基于慢特征分析的遥感影像变化检测方法

摘要

本发明公开了一种基于慢特征分析的遥感影像变化检测方法，本发明将慢特征分析方法引入多时相遥感影像的变化检测中，慢特征分析方法充分考虑多时相遥感影像间的空间对应关系，并将针对连续信号分析的原始理论发展成基于离散数据集的变化检测算法。慢特征分析方法能提取多时相遥感影像数据集中的不变特征作为特征空间，在此特征空间中，多时相遥感影像间的光谱差异得到了抑制，真实变化得到了突出，从而可提高变化检测精度。同时，慢特征分析方法实现简单、运算速度很快。

说明书

技术领域

本发明属于遥感影像处理技术领域，尤其涉及一种基于慢特征分析的遥感影 像变化检测方法。

背景技术

人类活动对地球表面环境造成了巨大影响，这种影响体现在环境变化、城市 发展等各个方面。因此，实时精确地获取地球地表覆盖的变化情况对于环境监测 和资源管理意义重大（Lu,D.,P.Mausel,E.Brondizio and E.Moran(2004). ″Change detection techniques.″International Journal of Remote Sensing25(12): 2365-2401.）。

变化检测是指通过对同一地区不同时间的地物分布情况进行观测来确定地 表变化（Singh,A.(1989).″Review Article Digital change detection techniques using  remotely-sensed data.″International Journal of Remote Sensing10(6):989-1003.）。遥 感影像能够长期提供较大范围的地表信息，在变化检测中有着重要的应用。目前 遥感影像变化检测已经广泛应用于土地覆盖\土地利用监测、城市发展研究、生 态系统变化监测和灾害监测等领域。如何快速准确地提取遥感影像中的变化区域 是遥感影像研究的难点问题。

虽然遥感影像可以以一定周期性覆盖同一区域，但是由于成像条件（太阳角 度、传感器、大气条件、土壤湿度、物候差异）不同，造成同一地物在不同时期 的遥感影像上显示出不同的亮度值，这种现象称为光谱差异。光谱差异是影响遥 感影像变化检测精度的最主要问题之一。如何有效地消除光谱差异，充分利用遥 感影像的多波段信息提取真实变化，是遥感影像变化检测亟需解决的难题。

针对遥感影像变化检测问题，国内外学者提出了很多方法，包括影像代数法、 分类后变化检测法、光谱混合分析法、物理特征提取法等。其中，影像变换法是 研究最为深入、应用最为广泛的一种方法。影像变换法是将多波段遥感影像看作 多维数据集，通过对多维数据进行投影的方法，将原始遥感影像投影到特征空间 中，在特征空间中进行变化检测。影像变换法能够充分利用多波段遥感影像的光 谱信息，通过数学原理寻找到最利于提取变化信息的特征空间，抑制噪声和光谱 差异，提高变化检测精度。主成分分析变换算法、GS变换算法和多元变化检测 变换算法是最主要的影像变换法。其中，主成分分析变换算法和GS变换算法由 于只考虑了数据的统计分布性质，没有考虑变化检测问题的特性，因此效果并不 能让人满意，大多只应用于预处理等步骤，不能直接获得地表覆盖变化信息。多 元变化检测算法由Nielsen提出（Nielsen,A.A.,K.Conradsen and J.J.Simpson  (1998).″Multivariate Alteration Detection(MAD)and MAF Postprocessing in  Multispectral,Bitemporal Image Data:New Approaches to Change Detection  Studies.″Remote Sensing of Environment 64(1):1-19.），该算法基于数学统计学中 的典型相关分析理论，针对两组数据之间的相关性进行分析，从而对两个数据集 分别进行变换，得到变换特征。多元变化检测算法得到了广泛的研究，但是该算 法没有考虑到同源遥感影像之间波段的对应关系，无法充分利用遥感影像的光谱 信息。

