t_tide 潮流调和分析实战：从水位数据到潮汐预报

简介基于Matlab的t_tide工具包R. Pawlowicz进行潮流调和分析的程序集合面向物理海洋学研究人员、海洋工程技术人员以及学习潮流要素计算的本科生和研究生。原始t_tide工具包被小幅修正并补充了部分功能同时配有操作说明文档可帮助读者快速完成潮流调和分析、潮汐要素计算与海洋数据处理。压缩包共31个文件大小约10.45MB其中包含18个m脚本文件、5个mat数据文件以及xls表格、dat数据、doc说明文档等覆盖主程序、误差估计、潮汐预测、调和常数提取等核心模块目录结构清晰便于按需调用。资源已有1683人学习下载适用范围覆盖科研与教学场景。借助该资源读者能够掌握t_tide的完整使用流程理解潮流调和分析的基本原理获得修正后的可运行代码、示例数据以及中英文操作指引节省调试时间提升分析效率。1. 潮流分析不靠猜t_tide 调和分析程序包能算出什么手里有一整年逐小时水位计数据想拆出潮汐成分、做预报却不知道该从哪里下手——这是很多物理海洋从业人员都撞过的墙。t_tide 潮流分析程序包是一个基于 MATLAB 的经典调和分析工具它把一段水位时间序列分解成若干天文分潮和浅水分潮的叠加输出每个分潮的振幅、迟角、信噪比以及合成潮位曲线。简单说只要给它一段均匀采样的水位序列和对应时间戳它就能告诉你这一带海域到底是半日潮主导还是全日潮主导M2 分潮有多大某个观测站位每天的潮差大概能到多少。这套工具适合三类人做物理海洋学研究需要计算潮流要素的学生做海洋工程需要设计基准面和水位设计值的工程师以及手里攒了大量水位原始数据、想系统提取潮汐信息的数据处理人员。它不需要你懂傅里叶变换的完整推导但需要你按它的规矩喂数据。接下来我按实际跑通流程把数据准备、参数配置、结果解读和几个容易翻车的地方一次讲清楚。2. 输入数据与配置时间序列规格、分潮选取与解析约束2.1 水位时间序列的三条硬规矩均匀采样、长度达标、基准面统一t_tide 的核心输入是两条向量时间t和水位x。时间必须是 MATLAB 的datenum格式单位为天水位单位为米。这三条规矩缺一条后面全白搭。第一采样必须均匀。t_tide 内部假设时间序列是等间隔的它通过参数dt以天为单位的采样间隔和序列长度来构造设计矩阵。如果你的原始数据是 10 分钟一条中间断了 3 小时又补了一段 20 分钟间隔的数据直接丢进去会怎样设计矩阵里的相位计算是按均匀间隔推进的实际数据间隔不均等于相位误差随时间累积出来的调和常数会系统性偏掉残差里也会残留明显的锯齿状信号。第二序列长度要能满足分潮分离的需要。这是整个调和分析里最容易被低估的约束。两个分潮的频率差决定分离它们需要的最短数据长度。以最常见的 M2周期 12.4206 小时和 S2周期 12.0000 小时为例两者频率差对应的拍频周期约为 14.77 天要稳定分离这两兄弟你至少需要约 15 天连续数据。要分离 K123.9345 小时和 P124.0659 小时则需要约 182 天。所以如果你只跑一个月的数据就不要指望 K1 和 P1 能同时给出可靠的独立解。第三水位基准面要统一。上游传来的数据可能混合了两种基准一部分以理论深度基准面为零点另一部分以平均海平面为零点。t_tide 不关心基准面的绝对定义它只对序列做调和拟合但如果序列中间换了基准面等于人为制造了一个阶跃这个阶跃会被当成一个极大周期的虚假信号吸收进拟合里污染邻近分潮的振幅。我一般在预处理阶段对整段序列做一次去平均再检查有没有明显跳变。常见做法是在喂给 t_tide 之前把时间先转成datenum并强制要求dt与数据实际时间差一致。下面这段预处理代码是每次跑分析的固定开场% 假设 raw_t 是 datetime 数组raw_x 是水位序列单位米 t_num datenum(raw_t); % 转成 datenum单位为天 dt0 median(diff(t_num)); % 计算实际中值采样间隔天 % 检查采样间隔是否稳定最大间隔与中值的偏差超过 1% 就要先插值 if max(abs(diff(t_num) - dt0)) 0.01 * dt0 t_reg t_num(1):dt0:t_num(end); % 生成均匀时间轴 x_reg interp1(t_num, raw_x, t_reg, linear); % 线性插值补全 else t_reg t_num; x_reg raw_x; end % 去平均并去掉明显超出物理范围的尖峰 x_reg x_reg - mean(x_reg); valid isfinite(x_reg); t_reg t_reg(valid); x_reg x_reg(valid);这里median(diff(t_num))取的是实际采样间隔的中值而不是直接相信设备配置里的采样率——设备偶尔跳秒、丢包是常态中值比标称值更接近真实间隔。插值用linear而不是spline原因在于 spline 会在缺测段附近制造过冲对潮汐这种窄带信号来说过冲会直接变成高频伪信号。去平均这一步不能省因为 t_tide 默认序列是零均值的常数项由拟合里的平均海平面项吸收但你提前去掉更稳。2.2 分潮不能想加就加奈奎斯特频率与 14.77 天分辨极限t_tide 不是把所有分潮一次性丢进拟合的。它内部有一套分潮列表默认情况下根据你的数据长度和采样间隔自动选出一组可解的分潮。这里有两个物理边界你必须心里有数。第一个边界是奈奎斯特频率。如果采样间隔是 1 小时那么你能分辨的最高频率周期是 2 小时任何周期短于 2 小时的分潮都会发生频率混叠被错误映射到低频段。潮汐分潮里常见的短周期成分是浅水分潮 M4周期 6.21 小时和 M6周期 4.14 小时1 小时采样完全没压力。但如果你用的是 6 小时间隔的验潮数据M6 就已经接近混叠区了这时候 t_tide 不会报错它只会默默把可疑分潮的误差放大输出里 SNR 会很低你要学会看这个信号。第二个边界是分潮之间的频率间隔。上面提过 M2 与 S2 需要约 14.77 天分离这是一个硬约束。程序包本身会按这个逻辑决定哪些分潮进拟合哪些分潮被合并或跳过。具体实现里它用 Rayleigh 准则判断当两个分潮的频率差大于1/TT 为总时长时才把它们当作独立分潮。如果数据不够长它会把次要分潮通过“分潮推断inference”方式处理而不是直接忽略。分潮推断是这个工具包里很实用的一个特性。比如你只有 30 天数据K1 和 P1 无法独立分离但你知道本海区 P1 的振幅大约是 K1 的 0.3 倍相位差也相对稳定就可以用inference, P1/K1让程序基于 K1 的解按比例推断 P1 的贡献而不是把 P1 直接扔掉。这个参数在短序列分析里几乎必用。% 若你确认本海区 P1/K1 振幅比约为 0.30相位差约为 0可显式指定推断关系 tidestruc t_tide(x_reg, start time, t_reg(1), dt, dt0, ... latitude, 30.5, inference, P1/K1, synthesis, 1);latitude传入观测站位纬度用于计算交点因子nodal factor和节点修正。你可能觉得这参数不重要但在高纬度海区全日分潮和半日分潮的节点因子随 18.61 年周期有明显摆动不传纬度会导致调和常数缺一块系统性误差。start time指定序列起点因为节点修正需要知道具体日期。synthesis设为 1 时程序会在输出结构体里附带一条拟合潮位曲线方便你直接和原序列对比。3. 跑通主程序预处理、核心调用与结果结构体3.1 预处理脚本清洗、均匀化与去平均在实际把序列交给 t_tide 之前我习惯先做一轮信号体检。所谓体检就三件事剔除尖峰、检查缺测比例、验证采样均匀性。尖峰的来源很多——设备打嗝、波浪瞬时越过传感器、生物附着后突然脱落。一个孤立尖峰对调和分析的影响比想象中大得多因为它包含全频段能量会让最小二乘拟合为了迁就它而扭曲附近分潮的振幅。我用的是滑动窗口加阈值法以 3 小时窗口计算局部中位数把偏离中位数超过 5 倍局部 MAD中位绝对偏差的点标记为尖峰然后剔除并线性插值补上。缺测比例的容忍度要分情况。如果缺测是零散的、单点缺失插值问题不大如果是连续缺测超过 3 小时而且总缺测比例超过 5%我会考虑分段分析或者干脆放弃这段数据。连续缺测在插值后看起来补上了但在频域上等于人为造了一段平滑过渡会压低高频分潮的振幅M4、M6 这类浅水分潮最吃亏。均匀化这一步在第 2.1 节的代码里已经做了这里再补一个细节插值前要先确认时间轴没有重复点。验潮站偶尔会因为时钟回拨产生两个完全相同的时间戳interp1遇到重复点会直接罢工所以要在插值前加一句% 去除时间重复点并保持原有顺序稳定 [~, idx_unique] unique(t_num, stable); t_unique t_num(idx_unique); x_unique raw_x(idx_unique);3.2 t_tide 调用方式与五个高频参数预处理完之后正式的调用代码如下% 核心调用返回 tidestruc 结构体 tidestruc t_tide(x_reg, ... start time, t_reg(1), ... % 序列起点用于节点因子 dt, dt0, ... % 采样间隔单位天 latitude, 30.5, ... % 站位纬度用于交点因子 inference, P1/K1, ... % 短序列下的分潮推断 synthesis, 1); % 输出合成潮位参数逐个说start time必须是 datenum 标量不要传字符串。传字符串虽然也能解析但多一道转换就多一处出错机会。dt的单位是天不是秒1 小时间隔就是1/24别填成 3600。latitude影响的是天文参数中的交点因子计算对中低纬度影响不算大但对高纬度海区比如 60°N 以上影响明显别漏传。inference的格式是P1/K1含义是用 K1 的结果推断 P1如果想同时推断多组可以写成元胞数组或者用逗号分隔的多组键值对。synthesis设为 1 后你可以在tidestruc.synthesis里拿到程序拟合的潮位曲线这个曲线是后续画残差图的基础。如果你用的是 Python 生态也有对应的移植版本可以用底层算法一致只是接口风格向 Python 侧靠拢。我大部分工作还是在 MATLAB 里做写 Python 版的人往往在配置文件和分潮列表上做了简化新手容易忽略latitude的作用导致结果和 MATLAB 版对不上。所以能用原版 t_tide 就用原版跨语言复核结果时再参考 Python 版。3.3 输出结构体怎么读振幅、迟角与信噪比调用完成后结果都存在tidestruc这个结构体里。最重要的字段是tidestruc.tidcon它是一个表格每一行对应一个分潮核心列包括频率、振幅、振幅误差、迟角、迟角误差和信噪比SNR。振幅的单位和你输入水位序列的单位一致输入米则振幅为米。有一点容易误解这里的振幅是分潮的半个潮差也就是在调和分析模型x(t) A*cos(ωt φ)中的系数 A不是该分潮贡献的峰峰值。如果你要描述“M2 分潮的半日潮差约 1.2 米”那指的是 2A这个坑几乎每个第一次用的人都会踩。迟角是相对于格林尼治天文时角的相位不是相对于本地时间的相位。所以不同经度站位之间的迟角不能直接比较要先按经度差把相位换算到同一参考。SNR 是振幅与振幅误差比值的平方是判断分潮是否可信的最直接指标。通常 SNR 大于 1 才认为该分潮被有效分离大于 2 才算比较可靠。如果某个分潮 SNR 小于 0.5基本可以判定它是噪声的产物不要拿来当真实潮流要素用。下面是一段绘制分潮振幅分布和合成潮位对比的常规代码% 读取 tidestruc 中的分潮表 tbl tidestruc.tidcon; % 只看 SNR 1.5 的分潮 idx_ok tbl.snr 1.5; sel tbl(idx_ok, :); % 画振幅误差棒图 figure; errorbar(sel.freq, sel.amp, sel.amp_err, o); xlabel(频率 (cycles/hour)); ylabel(振幅 (m)); title(主要分潮振幅分布); % 对比合成潮位与原始序列 figure; plot(t_reg, x_reg, k-); hold on; plot(t_reg, tidestruc.synthesis, r--); legend(观测, 调和拟合);tbl.snr这一列名在不同版本里可能略有差异如果字段名对不上用fieldnames(tidestruc)列一下所有字段名找对应关系。这个表就是后续所有分析的核心资产——有了分潮的调和常数你可以做水位预报、算设计潮位、反推潮汐类型也可以对比不同站位之间的潮汐差异。4. 常见问题与避坑时间格式、缺测与边界三条血泪记录4.1 时间格式搞错调和常数全发散现象t_tide 跑完输出的分潮振幅大得离谱M2 振幅算出来 8 米残差比原始序列还高一个量级SNR 全部小于 0.1。原因时间传错了。最常见的是把时间当成了从序列起点开始的秒数比如seconds(0:3600:36000)而 t_tide 期望的是绝对时间的 datenum。如果时间轴起点是 0它会把节点因子按 1800 年附近的日期去算交点因子完全偏离正常范围拟合当然发散。还有一种情况是把 datenum 传成了字符串数组或者把时间轴的单位搞混——datenum 的单位是天不是小时。解决在调用前强制校验。打印t_reg(1)和datestr(t_reg(1))确认起点日期合理打印dt0确认它约等于1/24而不是 1 或 3600。我自己养成一个习惯无论数据从哪来进 t_tide 之前必跑一行assert(abs(dt0 - 1/24) 1e-6, 时间间隔不是1小时)。这一行救过太多次了。4.2 缺测插值不当算出一整套假分潮现象分析一个缺测率只有 3% 的序列结果却很神秘——输出里多出若干个 SNR 挺高、但频率不在任何标准分潮表里的分潮而且残差里出现一种低频振荡。原因缺测段用线性插值补上后插值段在频谱上表现为一个低频平滑信号。这个平滑信号的能量会散布到邻近的潮频段如果恰好有一个待拟合分潮的相位和插值段的位置碰上了拟合就会把插值段的误差吸收进那个分潮。缺测越集中这种污染越严重。我遇到的一次案例里某观测站设备断电 5 小时插值后 M2 振幅被拉高了约 12%看起来很不合理最后查出来是插值段惹的祸。解决处理缺测前先看缺测分布。如果缺测是分散的单点线性插值没问题如果存在连续超过 4 个采样间隔的缺测段建议用潮汐信号特性做插值而不是纯数学插值。常见做法是先跑一次不含缺测段的初步调和分析得到分潮参数然后用这些调和常数预测缺测时段的水位并填回去再跑第二次正式分析。这种做法本质是用潮汐自身规律补数据比线性插值物理上更合理。4.3 序列截断边界引入 14.77 天伪周期现象序列本身是 30 天完整数据分析结果里 S2 分潮的 SNR 异常高但振幅明显比邻近站位小一圈。而残差序列里出现一个周期大约 14.77 天的波浪状起伏。原因大概率是你在数据两端做了截断处理或者从长序列里切了一段出来但切的时候没有保证首尾相位平滑。M2 与 S2 的差频周期恰好是 14.77 天当你用一段有限长序列做最小二乘拟合时频率分辨率是1/TM2 和 S2 两个频率在这个分辨率下可以被分离但会互相泄漏。如果数据段首尾恰好落在 M2 的高潮位附近泄漏就会集中在 S2 上导致 S2 振幅偏大或偏小。解决一是尽量让分析序列长度刚好是若干个完整天文日比如 30 天、60 天整。二是分析前用窗函数淡化边界。t_tide 本身不内置窗函数你可以在预处理阶段对序列乘一个 Hann 窗但要注意乘窗会改变振幅基准所以更实际的做法是做两次分析一次用原始序列一次用去首尾各 10% 的序列比较 M2 和 S2 的振幅变化。如果两次结果相差超过 3%说明边界泄漏严重需要回去查你的数据段是否有异常。4.4 深度不足时浅水分潮被当成噪声现象在近岸浅水区跑分析M4 和 MS4 分潮的 SNR 都很低只有 0.3 左右看起来像噪声但实测水位曲线在涨潮和落潮段明显不对称说明浅水分潮确实存在。原因浅水分潮振幅本来就小在强潮海区通常只占主分潮的 2%5%如果观测水位里有波浪噪声或者数据分辨率不足M4 就会被淹没。另一个原因是分潮列表不完整——如果在默认分潮表之外没有显式添加浅水分潮程序可能在拟合时把 M4 的能量分散到相邻频段。解决在调用 t_tide 时显式传shallow, M4 MS4 M6强制把这三个常见浅水分潮纳入拟合。加完之后再看 SNR通常能从 0.3 拉到 1 以上。这个过程也提醒我潮汐分析不是全自动的事你得知道自己所在的区域是半日潮型还是浅水分潮发育的海区然后主动告诉程序该关注谁。4.5 迟角符号与基准时区误读现象不同站位之间比较 K1 迟角时发现同一海区两个站位的迟角差接近 180°看起来像潮波反向传播但实测流场不支持这个结论。原因大多数资料里 K1 迟角写的是相对于格林尼治的相位但如果你参考的文献用的是本地标准时的相位两者会因为时区差异产生固定的相位偏移。对 K1 这种全日潮每差一个时区就偏 15 度。而且 t_tide 默认输出是相位用角度表示、包含符号习惯不同软件对于相位是超前还是滞后定义可能正好相反。跨软件比较时这个符号差会直接导致 180° 的误会。解决比较前先把所有数据统一换算。t_tide 输出角度是以余弦函数为基函数定义的换算公式这类细节在程序包自带的说明文档里有按它的约定来。我的做法是把不同站位的分析结果都换算到同一时区、同一基函数定义再画矢量图对比。如果发现符号反了通常不是数据问题而是软件约定问题直接乘一个负号就回来了。5. 进阶用法用调和常数做 30 天预报并验证残差5.1 t_predic 预报与置信区间拿到 tidestruc 之后的进阶操作是水位预报。t_tide 配套的t_predic函数可以基于分析得到的调和常数预测任意时间段的潮位曲线。预报长度超过原始序列长度后误差会逐渐增长但对 30 天量级的预报调和分析结果足够用了。% 预报未来 30 天间隔与原始数据一致 t_future (t_reg(end) dt0) : dt0 : (t_reg(end) 30); x_pred t_predic(t_future, tidestruc, synthesis, 1); % 用原始序列末段评估预报效果 idx_eval t_reg (t_reg(end) - 3); baseline x_reg(idx_eval) - tidestruc.synthesis(idx_eval); resid std(baseline); fprintf(历史残差标准差: %.3f m\n, resid);这里用历史残差的标准差作为未来预报精度的粗估计是工程上常用的做法。如果残差里还有明显周期性信号说明某个分潮没被充分提取需要回到第 4 章排查。t_predic的输出长度由你传的时间轴决定它不要求时间轴必须是等间隔但实际使用中保持等间隔能让你后面画图、算频谱都方便。预报开始时刻建议从整点开始避免跨过原始序列最后一段造成相位跳变。5.2 残差检验收敛没有看这两张图调和分析算完不能直接信残差是最后的照妖镜。第一张图残差时间序列观察有没有明显的正弦包络残留。如果有说明存在未纳入拟合的周期分潮常见的是遗漏了长周期分潮或浅水分潮。第二张图残差功率谱如果频谱在某个潮频段还有明显峰值说明对应分潮没被分离干净。判断收敛的量化标准也很简单残差的标准差应明显小于主要分潮振幅至少低于 M2 振幅的 10%。如果你的残差标准差接近 0.3 米而 M2 振幅只有 0.5 米那这个分析结果基本不可用不要拿去报设计潮位。做完这些验证我习惯把 tidestruc 存成二进制文件归档文件名里注明海区、时间范围和数据来源。从那以后我每次跑潮汐分析都强制走一遍流程先打印时间和 dt再看缺测分布分析后一定画残差频谱三个动作全部通过才敢用这套调和常数。这样做之后翻车的概率大幅下降希望帮到你。本文还有配套的精品资源点击获取