基于重力异常迭代延拓的南海海底地形反演

郭金运，魏志杰，祝程程，吴云龙，纪 兵

(1.山东科技大学 测绘与空间信息学院，山东 青岛 266590； 2.中国地震局地震研究所 中国地震局地震大地测量重点实验室，湖北 武汉 430071； 3.海军工程大学 导航工程系，湖北 武汉 430033)

在海洋科学中，海底地形数据分析可为海洋地质学、海洋地球物理学、海洋生态保护等提供可靠的研究资料。自20世纪70年代以来，卫星测高技术经过50多年的发展，其技术性能愈发成熟，测高精度和分辨率有了很大提升，海面高精度不断提高[1]，海洋重力场研究不断发展。1978年，Rummel等[2]利用逆Stokes公式通过大地水准面高确定了空间点的重力异常。最小二乘配置法(least squares collocation, LSC)作为大地测量学的经典方法[3]，得到不断改进和完善，现已成为利用卫星测高数据恢复重力场模型的重要方法。此外，逆Vening-Meinesz公式法[4]也是重力场反演的一种算法。随着发射卫星的增多以及卫星数据质量的提高，联合多卫星数据资料提高重力场反演精度[5-7]逐渐成为一种普遍的方法。

在重力场精度逐渐提高过程中，重力场反演海底地形的理论方法不断发展。目前常见的反演海底地形的方法有重力地质法[8]、最小二乘配置法[9]和导纳函数法[10-12]等，导纳函数法相对于前两种方法的优点是利用重力场和海底地形的导纳函数关系，可以很好地应用于有观测重力数据的海域，体现了卫星测高数据反演重力场的优势，大大提升了反演海底地形的效率。导纳函数不仅包含地形起伏信息，还包括岩石圈的均衡状态信息，为了提高海底地形反演的精度，需考虑岩石圈的均衡补偿效应[13]。地形反演过程对重力异常数据质量要求很高，一般需要将海面上的重力异常延拓到平均水深面。快速傅里叶变换(fast Fourier transform, FFT)是一种常用方法，向上延拓是稳定的过程，但向下延拓是误差放大的过程，结果不稳定。因迭代延拓法具有延拓深度范围大、结果稳定性强的特点，迭代法向下延拓开始得到应用，但对重力异常的迭代延拓法在地形反演研究的应用很少提及。

本研究以南海海域中沙群岛南部海域(114°E～118°E，12°N～15°N)为试验区，利用最小二乘配置法得到HY-2A重力异常，通过迭代延拓法得到延拓面上的重力异常，确定重力异常和海底地形之间相干性大于0.5的波长范围为20～115 km。采用CRUST 1.0全球地壳模型计算得到平均地壳厚度、地幔密度和有效弹性厚度等参数，进而求得试验区域均衡补偿模型。利用重力异常数据反演得到未加补偿的导纳水深模型和均衡补偿导纳水深模型，并将两个水深模型与ETOPO1水深模型和DTU10水深模型进行精度对比分析，讨论引起误差的原因。

1 原理与方法

1.1 迭代延拓法原理

重力异常向上延拓是收敛的过程，结果具有唯一性，而向下延拓是发散过程，结果稳定性差，误差较大。图1为海平面和延拓面的位置。

图1 海平面和延拓面的位置

假设海平面为S1=h，延拓面为S0=0，延拓高度为h=S1-S0。在S0面上重力异常的信息为p0(x,y,0),通过FFT得到P0(kx,ky,0)=F[P0(x,y,0)]，kx、ky是各自方向的波数。然后通过逆FFT延拓公式得到对应S1面的重力异常信息，公式为：

(1)

其中F-1为逆傅里叶变换。当h>0时，向上延拓；当h<0时，向下延拓。采用迭代法通过向上延拓计算得到向下延拓面的重力异常[14]，步骤如下：

1) 将海平面S1上已知点P的重力异常Pg0作为对应延拓面S0上重力异常的初始值；

2) 通过公式(1)计算海平面上的重力异常值Pgi，根据Pg0与Pgi的差值对Pg0进行改正；

3) 反复计算迭代，当计算值与初始值在误差允许范围时，即可得到延拓面S0上的重力异常。

1.2 导纳函数理论

“导纳”描述了海洋重力异常和海底地形在频率域上的频谱关系。Parker[15]在1973年将FFT引入位场理论，提高了密度界面反演的计算速度。在频率域内，重力异常和海底地形之间的关系为[16]:

(2)

