山区洪涝双极化SAR快速检测方法及系统与流程

本发明涉及一种山区洪涝双极化sar快速检测方法及系统，利用双极化sar影像、数字高程模型和水文地理数据，对山区洪涝灾害进行快速检测并抑制山体阴影误判，属于遥感灾害监测。

背景技术：

1、山区sar洪涝检测中，由于sar侧视成像机制影响，山体背坡、峡谷区域及局部高程遮挡区域容易形成雷达阴影。该类区域在sar影像中通常表现为低后向散射特征，与真实洪水区域具有相似的灰度响应，容易导致洪涝误判。

2、现有技术通常采用后向散射阈值、双极化特征和dem坡度约束进行洪水提取，但对于复杂山区地形，仅采用散射强度和简单坡度约束难以准确区分真实洪水与稳定地形阴影。

3、部分方法引入dem进行地形修正，但未充分结合sar成像几何关系、历史稳定低散射特征以及洪水空间水文连通规律，因此在山区环境下仍容易出现阴影误判率高、洪水边界不连续以及检测结果物理合理性不足的问题。

4、因此，需要一种能够结合sar成像几何、历史稳定低散射特征以及水文连通关系的山区洪涝检测方法，以降低山体阴影误判，提高山区洪涝检测精度。

技术实现思路

1、本发明的目的在于提供一种山区洪涝双极化sar快速检测方法，用于解决山区sar洪涝监测中山体阴影与真实洪水区域容易混淆的问题。

2、本发明采取的技术方案为：

3、山区洪涝双极化sar快速检测方法，包括步骤如下：

4、s1.获取目标山区多源基础数据并进行空间统一处理，多源基础数据包括历史水体数据、数字高程模型dem、水文地理数据、三期灾前sar时序影像、灾前vv及vh双极化sar影像、灾中vv及vh双极化sar影像；

5、s2.对灾前和灾中双极化sar影像进行预处理，提取灾前和灾中vv、vh双极化sar影像的vv同极化后向散射特征、vh交叉极化后向散射特征、vv-vh极化比值特征和vv-vh极化差异特征；

6、s3.利用灾前、灾中两期影像的后向散射强度差异、双极化比值差异和归一化差异特征，提取灾中疑似新增低散射区域，形成初始洪水候选区；

7、s4.构建雷达几何阴影概率模型：在获得初始洪水候选区后，基于数字高程模型dem和sar成像几何参数，对目标山区的雷达阴影、背坡弱照射和局部地形遮挡现象进行定量建模，计算每个像元的雷达几何阴影概率，用于识别低散射候选像元中由地形阴影引起的误判区域，将雷达几何阴影概率与初始洪水候选区进行空间叠加分析生成雷达几何阴影概率图；

8、s5.构建历史稳定低散射阴影先验模型：在获得雷达几何阴影概率图的基础上，进一步利用三期灾前vv、vh双极化sar时序影像，对目标山区长期稳定存在的低后向散射区域进行统计分析，构建历史稳定低散射先验图，用于区分灾中新增洪涝淹没区域与长期存在的山体阴影、深切沟谷、背坡遮蔽区、稳定水体及永久性低散射地貌区域；

9、s6.构建水文连通性约束模型：在获得初始洪水候选区、雷达几何阴影概率图以及历史稳定低散射阴影先验图的基础上，利用水文地理数据，对洪水候选像元的空间分布合理性进行约束，计算每个洪水候选像元的水文连通性得分，用于判断该低散射像元是否符合真实洪涝扩展的地形水文过程；

10、s7.融合多源约束因子计算洪水置信度：将双极化sar水体散射特征、sar成像几何约束、历史稳定低散射先验和地形水文连通约束进行综合融合，计算每个洪水候选像元的像元级洪水置信度，用于定量表征该像元属于真实洪涝淹没区域的可能性；

11、s8.阴影误判抑制与洪涝范围优化：对初始洪水候选区进行阴影误判抑制、连通域分析和空间一致性优化，得到空间连续性强、物理合理性高的山区洪涝灾害范围。

12、上述方法中，步骤s1所述的灾前vv及vh双极化sar影像和灾中vv及vh双极化sar影像用于提取洪涝发生前后地表后向散射变化信息；所述双极化sar影像限定为vv、vh组合，vv为同极化通道，vh为交叉极化通道；所述数字高程模型dem用于提取坡度、坡向、高程起伏和局部地形遮挡关系；所述水文地理数据包括河网、水库、沟谷、低洼区、汇流路径和流域边界；所述历史水体数据用于表征目标山区长期稳定水体的空间分布；所述三期灾前sar时序影像用于识别长期稳定存在的低散射区域。

