一种新的重力正演加速方法与流程

本发明涉及一种地球物理领域的重力正演方法，特别涉及大规模航空重力、卫星重力等的高精度、快速正演方法。

背景技术：

重力勘探是地球物理勘探的一种重要的方法，并被广泛应用到矿产资源勘探、地壳结构与莫霍界面探测和环境地球物理勘查中。测量的重力数据与地下密度分布体有着直接的关系，因此，通过对地下密度进行反演可以用来进行地下密度异常体的判别。随着重力勘探技术的发展，航空重力技术的逐渐成熟，对于给定密度分布的目标体，重力勘探往往含有较多的观测点，传统的重力正演方法计算耗时的特征越来越明显。因此，高效的正演显得尤为重要，目前，常用的重力计算加速算法主要包括两种：基于快速傅里叶变换(Fast Fourier Transformation,FFT)的重力正演和基于小波变换(wavelet transformation)的重力快速正演。以上两种方法计算效率快、精度高，但是都需要对三维地质体进行均匀的六面体网格剖分，这使得网格剖分不能很好的拟合复杂模型，方法不适合处理带地形的三维问题，同时该方法只能计算观测点在同一个平面上的重力正演。非结构化网格能够拟合复杂模型的优势日益突出，被广泛应用在三维数值模拟中。

因此，有必要设计一种高效、高精度的基于非结构化网格剖分的观测点任意分布的重力正演方法。

技术实现要素：

本发明所解决的技术问题是，针对现有技术的不足，提供了一种新的重力正演加速方法，能够基于非结构化网格剖分，对任意分布的观测点进行重力正演计算，大幅度地提高了重力正演计算的效率。

本发明的技术方案为：

步骤1、模型建立、区域离散与八叉树建立：

如图1所示，三维地质体Ω的密度分布为ρ(r′)，在观测点r处产生的重力位φ(r)、重力梯度g(r)和重力梯度张量T_g(r)可以表示为：

其中，为梯度运算符，R＝|r-r′|为观测点r到源点r′的距离；dv表示体积分，ρ(r′)为三维地质体Ω在源点r′处的密度，G为引力常数，G＝6.674×10^-11m³kg^-1s^-2；

模型建立：根据观测点与三维地质体(又称为源)的特征，设计三维地质体模型参数(模型坐标及密度信息)和观测点参数(观测点坐标)。

区域离散：对三维体地质体模型采用四面体单元进行网格离散，将三维地质体离散成多个四面体单元T_i,i＝1,2,...,N，N为四面体单元的个数。因此，三维地质体Ω可表示为：

八叉树建立：八叉树是一种用来描述三维空间的树状数据结构，本发明采用八叉树的数据构建方法，实现对观测点以及地下地质体建立相应的八叉树数据结构，即观测点八叉树和三维地质体单元(源)八叉树，如图3所示；

对源建立八叉树：构建一个立方体单元包含三维地质体，并将该立方体单元称第0层细胞；然后以第0层细胞为母细胞，将其细分成8个子细胞，该8个子细胞称为第1层细胞；仍然，按照上述剖分方式对每个子细胞进行细化，后续每剖分一层得到的子细胞依次记为第j(j＝2,3,4…,q)细胞层子细胞，直到每个子细胞里包含的四面体地质单元数少于设定的个数为止。最后一层的子细胞称为叶子细胞，同时删除不含有四面体地质单元的细胞。

同样地，构建观测点的八叉树数据结构：构建一个立方体单元包含所有的测点，称为第0层细胞；然后以第0层细胞为母细胞，将其细分成8个子细胞，该8个子细胞称为第1层细胞；仍然，按照上述剖分方式对每个子细胞进行细化，后续每剖分一层得到的细胞依次记为第i(i＝2,3,4…,p)层细胞直至每个细胞包含的观测点数目少于设定的个数为止，同时删除不含有观测点的细胞，最后一层的子细胞称为叶子细胞，。

步骤2、近区与远区的分离：

