機載LiDAR數(shù)據(jù)提取堤防工程特征信息圓環(huán)探測法

沈定濤，錢天陸，夏 煜，陳蓓青，張 煜，王結臣

1. 長江水利委員會長江科學院流域水資源與生態(tài)環(huán)境科學湖北省重點實驗室, 湖北 武漢 430010; 2. 南京大學地理與海洋科學學院江蘇省地理信息技術重點實驗室, 自然資源部國土衛(wèi)星遙感應用重點實驗室, 江蘇 南京 210023; 3. 長江水利委員會長江科學院空間信息技術應用研究所, 湖北 武漢 430010

洪水災害是最嚴重的自然災害之一，全球每年都因洪水事件造成大量的人員傷亡和財產(chǎn)損失[1]。堤防是抵御洪水危害的重要工程措施，近幾十年來，通過大規(guī)模的堤防建設，我國已建成堤防總長度達到20余萬km，在歷次抗洪中發(fā)揮了巨大的作用[2]。堤防通常由兩部分組成，即堤頂和堤坡，其中堤頂是指堤防頂部的平坦表面，其橫向?qū)挾缺鹊袒唐率侵傅谭赖膬蓚炔糠?，其外邊緣與地面持平，堤坡與堤頂組成堤身。堤防一側通常是臨水的，以抵御洪水侵襲，另一側為人類活動區(qū)域[3]。堤防的規(guī)格、形制如堤防最小高程、堤頂寬度、堤坡坡度等都有一套設計規(guī)范，我國于2013年發(fā)布的《堤防工程設計規(guī)范》即明確規(guī)定了不同等級堤防的建設標準[2]。

對堤防開展定期巡檢，提取堤防工程特征要素是加強洪水保護的一項重要工作。堤防要素包括堤防中心線、堤頂、堤坡、斷面線等參數(shù)[4]，這些數(shù)據(jù)是開展堤防工程信息化[5-6]、堤防安全評估[7-9]的重要數(shù)據(jù)源。傳統(tǒng)的使用GNSS技術開展堤防巡檢需要到堤防現(xiàn)場對堤防上各個特征點位的經(jīng)緯度和高程進行測繪，這項工作相當耗時費力，且對技術人員的工程經(jīng)驗要求很高[10-14]。機載LiDAR技術的快速發(fā)展為大范圍、高精度的堤防巡檢工作提供了便利。使用機載LiDAR點云生成的高分辨率DTM與GNSS測繪得到的點坐標不同，它提供的是連續(xù)的地表形態(tài)數(shù)據(jù)，基于DTM可以生成各種高精度數(shù)據(jù)產(chǎn)品[15-18]，如坡度圖、曲率圖、等值線、斷面線等。

研究人員嘗試將機載LiDAR數(shù)據(jù)用于堤防的研究，如堤防制圖[19]，堤頂斷面提取[20]，堤防中心線提取[21]，堤頂?shù)唐率噶窟吔缇€提取[22-23],堤防侵蝕及穩(wěn)定性分析[3]等，但這些方法大多對堤防數(shù)據(jù)完整性要求較高[3]，一些研究需要大量繁雜的人工輔助操作[21,24]。以堤防中心線提取為例，其基本思路是沿著堤防方向，找出一組相鄰DTM網(wǎng)格中心點連接而成的矢量線，且要求DTM網(wǎng)格數(shù)目盡可能少，而網(wǎng)格高程值盡可能大[21]。最小成本路徑分析法是較常見的一種堤防中心線提取方法[20,25]，該方法在堤防兩端設置起始點和終止點，以堤防高程為特征計算柵格成本距離，再沿著堤防方向在堤防起點和終止點之間基于成本路徑分析查找堤防最高高程點集合，得出起始點和終止點之間的最優(yōu)柵格路徑。最小成本路徑分析法較為簡便快捷，但其缺點同樣明顯：一是由于該方法是一種窮舉算法，計算效率較低，對計算機配置要求較高；二是一旦堤防存在破壞出現(xiàn)凹凸不平現(xiàn)象，將導致提取的柵格路徑出現(xiàn)“毛刺”現(xiàn)象，需要進行后處理平滑曲線。為此，文獻[21]又提出了矩陣法、垂線法和點匹配垂線法用于堤防中心線提取。在矩陣法中，通過橫向和縱向掃描DTM網(wǎng)格，提取行(列)高程值最大柵格單元，組成堤防中心線。該方法計算效率較高，只用完成兩次掃描即可提取所有柵格點，但是其仍存在諸多缺陷，一旦堤防存在蜿蜒曲折，一次橫向和縱向掃描會多次穿過堤防堤頂區(qū)域，導致漏掉部分柵格點，同時，對于近乎水平或垂直的堤防區(qū)域提取效果較差。在垂線法中，通過事先在堤防單側設置等間隔特征點，從特征點對堤防縱向畫垂直射線，統(tǒng)計射線穿過DTM中的最大高程所在的柵格，該方法仍難以避免射線多次穿過堤防堤頂區(qū)域的問題。點匹配垂線法對其進行了改進，通過人工在堤防兩側設置數(shù)目相同的等間距特征點，特征點之間兩兩組合構建穿過堤防的垂線，通過垂線與堤防相交統(tǒng)計高程值最高的柵格構建中心線。垂線法和點匹配垂線法需要開展大量的前處理操作，如特征點選取、垂線生成等，工作量較大。

