肖春蕾， 郭兆成， 张宗贵， 李迁， 尚博譞， 吴芳

(1.中国国土资源航空物探遥感中心，北京 100083; 2.国土资源部航空地球物理与遥感地质重点实验室，北京 100083)

0 引言

地裂缝是内、外地质营力或人类活动综合作用形成的具有一定长度、宽度和深度的地面裂隙[1]。地裂缝不仅对各类工程建筑、交通设施及土地资源造成破坏，而且会导致一系列的生态环境问题。地裂缝的形成机理主要有构造成因(如美国亚利桑那州南部的地震活动导致了已有断裂破裂面的重新复活)、地下水(或地下矿产)开采过量成因及构造与地下水(或地下矿产)开采复合成因等3种[2]。

近年来，遥感手段逐步被应用于地裂缝调查及其变化速率、微地貌变化的研究。利用差分干涉合成孔径雷达(InSAR)能够获取地表详细的位移信息，但植被变化会造成信号的去相关，使其对植被覆盖区活动地裂缝的识别与监测效果不佳[3]。同样，在植被高覆盖区使用高分辨率的光学图像提取地裂缝的分布非常困难[4]。航空数字影像虽然可以提供地表的纹理信息，但是通过数字摄影测量方法获取的数字地形模型(digital terrain model，DTM)往往不能消除植被的影响，难以获取真实的地裂缝微地貌信息[5]。机载激光扫描系统(airborne laser scanning，ALS)不受阴影及太阳高度角限制，能快速获取大面积、精细的地表三维数据。激光脉冲能够穿透植被，在植被覆盖区获取的高精度DTM，明显优于航摄图像和卫星图像。同时，LiDAR数据能够最大程度地减小地形切割带来的阴影影响[6]，通过其衍生的山体阴影图能够得到详细的地表参数，这使得在植被覆盖比较茂密的情况下，提取地裂缝等线性特征成为可能。

本文基于机载LiDAR激光点云数据，构建能刻画微地貌特征的高精度数字高程模型(digital elevation model,DEM)，设定定性及定量的识别参数，挖掘其在地裂缝提取及分析方面的能力，依据所提取的参数和微特征信息分析地裂缝的稳定性。

1 研究区概况及数据源

研究区位于湖南冷水江市，属于湖南省地质灾害比较典型的区域，地面塌陷、地裂缝现象比较常见。浪石滩位于冷水江市西南2 km处，是与新化县接壤的资江右畔丘陵区域。该区地裂缝、地表塌陷始于1987年，由于地下采煤造成坡脚和坡体深部采空，以及相应的矿区排水疏干，诱发出现地裂缝和地表塌陷等灾害。近年来，暴雨、洪水等因素加快了地表的变形发展，其地裂陷变形机制以沉陷拉裂为主。准确调查和评估地裂缝的分布及其发展状态是避让和防止次生灾害的前提，但由于该区域植被茂密，利用地面调查及其他光学影像解译等方法难以全面地获取区域地裂缝与地表塌陷的空间分布状况。

本研究利用的LiDAR数据由ALS50-II系统(详细参数见表1)所获取，点云密度6.8点/ m2，总数为6 039 290。该系统配备RCD中幅面像机，同时获取了测区的影像数据(空间分辨率10 cm)。

表1 获取实验数据采用的机载激光扫描参数

2 数据处理

2.1 激光点云滤波

激光点云滤波的目的是剔除地物点数据，得到裸露的地面点，以便于DEM的建立。机载LiDAR数据的滤波算法可以分为2大类： ①基于高程突变的滤波算法[7],该类算法的应用最为广泛，常用的有不规则三角网滤波算法、数学形态法和基于坡度的滤波方法，倘若地形高差起伏较大(如陡崖、地裂缝等)，该类算法容易将坡面点分离为地物点； ②基于激光脚点回波强度信息的滤波算法[8],可以在指定回波强度阈值范围内剔除或者分离出感兴趣的激光数据,但是只是阈值范围内激光离散点云数据的分离，需结合其他算法进一步的分类，以得到不同目标对象。

本文的目的是提取植被覆盖茂密区域地裂缝的信息，所以采用基于高程突变和激光回波强度信息相结合的滤波算法来保留地裂缝的坡面点，且分离低矮的植被点加入到地面点集中，以保证植被茂密地区地表形态的完整性。具体实现过程如下：

1)基于不规则三角网的滤波。利用不规则三角网滤波算法提取地面点，首先获取一些认为是地面点的低点构建一个初始的稀疏不规则三角格网，设定3个阈值：①最大地形坡度(terrain angle)，表达地形的起伏情况； ②最大内插距离值(interation distance)，定义内插点到三角网的距离； ③夹角(interation angle)，表达内插点和其最邻近地面点连线与三角网之间的夹角。对非地面点进行判断，每次将满足最大内插角度，到三角面的距离小于给定阈值的点纳入三角网中，同时高程阈值的选择要随着迭代次数的增加而适当减小，重复多次直到不再有新点加入为止。

