，

(武漢大學 遙感信息工程學院，武漢 430079)

1 研究背景

傳統(tǒng)的等高線提取主要基于攝影測量和地面測量。提取原理是首先通過空間后方交會或從星歷中解算出影像的外方位元素，然后進行影像匹配得到同名點的像坐標，最后通過前方交會解算出地面點空間三維坐標后內插得到格網(wǎng)DEM，根據(jù)DEM生成等高線。傳統(tǒng)方法效率低且成本高，由于需要大量人工操作，受人的經(jīng)驗因素影響較大。

隨著LiDAR(Light Detection and Ranging)技術的發(fā)展，利用LiDAR數(shù)據(jù)自動提取等高線的方法應運而生，彌補了傳統(tǒng)方法的不足。這類方法主要基于高程格網(wǎng)或TIN的等高線跟蹤方法。如Batcha等[1]提出離散點內插成規(guī)則格網(wǎng)或根據(jù)離散點建立三角網(wǎng)，再從規(guī)則格網(wǎng)或三角網(wǎng)上跟蹤提取等高線的方法；王濤等[2]提出一種優(yōu)化高程格網(wǎng)內插等高線的方法。由于LiDAR點云密度大，為避免內插的方法導致精度損失，保證地面精度，一般采用基于TIN的等高線跟蹤方法，如Boissonnat[3]的加密算法。根據(jù)國家成圖規(guī)范，等高線不能穿過建筑物、水系地物等所在的約束區(qū)域。傳統(tǒng)基于TIN的等高線跟蹤方法一般采用人工編輯處理約束區(qū)域的等高線，存在人工勞動強度過大的缺點。本文提出一種基于特征約束的等高線提取方法，先進行特征提取，然后自動生成特征約束的TIN三角網(wǎng)，最后在特征約束TIN的基礎上提取等高線。

2 LiDAR點云濾波

LiDAR點云數(shù)據(jù)是高精度的地表三維點云數(shù)據(jù)(DSM)，它除了包含地面點之外, 還包含地表的地物點，即非地面點。為了獲取反映真實地表的DEM，提高等高線精度，需要對DSM點云數(shù)據(jù)進行濾波，從而去除非地面點，得到裸露的地面點?，F(xiàn)有濾波的方法主要是基于激光點的幾何特性對地面點和非地面點進行區(qū)分，主要有以下3類：①基于形態(tài)濾波的方法[4-5]。這類方法的優(yōu)點是直觀，缺點是當?shù)匦巫兓瘎×視r，處理效果不夠理想；②基于曲面約束的方法。這類方法適用于平坦地形，對于地形起伏劇烈的山區(qū)效果不佳，且計算量大，運行效率較低；③基于內插的方法[6]。這類方法原理清晰，計算簡單，能較好處理地形起伏較大的山區(qū)，但要求預先剔除粗差。

本文數(shù)據(jù)采自長陽山區(qū)，地形較為復雜、建筑物較多，為滿足TIN算法對DEM精度的要求，本文采用改進的漸進加密三角網(wǎng)濾波算法[7]。為驗證改進的漸進加密三角網(wǎng)濾波的有效性，將試驗區(qū)點云數(shù)據(jù)分別利用形態(tài)濾波方法和基于曲面約束的濾波方法分別進行處理，并與改進的漸進加密三角網(wǎng)濾波法的濾波結果對比。通過比較幾種濾波方法產(chǎn)生的第1類錯誤即將地面點判為非地面點和第2類錯誤即將非地面點判為地面點，對濾波效果進行定量評價。為了更好地分析每種方法的濾波效果，本文分別對平地、山區(qū)、建筑物區(qū)和植被區(qū)這4種典型地形區(qū)域進行濾波處理和分析，其結果如表1。本文試驗所用LiDAR數(shù)據(jù)采自長陽山區(qū)，有平地、山地及建筑物，且植被覆蓋率較高、地形起伏較大、地物構成比較復雜。由表1可知，采用穩(wěn)健性較強的改進漸進加密三角網(wǎng)濾波法效果最優(yōu)，采用該方法得到的DEM三維渲染圖如圖1所示。

