CN111311559B - 一种无偏差的裂缝表征方法 - Google Patents

Abstract

本发明涉及一种无偏差的裂缝表征方法，属于油气地质裂缝评估领域。本裂缝表征方法包括裂缝密度表征和裂缝强度表征；表征方法包括以下步骤：获取露头、地质模型或地震裂缝属性数据，采用多阈值方法进行二值化处理，得到裂缝基础图像数据；采用扫描窗对裂缝基础图像数据部分位置扫描，获得计算区域图像；在计算区域图像上标记出裂缝的端点以及裂缝与扫描窗边界的交点，根据裂缝的线性方向连接对应端点和交点，获得裂缝标识线；统计端点的数量，得到总端点数；根据总端点数和扫描窗的面积，计算得到裂缝密度；根据总端点数和扫描窗的半径，计算得到裂缝强度。本裂缝表征方法适用范围广，不受选取扫描窗大小的限制，实现过程快捷。

Description

一种无偏差的裂缝表征方法

技术领域

本发明属于油气地质裂缝评估领域，具体涉及一种无偏差的裂缝表征方法。

背景技术

了解和系统分析裂缝发育模式及其几何属性在许多相关地球科学的学科中都是十分重要的。如结构地质学、构造学、岩石物理学、地球物理学和水文地质等。

目前，基于裂缝模型、岩石裂缝露头、含裂缝数字图像等资料，描述裂缝密度和强度的几何特征的方法主要是利用数学上的一些统计方法，如概率密度函数、累计分布函数等。裂缝密度常定义为单位面积内裂缝条数(1/m²)，通过统计裂缝中心点来计算裂缝条数。

目前有学者采用直线型扫描线方法统计与直线相交的裂缝，从而确定裂缝数量、方向、轨迹长度等。利用该方法结合测井资料表征了岩露头的裂缝密度。但是这类方法若扫描线不与裂缝垂直相交，将会出现表征错误。

有部分学者采用方形或圆形扫描窗的方法，统计任意裂缝的中心点数来计算裂缝密度，在图像处理变换方法的基础上，来实现裂缝密度和强度的表征，如图1所示。这类表征方法克服了直线扫描的缺点，但会出现假象。当裂缝被扫描窗截断时，原本不应该计入统计的裂缝被记录了，这种表征方法常导致表征的裂缝密度比实际的更高，造成偏差。

发明内容

本发明为了解决上述技术问题提供一种无偏差的裂缝表征方法，本裂缝表征方法适用范围广，不受选取扫描窗大小的限制，实现过程快捷。

本发明解决上述技术问题的技术方案如下：一种无偏差的裂缝表征方法，包括裂缝密度表征和裂缝强度表征；所述表征方法包括以下步骤：

获取露头、地质模型或地震裂缝属性数据，采用多阈值方法进行二值化处理，得到裂缝基础图像数据；

采用扫描窗对所述裂缝基础图像数据部分位置扫描，获得计算区域图像；

在所述计算区域图像上标记出裂缝的端点以及裂缝与扫描窗边界的交点，根据裂缝的线性方向连接对应所述端点和所述交点，获得裂缝标识线；

统计所述端点的数量，得到总端点数；

根据所述总端点数和扫描窗的面积，计算得到裂缝密度；

根据所述总端点数和所述扫描窗的半径，计算得到裂缝强度。

本发明的有益效果是：通过标记裂缝的端点，根据端点能够直接计算出裂缝密度和裂缝强度，能够非常直接的对裂缝进行表征，不会出现偏差；本表征方法，充分考虑了裂缝大小和空间展布的综合影响而设计，能够不受扫描窗大小的限制，不受裂缝的情况的限制，能够准确，高效的对裂缝进行表征，通过裂缝密度和裂缝强度来观察裂缝发育情况，对裂缝的表征准确不存在偏差。

在上述技术方案的基础上，本发明还可以做如下改进。

进一步，所述统计所述端点的数量，得到总端点数，包括以下步骤：

统计所述裂缝标识线的数量，得到总裂缝标识线数、扫描窗内裂缝标识线数和扫描窗截断裂缝标识线数；

根据所述总裂缝标识线数、所述扫描窗内裂缝标识线数和所述扫描窗截断裂缝标识线数换算，得到总端点数。

采用上述进一步方案的有益效果是：由于端点不易统计，将其转换为对所述总裂缝标识线数、所述扫描窗内裂缝标识线数和所述扫描窗截断裂缝标识线数的统计，使得计算更方便。

进一步，所述根据所述总裂缝标识线数、所述扫描窗内裂缝标识线数和所述扫描窗截断裂缝标识线数，计算得到总端点数，具体采用如下第二公式计算得到：

