一种球坐标系下重力场快速高精度正演方法与流程

本发明涉及地球物理技术领域，具体涉及一种大尺度球坐标系下重力场的快速高精度正演方法；该方法通过引入核矩阵等效性策略使得计算效率提高了两个数量级，同时大大降低了内存占用，该方法还通过引入自适应细剖分策略减少计算时的相对误差，使得正演精度大大提高。

背景技术：

重力勘探作为一项最基本的地球物理方法，已被广泛运用于地形校正、水文地质勘察、大地水准面测量、固体矿产资源勘探、岩石圈壳幔结构研究等技术领域。由于通过测量得到的重力场数据与地下密度异常体有着密切关系，因此该数据可被直接用于反演得到地下三维的密度分布，从而推断出地下构造运动、岩石圈结构构造等情况。

由于正演的速度和精度决定了反演的可行性及可信度，因此在一个高效的反演程序中，快速高精度正演方法占有极其重要的地位。对于小范围区域内的重力场勘探，可利用传统的直角坐标系并将地下场源剖分为规则的直立棱柱体，使用的快速高精度正演方法有fft算法、gauss-fft算法和快速多级展开法，对于这种情形下的重力场正演来说，计算效率和精度已得到保证。

然而，当研究至某一区域乃至全球尺度范围的重力场时，就不得不考虑到地球曲率的影响，相应地，所有的数据处理手段、正演和反演都将在球坐标系下进行。现有的球坐标系下三维重力场正演方法有二维高斯勒让德算法、三维高斯勒让德算法、泰勒级数展开法等；上述方法均存在正演效率低、计算精度差的缺点，严重制约着大规模重力场反演的广泛使用。

例如，传统的三维高斯勒让德算法非常耗时，尤其是当通过增加高斯节点来提高数值精度时，或当剖分的网格数和观测点数较多时，计算时间将呈指数增长，影响正演效率；另外，当观测点距离地下场源太近时，许多tesseroid单元体本身的大小与其中心点到观测点的距离不符合比例要求，导致正演精度很低。

在当前计算机硬件条件难以产生飞跃式突破的背景下，想要解决上述问题就只能寻求新的正演方法。

技术实现要素：

本发明的目的在于提供一种适用于大尺度球坐标系下且具有高精度和高效率的重力场正演方法，以解决背景技术中提出的问题。

为实现上述目的，本发明提供了一种球坐标系下重力场快速高精度正演方法，包括如下步骤：

步骤1)先将地下三维场源剖分成多个等间隔的规则tesseroid单元体，每个tesseroid单元体均由三对曲面围成，即一对同球心的曲面(r1＝const，r2＝const)、一对子午线断面(λ1＝const，λ2＝const)、一对同轴圆锥平面然后计算单个tesseroid单元体在场源体外一观测点p处的重力响应，其中，第q个tesseroid单元体在p点处产生的重力位v、重力加速度gα和重力梯度张量gαβ分别表示为如下形式：

其中α、β∈{x，y，z}，ρq是第q个tesseroid单元体的密度，g是重力常数，δαβ是克罗内克符号，如果α＝β，则δαβ＝1，如果α≠β，则δαβ＝0；并且

步骤2)在剖分tesseroid单元体时，将地下场源在径向、纬向、经向上分别剖分成段，剖分出的tesseroid单元体总数为在场源的正上方选取一个曲形的观测面，所述观测面在纬向和经向上分别剖分成和n'λ段且与地下场源的剖分一一对应，观测点位于观测面上且其位置与各tesseroid单元体的几何中心一一对应，观测点总数因此，公式(1)～(3)中的重力位v、重力加速度gα和重力梯度张量gαβ可以表示成矩阵相乘的形式，即：

kv·ρ＝v,(7)

kα·ρ＝gα,(8)

kαβ·ρ＝gαβ,(9)

其中，ρ为tesseroid单元体的密度向量且维度为ntess，v、gα和gαβ为正演结果向量且维度为nobs，kv、kα和kαβ为重力场正演核矩阵且维度为nobs×ntess；

步骤3)将公式(7)～(9)展开，可进一步表示为：

其中，p＝1、2……nobs和q＝1、2……ntess分别是观测点p(i，j)和场源点q(k，l，m)的索引，

