CN111736218B - 地层各向异性成因定量分析方法、设备及可读存储介质 - Google Patents
地层各向异性成因定量分析方法、设备及可读存储介质 Download PDFInfo
- Publication number
- CN111736218B CN111736218B CN202010476686.7A CN202010476686A CN111736218B CN 111736218 B CN111736218 B CN 111736218B CN 202010476686 A CN202010476686 A CN 202010476686A CN 111736218 B CN111736218 B CN 111736218B
- Authority
- CN
- China
- Prior art keywords
- wave
- transverse wave
- slowness
- fast
- slow
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 73
- 230000015572 biosynthetic process Effects 0.000 title claims abstract description 23
- 238000004445 quantitative analysis Methods 0.000 title claims abstract description 17
- 239000006185 dispersion Substances 0.000 claims description 110
- 230000005404 monopole Effects 0.000 claims description 37
- 238000004590 computer program Methods 0.000 claims description 14
- 238000004458 analytical method Methods 0.000 claims description 10
- 238000000605 extraction Methods 0.000 claims description 7
- 238000007781 pre-processing Methods 0.000 claims description 7
- 238000011084 recovery Methods 0.000 claims description 7
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000001914 filtration Methods 0.000 claims description 3
- 238000005094 computer simulation Methods 0.000 claims 1
- 238000010606 normalization Methods 0.000 claims 1
- 230000009466 transformation Effects 0.000 abstract description 5
- 230000000694 effects Effects 0.000 abstract description 4
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 abstract description 4
- 239000003345 natural gas Substances 0.000 abstract description 2
- 239000003209 petroleum derivative Substances 0.000 abstract description 2
- 238000005755 formation reaction Methods 0.000 description 20
- 230000009286 beneficial effect Effects 0.000 description 7
- 230000008859 change Effects 0.000 description 5
- 230000002159 abnormal effect Effects 0.000 description 3
- 238000012512 characterization method Methods 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 238000005553 drilling Methods 0.000 description 2
- 230000006870 function Effects 0.000 description 2
- 230000008569 process Effects 0.000 description 2
- 238000003491 array Methods 0.000 description 1
- 230000001351 cycling effect Effects 0.000 description 1
- 238000013500 data storage Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000003111 delayed effect Effects 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000009472 formulation Methods 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000035699 permeability Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
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/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/626—Physical property of subsurface with anisotropy
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
技术领域
本发明属于石油天然气勘探开发技术领域，涉及一种地层各向异性成因定量分析方法、设备及可读存储介质。
背景技术
复杂储层的储集空间多伴随裂缝发育，应力分布复杂。地层中存在裂缝或者应力分布有差异时，表征着地层存在各向异性，从测井资料中分析得到各向异性是由地层裂缝、孔洞还是地应力引起的，对后续的储层改造方案制定有直接影响作用。
目前利用阵列声波测井已经可以分析地层的各向异性，比如，发明专利US9140816，Jennifer Market等人，公开了2010-Apparatus and method for determingformation anisotropy；发明专利CN201210493762.0，伍东等人，公开了一种地层各向异性的评价方法以及Xiaoming Tang and Raghu K.Chunduru，Simultaneous inversion offormation shear-wave anisotropy parameters from cross-dipole acoustic-arraywaveform data[J].Geophysics,1999,64(5):1347-1656。但是，以上方法的主要思想是基于正交偶极四分量波形数据，结合测量仪器一号极板方位，利用Alford旋转和快慢波匹配等方法获取地层的各向异性大小及方位等信息，未对引起各向异性的原因进一步深入分析，进而导致后续的储层改造方案制定困难，不能针对性的改造储层，储层改造效果不佳。
发明内容
本发明的目的在于克服上述现有技术中未对引起储层各向异性的原因进行分析，进而导致储层改造效果不佳的缺点，提供一种地层各向异性成因定量分析方法、设备及可读存储介质。
为达到上述目的，本发明采用以下技术方案予以实现：
本发明第一方面，一种地层各向异性成因定量分析方法，包括以下步骤：
步骤1：获取储层连续深度的测井数据，包括单极全波和四分量偶极波；
步骤2：将储层当前深度的单极全波和四分量偶极波均依次进行增益恢复、延迟恢复和滤波处理，得到预处理单极全波和预处理四分量偶极波；
步骤3：根据预处理单极全波的波形由STC方法得到单极全波的横波慢度；并根据单极全波的横波慢度得到四分量偶极波的横波首波到时；
步骤4：结合当前深度采集的一号极板方位曲线，根据单极全波的横波慢度、四分量偶极波的横波首波到时和预处理四分量偶极波，反演得到快横波和慢横波的波形；
步骤5：根据快横波和慢横波的波形，通过STC方法得到快横波慢度和慢横波慢度；
步骤6：将快横波和慢横波的波形进行频散分析，得到快横波频散数据和慢横波频散数据；
步骤7：根据快横波频散数据和慢横波频散数据，拟合得到快横波频散曲线和慢横波频散曲线；
步骤8：提取频率在预设目标频率范围的快横波频散曲线和慢横波频散曲线，并进行位置判断，当快横波频散曲线和慢横波频散曲线交叉时，各向异性由地应力引起；否则，各向异性由裂缝引起；
步骤9：将预设目标频率范围的快横波频散曲线和慢横波频散曲线之间的差值求积分，用于表征当前深度的各向异性大小；
步骤10：重复步骤2至步骤9，至遍历步骤1中的储层连续深度。
本发明地层各向异性成因定量分析方法进一步的改进在于：
所述步骤3中根据预处理单极全波的波形由STC方法得到单极全波的横波慢度的具体方法为：
步骤3-1：选择单极全波的首波之后的2～3个波形周期为慢度提取时间区间；
步骤3-2：在慢度提取时间区间内设置周期窗长，并对每个周期窗长内的多组波形计算相关系数；
步骤3-3：预设慢度Si0，在慢度提取时间区间内时间窗内向后滑动周期窗长，得到一组波形相关系数数组Coef(si0，tj)；
步骤3-4：改变预设慢度，重复步骤3-3通过式(1)得到二维的相关系数数组Coef(si，tj)：
其中，M为接收站个数，Tw为模式波窗长，rm为第m个接收源接收的信号，s为慢度，Zm为第m个接收源到发射源的源距，Z1为第1个接收源到发射源的源距，τ为周期窗长起始位置；
步骤3-5：将二维的相关系数数组Coef(si，tj)进行慢度轴投影并在投影后进行归一化，然后通过式(2)得到不同慢度对应的时间慢度相关系数Slowness(j)′：
其中，Slowness(j)表示慢度为j对应的一组慢度时间相关系数，Slowness(jmax)表示慢度为j对应的一组慢度时间相关系数中的极大值，Slowness(j)′表示二维慢度时间相关系数提中得到全局极大值；
步骤3-6：Slowness(j)′的极大值对应的j即为单极全波的横波慢度。
所述步骤4中反演得到快横波和慢横波的波形以及快横波和慢横波的方位的具体方法为：通过Alford旋转法反演得到快横波和慢横波的波形以及快横波和慢横波的方位。
所述步骤6中将快横波和慢横波的波形进行频散分析时，采用的频散分析方法为Prony频散分析方法。
所述步骤7的具体方法为：
根据快横波频散数据和慢横波频散数据，通过多项式拟合方法或五点参数拟合方法拟合得到快横波频散曲线和慢横波频散曲线。
所述步骤8的具体方法为：
步骤8-1：提取频率在预设目标频率范围的快横波频散曲线和慢横波频散曲线；预设目标频率范围为0.4≤f≤2.5；
步骤8-2：确定频率为0.4、0.8、1.5和2.5的频率点，并通过式(4)计算各向异性成因指示参数index(i)：
其中，FWD(i，f)表示快横波频散曲线，SWD(i，f)表示慢横波频散曲线，f为频率；
步骤8-3：当index(i)>0时，表示快横波频散曲线和慢横波频散曲线分离，各向异性由裂缝引起；当index(i)<0时，表示快横波频散曲线和慢横波频散曲线交叉，各向异性由地应力引起。
所述步骤9的具体方法为：
步骤9-1：确定目标频率范围：0.4≤f≤2.5；
步骤9-2：通过式(5)建立基于快横波频散曲线和慢横波频散曲线的各向异性大小计算模型：
abs[(FWD(i，f)-SWD(i，f)] (5)
步骤9-3：在目标频率范围内对各向异性大小计算模型求积分，通过式(6)得到用于表征当前深度各向异性大小的Anidp(i)：
所述步骤4还包括：结合当前深度采集的一号极板方位曲线，根据单极全波的横波慢度、四分量偶极波的横波首波到时和预处理四分量偶极波，反演得到快横波和慢横波的方位；
所述步骤10还包括：通过快横波和慢横波的方位结合储层连续深度的各向异性大小得到储层最大水平地应力或者储层裂缝的延伸方向。
本发明第二方面，一种计算机设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现如上述地层各向异性成因定量分析方法的步骤。
本发明第三方面，一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现如上述地层各向异性成因定量分析方法的步骤。
与现有技术相比，本发明具有以下有益效果：
通过先对获取的波形数据预处理，再基于单极全波提取横波慢度，结合方位曲线和四分量偶极波形反演得到快横波和慢横波，将快横波和慢横波做频散处理，然后进行频散曲线拟合，分析快横波和慢横波的频散曲线的相互位置关系，进而判断引起各向异性的原因主要是地层裂缝还是地应力，最后对快慢波频散曲线差值求积分得到各向异性大小，用于定量表征由裂缝或地应力引起的各向异性大小。确定了各向异性是主要由地层裂缝引起还是地应力引起的，并定量表征了由裂缝或地应力引起的各向异性大小，进而为下一步进行储层改造提供有用信息，对储层改造措施提供指导，提高储层改造的效果。
附图说明
图1为本发明实施例的地层各向异性成因定量分析方法流程框图；
图2为本发明实施例的各向异性示意图；
图3为本发明实施例的各向异性裂缝成因单深度点的快、慢横波波形图；
图4为本发明实施例的各向异性裂缝成因单深度点的快、慢横波频散图；
图5为本发明实施例各向异性地应力成因单深度点的快、慢横波波形图；
图6为本发明实施例各向异性地应力成因单深度点的快、慢横波频散图；
图7为本发明实施例各向异性裂缝成因快横波和慢横波频散曲线拟合图；
图8为本发明实施例各向异性地应力成因快横波和慢横波频散曲线拟合图；
图9为本发明实施例的裂缝引起各向异性大小表征图；
图10为本发明实施例的地应力引起各向异性大小表征图。
具体实施方式
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。
需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。
下面结合附图对本发明做进一步详细描述：
本发明的目的在于在已经确定地层各向异性的基础上，提出一种地层各向异性成因定量分析方法，通过本发明的方法应用，进一步确定各向异性是由地层裂缝引起还是地应力引起的，为下一步储层改造提供有用信息。通过先对获取的波形数据预处理，再基于单极全波提取横波慢度，结合方位和四分量波形反演地层各向异性，对各向异性处理得到的快横波和慢横波做频散处理，然后进行频散曲线拟合，分析快横波和慢横波的频散曲线的相互位置关系，进而判断引起各向异性的原因主要是地层裂缝还是地应力，最后对快慢波频散曲线差值求积分得到各向异性大小，用于定量表征由裂缝或地应力引起的各向异性大小，具体如下：
参见图1，本发明地层各向异性成因定量分析方法，包括以下步骤：
步骤1)，获取储层连续深度的阵列声波或者交叉偶极声波的测井数据，包括四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)以及单极全波MWV(i,t)。
步骤1)的具体操作为：
11)进行阵列声波或者交叉偶极声波测井，用于处理分析的测井数据来自由一组垂直正交的X和Y两个偶极发射源和八组垂直正交的X和Y两个方向共十六个接收换能器组成的声波测井仪器。该声波测井仪器每个深度点采集的数据包括四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)，XX(i,t)表示X偶极发射源发射，X接收换能器接收；XY(i,t)表示X偶极发射源发射，Y接收换能器接收；YX(i,t)表示Y偶极发射源发射，X接收换能器接收；YY(i,t)表示Y偶极发射源发射，Y接收换能器接收。每组接收包括8组波形，对应不同源距位置上接收换能器接收到的波形。
步骤2)，对单极全波MWV(i,t)和四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)做预处理，得到预处理单极全波MWV(i,t)和预处理四分量偶极波XXO(i,t)、XYO(i,t)、YXO(i,t)和YYO(i,t)。
步骤2)的具体操作如下：
21)将对单极全波MWV(i,t)和四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)进行增益恢复。
22)将增益恢复后的单极全波MWV(i,t)和四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)进行延迟恢复。
23)将延迟恢复后的单极全波MWV(i,t)和四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)进行滤波，预设高频截止FreqH和低频截止FreqL，使用带通滤波器保留频率在FreqH和FreqL之间的单极全波MWV(i,t)和四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)，得到预处理单极全波MWV(i,t)和预处理四分量偶极波XXO(i,t)、XYO(i,t)、YXO(i,t)和YYO(i,t)。
步骤3)，对预处理单极全波MWV(i,t)利用STC方法即慢度时间相关法提取预处理单极全波MWV(i,t)的横波慢度。慢度时间相关法是一种计算多道信号之间相似性的方法，多道信号以同相轴方式排列时，通过对多道信号相似系数的计算，相似系数最高的位置连线形成的直线斜率表示模式波慢度的大小。
步骤3)的具体操作如下：
步骤31)，确定慢度提取时间区间，一般选择预处理单极全波MWV(i,t)中首波之后的2～3个波形周期。
步骤32)，对2～3个周期内的波形，对每个周期窗长WLength内的多组波形(一个接收站接收一组波形)计算相关系数。
步骤33)，固定慢度si0，在2～3个周期内时间窗内向后滑动周期窗长WLength，得到一组波形相关系数数组Coef(si0，tj)。
步骤34)，改变慢度，重复步骤33)，得到二维的相关系数数组Coef(si，tj)：
其中，M为接收站个数，Tw为窗长，rm为接收波列信号，s为慢度，Tw为信号长度，Zm为第m个接收的源距，Z1为第1个接收的源距，τ为时窗起始位置，Coef(si，tj)为相似相关系数。
步骤35)，对二维相关系数数组做慢度轴投影，并对投影后的数据归一化再乘以全局极大值：
步骤36)，在Slowness(j)′中找出极大值对应的j即为对应的预处理单极全波MWV(i,t)的横波慢度。
步骤4)，根据预处理单极全波MWV(i,t)的横波慢度计算四分量偶极波XX(i,t)、XY(i,t)、YX(i,t)和YY(i,t)的横波首波到时。
通常只处理XX(i,t)量波形的横波慢度与波至作为横波的慢度和波至：
其中，DTXX(j)表示深度j处的横波慢度，step表示仪器上提一个单位的距离，TTXX(i)深度i处的横波波至；TRmin表示发射源距离最近的接收源源距。
步骤5)，根据单个深度采集的一号极板方位曲线AZ(i)和四分量偶极波形数据XXO(i,t)、XYO(i,t)、YXO(i,t)和YYO(i,t)，利用Alford旋转等方法，反演得到快横波和慢横波的方位FAZ(i)和SWV(i,t)，以及快横波和慢横波的波形FWV(i,t)和SWV(i,t)。
步骤6)，基于STC方法分别通过快横波和慢横波的波形FWV(i,t)和SWV(i,t)计算对应的快横波慢度和慢横波慢度DTSF(i)和DTSS(i)。
步骤6)的具体操作类似步骤31)～37)。
步骤7)，基于快横波和慢横波的慢度，计算各向异性大小。
步骤7)的具体操作如下：
步骤71)，利用慢横波慢度与快横波慢度之差以及两者之和的比值，表征各向异性大小：
其中，Ani(i)表示单个深度各向异性大小。
步骤8)，根据步骤5)中各向异性反演的快慢横波FWV(i,t)和SWV(i,t)，计算快横波和慢横波的频散。
步骤8)的具体操作如下：
利用Prony方法频散分析方法，但不限于Prony方法，得到快横波和慢横波的频散离散结果。
步骤9)，基于快横波和慢横波的频散离散结果拟合，得到快横波和慢横波的连续频散曲线。
步骤9)的具体操作如下：
步骤91)，对快横波和慢横波的频散离散结果中的异常点进行剔除，如慢度异常高或异常低的值。
步骤92)，基于剔除异常点后的快横波的频散离散结果，利用多项式拟合或五点参数法等方法拟合得到快横波频散曲线FWD(i，f)；
步骤93)，基于剔除异常点后的慢横波的频散离散结果，利用多项式拟合或五点参数法等方法拟合得到慢横波频散曲线SWD(i，f)；
步骤10)对步骤9)处理得到的快横波和慢横波的连续频散曲线，判断快横波和慢横波的连续频散曲线是否交叉，给出各向异性成因指示。
步骤10)的具体实施如下：
步骤101)，在快横波和慢横波的连续频散曲线的频率-慢度坐标系内，确定目标信号主要能量频率范围，不同的仪器主要能量范围不同，如在0.4≤f≤2.5范围内，快横波和慢横波的连续频散曲线为非单调递增形态，因为低频探测深度更深，高频率信号探测深度更浅，由于钻井后，钻具破坏或储层渗透性变化，近井地层横波慢度会降低。
步骤102)，确定频率点a(f＝0.4)、b(f＝0.8)、c(f＝1.5)和d(f＝2.5)，计算各向异性成因指示参数index(i)：
步骤103)，如果index(i)>0，则各向异性主要由裂缝引起；如果index(i)<0，则各向异性主要由地应力引起。
步骤11)，基于步骤9)处理得到的快横波和慢横波的连续频散曲线定量计算各向异性大小。
步骤11)的具体操作如下：
步骤111)，确定计算频率范围：0.4≤f≤2.5；
步骤112)，确定基于快横波和慢横波的连续频散曲线的各向异性大小的计算模型：abs[(FWD(i，f)-SWD(i，f)]；
步骤113)，在频率范围0.4≤f≤2.5内对各向异性大小的计算模型求积分：
用Anidp(i)表征单个深度各向异性的大小。
步骤12)，循环到下一个深度，重复步骤2)到11)，得到连续深度的各向异性成因指示及大小。
本发明地层各向异性成因定量分析方法如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实现上述实施例方法中的全部或部分流程，也可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一计算机可读存储介质中，该计算机程序在被处理器执行时，可实现上述各个方法实施例的步骤。其中，所述计算机程序包括计算机程序代码，所述计算机程序代码可以为源代码形式、对象代码形式、可执行文件或某些中间形式等。计算机可读存储介质包括永久性和非永久性、可移动和非可移动媒体可以由任何方法或技术来实现信息存储。信息可以是计算机可读指令、数据结构、程序的模块或其他数据。需要说明的是，所述计算机可读介质包含的内容可以根据司法管辖区内立法和专利实践的要求进行适当的增减，例如在某些司法管辖区，根据立法和专利实践，计算机可读介质不包括电载波信号和电信信号。
在示例性实施例中，还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，所述计算机程序被处理器执行时实现所述地层各向异性成因定量分析方法的步骤。其中，所述计算机存储介质可以是计算机能够存取的任何可用介质或数据存储设备，包括但不限于磁性存储器(例如软盘、硬盘、磁带、磁光盘(MO)等)、光学存储器(例如CD、DVD、BD、HVD等)、以及半导体存储器(例如ROM、EPROM、EEPROM、非易失性存储器(NANDFLASH)、固态硬盘(SSD))等。
在示例性实施例中，还提供了一种计算机设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，所述处理器执行所述计算机程序时实现所述地层各向异性成因定量分析方法的步骤。处理器可能是中央处理单元(CentralProcessingUnit，CPU)，还可以是其他通用处理器、数字信号处理器(DigitalSignalProcessor，DSP)、专用集成电路(ApplicationSpecificIntegratedCircuit，ASIC)、现成可编程门阵列(Field-ProgrammableGateArray，FPGA)或者其他可编程逻辑器件、分立门或者晶体管逻辑器件、分立硬件组件等。
参见图2，各向异性处理成果图，从左至右依次第一道为深度道，第二道包含井径、自然伽马曲线，第三道为各向异性反演得到的传统慢度各向异性大小，相当于步骤4；从图2能得到反演的快慢横波、传统方法反演的各向异性大小；进而说明连续深度的常规各向异性表征参数及形式，包括井径、自然伽马常规曲线，慢度各向异性和快慢横波信息。
参见图3，各向异性较强的单个深度的8道快慢波波形；通过图3可得观察快慢横波首波波至和幅度的差异以及多道波列之间的波至和幅度变化，图中为某深度某井1390.6m的快慢波波形；通过图3说明单个深度的快慢波波形重叠显示，有利于直观观察快慢波的波至和幅度差异。
参见图4，为图3中单个深度8道快慢波处理得到的快慢波频散离散数据分布，横坐标为频率，纵坐标为慢度；如图所示，圆圈表示的为快横波频散数据分布，三角符号表示的为慢横波数据分布；从图4可以得到单个深度快慢波在频率和慢度域的数据分布，图中为某井1390.6m的快慢波频散数据分布；通过图4说明代表由裂缝引起的各向异性单个深度的快慢波频散数据在频率和慢度交会图中的分布，有利于观察快慢波频散数据的形态及相互关系。
参见图5，各向异性较强的另一个单深度的8道快慢波波形；通过图5可得观察快慢横波首波波至和幅度的差异以及多道波列之间的波至和幅度变化，图中为某深度某井1375.15m的快慢波波形；通过图5说明另一个单深度的快慢波波形重叠显示，有利于直观观察快慢波的波至和幅度差异。
参见图6，为图5中单个深度8道快慢波处理得到的快慢波频散离散数据分布，横坐标为频率，纵坐标为慢度；如图所示，圆圈表示的为快横波频散数据分布，三角符号表示的为慢横波数据分布；从图6可以得到单个深度快慢波在频率和慢度域的数据分布，图中为某井1375.15m的快慢波频散数据分布；通过图6说明代表由地应力引起的各向异性的另一个单深度的快慢波频散数据在频率和慢度交会图中的分布，有利于观察快慢波频散数据的形态及相互关系。
参见图7，为对图4中的快慢波离散频散数据拟合，分别得到快慢波频散数据拟合后的函数；从图7可以得到单个深度快慢波频散数据的拟合曲线，图中为某井1390.6m的快慢波频散数据拟合曲线；通过图7说明由裂缝引起的各向异性的一个单深度的快慢波频散数据拟合趋势，有利于预测慢度随频率的连续变化以及为定量分析快慢波频散曲线差异做前一步的准备工作。
参见图8，为对图5中的快慢波离散频散数据拟合,同样得到快慢波频散数据拟合后对应的函数；从图8可以得到单个深度快慢波频散数据的拟合曲线，图中为某井1375.15m的快慢波频散数据拟合曲线；通过图8说明由地应力引起的各向异性的另一个单深度的快慢波频散数据拟合趋势，有利于预测慢度随频率的连续变化以及为定量分析快慢波频散曲线差异做前一步的准备工作。
参见图9，为图7中快慢波拟合频散曲线在目标频率范围(一般为1.5～4kHz)之间的积分大小；从图9可以得到单个深度快慢波频散拟合曲线在特定频率范围的差异大小，用两者之间的积分大小表示，图中为某井1390.6m的快慢波频散数据拟合曲线；通过图9说明由裂缝引起的快慢波频散各向异性大小求取方式。
参见图10，为图8中中快慢波拟合频散曲线在目标频率范围(一般为1.5～4kHz)之间的积分大小；从图10可以得到单个深度快慢波频散拟合曲线在特定频率范围的差异大小，用两者之间的积分大小表示，图中为某井1375.15m的快慢波频散数据拟合曲线，通过图10说明由地应力引起的快慢波频散各向异性大小求取方式。
以上内容仅为说明本发明的技术思想，不能以此限定本发明的保护范围，凡是按照本发明提出的技术思想，在技术方案基础上所做的任何改动，均落入本发明权利要求书的保护范围之内。
Claims (3)
1.一种地层各向异性成因定量分析方法，其特征在于，包括以下步骤：
步骤1：获取储层连续深度的测井数据，包括单极全波和四分量偶极波；
步骤2：将储层当前深度的单极全波和四分量偶极波均依次进行增益恢复、延迟恢复和滤波处理，得到预处理单极全波和预处理四分量偶极波；
步骤3：根据预处理单极全波的波形由STC方法得到单极全波的横波慢度；并根据单极全波的横波慢度得到四分量偶极波的横波首波到时；
步骤4：结合当前深度采集的一号极板方位曲线，根据单极全波的横波慢度、四分量偶极波的横波首波到时和预处理四分量偶极波，反演得到快横波和慢横波的波形；
步骤5：根据快横波和慢横波的波形，通过STC方法得到快横波慢度和慢横波慢度；
步骤6：将快横波和慢横波的波形进行频散分析，得到快横波频散数据和慢横波频散数据；
步骤7：根据快横波频散数据和慢横波频散数据，拟合得到快横波频散曲线和慢横波频散曲线；
步骤8：提取频率在预设目标频率范围的快横波频散曲线和慢横波频散曲线，并进行位置判断，当快横波频散曲线和慢横波频散曲线交叉时，各向异性由地应力引起；否则，各向异性由裂缝引起；
步骤9：将预设目标频率范围的快横波频散曲线和慢横波频散曲线之间的差值求积分，用于表征当前深度的各向异性大小；
具体操作如下：
确定计算频率范围：0.4≤f≤2.5；
确定基于快横波和慢横波的连续频散曲线的各向异性大小的计算模型：abs[FWD(i，f)-SWD(i，f)]；
在频率范围0.4≤f≤2.5内对各向异性大小的计算模型求积分：
用Anidp(i)表征单个深度各向异性的大小；
步骤10：重复步骤2至步骤9，至遍历步骤1中的储层连续深度；
所述步骤3中根据预处理单极全波的波形由STC方法得到单极全波的横波慢度的具体方法为：
步骤3-1：选择单极全波的首波之后的2～3个波形周期为慢度提取时间区间；
步骤3-2：在慢度提取时间区间内设置周期窗长，并对每个周期窗长内的多组波形计算相关系数；
步骤3-3：预设慢度si0，在慢度提取时间区间内时间窗内向后滑动周期窗长，得到一组波形相关系数数组Coef(si0，tj)；
步骤3-4：改变预设慢度，重复步骤3-3通过式(1)得到二维的相关系数数组Coef(si，tj)：
其中，M为接收站个数，Tw为模式波窗长，rm为第m个接收源接收的信号，s为慢度，Zm为第m个接收源到发射源的源距，Z1为第1个接收源到发射源的源距，τ为周期窗长起始位置；
步骤3-5：将二维的相关系数数组Coef(si，tj)进行慢度轴投影并在投影后进行归一化，然后通过式(2)得到不同慢度对应的时间慢度相关系数SLowness(j)′：
其中，Slowness(j)表示慢度为j对应的一组慢度时间相关系数，Slowness(jmax)表示慢度为j对应的一组慢度时间相关系数中的极大值，Slowness(j)′表示二维慢度时间相关系数提中得到全局极大值；
步骤3-6：Slowness(j)′的极大值对应的j即为单极全波的横波慢度；
所述步骤4中反演得到快横波和慢横波的波形以及快横波和慢横波的方位的具体方法为：通过Alford旋转法反演得到快横波和慢横波的波形以及快横波和慢横波的方位；
所述步骤6中将快横波和慢横波的波形进行频散分析时，采用的频散分析方法为Prony频散分析方法；
所述步骤7的具体方法为：
根据快横波频散数据和慢横波频散数据，通过多项式拟合方法或五点参数拟合方法拟合得到快横波频散曲线和慢横波频散曲线；
所述步骤4还包括：结合当前深度采集的一号极板方位曲线，根据单极全波的横波慢度、四分量偶极波的横波首波到时和预处理四分量偶极波，反演得到快横波和慢横波的方位；
所述步骤10还包括：通过快横波和慢横波的方位结合储层连续深度的各向异性大小得到储层最大水平地应力或者储层裂缝的延伸方向。
2.一种计算机设备，包括存储器、处理器以及存储在所述存储器中并可在所述处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现如权利要求1所述方法的步骤。
3.一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现如权利要求1所述方法的步骤。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010476686.7A CN111736218B (zh) | 2020-05-29 | 2020-05-29 | 地层各向异性成因定量分析方法、设备及可读存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010476686.7A CN111736218B (zh) | 2020-05-29 | 2020-05-29 | 地层各向异性成因定量分析方法、设备及可读存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111736218A CN111736218A (zh) | 2020-10-02 |
| CN111736218B true CN111736218B (zh) | 2023-10-27 |
Family
ID=72646598
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010476686.7A Active CN111736218B (zh) | 2020-05-29 | 2020-05-29 | 地层各向异性成因定量分析方法、设备及可读存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111736218B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112360447B (zh) * | 2020-11-20 | 2024-05-28 | 中国石油天然气集团有限公司 | 一种评价储层射孔效果的方法 |
| CN112698407B (zh) * | 2020-11-23 | 2023-04-07 | 中国石油天然气集团有限公司 | 一种快速反演声波测井弯曲波频散曲线的方法 |
Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101553742A (zh) * | 2006-09-12 | 2009-10-07 | 普拉德研究及开发股份有限公司 | 使用图像和声波测井图的组合区分天然裂缝导致的声波各向异性和应力导致的声波各向异性 |
| CN110456418A (zh) * | 2019-09-12 | 2019-11-15 | 西南石油大学 | 阵列声波成像测井资料的处理和解释方法 |
| WO2020106287A1 (en) * | 2018-11-21 | 2020-05-28 | Halliburton Energy Services, Inc. | Enhanced anisotropy analysis with multi-component dipole sonic data |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7310285B2 (en) * | 2005-08-04 | 2007-12-18 | Schlumberger Technology Corporation | Method for characterizing shear wave formation anisotropy |
-
2020
- 2020-05-29 CN CN202010476686.7A patent/CN111736218B/zh active Active
Patent Citations (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101553742A (zh) * | 2006-09-12 | 2009-10-07 | 普拉德研究及开发股份有限公司 | 使用图像和声波测井图的组合区分天然裂缝导致的声波各向异性和应力导致的声波各向异性 |
| WO2020106287A1 (en) * | 2018-11-21 | 2020-05-28 | Halliburton Energy Services, Inc. | Enhanced anisotropy analysis with multi-component dipole sonic data |
| CN110456418A (zh) * | 2019-09-12 | 2019-11-15 | 西南石油大学 | 阵列声波成像测井资料的处理和解释方法 |
Non-Patent Citations (6)
| Title |
|---|
| CONTINUOUS STATISTICAL INDICATOR OF STRESS-INDUCED DIPOLE FLEXURAL WAVE ANISOTROPY;A. Syed,et al.;《SPWLA 52nd Annual Logging Symposium》;20110518;全文 * |
| Mechanical damage Detection and Anisotropy Evaluation using Dipole Sonic Dispersion Analysis;T.Plona,et al.;《SPWLA 43th Annual Logging Symposium》;20020605;全文 * |
| 交叉偶极声波测井资料在地层各向异性检测中的应用;梁晓东等;《石油仪器》;20101031(第05期);全文 * |
| 快弯曲波方位差异判断各向异性类型的方法及其应用;曾富强等;《石油学报》;20150430(第04期);全文 * |
| 横波各向异性分析方法及其在复杂储层中的应用;魏周拓等;《石油物探》;20110531;第50卷(第03期);第1-2节及附图1-8 * |
| 用正交偶极声波资料频散特性识别地层各向异性类型;苏远大等;《石油天然气学报》;20050630(第02期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111736218A (zh) | 2020-10-02 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN111736218B (zh) | 地层各向异性成因定量分析方法、设备及可读存储介质 | |
| US7423930B2 (en) | Methods and systems for detecting arrivals of interest | |
| CA2447245C (en) | Determination of dipole shear anisotropy of earth formations | |
| CN106154332B (zh) | 一种井中微地震纵横波事件初至识别方法 | |
| CN112593922B (zh) | 一种阵列声波测井评价固井二界面胶结质量的方法及装置 | |
| CN106814397B (zh) | 一种多参数联合反演计算岩石散射衰减的方法 | |
| CN102610227A (zh) | 声音信号处理设备、声音信号处理方法和程序 | |
| WO2008082828A2 (en) | Wave analysis using phase velocity processing | |
| CN107605470B (zh) | 一种纵横波径向速度变化成像方法 | |
| CN112487613B (zh) | 一种地层波走时的确定方法和装置 | |
| CN106930758B (zh) | 一种随钻声波测井方法 | |
| CN106643982A (zh) | 一种基于声波共振频率的液位测量方法 | |
| CN106908837A (zh) | 一种压裂裂缝形态和裂缝高度确定方法 | |
| CN109164492B (zh) | 一种提取套管井地层声波速度的方法 | |
| US20230026881A1 (en) | Improved Localization of an Acoustic Source | |
| CN107843919A (zh) | 一种微地震监测数据振幅均衡方法 | |
| Dong et al. | Arrival-time detection with multiscale wavelet analysis and source location of acoustic emission in rock | |
| CN107991707A (zh) | 一种基于shear let域内峰度特性的井中微地震初至波拾取方法 | |
| CN112595782B (zh) | 一种基于eemd算法的超声波横波起跳点识别方法及系统 | |
| CN112578435A (zh) | 岩石超声测试初至拾取方法及系统 | |
| EP2429214A2 (en) | Method for acoustic signal tracking | |
| CN114966823A (zh) | 一种地面微地震监测资料初至到时检测方法 | |
| US4933911A (en) | Method for determining seismic velocities | |
| CN114927141B (zh) | 异常水声信号的检测方法及系统 | |
| CN112649851A (zh) | 一种横波分裂垂直地震剖面裂缝预测方法及系统 |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |