CN113326470B - 遥感水深反演潮高改正方法 - Google Patents

Abstract

本发明提出了遥感水深反演潮高改正方法，属于遥感水深反演领域。该方法包括：(1)根据已有验潮站的潮高数据结合三次样条插值方法，计算卫星过境时刻采集数据时的潮汐高度；(2)根据潮汐改正方法，改正先验水深数据为先验卫星过境时刻水深数据；(3)基于先验卫星过境时刻水深数据结合半经验物理模型，反演卫星过境时刻水深；(4)根据潮高改正方法，改正卫星过境时刻水深为基于大地水准面的遥感反演水深。通过实际实验验证，该方法能提升遥感反演水深的精度，较早期的模型更高的水深反演精度。

Description

遥感水深反演潮高改正方法

技术领域

本发明涉及遥感水深反演领域，特别涉及遥感水深的潮高改正方法。具体而言，涉及基于反演遥感水深过程中水深数据需要进行潮高改正的方法。

背景技术

在海洋水深测量过程中，由于海洋受到太阳和月球的引力影响会产生周期性的潮涨和潮落现象，进而导致海水深度也会随着潮高的变化而变化。为保证水深反演结果归于以大地水准面为基准，需要从遥感反演的海洋水深中消除潮高的影响。

由于无法保证卫星过境时刻外部环境正好满足获取最佳卫星数据的时机，因此遥感反演水深潮高改正与多波束、机载LiDAR测量水深潮高改正截然不同。多波束、机载LiDAR测量水深潮高改正可以通过架设临时验潮站获取准确的潮高数据，应用于潮高改正。然而，由于卫星严重受外部环境影响，很难通过架设临时验潮站实现潮高改正。另外，目前遥感水深反演算法大部分都是半经验物理法，因此需要先验水深数据，而先验水深数据与卫星过境时刻水深数据间有潮高差，需要改正后才能用于水深反演。此外，用卫星过境时刻水深数据结合半经验物理法反演的水深需要消除潮高影响，改正到基于大地水准面的高程才能使用。

为解决上述问题，本发明提出了一种高精度三次样条插值潮高改正方法，可以用于先验水深数据和遥感反演水深的潮高改正，其要求是研究区或其附近有验潮站。

本发明与专利CN102221389A(CN102221389B)结合统计模型与动力模型的乘潮水位预报方法、CN105910587A(CN105910587B)一种基于潮高参数反演的潮汐预测方法、CN201010139189.4一种潮汐预报方法、CN201410741168.8潮汐预报方法、CN201610104433.0一种海上时移地震记录的潮汐校正方法、CN201610255994.0一种基于自适应变异粒子群优化的潮汐智能实时预报方法、CN201910438154.1一种基于FVCOM模型的可视化潮汐潮流预报方法、CN202010187806.1一种通过阶段拟合来预测内河任意点潮汐的方法和CN202010469480.1一种基于时空相关性的潮汐水位预报方法相比，本发明创新性在于充分利用三次样条插值方法和潮高随时间变化的光滑性，针对遥感反演水深提出的高精度三次样条插值潮高改正方法，可以保证误差更小，插值的潮高更精确，极大提升遥感水深反演精度。

本发明与专利CN202010711999.6基于改进GWR模型的多光谱遥感水深反演方法和CN201811623688.3一种基于深度学习的高光谱遥感水深反演方法相比，本发明创新性在于在反演水深的基础上考虑了潮高对水深反演的影响。

发明内容

本发明针遥感反演水深存在的潮高问题，提出了遥感水深反演潮高改正方法，以解决遥感反演水深过程中潮高对水深精度的影响。

遥感水深反演潮高改正方法包括如下步骤：

(1)根据已有验潮站的潮高数据结合三次样条插值方法，计算卫星过境时刻采集数据时的潮高高度。

(2)根据潮高改正方法，改正先验水深数据为先验卫星过境时刻水深数据。

(3)基于先验卫星过境时刻水深数据结合半经验物理模型，反演卫星过境时刻水深。

