基于构造恢复理论的含复杂断层三维地质建模方法

3D Geological Modeling Method Based on Tectonic Restoration Theory

摘要

在三维地质建模中，许多工作集中于断层本身的模拟，而忽略了断层时序性对地层构造的影响，然而考虑断层构造时序性及其导致的变形是十分必要和关键的.针对此问题，基于断裂恢复与演化影响域理论，提出了断裂矢量场建模方法用于含复杂断裂网络的三维地质建模，先使用断裂矢量场位移算子以断裂构造演化逆序恢复影响域范围内的地层和断层数据，再以演化正序逐步计算受断层影响后的地层和断层数据，获得含复杂断层网络的三维地质模型.通过建模实验和对比实验，验证了方法对数据具有更高的利用率以及在处理断裂接触时的能力，并且在建模过程会自动计算断裂导致的一切位移，提高了建模合理性和效率.方法充分考虑了断层构造时序性的影响，更适用于解决具有复杂接触关系的断裂网络模型构建问题.

关键词

三维建模 / 复杂断层网络 / 断裂构造恢复及演化 / 构造地质学

Key words

3D modeling / complex fault network / restoration and evolution of fault structure / structural geology

0 引言

研究断层及其三维结构建模技术具有重要的理论意义和实际意义（邰文星等，2023； 花卫华等，2023），目前构建三维地质模型技术较为成熟，但由于断层结构复杂、数据不连续等特征，含断层的建模方法一直是三维地质建模中的一个难点（李昌领等，2013； 王诏等，2015； 张建勇等，2021），国内外相关研究仍处于探索期.目前常用断层三维地质建模方法主要有4类：整体法、局部法、断层与地层的统一建模法和地层变形位移法.整体法基于地层恢复的思想，先恢复地层，再将断层两侧的地层看成一个整体，根据断层数据对三维地质模型进行恢复处理.沙浩（2016）依据整体法使用工程地质数据进行含断层的三维地质模型构建.该方法建模概念简单，但不便于进行局部修改，适用于构建断层数目较少、断距不大的地质体.局部法的核心思想是以断层面分割地层，分别构建两侧的三维地质模型，再合并为完整的三维地质模型.肖义等（2020）提出了一个扩展工具AHGW，该工具利用局部法将地质体沿断层进行划分，实现含断层的三维地质模型构建.该方法不需要用户过多干涉，但其适应性差，不适用于含复杂断层网络模型的构建.朱良峰和潘信（2008）提出了断层与地层的统一建模法，该方法基于断层数据特点对地层层面进行整体拟合，同时顾及断层线的约束作用，然后用断层面切割地层面，调整交线上的点，从而形成完整的实体模型.该方法避免了整体法和局部法的局限性，具有较强的柔韧性，但它需要足够多的控制数据才能构造任意复杂的断层系统，并且该方法未考虑断层时序性的影响.地层变形位移法用一个三维矢量来量化地层位移.基于该方法，Havana数值工具（Hollund et al.，2002）提供椭圆型断层模型和参数故障模型（PFM）两种模型，其中椭圆型断层模型主要适用于低于地震分辨率且断层面为椭圆的小断层，而参数故障模型通常适用于模拟地震数据中一致或部分一致且断层面通过一系列双线性平面进行表示的较大断层.该方法提供了完整的断层表示模型和地层移动变形算子，但滑动面两侧位移对称，与地质现象相违背，并且该方法适用于已经建立的油藏网络，难以在模型中添加几何断裂网络.

复杂断裂网络模型能够表达多条断裂间的几何形态和空间关系（Thore et al.，2002； Caumon et al.，2009），清晰地反映不同断裂在三维空间中的展布形态和拓扑结构.因此，许多学者也基于不同建模方法，实现了复杂断裂网络的表达方法.李培军（2000）、方燕等（2003）提出一种层状地质体构建断裂模型，首先建立三维断裂网络，再对断裂两侧插值地层面，但由于多值问题，难以建立复杂形态的断裂模型.武强和徐华（2005）利用一个或多个平面组合的方式构建断裂面，使用实测的走向、倾向和倾角等表示每条断裂的空间方位，并用一系列的断点数据构建断裂面的平面方程，但这种方法在处理断裂的相交相切时，很难发挥实际作用.主流商业建模软件大多采用Pillar的方法 （李兆亮等，2015），用几个断裂点来控制断裂面的空间展布，但对于一些大型断裂来说，几个断裂点难以描述其空间形态，所以Pillar方法难以表达断裂的真实形态.

