利用地形开度算子自动提取激光雷达点云断裂线

张 良，张 帆，姜晓琦，周佳雯，王杰栋

(1. 湖北大学资源环境学院，湖北 武汉 430062； 2. 区域开发与环境响应湖北省重点实验室，湖北 武汉 430062； 3. 武汉大学遥感信息工程学院，湖北 武汉 430079； 4. 浙江省第二测绘院，浙江 杭州 310012)

LiDAR系统是一种集激光、全球定位系统(GPS)和惯性导航系统(INS)3种技术于一身的新型传感器，无需大量地面控制点，便能够快速准确地获取地表高密度、高精度的三维点云，同时具备全天时观测能力、受天气影响小等优点，是高精度DEM的优质数据源之一[1-2]。地形表面通常不都是光滑的或均匀变化的，在陡坎、梯田、冲沟、堤岸等处，常出现转折或突变，被称为地形断裂线[3]。断裂线是建立DEM非常重要的依据，对数字高程模型的精度有相当大的影响，尤其在地形变化大的山区或水系多的地区[4]。因此，为了有效保留地表的关键地形特征，为用户提供精度更高的DEM，需要直接从LiDAR点云数据中提取地形断裂线作为构建DEM的约束条件。

近10年来，已有学者对基于LiDAR点云的断裂线提取进行了研究并取得了一定的成果。现有的断裂线自动提取算法包括断裂线候选点提取和断裂线矢量跟踪两大步骤，其中断裂线候选点的精度是断裂线提取精度和完整性的前提和关键，也是本文重点探讨的地方。文献[5—6]利用二阶导数等图像边缘算子处理LiDAR距离图像获取断裂线候选点并拟合断裂线。这类算法自动化程度高，但是自动提取的断裂线比较破碎且错提取率高。文献[7]通过局部平面相交的方法直接从离散点云提取较精确的断裂线候选点，该方法一般需要人工给出断裂线的起始位置和方向作为断裂线生长的初始值。文献[8]尝试利用三角网面片之间的法向差异提取候选断裂点，再基于方向优先策略追踪断裂线以提高算法的自动化程度。文献[9]首先利用双阈值从距离图像提取断裂线概略位置，再结合原始点云数据用分段构造法获取更精确的断裂线。双阈值法提取的断裂线完整性好，但大、小阈值的设置需要多次试验才能有效确定。

如前文所述，目前大部分断裂线提取算法都是基于地表局部邻域的变化开展，因此对地表突变和点云错分敏感，精度和完整性无法满足生产需求，因此目前有经验的处理人员更信任结合可视化、光照分析等工具判断地形的整体变化，进而人工绘制断裂线。目前大部分商业软件(如Qtmodeler)利用阴影(或与高程结合)渲染的方式表现三维地形细节。阴影渲染下，地表亮度与光源位置、地形起伏密切相关，地形起伏大的区域，部分地表细节因无法接受足够光照而被遗漏，而部分区域则可能曝光过度[10-11]。因此很难基于阴影渲染图提取完整的断裂线。地形开度(topoghaphic openness)算子利用均质虚拟辐射源代替单方向点阵直射光源，从根本上突破了基于阴影渲染采用平行直射光源导致的阴影和曝光过度问题。近几年来，作为一种新型的高分辨率地形可视化工具，在军事、考古、地质[12-15]等需要对微地貌进行详细解译的领域得到广泛应用。本文从人工解译的角度出发，将地形开度作为一种定量描述地表整体变化的地形特征，提出基于地形开度(topographic openness)的断裂线自动提取算法，以提高复杂地区高精度DEM的提取能力。

1 算法设计

1.1 地形开度算子