常规的重力计算方法是每个观测点对每个小的四面体地质单元进行一次积分，然后叠加求和得到该测点的重力值。快速多极展开算法的目的就是通过将积分区域划分成远区和近区，近区积分与常规的重力计算方法一致，对于远区的积分，通过选取扩展点，即每个细胞的中心点，将每个观测点与所有的四面体地质单元的积分转换成观测点所在细胞的中心点与四面体单元所在细胞的中心点的连接，从而减少了积分计算的次数，达到加速的目的。

为了达到这一加速的目的，将观测点与三维地质体之间的连接分为三种连接关系：第一种为三维地质体子细胞中心点与母细胞中心点的连接，通过M2M转换连接；第二种为三维地质体细胞中心点与观测点母细胞中心点的连接，通过M2L转换连接；第三种为观测点子细胞中心点与母细胞的中心点的连接，通过L2L转换连接。如图3所示。

从重力位的计算公式可以看出，积分核函数含有1/R项，当R趋近于0时，积分含有奇异性；当观测点r与三维地质体Ω内部积分点r′的距离趋近于0时，即R趋近于0，因此，将含有奇异性积分的区域称为近区，反之，其他区域称为远区。但是，实际计算过程中，对于每一个观测点r，若三维地质体Ω中某个四面体单元的中心点到观测点r的距离小于某一设定的数值ε【ε趋近于0，可取1-4个四面体单元的长度】，则认为该四面体单元所在区域为近区，记为T_near；其他区域为远区，记为T_far。如图2所示。那么，区域Ω可以表示为：

Ω＝T_near+T_far (2)

步骤3、分别计算近区与远区重力参数：

将重力正演的积分划分为对近区的积分和对远区的积分；基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区的积分，得到近区重力参数；并基于观测点八叉树和源八叉树，采用自适应快速多极展开算法计算远区的积分，得到远区重力参数；然后通过求和得到重力正演结果。

近区重力参数包括在观测点r处产生的近区重力位、近区重力梯度和近区重力梯度张量，分别记为远区重力参数包括在观测点r处产生的远区重力位、远区重力梯度和远区重力梯度张量，分别记为

步骤4、将步骤3中计算得到的近区重力参数与远区重力参数求和，得到观测点重力参数，包括在观测点r处产生的重力位φ(r)、重力梯度g(r)和重力梯度张量T_g(r)，其计算公式分别为：

以下对所述采用自适应快速多极展开算法计算远区的积分，得到远区重力参数的原理和步骤进行详细说明：

为了使用快速多极展开算法，需要对积分点进行转换，即对积分核函数1/|r-r′|进行展开，当满足条件时，积分核函数可以展开为：

其中，为源点r′所在叶子细胞的中心点(r′的扩展点)，R_nm和S_nm分别为常规的(R)和奇异的(S)球谐函数，为S_nm的复共轭。这两种球谐函数可以表示为：

其中，(r,θ,φ)为点x的球坐标表示，r>0，0≥θ≥π，0≥φ≥2π。为伴随勒让德函数，定义为：

其中，P_n为标准的n阶勒让德多项式。

此时，将积分核函数分成了两部分，一部分是关于观测点r和的函数另一部分为关于r′和的函数R_nm，目的是为了将观测点r从积分中分离出来。那么，重力位可以表示为：

其中，为重力位在处的多极距，叶子细胞的多极距直接通过多极距公式计算。多极距只与r′和有关，与观测点r无关，不随观测点的变化而变化，这就意味着，即使有很多个观测点，多极距也只需要计算一次。如果多极距一旦被计算好，当计算新的观测点的远区积分时，只需要进行一些简单的代数操作，即计算依赖于观测点的基础函数与不依赖于观测点的多极距的乘积。对于大规模重力计算，该方法可以大幅度的减少计算时间。如果所有的源细胞都在远区，那么计算M个观测点的φ(r)_far的时间复杂度将从O(NM)减少为O(bM+N)，其中，N为四面体单元的个数，b为细胞数。为了达到进一步减少计算时间的目的，需要借助自适应多层的方法，此时，时间复杂度可以减少为O(M log N)。

