矿集区深部结构、构造关系以及成矿环境是控制矿床形成的直接因素，因此，无论对于成矿学还是深部资源勘查，厘清深部三维形态和物性(包括断裂延伸、岩体的空间展布以及与地层的相互关系)都十分重要.一直以来，地质工作者试图利用各种已知信息建立地下三维结构模型创造出多种方法，其中，以重力观测数据为基础以其他先验地质信息为约束的三维(3 Dimensional，简称3D)重力反演，已经成为了解地球深部结构的一种重要手段(Oldenburg et al.，1997，2007； Fullagar et al.，2000，2007； Williams et al.，2004； Welford and Hall，2007).

按照反演单元划分，三维重力反演通常有网格节点(Voxels)和离散多面体(Discrete)两种方式(Li and Oldenburg，1998；姚长利，2002，2007).离散多面体反演因其既易于吸收先验地质、构造(包括地层倾向、断层和矿化体等)信息，又可使反演的理论场与观测场误差较小而受到欢迎.对于地质信息约束下的三维重力反演，不少学者做过工作(Roy and Clowes，2000；Jessell，2001；Goleby et al.，2002；Malehmir et al.，2006，2009；Martin et al.，2007；Williams，2008；Wang et al.，2011；吕庆田等，2014).综观过程，虽细节有所差别，但基本上可归纳为以下三步：①初始模型的建立，②反演的迭代与初始模型的升级，③模型的三维显示，其中步骤①和②至关重要.对于步骤①初始模型的建立，目前方法繁杂不一，尚未形成统一完美的解决方案，事实上，合理的初始模型能够减小后续重力反演的非唯一性，加快反演的收敛速度(McGaughey，2007；McInerney et al.，2007).步骤②则是对步骤①的进一步优化，以获得符合观测数据的物性参数空间分布.

贝叶斯方法(Bayes，1763)最初主要用于统计学领域，1987年法国Tarantola教授在其基础上创立了贝叶斯反演理论(Tarantola，1987)，之后，贝叶斯反演逐步应用于地球物理领域并得到广泛深化和应用(Green，1995；Scales and Snieder，1997；Gao，1997；Gouveia et al.，1998；Gr and is et al.，1999；Lasses et al.，2004；Tarantola，2005；Khodja et al.，2010；Stuart，2010；Fernández-Martínez et al.，2013).贝叶斯反演理论框架建立在概率论的基础上，它将数据信息与模型先验信息联系起来，通过模型的先验信息来约束模型的后验参数.由于贝叶斯推理把反演问题与观测数据信息、模型信息以及先验信息联系起来，使得求解约束泛函的极小点不需要计算梯度方向，因此，对约束项没有光滑性要求；同时，贝叶斯方法不只给出单一的反演值，而是给出解的后验分布采样，在此基础上得到解的条件期望、方差、置信区间等具有统计意义的结果，这使得贝叶斯方法能够很好地应用先验信息，得到丰富可靠的输出结果.将贝叶斯算法运用到三维重力反演中，应该有助于初始模型的改进.但是，在贝叶斯理论中，测量数据和先验信息包含在后验概率密度函数(PPD)中，它可以解释成模型的单点估计和不确定度等贝叶斯推断，这些信息的获取需要对反演问题进行优化求最优模型和在高维模型空间中对PPD进行采样积分；而且，从模型的后验概率密度函数中提取模型信息，需要计算它的模型估计、参数不确定度以及参数间的相关度等属性.所有这些，无疑会使贝叶斯方法受碍于维数，并且将为解决正演问题而付出高昂的计算成本.

为此，本文将开展以贝叶斯方法为基础的离散多面体三维重力反演研究，特别是在初始模型的建立中，吸收贝叶斯反演理论精华，将反演问题同观测数据、模型以及先验信息联系起来，计算中采用模型降阶技术克服反演的高维不确定性和复杂的正演成本，快速准确计算出地质体的密度和形态重现重力模型结构，再通过单位球体和任意形态几何体反演模拟实验以及安徽省泥河矿区应用实践，验证方法的正确性和实用性.

如前所述，贝叶斯反演是综合利用未知模型参数的先验信息、实测数据以及由物理规律得到的正演模型等信息来反演模型参数，最后得到未知模型参数的后验概率分布函数及其相应的特征量.设m表示未知的M维模型参数空间，d表示N维数据空 间，在贝叶斯反演中都看成是随机变量.实际应用中，数据d是已知的，模型参数化θ也是确定的，为了方便省略θ，得到贝叶斯公式(Carlin and Louis，2000)

