洪涝灾害是近几十年来最常见和最具破坏力的自然灾害之一^[1,2]。受气候变化的影响,洪涝灾害的规模和频率正呈持续上升趋势,未来洪涝事件可能会更加突发和严重^[3]。中国是遭受洪涝灾害影响最大的国家之一,洪涝灾害可对人民生命安全、经济状况以及基础设施造成重大损失^[2,4]。在2023年8月初,受台风“杜苏芮”北上的影响,我国华北、东北等地遭受强降雨,造成了严重的洪涝灾害。2023年8月11日台风“卡努”在辽宁省登陆,东北多地再次遭受强降雨,辽宁、黑龙江等地共8.8万人不同程度受灾,农作物受灾面积21.5千hm²,直接经济损失9.2亿元。因此,及时准确地监测洪涝灾害并绘制洪涝淹没地图对应急救援和灾后评估至关重要。

洪涝灾害期间恶劣的天气条件限制了光学遥感在洪涝监测中的有效性^[5]。合成孔径雷达(synthetic aperture radar, SAR)卫星具有全天时、全天候获取地面信息的能力,在云雾和黑夜条件下仍能获取地面洪涝信息,使其在时效性要求高的洪涝监测中发挥了重要作用^[6]。洪涝事件具有突发性,通常持续几天或十几天时间,期间的洪灾救援与决策支持需要高时效性的洪涝淹没动态监测结果,这就要求SAR影像获取间隔不宜过长^[4]。然而,SAR卫星通常具有固定的重访周期,以Sentinel-1为例,其在中国区域的重访周期为12 d,这导致其在洪涝事件期间通常仅能获取1幅SAR影像,有时甚至无法获取洪涝期间的SAR影像^[7]。随着在轨SAR卫星的增加,利用不同卫星获取的多源SAR影像进行洪涝淹没动态监测,成为解决洪涝期间高时效性洪涝淹没动态监测结果获取问题的可行思路。但是,不同SAR卫星影像的特征空间存在差异,导致多源SAR数据间不具有直接可比性。

从差异特征提取的角度,现有基于SAR数据的洪涝淹没制图方法主要可分为2类: 分类后比较方法和多时相SAR影像变化检测方法^[8]。其中,分类后比较法在洪涝淹没范围提取中被广泛使用,既可用于基于同源数据的洪涝监测,又可用于基于多源数据的洪涝监测。该类方法可看作一项二分类问题,即将单时相SAR影像分为水体和非水体。郭山川等^[9]依托GEE云计算平台,通过融合Sentinel-1 SAR和Sentinel-2光学影像提取水体,并在长江中下游地区洪涝淹没监测中进行了应用。张文璇等^[10]通过拟合Sentinel-1 SAR影像的概率密度分布,提出了一种大区域洪涝监测方法。由于水体在SAR影像上的反射特性不同于其他地物,SAR影像的直方图通常显示出双峰特性,因此,大多数研究基于OTSU,Kittler-Illingworth (KI)和MaxEntropy等全局阈值选择方法,利用单一阈值分割SAR影像得到水体和非水体分类结果。然而,针对大面积的水域,使用单阈值处理整幅SAR影像往往难以达到理想的效果,这是由于大面积水域的环境变化多样,导致不同区域的SAR影像特征差异较大,具有空间异质性。近年来,基于多时相SAR影像变化检测的方法在洪涝灾害监测中得到大量学者的关注,其把洪涝监测看作地表水体的变化问题^[11]。Twele等^[12]的研究结果表明,与基于SAR影像分类后比较的洪涝淹没提取方法相比,基于变化检测的方法在抑制与水体有相似后向散射特征的地物方面具有显著优势。冷英等^[13]把模糊聚类、最近距离聚类与变化检测的对数比值法相结合,进行淹没水体提取。黄平平等^[11]在差异图融合的基础上,提出了结合无监督聚类的洪涝淹没区域检测方法,并利用Sentinel-1数据提取了2020年鄱阳湖区域的洪涝淹没范围。然而,这些像素级的变化检测方法对斑点噪声的抑制能力较弱,而聚类使得洪涝淹没结果丢失了宝贵的边缘细节。Dong等^[14]评估了3种经典深度学习模型(BiT,STANets和SNUNet)在Sentinel-1 SAR影像中提取淹没水体的能力,结果表明BiT模型具有最佳性能。Zhao等^[15]提出了一种基于跨任务迁移学习策略的孪生网络(Siam-DWENet),实现了利用少量样本进行洪涝淹没检测。但是,基于深度学习的洪涝淹没检测方法往往需要高质量样本进行训练,在高时效性洪涝灾害监测中有一定的局限性。

