ADCP潮流调和分析与t_tide椭圆要素计算实战

简介：面向物理海洋学、海洋水文及相关专业学习者与科研人员的潮流调和分析工具包，基于Matlab的t_tide工具包整理，并对原工具作出小幅修正与补充，配有操作说明，可用于潮汐潮流要素计算、调和分析、结果检验及课程实验等场景。压缩包共31个文件，约10.45MB，以m脚本和mat数据为主，辅以xls表格、dat数据、zip子包以及doc、txt说明文档，类型较完整，兼顾程序运行、示例数据与使用指引。已有1685人学习下载，适合需要补足潮流调和分析实操材料、对照代码理解流程的读者与自学者。包内围绕潮流椭圆计算与绘制、潮汐预测、误差检查、调和常数读取等环节提供脚本和配套数据，并保留示例与说明文档，便于复用分析代码、核对计算结果、排查常见问题，从数据读取到结果呈现均可参考，也可作为科研入门或教学演示的资料。

1. 从一段 ADCP 流速序列说起：t_tide 能替你算清哪些潮流要素

手上有一份 ADCP 连续观测的流速记录，时间跨度从十几天到几个月不等，坐标是仪器坐标，采样间隔一小时。把这段时序拆成「余流 + 各分潮」两部分，再算出每个分潮的长轴、短轴、倾角、旋转方向，以及站点属于正规半日潮还是不正规全日潮——这就是潮流要素计算的标准流程，也是物理海洋学和海洋水文方向最常见的作业之一。t_tide 是 MATLAB 下一套成熟的调和分析工具箱，输入一列等间隔的时间序列，输出各分潮的振幅、相位、置信区间和合成序列，再配合一段几十行的椭圆换算脚本，就能把潮流调和常数表和椭圆要素表一次出齐。这份资源包里的程序，正是把「读数据—调和分析—算椭圆—判别潮型—出表」这条链路固化成脚本，适合手上已经有海洋数据获取与分析经验、不想每次从零写分析代码的人。

2. 潮流调和分析的底层账本：分潮表、瑞利准则与自由度

调和分析不是把一条曲线拟合得漂亮就算完，它本质是一次带先验频率的线性最小二乘。先验频率从哪来、你的记录长度能支撑几个频率、拟合出的系数可信到什么程度，这三件事决定了结果能不能拿去写报告。很多人第一次用 t_tide 时被一大堆参数吓住，其实参数背后就是这三件事。

2.1 分潮频率从哪来：t_constituents 表与常用八大分潮

平衡潮理论上，潮位可以展开成一串固定角速度的余弦项，每一项叫一个分潮，频率由天文引数的线性组合决定。t_tide 把这张表放在工具箱目录下的分潮表文件里，包含一百多个分潮的频率、平衡潮振幅和节点因子展开系数。默认调用时它并不会把一百多个全塞进去，而是按你的记录长度和瑞利准则筛出一批可分离的分潮，再逐个做信噪比检验。

实际工作中真正高频出现的是八个：M2、S2、N2、K2、K1、O1、P1、Q1。它们的频率是硬常数，可以直接背下来当检查表用。

拿到一份结果表，第一件事是对频率：M2 的振幅应该是半日潮里最大的，如果你的表里 S2 比 M2 还大，大概率是时间轴错了或者单位混了。这套自查习惯比任何后处理都管用。

2.2 瑞利准则：15 天和 30 天算出的 M2 为什么对不上

两个分潮能不能同时解出来，取决于它们的频率差在你的记录长度内能否分辨。判据就是瑞利准则：

T ≥ R / |f1 − f2|

T 是记录长度（小时），f 是频率（cycles/hour），R 是瑞利系数，t_tide 里由'rayleigh'参数控制，默认 1.0。R 取 1 是最低门槛，取 2 更保守。

把常用分潮对代进去，能得到一张很实用的门槛表：

这段判断可以直接写成脚本，跑数据前先打印一遍，比事后解释「为什么结果怪」省事得多。

% 瑞利准则门槛计算：输入分潮周期（小时），输出可分离所需最短记录长度 per = struct('M2',12.4206012,'S2',12.0,'N2',12.6583475,'K2',11.9672361, ... 'K1',23.9344721,'O1',25.8193387,'P1',24.0658889,'Q1',26.8683569); pairs = {'M2','S2'; 'M2','N2'; 'K1','O1'; 'K1','P1'; 'M2','K2'}; R = 1.0; % 瑞利系数，1 为最低门槛，2 更保守 for k = 1:size(pairs,1) f1 = 1/per.(pairs{k,1}); % cycles/hour f2 = 1/per.(pairs{k,2}); Tmin = R/abs(f1-f2); % 小时 fprintf('%-3s - %-3s : 频率差 %.7f cph, 最短记录 %.1f h (%.1f 天)\n', ... pairs{k,1}, pairs{k,2}, abs(f1-f2), Tmin, Tmin/24); end

