Matlab seawater 工具包实战基于 CTD 剖面数据计算 3 种地转流速与位势异常海洋科学研究中CTDConductivity-Temperature-Depth剖面数据是最基础也最重要的观测资料之一。如何从这些原始数据中提取出有物理意义的海洋动力学参数是每个海洋科研人员都会面临的挑战。本文将带你使用 Matlab 的 seawater 工具包一步步实现从 CTD 数据到位势异常和地转流速的完整计算流程。1. 环境准备与数据加载在开始之前我们需要确保 seawater 工具包已正确安装。这个工具包由 CSIRO 开发包含了计算海水各种物理性质的函数。安装完成后可以通过sw_check函数验证安装是否成功% 检查 seawater 工具包安装 sw_check接下来加载 CTD 剖面数据。假设我们有一个名为ctd_data.mat的文件包含了温度、盐度和深度信息% 加载 CTD数据 load(ctd_data.mat); % 查看数据结构 whos典型的 CTD 数据结构应包含以下字段depth: 深度向量单位米temp: 温度剖面单位°Csal: 盐度剖面单位PSUlat: 观测点纬度用于地转计算提示实际工作中CTD 数据可能来自不同格式如 .cnv, .asc 等需要先进行数据清洗和质量控制。2. 数据预处理与质量检查原始 CTD 数据通常需要进行以下预处理步骤% 去除NaN值 valid_idx ~isnan(depth) ~isnan(temp) ~isnan(sal); depth depth(valid_idx); temp temp(valid_idx); sal sal(valid_idx); % 按深度排序 [depth, sort_idx] sort(depth); temp temp(sort_idx); sal sal(sort_idx); % 计算压力db pressure sw_pres(depth, lat);数据质量检查是必不可少的环节。我们可以绘制基本剖面图进行视觉检查figure; subplot(1,3,1); plot(temp, -depth); ylabel(Depth (m)); xlabel(Temperature (°C)); grid on; subplot(1,3,2); plot(sal, -depth); xlabel(Salinity (PSU)); grid on; subplot(1,3,3); plot(sw_dens(sal,temp,pressure), -depth); xlabel(Density (kg/m^3)); grid on;3. 位势异常计算位势异常Geopotential Anomaly是计算地转流的基础物理量反映了水柱的势能差异。使用sw_gpan函数计算% 计算位势异常 gpan sw_gpan(sal, temp, pressure); % 可视化结果 figure; plot(gpan, -depth); xlabel(Geopotential Anomaly (m^2/s^2)); ylabel(Depth (m)); title(Geopotential Anomaly Profile); grid on;位势异常计算中需要注意以下几点输入的压力应为分压db盐度应为实用盐度PSU温度应为现场温度°C4. 三种地转流速计算方法地转流速的计算有多种方法下面介绍三种常用方法及其实现。4.1 相对地转流速常规方法% 计算相对地转流速 [vg, xg] sw_gvel(pressure, gpan, lat); % 可视化 figure; plot(vg, -depth); xlabel(Geostrophic Velocity (m/s)); ylabel(Depth (m)); title(Relative Geostrophic Velocity); grid on;4.2 绝对地转流速需参考层% 假设1500db为无运动层 ref_depth 1500; [~, ref_idx] min(abs(pressure - ref_depth)); % 计算绝对地转流速 vg_abs vg - vg(ref_idx); % 可视化 figure; plot(vg_abs, -depth); xlabel(Absolute Geostrophic Velocity (m/s)); ylabel(Depth (m)); title([Absolute Geostrophic Velocity (Reference: num2str(ref_depth) db)]); grid on;4.3 密度坐标下的地转流速% 计算位势密度 sigma_theta sw_pden(sal, temp, pressure, 0); % 在密度坐标下插值 sigma_levels linspace(min(sigma_theta), max(sigma_theta), 100); gpan_sigma interp1(sigma_theta, gpan, sigma_levels); pressure_sigma interp1(sigma_theta, pressure, sigma_levels); % 计算密度坐标下的地转流速 [vg_sigma, ~] sw_gvel(pressure_sigma, gpan_sigma, lat); % 可视化 figure; plot(vg_sigma, -pressure_sigma); xlabel(Geostrophic Velocity (m/s)); ylabel(Pressure (db)); title(Geostrophic Velocity in Density Coordinates); grid on;5. 结果可视化与分析将计算结果综合展示可以更全面地理解海洋动力结构figure; % 温度、盐度、密度剖面 subplot(1,4,1); plot(temp, -depth); ylabel(Depth (m)); xlabel(Temp (°C)); grid on; subplot(1,4,2); plot(sal, -depth); xlabel(Salinity (PSU)); grid on; subplot(1,4,3); plot(sw_dens(sal,temp,pressure), -depth); xlabel(Density (kg/m^3)); grid on; % 地转流速剖面 subplot(1,4,4); plot(vg, -depth, r, LineWidth, 1.5); hold on; plot(vg_abs, -depth, b--, LineWidth, 1.5); xlabel(Velocity (m/s)); legend(Relative, Absolute); grid on;对于更专业的分析可以绘制温盐图和地转流速矢量图% 温盐图 figure; plot(sal, temp, .); xlabel(Salinity (PSU)); ylabel(Temperature (°C)); title(TS Diagram); grid on; % 地转流速矢量图假设有多个站位 figure; quiver(lon, depth, vg, zeros(size(vg)), 0.5); xlabel(Longitude); ylabel(Depth (m)); title(Geostrophic Velocity Vectors);6. 实际应用中的注意事项在实际科研工作中使用 seawater 工具包计算地转流速时需要注意以下问题参考层选择无运动层level of no motion的选择直接影响绝对地转流速结果可通过水团分析或现场ADCP数据验证数据质量控制% 示例检测密度倒置 density sw_dens(sal, temp, pressure); if any(diff(density) 0) warning(Density inversion detected! Check data quality.); end计算效率优化对于大范围数据可预先计算并保存中间结果使用向量化操作代替循环与其他数据集的对比验证将计算结果与卫星高度计数据、漂流浮标数据对比使用不同方法交叉验证结果可靠性下表总结了三种地转流速计算方法的特点和适用场景方法类型优点局限性典型应用场景相对地转流速计算简单不依赖假设只能反映流速相对变化断面分析水团边界识别绝对地转流速提供绝对值估计依赖参考层假设流量估算动力高度计算密度坐标流速物理意义明确需要高质量密度数据水团分析等密度面研究7. 扩展应用与高级技巧掌握了基本计算方法后可以进一步探索以下高级应用多站位的断面分析% 假设有多个站位的CTD数据 num_stations 5; vg_matrix zeros(length(depth), num_stations); for i 1:num_stations % 加载每个站位数据 station_data load([station_ num2str(i) .mat]); % 计算地转流速 [vg_matrix(:,i), ~] sw_gvel(station_data.pressure, ... sw_gpan(station_data.sal, station_data.temp, station_data.pressure), ... station_data.lat); end % 绘制断面图 figure; contourf(1:num_stations, -depth, vg_matrix, 20, LineColor, none); colorbar; xlabel(Station Number); ylabel(Depth (m)); title(Geostrophic Velocity Section);时间序列分析% 假设有时间序列CTD数据 dates datetime(2020,1,1):days(7):datetime(2020,12,31); vg_time zeros(length(depth), length(dates)); for t 1:length(dates) % 加载每周数据 weekly_data load([ctd_ datestr(dates(t), yyyymmdd) .mat]); % 计算地转流速 vg_time(:,t) sw_gvel(weekly_data.pressure, ... sw_gpan(weekly_data.sal, weekly_data.temp, weekly_data.pressure), ... weekly_data.lat); end % 绘制Hovmöller图 figure; contourf(dates, -depth, vg_time, 20, LineColor, none); colorbar; xlabel(Date); ylabel(Depth (m)); title(Geostrophic Velocity Time Series);与卫星高度计数据融合% 假设已加载卫星数据 satellite_vg load(satellite_geostrophic.mat); % 插值到CTD剖面深度 vg_sat_interp interp1(satellite_vg.depth, satellite_vg.velocity, depth); % 比较 figure; plot(vg_abs, -depth, b, LineWidth, 2); hold on; plot(vg_sat_interp, -depth, r--, LineWidth, 2); xlabel(Velocity (m/s)); ylabel(Depth (m)); legend(CTD-based, Satellite-derived); grid on;