(4)最后，再根据潮高改正方法，改正卫星过境时刻水深数据为大地水准面的遥感反演水深。

附图说明

图1是潮高改正原理；

图2是实验流程图；

图3是涠洲岛附近水深、先验水深数据点及控制点线布设；

图4是未潮高改正的水深数据反演结果；

图5是潮高改正的水深数据反演结果；

图6是未潮高改正水深数据反演水深、潮高改正水深数据反演水深和真实水深的控制线剖面水深对比；

图7是未潮高改正水深数据反演水深、潮高改正水深数据反演水深和真实水深的直方图对比。

具体实施方式

下面结合附图对本发明的具体实施做进一步详细描述。

结合图1，说明本发明的潮高改正原理及其步骤。

潮汐每天会有高潮面和低潮面，潮高基准面则是距离低潮面D_b的平面，大地水准面则是某一地区的平均海平面，大地水准面到潮高基准面的距离为D_bg，某时刻潮高面到大地水准面距离为D_tg，某时刻潮高面到潮高基准面为D_t，潮高基准面到海底面的距离为D_bs，大地水准面到海地面的距离为D_w，某时刻潮高面到海底面的距离为D_rw，潮高改正就是将以某时刻潮高面起算(卫星过境时刻)的水深D_rw转成以大地水准面起算的水深D_w或将以大地水准面起算的水深D_w转成以某时刻潮高起算(卫星过境时刻)的水深D_rw。

步骤1：大地水准面水深D_w或卫星过境水深D_rw计算

大地水准面水深D_w计算具体步骤如下:

D_w＝-(D_rw-D_tg) (1)

式中，-代表水深是负数，D_rw为以某时刻从潮高面到海底面的距离，D_tg为潮高与大地水准面的距离。

卫星过境水深D_rw计算具体步骤如下:

D_rw＝-(D_w+D_tg) (2)

式中，-代表水深是负数，D_rw为以某时刻从潮高面到海底面的距离，D_tg为潮高与大地水准面的距离。

步骤2：潮高到大地水准面距离D_tg计算

D_tg＝D_t-D_bg (3)

式中，D_t为某时刻潮高，D_bg为大地水准面到潮高基准面的距离。

步骤3：某时刻潮高D_t计算

1、三次样条插值计算潮高

根据潮高数据的一致收敛、一阶可导、二阶可导等特点，潮高数据要满足以下条件方程组。

式中，S(x)是三次样条插值函数，x是以秒为单位的时间戳，x_j是j＝0,1,…,n时在区间[x₀,x_n]内的节点，D_j是D_i(i＝-12,…,-1,1,…,12)从新编号为j＝0,1,…,n后相应的潮高。S(x_j-0)是S(x)在节点x_j处的左极限，S(x_j+0)是S(x)在节点x_j处的右极限，S’(x_j-0)是S(x)在节点x_j处的一阶左导数，S’(x_j+0)是S(x)在节点x_j处的一阶右导数，S”(x_j-0)是S(x)在节点x_j处的二阶左导数，S”(x_j+0)是S(x)在节点x_j处的二阶右导数，S”(x₀)是S(x)在节点x₀处的二阶导数，D”(x₀)是D(x)在节点x₀处的二阶导数，S”(x_n)是S(x)在节点x_n处的二阶导数，D”(x_n)是D(x)在节点x_n处的二阶导数。

条件方程组构造包括三次样条函数S(x)及其一阶导函数S’(x)、二阶导函数S”(x)的构造。

由于S(x)在区间[x_j,x_j+1]上是三次多项式，故S”(x)在[x_j,x_j+1]上是线性函数，根据拉格朗日插值法可表示为

式中，S”(x)是S(x)的二阶导函数，M_j是j＝0,1,…,n时的参数，h_j＝x_j+1-x_j。

对公式(5)积分两次并利用S”(x_j)＝D_j及S”(x_j+1)＝D_j+1，可以确定积分常数，于是得三次样条表达式

