重力反演是当前全球海底地形模型构建的主要手段。国内外发布的全球海底地形数据资料中，绝大部分海域的数据均依赖重力数据恢复得到，如STO_IEU2020^[1]、BAT_WHU2020^[2]和GEOCO系列模型^[3]等。利用重力数据反演海底地形的方法主要包括重力地质法(gravity-geologic method，GGM)^[4-5]、导纳函数法^[6-7]、线性回归分析法^[8]、SA(simulated annealing)法^[9]、解析法^[10-11]等。应用重力反演海底地形的研究主要集中在规则矩形海域^[12-15]，鲜有涉及海陆交界区非规则形状海域。究其原因，一是导纳函数法和SA法一般在频率域处理重力数据和地形数据，涉及傅里叶正反变换操作，所以通常应用于规则矩形海域海底地形反演，对海陆交界区非规则形状海域海底地形恢复效果不佳；二是以上方法的实施一般涉及数据格网化，而海陆交界区非规则形状海域数据格网化结果势必受陆地区域数据干扰，导致常规海底地形反演方法恢复的数据结果鲁棒性较弱。可以说，利用重力数据反演海陆交界区非规则形状海域海底地形是目前重力反演海底地形研究的难点之一。然而，一方面海陆交界区海底地形是自主建立全球海底地形模型不容忽视的重要内容；另一方面，随着临海区域多源重力数据不断累积和数据处理技术不断进步，海陆交界区海洋部分相关可用数据质量将得到极大提升，可以预见，临海区域海底地形相关数据的开发与利用势在必行。

综上，针对利用重力反演海陆交界区非规则形状海域海底地形算法实施困难的问题，基于海底地形数据与重力数据间近似线性的关系，提出使用改进的GGM和回归分析方法，联合水深测量数据和重力异常数据构建海陆交界区非规则形状海域海底地形，并以中国南海部分海陆交界区海域为研究对象，融合SIO V31.1重力异常数据和声学水深数据，构建相应海陆交界区非规则形状海域海底地形模型。以研究海域实际船载声学水深数据作为外部检核条件，将重力反演海底地形结果与研究海域声学水深数据直接格网化结果、topo_23.1海深模型结果等进行比对，评估本文方法恢复的非规则形状海域海底地形数据质量。研究成果对融合重力数据的海陆交界区非规则形状海域海底地形资料获取，以及全球海底地形模型精细化构建有一定参考价值。

1 原理与方法 1.1 改进的重力地质方法

GGM反演海深的基本原理是将海面重力异常划分为长波参考重力异常和短波残余重力异常两个部分。长波参考重力异常与地壳下层物质质量有关，短波残余重力异常主要由海底地形变化引起。3种重力异常的关系为：

$ \Delta g=\Delta g_l+\Delta g_s $

(1)

式中，Δg为海面重力异常，Δg_l为长波参考重力异常，Δg_s为短波残余重力异常。研究点短波残余重力异常与海深的关系可表示为：

$ \Delta g_s^i=2 \pi G \Delta \rho\left(d_i-d\right) $

(2)

式中，G为万有引力常数，Δρ为研究点地壳密度和海水密度之差，Δg_sⁱ为研究点i处的短波重力异常，d_i为研究点i处水深，d为研究点参考水深(通常取研究区域最大水深)。依据式(2)可得研究点水深为：

$ d_i=\frac{\Delta g_s^i}{2 \pi G \Delta \rho}+d $

(3)

由此可以看出，GGM反演海底地形的质量与Δρ密切相关。通常以检核均方差等指标作为评价依据，不断调整、优化地壳和海水间密度差异常数，进而获取目标海域海底地形建模对应的最佳密度差异常数^[16]。范雕^[16]研究认为，密度差异常数的取值只是为了让重力异常的短波分量与海深保持一种线性关系，即让重力异常的短波分量与海深之间达到平衡。常规GGM通常仅可恢复规则矩形形状海域海底地形，而海陆交界区海底地形模型构建涉及不规则形状海域海底地形。因此，本文提出改进的GGM(称为iGGM)构建海底地形策略，以恢复海陆交界区海底地形数据，具体技术路线见图 1，主要步骤如下：

1) 将声学水深测量数据按照一定比例分为两类，一类为水深离散控制点，参与海底地形模型构建；另一类为水深离散检核点，用于海底地形模型质量评估，以确定最佳密度差参数。

2) 设置密度差参数初值，利用式(2)获得离散控制点处重力异常短波分量。

3) 基于目标海域海面重力异常和离散控制点处重力异常短波分量，采用式(1)计算得到离散控制点处重力异常长波分量。

4) 采用相应的格网化方法将控制点处重力异常长波分量格网化，获得目标区域格网化长波重力异常。

