一种基于数字图像鉴定纹层状页岩不同类型裂纹扩展的方法

本发明属于纹层状页岩裂纹扩展特征研究领域，具体为一种基于数字图像鉴定纹层状页岩不同类型裂纹扩展的方法。

背景技术：

1、随着国民经济的快速发展，常规油气资源难以满足我国的能源需求，而我国页岩油气资源丰富，加快页岩油气等非常规资源的勘探与开发是我国重要的能源发展战略。页岩是一种含有层理结构的多孔介质，而我国陆相页岩普遍发育含粉砂质、凝灰质或白云质等多类型的纹层，导致岩石的各向异性和非均质性更强。纹层状页岩内部富含天然微裂缝、层理、纹层等多类宏观非连续弱面结构，这些缺陷极大地影响了岩体的力学性质，导致岩体破坏时的裂纹形态复杂。不同裂纹起裂、扩展所需的应力不同，扩展后的路径差别也很大，对实际钻井工程中的井壁稳定性和压裂工程中的人工裂缝扩展等具有重要影响。

2、目前对页岩裂缝扩展规律的研究主要针对层理页岩而未考虑其中的纹层结构，且对裂纹扩展路径的描述停留在宏观尺度、定性描述阶段，对纹层状页岩油气勘探开发的指导具有局限性。本发明基于我国陆相页岩储层普遍分布多类型纹层的现状，提出一种基于数字图像鉴定纹层状页岩不同类型裂纹扩展的方法，从微观尺度阐明纹层状页岩裂纹扩展的动态特征，并依据裂纹扩展趋势对其进行准确地分类，分类结果体现了岩体同一介质内、不同介质间裂纹扩展的差异，对钻井过程中的井壁垮塌程度预测、压裂工程的人工裂缝扩展规模预测等实际工程问题具有重要意义。

技术实现思路

1、本发明的目的在于识别纹层状页岩不同类型的裂纹扩展情况，通过识别纹层状页岩裂纹不同类型的情况，进行数字图像处理，得到其全场应力应变分析结果。

2、本发明的技术方案如下：

3、一种基于数字图像鉴定纹层状页岩不同类型裂纹扩展的方法，其实施步骤如下：

4、步骤1、纹层识别：从尺寸为φ56.6mm×80mm的岩心切割成40mm×40mm×80mm的长方体岩样，并将各面磨平抛光，然后进行端面扫描，观测岩心端面，测试岩心参数。通过纹层识别的算法得到纹层类型识别结果并得到灰度图，确定样品的纹层分布特征，选择符合条件的岩样进行下一步实验。

5、步骤2、散斑制作：将上述的长方体岩样进行表面抛光打磨，并对岩样表面进行散斑制作，随后进行力学试验，并进行图像采集。

6、步骤3、图像采集：在进行力学测试实验同时，使用高速摄像机捕捉同步图像，得到一组包含位移的图像，将采集的图像数据以及力学测试实验得到的一系列实验数据进行分析。

7、步骤4、计算裂纹位移数据：在位移图像组中，选择感兴趣区域在区域内确定网格中的位移数据，采用了可靠性导向初值传递策略算法，获取裂纹位移数据。

8、步骤5、根据裂纹数据中两侧像素点位移趋势及位移量，构建三种裂纹类型标准，建立一种裂纹类型判别的新方法，并对不同的裂纹类型进行识别。

9、进一步，所述步骤1的具体步骤如下：

10、步骤1.1、扫描仪观测岩心端面，需要满足以下条件：如图6所示d代表相邻两个纹层之间页岩基质的厚度，d代表单个纹层的厚度，选取实验样品确保：max(d)＜max(d/2)，∑d小于0.3×∑d，当实际纹层与岩样不平行时，采用蒙特卡洛算法计算非规则纹层面积，如图2所示，需要满足以下条件：s表示相邻两个纹层之间页岩基质的面积，s表示单个纹层的面积。通过max(s)＜max(s/2)，∑s小于0.3×∑s。

11、步骤1.2、筛选并得到满足条件的岩样后，首先获得岩心灰度图像，其转换公式为：

12、gray＝0.299*r+0.587*g+0.114*b  (1)