基于上述原因，本文使用機載LiDAR生成的DTM，嘗試對堤防中心線、堤頂?shù)唐滦螒B(tài)和堤防斷面提取分別開展研究，并提出一整套的技術和方法。最后，將本文方法在洞庭湖區(qū)共雙茶垸蓄滯洪區(qū)堤防開展實際應用，驗證所提出方法的有效性，這對加強堤防保護，防范洪水風險具有很實際的應用意義。

1 研究區(qū)和數(shù)據(jù)

洞庭湖位于長江流域中游，是長江最重要的洪水調(diào)蓄湖泊(圖1)。洞庭湖區(qū)經(jīng)過20世紀50年代堤垸建設、60年代電排站建設以及70年代防洪配套治理工程建設，目前已形成大小堤垸226個，防洪大堤總長度5812 km。洞庭湖區(qū)堤防建有水閘、涵管5000多處，電排站3000多座，這些穿堤設施大多修建于20世紀50、60年代，老化嚴重，運行不良，對堤防安全構成嚴重威脅，防汛任務巨大。本文研究區(qū)共雙茶垸蓄滯洪區(qū)防洪干堤為一閉合環(huán)形結構(圖2)，總長度120.1 km。共雙茶垸蓄滯洪區(qū)位于沅江市北部，北臨草尾河，南臨黃土包河，東面南洞庭湖、西面隔水與赤山相望，四面環(huán)水，為洞庭湖區(qū)24個大型蓄滯洪區(qū)之一。

圖1 共雙茶垸蓄滯洪區(qū)在洞庭湖區(qū)的位置Fig.1 The location of the Gongshuangcha detention basin in the Dongting Lake area

本文使用德國TopoSys公司HARRIER 68i機載激光測量系統(tǒng)對共雙茶垸開展了航空攝影測量工作。機載激光測量系統(tǒng)中數(shù)碼相機使用的是Trimble公司的Rollei Metric AIC Pro(像素6000萬)；慣導系統(tǒng)型號為Applanix POS/AV系列，采樣頻率200 Hz；激光掃描儀的型號為Riegl LMS-Q680i，最大脈沖頻率80 KHz～400 KHz，掃描角度45°/60°。通過對采集的機載LiDAR點云進行處理，生成了1 m高分辨率堤防DTM數(shù)據(jù)(圖2)，平面位置中誤差為0.44 m，高程中誤差為0.04 m。

圖2 蓄滯洪區(qū)影像及堤防DTM數(shù)據(jù)Fig.2 Image of flood detention basin and DTM data of levee

2 特征信息提取方法

2.1 堤防中心線提取

2.1.1 圓環(huán)探測法基本原理

現(xiàn)有堤防中心線自動提取方法對堤防形態(tài)完整性要求較高，然而堤防通常會受到不同程度的破壞，使得堤防形態(tài)發(fā)生變化，堤防中心線提取難以獲得較滿意的效果。如圖3所示，以本文研究區(qū)為例，堤防受到洪水侵蝕、工程建設等的影響，很難將上述方法用于本文研究區(qū)堤防中心線的自動提取。

圖3 由于人類活動、侵蝕等造成堤防破壞Fig.3 Levee damage due to human activities, erosion, etc

