張連翀， 李國慶， 于文洋， 冉全

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

基于水平集的洪澇淹沒范圍時空模擬方法

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

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

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

洪澇淹沒范圍； 時空模擬； 水平集； 遙感； 黑龍江洪水

0 引言

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

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

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

1.1 研究區(qū)概況

2013年夏季汛期，黑龍江流域遭遇歷史罕見特大洪水。本文以淹沒范圍最大的黑龍江省同江市八岔段潰口淹沒區(qū)為研究區(qū)。該地區(qū)位于黑龍江干流下游南岸，地勢低平，黑龍江、松花江和烏蘇里江等多條河流匯合而來的洪水流經(jīng)該區(qū)域速度變緩，洪峰持續(xù)時間長。2013年8月23日該地區(qū)出現(xiàn)潰口，造成當(dāng)?shù)? 000余戶房屋受損，農(nóng)田成災(zāi)面積達(dá)7.6萬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)過相對輻射校正、幾何糾正等預(yù)處理后，首先計算歸一化差異水體指數(shù)(normalized difference water index,NDWI)以最大限度地突出水體與陸地、植被之間的差異； 然后采用最大類間方差法[10]確定分割閾值提取水體邊界； 再對研究區(qū)背景信息和洪水淹沒信息進(jìn)行二值化處理，最終得到不同時刻的洪澇淹沒范圍。

2 研究方法

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

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

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

(1)

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

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

(2)

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

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

3 實驗結(jié)果與分析

3.1 漲水過程模擬

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

圖1-1 洪澇淹沒范圍時空模擬結(jié)果(漲水過程)

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

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

3.2 退水過程模擬

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

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

圖2-2 洪澇淹沒范圍時空模擬結(jié)果(退水過程)

3.3 模擬結(jié)果精度評價

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

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

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

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

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

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

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

4 結(jié)論

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

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

[1] 俞云利,賴錫軍.二維平面非恒定流數(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ù)的洪水淹沒區(qū)快速生成算法[J].測繪學(xué)報,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的洪水淹沒模擬分析及虛擬現(xiàn)實表達(dá)[J].西南大學(xué)學(xué)報:自然科學(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)測應(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].國土資源遙感,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)國,孟令奎,張東映,等.資源一號02C星圖像水體信息提取方法[J].國土資源遙感,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] 周林滔,楊國范,趙福強(qiáng),等.EMD與分形相結(jié)合的遙感影像水體信息提取方法[J].國土資源遙感,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é)報:工學(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)測分析方法及其應(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)測分析[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

張連翀,李國慶,于文洋,等.基于水平集的洪澇淹沒范圍時空模擬方法[J].國土資源遙感,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

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

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

TP 751.1

A

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