一種面向海量數(shù)字高程模型數(shù)據(jù)的洪水淹沒(méi)區(qū)快速生成算法

沈定濤，王結(jié)臣，張 煜，黃 俊，胡承芳

1．江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210046；2．南京大學(xué) 地理信息科學(xué)系，江蘇 南京 210046；3．長(zhǎng)江水利委員會(huì)長(zhǎng)江科學(xué)院空間信息技術(shù)應(yīng)用研究所，湖北 武漢 430010

一種面向海量數(shù)字高程模型數(shù)據(jù)的洪水淹沒(méi)區(qū)快速生成算法

沈定濤1，2，3，王結(jié)臣1，2，張 煜3，黃 俊3，胡承芳3

1．江蘇省地理信息技術(shù)重點(diǎn)實(shí)驗(yàn)室，江蘇 南京 210046；2．南京大學(xué) 地理信息科學(xué)系，江蘇 南京 210046；3．長(zhǎng)江水利委員會(huì)長(zhǎng)江科學(xué)院空間信息技術(shù)應(yīng)用研究所，湖北 武漢 430010

常見(jiàn)種子點(diǎn)填充算法在實(shí)現(xiàn)數(shù)字高程模型（digital elevation model，DEM）數(shù)據(jù)下的洪水淹沒(méi)區(qū)生成時(shí)，具有難以處理大數(shù)據(jù)量、過(guò)多的遞歸計(jì)算易導(dǎo)致算法效率較低等缺點(diǎn)。針對(duì)此問(wèn)題，本文提出一種面向海量DEM數(shù)據(jù)的洪水淹沒(méi)區(qū)生成算法分塊壓縮追蹤法，該算法采用條帶分塊和實(shí)時(shí)柵格壓縮存儲(chǔ)技術(shù)，以解決海量地形數(shù)據(jù)下的淹沒(méi)分析計(jì)算問(wèn)題。最后，通過(guò)將本算法與常見(jiàn)種子點(diǎn)填充算法和分塊種子點(diǎn)填充算法進(jìn)行了對(duì)比測(cè)試，試驗(yàn)結(jié)果表明本算法不僅較好地解決了海量DEM數(shù)據(jù)下的洪水淹沒(méi)區(qū)生成問(wèn)題，與常規(guī)種子點(diǎn)填充算法和分塊種子點(diǎn)填充算法相比亦具有較高的計(jì)算效率。

海量地形數(shù)據(jù)；淹沒(méi)分析；種子填充；壓縮存儲(chǔ)；數(shù)字高程模型

1 前 言

洪水淹沒(méi)分析是進(jìn)行水文預(yù)測(cè)預(yù)報(bào)、洪水災(zāi)害評(píng)估的一項(xiàng)重要內(nèi)容，也是GIS地形三維仿真系統(tǒng)中的一個(gè)重要功能［1－3］。洪水淹沒(méi)分析主要基于DEM數(shù)據(jù)，結(jié)合水文水動(dòng)力學(xué)構(gòu)建的洪水演進(jìn)模型或根據(jù)設(shè)定的靜態(tài)“水位面”進(jìn)行計(jì)算。前者可以較精確地模擬洪水在不同時(shí)間段下的淹沒(méi)區(qū)變化過(guò)程，但是洪水演進(jìn)模型的構(gòu)建需要較多的水文水力參數(shù)條件，其處理較為復(fù)雜［4－8］。基于靜水位下的洪水淹沒(méi)分析是對(duì)洪水淹沒(méi)過(guò)程達(dá)到最終平衡狀態(tài)的近似模擬，盡管沒(méi)有洪水演進(jìn)模型方法精確，但其計(jì)算簡(jiǎn)單、實(shí)用，能夠快速地確定洪水淹沒(méi)區(qū)和水深分布，得到了廣泛的應(yīng)用。靜水位下洪水淹沒(méi)分析主要分為無(wú)源淹沒(méi)和有源淹沒(méi)兩種，其中無(wú)源淹沒(méi)相對(duì)簡(jiǎn)單，只需遍歷全部網(wǎng)格即可得出淹沒(méi)范圍，而有源淹沒(méi)需要考慮如環(huán)形山、堤壩等障礙，其計(jì)算相對(duì)復(fù)雜，是洪水淹沒(méi)分析算法研究的重點(diǎn)?，F(xiàn)有洪水淹沒(méi)分析算法的研究主要來(lái)源于圖像處理中的種子填充法及其變形，通過(guò)給定淹沒(méi)起算點(diǎn)，在DEM網(wǎng)格中按照四鄰域或八鄰域進(jìn)行擴(kuò)散，可以較好地模擬實(shí)際的地表徑流淹沒(méi)過(guò)程，以及自動(dòng)提取淹沒(méi)區(qū)。如文獻(xiàn)［9—10］提出的基于GIS的復(fù)雜地形淹沒(méi)分析及其災(zāi)害評(píng)估；文獻(xiàn)［11］提出的基于GIS格網(wǎng)模型的洪水淹沒(méi)分析；文獻(xiàn)［12］提出的基于DEM的“環(huán)形”洪水淹沒(méi)算法；文獻(xiàn)［13］提出的基于DEM的洪水淹沒(méi)分析等。