13、步骤s2中所述的预处理为：对灾前、灾中vv及vh双极化sar影像先通过轨道精校正提升影像地理定位精度，再经去热噪声、辐射定标处理，将影像灰度值转为分贝形式的后向散射系数；随后采用像元窗口的refined lee自适应滤波去除斑点噪声，兼顾影像细节与降噪效果；结合dem数据完成地形与正射校正，修正几何畸变且不剔除雷达阴影，通过初配准与精配准结合的方式实现影像精准匹配，最后开展辐射归一化处理，消除成像条件干扰。对于三期灾前sar影像，也采用相同的预处理流程。

14、步骤s2中分别计算vv、vh后向散射变化量、双极化比值变化量、双极化差异变化量和低散射变化幅度，各特征按下列公式组计算：

15、，

16、其中，表示极化通道p的灾前灾中后向散射变化量；表示灾前vv及vh双极化sar影像中极化通道p的后向散射系数；表示灾中vv及vh双极化sar影像中极化通道p的后向散射系数；p表示sar影像的极化通道，所述极化通道为vv和vh；

17、，

18、δr()=rafter()-rbefore()，

19、其中，t表示灾前时相before或灾中时相after；表示像元在时相t的vv/vh极化比值；表示像元在时相t的vv通道后向散射系数，表示像元在时相t的vh通道后向散射系数，δr()表示灾前、灾中极化比值变化量；ε为防止分母为零的极小正数。

20、dt(x)=σvv,t(x)-σvh,t(x)，

21、δd(x)=dafter(x)-dbefore(x)，

22、其中，dt(x)表示像元x在时相t的vv-vh极化差异；δd(x)表示灾前、灾中极化差异变化量。，

23、其中，alow(x)表示像元x的低散射变化幅度，其值越大，表示vv和vh通道在灾中相对于灾前的综合下降越明显。

24、同时，构建双极化归一化差异指数，用于表征同极化通道与交叉极化通道之间的相对散射差异，其计算公式为：

25、，

26、其中，表示像元x在时相t的双极化归一化差异指数；表示vv同极化通道后向散射系数；表示vh交叉极化通道后向散射系数；ε为防止分母为零的极小正数。

27、步骤s3中所述的洪水候选像元的判定表示为：

28、，

29、其中，cf(x)表示像元x的初始洪水候选标记；σp,after(x)、δσp(x)、δr(x)、δd(x)、alow(x)和ndpiafter(x)均为步骤s2的输出特征；、和分别表示空间分区z内的水体散射阈值、变化量阈值和低散射变化幅度阈值；、和分别表示水体对应的极化比值变化、极化差异变化和归一化差异指数取值范围，各阈值及取值范围在步骤s3中根据步骤s2输出特征的分区统计结果确定。

30、步骤s4中，基于目标像元x的坡度α(x)、坡向β(x)、sar入射角θ(x)、sar观测方位角φ(x)和雷达视线方向高程剖面，计算归一化坡度阴影因子s_α(x)、坡向背离因子s_β(x)、局部入射角因子s_θ(x)和局部高程遮挡因子s_h(x)，并按加权组合计算雷达几何阴影概率：

31、，

32、其中，pshadow(x)表示像元x的雷达几何阴影概率；wα、wβ、wθ、wh为非负权重，且wα+wβ+wθ+wh=1。

33、所述各归一化因子均由dem和sar成像几何参数直接计算，具体计算方式在实施例1的步骤s4中说明。

34、步骤s5中对多期灾前sar影像中每一像元的后向散射强度进行时序统计分析，设目标区域共有n期灾前sar影像，像元x在第i期灾前sar影像中极化通道的后向散射系数为，则该像元在极化通道上的历史后向散射均值可表示为：

35、 ，

36、其中，表示像元x在极化通道上的历史后向散射均值；(x)表示像元x在第i期灾前sar影像中极化通道的后向散射系数；n表示灾前sar影像的时相数量；表示vv和vh两个极化通道；

37、同时，计算该像元在多期灾前sar影像中的后向散射方差，用于评价其低散射状态的时序稳定性，其计算公式为：

38、 ，

39、其中，表示像元x在极化通道上的历史后向散射方差，若某一像元的历史后向散射均值较低，且方差较小，则说明该像元在多期灾前影像中长期保持稳定低散射状态，其低散射特征更可能来源于稳定地形阴影、深切沟谷或永久低散射地貌，而非灾中新增洪涝水体；

40、统计像元在多期灾前影像中被判定为低散射状态的出现频率，设低散射判定阈值为，则像元x的历史低散射出现频率表示为：

41、，

42、其中，i[·]为指示函数；括号内条件成立时i[·]=1，否则i[·]=0。

43、当某一像元同时满足历史后向散射均值低、历史后向散射方差小、低散射出现频率高，并且其空间位置与雷达几何阴影概率较高区域、陡坡背坡区域和局部高程遮挡区域一致时，将该像元判定为历史稳定低散射阴影像元，其判定条件表示为：