慢特征分析算法是一种从快速变化的输入信号中提取不变特征的非监督学 习方法（Wiskott,L.and T.J.Sejnowski(2002).″Slow feature analysis:Unsupervised  learning of invariances.″Neural computation14(4):715-770.）。研究发现从影像序列 中训练得到的慢特征函数与初级视觉皮层中的细胞种群有着定量上和定性上的 一致性。慢特征分析已经在空间变换的不变物体识别、非线性盲信号分解和人类 动作识别上得到了深入研究和应用。对于多时相遥感影像，变化像元存在着多种 变化类型，因此变化像元的差异是多种多样的，而未变化像元的差异只包含由于 外界条件变化所造成的光谱差异，情况比较单一。相比较突出变化信息，通过抑 制未变化像元的光谱差异来提高变化检测精度更加有效和实际。因此，我们可以 认为，在多时相遥感影像数据集上，时间维上的不变特征空间中，光谱差异得到 了最大程度的抑制，变化像元得到了足够的突出。

本文涉及如下参考文献：

[1] L.Bruzzone and D.F.Prieto,″Automatic analysis of the difference image for  unsupervised change detection,″IEEE Trans.Geosci.Remote Sens.,vol.38,pp. 1171-1182,May 2000.

[2] A.A.Nielsen,″The Regularized Iteratively Reweighted MAD Method for  Change Detection in Multi-and Hyperspectral Data,″IEEE Trans.Image Process., vol.16,pp.463-478,Feb2007.

[3] R.J.Radke,S.Andra,O.Al-Kofahi,and B.Roysam,″Image change  detection algorithms:a systematic survey,″IEEE Trans.Image Process.,vol.14,pp. 294-307,Mar 2005.

[4] R.E.Kennedy,P.A.Townsend,J.E.Gross,W.B.Cohen,P.Bolstad,Y.Q. Wang,and P.Adams,″Remote sensing change detection tools for natural resource  managers:Understanding concepts and tradeoffs in the design of landscape  monitoring projects,″Remote Sens.Environ.,vol.113,pp.1382-1396,Jul 2009.

[5] T.Celik,″Unsupervised Change Detection in Satellite Images Using  Principal Component Analysis and K-Means Clustering,″IEEE Geosci.Remote Sens. Lett.,vol.6,pp.772-776,Oct 2009.

[6] F.Bovolo and L.Bruzzone,″A Theoretical Framework for Unsupervised  Change Detection Based on Change Vector Analysis in the Polar Domain,″IEEE  Trans.Geosci.Remote Sens.,vol.45,pp.218-236,Jan 2007.

[7] P.Coppin,I.Jonckheere,K.Nackaerts,B.Muys,and E.Lambin,″Digital  change detection methods in ecosystem monitoring:a review,″Int.J.Remote Sens., vol.25,pp.1565-1596,May 2004.

发明内容

针对现有技术存在的不足，本发明提供了一种精度高、运算快、效率高的基 于慢特征分析的遥感影像变化检测方法。

为解决上述技术问题，本发明采用如下的技术方案：

一种基于慢特征分析的遥感影像变化检测方法，该方法将多时相遥感影像X 和Y分别读入大小为N×P的矩阵X和Y中，矩阵中各元素为各波段对应的像 素辐射值，N为多时相遥感影像的波段数，P为多时相遥感影像的像素数，基于 矩阵X和矩阵Y对多时相遥感影像X和Y进行如下操作：

步骤1，获取多时相遥感影像X和Y各波段的像素辐射值均值和像素辐射 值标准差；

步骤2，根据各波段的像素辐射值均值和像素辐射值标准差标准化多时相遥 感影像X和Y；

步骤3，基于标准化多时相遥感影像X和Y获取多时相遥感影像X和Y的 差值影像的协方差矩阵A、以及多时相遥感影像X和Y的协方差矩阵的和矩阵 B；

步骤4，求解矩阵A相对矩阵B的广义特征值及广义特征值对应的特征向 量，按照广义特征值大小对对应的特征向量进行排序，获得有序的特征向量矩阵， 其中，将广义特征值从小到大排序，排序为j的广义特征值为第j个波段的广义 特征值；

步骤5，采用有序的特征向量矩阵，将标准化多时相遥感影像X和Y投影 到特征空间，并获得多时相遥感影像X和Y各波段的特征差值；

步骤6，根据各波段的特征差值和广义特征值，获得各像素点的卡方距离， 并根据像素点的卡方距离获得多时相遥感影像X和Y的变化检测结果。

,步骤2具体为：