由于無(wú)法預(yù)先判斷洪水淹沒(méi)的大致范圍，常規(guī)種子點(diǎn)填充算法必須將DEM數(shù)據(jù)一次性載入內(nèi)存進(jìn)行計(jì)算，其所能載入的數(shù)據(jù)量依賴于計(jì)算機(jī)內(nèi)存配置。一般三維仿真系統(tǒng)中用于地形表達(dá)的DEM數(shù)據(jù)常常達(dá)到幾GB甚至幾十GB，這在目前個(gè)人計(jì)算機(jī)內(nèi)存配置條件下是難以全部讀入內(nèi)存的，當(dāng)進(jìn)行洪水淹沒(méi)分析時(shí)，容易受到內(nèi)存的限制導(dǎo)致計(jì)算失敗。另一方面，當(dāng)DEM數(shù)據(jù)較大，淹沒(méi)范圍較廣時(shí)，種子填充算法中大量的遞歸操作常導(dǎo)致算法效率急劇降低，同時(shí)容易出現(xiàn)堆棧溢出等問(wèn)題，這極大地限制了種子點(diǎn)填充淹沒(méi)分析算法的應(yīng)用范圍。因此，在進(jìn)行大規(guī)模DEM數(shù)據(jù)的洪水淹沒(méi)分析時(shí)，如何避免大量耗時(shí)的文件讀寫操作，同時(shí)又能在內(nèi)存中高效地處理大數(shù)據(jù)量DEM，是實(shí)現(xiàn)高效、大范圍、高精度和自動(dòng)化洪水淹沒(méi)分析的關(guān)鍵。

2 算法設(shè)計(jì)思路

2．1 分塊種子點(diǎn)填充算法設(shè)計(jì)思路

針對(duì)上述問(wèn)題，研究人員亦提出了一些優(yōu)化算法，如文獻(xiàn)［11］提出將DEM數(shù)據(jù)轉(zhuǎn)換為不規(guī)則多邊形格網(wǎng)和三角網(wǎng)，在多邊形和三角形格網(wǎng)的基礎(chǔ)上按照種子點(diǎn)擴(kuò)散方法進(jìn)行淹沒(méi)范圍計(jì)算。文獻(xiàn)［14］提出為DEM數(shù)據(jù)構(gòu)建瓦片金字塔，在大區(qū)域場(chǎng)景下用上層金字塔數(shù)據(jù)進(jìn)行淹沒(méi)分析，以實(shí)現(xiàn)快速三維實(shí)時(shí)顯示。這類方法通過(guò)對(duì)DEM數(shù)據(jù)進(jìn)行簡(jiǎn)化，以簡(jiǎn)化數(shù)據(jù)實(shí)現(xiàn)快速種子點(diǎn)填充計(jì)算，在一定程度上降低了算法對(duì)內(nèi)存的依賴，提升了計(jì)算速度，但是淹沒(méi)區(qū)范圍計(jì)算精度會(huì)直接受到簡(jiǎn)化數(shù)據(jù)影響，同時(shí)針對(duì)海量DEM數(shù)據(jù)的處理仍然有限。

一種較為有效的策略是對(duì)DEM數(shù)據(jù)進(jìn)行條帶分塊讀取，并實(shí)時(shí)計(jì)算淹沒(méi)區(qū)。以分塊種子點(diǎn)填充算法為例，其基本思路為：首先對(duì)原DEM數(shù)據(jù)進(jìn)行條帶劃分，在每個(gè)條帶上采用種子點(diǎn)填充算法計(jì)算所有潛在淹沒(méi)區(qū)并編號(hào)，然后對(duì)相鄰條帶間潛在淹沒(méi)區(qū)執(zhí)行連通域檢測(cè)，查找屬于同一連通域的潛在淹沒(méi)區(qū)并構(gòu)建關(guān)聯(lián)表，最后通過(guò)初始種子點(diǎn)所在淹沒(méi)區(qū)編號(hào)和關(guān)聯(lián)表提取整個(gè)淹沒(méi)區(qū)。如圖1所示，DEM數(shù)據(jù)被劃分為3個(gè)條帶，圖中灰色填充網(wǎng)格代表各條帶上的潛在淹沒(méi)區(qū)。通過(guò)對(duì)比相鄰條帶間對(duì)應(yīng)網(wǎng)格單元所在的淹沒(méi)區(qū)編號(hào)，執(zhí)行連通域檢測(cè)并構(gòu)建關(guān)聯(lián)表，如條帶2上的1、2號(hào)淹沒(méi)區(qū)和條帶3上3號(hào)淹沒(méi)區(qū)屬于同一個(gè)連通域，它們是相互關(guān)聯(lián)的。最后，根據(jù)種子點(diǎn)所對(duì)應(yīng)的1號(hào)淹沒(méi)區(qū)，逐個(gè)條帶查找關(guān)聯(lián)表信息并提取其他關(guān)聯(lián)淹沒(méi)區(qū)，構(gòu)建洪水淹沒(méi)范圍。

圖1 分塊種子點(diǎn)填充法示意圖Fig．1 Strip－divide seed filling algorithm