圖4是一段堤防示意圖，堤頂區(qū)域高程較堤坡區(qū)域高，堤防中心線即是要從堤頂區(qū)域提取高程值最大的系列柵格點。如圖5所示，堤防“中間高，兩側低”的形態(tài)特征使得DTM中間區(qū)域網(wǎng)格高程高于邊緣區(qū)。假定以一圓環(huán)覆蓋于堤防上，并探測圓環(huán)穿過的所有DTM網(wǎng)格高程，查找高程值最大的柵格，基于堤防形態(tài)特征，可以確定該柵格單元是屬于圖5位于堤頂區(qū)域的紅色柵格。

圖4 堤防示意圖Fig.4 Levee diagram

圖5 圓環(huán)與堤防DTM相交計算Fig.5 Intersection of ring and levee DTM

如圖6所示，若將查找到的DTM柵格中心(B點)作為圓心點繪制新的圓環(huán)，該圓環(huán)與圓環(huán)A相交，那么可以在位于圓環(huán)A外部的圓環(huán)B上，查到高程值最大的柵格，該網(wǎng)格中心為點C處于堤頂區(qū)域。以此類推，在堤防DTM上查找高程值最大的柵格，依次連接柵格中心點A、B、C、D、E等形成的矢量線，即是堤防中心線，可以將此方法稱為圓環(huán)探測法。

圖6 圓環(huán)探測堤防DTM高點Fig.6 Ring detection of high point on levee DTM

2.1.2 中心線提取方法

因為堤防平均寬度w可以看作是已知的，圓環(huán)半徑r的設定則相當重要。半徑r過小，會出現(xiàn)圓環(huán)完全處于堤防DTM內(nèi)部，或者圓環(huán)只與堤防一側相交，此時可能選不到堤防中心線上的點，也可能出現(xiàn)圓環(huán)之間形成“回路”，即新圓環(huán)與之前的圓環(huán)相交，導致計算失敗。若半徑r過大，在堤防大幅變向處，圓環(huán)探測會漏掉堤防拐彎處的“尖點”。因此，有必要分析堤防平均寬度w與圓環(huán)半徑r的關系，以盡量規(guī)避上述情形，進而實現(xiàn)堤防中心線的自動提取。

(1) 確定圓環(huán)半徑r。設定堤防平均寬度為w，則圓環(huán)半徑r必須滿足如下條件才能獲得滿意的計算效果

(1)

圖7 圓環(huán)半徑r設定Fig.7 Ring Radius r setting

(2) 確定圓環(huán)離散點n。假定當前圓環(huán)中心O點的坐標(xo,yo)，圓環(huán)半徑r，DTM柵格單元寬度為c。若將圓環(huán)離散化為一系列圓上點集，只要點足夠密集，就可以通過圓上點查找高程值最大的柵格。為避免漏掉圓環(huán)穿過的柵格，必須使得圓上相鄰離散點距離不超過一個柵格寬，則相鄰離散點的距離l需滿足l≤c，又已知圓環(huán)周長為2πr，若將圓環(huán)n等分，則圓環(huán)上相鄰離散點的距離l≈2πr/n，即必須滿足2πr/n≤c，進而可推導出n≥2πr/c。以本文研究區(qū)為例，DTM網(wǎng)格寬度c=1 m，圓環(huán)半徑r=50 m，則n≥314。為便于計算，將圓環(huán)360°等分，即n=360，此時圓環(huán)上相鄰離散點之間距離l≈2πr/n=0.87 m，完全滿足圓環(huán)探測需要。

如圖8所示，以圓環(huán)中心點O為坐標原點，以水平方向為x軸，垂直方向為y軸，從x軸正方向上作為圓環(huán)第一個點開始，依次按照1度遞增，形成圓環(huán)點集，共計360個點。那么，對點集中的A點而言，其與x軸正方向的角度α是已知的，設圓環(huán)中心O點的坐標為(xo,yo)，則點A的空間坐標xA，yA滿足式(2)

圖8 求堤防中心線上點A坐標Fig.8 Calculate the coordinate of point A on the central line of levee

(2)

由于DTM網(wǎng)格寬度為c，假設DTM范圍最小最大坐標為(xmin,ymin)、(xmax,ymax)，則點A對應的DTM網(wǎng)格行列號IA，JA滿足式(3)

