1 引言

合成孔径雷达(Synthetic Aperture Radar, SAR)具有全天时、全天候的工作能力，不受天气和光照的影响。随着SAR成像系统的发展，大量的高质量SAR图像为遥感应用提供了可靠的数据基础。其中基于多时相SAR数据的变化检测技术已经被广泛地应用于农业调查、环境监测、城市化进程研究以及灾害评估等方面^[1]。近年来，受气候变化和极端天气的影响，洪水已经成为全球频发性自然灾害，不仅严重影响局部生态系统，还会造成巨大的生命财产损失^[2]。尤其是2016年在超强厄尔尼诺现象的影响下，长江中下游地区发生持续性强暴雨，太湖、长江、淮河三大流域水位迅速上涨，导致全国26省、千余县区遭受洪涝灾害。洪水扩张范围的准确获取对规划救援和灾后评估都具有重要的指导作用，因此针对洪水的监测和变化检测研究具有重要意义。

变化检测技术是通过比较同一区域在不同时段的图像来识别该区域的变化情况。变化检测处理流程分为预处理和变化提取两个阶段。针对不同的检测目标对象，目前变化检测的方法可分为像素级的变化检测方法(Pixel-Based Change Detection, PBCD)、对象级的变化检测方法( Object-Based Change Detection, OBCD)和混合变化检测方法(Hybrid Change Detection, HCD)^[3,4]。其中，PBCD以单个像素为基本处理单元生成差异图，并对差异图进行分类。此类方法虽然简单、快速、有效，但受SAR图像固有的相干斑噪声影响较大，检测精度有限^[5]。而OBCD则是分别对变化前后的图像进行分类，然后再比较图像对象之间的差异。研究表明当SAR图像分辨率较高时，OBCD的检测性能要优于PBCD^[6]。然而针对不同的研究对象，分类算法的分割性能不同，另外该方法存在误差累积现象，因此需要根据具体的研究目标改进相应的处理流程。HCD融合了PBCD和OBCD两种方法来进一步降低噪声和杂散点对变化检测的影响^[7]。但是大多数HCD法都属于决策级融合，因此也会同时累积PBCD和OBCD两种方法固有的误差。Lu^[8]等人提出一种更为有效的基于算法级融合的HCD方法来检测洪水变化，该方法避免了未变化区域的误差累积效应，同时利用空间文本信息来提高检测精度。由于该方法在对象级变化检测阶段采用区域增长的方法提取水域区域，因此算法最终的提取精度取决于种子点的分布，当水域分布不连续时，容易漏掉局部区域，降低提取精度。

基于多时相SAR数据的变化检测所获取的地物变化信息可以被用来对洪水进行有效监测。Addabbo^[2]等人采用贝叶斯网络融合SAR强度信息、InSAR相干信息和其他参考图像，来检测布拉达诺河洪水区域。Pulvirenti^[9]等人通过采用形态学方法对COSMO-SkyMed图像进行自动分割来获取洪水演化图。Avendano^[10]等人采用无监督的分类方法对SAR图像进行分割，并融合平均比和对数比信息来获取最终的洪水监测结果。

为了提高多时相SAR数据变化检测的精度，本文提出一种改进的HCD方法。该方法在对象级变化检测阶段采用基于初始聚类中心选取的改进模糊C均值(Fuzzy Clustering Method, FCM)聚类算法来避免局部水域漏检的情况。另外本文采用Sentinel-1A卫星干涉宽测绘带模式多时相多极化SAR图像来监测洪水的变化情况。因此在像素级处理阶段，采用主成分分析法(Principal Component Analysis, PCA)融合分别由HV极化和VV极化SAR图像得到的差异图，进一步降低斑点噪声的影响。实验结果表明，本文的方法具有很强的鲁棒性，并且针对不同实验区域，本文方法的平均总误差率为3%，平均检测率为92%，而文献[8]中的HCD方法的平均总误差率为4%，平均检测率为87%。

2 算法描述

本文提出的算法的目标是融合PBCD和OBCD两种方法来提高变化检测的性能，算法的核心为通过PBCD处理结果来确定FCM的初始聚类中心和最近距离分类(Nearest Neighbor Clustering, NNC)的聚类中心。另外在进行HCD处理之前，对SAR图像做直方图均衡来保证图像具有相似强度的分布，避免图像强度差异带来的检测误差。算法流程如图1所示。