基于上述分析，大部分建模方法都没有考虑断层构造时序性及其对地层变形的影响，本文在整体法和地层变形位移法的基础上，引入断裂恢复与演化影响域理论，提出了断裂矢量场建模方法，将断裂构造的时序性及其对影响域范围内地层的变形考虑到三维地质建模中，从而有效提高了建模效率、数据利用率以及处理断裂接触时的能力，为含复杂断层的三维地质建模提供一种新的技术方法.

1 断裂矢量场与构造演化理论

1.1　断裂恢复与演化理论

断裂恢复与演化理论（Godefroy et al.， 2018）是基于年代地层学相关知识的启发提出的：从数学层面上定义一个新的空间，在这个空间里，所有的地层都是水平的，而断裂（如果有的话）则已经消失了.继而推测，无论这个新空间如何变化，无论后续使用何种插值方法，都更适合模拟地下的物理性质.该数学框架还能提供复杂问题的解决方案，如确定由构造事件引起的变形，可以将所有构造事件当作对初始地层的变形，随着地质演化，每次构造运动都会对上一次的构造运动产生影响，以便重现地质历史和相关变形.对于同一条断裂在多个时期进行构造运动的情况，本文实现的算法只计算一次该断裂的导致影响，即多次构造运动对地层的综合影响.

1.2　断裂矢量场

断裂矢量场定量描述了断裂的影响范围内地层的位移.在断裂矢量场中，矢量场位移强度由断距衰减函数计算得到，矢量场位移方向由广义径向基函数（RBF）（Carr et al.，2001； Cowan et al.，2003； Hillier et al.， 2014）插值断裂面梯度旋转得到，最后用矢量场位移算子定量表示地层变形.

1.2.1　矢量场断距衰减函数

本文中，断裂面被认为是一个零厚度的滑动面，并且构造的矢量场在这个滑动面的某个区域上是最大的，这个区域称为断裂中心.它还需要在某一边界之外为零，并在这两个区间内慢慢减小.因此，这个区域内的断距衰减函数由断裂中心、断裂走向、法向、倾向长度而确定.Barnett et al.（1987）及Walsh and Watterson（1987）最早利用在煤矿中获得的断裂数据，定义了孤立断裂的空间位移分布.断裂的影响范围为一个椭球体，在断裂中心附近具有最大的位移量，并且位移强度从断裂中心向四周逐渐减小，类似钟形曲线，该理论模型得到了在地球物理方法下获取的孤立断裂数值模型的证实（Maerten et al.，1999）.本文选取的就是该理论模型中的断距衰减函数，三维空间中的衰减函数如下：

\begin{array}{l}{f}_{\mathrm{0}}\left(x,y,z\right)=\mathrm{2}\cdot \left(\mathrm{1}-r\left(x,y,z\right)\right)\cdot \\ \sqrt[]{\frac{{\left(\mathrm{1}+r\left(x,y,z\right)\right)}^{\mathrm{2}}}{\mathrm{4}}-r{\left(x,y,z\right)}^{\mathrm{2}}},\end{array} (1)

式中： r \left(x,y,z\right) = \sqrt[]{{(\frac{x-x_{\mathrm{0}}}{l_x})}^{\mathrm{2}}+{(\frac{y-y_{\mathrm{0}}}{l_y})}^{\mathrm{2}}+{(\frac{z-z_{\mathrm{0}}}{l_z})}^{\mathrm{2}}}是到断裂中心点的径向距离，位于椭圆边界上 r \left(x,y,z\right) = \mathrm{1}以及该边界之外，衰减函数值为零. l_x、 l_y及 l_z为断裂走向、倾向及法向长度.

根据断裂上下盘相对运动的非对称特性，引入断裂上下盘位移比（ \gamma），这个比值可从断裂倾角估计得到，同时，这个比值的取值区间为［0，1］，当该比值为0和1时表示所有位移分别发生在断裂下盘和上盘.因此最终断距衰减函数可以定义为：

（1）上盘断距衰减函数：

f \left(x,y,z\right) = \mathrm{\gamma } \cdot d_{\mathrm{m}\mathrm{a}\mathrm{x}} \cdot f_{\mathrm{0}} \left(x,y,z\right) . (2)

（2）下盘断距衰减函数：

f \left(x,y,z\right) = \left(\mathrm{1}-\mathrm{\gamma }\right) \cdot d_{\mathrm{m}\mathrm{a}\mathrm{x}} \cdot f_{\mathrm{0}} \left(x,y,z\right) . (3)

