有限方位角条件下压噪叠前纵波各向异性裂缝检测方法

摘要

裂缝介质是各向异性介质；各向异性是指介质的物理性质(如地震波的速度，振幅等)随测量方向改变而变化。叠前地震资料可以用来检测裂缝。裂缝是油气的储集空间和运移通道，对裂缝的检测对寻找油气意义重大。常规裂缝检测方法往往采用宽方位地震资料。这是因为地震资料存在噪音，若应用有限(窄)方位资料，噪音的存在会使反演失效，而宽方位资料能够克服噪音的不利影响。但宽方位采集的成本非常高，可能是有限方位采集成本的好几倍，如何利用有限方位角地震资料进行裂缝检测具有重大经济意义。本发明在常规反演算法的基础上，考虑噪音因子，克服噪音对反演结果的不利影响，能够利用含有噪音的有限方位角地震资料，准确地反演出裂缝密度和走向。

说明书

技术领域

本发明属于石油勘探领域，涉及裂缝储层综合评价中应用叠前有限方位角地震数据对有效裂缝进行检测的方法。通过考虑噪音因子，有效地去除噪音对反演结果的不利影响。

背景技术

裂缝是岩石受力发生破裂、沿破裂面两侧的岩石没有发生明显位移的一种断裂构造现象。裂缝既是油、气、水等流体的储集空间，又是流体的运移通道。有效裂缝是指开启并充填流体的裂缝。有效裂缝是油、气储存的有效空间，为构造作用或成岩作用等形成。对有效裂缝的探测对于识别油气运移通道，寻找有效储集空间具有重要的指导意义。

目前国内比较成熟的裂缝预测技术是利用叠后的三维地震数据体进行裂缝检测，例如相干分析方法等。这些方法目前还存在很多不足之处。例如相干分析方法能够反映小断裂或断裂发育带，岩性突变带，但主要是用于断层边缘的检测，并且检测所有的断裂体系。当地下介质受多期次构造运动或成岩作用时，早期的裂缝可能被富含碳酸钙的流体所充填而成为无效裂缝。而这种无效裂缝在相干切片上仍能显示。相干分析技术很难预测出有效裂缝的发育特征。

应用地质信息对裂缝进行描述往往会受岩心和成像测井(FMI)资料数量的限制。而地震各向异性裂缝预测技术在横向大面积上具有预测意义，这种技术能够检测有效裂缝引起的各向异性特征。

我们将地震各向异性定义为介质物理性质(如地震波的速度、振幅等)的方位差异特征。裂缝介质模型是一种典型的各向异性介质模型。理论证明，当地下介质中含有一系列垂直排列且相互平行的裂缝时，地震波旅行时和振幅会随地震测线方位角的变化而变化。对于这一点许多学者做了大量的研究，如Ruger(1996，1998)、Tsvankin(2001)、Grechka等(1999)。当地震纵波在含有裂缝的各向异性介质中传播时，纵波在各个方位上传播的速度不同，沿着裂缝方向传播速度快，垂直裂缝方向传播速度慢。另外沿着裂缝方向振幅能量衰减小，垂直裂缝方向振幅能量衰减大。地震波在裂缝介质中传播时显示的方位差异特征，可以帮助我们探测裂缝的发育特征，包括裂缝走向和密度等。

世界著名的Ruger纵波反射系数公式是：

$>R(i,\phi )=\frac{1}{2}\frac{\mathrm{\Delta Z}}{\bar{Z}}+\frac{1}{2}\{\frac{\mathrm{\Delta \alpha }}{\bar{\alpha }}-{\left(\frac{2\bar{\beta }}{\bar{\alpha }}\right)}^2\frac{\mathrm{\Delta G}}{\bar{G}}+[\Delta {\delta }^V+2{\left(\frac{2\bar{\beta }}{\bar{\alpha }}\right)}^2\mathrm{\Delta \gamma }\left]{\cos }^2(\phi -{\phi }_s)\right\}{\sin }^{2}i$>

$>+\frac{1}{2}\{\frac{\mathrm{\Delta \alpha }}{\bar{\alpha }}+\Delta {ϵ}^V{\cos }^4(\phi -{\phi }_s)+{\mathrm{\Delta \delta }}^V{\sin }^2(\phi -{\phi }_s){\cos }^2(\phi -{\phi }_s)\}{\sin }^{2}i{\tan }^{2}i+...$>

