基于阈值分割与特征提取的支撑剂砂床形态分析方法

Morphological Analysis of Proppant Sand Beds Based on Threshold Segmentation and Feature Extraction

摘要: 裂缝内砂床形态预测是水力压裂领域研究的重点和难点。目前支撑剂运移实验数据处理仍以人工方式为主，数据数量有限且效率低下，砂床形态的形成和演化亟需深入研究。本文提出了一套用于砂床形态识别和特征参数提取的方法。该方法首先利用阈值分割算法进行精准的砂床轮廓提取。随后，通过计算机视觉方法识别和量化砂床的形态特征，进而应用特征参数提取技术获得六个关键的砂床形态特征参数。这些方法能够在支撑剂运移连续时刻下实现自动化的特征提取和定量分析，无需依赖人工测量。该方法能够准确且完整地捕捉支撑剂图像中的砂床形态特征及其变化规律，为水力裂缝内支撑剂砂床形态的研究提供了新的方法和思路。通过自动化的特征提取，可以更深入地理解砂床形态的演化过程，从而推动水力压裂技术的进步和优化。

Abstract: The prediction of proppant bed morphology within fractures is a key and challenging topic in the field of hydraulic fracturing. At present, the processing of proppant transport experimental data still relies primarily on manual methods, which are inefficient and limited in data volume. The formation and evolution of the proppant bed require further in-depth investigation. This study proposes a method for identifying proppant bed morphology and extracting characteristic parameters. The method first employs a threshold segmentation algorithm to accurately extract the bed contours. Then, computer vision techniques are used to identify and quantify morphological features of the bed, followed by the application of feature extraction algorithms to obtain six key morphological parameters. These methods enable automated feature extraction and quantitative analysis at continuous time steps of proppant transport, eliminating the need for manual measurement. The proposed approach can accurately and comprehensively capture the morphological features and dynamic changes of the proppant bed in images, offering new tools and perspectives for studying bed morphology within hydraulic fractures. Through automated feature extraction, the evolution of bed morphology can be better understood, thereby promoting the advancement and optimization of hydraulic fracturing technology.

1. 引言

水力压裂过程中，裂缝内支撑剂砂床形态对储层有效改造体积和增产效果具有重要影响[1]。砂床形态预测是水力压裂领域研究的重点和难点[2]。目前，主要利用透明平行板实验和数值模拟方法获取砂床演化数据，依靠研究人员提取部分时刻下砂床图像中形态参数[3] [4]。然而，支撑剂运移过程会产生庞大的数据集，其复杂性和规模使得仅凭人工方式难以确保数据的丰富性和准确性，在描述砂床形态和演化方面更多停留在表面。对于砂床特征参数的动态特性、各因素对砂床形态影响的量化表征有待深入研究。

近年来，数字图像处理技术取得了显著的进步和广泛的应用，已成为处理纹理复杂、规模庞大的信息分析和理解问题的极佳手段[5]-[7]，为支撑剂运移过程中砂床形态和演化分析提供了新的思路。阈值分割算法与特征参数提取技术已在图像分割、特征识别和物体检测等领域获得广泛应用[8]-[10]。Liu等利用计算机视觉技术与深度学习方法相结合，构建了颗粒分割、矿物学识别、孔隙类型智能识别等模型，实现了薄片智能识别[11]。为了提高特征识别和提取的准确性，通常需要对图像实施图像分割处理。杨明旺等研究了遥感影像中建筑物体提取方法，合适的图像分割方法可以依据颜色、纹理、形状等参数将建筑物体从复杂的遥感影像中分离出来，从而大幅提高特征提取的精度[12]。孙歧峰等提出了一种基于OpenCV计算机视觉库的井喷液柱高度测量方法，实时连续测量井喷液柱高度，第一时间获取与井喷相关的关键参数[13]。目前针对水力裂缝内砂床形态识别和特征提取方法还未见文献报道。

基于先进的图像处理技术，笔者首次提出了一套基于砂床形态识别和特征参数提取的方法。构建了Otsu阈值分割算法与特征参数提取技术相结合的砂床提取方法。该方法专门针对支撑剂运移影像数据进行处理，能精确定位并提取砂床轮廓。然后利用砂床形态特征参数提取方法，获得六个关键的砂床形态特征参数，进而得到特征参数随时间的变化规律，实现砂床运移过程中特征参数的准确和连续提取，这为深入分析砂床演化特性奠定了坚实的基础。

