復雜地形DSM的地面點識別及DEM提取

喻龍華，王雷光，吳 楠，廖聲熙，張懷清

(1. 中國林業科學研究院資源昆蟲研究所，云南 昆明 650224； 2. 西南林業大學林學院，云南 昆明 650224； 3. 中國林業科學研究院資源信息研究所，北京 100091)

DEM描述地面高程信息，它在工程勘測設計、土地利用、地理信息、林業等領域有著廣泛的應用。DEM精度的好壞直接關系到應用結果的可靠性[1-2]，因此，獲得高質量的DEM有著非常重要的意義。DSM包含地物的輪廓和高程信息[3]，因此可以由DSM通過分割和插值得到DEM。李忠等[4]利用附近最小值方法，從平坦地區的DSM中提取出DEM，通過10個檢查點計算出試驗結果的中誤差為3.51 m；其中10個檢查點中，結果與實際參考值相差小于2 m的點只有3個。Ismail等[5]探究了附近最小值法中，濾波窗口大小和DSM分辨率對DEM提取精度的影響，試驗結果表明：窗口尺寸太小則無法對大型建筑物等非地面點進行有效濾除，而尺寸太大會導致部分地面點被濾除，不容易確定最佳濾波尺寸；DSM的分辨率越高效果越好，當分辨率低于1 m時，隨著分辨率的繼續降低，得到的DEM的精度也隨之快速下降。由于該方法采用窗口中最小值作為判斷條件，當窗口位于山區、丘陵等地形復雜地區時，窗口中較高地面的高程值與最小值之間的差值很容易超過設定的閾值，致使產生錯誤的判斷，因此該方法只能用于較平坦的區域。

區域生長是一種重要的圖像串行區域分割算法，其關鍵因素包括初始種子點的選取、生長準測和終止條件[6-7]，并且直接影響圖像分割效果[8]。在生長算法的基礎上，許多研究者提出了擴展算法。Pohle等[9]把待分割圖像的像素值看作正態分布，先用原始區域生長算法估算出分布參數，再將該參數應用到第二遍生長過程中，從而獲得更好的結果。Zheng等[10]將區域生長算法和各向異性濾波技術結合，提出了一種無需種子點的自動分割算法，克服多數區域生長算法對初始種子點的選取順序和位置敏感的問題。無論是地形平坦地區，還是復雜地區，DSM數據中地面點之間都存在連續性和空間相關性[11]。區域生長算法具有能將連續相似性質的像素集進行合并的特性[12]，可以把連續的地面點結合成區，因此可以應用區域生長算法識別地面點和提取DEM。

本文采用區域生長算法和最大類間方差法實現自動選取種子點，以及生長準則和終止條件的自適應選擇，增強了算法的適用性，去除了地形條件對算法的限制，實現從平坦地區和地形復雜山區的DSM中識別地面點和提取DEM。

1 基于區域生長的DEM提取算法

不管是平坦地區還是地形復雜的山區，在全局上，DSM數據的地面點之間都表現出較好的連續性，建筑物、樹等其他地物與其周圍的地面點之間存在較大的跳躍。雖然在局部存在一定的階躍性，特別是在地形復雜的山區，但相鄰地面點之間還具有一定的連續性。因此，利用區域生長算法能夠把鄰近相似的像素點結合成區域的這一特性，可從DSM中有效識別地面點，如圖1所示。

圖1 算法流程

由圖1可以看出，先采用區域生長算法和最大類間方差法獲得DSM中的地面點，然后以八鄰域插值的方法對非地面點進行插值，得到完整的DEM。由于高架道路的起點是附著于地面，而后抬升，在DSM數據中與地面點之間存在很好的連續性。但高架路的主要部分又高于地面點，具有房屋等建筑的特征，若僅用生長運算1，會導致高架道路被錯誤識別到地面點集G中，最終產生錯誤的DEM。因此，當DSM中存在高架道路時需要調用模塊1；如不存在，則不需要。

1.1 種子點的自動選取

以正方形窗口W不重復地連續遍歷DSM數據，窗口W的大小取DSM中非地面地物的平均尺寸，如50×50。窗口W中滿足式(1)的像素點為種子點，但窗口W可能全部落在某一建筑物上或平坦的地面處。若是前者，則不能在其中取種子點。若是后者，則沒必要有種子點，因其是相對平坦的地面，它與其周圍窗口必然有一個平滑銜接，因此可以被周圍窗口中地面像素點的生長運算過程所識別出來。這兩種情況有一個共同特征，即當前窗口W的標準方差較小，因此算法中要求窗口W的標準方差大于整個DSM數據的標準方差的0.6倍時，才在W中尋找種子點，以避免前一種情況下產生的錯誤。

