在过去几十年里，遥感变化检测方法被广泛应用于各种遥感领域（如林业、农业、城市和灾害监测）^[1-7]。根据遥感变化检测的目的，可以将其分为2类：二元变化检测（目的是区分变化类和无变化类）和多元变化检测（目的是区分不同类型的变化）^[8]。与二元变化检测相比，多元变化检测在处理多地物影像时，变化类型数量的增加会提高计算复杂性和变化类型估计难度^[9]。

在当前的遥感变化检测方法中，分类后变化检测（post-classification comparison，PCC）方法通过对不同时相的遥感图像进行分类来识别变化信息，但易受到累积分类误差的影响。而变化向量分析法（change vector analysis，CVA）^[10-11]由于对图像的辐射校正要求严格，且易受到多种干扰因素的影响，应用受到了限制。为克服这些限制，Chen等^[12]提出了后验概率空间变化向量分析法（change vector analysis in posterior probability space，CVAPS），它放宽了对图像辐射校正的要求，并具有更高的变化检测性能。并且，CVAPS可以根据后验概率变化向量的方向，区分不同类型的变化，实现多元变化检测。然而，CVAPS利用支持向量机（support vector machine，SVM）^[13-14]来估计遥感图像像素的后验概率向量，易受到遥感图像中同物异谱、异物同谱、混合像元等因素的影响，从而难以准确估计复杂像元后验概率向量的强度和方向，并影响了其后多元变化检测的精度^[15]。为解决这一问题，李轶鲲等^[16-18]结合模糊C均值聚类（fuzzy C-means，FCM）和上下文敏感的贝叶斯网络（context-sensitive Bayesian network，CSBN），提出了FCM-CSBN-CVAPS方法。该方法通过建立像素与多种地物之间的随机链接，能够有效处理同物异谱、异物同谱及混合像元等问题，从而提高了变化检测的准确性和鲁棒性。因此本文基于CSBN估计后验概率变化向量，实现了更为精确的多元变化检测。在多元变化检测方法中，SVM-CVAPS根据后验概率变化向量与标准变化向量的最大相关程度（maximum correlation，MC）来判断变化类型^[15]。然而，较多的地物类型会生成过多的标准变化向量，这些标准变化向量互相干扰，大大降低了多元变化检测的精度。

因此，本文提出了一种多模式结合的多元变化检测算法FCM-CSBN-CVAPS-AT。该方法首先使用FCM捕获遥感影像中的信号类，再通过CSBN基于像素级空间关系建立每个像素和地物类型间多对多的随机链接，并利用CVAPS模型计算得到每个像素的后验概率变化向量。在此基础上，提出一种角度阈值（angle thresholding，AT）法。AT法通过计算像素后验概率变化向量与某一标准变化向量的夹角，当夹角小于一定阈值时，则判定该像素为该标准变化向量所代表的变化类型。由于AT法不需要在不同标准变化向量间进行比较，有效地避免了标准变化向量互相干扰的问题。实验证明，与比较算法相比，本文算法取得了更高的多元变化检测精度，且不易受算法参数的影响，具有一定的实用价值。

高分辨率遥感（high resolution remote sensing，HRRS）影像普遍存在的混合像元现象往往影响变化检测性能。因此，需要使用FCM算法将像素分解为多个信号类。与一般解混方法检测到的端元不同，信号类是具有相似光谱特征的像素簇。此外，引入空间信息可以有效地处理HRRS影像的光谱可变性，减轻不确定性。因此，CSBN通过引入像素级空间信息，建立信号对和地物类型之间的多对多随机链接。一方面，CSBN可以将不同的信号类（即不同的光谱）链接到同一种地物类型，以处理光谱可变性；另一方面，CSBN可以根据空间上下文信息将同一信号类链接到不同的地物类型，针对光谱和地物间一对多的关系建模。

