潮位分析中高低潮位數據的處理

• 基于高低潮數據的潮位插值方法
• • 三角波（分段線性）插值
  • • 原理
    • Matlab實現
  • 正弦波插值
  • • 原理
    • Matlab實現
  • 結果對比

基于高低潮數據的潮位插值方法

由于在傳統的調和分析中，需要輸入等間距的潮位數據集，但許多地區的長系列潮位資料仍是以高低潮的形式記錄的，即所有數據點為實測高/低潮位，以及對應的出現時間。
故我需要把這樣的高低潮數據經過處理（插值法），生成等間距的潮位時間序列。
以下三角波插值和正弦波插值方法參考楊鋒的碩士論文《基于高低潮數據的潮汐調和分析方法研究及應用》¹。
所用數據集來自中國南部某潮位測站的實測潮位記錄。

三角波（分段線性）插值

原理

對于已知節點 (x₀,y₀)和 (x₁,y₁)，在區間[x₀,x₁]內構造線性函數L(x)以逼近區間內的函數值。
$$L(x)=y_0+\frac{y_1?y_0}{x_1?x_0}(x?x_0)$$
在每一個區間內構建諸如上式的方程，即分段構造線性插值函數。

Matlab實現

選擇某測站1月份某9天的實測數據進行插值。

% data_time & data_z : 已知高低潮位數據(時間和潮位) x1 = linspace(0,9,9*24*60+1); y1 = interp1(data_time,data_z,x1); %分段線性插值 plot(data_time,data_z,'+',x1,y1), xlabel Time(d), ylabel Height(m); axis([-0.5 9.5,-1.5,1]); legend('實測高低潮','三角波插值','fontsize',14);

正弦波插值

原理

正弦波插值是以三角函數為插值函數的一種適用于周期函數的插值方法。潮位方程是一個以2k為周期的函數，令潮位方程為y = f(x)，則x₁,x₂,…,x_n對應函數值為f(x₁),f(x₂),…,f(x_n)。根據Lagrange插值原理，存在一個多項式p_n(x)，滿足
$$p(x_k)=f(x_k)$$
p_n(x)為三角函數插值多項式，其次數小于等于n。
取一個n階三角多項式
$$t_k(x)=\frac{sin\frac{x?x_0}{2}?sin\frac{x?x_{k?1}}{2}sin\frac{x?x_{k+1}}{2}?sin\frac{x?x_n}{2}}{sin\frac{x_k?x_0}{2}?sin\frac{x_k?x_{k?1}}{2}sin\frac{x_k?x_{k+1}}{2}?sin\frac{x_k?x_n}{2}},0\le k\le 2n$$
可簡寫成
$$t_k(x)=\sum _{l=0,l=k}^{2n}\frac{sin\frac{x?x_l}{2}}{sin\frac{x_k?x_l}{2}},0\le k\le 2n$$
可知，t_k(x)滿足
$$t_k(x)={\delta }_{kl}l,k=0,1,2\dots 2n$$
$$p_n(x)=\sum _{k=0}^{2n}f(x_k)t_k(x)$$
令x = x_k ，則有p(x_k)=f(x_k)。
由潮汐現象的特征可知，一年高低潮數據大概有1400多個，鑒于潮位所受影響因素太多，如果將所有節點都納入計算的話，函數復雜程度和計算量太大，擬合精度得不到保障，并且毫無必要。當只選取相鄰兩個節點進行正弦波插值時，插值效果和精度都較為適宜，故最終n取2。即只選擇相鄰點為端點，采用一階三角波插值函數p₁(x)作整個區間的插值函數。

也有，在區間[x_i,x_i+1]上
$$p_1(x)=\{\begin{array}{ll}\frac{y_{i+1}?y_i}{2}sin[\frac{\pi }{x_{i+1}?x_i}(x?\frac{x_{i+1}+x_i}{2})]+\frac{y_{i+1}+y_i}{2}& \text{if?}{y}_i\le y_{i+1}\\ \frac{y_i?y_{i+1}}{2}sin[\frac{\pi }{x_{i+1}?x_i}(x?(x_i?\frac{x_{i+1}?x_i}{2}))]+\frac{y_{i+1}+y_i}{2}& \text{if?}{y}_i>y_{i+1}\end{array}$$
推導過程如下所示。

或將上述兩個分段函數表達式整合成如下形式：
$$p_1(x)=\frac{abs(y_i?y_{i+1})}{2}sin[\pi \frac{x?x_i}{x_{i+1}?x_i}+(0.5+q)\pi ]+\frac{y_{i+1}+y_i}{2}$$
式中
$$q=\{\begin{array}{ll}1& \text{if?}{y}_i\le y_{i+1}\\ 0& \text{if?}{y}_i>y_{i+1}\end{array}$$

Matlab實現

%正弦波（一階）插值 % data_time & data_z : 已知高低潮位數據(時間和潮位) %三角函數插值d1 = size(data_time,1);x1 = linspace(0,9,9*24*60+1);% x1為待插值點的x坐標（時間點），數據集包含9天的數據，故從0~9day，每隔1min生成一個點 d2 = size(x1,2); y1 = zeros(1,d2);for i = 1:d2 flag = 1;for j = 1:d1 %判斷時間點 x1(i) 屬于數據集的哪個區間if(x1(i)<data_time(j))flag = j;break;endendxi = data_time(flag-1);xif = data_time(flag); %x1(i)屬于[xi , xif]dx = xif - xi;cx = 0.5 * (xi + xif);fxi = data_z(flag-1);fxif = data_z(flag);dfx = abs(fxif - fxi);cy = 0.5 * (fxi + fxif);%求解分段函數if(fxi<=fxif)y1(i) = 0.5*dfx*sin(pi*(x1(i)-cx)/dx) + cy;elsey1(i) = 0.5*dfx*sin(pi*(x1(i)-xi+0.5*dx)/dx) + cy;end endplot(data_time,data_z,'+',x1,y1,'r'), xlabel Time(d), ylabel Height(m); axis([-0.5 9.5,-1.5,1]); legend('實測高低潮','正弦波插值','fontsize',14);

結果對比

---

http://cdmd.cnki.com.cn/Article/CDMD-10319-1017280624.html ??

總結

以上是生活随笔為你收集整理的潮位分析中高低潮位数据的处理的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。