洪灾是最常见和最具破坏性的自然灾害之一，它严重损害了人类的财产和生命安全^［1］。对洪灾发生前后进行水体变化检测，确定洪水淹没范围，对灾情的评估和灾后的管理规划具有一定的指导作用^［2-4］。合成孔径雷达（Synthetic Aperture Radar， SAR）可以全天时、全天候地获取遥感数据，不受光照和天气的限制，弥补了光学和红外遥感的不足^［5-6］。随着遥感技术的发展，工作在不同波段的海量SAR影像为基于SAR图像的洪水监测提供了数据基础^［7］。

目前，基于SAR影像的水体变化检测根据是否需要先验知识可以分为监督和非监督两种检测方法^［8-10］。前者需要用已知的分类器进行样本训练，再利用其中已被训练的样本检测图像的变化区域^［11］。而后者不需要先验信息，可以直接对两时相图像间的变化进行检测，地面参考数据的缺失使得非监督变化检测方法在实际应用中更为广泛^［12-13］。使用非监督变化检测方法获取水体变化区域时，有两大关键步骤：差异图的获取、差异图的处理与分析^［14］。其中差异图的获取是非监督变化检测研究中很重要的一步，直接决定了检测精度^［15］。进行差异图的获取时，较为常用的方法有差值法和比值法^［16-18］。也有研究人员将主成分分析、对数比值法、邻域比值法等应用于差异图的构建^［19-20］，但类型单一，往往检测效果不够理想。为了降低变化检测结果的错误率，研究人员采用对数比值和均值比值差异图融合的方法^［21-22］，但是该算法忽略了差值差异图的变化信息。

基于以上分析，提出组合差异图和FCM 聚类的非监督洪水变化检测算法，对2021年7月河南省新乡市洪涝灾害进行动态监测。首先，使用差值法和均值比值法得到两幅强度差异图像；其次，为了充分利用两幅差异图像的变化信息和考虑邻域信息获取完整水灾信息，提出局部自适应加权和双边滤波融合规则计算组合差异图，得到准确反映水体变化信息的组合差异图；最后，利用FCM聚类算法将差异图分为变化区域和非变化区域，其中变化区域对应着洪水淹没区域。实验结果表明，与单一差异图和等权组合差异图检测结果相比，该算法大大地降低了检测结果的总错误数，提高了洪灾检测精度。

2021年7月中下旬，河南省出现了历史罕见的极端强降雨，新乡等市遭受特大洪涝灾害，造成重大人员伤亡和经济损失。新乡市位于河南省的北部（图1），总面积约为8 269 km²，是河南省北部最重要的交通枢纽之一。新乡市下辖4个市辖区（凤泉区、牧野区、卫滨区和红旗区）、2个县级市（长垣市、卫辉市）、6个县（辉县、获嘉县、原阳县、封丘县、延津县和新乡县）。新乡市属暖温带大陆性季风气候，受副热带高压影响十分显著，四季变化明显。新乡市南面靠近黄河，北部依靠太行山，自北向南海拔呈现从高到低的变化规律，地貌类型多样，北部为太行山石灰岩半干旱区，属于山地丘陵地貌。中部和南部为古黄河冲积平原，地形平坦，海拔坡度较低。新乡市南邻黄河，河流、湖泊众多，水资源丰富。其中，主要的河流包括黄河、金堤河、淇河、天然文岩河等，还有众多小的支流纵横交错于整个区域。

Sentinel-1是由A、B两颗同轨的极轨卫星组成的星座，是欧空局“哥白尼计划”的地球观测卫星系列之一，双星重访周期最短为6 d。数据影像有4种成像模式，分别为IW （Interferometric Wide）、EW （Extra Wideswath）、SM（Stripmap）和WV（Wave）。因其影像数据可免费下载，且影像质量和分辨率较高，Sentinel-1卫星数据已广泛应用于洪水监测和灾后损失评估。研究选择我国河南省新乡市洪涝灾害前后两时相（2021年7月15日和7月27日）的Sentinel-1 A IW模式SAR影像进行洪水监测处理，数据具体参数见表1，Sentinel-1 A SAR影像包含VV和VH两种极化方式，一般情况下，地物回波强度在交叉极化VH模式下比同极化VV低很多，为了更好增大水体和其他地物的差异，因此研究选择VV极化模式。两景VV极化模式数据经预处理（轨道修正、配准、滤波、地理编码、裁剪）后得到两时相的VV极化强度图像如图2所示，总像元数为18 808 547。结合研究区地理位置图（图1），对比前后时相强度图可以看出，受灾区域呈东西向分布，主要在卫河流域。