2. 支撑剂运移砂床形态表征

图1展示了40/70目石英砂在透明狭窄水槽内运移实验。将含有石英砂的压裂液体从左侧连续泵注到水槽内，颗粒在重力作用下沉降到裂缝底部，形成初始砂床(t₁阶段)。后续颗粒越过砂床连续充填前端区域(t₂阶段)。当砂床达到平衡高度后(t₃阶段)，高度恒定，保持固有的模式向裂缝深部推进(t₄阶段)。砂床特征参数主要包括高度H、长度L、前缘角θ、延伸角β、砂床面积A。引入裂缝覆盖率R，表示砂床在高度和长度方向上的面积与裂缝总面积的比值，用于评价裂缝有效覆盖面积与体积。考虑时间因素，得到相关特征参数的动态特性。

Figure 1. Schematic of 40/70 quartz sand bed morphology and evolution

图1. 40/70石英砂砂床形态与演化示意图

3. 砂床形态特征提取方法

3.1. 支撑剂运移实验

根据水力压裂实际施工原理，本文室内实验装置采用长1500 mm，高300 mm的可视化透明亚克力板代替实际储层压裂裂缝，裂缝宽度为4 mm，更接近实际储层裂缝宽度，与数值模型几何尺寸相同。左侧入口端均匀分布有五个入口，相邻间距为65 mm，入口直径约为10 mm，射孔密度约为16孔每米，入口设计是为了模拟实际情况，即射孔连接水力裂缝。将支撑剂按一定泵速从裂缝模型五个入口注入，将出口中下段封闭模拟半封闭裂缝，利用一台高清摄像机实时记录支撑剂颗粒在平直裂缝中运移及铺置形态，当砂床达到平衡状态时，实验结束。实验原理如图2所示。

Figure 2. Proppant migration experiments schematic diagram

图2. 支撑剂运移实验原理图

3.2. Otsu阈值分割算法

Otsu阈值分割算法是一种用于图像处理的自适应阈值分割方法[14]-[16]。首先计算实验图中所有像素点的灰度值，构建一个灰度值分布图，该图描述了图像中不同灰度级出现的频率。接着，遍历所有可能的灰度值，将其作为潜在的分割阈值。对于每个阈值，图像中的像素被分为两类：灰度值小于或等于阈值的像素被视为背景，而灰度值大于阈值的像素则被视为前景。通过计算这两类像素之间的类间方差，该方差反映了前景与背景在灰度值上的差异程度。在遍历所有可能的阈值并计算其对应的类间方差后，记录下使得类间方差达到最大的阈值。该阈值代表了图像中前景与背景差异最为显著的灰度值，即为最佳阈值。依据该阈值将实验图分割为二值图像，图像分割过程如图3所示。

Figure 3. Otsu threshold segmentation process

图3. Otsu阈值分割过程

设置砂床形态图像的最大灰度值与最小灰度值分别记为g_i与g_u，初始阈值为T⁰，根据初始阈值将图像分割为前景与背景，分别求出两者的平均灰度值A_b和A_f，进而根据类间方差公式计算出对应的g值，使g值最大的阈值T即为最佳分割阈值。上述计算过程公式为：

$T^0=\frac{g_i+g_u}{2}$ (1)

$A_b={\displaystyle \sum _{g=g_i}^{T^0}g}\times h\left(g\right)/{\displaystyle \sum _{g=g_i}^{T^0}h}\left(g\right)$ (2)

$A_f={\displaystyle \sum _{g=T^0+1}^{g_u}g}\times h\left(g\right)/{\displaystyle \sum _{g=T^0+1}^{g_u}g}\times h\left(g\right)$ (3)

$g={\omega }_0\cdot {\omega }_1\cdot {\left(A_f-A_b\right)}^2$ (4)

3.3. 砂床形态特征参数提取方法

通过分析像素连通性和追踪图像边界，提取砂床形态特征参数[17] [18]，如图4所示。在二值化图像中，利用深度优先搜索算法遍历所有像素，识别连通区域并分配标签，来标记砂床面积A [19]。选择初始边界点后，根据边界判别准则从该点出发，沿边界追踪至终止条件，记录所有边界点，通过图像中每个像素所代表的实际尺寸，将图像中的像素转换为实际度量单位[20]，得到砂床形态的长度L与高度H，上述计算过程公式如下。

$L=\max \left(x_i\right)-\min \left(x_i\right)$ (5)

$H=\max \left(y_i\right)-\min \left(y_i\right)$ (6)

提取前缘角θ的方法是先获取砂床底部最左侧边缘白色像素点坐标(x_l, y_l)。在该点的x坐标值基础上增加N个像素点，获取砂床左侧轮廓中对应的纵坐标，记为(x_{l + N}, y_{l + N})。计算两点连线的斜率，得到θ值。对于延伸角β，则获取砂床底部最右侧边缘的白色像素点坐标，记为(x_r, x_r)。在该点的x坐标值基础上减少M个像素点距离，获取砂床右侧轮廓中对应的纵坐标，记为(x_{r − M}, y_{r − M})。计算两个点之间连线的斜率，得到β值，上述计算过程公式如下：

$\theta =\arctan \left(\frac{y_{l+N}-y_l}{x_{l+N}-x_l}\right)$ (7)

$\beta =\arctan \left(\frac{y_r-y_{r-M}}{x_r-x_{r-M}}\right)$ (8)

Figure 4. Extracting sand bed morphology feature parameters

图4. 砂床形态特征参数提取过程

4. 实验结果分析

4.1. 砂床形态分割算法对比

由于水力裂缝内砂床形态多样且不规则，边界部分往往容易受到噪声、遮挡或低对比度等因素的影响，导致分割结果的边界不完整或不准确。针对该问题，对比了Otsu算法、迭代阈值分割算法、高斯自适应阈值分割算法的边界处理效果。图5(a1)~(c1)中，Otsu算法分割出的砂床整体轮廓清晰且完整，迭代阈值算法分割出的轮廓左下角有明显局部缺失，分割并不完整，高斯自适应算法分割出的轮廓出现大量模糊，无法将轮廓与背景分割出来；图5(a2)~(c2)中，Otsu算法与迭代阈值法都能够较好地保留砂床边缘细节，而高斯自适应法分割后的砂床边缘外围出现不规则噪点，且边缘比较模糊；图5(c3)中，高斯自适应算法在图像分割过程中出现噪声干扰，导致背景中出现多个黑色像素块，会影响特征识别效果和精度。因此，相比于另外两种，Otsu算法能够更准确地识别砂床轮廓。

为了进一步分析三种图像分割算法对形态特征参数提取准确性的影响，对一组时长400秒的支撑剂运移影像数据进行了形态识别和特征参数提取。选择特征参数L和θ，并比较了Otsu算法、迭代阈值分割算法、高斯自适应算法在四项评价指标下的表现。这些评价指标包括准确率(Acc)、精确度(Pre)、F1得分(F1-score)和召回率(Recall)。

从表1中的数据可以看出，在长度L的特征参数中，Otsu算法相比于迭代阈值算法和高斯自适应算法都有不同程度的提升。与高斯自适应分割模型相比，Otsu算法在Acc和Recall上提升较为显著，分别提高了8.1%和10%。虽然Otsu算法与迭代阈值法在Pre上的差异仅为0.41%，但在Acc上则提升了8.1%。这一差异表明，Otsu算法通过自动确定全局最佳阈值，有效提升了图像分割的效果。前缘角θ的特征参数中，Otsu算法的Pre和F1得分相较于高斯自适应分割法分别提升了12.78%和11.28%。此外，Otsu算法在Pre和F1得分方面，分别比迭代阈值法高出9.85%和8.2%。这些数据进一步证明了Otsu算法在图像分割能力和特征识别精度方面，相较于迭代阈值算法和高斯自适应算法，表现更为优异。

Figure 5. Comparison of three image segmentation algorithms

图5. 三种图像分割模型结果对比

Table 1. Comparison of three algorithmic evaluation metrics

表1. 三种模型评价指标对比

砂床形态特征	方法	Acc/%	Pre/%	F1-score/%	Recall/%
长度L	Otsu	98.21	96.22	93.84	91.58
	迭代阈值分割	90.11	95.81	89.86	84.61
	高斯自适应分割	89.30	88.77	85.02	81.58
前缘角θ	Otsu	98.24	95.91	93.07	90.41
	迭代阈值分割	90.36	86.06	84.87	83.72
	高斯自适应分割	88.04	83.13	81.79	80.50

4.2. 砂床形态几何特征参数提取

