面向全球DEM生產(chǎn)的點(diǎn)云智能濾波與DEM泊松編輯方法

胡 翰，丁雨淋，朱 慶，蔣 捷，文學(xué)虎，張 力，唐 偉，陽 俊，鐘若飛

1. 香港理工大學(xué)土地測(cè)量及地理資訊學(xué)系，香港； 2. 國(guó)土資源部城市土地資源監(jiān)測(cè)與仿真重點(diǎn)實(shí)驗(yàn)室，廣東 深圳 518000； 3. 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756； 4. 香港中文大學(xué)太空與地球信息科學(xué)研究所，香港； 5. 首都師范大學(xué)成像技術(shù)高精尖創(chuàng)新中心，北京 100037； 6. 國(guó)家基礎(chǔ)地理信息中心，北京 100083； 7. 四川省第二測(cè)繪地理信息工程院，四川 成都 610100； 8. 中國(guó)測(cè)繪科學(xué)研究院，北京 100083； 9. 黑龍江地理信息工程院,黑龍江 哈爾濱 150086； 10. 首都師范大學(xué)資源環(huán)境與旅游學(xué)院，北京 100037

“一帶一路”倡議和“走出去”國(guó)家戰(zhàn)略急需全球地理信息資源保障，由于全球地理信息資源的缺乏，我國(guó)急需開展“一帶一路”重點(diǎn)區(qū)域地理信息資源建設(shè)[1]。為實(shí)現(xiàn)全球重點(diǎn)區(qū)域的高精度數(shù)字高程模型(DEM)覆蓋，如何利用國(guó)產(chǎn)高分辨率衛(wèi)星影像高效生產(chǎn)全球DEM已經(jīng)成為我國(guó)全球地理信息資源建設(shè)工程的重大任務(wù)。

全球尺度的三維地理信息獲取通常采用光學(xué)立體影像或雷達(dá)干涉測(cè)量的方式，例如我國(guó)的資源三號(hào)(ZY-3)三線陣衛(wèi)星[2-3]與德國(guó)TanDEM-X雷達(dá)衛(wèi)星[4]。歐洲空中客車防務(wù)及航天公司(Airbus Defence and Space)，利用TanDEM-X雷達(dá)衛(wèi)星在2010年底至2014年采集的高分雷達(dá)生成DSM數(shù)據(jù)，通過對(duì)地形及水域等典型區(qū)域的人工交互式編輯，生產(chǎn)了覆蓋全球的高精度WorldDEM產(chǎn)品[5](格網(wǎng)分辨率12 m，垂直精度為2 m(相對(duì))/4 m(絕對(duì)))。我國(guó)國(guó)家測(cè)繪地理信息局于2015年啟動(dòng)全球地理信息資源建設(shè)工程，已利用資源三號(hào)衛(wèi)星(ZY-3)影像數(shù)據(jù)生產(chǎn)了境內(nèi)外DOM和DSM產(chǎn)品[6]。針對(duì)我國(guó)全球DEM生產(chǎn)任務(wù)，擬利用ZY-3影像的DSM產(chǎn)品，通過濾波編輯進(jìn)行高精度DEM生產(chǎn)。然而，我國(guó)各生產(chǎn)單位采用的傳統(tǒng)立體測(cè)圖的方式進(jìn)行DEM生產(chǎn)，一幅標(biāo)準(zhǔn)的1∶5萬DEM產(chǎn)品，耗時(shí)需1人/月，其中立體視圖下采集高程點(diǎn)需1～3 d，剩余90%的工作量則主要耗費(fèi)在植被、水域、城區(qū)等特殊地形地貌特征的交互式修復(fù)、修補(bǔ)等后處理工作。DEM生產(chǎn)的后處理操作量大，費(fèi)時(shí)費(fèi)力，且生產(chǎn)效率及產(chǎn)品精度完全依靠作業(yè)員的經(jīng)驗(yàn)和熟練程度，效率和精度難以保證。技術(shù)上仍面臨以下挑戰(zhàn)：

