[0020] 下面结合附图对本发明的具体实施做进一步详细描述。
[0021] 结合图1，说明本发明的潮高改正原理及其步骤。
[0022] 潮汐每天会有高潮面和低潮面，潮高基准面则是距离低潮面Db的平面，大地水准面则是某一地区的平均海平面，大地水准面到潮高基准面的距离为Dbg，某时刻潮高面到大地水准面距离为Dtg，某时刻潮高面到潮高基准面为Dt，潮高基准面到海底面的距离为Dbs，大地水准面到海地面的距离为Dw，某时刻潮高面到海底面的距离为Drw，潮高改正就是将以某时刻潮高面起算(卫星过境时刻)的水深Drw转成以大地水准面起算的水深Dw或将以大地水准面起算的水深Dw转成以某时刻潮高起算(卫星过境时刻)的水深Drw。
[0023] 步骤1：大地水准面水深Dw或卫星过境水深Drw计算
[0024] 大地水准面水深Dw计算具体步骤如下:
[0025] Dw＝‑(Drw‑Dtg) (1)
[0026] 式中，‑代表水深是负数，Drw为以某时刻从潮高面到海底面的距离，Dtg为潮高与大地水准面的距离。
[0027] 卫星过境水深Drw计算具体步骤如下:
[0028] Drw＝‑(Dw+Dtg) (2)
[0029] 式中，‑代表水深是负数，Drw为以某时刻从潮高面到海底面的距离，Dtg为潮高与大地水准面的距离。
[0030] 步骤2：潮高到大地水准面距离Dtg计算
[0031] Dtg＝Dt‑Dbg (3)
[0032] 式中，Dt为某时刻潮高，Dbg为大地水准面到潮高基准面的距离。
[0033] 步骤3：某时刻潮高Dt计算
[0034] 1、三次样条插值计算潮高
[0035] 根据潮高数据的一致收敛、一阶可导、二阶可导等特点，潮高数据要满足以下条件方程组。
[0036]
[0037] 式中，S(x)是三次样条插值函数，x是以秒为单位的时间戳，xj是j＝0,1,…,n时在区间[x0,xn]内的节点，Dj是Di(i＝‑12,…,‑1,1,…,12)从新编号为j＝0,1,…,n后相应的潮高。S(xj‑0)是S(x)在节点xj处的左极限，S(xj+0)是S(x)在节点xj处的右极限，S’(xj‑0)是S(x)在节点xj处的一阶左导数，S’(xj+0)是S(x)在节点xj处的一阶右导数，S”(xj‑0)是S(x)在节点xj处的二阶左导数，S”(xj+0)是S(x)在节点xj处的二阶右导数，S”(x0)是S(x)在节点x0处的二阶导数，D”(x0)是D(x)在节点x0处的二阶导数，S”(xn)是S(x)在节点xn处的二阶导数，D”(xn)是D(x)在节点xn处的二阶导数。
[0038] 条件方程组构造包括三次样条函数S(x)及其一阶导函数S’(x)、二阶导函数S”(x)的构造。
[0039] 由于S(x)在区间[xj,xj+1]上是三次多项式，故S”(x)在[xj,xj+1]上是线性函数，根据拉格朗日插值法可表示为
[0040]
[0041] 式中，S”(x)是S(x)的二阶导函数，Mj是j＝0,1,…,n时的参数，hj＝xj+1‑xj。
[0042] 对公式(5)积分两次并利用S”(xj)＝Dj及S”(xj+1)＝Dj+1，可以确定积分常数，于是得三次样条表达式
[0043]
[0044] 式中，S(x)是三次样条插值函数，Mj是j＝0,1,…,n‑1时的参数，hj＝xj+1‑xj(j＝0,1,…,n‑1)。
[0045] 为了确定S(x)，公式(6)中Mj是需要确定，因此对S(x)求导可得
[0046]
[0047] 式中，S’(x)是S(x)的导函数，Mj是j＝0,1,…,n时的参数(未知)，hj＝xj+1‑xj。
[0048] 据此，可知S(x)在区间[xj‑1,xj]上的一阶左导数和一阶右导数为
[0049]
[0050] 根据潮高曲线特征可知S(x)在区间[xj‑1,xj]上的一阶导数连续，可知[0051] S'(xj+0)＝S'(xj‑0) (9)
[0052] 将公式(8)代入公式(9)中得
[0053] μjMj‑1+2Mj+λjMj+1＝dj (10)
[0054] 式中，
[0055]
[0056] 条件方程组参数计算，主要是将初始参数xj和Di、公式(5)、(6)、(7)、(8)代入构造好的条件方程组(4)，且令λ0＝μn＝0,d0＝2D”0,dn＝2D”n，可得S(x)系数Mj的矩阵[0057]
[0058] 由已知条件可知，条件方程组(4)的系数矩阵元素λj,μj已确定，且满足下式[0059]
[0060] 因此，系数矩阵是对角占优矩阵，从而条件方程组(4)有唯一解。
[0061] 方程组(4)采用追赶法求解，具体如下。令
[0062]
[0063] 式中，αi,βi,γi为待定系数。由公式(13)得
[0064]
[0065] 式中，b1＝α1≠0，|b1|＞|c1|＞0,β1＝c1/b1，0＜|βi|＜1。公式(14)可写成矩阵方程形式为：
[0066] Am＝D (15)
[0067] 而求解矩阵方程组(15)等价于求解以下的三角方程组
[0068]
[0069] 求得m，从而得到解三对角线方程组的追赶公式。
[0070] ①计算的递推公式
[0071]
[0072] ②解Ly＝D
[0073]
[0074] ③解Um＝y
[0075]
[0076] 2、Di(i＝‑12,…,‑1,1,…,12)潮高计算
[0077] (1)获取卫星过境时间戳前后各12个小时整点潮高数据
[0078] 获取卫星过境时间戳前后各12个小时整点潮高数据公式
[0079] Di＝f(xi) (20)
[0080] 式中，Di为(i＝‑12,…,‑1,1,…,12)整点潮高，f是时间戳的潮高函数，xi为(i＝‑12,…,‑1,1,…,12)整点时间戳，x0为卫星过境时刻的时间戳，x‑1是卫星过境时刻前1小时内的整点时间戳，x1为卫星过境时刻后1小时内整点时间戳，其他以此类推。
[0081] (2)卫星过境时刻时间戳x计算方法见公式
[0082] x＝((X*24+h)*60+m)*60+s (21)
[0083] 式中，x是卫星过境时刻时间戳，X是从1970年1月1日至卫星过境时刻的天数，h是从1970年1月1日至卫星过境时刻减去整数天的小时数，m是从1970年1月1日至卫星过境时刻减去整数天和整数小时的分钟，s是从1970年1月1日至卫星过境时刻减去整数天、整数小时和整数分钟数的秒。
[0084] X的计算方法见公式
[0085] X＝Y+M+D‑719162 (22)
[0086] X是从1970年1月1日至卫星过境时刻的天数，Y是卫星过境时刻的年份，M是卫星过境时刻的月份，D是卫星过境时刻的当月第几天。
[0087] Y、M和D的计算方法见公式(23)
[0088]
[0089] 实验：
[0090] 使用涠洲岛西边的水深、2020年2月23日03:11:03获得的Landsat8 Oil无云卫星影像数据，结合Strumpf的对数比值水深反演方法，对遥感反演水深潮高改正方法的效果进行评价。按照图2的实验流程，具体实验步骤如下。
[0091] 步骤1：选择涠洲岛附近研究区，布设7个先验水深数据点、7个控制点、4条控制线，见图3。同时对遥感数据完成辐射定标、大气校正以及研究区裁剪。
[0092] 步骤2：将布设的7个先验水深数据点用公式(2)改正潮高，获取7个先验数据点卫星过境时刻水深，见表1。
[0093] 表1选择的7个先验水深数据点及其卫星过境时刻水深
[0094]
[0095] 步骤3：将7个先验水深数据点用于Strumpf的对数比值水深反演方法，反演结果为未潮高改正的水深数据反演结果，见图4；将7个先验数据点卫星过境时刻水深用于Strumpf的对数比值水深反演方法，反演研究区卫星过境时刻水深。
[0096] 步骤4：将反演的研究区卫星过境时刻水深用公式(1)改正潮高，获得以大地水准面(1985高程基准面)为起算面的水深数据，结果为潮高改正的水深数据反演结果，见图5。
[0097] 步骤5：对比7个控制点未潮高改正的水深数据反演水深、潮高改正的水深数据反演水深与真实水深见表2；对比4条控制线未潮高改正的水深、潮高改正水深与真实水深的剖面见图6；对比未潮高改正的水深数据反演水深、潮高改正的水深数据反演水深与真实水深直方图见图7。
[0098] 表2选择的7个控制点的未潮高改正的水深数据反演水深、潮高改正的水深数据反演水深与真实水深
[0099]
[0100] 从表2可知控制点1、2未潮高改正的水深数据反演水深更接近真实水深，而其他控制点中潮高改正的水深数据反演水深更接近真实水深是因为未潮高改正的水深数据反演水深集中在13‑15m内，参见图7。
[0101] 综上所述，遥感水深反演潮高改正方法在一定程度上能够提高水深反演的精度，使反演的水深更接近真实水深，也更有实际意义和价值。
[0102] 尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造概念，则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例以及落入本发明范围的所有变更和修改。
[0103] 显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包括这些改动和变型在内。