第一公式为：n＝(N-N_I-N_T)+2N_I＝N+N_I-N_T

式中，N表示总裂缝标识线数，N_I表示扫描窗内裂缝标识线数，N_T表示扫描窗截断裂缝标识线数。

采用上述进一步方案的有益效果是：从而方便的得出总端点数。

进一步，所述根据所述总端点数和扫描窗的面积，计算得到裂缝密度，采用如下第二公式计算得到：

第二公式为：

式中，n^u表示裂缝上端点数，n¹表示裂缝下端点数，n表示的是裂缝端点总数，F_den表示裂缝密度，S表示扫描窗的面积；

将第一公式带入第二公式中，得到如下第三公式：

第三公式为：

采用上述进一步方案的有益效果是：能够非常直接的得到裂缝密度，得到的裂缝密度与裂缝密度定义方式计算的结果相同，不会存在偏差。

进一步，所述根据所述总端点数和所述扫描窗的半径，计算得到裂缝强度，具体包括如下步骤：

在所述计算区域图像上建立坐标系(X，Z)；

将扫描窗包围的空间根据所述坐标系分为四个象限，建立每个象限中的概率密度函数，根据所述概率密度函数，得到夹角概率密度函数；

根据所述夹角概率密度函数采用积分求取，得到一个象限裂缝长度；

将每个所述象限裂缝长度加和，得到所述总裂缝长度值；

将所述总裂缝长度值带入裂缝强度定义公式中，换算成所述总端点数和所述扫描窗的半径，得到裂缝强度。

采用上述进一步方案的有益效果是：能够方便推导出根据总端点数计算得到裂缝强度，非常方便。

进一步，所述建立每个象限中的概率密度函数，根据所述概率密度函数，得到夹角概率密度函数，包括如下步骤：

设扫描窗为圆形，扫描窗的圆心与坐标系原点重合，则扫描窗的坐标为(0，R)；

建立概率密度函数p(z)为：p(z)＝1/R，0＜z≤R；

取扫描窗的边沿与所述X坐标轴的夹角为β，得到夹角概率密度函数f(β)为：

根据几何关系换算得出如下简化夹角概率密度函数f(β)：

f(β)＝cosβ，0＜β≤π/2。

采用上述进一步方案的有益效果是：由于象限内的任意裂缝是均匀分布的，从而引入概率密度函数能够方便推导出结果。

进一步，所述根据所述夹角概率密度函数采用积分求取，得到一个象限裂缝长度，采用如下第四公式计算得到：

第四公式为：

式中，

表示一个象限裂缝的长度，

表示一个象限内在扫描窗内的裂缝数值，β表示扫描窗的边沿与水平方向X轴的夹角，f(β)表示夹角概率密度函数；

将所述简化夹角概率密度函数f(β)带入第四公式中，换算得到如下第五公式：

所述将每个所述象限裂缝长度加和，得到所述总裂缝长度值，采用如下第六公式计算得到：

式中，

表示圆形扫描窗内裂缝的数目。

采用上述进一步方案的有益效果是：从而能够计算出裂缝总长度，用于带入公式换算。

进一步，所述裂缝强度定义公式，如下第七公式：

式中，I表示裂缝强度，r表示扫描窗的半径；

将第六公式带入第七公式中，得到如下第八公式：

将第八公式经积分换算得到如下第九公式：

第九公式为：

式中，n表示的是裂缝端点总数，r表示扫描窗的半径。

采用上述进一步方案的有益效果是：从而推导出，通过总端点数来计算裂缝强度，使得对裂缝强度的表征更方便。

进一步，所述裂缝强度表征还包括：

将裂缝在坐标系(0，π)分为S份不重叠的区域，每个区域记作Δθ_i，i＝1，2，3......s，令E[n(θ_i)]为在Δθ_i方向范围内的裂缝数量；

采用不同半径(r₁、r₂、r₃，…，r_s)的扫描窗扫描同一区域，得到扫描窗的裂缝数为n₁、n₂、n₃，…，n_s，得到各自裂缝强度为I₁、I₂、I₃，…，I_s，将其周长作为各自的权重，计算得到检对裂缝强度。

采用上述进一步方案的有益效果是：能够减小误差率。

进一步，所述将其周长作为各自的权重，计算得到检对裂缝强度，采用如下第九公式计算得到：

第九公式为：

采用上述进一步方案的有益效果是：提高对裂缝强度计算的准确性。

附图说明

图1为本发明现有的裂缝密度表征示意图；

图2为本发明本裂缝密度表征示意图；