(1) 傳統(tǒng)點(diǎn)云濾波方法的瓶頸問題。傳統(tǒng)點(diǎn)云濾波算法原理上需假設(shè)在局部鄰域內(nèi)的最低點(diǎn)為地面點(diǎn)，通過建立臨時(shí)的地面結(jié)構(gòu)，利用閾值參數(shù)對(duì)地面點(diǎn)進(jìn)行加密，區(qū)分地面與非地面點(diǎn)。全球地形起伏變化、地表結(jié)構(gòu)形態(tài)復(fù)雜(如圖1所示)，在點(diǎn)云密度和數(shù)據(jù)特性上均有較大差異?，F(xiàn)有濾波算法通常對(duì)同一區(qū)域采用一套參數(shù)，難以適應(yīng)錯(cuò)綜復(fù)雜的地面-非地面地形結(jié)構(gòu)，且對(duì)參數(shù)變化較為敏感，導(dǎo)致誤分類問題嚴(yán)重[7-8]。當(dāng)區(qū)域內(nèi)地形地貌種類較多時(shí)，無法批量處理整塊區(qū)域、濾波參數(shù)顧此失彼，需大量的交互式后處理操作。例如在處理圖1(c)所示的密集城區(qū)時(shí)，常出現(xiàn)大量房屋無法濾除或因過度濾波導(dǎo)致的山地地形特征缺失或河岸線消退等錯(cuò)誤結(jié)果。自動(dòng)、高效、可靠的點(diǎn)云濾波算法依舊是學(xué)術(shù)界尚未解決的難題，更是生產(chǎn)實(shí)踐普遍面臨的技術(shù)瓶頸[9-10]。

(2) 利用ZY-3衛(wèi)星影像的DSM數(shù)據(jù)進(jìn)行DEM生產(chǎn)的難題。光學(xué)衛(wèi)星影像的DSM，存在明顯噪聲及缺漏現(xiàn)象。采用光學(xué)立體衛(wèi)星影像，通過影像密集方式獲取密集匹配點(diǎn)云，并進(jìn)行DSM重建的過程，主要步驟包括：利用影像灰度信息，獲取影像像素點(diǎn)的坐標(biāo)對(duì)應(yīng)關(guān)系，并通過空間前方交會(huì)方式獲取地面/地表點(diǎn)的三維坐標(biāo)。不同于激光雷達(dá)(LiDAR)點(diǎn)云，影像密集匹配點(diǎn)云受限于影像紋理信息，在無紋理特征、紋理特征較弱、紋理特征重復(fù)出現(xiàn)的區(qū)域，由于像素灰度信息難以可靠恢復(fù)影像對(duì)應(yīng)關(guān)系，難以獲取有效的同名點(diǎn)，導(dǎo)致匹配噪聲多甚至匹配失敗。在利用DSM數(shù)據(jù)進(jìn)行DEM生產(chǎn)的過程中，此類噪聲以及缺失的地面信息，均需要后續(xù)交互式的平滑、修補(bǔ)處理等后處理操作。例如，無紋理特征的水域(如圖2所示)，傳統(tǒng)立體測(cè)圖方式下，DSM通常需進(jìn)行人工交互式編輯，在DSM自動(dòng)濾波或交互式編輯中，因水域高程低于河岸，常會(huì)出現(xiàn)河岸擴(kuò)張、湖中島缺失現(xiàn)象。因此，構(gòu)建顧及點(diǎn)云噪聲影響及保地形特征的點(diǎn)云濾波方法，是目前亟待解決的問題。

圖1 不同地形地表的DSM數(shù)據(jù)Fig.1 DSM data for different scenarios

圖2 水域處理效果Fig.2 Processing results for the water bodies

不同于LiDAR點(diǎn)云，光學(xué)影像密集匹配點(diǎn)云無多次回波效應(yīng)，無法穿透植被，因此在大范圍植被覆蓋區(qū)域，密集匹配點(diǎn)云中不包含地面點(diǎn)信息[11]。傳統(tǒng)DEM生產(chǎn)過程，對(duì)于大范圍林地覆蓋的區(qū)域，多是采用人工交互方式，將DSM降低固定的植被高度值后再做整體平滑操作，可靠性無法保證[12]。因此對(duì)于ZY-3 DSM，在林地等不包含地面信息的區(qū)域，需引入全球土地利用信息，如GLC30、正射影像分類信息，利用植被覆蓋信息，進(jìn)行定向精細(xì)濾波，降低植被高度。并且，由于林地通常出現(xiàn)于山區(qū)，在定向?yàn)V波的同時(shí)，需同時(shí)保留山地的典型地形特征，如山脊、山谷等特征線。

最后，DSM在高樓密集區(qū)域會(huì)因?yàn)檎趽鹾鸵暡顢嗔言驅(qū)е碌孛纥c(diǎn)缺失問題，導(dǎo)致在建筑物邊緣出現(xiàn)缺漏，使建筑物的直角特征缺失、變平滑。這類地形區(qū)域的自動(dòng)濾波結(jié)果通常呈現(xiàn)為某些或成片建筑無法過濾完整，或建筑物區(qū)域存在顯著高于周圍地表、卻低于建筑物高度的離散格網(wǎng)點(diǎn)。因此針對(duì)上述現(xiàn)象，需要針對(duì)性的能獲取平滑地表且可降低建筑物高度的后處理交互式編輯方法。

