基于颗粒识别分析系统的碎屑流堆积物颗粒识别和统计方法研究

陈 达，许 强，郑 光，彭双麒，王 卓，何 攀

（成都理工大学地质灾害防治与地质环境保护国家重点实验室，四川 成都 610059）

碎屑流是一种高速远程溃散性滑坡，具有非常高的运动速度和极大的位移，速度在30 m/s 以上[1]，具有极大的破坏力。为研究碎屑流运动特性及运动机理，有必要进行堆积物粒度定量分析，这也是碎屑流灾害研究的重点。传统颗粒的测量主要依靠人工，如需通过生成高分辨率数字地表模型（Digital Surface Model，DSM）对堆积物粒度进行人工解译[2-4]，存在效率低、精度不高以及人为因素误差大等不足。

计算机数字图像处理技术具有自动化、高效快捷等特点，被广泛运用于颗粒识别。吴义祥[5]研究了土细观结构图像的定量分析系统，可以获取细小颗粒面积、最长弦长度等基本要素；施斌等[6]使用Videolab图像处理系统直接分析SEM 等照片，定量处理土颗粒微细形状、大小和比例关系；涂新斌等[7]通过Qwin软件处理颗粒二值图像，获得了颗粒数、长度、宽度等颗粒数字特征；梁双华等[8]提出Mapinfo-Photoshop计算方法，利用Photoshop 图像处理功能和Mapinfo的统计分析功能，对土颗粒图像进行信息提取，定量表征颗粒、孔隙的个数以及周长、平均孔径等参数，极大地降低了图像分析费用，提高效率；Bai Baojun 等[9]通过Imagej 软件，直接将页岩原始孔隙图像进行二值化处理，进而定量统计孔隙、颗粒发育情况，得出孔隙率。可以看出，这些研究都是为了定量化分析图像的细观特征，实现图像自动识别，得出孔隙、颗粒的结构参数和数量特征。但是，这些方法在以下关键问题上存在一定的缺陷：（1）设定研究最小颗粒粒径及较小粒径颗粒识别方法；（2）一个图像中可能有上千个颗粒和孔隙，形状复杂，分布不均，宏观特征参数已经不能满足研究需要，为此需要一些特定的参数，分析每个颗粒的形状、大小，以此来反映颗粒系统特征。

孔隙（颗粒）及裂隙图像识别与分析系统（Pores and Cracks Analysis System，PCAS），可以定量分析碎屑流堆积物粒度特征，准确获得堆积物颗粒结构的几何参数和统计参数，且操作简单，具有自动化和可重复等优点。PCAS 在颗粒识别中，存在着三个参数的选取问题：阈值（T）、孔喉封闭半径（r）、最小孔隙面积（S0）。该方法在滑坡堆积物粒度分析中具有很好的效果，但实用范例很少，目前仅彭双麒等[10-11]通过PCAS 系统研究了碎屑流堆积体，得出碎屑流堆积物颗粒分布的一般规律，但没有具体阐述PCAS颗粒识别系统的工作原理、参数意义以及参数选取方法。为此，本文通过PCAS 软件对贵州纳雍崩塌碎屑流堆积物颗粒特征进行分析，得出图像识别系统中阈值、孔喉封闭半径、最小孔隙面积合适的参数取值，为PCAS 在崩塌堆积物颗粒的定量研究中提供了一个高效可行的方案。

1 方法介绍

PCAS 图像识别步骤见图1。

首先将图像导入PCAS 软件，调整阈值（T）进行二值化处理，区分出孔隙和颗粒；然后选取合适的孔喉封闭半径（r）和最小孔隙面积（S0），最后孔隙和颗粒的各种几何参数将汇集于数据表，包括个数、长度、宽度、面积、形状系数、定向性等，并统计得到颗粒和孔隙的分形维数、表观孔隙率、面积概率分布指数等统计参数，实现孔隙和颗粒结构的定量分析。本次主要研究颗粒个数和长度两个参数。

PCAS 得出的参数都是以像素为单位，可以通过换算得到实际的几何参数，如式（1）、（2）[12]：

式中：R—图像的分辨率；

S、C—像素面积和像素周长；

