一种获取地下三维密度结构的方法

摘要

本发明公开了一种获取地下三维密度结构的方法，涉及地球物理勘探领域，解决了现有技术无法准确获取地下三维密度结构的技术问题。本发明的主要技术方案为：确定建模区域的纵坐标、横坐标、深度及二维2D剖面间距；根据先验信息及建模区域，利用朴素贝叶斯分类方法进行重力反演建模，得到初始地质模型；采用离散体模拟方法将所述初始地质模型构建成2D地质模型；将所述2D地质模型拼接拟合成3D地质模型；对所述3D地质模型进行可视化和结构解译处理，得到地下三维密度结构。本发明主要基于贝叶斯分类反演原理，结合离散体模拟技术，实现重力密度模型的准确建立，提高所获取的地下三维密度结构的准确度。

权利要求

1.一种获取地下三维密度结构的方法，其特征在于，包括如下步骤：确定建模区域的纵坐标、横坐标、深度及二维2D剖面间距；

根据先验信息及建模区域，利用朴素贝叶斯分类方法进行重力反演建模，得到初始地质模型，具体为：根据所述先验信息确定出模型参数矢量；其中，所述模型参数矢量包括：地质体的埋深、厚度、密度及波速度；

根据所述建模区域、先验信息确定出核函数；其中，所述核函数表示观测数据与模型参数之间的函数关系；

根据所述模型参数矢量及核函数确定出观测向量数据，所述观测向量数据构成训练样本；

将所述观测向量数据划分为不同的类别，并计算出每个类别的观测向量数据在所述训练样本中的频率；

根据所述核函数、每个类别的观测向量数据在所述训练样本中的频率，计算出模型区域中地质体的最大似然解；

使用分类器鉴别满足所述地质体最大似然解的模型参数；

根据所述分类器的鉴别结果生成地质模型，并计算出所述地质模型的重力异常值和实际重力异常观测值的拟合方差；

其中，当所述拟合方差＜设定的正数ε时，由分类器的鉴别结果反演出的地质模型参数作为反演结果，得到初始地质模型；

当所述拟合方差≥设定的正数ε时，修正模型参数矢量，重复上述步骤；

采用离散体模拟方法将所述初始地质模型构建成2D地质模型；

将所述2D地质模型拼接拟合成3D地质模型；

对所述3D地质模型进行可视化和结构解译处理，得到地下三维密度结构。

2.根据权利要求1所述的方法，其特征在于，所述分类器为隐马尔科夫链分类器。

3.根据权利要求1所述的方法，其特征在于，所述设定的正数ε为0.001-0.03。

4.根据权利要求1所述的方法，其特征在于，所述先验信息包括地球物理信息和地质信息；

所述地球物理信息包括重力场信息、电磁信息、地震信息及岩石物理信息；

所述地质信息包括地层年代信息、钻孔资料信息及区域地质背景信息。

5.根据权利要求1所述的方法，其特征在于，所述采用离散体模拟方法将所述初始地质模型构建成2D地质模型的步骤，具体为：采用离散体模拟方法，结合先验信息对所述初始地质模型进行重力反演，得到2D地质模型。

6.根据权利要求1所述的方法，其特征在于，将所述2D地质模型拟合成3D地质模型的步骤，具体为：采用人机交互法对所述2D地质模型的每个2D剖面进行修改，叠加、拼接拟合出3D地质模型。

7.根据权利要求1所述的方法，其特征在于，采用GOCAD三维可视化模拟软件对3D地质模型进行可视化和结构解译处理，得到地下三维密度结构；或采用Encom PA地质资料分析、模拟软件对3D地质模型进行可视化和结构解译处理，得到地下三维密度结构。

说明书

一种获取地下三维密度结构的方法

技术领域

[0001] 本发明涉及地球物理勘探领域，尤其涉及一种获取地下三维密度结构的方法。

背景技术