采用公式标准化多时相遥感影像，其中，是多时相遥感影 像第i个像素第j个波段的标准化像素辐射值；是多时相遥感影像第i个像素 第j个波段的像素辐射值；是多时相遥感影像第j个波段的像素辐射值均值； 是多时相遥感影像第j个波段的像素辐射值标准差。

步骤3中所述多时相遥感影像X和Y的差值影像的协方差矩阵A为：

$A=\frac{1}{P}\underset{i=1}{\overset{P}{\Sigma }}({\hat{x}}^i-{\hat{y}}^i){({\hat{x}}^i-{\hat{y}}^i)}^T$

其中，和分别是多时相遥感影像X和Y第i个像素的标准化像素辐射 值向量，P是多时相遥感影像X的总像素数。

步骤3中所述的多时相遥感影像X和Y的协方差矩阵的和矩阵B为：

$B=\frac{1}{2P}[\underset{i=1}{\overset{P}{\Sigma }}\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T+\underset{i=1}{\overset{P}{\Sigma }}\left({\hat{y}}^i\right){\left({\hat{y}}^i\right)}^T]$

其中，和分别是多时相遥感影像X和Y第i个像素的标准化像素辐射 值向量，P是多时相遥感影像X的总像素数。

步骤5具体为：

将有序的特征向量矩阵分别与标准化多时相遥感影像X和Y相乘，并将相 乘后所得矩阵各对应元素分别相减，得到多时相遥感影像X和Y的特征差值影 像。

步骤6中所述的像素点的卡方距离其中，Tⁱ为第i个像素 的卡方距离；为第i个像素第j个波段的特征差值；λ_j为第j个波段的广 义特征值；N为多时相遥感影像的波段数。

步骤6中所述的根据像素点的卡方距离获得多时相遥感影像X和Y的变化 检测结果，具体为：

判断像素点的卡方距离与阈值的大小，卡方距离不小于阈值的像素点为变化 像素点，赋值1；卡方距离小于阈值的像素点为未变化像素点，赋值0；基于像 素点赋值分割多时相遥感影像X和Y并获得二值影像，即为多时相遥感影像X 和Y的变化检测结果。

另一种基于慢特征分析的遥感影像变化检测方法，包括步骤：

步骤1，以单位矩阵初始化像素权值矩阵；

步骤2，基于当前像素权值矩阵，获取多时相遥感影像X和Y各波段的像 素辐射值的加权平均值和加权标准差；

步骤3，根据各波段像素辐射值的加权平均值和加权标准差标准化多时相遥 感影像X和Y；

步骤4，基于标准化多时相遥感影像X和Y，并引入当前像素权值矩阵获取 多时相遥感影像X和Y的差值影像的协方差矩阵A、以及多时相遥感影像X和 Y的协方差矩阵的和矩阵B；

步骤5，求解矩阵A相对矩阵B的广义特征值及广义特征值对应的特征向 量，按照广义特征值大小对对应的特征向量进行排序，获得有序的特征向量矩阵， 其中，将广义特征值从小到大排序，排序为j的广义特征值为第j个波段的广义 特征值；

步骤6，采用有序的特征向量矩阵，将标准化多时相遥感影像X和Y投影 到特征空间，并获得多时相遥感影像X和Y各波段的特征差值；

步骤7，根据各波段的特征差值和广义特征值，获得各像素点的卡方距离， 并基于像素点的卡方距离和卡方分布更新像素点的权值，从而获得更新后的像素 权值矩阵，以更新后的像素权值矩阵为当前像素权值矩阵；

步骤8，判断本次迭代和上次迭代各波段的特征差值的最大差值，若最大差 值大于设定阈值或者本次迭代为首次迭代，重新执行步骤（2）～（7）；若最大差 值小于设定阈值，执行步骤（9）；所述的各波段的特征差值的最大差值为：针对 各波段，分别求取各波段在本次迭代和上次迭代中对应的特征差值的差值，所获 的差值的最大值即为各波段的特征差值的最大差值；

步骤9，根据像素点的卡方距离获得多时相遥感影像X和Y的变化检测结 果。

步骤4中所述的多时相遥感影像X和Y的差值影像的协方差矩阵A为：

$A=\frac{\underset{i=1}{\overset{P}{\Sigma }}\left[v^i({\hat{x}}^i-{\hat{y}}^i){({\hat{x}}^i-{\hat{y}}^i)}^T\right]}{\underset{i=1}{\overset{P}{\Sigma }}{v}^i}$

