海底地形模型在海洋学、固体地球物理学等研究领域具有重要意义。采用装载在船舶上的测深设备进行测量是获得海深的直接手段。但是，经过几十年的积累，当前船测海深数据在大洋上的覆盖依然非常稀疏，且大量数据是在卫星导航技术出现之前测量的，存在较大定位误差^[1]。当前，多波束测深系统等新技术虽然能够获得高精度、高分辨率的海底地形，但由于船舶航行速度很慢，要采用现代船测手段完成全球大洋海底地形的测绘，预计需要100~200 a时间和几十亿美元资金^{[2, 3]}，浅海区域地形测绘需要更长时间^[4]。因此，在短期内，不可能采用船测技术构建全球海底地形模型。

卫星测高技术的发展，为海底地形研究提供了全新技术手段。文献[5]首先证明了采用卫星测高数据预测海底地形的可行性。此后，海底地形反演技术渐趋成熟，所采用的数据一般为船测海深和测高重力异常，所采用的方法可分为解析法和统计法两大类^[6]。统计法由文献[7]首次采用，其后有研究者对该方法提出过一些改进^[8]，但由于海底地形与重力异常之间的互协方差函数难以确定，该方法的应用较少。解析法根据地壳均衡理论，直接建立海底地形与海面重力异常之间的函数关系，再进行海底地形反演计算，依据具体数据处理方法的不同，又可分为空域最小二乘法^{[9, 10]}和频域法^{[11, 12, 13]}。空域最小二乘法在计算过程中需涉及洋壳厚度、岩石圈挠曲强度、洋壳及地幔密度等地球物理参数，且计算速度慢，因而较少采用。

基于频域法，利用测高重力异常反演海底地形时，通常受地壳均衡、沉积层等影响较大，而垂直重力梯度异常数据可削弱这些异常物质的影响。文献[14]通过对高斯型海山模型的研究指出，垂直重力梯度异常具有放大短波信号、抑制长波信号的作用，岩石圈有效弹性厚度和莫霍面以下的密度异常对垂直重力梯度异常的影响，远小于对重力异常和大地水准面的影响，并据此研究了太平洋海山分布状况。文献[15]提出了根据重力垂直梯度异常反演海底地形的方法，推导了直角坐标系下的相应数学计算模型，并指出根据垂直重力梯度异常反演海底地形的优势在于能够抑制地壳均衡及沉积层密度异常等中长波信号误差的影响。此后，文献[16]根据这一原理，采用FFT方法反演计算了南中国海海底地形，但是由于垂直重力梯度异常数据高频噪声很大，反演结果不理想。文献[17]采用模拟数据，研究了根据垂直重力梯度异常数据反演海底地形时地壳密度异常、岩石圈有效弹性厚度等参数的影响。目前，在国内外，鲜有根据垂直重力梯度异常数据反演海底地形的成果发表。

本文研究联合船测海深、重力异常和垂直重力梯度异常数据，构建中国海及其周边1′×1′海底地形模型。利用三种数据的相对优势，采用稀疏分布的船测海深数据构建波长大于200 km的海底地形，采用垂直重力梯度异常数据反演100~200 km波段内的海底地形，以抑制地壳均衡、沉积层等干扰因素的影响，采用重力异常数据反演波长小于100 km的海底地形，克服垂直梯度异常反演时高频噪声过大的问题。以未参与海底地形计算的船测海深为检核数据，并与ETOPO1、DTU10、GEBCO和斯克里普斯海洋研究所(Scripps Institute of Oceanography,SIO) V15.1模型进行比较，分析了反演结果精度。

本文采用的数据有来自美国国家地球物理数据中心(National Geophysical Data Center,NGDC)的船测海深数据，共2 483 805个(见图 1)；来自SIO的测高重力异常数据(见图 2)和垂直重力梯度异常数据(见图 3)，版本均为V20.1^[18]。为尽量保证船测数据对长波地形的良好控制，将其中95％的数据用于海底地形模型构建，另5％用于反演结果的精度检核。

图 1 船测数据分布 Fig. 1 Distribution of Ship Soundings

图选项

图 2 测高重力异常 Fig. 2 Altimetric Gravity Anomalies

图选项