上述連通域標(biāo)記方法在執(zhí)行條帶種子點(diǎn)填充算法后可以及時(shí)清空分塊數(shù)據(jù)，降低了對(duì)內(nèi)存的依賴，使得處理海量柵格數(shù)據(jù)成為可能，在遙感及圖像處理領(lǐng)域得到了一定的應(yīng)用［15－18］。但將其用于淹沒(méi)分析時(shí)，主要弊端在于：①為進(jìn)行連通域判定，必須在條帶中執(zhí)行種子點(diǎn)填充算法提取所有潛在淹沒(méi)區(qū)，而大部分潛在淹沒(méi)區(qū)都不是實(shí)際淹沒(méi)范圍，這一處理過(guò)程較為低效；②條帶上種子填充計(jì)算結(jié)果無(wú)法保存在內(nèi)存中，因此需要預(yù)先生成一個(gè)與原DEM大小相同的空值柵格文件，及時(shí)寫入分塊條帶的計(jì)算結(jié)果，而在進(jìn)行連通域判定時(shí)則要多次讀取此淹沒(méi)結(jié)果文件；③當(dāng)淹沒(méi)區(qū)范圍遠(yuǎn)小于原DEM范圍時(shí)，大量條帶數(shù)據(jù)的讀取、計(jì)算和連通域判定是沒(méi)有必要的。

2．2 分塊壓縮追蹤算法設(shè)計(jì)思路

考慮到條帶中每個(gè)網(wǎng)格單元只有潛在淹沒(méi)和未淹沒(méi)兩種狀態(tài)，若在讀取條帶數(shù)據(jù)后直接將每行中高程值低于設(shè)定水位的相鄰網(wǎng)格單元按照游程編碼方法進(jìn)行壓縮，則可以較大程度上降低柵格數(shù)據(jù)量，進(jìn)而將潛在淹沒(méi)區(qū)存儲(chǔ)進(jìn)內(nèi)存。為避免盲目地讀取所有條帶數(shù)據(jù)，可以先讀取種子點(diǎn)所在條帶數(shù)據(jù)，實(shí)現(xiàn)壓縮存儲(chǔ)并以種子點(diǎn)所在壓縮單元進(jìn)行邊界追蹤標(biāo)記，然后向上和向下擴(kuò)展讀取新條帶并壓縮，每讀取一個(gè)新條帶都以其相鄰條帶頂柵格行或底柵格行中被追蹤過(guò)的壓縮單元兩側(cè)作為追蹤邊，在新條帶上執(zhí)行邊界追蹤，最后提取柵格場(chǎng)上被追蹤過(guò)的所有壓縮單元即可構(gòu)建淹沒(méi)區(qū)，可以稱之為分塊壓縮追蹤法。

如圖2所示，首先讀取條帶2并壓縮存儲(chǔ)，以種子點(diǎn)所在壓縮單元兩側(cè)追蹤得到淹沒(méi)區(qū)1，該條帶底柵格行和頂柵格行都有被追蹤過(guò)的壓縮單元，因此分別讀取條帶1和條帶3并進(jìn)行壓縮存儲(chǔ)，分別以條帶2底柵格行和頂柵格行上被追蹤過(guò)的壓縮單元兩側(cè)作為追蹤邊向下和向上追蹤得到條帶1的淹沒(méi)區(qū)1和條帶3的淹沒(méi)區(qū)3并進(jìn)而追蹤得到條帶2上的淹沒(méi)區(qū)2。因條帶1的底柵格行和條帶3的頂柵格行都沒(méi)有被追蹤過(guò)的壓縮單元，停止讀取新條帶，最后提取所有被追蹤過(guò)的壓縮單元生成洪水淹沒(méi)范圍。

與分塊種子點(diǎn)填充法相比，本算法的主要優(yōu)勢(shì)在于：①對(duì)分塊條帶直接進(jìn)行柵格壓縮存儲(chǔ)，避免了效率低下的種子點(diǎn)填充計(jì)算，而將潛在淹沒(méi)區(qū)存儲(chǔ)于內(nèi)存比存儲(chǔ)于文件效率更高；②條帶上潛在淹沒(méi)區(qū)以壓縮單元形式存儲(chǔ)于內(nèi)存，針對(duì)某壓縮單元的邊界追蹤可以跨越多個(gè)條帶，從而規(guī)避了復(fù)雜的連通域標(biāo)記過(guò)程；③按照向上和向下擴(kuò)展讀取新條帶并及時(shí)追蹤壓縮單元，讀取的條帶數(shù)目與實(shí)際淹沒(méi)范圍直接相關(guān)，沒(méi)有洪水淹沒(méi)的條帶則不會(huì)讀入內(nèi)存，避免了大量無(wú)效條帶數(shù)據(jù)的讀取與計(jì)算。

3 分塊壓縮追蹤法關(guān)鍵技術(shù)及其實(shí)現(xiàn)

3．1 算法總體流程

按照上述算法設(shè)計(jì)思路，分塊壓縮追蹤法的總體算法流程圖如圖3所示。

3．2 條帶壓縮與邊界追蹤標(biāo)記

條帶壓縮與邊界追蹤標(biāo)記包含了算法流程中的步驟（1）—步驟（3）。DEM數(shù)據(jù)條帶劃分示意圖如圖4所示，假定某DEM數(shù)據(jù)有24行、32列，設(shè)定條帶規(guī)格為8行，則DEM數(shù)據(jù)被劃分為3個(gè)條帶。為便于說(shuō)明，將網(wǎng)格上的高程值以整數(shù)進(jìn)行表達(dá)，以1—8的數(shù)字代表高程值，數(shù)字越大表示高程值越高。

圖4 柵格條帶分塊示意圖Fig．4 Raster strip－divide diagram

