地球物理学报  2012, Vol. 55 Issue (1): 304-309   PDF    
引用本文
孟小红, 刘国峰, 陈召曦, 等. 基于剩余异常相关成像的重磁物性反演方法. 地球物理学报, 2012, 55(1): 304-309, doi: 10.6038/j.issn.0001-5733.2012.01.030.  
MENG X H, LIU G F, CHEN Z X, et al. 3-D gravity and magnetic inversion for physical properties based on residual anomaly correlation. Chinese J. Geophys. (in Chinese), 2012, 55(1): 304-309, doi: 10.6038/j.issn.0001-5733.2012.01.030.  
基于剩余异常相关成像的重磁物性反演方法
孟小红, 刘国峰 , 陈召曦, 郭良辉     
中国地质大学(北京)地球探测与信息技术教育部重点实验室, 北京 100083
收稿日期 2011-01-21, 2011-12-14 收修定稿
基金项目: 国家公益性研究项目Sino-prob01及03课题,国家自然科学基金项目(41074095,41104083);中央高校基本科研业务费(2011YYL022)联合资助
作者简介: 孟小红,女,教授,博士生导师,主要从事计算地球物理,位场三维反演等研究.E-mailmxh@cugb.edu.cn
通讯作者: 刘国峰,博士,中国地质大学(北京)讲师.E-mail:liugf@cugb.edu.cn
摘要: 将场源区剖分成长方体单元,通过采集的重磁数据反演出这些单元的密度或者磁化率变化,勾画出场源的分布图像,这种方式是重磁三维反演的重要方向.重磁相关成像通过计算测量的重磁异常与地下各点在测区上的重磁异常的归一化相关,显示出异常地质体的空间赋存状态和等效剩余重磁物性.该方法计算速度快,方法简单、稳定,但是反演的结果只是在-1到+1之间的等效物性,不能够直接反演剩余密度或者磁化率,并且无法引入已知的地质约束.本文通过对物性模型的正演和实测结果的残差进行相关成像,迭代更新物性模型实现对物性参数的反演过程.模型实验证明该方法相对相关成像不仅能提高分辨率,还能够得到真正的物性参数.
关键词: 相关成像      剩余异常      物性模型      反演     
3-D gravity and magnetic inversion for physical properties based on residual anomaly correlation
MENG Xiao-Hong, LIU Guo-Feng, CHEN Zhao-Xi, GUO Liang-Hui     
Geo-detection Laboratory, Ministry of Education of China, Beijing 100083, China
Abstract: It is a main research direction to outline the source distribution from the imaging result of the gravity or magnetic inversion based on rectangular units. Correlation imaging of gravity and magnetic is a method that calculates the normalized correlation of the anomaly of rectangular unit and acquired anomaly. The result can outline the properties distribution and earth structure. The correlation imaging method is characterized by fast computation, stability, and simplicity, but it results in values between -1 and 1 but not the real physical properties, and also can't use the prior geology information. In this paper, we proposed a method based on the correlation imaging of the residuals of acquired data and forward data based on a physical model; after some iterations, we get the physical property model. The model test verified the effectiveness of our method.
Key words: Correlation imaging      Residual anomaly      Physical model      Inversion     
1 引 言

重磁数据三维反演是通过将地下反演空间剖分成长方体单元，应用采集的重磁数据，采用线性或者非线性反演算法来反演出这些单元的物性参数，勾画出场源的分布图像，从而确定矿产资源的赋存状态或者是地质构造形态[1-2].近几年来，很多专家学者对重磁三维反演问题进行了深入的研究，研究了诸如线性反演中的大矩阵求解问题，病态方程的约束问题，非线性反演的计算速度和存储问题，推动了三维反演算法的不断进步和走向应用[3-10].

重磁三维相关成像技术计算的是实测异常与地下各点产生的异常的归一化互相关系数.该类方法最早由Patella用于自然电位异常的解释，称为概率成像方法[11-12]，Mauriello等[13-14]将改方法推广到大地电磁，重力领域，Iuliano 等[15-16]在任意标量和矢量的地球物理场概率成像基础上发展了结合重力和磁力，重力和自然电位等的联合概率成像计算，国内许多专家学者都对概率成像方法在自然电场和大地电磁数据上的应用进行了深入的研究，取得了很多认识和进展.