一般地，岩石圈对海底地形载荷的均衡响应是通过区域补偿机制来实现的，岩石圈强度决定了挠曲量的大小。当海底地形载荷的波长小于岩石圈弹性厚度时，公式(2)等式右边第一项起主要作用，则公式(2)简化为：

G(k)=Z(k)H(k)，

(3)

其中：G(k)为重力异常的傅里叶变换；H(k)为海深的傅里叶变换；导纳函数为

Z(k)=2πG(ρc-ρw)exp(-2πkd)。

(4)

公式(4)为未考虑未补偿模型的导纳函数式。

1.3 挠曲均衡补偿导纳理论

实际问题中，Z(k)不仅包含海底地形信息，还包含岩石圈的均衡状态信息，只有确定岩石圈挠曲产生的重力异常，才能确定Z(k)。由于无法直接观测到岩石圈发生的挠曲情况，不能直接计算岩石圈的挠曲所产生的重力异常。为研究方便，可通过挠曲均衡模型来研究岩石圈的内部信息。当海底地形载荷的波长与岩石圈挠曲波长相近或者更长时，岩石圈在海底地形载荷作用下发生挠曲。根据挠曲均衡模型[13]，均衡补偿物质产生的重力异常与海底地形之间的关系为：

G(k)=2πG(ρc-ρw)Φ′e(k)exp(-2πk(d+t))H(k)，

(5)

(6)

其中：v为岩石圈的泊松比；E为弹性模量。结合公式(3)、(5)，整理得到在挠曲均衡模型下地形和补偿物质产生的重力异常为：

G(k)flex=2πG(ρc-ρw)e-2πkd(1-Φ′e(k)e-2πkt)H(k)=Z(k)flexH(k) ，

(7)

其中，重力异常的均衡补偿导纳函数

(8)

由公式(8)知，挠曲均衡补偿导纳函数与地壳平均厚度t、海底洋壳密度ρc、岩石圈挠曲刚度D以及弹性厚度Te等有关。

1.4 相干性分析

(9)

其中：*表示复数共轭；〈〉表示括号内功率谱密度乘积的实部在圆环区域内(ki-Δk≤k≤ki+Δk)的平均值。

2 试验数据及结果分析

2.1 研究区域

南海盆地的中央区域为深海平原，平均水深约为4 100 m，将114°E～118°E、12°N～15°N作为试验区域。在频率域上为消除计算过程中边缘效应的影响，实际计算时经纬度分别向外扩展2°，然后对计算结果进行处理，得到研究区域的海底地形。

2.2 数据来源

2.2.1 重力异常数据

各镇、各单位要利用多种形式，积极宣传农田林网建设与改善生态环境、提高人民群众生活质量的关系，与抵御自然灾害、改善农业生产条件的内在联系。积极宣传绿化先进典型，曝光毁林案件，通过一系列宣传教育，不断提高广大干部群众的绿化、美化意识，形成“植绿、爱绿、护绿、兴绿”和“保护生态、绿化家园”的文明行为，进一步增强农民群众建设农田林网的热情，形成全民动员、全民参与、全民动手，推进农田林网建设的浓烈氛围。

根据2011年8月发布的HY-2A卫星大地测量任务(geomatic mission, GM)的测高数据，获得南海海域(105°E～122°E, 2°N～26°N)1′×1′网格的重力异常。2016年3月30日进行变轨后开始GM，GM的周期从14天增加到168 天，工作范围为81°S～81°N。选取时间2016年3月30日—2018年8月22日，将采样频率1 Hz的HY-2A/GM 卫星高度计数据Level 2 Plus (L2P)产品[17]作为研究数据。首先，对HY-2A/GM 的海面高度进行校正和粗差消除预处理；然后，利用预处理的海面高来计算沿轨的大地水准面梯度，通过移去EGM2008的大地水准面梯度获得残余大地水准面梯度；最后，通过计算窗口半径为0.5°的最小二乘配置法，从残余大地水准面梯度求得1′×1′网格上的残余重力异常，通过恢复EGM2008的重力异常，从残余重力异常中获得最终的重力异常模型[18]，试验区的重力异常如图2(a)所示。

2.2.2 ETOPO1水深模型和DTU10水深模型

ETOPO1水深模型是由美国国家海洋与大气管理局(National Oceanic and Atmospheric Administration, NOAA)的国家环境信息中心(National Centers for Environmental Information, NCEI)提供的分辨率为1′的全球地形数据模型[19]，集成了来自多个国际组织的船测水深数据、测高重力数据、海岸线数据和陆地高程数据的评估和汇总。ETOPO1模型是海洋地球物理研究中使用最广泛的数据之一。

