Reconstruction of Discrete Gravity Data Using Three-Tier Equivalent Sources with Variable Sizes
-
摘要: 离散重力测量数据的重构是重力数据分析和处理中最重要的方法之一.等效源方法作为一种有效且稳定的重构方法，受到国内外学者的关注.基于位场理论，提出了三层变尺度等效源方法的技术方案，可有效恢复测量数据中长波长范围的信息，实现对离散重力数据的高精度重构，进而保障了在重力场重构基础上的重力梯度计算.详细地讨论了三层变尺度等效源的设置以及优化反演计算的过程，通过二维和三维理论模型试验，展示了方法的有效性及相对单层等效源方法的优势，并通过广西某地区实际重力测量数据的应用，验证了方法的抗噪能力和实用性.Abstract: Reconstruction of discrete gravity data is an important technique in data analysis and processing. As an effective and stable approach of reconstruction, Equivalent-source method has still been the interest of geoscientists. In this paper, the frame of multi-size equivalent-sources distributing in three layers is proposed based on the theory of potential field. This approach can effectively recover long-wavelength information in observed data and is able to realize the high-precision reconstruction, so that gravity gradients can be calculated accurately on the basis of the achieved equivalent-sources. Issues on settings of three-layer equivalent sources with variable sizes and the optimization of inversion algorithm are discussed in the work. The experiments of 2D and 3D models show the effectiveness and the advantage of the method relative to the reconstructed results using single-layer equivalent-source method. The real data example from a region of Guangxi, China, is applied and its result illustrates that the method has good noise suppression and practicability.
-
Key words:
- gravity /
- gravity gradient /
- reconstruction /
- three-tier equivalent sources /
- variable size
-
0. 引言
构建重力场模型是重力学研究的一项重要工作.利用地面、航空或海洋重力测量数据构建局部或区域重力场模型，可以为研究重力场空间分布特征、重力数据转换与分析及处理以及重力异常成因解译提供基础；另一方面，通过对离散重力观测数据的重构，可以得到描述重力场的连续函数，进而获得重力场在空间上连续的表达.
利用离散重力测量数据重构重力场模型的方法技术经历了几十年的发展，涌现出许多成功的方法(如马国庆等，2016)，其中基于等效源原理的方法是最常用的方法之一.等效源方法是利用位场的等效原理，通过在地下设置适当的虚拟场源(即“等效源”)，使其产生的场与观测数据达到充分一致，并以此代替实际场源来计算场源以外空间的场值及其变化率.Dampney(1969)最早应用等效源方法对重力场进行空间延拓，利用观测面下规则分布的单层质点构成等效场源，将复杂的空间延拓计算转换为简单的等效源正演.随后的国内外研究主要围绕两个方面的问题，即如何有效地设置与构造等效源和如何快速有效地求解等效源.
在构造等效源方面，许多学者将等效源构造成单层、偶层或多层，并采用不同的方案设置深度.例如：Needham(1970)在球谐域应用点质量等效源层；Bhattacharyya and Chan(1977)在曲面磁测数据处理中应用磁偶层等效源进行磁异常重构；Süenkel(1981)在Needham(1970)的方法基础上，讨论了单层与多层点质量层的球谐频谱特征并给出了快速计算等效源质量的方法；国内学者吴晓平(1984)最早利用不同大小区域平均重力数据建立多层点质量模型以求解物理大地测量边值问题；在大地测量领域，国内外学者对点质量模型进行了拓展和深化，如Vermeer(1991)、赵东明等(2001)、Klees et al.(2008)；Xia and Sprowl(1991)给出了依据1/2点距数据(0.5, 1.5, etc.)光滑度的单层等效源深度选取准则；Hansen and Miyazaki(1984)、王万银等(1991)、王万银和潘作枢(1996)等学者基于偶层等效源位实现曲线和曲面上的位场重构；徐世浙等(2004)联合边界单元和单层等效源进行重构且实现了航磁数据的向下延拓.在地球物理勘探领域，已有学者开始应用多层等效源，例如：庞旭林等(2010, 2011)和庞旭林(2014)在硕士论文和他的两篇会议论文中介绍了对航磁数据作航测高度调平(或“曲化平”)时利用3层等效源模拟不同深度场源的方法；魏华平(2014)和黎莎等(2015)则是应用多层等效源来简化重力三维物性反演.这些研究为优化等效源的设置、提高重构数据与测量数据的拟合精度提供了很好的借鉴.
无论以何种方式设置等效源，首先需要求取等效源的物性参数(如密度或视密度)，然而在实际应用中，这是一项基于海量数据的最优化问题，需要采用稳定且快捷的方法.前人的研究已提供了许多有效的途径.例如：Leão and Silva(1989)采用滑动窗口分块反演等效源以降低反演核矩阵维数；Mendonca and Silva(1994)采用数据稀疏压缩有效控制核矩阵维数；Li and Oldenburg(2000, 2010)、刘天佑等(2007)利用小波压缩技术压缩降低核矩阵维数；Oliveira and Barbosa(2013)利用“多项式等效源”取代等效源物性来减少参与反演计算的未知量等.此外，Li and Oldenburg(2001, 2010)、Barnes and Lumley(2011)在反演等效源的过程中引入了正则化用于压制数据包含的噪声.
利用等效源方法重构重力场的最终目标是使虚拟场源产生的重力场无限接近真实场.显然，若设置的虚拟场源与真实场源在空间上的分布十分接近，将会得到很好的重构效果，但是实现起来无异于地下空间的三维反演(如孙石达等，2016)，会带来诸多困难；等效源方法与三维密度反演的区别在于其目的是恢复空间场值，而不考虑对真实场源的准确重建，通过构建合理的等效场源模型，可以大幅减少计算困难，提高计算效率.重力场信号源于地表至地下深处均有密度分布，其波长范围广泛.根据离散等效源方法原理，等效源单元体(点、面或体)的密度均是常参数，所构造的等效源层通常被设置在近地表.如此，单层等效源难以保证所重构的重力场具有完整的波谱特性，随着高度(与等效源的距离)的增加，重构的重力场与真实场的差异也将增大，而设置多层等效源可有效弥补这方面的不足.此外，不同的等效源单元体(点、面、长方体)以及等效源设置深度的选取均会对重构效果产生影响.针对上述问题，本文提出了利用三层等效源重构重力场的技术方案，采用不同尺度的长方体作为等效源单元，并且根据重力场的波谱分析确定等效源的设置深度；通过构建三层变尺度等效源模型可以更有效地恢复数据中长波长信息，实现对离散重力数据高精度地重构，具有计算效率高、适用性强、重构精度高的优势.理论模型试验表明，在获得了准确重构重力场的基础上，还可得到高精度的重力梯度数据.
1. 方法原理
利用等效源方法重构重力场有3个主要技术环节，即等效源的设置、等效源物性反演和重构计算.
1.1 等效源设置方案
根据重力场的复杂程度，等效源在形式上可设置为单层或多层.理论上等效源层数越多构建的场源模型越接近于真实情况，重构的重力场精度也越高，但是必然使计算量成倍增长，而且等效源方法的精髓也在于等效性以及计算高效性.因此，为减小计算量及提高等效源模型的合理性(浅、中、深3个深度范围均有等效源分布)使其充分涵盖深部引起的所有频段成分，根据模型试验，建议设置三层等效源的方案，且各层等效源单元体几何尺度也随着深度不同而变化，其技术要点如下：
(1) 等效源单元体的选择.许多学者(如谢汝宽等，2015)认为若等效源深度方向延伸太小，将影响深部信息的重构，且在利用重构的重力场模型计算向上延拓时会产生较大误差.图 1展示了同等质量3种常用等效源单元体的重力效应随高度的衰减情况.图中单元体点质量(point mass)、面元质量(plate element)和长方柱体(rectangular prism)具有相同质量和中心深度，在与场源较近的高度范围内，长方柱体随高度增大衰减显然最慢.选择长方柱体作为等效源单元体有助于“吸收”更多的深部场源信息，故建议等效源单元体均采用长方柱体.另一方面，随着等效源设置深度的增加，等效源单元体的体积(尺度)也应加大，以保持等效源应有的敏感度，同时也可以减少等效源单元体数量，提高计算效率.
图 1 3种等效源单元体的重力效应衰减特征Fig. 1. Gravity attenuation characteristics of three equivalent sources(2) 第1层等效源参数选择.与前人设置单层等效源的方法类似，单元体水平尺寸可以依据测量数据的平均间距而定(MacLennan and Li, 2013).对于该层等效源深度，可根据Dampney(1969)及Xia and Sprowl(1991)的建议，将第1层等效源深度设置为1~3倍重力数据平均点距，这样既能较好地“吸收”观测重力数据中的短波长信息，又不易产生虚假高频噪声.第一层等效源可以水平展布，也可随观测面起伏(与观测面相似的等深面)展布，主要用于拟合重力场的短波长成分.
(3) 第2层等效源参数选择.设计该层等效源的目的是充分拟合重力场的中、长波长成分，弥补第1层等效源对重力场深源信息“吸收”的不足.前人的研究表明，近地表的单层等效源对重力场低频成分的表征是不足的，需要设置不同深度的等效源层进行组合(如吴晓平，1984).等效源的最佳深度与真实场源相同，但实际场源分布在不同深度上.因此，该层等效源深度设置是本方案的关键之一.估计场源深度的方法有很多，可以利用重力场波谱分析[如应用Syberg(1972)的匹配滤波方法]获得深部场源的平均深度，作为该层等效源深度选择的依据.当然，当观测面起伏较大时应用傅里叶变换计算的功率谱会有偏差，但依据其估计的场源平均深度可作为参考.根据模型试验，第2层等效源单元体几何尺寸可选择第1层的1.5~3.0倍，该层等效源可以沿水平面展布，也可随观测面起伏展布.
(4) 第3层等效源参数选择.加入第3层等效源旨在拟合重力场中“超长”波长成分，并压制边界效应.由于任何一个区域内的重力场都会包含诸如莫霍面起伏等深部因素引起的重力变化，即所谓的背景场.在应用等效源进行重构高空重力场时，这些成分会导致在边界附近产生边界效应(剧烈震荡)，背景场在边界上差异愈大，边界效应愈严重.此外，第1和第2层等效源也可能“泄漏”一些超长波长信息.可将第3层等效源置于第2层等效源之下，沿水平面展布.第3层等效源深度选取可参考区域内背景场“超长”波长特征，视具体情况而定.在没有其他资料参考的情况下，第3层等效源深度可选择3~10倍第2层的深度，第3层等效源单元体几何尺寸应比第2层更大(可选3~10倍)，同时等效源分布范围需大于重力数据的区域，具体可依据原始数据情况.根据模型试验，可选3~10倍于数据区域.
1.2 等效源物性反演算法
在等效源形状和空间分布确定之后，其重力效应与等效源物性(密度)之间为线性关系，利用矩阵方程表示可以写为：
$$ \boldsymbol{Gm} =\boldsymbol{d}, $$ (1) 式中：d为重力数据向量；m为等效源物性；G为正演核函数矩阵[公式(2)].本文选择的长方柱体作为等效源单元体，单个长方柱体重力效应正演计算采用Nagy(1966)的解析式，即：
$$ \begin{array}{l} d = Gm = \\ f\cdot m|\;\;|\;\;|A\cdot{\rm{ln}}\left({B + r} \right) + B\cdot{\rm{ln}}\left({A + r} \right) + C\cdot\\ {\rm{arctan}}\left({\frac{{Cr}}{{AB}}} \right)|_{{\xi _1}}^{{\xi _2}}\;|_{{\eta _1}}^{{\eta _2}}\;|_{{\zeta _1}}^{{\zeta _2}}, \\ A = \left({\xi - x} \right), B = \left({\eta - y} \right), C = \left({\zeta - z} \right), r = \\ \sqrt {{A^2} + {B^2} + {C^2}}, \end{array} $$ (2) 式中：d为柱体在(x, y, z)处重力效应值；G为核函数值；m表示柱体密度；f为万有引力常数；(ξ1, η1, ζ1)和(ξ2, η2, ζ2)分别表示柱体上、下角点坐标.
假设重力观测数据总数为M，第1~3层等效源的单元体总数分别为N1、N2和N3，等效源单元体总数为N，则式(1)可写成
$$ {\boldsymbol{G}_{M \times N}}{\boldsymbol{m}_N} = {\boldsymbol{d}_M}, $$ (3) 其中：
$$ \begin{array}{l} {\boldsymbol{m}_N} = [{m_1}, {m_2}, \cdots, {m_{{N_1}}}, {m_{{N_1} + 1}}, {m_{{N_1} + 2}}, \cdots, \\ {m_{{N_1} + {N_2}}}, {m_{{N_1} + {N_2} + 1}}, \cdots {m_{{N_1} + {N_2} + {N_3}}}{]^{\rm{T}}}. \end{array} $$ (4) 求解方程(3)即可得到等效源各个单元体的密度值.根据本方法等效源设置规则，观测数据数量一般会小于等效源单元体数量(M＜N)，因此方程欠定，需将方程(3)变换为：
$$ {\boldsymbol{G}^{\rm{T}}}\boldsymbol{Gm} = {\boldsymbol{G}^{\rm{T}}}\boldsymbol{d}, $$ (5) 进而采用共轭梯度法求解.
由于本方法中设置的等效源分布于不同深度，为避免场源分布过度集中于近地表，故在算法中引入“预优”技术.预优共轭梯度法亦称预优条件共扼梯度法，是求解线性方程组的有效方法之一.它是基于共扼梯度法的收敛速度且依赖于系数矩阵的特征值分布这一性态，事先对矩阵进行预处理，使其特征值分布较为集中，进而提高其收敛速度.预优共轭梯度法常被应用到地球物理反演中，如Pilkington(1997)、刘双等(2013)等.
将方程(5)两边同乘以预优矩阵P，使P(GTG)≈I，以改善矩阵(GTG)特征值分布.这里借鉴Pilkington(1997)的方法，取预优矩阵为：
$$ \boldsymbol{P} = {z^\beta }\boldsymbol{I}, $$ (6) 式中：z表示单元体中心到观测面的垂直距离(或平均垂直距离)；β为衰减系数，与场的衰减速率有关.加入预优矩阵，方程(5)变成：
$$ \boldsymbol{P}{\boldsymbol{G}^{\rm{T}}}\boldsymbol{Gm} = \boldsymbol{P}{\boldsymbol{G}^{\rm{T}}}\boldsymbol{d}. $$ (7) 关于预优因子zβ中β的取值，可参考Li and Oldenburg(1996, 1998)的“深度加权因子”中指数的选取方法.根据模型试验，β的取值范围在3~6之间为宜.事实上，这里的预优因子zβ还起到了深度加权的作用，可使反演得到的等效源更好地分布在不同深度上.
利用预优共轭梯度法迭代反演基本步骤如下：
(1) 首先设置反演精度(理论数据与实测数据的拟合差)，并令初始模型{mj}(0) (j=1, 2, …, N)=0；将观测数据dobs代入方程(7)进行反演，得到第一代密度模型m(1)，由此计算出测点处的理论数据并与实测数据比较，得到Δd(1)；
(2) 将Δd(1)代入方程(7)求解出下一代密度模型修改量Δm(0)，得到新一代模型m(1)=m(0)+Δm(0)；
(3) 如此重复，当经过第k次迭代后得到密度模型m(k)=m(k-1)+Δm(k-1)，能使拟合差Δd(k)满足精度要求，即终止迭代，m(k)便作为方程(7)的解minv.
1.3 重构重力场及重力梯度计算
将等效源密度向量minv各元素{mj}inv(j=1, 2, …, N)代入式(2)并累加求和，可得：
$$ d_i^{{\rm{rec}}} = \sum\limits_{j = 1}^N {{G_{ij}}m_j^{{\rm{inv}}}}, $$ (8) 式中：direc表示第i点处重构重力场值；Gij为第j柱体在第i点处的几何函数值.由式(8)可计算出空间各处的重力场值，即可获得重构重力场.
根据郭志宏等(2004)给出的长方体三坐标方向的重力梯度解析公式：
$$ {T_{xz}} = {G^{xz}}m = f\cdot m|\;\;|\;\;|\;{\rm{ln}}\left[ {\left({\eta - y} \right) + r} \right]{|_{{\xi _1}}}^{{\xi _2}}{|_{{\eta _1}}}^{{\eta _2}}{|_{{\zeta _1}}}^{{\zeta _2}}, $$ (9) $$ {T_{yz}} = {G^{yz}}m = f\cdot m|\;\;|\;\;|{\rm{ln}}\left[ {\left({\xi - x} \right) + r} \right]{|_{{\xi _1}}}^{{\xi _2}}{|_{{\eta _1}}}^{{\eta _2}}{|_{{\zeta _1}}}^{{\zeta _2}}, $$ (10) $$ \begin{array}{l} {T_{zz}} = - {G^{zz}}m = - f\cdot m|\;\;|\;\;|{\rm{arctan}}\\ \frac{{\left({\xi - x} \right)\left({\eta - y} \right)}}{{\left({\zeta - z} \right)r}}{|_{{\xi _1}}}^{{\xi _2}}{|_{{\eta _1}}}^{{\eta _2}}{|_{{\zeta _1}}}^{{\zeta _2}}, \end{array} $$ (11) 式中：Txz、Tyz、Tzz分别代表长方体三坐标方向重力梯度，Gxz、Gyz与Gzz是对应的几何函数值，m为柱体密度.将反演得到的等效源单元体密度{mj}inv(j=1, 2, …, N)分别代入式(9)、(10)和(11)，并累加求和可得：
$$ \left. \begin{array}{l} {({T_{xz}})_i}^{{\rm{rec}}} = \sum\limits_{j = 1}^N {} {G_{ij}}^{xz}{m_j}^{{\rm{inv}}}\\ {({T_{yz}})_i}^{{\rm{rec}}} = \sum\limits_{j = 1}^N {{G_{ij}}^{yz}{m_j}^{{\rm{inv}}}} \\ {({T_{zz}})_i}^{{\rm{rec}}} = \sum\limits_{j = 1}^N {} {G_{ij}}^{zz}{m_j}^{{\rm{inv}}}, \end{array} \right\}, $$ (12) 式中：(Txz)irec、(Tyz)irec与(Tzz)irec分别表示第i点处重构的三坐标方向重力梯度值，mj表示第j个等效源单元体密度，Gijxz、Gijyz与Gijzz为相应的几何函数值，由此可重构三方向重力梯度。此外，其他3个梯度(Txx)irec、(Tyy)irec与(Txy)irec重构，也可以代入相关的公式计算.
2. 理论模型试验
2.1 二维模型试验
为了验证方法效果，采用单元体为水平物质线和二维水平物质面两种单层等效源方法与本方法进行对比试验.如图 2a所示，本试验设计了包含4个不同尺度和深度二维密度异常体的组合模型，图 2b给出了它们在7 500 m长的观测剖面上的理论重力值(点距50 m)，其中模型密度均为(1 000 kg/m3).在试验中，单层等效源深度均为100 m；多层等效源深度分别为：第1层顶深50 m，单元体长宽为100 m×50 m；第2层顶深850 m(基于对数功率谱的估算)，单元体长宽为300 m×50 m；第3层顶深2 700 m，单元体长宽为500 m×50 m，分布于范围22.5 km，即分别向剖面两端延伸7.5 km.
图 2 二维组合模型及其产生的重力效应a.密度模型及多层等效源设置示意图，其中斜杠填充区域表示等效源分布区间；b.各模型理论重力异常及其总和Fig. 2. 2D synthetic example and theoretical gravity anomaly从图 3(左列)可见，本文的方法与单层等效源方法在观测面上重构的重力场(g)均与理论值拟合较好，但单层等效源重构得到的垂向重力梯度(Tzz)有显著误差.在利用重构方法计算向上延拓时，单层等效源方法的效果明显较差，无论重构重力场还是重力梯度，最大误差均接近或超过10%，且误差主要分布在边界附近，如图 3(右列)所示.可见，相对单层等效源方法而言，本文提出的重构方法在恢复重力场“全波段”信息和抑制边界效应方面都较单层等效源方法有明显改善.
图 3 观测剖面及750 m高度上重构结果及其与理论值的拟合差分布图中左列为各方法在观测面上重构重力和重力梯度及其与理论值的拟合差；右列为各方法在上延至750 m高度的重力和重力梯度及其与理论值的拟合差Fig. 3. Reconstruction results and their difference relative to theoretical values on observation surface and on the altitude of 750 m2.2 三维模型试验
为检验本方法在复杂条件下应用效果，试验中设置了起伏观测面、背景场和噪声(误差).三维模型由一组不同形状的模型体组成，如表 1及图 4a所示；观测面起伏的高程范围为0~411.6 m，如图 4b；理论总重力异常中加入了3%(均方差δ=±0.07 mGal)的高斯白噪声，如图 4e所示.
图 4 三维模型组合与观测面以及理论重力场a.三维模型示意图；b.起伏的观测面；c.观测面上的理论重力场；d.观测面上背景场；e.含噪声的模型重力场与背景场之和；f.对数功率谱曲线Fig. 4. 3D synthetic example, observation surface and theoretical gravity anomaly表 1 模型参数设置Table Supplementary Table Designed parameters of models模型编号 模型形状 模型中心距高程h=0
处平面深度(km)密度
(103 kg/m3)1 小方块 0.140 1.5 2 小方块 0.140 1.5 3 小方块 0.115 1.5 4 小方块 0.050 1.5 5 小方块 0.100 1.5 6 小方块 0.170 1.5 7 长方柱体 0.300 1.0 8 长方柱体 0.550 1.0 9 板状体 0.650 2.0 10 倾斜板状体 0.500 1.0 等效源设置方案：第1层等效源单元体水平尺度为100 m×100 m，与测点网间距一致，厚度50 m，其顶面距地面250 m，随地形起伏展布；第2层等效源单元体水平尺度为150 m×150 m，厚度200 m，根据由观测数据径向对数功率谱估算结果(图 4f)，设置其顶面在观测面之下380 m，随地形起伏展布；第3层等效源沿水平面展布，其单元体水平尺寸为第一层5倍，即500 m×500 m，厚度1 500 m，其顶面深度为海平面(h=0 m)之下1 600 m处，且分布范围5倍于数据区域.
依照上述等效源设置方案，经反演与重构得到的起伏观测面上重力场(图 5b)与理论值(图 5a)十分吻合，平均绝对误差和均方差分别为±0.033 mGal和±0.025 mGal，且噪声得到了有效压制.利用等效源计算向上延拓至h=1 000 m平面上所获得重力场与理论值也十分接近，平均绝对误差和均方差分别为±0.57 mGal和±0.46 mGal，重构结果中边界效应得到有效抑制.
图 5 三维模型及背景理论重力场及其重构结果a.无噪声理论模型与背景场；b.用图 4e重构的观测面重力场；c.h=1 000 m平面上的理论模型与背景场；d.用图 4e重构的h=1 000 m平面上结果；e.观测面上重构重力值与理论值拟合差；f.h=1 000 m平面上重构重力值与理论值拟合差Fig. 5. The gravity field from non-noise 3D models and background and reconstructive results利用离散重力观测数据计算重力梯度一般采用空间差分法或利用傅里叶变换进行转换.然而，对于在起伏观测面上测得的数据，转换的结果会产生畸变.高精度重力场重构结果为计算重力梯度提供了保障.笔者根据带噪声的三维模型及背景场数据(图 4e)进行了重力梯度重构试验，试验中未对数据中噪声施加任何处理，以此检验方法的抗干扰能力.图 6展示了起伏观测面上观测结果和“曲化平”延拓至420 m高度平面上两组结果.“曲化平”是将曲面数据延拓至地形最高点平面上的转换，由此能获得平面上场值或梯度，为后续数据分析和处理提供条件.
图 6 三维模型重构计算(Rec)的重力梯度结果与理论值比较左侧(1, 2)两列为起伏观测面上重构后计算(Rec)的重力梯度及其理论值；右侧(3, 4)两列为420 m高度上重构计算的(Rec)重力梯度及其理论值；a、b与c分别为Txz(北向水平梯度)，Tyz(东向水平梯度)和Tzz(垂向梯度)；图例单位为：重力值(mGal)Fig. 6. Comparison between the reconstruction gravity gradients and theoretical values of 3D models and background fields从观测面上重构的重力梯度(图 6a2, 6b2, 6c2)来看，Txz(北向水平梯度)、Tyz(东向水平梯度)和Tzz(垂向梯度)3个方向的重力梯度与无噪声的理论值(图 6a1, 6b1, 6c1)吻合较好，重构重力梯度较全面地反映了理论场梯度值的细节特征，但局部也存在小幅波动.Txz和Tyz与无噪声理论值之间均方差均小于±2 E，Tzz与无噪声理论值之间误差稍大，但仅为±8 E.分析表明，这些局部波动主要源于噪声干扰，但总体上对干扰的压制作用是明显的.图 6a4, 6b4, 6c4展示了重构重力梯度Txz、Tyz和Tzz在420 m高度平面上的分布.由于上延会使干扰影响减弱，从图上看局部波动有一定程度减小，重构的重力梯度与理论值更为接近，但中低频波动有所增大，但总体上与观测面上重构结果误差水平相当(误差统计见图 7).
图 7 三维模型重力梯度重构结果与理论值拟合差a、b分别表示观测面上与420 m高度上拟合差，(1~3)分别表Txz(北向水平梯度)、Tyz(东向水平梯度)和Tzz(垂向梯度)Fig. 7. Difference between the reconstruction gravity gradients and theoretical values of 3D models and background fields3. 实际数据应用
本文对广西某地区实测重力数据(已做必要处理)进行实际数据试验，重力数据采用自由网观测，平均测点间距750 m.研究区地形起伏较大，区内山谷沟壑纵横(见图 8a)，最高海拔898 m，最低15 m.研究区内重力异常总体表现为南西高、北东低，其中一条重力低异常带沿北东向展布，其中包含了多个重力低异常圈闭，如图 8b所示.区内局部重力异常主要反映了古生界沉积岩顶面起伏及其隐伏岩体的分布特征.根据在该地区开展的地质、地球物理工作所得到的认识，3层等效源设置参数如下：第1层和第2层等效源均设置为随观测面起伏，第1层等效源位于观测面下2.5倍点距(2.5×750 m)深度上，水平尺寸750 m×750 m，厚度500 m；第2层等效源顶面深度为观测面下3 000 m(图 8e)，水平尺寸取第1层的1.5倍，即1 125 m×1 125 m，厚度5 000 m；第3层等效源位于海平面以下9 000 m深处，水平尺寸10 km×10 km，厚度20 000 m.考虑到研究区背景场特征，对第1层等效源外延5倍点距(5×750 m)，第2层等效源外延7.5倍点距(7.5×750 m)，第3层等效源分布面积约为实际数据面积20倍.
图 8 研究区地形、实测重力异常与重构结果a.研究区地形等高线；b.实测重力异常；c.重构起伏地面上重力异常；d.重构海拔910 m平面上重力异常；e.对数功率谱曲线Fig. 8. Topography, observed gravity anomaly and reconstructive results in study area重构的观测面重力异常(图 8c)很好地还原了实测异常的区域和局部特征，并使干扰得到了有效压制，通过对重构数据拟合差的统计，测区绝大部分区域误差未超过1.5%.“曲化平”(延拓至海拔910 m)得到的平面重力异常(图 8d)较为真实地反映了基岩隆起及隐伏花岗岩体的分布；其整体受北东向、北西向及南北向地质构造的影响，多个显著重力异常低与出露花岗岩相吻合(冯建辉，2016).重构的起伏地面重力梯度(图 9a)及海拔910 m高处重力梯度(图 9b)反映更丰富的细节信息，受干扰影响较小，较好地反映了研究区的地质构造和地下岩性分布，尤其是研究区西部次级重力低圈闭的垂向重力梯度反映了更多构造信息，其中一些重力梯度异常成因已被其他地球物理方法所证实(陈家联等，2010; 朱国器等，2011).
图 9 重构的重力梯度a1、a2和a3分别为观测面上Txz(北向水平梯度)，Tyz(东向水平梯度)和Tzz(垂向梯度)；b1、b2和b3分别为上延至海拔910 m平面上的Txz, Tyz和TzzFig. 9. Reconstructive gravity gradients实际数据试验证明了本文提出的方法具有较好抗噪能力，重构所恢复的重力场信息可靠且波谱完整.
4. 结论
重构重力场的核心目标是尽可能地保全所有原始观测数据的有效信息.本文利用多层、变尺度等效源技术，尝试着在重构重力场过程中保全更多的信息，以提高离散重力数据重构精度，进而计算出高保真的重力梯度.通过理论模型试验和实际数据应用，得到以下几点认识：
(1) 三层、变尺度等效源在获取原始观测数据中全频段信息方面要优于单层等效源.本文采用的浅、中、深3层等效源技术方案，首先对等效源单元进行了优化选择，其次基于深部场源平均深度的估计确定第2层等效源深度，最后随深度逐级增大等效源单元体积，提高深部等效源敏感度.通过模型实验及实际数据的应用，证明本文提出的方案包含的信息频段更加完整，能减少重构中的频段“泄露”和边界效应，有效提高了计算精度，在重力场数据处理与解释中具有广泛的应用前景：常规重力场数据处理(网格化、曲化平、延拓、求导等)可以采用统一的等效场源正演实现，更符合位场规律；高精度重构结果，可以为重力场特征分析，重力异常处理与解释提供良好的基础，提升成果的可信度.
(2) 本文在重力数据反演(密度成像)方法基础上，加入“预优”技术，求解多层等效源.“预优”技术能显著提高求解方程的效率，而且预优因子可以显著提高深部场源的敏感度，同时配合等效单元体积随深度逐渐增大的设计.因此，本文的多层等效源联合反演，能够使得不同深度重力场信息在各等效层中得到更加合理的反映，并且基本不增加计算难度.
(3) 在模型实验与实际数据重构中，本文方案具有很好地压制噪声和消弱边界效应的能力，对噪声的压制得益于第1层和第2层等效源的合理设置，过浅或者过深都会对重构产生不利结果，而消弱边界效应主要基于第3层等效源对背景场信息较好的拟合.
(4) 大量模型试验表明，当背景场在边界差异较大时，利用本文提出的方法向上延拓且上延高度较大时，仍然会产生一定边界效应，其主要原因在于边界以外数据缺失，若能对边界以外的数据进行有效预测，则可进一步减小边界效应.
致谢: 中国地质大学(武汉)地球物理与空间信息学院杜劲松博士、孙石达博士以及冯建辉、孙凯等同学对本研究提供了宝贵的建议，刘双博士在反演工作方面提供了指导帮助；3位匿名审稿人给予了宝贵意见；在此一并致谢. -
图 1 3种等效源单元体的重力效应衰减特征
Fig. 1. Gravity attenuation characteristics of three equivalent sources
图 2 二维组合模型及其产生的重力效应
a.密度模型及多层等效源设置示意图，其中斜杠填充区域表示等效源分布区间；b.各模型理论重力异常及其总和
Fig. 2. 2D synthetic example and theoretical gravity anomaly
图 3 观测剖面及750 m高度上重构结果及其与理论值的拟合差分布
图中左列为各方法在观测面上重构重力和重力梯度及其与理论值的拟合差；右列为各方法在上延至750 m高度的重力和重力梯度及其与理论值的拟合差
Fig. 3. Reconstruction results and their difference relative to theoretical values on observation surface and on the altitude of 750 m
图 4 三维模型组合与观测面以及理论重力场
a.三维模型示意图；b.起伏的观测面；c.观测面上的理论重力场；d.观测面上背景场；e.含噪声的模型重力场与背景场之和；f.对数功率谱曲线
Fig. 4. 3D synthetic example, observation surface and theoretical gravity anomaly
图 5 三维模型及背景理论重力场及其重构结果
a.无噪声理论模型与背景场；b.用图 4e重构的观测面重力场；c.h=1 000 m平面上的理论模型与背景场；d.用图 4e重构的h=1 000 m平面上结果；e.观测面上重构重力值与理论值拟合差；f.h=1 000 m平面上重构重力值与理论值拟合差
Fig. 5. The gravity field from non-noise 3D models and background and reconstructive results
图 6 三维模型重构计算(Rec)的重力梯度结果与理论值比较
左侧(1, 2)两列为起伏观测面上重构后计算(Rec)的重力梯度及其理论值；右侧(3, 4)两列为420 m高度上重构计算的(Rec)重力梯度及其理论值；a、b与c分别为Txz(北向水平梯度)，Tyz(东向水平梯度)和Tzz(垂向梯度)；图例单位为：重力值(mGal)
Fig. 6. Comparison between the reconstruction gravity gradients and theoretical values of 3D models and background fields
图 7 三维模型重力梯度重构结果与理论值拟合差
a、b分别表示观测面上与420 m高度上拟合差，(1~3)分别表Txz(北向水平梯度)、Tyz(东向水平梯度)和Tzz(垂向梯度)
Fig. 7. Difference between the reconstruction gravity gradients and theoretical values of 3D models and background fields
图 8 研究区地形、实测重力异常与重构结果
a.研究区地形等高线；b.实测重力异常；c.重构起伏地面上重力异常；d.重构海拔910 m平面上重力异常；e.对数功率谱曲线
Fig. 8. Topography, observed gravity anomaly and reconstructive results in study area
图 9 重构的重力梯度
a1、a2和a3分别为观测面上Txz(北向水平梯度)，Tyz(东向水平梯度)和Tzz(垂向梯度)；b1、b2和b3分别为上延至海拔910 m平面上的Txz, Tyz和Tzz
Fig. 9. Reconstructive gravity gradients
表 1 模型参数设置
Table 1. Designed parameters of models
模型编号 模型形状 模型中心距高程h=0
处平面深度(km)密度
(103 kg/m3)1 小方块 0.140 1.5 2 小方块 0.140 1.5 3 小方块 0.115 1.5 4 小方块 0.050 1.5 5 小方块 0.100 1.5 6 小方块 0.170 1.5 7 长方柱体 0.300 1.0 8 长方柱体 0.550 1.0 9 板状体 0.650 2.0 10 倾斜板状体 0.500 1.0 
下载: 导出CSV
-
[1] Barnes, G., Lumley, J., 2011.Processing Gravity Gradient Data.Geophysics, 76(2):133-147. https://doi.org/10.1190/1.3548548 [2] Bhattacharyya, B.K., Chan, K.C., 1977.Reduction of Magnetic and Gravity Data on an Arbitrary Surface Acquired in a Region of High Topographic Relief.Geophysics, 42(7):1411-1430. https://doi.org/10.1190/1.1440802 [3] Chen, J.L., Zhang, Y., Chen, C., 2010.Satellite Magnetic Anomaly and Regional Geological Characteristics in Guangxi Area.Chinese Journal of Engineering Geophysics, 7(3):327-332(in Chinese with English abstract). [4] Dampney, C.N.G., 1969.The Equivalent Source Technique.Geophysics, 34(1):39-53. https://doi.org/10.1190/1.1439996 [5] Feng, J.F., 2016.Extraction and Explanation of Gravity and Magnetic Anomaly Features in Shedong Area of Guangxi(Dissertation).China University of Geosciences, Wuhan (in Chinese with English abstract). [6] Guo, Z.H., Guan, Z.N., Xiong, S.Q., 2004.Cuboid ΔT and Its Gradient Forward Theoretical Expressions without Analytic Odd Points.Chinese Journal of Geophysics, 47(6):1131-1138(in Chinese with English abstract). doi: 10.1007/978-3-642-21815-6_9/fulltext.html [7] Hansen, R.O., Miyazaki, Y., 1984.Continuation of Potential Fields between Arbitrary Surfaces.Geophysics, 49(6):787-795. https://doi.org/10.1190/1.1441707 [8] Klees, R., Tenzer, R., Prutkin, I., et al., 2008.A Data-Driven Approach to Local Gravity Field Modelling Using Spherical Radial Basis Functions.Journal of Geodesy, 82(8):457-471. https://doi.org/10.1007/s00190-007-0196-3 [9] Leão, J.W.D., Silva, J.B.C., 1989.Discrete Linear Transformations of Potential Field Data.Geophysics, 54(4):497-507. https://doi.org/10.1190/1.1442676 [10] Li, S., Liu, C.C., Jiang, Y.L., 2015.Application of 3D Equivalent Source Inversion and Imaging.Coal Geology & Exploration, 43(5):90-94 (in Chinese with English abstract). [11] Li, Y.G., Oldenburg, D.W., 1996.3-D Inversion of Magnetic Data.Geophysics, 61(2):394-408. https://doi.org/10.1190/1.1443968 [12] Li, Y.G., Oldenburg, D.W., 1998.3-D Inversion of Gravity Data.Geophysics, 63(1):109-119. https://doi.org/10.1190/1.1444302 [13] Li, Y.G., Oldenburg, D.W., 2000.Reduction to the Pole Using Equivalent Sources.SEG Technical Program Expanded Abstracts 2000, Technical Program Expanded Abstracts, 19(1):2484-2484. https://doi.org/10.1190/1.1816074 [14] Li, Y.G., Oldenburg, D.W., 2001.Stable Reduction to the Pole at the Magnetic Equator.Geophysics, 66(2):571-578. https://doi.org/10.1190/1.1444948 [15] Li, Y.G., Oldenburg, D.W., 2010.Rapid Construction of Equivalent Sources Using Wavelets.Geophysics, 75(3):L51-L59. https://doi.org/10.1190/1.3378764 [16] Liu, S., Feng, J., Gao, W.L., et al., 2013.2D Inversion for Borehole Magnetic Data in the Presence of Significant Remanence and Demagnetization.Chinese Journal of Geophysics, 56(12):4297-4309(in Chinese with English abstract). [17] Liu, T.Y., Yang, Y.S., Li, Y.Y., et al., 2007.The Order-Depression Solution for Large-Scale Integral Equation and Its Application in the Reduction of Gravity Data to a Horizontal Plane.Chinese Journal of Geophysics, 50(1):290-296(in Chinese with English abstract). http://www.oalib.com/paper/1567950 [18] Ma, G.Q., Ming, Y.B., He, Y., et al., 2016.Horizontal Derivative Iteration Method for Downward Continuation of Gravity and Magnetic Data.Earth Science, 41(7):1231-1237 (in Chinese with English abstract). [19] MacLennan, K., Li, Y.G., 2013.Denoising Multicomponent CSEM Data with Equivalent Source Processing Techniques.Geophysics, 78(3):E125-E135. https://doi.org/10.1190/geo2012-0226.1 [20] Mendonça, C.A., Silva, J.B.C., 1994.The Equivalent Data Concept Applied to the Interpolation of Potential Field Data.Geophysics, 59(5):722-732. https://doi.org/10.1190/1.1443630 [21] Nagy, D., 1966.The Gravitational Attraction of a Right Rectangular Prism.Geophysics, 31(2):362-371. https://doi.org/10.1190/1.1439779 [22] Needham, P.E., 1970.The Formation and Evaluation of Detailed Geopotential Models Based on Point Masses.Report 149.Department of Geodetic Science and Surveying, Ohio State University, Ohio. [23] Oliveira, Jr., V.C., Barbosa, V.C.F., 2013.Polynomial Equivalent Layer.Geophysics, 78(1):G1-G13. https://doi.org/10.1190/geo2012-0196.1 [24] Pang, X.L., 2012.Research on Reduction of Aeromagnetic Anomalies by Means of Equivalent Source Technology (Dissertation).China University of Geosciences, Beijing (in Chinese with English abstract). [25] Pang, X.L., Yao, C.L., Zheng, Y.M., et al., 2010.The Technique of the Altitude Correction on Aeromagnetic Data.In:Chinese Geophysical Society, Seismological Society of China, eds., The Chinese Geophysics 2010-Proceedings of the Twenty-Sixth Annual Conference of the Chinese Geophysical Society and the Thirteenth Academic Conference of the Chinese Seismological Society.Seismological Press, Beijing, 631 (in Chinese). [26] Pang, X.L., Yao, C.L., Zheng, Y.M., et al., 2011.The Improvements of the Technique of the Altitude Correction on Aeromagnetic Data.In:Chinese Geophysical Society, ed., Proceedings of the Twenty-Seventh Annual Conference of the Chinese Geophysical Society.University of Science and Technology of China Press, Hefei, 697 (in Chinese). [27] Pilkington, M., 1997.3-D Magnetic Imaging Using Conjugate Gradients.Geophysics, 62(4):1132-1142. https://doi.org/10.1190/1.1444214 [28] Sun, S.D., Chen, C., Du, J.S., et al., 2016.Magnetic Characteristics and Tectonic Implications of Crust in Junggar Basin and Its Surroundings.Earth Science, 41(7):1216-1224 (in Chinese with English abstract). [29] Sünkel, H., 1981.Point Mass Models and the Anomalous Gravitational Feld.Report 328.Department of Geodetic Sciences, Ohio State University, Columbus. [30] Syberg, F.J.R., 1972.A Fourier Method for the Regional-Residual Problem of Potential Fields.Geophysical Prospecting, 20(1):47-75.https://doi.org/10.1111/j.1365-2478.1972.tb00619.x doi: 10.1111/gpr.1972.20.issue-1 [31] Vermeer, M., 1991.Geoid Recovery at 0.5 Degree Resolution from Global Satellite Gradiometry Data Sets.International Association of Geodesy Symposia, 18-29. https://doi.org/10.1007/978-1-4612-3104-2_5 [32] Wang, W.Y., Pan, Z.S., 1996.Data Processing and Transformation Methods of Potential Fields on an Arbitray Surface with Dipole Layer Potential Technique.Journal of Xi'an College of Geology, 18(3):69-76 (in Chinese with English abstract). https://www.sciencedirect.com/science/article/pii/B9780444638908000347 [33] Wang, W.Y., Pan, Z.S., Li, J.K., 1991.Continuation Methods for Curved Surface of the Three-Dimensional High-Precision Gravity and Magnetic Potential Field.Geophysical & Geochemical Exploration, 15(6):415-422 (in Chinese with English abstract). [34] Wei, H.P., 2014.The Equivalent Source Method in Weak Gravity and Magnetic Information on Research and Application in Separation(Dissertation).Chengdu University of Technology, Chengdu (in Chinese with English abstract). [35] Wu, X.P., 1984.Point-Mass Model of Local Gravity Field.Acta Geodetica et Cartographica Sinica, 13(4):250-258 (in Chinese with English abstract). [36] Xia, J., Sprowl, D.R., 1991.Correction of Topographic Distortion in Gravity Data.Geophysics, 56(4):537-541. https://doi.org/10.1190/1.1443070 [37] Xie, R.K., Wang, P., Duan, S.L., et al., 2015.Analysis of the Reduction of Aeromagnetic Gradients Data to a Horizontal Plane.Progress in Geophysics, 30(6):2836-2840 (in Chinese with English abstract). [38] Xu, S.Z., Shen, X.H., Zou, L.J., et al., 2004.Downward Continuation of Aeromagnetic Anomaly from Flying Altitude to Terrain.Chinese Journal of Geophysics, 47(6):1127-1130(in Chinese with English abstract). [39] Zhao, D.M., 2001.Research on the Fast Approximation of Ballistic Missile Disturbing Gravity(Dissertation).Information Engineering University, Zhengzhou (in Chinese with English abstract). [40] Zhu, G.Q., Li, H.L., Wen, R.X., 2011.Prediction and Analysis of Deep Mineral Exploration in Guangxi.Chinese Journal of Engineering Geophysics, 8(6):713-722(in Chinese with English abstract). doi: 10.1007/s11771-009-0164-6 [41] 陈家联, 张毅, 陈超, 2010.广西地区卫星磁异常与区域地质特征.工程地球物理学报, 7(3):327-332. http://www.cnki.com.cn/Article/CJFDTOTAL-WTYH200305001.htm [42] 冯建辉, 2016. 广西社垌地区重磁异常特征提取与解释(硕士学位论文). 武汉: 中国地质大学. [43] 郭志宏, 管志宁, 熊盛青, 2004.长方体ΔT场及其梯度场无解析奇点理论表达式.地球物理学报, 47(6):1131-1138. http://www.cqvip.com/Main/Detail.aspx?id=11420770 [44] 黎莎, 刘春成, 江玉乐, 2015.三维等效源反演成像研究及应用.煤田地质与勘探, 43(5):90-94. http://industry.wanfangdata.com.cn/dl/Detail/Periodical?id=Periodical_mtdzykt201505021 [45] 刘双, 冯杰, 高文利, 等, 2013.强剩磁强退磁条件下的二维井中磁测反演.地球物理学报, 56(12):4297-4309. doi: 10.6038/cjg20131232 [46] 刘天佑, 杨宇山, 李媛媛, 等, 2007.大型积分方程降阶解法与重力资料曲面延拓.地球物理学报, 50(1):290-296. doi: 10.3321/j.issn:0001-5733.2007.01.036 [47] 马国庆, 明彦伯, 贺杨, 等, 2016.重磁数据稳定向下延拓的水平导数迭代法.地球科学, 41(7):1231-1237. http://www.earth-science.net/WebPage/Article.aspx?id=3332 [48] 庞旭林, 2014. 航磁异常数据曲面延拓等效源法技术研究(硕士学位论文). 北京: 中国地质大学. [49] 庞旭林, 姚长利, 郑元满, 等, 2010. 基于航磁测线数据高度改正方法研究. 见: 中国地球物理学会、中国地震学会编, 中国地球物理2010——中国地球物理学会第二十六届年会、中国地震学会第十三次学术大会论文集. 北京: 地震出版社, 631. [50] 庞旭林, 姚长利, 郑元满, 等, 2011. 基于航磁测线数据高度改正技术的改进. 见: 中国地球物理学会编, 中国地球物理学会第二十七届年会论文集. 合肥: 中国科学技术大学出版社, 697. [51] 孙石达, 陈超, 杜劲松, 等, 2016.准噶尔盆地及邻区地壳磁性特征及其构造意义.地球科学, 41(7):1216-1224. http://www.earth-science.net/WebPage/Article.aspx?id=3330 [52] 王万银, 潘作枢, 1996.偶层位曲面位场数据处理及转换方法.西安地质学院学报, 18(3):69-76. http://www.cqvip.com/qk/94561x/199603/2296336.html [53] 王万银, 潘作枢, 李家康, 1991.三维高精度重磁位场曲面延拓方法.物探与化探, 15(6):415-422. http://www.cnki.com.cn/Article/CJFDTOTAL-WTYH198802002.htm [54] 魏华平, 2014. 等效源法在重磁弱信息分离中的研究与应用(硕士学位论文). 成都: 成都理工大学. [55] 吴晓平, 1984.局部重力场的点质量模型.测绘学报, 13(4):249-258. http://kns.cnki.net/KCMS/detail/detail.aspx?filename=chxb198404001&dbname=CJFD&dbcode=CJFQ [56] 谢汝宽, 王平, 段树岭, 等, 2015.航磁梯度数据曲化平分析.地球物理学进展, 30(6):2836-2840. doi: 10.6038/pg20150650 [57] 徐世浙, 沈晓华, 邹乐君, 等, 2004.将航磁异常从飞行高度向下延拓至地形线.地球物理学报, 47(6):1127-1130. doi: 10.3321/j.issn:0001-5733.2004.06.028 [58] 赵东明, 2001. 弹道导弹扰动引力快速逼近的算法研究(硕士学位论文). 郑州: 信息工程大学. [59] 朱国器, 黎海龙, 温融湘, 2011.广西深部找矿特征分析与找矿预测.工程地球物理学报, 8(6):713-722. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=gcdqwlxb201106014 期刊类型引用(1)
1. 高宝龙，胡正旺，李端，杜劲松. 多层等效源方法在地面与航空磁异常数据融合中的应用. 地球科学. 2021(05): 1881-1895 . 
本站查看其他类型引用(0)
-
下载:
点击查看大图
计量
- 文章访问数: 4091
- HTML全文浏览量: 1555
- PDF下载量: 14
- 被引次数: 1