綜上所述，我國(guó)擬用ZY-3衛(wèi)星影像的DSM產(chǎn)品進(jìn)行全球DEM生產(chǎn)。ZY-3影像的DSM產(chǎn)品具有如下特點(diǎn)：①在低紋理、水面等區(qū)域存在明顯噪聲、缺漏現(xiàn)象；②在成片林地或其他大面積植被覆蓋區(qū)域，不包含地面信息，因此違背了濾波算法的局部鄰域內(nèi)最低高程為地面高程的基本假設(shè)，需要引入全球土地利用信息進(jìn)行定向?yàn)V波；③在建筑物覆蓋密集的城市區(qū)域，會(huì)由于遮擋、視差斷裂等因素，導(dǎo)致在建筑物邊緣出現(xiàn)缺漏，使建筑物的直角特征缺失、變平滑，該特性違背了點(diǎn)云濾波的地面平滑的基本假設(shè)，需要針對(duì)性的后處理交互式編輯。針對(duì)上述問題，筆者提出了一種面向全球DEM生產(chǎn)的點(diǎn)云智能濾波與泊松編輯方法，①顧及彎曲能量的點(diǎn)云自適應(yīng)濾波，以彎曲能量系統(tǒng)顯式定量刻畫復(fù)雜的地形地表結(jié)構(gòu)特征，并以此驅(qū)動(dòng)濾波模型參數(shù)的自適應(yīng)優(yōu)選；②構(gòu)建多邊界約束的地形泊松編輯框架，引入分類、地形地貌特征線等邊界信息，支撐并實(shí)現(xiàn)城區(qū)、林地、山地、水域等復(fù)雜地形的交互式保特征泊松編輯。

1 面向全球DEM生產(chǎn)的點(diǎn)云智能濾波與泊松編輯原理

1.1 顧及彎曲能量的點(diǎn)云自適應(yīng)濾波方法

顧及彎曲能量的點(diǎn)云自適應(yīng)濾波方法(圖3)，主要是針對(duì)DSM點(diǎn)云數(shù)據(jù)，在局部鄰域的最低點(diǎn)為地面點(diǎn)與地表平滑這兩個(gè)前提下，采用金字塔的濾波策略，逐級(jí)加密非地面點(diǎn)，從而逼近真實(shí)地表。在金字塔濾波的過程中，通過帶有抗噪性顧及彎曲能量正則化約束的DEM內(nèi)插算法與自適應(yīng)的參數(shù)優(yōu)選算法，提高自動(dòng)濾波算法對(duì)噪聲與參數(shù)選擇的穩(wěn)健性。

1.1.1 漸進(jìn)金字塔濾波策略

基于地面最低高程假設(shè)，在構(gòu)建點(diǎn)云金字塔時(shí)，通常采用局部鄰域內(nèi)的最低點(diǎn)作為金字塔對(duì)應(yīng)分辨率的格網(wǎng)數(shù)值，因此在從粗到精的金字塔濾波策略中，每層的點(diǎn)云屬于地面點(diǎn)的可能性將逐漸變小。若采用傳統(tǒng)的四叉樹的金字塔結(jié)構(gòu)，如圖4(a)所示，每層級(jí)之間的點(diǎn)云數(shù)量將為1∶4，待分類的點(diǎn)過多，易導(dǎo)致地面點(diǎn)誤分，而錯(cuò)誤點(diǎn)將在逐級(jí)濾波過程中被不斷放大。因此，本方法設(shè)計(jì)提出了漸進(jìn)金字塔的構(gòu)建方法，即金字塔的每層級(jí)之間分辨率漸進(jìn)加密，增加層級(jí)數(shù)目，提高分類可靠性[13]，如圖4(b)所示。

圖4 漸進(jìn)金字塔構(gòu)建Fig.4 Construction of the progressive pyramid

1.1.2 顧及彎曲能量正則化約束的局部抗噪性內(nèi)插方法

傳統(tǒng)的DEM內(nèi)插方法，通常是通過擬合地面控制點(diǎn)，P={pi=(xi,yi,zi)|i=1,2,3,…,n}，得到一張參數(shù)化的曲面，z=f(x,y)，在擬合過程中，由于地面控制點(diǎn)中不包含粗差，因此僅需要保持對(duì)地面控制點(diǎn)即可有較好的擬合度，即數(shù)據(jù)擬合度εdata(data term)具有較小的值，或完全經(jīng)過所有控制點(diǎn)，如式(1)所示

