CN107678063B - 一种基于等级相关分析的多分量转换波裂缝预测方法 - Google Patents

一种基于等级相关分析的多分量转换波裂缝预测方法 Download PDF

Info

Publication number
CN107678063B
CN107678063B CN201710872325.2A CN201710872325A CN107678063B CN 107678063 B CN107678063 B CN 107678063B CN 201710872325 A CN201710872325 A CN 201710872325A CN 107678063 B CN107678063 B CN 107678063B
Authority
CN
China
Prior art keywords
wave
component
sampled point
processed
rank correlation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710872325.2A
Other languages
English (en)
Other versions
CN107678063A (zh
Inventor
范晓
吴伟
陈小二
刘开元
王颀
唐浩
刘璞
赵尧
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China National Petroleum Corp
BGP Inc
Original Assignee
BGP Inc
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by BGP Inc filed Critical BGP Inc
Priority to CN201710872325.2A priority Critical patent/CN107678063B/zh
Publication of CN107678063A publication Critical patent/CN107678063A/zh
Application granted granted Critical
Publication of CN107678063B publication Critical patent/CN107678063B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/362—Effecting static or dynamic corrections; Stacking
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明提供了一种基于等级相关分析的多分量转换波裂缝预测方法，所述预测方法包括以下步骤：选取计算时窗，将待处理采样点的径向分量和横向分量按照计算时窗进行截取；对待处理采样点进行等级相关分析，得到最大等级相关系数；判断待处理采样点处是否具有裂缝发育，确定裂缝发育方位；对判定有裂缝发育的待处理采样点径向分量和横向分量进行分离，得到快横波和慢横波；根据快横波和慢横波时差计算各向异性系数；重复所述选取计算时窗步骤至所述计算各向异性系数步骤，循环至下一个待处理采样点直到地震数据计算完毕。根据本发明的方法采用了判别方法对计算结果进行判别，剔除裂缝明显不发育区域，在对裂缝发育进行统计时更为准确。

Description

一种基于等级相关分析的多分量转换波裂缝预测方法
技术领域
本发明涉及石油勘探领域，具体来讲，涉及一种基于等级相关分析的多分量转换波裂缝预测方法。
背景技术
裂缝性油气储层是当今石油地球物理勘探研究的重点之一，如何精确地判定裂缝在地下的存在状态又是重点中的难点。从世界上已开发油气田统计数字来看，裂缝型储集层在油气资源和生产能力方面，约占世界总量的一半。在油藏评价、区块的评估和提高油气采收率中，识别和刻画断层与裂缝系统是至关重要的步骤。在当前裂缝预测方法中，常用的多为利用叠后地震数据计算相干或曲率等属性或属性体预测裂缝，它们多适用于预测大尺度断裂带，无法充分利用叠前信息及地层中的各向异性信息，预测结果不够精细。因此，需要探索更为精细的裂缝预测技术。
随着多波研究的不断深入，多波技术越来越多地应用到油气勘探开发中。横波分裂作为天然地震、石油勘探领域的一项具有强大潜力的技术被地质学家和地球物理学家广泛关注。横波分裂理论即当横波穿过方位各向异性介质时会形成双折射，分裂成沿平行于裂隙方向传播的S1波和沿垂直裂隙方向传播的S2波，它们各自有对应的速度，前者传播的速度较快，叫快横波，后者传播的速度较慢，叫慢横波。由于传播的速度不同，在各向异性介质中形成快慢横波时差、振幅衰减等，这些特性都与各向异性性质有关。因此可根据快、慢波的时差、波形、振幅衰减等反过来研究裂缝系统的方位和密度。这个理论对预测与油气密切相关的裂隙发育带有着十分重要的作用，因此，利用横波分裂预测储层裂缝成为一种直接可靠的方法。
横波分裂的方位和时差是研究方位各向异性和预测裂缝发育带的主要两个参数，因此，准确地计算出裂缝发育方位和时差对于准确刻画裂缝系统在地下的分布形态具有重要意义。目前，常用的方法有互相关法和能量比法。互相关法是Naville(1986)提出的一套根据互相关原理确定横波分裂方位的方法，通过对不同分量数据旋转后求解互相关系数来确定裂缝方位。能量比法是Granger(2002)提出的一种方法，该方法认为当径向分量和横向分量旋转到裂缝发育的方位时，两者的能量比值最大。能量比检测方法当方位角间隔较大时，对裂缝走向的检测精度不足，尤其对于浅层的EDA介质，由于不规则、不均匀的分方位角道集分布，计算误差较大，因此只适合宽方位角三维观测系统。以上两种方法均对计算结果没有一个过滤机制，即一些裂缝不发育的地区也计算出了方位和时差，并绘制在裂缝发育图中。并且，由于横波震源比较昂贵加之横波震源激发的横波信噪比和分辨率都较低，利用纯横波震源来进行横波分裂勘探还难以实现大规模应用。然而，多分量转换波勘探克服了纯横波勘探激发难、成本高、静校正量大等缺陷，同时具有信息丰富、兼有纵波和横波的优势、资料信噪比较高等优点，使得多分量转换波勘探成为油气储层探测的有力工具。
发明内容
针对现有技术中存在的不足，本发明的目的之一在于解决上述现有技术中存在的一个或多个问题。
为了实现上述目的，本发明提供了一种基于等级相关分析的多分量转换波裂缝预测方法，所述预测方法可以包括以下步骤：选取一个包含目标深度范围的计算时窗，将待处理采样点的径向分量地震数据和横向分量地震数据按照计算时窗进行截取；对待处理采样点进行等级相关分析，得到等级相关系数最大值；判断待处理采样点处是否具有裂缝发育，若有裂缝发育，确定裂缝发育方位；对判定有裂缝发育的待处理采样点进行快横波和慢横波分离；根据快横波和慢横波时差计算各向异性系数；重复所述选取计算时窗的步骤至所述计算各向异性系数的步骤，循环至下一个待处理采样点直到地震数据计算完毕，其中，所述对待处理采样点进行等级相关分析的步骤可包括：将截取的待处理采样点的径向分量地震数据和横向分量地震数据分别记录到数组R和数组T中；对数组R和数组T进行旋转角度θj和时移τj的变换，得到变换后的分量LR(θj,t)和LT(θj,t+τj)，其中，旋转角度变换式为：
LT(θj,t+τj)＝LT(θj,t)*τj，其中，θj＝θk+Δθ，τj＝τk+Δτ，θk为起始旋转角度，τk为起始时移，Δθ为旋转角度循环步长，Δτ为时移循环步长，t为时间；
分别求取LR(θj,t)分量和LT(θj,t+τj)分量的秩，记为ZRi和ZTi；利用式1计算等级相关系数，得到等级相关系数最大值rmax＝max{r(θj,τj)}，其中，式1为其中，Di＝ZRi-ZTi，n为分量中数据个数。
与现有技术相比，本发明的方法采用了判别方法对计算结果进行判别，剔除裂缝明显不发育区域，剔除冗余的背景信息，在对裂缝发育进行统计时更为准确，进行裂缝发育图绘制时使裂缝带的分布和发育情况更加突出。同时，本发明方法提出了一套优化的多分量转换波裂缝预测流程，对裂缝预测过程中的多个步骤进行了方法的优选，包括将等级相关理论引入到多分量转换波裂缝预测中，替代传统的互相关计算，规避了互相关对数据分布的要求，提高了计算的准确度；采用分析-判别模式对裂缝预测结果做进一步精细判断；采用多种模式进行裂缝发育的综合绘制以及井周裂缝发育的统计功能。
附图说明
通过下面结合附图进行的描述，本发明的上述和其他目的和特点将会变得更加清楚，其中：
图1示出了根据本发明示例性实施例的基于等级相关分析的多分量转换波裂缝预测方法流程图。
图2示出了根据本发明示例性实施例的坐标旋转示意图。
图3示出了根据本发明示例性实施例的横波分裂与裂缝发育方位关系示意图。
具体实施方式
在下文中，将结合附图和示例性实施例详细地描述根据本发明的基于等级相关分析的多分量转换波裂缝预测方法。
具体来讲，转换横波是指野外勘探时采用纵波激发，地震波倾斜入射到弹性分界面时会产生反射横波，其中的上行横波为转换横波。转换横波拥有横波的全部性质，当上行转换横波与裂缝走向斜交时，同样会发生分裂现象，分裂为快横波和慢横波。分裂后的快、慢波传播到地表，按照地表观测系统坐标系分解并重新合成，被多分量检波器接收。因此，当测线与裂缝方向不平行时，地面X、Y分量记录中都含有快、慢波信息，然后利用多分量记录进行快、慢波分离，根据快波与慢波的时差、波形、振幅衰减等反过来研究裂缝系统的方位、密度和时差等信息，并根据以上信息预测裂缝发育。
根据横波分裂理论，当横波穿过方位各向异性介质时会形成双折射，分裂成快横波和慢横波。在本发明的裂缝预测的方法中，方位各向异性介质可以是扩容各向异性介质(EDA)，扩容各向异性介质可近似表示由于构造应力产生空间排列垂直裂缝群体引起的各向异性。EDA介质是一种典型的方位各向异性介质，由平行的垂直裂隙、微裂隙或定向的孔隙所引起。EDA介质具有以下特点：(1)裂缝走向与地层最小应力方向垂直，与最大应力方向平行；(2)在地下足够深处，由于地层压力，最小应力通常是水平的，所以裂缝往往是垂直的；(3)偏振波沿应力大的方向传播的速度比沿应力小的方向传播速度大。
图1示出了根据本发明示例性实施例的基于等级相关分析的多分量转换波裂缝预测方法流程图。图2示出了根据本发明示例性实施例的坐标旋转示意图。图3示出了根据本发明示例性实施例的横波分裂与裂缝发育方位关系示意图。
利用本发明的方法进行裂缝预测时，需要用到三种坐标系，分别是野外采集坐标、径向-横向坐标和自然坐标。
(1)野外采集坐标/X-Y坐标系
X轴沿Inline(主测线)方向，以地震排列站号增加方向为正；Y轴沿Crossline方向；Z方向垂直向下。
(2)径向-横向坐标系/R-T坐标系
径向-横向坐标系指将X-Y坐标系采集得到的地震资料中的每个接收点的X、Y分量旋转一个角度后得到的分量。径向分量(Radial，R)：以震源为起点，并指向接收点方向为R轴正方向，接收P-SV波，随具体的地震道变化。横向分量(Transverse，T)：为R分量旋转90度的方向为T轴正方向，接收P-SH波。
(3)裂缝方位坐标系/自然坐标系/S1-S2坐标系
指转换横波分裂后，快横波偏振方向和慢横波偏振方向。快横波(S1)方向：沿裂缝走向，以与炮检方位角的夹角小于90度的裂缝走向方位角为正向，即S1走向与裂缝走向一致。慢横波(S2)方向：裂缝走向的法线方向，即垂直于裂缝发育方向。
如图1所示，在本发明的一个示例性实施例中基于等级相关分析的多分量转换波裂缝预测方法可以包括：
步骤S01，对转换波三维三分量地震数据的水平分量X和Y进行旋转，得到旋转后的转换波径向分量地震数据和横向分量地震数据，并对转换波进行各向异性处理操作。
在本示例性实施例中，首先在野外进行三维三分量地震数据的采集。在三维三分量地震勘探中，地震采集方法是采集到高质量多分量原始资料的技术保障。由于转换波下行波和上行波具有不对称的射线路径，且转换波的转换点与纵横波速度比和深度有关，因此转换波三维三分量观测系统设计有其自身的特点。
以上，在三维三分量地震勘探中可以使用三分量接收器，三分量接收器分别为水平的X、Y分量接收器，以及垂直的Z分量接收器。由于Inline和Xline方向与横波极化方向不一致，造成X、Y方向上记录的横波既有P-SV波也有P-SH波，没有明确的偏振含义，各个偏振方向上的波场相互混杂，不利于转换波资料的处理及横波分裂信息的研究和提取。因此，为了获得一致的转换波场，在处理三维三分量数据时，需要将水平方向的两个分量能量重新分配，将转换波坐标旋转，实现由野外采集坐标系(X-Y坐标系)向径向-横向坐标系(R-T坐标系)旋转。通过坐标旋转处理后，转换波主要能量分布在径向分量上，横向分量主要代表各向异性的影响。X-Y坐标系向R-T坐标系旋转的示意图如图2所示，在图2中设炮点和检波点连线方向与Inline方向夹角为θ。当径向分量R与X分量方向的夹角0°＜θ＜90°时，坐标旋转式可以为：
其中，T为横向分量。
当径向分量R与X分量方向的夹角90°＜θ＜180°时，坐标旋转式可以为：
R＝-Xcos(180°-θ)+Ysin(180°-θ)，
T＝-Xsin(180°-θ)-Ycos(180°-θ)，
表达成矩阵形式为：
对获取的原始多波数据进行处理，转换波处理是多分量转换波地震勘探能否取得成功的关键。对转换波数据处理，为了得到高信噪比、高保真、无畸变的动校道集，这里，应当注意避免使用道均衡或类似于道均衡之类的均衡手段。多分量转换波数据处理可以包括以下几个关键技术：1、转换波静校正，包括转换波表层静校正和剩余静校正；2、去噪处理；3、地表一致性处理；4、速度分析；5、叠前时间偏移处理等。在数据处理过程中，由于转换横波在近零偏移距时能量较弱，而远偏移距振幅受入射角影响会发生畸变，因此要对偏移后的道集数据进行偏移距选择，以达到较高信噪比的叠加数据。
需要说明的是，转换波静校正、去噪处理、地表一致性处理、速度分析、剩余静校正、叠前时间偏移处理等均是在转换波坐标旋转后进行处理。为了获得高信噪比、高保真、无畸变的动校道集，对转换波数据处理的方法不限于此，其它常规对转换波进行各向异性处理的操作均可。这里的转换波静校正、去噪处理、地表一致性处理、叠前时间偏移处理等均属于常规操作。
步骤S02，选取一个包含目标深度范围的计算时窗，将待处理采样点径向分量地震数据和横向分量地震数据按照计算时窗范围进行截取。
在本示例性实施例中，计算时窗的选取需包含目标深度范围，即包含有裂缝性各向异性地层，可以通过层位控制或者时间段控制目标深度范围。从待处理采样点开始，当计算时窗确定后，将经步骤S01得到的所有采样点的径向分量地震数据和横向分量地震数据分别按照确定的计算时窗范围进行截取，截取的数据分别存入不同的数组中以便进行后续的计算。
以上，待处理采样点的径向分量地震数据和横向分量地震数据是经过处理后的多分量转换波数据。在同一工区中，所有采样点的计算时窗范围选取的层位范围是一致的，例如，通过层位控制，所有采样点的时窗范围均取包含层位1与层位2，但每个采样点对应的范围值在层位1和层位2中是不一样的，比如采样点C的范围可以取1200ms(层位1)～1300ms(层位2)，采样点D的范围可以取1300ms(层位1)～1400ms(层位2)。
步骤S03，对待处理采样点进行等级相关分析，得到等级相关系数最大值。
一方面，由于传统的分析方法多建立在数据变量为定量且服从正态分布的前提下，如果上述前提条件不成立，则会造成分析结果不可信或者错误。另一方面，传统的互相关法计算的是两变量的线性相关，用线性相关描述具有一定的局限性，并且传统的互相关法在抗噪性表现一般。因此，在本发明方法的地震数据分析过程中引入了等级相关分析进行计算，等级相关分析是一种非参数统计相关的分析方法，是一种计算的波形间的相似性是非线性的分析方法。
在本示例性实施例中，所述等级相关分析可以为斯皮尔曼等级相关分析。斯皮尔等级相关也可以称为斯皮尔等级相关系数，是测定变量之间等级相关程度的一种非参数统计相关分析方法。
斯皮尔曼等级相关的基本思路可以是：假设x，y是抽自两个不同总体X，Y的样本，其观察值分别为x1,x2,x3,...,xn和y1,y2,y3,...,yn，将他们配对形成(x1,y1)，(x2,y2)，……(xn,yn)。将xi和yi各自排序(同时为升序或者降序)，分别计算出它们的秩，计做Pi和Qi，得到n对秩(P1,Q1),(P2,Q2),……(Pn,Qn)。当X和Y完全相关时，∑(Pi-Qi)＝0，可以记做∑Di＝0，其中Di就可以度量X和Y的相关性，若Di越大，X和Y之间相关越不完全。然而Di的值可以为正值，也可以为负值，直接用∑Di来测量相关性会缩小Pi和Qi之间的差值，因此可以使用∑Di 2来表示Pi和Qi的差值大小。然而∑Di 2既受到Pi和Qi不一致程度的影响，同时也受到X和Y样本个数n的影响。鉴于此，为了更加准确地进行相关性分析，采用∑Di 2的最大值去除∑Di 2，得到一个相对的测量指标，叫做斯皮尔曼等级相关，用rs表示，其中，
其中和分别是Pi和Qi的平均值。
当每个数据的秩次均不相同时，上式可以写为：
其中，Di＝Pi-Qi，即为两个观察值秩之间的差距。
在上述式中，rs的取值范围在-1到1之间。当rs大于0时，两者为正相关；当rs小于0时，两者为负相关。rs值越接近于1，样本间的相关性越大，rs越接近于0，样本间的相关性越小。
例如，在本示例性实施例中，假设采样点A(A点的线号为X，道号为Y)的多分量地震数据的采样率为S，计算时窗范围采用层位控制方法，以已知层位hor1为时窗顶，hor2为时窗底。设置该采样点处的旋转角度循环步长为Δθ，时移循环步长为Δτ，阈值范围为θk≤θ≤π，τk≤τ≤τT；其中阈值范围作为循环的终止条件，θ为旋转角度，τ为时移，对采样点A进行等级相关分析，得到等级相关系数最大值的步骤可包括：
(A)根据选取的计算时窗，将截取的待处理采样点A的径向分量和横向分量地震数据分别记录到数组R和数组T中；
(B)对数组R和数组T进行旋转角度θj和时移τj的变换，得到变换后的分量LR(θj,t)和LT(θj,t+τj)，其中，旋转角度变换式为：
LT(θj,t+τj)＝LT(θj,t)*τj
其中，θj＝θk+Δθ，τj＝τk+Δτ，θk为起始旋转角度，τk为起始时移，t为时间。
一般而言，起始角度θk可以是0°，起始时移τk可以是0ms。
(C)针对不同的旋转角度和时移，分别求取LR(θj,t)分量和LT(θj,t+τj)分量的秩，记为ZRi和ZTi；
这里，首先对LR(θj,t)分量和LT(θj,t+τj)分量进行编秩操作，然后求两者的等级相关系数。不同的旋转角度和时移对应着不同的分量，因此，当分量不同时，求取的秩不相等，可以得到不同的ZRi和ZTi。
旋转角度和时移必须在设定的阈值范围内进行循环取值，当旋转角度或时移超出阈值范围后，将停止等级相关分析计算。
(D)利用式1计算等级相关系数，得到最大值rmax＝max{r(θj,τj)}，其中，式1为
其中，Di＝ZRi-ZTi，n为分量内数据的个数。
在本发明方法的等级相关计算过程中，r(θj,τj)值越大说明旋转后的径向分量与横向分量数据越相关，当求得r(θj,τj)的最大值rmax时，此时对应的旋转角度可以初步判定是裂缝的发育方位。
步骤S04，判断待处理采样点处是否具有裂缝发育，确定最终裂缝发育方位。
为了剔除裂缝不发育区域，确保对裂缝发育进行统计时更为准确，在进行裂缝发育图绘制时使裂缝带的分布和发育情况更加突出，需要对是否具有裂缝发育作进一步精细的判断，以确定最终的裂缝发育方位。
在本示例性实施例中，对步骤S03等级相关分析结果中的最大值rmax做进一步的判断，以确定待处理采样点处是否有裂缝发育。判断方法可以包括：
步骤S04.1，将待处理采样点的最大值rmax与判别值T1进行比较，若rmax＞T1，则进入下一步；
步骤S04.2，计算步骤S03中的待处理采样点所有r(θj,τj)值的标准差Dr，若Dr＞T2，则判定该采样点处有裂缝发育，并记录最大值rmax下的旋转角度θmax和时移τmax。
其中，判别值T1，T2可以是经验值，也可以是给定值。
以上，当判定该采样点处有裂缝发育时，此时的最大值rmax对应的旋转角度即为裂缝发育的方位，时移τmax即为分裂后快慢波的时差。当判定采样点处的裂缝不发育，则将不保存该采样点的数据，不绘制裂缝发育图。所述判别值T1，T2为本领域的常规数值。
因此，在本发明方法的裂缝预测过程中，要判别采样点处是否有裂缝发育，不但需要计算得到最大的等级相关系数，同时还需要满足以上的判别条件。
步骤S05，对判定有裂缝发育的待处理采样点进行快横波和慢横波分离。
在本示例性实施例中，径向分量地震数据和横向分量地震数据均包含有快、慢横波信息，因此，需要将快、慢横波从径向分量地震数据和横向分量地震数据中分离出来。
当确认采样点处有裂缝发育时，可以采用三角变换方法或Alford旋转方法。三角变换方法和Alford旋转方法均是通过坐标旋转的方式对裂缝处的横向分量地震数据以及径向分量地震数据进行分离得到快横波和慢横波。
以上，三角变换方法即是采用式得到快横波和慢横波，其中，S1代表快横波，S2代表慢横波。
Alford旋转方法的基本思路可以包括：设裂缝走向方位角β与径向分量正向方位角α的夹角为θ＝β-α，
定义正交旋转矩阵R为：
构造接受信号矩阵：
其中，R1、T1和R2、T2分别为两个不同方位角的径向量和横向量，两个方位角之间满足正交关系。
根据Alford旋转理论：V＝RSRT
其中，S1代表快横波，S2代表慢横波。
对式V＝RSRT进行正交旋转后，得到S＝RTVR。
步骤S06，根据快横波和慢横波时差计算方位各向异性系数。
具体来讲，裂缝越发育，转换波分裂时快慢波时差越大，因此，快慢波时差可以作为裂缝发育程度的一种重要指标，用来计算方位各向异性系数。
当判定采样点处有裂缝发育后，即可进行方位各向异性系数的计算。假设分离后的快波、慢波剖面反射波的存在时间分别为TS1和TS2，则快波、慢波时差ΔT＝TS1-TS2，由此可以确定各向异性系数为：
KC＝ΔT/TS1
这里，快波、慢波时差ΔT即为待采样点处的最大值τmax。
如图3所示，图3表示横波分裂与裂缝发育方位关系图。横波分裂后快横波沿裂缝发育方位，慢横波与裂缝发育方位垂直，根据波的叠加原理，径向分量和横向分量数据中均包含有快慢波分量。
步骤S07，重复步骤S02至步骤S06，循环至下一个待处理采样点至地震数据所有采样点计算完毕。
这里，还需要对待处理采样点的线号和道号进行判断，当满足start_inline＜X＜end_inline，start_crossline＜Y＜end_crossline时，则可以循环到下一个待处理采样点至地震数据计算完毕，其中，X为待处理采样点的线号，Y为待处理采样点的道号，start_inline是工区起始线号，end_inline是工区终止线号，start_crossline是工区起始道号，end_crossline是工区终止道号。
需要说明的是，在循环处理采样点时，采样点循环的顺序可以依据线号和道号进行循环。例如，在循环过程中，先保持某线号X不变，在道号范围start_crossline＜Y＜end_crossline内依照道号依次循环，道号循环完成后，线号变为X+1，再依次循环所有道号。
步骤S08，绘制裂缝发育预测图。
在本示例性实施例中，可以采用多种模式进行裂缝发育的综合绘制以及对井周围不同范围内的裂缝发育情况进行统计分析。例如，综合裂缝发育的方位以及各向异性系数绘制裂缝发育预测图，每一采样点的裂缝发育用短线表示，短线的指示方向代表该采样点的裂缝发育方位。裂缝发育的强度可以选择以下三种方式绘制：①用直线的长短来表示；②将直线赋予一种颜色，透明度高代表裂缝不发育，透明度越低裂缝越发育；③采用不同色彩代表不同的裂缝发育程度，即用长短、透明度、颜色三种不同方式表示。此外，可以对裂缝发育带的结果进行验证，在进行验证时，需要与测井技术中的裂缝发育情况进行对比，因此可以对井周围一定范围内的裂缝发育情况进行统计，制作玫瑰图，方便与测井数据进行对比分析。
综上所述，本发明的一方面将等级相关理论引入快慢波分离计算中，用非线性相关代替了传统的线性相关，规避了传统线性相关方法对地震数据分布的要求，提高了计算精度，更好地描述了快慢波间的相关性；另一方面采用了分析-判别的流程进一步提高了裂缝预测的准确性，通过设置判别值对裂缝发育计算结果做进一步判别，与传统方法相比，提高了裂缝预测的准确度，突出了裂缝发育较强的地方，能够有效地预测断裂带发育；再一方面采用了多种显示方式综合显示了裂缝防御带的方位和强度，更加直观的显示了裂缝发育情况。
尽管上面已经通过结合示例性实施例描述了本发明，但是本领域技术人员应该清楚，在不脱离权利要求所限定的精神和范围的情况下，可对本发明的示例性实施例进行各种修改和改变。

Claims (7)

1.一种基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述预测方法包括以下步骤：
选取一个包含目标深度范围的计算时窗，将待处理采样点的径向分量地震数据和横向分量地震数据按照计算时窗进行截取；
对待处理采样点进行等级相关分析，得到等级相关系数最大值；
判断待处理采样点处是否具有裂缝发育，若有裂缝发育，确定裂缝发育方位；
对判定有裂缝发育的待处理采样点进行快横波和慢横波分离；
根据快横波和慢横波时差计算各向异性系数；
重复所述选取计算时窗的步骤至所述计算各向异性系数的步骤，循环至下一个待处理采样点直到地震数据计算完毕，其中，
所述对待处理采样点进行等级相关分析的步骤包括：
将截取的待处理采样点的径向分量地震数据和横向分量地震数据分别记录到数组R和数组T中；
对数组R和数组T进行旋转角度θj和时移τj的变换，得到变换后的分量LR(θj,t)和LT(θj,t+τj)，其中，旋转角度变换式为：
LT(θj,t+τj)＝LT(θj,t)*τj，其中，θj＝θk+Δθ，τj＝τk+Δτ，θk为起始旋转角度，τk为起始时移，Δθ为旋转角度循环步长，Δτ为时移循环步长，t为时间；
分别求取LR(θj,t)分量和LT(θj,t+τj)分量的秩，记为ZRi和ZTi；
利用式1计算等级相关系数，得到等级相关系数最大值rmax＝max{r(θj,τj)}，其中，式1为
其中，Di＝ZRi-ZTi，n为分量中数据个数；
所述判断待处理采样点处是否具有裂缝发育的步骤包括：
将待处理采样点的最大值rmax与判别值T1进行比较，若rmax＞T1，则计算所述对待处理采样点进行等级相关分析步骤中的待处理采样点所有r(θj,τj)值的标准差Dr；若Dr＞T2，则判定该待处理采样点处有裂缝发育，并记录最大值rmax下对应的旋转角度和时移，其中，判别值T1，T2为经验值或给定值。
2.根据权利要求1所述的基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述径向分量地震数据和横向分量地震数据是对三维三分量地震数据的水平分量X和Y进行旋转而得到，所述水平分量X和Y进行旋转的坐标旋转式为：
其中，θ为径向分量与水平分量X的夹角，0°＜θ＜90°或90°＜θ＜180°。
3.根据权利要求1所述的基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述等级相关分析为斯皮尔曼等级相关。
4.根据权利要求1所述的基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述快横波和慢横波分离的方法包括Alford旋转方法或三角变换方法。
5.根据权利要求1所述的基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述计算各向异性系数的方法包括：
各向异性系数KC＝ΔT/TS1，其中，ΔT＝TS1-TS2，TS1为分离后的快横波剖面反射波存在时间，TS2为分离后的慢横波剖面反射波存在时间。
6.根据权利要求1所述的基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述预测方法还包括在所述地震数据计算完毕的步骤后绘制裂缝发育预测图。
7.根据权利要求6所述的基于等级相关分析的多分量转换波裂缝预测方法，其特征在于，所述绘制裂缝发育预测图的方法包括：每一采样点的裂缝发育用短线表示，裂缝发育方位用短线的指示方向表示，裂缝发育强度用短线的长短、短线的透明度或短线的不同颜色表示。
CN201710872325.2A 2017-09-25 2017-09-25 一种基于等级相关分析的多分量转换波裂缝预测方法 Active CN107678063B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710872325.2A CN107678063B (zh) 2017-09-25 2017-09-25 一种基于等级相关分析的多分量转换波裂缝预测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710872325.2A CN107678063B (zh) 2017-09-25 2017-09-25 一种基于等级相关分析的多分量转换波裂缝预测方法

Publications (2)

Publication Number Publication Date
CN107678063A CN107678063A (zh) 2018-02-09
CN107678063B true CN107678063B (zh) 2019-11-29

Family

ID=61136977

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710872325.2A Active CN107678063B (zh) 2017-09-25 2017-09-25 一种基于等级相关分析的多分量转换波裂缝预测方法

Country Status (1)

Country Link
CN (1) CN107678063B (zh)

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108957544B (zh) * 2018-08-16 2019-07-26 中国石油大学(北京) 近地表各向异性参数的测量方法、装置、地震计及介质
CN109799529A (zh) * 2019-01-31 2019-05-24 河海大学 一种基于互相关的横波分裂vsp裂缝预测方法
CN111830563B (zh) * 2019-04-17 2023-06-30 中国石油天然气集团有限公司 纯横波多层裂缝裂缝方向及快慢波时差确定方法及装置
CN111965701B (zh) * 2019-05-20 2023-09-26 中国石油天然气集团有限公司 近地表结构反演方法及系统
CN112946736A (zh) * 2019-11-26 2021-06-11 中国石油天然气集团有限公司 一种三维观测系统重建方法及系统
CN111856580B (zh) * 2020-07-27 2022-02-15 广州海洋地质调查局 一种超远偏移距obs数据初至波能量增强方法及处理终端
CN112068198B (zh) * 2020-08-24 2022-03-18 西南科技大学 基于地震波全波形特征的裂缝破裂尺度的描述方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101738639A (zh) * 2008-11-24 2010-06-16 中国石油天然气集团公司 提高岩石裂缝参数计算精度的方法
CN102053266A (zh) * 2009-11-09 2011-05-11 中国石油化工股份有限公司 地下裂缝预测方法
CN102053277A (zh) * 2009-10-30 2011-05-11 中国石油化工股份有限公司 一种利用地震资料进行检测储层裂缝发育方向的方法
CN104898162A (zh) * 2014-03-06 2015-09-09 中国石油化工股份有限公司 地质勘探中的裂缝检测方法
CN105116448A (zh) * 2015-08-11 2015-12-02 中国石油天然气集团公司 一种转换波方位各向异性校正方法及装置

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20120092958A1 (en) * 2010-10-18 2012-04-19 Geobiz Inc. Estimation of anisotropy from compressional waves from array sonic waveforms in well logging

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101738639A (zh) * 2008-11-24 2010-06-16 中国石油天然气集团公司 提高岩石裂缝参数计算精度的方法
CN102053277A (zh) * 2009-10-30 2011-05-11 中国石油化工股份有限公司 一种利用地震资料进行检测储层裂缝发育方向的方法
CN102053266A (zh) * 2009-11-09 2011-05-11 中国石油化工股份有限公司 地下裂缝预测方法
CN104898162A (zh) * 2014-03-06 2015-09-09 中国石油化工股份有限公司 地质勘探中的裂缝检测方法
CN105116448A (zh) * 2015-08-11 2015-12-02 中国石油天然气集团公司 一种转换波方位各向异性校正方法及装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
关于等级相关中斯皮尔曼公式的性质问题;徐唐先;《统计与决策》;19951231;22-23 *
基于Alford旋转的转换波各向异性校正技术;刘红爱 等;《天然气工业》;20111231;第31卷(第5期);41-43 *

Also Published As

Publication number Publication date
CN107678063A (zh) 2018-02-09

Similar Documents

Publication Publication Date Title
CN107678063B (zh) 一种基于等级相关分析的多分量转换波裂缝预测方法
Verdon et al. Microseismic monitoring using a fiber-optic distributed acoustic sensor array
Goertz-Allmann et al. Combining microseismic and geomechanical observations to interpret storage integrity at the In Salah CCS site
CN103076623B (zh) 一种基于叠前相干的裂缝检测方法
CN103645503B (zh) 一种三维时间域照明分析及振幅补偿方法
CN104360388B (zh) 一种三维地震观测系统评价方法
Lu et al. Prediction of coal seam details and mining safety using multicomponent seismic data: A case history from China
WO2012139082A1 (en) Event selection in the image domain
Pevzner et al. Estimation of azimuthal anisotropy from VSP data using multicomponent S-wave velocity analysis
CN106556861A (zh) 一种基于全方位地震资料的方位avo反演方法
CN108957548A (zh) 一种多波多分量联合观测地震页岩气富集区预测技术
CN105954797A (zh) 地震数据的断裂识别方法和装置
CN107894616A (zh) 多分量转换波裂缝预测方法
CN103869366B (zh) 一种确定裂隙裂缝走向的方法及装置
CN105445787B (zh) 一种最优方位子体相干的裂缝预测方法
Zhang et al. Automated microseismic event location by amplitude stacking and semblance
Ji et al. Observation of higher‐mode surface waves from an active source in the Hutubi Basin, Xinjiang, China
Bar et al. Receiver function analysis reveals layered anisotropy in the crust and upper mantle beneath southern Peru and northern Bolivia
Long et al. Automatic microseismic event detection with variance fractal dimension via multitrace envelope energy stacking
Copley Seismic characterization of Niobrara fluid and rock properties: A 4D study and multicomponent (3C) analysis
Guo et al. Multimodal dispersion curves extracted from deep seismic sounding profile to characterize the sedimentary structure
CN110174696B (zh) 一种介质对称轴与观测坐标轴互换的地震波采集方法
Andrews et al. Shear-wave splitting analysis of downhole microseismic data for reservoir characterization of the Montney Formation, Pouce Coupe, Alberta
CN109425893A (zh) 一种碳酸盐缝洞体系发育带预测方法及装置
Gajek Anisotropy estimation of Lower Paleozoic shales from northern Poland using microseismic data

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
CB02 Change of applicant information

Address after: No. 216, No. 216, Huayang Avenue, Tianfu New District, Sichuan, Sichuan

Applicant after: China National Petroleum Corporation Chuanqing Drilling Engineering Geophysical Exploration Company Ltd.

Address before: 610213 No. 1, No. 1, No. 1, Huayang Avenue, Huayang Town, Shuangliu County, Chengdu, Sichuan

Applicant before: China National Petroleum Corporation Chuanqing Drilling Engineering Geophysical Exploration Company Ltd.

CB02 Change of applicant information
TA01 Transfer of patent application right

Effective date of registration: 20180402

Address after: No. 189, fan Yangxi Road, Zhuozhou City, Baoding, Hebei

Applicant after: Dongfang Geophysical Exploration Co., Ltd., China Petrochemical Corp.

Address before: No. 216, No. 216, Huayang Avenue, Tianfu New District, Sichuan, Sichuan

Applicant before: China National Petroleum Corporation Chuanqing Drilling Engineering Geophysical Exploration Company Ltd.

TA01 Transfer of patent application right
GR01 Patent grant
GR01 Patent grant
TR01 Transfer of patent right

Effective date of registration: 20201113

Address after: 100007 Beijing, Dongzhimen, North Street, No. 9, No.

Patentee after: CHINA NATIONAL PETROLEUM Corp.

Patentee after: BGP Inc., China National Petroleum Corp.

Address before: No. 189, fan Yangxi Road, Zhuozhou City, Baoding, Hebei

Patentee before: BGP Inc., China National Petroleum Corp.

TR01 Transfer of patent right