周維娜，程曉光，嚴明，盧浩男

(1.常州市自然資源和規(guī)劃服務(wù)中心，江蘇 常州 213022；2.飛燕航空遙感技術(shù)有限公司，南京 210001)

0 引言

機載激光雷達(light detection and ranging，LiDAR)已經(jīng)成為陸地地形測繪的主要手段之一。在機載LiDAR測繪內(nèi)業(yè)生產(chǎn)中，水體斷裂線是主要的斷裂線類型之一[1-2]?？梢院唵握J為，水體斷裂線是水體在二維平面上的邊界。除數(shù)字高程模型(digital elevation model，DEM)生產(chǎn)外，水體邊界還可用于地圖制圖、河湖蓄水監(jiān)測等領(lǐng)域。對于水體較多較小的測區(qū)，手工勾繪水體斷裂線繁瑣，工作量大。如何自動、準確地提取水體邊界用作斷裂線，對提高機載LiDAR內(nèi)業(yè)生產(chǎn)效率很有意義。

用于陸地測繪的機載LiDAR一般采用波長在近紅外波段的激光脈沖[3]。由于純凈水體在近紅外波段反射率極低[4]，一般后向散射回波功率很弱，難以被激光接收器識別為信號，因而在點云中，水體一般呈現(xiàn)為空洞或點分布特別稀疏的情況。這使得基于機載LiDAR點云提取水體邊界具備獨特的優(yōu)勢。但是，由于風(fēng)、水生植物、水質(zhì)等因素的影響，水體內(nèi)部經(jīng)常存在一定數(shù)目的點，可能分布較密集、回波強度較高，甚至與岸邊相連。

目前基于機載LiDAR點云提取水體邊界的研究一般基于水面點稀少、回波強度弱等假設(shè)[5-7]，導(dǎo)致難以處理水體內(nèi)點較多或回波較強的情況。有研究人員使用光學(xué)影像輔助水體邊界提取[8-10]。但是在航攝光學(xué)影像上，不同水質(zhì)、水深、水草和陰影的水體的特征不一致，且相鄰圖幅間可能存在色差，而衛(wèi)星遙感影像的分辨率一般偏低。此外，提取建筑物屋頂邊界等的方法[11]不適用于提取水體邊界。

本文提出一種基于Delaunay三角網(wǎng)提取靜態(tài)水體邊界的方法。該方法不依賴影像輔助、點云預(yù)分類和回波強度假設(shè)，不需要已知水面高程，只需要水體存在面積夠大的平面空洞，既適用于水面點較少的情況，也適用于水面點較多的情況。首先，在點云空洞提取的基礎(chǔ)上估計水面高程。然后，基于點云高程沿三角網(wǎng)進行邊界擴張，直至邊界不再變化。在擴張中，動態(tài)更新對水面高程的估計。最后，修剪邊界。本文采用江蘇常州小黃山地區(qū)的點云進行了實驗，成功提取到所有主要靜態(tài)水體的邊界，且與數(shù)字正射影像(digital orthophoto map，DOM)中的水體邊界套合良好。

1 靜態(tài)水體水面點高程分析

靜態(tài)水體多為湖泊、池塘，或流動緩慢的河流(不考慮潮汐河口或水體與海水聯(lián)通的情況)，短時間內(nèi)水位變化小。對一片聯(lián)通水體來說，液態(tài)水作為一種流體，在無其他外力作用時，會在重力的作用下流動，直至水體表面沒有重力勢能差。對中小面積的靜態(tài)水體，如果沒有明顯局部重力異常，可以簡單認為，除出入水口外，此時水面的橢球體高程基本一致。

但是，波浪、水生植物、水質(zhì)等因素可能造成靜態(tài)水體水面高程不一致。波浪產(chǎn)生的水面點間的高程差一般在厘米級甚至分米級。在水生植物中，挺水植物可以形成明顯高出水面的點云，而浮藻、浮葉植物、漂浮植物等產(chǎn)生的點高程十分接近水面高程。如果沉水植物離水面很近且水質(zhì)不佳，則有可能水體不能完全吸收激光脈沖，沉水植物產(chǎn)生回波點，所以，點云中可能存在低于水面的點。不同類型的水生植物形成的點云差異較大，分析困難。但考慮到浮藻、浮葉植物、漂浮植物等的點高程接近水面高程，沉水植物的回波點少且容易被排除，而茂密的挺水植物經(jīng)常遮蓋水面，使人難以判斷水體邊界，所以在估計水面高程范圍時，本文不考慮水生植物的影響。