其中断距（最大位移长度） d_{\mathrm{m}\mathrm{a}\mathrm{x}}与椭球三轴长度由实际情况而确定，通过断裂认识的更新，可以灵活地修改各参数的数值，从而建立更加符合地质学家认知的断裂模型.

1.2.2　断裂矢量场位移方向

通过增加矢量方向约束，保证矢量场中变形地层运动方向的一致性，即断裂上下盘必须始终沿着断裂面运动，并且不能穿过断裂面.同时，为了实现断裂两侧运动方向相反，使用隐式插值函数表示一个给定的点p位于断裂的下盘还是上盘，这取决于点p的断裂势场值是小于还是大于断裂的等值面值.隐式构造建模过程中，利用断裂数据可以生成一个具有断裂表面几何形状的标量场，在该标量场中，具有以下特征：（1）断裂面对应标量场中的曲面，通常该曲面的属性值设为0；（2）上盘的值大于0，下盘的值小于0；（3）该标量场中的梯度总是垂直于断裂面.

断裂插值的方法与地质界面相同：使用断裂点作为数据点和断裂方向数据插值势场.因此，每个断裂都是由其自身的势场确定的，断裂的几何形状由断裂势场提取的等值面确定.利用广义径向基函数插值出的断裂标量场得到断裂梯度场 g ( x , y , z )，再根据不同运动类型的断裂对该梯度场进行相应的旋转操作，从而得到断裂矢量位移方向场 s ( x , y , z ).针对不同运动类型断裂计算断裂矢量位移场的方法如下：

（1）正断裂：利用断裂标量场得到断裂梯度场 g \left(x,y,z\right)，将方向矢量转化为该点的倾角 \theta和倾向 \alpha.上盘沿着断面向下运动，且 g ( x , y , z )垂直于断裂面，只需要将倾角 \theta加90°，下盘沿断裂面向上运动，运动方向与上盘相反，只需要对上盘位移向量取反.

\theta = a \cdot \mathrm{c} \mathrm{o} \mathrm{s} \left(g\left(z\right)\right) + \frac{\mathrm{\pi }}{\mathrm{2}} , \mathrm{\alpha } = \mathrm{\alpha } \cdot \mathrm{s} \mathrm{i} \mathrm{n} \left(\frac{g\left(x\right)}{\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)}\right)

\mathrm{上盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\\ \mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\\ \mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\end{array}\right] , (4)

\mathrm{下盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}-\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\\ -\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\end{array}\right] . (5)

（2）逆断裂：逆断裂上下盘的运动方向与正断裂相反.

\mathrm{上盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}-\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\\ -\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\end{array}\right] , (6)

\mathrm{下盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}\left(\alpha \right)\\ \mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}\left(\alpha \right)\\ \mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\end{array}\right] . (7)

（3）左移断裂：同理正断裂，上盘沿着断面向左运动，且 g ( x , y , z )垂直于断裂面，只需要对倾向 \alpha旋转90°.倾向的取值范围在0到 \mathrm{2} \mathrm{\pi }之间，需要对不同情况进行不同处理.

，

\theta = a \cdot \mathrm{c} \mathrm{o} \mathrm{s} \left(g\left(z\right)\right) , \alpha = a \cdot \mathrm{s} \mathrm{i} \mathrm{n} ( g ( x ) / \mathrm{s} \mathrm{i} \mathrm{n} ( \theta ) )，

\mathrm{上盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}(\alpha )\\ \mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}(\alpha )\\ \mathrm{c}\mathrm{o}\mathrm{s}(\theta )\end{array}\right] , (8)

\mathrm{下盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}-\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}(\alpha )\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}(\alpha )\\ -\mathrm{c}\mathrm{o}\mathrm{s}(\theta )\end{array}\right] . (9)

（4）右移断裂：右移断裂上下盘的运动方向与左移断裂正好相反.

\mathrm{上盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}-\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}(\alpha )\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}(\alpha )\\ -\mathrm{c}\mathrm{o}\mathrm{s}(\theta )\end{array}\right] , (10)

\mathrm{下盘} \mathrm{位移} \mathrm{向量} : \left[\begin{array}{c}{s}_x\\ s_y\\ s_z\end{array}\right] = \left[\begin{array}{c}\mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{s}\mathrm{i}\mathrm{n}(\alpha )\\ \mathrm{s}\mathrm{i}\mathrm{n}\left(\theta \right)\cdot \mathrm{c}\mathrm{o}\mathrm{s}(\alpha )\\ \mathrm{c}\mathrm{o}\mathrm{s}(\theta )\end{array}\right] . (11)