式中，S(x)是三次样条插值函数，M_j是j＝0,1,…,n-1时的参数，h_j＝x_j+1-x_j(j＝0,1,…,n-1)。

为了确定S(x)，公式(6)中M_j是需要确定，因此对S(x)求导可得

式中，S’(x)是S(x)的导函数，M_j是j＝0,1,…,n时的参数(未知)，h_j＝x_j+1-x_j。

据此，可知S(x)在区间[x_j-1,x_j]上的一阶左导数和一阶右导数为

根据潮高曲线特征可知S(x)在区间[x_j-1,x_j]上的一阶导数连续，可知

S'(x_j+0)＝S'(x_j-0) (9)

将公式(8)代入公式(9)中得

μ_jM_j-1+2M_j+λ_jM_j+1＝d_j (10)

式中，

条件方程组参数计算，主要是将初始参数x_j和D_i、公式(5)、(6)、(7)、(8)代入构造好的条件方程组(4)，且令λ₀＝μ_n＝0,d₀＝2D”₀,d_n＝2D”_n，可得S(x)系数M_j的矩阵

由已知条件可知，条件方程组(4)的系数矩阵元素λ_j,μ_j已确定，且满足下式

因此，系数矩阵是对角占优矩阵，从而条件方程组(4)有唯一解。

方程组(4)采用追赶法求解，具体如下。令

式中，α_i,β_i,γ_i为待定系数。由公式(13)得

式中，b₁＝α₁≠0，|b₁|＞|c₁|＞0,β₁＝c₁/b₁，0＜|β_i|＜1。公式(14)可写成矩阵方程形式为：

Am＝D (15)

而求解矩阵方程组(15)等价于求解以下的三角方程组

求得m，从而得到解三对角线方程组的追赶公式。

①计算的递推公式

②解Ly＝D

③解Um＝y

2、D_i(i＝-12,…,-1,1,…,12)潮高计算

(1)获取卫星过境时间戳前后各12个小时整点潮高数据

获取卫星过境时间戳前后各12个小时整点潮高数据公式

D_i＝f(x_i) (20)

式中，D_i为(i＝-12,…,-1,1,…,12)整点潮高，f是时间戳的潮高函数，x_i为(i＝-12,…,-1,1,…,12)整点时间戳，x₀为卫星过境时刻的时间戳，x_-1是卫星过境时刻前1小时内的整点时间戳，x₁为卫星过境时刻后1小时内整点时间戳，其他以此类推。

(2)卫星过境时刻时间戳x计算方法见公式

x＝((X*24+h)*60+m)*60+s (21)

式中，x是卫星过境时刻时间戳，X是从1970年1月1日至卫星过境时刻的天数，h是从1970年1月1日至卫星过境时刻减去整数天的小时数，m是从1970年1月1日至卫星过境时刻减去整数天和整数小时的分钟，s是从1970年1月1日至卫星过境时刻减去整数天、整数小时和整数分钟数的秒。

X的计算方法见公式

X＝Y+M+D-719162 (22)

X是从1970年1月1日至卫星过境时刻的天数，Y是卫星过境时刻的年份，M是卫星过境时刻的月份，D是卫星过境时刻的当月第几天。

Y、M和D的计算方法见公式(23)

实验：

使用涠洲岛西边的水深、2020年2月23日03:11:03获得的Landsat8 Oil无云卫星影像数据，结合Strumpf的对数比值水深反演方法，对遥感反演水深潮高改正方法的效果进行评价。按照图2的实验流程，具体实验步骤如下。

步骤1：选择涠洲岛附近研究区，布设7个先验水深数据点、7个控制点、4条控制线，见图3。同时对遥感数据完成辐射定标、大气校正以及研究区裁剪。

步骤2：将布设的7个先验水深数据点用公式(2)改正潮高，获取7个先验数据点卫星过境时刻水深，见表1。

表1选择的7个先验水深数据点及其卫星过境时刻水深

