基于三棱柱单元的青藏高原重力场地形效应计算

Calculation of Gravity Topographic Effect of Qinghai-Tibet Plateau Based on Triangular Prism Elements

摘要: 地球是一个梨形的不规则椭球体，其表面起伏的地形产生的重力场地形效应对地球物理勘探布格重力异常校正、大地测量学地形校正和确定大地水准面起伏都具有十分重要的意义。青藏高原是全球地质情况最复杂、地形变化最大的区域之一，其产生的地形效应是我国重力勘探普查、反演解释的研究热点。为了计算该地区的重力场地形效应，对该区域进行基础网格建模是最重要的步骤，其次还需进行高精度正演计算。本文提出了一种基于不规则三棱柱单元的带曲率大规模区域网格构建方法，并使用重力场闭合解析公式对其进行了正演计算，我们的结果与前人研究高度一致，证明了建模方法的正确性，同时该方法也可以被用于任意复杂的陆地区域建模，因此也具有应用指导意义。

Abstract: The earth is a pear-shaped irregular ellipsoid, and the topographic effect of the gravity field pro-duced by the undulating topography on its surface is of great significance to the correction of the Bouguer gravity anomaly in geophysical exploration, the topographic correction of geodesy and the determination of geoid fluctuations. The Qinghai-Tibet Plateau is one of the regions with the most complex geological conditions and the largest topographic changes in the world. The topo-graphic effect produced by this area is a research hotspot of gravity exploration census and inversion geological interpretation. In order to calculate the topographic effect of the gravity field the Qinghai-Tibet Plateau, the basic mesh construction of the area is the most important step, followed by high-precision forward calculation. In this paper, a method for constructing a large-scale regional grid with curvature based on irregular triangular prism elements is proposed, and the forward calculation is carried out using the closed analytical formula of the gravity field. The method which is verified correctly can also be used for arbitrarily complex land area modeling, so it also has the guiding significance.

1. 引言

地球重力场作为地下介质密度场源产生的天然物理场，在地球科学中有丰富的研究历史和重要价值。在物理大地测量学领域，通过地球重力场数据研究大地水准面起伏进而确定地球外部的形状。在地球物理学中，重力勘探则是最早在油气勘探中使用的方法，其应用历史可以追溯到19世纪 [1]。地球内部相对于层状球壳理论模型的密度差异是地球外部重力场异常的来源，故重力场数据中包含了地下密度结构的特征响应，如今常与电磁和地震勘探等其它地球物理方法进行联合或交叉梯度约束的反演和解释 [2]，因此确定地球重力场的精细结构及其时间相依变化将为现代地球科学解决人类资源、环境和灾害等紧迫课题提供地学信息。本世纪初，由德国地球科学研究中心和德国航空航天中心研制的CHAMP (Challenging Minisatellite Payload)卫星发射，开启了测量大规模全球重力场的序幕 [3]。后来随着GRACE (Gravity Recovery and Climate Experiment)和GOCE (Gravity field and steady-state Ocean Circulation Explorer)卫星的相继发射，研究人员获得了丰富的全球尺度重力场数据 [4] [5]，并基于球谐系数建立了多种地球重力模型 [6]，大大推动了物理大地测量学、地球物理学和海洋学等学科的发展。在勘探地球物理学领域，重力场数据常常是在空间域中进行处理与反演的，将反演结果结合实际的地质信息来解释地下异常体介质和结构的分布，从而揭示地球内部的演化规律与发展过程 [7]。