其中，和分别是多时相遥感影像X和Y第i个像素的标准化像素辐射 值向量；vⁱ是多时相遥感影像第i个像素的当前像素权值；P是多时相遥感影像 X的总像素数。

步骤4中所述的多时相遥感影像X和Y的协方差矩阵的和矩阵B为：

$B=\frac{1}{2\underset{i=1}{\overset{P}{\Sigma }}{v}^i}[\underset{i=1}{\overset{P}{\Sigma }}\left[v^i\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T\right]+\underset{i=1}{\overset{P}{\Sigma }}\left[v^i\left({\hat{y}}^i\right){\left({\hat{y}}^i\right)}^T\right]]$

其中，和分别是多时相遥感影像X和Y第i个像素的标准化像素辐射 值向量；vⁱ是多时相遥感影像第i个像素的当前像素权值；P是多时相遥感影像 X的总像素数。

步骤7中所述的基于像素点的卡方距离和卡方分布更新像素点的权值，具体 为：

遍历多时相遥感影像X和Y，计算以波段数N为参数的卡方分布中，各像 素点大于其卡方距离的概率，并将概率作为权值赋给对应的像素点。

与现有技术相比，本发明具有以下特点和有意效果：

1、基于慢特征分析，将多时相遥感影像投影到不变特征空间中，在此特征 空间中抑制了多时相遥感影像间的光谱差异，突出了真实变化。

2、通过迭代方法，基于特征差值计算卡方距离，通过卡方距离获取各像素 点未变化的概率，并将概率作为像素点的权值。在迭代过程中，将权值加入到数 据标准化和慢特征分析计算中，逐步将较大的权值赋给未变化像素，给变化像素 赋0权值。这样就使得慢特征分析计算中更多地考虑未变化像素，排除变化像素 的干扰，提高变化像素和未变化像素的分离程度，从而提高变化检测精度。

3、具有适应度高、自组织、自学习的特点，变化检测精度较高，运算速度 快，效率高，适合多时相遥感影像的数据特点，适用于多时相遥感影像变化检测。

附图说明

图1为本发明第一种技术方案的具体流程图；

图2为本发明第二种技术方案的具体流程图。

具体实施方式

本发明关键发明点为将慢特征分析方法引入多时相遥感影像的变化检测中， 慢特征分析方法充分考虑多时相遥感影像间的空间对应关系，并将针对连续信号 分析的原始理论发展成基于离散数据集的变化检测算法。慢特征分析方法能提取 多时相遥感影像数据集中的不变特征作为特征空间，在此特征空间中，多时相遥 感影像间的光谱差异得到了抑制，真实变化得到了突出，从而可提高变化检测精 度。同时，慢特征分析方法实现简单、运算速度很快。

以下结合具体实施方式详细说明本发明第一种技术方案。

本具体实施方式采用MATLAB平台实现，MATLAB遥感影像读写函数为实 施基础。调用遥感影像读取函数，输入待读取遥感影像文件名，遥感影像就被读 入大小为N×P的矩阵中，矩阵中各元素为各波段对应的像素辐射值，其中，N 为遥感影像的波段数，P为遥感影像的像素数。调用两次遥感影像读写函数，分 别将两幅多时相遥感影像X和Y读入矩阵X和矩阵Y。MATLAB遥感影像读写 函数为本技术领域的公知技术，在此不作赘述。

基于矩阵X和矩阵Y对多时相遥感影像X和Y进行变化检测，下述步骤均 对矩阵X和矩阵Y进行操作：

（1）分别计算多时相遥感影像X和Y各波段的像素辐射值均值，计算公式 如下：

${\mu }_{x_j}=\frac{1}{P}\underset{i=1}{\overset{P}{\Sigma }}{x}_j^i---\left(1\right)$

式（1）中，是多时相遥感影像X第j个波段的像素辐射值均值；是 多时相遥感影像X第i个像素第j个波段的像素辐射值；P是多时相遥感影像X 的总像素数。

基于公式（1）分别计算多时相遥感影像Y各波段的像素辐射值均值

（2）分别计算多时相遥感影像X和Y各波段的像素辐射值标准差，计算公 式如下：

