潮汐觀測資料調和分析及預報成圖

首先，需要說明兩點，
第一，本案例使用的是單月潮汐觀測資料，處理方法則是基于長期觀測資料的調和分析來進行處理，中期觀測資料的分析需要分別求主要分潮、隨從分潮，短期觀測資料分析則還需要計算不同觀測序列的權重，但是核心演算法與長期觀測資料分析是一致的，都是建立矛盾方程，然后使用最小二乘法建立法方程，求出法方程系數，再求出矩陣X、Y，
第二，本文主要介紹方法步驟，所用代碼大多為關鍵步驟實作，僅供參考，如需完整代碼，請關注博主的另一篇資源，
下面開始介紹本案例的處理，從理論上講，首先，我們需要選取分潮，確立我們所要分析的天文分潮，本案例用8個主要分潮——M2、S2、N2、K2、K1、O1、P1、Q1，四個半日潮，四個全日潮（其實去看潮汐相關研究的文獻就會發現，基本都是以這八個分潮為主的），相鄰資料時間間隔為一小時，以所有資料的中間數對應時刻作為時間原點，然后對觀測記錄資料進行排序，
從實際出發，在matlab中，首先清理空間、準備環境（這是一個良好習慣），然后需要匯入資料，對分潮進行排序，

%%匯入資料
clear all;close all;clc;
data=importdata('C:\Users\STAR\Desktop\TideData_01.txt');

%%分潮排序（八個）
M2=1;
S2=2;
N2=3;
K2=4;
K1=5;
O1=6;
P1=7;
Q1=8;

常量的準備，輸入杜德森數，以及天文元素隨時間的變化速度，

%%杜德森數（八個分潮各有七個杜德森數）
miu{M2}=[2,0,0,0,0,0,0];
miu{S2}=[2,2,-2,0,0,0,0];
miu{N2}=[2,-1,0,1,0,0,0];
miu{K2}=[2,2,0,0,0,0,0];
miu{K1}=[1,1,0,0,0,0,1];
miu{O1}=[1,-1,0,0,0,0,-1];
miu{P1}=[1,1,-2,0,0,0,-1];
miu{Q1}=[1,-2,0,1,0,0,-1];

%%天文元素隨時間的變化速度
rateOfChange=[14.49205211,0.54901653,0.04106864,0.00464183,0.00220641,0.00000196];

然后需要計算時間原點，個人習慣將代碼模塊化，本處使用自編函式TimeCalculation，

year=2003;
month=3;
day=1;
hour=0;
[year,month,day,hour]=TimeCalculation(2003,3,1,0,360);

計算分潮角速度

for i=M2:Q1
   sigma(i)=AngularVelocity(miu{i},rateOfChange);
end
sigma=deg2rad(sigma);   %注意，這里涉及到一個角度轉弧度的操作

然后是計算各天文元素

[tao,ss,hhh,pp,NNN,ppp]=AstronomicalElements(year,month,day,hour);
tao=deg2rad(tao);
ss=deg2rad(ss);
hhh=deg2rad(hhh);
pp=deg2rad(pp);
NNN=deg2rad(NNN);
ppp=deg2rad(ppp);
astronomicalElements=[tao,ss,hhh,pp,NNN,ppp];

計算分潮初始位相，同樣使用自編函式

for i=M2:Q1
   v0(i)=InitialPhase(miu{i},astronomicalElements);
end
v0(2) = 6.2832;%由于S2分潮的特殊性，直接賦值

計算交點因子及訂正角