[0002] 准确获取地下三维密度结构，对于了解和掌握地球深部地质体的分布形态非常重要。作为地球物理勘探的一种重要手段，重力三维建模技术可以通过测量区域的重力场信息计算出地下三维密度结构，进而推测地质体的分布形态。

[0003] 目前，重力三维建模技术主要是利用先验地质信息结合重力数据建立地下三维模型，进而获取地下三维密度结构。综观现有的这些先验地质信息约束下的重力三维建模技术，通常都包括三个部分：初始模型的建立、2D/3D重力反演模拟(2D指的是二维，3D指的是三维)及模型的3D显示与地质解释三个步骤。

[0004] 发明人发现重力三维建模技术中初始模型的准确建立能为后续的重力反演步骤奠定良好的基础；但是现有技术并未重视该步骤，只根据先验信息简单建立一个区域框架，准确度低，不仅给后续工作增加了不少工作量，还影响所获取地下三维密度结构的准确度。

发明内容

[0005] 有鉴于此，本发明实施例提供一种获取地下三维密度结构的方法，主要目的是通过准确建立初始模型，以提高所获取的地下三维密度结构的准确度。

[0006] 为达到上述目的，本发明主要提供如下技术方案：

[0007] 本发明实施例提供了一种获取地下三维密度结构的方法，该方法包括如下步骤：

[0008] 确定建模区域的纵坐标、横坐标、深度及二维2D剖面间距；

[0009] 根据先验信息及建模区域，利用朴素贝叶斯分类方法进行重力反演建模，得到初始地质模型；

[0010] 采用离散体模拟方法将所述初始地质模型构建成2D地质模型；

[0011] 将所述2D地质模型拼接拟合成3D地质模型；

[0012] 对所述3D地质模型进行可视化和结构解译处理，得到地下三维密度结构。

[0013] 前述的获取地下三维密度结构的方法，所述根据先验信息及建模区域，利用朴素贝叶斯分类方法进行重力反演，得到初始地质模型的步骤，具体为：

[0014] 根据所述先验信息确定出模型参数矢量；其中，所述模型参数矢量包括：地质体的埋深、厚度、密度及波速度；

[0015] 根据所述建模区域及先验信息确定观测数据与地下物理模型参数之间的函数关系，即参入反演的核函数；

[0016] 根据所述模型参数和核函数计算出观测向量数据，所述观测向量数据构成训练样本；

[0017] 将所述观测向量数据划分为不同的类别，并计算出每个类别的观测向量数据在所述训练样本中的频率；

[0018] 根据所述核函数、每个类别的观测向量数据在所述训练样本中的频率，计算出模型区域中地质体的最大似然解；

[0019] 使用分类器鉴别满足所述地质体最大似然解的模型参数；

[0020] 根据所述分类器的鉴别结果生成地质模型，并计算出所述地质模型的重力异常值和实际重力异常观测值的拟合方差；

[0021] 其中，当所述拟合方差＜设定的正数ε时，由分类器的鉴别结果反演出的地质模型参数作为反演结果，得到初始地质模型；

[0022] 当所述拟合方差≥设定的正数ε时，修正模型参数矢量，重复上述步骤。

[0023] 前述的获取地下三维密度结构的方法，所述模型区域中地质体的最大似然解的计算公式为：

[0024]

[0025] 其中，L(m)为模型区域中地质体的最大似然解；g(m)为正演算子，即上述核函数；d为观测数据矢量，m为模型参数矢量，T表示转置，下标G表示重力数据，下标S表示除重力数据外的其他已知地质数据。

[0026] 前述的获取地下三维密度结构的方法，所述分类器为隐马尔科夫链分类器。

[0027] 前述的获取地下三维密度结构的方法，所述设定的正数ε为0.001-0.03。

[0028] 前述的获取地下三维密度结构的方法，所述先验信息包括地球物理信息和地质信息；

[0029] 所述地球物理信息包括重力场信息、电磁信息、地震信息及岩石物理信息；

[0030] 所述地质信息包括地层年代信息、钻孔资料信息及区域地质背景信息。