综上所述,分类后比较法在大范围洪涝淹没监测中的精度不理想; 而基于变化检测的洪涝监测方法,通常使用同源SAR数据,无法直接用于多源SAR数据协同的洪涝淹没时序动态监测研究。为此,本文提出了一种基于多源时序SAR数据的洪涝淹没动态监测方法。首先,针对高分三号(GF-3)、Sentinel-1和涪城一号SAR影像特征差异明显问题,提出了基于后向散射特征的多源SAR数据特征空间对齐方法; 然后,引入渐进非局部理论(progressive non-local information,PNLI)从多时相SAR影像中提取洪水淹没差异图,进而通过自动阈值分割得到洪水淹没图; 最后,对2023年8月河北省涿州市及周边发生的洪涝灾害进行动态监测,并对涿州市洪涝淹没的动态变化特征进行分析。

涿州市位于河北省中部、保定市北部,其境内的拒马河是京冀两地重要的行洪河道。拒马河可分为南拒马河和北拒马河。其中,北拒马河自涿州开始向西南流经琉璃河后再入白沟河,而南拒马河则从涿州市南部起始,经涞水和定兴2县,最终汇入易水河。在台风“杜苏芮”和冷暖空气的共同影响下,自2023年7月29日起,河北涿州市遭遇连续强降雨,叠加上游洪水过境影响,该区域遭受了严重的洪涝灾害。本文研究区如图1所示,包括河北省的涿州市、固安县、高碑店市、定兴县、涞水县和北京市的房山区,图1显示了7月24日由Sentinel-1获取的洪水发生前SAR影像,选取红色矩形和蓝色矩形区域分别作为本文方法的验证集Ⅰ、验证集Ⅱ。

本研究使用由Sentinel-1,GF-3和涪城一号获取的5幅时间序列SAR影像进行涿州市洪涝淹没动态监测研究,数据基本情况如表1所示。

1)Sentinel-1数据。Sentinel-1是欧空局发射的C波段SAR卫星,具有4种成像模式,在中国区域的重访周期为12 d。本研究使用2023年7月24日、8月5日和8月17日在干涉宽条带(interferometric wide swath,IW)模式下采集的Sentinel-1 GRD数据。GRD产品提供具有10 m高空间分辨率的SAR数据,该产品相对于SLC级数据而言经过了多视处理。Sentinel-1的各级别产品免费提供,其中GRD产品被广泛应用于洪涝灾害监测。利用欧空局提供的软件SNAP可实现Sentinel-1数据的自动化预处理,主要包括6部分: ①轨道校正; ②热噪声去除; ③辐射定标; ④Refined Lee滤波; ⑤地形校正; ⑥分贝化。

2)GF-3数据。GF-3卫星是我国首颗高分辨率C波段多极化SAR卫星,共有12种成像模式^[16]。本文采用2023年8月1日在精细条带2模式(FSII)下采集的L1级SAR影像。通过ENVI SARscape5.6进行了预处理: ①数据导入; ②多视; ③Refined Lee滤波; ④辐射定标; ⑤地形校正; ⑥分贝化。GF-3号经过地形校正后的空间分辨率为10 m。

3)涪城一号数据。涪城一号是首颗“四川造”C波段轻小型SAR卫星,具备全天时全天候的差分干涉SAR数据获取能力,采用单星重轨干涉技术,能够稳定提供高频次的毫米级地面形变监测服务。本文使用在2023年8月3日在条带模式下采集,经过多视、辐射定标和地形矫正的L2级别的涪城一号SAR影像,并对其应用Refined Lee滤波消除斑点噪声,其空间分辨率为10 m。

为了验证本文方法的有效性和鲁棒性,选取图1中的验证集Ⅰ和验证集Ⅱ进行验证。2个验证集分别由不同时相的多源SAR影像对组成(图2)。为了进行定量验证,通过目视解译获取了参考洪涝淹没区域,如图2(c)和图2(f)所示,其中洪涝淹没区域标注为白色。

本文提出的洪涝淹没动态监测方法主要包括5个步骤: ①多源SAR数据预处理; ②基于后向散射特征对多源SAR数据进行特征空间对齐; ③基于PNLI的洪水淹没差异提取; ④自动阈值分割得到洪水淹没图; ⑤针对双时相淹没图进行逻辑运算得到洪涝动态监测结果。本文方法的基本过程如图3所示。