St、Ct—真实面积和真实周长。

PCAS 主要由三个参数确定：阈值（T）、孔喉封闭半径（r）、最小孔隙面积（S0）。本节主要通过介绍这三个参数的选取意义，分析PCAS的原理与方法。

1.1 阈值（T）

首先获取崩塌堆积物的正射影像图，圈定分析的范围，获得适当的图片大小，即开始图像识别。在图像识别前，需要进行二值化处理，获得二值图像，PCAS 采用阈值分割法进行图像二值化[13]，所谓阈值分割法，是基于影像图中各个组成要素的灰度值的不同，通过选择适当的灰度阈值将影像图转换为二值图像，根据图像中孔隙与颗粒的灰度值不同加以区分[14]。具体来说，图片由像素组成，而每一个像素是由红绿蓝（RGB）三原色分量组成，软件自动选取一个特征颜色值作为灰度值X（RGB），计算图像中各个部分的灰度值P（RGB）与这个特征灰度值的距离[12]：

与选取的阈值作对比，d小于阈值为黑色，识别为孔隙；大于阈值为白色，识别为颗粒；从而把图片识别成灰度图像（图2b）。

具体步骤为：起初设定较大的阈值，此时较多为颗粒的单元被识别为孔隙（黑色），再慢慢地减小阈值，直到黑色的细小颗粒单元被转化为白色，能清晰地识别出颗粒与孔隙，此时为较好的阈值。为了减少人为判断误差，通常对不同堆积区进行二值化处理，求取阈值平均值，作为选定的最终阈值。

1.2 孔喉封闭半径（r）

颗粒与颗粒间并不是完全分开，会有连接（图2b），在二值化过程中，颗粒通过细小连接，会被错误识别为一个颗粒。因此，PCAS 在识别过程中，会设定一个特定直径的孔喉封闭半径，定义为腐蚀结构元素的半径r，对图像做腐蚀运算[15]。当颗粒之间连接的直径小于2r时，则区分为两个独立区域，分开颗粒（图2c），再将剩余像素归并到颗粒上，实现颗粒与颗粒的自动区分[16]。

图2 PCAS 在T/r/S0=170/1/10时的颗粒识别结果Fig.2 PCAS particle recognition results at T/r/S0=170/1/10

1.3 最小孔隙面积（S0）

图片是由一系列像素点组成，当像素较低时，无法真实表征颗粒的真实形状，需要将这些像素点及噪声去除。因此，PCAS 识别程序中设定了可识别分析的最小孔隙面积，即设定一个面积识别下限，小于该值的颗粒面积都识别为杂点而被去除。在图2（b）可以看到二值图像中呈现为白色的细小杂点，设置最小孔隙面积，可以消除这些杂点（图2c），有效识别出颗粒。

PCAS 具有手动编辑的功能，对于颗粒不完整或较多明显杂点，可以通过界面的“Edit binary image”功能，手动修复，归并颗粒，区分出颗粒边界，减小系统误差。

2 PCAS 最优参数探索

以贵州纳雍普洒村崩塌-碎屑流为例，运用PCAS对其粒度进行分析，结合彭双麒等[2]的实测结果，研究PCAS 在崩塌碎屑流颗粒分析的应用，以及阈值、孔喉封闭半径和最小孔隙面积的参数选取。

2.1 研究区概况

彭双麒等[2]通过传统粒径统计方法，沿纳雍崩塌中部区域的主滑方向，将堆积区分为A1—A13 共13个区，每个区长100 m、宽50 m（图3）。通过生成地表数字模型（DSM），人工观测统计，得出每个区域颗粒粒径及对应数量。为了方便测量，粒径设定为颗粒对角线最大长度取值。

图3 崩塌堆积分区示意图[2]Fig.3 Partition diagram of the collapse accumulation body[2]

分别统计出每个区域0～1 m、1～2 m、2～3 m、3～4 m、4～5 m、5～6 m、6～7 m、7～8 m、8～9 m、9～10 m、10～17 m的粒径数量（识别最大粒径16.6 m），取值区间左开右闭。根据统计结果，绘制成累计级配曲线（图4），纵轴为小于某粒径颗粒的数量百分比，横轴为各粒径。