種子點(diǎn)所在網(wǎng)格位于DEM數(shù)據(jù)的第10行、第14列，屬于條帶2的第2行第14列，其高程值為3。假定淹沒(méi)水位為4，圖中以灰色填充的網(wǎng)格為種子點(diǎn)所在條帶的潛在淹沒(méi)范圍，其高程值均小于或等于設(shè)定水位值。

條帶數(shù)據(jù)的壓縮方法為：

（1）在內(nèi)存中生成一個(gè)以鏈表為元素的數(shù)組，數(shù)組大小為原DEM柵格行數(shù)目，鏈表的結(jié)點(diǎn)為壓縮單元，該鏈表結(jié)構(gòu)用于條帶數(shù)據(jù)壓縮后存儲(chǔ)對(duì)應(yīng)柵格行的所有壓縮單元。

（2）在條帶上按行遍歷，依次將柵格行上高程值小于或等于設(shè)定淹沒(méi)水位的相鄰網(wǎng)格采用游程編碼方法進(jìn)行壓縮，生成一個(gè)壓縮單元結(jié)構(gòu)，該壓縮單元標(biāo)記了潛在淹沒(méi)網(wǎng)格在此柵格行上跨越的起始列和終止列。

（3）每生成一個(gè)新的壓縮單元，則將該壓縮單元插入對(duì)應(yīng)行上壓縮單元鏈表尾部，最后完成整個(gè)條帶數(shù)據(jù)的壓縮存儲(chǔ)。圖5是種子點(diǎn)所在條帶壓縮存儲(chǔ)示意，以圖中第9行的壓縮單元為例，該壓縮單元標(biāo)記了1和18兩個(gè)值，表示該壓縮單元跨越了第1—18列。

圖5 條帶壓縮與邊界追蹤Fig．5 Strip compression and boundary tracing

條帶數(shù)據(jù)的壓縮存儲(chǔ)尚未標(biāo)記每個(gè)壓縮單元與淹沒(méi)源種子點(diǎn)的連通性，可以通過(guò)對(duì)種子點(diǎn)所在的壓縮單元進(jìn)行邊界追蹤以標(biāo)記與種子點(diǎn)連通的其他壓縮單元。主要步驟為：首先為每個(gè)鏈表設(shè)置一個(gè)標(biāo)記數(shù)組，數(shù)組元素的個(gè)數(shù)是鏈表中壓縮單元的兩倍，數(shù)組依次存儲(chǔ)了每個(gè)壓縮單元的左側(cè)和右側(cè)兩個(gè)標(biāo)記，默認(rèn)將所有的標(biāo)記數(shù)組元素設(shè)置為false，即表示所有壓縮單元兩側(cè)都未被追蹤過(guò)。其次，查找初始種子點(diǎn)所在的壓縮單元，分別以該壓縮單元左側(cè)和右側(cè)邊界作為起始追蹤邊，按照四方向追蹤法［19－20］進(jìn)行柵格邊界追蹤，每追蹤到新的壓縮單元左側(cè)或右側(cè)，就將該行上對(duì)應(yīng)的標(biāo)記數(shù)組元素值由false設(shè)置為true，當(dāng)追蹤到某側(cè)標(biāo)記為true的壓縮單元?jiǎng)t停止追蹤。最后，循環(huán)遍歷初始種子點(diǎn)所在柵格行上所有壓縮單元，若某壓縮單元一側(cè)被追蹤過(guò)而另一側(cè)未被追蹤，則以未被追蹤一側(cè)作為起始追蹤邊進(jìn)行邊界追蹤，直到該行上所有被追蹤過(guò)的壓縮單元滿足兩側(cè)都被追蹤。圖5顯示了壓縮單元邊界追蹤結(jié)果，追蹤邊界如圖中粗線所示。

3．3 依條帶擴(kuò)展存儲(chǔ)

依條帶擴(kuò)展存儲(chǔ)包含了算法流程中的步驟（5）—步驟（12），此步驟是一個(gè)循環(huán)過(guò)程，其循環(huán)結(jié)束條件為步驟（13）。因?yàn)榭梢灶A(yù)先知道當(dāng)前已讀取過(guò)哪些條帶數(shù)據(jù)，為便于描述，按照每個(gè)條帶所處的空間位置，將讀取的位于最上方的條帶稱為最上條帶，將位于最下方的條帶稱為最下條帶，將每個(gè)條帶中最上一柵格行稱為該條帶的頂柵格行，最下一柵格行稱為該條帶的底柵格行。在完成種子點(diǎn)所在條帶的壓縮存儲(chǔ)后，可以將該條帶同時(shí)設(shè)定為最上條帶和最下條帶。

此時(shí)可以向上和向下擴(kuò)展讀取新的條帶，并按照前述方法進(jìn)行條帶壓縮與邊界追蹤。其原理為：對(duì)最上條帶和最下條帶而言，需要繼續(xù)向上和向下讀取新條帶的唯一條件就是其頂柵格行或底柵格行上有被追蹤過(guò)的壓縮單元，表示淹沒(méi)區(qū)還未到達(dá)實(shí)際邊界。當(dāng)最上條帶頂柵格行上有被追蹤過(guò)的壓縮單元時(shí)，可以將該行所有被追蹤過(guò)壓縮單元的兩側(cè)標(biāo)記為起始追蹤邊，讀取最上條帶的上一條帶作為新的最上條帶，并進(jìn)行壓縮存儲(chǔ)，此時(shí)以標(biāo)記的起始追蹤邊向上追蹤即可得到當(dāng)前最上條帶的淹沒(méi)區(qū)，同理對(duì)最下條帶執(zhí)行此過(guò)程。不斷循環(huán)執(zhí)行最上條帶和最下條帶的擴(kuò)展存儲(chǔ)，當(dāng)最上條帶頂柵格行和最下條帶底柵格行都沒(méi)有被追蹤的壓縮單元，或都已到達(dá)DEM邊界則停止循環(huán)。

