張連翀， 李國(guó)慶， 于文洋， 冉全

(1．中國(guó)科學(xué)院遙感與數(shù)字地球研究所數(shù)字地球?qū)嶒?yàn)室，北京 100094； 2．中國(guó)科學(xué)院大學(xué)，北京 100049；3．海南省地球觀測(cè)重點(diǎn)實(shí)驗(yàn)室，海南 572023)

基于水平集的洪澇淹沒(méi)范圍時(shí)空模擬方法

張連翀1,2， 李國(guó)慶1,3， 于文洋1,3， 冉全1,2

(1．中國(guó)科學(xué)院遙感與數(shù)字地球研究所數(shù)字地球?qū)嶒?yàn)室，北京 100094； 2．中國(guó)科學(xué)院大學(xué)，北京 100049；3．海南省地球觀測(cè)重點(diǎn)實(shí)驗(yàn)室，海南 572023)

遙感技術(shù)能及時(shí)獲取洪水空間分布特征信息，已成為洪澇災(zāi)害監(jiān)測(cè)與損失評(píng)估的重要依據(jù)。然而受天氣和環(huán)境等因素影響，不能全天時(shí)接收遙感影像導(dǎo)致部分?jǐn)?shù)據(jù)缺失，無(wú)法提供動(dòng)態(tài)連續(xù)的洪澇淹沒(méi)過(guò)程資料。以2013年汛期黑龍江流域八岔段潰口淹沒(méi)區(qū)為例，基于多時(shí)相GF-1衛(wèi)星晴空遙感影像提取的洪澇淹沒(méi)范圍信息，將洪澇淹沒(méi)過(guò)程轉(zhuǎn)化為水平集函數(shù)的偏微分方程數(shù)值求解問(wèn)題，利用空間迎風(fēng)差分格式和時(shí)間歐拉差分格式模擬了從8月24日到10月8日洪水漲退過(guò)程的逐日淹沒(méi)范圍。精度評(píng)價(jià)結(jié)果表明，洪澇淹沒(méi)范圍的模擬結(jié)果與同時(shí)期遙感影像提取結(jié)果的Kappa系數(shù)分別為0.921 2和0.893 2； 與同時(shí)期洪澇淹沒(méi)范圍的統(tǒng)計(jì)數(shù)據(jù)相比略偏低，但相對(duì)誤差都小于10%。該方法的模擬結(jié)果與影像提取結(jié)果和實(shí)際統(tǒng)計(jì)數(shù)據(jù)都具有較好的時(shí)空一致性，為不依賴(lài)先驗(yàn)資料的洪澇災(zāi)害應(yīng)急響應(yīng)決策提供了科學(xué)依據(jù)。

洪澇淹沒(méi)范圍； 時(shí)空模擬； 水平集； 遙感； 黑龍江洪水

0 引言

準(zhǔn)確、科學(xué)地確定洪澇淹沒(méi)范圍，是防汛救災(zāi)決策和災(zāi)害損失評(píng)估的核心任務(wù)之一。水文學(xué)上多采用依賴(lài)先驗(yàn)資料的洪水演進(jìn)模型描述洪水運(yùn)動(dòng)過(guò)程，主要包括基于水文資料的洪水動(dòng)力學(xué)模型[1]和基于地形資料的格網(wǎng)模型[2]。前者通過(guò)求解水動(dòng)力學(xué)方程精確模擬水位、流量、流速及其隨時(shí)間的變化過(guò)程，但是輸入?yún)?shù)的地區(qū)差異性增加了模型結(jié)果和精度的不確定性； 后者通過(guò)種子蔓延算法[3]求取滿足條件的格網(wǎng)集合，但是遞歸判斷過(guò)多導(dǎo)致計(jì)算效率較低。而且，洪澇災(zāi)害具有突發(fā)性特點(diǎn)，很多受災(zāi)地區(qū)無(wú)法提供足夠的、實(shí)時(shí)的觀測(cè)資料，已成為防洪應(yīng)急指揮調(diào)度和搶險(xiǎn)救災(zāi)的薄弱環(huán)節(jié)。