重力场地形效应的计算是地球重力场研究的重要应用之一，随着大规模全球尺度、海量的重力场数据的采集，考虑带曲率的大规模区域的重力场地形效应成为了数据预处理的重要手段，也为后续重力场数据反演深部密度物质奠定了基础。此外，通过研究不同地形复杂度下重力场和地形的匹配理论与方法，能解释重力场和地形的变化特征，分析陆海质量迁移过程和洋壳均衡机制及地球圈层物质交换。计算重力场地形效应的方法主要有两种：基于球谐分析的谱域方法和基于球壳单元的空间域方法 [8]。基于球谐分析的谱域方法是根据不同阶次球谐系数建立的重力场模型进行计算，其缺点在于存在一定截断误差，计算时模型的分辨率有限。基于球壳单元的空间域方法近年来发展迅速 [9] [10]，已经成为计算重力场地形效应的主流方法，其核心在于将地形剖分成一系列由经纬度组成的球面六面体单元，称为球壳单元。因为球壳单元自带曲率从而能模拟真实大规模区域和地球模型，但因存在曲线与直线混合从而不能使用闭合解析方式进行重力场正演计算，只能通过数值方法替代求解。此外，近些年来使用多面体剖分方式进行重力场地形效应的研究也陆续出现，基于多面体单元(如四面体、三棱柱单元)计算的好处是可以使用闭合解析公式进行正演，从而得到准确的重力场地形效应。

为了让重力场地形效应结果更加明显，我们通常选取地形起伏大的区域进行计算验证。青藏高原作为全球地形起伏最剧烈的地区之一，包含了丰富的地学信息，也为人们研究欧亚–印度洋板块地质构造、地层变化关系、地震断裂带等提供了重要窗口 [11] [12]。因此计算青藏高原区域的重力场地形效应是十分有必要的，前人也进行许多了研究。球壳单元计算青藏高原地区重力场地形效应的文献较多 [13] [14] [15]；Casenave等基于不规则四面体单元构建的青藏高原区域模型并对地形效应进行了正演计算 [16]；Saraswati等同样基于非结构化四面体网格对青藏高原区域的地形效应进行了闭合解析公式的正演计算 [17]。本文提出了一种新的基于非结构化三角形网格、不规则三棱柱单元构建青藏高原区域模型的方式，该建模流程简便直观，同时考虑了青藏高原这类大规模区域的曲率。为了验证模型的正确性，我们使用Ren等 [18] [19] [20] 提出的闭合解析公式对重力场地形效应进行了正演计算，数值实验表明与前人结果吻合良好，并且我们的建模方法针对各类陆地区域均有指导意义，可以进一步推广应用。

2. 基于三棱柱单元的青藏高原区域建模

2.1. 大规模区域建模概述

大规模模型是相对于局部地球物理测量中的小规模区域而言。在局部地球物理应用中，大地可以看作是一个平面，各类数据处理手段都是在平面直角坐标系下进行，当我们的研究区域扩展到一定程度时，如果不考虑地球曲率的影响则会引入较大误差。相对于构建全球尺度的地球圈层模型，大规模区域模型可以看作是从中提取了一个板块，但从生成方式上来看，其原理与地球圈层模型完全不同。图1给出了一个简单带曲率、基于非结构化网格的大规模区域的轮廓图。

Figure 1. Large-scale regional model based on unstructured grid

图1. 基于非结构化网格的大规模区域模型轮廓图

2.2. 青藏高原二维曲面网格构建

为了构建青藏高原大规模区域的三维三棱柱网格，首先需要创建两个二维的非结构化三角形网格曲面，其中一个平面为高程为0的基准面，另一个则为加入青藏高原地形后的起伏界面。本文采用的地形数据为美国国家海洋和大气管理局下属的国家地球物理数据中心公布的1弧分全球高程数据集ETOPO1 [21]，它基于WGS84坐标系，图2给出了青藏高原建模区域的地形图，最高处为8266 m，最低处为2 m，平均高度2074 m。此次建模区域的经度范围为74˚~96˚，纬度范围为24˚~41˚，因此在实际建模时需要考虑地球曲率的影响，而不是简单当作平面看待。为了创建二维非结构化三角形网格，我们采用开源程序Gmsh进行非结构化剖分 [22]。

Figure 2. Topography in the Qinghai-Tibet Plateau

图2. 青藏高原地形图