表1 濾波錯誤率

圖1 DEM三維渲染圖

3 特征提取

特征約束區(qū)域提取是生成特征約束TIN三角網(wǎng)的基礎，通常通過提取典型地物的邊緣信息可得到典型地物的特征約束區(qū)域。在基礎測繪生產(chǎn)應用中，建筑物和水系這2類地物的邊緣特征最為典型，本文重點介紹建筑物及水系的邊緣提取方法。

3.1 建筑物

目前基于LiDAR數(shù)據(jù)進行建筑物提取的方法有很多,這些方法主要可以分為2類：第1類是基于原始點云再處理的建筑物激光點提取法(如崔建軍[8]，任自珍[9])，這類方法的缺點是在再處理的過程中無法避免地引入了不必要的誤差。第2類方法是基于濾波分類的建筑物激光點提取法。這類方法的代表是區(qū)域增長法(如Nicholas等[10])。此類方法以真實的數(shù)據(jù)為基礎，邊緣提取精度和可靠性較高，因此本文采用第2類方法。

(1) 在TIN形式的DEM中，建筑物區(qū)域常對應邊長較大的三角形。利用建筑物區(qū)域三角形的邊長較大的特性，設定適當?shù)拈撝?，可提取出建筑物區(qū)域，得到建筑物點云。

(2) 建筑物區(qū)域邊緣跟蹤。α-shapes[11]算法可以高效地提取各類建筑物邊緣。

3.2 水 系

水系所在區(qū)域的點云數(shù)據(jù)具有點密稀疏、高程相近、回波強度弱和水面低于周圍陸地的特點?；谒翟邳c云數(shù)據(jù)上的以上4個特征，采用直接利用LiDAR點云數(shù)據(jù)提取水系輪廓線的方法[12]。該方法首先利用單層粗格網(wǎng)提取水域。將滿足水系點云4個特征條件的格網(wǎng)塊判斷為水系的格網(wǎng)塊;再利用雙層格網(wǎng)提取較窄的水域。雙層格網(wǎng)的處理可以使得較窄的河流不易斷裂,在連通去噪后還能保留下來。單層及雙層格網(wǎng)的示意圖如圖2所示。

圖2 單層及雙層格網(wǎng)示意圖

在水體格網(wǎng)提取的基礎上,通過拉普拉斯邊緣提取算子獲得邊緣格網(wǎng)，求出各個邊緣格網(wǎng)中面向水系的邊緣點，試驗區(qū)數(shù)據(jù)覆蓋的主要水系為長陽的清江，其邊緣提取結果見本文4節(jié)中的圖3中深藍色曲線所示。

圖3 建筑物區(qū)域的約束TIN

4 基于特征約束的三角網(wǎng)TIN生成

本文第3節(jié)提取的建筑物及水系邊緣線所圍成的區(qū)域即為特征約束區(qū)域，將在約束區(qū)域內生成的三角網(wǎng)自動刪除，就能得到約束三角網(wǎng)。常用的建立約束三角網(wǎng)的方法很多，如：基于約束圖的CDT構造算法[13],非約束數(shù)據(jù)域Delaunay三角剖分的分割合并算法[14],加密算法[3],2步法[15-16]等。2步法處理方式更靈活且效率更高，因此本文采用2步法建立約束三角網(wǎng)。其結果如圖3所示。右上角圖為紅色框內區(qū)域的放大顯示。紫色線為建筑物邊緣線，藍色線為水系邊緣線，建筑物和水系邊緣線圍成區(qū)域內的三角網(wǎng)被剔除。

5 基于特征約束TIN的等高線生成

依據(jù)基礎測繪線畫圖產(chǎn)品規(guī)范，等高線不能穿過建筑物和水系。且在第4節(jié)中，已將建筑物及水系區(qū)域的三角網(wǎng)自動刪除，因此，本節(jié)將研究在帶有特征約束的TIN的基礎上，實現(xiàn)等高線的自動生成，其主要流程如下：

