一种基于极化雷达自身影像特征的水体提取方法及系统与流程

1.本发明涉及洪涝灾害监测领域，更具体地，涉及一种基于极化雷达自身影像特征的水体提取方法及系统。


背景技术：

2.星载合成孔径雷达系统(synthetic aperture radar,sar)成像时不依靠太阳光，使用波长较长的电磁波探测地物，受云、雨、雾等天气的影响较小，具有全天时、全天候获取地表信息的优势，能够弥补光学遥感的局限性，非常适用于对洪涝灾害的及时监测，是研究洪涝灾害的有效工具。
3.雷达遥感系统常用的极化电磁波有两种：水平极化(h)和垂直极化(v)，传统单极化雷达只能发射和接收一种极化波，如hh(水平发射、水平接收)或vv(垂直发射、垂直接收)。在有限的波谱信息下很容易将水体与裸地、广场，以及停车场等地物混淆。而极化雷达可以同时发射和接收多种极化波，即hh，hv(水平发射、垂直接收)，vh(垂直发射、水平接收)和vv，可以提供多个图像通道。因此，与传统单极化雷达相比，极化雷达可以获取更多的目标散射信息，实现更高精度的地物信息提取。
4.在水体提取方面，基于极化雷达影像的算法一般分为阈值分割、面向对象、机器学习三类。其中，面向对象考虑了地物的纹理、颜色特征，水体提取结果的连续性较强，不易受影像斑点噪声的影响；机器学习则因较强的特征提取能力，可以获得精度较高的提取结果。但是，成功应用上述方法的前提是高质量的训练数据和扎实的模型参数调整能力。因此，这两种方法的应用多见于学术论文中的小区域示范。
5.而阈值分割法因计算速度快、算法简易的优点，在诸如洪水淹没范围提取的应急响应中的具有较广泛的应用。例如，贾诗超等(2019)通过两个极化波段构建了sdwi指数来提取水体信息，李景刚等(2010)通过改进的otsu算法实现水体提取。但是，该方法易受山体阴影等低回波地物的干扰，提取精度不高。目前普遍的改进方法是借助dem(数字高程模型，digital elevation model)数据，通过阈值分割生成山体掩膜，滤去初始的水体指数提取结果。这种做法一方面需要额外的数据输入，计算量较大；另一方面引入了第二套数据的阈值分割结果，存在新的误差。


技术实现要素：

