本发明公开一种基于点云数据的露天矿台阶线提取和边坡监测的方法，该方法通过渐进形态学滤波来消除挡土墙等人工地物对露天矿台阶线提取的干扰，提高了精度；根据台阶线特征点的几何属性在提取特征点前先计算每点曲率，将曲率大于阈值的点作为潜在特征点，再对潜在特征点进行细化，大大减少了计算复杂度，提高了整个运行速度；将曲率指数结合AGPN算法提高了提取台阶线上的特征点的精度，本发明的方法从整体上提高了露天矿台阶线提取的自动化程度。此外通过对两期台阶线提取结果分析可以得出露天矿边坡坡度变化和边坡位移量，这对露天矿的安全生产具有重要意义。
1.一种基于点云数据的露天矿台阶线提取方法，其特征在于，包括如下步骤：步骤1：对露天矿点云数据使用渐进形态学滤波算法进行预处理，通过逐渐增加滤波器的窗口大小和不断增大的高程差阈值移除点云数据中的地物数据，提取地面点，过程如下：步骤1.1：渐进形态学滤波算法在进行第一次迭代时，将最小高程表面与初始滤波器的窗口大小作为输入，对三维激光点云数据作滤波运算，得到已滤波表面；
步骤1.2：比较已滤波表面在滤波前后的高程差是否小于高程差阈值，若小于则将窗口大小赋值给标记数组，来增加窗口大小；
所述高程差阈值根据研究区域的地形坡度确定，方法如下：假设坡度恒定，则地形的最大高程差，窗口大小和地形坡度之间存在下列公式：其中，dh0表示初始高程差阈值；s表示斜率；c表示格网大小；dhmax表示最大高程差阈值；
wk表示运行到第k次窗口大小，dhT,K表示高程差阈值，dhmax(t),k表示地形的最大高程差；
步骤1.3：将上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入，继续进行滤波运算迭代，直到窗口大小为最大窗口；
步骤2：基于顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取台阶线特征点，过程如下：步骤2.1：点P为点云中任意一点，对点云中点P及其邻域内的点使用最小二乘法拟合曲
面，并对其进行参数化：r(x,y)＝(x,y,ax+bxy+cy)，其中，x为点P的横坐标，y为点P的纵坐标，a、b、c为对曲面参数化后的参数；
步骤2.2：设r(x,y)的偏导数分别为rx，ry，rxx，rxy，ryy，然后根据以下公式计算平均曲率H和高斯曲率K：2
其中， 为曲面的单位法向量；E＝rx·rx＝1+(2ax+by) ；
F＝rx·ry＝(2ax+by)(2cy+bx)，G＝ry·ry＝1+(2cy+bx) ，步骤2.3：最后通过下列公式计算主曲率：其中，kmax为最大主曲率，kmin为最小主曲率；
步骤2.5：将曲率指数大于曲率指数阈值的点定义为潜在特征点，对潜在特征点使用顾及邻域几何属性的三维边缘检测算法AGPN细化特征点，将潜在特征点中的内点集进行法向量优化；将潜在特征点中的外点集直接定义为非特征点，过程如下：步骤2.5.1：对曲率指数大于曲率指数阈值的点其中任意一点P使用k‑d tree算法检索其邻域点集S；
步骤2.5.2：对邻域点集S使用RANSAC算法拟合平面，将点集S分为内点集和外点集；
步骤2.5.3：若点P在内点集中，则计算拟合平面的法向量，将其作为点P的法向量进行法向量优化；并将潜在特征点中的外点集直接定义为非特征点；
步骤2.6：根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系，设空间向量 在x轴上，向量 在y轴上，优化后的法向量为 计算角度阶跃Gθ，公式如下：Gθ＝max(θi+1‑θi)i＝1,…N‑1 (9)式中，空间向量 是由当前点P和其邻域内拟合平面提取的内点构成的，θi为向量与x轴的夹角；
步骤2.7：将角度阶跃大于阈值的点作为提取的台阶线特征点；
步骤3：采用移动最小二乘法拟合台阶线，得到提取后的露天矿台阶线，过程如下：步骤3.1：对提取出的台阶线特征点进行曲线拟合，在要进行曲线拟合的部分中，其目标子区域在拟合过程中使用的函数记作f(x)：T
α(x)＝[a(x1),a(x2),……,a(xm)] (11)T
p(x)＝[p(x1),p(c2),……,p(xm)] (12)T
其中，α(x)代表的为拟合曲线的一组系数；p (x)代表的为m次多项式空间的一组基，a(xi)为待定的第i个系数，p(xi)为基函数；
步骤3.2：若计算点x的邻域内有n个节点，近似函数f(x)在这些节点处的误差的加权平方和如下：其中，yi是x＝xi处的节点值，w(xi)是具有紧支特性的权函数；
步骤3.3：解出待定参数α(x)，让公式(13)中的J取得最小值；
α(x)＝A (x)B(x)y (15)B(x)＝[w(x1)p(x1),w(x2)p(x2),……,w(xn)p(xn)] (17)计算时需要用到形函数
2.一种基于台阶线提取的边坡监测方法，其特征在于，包括如下步骤：步骤1：采用权利要求1所述的基于点云数据的露天矿台阶线提取方法进行台阶线提取；
步骤2：通过对两个不同时期的露天矿台阶线提取结果进行分析，计算边坡曲面上每一部分的坡度值和位移变化量；
使用移动最小二乘法拟合台阶线后，将边坡曲面上每一部分对应的台阶线特征点坐标带入公式中计算坡度与位移变化量；
坡度i，即两点的高程差与其水平距离的百分比，其计算公式如下：h＝zi+1‑zi (20)式中：h为高程差，l为水平距离，x，y，z为台阶线特征点坐标；
一种基于点云数据的露天矿台阶线提取和边坡监测的方法
技术领域
[0001] 本发明涉及边坡数据采集与监测技术领域，尤其涉及一种基于点云数据的露天矿台阶线提取和边坡监测的方法。
背景技术
[0002] 露天矿在进行开采工作时，为了便于从上至下逐层开采，要把矿区内的矿物、岩石划分为具有一定厚度的水平分层，这些分层的工作面就是台阶。露天矿台阶线对露天矿的勘测、采矿设计、计划编制、区块设计、爆破设计、矿区划分等业务均具有重要的作用。
[0003] 现阶段在实际生产过程中提取露天矿上台阶线等特征线使用的主要方法是根据DSM进行人工绘制，这种方法工作量大、效率低。利用轻便的无人机，获取露天矿序列影像，进行三维建模，检测采场台阶线，自动精确解算剥离量，对加强矿山生产监控和成本管理具有重要意义。史红霞等在中国机械工程中发表的《基于法向量区域聚类分割的点云特征线提取》中提出采用自适应邻域的主成分分析法(PCA)估算模型法向，利用萤火虫算法(FA)优化的模糊C均值(FCM)算法对法向进行聚类实现模型的有效分割。构造点集剔除与合并准则从各分割块边界点集中析取候选特征点，再以局部邻域主轴方向为基准提取特征点的方法。柏宏强等在激光与光电子学进展中发表的《基于三维激光点云特征线提取的溶洞多分辨率三维重建方法研究》中提出采用改进邻近点几何特征提取溶洞特征值，增加法向量角作为检测特征点的参数；用社会粒子群(SPSO)算法与模糊C‑均值(FCM)聚类算法实现点云分类；采用折线生长方法将特征点连接成特征线的方法。上述方法中，存在边缘和特征线检测不全，受噪声点干扰大，且当场景复杂时，在非平面区域会产生现实中不存在的边缘或特征线等问题。
发明内容
[0004] 针对上述现有技术的不足，本发明提供一种基于点云数据的露天矿台阶线提取和边坡监测的方法。
[0005] 为解决上述技术问题，本发明所采取的技术方案是：一种基于点云数据的露天矿台阶线提取方法，包括如下步骤：
[0006] 步骤1：对露天矿点云数据使用渐进形态学滤波算法进行预处理，通过逐渐增加滤波器的窗口大小和不断增大的高程差阈值移除点云数据中的地物数据，提取地面点，过程如下：
[0007] 步骤1.1：渐进形态学滤波算法在进行第一次迭代时，将最小高程表面与初始滤波器的窗口大小作为输入，对三维激光点云数据作滤波运算，得到已滤波表面；
[0008] 步骤1.2：比较已滤波表面在滤波前后的高程差是否小于高程差阈值，若小于则将窗口大小赋值给标记数组，来增加窗口大小；
[0009] 所述高程差阈值根据研究区域的地形坡度确定，方法如下：
[0010] 假设坡度恒定，则地形的最大高程差，窗口大小和地形坡度之间存在下列公式：
[0011]
[0012]
[0013] 其中，dh0表示初始高程差阈值；s表示斜率；c表示格网大小；dhmax表示最大高程差阈值；wk表示运行到第k次窗口大小，dhT，K表示高程阈值，dhmax(t)，K表示地形的最大高程差。
[0014] 步骤1.3：将上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入，继续进行滤波运算迭代，直到窗口大小为最大窗口。
[0015] 步骤2：基于顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取台阶线特征点，过程如下：
[0016] 步骤2.1：点P为点云中任意一点，对点云中点P及其邻域内的点使用最小二乘法拟2 2
合曲面，并对其进行参数化：r(xy)＝(x，y，ax+bxy+cy)，其中，x为点P的横坐标，y为点P的纵坐标，a、b、c为对曲面参数化后的参数；
[0017] 步骤2.2：设r(x，y)的偏导数分别为rx，ry，rxx，rxy，ryy，然后根据以下公式计算平均曲率H和高斯曲率K：
[0018]
[0019]2
[0020] 其中， 为曲面的单位法向量；E＝rx·rx＝1+(2ax+by) ；
[0021] F＝rx·ry＝(2αx+by)(2cy+bx)，G＝ry·ry＝1+(2cy+bx)2，
[0022]
[0023] 步骤2.3：最后通过下列公式计算主曲率：
[0024]
[0025]
[0026] 其中，kmax为最大主曲率，kmin为最小主曲率；
[0027] 步骤2.4：计算点云中点P的曲率指数：
[0028]
[0029] 步骤2.5：将曲率指数大于曲率指数阈值的点定义为潜在特征点，对潜在特征点使用顾及邻域几何属性的三维边缘检测算法AGPN细化特征点，将潜在特征点中的内点集进行法向量优化；将潜在特征点中的外点集直接定义为非特征点，过程如下：
[0030] 步骤2.5.1：对曲率指数大于曲率指数阈值的点其中任意一点P使用k‑d tree算法检索其邻域点集S；
[0031] 步骤2.5.2：对邻域点集S使用RANSAC算法拟合平面，将点集S分为内点集和外点集；
[0032] 步骤2.5.3：若点P在内点集中，则计算拟合平面的法向量，将其作为点P的法向量进行法向量优化，将潜在特征点中的外点集直接定义为非特征点。
[0033] 步骤2.6：根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系，设空间向量 在x轴上，向量 在y轴上，优化后的法向量为 计算角度阶跃Gθ，公式如下：
[0034]
[0035] Gθ＝max(θi+1‑θi)(i＝1…N‑1) (9)
[0036] 式中，空间向量 是由当前点P和其邻域内拟合平面提取的内点构成的，θi为向量 与向量 夹角；
[0037] 步骤2.7：将角度阶跃大于阈值的点作为提取的台阶线特征点。
[0038] 步骤3：采用移动最小二乘法拟合台阶线，得到提取后的露天矿台阶线，过程如下：
[0039] 步骤3.1：对提取出的台阶线特征点进行曲线拟合，在要进行曲线拟合的部分中，其目标子区域在拟合过程中使用的函数记作f(x)：
[0040]
[0041] α(x)＝[a1(x)，ax(x)，……，am(x)]T (11)
[0042] pT(x)＝[p1(x)，p2(x)，......，pm(x)] (12)
[0043] 其中，α(x)代表的为待定的系数矢量；pT(x)代表的为m次多项式空间的一组基，ai(x)为待定的系数矢量，pi(x)为基函数；
[0044] 步骤3.2：若计算点x的邻域内有肝个节点，近似函数f(x)在这些节点处的误差的加权平方和如下：
[0045]
[0046] 其中，yi是x＝xi处的节点值，wi(x)是具有紧支特性的权函数；
[0047] 步骤3.3：解出待定参数α(x)，让公式(15)中的J取得最小值；
[0048] 对J取极值即对α(x)求导得
[0049]
[0050] α(x)＝A‑1(x)B(x)y (15)
[0051]
[0052] B(x)＝[w(x1)p(x1)，w(x2)p(x2)，......，w(xn)p(xn)] (17)
[0053] 计算时需要用到形函数
[0054]
[0055] 步骤3.4：获得最终的拟合函数：
[0056]
[0057] 另一方面，本发明提供一种基于台阶线提取的边坡监测方法，包括如下步骤：
[0058] 步骤1：上述基于点云数据的露天矿台阶线提取方法进行台阶线提取；
[0059] 步骤2：通过对两期露天矿台阶线提取结果进行分析，计算边坡曲面上每一部分的坡度值和位移变化量。
[0060] 使用移动最小二乘法拟合台阶线后，将边坡曲面上每一部分对应的台阶线特征点坐标带入公式中计算坡度与位移变化量；
[0061] 坡度，即两点的高程差与其水平距离的百分比，其计算公式如下：
[0062] h＝zi+1‑zi (20)
[0063]
[0064]
[0065] 式中：h为高程差，l为水平距离，x，y，z为台阶线特征点坐标；
[0066] 位移变化量计算公式如下：
[0067]
[0068] 式中，x，y，z为台阶线特征点坐标。
[0069] 采用上述技术方案所产生的有益效果在于：
[0070] 1、本发明提供的方法通过渐进形态学滤波来消除挡土墙等人工地物对露天矿台阶线提取的干扰，提高了精度。在未经过滤波处理时提取露天矿台阶线，挡土墙的脊点会被提取出来，该点并不是台阶线特征点。所以使用渐进形态学滤波消除干扰。
[0071] 2、本发明根据台阶线特征点的几何属性在提取特征点前先计算每点曲率，将曲率大于阈值的点作为潜在特征点，再对潜在特征点进行细化。这样大大减少了计算复杂度，提高了整个运行速度。
[0072] 3、本发明通过对两期台阶线提取结果分析可以得出露天矿边坡坡度变化和边坡位移量，这对露天矿的安全生产具有重要意义。
附图说明
[0073] 图1为本发明实施例中基于点云数据的露天矿台阶线提取和边坡监测的方法流程图；
[0074] 图2(a)为本发明实施例中点P位于折叠边缘时的空间向量在拟合平面上的投影情况示意图；
[0075] 图2(b)为本发明实施例中点P位于局部拟合平面内部时的空间向量在拟合平面上的投影情况示意图。
[0076] 图3为本发明实施例中台阶线提取效果图。
具体实施方式
[0077] 下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。
[0078] 本实施例中基于点云数据的露天矿台阶线提取和边坡监测的方法如图1所示，其中基于点云数据的露天矿台阶线提取方法，包括如下步骤：
[0079] 步骤1：对露天矿点云数据使用渐进形态学滤波算法进行预处理，通过逐渐增加滤波器的窗口大小和不断增大的高程差阈值移除点云数据中的地物数据，提取地面点，过程如下：
[0080] 步骤1.1：渐进形态学滤波算法在进行第一次迭代时，将最小高程表面与初始滤波器的窗口大小作为输入，对三维激光点云数据作滤波运算，得到已滤波表面；
[0081] 步骤1.2：比较已滤波表面在滤波前后的高程差是否小于高程差阈值，若小于则将窗口大小赋值给标记数组，来增加窗口大小；
[0082] 所述高程差阈值根据研究区域的地形坡度确定，方法如下：
[0083] 假设坡度恒定，则地形的最大高程差，窗口大小和地形坡度之间存在下列公式：
[0084]
[0085]
[0086] 其中，dh0表示初始高程差阈值；s表示斜率；c表示格网大小；dhmax表示最大高程差阈值；wk表示运行到第k次窗口大小，dhT，K表示高程阈值，dhmax(t)，K表示地形的最大高程差。
[0087] 步骤1.3：将上一迭代中获得的已滤波表面和增加的窗口大小作为滤波器的输入，继续进行滤波运算迭代，直到窗口大小为最大窗口。
[0088] 所述地物数据包括：车辆、植被、建筑、挡土墙等。消除地物对露天矿台阶线提取的干扰，提高了精度。
[0089] 步骤2：基于顾及邻域几何属性的三维边缘检测与曲率指数加权方法提取台阶线特征点，过程如下：
[0090] 台阶线上的特征点均是具有曲率极值的点，因此本发明根据曲率初步筛选潜在特征点。
[0091] 曲率是描述几何体弯曲程度的量，例如二维曲线偏离直线的程度，或者三维曲面偏离平面的程度，也可以用来确定曲面类型。
[0092] 平均曲率、主曲率和高斯曲率是曲率的三个基本要素，(1)平均曲率：是指曲面上过一点P有无数个不同方向的曲线，不同方向具有不同的曲率，任意取其中一对相互垂直的正交曲率，取平均值。其计算公式为：甘＝(K1+K2)/2，式中K1，K2为点P的一对正交曲率。(2)主曲率：过曲面上一点P上具有无穷个正交曲率，其中存在一条曲线使得该曲线的曲率为极大，这个曲率为极大值Kmax，垂直于极大曲率面的曲率为极小值Kmin，极值方向为主方向。(3)高斯曲率：其计算公式为K＝Kmax·Kmin，表示的是点P的弯曲程度。
[0093] 步骤2.1：点P为点云中任意一点，对点云中点P及其邻域内的点使用最小二乘法拟2 2
合曲面，并对其进行参数化：r(x，y)＝(x，y，ax+bxy+cy)，其中，x为点P的横坐标，y为点P的纵坐标，a、b、c为对曲面参数化后的参数；
[0094] 步骤2.2：设r(x，y)的偏导数分别为rx，ry，rxx，rxy，ryy，然后根据以下公式计算平均曲率H和高斯曲率K：
[0095]
[0096]2
[0097] 其中， 为曲面的单位法向量；E＝rx·rx＝1+(2ax+by) ；
[0098] F＝rx·ry＝(2ax+by)(2cy+bx)，G＝ry·ry＝1+(2cy+bx)2，
[0099]
[0100] 步骤2.3：最后通过下列公式计算主曲率：
[0101]
[0102]
[0103] 其中，kmax为最大主曲率，kmin为最小主曲率；
[0104] 步骤2.4：计算点云中点P的曲率指数：
[0105]
[0106] 曲率指数与坐标系无关，具有旋转和位移不变性，在三维特征点检测领域中被广泛使用。
[0107] 步骤2.5：将曲率指数大于曲率指数阈值的点定义为潜在特征点，对潜在特征点使用顾及邻域几何属性的三维边缘检测算法AGPN细化特征点，将潜在特征点中的内点集进行法向量优化；将潜在特征点中的外点集直接定义为非特征点；
[0108] AGPN算法引入了平面模型的RANSAC算法，RANSAC算法可以从一组包含“局外点”的观测数据集中通过迭代方式估计数学模型的参数，并且该数学模型能够适用于所有的“局内点”，它能从包含大量局外点的数据集中估计出高精度的参数，具有很好的鲁棒性。
[0109] 步骤2.5的过程如下：
[0110] 步骤2.5.1：对曲率指数大于曲率指数阈值的点其中任意一点P使用k‑d tree算法检索其邻域点集S；
[0111] 步骤2.5.2：对邻域点集S使用RANSAC算法拟合平面，将点集S分为内点集和外点集；
[0112] 露天矿地形复杂，少有规则的曲面，当检测不规则曲面相交构成的边缘时，由于曲面的微小局部可以看成是平面，平面模型可以有效地近似局部曲面的几何结构。如果潜在特征点中一点P是折叠边缘，则使用k‑d树搜索邻域后，空间向量 在点P邻域拟合平面上的投影会有较大的夹角，如图2(a)所示，该折叠边缘同时位于拟合平面和邻域内另一个相交表面上。当点P不是折叠边缘但也属于内点时，那么该点位于内点的拟合平面的内部，此时，空间向量 在拟合平面上的投影之间不存在较大的夹角，是连续分布的，如图2(b)所示。此外，若内点集的点数过少，则认为点P为噪声点。
[0113] 步骤2.5.3：若点P在内点集中，则计算拟合平面的法向量，将其作为点P的法向量进行法向量优化；并将潜在特征点中的外点集直接定义为非特征点。
[0114] 法向量优化过程为：对点云中某点P的法向量估计方法是对点P使用k‑d树搜索邻域，对其邻域采用最小二乘法拟合平面，构造协方差矩阵，求解协方差矩阵的特征值和特征向量，然后选择最小特征值对应的特征向量，并进行单位化，则该向量为点P法向量。当点P的邻域有多个平面时，优化后的法向量将会垂直于其中一个平面。
[0115] 步骤2.6：根据优化后的法向量和拟合平面上一对相互垂直的向量建立空间直角坐标系，设空间向量 在x轴上，向量 在y轴上，优化后的法向量为 计算角度阶跃Gθ，公式如下：
[0116]
[0117] Gθ＝max(θi+1‑θi)(i＝1…N‑1) (9)
[0118] 式中，空间向量 是由当前点P和其邻域内拟合平面提取的内点构成的，θi为向量 与向量 夹角；
[0119] 步骤2.7：将角度阶跃大于阈值的点作为提取的台阶线特征点。
[0120] 步骤3：采用移动最小二乘法拟合台阶线，得到提取后的露天矿台阶线；
[0121] 传统的最小二乘法拟合曲线是在计算偏差后，让平方和达到最小值，从而将问题转化为求解线性方程组，最终得到拟合曲线。该方法具有全局性，但是不能满足定位的需要，不适用于具有大量数据集的模型。与传统的最小二乘法相比，移动最小二乘法在最终拟合曲线的定位和光滑性方面有较好的效果，且拟合精度也比较高。移动最小二乘法有两个明显的优点：第一点是不同于传统方法中的多项式，拟合过程中使用的函数由基函数和系数向量组成；第二点是引入了紧支的概念，即点x处的值y只受x支撑域内节点的影响。
[0122] 步骤3的过程如下：
[0123] 步骤3.1：对提取出的台阶线特征点进行曲线拟合，在要进行曲线拟合的部分中，其目标子区域在拟合过程中使用的函数记作f(x)：
[0124]
[0125] α(x)＝[a1(x)，ax(x)，......，am(x)]T (11)
[0126] pT(x)＝[p1(x)，p2(x)，......，pm(x)] (12)
[0127] 其中，α(x)代表的为待定的系数矢量；pT(x)代表的为m次多项式空间的一组基，ai(x)为待定的系数矢量，pi(x)为基函数；
[0128] 步骤3.2：若计算点x的邻域内有n个节点，近似函数f(x)在这些节点处的误差的加权平方和如下：
[0129]
[0130] 其中，yi是x＝xi处的节点值，wi(x)是具有紧支特性的权函数；
[0131] 步骤3.3：解出待定参数α(x)，让公式(15)中的J取得最小值；
[0132] 对J取极值即对α(x)求导得：
[0133]
[0134] α(x)＝A‑1(x)B(x)y (15)
[0135]
[0136] B(x)＝[w(x1)p(x1)，w(x2)p(x2)，......，w(xn)p(xn)] (17)
[0137] 计算时需要用到形函数
[0138]
[0139] 步骤3.4：获得最终的拟合函数：
[0140]
[0141] 在本实施例中，某露天矿的台阶线提取效果如图3所示。
[0142] 另一方面，本发明提供一种基于台阶线提取的边坡监测方法，包括如下步骤：
[0143] 步骤1：上述基于点云数据的露天矿台阶线提取方法进行台阶线提取；
[0144] 步骤2：通过对两期露天矿台阶线提取结果进行分析，计算边坡曲面上每一部分的坡度值和位移变化量。
[0145] 使用移动最小二乘法拟合台阶线后，将边坡曲面上每一部分对应的台阶线特征点坐标带入公式中计算坡度与位移变化量；
[0146] 坡度是地形描述中比较常用的参数，是一个具有大小与方向的矢量，坡度是指曲面上某一点的垂直方向Z与法线方向N之间的夹角值。
[0147] 坡度，即两点的高程差与其水平距离的百分比，其计算公式如下：
[0148] h＝zi+1‑zi (20)
[0149]
[0150]
[0151] 式中：h为高程差，l为水平距离，x，y，z为台阶线特征点坐标；
[0152] 位移变化量计算公式如下：
[0153]
[0154] 式中，x，y，z为台阶线特征点坐标。