5) 使用相应的内插方法，将目标区域格网化长波重力异常内插，获得海陆交界区海洋部分网格长波重力异常。

6) 基于目标海域海面重力异常和海陆交界区海洋部分网格长波重力异常，解算得到海陆交界区海洋部分网格重力异常短波分量。

7) 利用海洋短波重力异常与海底地形间函数关系，反演得到海陆交界区海洋部分网格化海底地形，并使用事先预备的水深离散检核点对构建的海底地形模型进行精度评估。

8) 在设置的密度差异常数初值基础上，以一定步长调整密度差异，重复步骤2)~8)。

9) 以计算的海底地形模型检核统计数值结果为依据，将海底地形模型检核均方差最小作为确定最佳密度差的关键参考，相应密度差常数反演的海陆交界区海洋部分网格化海底地形即为海陆交界区域海底地形数据结果。

对比iGGM和常规GGM构建海底地形的策略可知，iGGM主要改进之处为：基于海陆交界非规则形状区海洋网格排布，使用相应内插方法获取非规则形状海洋区域短波重力异常，进而依据短波重力异常与海底地形间函数映射关系恢复非规则形状区海洋网格海底地形数据。

1.2 改进的回归分析方法

重力反演海底地形的导纳函数理论表明，经过滤波且向下延拓处理后的重力异常与滤波后海底地形具有较强的线性相关性^[17]，理论上，比例系数S=(2πGΔρ)^－1，即

$ \begin{gathered} F(b(x, y)) W\left(f_x, f_y\right)=\frac{1}{2 \pi G \Delta \rho} \cdot \\ \left\{F(\Delta g(x, y)) W\left(f_x, f_y\right) \mathrm{e}^{2 \pi f d}\right\} \end{gathered} $

(4)

式中，F(Δg(x, y))和F(b(x, y))分别为重力异常和海底地形的二维傅里叶变换，G为地球引力常数，Δρ为海水和地壳的密度差异，d为平均海深，f为频率，$f=\sqrt{f_x^2+f_y^2}$f_x和f_y分别为x和y方向的频率，W(f_x, f_y)表示为削弱重力数据高频噪声和地壳均衡补偿的低频影响所研制的带通滤波器，滤波器参数设置见文献[16]。

由式(4)可知，比例系数与相应频段重力异常的乘积为有限波段海底地形结果^[17]。回归分析方法构建海底地形结果可表示为：

$ h(x)=h_r(x)+S(x) y(x) $

(5)

式中，y(x)为经过延拓和滤波处理的重力异常，S(x)为海底地形与重力异常比例因子，h_r(x)为非反演波段海底地形，通常为船测水深输入数据通过预置的高通和低通滤波器获得的短波段海底地形与长波段海底地形之和^[17]。

全海洋覆盖区规则矩形海域海底地形模型构建一般将控制点处比例因子格网化^[18]或者使用移动窗口内比例因子^[19]。本文海陆交界区海底地形模型构建涉及不规则形状海域海底地形，为规避海陆交界区域比例因子格网化过程中陆地因素的影响，使用回归分析方法获取的是整个目标海域海底地形-重力数据比例因子，即整个目标海域共用同一个比例因子。海陆交界区域海底地形模型构建的改进回归分析方法(简称iRA)具体技术路线见图 2，实施步骤如下：

1) 将目标海域水深测量数据和海面重力数据进行格网化和滤波操作，分别获取非反演频段海底地形、反演频段海底地形和反演频段重力异常。

2) 基于海陆交界区海洋网格定位排布规律，以步骤1)获得的非反演频段海底地形和反演频段重力异常为输入，使用相应的内插方法获得海陆交界区海洋网格非反演频段海底地形和海陆交界区海洋网格反演频段重力异常。

3) 采用相应的内插方法，将反演频段海底地形和反演频段重力异常转换至海陆交界区海洋部分离散控制点处，获得离散控制点处反演频段海底地形和反演频段重力异常资料。

4) 以步骤3)得到的离散控制点处反演频段海底地形和反演频段重力异常资料为输入，使用回归分析技术得到目标海域海底地形-重力数据比例因子。

5) 步骤4)解算的海底地形-重力数据比例因子和步骤2)海陆交界区海洋网格反演频段重力异常的乘积即为海陆交界区海洋网格有限频段海底地形。

6) 步骤5)获得的海陆交界区海洋网格有限频段海底地形与步骤2)海陆交界区海洋网格非反演频段海底地形叠加，即为海陆交界区海洋网格化海底地形结果。

2 实验分析 2.1 实验区域概况与数据准备

