重磁三维反演技术(2010118)

重磁反演方法是一种地球物理勘探方法，用于研究地下的重力和磁力场。

它通过测量地球表面上的重力和磁力数据，推断地下的密度和磁性分布。

重磁反演方法的基本原理是根据地球物理学的基本方程，建立地下密度和磁性分布与地表重力和磁力场之间的关系。

然后，通过数学模型和计算方法，将地表观测数据转化为地下模型的参数。

在重磁反演方法中，常用的数学模型包括正演模型和反演模型。

正演模型是根据地下密度和磁性分布计算地表重力和磁力场的模型，而反演模型则是根据地表观测数据反推地下密度和磁性分布的模型。

重磁反演方法的应用范围广泛，可以用于研究地球内部的结构、地下矿产资源的勘探、地下水资源的调查等。

它在地质勘探、矿产勘探、环境地质等领域具有重要的应用价值。

需要注意的是，重磁反演方法是一种间接方法，其结果受到多种因素的影响，如观测误差、模型假设等。

因此，在实际应用中需要结合其他地球物理勘探方法和地质资料进行综合分析，以提高解释的准确性和可靠性。

重磁解释中健全的极大似然法及三维反演
方盛明;余钦范;楼海
【期刊名称】《物探与化探》
【年(卷),期】1994(18)2
【摘要】本文介绍了ＪＯＡＯ．Ｂ．Ｃ．ＳＩＬＶＡ提出的重磁解释中健全的（Ｒｏｂｕｓｔ）极大似然反演方法的基本原理，讨论了重力二度体非线性几何参数的反演情况，并将该方法扩展到三维重、磁非线性几何参数反演中。

重点讨论了场值中有较大随机干扰、有不同类型的地质噪音源影响下的反演结果。

理论反演结果表明，该方法有比最小二乘反演方法更接近实际的解，从而提高了解的可靠性，对于重磁局部异常或孤立异常有较好的反演效果。

【总页数】1页(P142)
【作者】方盛明;余钦范;楼海
【作者单位】不详;不详
【正文语种】中文
【中图分类】P631.24
【相关文献】
1.重磁联合反演及重磁与MT综合解释巴彦浩特盆地火成岩 [J], 邓荣来;李庆浩;宋桂桥;刘天佑
2.重磁联合反演解释在长白山天池深部构造中的应用 [J], 孙帮民;吴燕冈;管彦武
3.重磁电震联合建模正反演在火成岩解释中的应用 [J], 陶德强;赵文举;张嵘鑫;柴艳宁
4.大地电磁三维反演在土屋铜矿床探测中的应用研究--三维反演和解释 [J], 谭捍东;李聪嫔;王彩旗
5.重磁震联合反演及综合解释在海域构造解释中的应用 [J], 骆迪;张训华;蔡峰;刘展;张志珣;侯方辉;杜润林;魏龙
因版权原因，仅展示原文概要，查看原文内容请购买。