基于CSBN，本文提出了一种多模式结合的变化检测算法（FCM-CSBN-CVAPS-AT），其总体流程如图1所示。首先，通过FCM提取两时相遥感影像中的信号类，并通过CSBN估计双时相遥感图像的像素级后验概率向量，并在此基础上使用CVAPS生成后验概率变化强度图；其次，使用Otsu法^[19]计算变化阈值，生成变化二值图，再经过形态学处理，去除细小相斑和孔洞，生成优化后的变化二值图；最后，针对变化区域，利用AT法计算像素后验概率变化向量与某一标准变化向量的夹角，当夹角小于一定阈值，则判定该像素为标准变化向量所代表的变化类型，从而实现HRRS影像的多元变化检测。由于FCM-CSBN-CVAPS将空间信息融入后验概率向量的估计过程，并使用后验概率向量而不是使用分类图进行变化检测，因此，提出的算法对累积分类误差更具鲁棒性。

CSBN可以将空间信息与后验概率估计结合，因此本文使用CSBN来进行像素后验概率估计。然而，由于混合像元和累积聚类误差的影响，K-means无法有效捕获遥感图像中的信号类，大大地限制了变化检测性能。因此，本文改进了CSBN模型，将其与FCM相结合，在像素后验概率估计和模型训练阶段引入像素级空间信息。

本文构建了4层CSBN模型以计算像素p_i，j关于地物类型L_v的后验概率P（L_v|p_i，j）（图1）。第1层为遥感图像的像素p_i，j；第2层为像素对（p_i，j，p_g，h），其中p_g，h位于中心像素p_i，j的邻域N_i，j内；第3层为信号对（ω_m，ω_n），其中，m和n是信号类ω_m和ω_n的索引；第4层为遥感图像中地物类型。P（L_v|p_i，j）的计算公式为：

式中：C为信号类数量；P（m，n）是信号对（ω_m，ω_n）的先验概率，而条件概率P（（m，n）|L_v）则是信号对（ω_m，ω_n）和地物类型L_v间的随机链接。条件概率P（（m，n）|L_v）可以从训练样本中获得。本文假设先验概率P（L_v）服从一致分布。具体推导步骤请参见参考文献[16]。

通过相应的模糊隶属度可以计算概率P（（m，n）|（${{{p}_{i}}_{，}}_{j}$，p_g，h）），公式为：

式中P（m|p_i，j）和P（n|p_g，h）分别为像素p_i，j和p_g，h属于信号类ω_m和ω_n的条件概率，可基于模糊隶属度u_m（i，j）和u_n（g，h）估计。概率P（（p_i，j，p_g，h）|p_i，j）的大小反映了像素对（p_i，j，p_g，h）在估计后验概率P（L_v|p_i，j）时的重要程度。如果邻域N_i，j中的每个像素p_g，h对P（L_v|p_i，j）的估计同样重要，则P（（p_i，j，p_g，h）|p_i，j）=1/|N_i，j|，其中|N_i，j|为N_i，j中的像素数。

为了计算后验概率P（L_v|p_i，j），必须基于地物类型L_v的训练集T_v估计获得概率P（（m，n）|L_v）。估计公式基于信号对（ω_m，ω_n）出现在地物L_v中的频率SPF_v（m，n），其定义为：

由于模糊隶属度u_m（x，y）和u_n（g，h）的取值范围是0~1，因此需要对训练集T_v中所有像素对（${{{p}_{x}}_{，}}_{y}$，p_g，h）对应信号对（ω_m，ω_n）的模糊隶属度值求和，得到信号对频率SPF_v（m，n）。

基于各个信号对（ω_ϑ，ω_μ）（ϑ，μ = 1，2，…，C）的频率SPF_v（ϑ，μ），条件概率P（（m，n）|L_v）可以近似为：

进一步，为计算式（4）中的P（（m，n）|L_v），则先验概率P（（m，n））的计算采用以下总概率公式：

其中P（（m，n）|L_v）由式（3）和式（4）估计得到。

