贝叶斯推断下的胶西北金矿集区焦家断裂成矿构造不确定性建模

Uncertainty Modeling of Jiaojia Fault Metallogenic Structure Based on Bayes Inference in Northwestern Jiaodong Peninsula, Shandong Province

摘要: 焦家断裂带是胶西北金矿集区最重要的控矿构造之一，深部找矿潜力巨大。随着矿产勘查深度的增加，深部成矿构造研究已然成为深部找矿的基础前沿工作，是实现深部找矿三维结构分析与定量成矿预测的先行任务和前提条件。然而，焦家断裂深部地质信息获取难度大、数据信息量少，导致三维结构模型可能具有较大的不确定性，因此如何度量并且有效减小模型不确定性成为深部成矿构造研究需要解决的问题。本文提出将地质–地球物理信息融入三维结构建模，实现焦家断裂深部成矿构造模型的不确定性度量和可视化。研究表明，在地质约束的基础上进一步加入地球物理数据约束能够进一步减小三维结构模型不确定性，最终的不确定性模型为后续重建焦家深部成矿构造的三维精细结构提供理论基础，为焦家深部找矿提供引导。

Abstract: The Jiaojia fault is one of the most important ore-controlling faults in the Jiaodong Peninsula and has great potential for deep mineralization. With the increase of exploration depth, the study of deep mineralization structure has become the basic frontier work, which is the first task and prerequisite to the three-dimensional structure analysis and quantitative prediction for deep mineralization. However, the difficulty of obtaining geological information and lacking data information in depth in the Jiaojia fault leads to large uncertainty in the geological model. Thus, how to measure and effectively reduce the model uncertainty becomes the key point in the study of geological modeling. This paper proposes to integrate geological and geophysical information into 3D geological modeling for uncertainty measurement and visualization of the Jiaojia Fault. The result shows that the addition of geophysical constraints can further reduce the uncertainty of the 3D geological model. And the uncertainty model provides a theoretical basis for the subsequent reconstruction of the refined geological model of the Jiaojia fault and provides guidance for the Jiaojia mineralization.

1. 引言

通过矿床地质和找矿预测研究发现成矿流体容易在成矿构造交汇、拐弯、褶皱轴部、核部、岩体接触带等部位聚集 [1] [2] [3] [4] [5]，而随着预测深度增加，深部矿产预测更加依赖于对深部成矿构造的推断及其对矿化制约规律的认识 [6] [7] [8] [9]。虽然当前深部成矿构造的三维建模方法趋于成熟，但由于勘探深度的限制以及深部信息的缺失，深部成矿构造三维结构模型不确定性仍然存在，因此，深部成矿构造的三维不确定性建模分析成为开展深部成矿构造研究以及深部矿产资源预测工作中的重点。

焦家断裂带作为矿集区内主要控矿断裂之一，有很多大中型金矿如焦家金矿、望儿山金矿、新城金矿等沿断裂带分布 [10] [11] [12]。焦家断裂带累计探明金资源储量超过1200 t，并且还在不断增加，展现了深部重要的勘查和研究价值 [13]。为了构建焦家断裂深部成矿构造的不确定性模型，并且有效减小模型的不确定性，需要将多源、异构、多尺度的勘查数据融入建模过程，但是对于复杂多样且不同分辨率的输入数据，可能会导致建模结果的灵活性和再现性显著降低，并且难以考虑数据之间的分布概率 [14] [15]。由此发展起来的贝叶斯模型能够较好地顾及地质先验知识和观测数据之间的概率分布，补充在隐式建模插值步骤中无法考虑的额外地质信息，并对多源数据之间的关联进行自动评估，寻求先验知识和实测数据之间的最优拟合。此外，本文将采用马尔科夫链蒙特卡洛算法(MCMC)对贝叶斯模型的后验概率分布进行采样，以采样结果确定三维结构模型的后验概率分布，从而实现贝叶斯推断。基于MCMC算法的贝叶斯推断并不算新方法，前人在进行随机扰动及概率抽样的研究中也有广泛使用 [16] [17] [18]，蒙特卡洛积分从给定分布中抽样，构造样本均值来逼近期望，而在大型多维空间中，通常构建马尔科夫链产生样本的蒙特卡洛积分，其计算需要对高维概率分布进行积分，从而对模型参数作出推断或预测。马尔科夫链通过随机变量的在不同样本中的转移概率定义状态空间的随机探索形式，并且规定每一次迭代样本是否被接受只与上一次接受样本相关，不受之前接受样本的影响，且与其余的随机变量条件独立，以此获得相对独立的样本，充分提高采样效率，并以此较为直观、准确地获得三维结构的空间不确定性分布 [19] [20] [21]。

