本发明涉及水力学模拟的技术领域，特别是一种基于陆表水体产品的雷达数据水体信息提取方法及系统。
背景技术：
快速获取水体信息分布对于洪涝灾情评估、生态环境监测和水资源调查等方面具有重大意义。遥感技术是快速准确获取水体信息的重要手段。光学遥感由于自身波长较短，容易受到天气的影响，所以在应用上受到了一定的限制。而微波遥感中，星载合成孔径雷达（syntheticapertureradar,sar）具有全天时、全天候等特点，不受云、雨、雾的影响，在夜间也能成像，使得其成为洪涝灾害监测和湖泊动态监测的有力工具。欧空局于2014年4月发射的sentinel-1a系列卫星，单颗卫星的星轨周期为12天，两颗卫星为6天，宽幅最大为400km，分辨率最高为5m，具有单极化（hh/vv）和双极化（hh hv/vv vh）成像，而且数据免费，全球用户均可下载，获取方便。
目前sar水体信息提取主要有两种方法:一种是基于阈值分割的方法，一种是基于分类的方法。阈值法中应用较多的有单阈值法，多阈值法，双峰法和全局阈值法等。但在大范围的水体信息提取中，水体占比相对较低，阈值相对难于取值，而且由于受到雷达成像机理自身的影像，影像中会存在山地阴影，地形畸变，风浪导致水体表面有噪声，阈值法提取效果往往不好。而基于分类来识别大范围水体已经被征明是一种成功的算法。但基于分类的算法往往都是一种有监督的方法，需要足够的训练样本，这些样本传统上都是由人工进行选取，这极大的阻碍了水体信息提取的自动化效率。目前，已有利用现有的数据集，如modis和srtm衍生的水体掩膜来对模型进行训练，从而进行水体信息的自动提取,这种方法具有很大的潜力。但由于modis和srtm的空间分辨率较低，影响了水体提取的精度和效果。由欧盟联合研究中心jrc(jointresearchcentre)根据1984-2018的landsat数据制作的全球地表水体产品(分辨率为30m，有多种产品)，总体精度高达99.6%，可作为分类法中水体信息精细化提取的训练样本。
wenlihuang等人于2018年在remotesensing上提出了一种基于雷达卫星的水体信息自动提取方法，方法如下：(1)基于srtm水体数据集swbd和复合动态地表水范围cdswe生成水体标签，使用雷达后向散射系数、衍生指数和局部入射角作为特征，和水体标签组成训练样本；(2)通过水体样本和非水体样本的比例对样本数据进行随机选择部分训练样本；(3)使用随机采样的训练样本对随机森林进行训练，然后将训练的模型应用于所有像素，生成水体概率图，使用高概率水体、中概率水体、低概率水体和非水体标签对概率图进行分类。该方法提取速度较快，精度相对较高，但其未考虑训练样本中由于季节和时间等原因产生的显著错误样本，导致提取精度受到影响；同时其未考虑山地阴影的影响，提取水体中存在部分山地阴影。
技术实现要素：
为了解决上述的技术问题，本发明提出的一种基于陆表水体产品的雷达数据水体信息提取方法及系统，jrc的全球水体产品被用于自动生成训练数据集，并对数据集进行筛选，然后使用机器学习的分类方法对水体信息进行自动提取，并对提取的结果进行后处理，剔除小图斑和山地阴影的影像，从而得到精确的水体信息产品。
本发明的第一目的是提供一种基于陆表水体产品的雷达数据水体信息提取方法，包括获取卫星数据，还包括以下步骤：
步骤1：将所述卫星数据进行预处理,生成后向散射数据并计算所述后向散射数据的衍生系数；
步骤2：使用全球水体产品数据、所述后向散射数据和所述衍生系数创建训练样本数据集；
步骤3：基于随机森林模型对训练样本数据进行水体信息提取，得到初始水体提取结果；
步骤4：对所述初始水体提取结果进行后处理，得到最终的水体产品。
优选的是，所述预处理为对所述卫星数据经过基本处理后的1级地距产品的强度进行轨道矫正、辐射校正、斑点滤波、多视、地形矫正和db转换操作，最终生成后向散射数据vh和vv。
在上述任一方案中优选的是，述衍生系数包括极化比vhrvv，计算公式为
在上述任一方案中优选的是，所述衍生系数还包括归一化偏差极化指数ndpi，计算公式为
在上述任一方案中优选的是，所述衍生系数还包括归一化vh指数nvhi，计算公式为
在上述任一方案中优选的是，所述衍生系数还包括归一化vv指数nvvi，计算公式为
在上述任一方案中优选的是，所述步骤2包括将水体样本和非水体样本分别按所述vv值进行排序，选择所述水体样本中分位为n1%的值作为最大值阈值，删除大于最大值阈值的样本；非水体样本中分位为n2%的值作为最小值阈值，删除小于最小值阈值的样本。
在上述任一方案中优选的是，所述步骤3包括以下子步骤：
步骤31：随机选择训练数据集中的n3%的水体样本，当所述水体样本数量少于数量阈值时，选择全部水体样本；
步骤32：使用交叉验证网格gridsearch来确定随机森林分类器的最佳参数，生成训练模型；
步骤33：将所述训练模型应用于整个图像的所有像素，得到每个像素为水体的概率；
步骤34：根据阈值得到水体分布，得到初始水体提取结果。
在上述任一方案中优选的是，所述步骤4包括以下子步骤：
步骤41：利用srtm30的数字高程模型dem计算得到的坡度图；
步骤42：将所述初始水体提取结果中对应的坡度大于坡度阈值的提取结果赋值为非水体；
步骤43：通过种子点扩散的方法来删除小于面积阈值的水体；
步骤44：通过区域增长算法补充完整水体的边缘部分；
步骤45：生成最终的水体产品。
本发明的第二目的是提供一种基于陆表水体产品的雷达数据水体信息提取系统，包括用于获取卫星数据的数据获取模块，还包括以下模块：
预处理模块：用于将所述卫星数据进行预处理,生成后向散射数据并计算所述后向散射数据的衍生系数；
样本训练模块：用于使用全球水体产品数据、所述后向散射数据和所述衍生系数创建训练样本数据集；
信息提取模块：用于基于随机森林模型对训练样本数据进行水体信息提取，得到初始水体提取结果；
后处理模块：用于对所述初始水体提取结果进行后处理，得到最终的水体产品；
所述系统中的各个模块按照如第一目的所述的方法进行雷达数据水体信息提取。
优选的是，所述预处理为对所述卫星数据经过基本处理后的1级地距产品的强度进行轨道矫正、辐射校正、斑点滤波、多视、地形矫正和db转换操作，最终生成后向散射数据vh和vv。
在上述任一方案中优选的是，所述衍生系数包括极化比vhrvv，计算公式为
在上述任一方案中优选的是，所述衍生系数还包括归一化偏差极化指数ndpi，计算公式为
在上述任一方案中优选的是，所述衍生系数还包括归一化vh指数nvhi，计算公式为
在上述任一方案中优选的是，所述衍生系数还包括归一化vv指数nvvi，计算公式为
在上述任一方案中优选的是，所述样本训练模块用于将水体样本和非水体样本分别按所述vv值进行排序，选择所述水体样本中分位为n1%的值作为最大值阈值，删除大于最大值阈值的样本；非水体样本中分位为n2%的值作为最小值阈值，删除小于最小值阈值的样本。
在上述任一方案中优选的是，所述信息提取包括以下子步骤：
步骤31：随机选择训练数据集中的n3%的水体样本，当所述水体样本数量少于数量阈值时，选择全部水体样本；
步骤32：使用交叉验证网格gridsearch来确定随机森林分类器的最佳参数，生成训练模型；
步骤33：将所述训练模型应用于整个图像的所有像素，得到每个像素为水体的概率；
步骤34：根据阈值得到水体分布，得到初始水体提取结果。
在上述任一方案中优选的是，所述后处理包括以下子步骤：
步骤41：利用srtm30的数字高程模型dem计算得到的坡度图；
步骤42：将所述初始水体提取结果中对应的坡度大于坡度阈值的提取结果赋值为非水体；
步骤43：通过种子点扩散的方法来删除小于面积阈值的水体；
步骤44：通过区域增长算法补充完整水体的边缘部分；
步骤45：生成最终的水体产品。
本发明提出了一种基于陆表水体产品的雷达数据水体信息提取方法及系统，能够快速的自动提取大范围水体信息，而且不受山地阴影的影响，提取的水体较为完整和连通。
seasonality数据是指jrc全球水体产品的一种，它记录了最近一年内全球陆地区域每一个像元出现水体的域分次数，最高值为12，即该像元在一年内水体一直存在；最低为0，即非水体像元；1-11为季节性水体。最新的水体产品为基于2018年的landsat数据制作的seasonality2018。
附图说明
图1为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的一优选实施例的流程图。
图2为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的一优选实施例的流程图。
图3为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的一优选实施例的流程图。
图4为按照本发明的基于陆表水体产品的雷达数据水体信息提取系统的一优选实施例的模块图。
图5为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的另一优选实施例的技术路线图。
图6为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图3所示实施例的云阴影位移示意图。
图7为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的宽河道水体提取结果的一优选实施例的示意图。图8为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图7所示实施例的窄溪流水体提取结果示意图。
图9为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的洪水事件的一优选实施例的灾前卫星数据图。
图10为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图9所示实施例的灾前水体提取结果图。
图11为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图9所示实施例的灾中卫星数据图。
图12为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图9所示实施例的灾中水体提取结果图。
图13为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图9所示实施例的灾后卫星数据图。
图14为按照本发明的基于陆表水体产品的雷达数据水体信息提取方法的如图9所示实施例的灾后水体提取结果图。
具体实施方式
下面结合附图和具体的实施例对本发明做进一步的阐述。
实施例一
如图1所示，一种基于陆表水体产品的雷达数据水体信息提取方法，执行步骤100，获取卫星数据。卫星数据指的是sentinel-1a系列卫星返回的数据。
执行步骤110，将所述卫星数据进行预处理,生成后向散射数据并计算所述后向散射数据的衍生系数。预处理为对所述卫星数据经过基本处理后的1级地距产品的强度进行轨道矫正、辐射校正、斑点滤波、多视、地形矫正和db转换操作，最终生成后向散射数据vh和vv。衍生系数包括极化比vhrvv、归一化偏差极化指数ndpi、归一化vh指数nvhi和归一化vv指数nvvi。的极化比vhrvv计算公式为
执行步骤120，使用全球水体产品数据、所述后向散射数据和所述衍生系数创建训练样本数据集。将水体样本和非水体样本分别按所述vv值进行排序，选择所述水体样本中分位为n1%的值作为最大值阈值，删除大于最大值阈值的样本；非水体样本中分位为n2%的值作为最小值阈值，删除小于最小值阈值的样本。在本实施例中，n1=85，n2=15。
执行步骤130，基于随机森林模型对训练样本数据进行水体信息提取，得到初始水体提取结果。如图2所示，执行步骤131，随机选择训练数据集中的n3%的水体样本，当所述水体样本数量少于数量阈值时，选择全部水体样本，在本实施例中，n3=12.5。执行步骤132，使用交叉验证网格gridsearch来确定随机森林分类器的最佳参数，生成训练模型。执行步骤133，将所述训练模型应用于整个图像的所有像素，得到每个像素为水体的概率。执行步骤134，根据阈值得到水体分布，得到初始水体提取结果。
执行步骤140，对所述初始水体提取结果进行后处理，得到最终的水体产品。如图3所示，执行步骤141，利用srtm30的数字高程模型dem计算得到的坡度图。执行步骤142，将所述初始水体提取结果中对应的坡度大于坡度阈值的提取结果赋值为非水体。执行步骤143，通过种子点扩散的方法来删除小于面积阈值的水体。执行步骤144，通过区域增长算法补充完整水体的边缘部分。执行步骤145，生成最终的水体产品。
实施例二
如图4所示，一种基于陆表水体产品的雷达数据水体信息提取系统，包括数据获取模块200、预处理模块210、样本训练模块220、信息提取模块230和后处理模块240。
数据获取模块200：用于获取卫星数据。
预处理模块210：用于将所述卫星数据进行预处理,生成后向散射数据并计算所述后向散射数据的衍生系数。预处理为对卫星数据经过基本处理后的1级地距产品的强度进行轨道矫正、辐射校正、斑点滤波、多视、地形矫正和转换操作。衍生系数包括极化比vhrvv、归一化偏差极化指数ndpi、归一化vh指数nvhi和归一化vv指数nvvi。的极化比vhrvv计算公式为
样本训练模块220：用于使用全球水体产品数据、所述后向散射数据和所述衍生系数创建训练样本数据集。将水体样本和非水体样本分别按所述vv值进行排序，选择所述水体样本中分位为n1%的值作为最大值阈值，删除大于最大值阈值的样本；非水体样本中分位为n2%的值作为最小值阈值，删除小于最小值阈值的样本。
信息提取模块230：用于基于随机森林模型对训练样本数据进行水体信息提取，得到初始水体提取结果。信息提取包括以下子步骤：步骤31：随机选择训练数据集中的n3%的水体样本，当所述水体样本数量少于数量阈值时，选择全部水体样本；步骤32：使用交叉验证网格gridsearch来确定随机森林分类器的最佳参数，生成训练模型；步骤33：将所述训练模型应用于整个图像的所有像素，得到每个像素为水体的概率；步骤34：根据阈值得到水体分布，得到初始水体提取结果。
后处理模块240：用于对所述初始水体提取结果进行后处理，得到最终的水体产品。后处理包括以下子步骤：步骤41：利用srtm30的数字高程模型dem计算得到的坡度图；步骤42：将所述初始水体提取结果中对应的坡度大于坡度阈值的提取结果赋值为非水体；步骤43：通过种子点扩散的方法来删除小于面积阈值的水体；步骤44：通过区域增长算法补充完整水体的边缘部分；步骤45：生成最终的水体产品。
实施例三
本发明提出的水体信息自动提取主要包括四个步骤：
（1）将sentinel-1sar数据进行预处理,生成后向散射数据并计算其衍生系数；
（2）使用seasonality2018产品和后向散射数据及衍生系数创建训练样本数据集；
（3）基于随机森林模型进行水体信息提取；
（4）后处理：首先去除显著错分类水体;然后使用坡度数据去除山体阴影;再剔除小连通区域，使用区域增长补充水体，得到最终的水体产品；最后使用辅助图像对产品进行精度评估。
技术路线如图5所示。
1、sar数据预处理与指数计算
对高分辨率的sentinel-1的level-1地面检测产品grdh（high-resolutiongroundrangedetected,5mx20m）的强度进行轨道矫正，辐射矫正，斑点滤波，多视，地形矫正，转换为db。其中斑点滤波采用refine-lee滤波，多视中将图像转换为20m的正方形像素效果较好(10m的斑点噪声较多，30m则分辨率较低);使用srtm1弧秒（约30m）的数字高程模型进行地形矫正。
同时进行相关指数计算,并将其作为特征数据加入到分类中。其中包括极化比vhrvv，归一化偏差极化指数ndpi，归一化vh指数nvhi，归一化vv指数nvvi。指数计算如表1所示。
表1极化指数列表
2、水体训练样本创建
样本数据来源于jrcglobalsurfacewater的seasonality数据。这一数据记录的是最近一年内全球区域内陆地每一个水体像元出现的月份次数，最高值为12，即一年内水体一直存在；最低为0，即非水像元；1-11为季节性水体。最新的产品为seasonality2018。由于季节和精度的问题，水体样本和非水体样本中会出现较少的错分样本，所以本文结合雷达影像中水体的特性，首先将水体样本和非水体样本进行按vv值进行排序，然后选择水体样本中位于85%的分位的值作为最大值阈值，删除大于最大值阈值的样本；非水体样本中分位为15%的值作为最小值阈值，删除其中小于最小值阈值的样本。
3、区域水体信息提取
随机森林模型被用于从sentinel-1数据中分类出水体。随机森林在当前的所有算法中，具有极好的准确率。机器学习的分类方法会受到训练数据不平衡的情况影响，对于样本较少的类别，往往预测效果不太好。所以对于水体与非水体的比例大于1：20时，本文限制为1：20。为了简化训练样本，本文随机采样训练数据集中的1/8的水体样本，当水体样本少于10000时，选择全部水体样本，非水体按比例随机选取。同时使用交叉验证网格来确定随机森林分类器的最佳参数，例如使用的树的棵树和最大特征数。然后将训练出的模型应用于整个图像的所有像素，得到每个像素为水体的概率。最后根据阈值得到水体分布，本文选择大于0.5即为水体，从而得到初始水体提取结果。
4、后处理
经过信息提取后，有少部分异常的高值被错误分类为水体，为去这部分值的影响，将vv值大于水体最大阈值的水体赋值为非水体。在经过随机森林后的数据进行分类后，大部分水体已经被提取出来，但由于雷达图像的自身特征，山体阴影也呈现出和水体一样的特征，本文利用srtm30的数字高程模型dem计算得到的坡度图，对于所有的水体，当其对应的坡度大于5°时，赋值为非水体，即大于一定坡度时，我们认定其不能留住水体。此处生成的水体产品中有较多零散的像素点，我们通过种子点扩散的方法来删除小于一定面积的水体。由于dem的精度和误差问题，完整水体的边缘部分被去除，所以通过区域增长算法补充这部分水体，从而得到最终的水体产品。
实施例四
本实施例选择了2019年8月10日浙江省台州市临海洪水事件作为研究对象，浙江临海市包含3个区县，区域内平均海拔300m，地形较为复杂，既有山区又有平原地区;区域范围内水体类型多样，集水面积约254平方公里。
使用的数据有sentinel-1sar，landsat8，jrcglobalsurfacewater，srtmdem数据。
sentinel-1号卫星是欧空局发射的第一个哥白尼计划卫星星座，其由两颗卫星a、b组合，两个卫星的组合下，可以达到每6天对同一地点进行影像获取。其为c波段卫星，有四种成像模式：超精细模式（stripmapmode，sm）、干涉宽幅模式(interferometricwideswath，iw)、超宽幅模式(extra-wideswathmode，ew)、微波模式(wave-mode，wave)。其中对地球陆地进行覆盖的主要为iw模式，sm模式主要用于应急事件，ew和wave主要用于海洋监测。sentinel-1号卫星数据可以免费获取，iw模式的宽幅为250km，分辨率为5x20m，从卫星拍摄到数据分发至数据库只要3-6小时左右，对地观测能够穿透云雾，不受天气影像，这些特性使其非常适合应用于洪涝灾害遥感监测。本文sentinel-1数据从欧空局网站下载(https://scihub.copernicus.eu/)，也可从nasa转发的网站下载(https://search.asf.alaska.edu)。
jrcglobalsufacewater水体产品主要用于制作训练样本用于校准随机森林模型。该水体产品由欧洲联合研究中心制作，该产品包含1984年至2018年地表水的位置和时间分布图，并提供了水面的范围和统计数据，这些数据是使用1984年3月16日至2018年12月31日从landsat5、7和8采集的3,865,618个场景生成的。使用专家系统将每个像素分别分类为水/非水，并将结果整理为每月的历史记录整个时间段和两个时期（1984-1999、2000-2018）进行更改检测。本文主要使用的数据为seasonality2018，它主要记录了在2018年里水体出现了几个月，12个月即为永久性水体，少于12个月则为季节性水体。
实施例五
如图6所示，展示了实例四的水体提取结果，大部分水体区域已被提取，漏提现象较少，基本不受山地阴影的影响。在图7中局部区域1中宽河道水体被较为完整的提取出，且未受右上角的山地阴影影响；图8中局部区域2中的较窄溪流和水塘也被提取出来，河道较为完整。
实施例六
如图9-14所示，分别展示了本发明在浙江临海洪水事件中灾前、灾中和灾后的卫星图像和水体提取结果，结果展示水体提取较为完整，能后较好的应用于洪涝灾害遥感监测。
为了更好地理解本发明，以上结合本发明的具体实施例做了详细描述，但并非是对本发明的限制。凡是依据本发明的技术实质对以上实施例所做的任何简单修改，均仍属于本发明技术方案的范围。本说明书中每个实施例重点说明的都是与其它实施例的不同之处，各个实施例之间相同或相似的部分相互参见即可。对于系统实施例而言，由于其与方法实施例基本对应，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
1.一种基于陆表水体产品的雷达数据水体信息提取方法，包括获取卫星数据，其特征在于，还包括以下步骤：
步骤1：将所述卫星数据进行预处理,生成后向散射数据并计算所述后向散射数据的衍生系数；
步骤2：使用全球水体产品数据、所述后向散射数据和所述衍生系数创建训练样本数据集；
步骤3：基于随机森林模型对训练样本数据进行水体信息提取，得到初始水体提取结果；
步骤4：对所述初始水体提取结果进行后处理，得到最终的水体产品。
2.如权利要求1所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述预处理为对所述卫星数据经过基本处理后的1级地距产品的强度进行轨道矫正、辐射校正、斑点滤波、多视、地形矫正和db转换操作，最终生成后向散射数据vh和vv。
3.如权利要求2所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述衍生系数包括极化比vhrvv，计算公式为
4.如权利要求2所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述衍生系数还包括归一化偏差极化指数ndpi，计算公式为
5.如权利要求2所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述衍生系数还包括归一化vh指数nvhi，计算公式为
6.如权利要求2所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述衍生系数还包括归一化vv指数nvvi，计算公式为
7.如权利要求2所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述步骤2包括将水体样本和非水体样本分别按所述vv值进行排序，选择所述水体样本中分位为n1%的值作为最大值阈值，删除大于最大值阈值的样本；非水体样本中分位为n2%的值作为最小值阈值，删除小于最小值阈值的样本。
8.如权利要求7所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述步骤3包括以下子步骤：
步骤31：随机选择训练数据集中的n3%的水体样本，当所述水体样本数量少于数量阈值时，选择全部水体样本；
步骤32：使用交叉验证网格gridsearch来确定随机森林分类器的最佳参数，生成训练模型；
步骤33：将所述训练模型应用于整个图像的所有像素，得到每个像素为水体的概率；
步骤34：根据阈值得到水体分布，得到初始水体提取结果。
9.如权利要求8所述的基于陆表水体产品的雷达数据水体信息提取方法，其特征在于，所述步骤4包括以下子步骤：
步骤41：利用srtm30的数字高程模型dem计算得到的坡度图；
步骤42：将所述初始水体提取结果中对应的坡度大于坡度阈值的提取结果赋值为非水体；
步骤43：通过种子点扩散的方法来删除小于面积阈值的水体；
步骤44：通过区域增长算法补充完整水体的边缘部分；
步骤45：生成最终的水体产品。
10.一种基于陆表水体产品的雷达数据水体信息提取系统，包括用于获取卫星数据的数据获取模块，其特征在于，还包括以下模块：
预处理模块：用于将所述卫星数据进行预处理,生成后向散射数据并计算所述后向散射数据的衍生系数；
样本训练模块：用于使用全球水体产品数据、所述后向散射数据和所述衍生系数创建训练样本数据集；
信息提取模块：用于基于随机森林模型对训练样本数据进行水体信息提取，得到初始水体提取结果；
后处理模块：用于对所述初始水体提取结果进行后处理，得到最终的水体产品；
所述系统中的各个模块按照如权利要求1所述的方法进行雷达数据水体信息提取。
技术总结