这时P(d)是常量，P(d ｜ m)为已知数据 d 条件下随模型 m 变化的函数，称为似然函数，用L(m)表示.上式可以改写成

其中P(m ｜ d)表示后验概率分布函数(PPD).其原理可以简单概括如图 1.

图 1

Fig. 1

图 1 贝叶斯方法原理 Fig. 1 Bayesian method principle

似然函数L(m)(或P(d ｜ m))表达的是模型参数 m 与观测数据 d 之间的相对概率，在贝叶斯反演计算中，它的大小反映了模型响应与观测数据的适配程度，可见，似然函数 L(m)是贝叶斯反演计算的关键参量，模型 m 的似然函数可以表示为(Scales and Snieder，1997)

其中， g 为正演算子(又称核函数)， d 为观测数据矢量， m 为模型参数矢量，T表示转置， C _d为先验估计误差(以下公式中各参数皆如此).

引入目标函数E(m)，公式(3)似然函数L(m)可以改写为

其中

A是一个比例因子.

经典的贝叶斯反演算法就是寻求物理模型使得L(m)决定的概率尽可能大，等效于E(m)尽可能小(杨文采，1997；Yuji Mitsuhata，2004).采用最小二乘法求解，由，得到最大似然解

对应的协方差矩阵

可见，贝叶斯反演结果最终是由这个不确定度矩阵 C _m和最大似然解 m _ML共同定义，当数据误差服从Gaussian分布且是线性的，反演的结果也应满足Gaussian分布.实际反演通常存在欠参数化和超参数化两种情况：对于欠参数化情况，一般可采用贝叶斯信息准则(BIC)(Akaike，1980)来判断与数据信息一致的反演层数；对于超参数化情况，则采用倾向于简单结构的带先验信息反演，类似于正则化反演方法(Tikhonov and Arsenin，1977).本文带先验信息的重力反演属于后者.

通常，重力观测数据或其他先验地质信息是有噪音的，反演的模型结果也会有误差，假定这些噪音和误差满足Gaussian分布，在重力反演中似然函数L(m)、模型的先验概率P(m)、目标函数E(m)可以类推表达如下：

这里 m 代表待反演参数，d为重力观测数据， m _p为模型的先验估计值.当坐标原点位于地面，Z轴方向铅垂向下，X、Y轴位于水平面内，剩余质量为σ的地质体在观测点处产生的重力异常值为 g，根据万有引力定律， g 的表达式为

其中(ξ，η，ζ) 为地质体体积元的坐标，(x，y，z)为观 测点的坐标，G为常数，等于6.672×10^-11m³/(kg·s²).

计算中假设模型误差 C _m与先验估计误差 C _d都是相互独立且大小相等，σ_d、σ_m分别为第一个观测数据和第一个参数先验估计的方差， V _d、 V _m为单位矩阵，则有 C _d=σ^2d V ^d，C _m=σ^2m V _m.引入权重系数β，令β²= σ^2d/ σ^2m，目标函数E(m)可以改写为

方程(12)的求解实质是加权二乘拟合，关键是在重力观测数据拟合(第一项)和先验估计(第二项)之间寻求一种折衷，其中权重系数的选取非常重要.权重系数越大，反演解越接近先验估计；反之，反演解越接近于不稳定的初始模型.常用方法有L-曲线法(Hansen，1992)和广义交叉法(Wahba，1990)，详情 参考文献(Farquharson and Oldenburg，2000；Haber and Oldenburg，2000).

有了以上算法基础，以离散多面体为反演单元设计三维重力反演，包括：先验模型的确定、初始模型的建立、重力反演迭代与初始模型的升级、模型的三维显示，如图 2所示.详情以下展开.

图 2

Fig. 2

图 2 三维重力反演流程图 Fig. 2 3D gravity inversion flowchart

收集已有信息、简化地质单元、测量岩石物性、电磁或地震剖面以及钻测井资料的综合解释等，定义建模区域.