假設(shè)LiDAR測量和波浪波動相互獨立，則可以認為，水面高程的方差Vwater由兩部分構(gòu)成：LiDAR測量相對誤差的方差Vmeas和波浪造成的方差Vwave，見式(1)。

Vwater=Vmeas+Vwave

(1)

2 方法

本文提出了一種沿Delaunay三角網(wǎng)迭代擴張，動態(tài)估計水面高程的靜態(tài)水體邊界提取方法(圖1)?？紤]到水體內(nèi)的島嶼、設(shè)施等會使水體內(nèi)部存在空洞，將單個靜態(tài)水體抽象為OpenGIS簡單特征規(guī)范[12]中的多邊形(Polygon)，允許空洞存在。

圖1 本文方法流程圖

2.1 點云去噪

為避免高程異常點影響水面高程估計，對點云進行直通濾波，只保留高程位于合理范圍內(nèi)的點。

2.2 初始邊界計算

在計算水體初始邊界時，本文對文獻[13]的方法進行了修改，拋棄點云外邊界，只留下點云空洞邊界。設(shè)邊長閾值為Lthre，面積閾值為Athre。在得到點云空洞邊界后，構(gòu)建帶拓撲關(guān)系的水體初始邊界，規(guī)定外邊界點逆時針排列，空洞邊界點順時針排列。

此處定義：設(shè)邊界上有P0、P1、…、PN共N+1個頂點，P0=PN，則稱PiPi+1(0≤i≤N-1)為邊界段；如果PiPi+1是由ΔPiPkPi+1搜索到，則ΔPiPkPi+1的鄰接三角形中，Pk對著的那個三角形為搜索方向。

2.3 水面高程估計

假設(shè)機載LiDAR對同一個高程的觀測值服從正態(tài)分布。通過統(tǒng)計測區(qū)內(nèi)局部平坦地面的點云高程，可以大致估計Vmeas和對應(yīng)的標準差σmeas。Vwave難以實測，但與水面波動幅度正相關(guān)。根據(jù)式(1)，Vwater≥Vmeas，所以σwater≥σmeas，其中σwater是水面高程的標準差。

本文認為水體邊界上含有水面點，可以由外邊界和空洞邊界估計水面高程范圍。和落在水體周圍的植被、建筑物、地面上的點相比，水面點高程相對偏低。極少量來自沉水植物的點，則明顯低于水面點的正常高程范圍。為確定水面高程的可能范圍，采用聚集層次聚類[14]對邊界點高程進行分類。當(dāng)類間最小距離均大于ds時，聚類終止。對每個高程類，記錄高程最小值、最大值、高程數(shù)和點數(shù)。

在高程分類后，對所有高程類按照最小值從小到大的順序進行遍歷。最先滿足高程數(shù)大于等于Nzthre的高程類被認為含有水面點高程。如果沒有高程類滿足條件，可適當(dāng)增大ds，進行二次分類。如果二次分類還沒有滿足條件的高程類，則認為該邊界不是有效的水體邊界，拋棄該邊界。設(shè)滿足條件的高程類的最小值為Zmin，最大值為Zmax。在[Zmin，Zmax]內(nèi)，值越小，是水面高程測量值的可能性越大，所以設(shè)置水面高程范圍下限Zwmin=Zmin。為確定水面高程范圍上限Zwmax，采用如下方法。

步驟1：對[Zmin，Zmax]內(nèi)的每一個高程值Zt(受限于存儲精度，Zt數(shù)量有限)，計算高程在[Zmin，Zt]內(nèi)的邊界點的平均高程Zmean和高程標準差σt。

步驟2：取邊界點數(shù)大于5且σt≥σwater的最小Zt為Zwmax。如果沒有Zt滿足條件，則取Zwmax=Zmax。設(shè)Zwmax對應(yīng)的Zmean為Zwmean。