6.本发明的目的在于克服基于极化雷达影像提取水体时，阈值分割法易受山体阴影干扰的缺点，提出一种基于极化雷达自身影像特征的水体提取方法。该方法充分利用了极化雷达影像自身回波信号的统计特征，在不需要外部dem输入的情况下，即可消山体阴影的影响，实现影像图幅范围内的水体高精度提取。
7.为解决上述技术问题，本发明的技术方案如下：
8.本发明提供一种基于极化雷达自身影像特征的水体提取方法，包括以下步骤：
9.s1：根据需要开展水体提取的时间、地点，选择相应的极化雷达影像并完成预处
理，得到不同极化条件下的后向散射系数σ；
10.s2：绘制不同极化条件下的后向散射系数σ的统计直方图，根据统计直方图的双峰特征，依据预设的自动阈值提取算法，定位到双峰之间的谷底位置，提取局部最小值作为原极化雷达影像的水体阈值；
11.s3：计算不同极化雷达影像水体阈值之差的均值，取该均值的绝对值作为阈值偏移量t
δ
；
12.s4：针对双峰特征较弱的极化方式，在其后向散射系数σ上叠加阈值偏移量t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠；
13.s5：比较每个像元在不同极化方式中的后向散射系数大小，取最小值构成一景新影像σ
min
；
14.s6：构建水体提取指数，并根据步骤s2的方法计算水体阈值，提取图幅范围内的水体像元。
15.优选地，步骤s1具体包括：
16.s1.1：根据需要开展水体提取的时间、地点，下载相应的vv和vh极化雷达影像；
17.s1.2：对vv和vh极化雷达影像进行应用轨道文件校正处理；
18.s1.3：对vv和vh极化雷达影像进行热噪声去除处理；
19.s1.4：对vv和vh极化雷达影像进行斑点滤波处理；
20.s1.5：对vv和vh极化雷达影像进行地理编码处理；
21.s1.6：对vv和vh极化雷达影像进行分贝化处理，得到vv和vh两种极化方式的分贝化雷达后向散射系数，分别记为σ
vv
和σ
vh
。
22.优选地，步骤s2具体包括：
23.s2.1：对vv和vh极化雷达影像进行中值滤波，得到极化雷达影像后向散射系数σ
vv
和σ
vh
的统计直方图；
24.s2.2：按后向散射系数由小到大的顺序，取每n个区间(bin)的频数最大值，设该最大值为因变量y，对应bin所在的后向散射系数为自变量x，构建拟合函数f(x,y)，以达到平滑原始直方图数据中的频数曲线，降低局部bin之间的锯齿扰动的效果；
25.s2.3：根据拟合函数f(x,y)的一阶导数判断收敛性，找到频数的局部最小值y0，对应的bin即为直方图双峰之间的谷底，进而得bin所在的后向散射系数作为提取原始vv和vh极化雷达影像的水体阈值。
26.优选地，步骤s2.1中，中值滤波采用的滤波窗口大小为5
×
5像元。
27.优选地，步骤s2.2中，取n＝5。
28.优选地，步骤s3中，计算vv和vh极化雷达影像的阈值之差，取其绝对值作为阈值偏移量t
δ
。
29.优选地，步骤s4中，针对双峰特征较弱的vv极化雷达影像，每个像元均减去阈值偏移量t
δ
，即σ’vv
＝σ
vv-t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠。
30.优选地，步骤s5中，比较每个像元的σ’vv
和σ
vh
大小，取每个像元σ’vv
和σvh的最小值构成一景新影像，即全域最小值影像σ
min
。
31.优选地，步骤s6中，具体包括：
32.s6.1：构建水体提取指数pwi＝lg(σ
min
×
σ
min
)；
33.s6.2：按步骤s2的方法提取水体阈值t
wi
；
34.s6.1：按pwi》t
wi
提取水体范围。
35.本发明还提供一种基于极化雷达自身影像特征的水体提取系统，包括存储器和处理器，所述存储器中包括基于极化雷达自身影像特征的水体提取方法程序，所述基于极化雷达自身影像特征的水体提取方法程序被所述处理器执行时实现如下步骤：
36.s1：根据需要开展水体提取的时间、地点，选择相应的极化雷达影像并完成预处理，得到不同极化条件下的后向散射系数σ；
37.s2：绘制不同极化条件下的后向散射系数σ的统计直方图，根据统计直方图的双峰特征，依据预设的自动阈值提取算法，定位到双峰之间的谷底位置，提取局部最小值作为原极化雷达影像的水体阈值；
38.s3：计算不同极化雷达影像水体阈值之差的均值，取该均值的绝对值作为阈值偏移量t
δ
；
39.s4：针对双峰特征较弱的极化方式，在其后向散射系数σ上叠加阈值偏移量t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠；
40.s5：比较每个像元在不同极化方式中的后向散射系数大小，取最小值构成一景新影像σ
min
；
41.s6：构建水体提取指数，并根据步骤s2的方法计算水体阈值，提取图幅范围内的水体像元。
42.与现有技术相比，本发明技术方案的有益效果是：本发明提供一种基于极化雷达自身影像特征的水体提取方法及系统，根据需要开展水体提取的时间、地点，选择相应的极化雷达影像并完成预处理，得到不同极化条件下的后向散射系数σ；绘制不同极化条件下的后向散射系数σ的统计直方图，根据统计直方图的双峰特征，依据预设的自动阈值提取算法，定位到双峰之间的谷底位置，提取局部最小值作为原极化雷达影像的水体阈值；计算不同极化雷达影像水体阈值之差的均值，取该均值的绝对值作为阈值偏移量t
δ
；针对双峰特征较弱的极化方式，在其后向散射系数σ上叠加阈值偏移量t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠；比较每个像元在不同极化方式中的后向散射系数大小，取最小值构成一景新影像σ
min
；最后构建水体提取指数，并计算水体阈值，提取图幅范围内的水体像元。本发明克服基于极化雷达影像提取水体时，阈值分割法易受山体阴影干扰的缺点，充分利用了极化雷达影像自身回波信号的统计特征，在不需要外部dem输入的情况下，即可消山体阴影的影响，实现影像图幅范围内的水体高精度提取。
附图说明
43.图1为实施例1基于极化雷达自身影像特征的水体提取方法的流程图。
44.图2为哨兵1号grd数据vv极化雷达影像图。
45.图3为哨兵1号grd数据vh极化雷达影像图。
46.图4为σ
vv
的统计直方图。
47.图5为σ
vh
的统计直方图。
48.图6为修改后的vv和原始vv、vh的频数分布对比图。
49.图7为全域最小值影像σ
min
的示意图。
50.图8为水体提取指数pwi的直方图。
51.图9为通过本发明pwi方法和现有sdwi水体提取的效果对比图。
52.图10为图9中a、b、c、d区域的水体提取的效果对比图。
具体实施方式
53.为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式对本发明进行进一步的详细描述。需要说明的是，在不冲突的情况下，本技术的实施例及实施例中的特征可以相互组合。
54.在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用其他不同于在此描述的其他方式来实施，因此，本发明的保护范围并不受下面公开的具体实施例的限制。
55.实施例1
56.如图1所示，本实施例提供基于极化雷达自身影像特征的水体提取方法，包括以下步骤：
57.s1：根据需要开展水体提取的时间、地点，选择相应的极化雷达影像并完成预处理，得到不同极化条件下的后向散射系数σ(单位：db)；
58.在具体实施过程中，步骤s1包括：
59.s1.1：如图2-3所示，根据需要开展水体提取的时间、地点，下载相应的vv(垂直发射、垂直接收)和vh(垂直发射、水平接收)极化雷达影像；本实施例以哨兵1号卫星地距多视影像(ground range detected，grd)数据为例，水体提取的地点为为珠江河口地区，水体提取的时间(成像日期)为2021年12月4日，极化方式为vv和vh。影像通过哨兵系列卫星科研数据中心网站(sentinels scientific data hub，https://scihub.copernicus.eu/dhus/#/home)下载获取，原始影像文件名为s1a_iw_grdh_1sdv_20211204t103405_20211204t103430_040858_04d9e1_8879.zip。
60.s1.2：对vv和vh极化雷达影像进行应用轨道文件校正处理，替换元数据文件(.xml)中的哨兵1号卫星轨道状态信息；
61.s1.3：对vv和vh极化雷达影像进行热噪声去除处理，消除sar卫星装置(例如，发射机、功率放大器、接收机)内部的热损耗，提高sar影像的信噪比；
62.s1.4：对vv和vh极化雷达影像进行斑点滤波处理，使用目前较为常用的refined lee相干斑滤波器，去除影像相干斑噪声；
63.s1.5：对vv和vh极化雷达影像进行地理编码处理，使用range-dopplerterrain correction(距离多普勒法地形校正)，赋予影像实际的坐标信息；
64.s1.6：对vv和vh极化雷达影像进行分贝化处理，将后向散射系数σ执行对数变换，即σ(db)＝10*lg(σ)，得到vv和vh两种极化方式的分贝化雷达后向散射系数，分别记为σ
vv
和σ
vh
。
65.s2：绘制不同极化条件下的后向散射系数σ的统计直方图，根据统计直方图的双峰特征，依据预设的自动阈值提取算法，定位到双峰之间的谷底位置，提取局部最小值作为原极化雷达影像的水体阈值；
66.在具体实施过程中，步骤s2包括：
67.s2.1：对vv和vh极化雷达影像进行中值滤波，达到去除椒盐噪声、平滑影像的效果，得到极化雷达影像后向散射系数σ
vv
和σ
vh
的统计直方图，如图4-5所示；本实施例中，中值滤波的滤波窗口大小为5
×
5像元。
68.s2.2：按后向散射系数由小到大的顺序，取每5个区间(bin)的频数最大值，设该最大值为因变量y，对应bin所在的后向散射系数为自变量x，构建拟合函数f(x,y)，以达到平滑原始直方图数据中的频数曲线，降低局部bin之间的锯齿扰动的效果；
69.s2.3：根据拟合函数f(x,y)的一阶导数判断收敛性，找到频数的局部最小值y0，对应的bin即为直方图双峰之间的谷底，进而得bin所在的后向散射系数作为提取原始vv和vh极化雷达影像的水体阈值；本实施例中σ
vv
的水体阈值t
vv
为-18.42db，σ
vh
的水体阈值t
vh
为-23.62db。
70.s3：计算不同极化雷达影像水体阈值之差的均值，取该均值的绝对值作为阈值偏移量t
δ
；
71.在具体实施过程中，步骤s3中，计算vv和vh极化雷达影像的阈值之差，取其绝对值作为阈值偏移量t
δ
。本实施例为t
δ
＝5.2。
72.s4：针对双峰特征较弱的极化方式，在其后向散射系数σ上叠加阈值偏移量t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠；
73.在具体实施过程中，步骤s4中，针对双峰特征较弱的vv极化雷达影像，每个像元均减去阈值偏移量t
δ
，即σ’vv
＝σ
vv-t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠。修改后的vv和原始vv、vh的频数分布对比如图6所示。
74.s5：比较每个像元在不同极化方式中的后向散射系数大小，取最小值构成一景新影像σ
min
；
75.在具体实施过程中，步骤s5中，比较每个像元的σ’vv
和σ
vh
大小，取每个像元σ’vv
和σ
vh
的最小值构成一景新影像，即全域最小值影像σ
min
。全域最小值影像σ
min
如图7所示。
76.s6：构建水体提取指数，并根据步骤s2的方法计算水体阈值，提取图幅范围内的水体像元。
77.在具体实施过程中，步骤s6中，包括：
78.s6.1：构建水体提取指数pwi＝lg(σ
min
×
σ
min
)，其直方分布如图8所示。
79.s6.2：按步骤s2所述的方法提取水体阈值t
wi
；本实施例中得到的水体阈值tw＝6.384。
80.s6.1：按pwi》t
wi
提取水体范围，即对于每个像元，其pwi》t
wi
则该像元属于水体范围。该图9-10所示。可以看出，本方法在保留主干河流、水库水面提取效果的同时，可有效的消除山体阴影的干扰，提高极化雷达影像的水体提取精度。
81.实施例2
82.本实施例提供一种基于极化雷达自身影像特征的水体提取系统，包括存储器和处理器，所述存储器中包括基于极化雷达自身影像特征的水体提取方法程序，所述基于极化雷达自身影像特征的水体提取方法程序被所述处理器执行时实现如下步骤：
83.s1：根据需要开展水体提取的时间、地点，选择相应的极化雷达影像并完成预处理，得到不同极化条件下的后向散射系数σ(单位：db)；
84.s2：绘制不同极化条件下的后向散射系数σ的统计直方图，根据统计直方图的双峰特征，依据预设的自动阈值提取算法，定位到双峰之间的谷底位置，提取局部最小值作为原极化雷达影像的水体阈值；
85.s3：计算不同极化雷达影像水体阈值之差的均值，取该均值的绝对值作为阈值偏移量t
δ
；
86.s4：针对双峰特征较弱的极化方式，在其后向散射系数σ上叠加阈值偏移量t
δ
，使不同极化雷达影像的后向散射系数直方图在局部最小值处重叠；
87.s5：比较每个像元在不同极化方式中的后向散射系数大小，取最小值构成一景新影像σ
min
；
88.s6：构建水体提取指数，并根据步骤s2的方法计算水体阈值，提取图幅范围内的水体像元。
89.在本技术所提供的实施例中，应该理解到，所揭露的设备和方法，可以通过其它的方式实现。以上所描述的设备实施例仅仅是示意性的，例如，所述单元的划分，仅仅为一种逻辑功能划分，实际实现时可以有另外的划分方式，如：多个单元或组件可以结合，或可以集成到另一个系统，或一些特征可以忽略，或不执行。另外，所显示或讨论的各组成部分相互之间的耦合、或直接耦合、或通信连接可以是通过一些接口，设备或单元的间接耦合或通信连接，可以是电性的、机械的或其它形式的。
90.上述作为分离部件说明的单元可以是、或也可以不是物理上分开的，作为单元显示的部件可以是、或也可以不是物理单元；既可以位于一个地方，也可以分布到多个网络单元上；可以根据实际的需要选择其中的部分或全部单元来实现本实施例方案的目的。
91.另外，在本发明各实施例中的各功能单元可以全部集成在一个处理单元中，也可以是各单元分别单独作为一个单元，也可以两个或两个以上单元集成在一个单元中；上述集成的单元既可以采用硬件的形式实现，也可以采用硬件加软件功能单元的形式实现。
92.本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：移动存储设备、只读存储器(rom，read-only memory)、随机存取存储器(ram，random access memory)、磁碟或者光盘等各种可以存储程序代码的介质。
93.或者，本发明上述集成的单元如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明实施例的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机、服务器、或者网络设备等)执行本发明各个实施例所述方法的全部或部分。而前述的存储介质包括：移动存储设备、rom、ram、磁碟或者光盘等各种可以存储程序代码的介质。