一种基于地层相关结构约束的地质模型修正方法

Error Correction in Geological Model Based on Stratigraphic Interdependency

摘要

在地质勘探中，多数钻孔存在深部地层底板未采样的情况，不完整的采样信息限制了地质模型的准确性. 为了提高模型准确性，提出了一种基于地层相关结构约束的地层面模型修正方法. 由于地层的形成机制及后期构造运动，相邻地层的形态具有相似性. 基于此，利用Copula函数对相邻地层面的相关结构进行建模，构建相邻地层面的联合分布模型和待更新地层面的似然函数. 在贝叶斯框架中，利用似然函数对待修正的模型进行贝叶斯更新，得到地层面的后验分布，计算地层面的条件期望作为模型修正值. 利用该方法，对北海海岸带地区的钻孔数据和地层面模型进行了实验. 实验结果表明，模型修正后，地层面模型的误差降低，所提方法可以提高地质模型准确性.

关键词

地质模型 / 误差修正 / 地层相关结构 / 贝叶斯 / Copula

Key words

geologic models / error correction / stratigraphic interdependency / Bayes / Copula

0 引言

多数城市建设工程位于第四系沉积物覆盖区，为了保障工程建设和运营过程中的安全问题，有必要研究地质构造的几何形态和空间分布等信息，通过建立地层结构模型定量描述地质结构特征（谭飞等，2021）. 钻孔是建模过程中主要的一类约束信息，但是在地质调查过程中，多数钻孔的最深记录位置并不是地质界面，存在深部地层底板未采样的情况（见图1）. 基于这样的样本数据构建地质模型时，深部地层样本稀疏或接触位置测量不准确，导致地层的厚度和位置信息估计不准确，降低了模型的整体质量.

目前，多采用显式建模方法对第四系沉积层的结构和几何形态进行建模. 显式建模方法往往直接对地层面建模，主要思路为基于钻孔和剖面数据利用曲面拟合的方法构建地层结构的参数曲面（De Paor， 1996； de Kemp and Sprague， 2003； Sprague and de Kemp， 2005），或通过插值或者网格剖分的方法构建地层面的规则格网（Fremming， 2002）或不规则三角网（Lemon and Jones， 2003； Kessler et al.， 2009； Ming et al.， 2010； 张夏林等，2020）. 显式建模方法多基于三角剖分、参数剖分或空间插值，对各地层面单独建模，然后通过层与层之间的层序和接触关系，利用地层的拓扑关系将多层面模型组织在一起构成完整的地质体. 针对钻孔中深部地层底板位置不准确的问题，基于显式建模的方法通常采用以下几种方案构建底部地层面：（1）直接将钻孔的底界作为地层底板进行建模；（2）根据专家知识经验推断钻孔下方的地层位置信息，将钻孔延展到最深地层的底板位置；（3）舍弃钻孔中记录的最深位置信息，根据附近其他钻孔中记录的底板位置，利用插值方法计算当前孔位的地层面深度，该方法的准确性受附近样本信息的分布情况影响. 事实上，方法（1）得到的底板位置显然不准确，方法（2）虽然能保持地层的连贯性，但是会引入较大的主观不确定性，利用方法（3）创建的插值曲面与相邻地层面之间的关系往往与地层叠置原理不一致（Wellmann and Caumon， 2018），地质面会相互交叉或在域内留下间隙，这会违反地质模型的正确性条件（Caumon and Journel， 2004）.

对于存在问题的模型，还可以采用模型修正的方法使地质模型更接近真实情况.朱良峰等（2006， 2009）提出了基于地质剖面和虚拟钻孔的三维地层模型修正方法. 建模人员根据自身地质经验，通过增加地层面插值控制点的方式修改地质剖面，或增加包含地层位置信息的虚拟钻孔，利用修改后的剖面和虚拟钻孔约束，对模型中不符合要求的部分进行调整修正. 谢济仁等（2014）在建模过程中通过人工交互的方法编辑虚拟钻孔，调整地层的分层信息，迭代构建、调整地质模型直到模型符合认知. 修春华等（2015）通过巷道掘进过程中新布设的钻孔、工作面素描等，添加点、线、面数据修正克里金空间插值模型. 此类基于建模初始数据的模型误差修正方法往往需要在建模数据补充新信息（刘志锋等， 2012；侯征等， 2018），如新增钻孔、虚拟钻孔、辅助剖面、航磁等数据，再借助自动建模的方法实现模型的更新修正. 此外还有基于建模中间结果的模型误差修正方法，需要用户直接对地质模型的中间结果（吕迎红， 2007）进行人工修改. 许多研究使用地球物理正演模拟，然后根据正演模型和地球物理测量之间的视觉或定量比较手动调整地质模型（Gradmann et al.， 2013； Autin et al.， 2016； Haase et al.， 2017）. 此外还有一些基于信息整合的方法，例如Wellmann et al.（2017）在贝叶斯框架下将模型中的地层面点直接作为先验参数，使用地质知识和地球物理数据构建似然函数，估计地质模型参数的后验分布. 总体上，多数针对深部地层底板模型不准确的修正方法是利用新增数据（新增钻孔、剖面）或其他类型的数据（坡度、产状、重力、磁场、地震数据等）作为地质建模的约束项，或将专家知识经验（层厚、断层、褶皱的几何参数）转化为虚拟钻孔或剖面等数据约束，利用这些额外的信息约束重新计算模型，提升模型精度. 然而在很多时候，额外的数据难以获取，通过人工交互添加专家经验的方法耗时、费力、不易更新且会引入主观不确定性.

受地层沉积以及后期构造运动共同作用，新地层与老地层之间往往呈现出相似的空间形态特征（Tearpock and Bischke， 1991）. 因此，来自相邻地层的样本信息对地层面的估计能够起到约束. 针对深部地层底板位置信息不准确造成的模型误差问题，本文提出一种基于地层互相关的模型修正方法. 该方法首先利用不确定性分析方法得到待修正地层面模型的先验分布，然后在贝叶斯框架下，利用相邻地层间的互相关作为约束，构建模型更新的似然函数，对地层面先验分布进行贝叶斯更新，得到模型后验修正值. 利用该方法，本文以北海海岸带地区的钻孔数据和地层面模型为例进行实验，验证该方法的有效性.

1 方法

1.1　地层相关结构分析

1.1.1　地层面高程函数

