CN113189559B - 一种星载成像高度计遥感数据海底地形反演方法 - Google Patents

Abstract

本发明提供一种星载成像高度计遥感数据海底地形反演方法，实现宽刈幅、高分辨率和高精度海底地形的反演。本发明的方法反演的海底地形主要由三部分组成，分别为长波海底地形、中波海底地形和短波海底地形，如下式所示。中波海底地形的确定是海底地形反演重点，其通常由相关分析确定的中波段重力异常与比例因子相乘得到。本发明为解决传统高度计海底地形反演中SAS法和GGM法的不足，考虑到成像高度计数据特点；通过确定精确的中波段海底地形，最终得到高精度的海底地形。

Description

一种星载成像高度计遥感数据海底地形反演方法

技术领域

本发明属于海洋遥感的技术领域，具体涉及一种星载成像高度计遥感数据海底地形反演方法。

背景技术

海底地形是海洋测绘的主要目标一，根据海底地形数据绘制的海底地形图有着广泛应用，它能够为研究海底地质构造提供海洋观测数据，为海洋环境研究提供测绘保障，为海洋生物学研究提供基本数据。海底地形图也可以应用于航海及航运、渔业捕捞、海上划界、海洋工程和海洋资源开发等领域。

目前，利用遥感手段观测海底地形的传感器主要有SAR和雷达高度计，其中，SAR主要通过后向散射强度变化对浅海海底地形进行反演，雷达高度计则利用海面高度数据进行海底地形反演。雷达高度计相比于SAR，其主要应用于较深海区的海底地形反演，应用更广泛。但雷达高度计由于数据采集方法的先天不足，仅能获得一维的沿轨数据，数据的空间分辨率和时间分辨率低，这就导致最终反演得到的海底地形精度较低。

而卫星测高技术可以全天时、全天候进行观测，得到大范围、高精度的海面高数据，通过数学模型反演重力异常，再利用重力异常和海深之间的相关性反演得到海底地形，因此，卫星测高技术已成为观测海底地形的一种重要手段。

近年来，利用雷达高度计获取的海面高度进行海底地形反演方法时，学者通常采用GGM法(重力地质法)和SAS法。GGM法是利用迭代法确定出地壳与海水密度差异常数，进行海底地形反演的方法，GGM法无需先验模型、不考虑地壳均衡状态，计算速度更为便捷，且其在船测海深控制点分布均匀的大面积海域精度较高。而SAS法是由Smith and Sandwell提出的抗线性回归方法确定比例因子的海底地形反演方法，SAS法通过在频率域内确定中波段重力异常与残余海深之间的比例因子，从而进行海底地形的反演，其不受计算面积大小的限制，船测点分布的影响也较小，但比例因子值易受某些因素的影响(如存偏离较大或较小的点、重力异常与地形之间存在的高次项、滤波波段截取了部分有用的信号等)，从而降低了海底地形的反演精度。

为获取大范围的海面高度，大幅宽成像高度计应运而生。目前，国际上典型的成像高度计有两个，一是我国研制的并已经成功发射的天宫二号微波成像高度计，另一是美国将要发射的SWOT成像高度计。天宫二号微波成像高度计能够对海洋进行干涉成像观测，得到海面高度等观测参数，观测刈幅可达35km；SWOT卫星观测刈幅可达120km。

成像高度计与雷达高度计一样，都可以通过测量海面高度间接反演海底地形，两种传感器的海底地形反演方法具有较强的相似性。但是相较于雷达高度计，成像高度计获得的海面高度数据具有宽刈幅、高分辨率特点，因此，借鉴现有海底地形反演方法，进行适当的改造或优化，建立适用于成像高度计数据的高精度海底地形反演方法，具有重要的应用意义。

发明内容

本发明的目的是提供一种星载成像高度计遥感数据海底地形反演方法，实现宽刈幅、高分辨率和高精度海底地形的反演；从而弥补现有技术的不足。

本发明首先提供利用成像高度计的海面高数据反演得到重力异常数据的方法，包括如下的步骤：

S1:星载成像高度计遥感数据处理：

