崩积体粒径的图像识别与分析

2021-10-08徐今星杨根兰史文兵江兴元

科学技术与工程 2021年26期

徐今星，杨根兰，梁风，史文兵，江兴元

(贵州大学资源与环境工程学院/喀斯特地质资源与环境教育部重点实验室，贵阳 550025)

崩塌是一种常见的全球性泛生型山地灾害，其具有隐蔽性高、失稳破坏迅速、致灾后果严重的特点[1]。除崩塌体下落直接撞击带来危害外，堆积于高陡斜坡上崩塌堆积体也极易诱发二次崩滑灾害[2]，严重威胁人们的生命财产安全。

为减少崩塌堆积体产生的二次灾害，中外学者开展了众多研究。一部分学者从源头着手，研究崩塌堆积体成因机制[3-4]、运动特征[5]、落石范围[6]、冲击力[7]等。另一部分学者着重研究崩塌堆积体本身，关注崩塌堆积体的堆积特征[8-11]、地震[12]、降雨[13]工况下的稳定性、以及失稳后的变形破坏模式[14]等。

中外学者虽已对崩塌堆积体开展了众多研究，但目前针对崩塌堆积体块石粒径分布特征研究尚显不足。粒径分布特征直接影响堆积体内粒径迁移，最终影响堆积体稳定性，是崩塌堆积体重要的堆积特征。目前研究堆积体粒径主要依靠传统野外地质调查，需要地质工作者人工实地测量[15-16]。但山区地带自然地理条件复杂，对于处于高陡悬崖上的堆积体，大多难以到达，传统的地质调查工作难以开展。因此，如何不受地形条件制约、准确高效地识别、统计堆积体粒径已经成为迫切需要解决的问题。

因此，现以小茅坡崩塌堆积体为例，提出基于无人机航拍影像技术以及数字图像处理技术的堆积体粒径识别、统计方法。即通过无人机航拍影像三维重建堆积体；生成高分辨率数字地表模型(DSM)和生成数字正射影像(DOM)，通过孔隙(颗粒)与裂隙图像识别与分析系统[pores(particles) and cracks analysis system,PCAS]统计计算，获取堆积物表面信息，进而对堆积物粒度分布与堆积特点进行研究。研究结果对堆积体稳定性分析调查具有重要意义，为崩塌堆积体野外调查提供一种新的技术手段。

1 研究区域概况

小茅坡崩塌位于贵阳市开阳县金钟镇小茅坡村，地理坐标为106°48′47.01″E，27°4′2.79″N。研究区位于洋水背斜左翼，为侵蚀、剥蚀的低中山地貌类型，呈明显单面山的陡-缓-陡剥蚀台地地貌，研究区地质图如图1所示。单面山剥蚀台地的延伸方向近南北向，与背斜轴方向一致。危岩体所在山体的斜坡顶部高程1 400～1 450 m，高差300～350 m，坡向107°、坡度30°～51°。区内植被较发育，多为乔木和灌木，小茅坡全貌崩塌如图2所示。岩层主要为震旦系上统灯影组浅灰色厚层块状白云岩，岩层产状119°∠12°，主要发育两组共轭节理面，J1: 150°～170°∠60°～70°J2: 210°～220°∠65°～70°。

图1 崩塌区地质图Fig.1 Geological map of collapsed area

图2 小茅坡崩塌堆积区全貌图Fig.2 Overall view of Xiaomaopo collapsed accumulation area

2 无人机航空摄影测量工作方法

无人机倾斜摄影测量技术是以无人机为飞行平台搭载传感器设备，来获取地面信息的遥感方式。可以通过规划无人机航拍倾斜摄影，进而获得包含拍摄点地理位置信息的高清照片。

崩积体沿坡面堆积，而坡面倾斜且高差较大。为了提高获得正射影像图的分辨率，本次航拍采用大疆精灵4 RTK无人机最新的斜面航线技术(图3)。通过斜面航线功能，能够自动生成针对斜面或者立面场景的飞行航线，实现斜面或立面摄影测量的数据采集。斜面航线飞行操作流程如下：①飞手采用GS RTK App手动飞行打三个航点以定义被测目标平面；②设置外扩距离生成默认的平行四边形测区；③根据地面分辨率的需求调整任务高度。

图3 斜面航拍示意图Fig.3 Schematic diagram of aerial photography of slope

2.1 无人机航拍路线规划

