花岗岩裂隙储层各向异性高阶有限差分数值模拟

2012-10-10龚育龄

东华理工大学学报(自然科学版) 2012年4期

张华，龚育龄

(东华理工大学放射性地质与勘探技术国防重点学科实验室，江西抚州 344000)

高放射性核废物的安全处置是一个世界性难题，也是公众所关心的安全和环保问题。目前提出的方案是深地质处置，即在距地表深约500～1 000 m的地质体中建造“地质处置库”，地质处置库围岩之一为花岗岩岩体(郭永海等，2001;王驹等，2000)，然而在花岗岩这类硬岩中建设高放废物地质处置库，裂隙发育带的存在将是决定未来处置库成败的关键因素之一，因为在高放废物地质处置库中，裂隙带是放射性核素迁移的主要通道。随着时间的推移，核素随着地下水流在节理裂隙网络中的流动，极有可能迁移到生物圈，对人类环境造成极大的破坏;同时裂隙带发育程度直接指示着岩体的破碎程度，也在一定程度上预示着区域应力场分布及其方向(金远新等，2007;熊章强，2001;余运祥等，1994)，因此，基于处置库安全性与稳定性的考虑，必须查明处置库岩体的裂隙带发育情况，而地震勘探方法是有效预测和确定裂隙带分布的地球物理勘探方法(程纪星等，2002)。

然而由于各种条件的限制，地震勘探方法在花岗岩型地质处置库选址中的应用较少，因此，本文主要利用弹性波交错网格高阶有限差分算法，对花岗岩中裂隙带为各向同性介质模型和VTI型各向异性介质模型的弹性波传播特征进行数值模拟，获得地震波场快照和单炮记录，并分析其地震响应特征，为地震勘探在高放废物地质处置库选址中的应用提供理论基础和技术指导。

1 各向异性弹性波方程

二维一阶速度—应力VTI弹性波动方程组为:

其中，cij为弹性系数。

对于各向异性的表征，采用Thomsen(1986)提出弱各向异性介质理论。Thomsen各向异性参数有ε(纵波各向异性参数)，γ(横波各向异性参数)，δ(变异系数)。根据Thomsen原理，在VTI介质中，这三个各向异性参数与弹性系数转换关系为:

2 高阶有限差分算法

设函数u(x)连续，且具有2N+1阶导数，则交错网格一阶导数2N阶精度差分近似式表示为:

其中，x是实变量，x0是一实数，Δx为空间步长，m=1，2，…，N;cm为差分系数。据(裴正林，2006)的推导，当

在本文的数值模拟中采样交错网格高阶有限差分进行求解方程(1)。

3 边界条件

由于只能对有限区域进行数值模拟，必然会产生边界反射的影响。为消除这种边界，本次模拟根据完全匹配层吸收边界条件(PML)原理，采用变量分裂方法，对方程(1)中的一阶速度—应力VTI弹性波动方程进行分裂:

式中上标和代表该项只与相应的空间导数有关。PML吸收边界方程组为:

式中d(x)和d(z)分别是与x和z方向有关的衰减因子 (Bérenger，1994;钱进，2010)。

4 数值模拟与分析

4.1 地质模型

花岗岩为低孔隙率(孔隙度一般为0.3% ～0.7%)致密岩石，但是花岗岩体局部裂隙带构造发育，其厚度不均，它主要由于不同时期的构造应力或者地下水活动造成，裂隙之间被花岗岩脉、石英脉、方解石、伊利石等充填，同时也成为水储存、滞留形成地下水的空间。如果花岗岩裂隙带中裂隙不特别发育时，可以等效为各向同性介质，但当内部裂隙水平方向较为发育时，可以等效为VTI型各向异性介质。表1给出了VTI型裂隙带与花岗岩物性参数。模型的网格剖分采用1 000×1 000 m，网格间距为2 m，时间采样间隔=0.5 ms，记录长度为0.5 s，裂隙带埋深为500 m，裂隙带厚度分为5 m厚和20 m厚两种情况，顶底板围岩都为花岗岩，震源为主频为60 Hz的雷克子波，震源位于(500，40 m)位置。

表1 地质模型物性参数表Table1 Petrophysical parameters of geological model

4.2 三层均为各向同性介质的模型