(3)

按照式(2)和式(3)，查找圓O上所有離散點，找到高程值最大的點。如圖8所示，假定點A所在的DTM網(wǎng)格就是高程值最高網(wǎng)格，以點A為圓心，r為半徑做圓，因為圓心點A所在的圓和圓心點O所在的圓半徑都是r，則點A所在的圓必過點O。假定圓A和圓O相交于點P1和點P2，線段AO和線段P1P2相交于點P。要找到圓A上的某一堤防中心線上點B，則點B必須滿足式(4)

(4)

按照上述計算方法，依次得到所有堤防中心點，生成堤防中心線。

2.2 堤防形態(tài)參數(shù)提取

堤防監(jiān)測包括對堤防形態(tài)參數(shù)進行長期觀測，如堤防高度、堤頂、堤基寬度、坡度等，這些參數(shù)通常是使用GNSS測量獲得。堤防上每隔一定距離需要采集堤防斷面數(shù)據(jù)，斷面處通常是由斷面點構成，這些斷面點的間隔距離設定非常依賴于觀測人員的經(jīng)驗[12-14]。基于LiDAR數(shù)據(jù)提取的DTM數(shù)據(jù)有著較高的數(shù)據(jù)精度，用于堤防形態(tài)參數(shù)提取有著更好的效果，同時可以從堤頂和堤坡矢量面數(shù)據(jù)提取堤防輪廓線數(shù)據(jù)。

2.2.1 堤防堤頂、堤坡提取

獲取堤防每一區(qū)域的堤頂和堤坡參數(shù)信息是進行堤防破損檢測的重要基礎。按照《堤防工程設計規(guī)范》國家標準的規(guī)定，對土堤堤防而言，“1級堤防堤頂寬度不宜小于8 m，2級堤防堤頂寬度不宜小于6 m，3級堤防堤頂寬度不宜小于3 m；堤頂?shù)钠露葢蛞粋然騼蓚葍A斜，坡度宜采用2%～3%；1級、2級土堤的堤坡不宜陡于1∶3”(即坡度為18.43°)[2]。

基于堤防坡度區(qū)分堤頂和堤坡是較為常見的方式，如文獻[24]的研究中，堤防頂部坡度值范圍設定為[0°,8.43°]，堤坡坡度范圍設定為[8.43°,43.69°]。共雙茶垸蓄滯洪區(qū)主干堤防為3級堤防，由于該堤防建設時間較早，而《堤防工程設計規(guī)范》國家標準于2013年頒布，因此共雙茶垸蓄滯洪區(qū)干堤并不能完全滿足該標準對3級堤防的要求。按照國家標注對1級、2級堤防堤坡不陡于18.43°的規(guī)定，參照文獻[24]的研究將本堤防堤坡的坡度值容許加減±10°的誤差，即堤坡坡度為[8.43°,28.43°]，同理，國家標準對堤頂?shù)钠露纫鬄?%～3%，此處考慮堤防的建設現(xiàn)狀，將堤防堤頂坡度設定為[0°,8.43°]之間。

提取堤頂和堤坡形態(tài)參數(shù)的第1步即是將堤防DTM數(shù)據(jù)生成坡度圖，如圖9(a)所示。該坡度圖中最大坡度為68.87°，最小坡度為0°。按照上述對堤頂和堤坡的坡度范圍分析，將坡度圖中坡度分為3類：[0°,8.43°]、[8.43°,28.43°]和[28.43°,68.87°]，得到的坡度分類圖如圖9(b)所示。

最后，將坡度分類圖矢量化，得到坡度分類圖矢量數(shù)據(jù)，如圖9(c)所示。圖9中出現(xiàn)了較多的碎圖斑，這些碎圖斑面積一般小于100 m2(10 m×10 m)，多數(shù)是因為堤頂和堤坡上由于多年侵蝕和人工建設造成，使用ArcGIS的eliminate工具，將那些面積較大的圖斑內(nèi)部面積低于100 m2碎圖斑融合到大圖斑中，得到結果如圖9(d)所示。