对收集到的研究区内不同时间段的(沿轨)海面动力高度数据进行反距离加权平均处理，然后与平均海平面数据相加，从而确定出研究区域内的平均海面高度h；

所述的数据为成像高度计遥感观测的海面高度数据，海面地形数据、船测海深数据；

S2：利用处理后的海面高度数据求解残余垂线偏差，

S2-1:将处理后的海面高度数据减去海面地形得到大地水准面数据，再移去利用大地水准面模型计算得到的大地水准面数值，即可得到残余大地水准面数据

N＝h-δ-N_EGM2008

式中N为残余大地水准面高，h为海面高，δ为稳态海面地形，N_EGM2008为大地水准面模型；

S2-2：利用得到的残余大地水准面高数据得到研究区域的残余垂线偏差分量，

式中N、ξ、η均为矩阵数据，(i，j)为矩阵数据中每个元素的对应坐标，i是行，j是列；式中N为残余大地水准面高，

是纬度，λ是经度，单位为弧度；

S3：通过步骤S2计算得到的残余垂线偏差数据，进行残余重力异常的计算，

S3-1：利用逆Vening-Meinesz公式进行残余重力异常的计算，逆Vening-Meinesz公式具体表达如下：

式中，GM＝398600.4415km³/s²，为地球引力常数，R表示地球平均半径，p、q分别为计算点和流动点，Δ_qp表示流动点q到计算点p的方位角，ξ_q为流动点垂线偏差的子午分量，η_q为流动点垂线偏差的卯酉分量，H′为积分核函数的导数；

当计算点与流动点的球面距离为0，即S_qp＝0时，核函数H′奇异；当内圈带奇异时，重力异常的计算公式为：

s_o的计算式为：

ξ_y，η_x的表达式为：

S3-2：计算完残余重力异常与内圈带效应的重力异常后，加上重力异常模型的重力异常，即为最后所求的重力异常；

本发明再一个方面还提供一种星载成像高度计遥感数据海底地形反演方法，是使用上述方法获得的重力异常数据和船测海深数据进行海底地形反演，包括如下的步骤：

S1：基于重力异常与海深的相干性计算公式，进行相干性分析，确定出重力异常的反演波段；

式中，

为地形与重力异常的相干性，S_hh、S_gg为地形和重力异常的自功率谱密度函数，S_hg、S_gh为地形与重力异常的交叉功率密度函数；

S2：根据S1中确定的重力异常反演波段，对重力异常数据进行带通滤波并向下延拓到平均水面，获得反演波段内的重力异常；

S3：对下载的船测海深数据进行中误差限差处理，将剔除粗差后的船测数据2/3作为海深控制点，1/3作为检核点；对控制点数据进行格网化处理，并采用低通滤波得到长波海底地形模型；

S4：利用限定极值判定法求解比例因子，计算每个值对应的船测检核点观测海深与预测海深之间的标准差，再对所有计算的标准差进行线性拟合，对拟合曲线求导，标准差极小值点对应的系数，即为反演波段的重力异常与残余海深之间的比例因子；根据比例因子确定出中波段海底地形计算公式即为：

y＝S×G_中

式中，y为中波段海底地形，S为比例因子，G_中为反演波段重力异常值。

S5：比例因子与反演波段内的重力异常相乘即为中波海底地形模型；

S6：将船测水深数据格网化，采用低通滤波得到长波海底地形；长波和中波海底地形相加得到中长波海底地形，通过插值的方法，得到船测点处实测海深与中长波海深之差，将差值格网化后，生成短波海底地形；三者相加即为反演的海底地形。

本发明为解决传统高度计海底地形反演中SAS法和GGM法的不足，考虑到成像高度计数据特点，本发明的方法通过确定精确的中波段海底地形，最终得到高精度的海底地形。

附图说明

图1为本发明的海底地形反演技术路线图；

图2为本发明的成像高度计南海海面高度探测仿真图；

图3为本发明实施例的南海重力异常反演图；

图4为本发明实施例的重力异常与海底地形相关性估计图；

图5为本发明实施例的比例因子与标准差关系曲线图；