遙感數(shù)據(jù)作為洪澇災(zāi)害監(jiān)測(cè)與損失評(píng)估的重要信息來(lái)源[4]，既能夠直接用于提取洪澇淹沒(méi)范圍信息[5-7]，也能夠通過(guò)提取土地覆被、不透水面積比等下墊面特征信息推求水文模型參數(shù)[8]。然而受天氣和環(huán)境等因素影響，不能全天時(shí)接收遙感影像導(dǎo)致了部分?jǐn)?shù)據(jù)的缺失，無(wú)法提供動(dòng)態(tài)連續(xù)的洪澇淹沒(méi)過(guò)程資料[9]。水平集方法(level set method)作為一種通過(guò)極小化能量泛函追蹤界面移動(dòng)的數(shù)值技術(shù)，將洪澇淹沒(méi)過(guò)程轉(zhuǎn)化為淹沒(méi)范圍曲面間的拓?fù)渥冃?，從空間和時(shí)間維度進(jìn)行數(shù)值求解，以獲得具有時(shí)空一致性的洪澇淹沒(méi)范圍模擬結(jié)果。

1 研究區(qū)概況與數(shù)據(jù)源

1.1 研究區(qū)概況

2013年夏季汛期，黑龍江流域遭遇歷史罕見(jiàn)特大洪水。本文以淹沒(méi)范圍最大的黑龍江省同江市八岔段潰口淹沒(méi)區(qū)為研究區(qū)。該地區(qū)位于黑龍江干流下游南岸，地勢(shì)低平，黑龍江、松花江和烏蘇里江等多條河流匯合而來(lái)的洪水流經(jīng)該區(qū)域速度變緩，洪峰持續(xù)時(shí)間長(zhǎng)。2013年8月23日該地區(qū)出現(xiàn)潰口，造成當(dāng)?shù)? 000余戶房屋受損，農(nóng)田成災(zāi)面積達(dá)7.6萬(wàn)hm2。

1.2 數(shù)據(jù)預(yù)處理

選擇2013年(本文影像數(shù)據(jù)均為2013年，下文省略)汛期黑龍江流域的6景(7月10日、8月27日、8月28日、9月9日、9月18日、10月8日)GF-1衛(wèi)星多光譜晴空遙感影像，空間分辨率16 m。遙感影像經(jīng)過(guò)相對(duì)輻射校正、幾何糾正等預(yù)處理后，首先計(jì)算歸一化差異水體指數(shù)(normalized difference water index,NDWI)以最大限度地突出水體與陸地、植被之間的差異； 然后采用最大類(lèi)間方差法[10]確定分割閾值提取水體邊界； 再對(duì)研究區(qū)背景信息和洪水淹沒(méi)信息進(jìn)行二值化處理，最終得到不同時(shí)刻的洪澇淹沒(méi)范圍。

2 研究方法