图6展示了两个时刻下的砂床形态自动提取结果。红色矩形框界定了砂床在高度和长度方向最大范围，获得L和H值。β与θ可以自动计算并标注到轮廓图像特定位置处。A与R自动计算获得。

图7显示了一组时长为440秒的支撑剂运移数据。本文从模型中提取了四个特征数值及其变化速率，并与人工提取的结果进行了对比。结果表明，本文提出的模型计算得到的特征参数值与人工测量在特征参数和增长速率两组数据上匹配度很高，说明该模型在实际应用中可以有效替代人工测量，提高效率和准确性。在0~100秒内，四项砂床形态特征参数呈现快速增长趋势，增长速率均达到最高值。砂床的高度和长度在100秒时分别增长至189毫米和1037毫米，砂床面积在早期以平均600平方毫米/秒的速度快速增长，但在30秒后增长速率开始减缓。前缘角度在0~120秒内保持高速增长，120秒时达到35˚。在100~300秒内，各项特征参数的增长趋于平缓，砂床形态逐渐从砂丘向砂床转变，前缘角度以接近恒定的速率推进，增长速率较前期下降至0.2˚/秒。砂床的面积、长度和高度以较低速率持续增长，砂床面积增加较为明显，平均以260平方毫米/秒的速度稳定增长。在400秒后，特征参数不再变化，增长速率均降为零，表明砂床已达到平衡态。

Figure 6. Automated extraction and calibration of sand bed morphology features

图6. 砂床形态特征自动提取和标定

Figure 7. Four sand bed morphology parameters and their growth rates

图7. 四个砂床形态特征参数及增长速率

图8展示了在不同支撑剂粒径数据集(20/40目、40/70目、70/100目、100/140目石英砂)下，砂床达到平衡态时四项特征参数的误差对比情况。在这些数据集中，砂床面积的误差率相对较高，平均达到2.3%；而砂床高度和前缘角度的测量误差表现较为稳定，各自的平均误差率分别为1.75%和1.92%。总体而言，本文模型在不同支撑剂粒径数据集上表现良好，各项特征参数的误差率相对接近，显示出良好的泛化能力。

Figure 8. Comparison of errors in four key parameters of the sand bed at equilibrium under four different proppant particle sizes

图8. 四种支撑剂粒径下砂床平衡态时四项特征参数误差对比

5. 结论

(1) Otsu阈值分割算法展示了在支撑剂运移图像数据中高效准确地提取砂床轮廓的能力，不仅能够捕获砂床的边缘细节，还能精确地识别砂床的形态特征。为复杂形态裂缝中砂床形体的自动识别和特征参数提取提供了技术基础。

(2) 像素连通性分析和图像边界追踪算法在量化水力裂缝内砂床形态特征和区域分布方面发挥着关键作用。通过清晰定义的砂床轮廓，结合不同的特征参数计算方法，例如长度、面积、前缘角度等，可以有效地提取与砂床轮廓相关的特征参数。

(3) 本文提出的方法能够有效地提取连续时刻下实验图像中的砂床形态特征和相应的变化规律，误差值均在2.5%以内，表现出良好的准确性和泛化性。该方法为水力压裂过程支撑剂运移与形态展布研究提供了全新思路。

符号注释

T⁰——图像阈值分割时所设置的初始阈值；g_i——图像的最大灰度值；g_u——图像的最小灰度值；h(g) ——灰度值g出现的次数；A_b——图像前景的平均灰度值；A_f——图像背景的平均灰度值；ω₀——属于前景的像素点数占整幅图像的比例；ω₁——属于背景的像素点数占整幅图像的比例；g——经过类间方差公式计算得到的最佳阈值；max(x_i)——第i行在图像白色背景中x轴方向上最大值；min(x_i)——第i行在图像白色背景中x轴方向上最小值；max(y_i)——第i行在图像白色背景中y轴方向上最大值；min(y_i)——第i行在图像白色背景中y轴方向上最小值；(x_l, y_l)——砂床底部最左侧点；(x_r, y_r)——砂床底部最右侧点；(x_{l + N}, y_{l + N})——砂床轮廓前缘角度标记点；(x_{r − M}, y_{r − M})——砂床前进角标记点；Acc——准确率；Pre——精确度；F1-score——F1得分；Recall——召回率。

基金项目

重庆科技大学硕士研究生创新计划项目“基于数据驱动的支撑剂铺置形态预测方法研究”(YKJCX2420152)。

NOTES

^*通讯作者。