潮位分析中高低潮位數(shù)據(jù)的處理

• 基于高低潮數(shù)據(jù)的潮位插值方法
• • 三角波（分段線性）插值
  • • 原理
    • Matlab實(shí)現(xiàn)
  • 正弦波插值
  • • 原理
    • Matlab實(shí)現(xiàn)
  • 結(jié)果對(duì)比

基于高低潮數(shù)據(jù)的潮位插值方法

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

三角波（分段線性）插值

原理

對(duì)于已知節(jié)點(diǎn) (x₀,y₀)和 (x₁,y₁)，在區(qū)間[x₀,x₁]內(nèi)構(gòu)造線性函數(shù)L(x)以逼近區(qū)間內(nèi)的函數(shù)值。
$$L(x)=y_0+\frac{y_1?y_0}{x_1?x_0}(x?x_0)$$
在每一個(gè)區(qū)間內(nèi)構(gòu)建諸如上式的方程，即分段構(gòu)造線性插值函數(shù)。

Matlab實(shí)現(xiàn)

選擇某測站1月份某9天的實(shí)測數(shù)據(jù)進(jìn)行插值。

% data_time & data_z : 已知高低潮位數(shù)據(jù)(時(shí)間和潮位) 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('實(shí)測高低潮','三角波插值','fontsize',14);

正弦波插值

原理

正弦波插值是以三角函數(shù)為插值函數(shù)的一種適用于周期函數(shù)的插值方法。潮位方程是一個(gè)以2k為周期的函數(shù)，令潮位方程為y = f(x)，則x₁,x₂,…,x_n對(duì)應(yīng)函數(shù)值為f(x₁),f(x₂),…,f(x_n)。根據(jù)Lagrange插值原理，存在一個(gè)多項(xiàng)式p_n(x)，滿足
$$p(x_k)=f(x_k)$$
p_n(x)為三角函數(shù)插值多項(xiàng)式，其次數(shù)小于等于n。
取一個(gè)n階三角多項(xiàng)式
$$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)。
由潮汐現(xiàn)象的特征可知，一年高低潮數(shù)據(jù)大概有1400多個(gè)，鑒于潮位所受影響因素太多，如果將所有節(jié)點(diǎn)都納入計(jì)算的話，函數(shù)復(fù)雜程度和計(jì)算量太大，擬合精度得不到保障，并且毫無必要。當(dāng)只選取相鄰兩個(gè)節(jié)點(diǎn)進(jìn)行正弦波插值時(shí)，插值效果和精度都較為適宜，故最終n取2。即只選擇相鄰點(diǎn)為端點(diǎn)，采用一階三角波插值函數(shù)p₁(x)作整個(gè)區(qū)間的插值函數(shù)。

也有，在區(qū)間[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}$$
推導(dǎo)過程如下所示。

或?qū)⑸鲜鰞蓚€(gè)分段函數(shù)表達(dá)式整合成如下形式：
$$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實(shí)現(xiàn)

%正弦波（一階）插值 % data_time & data_z : 已知高低潮位數(shù)據(jù)(時(shí)間和潮位) %三角函數(shù)插值d1 = size(data_time,1);x1 = linspace(0,9,9*24*60+1);% x1為待插值點(diǎn)的x坐標(biāo)（時(shí)間點(diǎn)），數(shù)據(jù)集包含9天的數(shù)據(jù)，故從0~9day，每隔1min生成一個(gè)點(diǎn) d2 = size(x1,2); y1 = zeros(1,d2);for i = 1:d2 flag = 1;for j = 1:d1 %判斷時(shí)間點(diǎn) x1(i) 屬于數(shù)據(jù)集的哪個(gè)區(qū)間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);%求解分段函數(shù)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('實(shí)測高低潮','正弦波插值','fontsize',14);

結(jié)果對(duì)比

---

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

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。