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

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

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

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

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

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

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

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

0 引言

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

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

1 研究區概況與數據源

1.1 研究區概況

2013年夏季汛期，黑龍江流域遭遇歷史罕見特大洪水。本文以淹沒范圍最大的黑龍江省同江市八岔段潰口淹沒區為研究區。該地區位于黑龍江干流下游南岸，地勢低平，黑龍江、松花江和烏蘇里江等多條河流匯合而來的洪水流經該區域速度變緩，洪峰持續時間長。2013年8月23日該地區出現潰口，造成當地3 000余戶房屋受損，農田成災面積達7.6萬hm2。

1.2 數據預處理

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

2 研究方法

水平集方法首先由Osher和Sethian提出，是處理運動曲面隨時間拓撲變化的有效計算工具[11]。該方法將連續函數φ(x,y,t)∶R3→R描述為閉合演化曲線C(p,t)： 0≤p≤1在t時刻的隱式表達，即t時刻曲線C(p,t)對應于φ(x,y,t)的零水平集。本文設t1時刻的洪澇淹沒范圍為源曲面φ1，t2時刻的洪澇淹沒范圍為靶曲面φ2(t1

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

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

(1)

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

采用初始水平集函數描述φ1后，根據演化速度場求解Hamilton-Jacobi方程，實現水平集函數的動態演化，即

(2)

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

(3)

(4)

(5)

(6)

(7)

(8)

(9)

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

3 實驗結果與分析

3.1 漲水過程模擬

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

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

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

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

3.2 退水過程模擬

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

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

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

3.3 模擬結果精度評價

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

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

表1 模擬結果與提取結果的精度統計

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

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

表2 模擬結果與統計結果的相對誤差分析

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

4 結論

1)基于水平集的洪澇淹沒范圍時空模擬方法，將片段化的洪澇淹沒范圍數據擴展為連續動態的淹沒過程信息，能夠實現洪澇淹沒范圍的快速模擬。研究表明，模擬結果與影像提取結果具有較好的時空一致性，為洪澇淹沒歷時提取、淹沒程度和災害損失評估等后續工作提供了科學依據。相比于依賴先驗資料的洪水演進模型，該方法具有明顯的通用性優勢，提高了高空間分辨率遙感影像在洪澇災害監測中的應用前景。

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

[1] 俞云利,賴錫軍.二維平面非恒定流數學模型的遙感水位數據同化[J].水科學進展,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] 沈定濤,王結臣,張煜,等.一種面向海量數字高程模型數據的洪水淹沒區快速生成算法[J].測繪學報,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的洪水淹沒模擬分析及虛擬現實表達[J].西南大學學報:自然科學版,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] 孟令奎,郭善昕,李爽.遙感影像水體提取與洪水監測應用綜述[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衛星影像水體信息提取方法的適用性研究[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] 胡衛國,孟令奎,張東映,等.資源一號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] 周林滔,楊國范,趙福強,等.EMD與分形相結合的遙感影像水體信息提取方法[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] 袁迪,宋星原,張艷軍.基于遙感信息的新安江模型參數率定方法[J].武漢大學學報:工學版,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].水利水電科技進展,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年黑龍江大洪水遙感監測分析[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.

(責任編輯： 陳理)

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

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

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

TP 751.1

A

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