引言

致密油储层裂缝的精细描述和建模是当前石油地质界的世界性难题之一，目前仍处于探索阶段。实际岩层的裂缝受控于地层局域应力场、岩石属性及其空间分布^[1-3]。如果不分析构造演化历史就进行裂缝分布预测，则所建裂缝模型的可信度必然大大降低。Robinson 和Sahim 等提出了基于统计理论建立随机裂缝网络模型的方法^[4-5]，并在20 世纪80 年代和90 年代的岩石工程研究中得到了广泛应用^[6]。

裂缝随机网络模型的基本思路是将裂缝的位置、尺度、方向及其他一些与裂缝相关的特征视为随机变量，并赋予一定的概率分布函数。为了降低模型的复杂度，建模时通常利用简单的几何形状表示裂缝，如二维模拟时用直线表示，三维模拟时用椭圆或平面表示^[7]。由于能够采集到的裂缝统计数据，如露头、微电阻率成像测井、井下电视成像测井等，多为二维数据，且受到截断误差、边界效应等干扰，不能直接进行三维裂缝建模^[8-13]。在随机模拟时，裂缝的位置是裂缝随机建模中最重要的参数之一，通常独立于裂缝的其他属性而进行单独处理。Young 等提出基于地质统计学的裂缝位置模拟方法，选用裂缝密度作为模拟对象^[14]。该方法对于一个给定的裂缝数据集，通过逐点计算，首先确定裂缝位置空间相关函数；然后，利用克里金估计获取裂缝位置的无偏估计，这个过程可以协同地震等其他参数进行模拟^[15]；最后，基于特定的概率分布函数，再对裂缝的几何特征进行模拟。以上这些裂缝随机模拟方法，也称为标示点过程，即在进行裂缝的位置用一个点代替，而裂缝的其他属性，则当做该点上的各种标示^[16]。

裂缝的随机模拟方法将裂缝和断层视为简单几何形状的集合，但在实际应用中，并不容易获得该方法所需要的裂缝位置、长度、宽度等几何特征的真实概率分布函数。类似的随机建模方法，如基于像素的多点统计建模方法，不仅需要大量的训练图片，且可靠性无法确定，无法进行三维裂缝建模。从地质力学的角度考察，裂缝随机建模方法的主要弊端就是无法保证其符合地质事实，其中一个主要原因就是，裂缝的随机建模理论并没有考虑裂缝的成因，而仅依靠现今的裂缝系统统计量建模。

Egan 和Tanner 等提出了基于地质几何力学的裂缝建模方法，利用几何算法，模拟主断层的运动轨迹^[17-18]。在进行三维模拟时，通常利用3D Move软件对主断层进行断裂恢复，直到断层上下盘重合^[19]。断层的恢复方向通常选定为垂直褶皱轴线方向^[20-21]，并根据几何算法不同，略有差别。恢复后的层面模型，即被认定为一个可接受的变形前模型，而构造正演只需对恢复过程进行回放即可获得。在构造正演过程中，对比最终构造与当前模型之间的差值，可以获得绝对应变积分值^[22]。利用绝对应变积分可以求出一个绝对应力数值，就可以建立裂缝的位置和方向模型。从断裂力学的角度考察，基于几何力学的裂缝建模方法，并没有考虑裂缝尖端奇性场，因此，无法保证模拟满足断裂力学规律，也无法保证所得应变场的合理性^[23]。同时，在进行几何力学模拟时，没有考虑主断层的非均质扩展，也无法分析岩性非均质引起的应力集中所产生的小裂缝。

在传统均质裂缝研究中，一般假定裂缝为光滑曲面，然而越来越多的实验表明，事实并非如此。Golʹdshtejn 和Yavari 等研究了分形裂缝的力学机制，试图用分形维空间来解释裂缝的尖端奇性^[24-26]。Golʹdshtejn 等从裂缝的自相似角度出发，得到了同Yavari 同样的结论^[27]。Balankin 等也对这一理论进行了大量研究^[28-30]。目前，有些研究者开始认为，裂缝的分形维特征的本质在于介质本身的属性特征，如材料刚度。但由于边界条件设定带来的困难，基于分形维进行裂缝模拟，并没有取得比光滑裂缝模拟更好的结果。