重力数据的三维相关成像由郭良辉等[17-18]根据概率成像的概念发展而来，并将之推广到磁法数据相关成像，经过模型和实际数据的实验证明了该方法在一定程度上能够勾画出重磁异常的场源分布，并且具有计算效率高，计算过程稳定，容易实现等特点.但是，该方法求取的只是界于-1 和+1 之间的等效物性参数，而不是真实的密度或者是磁化率，并且在反演计算的过程中也无法引入已知地质认识的约束，一定程度上限制了该方法的广泛应用.本文根据重磁正演和相关成像计算的线性性质，提出了应用观测数据与模型正演结果残差进行相关成像更新物性模型的迭代反演过程，反演出地质体的重磁物性参数，该方法同时继承了相关成像的计算简单、稳定、易于实现等特点.

2 重磁相关成像基本理论

假设测区地下任意坐标为(xq，yq，zq)处，体积为vq的第q个长方体质量的剩余密度为Δσq，则它在测区上任意点(x，y，z)处的重力异常Δgq(x，y，z)可表示为

(1)

其中万有引力常数G=6.667×10-11m3/(kg·s2).

定义测区上实测重力异常与第q个长方体质量在测区上的重力异常的归一化互相关Cq为

(2)

其中，(xi，yi，zi)表示测区第i个观测点的坐标，Δg(xi，yi，zi)为该点的实测重力异常，Δgq(xi，yi，zi)为地下第q个点质量在该观测点上的重力异常，n为观测点数.

将公式(1)代入公式(2)可得到如下的表达：

(3)

其中，Bq(xi，yi，zi)为第q个质量点对观测点(xi，yi，zi)的重力异常基函数，表达式为

(5)

根据同样的推理，实测的ΔT磁异常与第q个磁偶极子的磁异常归一化相关Cq为

(5)

其中，ΔT(xi，yi，zi)为(xi，yi，zi)处的磁异常总强度，n为观测点数，Bq(xi，yi，zi)为第q个偶极子对观测点(xi，yi，zi)的磁力异常基函数，其具体表达如下：

(6)

其中，r是观测点到q点的距离，A，B，C，D，E，F是一组与磁化率、磁倾角等有关的参数，具体表达可详见文献[15].

三维相关成像的离散实现步骤为：

(1) 将下半空剖分成与实测网格一致的三维均匀网格；

(2) 计算每一网格结点单位剩余物性差所产生的异常与实测异常在一定窗口范围内的归一化互相关(称场源发生的概率)；

(3) 将计算结果置于网格点；

(4) 逐点移动扫描点使其遍历下半空间所有的网格点；

(5) 根据网格结点上场源出现的概率勾画出场源的分布情况.

重力和磁异常的相关系数Cq可理解为物性参数的等效表达，它的大小介于-1 和+1 之间，并且值的大小与物性参数的大小成线性对应关系，值为正则表示物性参数相对背景场的盈余，且越大说明盈余越多，负值表示物性参数相对背景场的亏损，且负的越多说明亏损越多，因此，相关成像在一定程度上能够表征地质体物性参数以及地质结构的变化.从上述公式可以看出，在相关系数的求解过程中，无大矩阵求逆等矩阵计算，其计算简单、稳定、可并行性强.但是相关成像的结果只是等效物性参数，并且计算过程中不能够引入已知地质条件的约束.

3 基于剩余异常相关成像的重磁物性反演

相关成像过程具有无大规模矩阵计算，可并行性强，计算效率和计算简单，但是也只能获得等效物性参数，无法引入已知地质信息约束.本文提出利用相关系数与物性参数成对应关系的性质，在相关系数的基础上反演物性参数.具体的做法是给定初始物性参数模型，进行正演计算，并以实测的数据与正演结果的均方误差作为物性反演结束的标准，如果达到标准，则反演完成，如果没有达到标准，则对实测与正演计算数据的残差进行相关成像，以一个物性参数的变化量乘以相关成像结果作为模型的扰动量来更新物性参数模型，反复进行上述计算，最终收敛得到反演结果.这样的物性参数反演过程为基于剩余异常相关成像的重磁物性反演方法，其流程可见图 1.

基于剩余异常相关成像的重磁物性反演方法是根据相关成像系数与物性参数的对应关系发展而来，整个反演过程继承了相关成像计算的计算过程简单，计算效率高的特点，并且在反演的过程中可以引入已知先验信息，深度加权等约束条件，这相比单纯的相关成像计算来说有了很大的灵活性，但基于相关成像的物性参数反演方法依然具有多解性问题，反演得到的物性参数也会受初始模型和引入的地质约束的影响.

