浅海传播损失图计算：简正波理论与Python实现

简介一份基于简正波理论的水声学实验程序包面向水声学、海洋声学方向的本科生或科研人员解决均匀水层与硬质海底波导模型中声传播损失的计算与分析问题。实验设定水层厚度100 m、声速1500 m/s海面海底平整点声源位于30 m深处发射25 Hz正弦信号程序分别计算10 m与90 m接收深度上的传播损失便于对比不同深度处的声场差异。压缩包共3个文件、约2 KB全部为MATLAB的.m脚本分别对应10 m、90 m以及二者联合计算的三份实现结构简洁适合直接运行、修改参数并观察简正波对传播损失的影响。目前已有509人学习下载适合用于课程实验、毕业设计或水声信道仿真入门。通过源码可了解简正波求和、深度采样及传播损失曲线绘制的基本流程并以此为基础扩展更多海底边界或声速剖面条件。1. 一百米水深、一百赫兹正弦信号声场里只装得下七个简正波传播损失图不是把声速剖面和海底参数丢进软件就能自动吐出来的东西它的质量取决于你到底在解哪个方程。对浅海几十米到两百米水深、频率几百赫兹以内的正弦信号射线声学要追几百条本征声线才能凑出像样的干涉结构而简正波通常只需要七八个模态就能把声场能量分布讲清楚。这个数量级上的差距让简正波成为浅海传播损失图计算绕不开的出发点。这篇文章要解决三个具体问题正弦信号作为单频声源时如何把频率换算成模态数量和截止条件砂质、泥质、岩质海底参数如何进入特征方程以及如何用几十行 Python 把传播损失图算出来并读懂干涉条纹、海底泄漏这些特征。适合做水声信道仿真、声呐性能预估和海底参数反演的一线工程师。2. 影响传播损失图的三个建模输入正弦信号、声速剖面与海底参数2.1 正弦信号频率先过两道门截止频率与模态数预判简正波理论里单频正弦激励等价于在某个固定角频率上求亥姆霍兹方程的特征展开。频率一旦定下来能传播的离散模态数就基本定死了。对最简单的 Pekeris 波导——等声速水层加均匀半空间海底第 n 阶模态的截止频率可以写成fc_n ≈ ((2n - 1) × c1) / (4H × sqrt(1 - c1² / c2²))其中 c1 是水层声速c2 是海底声速H 是水深。以 100 米水深、c11500 m/s、c21800 m/s 为例一阶模态截止频率约 6.8 Hz说明 10 Hz 以下信号在浅海几乎只剩下一个可传播模态干涉结构非常单调而 100 Hz 正弦信号对应的可传播模态数约为 Xmax/π 的整数部分Xmax H·sqrt(k1² - k2²)这里大约是 23.15除以 π 得 7.37也就是常说的“七个有效模态”。这个预判在做传播损失图之前就应该完成。它能帮你回答两个问题第一当前频率下用射线法还是简正波法更划算第二后续程序里最大模态数 N_max 该设成多少。实际工程中我一般会在算出的模态数基础上加 2 到 3 个用来覆盖海底折射带来的额外模态但不会盲目加到几十个——高阶模态接近截止角数值上很容易出现伪根。提示频率提高一倍模态数大致翻一倍但特征方程求根的数值难度也同步上升不要只看模态数来选方法。2.2 海底参数表声速、密度、衰减系数的工程取值与获取渠道海底在简正波方程里不是一条边界曲线而是三个物理量纵波声速 c2、密度 ρ2 和纵波衰减系数 α。这三个量决定特征方程的根落点进而决定传播损失图的全局形状。下面是一组工程常用的取值区间来源是 Hamilton 的地声模型和浅海实验反演结果的综合海底类型纵波声速 c2 (m/s)密度 ρ2 (kg/m³)衰减 α (dB/λ 1kHz)备注软泥1480–15501300–16000.1–0.2声速接近水层模态泄漏快粉砂质砂1600–17001700–19000.4–0.6常见大陆架浅海细砂1700–18501900–21000.5–0.8干涉条纹最清晰粗砂/砾石1850–20002000–22000.6–1.0高频衰减显著基岩2200–30002400–28000.1–0.3高声速强反射模态数偏多参数怎么拿分三种情况。有钻孔或地质采样时直接查柱状图有折射地震或走时层析结果时换算纵波声速什么都没有时用地声反演先估声速和密度再对衰减做敏感性扫描。我的建议是密度尽量别拍脑袋它直接影响特征方程的根位置反演时密度和声速存在强耦合拍错一个另一个必歪。衰减系数对近中场影响小对远场传播损失图的斜率影响大适合放在最后给误差带走。2.3 水层声速剖面Pekeris 均匀层假设什么时候成立真实海水声速随深度变化浅海常见的是春季等温层带来的负梯度、夏季温跃层带来的强负梯度以及混合层内的微弱正梯度。Pekeris 波导把水层简化为均匀声速这显然不严格但建模时得有取舍。判断标准很实用当声速剖面的相对变化率满足 (dc/dz)·H / c1 0.02 时用等效均匀声速代替传播损失图的误差通常能控制在 1 dB 以内如果温跃层强到百米水深内声速差超过 30 m/s就必须分层。分层的最简单做法是“阶梯近似”按 5 到 10 米一层把声速剖面板化成等声速层层间保留密度连续。这时特征方程不再是单个 tan 函数而是分层介质里逐层递推的阻抗关系。文章第 3 章先给均匀海底的完整 Python 实现第 5 章讲如何扩展到分层声速和分层海底——那是工程里最常见的升级路径。3. 用简正波计算正弦激励下的声场Python 最小实现3.1 特征方程的改写把复模态扫根问题变成有界分支求根Pekeris 波导的特征方程原始形式是 tan(X) -ρ1/ρ2 × k_z2 / k_z1其中 X k_z1 × H 是水层垂直波数的无量纲形式k_z2 是海底内的虚垂直波数。把 k_z2 用水平波数守恒关系消掉之后能得到一支只含 X 的实方程tan(X) -(ρ1/ρ2) × sqrt(Xmax² - X²) / XXmax 的含义前面已经算过它同时是模态总数的上限。方程左侧 tan(X) 在区间 ((n-1/2)π, nπ) 内从负无穷单调升到 0右侧在这个区间内是个负的有限值因此每一个 n 对应恰好一个根。这个“一区间一根”的结构很适合用 scipy.optimize.brentq 做有界求根不需要全局扫描也不怕漏根。import numpy as np from scipy.optimize import brentq f 100.0 # 正弦信号频率单位 Hz H 100.0 # 水深单位 m c1, c2 1500.0, 1800.0 # 水层/海底纵波声速 rho1, rho2 1000.0, 1800.0 # 密度单位 kg/m^3 alpha_b 0.5 # 海底衰减dB/lambda 1kHz w 2 * np.pi * f k1, k2 w / c1, w / c2 Xmax H * np.sqrt(k1**2 - k2**2) def find_root(n): L (n - 0.5) * np.pi R min(n * np.pi, Xmax - 1e-6) if L R: return None # 目标函数tan(X) (rho1/rho2) * sqrt(Xmax^2-X^2) / X 0 def F(X): return np.tan(X) (rho1 / rho2) * np.sqrt(Xmax**2 - X**2) / X try: return brentq(F, L 1e-6, R - 1e-6) except ValueError: return None roots [] n 1 while True: r find_root(n) if r is None: break roots.append(r) n 1这段代码的核心是把求根区间严格锁定在每个分支的内部避免 brentq 在 tan 的奇点附近抖动。注意 L 和 R 的取值左端点取 (n-1/2)π 偏移一个微小量右端点取 nπ 与 Xmax 的较小值两者都留了安全间隔。roots 列表里的每个值就是一个模态的 Xn它是接下来计算模态函数和水平波数的唯一输入。如果程序返回的模态数和你第 2.1 节预估的不一致先查 Xmax 计算有没有把频率和水深单位搞混。3.2 求解特征值并叠加声场归一化与源项处理的三个细节拿到 Xn 之后第 n 阶模态的水平波数是 k_rn sqrt(k1² - (Xn/H)²)。这里要留意k_rn 可能是个复数——当 Xn/H 超过 k1 时模态进入衰减区不过对 Pekeris 波导里正数分支的根这种情况很少发生真正的复数部分来自海底衰减 alpha_b 的等效吸收后面叠加声场时用指数衰减近似就够了不必把整个方程复数化那是全波数积分的事。模态函数的空间形态是 ψn(z) sin(Xn · z / H)归一化系数需要满足模态正交关系。严格做法是用包含密度权重的内积 ∫ψn²(z)/ρ(z) dz 1Pekeris 波导里密度在海底界面跳变这个积分要把水层和半空间分开算。下面的实现为保持可读性只在水平内积意义下做了归一化工程误差在可接受范围内如果要严格恢复声压绝对值需要把海底半空间项加回来krs np.sqrt(k1**2 - (np.array(roots) / H)**2).real # 模态函数在深度网格上的值z 取 0 到 H def mode_psi(z_grid, Xn): return np.sin(Xn * z_grid / H) # 归一化只考虑水柱部分教学用简化版本 def normalize(Xn): # 积分 sin^2(Xn*z/H) dz from 0 to H return np.sqrt(2.0 / (H - np.sin(2 * Xn) * H / (2 * Xn)))三个细节值得单独说。第一源深度 zs 不要取在水层边界上边界处模态函数可能恰好接近节点辐射效率算出来为零容易误判第二叠加声场时每个模态要乘以源深度的模态函数值 ψn(zs)这是点源激发系数的来源第三Hankel 函数的远场渐近式 1/sqrt(k_rn · r) 项决定了柱面扩散的几何衰减漏掉它传播损失图会整体偏移约 5 到 10 dB而且不同距离段偏移不一致。常见误用是直接拿模态数加权漏掉 k_rn 在分母上的距离项。3.3 生成传播损失图的完整脚本与输出判读这里给出可以完整运行的传播损失图生成脚本。声源是 100 Hz 正弦连续波源深度取 25 米接收距离从 500 米扫到 10 公里接收深度覆盖整个水柱。海底用细砂参数水层先假设等声速。水层体积衰减系数按 0.3 dB/km 1kHz 线性折算到当前频率。import matplotlib.pyplot as plt zs 25.0 rg np.linspace(500, 10000, 400) zg np.linspace(2.0, H, 200) alpha_w 0.3 * (f / 1000.0) # dB/km按频率线性折算 TL np.zeros((len(zg), len(rg))) for i, z in enumerate(zg): p np.zeros(len(rg), dtypecomplex) for Xn, kr in zip(roots, krs): ps np.sin(Xn * zs / H) pz np.sin(Xn * z / H) # 柱面波远场近似1/sqrt(kr*r) 为几何衰减指数项为相位与吸收 p ps * pz * np.exp(1j * kr * rg) p * np.exp(-alpha_w * rg / 1000.0) / np.sqrt(kr * rg 1e-6) TL[i, :] -20 * np.log10(np.abs(p) 1e-12) plt.figure(figsize(10, 6)) pc plt.pcolormesh(rg / 1000.0, zg, TL, shadingauto, cmapjet, vmin40, vmax90) plt.colorbar(pc, labelTransmission Loss (dB)) plt.xlabel(Range (km)); plt.ylabel(Depth (m)) plt.gca().invert_yaxis() plt.title(fTL map, {int(f)} Hz, source at {zs:g} m) plt.tight_layout(); plt.show()这段代码里p 的叠加顺序是“先相位后衰减”工程上有人先算衰减再乘相位结果画出来也没区别但前者更接近简正波理论表达式。vmin 和 vmax 的设置会让远场弱信号也能看出结构色不然靠近声源的位置会直接把整个色标拉爆。初看传播损失图时最优先确认三件事近场有没有以源深度为中心的对称扩张结构中场有没有明显的明暗相间条纹远场斜率是否稳定在每公里 3 到 5 dB 的量级。缺任何一条都意味着模态叠加或海底参数出了问题。4. 传播损失图读图与验证干涉条纹、海底泄漏和收敛性检查4.1 双模态干涉条纹在传播损失图上量出两种速度差传播损失图沿距离方向出现的明暗条纹本质是两阶主导模态的干涉拍。两阶模态水平波数差 Δk k_rm - k_rn对应的干涉周期是 Δr 2π / Δk。对 100 Hz、水深 100 米这个算例前两阶模态的 k_r 分别约为 0.4185 和 0.4175Δk 约 0.001干涉周期约 6 公里出头所以在 10 公里范围里只能看到一两个完整周期。要提高条纹可见度要么提高频率要么拉长接收距离段。实际操作里我一般会先在图上选一个接收深度把 TL(r) 沿距离提取出来做一次去趋势后找峰值间隔然后反推 Δk。这个反推值和特征方程求根的结果对不上时多半是某个主导模态被归一化过程吃掉了。另一种常见误读是把水平干涉条纹误当声影区边界。区分方法很简单——声影区边界在深度方向上有明显的突然跳变而干涉条纹从浅到深是连续倾斜的。4.2 海底泄漏与临界掠射角远离声源后为何声场贴底衰减更快海底声速比水层高时存在临界掠射角 θc arccos(c1/c2)本例是 33.6 度。实际海面声源激发的声波超过临界掠射角的能量每次碰海底都会向底下透射一部分这种损耗叫海底泄漏。在传播损失图上它表现为远离声源后声场能量逐渐向下沉近场的高阶模态携带大掠射角能量衰减特别快远场剩下的几乎是贴着临界角传播的低阶模态等效相位速度接近 c2。海底泄漏的量级可以直接从衰减 alpha_b 估算。远场情况下泄漏导致的额外衰减约等于 alpha_b × 反射次数 × 每反射路径长度反射次数粗算为 r / (2H·tanθ)其中θ取临界掠射角。对 10 公里传播距离、100 米水深这个值相当可观。这也是为什么传播损失图里最安全的接收器位置往往不是贴底而是放在声源深度或略低几米的地方——贴底位置收到的信号既不稳定又对海底衰减参数极其敏感。4.3 模态收敛性检查两种判据与 N_max 的“够用”标准计算完后不能直接拿图去写报告必须确认模态截断没有留下坑。两个判据缺一不可。第一个是能量判据最后一个参与计算的模态其等效声源强度 ψn(zs)² 相对所有模态总贡献要在 -40 dB 以下达不到就增加 N_max。第二个是残差判据把模态数从 N 加到 N5固定几个特征距离上的 TL 值变化应小于 0.1 dB。增加模态数对近场影响大还是远场影响大答案是近场。高阶模态在近场携带大量能量但衰减快远场只留低阶模态。如果收敛性检查只取远场 10 公里处的点近场可能早已失真。所以判据判定位置至少要取三个近场 500 米、中场 3 公里、远场 10 公里。有时扩大 N_max 后近场变好、远场不变那不是 bug是简正波方法的固有特性——高阶模态的泄漏衰减承担了近场到远场的平滑过渡。5. 从解析特征方程到分层海底相对层的扩展与最后调参技巧5.1 水层和海底做分层逐层递推求解与最小代码改动实际海域几乎没有严格均匀的海底。工程上最常见的扩展是把海底从半空间改成两层或三层每层有独立的声速、密度和衰减。特征方程的求解从 tan 函数变成逐层递推的阻抗关系但最小改动是保留原 Pekeris 求根框架把方程右侧替换成数值计算的有效阻抗 Z_b 除以水层特性阻抗def bottom_impedance(k_z1, layers): # layers: [(c, rho, alpha_att, h), ...]从最上层海底开始 Z None for c, rho, att, h in layers: k_z np.sqrt((w/c)**2 - (k1**2 - k_z1**2) 0j) rho_b rho * (1 1j * att / 8.686) if Z is None: Z rho_b * w / k_z else: # 分层介质递推式经多层反射校正 Z rho_b * w / k_z * (Z 1j * rho_b * w / k_z * np.tan(k_z * h)) / \ (rho_b * w / k_z 1j * Z * np.tan(k_z * h)) return Z这个递推式里的细节都在复数运算里每层的衰减直接折算进复数波速自然产生衰减虚部不再需要后置的指数吸收项。调试时看 Z 的实部虚部只是吸收修正实部才是决定模态截断的主导量。层数建议不超过四层再多层的特征值求根几乎必然出现伪根工作效率不升反降。5.2 海底声学参数的敏感性排序与误差条估算给传播损失图算误差条时不必做蒙特卡洛全空间采样先做敏感性排序。以细砂海底为基准逐项扰动声速 c2 变化 ±20 m/s远场 TL 变化约 ±1.5 dB密度 ρ2 变化 ±10%TL 变化约 ±0.5 dB衰减 alpha_b 从 0.5 调到 0.8 dB/λ10 公里处 TL 额外多衰减 3 到 4 dB。这个排序解释了为什么声速反演优先级最高——它对传播损失图的影响几乎是全局的。提示如果只想快速估算误差带把 c2 和 alpha_b 两个参数的扰动组合成四种极端情况渲染四张传播损失图逐点取包络线即可总量成本比蒙特卡洛低两个数量级。5.3 画传播损失图前的三处检查参考值、动态范围与网格尺度最后一件事三个小检查能省下大量返工。第一TL 参考值确认你用的参考距离是 1 米还是 1 公里不同参考差出 60 dB色标范围满盘皆输。第二动态范围不要拿全局 min-max 设色标要把 vmin 定在声源近场半经验衰减公式 20·log10(r) αr 的估计值附近否则近场强信号压缩掉所有远场细节。第三网格尺度距离步长要小于 2π/k_rn 的 1/20深度步长要小于 H/Xmax 即最小模态垂直半波长的 1/10不然传播损失图上的条纹是网格混叠画出来的不是声场算出来的。这三处检查做完图才能交出去当工程判读依据。本文还有配套的精品资源点击获取