图4 分区颗粒粒径累计曲线[2]Fig.4 Aggregate curve of particle size in subsection[2]

为了更直观显示A1—A13 每个统计区域内大粒径与小粒径含量之间的关系，分别作出各区域0～2 m、2～4 m、4～6 m、大于6 m颗粒所占比例柱状图（图5）。

图5 各区域不同粒径等级所占比例Fig.5 Proportion of different grain sizes in different regions

2.2 PCAS 参数讨论

运用PCAS 对每个区域颗粒进行定量分析。不断调整阈值、孔喉封闭半径、最小孔隙面积，分别统计不同参数下对应的颗粒分布情况，找到合适的参数取值，分析不同参数条件下对研究结果的影响，验证PCAS 软件的准确性。

（1）PCAS 参数选取

导入图像，不断调整阈值，对13个区域用PCAS软件进行二值化处理。通过实际分析本次纳雍普洒村崩塌的图像数据，得出阈值在170时可以清晰地区分出颗粒与孔隙（图6）。

探究孔喉封闭半径（r）、最小孔隙面积（S0）的取值对PCAS颗粒识别的影响。不断调整孔喉封闭半径、最小孔隙面积合适的取值，并进行相应的组合。如图7所示，当r/S0=0.5/5（这里为像素，通过与实际半径、面积换算，r/S0=0.5/5=0.05 m/0.05 m2，即r单位为0.1 m，S0单位为0.01 m2）时，许多杂点及不需要的微小颗粒被识别，且颗粒与颗粒间不能有效分隔开；当r/S0=4/40时，较大颗粒也小于最小孔隙面积，不能识别出来；对比图2（c），当r/S0=1/10时，识别所得颗粒分布结果更能反映堆积物颗粒各粒径分布情况。因此，选取孔喉封闭半径r为1，2，3（单位为0.1 m，下同），最小孔隙面积S0为10，20，30，40（单位为0.01 m2，下同），互相组合（表1），找到合适的参数取值（各参数均以像素为单位）。

（2）PCAS 粒径统计结果

图6 PCAS 在不同阈值下的二值化结果Fig.6 Binarization results of PCAS at different thresholds

图7 PCAS 在r/S0=0.5/5 与r/S0=4/40时的颗粒识别结果Fig.7 PCAS particle recognition results at r/S0=0.5/5 and r/S0=4/40

表1 三个参数的组合情况Table1 Combinations of the three parameters

由各个图像中颗粒的几何形态数据，统计出每个区域长轴在0～1 m、1～2 m、2～3 m、3～4 m、4～5 m、5～6 m、6～7 m、7～8 m、8～9 m、9～10 m、＞10 m的粒径数量（识别最大粒径11.8 m），取值区间为左开右闭，绘制出颗粒粒径累计曲线。并分别作出各区0～2 m、2～4 m、4～6 m、＞6 m 粒径所占比例柱状图。注意，PCAS 识别的是像素大小，每个区域图片像素为465×976，而实际尺寸为50 m ×100 m，由式（1）和式（2）换算，为方便统计，取1 像素等于0.1 m。

分别根据表1中的参数组合绘制颗粒粒径累曲线（图8）和各区域不同粒径等级所占比例柱状图（图9），图9中粒径百分比小于1%时没有标记。

由图8可知：

（1）A11、A12、A13 图像中粒径小于2 m的百分含量明显高于其他区域，说明运动距离越远，颗粒碰撞破碎越充分，小颗粒含量比重越大，分选越好。在阈值为170，孔喉封闭半径为2时，A1 图像中粒径小于2 m的百分含量也明显高于其他区域，曲线较陡，颗粒粒径级配差于其他区域。

（2）当孔喉封闭半径为1，2时，随着最小孔隙面积的变大，颗粒粒径小于2 m的百分比明显下降，但总的趋势不变；表明PCAS 在最小孔隙面积增大时，能去除较大的杂点，区分出有效颗粒，总的百分比就降低。但当孔喉封闭半径为3时，随着孔隙面积的变大，大粒径百分比无明显变化，小颗粒粒径百分比变小，小于1 m的粒径百分比大多在10%以下（除A11、A12、A13）；这表明当孔喉封闭半径增加到一定程度时，许多颗粒被进一步分割成更微小颗粒，当最小孔隙面积较小时（这里为10），大量微细颗粒就被当作杂点去除，此时即使最小孔隙面积变大，粒径百分比也无明显变化。