朴素贝叶斯分类(Naive Bayes Classifier，或NBC)是贝叶斯反演理论中的一种简单、稳定的分类算法，其思想是：对于给定的待分类项，求解在此项出现的条件下各个类别出现的概率，并按照概率最大划定待分类项的类别.建立朴素贝叶斯分类模型所需估计的参数很少，对缺失数据不太敏感，算法也比较简单，运用该方法对先验信息进行分类.考虑到通过朴素贝叶斯分类建立初始模型，无需了解模型参数之过程状态序列，只需计算模型参数最终状态的概率函数，因此，采用隐马尔可夫模型(Hidden Markov Model，或 HMM)更合适.设计如下主要流程：

第一步，确定模型参数特征属性及划分：根据先验信息p确定出模型参数矢量 m .其中，模型参数矢量包括模型中各地质体的位置、埋深、厚度、密度等物理参数及其可能变化的范围.

第二步，获取所建模型体对应的观测向量数据作为训练样本：根据先验地质、地球物理信息和建模区域确定观测数据与地下物理模型参数之间的函数关系，即参入反演的核函数 g .再由公式 d = gm，通过核函数 g (见公式11)与模型参数矢量 m 正演计算出所建模型体的观测向量数据 d .

第三步，计算训练样本中每个类别的频率：将模型体生成的观测向量数据划分为不同的类别.

这里是根据先验信息推算模型区域每个地质体的最大似然解 m _ML.对于重力贝叶斯反演，如公式(12)最大似然模型 m 的求解，涉及地质体坐标位置(ξ，η，ζ)、剩余质量(σ)和重力异常值(g)五个参数，计算维数较高，计算成本也高，而且影响反演解的收敛性增加不确定度.因此，借助模型降阶技术采用固定部分参数减小计算维数的方式，固定模型中地质体的坐标位置使得类别参数降到最少.这样，当体积确定时，根据“质量=密度×体积”，通过分段划分地质体与区域围岩的相对密度变化区间，等同于划分剩余质量的变化区间，再根据万有引力公式就可计算出模型地质体的重力异常数值范围.

第四步，计算这一位置固定的地质体的相对密度可能概率.

第五步，使用分类器鉴别所述模型地质体的最大似然解并将其分类.按照上述最大似然解算法，依据公式(12)计算模型，要求模型参数最大限度满足重力观测值与先验信息这一分类原则，确定模型中各个地质体(岩体)的坐标位置、埋深、厚度、密度等物理参数.

第六步，根据第五步的分类结果确定地质模型.经过分类器的鉴别生成地质模型，并计算模型重力与实际观测值的拟合方差.由万有引力公式正演算出地质模型的重力异常值，再计算它与实际重力异常观测值的方差.假设反演后的最大似然模型为 ，数据残差为，对于静态数据误差可通过以下方式计算反演数据协方差矩阵(Dosso and Wilmut，2006)

其中.对于Gaussian数据误差分布假设，归一化的误差应该随机分布在以0为中心、协方差为1的周围，理论值为G(0，1).计算到方差值小于设定的0~1之间的某个正数ε(ε为方差门槛值)时，将该模型参数作为反演结果.否则，修正模型地质体的特征属性划分，如：地质体与区域围岩的相对密度区间或地质体的坐标位置，再重复上述步骤，直到方差符合要求为止.方差门槛值应根据建模区域大小确定，一般地，10 km²面积区域，方差门槛值ε设为0.001-0.03.

计算过程如图 3所示，概括为：根据已知的先验信息选择初始模型；根据所选择的初始模型确定重力异常训练样本；计算训练样本中每个类别的频率(归结为计算模型中地质体剩余质量的变化区间)；通过贝叶斯分类器对观测向量数据进行分类鉴别(运用模型地质体的最大似然解)，确定初始模型参数的逼近区间.若逼近区间符合模型重力与观测重力的拟合方差要求，则输出反演结果；若该逼近区间不符合拟合方差要求，则根据先验信息重新选择初始模型.通过上述六步，就可针对性地给出模型中各地质体的密度值和分布范围，得到较为准确的初始地质模型.具体程序见“基于贝叶斯反演理论的重力建模软件”(刘彦等，2015).

图 3

Fig. 3

图 3 贝叶斯算法程序 Fig. 3 Flowchart based on bayesian classification method

对每一个赋予初始密度的模型体，采用人机交互法修改建立2D剖面模型，完成建模区域内所有2D剖面的重力模拟后，将每条2D剖面的模型走向方向长度延长同剖面长度，初始模型变为2.5D地质模型.再对2.5D模型进行重力曲线拟合、模型修改处理，直到达到满意的数据拟合和获得合理的地质模型.可采用如澳大利亚Encom Technology Pty Ltd.公司开发的专业软件Encom ModelVisionPro^TM等重力数据处理、反演软件完成(祁光等，2012).