1.2.3　矢量场位移算子

断裂位移算子用于表示地质地层变形，它不仅表示断裂几何形状，还能够表达位移.这种算子已被用于模拟具有地震断裂的地层，并在随机建模中用于修改地层模型，以及在模型构建过程中融入地质时序知识.本文使用的算子使法向上梯度在断裂算子极限处为零，因为位移远离断裂平面平滑地减小.当位移为零时，断裂消失，因此走向和倾向上在极限处为零.在最简单的（孤立断裂）情况下， D_s（走向方向上位移）和 D_d（倾向方向上位移）是连续的，最大值靠近断裂中心. D^f方向上（与断裂正交）在断裂中心处是不连续的，以说明断裂上变形场的不连续性. D^f方向上可以是不对称的，以解释脆性断裂破坏带中通常观察到的不对称.使用这些剖面对应的衰减函数和不同断裂对应的位移方向向量可以得到断裂矢量场的位移算子 D^f ( x , y , z )：

（1）上盘位移算子：

\left[\begin{array}{c}{D}_x^f\\ D_y^f\\ D_z^f\end{array}\right] = \left[\begin{array}{c}{D}_x\\ D_y\\ D_z\end{array}\right] + \left[\begin{array}{c}{S}_x\\ S_y\\ S_z\end{array}\right] \cdot f \left(x,y,z\right) . (12)

（2）下盘位移算子：

\left[\begin{array}{c}{D}_x^f\\ D_y^f\\ D_z^f\end{array}\right] = \left[\begin{array}{c}{D}_x\\ D_y\\ D_z\end{array}\right] - \left[\begin{array}{c}{S}_x\\ S_y\\ S_z\end{array}\right] \cdot f \left(x,y,z\right) , (13)

式中： D^f ( x , y , z )为矢量场位移后的地层点， D ( x , y , z )为矢量场中未发生位移的地层点， S ( x , y , z )为矢量场位移方向向量， f \left(x,y,z\right)为断距衰减函数， S ( x , y , z )及 f \left(x,y,z\right)与断裂本身几何形状及位置有关.

2 含复杂断层的三维地质模型构建

2.1　模型构建流程

基于断裂矢量场建模的方法本身只适用于孤立断裂建模，在建立复杂断裂时难以将断裂影响域中的数据点约束到建模过程中.针对该问题，在断裂矢量场的基础上，本文基于断裂恢复与演化影响域理论提出断裂矢量场建模方法用于含复杂断裂网络的三维地质模型构建.在构建三维地质模型前，首先采用隐式建模方法中的广义径向基插值法构建有断层的初始三维地质模型.然后采用本文提出的断裂矢量场建模方法构建含复杂断层的三维地质模型，构建流程如下：（1）断裂时序构建：结合地质资料、断裂构造图等，根据地层层序律和地质体间的切割律确定断裂发生的相对时间顺序；（2）断裂构造恢复：基于断裂构造恢复理论，使用矢量场位移算子以断裂构造演化逆序依次恢复断裂和地层数据，此时数据已经避免了断裂所造成的影响；（3）三维地质模型构建：按照断裂构造演化正序逐条加入断层，最后生成具有复杂断裂几何形状和接触类型的三维地质模型.

2.2　断裂时序构建

地质构造的连续性决定了断裂的特征及其相互关系，构造历史是有关断裂网络的重要信息来源，这一点必须在断裂网络建模中加以考虑.在断裂网络时序构建时，通常使用一颗断裂二叉树来表示断裂对建模空间的划分，从而来反映断裂之间的时间及空间关系，使得断裂网络模型的更新更加繁琐，并且该方法必须始终选择一个断裂作为主断裂，因此，多条依次为主辅关系断裂的情况不可能使用二叉树模型.在本文提出的方法中，只需要构建断裂构造的相对时间关系，而不需要表达断裂的空间关系，从而优化了断裂时序构建过程中的复杂度.

