王思雪,李英成,劉沛,耿中元,孫新博
(1.航空遙感技術(shù)國家測繪地理信息局重點(diǎn)實驗室,北京 100039;2.中測新圖(北京)遙感技術(shù)有限責(zé)任公司,北京 100039)
我國是一個地貌復(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)存占用的效果。
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)。
數(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
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。
查找連通塊的算法如圖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ū)域面積。
本文使用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處理效率對比
本文針對大區(qū)域DEM數(shù)據(jù)量較大而計算機(jī)性能有限的問題,使用塊碼壓縮和圖遍歷的方法進(jìn)行洪水淹沒范圍提取。實驗結(jié)果表明該算法能有效地提取出洪水淹沒范圍,對內(nèi)存的占用大大降低,DEM數(shù)據(jù)大小和個人計算機(jī)配置高低對算法的影響不大。本算法在輸出生成的淹沒區(qū)域DEM時,不需要預(yù)先生成與原DEM大小相同的空值柵格影像,只生成淹沒條帶間的區(qū)域,使得輸出的影像大小比原始影像小很多,大大節(jié)省了磁盤空間,提高了傳輸效率。