2)高程滤波。在不规则三角网滤波得到非地面点的基础上，通过设定高程阈值，获取一定低高度的低矮植被点作为地面点来进行DEM的构建。将数字地面模型(digital surface model,DSM)离散点云数据减去地面点数据可得到nDSM(normalized DSM)。在nDSM中，地表地物(如建筑物、车、树等)视为在同一水平面上，可以通过选取一定的高度阈值从nDSM中提取所需的低矮植被点。

3)利用激光点云回波强度信息分离坡面点。分析由不规则三角网及高程滤波得到的非地面点的回波强度信息，设定合理的回波强度阈值，提取陡峭坡度的激光脚点。

2.2 影像的正射镶嵌

航空数字影像经过正射纠正得到的正射镶嵌影像，已经应用到变形分析、土壤湿度分析以及裂隙构造分析等领域[9]。实验选取了17张影像进行正射镶嵌，初始外方位元素由POS系统获取，差分GPS获得像机航摄中心的瞬间位置(X，Y，Z)，惯导系统获得像片姿态(φ，ω，κ)。正射镶嵌的流程为： 首先，利用多项式进行像机畸变校正[10]； 然后，加入控制点信息，空中三角测量获取加密点，光束法区域网平差解算每张影像的外方位元素； 再由空三测量计算得到的地面离散点生成DEM，进行单片正射纠正； 采用OrthoVista软件进行影像的自动匀色，调节匹配相邻影像的颜色和亮度，计算得到每张片子的辐射校正参数[11]； 最后，利用OrthoVista软件的缝合线和自动羽化功能，得到无缝的正射镶嵌图。

3 基于LiDAR的地裂缝识别参数获取

通过滤波算法分离出的地面点数据，可构建不规则三角网(triangulated irregular network,TIN)，反距离加权内插可以得到DEM。本文以上述方法构建的TIN和DEM数据为基础，设定定性及定量识别参数，并挖掘这些参数对地裂缝信息提取及分析的能力。

3.1 山体阴影图

山体阴影图利用假想光源对地表进行照射，产生地形表面的阴影图，对实际地形特征进行“逼真”的立体模拟，增强地面的起伏感。它能够凸显很多微地貌特征，尤其是一些线性特征，比如地裂缝、线性构造、断层崖等。本文采用山体阴影图作为地裂缝识别的定性参数。

山体阴影图的计算有3个重要参数： ①太阳方位角，即假想光源发出的光线在水平面上的投影与正北方向的夹角，确定其值需分析试验区地裂缝的主体走向及河流山脉等的走向； ②太阳高度角，即假想光源发出的光线与水平面的夹角，一般选45°； ③地表灰度值，取值范围为0～255，或者依据高程变化情况选择最为合适的值。

3.2 坡度坡向图

地面上某点的坡度表示了地表在该点的倾斜程度。根据坡度起伏变化可以确定地裂缝、泥石流或严重的土壤侵蚀区。坡向(aspect)定义为坡面法线在水平面上的投影与正北方向的夹角，可确定地裂缝的走向。本文将坡度(slope)和坡向图作为地裂缝识别的定量参数。

地面上某点的坡度(S)和坡向(A)为地形曲率在东西(Y)、南北(X)方向上高程变化率的函数，即

,

(1)

,

(2)

式中：fx为X方向高程变化率；fy为Y方向高程变化率。通常，fx和fy的求解，是在3×3的移动窗口中(图1)，通过数值微分或局部曲面拟合方法进行。

图1 差分DEM示意图

坡度运用三阶反距离平方权差分算法计算，即

(3)

式中：d为格网分辨率；Zi(i=1,2,…，9)为中心点“5”周围i格网点的高程。

3.3 线性提取

目前，利用LiDAR点云数据自动或半自动提取地表断裂线的研究比较少。Brügelmann利用LiDAR点云高程值重采样为距离图像，以图像处理的方法计算像素的二阶导数来提取断裂点的候选点，并进一步拟合成断裂线[12]，这种方法不适应于有高差起伏的山区。Briese通过局部平面拟合相交的方法来提取断裂线[13]，该方法直接对离散LiDAR点云进行处理，较好地保持了断裂线的精度，但是需要人工给出断裂线的起始位置和方向，并且进行面片拟合相交逐步提取，计算量较大。本文采用基于曲率等值线的方法提取地裂缝矢量线，利用离散的三维地面点云数据构建TIN格网，追踪剖面最小曲率等值线，该方法即能较好地保持地裂缝的精度，又可减少相对于离散点云的计算量。