断裂时序构建就是要确定各个断裂的相对地质年代，即确定断裂发生的相对时间顺序，主要有以下两种方法来判断地质体（岩体、岩层、断裂）的新老关系：①基于地层层序律.未经强烈构造变动的地层，其沉积顺序总是上新下老，因此可以判断一般沉积地层的构建顺序为从下至上，从老至新.②基于地质体间的切割律.经过构造和岩浆活动，使得不同时代的地质体之间出现彼此切割关系，新形成的地质体总是切割老的地质体.岩体与岩层之间的沉积接触和侵入接触关系、岩层之间的不整合接触关系、岩体之间的切割关系，对于断裂构造其具体表现如图1所示，断裂A被断裂B所切穿，这就意味着断裂A先于断裂B产生.从地质历史演化角度出发，两次断裂的运动学相互作用是不同的，断裂A对矢量场空间中的岩体初次造成地层变形及位移，断裂B再在断裂A的基础上产生了一次平移运动，导致断裂A上下盘断裂成两部分，分别向不同方向运动，即地层发生两次矢量场位移变形.在基于断裂和地层连续插值的建模方法中，很难明确地考虑两次构造运动.因此，本文方法中需要记录断裂构造发生的相对时间顺序，方便后续使用矢量场位移算子进行运动学建模，按照断裂构造发生顺序模拟断裂网络，这一步也直接决定了断裂网络的拓扑关系是否正确.

2.3　断裂构造恢复

在断裂构造恢复过程中，基于断裂恢复与演化影响域理论和断裂矢量场，利用上文的公式求解矢量场位移算子，使用矢量场位移算子以断裂构造时序逆序依次恢复地层，即利用矢量场位移后的地层点 D^f ( x , y , z )求解矢量场中未发生位移的地层点 D ( x , y , z ).矢量场位移算子本质上是一个非线性方程组，不存在相应的求根公式，因此本文使用牛顿迭代法（成春蕾，2022），即把非线性方程组线性化的一种近似方法寻求方程的近似根，不断迭代用旧的 D^f ( x , y , z )推导新的 D ( x , y , z ).牛顿迭代法已经被证实：如果是连续的，并且待求的零点是孤立的，那么在零点周围存在一个区域，只要初始值位于这个邻近区域内，那么牛顿法必定收敛.本文使用C++GSL提供的Multidimensional Root-Finding组件求解非线性系统，该库为各种迭代求解器和收敛提供了低级组件，可以将它们组合起来以实现所需的解决方案，同时可以完全查询迭代的中间步骤.

如表1中的伪代码所示，完成断裂构造时序构建后，首先从最新的断裂开始，先计算新断裂的矢量场位移算子和标量场插值函数.由于断裂上下盘的位移函数不同，因此先使用标量场插值函数判断初始模型或经过上一条断裂恢复后的断裂和地层数据是位于新断裂的上盘还是下盘，再使用矢量场位移算子恢复断裂和地层数据.之后重复上述步骤，用较新断裂矢量位移场按照断裂构造时序逆序恢复前一步的断裂和地层数据，直至所有的断裂构建完毕.最后，所有的断裂和地层数据都已恢复到初始状态（未被断裂构造影响），即完成了断裂构造恢复.

2.4　三维地质模型构建

最后，在已经避免了断裂所造成影响的数据基础上，按照断裂时序正序利用矢量场位移算子调整断裂和地层数据.断裂一个接一个地产生，新产生的断裂是先前断裂的分支断裂，通过重复这一过程，断裂以不同的顺序绘制，并以不同的方式聚集，最终生成具有不同断裂几何形状和接触类型的断裂的三维地质模型.

这种基于构造演化的运动学模拟方法有非常强的应用前景，因为它允许考虑构造的不确定性，并从一个独特的数据集生成无限等价的一致构造模型，并且这种基于断裂构造演化的建模方法，建模并不是一次性完成，而是先建立初始地层模型，再在模型上初步叠加断裂构造事件，逐事件构建模型，这样就能保存建模的中间过程，有助于建模人员在建模过程对断裂参数进行调整.在其他方法中，受最近变形事件影响的断裂的每个部分都是独立表示的，需要手动调整断裂两侧模型运动的一致性，而在本文的方法无需手动调整断裂两侧位移，断裂导致的一切位移都是自动计算，从而提高了建模效率.

3 含复杂断层三维地质模型构建实例

3.1　建模数据

