基于多尺度窗口的DEM局部填洼方法

徐靜波，許捍衛，于艷超

基于多尺度窗口的DEM局部填洼方法

徐靜波1，許捍衛1，于艷超2

（1.河海大學 地球科學與工程學院，江蘇 南京 210098；2.貴州省水利水電勘測設計研究院，貴州 貴陽 550002）

為了去除DEM中的偽洼地，使水系提取結果更加精確，在M&V算法的基礎上提出了一種改進的填洼方法。用局部處理的思想將填洼步驟簡化，減小對DEM的計算范圍和迭代次數，從而達到優化算法的目的；同時借助Arcpy與GDAL庫，將其封裝為一個模塊集成到ArcGIS軟件中，以更高效地解決DEM填洼預處理問題。

填洼；算法；多尺度；局部；評估

DEM是以高程值為基礎的空間分布模型，用來反映區域的地形特征及趨勢[1]。隨著地理信息技術的發展，以及其與水文領域的密切結合，可借助一定算法自動提取研究范圍內的基礎水系。由于插值誤差或采樣誤差，原始DEM會形成許多偽洼地，與實際地形存在不符[2]。洼地是指地表的低洼處，在實際地形中真實洼地是極少的；偽洼地是局部地區的最低點或區域，在水流方向的計算中會導致水不能順利流出，從而對水系的提取造成障礙，影響結果的準確性。因此，在利用DEM提取水系前，必須進行填洼處理。

O'Callaghan J F[3]等提出用小窗口濾波平滑的方法處理DEM，該方法可填平較小的洼地但不適于大范圍的洼地填平。Jenson K[4]等提出J&D算法，主體思想是先填洼地為平地，再逐步墊高并找出潛在出口，確定流向。但這種方法的實現較復雜，對于低分辨率和數據量較小的DEM處理十分有效；對于高分辨率和大數據量的DEM處理就顯得效率低下。另外，在ArcGIS軟件中，J&D算法經過集成被應用于水文分析的填洼中[5]。Moran C J[6]等于1993年提出了一種快速填洼思想（M&V算法），2001年由Planchon O等編程實現[7]。該方法思路簡單，對高低分辨率的DEM普遍適用，且效率提高明顯，但由于算法的可移植性和發展的局限性，M&V算法作為一種先進算法尚未得到廣泛應用。本文借助ArcGIS平臺，利用Python語言和Arcpy站點包，在M&V算法的基礎上提出了一種改進算法，并將其封裝為模塊集成到ArcGIS軟件中，以此來改進填洼算法的應用，旨在更高效、更全面合理地解決DEM填洼預處理問題。

1 M&V算法

1.1 M&V算法結構

M&V算法的主要思想是先用一極大值水面淹沒原始DEM，再迭代移除DEM上多余的水，得到的結果就是填洼后的數據[8]。M&V算法的實現比較簡單，且適用于各種分辨率的DEM。其具體步驟為：

1）除邊界網格保留原值外，對剩余網格賦極大值水面，淹沒原始DEM。極大值默認選取網格最大水面高程值加1。

2）以3×3窗口逐行對所有網格（邊界除外）進行遍歷。以網格C為中間網格，遍歷其8鄰域：將網格C的當前值與鄰域網格值作比較，若鄰域網格值加上極小變量小于C的當前值，則將網格C重新賦值為鄰域網格值加上極小變量。由于填洼是將洼地水面逐漸抬高的過程，因此還需判斷網格原始高程與鄰域網格值加上極小變量的大小關系，若前者大于后者，則將網格新值重新賦值為原始高程。其中，極小變量值可在算法實現中酌情設定，默認值為0.1。

3）多次迭代，即重復步驟2），直到所有網格進行遍歷時，沒有可執行的單元為止，迭代結束。此時網格無限逼近于初始DEM，且洼地被填平。

1.2 算法實現

M&V算法充分考慮了水流流向問題，處理后的DEM滿足在不改變原始DEM結構的基礎上，對其中任意柵格點都能找到一條高程遞減的路徑將水順利排出。算法思想簡單完善，效率明顯高于當前ArcGIS軟件中的填洼工具，該工具使用J&D算法。實現M&V算法偽代碼為：