水体作为一种特殊的地物类型，其与电磁波的相互作用以镜面反射为主，所对应的SAR强度较弱，这也是对水体进行识别的主要依据之一。在本研究中，因两景影像获取时间间隔较短，同一位置所对应的两景影像后向散射强度若发生改变，则认为是因降水引起的，前后时相强度改变可以通过构建差异图来进行识别。进行差异图的获取时，较为常用的方法有差值法和比值法。

基于SAR强度信息进行洪水监测，可以根据不同时相图像间的差异来获取洪水淹没区域。图像间的差异通过构造差异图来获取，因SAR影像受乘性噪声的影响，因此基于SAR强度特性的差异图构造可以通过图像差值法和均值比值法进行构造。构造公式分别为：

其中：I2、I1分别表示变化前后的两景SAR影像强度值；D1表示差值差异图；μ1(i,j)和μ2(i,j)分别表示I1和I2在点（i，j）处对应像素的局部均值；D2表示均值比值差异图。

差值法可保留弱变化区域的细节信息，但受噪声影响较大；比值法可以在一定程度上抑制噪声，却会夸大不同时相影像中灰度值均较低区域的变化。单一类型的差异图存在检测精度低以及适用范围窄等问题，综合考虑上述因素，研究综合了差值图和比值图的信息优势，将其运用于组合差异图的构造中，在保留图像细节的同时，减少噪声的影响。

在对多个差异图进行组合时，各差异图的权值选择尤为重要。对各差异图的组合一般采取简单平均或者加权平均的方法，容易造成图像细节信息与能量的丢失，从而降低组合差异图的对比度。水体区域一般是连续且非孤立的点，像素之间具有很强的关联性，因此考虑像素邻域信息有利于提高水体变化区域的完整性。

在此基础上，首先采用局部能量加权算法最大限度保留各差异图的能量，具体步骤为：

（1）根据局部能量函数分别计算差异图D1、差异图D2每个像素点及邻域像素灰度值之和，将其定义为该中心像素点的能量值，能量值较大时，所对应的权值也相对较大。能量函数数学表达式为：

其中：Dμ分别表示差异图D1和 D2，Eμ表示差异图D1和 D2的能量值；w表示邻域窗口大小，一般取3。

（2）计算组合系数：

最终组合差异图表达式为：

为了最大限度获取源图像的细节信息，采用双边滤波规则。它通过空间距离和高斯核的组合过滤空间和距离域中的图像，是一种非线性边缘保留滤波器，其公式可表示为：

其中：DCo'表示原组合差异图像素点经双边滤波后新的灰度值；R(i,j)表示以点（i，j）为中心的灰度取值范围；DCo(k,l)表示原组合差异图像某一点的灰度值；q(i,j,k,l)表示使用叠加的高斯函数计算出的新灰度值。

对组合差异图DCo'提取变化信息时，有多种方法可对其进行水体变化检测，研究采用经典的FCM算法来完成洪水/非洪水分类，FCM算法可以很好的考虑图像的统计特性，并计算每个样本对所有类的隶属度^［13］。总体来说，本文方法有3个步骤：①均值差异图D1和差值差异图D2的生成；②利用局部能量加权和双边滤波算法组合规则得到最终组合差异图DCo'；③对DCo'进行聚类分析获得检测结果图。