图 3 重力垂直梯度异常 Fig. 3 Vertical Gravity Gradient Anomalies

图选项

从图 1看，船测数据主要密集分布于日本列岛周围、琉球海沟、马尼拉海沟等区域，在开阔海域，尤其是在菲律宾海南部，船测海深数据分布非常稀疏，不能仅依赖于船测数据构建较高分辨率的海底地形模型。图 2和图 3所示的由卫星测高技术获得的数据，在整个区域内具有较为密集而均匀的覆盖，可以作为船测数据的有效补充，用于构建分辨率和精度均较高的海底地形模型。

本文采用解析法中的频域法进行数据处理，其地球物理、数学基础是岩石圈挠曲均衡理论^[19]和Parker公式^[20]。根据岩石圈挠曲均衡理论，海面上的重力异常和重力垂直梯度异常的产生，主要源于海底地形起伏及其在莫霍面下的均衡补偿物质(洋壳挠曲)的质量异常。Parker则推导了由海平面以下的海底起伏，计算其在海面产生的重力异常的频域级数展开式：

式中，ΔG(k)为重力异常的傅里叶变换；G为万有引力常数，ρ_c、ρ_ω分别为洋壳、海水密度；k=2π/λ为波数；λ为地形波长；d为平均水深；H(k)为海底地形起伏的傅里叶变换；F表示傅里叶变换；h(x)为空域内海底地形起伏。

根据挠曲均衡模型，可计算均衡补偿物质产生的重力异常。首先，计算海底地形加载情况下的莫霍面起伏(洋壳挠曲量)：

式中，R(k)为地壳挠曲量的傅里叶变换；ρ_m为地幔密度；Φ_e(k)为岩石圈挠曲响应函数，有：

式中，g为平均重力加速度；D为岩石圈抗挠刚度；D=ET_e³/[12(1－υ²)]，E为杨氏模量，一般取10¹¹N/m²，T_e为岩石圈有效弹性厚度，υ为泊松比，一般取0.25。

当同时顾及海底地形起伏及洋壳挠曲时，根据式(1)，给出海底地形起伏与海面重力异常的解析函数关系为：

式中，t为平均洋壳厚度，其他符号意义同前。根据傅里叶变换导数定理，直接给出海底地形起伏与海面垂直重力梯度异常之间的函数关系为：

式中，ΔG_z(k)为重力垂直梯度异常的傅里叶变换，其他符号意义同前。

式(4)、(5)即分别是根据重力异常和垂直重力梯度异常反演海底地形的理论基础，公式中涉及洋壳平均厚度、岩石圈有效弹性厚度等未知参数，导致直接应用上述两式时存在一定困难。Smith和Sandwell指出，在海底地形构建过程中，为抑制均衡效应的影响，可以采用船测海深数据控制长波海底地形(>200 km)，而仅用重力异常反演短波段的海底地形，此时，式(4)、(5)中的均衡项忽略不计。Parker的研究表明，式(4)、(5)收敛速度很快，为简化计算过程，可以直接忽略两式中n=2以上高次项。因此，式(4)、(5)可以简化为：

当海面重力异常和重力异常垂直梯度已知时，可根据式(6)反算相应的海底地形：

从式(7)看，在短波反演波段内，海底地形与经向下延拓等处理后的重力异常(或垂直重力梯度异常)呈线性关系。实际数据处理过程中，由于存在沉积层等密度异常的影响，式(7)中的线性系数不根据理论密度参数直接计算，而是在船测海深数据约束下，分别计算不同区域的线性系数。

总之，联合船测海深、重力异常和垂直重力梯度异常构建海底地形模型的数据处理流程如图 4所示。

图 4 联合多源数据构建海底地形模型的数据处理流程图 Fig. 4 Data Processing Diagram for Bathymetry Predicting from Multi-source Data

图选项

