ARTICLE DETAIL

资讯详情

深耕网站建设、视觉设计与SEO优化的一线实战洞察。

Matlab潮汐调和分析实战:t_tide数据预处理、参数配置与预报验证

Matlab潮汐调和分析实战:t_tide数据预处理、参数配置与预报验证 简介这是一套基于 MATLAB 的潮汐调和分析实例工具包面向海洋学、地球科学等领域的科研人员与工程师用于将实测水位时间序列分解为不同频率的天文潮分量并开展潮汐预报与结果评估。资源共 18 个文件压缩包约 829KB其中包含 11 个 .m 脚本如 t_tide.m、t_predic.m、t_demo.m 等承担核心调和分析、预测、误差统计与示例运行功能5 个 .mat 数据文件存储潮汐常数、天文参数等基础数据集2 个 .dat 文件提供平衡潮与实测潮位的原始数据输入。目前已有 856 人学习使用。通过结合 t_demo.m 演示流程用户可快速掌握从数据导入、模型构建、最小二乘参数估计到残差检验、曲线绘制的完整分析方法并能基于 t_predic.m 进行后续潮位预测适用于课程研究、论文实验及工程潮汐分析场景。1. 潮汐调和分析为什么绕不开 Matlab 的 t_tide潮汐调和分析的日常任务是把一段实测水位拆成多个固定频率的余弦波每个波对应的振幅和迟角就是调和常数有了调和常数未来任何时刻的潮位都能直接推算。Matlab 生态里绕不开的实现就是 t_tide2002 年发布至今一直是海洋工程、港口设计和物理海洋数据处理里的默认选项。它的优点很突出接口简单、内置分潮表完整、还带误差估计但文档零散新手常在数据格式和时间基准上翻车。下面从数据约定、调用语法、参数调整到预报验证四个环节把 t_tide 的一整条链路说清楚。适合手头有水位或海流时间序列、想快速得到可用的调和常数和预报曲线的工程师和研究生。2. 潮汐调和分析的原理与 t_tide 数据约定2.1 分潮、振幅与迟角到底在算什么调和分析的数学结构不复杂把水位观测写成平均海面项加上一串余弦项的叠加每一项的频率由天体运动决定叫做分潮角速度。平时常说的 M2 是主太阴半日分潮S2 是主太阳半日分潮K1 和 O1 是主要全日分潮浅水区还有 M4、MS4 这些倍频分潮。t_tide 内置的分潮表远比这几个常用分潮大默认会根据数据长度自动筛选能独立求解的项。振幅的单位和输入一致输入米就得到以米为单位的振幅迟角是相对格林尼治天文相角的相位差单位是度。振幅和迟角合起来构成调和常数它是站点潮汐特征的定数同一站点不同年份分析出来的调和常数应当保持稳定。做长序列质量评估时连续两年的调和常数差超过误差范围就要警惕观测设备或基准面出了问题。迟角不是余弦函数里那个初始相位这一点特别容易混淆。t_tide 输出的迟角已经扣除了起始时刻的天文相角因此不同起始时间的数据算出的迟角不能直接比较。跨时段对比时必须用天文相角把两套结果换算到同一相位基准上否则会得到“滞后了好几十度”的假结论。这也解释了为什么两个团队用同一批数据只要 start 参数写法不同相位结果就可能对不上。2.2 数据长度决定能分辨哪些分潮分潮之间的频率间隔不同分辨它们需要的数据长度也不同。判断依据是最小分辨时长约为两个分潮差频的倒数。比如 M2 和 S2 的角速度差约为每小时 1.016 度对应周期约 14.8 天所以至少半个月的连续数据才能把这两个半日分潮分开。K1 和 P1 的差频很小要求数据长到 200 天以上这也是为什么短期验潮分析一般不给 P1 单独列结果。t_tide 自动做这个判断数据不够长就把无法分辨的分潮排除在最小二乘解之外。强行手动把分潮表拉长法方程接近奇异振幅和迟角会出现数值上虚高甚至负值但程序不会报错。想快速判断当前数据长度适合哪些分潮可以参考表 1 的经验值或者直接看 t_tide 输出里实际包含了哪些分潮。表 1 常用分潮对的最低分辨时长分潮对角速度差度/小时最低分辨时长O1 / K13.0504.9 天M2 / S21.01614.8 天N2 / M21.8818.0 天K1 / P10.072208 天K2 / S20.082183 天表格里是理论下限实际项目里建议按两倍时长准备数据否则误差带会宽到工程上没法用。另外要注意这里说的“数据长度”是有效连续观测长度中间有大量缺口时要折算成等效连续长度。2.3 输入数据格式与时间基准t_tide 的基本输入是等间隔采样序列配合采样间隔和起始时间。起始时间必须用 datenum 数值格式它参与天文相角计算填错或给成字符串程序虽然能跑但输出迟角会整体偏移。交叉验证里常见的“相位对不上”十有八九是起始时间基准没对齐。序列里的 NaN 会被 t_tide 跳过但缺口过大会导致法方程的数据量不足误差估计随之变大。常见做法是先把数据整理成干净、等间隔的序列再进入分析。下面这段代码从原始 CSV 读入并检查时间间隔% 读取原始观测并检查等间隔性 data readmatrix(tide_gauge.csv); % [datenum, 水位] 两列 t_num data(:,1); z data(:,2); dt_h diff(t_num) * 24; % 转成小时差 fprintf(间隔范围: %.4f ~ %.4f 小时\n, min(dt_h), max(dt_h)); if max(dt_h) 1.5 * median(dt_h) % 存在缺口按中位数间隔重采样 dt median(dt_h); t_new t_num(1) : dt/24 : t_num(end); z_new interp1(t_num, z, t_new, linear); t_num t_new(:); z z_new(:); else z z(:); end用 median 而不是 mean 作为采样间隔的基准是因为个别缺失点会把平均间隔拉大而中位数能反映“正常情况下的采样频率”。interp1 的线性插值对潮位这种平滑序列足够不会明显压低 M2 和 K1 的振幅。预处理结束后z 是列向量、单位是米、时间已经对齐到等间隔网格可以进入正式分析。3. 用 Matlab 跑通 t_tide 的完整分析流程3.1 安装 t_tide 并配置路径t_tide 不是 Matlab 官方工具箱要先拿到源码再添加路径。源码压缩包解压后会有 t_tide.m、t_predic.m 等核心文件和一些示例数据。把整个目录放到固定位置然后执行addpath(genpath(D:\toolboxes\t_tide)); savepath;genpath 递归添加目录下的所有子文件夹savepath 把当前路径写入 Matlab 的路径缓存之后不用每次启动重新添加。装完可以用which t_tide确认如果返回带完整路径的 .m 文件就算安装成功。这个工具箱代码写于本世纪初基本语法保持得很好最近几个大版本 Matlab 跑同一份代码没有出现兼容问题。需要注意 Matlab 后续版本对某些字符串语法更严格如果报错指向字符串拼接把单引号字符串改成双引号即可。安装阶段的另一个坑是路径里有中文空格genpath 处理不了建议放纯英文路径。3.2 水位和海流数据的预处理预处理有四个标准动作去野值、补缺口、去趋势可选、统一基准面。野值可以用滑动中位数找出来再用线性插值替换。对小时数据来说窗口 25 个点大约对应一天能识别出单点跳变又不会把真正的潮汐峰谷当野值。去趋势看用途。要做高程衔接、保留平均海面项时别去趋势让 t_tide 自己估计常数项如果主要关心分潮振幅和相位可以去掉 30 天滑动平均对应的长周期信号减少非潮汐水位变化对最小二乘的干扰。基准面统一是另一个常被忽略的环节。t_tide 只分析相对变化不同基准面不影响分潮振幅和迟角但影响常数项。如果要把多站点的调和分析结果作对比先确认各站高程基准一致否则常数项差异会被误读成潮差差异。预处理示例z_med movmedian(z, 25); bad abs(z - z_med) 3 * std(z - z_med); z(bad) NaN; z fillmissing(z, linear); % 可选去掉30天滑动平均趋势 z_detrend z - movmean(z, 24*30);先 movmedian 做中值滤波再用 3 倍标准差作阈值找野值。阈值选取要看站点潮差背景潮差大的海区可以放宽到 5 倍标准差避免大潮期间的正常高水位被误删。数据量大时也可以不插值保留有效段逐段分析最后对调和常数做加权平均。3.3 t_tide 核心调用与参数说明核心调用就一次[tidestruc, xout] t_tide(z, interval, 1, ... start, t_num(1), latitude, 30.5, ... infer, true, error, true);interval 是采样间隔单位小时10 分钟数据填 1/6一小时数据填 1start 是起始时刻的 datenum 值latitude 是观测纬度用于计算交点因子infer 开启分潮推断error 开启误差估计默认就是开的建议保持。参数用“参数名, 值”成对传入。表 2 整理了常见参数和推荐取值。参数默认值推荐设置说明interval1按实际采样填填错全盘皆错最容易翻车start无datenum 格式影响相位归算必须有latitude0实测纬度高纬度地区交点因子差异变大inferfalsetrue短序列时补 P1、K2 等相邻分潮shallowfalse浅水站 true增加 M4、MS4 等浅水分潮errortrue保持 true输出置信区间和 SNRoutputnone需要时改 full输出更多诊断信息漏掉 interval 是最常见的坑程序不会报错因为默认 1 小时如果你的序列是 10 分钟采样分潮频率全被解读错了。判断方法很简单跑完之后对比 xout 与原始序列如果相位漂移得很规律先检查 interval。3.4 读懂 t_tide 输出结构输出结构体 tidestruc 的核心字段是 tidecon、name、freq。tidecon 是一个 N 行乘 6 列的矩阵每行对应一个分潮列含义依次是振幅、振幅误差、迟角、迟角误差、SNR、是否纳入解算。name 是分潮名称列表freq 是角速度。表 3 给出常用字段速查。字段维度内容tidestruc.nameN×1 cell分潮名称如 M2、K1tidestruc.tideconN×6 double振幅、迟角与对应误差、SNR、标志tidestruc.freqN×1 double分潮角速度度/小时tidestruc.datum1×1平均海面项xout与输入等长拟合潮位序列拿到结果第一步不是看每个分潮振幅大小而是看 SNR。SNR 小于 1 的分潮基本不可信t_tide 默认不显示SNR 在 1 到 2 之间只能当参考。下面这段代码分潮结果整理成表格并按 SNR 排序tbl table(tidestruc.name, ... tidestruc.tidecon(:,1), ... tidestruc.tidecon(:,5), ... VariableNames, {分潮, 振幅_m, SNR}); tbl sortrows(tbl, SNR, descend); disp(tbl(1:10, :));输出表的前几行应该是 M2、K1、S2、N2、O1 这些主分潮SNR 至少在几十以上。如果 M2 的 SNR 很低或振幅与邻近站差一个数量级不要急着进入预报回头查数据质量。重点放在 M2 和 K1 上先确认它们没有异常再看其他分潮。4. t_tide 参数调整、推断分潮与常见坑4.1 值得细调的参数与阀门interval、start、latitude 属于必填项shallow、infer、error 属于按需调节。还有一个容易被忽视的 synthesis 参数它在 t_tide 和 t_predic 里含义略不同。在 t_tide 里synthesis 决定返回的 xout 是只用主分潮合成还是把推断分潮也包含进去。做预报时建议取 1让推断分潮参与合成曲线更平滑且接近真实潮型。纬度对结果的影响在高纬度更明显。交点因子与纬度相关北纬 60 度和赤道附近同一个分潮的振幅修正最多能差到 10%。如果项目站位靠近极区latitude 必须填实测值同时注意 t_tide 内部对交点的处理用的是与纬度相关的近似公式站在地磁异常区时误差会放大这一点在源码注释里有说明。shallow 参数决定是否包含浅水分潮。水深小于 20 米的河口和浅滩M4、MS4 不可忽略它们会使潮汐曲线出现明显的不对称表现为涨潮短、落潮长。做航道通航水深预报时浅水分潮必须保留否则低潮水位预测会系统性偏大。4.2 短序列用 inference 推断邻近分潮工程里最常见的情况是只有两个月水位数据拿不到 P1 和 K2。此时开启infer, truet_tide 会用 K1 推断 P1、用 S2 推断 K2按固定的振幅比和相位差折算。这些比例来自全球潮汐模型或地区统计经验对大多数海域精度足够。表 4 是 t_tide 里常见的推断分潮对。推断分潮依赖主分潮适用场景P1K1观测短于 200 天K2S2观测短于 180 天M4M2浅水分潮独立求解不稳定时MS4M2S2浅水区补充开启 inference 后推断分潮在 tidecon 里会有标志区分它们不参与最小二乘而是以固定关系出现在合成结果中。后续用 t_predic 预报时synthesis, 1会把它们包含进去。如果预报结果与实测有系统性偏差优先怀疑推断比例不适合本海区这时可以手动指定比例。手动指定用的是 inferap 和 infername 两个参数% 手动指定 P1 相对 K1 的振幅比和相位差 infer_ap [0.33, -2.5]; % 振幅比和相位差 infer_name {P1; K1}; [tidestruc, xout] t_tide(z, interval, 1, ... start, t0, latitude, lat, ... infer, true, inferap, infer_ap, ... infername, infer_name);手动指定的数值要来自附近长期站点的结果或文献不要随意填。不同海区的 P1/K1 振幅比差异不大但相位差与当地潮波传播路径相关复制别人的数值前先确认海域相近。4.3 排错路径与异常结果判断t_tide 报错集中在三点序列长度不足、有效数据点过少、interval 不匹配。长度不足时报错会提示某个分潮周期比序列还长。有效数据点过少多半是预处理插值没覆盖全部缺口回到预处理阶段检查 NaN 分布即可。interval 不匹配不报错只表现为 xout 与原始信号相位漂移排查方法前面已经说过。异常结果里最隐蔽的是台阶。仪器换电池后水位序列整体抬高或下降几十厘米t_tide 不会报警但残差会出现明显阶跃误差估计显著放大。正式分析前画一张全序列水位图扫一眼有没有水平台阶。有台阶时把数据拆成两段分别分析各出一套调和常数再按天数加权合并。这比强行平移数据更符合潮汐稳定性的实际。另一个经常误判的情况是 xout 与实测差异偏大。残差里有长周期波动属于正常调和分析拟合的是天文潮风暴潮、局地风涌、季节性水位变化都会留在残差里。区分“正常残差”和“分析错误”的方法是看残差里还有没有明显的半日或全日周期如果还能看出潮周期形状说明分潮表有遗漏或参数配置错误如果残差只有几天尺度的起伏则属于非潮汐信号不影响调和常数的可用性。5. 从调和常数做预报与验证的三个实用技巧5.1 t_predic 预测未来水位得到 tidestruc 后预测就是一句话的事t_future datenum(2025, 6, 1) : 1/24 : datenum(2025, 6, 8); z_pred t_predic(t_future, tidestruc, synthesis, 1); plot(t_future, z_pred); datetick(x); grid on;t_predic 的输入是目标时刻序列、tidestruc 和合成开关。synthesis 取 1预报使用全部有效分潮包括推断分潮曲线更平滑取 0 只输出主分潮结果适合做敏感性分析。需要批量预报多个站点时把每个站点的 tidestruc 存成结构数组循环调用即可。5.2 三种验证方式判断调和常数可用性第一种是回代检验用 tidestruc 反演分析时段内的水位计算 RMSE。这个值通常在几十厘米以内大潮期间误差偏大、小潮期间偏小。回代只能检验拟合质量不能代表预报能力。第二种是分割检验把序列前 70% 用于分析后 30% 独立用于预报对比。要求后 30% 至少覆盖一个大潮—小潮周期约 14 天否则反映不出天文潮的半月调制。这种方法最接近实际业务场景推荐优先使用。第三种是跨年检验用相邻年份的独立观测验证。如果调和常数稳定预报误差与回代误差在同一量级如果误差明显变大说明站点附近地形变化剧烈或观测环境不稳定。跨年检验适合连续观测两年以上的站结果可以直接回答“调和常数能不能用于长期预报”。5.3 海流观测的调和分析技巧海流的调和分析与水位流程一致区别是要对 U东向、V北向分别跑 t_tide。输出的振幅单位是 cm/s更适合用椭圆特征参数描述比如最大流速方向、椭圆率、旋转方向这些可以从 U、V 的振幅和迟角换算得到。海流比水位更容易受风驱动影响预处理时建议先低通滤波截止频率取 0.04 周期/小时约 25 小时周期以下能保留全日和半日潮信号同时去掉惯性振荡和大部分风驱高频能量滤波后的序列再做调和分析误差带通常能压缩到原来的三分之一。本文还有配套的精品资源点击获取
返回列表