简介:Matlab seawater工具包是一套专为海洋科学计算设计的扩展库,面向海洋科研人员、工程师及高校相关专业学生,用于在Matlab环境中精确计算海水物理性质与进行海洋动力学模拟。包内含40个m文件,整体压缩包约55KB,涵盖密度、声速、比热容、焓值、盐度-电导率转换、浮力频率、压力转换等核心函数,并附标准化与单位换算辅助工具,便于直接调用与二次开发。工具基于Tom McDougall和John Taylor算法实现,支持海水状态方程、潜热计算等底层物理过程,能够大幅简化海洋建模与数据分析中的复杂数学处理。目前已有1615人学习下载,适合需要处理海水实测数据、开展海洋声学或热力学研究、进行海洋环境数值模拟的中高级用户,是海洋物理计算中实用且轻量的基础工具包。 第一次意识到seawater工具包的价值,是在处理一组南海CTD剖面数据时。当时我要从原始的温盐深记录里算出位温和密度来做水团分析,手写状态方程算出来的结果总是跟文献对不上,后来才发现问题出在温度尺度和压力单位上。换成Matlab的seawater工具包之后,折腾了两天的计算变成了十几行代码,结果跟国际公开数据集完全一致。如果你也经常跟CTD、Argo浮标或海洋模式输出数据打交道,这个工具包几乎可以算作海洋学数据处理的标配。
1. 从CTD原始数据到科学变量:seawater为什么是海洋数据处理的标配
1.1 CTD到底测了什么:电导率、温度、压力与"盐度"的换算关系
很多人第一次接触CTD(Conductivity-Temperature-Depth,温盐深仪)时,会误以为盐度是直接测出来的。实际上CTD最核心的原始输出是电导率、温度和压力,其中电阻率传感器的电导率测量值需要经过换算才能得到实用盐度(Practical Salinity)。这个换算涉及PSS-78实用盐度标度,包含复杂的多项式拟合,手算不现实,正是seawater这类工具包存在的理由。
温度也不是我们平时说的"水温"那么简单。CTD探头测的是原位温度(in-situ temperature),但在物理海洋学分析中,更多时候需要的是位温(potential temperature),即水团绝热移动到海面或某个参考压力后应该具有的温度。深海4000米处原位温度可能是2°C,位温却可能接近-1°C,两者差值随压力增大而增大,这在深层水团分析中是不能忽略的。
压力同样有讲究。CTD输出的压力单位是分巴(dbar),1 dbar约等于1米水深,但严格来说这个换算依赖于纬度(重力加速度差异)。seawater中的sw_pres函数专门处理深度到压力的转换,输入深度和纬度就能得到准确的分巴值,省去了手工查表或者用近似公式的麻烦。
1.2 seawater工具包的来龙去脉:EOS-80标准与CSIRO的实现
seawater工具包由澳大利亚联邦科学与工业研究组织(CSIRO)海洋研究室的Phil Morgan等人开发,核心是实现了UNESCO 1980年海水状态方程(简称EOS-80)。这套方程是20世纪80年代以来物理海洋学计算密度、声速、比容等热力学性质的事实标准,绝大多数海洋学教材和经典文献中的数据都基于它。
工具包从2.0版本发展到现在的3.3.x系列,全部是纯MATLAB m文件实现,没有编译依赖,没有工具箱依赖,复制到本地就能用。3.x版本的一个重要更新是把温度基准从IPTS-68国际实用温标切换到了ITS-90国际温标,这使得它和现代CTD传感器的输出标准保持一致。老版本2.0在Matlab R2016a之后的版本上会出现一些警告提示,但核心计算函数仍然可以运行,只是温度尺度需要手动修正。
整套工具包的函数命名非常有规律,全部以sw_开头,后面跟着物理量的缩写,比如sw_ptmp表示potential temperature(位温),sw_dens表示density(密度),sw_vel表示sound velocity(声速)。这种命名约定让使用者不需要记忆大量散乱的函数名,看到函数名就能猜到用途。
1.3 常用函数一览:一张表看清sw_家族
以下是我在实际数据处理中最常用到的函数,它们覆盖了CTD数据后处理90%以上的需求。
| 函数 | 功能 | 典型调用 |
|---|---|---|
| sw_ptmp | 位温计算 | sw_ptmp(S,T,P,PR) |
| sw_dens | 原位密度 | sw_dens(S,T,P) |
| sw_pden | 位势密度(势密度) | sw_pden(S,T,P,PR) |
| sw_dens0 | 海表密度(0 dbar) | sw_dens0(S,T) |
| sw_vel | 声速 | sw_vel(S,T,P) |
| sw_salt | 电导率换算盐度 | sw_salt(C,T,P) |
| sw_cp | 比热容 | sw_cp(S,T,P) |
| sw_alpha | 热膨胀系数 | sw_alpha(S,T,P) |
| sw_beta | 盐收缩系数 | sw_beta(S,T,P) |
| sw_pres | 深度换算压力 | sw_pres(DEPTH,LAT) |
| sw_dpth | 压力换算深度 | sw_dpth(P,LAT) |
| sw_f | 科里奥利参数 | sw_f(LAT) |
其中sw_salt是处理原始CTD数据的第一步,因为它把电导率C(单位S/m)结合温度和压力换算成盐度。很多刚入门的同学不知道这一步,直接用CTD导出数据里的"盐度"字段,如果那个字段已经是仪器软件换算好的还好说,如果只是原始电导率,后面所有计算都会出问题。
2. 安装与路径配置:让老工具包在新版MATLAB里跑起来
2.1 三步完成安装:下载、解压、addpath
seawater工具包的安装是我见过最省事的一个,不需要install函数,不需要编译,不需要配置环境变量。下载方式主要有两个途径:CSIRO官网的seawater主页,以及GitHub上的seawater仓库。下载后得到一个zip压缩包,解压出来是一个名为seawater的文件夹。
打开MATLAB,在"主页"选项卡里找到"设置路径"(或"Set Path"),点击"添加并包含子文件夹",选中刚才解压出来的seawater文件夹,然后点击"保存"。如果想用命令行完成,直接在命令窗口输入:
addpath(genpath('你的路径/seawater')); savepath;注意genpath会把子文件夹也加进去,这个工具包内部有test等子目录,加进去不影响使用。验证是否安装成功很简单:
which sw_ptmp如果返回了seawater文件夹下的路径,说明安装成功。这个验证步骤我建议每次装完都做一次,因为我自己遇到过解压路径带了空格导致路径识别失败的情况,which命令能最快暴露问题。
2.2 版本差异与温度尺度:IPTS-68和ITS-90的坑
seawater 3.x版本与2.x版本最核心的差异在于温度基准。EOS-80本身是基于IPTS-68温标建立的,而现代CTD和温度传感器都使用ITS-90温标。如果输入的温度是ITS-90但工具包按IPTS-68计算,产生的误差虽然极小(海洋温度范围内约0.005°C),但在计算密度、位温时会被放大,尤其是在深层冷水区域。
判断你手上的seawater是哪个版本,可以在MATLAB里输入:
help sw_ptmp如果帮助文本明确写了"Temperature in ITS-90",就是3.x系列。如果是3.2之后的版本,还支持在函数调用时通过第四个参数设置温标。我个人建议直接下载3.3版本,因为现在绝大部分CTD数据已经是ITS-90了,没必要为了兼容旧数据自找麻烦。
如果你确实需要处理一批老文献或老航次数据,而这些数据记录的是IPTS-68温度,那就在输入seawater函数前先做一次转换:ITS-90温度约等于IPTS-68温度乘以0.99997再减去0.00004,这个近似在海洋温度范围内足够用了。但说实话,现在遇到这种情况的概率很低,更多时候是反过来的——新用户误以为seawater还在用旧温标,结果多此一举。
2.3 与gsw工具箱并存:函数前缀不同,互不干扰
很多人在装了seawater之后,又因为期刊要求而装了gsw工具箱(TEOS-10标准的官方实现)。这两个工具箱能不能共存?答案是肯定的。gsw的所有函数都以gsw_开头,与sw_开头的seawater函数完全不会重名,两者可以安全地同时出现在MATLAB路径里。
但我有一条经验要提醒:同一个脚本里不要混用两个标准下的函数。比如你用sw_ptmp算出了位温,然后把它传给gsw_rho去算密度,这在概念上是不对的,因为gsw是基于绝对盐度和保守温度的。跨标准混用会导致结果在近岸、河口等区域出现明显偏差。我的做法是:要么全部用seawater处理到底,要么全部用gsw处理到底,只在最后输出图件时做交叉验证对比。
3. 核心函数实战:从CTD剖面到位温、密度与T-S图
3.1 第一步:从深度计算压力,输入数据维度别搞错
拿到一次CTD下放的数据后,第一步通常是把深度换算成压力。这里要先弄清楚你的CTD数据里保存的是深度还是压力。大部分科研级CTD(如Seabird SBE 911plus)同时输出深度和压力两个通道,但有些近岸小型仪器只给深度。
如果只有深度,用sw_pres换算:
depth = (0:10:2000)'; % 深度,单位:米 lat = 20; % 站位纬度,单位:度 p = sw_pres(depth, lat); % 压力,单位:dbar输入维度的问题在这里就要注意:depth是201×1的列向量,lat是标量。在Matlab R2016b之后的版本里,标量会自动扩展,计算没问题。但如果你用的是旧版本,或者lat本身是一个与depth等长的向量(比如航次中纬度一直在变化),就必须保证两者形状一致,否则会报维度错误。
3.2 第二步:位温计算sw_ptmp,参考压力是关键参数
位温是水团分析里绕不开的量。它的物理含义很直观:取一团水,假设它被绝热地从深度P移动到参考压力PR,在移动过程中没有与外界发生热量交换,那么它在PR处应该具有的温度就是位温。绝热压缩和膨胀只会改变温度,所以深海的冷而高压的水团被移到海面时,会因为减压膨胀而进一步降温。
sw_ptmp的函数签名是:
theta = sw_ptmp(S, T, P, PR)其中S是盐度,T是原位温度(度,ITS-90),P是原位压力(dbar),PR是参考压力(dbar)。最常用的参考压力是0,表示把水团绝热移动到海面,得到的就是我们常说的potential temperature(海表位温)。如果你研究的是某个特定深度层的水团,也可以把PR设为那个层的压力,比如中层水分析常用1000 dbar。
我见过不少同学把PR参数漏掉,或者随意填一个数,这样算出来的位温完全不可比。位温必须说明参考压力才有意义,这是所有物理海洋学教材都会强调的。
3.3 第三步:密度、声速与T-S图
密度计算有两种:原位密度和位势密度。原位密度就是给定盐度、温度、压力状态下的实际密度,用sw_dens算:
rho = sw_dens(S, T, P); % 原位密度,单位kg/m^3 sigma_t = rho - 1000; % 密度异常sigma-t位势密度则是以某个参考压力为基准的密度,它消除了压缩效应,用于比较不同深度水团的"轻重"。常用的是参考0 dbar的位势密度异常,也就是sigma-theta:
sigma_theta = sw_pden(S, T, P, 0) - 1000;声速是声学海流计和声呐数据处理中的重要输入,seawater里一行代码就能算:
c = sw_vel(S, T, P); % 单位m/s密度和声速算出来后,最常见的数据可视化就是T-S图。T-S图的核心价值在于可以把一个剖面的水团特征压缩到一张二维图上:横轴是盐度,纵轴是温度,背景叠加等密度线。相同水团在T-S图上会聚成一团,不同水团的混合过程则表现为一条直线或曲线。
3.4 完整示例:一个2000米CTD剖面的全流程
我写了一个完整的处理脚本,从模拟CTD剖面到最终出图,你可以直接替换成自己的数据来跑:
% ===== 模拟一次CTD剖面 ===== % 深度序列:0~2000米,间隔10米 depth = (0:10:2000)'; lat = 20; % 北纬20度 % 由深度计算压力(单位:dbar) p = sw_pres(depth, lat); % 模拟温度剖面:海表约28.5°C,随深度递减 t = 28.5 - 25 * (1 - exp(-depth/300)); % 模拟盐度剖面:次表层盐度最大值,深层均匀 s = 34.6 + 0.5 * exp(-((depth - 200)/80).^2) + 0.2 * depth/2000; % ===== 计算位温 ===== % 参考压力0 dbar,即海表位温 theta = sw_ptmp(s, t, p, 0); % ===== 计算密度 ===== % 原位密度与sigma-t rho = sw_dens(s, t, p); sigma_t = rho - 1000; % 位势密度异常(参考0 dbar) sigma_theta = sw_pden(s, t, p, 0) - 1000; % ===== 计算声速 ===== c = sw_vel(s, t, p); % ===== 绘制T-S图,叠加等密度线 ===== figure('Color', 'w', 'Position', [100 100 750 500]); % 等密度线背景网格 s_grid = linspace(34, 36, 80); t_grid = linspace(0, 30, 80); [Sg, Tg] = meshgrid(s_grid, t_grid); sigma_grid = sw_dens0(Sg, Tg) - 1000; contour(Sg, Tg, sigma_grid, [21 22 23 24 25 26 27 28], 'k--', 'LineWidth', 0.6); hold on; % CTD剖面散点,颜色表示深度 scatter(s, t, 18, depth, 'filled'); colormap(jet); colorbar; xlabel('盐度 (PSS-78)'); ylabel('温度 (°C)'); axis([34 36 0 30]); title('T-S图:散点颜色代表深度'); grid on;运行之后你会看到一条从高温高盐端向低温方向延伸的曲线,散点颜色从红色逐渐变为蓝色,代表从表层到深层的过渡。这个图能直观看出次表层盐度最大值对应的温度范围,以及深层水团是否接近线性混合,这是做水团分析的基本功。
4. 最容易踩的坑:单位、参考压力与数据维度
4.1 分巴与米:什么时候能混用,什么时候不能
很多老海洋学家习惯说"800米处的密度",但代码里用的是800 dbar。这在大洋中上层问题不大,因为1 dbar约等于1.0197米水柱,误差不到2%。但在深海中这个误差会累积,如果用米代替分巴直接把2000米写成2000 dbar,密度计算在深层会有可见偏差,对精密研究(如地转流计算)是不可接受的。
正确的做法是从CTD压力通道直接读取压力,如果只有深度通道,务必用sw_pres换算。反过来也一样,如果你需要把等密度面深度画出来,先用sw_dpth把压力换回深度,别自己除以1.02。
4.2 盐度没有单位?PSS-78实用盐度的语义问题
PSS-78实用盐度在定义上是无量纲的,因为它本质上是电导率比值的函数,比值抵消了单位。但为了方便书写和交流,文献里普遍写成"PSU"(Practical Salinity Unit)或者直接写"Salinity"并标注PSS-78。我在图件坐标轴上一般写"盐度 (PSS-78)"或者"Salinity",不写PSU,因为严格来说PSS-78不是一个单位制。
这个细节在投稿时会被审稿人注意。如果你在论文里写"salinity in PSU",有些较真的审稿人会要求改成"PSS-78"或直接删除单位。我自己就因为这个被改过一次,后来养成了习惯:图上标注PSS-78,文中第一次出现时写清楚"practical salinity (PSS-78)",之后统一用S表示。
4.3 位温参考压力是哪个:说"位温"之前必须先说"参考多少"
位温这个名词如果不带参考压力,是没有明确指向的。默认情况下海表位温的参考压力是0 dbar,这也是大多数文献里"potential temperature"的含义。但如果你研究的是深层水团,有人会用1000 dbar或2000 dbar作为参考压力,这样算出来的位温数值会有差异,在对比文献时如果不注意参考压力,很容易得出错误结论。
同理,位势密度有各种命名:sigma-theta(参考0 dbar)、sigma-1(参考1000 dbar)、sigma-2(参考2000 dbar)、sigma-4(参考4000 dbar)。seawater的sw_pden通过PR参数区分:PR=0算出来的是sigma-theta,PR=1000是sigma-1,以此类推。引用别人数据时,先看方法部分写的是哪个sigma。
4.4 MATLAB旧版本与隐式扩展问题
seawater函数支持标量、向量、矩阵输入,但要求所有输入的维度可以广播。Matlab R2016b之前没有隐式扩展功能,如果你传入的盐度是201×1向量、温度是1×201矩阵,函数会报错。现在的新版本大多数情况能自动扩展,但为了代码健壮性,建议在调用前用size函数检查一下输入维度。
还有一个容易被忽略的点:sw_salt(电导率换算盐度)的输入电导率单位是S/m,也就是西门子每米。有些CTD厂商导出的电导率单位是mS/cm(毫西门子每厘米),两者相差10倍,如果不换算直接代入,算出来的盐度会完全离谱。我每次处理新航次数据时都会先随机抽几个值手工核对一下盐度范围,如果表层盐度算出35但实际应该在33左右,第一反应就查单位换算。
5. 要不要换TEOS-10:两个标准的比较与我的建议
5.1 EOS-80和TEOS-10到底差在哪
这里不展开复杂的数学推导,只说核心差异。TEOS-10是2010年国际海洋学组织推荐的新标准,它用Gibbs函数热力学势代替了EOS-80的经验多项式拟合。两者对盐度的定义不同:EOS-80用实用盐度SP,基于电导率;TEOS-10用绝对盐度SA,考虑了海水化学成分的微小变化。温度方面,TEOS-10引入了保守温度CT代替位温,理论上在能量守恒上更自洽。
差异有多大?对外海开阔海域的大多数工况,seawater和gsw工具箱计算出来的密度差在0.005到0.01 kg/m^3量级,位温差在0.01°C以内。这个误差对大多数物理海洋学分析来说可以忽略。但在近岸淡水注入区、冰川融水影响区、海底热液口附近,由于海水成分偏离标准海水,SA与SP的差异会明显变大,密度偏差可能达到0.1 kg/m^3以上,这时候必须用TEOS-10。
| 对比项 | seawater (EOS-80) | gsw (TEOS-10) |
|---|---|---|
| 盐度定义 | 实用盐度SP (PSS-78) | 绝对盐度SA |
| 温度变量 | 位温theta | 保守温度CT |
| 理论基础 | 经验多项式拟合 | Gibbs函数热力学势 |
| 计算成本 | 低 | 较高(依赖gsw工具箱) |
| 近岸/河口适用性 | 一般 | 更好 |
| 安装复杂度 | 轻量,纯m文件 | 依赖gsw工具箱 |
5.2 什么场景继续用seawater,什么场景该换gsw
我的经验是分三种情况。
第一种,日常快速处理和教学演示,用seawater。它轻量、直观、函数命名好记,学生们学起来没有负担。画T-S图、算个密度剖面,sw_ptmp和sw_dens两行代码搞定,效率很高。
第二种,正式发表论文且研究区域在开阔大洋,seawater的结果完全够用。很多高档次海洋学期刊并没有强制要求TEOS-10,只要你能在方法部分清楚说明自己用的状态方程版本,审稿人一般不会纠结。但要注意图件和表格里标清楚计算方法,避免被质疑。
第三种,研究区域有明显淡水影响、需要高精度能量收支分析、或者期刊明确要求TEOS-10标准,这时候直接上gsw工具箱。gsw的官方文档写得非常详细,每个函数都有理论背景和引用文献,写方法部分很方便。
5.3 实测对比:同一组CTD数据在两个标准下的差异
我曾经拿一组西北太平洋开阔海域的CTD数据分别用seawater和gsw算过密度剖面。在1000米以浅,两者sigma-t差异基本在0.005以内;2500米以深,差异略微增大但也没超过0.01。这个量级对地转流计算的影响几乎可以忽略。
但在一次处理长江口附近数据时,近岸低盐站位表层密度的差异就明显了,能达到0.08到0.1。原因很简单:长江冲淡水改变了局地海水离子组成,实用盐度SP无法反映这部分成分变化,只有绝对盐度SA能通过经纬度和深度信息做修正。那次之后我养成了一个习惯:先看站位离岸距离和水深,再决定用哪个工具包。
我现在的做法是:日常快速分析和教学用seawater,正式发表研究用gsw。如果你也经常跟CTD数据打交道,建议两个都装上,反正不冲突。最后提醒一句:不管用哪个工具包,拿到数据第一件事先确认单位和温度尺度,别急着跑函数。这个教训是我用无数次返工换来的,希望你能少走一次弯路。
本文还有配套的精品资源,点击获取