应用抗差估计方法构建日本海海底地形模型

范 雕，李姗姗，孟书宇，张金辉，单建晨，王傲明

（1.信息工程大学，郑州 450001；2.32022 部队，武汉 430074）

海洋测绘作为认识、开发、利用海洋的基础手段和系统工程，测绘成果可为维护国家海洋权益、发展海洋经济、保护海洋环境等国家重大发展战略顺利实施提供基础数据支撑。海底地形测量作为海洋测绘重要组成部分，数据及产品广泛应用于国民经济建设、国防建设以及科学研究等方方面面。目前海底地形数据获取方式仍主要通过船舶为载体搭载单波束或者多波束等声学仪器进行测量，面对占地球表面积大约71%的浩瀚海洋，声学测深数据覆盖范围十分有限，如最新版本的全球海底地形模型SRTM15+V2.0，其中收集的全球声学海深测量数据在空间格网大小为15"环境下采样率也仅为10.84%[1]。有学者指出，若望采用船基声学海底地形测量获取覆盖全球的深海水深数据，需要持续投入巨额资金，耗费大约200 船年（ship-year）时间方可完成[2]。因此面对海洋观测系统化、全球化发展趋势，探索低成本、高效率的海底地形数据获取手段显得尤为重要。

地壳均衡补偿理论表明海面重力数据与海底地形呈现较强相关关系，这种相关关系使得利用海面重力信息恢复海底地形成为可能[3]。特别是随着测高卫星、重力卫星等空间大地测量技术不断发展，依托星基平台高效快速地恢复了高精度全球海洋重力场信息[4,5]，从而国内外学者对采用海面重力数据反演海底地形开展了大量有益研究[6-9]，目前常用的依据海面重力数据反演海底地形算法主要有重力地质法（GGM 方法）[10,11]，回归分析技术[12-14]，导纳函数方法[15-17]，最小二乘配置技术[18,19]，模拟退火方法[20]等。其中利用海底地形和海面重力信息线性关系反演海底地形（称为回归分析技术）的关键核心是确定地形-重力比例因子。关于比例因子确定问题，国内外学者做了大量的研究探索工作，文献[12][21][22]采用最小二乘的方法获取比例因子，文献[23][24]直接将船测点残余海底地形和残余重力数据相除得到比例因子，然而若输入数据存在粗差，采用此两类比例因子获取方法会严重影响海底地形反演质量。文献[25]文献为了解决该问题，采用反演波段海底地形和反演波段重力数据的中位数之比获得比例因子，并强制地形-重力过原点，该方法具有一定的抗差性，但是忽略了地形-重力关系中常数项对海底地形模型构建结果的影响。

基于以上分析，兼顾输入数据中粗差对海底地形反演结果干扰和地形-重力关系中常数项对海底地形模型构建影响，本文提出采用抗差估计理论确定海底地形-重力比例因子，进而构建合理海底地形模型的新思路。然后以日本海部分海域作为算法验证区，分别利用试验海区卫星测高重力异常、卫星测高重力异常垂直梯度数据，结合稀疏船测水深数据，建立了相应的海底地形模型，并将本文海底地形模型反演结果与采用最小二乘方法、中位数方法构建的海底地形模型、DTU10 模型和ETOPO1 模型进行了模型横向比对和精度评估。

1 原理和方法

1.1 基本理论

Parker 公式推导了以海底地形为输入，海面重力异常为输出的频率域表达式：

式中：G为地球引力常数；Δρ为海水和地壳的密度差异；d为平均海深；F（Δg(x,y)）和F（h(x,y)）分别表示重力异常和海底地形的二维傅里叶变换；k为波数kx和ky分别表示x和y方向的波数（kx=2π/λx，ky=2π/λy）；λx、λy分别表示x和y方向的波长。当仅考虑式(1)中线性项时，海底地形和海面重力异常的频域线性关系可近似表示为

式中，Z(kx,ky)表征将海底地形转换为重力异常的能力，通常称为重力导纳函数。重力导纳函数Z(kx,ky)具有各向同性、空间不变性特点。研究表明，大尺度的海底地形特征可以用地壳均衡补偿理论来解释，那么当顾及地壳均衡补偿影响时，式(2)可表示为