因為在裁剪堤防DTM的過程中，堤防外側一些農(nóng)田等非堤防區(qū)域也包含在了堤防DTM中，在生成坡度圖時，這些區(qū)域的坡度也涵蓋在[0°,8.43°]、[8.43°,28.43°]坡度范圍中。如圖9所示，在局部放大圖中，位于堤防外側的品紅色圖斑與堤頂圖斑都屬于坡度范圍為[0°,8.43°]，該圖斑實際上不是堤頂，是處于堤防外部的農(nóng)田。因此，針對上述矢量化成果，需要剔除堤防外側的這類圖斑，以及坡度范圍為[28.43°,68.87°]的圖斑，得到最終包含堤頂和堤坡的堤防形態(tài)矢量數(shù)據(jù)如圖9(e)所示。

圖9 堤防堤頂、堤坡提取Fig.9 Extraction of levee crown and slope

2.2.2 堤防斷面線提取

堤防斷面數(shù)據(jù)是堤防工程最重要的特征信息，劃分堤防斷面，分析斷面處的堤防高程變化，可以對堤防橫切面形態(tài)有更加直觀的認識。本文在上述提取堤防中心線過程中按照50 m等間距對堤防劃分圓環(huán)，可以直接在每個圓環(huán)圓心處劃分斷面。如圖10所示，線段MN即是需要自動提取的斷面線，該線段是∠OAB的角平分線，M點和N點都位于以A為圓心的圓環(huán)上，即線段MN的長度為2r。只要求出點M和點N的空間坐標，即可得到穿過圓心A的斷面線MN。

如圖10，已知點O、點A和點B的空間坐標分別為(xO,yO)、(xA,yA)、(xB,yB)，因為點M位于圓環(huán)A上，且點M到點O和點B的距離相等，即點M滿足

圖10 堤防斷面線自動提取Fig.10 Automatic extraction of section lines of levee

(5)

按照前述圓環(huán)離散化方法，對圓環(huán)A上離散點進行遍歷，結合式(2)、式(3)、式(5)，必定可以求得點M的坐標，計算點N的坐標同理。對于堤防上任一斷面線，探測斷面線穿過的堤防DTM網(wǎng)格即可獲得斷面處堤防高程變化曲線。

3 結果與分析

3.1 堤防中心線結果分析

圖11和圖12是使用上述方法得到的堤防圓環(huán)和堤防中心線數(shù)據(jù)及其局部放大。與最小成本路徑法等逐柵格采樣方式相比，本文方法在數(shù)據(jù)采樣精度上存在一定精度損失，但本文方法對堤防形態(tài)變化等敏感度較低，針對圖3的各類堤防“異?！鼻樾我部梢缘玫捷^滿意的中心線數(shù)據(jù)。

圖11 圓環(huán)探測及中心線生成Fig.11 Ring detection and centerline generation

圖12 中心線局部放大Fig.12 Magnification of centerline

本文提出的圓環(huán)探測法有效規(guī)避了各種堤防異常問題，如圖13所示，堤防由于破損等情形的存在導致DTM中網(wǎng)格高程異常。其中的暗色斑點屬于異常高程點，這對如矩陣法、最小成本路徑法等逐柵格提取中心線的方法易使得計算結果效果較差。本文方法在堤防上間隔采樣，可以在一定程度上避開這些局部高程異常的網(wǎng)格，盡量保證堤防中心線處于堤防縱向中心位置。在一些堤防破損非常嚴重的區(qū)域，有少數(shù)這種高程異常凸起的網(wǎng)格點被采樣作為堤防中心線點，這些偏離堤防中心區(qū)域的異常點也表明了此處堤防破壞較為嚴重。

圖13 堤防局部高程異常Fig.13 Abnormal elevation of local dike

3.2 堤防堤頂、堤坡結果分析

盡管堤防在正常區(qū)域能夠提取出正確的堤頂、堤坡形態(tài)參數(shù)。但是，因為受到水土侵蝕、人類活動等影響，仍然存在部分堤防堤頂和堤坡形態(tài)被破壞的情形，如圖14所示。在圖14(a)中，因為在堤頂存在建筑物，堤防堤頂多邊形出現(xiàn)“拐彎”現(xiàn)象，且建筑物區(qū)域的坡度也處于堤頂坡度范圍，同時導致堤防中心線也出現(xiàn)“拐彎”現(xiàn)象；在圖14(b)中，因為人類活動和水土侵蝕，堤頂寬度出現(xiàn)劇烈變化，堤頂過寬或過細，這是堤防缺乏養(yǎng)護的典型表現(xiàn)；而在圖14(c)中，由于人類活動對堤防破壞太大，不僅出現(xiàn)堤頂寬度過寬，而且導致堤防中心線變道，堤防在縱向上，出現(xiàn)多條堤頂和堤坡區(qū)域。