图6为本发明实施例的南海中波段海底地形反演图；

图7为本发明实施例利用OSAS法反演的南海海底地形结果图。

具体实施方式

本发明的方法反演的海底地形主要由三部分组成，分别为长波海底地形、中波海底地形和短波海底地形，如下式所示。中波海底地形的确定是海底地形反演重点，其通常由相关分析确定的中波段重力异常与比例因子相乘得到。因此，比例因子的确定成为反演的关键。

h_OSAS＝h_短+h_长+SG_中

式中，h_OSAS为研究区海底地形，h_短、h_长为短波与长波部分的海底地形，S为中波重力异常与中波海底地形之间的系数，G_中为中波段重力异常。

中波段重力异常的获取分为两个步骤，首先利用带通滤波对计算得到的重力异常进行处理，然后再对处理后的重力异常进行向下延拓至平均海深面，公式如下：

ΔG_中(k)＝ΔG(k)W(k)e^2πkd

式中，ΔG(k)和ΔG_中(k)分别为测高重力异常与经带通滤波和向下延拓处理后的重力异常。W(k)为带通滤波，e^2πkd为延拓因子，k为径向频率，d为研究区内的平均海深。

为获取长波部分的海底地形，首先需要对船测海深数据进行格网化，再进行低通滤波处理，最后获取的对应波段海深值，即为长波海底地形。

H_长(k)＝H_船(k)×W₁(k)

式中，H_长(k)、H_船(k)分别为长波部分海深值与格网化的船测海深在波数域内的表达，W₁(k)为低通滤波。

分别将经过以上两个步骤处理的ΔG_中(k)和H_长(k)转换至空间域(分别为G_中和h_长)，再通过插值处理将格网的船测海深与中波段重力异常转换至船测点上，利用船测点上的船测海深值计算短波海深值。

h_短(x)＝h_船(x)-h_长(x)-SG_中(x)

h_短(x)、h_船(x)、h_长(x)和SG_中(x)分别为船测点上的短波海深值、船测海深值、长波海深值和中波海深值。最后将船测点上的短波海深值进行格网化处理即得到短波海底地形h_短。

本发明采用的中波段重力异常值是利用遥感获得的海面高度数据计算得到的，经船测重力异常数据验证表明，精度高于模型重力异常数据，说明重力异常精度较高。保障了比例因子值确定的准确性。

综合考虑多种因素的影响，本发明提出求解比例因子S，确定中波段海底地形计算公式的方法，其保证了反演结果的精度。具体步骤如下：

(1)首先确定比例因子的经验值范围，在比例因子取值范围内，以0.5为步长，计算不同比例因子下船测检核点观测海深与预测海深之间的标准差，再对计算得到的标准差和对应比例因子进行拟合，对拟合曲线求导，标准差极小值点对应的系数作为该区域比例因子。

(2)比例因子与中波段重力异常相乘即为中波段海底地形值，具体计算公式如下式所示:

y＝9.5×G_中

式中，y为中波段海底地形，G_中为中波段重力异常。

(3)将船测水深数据格网化，采用低通滤波得到长波海底地形；长波和中波海底地形相加得到中长波海底地形，通过插值的方法，得到船测点处实测海深与中长波海深之差，将差值格网化后，生成短波海底地形。

(4)中长波海底地形与短波海底地形叠加得到最终的OSAS海底地形。

下面结合实施例和附图，对本发明的具体实施方式作进一步详细描述。

实施例1：建立海底地形反演步骤

SAS法反演海底地形不受区域大小的限制，船测点分布的影响也较小，反演浅海区域海底地形的精度较高。相比于GGM法，SAS法更具有普适性，但现有比例因子选取方式得到的数值由于受到较多因素影响，会降低海底地形精度。

基于以上分析，如图1所示，本发明提出了一种成像高度计遥感数据海底地形反演方法。因此，比例因子的选取方法有待改进，以便确定出精确的中波段海底地形反演公式，其包括以下步骤：

S1:数据处理：

S1-1：确定研究区域，收集相应区域的成像高度计遥感观测的海面高度数据，海面地形数据、船测海深数据与重力异常数据。

