基于雷达影像的洪水监测技术方法详解-编程知识

洪水发生时候大多数是阴雨天气，光学影像基本上拍不到有效影像。雷达影像这时候就能发挥其不受天气影像的优点。现在星载的雷达卫星非常多，如高分三号、陆探一号、海丝一号（巢湖一号）、哨兵1号等。本文以哨兵1号L1地距(GRD)产品来介绍在洪水监测中的处理技术，其他雷达数据处理类似。

1 数据情况

因为要进行前后对比，因此需要两个时相哨兵1数据：

1、2023年8月2日：S1A_IW_GRDH_1SDV_20230802T214703_***.

2、2023年7月21日：S1A_IW_GRDH_1SDV_20230721T214702__***.

注：下载数据时候一般有三种格式下：SLC、GRD、GRD（COG），这里下载GRD即可。

2 处理流程

如下图为处理流程，可根据实际情况进行修改。

整个流程可分为两部分：SAR图像预处理和洪水信息提取。其中洪水信息提取采用的是分类后处理对比，即采用一种分类方法获取水灾前后的水面积信息，然后对两时相提取的水面积信息做对比得到水面积变化信息，增加部分就是洪水淹没部分。

3 SAR图像预处理

SARscape中提供流程化处理工具：/SARscape/Basic/Intensity Time Series Workflow。也提供分步式处理工具，本文使用分步处理工具来介绍。

在处理之前，为SARscape选择一套适合哨兵1数据的默认处理参数：Toolbox中单击/SARscape/Preferences，在打开的界面中，单击Load Preferences，选择Sentinel TOPSAR。

3.1 数据导入

SARscape5.6之后版本直接支持哨兵1的压缩包格式，因此导入数据的时候不需要解压缩。

（1）Toolbox中，打开SARscape/Import Data/SAR Spaceborne/Single Sensor/Sentinel-1。

（2）数据输入面板（Input Files）：选择两期压缩包文件。

（3）可选文件面板（Option file）：可以选择一个矢量文件对SAR数据进行裁剪，这里不选择。

（4）参数设置面板（Parameters）：按照默认。

（5）数据输出面板（Output Files）：选择一个输入目录。

（6）点击Exec执行处理。

得到两个时相、两个极化四个SAR图像。

3.2 图像配准

由于是两个时相SAR数据，使用自动SAR图像配准工具进行配准处理。

（1）在Toolbox中，选择/SARscape/Basic/Intensity Processing/Coregistration。

（2）在打开的Coregistration面板中：

数据输入（Input Files）面板：选择一个SAR数据作为基准（Input Reference File）。其余的数据作为待配准影像（Input File List）。
可选文件（Optional Files）面板：

DEM File：选择一个DEM数据可以提高配准精度，特别是高分辨率SAR数据和地形起伏较大区域，以及配准不同入射角图像的配准。由于该地区为平原，这里不选择DEM数据。

参数设置（Parameters）面板：按照默认参数。

注：设置Global->Generate Quick Look:False，不输出快视图，能提高一点处理速度。

数据输出（Output Files）面板：选择输出目录。

（3）单击Exec按钮执行配准过程。

执行完之后生成一个后缀为_meta的索引文件，每一景数据的配准结果作为一个波段在该文件中，类似于一个多波段数据，以及一个.series时间序列文件。

配准精度基本能达到亚像素，如下为配准效果图。

3.3 图像滤波

斑点类似噪音信号特征，一般产生于相干系统，比如SAR和激光雷达系统。表现为随机的“椒盐”散落在图像上。

SARscape提供两大类滤波，用于单波段雷达图像的滤波和多时相雷达图像滤波。本例子中使用多时相雷达图像滤波，效果比卷积滤波要好。

(1) 在Toolbox中，选择/SARscape/Basic/Intensity Processing/Filtering/De Grandi Spatio-Temporal Filtering。