(1)

在實(shí)際點(diǎn)云濾波過程中，地面點(diǎn)不可避免地會(huì)引入非地面粗差點(diǎn)，因此此時(shí)即使εdata為0，也并不意味著對(duì)地面有較好的擬合程度?，F(xiàn)有的研究已經(jīng)證明，采用地表一致性約束，可以一定程度上，降低噪聲敏感性，提高內(nèi)插算法的抗差性?，F(xiàn)有的聚類、平面擬合的一致性約束方法并不通用，且計(jì)算耗時(shí)。彎曲能量εsmooth見式(2)

(2)

式(2)可以表達(dá)一張參數(shù)化曲面的二階連續(xù)程度，因此若將εsmooth作為一個(gè)正則化約束引入地面內(nèi)插過程中，則可隱式地蘊(yùn)含著地表一致性約束[14]，抵抗噪聲點(diǎn)的干擾。即通過優(yōu)化如下能量函數(shù)，進(jìn)行DEM內(nèi)插，如式(3)所示

ε=εdata+λεsmooth

(3)

該能量函數(shù)可以很方便地采用帶正則化約束的薄板樣條函數(shù)來擬合求解[15]。

1.1.3 顧及彎曲能量的濾波參數(shù)自適應(yīng)優(yōu)選方法

在DEM內(nèi)插過程中，每個(gè)格網(wǎng)點(diǎn)對(duì)應(yīng)于一個(gè)參數(shù)化的薄板樣條函數(shù)曲面，因此也可以計(jì)算其彎曲能量值?？紤]到在山脊、地面斷裂處，在由粗到精的點(diǎn)云濾波過程中，內(nèi)插得到的地表通常會(huì)低于真實(shí)地表值，如圖5所示。一方面，由于這些區(qū)域的地形結(jié)構(gòu)彎曲起伏更為劇烈，通常表現(xiàn)為較大的彎曲能量，因此在對(duì)這些區(qū)域的點(diǎn)云進(jìn)行濾波處理時(shí)，依據(jù)彎曲能量的大小補(bǔ)償濾波參數(shù)閾值，以此顯式約束濾波參數(shù)閾值自適應(yīng)調(diào)整。另一方面，在金字塔濾波策略中，由于窗口大小不同，因此為了顧及坡度變化的影響，也需要考慮到尺度變化的影響，濾波閾值應(yīng)隨尺度的增大而增加。因此最終的濾波閾值Zt，如式(4)所示

zt=f(t,s,b)

(4)

式中，t為人為指定的初始閾值；s為尺度信息；b為彎曲能量信息。通過以上信息實(shí)現(xiàn)濾波閾值參數(shù)隨地形變化、尺度的自適應(yīng)優(yōu)選。

圖5 顧及彎曲能量的濾波參數(shù)自適應(yīng)優(yōu)選方法Fig.5 Automatic parameter tuning by bending energies

本文方法是在筆者的已有研究[13]基礎(chǔ)上進(jìn)行了擴(kuò)展。前期研究主要是針對(duì)Lidar點(diǎn)云濾波瓶頸，提出了顧及彎曲能量的Lidar點(diǎn)云自適應(yīng)濾波方法，該方法利用ISPRS提供的17種標(biāo)準(zhǔn)數(shù)據(jù)集進(jìn)行了多組試驗(yàn)分析，并設(shè)計(jì)了誤差評(píng)價(jià)體系，定量分析了彎曲能量方法的可靠性和自適應(yīng)能力。

1.2 多邊界約束的泊松編輯方法

1.2.1 泊松編輯方法的基本原理

DSM數(shù)據(jù)中可能由于遮擋、大范圍植被、水域等因素導(dǎo)致自動(dòng)濾波算法效率降低，交互式的后處理DEM編輯必不可少。在DEM編輯中需要同時(shí)顧及DEM的平滑特性與水域、山脊特征線等特征信息，傳統(tǒng)簡(jiǎn)單的DEM編輯方法并不能滿足上述需求，因此本方法引入一種泊松地形編輯(Poisson terrain editing)的策略[16]。泊松編輯的含義是，通過優(yōu)化域中的一個(gè)已知導(dǎo)向信息v，并給定邊界條件f|?Ω，通過下列積分方程，求解未知函數(shù)f。定義影像范圍S∈R2，待優(yōu)化域Ω∈S，影像上的一個(gè)閉合區(qū)域其邊界為?Ω，則f為優(yōu)化域Ω的一個(gè)未知函數(shù)。對(duì)本方法而言，函數(shù)f為DEM格網(wǎng)點(diǎn)的高程值

