1、引　言

随着近几年遥感技术的迅猛发展，遥感数据数量日益增长并广泛应用于环境监测、大气分析和城市规划等领域。其中SAR图像变化检测被应用到军事和民用领域，主要涉及灾害(如洪水、火灾和地震)区域的定位、城市的扩张情况分析以及军事应用中打击效果的评估等方面，研究SAR图像变化检测具有重要意义。

在遥感图像变化检测中，SAR图像利用合成孔径雷达成像，不受包括云、雨等天气气候变化等因素影响，可以提供稳定的数据源，而且星载SAR有着固定的运行周期，能够对地面进行定期的重复观测，因此它比光学遥感更适合对地面进行动态监视。由于SAR成像易受相干斑噪声(Speckle Noise)影响，因此SAR图像去噪(Salmon，2010；Soni 等，2013)问题获得了广泛重视。遥感图像变化检测本质上是一个二分类问题，即将差异图像分为变化类和非变化类。基于像素的方法(Coppin 等，2004；Lu 等，2004)通常利用差异像素灰度特征进行分类，并使用聚类方法获取尽量好的分类结果，CVA(Change Vector Analysis)(Bovolo和Bruzzone，2007)方法则利用差异图像的变化特征向量提取变化部分，Celik(2009)利用PCA降维并使用Kmeans聚类完成变化检测过程，Gong等人(2014)将马尔可夫随机场和模糊聚类相结合根据新的能量函数改进变化检测结果。然而像素特征不具备语义特征，存在语义鸿沟问题，因此利用图像分割分离异质区域从而获得半语义信息成为图像理解和识别的关键技术。常用的图像分割方法包括谱聚类分割算法(贾建华和焦李成，2010；Gou 等，2013)、主动轮廓线分割(Feng 等，2013；Gao 等，2011；Wang 等，2014)、MRF分割(宋晓峰 等，2010；Xiong 等，2012)和SRM(Nock和Nielsen，2004；Lang 等，2014)统计区域合并算法等。近年来，已有部分学者将图像分割方法应用于遥感图像变化检测中，Hichri等人(2013)将基于马尔可夫随机场MRF(Markov Random Filed)和水平集分割方法引入变化检测，并使用SVM(Supporting Vector Machine)进行分类，这是一种监督学习方法需要一定数量的训练数据。万红林等人(2012)将交互式分割和决策融合引入SAR图像变化检测中，但该方法对种子点的选取非常敏感。佃袁勇等人(2016)结合Meanshift算法对高分辨遥感图像进行多尺度分割获得变化检测结果，取得了一定效果。Wang等人(2016)根据提取差异图像中SIFT (Scale-invariant Feature Transform)点并根据特征点利用大金法进行分割，该方法分割效果取决于对SIFT点的有效提取，可以看出基于分割方法的变换检测方法正在逐渐增多，但大部分方法需要训练数据获得稳定的监督分割。SRM算法是一种快速区域生成分割算法，通过建立图像的统计模型利用统计推理方法控制像素及区域的合并，设定同质区域判决规则获得准确的分割结果，SRM方法不依赖于数据的概率分布假设并拥有较好的抗噪能力。本文利用SRM统计分割模型设计了一种新的非监督变化检测方法。首先针对SRM算法单次静态排序易导致的过分割现象，提出了一种基于多特征马氏距离的排序准则和动态排序模式的DSSRM算法；其次建立最小互信息准则构造差异数据集提高合并过程中的约束力，进一步提高检测性能，最后，利用SRM、DSSRM和简化SRM构成级联分割结构完成对差异图像的变化检测。

2、基于DSSRM级联分割的变化检测方法

本节首先介绍SRM分割算法的基本原理和流程，分析算法的不足之处，在此基础之上提出DSSRM算法，其次介绍基于最小互信息准则的差异数据集构造过程，最后提出一种级联分割变化检测框架。

    (2.1) SRM分割方法

经典SRM算法是2004年由Nock和Nielsen提出的一种图像统计模型，算法主要根据图像区域的同质性质在概率空间中寻找最优合并区域，能够感知并快速捕获图像的主要结构，具有较强的抗噪性能，针对不同的内容复杂度可以实现多尺度分割，是一种高效的图像分割算法。