重磁三维反演伊犁盆地中部密度和磁性结构索奎;张贵宾;梅岩辉;马勇胜【摘要】伊犁盆地中部区域地处多个构造单元的交汇区,构造复杂,缺乏高精度地质资料,制约了对盆地构造演化的进一步认识.本文联合使用多种先验信息作为约束条件,利用高精度重磁数据反演获得了研究区地下海拔-10 km以浅的三维密度和磁性结构,增强了反演结果的可靠性.结果揭示以白石墩次凸为中心发育了一个“北断南超”的南西-北东向不对称型凹陷,北部沉积厚度大于南部;区内主要断裂具有高角度特征,生烃中心受断裂控制,阿吾拉勒山前凹陷浅部可能发育逆掩断裂;火成岩以侵入的基性火成岩为主,伴有部分喷出火成岩和中酸性火成岩.【期刊名称】《地球物理学报》【年(卷),期】2018(061)008【总页数】10页(P3410-3419)【关键词】伊犁盆地;重磁;三维反演;断裂;火成岩【作者】索奎;张贵宾;梅岩辉;马勇胜【作者单位】华北水利水电大学资源与环境学院,郑州 450046;中国地质大学(北京)地球物理与信息技术学院,北京 100083;中国地质大学(北京)地球物理与信息技术学院,北京 100083;中国地质调查局油气资源调查中心,北京 100029;中国地质调查局油气资源调查中心,北京 100029【正文语种】中文【中图分类】P312;P6310 引言伊犁盆地是在天山造山带所夹持的微地块上发展演化而成的山间叠合盆地，属于伊犁—中天山微板块，矿产资源丰富，有较好的油气前景.多位学者对伊犁盆地及邻区的构造特征及形成演化等问题进行了研究，取得了丰富的成果.廖世南(1992)认为伊犁盆地是具有陆内裂谷和山间陆相的双层结构盆地，初步讨论了盆地的生成和发展；刘俊霞等(1997)利用航磁资料对盆地中部的基底结构、起伏变化、构造特征及盖层分布等问题进行了研究；张国伟等(1999)认为伊犁盆地分为中上元古界变质基底、中下石炭统裂谷火山岩系褶皱变形基底和盆地沉积岩系的三层结构，并将伊犁盆地主要构造划分为二山三盆；左国朝等(2008)基于野外调查结果和前人研究结论，研究了包括伊犁盆地在内的西天山地区的地层系统及演化过程.伊犁盆地阿吾拉勒地区出露了多处的火成岩，为研究盆地构造演化提供了丰富的证据(Yang et al.，2014；Ge et al.，2015；李永军等，2017).熊小林等(2001)研究了阿吾拉勒中酸性火成岩的地球化学特征，认为可能是晚古生代后碰撞阶段幔源岩浆底侵作用导致玄武质岩石熔融形成的；李注苍等(2006)基于大哈拉军组火山岩地球化学特征推断火山岩位于岛弧区，为板块俯冲造山带构造环境，这是目前的主流观点之一(李锦轶等，2006；Wang et al.，2007；Sun et al.，2008)；李晓英(2013)基于花岗岩的分析结果，认为阿吾拉勒地区在晚古生代不同时期存在岛弧、碰撞伸展和大陆裂谷环境；韩琼等(2015)认为阿吾拉勒大哈拉军山组时代跨度大，火山岩自西向东由老变新，铁矿床受火山活动控制.此外多人对伊犁盆地的油气前景进行了研究评价，钟红利(2011)基于二叠系以来古构造演化、烃源岩、储集层的沉积特征和油气的运移聚散特点，在伊宁凹陷划定了三个有利油气区；王亚军(2013)将伊犁盆地及邻区的构造演化划分为七个阶段，认为北部坳陷是演化条件最好的区域，上二叠统和中上三叠统是相对较好的油气层位.前人的成果为伊犁盆地的深入研究提供了良好的基础，但多数尺度较大，针对局部重点区域的精细研究有待加强.本文研究区域位于伊犁盆地中部多个次级构造单元交汇区域，以白石墩次凸为中心，包含了乌孙山隆起、阿吾拉勒凸起、伊宁凹陷的南部斜坡带和巩留次凹等构造单元的部分区域(图1).研究区是雅马渡—白石墩北西向断裂和赛里木湖—莫合尔北东向断裂这两条主要断裂系的交汇处，断裂带造成了白石墩凸起，分隔了伊宁凹陷和巩乃斯凹陷，并控制和影响了其他构造单元的构造分布形态和中新沉积(张国伟等，1999).由于研究区大部分被第四系地层覆盖，因此需要利用地球物理方法对本区域地下构造进行反演，获取精细构造，以期为盆地的构造演化等相关问题的深入研究提供地球物理证据.1 数据与方法1.1 数据本次重力测量设计网格为500 m×500 m，磁测量设计网格为500 m×250 m，在白石墩区域加密为250 m×250 m，实际重力测点数5004个，磁测点数10025个.对原始数据分别进行各项改正后得到了布格重力异常和ΔT磁异常，如图2所示.图1 伊犁盆地中部区域地质简图1 第四系砂土、砾石； 2 水西沟群含煤碎屑岩沉积； 3 石炭纪灰岩、碳酸盐岩、砂砾岩； 4 二叠纪花岗岩； 5 小泉沟群砂岩夹泥岩、砾岩； 6 二叠纪闪长玢岩； 7 二叠纪辉长岩； 8 晓山萨伊组、哈米斯特组、塔木其萨伊组、巴斯尔干组并层； 9 研究区范围； 10 断裂.Fig.1 Geological sketch map of central Yili Basin1 Quaternary System sandy soil，gravel; 2 Shuixigou Group; 3 Carboniferous Period limestone, carbonate rock, sandstone and conglomerate; 4 Permian Period; 5 Xiaoquangou Group; 6 Diorite porphyry; 7 Diorite gabbro; 8 Xiaoshansayi Formation, HamisiteFormation, Tamuqisayi Formation, Basiergan Formation; 9 Study area; 10 Fault.图2a显示研究区布格重力异常为负值，根据异常特征被分为五个区域.其中有三个高异常区域，分别是位于阿吾拉勒凸起的a、b和位于乌孙山隆起的c；有两个低异常区，分别是位于伊宁凹陷的d和巩乃斯凹陷的e；以及一个过渡交汇区f.与图1中阿吾拉勒凸起的范围相比，a范围向西北方向有所扩大，且a和b之间区域的低异常非常明显；c区与图1中乌孙山隆起相比向东北方向扩展较多，这可能是乌孙山隆起在地下的延伸的反映，也可能是伊宁凹陷和巩乃斯凹陷之间的次凸的反映；d和e为两个显著的低值区，与f区的交界处暗示了凹陷的边界位置.从图2b中可以看出，研究区磁异常总体较为平稳，其中有三个较为显著的高磁异常区g、h和k.g区与高重力异常区b吻合较好，磁异常略显杂乱，可能与阿吾拉勒凸起中的高磁性岩石分布有关；h是磁异常最高值区域，异常曲线较为平缓，闭合完整，可能是地下有一定埋深和体积的高磁性岩体反映；k区磁异常幅值相对较低，且较为散乱，反映了磁源可能埋深浅且体积较小.本次研究区岩石物性测量共采集标本321块，密度和磁化率统计结果如表1所示.结合东方地球物理公司在伊犁盆地进行的地质物性调查结果，总结研究区的主要物性特征如下：地层从上到下由新到老，密度逐渐增大；存在主要三个密度界面，第四系和新近系、三叠系和二叠系以及二叠系和石炭系；酸性岩浆岩为中低密度，中基性岩浆岩为高密度；新生界、中生界和古生界正常沉积地层为无磁性-微弱磁性地层，古生界二叠系火山岩地层为强磁性地层，岩浆岩由酸性到基性磁性逐渐增强.1.2 数据处理为了研究区域构造格架，本文使用三维反演方法获取地下密度和磁化率分布情况.基于实测数据分布范围和分辨率，确定反演垂向范围为地表到海拔-10 km，水平范围则为测区.由于研究区地表海拔高度从692 m到1593 m，差异较大，本文采用了带地形的三维物性反演方法，目的是尽量避免数据处理过程中引入的误差影响反演结果.图2 伊犁盆地中部重磁异常等值线图(a) 布格重力异常等值线图；(b) ΔT磁异常等值线图.Fig.2 Contour map of gravity and magnetic anomalies in central Yili Basin(a) Bouguer gravity anomaly contour map；(b) ΔT magnetic anomaly contour map.表1 伊犁盆地中部地层岩石物性参数统计表Table 1 Summary of densities and magnetic susceptibility measured from samples in central Yili Basin地层岩矿石名称标本数(件)密度(g·cm-3)磁化率(10-5SI)Q砂砾岩482.21217N砂岩、泥岩、砂砾岩102.37841J砂岩、泥岩、砾岩242.44333T砂岩、泥岩、砾岩402.39429P2tP2xP1砂岩、泥岩、页岩、砾岩、玄武岩、安山岩、安山玢岩252.54358712.59797182.621935C2C1d安山岩、安山玢岩、玄武岩、火山角砾岩、流纹岩、砂岩、泥岩、砾岩、灰岩222.649358522.678480岩浆岩花岗岩112.606123注：密度和磁化率均为几何平均结果.首先需要提取反映海拔-10 km到地表的密度不均匀体和磁性异常体产生的重磁异常场作为三维反演的依据.反演范围的顶面(地表)有起伏，侧面和底面则为直立或水平，为了避免插值运算引入的误差，尽可能保留原始数据的有效信息，经与多种先验信息对比后使用反演法(Li and Oldenburg，1998a；Li et al.，2012)获得了浅源的重力场和磁场.该方法的优点是提取到的浅源场数据均处于原始测点位置，避免了多次插值导致的误差，确保了高程数据准确，提高了三维反演数据的精度. 地形数据来自实测高程数据和STRM DEM数据融合.为了减少测量误差及浅层噪声对反演结果的影响，同时避免进行带地形反演时出现“观测点在地下”的情形，在异常分离时分别对浅源重磁场向上延拓了150 m，各个测点的高程也增加了相同高度.1.3 三维反演1.3.1 基本原理本文使用基于空间导数约束的重磁物性三维反演方法(Li and Oldenburg，1996，1998b)进行计算，该方法能够加入多种约束条件并实现大数据量的快速反演.反演目标函数如式(1)，通过最小化目标函数可以获得最佳解模型.(1)式中φd为数据拟合差函数；φm是模型目标函数；μ是正则化因子，用于平衡拟合差函数和模型目标函数的权重；最后一项为界限函数，λ为界限参数；aj是j单元的物性下限，bj为j单元的物性上限，mj为j单元的解模型，ln函数限制mj必须在aj和bj之间，形成了“硬约束”.数据拟合差函数φd为φd=‖Wd(G m-dobs)‖2，(2)式中Wd为各个单元的权重diag(1/σi)，σi为第i个单元的标准差；G是核函数矩阵；m是待求的物性模型；dobs为实际观测数据.式(1)中模型目标函数φm为(3)式中m0是参考物性模型；wr是距离加权函数；αs是模型长度系数；αx，αy和αz是模型光滑系数.从式(1)—(3)可以看出，该方法可以加入多种先验信息来对结果形成有效约束，包括根据先验信息建立的参考模型(m0)、参考模型的可靠程度(Wd)、解模型的数据范围(aj和bj)、解模型与参考模型的差异程度(αs)、解模型在各个方向的光滑程度(αx，αy和αz)和噪声水平调节(μ)等，从而获得最大程度逼近实际情况的解模型. 1.3.2 反演参数(1) 反演网格反演时地下剖分网格的设计需要基于实测数据分辨率和解释的需求，兼顾计算效率，同时考虑目标体的复杂程度和地形的起伏情况等因素.其中数据分辨率最为重要，实测数据的最大分辨率为数据的间隔，无法反映小于最小间隔的异常体，而且一般在浅部确定一个异常体的存在，该异常应该覆盖2～5个网格单元的范围.试验结果表明，反演网格尺寸小于数据采集间隔的1.2倍时都可以达到数据采集的最大分辨率(Boulanger and Chouteau，2001).基于上述因素，本文对重磁三维反演设计了相同的剖分网格：中心区域水平方向分辨率为300m ×300 m，四周各向外扩展3600 m；垂向分辨率为150 m(浅部)和200 m(中深部).反演网格范围是x：-3600～41400 m，y：-3600～39000 m，z：-11400～1600 m，共计1080576个网格单元.(2) 约束条件本文使用1 ∶20万的地质图基于物性统计结果建立了参考模型，该模型以地表信息为主，没有先验信息的深部区域设置参考模型的物性为0，且有先验信息的区域权重高于未知区域.基于物性统计结果设置解模型密度变化范围为-0.5～0.5 g·cm-3，磁化率变化范围为0～0.1 SI.考虑到研究区域以层状地层为主，设定重力反演模型的长宽高比值为10 ∶10 ∶1.磁反演的目标主要是区域内的火成岩体，经过初步反演后确定了区域内有两个主要的火成岩体，根据岩体的比例设定模型的长宽高比值为8 ∶10 ∶8. (3) 噪声水平模型试验结果表明，每个测点分别给定适当的噪声时，反演结果精度略优于全区域给定相同水平噪声的结果.本文根据实测数据的精度对噪声总体水平进行估计，再利用基于能量阈值的二维小波方法估计了每个测点的噪声，一般在异常变化平缓的区域噪声水平相对较低.最后根据反演结果对正则化因子进行微调，直至获得与先验信息拟合最佳的反演结果.2 反演结果及评价2.1 重力反演结果使用CPU为I7 2600，内存为8G DDR3，操作系统为Windows 7 64位专业版台式计算机进行计算，重力三维反演耗时36分33秒，获得的密度扰动结果如图3所示，显示出以下几个主要特征：密度扰动变化幅度较小，过渡平缓，但高低密度区域界限明显，为划分断裂位置提供了依据；高密度体主要分布在阿吾拉勒突起和乌孙山隆起区域，低密度区域则主要与凹陷带吻合；不同深度密度特征继承性较好，随着海拔高度的逐渐降低，横向密度不均匀性逐渐降低.2.2 磁反演结果使用与重力反演相同配置的计算机，磁三维反演耗时3小时6分54秒，获得的磁化率分布如图4所示，从图中可以看出：高磁化率区主要集中在东北部和东部两个区域，向下延伸较深且为整体；中部区域异常分散，幅值相对较低，埋深较浅；磁化率垂向变化较明显，差异随深度的增加逐步降低；伊宁凹陷区域在中深部有低幅值异常体存在而浅部没有.2.3 结果评价重力反演数据拟合差在-0.825～0.537 mGal之间，磁反演拟合差在-17.33～30.08 nT之间；在数据幅值较大的区域的归一化误差水平高，幅值小的区域归一化误差水平较低，与理论值相符合；两个解模型较为平滑，且与其他先验信息基本相符；这表明本次三维反演参数选择较为合适(Williams，2008).图3 伊犁盆地中部不同海拔高度密度扰动分布图(a) 海拔180 m； (b) 海拔-500m； (c) 海拔-1100 m； (d) 海拔-2100 m； (e) 海拔-5100 m； (f) 海拔-10100 m.F1、F2、F3、F4和F5为本文推断构造单元边界；ga、gb和gc为高密度异常区.Fig.3 Density distribution of different altitude in central Yili Basin(a) Altitude 180 m； (b) Altitude -500 m； (c) Altitude -1100 m； (d) Altitude -2100 m； (e) Altitude -5100 m； (f) Altitude -10100 m.F1, F2, F3, F4 and F5 are the inferred boundaries of different tectonic units. ga, gb and gc are high density abnormal areas.图4 伊犁盆地中部不同海拔高度磁化率异常分布图(a) 海拔180 m； (b) 海拔-500 m； (c) 海拔-1100 m； (d) 海拔-2100 m； (e) 海拔-5100 m； (f) 海拔-10100 m.ma、mb和mc为高磁化率异常区.Fig.4 Magnetic susceptibility distribution of different altitude in central Yili Basin(a) Altitude 180 m； (b) Altitude -500 m； (c) Altitude -1100 m； (d) Altitude -2100 m； (e) Altitude -5100 m； (f) Altitude -10100 m.ma, mb and mc are high magnetic susceptibility abnormal areas.3 讨论3.1 盆山构造和断裂从图3中可以看出研究区的主要构造格架，伊宁凹陷和巩乃斯凹陷呈现低密度特征且沉积层厚度较大，阿吾拉勒凸起和乌孙山隆起为高密度特征；随着深度的增加，密度的横向不均匀性逐渐降低，整体性增强；推断的主要断裂和构造单元边界基本位于高低密度异常过渡区.研究区北部边缘属于阿吾拉勒凸起，地表出露的是二叠系晓山萨依组(P2w)，由砂岩、砾岩和火成岩组成，密度范围变化较大.本次反演揭示了研究区内的高密度岩体，主要分布在ga、gb和gc三个区域.ga区域高密度体体积较小，埋深较浅，横截面积从海拔180 m到-5100 m几乎没有变化，再向下至-10100 m逐渐消失，这表明阿吾拉勒凸起在该处边界接触面倾向可能是北东向，F1为高角度逆断层.ga 东北部未进行数据采集，没有反演结果，根据反演的趋势结合地质图推断高密度体可能向北东向继续延伸.研究区东部阿吾拉勒凸起地表主要有二叠系乌郎组(P1w)的中基性火成岩出露，均为中高密度，反演结果显示gb区域存在一个整体性较好、体积较大的高密度体，在浅部(180～-1100 m)呈“C”字型，在深部则为一个整体，横截面积基本不变，结合磁反演结果，推测该区域浅部可能存在逆冲推覆构造，方向为从东北向西南，深度达到海拔-1100 m.该岩体一直向下延伸到-10100 m仍有显示，与白石墩次凸的接触面倾角较高，这表明F2为高角度正断层，倾向为南西向.该区域在本次测区的边缘位置，推测其在东北部和东南部应会继续延伸.gc 区域位于研究区西南部，范围最大，地表基本被新生界覆盖，其下方的高密度体体积大、埋深深，在-3100 m以浅高密度体为分散分布，随着深度的增加逐步融为一体，横截面积逐渐变大，向下延伸至-10100 m仍有明显的显示.该区域是乌孙山隆起与白石墩次凸的过渡区域，深部的高密度体可能是乌孙山隆起的延伸，据此本文认为乌孙山隆起的界线与前人划分的相比向北西向最远延伸了5 km.刘俊霞等(1997)认为该区域可能为逆掩推覆构造，乌孙山北缘断裂向北逆冲了15 km，本次反演结果显示为超覆构造，可能由于浅部的磁异常体存在干扰了推断结果.研究结果显示研究区内断裂均以中浅层断层为主，从反演结果可以看出区内几条主要断裂向下延伸的深度基本不超过11 km，属于盆地内部次级断裂，且各条断层角度较高，这与前人的结论基本一致.在前人确定构造单元边界位置的基础上，本文认为F1和F2断裂为阿吾拉勒凸起的边界；F3和F4为伊宁凹陷南部斜坡带和乌孙山隆起的边界；F5则可能是巩留次凹的边界.与前人划分的边界位置相比有变化，这是因为考虑了构造单元在地下延展情况.综上，本文认为研究区以白石墩次凸为主发育了一个南西—北东向不对称型凹陷，特点是北断南超，越靠近北部断层处沉积厚度越大，与阿吾拉勒凸起之间为高角度逆断层接触关系，且北部的山前凹陷浅部可能发育逆掩断裂；凹陷南部属于超覆沉积，厚度逐渐变小直至尖灭，超覆于乌孙山隆起之上；区域内沉积层厚度较大，生烃中心的主要层位是中生界和上古生界(王亚军，2013)，受到区内次级断裂的控制，其埋深从北向南逐渐变浅.3.2 火成岩研究区域内出露的中基性火成岩主要是位于阿吾拉勒凸起的二叠系火成岩，由安山玢岩、玄武岩、玄武玢岩、火山碎屑岩以及凝灰砂岩等组成；中酸性火成岩主要位于研究区域南部，主要是由安山玢岩、安山质凝灰岩及碎屑等中酸性喷发火成岩组成的下石炭统大哈拉军山组(C1d)，局部有酸性或基性的变化(王晓地，2001).物性资料统计结果显示，中基性火成岩呈现高磁化率、高密度的特征；中酸性火成岩的磁性变化范围较大，但一般是弱磁性显示且平均密度低于围岩的密度.本文以磁性反演结果为基础，结合重力反演结果圈定了6处中基性和1处中酸性火成岩分布区域，如图5所示.图5 伊犁盆地中部火成岩分布推断结果示意图1 已知基性火成岩； 2 已知中酸性火成岩； 3 推断基性火成岩； 4 推断中酸性火成岩.Fig.5 Deduction of igneous rocks distribution sketch map in central Yili Basin1 Known basic igneous rocks； 2 Deduction of intermediate-acidic igneous rocks； 3 Known basic igneous rocks； 4 Deduction of intermediate-acidic igneous rocks.从图5中可以看出，出露的中基性火成岩主要分布在阿吾拉勒凸起的范围内.图3和图4显示，图5中J1区域下方存在一个高密度、高磁化率岩体，主要形态特点是地表横截面范围较大，随着深度增加到海拔-500 m附近时横截面变小，再向下到-5100 m附近时横截面最大，深度再增加时横截面变小，在海拔-10100 m时有零星显示，推断为中基性火成岩，可能是底侵玄武质岩石(熊小林等，2001).J1的位置与已知火成岩出露区域基本相同，范围略大.J2区域地表全部为第四系覆盖，反演结果显示J2下方是与J1区域相似的高磁化率异常体，横截面积从地表到海拔0 m左右变化不大，再往深处逐渐变大，在海拔-2100 m横截面积最大，并在-5100 m附近分为南北两个岩体，南部岩体在-10100 m仍有显示.测区边界处为低重力异常，推测可能是由于低重力异常区位于巩留次凹的范围，盆地引起的重力低淹没了可能存在的高密度体引起的重力高，因此推测J2区域大部分为隐伏的中基性火成岩，且可能延伸到测区外.J3区域地表也为第四系覆盖，该区域以高重力异常和小范围不连续高磁异常为主，推断中基性火成岩分布较为零星，没有形成较大的岩体，且各个岩体埋深较浅，主要在海拔-1100 m以浅，推测可能为喷出岩，对深部构造影响有限.J4区域位于研究区西南部，有相对较弱的平滑磁异常场.在地表无明显高磁异常体存在，单纯根据磁异常数据处理结果难以确定J4下方深部是否存在火成岩，反演结果显示存在高磁性、高密度体，MT结果显示为高电阻率岩体，推断存在基性火成岩体，可能为侵入岩.J5范围较小，其所在区域在地质图上显示为中二叠统的晓山萨依组(P2x)，以碎屑岩—碳酸盐岩为主，但J5区域有较为明显的高磁异常和高重力异常显示，据此推断该区域有可能存在小规模的中基性火成岩，埋深适中.J6区域位于研究区域北部，属于阿吾拉勒凸起，根据地质图看该区域出露的为石炭系中酸性火成岩，而重磁异常(图2)和反演结果显示该区域为高密度高磁性岩体，研究表明本区域中的石炭系中酸性火成岩中也有部分中基性火成岩(王晓地，2001)，因此将其划分为中基性火成岩，范围较地质图上标明的范围有所扩大，向下延伸到海拔-1700 m左右逐渐消失，在其北部可能会继续延伸.研究区下部S1区域中部也有一个和J6类似的中基性火成岩体，磁三维反演结果显示为高磁异常，重力反演结果显示为高密度异常，MT结果显示为高电阻率异常，据此推断为中基性火成岩体.图1显示S1区域出露了石炭系大哈拉军山组，文献(白建科等，2015)显示乌孙山地区大哈拉军组有酸性-基性-中性-酸性的火山喷发旋回，本文划分的中酸性火成岩分布区域完全位于该范围之内，属于乌孙山隆起边缘位置，主要为弱磁异常和低密度异常显示；推断S1中部有中基性火成岩的主要原因是其显示了高磁性和相对高密度，这也与文献提到的局部基性变化相符(王晓地，2001).综上，本文认为区域内中深部的火成岩以侵入岩为主，这与前人的结论基本一致(李晓英，2013；熊小林等，2001；李凤鸣等，2015)，J3区域为喷出岩的可能性较大.4 结论本文利用高精度重磁数据在多种先验信息的约束下反演，揭示了伊犁盆地中部地下三维密度和磁性结构，为该区域构造演化的进一步研究提供了深部证据.基于反演结果本文认为：研究区为南西—北东向不对称型北断南超凹陷，沉积层厚度较大且北部沉积层厚度大于南部，南部沉积超覆于乌孙山隆起之上；区内主要断裂为高角度次级断裂，北部山前凹陷可能发育逆掩断裂；区内火成岩以侵入的基性火成岩为主，伴有部分喷出火成岩和中酸性火成岩.本次反演尽可能利用了已知的各种约束信息，限于数据精度和算法，所得解模型与实际情况会存在差异，随着勘探程度的不断加深，更多的先验信息的加入，会使得反演结果越来越接近实际情况.致谢感谢中国地质调查局油气资源调查中心提供了部分观测数据和相关资料，并且在研究过程中提供了诸多帮助.感谢两位审稿专家提出的宝贵意见. ReferencesBai J K, Li Z P, Xu X Y, et al. 2015. Detrital zircon U-Pb dating of Dahalajunshan Formation in Wusun Mountain region, western Tianshan, and its geological implications. Geology in China (in Chinese), 42(1): 85-95. Boulanger O, Chouteau M. 2001. Constraints in 3D gravity inversion. Geophysical Prospecting, 49(2): 265-280.Ge S S, Zhai M G, Safonova I, et al. 2015. Whole-rock geochemistry and Sr-。