For each grid cell c of the DEM '遍歷DEM中每個網格

if c is on the border： '邊界網格保持值不變

H(c)=Z(c)

else： '中間網格賦予極大值，默認為DEM所有網格最大值+1

H(c)=max(Z(n))+1

Do

loopFlag=False

For each grid cell c of the DEM

if H(c)!=Z(c)： '中間網格

For each existing neighbor n of c： '遍歷8鄰域

if H(c)＞(H(n)+min)：

H(c)=H(n)+min

loopFlag=True

if Z(c)＞(H(n)+min)：

H(c)=Z(c)

loopFlag=True

Loop While loopFlag=True

2 改進的M&V算法與算法實現

2.1 Python與ArcGIS

Python作為一種免費開源的通用腳本語言，其代碼簡潔明了，能與眾多第三方開源庫進行無縫銜接，編譯過程非常簡單，能直接對源代碼進行調試執行，且Python編寫的腳本程序能跨平臺運行。ArcGIS軟件是一款專門用于地理空間數據處理和分析的軟件，Esri公司在ArcGIS 10.0中引入了Arcpy。Arcpy是基于Python語言的站點包，提供了一種用于集成ArcGIS處理流程和開發工具腳本的豐富的動態環境。在基于Arcpy開發的ArcGIS應用程序和腳本中，可借助大量Python模塊和第三方開源庫，使功能實現更加便捷[9]。

2.2 基于多尺度窗口的局部填洼算法及實施步驟

M&V算法充分利用了邊界條件和DEM的初始值，先將洼地抬高，再通過多次迭代和與鄰域值進行比較，使柵格中的水都能順利流出，從而達到填洼的目的。迭代是M&V算法中最耗時的過程，其時間消耗主要包括兩個方面：迭代次數和每次迭代遍歷的柵格數。

本文從精簡迭代步驟的角度出發，先提取出洼地網格，找到該洼地的相對出口；再利用M&V算法，以半徑R為窗口對每個洼地網格進行計算，使洼地里的水可以順利流出，即以洼地局部迭代替代M&V算法中的全局迭代。為了最大程度地保留原始地形，需制定一個運算規則來確定洼地的迭代窗口R。基于多尺度窗口的局部填洼算法流程如圖1所示。

圖1 基于多尺度窗口的局部填洼算法流程圖

1）遍歷整個網格，標定洼地。洼地分為入流洼地和無向洼地兩種。入流洼地以單個網格形式存在，即其高程值小于其8鄰域網格的高程值；無向洼地通常以片區形式存在，在DEM中表現為網格高程值不低于該洼地所有鄰域網格高程值。

2）確定每個入流洼地的相對出口，此時的搜索半徑R即為該洼地M&V算法的迭代窗口大小。以半徑為R的窗口（初始值默認為5）遍歷洼地的周圍網格，若最小值小于洼地高程值，則標定該網格；若最小值仍高于洼地高程，則以步長為2逐步擴大搜索范圍，直到搜索到符合要求的相對出口為止。記錄找到相應出口時的R值，即為M&V算法迭代窗口的大小。最終洼地的流向就是標定的相對出口，相對出口的標定可作為算法檢驗的輔助參數。考慮到算法的可行性和合理性，R不可能無限增大，設定R的極限值為Rlim，當R增至Rlim仍未找到相對出口時，則增加洼地高程。

3）無向洼地迭代窗口大小的確定。根據無向洼地的分布，用一個最小外接矩形包住所有的洼地網格，迭代窗口即為最小外接矩形的長寬均增加一個長度后的窗口。

4）使用M&V算法對每個洼地進行迭代。

2.3 基于Python的ArcGIS模塊集成

本文開發的填洼算法工具界面如圖2所示。

圖2 基于多尺度窗口的局部填洼算法工具界面

3 性能評估

3.1 算法復雜度分析

