重力及其梯度异常正演的Moving-footprint大尺度模型分解方法

石泽玉，张志厚,2，刘鹏飞，范祥泰

(1.西南交通大学 地球科学与环境工程学院，四川 成都 611756; 2.西南交通大学 高速铁路线路工程教育部重点实验室，四川 成都 610031)

0 引言

重力及其梯度异常是由于地球局部质量分布不均匀而产生的[1]。重力及其梯度异常正演是根据已知地质体的形状、产状和剩余密度等来计算异常的分布规律，是重力勘探定量解释的基础，其对重力反演[2-4]及辅助导航[5]都具有非常重要的意义。大规模重力及其梯度异常的正演速度决定了反演的可行性[2]，也成为海量重力数据处理迫切需要解决的现实问题[6]。

为了提高重力及其梯度正演的计算效率，近年来众多学者提出了多种重力或其梯度的正演方法与技术，如Li等[7]应用小波变换和阈值小波系数组合压缩重力异常正演灵敏度矩阵来减少计算量，但仍需要对灵敏度矩阵进行计算；秦朋波等[8]发现了规则网络情况下核函数的对称性，并提出了一种快速计算灵敏度矩阵的计算方法；熊光楚[9]推导了长方体单元重力异常的傅里叶变换表达式；Shin等[10]提出了一种基于快速傅里叶变化的频率域方法，该方法通过改进已有算法加强了运算中的周期性，减少了运算中的数量级，在大型数据集进行计算时可以有效提升计算效率；Wu等[11]受偏移抽样技术的启发，提出了重力异常正演的高斯快速傅里叶变换(Gauss-FFT)计算方法，相比标准的快速傅里叶(FFT)正演方法，Gauss-FFT计算精度更高、速度更快；Ren等[12]引入了自适应快速多极子方法可对任意起伏地形的重力异常进行快速正演。在空间域，姚长利等[2]提出几何格架函数的方法，首先将场源划分成若干规则长方体单元，然后计算部分测点的格架函数并存储，其余测点格架函数可利用平移等效性和互换对称性直接调用，从而大大减少了计算量和存储量；陈召曦等[3]在几何格架函数方法的基础上提出了多核CPU加速并行计算的方法，以此提升正反演速度。张志厚等[4]提出了网格点几何格架函数的概念，该方法实质上是对几何格架函数方法进行加速。两者不同点在于几何格架函数的概念是相对于长方体单元，网格点几何格架函数的概念是相对于长方体的角点，其改进在于避免了网格点几何格架函数的多次重复计算。因此，当剖分单元足够大时，理论上网格点几何格架函数的计算效率能够提高近8倍。

但随着航空地球物理的发展，大面积海量重力数据面临着高精度快速处理的挑战。而以上方法只能对较少的数据量进行处理，如采用网格点几何格架函数的策略[4]对256×256×15的单元体进行正演至少需要3个多小时，当剖分单元扩大4倍到512×512×15时，其正演时间呈指数上涨，即在普通计算机上很难实现一次正演，难以满足实际生产的需求。

受航空电磁Moving-footprint大尺度模型分解正演技术的启发，本文提出了一种重力及其梯度异常正演的Moving-footprint大尺度模型分解计算方法。

文中将地下半空间规则剖分成若干长方体单元，某观测点重力及其梯度异常主要为其正下方一定范围内(子空间)的物性单元体产生，当观测点移动时，子空间跟随移动，即为“Moving-footprint”。文中划分了不同尺度的子空间进行正演计算，并将计算结果与理论结果进行对比，以此来验证方法的适用性。

1 重力及其梯度异常正演

重力及其梯度正演是将地下半空间剖分成若干个长方体单元(如图1所示)，然后计算每一个长方体单元在观测点观测到的异常值，再将每一个长方体的异常值求和，得到的结果即为地下半空间内的异常[13]。单个长方体单元在观测点产生的重力异常及重力梯度异常理论计算表达为[14]式(1)～(7)。

图1 地下半空间单元划分Fig.1 Underground half-space unit division

(1)

(2)

(3)

(4)

(5)

(6)

(7)

2 基于Moving-footprint技术的重力及其梯度异常正演方法