本文以焦家断裂带主断裂面为研究对象，在前人研究以及相关野外勘探数据资料的基础上，将地质信息以及地球物理重力数据融入贝叶斯推断过程，并结合马尔科夫链蒙特卡洛随机算法进行焦家断裂深部成矿构造的不确定性建模，为后续的精细建模以及深部找矿工作提供参考。

2. 研究区地质概况

胶西北金矿集区地处华北板块东南缘，北邻渤海湾，南为胶莱盆地，西部为郯庐断裂沂怵段，东部以五连－烟台为界与苏鲁超高压北延部分相隔(图1)。胶东地区断裂构造非常发育，而焦家断裂为胶西北金成矿最密集的构造带 [22] [23]。焦家断裂带自南向北大致分为三段，分别为寺庄段、马塘–新城段、新城–高家庄子段，长约60 km左右，走向30˚~50˚之间，倾角为30˚~50˚，局部可达80˚，最宽处可达1000 m。该断裂带主要走向NNE-NE向，但变化较大，形态不规则 [24]。断裂带在演化过程中形成了有序的断裂带结构，分带性明显，上下盘分带基本对称，由于强烈的热液蚀变和后期的脆性变形叠加，许多断裂早期韧性变形的组构已被破坏，仅见有糜棱岩角砾。由于焦家金矿带金矿体主要赋存于主裂面下盘蚀变程度较高的蚀变岩以及区域规模大并呈缓倾斜的韧、脆性叠加断裂带中 [10] [12]，因此进行对焦家断裂深部构造研究对于找矿工作具有重大意义。

Figure 1. Geological map of the Jiaodong peninsula (modified after Song et al., 2019) [11]

图1. 胶东金矿集区地质简图(改自宋明春等, 2019) [11]

焦家矿集区勘探工作中，3000 m的最深见矿钻孔是目前最佳研究载体，前人已开展初步的地质特征描述，取得一些基本的认识 [25] [26] [27]。本研究涉及的焦家断裂带成矿构造三维不确定性模型基本包含整个焦家断裂带主断裂面，根据前人研究 [28] [29] 以及相关勘察资料，选定的研究区上下盘均以玲珑二长花岗岩为主，以焦家断裂带主断裂面为界，东侧(下盘区域)除玲珑二长花岗岩外分布有部分郭家岭花岗闪长岩，西侧(上盘区域)除玲珑二长花岗岩之外分布有大片的变辉长岩和小部分变辉长岩质碎裂岩。构建的研究区岩性结构模型如图所示(图2)。

Figure 2. Schematic diagram of the initial model and lithology of the Jiaojia fault

图2. 焦家断裂带研究区域初始模型与岩性分布模型示意图

3. 方法

本文进行了研究区资料处理，整理并分析了地质及地球物理可用信息，并采用隐式建模方法表示深部成矿构造，随后结合马尔科夫链蒙特卡洛(Markov Chain Monte Carlo, MCMC)算法，对初始模型的关键点进行扰动，从而获得一系列随机模型，并基于贝叶斯框架融合平滑约束、重力约束以及产状约束生成后验概率，以此构建MCMC接受率对随机模型进行采样，获得一系列采样模型后，结合信息熵评估和度量模型不确定性。

3.1. 三维隐式建模

为更好模拟深部成矿构造形态及位置，首先需要将三维模型地质空间进行数字化处理，即将整个地下研究空间进行离散化。假设研究区域为规则立体空间D，该空间完全位于地下，并且包含整个成矿构造断裂带。将立体空间分割为若干个个大小一致的小立方体(图3)，每个小立方体包含的信息包括：1) 中心点x坐标(x, y, z)，2) 立方体边长R，以及3) 岩性密度属性ρ，其中岩性密度属性需要根据中心点坐标相对于成矿构造的位置来进行赋值。

