1 引言

相干变化检测(Coherent Change Detection, CCD)是一种可以检测场景中发生的微小变化的技术，它利用相位信息形成的相干性观测度量可以探测出亚波长量级的变化，并以相干性变化差异的形式将其表现出来^[1]。这是利用相干信息相对于非相干信息用于变化检测的优势之处。通过CCD处理形成的相干变化差异图具体表现为：场景中发生变化的区域呈现低相干特性(相干系数趋于0)，未发生变化的区域呈现高相干特性(相干系数趋于1)。然而，场景中并非所有呈现低相干特性的区域均为我们所关注的变化区域(称为目标变化区域)，也可能是低相干干扰造成的虚警。这些低相干干扰区域包括回波散射强度弱形成的低信噪比区域，极易受风吹雨淋造成体散射失相干的植被区域等。形成这些低相干干扰的根本原因是：CCD方法采用相位信息进行的相干性度量过于敏感，许多不是我们关注的区域也因发生足以用CCD方法检测出的相位变化而被SAR传感器敏锐地探测出来。因此，如何从低相干区域中排除引起虚警的干扰，是改善CCD性能的关键技术，最直观的解决方法为对低相干区域分类。传统的CCD方法只能区分高、低相干区域，不具备进一步区分能力。因此需要增加新的观测量实现进一步的分类。

目前，国内外的学者针对这一问题开展了许多研究工作。文献[2,3]从时间维度扩展，分别建立以时间为参量的模型对场景中的变化区域分类，提取出仅在特定时间发生变化的低相干区域；文献[4]考虑全极化情况，利用极化相干矩阵建立检验统计量，对不同极化散射类型的变化场景分类，提取出属于表面散射的冰雪消融变化；文献[5]通过估计场景中的噪声功率建立统计量，将真实变化区域与低信噪比区域区分，去除场景中低信噪比引起的低相干虚警干扰；还有的文献利用图像中目标的散射特征，分别采用聚类，半监督甚至深度学习的方式进行场景分类，实现变化检测^[6–8]。以上各方法或是侧重于观测方式的选取，或是关注于场景本身的散射属性，没有将二者有机地结合起来考虑，因而在适用范围上有一定的局限性。CCD方法如同使用一把具有精细刻度的尺子，能够丈量微小的变化，传统CCD方法之所以产生低相干虚警干扰，是因为这把尺子的尺度虽精细但单一，能够同时观测到目标变化与虚警干扰但无法对其进一步区分。将CCD方法在频率维度上扩展，采用不同频率电磁波对场景观测，利用场景中目标对多波段雷达信号的散射差异形成多尺度观测，使得对目标变化区域与虚警干扰区域形成不同的观测结果，从而将其有效区分。目前进行多波段CCD研究的文献较少。文献[9,10]采用极化变化检测的方法处理C, L两个波段的农作物场景SAR图像，文中所采用的方法将两个波段的数据结合起来，可检测出的变化区域为两个单独波段变化检测结果的并集。文献[11]使用X, C, S, L不同波段的SAR数据对北极地区分别进行CCD与非相干变化检测(NCCD)处理，通过对比实验结果得出结论：(1)CCD比NCCD能检测出更微小的变化；(2)NCCD处理结果与所使用的波段关系不大；(3)CCD中，观测电磁波波长越短(频率越高)，变化区域去相干程度越高。该文献只是对实验结果做初步定性分析，并指出：需要进一步对检测出的变化进行分类才具有应用意义。

在多波段观测模式下，不同波段的电磁波对场景探测，获得的目标散射特性存在较大差异，这与电磁波在不同频率下的穿透性能，可探测目标最小截面积，观测目标变化去相干程度等特性有关。不同类型的目标发生不同尺度的变化在各个波段探测下形成不同的相干变化检测结果。我们依此对各种变化分类，并从中提取我们感兴趣的目标变化区域。基于上述原理，本文充分分析不同波段电磁波的目标探测特性及其对不同尺度变化的去相干程度，据此提出一种多波段CCD方法。该方法首先获取场景在不同波段下的相干变化差异图像，然后根据目标在不同波段下的相干变化差异表征用改进的期望最大化(Expectation-Maximization, EM)算法进行分类，根据先验知识确定目标变化区域所属的类别，最后用Dempster-Shafer(DS)证据理论实现多波段相干变化差异图融合，提取目标变化区域。由于自然界中绝大多数场景对不同波段电磁波照射下呈现明显不同的散射差异^[12]，故采用本文提出的多波段CCD方法可适用于军事、农业、灾害监测等更广阔领域的场景中。通过实验数据处理结果表明，本文方法可以有效减小除目标变化区域外低相干干扰区域造成的影响。

