三维Tesseroid网格模型重力异常正演方法及并行算法

图1 Tesseroid单元体示意^[8]

Fig.1 Tesseroid model diagram^[8]

其中:

$\begin{array}{l} l = \sqrt[]{r_{o}^{2} + r_{s}^{2} - 2 r_{s} r_{o} \cos ψ} \\ \cos ψ = \sin φ_{o} \sin φ_{s} + \cos φ_{o} \cos φ_{s} \cos (φ_{s} - φ_{o}) \end{array}$

式中:G表示万有引力常量,G=6.67×10^-11 N·m²/kg²,ρ表示模型体密度,单位为g/cm³,观测点坐标Q₀(r₀,φ₀,λ₀),Tesseroid单元体中心点坐标Q(r_s,φ_s,λ_s),l为观测点Q₀与单元体中心点Q之间的距离,ψ为Q₀点与Q点相对于球心的夹角。

1.2 改进的高斯—勒让德积分法

将模型积分核函数积分区间化为勒让德多项式,随后使用近似积分核函数的勒让德多项式计算累加得到模型异常值。对球坐标系下重力异常积分表达式(1)进行GLQ积分分解,得到Tesseroid单元体的勒让德多项式表达式如下^[5]:

(2)

\begin{array}{l} g (γ, φ, λ) = \frac{(r_{2} - r_{1}) (φ_{2} - φ_{1}) (λ_{2} - λ_{1})}{8} \times \\ \overset{nλ}{\sum_{i = 1}} \overset{nφ}{\sum_{j = 1}} \overset{nr}{\sum_{k = 1}} ω_{i} ω_{j} ω_{k} \frac{Gρr'^{2} \cos φ' (r - r' \cos ψ)}{l^{3}}, \end{array}

式中:

$\begin{array}{l} λ_{si} = \frac{λ_{si}^{'} (λ_{2} - λ_{1}) + (λ_{2} + λ_{1})}{2}, \\ ω_{i} = \frac{2}{n_{λ} P_{nλ - 1} (λ_{si}^{'}) P_{nλ}^{'} (λ_{si}^{'})}, \\ φ_{sj} = \frac{φ_{sj}^{'} (φ_{2} - φ_{1}) + (φ_{2} + φ_{1})}{2}, \\ ω_{j} = \frac{2}{n_{φ} P_{nφ - 1} (φ_{sj}^{'}) P_{nφ}^{'} (φ_{sj}^{'})}, \\ r_{sk} = \frac{r_{sk}^{'} (r_{2} - r_{1}) + (r_{2} + r_{1})}{2}, \\ ω_{k} = \frac{2}{n_{r} P_{nr - 1} (r_{sk}^{'}) P_{nr}^{'} (r_{sk}^{'})} 。 \end{array}$