得到系列2.5D地质模型后，按照剖面的空间顺序依次将2.5D模型拼合成3D模型，计算3D模型的理论异常，并与实际异常对比，拟合误差较大的区域，需返回到2.5D剖面重新修改，如此调整得到3D地质模型.

采用三维可视化模拟软件(如：GOCAD或Encom PA)对3D地质模型进行可视化和结构解译处理，经过可视化及结构解译处理后的模型能够形象直观地反映地下三维密度结构.

为验证上述贝叶斯算法开展重力反演的有效性，选择单位球体和任意形态几何体两种模型体开展模拟实验.重力值可依据牛顿万有引力定律计算，球体的计算公式为(曾华霖，2005)

其中Δg为球体的重力异常值，其单位取μGal，1 μGal=10^-8 m·s^-2，G为万有引力常量，取值6.672×10^-11 m³/(kg·s²)，M为场源质量，h为场源中心深度或埋深，R为球体半径，X为观测平面上测点到场源中心的水平距离，这里指测线距离，h、R、X的单位均为m，σ为模型体相对围岩的密度差或相对密度，单位为kg·m^-3.对于任意形态几何体，可通过三角剖分法将其划分成若干个极小单元的四面体，根据这些最小单元四面体底面中心点的坐标(x，y，z)与三角形面积计算每个四面体的体积，再与密度乘积得到四面体的质量.为简化计算，任意形态几何体剖分至百个数量级单元体，剖分得到的小四面体近似于单位质量球体，重力异常值可采用公式(14)来计算，最后对所有剖分体的重力异常值求和，得到任意形态几何体的重力异常总值.

以单位球体为例，测点到场源中心的水平距离X和埋深h为固定值，采用上述贝叶斯算法，计算球体的密度范围.假设单位球体处于测点坐标X₀=1000 m处，埋深h₀为800 m(如图 4a所示)，球体与围岩的相对密度(即密度差)σ为30 kg·m^-3，测线总长2000 m，测点间距为10 m，加入正态分布的观测噪声，依据公式(14)正演计算得到重力观测值，如图 4b所示.

图 4

Fig. 4

图 4 单位球体数值模拟示意图 (a)单位球体；(b)加入正态分布噪声正演得到的重力观测值. Fig. 4Unit sphere numerical simulation modeling (a) Unit sphere sketch map; (b) Gravity observations by forward calculation with normal distribution noise.

采用设计的贝叶斯算法程序经过10次迭代运算误差小于0.02，满足精度要求，得到模型的相对密度值为31.7578125±5.2734375 kg·m^-3，运算过程如表 1所示.再采用隐马尔科夫链分类器鉴别 的反演计算，得到模型与围岩的相对密度值为 30.3107438614098 kg·m^-3，方差为1.87447932360261×10^-5，更接近实际密度差30 kg·m^-3.可见，隐马尔科夫链分类器鉴别的贝叶斯反演结果相对更准确.

表 1 (Table 1)

表 1 单位球体模拟贝叶斯计算结果 Table 1 Bayesian calculation results of unit sphere simulation 迭代 次数 相对密度 上限 (g·cm-3) 贝叶斯反演 概率 相对密度 下限 (g·cm-3) 贝叶斯反演 概率 计算 方差 0 1.83 0.334920234 -0.87 0.665079766 1 0.48 0.669840468 -0.87 0.330159532 2 0.48 0.339680937 -0.195 0.660319063 3 0.1425 0.679361873 -0.195 0.320638127 4 0.1425 0.358723746 -0.02625 0.641276254 5 0.058125 0.717174599 -0.02625 0.282825401 6 0.058125 0.436608921 0.0159375 0.563391079 7 0.03703125 0.579457701 0.0159375 0.420542299 8 0.03703125 0.512534703 0.026484375 0.487465297 ＜0.02 最终 密度 区间 0.0317578125±0.0052734375(g·cm-3)

表 1 单位球体模拟贝叶斯计算结果 Table 1 Bayesian calculation results of unit sphere simulation

选择与围岩的相对密度σ为30 kg·m^-3的任意形态几何体正演计算重力数据.为贴近实际，几何体定义为长轴a等于300 m、短轴b等于30 m的透镜体，三角剖分后如图 5a所示.透镜体位于测点坐 标X₀=1000 m处，埋深h₀为800 m，测线总长2000 m，测点间距为10 m.加入正态分布的噪声，正演计算得到的重力值如图 5b所示.