W(i,j)

(1)

式中，W(i,j)為窗口中位于(i,j)的像素點；min為當前窗口W中的最小值；level公式如下

level=[gt(W)-min]/t

(2)

式中，gt(W)為窗口W中最大類間方差法對應的閾值；t為常數系數，根據原始數據的標準方差取值，標準方差越大，對應取值也越大，取值范圍在10～20之間效果較好。level隨著不同窗口中gt(W)的變化發生變化。當窗口W中數據的值和值域較大時，表示地形起伏大，即當前窗口位于地形復雜的山區或丘陵，此時gt(W)的取值較大，level隨之取到較大值，在保證種子點質量的同時，獲得盡可能多的種子點，以提高算法的識別和插值效果。因此，在DSM數據處理中，地形平坦和地形復雜的部分都能自動找到較合適的閾值level，實現種子點的自適應選取。

1.2 生長準則和終止條件的自適應取值

種子點基礎上的生長過程是在整個DSM數據中完成的，生長準則和終止條件如下

(3)

式中，S(i,j)為種子點或生長過程中新識別出來的地面點；level與式(2)中相同；G為已經識別出的地面點集。種子點位于平坦位置時，與周圍地面點之間的差值相對較小，此時應該取相應較小的level作為生長準則和終止條件的限制閾值。當種子點位于起伏較大的復雜地形時，與周圍地面點之間的差值變大，此時level需要相應取較大的值。式(2)使level隨地形起伏程度自動取值，最終實現平坦和復雜地形下生長準則和終止條件的自適應取值。

1.3 城區高架道路識別

在已識別出來的點集G的基礎上再應用一次生長算法。區別在于上次以窗口中較小值為種子點，這次以窗口中較大值為種子點，且適當縮小窗口大小，當前窗口W′取前一次運算窗口大小的0.5～0.6倍，以使同一窗口中地形變化小，盡可能使高架道路和地面點之間的差異才是窗口中的主要差異。另外，此次生長過程不直接采用原始DSM，而采用對DSM進行中值濾波和開運算后的DSM′，避免異常值的影響，同時去掉細小突出部分，平滑各個地物輪廓，以增強識別結果的整體性。當窗口W′的標準方差大于整個DSM′數據的標準方差的0.6倍時，滿足式(4)，則判斷為種子點；在相對平坦的城區DSM中，滿足這兩個條件的點就是高架道路的種子點。

W′(i,j)>max-level′

(4)

式中，W′(i,j)為窗口中位于(i,j)的像素點；max為當前窗口W′中的最大值；level′公式如下

level′=[max-gt(W′)]/t′

(5)

式中，gt(W′)為窗口W′中最大類間方差法對應的閾值；t′為常數系數，與式(2)的含義一致，只是具體取值不一。

種子點基礎上的生長過程是在整個點集G中完成的，生長準則和終止條件如下

(6)

式中，S′(i,j)為種子點或生長過程中新識別出來的地面點；level′與式(5)中相同；G為已經識別出的地面點集；G′為識別出的高架道路點集；gt(W1)為窗口W1中最大類間方差法對應的閾值；W1是以當前待判斷點為中心的活動窗口，其大小是式(1)中窗口W的一半或更小，以使窗口W1在檢測高架道路過程中，盡可能只包含DSM′中的公路和地面信息。生長準則和終止條件(式(6)與式(3))的不同在于多了一個限制條件：G(i+m,j+n)≥gt(W1)。后一次生長過程是由大值開始，即從高架道路的種子點開始；而在生長過程中，當點(i+m,j+n)小于以點(i,j)為中心的窗口W1的最大類間方差法對應的閾值時，則說明點(i+m,j+n)不屬于高架道路一類，因此終止。

1.4 插 值

在地面點集G(如有高架道路，則G需要減去G′)的基礎上，對非地面待插值的點采用八鄰域中有效值點的均值，一定程度上保持DEM結果的連續性。對于個別八鄰域中沒有有效值點的點，用相應窗口W中地面點的平均值給其賦值。

2 試驗結果與分析

2.1 試驗數據