图3为本发明二值化结果示意图；

图4为本发明本裂缝密度计算的结果示意图；

图5为本发明本裂缝强度计算的结果示意图。

具体实施方式

下面通过具体实施例和对比例对本发明做进一步的详细描述，但是本领域技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。

实施例

本实施例提供一种无偏差的裂缝表征方法，包括裂缝密度表征和裂缝强度表征。表征方法包括以下步骤：

步骤101：将露头、地质模型或地震裂缝属性数据看作是一个个像素组成的图像。假设图像灰度级为K，某一像素点灰度级为i的像素有m_i(i＝0，1，2，...，K)个，则总像素有

步骤102：某一像素呈现的概率为P_i＝m_i/M，根据概率计算公式，获取概率分布直方图。

步骤103：根据直方图出现的峰值及拐点，将原始图像初始划分成N类，需要用N-1个分割阈值，设阈值集合t＝{t_k|k＝1，2，3，...，N-1}，同时令t₀＝0，t_N＝S，并给定多阈值的类内方差

具体如下第(1)公式：

公式(1)中ω_k代表第k类的概率，

代表第k类的类方差。其定义分别为：

步骤104：当类方差最小，即获取最优阈值，具体如下第(2)公式：

式中

表示最优阈值，等式右边下标0＜t₁＜t₂＜…＜t_N-1则表示初始阈值。即根据步骤1012和步骤1013获取的。

步骤105：根据上述步骤101、步骤102、步骤103和步骤104的描述，令N＝2便可实现最优二值化处理，得到裂缝基础图像数据。

步骤201：采用扫描窗对裂缝基础图像数据部分位置扫描，获得计算区域图像。

步骤301：在计算区域图像上标记出裂缝的端点以及裂缝与扫描窗边界的交点，根据裂缝的线性方向连接对应端点和交点，获得裂缝标识线。

步骤401：统计端点的数量，由于端点不易统计，将其转换为统计裂缝标识线的数量，得到总裂缝标识线数、扫描窗内裂缝标识线数和扫描窗截断裂缝标识线数。

步骤402：根据总裂缝标识线数、扫描窗内裂缝标识线数和扫描窗截断裂缝标识线数换算，得到总端点数，具体采用如下第(3)公式计算得到：

n＝(N-N_I-N_T)+2N_I＝N+N_I-N_T (3)

式中，N表示总裂缝标识线数，N_I表示扫描窗内裂缝标识线数，N_T表示扫描窗截断裂缝标识线数。

步骤501：根据公式(3)计算得到总端点数，以及扫描窗的面积，计算得到裂缝密度，具体采用如下第(4)公式计算得到：

式中，n^u表示裂缝上端点数，n¹表示裂缝下端点数，n表示的是裂缝端点总数，F_den表示裂缝密度，S表示扫描窗的面积。

步骤502：将公式(3)带入公式(4)中，得到如下第(5)公式：

从而，根据公式(5)计算得到裂缝密度F_den。

步骤601：在计算区域图像上建立坐标系(X，Z)。

步骤602：将扫描窗包围的空间根据坐标系分为四个象限，建立每个象限中的概率密度函数，根据概率密度函数，得到夹角概率密度函数，具体包括如下步骤：

步骤6021：设扫描窗为圆形，以圆形扫描窗进行扫描，扫描窗的圆心与坐标系原点重合，则扫描窗的坐标为(0，R)。

步骤6022：根据坐标(0，R)，建立概率密度函数p(z)为：

p(z)＝1/R，0＜z≤R。 (6)

步骤6023：取扫描窗的边沿与X坐标轴的夹角为β，根据概率密度函数p(z)得到夹角概率密度函数f(β)为：

步骤6024：根据几何关系z＝r sinβ，dz/dβ＝d(r sinβ)/dβ＝r cosβ，换算得出如下简化夹角概率密度函数f(β)：

f(β)＝cosβ，0＜β≤π/2。 (8)

步骤603：考虑到扫描窗与裂缝分布之间是相互独立的，而任意裂缝与扫描窗相交的概率是均匀分布的，根据夹角概率密度函数采用积分求取，得到一个象限裂缝长度，具体采用如下第(9)公式计算得到：

式中，

表示一个象限裂缝的长度，

表示一个象限内在扫描窗内的裂缝数值，β表示扫描窗的边沿与水平方向X轴的夹角，f(β)表示夹角概率密度函数。

步骤604：将简化夹角概率密度函数f(β)带入公式(9)中，换算得到如下第(10)公式：

步骤605：将每个象限裂缝长度加和，得到总裂缝长度值，具体采用如下第(11)公式计算得到：