地形开度算子或SVF(sky view factor)算子是近几年发展的利用均质辐射源代替单方向点阵直射光源的地形渲染算法。SVF算法从能量传播的角度出发，假设光源为一个虚拟的半圆形天球，待渲染的地表位于该天球的中心并受天球中均匀的散射光所照射，在该假设下，可以认为地表上某个点的亮度取决于该点及其邻域暴露在散射光源下的体积或表面积(如图1(b)所示)。图1(a)表示起伏地表的一个剖面截图，其中A为谷底待计算亮度值的像素，扇形区域代表A能够接收均质散射光辐射区域。因此，通过式(1)可估算出该点暴露在天球光线照射的体积Ω即该点的亮度。

(1)

式中，φ表示均值天球球体的纬度；λ表示经度。

图1 基于辐射度的地形渲染原理

地形开度[10]则相当于SVF的简化算子。地形开度包括正开度和负开度两个部分，正开度用于表示山脊、背脊等隆起地貌，负开度表示地形的整体凹陷程度，着重表达山谷、沟壑等地形。正开度通常取多个方位的天顶角的平均值(一般取8个方向)。图2表示方位角为90°和270°的天顶角剖面示意图。首先在一定半径内寻找点B，使得天顶角Ωi最小，即点A和点B之间的仰角γi最大。利用式(2)计算γi，再基于式(3)即可计算天顶角Ωi。通过同样的原理计算其余方向天顶角值之后，可根据式(4)计算出该点的正开度值Op。

(2)

Ωi=90-γi

(3)

(4)

式中，γi为点A和点B之间的仰角；HA、HB分别为点A和点B的高程；dAB为点A和点B之间的水平距离；Op为正开度值；n表示所取方向的数量，一般n值取8。

图2 天顶角计算示意图

负开度的计算原理和正开度类似，通常取8个方位的天底角的平均值。首先在搜索区域内找到点C，使得天底角Ψi最小，即点A和点C之间的俯视角δi最大，然后利用式(5)、式(6)计算该方向的天底角Ψi(图3为方位角90°和270°下天底角的示意图)。计算其余方向的天底角之后，利用式(7)统计负开度值On。

(5)

ΨiL=90+δi

(6)

(7)

式中，δi为点A和点B之间的俯视角；HA、HC分别为点A和点C的高程；dAC为点A和点C之间的水平距离；On为负开度值；n为所取方向的数量，一般n值取8。

图3 天底角计算示意图

1.2 基于地形开度算子的断裂线自动提取

大量的文献指出[12-14]，地形开度算子具有两大特点：①开度算子表现地表的主要变化，即地形结构主要变化得到加强，而缓慢的连续形变或小范围的突变则得到抑制，因此地形开度算子渲染后的地表具有极强的立体感和纵深感，也使得地形特征和周边地形的区别更加明显。②开度算子计算的地表亮度仅仅与搜索半径大小有关，不会因为光照条件的不同和地形起伏而遗漏地表细节。基于上述特点，本文在获取开度算子的基础上，实现断裂线的自动提取。算法具体步骤如下：

步骤1：利用不规则三角网(Tin)加密算法实现点云滤波并提取DEM。

步骤2：基于DEM分别提取正负开度Op和On，并且通过正负开度相减得到IDEM(式(8))。如图4(a)所示，若点在地表凹陷处，正开度小于负开度，IDEM取负值。如图4(b)所示，点落在起伏较为平缓的表面，其正、负开度值相近，IDEM值接近零。如图4(c)所示，若点落在断裂附近或山脊等突起地表附近，正开度大于负开度，IDEM取正值。因此，根据IDEM可区分凸出地表、较平缓地表(含倾斜地表)和凹陷地表。

I=(Op-On)/2

(8)

图4 正负开度和地表的关系

步骤3：通过边缘检测LOG算子对IDEM进行边缘提取。对该边缘提取结果进行二值化，利用面积阈值删除连通区域小的边缘像素，得到的高亮区域即为断裂线种子点及少量噪声。

步骤4：利用边缘检测LOG算子对原始DEM进行边缘提取。可发现地形断裂区域，地表起伏大的山坡、小的突起处都呈高亮状态，但是地形平坦区域的亮度明显低于其他区域，具有非常明显的区分度。基于该特点，通过亮度分割的方式提取地形平缓区域，对其进行亮度反转，然后利用形态学开运算得到较完整的平缓区域。

