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

Abstract

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

Description

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

技术领域

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

背景技术

人类活动对地球表面环境造成了巨大影响，这种影响体现在环境变化、城市发展等各个方面。因此，实时精确地获取地球地表覆盖的变化情况对于环境监测和资源管理意义重大（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 usingremotely-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 inMultispectral,Bitemporal Image Data:New Approaches to Change DetectionStudies.″Remote Sensing of Environment 64(1):1-19.），该算法基于数学统计学中的典型相关分析理论，针对两组数据之间的相关性进行分析，从而对两个数据集分别进行变换，得到变换特征。多元变化检测算法得到了广泛的研究，但是该算法没有考虑到同源遥感影像之间波段的对应关系，无法充分利用遥感影像的光谱信息。

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

本文涉及如下参考文献：

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

[2] A.A.Nielsen,″The Regularized Iteratively Reweighted MAD Method forChange 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 changedetection 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 resourcemanagers:Understanding concepts and tradeoffs in the design of landscapemonitoring projects,″Remote Sens.Environ.,vol.113,pp.1382-1396,Jul 2009.

[5] T.Celik,″Unsupervised Change Detection in Satellite Images UsingPrincipal 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 UnsupervisedChange Detection Based on Change Vector Analysis in the Polar Domain,″IEEETrans.Geosci.Remote Sens.,vol.45,pp.218-236,Jan 2007.

[7] P.Coppin,I.Jonckheere,K.Nackaerts,B.Muys,and E.Lambin,″Digitalchange 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}{\mathrm{\Σ}}}({\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}{\mathrm{\Σ}}}\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T+\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\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}{\mathrm{\Σ}}}\left[v^i({\hat{x}}^i-{\hat{y}}^i){({\hat{x}}^i-{\hat{y}}^i)}^T\right]}{\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{v}^i}$

其中，

和

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

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

$B=\frac{1}{2\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{v}^i}[\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\left[v^i\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T\right]+\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\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各波段的像素辐射值均值，计算公式如下：

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

式（1）中，

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

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

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

${\mathrm{\σ}}_{x_j}=\sqrt{\frac{1}{P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{(x_j^i-{\mathrm{\μ}}_{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-{\mathrm{\μ}}_{x_j}}{{\mathrm{\σ}}_{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}{\mathrm{\Σ}}}({\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}{\mathrm{\Σ}}}\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T+\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\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}{\mathrm{\Σ}}}\frac{{\left({\mathrm{SFA}}_j^i\right)}^2}{\sqrt{{\mathrm{\λ}}_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各波段的像素辐射值的加权平均值，计算公式如下：

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

式（9）中，

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

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

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

${\mathrm{\σ}}_{x_j}=\sqrt{\frac{\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\left[v^i{(x_j^i-{\mathrm{\μ}}_{x_j})}^2\right]}{\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{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}{\mathrm{\Σ}}}\left[v^i({\hat{x}}^i-{\hat{y}}^i){({\hat{x}}^i-{\hat{y}}^i)}^T\right]}{\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{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}{\mathrm{\Σ}}}{v}^i}[\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\left[v^i\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T\right]+\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\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}\×\mathrm{Npn}+\mathrm{Ncp}\×\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系数的评价更加客观。

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

Claims (9)

1.基于慢特征分析的遥感影像变化检测方法，其特征在于：

将多时相遥感影像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.如权利要求1所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

步骤2具体为：

采用公式

标准化多时相遥感影像，其中，

是多时相遥感影像第i个像素第j个波段的标准化像素辐射值；是多时相遥感影像第i个像素第j个波段的像素辐射值；

是多时相遥感影像第j个波段的像素辐射值均值；

是多时相遥感影像第j个波段的像素辐射值标准差。

3.如权利要求1所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

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

$A=\frac{1}{P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}({\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}{\mathrm{\Σ}}}\left({\hat{x}}^i\right){\left({\hat{x}}^i\right)}^T+\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\left({\hat{y}}^i\right){\left({\hat{y}}^i\right)}^T]$

其中，

和

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

4.如权利要求1所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

步骤5具体为：

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

5.如权利要求1所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

所述的像素点的卡方距离

其中，Tⁱ为第i个像素的卡方距离；

为第i个像素第j个波段的特征差值；λ_j为第j个波段的广义特征值；N为多时相遥感影像的波段数。

6.如权利要求1所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

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

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

7.基于慢特征分析的遥感影像变化检测方法，其特征在于，包括步骤：

步骤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的变化检测结果。

8.如权利要求6所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

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

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

其中，

和

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

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

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

其中，

和

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

9.如权利要求6所述的基于慢特征分析的遥感影像变化检测方法，其特征在于：

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

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

Images