式中，Zf(kx,ky)称为挠曲导纳函数；φ(kx,ky)为均衡响应函数：

式中，Tc为地壳厚度；Φe(kx,ky)为挠曲均衡响应函数（Flexural Response Function）。分析(3)式可知，均衡响应函数将由地形产生的重力异常修正为地形和均衡补偿物质对重力异常的贡献之和。挠曲均衡响应函数定义为：

其中，g为研究区域重力值；D为挠曲刚度，计算式如下：

式中，E为杨氏模量（Young’s Modulus）；Te为板块的有效弹性厚度；υ为泊松比。据式(3)可得利用重力异常反演海底地形算法模型为

式中，转换函数Qf(kx,ky)为导纳函数Zf(kx,ky)的倒数。为了抑制均衡补偿对海面重力异常的影响和克服重力异常向下延拓过程（转换函数Qf(kx,ky)包含向下延拓因子ekd）噪声放大的缺点以获得稳定的海底地形反演结果，引入带通滤波器W(kx,ky)限制转换函数Qf(kx,ky)，则

式中，h0为对应波段的海底地形。那么在未考虑地壳均衡补偿环境下，同时根据频率域求导法则，可得重力异常/重力异常垂直梯度（Δgz=-(∂Δg/ ∂z)）为系统输入，海底地形为系统输出的频率域表达式

由式(9)可知，重力异常/重力异常垂直梯度经滤波等处理后与海底地形呈线性比例关系，比例因子为S(x,y)=(1/2π GΔρ)。从而可以认为海底地形反演工作通常在中长波波段范围开展，而短波和长波部分海底地形信息可通过船测海深等方式进行补充，进而获得全波段的海底地形

式中，hl(x,y)为非反演波段海底地形；Δg0(x,y)为经滤波且向下延拓后的重力异常；Δgz0(x,y)为经滤波且向下延拓等处理后的重力异常垂直梯度相关量，其他参数意义同上。经分析可知，依据式(10)反演海底地形，关键在于确定比例因子。然而实际上由于海底地形地貌构造复杂和沉积层密度异常等因素影响，比例因子并非常数而是一个变化量[25]，不能根据理论密度参数直接计算，而是通过船测数据约束，分别计算不同区域的线性参数[22]。

图1 海底地形反演流程图Fig.1 Flow chart of seafloor topography inversion

基于以上分析，绘制海底地形的反演流程如图1所示。需要说明的是，传统方法首先获取控制点处离散比例因子，然后将比例因子格网化处理，进而构建海底地形模型[25]。文献[12]指出，依据移动窗口法获取比例因子构建海底地形质量高于传统方法。据此，本文采用移动窗口法求取比例因子参数。

1.2 比例因子获取方法

1.2.1 中位数方法（Median 方法）

假设反演海区有限频段海底地形和重力相关数据分别为b0(x,y) 和g0(x,y)，则二者比例因子可表示为：

式中，med(b0(x,y))和med(g0(x,y))分别表示地形和重力相关数据的中位数。

1.2.2 最小二乘方法（Ls 方法）

为获得海底地形-重力的比例因子，以重力数据为自变量，海底地形为因变量进行线性拟合获取比例因子和常数项参数结果。依据最小二乘方法解算参数结果为

1.2.3 抗差估计方法（Ro 方法）

抗差估计（RobustEstimation）也称为稳健估计或者鲁棒估计，源于统计学中的抗差性（Robustness）概念。当观测值中存在粗差的情况下，采用抗差估计可以得出正常模式下的最佳估值。抗差估计的原则是充分利用有效信息，限制利用有用信息，排除有害信息，在所假定的模型下，可以获得可靠、有效、具有实际意义的参数估值。

依据抗差估计理论，参数解算迭代策略为

式中，k为迭代次数；为拟合参数结果；A为重力数据构成的系数矩阵；L为海底地形输入向量；为等价权矩阵。中的元素为

式中，称为等价权；ω(vi)称为权因子；ψ(vi)是非线性函数；vik表示第k次迭代的残差分量。实际计算等价权的操作中，首先需将残差vi标准化

式中，miv是vi的中误差；Ai表示系数阵A的第i行向量；σ0是单位权中误差，采用标准化残差绝对值的中位数计算

当参数x 的迭代解满足式(17)，迭代停止，获得最终解。