为了定量地描述地质结构，首先需要对地质结构面进行数学化表达（Abrahamsen and More， 1994）. 地质建模领域通常利用钻孔等采样数据中包含的大量地层界面点信息，通过对界面点的插值拟合等方法构建地层界面的等值面描述地层的空间结构形态. 对于每个地层面 i，本文将位置 ( x , y )处的地层面 i的高程值定义为随机函数 Z_i ( x , y )，通过每个 ( x , y )位置的高程 Z ( x , y )描述地层面的空间几何形态. 这种基于 \left(X,Y\right)坐标的函数定义方法适用于沉积地层场景，需要指出，针对如倒转地层、透镜体、逆断层的地质界面在同一 \left(X,Y\right)位置会出现多个 Z值的复杂情况，可以采用多 Z值地质界面分块，满足每块区域的界面高程函数只对应一个 Z. 例如对于含有断层的区域，需要根据断层面将两侧划分成两个独立区域，对两侧的地质界面分开考虑.

1.1.2　地层界面的空间自相关和互相关

在地质领域中，地层作为一类地质实体，其空间分布也存在自相关和互相关的现象（如图2）. 正如美籍瑞士裔地理学家Waldo R. Tobler 于1970年提出的地理学第一定律（Tobler， 1970）所说，“Everything is related to everything else， but near things are more related to each other”，即地理事物或属性在空间分布上互为相关；距离越近，相关性越大；距离越远，相异性越大.

针对整合的沉积地层而言，地层的自相关和互相关特性主要源于地层的沉积过程. 丹麦地质学家Nicolaus Steno于1669年在《关于固体自然包裹于另一固体问题的初步探讨》一书中提出著名的地层层序律（Winter， 1968）：a.层序叠加律，在正常的岩层层序中，先形成的地层在下，后形成的地层在上；b.原始连续律，未经变动的地层在横向上连续延伸并逐渐尖灭；c.原始水平律，地层未经变动时时呈水平状.

假设 \left(x_j,y_j\right)位置上的地层面 i的高程函数为 Z_i \left(x_j,y_j\right)，根据叠加原理， Z_i \left(x_j,y_j\right)具有时空意义：相同的地层面形成于同一个时期，在不同位置 \left(x_j,y_j\right)间存在空间自相关现象；受沉积地层的形成机制影响，相邻的地层界面形成于相近的时期，在同一位置 \left(x_j,y_j\right)附近，不同时期的地层面之间形态存在相似性. 基于此原理，在钻孔数据较少的情况下，可以借助叠加原理通过这种新老地层间的相似性添加地质知识（Calcagno et al.， 2008； Zhu et al.， 2012； Perrin and Rainaud， 2013； De la Varga et al.， 2018），提高地质制图的准确性.

在不考虑断层的情况下，每一个地层面是一个连续的空间曲面. 由地层层序律可知，地层的空间自相关主要源于地层的连续性. 沉积地层在空间分布上具有连续性，服从原始连续律和原始水平律，但是地质过程中的随机因素导致局部变异，地层的自相关在结构性的基础上表现出一定的随机性. 这些空间曲面看似相互独立，但是在纵向上具有几何相似性，这种相似性正是源于地层的互相关. 地层间的互相关是新老地层在形成过程中的叠加规则的体现. 对于地层的自相关，可以通过地质统计学方法进行统计，并在平稳假设、内蕴假设等条件下对地层结构进行基于自相关规律的预测计算. 传统地质统计学方法多利用变量的空间自相关进行插值，对地层间的互相关考虑较少，但是相邻地层之间的信息可以通过互相关进行传递，互相关对地质变量预测和不确定性分析的影响不容忽视.

1.2　基于Copula的多维相关建模理论

为了同时利用两个地层的数据，需要一种能够对多个地质变量进行联合建模的工具，变量间的相关性建模是解决问题的关键. 很多空间统计学方法（如克里金、序贯高斯模拟等）基于多元高斯分布的假设构建多变量的联合分布，在错综复杂的地质作用下，很难保证地质测量数据都符合高斯分布. 有些方法（如对数正态克里金、指示克里金）虽然可以通过将非高斯分布的变量边际分布变换到高斯分布来处理，但是仍会面临相关结构不符合高斯假设的问题（Bárdossy and Li， 2008）.

Sklar定理（Sklar， 1959）指出，任何多元联合分布都可以用一元边际分布函数和一个描述变量之间相关关系的Copula函数来表示. Copula函数是利用变量间变化的一致性，将一维边际分布连接起来表达多维分布的连接函数，可以理解为一些边际是 [ \mathrm{0,1} ]上的均匀分布的随机变量构成的联合累积分布函数. 根据Sklar定理，令联合分布 F \left(z_{\mathrm{1}},z_{\mathrm{2}},...,z_n\right)的边际分布为 F_{\mathrm{1}} \left(z_{\mathrm{1}}\right) , F_{\mathrm{2}} \left(z_{\mathrm{2}}\right) , . . . , F_n \left(z_n\right)，则存在一个Copula函数 C \left(u_{\mathrm{1}},u_{\mathrm{2}},...,u_n\right)，使得

\begin{array}{l}F\left(z_{\mathrm{1}},z_{\mathrm{2}},...,z_n\right)=\\ C\left(F_{\mathrm{1}}\left(z_{\mathrm{1}}\right),F_{\mathrm{2}}\left(z_{\mathrm{2}}\right),...,F_n\left(z_n\right)\right)\end{array},(1)

其中： u_i = F_i \left(z_i\right).可以证明（Nelsen， 2007），对于连续的边际分布 F_{\mathrm{1}} \left(z_{\mathrm{1}}\right) , F_{\mathrm{2}} \left(z_{\mathrm{2}}\right) , . . . , F_n \left(z_n\right)，Copula函数 C \left(u_{\mathrm{1}},u_{\mathrm{2}},...,u_n\right)是唯一存在的. 根据Sklar定理可知，任意一个联合分布函数都可以表示其边际分布函数与连接函数的乘积，即：

\begin{array}{l}f\left(z_{\mathrm{1}},...,z_n\right)=\frac{{\partial }^{n}F\left(z_{\mathrm{1}},...,z_n\right)}{\partial z_{\mathrm{1}},...,\partial z_n}=\\ \frac{{\partial }^{n}C\left(F_{\mathrm{1}}\left(z_{\mathrm{1}}\right),...,F_n\left(z_n\right)\right)}{\partial z_{\mathrm{1}},...,\partial z_n}=\\ \frac{{\partial }^{n}C\left(u_{\mathrm{1}},...,u_n\right)}{\partial z_{\mathrm{1}},...,\partial z_n}\times {\prod }_{i=\mathrm{1}}^n\frac{\partial u_i}{\partial z_i}=\\ c\left(u_{\mathrm{1}},...,u_n\right)\times {\prod }_{i=\mathrm{1}}^n{f}_i\left(z_i\right)\end{array}，(2)

