0.   引言

准确获取结构面信息是岩质边坡稳定性分析的重要前提。传统的结构面接触式测量方法有测线法和测窗法，这2种方法操作简单，但对测量者的专业技术水平要求较高，易受测量者主观影响，对于复杂结构面需要进行多次测量和综合分析，工作量大且效率低；尤其对于高陡边坡，测量者无法直接进行人工编录，作业危险性大，难以获取全面的结构面信息^[1-2]。因此，高效准确地采集岩质高边坡结构面信息仍是一个关键难题。

针对传统接触式测量技术存在的各种问题，研究者提出了非接触式测量技术，包括三维激光扫描技术和无人机摄影测量技术^[3]。三维激光扫描技术对测量环境要求较高且相关设备价格昂贵。无人机摄影测量技术以无人机为云台，通过机身搭载高精度数码相机和传感器对被测物体拍摄高清影像，然后结合计算机图像处理软件进行后处理与分析^[4]，具有机动灵活、性价比高、维护简单和适用性强等优势，因而逐渐被应用于数字化地形测绘、建筑物测量等领域^[5]。

近年来，研究者将无人机摄影测量技术应用于岩质高边坡结构面识别中。赵明宇等^[6]通过无人机和全站仪对边坡进行联合测量，获取了边坡结构面迹线和产状信息。熊开治等^[7]开展了基于无人机摄影测量技术的随机裂隙网络模型构建方法研究。李水清等^[8]基于无人机和“三点法”提出了一种岩体结构面产状统计的半自动方法，克服了传统方法的局限性。虽然无人机摄影测量技术可以获取岩质高边坡表面的高清影像，但是缺少高效准确的影像后处理方法，导致结构面识别仍然存在准确性差和效率低等缺点。一些学者尝试采用不同数据处理算法进行结构面识别。沙鹏等^[9]在采用无人机航测获取边坡图像的基础上，结合运动恢复结构法（Structure from Motion，SfM）获取边坡岩体点云模型，进而建立包含离散裂隙网络的实际边坡三维离散元模型。贾曙光等^[10]使用无人机采集北京市某采石场地形数据，结合计算机视觉原理与SfM方法重建了三维模型测量结构面产状。张恺等^[11]采用无人机和SfM法对浙江某露天采石场边坡进行三维重构，提取岩体结构面产状。Liu Yongqiang等^[12]采用无人机拍摄高分辨率影像，开发了一种同时考虑产状、迹长和起伏度的结构面分组新方法。然而，当前算法在构建高密度三维点云模型、精准识别与拟合复杂结构面存在局限性，且没有充分考虑识别结构面产状信息统计特征的不确定性，亟需构建高效准确的融合高边坡摄影、后处理算法及统计分析的一体化结构面识别与产状信息采集方法。

本文以江西省南昌市某露天矿高边坡为研究背景，提出一种岩质高边坡结构面识别及产状统计信息采集方法。采用无人机摄影测量技术和计算机视觉原理，使用无人机对边坡进行航测并获取基础影像数据；采用Context Capture系列图像处理软件、Cloud Compare、DSE（Discontinuity Set Extractor）和多种算法进行边坡结构面识别，进而提取结构面产状信息，进行产状统计特征分析。

1.   基于无人机摄影测量的结构面信息采集

1.1   研究区概况

露天矿边坡现场如图1所示。通过现场踏勘得知，边坡长约200 m，坡顶坡底落差约100 m，边坡出露岩体节理裂隙发育。将研究区岩体划分为上层和下层2个部分，每层岩体高陡，难以直接采用传统的接触式测量方法进行结构面产状测量。选择地点开阔、人员稀少且不受磁场干扰影响的地区作为无人机飞行起降点，并通过云台设置覆盖整个研究区的飞行路线，进行无人机摄影测量工作。

图  1  露天矿边坡现场

Figure  1.  Open pit slope site

下载: 全尺寸图片 幻灯片

1.2   无人机影像处理

考虑到研究区地形复杂、障碍物多等因素，选择体积小、控制性能好、操作简单的Phantom 4 Pro V2.0四旋翼无人机进行测量。无人机系统主要包括飞行器、云台、视觉系统、红外感知系统、相机、遥控器、智能飞行电池和存储器等组件。

