王博，郭良辉,2，崔亚彤，王祥

(1.中国地质大学(北京)地球物理与信息技术学院，北京 100083；2.地质过程与矿产资源国家重点实验室 中国地质大学(北京)，北京 100083)

0 引言

三维网格模型的正反演是重力资料立体解释的重要手段。其中，正演是反演的基础，高精度、高效率的正演将促进高质量、高效率的反演解释。对于小范围、小尺度研究区，通常将地下三维空间剖分为众多规则排列的直立长方体组合，各长方体的密度各不相同，然后对直立长方体组合的三维网格模型进行正演和反演。然而，对于大面积、大尺度研究区，由于地球曲率的存在，常规的基于直立长方体组合的三维网格模型已不能适用，需采用球坐标系的三维网格模型方法。

当前球坐标系重力异常正演方法较多，主要有下面4类：① 将球坐标系下的单元体以点元、线元或者面元的方式对近似的规则形体进行处理[1]，这类方法主要的问题是计算精度非常低，难以满足日益增加的精度要求；② 将全球地形剖分为众多扇形环或球冠单元体等规则形体的组合，并给出严格的积分解析式，然后计算球坐标系下各规则模型体的正演重力效应，再累加作为全球地形改正值[2]，这类方法一般用于单点观测，且随着观测点的变化，相对应的球坐标系以及地形模型需要重新构建；③ 基于Tesseroid单元体的泰勒级数展开法，在正演模型体的几何中心点对模型积分核函数进行泰勒展开，计算球坐标系下各Tesseroid单元体的重力效应并叠加[1,3-4]，这种方法的计算精度与泰勒展开级数以及剖分网格的大小密切相关；④ 基于Tesseroid单元体的高斯—勒让德(GLQ)积分法，给出椭球形地球的精确剖分方案，然后将基于数值积分方法的GLQ积分法应用于椭球体积分计算中，同样通过叠加得到最终重力效应[5-7]，这种方法计算精度一般高于泰勒级数展开法。然而，不管是高斯—勒让德积分法还是泰勒级数展开法，它们在正演计算地表或近地表观测面的测点异常时，都会产生邻近效应，即误差会随着观测距离的减小而急剧增大。王祥等和Hao等就此对高斯—勒让德积分法作了改进,给出了新的单元体自适应剖分方案，极大提高该算法在地表观测面的计算精度，进而分别给出了基于改进高斯—勒让德积分法(GLQ2D)的单一密度层界面深度反演和视密度填图方法[8-9]。

本文研究将王祥等和Hao等[8-9]的单一密度层模型正演方法扩展到多个密度层的三维网格模型，为大尺度数值模拟和三维反演奠定基础。由于GLQ2D算法细化每个Tesseroid单元体剖分，导致剖分后的小单元体个数增加及正演计算量指数级增大，造成大规模三维Tesseroid网格模型的正演效率低下。

在重磁领域，国内外学者给出了重力异常正反演的一些快速算法，主要有频率域算法[10]，平移对称算法[11]和并行计算[12-16]等。其中，并行计算主要分为两大方向，主机端的CPU并行和设备端的GPU并行。CPU端的并行主要有OpenMP(open multi-processing)和MPI(message passing interface)等方式。OpenMP是一种用于共享内存并行系统的多线程程序设计的一套指导性注释，它可以使程序员把更多的精力投入到并行算法本身，而非其具体实现细节。而MPI则是信息传递接口，是独立于语言的通信协议，它需要程序员手动管理数据分配，实现进程通信以及维持同步。GPU端的并行方式主流的有NIVIDIA CUDA架构。按并行策略的不同，并行算法也分为两种：一种是任务并行，另一种则是数据并行。任务并行的并行思想是把同一批数据分给for循环中不同的循环体，进行处理。数据并行的思想则是不同的数据，用同一个程序处理。当前，国内外一些学者应用并行算法一定程度上提高了重磁正反演效率。比如，陈召曦等[12]给出了基于NIVIDIA CUDA并行计算的重力及重力梯度数据正演算法，Hou等[13-15]给出了基于多级混合并行的重力梯度三维密度反演算法，周雪等[16]给出了基于MPI和OpenMP并行计算的重力及重力梯度正演算法。但这些算法都是基于直角坐标系进行的。

本文针对大尺度、地表观测面研究区的高精度、高效率重力正演问题，研究球坐标系三维Tesseroid网格模型重力异常正演算法，采用GLQ2D算法提高大尺度、地表观测面的重力异常计算精度，采用基于OpenMP的MATLAB任务并行算法提高正演效率。最后通过理论模型数据和华东岩石圈三维密度网格模型作正演试验，评价本文方法的有效性。

1 方法原理

基于球坐标系的三维网格模型重力异常正演方法是将地下三维空间剖分为规则排列的Tesseroid单元体组合，各个单元体具有不同的密度值，再通过GLQ2D算法正演计算每个单元体在观测面每个测点的重力异常，观测面上任意测点的重力异常值即是所有Tesseroid单元体重力异常值在该测点处的累加。

1.1 Tesseroid单元体重力正演公式

在球坐标系中，对于任意Tesseroid单元体(图1)，其重力异常积分公式可由下式表示[1,6]:

图1 Tesseroid单元体示意[8]

(1)

其中：

cosψ=sinφosinφs+cosφocosφscos(φs-φo)

式中：G表示万有引力常量，G=6.67×10-11N·m2/kg2，ρ表示模型体密度，单位为g/cm3，观测点坐标Q0(r0,φ0,λ0)，Tesseroid单元体中心点坐标Q(rs,φs,λs)，l为观测点Q0与单元体中心点Q之间的距离，ψ为Q0点与Q点相对于球心的夹角。

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

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

(2)

式中：

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

(3)

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

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

(4)

其中：

L=max(Lφ,Lλ)

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

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

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

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

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

图3 三维Tesseroid网格模型

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

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

对图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/cm3。

本文设计了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 俯冲带三维模型示意

2.2 模型正演与效果评价

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

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

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

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

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

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

从表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]的速度Vp、Vs和密度经验式换算得到研究区岩石圈三维密度模型(图7)。该华东岩石圈三维密度模型网格大小111°×121°×95 km，经度、纬度方向网格间距为0.2°×0.2°，深度方向步长 2 km。最小密度值为2.26 g/cm3,最大密度值为 3.45 g/cm3，平均密度为3.144 g/cm3。针对该模型，传统的直角坐标系计算受地球曲率的影响会产生明显的误差，因此需要球坐标系计算。而该模型的球坐标系正演串行算法耗时较长(约15 h)，因此，本文采用并行算法来提高正演效率。

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

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

本文收集研究区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 重力正演计算效率对比

4 结论

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