铁路建设和运营中边坡形变监测尤为重要。以往的形变监测多是单点测量，受限于监测点个数，并不能真实反映边坡的整体形变情况，并且碍于危险一些区域也无法进行观测。随着软硬件的不断升级，一些新的对地观测手段如遥感技术、InSAR、三维激光扫描仪等不断应用于形变监测^[1-7]。研究表明，三维激光扫描仪能够对大规模的形变进行监测, 比传统测量方法更优越^[8-12]。Hesse等^[13]利用点云数据构建数字高程模型进行形变研究，克服了点云数据单点精度低的问题; Gosliga等^[14]对隧道点云数据进行拟合构建隧道横断面，将隧道形变问题转换为空间平面参数变化问题; 常明等^[15]利用反求拼接参数误差模型及格网化方式对点云数据配准过程中的误差进行消除，提高形变监测精度; Kang等^[16]将三维点云数据转化成二维的反射值影像，利用尺度不变特征变换(SIFT)进行特征点匹配，比较2期点云数据的变化来检测建筑模型的变化。

本文针对铁路边坡实验区域，利用LMZ620三维激光扫描仪进行多期数据观测，获取的点云数据在100 m范围的单点重复观测精度为1 cm。在对场地进行数据获取时，点云数据单点存在一定的偶然误差，为削弱偶然误差对形变监测的影响，需要进行体素化处理。针对第一次观测的参考数据，利用种子生长法对场景内满足一定空间参数的点云数据进行分类，将点云数据分割成很多细小平面P_i, 作为当前场地的基准数据; 以后每隔一段时间获取场地观测数据，通过比较体素化后的点云数据到平面P_i的距离判断当前边坡的情况。图 1是本文方法的流程图。

1 点云数据处理 1.1 点云数据体素化

三维激光扫描仪获取的点云数据存在大量的冗余，而且单点精度不高，存在一定的偶然误差。直接利用原始点云数据进行形变分析，微小的变形可能会被数据的偶然误差所覆盖，影响最后的分析结果。为提高处理效率并削弱偶然误差的影响，利用偶然误差数学期望为0的特点，在进行数据处理时对点云数据进行体素化^[17]处理。

首先需要确定点云数据的范围(x_min, x_max, y_min, y_max, z_min, z_max), 依据观测场地的实际对体素大小进行选择。体素过大时，对形变区域边界的判定及场景内边缘区域的坐标值有一定的影响; 体素过小时，无法有效削弱偶然误差，且影响点云数据的处理速度。假定体素大小选择为l, 则场景点云坐标系3个方向的体素个数为：