从图 4看，海底地形构建的数据源为船测海深、垂直重力梯度异常和重力异常。首先，联合船测数据和垂直重力梯度异常数据构建一个初步海底地形h_grad(x)，其简要数据处理流程包括两步：① 将船测海深格网化，并进行200 km低通滤波处理，获得长波海底地形模型h₂₀₀(x)；②对垂直重力梯度异常数据进行20~200 km的带通滤波处理(20 km的短波端截断是笔者经数值计算实验获得的较优值，如图 5所示)，然后进行海底地形反演计算，并将反演结果与h₂₀₀(x)相加，即可获得初步海底地形模型h_grad(x)。其次，以h_grad(x)为参考模型，对该模型进行100 km低通滤波，得到长波海底地形h₁₀₀(x)，再对重力异常数据进行100 km高通滤波处理，并反演短波海底地形h_grav(x)，将该结果与h₁₀₀(x)相加，即得最终反演结果hpredicted(x)：

图 5 低通滤波器设计 Fig. 5 Low-Pass Filter

图选项

式中，h₁₀₀(x)包含船测数据控制的长波海底地形h₂₀₀(x)和根据垂直重力梯度异常数据反演计算的100~200 km波段内海底地形。 为抑制重力异常数据的高频噪声，对短波端进行了低通维纳滤波处理：

式中，A为滤波器设计参数；SIO的学者在根据最新重力异常数据V20.1构建海底地形模型V15.1的过程中，取A=3 891 km⁴。此时，平均水深分别为2、4、6 km时对应的截断波长(使W(k)=0.5)分别为12、16、20 km，如图 5所示。

在短波端，对重力异常和垂直重力梯度异常数据进行低通滤波处理，主要是考虑到向上延拓效应的影响，使得观测数据在短波段信噪比较低。经图 5所示的低通滤波处理后，最终海底地形模型的真实分辨率即可认为是此滤波器的截断波长。

根据第§1给出的数据和计算原理，构建了中国海及周边的海底地形模型，大于200 km长波海底地形、100~200 km垂直重力梯度异常反演结果、小于100 km波段重力异常反演结果，以及最终模型结果，如图 6所示。

图 6 联合多源数据构建的中国海及邻区海底地形模型 Fig. 6 Bathymetry of China Sea and Its Adjacent Areas Predicted from Multi-source Data

图选项

根据文献[21]的建议，在船测海深数据分布密集的地区，可以对海底地形反演结果进行进一步校正，其方法是：根据反演模型内插船测点海深，再将此内差值减去船测值获得船测点上的残差值；然后，将残差值格网化获得残差模型；最后，在反演模型中减去此残差模型，即得最终反演结果(见图 6(d))。

采用未参与计算的船测海深数据对反演结果进行精度检核，并以ETOPO1、DTU10、GEBCO和V15.1模型为对比参考，评价反演结果的精度。在南中国海地区，检核点模型值与船测值之差统计见表 1，残差值的频率分布直方图见图 7。

图 7 南中国海检核点残差分布比率直方图 Fig. 7 Histogram of the Residual Depth on the Checking Points in South China Sea

图选项

从表 1看，在南中国海，本文计算的模型精度为136.802m，优于由纯船测数据构建的模型GEBCO以及目前被广泛采用的模型ETOPO1，略优于SIO的最新海底地形模型V15.1。表 1中，个别检核点的残差绝对值达千米以上，但从图 7看，残差值的分布非常集中，主要在±50m以内，占82.5％，在±200m以内的比率达到95.8%。

菲律宾海的精度检核见表 2和图 8。从表 2可知，本文构建的模型精度为120.85m，优于当前常见的海底地形模型GEBCO、ETOPO1和DTU10模型，并且较来自SIO的最新模型V15.1的精度提高了约8.2％。从图 8看，残差值的分布非常集中，±50m以内的点数占80.1 %，残差绝对值在200m以内的点数占95.7%。

图 8 菲律宾海检核点残差分布比率直方图 Fig. 8 Histogram of the Residual Depth on the Checking Points in Philippine Sea

图选项

从表 1和表 2的统计结果看，本文构建的海底地形模型精度在菲律宾海较南中国海高，这主要有两方面原因：一是菲律宾海北部船测数据分布较南中国海均匀、密集(见图 1)，二是开阔海域卫星测高数据的精度较边缘海高。参与精度对比的几个现有模型中，V15.1的精度最高，并且，它是当前国际上公认精度较高的海底地形模型。V15.1是联合船测海深、历史地形模型和V20.1重力异常数据构建的。与V15.1模型相比，本文模型的主要区别是，在100~200 km波段内采用了垂直重力梯度异常数据进行海底地形的反演计算。