选择中国南海海陆交界区域2°×2°(6°~8°N, 115°~117°E)范围为研究区。SIO UCSD(Scripps Institution of Oceanography, University of California, San Diego)于2021-11发布的分辨率为1′×1′的topo_23.1海深模型在实验区域的空间呈现见图 3，本文称为SCS_topo_23.1模型。实验海域开展数值分析的相关数据来源如下：1)海面重力异常数据来自SIO UCSD于2022-02发布的grav_31.1，为描述方便，研究区域重力异常称为SCS_grav_31.1(图 4)。2)船测多波束和单波束海深数据来自NOAA(National Oceanic and Atmospheric Administration)和NGDC(National Geophysical Data Center)发布的原始船测水深数据。

鉴于船测水深数据年代较为久远，本文选择SCS_topo_23.1海深模型为参照对其进行粗差处理，具体处理策略见文献[16]。剔除研究区域原始船测水深数据粗差后，得到15 424个“干净”的实测海深数据。均匀选择其中约20%的数据作为检核点，其余作为控制点，最终得到3 086个检核点和12 338个控制点，控制点和检核点的分布见图 5(黄色圆点为控制点，紫色三角形为检核点)。

2.2 海陆交界区海底地形模型构建

由于陆地因素影响，海陆交界区海洋部分并非规则矩形，利用GEBCO_2021海深模型获取实验海域海陆交界区海洋部分1′网格的定位排布规律，具体方法为：提取GEBCO_2021海深模型中海洋部分数据格网位置，因GEBCO_2021海深模型格网大小为15″，所以提取的海洋格网边长也为15″，采用滑动窗口方法将边长大小为15″的格网展布为边长为1′的离散网格形式，进而完成海陆交界区海洋部分边长为1′的离散网格定位，最终得到实验海域共计12 293个边长为1′的离散网格。以实验海域SCS_grav_31.1重力异常和船测水深控制点数据为建模输入数据，使用iGGM和iRA方法初步反演得到实验海域海底地形结果，分别称为SCS_iGGM0模型和SCS_iRA0模型。为分析重力数据在海陆交界区海底地形建模中的作用，仅以船测水深控制点数据为输入，使用格网化方法恢复海陆交界区海底地形资料，称为SCS_Grid0模型。SCS_iGGM0、SCS_iRA0和SCS_Grid0模型分别如图 6(a)、6(b)和6(c)所示，图中红色区域表示海底地形恢复失败(恢复的海深值为非负值)。

对比图 6和图 3可以看出，恢复失败的地方多为实验海域靠近陆地的部分。原因可能有两点：一是依据图 5可知，实验海域靠近陆地部分船测水深点稀少，存在大面积数据空白，从而该区域缺少必要的水深控制，导致恢复效果欠佳，特别是对于仅依靠船测水深构建的SCS_Grid0模型尤为明显；二是SCS_grav_31.1重力异常数据是通过卫星测高技术解算得到的，而卫星测高数据极易受陆地干扰，致使靠近陆地附近的卫星测高数据污染程度较高，所以基于卫星测高重力数据恢复的沿陆区海底地形效果欠佳。3种模型恢复的实验海域网格数统计结果见表 1。

表 1显示，SCS_iGGM0模型恢复10 958个网格数据，SCS_iRA0模型恢复10 440个网格数据，二者恢复的海底地形数据面积分别占实验海域总面积的89.14%和84.93%。说明融合重力数据反演海陆交界海域海底地形时，iGGM方法优于iRA方法。另外，仅依靠水深控制点的SCS_Grid0模型恢复的海底地形数据占实验海域总面积的83.74%，低于iGGM方法和iRA方法恢复的海底地形面积。综上，SCS_iGGM0模型对海陆交界区海底地形的恢复能力最强，其次是SCS_iRA0模型，SCS_Grid0模型最弱。

为完全恢复实验海域海底地形数据资料，使用最邻近内插/推估方法将SCS_iGGM0、SCS_iRA0和SCS_Grid0模型海深值扩充至整个海域，结果分别称为SCS_iGGM、SCS_iRA和SCS_Grid模型，分别见图 7(a)、7(b)和7(c)。引入DTU18海深模型和ETOPO1海深模型进行横向比对，二者在实验海域分别称为SCS_DTU18和SCS_ETOPO1模型，分别见图 7(d)和7(e)。从图 7(a)~7(e)可以看出，5种海深模型均能较好地描述研究海域海底地形地貌形态，但其中SCS_DTU18模型西北方向相较于其他4种模型存在4处明显的地形隆起。

2.3 海底地形模型精度分析评估

SCS_Grid、SCS_iGGM、SCS_iRA、SCS_DTU18、SCS_ETOPO1和SCS_topo_23.1模型数值特征统计见表 2(单位m)。

