利用垂直重力梯度异常反演海底地形的解析方法
徐焕1,2
, 于锦海1,2
, 安邦1,2,万晓云3
1. 中国科学院计算地球动力学重点实验室, 北京 100049;
2. 中国科学院大学地球与行星科学学院, 北京 100049;
3. 中国地质大学(北京)土地科学与技术学院, 北京 100083
基金项目：国家重点研发计划（2016YFB0501702）；国家自然科学基金（41774089；41674026）
关键词：海底地形 垂直重力梯度 边界效应 远区影响 正则化
引文格式：徐焕, 于锦海, 安邦, 等. 利用垂直重力梯度异常反演海底地形的解析方法[J]. 测绘学 ，2022，51(1)：53-62. DOI:
10.11947/j.AGCS.2022.20200578
XU Huan, YU Jinhai, AN Bang, et al. An analytical method for bathymetry inversion using vertical gravity gradient anomaly[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(1): 53-62. DOI: 10.11947/j.AGCS.2022.20200578
阅读全文：
http://xb.sinomaps.com/article/2022/1001-1595/2022-1-53.htm
引 言
海底地形和海底基岩密度可以反映出海底构造演化活动[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%的结论。
1 理论与方法
1.1 长方体产生的垂直重力梯度
假设有一个海底长方体海山Ω(图 1)，其长、宽、高分别为2a、2b、H–h，其中H为Ω底部至海平面的高度，h为Ω顶部至海平面的高度(即Ω的海深)。记R是Ω在海面对应的海域，显然确定Ω的形状仅需求解h。假设P是在海平面上任意一点，下面计算Ω在P处产生的VGG。
图 1 长方体海山Fig. 1 Schematic diagram of synthetic cylindrical seamount
图选项
建立图 1所示的坐标系，则Ω在P点产生的重力位为
(1)
式中，G是万有引力常数(6.672×10-11Nm2kg-2)；ρ=2700kg/m3是海山Ω的密度；(x1,y1,z1)为积分变量。因此，海山Ω在P点产生的VGG为
(2)
式中
(3)
积分区间
为矩形区间，进一步计算后可得IR(H,x,y)的表达式为
(4)
对于一般的矩形积分区间
;
，结合式(4)，IR(H,x,y)可写为
(5)
同理，式(2)中的IR(h,x,y)与IR(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足够小，在Rij中h(x,y)可以被看成是常数hij，于是，利用式(4)、式(5)和叠加原理可知，由Ω生成的VGG在海面P(x,y)处的值Γ(x,y)为
(6)
在式(6)中去掉海水的影响，则得Ω在海面上P处产生的VGG异常为
(7)
式中，Δρ=ρ–ρ0，而ρ0表示海水密度，在后面的计算中总是取ρ0=1030kg/m3。
如果事先能够给出海域R上的VGG异常的观测值，例如在R的所有格 点上VGG异常是已知的，可得
(8)
式中，(xp,yq)是海面上区域R的观测点，p,q=-N, …, 0, …N。这样VGG异常仅由Ω产生的前提下，建立了关于hij的观测方程组(8)，其中hij就是子格 Rij的平均海深。在式(8)中，方程的个数为(2N+1)2，未知数hij的个数为(2N)2，故满足可解条件。事实上，使用hij来替代h(x,y)的本质就是分块求出Ω的海深，即使用分片函数来表达Ω的海深。因此，当分割步长越小时，预测的Ω的地形就会越准确。由于式(8)关于hij是非线性的，因此，在给定初始值hij(0)后，可用迭代的方式来求解，即
(9)
式(9)就是线性化后的海底深度与VGG异常之间的观测方程组。
图 2 模拟20km×20km海底地形及对应格 Fig. 2 The simulated 20km×20km seafloor topography and corresponding grid diagram
图选项
1.3 观测方程组抗误差干扰分析
图 3 抗随机误差和系统误差的曲线图Fig. 3 The graph of the RMS error of the random error and the system error
图选项
2 边界效应的处理
上文建立的关于海深的观测方程式(8)或式(9)都是在海底仅有一个异常海山的前提下导出的，当处理实际海域时，显然海域下方的情况要复杂得多，特别是还存在着其他的异常海山。例如，图 4(a)模拟了40km×40km海域D下方的海山，如果仅需要求解中心20km×20km区域(称为内区R)下方的海底地形，显然在内区外的海山会对内区R上的VGG异常δΓ有贡献，因此，式(8)或式(9)需要进行调整。本节主要以图 4(a)模拟的海山为例来介绍求解方法和进行结果分析。
图 4 模拟40km×40km海底地形及对应格 化Fig. 4 The simulated 40km×40km seafloor topography and corresponding grid diagram
图选项
事实上，在式(7)中将区域取为D，而将P(x,y)点仅限于R中，则有
(10)
这里当i,j=-10, …, 9时Dij=Rij。因为式(10)中待求解的量hij共有402=1600个，所以在R中选择观测点
如下：
，
其中p,q=-20, …, 20。于是，可得观测方程组如下
(11)
式中，p,q=-20, …, 20。因此方程组(11)共含有412个方程。仿照式(9)，在(11)中对hij进行线性化，可得形如Ah=b的线性方程组，其中A是412×402矩阵，h是所有待解量hij排成的402阶向量，而b则是412阶常数向量。由于在(11)中观测点
均在R中，所以线性化后的方程组Ah=b可能会产生奇异性。为了保证可解性，引入正则化方程如下
(12)
这里E是单位矩阵，α>0是正则化参数(其单位为10-24m-2s-4)。
图 5 不同正则化参数α对应的反演结果的RMS误差Fig. 5 The RMS error of the results corresponding to different regularization coefficients α
图选项
由图 5可知，当正则化参数α=10-3时海深hij的RMS误差达到最小，此时在内区R上恢复的海深的误差分布由图 6给出。由图 6可知，最大误差为0.64m，RMS误差为0.48m。因此通过将内区R扩充到D，并引入正则化参数，就可以有效地处理D–R下方海山引起的边界效应。
图 6 考虑边界效应后，反演的内区下方海底地形的误差分布Fig. 6 The error distribution of the inverted seabed beneath Rafter introducing boundary effect
图选项
3 远区海山的影响
在上节扩展区域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异常趋于零。
图 7 沿y=0，从D的左边起算的Ω生成的VGG异常的图像Fig. 7 Alongy=0, the graph of VGG anomalies generated byΩfrom the left ofD
图选项
图 8 考虑远区影响Ω后的反演结果误差分布Fig. 8 The error distribution of the inverted seabed topography with the disturbance seamountΩ
图选项
由图 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)是不存在边界效应的。由此也可以得出这样的结论：有了边界效应后，反解内区下方海底地形的抗误差干扰性会降低。
(13)
然后在式(13)中同时对hij和ε进行求解。至于如何求解式(13)，仍然采用对hij进行线性化，再使用正则化参数的方法来迭代解算hij和ε。
按照上面的方法，对Ω的远区影响问题进行了计算，反演得到了内区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)作为研究对象，密度差Δρ=ρ–ρ0选取为1670kg/m3。此外，根据已有的海深数据，选择H=4.4km。
首先在72km×72km的海域D中选择内部52km×52km海域R作为内区，其上的VGG异常数据(图 9(a))来自GFZ。然后选择步长t=2km，由此便可以建立方程式(13)，对hij进行线性化，引入正则化参数α=10-3，迭代后便能解出R的子分割Rij的平均海深hij。最后将海深hij置于Rij的中心，便可生成R下方的海底地形如图 10所示。
声明：本站部分文章内容及图片转载于互联 、内容不代表本站观点，如有内容涉及侵权，请您立即联系本站处理，非常感谢！