(5)

對(duì)于DEM編輯而言，優(yōu)化域即DEM上的一個(gè)可帶洞的閉合多邊形區(qū)域，在區(qū)域內(nèi)，其導(dǎo)向域條件是依不同濾波算法給定的，該優(yōu)化域必須閉合，即被邊界條件完全包含。由于優(yōu)化域可為帶洞多邊形，因此在優(yōu)化域內(nèi)部也可以包含任意邊界條件，如圖6所示。一方面，由于邊界條件在優(yōu)化過程中保持不變，因此可以用于保留DEM中的地形結(jié)構(gòu)特征信息，如山脊線、山谷線、河岸線、河流區(qū)域等。另一方面，通過調(diào)整優(yōu)化域Ω中的導(dǎo)向信息的定義方式，實(shí)現(xiàn)不同類型的濾波效果，如平滑、恢復(fù)梯度信息等，如圖7所示。

圖6 泊松地形編輯的導(dǎo)向域與邊界條件Fig.6 The guided field and boundary conditions for Poisson terrain editing

1.2.2 泊松地形編輯方法

本方法中，泊松地形編輯方法包含3個(gè)主要步驟：①依據(jù)濾波需求創(chuàng)建多邊形區(qū)域，定義優(yōu)化域與包圍優(yōu)化域的邊界條件；②在優(yōu)化域中定義合適的導(dǎo)向信息以實(shí)現(xiàn)不同的濾波效果，并給定合適的邊界條件的初值；③求解方程，獲取泊松方程的解，即濾波后的DEM。

圖7 DEM編輯中的邊界條件Fig.7 Boundary conditions in DEM editing process

1.2.2.1 交互式編輯創(chuàng)建多邊形

如圖6和圖7所示，在二維影像域S中，待編輯的優(yōu)化域Ω及其邊界區(qū)域?Ω，都可用多邊形進(jìn)行表達(dá)，其中?Ω可為Ω的邊界線或Ω內(nèi)部的空洞，此在本方法中采用通用的帶洞多邊形來支持優(yōu)化域與邊界條件的創(chuàng)建與編輯。本方法主要支持兩種方式的多邊形創(chuàng)建方法：①交互式勾繪創(chuàng)建多邊形。②采用區(qū)域增長(zhǎng)的方式創(chuàng)建多邊形區(qū)域。對(duì)于前者，可以采用折線或流形的方式操作，其中折線為每次點(diǎn)擊創(chuàng)建一個(gè)點(diǎn)，最后閉合多邊形，流形通過拖動(dòng)鼠標(biāo)繪制平滑曲線并閉合。而區(qū)域增長(zhǎng)可根據(jù)DSM/DEM高程值或分類信息進(jìn)行區(qū)域增長(zhǎng)，并通過等值線提取算法獲取閉合多邊形區(qū)域。

1.2.2.2 濾波算法創(chuàng)建

在創(chuàng)建了優(yōu)化域Ω與邊界區(qū)域?Ω后，泊松方程的設(shè)計(jì)矩陣和結(jié)構(gòu)實(shí)質(zhì)上已經(jīng)可以確定，但為了獲取最終優(yōu)化域Ω中的未知函數(shù)f，仍然需要?jiǎng)?chuàng)建優(yōu)化域中的導(dǎo)向信息v與邊界區(qū)域的初值f*|?Ω，此信息決定了不同的濾波效果。目前本方法已集成實(shí)現(xiàn)了十幾種不同的濾波效果，在此列舉幾個(gè)說明導(dǎo)向信息與初值的創(chuàng)建方式。

(1) 復(fù)制DSM濾波：此濾波為恢復(fù)原始DSM的形狀，用于在某些濾波過程中由于濾波選擇不恰當(dāng)導(dǎo)致編輯DEM失真，此濾波可以保證在DEM邊界高程不變的情況下，恢復(fù)DSM的地形起伏特征。值得注意的是，此算法并非直接恢復(fù)DSM的高程值，因?yàn)檫@樣會(huì)導(dǎo)致編輯的邊緣產(chǎn)生突變。對(duì)此算法，導(dǎo)向域是直接從DSM中采用拉普拉斯算子計(jì)算所得的，而初值為DEM邊界區(qū)域中的高程值。