算法復雜度的計算涉及時間和存儲兩方面的需求，是一種事前估算，分為時間復雜度和空間復雜度。度量方式是指當問題的尺度以某種方式從1遞增到n時，解決這個問題的算法損耗的時間或占用的空間也以某種方式從1遞增到n，用函數表示為f (n)。本文采用時間復雜度進行度量。根據問題的不同，應考慮最壞的情況和最樂觀的情況，使用O表示：只有存在正整數c和n0使得T(n)≤cg(n)對所有的n≥n0成立，那么稱算法的漸進時間復雜度為T(n)= O(g(n))[7]。

對于一個大小為N的DEM，J&D算法需遍歷所有網格進行填洼，總的時間復雜度為O(N2)；M&V算法每次迭代都至少有一個網格被賦值，迭代最大長度為對角線長度（20.5N），時間復雜度為O(N1.414)；本文提出的先查找洼地再分別填充洼地的做法，在最壞的情況下，才會達到M&V算法的漸進時間復雜度，其時間復雜度接近O(NlogN)。

3.2 運行效率對比

采用3組數據對M&V算法和改進的M&V算法進行對比，90 m分辨率DEM是從地理空間數據云上下載的長江上游數據；30 m分辨率DEM為長江南京段數據，為水下地形測量數據；5 m分辨率DEM為珠江流域CGCS2000成果（表1）。

表1 算法運行效率

3.3 填洼量累計分析

由于M&V算法對高程值的改變有加有減，因此將增加和減少的絕對量之和作為衡量標準。本文算法的3 組數據累計改變高度分別為3 541、46 832和184 035，而M&V算法累計改變高度分別為3 967、50 236、254 793，在保留洼地周圍的地形特征上本文算法略優（表2）。

表2 填洼量累計分析結果

4 結 語

本文從DEM填洼算法需求的本源出發，在理解時下通用算法的基礎上，對M&V算法進行了改進，提出了一種新的基于多窗口尺度的局部填洼算法。算法考慮到每個洼地的地形特征，并針對區域性地形采用不同尺度的窗口對洼地作處理，減少了原始M&V算法的迭代次數和每一次迭代的數據量，并最大程度地保留了洼地周圍的地形特征。實驗證明，該方法比傳統填洼算法具有更高的填洼效率。

[1] 任立良,劉新仁. 基于DEM的水文物理過程模擬[J].地理研究, 2000,19(4)：369-376

[2] Martz L W, Garbrecht J Numerical Definition of Drainage Network and Subcatchment Areas from Digital Elevation Models[J]. Computers & Geosciences,1992,18(6)：741-761

[3] O'Callaghan J F, Mark D M. The Extraction of Drainage Networks from Digital Elevation Data[J].Computer Vision Graphics and Image Processing,1984,27(3)：323-344

[4] Jenson K, Domingue O. Extracting Topographic Structure from Digital Elevation Data for Geographic Information System Analysis[J]. Sensing,1988,54(11)：1 593-1 600

[5] 李輝,陳曉玲,張利華,等.基于三方向搜索的DEM中洼地處理方法[J].水利科學進展,2009,20(4)：473-479

[6] Moran C J, Vezina G. Visualizing Soil Surface and Crop Residues[J]. IEEE Computer Graphics and Applications, 1993,13(2)：40-47

[7] Planchon O, Darboux F. A Simple and Versatile Algorithm to Fill the Depressions of Digital Elevation Models[J]. Catena,2002, 46(2)：159-176

[8] 楊邦,任立良,賀穎慶.基于快速排序的數字高程模型分級填洼算法[J].計算機應用,2009,29(11)：3 161-3 164

[9] 彭海波,向洪普.基于Python 的空間數據批量處理方法[J].測繪與空間地理信息,2011,34(4)：81-82

[10] 殷人昆.數據結構(用面向對象方法與C++語言描述)[M].北京：清華大學出版社,1999

P208

B

1672-4623（2017）05-0098-03

10.3969/j.issn.1672-4623.2017.0053.0

徐靜波，碩士研究生，主要從事GIS開發與應用。

2015-08-31。

項目來源：國家自然科學基金資助項目（41101374、41101308）；水利部公益性行業科研專項經費資助項目（201201025）。