为了验证上述的相关理论与断裂矢量场建模方法的有效性与实用性，笔者基于C++开发了三维地质建模与可视化原型系统.在PC平台下（CPU：Intel（R） Core（TM） i7-7700 3.60GHz；内存：8GB；操作系统：Windows 10 专业版），以澳大利亚珀斯盆地的部分区域为研究区（Florian Wellmann and Guillaume Caumon， 2018），建立该区域的三维地质模型.珀斯盆地是北向西北走向的陆上和近海沉积盆地，沿澳大利亚大陆西南边缘延伸约1 300 km，是一个大型构造复杂的盆地，形成于二叠纪至白垩纪的澳大利亚板块和印度板块分离期间.它包括一个重要的陆上部分，并向近海延伸至大陆壳的边缘，珀斯盆地构造是二叠纪，三叠纪晚期至侏罗纪早期，中侏罗纪至早白垩期间的产物，叠加在已有的基底地形上，二叠纪期间的延伸沿西边缘产生了一系列深（长达15 km）的南北向裂谷盆地.如图2所示，本文选择该研究区内四个沉积地层（表2）共61个地层采样点和11个产状采样点，8条断裂（表3）共33个采样点和8个产状采样点进行建模.

3.2　建模实验

研究区地质演化如下：最先沉积形成4个初始地层，后期再经8次断裂构造运动的变形错动，形成了目前复杂的断裂网络布局.如表2，表3所示，地层的沉积序列及断裂的构造时序表明了断裂构造事件发生的时间顺序.

基于本文的断裂矢量场建模方法，主要分为以下两个步骤：（1）先计算每条断裂的矢量场位移算子，按照断裂构造时序演化逆序对每条断裂构造进行数据恢复.（2）按照断裂构造时序演化正序逐条加入断裂，直到重建完所有断裂，从断裂事件演化中解决断裂之间的拓扑关系（图3）.如图4a展示了在多个断裂影响下地层变形的一致性，所有地层均沿着断裂面方向位移，不会发生相交、穿模等拓扑错误.图4b展示了方法在处理断裂接触时的能力，南北向断裂将两条东西向平行断裂切割，并造成断裂两盘的相对位移，这种现象不仅符合本建模方法的理论基础，也满足该研究区的先验知识，表明方法可以构建出满足地学规律并且具有良好可视化效果的含复杂断层的三维地质模型.

3.3　对比实验

将本文的断裂矢量场位移方法与Bouziat（2012）提出的基于形状变形的矢量场方法（Vector Field based Shape Deformations）进行对比（图5）.由于VFSD方法只解决了孤立断裂建模问题，故该实验只选择研究区中的一条断裂参与建模实验.

在VFSD方法中，断裂影响域中的数据不参与建模，导致断裂影响域中的地层曲面与实际数据不吻合，如模型最上方的地层未进行错断，最下方地层在断层面附近形态失真，地层厚度不符合实际情况等（图5a）.在处理具有复杂接触关系的断裂网络建模中，该方法不仅没有涉及复杂断裂，而且还需要剔除位于断裂影响域中的数据才能够建立符合理论的模型.因此，导致模型的精度大大降低，甚至难以处理上述复杂的断裂模型.与VFSD方法相比，本文的方法利用断裂构造恢复理论，将断裂影响域中的数据约束到建模过程中，提高了数据的利用率，并且在构建的模型中断层对地层的错断更加明显，断层上下盘的地层形态更加合理（图5b）.本文的方法同样具有很强的灵活性，充分考虑了断裂构造时序演化，如上一节的建模实验中展示了基于矢量场的方法在处理各种类型断裂（正断裂、逆断裂、平移断裂、贯穿及非贯穿断裂）的能力，在含复杂断裂网络建模中具有更强的适用性.

4 结语

许多三维地质建模工作集中于断裂本身的模拟，而忽略了对断裂错断地层的模拟，为了模拟断裂对地层造成的位移变形，本文提出了断裂矢量场建模方法，利用矢量场位移算子使仅限于邻近断裂影响域中的地层产生变形，构建运动学一致的构造模型.该方法弥补了矢量场构建断裂时无法利用受影响数据的缺点，使得更多地质约束参与到建模过程中，提高了数据的利用率.而且该方法不仅考虑了断裂对周围岩层造成的位移变形，还在建模过程中考虑构造演化的时间顺序，从而能够解决具有复杂接触关系的断裂网络模型构建问题，更适用于含复杂断层的三维地质模型的构建.在其他方法中，受最近变形事件影响的断裂的每个部分都是独立表示的，因此，需要手动调整断裂两侧模型运动的一致性，而在本文提出的方法中，无需手动调整断裂两侧位移，断裂导致的一切位移都是自动计算，从而提高了建模效率.

