简介潮汐模型驱动程序TMD3.0版配套MATLAB代码包面向从事海洋潮汐模拟、水文预报及地球物理研究的高校师生与科研人员用于解决区域潮汐场建模、潮位预测与模型驱动计算等实际问题。压缩包共197个文件约93.55MB以113个png图示、38个m脚本、17个md说明文档和13个html页面为主另含nc数据文件、gif动画及xml、bib等配置与参考文献文件覆盖模型文档、格式说明与示例数据。已有345人学习下载。资源内含TMD模型文件格式说明、预测文档与EOT教程配合MATLAB脚本可直接复现潮汐驱动流程帮助读者理解模型输入输出结构、掌握潮汐预测实现方法并借助图示与动画直观核对结果适合作为潮汐建模入门与二次开发的参考素材。1. 潮汐模型驱动程序 TMD 3.0一个被低估的潮汐调和分析工具如果你手头有一年以上的水位观测序列想分离出八大分潮的振幅和相位又不想从零推导最小二乘的观测方程那 TMDTidal Model Driver大概率是你绕不开的一个工具。它最早由俄勒冈州立大学相关团队维护核心思路是把全球或区域潮汐模型的调和常数预先存成二进制格式再通过一套驱动程序按经纬度和时间插值出任意时刻的潮位。3.0 版把接口整理得更规整附带 MATLAB 调用代码意味着你不需要啃 Fortran 源码只要准备好模型文件和站位坐标就能在 MATLAB 里直接算出潮汐预报。它适合谁做近海工程水位校核的、跑海洋模式需要开边界潮汐驱动的、做验潮站数据质量控制的以及写论文需要批量算潮汐调和常数的。不适合谁只想画一张潮汐曲线图、对精度没要求的人用现成的在线预报就够了。这篇笔记按「TMD 3.0 到底怎么跑起来 → 模型数据怎么配 → 调和分析怎么落地 → 哪里容易翻车」的顺序讲代码能直接抄参数会逐个说明。2. 把 TMD 3.0 在 MATLAB 里跑通从解压到第一条潮位曲线2.1 先搞清楚 TMD 3.0 的目录结构和依赖关系拿到压缩包后别急着双击运行先看清楚里面有什么。TMD 3.0 的典型结构分三块一是核心驱动代码通常是 Fortran 或编译好的可执行文件负责读取模型数据并插值二是模型数据文件扩展名常见为.nc或自定义二进制格式存放各分潮的振幅和相位网格三是 MATLAB 封装层提供tmd_tide_pred、tmd_constit这类函数把底层调用包装成 MATLAB 能识别的接口。这里有个关键认知TMD 本身不是潮汐模型它只是「驱动程序」。模型数据才是真正的潮汐信息载体。所以你会看到两种文件要分开管理——驱动代码放一个目录模型数据放另一个目录MATLAB 路径只添加驱动代码所在目录模型数据路径通过参数传入。混在一起放不是不行但换模型时容易乱。依赖方面MATLAB 版本不需要太新R2016b 之后的版本基本都能跑因为用到的就是基础的矩阵运算和文件读写。如果你用的是 Windows 且核心是编译好的 exe注意系统架构要匹配32 位 MATLAB 调不动 64 位 exe这是血泪经验。Linux 下则通常需要自己编译 Fortran 源码gfortran是常见选择。2.2 最小可运行示例给定经纬度和时间算潮位下面这段代码是我一般会先跑的「冒烟测试」目的是确认驱动能读到模型、能返回数值、数值量级合理。假设模型数据文件叫TPXO9_atlas_v5.nc站位设在东海某点。% 添加 TMD 驱动代码路径按实际解压位置修改 addpath(genpath(D:\TMD3.0\driver)); % 模型数据文件路径 model_file D:\TMD3.0\data\TPXO9_atlas_v5.nc; % 站位经纬度东经、北纬为正 lon 122.5; lat 30.0; % 时间范围2024-01-01 00:00 到 2024-01-02 00:00每小时一个点 t_start datenum(2024,1,1,0,0,0); t_end datenum(2024,1,2,0,0,0); t_vec t_start:1/24:t_end; % 调用潮汐预报函数 [h, constituents] tmd_tide_pred(model_file, t_vec, lat, lon, z); % 输出前几个值看看量级 disp(h(1:5)); disp(constituents);逻辑说明tmd_tide_pred是封装好的预报函数第一个参数是模型文件第二个是 MATLAB datenum 格式的时间向量第三、四个是纬度和经度最后一个z表示要水位有些实现用h表示高度具体看函数签名。返回值h是潮位序列单位通常是米constituents返回参与计算的分潮列表。参数说明时间向量必须是 datenum 格式不是 datetime这是最容易翻车的地方。经纬度顺序是「先纬后经」和很多地理函数的习惯相反写反了会插值到错误网格甚至报错。如果你发现返回值全是 NaN先检查经纬度是否落在模型覆盖范围内——全球模型一般没问题区域模型就要看边界了。2.3 模型数据文件的加载与站位插值逻辑TMD 的插值不是简单的最近邻它会在模型网格上做双线性插值再对每个分潮的振幅和相位分别处理。振幅直接插值相位则要小心跨 0/360 度的跳变好的实现会做相位解缠。你不需要自己写这部分但要知道它存在因为当站位靠近网格边界或水深变化剧烈处插值误差会明显增大。模型文件加载是一次性的第一次调用会读进内存后续调用同一模型会快很多。如果你要批量算几百个站位建议把模型读一次然后循环站位而不是每个站位都重新加载。我见过有人写循环时把模型文件路径放在循环内反复读结果跑一晚上没跑完这就是没注意缓存机制。另外模型文件里的分潮数量决定了你能算哪些分潮。常见全球模型包含 M2、S2、N2、K2、K1、O1、P1、Q1 八个主要分潮有些还包含 M4、MS4 等浅水分潮。如果你的研究需要更多分潮得换模型或自己补充。MATLAB 封装层一般会返回实际参与计算的分潮名拿到后先扫一眼确认没有缺你需要的那个。3. 潮汐调和分析的落地从水位序列反算分潮参数3.1 调和分析的基本方程与 TMD 的角色区分调和分析的本质是已知水位序列反推各分潮的振幅和相位。方程是线性的每个时刻的水位等于各分潮贡献之和加上残差。TMD 在这里扮演两个角色——一是提供预报值用于对比二是它的模型数据本身就可以作为先验帮你判断反算结果是否合理。常见做法是用 TMD 预报一段时间的潮位再把这个预报值当作「真值」去检验你自己的调和分析代码是否正确。如果你从观测数据反算出的 M2 振幅和 TMD 预报的 M2 振幅差了一个量级那要么观测数据有问题要么你的分析代码有 bug。这个交叉验证思路非常实用比单纯看残差平方和直观得多。调和分析的观测方程可以写成矩阵形式A * x b其中A是设计矩阵每列对应一个分潮的 cos 和 sin 项x是待求系数b是观测水位。最小二乘解就是x (AA)^(-1) Ab。MATLAB 里用反斜杠运算符就能算不需要自己写求逆。3.2 用 MATLAB 实现最小二乘调和分析并对比 TMD 结果下面这段代码演示从合成水位序列反算分潮参数再和 TMD 预报对比。合成数据用 TMD 生成所以理论上反算结果应该非常接近。% 生成合成水位序列用 TMD 预报作为“观测” t_vec datenum(2024,1,1,0,0,0):1/24:datenum(2024,2,1,0,0,0); [h_obs, ~] tmd_tide_pred(model_file, t_vec, lat, lon, z); % 定义分潮频率周期单位小时 constits {M2,S2,N2,K2,K1,O1,P1,Q1}; periods [12.4206, 12.0000, 12.6583, 11.9672, ... 23.9345, 25.8193, 24.0659, 26.8684]; % 构建设计矩阵 N length(t_vec); A ones(N, 1); % 常数项平均水位 for k 1:length(constits) omega 2*pi / periods(k); % 角频率单位 1/小时 t_hours (t_vec - t_vec(1)) * 24; % 转成小时 A [A, cos(omega*t_hours(:)), sin(omega*t_hours(:))]; end % 最小二乘求解 coef A \ h_obs(:); % 提取振幅和相位 amp zeros(length(constits),1); pha zeros(length(constits),1); for k 1:length(constits) c coef(1 2*k - 1); s coef(1 2*k); amp(k) sqrt(c^2 s^2); pha(k) atan2(s, c) * 180/pi; % 转成度 end % 打印结果 for k 1:length(constits) fprintf(%s: 振幅%.4f m, 相位%.2f 度\n, constits{k}, amp(k), pha(k)); end逻辑说明设计矩阵第一列是常数项吸收平均水位之后每个分潮占两列分别是 cos 和 sin。时间转成「从起始时刻算起的小时数」是为了避免大数导致的数值精度问题。A \ h_obs直接给出最小二乘解MATLAB 内部用 QR 分解比显式求逆稳定。参数说明分潮周期必须用正确值M2 是 12.4206 小时不是 12 小时整差这一点会导致相位随时间漂移。相位用atan2而不是atan因为要区分象限。如果你发现某个分潮振幅异常大先检查该分潮周期是否写错或者观测序列长度是否足够——太短的序列无法分离相近频率的分潮比如 K1 和 P1 周期接近至少需要一个月以上数据才能分开。3.3 结果验证残差分析和与 TMD 预报的交叉对比算完调和常数后别急着出图。先做两件事一是看残差序列是否还有明显周期信号二是把反算的调和常数代回预报方程和原始序列对比。如果残差里还能看到明显的半日周期说明你漏了某个分潮或者分潮周期设错了。残差分析用 MATLAB 几行就能做h_fit A * coef; % 拟合值 resid h_obs(:) - h_fit; % 残差 fprintf(残差标准差: %.4f m\n, std(resid)); fprintf(残差最大值: %.4f m\n, max(abs(resid))); % 画图检查 figure; subplot(2,1,1); plot(t_vec, h_obs, b, t_vec, h_fit, r--); legend(原始,拟合); datetick(x,mm-dd); subplot(2,1,2); plot(t_vec, resid); datetick(x,mm-dd); title(残差序列);如果残差标准差在厘米量级说明拟合不错如果到了分米量级要么数据有噪声要么分潮没选全。交叉对比 TMD 预报值时注意 TMD 用的是模型网格插值而你反算用的是观测序列两者本身就有差异不要期望完全一致。合理的差异范围取决于站位处模型精度一般 M2 振幅差在几厘米内算正常。4. 避坑与排查TMD 3.0 跑不起来时先看这几条4.1 现象调用函数报「未定义函数或变量」原因MATLAB 路径没添加对或者添加了但函数名拼写和实际文件不一致。TMD 3.0 的 MATLAB 封装层函数名在不同版本里可能有细微差别比如tmd_tide_pred和tide_pred都出现过。解决先用which tmd_tide_pred确认 MATLAB 能不能找到这个函数。如果返回空检查addpath的路径是否指向了包含.m文件的目录而不是上一级目录。用genpath递归添加子目录更保险。如果函数名确实不对去解压目录里搜.m文件看实际名字是什么。4.2 现象返回的潮位全是 NaN 或零原因经纬度超出模型覆盖范围或者经纬度顺序写反了。区域模型尤其容易出这个问题全球模型一般不会。解决先确认模型覆盖范围通常模型文件里会有说明或属性。把经纬度限制在覆盖范围内再试。另外检查调用时的参数顺序TMD 的 MATLAB 封装一般是「先纬后经」和interp2的习惯相反。如果拿不准用两个明显不同的经纬度各跑一次看返回值有没有变化没变化就是没插值到正确位置。4.3 现象时间向量用 datetime 格式导致报错原因TMD 的 MATLAB 封装层大多按 datenum 设计datetime 是后来引入的格式不兼容。解决用datenum转换。如果原始数据是 datetime用datenum(dt)转一下。反过来如果函数返回的时间是 datenum你想用 datetime 画图用datetime(datenum_val, ConvertFrom, datenum)。这个转换在批量处理时容易忘建议在脚本开头统一转好。4.4 现象调和分析结果里某个分潮振幅特别大或特别小原因分潮周期设错或者观测序列太短导致分潮之间无法分离。K1 和 P1、S2 和 K2 这两对是重灾区。解决核对分潮周期表M212.4206、S212.0000、K123.9345、O125.8193这些值不能凭记忆写。序列长度至少覆盖 30 天最好 60 天以上。如果数据就是短考虑固定某些分潮的振幅相位关系比如用模型值约束而不是全部自由反算。4.5 现象Windows 下调用 exe 核心时闪退或无输出原因系统架构不匹配或者缺少运行时库。32 位 MATLAB 调 64 位 exe 会直接失败反之亦然。解决确认 MATLAB 和 exe 的架构一致用computer(arch)看 MATLAB 架构。如果 exe 是编译好的检查是否需要安装 Fortran 运行时库。实在搞不定就用纯 MATLAB 实现TMD 的插值逻辑不复杂自己写一个双线性插值加调和预报也就几十行可控性还更高。5. 进阶技巧用 TMD 做多站位批量预报与边界条件生成当你需要给海洋模式提供开边界潮汐驱动时单个站位预报就不够了得沿边界生成一系列站位的潮位时间序列。这时候效率是关键。我一般会先把模型读一次然后循环所有边界点每个点只做插值和调和求和不再重复读文件。具体做法是把边界点经纬度存成两列矩阵时间向量只建一次然后循环调用预报函数。如果 TMD 的封装层不支持批量传入多个站位就自己包一层循环。注意每次调用返回的潮位是列向量循环结束后拼成矩阵每列一个站位存成 NetCDF 或文本文件供模式读取。另一个技巧是分潮筛选。有些应用只关心 M2 和 S2那就在预报时只算这两个分潮速度会快很多。TMD 的封装层一般允许传入分潮列表参数查一下函数签名把不需要的分潮去掉。我见过有人为了算 M2 把八个分潮全算一遍白白多花七倍时间。验证批量结果时抽几个点单独跑一次对比批量输出是否一致。不一致通常是循环里变量没清干净或者经纬度索引写错了。这个检查花不了几分钟但能避免把错误边界条件喂给模式后跑几天才发现问题。最后说个习惯每次换模型数据文件先跑一遍单点冒烟测试确认量级和分潮列表正常再上批量。模型文件版本不同分潮命名和网格分辨率都可能变直接套用旧脚本容易翻车。希望帮到你。本文还有配套的精品资源点击获取