式中，ε为设置的迭代停止条件。本文选择的等价权函数为Huber 函数

由式(18)可知，当≤2σ0，Huber 函数等价于最小二乘方法，当＞2σ0，Huber 函数利用L1范数对观测值进行降权，这样可以有效抑制粗差的不利影响，同时保持良好的正态有效性。

值得说明的是，Ro 方法求解比例因子本质核心仍可认为是最小二乘估计。不同的是，Ls 方法获取比例因子仅仅是有船测海深数据约束的最小二乘方法；Median 方法使用排序信息获取比例因子；Ro 方法获取比例因子是具有外部Huber 约束地形信息的正则化最小二乘法。从而算法本质上是希望对最小二乘施加约束，进而获得更有效、合理的比例因子。

2 数值分析试验

2.1 数据准备和前期处理

选取日本海4 °×4 °（131 °E～135 °E,36 °N～40 °N）海域作为试验海区开展数值分析试验。该海域西北方向地形平坦，东北和东南方向存在海山等地形复杂地貌，从而认为该海域具备作为算法可靠性验证的典型区域条件。试验过程以船测点作为约束条件补充海底地形短波部分和长波部分地形信息。数据来源如下：①船测海深数据：来源于 NGDC（The National Geophysical Data Center）发布的研究海区实测水深数据。研究海域原始船测水深数据共39018 个，通过3σ准则对实测水深数据进行粗差剔除，最终得到37750个实测水深数据。选取其中约五分之四的数据（30877个）作为控制点，以期对构建的海底地形模型进行控制；剩余6873 个实测水深数据作为外部检核点供海底地形模型精度评估使用。控制点和检核点的分布如图2(a)所示，其中白色六边形为控制点，黑色三角形为检核点，图中背景为ETOPO1 海深模型。②卫星测高重力数据包括重力异常和重力异常垂直梯度数据均来自于SIO，UCSD（Scripps Institution of Oceanography,University of California,San Diego），V24.1 版本，分辨率为1’，研究海域重力异常和重力异常垂直梯度数据如图2(b)(c)所示。

图2 研究海域数据Fig.2 Study sea area data

2.2 海底地形模型构建

根据利用卫星测高重力数据反演海底地形理论可知，海底地形模型构建过程中，由于短波部分重力数据所含海底地形信噪比较低以及重力数据向下延拓存在噪声放大现象，从而采用低通滤波器对重力数据进行滤波处理以消除短波部分反演结果不稳定影响。为尽可能使用重力信息，同时考虑低通滤波器和转换函数特点以及反演分辨率等因素，结合文献[21]和文献[25]，最终设定短波截止波长为15 km。另外，为有效抑制地壳均衡补偿对海面重力信息的影响，充分利用卫星测高重力信息中的地形分量，依据文献[26]认为研究海域有效弹性厚度约为5 km，然后根据有效弹性厚度与波长的关系，长波截止波长最终设定为160 km。基于以上分析，日本海试验海域利用卫星测高重力数据反演海底地形的反演波段范围为15 km～160 km。

分别以重力异常和重力异常垂直梯度作为系统输入，采用Median 方法、Ls 方法和Ro 方法计算试验海区线性分析一次项和常数项参数结果分别如图3～图4所示。

图3 试验海区topo-vgg 线性分析结果（无粗差）Fig.3 Topo-vgg linear analysis results in experimental seaarea (no gross error)

图4 试验海区topo-dg 线性分析结果（无粗差）Fig.4 Topo-dg linear analysis results in experimental sea area(no gross error)

最终，以重力异常和重力异常垂直梯度为输入源，依据各方法获得的回归参数结果构建相应的海底地形模型如图5和图6所示。

观察图5和图6海深模型构建结果，六种海深模型均展现了研究海域地形地貌特点：海域西北部以海底平原为主要地形呈现；海域东北和东南部地貌形态丰富，其中海山大面积分布于该区域，海域东部边缘存在海底平原地貌。仔细对比图5和图6海深模型可以发现，使用Median 方法，以重力异常为输入数据构建的海深模型细节呈现明显丰富，如图5(a)西北部可以清晰看见3 处地形隆起，而图5(b)(c)、图6(a)～(c)目视仅1 处地形隆起，另外两处则较为模糊；又如在（39 ° N～40 ° N，133 ° E～134 ° E）区域，图5(a)中地形呈现细节明显好于其余模型。