图 5

Fig. 5

图 5 任意形态几何体数值模拟示意图 (a)任意形态几何体三角剖分；(b)加入正态分布噪声正演得到的重力观测值. Fig. 5 Any geometry numerical simulation modeling (a)Triangulation schemes of any geometry;(b)Gravity observations by forward calculation with normal distribution noise.

同样，采用上述贝叶斯算法程序计算，经过10次迭代运算误差小于0.02，满足精度要求，得到相对密度值为31.7578125±5.2734375 kg·m^-3，计算结果与单个球体相同.可见多个球体叠加运算对于贝叶斯算法影响并不大，同时也证明贝叶斯算法具有一定的稳定性，对梯度方向没有要求.采用隐马尔科夫链分类器鉴别反演计算，得到剩余密度值为 30.0626173207599 kg·m^-3，方差为1.72748273715972×10^-5，结果与实际给定的相对密度30 kg·m^-3非常接近，说明本文设计的贝叶斯反演方法对复杂形体同样有效.

泥河矿床位于安徽省庐枞火山岩盆地西北边缘，是长江中下游地区实施第二空间找矿以来在庐枞矿集区500~1500 m深度范围内最先发现的大型“玢岩型”铁矿床.矿区处于罗河—黄屯北东向成矿带上，南西距罗河铁矿床3 km，北东距龙桥铁矿床13 km(吴明安等，2011).泥河深部矿体的发现引起众多的关注，由此产生一系列关于“玢岩型”铁矿床成矿与控矿模式以及如何区分矿体与非矿体异常的思考.欲解决这些关乎找矿成败的关键问题，需要精细解剖地球物理异常性质，深入认识矿区深部结构，明晰区域地质状况.

早在2007年初，高锐等在庐枞矿集区及外围布设两个剖面总长达153.28 km的十字形深地震反射测线(高锐等，2010)，2008年吕庆田等部署两条穿过泥河、罗河、大包庄等主要矿体的10 km长的地震剖面测线(吕庆田等，2010)，2009—2010年吕庆田等又在庐枞地区完成5条总长达326 km的深地震反射剖面(吕庆田等，2014)，这些反射地震成像很好地显现出矿区所在区域的深部构造.不仅如此，针对泥河矿床更是实施完成不少具体工作：吴明安(2011)总结了泥河铁矿的发现过程及区域找矿方向，赵文广(2011)查明矿床的地质特征并分析了成因，认为辉石闪长玢岩是泥河铁矿的成矿母岩，正长斑岩穿切火山岩地层和矿体，是成矿期后形成的脉岩；刘彦(2012)、祁光(2012)以泥河矿床为例开展了基于重力异常分离方法寻找深部隐伏铁矿的实践以及地质信息约束下的三维重磁建模工作；张昆(2014)开展了音频大地电磁、可控源音频大地电磁、瞬变电磁和复电阻率四种电磁法探测试验，获得泥河铁矿区三维电性结构模型；范裕(2014)开展了矿床成岩成矿年代学研究，认为泥河铁矿床的成岩成矿作用几乎同时发生，盆地内130 Ma左右的辉石闪长玢岩侵入体是寻找泥河式玢岩型铁矿床的勘探靶区；张舒(2015)分析测试了泥河铁矿主成矿期矿石矿物的稀土元素、硫同位素及铅同位素，认为泥河铁矿是由次火山岩体演化产生的含矿高温热液在闪长玢岩穹窿顶部，通过交代充填作用形成的玢岩型铁硫矿床；张明明(2014)以钻孔资料为依据，采用离散光滑插值法建立成矿岩体的精细三维模型，重现闪长玢岩体的形态特征；周涛发(2014)则系统总结了泥河铁矿床的地质地球化学特征，提出泥河矿床的形成除了与辉石闪长玢岩有关外，还与三叠系膏盐层有着重要的成因关系，并建立起泥河铁矿床两期成矿模式.有了以上深入细致的工作，泥河矿区结构与成矿控矿特征渐趋明朗，积累了大量宝贵的资料，这些信息资料是三维重力反演试验成功的重要保证.

