简介:本资源是面向海洋科学、物理海洋学及水文工程领域科研人员与高校师生的专业计算工具包,提供基于国际标准TEOS-10(2010海水热力学方程)的GSW(Generalized Seawater)Matlab实现,用于高精度计算海水密度、比热容、声速、位势密度、地转流函数等关键物理属性。压缩包共345个文件,主体为337个Matlab函数(.m),涵盖核心算法(如gsw_gibbs、gsw_geo_strf_dyn_height)、数值稳定性处理(如gsw_stabilise_SA_CT)、完整性校验脚本(如gsw_check_functions)及模块化说明(Contents.m),辅以README和LICENSE等必要文档,总大小仅4.1MB,轻量易集成。已有362人学习下载,适用于海洋环流建模、CTD数据处理、气候模型参数化、潜水器浮力设计等实际科研与工程场景,开箱即用,无需额外依赖,显著提升海水物性计算的规范性与效率。 开头先讲个真实经历。我曾经在处理深海CTD数据时,用旧版SEAWATER工具箱算出来的位密断面图,和同事用TEOS-10标准算出来的结果在3500米以下总是有明显偏差。一开始以为是数据采集问题,后来才意识到,问题出在两个被海洋学界沿用了四十年的核心变量上——位温和实用盐度。于是我把整套数据处理流程迁到了TEOS-10,用起了官方发布的GSW工具箱。今天这篇围绕TEOS-10-GSW-Matlab.zip展开,聊聊这个压缩包怎么用、背后的变量体系为什么变了、以及我在真实航次数据处理中踩过的那些坑。
1. 为什么还在用EOS-80的人,该升级到TEOS-10了
先说结论:TEOS-10是用Gibbs函数描述海水热力学状态的新国际标准,它的Matlab实现就是GSW Oceanographic Toolbox,官方打包成zip形式发布,文件名类似gsw_matlab_v3_06.zip,也就是标题里那个TEOS-10-GSW-Matlab.zip。老一代的EOS-80连同SEAWATER工具箱(就是那一堆sw_开头的函数)在计算精度和物理一致性上,已经跟不上海洋学研究的需求了。
1.1 从“位温”到“保守温度”:一个看似微小却影响深远的变量替换
老工具箱里最常用的一个变量是位温(potential temperature),它的定义是海水从当前深度绝热位移到海面后所具有的温度,用于消除压缩效应带来的温度变化。这个变量在EOS-80体系里被当作保守量使用,但严格来说位温在等熵混合过程中并不是真正守恒的,尤其是在深层水形成、跨等密度面混合显著的区域,位温的“非保守性”会引入误差。
TEOS-10用保守温度(Conservative Temperature, CT)替代位温作为主要温度变量。保守温度的物理基础是焓,它被设计成在无相变、无热通量的等压混合过程中严格守恒。这意味着你做水团混合分析、端元分析时,用CT比用位温稳得多。
我举个实际例子。在北大西洋深层水与南极底层水的混合带上,两套水体温度差异不大,但位温和保守温度之间的偏差会在数千分之一的量级上影响混合比例估算。如果你做的是世纪尺度的水团变化研究,这个偏差足够让趋势信号失真。
1.2 绝对盐度与实用盐度:差的那零点几克哪里来的
实用盐度(Practical Salinity, SP)是PSS-78标准下的电导率盐度,它本质上是电导率比的一个换算结果,无量纲,数值上接近ppt但绝不是质量分数。TEOS-10改用绝对盐度(Absolute Salinity, SA),单位是g/kg,代表海水中溶解物质的总质量分数。
SA和SP之间存在一个校正量δSA = SA - SP。在大洋大部分海域,δSA的量级约为千分之几到百分之几g/kg,但在一些陆架海、半封闭海盆以及深层水形成区,δSA可以更大。这零点几克每千克的差异,对盐度绝对值的表述影响不大,但如果你在追踪0.01 g/kg量级的盐度异常,或者在做地转流计算,δSA的空间不均匀分布就会产生不能忽略的系统误差。
1.3 升级带来的实际收益
GSW工具箱不是简单地把密度公式改精确了一点,而是把整个热力学一致性问题解决了。它可以给出焓、熵、位焓、声速、热膨胀系数、盐收缩系数等一系列热力学变量,而且这些量之间严格满足热力学关系。举个例子,EOS-80里的热膨胀系数和密度是分开参数化的,两者在局部可能不一致;而TEOS-10里所有性质都从同一个Gibbs函数导出,物理上自洽。
对你手头的工作而言,最直接的收益包括:
- 密度、位密计算更精确,尤其是深层和大压力范围;
- 保守温度和绝对盐度做水团混合分析时不会引入假性的非保守信号;
- 能量、热通量、地转流计算有了统一规范的变量体系;
- 官方提供了波罗的海、北冰洋等区域专用盐度修正公式,可以覆盖更多特殊海区。
2. TEOS-10-GSW-Matlab.zip 的正确打开方式
这个压缩包本身并不复杂,但安装过程中的一些细节会直接影响后续使用。我见过不少人解压完就直接开跑,结果报错或者算出来的SA明显异常,最后发现是路径问题或者没有正确挂载子目录。
2.1 解压与目录规划
从TEOS-10官网下载zip后,先别急着双击解压到默认目录。建议把解压路径设成不含中文、不含空格的纯英文路径,比如D:\Ocean_Tools\gsw_matlab_v3_06或~/matlab/gsw。Matlab对路径里的中文和空格支持其实比前些年好多了,但海量数据脚本一跑起来,编码问题还是会冷不丁冒出来,没必要给自己埋雷。
Linux下用unzip命令解压很直接:
unzip gsw_matlab_v3_06.zip -d ~/matlab/解压后你会看到根目录下有一堆gsw_开头的.m函数文件,还可能有一个library子目录,里面是底层的核心函数,以及html文档目录和测试脚本test_gsw_matlab_v3_06.m。不同版本目录结构略有差异,但大框架一致。
这里有个容易被忽略的点:如果你下载的压缩包名字是TEOS-10-GSW-Matlab.zip,解压出来的文件夹名可能带版本号,也可能不带。为了后续升级方便,我习惯把文件夹改成类似gsw_v3_06这样能分辨版本的名称,同时保留原始函数文件不动。
2.2 在Matlab里挂载工具箱的两种方式
官方文档建议用addpath(genpath(...))把整个目录树加入Matlab路径。genpath会递归包含所有子目录,避免漏掉library等深层目录里的函数。
addpath(genpath('D:\Ocean_Tools\gsw_matlab_v3_06')); savepath;第一条命令将路径加入当前会话,savepath写入pathdef.m,让之后每次启动Matlab都自动加载。如果不希望写入全局路径,也可以不加savepath,每次新会话手动执行第一行。
另一种方式是用Matlab主界面上的Set Path按钮,GUI操作,点两下就能完成。但GUI方式在多版本切换时不方便,我建议直接用脚本。
2.3 安装验证
路径挂载完成后,先运行两个快速验证:
which gsw_SA_from_SP gsw_SA_from_SP(35, 100, -150, -10)which用于确认当前解析到的是你刚挂载的那个函数,而不是其他目录下的同名文件。第二个表达式调用绝对盐度计算函数,输入实用盐度35、压力100 dbar、经度-150度、纬度-10度,返回一个接近35的数值。只要结果是有限的、量级正确,基本可以认为安装成功。
更完整的验证是运行根目录下的测试脚本test_gsw_matlab_v3_06.m。这个脚本会执行一系列官方自检,运行时间不长,输出里会出现大量测试点的OK信息。如果某个点报错,优先检查Matlab版本和路径冲突。
提示:如果以前装过旧版本的GSW工具包,新版本解压后不要直接扔着不管。用
which gsw_rho看看当前调用的是哪个路径下的函数,确保没有新旧版本同时在路径列表里。
3. 从sw_到gsw_:变量体系和函数调用的底层差异
很多人迁移GSW时,第一反应是把sw_前缀改成gsw_,然后发现一堆报错。这很正常,因为两套工具箱不只是函数名不同,更重要的是输入输出变量体系变了。不理解背后的逻辑,只做文本替换,早晚会在数据上翻车。
3.1 核心输入输出变化:从SP到SA、从t到CT
最核心的差异是变量换血。老工具箱里几乎所有函数都吃实用盐度SP和位温(或者原位温t),而GSW工具箱的主流函数吃绝对盐度SA和保守温度CT。
我列一个常用函数对照表:
| 旧SEAWATER函数 | GSW函数 | 关键说明 |
|---|---|---|
| sw_dens(SP,t,p) | gsw_rho(SA,CT,p) | 输入必须为绝对盐度和保守温度 |
| sw_dens(SP,t,p) | gsw_rho_t_exact(SA,t,p) | 如果只有原位温,用这个便捷函数 |
| sw_ptmp(SP,t,p,pr) | gsw_pt_from_t(SA,t,p,pr) | 位温计算,在GSW中定义为从原位温算到参考压力 |
| sw_svel(SP,t,p) | gsw_sound_speed(SA,CT,p) | 声速计算 |
| sw_alpha(SP,t,p) | gsw_alpha(SA,CT,p) | 热膨胀系数 |
| sw_beta(SP,t,p) | gsw_beta(SA,CT,p) | 盐收缩系数 |
| sw_pden(...) | gsw_sigma0(SA,CT) | 位密,GSW直接用sigma系列函数 |
特别强调一点:gsw_rho的第二个参数是保守温度CT,不是原位温t,更不是位温pt。如果你拿旧数据里的原位温直接丢进去,算出来的密度会和sw_dens有偏差,而且这种偏差不是公式精度导致的,是输入变量类型错了。
正确的做法有二:要么先算CT再调gsw_rho:
SA = gsw_SA_from_SP(SP, p, lon, lat); CT = gsw_CT_from_t(SA, t, p); rho = gsw_rho(SA, CT, p);要么直接用gsw_rho_t_exact(SA, t, p)一步到位。我个人更推荐后者,少一次中间转换,调用语义也更清晰。
3.2 为什么经纬度坐标变得不可省略
老工具箱里sw_dens(SP,t,p)给三个参数就完事了。GSW里计算绝对盐度的标准函数gsw_SA_from_SP(SP, p, lon, lat)需要经度和纬度,这不是故意刁难,而是因为δSA的全球修正量依赖于地理位置。
TEOS-10文献里给出了全球海洋的δSA参考分布,这个分布是基于大量实测数据和海洋环流模式构建的。每个经纬度格点上都有一组对应的修正系数。如果没有经纬度,就无法查表插值,计算自然不可靠。
我在实际项目中见过有人图省事把lon和lat全填0,结果在赤道附近误差还能接受,一到高纬度区域SA偏差就会明显变大。如果你手里只有站位号没有经纬度,先根据航次记录把站位对应的经纬度补全,再进入计算流程,这是最稳妥的。
如果确实没有经纬度,且只做局地近似,可以关注GSW工具箱里针对特定海区提供的专用函数,比如波罗的海就有gsw_SA_from_SP_baltic。这些区域化公式牺牲了全球通用性,换来了局部精度。
3.3 压力与深度的互算
另一个高频问题是压力与深度互算。很多旧脚本习惯用深度/10近似压力,因为1 dbar约等于1米水深。但在精密计算中,压力还受纬度和实际海水压缩性影响。
GSW提供两个函数:
p = gsw_p_from_z(z, lat); z = gsw_z_from_p(p, lat);只需要深度(或压力)和纬度。为什么需要纬度?因为地球重力加速度随纬度变化,同样深度下的实际压力不同。在5000米深度,不同纬度之间的压力差异有几十dbar,对深层密度和位密的计算影响不可忽略。
4. 一次完整的实操:从CTD原始数据到断面密度图
理论说再多,不如跑一遍流程。下面用一个很典型的CTD数据处理场景串一遍GSW的典型用法。
4.1 数据准备与读取
假设你有一个航次断面的CTD数据文件ctd_cruise.csv,列包括:
station:站位号depth:深度(米)pressure:压力(dbar)temp:原位温(ITS-90,摄氏度)SP:实用盐度(PSS-78)lon、lat:经度、纬度
用Matlab读取:
data = readtable('ctd_cruise.csv'); p = data.pressure; t = data.temp; SP = data.SP; lon = data.lon; lat = data.lat;有些CTD数据文件里压力列缺失,只有深度,那就先转压力:
p = gsw_p_from_z(data.depth, lat);4.2 核心变量转换段
得到四个基本量之后,进入GSW的标准转换流程:
SA = gsw_SA_from_SP(SP, p, lon, lat); CT = gsw_CT_from_t(SA, t, p); sigma0 = gsw_sigma0(SA, CT);这里每一步都值得解释一下:
gsw_SA_from_SP完成实用盐度到绝对盐度的修正,需要经纬度;gsw_CT_from_t用SA、原位温和压力计算保守温度;gsw_sigma0算出相对参考压力0 dbar的位密,返回的是sigma值,即密度减1000 kg/m³。如果断面水深大,需要看深层等密度面,可以根据深度范围选gsw_sigma1、gsw_sigma2等。
4.3 绘制断面位密图
断面图目的通常是看水团结构。经典的画法是横轴为站位距离或经度、纵轴为压力(或者深度),填色加等值线:
figure; contourf(dist, p, reshape(sigma0, nstation, nlevel), 30, 'LineStyle', 'none'); set(gca, 'YDir', 'reverse'); xlabel('Distance (km)'); ylabel('Pressure (dbar)'); c = colorbar; c.Label.String = 'sigma0 (kg/m^3)';为什么用位密而不是原位密度?因为原位密度沿压力面变化太大,等值线基本是顺着深度走的,反映不出水团内部的分层信息。位密消除了压缩效应,能更清晰地展示水团的密跃层和等密度面结构。
4.4 一个容易出错的细节:盐度参数单位的理解
很多人在这一步栽跟头。CTD原始数据里的SP虽然数值上接近ppt,但它其实是无量纲的。GSW中SA的单位是g/kg。有些数据文件里列名写着PSU,其实和SP是同一个东西,可以直接作为SP输入。但如果你拿到的数据给的是绝对电导率(单位mS/cm),那就必须先转成电导率比,再计算SP。
标准代码如下:
% C是电导率值,单位mS/cm R = C / 42.914; % 42.914是C(35,15,0)的标准值 SP = gsw_SP_from_C(R, t, p);这个坑在近岸观测数据里特别常见,仪器输出的电导率往往不是盐度,需要自己转换。
4.5 与旧方法对比
为了直观感受TEOS-10带来的差异,可以在同一份数据上同时用sw_dens和gsw_sigma0计算位密,做差看看分布:
% 假设你还有旧SEAWATER工具箱在路径里 sigma0_old = sw_pden(SP, t, p, 0) - 1000; % 旧工具箱位密 delta = sigma0 - sigma0_old;在大深度处,两套标准计算出的位密差异可以达到0.01 kg/m³的量级,这在某些注重精密诊断的研究里是不可忽略的。这个差值主要来源于密度公式的更新与绝对盐度的修正。
5. 我在使用中踩过的坑(含解决思路)
下面这些坑,基本都是我在处理不同航次数据时真实遇到过的,写出来帮你省点时间。
5.1 经纬度缺失导致SA计算异常
有一年的航次数据只存了站位编号,站位对应的经纬度放在另一张表格里。我同事在跑GSW时没有关联经纬度,直接把lon和lat写成0,结果整个断面的SA都偏大,深层等密度面深度全部偏移。排查过程花了两个小时,最后定位到是经纬度缺失。
解决方案很简单:先用站位表把经纬度合并进数据表,再跑后续流程。如果有少量站位缺失经纬度,宁可用相邻站位的经纬度插值,也别直接填0。
5.2 把原位温直接输入gsw_rho导致密度算错
这个坑出现在一次跨海域的密度剖面对比中。我用gsw_rho_t_exact算了一个剖面,又用gsw_rho算另一个剖面,发现两边的密度差异不符合物理预期。后来检查发现,用gsw_rho的那个脚本里,第二个参数传的还是原位温,而不是保守温度。
这种问题在批量处理时非常隐蔽,因为数据不会报错,只是结果偏了。我的做法是写一个本地封装函数,把SA和CT的转换放在最前面:
function [rho] = my_rho_from_SP_t(SP, t, p, lon, lat) SA = gsw_SA_from_SP(SP, p, lon, lat); CT = gsw_CT_from_t(SA, t, p); rho = gsw_rho(SA, CT, p); end之后所有下游计算都用这个封装,避免在十几个脚本里反复出低级错误。
5.3 波罗的海、黑海、河口区等非开阔大洋区域的SA修正
GSW全球修正量参数化在开阔大洋的适用性较好,但在波罗的海、黑海、北冰洋陆架以及河口冲淡水区,δSA的参考修正量可能不够准确。对于波罗的海,官方提供了gsw_SA_from_SP_baltic专用函数,直接用它就行。
在长江口这类冲淡水和海水混合的高动态区域,我发现gsw_SA_from_SP算出的SA在低盐端可能变成负值,这显然没有物理意义。遇到这种情况,可以对比历史本海域的盐度-密度关系,判断是否需要引入局地修正公式,或者对SA做下限截断。
5.4 旧脚本批量替换时的隐藏问题
批量替换不是简单的sw_换成gsw_。参数顺序变了、单位变了、前置步骤变了,这些都必须一起处理。我整理过一个迁移检查清单:
- 输入变量是否已从SP转成SA?从t转成CT?
- 经纬度是否补齐?用的全球修正还是区域修正?
- 压力是否用
gsw_p_from_z重新计算? - 输出量纲是否匹配?sigma系列返回的是密度减1000后的值。
- 新旧工具包是否同时存在于路径中?用
which确认。
这个清单看起来简单,但只要其中一条没做到,最终结果就可能差之毫厘谬以千里。
5.5 版本差异:v3_05/v3_06中函数命名和返回参数的变化
GSW工具箱持续更新,不同版本中函数的命名和返回参数可能有调整。老版本里某些函数在v3_06之后增加了新版本,比如有些函数增加了异常值处理选项。我遇到过的情况是,用了某篇博文里的老函数名,在新版本中已经废弃,报错后去官网查release notes才明白。
建议下载最新版zip,安装后尽快用evalc('help gsw')之类的命令生成一份函数列表,存成文档备查。同时确认你没有把旧版和新版同时放进path,否则同名函数可能冲突。
6. 后续还能怎么用:GSW工具箱的高级功能与扩展
基础密度计算只是GSW的冰山一角。当数据处理链路稳定之后,你会发现自己还能顺手做很多过去要写一堆近似代码才能完成的计算。
6.1 动力高度与地转流
地转流计算在传统方法里常基于密度场积分,公式繁琐。GSW提供了gsw_geo_strf_dyn_height系列函数,可以直接从SA、CT和压力积分出动力高度。日常做法是确定一个参考面(比如1000 dbar或1500 dbar),然后计算两个站位之间动力高度的差值,进而得到地转流速。
dyn_height_ref = gsw_geo_strf_dyn_height(SA, CT, p, 1000);如果你之前用EOS-80手动处理,换成TEOS-10后最明显的感受是动力高度的物理定义更清晰,而且不用自己维护一堆积分系数。
6.2 冰-海相互作用
极地研究里经常要估计海冰融化对海水的稀释和冷却效应。GSW里有一组冰相关函数,比如gsw_SA_CT_ice_melting、gsw_T_ice_melting_poly等,可以算出给定质量的海冰融化后海水SA和CT的变化量。过去自己做这类计算往往要查海冰潜热、比热等参数,现在工具箱直接给。
6.3 与其他语言的互通
GSW官方不仅提供Matlab版本,还有Python版、Fortran版和C版。Python里用pip install gsw就能装,接口设计和Matlab版几乎一一对应。这意味着你用Matlab调通的流程,可以很轻松地迁移到Python数据处理管线里,方便团队里不同技术栈的人协作。
6.4 热力学诊断量:声速、位熵、焓
声速剖面是声学探测和海底地形反演的基础参数。GSW的gsw_sound_speed(SA, CT, p)可以替代旧的sw_svel,而且基于Gibbs函数得到的声速与其他热力学量严格一致。
位熵和焓也是水团分析里的强大工具。gsw_entropy_from_CT和gsw_pot_enthalpy_from_pt这些函数,在做等熵面分析和能量通量诊断时非常好用。简单说,只要你能想到的热力学量,在GSW里基本都有对应的函数,而且相互之间物理自洽。
最后分享一点个人体会。TEOS-10不是那种装上就能立刻让所有结果天翻地覆的标准,它的价值体现在你开始深挖物理机制、做精细混合分析、跨数据集对比的时候。我迁移到GSW之后,最明显的改变不是某个数值变准了,而是整个数据处理流程中不再需要担心“这个近似是否引入虚假信号”。如果你还在用旧工具箱,建议花半天时间把存量脚本迁过来,后面会省下很多冤枉时间。
本文还有配套的精品资源,点击获取