步骤4)根据核矩阵等效性，大小相等的tesseroid单元体在等距离观测点处产生的重力响应相等或者相差一个负号，从公式(4)和公式(10)～(12)可以看出，kv、kx、kz、kxx、kxz、kyy和kzz分量均是(λ'-λ)的偶函数，以kv为例，tesseroid单元体对于观测点的重力场核矩阵等效性可简化为如下公式：

kv(i,j；k,l,m)＝kv(i,|j-l|+1；k,1,m),(13)

而ky、kxy和kyz分量均是(λ'-λ)的奇函数，以ky为例，tesseroid单元体对于观测点的重力场核矩阵等效性可简化为如下公式：

ky(i,j；k,l)＝signum×ky(i,|j-l|+1；k,1),(14)

其中，

步骤5)根据重力场的叠加性对多个tesseroid单元体的重力响应进行累加求和，从而得到地下所有物质在观测点p处产生的重力场，首先算出由q(k，l，m)产生的核矩阵元素，然后通过等效关系式(13)～(15)得到其它核矩阵元素，利用上述等效性方法使原来的四维核矩阵降低为三维核矩阵，且核矩阵的存储量也减小至原来的1/n'λ，以提高计算效率。

所述核矩阵等效性的意思是：如果两个(或多个)tesseroid单元体的大小相等，并且这两个(或多个)tesseroid单元体到观测点的距离相等，那么这两个(或多个)单元体在观测点所产生的重力响应也相等或者相差一个负号。

优选地，为了确保正演精度，需将原始剖分后的tesseroid单元体细剖分至合适大小，首先计算观测点到tesseroid单元体几何中心的距离：

其中cosψt＝sinφsinφt+cosφcosφtcos(λ-λt)，

然后分别计算tesseroid单元体在径向、纬向和经向上的长度：

lr＝r2-r1,(17)

lφ＝r2arccos(sinφ2sinφ1+cosφ2cosφ1),(18)

lλ＝r2arccos(sin²φt+cos²φtcos(λ2-λ1)).(19)

最后判断观测点到tesseroid单元体中心点的距离与该tesseroid单元体在径向、纬向和经向上各长度间的关系：

对于每一个d称之为该方向上的长度尺度比，如果不等式(20)对于三个方向上的长度尺度比均满足均满足，则该单元体tesseroid就不需要再细剖分，否则继续细剖分该单元体tesseroid。

在本发明中，将场源体剖分成等间隔的规则tesseroid单元体，核矩阵等效性是运用于原始剖分后的tesseroid单元体上而非细剖分后的tesseroid单元体，如果某一原始tesseroid单元体被细剖分为许多个子单元体，参与核矩阵等效性的仍然是被原始剖分的tesseroid单元体。从上面的公式中可以看出，由于核矩阵等效性，在计算时只需计算第一个经度圈上所有单元体q(k，1，m)在所有观测点处的重力响应，而其它单元体在所有观测点处的重力响应不用再计算了。因此，自适应自剖分也只需在第一个经度圈上所有的原始单元体上进行，如果第一个经度圈上某一个单元体被细剖分成若干个子单元体，那就将这些子单元体在某一观测点上的重力响应都计算出来，再累加求和，就是这个原始单元体在该点处的重力响应，经过细剖分之后的计算精度比直接用使用原始单元体的计算精度要高很多。

本发明提供的技术方案至少具有如下有益效果：

1、本发明通过深度分析地下tesseroid单元体的剖分规则以及核矩阵表达式，提出了核矩阵等效性策略，利用核矩阵元素间的等效关系使原来的四维核矩阵降低为三维核矩阵，减少计算量，同时核矩阵的存储量也减小至原来的1/n'λ，减少计算量和内存占用，达到大幅度提高计算效率的目的。

2、本发明通过引入自适应细剖分策略，对观测点到tesseroid单元体中心点的距离与该tesseroid单元体在径向、纬向和经向上各长度间的比例关系进行判断，确保每个tesseroid单元体被细剖分至合适大小，当将两种策略结合在一起时，仅需要在第一列(即第一个经度圈)的原始tesseroid单元体上运用自适应细剖分策略，而其他核矩阵元素可利用核矩阵等效性快速推出，在提高正演精度的同时也提高了计算速率。

附图说明

为了更清楚地说明本发明技术方案，下面将对本发明实施例描述中所需使用到的附图作简单地介绍，显而易见地，下列附图仅仅用于帮助理解本发明中的部分实施例而非技术方案的全部，其中：

图1是本发明实施例中提到的球坐标系下tesseroid单元体的示意图；

图2是地下三维场源和观测面的剖分示意图，和δλ分别表示在径向、纬向和经向上的剖分间隔；

图3是核矩阵平移等效性的原理示意图，以第一层模型为例，图中位置关系相互平行的两条/三条虚线所表示的核矩阵元素相等；

图4是核矩阵互换对称性的原理示意图，以第一层模型为例，图中位置关系为相交的两条虚线所表示的核矩阵元素相等。

图5是本发明方法与传统三维高斯勒让德算法在计算时间及内存占用方面的对比图，图(a)是重力异常gz分量相对计算时间随观测高度的变化对比图，图(b)是重力梯度gzz分量相对计算时间随观测高度的变化对比图，图(c)是重力梯度gzz分量相对计算时间随网格剖分间隔的变化对比图，图(d)是内存占用随网格剖分间隔的变化对比图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明的部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其它实施例，都属于本发明保护的范围。

对比实施例1(采用传统三维高斯勒让德算法)：

如图1所示，先将地下三维场源剖分成多个等间隔的规则tesseroid单元体，每个tesseroid单元体均由三对曲面围成，即一对同球心的曲面(r1＝const，r2＝const)、一对子午线断面(λ1＝const，λ2＝const)、一对同轴圆锥平面然后计算单个tesseroid单元体在场源体外一观测点p处的重力响应，其中，第q个tesseroid单元体在p点处产生的重力位v、重力加速度gα和重力梯度张量gαβ分别表示为如下形式：

其中α、β∈{x，y，z}，ρq是第q个tesseroid单元体的密度，g是重力常数，δαβ是克罗内克符号，如果α＝β，则δαβ＝1，如果α≠β，则δαβ＝0；并且

需注意的是，由于公式(3)中包含对λ'和的椭圆积分，因此他们没有解析公式，需要通过数值求解的方法去离散，公式(1)～(3)可以采用如下三维高斯勒让德算法求解：

其中ω是积分区域，i、j、k分别为径向、经向、纬向上的高斯节点数，以及

是变换到每个tesseroid单元体上的高斯节点，和分别是区间[-1，1]上第i、j、k个高斯系数和高斯节点。

实施例1(采用本发明方法)

步骤1)如图1所示，先将地下三维场源剖分成多个等间隔的规则tesseroid单元体，每个tesseroid单元体均由三对曲面围成，即一对同球心的曲面(r1＝const，r2＝const)、一对子午线断面(λ1＝const，λ2＝const)、一对同轴圆锥平面然后计算单个tesseroid单元体在场源体外一观测点p处的重力响应，其中，第q个tesseroid单元体在p点处产生的重力位v、重力加速度gα和重力梯度张量gαβ分别表示为如下形式：