水平集方法首先由Osher和Sethian提出，是處理運(yùn)動(dòng)曲面隨時(shí)間拓?fù)渥兓挠行в?jì)算工具[11]。該方法將連續(xù)函數(shù)φ(x,y,t)∶R3→R描述為閉合演化曲線C(p,t)： 0≤p≤1在t時(shí)刻的隱式表達(dá)，即t時(shí)刻曲線C(p,t)對(duì)應(yīng)于φ(x,y,t)的零水平集。本文設(shè)t1時(shí)刻的洪澇淹沒(méi)范圍為源曲面φ1，t2時(shí)刻的洪澇淹沒(méi)范圍為靶曲面φ2(t1

在具體實(shí)現(xiàn)時(shí)，首先選取一個(gè)動(dòng)態(tài)演化的水平集函數(shù)φ(x,y,t)，并分別定義影像上的點(diǎn)X(i,j)到φ1和φ2的符號(hào)距離函數(shù)(signed distance function，SDF)d1和d2作為初始水平集函數(shù)，即

φ(x,y,t)=±d,

(1)

式中d是影像上的點(diǎn)X(i,j)到曲面網(wǎng)格的最短歐式距離。當(dāng)X(i,j)在曲面外時(shí)取正值，當(dāng)X(i,j)在曲面內(nèi)時(shí)取負(fù)值，當(dāng)X(i,j)在曲面上時(shí)值為0，該選擇的優(yōu)點(diǎn)是|φ|=1，有利于數(shù)值計(jì)算的穩(wěn)定性。

采用初始水平集函數(shù)描述φ1后，根據(jù)演化速度場(chǎng)求解Hamilton-Jacobi方程，實(shí)現(xiàn)水平集函數(shù)的動(dòng)態(tài)演化，即

(2)

Hamilton-Jacobi方程可以采用差分方法進(jìn)行數(shù)值求解。為了保證求解的精度，避免數(shù)值耗散，可以選取空間迎風(fēng)差分格式和時(shí)間歐拉差分格式進(jìn)行離散，即

(3)

(4)

(5)

(6)

(7)

(8)

(9)

利用隱式表達(dá)的水平集函數(shù)進(jìn)行曲面演化時(shí)，水平集函數(shù)的零等值面即為該時(shí)間點(diǎn)的洪澇淹沒(méi)范圍演化形狀。

3 實(shí)驗(yàn)結(jié)果與分析

3.1 漲水過(guò)程模擬

由于發(fā)生潰口前未獲得研究區(qū)晴空影像數(shù)據(jù)，所以選擇7月10日晴空影像提取河道本底水體信息，并將8月27日、9月9日洪澇淹沒(méi)范圍作為輸入數(shù)據(jù)，模擬該時(shí)間段內(nèi)洪澇淹沒(méi)過(guò)程(圖1)。

圖1-1 洪澇淹沒(méi)范圍時(shí)空模擬結(jié)果(漲水過(guò)程)

圖1-2 洪澇淹沒(méi)范圍時(shí)空模擬結(jié)果(漲水過(guò)程)

結(jié)果表明，自8月23日發(fā)生潰口后，河堤沿岸地區(qū)被迅速淹沒(méi)，洪水向東北部蔓延直至撫遠(yuǎn)縣，淹沒(méi)范圍迅速增大； 8月27日后洪水淹沒(méi)速度逐漸變緩，淹沒(méi)范圍向潰口處西南位置推進(jìn)，至9月9日達(dá)到最大面積。

3.2 退水過(guò)程模擬

選擇9月9日、10月8日洪澇淹沒(méi)范圍作為輸入數(shù)據(jù)，模擬該時(shí)間段內(nèi)洪澇淹沒(méi)過(guò)程。結(jié)果表明，自9月10日后淹沒(méi)區(qū)開(kāi)始退水，積水由西南部回流至黑龍江，至10月8日積水基本排空(圖2)。

圖2-1 洪澇淹沒(méi)范圍時(shí)空模擬結(jié)果(退水過(guò)程)

圖2-2 洪澇淹沒(méi)范圍時(shí)空模擬結(jié)果(退水過(guò)程)

3.3 模擬結(jié)果精度評(píng)價(jià)

通過(guò)求解水平集函數(shù)的偏微分方程，最終獲得八岔段潰口淹沒(méi)區(qū)從8月24日到10月8日的逐日洪澇淹沒(méi)范圍模擬結(jié)果。為了便于評(píng)價(jià)該方法的準(zhǔn)確性，分別利用2013年汛期不同時(shí)間的遙感影像提取結(jié)果和水利部統(tǒng)計(jì)數(shù)據(jù)，對(duì)洪澇淹沒(méi)范圍的模擬結(jié)果進(jìn)行了精度評(píng)價(jià)。

將8月28日和9月18日的模擬結(jié)果分別與基于同時(shí)期遙感影像的提取結(jié)果建立混淆矩陣，結(jié)果如表1所示。

表1 模擬結(jié)果與提取結(jié)果的精度統(tǒng)計(jì)

從表1中可以看出，2組數(shù)據(jù)的Kappa系數(shù)分別為0.921 2和0.893 2，說(shuō)明模擬結(jié)果與影像實(shí)際提取結(jié)果具有較好的時(shí)空一致性； 模擬結(jié)果與輸入數(shù)據(jù)的時(shí)間距離越遠(yuǎn)，模擬結(jié)果的正確率越低，但仍都保持在90%左右，模擬結(jié)果較為理想。

另外，依據(jù)水利部水利信息中心統(tǒng)計(jì)資料[12]，也對(duì)模擬結(jié)果進(jìn)行相對(duì)誤差分析(表2)。

表2 模擬結(jié)果與統(tǒng)計(jì)結(jié)果的相對(duì)誤差分析

從表2中可以看出，遙感影像提取的最大洪澇淹沒(méi)面積為773 km2，與統(tǒng)計(jì)數(shù)據(jù)的相對(duì)誤差僅為1.18%； 受無(wú)潰口發(fā)生時(shí)遙感影像的限制，同期模擬結(jié)果都略高于統(tǒng)計(jì)數(shù)據(jù)，但相對(duì)誤差均小于10%； 8月25—27日2 d淹沒(méi)區(qū)面積從293 km2增加到647 km2，面積增大速率為177 km2/d，與統(tǒng)計(jì)數(shù)據(jù)的面積增大速率(168.5 km2/d)相比較為一致。結(jié)果表明，本文方法能夠有效模擬洪澇淹沒(méi)區(qū)域的整體汛情狀況和洪水發(fā)展變化趨勢(shì)。

4 結(jié)論

1)基于水平集的洪澇淹沒(méi)范圍時(shí)空模擬方法，將片段化的洪澇淹沒(méi)范圍數(shù)據(jù)擴(kuò)展為連續(xù)動(dòng)態(tài)的淹沒(méi)過(guò)程信息，能夠?qū)崿F(xiàn)洪澇淹沒(méi)范圍的快速模擬。研究表明，模擬結(jié)果與影像提取結(jié)果具有較好的時(shí)空一致性，為洪澇淹沒(méi)歷時(shí)提取、淹沒(méi)程度和災(zāi)害損失評(píng)估等后續(xù)工作提供了科學(xué)依據(jù)。相比于依賴(lài)先驗(yàn)資料的洪水演進(jìn)模型，該方法具有明顯的通用性優(yōu)勢(shì)，提高了高空間分辨率遙感影像在洪澇災(zāi)害監(jiān)測(cè)中的應(yīng)用前景。

