一种基于电成像测井图像的页岩裂缝发育模式识别方法

1.本发明涉及油气勘探技术领域，具体涉及一种基于电成像测井图像的页岩裂缝发育模式识别方法。


背景技术：

2.页岩裂缝发育关系着页岩气的保存和运移，也是压裂改造及影响效益开发的重要因素之一。但页岩非均质性极强，发育裂缝的种类和组合样式复杂多样，由于尺寸原因，难以将单个或某些裂缝的特征与测井曲线响应完全对应起来。因此在复杂裂缝发育段中归纳并识别裂缝的主要发育模式及特征可更好的对压裂开发提供指导。目前可以通过电成像测井图像提取较为清晰明显的构造裂缝和层理信息，但仍有许多不明显的裂缝信息未能在图像上显示，给页岩气储层裂缝评价和分析带来了很大挑战。
3.目前页岩裂缝段的识别预测最主要的还是通过阵列声波等特殊测井的方式，大量地质工作者也从常规测井资料出发借助信号分析和深度学习等方法建立了裂缝预测模型。常规测井和阵列声波资料难以达到提取裂缝类型、数量、产状等信息，因此以上资料和方法对裂缝的评价局限于定性评价有无裂缝和裂缝发育程度。因电成像测井图像有着分辨率高的优势，陆续出现了通过电成像测井图像提取裂缝的产状信息。但众学者对信息复杂或不清晰的电成像测井图像上的裂缝信息，尤其是密集的层理信息的提取效果不明显。且成像上对裂缝的研究往往集中于图像处理方面，忽视了裂缝的类型和组合特征，而储层段不同的裂缝类型、数量和发育特征对后期压裂改造都有很大影响。
4.现有方法如专利申请号为2019110288605，一种细粒沉积陆相页岩裂缝识别评价方法，包括以下步骤：s1：提取测井属性，包括小波高频属性whf、小波奇异属性wsv和预测滤波误差属性pfe；s2：构建页岩裂缝评价子参数；s3：采用变异系数加权法，构建页岩裂缝评价综合子参数；s4：构建页岩裂缝评价综合参数；该方法实现单井页岩裂缝的识别和裂缝段内裂缝发育程度的评价，解决了利用常规测井资料识别裂缝段的问题。但是未涉及裂缝的类型和产状等信息，仅对某段是否有裂缝和裂缝发育程度进行定性的评价，对后期压裂改造的指导意义具有局限性。
5.如专利申请号为2019103754535，基于超声成像的油基泥浆井裂缝定量参数智能计算方法，在分析超声成像裂缝响应特征的基础上，通过滤波处理突显裂缝响应，提高信噪比，然后通过图像智能识别算法自动拾取裂缝轨迹并计算裂缝长度、密度和倾角，接着沿裂缝轨迹法向方向自动探索并统计计算裂缝宽度，再根据裂缝长度、宽度、密度参数计算裂缝面孔率。利用人工或计算机自动处理岩心滚扫描图片获取岩心裂缝参数。利用岩心裂缝参数对超声成像计算的裂缝参数进行刻度校正，可提高复杂非常规储层测井评价的准确性和可靠性，进而为制定合理高效的开发方案提供依据。该方法虽然实现了对裂缝密度和产状信息的提取，但没有涉及到不同裂缝类型发育的组合特征。


技术实现要素：