S1-2：对收集到的研究区内不同时间段的(沿轨)海面动力高度数据进行反距离加权平均处理，然后与平均海平面数据相加，从而确定出研究区域内的平均海面高度h。

S2：利用处理后的海面高度数据求解残余垂线偏差。因成像高度计得到的海面高数据为二维的网格化数据，因此在计算垂线偏差时可直接进行计算，不需进行雷达高度计计算垂线偏差时采用的格网化和间接平差等操作，一定程度上提高了计算的精度与速率。

垂线偏差的计算步骤具体如下：

S2-1:将处理后的海面高度数据减去海面地形得到大地水准面数据，再移去利用大地水准面模型计算得到的大地水准面数值，即可得到残余大地水准面数据。

N＝h-δ-N_EGM2008

式中，N为残余大地水准面高，h为海面高，δ为稳态海面地形，N_EGM2008为大地水准面模型。

S2-2：利用得到的残余大地水准面高数据，结合下式即可得到研究区域的残余垂线偏差分量。

式中N、ξ、η均为矩阵数据，(i，j)为矩阵数据中每个元素的对应坐标，i是行，j是列。式中N为残余大地水准面高，

是纬度，λ是经度，单位为弧度。

S3：通过步骤S2计算得到的残余垂线偏差数据，进行残余重力异常的计算。计算重力异常的过程主要利用了经典的“移去-恢复”法，在计算残余垂线偏差时，已移除模型的大地水准面数据，此为移去过程。重力异常模型与残余重力异常计算结果相加得到最终的重力异常即为恢复过程。具体步骤如下：

S3-1：首先，利用逆Vening-Meinesz公式进行残余重力异常的计算，逆Vening-Meinesz公式具体表达如下：

式中，GM＝398600.4415km³/s²，为地球引力常数，R表示地球平均半径，p、q分别为计算点和流动点，α_qp表示流动点q到计算点p的方位角，ξ_q为流动点垂线偏差的子午分量，η_q为流动点垂线偏差的卯酉分量，H′为积分核函数的导数。

当计算点与流动点的球面距离为0，即S_qp＝0时，核函数H′奇异。受卫星测高空间分辨率的限制，核函数发生奇异的面积从几平方千米到几十平方千米，称为奇异积分区，反演时必须考虑内区的奇异性效应。当内圈带奇异时，重力异常的计算公式为：

s_o的计算式为：

ξ_y，η_x的表达式为：

S3-2：计算完残余重力异常与内圈带效应的重力异常后，加上重力异常模型的重力异常，即为最后所求的重力异常。

S3-3：利用研究区内实测的重力异常数据以及EGM2008A模型对计算得到的重力异常数据进行精度检验，以保证海底地形反演的精度。

利用成像高度计的海面高数据反演得到的重力异常数据和船测海深数据即可进行海底地形反演，包括如下的步骤：

S4：基于重力异常与海深的相干性计算公式，进行相干性分析，确定出重力异常的反演波段。

式中，

为地形与重力异常的相干性，S_hh、S_gg为地形和重力异常的自功率谱密度函数，S_hg、S_gh为地形与重力异常的交叉功率密度函数。

S5：设计滤波器。考虑到Gibbs现象的影响，把带通滤波器取为高通滤波器与低通滤波器两个滤波器的乘积可以有效抑制长波和短波的奇异性。根据S4中确定的重力异常反演波段，对重力异常数据进行带通滤波并向下延拓到平均水面，获得反演波段内的重力异常。

S6：对下载的船测海深数据进行中误差限差处理，将剔除粗差后的船测数据2/3作为海深控制点，1/3作为检核点。对控制点数据进行格网化处理，并采用低通滤波得到长波海底地形模型。

S7：确定比例因子。利用限定极值判定法求解比例因子，具体步骤如下：

S7-1：在比例因子取值范围内，通过限定极值判定法，计算每个值对应的船测检核点观测海深与预测海深之间的标准差，再对所有计算的标准差进行线性拟合，对拟合曲线求导，标准差极小值点对应的系数，即为反演波段的重力异常与残余海深之间的比例因子。根据比例因子确定出中波段海底地形计算公式即为：