设点r_cⁱ为源八叉树第i层细胞的中心点，点为源八叉树第(i+1)层细胞的中心点，点r_cⁱ所在细胞为点所在细胞的母细胞；其中i＝0,1,2…,p-1；最高层细胞即第0层细胞的中心点为r_c⁰；对于母细胞，通过M2M转换公式，对该母细胞的所有子细胞累加求和得到母细胞的中心点处的多极距；

首先，通过点处的多极距与以下M2M转换公式计算点r_cⁱ处的多极距

然后，通过以下公式计算远区重力位：

(8)

母细胞的多极距通过其子细胞的多极距和M2M转换公式计算，M2M转换公式实现了子细胞与母细胞之间的联系，同样包含两部分，一部分为与和有关的另一部分为与r′和有关的即为子细胞的多极距。可以看出，通过M2M转换，点处的多极距只需要进行简单的几何运算即可求得，从而可以快速求得远区重力位。

类似地，如果点为观测点的局部扩展点，当满足条件时，积分核函数可以展开为：

那么，重力位在远区的积分可以表示为：

其中，为已知的基础函数，表示处的局部系数，为未知的局部系数，可以通过计算远区T_far的积分得到。但是，也可以通过已经计算好的多极距来计算局部系数从而达到减少计算量，加快计算速度的目的。

借助源细胞中心点的多极距来计算扩展点的局部系数具体计算方法为：

首先，对进行展开：

这里，此时，将展开成两部分，一部分为与r′和相关的另一部分为与和相关的可以看出，与观测点及其扩展点无关，很好的将观测点与源点分离开。由此得到M2L转换公式为：

由以上公式求得处的局部系数。

设点为观测点八叉树第(j+1)层细胞的中心点，点为观测点八叉树第j层细胞的中心点，点所在细胞为点所在细胞的母细胞，j＝0,1,2…,q-1；通过点处的局部系数和以下L2L转换公式计算点处的局部系数

那么，点的局部系数可以通过简单的代数求和直接计算得到。使用上述L2L转换，可以最大程度的减少计算时间。

然后，通过以下公式计算远区重力位：

逐层计算，从上往下，便可以计算得到每个观测点的重力位的远区积分。

基于同样的原理，可以快速计算远区重力场和远区重力梯度张量：

通过以下公式计算远区重力场和远区重力梯度张量：

用于计算远区重力场和远区重力梯度张量的局部系数与计算远区重力位的局部系数相同。这就意味着，如果计算远区重力位时，点处局部系数通过L2L，M2L和L2L转换公式计算得到，只需要再计算局部系数和局部基本的函数和的乘积，便可以计算得到重力场和重力梯度张量。

计算局部基础函数和需要使用坐标转换，这里x_i(i＝1,2,3)表示三个坐标分量，表示三个方向的单位向量。那么，可以写为：

可以写为：

其中，y_i(i＝1,2,3)表示为球坐标系下三个分量，y₁＝r>0，y₂＝θ和0≤θ≤π，

x_i和y_i之间的关系可以表示为：

所述步骤3中，采用基于近区所包含的四面体单元的重力积分无奇异性解析表达式直接计算近区重力参数的计算公式为：

其中，R＝1/|r-r′|。

本专利采用基于四面体单元的积分解析计算的方法对上述几种积分进行无奇异性、高精度的解析计算。

对于任意的四面体单元T，借助矢量恒等式，

标量形式可以重新写为：

其中，为四面体单元T的第i个三角面，d_i＝n_i·(r′-r)表示向量(r′-r)在面的法向n_i上的投影。当点r靠近点r′时，面积分k_i包含对1/R的含有弱奇异性的积分。

借助梯度定理，向量形式可以表示为：

其中，k_i定义和前面一致，因此，如果k_i计算了，那么体积分A和D便可以计算得到。

可以表示为：

下面，给出面积分和的解析表达式。下面使用K表示三角面用k和k^-1分别表示含有奇异性的面积分和这里，|r-r′|表示点r到点r′的距离。

如图5所示，点o为测点r在面K的投影点，并设点o为局部极坐标系的原点。为面K的面外法向，那么，有：