其中： f \left(z_{\mathrm{1}},...,z_n\right)是联合分布密度函数， f_i \left(z_i\right)是边际分布密度函数， c \left(u_{\mathrm{1}},...,u_n\right)是Copula密度函数. 同理，若已知 n维随机变量 \left(z_{\mathrm{1}},z_{\mathrm{2}},...,z_n\right)的边际分布 F_{\mathrm{1}} \left(z_{\mathrm{1}}\right) , F_{\mathrm{2}} \left(z_{\mathrm{2}}\right) , . . . , F_n \left(z_n\right)，则可以通过选取适当的Copula函数 C \left(u_{\mathrm{1}},u_{\mathrm{2}},...,u_n\right)得到它们的联合分布函数 F \left(z_{\mathrm{1}},z_{\mathrm{2}},...,z_n\right).

Copula的特殊性质提供了一种“尺度不变”的研究相关性度量的方法：随机变量在严格的单调变换下，其Copula函数是不变的，而边际可以任意改变. 即Copula函数取决于相关结构而不是边际分布，利用Copula表达的联合分布的特性在单调变换下是不变的（Wolff， 1981）. 因此，常用的数据转换（如自然对数或正态分数转换）对Copula函数没有影响（Bárdossy， 2006）.

为了模拟变量间各种特性的相关结构，不同的Copula函数被发展出来. 目前研究中采用较多的是椭圆族Copula（如Gaussian Copula，t Copula），阿基米德族Copula（如Frank Copula，Clayton Copula，Gumbel Copula）及其生存函数，此外还有很多其他类型的Copula函数可供选择（Sadegh et al.， 2017）. 各种Copula函数各自具有不同的特性，在实际应用中，可以根据不同的边际分布选择相应的Copula函数灵活构建合适的多元联合分布模型.

由于地质变量间的相关性可能是非线性的，基于Pearson线性相关系数和协方差的空间统计建模方法对非线性相关的问题难以给出令人满意的答案. Kendall \tau秩相关系数通过对变量观测值排序，利用变量观测值的秩计算变量间的相关系数. 因为变量的秩只与数据的排列顺序有关，与实际观测值无关，所在非线性相关和线性单调变换的情况下仍能有效地度量相关性. 由于Copula与Kendall \tau秩相关系数都具有在单调变换时保持不变的非线性优势，人们往往采用kendall \tau秩相关系数计算Copula函数，而不采用受到线性相关限制的Pearson线性相关系数.

Kendall \tau秩相关系数是衡量变量间变化一致性程度的指标，令 \left(a_{\mathrm{1}},b_{\mathrm{1}}\right)和 \left(a_{\mathrm{2}},b_{\mathrm{2}}\right)是随机变量 \left(a,b\right)的两组观测值，如果 且 ，或 a_{\mathrm{1}} > a_{\mathrm{2}}且 b_{\mathrm{1}} > b_{\mathrm{2}}，即 \left(a_{\mathrm{1}}-a_{\mathrm{2}}\right) \cdot \left(b_{\mathrm{1}}-b_{\mathrm{2}}\right) > \mathrm{0}，则称 \left(a_{\mathrm{1}},b_{\mathrm{1}}\right)和 \left(a_{\mathrm{2}},b_{\mathrm{2}}\right)是一致的；若 且 b_{\mathrm{1}} > b_{\mathrm{2}}，或 a_{\mathrm{1}} > a_{\mathrm{2}}且 ，即 ，则称 \left(a_{\mathrm{1}},b_{\mathrm{1}}\right)和 \left(a_{\mathrm{2}},b_{\mathrm{2}}\right)是不一致的. Kendall \tau定义为变量 a和 b变化一致性的概率和不一致性的概率之差：

.(3)

Kendall \tau秩相关系数与二维Copula函数 C \left(u,v|\theta \right)存在以下关系（Nelsen， 2000）：

\tau = \mathrm{4} {\int }_{\mathrm{0}}^{\mathrm{1}} {\int }_{\mathrm{0}}^{\mathrm{1}}C\left(u,v|\theta \right) \mathrm{d} C \left(u,v|\theta \right) - \mathrm{1}.(4)

当变量 u和 v的Kendall \tau秩相关系数已知，对应Copula函数的参数 \theta可以通过上式求解. 一些Copula族的参数甚至与Kendall \tau之间存在一对一的对应关系. 采用Kendall \tau秩相关系数估计的Copula函数参数 \theta与变量 u , v的边际分布函数无关.

1.3　贝叶斯框架下的地层面模型修正

针对深部地层底板位置信息不准确造成的模型误差问题，本节提出一种在贝叶斯框架下基于地层互相关的模型修正方法，无需新增数据和人工交互编辑，充分利用相邻地层数据隐含的互信息，借助相邻地层的相关结构为地质模型提供互相关约束，对深部地层面模型进行一定程度上的修正.

1.3.1　地层面模型先验分布

贝叶斯推断可以通过补充的信息对先验认识进行概率更新. 对于已建好的模型，贝叶斯方法可以通过似然函数引入条件约束，基于先验模型更新模型的参数. 相比于利用插值方法重新建模，在贝叶斯更新中对约束选择更加灵活.

在贝叶斯框架下，已建好的地层面模型被视为先验模型. 令模型中 \left(x,y\right)位置上待更新界面的高程值 Z作为模型参数，通过构建模型参数在给定辅助数据（相邻地层面的高程 Z_{\mathrm{a}\mathrm{d}\mathrm{j}}）时的似然函数 L \left(Z|Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right)，对模型位置参数 Z的先验分布 p \left(Z\right)利用贝叶斯推断进行更新，得到模型参数的后验分布 p \left(Z|Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right)：

p \left(Z|Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right) \propto p \left(Z\right) \cdot L \left(Z|Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right).(5)

似然函数 L \left(Z|Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right)可以基于Copula方法构建的联合分布 f \left(Z,Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right)和相邻地层的边际分布 f \left(Z_{\mathrm{a}\mathrm{d}\mathrm{j}}\right)得到. 式中， \left(x,y\right)位置上的地层面模型的先验分布 p \left(Z\right)按照 \left(x,y\right)位置处钻孔是否包含待修正地层的不完整观测数据分为两种情况（如图3）.