文中用到兩組試驗驗證數據，一組是地形復雜的山區DSM，原始數據由中國林業科學研究院資源信息所機載LiCHy系統于2014年4月獲取，LiCHy系統包括CCD相機、激光雷達傳感器和AISA Eagle Ⅱ高光譜傳感器[13]。試驗中所用DSM數據由其激光雷達的點云數據通過ENVI LiDAR處理而來，圖像大小為300×300像素，分辨率為1 m。影像區域地形復雜，起伏大，呈現中間和左上方地勢低，右方和下方地勢較高，主要有建筑物、植被、道路等地物。

另一組是2013年電氣和電子工程師協會地球科學與遙感學會(IEEE GRSS)數據融合大賽數據，原始數據包括機載高光譜影像、機載LiDAR數據和分類結果檢驗數據，獲取于2012年6月，覆蓋美國休斯敦大學及附近居民區，由IEEE GRSS提供[14]。高光譜影像的中心波長為364～1046 nm，波段數為144個，空間分辨率為2.5 m。機載LiDAR點云數據地面分辨率同為2.5 m，且兩者相互配準并有相同的覆蓋區域[15]。同樣通過ENVI LiDAR處理得到試驗中所用DSM數據，大小為1905×349像素，分辨率為2.5 m。分類結果檢驗數據共有15類地物。影像中地物種類較多，包括商用房屋、民用房屋、道路、樹等，但整體地勢平坦。本試驗以比賽結果檢驗數據中健康的草、假草、土壤、水、道路、地面停車場、網球場、跑道作為地面點的檢驗數據，以比賽結果檢驗數據中樹、居民區、商業區，以及高架高速公路作為非地面點的檢驗數據。

2.2 評價方法

采用目視解譯、混淆矩陣和制圖精度評價的方法，對本文算法與附近最小值法在識別地面點和提取DEM的精度及可靠性進行比較分析。

2.3 試驗結果

試驗結果如圖2和圖3所示。為了方便視覺比較，高光譜影像為真彩色顯示，并作5%拉伸，灰度圖像為線性拉伸顯示，試驗參數設置詳見表1。

圖2 復雜地形的DSM數據試驗

圖3 平坦城區的DSM數據試驗

試驗組別本文方法附近最小值法窗口W系數t窗口W'系數t'W1窗口W系數t數據一50×5012———15×152.5數據二50×502030×301225×2530×30 8

從目視效果來說，圖2中本文方法的識別結果明顯好于附近最小值法，基本實現了所有地面點的正確識別。對比提取的DEM結果，附近最小值法的結果只能反映整體地形走勢，而本文方法很好地保留了梯田等細節信息，效果明顯好于附近最小值法。同樣的，圖3中本文方法的識別結果明顯優于附近最小值法，本文方法提取的DEM也較好。但從結果(圖3(h))中還能較明顯看出原來建筑物的痕跡，其原因有兩個：①多數建筑地基部分會比周圍地面稍高，因此在DEM結果中，原建筑物位置處的值會大于周圍地面；特別是整個數據的地形平坦，并對DEM結果進行了線性拉伸顯示，擴大了這一細微差異，但不是錯誤。②如圖3(h)中Ⅲ處，雖然實現了正確識別，但由于插值過程采用的是八鄰域的均值，導致在微地貌區域插值效果不理想[16]，存在與周圍有明顯差異的人為痕跡。從圖3(e)可以看出，本文方法很好地識別出了高架道路，防止了因采用區域生長算法而把高架道路分割到地面，并且避免了Liao等[14]在2013年數據融合大賽把圖中Ⅰ的高架高速公路錯分到商業區的錯誤。結果(圖3(e))中區域Ⅱ存在一個明顯錯誤。根據本文算法原理及通過觀察對應位置DSM可知，這一錯誤是由于本文算法采用的DSM數據只反映地物高程信息，而圖中Ⅱ的建筑物(高架高速公路)、植物及地面三者之間在

DSM上存在較好連續性，難以區分。如果結合植被指數(NDVI)作為輔助分析數據，則可避免這一錯誤。試驗效果如圖4所示，基本正確識別出了DSM中全部高架道路。

圖4 結合NDVI后的效果

2.4 定量評價

利用混淆矩陣和制圖精度質量評價指標進行定量評價，統計結果見表2和表3。由于識別的結果只有地面點和非地面點兩類，因此不能單從地面點或非地面點的制圖精度來評價方法優劣，而應該兼顧兩類的精度。對于地形復雜的山區DSM數據，獲取的地面點數量越多，質量越高，得到的DEM結果越好。從表3中可以看出，本文方法的精度都優于附近最小值法，精度都超過90%，是本文方法識別地面點和提取DEM效果好的直接體現。

