基于预优共轭梯度的磁化强度反演方法及应用

曾琴琴，刘 双，王永华

(1. 成都地质矿床研究所，四川成都 610081；2. 中国地质大学地球物理与空间信息学院，湖北武汉 430074)

基于预优共轭梯度的磁化强度反演方法及应用

曾琴琴1，刘 双2，王永华1

(1. 成都地质矿床研究所，四川成都 610081；2. 中国地质大学地球物理与空间信息学院，湖北武汉 430074)

基于预优共轭梯度算法的反演方法在求解磁化强度大、规模欠定线性方程组时，可改善方程组的条件数，并提高算法的纵向分辨率。该方法对磁性体形状无要求，且收敛速度快、稳定性好，能较直观地反映磁性体特征。本文通过理论模型试验结果分析证明了基于预优共轭梯度的磁化强度反演方法的稳定可靠性，并对云南北衙铁金矿区万硐山矿段L32线磁异常进行了反演，结果表明了该方法应用于实际资料处理与解释的有效性和可靠性。

预优共轭梯度 深度加权预优因子 磁化强度反演 北衙铁金矿区

Zeng Qin-qin, Liu Shuang, Wang Yong-hua. The magnetization inversion method based on preconditioned conjugate gradient algorithm and its application[J].Geology and Exploration,2014,50(4):0756-0762.

1 引言

岩、矿石的密度与磁性是研究构造类型、岩石年代及类型的重要参数，也是位场异常解释、地质填图、矿产分布等研究的重要依据。由于野外采集标本有限，地下大量岩、矿石物性信息采集不全，因而利用观测的位场异常反演地下岩、矿体位置对位场资料的解释具有重要参考价值。依据测量异常数据，结合地质和其它资料，求解地下场源体的空间位置、形状及物性特征参数就是地球物理反演所要解决的问题(刘璎等，2011)。

磁化强度反演技术以大量不同磁化强度的网格剖分单元的组合来模拟地下复杂的磁性体，每个网格单元磁化强度与磁异常数据之间为线性关系，通过在空间域建立目标函数，利用比较成熟的线性反演理论来求解各个单元格的磁化强度值(Lietal.,1996)，以得到地下介质的磁性分布特征。从刻画地下复杂磁性体的灵活性和对观测数据的拟合能力来说，磁化强度反演技术比传统解释方法更具有优势(杨波，2010)。预优共轭梯度算法是近年来发展起来的一种快速而精确的计算方法，它具有迭代次数少、收敛迅速等优点(杨振威等，2012)，从而被引入到位场资料处理解释中。张世晖(2003)将预优共轭梯度算法与密度成像相结合，并采用小波压缩方法对建立的方程进行压缩；Zhdanovetal.,(2004)利用重加权正则化共轭梯度算法，采用聚焦正则化因子约束目标函数较好地反映了尖锐边界特征；Wilson等(2011)针对巨大的航空重力梯度测量数据，利用共轭梯度法反演，并结合并行算法提高了算法的计算效率。

本文将深度加权预优因子引入到共轭梯度反演中，使得预优共轭梯度矩阵核函数的特征值集中分布，由此提高了磁化强度反演算法的收敛速度和纵向分辨率，并通过理论模型试验和北衙铁金矿区L32线解释结果证明了该方法的有效性和可靠性。

2 方法原理

2.1 二度矩形截面棱柱正演

如图1所示，斜磁化二度板状体模型沿走向无限延长，中心点坐标为(x0,z0)，宽度为2b，延深长度为2l，倾角为α，有效磁化强度为Ms，磁化倾角为i。设：

图1 二度板状体坐标示意图Fig. 1 Sketch of coordinate of a 2D tabular body

则该板状体在P(x,z)点产生的磁异常可表示为：

(1)

(2)

上式中，Ha,Za分别表示磁场强度的水平和垂直分量。若地磁场的倾角和测线的磁方位角分别为I和A，则ΔT可表示为(管志宁，2005)：

(3)

2.2 预优共轭梯度反演方法

将地下空间划分为N个致密排列大小固定的块体单元，可以建立磁化强度参数与磁测数据的线性方程组，通过正则化方法反演光滑解(刘天佑，2007；刘展等，1999)。若测线上有M个观测点，则第j个块体单元在第i个观测点的磁异常为

式中，mj表示第j个棱柱体单元的磁化强度(或磁化率)；Gij表示单位大小的磁化强度的第j个块体单元在第i个观测点产生的异常。根据位场的叠加原理，第i个观测点的磁异常ΔTi是不同深度的块体单元在该点处的磁异常之和(Boulangeretal.,2001; Oldenburgetal.,1994)，即

即

(4)

式中，ΔT为M×1维向量，表示M个观测点上的磁异常；m为N×1维向量，表示N个块体单元的磁化强度；G为M×N维矩阵，称为核矩阵。通常，观测数据的个数远远小于地下模型参数的个数，即M<

