简介seawater工具箱是Matlab环境下专为物理海洋学设计的函数库采用C语言核心计算与Matlab封装形式提供海水密度、声速、盐度、电导率等关键参数的高精度求解方案。资源共41个文件包含40个.m源码文件及1份readme说明压缩包仅56KB轻量易部署。其中sw_dens、sw_svel、sw_salt等核心函数覆盖海水状态方程、声速剖面、盐度换算等常用计算sw_dist、sw_gvel、sw_ptmp等辅助函数可完成数据插值、地转流分析与位温订正形成从参数求解到结果可视化的完整链路。源码严格遵循EOS-80等国际海水方程细致处理压力、温度、盐度三变量的耦合关系并注重数值稳定性与计算效率的平衡。已有679人学习使用研究者可参照readme快速上手并通过修改源码集成自定义物理模型适用于海洋工程、气候建模、海洋生态等多领域深度开发是理解海洋物理计算细节的实用参考。1. 海洋物理研究里的 seawater 源码它到底是什么为什么值得你读先给结论seawater 不是一个海洋数值模式也不是观测系统它是一套为海洋物理研究服务的开源计算库核心任务是把海洋学里最常用的状态量——盐度、温度、压力、密度、声速、位势高度——用统一的代码接口算出来并且能直接对接 CTD、Argo 浮标、走航观测这些实测数据。做海洋物理的人都知道海水状态方程里每个量都牵一发动全身密度算错千分之一地转流流速就能偏出好几厘米每秒这对研究大洋环流、水团分析、锋面识别来说是不可接受的误差。我说几个最典型的场景你就明白了用 Argo 剖面数据算位势高度异常、把 CTD 航次数据插值到标准层之后计算密度、用温盐深数据诊断水团是否发生混合。这些活儿如果自己从零写公式不仅要反复查 UNESCO 方程系数还要处理单位、纬度修正、压力单位这些容易被忽略的细节。seawater 这类库把这些公式固化成经过验证的函数你要做的就是理解每个参数的物理意义把数据喂进去然后检查输出量级是否合理。适合的读者也很明确正在处理海洋实测数据的硕士博士、做物理海洋或者海洋环境工程的一线从业者以及刚入门想用标准工具替代手写公式的开发者。下面我按实际工作中最常用到的计算路径把这份源码的用法和坑一次讲透。2. 从温盐深数据到密度场seawater 核心函数与调用约定2.1 源码里最常用的函数族密度、声速、位势高度的计算入口以常见的 seawater 实现Python 版本和 MATLAB 版本都遵循 UNESCO 1983 年方程为例日常使用频率最高的是这样几组函数盐度计算从电导率比出发salt或sw_salt接受电导率比和温压条件返回到实用盐度密度计算通常用density或sw_dens输入温度、盐度、压力后输出密度。声速剖面在声学定位和海洋探测中很常用对应sound或sw_svel。位势高度则依赖geostrophic或sw_geostr它把密度剖面沿垂向积分得到相对于参考面的动力高度。我一般建议把调用约定记成温度单位是 ITS-90、盐度单位是 PSU、压力单位是 dbar这一条主线。温度必须转换成摄氏温度绝对温标在接口里通常不被接受压力字段如果来自 CTD 原始数据往往是 db 为单位的海水静压而不是大气压加海水压的绝对压。如果你直接把压力值传进去而忘记去掉大气压密度在表层会有一个系统偏差这会直接污染后续地转流计算结果。import seawater as sw # 典型 CTD 单站数据: 温度(degC), 盐度(PSU), 压力(dbar) # 注意这里压力是 CTD 原始输出, 包含约 10 dbar 的大气压影响 temp [24.5, 22.3, 18.7, 14.2, 9.8] # 从表层往下 salt [35.2, 35.4, 35.6, 35.1, 34.8] pres [2.0, 50.0, 100.0, 200.0, 400.0] # 计算海水密度, 默认使用 UNESCO 1983 状态方程 dens sw.dens(salt, temp, pres) print(密度输出:, dens) # 计算声速剖面, 用于声学设备校准 svel sw.svel(salt, temp, pres) print(声速输出:, svel)这段代码背后的逻辑并不复杂sw.dens按 UNESCO 方程先计算某个温盐压状态下的比容再换算成密度返回值的单位是 kg/m³典型海水密度约在 1020 到 1030 之间。sw.svel则直接套用声速经验公式输出单位是 m/s约在 1450 到 1550 范围。如果你看到密度输出在 1000 以下或者声速低于 1400先不要怀疑函数写错大概率是温度盐度压力三个数组没有对齐比如压力传成了负值或者温度传成了开尔文。2.2 为什么密度计算要区分位势深度与实际深度seawater 的函数在设计上还有一个很多人会忽略的地方垂向坐标可以用实际深度也可以用位势深度。实际深度是仪器测出来的几何深度位势深度则把地球重力场的变化考虑进去用动力米表示。地转流计算和位势高度积分必须使用位势深度否则在深层大尺度计算里会累积误差。海水密度分布不均匀同样的几何深度在不同海域对应的位势面并不相同直接拿几何深度去做位势积分等于默认海水正压这在锋面强烈的区域是站不住脚的。import seawater as sw import numpy as np # 假设有连续剖面数据 z np.arange(0, 500, 10) # 实际深度, 单位米 lat 28.5 # 所在纬度, 度 # 把实际深度转换成位势深度, 用于位势高度/地转流计算 dz sw.z2p(z, lat) # 深度转压力(dbar), 同时考虑纬度重力修正 print(不同深度的压力值(dbar):, dz) # 再计算位势高度异常, 参考面取 500 dbar # 先构造盐度温度剖面, 这里用简单插值示例 temp_profile 25 - z * 0.03 salt_profile 35.5 z * 0.001 # 用位势压力坐标计算各层密度 dens_profile sw.dens(salt_profile, temp_profile, dz)sw.z2p做的事情是把几何深度转换成标准海水压力换算时引入纬度修正因子因为地球不是完美球体重力加速度随纬度变化。你会发现输出压力值和 z 不是简单的一一对应深层差异会到几十 dbar这就是位势修正的体现。之后再算位势高度时seawater 内部积分用的就是这套修正后的压力坐标。这里要强调一点不要自己写一个压力深度×1.02的粗略公式来替代尤其在深水区或者高纬度海区这个近似带来的误差会直接影响地转流方向判断。2.3 实测数据处理里必须做的单位对齐与纬度修正实测数据进入 seawater 之前单位对齐是翻车率最高的环节。Argo 浮标下载的 NetCDF 文件里温度通常是摄氏度但盐度字段有两种可能有的已经转换成实用盐度有的仍保留为电导率比或原始电导率CTD 数据中压力是 dbar但某些老式仪器以米为单位存储深度字段走航观测的盐度可能来自实验室盐度计分析数值是绝对盐度而不是实用盐度。这些都是我在实际处理数据时踩过的坑。import xarray as xr import seawater as sw # 读取 Argo 剖面 NetCDF 文件示例 ds xr.open_dataset(argo_profile.nc) temp ds.TEMP.values # 摄氏度, ITS-90 pres ds.PRES.values # dbar psal ds.PSAL.values # 实用盐度, 已经是 PSU # 如果文件中盐度字段名是 PSAL_QC 或原始电导率, 需要先检查单位 if PSAL in ds.variables: salt psal else: # 电导率比转为盐度需要温度和压力, 这是最容易被忽视的地方 salt sw.salt(ds.CNDC.values / 42.9, temp, pres) # 计算密度 dens sw.dens(salt, temp, pres) # 对压力做质量控制: 负压或超过仪器量程的值直接剔除 mask (pres 0) (pres 2000) dens_clean dens[mask]如果你遇到的是电导率比数据转换成实用盐度时的关键参数是电导率比的定义标准海水在 15 度、1 个大气压下的电导率是 42.9 mS/cm所以要把原始电导率除以 42.9 得到比值r然后才能调用sw.salt(r, temp, pres)。很多人在这里直接用原始电导率传进去导致盐度全部偏大再往下算密度就完全失真。此外纬度修正主要作用于压力转换和位势积分sw.dens本身不要求纬度参数但如果你调用的是地转流函数纬度就必须传对否则科氏参数就错了。3. 用 seawater 从零构建一个水团分析流程从 Argo 数据到 T-S 图3.1 数据准备多剖面数据如何组织成 seawater 可用的输入格式水团分析是物理海洋里最常用的诊断手段之一做法是把每个观测点的温盐数据画在 T-S 图上等密度线叠加上去就能判断水团的来源、混合程度和变化趋势。做这件事之前需要把 Argo 或 CTD 数据按照剖面号 深度层的结构组织好。常见做法是把所有剖面拉平成一个大表每条记录包含温度、盐度、压力、纬度、经度、时间。我建议用 Pandas 的 DataFrame 来组织好处是可以直接做掩膜筛选比如去掉压力小于 10 dbar 的表层数据避免日晒导致的温盐异常干扰水团分析。再把 DataFrame 里的列直接转成 numpy 数组传入 seawater。注意不要把 DataFrame 本身传给 seawater 函数大部分接口只接受 numpy 数组或列表传入 DataFrame 会触发类型异常而且输出结果不再保留索引后续回溯剖面会麻烦。import pandas as pd import numpy as np import seawater as sw # 构造示例: 三个 Argo 剖面的温盐数据 data { profile_id: [1, 1, 1, 2, 2, 2, 3, 3, 3], pressure: [10, 100, 300, 15, 80, 250, 5, 120, 400], temperature: [26.1, 19.5, 11.2, 25.8, 18.9, 10.5, 27.0, 17.8, 8.9], salinity: [35.1, 35.6, 34.9, 35.0, 35.5, 34.8, 34.9, 35.3, 34.6], latitude: [28.0, 28.0, 28.0, 29.5, 29.5, 29.5, 27.2, 27.2, 27.2], } df pd.DataFrame(data) # 去掉表层日晒影响的浅层数据 df df[df[pressure] 10] # 转成数组传给 seawater salt df[salinity].values.astype(np.float64) temp df[temperature].values.astype(np.float64) pres df[pressure].values.astype(np.float64) # 计算每条记录的密度 dens sw.dens(salt, temp, pres) df[density] dens # 查看密度范围, 验证是否处于海水合理区间 print(df[[profile_id, pressure, density]].describe())这段代码的逻辑是先把数据规整成每一行是一条观测记录的长表格式然后按压力筛选再批量计算密度。astype(np.float64)这一步不能省因为 Pandas 列可能有NaN或者整型seawater 在混入整型时会出现广播错误而且浮点精度不足会在深层密度积分时产生可见误差。输出密度描述统计之后你应该看到密度范围在 1020 到 1028 之间如果出现低于 1018 的异常值说明盐度或者温度数据有坏值需要在进入 T-S 分析前处理。3.2 T-S 图绘制与等密度线叠加的完整代码T-S 图的视觉核心是散点图叠加等密度线。等密度线不用自己解方程去画seawater 或配套工具提供了直接在图上生成等值线的辅助函数。标准做法是先生成一个盐度-温度网格在这个网格上计算密度再用contour画出密度等值线最后把实测的温盐散点叠加在图上。import matplotlib.pyplot as plt import numpy as np import seawater as sw # 生成盐度-温度网格, 覆盖亚热带水团范围 s_grid np.linspace(33.5, 36.5, 100) t_grid np.linspace(5, 30, 100) S, T np.meshgrid(s_grid, t_grid) # 固定压力为 0 (表层) 或取 200 dbar 代表次表层 p_fixed 200 RHO sw.dens(S, T, np.full_like(S, p_fixed)) plt.figure(figsize(8, 6)) # 画等密度线, 每 0.5 kg/m3 一条 cs plt.contour(S, T, RHO, levelsnp.arange(1020, 1029, 0.5), colorsgray, linewidths0.8) plt.clabel(cs, inlineTrue, fontsize8) # 叠加实测散点, 用剖面/站位做颜色区分 sc plt.scatter(df[salinity], df[temperature], cdf[profile_id], cmapviridis, s30, edgecolork) plt.colorbar(sc, labelprofile ID) plt.xlabel(Salinity (PSU)) plt.ylabel(Temperature (degC)) plt.title(T-S Diagram with Isopycnals at 200 dbar) plt.show()这里需要注意等密度线是在固定压力平面上画的严格来说应该用中性密度或者位势密度来对比不同深度的水团但大多数人做初步诊断时直接用某一压力层的等密度线作为参考。要画位势密度即把水团绝热移到参考面后的密度seawater 提供sw.pden它接受盐度、温度、压力、参考压力参数建议在 1000 dbar 参考面上使用这样能更好地区分深层水团。散点的颜色映射到剖面编号可以直观看出不同空间位置的水团是否落在同一条等密度线上如果不同颜色聚集在不同密度区间说明存在明显的水团边界。3.3 水团边界判断与混合诊断的常用参数把 T-S 图做出来之后真正有价值的判断是看散点走向。如果散点沿等密度线分布说明水团之间发生的是等密度混合温度和盐度呈补偿变化如果散点穿越等密度线说明存在越密混合可能是双扩散或强风搅拌导致的。seawater 提供的密度计算可以帮你量化混合强度比如计算每个散点到给定等密度线的垂直距离用这个距离作为混合强度指标。# 对每条观测记录, 计算它在密度空间里到参考水团的距离 ref_salt 35.0 ref_temp 18.0 ref_dens sw.dens(ref_salt, ref_temp, 200) # 计算每条记录的密度差 dens_diff df[density].values - ref_dens # 用密度差大小判断水团是否接近参考水团 df[dens_diff] dens_diff print(密度差小于 0.2 的记录数:, (np.abs(dens_diff) 0.2).sum()) print(密度差大于 0.5 的记录数:, (np.abs(dens_diff) 0.5).sum())这段代码把参考水团的密度算出来然后逐点计算差值。密度差绝对值小于 0.2 可以认为该观测点与参考水团性质接近大于 0.5 则说明水团性质已经明显不同。这个阈值不是固定的在分析不同海域时应该根据研究区密度层结强度调整——层结强的区域密度差阈值可以放小到 0.1层结弱的区域则放大到 0.5。我一般先看等密度线间隔和散点分布范围再决定阈值避免用统一标准导致漏判。4. 位势高度与地转流计算seawater 里最容易被误用的高阶功能4.1 为什么位势高度要从密度剖面积分而来地转流是物理海洋最核心的动力学量之一它假设流动在水平压力梯度与科氏力之间达到平衡。实际计算中我不能直接测压力梯度而是从密度场入手先算出位势高度动力高度再沿水平方向求梯度就得到地转流速。位势高度的计算本质上是对比容的垂向积分起点是参考面向上积分到目标深度。seawater 的sw.geostr函数正是为这个目的设计的它接受剖面深度、温度、盐度和纬度输出各层的位势高度异常值。但这里有两个前提条件必须满足一是剖面必须是等间距或者至少单调排列的深度序列二是参考面的选择要合理。理想的参考面是流速为零或近似为零的深层比如 1000 dbar 或 2000 dbar如果你用表层做参考面地转流会完全失真因为表层本身流速并不为零。import seawater as sw import numpy as np # 构造一个典型深海剖面 depth np.arange(0, 1200, 50) # 单位米, 每 50m 一层的规则网格 temp 20 * np.exp(-depth / 300) 4 # 模拟温度随深度衰减 salt np.full_like(depth, 35.0) # 简化, 假设盐度恒定 lat 30.0 # 把几何深度转换为压力(dbar) p sw.z2p(depth, lat) # 计算相对于 1200 dbar 参考面的位势高度异常 # 注意: 返回两个值, 第一个是各层位势高度异常, 第二个是各层深度 dyn_height, p_out sw.geostr(p, temp, salt, lat, axis0)sw.geostr的返回结构要特别注意dyn_height是位势高度异常序列不同版本的 seawater 接口可能还返回参考面信息或实际深度序列建议执行之后先打印形状再继续用。位势高度异常的单位是动力米数值通常在零点几到几之间深层接近零、表层正或负取决于相对于参考面的密度分布。如果把dyn_height当成绝对位势高度去画图会发现量级很小不对注意它本身就是异常要叠加参考面的绝对位势才完整。4.2 地转流速的差分计算和纬度适用范围得到位势高度异常场之后地转流由相邻两个站位之间的位势高度差决定。标准公式是地转流速正比于位势高度梯度除以科氏参数所以纬度决定科氏参数低纬度地区科氏参数小同一水位梯度产生的流速更大。在赤道附近科氏参数趋近于零地转近似不成立seawater 函数在纬度绝对值小于 5 度时会给出不合理的极大流速这时候需要切换成赤道波或边界层理论来分析。# 假设有两个相邻站位的位势高度异常数组 dh1, dh2 dh1 np.array([0.05, 0.12, 0.18, 0.21, 0.23]) dh2 np.array([-0.02, 0.04, 0.09, 0.13, 0.16]) # 站位间距 20 km dx 20_000 # 米 # 位势高度差 delta_h dh1 - dh2 # 科氏参数 f 2 * 7.2921e-5 * np.sin(np.radians(lat)) # 地转流速: v (g/f) * (dh/dx), g 约 9.8 m/s2 g 9.8 v_geostrophic (g / f) * (delta_h / dx)地转流速计算里的单位是一个经典陷阱位势高度单位是动力米动力米乘以 g 就还原成几何位势的米但在差分公式里 g/f 的系数已经包含了单位转换如果你把位势高度先乘 9.8 再差分流速会大 9.8 倍。实际工作中大部分人会直接用 seawater 自带的地转流函数内部已经处理好单位。但如果你需要半地理流或者斜压切变就绕不开手写差分式。还有一个容易被忽视的点站位间距不是直线距离而是沿等深线方向的间距计算前要用球面距离公式算清楚否则在地形复杂的陆架区会出现虚假强流。4.3 参考面选取的三种策略与对应误差参考面怎么选直接影响地转流绝对值。有人说选取越深越好但深层流速并非绝对为零只是相对于表层流较小。我常用三种策略固定深度参考面、基于密度层结的等密度参考面、以及中性密度参考面。固定深度最简单但在陆架区可能落到海底之下等密度参考面需要先计算密度层结的稳定度再挑一个密度梯度最小的层作为参考层这种做法在边界流区域更合理。# 策略一: 固定深度参考面 1000 m ref_depth_idx np.argmin(np.abs(depth - 1000)) dh_ref_1000 dyn_height - dyn_height[ref_depth_idx] # 策略二: 选取密度梯度最小的层作为参考面 rho sw.dens(salt, temp, p) rho_gradient np.abs(np.diff(rho)) # 密度梯度最小值对应的层, 加 1 是因为 diff 后长度少 1 min_grad_idx np.argmin(rho_gradient) 1 dh_ref_densgrad dyn_height - dyn_height[min_grad_idx] # 比较两种参考面得到的表层位势高度差异 print(固定参考面表层位势高度:, dh_ref_1000[0]) print(密度梯度最小参考面表层位势高度:, dh_ref_densgrad[0])这段代码展示了参考面选取对结果的直接影响。你会发现两种方法给出的表层位势高度有差距在陆架-深海过渡区域差距甚至可能超过 0.05 动力米换算成地转流就是每秒几厘米的差异。没有哪个参考面是绝对正确的关键是把参考面选取方式写进方法说明里并在对比不同时期或不同海域数据时保持统一。如果研究区存在强深层流比如南极绕极流或深水溢流固定深度参考面会引入系统性偏差这时候要用中性密度参考面或直接做逆方法。5. 常见坑与排查方法seawater 计算中反复翻车的五个点5.1 压力单位与大气压混用导致密度偏高现象计算结果密度普遍比同区域气候态数据高出 0.5 到 1 kg/m³。原因CTD 的原始压力包含了海面大气压大约 10 dbar计算时未扣除导致整体压力偏大深层尤其明显。解决如果数据来自 CTD 原始文件且压力没有经过海面校正应统一减去 10 dbar 或准确的海面大气压值Argo 浮标数据一般已经做好气压校正不需要再减。日常处理时要把压力字段的最小值打印出来检查如果最小压力在 10 附近说明很可能还没扣大气压。5.2 盐度字段实际是电导率比而不是实用盐度现象T-S 图上的散点整体偏离气候态等值线而且密度偏小或偏大没有规律。原因部分数据库中盐度字段存的是电导率比没有经过转换直接用这个值参与密度计算导致盐度量纲错误。解决查看数据文档确认字段定义如果是电导率比用sw.salt转换转换时注意温度必须用同一个剖面的实测温度不能用固定值。转换前后盐度值的差异通常在小数点后两位到三位之间但对密度的影响能到 0.1 kg/m³ 量级足以改变水团归类。5.3 深度与压力坐标混淆导致位势积分错位现象位势高度剖面出现不正常的锯齿形或者地转流方向与已知环流相反。原因把一个轴是深度米、另一个轴是压力 dbar 的两组数据混在一起传给sw.geostr函数内部默认输入为压力深度值被当成了压力。解决统一用sw.z2p把深度转换成压力或者直接用 CTD 输出的压力数组不要混用。判断方法很简单打印输入的p数组最大值如果超过 2000 而研究区水深只有几百米多半是把深度当成压力了。5.4 参考面选取不同导致地转流结果无法重复现象同一批数据在不同文章或报告里地转流方向一致但流速差异很大互相无法验证。原因有人用 1000 dbar 做参考面有人用 2000 dbar还有人用海底最深层导致绝对流速不同。解决在数据处理流程里把参考面定义为变量并输出到报告里方便复现。更规范的做法是同时输出不同参考面下的结果供下游研究者和读者对比如果参考面深度处的观测有缺失要用插值补齐而不是直接填 0否则会产生虚假的压力梯度。5.5 剖面数据存在倒转或非单调深度序列现象计算密度时程序不报错但结果出现莫名其妙的负值或者密度在密度层结图上交叉。原因剖面数据没有按压力单调排序或者存在重复站位混入。解决在传给 seawater 前用np.argsort排序同时检查是否存在重复剖面号混合多剖面数据时更要谨慎。排序前要保留原始索引方便结果回溯。这个坑最隐蔽因为错的不是公式而是输入顺序且只影响部分层位肉眼不容易发现。6. 进阶用法把 seawater 接入业务化数据处理流水线如果只是做单次研究分析直接调函数就够了。但实际工程里我通常会面临这样的场景每个月有新的 Argo 或船测数据进来需要自动完成质控—标准层插值—密度/位势高度计算—T-S 图更新—异常报警这一整套流程。这个需求催生了把 seawater 封装进自动化流水线的做法。我的习惯是写一个数据处理的 Pipeline 类把前面说的所有坑用代码卡死在入口处这样后面的人接手不会再把电导率比当盐度用。import seawater as sw import numpy as np import pandas as pd from dataclasses import dataclass dataclass class SeawaterProfile: pressure: np.ndarray temperature: np.ndarray salinity: np.ndarray latitude: float def validate(self): # 入口质控: 单位检查与排序 assert self.pressure.min() 0, 压力存在负值 idx np.argsort(self.pressure) self.pressure self.pressure[idx] self.temperature self.temperature[idx] self.salinity self.salinity[idx] return self def density(self): return sw.dens(self.salinity, self.temperature, self.pressure) def dynamic_height(self, ref_pressure1000): ref_idx np.argmin(np.abs(self.pressure - ref_pressure)) dh, _ sw.geostr(self.pressure, self.temperature, self.salinity, self.latitude) return dh - dh[ref_idx]这个 Pipeline 类把单位检查和排序放在构造函数之后、正式计算之前保证下游函数永远不会吃到脏数据。dynamic_height方法重定义了参考面外部调用时可以灵活指定。实际流水线里我会再加一个验证步骤把计算结果与 WOA 气候态数据做对比如果密度偏离超过 1.5 个标准差就触发警告通过邮件或消息推送给数据管理员。这样遇到传感器漂移或者生物附着影响盐度时不会在几周之后才发现数据质量问题。最后分享一个测量习惯每次用 seawater 算完一批数据我会把典型剖面的密度最大值、位势高度表层值、最大地转流速三个指标存成一个小的 JSON 日志文件积累几周之后回头看趋势。如果某个指标突然跳变往往是输入数据有问题而不是代码问题。这种方法不依赖任何复杂平台只用标准库就能实现但对从业者来说省去了大量回头排查的体力活。希望你也能在自己的数据处理链条里建立类似的复核机制——花十分钟加一层检查能省掉后面几天重新算数据的代价这是我在这个方向上最想强调的一条经验。本文还有配套的精品资源点击获取