利用收集到的资料和数据如上述设计程序开展矿区三维重力反演，通过试验效果检验方法的效力.资料包括：1 ∶ 10000地形地质图，1 ∶ 50000高精度重力观测数据，74个钻孔累计8.33万米的钻井资 料(含69个钻孔测井数据)，882块钻孔岩芯物性测 试数据，以及庐枞矿集区5条电磁剖面和9条地震剖面.

首先简化地质单元，提取重要信息.将1 ∶ 10000地形地质图简化(见图 6)，泥河矿区存在第四系(Q)砂砾岩、白垩统杨湾组(K₁y)砂岩、白垩统双庙组(K₁sh)和白垩统砖桥组(K₁zh)火山岩4套地层.地层产状平缓，褶皱不发育，往深部地层产状略有起伏变化.构造简单，仅有火山岩地层的北西倾斜单斜构造和成矿期后的浅层断裂.需要重点关注的是，泥河铁矿体主要赋存在闪长玢岩体的顶部与砖桥组下段的火山碎屑岩中，闪长玢岩体的穹状隆起部位是赋矿的有利构造(吴明安等，2011；赵文广等，2011).

图 6

Fig. 6

图 6 泥河矿区地质简化图 Fig. 6 Simplified geologic map of Nihe mining area

接着备好数据，确定模型空间.重力资料采用1 ∶ 50000高精度观测数据(源自安徽省地质调查院)，需开展位场分离工作保留剩余重力值(刘彦，2012).反演区域是在综合矿区地层、岩体及整体结构并充分考虑到钻孔资料可明确的1200 m岩性深 度这些已知信息基础上划定的.建模2D剖面部署 如图 6方框内线条所示，平面面积为5.6 km²(2.8 km× 2.0 km)，剖面垂直构造走向取北偏西40°方向，间距 为100 m，共划分28条平行剖面，区域重力场的背景密度设为2.62 g·cm^-3.

运用贝叶斯反演程序建立初始模型.采用离散多面体建模需要依次对每条剖面开展重力反演工作.以矿区9线剖面为例，依据已有地质信息、钻孔资料以及电磁地震等认识，勾画出先验模型如图 7a所示，该先验模型基本界定地质体的坐标位置和埋深，这等同于固定地质体的体积.再计算各地质体的密度区间，采用刘彦等(2015)设计的贝叶斯反演程序.结合测试的882块物性数据，给出各地层和岩体的密度区间范围，利用程序快速算出初始模型中各 地质体的密度(见表 2)，省却常规离散体反演建模 过程中反复调试密度的困顿，可大大提高建模效率. 从地表起由上至下，第四系砂砾岩(Q)密度为2.10± 0.03 g·cm^-3，白垩统杨湾组(K₁y)地层密度为2.30±0.05 g·cm^-3，白垩统双庙组上段(K₁sh₂)地层密度为2.55±0.04 g·cm^-3，白垩统双庙组下段(K₁sh₁)地层密度为2.57±0.03 g·cm^-3，白垩统砖桥组(K₁zh)地层密度为2.65±0.03 g·cm^-3； 其它岩体及矿体的密度范围：闪长玢岩(δμ)岩体密 度为2.90±0.04 g·cm^-3，磁铁矿(Mt)、黄铁矿(Py)以及石膏矿(Ah)密度分别是3.48±0.09 g·cm^-3、3.35±0.09 g·cm^-3、2.93±0.07 g·cm^-3.

图 7

Fig. 7

图 7 泥河矿区9线剖面重力反演示意图 (a)剖面地层和岩矿体初始模型；(b)重力离散体反演2.5D拟合剖面. Fig. 7 9 line profile gravity inversionof Nihe mining area in Anhui province (a)Initial modelof profile stratum and rock ore body;(b)2.5D fitting section of gravity discrete simulation inversion.

表 2 (Table 2)

表 2 采用贝叶斯计算的泥河矿区初始模型中各地质体的密度 Table 2 Initial model density of Nihe mining area by bayesian calculation 地层/地质体名称 符号 密度区间(g·cm-3) 第四系砂砾岩 Q 2.10±0.03 白垩统杨湾组 K1y 2.30±0.05 白垩统双庙组上段 K1sh2 2.55±0.04 白垩统双庙组下段 K1sh1 2.57±0.03 白垩统砖桥组 K1zh 2.65±0.03 闪长玢岩 δμ 2.90±0.04 磁铁矿 Mt 3.48±0.09 黄铁矿 Py 3.35±0.09 石膏矿 Ah 2.93±0.07

