裂缝型非均匀介质地震波正演模拟研究

2019-11-19杨小江赵红娟赵伟超

物探化探计算技术 2019年5期

侯凯，李黎，罗泽，杨小江，赵红娟，赵伟超

(中海石油(中国)有限公司深圳分公司研究院，深圳南山 518000)

0 引言

裂缝性油气藏作为一种非常复杂的隐蔽油气藏，广泛发育在礁灰岩、泥页岩以及火成岩中。据统计，世界上已探明的总储量的一半为裂缝性油气藏[1]。裂缝性油气藏研究最突出的问题是裂缝的尺度非常小，远远低于其他地质体的地震响应，地震采集的数据中虽然包含裂缝型储层的响应信息，但是裂缝能量的响应很弱，识别具有很大的困难。在20世纪80年代以前，裂缝性储层的研究依据野外露头、岩芯和测井资料等，从构造地质学角度定性研究裂缝的类型、成因机制和分布特征。后来根据应力场数值模拟和岩石破裂准则，进行半定量研究。随着计算机技术的快速发展，裂缝研究逐渐定量化[2]。Wu Huang等[3]利用De Wolf近似和相位屏算子计算三维背景下的背向散射场，符力耘等[4]建立了均匀背景介质的Green函数基本解体积分方程，尹军杰等[11]对于大量的地质模型进行散射波数值模拟，系统分析了散射波与反射波的差异，蕾蕾[5]根据散射积分算法，对于缝洞模型进行了可变网格的积分模拟，分析了散射波长的传播特征。裂缝储层的正演模拟虽然取得了很大地进步，但是仍然无法精确刻画地下厘米和毫米级别的裂缝响应特征。

本次研究从散射理论出发，采用多级可变网格的正演模拟方法对于裂缝进行了大量的模拟，实现了厘米甚至毫米级别裂缝的地震正演模拟，并统计和定量分析了裂缝群的发育密度、速度、倾角和分布范围等不同时，地震波场的传播特征，准确地找出裂缝的各项参数与地震接收记录的波场对应关系，大大提高了裂缝模拟的计算效率，为裂缝型油气藏的准确预测提供了有力根据。

1 地震散射理论

散射理论最早是在量子力学中提出的，定义为由传播介质的不均匀性引起的光线向四周射去的现象。假设介质是完全均匀的，地震波从震源出发将按照球面扩散的形式以不变的速度传播。若介质中存在与均匀背景不一样的非均质体时，地震波的传播发生改变，此时的非均质体可以看作二次震源，我们称这种现象为散射现象[6-8]。根据散射理论，我们将介质分为背景介质和扰动介质，在背景介质中传播的波场，称之为背景波场。如果介质的性质发生变化，那么这部分介质称为扰动介质。由介质的扰动部分产生的波场称为散射波场。

地震波在均匀介质中声波的传播符合以下规律：

(1)

其中：C(r)为均匀介质中的传播速度；u0(r,t)为产生的波场，将式(1)转换到频率域得到式(2)。

(2)

在索莫菲辐射边界条件的假设下：

r→∞,r=|r|

(3)

ε(r)为相对于C0(r)的扰动波场，并满足C0(r)横向不变，C0(r)可以表示为C0(z)。

(4)

根据扰动理论，地震波场u可以分解为均匀介质中的波场u0与扰动介质产生的波场的叠加，即：

u(x,y,z;ω)=u0(x,y,z;ω)+U(x,y,z;ω)

(5)

均匀介质中的波场可以表示成式(6)。

(6)

将式(5)、式(6)带入到式(2)中得式(7)。

u(x,y,z;ω)

(7)

用格林函数法求解求解方程(7)，最终得到波场的积分表达行式如下：

u(x,y,z;ω)=u0(x,y,z;ω)+us(x,y,z;ω)=

u0(x,y,z;ω)+

(8)

方程(8)就是Lippmann-Swinger方程[9]，此方程是一个非线性方程，在实际中需要使用born公式简化。

2 多级可变网格

常规的有限差分正演模拟方法考虑到模拟的精度，又要考虑计算的效率，空间采样间隔通常是5 m～10 m。Moczo[16]首次提出了可变网格的思想，根据非均匀体的尺度，对空间步长调整，可以在一定程度上，保证模拟效率的情况下提高精度。然而对裂缝储层来说，仍然不能准确刻画，为了克服可变网格倍数的限制，笔者首次提出并采用了多级变网格的方法对于裂缝进行模拟，成功的将变网格区域继续放大。

如图1所示，常规网格无法对于小尺度模型准确刻画，经过可变网格加密后，可以提高大尺度溶洞的刻画程度，但仍然无法准确地模拟裂缝。在这种情况下，在可变网格的基础上，对于微裂缝继续空间加密，经过两次的加密可以达到了模型放大100倍的效果，准确地刻画微裂缝，这种局部加密的方法提高了模拟的精度，又保证了正演的效率。