本文的结构安排如下：第2节介绍不同波段电磁波探测目标特性的差异，第3节为多波段CCD算法原理与处理流程，第4节为实验数据验证，第5节对本文内容作总结。

2 目标的多波段相干性表征

场景中目标受电磁波照射，产生的响应为场景中各个目标散射中心响应之和，即目标反映在SAR图像上的复投影结果。我们参考文献[13,14]建立的参数模型，获取的SAR图像 $I$ 表示为：

其中， $E$ 为场景中目标总响应， $N_{\rm noi}$ 为噪声项。 ${E_i}$ 为场景中每个目标的散射中心响应，包括复散射系数项和相位项。其中，前者与电磁波频率 $f$ ，目标观测角度信息(俯仰角 $\theta $ 与偏航角 $\phi $ )及目标结构参数 $T\,$ 相关，后者与电磁波频率 $f$ 及目标散射中心距离参数 $R$ 有关。为了避免空间去相干对CCD结果的干扰，在数据采集过程中，我们保持场景变化前后的两次观测俯仰角 $\theta $ 与偏航角 $\phi $ 固定，使两次观测之间的复相干系数 $\gamma $ 不包含空间去相干。根据噪声去相干与变化去相干理论可知，影响复相干系数 $\gamma $ 的重要因素为SAR图像信噪比与发生变化的程度^[15]，这两点分别对应于场景中每个目标的复散射系数项与相位项，因此对于每个目标，其相干系数统计值可表示为 ${\hat \gamma _i}$ $ = {\hat \gamma _i}$ $\left( {f,{T_i},\Delta {R_i}} \right)$ 。

关于场景中的不同结构目标在不同频率电磁波探测下对相干系数的影响，已有一些文献从不同角度直接或间接地通过理论分析与实验验证得出结论^[15–19]，可总结如下：

(1) 高波段雷达可通过CCD方法检测小尺寸或表面粗糙程度较弱的目标造成的微小尺度变化，同时对大尺寸或表面粗糙程度强的目标造成的变化更具有良好检测性能；

(2) 低波段雷达通过CCD方法可检测大尺寸或表面粗糙程度较强的目标造成的较大尺度变化，对于小尺寸或表面粗糙程度弱目标的幅值响应弱，使得图像信噪比低，导致相干性差，易引起虚警；

(3) 中间波段雷达对目标的检测效果介于高、低波段之间，使得在目标尺寸结构或表面粗糙程度处于一定范围内的情况下，对于一定尺度目标变化的检测结果具有区分度；

(4) 低波段雷达独有穿透特性，可以检测深层或植被冠层以下目标发生的较大尺度变化。

根据上述结论可知，在场景中存在尺寸和表面粗糙度有较大差异的目标，且发生不同尺度变化的时候，得到的CCD结果会有许多值，可据此把场景中的目标分成许多类；每种目标在不同波段电磁波探测下得到不同大小的CCD结果，分别分属于每个波段的不同分类。我们利用这种分类及不同波段结果中分类的差异性设计多波段CCD方法，具体在第3节中阐述。

3 多波段CCD方法

根据上一节分析的目标在多波段探测下的相干性表征特性，本文提出一种多波段CCD方法。该方法能够充分利用不同波段电磁波探测目标及目标变化造成的散射差异进行相干变化检测，通过对多波段相干变化差异图分类，结合目标在不同波段探测下的散射特性构成的先验知识提取目标变化区域，排除CCD结果中的虚警。多波段CCD方法的主要步骤如下：

步骤1 对场景多波段SAR图像进行基于SIFT特征的图像配准与辐射校正，使之满足CCD处理需求^[20]。

步骤2 分别计算配准后各个波段变化前后SAR图像的相干变化统计量(一般为相干系数)，形成各波段的相干变化差异图。

步骤3 对每个波段分别处理：以该波段的相干变化差异图为数据，建立基于相干变化检验统计量概率密度函数(probability density function, pdf)的有限混合模型(Finite Mixture Modem, FMM)，用改进的EM算法自适应求解FMM中分模型数量 $n$ 及各个分模型pdf参数，据此在该波段对检测场景按照变化情况分类，形成该波段观测下的分类结果，并确定场景在该波段观测下对每种分类的响应度。

步骤4 在各个波段中，根据场景待检测区域中感兴趣地物的少量监督样本在该波段相干变化差异图的直方统计结果确定目标变化区域与非变化区域所属分类，然后根据响应度结果采用DS证据理论进行多波段相干变化差异图像融合，获取最终的多波段CCD结果。