采用无人机获取高分辨率边坡影像后，需要选择合适的计算机图像处理软件基于图像二维坐标构建三维点云模型。目前，常用的无人机航测影像后处理软件主要包括Pix4D、大疆智图、Photo Scan、Context Capture、Open Drone Map等。其中， Context Capture具有快速、简单、全自动等优势，能够在无人工干预情况下，从简单连续的二维影像还原出真实的三维实景模型和超高密度点云模型。此外，该软件还具有良好的输入数据兼容性和广泛的输出数据格式，能够满足绝大多数的工作需要。因此，本文选择Context Capture软件进行无人机拍摄影像的后处理。

由于无人机配备有惯性测量装置传感器及全球定位系统（Global Positioning System，GPS）与全球轨道导航卫星系统（Global Orbiting Navigation Satellite System，GLONASS），在飞行作业中拍摄的影像均保留了传感器信息、相机参数、照片姿态参数、照片位置参数等元数据，从而保证能够输出色彩真实、几何形态及细节精细的三维模型。Context Capture构建三维点云模型的流程如图2所示。生成的边坡上层和下层岩体高密度三维点云模型分别如图3和图4所示，可看出模型与边坡实物图吻合度较高。

图  2  三维点云模型构建流程

Figure  2.  Construction process of 3D point cloud model

下载: 全尺寸图片 幻灯片

图  3  上层岩体三维点云模型

Figure  3.  3D point cloud model of rock mass at the upper part of the slope

下载: 全尺寸图片 幻灯片

图  4  下层岩体三维点云模型

Figure  4.  3D point cloud model of rock mass at the lower part of the slope

下载: 全尺寸图片 幻灯片

2.   岩质高边坡结构面识别

2.1   识别流程

在完成原始点云数据重采样后，采用A. J. Riquelme等^[13]基于Matlab编写的开源软件DSE进行结构面识别，流程如图5所示。

图  5  结构面识别流程

Figure  5.  Identification process of structural plane

下载: 全尺寸图片 幻灯片

1） 将带有空间点坐标（X，Y，Z）和颜色信息的重采样点云数据输入到结构面提取程序中。

2） 设置K近邻（K-Nearest Neighbor，KNN）算法中需要的参数−近邻点数量K和偏差参数${\eta _{\max }}$，并根据确定好的参数寻找共面点。

3） 在搜寻到所有共面点子集后，分别进行最优平面拟合，计算拟合平面的法向量并作为对应目标点的法向量。

4） 在法向量密度图的基础上，采用Riquelme算法计算极点密度，并确定主极点。

5） 采用基于密度的聚类（Density-Based Spatial Clustering of Applications with Noise，DBSCAN）算法进行聚类分析和结构面识别，并通过最小二乘法计算结构面的平面方程。

6） 输出三维点云及结构面数据。

2.2   关键技术

2.2.1   点法向量计算和极点统计分析

在边坡结构面识别过程中，寻找特征向量是进行后续结构面分组与提取的前提。通常选择点法向量作为特征向量，但是每个单独的点无法计算法向量，需要利用每个点和其相似点构成平面的法向量作为该点的法向量。KNN是一种非参数、有监督的学习分类器，通过邻近度指标对单个数据点进行分类或预测^[14]。KNN算法主要包括近邻点距离度量法和确定近邻点数量法。① 近邻点距离度量法：为了确定最接近给定目标点的近邻点，需要计算该点与其他数据点之间的距离，以该点为球心，将所有满足给定距离的相邻点都归属于同一个点集。② 确定近邻点数量法：规定目标点需要寻找的近邻点数量为一定值，从目标点出发，将距离目标点最近的K个点归属于同一个点集。

由于点云密度存在不均匀性，采用近邻点距离度量法进行搜索时，可能会因为边界处包含错误点而形成错误平面，导致后续计算流程出现错误^[15]。而选择确定近邻点数量法进行搜索时，可通过适当缩小搜索半径降低结果出错率。因此，本文选择确定近邻点数量法构建相似点集。