其中Δε^V、Δδ^V、Δγ是Thomsen参数，α、β分别为纵波速度和横波速度，φ为采集方位角，φ_s为裂缝走向，i为入射角。Z为纵波垂直入射时的波阻抗，G为横波的切向模量。由Ruger公式进一步推导可得：

$>\left(\begin{array}{c}R(i,\phi )=C_1+C_2{\cos }^2\phi {\sin }^{2}i+2C_3\sin \phi \cos \phi {\sin }^{2}i+C_4{\sin }^2\phi {\sin }^{2}i\\ C_2=B^{\mathrm{ani}}{\cos }^2{\phi }_s\\ C_3=B^{\mathrm{ani}}\sin 2{\phi }_s\\ C_4=B^{\mathrm{ani}}{\sin }^2{\phi }_s\end{array}\right)$>

B^ani是各向异性参数，也即我们要反演的裂缝密度；C₁，C₂，C₃，C₄是四个未知量，假设共反射点(CRP)点的覆盖次数为n，，则上式可以简化为：

AX＝B

A是n行4列的矩阵；X是C₁，C₂，C₃，C₄四个未知量构成的列向量；B是不含噪音的地震反射振幅。上式是超定线性方程组，一般算法的裂缝反演问题就是求解上述线性方程组中未知量X的问题。很多学者应用上述算法对许多地区宽方位地震资料进行了裂缝预测，如David Y.Wang等(2006)、Gregg Zelewski等(2009)。

应用以上算法进行裂缝反演的不足之处在于，要应用宽方位采集的地震资料进行反演，才能克服噪音对反演结果的不利影响。因为当使用的资料为有限(窄)方位角时，含有噪音的资料会使反演失效或不准确。但是宽方位地震数据所需要的成本非常高，其成本可能是有限方位角资料的好几倍。

发明内容

针对有限方位角下常规裂缝探测方法的不足，本发明提供一种利用有限方位角资料进行压噪叠前裂缝检测的方法。通过本发明的算法改进，能够利用有限方位角资料较为准确地进行裂缝探测，可以极大地减少因采集为宽方位而造成的采集成本，具有重大的经济意义。

本发明实现上述目的的具体实施方案如下：

步骤1：对地震资料进行分方位处理(包括分方位保幅偏移)，并抽取分方位共反射点(CRP)道集。

步骤2：将分方位的CRP道集进行数据合并，并为下一步做准备；

步骤3：用本发明所采用的算法进行裂缝反演，算出每个CRP点上，每个时间采样点上的裂缝密度和裂缝方位角值；

步骤4：对反演数据进行质量控制，其目的在于观察本发明算法中方程求解的准确程度。质量因子代表观测的地震反射振幅和反演的反射振幅(矩阵A与求解的X的乘积AX)的差值，差值越小，则质量因子越大，表明求解的准确程度越高；

步骤5：输出裂缝反演结果即裂缝密度数据体和方位角数据体；

步骤6：对输出后的数据体提取层位数据，并进行地质意义分析。

本发明所提供的有限方位角条件下压噪裂缝探测方法的基本原理如下：在Ruger给出的裂缝模型介质的纵波反射系数公式基础上进行改进，达到消除噪音的影响的目的。本发明的改进之处是：在Ruger反射系数方程AX＝B基础上，考虑噪音因子N，使要求解的方程转化为：

AX＝B+N

通过奇异值分解的方法求解上方程方法如下：

将矩阵A进行奇异值分解可得：A＝UWV

式中U是n行4列正交矩阵；W是4行4列对角矩阵，对角元素是正值或0；V是4行4列正交矩阵。根据奇异值分解的方法求解A的广义逆，进而解上述方程可得：

X＝(A^TA)^-1A^T(B+N)

 ＝V^TW^-1U^TB+V^TW^-1U^TN

通过上式，可以将地震数据的有效信号项和噪音项分开。具体方式是将N乘上一个很小的数或除以一个很大的数，从而使噪音项趋近于0。

本发明的特色是：可以对含有噪音的有限方位角资料，准确地反演出裂缝的各向异性特征，极大地节约因采集方式为宽方位造成的采集成本。