GTGm=GTΔ

即

Am=b

(5)

式中，A为N阶对称正定矩阵，b为N维向量。而求解以上方程组等价于求解二次泛函的极小点。

(6)

因此，可以利用泛函φ(m)的二次终止性，沿共轭梯度方向搜索可快速的求解，这种方法即为共轭梯度法。目前，共轭梯度法已经成为求解大型线性方程组最有效的方法之一，在地球物理反演中得到广泛的应用(吴小平等，1998，2000；朱培民等，2000；李冰等，2002；翁爱华，2007；翁爱华等，2012)。

为避免迭代次数太大，提高共轭梯度法的收敛速度，VanDecar等(1994)提出采用预优矩阵改善方程组条件数，并导出求解类似方程组的预优共轭梯度算法，即

PATAm=PATΔb

(7)

式中，P为预优矩阵，近似等于(ATA)-1，即P(ATA)≈I(I为单位阵)。P(ATA)的奇异值集中分布在对角线上，且非常接近1，由此改善了方程组的条件数。

由于构造模型的核函数是线性的，位场异常与场源到观测点的距离呈指数衰减，导致核矩阵中数值随深度增加而急剧减小，相同的棱柱体单元，深部的异常响应比浅部的要弱得多，容易出现“趋肤效应”。Li and Oldenburg(1996)在3D磁化率成像反演中将模型范数修改为目标函数的深度加权函数，随后又将磁化率反演算法应用于密度反演(Lietal.,1998)。Pilkington引入了简单可靠的预优因子(Pilkington,1997)，即P=z3I，z为块体单元中心埋深，I为单位阵，本文即引用该预优因子提高算法的纵向分辨率。

3 理论模型试验

理论模型由有限延深板状体组成，已知测线方位角为90°，地磁倾角为45°，图2为不同板状体的模型反演结果对比。其中，图中上部曲线为模型正演磁异常；中部红色块体为理论模型体，已知磁化强度为10000×10-3A/m；下部为利用预优共轭梯度反演方法得到的磁性体。可见，基于预优共轭梯度的反演方法反映的不同埋深、产状的磁性体特征结果与理论模型较为一致。但该方法是一种光滑成像，其反映的磁性体边界要比实际模型更圆滑，且随磁性体埋深增大，反演效果逐渐减弱。

图2 理论模型磁化强度反演结果(A=90°，I=45°)Fig.2 Inversion results of theoretical models(A=90°,I=45°) a-模型试验1；b-模型试验2；c-模型试验3 a-model test 1; b-model test 2; c-model test 3

4 应用实例

北衙铁金矿区位于云南省大理州鹤庆县北衙乡，大地构造上处于扬子板块的西部边缘(葛良胜等，2002)。区内发育古生界志留系、泥盆系及二叠系，中生界三叠系，新生界始新统及上新统和第四纪地层。以三叠系最为发育，早古生代地层仅见于研究区的西北角，且零星出露。其南部及东南部大面积出露峨眉山玄武岩，另见大面积的古近纪斑岩(石英斑岩、正长斑岩)侵入到峨眉山玄武岩和三叠纪地层中，沿接触带矽卡岩化、大理岩化、角岩化、褐铁矿化发育，并有少量的基性-超基性脉岩侵入到三叠系中。矿区铁金矿体围岩为三叠系中统北衙组不纯碳酸盐岩，该套地层中矿化元素含量普遍较高，受岩浆后期含矿热液作用，地层普遍出露热变质及交代蚀变形成大理岩、矽卡岩及含金的磁铁矿-硫化物矿体(原生矿石)。矿区岩(矿)石磁性统计结果(见表1)表明，区内铁矿石、含金铁矿石磁性最强，与围岩及其它岩石磁性差异较大，磁法勘探具有前提条件。

图3为北衙铁金矿区磁异常平面图。研究区东部玄武岩广泛出露地带表现为一系列沿南北向断续相接的条带状磁异常；西部磁异常主要由万硐山和红泥塘两个矿段组成，其中，南部红泥塘矿段磁异常规模大，呈宽缓的弧形异常带延伸至金钩坝以南；北部万硐山矿段磁异常规模小、幅值大，呈规则的南正北负等轴状椭圆形，表现为浅埋磁性体异常特征。

表1 北衙矿区岩(矿)石磁性统计表(据尹潘等，2005)

图3 北衙铁金矿区磁异常平面图Fig. 3 Map showing magnetic anomalies of the Beiya iron-gold deposit