DTU10水深模型来自丹麦科技大学(Technical University of Denmark, DTU)发布的分辨率为1′的全球地形数据模型[20]，通过ERS-1和GEOSAT大地测量任务的高度计数据计算得到全球地形模型数据。

2.2.3 船测水深数据

船载水深数据从NCEI官网下载得到[21]，时间跨度为1960—2016年。NCEI控制数据质量并将合格的数据收集到数据库中。尽管总体数据质量准确可靠，但一些早期数据存在较大的测量误差，因此有必要剔除这些误差数据。将船测水深数据减去ETOPO1数据，差值大于船测水深数据的点将被删除。图2(b)为试验区的船测轨迹，黑色点表示试验控制点(13 934控制点)，白色点表示水深反演结果的检核点(9 464个)，背景为ETOPO1海深模型。

图2 重力异常图和船测水深轨迹分布图

2.2.4 CRUST 1.0全球地壳模型

CRUST 1.0为全球三维地壳模型，其中的海水深度和地形数据是从ETOPO1地形模型中提取的，并通过平均处理得到间隔1°的网格数据。官方下载网址为：http://igppweb.ucsd.edu/～gabi/crust1.html[22]。

2.3 迭代延拓法延拓结果

将HY-2A重力异常数据向下延拓到平均水深面得到迭代次数和误差的关系，如图3所示。

图3 迭代次数和延拓误差的关系

2.4 确定海底地形反演波段范围

根据公式(9)对重力数据和海底地形在频率域进行相干性分析。信噪比(signal noise ratio, SNR)描述信号与噪声的比例关系：

(10)

其中：当相干coh=0.5时，SNR为1；当coh>0.5时，信噪比大于1。

试验区内重力异常数据和海底地形的相干性如图4所示。由图4看出，在波长范围20～115 km内，重力异常与海底地形的相干性>0.5，选取该波段进行海底地形反演。

图4 重力异常和海深的相干性

2.5 挠曲均衡补偿模型参数的确定

在试验区域内通过均衡补偿模型计算得到地球物理参数如表1所示。

表1 挠曲均衡补偿模型所需要的参数

通过CRUST 1.0全球地壳模型获取海水密度ρw，海底洋壳密度ρc，地幔密度ρm，平均地壳厚度t，岩石圈泊松比v和弹性模量E。平均海水深度由NCEI发布的船测水深数据获得。

在实际计算时，通过分析观测导纳和理论导纳曲线的走势来确定岩石圈的有效弹性厚度Te，并通过重力异常和海深分析二者的“观测导纳”关系，公式如下：

(11)

当Te分别取2、6、10 km时，导纳函数Z(k)flex的曲线变化如图5所示。低频部分曲线主要受有效弹性厚度Te的影响，反映了深层地壳物质对重力场的影响；当Te取6 km时，低频部分“观测导纳”和“理论导纳”曲线走势最接近，因此确定研究区域的弹性厚度为6 km。

图5 不同有效弹性厚度下“理论导纳”和“观测导纳”的曲线

2.6 海底地形反演

采用移去恢复法构建海底地形，在未补偿条件和均衡补偿条件下，通过导纳函数法反演海底地形的步骤如下：

1) 将试验区的船测水深数据进行网格化，得到网格化水深模型；

2) 将网格化的水深进行二维傅里叶变换，通过滤波处理得到短波段(<20 km)海深模型和长波段(>115 km)海深模型；

3) 重力异常数据经带通滤波处理后，在未补偿条件和均衡补偿条件下，分别利用迭代延拓法向下延拓得到延拓面的重力异常，在波段20～115 km波长范围内反演海深；

4) 将各个波段的海深相加，得到未补偿条件下的海深模型和均衡补偿条件下的海深模型。将未补偿导纳水深模型记为Model1，均衡补偿导纳水深模型记为Model2。

2.7 试验结果

图6显示了试验区域反演的两个海深模型Model1、Model2与ETOPO1水深模型和DTU10水深模型的结果，可以看出试验区域整体的地形分布，4个模型都描绘了试验区域的主要地形。表2为格网化的NCEI船测水深数据、Model1、Model2、ETOPO1水深模型和DTU10水深模型的数据统计。NCEI船测水深数据、Model1、Model2和ETOPO1水深模型的平均水深相差不大，Model2相对于Model1精度有所提升，而DTU10水深模型的海深平均值相对于其他数据海深平均值偏差较大，误差超过100 m。