[0031] 前述的获取地下三维密度结构的方法，所述采用离散体模拟方法将所述初始地质模型构建成2D地质模型的步骤，具体为：

[0032] 采用离散体模拟方法，结合先验信息对所述初始模型进行重力反演，得到2D地质模型。

[0033] 前述的获取地下三维密度结构的方法，将所述2D地质模型拼接拟合成3D地质模型的步骤，具体为：

[0034] 采用人机交互法对所述2D地质模型的每个2D剖面进行修改，叠加拼接拟合出3D地质模型。

[0035] 前述的获取地下三维密度结构的方法，采用GOCAD三维可视化模拟软件对3D地质模型进行可视化和结构解译处理，得到地下三维密度结构；或

[0036] 采用Encom PA地质资料分析、模拟软件对3D地质模型进行可视化和结构解译处理，得到地下三维密度结构。

[0037] 与现有技术相比，本发明实施例提出的获取地下三维密度结构至少具有如下有益效果：

[0038] (1)本发明实施例提出的获取地下三维密度结构的方法通过贝叶斯分类方法开展重力反演，准确地给出模型中各地质体的密度值和分布范围，得到准确的初始地质模型；而初始模型的准确建立，不仅使得后续的重力位场拟合更省时省力，关键是提高所获取的地下三维密度结构的准确度。

[0039] (2)本发明实施例提出的获取地下三维密度结构的方法，是基于离散体模拟的重力3D反演建模技术，改进了传统的离散体建模方法，吸收贝叶斯反演理论的优点，综合利用模型的先验信息、实测数据及物理规律等已有信息来反演模型参数，采取最大似然函数算法求解，得到模型的后验概率分布及其相应的特征量，运算中将隐马尔科夫链应用于重力反演分类器的鉴别中，准确地计算地质体的密度和分布范围，实现重力模型结构的精细刻画，提高所获取的地下三维密度结构的准确度。

附图说明

[0040] 图1为本发明实施例1提供的获取地下三维密度结构的流程图；

[0041] 图2为本发明实施例1提供的基于贝叶斯分类原理的重力反演流程图；

[0042] 图3为本发明实施例2提供的任意形态地质体三角剖分示意图；

[0043] 图4为本发明实施例2中任意形态地质体的重力正演计算结果示意图；

[0044] 图5为利用本发明实施例1的方法得到的泥河矿区矿体与地层的3D模型图。

具体实施方式

[0045] 为更进一步阐述本发明为达成预定发明目的所采取的技术手段及功效，以下结合较佳实施例，对依据本发明提出的一种获取地下三维密度结构的方法具体实施方式、特征及其功效，详细说明如下。

[0046] 重力三维建模技术可以通过测量区域的重力场信息，推算出地下三维密度结构，进而推测出地质体的分布形态。现今比较有效的重力三维建模技术主要是利用重力数据结合先验地质信息建立地下三维模型，进而获取地下三维密度结构。综观现有的这些先验地质信息约束下的重力三维建模技术，通常都包括三个部分：初始模型的建立、2D/3D重力反演模拟(2D指的是二维，3D指的是三维)及模型的3D显示与地质解释三个步骤。其中，初始模型的建立步骤、2D/3D重力反演模拟步骤在建模过程中起着至关重要的作用。发明人发现初始模型的建立这一步骤，目前尚未引起重视，还比较粗放。目前主要是根据先验信息简单建立一个区域框架，准确度低，不仅给后续工作增加了不少工作量，使所获取地下三维密度结构的准确度低，更影响地质体分布形态的推测效果。2D/3D重力反演模拟是对初始模型的进一步优化，并提供构建模型的物理和几何参数的空间分布。目前，2D/3D重力反演模拟主要有表面体模拟、网格反演模拟以及离散体模拟三种方法。发明人发现，表面体模拟方法可以较好地模拟地质体的构造信息，但重力场拟合误差较大；网格反演模拟方法容易获得模型的物性分布，模型产生的理论异常与实测异常的对比度好，但是网格反演很难加入先验地质信息的约束，模拟结果反映地质体的宏观分布尚可，局部细节往往与实际差距较大；而离散体模拟方法则具有上述两种方法的综合优势，既可以很好地反映地质信息，又可使观测场与理论场的拟合误差较小，它可以充分吸收先验地质信息(如地层倾向、断层和矿化体等)，还可以最大限度地发挥地质学家的经验和对区域地质的理解。因此，在2D/3D重力反演模拟步骤中，发明人认为离散体模拟方法比较实用。