13、式中，gray为岩心灰度图像的灰度值；r为岩心灰度图像中的红色分量；g为岩心灰度图像中的绿色分量；b为岩心灰度图像中的蓝色分量。

14、然后利用二维卷积模型对灰度图像进行强化并提取纹层特征，在得到岩心图片的四分位处提取灰度曲线，利用灰度曲线的平均灰度值构建平均灰度曲线后，对平均灰度曲线进行小波变换去除噪声信号，得到小波变换处理后的平均灰度曲线。计算小波变换前后平均灰度曲线之间的灰度差值得到灰度差值曲线。再以小波变换前后的平均灰度值和灰度差值作为特征值进行聚类分析，通过首轮识别深色纹层、次轮识别亮色纹层和暗色纹层，得到纹层类型识别结果。

15、步骤1.3、利用python数据库计算岩石图片各行四分位处的灰度值并转换为整数类型得到四分位数，将岩石图片的深度和四分位数存储于python数据库中，生成三条灰度曲线，对平均灰度曲线和灰度差值曲线进行归一化处理，设置k-means算法的聚类簇数，利用k-means算法对标准化后的特征值进行聚类分析，计算相邻深度差值表示纹层厚度，并将聚类分析得到聚类结果、深度差值和特征值保存至python数据库中，最后根据聚类结果绘制差值频率直方图。通过算法得到纹层类型识别结果并得到灰度图。

16、进一步，所述步骤2的具体步骤如下：

17、步骤2.1、散斑处理：对纹层状页岩表面进行散斑制作，使用黑色散斑，散斑密排率大致为50％，即任意子区内亮暗像素所占的面积大致相同。

18、步骤2.2、实验图像采集：进行纹层状页岩单轴压缩实验时，采用两个白灯光源间隔60°对岩石正对面进行持续照射，试样表面应与相机光轴垂直，环境温度控制在15℃～25℃，拍摄环境无窗或观察窗。

19、进一步，所述步骤3的具体步骤如下：

20、步骤3.1、图像分析需使用光学图像(不包括扫描电子显微镜、原子力显微镜或x射线拍摄的图像)，调整高速相机设置(相机分辨率需要大于1080p hd/30fps)和捕捉同步图像。

21、步骤3.2、在进行力学实验时，需控制加载速率控制在0.5～1.0mpa/s。

22、步骤3.3、图像采集中的帧率的选择：为防止帧间位移过大，需要满足以下条件：选择子区大小是25pixel且图像尺寸为20pixel·mm-1，相邻两幅图像间的最大位移应小于1.25mm。试样运动速度小于1mm·s-1，最小帧率设为0.8hz。

23、步骤3.4、建立参考子集和当前子集之间相似性的度量，这是通过将最终参考子集点的灰度值与最终当前子集点的灰度值进行比较来实现的，以下两个指标是数字图像分析中两个指标，计算公式如下：

24、

25、

26、式中f和g分别是参考图像函数和当前图像函数，并返回与指定(x，y)点对应的灰度值。fm和gm分别是最终参考和当前子集的平均灰度值，定义如下：

27、

28、式中n(s)是s中的元素数。

29、步骤3.5、非线性优化，主要分为三个部分：提供初始猜测、高斯-牛顿(gn)迭代优化方案和拉格朗日插值。初始猜测，因为迭代优化方案收敛到局部最大值/最小值；因此，需要在全局最大/最小附件进行初始猜测。在gn迭代优化方案部分，讨论了两种不同的方案，正向加法高斯-牛顿法(fa-gn)和逆组合高斯-牛顿法(ic-gn)。首先讨论fa-gn方法，因为这是数字图像分析中常用的算法(以及一般的优化)。另一方面，ic-gn方法是非线性优化(通常用于图像对齐算法)的特例，它比标准的fa-gn方法更快，最近被应用于数字图像分析。最后，讨论了插值，因为需要亚像素位置的灰度值和梯度来计算fa-gn和ic-gn方法中使用的量。

30、步骤3.6、高斯-牛顿非线性迭代最小二乘法，此方案可以通过查找函数导数的根来推断为优化。更进一步，它可以推广为多元优化，方法是使用梯度代替导数，然后确定梯度范数收敛到零的位置。

