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

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

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

全球尺度的三維地理信息獲取通常采用光學(xué)立體影像或雷達(dá)干涉測量的方式，例如我國的資源三號(ZY-3)三線陣衛(wèi)星[2-3]與德國TanDEM-X雷達(dá)衛(wèi)星[4]。歐洲空中客車防務(wù)及航天公司(Airbus Defence and Space)，利用TanDEM-X雷達(dá)衛(wèi)星在2010年底至2014年采集的高分雷達(dá)生成DSM數(shù)據(jù)，通過對地形及水域等典型區(qū)域的人工交互式編輯，生產(chǎn)了覆蓋全球的高精度WorldDEM產(chǎn)品[5](格網(wǎng)分辨率12 m，垂直精度為2 m(相對)/4 m(絕對))。我國國家測繪地理信息局于2015年啟動全球地理信息資源建設(shè)工程，已利用資源三號衛(wèi)星(ZY-3)影像數(shù)據(jù)生產(chǎn)了境內(nèi)外DOM和DSM產(chǎn)品[6]。針對我國全球DEM生產(chǎn)任務(wù)，擬利用ZY-3影像的DSM產(chǎn)品，通過濾波編輯進(jìn)行高精度DEM生產(chǎn)。然而，我國各生產(chǎn)單位采用的傳統(tǒng)立體測圖的方式進(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ù)對地面點(diǎn)進(jìn)行加密，區(qū)分地面與非地面點(diǎn)。全球地形起伏變化、地表結(jié)構(gòu)形態(tài)復(fù)雜(如圖1所示)，在點(diǎn)云密度和數(shù)據(jù)特性上均有較大差異?，F(xiàn)有濾波算法通常對同一區(qū)域采用一套參數(shù)，難以適應(yīng)錯綜復(fù)雜的地面-非地面地形結(jié)構(gòu)，且對參數(shù)變化較為敏感，導(dǎo)致誤分類問題嚴(yán)重[7-8]。當(dāng)區(qū)域內(nèi)地形地貌種類較多時(shí)，無法批量處理整塊區(qū)域、濾波參數(shù)顧此失彼，需大量的交互式后處理操作。例如在處理圖1(c)所示的密集城區(qū)時(shí)，常出現(xiàn)大量房屋無法濾除或因過度濾波導(dǎo)致的山地地形特征缺失或河岸線消退等錯誤結(jié)果。自動、高效、可靠的點(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)對應(yīng)關(guān)系，并通過空間前方交會方式獲取地面/地表點(diǎn)的三維坐標(biāo)。不同于激光雷達(dá)(LiDAR)點(diǎn)云，影像密集匹配點(diǎn)云受限于影像紋理信息，在無紋理特征、紋理特征較弱、紋理特征重復(fù)出現(xiàn)的區(qū)域，由于像素灰度信息難以可靠恢復(fù)影像對應(yīng)關(guān)系，難以獲取有效的同名點(diǎn)，導(dǎo)致匹配噪聲多甚至匹配失敗。在利用DSM數(shù)據(jù)進(jìn)行DEM生產(chǎn)的過程中，此類噪聲以及缺失的地面信息，均需要后續(xù)交互式的平滑、修補(bǔ)處理等后處理操作。例如，無紋理特征的水域(如圖2所示)，傳統(tǒng)立體測圖方式下，DSM通常需進(jìn)行人工交互式編輯，在DSM自動濾波或交互式編輯中，因水域高程低于河岸，常會出現(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)過程，對于大范圍林地覆蓋的區(qū)域，多是采用人工交互方式，將DSM降低固定的植被高度值后再做整體平滑操作，可靠性無法保證[12]。因此對于ZY-3 DSM，在林地等不包含地面信息的區(qū)域，需引入全球土地利用信息，如GLC30、正射影像分類信息，利用植被覆蓋信息，進(jìn)行定向精細(xì)濾波，降低植被高度。并且，由于林地通常出現(xiàn)于山區(qū)，在定向?yàn)V波的同時(shí)，需同時(shí)保留山地的典型地形特征，如山脊、山谷等特征線。

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

綜上所述，我國擬用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ū)域，會由于遮擋、視差斷裂等因素，導(dǎo)致在建筑物邊緣出現(xiàn)缺漏，使建筑物的直角特征缺失、變平滑，該特性違背了點(diǎn)云濾波的地面平滑的基本假設(shè)，需要針對性的后處理交互式編輯。針對上述問題，筆者提出了一種面向全球DEM生產(chǎn)的點(diǎn)云智能濾波與泊松編輯方法，①顧及彎曲能量的點(diǎn)云自適應(yīng)濾波，以彎曲能量系統(tǒng)顯式定量刻畫復(fù)雜的地形地表結(jié)構(gòu)特征，并以此驅(qū)動濾波模型參數(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)，主要是針對DSM點(diǎn)云數(shù)據(jù)，在局部鄰域的最低點(diǎn)為地面點(diǎn)與地表平滑這兩個前提下，采用金字塔的濾波策略，逐級加密非地面點(diǎn)，從而逼近真實(shí)地表。在金字塔濾波的過程中，通過帶有抗噪性顧及彎曲能量正則化約束的DEM內(nèi)插算法與自適應(yīng)的參數(shù)優(yōu)選算法，提高自動濾波算法對噪聲與參數(shù)選擇的穩(wěn)健性。

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