${\sigma }_{x_j}=\sqrt{\frac{1}{P}\underset{i=1}{\overset{P}{\Sigma }}{(x_j^i-{\mu }_{x_j})}^2}---\left(2\right)$

式（2）中，是多时相遥感影像X第j个波段的像素辐射值标准差；是多时相遥感影像X第j个波段的像素辐射值均值；是第i个像素第j个波段 的像素辐射值；P是多时相遥感影像X的总像素数。

基于公式（2）分别计算多时相遥感影像Y各波段的像素辐射值标准差

（3）根据步骤（1）获得的像素辐射值均值和步骤（2）获得的像素辐射值 标准差对多时相遥感影像X和Y进行标准化。

对多时相遥感影像各波段的像素辐射值进行标准化，具体为：

将各波段的像素辐射值分别减去该波段的像素辐射值均值，然后除以该波段 的像素辐射值标准差，得到各波段的标准化像素辐射值。

上述标准化过程见下式：

${\hat{x}}_j^i=\frac{x_j^i-{\mu }_{x_j}}{{\sigma }_{x_j}}---\left(3\right)$

式（3）中，是多时相遥感影像X第i个像素第j个波段的标准化像素辐 射值；是多时相遥感影像X第i个像素第j个波段的像素辐射值；是多时 相遥感影像X第j个波段的像素辐射值均值；是多时相遥感影像X第j个波 段的像素辐射值标准差。

基于公式（3）计算多时相遥感影像Y的标准化像素辐射值

（4）根据慢特征分析理论，基于标准化的多时相遥感影像X和Y计算多时 相遥感影像X和Y的差值影像的协方差矩阵A。

协方差矩阵A的计算公式如下；

$A=\frac{1}{P}\underset{i=1}{\overset{P}{\Sigma }}({\hat{x}}^i-{\hat{y}}^i){({\hat{x}}^i-{\hat{y}}^i)}^T---\left(4\right)$

式（4）中：

和分别是多时相遥感影像X和Y第i个像素的标准化像素辐射值向量， ${\hat{x}}^i=\left(\begin{array}{cccccc}{\hat{x}}_1^i& {\hat{x}}_2^i& ...& {\hat{x}}_j^i& ...& {\hat{x}}_N^i\end{array}\right),$${\hat{y}}^i=\left(\begin{array}{cccccc}{\hat{y}}_1^i& {\hat{y}}_2^i& ...& {\hat{y}}_j^i& ...& {\hat{y}}_N^i\end{array}\right),$是多时相遥感影像X第i个像素第j个波段的标准化像素辐射值，是多时相 遥感影像Y第i个像素第j个波段的标准化像素辐射值，N为波段数；

P是多时相遥感影像X的总像素数。

所得协方差矩阵A大小为N×N。

（5）根据慢特征分析理论，基于标准化的多时相遥感影像X和Y计算多时 相遥感影像X和Y的协方差矩阵之和，即为矩阵B。

矩阵B的计算公式如下：

$B=\frac{1}{2P}[\underset{i=1}{\overset{P}{\Sigma }}\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T+\underset{i=1}{\overset{P}{\Sigma }}\left({\hat{y}}^i\right){\left({\hat{y}}^i\right)}^T]---\left(5\right)$

式（5）中：

和分别是多时相遥感影像X和Y第i个像素的标准化像素辐射值向量， ${\hat{x}}^i=\left(\begin{array}{cccccc}{\hat{x}}_1^i& {\hat{x}}_2^i& ...& {\hat{x}}_j^i& ...& {\hat{x}}_N^i\end{array}\right),$${\hat{y}}^i=\left(\begin{array}{cccccc}{\hat{y}}_1^i& {\hat{y}}_2^i& ...& {\hat{y}}_j^i& ...& {\hat{y}}_N^i\end{array}\right),$是多时相遥感影像X第i个像素第j个波段的标准化像素辐射值，是多时相 遥感影像Y第i个像素第j个波段的标准化像素辐射值，N为波段数；

P是多时相遥感影像X的总像素数。

所得矩阵B大小为N×N。

（6）求解矩阵A相对矩阵B的广义特征值及广义特征值对应的特征向量， 并根据广义特征值大小对对应的特征向量排序，获得有序的特征向量矩阵。