如圖5所示，鏈表數(shù)組中只存儲(chǔ)了一個(gè)條帶的壓縮數(shù)據(jù)，其頂柵格行（第16行）有兩個(gè)壓縮單元且都被追蹤過(guò)，因此需要從原DEM中讀取當(dāng)前條帶的上一條帶，并執(zhí)行游程壓縮存儲(chǔ)。圖6顯示了讀取上一條帶后的鏈表數(shù)組變化，在第17—24行上新增加了壓縮單元。此時(shí)分別以第16行被追蹤過(guò)邊界的兩個(gè)壓縮單元的左側(cè)和右側(cè)為起始追蹤邊向上進(jìn)行邊界追蹤，圖中以粗線條顯示的為邊界追蹤結(jié)果。當(dāng)前最下條帶底柵格行位于第9行，有一個(gè)壓縮單元，且兩側(cè)已經(jīng)被追蹤過(guò)，需要讀取DEM中對(duì)應(yīng)的下一條帶并進(jìn)行壓縮存儲(chǔ)，結(jié)果如圖7所示。此時(shí)以第9行壓縮單元兩側(cè)作為追蹤邊，執(zhí)行向下追蹤，由于第10行上壓縮單元與第9行壓縮單元并不相連，追蹤停止。

圖6 上一條帶壓縮存儲(chǔ)及邊界追蹤Fig．6 Upper strip compression and boundary tracing

圖7 下一條帶壓縮存儲(chǔ)及邊界追蹤Fig．7 Lower strip compression and boundary tracing

3．4 淹沒(méi)區(qū)孤島處理

淹沒(méi)區(qū)孤島處理對(duì)應(yīng)算法流程圖中步驟（14）。在完成擴(kuò)展條帶壓縮存儲(chǔ)及邊界追蹤標(biāo)記后，需要追蹤部分未被標(biāo)記的孤島。如圖7所示，淹沒(méi)區(qū)中仍然有兩個(gè)孤島未被追蹤出來(lái)，其中一個(gè)孤島橫跨第12—15行，另一個(gè)孤島橫跨第19—22行?？梢灾鹦斜闅v追蹤標(biāo)記數(shù)組，并兩兩取值對(duì)比，該值代表對(duì)應(yīng)壓縮單元的左側(cè)和右側(cè)追蹤情況。若兩個(gè)標(biāo)記值分別為true和false，則需對(duì)標(biāo)記值為false的一側(cè)進(jìn)行追蹤。如圖7中第12行第一個(gè)壓縮單元，其左側(cè)已經(jīng)被追蹤，但是右側(cè)未被追蹤過(guò)，此時(shí)以該側(cè)為起始追蹤邊進(jìn)行邊界追蹤標(biāo)記，將孤島邊界追蹤出來(lái)。按照此方法，最后追蹤結(jié)果如圖8所示，此時(shí)整個(gè)淹沒(méi)區(qū)外側(cè)邊界和內(nèi)側(cè)邊界已經(jīng)全部追蹤出來(lái)。

圖8 淹沒(méi)區(qū)孤島提取Fig．8 Island extracting in flood region

4 試驗(yàn)與分析

基于上述算法設(shè)計(jì)原理，筆者采用C＃語(yǔ)言實(shí)現(xiàn)了洪水淹沒(méi)區(qū)生成算法，并在“金沙江水文泥沙信息三維可視化系統(tǒng)”中得到應(yīng)用。圖9是金沙江流域研究區(qū)DEM數(shù)據(jù)，該數(shù)據(jù)包含的網(wǎng)格行列數(shù)為19 719行×19 454列，數(shù)據(jù)文件大小為1．43 GB，網(wǎng)格寬度25 m，圖10顯示了淹沒(méi)源水位上升50米時(shí)的洪水淹沒(méi)區(qū)范圍。

圖9 研究區(qū)DEMFig．9 Study area of DEM data

圖10 淹沒(méi)范圍示意圖Fig．10 Flood region diagram