$ \left\{ \begin{array}{l} {\rm{Nu}}{{\rm{m}}_x} = {\rm{ floor }}\left( {\frac{{{x_{{\rm{max}}}} - {x_{{\rm{min}}}}}}{l}} \right)\\ {\rm{Nu}}{{\rm{m}}_y} = {\rm{ floor }}\left( {\frac{{{y_{{\rm{max}}}} - {y_{{\rm{min}}}}}}{l}} \right)\\ {\rm{Nu}}{{\rm{m}}_z} = {\rm{floor}}\left( {\frac{{{z_{{\rm{max}}}} - {z_{{\rm{min}}}}}}{l}} \right) \end{array} \right. $

(1)

用体素内所有点集p_i的坐标平均值作为当前体素的坐标值，完成体素化：

$ \left\{ \begin{array}{l} {T_{ix}} = {\rm{floor}}\left( {\frac{{{p_{ix}} - {x_{\min }}}}{l}} \right)\\ {T_{iy}} = {\rm{floor}}\left( {\frac{{{p_{iy}} - {y_{\min }}}}{l}} \right)\\ {T_{iz}} = {\rm{floor}}\left( {\frac{{{p_{iz}} - {z_{\min }}}}{l}} \right) \end{array} \right. $

(2)

1.2 形变量提取

边坡的点云数据是单点的离散数据，不具备连续的空间特性。为了能够更好地对目标数据进行表述，本文利用种子生长法^[18]对点云数据进行分类，对满足同一空间特性的点云数据进行提取，作为形变监测的参考基准面。不同的生长结果具有不同的空间属性，为了能够更好地显示分类结果，对具有不同空间属性的点云数据赋予不同的RGB值。原始点云数据及分类结果如图 2所示。

在进行观测数据与目标数据的形变监测之前，需要对观测数据和目标数据进行配准，将2期点云数据统一在同一坐标系下。本文采用文献[15]提出的反算拼接参数误差模型提高配准精度，将多次获取的观测数据进行体素化处理，找到每个体素点坐标对应的种子生长法得到的分类面，利用体素点到分类面的距离作为当前体素的形变量。假设当前观测数据的体素点为(GP_x, GP_y, GP_z), 则依据式(2)，当前体素点在目标数据中的对应体素点为(TP_x, TP_y, TP_z)，利用体素点的空间坐标找到对应的分类面及其空间参数：

$ ax + by + cz + d = 0 $

(3)

假设点i的形变量、速率、加速率分别为ξ_i、λ_i、Ω_i，则：

$ \left\{ \begin{array}{l} {\xi _i} = \frac{{\left| {a{\rm{G}}{{\rm{P}}_{ix}} + b{\rm{G}}{{\rm{P}}_{iy}} + c{\rm{G}}{{\rm{P}}_{iz}} + d} \right|}}{{\sqrt {{a^2} + {b^2} + {c^2}} }}\\ {\lambda _i} = \frac{{{\xi _i} - {\xi _{i - 1}}}}{t}\\ {\mathit{\Omega }_i} = \frac{{{\lambda _i} - {\lambda _{i - 1}}}}{t} \end{array} \right. $

(4)

2 形变监测与分析

本文利用扩展卡尔曼滤波算法^[19]，将边坡的形变量、速率、加速率作为状态向量，通过前一期的数据结果对下一次的状态向量进行预测，完成预测更新; 利用当前获取的数据结果对预测状态向量进行修正，完成状态更新，实现对边坡形变情况的求解。

2.1 系统空间和状态向量

形变量的观测方程为：

$ \mathit{\boldsymbol{L}} = \mathit{\boldsymbol{BX}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} $

(5)

式中，L为观测值，B为已知系数矩阵; Δ为观测误差向量，即观测噪声向量。边坡形变过程可以理解为由处于运动状态的地面点构成，是关于时间t的函数。扩展卡尔曼滤波的状态方程可以表述为：

$ {\mathit{\boldsymbol{X}}_{k + 1}} = {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{k + 1,k}}\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{\xi }}_k} + {\mathit{\boldsymbol{ \boldsymbol{\varGamma} }}_{k + 1}}{\mathit{\boldsymbol{\omega }}_k} $

(6)

随机模型为：

$ \left\{ \begin{array}{l} E\left( {{\omega _k}} \right) = 0,E\left( {{\Delta _k}} \right) = 0\\ {\mathop{\rm cov}} \left( {{\omega _k},{\omega _j}} \right) = {\mathit{\boldsymbol{D}}_k}{\xi _{kj}}\\ {\mathop{\rm cov}} \left( {{\Delta _k},{\Delta _j}} \right) = {\mathit{\boldsymbol{D}}_{{\Delta _k}}}{\xi _{kj}}\\ {\mathop{\rm cov}} \left( {{\omega _k},{\omega _j}} \right) = 0\\ E\left( {{\mathit{\boldsymbol{X}}_0}} \right) = {\mathit{\boldsymbol{\mu }}_x}\left( 0 \right) = \mathit{\boldsymbol{\hat X}}\left( {0/0} \right)\\ {\mathop{\rm var}} \left( {{\mathit{\boldsymbol{X}}_0}} \right) = {\mathit{\boldsymbol{D}}_{{x_0}}}\\ {\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{X}}_0},{\omega _k}} \right) = 0,{\mathop{\rm cov}} \left( {{\mathit{\boldsymbol{X}}_0},{\Delta _k}} \right) = 0 \end{array} \right. $

(7)