对于不同SAR传感器在不同时间所获取SAR影像,其后向散射系数存在较大的时空异质性,通常难以进行直接比较。多源光学遥感影像处理中也存在类似问题,现有研究通常采用基于色彩空间变换的方法来减小成像差异^[17]。然而,光学与SAR影像的成像机制不同,无法将光学影像的处理方法简单应用于SAR强度影像。为此,本文顾及SAR影像成像机制,首先把多源SAR卫星获取的复数影像转换到后向散射特征空间,进而将前一时相SAR影像的后向散射特征空间迁移到后一时相SAR影像中,以减少洪涝监测中多源SAR影像的成像条件影响。具体来说,首先,将SAR影像强度特征转换为以分贝表示的后向散射特征,公式为:

其中,σ₀(dB)为以分贝为单位的后向散射系数影像; I₀为待转换的强度影像。然后,计算某t时相SAR影像I_t的方差S_t和平均值Gt¯,其中,t=1,2(t为1时表示前一时相,t为2时表示后一时相)。转换Z_t和缩放H_t参数计算公式如下:

最后利用这些参数实现前一时相SAR影像和后一时相SAR影像特征空间对齐,对齐的影像I的计算公式为:

基于后向散射特征的多源SAR数据特征空间对齐结果如表2所示。

基于分类比较的洪水淹没检测方法过度依赖于单时相水体的提取精度,直接从多时相SAR影像中提取洪水淹没信息可以有效提高数据处理的自动化^[4]。本文引入从多时相SAR影像的时空邻域中提取PNLI方法,生成洪水差异影像^[18]。PNLI利用基于比值的高斯核函数来表征非局部邻域内不同像元的相似性,可以有效减少SAR影像固有的相干斑噪声。此外,PNLI利用渐进优化的策略来进一步提高洪涝SAR影像中时空非局部信息提取的可靠性。这些特点使得PNLI在洪涝监测中可以有效减少外界环境对SAR影像的影响。PNLI方法主要分为2个阶段: ①利用渐进优化策略得到匹配窗口中的最优同质区域,提高非局部邻域内像元相似性表征的可靠性; ②提取时空非局部信息并生成表征洪水淹没的差异图像。

第一阶段: 首先,设定初始搜索窗口的大小与匹配窗口相同; 其次,采用标准差为σ的各向同性高斯模型来确定匹配窗口内位置i处像元的权重M(i); 最后,计算初始搜索窗口(即匹配窗口)内像元的权重w'_t,公式为:

式中: w'_t(i,x)为某一时相t的影像中位置i处的像元与位置x处像元的相似性权重; R为V(i)和V(x)之间的比值向量; V(i)为中心像元位置为i、大小等于匹配窗口的邻域内像元组成的向量; V(x)为中心像元位置为x、大小等于匹配窗口的邻域内像元组成的向量,可以参考文献[18]计算; a为高斯核的标准差; Z'_t (x)为归一化参数,其中Ω'_x为影像I_t中位置x处像元的初始搜索窗口组成的邻域。

PNLI基于w'_t得到同质性区域不规则分布的匹配窗口权重w_t,公式为:

式中: w_t(i)为影像I_t中匹配窗口内位置i处像元用于计算相似性的权重; Z'为归一化参数; w″_t为相似性较大的同质性区域二值掩模图像,可参考文献[18]计算; M(i)为式(5)中匹配窗口内位置i处像元用于计算相似性的权重。

第二阶段: 首先,利用优化后的匹配窗口权重影像和确定的搜索窗口计算差异影像,公式为:

式中: D为使用PNLI得到的差异影像; u₁(x),u₂(x)分别为影像I₁与I₂中位置x处的PNLI。在差异影像D中,非淹没区像元具有较小的强度值,淹没区具有较大的强度值。其中u₁(x)计算公式为:

式中: I(t,i)为SAR影像I_t中位置i处像元的强度值; Ω″为搜索窗口组成的时空邻域; w‴(t,i,1,x)为表征影像I_t中位置i处的像元与影像I₁中位置x处像元相似性的权重^[18]。

PNLI主要有噪声水平、匹配窗口和搜索窗口3个参数。为了避免用户经验输入噪声的水平的影响,本文将使用文献中提出的噪声水平自适应估计方法来确定噪声水平^[19]。搜索窗口大小使用现有研究推荐的设置: 21×21^[20]。匹配窗口采用基于邻域异质性的窗口优选方法自适应选择^[18]。

本文使用MaxEntropy自动化阈值分割方法获取洪水淹没图,其通过使分割后的前景和背景具有最大的不确定性和信息量来自动确定阈值^[21]。MaxEntropy分割方法的目的是找到一个合适的阈值,使得洪涝淹没范围和非淹没范围熵的和最大化。阈值为q时的熵计算公式为:

