2. 中国科学院大学,北京市玉泉路甲19号,100049
滑坡是非常严重的自然灾害,它的触发因素有很多,例如降雨、地下水位变化、地震火山活动以及人类活动等。目前有很多技术可以监测滑坡形变,主要分成两类:一是传统的点测量方法[1],主要有全站仪、全球定位系统(GPS)等;二是最新的面监测技术,有合成孔径雷达干涉测量(interferometric synthetic aperture radar, InSAR)[2]以及地面三维激光扫描测量(terrestrial laser scanner,TLS)[3]等。地面三维激光扫描(TLS)可以快速(每秒百万个点)、高空间分辨率(可达1 mm一个点)、高精度(可达1~2 mm)[4]地获取地表三维信息,这些信息以点云形式保存下来。TLS在国内外应用非常广泛。国外学者Masi[5]将其用于古文物数字化,Fanti[6]将TLS技术应用到岩滑当中,Travelletti[7]、Pesci[8]、Dunning[9]等也成功地将TLS技术应用到滑坡监测当中。国内起步较晚,但也有很多成功案例,例如邱俊玲[10]将TLS应用到矿山边坡监测,刘绍堂[11]利用TLS来监测隧道形变,徐进军[12]成功地利用TLS来监测分析桥梁的变形情况。滑坡区域的地形表面在滑动之后会发生变化,可以通过由TLS获取地表点云生成高精度的数字表面模型(DSM),对多期DSM进行比较,得出形变随时间的变化,并且可以估计位移量、表面变化、体积变化。本文以四川理县滑坡为例,探讨了点云数据获取、数据预处理流程的关键技术,给出了一套非常全面、直观的滑坡形变分析方法。为了充分利用非地面点云数据,提出了一种新方法来分析滑坡体的水平变化。
1 三维激光扫描测量原理和扫描精度分析 1.1 三维激光扫描测量原理本文以FARO Laser Scanner Focus 3D系列扫描仪为例,介绍三维激光扫描的测量原理。首先,扫描仪内部的激光发射装置以每秒百万激光束的速度向目标表面发射激光脉冲,再由扫描仪内部脉冲接收装置接受反射回来的激光脉冲,同时记录脉冲传播时间、脉冲横向扫描角度α和纵向扫描角度β。根据脉冲在介质中的传播速度和α、β便可计算出目标点P在仪器自定义坐标系下的坐标,如图 1所示。则目标点P的坐标计算公式为:
(1) |
虽然扫描仪是面测量方式,但其初始获得的依然是单点坐标。因此,我们可将扫描点空间位置的计算公式线性化,从而推导出扫描点空间位置的理论精度。
对式(1)全微分,获得系数矩阵:
根据协方差传播定律,便可获得三维激光扫描点空间位置精度的理论值为:
(2) |
式中, σα、σβ、σS由仪器厂家给定,故单个扫描点的空间位置理论精度最终由扫描点P到扫描仪的距离S和纵向扫描角度β共同决定。此外,物体本身的反射率、外界环境对点云数据精度的影响同样很大,应尽量避免发射率较低及雨、雪、沙尘暴等异常天气的情况。
2 点云数据采集本文以四川理县一地质滑坡为例采集点云。滑坡体为土质,滑坡区域大约为200 m×200 m,且滑坡体前沿为当地小学教学楼,因此滑坡监测和防护措施尤为重要。数据采集时间分别为2015-01-23和2015-12-28,采用美国FARO公司生产的FARO Laser Scanner Focus 3D X 120三维激光扫描仪,测量精度可达±2 mm/25 m,扫描速度为97.6万点/s,视角范围为水平360°、垂直270°,最高分辨率为10 m处点间距0.9 mm。扫描过程如下:
1) 布置球形标靶。首先对滑坡体进行踏勘,确定扫描站数和仪器放置位置,并合理设计球形标靶的放置位置,然后布置球形标靶。为了保证相邻测站点云能够利用布尔沙坐标转换模型拼接,两测站至少应有3个以上可通视的重合标靶。
2) 扫描滑坡体。将扫描仪安置在测站上,调整仪器的姿态。然后打开扫描仪,设置扫描质量和分辨率等开始扫描。扫描质量和分辨率越高,扫描时间越长,可根据实际精度要求合理选择。如图 2,三角形为扫描仪位置,圆形为标靶位置,线条内为滑坡区域。
点云过滤主要目的有两个。一是删除与地形无关的建筑物、树木、人等点云,以及雾、施工浮尘等空气中杂质反射而成的离群噪声点。对于与建模或者地形无关的建筑物、树木、人等点云,采取目测手动的方式剔除,而对于离群噪声点,采用“离群”方法自动删除。其基本思想如下:首先,设定栅格尺寸、距离阈值、分配阈值的大小;然后,获取在由栅格尺寸设定的周围区域内的有效点,并计算其中所有点到扫描仪的距离与该扫描点到扫描仪的距离的差值;最后,计算差值小于距离阈值的扫描点的百分比,如果百分比也在分配阈值内,则保留此扫描点,反之则删除。二是降低点云密度,提高数据处理效率。高斯滤波在指定区域内权重分布为高斯分布,其滤波之后的结果对原始数据损害较小,能够很好地保持原始数据的形态,故在本文选择高斯滤波器来降低点云密度(图 3)。
对于一个范围较大的滑坡而言,由于受到视线遮挡、扫描仪测量距离有限等影响,需要在不同测站从不同视角对滑坡体进行扫描,从而得到了不同仪器坐标系的点云数据。为了构建滑坡体的整体模型,就必须将不同测站的点云转换到同一坐标系下,这一过程叫作点云拼接。点云拼接的方法主要分为两种:基于几何特征的拼接和基于标靶的拼接。由于在本滑坡场景中几何特征不明显,本文采取基于标靶的拼接方法。其拼接步骤如下:首先,利用拟合法获取至少3个以上球形标靶的球心坐标;然后,根据布尔沙坐标转换模型[13-14],按照最小二乘法求解转换参数;最后,根据转换参数,将不同视角的点云转换到同一坐标系下。其流程如图 4所示,拼接结果见图 5。
将两期或者多期点云转换至同一绝对坐标系下叫作地理参考。由于两期布置的标靶位置很难一致,而且为了方便之后的形变分析,需要将拼接之后的点云转换到同一绝对坐标系下[15]。首先利用GPS获取标靶中心的GPS坐标,并拟合出至少3个以上标靶中心在扫描仪坐标系下的坐标;然后利用布尔沙坐标转换模型根据最小二乘原理求解转换参数;最后将扫描仪坐标系下的两期或多期点云转换到绝对坐标系下。
4 滑坡形变分析 4.1 基于点云比较的形变分析利用单点进行滑坡监测的准确性相对较高,但只能监测某一部位的形变情况,无法监测整个区域的形变。而利用两期点云进行滑坡形变分析,能够全方位地反映滑坡体的整体变形以及滑坡体各个部位的变形大小。基于点云比较的滑坡整体形变分析步骤如下:1)读取已经经过点云预处理的两期点云20150123.xyz和20151218.xyz;2)对第一期点云进行封装,将点云数据转换成多边形,并对多边形进行修复,消除钉状物、空洞等缺陷多边形;3)以第一期点云生成的多边形为参考,将第二期的点云与之比较,自定义一个偏差方向,生成两期观测值在该方向上的差值,即形变值。由图 6(单位m)可知,两期点云在Z方向的变形量在2 dm左右,这与GNSS实时远程变形监测结果基本一致。
利用两期点云生成的不规则三角网(TIN)进行滑坡形变分析,也能够全方位反映滑坡体的整体变形以及滑坡体各个部位的变形大小。其步骤如下:首先读取两期点云,根据Delaunay法则生成两期点云的不规则三角网TIN;然后以第1期TIN为参考,将第2期TIN与之比较,自定义一个偏差方向,生成两期TIN在该方向上的差值,即形变量。由图 7(单位m)可知,TIN比较结果与点云比较结果基本一致。
扫描仪获取的点云数据中包括地面数据和非地面数据。利用地面数据可以分析出整个滑坡体的形变。除此之外还有大量的非地面数据,如电线杆、树干等,这些地物有个共同点就是不易发生移动,因此可以提取这些地物的几何特征并应用到多期观测的同名点比较,进而可获得一些特征点的形变值。为了充分利用非地面数据,本文提出一种利用滑坡区域的电线杆或者树干来分析滑坡体水平位移的方法。其基本原理如下:首先在固定高度处截取电线杆的一个截面(图 8(a)),然后利用最小二乘拟合出截面的圆心坐标(图 8(b)),最后利用两期同一高度截面的圆心坐标得出电线杆所在滑坡点位的水平位移(图 8(c))。此外,也可以在滑坡的关键位置上放置一些固定标靶,通过拟合其几何中心来分析滑坡的形变(表 1)。
另外,在固定标靶处安装有GNSS实时远程变形监测系统,其在2015-01-23~2015-12-28期间的连续监测结果如图 9。通过对比可以看出,基于固定标靶的形变分析结果与实时GNSS实时监测结果基本一致。
1) Kriging插值法生成两期DEM。基于Kriging插值法能获得高精度DEM,并且基于DEM的土方量计算是当前的主要方法之一[16]。Kriging插值法与其他插值法的一个重要不同之处在于Kriging插值法不仅考虑了待插值点与邻近点之间的空间位置关系,也考虑了邻近点之间的位置关系,内插出的结果更符合真实情况。为了便于分析两期观测的滑坡形变,本文以1 dm的栅格大小,利用Kriging插值法生成两期DEM。
2) 滑坡土方量分析。将相同大小的两期DEM求差便可得到差值栅格文件,利用差值栅格文件,计算出2015-01-23~2015-12-18滑坡的沉降土方量约为7 034 m3。结合上文生成的3D比较结果,可以非常直观地得出滑坡引起的体积变化和体积变化量级在整个滑坡体的分布情况,这大大提高了滑坡灾害评价和治理工作的效率。
3) 滑坡剖面分析。通过事先设计好的剖面线矢量文件和已经生成的两期DEM,可以提取出滑坡区域的多条剖面线或者等高线。本文利用前后两次扫描生成的DEM和同一剖面线矢量文件,绘出滑坡前后的剖面,通过比较,可以得出滑坡的变形趋势(图 10)。
本文利用地面三维激光扫描仪获取了四川理县一地质滑坡的点云数据,探讨了点云数据预处理的关键技术。通过点云比较、TIN比较、特征点比较、剖面线比较分析了滑坡的形变量和变形趋势。为了充分利用非地面点云,本文提出一种提取滑坡区域树干或者电线杆等地物的几何特征来分析坡体水平变形的方法。实验结果表明,三维激光扫描仪以其独特的面测量方式,能够获取滑坡体表面丰富的点云信息。利用这些点云可以全面而又直观地分析滑坡体的表面变化、体积变化和变形趋势,为地质灾害防护提供宝贵的数据资料。
[1] |
Kasperski J, Delacourt C. Application of a Terrestrial Laser Scanner (TLS) to the Study of the Séchilienne Landslide (Isère, France)[J]. Remote Sensing, 2010(22): 2785-2802
(0) |
[2] |
Squarzoni C, Delacourt C, Allemand P. Nine Years of Spatial and Temporal Evolution of the La Valette Landslide Observed by SAR Interferometry[J]. Engineering Geology, 2003, 68: 53-66 DOI:10.1016/S0013-7952(02)00198-9
(0) |
[3] |
Lemmens M. Geo-Information[M]. Netherlands: Springer Netherlands, 2011
(0) |
[4] |
Shan J, Charles K. Topographic Laser Ranging and Scanning: Principles and Processing[M]. London: CRC Press, Taylor & Francis Group, LLC, 2008
(0) |
[5] |
Masi A D. Engineering Geology for Society and Territory-Volume 8[M]. Switzerland: Springer International Publishing, 2014
(0) |
[6] |
Fanti R, Gigli G. Terrestrial Laser Scanning for Rockfall Stability Analysis in the Cultural Heritage Site of Pitigliano (Italy)[J]. Landslides, 2013, 10: 409-420 DOI:10.1007/s10346-012-0329-5
(0) |
[7] |
Travelletti J, Malet J P, Delacourtc C. Image-Based Correlation of Laser Scanning Point Cloud Time Series for Landslide Monitoring[J]. International Journal of Applied Earth Observation and Geoinformation, 2014, 32: 1-18 DOI:10.1016/j.jag.2014.03.022
(0) |
[8] |
Pesci A, Teza G. Multitemporal Laser Scanner-Based Observation of the Mt. Vesuvius Crater: Characterization of Overall Geometry and Recognition of Landslide Events[J]. ISPRS Journal of Photogrammetry and Remote Sensing, 2011, 66: 327-336 DOI:10.1016/j.isprsjprs.2010.12.002
(0) |
[9] |
Dunning S A, Massey C I, Rosser N J. Structural and Geomorphological Features of Landslides in the Bhutan Himalaya Derived from Terrestrial Laser Scanning[J]. Geomorphology, 2009, 103: 17-29 DOI:10.1016/j.geomorph.2008.04.013
(0) |
[10] |
邱俊玲, 夏庆霖, 姚凌青, 等. 基于三维激光扫描技术的矿山地质建模与应用[J]. 地球科学, 2012, 37(6): 1209-1215 (Qiu Junling, Xia Qinglin, Yao Lingqing, et al. Mine Geological Modeling and Application Based on the Three-Dimensional Laser Scanner Technology[J]. Earth Science, 2012, 37(6): 1209-1215)
(0) |
[11] |
刘绍堂, 刘文锴, 周跃寅. 一种隧道整体变形监测方法及其应用[J]. 武汉大学学报:信息科学版, 2014, 39(8): 981-986 (Liu Shaotang, Liu Wenkai, Zhou Yueyin. A Tunnel Overall Deformation Monitoring Method and Its Application[J]. Geomatics and Information Science of Wuhan University, 2014, 39(8): 981-986)
(0) |
[12] |
徐进军, 廖骅, 韩达光, 等. 大跨度桥梁桥面线形测量新方法[J]. 测绘通报, 2016(1): 91-94 (Xu Jinjun, Liao Hua, Han Daguang, et al. New Methods for Deflection Measurement of Bridge with Large Span[J]. Bulletin of Surveying and Mapping, 2016(1): 91-94)
(0) |
[13] |
段鹏硕, 刘根友, 龚有亮, 等. 空间坐标系变换的函数梯度描述方法[J]. 测绘学报, 2014, 43(10): 1005-1012 (Duan Pengshuo, Liu Genyou, Gong Youliang, et al. The Functional Gradient Description Method of Space Coordinate Transformation[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(10): 1005-1012)
(0) |
[14] |
杨开伟, 李娟娟, 甘兴利, 等. 基于布尔萨模型的大旋转角坐标转换参数算法[J]. 全球定位系统, 2011(6): 45-49 (Yang Kaiwei, Li Juanjuan, Gan Xingli, et al. Arithmetic of Coordination Transformation Parameters of Big Rotation Angle Based on Bursa Model[J]. GNSS World of China, 2011(6): 45-49 DOI:10.3969/j.issn.1008-9268.2011.06.013)
(0) |
[15] |
Olsen M J, Johnstone E, Driscoll N, et al. Terrestrial Laser Scanning of Extended Cliff Sections in Dynamic Environments: Parameter Analysis[J]. Journal of Surveying Engineering, 2009, 135(4): 161-169 DOI:10.1061/(ASCE)0733-9453(2009)135:4(161)
(0) |
[16] |
曹俊茹, 刘强, 姚吉利, 等. 基于Kriging插值DEM的计算土方量方法的研究[J]. 测绘科学, 2011, 36(3): 98-99 (Cao Junru, Liu Qiang, Yao Jili, et al. Method of Calculating Earthwork Based on the DEM of Kriging Interpolation[J]. Science of Surveying and Mapping, 2011, 36(3): 98-99)
(0) |
2. University of Chinese Academy of Sciences, A19 Yuquan Road, Beijing 100049, China