大地电磁测深数据和重力数据三维联合反演汇报人：日期:•引言•大地电磁测深和重力数据采集与处理目录•三维联合反演的理论和方法•实验和结果分析•结论和展望•参考文献01引言大地电磁测深和重力数据在地球科学领域的应用大地电磁测深和重力数据是地球科学领域重要的数据来源，对于研究地球内部结构、地壳厚度、地幔流动、地核状态等具有重要意义。

三维联合反演的必要性传统的二维反演方法在处理复杂地球内部结构和多参数反演时存在一定的局限性，因此需要采用三维联合反演方法以提高反演精度和可靠性。

国内外研究现状目前，国内外学者在大地电磁测深和重力数据联合反演方面已经取得了一定的研究成果，但仍存在一些问题和挑战，如多参数反演的复杂性、数据分辨率和信噪比等问题。

研究背景和意义•国内外研究现状：目前，国内外学者在大地电磁测深和重力数据联合反演方面已经取得了一定的研究成果，如基于波动方程和射线理论的反演方法、全波形反演方法等。

同时，随着计算机技术和数值计算方法的发展，越来越多的学者开始关注三维联合反演方法的研究和应用。

发展趋势：未来，大地电磁测深和重力数据联合反演将朝着以下几个方向发展 1. 高分辨率和高质量的数据采集技术；2. 更加精确和可靠的反演算法和技术；3. 多参数、多尺度和多角度的综合反演方法；4. 人工智能和机器学习等新技术的应用。