采用MATLAB工具，调用特征值函数eig（）获得矩阵A相对矩阵B的广 义特征值矩阵和特征向量矩阵。计算公式为：

AW＝BWΛ  （6）

式（6）中，W是特征向量矩阵，Λ是广义特征值矩阵。

按照广义特征值大小从小到大排序，对广义特征值对应的特征向量进行相应 排序，得到有序的特征向量矩阵。排序后的广义特征值和波段一一对应，排序为 j的广义特征值为第j个波段对应的广义特征值。

（7）采用有序的特征向量矩阵，将标准化后的多时相遥感影像X和Y投影 到特征空间，并获得多时相遥感影像X和Y各波段的特征差值。

具体实施方式为：

将有序的特征向量矩阵分别乘以标准化后的多时相遥感影像X和Y，并将 相乘后所得矩阵各对应元素相减，得到多时相遥感影像X和Y的特征差值矩阵 SFA，即，特征差值影像。特征差值矩阵SFA是由各波段的特征差值SFA_j构成 的矩阵向量，特征差值SFA_j的计算如下：

${\mathrm{SFA}}_j=w_j\hat{x}-w_j\hat{y}---\left(7\right)$

式（7）中，SFA_j表示第j个波段的特征差值；w_j是第j个特征向量；和 分别为多时相遥感影像X和Y第j个波段的标准化辐射值向量矩阵，即，为 多时相遥感影像X各像素的标准化辐射向量构成的矩阵，为多时相遥感影 像Y各像素的标准化辐射向量构成的矩阵。

（8）根据特征差值和各波段对应的广义特征值，计算各像素点的卡方距离 T。

具体实施为：

分别计算各波段的像素点特征差值的平方并除以该波段对应的广义特征值 的开方，然后对所有波段求和得到像素点的卡方距离T，该距离满足卡方分布。

卡方距离T的计算公式如下：

$T^i=\underset{j=1}{\overset{N}{\Sigma }}\frac{{\left({\mathrm{SFA}}_j^i\right)}^2}{\sqrt{{\lambda }_j}}---\left(8\right)$

式（8）中，Tⁱ为第i个像素的卡方距离；为第i个像素第j个波段的 特征差值；λ_j为第j个波段的广义特征值；N为多时相遥感影像的波段数。

（9）基于像素点的卡方距离获得多时相遥感影像X和Y的变化检测结果。

判断像素点的卡方距离与阈值的大小，卡方距离不小于阈值的像素点为变化 像素点，赋值1；卡方距离小二阈值的像素点为未变化像素点，赋值0。基于各 像素点赋值分别分割多时相遥感影像X和Y并获得二值影像，保存二值影像作 为多时相遥感影像X和Y的变化检测结果图。

上述技术方案是根据慢特征分析原理从未变化数据集中提取不变特征，并认 为变化像素不会影响整个影像数据集的分布特性。如果影像区域内变化像素过 多，就会造成数据分布变化，从而影响慢特征分析变化检测效果。

为了克服慢特征分析方法存在的上述问题，获得更优的变化检测效果，本发 明的第二种技术方案将迭代定权方法引入上述多时相遥感影像的变化检测中。该 第二种方案根据卡方距离获得像素点为未变化像素点的概率，并将该概率作为该 像素点权值，并进行迭代。在迭代中，更多考虑未变化像素点，并逐步提出变化 像素点的影响。收敛时，变化像素点和未变化像素点可得到充分的分离。本发明 第二技术方案具体如下：

（1）初始化像素权值矩阵v，此时权值矩阵v为单位矩阵，权值矩阵v大小 为1×P，P为遥感影像的像素数。

采用单位权值矩阵v为多时相遥感影像X和Y中各像素点赋初始权值1。 权值矩阵v=[v¹，v²，….，vⁱ，…，v^P]，vⁱ表示影像第i个像素的权值，初始 化vⁱ=1，i＝1,2,...,P。

（2）根据当前权值矩阵，分别计算多时相遥感影像X和Y各波段的像素辐 射值的加权平均值，计算公式如下：

${\mu }_{x_j}=\frac{\underset{i=1}{\overset{P}{\Sigma }}\left(v^i{x}_j^i\right)}{\underset{i=1}{\overset{P}{\Sigma }}{v}^i}---\left(9\right)$