航拍路线规划要充分考虑研究区地形地貌特征，确保无人机航拍安全以及相片重叠率。由于陡崖面、堆积体斜面高差较大且部分堆积体粒径较小，为最大程度还原崩塌危岩体与堆积体特征，航拍规划为两架次。第一架采用斜面低空航拍飞行拍摄堆积体;航向重叠率不低于80%，旁向重叠率不低于80%;规划4条航线，航线间距为15 m，航拍高度根据预设斜面变化。第二架次手动控制拍摄整个崩塌区域，保证崩塌危岩体、陡崖面以及堆积体细节信息获取充足。第一架次、第二架次分别获取781、400张照片。

2.2 高清数字地表模型建立

采用Agisoft PhotoScan对带有地理位置信息的航拍影像照片进行拼接处理，生成高清数字地表模型(DSM)，与高清正射影像图(DOM)。所生成的高清DSM可以测量崩塌危岩体结构面产状、人工测量并统计堆积体粒径大小，具有安全、高效、便捷等特点。Agisoft PhotoScan通过导入挑选好的照片、对齐照片可以自动提取稀疏点云，并依次加密点云、生成空间白模、网格、纹理，输出正射影像图。

2.3 精度对比与验证

分辨率是指图像中一个像素所代表的实际地物的面积。分辨率是后文识别图像中岩块实际几何尺寸的依据。虽然通过三维重建后建模软件会自动计算图像分辨率，但有时有较大误差，因此需要设置参照物。为了更好地区别出参照物，使用与岩块、环境颜色差异明显红色油漆对参照物表面进行涂抹，并测量参照物的实际长度与宽度，如图4所示。

图4 参照物图Fig.4 Reference object map

3 PCAS图像处理方法

PCAS是由南京大学开发的一种通过数字图像处理，常用于微观颗粒、孔隙以及裂隙的识别分析的专业软件[17]，具有高效便捷的工作特点，用于岩块粒度识别有较好的效果。软件通过将图像二值化、矢量化处理，可以计算颗粒、孔隙以及裂隙的数目、长宽、方向、形状系数等统计参数[18-19]。PCAS软件处理流程主要包括以下几个部分：①正射影像图获取；②图像预处理；③图像二值化处理；④二值图加工处理；⑤堆积体粒径识别与统计。PCAS图像处理工作流程图如图5所示。

图5 PCAS图像处理过程Fig.5 PCAS image processing process

(1)正射影像图获取。采用2.2节中Agisoft PhotoScan拼接的整体崩塌区域高清正射影像图。

(2)图像预处理。整体崩塌区域高清正射影像图像素过大，无法直接导入PCAS进行处理，需按照分区图分别裁剪成10小段方能导入。此外，为了保证良好的识别效果，应根据崩塌堆积体分布范围通过二次裁剪尽可能将细小粒径与大粒径区分开，如图6(a)所示。此外，当细小块石边界模糊时，可采用Photoshop滤镜的UnsharpenMask(USM)锐化进行锐化处理，选中模糊边界通过设置锐化半径、数量以及阈值，使边界变得更加清晰。

图6(a)中白色区域为PCAS识别所得块石，黑色为块石之间的空隙；图6(d)中彩色区域为PCAS识别所得块石，黑色为空隙图6 PCAS图像处理过程Fig.6 PCAS image processing process

(3)图像二值化处理。正射影像图为真彩色图像，其利用3个大小相同的二维数组R(红色)、G(绿色)、B(蓝色)表示一个像素，其中每个数组的范围为0～255，通过这三种颜色可以合成任意颜色。由于堆积体块石与堆积环境以及块石间孔隙颜色具有明显差异，二值化处理时以某一需要识别的崩积块石的某一点像素点的RGB值(P)作为基础，图像中其余像素点RGB值(X)与其有差值(D)，处理前需设置一个合适的阈值(T)，当D>T时，像素被识别为空隙，当D

D=|X-P|

(1)

(4)二值图加工处理。由于堆积体块石不会完全是同一颜色，外加环境中植被等因素的干扰，二值图中往往会出现岩石识别不全，导致出现空洞或者边缘识别不清，如图6(b)所示，图6(b)中红色圈出的块石出现了识别不全的现象。此外，由于块石粒径过小或者块石堆叠导致的块石边界模糊，PCAS软件识别时，会因块石边界连接错误地将多个小块石识别为一个大的块石。为解决此类问题，可以手动处理二值图，即删除像素点块石边界清晰，或增加像素点填补块石边界及块石中的空洞。图6(c)为人工处理后的二值图。

