海洋要素计算：基于EOS-80的MATLAB工具箱全流程解析

简介这份资源是一套面向海洋科学研究者的海洋要素计算工具箱围绕海水温度、盐度、密度、压力及流速流向等关键参数提供了一系列MATLAB计算函数适合开展海洋数据处理、环流分析和模型构建的科研人员使用。包内共36个文件以34个m函数为主涵盖密度、位温、声速、混合层、地转流等常用算法另含1个mat数据文件与1个readme说明整体仅52KB轻量实用。目前已有352人学习下载。借助这些函数使用者无需从底层编写全部计算逻辑即可快速完成海洋观测数据的基础性物理量换算与剖面分析同时还能结合自带数据样例验证算法流程为后续研究、教学或项目开发提供直接支持。1. 海洋要素计算先过状态方程这道坎在南海投放一次 XCTD 剖面拿回的是温度、电导率和压力三组原始序列但论文里要的却是位温、位密、浮力频率和地转流速。这些量没有一个是传感器直接读出来的全部要靠海洋要素计算里的状态方程在温盐之间做换算。sea_water 这套工具包把 EOS-80 标准下三十多个计算函数打包在一起从纯水密度到声速、从实用盐度到位密一条命令就能从原始观测落到物理量。对做物理海洋、走航观测和数值模式后处理的工程师来说这是绕不开的基础工具包。本文从 sw_data.mat 里的实测剖面出发把密度、盐度、浮力频率和地转流计算完整跑一遍最后用 sw_new 示范如何扩展自定义诊断函数。2. 从 sw_data.mat 到 sw_dens一条完整的水柱计算链路2.1 先看清 sw_data.mat 里的数据结构工具箱自带的 sw_data.mat 是一份标准示例剖面包含了 14 层标准压力上的温度、盐度以及观测纬度。在 MATLAB 里加载它只需要两行load sw_data.mat whoswhos会在命令行列出工作区里所有变量的名称、大小和类型。拿到这组变量后先别急着算确认一下压力和深度的对应关系disp([P, sw_dpth(P, LAT)])sw_dpth做的是压力到深度的换算输入压力 P 和纬度 LAT输出该压力对应的几何深度单位是米。为什么需要纬度因为重力加速度 g 随纬度和压力是变化的深度换算必须考虑这个因素。disp把压力列和深度列并列打印出来我一般先看一眼这个数组确认剖面没有出现压力倒置或者深度跳变再做下一步计算。这个工具箱的数据组织方式非常直接所有函数都接受变量数组作为输入数组的每一列代表一个观测站位的垂直剖面每一行代表一层标准压力。后面的密度、位温、位密计算全都遵循这个约定所以 sw_data.mat 里存的就是一个可以被所有 sw_ 前缀函数直接消费的标准输入格式。2.2 用 sw_dens 算原位密度先搞清楚 EOS-80 的输入约定拿到 S、T、P 三个剖面后最自然的第一个计算就是原位密度rho sw_dens(S, T, P);sw_dens输入是实用盐度 Spsu、原位温度 T摄氏度和原位压力 Pdbar输出是原位密度 rho单位是 kg/m³。这里有一个非常容易踩的坑S 必须是实用盐度而不是电导率原始读数T 必须是原位温度而不是位温P 必须是压力而不是深度。很多刚接触海洋数据的工程师直接把 CTD 输出的电导率塞进去算出来的密度差得离谱。EOS-80 的密度公式并不是摘要里写的 ρ ρ₀(1 αΔT βS) 那么简单那个线性近似只能用于教学演示。真实的 EOS-80 是在 0–40°C、0–42psu、0–10000dbar 范围内用十几项多项式拟合出来的温度、盐度、压力之间还有交叉项。工具箱先用sw_dens0对 35psu、0°C 的状态做参考态归一化再在sw_dens内部完成完整的非线性展开。α 和 β 也不是常系数而是随着温压状态变化的偏导数分别对应sw_alpha和sw_beta两个函数。2.3 位温和位密为什么不能拿原位密度做层结分析原位密度受到绝热压缩的影响一个 500m 深的水团被抬到海面压力释放之后密度变小但它的位密几乎不变。所以做层结稳定性分析时必须先把水团绝热移到参考压力面上再比较密度。工具箱里对应的是两个函数pt sw_ptmp(S, T, P, 0); % 位温参考压力取 0 dbar pden sw_pden(S, T, P, 0); % 位密参考压力取 0 dbarsw_ptmp的第四个参数是参考压力 PR传 0 表示把水团绝热移到海面sw_pden用同样的参考压力做位密计算。在物理海洋论文里pden通常写作 σ_θ也就是位密减掉 1000 之后的数值比如 1026.4 kg/m³ 写成 26.4 σ_θ。参考压力选多少取决于研究区域。深海环流研究习惯把参考压力选在 2000 dbar 甚至 4000 dbar用两个不同参考压力算出来的位密剖面差异可以在 0.1 kg/m³ 量级这个量级足以改变层结稳定性判断。下面这张表列出这一节用到的核心函数方便对照参数函数函数签名返回物理量单位sw_dens0sw_dens0(S)参考态密度kg/m³sw_pdensw_pden(S,T,P,PR)位密kg/m³sw_ptmpsw_ptmp(S,T,P,PR)位温°Csw_dpthsw_dpth(P,LAT)深度msw_pressw_pres(D,LAT)压力dbar表中sw_pres是sw_dpth的逆函数做单位换算时经常两个一起用。水深和压力在 1000m 以上相差不大1000 dbar 约等于 1000m但到了 5000m 深度压力数值比深度数值大 2% 左右不能直接画等号。3. sw_bfrq 与 sw_gvel把密度剖面换算成动力诊断量3.1 浮力频率稳定性分析的入口拿到位密剖面之后第一步动力学诊断通常是浮力频率。海洋学里用 N² 表示层结强度N 是 Brunt–Väisälä 频率单位是 s⁻¹工具箱里的 sw_bfrq 可以直接算[N2, midP] sw_bfrq(S, T, P, LAT); N sqrt(abs(N2)); % 取绝对值再开方避免负值报错sw_bfrq的输入是完整的三条剖面 S、T、P 和纬度 LAT纬度在这里用于计算当地重力加速度。输出N2是相邻压力层中点处的浮力频率平方单位是 s⁻²midP是对应的中间压力值。N2出现负值说明那一段水柱是静力不稳定的实际海洋里实时观测出现负值通常意味着发生了对流混合或者传感器噪声分析的时候要把负值区间单独标出来看。N² 的单位换算在工程里经常搞混。周期是 N 的倒数乘以 2π比如 N 0.01 s⁻¹ 对应周期约 10 分钟内波数值模式里通常直接用 N² 参与计算。我在处理潜标数据时习惯把 N² 画成对数坐标剖面因为海洋上层 N² 可以到 10⁻³ 量级深水区掉到 10⁻⁶ 量级线性坐标根本看不出深水层的微弱层结。3.2 地转速度从密度场到流场的关键一步浮力频率描述的是水柱的垂向稳定性地转流则利用密度场的水平梯度反演流速。sw_gvel做的是经典动力计算方法先沿等压面计算位势再做水平差分最后除以科氏参数得到地转速度剖面。调用方式是gvel sw_gvel(S, T, P, LAT);这里的 S、T、P 是二维数组每一列是一条观测剖面的垂直数据每一行代表同一个压力面。输出gvel是与输入同形状的速度矩阵单位是 m/s正值代表与差分方向一致的地转流动。实际使用中有一个限制必须讲清楚sw_gvel给出的只是相对流速。动力方法默认某一深层为零速度面如果两条站位之间没有直接测流数据这个参考面只能靠经验取。在大洋中上层研究里有人取 1000 dbar 作为零面有人取 2000 dbar取法不同速度剖面整体会平移但剪切结构不变。这也是为什么论文里地转流结果通常配一句相对 XX dbar 的流速。3.3 声速剖面多波束和定位里的实际用途密度剖面除了做层结和地转流还能直接换算成声速剖面用于多波束测深、水声定位和声学释放器的距离校正。工具箱里的函数是c sw_seck(S, T, P);输出 c 的单位是 m/s典型值在 1450–1550 m/s 之间。sw_seck用的是 Chen 与 Millero 1977 年的经验公式对温度和压力的灵敏度最高盐度的影响相对小但深海环境下盐度的累积误差同样会让声速产生 0.1 m/s 量级的偏差对长基线定位来说这个量级不能忽略。工具箱里同时存在sw_svel这个函数名它是旧脚本兼容入口新代码统一调sw_seck即可。用声速剖面做声学计算时常见做法是把 c 按深度积分得到声学路径距离再和直读距离做差校验收敛指标。4. 盐度从电导率到实用盐度sw_salt / sw_cndr / sw_salrp 的换算细节4.1 实用盐度不是一个浓度而是一个电导率比很多刚接触海洋数据分析的工程师会被盐度这个名字带偏以为它是海水中溶解盐类的质量浓度。实际上 PSS-78 实用盐度定义的是 15°C、1 个标准大气压下海水样品电导率与 KCl 标准溶液电导率的比值 K₁₅再经过一个非线性函数转换得到 0–42 之间的无量纲数值习惯上写成 psu 但严格意义上没有单位。工具箱里sw_salt就是电导率比到实用盐度的直接入口CNDR 0.85; % 仪器给出的电导率比原值 S sw_salt(CNDR, 25, 100);入参CNDR是电导率比25 是原位温度100 是原位压力函数内部完成温度补偿和压力补偿后返回实用盐度。问题在于CTD 探头给出的通常不是直接的电导率比而是经过标定系数换算过的电导率单位是 S/m。要得到CNDR必须先用海水电导率除以参考电导率这个参考值就是sw_c3515返回的 42.914 mS/cm。4.2 sw_cndr 与 sw_salrp反算与压力修正如果手里是盐度而不是电导率比又需要电导率比做质量控制可以用sw_cndr反算R sw_cndr(35, 25, 100);这里 R 是 25°C、100 dbar 条件下的电导率比函数内部先通过 PSS-78 的定义式反解出 15°C 参考条件下的比值再修正到目标温度和压力。sw_salrp则是单独拿出来做压力修正的中间函数它接受任意条件的电导率比 R、温度 T 和压力 P返回消除压力影响后的修正比值 Rx。实际数据处理里我一般在两个地方用到sw_cndr一个是把 CTD 原始电导率记录转成标准电导率比另一个是交叉验证不同传感器的盐度漂移。把历史盐度转成电导率比之后和原始电导率对比曲线的偏移量可以直接映射为探头污染导致的盐度漂移速率。4.3 盐度换算链路的核心函数表下面这张表把盐度相关的几个函数串起来标注了输入输出方向方便排查数据流函数输入 → 输出典型用途sw_c3515无 → 42.914取参考电导率常数sw_cndr(S,T,P) → R盐度反算电导率比sw_salrp(R,T,P) → Rx电导率比压力修正sw_salrt(R,T) → Rt电导率比温度修正sw_salt(CNDR,T,P) → S电导率比转实用盐度sw_salds(S,T,P) → dS盐度随压力变化率sw_salrt在工具包里的角色是把任意温度下的电导率比折到 15°C 参考条件下sw_salds则给出盐度对压力的导数用于计算等熵面梯度。搞懂了这几个函数的关系也就明白了为什么 CTD 数据手册里总强调盐度计算必须同时记录温度和压力电导率的温度系数高达每度 2% 左右压力修正虽然小但深水站的累积误差完全可以到 0.01 psu 量级。5. 用 sw_new 扩展海洋要素计算一个混合层深度诊断示例5.1 用 sw_new 生成标准文件头sw_new 是工具箱提供的扩展脚手架在 MATLAB 命令行运行sw_new按提示输入函数名sw_mld它会生成一个带标准注释段的新函数模板。自定义函数放进工具箱目录后要遵守工具箱的规则文件名以小写sw_开头输入参数顺序固定。生成的模板长这样function mld sw_mld(S, T, P) % SW_MLD Mixed layer depth based on density threshold. % The reference density is computed at 10 dbar. % Syntax: mld sw_mld(S, T, P) % Inputs: S salinity profile (psu) % T in-situ temperature profile (deg C) % P pressure profile (dbar) % Output: mld mixed layer depth (m) % % See also SW_DENS, SW_PDEN % Add your code below end这个模板的核心是保留了工具箱的统一调用约定后续所有脚本都能用和sw_dens相同的方式调用它。5.2 用密度阈值法计算混合层深度混合层深度的一个常见定义是从 10 dbar 参考密度出发向下找密度增量超过 0.03 kg/m³ 的第一个深度跨越点。基于这个定义填好的函数如下function mld sw_mld(S, T, P) rho0 sw_dens(S(1), T(1), P(1)); % 表层参考密度 rho sw_dens(S, T, P); % 原位密度剖面 dP abs(P - P(1)); target rho0 0.03; % 密度阈值 idx find(rho target dP 10, 1, first); if isempty(idx) mld NaN; return; end % 线性插值到阈值穿越点 mld interp1(rho(idx-1:idx), P(idx-1:idx), target); endrho0取表层第一个有效压力层的密度作为参考target在参考密度上加 0.03 kg/m³ 的密度差阈值。find找到第一个同时满足密度超过阈值且离表层至少 10 dbar 的层位interp1在这个层位和上一层之间做线性插值把阈值穿越的精确压力求出来。这里用压力近似水深输出混合层深度严格做法是把sw_dens换成sw_pden并用sw_dpth把压力转成深度。5.3 边界条件与验证方法混合层深度的阈值没有全球统一标准0.03 kg/m³ 是热带和副热带常见配置高纬度海域由于密度层结弱很多人改用 0.01–0.02 kg/m³。验证自定义函数是否可靠可以把它和sw_bfrq的结果对照混合层底的位置应该对应 N² 第一次明显增强的深度面两者相差超过 20% 时优先检查参考面选取是否合理。这个例子说明sw_new 生成的骨架可以快速把一套工程判断固化成工具箱的标准函数后续批量处理历史 CTD 数据时直接调用sw_mld(S, T, P)就能对每个站位输出统一的混合层深度不再需要复制粘贴脚本里的阈值逻辑。本文还有配套的精品资源点击获取