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

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

三角波（分段线性）插值

原理

对于已知节点 (x₀,y₀)和 (x₁,y₁)，在区间[x₀,x₁]内构造线性函数L(x)以逼近区间内的函数值。
L(x)=y0+y1−y0x1−x0(x−x0)L(x) = y_{0} + \frac{y_{1}-y_{0}}{x_{1}-x_{0}} (x-x_{0}) L(x)=y0+x1−x0y1−y0(x−x0)
在每一个区间内构建诸如上式的方程，即分段构造线性插值函数。

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(xk)=f(xk)p(x_{k}) = f(x_{k}) p(xk)=f(xk)
p_n(x)为三角函数插值多项式，其次数小于等于n。
取一个n阶三角多项式
tk(x)=sinx−x02⋯sinx−xk−12sinx−xk+12⋯sinx−xn2sinxk−x02⋯sinxk−xk−12sinxk−xk+12⋯sinxk−xn2,0≤k≤2nt_{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 tk(x)=sin2xk−x0⋯sin2xk−xk−1sin2xk−xk+1⋯sin2xk−xnsin2x−x0⋯sin2x−xk−1sin2x−xk+1⋯sin2x−xn,0≤k≤2n
可简写成
tk(x)=∑l=0,l≠k2nsinx−xl2sinxk−xl2,0≤k≤2nt_{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 tk(x)=l=0,l=k∑2nsin2xk−xlsin2x−xl,0≤k≤2n
可知，t_k(x)满足
tk(x)=δkll,k=0,1,2…2nt_{k}(x) = \delta _{kl} \space \space l,k=0,1,2 \dots 2n tk(x)=δkl l,k=0,1,2…2n
pn(x)=∑k=02nf(xk)tk(x)p_{n}(x) = \sum_{k=0}^{2n} f(x_{k})t_{k}(x) pn(x)=k=0∑2nf(xk)tk(x)
令x = x_k ，则有p(x_k)=f(x_k)。
由潮汐现象的特征可知，一年高低潮数据大概有1400多个，鉴于潮位所受影响因素太多，如果将所有节点都纳入计算的话，函数复杂程度和计算量太大，拟合精度得不到保障，并且毫无必要。当只选取相邻两个节点进行正弦波插值时，插值效果和精度都较为适宜，故最终n取2。即只选择相邻点为端点，采用一阶三角波插值函数p₁(x)作整个区间的插值函数。

也有，在区间[x_i,x_i+1]上
p1(x)={yi+1−yi2sin[πxi+1−xi(x−xi+1+xi2)]+yi+1+yi2if yi≤yi+1yi−yi+12sin[πxi+1−xi(x−(xi−xi+1−xi2))]+yi+1+yi2if yi>yi+1p_{1} (x)= \begin{cases} {\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{cases} p1(x)={2yi+1−yisin[xi+1−xiπ(x−2xi+1+xi)]+2yi+1+yi2yi−yi+1sin[xi+1−xiπ(x−(xi−2xi+1−xi))]+2yi+1+yiif yi≤yi+1if yi>yi+1
推导过程如下所示。

或将上述两个分段函数表达式整合成如下形式：
p1(x)=abs(yi−yi+1)2sin[πx−xixi+1−xi+(0.5+q)π]+yi+1+yi2p_{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} p1(x)=2abs(yi−yi+1)sin[πxi+1−xix−xi+(0.5+q)π]+2yi+1+yi
式中
q={1if yi≤yi+10if yi>yi+1q = \begin{cases} {1} &\text{if } y_i \le y_{i+1} \\ {0} &\text{if } y_i > y_{i+1} \end{cases} q={10if yi≤yi+1if yi>yi+1

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 ↩︎

潮位分析中高低潮位数据的处理相关推荐

使用ArchR分析单细胞ATAC-seq数据(第四章)
本文首发于我的个人博客, http://xuzhougeng.top/ 往期回顾: 使用ArchR分析单细胞ATAC-seq数据(第一章) 使用ArchR分析单细胞ATAC-seq数据(第二章) 使用 ...
高级转录组分析和R数据可视化第11期(报名线上课还可免费参加线下课2020.6)
福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现决定安排扩增子16S分析.宏基因组.Python课程和转录组的线上直播课.报名参加线上直播课的老师可在1年内选择参加同课程的一次线下课 ...
高级转录组分析和R数据可视化第11期(课程推迟，可先报名，时间另行告知)
封面来源:https://www.zhihu.com/question/304747766 常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后 ...
高级转录组分析和R数据可视化(2020.2,课程推迟,可先报名,时间另行告知)
封面来源:https://www.zhihu.com/question/304747766 常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后 ...
这是入门生信，学习生信分析思路和数据可视化的首选？
封面来源:https://www.zhihu.com/question/304747766 常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后 ...
高级转录组分析和R数据可视化专题研讨会(2019.12)
常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后应该就如做个PCR一样常见.而且分析思路简洁清晰,是入门生信,学习生信分析思路和数据可视化的 ...
转录组分析_高级转录组分析和R数据可视化
封面来源:https://www.zhihu.com/question/304747766 常规转录组是我们最常接触到的一种高通量测序数据类型,其实验方法成熟,花费较低,是大部分CNS必备的技术,以后 ...
计算机网络数据分析报告,贵州大学计算机网络实验报告-实验四-分析IP协议数据包格式...
贵州大学计算机网络实验报告-实验四-分析IP协议数据包格式 (7页) 本资源提供全文预览,点击全文预览即可全文预览,如果喜欢文档就下载吧,查找使用更方便哦! 9.9 积分贵州大学GUIZHOU UN ...
如何在阿里云上使用Data Lake Analytics分析Table Store数据
0. Data Lake Analytics(简称DLA)介绍数据湖(Data Lake)是时下热门的概念,更多阅读可以参考: https://en.wikipedia.org/wiki/Data_ ...

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

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

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

三角波（分段线性）插值

原理

Matlab实现

正弦波插值

原理

Matlab实现

结果对比

潮位分析中高低潮位数据的处理相关推荐

最新文章

热门文章