Gmsh剖分后的二维表面如图1中的上表面所示，为了更好的模拟地形起伏，在表面处可以采用较小的分辨率以获得精细的结果。而ETOPO1的网格文件是基于地理坐标系的经纬度网格，其所给的高程值只在经纬度网格的节点处有高程值，但我们生成的非结构化网格是在整个表面不规则分布的，显然不会只出现在经纬度网格节点，因此插入地形的整体思路是搜索和插值：在我们现有的表面网格节点上查找离该节点最近的ETOPO1模型的高程点，如果这两个点的重合，则直接将高程赋值给表面节点(但一般不会重合)，对于不重合的点就看表面节点落在ETOPO1数据网格的哪一个栅格中，找到该栅格的东南西北四个节点高程值进行插值得到该点的高程值。这又涉及到插值函数的选择，为了符合带曲率网格的情况，我们选择采用Wang H [23] 等人提出的球面双线性插值法，只要得到经纬度坐标下网格西南、西北、东北和东南四个角的坐标及其高程，就可以插值得到网格内部任意点的高程，该方法已经被证明插值误差几乎可以忽略。图3给出了加入ETOPO1高程地形后的青藏高原非结构化二维三角形网格曲面。

从图3中可以看出，青藏高原的地形还原的较为拟合，同时该二维网格自带曲率，不存在直接使用平面正演计算带来的拟合误差。因为本文提出的建模方法是针对陆地区域的，因此加入地形后的网格高程一般均为正值，这样上表面的网格曲面始终在未加入高程网格的上方，这样便形成了两个曲面的高程差，这也是三维三棱柱单元形成的基础，如图4所示。

Figure 3. 2D surface mesh of Qinghai-Tibet Plateau with topography

图3. 加入地形后的青藏高原区域二维曲面网格

Figure 4. The upper surface with topography (yellow) and the lower surface without topography (green) of the Qinghai-Tibet Plateau region

图4. 加入地形后的青藏高原区域上曲面(黄)和未加入地形的下曲面(绿)

2.3. 青藏高原三维三棱柱网格构建

为了构建三棱柱网格单元，在两个二维非结构化三角形曲面网格的基础上将对应的节点连接即可形成上下表面均为三角形的三棱柱单元。为了保证三棱柱单元的一致性，需要保证上下两个表面三角形节点数量和位置完全一致。这是我们在构建时就已经考虑到的，因为我们首先生成下表面的节点，仅仅通过插值便生成了上表面，此时便只需要连接即可，如图5所示。图5给出了上下表面连接生成三棱柱单元的放大图，可以明显看出上表面为带有起伏地形的曲面，下表面为高程为0的基准面，中间的区域便是真实的地形起伏。图6给出了使用Gmsh可视化后的三维网格切面图，本次构建的网格数量信息为：73,768个节点、146,120个三角形单元、73,060个三棱柱单元。上表面地形起伏与图2给出的吻合良好，为了验证模型的准确型，下一节中我们将使用闭合解析公式进行重力场地形效应计算。

Figure 5. Schematic diagram of triangular prism element

图5. 三棱柱单元生成示意图

Figure 6. The Qinghai-Tibet Plateau constructed by triangular prism elements

图6. 青藏高原三棱柱网格

3. 青藏高原重力场地形效应计算

3.1. 重力场闭合解析正演方法

为了保证获得高精度的重力场地形效应，本次重力场正演计算采用任意多面体闭合解析公式法，每一个三棱柱单元的密度均为常量，取地壳的平均密度2700 kg/m³。下面简单给出正演公式的推导过程，详细推导可参考文献 [18] [19] [20]。对于任意三棱柱单元H，它由五个面 $\partial H$ 组成，即：

$H={\displaystyle \underset{i=1}{\overset{N}{\sum }}\partial H_i}$ (1)

其中 $N=5$ ，有两个面是三角形，三个面是四边形。根据万有引力定律及其重力场理论 [24]，任意三棱柱单元H在观测点产生的重力位 $\varphi \left(r^{\prime }\right)$ 为：

