CN113129318B - 一种利用sar影像计算蓄滞洪区蓄水量的方法 - Google Patents

Abstract

本发明公开了一种利用SAR影像计算蓄滞洪区蓄水量的方法，包括以下步骤：S1：对SAR影像进行预处理，利用阈值分割方法，提取水体，得到水体二值图像；S2：根据洪水淹没范围，分别对水体二值图像和数字高程模型DEM数据进行影像裁剪，对裁剪后的水体二值图像和DEM图像进行平面投影，并完成波段合成；S3：利用DEM图像，计算水体二值图像中每一个水体像元的水深，结合水体像元的面积，得到每一个水体像元对应立方体的蓄水量；S4：计算水体二值图像中所有水体像元对应立方体的蓄水量之和，即蓄滞洪区的蓄水总量；计算蓄滞洪区内的蓄水量，为蓄滞洪区的科学规划、合理运用提供决策支持。

Description

一种利用SAR影像计算蓄滞洪区蓄水量的方法

技术领域

本发明涉及领域，特别是一种利用SAR影像计算蓄滞洪区蓄水量的方法。

背景技术

我国是一个洪涝灾害频发的国家，为了提高防洪能力，在沿江河低洼地区和湖泊等设立蓄滞洪区。蓄滞洪区是分洪区、蓄洪区、滞洪区的统称，是江河防洪体系的重要组成部分，是我国防御流域性洪水的一项重要水利工程。它具有应急分洪蓄水能力，当发生大洪水时，可以分洪、滞蓄洪水，实现蓄洪削峰，保障重点地区的防洪安全。目前，我国设立了98处国家级蓄滞洪区、50处地方蓄滞洪区以及黄河滩区2处，用于分蓄流域超额洪水。蓄滞洪区开启后，及时确定蓄滞洪区蓄水量对科学规划、合理运用蓄滞洪区，具有重要的经济和社会意义。

遥感技术具有宏观、客观、高效、便捷等特点。洪涝发生时，往往多云、多雨，光学卫星影像难以发挥作用。而合成孔径雷达SAR，不受云、雾、雨等影响，可以全天时、全天候获取地面观测资料，是开展洪涝等灾害监测的重要技术手段。2016年，我国发射第一颗空间分辨率达到1米的C频段多极化合成孔径雷达SAR成像卫星，具备12种成像模式，成像幅宽大，与高空间分辨率优势相结合，既能实现大范围普查，也能详查特定区域，已经成为我国遥感减灾救灾的重要数据源之一。

蓄水量计算的传统方法都依赖地面观测。主要包括两种技术，一是利用蓄滞洪区进水处的计量水闸，对水闸分洪流量过程进行积分；二是根据水文站监测的水位，结合水位库容特性曲线计算蓄水量。

利用蓄滞洪区进水水闸，计算泄洪水量，即为蓄水量；

进入蓄滞洪区的水，部分下渗进入土壤，土壤含水量饱和后，才会形成地表积水；部分水通过蒸发进入大气。另外，洪水流经蓄滞洪区时，还会有沿途损失。因此，蓄水量并不完全等于泄洪水量。

根据水文站监测的水位，结合水位库容特性曲线计算蓄水量。

一是蓄滞洪区没有水文站时，需要人工测量水位；二是水闸泄洪时，水面并不是平的，而有限的测量点不能反映真实的水位情况，不能得到准确的蓄水量；三是随着人类活动的增加，蓄滞洪区内的土地利用已发生变化，因此，设计水位-容积曲线不能反映真实的水位与蓄水量关系，无法计算正确的蓄水量。

发明内容

针对现有依赖地面观测数据获取蓄滞洪区蓄水量的弊端，本发明提出了一种利用SAR影像计算蓄滞洪区蓄水量的方法，利用SAR影像宏观、高效，以及全天时、全天候的优势，不依赖地面观测数据，获取蓄滞洪区水面，结合数字高程模型DEM数据，计算蓄滞洪区内的蓄水量，为蓄滞洪区的科学规划、合理运用提供决策支持。

一种利用SAR影像计算蓄滞洪区蓄水量的方法，包括以下步骤：

S1：对SAR影像进行预处理，利用阈值分割方法，提取水体，得到水体二值图像；

S2：根据洪水淹没范围，分别对水体二值图像和数字高程模型DEM数据进行影像裁剪，对裁剪后的水体二值图像和DEM图像进行平面投影，并完成波段合成；

S3：利用DEM图像，计算水体二值图像中每一个水体像元的水深，结合水体像元的面积，得到每一个水体像元对应立方体的蓄水量；

S4：计算水体二值图像中所有水体像元对应立方体的蓄水量之和，即蓄滞洪区的蓄水总量。

优选地，S1包括以下子步骤：

S11：对SAR影像预处理，包括辐射校正、多视、滤波和地理编码；

S12：对经过预处理后的SAR影像，利用阈值分割方法，提取出水体；

S13：生成水体二值图像，其中，水体像元值为1，非水体像元值为0。

