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

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

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

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

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

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

填洼；算法；多尺度；局部；評(píng)估

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

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

1 M&V算法

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

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

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

2）以3×3窗口逐行對(duì)所有網(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)格值加上極小變量。由于填洼是將洼地水面逐漸抬高的過(guò)程，因此還需判斷網(wǎng)格原始高程與鄰域網(wǎng)格值加上極小變量的大小關(guān)系，若前者大于后者，則將網(wǎng)格新值重新賦值為原始高程。其中，極小變量值可在算法實(shí)現(xiàn)中酌情設(shè)定，默認(rèn)值為0.1。

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

1.2 算法實(shí)現(xiàn)

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

For each grid cell c of the DEM '遍歷DEM中每個(gè)網(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算法與算法實(shí)現(xiàn)

2.1 Python與ArcGIS

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

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

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

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

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

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

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

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

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

2.3 基于Python的ArcGIS模塊集成

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

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

3 性能評(píng)估

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

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

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

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

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

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

3.3 填洼量累計(jì)分析

由于M&V算法對(duì)高程值的改變有加有減，因此將增加和減少的絕對(duì)量之和作為衡量標(biāo)準(zhǔn)。本文算法的3 組數(shù)據(jù)累計(jì)改變高度分別為3 541、46 832和184 035，而M&V算法累計(jì)改變高度分別為3 967、50 236、254 793，在保留洼地周?chē)牡匦翁卣魃媳疚乃惴詢(xún)?yōu)（表2）。

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

4 結(jié) 語(yǔ)

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

[1] 任立良,劉新仁. 基于DEM的水文物理過(guò)程模擬[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].計(jì)算機(jī)應(yīng)用,2009,29(11)：3 161-3 164

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

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

P208

B

1672-4623（2017）05-0098-03

10.3969/j.issn.1672-4623.2017.0053.0

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

2015-08-31。

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