$\varphi \left(r^{\prime }\right)=G{\displaystyle \underset{H}{\iiint }\frac{\lambda \left(r\right)}{R}\text{d}v}$ (2)

其中 $r^{\prime }\left(x^{\prime },y^{\prime },z^{\prime }\right)$ 为观测点坐标， $r\left(x,y,z\right)$ 为产生重力场响应的密度体源点坐标， $\lambda \left(r\right)$ 为密度体在其内部源点所在的密度函数，G为万有引力常量，取为 $6.673\times {\text{10}}^{-11}{\text{m}}^{\text{3}}\cdot {\text{kg}}^{-\text{1}}\cdot {\text{s}}^{-\text{2}}$ ， $R=\left|r^{\prime }-r\right|$ 为观测点与源点之间的距离。根据重力场与重力位的关系，观测点处的重力场 $g\left(r^{\prime }\right)$ 为：

$g\left(r^{\prime }\right)={\nabla }_{r^{\prime }}\varphi \left(r^{\prime }\right)=-G{\displaystyle \underset{H}{\iiint }\lambda \left(r\right){\nabla }_r\left(\frac{\text{1}}{R}\right)\text{d}v}$ (3)

其中 $\nabla$ 为梯度算子， ${\nabla }_{r^{\prime }}$ 代表在观测点处求梯度， $\text{d}v$ 为密度体体积元。上式基于以下矢量恒等式：

${\nabla }_r\left(\frac{\text{1}}{R}\right)=-{\nabla }_{r^{\prime }}\left(\frac{\text{1}}{R}\right)$ (4)

再根据另一个矢量恒等式：

${\nabla }_r\left(\phi \psi \right)=\phi {\nabla }_r\psi +\psi {\nabla }_r\phi$ (5)

其中 $\phi$ 和 $\psi$ 为任意标量函数，带入式(3)中可以得到重力场 $g\left(r^{\prime }\right)$ 展开式：

$g\left(r^{\prime }\right)\text{=}-G{\displaystyle \underset{H}{\iiint }\left[{\nabla }_r\frac{\lambda \left(r\right)}{R}-\frac{\text{1}}{R}{\nabla }_r\lambda \left(r\right)\right]\text{d}v}$ (6)

根据散度定理和公式(1)，将上式改写成：

$g\left(r^{\prime }\right)\text{=}-G\left[{\displaystyle \underset{i=1}{\overset{N}{\sum }}{\hat{n}}_i}{\displaystyle \underset{\partial H_i}{\iint }\frac{\lambda \left(r\right)}{R}\text{d}s}-{\displaystyle \underset{H}{\iiint }\frac{{\nabla }_r\lambda \left(r\right)}{R}\text{d}v}\right]$ (7)

这里 ${\hat{n}}_i$ 为密度体边界 $\text{d}s$ 的外法线方向单位矢量，这样原先的体积分转化成了一个面积分和体积分求和的形式。对于密度函数 $\lambda \left(r\right)$ ，本文使用常量密度函数 ${\lambda }_0\left(r\right)$ ，所以有：

${\nabla }_r\lambda \left(r\right)={\nabla }_r{\lambda }_0\left(r\right)=0$ (8)

代入式(6)中右端第二项变为0，故常量密度函数的三棱柱单元的重力解析解表达式：

$g_0\left(r^{\prime }\right)\text{=}-G{\lambda }_0\left(r\right){\displaystyle \underset{i=1}{\overset{N}{\sum }}{\hat{n}}_i}{\displaystyle \underset{\partial H_i}{\iint }\frac{1}{R}\text{d}s}$ (9)

令 $N=5$ ， ${\lambda }_0\left(r\right)=2700$ 并且给出三棱柱的五个面的信息以及六个点的坐标，式(9)便可以直接算出，任意观测点处的重力便可以得到。由于重力场是位场，满足叠加原理，所以计算青藏高原处测点的重力场总值时只需将组成青藏高原地形的所有三棱柱单元的重力场值叠加即可。

3.2. 模型验证

