基于三维激光扫描的矿区建筑物变形监测与采动损害评估

王 磊，朱尚军，2,蒋 创,江克贵,郭庆彪

(1.安徽理工大学 空间信息与测绘工程学院，安徽 淮南 232001;2.安徽省地质测绘技术院，安徽 合肥 230000)

0 引 言

在煤炭资源地下开采过程中，岩层移动随着工作面的推进与开采强度的提高从下向上逐渐发展到地表，并表现为地表下沉、倾斜、曲率、水平移动与水平变形等移动变形。上述移动变形将导致各种地质环境灾害，如滑坡、地表塌陷，房屋桥梁垮塌和地表裂缝等[1-4]，对开采沉陷盆地内的建筑物、铁路、公路和高压线塔造成严重影响。因此，分析煤炭开采过程中地表建筑物的稳定性并对建筑物的损坏等级进行评定与防治对矿区的生产建设及安全有着重要的意义[5-7]。

常见的测量及评估矿区地表建筑物破坏程度方法有：①钢尺量距法，直接量取房屋的裂缝状态；②在建筑物周围及内部布设观测站，利用高精度的水准仪与全站仪量取房屋的下沉与水平移动，经过后续数据的处理与分析评估房屋的损坏等级[8-10]。然而，上述2种传统方法获得的数据都是离散型、低密度的建筑物局部信息，不能完整有效地对建筑物损坏等级进行评价。

三维激光扫描技术(3DLS)作为一种新型测量技术，具有自动化、连续性、非接触式、高密度等优势弥补了传统技术监测矿区沉陷的不足。三维激光扫描仪通过连续布置测站，可以获得建筑物整体的、连续的点云数据，经过点云数据的预处理、配准后可以连续地、完整地描述建筑物的倾斜与水平变形，达到能够细致准确地评估煤矿区地表建筑物损坏等级的目的[11-13]。戴华阳等[14]提出利用三维激光扫描技术提取采动区房屋特征点信息获得采动区地表建筑物形变量，并经过误差分析理论，评价直接获取数据与间接计算数据的精度，得出其精度结果能满足采动区房屋移动变形监测标准。吴侃等[15]分析了三维激光扫描技术的单点定位的精度，数据采集及处理方法，提出利用建筑物的特征线来判断房屋是否变形的思路。李得军等[16]提出三维激光扫描技术提取高压线塔的变形量，并利用最小二乘平面拟合方法分析了其精度。张豪杰[17]研究了三维激光扫描技术应用于采动区高等级公路下变形，提出了一种基于遗传算法的提取圆台中心线的算法，为高等级公路下压煤开采提供了理论与数据的支持。

高精度的三维激光技术应用于矿区建筑物变形监测不仅依赖于精密的模型和误差分析理论，同时有效的标靶球拟合方法、合理的矿区建筑物变形监测设计方法对采动损害评估也有着至关重要的意义；鲁铁定等[18]提出了一种地面激光扫描标靶球定位的整体最小二乘定位方法，并推导出和比较2种整体最小二乘解算法；陶武勇等[19]提出求解球面拟合的改进总体最小二乘法，该算法考虑到观测向量中的非线性形式，得到了更加准确地拟合球心。通过上述文献发现，现阶段的标靶球拟合方法在求取球心坐标过程中存在计算量大、线性化过程复杂等问题；异于传统的矿区建筑物变形监测设计方法，三维激光技术实地观测存在设计困难、精度难以控制等问题，难以有效地对于矿区的建筑物进行采动损害评估。

针对上述问题，首先介绍了基于三维激光的建筑物变形监测精度设计方法，建筑物移动与变形计算与采动损害评估方法，其次利用QPSO算法对标靶球球心拟合精度进行了分析，将多个三维激光扫描仪测站转换到同一坐标系统下，并根据三维激光扫描仪合理的观测站布设方案对淮南某矿区地表建筑物进行了实时观测，最后基于三维激光扫描技术提取了开采沉陷移动变形盆地内地表建筑物的移动变形值，并总结了矿区地表建筑物的损坏等级规律分布。为附近矿区的地质环境灾害评估与防治提供了可靠的依据。