per里放的是各分潮周期，取倒数得频率；R改成 2 会把所有门槛翻倍，适合在报告里给出保守估计。这段代码不改任何数据，只用来做可行性预判：记录只有 12 天，就别指望 N2 的振幅能当真；只有 20 天，P1 一定会被 t_tide 自动并进 K1 或者干脆不出现。

还有一层容易被忽略：瑞利准则只回答「能不能分开」，不回答「分得准不准」。序列越长，参与拟合的方程越多，系数的形式误差越小，但前提是这期间潮流本身没有明显非平稳。实测数据里常见的情况是前 20 天 M2 很稳，后 10 天来了一个天气过程，拟合出的 M2 振幅被拉偏。所以长序列不一定更好，分段分析再对比反而更能说明问题。

2.3 信噪比与置信区间：tidecon 四列到底怎么读

t_tide 的核心输出之一是tidecon，一个 N×4 的矩阵，四列依次是振幅、振幅的标准差、相位（度）、相位的标准差。注意这里的「标准差」不是振幅本身的波动范围，而是最小二乘拟合给出的形式误差，它随记录长度增加而减小，随残差方差增大而增大。

选用'error','linear'时，误差由线性化协方差矩阵直接算出，速度快；选用'error','cboot'时走 bootstrap 重采样，结果更稳健，尤其适合残差明显非高斯的流速数据，代价是计算量按重采样次数线性增长。当振幅与误差之比小于 2，也就是信噪比不足时，这个分潮的调和常数不该单独拿出来讨论，比较稳妥的做法是把它归入残差或者用推断分潮的方式约束。

读数的时候我一般按这个顺序过一遍：先看 M2 的振幅误差是否小于振幅的 20%，再看半日潮组的相位是否连续，最后看全日潮组有没有出现振幅接近误差的项。出现过不了的项，回头查数据缺口或者缩短分析窗口，通常能定位到问题。

3. 把数据喂进 t_tide：时间轴、缺测与调用参数

工具箱本身不挑数据来源，但它对输入有三条硬性假设：等间隔、数值型时间轴、时间基准一致。绝大多数「跑出来不对」的案例，问题都出在这三条上，而不是算法本身。

3.1 时间轴的三个硬约束：等间隔、UTC、起点

t_tide 用 MATLAB 的 datenum 作为时间轴，这个数值没有时区概念。观测记录如果是北京时间，直接转 datenum 再喂进去，相位就会整体偏 8 小时对应的角度——半日潮偏 8 小时相当于偏了 232 度，结果看上去「有数」，但和其他站完全对不上。所以读数据的第一步就应该是统一到 UTC，并在脚本里留一行注释写清楚。

% 读入：第 1 列是 'yyyy-MM-dd HH:mm:ss'（UTC），第 2 列是水位或流速 opts = detectImportOptions('wl_sample.csv'); opts.VariableNames = {'tstr','val'}; opts = setvartype(opts, {'datetime','double'}); T = readtable('wl_sample.csv', opts); % 字符串 -> datetime -> datenum，显式声明 UTC，避免受本机时区影响 tdt = datetime(T.tstr, 'InputFormat','yyyy-MM-dd HH:mm:ss', 'TimeZone','UTC'); t = datenum(tdt); % MATLAB datenum，无时区含义 h = T.val(:); % 列向量，t_tide 对行向量不友好 % 等间隔检查 dt_h = median(diff(t))*24; % 采样间隔，小时 assert(abs(dt_h - 1) < 1e-6, '采样间隔不是 1 小时，实测 %.4f', dt_h); % 缺口检查 gap = diff(t)*24; fprintf('采样间隔 %.3f h，最大缺口 %.2f h，缺口数 %d\n', ... dt_h, max(gap), sum(gap > 1.5*dt_h));

TimeZone一定要写，否则datetime会按本机时区解释字符串，换一台机器结果就变。h强制转列向量是血泪经验，t_tide 内部按列做矩阵运算，行向量进去有时不报错但结果错位。缺口的统计不是可选项，后面决定要不要分段就靠这个数。

