这里我想分享一个自己写的潮汐调和分析的MATLAB小程序。当年大三下学期学海洋要素计算时，老师会布置一道大作业就是自己编写潮汐调和分析的FORTRAN程序，记得当时花了好几周写了几百行代码，虽然最后代码跑通了，运行结果看上去也不错，但是觉得自己还是没搞明白。读研后，由于我的研究方向就是潮汐，所以我又重新研究了一遍调和分析方法，我发现其实调和分析基本原理非常简单（其实就是最小二乘，具体可见上面的照片），之所以写程序很麻烦，是因为考虑了太多不必要的东西，比如得自己写程序计算交点因子和交点订正角，甚至得自己编程序求矩阵的逆。我修改了t_tide源程序，将这些非核心的东西剥离，重新写潮汐调和分析的程序（只用了20分钟），代码如下所示，代码也就43行，其中核心代码10行左右。输入参数xin是观测水位，tides是使用哪些分潮，ntides是分潮的个数，dt是观测的时间间隔。输出参数xout是回报的水位，S是平均海平面，H是分潮振幅，G是分潮迟角。

function [S,H,G,coef,xout,ju,coefint]=mytide2(xin,tides,ntides,dt)

%潘海东（中国海洋大学）   Haidong Pan(Ocean University of China) 2017/05/06

%email:1027128116@qq.com

%经典调和分析 classical harmonic analysis

disp(['********* Classical Harmonic Analysis **********  '])

disp(['************* written by Haidong Pan(OUC) ********************  '])

 

load t_constituents

for i=1:ntides

    ju(i)=strmatch(upper(tides(i)),const.name);

    fu(i)=const.freq(ju(i));

end

nobs=length(xin);

gd=find(isfinite(xin(1:nobs)));

ngood=length(gd);

fprintf('   Points used: %d of %d\n',ngood,nobs)

 

nobsu=nobs-rem(nobs-1,2);

t=dt*([1:nobs]'-ceil(nobsu/2));  % 平移到中间时刻

JJ=length(fu);

 

tc=zeros(nobs,1+2*JJ,'double');

tc(:,1)=1;

for i=1:nobs

   for j=1:JJ

      tc(i,1+j)=cos((2*pi)*t(i)*fu(j));

      tc(i,1+JJ+j)=sin((2*pi)*t(i)*fu(j));

   end

end

%coef=tc(gd,:)\xin(gd);

[coef,coefint]=regress(xin(gd),tc(gd,:),0.05);

disp([' Least squares soln:use regress '])

S=coef(1);

 

aa=coef(2:JJ+1);bb=coef(JJ+2:2*JJ+1);

H=sqrt(aa.^2+bb.^2);G=atand(bb./aa);

G(aa<0)=G(aa<0)+180;

xout=tc*coef;

varx=cov(xin(gd));varxp=cov(xout(gd));

disp(['var(x)= ',num2str(varx),'   var(xp)= ',num2str(varxp)]);

disp(['percent var predicted/var original=',num2str(100*varxp/varx),'%']);

disp(['RMSE=',num2str(sqrt(nanmean((xout(gd)-xin(gd)).^2)))]);

disp(['Max error=',num2str(max(abs(xout(gd)-xin(gd))))]);

 

 

大家可以去对比我写的程序和T_TIDE(代码如下，结果几乎完全一样) 。哥伦比亚河口Astoria站的水位数据下载网址（在Water_level_records_Columbia_river. mat文件中）：

[S,H,G,coef,xout,ju,coefint]=mytide2(astoria(1:720,7),{'M2';'S2';'K1';'O1'},4,1);

[NAME,FREQ,TIDECON,XOUT2]=t_tide(astoria(1:720,7),'rayleigh',['M2';'S2';'K1';'O1'] ,'interval',1);

最后简单介绍一下本人开发的新型潮汐调和分析工具包S_TIDE, 可以处理非均匀采样数据，可以处理平稳潮，可以计算多个潮流椭圆参数，欢迎使用，简易教程请见：

工具包下载链接：