[0047] 基于上述问题，本发明的发明人提出利用离散体模拟的重力3D反演建模技术获取地下三维密度结构的方法，该方法的核心是在建模区域内根据地质认识、建造若干个赋予密度值的二维地质剖面，即初始模型，再采用人机交互法逐个修改这些2D模型剖面，直到获得满意的数据拟合和合理的地质模型为止。其中，初始模型的准确建立是本发明的改进重点；本发明的发明人基于贝叶斯反演理论、结合离散体模拟技术，采用隐马尔科夫链和最大似然函数算法，精细判断出每个地质体的密度和分布范围，实现重力密度模型的准确建立，使得后续的2.5D/3D重力位场拟合更加省时省力。

[0048] 下面结合实施例对本发明作进一步的详细说明。

[0049] 实施例1

[0050] 本实施例提供一种获取地下三维密度结构的方法，如图1所示，具体包括如下步骤：

[0051] 1、定义建模区域：确定建模区域的纵坐标、横坐标、深度及二维2D剖面间距。

[0052] 该步骤主要根据所要研究的目标地质区域确定出建模区域。

[0053] 2、先验信息的预处理：

[0054] 如图1所示，该步骤中的先验信息包括地球物理信息和地质信息。地球物理信息包括重力场信息、电磁信息、地震信息及岩石物理信息。地质信息包括地层年代信息、钻孔资料信息及区域地质信息。

[0055] 该步骤主要收集目标地质区域的已有信息，简化地质单元，测量目标地质区域内的岩石物性，重力数据的编辑与处理，电磁或地震剖面的解释，以及钻井测井资料的解译等。

[0056] 3、初始模型的建立：根据先验信息、建模区域，利用朴素贝叶斯分类方法进行重力反演建模，得到初始地质模型。

[0057] 该步骤主要依据贝叶斯分类算法，准确给出模型中各地质体的密度值和分布范围。具体地，采用朴素贝叶斯分类方法开展重力反演，设计流程如图2所示，主要步骤如下：

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

[0059] 第二步，获取训练样本：根据先验地质、地球物理信息和建模区域确定观测数据与地下物理模型参数之间的函数关系，即参入反演的核函数gG和gS。再由公式d＝gm，通过核函数g与模型参数矢量m正演计算出表征所建模型体的观测向量数据d。

[0060] 这里，重力异常值与模型之间的函数关系有牛顿万有引力公式。当坐标原点位于地面，Z轴方向铅垂向下，X、Y轴位于水平面内，剩余质量为σ的地质体在观测点处产生的重力异常，根据万有引力定律，计算公式为：

[0061]

[0062] 式中，(ξ,η,ζ)为地质体体积元的坐标，(x,y,z)为观测点的坐标，G为常数，等于6.67×10-11m3/(Kg·s2)。

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

[0064] 该步骤中的类别归结为模型中每个地质体剩余质量的变化范围。即，根据模型地质体的分布范围将其划分为不同的变化区域段。该步骤具体为：对于模型中坐标位置固定的地质体，体积确定时，通过分段划分地质体与区域围岩的相对密度变化范围，根据“质量＝密度×体积”，等同于划分剩余质量的变化范围，再根据上述万有引力公式，也就计算出建模地质体的重力异常范围。与此同时，计算这一坐标位置固定的地质体的相对密度的可能概率，即计算每个类别在训练样本中的频率。