图4为万硐山矿段L32线综合解释结果，其中，图(a)为剖面磁异常曲线，(b)为基于共轭梯度的反演方法得到的磁性体分布特征，(c)为对应的地质勘探剖面。由图(c)可知，在剖面内200 m深处，钻孔32ZK20钻遇铁金矿体KT4中心，其埋深约100 m，厚度达30 m，矿体沿两端逐渐减薄；图(b)显示，剖面内主磁性体分布于中部100～300 m段，磁化强度≥3300×10-3A/m的地质体埋深约85 m，中心厚度达42 m，呈顶部宽、底部窄的“▽”形，与矿体KT4特征基本吻合，结合区内磁性统计结果，推测该磁性体主要反映了矿体KT4的平均效应。受矿体厚度和磁性分布的不均匀分布影响，局部位置反演效果减弱。

图4 万硐山矿段L32勘探线反演与解释结果Fig. 4 Inversion and interpretation results of exploration line L32 in the Wandongshan ore section a-剖面磁异常；b-磁化强度反演结果；c-地质勘探剖面 a-profile magnetic anomaly; b-magnetization inversion result; c-geological exploration profile

5 结论

基于预优共轭梯度算法的磁化强度反演方法收敛速度快，耗时短、稳定性好，且算法直接从目标函数出发，无需形成显式线性方程组，运算成本低。该算法对磁性体个数和形状无要求，在求解地球物理大规模反演问题中具有优势，尤其在研究浅埋深磁性体时，能够提供较为直观可靠的参考资料。

北衙铁金矿区万硐山矿段L32线解释结果表明，基于预优共轭梯度法的磁化强度反演方法解释结果具有较高的可靠性，该方法可广泛应用于万硐山矿段的磁异常分析与解释，进一步指导该矿段矿产勘探与开发工作。

Boulanger O., Chouteau M. 2001. Constraints in 3D gravity inversion[J]. Geophysics Prospecting, 49:265-280

Ge Liang-sheng, Guo Xiao-dong, Zou Yi-lin, Li Zhen-hua, Zhang Xiao-hui. 2002. Geological characteristics and genesis of Beiya gold depsit, Yunnan Province[J]. Contributions to Geology and Mineral Resources Research, 17(1):32-41(in Chinese with English abstract)

Guan Zhi-ning. 2005. Geomagnetic field and magnetic exploration[M]. Beijing: Geological Publishing House: 109-112(in Chinese)

Li Bing, Liu Hong, Li You-ming. 2002. 3-D seismic data discrete smooth interpolation using conjugate gradient[J]. Chinese J. Geophys, 45(5): 691-699(in Chinese with English abstract)

Li Y, Oldenburg D. W. 1996. 3-D inversion of magnetic data [J]. Geophysics,61(2):394-408

Li Y, Oldenburg D.W.1998. 3-D inversion of gravity data[J]. Geophysics, 63:109-119

Liu Ying, Meng Gui-xiang, Yan Jia-yong, Deng Gang.2011.Application of 3D property inversion for gravity and magnetic data to metal mineral exploration[J].Geology and Exploration,47(3): 448-455(in Chinese with English abstract)

Liu Tian-you. 2007. New data processing methods for potential field exploration[M]. Beijing: Science Press:88-113(in Chinese)

Liu Zhan, Wang Wan-yin. 1999. A method of regularized linear programming applying to directly inversion of potential fields on a curved surface[J]. Geology and Exploration,35(6):62-66(in Chinese with English abstract)

Oldenburg D. W., Li Y.G. 1994. Subspace linear inversion methods[J].Inverse Problems,10: 915-935

Pilkington M. 1997. 3-D magnetic imaging using conjugate gradients[J]. Geophysics,62:1132-1142

Vandecar J.C., Snieder R. 1994. Obtaining smooth solutions to large linear inverse problems[J]. Geophysics,59:818-829

Weng Ai-hua. 2007. Occam’s inversion and its application to transient electromagnetic method[J]. Geology and Exploration,42(5):74-76(in Chinese with English abstract)

Weng Ai-hua, Liu Yun-he, Jia Ding-yu, Liao Xiang-dong, Yin Chang-chun. 2012. Three-dimensional controlled source electromagnetic inversion using non-linear conjugate gradients[J]. Chinese J. Geophys, 55(10):3506-3515(in Chinese with English abstract)

Wilson G, Cuma M, Zhdanov M S. 2011. Massively parallel 3D inversion of gravity and gravity gradiometry data[J]. Preview, 152(6): 29-34

Wu Xiao-ping, Xu Guo-ming. 1998. Occam inversion improveness of magnetotelluria data[J].Chinese J.Geophys, 41(4): 547-554(in Chinese with English abstract)

Wu Xiao-ping, Xu Guo-ming. 2000. Study on 3-D resistivity inversion using conjugate gradient method[J]. Chinese J. Geophys, 43(3):420-427(in Chinese with English abstract)

Yang Bo. 2010. 3D smooth inversion of magnetic data and its application to the Daye iron mine[J]. Progress in Geophys,25(4):1433-1441 (in Chinese with English abstract)