本文提出的基于随机扩展有限元模拟的地质裂缝建模方法，将地质统计的概率观点和断裂模拟的动力学原理及其微观断裂机理结合在一起，构建更符合地质事实的裂缝模型。该方法共包含4 个步骤：（1）收集与裂缝有关的统计数据，包括区域构造演变历史、构造解释数据、钻井取芯数据及各种成像测井数据等。（2）建立属性模型，基于地质统计学建立属性模型，为裂缝模拟提供基础模型，要求基础模型要大于目标研究区的尺度，以便减小边界效应对研究区裂缝建模的影响。由于储层岩石高度非均质，其应力强度因子的解析形式不易获得，通常只能根据大量不同孔隙度和几何形状的裂缝进行裂缝扩展实验结果。（3）根据属性模型，基于扩展有限元方法，进行主断层扩展历史拟合。（4）利用多场趋势的克里金估计方法，建立裂缝特征最优估计模型。

1 裂缝扩展反演模型

裂缝演化系统是一个复杂的非线性动力系统^{[20-23, 28-32]}。为了考察裂缝的反演分析及分布预测问题，建立储层裂缝扩展正向演化模型

$\dot{X}=G(X,{{X}_{0}},{{P}_{0}})$

(1)

式中：

X—裂缝状态量；

X₀—裂缝扩展初始状态量；

P₀—地质系统的属性特征、初始条件、边界条件和断裂力学特征等裂缝演化系统参量。

对于单一裂缝演化系统，采用统一的裂缝扩展规律G，则裂缝演化预测系统的差分形式可以转化为裂缝扩展的反演问题，并可定量地写为以下最小化问题

${{J}^{i}}=\frac{1}{N}\sum\limits_{k=1}^{N}{||X_{k}^{i}-{{Y}_{k}}||},i\in [1,M]$

(2)

$P=\min {{J}^{M}}$

(3)

式中：

M—模拟步总数；

N—裂缝系统参量估计时所用裂缝属性个数；

Y_k—通过岩芯、测井和地震等资料获得的第k种裂缝状态量；

X_kⁱ—第i模拟步时预测的第k种裂缝状态量；

P—模型的最优解。

由于裂缝扩展演化系统属于非线性系统，反演会出现多个在局部最优解。为了搜索全局最优解，本文采用模拟退火技术^[33-34]。模拟退火方法是通过模拟熔融金属冷却系统，寻找优化问题全局最优解的方法。

在裂缝演化系统中，任意演化时刻的系统能量记为$E({t_i})= {J^i}$，下一个时刻的系统能量依赖于两个时间点间的能量差($\Delta E({t_i},{t_{i + 1}})= E({t_{i + 1}})-E({t_i})$)。 裂缝反演系统各个状态出现的概率$P_{E({t_{i + 1}})}= \min(1,-{\Delta E({t_{i + 1}},{t_i})}/{kT})$。如果$\Delta {J^i}$＜$0$，则接收第i步的裂缝扩展模拟结果，否则以概率$\exp(-\Delta {J^i})$接收模拟结果，其中$\Delta {J^i} = {J^{i-1}}-{J^i}$。

2 裂缝扩展正向演化模型 2.1 基本控制方程

模型边界条件采用准静态加载方式，将岩层应变控制方程写成虚功原理的形式

$\delta w=0$

(4)

式中：

w—总能量。

对于线弹性材料，w 可写为

$w=\int{\left\{ \frac{1}{2}\left[ -\varepsilon \left( {{\sigma }_{k}}-{{p}_{\text{w}}} \right) \right] \right\}}\text{d}s+\int\limits_{L}{{{u}_{l}}\cdot {{f}_{l}}}\text{d}l+\sum\limits_{i}{{{u}_{i}}}\cdot {{f}_{i}}$

(5)

式中：

${{\sigma }_{k}}$—应力张量；

${{\varepsilon }_{k}}$—应变张量；

${{p}_{\text{w}}}$—地层孔隙水压力；

${{f}_{l}}$—边界线应力；

${{f}_{i}}$和${{u}_{i}}$—节点位移和节点作用力；

k—下标，k=x，y。

2.2 裂缝扩展正向演化数值模型

为了进行裂缝扩展分析，在有限元方法的基础上，引入储层裂缝扩展有限元^[35-38]。扩展有限元方法与非扩展有限元方法的区别在于，扩展有限元方法在裂缝尖端增加了扩展自由度变量^[35-36]，从而将传统有限元模拟方法中的位移u改写为

$\begin{align} & u={{N}_{i}}(x)[{{u}_{i}}+H(x){{a}_{i}}+{{F}_{j}}(x)b_{i}^{j}], \\ & i=1,2,3,\cdots ,N,j=1,2,3,4 \\ \end{align}$

(6)

式中：

${{N}_{i}}(x)$—一般的节点形函数；

${{u}_{i}}$—连续介质的节点位移量；

${{a}_{i}}$—引入的节点扩展自由度变量；

$b_{i}^{j}$—节点扩展自由度变量；

