海底地形和海底基岩密度可以反映出海底构造演化活动^[1]，可帮助探明海底矿产资源^[2]，此外，海底地形对于海事搜救等活动也是非常重要的，因此研究海底地形不仅有理论意义，而且有实用价值。目前，利用多波束回声测量装置直接测量海深是最好的方法，但是由于成本高、耗时长，因此，测得的海深数据比较稀疏，难以覆盖整个海洋^[3]。近年来，利用海洋重力数据反演海底地形的工作日益受到重视，其原因是成本低、覆盖范围广。从目前已有的研究来看，反演海底地形主要利用重力异常和垂直重力梯度异常。

重力地质法(gravity-geologic method，GGM)^[4]，其原理是将海面上重力异常中的短波部分与海底地形存在着线性关系，然后在实际使用时，需要结合船测的海深数据，拟合出该线性关系，从而得到海底地形。GGM方法目前已经被广泛地应用于海底地形的反演^[5-8]。事实上，GGM方法在海底较为平坦的区域是非常有效的，但对起伏较大的海底来说，反演海底地形的精度较差^[9-10]，因此GGM方法在应用上有一定的局限性。

频域方法源自Parker将引力位表达式进行傅里叶展开^[11]，给出了海面重力数据与海深之间的频率域表达式，略去海深的高阶小量后，即得频域内的线性关系，这为利用重力异常直接反演海底地形提供了理论依据。当考虑地壳的均衡补偿，Parker理论得到进一步完善，基于该理论，文献[12]反演皇帝海山链的海深，由于地形变化剧烈，海山链沿线附近精度不甚理想。文献[13]联合多源重力数据反演菲律宾海域的海底地形，精度提高30%左右。此外，文献[14]证实了在50~300km波段范围内，海底地形与卫星测高数据之间存在高度相关性。在文献[12, 14]的研究基础上，文献[15]说明了在15~160km的波段内，海底地形和重力异常间具有线性相关的特征。因为在应用Parker理论反演海底地形时都会舍去海深二阶以上的量，所以精确度存在不确定性^[16]。例如，文献[17]对比了GGM和Parker理论反演海底地形的结果后指出，GGM方法在可靠性和精度方面具有明显的优势。

作为重力数据的延伸，垂直重力梯度(vertical gravity gradient, VGG)对海底地形高频部分的敏感度超过重力本身^[18-19]。文献[20]推导出长方体产生的VGG公式，该公式已在地球物理勘探方面得到了许多应用。例如，文献[21]利用重力梯度全张量数据，采用共轭梯度法反演地质体与密度，但其仅仅限于模拟计算。文献[22]采用聚焦反演方法来预测异常地质体的埋深信息。此外，文献[23]也利用该公式，计算了不规则地质异常体产生的垂直重力梯度，并将其用于评估西墨西哥Jasper海山的密度。在将VGG应用于预测海深方面，文献[24]利用SIO(scrioos institude of oceangraphy)的VGG数据，并结合船测数据来反演海底地形，例如，在南印度洋海域的反演结果与船测深度的均方根误差约为188m。文献[25]利用VGG异常数据，顾及了Park理论中的二阶量影响和采用模拟退火法反演了海底地形，并给出其精度提高了22%的结论。

综合上述各种方法，其本质上还是利用经验公式来建立起重力数据与海底地形之间的线性关系，而这样的线性关系是随着海域的不同而变化的。为避开使用经验公式带来的影响，本文建立了一种解析的利用VGG数据来反演海底地形的方法，其原理是使用矩形分割将研究区域进行格网化，然后利用各个子矩形柱体产生的VGG表达式建立起各个分割子域的海深与该海域上VGG异常之间的函数关系，即导出了关于各个分割子域海深的观测方程组。但是在求解该观测方程组时，研究区域边界附近海深的求解易于受到区域外VGG异常的影响，即会产生边界效应，因此，如何评估边界效应的影响以及分离该影响的方法也就成了十分重要的研究课题。