步骤5：通过叠加求并运算，利用步骤4提取的地形平坦区域进一步消除步骤3提取的噪声，提取结果即为断裂线种子点。

步骤6：利用文献[8]的方法，结合原始点云采用基于平面对的分段构造法和基于三角网的直接构造法提取完整的断裂线矢量。

2 试验和精度评价

2.1 试验数据和试验工具

为了检验算法的有效性，笔者选取一块中国陕西某地区的点云数据(如图5(a)所示)进行试验。试验数据位于黄土高原，由于植被覆盖较少，土质疏松，气候干旱，导致水土流失严重，侵蚀切割强烈，断裂和裂隙纵横交错。试验区域范围为0.6 km×0.6 km，平均点间距为0.6 m，激光脚点数量为1 418 228，试验区内包含大量的不同种类的断裂、村庄和低矮植被。利用LiDAR点云数据处理平台LiDAR_Suite对试验数据进行滤波和人工分类并基于LiDAR_Suite软件进行二次开发，利用C++语言实现了本文算法。

2.2 试验过程和结果分析

(1) 在LiDAR点云处理平台LiDAR_Suite上实现点云精细分类，然后基于线性三角网插值算法获取1 m分辨率的DEM，如图5(b)所示。

(2) 基于DEM提取开度图像IDEM。根据文献[11]的建议，开度半径设置为50像素。如图5(c)所示，开度图像中，突起地形(如断裂线或道路桥梁边缘)呈明显的高亮状态，山坡、山谷则亮度值非常低，具有非常明显的区分度。

(3) 利用LOG算子对IDEM进行边缘提取，其中高斯拉普拉斯卷积核的大小为25×25像素，高斯方差设为2.5，并且利用灰度直方图自动进行二值分割，提取结果如图5(d)所示，二值化后的结果基本上涵盖了全部地形断裂，以及少量由于地形突变和错分类导致的噪声。

(4) 利用LOG算子直接对DEM进行边缘提取，LOG算子参数与步骤3相同，处理结果如图5(e)所示，可见在黄土高原这种特殊地形条件下，直接基于LOG算子也能取得较好的效果，但是断裂线完整度不如基于IDEM的提取结果，噪声数量也明显更大。另一方面，如图5(e)所示，利用LOG算子处理后的DEM，地形平坦区域的亮度明显低于其他区域，利用该特征对边缘提取结果进行二值化处理，并且通过形态学算子提取较完整的平缓区域，结果如图5(f)所示。

(5) 基于步骤3和步骤4的提取结果进行叠加分析，消除噪声，得到完整的断裂线种子点。在此基础上，结合原始点云数据采用基于平面对的分段构造法和基于三角网的直接构造法提取断裂线矢量。

图5 基于开度的断裂线提取

(6) 从两个方面来评价本文方法自动提取断裂线的效果：①将最终提取的断裂线矢量和开度图叠加，结果如图5(g)所示，可见本文所提取的断裂线具备非常好的完整性，同时与实际的断裂边缘吻合度也较好；另一方面，地形结构变化不够明显的区域，如图5(g)左下角的桥梁边缘没有识别，主要原因在于为了更好地表达地表整体的变化，开度半径设置偏大，因此小的结构变化不够突出。②在LiDAR_Suite平台通过人机交互的方式提取断裂线，与本文算法进行对比，对比结果如图5(h)所示，其中左侧为本文算法提取的断裂线，右侧为人工提取的断裂线，通过对比发现，两种方式提取的断裂线基本一致，通过本文提取的断裂线细节更加丰富，保真度更高。

3 总结和展望