优选地，S2包括以下子步骤：

S21：根据蓄滞洪区内洪水淹没范围，分别对分别对水体二值图像和DEM数据进行影像裁剪，获取淹没范围的水体二值图像和DEM图像；

S22：对淹没范围的水体二值图像和DEM图像进行平面投影，投影后图像空间分辨率为R米；

S23：利用ENVI软件的Layer Stacking模块，对完成平面投影的水体二值图像和DEM图像进行波段合成。

优选地，S3包括以下子步骤：

S31：对水体二值图像Water中某一个水体像元(i，j)，求算水位值Z；根据DEM图像，该像元垂直对应的蓄滞洪区高程值为H，那么该像元的水深L＝Z-H；

S32：水体像元面积S＝R²，水体像元对应立方体的蓄水量V＝S·L＝R²·(Z-H)；

S33：水体像元总数为N，将所有水体像元对应立方体的蓄水量相加，得到蓄滞洪区总蓄水量

优选地，S31包括以下子步骤：

S311：设置滑动窗口，以(i，j)为起点，向左一步步滑动，步长为1，移动步数记为left，直到Water(i，j-left)＝1，且Water(i，j-left-1)＝0，即可判断找到水体左边界像元，左边界像元的高程值DSM(i，j-left)为Z_left；

S312：同S311，依次找到右边界像元(i，j+right)、上边界像元(i+up，j)和下边界像元(i-down，j)，得到Z_right、Z_up和Z_down；

S313：水体像元(i，j)的水位为Z_left、Z_right、Z_up和Z_down的最大值，即Z＝max{Z_left,Z_right,Z_up,Z_down}。

本发明利用SAR卫影像计算蓄滞洪区蓄水量的方法的有益成果：

不依赖地面观测数据，利用雷达全天时、全天候观测的优势，结合数字高程数据产品，计算蓄滞洪区蓄水量，为科学规划、合理运用蓄滞洪区提供数据数据支撑，避免洪涝造成的巨大经济损失，具有重要的经济和社会价值。

附图说明

图1为总流程图。

图2(a)为步骤3.1的平面图，图2(b)为步骤3.1的剖面图。

图3为步骤1.1预处理后的影像。

图4(a)为淹没区水体二值图像；图4(b)为淹没区DEM图像。

图5为蒙洼蓄洪区的总蓄水量与水闸测量的泄洪量对比图。

具体实施方式

下面结合附图和具体实施方式，对本发明作进一步说明。

下面对本发明的具体实施方式进行描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

本发明提出了一种利用SAR影像计算蓄滞洪区蓄水量的方法，总体技术方案包括：

步骤1：对SAR影像进行预处理，利用阈值分割方法，提取水体，得到水体二值图像；

步骤2：根据洪水淹没范围，分别对水体二值图像和数字高程模型DEM数据进行影像裁剪，对裁剪后的水体二值图像和DEM图像进行平面投影，并完成波段合成；

步骤3：利用DEM图像，计算水体二值图像中每一个水体像元的水深，结合水体像元的面积，得到每一个水体像元对应立方体的蓄水量；

步骤4：计算水体二值图像中所有水体像元对应立方体的蓄水量之和，即蓄滞洪区的蓄水总量。

步骤1包括：

(1.1)对SAR影像预处理，包括辐射校正、多视、滤波和地理编码等；

(1.2)对经过预处理后的SAR影像，利用阈值分割方法，提取出水体；

(1.3)生成水体二值图像，其中，水体像元值为1，非水体像元值为0。

步骤2包括：

(2.1)根据蓄滞洪区内洪水淹没范围，分别对分别对水体二值图像和DEM数据进行影像裁剪，获取淹没范围的水体二值图像和DEM图像；

(2.2)对淹没范围的水体二值图像和DEM图像进行平面投影，投影后图像空间分辨率为R米；

(2.3)利用ENVI软件的Layer Stacking模块，对完成平面投影的水体二值图像和DEM图像进行波段合成。

步骤3包括：

(3.1)对水体二值图像Water中某一个水体像元(i，j)，求算水位值Z；根据DEM图像，该像元垂直对应的蓄滞洪区高程值为H，那么该像元的水深L＝Z-H，示意图见图2，其中，图2(a)为平面图，图2(b)为剖面图；

包括：

(3.1.1)设置滑动窗口，以(i，j)为起点，向左一步步滑动，步长为1，移动步数记为left，直到Water(i，j-left)＝1，且Water(i，j-left-1)＝0，即可判断找到水体左边界像元，左边界像元的高程值DSM(i，j-left)为Z_left；

(3.1.2)同(3.1.1)步骤，依次找到右边界像元(i，j+right)、上边界像元(i+up，j)和下边界像元(i-down，j)，得到Z_right、Z_up和Z_down；

(3.1.3)水体像元(i，j)的水位为Z_left、Z_right、Z_up和Z_down的最大值，即Z＝max{Z_left,Z_right,Z_up,Z_down}；