堆积体块石识别与统计。采用PCAS软件识别块石时，基于种子算法，通过设置最小识别面积阈值，除去像素小于阈值的细小噪声，并识别像素大于面积阈值块石；基于腐蚀运算，通过设置封闭半径(element radius)，腐蚀掉块石重叠部分后，通过种子算法识别独立块石，并将腐蚀掉的像素归还到原块石中得到块石参数。识别完毕后，可得到块石的数目、长宽、方向、形状系数等统计参数。其中，块石的长、宽分别表示最大和最小费雷特直径。费雷特直径(Feret diameter，dF)定义为相对于特征的一对平行切线与特征之间的正交距离，即块石长轴投影到统一坐标方向上的最大距离。

4 崩塌体块石统计结果

4.1 PCAS统计

由于整体堆积体区域高清正射影像图像素过大，无法直接导入PCAS进行处理，为了更好地统计堆积体粒度组成，需将主要堆积体区进行分区。沿主崩塌方向将堆积体由崩塌后缘的堆积区开始，每20 m划分为一个的矩形区域，共10个区域，编号从后缘开始，依次为A1～A10，主要堆积区分区如图7所示。将分区后的图像裁剪出来分别处理，取块石的长为块石的直径，得到每一个区域的块石粒径的像素大小。根据设置最小识别面积阈值50像素、封闭半径2像素可知，PCAS软件可以识别的最小粒径为0.2 m。将长度小于1 m块石每隔0.3 m划分为一个统计范围，将大于等于1 m块石每隔1.0 m划分为一个统计范围，共划分为12个统计范围。统计时，取左开右闭区间，PCAS统计结果如表1所示。

A1～A10、B1～B3、C1～C3为分区编号图7 主要堆积体分区图Fig.7 Main integral partition figure

表1 PCAS纵向统计结果表Table 1 PCAS longitudinal statistical results table

由于崩塌后落石运动的轨迹与坡面有关，为了更好地将崩积体的崩积特点与崩塌的运动过程联系起来，选取堆积主轴方向切剖面图，如图8所示。

Zbdn为震旦系上统灯影组图8 堆积方向剖面图Fig.8 Stacking direction profile figure

4.2 数字地表模型块石人工统计

考虑到A10区整体粒径偏小，仅统计A1～A9可见块石粒径。由于崩积体中上部高陡调查人员无法前往，块石人工统计分为模型统计与野外统计两个部分。

(1)野外统计：调查人员可以到达的区域采用野外实地测量的方法进行统计，统计结果如表2所示。

表2 野外纵向统计结果表Table 2 Field longitudinal statistical results table

(2)模型统计:参照彭双麟等[20]，采用对DSM图中崩塌堆积岩块的人工测量统计方法进行验证。统计较多的A4、A5、A6区数量在600块左右，而统计较少的A1、A9区数量在100块左右，其余分区统计数量在100～200块。由于肉眼无法辨别细小块石，所以模型测量统计总量为PCAS统计的1/10～1/50，模型统计结果如表3所示。

表3 模型纵向统计结果表Table 3 Model longitudinal statistical result table

5 统计结果验证与分析

5.1 PCAS统计可行性验证

三种方法统计的结果总量上有较大差异，因为统计的是同一个区域，所以两种统计的结果应服从相同的分布规律。据此，为了验证方法的可行性，对三种统计结果采用IBM SPSS Statistics软件进行回归曲线估算。所有估算结果表明方程为Y=A+B/t的逆模型决定系数R2最高，表明统计结果最符合逆分布。

根据SPSS估算结果，采用方程为Y=A+B/t的逆模型在origin中拟合两种统计结果，如图9所示。野外统计R2分别为0.87、0.89、0.89、0.88；模型统计R2分别为0.92、0.90、0.91、0.91、0.90、0.91、0.91、0.94、0.96；PCAS统计R2分别为0.84、0.87、0.89、0.81、0.81、0.82、0.82、0.84、0.85。由于所有区决定系数R2均大于0.8，可以认为统计结果均服从逆分布，即PCAS统计结果真实可靠。A1～A3区模型统计与PCAS统计拟合曲线高度相似，而由于A4～A9区统计总量较大，小颗粒较多，导致模型统计与PCAS统计的拟合曲线有一定差异；A6～A7野外统计与模型统计拟合曲线高度相似，而A8～A9区野外统计与PCAS统计拟合曲线高度相似。从整体上看，PCAS统计数量>野外统计>模型统计。