为了减少边界效应的影响，本文将研究海域适当进行扩充，然后将扩充海域内的海深同时进行求解。为避免观测方程组出现奇异性，在求解观测方程组时引入了正则化参数，通过仿真计算后，确认正则化参数选取为10^-3时，反演的海域下方海深具有较好的精度。在对实际海域的海底地形进行反演计算时，除了研究海域下方海底地形会对VGG产生影响外，其他的海底地形也会对VGG异常产生影响，本文将扩充海域外的海底地形对VGG的影响统称为远区影响。为了消除远区影响，本文将其视为干扰常数ε，并在求解观测方程组时对其同步进行求解。

1 理论与方法 1.1 长方体产生的垂直重力梯度

假设有一个海底长方体海山Ω(图 1)，其长、宽、高分别为2a、2b、H-h，其中H为Ω底部至海平面的高度，h为Ω顶部至海平面的高度(即Ω的海深)。记R是Ω在海面对应的海域，显然确定Ω的形状仅需求解h。假设P是在海平面上任意一点，下面计算Ω在P处产生的VGG。

建立图 1所示的坐标系，则Ω在P点产生的重力位为

式中，G是万有引力常数(6.672×10^-11Nm²kg^-2)；ρ=2700kg/m³是海山Ω的密度；(x₁, y₁, z₁)为积分变量。因此，海山Ω在P点产生的VGG为

式中

积分区间为矩形区间，进一步计算后可得I_R(H, x, y)的表达式为

对于一般的矩形积分区间; ，结合式(4)，I_R(H, x, y)可写为

同理，式(2)中的I_R(h, x, y)与I_R(H, x, y)有相同的表达形式。因此，得到了矩形海山Ω在海平面上任何一点P产生垂直重力梯度的表达式。

1.2 观测方程与可解性分析

为了讨论方便，下面恒设R={(x, y); -a≤x, y≤a}是海面上一个正方形区域，Ω是R下方对应的海山，其密度为常数ρ，形状为Ω={(x, y, z); (x, y)∈R, H-h(x, y)≤z≤H}，其中h(x, y)是Ω的顶部至海面的高，因此，要确定R下方的海底地形就是确定函数h(x, y)。本节将在R上VGG异常仅是由Ω产生的前提下来研究如何求解h(x, y)。

选择一个步长t，将区间[-a, a]进行2N等分，即a=N·t，这样就能将区域R进行格网化，其中典型的子格网区域是。若步长t足够小，在R_ij中h(x, y)可以被看成是常数h_ij，于是，利用式(4)、式(5)和叠加原理可知，由Ω生成的VGG在海面P(x, y)处的值Γ(x, y)为

在式(6)中去掉海水的影响，则得Ω在海面上P处产生的VGG异常为

式中，Δρ=ρ-ρ₀，而ρ₀表示海水密度，在后面的计算中总是取ρ₀=1030kg/m³。

如果事先能够给出海域R上的VGG异常的观测值，例如在R的所有格网点上VGG异常是已知的，可得

式中，(x_p, y_q)是海面上区域R的观测点，p, q=-N, …, 0, …N。这样VGG异常仅由Ω产生的前提下，建立了关于h_ij的观测方程组(8)，其中h_ij就是子格网R_ij的平均海深。在式(8)中，方程的个数为(2N+1)²，未知数h_ij的个数为(2N)²，故满足可解条件。事实上，使用h_ij来替代h(x, y)的本质就是分块求出Ω的海深，即使用分片函数来表达Ω的海深。因此，当分割步长越小时，预测的Ω的地形就会越准确。由于式(8)关于h_ij是非线性的，因此，在给定初始值h_ij⁽⁰⁾后，可用迭代的方式来求解，即

式(9)就是线性化后的海底深度与VGG异常之间的观测方程组。