其中，ω₀可以为任意实数，ω₀＝0表示点r在面K上。为极坐标系下的向量，从点o指向积分点r′，可以写为其中单单位向量，ρ≥0为向量的长度。此时，R＝|r-r′|可以表示为的奇异性取决于ω₀和ρ的数值，当ω₀＝0时，在点o附近的无穷小得区域里积分存在奇异性，因此，将积分区域K分为两个部分：可能含有奇异性的无穷小区域C(ε)和其他不含奇异性的区域K-C(ε)，那么有：

K＝[K-C(ε)]+C(ε) (26)

无穷小区域C(ε)为一个角度为α半径ε→0的圆形区域。那么α可能为：

借助如下两个向量恒等式(Wilton,et al,1984)：

其中，为对点r的梯度，表示对积分点r′的面梯度。那么，有：

在积分区域K-C(ε)，R^-1无奇异性，借助梯度定理：

其中，为在边上的投影在外法向上的投影，在边上，为常数。类似地，

其中，积分是与点o相关的奇函数。积分只有当观测点在三角面上时才消失，此时a(o)＝2π。当观测点在三角面的边上或者点上时，积分存在奇异性。因此，有：

现在，考虑两个一维线积分。建立一维局部坐标系，定义：

其中，为边的切向方向，为面K的面外法向，为边的外方向。s＝0表示向量与边垂直，其值可以表示为：

这里，和可以为正数，也可以是复数，表示点o在边的延长线上，表示点o在边上。现在，距离R可以表示为这里等于向量的方向与方向的点积。和为点r到点和点的距离。和为点r到边的距离。

现在，一维线积分可以表示为：

那么，

这里，角度a(o)可以表示为：

当观测点在边上时，将使数值计算不稳定，为了避免该问题，引入此时，有：

其中，此时，面积分没有奇异性。

从上可知，其中的第二个积分有如下结果：

当时，第一部分可以写为(Gradshteyn,et al,2014)：

当时，如果积分存在奇异性，从新计算积分

将近区与远区对观测点r处产生的重力位、重力梯度和重力梯度张量求和得到整个区域对观测点r处产生的重力位、重力梯度和重力梯度张量，其计算公式分别为：

有益效果：

常规的重力计算方法是每个观测点对每个小的四面体地质单元进行一次积分，然后叠加求和得到该测点的重力值。本发明采用快速多极展开算法，通过将积分区域划分成远区和近区，近区积分与常规的重力计算方法一致，对于远区的积分，将每个观测点与所有的四面体地质单元的积分转换成观测点所在细胞的中心点与四面体单元所在细胞的中心点的连接，从而减少了积分计算的次数，实现加速的目的。为了达到这一加速的目的，分别建立观测点与三维地质体八叉树结构，将观测点与三维地质体之间的连接分为三种连接关系：第一种为三维地质体子细胞中心点与母细胞中心点的连接，通过M2M转换连接；第二种为三维地质体细胞中心点与观测点母细胞中心点的连接，通过M2L转换连接；第三种为观测点子细胞中心点与母细胞的中心点的连接，通过L2L转换连接。计算三维地质体八叉树中各个细胞的多极距或观测点八叉树中各个细胞的局部系数，从而计算远区重力位；三维地质体八叉树中叶子细胞的多极距直接通过多极距公式计算，母细胞的多极距通过M2M转换公式计算；多极距与观测点r无关，不随观测点的变化而变化，这就意味着，即使有很多个观测点，多极距也只需要计算一次。如果多极距一旦被计算好，当计算新的观测点的远区积分时，只需要进行一些简单的代数操作。对于大规模重力计算，该方法可以大幅度的减少计算时间。进一步地，本发明也可以通过计算观测点八叉树中各个细胞的局部系数，从而计算远区重力位；观测点八叉树中最高层细胞的局部系数通过三维地质体八叉树中最高层细胞的多极距和M2L转换公式计算，再通过L2L转换公式计算子细胞的局部系数。

本发明可以进行大规模重力高精度、快速正演计算与大规模重力地形改正计算，大幅度地提高了计算的效率。

附图说明