2.1 PBCD检测

PBCD检测主要分为两步：差异图生成、差异图分割。若在PBCD检测阶段，图像对象相关的参数估计得较准确，则在OBCD检测阶段更容易提取出精度较高的图像对象^[8]。因此差异图分割的目的就是获取准确的图像对象。

2.1.1 差异图生成 由于SAR相干斑噪声属于乘性噪声，针对多时相SAR图像的残差计算，比值运算更适合SAR图像的统计特征。因此本文采用局部平均似然比变化检测算子进行差异图像的生成，并利用小矩形窗来代替局部同质区域，则平均似然比的简化计算公式可表示为^[11]：

$\eta \left( {i,j} \right) = \frac{{\sum\limits_{\left( {m,n} \right) \in {V_{i,j}}} {R\left( {m,n} \right)} }}{{\sum\limits_{\left( {m,n} \right) \in {V_{i,j}}} {F\left( {m,n} \right)} }}{\rm{ + }}\frac{{\sum\limits_{\left( {m,n} \right) \in {V_{i,j}}} {F\left( {m,n} \right)} }}{{\sum\limits_{\left( {m,n} \right) \in {V_{i,j}}} {R\left( {m,n} \right)} }}$

(1)

其中，R为参考强度图，F为洪水强度图，V表示像素(i, j)在w×w窗内的局部邻域像素集(w通常取值为3)。随后将似然比图像的取值范围调整为[0, 255]。

本文采用的Sentinel-1A SAR数据包含HV极化和VV极化两种图像。在图像处理技术中，图像融合通过联合不同来源的信息产生相对于单幅图像来说质量更高的图像^[12]。因此在本文提出的算法中，PCA被用来融合分别由HV极化时相图像对和VV极化时相图像对生成的差异图。PCA融合的具体实施步骤如下所示^[13]：

(1) 由P组差异图中的所有数值构成一个2维矩阵 ${X} = {\left[ {{{x}_i} \cdots {{x}_P}} \right]^{\rm{T}}}$ ，其中 x_i 表示将第i 幅图像矩阵转化为行向量。则矩阵 X 的行表示以单幅图像作为变量，矩阵 X 的列代表这些变量的观测值。

(2) 计算 X 的协方差矩阵：

$K = \frac{1}{{N - 1}}{X}{{X}^{\rm{T}}}$

(2)

其中N表示单幅图像中所有像素的个数，K表示对应的观测对之间的相关性。

(3) 独立测量值之间的协方差为0，才能够实现通过冗余观测值消除噪声。因此对K进行特征值分解：

$K = {ED}{{E}^{\rm{T}}}$

(3)

其中 E 表示特征向量， D 表示由特征值构成的对角矩阵。

(4) 将特征向量按照对应的特征值从大到小重新排列，则差异图对应的主成分可表示为：

${\tilde {{X}}}= \left( {{{{E}}^{\rm{T}}}} \right){{X}}$

(4)

(5) 其中 ${\tilde {{X}}}$ 的大小为P×N, ${\tilde {{X}}}$ 中每一行都代表差异图数据集的主成分，图像中的信息主要包含在第1行(第一主成分)中，其他行则以噪声为主^[13]，因此本文只保留第一主成分，并将其重新转化为2维矩阵，作为最终的差异图。

2.1.2 图像对象提取 在洪水监测应用中，图像对象指水域对象，初始分割得到的变化区域直接对应洪水SAR图像中的核心水域，另外该区域也对应着参考SAR图像中的植被或岸滩区域。初始分割方法包括阈值分割、主动轮廓模型、马尔科夫随机场等技术^[14–16]。阈值法简单、快速，因此应用比较广泛。本文利用差异图的相邻直方图比值方法进行阈值选取，从而获取初始分割。差异矩阵直方图的峰值代表无变化区域，浮动区域代表变化区域，因此这两个区域之间的过渡点可以设为检测门限^[11]。过渡点的查找可以通过计算差异直方图h _{$\eta $} 相邻灰度级之间的比值来实现，计算公式可表示为：

$r\left( k \right) = {h_\eta }\left( k \right)/{h_\eta }\left( {k{\rm{ + }}1} \right)$

(5)