为了验证式(8)关于h_ij(或式(9)关于h_ij^(k+1))是可解的，本文通过一个算例来进行说明。假设在20km×20km海域R下有形如图 2(a)所示的海山Ω，海山的最大和最低高度分为0.3km和0.1km，以及H=2km。取步长t=2km，则N=10，海域R的格网化的结果如图 2(b)所示。设定容许误差是0.01m，而将初始值h_ij⁽⁰⁾限定在0与真值之间，采用最小二乘法来求解式(9)，经计算后可知，通过6次迭代计算就能达到设定的精度。由此说明式(9)具有很好的可解性。

1.3 观测方程组抗误差干扰分析

为了分析式(8)或式(9)的抗误差干扰特性，本文对上节算例中的观测值VGG异常分别加入随机误差和系统性误差：随机误差服从均值为0，标准差为δ Eötvös(1 Eötvös=10^-9s^-2)的正态分布；而系统误差是直接给观测值VGG异常加入相同的异常值。然后再利用式(9)来迭代求解h_ij。图 3(a)是针对随机误差的结果，其图像是反演结果的均方根(RMS)误差随标准差δ的变化；而图 3(b)则是系统误差的图像。由图 3(a)可知，随机误差的标准差增大时，反演结果的RMS误差也将变大，特别地当标准差δ为10 Eötvös时，反演结果的RMS误差为13.14m。而图 3(b)显示，反演结果的RMS误差与系统误差呈线性正比关系，特别当系统误差为10 Eötvös时，反演结果的RMS误差达到了76.3m。通过上述分析可知，式(9)具有很好的抗随机误差干扰性质。

2 边界效应的处理

上文建立的关于海深的观测方程式(8)或式(9)都是在海底仅有一个异常海山的前提下导出的，当处理实际海域时，显然海域下方的情况要复杂得多，特别是还存在着其他的异常海山。例如，图 4(a)模拟了40km×40km海域D下方的海山，如果仅需要求解中心20km×20km区域(称为内区R)下方的海底地形，显然在内区外的海山会对内区R上的VGG异常δΓ有贡献，因此，式(8)或式(9)需要进行调整。本节主要以图 4(a)模拟的海山为例来介绍求解方法和进行结果分析。

以后称40km×40km区域D为内区R的扩展区域，选取步长t=2km对D进行格网化(图 4(b))。称D-R为内区R的边界区域，其下方海山对内区R上VGG异常的贡献称为R的边界效应。假设在内区R上给出了VGG异常值(包含了边界效应)，本文的做法是：使用内区R上给出的VGG异常值来求解出整个扩展海域D下方的海底地形，然后剔除D-R下方的海山，仅取R下方海山作为最终的反演结果。注意：采用的已知数据仅是内区R上给出了VGG异常值。

事实上，在式(7)中将区域取为D，而将P(x, y)点仅限于R中，则有

这里当i, j=-10, …, 9时D_ij=R_ij。因为式(10)中待求解的量h_ij共有40²=1600个，所以在R中选择观测点如下：，其中p, q=-20, …, 20。于是，可得观测方程组如下

式中，p, q=-20, …, 20。因此方程组(11)共含有41²个方程。仿照式(9)，在(11)中对h_ij进行线性化，可得形如Ah=b的线性方程组，其中A是41²×40²矩阵，h是所有待解量h_ij排成的40²阶向量，而b则是41²阶常数向量。由于在(11)中观测点均在R中，所以线性化后的方程组Ah=b可能会产生奇异性。为了保证可解性，引入正则化方程如下

这里E是单位矩阵，α>0是正则化参数(其单位为10^-24m^-2s^-4)。

因为α>0，所以形如式(12)的方程关于h总是存在唯一解。因此通过迭代求解式(12)的方程组后，便能从式(11)解得所有的h_ij(i, j=-20, -19，…，18, 19)。从中选择出h_ij(i, j=-10, -9，…，8, 9)，便可得R的子分割R_ij的平均海深h_ij。关于正则化参数α的选择，本文以得到的R_ij下方平均海深h_ij的RMS误差作为标准，绘制了相应的曲线(图 5)。