${{F}_{j}}(x)$—裂缝尖端非对称奇性函数，用于展开裂缝尖端奇性场；

$H(x)$—裂缝的界面阶跃函数，一般选用阶梯函数；

$a$—裂缝尖端极坐标下的扩展方位角。

位移公式（6）中，右端第1 项对模型中全部节点均有效，第2 项仅对其形函数被裂缝切断的节点有效，第3 项则仅对其形函数被裂缝尖端切断的节点有效。裂缝扩展正向演化数值网格如图 1。

在非各向同性的储层岩石介质中，裂缝尖端奇性场依赖于裂缝尖端的具体位置，因此要准确地模拟裂缝尖端奇性场，需要不断追踪裂缝的扩展轨迹，降低模拟速度。一种替代方案是基于裂缝扩展沿断裂径迹扩展行为的扩展有限元模拟方法。这种模拟方法可以广泛应用于脆性和韧性裂缝的扩展模拟，且不需要提前指定裂缝的断裂轨迹。由于这种模拟方法不再需要跟踪处理裂缝尖端奇性场，而只保留了裂缝阶跃函数，因此可以提高模拟速度。为了考察裂缝自由扩展路径和尖端奇性，本文并不提前制定裂缝的断裂轨迹。

2.3 多场趋势克里金系统模型

基于地质统计学的裂缝建模模拟方法是实际裂缝建模中应用较为广泛的方法之一^{[1-9, 13-16]}，但现有的地质统计建模方法^[39]，包括普通克里金（OK）、指示克里金（IK）和趋势克里金（KT）等，没有考虑裂缝尖端断裂扩展动力学因素影响。为了考察裂缝尖端扩展动力学因素对裂缝属性随机分布的影响，在趋势克里金估计系统的基础上，建立多场趋势克里金系统模型，该模型将裂缝特征属性Z(X) 的分布估计模型写作两项之和

$Z(X)=m(X)+Y(X)$

(7)

式中：

$m(X)$—确定部分，也就是根据裂缝反演系统获得的裂缝分布趋势；

$Y(X)$—具有空间相关性的随机变量。

含场趋势的克里金估计方法与泛克里金估计方法很相似，但是需要将泛克里金估计中的确定部分$m(X)$改写为1组辅助模拟场量$\left( \left\{ {{\Psi }_{0}}\left( X \right),{{\Psi }_{1}}\left( X \right),\cdots ,{{\Psi }_{p}}\left( X \right) \right\} \right)$的函数形式

$m(X)={{c}_{0}}{{\Psi }_{0}}(X)+{{c}_{1}}{{\Psi }_{1}}(X)+\cdots +{{c}_{p}}{{\Psi }_{p}}(X)$

(8)

式中：

${{\Psi }_{0}}\left( X \right)=1$；

${c_i}$—待定系数，$i=0,\cdots ,p$；

p—辅助模拟场量的数目。

根据克里金估计理论，辅助场量$\left( \left\{ {{\Psi }_{0}}\left( X \right),{{\Psi }_{1}}\left( X \right),\cdots ,{{\Psi }_{p}}\left( X \right) \right\} \right)$的无偏估计需满足如下方程^[39]

$\begin{align} & \sum\limits_{i=1}^{n}{{{\lambda }_{i}}{{\Psi }_{j}}({{X}_{i}})}={{\Psi }_{j}}({{X}_{0}}) \\ & \forall j=0\cdots p \\ \end{align}$

(9)

在保证无偏性的前提下，若令估计方差最小，可以使用经典的拉格朗日方法求解。此时，克里金估计系统可以写为如下矩阵形式

$\left[ \begin{matrix} \Sigma & \text{ }F \\ {{F}^{T}} & 0 \\ \end{matrix} \right]\left[ \begin{matrix} \lambda \\ \mu \\ \end{matrix} \right]=\left[ \begin{matrix} {{\sigma }_{0}} \\ {{f}_{0}} \\ \end{matrix} \right]$

(10)

式中：

$\Sigma -Z$的协方差矩阵；

F—观测点处的辅助场量趋势；

$\lambda $—克里金权重矢量；

$\mu $—拉格朗日乘子矢量；

${{\sigma }_{0}}$—待估点和观测点处的属性协方差矢量；

${{f}_{0}}$—待估点处的场量模拟值。

裂缝特征属性估计值可以写为

$Z_{0}^{*}={{\lambda }^{\text{T}}}Z=[{{Z}^{\text{T}}}\text{ }0]\left[ \begin{matrix} \lambda \\ \mu \\ \end{matrix} \right]$

(11)