设p_i，j表示给定图像的第i行和第j列中的像素。当混合像元p_i，j由地物类型L₁和L₂组成时，遥感图像中像元在t₁和t₂时刻的光谱变化较小，导致其在两时相图像中估计的后验概率向量变化不显著。假设在t₁时刻，后验概率P（L₁|p_i，j）和P（L₂|p_i，j）分别为51%和49%，根据最大后验概率（maximum a posteriori，MAP）原则，将像素p_i，j分类为L₁；假设在t₂时刻，后验概率P（L₁|p_i，j）和P（L₂|p_i，j）分别为49%和51%，像素p_i，j被归类为L₂。因此，PCC方法将像素p_i，j识别为变化的像素，而不考虑其在两时相图像中所对应的微小光谱差异。因此MAP分类错误会引起PCC的误判。

为了解决这个问题，Chen等^[12]提出了CVAPS框架，分别用ρ¹=（${\rho }_{1}^{1}$，…，${\rho }_{v}^{1}$，…，${\rho }_{s}^{1}$）和ρ²=（${\rho }_{1}^{2}$，…，${\rho }_{v}^{2}$，…，${\rho }_{s}^{2}$）表示像素在t₁和t₂时刻的后验概率向量。${\rho }_{v}^{1}$和${\rho }_{v}^{2}$分别表示t₁和t₂时刻的后验概率P（L_v|p_i，j），s表示两时相图像中地物类型的数量。因此，像素p_i，j的变化向量Δρ定义为：

像素p_i，j的变化幅度‖Δρ‖定义为：

本文根据式（7）生成后验概率变化幅度图，经过自动阈值算法处理，生成变化二值图。由于CVAPS比较两时相后验概率向量而不是分类结果，因此显著降低了累积分类误差。此外，与CVA相比，CVAPS将类内和类间距离进行归一化处理，使得不同类型的变化幅度‖Δρ‖在同一尺度内，因而适合使用单一阈值有效确定变化与非变化像素，无需采用复杂的多阈值算法。

FCM算法是一种基于目标函数的模糊聚类算法，能将单一像元以不同隶属度分解成不同信号类，从而建立起像元与信号类间的随机链接，分解混合像元。假设存在t时相M行N列像素构成的遥感影像I，将影像I划分为C个信号类，u_k（i，j）（1≤k≤C）是像素${{{p}_{i}}_{，}}_{j}$关于信号类ω_k的隶属度。FCM聚类算法的目标函数定义为：

式中：J为目标函数；U为由U_k组成的模糊矩阵；ψ={Ω₁，…，Ω_k，…，Ω_C}是信号类中心的集合；Ω_k为信号类ω_k的中心；q为模糊度参数。

多变化类型问题可以表现为如下形式：设Ω={Ω_n，Ω_c}为所有类的集合，其中Ω_n为不变类，Ω_c={Ω_c1，Ω_c2，…，Ω_ck}是k种可能变化类型的集合。多变化类型任务的目标即为检测变化像元Ω_c，并识别它们在{Ω_c1，Ω_c2，…，Ω_ck}中的类型。

首先，使用E_d来表示仅属于类d的标准像素后验概率向量，表达式为：

如果一个标准像素在一段时间内从a类变为b类，则其标准变化向量可表示为：

设t₁和t₂时刻在同一像素上获取的后验概率变化向量为Δρ。在多变化类型检测中，SVM-CVAPS根据后验概率变化向量与标准变化向量的MC来判断变化类型。MC法使用标准变化向量集SE={ΔE_ab|a，b=1，2，…，m，a≠b}确定Δρ的变化类型。假设标准变化向量集SE中的ΔE_qr与Δρ的夹角最小，表明ΔE_qr与Δρ的变化类型相似。因此，Δρ的变化类型被确定为“从q到r”。所有标准变化向量ΔE_ab具有相同的范数值，即最小夹角相当于内积的最大值，那么，就可以通过式（12）来确定变化类型：

假设在两时相遥感影像中共有N类地物，每一个像素的后验概率变化向量需要与N×（N-1）个标准变化向量比较，并根据MC确定变化类型。然而，较多的地物类型会生成过多的标准变化向量，这些标准变化向量互相干扰，大大降低了多元变化检测的精度。

针对上述MC法的问题，本文在MC法提出的标准变化向量基础上，根据夹角阈值确定变化类型。由于后验概率变化向量Δρ和标准变化向量ΔE_ab间夹角θ可以有效反映不同类型的变化信息，因此本文提出的AT法使用夹角θ来区分多种变化类型。夹角θ的计算公式为：

