简介这份资源是面向海洋科学计算与工程应用场景的 TEOS-10-GSW Matlab 工具包适合从事物理海洋学、气候建模、浮标数据处理及潜水器设计等方向的科研人员与工程师使用。它把国际通用的 TEOS-10 海水热力学方程与 GSW 广义海水算法封装为 Matlab 接口让使用者无需深究底层数学即可计算温度、盐度、压力、密度、比热容等海水属性。压缩包共 345 个文件以 337 个 .m 脚本与函数为主体另含少量 .p、.mat、.xml、readme、md 与 txt 文件整体约 4.1MB结构紧凑、便于按模块检索。内容涵盖 GSW 函数接口、校验脚本、示例调用、测试数据与说明文档可帮助读者快速在 Matlab 中完成海洋水文物理计算与结果验证。目前已有 362 人学习下载适合需要系统掌握 TEOS-10 标准计算流程的中高级用户参考。1. 海水物性计算为什么绕不开 TEOS-10 与 GSW从一份 Matlab 工具箱说起做海洋观测、温盐深数据处理或者数值模式后处理的人迟早会撞上一件事手里有一堆温度、盐度、压力剖面想算密度、声速、热膨胀系数结果发现不同文献给的公式对不上同一组数据算出来的密度能差出零点零几。这不是你算错了而是海洋学在 2010 年前后换过一次「度量衡」。老的 EOS-80 用实用盐度 PSU 和 IPTS-68 温度新的 TEOS-10 改用绝对盐度 SA 和保守温度 CT两者在深层水和极地海域的偏差会明显放大。GSW 就是 TEOS-10 的官方参考实现全称 Gibbs SeaWater Oceanographic ToolboxMatlab 版本通常打包成 TEOS-10-GSW-Matlab.zip 这样的形式分发。它解决的核心问题只有一个用一套热力学自洽的公式把海水的基本状态量算准并且让全球做海洋的人算出来是同一个数。适合谁凡是碰 CTD 剖面、Argo 浮标数据、海洋模式输出、声场建模的人都该把这套东西装进自己的工具链而不是继续抄十几年前的旧函数。2. TEOS-10 与 GSW 的底层逻辑为什么不能再用 EOS-80 凑合2.1 从实用盐度到绝对盐度差的不只是一个系数EOS-80 时代的盐度是「实用盐度」靠电导率比值定义本质是一个经验标定值。它的问题在于海水不是单纯的 NaCl 溶液不同海域的离子比例不同电导率换算出来的盐度并不等于真实的质量分数。TEOS-10 引入绝对盐度 SA单位是 g/kg直接对应溶液的质量分数并且通过一个叫 SRReference Salinity的中间量把电导率盐度、营养盐修正、质量守恒串起来。你在 GSW 里会看到SA_from_SP这个函数它做的就是这件事输入实用盐度 SP、经度、纬度、压力输出绝对盐度 SA。经度和纬度不是摆设因为 SR 与 SP 的换算系数随地理位置变化这是 EOS-80 完全没有考虑的因素。保守温度 CT 同理。老的温度是「现场温度」受压力影响同一水团在不同深度测出来的温度不一样做水团分析时非常别扭。CT 是把现场温度绝热修正到参考压力后的值相当于给温度也做了「位势化」。GSW 里CT_from_t负责这个转换。你如果直接拿现场温度去算密度在 1000 米以下误差会累积做深层水团识别时可能把两个水团混成一个。2.2 GSW 工具箱的模块划分与调用习惯GSW 的 Matlab 版本不是一个大函数而是一组按功能分类的 m 文件加一个 MEX 编译的底层库。常见目录结构里你会看到gsw_*前缀的函数比如gsw_rho、gsw_CT_from_t、gsw_sound_speed。这些函数分几类转换类SP↔SA、t↔CT、物性类密度、声速、热容、膨胀系数、稳定性类位势密度、中性密度、浮力频率、以及少量专门用于冰-海-气界面的函数。调用习惯上有一个坑GSW 的函数大多要求输入是绝对盐度 SA 和保守温度 CT而不是你从 CTD 原始文件里直接读出来的 SP 和 t。所以标准流程永远是「先转换再计算」。我一般会写一个包装函数把原始剖面一次性转成 SA、CT后面所有物性计算都基于这两个量避免在几十个地方重复转换、重复出错。% 假设从 CTD 文件读入的列向量SP(实用盐度), t(现场温度), p(压力 dbar) % lon, lat 为剖面所在位置的经纬度度 % 第一步实用盐度转绝对盐度 SA gsw_SA_from_SP(SP, p, lon, lat); % 第二步现场温度转保守温度 CT gsw_CT_from_t(SA, t, p); % 第三步基于 SA、CT、p 计算密度 rho gsw_rho(SA, CT, p); % 第四步计算声速用于声场建模 c gsw_sound_speed(SA, CT, p);这段代码的逻辑很直白GSW 的所有物性函数签名基本是gsw_xxx(SA, CT, p)顺序固定。参数说明上SP 和 t 是原始观测值p 是压力lon/lat 是剖面位置缺一不可。注意gsw_SA_from_SP的 lon/lat 可以是标量也可以是数组如果整条剖面在同一位置传标量即可如果是走航观测每个点位置不同就要传同长度的向量。密度输出单位是 kg/m³声速是 m/s这些在函数头注释里都有但很多人不看直接拿去用单位错了都不知道。2.3 安装与路径配置别让 MEX 编译卡住你TEOS-10-GSW-Matlab.zip 解压后通常包含gsw_matlab和gsw_matlab_v3_06_xx之类的目录里面既有 m 文件也有 C 源码。Matlab 版本要能调用底层 C 库需要先编译 MEX。常见做法是在 Matlab 命令行里进入工具箱目录运行gsw_install或者手动mex编译。如果你机器上没有配置好 C 编译器这一步会报错很多人就卡在这里。我一般会先确认 Matlab 的mex -setup能识别到编译器。Windows 上常见的是 MinGW-w64Linux 上就是 gcc。如果编译失败先看报错里有没有「找不到头文件」或者「架构不匹配」前者通常是编译器路径没配好后者是 32/64 位混用。编译成功后用which gsw_rho确认路径已经加到 Matlab 搜索路径里。如果没加手动addpath(genpath(你的工具箱路径))并savepath。提示不同 Matlab 版本对 MEX 的支持有差异R2018b 之后对 C99 的支持更完整。如果你用的是很老的版本建议先升级 Matlab而不是硬改源码。3. 用 GSW 算密度、声速与稳定性一份可复现的剖面处理脚本3.1 读数据与预处理的三个检查点拿到一份 CTD 剖面别急着算。我一般先做三件事检查压力是否单调递增、检查有没有重复深度、检查盐度和温度有没有明显的尖峰。压力不单调会导致后续插值出错重复深度会让某些函数返回 NaN。尖峰可能是传感器故障也可能是真实的小尺度信号需要结合原始记录判断。% 读入数据假设是三列文本p, t, SP data load(ctd_profile.txt); p data(:,1); t data(:,2); SP data(:,3); % 检查压力单调性 if any(diff(p) 0) warning(压力序列非单调可能存在重复或倒转); % 常见处理排序并去重 [p, idx] sort(p); t t(idx); SP SP(idx); [p, ia] unique(p, stable); t t(ia); SP SP(ia); end % 检查 NaN valid ~isnan(p) ~isnan(t) ~isnan(SP); p p(valid); t t(valid); SP SP(valid);这段代码的关键在unique的stable选项它保留第一次出现的顺序避免打乱剖面结构。参数上diff(p) 0判断非递增sort后重新索引保证后续计算不会因为顺序问题出错。很多人跳过这一步结果gsw_rho返回一堆 NaN回头查半天。3.2 计算位势密度与浮力频率位势密度是海洋学里判断水团稳定性的核心量。GSW 提供gsw_rho算现场密度gsw_CT_from_t和gsw_SA_from_SP转换后再用gsw_rho算参考压力下的密度就是位势密度。浮力频率 N² 则用gsw_Nsquared它直接接受 SA、CT、p。% 转换 SA gsw_SA_from_SP(SP, p, lon, lat); CT gsw_CT_from_t(SA, t, p); % 现场密度 rho gsw_rho(SA, CT, p); % 位势密度参考压力 0 dbar rho_pot gsw_rho(SA, CT, zeros(size(p))); % 浮力频率平方 [N2, p_mid] gsw_Nsquared(SA, CT, p, lat);gsw_Nsquared返回的 N2 长度比输入少 1因为它算的是层间值对应压力是p_mid。这个细节在画图时很重要否则维度对不上。参数lat用于科氏参数计算不能省略。如果你做的是赤道区域lat 接近 0N2 的计算会涉及 f 的平方数值上仍然稳定但物理意义要小心解释。3.3 声速计算与声场建模的衔接声速是 GSW 里少数几个对渔业、水声、海洋哺乳动物研究都很实用的输出。gsw_sound_speed基于 TEOS-10 的热力学框架比老的经验公式比如 Mackenzie在深层和低温环境下更准。% 声速 c gsw_sound_speed(SA, CT, p); % 画剖面 figure; subplot(1,3,1); plot(rho, -p); xlabel(密度 kg/m^3); ylabel(压力 dbar); subplot(1,3,2); plot(c, -p); xlabel(声速 m/s); ylabel(压力 dbar); subplot(1,3,3); plot(N2, -p_mid); xlabel(N^2 s^{-2}); ylabel(压力 dbar);画图时把压力取负是为了让深度向下增加符合海洋学惯例。这一步没有技术含量但图看起来专业很多。声速剖面的拐点往往对应声道轴做声场建模时这个轴的位置直接决定声波能不能远距离传播。注意GSW 的声速是绝热声速如果你做的是低频、长距离传播还需要考虑频率修正那是另一个模型的事别把 GSW 的输出直接当成全频段声速。4. 避坑与排查GSW 在 Matlab 里最容易翻车的五个地方4.1 现象所有输出都是 NaN。原因输入里混了 NaN 或 Inf。解决在转换前统一清洗用valid索引过滤别指望 GSW 内部帮你处理。4.2 现象密度值比预期大很多。原因把实用盐度 SP 直接传给了gsw_rho而gsw_rho要的是绝对盐度 SA。解决先调gsw_SA_from_SP确认 SA 和 SP 的差值在合理范围一般 0.10.3 g/kg 量级。4.3 现象MEX 编译报错提示找不到gsw_teos10.h。原因编译时工作目录不对或者头文件路径没包含。解决在工具箱根目录下运行安装脚本不要在其他目录里手动mex。4.4 现象gsw_Nsquared返回的 N2 全是负值。原因剖面本身存在密度倒转可能是真实的对流不稳定也可能是数据错误。解决先检查温度和盐度剖面有没有明显的倒置如果是真实信号负 N2 就是物理结果不用改如果是传感器问题回到原始数据排查。4.5 现象同一组数据别人算的密度和我差 0.01。原因对方可能用了不同的参考压力算位势密度或者用了 EOS-80。解决统一用 TEOS-10 的 SA、CT并明确参考压力。做对比时把对方的公式和参数要过来逐项核对。5. 进阶用法把 GSW 嵌进批量处理与模式后处理流程5.1 批量处理 Argo 剖面的向量化写法Argo 数据动辄几千条剖面逐条循环调用 GSW 会慢得让人想砸键盘。GSW 的函数本身支持向量输入只要你把同一剖面的 SA、CT、p 组织成列向量一次调用就能算完整条剖面。更进一步如果你把多条剖面拼成矩阵每列一条剖面GSW 也能处理但要注意 lon/lat 的维度匹配。% 假设 SA_all, CT_all, p_all 是 n×m 矩阵每列一条剖面 % lon_vec, lat_vec 是 1×m 向量 rho_all gsw_rho(SA_all, CT_all, p_all); % 如果 lon/lat 每条剖面不同转换时要传对应向量 for j 1:m SA_all(:,j) gsw_SA_from_SP(SP_all(:,j), p_all(:,j), lon_vec(j), lat_vec(j)); end向量化的关键是提前把内存分配好别在循环里动态增长数组。我一般会预估剖面数量用zeros预分配再填充。这样处理一万条剖面比逐条循环快一个数量级。5.2 与海洋模式输出的对接海洋模式比如 ROMS、MITgcm输出的盐度通常是实用盐度或模式自定义盐度温度是位温或现场温度。对接 GSW 时第一步永远是确认模式输出的变量定义。如果模式用的是位温那gsw_CT_from_t的输入就要相应调整不能直接套。常见做法是写一个适配层把模式变量名映射到 GSW 的输入再统一转换。5.3 验证方法用官方测试数据对拍GSW 工具箱里通常带一组测试数据包含已知输入和期望输出。安装完成后运行测试脚本看所有断言是否通过。如果通过说明编译和路径都没问题。如果失败先看是哪个函数再对照官方文档检查输入维度。我习惯在每次升级 Matlab 或换机器后跑一遍测试花两分钟省掉后面几小时的排查。检查项期望结果常见异常MEX 编译无报错生成 mex 文件编译器未配置测试脚本全部通过路径未添加SA 与 SP 差值0.10.3 g/kg经纬度传错N2 符号稳定层结为正数据倒转这张表是我自己排查时用的清单每次换环境先过一遍。最后说个习惯我从不把 GSW 的输出直接写进论文或报告一定先画剖面图肉眼确认没有异常跳变。图对了数才敢用。希望帮到你。本文还有配套的精品资源点击获取