图9 PCAS统计与人工统计拟合结果图Fig.9 PCAS statistics and manual statistics fitting results

5.2 崩积体的纵向分布规律

为了能够直接对比分析两种统计结果，并减少由于统计总量差异导致的误差，将两种统计结果分别制成粒径累计曲线，如图10和图11所示。

图10 人工统计粒径累计曲线Fig.10 The cumulative curve of artificial particle size

图11 PCAS粒径累计曲线Fig.11 Cumulative curve of PCAS particle size

为了直观地了解各个分区粒径的分布特点，将PCAS统计结果按照分区制成百分比堆积柱状图，如图12所示。崩积体的纵向分布规律如下。

由图10和图11可知：总体上来看，PCAS统计与人工统计的每个区域累计曲线所表达出累计特点、趋势基本相同。由于PCAS统计的数据是人工统计的10～50倍，而且可以统计肉眼无法分辨的小块石，PCAS累计曲线连续、较为光滑可靠，并且2 m及以下的区域与区域之间差异更加明显。由人工统计小于1.0 m块石所占的含量为68%、PCAS统计小于1.0 m块石所占的含量为81%可知，小块石为堆积体的主要组成部分，表明危岩体崩落后碎裂化程度较大。

图12 各区不同范围粒径所占的百分比柱状图Fig.12 Histogram of the percentage of particle size in different areas in each area

(1)粒径小于2.0 m块石。粒径小于2.0 m块石所占的百分比含量都呈现先减小后增大在减小的趋势，这与通常情况不同。PCAS统计结果显示小块石在主要堆积区的中部A3～A7区、下部A9区含量较高，所有区小块石百分比含量均超过95.0%，且各区间差异小于3%。产生这种现象的原因有两个：①A1～A3区小于0.2 m的细小块石占一定的比例(图7)，但由于测量精度不足导致的统计缺失；②崩塌发生时，落石势能转化为动能，落石在坡面滚动、弹跳，坡面突然发生改变，由A1～A3的陡倾坡面转为A4～A7较缓坡面，使落石在坡面运动过程中相互碰撞、撞击坡面，这导致落石碎裂化加大，进而导致更多小块石的产生。

(2)粒径在2.0～3.0 m块石。粒径在2.0～3.0 m块石所占的百分比含量出现了前大后小两个峰值的双峰情况，变化曲线为一不规则M形，与郑光在碎屑流中发现规律一致[21]，呈现先增大后减小再增大最后在减小的趋势。第一峰值在崩源区的A1～A2区附近，第二个峰值同时位于远离崩源区的A8区，这表明2.0～3.0 m粒径中等的块石主要分布于A1～A2与A8区。

(3)粒径大于3.0 m块石。粒径大于3.0 m块石所占的百分比含量块石所占的百分比含量都呈现先减小后增大在减小的趋势，第一峰值在崩源区的A1区附近，第二个峰值同时位于远离崩源区的A8区，即大块石主要分布于A1以及A8区。这表明先崩落的大块石，由于能量较大，运动距离较远，但与A1～A6区相比，A7～A9区坡度较缓，大块石受到坡面的阻力较大，大块石能量被迅速耗散，最终停积于此；而后崩落的大块石，由于坡面先崩落块石的缓冲作用以及植被的阻挡，并没有发生滚动，而是直接积于此。

5.3 崩积体的横向分布规律

为研究堆积体粒径的横向分布规律，将位于堆积体前、中、后部位的A3、A6、A7三个区分解为B1～B3、C1～C3、D1～D3。采用PCAS分别将以上各区粒径进行统计，统计结果如表4所示。

表4 PCAS横向统计结果表Table 4 PCAS horizontal statistical result table

由图13、图14可知，距离崩塌源区距离相同的横向各区粒径组分的含量基本相同，中粒径块石与大粒径块石在横向上呈现先增大后减小的趋势，这与彭双麟等[20]对崩塌-碎屑统计结果相似。其中，位于堆积线上的各区(B2、C2、D2)中粒径块石与大粒径块石所占比例最大。结合现场调查，原始坡面全部为相同的残坡积土，无耕植土层且坡面横向变化较小。崩塌块石崩落后，块石产生弹跳、滚动，沿142°堆积，并向两侧扩散。因此，位于堆积线上的各区块石数目最多，且大、中粒径块石含量最大。此外，随着落石运动距离的增加，横向各区粒径差异在不断增大。与5.2节相同，横向各区小粒块石呈现先减小后增大在减小的趋势，其他粒径则与之相反。