2)遙感數(shù)據(jù)中包含了豐富的洪澇淹沒(méi)時(shí)空信息。受高空間分辨率遙感影像獲取時(shí)間的限制，模擬結(jié)果與輸入影像的時(shí)間距離越遠(yuǎn)，模擬結(jié)果精度越低。后續(xù)工作中擬主要考慮融合多源、多類(lèi)型、多分辨率的遙感影像，高空間分辨率和高時(shí)間分辨率的數(shù)據(jù)相結(jié)合以提高水體邊界提取精度，進(jìn)一步改進(jìn)本文方法及其模擬結(jié)果。

[1] 俞云利,賴(lài)錫軍.二維平面非恒定流數(shù)學(xué)模型的遙感水位數(shù)據(jù)同化[J].水科學(xué)進(jìn)展,2008,19(2):224-231. Yu Y L,Lai X J.2D horizontal unsteady flow model for assimilating remote sensing water levels[J].Advances in Water Science,2008,19(2):224-231.

[2] 沈定濤,王結(jié)臣,張煜,等.一種面向海量數(shù)字高程模型數(shù)據(jù)的洪水淹沒(méi)區(qū)快速生成算法[J].測(cè)繪學(xué)報(bào),2014,43(6):645-652. Shen D T,Wang J C,Zhang Y,et al.A quick flood inundation algorithm based on massive DEM data[J].Acta Geodaetica et Cartographica Sinica,2014,43(6):645-652.

[3] 楊軍,賈鵬,周廷剛,等.基于DEM的洪水淹沒(méi)模擬分析及虛擬現(xiàn)實(shí)表達(dá)[J].西南大學(xué)學(xué)報(bào):自然科學(xué)版,2011,33(10):143-148. Yang J,Jia P,Zhou T G,et al.DEM-based simulation analysis and virtual reality expression of flood submergence[J].Journal of Southwest University:Natural Science Edition,2011,33(10):143-148.

[4] 孟令奎,郭善昕,李爽.遙感影像水體提取與洪水監(jiān)測(cè)應(yīng)用綜述[J].水利信息化,2012(3):18-25. Meng L K,Guo S X,Li S.Summary on extraction of water body from remote sensing image and flood monitoring[J].Water Resources Informatization,2012(3):18-25.

[5] 段秋亞,孟令奎,樊志偉,等.GF-1衛(wèi)星影像水體信息提取方法的適用性研究[J].國(guó)土資源遙感,2015,27(4):79-84.doi:10.6046/gtzyyg.2015.04.13. Duan Q Y,Meng L K,Fan Z W,et al.Applicability of the water information extraction method based on GF-1 image[J].Remote Sensing for Land and Resources,2015,27(4):79-84.doi:10.6046/gtzyyg.2015.04.13.

[6] 胡衛(wèi)國(guó),孟令奎,張東映,等.資源一號(hào)02C星圖像水體信息提取方法[J].國(guó)土資源遙感,2014,26(2):43-47.doi:10.6046/gtzyyg.2014.02.08. Hu W G,Meng L K,Zhang D Y,et al.Methods of water extraction from ZY-1 02C satellite imagery[J].Remote Sensing for Land and Resources,2014,26(2):43-47.doi:10.6046/gtzyyg.2014.02.08.