图5 重力异常构建海深模型Fig.5 Bathymetry model derived by gravity anomaly

图6 重力异常垂直梯度构建海深模型Fig.6 Bathymetry model derived by vertical gravity gradient anomaly

分析沉积层对比例因子和相关系数结果影响，依据CRUST 1.0 绘制试验海域沉积层厚度和密度分布图如图7所示。比对图3、图4和图7可大致发现，海底地形-重力比为适应沉积层厚度和密度变化而变化：在沉积层厚度较薄区域，沉积层密度较小，海水和地壳密度差异较大，从而海底地形-重力呈现出较强相关性，如研究海域东北部区域；沉积层厚度较厚区域，沉积层密度较大，海底地形起伏较小，覆盖在地形表面的沉积层也会对海面重力异常产生影响，该区域地形-重力相关性较弱，如西南部和西北部研究海域。另外，对于海底地形和重力数据比例因子较大的原因，笔者认为可能是，因地壳均衡补偿作用对海面重力影响，导致该海区重力值（重力异常/重力异常垂直梯度）较小，从而海底地形和重力异常比例因子较大。虽然实验开始阶段输入数据经过滤波处理，然而因地形挠曲波长的选择为试验海区大范围内的平均值，当部分海域地形波长可能小于挠曲波长时，该区域可能发生局部均衡补偿现象，从而导致海面重力受到均衡补偿影响。此时可通过缩小滤波窗口，减小挠曲波长的方法抑制均衡补偿影响，然而缩小滤波窗口将引发其他问题。

图7 试验海域沉积层Fig.7 Sediment in study area

2.3 精度评价分析

以检核点作为外部检核条件，引入国际广泛使用的ETOPO1 海深模型和DTU10 海深模型进行比较，评估各海深模型精度和效能。采用双线性插值方法将模型海深值内插到外部检核点处，并与检核点处实测水深值作差以评估海深模型精度。

统计结果表明，以重力异常为数据源，采用Median 方法建立的BGMed 海深模型检核差值最大值和最小值绝对值均超过 1000 m，检核均方差为243.37 m，检核精度远低于基于Ls 方法和Ro 方法构建的BGLs 模型和BGRo 模型。其中采用Ro 方法反演得到的BGRo 海深模型的检核均方差为90.03 m，Ls 方法反演海深模型的检核均方差为94.04 m，可以认为利用Ro 方法与Ls 方法在研究海域构建的海底地形模型整体效果相当，因此，从整体上看，Ro 方法具有良好的适用性。ETOPO1 海深模型和DTU10 海深模型的检核均方差为99.94 m 和125.42 m，检核精度高于BGMed 模型，低于BGLs 模型和BGRo 模型。采用重力异常垂直梯度为数据源建立的BVGMed 模型、BVGLs 模型和BVGRo 模型中，BVGRo 海深模型检核精度最高，为89.51 m，BVGMed 海深模型检核精度最低，为264.77 m。BVGRo 模型相较于ETOPO1海深模型和DTU10 海深模型精度分别提高了10.44%和28.63%左右。对比不同卫星测高重力数据构建的海深模型，BVGLs 海深模型精度优于BGLs 海深模型，BGRo 海深模型检核精度低于BVGRo 海深模型。由于重力异常和重力异常垂直梯度随距离的衰减特性不同，为探究二者面对浅海域和深海域表现状态，绘制BGRo 检核差值绝对值和BVGRo 检核差值绝对值之差与海深关系如图8所示。考察图8可以发现，本文试验区域水深小于1500 m 左右深度，BVGRo 检核差值绝对值比较明显小于BGRo 检核差值绝对值，即该水深范围内重力异常垂直梯度反演效果较好；水深大于1500 m 左右范围，BGRo 检核差值绝对值比较明显小于BVGRo 检核差值绝对值，即该水深范围内重力异常反演效果较佳。笔者认为从重力异常和重力异常垂直梯度数据的特性角度考虑：重力异常反映的是相对低频的信息，而重力异常垂直梯度反映的是相对高频的信息，因此海域较深，重力异常数据的反演结果相对较好，而浅部则是重力异常垂直梯度较好。