1 2 3研究内容：本研究旨在利用大地电磁测深和重力数据，开展三维联合反演方法的研究和应用。

具体研究内容包括1. 大地电磁测深和重力数据的预处理和分析；2. 三维联合反演算法的建立和优化；0102034. 实例应用和效果评估。

创新点：本研究具有以下创新点3. 反演结果的分析和解释；011. 提出了一种基于波动方程和射线理论的联合反演方法，提高了反演精度和可靠性；022. 开发了一套完整的三维联合反演软件系统，实现了自动化和批量化处理；033. 对多种地球内部结构和参数进行了反演实验和分析，验证了方法的可行性和有效性；044. 将研究成果应用于实际地球探测任务中，取得了良好的应用效果。

重磁遗传算法三维反演中高速计算及有效存储方法技术
姚长利;郝天珧;管志宁;张聿文
【期刊名称】《地球物理学报》
【年(卷),期】2003(046)002
【摘要】将地下场源区域规则划分成很多小长方体单元,并且通过反演确定这些单元的物性变化,勾画出场源的分布图像,这种方式逐步成为重磁反演,特别是三维反演的重要方向;遗传算法等非线性技术进行该类反演将逐步成为发展趋势. 本文指出,在应用遗传算法进行该类反演过程中,隐含着数据量较大时超常规的计算量,它已成为制约该类反演充分发挥作用的瓶颈问题;同时,本文提出了针对性的分离并存储几何格架的计算策略、以及独特的几何格架等效压缩存储技术,可以从根本上提高非线性反演计算速度,为该类反演的有效应用奠定了坚实的基础.
【总页数】7页(P252-258)
【作者】姚长利;郝天珧;管志宁;张聿文
【作者单位】中国地质大学地球物理学院,北京,100083;中国科学院地质与地球物理研究所,北京,100101;中国地质大学地球物理学院,北京,100083;中国地质大学地球物理学院,北京,100083
【正文语种】中文
【中图分类】P631
【相关文献】
1.重磁异常三维物性反演随机子域法方法技术 [J], 姚长利;郑元满;张聿文
2.重磁遗传算法三维反演中动态数组优化方法 [J], 姚长利;郑元满
3.重磁反演约束条件及三维物性反演技术策略 [J], 姚长利;郝天珧;管志宁
4.重磁异常正反演可视化实时方法技术改进 [J], 姚长利
5.重磁电三维反演在下嘎来奥伊河上游矿区中的应用 [J], 李少朋; 李桐林; 郑军; 陈汉波; 张忠禹
因版权原因，仅展示原文概要，查看原文内容请购买。