3.2 缺测、尖峰与去趋势的处理顺序

正确的处理顺序是：先剔除明显尖峰，再补短缺口，最后决定要不要去趋势。顺序反了会把插值出来的假点当成真值再去做尖峰检测，越补越乱。

短缺口（连续缺测不超过 2 到 3 个采样点）用线性插值补上，对半日潮的相位影响可以忽略。长缺口不要硬补，按缺口切成若干连续段，每段单独做调和分析，再比较各段的调和常数是否一致。如果各段差异明显，说明这段观测本身不满足平稳假设，硬拼成一条长序列得到的「平均调和常数」没有物理意义。

尖峰用滑动中值和绝对偏差判据处理，比简单的 3σ 稳，因为流速序列本身有很强的周期性，直接算标准差会被周期项撑大。

% 1) 尖峰剔除：滑动中值 + 局部 MAD 判据 win = 25; % 约一个 M2 周期的窗口 med = movmedian(h, win, 'omitnan'); mad_ = movmad(h, win); % 局部绝对偏差 bad = abs(h - med) > 6 * 1.4826 * mad_; % 1.4826 把 MAD 折算成标准差 h(bad) = NaN; fprintf('剔除尖峰 %d 个（占 %.2f%%）\n', sum(bad), 100*sum(bad)/numel(h)); % 2) 短缺口线性插值，长缺口打标 idx = (1:numel(h))'; h = interp1(idx(~isnan(h)), h(~isnan(h)), idx, 'linear', NaN); run = cumsum(isnan(h)); % 连续缺测段编号，插值后此处仅作示意

实际脚本里我会先用bwlabel或差分累计把连续缺测段标出来，长度超过 3 的段落记下起止时间，分析时按段切分。代码里第二步的插值只对短缺口生效，这一步的判据要放在插值之前算，不能反过来。

关于去趋势：t_tide 本身会拟合一个常数项（相当于平均水位或余流），'secular','linear'会额外拟合线性趋势并从序列中扣除。对水位数据，线性趋势可能来自地面沉降或海平面变化，扣不扣要看你关心什么；对流速数据，线性趋势往往就是季节性环流的一部分，扣掉之后余流就没了。我一般默认用'secular','mean'，只去均值，保留线性项在残差里，事后单独评估。

3.3 主调用与返回值的逐参数说明

把数据整理干净之后，主调用其实只有一行。但每个参数都影响结果，值得逐条说清楚。

lat = 30.5; % 站点纬度，单位度，北纬为正；按实测填，此处仅作占位 dt_h = 1.0; % 采样间隔，小时 t0 = t(1); % 起始时间（datenum），必须显式给出 [tname, tfreq, tcon, hsyn] = t_tide(h, ... 'interval', dt_h, ... % 采样间隔，小时 'start', t0, ... % 相位基准时刻 'latitude', lat, ... % 纬度，用于节点因子修正 'rayleigh', 1.0, ... % 瑞利系数 'error', 'linear', ... % 误差估计方式：linear / cboot 'synthesis', 1, ... % 输出合成序列 'secular', 'mean', ... % 只去均值，不扣线性趋势 'output', 'none'); % 不在命令行打印长表

'interval'与'start'决定相位基准：t_tide 输出的相位是相对t0的滞后角，单位度，含义是A*cos(ω*(t−t0) − g)。换一个起点，相位就整体平移，振幅不变。所以两批数据要对比相位，t0必须一致，或者事后按角速度换算到同一基准。

'latitude'影响节点因子 f 和天文辐角修正 u，主要作用在 K1、O1、M2 这些强分潮上。纬度填 0 相当于关掉大部分节点修正，短序列影响不大，长序列会体现出来，别图省事填 0。

'synthesis', 1让函数多返回一个合成序列hsyn，长度与输入一致，只包含被判定为显著的分潮。这个返回值是后面校验的核心素材，别省这一步。

tname在多数打包版本里是字符矩阵，每行一个分潮名，长度不足处补空格；也有版本返回 cell 数组。写后处理代码时用cellstr(tname)包一层，两种都能吃。

tfreq是各分潮频率（cycles/hour），tcon是前面说的 N×4 矩阵。三者的行序一一对应，做椭圆换算时必须按名字对齐，不能按行号硬对——u 和 v 两次调用筛出的分潮集合经常不一样。

3.4 流速分量：仪器坐标到真北的旋转