其中α、β∈{x，y，z}，ρq是第q个tesseroid单元体的密度，g是重力常数，δαβ是克罗内克符号，如果α＝β，则δαβ＝1，如果α≠β，则δαβ＝0；并且

步骤2)如图2所示，在剖分tesseroid单元体时，将地下场源在径向、纬向、经向上分别剖分成段，剖分出的tesseroid单元体总数为在场源的正上方选取一个曲形的观测面，所述观测面在纬向和经向上分别剖分成和n'λ段且与地下场源的剖分一一对应，观测点位于观测面上且其位置与各tesseroid单元体的几何中心一一对应，观测点总数因此，公式(1)～(3)中的重力位v、重力加速度gα和重力梯度张量gαβ可以表示成矩阵相乘的形式，即：

kv·ρ＝v,(7)

kα·ρ＝gα,(8)

kαβ·ρ＝gαβ,(9)

其中，ρ为tesseroid单元体的密度向量且维度为ntess，v、gα和gαβ为正演结果向量且维度为nobs，kv、kα和kαβ为重力场正演核矩阵，维度为nobs×ntess；

步骤3)将公式(7)～(9)展开，可进一步表示为：

其中，p＝1、2……nobs和q＝1、2……ntess分别是观测点p(i，j)和场源点q(k，l，m)的索引，