1 基于3DLS的建筑物变形监测与采动损害评估原理

1.1 三维激光扫描基本原理

三维激光扫描仪在扫描获取点云信息的过程中，是利用测距的原理，仪器依据接收回来的脉冲信号来识别物体信息，三维激光扫描仪架设位置确定，则该测站的坐标轴即确定，扫描仪架设点为三维坐标原点，横向扫描为X轴，Y轴所确定的平面，竖向扫描物体包含Z轴，Z轴方向向上为正。其三维坐标系统如图1所示。

图1 三维激光扫描仪坐标系原理

由图1可知，点云中P点的坐标计算公式为：

(1)

式中：S为P点至坐标原点O点的距离，为距离观测值;α为纵向扫描角度;β为横向扫描角度。

1.2 基于3DLS建筑物变形监测精度设计方法

对于一个建成的建筑物变形观测网，变形观测点所能够达到的精度，除取决于监测网对观测点的测量精度，也受监测网本身的精度影响；三维激光扫描仪对于矿区建筑物变形监测也需确定建筑物变形测量等级及精度，表1为建筑物变形测量等级及精度。三维激光扫描仪直接获取观测点的三维坐标值，依据三维激光扫描仪测量原理，根据文献[21],依据式(1)推导得激光点位坐标Xi=[X,Y,Z]的协方差为

表1 建筑物变形测量等级及精度[20]

(2)

(3)

(4)

(5)

式中：σS为测距误差；σα、σβ为纵向、横向扫描角度误差。

通过以上模型,依据表1建筑物变形测量等级及精度要求，通过误差传播定律可求三维激光扫描仪的最大扫描测量距离，并由此可设计出实地勘察矿区建筑物三维激光扫描仪用于变形监测的精度设计方案；在监测矿区建筑物形变时会出现多站观测需进行点云数据拼接，点云数据进行拼接时也存在误差，故笔者将每一站测量控制在四等变形测量等级，直接将多个观测站进行坐标转换，将各观测站数据统一到已建立的坐标系统中，避免了点云数据的拼接误差。

1.3 建筑物移动与变形计算与采动损害评估方法

笔者采用手动选取同名特征点。在三维激光扫描观测过程中，加大了特征点区域的点云密度，并且在计算地表建筑物的形变过程中多次计算取均值。以房屋的角点作为特征点，计算了特征点n+1至特征点n的倾斜值i′、水平变形ε，用此方法沿着建筑物的长轴方向计算了若干房屋的倾斜值i′、水平变形ε，分析各个房屋的倾斜值i′、水平变形ε得出结论。

各个房屋2个特征点之间的倾斜值i′、水平变形ε，用以下公式获得：

相邻2点的倾斜为:

(6)

式中:ln～n+1为n号点至n+1号点的水平距离，可以由直接观测所得，也可以由两点的坐标反算求得；Wn+1、Wn为n+1号点、n号点的下沉值，mm。

其中正负号为在倾斜断面上，指向上山方向为正，指向下山方向为负；走向断面上，指向右侧为正，指向左侧为负。

n～n+1号点之间的水平变形为:

(7)

式中：(ln+1-n)0、(ln+1-n)m为n+1～n号点在首次观测时和m次观测时的水平距离，可以由直接观测所得，也可以由两点的坐标反算求得。

地表曲率Km～n～p是两相邻线段的倾斜差与两线段中点间的水平距离的比值，它反映了观测线断面上的弯曲程度，计算公式如下：

(8)

依据建筑物下、铁路下和水体下采煤规范并结合上述建筑物移动与变形计算方法对房屋的损坏等级进行分析，由于实地观测房屋大部分为砖混结构，因此本文实测矿区的地表建筑物满足表2所设定。

表2 砖混结构建筑物损坏等级[22]

2 基于QPSO的三维激光标靶球中心拟合方法