在边坡形变测量中，假设地面点i在t时刻的位置向量为ξ_i(t)，即$({{x}_{i}}, {{y}_{i}}, {{z}_{i}}) $, 其沉降瞬时速率为λ_i(t)，地面沉降的瞬时加速率为Ω_i(t)，则：

$ {\dot \xi _i}\left( t \right) = {\lambda _i}\left( t \right) $

(8)

$ {\dot \lambda _i}\left( t \right) = {\mathit{\Omega }_i}\left( t \right) $

(9)

地面点状态向量的描述可认为是位置向量、瞬时速率以及瞬时加速率的函数，记为f(ξ_i(t), λ_i(t), Ω_i(t))，k时刻的状态向量为X_k(t)，所有时刻的状态向量为X(t)，则边坡形变的状态模型可以表述为：

$ {\mathit{\boldsymbol{X}}_k}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {{\xi _k}\left( t \right)}\\ {{\lambda _k}\left( t \right)}\\ {{\mathit{\Omega }_k}\left( t \right)} \end{array}} \right],\mathit{\boldsymbol{X}}\left( t \right) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{X}}_1}\left( t \right)}\\ {{\mathit{\boldsymbol{X}}_2}\left( t \right)}\\ \vdots \\ {{\mathit{\boldsymbol{X}}_k}\left( t \right)} \end{array}} \right] $

(10)

2.2 预测更新

在完成观测数据的更新时，利用上一次观测结果对当前结果进行预测，预测状态方程和观测方程为：

$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\bar X}}}_{k + 1,k}} = {\mathit{\boldsymbol{\varphi }}_{k + 1,k}}{\mathit{\boldsymbol{X}}_k} + {\mathit{\boldsymbol{\omega }}_{k + 1}} = \\ \left[ {\begin{array}{*{20}{c}} \mathit{\boldsymbol{I}}&{\Delta {t_k}\mathit{\boldsymbol{I}}}&{\frac{1}{2}\Delta {t_k}\mathit{\boldsymbol{I}}}\\ 0&\mathit{\boldsymbol{I}}&{\Delta {t_k}\mathit{\boldsymbol{I}}}\\ 0&0&\mathit{\boldsymbol{I}} \end{array}} \right]{\mathit{\boldsymbol{X}}_k} + {\mathit{\boldsymbol{\omega }}_{k + 1}}\\ {\mathit{\boldsymbol{L}}_{k + 1}} = {\mathit{\boldsymbol{B}}_{k + 1,k}}{\mathit{\boldsymbol{X}}_{k + 1}} + \Delta {\mathit{\boldsymbol{v}}_{k + 1}} \end{array} \right. $

(11)

式中，ω_k+1为状态误差，Δv_k+1为观测误差。

利用上一次观测结果的协方差矩阵对当前时刻的最优预测误差协方差矩阵进行求解：

$ {\mathit{\boldsymbol{D}}_{k + 1,k}} = {\mathit{\boldsymbol{\varphi }}_{k + 1,k}}{\mathit{\boldsymbol{D}}_k}\mathit{\boldsymbol{\varphi }}_{k + 1,k}^{\rm{T}} + {\mathit{\boldsymbol{R}}_k} $

(12)

式中，φ_{k+1, k}为状态转移矩阵，R_k为噪声误差。

2.3 状态更新

在完成状态预测后，需要对卡尔曼增益矩阵J_k进行求解，进一步修正状态向量预测值$\mathit{\boldsymbol{\bar X}}{_{k + 1, k}} $：

$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{J}}_{k + 1}} = }\\ {{\mathit{\boldsymbol{D}}_{k + 1,k}}\mathit{\boldsymbol{B}}_{k + 1}^{\rm{T}}{{\left[ {{\mathit{\boldsymbol{B}}_{k + 1}}{\mathit{\boldsymbol{D}}_{k + 1,k}}\mathit{\boldsymbol{B}}_{k + 1}^{\rm{T}} + {\mathit{\boldsymbol{D}}_{\Delta \left( {k + 1} \right)}}} \right]}^{ - 1}}} \end{array} $