表 2 采用贝叶斯计算的泥河矿区初始模型中各地质体的密度 Table 2 Initial model density of Nihe mining area by bayesian calculation

2.5 D模型的构建.有了确定的密度区间再调整修改地质体的形态，即进行重力的反演迭代与初始模型的升级，直到模型曲线和实际重力曲线拟合完好，此时模型生成的重力值与实际重力观测值方差满足允许范围，由此建立剖面的2.5D初始模型(如图 7b所示).同样方式，依次做好全部28条剖面的重力反演迭代与初始模型的升级工作.

3D模型的生成与显示.选用离散多面体反演建模专业软件Encom ModelVision Pro^TM 拼接拟合，完成整个泥河矿区的三维重力反演.反演最后得到的模型加载到三维可视化专业软件Encom PA，实现多方位成像和结构解译，最终模型如图 8所示.

图 8

Fig. 8

图 8 安徽省泥河矿区三维重力反演模型 Fig. 8 3D gravity inversion model of Nihe mining area in Anhui province

图 8所示的3D模型显示：泥河矿区内的主要矿体为磁铁矿、黄铁矿以及石膏矿；矿体整体呈北东向展布，延伸至东北部时矿体稍有抬升；矿区磁铁矿、黄铁矿较多，石膏矿相对较少；磁铁矿主要位于矿区的西南部，呈透镜状；黄铁矿主要集中在矿区的东北部，中部少量的石膏矿；黄铁矿和石膏矿埋深相 对于磁铁矿较浅，矿体深度大致在地下600~1100 m 范围内.矿区东北部矿体主要是垂直方向上的两层黄铁矿，上层矿体体积较小，呈层状分布，矿体在西南部较小，东北部较大，平均宽度约为245 m，埋深约为600 m，厚度约40 m；下层矿体体积较大，呈透镜状，埋深在800~1050 m，最大宽度约为680 m.三维重力反演结果还揭示出泥河矿床的控矿特征：矿区内地层由浅到深主要为第四系盖层(Q)、白垩统杨湾组(K₁y)、白垩统双庙组(K₁sh)和白垩统砖桥组(K₁zh)；较浅层的黄铁矿和石膏矿多位于砖桥组上段地层，体积较小，呈层状似层状分布；矿体主要集中在砖桥组下段侵入岩中，侵入岩主要为辉石闪长玢岩和脉岩；侵入体顶面形成隆起，矿区西南部隆起陡立，矿区东北部隆起宽缓.该模型刻画展示的矿体和岩体的空间关系，与吴明安(2011)、赵文广(2011)等运用地质学，范裕(2014)等运用成矿年代学，周涛发(2014)、张舒(2015)等运用地球化学，张明明(2014)等运用离散光滑插值法，张昆(2014)等运用电磁法等推论结果吻合一致，说明本方法开展三维重力反演是实用有效的.

(1)本文发挥离散多面体反演技术和贝叶斯方法优势，通过理论分析、综合研究，设计程序完成模型试验和矿区实践，改进了传统的三维重力反演方式，得到如下认识：将贝叶斯方法引入重力反演初始模型的建立中，可以定量地利用已知先验信息、观测数据以及物理规律来确定模型参数具体数值，节省模型选择时间，提高反演效率；最大似然函数算法应用于分类器的鉴别中改善了传统朴素贝叶斯分类方法，采取模型降阶技术固定所建模型中几何体的形态或密度，达到在几何体形态、密度和重力值五个参数中降低维数，可以减小高维不确定性和正演的计算量，有利于降低多解性，提高反演的准确度；任意形态几何体模拟试验以及泥河矿区反演实践证明，这种基于贝叶斯算法的三维重力反演方式，能够充分吸收先验信息，改善反演精度，提高建模效率.

(2)重点主要集中于改善反演初始模型的建立方式，三维重力反演的实现是通过本次设计的贝叶斯方法模块与其它成熟反演建模软件(如：Encom ModelVisionPro^TM)整合完成的.在贝叶斯反演计算中，可以通过固定密度、变化体积这种方式准确圈定出异常体的展布形态，这对于实际勘探生产将更具有应用价值.

致谢 本文得到中国地质大学(北京)薛涛老师及学生张永强的大力支持和帮助，审稿专家富有建设性的意见对于本文的完善起了很好的作用，稿件修改期间还得到中国地质科学院矿产资源研究所史大年研究员的指导以及与四川大学周永道博士的有益探讨，所有这些一并致谢！