圖14 堤防堤頂、堤坡形態(tài)受到破壞Fig.14 Damaged levee crown and slope

3.3 堤防斷面結果分析

基于上述求斷面線端點坐標式(4)，得到的堤防斷面線如圖15所示，研究區(qū)堤防共生成斷面2401個。

圖15 堤防斷面線提取結果Fig.15 Extraction results of levee section lines

堤防破損對斷面高程的影響也不容忽視，在堤防破損比較嚴重的區(qū)域堤頂寬度和堤坡坡度都發(fā)生了顯著的變化。如圖16所示，圖(a)為正常堤防區(qū)域某斷面的高程變化，圖(b)為異常堤防區(qū)域某斷面的高程變化。本文研究區(qū)堤防堤頂高程在35 m左右，在正常堤防區(qū)域堤頂部分高程值與堤坡高程值有很明顯的對比，堤頂寬度也符合堤防設計規(guī)范，而在異常堤防區(qū)域，由于堤防受到較嚴重的破壞，導致堤防堤頂區(qū)域不明顯，整個橫斷面有一大片區(qū)域高程值都達到35 m，導致堤坡窄而陡。

圖16 堤防斷面高程變化對比Fig.16 Comparison of elevation changes of levee cross section

4 討 論

4.1 與幾類堤防中心線生成方法的對比分析

圖17是文獻[21]提出的堤防中心線生成方法，其中矩陣法是自動算法，垂線法和點匹配垂線法是半自動算法。圖17(a)為矩陣法示意圖，掃描線EF實際跨越了兩個堤防中心線點，但是算法只能提取其中一點，因此對于堤防彎曲較大的區(qū)域存在漏掉關鍵點的可能。同時，當局部堤防區(qū)域過于水平或垂直，掃描線很難提取到所有中心線點。如圖中掃描線CD幾乎與堤防中心線重合，再加上此區(qū)域橫向掃描存在跨越兩次堤防中心線的情形，使得中心線點提取存在困難。

如圖17(b)和圖17(c)所示，垂線法是在堤防邊緣單側取點，點匹配垂線法是在堤防邊緣雙側取點，然后手工繪制垂線，通過垂線與堤防相交查找中心線點。該方法較為煩瑣，需要大量的輔助操作，以點匹配垂線法為例，若用于本文研究區(qū)堤防中心線提取，假定按照本文研究設定的50 m取點，需要在120 km長的堤防雙側各選取2400個點構建垂線段。

圖17 3種堤防中心線生成方法[21]Fig.17 Three generating methods of levee center line[21]

圖18和圖19分別是使用矩陣法和最小成本路徑法對本文研究區(qū)堤防數(shù)據(jù)提取中心線。采用矩陣法提取的最大高程值柵格在空間上很難做到相互連接，若直接執(zhí)行柵格轉(zhuǎn)換矢量會生成大量斷線，因此該方法仍需要大量的后處理操作才能得到滿意的堤防中心線數(shù)據(jù)。使用最小成本路徑法獲取的堤防中心線雖然沒有矩陣法出現(xiàn)的散亂斷線問題，但是存在大量的“毛刺”現(xiàn)象，需要進行平滑后處理操作。

圖18 矩陣法生成堤防中心線Fig.18 Generation of levee center line by matrix method

圖19 最小成本路徑法生成堤防中心線Fig.19 Generation of levee center line by cost path analysis