需要指出的是，设置合适的近邻点数量K是保证KNN算法性能的关键。如果K值过小（如K < 15），那么在计算极点密度过程中会保留噪声；若K值过大（如K > 30），则会降低结果的准确性。同样，对于偏差参数${\eta _{\max }}$，如果${\eta _{\max }} < 15{\text{%}} $，会导致符合的共面点被丢弃；而当${\eta _{\max }} > 25{\text{%}} $时，则会产生边缘点。根据A. C. Rencher和W. F. Christensen^[16]的研究，如果1组主成分的方差贡献率不小于80%，则可认为其能够正确地表示数据。因此，本文中取K = 30，${\eta _{\max }}{\text{ = }}20{\text{%} }$。

确定点法向量后，采用Riquelme算法计算极点的密度函数。该方法为点云中每个点指定主方向，但由于曲面的读取误差和极点分散等问题，使得局部极大值极点中只有少数点才是主极点。为此，通过设置参数${n_{\mathrm{p}}}$和$\gamma $来寻找局部最大值点并作为主极点，${n_{\mathrm{p}}}$为可建立结构面组的最大数量，$\gamma $为点法向量与指定主极点法向量之间的最大允许角度。本文设置${n_{\mathrm{p}}}{\text{ = }}10$，$\gamma $取20～30°。

2.2.2   聚类分析与平面拟合

为了找到每个结构面组的数据子集并实现结构面识别，本文采用文献[17]提出的DBSCAN算法进行聚类分析。文献[18]证明了该算法在点云数据处理中的实用性。该算法需要设置2个输入参数：ε和 min-pts。ε 表示可视为相邻点的两点间最大距离，一般根据点云密度确定；min-pts为DBSCAN算法中的阈值参数，即在1个目标点的邻域内至少包含的点数。min-pts过大会导致碎片化，min-pts过小会导致子集连接过多。本文设置ε为4，min-pts为5。

聚类分析原理：在数据集中任选1个初始点，如果在距离该点ε 距离内有超过min-pts个点，则认为这些点都属于同一个聚类。以所有新添加的点为起始点，搜索距离该点ε 距离内是否也有超过min-pts个点。如果有超过min-pts个点，则把这些点也增加到与初始点同一聚类中，依此循环，直至没有符合条件的可添加点为止，便可得到1个聚类的所有子集。以该聚类外的任一点为初始点重复上述过程，直到所有的聚类分组都被找到。

完成聚类分析后，需要计算结构面的平面方程。设结构面组i的点集为${R_i}$，属于该组结构面聚类j的点集为${R_{ij}}$，寻找${R_{ij}}$的最佳拟合平面，其平面方程为

式中：（$x$,$y$,${\textit{z}}$）为点数据的坐标信息；${A_{ij}}$，${B_{ij}}$，${C_{ij}}$分别为平面单位法向量在x，y，z方向的分量；${D_{ij}}$为原点到平面的垂直距离。

2.2.3   点云数据前处理

经过Context Capture软件处理获得的点云密度通常很大，直接用于结构面识别需要大量计算，时间成本大，效率低。为此，对原始点云数据进行前处理，在保证计算结果准确性的前提下，最大限度地提高计算效率。Cloud Compare软件是一款开源的三维点云编辑与处理软件，支持许多点云处理算法（包括点云配准、重采样、统计计算、自动分割等），并包含大量插件可供选择。因此，本文选择Cloud Compare软件进行点云数据重采样处理。重采样方法包括空间距离重采样、随机重采样、八叉树重采样3种。综合比较，本文选择空间距离重采样方法。由于研究区域范围较大，原始点云数据过于庞大，为提高计算效率，在Cloud Compare软件中进行重采样时，设置两点间的最小距离为5 cm。

2.3   识别结果验证

为了验证边坡结构面识别结果的准确性，将边坡结构面识别结果与野外实地勘测的结构面信息进行对比分析。

图  6  上层岩体法向量密度

Figure  6.  Normal vector density of rock mass at the upper part of the slope

下载: 全尺寸图片 幻灯片

分别对图1中高边坡的上层岩体和下层岩体进行结构面识别，得到岩体法向量密度图和三维点云，如图6—图9所示。边坡上层岩体共找到6个主极点J₁—J₆（倾向/倾角）：J₁（156.26°/66.07°），J₂（177.52°/57.46°），J₃（147.54°/34.49°），J₄（197.45°/82.27°），J₅（127.88°/78.31°），J₆（221.05°/88.77°）。边坡下层岩体共找到5个主极点J₇—J₁₁：J₇（155.00°/43.01°），J₈（177.51°/57.45°），J₉（135.02°/75.48°），J₁₀（209.04°/72.65°），J₁₁（354.48°/73.11°）。