式中，

表示圆形扫描窗内裂缝的数目。。

步骤606：将总裂缝长度值带入裂缝强度定义公式中，裂缝强度定义公式，如下第(12)公式：

式中，I表示裂缝强度，r表示扫描窗的半径。

步骤607：将公式(11)带入公式(12)中，得到如下第(13)公式：

步骤608：由于一定区域内的裂缝强度，不会随扫描窗半径大小的改变而改变，即扫描窗大小增加，其包含裂缝的长度也会相应的增加，根据其定义，裂缝强度应该保持不变，因此

可以表示为如下第(14)公式：

式中，c表示为常数。

步骤609：将公式(14)带入到公式(13)中，得到如下第(15)公式：

c＝4I。 (15)

步骤610：再将公式(15)带入到公式(14)中，用实际裂缝数n代替

将总裂缝长度值换算成总端点数和扫描窗的半径，根据总端点数和扫描窗的半径得到裂缝强度，具体如下第(16)公式：

从而通过总端点数能够直接计算出裂缝密度和裂缝强度，能够对裂缝进行表征。

步骤701：将裂缝在坐标系(0，π)分为S份不重叠的区域，每个区域记作Δθ_i，i＝1，2，3......s，令E[n(θ_i)]为在Δθ_i方向范围内的裂缝数量。

步骤702：将其带入到公式(16)中，得到如下第(17)公式：

当s趋于无穷大时，对于整个区域，公式(17)可换算为如下第(18)公式：

式中

表示裂缝与扫描窗相交的总数平均值。

步骤703：采用不同半径(r₁、r₂、r₃，…，r_s)的扫描窗扫描同一区域，得到扫描窗的裂缝数为n₁、n₂、n₃，…，n_s，得到各自裂缝强度为I₁、I₂、I₃，…，I_s，将其周长作为各自的权重，计算得到检对裂缝强度，具体采用如下第(19)公式计算得到：

将公式(18)带入到公式(19)中，得到如下第(20)公式：

即计算得出裂缝强度。

通过计算出的裂缝密度和裂缝强度能够对裂缝进行表征。

以下采用实际数据进行表征示例

如图1-图4所示，假设得到的二值化图是长宽分别为16米、20米的裂缝区域带，如图1所示，采用半径为9米的圆形扫描窗进行扫描，扫描窗中黑色实线表示裂缝。

S1：将裂缝的端点标记出来，以及裂缝与扫描窗边界的交点标记出来，根据裂缝的线性方向连接对应端点和交点，获得裂缝标识线。

S2：统计裂缝标识线的数量，得到总裂缝标识线数、扫描窗内裂缝标识线数和扫描窗截断裂缝标识线数，并带入公式(5)中，得到如下计算过程：

从而得到裂缝密度数据为0.016。

本计算结果与真实情况完全相符，不存在偏差或误差。

S3：根据统计的裂缝标识线的数量，得到总裂缝标识线数、扫描窗内裂缝标识线数和扫描窗截断裂缝标识线数，带入到公式(3)中，得到n＝4，将其带入到公式(16)中，具体如下计算过程：

得到裂缝强度为0.111。与根据定义得到的裂缝强度数据基本无差异。

图1为现有的裂缝密度表征示意图；图2为本发明本裂缝密度表征示意图；图1与图2简要说明了常规方法与本发明的比较，本发明计算得到的结果较常规方法更准确。图3是一个二值化后的裂缝模型，图4是根据图3采用本发明方法计算出的裂缝密度，图5是根据图3采用本发明方法计算出的裂缝强度，三种结合在图3中圆圈及箭头标志处，可以明显的看出，无裂缝的区域，本发明计算结果为零，在裂缝密集的地方，裂缝密度明显增大，在裂缝宽度大的地方，裂缝强度明显强，而裂缝密度大的地方裂缝强度不一定大，说明裂缝强度与裂缝密度并不一致。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (9)

1.一种无偏差的裂缝表征方法，其特征在于，包括裂缝密度表征和裂缝强度表征；所述表征方法包括以下步骤：

获取地质露头、地质模型或地震裂缝属性数据，采用多阈值方法进行二值化处理，得到裂缝基础图像数据；

采用扫描窗对所述裂缝基础图像数据部分位置扫描，获得计算区域图像；

在所述计算区域图像上标记出裂缝的端点以及裂缝与扫描窗边界的交点，根据裂缝的线性方向连接对应所述端点和所述交点，获得裂缝标识线；

统计所述裂缝标识线的数量，得到总裂缝标识线数、扫描窗内裂缝标识线数和扫描窗截断裂缝标识线数；