水体表面相对光滑，以镜面反射为主，回波信号较弱，在前、后时相影像的相干图上表现为低相干，对应的相干系数较小。根据水体的相干特性，在相位解缠时，设置适当的相位解缠阈值和解缠等级，可以使水体在干涉图中得到凸显。一般来说，我们可以将洪水前后时相的SAR影像进行复共轭相乘得到干涉图，干涉相位可表示为：

其中：Real（I）和Imag（I）代表干涉图的实部和虚部。以Phase、Real（I）和Imag（I）进行RGB假彩色合成，可以获取水域分布，灾后相对灾前的扩张水体即为InSAR处理获取的洪水淹没区。已有研究证明在分别中高分辨率SAR影像进行多时相地表变化检测时，相位信息的检测精度优于强度信息^［14］。因此，将通过彩色合成提取的洪水淹没区作为参考真值图，检验本文方法的优劣。研究的整体技术流程图如图3所示。

将2021年7月15日和7月27日所采集影像的强度图按照公式（1）和公式（2）构造差值差异图和均值差异图，结果如图4所示，图中像元值越接近于1，表示地表发生变化的可能性越大，即可表示洪水区域，相反，越接近于0，表示非洪水区域。图中大部分像元呈蓝色调，差异图的值较小，即多数区域为非洪水区域。在两幅差异图中洪水淹没区域表现较为明显，分布集中新乡市中部地区，呈东西走向，沿卫河流域，在新乡市南部黄河流域，也有部分区域在差值差异图上像元值较大，接近于1，表现为受灾区域。图4（a）显示的局部放大图可以看出，差值差异图受噪声影响较为严重，图像细节信息损失严重。统计差值差异图的直方图如图4（a）中所示，差值差异图的像元值主要集中在0.2附近，均值为0.244 1，方差为0.020 4。而均值差异图的像元值集中在0.4附近， 均值为0.377 7，方差为0.041 3。对比图4中的局部放大图可以看出，受灾区域均可与背景区分开。相较于差值差异图，均值差异图的优点是其变化信息内部均一性更好，缺点是变化信息的强度被抑制，且洪水/非洪水在均值差异图上对比度不高，对受灾区域的提取增加难度。

为真实反映受灾信息，提出一种局部能量加权和双边滤波组合规则进行差异图的融合，将差值差异图和均值差异图按照3.2节算法描述进行组合，组合效果如图5（a）所示。根据组合差异图的统计直方图可以看出，图中绝大多数像元值集中在0值附近，对应该研究区范围内大范围的未受灾区域，而在卫河和黄河流域范围内，存在像元值较大的区域，对应受灾区。研究提出的组合差异图方差为0.015 2，相较于单一差异图来说，方差最小，抑制了噪声的影响，拉大了受灾/非受灾区域图像对比度，受灾区域得到有效增强。

为与本文方法对比，同时生成了等权组合差异图，即将差值差异图和均值差异图的组合权重因子都设为0.5，结果如图5（b）所示。等权方式得到的组合差异图方差为0.029 8，大于图5（a）的方差，且噪声明显，受灾/非受灾区域图像对比度低于图5（a）所示的组合差异图。

利用研究提出的局部能量加权和双边滤波规则得到的组合差异图经FCM聚类后，洪水检测结果如图7（a）所示，图中红色区域表示受灾区域。同时，为比较本文算法的优势，对预处理后的Sentinel-1 SAR影像数据，分别采用单一差异图法、等权组合差异图法来进行洪水信息提取，提取结果如图7（a）~图7（d）所示，图7（e）表示的是由相干法得到的参考真值图，其洪水像元数为360 544。与图7（e）相比较，4种方法都可以有效检测出受灾区域，主要分布在卫河流域和黄河流域。其中，单一差值差异图提取的受灾像元数为186 640（图7（b）），单一均值差异图提取的受灾像元数为246 151个（图7（c））。本文方法提取的洪水淹没像元数为284 671，提取的水体边缘连续性和区域完整性相对较好，受噪声影响较小；等权组合差异图提取的洪水淹没像元数为238 277个。为了定量对比4种水体提取方法的优劣，将4种方法得到的水体信息与相干法结果进行比较，采用虚检率FP（probability of false）、总错误率OE（total error）、正确率PCC（Probability of Correct Classification）和Kappa系数（K）等多种评价指标来定量评价洪灾信息提取的精度，计算公式如下所示：

