简介:面向物理海洋学、水利工程与地理信息系统研究者的TMD2.05潮汐分析MATLAB工具箱,提供从潮汐数据解析、主要分潮识别到长期预测的完整计算流程,适合需要开展潮汐模拟、工程水位推算或风暴潮预警研究的科研与工程人员。压缩包共82个文件,以71个m脚本为主体,涵盖主程序TMD.m、子模型加载、谐波分析、椭圆参数提取及多种插值函数;另含lat_lon系列坐标文件、tpxo8_atlas模型文件与临时数据文件,整体约379KB。已有1012人学习下载。解压后在MATLAB中加载即可使用,内含GUI界面、tmd_tide_pred等预测函数及模型参数处理模块,可帮助快速建立潮汐模型、输出分潮调和常数,并结合tpxo8全球数据实现多站点预报,对海洋工程设计和潮汐灾害风险评估具有直接参考价值。
1. 聊 TMD2.05.zip,先聊你手里的潮汐数据为什么算不准
如果你写过海洋、港口、海岸工程或海上风电相关的代码,大概率会遇到这个场景:甲方只给了你一个 GPS 坐标,却要求你给出未来一年每小时的潮位预报,甚至要画出某海域的潮汐流矢量场。MATLAB 用户圈子里,TMD(Tidal Model Driver)2.05 一直是个口口相传的工具包——它已从创始人 R. S. Pawlowicz 那里以TMD2.05.zip的形式流传了十几年,专门用于读取全球海潮模型(最常见的是 TPXO 系列),把其中有调和常数提取出来,再对任意时空点做潮位、潮流、潮汐分潮的复现与预测。它不靠插值网格点上的经验值,而是按球坐标做空间插值后调用格林函数或直接使用数据集内的分潮常数,精度和计算效率在工程场景里都优于自己手写调和分析代码。它能解决的是「给定坐标和时间,快速给出逐时潮位与潮向」这类问题,适合水动力工程师、海事算法工程师、避航研究和数值预报组的同学。标题里的 2.05 指代的是一个稳定发行版,网上流传的 zip 包结构基本一致,但解压方式、路径配置和模型数据存放位置各有说法——这篇我们就从实际使用者的角度把这一整套流程讲透。
2. TMD 潮汐计算的理论骨架与数据集关系
2.1 TMD 2.05 到底做了什么:从调和常数到潮位合成
TMD(Tidal Model Driver)的核心并不是去解 Navier-Stokes 方程,它是把前人计算好的潮汐模型数据作为「常数库」读取,然后利用潮汐调和分析的公式综合出任意时刻的潮位。基本公式是(对你没看错,这是一个叠加式子):
h(t) = amp0 + sum_n( f_n * A_n * cos(omega_n * t + V_n + u_n - g_n) )其中A_n是第 n 个分潮的振幅,g_n是格林尼治相位滞后,omega_n是分潮角速度,V_n是天文初相位,f_n和u_n是交点因子和交点订正角,amp0是平均海平面。TMD 在运行时会先从网格数据中读取对应经纬度上的A_n、g_n、amp0,然后基于你给定的时间序列计算V_n、f_n、u_n,最终把多个分潮结果叠加起来。
注意:TMD 不负责做调和分析(就是把实测水位拆成多个分潮),它负责的是用已有的调和常数组去做正演计算。如果你的需求是「自己分析实测水位」,请搜 t_tide 工具包,TMD 与 t_tide 是配套使用的:t_tide 做逆向分析,TMD 做正向预测。
2.2 为什么是 2.05 版本而不是新版本
在 GitHub 上你能找到更新到了 2021 年左右的 GitHub 镜像仓库(如 tmd 的镜像),但TMD2.05.zip是流传度最高的一个打包版本。原因有两点:第一,2.05 之后的版本改动主要是数据结构适配新的海洋潮汐模型 (如 TPXO9),但很多老项目中模型文件是 TPXO7.2 或 TPXO8.0,2.05 版本对这些老模型兼容得更好;第二,网上大多数博客、教学视频、开源代码仓库引用的路径和函数签名(如tmd_predict、tmd_interp、tmd_extract_HC)都是在 2.05 上验证过的,对新人有实打实的参考价值。
这里我不建议你盲目追求新版本,先把你拿到的模型数据格式确定下来:如果是 NetCDF 格式的辐合数据,且模型自带h、u、v、hZ、uZ等变量,2.05 足够用。如果你要读 TPXO9-Atlas 或 FES2014 这类新型模型,你可能需要 smedigau 的TMDfork(即 GitHub 上的smedigau/tmd),那个版本改了数据读取层,但函数名基本沿用了 2.05。
2.3 TMD2.05.zip 核心文件清单
解压后你会看到这样一个结构(常见做法是丢进 MATLAB 搜索路径):
TMD_2.05/ tmd_predict.m tmd_interp.m tmd_extract_HC.m tmd_ellipse.m tmd_conlist.m tmd_equilibria.m tmd_plot.m tmd_convert_ll.m tmd_coastline.m tmd_constant.h load_data.m DATA/ Docs/关键函数角色如下:
| 函数名 | 作用 | 典型调用方式 |
|---|---|---|
tmd_predict | 根据模型和测站坐标、时间,输出潮位 | tmd_predict(model, lat, lon, time) |
tmd_interp | 在模型网格上插值出某点调和常数 | tmd_interp(model, lat, lon, 'h') |
tmd_extract_HC | 提取某点所有分潮的振幅与相位 | hc = tmd_extract_HC(model, lat, lon) |
tmd_ellipse | 绘制/计算潮流椭圆参数 | tmd_ellipse(hc, 'h') |
tmd_conlist | 列出模型里包含哪些分潮 | tmd_conlist(model) |
load_data.m是一个交互式脚本,它的作用是从美国俄勒冈州立大学(OSU)的服务器上下载特定区域的 TPXO 模型数据。但由于官方服务器和网络环境的未知性,你经常会遇到下载超时或只能拿到部分网格的情况,所以我建议你在解压 zip 后另外手动找到你所在海域的范围文件,常见做法是去tmd的仓库 Release 页面或者 OSU 的 TPXO 数据分发中心手动下载一个h_*_*.nc文件,再配合模型的grid元数据处理。
3. 从 zip 到第一张潮汐图:完整复现步骤
3.1 安装与路径配置的三种方式
先说结论:TMD2.05.zip实际上不是自动安装程序,而是一个源码打包,你没有必要跑install脚本,只需保证 MATLAB 能索引到目录即可。三种方式任选:
%% 方式一:一次性添加路径 addpath(genpath('D:/tools/TMD_2.05'));%% 方式二:保存为默认路径 pathtool; % 手动添加 TMD_2.05 与 DATA 文件夹,然后 save%% 方式三:启动脚本统一加载 % 在你的 startup.m 中添加以下行 addpath(genpath('D:/tools/TMD_2.05')); tmdver = tmd_conlist('D:/DATA/h_global_tpxo7.2.nc'); % 测试读取这里的tmd_conlist会读取模型文件并列出内部包含的分潮名列表,若能正常返回,说明你的路径和数据文件都OK。genpath会把子目录全部加入,包括DATA和Docs,但你后续真正读取的模型文件不一定放在DATA下,建议单独用变量管理。
提示:避免使用
addpath(genpath(pwd))这类写法,因为你一旦切换当前文件夹,工具箱路径会失效,且可能出现两个同名函数冲突。
3.2 模型数据文件:z 文件、netcdf 与 struct 三种格式的取舍
TMD 兼容三种模型数据形式:
- OSU 原始格式:通常以一堆无扩展名或
.nc结尾的文件组成,每月一个分潮一套文件(比如h_m2.out、h_s2.out)。这种格式年代较老,但读取速度快。 - NetCDF 格式:对于 TPXO7.2、TPXO8.0 和 TPXO9 来说,常用一个
h_global_tpxo7.2.nc同时装下几十个分潮,TMD 2.05 需要 MATLAB 的netcdf函数支持。老版本 MATLAB 用ncgeodataset(来自 NCTOOLBOX),新版本用内置ncread系列即可。 - 已加载的 struct:用
tmd_interp之后返回的model是一个 MATLAB struct,里面存有lat、lon、h、u、v、Z等网格数组。你可以把 struct 通过save变成.mat文件,之后直接用load读取,这是速度最快的方式,适合你反复做不同时间段的预测。
实际操作时,我强烈建议你拿到.nc模型后先转换成.mat再复用:
% 将 netcdf 模型封装为 mat 中的 struct model = tmd_interp('D:/DATA/h_global_tpxo7.2.nc', lat, lon, 'h'); save('D:/DATA/tpxo7.2_model.mat', 'model', '-v7.3');注意-v7.3参数:TPXO 模型的网格是 1440x721 甚至更细,数据量可能会超过 2GB,-v7.3使用 HDF5 底层格式,支持超过 2GB 的变量。如果你不写这个选项,保存过程很可能直接报错或让你的内存爆炸。
3.3 用 tmd_predict 复现某站逐时潮位
这是模拟潮汐最核心的一步。假设你要计算东经 122.2 度、北纬 37.3 度(威海附近),2024 年 1 月 1 日 0 点到 1 月 3 日 0 点时段的逐时潮位:
% 加载模型 struct(假设已经转成 .mat) load('D:/DATA/tpxo7.2_model.mat', 'model'); % 构建时间序列:每 1 小时一个点 t0 = datenum(2024, 1, 1, 0, 0, 0); % MATLAB 时间格式 t1 = datenum(2024, 1, 3, 0, 0, 0); time = t0:1/24:t1; % 步长为 1 小时 % 预测潮位,单位:米(相对于模型海平面,通常为平均海平面) h = tmd_predict(model, 37.3, 122.2, time); % 画图 plot(time, h); datetick('x', 'yyyy-mm-dd HH:MM'); xlabel('时间'); ylabel('潮位(m)'); grid on;这段代码有三个值得关注的参数:
- 纬度在前、经度在后,别传反了,
tmd_predict(model, lat, lon, time)的顺序是固定好的。用tmd_convert_ll可以帮你把经纬度在+-180和0~360之间转换,但大多数情况下直接用122.2(东经为正)是没有问题的。 tmd_predict内部会识别time为 MATLAB 的datenum形式,不是datetime。如果你在 R2019b 之后习惯了datetime,务必先用datenum(...)显式转换,或者在调用前检查class(time)。- 输出的
h是一个列向量,行列方向由输入的time决定:若time是行向量,则h也是行向量。这点常被忽视,导致后续拼接数据维度报错。
运行成功后,你会看到一条明显的正弦波动曲线。如果画图结果是一条几乎水平线,请检查模型数据是否覆盖该坐标点(在海里才算有效,陆地上 TMD 返回的插值结果会是NaN),或者你的时间序列长度是否不足一整个潮汐周期。
4. 潮汐模拟进阶:调和常数提取、分潮订正与效应物理解读
4.1 tmd_extract_HC 返回的 12 个分潮如何解读
TMD 可以使用 TPXO 数据里全部分潮,也可以只算前几个。默认情况下,tmd_predict会根据你给的模型自动决定分潮数目。如果你想查看具体每个分潮的贡献,就必须显式提取调和常数:
% 提取威海附近的分潮调和常数 hc = tmd_extract_HC(model, 37.3, 122.2); % hc 里有哪些字段 fields(hc) % 输出通常包括: amp 和 phase % hc.amp 是 [n_freq x 1],hc.phase 是 [n_freq x 1] % 对应的频率名在 hc.name 里常见输出示意表(数值只是示例,不要当作真实数据):
| hc.name | hc.amp (m) | hc.phase (deg) | 意义 |
|---|---|---|---|
| M2 | 0.85 | 127.3 | 主要半日分潮 |
| S2 | 0.32 | 141.8 | 太阳主要半日分潮 |
| N2 | 0.18 | 112.9 | 椭圆摄动引起的半日潮 |
| K1 | 0.42 | 88.5 | 日月合成的日潮 |
| O1 | 0.25 | 74.2 | 主太阴日潮 |
| M4 | 0.03 | 120.0 | 浅水分潮(非线性效应) |
hc.amp单位是米,相位单位是度(格林尼治相位迟角)。这个数据结构对后续分析很关键:tmd_predict是「外部封装的合成器」,而tmd_extract_HC能给你更细粒度的调参入口。如果你想做自定义的潮汐计算——比如只算 M2 + S2 + K1 + O1 四个主分潮,那么你可以这样做:
idx = ismember(hc.name, {'M2', 'S2', 'K1', 'O1'}); custom_amp = hc.amp(idx); custom_phase = hc.phase(idx); % 再结合你自己的格林尼治相位累计逻辑进行正演4.2 调和常数的空间分布:从点到位
只拿一个点做潮位线还不够,潮汐模拟的核心优势在于空间全场。TMD 提供tmd_interp把模型网格插值到任意坐标点,你可以对一个区域做网格循环,得到该区域的分布场。常见做法是生成网格坐标矩阵再向量化插值:
lons = 120:0.2:124; % 经度范围 120E ~ 124E lats = 35:0.2:39; % 纬度范围 35N ~ 39N [LON, LAT] = meshgrid(lons, lats); % tmd_interp 需要 reshape 成列 lonv = LON(:); latv = LAT(:); % 插值 M2 分潮的振幅和相位 hM2 = tmd_interp(model, latv, lonv, 'h'); % 第五参数缺省时返回振幅 % 更常用:直接返回所有变量结构 [s, ~, ~] = tmd_extract_HC(model, latv, lonv); % 不推荐在大网格上循环这里有两个性能陷阱:
- 不要用 for 循环逐点调用
tmd_extract_HC。TPXO 网格是 0.25 度分辨率,一个渤海区域就有上千个格点,逐点调用会重复打开 netcdf 文件句柄,耗时可能从秒级退化到分钟级。正确方法是调整model里的网格向量,利用tmd_interp一次传入全部坐标点。 model里的lat、lon本来就是一维网格坐标,tmd_interp内部会对输入纬度和经度做双线性插值。你传入的latv、lonv必须是一维向量,tmd_interp会自动把它们当作散点列表处理,而不是要求它们构成单调序列。
插值完成后,你可以用pcolor或contourf画出 M2 振幅等值线,这对雷雨天气下的风暴潮预报叠加非常有意义。
4.3 tmd_ellipse 绘制潮流椭圆与旋转方向
如果你做的是潮流模拟而不只是潮位模拟,就需要关注当前测点的潮流椭圆参数。tmd_ellipse根据当前点的u(东向分量)与v(北向分量)分潮调和常数,计算出椭圆长半轴、短半轴、旋转方向和分潮起始时刻的方向:
% 计算潮流椭圆 ell = tmd_ellipse(hc_u_v); % hc_u_v 需从 tmd_extract_HC 中获得 % ell 输出字段:semi_major, semi_minor, incl, phase如果你是从tmd_extract_HC(model, lat, lon, 'u')或'v'分别拿到流速分量,就还需要把它们组合成复速度形式。实测中更快的路径是直接把model中某点u、v插值出来,写成这样:
uM2 = tmd_interp(model, lat, lon, 'u'); vM2 = tmd_interp(model, lat, lon, 'v'); uv_complex = uM2 + 1i * vM2; % 复数方式表示速度这种复数表示在 TMD 中非常常见,因为潮汐流的调和常数据天然带有方向信息,复数运算可以直接体现顺时针或逆时针旋转。
5. 潮汐预测实战:未来多日预报与区域批量计算
5.1 构建任意时段的预测时间轴
潮汐预测的数学基础是固定的,时间轴是否正确决定了预报质量。你要注意三件事:
- 时间步长要与需求一致:工程上常用每分钟或每小时一个点,但潮汐谱的高频分潮只有 M4、M6 等少数浅水分潮,你用 6 分钟步长并不会让结果更准确,反而使计算量增大 10 倍。我一般做港口作业窗口分析用 10 分钟步长,做航海分析用 1 小时步长。
- 起始时刻从整点开始,时间步的单位是「天」,而不是秒。
1/24表示 1 小时,1/1440表示 1 分钟,这个很容易写错。 - 跨年的时候注意
datenum对闰年判断是准确的,你只需要关心 MATLAB 的日期边界,不存在 UTC+8 导致的时区问题——TMD 是基于静态时间坐标的,不涉及时区转换。
% 预测:2025-03-01 00:00 至 2025-03-08 00:00,10 分钟间隔 step = 1/1440 * 10; % 10 分钟对应的天数 t_start = datenum(2025, 3, 1); t_end = datenum(2025, 3, 8); time = t_start:step:t_end; h_pred = tmd_predict(model, lat, lon, time);运行这段代码前先检查内存:1441 * 7 = 10080 个点,完全没问题;但如果结合区域批量,比如 200 个站,那么就是 200 万数据,输出h_pred是 200 行 10080 列的双精度矩阵,占内存 16MB,也不大,放心算。
5.2 批量站点数据组织与并行化
真实项目中,测站往往有多个且坐标存在于 CSV、Excel 或数据库中。推荐做法是先把站点表读成表格:
% 读取站点表:站点ID,纬度,经度 stations = readtable('stations.csv'); lats = stations.lat; lons = stations.lon; ids = stations.id; % 预分配 n_station = height(stations); n_time = length(time); h_all = zeros(n_station, n_time); for i = 1:n_station h_all(i, :) = tmd_predict(model, lats(i), lons(i), time); end如果你的 MATLAB 版本支持parfor,且你有多核 CPU,可以改写成:
parfor i = 1:n_station h_all(i, :) = tmd_predict(model, lats(i), lons(i), time); end但注意:parfor需要提前开启并行池(parpool),并且每个 worker 必须能访问到model。model是只读数据,不会有数据竞争。如果你的站点数少于 50,并行加速效果不明显,反而要付出约 1-2 分钟启动池的时间,自己掂量。
5.3 预测结果的验证、赢点与易错点
预测值和实测值对比时,常见误差来源不是 TMD 本身,而是模型空间分辨率不够。TPXO7.2 的网格约 0.25 度(大约 25 公里),在近岸地形急剧变化的海域(比如长江口、胶州湾、狭窄水道),25 公里分辨不出地形效应,潮位偏差可能达到 20-30 厘米。TMD2.05.zip本身只是驱动,不包含高分辨率局部模型。这种情况,常见做法是先用 TMD 跑出开边界条件,再用 DELFT3D 或 FVCOM 做局部加密数值模拟。
另一大易错点:我们前面提到tmd_predict返回的是相对于平均海平面的潮高,如果你拿到的现场潮位基准是当地理论最低潮面(chart datum),那两者之间会有一个固定的基准差。你需要从当地潮汐表查这个基准差,然后加到预测结果上:
% 预测结果从平均海平面转移到理论最低潮面(示意) h_chart_datum = h_pred - constant_offset;constant_offset具体是多少、符号是正是负,完全取决于当地基准定义,不要靠猜。一般港口局提供的潮高基准为「理论上可能最低潮面」,这个值小于平均海平面,所以h_pred减去该常数才会和实际潮高表对得上。你可以用当地的潮汐表当天的潮高和 TMD 预测值做回归,估算这个常数。
提示:如果校验结果中偏差存在明显的周期特征(比如每天有两个峰值),说明你的基准差不是常数,而是你还漏掉了日潮或半日潮的相位订正——建议检查
tmd_conlist输出的分潮列表中是否有K2或T2等微小分潮,TPXO7.2 默认可能不包含它们。
6. 一个小技巧:用tmd_predict与tmd_interp组合,反推任意时刻潮位梯度场
做水动力数据同化或海面地形分析时,你有没有试过把 TMD 当成一个「纯函数」连续调用?我常用的一个技巧是:并不直接等距离散站点,而是先用tmd_interp插值出一个大区域的网格调和常数,然后在每个时刻对调和常数做加权叠加。这样可以避免每个站点、每个时间戳都重复做一次tmd_predict的沉重插值计算,单位时间从O(Nsite)降到O(Ngrid),在绘制潮位场动画时非常划算。
% 1. 构建目标网格(统一分辨率,比如 0.1 度) lon_g = 120:0.1:124; lat_g = 36:0.1:40; [lon_mesh, lat_mesh] = meshgrid(lon_g, lat_g); % 2. 插值所有分潮的调和常数(振幅和相位都存下) hc_grid = tmd_extract_HC(model, lat_mesh(:), lon_mesh(:)); % 该函数返回的 hc_grid.amp 尺寸为 [N * nfreq] % 3. 对每个分潮做时间合成 time_vec = t_start:1/24:t_end; for it = 1:length(time_vec) % 时间正演项:需要了解 tmd_equilibria 计算每时刻的均衡参数 % 这里简化处理,示意思路: % 真实情况请读取 hc_grid.amp 与相位后自行合成 h_field = zeros(size(lat_mesh)); % 循环 nfreq 个分潮,叠加到 h_field % h_field(:) = h_field(:) + hc_grid.amp(:, ifreq) * cos(...) ... % 再显示或保存 end完整的精密合成需要处理f、u、V等天文参数,这其实被 TMD 封装在tmd_equilibria.m里。你可以直接调用它拿到各分潮在指定时刻的V0+u值:
[V0, u, f] = tmd_equilibria(hc_grid.con, time_vec);这样你完全绕开了tmd_predict对坐标点循环的约束,把计算变成了两个阶段:一次网格插值、N 次时间合成。这个技巧特别适合你后期做数据可视化或向 Web 端口提供实时预报接口——把调和常数缓存为.mat或数据库行,前端只做时序计算。如果哪天你发现某个坐标附近预测结果与实测始终差了一个固定角度,记得回头检查从tmd_interp取出的 phase 单位到底是度还是弧度,还有模型里con列表的分潮排列顺序是否与tmd_equilibria返回的数组严格对应——这两处是最隐蔽的系统性误差来源。
本文还有配套的精品资源,点击获取