基于上述处理步骤形成的多波段CCD方法流程图如图1所示。

上述处理过程中，步骤1和步骤2已有较多文献进行相关方法研究，可以分别参考文献[21]与文献[1]中的算法进行处理。下面针对能够体现多波段CCD方法特点的步骤3和步骤4进行详细讨论。

3.1 变化检测场景分类

根据第2节的理论，受雷达频率，目标尺寸，表面粗糙程度，变化尺度等因素影响，使用不同波段探测场景中的目标获取的CCD检测结果中存在相干程度的差异。因此，需要根据每个波段获取的相干变化差异图对检测场景的变化情况依相干度大小进行分类，为下一步根据分类结果实现多波段融合提供可靠的数据。

在每个波段数据中，我们需要借助步骤2采用相干变化检验统计量的pdf进行场景分类。这里我们不妨采用传统的相干系数检验统计量，其pdf表示为^[22]：

式(2)中， $\left| \gamma \right|$ 表示相干系数的真实值， $N$ 表示像元统计个数(统计窗口大小)， ${}_m{F_n}\left( { \bullet ; \bullet ; \bullet } \right)$ 为超几何分布， $\left( {\left| \gamma \right|,N} \right)$ 可作为式(2)所示pdf的一组参数。传统CCD只是根据检验结果相干性的高低将场景分成变化/非变化两类，本文描述的多波段CCD中，分类情况更为复杂：受目标结构影响，场景未发生变化时 $\left| \gamma \right|$ 值并非固定在1附近，可能因信噪比差异而在 $\left[ {0,1} \right]$ 区间内取值；受目标变化程度影响，场景发生变化时 $\left| \gamma \right|$ 值并非固定趋近于0，也可能在 $\left[ {0,1} \right]$ 区间内取值。因此根据不同目标散射特性与变化程度，依据相干统计结果可将场景分为许多类。每一类受目标纹理信息影响，统计像元数据之间的相关性可导致 $N$ 的实际值与名义值不符，且场景中不同分类的 $N$ 值也可能不同。据此，我们将某一波段观测下检测场景的相干变化差异图看作是基于式(2) pdf分布的FMM，将每个单一分布看作一个分类(分模型)，总体分布表示为：

其中， ${a_i}$ 为第 $i$ 分类的加权系数，满足 $\displaystyle\sum\nolimits_{i = 1}^n \!\!{{a_i}} = 1$ ，与第 $i$ 个分布的形状参数 $\left( {\left| {{\gamma _i}} \right|,{N_i}} \right)$ 一起构成第 $i$ 个分布的分模型参数。 $n$ 为FMM分类个数。本文以EM算法为基础实现基于FMM的变化场景分类。

EM算法是一种极大似然估计算法，它先构建初始分类，通过迭代的方式用最大似然准则调整样本的类别归属，其中E步计算场景中每个像元对当前模型参数中的每个分模型的响应度，M步根据E步计算的响应度估计每个分模型新的模型参数。这样循环迭代，直至模型参数满足收敛条件，实现混合分布模型中各个分模型的参数估计。将传统EM算法用于本文应用中存在以下两个问题：一是分模型个数 $n$ 固定，而本文应用中，分模型的个数未知；二是在EM算法的M步对分模型参数估计时，无法像高斯混合分布模型那样建立关于分布参数最大似然估计的显式表达式，也无法获取本文分布的对数累积量显式表达式^[23]。对此，本文方法做出如下改进。

针对第1个问题，本文采用自适应估计方法。首先，建立一个足够多分模型个数的混合分布模型。然后在迭代过程中，如果某一个分模型的加权系数 ${a_i}$ 小于某个阈值，则可认为该分类无意义予以去除，并重新计算剩下分模型的加权系数；如果有两个分模型参数 $\left| {{\gamma _i}} \right|$ 之差小于某个范围时，可认为它们属于同一分类，将其合并，具体操作为：将两个分模型加权系数相加，取加权系数较大的参数 $\left| {{\gamma _i}} \right|,{N_i}$ 共同作为合并后的分模型参数。这样实现混合模型个数的自适应估计。

针对第2个问题，本文在M步中采用中心矩估计量估计分模型参数。对于式(2)所示的分布，其各阶原点矩表达式^[22]为：

对于统计数据，我们用样本加权平均方法计算每个分模型的样本 $k$ 阶中心矩 $\bar \beta _i^k$ ，计算公式为：

其中， $\xi \left( {i,j} \right)$ 为第 $j$ 个样本像元在第 $i$ 个分模型处的响应度。我们通过式(6)所示的极值计算方法获取分模型参数。

