基于圖遍歷的計算DEM數(shù)據(jù)洪水淹沒范圍的算法

王思雪，李英成，劉沛，耿中元，孫新博

(1.航空遙感技術(shù)國家測繪地理信息局重點(diǎn)實驗室，北京 100039；2.中測新圖(北京)遙感技術(shù)有限責(zé)任公司，北京 100039)

0 引言

我國是一個地貌復(fù)雜多樣的國家，山川河流縱橫交錯，洪水災(zāi)害常伴隨強(qiáng)降雨發(fā)生，尤其是南方近河流地區(qū)發(fā)生比較頻繁。無論是災(zāi)害預(yù)警還是應(yīng)急救災(zāi)，及時有效地獲取洪水淹沒地區(qū)范圍，對保障人民生命和財產(chǎn)安全有重大意義。近年來，國內(nèi)外已廣泛利用GIS、RS、LiDAR、SAR等手段對洪水淹沒范圍進(jìn)行預(yù)測和評估[1-4]。在以通過各種手段獲取并經(jīng)過處理得到的DEM為柵格影像數(shù)據(jù)進(jìn)行洪水淹沒分析時，有靜態(tài)水平面[5]和建立水動力學(xué)模型[6]2種方法，根據(jù)靜態(tài)水平面計算淹沒范圍相比于基于水文、水動力學(xué)構(gòu)建的洪水演進(jìn)模型的方法來說，精確度不如后者，但是實用性強(qiáng)，在防洪抗洪工作中有較大的實際意義。

靜態(tài)水平面法按給定水位的“水平面”計算淹沒范圍，分為無源式淹沒和有源式淹沒2種方式。無源式水淹分析算法通過將DEM中高程點(diǎn)低于設(shè)定水位的視為淹沒點(diǎn)。有源式淹沒則需考慮水淹區(qū)域的連通性，將高程低于設(shè)定水位并與指定位置有連通的DEM格網(wǎng)點(diǎn)視為淹沒點(diǎn)，我國近河流地形條件的復(fù)雜性決定有源淹沒更具有適用性。劉仁義等[7-8]采用種子蔓延法通過給定的種子點(diǎn)(即指定淹沒源)按照四鄰域或八鄰域進(jìn)行擴(kuò)散找到鄰接的低于給定水位高度的連通點(diǎn)，所有連通點(diǎn)的集合就是淹沒區(qū)范圍。楊啟貴等[9]采用分塊種子蔓延算法對DEM數(shù)據(jù)進(jìn)行將研究區(qū)域分成預(yù)先設(shè)定大小的塊，對所有的正方形塊內(nèi)使用種子蔓延法計算，對相鄰塊間判斷連通性，最后通過包含初始種子點(diǎn)的淹沒區(qū)以及其所對應(yīng)的關(guān)聯(lián)表來提取出整個淹沒區(qū)。沈定濤等[10]采用分塊壓縮追蹤法將研究區(qū)域劃分為條帶，將條帶中每行中高程值低于設(shè)定水位的相鄰網(wǎng)格單元進(jìn)行游程編碼壓縮，在較大程度上降低了柵格數(shù)據(jù)量，由種子點(diǎn)壓縮單元開始，以壓縮單元兩側(cè)作為追蹤邊，在不同條帶上進(jìn)行追蹤，最后提取柵格場上被追蹤過的所有壓縮單元即可構(gòu)建淹沒區(qū)。

在分析了幾種洪水淹沒算法后，考慮到處理系統(tǒng)可靠性和算法安全性，對算法進(jìn)行了如下設(shè)計：在DEM數(shù)據(jù)讀入時進(jìn)行分條帶讀入計算機(jī)內(nèi)存，在數(shù)據(jù)壓縮方式上采用了塊碼壓縮將潛在連通塊標(biāo)記并將條帶序列化為二進(jìn)制流存入磁盤，在搜索連通區(qū)域時，將遍歷到的條帶從磁盤反序列化釋放到內(nèi)存，利用求解圖的連通性問題的思路，采用圖遍歷的廣度優(yōu)先搜索法進(jìn)行連通區(qū)域的遞歸查找。程序上通過設(shè)定的洪水水位、選定種子點(diǎn)計算某個區(qū)域內(nèi)的淹沒范圍。選取北京市、四川省的DEM數(shù)據(jù)進(jìn)行實驗，與種子蔓延算法和分塊種子蔓延法進(jìn)行比較，最后輸出淹沒區(qū)所在條帶的DEM，實現(xiàn)了快速計算和減小內(nèi)存占用的效果。