由图 5可知，当正则化参数α=10^-3时海深h_ij的RMS误差达到最小，此时在内区R上恢复的海深的误差分布由图 6给出。由图 6可知，最大误差为0.64m，RMS误差为0.48m。因此通过将内区R扩充到D，并引入正则化参数，就可以有效地处理D-R下方海山引起的边界效应。

3 远区海山的影响

如果选定了内区R，通过适当扩充R得到扩展区域D，就可以求解R下方海底地形的边界效应。事实上，在扩展区域D外部的下方仍有海山会对R上的VGG值产生影响，本文称这样的影响为远区影响。尽管远区影响总体来看在R上是低频的，但其对R下方海底地形的反演到底影响到何种程度是需要探讨的。本节将在上节基础上，通过引入远区海山来进一步地分析远区影响及处理方法。

在上节扩展区域D的正左方再放置一个高为0.5km，长与宽则分别是40km和20km的海山Ω，其底部至海面深度仍取H=2km。利用式(4)极易算出Ω生成的VGG异常。图 7绘制了沿y=0，从D的左边起算的Ω生成的VGG异常图像。由图 7可知，在D的左边Ω生成的VGG异常为8.3 Eötvös，进入内区R后该VGG异常从-5.8 Eötvös变化至-0.75 Eötvös，之后该VGG异常趋于零。

如果将远区影响看成是R上VGG异常的系统误差，则可直接使用本文处理边界效应的方法来反演R下方海底地形。增加远区影响Ω后，其反演结果的误差分布由图 8(a)给出。从图 8(a)可知，受Ω影响，反演的内区R下方的海底地形的最大与RMS误差分别是71.8m与39.4m。

由图 7可知，由Ω生成的VGG异常在R上的绝对值是介于0.73与5.8 Eötvös之间的。事实上，若将远区Ω在R上生成的VGG异常看成是R上VGG异常的系统误差的话，则该系统误差是介于0.73~5.8 Eötvös之间。进一步地将图 8(a)与图 3(b)进行比较，在图 3(b)中，对应于系统误差为5.8 Eötvös时，其反演海深的RMS误差为35m；而图 8(a)说明系统误差介于0.73至5.8 Eötvös之间时，反演海深的RMS误差达到了39.4m。原因就是所解算的问题是不一样的，图 8(a)是基于边界效应求解的，而图 3(b)是不存在边界效应的。由此也可以得出这样的结论：有了边界效应后，反解内区下方海底地形的抗误差干扰性会降低。

在处理实际海域问题时，一旦选定了内区，其远区影响是非常大的，这是因为实际的远区影响是由多种因素产生的。例如：除了本文描述的扩展区域外海山的影响外，还存在着地球深部质量分布不均匀等因素的影响。如果相对于内区R来说，远区影响太大，则观测方程组(11)中会出现很大的误差项，这将会导致反演结果完全不正确。因此，必须选择某种途径去处理远区影响。

尽管远区影响的量级会很大，但是在R上VGG异常会表现出低频的特征，其理由是远区影响距内区R相对较远，而VGG自身在离开异常体后会表现出衰减很快。因此，本文将引入一个常数ε来表示R的远区影响，将式(11)改写为

然后在式(13)中同时对h_ij和ε进行求解。至于如何求解式(13)，仍然采用对h_ij进行线性化，再使用正则化参数的方法来迭代解算h_ij和ε。

按照上面的方法，对Ω的远区影响问题进行了计算，反演得到了内区R下方的海底地形，其误差分布由图 8(b)给出，其中最大与均方根误差分别是58.7m与32.5m。相较于图 8(a)的结果，最大误差和均方根误差分别降低18.2%和17.5%。这说明加入远区影响项ε后，确实是可以消除部分远区影响的。

4 实际算例及分析

