潮位分析中高低潮位数据的处理

基于高低潮数据的潮位插值方法

由于在传统的调和分析中，需要输入等间距的潮位数据集，但许多地区的长系列潮位资料仍是以高低潮的形式记录的，即所有数据点为实测高/低潮位，以及对应的出现时间。
故我需要把这样的高低潮数据经过处理（插值法），生成等间距的潮位时间序列。
以下三角波插值和正弦波插值方法参考杨锋的硕士论文《基于高低潮数据的潮汐调和分析方法研究及应用》¹。
所用数据集来自中国南部某潮位测站的实测潮位记录。

三角波（分段线性）插值

原理

对于已知节点 (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}\cdots sin\frac{x-x_{k-1}}{2}sin\frac{x-x_{k+1}}{2}\cdots sin\frac{x-x_n}{2}}{sin\frac{x_k-x_0}{2}\cdots sin\frac{x_k-x_{k-1}}{2}sin\frac{x_k-x_{k+1}}{2}\cdots sin\frac{x_k-x_n}{2}},0\le k\le 2n$$
可简写成
$$t_k(x)=\sum _{l=0,l\ne 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;
       end
    end
    xi = 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;
    else
        y1(i) = 0.5*dfx*sin(pi*(x1(i)-xi+0.5*dx)/dx) + cy;
    end
end

plot(data_time,data_z,'+',x1,y1,'r'), xlabel Time(d), ylabel Height(m);
axis([-0.5 9.5,-1.5,1]);
legend('实测高低潮','正弦波插值','fontsize',14);

结果对比

---

1. http://cdmd.cnki.com.cn/Article/CDMD-10319-1017280624.html ↩︎