基于上述离散化研究空间，焦家主断裂面能够通过符号距离函数进行表示。对研究区域的地质剖面图进行矢量化，提取成矿构造主断裂面的线串模型，以此圈定断裂面的边界，构建深部成矿构造的显式地质模型，并计算该空间中每个立体单元的中心点到线串的最短距离，以此构建隐函数，该隐函数 $f\left(x\right)$ 具有如下特征：

$f\left(x\right)\{\begin{array}{ll}>0,\hfill & x\in 上盘区域\hfill \\ =0,\hfill & x\in 断裂面内\hfill \\ <0,\hfill & x\in 下盘区域\hfill \end{array}$ (1)

即在整个地下研究空间中，若离散化立方体x位于断裂面上盘区域，则符号为正，反之，点x位于下盘区域，符号为负，以此实现深部成矿构造的数字化表达。用上述隐式函数表达作为初始模型，以便后续进行迭代计算。

Figure 3. Schematic diagram of the discrete study area

图3. 离散空间划分示意图

3.2. 多源数据同化框架

贝叶斯理论是一种经典的统计推断方法，目前已有大量研究学者将此类统计理论应用于融合多源数据的复杂分布概率 [30] [31]。在本文案例的研究中，使用贝叶斯框架下的地质先验信息和地球物理重力似然函数。地质先验信息使构建的地质模型更加符合地质认知，从而得出合理的地质解释，重力似然函数考虑地球物理勘探中的数据样本信息，使模型更加具有可靠性，保证深部成矿构造模型更加接近地下真实形态。

在本文中，对于给定的地球物理观测重力异常数据集合 $\epsilon =\left\{{\epsilon }_1,{\epsilon }_2,\cdots ,{\epsilon }_N\right\}$ 和推断的目标地质界面f，贝叶斯后验概率 $P\left(f|\epsilon \right)$ 表示为：

$P\left(f|\epsilon \right)=\frac{P\left(f\right)P\left(\epsilon |f\right)}{P\left(\epsilon \right)}$ (2)

其中， $P\left(f\right)$ 为先验项，用以表达深部成矿构造的先验信息，从地质的角度出发约束迭代模型的形态及产状。 $P\left(\epsilon |f\right)$ 为似然项，用以表达在迭代过程中重力计算异常值与实测值之间的相关性， $P\left(\epsilon \right)$ 为一恒定的常数。

结合重力观测数据的似然函数 $P\left(\epsilon |f\right)$ 可以表示为：

$P\left(\epsilon |f\right)={\displaystyle {\prod }_{i=1}^{N}P\left({\epsilon }_i|f,\theta \right)}$ (3)

其中， $P\left({\epsilon }_i|f,\theta \right)=\frac{1}{\sqrt{2\text{π}}\sigma }\exp \left(-\frac{{\left({\epsilon }_i-fw{d}_{f,\theta }\right)}^2}{2{\sigma }^2}\right)$，体现了重力正演值与观测值之间的差异。 $\sigma$ 为常数， ${\epsilon }_i$ 表示重力异常观测值， $fw{d}_{f,\theta }$ 表示对应观测点的重力正演值。由于断裂面上下盘可能具有不同的岩性，因此剩余密度并非一个恒定的值，并且同一种岩性内部密度也有所变化，此时需对断裂面研究区进行分区最优密度反演。

地质先验知识完全是根据人们的经验和认知来进行判断，对先验信息的认知程度会对后验概率的计算产生一定的影响。这样可以使得通过该方法得到的三维模型能够更好的符合专业的地质认知，并可以得出合理的地质解释。本文中，先验项主要融合了地质约束，如模型一致性检验、断裂面法向量先验以及断裂面平滑项，对于所有先验项 $P\left(f\right)$ 都可以以下形式进行表达：

$P\left(f\right)=\frac{1}{Z}\exp \left(-E\left(f\right)\right)$ (4)

其中，Z为常数，用于表达先验项的权重， $E\left(f\right)$ 为不同先验项的能量项，在本文中使用高斯分布进行表达。

通过上述先验信息以及重力约束似然项的计算，可以构建出贝叶斯框架下的后验概率，利用概率形式表现模型的附加约束信息，成功实现地质–地球物理数据同化，确保模型的准确性和可靠性，启发后续工作在模型模拟数据与实测观测之间寻找最优解过程中的自动化，为深部成矿构造结构推断及不确定性度量提供理论支撑。

3.3. 基于MCMC采样的贝叶斯推断

由于贝叶斯理论在高维离散型问题中计算复杂且积分求解困难，难以使用解析方法获得后验分布，本文结合MCMC算法 [32]，使初始模型在一定约束下进行随机扰动，实现在希尔伯特空间中寻找深部成矿构造模型的概率分布，并且利用贝叶斯框架下结合了地–物信息的后验概率构建MCMC接受率，从而在扰动模型中进行采样，探索所有可能的断裂面模型，由此实现贝叶斯推断过程(见图4)。

MCMC采样的基本思路为对给定工程控制线串设置一个允许范围，使线串上的每一个关键点进行随机扰动，并结合贝叶斯框架中的后验概率采样出下一个符合要求的样本，且每一次样本只与前一次样本有关。通过每次采出的样本，设置接受率，再根据给定的接受率来判断当前样本是否被接受，如果是，则作为下一次采样的起始样本。利用MCMC算法对扰动后的模型构建接受率 $\alpha$ 表示为：

$\alpha =\min \left(1,\frac{\text{π}\left(z^{\ast }\right)q\left(z|z^{\ast }\right)}{\text{π}\left(z\right)q\left(z^{\ast }|z\right)}\right)$ (5)

其中，z为前一次迭代模型， $z^{\ast }$ 为当前扰动获得的模型， $\text{π}(\cdot )$ 为贝叶斯推断的后验概率， $q(\cdot )$ 为当前采样与前一次采样的转移概率，且服从高斯概率分布。此外，由于马尔科夫链需要经过多次状态转移才能达到稳定，因此对之前的采样进行舍弃，使采样结果更加接近真实结果，舍弃的样本数需根据不同约束条件进行确定。

Figure 4. Schematic diagram of the disturbed process

图4. 扰动示意图

3.4. 模型不确定性结果可视化

本文将使用信息熵结果对焦家断裂带模型的不确定性结果进行可视化。信息熵概念最初是在信息理论领域发展的，可对信息进行一定程度的量化，后来又被应用于地质建模的不确定性定量评估 [33]。在地质领域，信息熵可以衡量地质模型区域的离散化属性有关的信息变化量，如果将模型研究区域离散化为若干个小立方体，则每个体元的熵值 $H\left(x\right)$ 可以写为：

$H\left(x\right)=-{\displaystyle {\sum }_{i=1}^{N}p\left(x_i\right)\log \left(p\left(x_i\right)\right)}$ (6)

其中，N表示区域中体元的总数， $p\left(x_i\right)$ 表示第i个体元 $x_i$ 相关信息变化的概率。

为计算整个断裂面研究空间的每个离散立方体单元 $x_i$ 处的信息熵值，设置位置指示函数为 $I_f$ 为：

$I_{f_{up}}\left(x_i\right)=\{\begin{array}{ll}1\hfill & x_i\in f_{up}\hfill \\ 0\hfill & x_i\in f_{down}\hfill \end{array}$ (7)

$I_{f_{down}}\left(x_i\right)=\{\begin{array}{ll}1\hfill & x_i\in f_{down}\hfill \\ 0\hfill & x_i\in f_{up}\hfill \end{array}$ (8)

其中 $I_{f_{up}}$ 统计了离散单元在多次采样中位于上盘的次数， $I_{f_{down}}$ 统计了离散单元位于下盘的次数。根据每个小立方体单元隐函数符号判断其位置，并对每个单元的位置指示函数进行赋值。

由MCMC算法得到的n个断裂面线串模型，可计算每个立体单元的位于上盘( $p_{f_{up}}$ )或下盘( $p_{f_{down}}$ )的信息概率：

$p_{f_{up}}\left(x_i\right)={{\displaystyle \sum }}^{\text{}}\frac{I_{f_{up}}\left(x_i\right)}{n}$ (9)

$p_{f_{down}}\left(x_i\right)={{\displaystyle \sum }}^{\text{}}\frac{I_{f_{down}}\left(x_i\right)}{n}$ (10)

位置信息概率表示每个小立方体单元位于上下盘的出现概率，由此计算出每个小立方体单元的信息熵值，根据计算值，信息熵较大的区域不确定性较高，而熵值较小的区域不确定性较低(图5)，由此可视化整个断裂面模型的不确定性。

Figure 5. The relationship of the fault surface and the information entropy of cubes

图5. 断裂面位置与离散单元信息熵关系示意图

4. 不确定性建模结果

4.1. 仅平滑约束下的不确定性模型

在第一个示例模型中，我们将评估地质先验知识对焦家断裂带空间不确定性分布的影响(见图6)。在此只使用平滑数据，包括初始模型一致性信息以及几何平滑信息，对三维地质模型进行贝叶斯推断。

这个例子说明了焦家主断裂面在地质平滑先验认知约束下，将其三维结构模型的不确定性圈定于一个范围，为后续的贝叶斯推断提供基础。

Figure 6. Uncertainty model of the Jiaojiafault surface with smoothing constraint only: (a) spatial information entropy distribution; (b) profile of non-zero information entropy; (c) information entropy of six different segments in the study area

图6. 仅平滑约束下的焦家主断裂面空间不确定性分布模型：(a) 空间信息熵分布；(b) 信息熵非零部分的空间剖面分布；(c) 研究区内六个不同分段的信息熵分布

4.2. 地质约束与重力约束下的不确定性模型

在上述地质先验知识的基础上加入研究区范围内1:20万地球物理布格重力异常数据，以观察焦家断裂带主断裂面空间不确定性分布的变化趋势。

在地质平滑约束的基础上加入重力约束后，通过焦家断裂带主断裂面三维结构模型的信息熵分布，可以看出在加入布格重力异常数据后近地表区域信息熵值变化不大，但在断裂面深部信息熵深部有较大的变化，且均有变小的趋势(图7)。由此可总结出在贝叶斯推断中加入地球物理重力勘探数据能够在一定程度上减少模型不确定性，寻找与实际观测最优拟合的三维结构模型，获得的后验概率分布更具准确性和可靠性。

Figure 7. Uncertainty model of the Jiaojiafault surface with smoothing and gravity constraints: (a) spatial information entropy distribution; (b) profile of non-zero information entropy; (c) information entropy of six different segments in the study area

图7. 地质约束与重力约束下的焦家主断裂面空间不确定性分布模型：(a) 空间信息熵分布；(b) 信息熵非零部分的空间剖面分布；(c) 研究区内六个不同分段的信息熵分布

4.3. 附加产状约束下的不确定性模型

现在我们评估附加地质约束数据(如产状数据)将如何改变三维结构模型的后验概率分布。在本文中选择法向量数据作为附加地质约束信息，作为贝叶斯先验分布以约束模型推断。法向量是通过三维隐式模型结合索贝尔算子计算得到，在一定程度上能够代表焦家主断裂面的法向方向。图8展示了焦家断裂带主断裂面在附加地质约束后不确定性的空间分布(图8)。

将附加地质约束进一步融入贝叶斯推断过程后，焦家断裂带主断裂面的不确定性分布进一步减小，而由于输入数据在近地表区域的精度较高，因此在不断加入约束数据后，贝叶斯推断的不确定性结果无较大的变化。

Figure 8. Uncertainty model of the Jiaojia fault surface with the additional dipping constraints: (a) spatial information entropy distribution; (b) profile of non-zero information entropy; (c) information entropy of six different segments in the study area

图8. 附加法向量约束下的主断裂面空间不确定性分布模型：(a) 空间信息熵分布；(b) 信息熵非零部分的空间剖面分布；(c) 研究区内六个不同分段的信息熵分布

5. 结论

1) 本文建立了胶西北金矿集区焦家断裂带主断裂面的完整三维结构不确定性模型，为该地区三维精细结构建模提供理论支撑。

2) 焦家断裂带主断裂面不确定性模型能够较好地观察到焦家断裂带深部走向以及陡缓趋势变化，并且推断过程中融合了地质–地球物理信息，也进一步保证了模型的准确性，对焦家区域的深部成矿预测与找矿勘探工作提供参考。

3) 本文所使用的贝叶斯推断方法不仅结合多源数据，还基于MCMC采样对断裂面形态的后验概率分布进行推断，因此该方法未来也有望在加入更多约束数据的基础上，推广应用于其他深部成矿构造的三维不确定性建模工作，并且获得更加不确定更小的深部成矿构造三维结构模型。