前面各节分别讨论了单个海山、周边海山(边界效应)、远区海山(远区影响)等情况下如何反演某研究海域R下方海底地形的反演方法。本节将选取某个实际海域作为研究对象，来反演其下方的海底地形。这里选取南中国海72km×72km的海域(112.6°E~113.25°E, 11.0°N~11.65°N)作为研究对象，密度差Δρ=ρ-ρ₀选取为1670kg/m³。此外，根据已有的海深数据，选择H=4.4km。

首先在72km×72km的海域D中选择内部52km×52km海域R作为内区，其上的VGG异常数据(图 9(a))来自GFZ。然后选择步长t=2km，由此便可以建立方程式(13)，对h_ij进行线性化，引入正则化参数α=10^-3，迭代后便能解出R的子分割R_ij的平均海深h_ij。最后将海深h_ij置于R_ij的中心，便可生成R下方的海底地形如图 10所示。

为了检验反演出的R下方的海底地形的精度，本文选取了该海域全部289个来自NGDC(National Geophysical Data Center)船测点的海深数据，其分布由图 9(b)给出。将本文反演结果与NGDC数据进行比较后可知：其最大、最小和均方根误差分别为339、-157和109m。统计反演结果的平均深度为4173m，均方根误差相当于平均深度的2.6%(109m/4173m)。另外，289个点上的误差绝对值统计见表 1，结果表明：误差主要集中在150m以下，占比为84.1%，误差大于150m的点，占比为15.9%，且这部分点主要分布于R的边缘处。如若将最终的反演结果再向内收缩一个格网，即2km，反演结果的最大误差和均方根误差将分别下降为207m和99m，同比分别下降38.9%和9.2%。

5 总结与讨论

目前，研究重力场与海底地形的关系主要还是使用经验公式，即，通过部分实测的海深数据来拟合出海底地形与测高、重力和重力梯度等数据之间的线性关系，以此来进一步地推算海域的海底地形。显然，这样的算法是随着海域的不同而变化的。本文工作则是依靠建立起的海底地形与VGG异常之间的严格函数关系(或观测方程组)而展开的，理论上看，本文工作是解析的方法。本文：①通过在观测方程组中加入随机误差和系统误差后，进行模拟计算，结果表明观测方程组具有极强的抗误差干扰特征；②为了处理边界效应，引入了扩充区域，并在求解线性化方程组时，引入正则化参数，确保了线性化方程组的可解性，同时经模拟计算表明由此解算的海底地形具有很高的精度，统计出的均方根误差为0.48m；③分析了远区影响对反演结果的影响，使用某个常数来吸收远区影响，这样就能确保远区影响较大时，观测方程组依然有很好的可解性；④将本文的理论方法用于实际海域中，并得到统计的均方根误差为109m的海底地形反演结果，从而验证了本文方法的有效性。

本文对海域格网化时，采用的步长为2km，即在2km×2km范围内，仅使用平均海深来表示海山，这意味着在2km范围内海底地形的起伏是难以分辨出来的。如果将步长取得更小，理论上讲对海山的解算精度也会增加，但是重力场数据的分辨率会制约步长的选择。例如，EGM2008重力场模型的阶次是2160^[26]，其空间分辨率约为8km，如果步长小于4km，那么来自EGM2008模型的重力数据是不能反映出4km×4km内的地形起伏的。因此，处理实际问题时，步长的选择与重力数据的分辨率是相关的。基于此，若要提高海底地形的反演精度，不仅要提高重力数据的精度，而且还需要提高数据的分辨率。

本文主要的研究内容还是侧重于理论与方法上，即针对建立的从垂直重力梯度数据反演海底深度的观测方程组，进行抗误差性分析以及处理边界效应等分析、验证。因此，在对实际海底地形进行预测时，都将海山和海水密度取成了常用的推荐值。事实上，在实际预测海底地形时，应该尽可能地参考当地的地质与海洋资料来选择海山与海水的密度分布。

致谢: 感谢GFZ和NGDC提供的计算数据。