2.4 水體空洞與建筑物空洞的區(qū)分

在機載LiDAR點云中，除了水體造成的空洞外，還經(jīng)常存在建筑物造成的空洞(圖2)。后者往往一側(cè)為屋頂，另一側(cè)為更低的屋頂或地面。上述特征使得建筑物空洞的邊界點高程分布直方圖呈現(xiàn)明顯的雙峰特點。

圖2 建筑物空洞示例

使用如下方法區(qū)分建筑物空洞和水體空洞。

步驟1：采用2.3節(jié)的方法對空洞邊界點進行高程分類。

步驟2：尋找高程數(shù)最多的2個高程類。

步驟3：如果2個高程類的點數(shù)均超過空洞邊界點總數(shù)的1/3，且較低高程類的最大值低于較高高程類的最小值2.5 m以上，則認為該空洞是建筑物空洞，否則是水體空洞。

2.5 基于高程的邊界擴張

在得到水體初始邊界后，對每個水體沿外邊界和空洞邊界進行基于高程的邊界擴張。首先定義，如果某三角形三個頂點高程均在[Zwmin，Zwmax]內(nèi)，則該三角形為水平三角形，否則為非水平三角形。在單次擴張中，沿著搜索方向試圖擴張一個三角形的距離，將水平三角形擴張到邊界內(nèi)部。設(shè)當(dāng)前對邊界段PiPi+1進行擴張，方法如下(參考圖3)。

步驟1：刪除已擴張邊界的最后一個頂點。

步驟2：計算頂點含Pi的三角形的集合SΔ，以Pi為中心逆時針排列。

步驟3：篩選SΔ中位于PiPi+1和已擴張邊界的最后一段PlPi之間，而且包含搜索方向的三角形集合St，以Pi為中心順時針排列。

注：白箭頭表示搜索方向。在(b)～(i)中，黑箭頭表示擴張后的有序邊界；白三角形為水平三角形；灰三角形為非水平三角形。圖3 基于高程的邊界擴張示意圖

圖3(b)至圖3(i)顯示了在圖3(a)St中的3個三角形的8種組合(是否水平三角形)對應(yīng)的擴張邊界。對所有邊界段擴張一遍，并排除無效邊界，稱為一輪擴張。對邊界進行多輪擴張，至邊界不再變化，則基于高程的邊界擴張終止。在每輪擴張后，都需要動態(tài)更新對水面高程的估計。

2.6 交叉多邊形處理

在邊界擴張的過程中，擴張后的多個邊界之間可能發(fā)生交叉。例如，同一片水體被淺灘分為多個部分，每個部分的邊界擴張到越過淺灘導(dǎo)致邊界間發(fā)生交叉。單個邊界可能發(fā)生自交叉，而自交叉多邊形不符合OpenGIS簡單特征規(guī)范。為此，需要對交叉多邊形(包括多個多邊形間交叉和單個多邊形自交叉)求并集。圖4為示意圖。圖中，粗閉合實線為一個凹多邊形的邊界(粗花箭頭表示搜索方向，實粗箭頭表示邊界點順序，虛粗箭頭表示邊界點逆序)；細閉合實線為一個帶洞多邊形的邊界(細花箭頭表示搜索方向，實細箭頭表示邊界點順序，虛細箭頭表示邊界點逆序)；1～6是候選多邊形編號。首先，對邊界采用文獻[15]的方法進行分割。

圖4 交叉多邊形處理示意圖

步驟1：對交叉邊界按照交點分割為正序邊界段。

步驟2：對正序邊界段生成對應(yīng)的逆序邊界段。

步驟3：使用最小角準則連接所有邊界段(包括正序和逆序邊界段)，生成的多邊形稱為候選多邊形。

由圖4可以發(fā)現(xiàn)，1的邊界是并集外邊界，搜索方向朝外，逆時針排列；4和7的邊界是并集空洞邊界，搜索方向朝內(nèi)，順時針排列。本文規(guī)定：邊界點逆時針排列且搜索方向朝外的候選多邊形的邊界為并集外邊界，邊界點順時針排列且搜索方向朝內(nèi)的候選多邊形的邊界為并集的空洞邊界。此處只保留這兩種候選多邊形的邊界。