式中: p(m)为差异影像D中灰度值为m的概率; p₀(q),p₁(q)分别为q阈值下非淹没区和淹没区的累积概率; H₀(q)和H₁(q)分别为非淹没范围和淹没范围对应的熵; K为差异影像中灰度值的理论最大值,对于8位存储影像,K=255; H(q)为图像的总熵。

随着SAR卫星数量的增加,在发生洪涝等自然灾害时,可以获得更多的SAR影像数据。利用多源SAR影像可以提供高空间分辨率和逐天级别的观测能力,这对于监测洪涝的快速演变以及评估洪涝的影响范围和程度至关重要。因此,本文提出了一种基于变化检测结果图的洪涝动态监测方法。首先基于PNLI得到洪水淹没图; 然后将洪涝灾害中的变化像元赋值为1,不变化像元赋值为0; 最后将前一阶段的变化图和后一阶段的变化图通过以下逻辑运算得到洪涝动态变化结果:

式中: X₁和X₂分别为前后时间段的洪涝变化图; W={1,2,3,4}为洪涝淹没变化范围,1为退水区域,2为洪涝淹没区域,3为新增洪涝区域,4为其他区域。

为了验证本文方法的有效性,与其他7种先进的洪水淹没制图方法进行了定量比较。实验中采用的对比方法既包括对数比值法(log-ratio, LR)^[22]、均值比值法(mean-ratio, MR)^[23],PNLI和深度学习方法DDNet^[24]等基于变化检测的方法,也包括K-means^[25],OTSU^[26]和Kittler and Illingworth (KI)^[27]等在洪水淹没制图中广泛使用的分类后比较方法。为了定量分析本文方法的有效性,本文对验证集上的洪涝淹没检测结果计算了如下5个指标: Kappa系数、准确度(overall accuracy,OA),F1分数、精确度(precision,P)和召回率(recall,R)^[28]。

本文在验证集Ⅰ与Ⅱ中测试了上述方法的性能,得到了2个验证集的洪水淹没结果图(图4、图5),其中参考变化图为目视解译结果,图中白色为淹没区,黑色为未淹没区。从图4可以明显看出,直接从双时相SAR影像中提取洪水差异的结果优于分类后比较方法。图4(a)-(c)的结果与其他结果相比具有显著的差异,这是由于图2(a)所示的影像双峰特征不明显,分类结果中很多非水体区域被误分为水体,进而在分类后比较时把误差传递给洪水淹没结果。此外,LR的结果受斑点噪声的影响较大,MR虽然能较好地抑制斑点噪声的影响,但其对洪水淹没边缘信息丢失较为严重。从图4和图5中可以直观看出,相比本文方法,LR,MR和DDNet在蓝色矩形区域存在较多的误检,在绿色矩形区域存在较多的漏检。图4(h)和(i)与图5(h)和(i)中显示了从不同方向的特征空间对齐方法获取的洪水淹没图,可以直观看出不同方向的结果相似。相比于未进行特征空间对齐的PNLI方法,本文方法更好的抑制了成像差异和斑点噪声的影响。因此,本文方法既可以抑制SAR斑点噪声的影响,也可以保留更好的洪水淹没边缘细节。

从定量评估结果来看(表3),本文方法在整体评价指标(Kappa,F1和OA)中均优于其他方法。验证集Ⅰ主要包括建设用地,可以看出分割聚簇方法对建设用地的洪水淹没检测较为困难,K-means,OTSU和KI的F1均小于0.5; MR和DDNet虽然具有高的P(MR为0.912 9,DDNet为0.908 4),但它们具有较低的R(MR为0.724 6,DDNet为0.604 4),这导致它们遗漏了更多的洪水淹没区域。这是因为在洪涝应急制图中,P值越高,代表更少的非洪水淹没像素误判为洪水淹没像素; R值越高,代表更少的洪水淹没像素漏判为非洪水像素,因此确保更高的R可能比更高的P更为关键^[1]。另外,验证集Ⅱ中有大量植被区域,MR可以检测到大部分淹没区域,这是因为有效利用邻域信息抑制了斑点噪声的影响。此外,与未进行特征空间对齐的PNLI方法相比,本文方法具有更高的鲁棒性,在验证集Ⅰ上F1平均提升了2.1%,在验证集Ⅱ上F1平均提升了5.2%。不同方向特征空间对齐的定量评估结果相差不大,表明该方法对转换方向的灵敏度较低,即“T1到T2”和“T2到T1”的不同方向对齐对结果影响不大。综上所述,本文方法具有最优的性能,在2个验证集中F1均大于0.86,Kappa均大于0.84。性能优异的原因主要有2个方面: 首先是本文方法基于后向散射特征的多源SAR数据特征空间对齐,使它们具有相似的特征空间从而减小由于成像卫星不同而引起的差异; 其次是本文方法充分利用多时相SAR影像的非局部时空信息,减小了时空异质性的影响,提高了洪涝淹没检测的准确性和鲁棒性。

