张连翀， 李国庆， 于文洋， 冉全

(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