DEM的制作一直是测绘地理信息系统领域的重点之一，提取地形断裂线并作为构建DEM的约束条件，可以有效减少地貌的失真，提供更高精度的DEM。本文从人工解译的角度出发，引入地形开度作为一种定量描述地表整体变化的地形特征，提出了基于开度算子的断裂线自动提取算法，突破了传统方法仅仅基于地形局部变化考虑的局限，结合边缘提取算子和形态学算子，能够快速提取完整准确的断裂线。算法简洁有效，不需要人工干预，具备较好的推广价值。

开度算子表示的地形细节的丰富程度与搜索半径相关，为了提取更高精度的断裂线，需要研究基于不同地形地貌自适应调整开度半径的方法；同时，本文着重探讨了断裂线候选点的提取方法，如何结合开度算子，基于断裂线候选点提取更精确的断裂线矢量，还需要进一步研究。另外，目前对断裂线提取精度的评价仅仅是与影像和手工提取的结果进行比较，还难以进行定量分析，这也是今后需要深入研究的问题。

参考文献：

[1] AXELSSON P.Processing of Laser Scanner Data-Algorithms and Applications[J].Isprs Journal of Photogrammetry & Remote Sensing，1999，54(2-3):138-147.

[2] 郑辑涛，张涛.基于可变半径圆环和B样条拟合的机载LiDAR点云滤波[J].测绘学报，2015，44(12):1359-1366.

[3] YANG B，HUANG R，DONG Z，et al.Two-step Adaptive Extraction Method for Ground Points and Breaklines from Lidar Point Clouds[J].Isprs Journal of Photogrammetry & Remote Sensing，2016，119:373-389.

[4] 易俐娜，刘鹏飞，乔小军,等.结合遥感影像和DEM的线性体特征提取[J].测绘通报，2014(10):19-22.

[5] UGELMANN R B.Automatic Breakline Detection from Airborne Laser Range Data[J].International Archives of Photogrammetry and Remotes Sensing,2000,33(B3):109-115.

[6] TOSCANO G J，GOPALAM U，DEVARAJAN V.A Novel Method for Automation of 3D Hydro Break Line Generation from LIDAR Data Using Matlab[J].ISPRS-International Archives of the Photogrammetry，Remote Sensing and Spatial Information Sciences，2013，XL-2/W2(1):99-104.

[7] BRIESE C.Three-dimensional Modelling of Breaklines from Airborne Laser Scanner Data[J].International Archives of Photogrammetry Remote Sensing & Spatial Information Sciences，2011，35:1097-1102.

[8] 徐景中，万幼川，张圣望.基于机载激光雷达点云的断裂线自动提取方法[J].计算机应用，2008，28(5):1214-1216.

[9] 彭检贵，马洪超，王宗跃,等.机载LiDAR点云的双阈值自动提取断裂线方法[J].测绘科学技术学报，2010，27(4):275-279.

[10] SHIRASAWA M，YOKOYAMA R.Visualizing Topography by Openness: A New Application of Image Processing to Digital Elevation Models[J].Photogrammetric Engineering & Remote Sensing，2002，68(3):257-266.

[11] NGUYEN H T，PEARCE J M.Incorporating Shading Losses in Solar Photovoltaic Potential Assessment at the Municipal Scale[J].Solar Energy，2012，86(5):1245-1260.

[12] DONEUS M.Openness as Visualization Technique for Interpretative Mapping of Airborne Lidar Derived Digital Terrain Models[J].Remote Sensing，2013，5(12):6427-6442.

[13] ELMAHDY S I.Topographic Openness Algorithm for Characterizing Geologic Fractures of Kuala Lumpur Limestone Bedrock Using DEM[J].Journal of Geomatics，2010，4(2):61-68.

[14] HODUL M，KNUDBY A，HO H.Estimation of Continuous Urban Sky View Factor from Landsat Data Using Shadow Detection[J].Remote Sensing，2016，8(7):568.

[15] FAVALLI M，FORNACIAI A.Visualization and Comparison of DEM-derived Parameters Application to Volcanic Areas[J].Geomorphology，2017，290: 69-84.