首先模拟三层均为各向同性介质的波场快照和单炮记录，以便与各向异性层状介质的模拟结果进行对比，应用空间交错网格高阶有限差分算法，得到t=180 ms时的水平分量和垂直分量的波场快照和单炮记录，其中本文所有波场快照纵横波坐标都为记录道数，单位为相对值，单炮记录纵坐标表示时间(s)，横坐标表示记录道数，如图1所示，该模型的波场快照经过界面反射和透射后，明显存在的弹性波有:①反射P波、②反射转换S波、③透射P波和④透射转换S波⑤直达波五种，见图1中详细标注。同时结合单炮记录可以看出，在水平分量上，反射P波能量相对较弱，其振幅小于反射S波和透射P波。在垂直分量上，反射P波能量较强，转换波次之，反射S波投影相对很弱，需要强调的是波场快照中由于时间选取较大，因此没有出现直达波。

同时可以看到底板反射的折射波并没有在波场快照中出现，这说明对于薄层花岗岩裂隙带来说，并不会产生底板反射的折射波。此外，从图中还可以看出由于裂隙带为薄层，只有一个反射界面有强反射产生，所以该层反射波应为裂隙带顶底板反射的复合波，无法将其分离。

图1 各向同性裂隙带单炮记录和波场快照(5 m)Fig.1 Single shot record and wave field snapshot of isotropic fracture zones(5 m)

图2 各向同性裂隙带波场快照和单炮记录(20 m)Fig.2 Single shot record and wave field snapshot of isotropic fracture zone(20 m)

图2为裂隙带厚度为20 m时的波场快照和单炮记录，由图2可以看出，随着裂隙带厚度的增加，在该模型波场快照中存在的弹性波除了与图1相同特征的波场外，且在上下两个反射界面上都有相对较强的反射波(⑥，⑦)和透射波(⑧，⑨)产生，使得裂隙带顶底板反射波可以分离，说明该裂隙带的厚度可以进行预测。此外，反射P波和反射转换S波下方有一些裂隙带顶底板反射的多次波出现(⑩)，其中水平分量的多次波要弱于垂直分量，且这一现象在反射S波上比反射P波更为明显，因此在实际资料处理与解释过程中时应加以注意。

4.3 均匀VTI型裂隙带模型

根据表1中的岩石物性参数分别建立均匀VTI裂隙带模型，模型大小、空间和时间采样间隔、记录长度和震源子波主频以及有限差分精度与上述各向同性介质模型相同，震源位于(500 m，500 m)处。图3为VTI型裂隙带在220 ms两个时刻的波场快照。从图3中可以看出，在均匀VTI型裂隙带中存在着两类波:①拟P波和②拟SV波。拟SV波在对角线方向上出现了分叉现象。而且，VTI型裂隙带在沿垂直对称轴方向传播的速度明显要大于沿水平对称轴方向，即弹性波沿裂隙方向的传播速度要明显大于沿垂直裂隙方向的速度。从而使得运用地震各向异性理论预测裂隙发育方向成为可能，同时也验证了VTI等效介质理论的准确性。并且从图3波场快照可以看出，本文所采用的PML吸收边界条件，能较好地压制人工边界反射。

4.4 中间层为含VTI型裂隙带模型

利用表1中花岗岩裂隙带物性参数建立二维层状模型来观察VTI型裂隙带反射波特征。模型参数与上述各向同性介质模型相同。考虑到算法和边界吸收条件的应用，将震源置于(500 m，100 m)处，接收排列置于深度80 m处进行水平布置。分别对裂隙带厚度为5 m和20 m两种情况进行数值模拟，模拟的单炮记录如图4所示。同样，地震单炮记录也分水平分量和垂直分量。

从图4中可以看出，地震单炮记录中共包含四类波:①直达波，②裂隙带反射纵波(PP)，③裂隙带反射横波(PS);④为多次波。对于薄层裂隙带，PP与PS均为裂隙带顶底板反射的复合波，无法进行分离，见图4a和图4b，当裂隙带厚度为20 m时，厚裂隙带层顶底界面的反射波可以分清，从能量角度来看，垂直分量中PP能量要强于PS，这取决于弹性界面的纵横波反射系数;而水平分量则与之相反，这主要是因为弹性波传播过程中存在一定的模式转换和能量转换。同时，从图4也可以看出，反射波的振幅随偏移距和厚度有较为明显的变化。而且随着裂隙带厚度的增加，裂隙带反射波会在顶底板发生连续反射，从而导致了裂隙层中出现多次波，这一现象在PS波上反映更为明显，见图4c和图4d，而这点与各向同性介质模型类似。