1 算法設(shè)計

1.1 算法基本思路

DEM是以柵格形式存儲的地形數(shù)據(jù)，其表示形式包括規(guī)則矩陣格網(wǎng)和不規(guī)則三角網(wǎng)，本文研究的DEM是規(guī)則格網(wǎng)形式的，由一系列地面點(diǎn)及其所具有的高程組成，在計算機(jī)高級語言中，它就是一個二維數(shù)組或者數(shù)學(xué)上的一個二維矩陣{Zij}[11-12]。

當(dāng)大范圍DEM的數(shù)據(jù)量較大時，一次讀入和輸出通常會造成對內(nèi)存很大的占用，容易造成系統(tǒng)崩潰?？紤]到系統(tǒng)安全性，采取對DEM數(shù)據(jù)進(jìn)行劃分，分次讀入和分次輸出。劃分方式可以按行列劃分為規(guī)則矩形，可以按行劃分為條帶。考慮到本文需要將DEM數(shù)據(jù)分批次由磁盤讀入內(nèi)存，再經(jīng)過將潛在淹沒柵格存入磁盤待調(diào)用，故使用將DEM劃分為條帶的方法，使其能更方便地存儲和調(diào)用，且按條帶比按規(guī)則矩形更容易實現(xiàn)辨別水流橫向上連通關(guān)系，減少循環(huán)迭代的判斷次數(shù)。

算法主要分為數(shù)據(jù)壓縮存儲階段和查找連通區(qū)階段。矩陣格網(wǎng)DEM具有可進(jìn)行壓縮存儲的特點(diǎn)，基于這個特點(diǎn)，本文在數(shù)據(jù)壓縮存儲階段，對每條帶內(nèi)不同行和列的低于某個高程值的相鄰像元(具有潛在連通性的淹沒柵格)形成的面實體以塊碼壓縮存儲，即采用方形區(qū)域作為記錄單元，每個記錄單元包括相鄰的若干柵格，數(shù)據(jù)結(jié)構(gòu)由初始位置(行、列號)和半徑，再加上記錄單位的代碼組成[13]。塊碼在合并、插入、檢查延伸性、計算面積等操作時有明顯的優(yōu)越性。

圖遍歷方法是數(shù)據(jù)結(jié)構(gòu)中求解圖的連通性問題、拓?fù)渑判蚝颓箨P(guān)鍵路徑等算法的基礎(chǔ)。圖的遍歷方法目前有深度優(yōu)先搜索法和廣度優(yōu)先搜索法2種算法[14]。廣度優(yōu)先搜索是從圖的某個頂點(diǎn)出發(fā)，訪問該頂點(diǎn)后依次訪問其相鄰的未被訪問的點(diǎn)，然后分別從這些鄰接點(diǎn)出發(fā)依次訪問它們的鄰接點(diǎn)，直至圖中所有已被訪問的頂點(diǎn)的鄰接點(diǎn)都被訪問到。若此時圖中還有頂點(diǎn)未被訪問，則另選圖中一個未被訪問的低點(diǎn)作為起點(diǎn)，重復(fù)上述過程，直到圖中所有頂點(diǎn)都被訪問到為止。深度優(yōu)先搜索與廣度優(yōu)先搜索類似，只是在遍歷順序上是以頂點(diǎn)訪問后的鄰接點(diǎn)出發(fā)，使用深度優(yōu)先遍歷圖。圖遍歷的實質(zhì)就是對每個頂點(diǎn)查找鄰接點(diǎn)的一個遞歸的過程。本文在查找連通區(qū)階段使用廣度優(yōu)先搜索法由種子點(diǎn)所在塊(簡稱種子塊)作為頂點(diǎn)出發(fā)找到鄰接塊，再由訪問后的塊找到鄰接的未訪問到的塊，所有被訪問的塊構(gòu)成的集合就是有源淹沒的淹沒范圍。本文算法的基本思路如圖1所示。