步骤4)如图3和图4中所示的核矩阵等效性，大小相等的tesseroid单元体在等距离观测点处产生的重力响应相等或者相差一个负号，从公式(4)和公式(10)～(12)可以看出，kv、kx、kz、kxx、kxz、kyy和kzz分量均是(λ'-λ)的偶函数，以kv为例，tesseroid单元体对于观测点的重力场核矩阵等效性可简化为如下公式：

kv(i,j；k,l,m)＝kv(i,|j-l|+1；k,1,m),(13)

而ky、kxy和kyz分量均是(λ'-λ)的奇函数，以ky为例，tesseroid单元体对于观测点的重力场核矩阵等效性可简化为如下公式：

ky(i,j；k,l)＝signum×ky(i,|j-l|+1；k,1),(14)

其中，

步骤5)根据重力场的叠加性对多个tesseroid单元体的重力响应进行累加求和，从而得到地下所有物质在观测点p处产生的重力场，首先算出由q(k，l，m)产生的核矩阵元素，然后通过等效关系式(13)～(15)得到其它核矩阵元素，利用上述等效性方法使原来的四维核矩阵降低为三维核矩阵，且核矩阵的存储量也减小至原来的1/n'λ，以提高计算效率。

在本实施例中，为了确保正演精度，需将原始剖分后的tesseroid单元体细剖分至合适大小，首先计算观测点到tesseroid单元体几何中心的距离：

其中cosψt＝sinφsinφt+cosφcosφtcos(λ-λt)，

然后分别计算tesseroid单元体在径向、纬向和经向上的长度：

lr＝r2-r1,(17)

lφ＝r2arccos(sinφ2sinφ1+cosφ2cosφ1),(18)

lλ＝r2arccos(sin²φt+cos²φtcos(λ2-λ1)).(19)

最后判断观测点到tesseroid单元体中心点的距离与该tesseroid单元体在径向、纬向和经向上各长度间的关系：

对于每一个d称之为该方向上的长度尺度比，如果不等式(20)对于三个方向上的长度尺度比均满足均满足，则该单元体tesseroid就不需要再细剖分，否则继续细剖分该单元体tesseroid。

上述计算过程可通过计算机编程手段实现，用于正演模拟实例时，先在程序中设置好相关的正演参数，包括观测高度的确定、网格剖分的间隔和段数等，然后确定所需计算的重力场分量并给定长度尺度比、高斯节点和高斯系数，接着启用上述自适应细剖分策略和核矩阵等效性策略，运行程序后可快速且精准地得出结果。

一方面，根据自适应细剖分策略，当观测点距离场源较近时，原始的tesseroid单元体将会被进一步剖分成更小的单元体，而tesseroid单元体数量的增加将导致计算效率的降低，因此观测高度将不可避免对计算效率产生影响；另一方面，在球坐标系下的tesseroid单元体没有解析解，因此，为了对比本发明方法的高效性，需用一个均质球壳层作为标准模型，其在场源外一点的重力场可以通过解析公式求得。

分别采用对比实施例1和实施例1中的方法在全球尺度范围(卫星观测)内进行重力场正演模拟，球壳层厚度为100km，起始范围分别为6271到6371km(即地球实际半径)，球壳密度设置为1000kg/m³，将该球壳层在径向和纬向上均剖分成间隔为1°的tesseroid单元体(共计180×360个)，观测高度范围为地面上距离1km到250km。所得的数据结果见图5；其中，图(a)是重力异常gz分量相对计算时间随观测高度的变化对比图，图(b)是重力梯度gzz分量相对计算时间随观测高度的变化对比图，图(c)是重力梯度gzz分量相对计算时间随网格剖分间隔的变化对比图，图(d)是内存占用随网格剖分间隔的变化对比图。

首先测试计算效率随观测高度的变化。从图(a)和图(b)中可以看出，本发明方法在计算效率方面相比传统方法提高了近两个数量级，且离场源越近，提高越明显。此外，gz分量的计算时间几乎不随观测高度的变化而变化，而gzz分量计算时间随着观测高度的降低呈指数增长。

然后测试gzz分量相对计算时间及内存占用与网格剖分间隔的变化关系，在这个测试中，观测高度固定为50km不变，改变网格的剖分间隔(也即网格个数)且剖分间隔从0.25°到20°变化，测试全程在一台i7-4790kcpu、16gb的台式机上进行且无多核并行计算。所得关于gzz分量绝对计算时间和内存占用与网格剖分间隔间的关系见表1。

需要注意的是，所述内存占用指的是核矩阵的大小，即反演中雅克比矩阵的大小。传统方法在这里是不需要存储核矩阵的，但是在反演计算中，核矩阵(雅克比矩阵)的存储是必须的。

表1

从图(c)、图(d)和表1中可以看出，在本发明方法和传统方法中，相对计算时间都随网格剖分间隔的减小而增大，但存在明显区别的是，本发明方法相比传统方法在计算效率上提高了两个数量级，内存占用反而减小了两个数量级，且剖分间隔越小，网格数越多，本发明方法的优势越明显，尤其对于网格剖分间隔小于1°的实施例来说，本发明方法具有绝对优势，这种优势在重力场反演中也更加明显，大大提高了正演精度，保证最大相对误差低于0.1％。

以上所述仅为本发明的优选实施例，并非因此限制本发明的专利保护范围，对于本领域的技术人员来说，本发明可以有各种更改和变化。在本发明的精神和原则之内，凡是利用本发明说明书及附图内容所作的任何改进或等同替换，直接或间接运用在其它相关的技术领域，均应包括在本发明的专利保护范围内。