简介一份基于Matlab t_tide工具包的潮位调和分析脚本聚焦潮汐验证与潮流分析适用于海洋工程、港口航道规划、海洋环境监测等领域的科研人员、工程师及高校相关专业学生。压缩包内仅含1个tiaohe.m源文件体积约2KB代码结构紧凑重点演示了利用t_tide进行潮位时间序列调和分析的核心流程从原始观测数据导入、缺失与异常值预处理、主要分潮如M2、S2、N2等的选取到最小二乘参数估计、残差检验以及潮位拟合与预测的可视化输出并给出常用精度评价指标。脚本可直接替换数据运行也可作为二次开发模板便于使用者理解调和常数计算逻辑、潮汐动力机制及模型精度验证方法。已有581人学习下载对于希望快速上手潮汐调和分析、减少重复编程工作的入门者与进阶者这份脚本提供了清晰、可复用的代码参照与实际范例。1. 从一段实测潮位到一张可用的调和分析结果表港口设计、航道疏浚、海洋环境监测甚至 offshore 工程的作业窗口规划都绕不开一个问题这段潮位或海流数据里到底藏着哪些周期性分量每个分量的振幅和相位是多少。t_tide 这套在海洋学里用了二十多年的 Matlab 工具包配合 tiaohe.m 这样的调用脚本解决的就是这件事。它不靠“猜趋势”而是把非平稳的潮位时间序列强制分解成一组固定频率的分潮M2、S2、K1、O1 等再用最小二乘拟合出每个分潮的振幅和相位最后给出预测序列、残差和置信区间。适合手里有实测水位或流速数据、需要出调和常数或者做潮汐预报的从业者也适合给论文里的“潮位验证”章节提供可复现的处理流程。下面这套流程里有几个地方容易踩坑比如数据长度不够导致分潮混淆、缺失值直接喂给 t_tide 导致 NaN 传导到所有输出这些问题我会在代码段里一并处理掉。2. 调和分析的数学基础与 t_tide 的输入输出约定2.1 潮位序列为什么能拆成有限项余弦波叠加潮位观测值 h(t) 在调和分析框架下被表示成平均海平面加上若干分潮余弦项的叠加再加一个非潮汐残差项。每个分潮都对应一个由天体引潮势推导出的固定角速度 ωi这个角速度是天文常数不随地点改变。模型写成h(t) Z0 Σ fi·Hi·cos(ωi·t (V0u)i − Gi) R(t)其中 Z0 是平均海平面Hi 和 Gi 是第 i 个分潮的振幅和迟角需要求解的核心参数fi 和 (V0u)i 是交点因子和交点订正角由天文参数推算t_tide 内部会自动计算。R(t) 是残差包含气象强迫、非线性效应和观测噪声。所有分潮频率两两之间的最小频率差决定了理论上能分辨的最短数据长度这就是 Rayleigh 准则两个分潮要想被分离至少需要 1/(ωi−ωj) 长度的数据。例如 K1 和 P1 的频率差约为 0.0000000729 rad/s对应约 182 天的数据才能完全分开。实际工作中如果只有一个月数据t_tide 的 default 配置会自动合并或剔除这类“混淆”的分潮这也是为什么同一组数据用不同 length 参数会得到不同分潮数量。2.2 t_tide 的函数族与核心参数解法t_tide 工具包的核心函数是t_tide它用最小二乘求解上面的线性方程。设计矩阵的每一列是 cos(ωi·t) 和 sin(ωi·t) 的节点因子修正形式未知量是每个分潮的余弦系数和正弦系数最后转换为振幅和迟角。求解完成后返回椭圆参数输出如果是流速分量输入、预测序列和误差估计。t_tide的输入是单列或双列时间序列双列对应 u、v 流速分量关键参数包括start指定起始时间、latitude传入测站纬度、output控制输出文件、rayleigh自定义 Rayleigh 准则倍数、synthesis用已有调和常数合成时间序列。t_tide内部用到一个核心子程序来处理从 T_xtide 全球潮汐模型插值得到的初始分潮列表如果不调用这个初始化函数比如直接运行旧版脚本分潮表会少很多浅水分潮结果差异明显。下面用表格列出 t_tide 中最常用的几个参数及其职责这些参数在 tiaohe.m 的改造过程中会反复用到。参数名取值示例作用典型使用场景start[2023 6 1 0 0 0]指定序列起始日期时间保证与真实时间对齐节点因子计算正确latitude30.5测站纬度用于交点因子精细修正高纬度测站必须设置否则分潮参数不准rayleigh1.0控制分潮分辨严格度短数据可放宽到 0.9长数据设 1.0 以上outputtest.out输出文本结果文件归档、批处理时常用synthesis1用已有常数合成潮位序列预测未来潮位、回代比较inferK1对弱分潮做潮汐推断数据不足时把 P1 从 K1 推断得到2.3 tiaohe.m 脚本应具备的数据流向tiaohe.m 作为调用 t_tide 的驱动脚本数据流向是原始时间序列 → 时间向量构造 → 缺失值处理 → 分潮配置 → t_tide 求解 → 结果提取与误差计算 → 绘图。每一步都有独立的检查点而不是把所有逻辑塞进一个函数。时间向量用datenum构造需要注意 Matlab 的 datenum 单位是天分潮角速度单位是 rad/st_tide 内部会做单位换算但如果自己写了合成序列的代码单位转换是一个常见 bug 来源。下面我给出一段 tiaohe.m 的骨架代码涵盖数据加载、预处理和 t_tide 调用。% tiaohe.m - 基于 t_tide 的潮位调和分析骨架 % 输入: 水位观测 CSV, 第一列为时间, 第二列为水位(单位: m) data readmatrix(tide_obs.csv); tnum datenum(data(:,1), data(:,2), data(:,3), ... data(:,4), data(:,5), data(:,6)); h data(:,7); % 缺失值线性插值, 超过10%连续缺失直接报错 if any(isnan(h)) if sum(isnan(h))/length(h) 0.1 error(缺失超过10%, 请重新检查原始数据); end h fillmissing(h, linear); end % 调用 t_tide 做调和分析, 返回主要结果 [coef, pred, err] t_tide(h, start, tnum(1), ... latitude, 30.5, ... output, tide_res.out, ... rayleigh, 1.0); fprintf(分潮数量: %d\n, length(coef.tidecon(:,1)));这段代码里coef.tidecon是一个四列矩阵第一列是分潮角速率、第二列是振幅、第三列是迟角、第四列是信噪比逐行对应一个分潮。pred是模型预测的潮位序列长度与输入一致err是每个分潮参数的标准差估计。readmatrix读取 CSV 时假设时间是分裂的六列数值如果原始数据是字符串时间应该改用datetime解析然后调用datenum转换。fillmissing做线性插值是快速方案如果数据有明确的天文潮周期成分可以先用低通滤波器分离出趋势再做插值效果会更好。3. 用 tiaohe.m 完成单站潮位调和分析与质量评估3.1 构造干净的时间序列并处理闰秒、时区和基准面问题单站潮位调和分析的第一步不是调用任何分析函数而是想清楚时间基准面。观测数据如果用的是东八区时间而 t_tide 默认按 UTC 处理结果差 8 小时对应的相位偏差K1 分潮的迟角会差 120 度。因此要把本地时间显式转成 UTC或者在调用 t_tide 时不做转换但把所有分潮迟角加上一个区域常数偏移。常见做法是在脚本开头统一转 UTC用datetime(..., TimeZone, Asia/Shanghai)构造本地时间再用TimeZoneUTCLeapSeconds转基准最后取datenum。这一段处理虽然和调和分析本身无关但直接决定迟角结果能不能与沿岸潮汐表对比。潮位基准面也需要对齐。t_tide 对绝对高程没有要求它拟合的是相对变化所以只要序列是等间距采样即可。但如果要输出“预报潮位”用于工程需要把平均海平面、理论最低潮面的偏移量叠加回去。具体做法是在调和分析完成后把pred序列加上原始序列的均值差pred_abs pred (mean_obs - mean_harmonic)其中mean_harmonic是分潮常数中平均海平面项的拟合值。这套流程下来预测序列的均值才能与观测一致否则 RMS 误差会比真实水平大一个常数项。3.2 分潮选取策略与 Rayleigh 准则冲突的处理t_tide 默认的分潮列表有 68 个分潮实际参与拟合的分潮个数由数据长度决定。一个常见误区是直接把rayleigh设为很小的值比如 0.5强行把本不该分离的分潮分开这样拟合的振幅和相位会出现高度相关参数误差会暴涨。我的经验是一个月数据用默认配置三个月数据可以把rayleigh调到 1.1 并加上浅水分潮列表半年以上数据再考虑 M4、MS4 等浅水分潮纳入主列表。下面是一段根据数据长度自动配置分潮的示例代码。% 根据数据时长自动配置分潮选取策略 dt tnum(2) - tnum(1); % 采样间隔, 单位: 天 N length(tnum); span_days N * dt; if span_days 35 opts {rayleigh, 0.8}; % 短数据放宽标准 elseif span_days 100 opts {rayleigh, 1.0}; % 标准配置 else opts {rayleigh, 1.2, shallow}; % 长数据增加浅水分潮 end [coef, pred, err] t_tide(h, start, tnum(1), ... latitude, 30.5, opts{:});dt是用相邻时间差计算的实际采样间隔不直接使用名义采样率这样可以发现数据里是否混入了重采样或者缺失插值造成的非等间隔问题。span_days判断数据跨度35 天以内允许rayleigh0.8是实践中的折中——此时 P1 和 K1 不会被强行分离但 S2 和 K2 能勉强分开。100 天以上加上shallow选项后t_tide 会在分潮列表中追加 35 个浅水分潮M4、MS4、MN4 等适用于近岸或河口站点。数据起点不是整点小时也没关系t_tide 会利用start参数里的具体时间计算天文初位相关键是不要对时间做“四舍五入到整点”这类预处理这会人为改变分潮初相位。3.3 结果验证残差谱密度、RMSE 与置信区间调和分析完成后验证工作主要看三个指标模型预测与观测的 RMSE、残差的频谱中是否还有显著潮汐能量、分潮参数的信噪比和误差区间。t_tide 直接给出err参数其中第四列是每个分潮的 SNR。按经验SNR 小于 1 的分潮基本不可信适合从结果表里剔除SNR 在 1 到 2 之间的分潮可以保留但要在报告中注明SNR 大于 5 的分潮是可靠的。下面用一段代码计算残差功率谱并检查分潮频率处是否还有“残余峰”。residual h - pred; % 计算残差功率谱 fs 1 / (dt * 86400); % 采样频率, 单位: Hz [pxx, f] pwelch(residual, [], [], [], fs); % 检查 M2 和 S2 频率处是否有残余谱峰 m2_freq 1.9322736e-5; % M2 角频率对应的 Hz idx find(abs(f - m2_freq) 1e-7); if ~isempty(idx) peak_m2 10*log10(pxx(idx)); fprintf(M2处残差谱密度: %.2f dB\n, peak_m2); end rmse sqrt(mean(residual.^2)); fprintf(RMSE: %.3f m\n, rmse);pwelch用的是 Welch 平均周期图法默认加汉宁窗频率分辨率由窗长决定。m2_freq的值是 1.9322736e-5 Hz对应 M2 分潮 12.4206 小时的周期。残差谱在 M2 处如果没有明显尖峰说明潮汐能量已经被模型充分提取如果有尖峰意味着该分潮没有被正确拟合常见原因是数据长度不足导致 M2 和 N2 混淆或者振幅变化被截断。rmse是整体精度的直接指标一般验潮站 30 天调和分析的 RMSE 在 0.1 到 0.2 米之间如果超过 0.3 米应优先检查是否有非潮汐水位变化风暴增水、径流影响存在于原始序列中。下面给出一个典型的分潮结果表展示 t_tide 输出的主要列和验证结论。分潮名称振幅 (m)迟角 (deg)SNR使用建议M2主要太阴半日分潮1.42218.5156.2主分潮必须保留S2主要太阳半日分潮0.38245.142.8保留K1太阴日分潮0.2297.35.6保留注意与 P1 的混淆O1主要太阴日分潮0.18112.44.1保留P1太阳日分潮0.0788.60.8SNR 低建议剔除或推断N2椭圆率效应分潮0.29205.72.3保留但误差区间较大4. 多站联测数据与准调和潮流椭圆分析4.1 流速分量进 t_tide 后输出的椭圆参数含义当输入是海流观测的 u、v 分量时t_tide 不再输出简单的振幅和迟角而是输出潮流椭圆参数长半轴长度最大流速、短半轴长度最小流速、长轴方向、相位和旋转方向。这些参数描述的是一段时间内潮流矢量的端点轨迹——在开阔海域通常是旋转流轨迹近似椭圆在近岸水道则是往复流短半轴趋近于零。t_tide 内的处理逻辑是对 u、v 序列分别做调和分析再把同频率的余弦与正弦系数组合成椭圆参数。因此多站联测数据的核心问题是保证各站时间基准一致否则不同站点的相位无法对比潮流场的空间分布会失真。实际工作中多站联测数据的来源可能是走航 ADCP 按断面分层采集也可能是多个底式海流计同步投放。前者每个深度单元的时间长度通常只有 1 到 3 天做准调和分析比全调和更现实后者如果有 15 天以上的连续记录可以做完整调和分析。准调和分析的做法是固定分潮频率通常是 M2、S2、K1、O1再加 M4 和 MS4不求解天文节点因子直接用观测序列拟合复振幅。t_tide 不直接提供准调和模式但可以通过设置rayleigh为极大值并手动指定分潮列表来近似实现。4.2 批量处理多站数据的结构化脚本设计多站数据的处理不能逐站手点脚本应该把所有文件统一命名、统一格式然后用一个循环批量处理。这里给出一段批量处理框架输出每个站的调和常数并汇总到一个结构体数组中方便后续绘制潮汐图或计算相位差。% 批量处理多站 ADCP 数据 files dir(station_*.mat); results []; % 用于保存所有站的调和常数 for k 1:length(files) load(fullfile(files(k).folder, files(k).name), u, v, t); % u, v 是流速分量, t 是时间向量, 单位: datenum % 剔除流速低于阈值的“静止”记录, 避免噪声主导 speed sqrt(u.^2 v.^2); valid speed 0.02; % 阈值 2 cm/s u u(valid); v v(valid); t t(valid); if length(t) 660 % 少于 660 点(约 11 小时1min采样) 则跳过 warning(站 %s 数据过短, 已跳过, files(k).name); continue; end % 调用 t_tide 椭圆模式, 输入两列矩阵 [coef, pred, err] t_tide([u(:) v(:)], ... start, t(1), latitude, 30.5, ... rayleigh, 1.0, output, [files(k).name .out]); % 提取 M2 椭圆参数 m2idx find(coef.freq 1.405189e-4); % M2 频率 if ~isempty(m2idx) results(k).station files(k).name; results(k).m2_major coef.ellipse(m2idx, 1); results(k).m2_minor coef.ellipse(m2idx, 2); results(k).m2_theta coef.ellipse(m2idx, 3); results(k).m2_phase coef.ellipse(m2idx, 4); end end % 保存汇总结果 save(harmonic_summary.mat, results);这段代码中coef.freq是分潮频率数组坐标系约定可以由t_tide的coor参数指定默认是“向东为正 u、向北为正 v”的海流坐标系。coef.ellipse的每一行对应一个分潮四列分别是长半轴、短半轴、长轴方位角从北顺时针和相位。speed 0.02的低速滤波是一种经验做法因为 ADCP 在近底层测得的微小流速常常是湍流伪影不参与调和拟合反而能减少噪声但这也意味着要保证过滤后的数据长度仍然覆盖至少一个完整潮周期。4.3 短时段观测的准调和分析与 t_tide 的限制边界如果观测窗口只有 25 小时或 48 小时完整调和分析无法分离 S2 与 K2、P1 与 K1 这些相近频率的分潮。此时更稳妥的做法是采用“固定差比法”利用邻近长期验潮站给出的振幅比和迟角差把次要分潮从主分潮中推断出来。t_tide 提供infer参数和t_tide内部支持的推断机制。不过对于流数据分析短时段的强流过程比如风暴潮叠加往往会让调和常数产生显著偏差所以短时段结果应当作为参考而不是最终工程结论这和潮位分析的标准做法是保持一致的。5. 直接从 t_tide 结果跳到数据对比和异常诊断的实用技巧调和分析的最后一步不是“跑通”而是“审结果”。这里给出一个我常用的诊断流程先用 t_tide 自带参数合成预测序列再计算逐时残差然后把残差做 24 小时滑动平均这样可以快速识别非潮汐水位变化事件。如果残差滑动平均序列在某个时段明显偏离零线对应的时间点大概率有 storm surge 或径流峰值这些时段的原始数据可以考虑做标记但不删除——因为调和分析要求连续等间隔采样删除会造成频谱泄漏。正确做法是保留数据在结果报告中单独说明这段时间的残差特征。另一个容易忽略的验证手段是分潮振幅比和迟角差的合理性检查。对同一测站M2/S2 的振幅比通常在 2.0 到 4.0 之间如果超出这个范围往往意味着分潮混淆或数据质量问题。将 M2 与 S2 的迟角差G_S2 − G_M2换算成时间差应约等于大潮周期14.77 天内的一个固定值。这个检查不需要额外数据用小脚本直接对比 t_tide 输出即可。下面给出一段迟角差检查代码。% 检查 M2 和 S2 分潮的迟角差是否合理 m2_idx find(coef.freq coef.freq(1)); % 需按实际频率定位 % 实际应以频率匹配为准, 这里示意提取 freqs coef.freq; m2f freqs(1); % 示意赋值, 实际用 1.405189e-4 s2f freqs(2); % 实际用 1.454441e-4 f_m2 find(abs(freqs - 1.405189e-4) 1e-6); f_s2 find(abs(freqs - 1.454441e-4) 1e-6); if ~isempty(f_m2) ~isempty(f_s2) phase_diff coef.tidecon(f_s2, 3) - coef.tidecon(f_m2, 3); phase_diff mod(phase_diff, 360); time_diff phase_diff / (360 / 12.4206); % 小时 fprintf(S2与M2相位差换算时间: %.2f h\n, time_diff); end迟角差换算成小时后的合理参考范围可以由理论推算但更实用的是与相邻验潮站多年调和常数对比偏差应在正负一小时以内。mod(phase_diff, 360)这步很重要因为 t_tide 输出的迟角范围是 0 到 360直接做差会出现负值或超过 360 的情况必须先做周期折叠再比较。最终交付的调和常数表应该附上置信区间来自err输出、数据时段、采样间隔滤波方式说明这样下游做工程设计的人才能判断数据边界在哪里。本文还有配套的精品资源点击获取 基于t_tide的Matlab潮汐调和分析全流程指南
简介一份基于Matlab t_tide工具包的潮位调和分析脚本聚焦潮汐验证与潮流分析适用于海洋工程、港口航道规划、海洋环境监测等领域的科研人员、工程师及高校相关专业学生。压缩包内仅含1个tiaohe.m源文件体积约2KB代码结构紧凑重点演示了利用t_tide进行潮位时间序列调和分析的核心流程从原始观测数据导入、缺失与异常值预处理、主要分潮如M2、S2、N2等的选取到最小二乘参数估计、残差检验以及潮位拟合与预测的可视化输出并给出常用精度评价指标。脚本可直接替换数据运行也可作为二次开发模板便于使用者理解调和常数计算逻辑、潮汐动力机制及模型精度验证方法。已有581人学习下载对于希望快速上手潮汐调和分析、减少重复编程工作的入门者与进阶者这份脚本提供了清晰、可复用的代码参照与实际范例。1. 从一段实测潮位到一张可用的调和分析结果表港口设计、航道疏浚、海洋环境监测甚至 offshore 工程的作业窗口规划都绕不开一个问题这段潮位或海流数据里到底藏着哪些周期性分量每个分量的振幅和相位是多少。t_tide 这套在海洋学里用了二十多年的 Matlab 工具包配合 tiaohe.m 这样的调用脚本解决的就是这件事。它不靠“猜趋势”而是把非平稳的潮位时间序列强制分解成一组固定频率的分潮M2、S2、K1、O1 等再用最小二乘拟合出每个分潮的振幅和相位最后给出预测序列、残差和置信区间。适合手里有实测水位或流速数据、需要出调和常数或者做潮汐预报的从业者也适合给论文里的“潮位验证”章节提供可复现的处理流程。下面这套流程里有几个地方容易踩坑比如数据长度不够导致分潮混淆、缺失值直接喂给 t_tide 导致 NaN 传导到所有输出这些问题我会在代码段里一并处理掉。2. 调和分析的数学基础与 t_tide 的输入输出约定2.1 潮位序列为什么能拆成有限项余弦波叠加潮位观测值 h(t) 在调和分析框架下被表示成平均海平面加上若干分潮余弦项的叠加再加一个非潮汐残差项。每个分潮都对应一个由天体引潮势推导出的固定角速度 ωi这个角速度是天文常数不随地点改变。模型写成h(t) Z0 Σ fi·Hi·cos(ωi·t (V0u)i − Gi) R(t)其中 Z0 是平均海平面Hi 和 Gi 是第 i 个分潮的振幅和迟角需要求解的核心参数fi 和 (V0u)i 是交点因子和交点订正角由天文参数推算t_tide 内部会自动计算。R(t) 是残差包含气象强迫、非线性效应和观测噪声。所有分潮频率两两之间的最小频率差决定了理论上能分辨的最短数据长度这就是 Rayleigh 准则两个分潮要想被分离至少需要 1/(ωi−ωj) 长度的数据。例如 K1 和 P1 的频率差约为 0.0000000729 rad/s对应约 182 天的数据才能完全分开。实际工作中如果只有一个月数据t_tide 的 default 配置会自动合并或剔除这类“混淆”的分潮这也是为什么同一组数据用不同 length 参数会得到不同分潮数量。2.2 t_tide 的函数族与核心参数解法t_tide 工具包的核心函数是t_tide它用最小二乘求解上面的线性方程。设计矩阵的每一列是 cos(ωi·t) 和 sin(ωi·t) 的节点因子修正形式未知量是每个分潮的余弦系数和正弦系数最后转换为振幅和迟角。求解完成后返回椭圆参数输出如果是流速分量输入、预测序列和误差估计。t_tide的输入是单列或双列时间序列双列对应 u、v 流速分量关键参数包括start指定起始时间、latitude传入测站纬度、output控制输出文件、rayleigh自定义 Rayleigh 准则倍数、synthesis用已有调和常数合成时间序列。t_tide内部用到一个核心子程序来处理从 T_xtide 全球潮汐模型插值得到的初始分潮列表如果不调用这个初始化函数比如直接运行旧版脚本分潮表会少很多浅水分潮结果差异明显。下面用表格列出 t_tide 中最常用的几个参数及其职责这些参数在 tiaohe.m 的改造过程中会反复用到。参数名取值示例作用典型使用场景start[2023 6 1 0 0 0]指定序列起始日期时间保证与真实时间对齐节点因子计算正确latitude30.5测站纬度用于交点因子精细修正高纬度测站必须设置否则分潮参数不准rayleigh1.0控制分潮分辨严格度短数据可放宽到 0.9长数据设 1.0 以上outputtest.out输出文本结果文件归档、批处理时常用synthesis1用已有常数合成潮位序列预测未来潮位、回代比较inferK1对弱分潮做潮汐推断数据不足时把 P1 从 K1 推断得到2.3 tiaohe.m 脚本应具备的数据流向tiaohe.m 作为调用 t_tide 的驱动脚本数据流向是原始时间序列 → 时间向量构造 → 缺失值处理 → 分潮配置 → t_tide 求解 → 结果提取与误差计算 → 绘图。每一步都有独立的检查点而不是把所有逻辑塞进一个函数。时间向量用datenum构造需要注意 Matlab 的 datenum 单位是天分潮角速度单位是 rad/st_tide 内部会做单位换算但如果自己写了合成序列的代码单位转换是一个常见 bug 来源。下面我给出一段 tiaohe.m 的骨架代码涵盖数据加载、预处理和 t_tide 调用。% tiaohe.m - 基于 t_tide 的潮位调和分析骨架 % 输入: 水位观测 CSV, 第一列为时间, 第二列为水位(单位: m) data readmatrix(tide_obs.csv); tnum datenum(data(:,1), data(:,2), data(:,3), ... data(:,4), data(:,5), data(:,6)); h data(:,7); % 缺失值线性插值, 超过10%连续缺失直接报错 if any(isnan(h)) if sum(isnan(h))/length(h) 0.1 error(缺失超过10%, 请重新检查原始数据); end h fillmissing(h, linear); end % 调用 t_tide 做调和分析, 返回主要结果 [coef, pred, err] t_tide(h, start, tnum(1), ... latitude, 30.5, ... output, tide_res.out, ... rayleigh, 1.0); fprintf(分潮数量: %d\n, length(coef.tidecon(:,1)));这段代码里coef.tidecon是一个四列矩阵第一列是分潮角速率、第二列是振幅、第三列是迟角、第四列是信噪比逐行对应一个分潮。pred是模型预测的潮位序列长度与输入一致err是每个分潮参数的标准差估计。readmatrix读取 CSV 时假设时间是分裂的六列数值如果原始数据是字符串时间应该改用datetime解析然后调用datenum转换。fillmissing做线性插值是快速方案如果数据有明确的天文潮周期成分可以先用低通滤波器分离出趋势再做插值效果会更好。3. 用 tiaohe.m 完成单站潮位调和分析与质量评估3.1 构造干净的时间序列并处理闰秒、时区和基准面问题单站潮位调和分析的第一步不是调用任何分析函数而是想清楚时间基准面。观测数据如果用的是东八区时间而 t_tide 默认按 UTC 处理结果差 8 小时对应的相位偏差K1 分潮的迟角会差 120 度。因此要把本地时间显式转成 UTC或者在调用 t_tide 时不做转换但把所有分潮迟角加上一个区域常数偏移。常见做法是在脚本开头统一转 UTC用datetime(..., TimeZone, Asia/Shanghai)构造本地时间再用TimeZoneUTCLeapSeconds转基准最后取datenum。这一段处理虽然和调和分析本身无关但直接决定迟角结果能不能与沿岸潮汐表对比。潮位基准面也需要对齐。t_tide 对绝对高程没有要求它拟合的是相对变化所以只要序列是等间距采样即可。但如果要输出“预报潮位”用于工程需要把平均海平面、理论最低潮面的偏移量叠加回去。具体做法是在调和分析完成后把pred序列加上原始序列的均值差pred_abs pred (mean_obs - mean_harmonic)其中mean_harmonic是分潮常数中平均海平面项的拟合值。这套流程下来预测序列的均值才能与观测一致否则 RMS 误差会比真实水平大一个常数项。3.2 分潮选取策略与 Rayleigh 准则冲突的处理t_tide 默认的分潮列表有 68 个分潮实际参与拟合的分潮个数由数据长度决定。一个常见误区是直接把rayleigh设为很小的值比如 0.5强行把本不该分离的分潮分开这样拟合的振幅和相位会出现高度相关参数误差会暴涨。我的经验是一个月数据用默认配置三个月数据可以把rayleigh调到 1.1 并加上浅水分潮列表半年以上数据再考虑 M4、MS4 等浅水分潮纳入主列表。下面是一段根据数据长度自动配置分潮的示例代码。% 根据数据时长自动配置分潮选取策略 dt tnum(2) - tnum(1); % 采样间隔, 单位: 天 N length(tnum); span_days N * dt; if span_days 35 opts {rayleigh, 0.8}; % 短数据放宽标准 elseif span_days 100 opts {rayleigh, 1.0}; % 标准配置 else opts {rayleigh, 1.2, shallow}; % 长数据增加浅水分潮 end [coef, pred, err] t_tide(h, start, tnum(1), ... latitude, 30.5, opts{:});dt是用相邻时间差计算的实际采样间隔不直接使用名义采样率这样可以发现数据里是否混入了重采样或者缺失插值造成的非等间隔问题。span_days判断数据跨度35 天以内允许rayleigh0.8是实践中的折中——此时 P1 和 K1 不会被强行分离但 S2 和 K2 能勉强分开。100 天以上加上shallow选项后t_tide 会在分潮列表中追加 35 个浅水分潮M4、MS4、MN4 等适用于近岸或河口站点。数据起点不是整点小时也没关系t_tide 会利用start参数里的具体时间计算天文初位相关键是不要对时间做“四舍五入到整点”这类预处理这会人为改变分潮初相位。3.3 结果验证残差谱密度、RMSE 与置信区间调和分析完成后验证工作主要看三个指标模型预测与观测的 RMSE、残差的频谱中是否还有显著潮汐能量、分潮参数的信噪比和误差区间。t_tide 直接给出err参数其中第四列是每个分潮的 SNR。按经验SNR 小于 1 的分潮基本不可信适合从结果表里剔除SNR 在 1 到 2 之间的分潮可以保留但要在报告中注明SNR 大于 5 的分潮是可靠的。下面用一段代码计算残差功率谱并检查分潮频率处是否还有“残余峰”。residual h - pred; % 计算残差功率谱 fs 1 / (dt * 86400); % 采样频率, 单位: Hz [pxx, f] pwelch(residual, [], [], [], fs); % 检查 M2 和 S2 频率处是否有残余谱峰 m2_freq 1.9322736e-5; % M2 角频率对应的 Hz idx find(abs(f - m2_freq) 1e-7); if ~isempty(idx) peak_m2 10*log10(pxx(idx)); fprintf(M2处残差谱密度: %.2f dB\n, peak_m2); end rmse sqrt(mean(residual.^2)); fprintf(RMSE: %.3f m\n, rmse);pwelch用的是 Welch 平均周期图法默认加汉宁窗频率分辨率由窗长决定。m2_freq的值是 1.9322736e-5 Hz对应 M2 分潮 12.4206 小时的周期。残差谱在 M2 处如果没有明显尖峰说明潮汐能量已经被模型充分提取如果有尖峰意味着该分潮没有被正确拟合常见原因是数据长度不足导致 M2 和 N2 混淆或者振幅变化被截断。rmse是整体精度的直接指标一般验潮站 30 天调和分析的 RMSE 在 0.1 到 0.2 米之间如果超过 0.3 米应优先检查是否有非潮汐水位变化风暴增水、径流影响存在于原始序列中。下面给出一个典型的分潮结果表展示 t_tide 输出的主要列和验证结论。分潮名称振幅 (m)迟角 (deg)SNR使用建议M2主要太阴半日分潮1.42218.5156.2主分潮必须保留S2主要太阳半日分潮0.38245.142.8保留K1太阴日分潮0.2297.35.6保留注意与 P1 的混淆O1主要太阴日分潮0.18112.44.1保留P1太阳日分潮0.0788.60.8SNR 低建议剔除或推断N2椭圆率效应分潮0.29205.72.3保留但误差区间较大4. 多站联测数据与准调和潮流椭圆分析4.1 流速分量进 t_tide 后输出的椭圆参数含义当输入是海流观测的 u、v 分量时t_tide 不再输出简单的振幅和迟角而是输出潮流椭圆参数长半轴长度最大流速、短半轴长度最小流速、长轴方向、相位和旋转方向。这些参数描述的是一段时间内潮流矢量的端点轨迹——在开阔海域通常是旋转流轨迹近似椭圆在近岸水道则是往复流短半轴趋近于零。t_tide 内的处理逻辑是对 u、v 序列分别做调和分析再把同频率的余弦与正弦系数组合成椭圆参数。因此多站联测数据的核心问题是保证各站时间基准一致否则不同站点的相位无法对比潮流场的空间分布会失真。实际工作中多站联测数据的来源可能是走航 ADCP 按断面分层采集也可能是多个底式海流计同步投放。前者每个深度单元的时间长度通常只有 1 到 3 天做准调和分析比全调和更现实后者如果有 15 天以上的连续记录可以做完整调和分析。准调和分析的做法是固定分潮频率通常是 M2、S2、K1、O1再加 M4 和 MS4不求解天文节点因子直接用观测序列拟合复振幅。t_tide 不直接提供准调和模式但可以通过设置rayleigh为极大值并手动指定分潮列表来近似实现。4.2 批量处理多站数据的结构化脚本设计多站数据的处理不能逐站手点脚本应该把所有文件统一命名、统一格式然后用一个循环批量处理。这里给出一段批量处理框架输出每个站的调和常数并汇总到一个结构体数组中方便后续绘制潮汐图或计算相位差。% 批量处理多站 ADCP 数据 files dir(station_*.mat); results []; % 用于保存所有站的调和常数 for k 1:length(files) load(fullfile(files(k).folder, files(k).name), u, v, t); % u, v 是流速分量, t 是时间向量, 单位: datenum % 剔除流速低于阈值的“静止”记录, 避免噪声主导 speed sqrt(u.^2 v.^2); valid speed 0.02; % 阈值 2 cm/s u u(valid); v v(valid); t t(valid); if length(t) 660 % 少于 660 点(约 11 小时1min采样) 则跳过 warning(站 %s 数据过短, 已跳过, files(k).name); continue; end % 调用 t_tide 椭圆模式, 输入两列矩阵 [coef, pred, err] t_tide([u(:) v(:)], ... start, t(1), latitude, 30.5, ... rayleigh, 1.0, output, [files(k).name .out]); % 提取 M2 椭圆参数 m2idx find(coef.freq 1.405189e-4); % M2 频率 if ~isempty(m2idx) results(k).station files(k).name; results(k).m2_major coef.ellipse(m2idx, 1); results(k).m2_minor coef.ellipse(m2idx, 2); results(k).m2_theta coef.ellipse(m2idx, 3); results(k).m2_phase coef.ellipse(m2idx, 4); end end % 保存汇总结果 save(harmonic_summary.mat, results);这段代码中coef.freq是分潮频率数组坐标系约定可以由t_tide的coor参数指定默认是“向东为正 u、向北为正 v”的海流坐标系。coef.ellipse的每一行对应一个分潮四列分别是长半轴、短半轴、长轴方位角从北顺时针和相位。speed 0.02的低速滤波是一种经验做法因为 ADCP 在近底层测得的微小流速常常是湍流伪影不参与调和拟合反而能减少噪声但这也意味着要保证过滤后的数据长度仍然覆盖至少一个完整潮周期。4.3 短时段观测的准调和分析与 t_tide 的限制边界如果观测窗口只有 25 小时或 48 小时完整调和分析无法分离 S2 与 K2、P1 与 K1 这些相近频率的分潮。此时更稳妥的做法是采用“固定差比法”利用邻近长期验潮站给出的振幅比和迟角差把次要分潮从主分潮中推断出来。t_tide 提供infer参数和t_tide内部支持的推断机制。不过对于流数据分析短时段的强流过程比如风暴潮叠加往往会让调和常数产生显著偏差所以短时段结果应当作为参考而不是最终工程结论这和潮位分析的标准做法是保持一致的。5. 直接从 t_tide 结果跳到数据对比和异常诊断的实用技巧调和分析的最后一步不是“跑通”而是“审结果”。这里给出一个我常用的诊断流程先用 t_tide 自带参数合成预测序列再计算逐时残差然后把残差做 24 小时滑动平均这样可以快速识别非潮汐水位变化事件。如果残差滑动平均序列在某个时段明显偏离零线对应的时间点大概率有 storm surge 或径流峰值这些时段的原始数据可以考虑做标记但不删除——因为调和分析要求连续等间隔采样删除会造成频谱泄漏。正确做法是保留数据在结果报告中单独说明这段时间的残差特征。另一个容易忽略的验证手段是分潮振幅比和迟角差的合理性检查。对同一测站M2/S2 的振幅比通常在 2.0 到 4.0 之间如果超出这个范围往往意味着分潮混淆或数据质量问题。将 M2 与 S2 的迟角差G_S2 − G_M2换算成时间差应约等于大潮周期14.77 天内的一个固定值。这个检查不需要额外数据用小脚本直接对比 t_tide 输出即可。下面给出一段迟角差检查代码。% 检查 M2 和 S2 分潮的迟角差是否合理 m2_idx find(coef.freq coef.freq(1)); % 需按实际频率定位 % 实际应以频率匹配为准, 这里示意提取 freqs coef.freq; m2f freqs(1); % 示意赋值, 实际用 1.405189e-4 s2f freqs(2); % 实际用 1.454441e-4 f_m2 find(abs(freqs - 1.405189e-4) 1e-6); f_s2 find(abs(freqs - 1.454441e-4) 1e-6); if ~isempty(f_m2) ~isempty(f_s2) phase_diff coef.tidecon(f_s2, 3) - coef.tidecon(f_m2, 3); phase_diff mod(phase_diff, 360); time_diff phase_diff / (360 / 12.4206); % 小时 fprintf(S2与M2相位差换算时间: %.2f h\n, time_diff); end迟角差换算成小时后的合理参考范围可以由理论推算但更实用的是与相邻验潮站多年调和常数对比偏差应在正负一小时以内。mod(phase_diff, 360)这步很重要因为 t_tide 输出的迟角范围是 0 到 360直接做差会出现负值或超过 360 的情况必须先做周期折叠再比较。最终交付的调和常数表应该附上置信区间来自err输出、数据时段、采样间隔滤波方式说明这样下游做工程设计的人才能判断数据边界在哪里。本文还有配套的精品资源点击获取