式（9）中，是多时相遥感影像X第j个波段的像素辐射值的加权平均 值；vⁱ是多时相遥感影像X第i个像素的当前权值；是多时相遥感影像X第 i个像素第j个波段的像素辐射值；P是多时相遥感影像X的总像素数。

基于公式（9）分别计算多时相遥感影像Y各波段的像素辐射值的加权平均 值

（3）分别计算多时相遥感影像X和Y各波段的像素辐射值的加权标准差， 计算公式如下：

${\sigma }_{x_j}=\sqrt{\frac{\underset{i=1}{\overset{P}{\Sigma }}\left[v^i{(x_j^i-{\mu }_{x_j})}^2\right]}{\underset{i=1}{\overset{P}{\Sigma }}{v}^i}}---\left(10\right)$

式（10）中，是多时相遥感影像X第j个波段的像素辐射值的加权标准 差；是多时相遥感影像X第j个波段的像素辐射值的加权平均值；是多时 相遥感影像X第i个像素第j个波段的像素辐射值；vⁱ是多时相遥感影像X第i 个像素的当前权值；P是多时相遥感影像X的总像素数。

基于公式（10）分别计算多时相遥感影像Y各波段的像素辐射值的加权标 准差

（4）根据像素辐射值的加权平均值和像素辐射值的加权标准差对多时相遥 感影像X和Y进行标准化。

标准化方法同第一种技术方案的步骤（3），将各波段的像素辐射值分别减去 该波段的像素辐射值的加权平均值，然后除以该波段的像素辐射值的加权标准 差，得到各波段的标准化像素辐射值分别表示多时相遥感影像 X和Y第j个波段第i个像素的标准化辐射值。

（5）根据慢特征分析理论，基于标准化的多时相遥感影像X和Y计算多时 相遥感影像X和Y的差值影像的协方差矩阵A。

和第一种技术方案中计算多时相遥感影像X和Y的差值影像的协方差矩阵 A不同，本方案在计算矩阵A时考虑了像素权值，计算公式如下：

$A=\frac{\underset{i=1}{\overset{P}{\Sigma }}\left[v^i({\hat{x}}^i-{\hat{y}}^i){({\hat{x}}^i-{\hat{y}}^i)}^T\right]}{\underset{i=1}{\overset{P}{\Sigma }}{v}^i}---\left(11\right)$

式（11）中，和分别是多时相遥感影像X和Y第i个像素的标准化像 素辐射值向量；vⁱ是多时相遥感影像第i个像素的当前权值；P是多时相遥感影 像X的总像素数。

（6）根据慢特征分析理论，基于标准化的多时相遥感影像X和Y计算多时 相遥感影像X和Y的协方差矩阵之和，即为矩阵B。

和第一种技术方案中计算多时相遥感影像X和Y的协方差矩阵之和不同， 本方案在计算矩阵B时考虑了像素权值，计算公式如下：

$B=\frac{1}{2\underset{i=1}{\overset{P}{\Sigma }}{v}^i}[\underset{i=1}{\overset{P}{\Sigma }}\left[v^i\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T\right]+\underset{i=1}{\overset{P}{\Sigma }}\left[v^i\left({\hat{y}}^i\right){\left({\hat{y}}^i\right)}^T\right]]---\left(12\right)$

式（12）中，和分别是多时相遥感影像X和Y第i个像素的标准化像 素辐射值向量；vⁱ是多时相遥感影像第i个像素的当前权值；P是多时相遥感影 像X的总像素数。

（7）求解矩阵A相对矩阵B的广义特征值及广义特征值对应的特征向量， 并根据广义特征值大小对对应的特征向量进行排序，获得有序的特征向量矩阵。

本步骤的具体实施方式参见第一种技术方案的步骤（6），在此不作赘述。

（8）采用有序的特征向量矩阵，将标准化后的多时相遥感影像X和Y投影 到特征空间，并获得多时相遥感影像X和Y各波段的特征差值。

本步骤的具体实施方式参见第一种技术方案的步骤（7），在此不作赘述。

（9）根据各波段的特征差值和广义特征值，计算各像素点的卡方距离T， 并基于像素点的卡方距离和卡方分布更新像素点的权值，从而获得更新后的像素 权值矩阵，以更新后的权值矩阵为当前权值矩阵。