[7] 周林滔,楊國(guó)范,趙福強(qiáng),等.EMD與分形相結(jié)合的遙感影像水體信息提取方法[J].國(guó)土資源遙感,2014,26(4):41-45.doi:10.6046/gtzyyg.2014.04.07. Zhou L T,Yang G F,Zhao F Q,et al.Water information extraction from remote sensing image using EMD and fraction method[J].Remote Sensing for Land and Resources,2014,26(4):41-45.doi:10.6046/gtzyyg.2014.04.07.

[8] 袁迪,宋星原,張艷軍.基于遙感信息的新安江模型參數(shù)率定方法[J].武漢大學(xué)學(xué)報(bào):工學(xué)版,2014,47(1):23-27. Yuan D,Song X Y,Zhang Y J.Research on parameter calibration method of Xinanjiang model based on remote sensing information[J].Engineering Journal of Wuhan University,2014,47(1):23-27.

[9] 丁志雄,李紀(jì)人.流域洪水汛情的遙感監(jiān)測(cè)分析方法及其應(yīng)用[J].水利水電科技進(jìn)展,2004,24(3):8-11. Ding Z X,Li J R.Remote sensing monitoring analysis methods and their application to flood monitoring in river basins[J].Advances in Science and Technology of Water Resources,2004,24(3):8-11.

[10]Otsu N.A threshold selection method from gray-level histograms[J].IEEE Transactions on Systems,Man,and Cybernetics,1979,9(1):62-66.

[11]Osher S,Sethian J A.Fronts propagating with curvature-dependent speed:Algorithms based on Hamilton-Jacobi formulations[J].Journal of Computational Physics,1988,79(1):12-49.

[12]王伶俐,陳德清.2013年黑龍江大洪水遙感監(jiān)測(cè)分析[J].水文,2014,34(5):31-35,93. Wang L L,Chen D Q.Satellite remote sensing monitoring and analysis of Heilongjiang River flood in 2013[J].Journal of China Hydrology,2014,34(5):31-35,93.

(責(zé)任編輯： 陳理)

Approach to simulating the spatial-temporal process of flood inundation area

ZHANG Lianchong1,2, LI Guoqing1,3, YU Wenyang1,3, RAN Quan1,2

(1．KeyLabofDigitalEarthScience,InstituteofRemoteSensingandDigitalEarth,ChineseAcademyofSciences,Beijing100094,China; 2.UniversityofChineseAcademyofSciences,Beijing100049,China; 3.HainanLabofEarthObservation,Hainan572023,China)

Remote sensing data, as important information for flood disaster monitoring and loss assessment, can timely obtain the spatial-temporal distribution characteristics of flood. However, as it is restricted by weather conditions, it cannot form a dynamic and continuous process data. In this study, multi-temporal GF-1 satellite remote sensing clear images were used to extract the flood extent area based on bacha breach on the Heilong River in 2013. The flood inundation process was transformed into a numerical problem of partially differential equations by level set function. Finite difference method both in space and time was used to simulate the results of daily flood inundation area from August 24 to October 8. The results show that，compared with remote sensing data, the spatial-temporal consistency and the Kappa coefficients are 0.921 2 and 0.893 2; Compared with statistic data，the relatively error is less than 10%. This method has provided a scientific basis for the decision of flood disaster emergency response without prior information.

flood inundation area; spatial-temporal simulation; level set; remote sensing; Heilong River flood

10.6046/gtzyyg.2017.01.14

張連翀,李國(guó)慶,于文洋,等.基于水平集的洪澇淹沒(méi)范圍時(shí)空模擬方法[J].國(guó)土資源遙感,2017,29(1)：92-96.(Zhang L C,Li G Q,Yu W Y,et al.Approach to simulating the spatial-temporal process of flood inundation area[J].Remote Sensing for Land and Resources,2017,29(1):92-96.)

2015-09-17；

2015-11-19

國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目“多源遙感監(jiān)測(cè)數(shù)據(jù)在線融合及協(xié)同分析云平臺(tái)”(編號(hào)： 2016YFB0501504)和中國(guó)科學(xué)院數(shù)字地球重點(diǎn)實(shí)驗(yàn)室開(kāi)放基金項(xiàng)目“面向按需處理的遙感信息模型自動(dòng)化計(jì)算方法”(編號(hào)： 2015LDE005)共同資助。

張連翀(1985- ),男,博士研究生,主要從事高性能地學(xué)計(jì)算方面的研究。Email： zhanglc@radi.ac.cn。

TP 751.1

A

1001-070X(2017)01-0092-05