(1) 基于約束三角網(wǎng)跟蹤等高點。依次判斷每條等高線是否與三角形的邊相交，將與等高線相交的三角形內插得到等高點坐標，內插公式為

(1)

(2)等高線平滑。將所有等高點連接起來得到的是不光滑的折線，需使用光滑算法進行處理以得到光滑等高線。本文采用精度較高的張力樣條函數(shù)插補法[17]，主要原理：設(x1,y1),(x2,y2),…(xn,yn)是一組已知的數(shù)據(jù)點，且x1

f″-σ2f(x)=[f″(x)-σ2yi](xi+1-x)/Hi+

[f″(xi+1)-σ2yi+1](x-xi)/Hi。

(2)

張力樣條數(shù)為單值情況下公式(2)的解要得到更精確和美觀的曲線，需要設定合適的張力系數(shù)值。整個試驗區(qū)等高線如圖4所示，建筑物區(qū)域等高線局部放大圖右上角所示，可見等高線未穿過建筑物和水系。

圖4 建筑物區(qū)域等高線

6 等高線精度評定

等高線精度檢查利用導線測量得到的檢查點的高程值與等高線內插出的對應檢查點的高程值進行比較，從而得到等高線的高程精度。實驗選取52個檢查點，以等高線高程Z為橫軸，誤差DeltaZ為縱軸繪制誤差分布散點圖，等高線誤差分布圖如圖7所示。其中平均高差為+0.022 m，最小高差為-1.221 m，最大高差為+1.570 m，平均誤差為0.145 m，最小二乘殘差為0.326 m，方差為0.328 m，滿足1∶2 000 DLG成圖規(guī)范的要求。

圖5 等高線誤差分布

7 結 語

本文以特征約束的TIN三角網(wǎng)為基礎，利用LiDAR點云數(shù)據(jù)實現(xiàn)長陽山區(qū)等高線自動提取。

針對長陽測區(qū)數(shù)據(jù)地形起伏較大，地物構成復雜的特點，采用穩(wěn)健性較高的改進的漸進三角網(wǎng)濾波方法生成高精度的DEM點云數(shù)據(jù)，實現(xiàn)地面高程模型TIN三角網(wǎng)的自動生成，并以每個三角網(wǎng)邊的長度為約束條件，實現(xiàn)建筑物和水系這2類典型地物邊緣的自動提取，同時自動剔除特征區(qū)域內的三角網(wǎng)，生成帶有特征約束的TIN三角網(wǎng)。在此基礎上，判定等高線是否穿越當前三角網(wǎng)，利用三角網(wǎng)頂點高程值內插，自動提取等高線。本文在試驗區(qū)域選取52個檢查點，對生成的等高線進行了精度評價，評價結果表明該方法流程滿足1∶2 000 DLG成圖規(guī)范的要求，與傳統(tǒng)的等高線提取方法相比，該方法更具科學性、可靠性，節(jié)省了大量人工操作，具有效率快和作業(yè)精度高的特點。

但針對以下2種情況：①地形起伏變化劇烈，如懸崖、陡坎等比較豐富的區(qū)域；②經(jīng)過人工改造的平坦區(qū)域，如農田、果園、梯田等區(qū)域，特征提取需要進一步完善，TIN三角網(wǎng)的生成條件會更加嚴格，因此等高線的生成算法，有待于進一步研究和完善。

參考文獻：

[1] BATCHA J, REESE J. Surface Determination and Automatic Contouring for Mineral Exploration, Extraction, and Processing[J]. Colorado School of Mines Quarterly, 1964, 59: 1-4.

[2] 王 濤,毋河海.一種從高程格網(wǎng)中提取等高線的算法[J]. 測繪科學,2006, 31(2):108-110. (WANG Tao, WU He-hai. An Algorithm for Extracting Contour from Elevation Grid[J]. Science of Surveying and Mapping, 2006, 31(2): 108-110. (in Chinese))

[3] BOISSONNAT J. Shape Reconstruction from Planar Cross Sections[J]. Computer Vision, Graphics, and Image Processing, 1988 ,(1): 1-29.