假设ΔE_qr（ΔE_qr∈SE）与Δρ的夹角θ小于阈值角度T（T=30°），则表示ΔE_qr与Δρ的变化方向相似，因此确定Δρ的变化类型为“从q到r”。由于AT法仅仅根据单一标准变化向量ΔE_qr确定像素的变化类型，避免了过多标准变化向量互相干扰的问题，显著提高了变化类型判断精度。

本文使用的影像数据来自SenseEarth平台数据集（地址为： https://rs.sensetime.com/），空间分辨率为3 m，影像尺寸为512像素×512像素。遥感影像经过了辐射校正、几何纠正和图像配准预处理操作。研究区1实验影像共包括1类地物变化类型，为林地到荒地；研究区2实验影像共包括2类地物变化类型，分别为建筑物到荒地、建筑物到林地；研究区3实验影像共包括3类地物变化类型，分别为荒地到建筑物、林地到建筑物、建筑物到荒地。各研究区总地物类型数量和变化类型数量如表1所示。算法运行环境为 13th Gen Intel Core i7-13700KF 3.42 GHz，32 GB内存，16核24线程处理器。为了验证方法的有效性，本文选用FCM-SBN-CVAPS-AT，SVM-CVAPS-AT，FCM-CSBN-CVAPS-MC，FCM-SBN-CVAPS-MC，SVM-CVAPS-MC，FCM-CSBN-PCC和DCVA^[20]算法作为对比算法，同本文提出的FCM-CSBN-CVAPS-AT算法进行比较。

本节实验中用到的FCM-CSBN-CVAPS-AT，FCM-SBN-CVAPS-AT，SVM-CVAPS-AT，FCM-CSBN-CVAPS-MC，FCM-SBN-CVAPS-MC，SVM-CVAPS-MC和FCM-CSBN-PCC算法的运行参数为：聚类数为10，模糊参数q=3.5，每类训练样本数为1 000；DCVA算法无运行参数；AT法的角度阈值T=30°。

综合考虑研究区内地物变化，本文关注的地物变化类型主要有以下5种：建筑物到荒地，建筑物到林地，林地到荒地，林地到建筑物，荒地到建筑物。图2—4分别为研究区1、研究区2和研究区3的遥感影像及其检测变化结果。通过对比可以发现：① FCM-CSBN-CVAPS-AT算法性能明显优于SVM-CVAPS-AT，FCM-SBN-CVAPS-AT，FCM-CSBN-PCC和DCVA算法；② FCM-CSBN-CVAPS-MC算法性能明显优于SVM-CVAPS-MC，FCM-SBN-CVAPS-MC和FCM-CSBN-PCC算法。这证明了FCM-CSBN-CVAPS算法在多元变化检测中的优越性。并且，在绝大多数情况下，各个算法的AT版本精度均优于MC版本，证明了AT法的优势。

研究区1的地物变化类型为林地到荒地。从图2（d）和（e）可以看出，SVM-CVAPS-AT算法无明显漏检区域，但边界较为模糊，漏检区域较大；SVM-CVAPS-MC算法漏检严重，边界更为模糊。从图2（f）和（g）可以看出，FCM-SBN-CVAPS-MC算法虽检测出较为清晰的边界区域，但存在一定面积错检区域；FCM-SBN-CVAPS-AT算法虽也存在错检区域，但边界区域较FCM-SBN-CVAPS-MC法更为清晰。从图2（h）可以看出，FCM-CSBN-PCC算法相较于其他算法漏检范围小，但存在大量错检区域。从图2（i）可以看出，DCVA算法检测结果边界较为模糊，存在大量错检区域。从图2（j）和（k）可以看出，与MC法相比，FCM-CSBN-CVAPS-AT算法的错检、漏检范围小，较好地保留了变化区域的边界信息，取得了最好的变化检测结果。