Moving-footprint即移动脚印技术，该技术广泛应用于航空电磁探测领域，实现了大规模航空电磁数据的快速正反演[15-16]。Yin等[16]将航空电磁系统的Moving-footprint定义为地下有限导电半空间中的某一子空间，认为该子空间整体电磁响应约为地下半空间总电磁响应的90%，当发射接收系统移动，该子空间随之移动，即为航空电磁系统的Moving-footprint大尺度模型分解技术。在重力异常正演计算中应用Moving-footprint技术，即只考虑在距离观测点一定的范围内的网格点(子空间内)在观测点产生的异常影响(如图2所示)，而超出选定的子空间范围的网格点由于距离太远，在观测点产生的异常影响十分微弱，在计算中可以忽略不计。文中方法与航空电磁领域的Moving-footprint不同之处是选取的子空间在深度上与地下半空间整体深度一致，而不是航空电磁的浅层部分深度。

注：红线包围区域为子空间note：the area enclosed by the red line is the subspace图2 地下网格体子空间划分示意Fig.2 Underground grid subspace division schematic diagram

本文所提Moving-footprint的重力及其梯度异常正演方法是在网格点几何格架函数技术的基础上进行了改进。因此，首先在全空间内选取子空间的大小；然后计算子空间内部分测点的网格点的几何格架函数，并存储以备调用；随后计算某观测点异常时，先判断该观测点所在的子空间，再调用网格点格架函数和该子空间单元体的剩余密度，求和后获得该观测点异常；最后，利用Moving-footprint完成观测区域所有点的重力及其梯度异常正演，主要计算步骤：

1)将地下半空间剖分为规则的长方体单元，并对单元体的剩余密度进行赋值；

2)确定子空间大小，计算子空间的网格点格架函数并进行存储以备调用；

3)根据观测点确定子空间的相对位置，利用平移等效性和对称互换性调用网格点几何格架函数，并利用式(1)～(7)完成该观测点重力及其梯度异常的正演计算；

4)观测点移动，子空间随之移动，利用步骤3完成移动观测点重力及其梯度异常的正演计算；最终完成所有观测点的正演计算。

基于Moving-footprint技术，只计算对观测点起主要贡献的长方体单元产生的异常，从而大大减少了总运算量，提高了计算效率。

3 模型实验

为了检验本文所提方法的计算效果，将地下半空间剖分为256×256×15个长方体单元，长方体单元大小为0.1 km×0.1 km×0.1 km。采用4个长方体组合模型进行检验，长方体组合模型的剩余密度都为0.5 g/cm3，长方体组合模型大小为4.0 km×4.0 km×0.5 km，其中心点坐标分别为(13.0 km, 13.0 km, 0.75 km)、(20.0 km, 6.0 km, 0.75 km)、(6.0 km, 6.0 km, 0.75 km)及(20.0 km, 20.0 km, 0.75 km)，模型示意如图3所示，其中，正演计算点距为0.1 km×0.1 km。

注：红线包围区域为计算子空间，蓝色立方体区域为异常体区域note：the area surrounded by the red line is the calculation subspace, and the blue cube is the gravity anomaly area图3 模型示意Fig.3 Model diagram

通过式(1)～(7)利用网格点几何格架函数方法[4]计算获得重力及其梯度异常。图4所示为选取256×256全空间计算所得准确的地下异常体正演结果，其与地下异常体一一对应。采用32×32、24×24及16×16的子空间分别计算，获得重力及其梯度异常的结果分别如图5～图7所示，各子空间的计算时间如表1所示，以及全空间与各子空间的计算时间比如图8所示。

由图4～图8可以看出：①随着子空间范围的缩小，运算时间随之缩短，大大提高了计算效率；②子空间计算范围缩小，计算精度下降。

为了评价本文所提方法的精度，将子空间为32×32的计算结果(图5)与理论值(图4)相减，其结果如图9所示。可以看出误差值基本上在零值附近。

为了定量评价误差的大小，文中统计了重力及其梯度异常的最大、最小值，以及计算结果与理论值误差的均值和均方差(如表2所示)。均值公式为：

图4 256×256全空间运行结果Fig.4 256×256 full space operation result

图5 32×32子空间运行结果Fig.5 32×32 subspace operation result

图6 24×24子空间运行结果Fig.6 24×24 subspace operation result