其中， $k \in \left[ {{k_0},255} \right)$ , k₀为h_h最大值对应的灰度级。一旦r(k)小于1，初始分割门限t_Init可设置为：

${t_{{\rm{Init}}}} = \left\{ {\left. k \right|r\left( k \right) \le 1} \right\}$

(6)

采用t_Init对差异图进行初始分割，那么核心水域的掩膜可定义为：

${M_{{\rm{Init}}}} = \left\{ {\begin{array}{*{20}{c}}{\begin{array}{*{20}{c}}\!\!\!\!\! {1,} & {\eta \left( {i,j} \right) \ge {t_{{\rm{Init}}}}}\end{array}}\\{\begin{array}{*{20}{c}}\!\!\!\!\! {0,} & {\eta \left( {i,j} \right) < {t_{{\rm{Init}}}}}\end{array}}\end{array}} \right.$

(7)

2.2 OBCD检测

OBCD检测的关键步骤就是分别识别多时相图像中的图像对象，然后通过比较图像对象之间的差异来实现变化检测，因此图像分割的质量尤为重要。在本文的应用中，水域为待研究的图像对象。针对水域对象提取，相关研究方法包括阈值分割、聚类、主动轮廓模型等^[17–19]。鉴于SAR图像中水域表现为低灰度同质区域，文献[8]选取核心水域中的部分像素作为种子点进行区域增长，从而得到洪水区域。结合文献[8,17–19]可知，相同地物在SAR图像中都具有一定的同质性，因此可以选取不同地物的部分像素作为种子点进行分类。本文提出的OBCD检测算法分为3步：估计初始聚类中心，FCM聚类，最近距离分类(NNC)。

2.2.1 估计初始聚类中心 本文所提出的算法关注点在于水域和非水域之间的区分，无需较多分类，因此文中将整个图像分为3类：水域、背景区域、过渡区域。各个类别的聚类中心可通过种子点的选取来确定。

首先采用掩膜M_Init提取洪水SAR图像中的核心水域，并选取该区域强度直方图峰值对应的强度值作为门限，提取水域种子点，则水域种子点门限 $\tau $ _w 为：

${\tau _w} = \left\{ {\left. l \right|{h_w}\left( l \right) = \max \left( {{h_w}\left( l \right)} \right)} \right\}$

(8)

其中，h_w 为核心水域强度直方图。水域的初始聚类中心即为水域种子点的强度均值。

接着采用掩膜M_Init提取参考SAR图像中的对应区域，同样选取该区域强度直方图峰值对应的强度值作为门限，提取背景区域种子点，则背景区域种子点门限 $\tau $ _b为：

${\tau _{\rm{b}}} = \left\{ {\left. l \right|{h_{\rm{b}}}\left( l \right) = \max \left( {{h_{\rm{b}}}\left( l \right)} \right)} \right\}$

(9)

其中，h_b为背景区域强度直方图。背景区域聚类中心即为背景区域种子点的强度均值。

对于FCM聚类，将水域种子点的强度均值作为第一个聚类中心，计算每个像素点与最近聚类中心间的最短距离dis(x)，根据概率分布 ${\rm{dis}}{\left( x \right)^2}/\sum\nolimits_{x \in {X}} {{\rm{dis}}{{\left( x \right)}^2}} $ 来选取新的像素点作为新的聚类中心，直至选取了3个聚类中心点，其中x为像素坐标。

NNC分类需要对过渡区域进行二次分割，而植被密集区域灰度级较低，因此将上述水域的初始聚类中心和背景区域聚类中心的均值作为背景二次分类中心，同时仍将水域的初始聚类中心作为水体二次分类中心。

2.2.2 FCM聚类 FCM是一种无监督聚类，该算法对初始值敏感，迭代容易陷入局部极值, 难以取得全局最优。文献[20]通过设定聚类中心之间的最小距离来保证各类之间的分离性，避免随机求取初始聚类中心时容易使算法收敛到局部极小的情况。本文提出的算法根据概率分布来确定初始聚类中心，同样保证了各类之间的分离性。FCM聚类通过迭代求解目标函数最小化的过程来确定每一个样本隶属于某一类的隶属度。算法步骤如下所示：

(1) 给定类别数c，模糊指数m，容许误差$\xi $。其中 $m \in \left[ {{\rm{1}},{\rm{ + }}\infty } \right)$ ，当m=1时，该算法退化为k均值聚类，一般情况下，m取值为2。

(2) 按前述方法设置初始聚类中心，并设置循环次数k=1。

(3) 计算样本的隶属度，即SAR图像中所有像素的隶属度，其计算公式为：

$\begin{array}{l}{\mu _{ij}}\left( k \right) = \frac{1}{{{{\sum\limits_{k = 1}^c {\left[ {d_{ij}^2\left( {{x_j},{\omega _i}} \right)/d_{kj}^2\left( {{x_j},{\omega _k}} \right)} \right]} }^{1/\left( {m - 1} \right)}}}},\\\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;i = 1, \cdots ,c; \ j = 1, \cdots ,N\end{array}$

(10)

其中，i代表类别编号，j代表样本编号，$\omega $为聚类中心，d为欧式距离，N为图像所有像素个数，则μ_ij 为像素样本x_j 对${\omega _i}$的隶属度。

(4) 修正所有的聚类中心${\omega _i}$(k+1)，修正公式为：

${\omega _i} = \frac{{\sum\limits_{j = 1}^N {{{\left( {{\mu _{ij}}} \right)}^m}{x_j}} }}{{\sum\limits_{j = 1}^N {{{\left( {{\mu _{ij}}} \right)}^m}} }}$

(11)

(5) 计算误差 $e = {\sum\limits_{i = 1}^c {\left\| {{\omega _i}\left( {k + 1} \right) - {\omega _i}\left( k \right)} \right\|} ^2}$ ，如果e ≤ $\xi $，则停止迭代，否则重复(3), (4), (5)步。

(6) 迭代结束后，如果 $d_{ij}^2\left( {{x_j},{\omega _i}} \right) < d_{kj}^2\left( {{x_j},{\omega _k}} \right),$ $k = 1, \cdots ,c,k \ne i$ ，那么可以将x_j 归为第i类。

2.2.3 NNC分类 在本文的OBCD检测阶段，先采用FCM将图像分为3类，即水体目标、背景、过渡区域。受灰度级较低的植被区域或相干斑噪声的影响，水体与背景有重叠部分，这部分像素主要包含在过渡区域中，因此本文通过计算最近距离对过渡区域进行二次划分。计算公式为：

${\mathop{\rm label}\nolimits} \left( {{I^{l \in {{{Ω}} _{{{ω}_i}}}}}} \right) \!\!=\!\! \left\{ \begin{array}{l}\!\!\omega _w^2,{\left\| {{x_l} - \omega _w^1} \right\|^2} < {\left\| {{x_l} - \omega _b^1} \right\|^2}\\\!\!\omega _b^2,{\left\| {{x_l} - \omega _w^1} \right\|^2} \ge {\left\| {{x_l} - \omega _b^1} \right\|^2}\end{array} \right.$

(12)

其中，I为坐标空间， ${{{Ω}} _{\omega_i}}$ 为过渡区域像素空间，x_l 为像素样本， $\omega _w^1$ 为水体二次分类中心， $\omega _b^1$ 为背景二次分类中心。将标记为 $\omega _b^2$ 的像素与FCM划分的背景部分合并，同时将标记为 $\omega _w^2$ 的像素与FCM划分的水体部分合并，从而得到最终的水体分布。

3 实验结果及分析 3.1 实验数据说明

本文采用Sentinel-1A获取的多时相多极化SAR图像来进行算法验证，实验数据具体参数如表1所示。2016年入汛以来，我国已有多条江河水位超过1998年同期，6月下旬，我国南方降雨集中、强度加大，且强降雨落区与前期降雨叠加，与常年同期相比，6月以来，长江、淮河及太湖流域降雨量偏多3成到1倍。气象卫星监测显示，6月29日淮河蒙洼至颍上段部分河道水体面积较6月16日有所增大。经估算，淮河蒙洼段约25 km河道水体宽度增加0.5–1.5 km；淮河颍上段约10 km河道水体宽度增加约1 km。另一处具有代表意义的洪水场景为鄱阳湖水域，受强降雨和长江江水倒灌双重影响，鄱阳湖水面明显增大。气象卫星对鄱阳湖水面的遥感监测显示，7月7日鄱阳湖主体及附近水域面积为3966 km²，较今年6月5日监测的水域面积(3559 km²)增大407 km²，比历史同期(平均为3287 km²)偏大679 km²。因此本文选取相近时间段的SAR图像做洪水变化检测研究，如图2和图3所示：图像场景1为淮河蒙洼至颍上段SAR图像，数据获取日期分别为2016年6月11日和2016年7月5日；图像场景2为鄱阳湖水域SAR图像，数据获取日期分别为2016年6月11日和2016年7月5日。其中图2(c)和图3(c)中蓝色区域表示未变化水体，红色区域为洪水扩张区域。为了更好地进行算法的细节展示，文中分别选取典型区域A, B进行算法各个步骤的分析。A区域的图像尺寸为757×1390, B区域的图像尺寸为2048×1536，如图4所示。

3.2 变化检测步骤分析

通常情况下，使用局部平均似然比来获取多时相SAR数据的变化差异图，这种方法可以在一定程度上降低相干斑噪声的影响。然而在HV极化方式下，有些植被区域的后向散射强度变化较明显，因此会对洪水区域的变化产生干扰，如图5(a)所示。与HV极化图像相比，VV极化所获取的似然比差异图中，洪水区域和背景区域的区分比较清晰，如图5(b)所示。针对多源数据，PCA被用来融合不同极化方式下获取的似然比图。相对于单一数据，该方法可以得到更可靠的数值，有利于后续同质区域的划分。由图5中的数值可以看出，在经过PCA融合后得到的图像中，变化区域和不变区域之间的差值变大。

将最终得到的差异图转换为灰度图，并计算得到其归一化直方图，如图6(a)所示。根据公式(5)计算相邻灰度值之间的比值，得到图6(b)所示的比值曲线，选取首个小于1的点对应的灰度级作为PBCD检测的阈值。则PBCD提取的变化区域如图6(c)所示。大部分变化区域都能提取出来，但是细节上有缺失。

3.3 变化检测定量分析

为了进一步验证算法的有效性，将本文算法与其他方法的检测结果进行对比。这些方法包括Lu等人提出的HCD法^[8]，Li等人提出的基于两级聚类的无监督变化检测法^[21]，试差法(Try and error)^[11]以及基于FCM的对象级变化检测方法。其中文献[21]的方法和试差法都是直接对变化差异图进行处理，而基于FCM的对象级变化检测方法则先基于Gabor特征对变化前后的SAR图像进行分类，并对提取出水体分布做差值计算得到最终的变化结果。本文选取检测率、虚警率、总误差率和Kappa系数4个参量进行性能评估，其中检测率为正确检测出的变化像素占总的变化像素的百分比；虚警率为被误分为变化区域的背景像素占背景像素总数的百分比；总误差率为误检像素与漏检像素之和占总像素个数的百分比。

为了便于后续分析，图7(a)和图8(a)给出了手动提取的洪水变化分布作为真实变化区域。由于变化差异图的计算引入局部平均，因此图像中部分细小变化丢失，因此直接对变化差异图进行处理，不能得到完整的变化信息，如图7(d)、图7(e)所示。而文献[8]中的方法在图像中存在局部封闭区域的情况下，检测结果取决于种子点分布的均匀性，如图7(c)所示，部分封闭区域因为没有种子点而被遗漏，影响了最终的检测结果。从图7(f)和图8(f)可以看出，本文方法的检测效果要优于直接采用FCM进行聚类的检测效果。B区域场景较复杂，图8所示的直观对比效果有限，因此本文通过表2的变化检测精度量化对比进一步分析算法的有效性，针对B区域这样的复杂场景，相比于直接由变化差异图提取变化区域的方法，本文所提出的方法检测率大幅提高，虚警率略高，总体误差率较低，Kappa系数较高。

4 结论

本文在算法级混合变化检测的基础上，提出综合考虑像素级变化检测与聚类中心之间关系的变化区域提取方法。通过初始变化区域提取结果分别估计FCM和NNC这两种串联分类方法的聚类中心，并分别对洪水发生前后的两组SAR图像进行分类，将最终估计得到的水体区域做差值计算。通过对Sentinel-1A多时相双极化SAR图像变化检测进行验证，该方法的检测率较高，总体误差率相对较低，可以在雷达遥感图像的洪水监测中进行应用。