地形表面曲率是局部地形曲面在各个界面方向上的形状、凹凸变化的反映，显示出地形曲面在不同方向上的结构和形态特征。Wood提出用纵向曲率(longitudinal curvature)和断面曲率(cross section curvature)来进行地形特征的识别和提取[14]。地形剖面曲率的实质为地面坡度的变化率，即

(4)

4 实验结果与分析

4.1 正射镶嵌及点云滤波精度

研究区位于资江右畔的丘陵区域，由正射镶嵌后的影像(图2)可以看出，研究区植被覆盖茂密，建筑物、农田较多。在研究区布设了6个控制点来评定正射影像的精度，X，Y，Z这3个方向的均方根误差RMSE分别为0.282 m，0.259 m和0.883 m。

图2 正射镶嵌影像

首先，对激光点云数据进行滤波处理，设定3个阈值参数分别为：Terrainangle=88°，Interationdistance=2 m ，Interationangle= 8°，进行不规则三角网滤波算法获取初始地面点； 然后，设定0.6 m的高程阈值参数，提取低矮植被加入到地面点集中； 最后，分离灰度反射强度为22～30、高程为0.6～1.2 m的植被点到地面点集，以保留陡峭的坡面点。采用目前应用最广泛的商业软件Terrasoild提取的地面点仅占点云总数的13.1%(表2)，且有地裂缝的坡面点被误分为植被点的现象(图3)。

表2 不同方式获取的地面点对比

(a) TerraSoild软件提取的地面点和植被点(b) 激光点云的空间分布

为了更好地识别与提取地裂缝，本文在基于高程突变滤波的基础上，分析裸地反射强度信息，保留了地裂缝的坡面点(图4)，分离的地面点数目占点云总数的43.2%(表2)。

(a) 本文方法提取的地面点和植被点(b) 激光点云的空间分布

利用滤波分离出的地面点构建DEM，为后续地裂缝的提取与识别提供数据基础。原始激光点云数据内插生成的DSM见图5(a)。由于地表地物(植被、房屋等)的影响，很难辨识地裂缝的位置。在植被覆盖的区域，利用数字摄影测量手段提取的地面点较少，构建的DEM(图5(b))比较粗糙，难以刻画地表真实形态，尤其是地裂缝等地表微特征。利用TerraSoild软件滤波分离地面点时，地裂缝的坡面点被误分为植被点，因此生成的DEM(图5(c))虽然在视觉效果上更能突显地裂缝的位置，但是却影响后续地裂缝提取的精度。图5(d)为采用本文提出的滤波方法分离地面点派生的DEM，能够更细致地刻画地表形态，有效地突显地裂缝等地表微特征。

(a) 原始激光点云派生出的DSM(b) 传统数字摄影测量方法提取的DEM

(c) Terrasoild软件滤波提取地面点所派生出的DEM(d) 本文方法提取地面点所派生的DEM

4.2 识别参数及地裂缝微特征分析

基于生成的DEM，设定太阳方位角315°及太阳高度角45°，得到的定性参数山体阴影图(图6(a))，其表达地表形态的立体能力比较强。对于植被覆盖茂密的区域，较正射镶嵌影像(图6(b))而言，山体阴影图(图6(c))更能清晰地反映地裂缝的起伏及破碎特征。

(a) 实验区山体阴影图

(b) 局部有分布地裂缝处的正射影像(c) 局部有分布地裂缝处的山体阴影图((a)图黑色方框范围内)((a)图黑色方框范围内)

图7为研究区的坡度和坡向图。地裂缝分布的区域为坡度比较大(约在30°～50°之间)的地区，分布规律性较强，地裂缝处的坡度最小为40°，最大为65°，大部分位于50°～60°之间。坡向表示斜坡体表面的任一高程值的变化量的最大变化方向，由坡向图可以看出，地裂缝的坡向为NW方向，位于南坡区域，即阳坡，而坡向决定斜坡的稳定性，阳坡的侵蚀强度往往远高于阴坡。

对于地裂缝线性提取，首先用获取的地面点数据构建TIN，追踪曲率的最小值Kmin，得到地裂缝矢量线(图8)，长度分别为(从左到右)24.8 m，43.2 m和81.8 m，总长149.8 m。将矢量线与正射影像图相叠合，如图8(b)所示，可以看出地裂缝提取的位置比较准确。

(a) 地裂缝矢量线与TIN格网叠合图(b) 地裂缝矢量线与DOM叠合图

图8地裂缝矢量线与TIN格网及DOM叠合图

Fig.8CombinationofgroundfissurecontourandTIN(left)andDOM(right)