31、步骤3.7、获取全场位移，能够计算位于初始参考子集中心的单个材料点的位移数据。下一步是获取位移值的完整数组。通常，这是通过首先选择感兴趣区域，然后在感兴趣区域内确定网格中的位移数据来实现的。然后，对位移进行减小或插值，形成一个“连续”的位移场。图7显示了具有复杂感兴趣区域的样品在本例中为带有裂缝的样品。感兴趣区域是白色突出显示的区域。示例子集叠加在感兴趣区域上。裂缝附近的子集是一个特例，它表明子集需要在感兴趣区域内是连续的。

32、步骤3.8、通过使用四个位移梯度获得格林-拉格朗日应变，如下所示：

33、

34、位移梯度实际上是通过ic-gn方案直接获得的。但是，这些值是嘈杂的，因此在计算应变场之前必须以某种方式进行“平滑”。一种可能的方法是在使用上述方程之前使用低通滤波器来平滑位移梯度场。其中使用的方法，它使用位移数据子集(u和v)的最小二乘平面拟合来查找以下方程中的平面参数。

35、

36、这里的想法是形成一个过约束的方程组，其形式为：

37、

38、进一步，所述步骤4的具体步骤如下：

39、步骤4.1、可靠性导向初值传递策略算法(如图9)是通过首先选择一个种子点，计算以此点为中心的周围的子集的对应变形参数，然后标记每个已经完成遍历的点，将其停用。

40、步骤4.2、子集形成一个队列，然后将变形参数和子集中心的位置储存在队列中。程序会进入一个while循环，在每次迭代时，队列中的顶部元素将会被复制并且从队列中删除，然后，将加载的队列中的数据添加到数据图中。

41、步骤4.3、接下来，使用来自加载队列点的位移数据作为非线性优化方案的初始猜测，分析周围四个像素点(如果它们还未被标记或在感兴趣范围内)。重复该程序，直到队列为空，此时将分析所有连续的感兴趣点。

42、步骤4.4、在相邻计算阶段，中间的白色十字是从队列中加载的材料点。黑色十字是添加到队列中的未分析点。灰叉已经过分析，并被跳过。这个过程具有较高准确度，首先处理不佳数据点，比如具有较低的cls的数据点(cls为变形参数)，从而防止该数据被用作相邻点的初始猜测。此外，相邻子集提供的初始猜测通常效果相当好(只要子集间距相对较小)，因为位移场通常相对平滑。此算法的唯一问题是它必须串行完成，为解决此问题的一种方法是对感兴趣区域进行分区，并选择并行处理多个种子点。图9是一个可靠性导向初值传递策略算法流程图。

43、进一步，所述步骤5的具体步骤如下：

44、步骤5.1、在裂纹附件选择几组关于中间裂纹的对称点h0和h1，将对称点在u-v坐标下的位移转化到新建立的u’-v’化坐标系(如图10所示)中，通过点h0和h1的u’和v’将裂纹分为三种类型(如图11所示)：当仅有对称点的u’为一正一负时，说明裂纹两侧的点在裂纹垂直方向做背离裂纹运动，该裂纹为拉伸裂纹；当对称点h0和h1的v’的差值>0.15时，说明裂纹两侧的点在平行与裂纹方向做相对反向运动，说明裂纹两侧的点在平行于裂纹方向做相对反向运动，该裂纹为剪切裂纹；当两对称点h0和h1的u’均为一正一负且v’的差值>0.15时，裂纹为拉剪复合裂纹。

45、u'＝u·sinβ+v·cosβ  (12)

46、v'＝u·sinβ+v·sinβ  (13)

47、综上所述，由于采用了上述技术方案，本发明的有益效果是：

48、通过识别不同类型的纹层以及不同裂纹类型来更加精准地对裂纹扩展进行应力应变分析，可以根据不同纹层来计算其应力应变，通过识别不同裂纹类型，预测裂纹扩展。通过最开始的初始猜测，进而非线性优化对算法进行更加进一步优化，采用高斯-牛顿非线性迭代最小二乘法、正向加法、逆向组合法等一系列方法最终获得不同介质的全场位移和全场应变，得到应变图像、二位云图等。最终为岩石纹层及裂纹类型进行识别达到更精准地应力应变分析，对钻井过程中的井壁垮塌程度、压裂工程的人工裂缝扩展规模预测等实际工程问题具有重要意义。