44、，

45、其中，表示像元x的历史稳定低散射阴影标记；当=1时，表示该像元属于历史稳定低散射阴影像元；当=0时，表示该像元不属于历史稳定低散射阴影像元；表示历史后向散射均值阈值；表示历史后向散射方差阈值；表示历史低散射出现频率阈值；表示像元x的雷达几何阴影概率；表示阴影概率判定阈值。

46、为获得所需的连续历史稳定低散射阴影概率，对历史低均值、低方差和低散射频率进行归一化，并与雷达几何阴影概率融合：

47、，

48、其中，表示历史低均值得分；表示历史低方差得分；表示vv、vh通道历史低散射出现频率的平均值；、、、为非负权重且其和为1。像元x的历史稳定低散射阴影概率经截断后取值范围为[0，1]；历史水体像元不作为历史阴影像元，并将其)置为0。

49、步骤s6水文连通性得分的计算过程为：

50、首先将河网、水库、历史水体、沟谷线、低洼区、流向数据、汇流累积量数据和子流域边界数据统一至与sar影像一致的空间参考系统、像元尺度和研究区范围，对于矢量形式的河网、水库、沟谷线和流域边界数据，将其栅格化至sar影像对应的空间分辨率；对于由dem计算得到的流向、汇流累积量和低洼区数据，则与sar影像进行空间配准和尺度统一，保证水文因子与洪水候选像元之间具有一一对应的空间关系；

51、然后，基于统一后的水文地理数据，分别计算洪水候选像元与水文通道之间的距离关系、与历史水体和水库之间的连通关系、与低洼区之间的重合关系、与汇流路径之间的一致性以及其所在位置的汇流累积强度；

52、其中，洪水候选像元与河网之间的距离用于衡量该像元是否靠近主要水文通道；设洪水候选像元x到最近河网和水体边界的距离为，则距离越小，表示该像元与河流、水库和历史水体的空间关联越强，其属于洪涝淹没区域的可能性越高，采用距离衰减函数将其转换为河网邻近得分：

53、，

54、其中，表示像元x的河网邻近得分；表示像元x至最近河网、水库和历史水体边界的距离；表示距离衰减尺度参数；

55、计算洪水候选像元与低洼区之间的空间重合关系，低洼区通常是洪水汇集和滞留的重要区域，若某一低散射候选像元位于低洼区内，则其符合洪涝淹没的地形条件，设表示像元x的低洼区重合因子，当像元x位于低洼区和低高程汇流区域内时，取高值；当像元x位于陡坡、山脊和非汇流区域时，取低值；

56、计算洪水候选像元与汇流路径之间的方向一致性，由dem计算地表流向和汇流路径，真实洪水扩散方向通常与地表汇流方向和河谷延伸方向具有较高一致性，设像元x处的局部洪水候选斑块延伸方向为，dem推导的汇流方向为，则方向一致性得分表示为：

57、，

58、其中，表示像元x的汇流方向一致性得分；表示洪水候选斑块在像元x附近的空间延伸方向；表示dem推导的地表汇流方向；

59、同时，利用汇流累积量表征像元所在位置的集水能力，汇流累积量越高，说明该位置接收上游径流汇入的能力越强，发生洪涝积水和漫溢的可能性越高，设像元x的汇流累积量为，对其进行归一化处理后得到汇流累积得分：

60、，

61、其中，表示像元x的汇流累积得分；表示像元x的汇流累积量；和分别表示研究区内汇流累积量的最小值和最大值；为防止分母为零的极小正数；

62、此外，计算洪水候选像元与历史水体、水库、河网和沟谷之间的连通关系，设c(x)表示像元x所在候选斑块的水文连通因子，当该斑块与河网、水库、历史水体和沟谷连通时，c(x)取高值；当该斑块孤立分布且不具备水文连通关系时，c(x)取低值；

63、综合上述水文因子，构建像元级水文连通性得分，所述水文连通性得分表示为：

64、，

65、其中，表示像元x的水文连通性得分；表示河网邻近得分；表示低洼区重合因子；表示汇流方向一致性得分；表示汇流累积得分；表示水文连通因子；、、、、分别表示各水文因子的权重系数，并满足：

66、，

67、其中，各权重系数根据研究区地形类型、洪涝灾害类型、河网密度、dem精度和历史灾害样本进行设定，并通过历史洪涝样本和少量人工核查样本进行自适应校准。

68、步骤s7中像元级洪水置信度表示为：

69、，

70、其中，pflood(x)表示像元x的洪水置信度；pwater(x)表示水体散射得分；phydro(x)表示水文连通性得分；pshadow(x)表示步骤s4得到的雷达几何阴影概率；phistory(x)表示步骤s5得到的历史稳定低散射阴影概率；λ1、λ2、λ3、λ4分别表示上述四项的权重系数。