y＝S×G_中

式中，y为中波段海底地形，S为比例因子，G_中为反演波段重力异常值。

S8：比例因子与反演波段内的重力异常相乘即为中波海底地形模型。将船测水深数据格网化，采用低通滤波得到长波海底地形；长波和中波海底地形相加得到中长波海底地形，通过插值的方法，得到船测点处实测海深与中长波海深之差，将差值格网化后，生成短波海底地形。三者相加即为最后求得的海底地形。

S9：对S8得到的海底地形模型进行精度验证。将海底地形模型数据内插到船测检核点处与检核点处的海深值作比较，统计相对偏差，从而进行精度评价。

实施例2：海底地形反演方法的应用

下面以中国南海(112°E～119°E，12°N～19°N)范围内，2020年1月1日到7日，连续7天分辨率为1/24°的成像高度计海面高度探测仿真数据为例，说明本发明方法海底地形反演的具体过程。

1)海面高数据预处理

海面高数据预处理主要是通过对获取的研究区域内多幅7天为一组的海面动力高度数据进行反距离加权平均，再与平均海平面数据相加，从而确定出研究区域内最终的海面高度，如图2所示。此过程能够削弱海面时变信息对海面高的影响。

2)计算垂线偏差

利用处理后的海面高度数据减去DOT2008A海面地形，得到大地水准面高数据，再移去利用EGM2008A大地水准面模型，得到垂线偏差计算所需的残余大地水准面数据。根据垂线偏差计算公式，计算出格网化的残余垂线偏差数据分辨率为1/24°。

3)计算重力异常

重力异常的计算采用了“移去-恢复”法，移去的过程为上述移去大地水准面模型的过程。基于残余垂线偏差数据，利用逆Vening-Meinesz公式计算残余重力异常与内圈带内的重力异常，最后，加上EGM2008A重力异常模型数据，得到研究区内的重力异常数据，此过程即为恢复的过程。重力异常的计算结果如图3所示。利用实测船测重力异常数据(来自NGDC的船测数据)并结合EGM2008A模型对计算得到的重力异常数据进行精度检验，如表1统计结果显示，本方法得到的重力异常反演结果精度较高，也进一步说明了成像高度计的二维平面数据能够提高重力异常的计算精度。

表1：试验重力异常与船测数据插值统计表

4)确定重力异常反演波段

对计算得到的重力异常数据与研究区内的船测海深数据进行相干性分析。确定出研究区域的重力异常与海深的相干性波段。如图4所示，当波长在10-120km波段范围内时，相干性较强。

5)船测数据处理

船测海深为测线数据，来自NGDC(美国地球物理中心)。对船测海深数据进行中误差限差处理，剔除粗差数据，将处理后的船测数据的2/3作为海深控制点，1/3作为检核点。对控制点数据进行插值处理，得到格网化海深数据。

6)滤波处理

根据确定的重力异常反演波段，对重力异常数据进行带通滤波并向下延拓至平均水面，获得反演波段内的重力异常；对格网化船测海深数据进行低通滤波处理得到长波海底地形。

7)确定比例因子

首先确定该系数的经验值范围，在取值范围内，通过限定极值判定法，计算每个值对应的船测检核点观测海深与预测海深之间的标准差，再对所有计算的标准差进行线性拟合，对拟合曲线求导，标准差极小值点对应的系数作为二者的比例因子，即为反演波段的重力异常与残余海深之间的比例因子。如图5所示为比例因子与标准差的关系曲线图。本发明确定该值大小为9.5m/mGal。因此，该区域中波段海底地形计算公式为y＝9.5×G_中。

8)海底地形反演

如图6所示，比例因子与反演波段内的重力异常相乘即为中波海底地形模型。将长波和中波海底地形相加得到中长波海底地形，通过插值方法，得到船测点上实测海深与中长波海深之差，将差值格网化后，生成短波海底地形模型。中长波地形与短波地形叠加得到海底地形模型，如图7所示为反演结果图。

9)精度检验