[0065] 第四步，计算每个类别条件下各个特征属性划分频率。即，根据先验信息推算建模区域每个地质体的最大似然解L(m)。

[0066] 该步骤中，模型区域中每个地质体的最大似然解的计算公式为：

[0067]

[0068] 其中，L(m)为模型地质体的最大似然解(即，模型区域中地质体的最大似然解)；g(m)为正演算子，d为观测数据矢量，m为模型参数矢量，T表示转置，下标G表示重力数据，下标S表示除重力数据外的其他已知地质数据。

[0069] 第五步，使用分类器鉴别所述模型地质体的最大似然解，并将其分类。

[0070] 根据上述最大似然解计算公式，要求模型参数最大限度满足重力观测值与先验地质信息这一分类原则，准确建立模型中的各个地质体(岩体)。

[0071] 较佳地，该步骤使用隐马尔科夫链分类器开展重力反演建模。即，在运算中将隐马尔科夫链应用于重力反演分类器的鉴别中，不需要了解模型地质体生成的状态过程，只需掌控转换的随机函数，由此设计程序完成反演计算，以有效改善重力三维建模方式。

[0072] 第六步，根据第五步的分类结果确定地质模型：根据所述分类器的鉴别结果生成地质模型，并计算出拟合方差。

[0073] 依据万有引力公式，由上述核函数，正演计算地质模型的重力异常值，再计算它与先验的实际重力异常观测值的方差(即，上述的拟合方差)。当这个方差值小于设定的某一个正数ε，则将该模型参数作为反演结果，进而得到初始模型。重力建模中，方差范围应根据建模区域大小确定，一般10km2地表面积建模区域，方差值ε设为0.001-0.03，此时重力异常拟合较好。否则，修正建模地质体的特征属性划分，如：相对密度范围或坐标位置，再重复上述步骤，直到符合要求为止。

[0074] 上述六步的实施可以参考流程图2，根据准备的数据、确定的模型空间、已知的先验信息，选择初始模型，根据所选择的初始模型确定重力异常训练样本(即：通过初始模型与核函数正演计算获取观测向量数据)，计算训练样本中每个类别的频率(归结为计算模型地质体剩余质量的变化范围)，采用贝叶斯分类器对观测向量数据进行分类鉴别(通过求解模型地质体的最大似然解)，计算初始模型参数的逼近区间，若该逼近区间符合拟合方差要求，则输出反演结果，得到初始地质模型。若该逼近区间不符合拟合方差要求，则根据先验信息，重新选择初始模型。

[0075] 通过上述六步，能准确地给出模型中各地质体的密度值和分布范围，得到准确的初始地质模型，这样不仅使得后续的重力位场拟合更省时省力，关键是提高了所获取的地下三维密度结构的准确度。

[0076] 4、采用离散体模拟方法将所述初始地质模型构建成2D地质模型。

[0077] 该步骤主要是2D模型的构建，实际就是带先验信息的重力反演，以对初始模型进一步优化，并提供构建模型的物理、几何参数的空间分布。

[0078] 该步骤可以采用Tikhonov正则化方法运算，具体如下：

[0079] 由于重力观测数据和建立的模型之间存在线性泛函关系，并可以表示为第一类弗雷德霍姆方程

[0080]

[0081] 这里Ek表示是第k个观测数据，Gk(x)是万有引力公式关系的核函数，m(x)是定义域为a≤x≤b的待求模型连续函数，N是观测数据的个数。用矩阵形式表示为：

[0082] Γα＝E

[0083] 其中，α＝(α1,α2,L,αN)T，E＝(E1,E2,L,EN)T，Γ是由核函数构成的内积对称正定矩阵，Γkn表示该矩阵的元素，由此可计算出α。

[0084] 根据先验信息通过引入正的加权函数ω(x)对部分模型m(x)进行约束，使得各处的m(x)拥有相等的权；构制模型，使解的欧几里德长度取极小，得到模型的解m(x)[0085]