deltaMiu4{M2}=[0,0,0,2,2];
deltaMiu5{M2}=[-2,-1,0,0,1];
rho{M2}=[0.0005,-0.0373,1,0.0006,0.0002];
deltaMiu4{K2}=[0,0,0,0];
deltaMiu5{K2}=[-1,0,1,2];
rho{K2}=[-0.0128,1,0.2980,0.0324];
deltaMiu4{K1}=[-2,0,0,0,0,0];
deltaMiu5{K1}=[-1,-2,-1,0,1,2];
rho{K1}=[0.0002,0.0001,-0.0198,1,0.1356,-0.0029];
deltaMiu4{P1}=[0,0,0,2,2];
deltaMiu5{P1}=[-2,-1,0,0,1];
rho{P1}=[0.0008,-0.0112,1,-0.0015,-0.0003];
deltaMiu4{O1}=[0,0,0,2,2,2];
deltaMiu5{O1}=[-2,-1,0,-1,0,1];
rho{O1}=[-0.0058,0.1885,1,0.0002,-0.0064,-0.0010];
for i=[M2,K2,K1,P1,O1]
    [f(i),u(i)]=IntersectionFactorAndCorrectionAngle(deltaMiu4{i},deltaMiu5{i},rho{i},pp,NNN);
end
f(Q1)=f(O1);
u(Q1)=u(O1);
f(N2)=f(M2);
u(N2)=u(M2);
f(S2)=1;
u(S2)=0;

建立法方程并求x和y

deltaT=1;
N=data(end,1);
NN=(N-1)/2;
A(0+1,0+1)=N;
for i=M2:Q1
    A(0+1,i+1)=(sin(N*sigma(i)*deltaT/2))/(sin(sigma(i)*deltaT/2));
    A(i+1,0+1)=(sin(N*sigma(i)*deltaT/2))/(sin(sigma(i)*deltaT/2));
    A(i+1,i+1)=(N+(sin(N*sigma(i)*deltaT)/sin(sigma(i)*deltaT)))/2;
    B(i,i)=(N-(sin(N*sigma(i)*deltaT)/sin(sigma(i)*deltaT)))/2;
end
for i=M2:Q1
   for j=M2:Q1
      if ~(i==j)
           A(i+1,j+1)=(((sin(N/2*(sigma(i)-sigma(j))*deltaT))/(sin(1/2*(sigma(i)-sigma(j))*deltaT)))+(sin(N/2*(sigma(i)+sigma(j))*deltaT))/(sin(1/2*(sigma(i)+sigma(j))*deltaT)))/2 ;
           B(i,j)=(((sin(N/2*(sigma(i)-sigma(j))*deltaT))/(sin(1/2*(sigma(i)-sigma(j))*deltaT)))-(sin(N/2*(sigma(i)+sigma(j))*deltaT))/(sin(1/2*(sigma(i)+sigma(j))*deltaT)))/2 ;
       end 
   end
end
F1(0+1)=sum(data(:,2));
for i=M2:Q1
    F1(i+1)=sum(data(:,2).*cos((data(:,1)-361)*sigma(i)*deltaT));
    F2(i)=sum(data(:,2).*sin((data(:,1)-361)*sigma(i)*deltaT));
end

X=F1/A(:,1:end);
Y=F2/B(:,1:end);

計算準調和振幅R和位相theta

R=(X(2:end).^2+Y.^2).^0.5;
for i=M2:Q1
    theta(i)=CalculatedPhase(R(i),X(i+1),Y(i));
end

計算分潮的調和常數

H=R./f;
g=theta+v0+u;
for i=M2:Q1
   g(i)=rad2pi(g(i));
end

潮汐預報、計算自報余差

S0=X(0+1);
for i=1:721
    h(i)=S0+sum(f.*H.*cos(sigma*(i-361)+v0+u-g));
end
r=data(1:end,2)'-h(1:end);
delta=sum(r.^2)^0.5/data(end,1);

計算預報潮位

starttime=(31-1+30+1)*24+1;
endtime=starttime+31*24+1;
forecastTime=1:31*24+1+1;
for i=starttime:endtime
    forecastTide(i-starttime+1)=S0+sum(f.*H.*cos(sigma*(i-361)+v0+u-g));
end

繪制預報圖

figure(2);
plot(forecastTime,forecastTide);
title("Forecast(May 1 to June 1)");
xlabel("Serial number");
ylabel("Stage");
legend("Forecast",'Location','Best');

寫在文末，潮汐體現的是海洋動力及水文要素的變化規律和控制機制，從古至今，關于潮汐的研究從未停止，遠有沈括著作《夢溪筆談》，近有牛頓平衡潮理論之廣泛應用，希望本文能為有需要的人提供一點思路及方法，