图7 16×16子空间运行结果Fig.7 16×16 subspace operation result

表1 256×256全空间不同子空间运行时间

图8 全/子空间运算时间比值Fig.8 Full/subspace operation time ratio

图9 256×256全空间与选取32×32子空间计算误差Fig.9 256×256 full space and selected calculated 32×32 subspace error

(8)

式中：a1,a2,…,an为计算结果与理论值的误差矩阵元素；n为矩阵所包含的元素的数量。

均方差的公式为：

(9)

由表2可得，重力异常计算值与理论值误差的均值与均方差分别为1.326 5 g.u.、0.717 5 g.u.，相比其理论最大值(72.127 g.u.)、最小值(0.089 6 g.u.)的范围，误差相对较小。

为了进一步定量衡量计算结果的精度，文中同时也统计了子空间32×32计算结果的均方差和平均相对误差(如表3所示)，计算结果的均方差公式为：

(10)

式中：xi为子空间各观测点计算结果；(x*)i为全空间各观测点计算结果；n为观测点数量。理论值与其平均值的均方偏差为：

表2 全空间重力及其梯度异常最大、最小值及计算值与理论值误差的均值和均方差

(11)

(12)

表3 32×32子空间计算结果的均方差和平均相对误差

由表3可得：32×32子空间计算所得结果中重力异常和Uxz、Uyz两个梯度异常所得结果的精度较高。重力异常的平均相对误差为9.01%，Uxz、Uyz的平均相对误差数值在5%以下，另Uzz、Uxx、Uyy的平均相对误差数值在20%以下，Uxy的平均相对误差数值为29.09%。通过表3，并结合表2中理论最大值和最小值，可以看出，重力及部分重力梯度异常计算结果误差较小。

4 反演结果对比

有限内存拟牛顿方法已被证明在解决重力及其梯度正演中具有一定的优势性[8]。因此应用有限内存拟牛顿方法对Moving-footprint方法的计算效果进行检验。为验证本文所提正演方法在计算中的优势，采用文中方法与文献[4]所提正演方法在反演中的效果进行对比。首先选择较小空间的数据集进行运算。全空间网格剖分为16×16×9，子空间为8×8×9。模型为一长方体模型，大小为500 m×500 m×500 m，顶面埋深为500 m，剩余密度为1 g/cm3。分别应用现有计算方法与本文所提Moving-footprint方法得出的计算结果如图10所示。

正演方法和本文方法这两种计算方式所需的时间分别为：29.222 s、26.602 s。由此可知应用Moving-footprint方法可以降低运算时间。当选择32×32×9的全空间进行计算时，所得结果如图11所示，现有方法计算所需时间为1 543.53 s，而在应用Moving-footprint方法后，计算时间为703.25 s，因此应用Moving-footprint技术对计算效率进行了提升，且随着模型空间剖分数量的增加，反演的计算效率有显著提升。

a—原始方法；b—Moving-footprint方法a—the result obtained by the original method；b—the result obtained by applying the Moving-footprint method图10 全空间16×16×9不同方法反演结果Fig.10 Inversion results of different methods in full space 16×16×9

当选择更大模型剖分空间进行反演时，如256×256×9，嵌套已有的正演算法无法完成迭代过程。而嵌套本文所提Moving-footprint技术可以有效完成对大尺度模型的反演计算，迭代1次所需时间约为30 min。

a—原始方法；b—Moving-footprint方法a—the result obtained by the original method；b—the result obtained by applying the Moving-footprint method图11 全空间32×32×9不同方法反演结果Fig.11 Inversion results of different methods in full space 32×32×9

5 结论与建议

本文借鉴航空电磁正演计算的“Moving-footprint”技术，提出了基于“Moving-footprint”重力及其梯度异常的正演计算方法。该方法选择全空间范围内的一定子空间，只计算存储在子空间范围内的网格点的格架函数，即只考虑子空间范围内的计算点在观测点产生的重力及其梯度异常。在保证一定的计算精度的同时，缩短了运算时间，提升了计算效率。

海量重力数据线性迭代反演时，初步迭代应用“Moving-footprint”技术选取合适的子空间大小，提升了计算效率。从而使得大范围全空间重力及其梯度异常正反演计算成为可能。