[0086] 将上面矩阵方程解出的α代入，即可得到模型的解m(x)。

[0087] 5、将2D地质模型拟合成3D地质模型。

[0088] 该步骤主要是2.5D/3D重力位场的拟合，对每一个赋予初始密度的模型体，采用人机交互法对2D地质模型的每个2D剖面进行修改，直到获得合理的地质模型和满意的数据拟合为止。

[0089] 具体地，如图1所示，该步骤具体为，首先完成建模区域内所有2D剖面的重力模拟，根据Modelvision Pro重磁数据处理、反演软件，将每条2D剖面的模型走向长度Y方向缩短为剖面间距，使初始模型变为2.5D地质模型；再对2.5D模型进行重力曲线拟合、模型修改处理，得到系列2.5D地质模型；最后，按照剖面的空间顺序依次将2.5D模型拼合成3D模型，计算3D模型的理论异常，并与实际异常对比，拟合误差较大的地方，返回到2D剖面进行修改，如此将2.5D地质模型在3D环境中进行拟合、整合处理，得到3D地质模型。

[0090] 6、对所述3D地质模型进行可视化和结构解译处理，得到地下三维密度结构。

[0091] 该步骤主要采用GOCAD三维可视化模拟软件(或Encom PA地质资料分析、模拟软件)对3D地质模型进行可视化和结构解译处理，经可视化及结构解译处理后的数据能够直观地形象地获得地下三维密度结构。

[0092] 实施例2

[0093] 本实施例主要用于验证贝叶斯反演方法建立初始模型的可行性。具体方法如下：选择与围岩的相对密度σ为30kg·m-3的任意形态地质体开展正演模拟实验。

[0094] 为方便计算，对于任意形态地质体，采用三角剖分方法划分成若干个极小单元的四面体，将剖分到最终的每个四面体近似看作单位质量球体，根据这些最小单元四面体底面中心点的坐标(x，y，z)与三角形面积计算每个四面体的体积，再与密度乘积。每个四面体重力值的计算方法等效于球体计算，根据牛顿万有引力定律，重力值计算公式为：

[0095]

[0096] 其中，G＝6.67×10-11m3/(Kg·s2),M为场源质量，h为场源中心深度或埋深，R为球体半径，X为观测平面上测点到场源中心的水平距离，这里指测线距离，h、R、X的单位均为m，σ为模型体相对围岩的密度差或相对密度，单位为kg·m-3，Δg单位为μGal；1μGal＝10-9m·s-2。最后对所有四面体的重力值求和，由此求出任意形态地质体的重力值。

[0097] 为贴近实际，地质体定义为长轴a等于300m、短轴b等于30m的透镜体，地质体经三角剖分后如图3所示。假设地质体位于测点坐标X＝1000m处，埋深h为800m，测线总长为2000m，测点间距为10m。加入正态分布的噪声，正演计算得到的重力观测值如图4所示。先采用朴素贝叶斯反演程序计算，经过10次迭代运算，误差小于0.02，满足精度要求，得到相对密度值为31.7578125±5.2734375Kg·m-3。可见，采用三角剖分方法叠加计算复杂形态地质体对于贝叶斯反演运算影响并不大，同时也说明了贝叶斯反演解具有一定的稳定性，对梯度方向没有要求。采用隐马尔科夫链分类器鉴别再进行反演计算，得到剩余密度值为-3 -5
30.0626173207599Kg·m ，拟合方差为1.72748273715972×10 ，结果与实际给定的相对密度σ为30kg·m-3非常接近，说明设计的贝叶斯反演方法对复杂形体同样有效。

[0098] 实施例3

[0099] 本实施例为应用实例，具体将实施例1提供的获取地下三维密度结构的方法应用到实际地质勘查中。