(13)

更新后的状态向量及协方差矩阵为：

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{X}}_{k + 1}} = {\left[ {\mathit{\boldsymbol{D}}_{k + 1,k}^{ - 1} + \mathit{\boldsymbol{B}}_{k + 1}^{\rm{T}}\mathit{\boldsymbol{D}}_{{\Delta _{\left( {k + 1} \right)}}}^{ - 1}{\mathit{\boldsymbol{B}}_{k + 1}}} \right]^{ - 1}}\\ \left[ {\mathit{\boldsymbol{D}}_{k + 1,k}^{ - 1}{{\mathit{\boldsymbol{\bar X}}}_{k + 1,k}} + \mathit{\boldsymbol{B}}_{k + 1}^{\rm{T}}\mathit{\boldsymbol{D}}_{{\Delta _{\left( {k + 1} \right)}}}^{ - 1}{\mathit{\boldsymbol{L}}_{k + 1}}} \right]\\ {\mathit{\boldsymbol{D}}_{k + 1}} = {\left[ {\mathit{\boldsymbol{D}}_{k + 1,k}^{ - 1} + \mathit{\boldsymbol{B}}_{k + 1}^{\rm{T}}\mathit{\boldsymbol{D}}_{{\Delta _{\left( {k + 1} \right)}}}^{ - 1}{\mathit{\boldsymbol{B}}_{k + 1}}} \right]^{ - 1}} \end{array} \right. $

(14)

完成对状态向量的更新，得到具体的形变量、速率及加速率。

3 边坡实验

选取铁路旁的边坡为观测对象。该场地之前发生过滑坡，为避免再次滑坡，利用水泥进行加固。实验场地如图 3所示，红色框为实验区域。

对实验场进行为期1 a的监测，利用三维激光扫描仪每2个月进行一次数据获取，以首次观测作为目标数据。假定目标数据稳定且不存在形变，针对获取的数据进行体素化处理并分类(图 2)，获取每一类点云数据的空间参数。每隔一定时间对场地进行观测，求解场地的形变量。每次获取的数据都是以扫描仪为坐标中心的局部三维坐标系，不需要严格对中整平，只需仪器所在位置大致统一即可。

选择实验场地中稳定建筑物的点云数据进行配准，并通过人工选点的方式进行精度检验，保证误差控制在0.5 mm以内，如图 4所示。为了视觉效果更好，基础数据和观测数据分别用不同的颜色显示，其中白色为基础数据，绿色为观测数据。

完成配准后，按照本文的方法对整个场内的点云数据进行处理，得到每个体素点经过滤波后的形变量、速率、加速率。本次实验共进行7次数据获取，其中第一次为基础数据，其余6次为观测数据，如图 5所示。

在进行数据获取时，碍于视角的影响，实验场的点云数据存在空白部分，并不影响整体观测的结果。结合监测区域数据处理结果分析，第1~3次观测结果显示，边坡整体不存在明显形变量，局部变化为场地误差。第4次观测结果显示，右上角和左下角存在一定的区域性形变。结合场地影像发现，左下角是由于场地加固部分生长出一部分杂草，与边坡形变无关; 右上角存在约3 cm的形变量。第5次观测结果显示，去除杂草的影响，场中右上角仍存在约5~7 cm的形变量，需要引起重视。第6次监测结果显示，场地右上角发生15 cm的形变，结合场地影像(图 6)发现，实验场地右上角土坡发生滑落。

4 结语

利用三维激光扫描仪对边坡实验场地进行数据获取，采用本文提出的处理模式，即数据体素化预处理-体素点分类-分类结果拟合-形变量求解，对实验场地的6次观测数据进行处理分析。结果显示，在第4次观测中可以明显看出场地右上角存在的形变量，第5次观测中形变量继续扩大，直至第6次发生滑坡。证明在边坡形变监测中，利用非接触式的三维激光扫描仪对数据进行获取，能够在保障人身安全的同时对实验区域内的全部场景进行有效的形变量解算，从而对可能存在滑坡的情况进行监测、预警、预报。