兩幫劈裂示意類似于圖1a與圖1b，在巷道開挖前，兩幫內(nèi)有原生節(jié)理裂隙及采動產(chǎn)生的裂隙，開挖后在強卸載作用下節(jié)理裂隙很快受剪切、拉伸破壞并貫通，發(fā)生剪脹變形，如果受到動載沖擊波沖擊，則這一過程加劇。貫通后的節(jié)理裂隙繼煤體深部延伸，剪脹變形加劇，直發(fā)生大變形甚至失穩(wěn)。

2.7 孤立點邊界移除

由于水面高程估計不一定非常精確，所以可能有少量水面點的高程不在[Zwmin，Zwmax]內(nèi)。這導(dǎo)致在基于高程的邊界擴張后，可能出現(xiàn)一些孤立點造成的邊界(圖5)。對這些孤立點產(chǎn)生的邊界，需要將它們排除。這些邊界的頂點相對較少，邊界段搜索方向朝內(nèi)，且邊界內(nèi)部由環(huán)繞孤立點的三角形構(gòu)成。基于以上特征，使用如下方法進行處理。

注：白點為高程在[Zwmin，Zwmax]內(nèi)的點；黑點為高程不在[Zwmin，Zwmax]內(nèi)的點；粗黑線為基于高程擴張后的邊界。圖5 孤立點邊界示意圖

步驟1：如果某邊界的邊界段數(shù)目小于等于Eiso，且搜索方向朝內(nèi)，則遍歷邊界段，進入步驟2。

步驟2：對邊界段PiPi+1，設(shè)代表其搜索方向的相鄰三角形為Δneib，查找在Δneib中除Pi和Pi+1外的另一個頂點Pv。如果Pv不在邊界上，則將Pv放入集合SV。

步驟3：遍歷完后，如果SV中不重復(fù)頂點的數(shù)目小于等于Niso，則認為邊界由孤立點造成，拋棄該邊界，否則保留該邊界。

2.8 邊界修剪

上述步驟得到的水體邊界可能存在重合邊界點(非首尾點重合)造成的突出部分(簡稱突部)和內(nèi)陷部分(簡稱陷部)。為了讓邊界平滑，最好從水體主邊界上去除細碎的突部和陷部。在圖6中，水體0是外邊界包圍部分1減去陷部對應(yīng)的空洞2，而1由水體主體3和突部4構(gòu)成。簡單地通過面積閾值即可判斷是否該保留陷部和突部。

圖6 邊界修剪示意圖

方法如下。

步驟1：對待修剪邊界計算外邊界和空洞邊界。

步驟2：如果空洞面積小于Athre，則拋棄空洞。

步驟3：對外邊界利用最小角準則進行連接，得到1個或多個不重疊的多邊形，稱為子多邊形。

步驟4：如果子多邊形面積小于Athre，則拋棄子多邊形，否則保留。

步驟5：對剩余空洞和子多邊形重構(gòu)拓撲關(guān)系。

3 實驗

本文利用2019年9月由RIEGL VQ-1560i機載LiDAR采集的江蘇省常州市小黃山地區(qū)的數(shù)據(jù)進行實驗。小黃山的地類以森林、農(nóng)田、居民區(qū)、道路和水體為主。在測區(qū)內(nèi)分布著200多個水體，絕大部分為靜態(tài)水體，包括天然湖、魚塘、沼澤等，最大面積約62 000 m2，最小面積約120 m2。數(shù)據(jù)采集時設(shè)計地面點密度為16點/m2。除點云外，還同時獲取了可見光影像，生成了0.1 m分辨率的DOM。在DOM輔助下，選取了220份水體數(shù)據(jù)，每份數(shù)據(jù)包括一個或多個鄰近水體。采用C++編程實現(xiàn)了本文方法。取Lthre=1.0 m，Athre=1.0 m2。

3.1 LiDAR測量誤差分析

