徐靜波，許捍衛(wèi)，于艷超

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

徐靜波1，許捍衛(wèi)1，于艷超2

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

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

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

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

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

1 M&V算法

1.1 M&V算法結(jié)構(gòu)

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

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

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

3）多次迭代，即重復(fù)步驟2），直到所有網(wǎng)格進(jìn)行遍歷時，沒有可執(zhí)行的單元為止，迭代結(jié)束。此時網(wǎng)格無限逼近于初始DEM，且洼地被填平。

1.2 算法實現(xiàn)

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

For each grid cell c of the DEM '遍歷DEM中每個網(wǎng)格

if c is on the border： '邊界網(wǎng)格保持值不變

H(c)=Z(c)

else： '中間網(wǎng)格賦予極大值，默認(rèn)為DEM所有網(wǎng)格最大值+1

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

Do

loopFlag=False

For each grid cell c of the DEM

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

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 改進(jìn)的M&V算法與算法實現(xiàn)

2.1 Python與ArcGIS

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

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

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

本文從精簡迭代步驟的角度出發(fā)，先提取出洼地網(wǎng)格，找到該洼地的相對出口；再利用M&V算法，以半徑R為窗口對每個洼地網(wǎng)格進(jìn)行計算，使洼地里的水可以順利流出，即以洼地局部迭代替代M&V算法中的全局迭代。為了最大程度地保留原始地形，需制定一個運(yùn)算規(guī)則來確定洼地的迭代窗口R?；诙喑叨却翱诘木植刻钔菟惴鞒倘鐖D1所示。

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

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

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

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

4）使用M&V算法對每個洼地進(jìn)行迭代。

2.3 基于Python的ArcGIS模塊集成

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

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

3 性能評估

3.1 算法復(fù)雜度分析

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

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

3.2 運(yùn)行效率對比

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

表1 算法運(yùn)行效率

3.3 填洼量累計分析

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

表2 填洼量累計分析結(jié)果

4 結(jié) 語

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

[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].水利科學(xué)進(jìn)展,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] 楊邦,任立良,賀穎慶.基于快速排序的數(shù)字高程模型分級填洼算法[J].計算機(jī)應(yīng)用,2009,29(11)：3 161-3 164

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

[10] 殷人昆.數(shù)據(jù)結(jié)構(gòu)(用面向?qū)ο蠓椒ㄅcC++語言描述)[M].北京：清華大學(xué)出版社,1999

P208

B

1672-4623（2017）05-0098-03

10.3969/j.issn.1672-4623.2017.0053.0

徐靜波，碩士研究生，主要從事GIS開發(fā)與應(yīng)用。

2015-08-31。

項目來源：國家自然科學(xué)基金資助項目（41101374、41101308）；水利部公益性行業(yè)科研專項經(jīng)費(fèi)資助項目（201201025）。