图1 多级变网格的模型剖分示意图Fig.1 The schematic diagram of multistage variable mesh model

图2 过渡区域空间导数网格计算点Fig.2 Spatial derivative grid computing point of transition region

图3 含裂缝的断层模型Fig.3 A fractured fault model

1)多级可变网格空间过渡。为了保证变网格区域的波场稳定过渡，需要重新定义过渡区域网格差分格式，让加密区域的波场得到耦合。如图2所示，变网格区域参与导数计算点是不对称的，为了解决此问题，小网格区域取点不是等距的。对差分格式和差分系数做如下修正：

u(x-mNdx,z)-2u(x,z) ]

(9)

u(x,z-mNdz)-2u(x,z) ]

(10)

式(9)和式(10)分别给出了加密区域x方向和z方向的网格空间二阶导数M阶精度的计算公式。dx表示加密网格步长，m表示空间导数的点到参与计算点的距离。

图4 含裂缝断层模型的接收记录Fig.4 The receiving record of the fractured model(a)单程波的接收记录及裂缝响应；(b)双程波的接收记录及裂缝响应

2)多级变网格时间过渡方法。在使用多级可变网格时，模型的时间采样步长同样也发生了变化。波动方程的空间采用间隔与时间采样间隔必须满足稳定性条件，才能保证算法是稳定性的。当空间采样间距减小时，时间采样的间距必须相应的减小。多级变网格的时间采样步长与空间采样步长倍数的变化保持一致。设变网格的倍数是M时，t时刻到t+1时刻之间的波场值用线性插值算法计算M-1时刻的波场值，在小网格步长内计算t时刻到t+MΔt的波场(Δt是小网格的时间步长)，如果大网格的时间步长为ΔT，则Δt=1/M(ΔT)。最后将t+MΔt时刻的波场值回归到大网格的网格上。若在一阶变网格区域内需要进行二级变网格时，需要对于二级变网格的时间步长和空间步长做与同样方法的改变即可。经过测试可知，改变时间采样步长不会引起明显的人为反射，因此不做特殊的处理。

2.1 积分法与有限差分微分法的波场特征对比

图3为一个断层模型，在断嵌位置处发育垂直裂缝群。模型设置的参数如下：裂缝的宽度为0.5 m，断开深度为50 m，模型大小设置为600×500，间隔为5 m。激发主频30 Hz的雷克子波，时间间隔是0.5 ms。震源位于地面的中心位置。可变网格积分法使用高斯衰减边界条件，有限差分微分法使用PML边界条件。

图4为两种裂缝模拟方法的接受记录。从接受记录情况看，波动方程微分法可以接收到丰富的波场信息[10]，同时也存在裂缝响应弱，容易被其他信息掩盖的缺点。波动方程积分法得到的地震记录不含直达波、多次反射波等信息，如图4所示，相对于双程波的模拟来说，单程波模拟的裂缝响应更清晰，且噪音成分少，可以直观地找到裂缝非均质体的响应特征。综合分析，对于含裂缝的储层使用单程波进行裂缝响应的研究，更能观测到裂缝的响应，满足定量分析的要求，同时也能满足微裂缝模拟对计算机的要求[11-13]。

图5 模型示意图Fig.5 A fractured model diagram(a)常规网格;(b)一级变网格;(c)二级变网格

图6 地面记录的接收结果Fig.6 Seismic record(a)常规网格的模拟结果；(b)多极变网格的模拟区域

图7 接收记录的波形分析Fig.7 Wave figure from cmp

图8 背斜模型Fig.8 The anticline model(a)无裂缝模型;(b)含裂缝模型

2.2 多级可变网格积分法与常规网格积分法对比

采用同一模型，分别使用常规网格积分和多级可变网格的方法进行模拟(图5)。震源子波是雷克子波，主频为40 Hz，采样间距是1 ms。为了测试多级可变网格，我们将模型中的第二个矩形溶洞进行二级变网格的加密，第一次选取模型中的虚形框部分。横向加密10倍得到一级变网格的模型，变网格区域内的空间采样间隔变为1 m，时间采样间距为0.1 ms。然后对虚形框部分继续加密10倍，得到二级变网格模型，这时候的网格间距变为0.1 m，时间采样间距为0.01 ms。

对比图6(a)、图6(b)可看出，多级变网格的模拟区域的记录与常规网格的记录一致，多级变网格没有带来明显的边界反射。抽取如图7所示位置出的波形图做对比分析，可以清晰地看出,常规网格的模拟方法与多级变网格的模拟方法得到波形基本一致。可以看到多级可变网格可以实现高倍数的变网格模拟，而且稳定性和精度很高。

2.3 网格积分法与网格微分法计算效率对比