重磁三维可视化反演系统的设计与实现
陈建国;肖敦辉;梁玉辉
【期刊名称】《地质学刊》
【年(卷),期】2012(36)3
【摘要】重磁三维可视化反演解释系统一直是国内外重磁勘探领域的研究重点之一.介绍了基于Visual C++与OpenGL环境研发的重磁三维可视化反演系统,详细阐述了系统的设计思想与实现的关键技术,包括可视化技术、面向对象技术、图形拾取技术、碰撞检测技术及反演约束技术等.系统具有人机交互几何反演与最优化物性自动反演两种功能,可满足目标异常、区域模拟和盆地建模.
【总页数】6页(P250-255)
【作者】陈建国;肖敦辉;梁玉辉
【作者单位】中国地质大学地质过程与矿产资源国家重点实验室,湖北武汉430074;中国地质大学(武汉)资源学院,湖北武汉430074;中国地质大学地质过程与矿产资源国家重点实验室,湖北武汉430074;中国地质大学(武汉)资源学院,湖北武汉430074;中国地质大学地质过程与矿产资源国家重点实验室,湖北武汉430074;中国地质大学(武汉)资源学院,湖北武汉430074
【正文语种】中文
【中图分类】P631.2;TP391
【相关文献】
1.合肥盆地重磁资料处理及重磁震联合反演 [J], 李云平;刘金连;林治模;胡加山;阎红;夏吉庄
2.重磁联合反演及重磁与MT综合解释巴彦浩特盆地火成岩 [J], 邓荣来;李庆浩;宋桂桥;刘天佑
3.重磁反演约束条件及三维物性反演技术策略 [J], 姚长利;郝天珧;管志宁
4.重磁三维自动反演软件系统的分析与设计 [J], 郑元满;姚长利;张晨
5.任意形状地质体数值积分法重磁场三维可视化反演 [J], 杨宇山;刘天佑;李媛媛因版权原因，仅展示原文概要，查看原文内容请购买。

火山盆地三维重、磁数据反演一、三维重、磁反演的可行性分析大量的研究工作表明火山盆地的变质基底和火山岩盖层间之间存在明显的密度和磁化率差异：变质岩基底为高密度、低磁性，而火山岩则为低密度、相对高磁性，因此可以利用重力和磁法数据圈定火山盆地的变质基底、火山通道等区域构造格架。

二、拟开展的主要工作1、岩石标本采集与测试根据研究区的具体情况，争取本着钻孔岩芯标本为主，地表岩石标本为辅的原则，每种地质单元（岩性至少采集50块标示进行密度、磁化率和电阻率测试。

2、重磁资料的收集整理为开展区域重、磁数据三维反演，必须收集研究区1：5万地面高精度重力和磁力测量原始数据或经数据预处理后的网格化数据和数字化的研究区地形资料，以及深钻孔编录资料和研究区主要岩石的物性资料。

3、三维重、磁数据反演（1）、重、磁场特征分析（2）、重磁三维物性反演三维物性反演拟采用UBC-GIF反演软件。

UBC-GIF软件是由英国哥伦比亚大学地球物理反演工作室研制的带稀疏先验信息约束的3D 重磁物性反演软件，分为MAG3D和GRAV3D两个模块。

该软件灵活性比较好，可以将不同的地质勘探阶段所获得的地质信息加入反演中，进行约束反演，从而使反演结果能够拟合观测数据，而且与已知的地质信息相吻合，更好地反映地下真实的地质信息。

具体反演流程：（3）、三维地质-地球物理模型构建三维地质-地球物理模型的构建的思路是利用三维重磁物性反演获得的物性属性数据体，获取各个方向的切片图，以切片图为基础，综合各种地质地球物理资料，构建合理的二维地质地球物理初始模型，然后将所有的2D模型组合成3D地质地球物理模型，以此模型作为约束属性模型进行再反演计算，直止反演精度满足要求，且模型与地质认识相一致。

具体建模型流程为：三、经费建议1、岩石物性标本采集与测试300块*500元/块=15万元。

2、重、磁三维反演中国地质调查局没有预算标准，我们承担的中国地质调查局相山铀矿田三维地质填图试点项目，盆地面积为582m2，预算经费为80万元。