尽管文本提出的方法能够实现部分断裂的建模，但仍然存在一定的不足：①只考虑了断裂构造造成的位移，但实际上其他构造（如褶皱、底辟构造等）也会造成地层位移，因此后续工作中需要考虑更多构造产生的位移.②未考虑断裂本身的复杂性，本文提出的方法建立的断裂模型的位移方向始终平行于断裂面，虽然保证了地层在位移变形的一致性不会发生拓扑错误，但断裂造成的运动方向并不总是平行于断裂面，因此下一步工作中需要更加精确的位移矢量来约束模型的运动.③本文未对断裂模型进行不确定性评价、分析不确定性产生的原因，为模型进一步优化提供参考.

参考文献

[1] Barnett, J., Mortimer, J., Rippon, J., et al., 1987. Displacement Geometry in the Volume Containing a Single Normal Fault. AAPG Bulletin, 71(8):925-937. https://doi.org/10.1306/948878ed-1704-11d7-8645000102c1865d

[2] Bouziat, A., 2012. Vector Field Based Fault Modelling and Stratigraphic Horizons Deformation. In: Proceedings of the 31st Gocad Meeting, San Antonio, Texas, USA. SPE. https://doi.org/10.2118/160905-stu

[3] Carr, J. C., Beatson, R. K., Cherrie, J. B., et al., 2001. Reconstruction and Representation of 3D Objects with Radial Basis Functions. In: Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques. ACM, 67-76. https://doi.org/10.1145/383259.383266

[4] Caumon, G., Collon-Drouaillet, P., Le Carlier de Veslud, C., et al., 2009. Surface-Based 3D Modeling of Geological Structures. Mathematical Geosciences, 41(8): 927-945. https://doi.org/10.1007/s11004-009-9244-2

[5] Cheng, C.L., 2022. Newton Iteration Method and Its Implementation by MATLAB. China BroadBand, (3): 40-42 (in Chinese).

[6] Cowan, J., Beatson, R., Ross, H.J., et al., 2003. Practical Implicit Geological Modelling. AusIMM Publication Series, Bendigo, Victoria, Australia, 89-99.

[7] Fang, Y., Li, M., Hu, Y.J., 2003. Three-Dimensional Visualization of Stratified Geological Objects. Journal of Jiaozuo Institute of Technology, 22(6): 441-444 (in Chinese with English abstract).

[8] Godefroy, G., Caumon, G., Ford, M., et al., 2018. A Parametric Fault Displacement Model to Introduce Kinematic Control into Modeling Faults from Sparse Data. Interpretation, 6(2): B1-B13. https://doi.org/10.1190/INT-2017-0059.1

[9] Hillier, M. J., Schetselaar, E. M., de Kemp, E. A., et al., 2014. Three-Dimensional Modelling of Geological Surfaces Using Generalized Interpolation with Radial Basis Functions. Mathematical Geosciences, 46(8): 931-953. https://doi.org/10.1007/s11004-014-9540-3

[10] Hollund, K., Mostad, P., Nielsen, B.F., et al., 2002. Havana; A Fault Modeling Tool. In: Koestler, A.G., Hunsdale, R., eds., Hydrocarbon Seal Quantification. Norwegian Petroleum Society Conference, Norway, 157-171. https://doi.org/10.1016/S0928-8937(02)80013-3

[11] Hua, W.H., Guo, D.Y., Liu, X.G., et al., 2023. Unified Correction and Connection Method of Stratigraphic Sequence with Complex Inversion. Earth Science, 48(4): 1532-1542 (in Chinese with English abstract).

[12] Li, C.L., Zhang, H., Zhu, L.F., 2013. Algorithm for True-3D Modeling of Geological Body with Single-Fault. Computer Engineering and Design, 34(10): 3590-3594 (in Chinese with English abstract).

[13] Li, P.J., 2000. Three-Dimensional Simulation and Visualization of Layered Geological Bodies. Earth Science Frontiers, 7(S2): 271-277(in Chinese with English abstract).

[14] Li, Z.L., Pan, M., Yang, Y., et al., 2015. Research and Application of the Three-Dimensional Complex Fault Network Modeling. Acta Scientiarum Naturalium Universitatis Pekinensis, 51(1): 79-85 (in Chinese with English abstract).

[15] Maerten, L., Willemse, E. J. M., Pollard, D. D., et al., 1999. Slip Distributions on Intersecting Normal Faults. Journal of Structural Geology, 21(3): 259-272. https://doi.org/10.1016/s0191-8141(98)00122-9

