2018, Vol. 61
2. 北京派特森科技股份有限公司, 北京 100015;
3. 东北大学资源与土木工程学院, 沈阳 110819;
4. 国防科技大学机电工程与自动化学院, 长沙 410000
2. Beijing Petrosound Geoservices Co., Ltd, Beijing 100015, China;
3. School of Resources and Civil Engineering, Northeastern University, Shenyang 110819, China;
4. Department of Mechatronic Engineering and Automation, National University of Defense Technology, Changsha 410000, China
重力场与相关的梯度数据联合处理和解释中的优势在于能够迅速和准确的推断地下地质构造的场源边界、相互关系、深度变化、分布规模、场源性质等地球物理现象(马国庆等, 2012a, 2014；周帅等，2016).尤其是，目前广泛开展的不同空间探测技术的广泛应用，迫切需要一种能够将重力场异常和梯度数据同时进行自动反演和解释的计算方法与技术(Blakely, 1995)，提升发现目标的能力.一种好的位场反演方法应该是在较少先验地质信息条件下，解决一定的实际地质问题，这样就能最大限度地减少人为推断误差，而欧拉反褶积法就是这样一种方法.
欧拉反褶积方法作为一种能够同时利用位场和梯度数据完成自动估算场源位置的计算分析方法，近十多年来被广泛应用，取得了良好效果.理论上，它以欧拉齐次方程为基础，运用场源数据和三个方向空间导数，以及各种不同地质体特有的与场源性质有关的构造指数来确定观测数据与场源位置和性质的对应关系.它的优势在于能够较灵活的适应地质对象变化情况，不需要精确假设地质模型形态，不需要考虑磁场的磁化方向，对同源重磁场具有联合处理解释的融合能力，能在较少地质先验信息情况下独立完成自动解释(Reid，2007).欧拉法一开始提出就是根据齐次的欧拉方程用构造指数来求出场源的水平位置和垂直埋深.提出时只能用于以点源和偶极源的反演.经过Thompson(1982)的扩展，推算，得到了二维的反演方法.后来，Reid等(1990)又把欧拉反褶积方法的应用拓展到了三维.Zhang(2000)等提出张量欧拉反褶积法，并很快得到进一步的扩展和推广(Mushayandebvu et al., 2001; Ravat et al., 2002)，对于更加多样化的场源类型都得到很好的应用(Stavrev and Reid, 2007, 2010).Salem等(2008)将倾斜角和欧拉反褶积联合用于磁数据的自动解释.Reid和其他作者相继做了很多关于欧拉反褶积的工作(Reid et al., 2012; Reid and Thurston, 2014), 并对构造指数的选取问题进行了详细讨论.近年来，Roth(2009)对海洋全张量重力数据进行分析研究取得了较好的效果.周文月等(2017)将航空重力全张量数据欧拉反褶积方法用于文顿盐丘地区的研究，有效的圈定地下异常体范围.Fedi等提出一种多尺度的分析方法并将其应用于欧拉反褶积(Fedi et al., 2009, 2015; Florio and Fedi, 2014)，这种新的方法可以很好的获取欧拉解，在一定程度上提高了欧拉反演的精度.
针对以上研究，本文提出多个不同高度数据联合欧拉反褶积方法.基于不同高度的观测数据可同时对应同一场源并用于反演场源位置的分析原理，拓展不同高度场数据在欧拉反褶积方法中的应用范围.不同高度数据可以通过延拓或者航空、潜航数据获得.向上延拓数据和航空数据，可以压制噪声，滤除不需要的小局部异常体，适合大面积的重磁勘探；向下延拓数据和潜航数据，在深度上更容易靠近异常体，适合寻找发现某些特定的目标体.此外，欧拉反褶积方程在不同的高度表现形式相似，不同高度数据反应同一个场源的位置信息，因此可将不同高度数据融合并实现联合反演.在实际应用中，该方法不仅可以解决场源的定性问题、同时可以解决场源位置变化的定量计算问题，可有效圈出异常源的基本轮廓，尤其适合于大面积网格重力数据计算和大深度分析和解释.
1 不同高度数据欧拉公式的提出与系统化研究三维欧拉方程由Thompson(1982)提出，它运用场源数据和其三个方向空间导数，以及各种不同地质体特有的与场源性质有关的构造指数来估算场源位置和深度.公式如下：
|
(1) |
式中，(x0, y0, z0)是场源位置，(x, y, z)为观测点坐标，Bz和N分别是背景场和构造指数，fz是重力异常，fzx, fzy和fzz分别是重力异常在x, y和z三个方向的导数.移向后上式可写为
|
(2) |
通过改变高度因子h对(2)式进行改进，设另外两个不同的观测高度分别为z1, z2, 观测面z1上的观测点坐标为(x, y, z1), 其中z1=z+h1, 重力及其梯度异常值分别为：fz1, fxz1, fyz1, fzz1; 观测面z2上的观测点坐标为(x, y, z2), 其中z2=z+h2, 重力及其梯度异常值分别为：fz2, fxz2, fyz2, fzz2.则由公式(2)可以得到观测高度z1平面的欧拉公式(3)：
|
(3) |
对比分析(2)和(3)式可以看出，欧拉反褶积方程在不同的高度其表现形式相似，且两个不同高度反应同一个场源的位置信息，因此将不同高度数据融合来实现多种高度数据联合反演计算.两个不同高度(z和z1观测面)数据联合欧拉求解矩阵表达式为
|
(4) |
同理，对z2观测面而言，可以得到类似公式(3)表达式，三个观测面z, z1, z2可形成联合表达式，即三个不同高度数据联合欧拉求解矩阵表达式为
|
(5) |
对于m+1个不同高度的数据联合欧拉求解矩阵表达式可表示为
|
(6) |
公式(6)构成在某一水平测量节点(x，y)上的m+1个不同高度重力场和梯度数据联合反演的欧拉反褶积求解矩阵方程.
2 理论模型试验 2.1 单一长方体加小局部异常在笛卡尔直角坐标系下，以x轴正方向为正北方向，y轴正方向为正东方向，z轴垂直向下.选择测区范围为x方向0~4000 m，y方向0~4000 m.采样间隔为40 m.建立如下的长方体模型：棱柱体长1000 m，宽1000 m，高200 m，中心埋深(2000, 2000, 200)m，异常体密度为1000 kg·m-3.在上述模型中加入一小的局部异常，局部异常体为边长为200 m的正方体，中心点坐标为(3000, 1000, 270)m，异常体密度为1000 kg·m-3.观测面在z=0时，欧拉反褶积平面结果能显示出局部异常的存在(图 1a).向下延拓50 m时，反演结果成功显示出局部异常，且结果比较收敛(图 1b).这是因为随着向下延拓深度的增加，观测面离异常体位置变近，可测到的重力及其梯度异常值增大.说明向下延拓能够突出有效局部小异常.
|
图 1 单一观测面有局部异常体存在时欧拉反演结果 (a)观测面在z=0 m时，欧拉反褶积平面结果；(b)观测面在z1=50 m时，欧拉反褶积平面结果. Fig. 1 Euler result with local anomaly of one single observation data (a) Euler result when the observation surface is z=0 m; (b) Euler result when the observation surface is z1=50 m. |
图(2a，2b)为有噪声存在时的欧拉反演结果.从图 2a可以看到，观测面为地面时，噪声和小干扰异常都存在，反演结果不是很好；当向上延拓50 m后，成功滤掉随机噪声和小干扰异常.这是因为，高度增加，观测面离异常体位置变远，可测到的重力及梯度异常值减小，噪声和小局部异常体产生的干扰在向上延拓到一定高度后弱化或消失.说明向上延拓能够压制噪声和小的局部干扰异常，降低局部干扰异常对反演结果聚焦准确度的干扰力度.
|
图 2 单一观测面有噪声存在时欧拉反演结果 (a)观测面在z=0 m时，欧拉反褶积平面结果；(b)观测面在z2=-50 m时，欧拉反褶积结果. Fig. 2 Euler result with noise of one single observation data (a) Euler result when the observation surface is z=0 m; (b) Euler result when the observation surface is z2=-50 m. |
在笛卡尔直角坐标系下，x轴正方向设为正北方向，y轴正方向为正东方向，z轴垂直向下.选择测区范围为x方向0~4000 m，y方向0~4000 m.采样间隔为40 m.建立如下的长方体模型：两个棱柱体边长分别为长1000 m，宽1000 m，高200 m，中心点坐标分别为(1500, 1500, 200)m，(2500, 2500, 300)m，异常体密度为1000 kg·m-3.
2.2.1 顶点接触图 3a为两个埋深不同的长方体重力异常数据，图 3b为窗口选用7个采样间隔时单一观测面欧拉反演平面结果，可以看到埋深较浅的地质体识别效果较好，而深部地质体边界识别效果一般.考虑到高度增加，观测面离异常体位置变远，可测到的有效重力及其梯度异常值较少，采用加大窗口的方法增加有效解.图 3c为窗口大小增大为9个采样间隔时，欧拉反演的结果.增大窗口后，解的数量明显增加，两个异常体的位置均能有效识别.
|
图 3 模型重力异常及欧拉反演平面结果 (a)重力异常；(b)窗口大小7×7时单一观测面数据反演结果；(c)窗口大小9×9时单一观测面数据反演结果；(d)窗口大小7×7，两个高度数据联合反演结果. Fig. 3 Gravity anomaly of model and Euler deconvolution result (a) Gravity anomaly; (b) Euler result when the window size is 7×7; (c) Euler result when the window size is 9×9; (d) Joint Euler result of two different heights observation data when the window size is 7×7. |
单一高度数据加大窗口后，能大致识别出边界信息，但欧拉反演结果比较发散.而从图 3d中可以看到，两个高度数据同时反演时，不用过多加大窗口的情况下，就能得到两个地质体的边界信息.不同高度数据联合欧拉反演，不但能识别出两个地质体边界接触的角点信息，而且反演结果明显比单一观测面数据反演结果收敛程度好.
2.2.2 相互接触模型将两个棱柱体的中心点坐标分别设为(1500, 1500, 200)m和(2500, 2000, 300)m，其他参数均不变.模型重力异常及欧拉反演结果如图 4所示.不同高度观测面z=0 m, z2=-50 m以及z1=50 m数据单独反演时，如图 4b、图 4c和图 4e所示，反演结果能够大致识别出较浅的异常体，对于埋深较深的异常体只能识别出一部分边界，对两个异常接触附近的边界不能识别.这是因为单独不同高度数据反演时，不同场源叠加到同一窗口，移动窗口把某些深部有效异常作为干扰处理，使得欧拉反演结果不准确.因此采用不同高度数据的联合来完成欧拉反演.当用几个不同高度数据联合反演时，结果如图 4d和图 4f所示：不同高度数据融合联合欧拉反褶积结果能够圈定出两个接触长方体的边界，并且能识别出相接触的两个角点位置(2000, 1500)m和(2000, 2000)m信息.这是任何一个单一观测面数据或者合成几个观测面单独计算结果所不能达到的.因此，不同高度数据融合来进行联合欧拉反褶积能够有效地改进单一观测面数据结果的准确性，提高反演结果的准确度.
|
图 4 模型重力异常及欧拉反演结果 (a)重力异常；(b)观测面z=0 m欧拉反演结果；(c)观测面z2=-50 m欧拉反演结果；(d)两个观测面数据联合反演(窗口7×7)；(e)观测面z1=50 m欧拉反演结果；(f)三个不同高度数据联合欧拉反演结果(窗口7×7). Fig. 4 Gravity anomaly of model and Euler deconvolution result (a) Gravity anomaly; (b) Euler result when the observation surface is z=0 m; (c) Euler result when the observation surface is z2=-50 m; (d) Joint Euler result of two different heights observation data when the window size is 7×7; (e) Euler result when the observation surface is z1=50 m; (f) Joint Euler result of three different heights observation data when the window size is 7×7. |
为了进一步描述深度结果聚焦程度，对不同高度欧拉反演结果进行深度统计.根据模型深度范围进行结果统计，统计解的个数及反演出的异常体的平均埋深并对解得收敛程度进行评估，最后对反演结果的质量做出评价.
对解的收敛程度进行评估(参考Marson and Kingele, 1993)，和方差的定义相似，公式如下：
|
(7) |
式中，X代表欧拉反演结果，μ是模型的真实深度，E[]代表X和μ最小二乘的期望值.方差越小表明越多的欧拉反褶积的解接近于真实深度.
表 1给出了构造指数为0，窗口大小为5个采样间隔，模型中心埋深为200 m时，不同高度数据欧拉反褶积计算结果.可以得出，不同高度观测数据都可以很好的反应出场源信息，但延拓到一定高度后，深度方向方差随延拓高度的增加而增大；几个不同高度数据联合反演后方差介于单独高度反演方差的最大值和最小值之间.综上所述，不同高度联合欧拉反褶积在一定程度上可以使深度反演结果更接近于真实的深度.
|
表 1 欧拉反演结果评价 Table 1 The evaluation of Euler deconvolution result |
为了说明不同高度数据融合联合欧拉反褶积反演方法应用的有效性.结合来自龙门山断裂地区的实测重力数据进行分析研究.通过广泛收集整理龙门山地区地球物理探测资料，分析该地区重力异常分布特征，然后利用欧拉反褶积法划分龙门山地区及邻区断裂分布特征及参考深度.
3.1 研究区区域重力场特征根据龙门山地区1:250万布格重力异常图，对该地区重力异常特征进行了描述；并通过参考前人研究的地质资料和断裂划分的位置，对研究区进行了区域划分.
图 5是研究区的布格重力异常场，可以看到研究区异常值均为负值，这是因为布格改正、地形改正后，等同于排除了大地水准面以上的物质质量，引起地壳质量不足；异常值整体上自北西向南东逐渐升高；东西两侧异常形态变化比较平缓，而中间是重力异常梯级带.根据重力异常值的变化，对研究区进行了划分(图 5).
|
图 5 研究区布格重力异常图及区域的划分 Fig. 5 The Bouguer gravity anomaly and region division |
主要分了三个单元的区域：一级单元界限(粗线)位于龙门山断裂带上，由于在此处重力异常呈现近北东走向的高梯度带，是布格重力异常最为密集的位置.汶川大地震的几个重大发生地点如汶川、茂县、北川等都是在这条界限上.在一级单元界限的两侧，划分了马尔康块体(Ⅱ1)和茂县—丹巴块体(Ⅱ2).一级单元界限的东南方向内，又将其进行详细的划分.首先，按照重力异常的疏密程度划分出了另一条二级单元界限，界限的西部是近于南北向的条带状重力异常，变化幅度较大且变化较快，大概50 km的距离变化了60×10-5 m·s-2, 为龙门山断裂带(Ⅱ3)；将二级界限的东部划为川西台陷(Ⅲ1)，川中台拱(Ⅲ2)和川北台陷(Ⅲ3)三部分(马国庆等, 2012b, c).
3.2 断裂分布特征用欧拉反褶积法对龙门山地区重力数据进行反演，由表 2可知，构造指数选0.根据反演结果可以划分得到龙门山及邻区断裂的分布特征.如图 6所示，可以看到该地区存在深大断裂.
|
图 6 重力异常欧拉反褶积平面结果 Fig. 6 The Euler deconvolution result of gravity anomaly |
|
|
表 2 欧拉构造指数表 Table 2 The structure index of Euler |
为进一步研究龙门山地区深大断裂，获得断裂大致深度，对布格重力异常向上延拓不同高度，得到不同高度重力数据进行欧拉反演.将地面观测数据、向上延拓5000 m的数据以及向上延拓10000 m的数据融合.上延后仍存在重力异常说明有断裂存在且断裂深度大于此上延高度.三个高度联合欧拉反褶积结果如图 7所示.可知龙门山地区有深大断裂，走向近于北东，深达20 km以上.
|
图 7 龙门山地区三个不同高度数据联合欧拉反褶积平面结果 Fig. 7 Joint Euler result of three different heights observation data of Longmenshan |
根据不同高度数据联合欧拉反褶积法反演结果，对该地区断裂进行了划分，如图 8所示.依据平面欧拉结果可以划出研究区内断裂的分布规律.从图 8中可以看出，研究区共存在14条断裂，有8条为北东走向，北西向3条，近东西向3条.龙门山三条断裂分别是F1、F2、F3，以这三条主断裂为例，描述研究区断裂特征并与前人研究地质资料对比.如表 3所示，可以得到不同高度数据联合欧拉反演结果与地质资料对应较好.从联合后重力异常反演结果可知，龙门山地区三条断裂带延深较大，由于物质运移是发生在上地幔，因此这三条断裂带为物质的运移通道.
|
图 8 断裂划分示意图 F1汶川—茂县断裂；F2映秀—北川断裂；F3都江堰—江油断裂. Fig. 8 The characteristic of faults F1 Wenchuan-Maoxian Fault; F2 Yingxiu-Beichuan Fault; F3 Dujiangyan-Jiangyou Fault. |
|
|
表 3 不同高度数据联合欧拉反演得到的断裂特征与地质资料的对比 Table 3 Comparison between actual geological data and fault recognition results of different heights Euler deconvolution |
本文讨论了不同高度数据在欧拉反褶积法解决场源问题中的应用，研究发现了增加不同高度数据联合计算解释能够改进移动窗口内数据量对结果的影响.分析结果表明，数据量越大越能满足统计规律，收敛过程由此得到改善.针对不同高度数据融合并实现联合反演问题，将不同高度异常范围变化作为窗口选取的参考依据，从中确保了窗口内数据使用的统计前提条件，也避免了由于窗口选择过大对反演结果的精细程度产生的影响.并结合来自龙门山断裂地区的重力数据进行分析研究，划分了断裂，取得了良好的结果.
致谢感谢导师黄大年教授在论文选题上的指导，感谢吉林大学地球探测科学与技术学院马国庆副教授提供的数据支持，感谢审稿专家及编辑对本文提出的宝贵的修改意见.
Blakely R J. 1995. Potential Theory in Gravity and Magnetic Applications. Cambridge: Cambridge University Press.
|
Fe di, M, Florio G, Quarta T A M. 2009. Multiridge analysis of potential fields:Geometric method and reduced Euler deconvolution. Geophysics, 74(4): L53-L65. DOI:10.1190/1.3142722 |
Fedi M, Florio G, Paoletti V. 2015. MHODE:A local homogeneity theory for improved source-parameter estimation of potential fields. Geophysical Journal International, 202(2): 887-900. DOI:10.1093/gji/ggv185 |
Florio G, Fedi M. 2014. Multiridge Euler deconvolution. Geophysical Prospecting, 62(2): 333-351. DOI:10.1111/gpr.2014.62.issue-2 |
Ma G Q, Du X J, Li L L. 2012. Comparison of the tensor local wavenumber method with the conventional local wavenumber method for interpretation of total tensor data of potential fields. Chinese Journal of Geophysics (in Chinese), 55(7): 2450-2461. DOI:10.6038/j.issn.0001-5733.2012.07.029 |
Ma G Q, Huang D N, Yu P, et al. 2012. Application of improved balancing filters to edge identification of potential field data. Chinese Journal of Geophysics (in Chinese), 55(12): 4288-4295. DOI:10.6038/j.issn.0001-5733.2012.12.040 |
Ma G Q, Meng L S, Li L L. 2012. Fault distribution of Longmenshan and adjacent regions and fault morphological differences before and after the earthquake. Journal of Jilin University (Earth Science Edition) (in Chinese), 42(2): 519-525. |
Ma G Q, Huang D N, Li L L, et al. 2014. A normalized local wavenumber method for interpretation of gravity and magnetic anomalies. Chinese Journal of Geophysics, 57(4): 1300-1309. DOI:10.6038/cjg20140427 |
Marson I, Kingele E E. 1993. Advantages of using the vertical gradient of gravity for 3-D interpretation. Geophysics, 58(11): 1588-1595. DOI:10.1190/1.1443374 |
Mushayandebvu M F, Van Driel P, Reid A B, et al. 2001. Magnetic source parameters of two-dimensional structures using extended Euler deconvolution. Geophysics, 66(3): 814-823. DOI:10.1190/1.1444971 |
Ravat D, Wang B, Wildermuth A, et al. 2002. Gradients in the interpretation of satellite-altitude magnetic data:An example from central Africa. Journal of Geodynamics, 33(1-2): 131-142. DOI:10.1016/S0264-3707(01)00059-X |
Reid A B, Allsop J M, Granser H, et al. 1990. Magnetic interpretation in three dimensions using Euler deconvolution. Geophysics, 55(1): 80-91. DOI:10.1190/1.1442774 |
Reid A B. 2007. Euler deconvolution. //Gubbins D, Herrero-Bervera E eds. Encyclopedia of Geomagnetism and Paleomagnetism. Netherlands: Springer, 263-265.
|
Reid A B, Ebbing J, Webb S J. 2012. Comment on 'A crustal thickness map of Africa derived from a global gravity field model using Euler deconvolution' by Getachew E. Tedla, M. Van Der Meijde, A. A. Nyblade and F. D. Van Der Meer. Geophysical Journal International, 189(3): 1217-1222. |
Reid A B, Thurston J B. 2014. The structural index in gravity and magnetic interpretation:Errors, uses, and abuses. Geophysics, 79(4): J61-J66. DOI:10.1190/geo2013-0235.1 |
Roth M. 2009. Marine full tensor gravity gradiometry data analysis and Euler deconvolution[Ph. D. Thesis]. Stuttgart: University of Stuttgart.
|
Salem A, Williams S, Fairhead D, et al. 2008. Interpretation of magnetic data using tilt-angle derivatives. Geophysics, 73(1): L1-L10. DOI:10.1190/1.2799992 |
Stavrev P, Reid A. 2007. Degrees of homogeneity of potential fields and structural indices of Euler deconvolution. Geophysics, 72(1): L1-L12. DOI:10.1190/1.2400010 |
Stavrev P, Reid A. 2010. Euler deconvolution of gravity anomalies from thick contact/fault structures with extended negative structural index. Geophysics, 75(6): 151-158. DOI:10.1190/1.3496902 |
Thompson D T. 1982. A new technique for making computer-assisted depth estimates from magnetic data. Geophysics, 47(1): 31-37. DOI:10.1190/1.1441278 |
Zhang C Y, Mushayandebvu M F, Reid A B, et al. 2000. Euler deconvolution of gravity tensor gradient data. Geophysics, 65(2): 512-520. DOI:10.1190/1.1444745 |
Zhou S, Huang D N, Jiao J. 2016. A filter to detect edge of potential field data based on three-dimensional structural tensors. Chinese Journal of Geophysics (in Chinese), 59(10): 3847-3858. DOI:10.6038/cjg20161028 |
ZHOU Wen-Yue, MA Guo-Qing, HOU Zhen-Long, et al. 2017. The study on the joint Euler deconvolution method of full tensor gravity data. Chinese Journal of Geophysics (in Chinese), 60(12): 4855-4865. DOI:10.6038/cjg20171225 |
马国庆, 杜晓娟, 李丽丽. 2012a. 解释位场全张量数据的张量局部波数法及其与常规局部波数法的比较. 地球物理学报, 55(7): 2450-2461. DOI:10.6038/j.issn.0001-5733.2012.07.029 |
马国庆, 黄大年, 于平, 等. 2012b. 改进的均衡滤波器在位场数据边界识别中的应用. 地球物理学报, 55(12): 4288-4295. DOI:10.6038/j.issn.0001-5733.2012.12.040 |
马国庆, 孟令顺, 李丽丽. 2012c. 龙门山及邻区断裂分布及地震前后断裂形态差异. 吉林大学学报(地球科学版), 42(2): 519-525, 553. |
马国庆, 黄大年, 李丽丽, 等. 2014. 重磁异常解释的归一化局部波数法. 地球物理学报, 57(4): 1300-1309. DOI:10.6038/cjg20140427 |
周帅, 黄大年, 焦健. 2016. 基于三维构造张量的位场边界识别滤波器. 地球物理学报, 59(10): 3847-3858. DOI:10.6038/cjg20161028 |
周文月, 马国庆, 侯振隆, 等. 2017. 重力全张量数据联合欧拉反褶积法研究及应用. 地球物理学报, 60(12): 4855-4865. DOI:10.6038/cjg20171225 |