三维激光点云数据预处理过程中包括坐标转换与点云数据拼接，在坐标转换过程中需要利用相对坐标系下的标靶球球心坐标与北京54坐标系下的标靶球球心坐标进行点云数据的整体坐标转换；点云数据拼接过程中高精度的标靶球球心坐标有助于提高点云数据拼接的精度，因此得到高精度的相对坐标系下的标靶球球心坐标对点云数据的预处理过程尤为重要。

2.1 方法构建

根据球心方程：

(9)

式中:(xm,ym,zm)为标靶球球心拟合数据点m坐标；(x0,y0,z0)为标靶球球心坐标。

由式(9)知该方程是多参数非线性化方程，目前，标靶球球心拟合的方法一般有最小二乘法、总体最小二乘法、混合总体最小二乘法等[18-19]，上述方法在求取球心坐标过程中存在计算量大、线性化过程复杂等问题，文献[23-24]表明QPSO算法可以弥补这些不足，并且在水利领域与计算机算法领域都有着实际应用价值，而在标靶球球心拟合方面取得的应用成果较少。故此本节讨论了基于QPSO算法的标靶球拟合过程及精度分析。

R为标靶球半径，则该点到球心的残差Vm为

vm=vm∑|Rm-R|

(10)

根据误差绝对值之和最小原则，构建出球心求参准则如下：

(11)

式中:M为测点数；可采用量子粒子群算法对式(11)求解球心坐标参数。

上文所述，基于MATLAB编码实现量子粒子群优化算法，其流程总结如下：

在解空间有1 000个粒子，每个粒子有三维,在第t次迭代中粒子的当前位置用xi3=[xi1(t),xi2(t),xi3],i=1,2,…,1 000表示；粒子的历史最优位置用pi3(t)=[pi1(t),pi2(t),pi3(t)]表示，群体最优用pg3(t)=[pg1(t),pg2(t),pg3(t)]表示；将球面上加入误差的点求得的球半径Rm与真实半径R之差的绝对值累加作为本文的适应度函数f。公式如下：

f=∑|Rm-R|

(12)

适应度函数越小说明粒子位置越好，即求得球心坐标结果越接近于真实值。

具体操作过程如下：

1)初始化种群，将粒子当前位置初始化为个体历史最优位置，计算适应度找到群体最优位置。

2)确定粒子i(1≤i≤1 000)的介于个体历史最优与群体最优之间的位置pgi。

3)计算所有粒子个体历史最优位置的均值Mbest。

4)更新粒子的位置。

5)计算当前迭代下的粒子适应度值与前一次迭代比较，如果适应度值小，将粒子历史最优位置更换为当前粒子的位置，否则不变；找到群体最优位置与前一次迭代比较，适应度值小则进行替换，否则不变。

6)重复步骤2)～5),若循环达到最大迭代次数或满足精度，跳出循环，输出最终所得的3个球心参数。

基于量子粒子群优化算法的球心坐标求取程序实现如图2所示：

图2 QPSO标靶球中心拟合程序流程

2.2 模拟试验

为高度模拟标靶球球面点的情况，选择与标靶球相同半径R=0.075 m的模拟球面，将在球面上随机生成1 000个点即模拟点，将每个模拟点三维坐标(x,y,z)方向上均加入随机误差v～N(0，0.003)，v～N(0，0.005)，v～N(0，0.007)(v为随机误差，N为正态分布，单位m)，分3次模拟，误差依次增大，QPSO算法中每次模拟粒子数目皆为 1 000，迭代100次，3次模拟球球心坐标皆为(0.1,0.1,0.1)，每次误差条件下算法各运行10次，取均值，其结果如下：

考虑到球面上点的密度也会对标靶球球心的拟合精度有影响，故此将每个模拟点三维坐标(x,y,z)方向上随机误差控制为v～N(0，0.007)，将球面点数量降至为600、300、100；再次讨论标靶球球心拟合结果精度。