将海底地形模型数据内插到船测检核点处与检核点处的海深值作比较，统计相对偏差，确定模型的精度。结果如表2所示，相对偏差为2.71％，在5％以内精度较高。

表2：模型检核统计结果表

本发明方法确定出中国南海(12°N～19°N，112°E～119°E)在海底地形反演过程的中波段重力异常与中波段海深之间的比例因子为9.5m/mGal，为该区域海底地形的精确反演提供有力支撑。

Claims (2)

1.一种星载成像高度计遥感数据海底地形反演方法，其特征在于，所述的方法是利用成像高度计的海面高数据反演得到重力异常数据的方法获得的重力异常数据和船测海深数据进行海底地形反演；

所述得到重力异常数据的方法包括以下步骤：

星载成像高度计遥感数据处理：

对收集到的研究区内不同时间段的海面动力高度数据进行反距离加权平均处理，然后与平均海平面数据相加，从而确定出研究区域内的平均海面高度h；

利用处理后的海面高度数据求解残余垂线偏差：

将处理后的海面高度数据减去海面地形得到大地水准面数据，再移去利用大地水准面模型计算得到的大地水准面数值，得到残余大地水准面数据：

N＝h-δ-N_EGM2008

式中，N为残余大地水准面高，h为海面高，δ为稳态海面地形，N_EGM2008为大地水准面模型；

利用得到的残余大地水准面高数据得到研究区域的残余垂线偏差分量，

式中N、ξ、η均为矩阵数据，(i，j)为矩阵数据中每个元素的对应坐标，i是行，j是列；式中N为残余大地水准面高，

是纬度，λ是经度，单位为弧度；

通过所述残余垂线偏差分量，进行残余重力异常的计算：

利用逆Vening-Meinesz公式进行残余重力异常的计算，逆Vening-Meinesz公式具体表达如下：

式中，GM＝398600.4415km³/s²，为地球引力常数，R表示地球平均半径，p、q分别为计算点和流动点，α_qp表示流动点q到计算点p的方位角，ξ_q为流动点垂线偏差的子午分量，η_q为流动点垂线偏差的卯酉分量，H′为积分核函数的导数，其中，当计算点与流动点的球面距离为0，即S_qp＝0时，核函数H′奇异；当内圈带奇异时，重力异常的计算公式为：

s_o的计算式为：

ξ_y，η_x的表达式为：

计算完残余重力异常与内圈带效应的重力异常后，加上重力异常模型的重力异常，获取所述重力异常数据；

所述反演方法包括如下的步骤：

S1：基于重力异常与海深的相干性计算公式，进行相干性分析，确定出重力异常的反演波段：

式中

为地形与重力异常的相干性，S_hh、S_gg为地形和重力异常的自功率谱密度函数，S_hg、S_gh为地形与重力异常的交叉功率密度函数；

S2：根据S1中确定的重力异常反演波段，对重力异常数据进行带通滤波并向下延拓到平均水面，获得反演波段内的重力异常；

S3：对下载的船测海深数据进行中误差限差处理，将剔除粗差后的船测数据2/3作为海深控制点，1/3作为检核点；对控制点数据进行格网化处理，并采用低通滤波得到长波海底地形模型；

S4：利用限定极值判定法求解比例因子，计算每个值对应的船测检核点观测海深与预测海深之间的标准差，再对所有计算的标准差进行线性拟合，对拟合曲线求导，标准差极小值点对应的系数即为反演波段的重力异常与残余海深之间的比例因子；根据比例因子确定出中波段海底地形计算公式即为：

y＝S×G_中

式中，y为中波段海底地形，S为比例因子，G_中为反演波段重力异常值；

S5：比例因子与反演波段内的重力异常相乘即为中波海底地形模型；

S6：将船测水深数据格网化，采用低通滤波得到长波海底地形；长波和中波海底地形相加得到中长波海底地形，通过插值的方法，得到船测点处实测海深与中长波海深之差，将差值格网化后，生成短波海底地形；三者相加即为反演的海底地形。

2.如权利要求1所述的方法，其特征在于，所述星载成像高度计遥感数据包括海面高度数据，海面地形数据、船测海深数据。

Images