遍历多时相遥感影像X和Y，计算以波段数N为参数的卡方分布中，各像 素点大于卡方距离的概率，并将概率作为权值赋给对应的像素。

各像素点权值的计算公式如下：

vⁱ＝P{χ²(N)＞Tⁱ}  （13）

式（13）中，vⁱ表示第i个像素的当前权值；χ²(N)表示以波段数N为参 数的卡方分布；P{χ²(N)＞Tⁱ}表示以波段数N为参数的卡方分布中、第i个像 素大于其卡方距离Tⁱ的概率。

（10）判断本次迭代和上次迭代各波段的特征差值的最大差值，如果最大差 值大于设定阈值或者本次迭代为首次迭代，重新执行步骤（2）～（9）；如果最大 差值小于设定阈值，执行步骤（11）。

设定阈值用于迭代收敛，可根据经验进行设定，本具体实施中，该阈值设定 为10^-6。

上述各波段的特征差值的最大差值为：

针对各波段，分别求取各波段在本次迭代和上次迭代中对应的特征差值的差 值，所获的差值的最大值即为各波段的特征差值的最大差值。

（11）基于像素点的卡方距离获得多时相遥感影像X和Y的变化检测结果。 具体可参见第一技术方案的步骤（9）。

以下通过对比试验来验证本发明的有益效果。

本试验采用的数据为Landsat ETM+数据，分辨率30m，共6个波段，影像 尺寸400像素×400像素。分别采用经典的变化向量分析法（方法1）、主成分分 析法（方法2）、多元变化检测法（方法3）和本发明方法进行变化检测。

变化检测评价指标：采用定量评价方法，参考样本通过目视解译标注，总共 选择了3724像素的变化样本和10888像素的未变化样本。评价指标采用如下两 个指标：

1)Kappa系数：

Kappa系数是用来评价分类问题的权威评价指标。Kappa系数越大，精度越 高。在变化检测中，可以将变化检测结果看做二分类问题（变化和未变化）。本 试验中，选取的是方法1～3所能获得的最高Kappa系数来评价方法1～3的检测 能力。

Kappa系数计算方法为：

根据样本获取混淆矩阵，见表1。

表1 混淆矩阵

表1中，Nnn表示未变化样本被检测成为未变化的数量；Nnc表示未变化样 本被检测成为变化的数量；Ncn表示变化样本被检测成为未变化的数量；Ncc表 示变化样本被检测成为变化的数量；Nnp为Nnc和Nnc之和，Ncp为Ncn和Ncc 之和，Npn为Nnn和Ncp之和，Npc为Nnc和Ncc之和，N为总样本数。

Kappa系数的计算公式为：

$\mathrm{Kappa}=\frac{(\mathrm{Nnp}\times \mathrm{Npn}+\mathrm{Ncp}\times \mathrm{Npc})/N^2-(\mathrm{Nnn}+\mathrm{Ncc})/N}{1-(\mathrm{Nnn}+\mathrm{Ncc})/N}---\left(10\right)$

2)整体精度：

整体精度（Overall Accuracy，OA）是用来评价分类问题的权威评价指标。 整体精度越高，检测精度越高。OA的计算方法同样基于表1所示的混淆矩阵， 整体精度OA计算公式为：

$\mathrm{OA}=\frac{\mathrm{Nnn}+\mathrm{Ncc}}{N}---\left(11\right)$

采用Kappa系数和整体精度评价方法1～3和本发明方法的变化检测能力， 评价指标见表2。

表2 对比试验结果

  本发明方法 方法1 方法2 方法3 整体精度 0.9885 0.9241 0.9245 0.9535 Kappa系数 0.9699 0.7859 0.7853 0.8944

从表2可见，本发明方法能获得更高的整体精度和Kappa系数，表明本发 明方法具有更强的变化检测能力。本发明方法和方法1～3的Kappa系数的差异， 比整体精度差异更加明显，因为整体精度未考虑变化和未变化样本的数量问题， 而Kappa系数考虑到了该问题，因此，Kappa系数的评价更加客观。

由此可得出结论，与传统变化检测方法相比，本发明方法拥有更高的变化检 测精度。本发明充分考虑了两幅遥感影像之间的联系，突出变化区域，抑制未变 化区域，提高变化检测精度。