71、步骤s8中所述阴影误判剔除规则表示为：

72、，

73、其中，mremove(x)表示像元x的剔除标记；当mremove(x)=1时，将该像元判定为阴影误判像元并从洪水候选区中剔除；表示归一化后的洪水置信度；pshadow(x)、phistory(x)和phydro(x)分别表示雷达几何阴影概率、历史稳定低散射阴影概率和水文连通性得分；tflood、tshadow、thistory和thydro分别表示洪水判定阈值、阴影概率阈值、历史阴影阈值和水文连通性阈值。

74、山区洪涝双极化sar快速检测系统，包括数据获取模块、sar预处理模块、洪水候选区提取模块、雷达几何阴影概率计算模块、历史稳定阴影先验构建模块、水文连通性计算模块、洪水置信度计算模块、阴影误判抑制模块、结果输出模块，

75、数据获取模块，用于获取目标山区多源基础数据并进行空间统一处理，多源基础数据包括历史水体数据、数字高程模型dem、水文地理数据、三期灾前sar时序影像、灾前vv及vh双极化sar影像、灾中vv及vh双极化sar影像；

76、sar预处理模块，用于对灾前和灾中双极化sar影像进行预处理，提取灾前和灾中vv、vh双极化sar影像的vv同极化后向散射特征、vh交叉极化后向散射特征、vv-vh极化比值特征和vv-vh极化差异特征；

77、洪水候选区提取模块，用于利用灾前、灾中两期影像的后向散射强度差异、双极化比值差异和归一化差异特征，提取灾中疑似新增低散射区域，形成初始洪水候选区；

78、雷达几何阴影概率计算模块，用于在获得初始洪水候选区后，基于数字高程模型dem和sar成像几何参数，对目标山区的雷达阴影、背坡弱照射和局部地形遮挡现象进行定量建模，计算每个像元的雷达几何阴影概率，用于识别低散射候选像元中由地形阴影引起的误判区域，将雷达几何阴影概率与初始洪水候选区进行空间叠加分析生成雷达几何阴影概率图；

79、历史稳定阴影先验构建模块，用于在获得雷达几何阴影概率图的基础上，进一步利用三期灾前vv、vh双极化sar时序影像，对目标山区长期稳定存在的低后向散射区域进行统计分析，构建历史稳定低散射先验图，区分灾中新增洪涝淹没区域与长期存在的山体阴影、深切沟谷、背坡遮蔽区、稳定水体及永久性低散射地貌区域；

80、水文连通性计算模块，用于在获得初始洪水候选区、雷达几何阴影概率图以及历史稳定低散射阴影先验图的基础上，利用水文地理数据，对洪水候选像元的空间分布合理性进行约束，计算每个洪水候选像元的水文连通性得分，判断该低散射像元是否符合真实洪涝扩展的地形水文过程；

81、洪水置信度计算模块，用于将双极化sar水体散射特征、sar成像几何约束、历史稳定低散射先验和地形水文连通约束进行综合融合，计算每个洪水候选像元的像元级洪水置信度，定量表征该像元属于真实洪涝淹没区域的可能性；

82、阴影误判抑制模块，用于对初始洪水候选区进行阴影误判抑制、连通域分析和空间一致性优化，得到空间连续性强、物理合理性高的山区洪涝灾害范围；

83、结果输出模块，用于输出山区洪涝灾害检测结果及其不确定性表达结果。

84、本发明的有益效果是：

85、（1）通过灾前、灾中vv和vh双极化后向散射变化、极化比值变化、极化差异变化及归一化差异指数联合构建初始洪水候选区，相较于单一后向散射阈值方法，能够增强灾中新增水体与其他低散射地物的可分性。

86、（2）将坡度、坡向背离程度、局部入射角和雷达视线高程遮挡关系转换为可直接计算的雷达几何阴影概率，并结合灾前时序影像形成连续的历史稳定低散射阴影概率，能够从成像机理和时间稳定性两个维度识别山体阴影、深切沟谷及永久低散射地貌。

87、（3）利用河网邻近、低洼区重合、汇流方向一致性、汇流累积量及水文连通关系形成正向水文约束，再与水体散射得分、几何阴影概率和历史阴影概率融合，可避免将水文上不合理的孤立暗斑直接判为洪水，提高洪涝范围的物理合理性和空间连续性。

88、（4）上述双极化变化特征提供灾中新增水体证据，几何阴影概率和历史阴影概率提供互补的反向证据，水文连通性提供正向过程约束；各创新点在像元级置信度融合及斑块级空间优化中协同作用，在不依赖光学影像的条件下兼顾快速检测、阴影误判抑制、边界连续性和结果可解释性。