根据所述总裂缝标识线数、所述扫描窗内裂缝标识线数和所述扫描窗截断裂缝标识线数换算，得到总端点数；

根据所述总端点数和扫描窗的面积，计算得到裂缝密度；

根据所述总端点数和所述扫描窗的半径，计算得到裂缝强度。

2.根据权利要求1所述的无偏差的裂缝表征方法，其特征在于，所述根据所述总裂缝标识线数、所述扫描窗内裂缝标识线数和所述扫描窗截断裂缝标识线数，计算得到总端点数，具体采用如下第一公式计算得到：

第一公式为：n＝(N-N_I-N_T)+2N_I＝N+N_I-N_T

式中，N表示总裂缝标识线数，N_I表示扫描窗内裂缝标识线数，N_T表示扫描窗截断裂缝标识线数。

3.根据权利要求2所述的无偏差的裂缝表征方法，其特征在于，所述根据所述总端点数和扫描窗的面积，计算得到裂缝密度，采用如下第二公式计算得到：

第二公式为：

式中，n^u表示裂缝上端点数，n^l表示裂缝下端点数，n表示的是裂缝端点总数，F_den表示裂缝密度，S表示扫描窗的面积；

将第一公式带入第二公式中，得到如下第三公式：

第三公式为：

4.根据权利要求1所述的无偏差的裂缝表征方法，其特征在于，所述根据所述总端点数和所述扫描窗的半径，计算得到裂缝强度，具体包括如下步骤：

在所述计算区域图像上建立坐标系(X，Z)；

将扫描窗包围的空间根据所述坐标系分为四个象限，建立每个象限中的概率密度函数，根据所述概率密度函数，得到夹角概率密度函数；

根据所述夹角概率密度函数采用积分求取，得到一个象限裂缝长度；

将每个所述象限裂缝长度加和，得到总裂缝长度值；

将所述总裂缝长度值带入裂缝强度定义公式中，换算成所述总端点数和所述扫描窗的半径，得到裂缝强度。

5.根据权利要求4所述的无偏差的裂缝表征方法，其特征在于，所述建立每个象限中的概率密度函数，根据所述概率密度函数，得到夹角概率密度函数，包括如下步骤：

设扫描窗为圆形，扫描窗的圆心与坐标系原点重合，则扫描窗的半径为R；

建立概率密度函数p(z)为：p(z)＝1/R,0＜z≤R；

取扫描窗的边沿与X坐标轴的夹角为β，得到夹角概率密度函数f(β)为：

根据几何关系换算得出如下简化夹角概率密度函数f(β)：

f(β)＝cosβ,0＜β≤π/2。

6.根据权利要求5所述的无偏差的裂缝表征方法，其特征在于，所述根据所述夹角概率密度函数采用积分求取，得到一个象限裂缝长度，采用如下第四公式计算得到：

第四公式为：

式中，

表示一个象限裂缝的长度，

表示一个象限内在扫描窗内的裂缝数值，β表示扫描窗的边沿与水平方向X轴的夹角，f(β)表示夹角概率密度函数；

将所述简化夹角概率密度函数f(β)带入第四公式中，换算得到如下第五公式：

所述将每个所述象限裂缝长度加和，得到所述总裂缝长度值，采用如下第六公式计算得到：

式中，

表示圆形扫描窗内裂缝的数目。

7.根据权利要求6所述的无偏差的裂缝表征方法，其特征在于，所述裂缝强度定义公式，如下第七公式：

式中，I表示裂缝强度,r表示扫描窗的半径；

将第六公式带入第七公式中，得到如下第八公式：

将第八公式经积分换算得到如下第九公式：

第九公式为：

式中，n表示的是裂缝端点总数，r表示扫描窗的半径。

8.根据权利要求1所述的无偏差的裂缝表征方法，其特征在于，所述裂缝强度表征还包括：

将裂缝在坐标系(0，π)分为S份不重叠的区域，每个区域记作Δθ_i,i＝1,2,3......s，令E[n(θ_i)]为在Δθ_i方向范围内的裂缝数量；

采用不同半径(r₁、r₂、r₃，…，r_s)的扫描窗扫描同一区域，得到扫描窗的裂缝数为n₁、n₂、n₃，…，n_s，得到各自裂缝强度为I₁、I₂、I₃，…，I_s，将其周长作为各自的权重，计算得到模型整体的裂缝强度。

9.根据权利要求8所述的无偏差的裂缝表征方法，其特征在于，所述将其周长作为各自的权重，计算得到检对裂缝强度，采用如下第九公式计算得到：

第九公式为：

Images