表2 識別結果的混淆矩陣

本文算法對地形平坦地區和復雜地形的山區都能提取到質量較好的DEM，但還存在兩個局限性：①本算法不適合含有郁閉度非常高的林區的DSM數據。在算法運算中，針對林區部分的DSM數據難以識別出足夠多的地面點，插值結果存在明顯人為痕跡，增加了DEM結果的誤差。尤其是當郁閉度高的林區位于地形復雜的山區時，不能提取出足夠多有效反映地形變化的地面點，致使DEM結果誤差更大。②本算法要求DSM數據中必須包含高架道路附著于地面的起點部分，才能正確識別高架道路。原因在于識別高架道路是基于其的兩個屬性：起點部分附著于地面，中間主要部分高于地面。

表3 識別結果的制圖精度 (%)

3 結 語

本文將最大類間方差法和區域生長算法結合，運用中值濾波和開運算減少噪聲和局部極值對結果的影響，增強分割結果的整體性，實現了圖像分割過程中的種子點、生長準則和終止條件的自適應選擇，從平坦地區和地形復雜的山區DSM中識別地面點和提取DEM。此外，本算法使用兩次區域生長算法實現高架道路的識別。試驗結果及評價指標表明，本算法能夠有效提高DEM的提取精度。

在后期工作中，針對本文算法局限性①，采用連續性更好的插值方法替代八鄰域插值，如雙三次插值法，改善林區部分的插值效果；同時，嘗試對式②的常系數t做到自適應取值，從而進一步增強算法的自適應能力和適用范圍。

參考文獻：

[1] 雷蓉，邱振戈，張士濤.基于遙感影像生成DEM的質量檢查[J].測繪通報,2005(4):36-39.

[2] 湯博麟，蔡海峰，牟玉興.DEM制作檢查及精度問題的解決方法[J].地理空間信息,2016,10(6):76-78.

[3] 馬紅.利用高精度DSM數據提取建筑物輪廓算法研究[J].測繪通報,2015(4):111-124.

[4] 李忠，李梅，杜緒奎.平坦地區DSM到DEM的試驗[J].測繪與空間地理信息,2008,31(2):29-31.

[5] ISMAIL Z，JAAFAR J.DEM Derived from Photogram-metric Generated DSM Using Morphological Filter[C]∥Control and System Graduate Research Colloquium.[S.l.]：IEEE，2013.

[6] 徐海榮，田聯房，陳萍，等.改進的區域生長算法在醫學圖像分割中的應用[J].生物醫學工程研究,2005,24(3):187-190.

[7] 林瑤，田捷.醫學圖像分割方法綜述[J].模式識別與人工智能,2002,15(2):192-204.

[8] KAMDI S，KRISHNA R K.Image Segmentation and Region Growing Algorithm[J].International Journal of Computer Technology & Electronics Engineering,2012,1(2):103-107.

[9] POHLE R，TOENNIES KD.A New Approach for Model-based Adaptive Region Growing in Medical Image Analysis[M].Berlin Heidelberg：Springer，2001：238-246.

[10] LIN Z，JIN J，TALBOT H.Unseeded Region Growing for 3D Image Segmentation[J].Journal of Trauma,2000,24(11):938-945.

[11] 朱紹攀，張書畢.分形地形復雜度研究[J].地理空間信息,2011,9(3):116-119.

[12] 蘇明超，郭書軍.基于區域生長算法的CT圖像肺組織分割[C]∥中國電子學會全國網絡編碼學術年會論文集.北京：國防工業出版社，2008.

[13] 劉怡君，龐勇，廖聲熙，等.機載LiDAR和高光譜融合實現普洱山區樹種分類[J].林業科學研究,2016,29(3):407-412.

[14] LIAO W,PI?URICA A,BELLENS R,et al.Generalized Graph-based Fusion of Hyperspectral and LiDAR Data Using Morphological Features[J].IEEE Geoscience & Remote Sensing Letters,2015,12(3):552-556.

[15] 董彥芳，龐勇，許麗娜，等.高光譜遙感影像與機載 LiDAR 數據融合的地物提取方法研究[J].遙感信息,2014(6):73-76.

[16] 彭儀普，劉文熙.實時交互的多分辨率地形模型[J].武漢大學學報(信息科學版),2003,28(4):480-483.