(2) 林地高程降低濾波：在林地，可通過分類信息創(chuàng)建優(yōu)化域與邊界區(qū)域，并給定一個(gè)待降低的高程值，整體降低林地高度。然而若直接整體降低高程，勢(shì)必會(huì)在邊緣區(qū)域產(chǎn)生突變效應(yīng)，因此對(duì)于此濾波，導(dǎo)向域?yàn)楫?dāng)前DEM的拉普拉斯算子計(jì)算所得，即保證當(dāng)前DEM的起伏特性不變，邊界區(qū)域的初值，在編輯區(qū)域的邊緣采用DEM的高程值，而在內(nèi)部隨機(jī)創(chuàng)建一系列點(diǎn)狀的邊界區(qū)域，其初值為DEM減去降低高程值。采用此方法，可保證編輯后植被高度正常降低且不會(huì)產(chǎn)生顯著突變效果。

(3) 均值/中值/最低值平滑：同樣，若簡(jiǎn)單地采用高斯平滑，難免會(huì)在編輯的邊緣處產(chǎn)生突變效果，因此本軟件采用泊松編輯的方式實(shí)現(xiàn)平滑，其核心思路是采用0值作為導(dǎo)向域信息，而采用均值/中值/最低值作為內(nèi)部邊界區(qū)域的初值。該方法等價(jià)于薄膜樣條函數(shù)內(nèi)插(membrane interpolation)[17]。

1.2.2.3 泊松方程求解

在構(gòu)建了泊松方程后，在離散的二維影像域中可將式(5)確定的初值積分方程，轉(zhuǎn)化為一個(gè)大規(guī)模的稀疏方程AX=B，其中A為可逆方陣，而X為最終優(yōu)化域中每個(gè)格網(wǎng)點(diǎn)的值，該方程可采用稀疏矩陣的LU分解[18]或共軛梯度優(yōu)化方法進(jìn)行求解(conjugate gradient)[16,19]。由于傳統(tǒng)基于CPU的優(yōu)化算法時(shí)間較長(zhǎng)，難以滿足交互式的實(shí)時(shí)需求，因此本文基于通用的并行GPU加速計(jì)算框架CUDA與OpenCL實(shí)現(xiàn)了共軛梯度法的并行算法，速率較傳統(tǒng)CPU算法提高10倍，可滿足交互式的實(shí)時(shí)需求。

2 試驗(yàn)分析

基于本文所述兩項(xiàng)關(guān)鍵技術(shù)，自主研制了點(diǎn)云智能濾波與泊松編輯軟件LINK。在試驗(yàn)分析部分，本文將著重介紹利用不同類型不同地區(qū)的全球DEM生產(chǎn)性試驗(yàn)情況，并將試驗(yàn)結(jié)果與當(dāng)前我國(guó)其他兩款主流軟件GEOWAY CIPS和PixelGrid進(jìn)行對(duì)比分析。

2.1 數(shù)據(jù)

試驗(yàn)數(shù)據(jù)包括全球地理信息資源建設(shè)工程前期生產(chǎn)的DSM、DOM成果，均采用WGS-84坐標(biāo)系，以1∶5萬標(biāo)準(zhǔn)分幅為存儲(chǔ)單元?；贒SM生產(chǎn)DEM的難點(diǎn)主要在于大面積植被覆蓋區(qū)域和建筑物覆蓋區(qū)域，因此試驗(yàn)數(shù)據(jù)涵蓋了建筑區(qū)、森林、水域等典型數(shù)據(jù)具體見表1。

表1 試驗(yàn)數(shù)據(jù)描述

試驗(yàn)數(shù)據(jù)分布情況如圖8所示。

圖8 試驗(yàn)數(shù)據(jù)分布Fig.8 The distribution of test data

2.2 試驗(yàn)結(jié)果質(zhì)量定量評(píng)定與分析

本次試驗(yàn)進(jìn)行了多種檢查方法對(duì)DEM成果數(shù)據(jù)進(jìn)行全面地比較、質(zhì)量評(píng)定與分析，其中包括：①采用人工采集的檢查點(diǎn)進(jìn)行質(zhì)量評(píng)定；②利用原始DSM與DEM成果數(shù)據(jù)的差值檢查濾波可靠性；③在立體環(huán)境下檢查DEM成果與空三加密成果符合程度；④以傳統(tǒng)的立體采集地貌方式生產(chǎn)的DEM數(shù)據(jù)為參考檢查濾波成果。以下重點(diǎn)介紹采用人工采集檢查點(diǎn)進(jìn)行質(zhì)量評(píng)定與分析的結(jié)果。