本发明中裂缝探测的适用范围如下：本发明中探测的裂缝密度，是指有效裂缝引起的各向异性强度。所谓有效裂缝是指开启并充填了流体的裂缝，是油、气、水的有效储集空间。各向异性强度是多组裂缝共同作用的结果：当一定深度范围内有效裂缝走向大致相同时，各向异性程度高；反之，若在一定的深度范围内，裂缝的走向呈现出不同的方向，地震波的各向异性特征会被削弱，甚至会显示出各向同性的特征，这就无法通过反演获得地下裂缝的信息。

附图说明

图1(a)是地震波在裂缝介质中的传播示意图；

图1(b)是纵波方位振幅差异示意图。

图2(a)～图2(d)是对本算法优越性的说明。

图3是裂缝反演流程图。

图4是新疆某区块资料采集方位角分布图。

图5是新疆该区块颗粒灰岩层段统计的裂缝密度和充填特征(有效裂缝地质统计)。

图6是新疆该区块颗粒灰岩层段反演的裂缝密度和方位角叠合图。

图7是新疆该区块颗粒灰岩层段反演结果质量控制因子。

具体实施方式

以下结合附图说明本发明具体实施方式。

图1是纵波方位振幅差异原理图。如图1(a)所示，各向异性裂缝介质模型中含有一系列垂直排列且相互平行的裂缝。这里将测线方位角定义为地震波的激发点和接收点的连线与裂缝走向的夹角。当地震波以不同的测线方位角激发和接收时，会存在随着方位角而变化的振幅差异。图1(b)中，横坐标代表测线方位角，纵坐标代表时间。如图所示，在时间为2.0秒时，纵波振幅随着测线方位角而变化。本算法就是基于对纵波方位振幅差异属性的分析，反演出地下介质裂缝的发育情况。

图2(a)～图2(d)是对本算法优越性的说明。图2(a)左侧是全方位采集示意图，采集时炮点和检波点均匀分布在360范围内；右侧是有限方位采集示意图，采集时炮点和检波点在有限的角度范围内，为局部方位分布。图2(b)为对不含噪音的合成地震记录共反射点(CRP)道集进行反演。分为左、中、右三个图进行说明：左图和右图分别为全方位采集和有限方位采集合成的CRP道集，对这两种道集分别采用不加噪音因子的算法和加过噪音因子的算法进行反演计算，反演结果相同，如中图所示，都能显示从左到右依次变化为弱、强、弱、强的各向异性特征。图2(c)为对全方位采集且资料含噪音的CRP道集进行反演，右侧显示用一般算法和本算法都能够反演出从左到右依次变化为弱、强、弱、强的各向异性特征。图2(d)为对有限方位采集且含有噪音的CRP道集进行反演，分左、中、右图进行说明：中图为用不加噪音因子的算法反演后的结果，椭圆圈内表明各向异性无显示；右图为应用本发明所采用的算法，即加过噪音因子后的算法进行反演，椭圆圈内显示了反演出的从左到右较为均匀的各向异性特征。通过以上说明，可以得到如下结论：本算法的优势在于，可以对含有噪音的有限方位角资料，准确地反演出裂缝的各向异性特征。而这正是实际生产中需要解决的问题，通过本算法，可以极大地节省因资料采集为全方位造成的成本。

图3为本发明的裂缝反演流程图：

步骤1：对地震资料进行分方位处理(包括分方位保幅偏移)，并抽取分方位共反射点(CRP)道集。

步骤2：将分方位的CRP道集进行数据合并，并为下一步做准备；

步骤3：用本发明所采用的算法进行裂缝反演，算出每个CRP点上，每个时间采样点上的裂缝密度和裂缝方位角值；

步骤4：对反演数据进行质量控制，其目的在于观察本发明算法中方程AX＝B+N求解的准确程度。质量因子代表观测的地震反射振幅和反演的反射振幅(矩阵A与求解的X的乘积AX)的差值，差值越小，则质量因子越大，表明求解的准确程度越高；

步骤5：输出裂缝反演结果即裂缝密度数据体和方位角数据体；

步骤6：对输出后的数据体提取层位数据，并进行地质意义分析。

以下以该方法在新疆某区块实际资料上的应用为例子具体说明该方法的实际应用效果。

图4是该地区资料采集时的分方位示意图，如图所示，黑粗线分别为施工方向和垂直施工方向。所分方位角，对称分布在施工方向两侧，按照15度一个方位角，逆时针方向分别为方位角1，2，3，4。方位角范围未达到360°，因此属于局部方位覆盖。