从表3可看出：① 每个随机误差条件下求得的球心坐标取均值后，均值的相对中误差在三种随机误差条件下的数值皆小于0.7%；② 从各个误差条件下10次试验结果的参数拟合中误差可看出随着随机误差的提高，x,y,z的参数拟合中误差没有发生波动性改变，参数拟合中误差最大为0.7 mm；③ 从点位中误差可看出随着随机误差大幅度地增加，球心均值点位中误差增加幅度较小且球心的点位中误差皆小于1 mm。

表3 同点云密度、不同方差条件下标靶球球心拟合结果精度分析

从表4可看出随着点云密度的降低，球心坐标均值偏差有所增加但各个轴向误差皆小于3 mm，球心均值点位中误差皆小于4 mm。图3为标靶球球面点云图。

表4 不同点云密度、同方差条件下标靶球球心拟合结果精度分析

图3 标靶球球面点云图

3 工程实例

3.1 矿区建筑物变形监测概况

淮南顾北矿南—煤采区1312(1)是顾北矿南的首采工作面，整个工作面采用后退式采煤,机械化掘进，全部垮落法管理顶板。 1312(1)工作面为正在开采中的工作面，回采速度3.4 m/d。工作面沿煤层走向布置，工作面自2019-07-10开采，至2019-11-02，工作面运输巷退尺长430.9 m，轨道巷退尺长429.6 m，平均退尺430.3 m，宽度204.6 m，工作面尺寸满足D1/H=0.396且小于1.2～1.4，D3/H=0.833且小于1.2～1.4条件(D1为工作面走向长度；D3为工作面倾向长度；H为埋深)，故工作面走向、倾向皆为非充分采动，总体为非充分采动。平均采高为4.29 m；煤层倾角平均为5°,为近水平煤层。工作面平均深度为518 m。本次试验使用的三维激光扫描仪为中海达HS-650扫描仪，测距方式采用脉冲式，全波形测量技术，可多次回波输出，最大测距650 m，最小测距1.5 m，100 m测距精度可达5 mm，扫描频率室内最大300 kHz，室外最大100 kHz，垂直角与水平角分辨率皆为0.001°,水平视场角0°～360°，垂直视场角-40°～60°，可全天候无差别进行测量试验。本次三维激光实地监测房屋A、B、C为砖混结构与工作面走向呈平行分布，房屋A、B、C东墙面距离开切眼依次为26、45、69 m；距离房屋各个特征点20 m处架设三维激光观测点，共布置了3个观测站，各站之间距离35 m。

3.2 三维激光扫描仪观测方案设计

利用三维激光扫描仪对矿区建筑物监测的方法如下;沿着房屋周围设立多个观测站，每一期每一站测量时在扫描仪周围利用移动站测量出四个当地北京54坐标系统下的点坐标信息(矿上测量时有4个标靶球，故测量出4个点的北京54坐标)，该4个点坐标用于坐标转换，在4个已知点上利用光学对中将标靶球球心与已知点对应，测量出仪器高，将已知点高程加上仪器高可得到标靶球的球心在当地北京54坐标系统下的坐标信息，通过观测的点云信息利用上述的QPSO算法进行标靶球拟合求得相对坐标系下的标靶球球心坐标，将所有的点云信息转换到统一的坐标系统下，即北京54坐标系统；不需要设立公共标靶点，直接将多个观测站进行坐标转换，将各观测站数据统一到已建立的坐标系统中，直接通过点云数据分析提取出矿区建筑物的下沉、倾斜、曲率、水平移动及水平变形，最后依据建筑物下、铁路下和水体下采煤规范对矿区建筑物的损坏等级进行合理的评估。

由于实地勘察矿区建筑物三维激光扫描仪用于变形监测满足四等变形测量等级即可，并且HS-650三维激光扫描仪垂直视场角-40°～60°，扫描测量距离太近会导致建筑物扫描不完整，故此在矿区实地观测时具体操作如下：

1)实地勘察确定建筑物变形测量等级及精度，建立建筑物变形监测精度设计方法、建筑物移动与变形计算与采动损害评估方法。

2)三维激光扫描仪观测方案设计，三维激光观测站与房屋特征点距离控制在15～40 m；观测站之间距离控制在35 m。