在实际处理中，我们建立关于参数 $\left| \gamma \right|$ 和 $N$ 的 $k$ 阶中心矩2维索引图代入式(6)，从而避免了大量复杂耗时的计算。

3.2 多波段融合

我们将3.1节获取的变化检测中第 $f$ 个波段的第 $i$ 个分类结果记作 ${\rm{Cl}}_i^f\left( {a_i^f,\left| {\gamma _i^f} \right|,N_i^f} \right)$ ，为便于下文分析，在各个波段中，不妨将每个分类结果按照 $\left| {\gamma _i^f} \right|$ 由小到大顺序排列，即令 $0 \le \left| {\gamma _1^f} \right| < \left| {\gamma _2^f} \right| < ·\!·\!· \le 1$ 。接下来采用图像融合方法获取多波段融合后的相干变化差异图。

多波段图像融合方法主要有基于DS证据理论^[24]和基于模糊集理论^[25]的融合方法，它们都具有不确定信息的处理能力。针对本文对多波段差异图像融合的具体应用，我们采用没有对分类状态响应度函数作近似的DS证据理论融合方法。

在进行多波段融合之前，我们先根据DS证据理论框架构建变化状态空间，并需要根据一定的先验知识确定各个波段中每个分类隶属的变化状态。根据3.1节描述，每个分类所属变化状态与该波段探测下目标的信噪比以及变化程度有关，在这些参数难以预估的情况下，选择少量监督样本可以辅助确定分类所属的变化类型^[3]。监督样本的选择有两个准则：一是选择的样本区域在场景中为具有一定代表性的地物，二是所选择区域具有明确的变化状态(变化/不变化)，且这两种状态均有与之相符的样本。根据监督样本所属的变化状态及其相干变化差异图的直方图统计结果确定变化类集合 $\left\{ {{\rm{Cl}}_i^f\Bigr| {{c_l} \le i \le {c_h}} } \right\} \subseteq C\,$ 与非变化类集合 $ \left\{ {\rm{Cl}}_i^f\left| {u_l}\right. \right.$ $ \le i \le {u_h} \Bigr\} \subseteq U\,$ 。对于可能存在的 $\left\{ {{\rm{Cl}}_i^f\Bigr| {{c_h} < i < {h_l}}} \right\}$ 分类集合，其 $\left| \gamma \right|$ 值介于变化与非变化之间，它们可能属于低信噪比的非变化情况，也可能属于一定程度变化去相干的微小变化情况，在没有先验知识的情况下，令 $\left\{ {{\rm{Cl}}_i^f\left| {{c_h} < i < {h_l}} \right.} \right\} \subseteq \varTheta = \left\{ {C,U} \right\}$ 表示不确定情况。剩下其它类集合表示为 $\varnothing $ 。这样构成了关于变化状态的完备空间 ${2^\varTheta } = \left\{ {C,U,\varTheta ,\varnothing } \right\}$ 。令场景中的每个像元对状态空间中每种状态的概率分配函数等于对该状态类别的响应度，这样可保证每个像元中对所有状态的概率分配函数之和为1。根据证据合成方法，变化状态的多波段融合结果表示为：

其中， ${m_f}\left( \cdot \right)$ 表示第 $f$ 个波段的概率分配函数， $ \oplus $ 表示DS证据融合操作。这样对场景中逐个像元求解，即可获取多波段融合相干变化差异图。

需要说明的是，所有波段电磁波对水面等几乎完全镜面反射区域的散射强度都比较弱，一般均低于噪声水平，因此在各个波段的相干变化差异图中均呈现低相干特性。对这类区域的变化检测超出了本文多波段CCD方法的检测范围，因此我们假定该区域没有发生变化，将其看作虚警干扰从SAR幅度图像上做掩模处理将其滤除，即：在所有波段变化前后SAR图像上，幅度均低于某一阈值(噪声水平)的像元，所在区域认为未发生变化，在最终的多波段CCD结果中将该位置数据置为1。

4 实验数据验证