圖1 算法基本思路鏈圖

由于水流本身蜿蜒曲折的特點(diǎn)，單個條帶肯定會被多次重復(fù)搜索淹沒柵格，故本算法首先由種子點(diǎn)所在條帶(簡稱種子條帶)開始，確定每個條帶內(nèi)所有具有潛在連通性的淹沒柵格，壓縮成塊，塊作為判斷條帶間連通性的基本形式。條帶以塊為組成，塊由壓縮單元組成。

本算法定義的條帶間連通條件：塊頂端(塊頭)存在行號等于所在條帶首行行號的壓縮單元，或者塊末端(塊尾)存在行號等于所在條帶末行行號的壓縮單元。

本算法定義的條帶上下追蹤條件：條帶內(nèi)存在行數(shù)等于條帶柵格行數(shù)的塊，可以向上一條帶或者下一條帶繼續(xù)追蹤滿足條帶間連通條件的塊。

將滿足條帶上下追蹤條件的條帶由內(nèi)存寫入磁盤，再由種子塊開始在條帶間循環(huán)遍歷，找到不同條帶間連通的淹沒壓縮塊，從磁盤將數(shù)據(jù)調(diào)出來，標(biāo)記條帶內(nèi)滿足追蹤條件的塊頭和/或塊尾并計算塊內(nèi)柵格個數(shù)，經(jīng)過多次循環(huán)迭代最終計算得到整個淹沒連通區(qū)的范圍大小。算法分為數(shù)據(jù)壓縮存儲和查找連通塊2個階段執(zhí)行，如圖2所示。

圖2 算法示意圖

數(shù)據(jù)壓縮階段：

①讀取種子點(diǎn)的條帶1，對種子點(diǎn)所在柵格進(jìn)行數(shù)據(jù)壓縮得到塊0。如果塊0符合連通條件，即種子點(diǎn)選取合適，在該條帶從左到右遍歷找到符合連通條件的塊1、塊2，然后該條帶存入磁盤。

②對條帶1上條帶的柵格執(zhí)行數(shù)據(jù)壓縮和存盤，并依次向上方條帶執(zhí)行，直到條帶的頂柵格行中無壓縮單元。

③對條帶1下一條帶執(zhí)行數(shù)據(jù)壓縮和存盤，并依次向下方條帶執(zhí)行，直到條帶的底柵格行中無壓縮單元。

查找連通塊階段：

①將條帶1從磁盤取到內(nèi)存，由塊0的頂部壓縮單元向上一條帶搜索(條帶0取到內(nèi)存)，找到條帶0中塊1、塊2，標(biāo)記塊底壓縮單元，并向下一條帶搜索，找到條帶1中滿足連通條件且未標(biāo)記的塊1，標(biāo)記其塊頂壓縮單元。由條帶1中塊1向上一條帶搜索，無未標(biāo)記連通塊，結(jié)束查找。

②由條帶1中塊0的底部2個壓縮單元向下一條帶搜索(條帶2取到內(nèi)存)，找到條帶2中塊0，標(biāo)記塊頂壓縮單元，因其塊底壓縮單元不在條帶底部，不滿足條帶上下追蹤條件，結(jié)束查找。

③對條帶0～2中所有滿足條帶上下追蹤條件的被標(biāo)記的塊，通過讀取塊結(jié)構(gòu)內(nèi)柵格總數(shù)計算得到淹沒面積，并記錄有效條帶號(0～2)。

1.2 數(shù)據(jù)壓縮存儲

數(shù)據(jù)壓縮存儲的算法流程圖如圖3所示。由種子條帶開始，讀取種子條帶的DEM柵格數(shù)據(jù)，開始遍歷所有淹沒柵格，并將每行中連續(xù)淹沒柵格記錄成一個壓縮單元DataCell，壓縮單元內(nèi)記錄起始列號、終止列號、行號、編號。

