本发明属于工程技术领域，涉及一种露天矿山边坡岩体节理产状快速精细取值方法，特别的是本发明是采用三维激光扫描测量方法获得了的结构面产状数据，利用k均值聚类方法进行了结构面分组，在各组结构面产状统计分析的基础上，结合野外现场测量结果，确定边坡各组岩体结构面产状的统计分布规律与统计特征值，野外现场调查与室内统计分析相结合，显著提高了露天矿山边坡岩体结构面产状确定的可靠性。
背景技术：
结构面是边坡岩体中由力学强度较低的部位或岩性相对软弱的夹层构成的不连续面，岩体的变形与稳定性主要取决于结构面的发育条件。露天矿山边坡的稳定性受岩体结构面控制，结构面的产状、规模、组合型式及其与边坡面的空间组合关系决定了矿山边坡潜在的破坏模式。同时，结构面的力学特性决定了矿山边坡潜在滑移破坏的可能性。为了实现露天矿山边坡稳定性精细评价，必须建立正确的计算模型，调查清楚边坡岩体结构面的发育特征(如产状、迹长、间距等)是边坡稳定性评价的基础，其中结构面产状精细确定是判断边坡潜在破坏模式、构建计算模型的关键。
目前，岩体结构面产状测量方法主要包括：测线法、窗口统计法；、钻孔定向取芯法；摄影测量法；三维激光扫描法。近年来，一些代表性的岩体结构面产状研究如下：
吴乐文(2014)公开了一种基于激光或红外测距仪的空间测距原理的种非接触式结构面产状测量方法，突破了传统地质罗盘仪接触式测量的应用限制。
朱爱玺(2014)公开了一种基于图像的地下洞室掌子面快速安全预测方法，以红外相机拍照为测量手段，对掌子面图像进行图像预处理与聚类分析，进行边界提取，得到掌子面层理节理信息，最终根据提取到的结构面信息建立三维地质模型。
许振浩等(2015)公开了一种车载式隧道全空间裂隙网络检测成像与预警系统及方法，该系统集激光扫描与数码扫描于一体，能够获取隧道全空间裂隙网络图像、掌子面多元裂隙网络信息。
侯清波(2016)公开了一种地质结构面产状测量方法，在野外现场采用非接触式测量方法和接触式手持android采集设备测量结构面产状，实现了对远距离或难以接触的结构面进行产状测量。
张文等(2017)公开了一种基于夜视仪和微型三维激光扫描仪的岩体结构面原位信息采集自动化装置，通过主机的图像处理模块能自动解译结构面的位置、产状、迹长、密度与宽度等参数。
姜清辉等(2017)公开了一种开挖边坡岩块信息提取与三维重构方法，运用坡面三维激光扫描点云数据，实现了岩石块体信息的自动化提取，通过岩块统计信息来评价岩体完整程度。
胡建平(2018)公开了一种基于无人机的结构面产状及性状的多角度识别方法，利用无人机对确定范围的结构面进行拍照，采用smart3dcapture软件形成三维实体影像图，对三维实体影像图进行地质解译，确定结构面空间出露位置、长度、产状及性状。
这些研究成果为岩体结构面产状的智能化提取提供了宝贵的研究基础与经验。然而，在前人研究中并没有给出详细、可靠的岩体结构面产状快速精细取值的方法。岩体结构面产状测量仅仅做到每一条结构面的精确测量是不充分的，必须以此为基础进行结构面的分组，依据分组后得到各组结构面的统计测量结果，确定出具有代表性的结构面产状。另一方面，目前有学者已经在结构面产状的聚类分析方法作了大量工作，其中在k均值聚类算法是目前广泛应用的方法。然而，在k均值聚类算法中，聚类数目k需要事先指定，k值的确定非常难以估计，事先无法预知给定的数据集应该分成多少类别最为合适。此外，该方法对初值敏感，对不同的初值可能会导致不同的聚类结果，容易陷入局部最优。因此，针对目前技术方法的局限性，亟需提出一套完整的露天矿山边坡岩体结构面产状快速精细取值方法。
技术实现要素：
为了克服现有露天矿山边坡岩体结构面产状精细快速取值问题，本发明提供一种基于三维激光扫描的露天矿山边坡岩体节理产状快速精细取值方法。
本发明解决其技术问题所采用的技术方案是：
一种基于三维激光扫描的露天矿山边坡岩体节理产状快速精细取值方法，所述方法包括以下步骤：
(1)现场结构面产状接触式测量；
(2)结构面三维激光扫描精细测量，过程如下：
2.1：根据观测边坡岩体范围与空间位置，选取岩体表面的某个测量点作为初始测量点；
2.2：根据各个测量点与三维激光扫描仪的可旋转测量头之间的距离，设定三维激光扫描仪的主机的测点空间位置和姿态；
2.3：等分辨率扫描边坡岩体结构面空间点云数据；
2.4：检查点云数据记录结果，对扫描过程中产生的噪声、孔洞进行修补，对局部稀疏数据增补处理；
(3)结构面点云数据分析与处理，过程如下：
3.1：计算拓扑构造后的点云中当前点与相邻点的距离与距离均值，通过距离阈值对点云数据中噪声点进行识别和剔除；
3.2：根据设定的测量装置自身空间坐标初始位置与现场结构面产状接触式测量步骤(2)中所测得的参考结构面产状，确定点云数据的空间三维坐标；
3.3：采用快速k近邻搜索方法，估计点云表面法线；
3.4：采用区域生长算法进行结构面识别，并对识别出的结构面进行分割提取；
3.5：对提取出的结构面进行平面拟合，获取所以结构面的产状数据，即倾向a1，倾角b1；
(4)结构面产状k均值聚类与统计分析。
进一步，所述步骤(4)的处理过程如下：
4.1：将现场结构面产状接触式测量确定的k组结构面倾向a0*，倾角b0*，作为每个聚类确定一个初始聚类中心，得到k个初始聚类中心；
4.2：根据相似性度量准则计算每个样本数据距聚类中心的距离；将每个样本数据分配到距它最近的聚类中心，得到k组数据；
4.3：对于每组结构面，采用特征模量分析方法求解每组数据的聚类中心，假设某组内存在l个数据，求解它们的聚类中心采用如下：
首先，采用如下公式计算矩阵s
式中(xi，yi，zi)(i＝1,2，…，l)为任意结构面的单位法向量；
然后，求解矩阵s的特征值(τ1，τ2，τ3)及其对应的特征向量(ξ1，ξ2，ξ3)，τ1＜τ2＜τ3，最大特征值对应的特征向量ξ3为组内l个向量的平均向量，该平均向量作为新的聚类中心；
4.4：按照4.2～4.3重复计算，直到所有聚类中心的位置都固定，结构面数据的组别分配也随之固定；
4.5：将4.4得到的以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据；
4.6：将4.5得到的每组结构面的倾角数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ,m+σ]；
4.7：判断k组结构面的初始聚类中心的倾角b0*，是否落在稳健区间[m-σ,m+σ]范围内，若在该范围内，结构面聚类分析完成；若不在该范围内，说明初始聚类数目k估计不正确，需要对k值进行修改，将新的聚类数目k1，k1＝k+1，带入4.1，重新按照4.1～4.6进行聚类，直到初始聚类中心的倾角b0*，落在稳健区间[m-σ,m+σ]范围内；
4.8：根据4.7得到的最后结构面分组情况，绘制各组结构面倾向、倾角的统计直方图，计算各组结构面的倾向、倾角平均值。
再进一步，所述步骤(1)的处理过程为：
1.1：通过露天矿山边坡分级分析确定所需研究的边坡研究范围，观察研究范围内岩体结构面发育的总体特征；
1.2：依据观察到的结构面发育总体特征，初步判断边坡岩体由k0组结构面构成，采用地质罗盘对各组结构面中露头完整、测量条件较好的结构面进行测量，作为参考结构面产状，照片记录相应结构面位置；
1.3：采用地质罗盘对每组结构面产状(倾向a0，倾角b0)进行测量，至少获得9个测点数据；
1.4：对于倾角小于等于80°的结构面，每组结构面产状数据依据结构面倾角大小，剔除最大2个值、最小2个值，计算其余5个产状数据的算术平均值倾向a0*和倾角b0*，倾向a0*，倾角b0*视为代表该组结构面总体规律的产状；
对于倾角大于80°的结构面，每组结构面产状数据依据结构面倾向，按产状数据集中原则归类剔除异常数据后，计算其余5个产状数据的算术平均值，倾向a0*，倾角b0*视为代表该组结构面总体规律的产状。
本发明采用三维激光扫描测量方法获得结构面产状数据，利用k均值聚类方法进行了结构面分组，并在各组结构面产状统计分析的基础上，结合野外现场测量结果，确定了边坡各组岩体结构面产状的统计分布规律与统计特征值，做到了野外现场调查与室内统计分析相结合。
本发明的有益效果主要表现在：采用三维激光扫描测量方法快速获得结构面产状数据，将现场结构面产状接触式测量结果作为初始聚类中心，利用k均值聚类方法进行了结构面分组，并在各组结构面产状统计分析的基础上，结合野外现场测量结果对分组结果进行检验，确定了边坡各组岩体结构面产状的统计分布规律与统计特征值。
附图说明
图1是露天矿山边坡结构面产状极点图。
图2是第一组结构产状统计分布规律。
图3是第二组结构产状统计分布规律。
具体实施方式
下面结合附图对本发明作进一步描述。
参照图1～图3，一种基于三维激光扫描的露天矿山边坡岩体节理产状快速精细取值方法，包括现场代表性结构面产状测量、结构面三维激光扫描精细测量、结构面点云数据分析与处理、结构面产状k均值聚类与统计分析四个方面内容，各个部分内容具体介绍如下：
(1)现场结构面产状接触式测量：
1.1：通过露天矿山边坡分级分析确定所需研究的边坡研究范围，观察研究范围内岩体结构面发育的总体特征；
1.2：依据观察到的结构面发育总体特征，初步判断边坡岩体由k0组结构面构成，采用地质罗盘对各组结构面中露头完整、测量条件较好的结构面进行测量，作为参考结构面产状，照片记录相应结构面位置；
1.3：采用地质罗盘对每组结构面产状(倾向a0，倾角b0)进行测量，至少获得9个测点数据；
1.4：对于倾角小于等于80°的结构面，每组结构面产状数据依据结构面倾角大小，剔除最大2个值、最小2个值，计算其余5个产状数据的算术平均值(倾向a0*，倾角b0*)，倾向a0*，倾角b0*可视为代表该组结构面总体规律的产状。
对于倾角大于80°的结构面，每组结构面产状数据依据结构面倾向，按产状数据集中原则归类剔除异常数据后，计算其余5个产状数据的算术平均值(倾向a0*，倾角b0*)，倾向a0*，倾角b0*可视为代表该组结构面总体规律的产状。
(2)结构面三维激光扫描精细测量
2.1：根据观测边坡岩体范围与空间位置，选取岩体表面的某个测量点作为初始测量点；
2.2：根据各个测量点与三维激光扫描仪的可旋转测量头之间的距离，设定三维激光扫描仪的主机的测点空间位置和姿态；
2.3：等分辨率扫描边坡岩体结构面空间点云数据；
2.4：检查点云数据记录结果，对扫描过程中产生的噪声、孔洞进行修补，对局部稀疏数据增补处理。
(3)结构面点云数据分析与处理
3.1：计算拓扑构造后的点云中当前点与相邻点的距离与距离均值，通过距离阈值对点云数据中噪声点进行识别和剔除；
3.2：根据设定的测量装置自身空间坐标初始位置与现场结构面产状接触式测量步骤(2)中所测得的参考结构面产状，确定点云数据的空间三维坐标；
3.3：采用快速k近邻搜索方法，估计点云表面法线；
3.4：采用区域生长算法进行结构面识别，并对识别出的结构面进行分割提取；
3.5：对提取出的结构面进行平面拟合，获取所以结构面的产状数据(倾向a1，倾角b1)。
(4)结构面产状k均值聚类与统计分析
4.1：将现场结构面产状接触式测量确定的k组结构面倾向a0*，倾角b0*，作为每个聚类确定一个初始聚类中心，得到k个初始聚类中心。
4.2：根据相似性度量准则计算每个样本数据距聚类中心的距离；将每个样本数据分配到距它最近的聚类中心，得到k组数据。
4.3：对于每组结构面，采用特征模量分析方法求解每组数据的聚类中心。假设某组内存在l个数据，求解它们的聚类中心可采用如下：
首先，采用如下公式计算矩阵s
式中(xi，yi，zi)(i＝1,2，…，l)为任意结构面的单位法向量。
然后，求解矩阵s的特征值(τ1，τ2，τ3)(τ1＜τ2＜τ3)及其对应的特征向量(ξ1，ξ2，ξ3)，最大特征值对应的特征向量ξ3为组内l个向量的平均向量，该平均向量作为新的聚类中心。
4.4：按照4.2～4.3重复计算，直到所有聚类中心的位置都固定，结构面数据的组别分配也随之固定。
4.5：将4.4得到的以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据。
4.6：将4.5得到的每组结构面的倾角数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算倾角数据的稳健区间[m-σ,m+σ]。
4.7：判断k组结构面的初始聚类中心的倾角b0*，是否落在稳健区间[m-σ,m+σ]范围内，若在该范围内，结构面聚类分析完成；若不在该范围内，说明初始聚类数目k估计不正确，需要对k值进行修改，将新的聚类数目k1(k1＝k+1)，带入4.1，重新按照4.1～4.6进行聚类，直到初始聚类中心的倾角b0*，落在稳健区间[m-σ,m+σ]范围内。
4.8：根据4.7得到的最后结构面分组情况，绘制各组结构面倾向、倾角的统计直方图，计算各组结构面的倾向、倾角平均值。
以浙江省绍兴市上虞区某矿山边坡岩体作为研究对象，介绍本发明的具体实施方式：
(1)现场结构面产状接触式测量
1.1：通过对某矿山边坡进行分级分析确定了矿山a边坡和矿山b边坡，以矿山a边坡作为研究对象，测量该边坡的高度与宽度，观察研究范围内岩体结构面发育的总体特征；
1.2：依据观察到的结构面发育总体特征，初步判断边坡岩体由k＝2组结构面构成，对各组结构面中露头完整、测量条件较好的结构面进行测量，作为参考结构面产状，照片记录相应结构面位置；
1.3：采用地质罗盘对每组结构面产状进行测量，每组获得9个测点数据；
1.4：结构面倾角均小于等于80°的结构面，每组结构面产状数据依据结构面倾角大小，剔除最大2个值、最小2个值，计算其余5个产状数据的算术平均值，代表第一组和第二组结构面总体规律的产状分别为59°∠121°、75°∠220°。
(2)结构面三维激光扫描精细测量
2.1：根据观测边坡岩体范围与空间位置，选取岩体表面的某个测量点作为初始测量点；
2.2：根据各个测量点与三维激光扫描仪的可旋转测量头之间的距离，设定三维激光扫描仪的主机的测点空间位置和姿态；
2.3：等分辨率扫描边坡岩体结构面空间点云数据；
2.4：检查点云数据记录结果，对扫描过程中产生的噪声、孔洞进行修补，对局部稀疏数据增补处理。
(3)结构面点云数据分析与处理
3.1：计算拓扑构造后的点云中当前点与相邻点的距离与距离均值，通过距离阈值对点云数据中噪声点进行识别和剔除；
3.2：根据设定的测量装置自身空间坐标初始位置与现场结构面产状接触式测量步骤(2)中所测得的参考结构面产状，确定点云数据的空间三维坐标；
3.3：采用快速k近邻搜索方法，估计点云表面法线；
3.4：采用区域生长算法进行结构面识别，并对识别出的结构面进行分割提取；
3.5：对提取出的结构面进行平面拟合，获取所以结构面的产状数据(倾向a1，倾角b1)，结构面产状的赤平投影极点图如图1所示。
(4)结构面产状k均值聚类与统计分析
4.1：将现场结构面产状接触式测量确定的k组结构面倾向a0*，倾角b0*，作为每个聚类确定一个初始聚类中心，得到k个初始聚类中心。
4.2：根据相似性度量准则计算每个样本数据距聚类中心的距离；将每个样本数据分配到距它最近的聚类中心，得到k组数据。
4.3：对于每组结构面，采用特征模量分析方法求解每组数据的聚类中心。假设某组内存在l个数据，求解它们的聚类中心可采用如下：
首先，采用如下公式计算矩阵s
式中(xi，yi，zi)(i＝1,2，…，l)为任意结构面的单位法向量。
然后，求解矩阵s的特征值(τ1，τ2，τ3)(τ1＜τ2＜τ3)及其对应的特征向量(ξ1，ξ2，ξ3)，最大特征值对应的特征向量ξ3为组内l个向量的平均向量，该平均向量作为新的聚类中心。
4.4：按照4.2～4.3重复计算，直到所有聚类中心的位置都固定，结构面数据的组别分配也随之固定。
4.5：将4.4得到的以单位法向量表示的结构面产状数据转换为以倾向、倾角表示的结构面产状数据。
4.6：将4.5得到的每组结构面的倾角数据进行统计分析，计算结构面倾角的平均值m与标准差σ，计算得到第一组结构面的倾角数据的稳健区间[54,82]，第二组结构面的倾角数据的稳健区间[64,81]。
4.7：判断可知两组结构面的初始聚类中心的倾角均落在稳健区间范围内。
4.8：根据4.7得到的最后结构面分组情况，分别绘制第一组、第二组结构面倾向、倾角的统计直方图如图2、3所示。第一组结构面的倾向的平均值为150.3°、倾角平均值为68.2°，第一组结构面的倾向的平均值为269.3°、倾角平均值为72.9°。