Yang Zhen-wei, Yan Jia-yong, Liu Yan, Wang Hua-feng. 2012. Research progresses of the high-density resistivity method[J]. Geology and Exploration,48(5): 969-978 (in Chinese with English abstract)

Yin Pan, Song Li-jun. 2005. Alkali rich porphyry metallogenic regularity research and evaluation Beiya gold deposit in Yunnan[J]. Journal of nonferrous metal(Mine section),57(3):22-24(in Chinese with English abstract)

Zhang Shi-hui. 2003. Marine gravity data from Satellite altimetry processing method research and its application in Okinawa trough[D].Wuhan: China University of Geosciences:1-60(in Chinese)

Zhdanov M S, Ellis R, Mukherjee S. 2004.Three-dimensional regularized focusing inversion of gravity gradient tensor component data[J]. Geophysics, 69(4): 925-937

Zhu Pei-ming, Wang Jia-ying, Zhan Zhen-gbin, Gu Han-ming, Zhu Guang-ming. 2000. Stochasitic conjugate gradient inversion[J]. Oil Geophysical Prospecting, 35(2):208-213(in Chinese with English abstract)

[附中文参考文献]

葛良胜，郭晓东，邹依林，李振华，张晓辉．2002．云南北衙金矿床地质特征及成因研究[J]．地质找矿论丛，17(1)：32-41

管志宁．2005．地磁场与磁力勘探[M]．北京：地质出版社：109-112

李冰，刘洪，李幼铭．2002．三维地震数据离散光滑插值的共轭梯度法[J]．地球物理学报，45(5)：691-699

刘璎，孟贵祥，严加永，邓刚．2011．重磁3D物性反演技术在金属勘探中的应用[J]．地质与勘探，47(3)：448-455

刘天佑．2007．位场勘探数据处理新方法[M]．北京：科学出版社：88-113

刘展，王万银．1999．曲面位场的正则化线性规划法直接反演技术[J]．地质与勘探，35(6)：62-66

翁爱华．2007．Occam反演及其在瞬变电磁测深中的应用[J]．地质与勘探，42(5)：74-76

翁爱华，刘云鹤，贾定宇，廖祥东，殷长春．2012．地面可控源频率测深非线性共轭梯度反演[J]．地球物理学报，55(10)：3506-3515

吴小平，徐果明．1998.大地电磁数据的Occam反演改进[J]．地球物理学报，41(4)：547-554

吴小平，徐果明．2000．利用共轭梯度法的电阻率三维反演研究[J]．地球物理学报，43(3)：420-427

杨波．2010．三维光滑磁化率成像反演及在大冶铁矿的应用[J]．地球物理学进展，25(4)：1433-1441

杨振威，严加永，刘彦，王华峰．2012．高密度电阻率法研究进展[J]．地质与勘探，48(5)：969-978

尹潘，宋立军．2005．云南北衙金矿富碱斑岩成矿规律研究及评价[J].有色金属(矿山部分)，57(3)：22-24

张世晖．2003．海洋卫星测高重力数据处理方法研究及在冲绳海槽的应用[D].武汉:中国地质大学:1-60

朱培民，王家映，詹正彬，顾汉明，朱光明．2000．随机共轭梯度反演法[J]．石油地球物理勘探，35(2)：208-213

The Magnetization Inversion Method Based on PreconditionedConjugate Gradient Algorithm and its Application

ZENG Qin-qin1, LIU Shuang2, WANG Yong-hua1

(1.ChengduInstituteofGeologyandMineralResources,Chengdu,Sichuan610081;2.InstituteofGeophysics&Geomatics,ChinaUniversityofGeosciences,Wuhan,Hubei430074)

The inversion method based on preconditioned conjugate gradient algorithm is used to solve large-scale underdetermined linear magnetization equations, which can improve the condition number of equations and vertical resolution of the algorithm. This method not only has no requirements on magnetic body shape, but also has a fast convergence speed and good stability, with which the characteristics of magnetic bodies can be reflected more intuitively. The tests on theoretical models prove this method is stable, reliable and fast. And its use to and the inversion of real magnetic data of exploration line L32 in the Wandongshan ore section of the Beiya iron-gold deposit also indicates the effectiveness and reliability of this method.

preconditioned conjugate gradient, depth weighted preconditioned factor, magnetization inversion, Beiya iron-gold deposit

2013-12-27；

2014-05-31；[责任编辑]郝情情。

中国地质调查局地质调查工作项目(编号1212011120918、1212011220249和12120113051500)联合资助。

曾琴琴(1982年-)，女，2011年毕业于中国地质大学(武汉)，获博士学位，工程师，长期从事地球物理数据处理解释及方法技术的研究工作。E-mail：zqqcgs@163.com。

P631

A

0495-5331(2014)04-0756-7