潮流椭圆对坐标系极其敏感。ADCP 默认输出的是仪器坐标（beam 或 xyz），如果没转到东-北-天坐标，算出来的倾角就是相对仪器朝向的角度，看上去数值正常，实际上整体偏了一个航向角。这种错误在图上很难发现，因为椭圆形状是对的，只是方向系统性偏了。

常见做法是先转到东分量 u 和北分量 v，再分别做调和分析。如果手上只有仪器坐标和磁罗经朝向，至少要做一次水平旋转：

% 仪器坐标 (x, y) 转真北坐标 (u, v) % heading_deg: 仪器 x 轴相对真北的方位角，顺时针为正 th = deg2rad(heading_deg); u = x*cos(th) + y*sin(th); % 东分量 v = -x*sin(th) + y*cos(th); % 北分量

heading_deg的来源要交代清楚：是磁罗经读数就得先加磁偏角，是 GPS 罗经就直接用。旋转只改方向不改大小，sqrt(u.^2+v.^2)应该和sqrt(x.^2+y.^2)逐点相等，写完拿这个恒等式验一遍，能挡住绝大多数符号错误。

4. 从调和常数到潮流椭圆：四个要素的算法与实现

潮位分析到调和常数就结束了，潮流还多一步：每个分潮在水平面上是一个椭圆，需要把东、北两个分量的振幅和相位合成出长轴、短轴、倾角和旋转方向。这一步没有现成函数，得自己写。

4.1 椭圆要素公式与旋转方向的符号约定

设某个分潮的东分量和北分量分别为

u(t) = Au·cos(ωt − gu)，v(t) = Av·cos(ωt − gv)

在 (u, v) 平面上消去 t，轨迹是一个椭圆。记 Δ = gu − gv，则四个要素为：

• 长半轴 W = sqrt( 0.5·(S + R) )，其中 S = Au² + Av²，R = sqrt( (Au² − Av²)² + 4·Au²·Av²·cos²Δ )
• 短半轴 w = sqrt( 0.5·(S − R) )
• 倾角 θ = 0.5·atan2( 2·Au·Av·cosΔ, Au² − Av² )，从东向逆时针量
• 旋转方向由 sinΔ 的符号决定

这里有个必须自己验一次的坑：不同教材对相位差的定义是 gu − gv 还是 gv − gu，对「逆时针为正」还是「顺时针为正」的约定也不统一。符号一旦搞反，整张表的旋转方向全部颠倒，而数值大小完全正常，肉眼看不出来。

我的做法是用合成序列反算一次：取一个分潮，用它自己的调和常数合成 u、v 各一个周期，画出轨迹，看它是顺时针还是逆时针，再和公式给出的符号对照。跑通一次以后把这个自检固化成脚本，换数据时重跑。

4.2 批量分潮的椭圆计算脚本

下面这段是核心换算，输入是两次t_tide调用的返回值，输出一张按分潮排列的椭圆要素表。

% 输入：uname/ucon 来自 u 分量的 t_tide，vname/vcon 来自 v 分量 % ucon/vcon 均为 N×4：[振幅, 振幅误差, 相位(度), 相位误差(度)] un = cellstr(uname); vn = cellstr(vname); [common, iu, iv] = intersect(un, vn); % 按名字对齐，不能用行号 Au = ucon(iu,1); gu = ucon(iu,3); Av = vcon(iv,1); gv = vcon(iv,3); dphi = deg2rad(gu - gv); % 相位差，弧度 S = Au.^2 + Av.^2; D = Au.^2 - Av.^2; R = sqrt(D.^2 + 4*Au.^2.*Av.^2.*cos(dphi).^2); W = sqrt(0.5*(S + R)); % 长半轴 w = sqrt(0.5*(S - R)); % 短半轴（取正） theta = 0.5*atan2(2*Au.*Av.*cos(dphi), D); % 倾角，弧度 theta_deg = mod(rad2deg(theta), 180); % 折算到 [0,180) rot = sign(sin(dphi)); % +1 逆时针，-1 顺时针 T = table(common, Au, Av, W, w, theta_deg, rot, ... 'VariableNames', {'tide','Au','Av','major','minor','incl_deg','rot'}); disp(T)

intersect那行是关键：u、v 两次调和分析筛出的分潮集合经常不同，比如 u 分量里 N2 过了信噪比检验而 v 分量没过，直接按行号配对就会把 N2 的 u 和 K2 的 v 拼在一起，算出来的椭圆毫无意义。intersect返回公共分潮名和各自的索引，iu/iv就是对齐用的。