在图3中，将地层A的底面简记作地层面A，将地层B的底面简记作地层面B. 以待修正地层面B和相邻地层面A为例，在2号钻孔处，假设由于钻孔2未穿过地层面B，地层面B的样本数据未能准确获取，地层面B的实际位置 {\mathrm{B}}_{\mathrm{2}}未知，此时有：情况a. 2号钻孔的数据未记录地层面B的任何位置信息；情况b. 2号钻孔的数据记录的地层B的最深位置为 {{\mathrm{B}}_{\mathrm{2}}}^{\mathrm{\text{'}}}，但是实际地层底面位置 {\mathrm{B}}_{\mathrm{2}}仍未知.

对于情况a，2号钻孔位置 \left(x,y\right)处，地层面B的先验分布 p \left(Z_{\mathrm{B}}\right)可以采用Tacher等人提出的地层面模型不确定性评价方法（Tacher et al.， 2006）进行估计，即：将已有地层面模型在 \left(x,y\right)位置上的高程值 z_{\mathrm{m}} \left(x,y\right)作为模型先验分布的期望，利用样本数据计算 \left(x,y\right)位置上的克里金估值方差 {\sigma }_{\mathrm{k}}^{\mathrm{2}} \left(x,y\right)作为分布的方差，构建一个正态分布 N ( z_{\mathrm{m}} , {\sigma }_{\mathrm{k}}^{\mathrm{2}} )作为 \left(x,y\right)位置上待修正地层面模型的先验分布.

对于情况b，2号钻孔位置处地层面B的真实高程值 Z_{\mathrm{B}\mathrm{2}}应小于2号钻孔数据记录的地层B的最深位置 {{\mathrm{B}}_{\mathrm{2}}}^{\mathrm{\text{'}}}处的高程值 {Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}，此时有不等式约束 Z_{\mathrm{B}\mathrm{2}} \le {Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}. 2号钻孔位置处 Z_{\mathrm{B}\mathrm{2}}的概率分布 f \left(Z_{\mathrm{B}\mathrm{2}}\right)可以用Tacher等人的方法根据邻域内周围钻孔中地层面B的样本点计算，同时考虑不等式约束 Z_{\mathrm{B}\mathrm{2}} \le {Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}，将分布 f \left(Z_{\mathrm{B}\mathrm{2}}\right)转化为 Z_{\mathrm{B}\mathrm{2}} \le {Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}的截尾分布（如图4）. Z_{\mathrm{B}\mathrm{2}}的截尾分布概率密度函数 f \left(Z_{\mathrm{B}\mathrm{2}};{Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right)可以用下面公式计算：

\left\{\begin{array}{c}f\left(Z_{\mathrm{B}\mathrm{2}};{Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right)=\frac{f\left(Z_{\mathrm{B}\mathrm{2}}\right)}{F\left({Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right)},Z_{\mathrm{B}\mathrm{2}}\in \left.\left(-\mathrm{\infty },{Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right.\right]\\ f\left(Z_{\mathrm{B}\mathrm{2}};{Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right)=\mathrm{0},Z_{\mathrm{B}\mathrm{2}}\in \left({Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}},+\mathrm{\infty }\right)\end{array},\right. (5,6)

式中： F \left({Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right)是对应 Z_{\mathrm{B}\mathrm{2}} \le {Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}的累积概率密度函数. 在2号钻孔位置处，截尾分布 f \left(Z_{\mathrm{B}\mathrm{2}};{Z_{\mathrm{B}\mathrm{2}}}^{\mathrm{\text{'}}}\right)即为地层面B的先验分布 p \left(Z_{\mathrm{B}}\right).

1.3.2　基于地层相关结构约束的模型修正

当待修正的地层面上侧有一个位置准确的相邻地层面能够作为条件约束时，可以利用这两个地层面构建双地层面高程函数联合分布，采用双变量Copula函数对地层面的相关结构建模，然后根据联合分布计算似然函数，对先验模型进行贝叶斯更新修正. 下面介绍针对两个相邻地层组成的场景，如何利用双变量Copula和贝叶斯推断方法实现对先验模型的修正.

假设AB两个地层面为沉积地层中两个相邻的地层面，A和B的高程变量间存在稳定的互相关结构，AB两个地层面的模型是仅使用各自地层的采样数据独立建模得到的，B地层为待修正地层. 在空间位置 s_{\mathrm{0}} = \left(x,y\right)处，模型中地层面B的高程概率密度函数为 f_{\mathrm{B}} \left(Z_{\mathrm{B}}\right)，利用双变量条件Copula和贝叶斯更新方法通过A的数据 V_{\mathrm{A}}修正B的先验模型值 Z_{\mathrm{B}\mathrm{U}}的主要步骤为：

（1）根据模型不确定性评价方法估计B地层面的模型值 Z_{\mathrm{B}\mathrm{U}}的概率分布 f \left(Z_{\mathrm{B}}\right)作为地层面模型参数 Z_{\mathrm{B}}的先验概率分布，根据不完整的测量数据t将先验的正态分布 f_{\mathrm{B}} \left(Z_{\mathrm{B}}\right)转化为截尾分布 f \left(Z_{\mathrm{B}};t\right)（图4）.

（2）统计每个地层面的样本数据，估计AB两个地层面高程变量的边际分布 F_{\mathrm{A}} , F_{\mathrm{B}}，将AB界面高程边际分布变换为 \left[\mathrm{0,1}\right]均匀分布下的分位数 u_{\mathrm{A}} , u_{\mathrm{B}}，估计双变量Copula函数 C \left(u_{\mathrm{A}},u_{\mathrm{B}}\right)的最佳函数族和对应参数；

（3）根据双变量Copula密度函数 c_{\mathrm{A},\mathrm{B}} \left(u_{\mathrm{A}},u_{\mathrm{B}}\right)得到联合概率密度函数 f_{\mathrm{A},\mathrm{B}} \left(Z_{\mathrm{A}},Z_{\mathrm{B}}\right)，由A地层面的高程值概率密度函数 f_{\mathrm{A}} \left(Z_{\mathrm{A}}\right)和 s_{\mathrm{0}} = \left(x,y\right)处的高程值 V_{\mathrm{A}}，计算在 V_{\mathrm{A}}约束下的B地层面的似然函数：

L \left(Z_{\mathrm{B}}|Z_{\mathrm{A}}=V_{\mathrm{A}}\right) = f_{\mathrm{A}|\mathrm{B}} \left(Z_{\mathrm{A}}=V_{\mathrm{A}}|Z_{\mathrm{B}}\right).(7)

（4）对先验截尾分布 f \left(Z_{\mathrm{B}};t\right)进行贝叶斯更新：

f \left(Z_{\mathrm{B}}|Z_{\mathrm{A}}=V_{\mathrm{A}}\right) \propto f \left(Z_{\mathrm{B}};t\right) \cdot L \left(Z_{\mathrm{B}}|Z_{\mathrm{A}}=V_{\mathrm{A}}\right).(8)

归一化后即可得到B地层面的后验概率分布 f_{\mathrm{B}|\mathrm{A}} \left(Z_{\mathrm{B}}|Z_{\mathrm{A}}=V_{\mathrm{A}}\right)，计算条件期望值或中位数作为模型修正值. 其中，条件期望的计算公式为：

E_{\mathrm{B}|\mathrm{A}} = \underset{R}{\int } f\left(Z_{\mathrm{B}}|Z_{\mathrm{A}}=V_{\mathrm{A}}\right)\mathrm{d}{Z}_{\mathrm{B}}.(9)

重复上述步骤，直到所有需要被修正的位置 s_{\mathrm{0}}被遍历，完成界面模型的修正. 基于双地层Copula的模型修正方法流程如图5所示.

需要注意的是，根据样本数据计算相关结构模型的参数（Copula族和Copula参数）时，只考虑待修正地层和相邻地层数据同时存在的区域. 对待修正地层模型进行修正时也只涉及该区域. 此外，模型修正处理过程中也存在不确定性，修正后的模型值不仅受到先验模型准确性的影响，而且在修正过程中需要进行基于Copula函数的相关结构建模，该相关结构模型的准确性也会影响最终的模型修正效果. 理论上，作为条件约束的辅助地层面的高程数据准确性需要高于被修正地层面的高程数据准确性. 当辅助地层条件约束直接采用测量得到的硬数据时，修正效果最佳，当辅助地层的模型准确性明显高于待修正地层时，也可以使用准确性较高的辅助地层模型估计值去修正准确性较低的待修正地层模型估计值. 当辅助地层的采样或建模准确性较低时，不再适合作为待修正地层的约束条件.

2 实验

2.1　研究区概况

本节以北海海岸带地区地质数据为例，对基于相邻地层互相关约束的模型修正方法进行实验验证. 研究区位于广西北海市沿海地区，南及南西面临海，北面靠陆地，地势北高南低及北东高南西低. 大部分为平原区，海拔一般10~30 m，在海滩地区，海拔一般为1~5 m. 北东及南西侧冠头领有少量垄状低丘分布，一般海拔为45~90 m，最高是位于东面的采花岭117.80 m，垄状低丘呈南西-北东向起伏延伸. 工作区内广泛分布第四系，局部出露有第三系、泥盆系含石英斑岩脉. 第四系、第三系以松散沉积土体为主，岩性为砂、砾砂、含黏土砂砾、粉质黏土、砂质黏土、黏土等，泥盆系岩性主要是粉砂岩、泥质粉砂岩与页岩互层、细砂岩，局部含石英斑岩脉. 研究区钻孔的钻探结果显示该区域勘探地层主要有：第四系人工填筑层（Q^ml）、第四系冲积层（Q^al）、第四系海积层（Q^m）、第四系海陆混合堆积层（Q^mc）、第四系北海组（Q_p ^b）、第四系湛江组（Q_p ^z）、第三系南康群上段（N_n）. 土层按其性质和状态自上而下分为8层，分别为填土（①）、耕（表）土（②）、中粗砂（③）、於泥质黏土（④）、粉质黏土（⑤）、黏土（⑥）、砾砂（⑦）、粉细砂（⑧）等. 其中，⑤和⑥又细分为：⑤-1可塑状粉质黏土、⑤-2砾砂和⑥-1可塑状粘土、⑥-2粉细砂. 对研究区建模得到的三维地质模型如图6所示.

2.2　地层面模型修正

研究区共89个钻孔，其中有80个钻孔的样本数据中均含有⑤-2（砂砾）和⑥-1（可塑状粘土）地层面的位置信息. 在这80个钻孔中，有7个钻孔存在⑥-1地层未钻透，地层底板记录不准确的情况，剩余73个钻孔的⑥-1地层面位置信息是准确的. 下面实验利用⑤-2地层面样本数据修正⑥-1地层面模型. 将位于上侧的⑤-2地层下底面作为辅助地层面，位于下侧的⑥-1地层下底面作为待修正地层面，取⑥-1地层下底面位置准确的样本（73个）作为验证集，利用留一法交叉验证检验修正前后的模型精度. 为了展示所提方法对1.3.1节中a和b两种情况的差别，本实验对验证集中的样本添加了0.5倍地层厚度作为系统误差，即在真实界面位置的基础上升高当前位置地层厚度值的50%作为假设的钻孔最大深度，在实验中通过忽视不准确测量结果直接计算先验分布模拟情况a，通过构建截尾分布模拟情况b.

实验对验证集中每一个验证位置，以除该位置之外的验证集样本作为已知数据，利用反距离加权方法（inverse distance weighted， IDW）对验证位置上待修正地层面的高程进行预测，利用克里金预测方差作为先验模型不确定性的估计，计算预测结果与验证集的真实采样数据之间的偏差，统计得到每个地层面模型的绝对值平均误差（mean absolute error， MAE）、平均误差（mean error， ME）和均方根误差（root mean square error， RMSE）作为模型准确性指标. 平均误差（ME）能反映预测结果的系统偏差，平均绝对误差（MAE）可以反映预测精度，均方根误差（RMSE）通过标准差的形式反映了预测结果中的极值情况（Daly， 2006），总体上，误差越小，预测精确度越高.

本实验中IDW和克里金的搜索邻域均设置为最近的3个点，IDW权重采用距离平方的倒数. 实验将验证位置上的IDW插值结果作为先验模型，利用IDW预测值和该位置上的克里金预测方差构建模型先验分布，并根据不准确的钻孔位置信息计算截尾分布，借助相邻层的样本数据，通过Copula方法建立相邻地层的相关结构模型，基于贝叶斯推断方法对带有误差的模型值进行更新修正. 对修正后的模型值再次与真实采样数据比较，计算MAE、ME和RMSE这3个误差指标，进行修正前后的误差对比. 为了说明在模型修正中，利用不准确测量数据做约束（情况b）相比于忽略该数据（情况a）的效果差别，实验分别计算了两种情况下的修正后模型误差. 为了便于标注，在下文图中分别以 Z_{\mathrm{1}} , Z_{\mathrm{2}}指代⑤-2和⑥-1地层的下底面高程变量，以 u_{\mathrm{1}} , u_{\mathrm{2}}指代⑤-2和⑥-1地层的下底面高程分位数. ⑤-2和⑥-1地层面高程值的数值分布图如图7所示. 根据钻孔数据计算的⑤-2和⑥-1地层面高程值统计直方图如图8所示，图8a为⑤-2地层面高程统计直方图（红线为根据样本数据拟合的概率密度函数曲线）；图8b为⑤-2地层面高程样本数据拟合累积概率密度函数；图8c为⑥-1地层面高程统计直方图（红线为根据样本数据拟合的概率密度函数曲线）；图8d为⑥-1地层面高程样本数据拟合累积概率密度函数.

对⑤-2和⑥-1地层的准确测量样本数据进行统计，计算得到两个地层底面的秩相关系数kendall \tau = \mathrm{0.56}. 最佳拟合Copula函数为Survival Gumbel Copula （ \theta = \mathrm{2.72}）：

\begin{array}{l}{C}_{\mathrm{S}\mathrm{G}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{l}}\left(u,v\right)=\\ \mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\left[{\left(-\mathrm{l}\mathrm{n}\left(\mathrm{1}-u\right)\right)}^{\theta }+{\left(-\mathrm{l}\mathrm{n}\left(\mathrm{1}-v\right)\right)}^{\theta }\right]}^{\frac{\mathrm{1}}{\theta }}\right\},\theta \in \left.\left[\mathrm{1},+\mathrm{\infty }\right.\right)\end{array},(10)

⑤-2和⑥-1地层面高程Copula如图9所示，图9a为⑤-2地层面（u₁）与⑥-1地层面（u₂）Copula密度立体图；图9b为⑤-2地层面（u₁）与⑥-1地层面（u₂）Copula密度等值线图，图中变量边际分布均为 \left[\mathrm{0,1}\right]均匀分布. 利用该Copula函数计算⑤-2和⑥-1地层面联合分布和⑥-1地层面高程似然函数，利用似然函数更新⑥-1地层面高程先验分布，得到后验分布.

图10以验证位置#8钻孔为例，展示通过⑤-2地层面样本数据修正⑥-1地层面时，⑥-1地层面模型概率分布的更新过程. #8钻孔位置上⑤-2地层面样本高程测量值 V_{\mathrm{1}} = - \mathrm{24.60} m，⑥-1地层面真实样本高程测量值 V_{\mathrm{2}} = - \mathrm{34.60} m，模拟的不准确高程测量值 {V_{\mathrm{2}}}^{\mathrm{\text{'}}} = - \mathrm{29.60} m，利用周围钻孔的IDW插值得到的⑥-1地层面模型高程值 V_{\mathrm{m}} = - \mathrm{15.08} m，克里金估计方差 {\sigma }_{\mathrm{k}}^{\mathrm{2}} = \mathrm{22.77}. 不考虑不准确的测量数据时，#8钻孔位置上⑥-1地层面模型先验分布 f \left(z_{\mathrm{2}}\right) \sim N \left(-\mathrm{15.08,22.77}\right)，如图10中绿线所示. #8钻孔位置上⑥-1地层面似然函数 f \left(z_{\mathrm{1}}=V_{\mathrm{1}}|z_{\mathrm{2}}\right)如图11所示. 不考虑#8钻孔位置上不准确的测量数据 {V_{\mathrm{2}}}^{\mathrm{\text{'}}}时，#8钻孔位置上⑥-1地层面后验分布 f \left(z_{\mathrm{2}}|z_{\mathrm{1}}=V_{\mathrm{1}}\right)如图10中红色曲线所示，后验期望值 E \left(z_{\mathrm{2}}\right) = - \mathrm{22.61} m. 考虑不准确测量数据 {V_{\mathrm{2}}}^{\mathrm{\text{'}}}的不等式约束， #8钻孔位置上⑥-1地层面模型的截尾分布 f \left(z_{\mathrm{2}}|z_{\mathrm{1}}=V_{\mathrm{1}};{V_{\mathrm{2}}}^{\mathrm{\text{'}}}\right)如图10中紫线所示，后验期望值 E {\left(z_{\mathrm{2}}\right)}_{\mathrm{T}} = - \mathrm{31.41} m. 可以看出，相比于利用周围钻孔插值得到的IDW先验模型，贝叶斯更新后的模型后验分布的期望更接近真实样本值；考虑不准确测量值 {V_{\mathrm{2}}}^{\mathrm{\text{'}}}的后验分布 f \left(z_{\mathrm{2}}|z_{\mathrm{1}}=V_{\mathrm{1}};{V_{\mathrm{2}}}^{\mathrm{\text{'}}}\right)比忽略 {V_{\mathrm{2}}}^{\mathrm{\text{'}}}的分布 f \left(z_{\mathrm{2}}|z_{\mathrm{1}}=V_{\mathrm{1}}\right)更接近真实情况（图11）.

类似地，对73个验证位置上的模型值进行修正，修正前后的模型值结果（部分）如表1所示，误差统计结果（MAE， ME， RMSE）如表2所示. 表中“未截尾修正值”表示修正过程中忽略不准确的样本的结果，“截尾修正值”表示修正过程中利用不准确的样本进行不等式约束的结果.

针对⑥-1地层未钻透的7个钻孔位置上地层底板记录不准确的情况，本文也进行了修正实验，修正实验结果如表3所示.

2.3　结果与讨论

从实验结果可以看出，虽然修正后的⑥-1地层面模型值与验证集中的钻孔真实样本值存在一定偏差，但在大多数位置，修正后的模型值比之前的IDW预测值更接近真实样本值. 总体上看，模型修正后，地层面模型的误差降低. 需要指出，利用Copula方法的模型修正可以在统计意义上改善模型，但是不能保证对每个位置的修正都会减小误差. 相对于未截尾的先验分布，从截尾后的先验分布得到的后验修正值绝对值平均误差（MAE）和均方根误差（RMSE）降低，平均误差（ME）显著向负数方向移动. 这是由于由不准确的样本信息提供了不等式约束，不等式所包含的信息降低了不确定性，提高了建模精度，平均误差（MAE）和均方根误差（RMSE）降低. 同时地层面的期望值应位于不准确样本的下方，不等式造成先验分布变成右截尾分布，先验分布的变化导致后验分布的期望也整体向负数方向移动，而平均误差是所有真误差的平均值，统计上看，整体真误差的均值也向负数方向移动，平均误差（ME）变为负值.

从信息论的角度来看，事件发生的不确定性来源于未获得的信息. 地层结构通过同一地层的空间自相关和相邻地层的互相关影响着信息的传递. 由于相关性的存在，从一个信源获得的信息包含了与其他信源的互信息. 互信息可以看作是一个随机变量由于知道另一个随机变量而减少的不确定性. 根据样本值的统计，在不考虑边际分布差异的情况下（均转换到 [ \mathrm{0,1} ]上的均匀分布），⑤-2地层样本数据信息量为1.384 bit，⑥-1地层样本数据的信息量为1.385 bit. 两地层数据的互信息为0.395 bit. 这说明相邻界面间存在这相关性，相邻地层的数据中隐藏着可用的信息，有效利用这部分互信息可以增进人们对地质规律的认识，更准确地预测地质现象. 基于相邻地层互相关约束的模型修正方法正是利用这种相邻地层面的互信息来降低模型的不确定性，实现模型的修正. 由于利用Copula函数建立的相关结构模型中同样存在不确定性，因此修正后的结果也难免与真实值存在误差.

3 结论

为了提高地质模型中地层底板位置的准确性，本文提出了一种基于地层相关结构约束的地质模型修正方法. 该方法首先基于Copula理论，建立相邻地层面的互相关结构模型，然后构造待修正地层面的似然函数，最后在贝叶斯框架中，利用不准确的钻孔测量数据和相邻地层的数据约束对先验界面模型进行更新修正. 本文利用北海海岸带地质数据对IDW方法得到的地层面插值模型进行了修正实验，通过交叉验证对比了修正前后的模型误差. 实验结果表明基于地层相关结构约束的模型修正方法无需新增其他类型辅助数据，方法可以利用相邻地层的数据约束降低模型误差，提高地层面模型的建模准确性.

参考文献

[1] Abrahamsen, P., Omre, H., 1994. Random Functions and Geological Subsurfaces. Ecmor IV-European Conference on the Mathematics of Oil Recovery, Røros, Norway, 21. https://doi.org/10.3997/2214-4609.201411142

[2] Autin, J., Scheck-Wenderoth, M., Götze, H. J., et al., 2016. Deep Structure of the Argentine Margin Inferred from 3D Gravity and Temperature Modelling, Colorado Basin. Tectonophysics, 676(47): 198-210. https://doi.org/10.1016/j.tecto.2015.11.023

[3] Bárdossy, A., 2006. Copula-Based Geostatistical Models for Groundwater Quality Parameters. Water Resources Research, 42(11). https://doi.org/10.1029/2005wr004754

[4] Bárdossy, A., Li, J., 2008. Geostatistical Interpolation Using Copulas. Water Resources Research, 44(7). https://doi.org/10.1029/2007wr006115

[5] Calcagno, P., Chiles, J. P., Courrioux, G., at al., 2008. Geological Modelling from Field Data and Geological Knowledge: Part I. Modelling Method Coupling 3D Potential-Field Interpolation and Geological Rules. Physics of the Earth and Planetary Interiors, 171 (1-4): 147-157. https://doi.org/10.1016/j.pepi.2008.06.013

[6] Caumon, G., Journel, A. G., 2004. Early Uncertainty Assessment: Application to A Hydrocarbon Reservoir Appraisal. Geostatistics Banff 2004, Springer, Dordrecht, 551-557. https://doi.org/10.1007/978-1-4020-3610-1_56

[7] Daly, C., 2006. Guidelines for Assessing the Suitability of Spatial Climate Data Sets. International Journal of Climatology, 26(6): 707-721. https://doi.org/10.1002/joc.1322

[8] de Kemp, E. A., Sprague, K. B., 2003. Interpretive Tools for 3-D Structural Geological Modeling Part I: B´ezier-Based Curves, Ribbons and Grip Frames. Geoinformatica, 7 (1): 55-71. https://doi.org/10.1023/a:1022822227691

[9] De la Varga, M., Schaaf, A., Wellmann J. F., 2018. GemPy 1.0: Open-Source Stochastic Geological Modeling and Inversion. Geoscientific Model Development Discussions, 1-50. https://doi.org/10.5194/gmd-12-1-2019

[10] De Paor, D. G., 1996. B´ezier Curves and Geological Design. In: De Paor, D. G., ed., Computer Methods in the Geosciences. Pergamon, 389-417. https://doi.org/10.1016/s1874-561x(96)80031-9

[11] Fremming, N. P., 2002. 3D Geological Model Construction Using a 3D Grid. ECMOR VIII-8th European Conference on the Mathematics of Oil Recovery. https://doi.org/10.3997/2214-4609.201405917

[12] Gradmann, S., Ebbing, J., Fullea, J., 2013. Integrated Geophysical Modelling of a Lateral Transition Zone in the Lithospheric Mantle under Norway and Sweden. Geophysical Journal International, 194 (3): 1358-1373. https://doi.org/10.1093/gji/ggt213

[13] Haase, C., Ebbing, J., Funck, T., 2017. A 3D Regional Crustal Model of the NE Atlantic Based on Seismic and Gravity Data. Geological Society, London, Special Publications, 447 (1): 233-247. https://doi.org/10.1144/sp447.8

[14] Hou, Z., Wang, T. Y., Yu, C. C., et al., 2018. Study of 3d Geological Modeling Based on Aeromagnetic Data. Advances in Earth Science, 33(3): 257-269 (in Chinese with English abstract).

[15] Kessler, H., Mathers, S., Sobisch, H. G., 2009. The Capture and Dissemination of Integrated 3D Geospatial Knowledge at the British Geological Survey Using GSI3D Software and Methodology. Computers & Geosciences, 35 (6): 1311-1321. https://doi.org/10.1016/j.cageo.2008.04.005

[16] Lemon, A. M., Jones, N. L., 2003. Building Solid Models from Boreholes and User-Defined Cross-Sections. Computers & Geosciences, 29(5): 547-555. https://doi.org/10.1016/s0098-3004(03)00051-7

[17] Liu, Z. F., Wei, Z. H., Huang, X. J., et al., 2012. Research on Dynamic Update Method of 3d Geological Model: A Case Study of Water Resources and Hydropower Projects. China Rural Water and Hydropower, 11:93-96 (in Chinese with English abstract).

[18] Lyu, Y. H., Du, Y., Zou, C. Y., 2007. Improved Numerical Simulation Method of Complex Fault Block Oil Reservoir. Fault-Block Oil & Gas Field, 14(6): 21-22, 94 (in Chinese with English abstract).

[19] Ming, J., Pan, M., Qu, H. G., et al., 2010. GSIS: A 3D Geological Multi-Body Modeling System from Netty Cross-Sections with Topology. Computers & Geosciences, 36(6): 756-767. https://doi.org/10.1016/j.cageo.2009.11.003

[20] Nelsen, R. B., 2007. An Introduction to Copulas. Springer, New York.

[21] Perrin, M., Rainaud, J. F., 2013. Shared Earth Modeling: Knowledge Driven Solutions for Building and Managing Subsurface 3D Geological Models, Editions Technip.

[22] Sadegh, M., Ragno, E., Aghakouchak, A., 2017. Multivariate Copula Analysis Toolbox (MvCAT): Describing Dependence and Underlying Uncertainty Using a Bayesian Framework. Water Resources Research, 53. https://doi.org/10.1002/2016WR020242.

[23] Sklar, A, 1959. Fonctions de Repartition An Dimensions et Leurs Marges. Publications de l'Insitut de Statistique de Paris, 229-231

[24] Sprague, K. B., de Kemp, E. A., 2005. Interpretive Tools for 3-D Structural Geological Modelling Part II: Surface Design from Sparse Spatial Data. GeoInformatica, 9(1): 5-32. https://doi.org/10.1007/s10707-004-5620-8

[25] Tacher, L., Pomian-Srzednicki, I., Parriaux, A., 2006. Geological Uncertainties Associated with 3-D Subsurface Models. Computers & Geosciences, 32(2): 212-221. https://doi.org/10.1016/j.cageo.2005.06.010

[26] Tan, F., Wang, J., Jiao, Y.Y., et al., 2021. Current Situation and Development of Urban Underground Space Suitability Evaluation. Earth Science, 46(5): 1896-1908 (in Chinese with English abstract).

[27] Tearpock, D. J., Bischke, R. E., 1991. Applied Subsurface Geological Mapping. Prentice Hall.

[28] Tobler, W. R., 1970. A Computer Movie Simulating Urban Growth in the Detroit Region. Economic Geography, 46: 234. https://doi.org/10.2307/143141

[29] Wellmann, J. F., de la Varga, M., Murdie, R. E., et al., 2017. Uncertainty Estimation for a Geological Model of the Sandstone Greenstone Belt, Western Australia: Insights from Integrated Geological and Geophysical Inversion in a Bayesian Inference Framework. Geological Society, London, Special Publications, 453(1): 41-56. https://doi.org/10.1144/sp453.12

[30] Wellmann, J. F., Caumon, G., 2018. 3-D Structural Geological Models: Concepts, Methods, and Uncertainties, Advances in Geophysics, 59: 1-121. https://doi.org/10.1016/bs.agph.2018.09.001

[31] Winter, J. G., 1968.The Prodromus of Nicolaus Steno's Dissertation Concerning A Solid Body Enclosed by Process of Nature within A Solid, Hafner.

[32] Wolff, B, S. F., 1981. On Nonparametric Measures of Dependence for Random Variables. Annals of Statistics, 9(4):879-885. https://doi.org/10.1214/aos/1176345528

[33] Xie, J. R., Qiao, S. F., Qian, H., et al., 2014. The Application of Virtual Drilling Technology in the Three-Dimensional Geological Modeling of Water Resources and Hydropower. Journal of Railway Science and Engineering, (3):123-128 (in Chinese with English abstract).

[34] Xiu, C. H., Che, D. F., Jia, G. B., 2015. A 3D Dynamic Modeling Method for Coal Seams with Complex Geological Structures. Mine Surveying, 6: 52-55 (in Chinese with English abstract).

[35] Zhang, X.L., Wu, C.L., Zhou, Q., et al., 2020. Multi-Scale 3D Modeling and Visualization of Super Large Manganese Ore Gathering Area in Guizhou China. Earth Science, 45(2): 634-644 (in Chinese with English abstract).

[36] Zhu, L. F., Wu, X. C., Pan, X., 2006. Mechan6ism and Implementation of Error Correction for 3D Strata Model. Rock and Soil Mechanics, 27(2):268-271 (in Chinese with English abstract).

[37] Zhu, L. F., Wu, X. C., Pan, X., 2009. Theory of Accuracy Assessment and Methods for Error Correction in 3D Geological Structure Models. Earth Science Frontiers, 16(4): 363-371 (in Chinese with English abstract).

[38] Zhu, L. F., Zhang, C. J., Li, M. J., et al., 2012. Building 3D Solid Models of Sedimentary Stratigraphic Systems from Borehole Data: An Automatic Method and Case Studies. Engineering Geology, 127(8): 1-13. https://doi.org/10.1016/j.enggeo.2011.12.001

[39] 侯征, 王天意, 于长春, 等, 2018. 基于航磁数据的三维地质建模研究. 地球科学进展, 33(3): 257-269.

[40] 刘志锋, 魏振华, 黄笑鹃, 等, 2012. 三维地质模型动态更新方法研究——以水利水电工程为例. 中国农村水利水电, 11:93-96.

[41] 吕迎红, 杜燕, 邹存友, 等, 2007. 复杂断块油藏地质模型修正技术探讨. 断块油气田, 14(6): 21-22, 94.

[42] 谭飞, 汪君, 焦玉勇, 等, 2021. 城市地下空间适宜性评价研究国内外现状及趋势. 地球科学, 46(5): 1896-1908.

[43] 谢济仁, 乔世范, 钱骅, 等, 2014. 虚拟钻孔技术在水利水电三维地质建模中的应用. 铁道科学与工程学报, 3: 123-128.

[44] 修春华, 车德福, 贾国兵, 2015. 含复杂地质构造的三维煤层动态建模方法. 矿山测量, 6:52-55.

[45] 张夏林, 吴冲龙, 周琦, 等, 2020. 贵州超大型锰矿集区的多尺度三维地质建模. 地球科学, 45(2): 634-644.

[46] 朱良峰, 吴信才, 潘信, 2006. 三维地层模型误差修正机制及其实现技术.岩土力学, 27(2):268-271.

[47] 朱良峰, 吴信才, 潘信, 2009. 三维地质结构模型精度评估理论与误差修正方法研究. 地学前缘, 16(4):363-371.

基金资助

国家重点研发计划课题“城市地下全要素信息集成与智能建模技术”(2019YFC0605102)