联立式(10)和式(11)，可以得出裂缝特征属性估计值

$\begin{array}{*{35}{l}} Z_{0}^{*}=[{{Z}^{\text{T}}}0]{{\left[ \begin{matrix} \Sigma & F \\ {{F}^{\text{T}}} & 0 \\ \end{matrix} \right]}^{-1}}\left[ \begin{matrix} {{\sigma }_{0}} \\ {{f}_{0}} \\ \end{matrix} \right] \\ =[{{b}^{\text{T}}}{{c}^{\text{T}}}]\left[ \begin{matrix} {{\sigma }_{0}} \\ {{f}_{0}} \\ \end{matrix} \right] \\ \end{array}$

(12)

式中：

$\begin{align} & {{b}^{\text{T}}}=[1-F{{({{F}^{\text{T}}}{{\Sigma }^{-1}}F)}^{-1}}{{F}^{\text{T}}}{{\Sigma }^{-1}}]{{\Sigma }^{-1}}{{Z}^{\text{T}}}; \\ & {{c}^{\text{T}}}=F{{({{F}^{\text{T}}}{{\Sigma }^{-1}}F)}^{-1}}{{\Sigma }^{-1}}{{Z}^{\text{T}}}。 \\ \end{align}$

式(12)显示，$Z_{0}^{*}$可以表示为随机项${{b}^{\text{T}}}{{\sigma }_{0}}$和确定项${{c}^{\text{T}}}{{f}_{0}}$之和。

3 数值模拟 3.1 研究区域

大情字井地区地处松辽盆地中央拗陷长垣长岭凹陷构造带内，受孙吴双辽基底断裂活动影响较强烈，是中央拗陷构造带内构造变动相对发育地区^[40-41]。早白垩世末，由于研究区经历了反时针向南北直扭的构造变动，形成了区域性的、并受基底断裂活动所控制的裂缝体系^[41]。早新生代晚期，即中新世，由于研究区又经历了近东西向的挤压构造变动^[42]，形成了新的区域性的、并与断层活动伴生的局部性的裂缝体系。这两个构造期所形成的裂缝体系的叠加或复合，构成了研究区现今所展现的裂缝网络系统。本文研究对象为大情字井地区青山口扶余致密油层（图 2）。从目前该区完钻井钻井液漏失情况统计可以看出，受裂缝影响该区探井85% 以上都有漏失现象。

3.2 断层扩展模拟结果

裂缝反演共包含3 个步骤：（1）收集与裂缝有关的统计数据，包括区域构造演变历史、构造解释数据等，另外还包括钻井取芯以及各种成像测井数据等。（2）建立属性模型，基于地质统计学建立属性模型，为裂缝模拟提供基础模型，要求基础模型要大于目标研究区的尺度，以便减小边界效应对研究区裂缝建模的影响。由于储层岩石高度非均质，其应力强度因子的解析形式不易获得，通常只能根据大量不同孔隙度和几何形状的裂缝进行裂缝扩展实验。（3）根据属性模型，基于扩展有限元方法，进行主断层扩展历史拟合，拟合采用退火方式。

进行退火模拟时需要不断修正边界加载条件和岩层断裂参量，因此会产生众多的中间反演模型，限于篇幅，这里仅给出最终的裂缝反演模型，如图 3。图 3a 为初始模拟步有效应力分布图；图 3b 为第50 步裂缝扩展图和有效应力分布图；图 3c 为第100 步裂缝扩展图和有效应力分布图；图 3d 终止步裂缝扩展图和有效应力分布图。

3.3 裂缝断裂准则

实验表明，剪性裂缝形成于压应力条件，根据Mohr-Coulomb 失效原理，对于局部各向同性的低孔隙度岩石，在三轴应力（主应力，σ₁ > σ₂ > σ₃）条件下，岩石的破裂面平行于σ₂，且与σ₁ 的夹角为θ_i = 45 − φ_i/2，其中φ_i 为岩石内摩擦角；张性裂缝形成于拉应力条件下，裂缝面平行于σ1，而垂直于σ₃。在剪性裂缝和张性裂缝之间的过渡带，还存在着混合裂缝^[43]。混合裂缝的裂缝面方向可以由Griffith 的裂缝扩展理论给出。将Griffith 理论和Mohr-Coulomb 原理结合，可以得到一个混合失效判据^[43] （图 4）。