(a)Model1；(b)Model2；(c)ETOPO1水深模型；(d)DTU10水深模型图6 海底地形模型

表2 NCEI、Model1、Model2、ETOPO1和DTU10水深数据统计表Tab. 2 Data statistics of NCEI, Model1, Model2, ETOPO1 and DTU10 bathymetry m

2.8 精度分析

为进一步分析各模型精度，将4个海深模型在NCEI船测水深检核点上进行模型插值计算得到相应的海深值，通过与NCEI船测水深数据作差比较各自的模型精度。表3给出了各海深模型在检核点的差异统计结果。本研究反演的两个水深模型差异的平均值分别为2.5、-1.4 m，ETOPO1水深模型与反演的两个水深模型平均海深差异接近，但DTU10水深模型差异的平均值则相对较大。Model1、Model2和ETOPO1水深模型的相对误差接近，而DTU10水深模型相对误差达4.48%。按照标准差质量统计来看，水深模型精度从高到低依次为Model2、Model1、ETOPO1水深模型和DTU10水深模型(标准差分别为150.8、 159.7、164.3、196.1 m)。由表3可知，通过迭代延拓法处理得到的重力异常，经反演计算得到的两个水深模型精度均优于ETOPO1水深模型和DTU10水深模型。

表3 Model1、Model2、ETOPO1和DTU10数据在检核点上误差的统计结果Tab. 3 Data statistical results of error of Model1, Model2, ETOPO1 and DTU10 at check points

图7分别显示了Model1、Model2水深差分布直方图，海深差在50 m范围内分别占51.7%和62.3%，海深差在200 m范围内分别占90.6%和91.9%。图8分别显示了两个反演海深模型与NCEI船测水深数据的相对误差和海深差分布关系，结果没有明显的差异特征，但在大于-3 000 m的深海区域反演海深精度相对较高。图9显示了Model1、Model2与NCEI船测水深数据相对误差的空间位置分布。

图7 海深差统计直方图

图8 相对误差和海深差随海深变化的情况

图9 相对误差的空间位置分布情况

通过上述统计结果，结合图6和图9海底地形模型分析可知，在中南暗沙附近(114.75°E～115.75°E，13.25°N～14.5°N)海底地形起伏较大，反演的水深精度相对较差，最大海深差可达2 000 m以上，ETOPO1水深模型和DTU10水深模型也存在相同的情况，这是地形变化剧烈影响的结果，是目前海底地形模型存在的问题。Model1在整体精度上优于Model2和ETOPO1模型，但Model1和Model2水深模型相对于ETOPO1水深模型在水深刻画上存在不足。

通过上述分析，认为影响水深反演的可能原因为：①船测水深数据时间跨度大，由于时效性问题前期测量的数据精度有待提高，粗差很难完全剔除，可能导致在分析海底地形与重力异常的导纳性上出现误差；②理论上海底地形与重力数据在频率域上存在良好的线性关系，但实际上海底地形复杂，并且在频率域上求二者相干性时用到的圆环内的平均值也会产生误差，导致反演水深时产生误差；③将船测水深数据和离散数据网格化导致格网数据精度不高，得到的水深模型也会受到影响；数据在滤波过程中有信号损失，是导致精度损失的原因之一。

3 结论

通过重力导纳理论，利用HY-2A重力异常延拓数据反演得到未补偿导纳水深模型Model1和均衡补偿导纳水深模型Model2，并将反演结果与目前公认的ETOPO1水深模型和DTU10水深模型进行外部检核比较，得到如下结论：

1) 利用迭代延拓法将HY-2A重力数据向下延拓，得到平均水深面上的重力异常，其精度误差在0.1 mGal以内。

2) 通过实际计算得到的“观测导纳”和均衡补偿理论计算的“理论导纳”进行比较，可以看出低频部分主要受地壳有效弹性厚度的影响，最后确定试验区域内的有效弹性厚度为6 km。

3) 虽然整体上Model2精度略高于Model1，标准差提升约9 m，但精度提升效果不明显，说明有效弹性厚度为6 km时，均衡补偿模型对试验区地形的精度改善效果有限。

4) 利用迭代延拓法求得延拓面上的重力异常，经反演计算得到的两个水深模型精度均优于ETOPO1水深模型和DTU10水深模型，并且Model2优于Model1。通过误差分析，在地形变化剧烈的区域Model1、Model2、ETOPO1水深模型和DTU10水深模型相对于船测水深数据精度较差，这也是目前水深模型存在的问题。