本科毕业论文-—磁化强度成像及其在航磁资料解释中的应用

桂林理工大学硕士研究生学位论文分类号： 密级： 编号：212011479 桂林理工大学 硕 士 研 究 生 学 位 论 文磁化强度成像及其在航磁资料解释中的应用专 业： 地质工程研究方向： 应用地球物理研 究 生： 周海滨指导教师： 王有学 教授论文起止日期：2012年4月至2013年4月Magnetic imaging and its application in Aero-magnetic Data interpretationMajor：Geological EngineeringDirection of Study：Application of geophysical Graduate Student： Haibin ZhouSupervisor： Prof. Youxue WangCollege of Earth ScienceGuilin University of TechnologyApril, 2012 to April, 2013 研究生学位论文独创性声明和版权使用授权书独 创 性 声 明本人声明：所呈交的论文是本人在导师指导下进行的研究工作及取得的研究成果。据我所知，除了文中特别加以标注和致谢的地方外，论文中不包含他人已经发表或撰写过的研究成果，也不包含为获得其它教育机构的学位或证书而使用过的材料。对论文的完成提供过帮助的有关人员已在论文中作了明确的说明并表示谢意。 学位论文作者（签字）： 签字日期： 学位论文版权使用授权书 本学位论文作者完全了解(学校)有关保留、使用学位论文的规定，有权保留并向国家有关部门或机构送交论文的印刷本和电子版本，允许论文被查阅和借阅。本人授权(学校)可以将学位论文的全部或部分内容编入有关数据库进行检索，可以采用影印、缩印或扫描等复制手段保存、汇编学位论文。同时授权中国科学技术信息研究所将本学位论文收录到中国学位论文全文数据库，并通过网络向社会公众提供信息服务。（保密的学位论文在解密后适用本授权书）本论文是否保密： 是 否 如需保密，保密期限为：学位论文作者签名： 导师签字：签字日期： 年 月 日 签字日期： 年 月 日摘 要地下各种具有不同磁性差异的岩（矿）石会在其周围产生不同的磁性异常场，磁力勘探是通过分析解释观测所得到的磁异常数据，用以勘察有用矿产及寻找地下地质构造的一种较为常用地球物理方法，具有广泛的应用。磁测数据的采集简单易行、数据采集效率高，且不受地域的影响。磁异常反映了地下磁性特征的变化，包括磁性异常体的分布和磁化强度。通过有效的反演方法，结合已知的地质资料，就可以对磁性异常体的大小、形状、磁化强度等信息进行估算。本文在不考虑磁性异常体的剩磁及退磁影响的情况下，对二度磁性异常体进行了磁化强度成像研究。在实际工作中，由于二度体和三度体在剖面上可能有相同的异常形态，而本文所介绍的方法中，磁异常数据为一维(剖面)数据，所有地质体都当作二度体的来解释，部分解释结果与实际情况不一致是在所难免的，若要进行更为细致和准确的解释，那么建议采用平面异常数据，而对于位场反演问题来说，三维反演与二维反演相比，并没有本质上的差异，实现三维反演仅需要做细节上的调整便可，反演方法并没有太大的改变。在研究了二度磁异常体正演公式的基础上，本文给出了二维磁化强度成像的基本公式。在构造目标函数时，基于Occams反演方法，在目标函数中加入了光滑度项及最小模型项来对模型进行约束，同时为了克服反演结果富集在地表附近的现象，采用深度加权函数对空间单元赋予不同的权重，得到比较符合实际的结果。本文还讨论了如何选择最优的正则化参数，通过选取合理的正则化参数，使反演的结果较为合理。在求解方程组时，采用了计算速度快、节省内存且收敛速度快的共轭梯度算法。通过对多组理论模型的正演计算及其磁化强度成像结果的对比分析，讨论了反演算法的分辨率和可靠性，并对在反演中所选择的参数进行了总结和归纳。本文的反演算法不依赖于初始模型设置，成像效果较好。将该反演方法应用于航磁剖面，反演的结果为构造的解释提供了一定依据。关键词：二维光滑磁化强度成像；Occams反演；深度加权AbstractA variety of underground rock (ore) with different magnetic differences will produce different magnetic anomaly field. Magnetic exploration is a commonly used geophysical method which through the analysis of magnetic anomaly data to investigate valuable mineral and find underground geological formations with a wide range of applications. Collecting magnetic data is simple, high efficient and without geographical impact. The magnetic anomaly reflects the changes of the underground magnetic characteristics, including the distribution and magnetization of the magnetic anomalies. Through effective inversion method which combined with known geological data, we can estimate the information about the size, shape and the intensity of magnetization of magnetic anomaly.Inthis paper, without considering about the remanent magnetization and demagnetization. In actual work, 2D-bodys magnetic anomaly and 3D-bodys magnetic anomaly usually have the same abnormal morphology, but the method described in this article, magnetic anomaly data is one-dimensional (profile) data，all the Geological body will explain as the 2D-bodys magnetic anomaly, therefore, partial interpretation results will different with the actual condition can hardly be avoided.Based upon the research of 2D-bodys magnetic anomaly forward formulas, this paper presents the magnetization imaging inversion formulas of 2D-bodys magnetic anomaly. When constructing the objective function, based on the Occams inversion method, theoretical aspects in 2D inversion of magnetic data are mentioned, such as the addition of smoothness matrix will be used to constraint the model. At the same time, I formulate the inversion using Tikhonov regularization. Further more, the experiments of some theoretical model indicated that the inversion method seemed to be reliable and practical. In this paper, I apply the conjugate gradient (CG) method to solve multi-dimensional linear equation. By comparing the results of the forward calculation of the theoretical model and its magnetization inversion, it will discuss the resolution and reliability of the inversion algorithm, and summarize how to choose the parameters in the inversion. The inversion algorithm does not depend on the initial model set. The method is applied to the inversion of aeromagnetic profile; the inversion results will provide the foundations of structural interpretation.Keywords：2D smooth magnetization imaging; Occam inversion; depth weighted目 录摘 要IAbstractII目 录III第1章 绪论11.1 问题的提出及研究背景11.1.1 研究背景11.1.2 问题的提出21.2 研究现状21.3 本文研究内容3第2章 磁场理论基础52.1 地磁场52.1.1 磁位52.1.2 磁位的微分方程72.1.3 二维直立板状体正演公式82.2 岩石的岩石磁性特征及场源磁模型82.3 磁异常11第3章 二维介质中的磁化强度成像123.1 引言123.2 2D磁化强度成像的基本原理及正演计算143.2.1 模型构建及正演计算143.2.2 模型网格剖分问题163.3 磁化强度成像的反演技术173.3.1 共轭梯度法173.3.2 反演模型设计183.3.3 无约束的剖面反演计算193.3.4 深度加权函数203.3.5 物性绝对约束233.3.5 光滑磁化强度成像反演方法26第4章 数值模拟314.1 模型试算一324.2 模型试算二334.3 模型试算三344.5 模型试算四36第5章 二维磁化强度成像在航磁数据解释中的应用415.1 位置概况425.2 航磁资料与磁化强度成像435.3 航磁异常解释44第6章 结论与建议476.1 结论476.2 几点建议47参考文献48个人简历51致谢52IV第1章 绪论1.1 问题的提出及研究背景1.1.1 研究背景磁力勘探通过研究地下磁异常体在地表产生的磁异常数据来分析研究磁异常体的空间位置、大小、形状及产状等情况，从而达到找矿和解决地质问题的一种地球物理勘探方法。其具有轻便、快速、效率高及效果明显等特点，因此被广泛的应用于大面积的扫面工作，是一种必不可少的地球物理勘探手段1,2。随着地球科学技术的发展，人类对地球浅层沉积建造和性质的研究已经不能满足生产、经济和社会发展，特别是可持续发展的要求。并且随着人类对地球的研究不断地深入，浅层地壳已经无法满足现今的研究需求，这迫使地球科学研究向纵深发展成为一个趋势，因为岩石圈是人类居住与获取各种资源、能源、改造和利用场所。在二十一世纪，地质构造作为岩石圈的重要组成部分与地球物理学息息相关。自二十世纪八十年代初以来，我国对大陆岩石圈的研究开始引起了地学界的广泛关注，国家也成立了相关部门3，实施了一系列的地学断面的研究工作，采用了多种地球物理方法，开展了对中国大陆岩石圈深部结构的研究，并产生了广泛的国际影响。在众多的地球物理方法中，唯有航磁调查覆盖了整个中国大陆及其海域，从而成为我国地学研究的一个亮点，吸引了众多国内外科学家的关注，为我国的地学研究提供了不可多得的地球物理资料。航磁异常与中国大陆的大地构造格局之间的存在着一定的内在联系，因此通过探讨航磁异常的反演结果，用以研究地壳深部构造。磁测数据资料中包含着地下丰富的信息，但目前的数据解析并没有充分利用好这些信息，特别是在对磁测数据的反演上。反演是磁测资料定量解释中的重要环节之一， 和其他地球物理方法一样，磁测资料解释的目的在于通过地面或航空等实测数据利用某种手段推断出地下磁化率（磁化强度）分布规律，从而达到寻找目标地质体的目的4-9。在这样的背景下，本文对磁测异常进行了磁化强度成像研究，通过采用有效的约束方法减少反演结果的多解性。当不考虑剩磁时，磁化强度与磁化率是正比例关系，比例系数仅与地磁场强度有关，因此，本文对磁化强度进行反演而非磁化率。1.1.2 问题的提出面对日益复杂的地质问题及高质量的磁测数据，传统的定量、半定量的反演并不能满足勘探需求。如何能更加有效的对实测磁数据进行定量反演解释，这是本文所要解决的难题。当前人们对磁化强度成像技术还有抱有一定的怀疑态度，其原因是位场反演固有的多解性，影响了反演的结果。磁化强度成像是一种通过将地表观测到的磁异常数据直接反演地下场源的方法，其工作方式为：将地下场源一定范围的区域剖分成为一系列规则组合单元体通常为小立方体或长方体，通过反演算法获得这些单元的磁化强度大小，由这些网格的强化强度变化推测场源的分布范围、形态特征等。所以，磁化强度成像结果能比较直观地重构出磁性异常体及构造的分布情况。在反演过程中，虽然观测值与理论模型产生的理论值很容易拟合，得到一个数据拟合空间的解模型，但是由于磁场观测数据的有限性及不准确性等因素的存在，导致反演结果具有非唯一性，以至于难以得到一个与实际情况相符合的解。解的非唯一性是对于所有反问题来说都是普遍存在的，这种非唯一性与反演方法、数据处理的好坏无关，这是地球物理反演所要解决的难题。磁化强度成像属于不适定的问题，对于不适定问题应当采用正则化方法解决，或者说在对解的估计时应当在使结果获得较高的分辨率与拟合方差最小之间取折衷。这就需要在反演的过程中利用正则化因子添加一些约束条件及先验信息，对模型进行限制。如何得到一个相对可靠的解，本文针对这个问题作了详细的研究。1.2 研究现状 对于位场层析成像方法的研究开始于二十世纪七十年代10-13,国内在八十年代也开始了相关研究14-16，九十年代后，对位场层析成像的研究已经取得了较大的进展，不再仅限于理论上的研究，在处理实际数据中也得到了应用17-21,23，并取得不错的效果。对于磁化率（磁化强度）成像问题的研究，Yaoguo Li & Oldenberg20，21，给出了基于深度加权的3D磁化率和密度成像方法，其在核函数中采用的深度加权的方法，对于提高纵向分辨率及恢复异常体的真实位置都具有很好的理论及实际物理意义，同时他们还提出了最光滑反演方法，加入了模型梯度模量最小的约束条件，反演结果在边缘有一定的过渡，又能够很好地反映理论模型的位置，比较符合实际情况；为了提高纵向上的分辨率，Yaoguo Li & Oldenberg22，将井中磁测数据和地面磁测数据进行了联合反演，取得了一定的效果；由于3D成像需要求解规模很大的线性方程组，为了减少对计算时间和存储空间的需求Pilkington23利用预优共轭梯度法进行3D磁化率成像，得到了比较理想的结果；在地球物理反演的过程中，多解性问题是不可避免的，为了减少多解性，Yaoguo Li & Oldenberg24、Oleg Portniaguine 和 Micheal Zhdanov26引入了数据压缩技术来减小计算量、加快反演速度；Oleg Portniaguine 和 Micheal Zhdanov25，26提出了聚焦反演算法，加入了异常边界锐化的约束条件；Pilkington27引入柯西范数来度量模型目标函数，从而使得到的解能够用最少的非零值来拟合观测数据；Constable28将Occam反演方法应用于大地电磁反演中求解多层地球物理模型的光滑解。 国内对于磁化率（磁化强度）成像的研究也取得了一些成果。郭良辉，孟小红等29提出一种三维相关成像方法并应用于磁异常三维成像，具有良好的成像效果；杨波30研究了三维光滑磁化率成像的方法，并将磁化率成像反演应用于大冶铁矿的反演解释中，反演效果良好；刘圣博31对具有强剩磁的磁性异常体三维磁化率成像进行了研究，采用磁异常模量进行成像反演，磁异常模量反演减弱了剩磁对反演结果的影响；王若等32利用磁性标量反演方法对航磁数据进行了反演成像；姚长利,郝天珧等33、李军,李才明等34对重磁反演存在的问题以及解决的方法进行了分析讨论，并且就重磁的三维物性反演方法提出了一些方向性的意见；姚长利等35就重磁三维物性反演中的快速正演计算问题进行了研究；姚长利，郑元满等35，36采用遗传算法对三维物性进行了反演，并对其中的计算速度和有效存储提出了一些方法；姚长利等37三维物性反演中的随机子域方法进行了详细的分析；王妙月等38利用磁化强度矢量进行了磁化强度层析成像的研究，通过模型反演，证明了该方法的可行性。磁化强度成像问题是一个不适定问题，必须进行约束反演。本文详细研究了磁化强度成像的方法，并亚东-格尔木地学断面的航磁数据进行了磁化强度成像，通过对成像结果的分析，用以了解地下的构造情况，对研究地壳深部情况有一定的意义。1.3 本文研究内容本文研究二维的磁化强度成像的方法。首先将观测区域地下空间剖分为规则的网格单元，网格的大小固定，采用光滑磁化强度成像算法求解网格单元的磁化强度，从而求出地下磁性体的分布情况，通过研究多组反演模型的磁化强度成像的规律和特点，最后通过实例反演研究磁性体的分布与大地构造的关系。在磁化强度成像问题中，目标函数是欠定的，没有唯一的解，为了求得反问题的一个相对准确的解，必须对目标方程进行约束，本文对如何获得唯一解作了相关研究。反演中，必须拟合数据，否则实测数据将失去意义，但仅仅拟合实测数据并不能获得唯一解，因而通过正则化方法引入模型目标函数进行约束，使得在无穷多个解中，找到一个既能满足数据拟合目标函数也能满足模型目标函数的解。模型目标函数的作用是降低反演中解的自由度，从而获得唯一解。因此，在建立的目标方程中包含了数据拟合函数目标函数和模型目标函数。在目标函数中通过选择合适的正则化参数权衡数据拟合目标函数和模型目标函数之间的权重，使反演结果既不过分拟合观测数据，也不过多地偏向模型目标函数，导致模型目标函数的约束空间变小，数据拟合不够，丢失观测数据中的有用信号。通常，在磁化强度成像的过程中，其成像结果会出现明显的地表富集现象，这不符合实际，针对这个问题，目标函数中加入了深度加权约束，使得成像重构的结果处于合适的位置。在传统的反演中，大多数都是采用了物性范围的绝对约束，在磁化强度成像中也采用了绝对约束，在一个合理且固定的磁化强度范围内取值，并获得了理想的反演结果。本文采用收敛快、效率高的共轭梯度法对目标方程进行求解。第2章 磁场理论基础2.1 地磁场地球外部的地磁场是由来源于地球内部和地球大气层中的不同的磁场组成，在地表可以观测到地磁场的各个分量 39、40，通常主要将地磁场分为两个部分，一是来源于地球内部的稳定磁场，这是地磁场中的主要部分，二是由地球外部引起的变化磁场，其大小比稳定磁场弱得多，由于其影响，使地磁场产生短周期的变化。本文主要研究的是来源于固体地球内部的稳定磁场，一般认为磁测数据经过日变矫正之后便是稳定磁场。 图 2.1 地磁三要素示意图地面任意点的地磁场总强度矢量可用三个参数对其进行描述：地磁场的倾角,地磁场的偏角,以及地磁场的大小,地磁场三要素如图2.1所示：为地磁场总强度矢量与水平面的夹角；为磁北方向与地理北之间的夹角；为地磁场场值的大小，与方向无关。本章主要研究磁场的基本理论，包括磁位、磁场的积分表达式及微分表达式、泊松关系式等，这些基本理论公式都是磁性体正演的基础。2.1.1 磁位关于地磁场基本公式的推导，一般认为地下介质是无传导电流的，因此对于地磁场的而言，只研究无电流条件下的稳定磁场。在磁性体的外部空间，稳定磁场是保守场（即），所以稳定磁场是一个无旋场，因为标量场的梯度必为无旋场，所以无旋场可以表为一个标量场的梯度，即可引入一个标量磁位，使得： (2.1)这里,为磁场强度；为磁位；式(2.1)中负号表示磁场强度的方向始终指向磁位减小最快的方向。磁场强度的分量为 (2.2)在无磁性介质的空间中， (2.3)式(2.3)中，为磁感应强度，为真空磁导率，综合以上分析，在直角坐标轴上的分量分别为 (2.4)所以，磁位方向导数的负数就是磁性体的磁场，通常用、分别表示磁场的垂直分量、沿轴的水平分量、沿轴的水平分量；用表示磁性体总磁场模，它与磁场的三个分量的关系为：；如图2.2所示，一个任意形状、体积为、密度为而磁化强度为的物体，取地面上点为坐标原点，轴垂直向下，轴在水平面内。位于点的体积元（可视为磁偶极子）在空间中任意点的磁位可表示为 (2.5) 图2.2 任意磁性异常体在地表引起的磁异常整个物体的磁位应为 (2.6)当为常矢量（均匀磁化）时，可移至积分号之外，根据互易定理，上式可写成 (2.7)对于任何三度磁性体产生的磁场强度都可以根据(2.7)求得其磁位的值，然后代入式(2.4)即可求出。对于二度体而言，只在平面内考虑引力位，与方向无关。2.1.2 磁位的微分方程已知磁介质中的稳定磁场与和的有如下关系 (2.8)式中，为感应磁化强度，为剩余磁化强度。在实际中，剩余磁化强度一般比感应磁化强度小得多，因此常被忽略不计。因为，所以上式变为： (2.9)式(2.9)中，而，、分别为介质的磁导率和相对磁导率。稳定磁感应强度是无散场，即 (2.10)上式可解得 (2.11)式中，分别为磁荷密度、感应磁荷密度和剩余磁荷密度。将(2.1)式代入(2.11)式中，可得 (2.12)上式为无传导电流磁介质中稳定磁场的磁位微分方程，即泊松方程。在磁介质之外或在均匀磁介质中，由于（即），磁位满足拉普拉斯方程 (2.13)在无传导电流的磁介质中，稳定磁场是有旋场，即 (2.14)在磁介质之外或在均匀磁介质中，由于，故上式变为0 (2.15)由以上讨论可知，在无传导电流的非磁性介质中（如空气或水中）或在均匀磁介质内，稳定磁场（和）均是无旋无散场。该结论是将稳定磁场与统一起来研究的基础。2.1.3 二维直立板状体正演公式本文网格划分的基本单元是矩形板状体，地表观测到的磁异常是由各个网格单元在观测点所产生的磁异常的叠加。我们假定矩形板状体单元均匀磁化，忽略剩磁及退磁的影响，只有磁荷分布，并且磁荷密度在同一磁荷面分布均匀，则单个板状体单元在地面任意点产生的磁场，是板状体所有磁荷面在该点产生磁场的总和。因此，计算板状体的磁场可先计算磁荷面的磁场，并将各个面求和，并且通过上面的分析，磁场的表达式还可由引力位与磁位的泊松公式导出，这里直接给出二度均匀直立板状体的磁异常解析表达式： (2.16) (2.17)式(2.16)、(2.17)中各个参数的物理意义下文将详细介绍。2.2 岩石的岩石磁性特征及场源磁模型通过对磁异常的解释及分析，用以解决地质找矿及构造分析问题，必须对岩石的磁性（磁源）特征有所了解。岩石和矿石磁性的强弱可用磁化强度表示，组成岩石和矿石的矿物质有顺磁性、反磁性及铁磁性。岩石在地壳中形成后，经历了较长的时期，使它磁化的地球磁场无论在强度或方向上都发生过变化。岩石还曾经受过高温高压及机械力的作用。因此，岩石中如果含有铁磁性的矿物，它就可能有一定的剩余磁化强度或称永久磁化强度，剩余磁化强度与岩石的经历及其铁磁性矿物含量有关，与现今地磁场大小及其方向无关。通常，在磁异常资料解释工作中，剩余磁化强度一般情况下都很弱，磁性源的磁化方向通常不会受其影响，所以不考虑剩余磁化的存在。此外，岩石还在现今的磁场中被磁化而具备一部分感应磁化强度，它与现今的正常地磁场成正比，方向也与地磁场的方向一致39。在不考虑剩余磁化的情况下，岩石的感应磁化强度可表示为： (2.18)式中，是现今地磁场总强度（磁感强度）；是岩矿石的磁化率，它取决于岩矿石的性质，决定了岩矿石被磁化的难易程度；为真空磁导率，在制中的大小为。地磁场总强度在地球表面上一般情况下是已知的，且相当大范围内可认为是不变的，是一已知常数。磁化率是一个无量纲的量，因此可知感应磁化强度的方向与地磁场的方向一致。对于地学断面的岩石磁性资料不多，刘蓓莉44和余钦范等45分别对阿尔泰-台湾地学断面和格尔木-额济纳旗地学断面的岩石磁性进行收集整理和测试，通过对比分析，两条断面上的岩石磁性参数均有共同特征（表2.1）： 表2.1岩石的磁性(刘蓓莉, 1997)岩类岩石名称磁化率/10-5(SI)沉积岩碎屑岩9.37碳酸盐岩5.20侵入岩花岗岩107闪长岩2030辉长岩2207火山岩流纹质凝灰岩01800流纹岩06600英安质凝灰岩08200变质岩片岩12.2千枚岩11.5（1） 沉积岩盖层和一些浅变质岩(如大理岩、片岩、砂岩、灰岩和硅质岩等)磁性都很低。（2） 各类中酸性岩浆岩，包括安山岩、闪长岩、花岗岩、花岗闪长岩、闪长玢岩和中酸性火山岩及凝灰岩等，具有中等磁性。（3） 基性和超基性岩浆岩，包括橄榄岩、辉长岩、玄武岩等，磁性较强，磁化率为。这些岩石的磁性变化大，不稳定。（4） 深变质岩也具有较强的磁性，磁化率为。当这类岩石大面积出露时，会引起大范围的磁异常。（5） 太古界(结晶基底)具有强磁性，地表岩石的磁化率为。此外，前人研究表明，结晶基底的磁性并不是向下一直延深至地核，而是在一定深度以后，由于地温随深度增加，结晶基底中的铁磁性矿物因达到居里点而变为顺磁性矿物，由于顺磁性矿物较弱，一般认为是无磁性物质。达到居里点时的界面成为居里面，居里面可以认为是磁性底界面或者太古界(或结晶基底)底界面。目前还有一种观点认为，莫霍面是地壳磁性层底界面46。根据上述岩石磁性特征，不管采用哪一种观点作为磁性底界面均可建立如图2.3所示的地壳磁模型。该模型主要有四部分组成:无磁性沉积层、局部磁性体、磁性层以及磁性层下界面以下的无磁性岩层。图 2.3 地壳场源模型无磁性沉积层实际上包括沉积岩和浅层变质岩，具有一定磁性，但相对岩浆岩或者铁磁性物质而言要小得多，因此认为是无磁性的。局部磁性体主要指铁磁性矿物、岩浆岩的侵入体(主要指由岩浆岩充填的断裂带)等。磁性层主要是由深成变质岩和结晶基底组成，其上界面可能为结晶基底的顶界面，下界面可能为居里面，也可能为莫霍面。磁性层上界面深，下界面深度变化较大(如表2.2所示)。由于磁性层的位于地壳深部，埋深较大，我们可以认为其磁性变化不大，他所产生的异常只是由于顶底界面的起伏引起，在定量解释过程中将磁性层假定为均匀磁性层。表2.2 根据大陆磁模拟计算不同地区居里面深度变化范围(管志宁等, 1997)北 美苏 联中国英 国20-48km15-55km20-50km15-26km通常，我们把磁性层因顶底界面起伏产生的异常称作为区域异常；深部岩浆岩一定通道(比如深大断裂)侵入到地表，或者规模较大的磁铁矿也会产生规模较小航磁异常，我们称之为局部异常。航磁的区域异常特征为：磁异常横向变化较为平缓，频率较低；局部异常特征为：横向变化较为剧烈，频率较高。通过以上分析，地壳岩石可能的磁化强度一般为，但是可能个别地区磁化强度会偏高，达到39。2.3 磁异常地表观测到的磁场是不同磁场叠加的综合，去除短周期变化的磁场后，磁异常就是观测地磁场与正常背景场之间的差异。磁异常是地下空间各种磁性异常体受到地球磁场的磁化作用后，在其附近区域产生并添加在正常背景场上的次生场。它属于固体地球内部的稳定场的一部分。磁测数据中，其内包含了两部分的磁异常数据：大范围的区域异常和小范围的局部异常。它们是由分布大小不同的磁性体产生的。磁力勘探中，可对不同的地磁异常分量进行观测，磁异常分量分为垂直异常（用或表示）分量、水平磁异常（用或表示）分量以及地磁场总强度模量差（用表示）。其中，垂直磁异常分量是实测地磁场垂直分量（）与正常磁场垂直分量（）的模量差。水平磁异常分量是实测地磁场垂直分量（）与正常磁场水平分量（）模量差。在航空磁测和海洋磁测以及目前的野外地面磁测中，基本都是测定总磁场磁感应强度的相对变化值,称为总磁场模量异常。在数值上等于实测总磁场感应强度与正常场总磁场感应强度之差。在实际工作中，之所以选择测定总磁场的模量异常，是因为它更容易精确地确定， 是与的模量差，即： (2.19)通常情况下，当总磁异常矢量值小于为正常背景场值的数十分之一时， 就可以近似把看成是方向的分量且的方向在区域不大的情况下一般认为是不变的（一万平方公里内约变化一度），所以，可把看作在固定方向的分量，即相当于在固定方向的分量。这样的物理意义与、类似，都是在固定方向的分量。在航空磁测中，磁异常通常在数百至一千纳特之间，以上关系的近似程度比较高。因此，对于二度体而言，的表达式为： (2.20)式中，、分别为在和上的分量，、分别为正常场的倾角和偏角。第3章 二维介质中的磁化强度成像3.1 引言磁化强度成像的方法类似于密度成像、地震层析成像方法，它是在剖面范围内，将地下介质剖分成大小相同、物性均匀的矩形截面的二度磁性体模型，通过最优化反演方法，在空间域对所建立的目标函数进行求解，最终得到地下二维空间的磁性源的分布情况。地表测量的磁异常是由地下磁性体在地磁场的作用下引起的，而地下介质是连续分布的，地下空间中的磁性体可视为由一系列较高或者较低磁化强度块体组成的磁异常区域（场源）。根据地面观测的磁异常反演地下场源体分布时，可将观测区域地下空间介质划分为一系列块体单元，当块体单元足够小时，假设每个块体单元内部磁化强度均匀分布，允许不同块体单元具有不同的磁化强度。根据位场理论，当不考虑剩磁和退磁影响时，磁化强度均匀的块体单元所产生的磁异常与磁化强度成正比，比例系数仅仅与块体单元形状及空间坐标有关。当块体单元大小位置固定时，比例系数可以预先确定出来，在反演的过程中，磁化强度值是唯一描述地质体和需要反演的参数。这样，根据磁异常反演磁性体的问题就转变为求解磁化强度在地下二维介质中的分布问题，这种反演方法称为磁化强度成像。前面我们已经推导出磁异常的解析公式，其解为精确解。地球物理正演是是一个由源求场的过程，即已知场源和地下物性分布，计算地球物理场的空间和时间分布。已知地质体的物性及其赋存状态可称为模型。 图3.1为磁异常正演的示意图。图3.1左边是地下二维磁性体的正演模型，通过解析公式得到地下磁性异常体在地面的磁异常，地面观测的磁异常特征如图3.1右边所示。其中为磁场的正演算子。图 3.1 模型正演示意图地球物理反演问题则是由场求源的过程，采用某种技术方法对实测数据进行处理，从而判断地下目标体的物性分布。解决反演问题的前提是先解决实测数据与地球物理模型的联系问题，确定它们之间的函数关系，所以正演问题的解决是进行反演的基础，无论何种形式的反演方法都必须要依赖于正演问题的解决。在没有解决正演问题时，实测数据与地球物理模型参数就无法联系起来，反演就会出现盲目性。只有解决了正演问题才能通过线性或非线性反演方法解决反演问题。地球物理反演属于一种定量的解释工具，通过有效的反演方法，能较可靠的重构地下未知场源的特征。图3.2为磁异常反演示意图。图3.2左边为观测数据，通常在观测数据中不但包含了有用的地下场源信息，并且还具有高频的噪声以及人为或仪器的一些干扰因素等。根据处理后的实测观测数据，通过反演算子重构地下未知场源的特征，图3.2右边的二维图为反演重构的磁性异常体。图 3.2 模型反演示意图在反演过程中，反演算子的形式可能会很复杂。一般有以下两点导致了反演算子的复杂性：1.物理模型过于复杂且一般是非线性问题；2.通常反演具有多解性即非唯一性，即不同模型都能够拟合同一组实测数据。所以，在实际操作中，为了减少非唯一性，使得结果较为可信，一般要采用正则化方法添加一些其他合理、有用信息进行约束反演，这将增加反演算子的复杂性。重磁反演固有的多解性主要由以下三个方面的因素影响。1.通常重磁反演属于不适定问题，所得的方程是欠定的，属于欠定问题，并且数据也是有限的,而地下模型则是连续的,对于离散反演而言，反演模型个数多于测量的数据点个数；2.地表采集的数据不可避免地存在高频噪声和其他非人为因素的干扰；3.地表采集的数据是地下所有场源的综合效应,位场反演存在固有的多解性。通过以上分析可知，要想获得相对可信的解，必须进行约束反演。3.2 2D磁化强度成像的基本原理及正演计算3.2.1 模型构建及正演计算 本文讨论的是二度磁性矩形体的正反演问题，所谓二度体即磁性异常体沿走向无限伸长，且沿走向，矩形柱体的埋深、截面形状、磁性变化特征稳定的情况。这样我们只需要计算剖面上的磁场就可以了。二维的模型构制示意图及其异常曲线如图3.3所示，将地下介质网格化，这样磁性体就是横截面为矩形的二度体，具体的网格划分大小与数据采集的间隔有关，下文将作详细介绍。在地面可测得水平分量、垂直分量以及总强度异常。在实际工作中，总强度异常一般比较容易精确地测得，本文模型反演所采用的数据均为通过正演所得的总强度异常。图 3.3 二维模型构制及其异常示意图矩形单元内的磁化强度均为常数，允许不同矩形单元存在磁性差异。图3.3中矩形单元的大小为，红色区域内的磁化强度为，其他区域为。地磁倾角，二度体走向与磁北夹角,即有效磁化倾角。异常曲线为、模型的正演曲线。从图3.3中可以看出，在反演的过程中，我们将地下空间网格剖分为致密排列的二度直立板状磁性单元后，根据地表所采集的数据，通过磁化强度成像方法，得到单元的磁化强度，从而确定地下场源的分布规律。单个直立板磁性单元如图3.4所示。图 3.4 单个板状体单元截面示意图 二度体的正演公式如式2.16、式2.17及式2.20，则第个直立板状体在观测点处的磁场为： (3.1) (3.2) (3.2) 为第个板状体的磁化强度；为地磁场倾角；为地磁场在地面上的投影与测线方向的夹角；为第个板状体的有效磁化强度倾角，； 、分别为第个板状体截面角点、到观测点的距离，、分别为、与轴正方向的夹角，由轴顺时针算起。在直立板状体的水平半宽度、垂直半宽度、中心位置及观测位置已知的情况下，这些参数可以表示为：那么整个地下磁性体在地面观测点的总磁异常为： (3.4) (3.5) (3.6)其中，直立板状体模型的个数。 通过以上分析，第个直立板状体在观测点处的总强度异常磁场可表示为： (3.7)为第个板状体的磁化强度，由（3.1）、（3.2）、（3.3）式可知，是由第个直立板状体的几何参数及其中心与观测点的距离所决定的核系数矩阵。我们知道，在观测点处所测得的磁场为地下个直立板状体在该点的总和，即 (3.8)我们假设观测点的个数,(3.8)式可写成矩阵形式： (3.9)即 (3.10)这里，为观测剖面上的维磁异常向量；表示维核系数矩阵；表示各个直立板状体的维磁化强度矩阵。3.2.2 模型网格剖分问题在磁化强度成像中，对地下模型网格剖分后，其大小及位置保持均不变，并且与观测点之间的位置关系也不变。网格剖分的大小应根据测点间距大小、最小目标体的大小以及计算速度为参考设置，而测点间距的设置与需要观测的磁性地质体的规模大小有关。对于二维模型而言，矩形是任意二维地质体的基本单元，模型的离散化主要有两种方式：第一种是矩形横向长短固定，整个模型空间由同层相等、各层不等的矩形水平分层排满，且随深度曾大各层模型的尺度也增加，从而使地下不同深度的矩形单元在地面上的磁异常位于相同数量级上，因而能克服随深度衰减的影响；第二种划分方式与第一种类似，只是将所有的网格单元的大小都划分为相同大小的矩形，采用合适的深度加权函数来消除随深度衰减的影响。第一种网格划分方式反演时的分辨率随深度加大而减小，因此第二种方式被广泛采用。 在网格剖分大小时，受到两方面的限制，一是分辨率的限制；二是计算速度的限制。剖分单元太大分辨率不足、不足以客观地反映地下磁化强度分布情况，而剖分单元太小则单元数量太多，计算量大，影响效率。Boulanger和Chouteau(2001)讨论了块体单元大小设置的两个准则：（1） 单元划分时应尽量小，使得观测数据中高频成分的信息能在反演结果中得到重构；（2） 单元划分时应足够大，这样节省反演时间，提高效率。事实上，观测数据的波长决定了单元大小的划分。浅层地质体影响观测数据中的短波长，而观测数据中的长波长则是深部地质信息的反应。Bhattacayya和Leu根据采用功率谱，分析了三维模型体顶层的能够对观测数据进行拟合的奈奎斯特频率的最大块体，指出单元划分大小应满足 (3.11)其中，为观测数据在方向的采样间距，为矩形单元在方向上的边长。这条准则同样适用于二维模型，则为矩形单元横向上的长度，为观测数据的点距。本文所剖分的正演模型大小都满足这个原则。3.3 磁化强度成像的反演技术3.3.1 共轭梯度法 共轭梯度法也叫共轭斜量法，最早是为了求解线性方程组而提出来的。当为对称正定矩阵时，求解上述方程组的解的问题与求解二次泛函的极小问题是等价的。共轭梯度法就是合理挑选一个方向向量，经过有限次迭代就能达到一个极小值。如果不是对称正定矩阵，则先将方程变换成。下面给出不加约束条件的共轭梯度法计算步骤：（1） 给定观测数据，初始模型；（2） 令,为迭代的最大次数；（3） 计算观测值与初始模型正演理论值之间的残差；（4） 令，；（5） ,共轭梯度（CG）迭代开始，为最大迭代次数； ； ； ； ； ； 如果，则迭代满足精度要求，跳出循环。这里为给定的误差限。 ； （6） 共轭迭代结束，得到，即 （7）令，更新初始模型，回到第三步重新迭代，直到满足给定的误差限。 理论上，共轭梯度法经过有限次的迭代应得到精确的结果。但实际上由于舍入的影响，结果并不如此。3.3.2 反演模型设计为了检验反演算法的有效性和可靠性，本小节设置了一个基础模型，假设地下截面为矩形的板状体走向为南北走向，其截面大小为，沿走向无无限延伸，观测剖面沿东西方向，长。数据的采集间距为，共采集41个数据点，即从剖面最左端的处开始，每隔采集一个点，一直到剖面右端的处。矩形柱体位于观测剖面的中间（图3.5红色区域），其顶端距地面。矩形柱体的磁化强度为，其他区域为，有效磁化倾角。理论模型及其正演总强度异常曲线如图3.5所示。通过对模型进行正演计算，然后利用正演的数据进行反演。在进行反演时，需要对地下进行网格划分，上面对网格划分大小的讨论可知，划分单元不能太大，也不能太小，太大则分辨率不够，太小则待求量过多，计算量过大，影响效率，所以按照上文所提到的网格划分规则进行划分。网格的大小为，则纵向划分为20份，横向为40份，共有个网格单元，如图3.5所示。图 3.5 基础模型及其正演异常曲线通过反演算法，得到模型各个单元的磁性分布，如果反演比较理想，则红色区域内的磁化强度应为，其他区域为。3.3.3 无约束的剖面反演计算 通过上面的分析我们已经得到网格化模型的正演计算公式3.10，在无约束的情况下，利用正演所得到的数据进行反演，采用共轭梯度法求解目标函数3.10式，所得的没有加入任何约束条件的反演结果如图3.6所示。 图 3.6 无约束反演结果及其异常曲线上图的磁化强度成像结果中，物性分布趋于地表。在该模型的试验中，地表测点数据有41个，而地下网格单元有800个，未知数的数量大于方程的个数，属于欠定问题，所以理论上求解目标方程时会有无数个解。正演的异常曲线与反演重构后的异常曲线非常吻合，肉眼无法区分，均方差达到了。可见，反演时仅仅拟合观测曲线是不够的，有时虽然观测曲线拟合得很好，但是重构的结果却与实际不符。反演结果富集地表，与实际情况相差甚远，这并不说明计算是错误的，而是说明了该结果是欠定问题的众多解的一种解释而已。所以，要想使成像的位置与模型异常体的位置较好地对应，获得可靠的反演结果，约束条件的添加必不可少。3.3.4 深度加权函数在空间上，地面或空中观测的磁异常数据是地下不同深度所有磁性地质体所产生的叠加异常，一般情况下，浅部小的异常体与深部大的异常体产生的叠加场难以区分，磁场数据在深度方向上的分辨率低。例如，对于最小模型目标函数，反演结果将使“能量”倾向于集中在地表附近，这显然与实际异常体位置不符。这是因为构造模型的核函数是线性的，磁异常值与场源到观测点的距离呈指数衰减，导致核矩阵中数值随深度增加而急剧减小，相同的网格单元，深部的异常响应应比浅部的要弱得多，对观测数据的贡献相对很小，故容易出现这种现象。由于这种现象的影响，虽然反演得到磁化强度模型可以很好地拟合观测数据，但缺乏深度方向上的分辨率。针对这种情况，我们采用深度加权函数来补偿这种衰减，使反演重构的结果有一定的纵向分辨率。其基本思想是将核函数中由深度引起而衰减的部分抵消。为了说明深度加权的有效性及必要性，本文采用一维的数学模型对其作了验证。假设有一组数据,其正演公式为公式(3.12)。 ， (3.12)式中，核函数为 (3.13)在核函数中，由于衰减因子的存在，使得核函数分量随着的增大而不断减少。因此，重构模型时，非零值大多集中在取小值的区域，如图3.7直接反演结果所示。可知，模型参数是核函数的线性组合，即 (3.14)公式(3.14)中，随着的变大而不断衰减，我们采用加权的目的就是去掉这部分的衰减。为了消除衰减系数的影响，改造公式（3.12），则有 ， (3.15)则加权后的模型参数为 (3.16)结合公式(3.14)，可得到经加