(注：圖3默認(rèn)選取的種子點(diǎn)是淹沒點(diǎn)，種子塊滿足條帶間連通條件。)圖3 數(shù)據(jù)壓縮存儲流程圖

柵格點(diǎn)結(jié)構(gòu)如下：

public class Point

{

public int x； //列號

public int y； //行號

public double elevation；//像素值(高程值)

}

壓縮單元結(jié)構(gòu)如下：

public class DataCell

{

public int columnStartId；//起始列號

public int columnEndId；//終止列號

public int rowId；//行號

public int cellId；//壓縮單元編號

public bool isFlood；//是否已標(biāo)記

}

以種子點(diǎn)所在的壓縮單元為頂點(diǎn)，用深度優(yōu)先搜索法將不同行之間的相鄰壓縮單元壓縮成塊，塊內(nèi)以鏈表形式存儲塊頭和塊尾的壓縮單元，并記錄塊內(nèi)柵格行數(shù)、塊號、條帶號、柵格總數(shù)。

塊結(jié)構(gòu)如下：

public class Lump

{

public Queue m_CompressCell；//包含塊頭和塊尾的壓縮單元的隊列

public int rowsNum；//塊內(nèi)行數(shù)

public int lumpId；//塊的編號

public int stripeId；//所在條帶編號

public int pointsNum；//塊內(nèi)柵格總數(shù)

}

對條帶內(nèi)滿足條帶間連通條件的塊按從左到右、從上到下的順序進(jìn)行編號，將不滿足條帶間連通條件的塊舍棄，條帶以鏈表形式存儲已編號的塊。

條帶需滿足條帶上下追蹤條件。具體地，如果種子條帶的頂柵格行有被編號的塊，向種子條帶的上一條帶執(zhí)行同樣方式的壓縮存儲，并以此一直向上執(zhí)行，直到向上的條帶內(nèi)沒有滿足條帶上下追蹤條件的塊，表示已不能再繼續(xù)向上追蹤，向上搜索結(jié)束，記錄最小截止條帶號TStripeId。

同樣地，如果種子條帶的底柵格行有被編號的塊，向種子條帶的下一條帶執(zhí)行壓縮存儲，并以此一直向下執(zhí)行，直到向下的條帶內(nèi)沒有滿足條帶上下追蹤條件的塊，結(jié)束向下搜索，記錄最大截止條帶號EStripeId。

1.3 查找連通塊

查找連通塊的算法如圖4所示。為了在遍歷過程中標(biāo)記塊頭和塊尾壓縮單元是否是連通的，需設(shè)置布爾對象isFlood，初值設(shè)為“false”，當(dāng)訪問到該壓縮單元時值為“true”，表示是已標(biāo)記的連通的壓縮單元，避免重復(fù)運(yùn)算。若塊頭或者塊尾存在已標(biāo)記的壓縮單元，代表該塊為淹沒塊，若已標(biāo)記塊的行數(shù)等于所在條帶行數(shù)，則表示該塊滿足條帶上下追蹤條件。

圖4 查找連通塊流程圖

根據(jù)條帶號從硬盤將種子條帶反序列化到內(nèi)存，在上下截止條帶號區(qū)間內(nèi)，由種子塊開始，取種子塊頭壓縮單元(多個DataCell)加入隊列，取種子條帶頂柵格行上所有壓縮單元和其上一條帶(反序列化從磁盤讀取出)底柵格行上的所有壓縮單元放入搜索隊列，用廣度優(yōu)先搜索找到與種子塊(已標(biāo)記)相鄰接的壓縮單元，進(jìn)而找到對應(yīng)的連通塊。

同樣地，取種子塊尾壓縮單元(多個DataCell)加入隊列，取種子條帶底柵格行上所有壓縮單元和其下一條帶(反序列化從磁盤讀取出)的所有底柵格行上的壓縮單元，通過條帶之間鄰接的壓縮單元到對應(yīng)連通塊。

從種子條帶開始，將所有相鄰2個條帶間的鄰接兩行(頂/底)滿足追蹤條件的塊頭和塊尾壓縮單元放入搜索隊列。同時將滿足追蹤條件的該條帶帶號存入條帶號隊列，用于存儲有效條帶號。最后輸出和顯示有效條帶間的淹沒區(qū)域，減小輸出DEM的大小。