洪涝事件具有突发性,通常持续几天或十几天时间,期间的洪灾救援与决策支持需要高时效性的洪水动态监测地图而进行动态监测。利用不同卫星获取的多源SAR影像进行洪涝淹没动态监测,成为解决洪水期间高时效性洪涝淹没动态监测结果获取问题的可行思路。基于涿州市及周边区域发生洪涝灾害时的多源时序SAR影像,采用本文方法进行洪涝监测,获得了从2023年7月24日至2023年8月17日的洪涝淹没变化图(图6),其中灾前的水体分布情况是由专家对7月24日SAR影像进行目视解译而获取的。结果表明,截至8月1日(图6(a)),由于上游来水增多,拒马河和小清河周边大量区域被淹没,涿州市主城区西北部出现明显的洪涝区域。同时,白沟河下游也出现大量新增水体。截至8月3日(图6(b)),涿州市城区西部的水体逐渐退去,城区北部的水体略有褪去。而大量的洪涝向白沟河涌去,白沟河下游出现大量的新增洪涝,大面积的农田被淹没,多座村庄被洪水包围。截至8月5日(图6(c)),涿州市城区的洪涝基本退去,白沟河上游的洪涝也开始逐渐消退,下游的农田和村庄仍处于洪水之中。截至8月17日(图6(d)),涿州市和其周边区域的大面积洪涝已经退去。以上结果表明,基于本文方法可精确地描绘此次洪水事件的动态变化,表明该方法可为应急救援、决策评估和灾后评估及时提供宝贵信息。

洪涝灾害是威胁人类生命和财产的主要灾害之一。利用多种SAR影像及时进行淹没范围及其动态变化监测,在减少洪灾损失中发挥着重要作用。但是受不同卫星的传感器参数和成像条件等限制,非同源SAR影像成像方面存在差异,这是洪涝灾害动态监测面临的挑战之一。不同时间和区域的差异可导致检测到大量的伪淹没区,本文使用的后向散射特征空间校正方法,可以有效减轻时间序列图像之间的数据分布差异。此外,洪涝灾害往往以突发性的方式发生,极端情况下,可在短时间内快速形成。这种突发性使得变化检测方法能够将注意力集中在水体变化上,而无需过度关注其他类别的变化。因此,在洪涝灾害区域中,由于洪涝迅速泛滥和水体广泛扩展,可以将变化区域定义为水体的变化。通常,洪涝测绘使用比值法和对数比值法等像素级变化检测算法,但无法抑制斑点噪声影响。相较于仅使用单个图像进行水体识别并进行洪涝测绘,本文方法联合多时相SAR影像的时空非局部邻域提取洪涝信息,在复杂地区洪涝监测中具有明显的优势。

本研究针对突发洪涝灾害淹没动态监测的迫切需求,提出了一种基于多源时序SAR数据的洪涝淹没动态监测方法。一方面,不同于现有广泛使用的分类后比较洪水淹没制图方法,本文方法可从多时相异源SAR数据中直接提取洪水淹没信息。另一方面,不同于现有基于局部特征的洪水淹没信息提取方法,本文利用多时相SAR影像的非局部时空信息进行洪水淹没测绘,在抑制噪声影响的同时,可以获取边缘信息可靠的洪水淹没范围测绘结果。为验证本文方法的有效性,以2023年7月底8月初河北省涿州市及周边发生的洪涝灾害为例,利用Sentinel-1,GF-3和涪城一号获取的多源SAR数据,开展了洪涝淹没动态监测与分析实验。结果显示: ①8月1日,涿州市主城区大部分被淹没,随后8月3日,白沟河聚集了大量来自上游的洪水,而到了8月5日,涿州主城区洪水基本退去; ②利用多源SAR数据对洪涝灾害进行时序监测,能够有效地追踪洪涝的增减变化,为灾害救援和灾后评估提供有力支持。

本文的主要贡献为: ①提出了基于后向散射特征的多源SAR数据特征空间对齐方法,有效抑制了传感器参数差异和外界环境造成的多源SAR影像时空异质性; ②基于时空非局部信息,开展了2023年7月底8月初涿州市及周边洪涝灾害的动态监测,分析了7月28日至8月17日涿州洪涝淹没范围的时空变化特征。