由表 2可见，采用iGGM和iRA构建的海深模型水深最浅处近乎0，最深处约3 000 m；水深平均值和标准差分别约为930 m和1 000 m。SCS_iGGM和SCS_iRA模型最大值、平均值和标准差等数值特征与其他4种模型接近。由此可以看出，采用本文提出的iGGM和iRA构建的海陆交界区非规则形状海域海底地形模型结果可靠。比较各海底地形模型间的差异，统计结果见表 3(单位m)。

由表 3可见，SCS_iGGM、SCS_iRA、SCS_Grid模型与SCS_DTU18模型间差值标准差约为145 m，与SCS_ETOPO1模型差值标准差约为105 m，与SCS_topo_23.1模型差值标准差约为115 m。相较而言，SCS_iGGM、SCS_iRA、SCS_Grid模型与SCS_ETOPO1模型更为接近。模型间差值平均值表明，SCS_ETOPO1模型与本文构建的3种海底地形模型间系统误差最大，约54 m；SCS_topo_23.1模型与本文构建的3种海底地形模型间系统误差最小，为23 m左右。以SCS_topo_23.1模型为参照，SCS_iGGM、SCS_iRA、SCS_Grid模型与其差值标准差分别为111.48 m、113.66 m、120.87 m。由此可知，本文构建的3种模型中，SCS_iGGM与SCS_topo_23.1模型更为接近，SCS_Grid与SCS_topo_23.1模型符合程度较弱。以SCS_ETOPO1和SCS_DTU18模型为参照，所得结论相同。综上表明，融合重力数据的海陆交界实验区海底地形模型构建效果相较于船测水深数据直接格网化建模更优。

以3 086个检核点作为外部检核条件，将上述6种海深模型内插到检核点，内插结果与检核点处船测海深值作差，统计结果见表 4(单位m)。

由表 4可见，SCS_iGGM、SCS_iRA、SCS_Grid模型检核差值RMS和标准差均为40 m左右；SCS_DTU18、SCS_ETOPO1、SCS_topo_23.1模型检核差值RMS分别为183.79 m、138.15 m、87.04 m，检核差值标准差分别为174.04 m、121.45 m、84.60 m。就检核差值RMS或者检核差值标准差而言，SCS_iGGM、SCS_iRA、SCS_Grid模型检核精度相比SCS_DTU18、SCS_ETOPO1、SCS_topo_23.1模型分别提高约76%、70%、53%。另外，SCS_DTU18、SCS_ETOPO1、SCS_topo_23.1模型检核差值平均值分别为59.07 m、65.84 m、20.49 m，说明3种模型分别存在不同程度的系统误差；而SCS_iGGM、SCS_iRA、SCS_Grid模型系统误差近乎0。同时，本文海底地形模型检核差值最值统计结果也优于SCS_DTU18、SCS_ETOPO1、SCS_topo_23.1模型。比较SCS_iGGM、SCS_iRA、SCS_Grid模型检核差值统计结果可知，SCS_iGGM模型略好于SCS_iRA模型，二者均优于SCS_Grid模型。也就是说，实验海陆交界海域顾及重力贡献的海底地形建模相较于仅依靠船测水深构建的海底地形模型效果更佳。SCS_iGGM模型外符合精度优于SCS_iRA模型，原因在于，SCS_iGGM和SCS_iRA模型虽均基于海底地形与重力数据间线性关系构建，但SCS_iGGM模型通过持续调整和优化实验海域密度差常数使得海底地形与重力数据达到最佳平衡状态，而SCS_iRA模型则依据回归分析方法一次性拟合实验海域海底地形与重力数据比例因子，并未对比例因子这一关键因素再次进行优化。

3 结语

本文提出使用改进的GGM和回归分析，联合水深测量数据和重力异常数据的海陆交界非规则形状海域海底地形构建方法。选择中国南海部分海陆交界区海域作为实验区，利用SIO V31.1重力异常分别构建SCS_iGGM和SCS_iRA模型，并与船测水深数据直接格网化结果、SCS_topo_23.1等海深模型进行对比分析，主要得到以下结论：

1) 外部检核结果表明，SCS_iGGM和SCS_iRA模型精度相当，明显优于SCS_DTU18、SCS_ETOPO1和SCS_topo_23.1模型；相较于SCS_DTU18、SCS_ETOPO1和SCS_topo_23.1模型，基于改进的GGM方法和回归分析方法构建的海底地形模型精度分别提高约76%、70%和53%。

2) 相比仅依靠船载水深数据恢复的海底地形SCS_Grid，联合重力异常和船载水深数据建立的海底地形模型SCS_iGGM和SCS_iRA的数据恢复能力和模型检核质量均更好。也就是说，针对船测水深数据稀疏的海域，联合重力数据的海底地形模型效果更佳。

值得说明的是，目标海域水深控制点分布的均匀程度也会影响SCS_iGGM、SCS_iRA、SCS_Grid三种模型的效能，其具体影响有待进一步研究。