λ_si,φ_sj,r_sk是球坐标系下Tesseroid单元体经度、纬度以及半径方向坐标积分区间的高斯节点, $λ_{si}^{'}$ , $φ_{sj}^{'}$ , $r_{sk}^{'}$ 为单元体各坐标方向(-1,1)区间内的高斯节点。P_nr_-1( $r_{sk}^{'}$ )是单元体高斯节点 $r_{sk}^{'}$ 处(nk-1)阶高斯—勒让德多项式的值。 $P_{nr}^{'}$ ( $r_{sk}^{'}$ )是单元体高斯节点 $r_{sk}^{'}$ 处(nk)阶高斯—勒让德多项式的一阶导数值。根据推导出来的表达式,可以发现上述GLQ积分法可以将模型体重力异常的计算从表面转移到模型体内部,有效地避免奇异值的产生,美中不足的是,该方法同时也带来了较大的计算量。

Heck和Seitz^[1]以及王祥等^[8]针对这一问题对GLQ积分法进行了改进:

(3)

\begin{array}{r} g (γ, φ, λ) & = Gρ \int_{λ_{1}}^{λ_{2}} \int_{φ_{1}}^{φ_{2}} \int_{r_{1}}^{r_{2}} \frac{r_{s} λ_{s}^{2} \cos φ_{s} (r_{o} - r_{s} \cos ψ)}{l^{3}} d r_{s} d φ_{s} d λ_{s} \\ = Gρ \int_{φ'}^{=} \int_{λ'}^{=} f (r, λ, φ, λ', φ', r_{1}, r_{2}) d φ'd λ' \\ = \frac{(λ_{2} - λ_{1}) (φ_{2} - φ_{1})}{4} \\ \overset{J}{\sum_{nj = 1}} \overset{I}{\sum_{ni = 1}} ω_{ni} ω_{nj} f (r, λ, φ, {\dot{λ}}_{ni}^{'}, {\dot{φ}}_{nj}^{'}, r_{1}, r_{2}) \\ = \frac{(λ_{2} - λ_{1}) (φ_{2} - φ_{1})}{4} \\ \overset{J}{\sum_{nj = 1}} \overset{I}{\sum_{ni = 1}} ω_{ni} ω_{nj} \frac{Gρ \cos φ'}{r} \times \frac{r'^{3}}{l} - l (r' + 3 r \cos ψ) - \\ r^{2} (3 co s^{2} ψ - 1) \ln (l + r' - \cos ψ) 。 \end{array}

与传统的GLQ积分法对比,Heck和Seitz^[1]改进后的GLQ积分法先对球坐标系下重力异常积分表达式r方向积分,王祥等^[8]进一步作GLQ积分分解,实现了用较少的权重控制点来计算模型重力异常(即改进后不在r方向展开),既提升了模型重力计算的精度,又减少了模型计算的时间。这种改进后的GLQ积分法称之为改进的高斯—勒让德积分,即GLQ2D算法。

为了进一步提高地表观测面的正演精度,Hao等^[9]在王祥等^[8]的基础上对前人的自适应剖分方案进行简化改进,引入根据不同计算精度要求而给定的距离判定值W。对一个径向、纬向和经向的计算范围分别为[r₁,r₂]、[φ₁,φ₂]、[λ₁,λ₂]的Tesseroid单元体是否进行进一步的剖分,可以由判别式(4)来进行判定。

(4)

\frac{d}{L} \geq W,

其中:

$\begin{array}{l} d = \sqrt[]{r^{2} + r_{c}^{2} - 2 r r_{c} \cos ψ_{c}} \\ L = \max (L_{φ}, L_{λ}) \end{array}$

式中:L_φ,L_λ根据Tesseroid单元体的弧长与周长的比例得来。根据改进后的自适应剖分方案,有效的解决了近地表观测时邻近效应带来的精度不足的问题。

1.3 三维Tesseroid网格模型重力正演流程及并行算法

本文球坐标系下三维网格模型正演流程与直角坐标系下三维网格模型正演的类似。图2显示了球坐标系三维网格模型重力正演串行算法流程图。

图2

图2 三维Tesseroid网格模型重力正演串行流程

Fig.2 Serial flow chart of gravity forward modeling based on a 3-D Tesseroid mesh model

本文正演串行算法步骤如下:① 输入三维Tesseroid网格模型相关数据以及观测面相关数据,模型数据包括模型经度、纬度、深度方向的坐标最小值、纵横向步长、Tesseroid网格数量(nx,ny,nz)以及各单元体密度值,观测面数据包括测点总数(ns)及每个测点的经度、纬度、高度。② 根据模型深度轴层数nz逐层遍历;根据观测面测点数ns逐点遍历;进一步,每一层根据模型规则网nx×ny个单元体按照先y轴方向、再x轴方向依次遍历计算。按照上述遍历顺序,利用GLQ2D积分法计算单个单元体对单个观测点的重力效应,然后叠加单层所有单元体对单点重力效应,最后叠加所有层单元体对单点重力效应,也即得到模型在观测面各测点的重力异常。③ 输出该模型对观测面重力正演结果。图3是球坐标系的多层Tesseroid单元体模型剖分示意。

图3

图3 三维Tesseroid网格模型

Fig.3 3-D Tesseroid mesh model

本文的并行思路采用基于OpenMP的MATLAB任务并行算法,即利用parfor对for循环进行并行(图4)。MATLAB进行parfor循环时采用client和worker模式。其中client为编写和启动该代码的MATLAB端,而worker指运行该代码的MATLAB端。电脑中的MATLAB软件可视为一个进程,同一台电脑可以同时运行多个MATLAB进程,每个worker对应的物理单元为处理器或处理器核。每个MATLAB进程之间可以通过一定方式开展数据传输。用户首先在client端编写所要运行的代码,client端在运行该代码的过程中,将需要并行的代码段分配到其他MATLAB进程运行。

图4

图4 MATLAB parfor任务并行算法示意

Fig.4 MATLAB parfor task parallel algorithm schematic

对图2所示的球坐标系三维网格模型重力正演计算的并行策略可以有很多种,但效果基本相当。本文根据最外层(模型深度层)循环并行加速最优的原则,提出并行策略如下:首先明确需要并行的代码段,即模型分层的循环结构,然后通过任务并行加速在client端实现代码设计,MATLAB client端自动将模型不同层使用算法分配给不同worker处理单元,各处理单元分别计算单层重力效应,最终叠加得到模型正演结果。

2 SLAB模型试验

2.1 模型参数

该系列模型为俯冲带简单模型模拟,经度、纬度范围分别为100°E~110.2°E,20°N~30.2°N,深度方向范围为0~200 km,俯冲带用以模拟板块汇聚边缘,其顶部距观测面的垂直距离为30 km,并以东倾60°的方向一直延伸至200 km,剩余密度为0.08 g/cm³。

本文设计了3个不同数量级的模型:① 模型大小51°×51°×20 km,经度、纬度方向网格间距为0.2°×0.2°,深度方向步长10 km;② 模型大小51°×51°×40 km,经度、纬度方向网格间距为0.2°×0.2°,深度方向步长5 km;③ 模型大小101°×101°×40 km(图5),经度、纬度方向网格间距为0.1°×0.1°,深度方向步长5 km。3个模型观测面高度统一设置为0 m。

图5

图5 俯冲带三维模型示意

Fig.5 Three dimensional model of subduction zone

2.2 模型正演与效果评价

本文在集群机上实现基于CPU的并行计算,集群CPU信息:Intel(R) Xeon(R) CPU E5-2680 v2 @2.80GHz,并行使用软件为:MATLAB R2018a,模型并行计算单元数即并行使用CPU核数:12。图6为GLQ2D串行正演结果,并行正演结果,以及模型24°N剖面结果对比。通过对比可见并行算法与串行算法正演结果完全一致。

图6

图6 GLQ2D串行和并行正演结果对比

注:白色虚线方框为俯冲带沿倾角方向投影到观测面上位置;黑色虚线为对比剖面位置

Fig.6 The comparison results between the GLQ2D serial and parallel forward modeling

Note: the white dotted box shows the position of the subduction zone projected onto the observation plane along the dip angle;the back dotted line is the position of the comparison profile

对比图6a、b,从数值来看,其异常结果均呈现中间异常高,两侧异常低,且俯冲带右侧异常降低较左侧慢的特征。然后进一步对比图中24°N剖面正演曲线,如图6c所示,重力异常中间高两侧低,俯冲带左侧异常降低较右侧快,重力异常串行结果和并行结果完全吻合,验证了对于该数据模型使用并行计算不会改变它的计算结果和计算精度。

表1与表2分别对比了不同核数并行俯冲带模型重力正演计算效率以及不同数量级俯冲带模型的重力正演计算效率对比。在这里为了解释方便,引入加速比的定义,即串行运行用时比并行运行时间^[17]。

$S = \frac{T_{串行}}{T_{并行}}$

表1 不同并行核数重力正演计算效率对比

Table 1 The efficiency comparison results of gravity forward modeling with different parallel cores

模型大小	核数	时间/s	加速比
101×101×40	串行	11432.970
	2核	5871.609	1.947
	4核	2882.673	3.966
	6核	2033.522	5.622
	8核	1704.662	6.707
	10核	1441.482	7.931
	12核	1193.787	9.577

新窗口打开| 下载CSV

表2 不同数据量模型重力正演计算效率对比

Table 2 The efficiency comparison results of gravity forward modeling with different data volume models

数据量	并行时间/s	串行时间/s	加速比
51×51×20	86.243	529.506	6.159
51×51×40	147.773	1234.136	8.372
101×101×40	1193.787	11432.970	9.580

新窗口打开| 下载CSV

从表2可以看出并行计算相比串行计算有明显且稳定的速度优势,且随并行核数的增加,加速比不断增加,在12核并行计算时得到最大平均加速比9.580。因此本文中其他正演算法均采用12核并行计算。

对一个并行系统(包含但不限于算法,程序)来说,如果其性能(如加速比)可以随处理单元数目的增加而按比例增加,我们称该并行系统具有可扩放性^[17]。对于球坐标系MATLAB重力异常正演并行算法,保持程序的可扩放性是不可缺少的。对比表1中不同并行核数模型计算运行时间及其加速比以及表2中不同数据量模型计算运行时间及其加速比,可以得知该并行策略稳定可行,即具有良好的可扩放性,且随计算数据量的增加,计算加速比不断增加并趋于处理单元数。

3 华东三维岩石圈模型试验

本次正演试验所选取的研究区为中国东部,经度范围99.75°E~122.25°E,纬度范围20.75°N~45.25°N。华东是我国地质构造和动力学研究及资源能源勘查的重要区域,其区域地球物理场数值模拟和壳幔结构成像具有重要的科学意义。本文从Shen等^[18] 中提取出研究区岩石圈三维横波速度结构模型,根据Brocher^[19]的速度V_p、V_s和密度经验式换算得到研究区岩石圈三维密度模型(图7)。该华东岩石圈三维密度模型网格大小111°×121°×95 km,经度、纬度方向网格间距为0.2°×0.2°,深度方向步长 2 km。最小密度值为2.26 g/cm³,最大密度值为 3.45 g/cm³,平均密度为3.144 g/cm³。针对该模型,传统的直角坐标系计算受地球曲率的影响会产生明显的误差,因此需要球坐标系计算。而该模型的球坐标系正演串行算法耗时较长(约15 h),因此,本文采用并行算法来提高正演效率。

图7

图7 华东岩石圈三维密度模型

Fig.7 Three dimensional density model of Eastern China lithosphere

应用本文方法对华东岩石圈三维模型进行了正演计算,并截取模型30°N剖面进行效果对比。图8显示了华东模型GLQ2D串行正演结果和12核并行正演结果,以及其30°N剖面的正演效果对比。通过对比可见模型并行正演与串行正演结果完全一致。研究区岩石圈密度模型的理论重力异常变化整体较为平缓,幅值范围为 -450~650 mGal。其中,内蒙古中东部和河北北部及青藏高原一带呈现大范围的低值分布,扬子克拉通呈现大范围的高值分布,鄂尔多斯地块、松辽盆地呈现中等范围的高值分布,东南沿海一带异常幅值整体在-100~200 mGal之间。

图8

图8 华东岩石圈密度模型正演结果对比

注:黑色虚线为对比剖面位置

Fig.8 The comparison results of the Eastern China lithosphere density model

Note:the back dotted line is the position of the comparison profile

本文收集研究区EGM2008^[20]的实际布格重力异常数据(图8d)用于异常比较,数据网格间距为0.2°×0.2°。研究区EGM2008实际布格重力异常总体呈东高西低的趋势,幅值范围为-350~200 mGal。本文正演结果的整体异常特征与实际布格重力异常差别较大,究其原因主要有两点:① 本文正演所依据的模型为0~200 km深度范围内的华东岩石圈三维密度模型,不包括周边区域及200 km深度以下的介质模型;② 本文所用的岩石圈三维密度模型是由背景噪声成像获得的岩石圈三维横波速度模型^[18]通过速度—密度经验式^[19]换算而来的,该模型没有得到重力数据约束,而且所用的经验式并不能准确反映华东的实际速度和密度关系,不能反映研究区真实壳幔密度分布,从而导致上述的异常差别。

表3为华东模型运行时间表,从表3可以看出,并行计算相比串行计算有明显且稳定的速度优势,且对比理论模型,该模型串行运行时间均值比上12核并行计算运行时间均值得到平均加速比为 11.206,远大于前者的9.580。这进一步验证了该并行系统的可扩放性,即随计算数据量的增大,并行加速比随之增大且趋于处理单元数。

表3 重力正演计算效率对比

Table 3 The efficiency comparison results of gravity forward modeling

数据量	并行时间/s	串行时间/s	加速比
111×121×95	4841.807	54255.090	11.206

新窗口打开| 下载CSV

4 结论

本文针对大尺度、地表观测面的高精度、高效率重力正演问题,给出了基于改进的高斯—勒让德积分法的并行计算方案,并通过理论模型和华东岩石圈三维模型数据试验验证了本文方法的有效性。数据试验表明本文方法可以实现大尺度、地表观测面重力异常的高精度、高效率正演,基于OpenMP的MATLAB任务并行算法具有可扩放性及稳定性,即随数据量的增大,加速比会趋于并行单元数;数据量不变,并行单元增大,加速比也会增大。本文方法为高效的大尺度重力场模拟和三维反演奠定技术基础。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Heck

, Seitz

A comparison of the Tesseroid, prism and point-mass approaches for mass reductions in gravity field modelling

[J]. J. Geodesy, 2007, 81(2):121-136.

DOI:10.1007/s00190-006-0094-0 URL [本文引用: 5]

[2]

杨学祥

布格改正和地形改正的误差——关于区域重力测量中地形改正最大半径的讨论

[J]. 地壳形变与地震, 1992(2):1-6.

Yang X

Errors in Bouguer and topographic correction—Discussion on the maximum radius of topographic correction in regional gravity measurement

[J]. Crustal Deformation and Earthquake, 1992(2):1-6.

[3]

梁青

月球重力异常特征与三维密度成像研究

[D]. 武汉:中国地质大学(武汉), 2010.

Liang

Gravity anomaly features and 3D density imaging of the moon

[D]. Wuhan: China University of Geosciences (Wuhan), 2010.

[4]

Grombein

, Seitz

, Heck

Optimized formulas for the gravitational field of a Tesseroid

[J]. J.Geodesy, 2013, 87(7):645-660.

[5]

Uieda

, Barbosa V C

Fast nonlinear gravity inversion in spherical coordinates with application to the South American Moho

[J]. Geophys. J. Int., 2017, 208(1):162-176.

DOI:10.1093/gji/ggw390 URL [本文引用: 2]

[6]

Asgharzadeh M

Von Frese R R

Kim H

Spherical prism magnetic effects by Gauss-Legendre quadrature integration

[J]. Geophys. J. Int., 2008, 173(1):315-333.

DOI:10.1111/gji.2008.173.issue-1 URL [本文引用: 2]

[7]

Zhang

, Wu Y

, Yan J

, et al.

3D inversion of full gravity gradient tensor data in spherical coordinate system using local north oriented frame

[J]. Earth, Planets and Space, 2018, 70(58):1-23.

DOI:10.1186/s40623-017-0766-4 URL [本文引用: 1]

[8]

王祥, 郭良辉.

球坐标系密度界面反演方法及在华南大陆的应用

[J]. 物探与化探, 2020, 44(5):1161-1171.

[本文引用: 7]

Wang

, Guo L

Density interface inversion method in spherical coordinate and its application in South China mainland

[J]. Geophysical and Geochemical Exploration, 2020, 44(5):1161-1171.

[本文引用: 7]

[9]

Hao A

, Guo L

, Wang

The apparent density mapping approach in spherical coordinates and the crustal density distribution of Chinese mainland

[J]. IEEE Access, 2019, 7(1):160705-160717.

DOI:10.1109/Access.6287639 URL [本文引用: 3]

[10]

Cui Y

, Guo L

A wavenumber-domain iterative approach for 3D imaging of magnetic anomalies and gradients with depth constraints

[J]. Journal of Geophysics and Engineering, 2019, 16(6):1032-1047.

DOI:10.1093/jge/gxz064 URL [本文引用: 1]

[11]

Zhao G

, Chen

, Uieda

, et al.

Efficient 3D large-scale forward-modeling and inversion of gravitational fields in spherical coordinates with application to lunar mascons

[J]. Journal of Geophysical Research: Solid Earth, 2019, 124(4):4157-4173.

DOI:10.1029/2019JB017691 URL [本文引用: 1]

[12]

陈召曦, 孟小红, 郭良辉, 等.

基于GPU 并行的重力、重力梯度三维正演快速计算及反演策略

[J]. 地球物理学报, 2012, 55(12):4069-4077.

Chen Z

, Meng X

, Guo L

, et al.

Three-dimensional fast forward modeling and the inversion strategy for large scale gravity and gravimetry data based on GPU

[J]. Chinese Journal of Geophysics, 2012, 55(12):4069-4077.

[13]

Hou Z

, Wei X

, Huang D

, et al.

Full tensor gravity gradiometry data inversion: performance analysis of parallel computing algorithms

[J]. Applied Geophysics, 2015, 12(3):292-302.

DOI:10.1007/s11770-015-0495-z URL [本文引用: 2]

[14]

Hou Z

, Huang D

Multi-GPU parallel algorithm design and analysis for improved inversion of probability tomography with gravity gradiometry data

[J]. Journal of Applied Geophysics, 2017, 144:18-27.

DOI:10.1016/j.jappgeo.2017.06.009 URL [本文引用: 2]

[15]

Hou Z

, Huang D

, Wang E

, et al.

3D density inversion of gravity gradiometry data with a multilevel hybrid parallel algorithm

[J]. Applied Geophysics, 2019, 16(2):141-152.

DOI:10.1007/s11770-019-0763-4 URL [本文引用: 2]

[16]

周雪, 于平, 翁爱华, 等.

基于MPI和OpenMP的重力及重力梯度数据并行正演算法研究

[J]. 世界地质, 2018, 37(3):897-904.

Zhou

, Yu

, Weng A

, et al.

Parallel forward modelling algorithm with gravity and gravity gradient data based on MPI and OpenMP

[J]. Global Geology, 2018, 37(3):897-904.

[17]

陈国良

. 并行计算:结构·算法·编程[M]. 北京: 高等教育出版社, 2003:83-88.

Chen G

. Parallel: construction· algorithms· programming [M]. Beijing: Higher Education Press, 2003:83-88.

[18]

Shen W

, Michael H

, Kang

, et al.

A seismic reference model for the crust and uppermost mantle beneath China from surface wave dispersion

[J]. Geophys. J. Int., 2016, 206(2):954-979.

DOI:10.1093/gji/ggw175 URL [本文引用: 2]

[19]

Brocher T

Empirical relations between elastic wavespeeds and density in the Earth’s crust

[J]. Bull. Seismol. Soc. Am., 2005, 95(6):2081-2092.

DOI:10.1785/0120050077 URL [本文引用: 2]

[20]

Pavlis N

, Holmes S

, Kenyon S

, et al.

The development and evaluation of the Earth Gravitational Model 2008 (EGM2008)

[J]. J. Geophys. Res., 2012, 117(4):B04406.