3)数据采集，不同监测时段三维激光观测站的观测对象要相同；由于将比较不同时期地表同一建筑物相同特征点位置的变形量，故此要求每次测量对象位置要相同。

4)数据预处理，利用QPSO算法求得三维激光标靶球球心坐标，将所有的点云信息转换到统一的坐标系统下。

5)建筑物变形提取及采动损害评估。

6)总结矿区建筑物的损坏程度与开切眼的距离及工作面开采强度的关系。

图4为三维激光扫描仪观测方案设计流程。图5为2019-07-09的矿区建筑物点云信息与工作面井上下对照图。图5中红色点云为矿区实测数据经过数据预处理，及坐标转换后的点云分布情况，从图中可看出数据处理结果与实地观测区域相符合，部分点云不呈现圆周状态是由于观测区域有房屋与树木的阻挡。图6为数据预处理过程中提取矿区建筑物的倾斜与水平移动房屋B的三维点云图。

图4 三维激光扫描仪观测方案设计流程

图5 矿区建筑物点云信息与工作面井上下对照

图6 矿区房屋B三维激光扫描点云图

3.3 建筑物变形提取与采动损害评估

基于上述矿区概况沿工作面走向平行线设立多个三维激光扫描仪观测站，对走向线上建筑物进行了多时段观测，2019-07-09第一次观测值作为首次观测值(工作面未进行开采)，随后每隔固定时间观测矿区地表建筑物。

沿着建筑物的长轴方向依据式(6)—式(7)计算多时段房屋的倾斜值i、水平变形ε见表5,曲率变形远小于0.1 mm/m2，因此未提取房屋的曲率变形。

表5中房屋A、B、C沿平行于走向线分布，房屋A最靠近开切眼，房屋C离开切眼最远；从表5中可看出房屋的损坏程度与开切眼的距离有关，房屋的损坏程度与其距离开切眼的距离成正比，距离开切眼越近房屋的损坏程度越高；房屋的损坏程度与工作面的开采进度也有关联，工作面开采强度越大，房屋的损坏程度越明显。

表5 矿区建筑物倾斜、水平变形提取

从2019-08-02变形差值可看出房屋A、B皆达到了Ⅳ级损坏程度，房屋C达到了Ⅲ级损坏程度，随着工作面的掘进，房屋C的倾斜值i、水平变形ε皆增大，故此需加强矿区地表建筑物的搬迁和维护工作。图7为矿区实拍地表房屋C采动损害图。

图7 矿区地表房屋C采动损害实拍

4 结 论

1)构建了基于QPSO的标靶球拟合算法，模拟试验结果：①同点云密度不同方差条件下，球心坐标取均值后，各轴向均值相对中误差数值皆小于0.7%;由均值拟合中误差分析可得，随着随机误差的提高，X、Y、Z的参数拟合中误差没有发生波动性改变，参数拟合中误差最大为0.7 mm，球心的点位中误差皆小于1 mm。②不同点云密度同方差条件下,球心坐标均值偏差有所增加，但各个轴向误差皆小于3 mm，球心均值点位中误差皆小于4 mm。试验结果表明QPSO算法应用于标靶球拟合有助于提高坐标转换精度和监测的结果准确性。

2)介绍了基于三维激光的建筑物变形监测精度设计方法及采动损害评估方法；开展了基于三维激光扫描的矿区建筑物变形监测及采动损害案例研究。案例结果：地表建筑物的损坏等级与距离开切眼的距离成正比，距离开切眼越近房屋的损坏等级越高，同时地表建筑物的损坏等级与工作面的掘进进度有关，随着工作面开采强度的提高地表建筑物的损坏等级也随之提高，从2019-08-02变形和2019-10-25变形数据对比上可看出随着工作面的掘进，房屋A、B、C的倾斜、水平变形皆在增加。

3)提出的基于三维激光扫描技术的矿区建筑物变形监测与采动损害评估方法对矿区地表建筑物的形变监测与采动损害评估有着重要的意义。