识别的三维点云数据与实际边坡结构面的对比如图10所示。可看出，用本文方法识别的结构面与边坡工程现场实际情况基本吻合。

图  7  下层岩体法向量密度

Figure  7.  Normal vector density of rock mass at the lower part of the slope

下载: 全尺寸图片 幻灯片

图  8  上层岩体三维点云

Figure  8.  3D point cloud of rock mass at the upper part of the slope

下载: 全尺寸图片 幻灯片

图  9  下层岩体三维点云

Figure  9.  3D point cloud map of rock mass at the lower part of the slope

下载: 全尺寸图片 幻灯片

图  10  结构面识别结果与实际边坡结构面的对比

Figure  10.  Comparison between structural plane recognition results and actual slope structural plane

下载: 全尺寸图片 幻灯片

3.   边坡结构面产状统计信息采集

在聚类分析过程中对6组结构面的每个子聚类进行平面拟合，分别得到其平面方程参数，利用这些参数可计算得到每个结构面的倾角和走向，进而得到其倾向^[19-20]。因为聚类结果中存在很多小的聚类，难以形成有效的结构面，所以设定一个阈值500，每个子聚类中的所有点数量不小于500时才在输出结果中显示。利用Python编写的计算程序得到边坡稳定性定量评价需要的结构面产状信息。

对边坡上层岩体进行分析，得到拟合平面参数及结构面产状信息（表1），J₁—J₆所在结构面组中出露的结构面数量依次为70，66，103，49，55，20。通过统计分析，得到各组结构面倾角和倾向的概率分布直方图及函数拟合曲线，如图11所示，μ为均值，σ为标准差。

表  1  上层岩体拟合平面参数及结构面产状信息

Table  1.  Parameters of fitting plane and orientation information of structural plane of rock mass at the upper part of the slope

主极点	聚类号	点数量	拟合平面方程参数				倾角/（°）	倾向/（°）
			Aij	Bij	Cij	Dij
J1	1	62 356	−0.407	0.813	−0.416	−49.929	65.43	153.40
	2	34 681	0.268	−0.877	0.399	42.268	66.46	163.01
J2	1	54 135	−0.010	0.860	−0.510	−33.920	59.33	179.36
	2	37 259	0.253	0.810	−0.529	−39.055	58.04	197.33
J3	1	62 073	0.334	−0.288	0.898	−58.570	26.17	130.82
	2	32 684	0.323	−0.492	0.808	−10.036	36.07	146.75
J4	1	88 328	0.269	0.947	−0.177	−44.659	79.83	195.89
	2	7 973	0.241	0.961	−0.136	−47.637	82.17	194.10
J5	1	37 880	−0.902	0.420	−0.100	−11.179	84.26	114.98
	2	27 762	0.748	−0.620	0.237	33.931	76.31	129.69
J6	1	20 051	0.643	0.766	−0.020	−31.941	88.88	220.03
	2	6 660	0.600	0.792	−0.116	0.267	83.32	217.13

下载: 导出CSV 

| 显示表格

图  11  结构面组J₁—J₆倾角和倾向概率分布直方图与拟合曲线的比较

Figure  11.  Comparison of probability distribution histograms and fitting curves for dip and dip angles of structural plane groups J₁-J₆

下载: 全尺寸图片 幻灯片

边坡下层岩体拟合平面参数及结构面产状信息见表2。J₇—J₁₁所在结构面组中出露的结构面数量依次为102，92，41，42，12。同理，通过统计分析可得到各组结构面倾角和倾向的概率分布直方图、函数拟合曲线和结构面产状分布特征。结构面产状统计特征见表3。边坡结构面产状概率统计分布特征可以为高边坡概率稳定性的定量评价提供重要的数据来源。

表  2  下层岩体拟合平面参数及结构面产状信息

Table  2.  Parameters of fitting plane and orientation information of structural plane of rock mass at the lower part of the slope

主极点	聚类号	点数	拟合平面方程参数				倾角/（°）	倾向/（°）
			Aij	Bij	Cij	Dij
J7	1	52 786	−0.297	0.687	−0.664	−13.625	48.42	156.60
	2	23 217	0.317	−0.645	0.696	22.360	45.93	153.80