图13 PCAS横向粒径累计曲线Fig.13 Horizontal cumulative curve of PCAS particle

图14 横向各区不同范围粒径所占的百分比柱状图Fig.14 Horizontal histogram of the percentage of particle size in different areas in each area

5.4 崩积体特征

崩积体以白、灰白色白云岩块体为主，呈条带状向冲沟内堆积，处于基本稳定状态，堆积方向142°。崩积体上部靠近崩塌源处高程1 290 m，下部最远端高程1 010 m。该崩积体纵向长度约480 m，横向平均宽度约86 m(图13)。崩积体下部原地表植被主要为低矮灌木，崩塌落石自崩塌源而下，沿崩落路径翻滚、跳跃、滑动，在地表形成崩落坑、倒石堆以及刮擦痕迹，原坡面被改造，植被遭到破坏，尤其是近崩塌源的陡坡处，虽然坡面刨蚀深度较小，但植被完全被破坏。根据统计区域的数目特征以及野外调查结果，将崩积体划分为3个区，A1～A7为主积区，A8～A10为散落区，A10下侧至最远滚落的崩积边缘为影响区。

(1)主积区：高程1 230～1 290 m，为落石的主要堆积区，堆积块石共27 186个，占全部的95.81%，区内刮铲现象明显、大部分植被被破坏，破坏面积占区内总面积的90%左右。上部A1～A3区坡度达到40°～51°，平均坡度45°。区内崩积体粒径主要分布范围在0.2～2 m，小块石主要堆积在区内上部靠近崩塌源处的A1～A3区，A1～A3区中小块石占主积区内崩积体的27.43%，区内最大块石为8.5 m，位于A2区；自中部向下的A4～A7区落石开始大量堆积，A4～A7区坡度达到33°～44°，平均坡度37°。与其他区相比，区内中等块石数目较多共402个，占该区的2.05%，区内堆积崩积体粒径分布范围在0.2～4 m，区内最大块石为9.6 m，位于A4区。块石间为碎石填充，崩积体中出现架空现象。

(2)散落区：高程1 210～1 230 m，为少量落石散落堆积区，堆积块石共1 174个，占全部的4.14%，区内坡度为32°，少量植被被破坏，崩积体粒径主要分布范围在0.2～3 m，区内最大块石为7.4 m，位于A8区。散落区内主要影响简易公路，但由于人为扰动，目前简易公路上所有块石均被搬运至简易公路两侧。

(3)影响区：高程1 151～1 210 m，为落石运动最远的区域，该区范围较大，堆积块石共14个。区内坡度为32°，堆积崩积体粒径主要分布范围在0.2～4 m。影响区内主要影响简易公路。

6 讨论

通过对比分析人工与PCAS统计结果的异同，可知无人机倾斜摄影-PCAS统计的优劣势。

6.1 优势

(1)与传统人工统计相比，无人机倾斜摄影-PCAS统计不需要现场实地统计，因此具有便捷、安全、成本低、高效等工作特点。

(2)与传统人工统计相比，PCAS统计通过数字图像处理技术对图形中的岩石粒径大小进行直接识别，具有高精度性的特点。这使得PCAS在统计小粒径块石时远远优于人工统计。

(3)与传统人工统计相比，PCAS统计结果数据量大，不会出现漏记、错记、重复记录情况，这使得统计曲线更加平滑可靠。

6.2 劣势

(1)PCAS是通过图像处理识别统计块石粒径大小，因此其统计结果是二维的表观块石。与碎屑流相比，崩塌堆积体碎裂化并不彻底，实际的块石并不是长度、宽度、高度相近的几何体。由于采用的是无人机航拍正射影像图，当块石的高度远远大于长度与宽度时PCAS统计的结果会偏小。此外，当多个块石三维重叠时，重叠面积大于腐蚀运算的封闭半径的多个块石会被错误的识别为一个大的块石。

(2)由于PCAS采用二值化的方法减小计算量，但二值化过程中会产生噪声。为了去除细小的噪声，采用最小识别面积阈值的方法，这导致真实面积小于阈值的块石被一并当作噪声除去。

(3)PCAS识别的本质是识别块石颜色与周围环境的差异，当块石颜色与周围环境的差异较小时，识别效果较差。此外，环境中植被等因素的干扰较大时，二值图中往往会出现岩石识别不全，导致出现空洞或者边缘识别不清，需要人工干预量较大。