就外部检核结果统计指标整体而言，Ro 方法和Ls 方法在试验海域构建模型精度相当，优于Median方法；而Median 方法对于试验海域地形呈现丰富度又好于Ro 方法和Ls 方法。可见，三种方法各自具有不同的海底地形构建优势。

图8 BGRo 检核差绝对值和BVGRo 检核差绝对值之差与海深关系Fig.8 Relationship between the difference between checking difference absolute value of BGRo and BVGRo and sea depth

分别统计以重力异常和重力异常垂直梯度为数据源，采用Median 方法、Ls 方法和Ro 方法构建的海深模型不同检核差值范围内检核点数量，结果如图9所示。图9显示Median 方法构建的海底地形模型检核差值在不同差值区间数量稳定；Ro 方法构建的海深模型检核差值小于50 m 左右范围内检核点数量明显多于Ls 方法和Median 方法，检核差值大于50 m 范围内检核点数量明显少于Ls 方法和Median 方法，进一步说明了Ro 方法良好的估计性质和适用性。进一步验证Ro 方法的优越性，首先以试验海域反演波段海底地形和重力异常/重力异常垂直梯度为研究对象，分别选取其中3025 个数据点进行分析研究，包括反演波段海底地形、反演波段重力异常和重力异常垂直梯度数据。试验过程中认为原始反演波段输入数据中不含粗差，原始反演波段海底地形、反演波段重力异常和重力异常垂直梯度输入数据（分别简称为topo、dg 和vgg）统计结果如表1。实验过程中认为重力数据质量良好无粗差存在[25]，分别在与dg 和vgg相拟合的topo 中随机加入粗差（在原反演波段地形数据基础上增加或减小5～10 倍中误差作为粗差），最终与dg 相拟合的topo 数据中加入粗差点数为24 个，约占总点数的7.93‰，得到含粗差的地形数据topoErrDg；与vgg 相拟合的topo 数据加入粗差的点数为12 个，约占总点数的 3.97‰，获得含粗差的地形数据topoErrVgg。添加粗差的地形数据统计结果见表1。相较于添加粗差前的topo 数据，添加粗差后的地形数据最值变化明显，数据整体波动剧烈。

图9 模型比较Fig.9 Bathymetry model comparison

表1 输入数据统计结果Tab.1 Statistical results of input data

分别采用Median 方法、Ls 方法和Ro 方法确定地形-重力比例因子，其中反演波段海底地形和重力异常/重力异常垂直梯度线性回归分别简称为topo-dg 和topo-vgg。未加粗差的情况下，采用三种线性分析技术获得海底地形和重力异常/重力异常垂直梯度的拟合结果如图10(a)(c)所示，图10(a)(c)显示海底地形和重力表现良好的线性性质。在地形数据随机加入粗差后与相应的重力数据（dg 和vgg）采用相同的方法获得的线性分析结果如图10(b)(d)所示，图10(b)(d)中蓝色圆点表示加入粗差后的地形数据结果。图10(b)(d)显示，加入粗差对抗差估计结果影响较小。

图10 海底地形和重力数据Fig.10 Seafloor topography and gravity data

三种方法拟合参数统计结果表明，反演波段海底地形数据加入粗差后，海底地形和重力数据相关性明显降低，如未加粗差的海底地形和重力异常垂直梯度的相关系数为0.76，加入粗差的海底地形和重力异常垂直梯度相关系数为0.65。因此海底地形和重力回归分析过程中应充分利用有效信息，限制利用有用信息，排除有害信息，进而尽量避免粗差对结果的影响。经计算，不含粗差情况下，Median 法、Ls 法和Ro 法获得的比例因子分别为28.54、14.95 和16.29；加入粗差后三种方法获得的比例因子分别为28.59、13.78 和16.32，与不含粗差环境相比，分别在原值基础上变动1.75‰、7.83%和1.84‰，由此可见粗差对最小二乘影响较大，对抗差估计结果影响较小。同理可得加入粗差后，三种线性分析方法解算的海底地形和重力异常垂直梯度的比例因子较不含粗差情况分别变动0.00%、2.42%和0.00%，体现了抗差估计在粗差环境中优良的线性分析性质。