6.本发明针对现有技术存在的问题提供一种基于电成像测井图像的页岩裂缝发育模式识别方法。
7.本发明采用的技术方案是：
8.一种基于电成像测井图像的页岩裂缝发育模式识别方法，包括以下步骤：
9.步骤1：获取电成像，进行双边滤波处理；
10.步骤2：电成像二值化，然后通过改进霍夫变换提取构造裂缝的产状信息，统计每米的构造裂缝条数和开度；
11.步骤3：将步骤1滤波后的图像进行灰度化和锐化处理，统计每米层理条数n；
12.步骤4：根据步骤2得到的开度计算每米构造裂缝开度之和，开度之和除以每段层理条数n得到σ；σ乘以n最大值记作k；
13.以n为横坐标，以σ为纵坐标构造直角坐标系；在该直角坐标系下做曲线σ＝k/n；将直角坐标系的横坐标和纵坐标及曲线形成的区域划分为六种区域，对应六种裂缝发育模式；将样本点投影到直角坐标系中，根据其投影的位置进行裂缝发育模式的识别。
14.进一步的，所述步骤1中的双边滤波处理过程如下：
[0015][0016]
式中：i为输入图像，bf为双边滤波后的图像，p为滤波核窗口的中心坐标，q为过滤核窗的非中心坐标，s为空间域，w
p
为归一化因子，为空间域权重系数，为值域权重系数，iq为对应q坐标的像素值。
[0017]
进一步的，所述步骤2中提取构造裂缝的产状信息过程如下：通过改进霍夫变换提取固定周期正弦线基线来提取构造裂缝的产状信息；
[0018]
基线的确定方法如下：
[0019]
s1：将图像分为左右两等分i1、i2；
[0020]
s2：设置一个一维累加器数组acc[y0]，将其各元素初始化为0；其中y0为正弦曲线的基线位置；
[0021]
s3：对于i1中的任意边界点p(x
p
，y
p
)，执行步骤s4～s7；
[0022]
s4：对于i2中垂直扫描线x＝x
p
上任意边缘点p
′
(x
p
′
,y
p
′
)，将p和p
′
配成点对，执行步骤 s5～s7；
[0023]
s5：计算点对p
′
的终点纵坐标yc；
[0024]
s6：与yc对应的累加器数组元素的数值加1，acc[y0]＝acc[y0] 1；
[0025]
s7：循环直到图像中所有的点对均被访问；
[0026]
s8：累加器数组的峰值位置即为正弦线基线位置。
[0027]
进一步的，所述正弦曲线如下：
[0028]
y0＝y-a sin(ωx-β)
[0029]
式中：x，y为对应图像空间中的像素点坐标，a为正弦线的振幅，ω为角速度，β为初相位，y0为正弦曲线的基线位置。
[0030]
进一步的，所述裂缝开度计算方法如下：
[0031][0032]
式中：wi为纵向上该条裂缝所占像素点个数，m为等间隔位置个数，为等间隔在m 个位置上求取该裂缝纵向上覆盖像素点的平均个数，wd为单个像素点的实际高度。
[0033]
进一步的，所述灰度化和锐化过程如下：
[0034]
灰度化采用matlab中的rgb2gray函数进行处理；锐化采用matlab中的sobel算子和imfilter函数进行处理。
[0035]
进一步的，所述步骤4中六种裂缝发育模式分别为：裂缝与层理混合发育型、中密度层理型、密集层理型、强构造裂缝弱层理型、大尺度构造裂缝型、致密段。
[0036]
本发明采用的技术方案是：
[0037]
(1)本发明裂缝模式识别方法能够精细地拾取成像上的裂缝信息并准确的识别该段的裂缝发育模式，为成像裂缝评价提供了新的思路，为后期压裂方案优化提供了依据；
[0038]
(2)本发明通过双边滤波，能够实现降噪且能保留边缘信息，使得裂缝信息的提取更加准确；
[0039]
(3)本发明通过图像灰度和锐化处理，突出了成像上的层理信息，使人机交互统计层理信息更加便捷与准确；
[0040]
(4)本发明采用1m段为一样本点研究，重点研究和突出该段裂缝发育的模式与特征，淡化了裂缝信息提取方面误差的影响。
附图说明
[0041]
图1为本发明方法流程示意图。
[0042]
图2为本发明实施例1中的改进霍夫变换提取固定周期正弦线的原理示意图。
[0043]
图3为本发明实施例1中的改进霍夫变换提取固定周期正弦线的原理应用在简单图像(图3a构造裂缝)和复杂图像(图3b层理信息较多)上的效果图。
[0044]
图4为本发明实施例1中的电成像测井图(a)、双边滤波(b)、灰度图锐化后(c)和二值图(d)的效果图。
[0045]
图5为本发明实施例1中的裂缝发育模式类型的示意图。
[0046]
图6为本发明实施例1中的不同裂缝发育模式在σ和n坐标系中的区域划分概念图。
[0047]
图7为本发明实施例1中将样本点投影到σ和n坐标系中，确定不同裂缝发育模式的划分界限图。
[0048]
图8为本发明实施例1中深层页岩裂缝发育模式识别方法应用效果图。
具体实施方式
[0049]
下面结合附图和具体实施例对本发明做进一步说明。
[0050]
如图1所示，获取电成像图，对电成像图进行双边滤波，双边滤波过程如下：
[0051][0052]
式中：i为输入图像，bf为双边滤波后的图像，p为滤波核窗口的中心坐标，q为过滤核窗的非中心坐标，s为空间域，w
p
为归一化因子，为空间域权重系数，为值域权重系数，iq为对应于q坐标的像素值。
[0053]
当窗口出现边缘情况，权值变大，使得边缘信息保留。当窗口中像素值平滑时，权重变小，使得区域内被降噪处理。电成像图和双边滤波后的图像如图4b所示。
[0054]
步骤2：电成像二值化，然后通过改进霍夫变换提取构造裂缝的产状信息，统计每米的构造裂缝条数和开度；
[0055]
改进的hough变换通过图像空间域参数空间的点线对偶性，将直线方程改为正弦曲线方程，其正弦曲线为：
[0056]
y0＝y-a sin(ωx-β)
[0057]
式中：x、y为对应图像空间中的像素点，a为正弦线的振幅，ω为角速度，β为初相位， y0为正弦曲线的基线位置，基线的确定方法如下：
[0058]
基线的确定方法如下：
[0059]
s1：将图像分为左右两等分i1、i2；
[0060]
s2：设置一个一维累加器数组acc[y0]，将其各元素初始化为0；其中y0为正弦曲线的基线位置；
[0061]
s3：对于i1中的任意边界点p(x
p
，y
p
)，执行步骤s4～s7；
[0062]
s4：对于i2中垂直扫描线x＝x
p
上任意边缘点p
′
(x
p
′
,y
p
′
)，将p和p
′
配成点对，执行步骤 s5～s7；
[0063]
s5：计算点对p
′
的终点纵坐标yc；
[0064]
s6：与yc对应的累加器数组元素的数值加1，acc[y0]＝acc[y0] 1；
[0065]
s7：循环直到图像中所有的点对均被访问；
[0066]
s8：累加器数组的峰值位置即为正弦线基线位置。
[0067]
上述过程能够快速且有效的识别构造裂缝在电成像测井图像上正弦线的基线(图3a)，但对于层理较多的图片，该方法不能有效的提取层理信息(如图3b所示)。
[0068]
构造裂缝的开度计算公式如下：
[0069][0070]
式中：wi为纵向上该条裂缝所占像素点个数，m为等间隔位置个数，为等间隔在m 个位置上求取该裂缝纵向上覆盖像素点的平均个数，wd为单个像素点的实际高度。本实施例中m取8。
[0071]
步骤3：将步骤1滤波后的图像进行灰度化和锐化处理，通过人机交互统计每米层理条数n；
[0072]
对于图像不清晰且较为复杂、层理较多的图片，经过灰度化和锐化处理后可突出层理信息(如图4c所示)。通过人机交互方式统计每米层理条数更加准确。
[0073]
灰度化采用matlab中的rgb2gray函数进行处理；
[0074]
调用函数方式如下：
[0075]
a_gray＝rgb2gray(bf)
[0076]
式中：bf为双边滤波后的图像，a_gray为灰度化后的图像。
[0077]
锐化采用matlab中的sobel算子和imfilter函数进行处理，锐化过程如下：
[0078]
f_sobel＝fspecial(
′
sobel
′
)
[0079]
b＝imfilter(a_gray,f_sobel)
[0080]
a_sobel＝a_gray-b
[0081]
式中：f_sobel为生成的滤波函数，a_gray为灰度化后的图像，b为滤除的信息，a_sobel 为锐化后突出细节的灰度图。
[0082]
根据步骤2得到的开度计算每米构造裂缝开度之和，开度之和除以每段层理条数n得到σ；σ乘以n最大值记作k；
[0083]
以n为横坐标，以σ为纵坐标构造直角坐标系；在该直角坐标系下做曲线σ＝k/n；将直角坐标系的横坐标和纵坐标及曲线形成的区域划分为六种区域，对应六种裂缝发育模式；将样本点投影到直角坐标系中，根据其投影的位置进行裂缝发育模式的识别。
[0084]
通过大量调研和岩心、成像观察，结合构造裂缝和层理发育和组合的特征将深层页岩裂缝按照发育模式分为了六大类(如图5所示)；分别为：构造裂缝与层理混合发育型(i类)、中密度层理型(ⅱ类)、密集层理型(ⅲ类)、强构造裂缝弱层理型(ⅳ类)、大尺度构造裂缝型(
ⅴ
类)、致密段(vi类)。划分以上六种裂缝发育模式的方法具体过程如下：
[0085]
s11：将电成像测井图像分为以1m为一段的图像，如步骤1～3方法统计每小段(1m) 的构造裂缝开度和(假如该1m段有j条构造裂缝，则该段构造裂缝开度和为)和每小段的层理条数n。
[0086]
s12：计算每段构造裂缝开度和除以层理条数的值(σ)：
[0087][0088]
s13：计算每段σ乘n的最大值，记作k；
[0089]
s14：根据六种裂缝发育模式都有各自的σ和n的特征，建立以σ为纵坐标、n为横坐标的直角坐标系中，并将区域(n》0,σ》0,σ＝k/n)划分为六个分区域(如图6所示)。每个区域对应一种裂缝发育模式(六种裂缝发育模式分别为：构造裂缝与层理混合发育型(i类)、中密度层理型(ⅱ类)、密度层理型(ⅲ类)、强构造裂缝弱层理型(ⅳ类)、大尺度构造裂缝型(
ⅴ
类)、致密段(vi类))。其中i类n、σ值较大，ⅱ类n值较大，σ值较小，ⅲ类n值最大，σ值很小，ⅳ类n值较小，σ值很大，
ⅴ
类裂缝n值最小，σ值最大，vi类裂缝欠发育n值和σ值都最小。
[0090]
s15：将每小段作为一个样本点并投影在以σ为纵坐标、n为横坐标的直角坐标系中(如图7所示)。每个样本点的投影都应该在σ＝k/n曲线的限制下(n》0,σ》0)。根据六种裂缝发育模式都有各自的σ和n的特征结合电成像测井图像确定了每种裂缝发育模式的σ和n的
界限 (如表1所示)。完成了成像裂缝发育模式识别方法的建立。
[0091]
表1.不同裂缝发育模式类型界限
[0092][0093]
本发明公开的基于电成像测井图像的页岩裂缝发育模式识别方法，能够精细的拾取成像上的裂缝信息并准确的识别该段的裂缝发育模式，为成像裂缝评价提供了新的思路；本发明方法通过双边滤波，能够实现降噪且能保留边缘信息，使得裂缝信息的提取更加准确。通过将图像灰度和锐化处理，突出了成像上的层理信息，使人机交互统计层理信息更加便捷与准确。以1m段为一样本点研究，重点研究和突出该段裂缝发育的模式与特征，淡化了复杂裂缝和层理信息提取方面误差的影响，为后期压裂方案优化提供了依据。