图8 不同参数取值下各区域颗粒粒径累计曲线Fig.8 Accumulation curve of particle size in each region under different parameters

（3）以N2、N6、N10为例，当最小孔隙面积一定时，随着孔喉封闭半径的增大，颗粒粒径小于2 m的百分比明显下降，且级配变差。N10 条件下，每个区域粒径小于2 m的百分比分布在50%～90%，而N2、N6 小于2 m的百分比则在80%～95%；这表明随着孔喉封闭半径增大，形状不规则的大颗粒被腐蚀，分割为小颗粒，当小于最小孔隙面积时，被作为杂点去除，所占比重减小。

（4）小于3 m的颗粒百分比在每个区域的占比都在80%以上，而图4小于3 m的为70%以上，说明PCAS 相比人工计数，颗粒计数更精细，识别更完整，但总的效果相差不大，粒径级配曲线与实际情况相符，验证了PCAS 在碎屑流堆积物颗粒识别上的可靠性。

由图9可知：

（1）由图N1、N2、N5、N6可得，当孔喉封闭半径为1，2，最小孔隙面积为10，20时，各区域0～2 m 小颗粒所占百分比都在80%以上，大于4 m颗粒所占百分比都在4%以内；表明当孔喉封闭半径与最小孔隙面积较小时，各区域较多的小颗粒与杂点被识别，且不能被错误识别为单个颗粒。对比N3、N4、N7、N8，当孔喉封闭半径或最小孔隙面积增大时，0～2 m 小颗粒占比降低，大于4 m的大颗粒比重增加。

（2）由图N9、N10、N11、N12可得，当孔喉封闭半径为3时，各区域小于2 m的颗粒占比都具有先减小后增大的趋势，拐点在A3区；有3个明显波峰，极点在A2、A9、A12区，有三个明显波谷，极点在A3、A8、A10区。各区域4～6 m颗粒占比都在2%～4%附近，小于6 m的颗粒占比具有2个波谷，极点在A4 与A8区。各粒径占比分布情况与图5对比，基本吻合，接近真实值。

（3）对比图N9、N10、N11、N12 与图5，大于6 m颗粒占比在A3、A4、A5 三个区域与实际情况相符。图5中，0～2 m颗粒占比在A3、A4、A5区呈现波谷形态，但在PCAS 粒径识别占比中，A1～A10区0～2 m颗粒占比都在52%～69%，表明在PCAS颗粒识别系统中，小粒径颗粒识别效率高，各区占比较接近。

2.3 PCAS 优劣及可行性讨论

根据PCAS 在阈值为170、孔喉封闭半径为3、最小孔隙面积为30 三个参数下取得的各区粒径分布结果，结合人工统计数据，作各区粒径块数对比柱状图（图10）。结合图8、图9可得如下结果：

图10 PCAS 统计方法在T/r/S0=170/3/30时与人工统计各区颗粒块数对比图Fig.10 Comparison chart of the PCAS statistical method at T/r/S0=170/3/30 and artificial statistics

2.3.1 优势

(1) 高效性：PCAS 统计的细颗粒数量相比人工统计的多，但颗粒数量变化趋势基本一致，各区数量变化情况基本对应；0～2 m 粒径块数在A2、A4、A7 附近出现极值，2～4 m 粒径数量在A4、A7 附近出现极值。表明PCAS可以有效识别肉眼难以计数的细小颗粒，可以极大提高对小粒径颗粒（小于4 m）的识别效率，相比人工细小颗粒统计更高效。

(2) 准确性：当颗粒粒径大于4 m时，人工统计与PCAS 统计结果基本吻合（A8区异常）。人工对大颗粒的识别相比小颗粒准确性较高，这也证实了PCAS对堆积体识别的准确性，人工可以识别的，PCAS 也能自动准确识别。