本文综合利用船测海深、垂直重力梯度异常和重力异常数据构建海底地形模型。船测海深数据通常用于构建长波海底地形模型，测高重力异常数据一般用于反演约15~200 km波段内的海底地形，本文用于反演小于100 km波段的海底地形。垂直重力梯度异常数据，首次被用于反演100~200 km波段内海底地形，且取得了较好效果，这有其内在必然性。

首先，垂直重力梯度异常数据可以用于反演20~200 km波段内海底地形。在20~200 km波段内，海底地形与垂直重力梯度异常之间具有很强的相干性，如图 9所示，且Parker式收敛速度很快，式(6)的线性关系成立。图 9采用菲律宾海北部地区船测海深和垂直重力梯度异常数据，数据范围是(20°～30°N,130°～140°E)，用GMT软件直接计算两者间相干性。

图 9 海底地形与垂直重力梯度异常之间的相干性 Fig. 9 Coherency Between Bathymetry and Vertical Gravity Gradient Anomalies

图选项

其次，重力异常之所以仅用于反演约15~200 km波段海底地形，一是因为在长波端，地壳均衡现象、沉积层等可能带来较大误差；二是因为在短波端，受海深向上延拓效应影响，观测数据信噪比较低。文献[18]通过高斯型海山的模拟计算，发现岩石圈有效弹性厚度对垂直梯度异常的影响比对重力异常的影响小。文献[19]也指出，垂直梯度异常数据反演海底地形的优势，主要是它能够抑制沉积层、地壳均衡现象等引起的较长波长误差。理论上，垂直梯度异常较重力异常更能反映短波海底地形起伏，但是垂直梯度异常数据在短波段的信噪比很低，不太实用。因而，本文利用垂直梯度异常数据反演100~200 km波段海底地形模型，而小于100 km的短波海底地形依然采用重力异常反演计算。

在100~200 km波段范围内，海底地形反演时，垂直梯度异常较重力异常的优势，可以通过数值计算实验稍加说明。在菲律宾海北部地区，假设V15.1地形为真值，参考CRUST2.0地壳模型，建立洋壳密度为2 800 kg/m³，地幔密度为3 350 kg/m³，平均水深为-3 400m，洋壳平均厚度为6.5 km的单层洋壳模型。采用三维导纳分析法，联合V15.1模型和测高重力异常，分析得到该地区岩石圈有效弹性厚度约为3 km。基于此模型，分别正演计算100~200 km波段内，海底地形起伏及其挠曲均衡补偿物质对海面重力异常和垂直重力梯度异常的影响。计算结果表明，地形起伏产生的重力异常均方根为102.850 mGal，挠曲均衡补偿物质产生的重力异常均方根为99.733 mGal，与地形影响的均方根之比为96.970％；地形起伏产生的垂直重力梯度异常均方根为12.846 E(1 E=10^-9/s²)，挠曲均衡补偿物质产生的垂直梯度均方根为10.489 E，与地形影响的均方根之比为81.652％。这说明在100~200 km波段内，地壳均衡现象对垂直梯度的影响比对重力异常的影响小。

根据上述洋壳模型和正演计算的重力异常和垂直重力梯度异常数据，分别反演计算海底地形，在100~200 km波段内，重力异常反演结果与V15.1之差的标准差为100.367m，垂直梯度反演结果与V15.1之差的标准差为94.744m，后者精度较前者高5.623％。

综上所述，结合前人研究结论和本文计算结果，可以得出两个结论：

1) 垂直重力梯度异常数据可以被用于反演海底地形模型，且在100~200 km波段内较重力异常反演稍有优势。

2) 在中国海及其周边地区，联合船测海深、垂直重力梯度异常和重力异常数据构建的海底地形模型，其精度较GEBCO、DTU10、ETOPO1和V15.1模型高，在菲律宾海地区，本文计算模型精度较V15.1模型提高了约8.2％。

本文模型以及GEBCO、DTU10、ETOPO1和V15.1模型给出的均是1′×1′的海底地形数据，但模型真实分辨率与采用的数据源有关，并受海深向上延拓效应等物理现象制约。因此，可以认为，式(8)的截断波长即为模型真实分辨率。