為驗(yàn)證算法效率，筆者分別采用常規(guī)種子填充算法、分塊種子點(diǎn)填充算法和本文提出的分塊壓縮追蹤法進(jìn)行了測(cè)試，其中常規(guī)種子填充算法因?yàn)閮?nèi)存不足導(dǎo)致計(jì)算失敗。表1給出了分塊壓縮追蹤法與分塊種子點(diǎn)填充法的算法耗時(shí)對(duì)比。測(cè)試結(jié)果表明，在相同的條帶規(guī)格下，分塊壓縮追蹤法耗時(shí)遠(yuǎn)小于分塊種子填充法耗時(shí)。試驗(yàn)中設(shè)定淹沒(méi)水深為定值，將單個(gè)條帶的行數(shù)從50逐步擴(kuò)展到2000，當(dāng)一次讀入的條帶范圍變大時(shí)，由于減少了分塊數(shù)量，同時(shí)降低了分塊間連通性判定的次數(shù)，因此分塊種子點(diǎn)填充法耗時(shí)是隨著分塊條帶數(shù)據(jù)的逐步變大而降低的，其最佳算法運(yùn)行條件為分塊條帶能剛好被計(jì)算機(jī)一次性讀入內(nèi)存。分塊壓縮追蹤法是依條帶進(jìn)行擴(kuò)展讀取的，其讀取的條帶總體范圍與實(shí)際淹沒(méi)范圍直接相關(guān)，盡管在條帶分塊更小時(shí)讀取的條帶數(shù)目更多，但是其實(shí)際讀取的數(shù)據(jù)范圍可能更小，即文件I／O耗時(shí)更少，在條帶分塊較大時(shí)這種差別更明顯。而壓縮單元邊界追蹤在不同分塊情形下的差別并不大，因?yàn)檫吔缱粉櫺枰阉鞯膲嚎s單元數(shù)目是一致的。隨著條帶行數(shù)的逐步增大，算法總耗時(shí)增加并不顯著，相反在條帶行數(shù)較少時(shí)算法耗時(shí)更少，條帶行數(shù)從50擴(kuò)大到2000時(shí)算法耗時(shí)僅增加了26 s，而讀取每個(gè)分塊條帶數(shù)據(jù)所占據(jù)的內(nèi)存則從8 MB（50行×2萬(wàn)列×8字節(jié)）增加到320 MB（2000行×2萬(wàn)列×8字節(jié)）。在條帶分塊策略方面，只要分塊數(shù)據(jù)能夠一次性讀入內(nèi)存都是可以滿足算法要求的，但分塊條帶更小時(shí)可減少文件I／O耗時(shí)，且在一定程度上降低了內(nèi)存峰值，是較優(yōu)的選擇。

表1 分塊壓縮追蹤法與分塊種子點(diǎn)填充法算法耗時(shí)對(duì)比Tab．1 Time comparison of divide compression tracing and stripe－divide seed filling algorithm

在內(nèi)存占用方面，常規(guī)種子點(diǎn)填充算法需要讀取整個(gè)DEM數(shù)據(jù)，以本算法測(cè)試數(shù)據(jù)為例，其需要占用內(nèi)存約3．2 GB（2萬(wàn)行×2萬(wàn)列×8字節(jié)）的數(shù)據(jù)量，這在目前個(gè)人計(jì)算機(jī)上是難以做到的。分塊種子點(diǎn)填充法與分塊壓縮追蹤法較好地解決了這一問(wèn)題，在讀取條帶數(shù)據(jù)后通過(guò)計(jì)算潛在淹沒(méi)區(qū)，可以及時(shí)清空分塊條帶數(shù)據(jù)所占用的內(nèi)存，其內(nèi)存占用的峰值主要來(lái)自于條帶數(shù)據(jù)的讀入，如條帶行數(shù)為1000時(shí)占用的內(nèi)存也僅160 MB。分塊壓縮追蹤法將潛在淹沒(méi)區(qū)都存儲(chǔ)在內(nèi)存中，測(cè)試中整個(gè)柵格場(chǎng)共生成了258 557個(gè)壓縮單元，按照每個(gè)壓縮單元存儲(chǔ)起始列和終止列兩個(gè)字段計(jì)算，需要8字節(jié)的數(shù)據(jù)量，則總數(shù)據(jù)量?jī)H為2 MB。若按此柵格數(shù)據(jù)壓縮比進(jìn)行計(jì)算，當(dāng)柵格數(shù)據(jù)量達(dá)到3．2 TB（20萬(wàn)行×200萬(wàn)列×8字節(jié)）時(shí)，需存儲(chǔ)的壓縮單元數(shù)據(jù)量約為2 GB。隨著64位操作系統(tǒng)的普及，個(gè)人計(jì)算機(jī)內(nèi)存配置將逐步過(guò)渡到4—8 GB，本算法將具備處理TB級(jí)數(shù)據(jù)的應(yīng)用潛力。

5 結(jié) 論

本文針對(duì)常規(guī)種子點(diǎn)填充算法難以用于海量DEM數(shù)據(jù)淹沒(méi)分析問(wèn)題，提出采用分塊壓縮追蹤法進(jìn)行淹沒(méi)區(qū)生成。通過(guò)將DEM數(shù)據(jù)進(jìn)行條帶劃分，并將條帶中柵格行上連續(xù)多個(gè)淹沒(méi)單元進(jìn)行壓縮存儲(chǔ)，以降低數(shù)據(jù)量，最后采用壓縮單元邊界追蹤標(biāo)記方法提取淹沒(méi)范圍，從而實(shí)現(xiàn)了復(fù)雜地形條件下的淹沒(méi)區(qū)生成。與常見(jiàn)種子點(diǎn)填充方法相比，本文方法使用較小的內(nèi)存配置處理了較大的柵格數(shù)據(jù)量，同時(shí)避免了種子填充法中的大量遞歸判斷，提升了計(jì)算速度，具有一定的實(shí)用性。