(1) 采用人工采集的檢查點(diǎn)進(jìn)行質(zhì)量評(píng)定。人工采集的檢查點(diǎn)是在DOM上選取位于建筑物頂部、植被頂部等區(qū)域的點(diǎn)位，其X、Y值從DOM獲取，Z值在DSM的高程值基礎(chǔ)上，立體環(huán)境下讀取建筑物、植被的高度，對(duì)DSM高程值進(jìn)行修正，從而得到最終的Z值。在建筑和植被區(qū)域，均勻采集了大量檢查點(diǎn)；在裸地、水域等區(qū)域，主要采集在能反映地形特征位置的檢查點(diǎn)，具體如圖9所示。

圖9 檢查點(diǎn)分布Fig.9 Check points distributions

(2) 精度檢查。利用檢查點(diǎn)數(shù)據(jù)對(duì)3種軟件的濾波結(jié)果進(jìn)行精度檢查，結(jié)果見表2。

表2 精度分析

從表2可以看出，針對(duì)建筑、植被及裸地區(qū)域，本方法的濾波精度明顯優(yōu)于其他兩款軟件，對(duì)于植被區(qū)域的濾波，泊松濾波方法優(yōu)勢(shì)更突出。

2.3 與傳統(tǒng)人工立體采集生產(chǎn)的DEM對(duì)比

為了將DSM濾波生產(chǎn)的DEM與傳統(tǒng)方式生產(chǎn)的DEM進(jìn)行效果比較，本次試驗(yàn)采用立體采集地貌的方式生產(chǎn)了6幅試驗(yàn)數(shù)據(jù)的DEM成果，并對(duì)兩種生產(chǎn)方式的DEM效果、精度進(jìn)行了詳細(xì)的對(duì)比分析，如圖10所示。

從表3可以看出，利用DSM濾波生產(chǎn)DEM的方式相比傳統(tǒng)生產(chǎn)方式而言，精度明顯提高。證明本文方法在進(jìn)行DEM生產(chǎn)時(shí)，避免了人工采集誤差、格網(wǎng)點(diǎn)內(nèi)插誤差等影響因素，可以盡量保持實(shí)際地貌細(xì)節(jié)。本方法相比傳統(tǒng)生產(chǎn)方式而言，生產(chǎn)效率大大提高，約提高至少2～3倍。尤其是在山地、高山地等DEM生產(chǎn)困難地區(qū)，本方法的優(yōu)勢(shì)更加明顯。

圖10 成果DEM對(duì)比Fig.10 Comparison between different DEM results

圖幅地形檢查點(diǎn)個(gè)數(shù)中誤差/m傳統(tǒng)方式DEM成果LINK生產(chǎn)DEM成果NE47E006011平地981.8711.573NE47E019001平地3644.5743.430NE47E017002丘陵地2905.1052.652NH41E005013丘陵地4234.2651.618NE47E015008山地2586.9842.646NJ42E020023高山地435.1414.007

2.4 與其他主流軟件的試驗(yàn)結(jié)果對(duì)比分析

在對(duì)比分析LINK與當(dāng)前我國(guó)其他兩款主流軟件GEOWAY CIPS和PixelGrid的濾波結(jié)果時(shí)，將重點(diǎn)分析建筑、植被和水域的處理效果：建筑、植被是否有效濾掉，是否將非地面點(diǎn)有效降至地面；在水域等地貌突變區(qū)域，DEM數(shù)據(jù)是否準(zhǔn)確保留了地貌細(xì)節(jié)、地形結(jié)構(gòu)特征，比如湖中小島，河流邊緣紋理，細(xì)窄的溝壑，紋理變化復(fù)雜的山地區(qū)域等。

2.4.1 建筑區(qū)濾波效果對(duì)比

建筑區(qū)濾波效果對(duì)比，如圖11所示，可見其他兩款軟件均存在大量建筑物未濾除，LINK對(duì)建筑物區(qū)域的濾波效果較好。LINK軟件設(shè)計(jì)了建筑物區(qū)域的泊松濾波模塊，可根據(jù)自動(dòng)分類提取出的建筑區(qū)域多邊形，實(shí)現(xiàn)建筑區(qū)域的局部自適應(yīng)濾波；根據(jù)地形特征可將建筑區(qū)分為平地區(qū)域、丘陵區(qū)域及山地高山地區(qū)域，通過調(diào)整鄰域大小及坡度閾值達(dá)到不同地形建筑區(qū)的最佳濾波效果。

圖11 建筑物濾波效果Fig.11 Filtering results for building areas

2.4.2 植被區(qū)域?yàn)V波效果對(duì)比