基于以上分析，在每个移动窗口内的初始反演波段海底地形输入数据中随机加入粗差（粗差点占总数的3%～6%），然后分别依据Median 方法、Ls 方法和Ro 方法构建输入海底地形数据含粗差的海深模型，构建的海深模型检核统计结果如表2。

表2 海深模型（含粗差）差值检核统计结果（单位：米）Tab.2 Statistical checking results of difference bathymetry model (including gross errors) (Unit:m)

为描述方便，以重力异常为输入源，采用Median方法、Ls 方法和Ro 方法建立的海深模型分别称为BGMedErr 模型、BGLsErr 模型和BGRoErr 模型；以重力异常垂直梯度为输入数据，采用Median 方法、Ls 方法和 Ro 方法构建的海深模型分别称为BVGMedErr 模型、BVGLsErr 模型和BVGRoErr 模型。表2检核统计结果显示，在输入海底地形数据中添加粗差，将不同程度地影响海深模型检核精度。如以重力异常为系统输入，BGMedErr 模型检核差值均方差为244.21 m，BGMed 模型与检核点差值均方差为243.37 m，反映了Median 方法具有一定的抵抗粗差能力。但是采用Median 方法获得的海深模型，由于反演过程强制拟合过原点，忽略了常数项对结果的影响，最终的海深模型反演精度低于 BGLsErr 模型和BGRoErr 模型。BGLsErr 模型相较于BGLs 模型检核均方差变化了16.88 m，而BGRoErr 模型相较于BGRo模型检核均方差仅变化了0.02 m，证明了采用Ro 方法反演海深模型能够充分利用有效信息，排除有害信息，获得稳定、有效的解算结果。以重力异常垂直梯度构建的海深模型检核统计结果也表明，Median 方法获得海深模型精度最低，Ro 方法获取的海深模型检核精度最高。同时BVGLsErr 模型相较于BVGLs 模型检核均方差变化了11.73 m，而BVGRoErr 模型相较于BVGRo 模型检核均方差仅变化了0.25 m。

各个海深模型加入粗差前后不同检核差值区间检核点变化情况如图11所示。可发现加入粗差后对Ls 方法构建的海深模型影响相对较大：不同检核差值区间检核点数量变化明显，检核差值小的检核点数量减小，检核差值大的检核点数量增多；而加入粗差后Ro 方法构建的海深模型不同检核差值区间检核点数量并未发生较大变化，进一步证明了Ro 方法对待含粗差输入数据具备较高稳定性。

图11 加入粗差后模型比较Fig.11 Bathymetry model comparison after adding gross errors

3 结 论

针对当前利用回归分析技术恢复海底地形算法中有限波段重力数据与海底地形之间比例因子获取忽略常数项或者抵抗粗差能力弱缺陷，提出了可采用抗差估计理论，实现兼顾输入数据粗差影响和常数项的比例因子参数求解新思路。选取日本海部分海域作为算法验证区，以卫星测高重力异常和重力异常垂直梯度为重力输入，分别使用最小二乘法、中位数法和抗差估计法恢复了试验海域海底地形；同时讨论了初始海底地形存在粗差环境下，三种比例因子获取方法抵抗粗差影响能力强弱，引入ETOPO1 和DTU10 海深模型对本文模型进行了横向比对分析。结果表明，由于重力异常垂直梯度在较高频段所含“地形”信息相比重力异常更加丰富，从而研究海域以重力异常垂直梯度作为输入源构建的海底地形模型在浅海域整体质量优于以重力异常作为输入源构建的海底地形模型；同理，深海域重力异常建立海底地形效果较好。中位数方法、最小二乘方法和抗差估计方法构建的试验海域海底地形模型中，由于中位数方法强制回归分析过原点，忽略常数项的影响，检核精度最低；抗差估计构建的海底地形模型检核精度最高。最小二乘方法和抗差估计方法构建海底地形模型地貌呈现整体效果相当，而中位数方法构建的海底地形模型地貌细节呈现更为丰富。输入海底地形数据含粗差的情况下，抗差估计方法构建的海底地形模型检核精度变化最小，反映了利用抗差估计构建海底地形模型具有可靠性高、实用性强特点。

综上，依托回归分析方法恢复试验海域海底地形时，建议可根据实际需求，重点考虑使用重力异常垂直梯度，采用抗差估计方法构建海底地形模型。