简介一套面向海洋物理研究者的Matlab工具箱源码包聚焦seawater工具箱的核心实现帮助理解海水密度、声速、盐度、温度、压力等参数的物理计算过程。压缩包共含41个文件其中40个为.m格式脚本覆盖密度、声速、盐度、位温、混合层深度、地转流等常用计算模块另附1个README说明文件整体体积仅56KB结构紧凑便于按需查阅和逐文件研读。目前已有679人学习使用。借助该源码可深入掌握EOS-80等国际海水方程的实现细节理清压力、温度、盐度三者的耦合关系并在此基础上修改或新增函数扩展适用于特定海区或研究场景的水文模型对于从事海洋工程、气候建模、海洋生态研究的读者也能从中借鉴数值稳定性和算法优化思路。1. 为什么海洋物理计算要自己啃源码从一枚 CTD 剖面说起海洋物理研究中CTD温盐深仪下放一次会带回数千组温度、电导率和压力读数但你要的不是这三条原始曲线而是位温、位密、声速、混合层深度这些派生量。以前大家习惯直接调某个桌面软件点几下导出可一旦涉及批量处理 200 个站位、或者要改换状态方程版本就会意识到手里没有一份能改的 seawater 源码工作是做不下去的。这篇笔记就是讲清楚一套常见做法拿着 seawater 这套海洋物理源码怎么在本地把密度、声速、位温算明白参数怎么标定哪些地方容易翻车以及最终怎么把它封装进自己的数据处理流程。适合正在做水文数据处理、写毕业论文或者给浮标数据做质量控制的人新手能跟着一步步跑通熟手可以跳过安装直接看避坑章节。2. 拆开 seawater 源码状态方程、温盐深换算与声速这三个核心模块2.1 海水密度不是查表查出来的状态方程的迭代逻辑很多人第一次打开 seawater 源码时会下意识去找一个查表函数但实际上海水密度计算是一个状态方程迭代过程。源码里最核心的是密度模块它输入盐度、温度和压力输出密度公式采用国际常用的 EOS-80 标准。EOS-80 把海水看作一个压缩流体密度由纯水密度、盐度贡献和压力压缩项三部分叠加其中纯水密度本身就是一个关于温度的多项式源码里一般写成一个独立的函数比如dens0。真正的迭代逻辑在深水计算里压力影响密度密度又影响压力梯度所以当你要算某一深度的原位密度时不能直接带一个静压值进去而是要分层积分。源码里通常提供两种密度接口一种是直接给盐度、温度、压力返回原位密度另一种是给盐度、温度、深度返回密度——后者内部会先把深度转成压力再调用前者。这地方就是很多误用的根源后面避坑章节会专门讲。参数说明上盐度单位是 PSU实用盐度1967 或 PSS-78 定义温度是摄氏度压力是 dbar分巴注意不是 bar更不是 Pa。dbar 和深度的关系大约是 1 dbar ≈ 1 米水深所以在近海浅层很多人偷懒直接拿深度数值当 dbar 用这是可以接受的工程近似但在深水会偏差很大。源码里压力参数如果传错了单位密度第三位小数就开始错声速偏差可达每秒几米。2.2 盐度、电导率、位温你可能用错的那组换算函数seawater 源码里不只有密度还有一组专门做盐度换算的函数。CTD 实测的是电导率不是盐度所以第一道换算就是由电导率比到实用盐度。这个换算关系在源码里不是简单的多项式而是分段函数并且依赖温度和压力因为电导率对温度极其敏感温差 1 度可能让盐度偏差 0.01 以上。源码中典型做法是传入电导率比其实就是实测电导率除以标准 KCl 溶液电导率、温度和压力最后返回 PSU。位温的概念容易被新手忽略。海水的绝热温度梯度意味着一个水团从深处抬升到海面时温度会因为膨胀而降低如果你拿原位温度去做水团分析会以为两处水温度不同其实它们可能来自同一个源。源码里的位温函数就是做这个事输入原位温度、盐度、压力通过绝热温度梯度做积分输出参考压力通常是 0 dbar下的位温。注意这里又出现压力单位和密度模块一致都要用 dbar。实际使用中还有一个容易混淆的有些源码版本会把盐度换算函数命名为salinity有些叫salt返回的还有可能是绝对盐度而非实用盐度这直接关系到 TEOS-10 和 EOS-80 两个标准的选择。你手里的 CTD 数据如果是近十年的厂家通常已经按 TEOS-10 输出绝对盐度了但许多开源 seawater 源码还停留在 EOS-80这时候不能直接拿绝对盐度去喂给旧的密度公式得做换算或者改用新版源码。2.3 声速剖面计算从 UNESCO 公式到源码实现声速剖面是水声通信、声呐性能预报里最基础的环境输入也是海洋物理源码里实现得最直接的部分。常用的 UNESCO 声速经验公式是一个三项多项式第一项只跟温度有关第二项是盐度的线性项第三项是压力项里面还嵌套温度和压力的交叉项。源码里通常就一个函数入参是盐度、温度、压力出参是声速单位米每秒。这个公式本身并不复杂但工程上有一个坑声速公式分好几个版本有的适用温度范围是 -2 到 35 度盐度 0 到 42 PSU压力 0 到 10000 dbar超出范围公式照样能算但结果已经不可信。源码里一般只做很少的边界检查所以如果你用极地海底高温热液口的数据算出来的声速可能比实际偏大十几米每秒。这里建议不要迷信源码给出的值要结合实测声速仪校验或者至少做范围限制。另外有些 seawater 源码还会把声速的梯度计算也放进来目的是做射线追踪。如果你要做声线弯曲分析就得检查源码里是否提供了对温度和压力的一阶偏导函数。没有的话只能用数值差分自己算这涉及步长的选择温度步长 0.1 度、压力步长 1 dbar 通常就够太小会被浮点误差吃掉。3. 在本地跑通 seawater 源码最小调用示例与参数标定3.1 安装与依赖一份能跑的最小环境常见做法是把 seawater 源码包放进自己的项目目录而不是全局安装这样你能直接改源码看内部逻辑也方便固定某个版本。我一般会建一个干净的虚拟环境只装 numpy 和 scipy因为 seawater 源码里插值部分通常依赖scipy.interpolate你如果不想装 scipy也可以把插值函数替换成自己写的一维线性插值但没必要scipy 很常见。先确认源码目录结构。典型的 seawater 源码会有这些文件seawater/ __init__.py density.py salinity.py sound.py temperature.py constants.py把项目目录放到你的工程根目录下然后用相对导入就能用。下面是最小验证代码import sys sys.path.append(./) # 确保项目根目录在模块搜索路径里 import seawater as sw # 导入源码包 # 给一个标准海水参数 s 35.0 # 盐度PSU t 20.0 # 温度摄氏度 p 100.0 # 压力dbar rho sw.dens(s, t, p) print(密度: {:.4f} kg/m^3.format(rho))这段代码的逻辑是调用dens函数传入标量参数得到的密度大概在 1021 kg/m³ 左右。注意dens内部可能还需要一个函数来换算深度到压力你不需要手动调。参数说明盐度 35 是海洋学标准海水盐度温度 20 是常见表层温度压力 100 dbar 约相当于 100 米水深。如果你看到返回的密度小于 1000说明盐度或温度单位传错了。3.2 用源码计算标准海水密度10 行代码复现教科书记录值为了确认源码没问题我会用 UNESCO 教科书上的经典值来做冒烟测试。一个著名的参考点是盐度 35温度 20压力 0 时密度应该是 1024.77 kg/m³ 附近具体看用的是哪个状态方程EOS-80 的值是 1024.77。下面用一个二维网格来计算顺便验证数组输入是否正常import numpy as np import seawater as sw # 构造一个温度范围 0~30盐度范围 30~40 的网格 temp np.array([0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0]) salt np.array([30.0, 32.0, 34.0, 35.0, 36.0, 38.0, 40.0]) # 广播成二维网格 T, S np.meshgrid(temp, salt) # 计算海面压力下的密度 rho0 sw.dens(S, T, 0.0) # 找盐度35、温度20的那个值 ref rho0[3, 4] # 盐度36的行温度20的列 print(盐度35温度20密度:, ref)这里有一个非常关键的工程细节sw.dens函数内部如果支持 numpy 数组它会把数组当成逐元素计算返回同维度数组。但有些源码版本用了if type float之类判断数组进去会报错或者按列表处理所以这一步能帮你测出版本行为。上面代码里我在盐度数组中故意没有写 35因为 35 恰好落在索引 3 的位置30,32,34,35所以rho0[3,4]就是盐度 35、温度 20 的密度。实际输出应该在 1024.7 这个量级。如果差得很远优先检查温标很多老源码默认用 IPTS-68而你输入的是 ITS-90 温度两者相差 0.005 度对密度影响可忽略但对声速会有可察觉的影响。3.3 参数单位陷阱PSS-78、EOS-80 与 TEOS-10 的边界这是最容易踩坑的地方。seawater 旧版源码基于 EOS-80输入盐度是 PSS-78 实用盐度定义域是 2~42 PSU。但近年来 CTD 默认输出的是 TEOS-10 绝对盐度绝对盐度数值比实用盐度通常大 0.01~0.02在河口还会差更多。如果你直接把 CTD 导出的绝对盐度塞进旧源码密度误差约 0.01 kg/m³声速误差约 0.01 m/s单独看可以忽略但在做浮力频率Brunt-Väisälä 频率计算时密度微弱偏差会被求导放大十倍以上。我建议在使用任何 seawater 源码之前先核对一下输入数据是由什么标准导出的。海鸟 CTD 的软件里可以设置输出为 Practical Salinity 还是 Absolute Salinity如果你不确定可以看数据文件头。如果是绝对值先做减法转换# 绝对盐度 SA 到实用盐度 SP 的粗算非官方公式仅用于调试 SA 35.0123 SP (SA - 0.01) / (1 - 0.005) # 简易线性近似切勿用于正式研究更好的做法是直接用支持 TEOS-10 的gsw库做交叉验证这个后面会讲。压力单位更是重灾区。看源码里的常量表P_REF 0.0是参考压力单位是 dbar。有些函数名带p的需要的是压力而不是深度比如sw.dens(s, t, p)第三个参数是压力。如果你手头只有深度需要把深度转压力常见近似是p depth / 0.1即深度米数乘以 10 得到 dbar不对1 dbar 约等于 1 米水深所以深度 100 米压力约 100 dbar因此p depth就行但严格要加上大气压约 10 dbar也就是p depth 10。深水时这个差异不明显但表层的静压与深度只差 10 dbar影响密度第四位小数声速约 0.2 m/s。所以如果做高精度声速预报压力必须用depth * 1.004 10之类更精细的换算具体系数看源码里有没有提供pressure函数。4. 源码调试与验证用历史 CTD 数据反推计算精度的三条路径4.1 拿 WOA 数据集做基准比对拿到一套不是自己写的 seawater 源码第一件事不是相信它而是验证。最省力的基准数据是全球海洋 atlas 类的格点数据比如 WOAWorld Ocean Atlas它给出了各个标准层的温盐场以及根据 EOS-80 计算好的密度、声速等量。你可以从官方下载一个 0.25 度格点的剖面然后截取你需要的经纬度位置把温盐输入你的 seawater 源码对比它算出的密度与 WOA 自带的密度。具体做法是先读取 WOA 的 netCDF 文件取出温度、盐度、压力WOA 给的是深度需要转换为压力然后逐层调用源码函数算完之后与数据文件里的density变量做相减统计残差。残差的合理范围是 ±0.03 kg/m³如果超过这个量级说明你的输入单位或标准不对。这一步一定要写在自动测试脚本里以后每次改源码都能跑一遍回归测试。下面是个伪代码结构的测试示意import netCDF4 as nc import numpy as np import seawater as sw ds nc.Dataset(woa18_grid.nc) t ds.variables[temperature][0, 20, 10, :] # 某站位的温度剖面 s ds.variables[salinity][0, 20, 10, :] depth ds.variables[depth][:] # 深度转压力 p depth 10.0 # 约等于 dbar rho_calc sw.dens(s, t, p) rho_woa ds.variables[density][0, 20, 10, :] diff rho_calc - rho_woa print(最大偏差:, np.max(np.abs(diff)))如果diff在浅层有明显的锯齿状大概率是 WOA 的深度网格不均匀而你的源码内部做了插值插值方法引入误差。这时要检查源码里密度计算是否对压力做了线性插值如果只是简单的一维插值建议改用更高阶插值或者直接把 WOA 标准层的压力作为输入不自行插值。4.2 敏感性分析盐度 ±0.01 对密度和声速的影响在把源码交给别人或用于业务前我习惯做一次敏感性分析搞清楚输入误差的传播。这个目标不是质疑源码而是帮自己确定数据的精度要求。比如目标是要让声速误差小于 0.5 m/s那么需要温盐的精确度达到多少以下代码用有限差分计算敏感性代码量很小但信息量很大import seawater as sw s 35.0 t 20.0 p 100.0 # 盐度扰动 0.01 PSU ds 0.01 rho_base sw.dens(s, t, p) rho_up sw.dens(s ds, t, p) sound_base sw.soundvel(s, t, p) sound_up sw.soundvel(s ds, t, p) print(盐度0.01 密度变化:, rho_up - rho_base) print(盐度0.01 声速变化:, sound_up - sound_base)参数说明这里之所以把盐度扰动设在 0.01因为这是大多数 CTD 在实际海上的盐度精度极限。你会发现密度变化大约在千分之几 kg/m³声速变化约 0.06 m/s。个位数的盐度误差比如 0.1 PSU会导致声速 0.6 m/s 的偏差这对水声定位来说已经不可接受了。同理可以做温度扰动 0.001 度的敏感性。这个结果能指导你在数据质控时对温度、盐度的异常值筛选阈值做合理化设置而不是拍脑袋。4.3 与国外开源库 gsw 的交叉验证如果你的进程里只依赖一个 seawater 源码出了错很难判断是算法问题还是参数问题所以我会用gswGibbs SeaWater库作为独立参照。gsw 是 TEOS-10 标准库函数接口和 seawater 完全不同但计算物理量相同两者交叉验证能互相暴露问题。用 gsw 计算同一温盐压下的密度和声速import gsw SA 35.0 # 绝对盐度 CT 20.0 # 保守温度近似等于位温 p 100.0 # 压力 dbar lon 120.0 lat 30.0 # 绝对盐度到实用盐度 SP gsw.SP_from_SA(SA, p, lon, lat) rho gsw.rho(SA, CT, p) sound gsw.sound_speed(SA, CT, p) print(gsw密度:, rho, 声速:, sound)注意 gsw 中的温度参数是保守温度Conservative Temperature和 seawater 里的原位温度差一个量级很小的绝热效应。在做交叉验证时你要先用gsw.CT_from_t(SA, t, p)把原位温度转换成保守温度否则会引入压力项误差。理想情况下同一个剖面两个库算出的声速差别应小于 0.2 m/s。如果差别大检查 seawater 源码里是否用了旧版声速公式有的源码还在用 Wilson 1960 公式或者温标不一致。这一步很值得自动化但别把 gsw 的依赖写进生产程序只是在调试阶段用因为 gsw 的包很大而且启动慢。5. 常见问题与避坑seawater 源码使用中容易翻车的五个细节5.1 输入盐度不是 PSU 而是绝对盐度为什么结果偏了现象用 CTD 实测数据算密度结果总是比海区历史值低 0.03~0.05 kg/m³而且越往深海偏差越大。原因CTD 后处理软件默认输出绝对盐度而 seawater 源码按实用盐度设计。绝对盐度的数值偏高尤其在大洋深层两者差可达 0.05 PSU 量级。解决在数据读取模块中明确区分两个字段。如果数据文件里带有flag参数优先看图没有的话用 gsw 库把绝对盐度转成实用盐度或者用一条非常见盐度的水样校准。注意不要直接对所有数据做常数减法因为转换关系跟区域有关系不同海域的构成离子比例不同偏差并非恒定。5.2 压力参数用的是 dbar 还是 bar一个单位引发的 10% 偏差现象声速计算结果整体偏大或偏小且随深度线性增长到 5000 米处偏差超过 20 m/s。原因源码里的压力参数单位是 dbar但有人把 CTD 输出的压力单位通常是 dbar误以为是 bar或者反过来把 bar 数据直接喂进去。1 bar 1000 dbar但压力数值在 0~1000 范围内时两者差一个量级。解决写一个强制校验读入压力数组后先打印它的取值范围。如果最大压力超过 5000而最大深度只有 1000 米那单位一定有问题。另外注意源码里pressure函数可能会对压力做加 10 dbar 处理你要确认你的输入是否已经包含大气压避免重复加。5.3 数组输入与标量输入的维度陷阱现象单点计算正确但传入一维数组时结果变成 NaN或者返回一个嵌套列表。原因源码内部函数有的用math.sqrt而不是numpy.sqrtmath只接受标量遇到数组就抛出 TypeError有的在计算中用了if判断数组整体导致布尔歧义。解决在使用源码前先把测试输入转成 numpy 数组再传或者检查源码里是否写了np.asarray。如果你不想大改源码可以在调用处统一加一个包装函数把所有输入转成 float 或者按元素循环计算这样牺牲一点性能但绝对稳定。对于要批处理的剖面循环 5000 个点也就几秒钟完全能接受。5.4 温度是 ITS-90 还是 IPTS-68老版本源码的温标兼容问题现象和同事的结果对比声速总是差 0.3 m/s且不随深度变化像个系统偏移。原因很多老 seaswater 源码采用的是 IPTS-68 温标而现代 CTD 输出的温度是 ITS-90。两个温标之间相差约 0.005~0.01 度虽然不大但声速对温度很敏感0.01 度足以造成 0.03 m/s 的偏差少数公式更敏感。解决在源码的配置区或常量文件里查看是否有TEMP_SCALE字段。如果没有就要用ITS90 IPTS68 - 0.0003 * IPTS68这种近似换算注意这是经验式不同温度段系数略有差别。更稳妥的做法是升级到支持 ITS-90 的版本或者用 gsw 库替代。判断是温标问题的最简单方法把输入温度整体加 0.01 度看结果是否逼近你的参考值。5.5 源码里用魔数改参数时别动这几个常量现象想调整源码适应高盐环境直接改了密度公式里的某个系数结果整个海域剖面全部失真。原因seawater 源码中会有一批没有注释的魔数比如密度公式里的4.8314e-4、声速公式里的1.389等这些不是随便的比例系数而是从国际标准回归出来的牵一发动全身。有人调试时看到某个数字不顺眼就改了导致多项交叉检验全失败。解决不要修改源码里的基本常量。如果确实需要适配特定海域应该做的是在调用端做数据校正而不是改算法。源码里的常量表有什么用它至少能告诉你公式的地球参数设定比如重力加速度g 9.7803用的是 45 度纬度海平面值如果你在北极圈做计算重力差异对位密度的影响约 0.005 kg/m³需要单独修正但改常量会让别人无法复核。正确做法是在自己的代码里覆盖g参数而不是改原包。6. 把源码接进自己的水文处理流程封装成函数并验证输出6.1 封装一个 seawater_props 函数返回密度、声速、位温单独调用散落的源码函数不利于维护我会把它们统一封装成一个数据类或者字典返回式函数这样业务代码只需要引入一个入口。下面是建议的封装模板import numpy as np import seawater as sw def seawater_props(sp, t, p): 输入: 盐度 sp (PSU), 温度 t (ITS-90, °C), 压力 p (dbar) 输出: 密度、声速、位温 t_68 t - 0.0003 * t # ITS-90 转 IPTS-68供老源码用 rho sw.dens(sp, t_68, p) sound sw.soundvel(sp, t_68, p) pt sw.ptmp(sp, t_68, p) # 位温参考压力 0 dbar return { density: rho, sound_speed: sound, potential_temperature: pt } # 批量处理一条剖面 depths np.arange(0, 1000, 10) s_profile np.full_like(depths, 35.0) t_profile 20.0 - 0.01 * depths result [seawater_props(s, t, p) for s, t, p in zip(s_profile, t_profile, depths)] # 取声速列 sound_curve [r[sound_speed] for r in result]逻辑说明这里刻意做了温标换算因为老源码默认 IPTS-68。np.full_like用来生成盐度剖面很方便。输出用字典业务代码可以按需取键避免记错函数名。这个封装函数后续如果要换用新标准只需要改内部实现调用方完全不用动。参数说明如果你的数据是深度而不是压力可以在封装函数内部加一个depth_to_pressure的处理但建议在数据准备阶段就转好不要让封装函数承担过多职责。6.2 批量处理 CTD 剖面时的性能与内存注意当你有 200 个站位、每个站位 5000 个深度层时用上面那个列表推导式会生成大量小字典内存开销不小。更稳的做法是用 numpy 向量化但要先确认源码函数本身支持向量化。如果不支持可以用numpy.frompyfunc把它包装成能接受数组的函数import numpy as np import seawater as sw # 将标量函数包装成 ufunc vector_dens np.frompyfunc(sw.dens, 3, 1) # 准备整条剖面的数组 s_arr np.full(5000, 35.0) t_arr np.linspace(20, 2, 5000) p_arr np.linspace(0, 100, 5000) rho_arr vector_dens(s_arr, t_arr, p_arr).astype(float)用frompyfunc的好处是保留 numpy 的广播机制也避免 Python 循环过慢。但要注意返回的是 object 数组所以要.astype(float)。实测 5000 个点的剖面计算在这个方式下通常小于 0.1 秒瓶颈反而是 IO。另外如果源码内部用了scipy.interpolate它会自动处理向量化就不用这个包装。6.3 用科学计算容器固化运行环境源码包最容易被破坏的方式就是依赖库升级。我早期吃过这个亏跑完一组实验过两个月再跑发现结果和之前记录对不上查了半天才知道是 numpy 的行为变化影响了插值。后来我习惯在项目根目录放一个环境描述文件把 numpy、scipy 的版本锁死。更彻底一点用容器技术把整个运行时固化下来这样换机器也不用重新折腾。最后分享一个习惯每一批数据计算前我都会跑一次标准海水点温度 20、盐度 35、压力 0作为冒烟测试如果这个点的密度不在 1024.7 到 1024.8 范围内说明环境或源码已经被改坏我会立刻停止任务去排查而不是继续处理。这个习惯帮我省下过很多次返工。希望帮到你。本文还有配套的精品资源点击获取