Elfouhaily海浪谱详解：MATLAB实现海面电磁散射仿真统一谱模型

做海面电磁散射仿真那阵子我一直在找一个能同时覆盖重力波和毛细波的海浪谱模型。PM谱、JONSWAP谱虽然很成熟但都主要描述重力波段雷达后向散射关心的布拉格波数通常在几十到几百 rad/m恰好落在毛细波到短重力波的过渡段。后来换成 Elfouhaily 谱问题一下子理顺了它把长波和短波统一到同一个曲率谱框架里还自带方向分布非常适合海面建模和遥感仿真。这篇文章把我用 MATLAB 绘制 Elfouhaily 海浪谱的完整过程写出来包括公式拆解、可直接运行的代码、画图时容易踩的坑以及参数调试的体会。1. 为什么选 Elfouhaily 谱长波短波一把抓的“统一谱”1.1 PM、JONSWAP 这些经典谱卡在哪先说结论不是经典谱不好而是它们的设计目标决定了适用范围。Pierson-Moskowitz 谱描述的是充分成长状态下的重力波形式简单只有一个风速参数JONSWAP 谱在 PM 基础上加了峰增强因子更适合有限风区成长中的海浪。但这两个谱都有一个共同边界它们默认只考虑重力波表面张力那部分完全没有。这个缺陷在波浪统计和港口工程里影响不大但放到电磁散射、海面遥感、雷达杂波仿真里就麻烦了。微波散射对海面的敏感尺度是厘米到分米级对应波数大致在 10~1000 rad/m。这个区间正好是短重力波向毛细波过渡的地带表面张力开始主导色散关系仅仅靠重力波谱外推是不行的。我最早用过普通的截断处理把 PM 谱向高波数延伸再加个指数衰减。结果谱形虽然“像那么回事”但物理上没有依据换成不同风速时高频段的能量变化不可控仿真出来的海面偏“软”散射系数怎么调都不对。后来换到 Elfouhaily 谱才算是把长波和短波从同一个公式里自然长出来。1.2 Elfouhaily 谱的工程定位Elfouhaily 等人在 1997 年提出的这个谱目标很明确给短波和长波一个统一描述并且能直接用于电磁散射计算。它有几个特点是我在实际项目里最看重的波数覆盖范围宽从几十米长的涌浪到厘米级毛细波都在一个公式里长波部分和短波部分各有独立的饱和系数物理上对应不同的平衡机制自带方向扩展函数能描述风浪沿风向的分布不是简单的无方向谱公式里的参数都可以由 10m 风速和逆波龄反演输入参数少工程上方便。当然它也是半经验模型不是万能的。它对“单一风场、充分发育海面”这类场景拟合得比较好如果是复杂涌浪叠加或者风场剧烈变化的海况还是要配合其他谱或者数值模型来做。但作为默认的海面谱它在仿真和遥感反演里的使用频率非常高值得把实现细节摸透。2. 画图前先把公式和量纲理顺2.1 波数谱 S(k) 和曲率谱 B(k) 别混新手最容易搞混的就是这两个量。Elfouhaily 谱的标准形式是波数谱 (S(k))它满足海面高度方差积分[ \int_0^\infty S(k),dk \langle \eta^2 \rangle ]这里 (k) 是角波数单位是 rad/m所以 (S(k)) 的单位是 m³/rad。经常在文献里看到的曲率谱 (B(k)) 或者饱和谱定义是[ B(k) k^3 S(k) ]这个量是无量纲的它描述的是海面斜率谱的饱和程度。画图时如果直接画 (S(k))动态范围太大峰值和高频段差了七八个数量级画 (B(k)) 就能更清楚地看到重力波平衡段和毛细波平衡段。所以我一般在一个图里同时画这两个量左边看绝对能量右边看谱形结构。2.2 输入参数怎么取这个谱的核心输入不多我整理成一张表参数符号含义典型取值10m 风速(U_{10})海面上方 10m 处的风速5~20 m/s逆波龄(\Omega)表征风浪成长状态0.84 为充分成长1~3 为年轻浪风向(\phi_w)风浪传播主方向0 rad重力加速度(g)物理常数9.81 m/s²毛细波过渡波数(k_m)表面张力开始显著作用的波数370 rad/m最小相速度(c_m)毛细波过渡区的最小相速度0.23 m/s逆波龄 (\Omega) 是个容易被忽略的参数。它的定义是[ \Omega \frac{U_{10}}{c_p} ]其中 (c_p) 是谱峰值处的相速度。(\Omega0.84) 对应充分成长的风浪(\Omega) 越大说明波浪越“年轻”峰值波数越高。同一个风速下(\Omega) 不同谱的形状差别会很直观。峰值波数 (k_p) 的计算公式很简洁[ k_p \frac{g \Omega^2}{U_{10}^2} ]这个公式建议直接记住后面做参数敏感性分析时很有用。举个例子充分成长状态下 (\Omega0.84)不同风速对应的峰值波数和波长如下(U_{10}) (m/s)(k_p) (rad/m)谱峰波长 (m)50.27722.7100.069290.8150.0308204.3200.0173363.0风向 (\phi_w) 默认取 0 就行但方向谱画图时一定要把它作为参数传进去否则后面调风向时容易晕。2.3 核心公式里的关键物理量在写代码之前把公式完整过一遍。色散关系用的是包含表面张力的形式[ c(k) \sqrt{\frac{g}{k}\left(1\frac{k^2}{k_m^2}\right)} ]峰值相速度直接用[ c_p \frac{U_{10}}{\Omega} ]摩擦风速用一个简化拖曳系数计算[ u_* U_{10} \sqrt{0.00144} ]一维谱分成两部分长波项 (B_L(k)) 和短波项 (B_H(k))[ S(k) \frac{B_L(k) B_H(k)}{k^3} ]长波项[ B_L(k) 0.5 \alpha_p \frac{c_p}{c} F_p(k) ][ F_p(k) \exp\left[-\frac{5}{4}\left(\frac{k_p}{k}\right)^2\right] \exp\left[-\frac{\Omega}{\sqrt{10}}\left(\sqrt{\frac{k}{k_p}}-1\right)\right] ]短波项[ B_H(k) 0.5 \alpha_m \frac{c_m}{c} F_m(k) ][ F_m(k) \exp\left[-\frac{1}{4}\left(\frac{k}{k_m}-1\right)^2\right] ]系数 (\alpha_p) 和 (\alpha_m) 分别控制长波、短波的饱和水平[ \alpha_p 0.006\sqrt{\Omega} ][ \alpha_m \begin{cases} 0.01\left(1\ln\frac{u_}{c_m}\right), u_ c_m \ 0.01\left(13\ln\frac{u_}{c_m}\right), u_\ge c_m \end{cases} ]方向分布函数[ D(k,\phi) \frac{1}{2\pi}\left[1\Delta(k)\cos\left(2(\phi-\phi_w)\right)\right] ][ \Delta(k) \tanh\left[\frac{\ln 2}{4}\frac{u_*}{c}\frac{c_p}{c}\frac{1}{\alpha_p}\right] ]最终二维方向谱就是[ \Psi(k,\phi) S(k) D(k,\phi) ]我第一次看到这一串公式时觉得头大但实际上写成代码很机械关键是别把单位搞错尤其是ln和log10。3. 核心代码一个函数搞定一维谱和方向分布3.1 可直接运行的 MATLAB 函数我按上面的公式写了一个函数既返回一维谱 (S(k))也返回二维方向谱 (\Psi(k,\phi))。把它存成elfouhaily_spectrum.m就能直接用。function [Sk, Psi] elfouhaily_spectrum(U10, k, phi, Omega, phiw) % elfouhaily_spectrum 计算Elfouhaily海浪谱 % 输入: % U10 : 10m风速, m/s % k : 角波数, rad/m, 可为向量 % phi : 方位角, rad, 可为向量 % Omega : 逆波龄, 默认0.84 % phiw : 风向, rad, 默认0 % 输出: % Sk : 一维波数谱 S(k), 1 x length(k) % Psi : 二维方向谱 Psi(k, phi), length(phi) x length(k) if nargin 4 || isempty(Omega) Omega 0.84; end if nargin 5 || isempty(phiw) phiw 0; end g 9.81; km 370; % 重力-毛细波过渡波数, rad/m cm 0.23; % 最小相速度, m/s Cd 0.00144; % 简化拖曳系数 u_star U10 * sqrt(Cd); kp g * Omega^2 / U10^2; cp U10 / Omega; % 包含表面张力的色散关系 c sqrt( g ./ k .* (1 (k / km).^2) ); % 长波项 alpha_p 0.006 * sqrt(Omega); Lpm exp(-1.25 * (kp ./ k).^2); Gamma exp(-(Omega / sqrt(10)) .* (sqrt(k ./ kp) - 1)); Fp Lpm .* Gamma; BL 0.5 * alpha_p .* (cp ./ c) .* Fp; % 短波项 if u_star cm alpha_m 0.01 * (1 log(u_star / cm)); else alpha_m 0.01 * (1 3 * log(u_star / cm)); end Fm exp(-0.25 * (k / km - 1).^2); BH 0.5 * alpha_m .* (cm ./ c) .* Fm; % 一维谱 Sk reshape((BL BH) ./ k.^3, 1, []); % 方向扩展项 Delta tanh( log(2) / 4 .* (u_star ./ c) .* (cp ./ c) / alpha_p ); Delta reshape(Delta, 1, []); % 生成二维网格 [~, Phigrid] meshgrid(k, phi); Dphi (1 / (2 * pi)) * (1 Delta .* cos(2 * (Phigrid - phiw))); % 二维方向谱 Psi Dphi .* Sk; end这个函数符合 MATLAB 的向量化习惯k和phi传向量就行。有一点要提醒log在 MATLAB 里默认是自然对数千万别为了“保险”改成log10算出来的谱会完全不对。3.2 方向扩展项为什么这样写方向扩展项是我最初实现时最没把握的地方。如果只看一些简化的博客代码很多人会把方向分布写成固定的 (\cos^2) 分布这样虽然能出图但丢失了 Elfouhaily 谱的一个重要特征方向扩展随波数变化。长波段的峰值波数附近风浪的方向性比较强谱能量集中在风向上到了毛细波段方向性逐渐变弱能量分布向各向同性过渡。 (\Delta(k)) 这个参数就负责控制这个变化过程。公式里cp ./ c和1 / alpha_p这两个因子不能省。cp/c描述相速度比相当于把长波峰的“记忆”传递到当前波数1/alpha_p引入逆波龄的影响让方向扩展程度和波浪成长状态挂钩。我第一次偷懒去掉之后画出来的方向图几乎是个圆跟实际风浪差太远了所以这个地方一定要按公式来。3.3 波数网格怎么布置画谱之前要先生成波数向量。我的经验是k logspace(-3, 3, 512);也就是从 (10^{-3}) 到 (10^3) rad/m对数均匀取 512 个点。这个范围基本覆盖了涌浪到毛细波的主要能量区间点数也足够让峰值的形状光滑。如果你对尺度有特定要求比如只看重力波段可以把范围缩到logspace(-2, 1, 256)但如果要做雷达散射仿真高频段一定要保留到 1000 rad/m 以上否则布拉格波数附近没有谱值。方位角网格我一般用phi linspace(-pi, pi, 361);如果做极坐标方向图361 个点足够平滑。二维直角坐标画图时phi 的数量可以降到 181因为pcolor的网格显示天然有平滑效果点数太多反而增大内存。4. 从一维到二维三张图把谱画明白4.1 一维图S(k) 与 B(k) 双轴对比画一维谱是最直接的验证方式。我习惯把 (S(k)) 和 (B(k)) 放在同一张图用双纵轴方便比较。U10 10; Omega 0.84; phiw 0; k logspace(-3, 3, 512); [Sk, ~] elfouhaily_spectrum(U10, k, [], Omega, phiw); Bk k.^3 .* Sk; figure; yyaxis left; loglog(k, Sk, LineWidth, 1.5); ylabel(S(k) (m^3/rad)); yyaxis right; loglog(k, Bk, LineWidth, 1.5); ylabel(B(k)); xlabel(wavenumber k (rad/m)); legend({S(k), B(k)}, Location, southwest); set(gca, XScale, log, YScale, log);画出来能看到几个特征(S(k)) 在 (k_p) 附近有明显的峰然后向高频段衰减(B(k)) 在高频段会有一个隆起对应的就是毛细波平衡段谱峰的位置和风速、逆波龄有直接对应关系。这张图也是检查代码是否写对的第一步。如果 (B(k)) 在高频段没有隆起或者峰值位置明显不合理多半是系数或者单位出了问题。4.2 方向分布极坐标图方向分布最能直观看出“风浪沿风向展开”的效果。我用若干特征波数画极坐标图把不同尺度波浪的方向性并列展示。phi linspace(-pi, pi, 361); k_test [0.05, 0.1, 1, 10, 370]; figure; for i 1:numel(k_test) [~, Psi_i] elfouhaily_spectrum(U10, k_test(i), phi, Omega, phiw); polarplot(phi, Psi_i ./ max(Psi_i), LineWidth, 1.5); hold on; end legend(cellfun((x) sprintf(k%.1f rad/m, x), ... num2cell(k_test), UniformOutput, false), Location, northeast);用polarplot的好处是不用手动把角度转成直角坐标。归一化之后能清楚看到波数越大花瓣越宽也就是方向性越弱。这个图放进论文或者汇报材料里很直观。4.3 二维方向谱的直角坐标渲染二维谱的正确定义是在极坐标下的 ( \Psi(k,\phi) )但展示时通常要转到 (k_x,k_y) 直角坐标平面。这里有个容易忽视的细节在 (k_x,k_y) 平面上谱密度需要除以 (k)因为直角坐标面积元和极坐标面积元的关系是[ dk_x dk_y k,dk,d\phi ]所以如果直接画Psi2积分关系会不对。正确做法是先将Psi2除以K得到直角坐标谱密度kx linspace(-1.5, 1.5, 401); ky linspace(-1.5, 1.5, 401); [KX, KY] meshgrid(kx, ky); K sqrt(KX.^2 KY.^2); PHI atan2(KY, KX); % 避免零波数处除零 K(K 1e-3) 1e-3; [~, Psi2] elfouhaily_spectrum(U10, K(:), PHI(:), Omega, phiw); Psi2 reshape(Psi2, size(K)); % 转换到 kx-ky 平面的谱密度 Fcart Psi2 ./ K; figure; pcolor(KX, KY, log10(Fcart)); shading interp; axis equal; colorbar; clim([-4 2]); xlabel(k_x (rad/m)); ylabel(k_y (rad/m));画出来的二维图能很清晰地看到风浪能量在风向上的“长条形”分布。要注意clim的取值范围会因为风速和波数范围变化第一次画的时候可以先用colorbar看一下数据范围再手动调整。4.4 顺手做个海面高度场可选谱画完之后很多人的下一步是生成随机海面高度场。思路并不复杂对直角坐标谱密度开根号乘上随机相位再做逆傅里叶变换。核心代码大致是dk kx(2) - kx(1); A sqrt(2 * Fcart * dk^2) .* exp(1i * 2 * pi * rand(size(Fcart))); eta real(ifft2(ifftshift(A)) * numel(A));这里要注意随机相位的单位以及ifftshift的用法不同 MATLAB 版本对频率排列的要求略有差异。这个方向展开去写又是一篇文章本文先不展开但谱本身已经为它提供了完整输入。5. 参数调一调谱就变脸风速和逆波龄的影响5.1 风速 U10峰值波数整体平移风速对谱的影响最容易理解。从 (k_p g\Omega^2/U_{10}^2) 可以看出风速越大峰值波数越小也就是海浪的主波长越长。能量水平也跟着整体抬升。实际计算时你会发现(U_{10}5) m/s 和 (U_{10}15) m/s 的谱峰值位置差了将近一个数量级。如果海面仿真的风场是空间变化的意味着不同区域的谱峰波数差异很大网格分辨率要提前留足余量。还有一个不太直观的点风速增大后(u_*) 会超过 (c_m)此时 (\alpha_m) 的计算会从1 log切换到1 3*log短波段的饱和水平会明显上升。这对应着强风下毛细波能量增强雷达后向散射也会跟着增强。这个切换点大约在 (U_{10} \approx 6) m/s 附近低风速仿真时要留意。5.2 逆波龄 Ω从“成熟海”到“年轻海”逆波龄 (\Omega) 是我觉得这个谱里最有意思的参数。(\Omega0.84) 表示充分成长的风浪实际海面往往达不到这个状态尤其是有限风区波浪还没长起来(\Omega) 可能到 1.5 甚至 2 以上。当 (\Omega) 增大时峰值波数 (k_p) 变高主波长变短(\alpha_p) 变大长波部分饱和水平上升方向扩展项里的 (1/\alpha_p) 变小整体方向扩展变宽。翻译成物理图像就是年轻浪更“碎”、能量分布更宽、方向性更弱。这个参数对雷达海杂波仿真的影响很大因为不同波龄下的短波能量分布差别明显。5.3 自检积分还原和数值稳定写完函数之后我强烈建议做一次自检把二维方向谱对 (\phi) 积分应该能还原出一维谱 (S(k))。phi linspace(-pi, pi, 361); [Sk, Psi] elfouhaily_spectrum(U10, k, phi, Omega, phiw); % 积分还原一维谱 Sk_check trapz(phi, Psi, 1); figure; loglog(k, Sk, LineWidth, 1.5); hold on; loglog(k, Sk_check, --, LineWidth, 1.5); legend({S(k), integral of Psi}, Location, southwest);如果两条线重合说明方向分布函数的归一化正确。我一开始画极坐标图时发现方向图面积不对劲就是因为忘了方向分布前面除以 (2\pi)。这个自检能省掉很多排查时间。数值稳定性方面注意k不能取 0对数网格的下限本来就大于 0问题不大。但如果有人手动构造带 0 的波数向量函数里会出现除零建议在外部统一做k(k 1e-3) 1e-3。6. 我踩过的坑和调试建议6.1 单位换算和 log 函数这个坑我踩得最痛。Elfouhaily 谱的公式全部基于角波数 (k)单位是 rad/m。如果从频率谱或者波长谱转换过来要记得雅可比变换。例如波长谱 (S(\lambda)) 和波数谱 (S(k)) 的关系是[ S(\lambda) \frac{2\pi}{\lambda^2} S(k), \quad k \frac{2\pi}{\lambda} ]如果直接拿波长网格代入波数公式不换算谱峰位置会完全错位。别问我怎么知道的。另外再次强调MATLAB 的log是自然对数。公式里所有 ( \ln ) 和 ( \log ) 都用log实现即可不要改成log10。6.2 低频截止的陷阱低频段虽然没有峰值能量但涌浪部分对长波段仍有贡献。如果 (k) 的下限取得太大比如从 (10^{-1}) 开始谱峰左侧的形状会被截掉仿真海面会缺失长波起伏。我做海面高度场时kmin至少要比 (k_p) 小两个数量级。比如 (U_{10}10) m/s 时 (k_p\approx0.07)那kmin取 (10^{-3}) 比较稳妥。同样的道理kmax要覆盖到km的数倍以上否则毛细波段看不到饱和形态。6.3 绘图性能与样式细节如果只是画一维谱512 个点足够。二维直角坐标图要注意网格点数401 x 401已经能画出很光滑的图点数再翻一倍速度会明显变慢而且视觉上差别不大。绘图样式上我一般会把LineWidth设为 1.5坐标字号用Times New Roman图例放在不遮挡曲线的地方。科研绘图最大的原则是“看图的人第一眼就能抓到谱峰位置”所以对数坐标轴一定要标清楚物理量单位。6.4 最后的调试心得如果画出来的谱形不对先别急着调代码逻辑按下面顺序排查检查kp的数量级是否合理(U_{10}10) m/s 时 kp 应该在 0.07 附近检查alpha_m是否落在 0.005 到 0.05 之间如果出现负数多半是u_star/cm的对数取值不对检查一维谱积分还原二维谱是否通过这是最快定位方向分布问题的方法最后再看高频段 (B(k)) 是否在 (k370) 附近有突起这是毛细波段的标志没有突起说明短波项没生效。我个人在实际操作中的体会是Elfouhaily 谱的代码实现本身不复杂真正的难点在于理解每个参数背后的物理含义并且始终带着“积分、单位、谱形”三个校验去工作。只要你把这几点控制住这套代码不管是做海面仿真还是雷达散射系数计算都能当成一个稳定的底子继续往上层搭。