1. 为什么是Elfouhaily谱从涌浪到毛细波的全尺度建模做海面电磁散射仿真、SAR海洋遥感或者海杂波建模的时候手里没有一套靠谱的海浪谱模型基本等于巧妇难为无米之炊。Elfouhaily海浪谱是我这几年用得最多的方向谱模型之一也是目前海洋遥感领域事实上的标准工具。它最打动我的一点是能把几毫米的毛细波和几百米的涌浪放进同一个公式框架里从大尺度重力波到小尺度毛细波平滑过渡中间不存在模型拼接的断点。很多刚接触海浪谱的同学会问PM谱和JONSWAP谱我还会用为什么要额外学一个Elfouhaily谱答案很直接PM谱描述充分发展的风浪JONSWAP谱描述有限风区的成长过程但这两个经典模型都只覆盖重力波范围对波数大于几十弧度每米的毛细波区间基本无能为力。而海面电磁散射的敏感性恰恰集中在小尺度波雷达后向散射的主要贡献来自布拉格共振波波长通常在厘米量级。用PM谱或者JONSWAP谱去算后向散射系数高频段完全就是一片空白。Elfouhaily谱是1997年由法国海洋开发研究院的几位学者在论文A unified directional spectrum for short and long wind-driven waves中提出的核心思想是把长波段的曲率谱和短波段的曲率谱统一到一个表达式中。这个谱以波数k为自变量同时包含波向分布天然适合电磁散射仿真和SAR图像模拟。如果你是做海杂波统计特性分析、海面双站散射计算或者多普勒谱仿真Elfouhaily谱几乎是绕不开的起点。这篇博文我会从模型原理讲到MATLAB代码实现再讲到绘图和工程应用中的一些坑。看完你不仅能画出好看的谱图还能理解每个参数背后的物理含义包括逆波龄怎么取值、截断波数怎么选、算出来怎么验证。2. 公式拆解把Elfouhaily谱的每个参数都弄明白2.1 从曲率谱出发长波与短波的统一Elfouhaily谱最核心的变量是波数k单位是弧度每米。频谱密度S(k)与曲率谱B(k)的关系很简单S(k) B(k) / k³为什么要用曲率谱因为海浪谱的动态范围实在太大如果直接画S(k)从波数0.001到10000跨越七个量级谱密度值变化几十个数量级很难看出结构。而曲率谱B(k) k³S(k)把高频段的衰减趋势抵消掉一部分在loglog坐标下能清楚地看到两个特征区域低频段的重力波峰和高频段的毛细波平台。曲率谱由长波分量和短波分量相加得到B(k) B_l(k) B_h(k)其中长波分量描写重力波段的能量分布短波分量描写毛细波段。两个分量的表达式分别是B_l(k) 0.5 * α_p * (c_p / c) * F_p B_h(k) 0.5 * α_m * (c_m / c) * F_mα_p和α_m是平衡水平常数c_p是谱峰处的相速度c_m是最小相速度c c(k)是包含表面张力修正的相速度。这两个式子看起来简单但每个量都牵涉到一系列物理参数下面逐一拆开讲。2.2 风速、逆波龄、摩擦速度三个核心输入参数Elfouhaily谱的主要输入有三个10米高度风速U10、逆波龄Ω、风向。其中逆波龄Ω是最容易被忽视但在物理上极其重要的参数。逆波龄的定义是Ω U10 / c_p其中c_p是谱峰相速度。Ω越大说明风浪越年轻谱峰向高波数移动。完全发展的风浪Ω ≈ 0.84这时候谱峰相速度c_p ≈ U10/0.84 ≈ 1.19U10谱峰对应的波长通常在实际海况中能看到。年轻风浪的Ω可以从1.5一直到5谱峰波长大幅缩短。由Ω可以推出谱峰波数k_p g * Ω² / U10²这个公式非常常用我一般会先算k_p来判断谱峰应该出现在哪个波数位置。举个例子U1010 m/s、Ω0.84时k_p 9.81*0.84²/100 ≈ 0.069 rad/m对应波长2π/k_p ≈ 91米是典型的中等尺度风浪。如果Ω变成3k_p ≈ 0.88 rad/m对应波长只有7米左右这就是一个典型的年轻海况。摩擦速度u*需要从风速和拖曳系数算C_d (0.8 0.065 * U10) * 1e-3 u* U10 * √C_d这里的拖曳系数公式是Elfouhaily论文里采用的中等风速经验关系。U1010 m/s时C_d 1.45e-3u* ≈ 0.38 m/s。2.3 方向分布波向怎么展开方向谱的构造方式是S(k, φ) S(k) * Δ(k, φ)其中Δ(k, φ)是方向分布函数。完整的Elfouhaily方向谱包含多个傅里叶谐波项来描写上风向/下风向的不对称性但工程实现中绝大多数场景用二阶谐波就足够了也就是Δ(k, φ) (1/2π) * [1 Δ_2(k) * cos(2(φ - φ_wind))]这里的Δ_2(k)是方向集中度参数它同样是波数的函数。物理图像上长波方向性通常比较集中短波因为受局地风应力驱动和破碎过程影响方向分布相对更弥散。实际实现时Δ_2(k)由长波分量和短波分量加权合成Δ_2(k) (B_l * Δ_l B_h * Δ_h) / (B_l B_h)权重就是两个曲率谱分量的相对大小这保证在重力波主导的波段和毛细波主导的波段方向谱能平滑过渡。2.4 单位约定与谱积分这节内容虽然枯燥但必须搞清楚否则后续计算有义波高、后向散射系数的时候会出乱子。波数谱S(k)的单位是m³表示单位波数间隔内的海浪方差贡献。方向谱S(k,φ)的单位是m³/rad方向分布函数Δ量纲为一。总方差和方向谱之间的关系是m0 ∫∫ S(k, φ) * k * dφ * dk注意这里有个k的雅可比因子因为极坐标下的面积元是kdkdφ。很多人第一次算有义波高就对不上原因就是漏了这个k。如果只用一维谱算那就是m0 ∫ S(k) dk两种方式结果应该一致前提是方向分布函数满足归一化条件∫Δ(k,φ)dφ 1。有义波高H_s 4√m0。3. MATLAB实现从零写一个干净可复用的函数3.1 代码结构设计我习惯把代码拆成两个函数加一个演示脚本。第一个函数计算一维频率谱第二个函数在它的基础上叠加方向分布生成二维方向谱。这样划分的好处是单向用时只调一维函数逻辑清晰需要方向信息时再调二维函数不会被迫生成用不到的大矩阵。代码里有几个关键设计点需要提前说清楚。第一波数向量k必须是列向量角度向量phi必须是行向量这样才能利用MATLAB R2016b引入的隐式扩展机制直接生成二维矩阵不用写meshgrid循环。第二所有物理常数重力加速度g、最小相速度波数k_m提前定义好避免在循环里反复赋值。第三输入输出加注释明确每个变量的物理单位和含义这个习惯在工程协作时能省大量沟通成本。function [S, B, kp, ustar] elfouhaily_omni(U10, Omega, k) % elfouhaily_omni Elfouhaily一维全向波数谱 % 输入: % U10 - 10米高度风速 [m/s] % Omega - 逆波龄, 完全发展海况约0.84 % k - 波数列向量 [rad/m], 必须为正且不含0 % 输出: % S - 全向频谱密度 S(k) [m^3] % B - 曲率谱 B(k) k^3 * S(k) % kp - 谱峰波数 [rad/m] % ustar - 摩擦速度 [m/s] g 9.81; km 370; % 最小相速度对应的波数 [rad/m] % 拖曳系数与摩擦速度 Cd (0.8 0.065 * U10) * 1e-3; ustar U10 * sqrt(Cd); % 谱峰参数 kp g * Omega^2 / U10^2; cp U10 / Omega; % 考虑表面张力的相速度 c sqrt(g ./ k .* (1 (k / km).^2)); cm sqrt(2 * g / km); % 平衡水平常数 if Omega 0.35 alpha_p 0.006 * Omega^0.55; else alpha_p 0.006; end if ustar / cm 1 alpha_m 0.008; else alpha_m 0.008 * (1 log(ustar / cm)); end % 谱形函数 Fp exp(-1.25 * (kp ./ k).^2) .* ... exp(-Omega / sqrt(10) * (sqrt(k ./ kp) - 1)); Fm exp(-0.25 * (k / km - 1).^2); % 曲率谱分量 Bl 0.5 * alpha_p * (cp ./ c) .* Fp; Bh 0.5 * alpha_m * (cm ./ c) .* Fm; B Bl Bh; % 全向频谱 S B ./ k.^3; end这段代码有几点值得解释。相速度公式c √[g/k * (1(k/k_m)²)]里的表面张力项是关键k_m 370 rad/m对应约1.7厘米的波长正是毛细波和重力波的过渡区。如果去掉这一项高频段谱形会完全错误。α_p的表达式根据Ω是否大于0.35分为两段这是论文里的原始设定目的是防止逆波龄过小时出现不合理的谱峰隆起。α_m的对数修正项来自实测数据拟合物理背景是强风条件下短波平衡水平会随摩擦速度增大。3.2 方向谱函数方向谱函数在一维谱的基础上叠加方向分布。这里的核心是方向集中度参数Δ_2(k)的计算以及矩阵广播的维度匹配。function [S_dir, D] elfouhaily_dir(U10, Omega, k, phi, wind_dir) % elfouhaily_dir Elfouhaily二维方向谱 % 输入: % U10 - 10米高度风速 [m/s] % Omega - 逆波龄 % k - 波数列向量 [rad/m], Nx1 % phi - 方位角行向量 [rad], 1xM % wind_dir - 风向方位角 [rad] % 输出: % S_dir - 方向谱 S(k,phi) [m^3/rad], NxM % D - 方向分布函数, NxM % 先算一维分量 [S_omni, B, kp, ustar] elfouhaily_omni(U10, Omega, k); g 9.81; km 370; c sqrt(g ./ k .* (1 (k / km).^2)); cm sqrt(2 * g / km); cp U10 / Omega; if Omega 0.35 alpha_p 0.006 * Omega^0.55; else alpha_p 0.006; end if ustar / cm 1 alpha_m 0.008; else alpha_m 0.008 * (1 log(ustar / cm)); end Fp exp(-1.25 * (kp ./ k).^2) .* ... exp(-Omega / sqrt(10) * (sqrt(k ./ kp) - 1)); Fm exp(-0.25 * (k / km - 1).^2); Bl 0.5 * alpha_p * (cp ./ c) .* Fp; Bh 0.5 * alpha_m * (cm ./ c) .* Fm; % 方向集中度参数 Delta_l 0.35 * exp(-Omega / sqrt(10) * (sqrt(k ./ kp) - 1)); Delta_h 0.35 * exp(-0.25 * (k / km - 1).^2); Delta2 (Bl .* Delta_l Bh .* Delta_h) ./ (Bl Bh); % 方向分布函数: 利用隐式扩展生成 N x M 矩阵 dphi phi - wind_dir; % 1 x M D (1 / (2 * pi)) * (1 Delta2 .* cos(2 * dphi)); % N x M % 方向谱 S_dir S_omni .* D; % N x M end关于方向集中度参数我在这里用了工程上最常见的简化形式长波部分的集中度随波数衰减短波部分在k_m附近有一个高斯形状的分布。更严格的Elfouhaily论文版本还包含上风向/下风向不对称的奇次谐波项对多普勒谱和极化散射计算有影响但对绝大多数谱图绘制和应用场景来说这个简化版本已经足够。后续需要精确复现论文结果时再把高阶项补上即可。3.3 验证代码正确性写完代码第一件事不是画图而是验证。我通常做三个检查第一个检查是看谱峰位置。把k_p g*Ω²/U10²算出来然后画出谱找峰值对应的波数两者应该基本一致。第二个检查是算有义波高。完整发展海况下U1010 m/s时有义波高大约2米左右。如果算出来是0.2米或者20米那一定是什么地方出了问题。% 有义波高检查 m0 trapz(k, S_omni); Hs 4 * sqrt(m0); fprintf(U10 %.1f m/s, Omega %.2f, Hs %.2f m\n, U10, Omega, Hs);第三个检查是方向谱的归一化。对任意固定波数k_0把方向谱沿所有角度积分结果应该等于对应的一维谱值∫ S(k_0, φ) dφ S(k_0)这个检查能验证方向分布函数的归一化系数是否正确。4. 谱图绘制把物理特征画到读者眼前4.1 一维频谱的绘制一维谱用loglog坐标画是最合适的因为波数从0.001到10000跨越七个数量级线性坐标根本没法看。一个典型的绘制代码如下% 演示脚本: 绘制Elfouhaily一维谱 clear; clc; U10 10; Omega 0.84; k logspace(-3, 4, 800); [S, B, kp] elfouhaily_omni(U10, Omega, k); figure(Color, w, Position, [100 100 900 380]); subplot(1, 2, 1); loglog(k, S, b-, LineWidth, 2); hold on; xline(kp, r--, LineWidth, 1.5); grid on; xlabel(波数 k (rad/m)); ylabel(S(k) (m^3)); title([U10 num2str(U10) m/s, \Omega num2str(Omega)]); xlim([1e-3 1e4]); legend(Elfouhaily, 谱峰位置, Location, southwest); subplot(1, 2, 2); loglog(k, B, r-, LineWidth, 2); grid on; xlabel(波数 k (rad/m)); ylabel(曲率谱 B(k) k^3 S(k)); title(曲率谱); xlim([1e-3 1e4]);画出来的图你应该能看到三个特征。第一低频段谱线快速上升到达谱峰谱峰位置与k_p一致。第二中频段谱线以k⁻³的斜率缓慢下降对应的曲率谱近似平台。第三在k 370 rad/m附近出现第二个凸起这是毛细波的能量聚集区也是曲率谱上最明显的短波特征。4.2 二维方向谱的绘制二维方向谱的绘制比一维复杂要把极坐标下的(k, φ)网格转换到笛卡尔坐标系。转换关系是k_x kcos(φ)k_y ksin(φ)。由于谱值的动态范围很大直接画线性值会丢失细节我习惯用分贝值来显示。% 二维方向谱绘制 phi linspace(-pi, pi, 360); [S_dir] elfouhaily_dir(U10, Omega, k, phi, wind_dir); % 转换到笛卡尔波数空间 [K, PHI] meshgrid(k, phi); X K .* cos(PHI); Y K .* sin(PHI); % 分贝显示, 以最大值为参考 SdB 10 * log10(S_dir / max(S_dir(:))); figure(Color, w, Position, [100 100 800 650]); surf(X, Y, SdB, EdgeColor, none); view(2); axis equal tight; colormap(jet); caxis([-40 0]); colorbar; xlabel(k_x (rad/m)); ylabel(k_y (rad/m)); title(二维谱 S(k,\phi) 归一化 dB 显示);这图画出来应该是一个沿风向拉长的椭圆形能量分布中心是重力波主峰外围扩散出去方向性比较集中。同时你会在k≈370附近看到第二圈能量对应毛细波这圈能量的方向分布比较宽整个图形看起来有点像嵌套的椭圆环。4.3 多工况对比单看一条谱线很难感受到参数的影响我建议做一个多工况叠加图% 不同风速对比 U10_list [5, 10, 15]; Omega 0.84; figure(Color, w, Position, [100 100 900 400]); hold on; colors lines(3); for i 1:3 [S_i, ~, kp_i] elfouhaily_omni(U10_list(i), Omega, k); loglog(k, S_i, Color, colors(i, :), LineWidth, 2); end grid on; xlabel(波数 k (rad/m)); ylabel(S(k) (m^3)); xlim([1e-3 1e4]); legend(U105 m/s, U1010 m/s, U1015 m/s, Location, southwest);风速的影响非常直观风速增大谱峰向低波数移动波长变长谱密度整体抬升有义波高显著增大。逆波龄的影响则不同% 不同逆波龄对比 Omega_list [0.84, 2, 3.5]; U10 10; figure(Color, w, Position, [100 100 900 400]); hold on; for i 1:3 [S_i, ~, kp_i] elfouhaily_omni(U10, Omega_list(i), k); loglog(k, S_i, Color, colors(i, :), LineWidth, 2); end grid on; xlabel(波数 k (rad/m)); ylabel(S(k) (m^3)); xlim([1e-3 1e4]); legend(\Omega0.84 成熟海, \Omega2.0, \Omega3.5 年轻海, ... Location, southwest);逆波龄增大时谱峰向高波数明显移动谱峰变尖锐总能量下降。这对应一个物理事实同样的风速下年轻海浪没有充分发展有效波高比成熟海况低。5. 参数分析与数值计算要点5.1 逆波龄Ω怎么取很多初学者最纠结的就是Ω取值。这里有三个典型场景完全发展海况取0.84一般工程默认值有限风区或存在强风急流的海况取1到2台风或飓风年轻浪况取3到5。如果你在做电磁散射仿真需要覆盖从成熟到年轻海面的不同风况我建议把Ω作为一个扫描参数而不是固定值这样能同时看到谱形状变化对后向散射系数的影响。需要注意的是Ω和U10不是独立的。给定风速下Ω由风区长度和风时共同决定。严格的海洋学计算要用风区模型反推但工程仿真中直接指定Ω是常见做法。有一点必须记住Ω过小时k_p会非常小可能低于你设定的波数下限。比如U105、Ω0.5时k_p 9.81*0.25/25 ≈ 0.098 rad/m如果波数下限设成0.1谱峰就落在计算范围外了。5.2 波数截断怎么选波数范围的选取取决于应用目的。电磁散射仿真需要覆盖雷达 Bragg 波数比如C波段5.3 GHz的Bragg波数大约在110 rad/m左右X波段9.6 GHz大约在200 rad/m左右Ku波段更高。这意味着你的波数上限至少要达到几百甚至上千才能正确计算高频段的后向散射。下限的选取取决于你是否关心涌浪。如果只关心风浪本身从0.01 rad/m开始就够了如果要模拟海面大尺度起伏对导航雷达的干扰下限要延伸到0.001甚至更低。波数向量用logspace生成是最合理的我之前用线性间隔试过低波数段完全没分辨率高波数段又浪费大量内存。5.3 从谱到有义波高的数值积分算有义波高的关键是数值积分精度。谱在谱峰附近变化很快如果波数采样太粗trapz积分会显著低估总能量。我的经验是每十倍频程至少100个采样点800到1000个点覆盖0.001到10000这个范围比较稳妥。积分代码如下m0 trapz(k, S); Hs 4 * sqrt(m0);验证一下量级U1010、Ω0.84时算出来的Hs应该在2米左右。如果偏差超过20%检查两个地方一是波数下限是否覆盖了谱峰左侧的上升段二是波数上限是否够高。低频段的能量对m0贡献很大因为谱密度在低频段虽然数值小但波数分辨率跨度大。5.4 方向集中度的物理验证画完二维谱后建议做一个快速物理检查在谱峰波数处方向分布函数D(k_p, φ)的半功率宽度应该在30到40度左右这是成熟海况的典型值。如果宽度只有5度说明Δ_2参数设置过大如果超过60度说明集中度不足。短波端的半功率宽度应该明显更宽通常在50到70度。这种直观检查比我之前见过的任何单元测试都管用。6. 常见问题与排查实录6.1 谱高频段出现不正常的振荡这是我遇到过最多的问题症状是谱线在波数超过1000 rad/m之后出现波浪形起伏甚至变成负值。原因通常是两种一是波数向量中包含0或者负值导致相速度公式出现除零和虚数二是矩阵计算中隐式扩展的维度搞错比如k不小心变成了行向量导致运算结果形状不对。检查方法很简单在函数入口加一句assert确认k的所有元素大于零且单调递增。再disp一下size(k)和size(S)看维度是否一致。我见过有人把k写成1x800的行向量谱线照样画出来了但方向谱全错。6.2 方向谱中心出现凹陷或者十字形伪影方向谱画出来中心凹陷通常是因为绘图时surf网格太粗糙或者caxis动态范围设得太大导致峰值区域饱和。另一个常见问题是meshgrid的维度顺序搞反了。在我上面的代码里meshgrid(k, phi)生成的是MxN矩阵所以S_dir需要转置才能匹配。如果忘了转置画出来的方向谱会相对于风向旋转90度。十字形伪影多半是数值问题。当波数k非常小或非常大时谱密度值可能溢出在loglog坐标下产生Inf或NaN。解决办法是把波数上下限设得保守一点并在谱密度计算后加一个有限值检查S(~isfinite(S)) 0;6.3 算出的有义波高和实测差很远有义波高偏大通常是因为波数下限太低把不存在的涌浪能量也算进去了。偏小则是因为波数上限不够或者采样太粗。另外要注意Elfouhaily谱描述的是风浪谱不包含涌浪成分。如果你的仿真场景是风浪和涌浪混合的海况直接在Elfouhaily谱上叠加一个窄带涌浪谱是常见的做法。6.4 老版本MATLAB的兼容问题隐式扩展是R2016b才有的特性如果实验室还有人在用老版本代码会报矩阵维度不匹配。解决办法是手动用bsxfun或者repmat。举一个例子% R2016a及更早版本 D (1 / (2 * pi)) * (1 bsxfun(times, Delta2, cos(2 * dphi)));这个改动只影响方向谱函数一维函数完全不受影响。6.5 性能优化二维谱在精细网格下会占用大量内存。比如k取1500个点、phi取720个点S_dir就是1500×720的双精度矩阵大约8.6 MB看起来不大。但如果做蒙特卡洛仿真要生成几千个海面实现这个矩阵会被反复计算时间成本就上来了。我常用的优化手段有三个把k和phi的采样点数降到合理的最低值把双精度改成单精度存储方向谱矩阵以及在循环外部预计算所有不依赖随机变量的谱分量。7. 应用拓展Elfouhaily谱能带你去哪7.1 海杂波仿真拿到方向谱之后最常见的下一步就是生成随机海面。方法是对谱做傅里叶逆变换在频域构造一个符合谱密度分布的随机复幅值然后逆变换到空间域。具体操作是取sqrt(S(k,φ)kdk*dφ)作为幅度相位用均匀随机数然后做二维IFFT。这一步做出来就是带方向性的随机海面高度场可以直接叠加到电磁散射计算当中。用Elfouhaily谱生成的海面有一个优点就是小尺度波纹的统计特性是自洽的。我曾经用不同截断波数生成海面并计算后向散射系数对比发现只要截断波数覆盖了布拉格共振区结果就非常稳定。这是用PM谱或者高斯谱做不到的。7.2 海面电磁散射计算的衔接在双站散射和SAR成像仿真中Elfouhaily谱提供了每个波数-方向分量上的能量密度结合微扰法或者小斜率近似就能算出归一化雷达散射截面。具体而言Bragg散射的核心是把共振波数处的谱密度提取出来代入到散射系数公式里。方向谱的优势在这里体现得特别明显它可以直接给出任意雷达入射方位角下的散射系数变化也就是方位向的调制效应。7.3 进一步扩展的方向如果后续研究需要更精细的谱模型可以在Elfouhaily基础上做三件事叠加涌浪谱和风浪谱形成混合谱引入时变风场让谱随空间位置渐变将谱离散化成海面网格时加入各向异性的调制。这些扩展在工程中都有现成需求理解Elfouhaily谱的实现细节是所有这些工作的共同基础。回到这篇博文的起点我当时自己动手写Elfouhaily谱的MATLAB代码本意只是做一个海面生成的前置步骤没想到在核对公式和排查问题的过程中把海浪谱模型的前因后果摸了个遍。前前后后踩过的坑里最想提醒后来人的一点是拿到公式不要急着抄代码先把每个参数的物理含义和量纲搞清楚。我见过太多人把α_m的临界判断写成大于号或者忘记相速度里的表面张力修正最后画出来的谱形不对还反过来怀疑模型本身。代码本身不难难的是理解每一个数值背后的海洋学意义。