研究区2的地物变化类型为：建筑物到林地、建筑物到荒地。从图3（d）—（i）可以看出，与研究区1的结果类似，各个对比算法均未取得满意的检测结果。从图3（j）和（k）可以看出，FCM-CSBN-CVAPS-MC算法虽可以检测出边界信息，但仍存在少量漏检区域；而FCM-CSBN-CVAPS-AT算法能较好地检测和保留变化区域边界细节，有效检测到2种变化类型，具有最高的精度。

研究区3的地物变化类型为：荒地到建筑物、林地到荒地和建筑物到荒地。从图4（d）—（i）可以看出，与研究区1和2类似，各个对比算法无法取得满意的检测结果。从图4（j）可以看出，FCM-CSBN-CVAPS-MC算法存在明显的错检和漏检区域，识别建筑物的变化类型时，将大块区域的荒地到建筑物检测为建筑物到荒地，将未变化的水体检测为林地到建筑物，未能正确识别出变化类型。从图4（k）可以看出，FCM-CSBN-CVAPS-AT算法只存在少量漏检和错检区域，可以较完整地识别出变化区域和变化类型，检测效果最好。

本实验以FCM-CSBN-CVAPS-AT算法检测研究区2的变化类型为例，测试角度阈值对本文算法的影响。角度T取6个值，即10°，20°，30°，40°，50°和60°，如图5所示，在T为30°时取得最高Kappa系数，当T大于或小于30°时，Kappa值均小于0.771 2。此外，T为10°时未检测出变化区域。由此可见，过高或过低的角度阈值都会导致Kappa值降低。

本文使用错检率、漏检率、总体精度和Kappa系数4个指标对本文方法和对比方法的多变化检测性能进行比较分析。

从表2—3的实验结果可以看出，本文算法虽未完全取得最低错检率和漏检率，但总体精度和Kappa系数均高于其他对比算法。以研究区2的建筑物到林地变化类型为例（表2），本文算法比FCM-SBN-CVAPS-AT算法的总体精度和Kappa系数分别高出3.46%和0.230 7，比SVM-CVAPS-AT算法的总体精度和Kappa系数分别高出2.08%和0.281 7。

由于DCVA算法预训练网络精度不适用于本文实验影像，所以本文算法的总体精度和Kappa系数分别比DCVA算法高出28.36%和0.656 4。

如表2所示，本文算法的Kappa系数比FCM-CSBN-PCC算法高出0.244 2，证明CVAPS比PCC更不易受累积分类误差的影响。此外，FCM-CSBN-CVAPS-AT和FCM-SBN-CVAPS-AT算法分别比FCM-CSBN-CVAPS-MC和FCM-SBN-CVAPS-AT算法Kappa系数高出0.062 9和0.037 9，由此证明AT法比MC法具有更高的检测精度，在多元变化类型检测中更具优势。

SVM-CVAPS算法估计的后验概率向量具有较低不确定性，因而大大高估了变化区域。从研究区2和3的实验结果中可以看出，SVM-CVAPS算法具有较高的漏检率和错检率，并且其在多类地物变化检测中无法同时正确识别多种地物类型变化情况。相比之下，FCM-CSBN-CVAPS使用CSBN合理引入了空间信息，其估计后验概率向量虽然也具有较低的不确定性，但没有明显高估变化区域，检测的变化区域最为准确且边界较为平滑。

由于地物光谱的低可分性和地物结构的高度复杂性，遥感图像的不确定性常常限制了基于遥感的多变化检测性能。本文所使用的FCM-CSBN-CVAPS-AT算法通过CSBN引入空间信息估计后验概率变化向量并根据后验概率变化向量的方向来判断变化类型，实现针对HRRS影像多变化类型检测。实验证明，与对比算法相比，本算法精度更高且不易受参数影响，具有更高鲁棒性。

尽管本文所提出的算法取得了较好的多变化检测性能，但仍存在以下缺点：首先，确定模糊度q、聚类数和CSBN的窗口大小需要通过实验来完成；其次，CSBN的计算复杂度较高，计算时间明显长于对比算法。因此，未来工作将着重于开发一种自动确定最佳模糊度q、信号类数量和窗口大小的方法，并通过并行计算降低计算时间，以提高算法的实用性。