［1］ LI Yun，F(xiàn)AN Ziwu，WU Shiqiang，et al．Numerical Simulation and 3D Visualization of Flood Propagation in Large－scale Detection Basins［J］．Journal of Hydraulic Engineering，2005（10）：1158－1164．（李云，范子武，吳時(shí)強(qiáng)，等．大型蓄滯洪區(qū)洪水演進(jìn)數(shù)值模擬與三維可視化技術(shù)［J］．水利學(xué)報(bào)，2005（10）：1158－1164．）

［2］ DONG Wenfeng，YUAN Yanbin，DU Yingze，et al．Simulation of 3D Terrain and Dynamic Simulation of Flood Routing Method［J］．Water Resources and Power，2001，19（3）：37－39．（董文鋒，袁艷斌，杜迎澤，等．流域三維地形仿真及洪水演進(jìn)動(dòng)態(tài)模擬［J］．水電能源科學(xué)，2001，19（3）：37－39．）

［3］ JIANG Jie，WU Lingda，XU Jiangbin．Digital Earth Based Flood Routing Model Visualization［J］．Computer Engineering and Applications，2009，45（36）：1－4．（蔣杰，吳玲達(dá)，徐江斌．基于數(shù)字地球的洪災(zāi)演進(jìn)模型表現(xiàn)［J］．計(jì)算機(jī)工程與應(yīng)用，2009，45（36）：1－4．）

［4］ WANG Zhili，GENG Yanfen，JIN Sheng．The Two－dimen－sional Flood Routing Simulation［J］．Chinese Journal of Computational Mechanics，2007，24（4）：533－538．（王志力，耿艷芬，金生．二維洪水演進(jìn)數(shù)值模擬［J］．計(jì)算力學(xué)學(xué)報(bào)，2007，24（4）：533－538．）

［5］ LI Daming，LIN Yi，XU Yanan，et al．Numerical Model of Flood Propagation［J］．Journal of Tianjin University，2009，42（1）：47－55．（李大鳴，林毅，徐亞男，等．河道、滯洪區(qū)洪水演進(jìn)數(shù)學(xué)模型［J］．天津大學(xué)學(xué)報(bào)，2009，42（1）：47－55．）

［6］ WEI Wenli，SHEN Yongming，SUN Guangcai，et al．Numerical Simulation of 2D Dam－break Flood Wave［J］．Journal of Hydraulic Engineering，2003（9）：43－47．（魏文禮，沈永明，孫廣才，等．二維潰壩洪水波演進(jìn)的數(shù)值模擬［J］．水利學(xué)報(bào)，2003（9）：43－47．）

［7］ ZHANG Xibing，LU Jinyou，F(xiàn)AN Beilin，et al．Dam Break Flood Routing and Drainage Process Modeling for the Tangjiashan Barrier Lake［J］．Yangtze River，2008，39（22）：21－22．（張細(xì)兵，盧金友，范北林，等．唐家山堰塞湖潰壩洪水演進(jìn)及下泄過(guò)程復(fù)演［J］．人民長(zhǎng)江，2008，39（22）：21－22．）

［8］ LI Guangchi，WANG Chuanhai．Universal Algorithm for River Basin Flood Routing Model［J］．Journal of Hohai University，2005，33（6）：624－628．（李光熾，王船海．流域洪水演進(jìn)模型通用算法研究［J］．河海大學(xué)學(xué)報(bào)：自然科學(xué)版，2005，33（6）：624－628．）

［9］ LIU Renyi，LIU Nan．A GIS Based Model for Calculating of Flood Area［J］．Acta Geographica Sinica，2001，56（1）：1－6．（劉仁義，劉南．基于GIS的復(fù)雜地形洪水淹沒(méi)區(qū)計(jì)算方法［J］．地理學(xué)報(bào)，2001，56（1）：1－6．）

［10］ LIU Renyi，LIU Nan．A New DEM－based Method for Flood Area Calculation and Damage Evaluation［J］．Journal of Image and Graphics，2001，6（2）：118－122．（劉仁義，劉南．一種基于數(shù)字高程模型DEM的淹沒(méi)區(qū)災(zāi)害評(píng)估方法［J］．中國(guó)圖像圖形學(xué)報(bào)，2001，6（2）：118－122．）

［11］ DING Zhixiong，LI Jiren，LI Lin．Method for Flood Submergence Analysis Based on GIS Grid Model［J］．Journal of Hydraulic Engineering，2004（6）：56－67．（丁志雄，李紀(jì)仁，李琳．基于GIS格網(wǎng)模型的洪水淹沒(méi)分析方法［J］．水利學(xué)報(bào)，2004（6）：56－67．）

［12］ SUN Hai，WANG Cheng．A“Ring”Method for Flood Submergence Based on DEM［J］．Geomatics and Information Science of Wuhan University，2009，34（8）：948－951．（孫海，王乘．利用DEM的“環(huán)形”洪水淹沒(méi)算法研究［J］．武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2009，34（8）：948－951．）

［13］ GUO Lihua，LONG Yi．Analysis of Flood Submerging Based on DEM［J］．Bulletin of Surveying and Mapping，2002（11）：25－30．（郭利華，龍毅．基于DEM的洪水淹沒(méi)分析［J］．測(cè)繪通報(bào)，2002（11）：25－30．）Binary Pictures［J］．Journal of Image and Graphics，2003，8（2）：198－202．（張修軍，郭霞，金心宇．帶標(biāo)記矯正的二值圖像連通域像素標(biāo)記算法［J］．中國(guó)圖像圖形學(xué)報(bào)，2003，8（2）：198－202．）

