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

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

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

0 引言

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

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

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

1 算法設(shè)計(jì)

1.1 算法基本思路

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

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

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

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

圖1 算法基本思路鏈圖

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

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

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

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

圖2 算法示意圖

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

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

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

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

查找連通塊階段：

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

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

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

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

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

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

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

public class Point

{

public int x； //列號(hào)

public int y； //行號(hào)

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

}

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

public class DataCell

{

public int columnStartId；//起始列號(hào)

public int columnEndId；//終止列號(hào)

public int rowId；//行號(hào)

public int cellId；//壓縮單元編號(hào)

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

}

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

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

public class Lump

{

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

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

public int lumpId；//塊的編號(hào)

public int stripeId；//所在條帶編號(hào)

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

}

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

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

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

1.3 查找連通塊

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

圖4 查找連通塊流程圖

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

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

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

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

2 實(shí)驗(yàn)與分析

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

實(shí)驗(yàn)所用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)做測(cè)試。劃分條帶數(shù)量過(guò)多會(huì)造成迭代次數(shù)多、計(jì)算機(jī)運(yùn)行效率低。劃分條帶數(shù)量過(guò)少會(huì)使計(jì)算機(jī)一次讀入內(nèi)存的數(shù)據(jù)量大。經(jīng)過(guò)多次實(shí)驗(yàn)，使用取素?cái)?shù)的方法取得符合512×512大小內(nèi)最大的行數(shù)，以整除為優(yōu)先。設(shè)定北京地區(qū)洪水水位為250 m，四川省洪水水位為800 m，使用軟件ERDAS IMAGINE 2013展示淹沒(méi)區(qū)域。如圖5所示，右側(cè)疊加在原始DEM圖上的白色區(qū)域，即為算法輸出的淹沒(méi)區(qū)域DEM。

圖5 原始DEM和提取淹沒(méi)區(qū)域后DEM對(duì)比

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

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

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

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

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

3 結(jié)束語(yǔ)

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