由于DSM成果的格網(wǎng)間距為10 m，地物地貌細(xì)節(jié)體現(xiàn)不足，常規(guī)統(tǒng)一濾波算法難以將大面積植被區(qū)域?qū)?yīng)的非地面點(diǎn)降至地面。LINK針對(duì)性設(shè)計(jì)了面向植被區(qū)域的泊松濾波模塊，即根據(jù)自動(dòng)分類提取出的植被區(qū)域多邊形及其植被高度屬性值，結(jié)合林地區(qū)域所在的地形特征，通過調(diào)整平滑大小和降低高程兩個(gè)參數(shù)值，在保持地形特征的同時(shí)將植被區(qū)域非地面點(diǎn)高程降至地面。在大面積植被覆蓋區(qū)域，通過濾波處理后，濾波結(jié)果的高程值普遍小于原始DSM高程值。因此試驗(yàn)將原始DSM數(shù)據(jù)減去3種軟件的濾波結(jié)果，對(duì)柵格相減的值進(jìn)行比較，效果如下(白色為DSM-DEM≤0，黑色為DSM-DEM≥0)，如圖12所示。

通過對(duì)比3種軟件DSM與DEM差值柵格圖可以看出，LINK軟件可以根據(jù)圈定的植被范圍，有針對(duì)性地將植被區(qū)域的非地面點(diǎn)降至地面，而CIPS和PixelGrid采用的整體濾波方法，難以識(shí)別大片植被覆蓋區(qū)域，濾波效果不理想。

2.4.3 水域?yàn)V波效果對(duì)比

對(duì)于水域區(qū)域，另兩款軟件均無特定的處理模塊，濾波結(jié)果均需要進(jìn)行大量的交互式后處理操作。LINK軟件針對(duì)性的設(shè)計(jì)了水域泊松濾波模塊。DSM成果數(shù)據(jù)的湖泊等靜止水域范圍是整體置平的，河流等流動(dòng)水域是分段置平的，針對(duì)這一特點(diǎn)，LINK提供了水域自動(dòng)提取功能，基于自動(dòng)提取的水域范圍設(shè)置緩沖區(qū)域，并將水域范圍在后續(xù)處理中會(huì)作為泊松濾波的控制多邊形存在，從而使河岸線等水域特征不受周圍數(shù)據(jù)編輯處理的影響。

圖12 大面積植被覆蓋區(qū)域DSM-DEM結(jié)果Fig.12 Filtering results for vegetation areas

(1) 湖心島濾波效果對(duì)比。湖心島濾波結(jié)果如圖13所示。

圖13 湖心島濾波結(jié)果Fig.13 Filtering results for island

(2) 河岸線濾波效果對(duì)比。用原始DSM數(shù)據(jù)與3種軟件的濾波結(jié)果進(jìn)行柵格減運(yùn)算，通過河岸線區(qū)域的變化值來說明3種軟件的濾波方法對(duì)河岸線的保持效果。由圖14可見，LINK能夠完整保持河岸線特征，而其他兩款軟件會(huì)出現(xiàn)河岸線消退的情況。

3 結(jié) 論

針對(duì)全球地形地表結(jié)構(gòu)具有多樣性和復(fù)雜性特點(diǎn)，現(xiàn)有依靠單一濾波模型或有限濾波規(guī)則的點(diǎn)云濾波方法均難以根據(jù)復(fù)雜多樣的地表結(jié)構(gòu)特征進(jìn)行自適應(yīng)調(diào)整，濾波結(jié)果的準(zhǔn)確性和可靠性問題十分突出，導(dǎo)致大量人工交互后處理工作量，DEM的質(zhì)量和精度難以保證。針對(duì)上述問題，本文提出了一種智能濾波與泊松編輯方法：以彎曲能量定量顯式刻畫地表形狀特征，以此驅(qū)動(dòng)濾波模型的自適應(yīng)調(diào)整，提高濾波結(jié)果的可靠性；設(shè)計(jì)提出了多邊界約束的泊松地形編輯框架，引入分類、地形地貌特征線等邊界信息，實(shí)現(xiàn)建筑區(qū)、森林和水域等區(qū)域的定向智能精準(zhǔn)編輯。經(jīng)過大量不同類型不同地區(qū)的全球DEM生產(chǎn)性試驗(yàn)，證明了本文方法的可靠性和有效性，為全球DEM大規(guī)模生產(chǎn)提供了有力的技術(shù)方法支撐。

圖14 河岸線特征保持情況(DSM-DEM)Fig.14 Filtering results for river bank