如果已標(biāo)記過一端的塊滿足條帶上下追蹤條件，將塊的另一端未標(biāo)記壓縮單元也標(biāo)記，遞歸使用廣度優(yōu)先搜索法對[TStripeId，EStripeId]條帶間滿足條帶間連通條件的塊進(jìn)行搜索，最終得到連通區(qū)間有效截止條帶號，并通過統(tǒng)計連通塊的面積得到整個淹沒區(qū)域面積。

2 實驗與分析

本文使用VS2010.NET3.5平臺，C#編程語言開發(fā)，為驗證算法運(yùn)行效率和內(nèi)存占用，使用運(yùn)行環(huán)境：Windows 64位，CPU為Intel 酷睿i5 5200u處理器，4 GB內(nèi)存空間的筆記本做測試。

實驗所用DEM數(shù)據(jù)選取北京地區(qū)DEM(數(shù)據(jù)包含的網(wǎng)格數(shù)8 621列×11 102行，數(shù)據(jù)文件大小100 MB)和四川省DEM(數(shù)據(jù)包含的網(wǎng)格數(shù)45 210列×25 201行，數(shù)據(jù)文件大小2.57 GB)做測試。劃分條帶數(shù)量過多會造成迭代次數(shù)多、計算機(jī)運(yùn)行效率低。劃分條帶數(shù)量過少會使計算機(jī)一次讀入內(nèi)存的數(shù)據(jù)量大。經(jīng)過多次實驗，使用取素數(shù)的方法取得符合512×512大小內(nèi)最大的行數(shù)，以整除為優(yōu)先。設(shè)定北京地區(qū)洪水水位為250 m，四川省洪水水位為800 m，使用軟件ERDAS IMAGINE 2013展示淹沒區(qū)域。如圖5所示，右側(cè)疊加在原始DEM圖上的白色區(qū)域，即為算法輸出的淹沒區(qū)域DEM。

圖5 原始DEM和提取淹沒區(qū)域后DEM對比

本次水淹分析的算法使用北京地區(qū)DEM和四川省DEM作為數(shù)據(jù)源，對種子蔓延法、分塊種子蔓延算法與本文算法進(jìn)行比較(其中運(yùn)行平均時間包括計算淹沒水面面積和輸出淹沒區(qū)域DEM)。

由表1可見，對柵格數(shù)據(jù)壓縮成壓縮單元明顯提高了運(yùn)行效率，內(nèi)存占用最小，且采用本文只對有效影像條帶區(qū)域輸出的方法(圖5)，輸出的淹沒區(qū)域DEM數(shù)據(jù)相比于原始DEM在數(shù)據(jù)量上很大程度地縮小，有利于影像快速輸出和傳輸。

表1 3種算法運(yùn)算效率和內(nèi)存占用對比

由表2可見，本算法在對于不同數(shù)據(jù)大小的DEM數(shù)據(jù)進(jìn)行洪水淹沒范圍提取中兼顧了運(yùn)行效率和防止內(nèi)存溢出兩方面的優(yōu)勢，內(nèi)存占用幾乎不受原始DEM數(shù)據(jù)大小限制。

表2 算法對不同數(shù)據(jù)大小的DEM處理效率對比

3 結(jié)束語

本文針對大區(qū)域DEM數(shù)據(jù)量較大而計算機(jī)性能有限的問題，使用塊碼壓縮和圖遍歷的方法進(jìn)行洪水淹沒范圍提取。實驗結(jié)果表明該算法能有效地提取出洪水淹沒范圍，對內(nèi)存的占用大大降低，DEM數(shù)據(jù)大小和個人計算機(jī)配置高低對算法的影響不大。本算法在輸出生成的淹沒區(qū)域DEM時，不需要預(yù)先生成與原DEM大小相同的空值柵格影像，只生成淹沒條帶間的區(qū)域，使得輸出的影像大小比原始影像小很多，大大節(jié)省了磁盤空間，提高了傳輸效率。