海洋潮汐调和分析实战：t_tide工具箱提取分潮常数与潮汐预报

简介海洋潮汐分析工具箱 t-tide v1.3beta 是一套面向海洋科学、港口航道及水利工程从业者的 MATLAB 潮汐分析工具集主要解决潮汐观测数据的调和分析、分潮提取、误差评估与未来潮位预测等问题。压缩包共 18 个文件包含 11 个 m 脚本、5 个 mat 数据文件和 2 个 dat 原始潮汐数据整体大小约 887KB。核心脚本 t_tide.m 基于谐波分析方法分解主要分潮t_predic.m 可依据模型生成预测序列t_errors.m 负责计算拟合误差与置信区间t_synth.m 支持合成潮汐过程以验证模型配套的 tide3.dat 与 t_equilib.dat 则提供历史观测和平衡潮汐参考数据。资源虽小巧但功能链完整适合具备一定 MATLAB 基础的科研人员与工程技术人员直接调用、修改或集成到业务系统中。该资源已有 178 人学习下载可作为潮汐数据预处理、特征分析和工程预报的实用工具箱。1. 海洋潮汐分析工具箱拿原始水位曲线换一套可信的分潮常数拿到某海湾一整年的验潮资料先别急着画水位曲线——用海洋潮汐分析工具箱 t_tide v1.3beta 把时间序列拆成几十个天文分潮才能把“这天的水位为什么高”从玄学变成可计算的迟角。这套工具箱解决的核心问题是从原始水位波动里分离出天文潮的确定性部分输出每个分潮的振幅和迟角并给出预报序列与拟合残差。它适合做港口工程、海岸带环境评价、验潮站数据质控和短期潮汐预报的人新手能靠它把一段杂乱的水位记录拆成一张分潮清单熟手则拿它做风暴潮余水位的预处理。如果手上的数据只有几天长度也能跑通只是分潮数量和置信区间都要降级处理这一点后面会反复强调。2. 主程序 t_tide 的调用与输出解读三个参数决定调和分析成败2.1 调和分析的最小二乘背景——为什么选 t_tide潮汐信号不是一堆随机波动而是若干频率已知的分潮叠加。每个分潮的角频率由天文参数展开确定水位序列可以写成若干个余弦函数的线性组合。调和分析要做的事情就是在一堆已知频率的余弦基底下用最小二乘估计每个分潮的振幅和初相。理论上只要频率已知剩下就是一个线性回归问题但实际操作里最卡人的是两个分潮频率太接近资料长度不够时它们会互相“抢能量”解出来的振幅和迟角让下游预报直接翻车。t_tide 的聪明之处在于把天文参数内置。分潮的振幅并不是恒定不变的月球轨道交点以 18.61 年为周期缓慢变化交点因子和节点订正值随之起伏t_tide 在每次运行时会根据起始时间重新计算这些因子不需要你手动查天文年表。同时它用 rayleigh 参数控制分潮分离的松紧程度用 inference 参数让弱分潮从强分潮那里“借”振幅和迟角这些都是短资料场景下的后悔药。选型理由很简单同样功能我过去也尝试自己写最小二乘求解器但潮汐频率表、节点修正、浅水分潮组合这三块工作量完全不在一个量级直接用现成工具箱再核实边界条件是性价比最高的路径。2.2 第一次完整调用输入、返回结构与最小可用例子调用 t_tide 的输入不复杂但有一个细节必须说清楚第二个参数是采样间隔单位是小时。整点每小时采样的数据就传 1半小时采样的数据传 0.5而不是传频率也不是传时间向量。% 读取两列数据[时间, 水位] data load(tide_obs.txt); tday data(:,1); % 单位天相对某个起点 z data(:,2); % 单位m dt_h 1; % 采样间隔1小时 starttime datenum(2023,1,1,0,0,0); % 序列第一点对应时刻 % 调用主程序 [tidestruc, xout] t_tide(z, dt_h, ... starttime, starttime, ... rayleigh, 0.95);这段代码里z 必须是单列水位序列t_tide 不接受带缺测的输入有 NaN 需要先补。starttime 必须和序列第一个点物理上是同一个时刻它决定了分潮初始相位的计算基准差一个小时后面所有分潮的迟角都会系统性漂移。xout 是原始水位减去潮汐拟合后的残差序列长度和 z 一样平时做余水位分析经常要用它。返回的 tidestruc 是一个结构体最常用的字段是 name、tidecon 和 fu。name 是本次求解包含的分潮名列表比如 M2、S2、K1、O1tidecon 是对应的参数矩阵fu 存的是交点因子。以下是几个关键参数的经验值参数作用备注interval采样间隔单位小时均匀采样缺测先插补starttime序列起点绝对时间用 datenum 精度到秒rayleigh分潮分离判据默认 1短资料可放宽到 0.9inference弱分潮推断如 P1/K1 表示 P1 由 K1 推断shallow是否加入浅水分潮河口等浅水区域建议开 1synthesis预报时是否包含所有分潮配合 t_predic 使用rayleigh 这个参数值得多说一句。它类似频率分辨的门槛数值越小允许被一起解出的分潮对就越多但代价是解的不稳定性上升。默认值 1 是资料长度比较充裕时的选择只有几天到十几天数据时我一般会手动放宽到 0.9同时配合 inference 减少独立未知数。不要为了把某个分潮强行解出来把值压到 0.8 以下那样看到的振幅基本都是虚高。2.3 tidecon 矩阵怎么读振幅、迟角与信噪比跑完 t_tide真正的产出在 tidestruc.tidecon 这个六列矩阵里。每一行对应 name 里的一个分潮列含义分别是频率、振幅、迟角、振幅中误差、迟角中误差、信噪比。振幅的单位和输入水位保持一致输入米就是米输入厘米就是厘米迟角单位是度。tide tidestruc.tidecon; % 用分潮名定位 M2避免行号记错 idx strcmp(tidestruc.name, M2); fprintf(M2: 振幅%.3f m, 迟角%.2f°, SNR%.1f\n, ... tide(idx, 2), tide(idx, 3), tide(idx, 6));第六列信噪比是我每次先看的值。信噪比等于振幅除以振幅中误差如果小于 2基本说明这个分潮没有被可靠分离振幅和迟角都不可信更别说拿去预报。常见的情况是短资料下 O1 和 K1 分不开两者振幅差不多但相位互相干扰SNR 一片稀烂。这时候不要硬取数值老老实实把资料凑够一个月以上或者用 inference 把弱分潮挂到强分潮上才是正确方向。还有一种我踩过的坑直接按前面输出的行号去手动取值后来换了数据长度后行号移动取到了隔壁分潮。从那以后我每次都用 strcmp 按分潮名定位绝不用行号索引这个习惯也建议你直接抄走。3. 从调和常数到预报序列t_predic 的相位起点与置信区间3.1 交点因子与节点修正——为什么不能拿常数直接外推很多刚接触调和分析的人犯一个共同的错把 tidestruc 里输出的振幅和迟角当成固定值直接代入余弦函数去推算未来几个月的水位。这个做法在十天半个月内可能看起来还行一旦拉长到半年以上误差会越来越大。原因是分潮振幅其实随时间缓慢变化最典型的是月球交点 18.61 年周期造成的幅度调制另外一个分潮的迟角也会因为天文参数的变化产生偏移。t_tide 在分析时段内部已经把这一时刻对应的交点因子 f 和节点订正值 u 算进了解的结果也就是说你拿到的分潮常数是“这个观测时段内有效”的等效常数。要做预报必须让预报程序重新计算每一天的交点因子和节点订正而不是机械地复制分析结果。这就是为什么 t_tide 后面必须接 t_predic而不是自己在脚本里写一行 cos 循环。t_predic 做的事情本质上是把同一套天文参数框架重新演算一遍只是把时间轴从观测段平移到了预报段。3.2 t_predic 的调用与时间换算t_predic 的第一个参数是相对小时数是从分析段起点开始往后的小时偏移量。这个参数不是绝对时间不能直接把 datenum 传进去这是新手在时间轴上报错最多的地方。假设分析起点是 2023 年 1 月 1 日 0 点要预报起点后 0 到 30 天、步长半小时的潮位正确写法是这样dt_pred 0.5; % 预报步长0.5小时 hours (0:dt_pred:30*24); % 从起点之后30天单位小时 [pred, ~] t_predic(hours, tidestruc, synthesis, 1); % 自己换算绝对时间用于画图 abs_time starttime hours / 24;synthesis 参数设为 1 表示把参与分析的所有分潮、包括推断出来的分潮都合成进预报值。如果不加这个参数或者设为 0只输出一部分分潮的叠加结果用途不大。返回的第二个参数在很多版本里是相对时间或者别的参考系我一般直接无视它自己用 starttime 加上偏移量换算绝对时间这样可以完全规避版本间返回值含义不一致的问题。预报步长的选择也有讲究。半小时步长足够给绝大多数港口工程使用但如果你需要精确到高低潮时刻建议把 dt_pred 再加密到 0.25 小时也就是 15 分钟。太粗的步长会让高潮峰顶被“压扁”和实测资料比较时误以为预报偏低。3.3 长期预报与置信区间t_predic 本身不直接返回误差带但 tidecon 里的振幅中误差给了我们一个快速估算置信区间的途径。做法是对振幅做上下偏移生成两组预报序列包出来的区间就是近似的预报可信带amp tidestruc.tidecon(:, 2); amp_se tidestruc.tidecon(:, 4); % 造一份上下1.96倍中误差的副本其余参数不动 td_low tidestruc; td_low.tidecon(:, 2) amp - 1.96 * amp_se; pred_low t_predic(hours, td_low, synthesis, 1); td_high tidestruc; td_high.tidecon(:, 2) amp 1.96 * amp_se; pred_high t_predic(hours, td_high, synthesis, 1);这个办法只考虑了振幅误差没有把迟角误差和分潮间协方差算进去得到的区间会比严格置信区间窄一些。做工程初判够用写进正式报告前需要从最小二乘协方差矩阵严格推导。更实际的经验是如果这套分析用的观测资料不足一个月我基本不会拿它做超过三个月的预报因为置信带很快就会宽到失去意义只有凑够了完整的一个月以上资料置信带才会明显收窄这时候预报结果才敢交给下游做航道水深复核。4. 避坑指南从矩阵求逆失败到预报相位漂移的排查记录4.1 报错 Matrix is singular or badly scaled现象t_tide 运行到一半直接报矩阵奇异或者给出一个极其夸张的振幅比如 M2 振幅算出 5 米。原因输入水位序列里有 NaN 缺测或者整段序列全是同一个常数最小二乘正规方程无法求逆。解决先执行 isnan 检查所有缺失点用 fillmissing 做线性插补。注意不要用 0 填空那等同于引入一段人为的低频信号。插补后如果序列长度还不够一次完整潮周期那就换数据源强行分析没有意义。4.2 解出的 M2 振幅比邻近站点大一圈现象某段短资料算出的 M2 振幅明显高于当地长期潮汐表而且 SNR 还不低。原因资料长度太短M2 和它旁边的浅水分潮或 N2 没有彻底分开部分能量被错误分配给了 M2。表面看拟合残差不大实际是分潮之间互相“抢”振幅。解决先看 tidecon 里 N2、M4 的振幅和迟角是否异常再用 rayleigh0.9 重跑一遍对比如果依然异常用 inference 把 N2 挂到 M2 上减少独立未知数。千万不要只盯着 M2 一个分潮调和分析是整体解单个分潮好看不代表全部可靠。4.3 预报与实测高低潮时刻系统性偏移几小时现象t_predic 报出的高潮时刻比实测稳定提前 3 小时振幅看起来又差不多。原因starttime 和数据实际起点错位了。最常见的是数据文件里记录的是当地时而 starttime 用了 UTC或者反过来另一个常见错误是数据第一行对应的物理时刻不是整点但 starttime 写成了整点。解决把数据文件里的第一行时间和 starttime 严格对齐精度到秒。处理陌生站点数据时我在跑正式分析前会先画一张观测和预报的对比图如果发现相位存在固定平移第一件事就是回去核对时间基准而不是怀疑分潮常数本身。4.4 K1、P1 振幅一正一负互相“吸”现象同一个站用不同长度资料重复分析K1 和 P1 的振幅每次差异都很大符号都有翻车的可能。原因K1 和 P1 的频率差非常小需要约 14 天以上的资料才能分开。资料短于这个长度时两者几乎完全混淆最小二乘解极易振荡。解决短资料场景下显式添加 inference 参数把 P1 的振幅和迟角由 K1 推断例如[tidestruc, xout] t_tide(z, dt_h, ... starttime, starttime, ... inference, P1/K1);这样 P1 就从独立求解变成了“K1 的附属”整体解的稳定性会显著提升。记住一句话宁可少解几个分潮也不要解出一堆符号乱跳的虚假常数。4.5 残差 xout 里明显残留正弦波动现象t_tide 跑完后xout 仍然呈现明显的周期性起伏而不是噪声样子的短周期波动。原因某个重要分潮没有进入求解列表。可能是 shallow 参数没开漏掉了浅水分潮也可能是资料长度不够频率接近的分潮被 rayleigh 判据直接剔除了。解决打开 shallow1 试一次对比 xout 的标准差是否显著下降如果依然有长周期残留考虑把资料段做带通滤波后再分析。我一般拿 xout 的标准差和原始水位信号的标准差做比值如果残留能量超过总能量的 10%这个调和分析结果我不会拿去出预报产品。5. 把输入数据做对文件清单、预处理流程与 M2 分潮提取实例5.1 工具箱文件清单与职责分工拿到 t_tide v1.3beta 的压缩包解压后不是把所有文件一股脑丢进 work 目录就完事先认清楚每个文件干什么文件名职责通常是否修改t_tide.m主调和分析程序否t_predic.m由分析结果生成预报序列否t_vuf.m计算交点因子和节点订正值否t_astro.m计算天文参数否readme.txt使用说明与版本差异说明否这个工具箱的依赖关系是 t_tide 和 t_predic 调用 t_vuft_vuf 再调用 t_astro所以四个 .m 文件必须放在同一个可访问目录里。安装方式就是在 MATLAB 里把解压目录加到搜索路径addpath(你的路径/t_tide_v1.3beta);有些版本还带一个名为 uttide 的子目录里面是程序内部使用的辅助函数不要单独把它移到别处保持目录结构完整即可。如果已经有更老版本的 t_tide 在环境中建议先清理掉旧路径否则两个版本同名函数会互相覆盖这种冲突在 MATLAB 里很难排查。5.2 输入数据格式与预处理建议数据格式上最常见的可复现方案是 ASCII 文本文件列顺序为年、月、日、时、分、秒、水位。水位的单位我在文件头注释里写清楚是米还是厘米防止几个月后自己回读数据时犯迷糊。d load(obs.dat); % 假设列顺序YYYY MM DD HH MM SS 水位(m) tnum datenum(d(:,1), d(:,2), d(:,3), d(:,4), d(:,5), d(:,6)); z d(:,7); % 1) 缺测线性插补 z fillmissing(z, linear); % 2) 尖峰压制与5点中值滤波结果比较偏差超过阈值的点替换 z_med medfilt1(z, 5, truncate); bad abs(z - z_med) 0.5; % 阈值按当地潮差调整 z(bad) z_med(bad);中值滤波的窗口大小和阈值不是通用参数。某站半日潮潮差 3 米时我用 0.5 米做尖峰判定是合理的如果一个站潮差只有 0.3 米同样阈值会把正常潮汐变化也压掉必须把阈值降到 0.1 米左右。这一点靠脚本自适应比较难我的习惯是先用一年数据算出振幅标准误以 5 倍标准误作为尖峰判定阈值。5.3 从原始水位到 M2 分潮常数的完整流程数据预处理干净后完整的 M2 分潮提取流程可以一把跑完。这里我习惯顺手把结果写到 CSV方便后续和官方潮汐表对比% 主分析 [tidestruc, xout] t_tide(z, 1, ... starttime, tnum(1), ... rayleigh, 0.95, ... shallow, 1); % 按分潮名定位输出M2相关常数 names tidestruc.name; tidecon tidestruc.tidecon; fid fopen(m2_result.csv, w); fprintf(fid, name,amp_m,phase_deg,amp_se_m,phase_se_deg,snr\n); for k 1:numel(names) if contains(names{k}, M2) fprintf(fid, %s,%.4f,%.2f,%.4f,%.2f,%.1f\n, ... names{k}, tidecon(k,2), tidecon(k,3), ... tidecon(k,4), tidecon(k,5), tidecon(k,6)); end end fclose(fid); % 画出观测与预报对比快速目检 hours_obs (0:numel(z)-1); pred_full t_predic(hours_obs, tidestruc, synthesis, 1); plot(tnum, z, k); hold on; plot(tnum, pred_full, r--);用 contains 匹配 M2 可以同时把 M2、M2 的浅水分潮同名变体尽量覆盖到但要注意可能匹配到 M2 以外的含 M2 字样的分潮名。保险起见输出前看一眼 names 列表。绘图对比这一步不能省它能把相位错位、振幅异常这类问题暴露得比任何统计指标都直观。5.4 与官方潮汐表的比对验证拿到 M2 振幅和迟角后要和附近长期验潮站或官方潮汐表的同分潮常数做交叉验证。工程上我一般按下面的口径判断振幅相对偏差在 5% 以内算正常迟角差在 10 度以内算吻合。迟角差经常包含一个固定常数那不是程序问题是基准问题——官方表用的时间标准可能和你传入的 starttime 基准不一致需要先把两边统一到同一时间体系再比较。如果振幅差在 10% 到 20% 之间先检查是不是资料长度不足导致的混淆再检查浅水分潮是否参与求解。如果差超过 20%多半是数据源头有问题比如水位计零点漂移、记录时间跳变这种时候不要急着调参数先回去看原始数据的连续性和水准联测记录。6. 进阶用法余水位分离、残差诊断与长期资料的业务化处理6.1 从观测减预报分离余水位潮汐分析做完以后把整段观测序列减去同长度的潮汐预报序列得到的残差就是余水位。余水位代表气象扰动、径流顶托、地形非线性等非天文因素造成的水位偏离是风暴潮研究和航道回淤分析的重要输入。pred_full t_predic((0:numel(z)-1) * dt_h, tidestruc, synthesis, 1); resid z(:) - pred_full(:); % 检查残差是否还有明显的周期成分 plot(tnum, resid);如果你发现 resid 里依然存在明显的半日周期波动说明调和分析没有把主天文潮完全分离干净问题通常出在分潮数量不足或浅水分潮没开而不是余水位本身如此。6.2 残差诊断的量化指标判断一次调和分析是否到了可以交付的状态我会计算残差的均方根误差和纳什系数rmse sqrt(mean(resid.^2)); nse 1 - sum(resid.^2) / sum((z - mean(z)).^2);RMSE 要结合站点潮差理解一般控制在潮差的 10% 以内算合格NSE 大于 0.9 说明潮汐预报抓住了主要方差。这两个指标每跑完一次分析都输出一次如果 NSE 低于 0.8我不会继续往下游分析走先回头补数据或调整分潮配置。6.3 长期资料的分段处理与分潮推断策略遇到多年连续资料时我的处理习惯是分段分析而不是整段一次塞进去。每一年单独跑一次 t_tide保留每年的 M2、K1、O1 等主分潮常数然后看年际变化。分段的好处是能直观看到交点因子变化对有效振幅的影响也能及时发现某一年水位计故障造成的异常段。对资料长度不足的分潮统一用 inference 策略例如把 P1 挂到 K1、把 N2 挂到 M2 上避免每年解出来的弱分潮常数互相乱跳。前两年我在某海湾处理逐时水位资料时急着出结果用固定振幅外推了三个月的预报最后和实测差了一截。从那以后每个新站点的数据我都会强制走一遍同样的流程先看 SNR再画观测减预报残差曲线确认没有周期残留后才把调和常数交出去。希望帮到你。本文还有配套的精品资源点击获取