图5为通过岩心资料和成像测井(FMI)资料对该层段的裂缝条数和充填特征分析(有效裂缝地质统计)。当FMI识别出的充填特征弱时(FMI显示的电阻率低)，裂缝为有效裂缝，各向异性特征明显；另外当一定深度范围内FMI识别出的裂缝条数多时，表示该深度范围内裂缝密度高，各向异性特征明显。分析该区域的裂缝条数及充填特征，可以得到如下结论：

(1)由I井向M井裂缝线密度逐渐增加。I井FMI识别出3条裂缝，K井识别出4条，M井识别29条，裂缝线密度由0.06条/m增至0.07条/m和0.17条/m。岩心裂缝线密度受取心井段和取心长度的影响，对裂缝密度变化规律表现不明显，但整体上与FMI反映的裂缝变化规律趋于一致。

(2)有效裂缝具有由东南向北西方向逐渐增加的趋势。I井FMI识别的有效裂缝占总裂缝的43％，K占78％，M井占86％。

图6是该地区颗粒灰岩段反演的裂缝密度和方位角叠合图。提取该层段对应时窗上裂缝密度的均方根属性，提取该层段对应时窗顶部的方位角值进行叠合。叠合时采用MATLAB绘图，箭头方向代表裂缝方位角(即裂缝走向)，其长度与裂缝密度相对应。同时底色深浅代表裂缝密度大小。白色线条为该地区的断层平面图。通过FMI(成像测井)资料统计这一层段上的有效裂缝条数和走向，并用玫瑰图的形式显示。玫瑰图的轴向角度指示FMI统计的裂缝方位角，其径向长度表示在该层段该方位上的裂缝条数。

各向异性特征是多组裂缝共同作用的结果，当裂缝走向大致相同时，各向异性强；反之，对于多组不同走向的裂缝，各向异性效应弱。为验证反演结果的正确性，将反演结果与井处的FMI资料进行对比。将反演结果中D井附近裂缝密度和走向叠合图进行放大，并与玫瑰图进行对比，如左下角小图所示：叠合图上箭头方向与玫瑰图上裂缝走向一致，而且叠合图上D井裂缝密度较高，说明反演正确。K井和J井玫瑰图上显示裂缝条数较多，但走向呈不同的方向，这可能造成各向异性效应变弱，而叠合图上K井和J井反演的裂缝密度值较低，这说明反演的结果反应的是多组有效裂缝造成的各向异性特征。

通过对比图6和图5，图6各向异性反演的裂缝密度值G井较低，C井和D井裂缝密度值较高，说明反演的裂缝密度具有从东南向北西方向逐渐增加的趋势，与图5的有效裂缝地质统计规律一致。进一步验证了反演的结果与地质特征的一致性。

图7是对本算法进行的质量控制，将图6和图7进行对比，反演裂缝密度高的地方质量控制因子呈高值，表明通过本发明算法求解方程准确。

以上具体实施方式仅用于说明本发明，而非用于限定本发明。

参考文献

王永刚，乐友喜，张军华.地震属性分析技术[M].东营：中国石油大学出版社，2006.

Tsvankin，I.2001.Seismic Signatures and Analysis of Reflection Data in Anisotropic Media.New York：PergamonPress.

Ruger，A.1996.Reflection Coefficients and Azimuthal AVO Analysis in Anisotropic Media.PhD thesis，ColoradoSchool of Mines，Leadville，Colorado.

Ruger，A.Variation of P-wave reflectivity with offset and azimuth in anisotropic media[J].Geophysics，1998，63(3)：935-947.

Grechka，V.Tsvankin，I.and Cohen，J.K.1999.Generalized Dix Equation and Analytic Treatment ofNormal-Moveout Velocity for Anisotropic Media.Geophysical Prospection47(2)：117-148.DOI：10.1046/j.1365-2478.1999.00120.x.

David Y.Wang.et al.Fracture Analysis for Carbonate Reservoirs Using 3D Seismic P-Wave Data：Middle EastCase Study.Society of Petroleum Engineers，SPE 101596.

Gregg Zelewski.et al.P-wave Seismic Azimuthal Anisotropy for Detection and Prediction of Fractures in aMiddle Eastern Carbonate Reservoir.Society of Petroleum Engineers，SPE 13903.