［19］ XIE Shunping，DU Jinkang，WANG Jiechen．Method for Tracing Run－length Outline to Implement Vectorization of Raster Graphics and Image Data［J］．Journal of Remote Sensing，2004，8（5）：465－470．（謝順平，都金康，王結(jié)臣．實(shí)現(xiàn)柵格圖形和圖像數(shù)據(jù)矢量化提取的游程輪廓追蹤法．遙感學(xué)報(bào)，2004，8（5）：465－470．）

［20］ XIE Shunping，DU Jinkang，WANG Lachun，et al．Approach of Vectorization for GIS Raster Data Based on Run－length Encoding System［J］．Acta Geodaetica et Cartographica Sinica，2004，33（4）：323－327．（謝順平，都金康，王臘春，等．基于游程編碼的GIS柵格數(shù)據(jù)矢量化方法［J］．測(cè)繪學(xué)報(bào)，2004，33（4）：323－327．）

（責(zé)任編輯：陳品馨）

［14］ JIANG Rengui，XIE Jiancang，LI Jianxun，et al．Analysis and Simulation of Flood Inundation Based on Digital Earth［J］．Engineering and Application，2011，47（13）：219－222．（姜仁貴，解建倉(cāng)，李建勛，等．基于數(shù)字地球的洪水淹沒(méi)分析及仿真研究［J］．計(jì)算機(jī)工程與應(yīng)用，2011，47（13）：219－222．）

［15］ MA Jianglin，ZHAO Zhongming，MENG Yu，et al．Connected Component Labelling for Massive Remote Sensing Classification Image［J］．Computer Engineering，2008，34（1）：262－264．（馬江林，趙忠明，孟瑜，等．海量遙感分類圖連通域標(biāo)記方法［J］．計(jì)算機(jī)工程，2008，34（1）：262－264．）

［16］ WANG Jing，ZHANG Yanning，LUO Jiancheng，et al．Improved Connected Component Labeling on High－resolution Remote Sensing Image［J］．Computer Engineering and Applications，2005，41（10）：37－39．（王晶，張艷寧，駱劍承，等．針對(duì)高分辨率遙感影像分割的改進(jìn)連通域標(biāo)記方法［J］．計(jì)算機(jī)工程與應(yīng)用，2005，41（10）：37－39．）

［17］ ZHAO Fei，ZHANG Lu，ZHANG Zhiyong，et al．A Hardware Acceleration Based Algorithm for Real－time Binary Image Connected－component Labeling［J］．Journal of Electronics＆Information Technology，2011，33（5）：1069－1075．（趙菲，張路，張志勇，等．基于硬件加速的實(shí)時(shí)二值圖像連通域標(biāo)記算法［J］，電子與信息學(xué)報(bào)，2011，33（5）：1069－1075．）

［18］ ZHANG Xiujun，GUO Xia，JIN Xinyu．The Pixel Labeled Algorithm with Label Rectified of Connecting Area in

A Quick Flood Inundation Algorithm Based on Massive DEM Data

SHEN Dingtao1，2，3，WANG Jiechen1，2，ZHANG Yu3，HUANG Jun3，HU Chengfang3
1．Jiangsu Provincial Key Laboratory of Geographic Information Science and Technology，Nanjing，210046，China；2．Department of Geographic Information Science，Nanjing University，Nanjing 210046，China；3．Spatial Information Technology Application Institute，Changjiang River Scientific Research Institute，Changjiang Water Resources Commission，Wuhan 430010，China

The conventional flood inundation methods based on DEM data usually have some disadvantages．For example，the seed filling method which is a popular flood inundation method can not get better result when the data amount is huge，and also it reflects low computational efficiency due to too many recursive operations．To resolve this problem，this paper proposes a quick flood inundation algorithm based on massive DEM data．It focuses on strip－divide method and real time raster compression storage technology．The paper makes the comparison between the common seed filling algorithm and strip－divide seed filling algorithm，the results show that this algorithm greatly improves the computational efficiency．At the same time，it resolves the problem of massive DEM data analysis．

massive DEM data；flood inundation；seed fill；compressed storage；digital elevation model

SHEN Dingtao（1983—），male，engineer，PhD candidate，majors in hydroinformatics and geographical information science．

P208

A

1001－1595（2014）06－0645－08

國(guó)家自然科學(xué)基金（41101408；41301435）；中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)項(xiàng)目（CKSF2013016／KJ；CKSF2014031／KJ）；水利部公益性行業(yè)科研專項(xiàng)（201201051）；教育部新世紀(jì)優(yōu)秀人才支持計(jì)劃（NCET－13－0280）

2013－03－20

沈定濤（1983—），男，工程師，博士生，主要研究方向?yàn)樗畔W(xué)和GIS。

SHEN Dingtao，WANG Jiechen，ZHANG Yu，et al．A Quick Flood Inundation Algorithm Based on Massive DEM Data［J］．Acta Geodaetica et Cartographica Sinica，2014，43（6）：645－652．（沈定濤，王結(jié)臣，張煜，等．一種面向海量數(shù)字高程模型數(shù)據(jù)的洪水淹沒(méi)區(qū)快速生成算法［J］．測(cè)繪學(xué)報(bào)，2014，43（6）：645－652．）

10．13485／j．cnki．11－2089．2014．0104

修回日期：2014－02－16