在DOM輔助下，選取了5塊位于局部近似平面的點云，去噪，計算高程標準差作為點云測量的高程相對誤差σmeas，結(jié)果見表1。可以看出，σmeas分布在[0.017，0.025]內(nèi)。因此，σwater需要大于等于0.025 m。

表1 局部近似平面的高程標準差

3.2 水面高程范圍估計

圖7是典型的水體邊界點高程分布直方圖。在11.47、11.52和11.54 m處存在3個離群低點。而從11.56 m開始到11.98 m，每個厘米值上均有點分布。設(shè)ds=0.01 m，則高程在[11.47，11.98]內(nèi)的點共分為4個高程類。取Nzthre=6，則Zmin=11.56 m，Zmax=11.98 m。

圖7 典型的水體邊界點高程分布直方圖

圖8展示了Zwmean和σwater對Zwmax的依賴關(guān)系?？梢钥闯觯琙wmean和σwater是Zwmax的單調(diào)遞增函數(shù)。σwater越大，對應(yīng)的Zwmax越大。

圖8 Zwmean和σwater隨Zwmax變化的折線圖

3.3 水面高程標準差對邊界提取的影響

σwater取值影響水面高程范圍估計，進而影響最終結(jié)果。為此，實驗了0.026、0.028、0.030、0.032和0.034 m共5個σwater值用于水體邊界提取的效果，取Eiso=12，Niso=2?；诓煌襴ater得到的水體邊界一般存在差異。對結(jié)果準確度進行定量評估非常困難。但是，通過與測區(qū)DOM疊加顯示，可以目視判斷由點云提取的水體邊界和影像中邊界的套合程度，作為評估結(jié)果的依據(jù)。

結(jié)果表明，對所有220個數(shù)據(jù)來說，這5個σwater值均100%提取到主要水體。在水體內(nèi)沒有很多點的情況下，即使σwater較小也能取得良好的效果，使用較大的值并不會明顯改進效果。

圖9的魚塘中，σwater取0.026 m和0.034 m得到的邊界基本一致，只有細微差別。

但是，對于水體內(nèi)存在大量點的情況，σwater可能明顯影響邊界提取的效果。

圖10中的4個完整魚塘，除了中間的小魚塘外，其他三個大魚塘在中南部均存在大量點，其中極少量點位于增氧機，其余點位于水面。這可能是在采集數(shù)據(jù)時魚塘南部刮南風(fēng)造成的，未影響到魚塘北部。在σwater=0.026 m時，三個大魚塘內(nèi)部都存在大量空洞。隨著σwater的增大，空洞越來越少。到σwater=0.034 m時，除了魚塘內(nèi)的增氧機外，已經(jīng)沒有任何水面點造成的空洞，而魚塘外邊界、增氧機邊界與DOM中的邊界套合良好。

注：在(b)、(c)內(nèi)，黑色實線為提取到的邊界。圖9 水體內(nèi)點較少情況時的邊界提取效果示例

注：在(b)至(f)內(nèi)，黑色實線為提取到的邊界。圖10 水體內(nèi)點較多情況時的邊界提取效果示例

4 結(jié)束語

水面高程標準差σwater的取值對水體內(nèi)點少的情況影響小，但對水體內(nèi)點較多的情況影響大，較大的σwater傾向于給出較好的效果。對后者來說，增大σwater會增大Zwmax，使更多波浪點的高程小于等于Zwmax，從而在基于高程的邊界擴張中被擴張到。對一些動態(tài)水體，比如噴泉，采用較大的σwater也可以取得良好效果。但是，對于建筑物和水體相鄰的情況，如果建筑物空洞與水體空洞相連，則提取到的水體邊界可能不準。本文方法不適用于沒有足夠大面積的空洞存在，以及挺水植物較多的水體。

本文方法提取的水體邊界可以輔助測繪內(nèi)業(yè)中DEM水體斷裂線的生產(chǎn)，一般經(jīng)少量人工編輯即可獲得較常規(guī)生產(chǎn)手段更真實準確的水體邊界，大大減少人工工作量。

水體邊界提取是空間不確定性研究中的典型問題[16]。由于水體自身特點和觀測手段的限制，難以精確獲取某時刻水體的邊界。如何定量評估水體邊界的準確度需要進一步研究。