步骤3：将7个先验水深数据点用于Strumpf的对数比值水深反演方法，反演结果为未潮高改正的水深数据反演结果，见图4；将7个先验数据点卫星过境时刻水深用于Strumpf的对数比值水深反演方法，反演研究区卫星过境时刻水深。

步骤4：将反演的研究区卫星过境时刻水深用公式(1)改正潮高，获得以大地水准面(1985高程基准面)为起算面的水深数据，结果为潮高改正的水深数据反演结果，见图5。

步骤5：对比7个控制点未潮高改正的水深数据反演水深、潮高改正的水深数据反演水深与真实水深见表2；对比4条控制线未潮高改正的水深、潮高改正水深与真实水深的剖面见图6；对比未潮高改正的水深数据反演水深、潮高改正的水深数据反演水深与真实水深直方图见图7。

表2选择的7个控制点的未潮高改正的水深数据反演水深、潮高改正的水深数据反演水深与真实水深

从表2可知控制点1、2未潮高改正的水深数据反演水深更接近真实水深，而其他控制点中潮高改正的水深数据反演水深更接近真实水深是因为未潮高改正的水深数据反演水深集中在13-15m内，参见图7。

综上所述，遥感水深反演潮高改正方法在一定程度上能够提高水深反演的精度，使反演的水深更接近真实水深，也更有实际意义和价值。

尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包括这些改动和变型在内。

Claims (1)

1.遥感水深反演潮高改正方法，其特征在于包括如下步骤：

(1)根据已有验潮站的潮高数据结合三次样条插值方法，计算卫星过境时刻采集数据时的潮汐高度；

(2)根据潮高改正方法，改正先验水深数据为先验卫星过境时刻水深数据；

(3)基于先验卫星过境时刻水深数据结合半经验物理模型，反演卫星过境时刻水深；

(4)最后，再根据潮高改正方法，改正卫星过境时刻水深数据为大地水准面的遥感反演水深；

其中，潮高改正原理及其步骤如下：潮汐每天会有高潮面和低潮面，潮高基准面则是距离低潮面D_b的平面，大地水准面则是某一地区的平均海平面，大地水准面到潮高基准面的距离为D_bg，卫星过境时刻潮高面到大地水准面距离为D_tg，卫星过境时刻潮高面到潮高基准面为D_t，潮高基准面到海底面的距离为D_bs，大地水准面到海底面的距离为D_w，卫星过境时刻潮高面到海底面的距离为D_rw，潮高改正就是将D_rw转成D_w或将D_w转成D_rw；

步骤1：大地水准面到海底面的距离D_w或卫星过境时刻潮高面到海底面的距离D_rw计算

大地水准面到海底面的距离D_w计算具体步骤如下:

D_w＝-(D_rw-D_tg) (1)

式中，-代表水深是负数；

卫星过境时刻潮高面到海底面的距离D_rw计算具体步骤如下:

D_rw＝-(D_w+D_tg) (2)

步骤2：卫星过境时刻潮高面到大地水准面距离D_tg计算

D_tg＝D_t-D_bg (3)

步骤3：卫星过境时刻潮高面到潮高基准面D_t计算

第一、三次样条插值计算潮高

根据潮高数据的一致收敛、一阶可导、二阶可导特点，潮高数据要满足以下条件方程组：

式中，S(x)是三次样条插值函数，x是以秒为单位的时间戳，x_j是j＝0,1,…,n时在区间[x₀,x_n]内的节点，D_j是D_i，i＝-12,…,-1,1,…,12重新编号为j＝0,1,…,n后相应的潮高；S(x_j-0)是S(x)在节点x_j处的左极限，S(x_j+0)是S(x)在节点x_j处的右极限，S’(x_j-0)是S(x)在节点x_j处的一阶左导数，S’(x_j+0)是S(x)在节点x_j处的一阶右导数，S”(x_j-0)是S(x)在节点x_j处的二阶左导数，S”(x_j+0)是S(x)在节点x_j处的二阶右导数，S”(x₀)是S(x)在节点x₀处的二阶导数，D”(x₀)是D(x)在节点x₀处的二阶导数，S”(x_n)是S(x)在节点x_n处的二阶导数，D”(x_n)是D(x)在节点x_n处的二阶导数；