采用ESAR机载SAR数据对本文方法进行实验验证。该数据变化前后采集时间间隔为1个月(2007年2月－2007年3月)，地点为瑞典Eggby小镇郊外的一片场景。对该场景分别用L波段与P波段HH极化雷达进行变化前后数据采集，所得SAR图像及对照光学图像如图2所示。在变化前后的两次数据采集时间间隔内，场景中发生土壤翻修，草地修剪，冰雪消融等类型的目标变化，主要发生在场景中草地，裸地的部分区域以及湖面区域，在图2(e)的光学图像中将部分变化区域予以简要标注。经过图像配准后对每个波段SAR图像进行CCD处理，得到的相干变化差异图如图3所示。比较图3中两幅图可发现，场景中左侧大面积植被区域在L波段差异图中呈现低相干特性，在P波段差异图中为高相干特性，这是因为相比于L波段，P波段电磁波具有良好的穿透特性，没有受到植被冠层受风吹雨淋发生的变化去相干影响；场景中局部裸地、草地区域(红色圆圈部分)在L波段中呈现高相干特性，而在P波段中呈现低相干特性，这是因为该处场景表面粗糙度较弱，导致波长较长的P波段对该处的目标散射强度较弱，造成低信噪比，从而受到噪声去相干的影响。这两部分区域均不是我们感兴趣的目标变化区域，在各自波段的检测结果中形成虚警干扰。

基于图3所示两个波段的相干变化差异图，我们用3.1节提出的FMM及EM算法对场景分类。结合SAR图像与差异图信息，我们初始制定分模型数目均为5, M步迭代中取式(6)参数 $k = 2$ 。经过迭代后，L波段与P波段结果中分模型数目分别为4和3。图4分别给出了两个波段差异图的统计直方图与建立的FMM以及分模型的最终迭代结果，可见FMM对统计直方图的拟合效果良好。

在图3的两幅相干变化差异图中，我们分别选取黄色和绿色方框区域作为变化区域与非变化区域的样本监督数据，并将各自波段的直方图与各自的拟合分布模型比较，结果如图5所示。为方便比较，图5中将样本监督数据直方图乘以相应的倍数，使其与拟合分布分模型的概率密度大小相当。通过比较我们得知，L波段数据中，变化类与非变化类分别对应第1分模型与第3、第4分模型；P波段数据中，变化类与非变化类分别对应第1分模型与第3分模型。根据3.2节处理方法可得，两个波段的第2分模型分别对应各自的不确定类。然后计算像元在各个分类中的响应度(即计算概率分配函数)，代入式(7)计算，求得多波段融合相干变化差异图结果 $\left( {{m_P} \oplus {m_L}} \right)\left( C \right)$ 如图6(c)所示(场景中无因散射强度弱以致掩模去除区域)。图6(a)与图6(b)为使用模糊集理论的融合方法^[25]，前者使用加权算术平均的方式融合，后者使用加权几何平均方式融合。结合图3与图6可见，相比于单一波段CCD结果，多波段融合结果中，既排除了L波段中植被受风吹扰动造成的低相干干扰，又免受P波段中低散射强度导致低信噪比造成的低相干干扰影响，使得结果中只剩下土壤翻动，草地修剪，湖面冰雪消融造成的目标变化区域，更加接近真实值。在模糊融合结果与证据融合结果比较中，比较图6(a)–图6(c) 3种多波段融合结果，可见加权算术平均的模糊融合结果在检测阈值大于0.5情况下存在较多虚警，其它两种融合结果均具有较大的阈值选取动态范围，具体检测性能需要通过详细指标分析。

为了进一步说明本文使用的多波段CCD方法的性能，采用受试者工作特征(Receiver Operating Characteristic, ROC)曲线与最优Kappa系数两个指标进行分析。ROC曲线越靠近左上角，最优Kappa系数值越大，说明方法性能越好。图7(a)，图7(b)分别展示了实验采用的两个单波段CCD与多波段CCD方法的ROC曲线和选取不同检测门限的Kappa系数。从图中可以看出，本文使用的DS证据融合多波段CCD方法的ROC曲线更接近左上方，最优Kappa系数值为0.81，大于两个单波段CCD方法(L波段：0.63，P波段：0.54)和两种模糊融合方法(均为0.79)。因此，以上分析验证了本文提出的多波段CCD方法对去除单波段CCD方法中低相干干扰，降低虚警率方面的正确性与有效性。

5 结束语

本文根据场景中目标及其发生的变化在多波段电磁波观测情况下表征出的相干性差异，提出一种多波段CCD方法。该方法先将待检测场景根据各波段的相干变化差异图用改进EM算法进行各自的自适应分类处理，然后依据少量监督样本分别在各个波段中确定分类的变化状态，最后使用DS证据理论方法处理，获取多波段融合CCD结果，去除了每个波段各自产生的低相干干扰，达到降低虚警率的目的。由于用于多波段CCD处理的重轨SAR数据较少，本文中仅用一组双波段数据验证该方法的有效性与正确性。接下来的研究应以更多波段对不同场景探测的重轨SAR数据为基础，进一步分析波段种类、波段数目、场景地物类型等因素对多波段CCD方法实现效果的影响。

致谢　本文作者感谢欧空局提供ESAR机载L波段与P波段重轨SAR数据。