SRM算法主要针对RGB三通道彩色图像，假设图像I含有个像素，每个像素包含(R，G，B)3个颜色通道，每个通道取值范围为{1，2， $ \cdots $ ，g}，g=2^q，q为图像每个通道的量化阶数，每个颜色通道采用Q个独立的随机变量表示，取值范围为[0，g/Q]。图像I是真实场景I^*的一个观测，I^*中每个统计区域代表一个真实对象，因此统计区域应该满足同质特性：

(1) 对于任意统计区域内部的任意颜色通道(R，G，B)，统计像素在该通道内都应该具有同样的均值；

(2) 任意相邻统计区域中至少有一个颜色通道的均值与其他通道不相同。

根据上述统计区域的同质特性，推导出待合并区域对 $(R,\;\; R')$ 的合并准则如式(1)(2)：

$P(R,R') = \left\{ {\begin{array}{*{20}{l}}{\begin{array}{*{20}{c}}\!\!\!\!\!\! 1 \begin{array}{l}\forall a \in (R,G,B),\mathop {\max }\limits_{a \in \left\{ {R,G,B} \right\}} \\[8pt]\left| {\overline {{R'_{a}}} - \overline {{R_a}} } \right| \leqslant \sqrt {{b^2}(R) - {b^2}(R')} \end{array}\end{array}}\\[17pt]\!\!\!\!\!\! {\begin{array}{*{20}{c}}0 \quad\quad{{\text{其他}}}\end{array}}\end{array}} \right.$

(1)

$b(R) = g\sqrt {\frac{1}{{2Q\left| R \right|}}\ln \frac{{\left\| {{R_{\left| R \right|}}} \right\|}}{\delta }} $

(2)

式中， $\overline {R_{{a}}'} $ 和 $\overline {{R_{{a}}}} $ 分别表示统计区域 $R'$ 和 $R$ 在通道a中的平均观测值， ${R_{\left| R \right|}}$ 表示有 $\left| R \right|$ 个像素的区域集合，且有 $\left\| {{R_{\left| R \right|}}} \right\| \leqslant {\left( {n + 1} \right)^{\min \left( {\left| R \right|,g} \right)}}$ ，参数Q代表真实场景I^*统计复杂度，通过调节Q可以获得不同尺度的分割结果，常数 $\delta = \frac{1}{{6{{\left| I \right|}^2}}}$ ， $\left| I \right|$ 表示图像含有的像素个数。

SRM算法的核心技术包括两部分：

(1) 排序算法。在判断是否合并之前首先要确定待合并的区域对 $(R,\;\; R')$ ，SRM算法按照像素四邻域方向梯度幅值大小进行一次静态排序，获得合并优先级列表。排序函数如式(3)：

$f(p,{p'}) = \mathop {\max }\limits_{a \in \left\{ {R,G,B} \right\}} \left| {p_{_a}' - {p_a}} \right|$

(3)

式中，(p，p’)为相邻像素对， $\left| {{p_a}' - {p_a}} \right|$ 为通道a梯度幅值，式(3)表示取R、G、B 3通道中最大值作为最终的排序依据。排列顺序的不同将对最终合并结果产生极大影响，如图1(a)中有标记为1、2、3的3块待合并区域，共有6种排列结果即{[1，2，3][2，1，3]；[1，3，2]，[3，1，2]；[2，3，1]，[3，2，1]}，这6种排序将产生如图1(b)、(c)、(d)所示的3种不同合并结果。

(2) 合并准则。判断满足什么样的条件图像区域对象 $(R,\;\; R')$ 才能合并。SRM算法的合并准则主要考虑图像区域的幅值均值特征 $\left| {\overline {R_a'} - \overline {{R_a}} } \right|$ ，如式(1)(2)所示。

    (2.2) DSSRM方法

统计区域合并SRM算法存在不足之处：

(1)静态排序方式：SRM算法采用一次静态排序方式，排列顺序一次生成，在合并过程中并不更新，由于区域合并是一个动态过程，每次合并都可能生成新的区域，该区域的特征很可能发生改变，从而造成排列顺序的重新分布，显然一次排序无法适应合并区域特征的动态变化过程；

(2)单特征排序准则：SRM算法的排序准则主要基于像素对和区域对特征相关性，仅仅考虑灰度均值特征，初始合并的主要对象为像素，灰度均值能够准确代表像素特征，随着合并区域面积逐步增大，区域形状和纹理逐渐复杂，依赖单特征的灰度均值易造成图像局部过分割现象，如图2被圈中部分所示，Bern数据差异图像变化区域被腐蚀，造成检测漏警增大。

针对所述SRM算法的不足，提出一种多特征判决准则的动态排序算法DSSRM，考虑引入区域均值、区域面积和区域灰度直方图等特征，并采用无量纲的马氏距离构造排序函数，定义最优待合并区域对 $(R,\;\; R')$ 为

$(R,{R'}) = \arg \mathop {\max }\limits_{a \in \left\{ {{D_1},\cdots,{D_n}} \right\}} f({R_a},R_a^{'})$

(4)

式中，f $(R,\;\; R')$ 为排序函数， ${D_1},\cdots,{D_n}$ 表示n个数据通道，自然图像中n=3，表示R、G、B 3个通道，多光谱和高光谱图像中n对应光谱段数据，对于SAR图像而言，只有灰度通道，为了提高合并约束力，在2.3节构造差异数据通道并给出相关构造依据。

排序函数f $(R,\;\; R')$ 的定义如式(5)所示：

$\left\{ \begin{array}{l}f(R,{R'}) = \sqrt {\left[ {{{H}_R} - {{H}_{{R'}}}} \right] \times {{S}^{ - 1}} \times {{\left[ {{{H}_R} - {{H}_{{R'}}}} \right]}^{\rm{T}}}} \\r(R,{R'}) = 1\\\max \left[ {\left| {\displaystyle\frac{{{A_R}}}{{{A_{{R'}}}}}} \right|,\left| {\displaystyle\frac{{{A_{{R'}}}}}{{{A_R}}}} \right|} \right] < \delta \end{array} \right.$

(5)

式中， ${{{H}_R}}$ 和 ${{{H}_{{R'}}}}$ 分别表示区域 $R$ 和 $R'$ 对应的特征向量， ${{H}_{{R'}}} = \left\{ {\mathop {R_{}'}\limits^ - ,A\left( {R_{}'} \right),\left| {{R'}} \right|} \right\}$ ， ${{H}_R} = \left\{ {\mathop R\limits^ - ,A(R),\left| R \right|} \right\}$ ，灰度均值 $\overline R $ 、统计直方图A(R)和区域面积 $\left| R \right|$ 构成特征向量，其中 $\overline R $ 表达了区域平均灰度相似性，A(R)表达了区域灰概率度分布的相似性，而 $\left| R \right|$ 则表达了区域大小相似性，采用无量纲的马氏距离 $\sqrt {\left[ {{{H}_R} - {{H}_{{R'}}}} \right] \times {{S}^{ - 1}} \times {{\left[ {{{H}_R} - {{H}_{{R'}}}} \right]}^{\rm{T}}}} $ 定义待合并区域对特征向量 ${{{H}_R}}$ 和 ${{{H}_{{R'}}}}$ 的相似度，其中S为所有待合并区域特征向量的协方差矩阵，S^–1为其逆矩阵。

$r(R,\;\; R')$ 函数表示区域对 $(R,\;\; R')$ 的相邻关系，1表示相邻，0表示不相邻，根据邻接矩阵元素值确定满足相邻关系的区域对 $(R,\;\; R')$ ，邻接矩阵随着合并进行动态更新。A_R表示区域R中的像素个数， $\max \left[ {\left| {\frac{{{A_R}}}{{{A_{{R'}}}}}} \right|,\left| {\frac{{{A_{{R'}}}}}{{{A_R}}}} \right|} \right] < \delta $ 表示待合并的区域对面积比不能过大以防止图像细节丢失，δ为经验常数。

最后，为了适应区域合并过程中区域特征的动态变化特性，本文采用动态排序方式，每完成一次合并就重新计算新区域的特征向量，更新相邻矩阵，重新计算 $f(R,\;\; R')$ 并更新排序列表，保证下一次合并区域最相似。

综上所述，DSSRM算法流程如下：

(1) 遍历待分割图像I中相邻的区域，确定相邻区域对个数M，生成邻接矩阵B；

(2) 计算每个区域对的相似度 $f(R',\;\; R)$ ；

根据式(4)(5)，确定候选合并区域对 $(R,\;\; R')$ ；

(3) 根据式(1)(2)判断 $(R,\;\; R')$ 是否合并，如果合并，则更新邻域矩阵B及相邻区域对个数M，跳到(2)；否则，遍历所有待合并区域，如果未遍历完，跳到(3)；

结束。

    (2.3) 基于最小互信息准则的的差异数据集构造

根据式(4)计算的排序函数中， ${D_1},\cdots,{D_n}$ 表示n个数据通道，对于自然图像n=3，数据通道分别为R、G、B。SAR图像是由雷达回波构成的合成图像，在空间分辨率上只有雷达回波强度一个通道，而自然图像、高光谱图像有多个数据通道，通道约束能力更强。为了提高SRM算法的变换检测效果，必须构造SAR图像差异数据集，提高SRM算法合并的约束能力。定义输入两时相SAR图像分别为I₁、I₂，对数比值差异图D₁、均值比值差异图D₂和高斯比值差异图D₃定义如式(6)所示，其中 ${\mu _{{I_1}}}$ 和 ${\mu _{{I_2}}}$ 分别是输入时相图像I₁、I₂的邻域灰度均值，g₁和g₂分别是I₁、I₂的邻域高斯滤波值。

$\left\{ \begin{array}{l}{D_1} = \left| {\log \displaystyle\frac{{{I_1}}}{{{I_2}}}} \right|\\{D_2} = 1 - {\rm{min}}\left( {\displaystyle\frac{{\mathop {{\mu _{{I_1}}}}\limits^ - }}{{\mathop {{\mu _{{I_2}}}}\limits^ - }},\displaystyle\frac{{\mathop {{\mu _{{I_2}}}}\limits^ - }}{{\mathop {{\mu _{{I_1}}}}\limits^ - }}} \right)\\[15pt]{D_3} = \log \left[ {\max \left( {\displaystyle\frac{{{g_{_1}}}}{{{g_{_2}}}},\displaystyle\frac{{{g_{_2}}}}{{{g_{_1}}}}} \right)} \right]\end{array} \right.$

(6)

Bern数据集两时相、对数比值差异图像、均值比值差异图像和高斯比值差异图像分别如图3(a)—(e)所示。

建立图像互信息最小准则构造差异数据集，对于输入图像A、B，互信息定义为式(7)：

$\left\{ \begin{array}{l}I(A,B) = H(A) + H(B) - H(A,B)\\H(A) = - \sum\limits_{i = 1}^n {{P_A}({a_i})\log } {P_A}({a_i})\\H(B) = - \sum\limits_{i = 1}^n {{P_B}({b_i})\log } {P_B}({b_i})\\H(A,B) = - \sum\limits_A {\sum\limits_B {{P_{AB}}({a_i},{b_i})} } \log {P_{AB}}({a_i},{b_i})\end{array} \right.$

(7)

式中，H(A)为图像A的熵， ${P_A}({a_i})$ 为图像A的灰度级概率，H(B)为图像B的熵， ${P_B}({b_i})$ 为图像B的灰度级概率，H(A，B)为图像A与B的联合熵， ${P_{AB}}({a_i},{b_i})$ 为图像A和B的联合灰度级概率。通过计算图像A与B的灰度直方图统计H(A)和H(B)，计算图像A与B的联合直方图统计联合熵H(A，B)。

统计常用的SAR图像变化检测Berne、Otawwa以及Yellow River数据的4种差异图像的互信息形成了列表，结果如表1所示，其中L表示对数比差异图，R表示均值比差异图，G表示高斯均值比差异图，符号‘&’连接计算互信息的对象。从表1中可以看出A&B的互信息最小，而B&C最大，即均值比值差异图和对数比值差异图之间的差异性最大，作为SRM最终的差异数据集。

    (2.4) 基于DSSRM算法的级联分割检测框架

在DSSRM分割方法基础之上提出一种混合级联合结构的SAR图像变化检测框架，该框架结构如图4所示。首先利用非局部均值滤波去除SAR相干斑噪声，然后采用差异性较大的对数比和均值比方法构建双通道差异图像，最后利用混合级联结构对双通道差异图像提取变化区域，其中第一级采用原始SRM算法将差异图从像素空间转换为超像素空间，第二级利用DSSRM算法完成差异图像合并，最后使用简化的SRM算法获得最终的变化检测结果。

采用非局部均值滤波滤除SAR图像相干斑噪声，非局部均值通过使用图像中所有像素灰度值的加权平均值来得到待处理像素点的灰度估值，对于SAR图像的相干斑噪声有较好的效果。

已知一幅图像 $v = \left\{ {v\left( x \right),x \in {I}} \right\}$ ，I为图像像素域，在经过非局部均值滤波器处理之后，图像中任意像素点v(x)的灰度估计值为

$NL\left( v \right)\left( i \right) = \sum\limits_{j \in {I}} {\omega \left( {i,j} \right)} v\left( i \right)$

(8)

式中，权值 $\omega \left( {i,j} \right)$ 表示像素v(i)与v(j)之间的相似度，其满足条件 $0 \leqslant \omega \left( {i,j} \right) \leqslant 1$ 且 $\sum\limits_{j \in {I}} {\omega \left( {i,j} \right)} = 1$ 。Nock和Nielsen(2004)提出的SRM方法使用一次排序完成图像合并过程，易产生过分割现象，尤其对SAR差异图过分割现象较为严重，因此采用DSSRM方法进行二次合并过程使分割结果稳定收敛，并采用简化SRM将差异图合并成变化和不变化两部分，这就是级联分割结构。DSSRM使用多特征组合进行分割，但是首先要获得差异图像的超像素初始分割图，简化的SRM算法是为了将DSSRM分割结果继续合并成变化和非变化两个区域并进行二值化，获得最终的检测结果。简化SRM算法放松了原始SRM算法中合并区域必须相邻的条件，对第二次合并结果中每个区域的灰度均值进行升序排序，并按照合并规则进行合并，直至将图像划分为变化区域和非变化两种区域。级联合并的三次合并结果分别如图5(a)—(c)所示。

3、实验结果与分析

为了验证基于DSSRM算法的混合级联框架在SAR图像遥感图像变化检测中的有效性，设计了3种仿真实验加以验证：(1)对比4种变化检测方法在SAR图像变换检测中的性能指标，包括SRM+Kmeans方法、PCA_Kmeans方法、Fuzzy+MRF以及基于双通道差异图及DSSRM算法的级联分割方法简称DSSRM；(2)对比在级联框架下单通道差异图方法、双通道差异图以及三通道差异图条件下的变化检测效果，验证差异图互信息准则的有效性；仿真数据集采用Bern数据集、Ottawa数据集和Yellow River数据集。数据集的分辨率分别为：301×301、290×350、291×306个像素大小。仿真实验的计算环境如下：CPU Intel(R) Core i5-34703.2 GHZ、内存4 G、WINDOWS 7操作系统，仿真软件版本MATLAB2013b。

4种方法在Bern数据集上的检测结果如图6(a)—(d)所示，表2列出了4种方法在3种数据集上的检测性能指标对比，主要包括虚警像素个数F_p、漏警像素个数F_n，正确分类概率PCC以及一致性系数Kappa，并在表中用加粗方式标出了最佳检测指标。从表2可以看出，DSSRM方法相对于其他方法获得了更好的检测结果，尤其在虚警错误数上及总错误数上获得了较好的控制，使得虚警和漏警相对均衡，因此在Kappa系数上也获得了较好的性能指标。这主要得益于级联分割框架中SRM算法为初始分割将像素空间转换为超像素空间，使待合并区域具有部分语义信息，DSSRM利用区域特征使得最相似的区域首先合并，避免过分割现象，双通道数据也对进一步提高了合并的有效性。

为了验证互信息选择准则的正确性，对比在级联分割框架下单通道差异图(以对数比值差异图为例)、双通道差异图以及三通道差异图的变化检测结果分别如图7(a)(b)(c)，各项性能指标如表3所示。双通道差异图获得了最好的检测结果，主要是因为对数比值和均值比值差异性最大，由于均值比值差异图和高斯比值差异图相似性较大，在合并过程降低了合并约束能力，导致了过分割现象，如图7(c)中的红圈部分，因此三通道差异图效果反而变差。仿真实验说明检测效果主要取决于通道差异性。

4、结　论

本文提出了一种基于DSSRM级联分割的SAR图像变化检测方法。DSSRM算法在半语义的图像超像素空间完成收敛合并，简化SRM则获得最终的变化检测结果。通过构造合并区域多特征的马氏距离空间作为排序准则，用动态更新方式替代一次静态排序过程，并利用最小互信息准则构造差异数据集提高合并的约束能力。实验表明，DSSRM的虚警错误比基于Kmeans和MRF的检测方法低，虚警概率和漏警概率更均衡，因此能够获得更高的检测性能。图像合并过程中尺度因子应与图像尺寸相关，如何根据SAR图像特征获得理想的尺度并探索自适应式尺度合并是本文下一步研究的方向。