J8	1	43 880	0.001	−0.825	0.565	30.001	55.62	179.93
	2	28 671	−0.197	0.846	−0.495	−26.153	60.35	166.89
J9	1	90 517	0.742	−0.622	0.251	15.808	75.44	129.97
	2	60 842	0.723	−0.691	0.001	56.813	89.96	133.68
J10	1	27 216	0.482	0.793	−0.372	15.417	68.19	211.28
	2	10 396	0.252	0.928	−0.274	−42.152	74.07	195.16
J11	1	7 602	0.158	−0.987	0.024	50.753	88.61	350.88
	2	1 690	0.137	−0.989	0.049	53.162	87.16	352.14

下载: 导出CSV 

| 显示表格

表  3  边坡上层和下层岩体结构面产状统计特征

Table  3.  Statistical characteristics of rock mass structural plane orientation at the upper and lower part of the slope

上层岩体						下层岩体
主极点	产状	分布特征	均值$\mu $	标准差$\sigma $		主极点	产状	分布特征	均值$\mu $	标准差$\sigma $
J1	倾角	对数正 态分布	68.32	2.92		J7	倾角	正态分布	43.39	4.42
	倾向	对数正 态分布	156.21	6.78			倾向	正态分布	154.31	8.50
J2	倾角	正态分布	56.97	5.31		J8	倾角	对数正 态分布	59.33	5.36
	倾向	正态分布	184.40	8.40			倾向	对数正 态分布	170.74	5.73
J3	倾角	正态分布	34.95	6.64		J9	倾角	正态分布	72.74	7.67
	倾向	正态分布	142.94	15.87			倾向	对数正 态分布	138.26	12.24
J4	倾角	正态分布	81.39	2.93		J10	倾角	正态分布	72.09	16.48
	倾向	对数正 态分布	201.24	3.16			倾向	—	—	—
J5	倾角	正态分布	77.58	4.47		J11	倾角	—	—	—
	倾向	正态分布	127.87	6.79			倾向	—	—	—
J6	倾角	正态分布	83.22	4.97
	倾向	对数正 态分布	216.16	3.17

下载: 导出CSV 

| 显示表格

4.   讨论

对比分析研究区高边坡上层和下层岩体结构面的概率分布直方图及拟合结果可知，大部分结构面的函数拟合效果较好，仅有少数结构面的拟合效果较差。其可能原因主要如下：① 无人机获取的图像数量不够多，导致部分结构面细节丢失，进而增大了点云数据提取过程中的误差。② 研究区边坡下层岩体表面局部存在植被，干扰点云数据聚类分析，影响结构面产状的计算结果。③ 一些结构面组出露的结构面数量很少，无法对其进行准确拟合。

相较于当前结构面识别研究^[21-22]，本文方法具有以下优势。

1） 相较于三维激光扫描方法，本文方法具有效率高、操作简单和成本低等优势，尤其在一些人力无法到达和不利于架设观测站的情况（如水电站库岸边坡等），更能突出无人机机动灵活的优势。

2） 相较于其他数据处理算法（K-Means算法等），DBSCAN算法不需要事先指定聚类个数，能够自动识别簇的数量和形状，同时能识别离群点，有效过滤噪声数据，使得结构面识别过程更高效、结果更准确，且无需过多人工干预。

5.   结论

1） 提出了基于无人机摄影测量技术和Context Capture系列软件的边坡结构面点云数据获取方法，具有效率高、操作简单和成本低等优势，为解决高边坡结构面勘测难题提供了有效途径。

2） 采用KNN算法和DBSCAN算法进行边坡结构面识别，实现了复杂边坡结构面的准确高效识别和产状信息快速提取，结构面识别结果与工程实际吻合。

3） 岩质高边坡结构面识别及产状统计信息采集方法可获取高边坡结构面数量、产状信息及其统计特征，大部分结构面倾角和倾向概率分布与实测数据拟合较好，为高边坡裂隙网络模型构建及稳定性分析提供了重要数据来源。

4） 下一步将继续开展岩质高边坡结构面网络模型构建和边坡稳定性与可靠度分析，并揭示结构面分布对边坡稳定性的影响机制，最终实现科学指导水利水电工程、公路和露天矿边坡的防护与安全运行，降低边坡支护成本和失稳破坏造成的潜在损失。