其中：P12表示被算法判断为洪水的非洪水像元数；P21表示被算法判断为非洪水的洪水像元数；P11表示算法提取到的洪水像元数；P22表示算法提取到的非洪水像元数；P1表示参考真值图中洪水像元数；P2表示参考真值图中非洪水像元数；P0表示总分类精度。PCC越高表示变化检测的性能越好，K通常用来衡量变化结果图与参考图像之间的相似性，评价结果如表2所示。本文算法的OE为0.35%，是4种方法中最低的，大大降低了检测结果的虚警率。而单一差值差异图检测的OE最大，为1.19%。以组合差异图的方式监测洪水，OE是要明显低于单一差异图，PCC是要高于单一差异图；以等权方式组合差异图，其Kappa系数为0.785 5，低于本文方法的0.899 6。而且，本文算法检测的FP显著低于其他3种方法，为0.84%。总体而言，与其他3种方法相比，以本文方法检测洪水虚警率最低，检测的像元总错误数最少，Kappa系数最大，表明本文方法不仅在很大程度上抑制了噪声对检测结果的影响，而且检测精度要明显优于其他方法。最终经计算，基于本文方法监测到的新乡市洪水面积为77.05 km²，新乡市政府在2021年7月22日公布的成灾面积110 km²，本文计算洪水面积是利用7月12日和7月27日的影像，考虑到洪水部分退去等因素影响，本文检测的洪水面积基本与政府公布数据吻合，进一步验证了本文算法的可靠性。

遥感技术是洪灾监测的主要手段，而SAR作为一种主动式遥感，其全天时全天候的优势在灾害评估中发挥着巨大作用。借助SAR的独特优势对新乡市的洪灾进行监测，并提出了组合差异图和双边滤波组合的算法，该算法首先根据局部能量加权函数分别计算差值差异图和均值差异图的每个像素点及邻域像素灰度值之和，确定原像素的组合权值，实现局部自适应权值计算，然后，根据两幅差异图每个像素的权值系数获得组合差异图，最后，根据双边滤波算法对组合差异图进行滤波，得到最终的差异图。实验结果证实该算法的检测结果精度优于单一差异图和等权组合差异图，另外还有其他几个方面在后续研究中可进一步提高，主要包括：①研究受数据获取时间的影响，其结果与政府公布的实际受灾面积稍有差异，洪灾检测往往要求前后时相影像间隔时间较短，但目前的研究受限于传感器的重访周期；②对于Sentinel-1 SAR影像，我们考虑到地面物体在VH回波强度远低于VV极化，本研究仅使用极化方式为VV的强度影像，但是对于SAR的相位信息未得到充分利用，因此，在后续研究中，将继续分析不同极化条件下地物相干性的差异，以期获得更好的检测结果。③对于聚类方法，本文采用传统FCM方法，后续工作中将进一步考虑聚类算法的优化。④另外，还需进一步探索所提出的方法是否可以与其他方法相结合，例如基于深度学习的方法。

（1）利用单一差值差异图、单一均值差异图、等权组合差异图和加权组合差异图，四种方法的检测结果正确率均超过98%，均可以对新乡市2021年洪灾的灾情进行有效评估，揭示洪水淹没范围主要集中在卫河流域。

（2）单一差值图检测结果的虚警率为55.14%，远高于单一均值差异图和组合差异图，且检测的错误像元数最高（OE为1.19%），Kappa系数最小（0.585 8），检测结果精度最差；

（3）无论是等权还是局部能量加权，组合差异图检测结果的虚警率都要低于单一差异图结果，本文方法虚警率最低，为0.84%；

（4）局部能量加权和双边滤波组合规则大大提高了洪水检测的Kappa系数（0.899 6），使检测精度得到提升，用本文方法检测到新乡市2021年7月份洪水受灾面积约为77.05 km²。