提取地裂缝的准确位置之后，可通过点云的剖面信息浏览，来分析地裂缝的微细节，如深度、宽度及几何形态等(图9)。

(a) 地裂缝A处影像(左)及其点云剖面(右)(b) 地裂缝B处影像(左)及其点云剖面(右)

(c) 地裂缝C处影像(左)及其点云剖面(右)(d) 地裂缝D处影像(左)及其点云剖面(右)

地裂缝的形态为规则的V字形(图9(a)(c)右)； 由图9(b)(c)可以发现，这2处地裂缝的坡面有变缓的现象，而根据地裂缝提取辅助识别参数判别，2处坡面都为坡度陡峭、雨水侵蚀较为严重的阳坡，坡面的土壤有塌陷，地裂缝有继续发育的可能。

5 结论

本文以机载LiDAR(ALS)数据为基础，开展激光点云数据对地裂缝识别提取能力的研究。主要结论如下：

1)本文依次基于不规则三角网滤波、高程滤波及回波信息强度滤波提取地面点，很好地保留了地裂缝陡峭边的坡面点，保证了地表微特征的完整性。

2)提取了地裂缝辅助识别参数，定性参数山体阴影图能突显地裂缝的地貌特征，以确定地裂缝分布区域，由定量识别参数坡度坡向图可确定地裂缝发生在坡度为30°～50°的区域，坡向为NW方向，分布规律性较强。

3)为了提取地裂缝的矢量线，得到地裂缝的长度及准确位置信息，本文追踪TIN剖面最小曲率等值线，获取地裂缝矢量线，再与正射镶嵌影像相叠合，定性分析确定其位置，结果比较精确。并通过对地裂缝剖面信息与辅助识别参数的综合分析，得出了研究区地裂缝有继续发育可能的结论。

参考文献(References)：

[1] Jachens R C,Holzer T L.Differential compaction mechanism for earth fissures near Casa Grande,Arizona[J].Geological Society of America Bulletin,1982,93(10):998-1012.

[2] Holzer T L,Davis S N,Lofgren B E.Faulting caused by groundwater extraction in south-central Arizona[J].Journal of Geological Research,1979,84(B2):603-612.

[3] Belardinelli M E,Sandri L,Baldi P.The major event of the 1997 Umbria-Marche(Italy)sequence:What could we learn from DInSAR and GPS data[J].Geophysical Journal International,2003,153(1):242-252.

[4] Niebergall S,Loew A,Mauser W.Object-orientated analysis of very high resolution QuickBird data for mega city research in Delphi/India[C]//Proceedings of the Urban Remote Sensing Joint Event,Paris.ISBN:1-4244-0712-5.2007:8-13,IEEE07EX1577.

[5] Baltsavias E P.A comparison between photogrammetry and laser scanning[J].ISPRS J Photogrammetry and Remote Sensing,1999,54(2):83-94.

[6] Carter W E,Shrestha R L,Slatton K C.Geodetic laser scanning[J].Phys Today,2007:60(12):41-47.

[7] 张小红.机载激光雷达测量技术理论与方法[M].武汉:武汉大学出版社,2007:93-116。

Zhang X H.The Theory and Technical Method of Airborne Laser Radar Measurement[M].Wuhan:Wuhan University Process,2007:93-116.

[8] 路兴昌,张雪霞.基于回波强度和采样点距离的点云滤波研究[J].测绘科学,2009,34(6):196-197.

Lu X C,Zhang X X.Study on points cloud filtering based on reflectance intensity and range[J].Science of Surveying and Mapping,2009,34(6):196-197.

[9] Niethammer U,Rothmund S,Joswig M.UAV-based remote sensing of the slow moving landslide Super-Sauze[C]//Malet J P,Rema tre A,Boogard T.Proceedings of the International Conference on Landslide Processes:From Geomorphologic Mapping to Dynamic Modelling.Strasbourg:CERG Press,2009:69-74.

[10]Niethammer U,Rothmund S,James M R,et al.UAV-based remote sensing of landslides[C]//International Archives of Photogrammetry,Remote Sensing and Spatial Information Sciences,Vol.XXXVIII,Part 5 Commission V Symposium,Newcastle upon Tyne,UK.2010:496-501.

[11]OrthoVista.Official OrthoVista software homepage[EB/OL].http://www.orthovista.com.(accessed 1 August 2010).

[12]Brügelmann R.Automatic breakline detection from airborne laser range data[J].International Archives of Photogrammetry and Remote Sensing,2000,33(B3):109-115.

[13]Briese C.Three-dimensional modelling of breaklines from airborne laser scanner data[J].International Archives of Photogrammetry Remote Sensing,2004,35(B3):109-I102.

[14]Wood J.The geomorphological characterization of digital elevation models[D].Leicester,UK:University of Leicester,1996.