简介面向物理海洋学、水利工程与地理信息系统研究者的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 用户圈子里TMDTidal Model Driver2.05 一直是个口口相传的工具包——它已从创始人 R. S. Pawlowicz 那里以TMD2.05.zip的形式流传了十几年专门用于读取全球海潮模型最常见的是 TPXO 系列把其中有调和常数提取出来再对任意时空点做潮位、潮流、潮汐分潮的复现与预测。它不靠插值网格点上的经验值而是按球坐标做空间插值后调用格林函数或直接使用数据集内的分潮常数精度和计算效率在工程场景里都优于自己手写调和分析代码。它能解决的是「给定坐标和时间快速给出逐时潮位与潮向」这类问题适合水动力工程师、海事算法工程师、避航研究和数值预报组的同学。标题里的 2.05 指代的是一个稳定发行版网上流传的 zip 包结构基本一致但解压方式、路径配置和模型数据存放位置各有说法——这篇我们就从实际使用者的角度把这一整套流程讲透。2. TMD 潮汐计算的理论骨架与数据集关系2.1 TMD 2.05 到底做了什么从调和常数到潮位合成TMDTidal 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.02.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.namehc.amp (m)hc.phase (deg)意义M20.85127.3主要半日分潮S20.32141.8太阳主要半日分潮N20.18112.9椭圆摄动引起的半日潮K10.4288.5日月合成的日潮O10.2574.2主太阴日潮M40.03120.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 的日期边界不存在 UTC8 导致的时区问题——TMD 是基于静态时间坐标的不涉及时区转换。% 预测2025-03-01 00:00 至 2025-03-08 00:0010 分钟间隔 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里。你可以直接调用它拿到各分潮在指定时刻的V0u值[V0, u, f] tmd_equilibria(hc_grid.con, time_vec);这样你完全绕开了tmd_predict对坐标点循环的约束把计算变成了两个阶段一次网格插值、N 次时间合成。这个技巧特别适合你后期做数据可视化或向 Web 端口提供实时预报接口——把调和常数缓存为.mat或数据库行前端只做时序计算。如果哪天你发现某个坐标附近预测结果与实测始终差了一个固定角度记得回头检查从tmd_interp取出的 phase 单位到底是度还是弧度还有模型里con列表的分潮排列顺序是否与tmd_equilibria返回的数组严格对应——这两处是最隐蔽的系统性误差来源。本文还有配套的精品资源点击获取