[4] KILIAN J, HAALA N, ENGLICH M. Capture and Evaluation of Airborne Laser Scanner Data[J] . IAPRS,1996,31(B3) : 383-388.

[5] VOSSELMAN G. Slope Based Filtering of Laser Altimetry Data[J]. IAPRS, 2000，(B3) : 935-942.

[6] AXELSSO N P. DEM Generation from Laser Scanner Data Using Adaptive TIN Models [J]. IAPRS,2000,33(B4/1) : 110-117.

[7] 李 卉,李德仁,黃先鋒,等.一種漸進加密三角網(wǎng)LIDAR點云濾波的改進算法[J]，測繪科學,2009,34(3):39-41. (LI Hui, LI De-ren, HUANG Xian-feng,etal. An Algorithm for LIDAR Point Cloud Filtering Based on Progressive Encrypting Triangulation[J]. Science of Surveying and Mapping, 2009, 34(3):39-41.(in Chinese))

[8] 崔建軍,隋立春,徐花芝,等.基于邊緣檢測算法的LiDAR數(shù)據(jù)建筑物提取[J].測繪科學技術學報，2008,25(2):98-100. (CUI Jian-jun, SUI Li-chun, XU Hua-zhi,etal. Building Extraction from the LiDAR Data Based on the Edge Detection Algorithm[J]. Journal of Zhengzhou Institute of Surveying and Mapping, 2008, 25(2): 98-100. (in Chinese))

[9] 任自珍,岑敏儀,張同剛,等.基于等高線形狀分析的LiDAR建筑物提取[J].西南交通大學學報,2009,44(1):83-88. (REN Zi-zhen, CEN Min-yi, ZHANG Tong-gang,etal. Building Extraction from the LiDAR Data Based on Contour Shape Analysis[J]. Journal of Southwest Jiaotong University, 2009,44(1):83-88.(in Chinese))

[10] NICHOLAS S,TAKIS K.Triangulated,Connected Sets for Building Detection from Irregularly Spaced LiDAR[C]∥ Proceedings of the 3rd International Symposium on Communications, Control and Signal Processing, Malta, March 12-14, 2008: 560-565.

[11] 沈 蔚,李 京,陳云浩,等. 基于LIDAR數(shù)據(jù)的建筑輪廓線提取及規(guī)則化算法研究[J].遙感學報, 2008,(5):692-698.(SHEN Wei, LI Jing, CHEN Yun-hao,etal. Algorithm Research Based on Buildings Contour Extraction and Regularization from LIDAR Data[J]. Journal of Remote Sensing, 2008,(5):692-698.(in Chinese))

[12] 王宗躍,馬洪超,徐宏根,等.基于LiDAR 點云數(shù)據(jù)的水體輪廓線提取方法研究[J]， 武漢大學學報·信息科學版,2010,35(4) :432-435. (WANG Zong-yue, MA Hong-chao, XU Hong-gen,etal. Research of Water Contour Extraction Based on LiDAR Point Cloud[J]. Geomatics and Information Science of Wuhan University, 2010, 35(4) :432-435. (in Chinese))

[13] LEE D, LIN A. Generalized Delaunay Triangulation for Planar Graphs[J]. Discrete and Computational Geometry, 1986, (1): 201-217.

[14] LEE D, SCHACHTER B. Two Algorithms for Constructing a Delaunay Traingulation[J]. International Journal of Parallel Programming, 1980，(3): 219-242.

[15] SLOAN S. A Fast Algorithm for Constructing Delaunay Triangulations in the Plane[J]. Advances in Engineering Software, 1987, (1): 34-55.

[16] DE FLORIANI L, PUPPO E. An On-line Algorithm for Constrained Delaunay Triangulation[J]. CVGIP: Graphical Models and Image Processing, 1992, (4): 290-300.

[17] SCHWEIKERT D. An Interpolation Curve Using a Spline in Tension[R]. USA: Department of Applied Mathematics of Brown University, 1965.