简介这是一款基于t-tide v1.3beta的海洋潮汐分析工具箱面向海洋科学研究者、港口工程师及物理海洋学相关专业学生用于潮汐调和分析、分潮参数提取与未来潮位预测。压缩包共18个文件以11个Matlab脚本为核心覆盖潮汐成分分解与重构、合成时间序列、误差评估、天文常数读取、潮汐预测以及流速向量分析等完整流程另含5个mat数据文件为分潮参数与示例数据2个dat文件记录历史观测和平衡潮汐数据整体大小约887KB。资源已有178人学习适合具备Matlab基础、需要处理长周期潮位观测数据或开展潮汐动力学研究的进阶用户。借助该工具箱可快速求解主要分潮的振幅与相位生成预测曲线检验分析误差并可根据研究区域定制模型相比从零编写谐波分析代码这套工具能显著提升潮汐数据处理与研究效率。1. 海洋潮汐分析工具箱先别急着解压先搞清楚它解决什么问题接到一个港湾工程的水位分析任务手头是某临时验潮站一年的逐时水位目标是给出主要分潮的调和常数、预报未来一个月的逐时潮位还要估算误差。如果手动按潮汐表一个个分潮去掰周期和相位牵扯到几十个正弦分量根本算不完。海洋潮汐分析工具箱 t-tide v1.3beta 就是干这个的它在 MATLAB 环境下把一整段水位序列拆解成多个天文分潮通过最小二乘拟合得到振幅和相位再同步给出误差条和预报潮位。适合港航设计、岸线工程、环境监测需要长期复核潮位规律的技术人员也适合刚接触潮汐分析但不想从零写谐波算法的研究生。它不是一个“点一下出图”的傻瓜工具而是一套需要理解输入输出和参数控制的底层函数包这篇笔记打算把它讲透。2. 潮汐调和分析的底子为什么这类工具箱能“预报”潮汐2.1 天文潮的本质是一组合成正弦波潮汐产生的直接驱动力是月球和太阳的引潮力对应的海平面响应并不是一个正弦波而是若干频率不同、振幅不同、初相位不同的正弦分潮叠加。每个分潮有固定的天文角速度也都有自己的工程名称。最常见的是 M2主太阴半日分潮、S2主太阳半日分潮、K1日月合成日分潮、O1主太阳日分潮。在近岸河口还会出现 M4、MS4 这类浅水分潮它们不是直接来自引潮力而是主要分潮在浅水中非线性变形产生的倍频分量。观测上任意时刻的潮位可以近似写成z(t) Σ Ai·cos(σi·t φi) R(t)其中 σi 是第 i 个分潮的角速度Ai 是振幅φi 是初相位R(t) 是残差项。工程上最关心的是 Ai 和 φi因为它们合起来决定了当地的潮汐类型和预报精度。潮汐类型的判断也依赖这个结果M2 振幅占比很大的区域往往属于规则半日潮O1 和 K1 为主的区域则偏向日潮。教科书经常把这个问题讲成“从一片混沌的曲线里分离规律成分”听起来很玄但到了数值计算层面就是一个求解线性模型参数的问题。每个分潮的角速度可以由天文参数精确算出来并不需要从数据里去猜。比如 M2 的角速度约是 1.40519 度/小时S2 约是 1.45445 度/小时两者非常接近这就是为什么短周期数据很难把这两个分潮分开的原因。频率相近的分潮需要足够长的观测时间才能被区分这个时间边界在调和分析里有一个专门的名字叫 Rayleigh 判据。2.2 最小二乘拟合t_tide v1.3beta 的求解核心任何潮汐调和分析工具都不会逐个分潮去猜而是把所有选定分潮的角速度代入同一个线性模型中一次性解出全部振幅和相位。把方程写成矩阵形式目标是用观测值 y 找到未知向量 x使 A·x 尽可能等于 y。% 以下代码展示最小二乘求解思路不是 t_tide 的完整调用 N length(y); A zeros(N, 2 * num_components); for k 1:num_components A(:, 2*k-1) cos(sigma(k) * t); A(:, 2*k) sin(sigma(k) * t); end x A \ y; % 最小二乘解 Amp hypot(x(1:2:end), x(2:2:end)); Pha atan2(x(2:2:end), x(1:2:end)) * 180/pi;这段代码的逻辑是把每个分潮拆成一个余弦项加一个正弦项两个待定系数合成一个振幅和一个相位把所有分潮的基函数水平堆在一起形成设计矩阵 A然后用 MATLAB 的反斜杠做最小二乘求解。这里的t是相对某个起算时刻的时间向量单位可以是小时或秒但必须和sigma的单位一致。num_components是参与拟合的分潮数如果数据长度小于分潮数两倍左右方程就开始变得病态解出的振幅和相位不再可靠。t_tide v1.3beta 内部就是在最小二乘基础上加入了很多工程化细节包括逐项输出误差估计、根据数据长度和分潮角速度自动判断可分辨性、支持节点因子修正和倾斜项。所以拿到结果后不需要手动算置信区间它直接给出一列误差值。这也是为什么在同一个场景下用它比从零写脚本要稳妥得多。你可能会问既然最小二乘原理这么简单为什么还有人不直接写脚本原因是写一个能跑的概率高但要写出面对缺测、长短序列、约束条件、误差传播与异常谱峰都不翻车的脚本就远不是一两百行能解决的事。2.3 为什么基于 MATLAB 的经典工具箱还不过时这类工具箱虽然在工程圈流传已久但在实际项目里仍是默认选择。原因有几个一是 MATLAB 在时间序列处理、绘图和批量调用上足够顺手不需要另搭编译环境二是潮汐分析本质上是数值线性代数MATLAB 的矩阵运算和线性方程求解成熟稳定三是这个版本的函数接口长期保持稳定很多业务代码和内部规范都是围绕它写的。相比之下临时写的自定义脚本看着只处理一个站点时很简单但遇到数据空缺、异常值或分潮混叠时就很难收场。选择一个久经考验的工具箱相当于把校验、误差估计和分潮挑选的工作交给一套成熟逻辑自己只需要关注数据质量和后处理。判断“能不能用”的简单标准是把函数目录放进 MATLAB 路径先跑通自带示例数据再跑自己的真实数据看残余序列是否稳定、误差条是否合理。2.4 观测序列里不止有天文潮残差也是一部分结论很多第一次上手的人容易把“调和分析”理解成“从实测水位里把潮位分离出来”这是不够准确的。天文潮的确定性成分只占总水位变化的一部分剩下的残差里还有气象增水、径流变化、波浪非线性、仪器漂移等因素。t_tide 的作用不是把所有波动都解释掉而是通过一组分潮把确定性最强的天文潮提取出来让残差尽可能不包含周期性成分。这一点在工程上的意义很大。例如在做风暴潮增水分析时通常先用 t_tide 或类似工具算出天文潮再从实测水位里减去天文潮得到的“增水序列”才更接近真实的气象扰动水位。如果分潮选取不完整残差里可能残留明显的周期信号后续增水分析就会被污染。所以要养成一个习惯每次跑完调和分析先看残差序列里是否还有周期性波形而不是只看拟合阶段的贴合度。3. 从压缩包到第一次跑通安装、目录结构与一次完整流程3.1 解压后的目录结构与路径设置拿到压缩包后先不要急着双击运行某个脚本应该先把文件夹完整解压然后看一眼里面有哪些函数文件。常见做法是目录下至少包含t_tide.m主分析函数、t_predic.m预报函数以及若干用于读取、画图或误差计算的子函数。如果压缩包还带一个样例数据子目录那是最理想的因为可以先拿示例数据跑一遍验证环境避免一开始就把自己手头数据的误差带进去。在 MATLAB 中设置路径有两种方式。如果只准备在当前会话使用直接addpath足够如果希望以后打开 MATLAB 自动加载则需要同时执行savepath。需要注意路径含有中文字符或 MATLAB 不支持的符号时函数可能找不到最好把工具箱放到一个纯英文路径下。cd /data/tools/ unzip t-tide-v1.3beta.zip ls -la t-tide-v1.3beta/% 把工具箱加入当前MATLAB会话并保存路径配置 addpath(/data/tools/t-tide-v1.3beta); savepath;代码逻辑很简单先在系统命令行解压再用addpath将解压后的文件目录加进 MATLAB 搜索路径。savepath会把当前路径配置写入配置文件这样下次启动时依然有效。最容易被忽略的是路径分隔符在 Windows 平台上的写法写addpath(C:\tools\t-tide-v1.3beta)时字符串里的反斜杠需要转义成双反斜杠否则 MATLAB 会把它当成非法字符。3.2 准备观测数据从验潮站原始记录到可分析序列t_tide 的主函数要求输入一个单列等间距水位序列并附带起始时刻和采样间隔。如果你的原始数据是逐分钟记录不要直接丢进去因为奇异点太多。常见做法是把原始数据处理成整点小时级序列去掉明显缺测再用插值补齐个别空洞。数据列的顺序、时间格式、水位单位都必须提前核对好。raw readmatrix(sta_level.csv); % 假定的两列CSV raw(raw(:,2) -9999, 2) NaN; % 剔除缺测标志值 tnum raw(:,1); wl raw(:,2); dt diff(tnum); if max(abs(dt - median(dt))) 1e-5 warning(时间间隔不均匀建议先重采样); end wl fillmissing(wl, linear); % 缺测较少时线性插值这段代码先读取一个两列文件第一列是 MATLAB 的datenum时间戳第二列是水位。数据里常约定某个特异值表示缺测比如 -9999因此先把它转成NaN。diff用于检查相邻时间差是否足够均匀如果某一段缺测太久直接线性插值会引入虚假低频能量这种情况其实应该把那段数据切掉或按更精细的办法处理。补充一点如果原始数据是字符串格式的年月日必须先用datenum(text, yyyy-mm-dd HH:MM)转成数值时间绝不能把字符串直接传进去。这个预处理环节看似简单却是后面所有结果是否可信的地基。潮汐调和分析对时间基准非常敏感时间起点差一个小时相位结果就会整体偏一个角度水位单位从米变厘米振幅就差两个数量级。所以每次换新站点数据我都会强制把第一列、第二列的采样间隔、数据起止时间打印出来确认没有量纲混淆。3.3 第一次调用 t_tide最小可运行示例数据准备好之后直接运行主函数。最常见的最小调用如下[tcon, pred] t_tide(wl, start, tnum(1), interval, 1); figure plot(tnum, wl, Color, [0.7 0.7 0.7]); hold on plot(tnum, pred, r-, LineWidth, 1.2) legend({实测, 合成潮位}, Location, best)代码逻辑wl必须是列向量t_tide会默认选择一组适合数据长度的主分潮用最小二乘拟合出调和常数。输出tcon是结构体里面有分潮名称、频率、振幅、相位和误差pred则是和输入序列等长的合成潮位序列。画图时可以看到如果数据干净预报潮位和实测值在大部分时段能够贴合偶尔会有风暴增水造成的偏移。参数说明start必须是准确的时间起点interval表示观测数据的时间间隔单位是小时。如果你的数据是 6 分钟间隔就不能写interval1而应写成0.1十分钟一个点则写1/6。第一次跑通时不要急着加一堆优化参数先跑默认值通过残差判断数据有没有明显问题。默认输出里通常会有一些警告比如某些分潮因数据长度不足没有被解析这些是正常提示具体影响放到后面调优时处理。3.4 用自带样例或模拟数据验证环境完整性如果压缩包里带示例数据建议先跑通。没有样例数据时也可以自己生成一段合成潮位测试工具箱是否工作正常。这个步骤不复杂但对后续排查非常有价值。你可以用两个已知分潮合成一段一小时间隔的序列再调用 t_tide 反解看看振幅相位能不能还原。% 构造24小时模拟序列包含M2和K1 fs 1; % 每小时一个点 t (0:24*7-1)/24; % 一周时间单位天 wl_sim 1.2 * cos(2*pi/12.42*t) 0.5 * cos(2*pi/23.93*t 0.3); [tcon_sim, pred_sim] t_tide(wl_sim, start, datenum(2024,1,1), interval, 1); disp(tcon_sim.name) disp(tcon_sim.amp)这段代码生成了一周逐时水位由 M2 和 K1 两个分潮叠加而成。通过 t_tide 反解后应该能在tcon_sim.amp里看到接近 1.2 和 0.5 的两个振幅。如果这个自检都通不过问题大概率出在路径设置或函数调用版本上而不是在手头数据上。这个习惯我一直保留遇到异常结果会先用合成数据做对照实验快速定位问题出在工具箱还是数据本身。4. 参数调优与输出解读让预报结果更贴近工程实际4.1 分潮选取rayleigh 准则、自定义分潮与浅水分潮默认情况下t_tide 会根据数据长度自动选择可分辨的分潮。但工程项目里往往有明确要求比如只需要给出主要半日和全日分潮或者需要关注某几个浅水分潮。这时可以通过constituents参数手工指定。建议使用前先看一组默认结果再用 Rayleigh 准则检查分潮是否分得开。my_con {M2, S2, N2, K1, O1, M4, MS4}; [tcon, pred] t_tide(wl, start, tnum(1), interval, 1, ... constituents, my_con, rayleigh, 0.9);这段代码指定了七个分潮其中 M4 和 MS4 是浅水分潮。rayleigh参数控制两个分潮角速度差与数据长度的匹配程度常用值在 0.9 到 1.0 之间。如果序列只有 30 天M2 和 S2 之间的频率差很小虽然不一定能彻底分开但用rayleigh小于 1 仍能拟合出结果代价是误差条变大。对于要求严格的潮位预报最好使用 1.0 以上并准备至少 60 天以上连续数据。为什么需要手工指定而不是全自动因为某些浅水分潮在河口和近岸非常显著默认分组可能忽略它们反过来如果数据长度不够而强行塞进太多分潮最小二乘方程接近奇异解出来的振幅可能双双异常。这个区域是踩坑高频区也是最容易出现“两个分潮互相抢振幅”的地方。建议每次调组分时保留一份tcon副本记录组分和时间范围方便复盘对比。4.2 读懂 tcon 结构体从振幅、相位到误差条第一次运行后很多人习惯直接看pred却不看tcon这是不完整的。tcon里的每个字段都对应一个分潮的估计结果。在工程报告里最常提取的是名称、振幅、迟角、标准误差。tcon.name tcon.freq tcon.amp tcon.pha tcon.amp_err tcon.pha_err代码说明name是字符串元胞数组amp和pha分别是振幅和相位amp_err和pha_err是标准误差。如果amp_err大于对应振幅的 20%说明这个分潮的可信度很低要么数据长度不够要么附近有更强的分潮在干扰。把这一整组结果输出成表格可以进一步做基准面对比或不同测站间的空间插值。result_table table(tcon.name, tcon.amp, tcon.pha, ... tcon.amp_err, tcon.pha_err, ... VariableNames, {name,amp,phase,amp_err,phase_err}); writetable(result_table, harmonics_result.csv);这里的phase单位是度参考时间由start决定。如果你换了起点时间相位值会系统性变化所以在跨测站横向对比时要确认大家使用同一个时间基准。不少新手直接拿不同起点的相位数字做对比结果差了十几度还以为站点存在显著差异实际只是时间基准没对齐。4.3 用 t_predic 重建任意时段的潮位拿到tcon之后可以用它做两件事第一是对已有观测区间做合成对比第二是对没有观测的区间做预报。t_predic的输入是目标时间向量与调和常数结构体输出是对应潮位序列。t_new datenum(2025,1,1):1/24:datenum(2025,2,1); pred_new t_predic(t_new, tcon, synthesis, 1); plot(t_new, pred_new, b-);代码逻辑先构造富 2025 年 1 月 1 日到 2025 年 2 月 1 日的逐小时时间向量1/24表示一小时一个点。接着调用t_predic传入tcon结构体和目标时间。synthesis参数控制是否进行节点修正一般填 1 代表使用完整的合成因子。输出序列长度与时间向量一致逐点就是预报潮位。这个函数在业务里的常见用途是制作潮位预报表或者为工程施工窗口寻找低潮时段。但要牢牢记住它只预报天文潮部分台风、寒潮、径流影响不会被包含。当实测序列里出现过强增水时tcon的振幅和相位已经受了影响预报偏差也会比较大。解决办法是把极端事件剔除出拟合窗口或者用长序列稀释事件影响。4.4 多基准面对比与节点修正工程上经常会遇到“平均海平面”“理论深度基准面”“当地潮高基准”等不同基准面t_tide 输出的振幅和相位是相对输入水位序列的输入序列本身是什么基准输出的调和常数就躺在什么基准上。跨项目对比时必须先把输入序列统一到同一高程基准。% 若实测水位是相对当地潮高基准想转换到平均海平面 datum_offset 0.35; % 当地验潮零点与平均海平面差单位米 wl_corrected wl - datum_offset; [tcon_corrected, pred_corrected] t_tide(wl_corrected, ... start, tnum(1), interval, 1, ... constituents, {M2,S2,K1,O1});这里引入一个常数偏移量datum_offset把整段水位平移后重新拟合。由于调和常数只描述周期成分基线平移一般不会改变振幅但会影响后续水位预报的绝对高程。实际工程中这个环节很容易被漏掉最后做频率分析时发现所有预报潮位都比实测低 30 厘米返工成本很高。我的建议是在数据导入阶段就确定好基准并在脚本注释里写清楚偏移量来源不要等出图以后再去猜。5. 避坑清单潮汐分析里最容易翻车的五个现场在拆项目过程中有几个坑反复出现下面按“现象—原因—解决”的方式记录方便你对照排查。5.1 现象结果振幅比周围测站明显偏大或偏小原因输入数据单位不统一。t_tide的振幅单位严格跟随输入序列单位。如果你传入的是“厘米”单位但后续比较时用了“米”振幅自然会差 100 倍另一种情况是inorm或synthesis参数被修改过导致输出不是原始观测幅度。解决先检查原始数据列到底是不是米再检查调用时是否误写了inorm或synthesis。如果是厘米先除以 100。如果脚本里加了某个参数但不确定含义就把该参数删掉用默认值跑一遍对比。还可以打印前三个主要分潮振幅与区域内的公开参考值做量级核对M2 振幅一般能到几十厘米以上不会出现 0.0002 这种数字。若不正常优先怀疑量纲与参数而不是分析算法本身。5.2 现象程序报错提示输入数据含有时间间隔不均匀或 NaN原因预处理把大量缺测用线性插值补齐了但缺测段太长时插值序列并不均匀t_tide 内部检测到时间轴断裂后会直接拒绝执行。另一个常见问题是原始数据里已有重复时间戳或倒序时间戳。解决不要对一大段缺测做线性插值应该把缺测段切割开或舍弃。如果缺测零星分散可以先检查最大连续缺测时长再决定采用线性插值还是删除缺失段。时间检查这一步可以用相邻diff来定位异常点并把异常区间打印出来看是否属于仪器停机。遇到倒序时先用sortrows按时间升序重排再进入下一步。[~, sort_idx] sort(tnum); tnum tnum(sort_idx); wl wl(sort_idx); dt_hours diff(tnum) * 24; bad_idx find(dt_hours 1.05 | dt_hours 0.95); if ~isempty(bad_idx) fprintf(发现 %d 个非整点间隔最大位置在 %d\n, length(bad_idx), bad_idx(1)); end这段代码先排序再计算相邻间隔小时数筛出明显不是 1 小时的间隔。如果数目很少可以局部插值如果大面积不是整点就需要回到原始数据确认采样策略。t_tide 对等间隔要求比较严前提工作做扎实能省掉很多后期麻烦。5.3 现象某个主要分潮的误差条巨大且预报潮位和实测总是差 0.2 米以上原因分潮选取过多或过少。分潮过多时最小二乘矩阵接近病态部分分潮振幅互相“抢能量”分潮过少时未被解释的周期成分会混入残差拟合潮位出现系统低估。两者都会导致误差条变大。解决用默认分潮先跑一遍把误差超过振幅五分之一的那些分潮从列表里剔除再进行下一次拟合。对比两次结果的pred与实测之间的均方根误差如果误差没有明显下降说明继续增加分潮意义不大应该把注意力转到残差序列上。keep_idx strcmp(tcon.name, my_con); % 或使用 all(name, ...) valid tcon.amp_err ./ tcon.amp 0.2; disp(tcon.name(~valid)); % 显示误差偏大的分潮这段逻辑只是展示一种筛选方式用小误差阈值找出可信度低的分潮。实际项目中可能还需要结合分潮名和物理意义来判断不能只按比例盲目剔除因为某些浅水分潮本来振幅就小误差比例偏高也是正常的。5.4 现象相位结果从某一天开始突然偏移 180 度原因start参数设置错误。t_tide的相位参考点取决于start如果把它错给了另一个站点时段或错误年份而角速度没有变相位就会整体旋转。另外一个常见原因是把datetime类型直接当作start没有转成datenum。解决重新确认传入的tnum(1)是实际观测开始的绝对时间且类型是数值。用datenum(dt)显式转换。这个错很阴险因为振幅完全不受影响只有相位偏移单看振幅曲线很难发现。每次换数据时都检查起始时间打印值可以提前拦截这一类问题。fprintf(数据起始时间: %s\n, datestr(tnum(1), yyyy-mm-dd HH:MM));这段代码只是把起始时间打印出来。看起来不起眼但能帮你快速发现“是不是已经提前被插值到了某个新日期”或“时间格式被 MATLAB 自动转成了天数”。遇到相位批量偏移时先跑这一行比盲调参数有效得多。5.5 现象调和分析的残差在风暴季节呈现明显低频摆动原因天文潮是确定性分量但风暴增水和河流径流不是正弦过程。把强增水事件直接放进拟合序列会让低频部分受到明显扭曲也会影响浅水分潮结果。残差里看到的低频摆动正是这些气象和径流信号没被分离干净。解决识别出大风过程或极端水位时段将该段从拟合窗口中剔除或替换为随后一周的天文潮预报值。更稳妥的做法是采用一个月以上的逐时数据避开强增水事件使最小二乘结果接近纯天文潮。如果项目需要分析风暴潮则恰恰可以利用这个残差序列做后续的气象水位分离不要把残差简单看成噪声丢掉。6. 进阶把 t_tide 接进业务化流程——批量处理长序列与质量控制6.1 多年数据滚动分析逐年提取调和常数如果拥有多年连续观测不要一次性把全部数据丢进去算一个调和常数。潮汐调和常数会随时间缓慢变化工程上常按逐月或逐年窗口滑动观察振幅和相位的稳定性。滚动分析还能有效暴露单一年份缺失数据或强风暴污染的问题。amp_m2 []; years 2020:2024; for yr years idx tnum datenum(yr,1,1) tnum datenum(yr1,1,1); [tcon_y, ~] t_tide(wl(idx), start, tnum(find(idx,1)), interval, 1); m2_pos strcmp(tcon_y.name, M2); amp_m2(end1) tcon_y.amp(m2_pos); %#okSAGROW end这段代码把整段序列按年份切分每一年单独调用t_tide从返回的tcon_y里提取 M2 振幅存进数组。循环结束后可以看到逐年振幅没有明显跳变。正常情况下相邻年份间差应小于 5%。如果某年偏大一般说明这一年的数据包含了强增水事件或较多缺测。这种滚动式分析比直接算全时段平均值更健壮也更容易定位哪一段数据需要返工。6.2 质量控制的最后一步残差分析和结果记录调和分析的输出pred不应该直接作为成果还要与实测做差形成residual。分析残差的目的不只是看拟合效果更是为后续风暴潮或工程水位统计分离出干净的潮汐信号。残差的均方根误差和平台偏移直接反映调和常数提取的可靠程度。residual wl - pred; residual(residual 0) NaN; rms_res rms(residual(~isnan(residual))); fprintf(残差均方根误差%.3f m\n, rms_res); figure plot(tnum, residual, ., MarkerSize, 3) xlabel(时间); ylabel(残差 (m)) title(实测水位与天文潮合成后的残差序列)代码逻辑实测序列减去合成序列得到每个时间点残差rms计算残余序列的均方根打印出来。如果rms_res明显偏大说明还有未解释的周期成分或异常事件需要回到数据清洗阶段。画残差曲线时一条围绕零线窄幅波动的带宽说明天文潮提取充分如果残差里有波浪状起伏说明还有未分离的浅水分潮或增水影响。我在实际项目里养成的习惯是无论时间多紧都要先画出水位总览、时间间隔检查、缺测统计最后看残差。只有这三步都过了才会把tcon写进报告。从那以后我每次处理新站点潮汐分析都强制走一遍这个流程解压工具箱、检查路径、确认时间基准、跑默认结果、看残差再进入参数调优。这个流程能在两小时内排除八成数据坑。希望帮到你。本文还有配套的精品资源点击获取