姚长利,郑元满,张聿文.重磁异常三维物性反演随机子域法方法技术.地球物理学报,2007,50(5):1576～1583Y ao C L ,Zheng Y M ,Zhang Y W.32D gravity and magnetic inversion for physical properties using stochastic subspaces.Chinese J .G eophys .(in Chinese ),2007,50(5):1576～1583重磁异常三维物性反演随机子域法方法技术姚长利1,2,郑元满1,2,张聿文1,21地质过程与矿产资源国家重点实验室和地下信息探测技术与仪器教育部重点实验室,北京　1000832中国地质大学(北京)地球物理与信息技术学院,北京　100083摘　要　本研究针对三维物性反演中存在的困难和问题,提出三维物性反演的随机子域方法技术,首先是将正反演中保持不变的几何格架分离计算并存储,避免重复计算,从而提高正反演计算速度;其次是利用对称性等实现等效计算,明显降低格架计算和存储要求;再通过随机子域方式,降低反演的维数问题;另外,通过概率方式控制子域生成的分布,实现约束新机制.模型和实例计算表明了方法技术的效果,为大面积重磁数据的三维反演提供了有效的途径.关键词　重磁数据,三维反演,随机子域,约束条件文章编号　0001-5733(2007)05-1576-08中图分类号　P631收稿日期2006-11-02,2007-04-30收修定稿基金项目　国家自然科学基金项目(40374039)和新世纪优秀人才支持计划(NCET 20420726)资助.作者简介　姚长利,男,1965年生,教授,主要从事重磁勘探研究与教学工作.E 2mail :clyao @32D gravity and m agnetic inversion for physical propertiesusing stochastic subspacesY AO Chang 2Li 1,2,ZHE NG Y uan 2Man 1,2,ZHANG Y u 2Wen 1,21State K ey Laboratory o f G eological Processes and Mineral Resources ,G eo 2detection Laboratory ,Ministry o f Education o f China ,Beijing 100083,China2School o f G eophysics and In formation Technology ,China Univer sity o f G eo sciences ,Beijing 100083,ChinaAbstract F ocused on the existing difficulties and problems in 32D inversion for physical properties ,the paper brings forward stochastic subspace methodology for 32D inversion for physical properties.Firstly ,it com putes separately and saves the geometric trellis which keeps unchanged during forward simulation and inversion to av oid repetitive com putation s o as to increase the speed of forward simulation and inversion com putation.Secondly ,it uses symmetry to realize equivalent com putation ,which distinctly lowers the requirement of trellis com putation and storage.And thirdly with stochastic subspace inversion method it reduces the number of dimensions of inversion.In addition ,it controls the distribution of the subspaces generated through probability method to realize the new mechanism of constraint.The com putations of m odel and field data dem onstrate the effect of the methodology which is hopeful to be of practicality 2oriented 32D inversion for physical properties of large scale gravity and magnetic data and meets the requirement of explanation of 32D inversion.K eyw ords G ravity and magnetic data ,32D inversion ,Stochastic subspace ,C onstraints第50卷第5期2007年9月地　球　物　理　学　报CHI NESE JOURNA L OF GE OPHY SICSV ol.50,N o.5Sep.,20071　引　言 随着对重磁反演及定量推断解释要求的逐步提高,基于二维的重磁反演已经难以满足全方位确定地质构造精细展布图像的高要求,所以重磁反演已逐步发展到三维反演阶段.重磁异常反演的模型主要为形态模型和物性模型两类,在三维反演中,形态模型以多面体模拟地质体为主,通过反演实现多面体形态的变化,以达到逼近地质体的目的,多面体的物性通常是给定不变的.而物性模型则以将场源区域划分成小的单元组合(主要是长方体或立方体单元),在反演过程中,单元的形态不变,物性发生变化,通过物性变化勾画场源范围.比较而言,由于三维物性反演具有模型物性易于操作、能模拟任意复杂地质体的能力,反演方法技术受限制条件少,不用涉及到复杂的形态变化,因而这种反演方式已成为重磁反演尤其是三维反演近年的一个主要方向[1～11].但是,重磁反演通常是一项困难的工作,并且在三维物性反演中表现得更加突出.首先就是反演的多解性更加严重,因为三维反演对象是面积性测量数据,数据量远远大于二维情况的剖面数据.另一方面,三维反演的目的也是为了更精细地刻画地质构造研究对象,其结果是剖分的组合模型数量往往需要很大,造成计算中解空间的维数非常大.这样,重磁反演的多解性这个本质问题在这类反演中暴露得更加明显.已有的研究表明,增加约束条件可能是减少多解性最好的方法.从数学角度看,约束的作用是给数据提供额外信息,从而有效求解多数情况下是病态的反演问题.如果从反演效果上看,约束的作用则是引导反演的进行过程,在多解的“十字路口”选择合理的方向.除此之外,目前物性反演所采取的方法主要还是局部最优化方法,面对三维反演,其局限性已严重制约了物性反演的效果.因为当反演面积数据量大,模型单元多时,超大规模线性方程组求解变得非常困难.举例说明:如一个测区数据为100×100的网格(这只是一般大小的面积数据量),作物性反演时取简单的模型对应剖分关系:100×100×50,即深度上有50层剖分模型,水平每一层为100×100的组合单元体,模型单元水平位置与测点位置对应,总的物性单元数则为500000个.采取如广义逆求解最小二乘法反演方法,将反演目标函数转换成线性方程组,其系数矩阵的数组大小将为500000×500000,该超大规模线性方程组的数学求解实际上几乎无法进行(一些针对性的压缩技术无法从根本上解决这个问题).除此算法问题之外,简单换算可知,该系数矩阵对计算机的内存要求是极其巨大的(25×1010,即单精度数据会达到1000G B!),在很长的时间内这几乎也很难实现.所以,现有的三维物性反演研究中,异常数据和剖分的模型单元数都相当少,这势必严重影响三维反演的实际应用.局部最优化方法的另一个问题是与约束结合时存在的困难.带约束的反演目标函数形式如F=Fd(m)+αF(1)m(m)+βF(2)m(m)+…,其中m为模型物性参数向量,Fd(m)表示数据的吻合程度,其单位为nT (磁)或mG al(重力),F(1)m(m)、F(2)m(m)等则为模型参数的各种约束评价函数,它们可能会包括模型的空间关系信息或(和)物性信息,显然其物理量纲是复杂的.由于Fd、F(1)m和F(2)m等物理量纲的不同,如何调和难以把握,所以正则系数α,β等的数值往往是经验的.除此之外,复杂的导数运算9F9m i= 9F d(m)9m i+α9F(1)m(m)9m i+β9F(2)m(m)9m i+…又增添了不少麻烦,一方面其表达式繁琐、计算量很大,另一方面,为了进行导数计算,一些地质约束条件的提取和数学表达变得很困难,所以造成实际中一些约束无法清晰描述,并难以与反演有效融合.比较而言,非线性算法反演则具有明显的优势,如遗传算法和控制随机算法这类非线性方法,是按照一定的随机原则产生解模型空间的采样点并进行评价的,理论上具有全局寻优的特点;另外,由于是通过正演计算实现反演的,非线性全局优化方法不需要对高维目标函数进行复杂的求导计算,从而大大降低了与各种约束条件结合的技术“门槛”,必然会简化地质约束条件的数学描述,更加容易地结合一些过去难以结合的地质、地球物理约束条件.但到目前为止,非线性技术在地球物理反演中实际应用效果并不很理想,在重磁三维物性反演中的研究成果还很少,究其根本原因是其普遍存在的计算瓶颈问题.通过分析研究可以看出,非线性反演方法是要通过巨大的正演计算量来避免导数计算并且实现对解空间的访问搜索的.但在三维物性反演中,解空间的维数是如此之高(例如前面例子中所面对的反演问题,解空间就是n500000维的,其中n为一维参数的取值数),对解空间状况要达到相当的了7751　5期姚长利等:重磁异常三维物性反演随机子域法方法技术解,其访问搜索量之巨是无法承受的.所以,要想从根本上取得突破,只有大幅度降低反演解空间的维数.基于以上分析,我们根据物性反演的特点,针对性地提出新的重磁三维物性反演思路,称之为“随机子域法物性反演”方法.这里的子域是指整个反演模型区域中的某个小区域,它通常包括一定数量的剖分模型单元.2　方法原理 根据非线性反演的物性反演方案,在计算中不再有导数运算,只涉及到规则长方体的正演计算.在教科书上就可以找到其计算公式.以三维密度模型的正演公式说明剖分模型的正演情况,如图1所示的组合三维模型区域,其中任一地质体单元j 在坐标观测点p (x ,y ,z )的重力异常为Δg j (x ,y ,z )=σj S j (x ,y ,z ),(1)其中,S j (x ,y ,z )=G∑2l =1∑2m =1∑2n =1(-1)l +m +n×(x l -x )ln[(y m -y )+R lmn ]+(y m -y )ln[(x l -x )+R lmn ]+(z n -z )arctan(z n -z )R lmn(x l -x )(y m -y),G 为万有引力常数,σj 为第j 个单元模型的密度,R lmn =(x l -x )2+(y m -y )2+(z n -z )2.图1　反演模型示意图(a )组合模型;(b )任一模型单元.Fig.1　The m odels(a )The m odel consisting of a set of juxtaposed cells ;(b )One of the cells.211　降低计算和存储要求的等效计算方法在三维物性反演中,每一种变化的模型单元物性分布,都需要计算其重磁异常效应,长方体单元的重磁场正演计算量是相当大的.由于非线性反演中时刻进行着大量的目标函数采样正演计算与评价,加快正演计算的措施,必然提高反演计算速度.仔细分析重磁物性反演的特点,分析上述正演计算表达式,其中S j (x ,y ,z )由场源的几何参数与计算点坐标组成的特点,我们称之为几何格架.不难发现,在物性反演中,一旦模型的剖分关系确定下来后,其几何形态及与测点的相对关系将始终保持不变.针对这个特点,采取存储模型几何格架的策略解决重复计算问题,即将所有模型的几何格架计算一次并存储起来,然后,后续的、反复的正演计算就变成了几何格架与对应的物性非常简单的一个乘积运算,模型的正演计算量几乎消失,这样,简单的存储策略就达到了极大地提高计算速度的目的,从而极大地提高了非线性算法物性反演的能力.上述几何格架分离技术,合理利用了物性反演的特点,从根本上提高了后续反演计算时的速度,从而为非线性方法在三维物性反演这个特定场合下的应用奠定了一个关键条件.但是,在实际三维物性反演中,几何格架即使只计算一遍,其计算量往往也是非常大的.同时,该策略在三维物性反演中却引发了新的困难———巨大的几何格架存储问题.以前面三维模型剖分为例,如测区网格数据为100×100规模,模型剖分仍为100×100×50的三维网格状,在这里提出的格架分离存储方案中,每一个单元体有100×100个几何格架计算及存储点,所有模型总的几何构架存储量将是5×109.如果以单精度(4Bytes )保存,数据量达20G B ,双精度数据(8Bytes )则高达40G B ,仅几何格架就需要这么大的存储量,最近若干年的计算机内存还无法承受.另外,我们在测试该组合模型几何格架计算时(应用的是P 2III (1G )、256M B 内存的计算机),仅几何格架计算一遍就需要近16h !如果要减少计算量和存储量,只能极大地减少模型的剖分数,这样必然会使反演模型的刻画能力明显降低.因此,在三维反演情况下,分离的几何格架直接存储方法是不可取的.但仔细分析,如果面积数据是水平分布的(否则需先作重磁异常曲化平转换),并且使模型的剖分与数据网格采取某种对应关系(如取一一对应),则会存在同一层模型各单元之间几何构架具有特定的等价性.利用这个等价性进行等效存储,使每一层的构架存储量减少到只相当于一个模型单元的存储量.下面对其进行具体解释分析.8751地球物理学报(Chinese J.G eophys.)50卷　取模型为均匀划分,且与规则网格的测点对应,m 、n 为测区网格行列大小,设水平观测面上任一测点p (k ,l )(k =1,2,…,m ;l =1,2,…,n ),与地下模型单元σ(k ,l )对应(见图2所示,这里只画出一层),任意单元σ(i ,j )(i =1,2,…,m ;j =1,2,…,n )在网格点p (k ,l )的几何格架为S i ,j (k ,l )(i =1,2,…,m ;j =1,2,…,n ;k =1,2,…,m ;l =1,2,…,n ).显而易见,该几何构架是一个四维数组(再考虑到深度上分层,则变为五维数组),其数据量往往极其巨大.图2　模型单元与测点网格一一对应关系Fig.2　Relationship between observation stations andcells on a layer 分析上述均匀剖分模型重磁场的计算公式,可以发现以下平移等效性(示意关系见图3所示):S 1,1(k 0,l 0)=…=S i ,j (i +k 0-1,j +l 0-1)=…=S m -k 0+1,n -l 0+1(m ,n ),(2)其中,i =1,2,…,m ;j =1,2,…,n ;k 0=1,2,…,m ;l 0=1,2,…,n .利用(2)式,可以避免重复计算和存储,从而极大地减少存储量.另外,对于重力异常,以及垂直磁化的磁异常垂直分量,可以进一步利用简单的对称性(见图4所示),因而有以下等价关系式:S i ,j (i +k 0-1,j +l 0-1)=S i +k 0-1,j +l 0-1(i ,j ).(3)这样,综合(2)、(3)两式,任何几何格架的计算都可以通过以下等价计算公式得到S i ,j (k ,l )=S 1,1(|k -i |+1,|l -j |+1),(4)其中i =1,2,…,m ;j =1,2,…,n ;k =1,2,…,m ;l =1,2,…,n .也就是说,任意单元σ(i ,j )的几何格架S i ,j (k ,l )与特定单元σ(1,1)产生的几何格架S 1,1(k ′,l ′)(k ′=|k -i |+1=1,2,…,m ;l ′=|l -j |+1=1,2,…,n )完全相同.因此,事先计算好单元σ(1,1)产生的几何格架S 1,1(k ,l )(k =1,2,…,m ;l =1,2,…,n ),像表格一样保存起来,要计算任意单元σ(i ,j )的S i ,j (k ,l ),只要查表,按确定的关系式(4)找到正确的位置即可.经过这样巧妙的技术措施,将四维数组S i ,j (k ,l )变成等价的二维数组S 1,1(k ,l ),从而节省了大量的计算量和存储量.图3　平移等效性Fig.3　The equivalence indisplacement图4　互换对称性Fig.4　The equivalence in reciprocation(4)式证明,利用上述等效性,可以使每一层的几何构架存储量减少到只相当于一个模型单元的几何格架,因而大大减少了对内存的需求.对于上述100×100×50剖分的组合模型例子,双精度数据格式的格架存储量将由40G B 减少到4M B ,即原来的1Π10000,可见其减少量是惊人的.同时,由于可以利用等效单元的几何格架,没有必要再计算所有单元的几何格架,这样,格架分离技术中格架的初始计算也得到极大压缩(对于本例,为原来的1Π10000,计算时间将由16h 减少到约6s ).此项关键技术保证了物性反演中格架分离措施大幅度提高计算速度的顺利执行.212　降低反演维数的随机子域选择反演方式尽管有了上述大幅度提高计算速度的方法技术,但是在三维物性反演中,由于划分的物性小单元9751　5期姚长利等:重磁异常三维物性反演随机子域法方法技术数量往往很多,这样,反演的维数很高,解模型空间极大,前面已指出,必须降维,否则无法实现有效的反演,为此我们采取子域方案.例如将反演区域分成若干个小块———子域,分别对各子域进行独立反演,依此完成对全区的反演计算.但这样的子域反演是有缺点的,例如其确定的子域排列顺序反演就会造成场源分布主观上的倾向性,这是不公平的.为此,我们采取动态的随机子域法,即随机产生某个位置、某个尺度大小的子域,反演其物性变化,完成后再产生新的子域并继续进行反演,这样迭代进行.由于子域的位置和大小都是随机产生的,希望使整个反演区域都有相等的反演机会.针对每一个生成的子域,由于其维数已大大降低,另外,迭代方式使每个物性单元都会多次被包含,单个子域的具体反演方式不再是决定性的,所以可以采取非线性方法或线性方法进行反演.实际上为使计算更加简单,针对随机生成的子域,我们可以将其当成一个整体的均匀单一模型来反演,再将反演结果分配到所包含的组合单元中.这样,在计算上随机迭代子域反演将变得非常简单.213　格架权约束异常分离技术在应用之前,还有一个关键问题需要解决———子域反演所对应的重磁异常选择问题.在子域的具体反演时,如果不采取措施,先被选中的子域将会反演掉目标异常的绝大部分,后续选择的新子域将只对所剩很少的剩余异常进行反演,我们称其为子域的“贪婪现象”,这是必须防止的.如何从总的异常中选取与子域对应的部分?这是蕴含在反演中的事先无法回答的问题.但显然有一定的规律在里面,子域所对应的异常份额主要由三个因素组成:与子域的尺度大小有关,与子域的深度有关,另外,还与子域的物性大小分布有关.因为物性大小是未知的待求量,无法加以考虑,而子域的格架则包含了子域大小和子域深度信息,所以我们选择子域的格架在反演区域总格架中所占的比例作为加权系数,从总异常中选择相应份额的异常作为子域的反演目标,进行加权反演,克服子域的贪婪现象,弱化随机子域生成次序上先后的影响,保证子域公平竞争的机会,实现子域反演方法.具体加权系数为W s (x ,y ,z )=∑n2j =n1S j(x ,y ,z )∑n 0i =1Si(x ,y ,z ),(5)其中,W s 为任意计算点的权系数,n 0为所有单元数,n 1,n 2为当前子域所包含单元的编号.由于模型产生的主要异常在其正上方附近,实际计算中为使问题简化,上述权重可以只取一个比值,即正上方格架数值之比代替随点变化的精确权重,这样会使问题简单,其误差会在后续的迭代中不断修正,迭代过程也就是不断修正的过程.根据子域格架提取相应的权重份额异常,可以理解为是一种新的约束,但不是对场源的约束,这一点与传统的反演约束机制不同.从约束的角度来评价,它不会带来额外的副作用,就是前面指出的常规约束存在物理量纲混乱以及相关的技术困难等.214　概率约束子域选取机制在随机子域物性反演中,如何结合约束?这是要加以研究的.约束的作用简单说就是希望产生什么样的结果,为此,我们采取概率子域生成方式体现约束.这里我们将约束倾向转化成子域的选取机会,以概率的方式体现出来.例如,Li (1996,1998)以及Pilkington (1997)在他们研究的方法技术中,为了避免反演结果中重磁场源集中在地表附近,在反演中结合了物性深度加权约束措施.此处概率约束方式就可以通过使深部具有更大的生成权,“倾向性”地让深部模型单元有更多的机会被选中参加反演计算.再比如,地质上或其他地球物理方法推断出某种倾向的构造分布(可称其为构造倾向性约束),这是属于外部附加的约束条件,在子域反演中就可以方便地将其转换成对应范围的相对概率分布,随机子域的生成相应地受此概率分布控制,选中子域后的具体反演过程与先前相同.可以看出,此处约束的作用体现在以机会的形式影响模型的生成上,至于所希望的模型能否得到支持并“站得住脚”,则要依靠目标函数评价,此处的约束只是起“引见”的作用.传统约束方式同时也可以保留,两者并不矛盾,只是作用的方式更加简单.例如,有些强制性约束可以通过对反演结果的滤波实现,如物性变化范围约束,以及最小构造约束等.3　模型计算 为了检验方法的计算效果,需要进行实际模型检验.由于物性反演方面的研究工作很多,也是近年来重磁反演的重点,为便于对比,我们设计了Li (1996,1998),P ortniaguine (2002)等多位研究者在论文中都使用的模型,受篇幅所限,下面选取一个0851地球物理学报(Chinese J.G eophys.)50卷　例子.图5为组合岩脉例子,图5a是模型产生的重力异常(其中附加了随机干扰),图5b为密度模型的一个垂直剖面图,两个岩脉互为反方向倾斜,图中红颜色延伸大的密度为110g・cm-3,黄色小的模型密度为018g・cm-3.图5c为该组合模型的一个水平切片图(深度为75m),可以看出大的模型水平走向也较长一些.测区平面数据网格大小为41×41,将反演区域剖分成41×41×20=33620个网格立方体小单元,在深度方向上为20层(由位场的特点,深部场源在地表产生的异常宽缓,浅部场源产生的异常窄陡.通常,测区范围内的异常主要由一定深度以上的场源产生,故一般取反演深度为水平范围的一半以内即可,更深的场源产生的异常可作为区域场去除).在反演中,为了使整个模型空间都能被随机子域覆盖到,需要有一定量的随机子域抽取次数,也就是迭代的次数.本例反演时的主要参数设置为:随机提取子域数量5×(41×41×20)=5×33620个,相当于模型总的单元数的5倍,子域的尺度在给定的范围内任意随机选取,该例中,子域中的单元大小为从1×1×1即只包含1个单元,到12×12×12即包含1728个单元;子域在水平位置为均匀概率分布,而深度分布上则取为线性增加概率分布,取从顶层的相对概率011到最底层的110.另外,附加强约束密度范围为010～210g・cm-3.迭代反演均方误差为01035mG al,反演约耗时20min(使用Pentium2III笔记本计算机,主频1G H z, 256M内存).另外,对采取快速算法前后计算速度进行了比较,针对41×41×20=33620个模型单元,计算一遍41×41=1681个网格计算点所需时间,采取等效几何格架快速计算技术后要65s,而没有采用时则要62min,速度提高很大.图6a～6c是反演结果的图示,其中图6a相当于y=300位置的东西向剖面,反映出的场源分布与图5c中y=300剖面位置吻合;图6b和6c分别相当于y=400位置和y=500位置的东西向剖面,将场源分布中两个场源的空间主要特点显露出来,图6c中左边场源的影响较大,接近与右边场源相连.比较表明,取得的结果与理论模型吻合得相当好.4　实际资料计算 图7是中国内蒙古布敦花地区某处航磁异常的反演情况图,图7a为实测磁异常图.反演中取该地区地磁倾角约60°,偏角约0°,反演时不考虑剩磁的影响,故磁化方向取地磁场方向.测区范围为12175 km×14175km,测区平面数据网格大小为60×53= 3180测点.将反演区域剖分成53×60×26=82680个网格立方体小单元,在深度方向上为26层.在反演中,参数基本上和前面模型计算时的一致,为提高计算速度,随机子域的尺度限制在不大于8×8×8个组合单元.另外,根据该地区的地质资料,磁化强度限制在(0～200)×10-2A・m-1范围.图7b是反演结果模型产生的磁异常,图7c是反演结果三维显示图(其中红色为磁性体中大于100×10-2A・m-1的磁化强度等值面分布状态,小于该值设为蓝色).图7d～7f为几个不同位置的垂直剖面图,为磁化强度分布情况(单位为10-2A・m-1).综合上述图件可以看出,地下磁性场源的分布及变化,从整体把握场源的空间展布关系.该岩体呈弯月形分布与地质情况及其他方法的推断相当吻合,其空间赋存状态细节有助于更进一步的详细地质研究.5　结　论 从以上模型计算和实际资料反演计算可以看出,我们提出的随机子域加权反演方法,具有方法简单、计算稳定、实用能力强的特点.其中快速计算和有效存储技术可以作为各种三维物性正反演方法的基础.另外,这里提出了两种新的约束方式:(1)根据格架权进行子域异常分离,是对重磁异常进行约束性选择,这是内部固有的约束;(2)子域选取生成的概率分布加权方式,这是外加的约束,如深度加权约束将深度补偿转化成模型生成概率补偿.以上两项约束,都没有影响目标函数的结构,较传统的约束简单,避免了反演目标函数中附加约束项造成的数据量纲混乱问题.1851　5期姚长利等:重磁异常三维物性反演随机子域法方法技术图5　重力模型例子Δg 为重力异常,h 为深度,ρ为密度.(a )加干扰的模型重力异常;(b )组合岩脉模型垂直剖面;(c )原始模型水平切片(深度75m ).Fig.5　Synthetic exam pleΔg is the gravity anomaly ,h is the depth ,ρis the density.(a )N oise 2corrupted gravity anomalies ;(b )The vertical profile view of the m odels ;(c )The horiz ontal section of the m odels at depth 75m.图6　重力异常反演结果的不同剖面显示(a )y =300m ;(b )y =400m ;(c )y =500m.Fig.6　The results of inversion ,which shown in the form of vertical profiles at different positions ,am ong whichthe y coordinate is equal to 300m (a ),400m (b )and 500m (c ),respectively图7　实测数据反演ΔT 为磁异常,M 为磁化强度.(a )原始磁异常;(b )反演结果模型磁异常;(c )反演结果三维立体透视图,及其三个剖面图(d )y =510km ,(e )y =710km ,(f )y =910km.Fig.7　The inversion of field dataΔT is the m agnetic an om aly ,M is the m agnetization.(a )T he observed m agnetic an om aly ;(b )T he calculated an om aly caused by the resulted m odels of inversion ;(c )A 32D view of the m odels ,and three vertical profiles at different y positions ,wh ose coordinates are 510km (d ),710km (e )and 910km (f ),res pectively.2851地球物理学报(Chinese J.G eophys.)50卷　。