[0100] 泥河铁矿是安徽省地质调查院在玢岩铁矿成矿模式和大型矿集区成矿理论指导下，利用钻探对重磁异常进行验证，于2007年5月首次发现的。泥河铁矿床位于安徽省庐枞火山岩盆地西北边缘，处于罗河-黄屯北东向成矿带上。南西距罗河铁矿床3km，北东距龙桥铁矿床13km，与罗河铁矿有着相同的成矿地质条件。泥河矿区内地层主要为白垩统砖桥组(K1zh)和白垩统双庙组(K1sh)火山岩，白垩统杨湾组(K1y)砂岩以及第四系(Q)砂砾土。矿区构造较为简单，主要为火山岩地层的北西倾斜单斜构造和成矿期后的浅层断裂。地层产状平缓，褶皱不发育，往深部地层产状略有起伏变化。据前人研究资料，泥河铁矿主要赋存在闪长玢岩(次火山岩)体的顶部与砖桥组下段的火山碎屑岩中，砖桥组下段的火山碎屑岩是寻找铁矿的地层标志。赋矿的有利构造为闪长玢岩体的穹状隆起部位，矿床成因类型为广义玢岩型铁矿。

[0101] 泥河铁矿床的发现绝非偶然，它的结构特征具体如何、成矿机理和找矿线索又是什么，对其它矿集区有着怎样的找矿启示？矿区除了钻探发现的700-1200m深度的铁矿体外，在其深部和外缘是否还有新的矿体存在？目前还存在的这些疑问，究其原因，在于矿区结构研究还不够精细。为了解决上述疑问，本实施例以1:1万重力地面测量资料为基础，采用实施例1提供的方法对泥河矿区开展重力三维建模，刻画该矿区地下三维密度结构，为解决上述疑问提供坚实依据。根据实施例1提出的方法得到的建模结果如图5所示。

[0102] 图5中的泥河矿区3D模型图清晰地刻画出岩矿体的三维形态和展布特征。如图5所示，泥河矿区内的主要矿体为磁铁矿、黄铁矿以及石膏矿；由图5所示的3D模型还可以看出：矿体整体呈北东向展布，延伸至东北部时矿体稍有抬升；磁铁矿以及黄铁矿的含量较高，石膏矿含量较少；磁铁矿主要位于研究区的西南部，呈透镜状；黄铁矿主要集中在矿区的东北部，中部少量的石膏矿；黄铁矿和石膏矿埋深相对于磁铁矿较浅，矿体赋存深度大致在地下
600m-1100m范围内。研究区东北部矿体以垂直方向上的两层黄铁矿为主，上层矿体体积较小，呈层状分布，矿体西南部宽度较小，东北部较大，平均宽度约为245m，埋深约为600m，平均厚度40m。下层矿体体积较大，呈透镜状，埋深在800m-1050m，最大宽度约为680m。

[0103] 图5所示的3D地质模型还揭示出泥河矿区的控矿特征：研究区内地层由上至下主要为第四系盖层(Q)、白垩统杨湾组(K1y)、白垩统双庙组(K1sh)和白垩统砖桥组(K1zh)。较浅层的黄铁矿和石膏矿多位于砖桥组上段地层中，体积较小，呈层状似层状分布；主要矿体集中在砖桥组下段侵入岩中，侵入岩主要为辉石闪长玢岩和脉岩。侵入体顶面形成隆起，矿区南西部隆起陡立，矿区北东部隆起宽缓。

[0104] 这些与其它方法(成矿年代学、地球化学等)推论结果一致，说明本实施例1提出的方法获取地下三维密度结构是有效的。

[0105] 综上，本发明实施例提供的获取地下三维密度结构的方法，是基于离散体模拟的重力3D反演建模技术，改进了传统的离散体建模方法，吸收贝叶斯反演理论的优点，综合利用模型的先验信息、实测数据及物理规律等已有信息来反演模型参数，采取最大似然函数算法求解，得到模型的后验概率分布及其相应的特征量，运算中将隐马尔科夫链应用于重力反演分类器的鉴别中，准确地计算地质体的密度和分布范围，实现重力模型结构的精细刻画，提高所获取的地下三维密度结构的准确度。

[0106] 以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。