图1为典型的任意复杂的重力勘探示意图，r表示观测点，r'表示三维地质体Ω内的任意一点，三维地质体的密度分布为ρ(r′)，h表示观测点相对地面的高度。

图2为三维体Ω分成近场区T_near和远场区T_far示意图。

图3为观测点八叉树和源八叉树中M2M、M2L和L2L转换示意图。

图4为基于快速多极展开算法的重力场计算方法的流程图。总场通过对近区场和远区场求和得到，近区场直接采用解析公式计算，远区场通过快速多极展开算法计算。

图5为基于四面体单元无奇异性积分局部坐标系示意图；

图6为六面体模型示意图。六面体的密度为2000kg/m³，并被剖分成一系列的四面体单元，观测平面距离六面体的上表面的距离为h。

图7为本发明的基于四面体单元的解析表示式计算结果与基于六面体单元的解析表达式计算结果的对比图。图(a)为基于六面体单元的解析表达式计算结果；图(b)为本发明的基于四面体单元的解析表示式计算结果；图(c)为二者的相对误差图。

图8为本发明的基于快速多极展开算法的重力场计算结果与基于四面体单元的解析表达式计算结果的对比图。图(a)为本发明的基于四面体单元的解析表达式计算结果；图(b)为本发明的基于快速多极展开算法的重力场计算结果；(c)为二者的相对误差图。

具体实施方式

以下结合附图和具体实施方式对本发明作进一步的说明。

本发明涉及的重力计算方法包括以下步骤：

步骤1、测点和源文件的设计：确定测点的个数与坐标、密度体位置坐标文件。并进行非结构化网格剖分。分别对观测点与密度体建立对应的八叉树结构，将积分区域分为近区与远区。

步骤2、远区场计算：通过自适应快速多极展开加速算法计算远区场积分。

步骤3、近区场计算：通过解析表达式计算高精度、无奇异性的计算近区积分。

步骤4、计算观测点重力场：通过对近区场积分和远区场积分求解计算得到观测点重力场值。

以下为本发明计算一个六面体重力场的实例。

如图6所示，设单一六面体的密度为2000kg/m³，长100m，宽100m，高100m。令六面体的上表面中点为笛卡尔坐标系的原点，向下方向为z轴的正方向。M＝210201个观测点分布在六面体上表面上方高度为h＝100m的平面上，观测平面为200m×200m。采用TetGen将六面体剖分成N＝20415个四面体单元。

为了验证本发明方法的正确性，进行一下两个验证：四面体重力场的解析表达式的正确性验证和采用快速多极展开的重力场计算的正确性验证。首先验证四面体重力场的解析表达式的正确性验证，将本发明的解析解计算的结果与六面体单元的重力解析解计算结果(Li,et al,1998)进行对比，只对比了重力位与垂直方向的重力梯度两个结果，如图7所示，两者的结果十分吻合，重力位的最大误差为6.0×10^-11％，垂直方向的重力梯度的最大误差为4.0×10^-9％。说明本发明推导的四面体重力场解析表达式的正确性。下面验证采用自适应快速多极展开算法的重力场计算的正确性，将结果与采用四面体解析表达式计算的结果作对比，同样只比较重力位与垂直方向的重力梯度两个分量，如图8所示，重力位的最大误差为1.0×10^-5％，垂直方向的重力梯度的最大误差为5.0×10^-9％，二者的结果十分吻合，说明本发明基于快速多极展开算法的重力场计算方法的正确性。使用基于四面体单元的解析表达式进行计算的方法计算时间为63.6s，使用本发明基于快速多极展开算法的重力场计算方法的计算时间为7.1s，加速比为8.95倍。并且，观测点数越多，加速比越大。同时采用OpenMP并行计算加快计算，当M＝1002001时，加速比达到8661倍。

参考文献

Li X,Chouteau M.Three-dimensional gravity modeling in all space[J].Surveys in Geophysics,1998,19(4):339-368.

Wilton D R,Glisson A W,Schaubert D H,et al.Potentials integrals for uniform and linear source distributions on polygonal and polyhedral domains[J].IEEE Transactions on Antennas and Propagation,1984,32(3):276-281。