theta_deg用mod(...,180)是因为椭圆的长轴方向在 0 到 180 度之间就有定义，超过 180 度等价于减去 180。这一步不做的话，atan2 会返回 −90 到 90 之间的值，直接拿去画图会出现负角度，看着像出错。

rot用sign(sin(dphi))，返回 ±1，画图时映射成「逆时针/顺时针」两个标记。前面说的符号自检就是验这一列。

4.3 潮流类型判别、余流分离与结果落表

椭圆算完，还有两件事要一起做：判潮型和分离余流。

潮型用 Form Factor 判别，公式是

F = (A_O1 + A_K1) / (A_M2 + A_S2)

注意这个判别式用的是潮位振幅，用在流速上时要用同一分潮的长半轴代替振幅，量纲才一致。有的文献用椭圆长轴做 F 因子，结论和用潮位算出来的可能不同，报告里要写清楚用的是哪一种。

余流是原始序列减去所有显著分潮的合成序列之后剩下的平均。t_tide 返回的合成序列hsyn只包含显著分潮，所以余流直接取mean(h - hsyn)。

% 潮型判别（以流速长半轴代替振幅） getA = @(nm) T.major(strcmp(strtrim(T.tide), nm)); % 按名取长半轴 F = (getA('O1') + getA('K1')) / (getA('M2') + getA('S2')); if F < 0.5, type = '正规半日潮'; elseif F < 2.0, type = '不正规半日潮'; elseif F < 4.0, type = '不正规全日潮'; else, type = '正规全日潮'; end fprintf('Form Factor = %.3f，潮型：%s\n', F, type); % 余流：原始序列减去显著分潮合成序列 res_u = u - usyn; res_v = v - vsyn; fprintf('余流 U = %.3f, V = %.3f\n', mean(res_u,'omitnan'), mean(res_v,'omitnan')); % 落表 writetable(T, 'ellipse_result.csv');

getA用匿名函数按名字取长半轴，比手写索引稳，换数据不用改。潮型判别这一段的边界条件写的是<，也就是 F 正好等于 0.5 时归入不正规半日潮，这个约定要和报告里的定义对齐，不同规范略有差异。

res_u/res_v除了取均值算余流，还应该看一下它的标准差——残差标准差过大，说明还有未被模型解释的信号，可能是浅水分潮（M4、M6）或者非潮汐过程，值得单独加项再分析一次。

5. 排查与避坑：五类翻车现场的现象、成因与处置

这一节里的每一条，都是实际跑数据时真出现过的问题，按「看到什么—为什么—怎么办」写，遇到了直接对号入座。

5.1 现象：同一段数据两次跑，相位对不上

原因几乎都是'start'没显式给。t_tide 在不给起点时会用序列的第一个时刻，但如果你的数据是分段读入的，两次调用起点不同，相位基准就不同，振幅一模一样、相位差一大截，看上去像算法不稳定。

处置方式是把起点固定成一个常量写进配置，比如统一用整点日期datenum(2023,1,1,0,0,0)，所有批次都相对这个基准输出相位。事后换算也简单：相位差等于角速度乘以时间差，按-freq*360*dt修正即可。

5.2 现象：结果表里死活找不到 P1，或者 K1 振幅明显偏大

原因是记录长度不够，K1 和 P1 的频率差只有 0.0002281 cph，需要约 183 天才能分离。短序列里 t_tide 会按瑞利准则把 P1 剔除，或者用推断分潮的方式把它并入 K1，导致 K1 的振幅包含了 P1 的贡献。

处置方式有两个方向：如果只需要 K1 的总量，就接受这个合并结果，但在报告里注明；如果需要区分，用'inference'参数指定推断关系，或者干脆把记录截到能分离的长度再分析。不要为了「让 P1 出现」去调小瑞利系数，那只是把不可分辨的两个频率强行分开，结果不可信。

5.3 现象：函数报错或者结果里出现 NaN

原因是序列里有缺测值没有处理。t_tide 内部按矩阵求解，NaN 会污染整个最小二乘，轻则某个分潮的误差变成 NaN，重则直接报错。

处置方式是进函数之前先把 NaN 处理干净：短缺口线性插值，长缺口按段切分。处理完再检查一遍sum(isnan(h))是否为 0，这一步别省。

5.4 现象：椭圆倾角整体偏了一个固定角度