基于地面最低高程假設(shè)，在構(gòu)建點(diǎn)云金字塔時(shí)，通常采用局部鄰域內(nèi)的最低點(diǎn)作為金字塔對應(yīng)分辨率的格網(wǎng)數(shù)值，因此在從粗到精的金字塔濾波策略中，每層的點(diǎn)云屬于地面點(diǎn)的可能性將逐漸變小。若采用傳統(tǒng)的四叉樹的金字塔結(jié)構(gòu)，如圖4(a)所示，每層級之間的點(diǎn)云數(shù)量將為1∶4，待分類的點(diǎn)過多，易導(dǎo)致地面點(diǎn)誤分，而錯誤點(diǎn)將在逐級濾波過程中被不斷放大。因此，本方法設(shè)計(jì)提出了漸進(jìn)金字塔的構(gòu)建方法，即金字塔的每層級之間分辨率漸進(jìn)加密，增加層級數(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)中不包含粗差，因此僅需要保持對地面控制點(diǎn)即可有較好的擬合度，即數(shù)據(jù)擬合度εdata(data term)具有較小的值，或完全經(jīng)過所有控制點(diǎn)，如式(1)所示

(1)

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

(2)

式(2)可以表達(dá)一張參數(shù)化曲面的二階連續(xù)程度，因此若將εsmooth作為一個正則化約束引入地面內(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)插過程中，每個格網(wǎng)點(diǎn)對應(yīng)于一個參數(shù)化的薄板樣條函數(shù)曲面，因此也可以計(jì)算其彎曲能量值?？紤]到在山脊、地面斷裂處，在由粗到精的點(diǎn)云濾波過程中，內(nèi)插得到的地表通常會低于真實(shí)地表值，如圖5所示。一方面，由于這些區(qū)域的地形結(jié)構(gòu)彎曲起伏更為劇烈，通常表現(xiàn)為較大的彎曲能量，因此在對這些區(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ò)展。前期研究主要是針對Lidar點(diǎn)云濾波瓶頸，提出了顧及彎曲能量的Lidar點(diǎn)云自適應(yīng)濾波方法，該方法利用ISPRS提供的17種標(biāo)準(zhǔn)數(shù)據(jù)集進(jìn)行了多組試驗(yàn)分析，并設(shè)計(jì)了誤差評價(jià)體系，定量分析了彎曲能量方法的可靠性和自適應(yīng)能力。

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

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

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

(5)

對于DEM編輯而言，優(yōu)化域即DEM上的一個可帶洞的閉合多邊形區(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個主要步驟：①依據(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ū)域增長的方式創(chuàng)建多邊形區(qū)域。對于前者，可以采用折線或流形的方式操作，其中折線為每次點(diǎn)擊創(chuàng)建一個點(diǎn)，最后閉合多邊形，流形通過拖動鼠標(biāo)繪制平滑曲線并閉合。而區(qū)域增長可根據(jù)DSM/DEM高程值或分類信息進(jìn)行區(qū)域增長，并通過等值線提取算法獲取閉合多邊形區(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，仍然需要創(chuàng)建優(yōu)化域中的導(dǎo)向信息v與邊界區(qū)域的初值f*|?Ω，此信息決定了不同的濾波效果。目前本方法已集成實(shí)現(xiàn)了十幾種不同的濾波效果，在此列舉幾個說明導(dǎo)向信息與初值的創(chuàng)建方式。

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

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

(3) 均值/中值/最低值平滑：同樣，若簡單地采用高斯平滑，難免會在編輯的邊緣處產(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)化為一個大規(guī)模的稀疏方程AX=B，其中A為可逆方陣，而X為最終優(yōu)化域中每個格網(wǎng)點(diǎn)的值，該方程可采用稀疏矩陣的LU分解[18]或共軛梯度優(yōu)化方法進(jìn)行求解(conjugate gradient)[16,19]。由于傳統(tǒng)基于CPU的優(yōu)化算法時(shí)間較長，難以滿足交互式的實(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)前我國其他兩款主流軟件GEOWAY CIPS和PixelGrid進(jìn)行對比分析。

2.1 數(shù)據(jù)

試驗(yàn)數(shù)據(jù)包括全球地理信息資源建設(shè)工程前期生產(chǎn)的DSM、DOM成果，均采用WGS-84坐標(biāo)系，以1∶5萬標(biāo)準(zhǔn)分幅為存儲單元?；贒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ì)量定量評定與分析

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

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

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

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

表2 精度分析

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

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

為了將DSM濾波生產(chǎn)的DEM與傳統(tǒng)方式生產(chǎn)的DEM進(jìn)行效果比較，本次試驗(yàn)采用立體采集地貌的方式生產(chǎn)了6幅試驗(yàn)數(shù)據(jù)的DEM成果，并對兩種生產(chǎn)方式的DEM效果、精度進(jìn)行了詳細(xì)的對比分析，如圖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)勢更加明顯。

圖10 成果DEM對比Fig.10 Comparison between different DEM results

圖幅地形檢查點(diǎn)個數(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é)果對比分析

在對比分析LINK與當(dāng)前我國其他兩款主流軟件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ū)濾波效果對比

建筑區(qū)濾波效果對比，如圖11所示，可見其他兩款軟件均存在大量建筑物未濾除，LINK對建筑物區(qū)域的濾波效果較好。LINK軟件設(shè)計(jì)了建筑物區(qū)域的泊松濾波模塊，可根據(jù)自動分類提取出的建筑區(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波效果對比

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

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

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

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

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

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

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

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

3 結(jié) 論

針對全球地形地表結(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ì)量和精度難以保證。針對上述問題，本文提出了一種智能濾波與泊松編輯方法：以彎曲能量定量顯式刻畫地表形狀特征，以此驅(qū)動濾波模型的自適應(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