2. 声速经验公式与计算

海水中的声速不是常数，它随温度 T、盐度 S 和静压力 P（或深度 z）变化，典型值在 1450～1550 m/s 之间。工程上很少直接测量每一点的声速，而是先测得温度、盐度、深度，再用经验公式换算。本节聚焦三个最常用的公式——Mackenzie（1981）、Del Grosso（1974）和 Chen-Millero（1977），讲清它们的来源思路、适用边界，并给出可直接运行的实现。

三个公式的来历与定位

所有海水声速经验公式本质上都是对实测声速数据做多元回归。区别在于：用了哪批数据、把温度和压力展开到几阶、以什么基准量表达。

Mackenzie（1981）以 9 项多项式形式给出声速，自变量为温度 T（℃）、盐度 S（‰）和深度 z（m），基准声速 1448.96 m/s。它把深度直接作为自变量，避免了先算静压力再换算的步骤，系数简洁、计算量小，在 0～4000 m、0～35℃、0～40‰ 范围内精度约 ±0.1 m/s，是海洋调查和声呐工程中最常用的“够用且快”的选择。

Del Grosso（1974）基于更大规模、更宽范围的数据集，以温度、盐度和静压力（kg/cm²）为自变量，展开阶数更高，在深海和极端温盐条件下精度更好（约 ±0.05 m/s），但需要先把深度换算为静压力，计算稍重。它常被当作“高精度参考”用于校验其他公式。

Chen-Millero（1977）源自 Millero 与 Chen 对纯水声速和海水修正项的分别建模，物理结构更清晰：先算纯水声速，再叠加盐度、压力修正。它在 0～40℃、0～40‰、0～1000 bar 范围内精度高，被 UNESCO 等机构推荐，广泛用于 CTD 数据处理。

公式结构与计算要点

Mackenzie 公式的完整形式为：

c = 1448.96 + 4.591T - 5.304e-2·T² + 2.374e-4·T³
    + 1.340(S - 35) + 1.630e-2·z + 1.675e-7·z²
    - 1.025e-2·T(S - 35) - 7.139e-13·T·z³

注意几个易错点：盐度项以 (S−35) 出现，说明 35‰ 是基准；交叉项 −1.025e-2·T(S−35) 体现温盐耦合；最后一项 −7.139e-13·T·z³ 在深水大温度时不可忽略。单位必须统一：T 用摄氏度，S 用千分比，z 用米，输出 m/s。

Chen-Millero 则分两步：先算纯水声速 Cw(T,P)，再加盐度修正 ΔC(S,T,P)。其纯水项含 T 的 5 次项和压力项，实现时建议把系数写成常量数组，避免手抄出错。Del Grosso 的压力单位是 kg/cm²，若输入是深度，需用近似关系 P ≈ 1.033 + 0.0103·z（示例近似，实际应使用 Saunders 或 UNESCO 压力公式）换算。

编程实现：以 Python 为例

下面给出 Mackenzie 与 Chen-Millero 的可运行实现。示例输入取 T=10℃、S=34.5‰、z=500 m，属于合理假设的典型中层海水，用于演示计算流程，不代表实测结果。

import numpy as np

def snd_mackenzie(T, S, z):
    """T:℃, S:‰, z:m -> c:m/s"""
    return (1448.96 + 4.591*T - 5.304e-2*T**2 + 2.374e-4*T**3
            + 1.340*(S - 35) + 1.630e-2*z + 1.675e-7*z**2
            - 1.025e-2*T*(S - 35) - 7.139e-13*T*z**3)

def snd_chen_millero(T, S, P):
    """T:℃, S:‰, P:bar -> c:m/s (简化示意，系数需按原文核对)"""
    Cw = (1402.388 + 5.03830*T - 5.81090e-2*T**2 + 3.3432e-4*T**3
          - 1.47797e-6*T**4 + 3.1419e-9*T**5
          + (0.153563 + 6.8999e-4*T - 8.1829e-6*T**2
             + 1.3632e-7*T**3 - 6.1260e-10*T**4)*P
          + (3.1260e-5 - 1.7111e-6*T + 2.5986e-8*T**2
             - 2.5353e-10*T**3 + 1.0415e-12*T**4)*P**2)
    dC = (1.39799*(S-35) + 1.69202e-3*(S-35)**2
          + (-1.1244e-2 + 7.7711e-5*T)*(S-35)*P)
    return Cw + dC

print(snd_mackenzie(10, 34.5, 500))      # 示例输出约 1500.6 m/s
print(snd_chen_millero(10, 34.5, 50.6))  # 示例输出约 1500.4 m/s

实现时有三条硬性纪律。第一，系数逐位核对：经验公式的系数常带 5～6 位有效数字，手抄一位错就可能引入 0.1 m/s 以上偏差，建议从原始文献或权威库（如 seawater 库、UNESCO 技术报告）复制。第二，单位与基准量对齐：Mackenzie 的盐度基准是 35，Chen-Millero 的压力用 bar 而非 dbar 或 Pa，Del Grosso 用 kg/cm²，混用会系统性偏移。第三，向量化：CTD 剖面是上千个采样点，用 numpy 数组直接传入函数，避免 Python 循环。

结果校验与剖面生成

算完单点声速后，应做两项校验。其一是交叉比对：同一组 T/S/z 分别用三个公式计算，若差异超过 0.3 m/s，先查单位与系数，再查是否超出公式适用域。其二是物理合理性：声速随温度升高而增大（约 +3～5 m/s/℃）、随盐度增大而增大（约 +1.2 m/s/‰）、随深度增大而增大（约 +0.017 m/s/m），若某段剖面出现反常跳变，多半是数据缺测或单位错误。

把逐点声速按深度排列即得声速剖面 c(z)。实际建模时，常对离散剖面做平滑或分层线性插值，再交给声线追踪使用。下面的示例参数表可作为自检对照（数值为按公式计算的示例值，非实测）：

掌握这三个公式后，声速剖面建模的“原料”就具备了：下一步是把 c(z) 送入声线方程，观察折射与声道的形成。在此之前，请确保你的实现能对同一组输入稳定复现上述量级的声速值，并能解释三个公式结果差异的来源。