(3) 规律性：从各区颗粒数量变化趋势可以看出，PCAS 统计颗粒数量变化呈相当好的规律，曲线圆润平滑，“波峰波谷”的数量极值与贵州纳雍等碎屑流崩塌实际情况对应较好[17-18]，且符合碎屑流堆积体的一般规律[19]。而人工统计由于肉眼的局限性，在颗粒数量、颗粒大小上总会存在判断失误，与实际情况会有一定的差异。

(4)可操作性：PCAS颗粒识别系统操作简单，只需设定阈值（T）、孔喉封闭半径（r）、最小孔隙面积

（S0）三个参数，即可得到堆积体颗粒分布情况，可节省大量人力物力。

2.3.2 劣势

PCAS 劣势主要体现在三个参数自身的特性上，分析如下：

(1)A2、A8区颗粒数量出现奇点；图10（a）、（b）中，A2区人工统计粒径小于4 m颗粒数量比PCAS 计数高，主要是因为当颗粒小于PCAS 设置的最小孔隙面积时，就被自动删除，导致A2区粒径小于4 m的颗粒数量较低；图10（c）、（d）中，A8区人工统计粒径4 m以上颗粒数量较高，主要是因为PCAS 设置的孔喉封闭半径过大，许多形状不规则的大颗粒被错误分割成多个小颗粒，导致大颗粒数量变少，并且这些误分割成的小颗粒一部分大于最小孔隙面积，被识别成单个颗粒，导致A8区2～4 m 粒径颗粒数量较人工统计的多。

(2)PCAS 图像二值化，是根据图像中孔隙与颗粒的灰度值不同加以区分的，主要取决于所得堆积体图像。如果堆积体图像的分辨率、对比度、饱和度以及拍摄角度等不同时，选取的阈值可能会有差异，图像二值化结果也可能不同。

2.3.3 可行性

(1) 综合图8～10，当PCAS颗粒识别系统中采用不同的参数取值时，对分析崩塌碎屑流颗粒分布情况会有一定的影响。当所得图像与本文纳雍崩塌图像相似，阈值170、孔喉封闭半径3、最小孔隙面积30为最优参数选择，各区不同粒径数量与人工统计所得结果基本对应，粒径累计曲线与各区域不同粒径等级所占比例最接近实际情况；反之，可根据本文所介绍的PCAS 识别方法、步骤以及参数选取的判定依据，依据具体情况以“170、3、30”为基准进行合理的参数调整。

(2) PCAS颗粒识别系统已在贵州纳雍崩塌[10]和金沙江白格滑坡[11]中得到应用，并且效果显著。得到滑坡堆积体图像后，只需设定三个参数，即可对堆积体进行颗粒识别分析，操作简单，并且PCAS 具有手动编辑的功能，可以通过界面的“Edit binary image”功能，根据实际情况进行颗粒的调整。在新技术不断发展的今天，PCAS 不失为崩塌碎屑流粒径研究的一个高效、可行的途径，具有应用价值。

3 结论

（1）PCAS 能自动准确地识别碎屑流堆积物颗粒与孔隙，具有高效、准确、可操作性强等特点，可用于堆积物颗粒数字图像的识别、量化、分析；确定三个关键参数：阈值（T）、孔喉封闭半径（r）、最小孔隙面积（S0）的合适取值后，即可得出粒径分布情况，节省大量人力物力。

（2）针对颗粒粒径较大、二值化后孔隙与颗粒区分不明显、噪点较多的图像，采用较大的r/S0值更能反应颗粒分布宏观情况；反之，宜适当减小r/S0进行分析。当所得图像与本文纳雍崩塌图像相似，阈值170、孔喉封闭半径3、最小孔隙面积30为最优的参数选择，可得到有效的碎屑流颗粒分布情况；反之，可以“170、3、30”为基准进行合理的参数调整。

（3）PCAS 具有较高的可行性，所得贵州纳雍崩塌堆积物粒径识别和统计结果与实测值接近，粒径占比、分布规律基本吻合，即堆积体小粒径占比较多，0～2 m颗粒粒径各区占比都在50%以上，且运动距离越远小粒径比例越大，各区“波峰波谷”变化趋势基本对应。