原因是流速没从仪器坐标转到真北坐标，或者旋转时用错了航向角。这个错误最隐蔽的地方在于椭圆形状、长短轴之比、旋转方向全都正确，只有倾角整体平移。

处置方式是回到原始数据，确认 ADCP 输出的坐标系类型，找到航向角来源（磁罗经要加磁偏角），做一次旋转，并用sqrt(u.^2+v.^2)与sqrt(x.^2+y.^2)逐点相等来验证旋转没引入缩放。

5.5 现象：用'error','cboot'跑了几个小时没出结果

原因是 bootstrap 重采样默认次数不低，每增加一个分潮、每增加一倍数据量，耗时都往上走。一段一年的小时序列，在普通笔记本上跑一晚是常事。

处置方式是按需选择：正式出报告前用'linear'快速过一遍，确认分潮选择和量级没问题，最后定稿时再对关键站跑一次cboot。或者用'error','linear'加上分段分析，用多段结果的一致性来代替 bootstrap 的稳定性评估，工程上更划算。

6. 校验与进阶：用预报残差判断调和常数能不能用

调和常数算出来不等于能用。判断依据不是「看起来合理」，而是拿这套常数去预报一段没参与拟合的数据，看预报值和实测值差多少。这一步做了，报告里的结论才站得住。

6.1 三个校验指标：RMS、相关系数、方差解释率

把序列按时间切成两段，前 70% 做拟合，后 30% 做校验。用t_predic按拟合出的调和常数合成后半段，和实测比三个数：

• 残差 RMS：单位与原始数据一致，直接反映绝对误差量级
• 相关系数：反映相位和波形是否对得上
• 方差解释率 = 1 − var(残差)/var(实测)：反映潮流项占了多少总方差

% tq: 校验段的时间（datenum），hq: 校验段实测值 hq_pred = t_predic(tq, tname, tfreq, tcon(:,1:2), 'latitude', lat); res = hq(:) - hq_pred(:); rms_ = sqrt(mean(res.^2, 'omitnan')); cc = corr(hq(:), hq_pred(:), 'rows','complete'); ve = 1 - var(res,'omitnan')/var(hq,'omitnan'); fprintf('RMS=%.3f, R=%.3f, 方差解释率=%.1f%%\n', rms_, cc, 100*ve);

tcon(:,1:2)只取振幅和相位两列，t_predic不需要误差列。方差解释率低于 70% 时，先别怀疑算法，回头查浅水分潮——半日潮占优的站位，M4 和 M6 往往能再解释掉 5 到 10 个百分点，把它们加进分潮表再拟合一次，指标会明显改善。

6.2 节点因子与长序列的 18.6 年修正

节点因子 f 和天文辐角 u 描述的是月球升交点西退对分潮振幅和相位的影响，主周期 18.6 年。短于一年的观测，节点因子的变化很小，做不做修正差别不大；但如果你要做多年尺度对比，或者把不同年份的调和常数放一起比较，就必须统一到同一参考时刻。

t_tide 在给定'latitude'时会自动算节点因子，但参考时刻默认是序列中点。跨年份对比时，我一般显式指定参考时刻，把所有批次的调和常数都归算到同一个历元再比。这一步没有捷径，只能靠脚本里一行注释把参考时刻记死，不然过两个月自己都不记得当时用的是哪个基准。

6.3 批量跑站的脚本骨架

单站跑通之后，剩下的是把它包成循环。骨架大致是这样：读配置表（站号、文件路径、纬度、采样间隔）→ 逐站读数据并清洗 → 调和分析 → 椭圆换算 → 潮型判别 → 落表 → 汇总。每一站的中间结果单独存一个 csv，汇总表只放关键列，出问题能快速定位到站。

跑批量时最容易忽略的是异常捕获。某一站数据缺口太长，函数直接抛错，整个循环中断，前面跑完的结果也白跑。在循环体里包一层try/catch，把失败站号记进日志，跑完统一处理，比逐站守着强得多。

最后说个我自己的习惯。早些年我图快，拿到数据直接进t_tide，出来一张表就往下写，结果有一次把一段只有 11 天的序列当月的资料用，M2 和 S2 根本没分开，两个振幅互相污染，画出来的椭圆长轴偏了将近 20%，直到评审时被问「S2 怎么比 M2 还大」才发现。从那以后我每次跑分析前都强制走一遍三件事：打印瑞利门槛表、检查缺口统计、留 30% 数据做预报校验。三步加起来不到五分钟，但能挡掉绝大多数返工。希望帮到你。

本文还有配套的精品资源，点击获取