(3.2)水体像元面积S＝R²，水体像元对应立方体的蓄水量V＝S·L＝R²·(Z-H)。

(3.3)水体像元总数为N，将所有水体像元对应立方体的蓄水量相加，得到蓄滞洪区总蓄水量

本发明在实施时，本发明实例选取淮河中部的蒙洼蓄洪区。蒙洼蓄洪区位于安徽省阜阳市阜南、颍上两县境内，是淮河流域于1953年设立的第一座行蓄洪区，也是我国重要的国家级蓄洪区之一。2020年7月，淮河流域发生区域性大洪水，7月20日上午8时31分，淮河王家坝闸开闸放水，蒙洼蓄洪区启用蓄洪。收集了7月20日13时和18时的SAR影像，利用本发明提出的方法，分别计算了泄洪时蒙洼蓄洪区的蓄水量，并与同一时间王家坝水闸实际泄洪量进行了对比。

步骤(1)：获取蒙洼蓄洪区SAR影像，完成SAR影像预处理，提取水体，生成水体二值图像；

步骤(11)：获取覆盖蒙洼蓄洪区SAR影像，完成辐射校正、多视、滤波和地理编码等预处理，预处理后的影像如图3所示；

步骤(12)：利用阈值分割方法，提取水体；

步骤(13)：生成水体二值图像，其中，水体像元值为1，非水体像元值为0。

步骤(2)：收集蒙洼蓄洪区的DEM数据，分别对水体二值图像和DEM数据进行影像裁剪，获取淹没区的水体二值图像和DEM图像，对裁剪后的水体二值图像和DEM图像进行平面投影，平面投影坐标为UTM Zone 50N，空间分辨率为25米，如图4所示，其中，图4(a)为淹没区水体二值图像，图4(b)为淹没区DEM图像；

步骤(3)：利用DEM图像，计算水体二值图像中每一个水体像元的水深Z，水体像元的面积S＝25m*25m＝625m²，得到每一个水体像元对应立方体的蓄水量V；

步骤(4)：计算淹没区所有水体像元对应立方体的蓄水量之和，即蒙洼蓄洪区的总蓄水量，与水闸测量的泄洪量对比如图5。利用本发明提出的方法估算的蓄水量，与水闸泄洪量相比，偏差在5％左右，偏差主要是水体下渗和蒸发导致。

Claims (1)

1.一种利用SAR影像计算蓄滞洪区蓄水量的方法，其特征在于，包括以下步骤：

S1：对SAR影像进行预处理，利用阈值分割方法，提取水体，得到水体二值图像；

S2：根据洪水淹没范围，分别对水体二值图像和数字高程模型DEM数据进行影像裁剪，对裁剪后的水体二值图像和DEM图像进行平面投影，并完成波段合成；

S3：利用DEM图像，计算水体二值图像中每一个水体像元的水深，结合水体像元的面积，得到每一个水体像元对应立方体的蓄水量；

S4：计算水体二值图像中所有水体像元对应立方体的蓄水量之和，即蓄滞洪区的蓄水总量；

所述S1包括以下子步骤：

S11：对SAR影像预处理，包括辐射校正、多视、滤波和地理编码；

S12：对经过预处理后的SAR影像，利用阈值分割方法，提取出水体；

S13：生成水体二值图像，

其中，水体像元值为1，非水体像元值为0；

所述S2包括以下子步骤：

S21：根据蓄滞洪区内洪水淹没范围，分别对水体二值图像和DEM数据进行影像裁剪，获取淹没范围的水体二值图像和DEM图像；

S22：对淹没范围的水体二值图像和DEM图像进行平面投影，投影后图像空间分辨率为R米；

S23：利用ENVI软件的Layer Stacking模块，对完成平面投影的水体二值图像和DEM图像进行波段合成；

所述S3包括以下子步骤：

S31：对水体二值图像Water中某一个水体像元(i，j)，求算水位值Z；根据DEM图像，该像元垂直对应的蓄滞洪区高程值为H，该像元的水深L＝Z-H；

S32：水体像元面积S＝R²，水体像元对应立方体的蓄水量V＝S·L＝R²·(Z-H)；

S33：水体像元总数为N，将所有水体像元对应立方体的蓄水量相加，得到蓄滞洪区总蓄水量

所述S31包括以下子步骤：

S311：设置滑动窗口，以(i，j)为起点，向左一步步滑动，步长为1，移动步数记为left，直到Water(i，j-left)＝1，且Water(i，j-left-1)＝0，即可判断找到水体左边界像元，左边界像元的高程值DEM(i，j-left)为Z_left；

S312：同S311，依次找到右边界像元(i，j+right)、上边界像元(i+up，j)和下边界像元(i-down，j)，得到Z_right、Z_up和Z_down；

S313：水体像元(i，j)的水位为Z_left、Z_right、Z_up和Z_down的最大值，即Z＝max{Z_left,Z_right,Z_up,Z_down}。

Images