图8给出的背斜模型，模型的速度在图中已经标出。在图(b)中发育一组垂直的裂缝群，裂缝的速度为1 500 m/s。模型的网格大小为120×100，空间采样间距为10 m，震源使用雷克子波，主频为30 Hz，时间采样步长是0.5 ms，采样点是1 500 个。对这两个模型分别使用多级变网格积分法和多级变网格微分法正演模拟，两种方法给定的模型参数以及震源参数是完全相同的。

图9 地面接收记录Fig.9 Seismic record(a)单程波无裂缝模拟结果;(b)单程波有裂缝模拟结果;(c)裂缝响应;(d)双程波无裂缝模拟结果;(e)双程波有裂缝模拟结果;(f)裂缝响应

由图9可以看出，双程波模拟结果中有直达波的成分，而单程波模拟仅含有一次反射波。两种方法都需要进行高倍数的变网格来模拟裂缝。多级可变网格采用的10倍的加密，采用微分法正演模拟的计算时间为6 h。采用格林函数积分法的单程波模拟，进行正演模拟的计算时间为20 min。因此可以得出，对于含有裂缝非均质体的模型，格林函数积分法要比多级变网格的双程波模拟效率高20倍左右。

3 微裂缝地震响应定量分析

3.1 微裂缝响应特征

通过图10中给定2个模型分析裂缝的接收记录响应特征，给定裂缝的宽度为5 mm，倾角为60°，裂缝间的横向间距为0.1 m,图10(a)中发育3组裂缝，裂缝的切深为10 m，范围是35 m。图10(b)中不发育裂缝。

通过对比图10可以看出，裂缝在地层中的能量很弱，裂缝群能量会对界面反射有干涉作用，在地震记录中以尾波的形式出现。

3.2 震源主频度裂缝能量影响

图11模型范围是3 000 m*3 000 m，空间采样间隔为5 m，时间采样率是1 ms，倾斜裂缝发育的位置是1 000 m，裂缝横线的断开大小只有5 mm，裂缝速度为2 000 m/s，裂缝间的距离是0.25 m，角度为60°，按照表1改变震源的主频，定量研究了裂缝的响应特征。

通过抽取地震接受记录，对地震波形和能量变化定量分析(图12)。由图12可以看出，裂缝对于不同主频的震源响应不同，主频越低地震波形越分散。裂缝能量并不是随着主频的升高而升高，而是呈现抛物线的变化特点，当裂缝角度为60°时，主频为50 Hz的裂缝能量最强，然后随着主频的升高而能量降低。

3.3 裂缝发育密度与能量响应的定量分析

裂缝发育密度越大，裂缝响应的能量越强，裂缝发育密度与能量间的响应存在一定的规律，通过设定模型参数来定量分析。模型的大小为2 500 m*2 500 m，空间采样间隔5 m，时间采样间隔为1 ms，给定裂缝的尺度为0.005 m，裂缝速度为2 000 m/s，纵向切深为10 m，裂缝的倾角为60°。震源主频为50 Hz。按照表2改变裂缝的发育密度，获得裂缝的正演地震记录。

图10 微裂缝地震响应Fig.10 The seismic record of micro-fractuew (a)模型1及对应的炮记录；(b)模型2及对应的炮记录；(c)中心道炮记录波形图

图11 倾斜模型及对应炮记录Fig.11 Tilt model and corresponding seismic record

20304050607080

表2 裂缝间距(单位：m)

图14为正演地震记录，分析得到裂缝密度越大，裂缝的响应也越强。图15为裂缝能量与裂缝密度的统计关系，可以看出裂缝发育密度与裂缝能量响应之间并不是单纯的线性关系，而是呈现一定的对数关系，通过分析裂缝的能量，反过来可以得到裂缝的发育密度，对于裂缝地识别和定量描述有着十分重要的意义。

3.4 实际模型测试

经过理论模型的验证，该算法可以快速模型微裂缝的地震特征。图16为某工区低幅度构造的速度模型，模型大小为2 000 m*2 500 m，在模型深度1 850 m～1 870 m处发育有6组裂缝群，裂缝断开大小为0.005 m，裂缝速度2 000 m/s，裂缝间的横向距离是0.1 m，裂缝断开深度为10 m。震源为30 Hz的雷克子波，分别在模型的左边、中间和右边放炮，图16(b)地震接收到的信号。

图17对比了有裂缝和无裂缝发育的地震接收信差别，可以看出，由于裂缝的干涉作用地层界面的反射能量明显地增强了。同时也验证了该算法对于毫米级别的裂缝模拟有着很好的适用性。

4 结论

地震波正演模拟在地震资料采集、处理和解释过程中具有十分重要的作用，地震波正演的热点逐渐集中于对复杂微小构造的适用性和模拟效率的高效性[14-15]。对于含有裂缝的非常规油气藏，由于裂缝的发育尺度只有厘米和毫米级别，要实现裂缝的准确模拟，常规的地震正演模拟方法存在模拟效率和精度低、对计算机要求高等不足。