(2) 打开De Grandi Spatio-Temporal Filtering面板

数据输入（Input Files）项：选择配准好的四个数据。
参数设置（Parameters）项：按照默认。

注：设置Global->Generate Quick Look:False，不输出快视图，能提高一点处理速度。

Output Files面板，数据输出路径和文件名按照默认。文件名自动添加_fil后缀。

(3) 单击Exec执行。

多时相滤波阈值斑点噪声的效果还是非常好的，如下图所示。

3.4 地理编码&辐射定标

这一步将前面配准、滤波处理得到的结果进行地理编码和辐射定标，得到各个时相数据的雷达强度图像或者后向散射系数图像。

该步骤可选择DEM数据参与地理编码提高几何精度，这里使用/SARscape/Import Data/DEM Extraction/ALOS World 3D 30m 工具自动下载SAR图像范围内的DEM，操作比较简单这里不详细介绍。

（1）在Toolbox中，选择/SARscape/Basic/Intensity Processing/Geocoding/Geocoding and Radiometric Calibration。

（2）打开地理编码和辐射定标面板（Geocoding and Radiometric Calibration）。

Input file：鼠标单击Brower按钮，按住shift键，将上一步滤波得到的4景数据全选，点击打开，数据列在Input File List列表下。
Optional Files面板，Geometry GCP File和Area File（应用于选择"True area"后向散射系数反演方法时候）这两个文件是可选项，这里不使用这两个文件。
DEM/Cartographic System面板，输入DEM文件或投影信息。若是输入DEM数据，最后输出结果默认以DEM投影参数为准。如果不输入DEM数据，则设置Output Projection。这里输入DEM文件
Parameters面板，主要参数（Principal Parameters）
- 像元大小（X Grid Size）：15
- 像元大小（Y Grid Size）：15
- 辐射定标（Radiometric Calibration）：Ture

注：这里辐射定标为后向散射系数，如果采用雷达水指数方式提取水体，这里建议定标为后向散射系数。其他水体提取方法看情况来定。

散射面积（Scattering Area）：Local Incidence Angle
辐射归一化（Radiometric Normalization）：False
辐射归一化方法（Normalization Method）：True
局部入射角（Local Incidence Angle）：False
输出类型（Output type）：dB

注:输出类型linear和Db的关系为： =10*log

注：其他参数按照默认，绝大多数情况下使用默认参数就能得到较好的结果。

Output Files面板，输出路径和文件名按照默认，自动添加了_geo后缀。

（3）单击Exec执行。

在生成结果中，默认会生成带有坐标、tif格式的快视图（* _ql.tif），如果后面采用目视解译的方法提取水体信息，可以直接使用该快视图图像。如下图为生成的快视图。

注：背景透明显示了。

4洪水信息提取与成果整理

在SAR图像上，有水体的地方一般显示黑色。可选择很多方法来提取水体信息。

在应急情况下，其实手动提取水面积是最快，精度也最高的方式。当然如果事先有利用大量样本训练好的深度学习模型，深度学习方法也是效率很高的一种方法。如下为几种方法的总结：

方法	说明
阈值分割法	在平原区域效果良好，几乎不需要人工编辑，但是在山区，阴影同时都被提取出来了。需要使用分类后处理的工具进行人机交互处理。
决策树分类方法	加入DEM，可有效避免山区阴影的混淆，对于坡度大的区域效果显著。后期仍然需要一定量的人工编辑，但工作量较阈值法小。
深度学习方法	效果良好，需要前期进行充分的影像水体特征分析，表现为选择若干有代表性的子区域，可有效区分山体阴影和水体，前期工作量较大，优势是训练好的水体提取模型可以用于同类型多幅影像，多个地区的水体提取，适用于批量处理。
面向对象方法	可直接获取水域边界，加入DEM作为辅助数据，对山体阴影的抑制也有作用。
人工+魔术棒	精度最高，需要人工参与。