条件方程组构造包括三次样条插值函数S(x)及其一阶导函数S’(x)、二阶导函数S”(x)的构造：

由于S(x)在区间[x_j,x_j+1]上是三次多项式，故S”(x)在[x_j,x_j+1]上是线性函数，根据拉格朗日插值法可表示为

式中，S”(x)是S(x)的二阶导函数，M_j是j＝0,1,…,n时的参数，h_j＝x_j+1-x_j；

对公式(5)积分两次并利用S”(x_j)＝D_j及S”(x_j+1)＝D_j+1，确定积分常数，于是得三次样条插值函数表达式

式中，S(x)是三次样条插值函数，M_j是j＝0,1,…,n-1时的参数，h_j＝x_j+1-x_j；

为了确定S(x)，公式(6)中M_j是需要确定，因此对S(x)求导可得

式中，S’(x)是S(x)的一阶导函数，M_j是j＝0,1,…,n时的参数，h_j＝x_j+1-x_j；

据此，可知S(x)在区间[x_j-1,x_j]上的一阶左导数和一阶右导数为

根据潮高曲线特征可知S(x)在区间[x_j-1,x_j]上的一阶导数连续，可知

S'(x_j+0)＝S'(x_j-0) (9)

将公式(8)代入公式(9)中得

μ_jM_j-1+2M_j+λ_jM_j+1＝d_j (10)

式中，

条件方程组参数计算，主要是将初始参数x_j和D_i、公式(5)、(6)、(7)、(8)代入构造好的条件方程组(4)，且令λ₀＝μ_n＝0,d₀＝2D'₀',d_n＝2D'_n'，可得S(x)参数M_j的矩阵

由已知条件可知，条件方程组(4)的系数矩阵元素λ_j,μ_j已确定，且满足下式

因此，系数矩阵是对角占优矩阵，从而条件方程组(4)有唯一解；方程组(4)采用追赶法求解，具体如下：

令

式中，α_i,β_i,γ_i为待定系数；由公式(13)得

式中，b₁＝α₁≠0，|b₁|＞|c₁|＞0,β₁＝c₁/b₁，0＜|β_i|＜1；

公式(14)可写成矩阵方程形式为：

Am＝D (15)

而求解矩阵方程组(15)等价于求解以下的三角方程组

求得m，从而得到解三对角线方程组的追赶公式：

①计算的递推公式

②解Ly＝D

③解Um＝y

第二、D_i潮高计算

(1)获取卫星过境时间戳前后各12个小时整点潮高数据

获取卫星过境时间戳前后各12个小时整点潮高数据公式

D_i＝f(x_i) (20)

式中，D_i为i＝-12,…,-1,1,…,12整点潮高，f是时间戳的潮高函数，x_i为i＝-12,…,-1,1,…,12整点时间戳，x_-1是卫星过境时刻前1小时内的整点时间戳，x₁为卫星过境时刻后1小时内整点时间戳，其他以此类推；

(2)卫星过境时刻时间戳x计算方法见公式

x＝((X*24+h)*60+m)*60+s (21)

式中，x是卫星过境时刻时间戳，X是从1970年1月1日至卫星过境时刻的天数，h是从1970年1月1日至卫星过境时刻减去整数天的小时数，m是从1970年1月1日至卫星过境时刻减去整数天和整数小时的分钟，s是从1970年1月1日至卫星过境时刻减去整数天、整数小时和整数分钟数的秒；

X的计算方法见公式

X＝Y+M+D-719162 (22)

Y是卫星过境时刻的年份，M是卫星过境时刻的月份，D是卫星过境时刻的当月第几天；

Y、M和D的计算方法见公式(23)

Images