图 1 基于剩余异常相关成像的重磁物性反演流程 Fig. 1 Flowchart of imaging based on residual anomaly correlation
4 模型数据实验

为验证本文反演方法的有效性，设计了包含了两个异常磁化率的三维模型进行测试，异常体的长宽各200m，厚度为100 m;反演数据体的长宽各1200m，厚度为500 m，反演网格间距为10 m×10m×10 m，异常体的具体位置参数如表 1.图 2是模型正演获得的ΔT异常，单位为nT.

表 1 模型的异常体参数 Table 1 Parameter of the model
图 2 模型异常的等值线分布图 Fig. 2 Magnetic anomaly contour of the model

采用过两个异常体中心的线作质控，进行模型相关成像和本文反演的对比，为简化称这条线为A线(图 3).

从图 3相关成像的结果上基本能够刻画出两个磁异常区域，但整体的分辨率很差.利用模型相关成像的结果乘以0.01作为初始磁化率模型，应用本文方法对该模型进行物性参数反演，整个过程暂不设置任何约束.图 4是模型反演结果在A 处的剖面，其中反演过程中迭代次数为12次，反演残差为0.024.

从图 4的形态上看，基于本文方法的反演结果对比相关成像结果分辨率明显提高，基本上能够勾画出异常体的范围，并且得到了真正的磁化率参数.

本文方法是物性参数反演方法，相对于相关成像而言，能够引入已知的地质信息作为约束，例如在反演模型更新过程中，对已知地质信息的区域设置反演参数的最大值和最小值.对于本例的模型，如果在反演过程中设置异常体外的值在背景异常附近小范围波动，以图 4结果为初始模型，迭代6次后得到如图 5的反演结果.而这个结果也间接证明了引入模型约束对反演结果可靠性的提高至关重要.

图 3 A线的异常（a)，过A线的模型（b)和相关成像结果（c) Fig. 3 Aanomaly of profile A (a)，the model (b) and correlation imaging resul to fprofile A (c)
图 4 本文方法反演结果 (a)反演曲线和采集曲线；（b)过A线的反演的磁化率(SI)结果. Fig. 4 The result of inversion using this paper method (a) Inversion curve and observed curve，(b) Inversed magnetic susceptibility (SI) crossing A of the model.
图 5 加上已知信息的本文方法反演结果 (a)反演曲线和采集曲线；（b)过A线的反演剖面. Fig. 5 The result of inversion using this paper method with geology constrain(a) Inversion curve and observed curve，(b) Inversion result cross A of the model.
5 影响基于剩余异常相关成像的重磁物性反演的因素分析

本文介绍的反演方法具有计算简单，计算效率高等特点，但是同时他也具有其他反演方法所具有的多解性，对初始模型依赖性强等特点，从反演的流程上看，本文方法主要是依靠相关成像的结果乘以一个给定的物性参数的变化值来更新物性参数模型，因此物性参数变化值的选择对反演过程的收敛速度有很大的影响，本文建议在迭代计算的开始，选择一个相对较大的物性参数变化值，等即将收敛时，物性参数的变化值也同时减小.

6 结 论

本文提出一种基于剩余异常相关成像的重磁物性反演方法，该方法相对于单独的相关成像来说，提高分辨率的同时能够反演真实的物性参数，而不是等效参数，并且可以引入已知的地质先验信息的约束提高反演的分辨率.另外，相对传统的物性参数反演方法而言，没有大矩阵计算所带来的奇异性等问题，因此具有计算简单，计算效率高等特点.在深部矿产勘探和地质调查上有广泛的应用前景.