本次测试模型的测点安排青藏高原区域上表面所有地形节点高1 cm处，可以认为是实际重力勘探在地表处测量得到的结果。为了更加直观的表现重力场分量特征，我们通过坐标转换将在直角坐标系下的结果转换为了球坐标系下。球坐标系的球心放置在地球球心，r为径向方向，正方向指向外； $\theta$ 为纬度方向，范围为−90˚~90˚； $\phi$ 为经度方向，范围为−180˚~180˚。图7~9分别给出了青藏高原区域重力场地形效应分量 $g_r$ 、 $g_{\theta }$ 和 $g_{\phi }$ 等值线图。

Figure 7. $g_r$ component of gravity topographic effect of the Qinghai-Tibet Plateau

图7. 青藏高原重力场地形效应 $g_r$ 分量等值线图

Figure 8. $g_{\theta }$ component of gravity topographic effect of the Qinghai-Tibet Plateau

图8. 青藏高原重力场地形效应 $g_{\theta }$ 分量等值线图

Figure 9. $g_{\phi }$ component of gravity topographic effect of the Qinghai-Tibet Plateau

图9. 青藏高原地形效应 $g_{\phi }$ 分量等值线图

从图7可以看出， $g_r$ 分量对地形的反映较为清晰，在地形高的地方(如喜马拉雅山脉)出现了异常极大值，在平坦的区域(如印度东北部)地形效应接近于0。出现异常负值是因为r方向的正方形指向外，而 $g_r$ 分量是指向内的。该结果与前人使用四面体单元、球壳单元正演的结果相比具有良好的一致性，证明了我们建模的准确性。图8和图9分别给出了 $g_{\theta }$ 纬度和 $g_{\phi }$ 经度分量的等值线图，可以看出纬度分量对上下边界勾勒明显，而经度方向对东西方向边界勾勒明显，大部分区域的两者地形效应也基本为0。

3.3. 对比分析

传统的球壳单元模拟地形起伏时上下表面为四边形网格，而三棱柱单元上下表面由非结构化三角形组成，相比于球壳单元能够更好的模拟地形起伏情况，更接近真实地形；在正演计算时，球壳单元正演地形效应只能采用基于高斯–勒让德积分或者泰勒展开等数值方法，虽然精度在诸多文献改进后有所提高但仍比闭合解析公式精度低。故使用三棱柱单元计算的重力场地形效应更准确，但当测点很多时解析方法耗时较长，同时非结构化三角形网格的构建过程很依赖于剖分程序，不如球壳单元生成那样更利于控制。

相比于使用其它网格单元进行的闭合解析公式正演，三棱柱单元的优势在于生成方式简单，例如相比于Casenave等 [16] 使用的四面体区域构建方法就容易许多。但劣势在于目前三棱柱单元建模方式只适用于陆地区域，即只能处理高程均大于(或小于) 0的区域，对于带有海岸线的海陆交叉区域还无法构建。

4. 结论

通过构建基于非结构化三角形网格的两个曲面构建了青藏高原的三棱柱单元地形模型，再通过重力场闭合解析公式对青藏高原重力场地形效应进行了正演，得出了以下结论：

1) 本文提出的针对任意陆地区域的三棱柱网格建模得到了验证，并以青藏高原大规模区域为例得到了精确的结果。

2) 青藏高原重力场地形效应显示，对于球坐标系下的重力场分量 $g_r$ 对地形反映最为清晰，而 $g_{\theta }$ 和 $g_{\phi }$ 分量对边界轮廓勾勒的明显。如果要对该地区的地形效应进行更进一步的研究，计算球坐标系下的全张量重力场梯度分量是将来的研究重点。

3) 三棱柱单元生成方式简单易懂，计算的重力场地形效应相比基于球壳单元的空间域正演方法更准确，劣势在于当测点很多时解析方法耗时较长，同时三棱柱单元建模方式只适用于陆地区域，生成时很依赖二维的非结构化三角形网格，不利于后续人为控制。