应用混合失效判据可以判断从剪性裂缝到张性裂缝区间内任意裂缝的开裂方向。首先，基于应力场模拟结果，结合趋势克里金方法，得到裂缝发育的估测位置。然后，根据应力模拟得到的主应力，计算裂缝发育位置处的Mohr 圆。之后，根据Mohr圆与混合失效判据线的相对关系判断裂缝的有效性：当Mohr 圆小于混合失效判据线时，则未达到破裂条件，裂缝无效；当Mohr 圆与混合失效判据线相切时，裂缝达到破裂条件，为有效裂缝。最后，计算切点—Mohr圆心连线的倾角2θ_i（图 4），进而确定破裂面与σ₁ 的夹角θ_i。

3.4 多场趋势克里金系统模拟结果

利用成像测井解释数据，获得了研究区及周边地区5 口井的裂缝倾角和倾向数据，如表 1。作为多场趋势克里金模拟的基本输入数据，利用裂缝的成像测井解释数据可以对裂缝进行分组，然后对每一组裂缝单独进行建模^[44-46]，因此，模拟产生的裂缝分组数通常等于或者大于成像测井分析得到的裂缝分组数。

图 5a 为基于随机扩展有限元裂缝建模方法获得的大情字井地区黑43 块青山口组的裂缝分布图。

基于裂缝反演系统获得的主应力场（σ₁，σ₂，σ₃）和主应变场模拟（ε₁，ε₂，ε₃）结果，以及趋势克里金方法得到的裂缝的走向向量r和法向向量n，可以计算出每个裂缝所受的剪切应力（τ = r^TDiag(σ₁,σ₂,σ₃)n）和法向应力（σ_n = n^TDiag(σ₁,σ₂,σ₃)n），从而获得每个裂缝的开闭指数：Ω = σ₁ − σ_n/(σ₁ − σ₃)，滑动指数：S = τ/σ_n。图 5b 为裂缝开闭指数分布图，图 5c 为裂缝滑动指数分布图。图 5b、图 5c 中的开闭指数和滑动指数均进行了归一化处理，归一化后的开闭指数${\tilde{\Omega }}$= (Ω − min(Ω))/(max(Ω) − min(Ω))，归一化后的滑动指数${\tilde{S}}$= (S − min(S ))/(max(S ) − min(S ))。

4 讨论

模拟结果显示，裂缝正演和反演系统可以应用于裂缝的属性建模，且可以给出相对严格的裂缝尖端扩展结果。多场克里金方法不仅能够充分利用裂缝反演系统的模拟成果，还添加了岩芯、测井等约束数据，得到的裂缝分布模拟结果与成像测井解释获得的裂缝属性数据符合良好。裂缝反演系统的裂缝扩展方向和裂缝轨迹结果与真实裂缝轨迹总是存在一定的误差，最终得到单个井点的裂缝分布也与成像测井解释获得的裂缝属性分布存在一定的误差。裂缝扩展模拟和多场趋势克里金模拟结果的误差说明了储层裂缝具有复杂的空间结构。如果能够找到正确的模拟参数，减小边界条件和属性的不确定性，必然能够大幅度减小模拟误差。

本文模型考虑了模型参数的不确定性，因而引入了模拟退火方法，每个模型均经过大量迭代运算，然而裂缝的反演结果仍然存在误差。这主要有3 个方面的原因：（1）反演系统的目标搜索空间维度很多，使裂缝反演出现多解性；（2）裂缝反演系统是个非线性系统，某一步微小的模拟误差，在传播多步以后，被非线性放大；（3）储层岩石的弹性模量、泊松比和断裂力学参量，存在小幅度空间涨落，导致模拟出现误差。

5 结论

（1）建立了随机扩展有限元裂缝建模方法，该方法包括裂缝扩展的正反演系统和多场克里金模拟系统。裂缝扩展正演系统采用扩展有限元分析方法，裂缝扩展反演系统采用模拟退火方法自动寻找全局最优解。

（2）建立的裂缝扩展正演系统在不设定裂缝开裂轨迹的情况下，模拟剪切裂缝、张性裂缝和混合裂缝的尖端扩展历史；裂缝扩展反演系统则可以自动调整应力边界条件，得到的裂缝模拟结果与真实裂缝符合良好。利用多场克里金模拟系统，综合成像测井裂缝解释资料和裂缝扩展正反演系统的模拟结果，最终建立了大情字井地区上白垩统青山口组的裂缝模型。

（3）基于随机扩展有限元的裂缝建模方法综合了随机建模方法和几何动力学裂缝建模方法的优点，对于裂缝扩展历史复杂，成像测井裂缝解释资料不足的地区，可自动拟合出区域裂缝扩展历史，根据断裂判据，给出裂缝分布趋势；对于成像测井解释资料丰富的地区，在裂缝扩展自动拟合的基础上，采用多场克里金估计，可以给出统计最优的裂缝分布模型。