参考文献
[1] 郝天珧, 徐亚, 赵百民, 等. 南海磁性基底分布特征的地球物理研究. 地球物理学报 , 2009, 52(11): 2763–2774. Hao T Y, Xu Y, Zhao B M, et al. Geophysical research on distribution features of magnetic basements in the South China Sea. Chinese J. Geophys. (in Chinese) , 2009, 52(11): 2763-2774.
[2] 郝天珧, 黄松, 徐亚, 等. 关于黄海深部构造的地球物理认识. 地球物理学报 , 2010, 53(6): 1315–1326. Hao T Y, Huang S, Xu Y, et al. Geophysical understandings on deep structure in Yellow Sea. Chinese J. Geophys. (in Chinese) , 2010, 53(6): 1315-1326.
[3] 徐世浙, 曹洛华, 姚敬金. 重力异常三维反演-视密度成像方法技术的应用. 物探与化探 , 2007, 31(1): 25–28. Xu S Z, Cao L H, Yao J J. 3D inversion of gravity anomaly: an application of the apparent density imagery technology. Geophysical & Geochemical Exploration (in Chinese) , 2007, 31(1): 25-28.
[4] 姚长利, 郝天姚, 管志宁, 等. 重磁遗传算法三维反演中的高速计算及有效存储方法技术. 地球物理学报 , 2003, 46(2): 252–258. Yao C L, Hao T Y, Guan Z N, et al. High speed computation and efficient storage in 3-D gravity and magnetic inversion based on genetic algorithms. Chinese J. Geophys. (in Chinese) , 2003, 46(2): 252-258.
[5] 姚长利, 郑元满, 张聿文. 重磁异常三维物性反演随机子域法方法技术. 地球物理学报 , 2007, 50(5): 1576–1583. Yao C L, Zheng Y M, Zhang Y W. 3-D gravity and magnetic inversion for physical properties using stochastic subspaces. Chinese J. Geophys. (in Chinese) , 2007, 50(5): 1576-1583.
[6] 杨文采, 余长青. 根据地球物理资料分析大别—苏鲁超高压变质带演化的运动学与动力学. 地球物理学报 , 2001, 44(3): 346–359. Yang W C, Yu C Q. Kinetics and dynamics of development of the Dabie-Sulu UHPM Terranes based on geophysical evidences. Chinese J. Geophys. (in Chinese) , 2001, 44(3): 346-359.
[7] Barbosa V C F, Silva J B C. Generalized compact gravity inversion. Geophysics , 1994, 59(1): 57-68. DOI:10.1190/1.1443534
[8] Li Y G, Oldenburg D W. 3-D inversion of magnetic data. Geophysics , 1996, 61(2): 394-408. DOI:10.1190/1.1443968
[9] Li Y G, Oldenburg D W. 3-D inversion of gravity data. Geophysics , 1998, 63(1): 109-119. DOI:10.1190/1.1444302
[10] Silva J B C, Hohmann G W. Nolinear magnetic inversion using a random search method. Geophysics , 1983, 48(12): 1645-1658. DOI:10.1190/1.1441445
[11] Rene R M. Gravity inversion using open, reject, and "shape to the anomaly" fill criteria. Geophysics , 1986, 51(4): 988-994. DOI:10.1190/1.1442157
[12] Patella D. Introduction to ground surface self-potential tomography. Geophysical Prospecting , 1997, 45(4): 653-681. DOI:10.1046/j.1365-2478.1997.430277.x
[13] Mauriello P, Patella D. Principles of Probability tomography for natural-source electromagnetic induction fields. Geophysics , 1999, 64(5): 1404-1417.
[14] Mauriello P, Patella D. Localization of maximum-depth gravity anomaly sources by a distribution of equivalent point masses. Geophysics , 2001, 66(5): 1431-1437. DOI:10.1190/1.1487088
[15] 毛立峰, 王绪本, 高永才. 大地电磁概率成像的效果评价. 地球物理学报 , 2005, 48(2): 429–433. Mao L F, Wang X B, Gao Y C. Appraisement on the MT probability tomography. Chinese J. Geophys. (in Chinese) , 2005, 48(2): 429-433.
[16] 王绪本, 毛立峰, 高永才. 电磁导数场概率成像方法研究. 成都理工大学学报 (自然科学版) , 2004, 31(6): 679–684. Wang X B, Mao L F, Gao Y C. Probability tomography of electromagnetic field-derivative method. Journal of Chengdu University of Technology (Sciences & Technology Edition) (in Chinese) , 2004, 31(6): 679-684.
[17] 郭良辉, 孟小红, 石磊, 等. 重力和重力梯度数据三维相关成像. 地球物理学报 , 2009, 52(4): 1098–1106. Guo L H, Meng X H, Shi L, et al. 3-D correlation imaging for gravity and gravity gradiometry data. Chinese J. Geophys. (in Chinese) , 2009, 52(4): 1098-1106.
[18] 郭良辉, 孟小红, 石磊. 磁异常ΔT三维相关成像. 地球物理学报 , 2010, 53(2): 435–441. Guo L H, Meng X H, Shi L. 3D correlation imaging for magnetic anomaly ΔT data. Chinese J. Geophys. (in Chinese) , 2010, 53(2): 435-441.