新浪博客
，运行结果看上去也不错，但是觉得自己还是没搞明白。读研后，由于我的研究方向就是潮汐，所以我又重新研究了一遍调和分析方法，我发现其实调和分析基本原理非常简单（其实就是最小二乘，具体可见上面的照片），之所以写程序很麻烦，是因为考虑了太多不必要的东西，比如得自己写程序计算交点因子和交点订正角，甚至得自己编程序求矩阵的逆。我修改了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)
%潘海东（中国海洋大学）
%email:1027128116@qq.com
%经典调和分析 classical harmonic analysis
disp(['********* Classical Harmonic Analysis **********
disp(['************* written by Haidong Pan(OUC) ********************
load t_constituents
for i=1:ntides
end
nobs=length(xin);
gd=find(isfinite(xin(1:nobs)));
ngood=length(gd);
fprintf('
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
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),'
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文件中）：
https://www.researchgate.net/publication/323880017_Water_level_records_Columbia_river
[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, 可以处理非均匀采样数据，可以处理平稳潮，可以计算多个潮流椭圆参数，欢迎使用，简易教程请见：
https://www.bilibili.com/read/cv7785513
工具包下载链接：
https://www.researchgate.net/project/A-non-stationary-tidal-analysis-toolbox-S-TIDE
function [S,H,G,coef,xout,ju,coefint]=mytide2(xin,tides,ntides,dt)
%潘海东（中国海洋大学）
%email:1027128116@qq.com
%经典调和分析 classical harmonic analysis
disp(['********* Classical Harmonic Analysis **********
disp(['************* written by Haidong Pan(OUC) ********************
load t_constituents
for i=1:ntides
end
nobs=length(xin);
gd=find(isfinite(xin(1:nobs)));
ngood=length(gd);
fprintf('
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
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),'
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文件中）：
https://www.researchgate.net/publication/323880017_Water_level_records_Columbia_river
[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, 可以处理非均匀采样数据，可以处理平稳潮，可以计算多个潮流椭圆参数，欢迎使用，简易教程请见：
https://www.bilibili.com/read/cv7785513
工具包下载链接：
https://www.researchgate.net/project/A-non-stationary-tidal-analysis-toolbox-S-TIDE