[16] Sha, H., 2016. Fault 3D Geological Modeling Based on Integral Method. Coal and Chemical Industry, 39(8): 24-26 (in Chinese with English abstract).

[17] Tai, W.X., Zhou, Q., Yang, C.F., et al., 2023.3D Geological Visualization Modeling and Its Application in Zhexiang Gold Deposit, Southwest Guizhou Province. Earth Science, 48(11): 4017-4033 (in Chinese with English abstract).

[18] Thore, P., Shtuka, A., Lecour, M., et al., 2002. Structural Uncertainties: Determination, Management, and Applications. Geophysics, 67(3): 840-852. https://doi.org/10.1190/1.1484528

[19] Walsh, J. J., Watterson, J., 1987. Distributions of Cumulative Displacement and Seismic Slip on a Single Normal Fault Surface. Journal of Structural Geology, 9(8): 1039-1046. https://doi.org/10.1016/0191-8141(87)90012-5

[20] Wang, Z., Liu, Z., An, C.R., et al., 2015.3D Fault Modeling Method Based on Profile. Journal of Xi’an Shiyou University (Natural Science Edition), 30(6): 50-54, 61, 9(in Chinese with English abstract).

[21] Wellmann, F., Caumon, G., 2018.3-D Structural Geological Models: Concepts, Methods, and Uncertainties. Advances in Geophysics. Elsevier,Amsterdam, 1-121. https://doi.org/10.1016/bs.agph.2018.09.001

[22] Wu, Q., Xu, H., 2005. Technique for 3D Fault Simulation in Virtual Mine System. Journal of Liaoning Technical University, 24(3): 316-319 (in Chinese with English abstract).

[23] Xiao, Y., Guo, T.T., Xu, S.G., et al., 2020. Three-Dimensional Geological Modeling of Faults Based on AHGW. China Water Transport, 20(6): 175-177 (in Chinese with English abstract).

[24] Zhang, J.Y., Cui, Z.D., Han, W.G., et al., 2021. Research Progress of 3D Fault Modeling Technology Considering the Process of Geological Evolution. Xinjiang Geology, 39(2): 336-339 (in Chinese with English abstract).

[25] Zhu, L.F., Pan, X., 2008. Study on Computer Modeling of Geological Faults in 3D. Rock and Soil Mechanics, 29(1): 274-278 (in Chinese with English abstract).

[26] 成春蕾, 2022. 牛顿迭代法及其MATLAB实现. 中国宽带, (3): 40-42.

[27] 方燕, 李梅, 胡友健, 2003. 层状地质体的三维可视化研究. 焦作工学院学报(自然科学版), 22(6): 441-444.

[28] 花卫华, 郭丹阳, 刘修国, 等, 2023. 含复杂倒转的地层层序统一修正与连接方法. 地球科学, 48(4): 1532-1542.

[29] 李昌领, 张虹, 朱良峰, 2013. 单断层地质体真三维建模算法. 计算机工程与设计, 34(10): 3590-3594.

[30] 李培军, 2000. 层状地质体的三维模拟与可视化. 地学前缘, 7(增刊2): 271-277.

[31] 李兆亮, 潘懋, 杨洋, 等, 2015. 三维复杂断层网建模方法及应用. 北京大学学报(自然科学版), 51(1): 79-85.

[32] 沙浩, 2016. 基于整体法的断层三维地质建模. 煤炭与化工, 39(8): 24-26.

[33] 邰文星, 周琦, 杨成富, 等, 2023. 黔西南者相金矿床三维地质可视化建模及应用. 地球科学, 48(11): 4017-4033.

[34] 王诏, 刘展, 安聪荣, 等, 2015. 基于剖面的三维断层建模方法. 西安石油大学学报(自然科学版), 30(6): 50-54, 61, 9.

[35] 武强, 徐华, 2005. 虚拟矿山系统中三维断层模拟技术. 辽宁工程技术大学学报, 24(3): 316-319.

[36] 肖义, 郭婷婷, 徐世光, 等, 2020. 基于AHGW的含断层三维地质建模. 中国水运(下半月), 20(6): 175-177.

[37] 张建勇, 崔振东, 韩伟歌, 等, 2021. 考虑地质演化过程的断层三维建模技术研究进展. 新疆地质, 39(2): 336-339.

[38] 朱良峰, 潘信, 2008. 地质断层三维构模技术研究. 岩土力学, 29(1): 274-278.

基金资助

国家重点研发计划课题(2019YFC0605102)

国家自然科学基金项目(41972307)