圖20為采用最小成本路徑法、矩陣法和圓環(huán)探測法得到的堤防中心線結果對比。由于最小成本路徑法需要計算路徑最短、成本最小(高程最高)的柵格集合，因此在堤防彎曲較大的區(qū)域，最小成本路徑法得到的堤防中心線出現(xiàn)“抄近路”現(xiàn)象而漏掉堤防中心線關鍵點。由于圓環(huán)探測法中最大高程點探測是在圓環(huán)上，一旦堤防“尖點”位于圓環(huán)內(nèi)部，則同樣存在漏掉關鍵點的可能。一種有效的解決方法是對圓環(huán)進行加密探測，如圖21所示，A、B、C為堤防中心線點，由于點A和點C位于圓環(huán)B的360個環(huán)上點次序是可知的，根據(jù)對稱性，可以計算得出平分線段AC且穿過點B的圓環(huán)上點E、H，進而分別計算得出弧AH、HC、AE和EC的二等分點G、I、D、F。可以類似點匹配垂線法，計算線段GD和線段FI穿過堤防區(qū)域的最大高程點，分別統(tǒng)計最大高程點與線段AB和BC的垂直距離，通過預先設定的距離閾值確定是否保留該最大高程點作為堤防中心線，從而可以避免漏掉堤防關鍵點。

圖20 堤防中心線結果對比Fig.20 Comparison of levee center lines

圖21 圓環(huán)加密探測示意圖Fig.21 Schematic diagram of ring encryption detection

4.2 與Bresenham繪圓算法的對比分析

Bresenham繪圓算法是計算機圖形學領域的經(jīng)典算法，該算法是為解決幾何圖形的計算機顯示問題而設計的。本文嘗試使用Bresenham繪圓算法提取堤防中心線，并與本文所提出的方法進行對比分析。由于Bresenham繪圓算法是基于“點陣逼近”模擬圓形，確定點陣間距非常重要。為了避免圓環(huán)在DTM網(wǎng)格上探測時漏掉柵格，必須保證圓環(huán)上相鄰點間隔不超過1個網(wǎng)格。本文中DTM網(wǎng)格空間分辨率為1 m，為此將Bresenham繪圓算法點陣間距設置為1 m可以滿足上述要求。圖22是使用Bresenham繪圓算法繪制的半徑為50 m的圓形，共288個點。與Bresenham繪圓算法不同，本文所使用的圓環(huán)繪制方法屬于“多邊形逼近”方式繪圓，多邊形有360個點均勻分布在圓上。同時，Bresenham繪圓算法特點使得各點位于圓環(huán)內(nèi)側或外側，并不完全位于圓上，如圖22紅色圓弧所示，“多邊形逼近”繪圓效果明顯好于“點陣逼近”繪圓效果。

圖22 兩種繪圓算法對比Fig.22 Comparison of two circle drawing algorithms

設定DTM網(wǎng)格寬度為c，堤防長度為L，圓環(huán)半徑為r，可知圓環(huán)探測法所需圓環(huán)數(shù)量約為L/r，按照前述分析，圓環(huán)離散點n=2πr/c，則本文研究所采用的圓環(huán)探測法算法時間復雜度為T(n)=nL/r=2πL/c，即算法時間復雜度是關于堤防長度與DTM網(wǎng)格寬度比值的一個線性函數(shù)。受圓環(huán)相鄰點間隔必須控制在一個網(wǎng)格距離的約束，本文所采用的繪圓方法和Bresenham繪圓算法二者所需的圓環(huán)點數(shù)量比可以近似看作是一個常數(shù)，因此，兩種繪圓方法的時間復雜度是一致的。

圖23是兩種繪圓方法提取堤防中心線的局部結果對比，Bresenham繪圓鋸齒狀明顯。由于兩種方法的圓環(huán)采樣點選取原理不同，使得最大高程點柵格選擇會出現(xiàn)一定偏差，導致兩種圓環(huán)探測結果并不完全重合。本文所采用的圓環(huán)繪制方法僅基于圓環(huán)中心坐標和半徑即可快速計算環(huán)上任一點的坐標，且圓環(huán)點對稱性更好，用于圓環(huán)加密探測(圖21)和圓環(huán)斷面提取(圖10)時簡便快捷，比Bresenham繪圓更便于程序?qū)崿F(xiàn)。

圖23 兩種繪圓算法提取堤防中心線對比Fig.23 Comparison of two circle drawing algorithms for extracting levee centerline

5 結 論

利用機載LiDAR生成的高分辨率DTM提取各種堤防工程特征信息，有效克服了傳統(tǒng)測量方法的缺點，為堤防結構的特征監(jiān)測提供了更為精確和精細的產(chǎn)品。這些數(shù)據(jù)產(chǎn)品在堤防數(shù)據(jù)建庫、堤防管理信息系統(tǒng)和堤防穩(wěn)定性和安全評估方面都有著很好的應用潛力。