CN103903258A - 基于次序统计量谱聚类的遥感图像变化检测方法 - Google Patents

Abstract

本发明公开了一种基于次序统计量谱聚类的遥感图像变化检测方法，它属于图像处理技术领域，主要解决现有技术检测效果不理想，检测精度低、适用范围窄的问题。包括图像读入、差异图处理、样本选择、计算矩阵、估算特征、选取矩阵、聚类和图像输出等步骤，提高了变化检测的精度，减少了错分率，对不同类型的遥感图像均可获得较好的效果，可应用于环境监测、灾害评估等领域。

Description

基于次序统计量谱聚类的遥感图像变化检测方法

技术领域

本发明属于图像处理技术领域，涉及遥感图像的变化检测，特别涉及一种基于次序统计量谱聚类的遥感图像变化检测方法，该方法可用于植被覆盖、土地利用、自然灾害评估和城区规划等领域。

背景技术

变化检测是遥感技术的重要研究方向之一，它的目标是通过分析在同一地区不同时段获得的两幅遥感图像来辨识其变化信息。它在土地利用、火灾面积调查、开展资源调查、环境检测、军事侦察等方面有广泛的应用领域。变化检测方法大体可以分为两类：有监督法和无监督法。但由于有监督法需要地面真实数据来获取变化区域的训练样区，从而进行变化检测，然而想要得到地面的真实信息是比较困难的，所以无监督变化检测方法是常用的变化检测方法。

无监督的遥感图像变化检测过程通常可以分为三个步骤：(1)预处理；(2)图像比较；(3)差异图分析。

传统的变化检测算法是在差异图上进行阈值判断来提取变化信息，但传统构造差异图方法存在许多缺陷，例如没有考虑图像像素的邻域信息，简单计算可能出现虚警率，变化阈值难以确定等问题。另外传统的聚类算法，如Kmeans、EM算法等都是建立在凸球形的样本空间上，当样本空间非凸时，算法易陷入局部最优。

谱聚类算法是一种比较流行的聚类算法，它已广泛应用于计算机视觉、语音识别、文本挖掘等领域，该算法可以在任意形状的样本空间上聚类且收敛于全局最优。与传统聚类算法相比，谱聚类算法的优点在于：(1)谱聚类是直接对拉普拉斯矩阵的特征向量求解，然后进行聚类划分，所以谱聚类可以对样本空间是非凸的数据进行聚类求解；(2)谱聚类仅与数据点的数目有关，而与维数无关，因而可以避免求解高维特征向量时造成奇异性问题；(3)谱聚类是根据样本间的相似关系建立矩阵，通过计算特征向量找出样本间的内在联系，实质上谱聚类是通过把低维数据映射到高维空间，然后在高维空间进行聚类，所以谱聚类的性能通常优于一般的聚类方法。

虽然谱聚类有很多优点，但是当聚类样本数目比较大时，谱聚类在计算和存储相似度矩阵方面是十分困难的。比如在n个样本之间构造相似度矩阵需要内存是O(n²)，计算特征向量分解需要的时间复杂度为O(n³)，内存使用复杂度是O(n²)。Fowlkes等人提出使用Nystrom逼近来避免计算整个相似度矩阵，此算法很好的降低了算法时间和空间的复杂度。基于Nystrom的谱聚类算法可简述为：在样本中随机选取少量样本作为采样样本，计算采样样本之间的相似度矩阵A，并计算采样样本与剩余样本之间的相似度矩阵B；再使用Nystrom方法估计全体样本之间的相似度矩阵W及其主特征向量；并从这些主特征向量中选取k个特征向量；最后采用K-means等传统聚类算法完成最终的聚类划分。但由于在构造相似矩阵时仅简单使用了样本间的欧式距离来衡量它们之间的相似性，不能很好的抑制噪声点的影响，从而会影响变化区域的检测，造成检测精度不够高的问题。

发明内容

本发明的目的在于提出了一种基于次序统计量谱聚类的遥感图像变化检测方法，针对现有技术检测效果不理想，适用范围窄等问题，能更好的检测出图像的变化区域。

为实现上述目的，本发明的技术方案是：基于次序统计量谱聚类的遥感图像变化检测方法，包括如下步骤：

（1）输入在两个时相获取的同一地区的两幅大小为M×N的遥感图像T1和T2，其中M、N为自然数，代表图像的行数和列数，分别表示图像的长和宽；

（2）利用对数比值算子得到差异矩阵C，并对差异矩阵C进行归一化处理得到差异图DI；

（3）样本选择：假设差异图DI有n个像素点，从中随机采样m个像素点作为样本点，则剩余为n-m个像素点；

（4）利用基于次序统计量的像素灰度相似度函数来分别计算采样点间相似矩阵A，采样点与剩余点间相似矩阵B；

（5）使用Nystrom估计计算特征值Λ_R和相应的特征向量U_R；

（6）从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，并对矩阵V进行行归一化处理，得到归一化矩阵

（7）将

的每行视为一个点，用K-means算法将它们聚为k类；

（8）若

的第i行数据属于第j类，则将对应像素点也归为第j类，再将结果转化为最终的变化检测结果。

上述步骤(2)所述的利用对数比值算子得到差异图，并对差异图进行归一化处理，按照如下步骤进行：

a)对图像T1在坐标(i,j)处的灰度值T1(i,j)与图像T2在坐标(i,j)处的灰度值T2(i,j)进行对数比值运算，得到差异矩阵C在坐标(i,j)处的值C(i,j)：

进而得到差异矩阵C＝{C(i,j)}；

其中log表示自然对数,0＜i≤M,0＜j≤N；

b）归一化差异矩阵，得到差异图DI在坐标(i,j)处的值DI(i,j)：

$\mathrm{DI}(i,j)=\frac{C(i,j)-\min \left(C\right)}{\max \left(C\right)-\min \left(C\right)},$ 进而得到差异图：DI＝{DI(i,j)}。

其中min为求最小值函数，max为求最大值函数，C为差异矩阵，0＜i≤M，0＜j≤N；

上述步骤(4)所述的利用基于次序统计量的像素灰度相似度函数来分别计算采样点间相似矩阵A，采样点与剩余点间相似矩阵B，按照如下公式进行：

$S_{\mathrm{ij}}=\exp (-\frac{{\left|\right|x\left(i\right)-x\left(j\right)\left|\right|}^2}{{2\mathrm{\σ}}^2})$

其中，x(i)和x(j)分别为像素i，j灰度的n次序统计量，即通过将像素i和j各自的1环邻域像素灰度近似的视为i，j灰度的n次观测值，且将其分别由小到大排序得到，σ为尺度参数，一般通过经验取得，exp为指数函数，S_ij为像素i，j的相似度，在计算相似矩阵A时，0＜i≤m，0＜j≤m，在计算相似矩阵B时，0＜i≤m，0＜j≤n-m，从而得到相似矩阵A∈R^m×m和相似矩阵B∈R^m×(n-m),R为实数集。

上述步骤(5)所述的使用Nystrom估计计算特征值和相应的特征向量，按照如下步骤进行：

a)对角矩阵D按如下公式计算，

$D=\mathrm{diag}(\left[\begin{array}{c}{\mathrm{Al}}_m+{B1}_{n-m}\\ B^T{1}_m+B^T{A}^{-1}{B1}_{n-m}\end{array}\right])$

其中1代表元素为1的列向量，diag为求对角矩阵函数，A为采样点间相似矩阵，B为采样点与剩余点间相似矩阵，B^T为矩阵B的转置矩阵，A^-1为矩阵A的逆矩阵；

b)分别计算相似矩阵A和B的共轭矩阵

和

$\stackrel{\‾}{A}=D^{-\frac{1}{2}}{\mathrm{AD}}^{-\frac{1}{2}},\stackrel{\‾}{B}=D^{-\frac{1}{2}}{\mathrm{BD}}^{-\frac{1}{2}}$

其中A为采样点间相似矩阵，B为采样点与剩余点间相似矩阵，D为对角矩阵。

c)构造矩阵 $R=\stackrel{\‾}{A}+{\stackrel{\‾}{A}}^{-\frac{1}{2}}{\stackrel{\‾}{\mathrm{BB}}}^T{\stackrel{\‾}{A}}^{-\frac{1}{2}};$

其中

为相似矩阵A的共轭矩阵，为相似矩阵B的共轭矩阵，为矩阵

的转置矩阵。

d)对矩阵R进行特征分解，

其中Λ_R为矩阵R的特征值，U_R为矩阵R的特征向量,

为U_R的转置矩阵。

上述步骤(6)所述的从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，并对矩阵V进行行归一化,得到归一化矩阵

按照如下步骤进行：

a)从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，按照如下公式进行：

$V=\left[\frac{\stackrel{\‾}{A}}{{\stackrel{\‾}{B}}^T}\right]{\stackrel{\‾}{A}}^{-\frac{1}{2}}{\left(U_R\right)}_{:,1:K}{\left({{\mathrm{\Λ}}_R}^{-\frac{1}{2}}\right)}_{1:K,1:K}$

其中

为相似矩阵A的共轭矩阵，

为相似矩阵B的共轭矩阵，为矩阵

的转置矩阵,Λ_R为矩阵R的特征值，U_R为矩阵R的特征向量，K为选取的特征值个数。

b)对矩阵V进行行归一化处理，得到归一化矩阵

${\stackrel{\‾}{V}}_{\mathrm{ir}}=\frac{V_{\mathrm{ir}}}{\sqrt{\underset{r=1}{\overset{K}{\mathrm{\Σ}}}{V_{\mathrm{ir}}}^2}},i=1,...,n,r=1,...K$

其中n为像素个数，i和r为矩阵V中元素的坐标，K为选取的特征值个数。

本发明与现有的技术相比具有以下优点：

本发明采用基于次序统计量的像素灰度相似度函数来计算相似矩阵，很好的利用了样本的邻域信息，有效的抑制了噪声点的影响，提高了分类精确率，具有很好的鲁棒性。

附图说明

图1是本发明的流程图；

图2是Bern地区SAR图像数据集和该数据集的标准变化检测图；

图3是用本发明和对比方法对Bern地区SAR图像数据集的变化检测结果图；

图4是Ottawa地区SAR图像数据集和该数据集的标准变化检测图；

图5是用本发明和对比方法对Ottawa地区SAR图像数据集的变化检测结果图。

具体实施方式

参照图1，本发明的具体实施过程如下：

步骤1：输入在两个时相获取的同一地区的两幅大小为M×N的遥感图像T1和T2，其中M、N为自然数，代表图像的行数和列数，分别表示图像的长和宽。

步骤2：利用对数比值算子得到差异矩阵C，并对差异矩阵C进行归一化处理得到差异图DI。

2a)对图像T1在坐标(i,j)处的灰度值T1(i,j)与图像T2在坐标(i,j)处的灰度值T2(i,j)进行对数比值运算，得到差异矩阵C在坐标(i,j)处的值C(i,j)：

进而得到差异矩阵C＝{C(i,j)}；

其中log表示自然对数,0＜i≤M,0＜j≤N；

2b)归一化差异矩阵，得到差异图DI在坐标(i,j)处的值DI(i,j)：

$\mathrm{DI}(i,j)=\frac{C(i,j)-\min \left(C\right)}{\max \left(C\right)-\min \left(C\right)},$ 进而得到差异图：DI＝{DI(i,j)}；

其中min为求最小值函数，max为求最大值函数，C为差异矩阵，0＜i≤M，0＜j≤N。

步骤3：样本选择：假设差异图DI有n个像素点，从中随机采样m个像素点作为样本点，则剩余为n-m个像素点；

3a)将差异图DI中的像素点排为一列，即为x₁,...,x_n；

3b)从中随机采样m个样本点，即为x₁,...,x_m；

3c)则剩余为n-m个像素点，即为x_m+1,...,x_n。

步骤4：利用基于次序统计量的像素灰度相似度函数来分别计算采样点间相似矩阵A，采样点与剩余点间相似矩阵B；

按照如下公式进行：

$S_{\mathrm{ij}}=\exp (-\frac{{\left|\right|x\left(i\right)-x\left(j\right)\left|\right|}^2}{{2\mathrm{\σ}}^2})$

其中，x(i)和x(j)分别为像素i,j灰度的n次序统计量，即通过将像素i和j各自的1环邻域像素灰度近似的视为i,j灰度的n次观测值，且将其分别由小到大排序得到，σ为尺度参数，一般通过经验取得，exp为指数函数，S_ij为像素i,j的相似度，在计算相似矩阵A时，0＜i≤m，0＜j≤m，在计算相似矩阵B时，0＜i≤m，0＜j≤n-m，从而得到相似矩阵A∈R^m×m和相似矩阵B∈R^m×(n-m),R为实数集。

步骤5：使用Nystrom估计计算特征值Λ_R和相应的特征向量U_R；

5a)对角矩阵D按如下公式计算，

$D=\mathrm{diag}(\left[\begin{array}{c}{\mathrm{Al}}_m+{B1}_{n-m}\\ B^T{1}_m+B^T{A}^{-1}{B1}_{n-m}\end{array}\right])$

其中1代表元素为1的列向量，diag为求对角矩阵函数，A为采样点间相似矩阵，B为采样点与剩余点间相似矩阵，B^T为矩阵B的转置矩阵，A^-1为矩阵A的逆矩阵；

5b)分别计算相似矩阵A和B的共轭矩阵和

$\stackrel{\‾}{A}=D^{-\frac{1}{2}}{\mathrm{AD}}^{-\frac{1}{2}},\stackrel{\‾}{B}=D^{-\frac{1}{2}}{\mathrm{BD}}^{-\frac{1}{2}}$

其中A为采样点间相似矩阵，B为采样点与剩余点间相似矩阵，D为对角矩阵。

5c)构造矩阵 $R=\stackrel{\‾}{A}+{\stackrel{\‾}{A}}^{-\frac{1}{2}}{\stackrel{\‾}{\mathrm{BB}}}^T{\stackrel{\‾}{A}}^{-\frac{1}{2}};$

其中

为相似矩阵A的共轭矩阵，为相似矩阵B的共轭矩阵，

为矩阵

的转置矩阵。

5d)对矩阵R进行特征分解，

其中Λ_R为矩阵R的特征值，U_R为矩阵R的特征向量,

为U_R的转置矩阵。

步骤6：从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，并对矩阵V进行行归一化处理，得到归一化矩阵

6a)从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，按照如下公式进行：

$V=\left[\frac{\stackrel{\‾}{A}}{{\stackrel{\‾}{B}}^T}\right]{\stackrel{\‾}{A}}^{-\frac{1}{2}}{\left(U_R\right)}_{:,1:K}{\left({{\mathrm{\Λ}}_R}^{-\frac{1}{2}}\right)}_{1:K,1:K}$

其中为相似矩阵A的共轭矩阵，为相似矩阵B的共轭矩阵，

为矩阵

的转置矩阵,Λ_R为矩阵R的特征值，U_R为矩阵R的特征向量，K为选取的特征值个数。

6b)对矩阵V进行行归一化处理，得到归一化矩阵

${\stackrel{\‾}{V}}_{\mathrm{ir}}=\frac{V_{\mathrm{ir}}}{\sqrt{\underset{r=1}{\overset{K}{\mathrm{\Σ}}}{V_{\mathrm{ir}}}^2}},i=1,...,n,r=1,...K$

其中n为像素个数，i和r为矩阵V中元素的坐标，K为选取的特征值个数。

步骤7：将

的每行视为一个点，用K-means算法将它们聚为k类；

步骤8：若

的第i行数据属于第j类，则将对应像素点也归为第j类，再将结果转化为最终的变化检测结果；

8a)给标签为1的点赋予灰度值0，给标签为2的点赋予灰度值255。则得到变化检测结果图；

8b)输出变化检测结果图。

本发明的效果可通过以下仿真实验进一步说明：

1.实验条件：

仿真实验环境为：Windows7,Dual-Core AMD Opteron(tm)Processor2220,基本频率2.8GHz,软件平台为MATLAB7.8.0(R2009a)。

第一个数据集是Bern地区SAR图像数据集，数据集的原始图像都是通过ERS-2获得的瑞士Bern地区的前后两时相SAR图像组成，前一时相的SAR图像是在洪水灾害刚刚发生后获得的，图像中阴暗部分为受洪水影响的区域，如图2(a)所示，后一时相的SAR图像是洪水几乎完全消失的时获得，如图2(b)所示。两幅图像的尺寸均为301×301像素，灰度级为256，图2(c)为该数据集的标准变化检测图。

第二个数据集是Ottawa地区SAR图像数据集，它由两幅在不同时刻拍摄的Radarsat SAR图像组成，变化区域主要是由洪水灾情引起的，如图4(a)和4(b)所示，其中图4(a)为1997年5月图像，图4(b)为1997年8月的图像，图像大小均为290×350像素，灰度级为256。图4(c)为该数据集的标准变化检测图。

2.实验评价指标

实验使用的评价指标是漏检数，错检数、总错误数和正确检测率。其中，漏检数为没有检测出来的实际发生变化的像素的总和，错检数为实际没有发生变化但被当作变化检测出来的像素的总和，总错误数是漏检数和错检数之和，正确检测率为像素总数减去总错误数，然后再除以像素总数。我们一般认为漏检数、错检数和总错误数越小，并且正确检测率越大，则变化检测效果越好。

3.实验内容和实验结果

用本发明方法和现有的3种变化检测方法对2个图像数据集进行变化检测，该3种对比方法分别为：采用对数比值算子构造差异图然后用K-means聚类的方法记作Kmeans法；采用对数比值算子构造差异图然后用基于Nystrom逼近的谱聚类方法记做SC法；采用对数比值算子构造差异图然后用Krinidis在文章“Arobust fuzzy local information C-means clustering algorithm”中提出的方法聚类，记作FLICM法。

实验1.

用本发明方法和现有的3种变化检测方法对第一个数据集进行变化检测，参数设置为，核参数σ为19，模糊因子m为2，实验结果如图3和表1所示，其中3(a)是Kmeans法的变化检测结果图，3(b)是SC法的变化检测结果图，3(c)是FLICM法的变化检测结果图，3(d)是本发明的变化检测结果图。

从图3可以看出，3(a)和3(b)两幅图中存在较多的伪变化区域，造成过多的误检数，从而使得总错误数过大；图3(c)虽然直观效果看起来很好，但从表1的数值统计结果可以看出，漏检数很大，从而总错误数也较大，本发明的方法与FLICM法相比具有更少的错误数，且能获得很好的直观效果。

实验2.

用本发明方法和现有的3种变化检测方法对第二个数据集进行变化检测，参数设置为，核参数σ为34，模糊因子m为2，实验结果如图5和表1所示，其中5(a)是Kmeans法的变化检测结果图，5(b)是SC法的变化检测结果图，5(c)是FLICM法的变化检测结果图，5(d)是本发明的变化检测结果图。

从图5可以看出，图5(a)和5(b)中伪变化信息较多，误检严重，且噪声影响较大；图5(c)则漏检严重。本发明方法在降低噪声的同时有效地提取出了变化信息。

用本发明方法和现有的3种变化检测方法对2个数据集进行变化检测的数值统计结果如表1所示。

表1.变化检测的数值统计结果

由表1可以看出，Kmeans法的错检数较多，FLICM法的漏检数较多，而传统的谱聚类SC法虽有所中和，但仍存在较多的总错误数；本发明通过利用邻域信息来构造相似矩阵，很好的抑制了噪声点的影响，大幅度减少了总错误数，提高了分类精确率。相比于其它算法的聚类结果，本发明更准确的检测出图像的变化区域。

Claims (5)

1.基于次序统计量谱聚类的遥感图像变化检测方法，其特征在于：包括如下步骤：

（1）输入在两个时相获取的同一地区的两幅大小为M×N的遥感图像T1和T2，其中M、N为自然数，代表图像的行数和列数，分别表示图像的长和宽；

（2）利用对数比值算子得到差异矩阵C，并对差异矩阵C进行归一化处理得到差异图DI；

（3）样本选择：假设差异图DI有n个像素点，从中随机采样m个像素点作为样本点，则剩余为n-m个像素点；

（4）利用基于次序统计量的像素灰度相似度函数来分别计算采样点间相似矩阵A，采样点与剩余点间相似矩阵B；

（5）使用Nystrom估计计算特征值Λ_R和相应的特征向量U_R；

（6）从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，并对矩阵V进行行归一化处理，得到归一化矩阵

（7）将

的每行视为一个点，用K-means算法将它们聚为k类；

（8）若

的第i行数据属于第j类，则将对应像素点也归为第j类，再将结果转化为最终的变化检测结果。

2.根据权利要求1所述的基于次序统计量谱聚类的遥感图像变化检测方法，其特征在于：其中步骤(2)所述的利用对数比值算子得到差异图，并对差异图进行归一化处理，按照如下步骤进行：

a)对图像T1在坐标(i,j)处的灰度值T1(i,j)与图像T2在坐标(i,j)处的灰度值T2(i,j)进行对数比值运算，得到差异矩阵C在坐标(i,j)处的值C(i,j)：

进而得到差异矩阵C＝{C(i,j)}；

其中log表示自然对数,0＜i≤M,0＜j≤N；

b）归一化差异矩阵，得到差异图DI在坐标(i,j)处的值DI(i,j)：

$\mathrm{DI}(i,j)=\frac{C(i,j)-\min \left(C\right)}{\max \left(C\right)-\min \left(C\right)},$ 进而得到差异图：DI＝{DI(i,j)}；

其中min为求最小值函数，max为求最大值函数，C为差异矩阵，0＜i≤M，0＜j≤N。

3.根据权利要求1所述的基于次序统计量谱聚类的遥感图像变化检测方法，其特征在于：其中步骤(4)所述的利用基于次序统计量的像素灰度相似度函数来分别计算采样点间相似矩阵A，采样点与剩余点间相似矩阵B，按照如下公式进行：

$S_{\mathrm{ij}}=\exp (-\frac{{\left|\right|x\left(i\right)-x\left(j\right)\left|\right|}^2}{{2\mathrm{\σ}}^2})$

其中，x(i)和x(j)分别为像素i，j灰度的n次序统计量，即通过将像素i和j各自的1环邻域像素灰度近似的视为i，j灰度的n次观测值，且将其分别由小到大排序得到，σ为尺度参数，一般通过经验取得，exp为指数函数，S_ij为像素i，j的相似度，在计算相似矩阵A时，0＜i≤m，0＜j≤m，在计算相似矩阵B时，0＜i≤m，0＜j≤n-m，从而得到相似矩阵A∈R^m×m和相似矩阵B∈R^m×(n-m),R为实数集。

4.根据权利要求1所述的基于次序统计量谱聚类的遥感图像变化检测方法，其特征在于：其中步骤(5)所述的使用Nystrom估计计算特征值和相应的特征向量，按照如下步骤进行：

a)对角矩阵D按如下公式计算，

$D=\mathrm{diag}(\left[\begin{array}{c}{\mathrm{Al}}_m+{B1}_{n-m}\\ B^T{1}_m+B^T{A}^{-1}{B1}_{n-m}\end{array}\right])$

其中1代表元素为1的列向量，diag为求对角矩阵函数，A为采样点间相似矩阵，B为采样点与剩余点间相似矩阵，B^T为矩阵B的转置矩阵，A^-1为矩阵A的逆矩阵；

b)分别计算相似矩阵A和B的共轭矩阵

和

$\stackrel{\‾}{A}=D^{-\frac{1}{2}}{\mathrm{AD}}^{-\frac{1}{2}},\stackrel{\‾}{B}=D^{-\frac{1}{2}}{\mathrm{BD}}^{-\frac{1}{2}}$

其中A为采样点间相似矩阵，B为采样点与剩余点间相似矩阵，D为对角矩阵；

c)构造矩阵 $R=\stackrel{\‾}{A}+{\stackrel{\‾}{A}}^{-\frac{1}{2}}{\stackrel{\‾}{\mathrm{BB}}}^T{\stackrel{\‾}{A}}^{-\frac{1}{2}};$

其中

为相似矩阵A的共轭矩阵，

为相似矩阵B的共轭矩阵，

为矩阵

的转置矩阵；

d)对矩阵R进行特征分解，

其中Λ_R为矩阵R的特征值，U_R为矩阵R的特征向量,

为U_R的转置矩阵。

5.根据权利要求1所述的基于次序统计量谱聚类的遥感图像变化检测方法，其特征在于：其中步骤(6)所述的从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，并对矩阵V进行行归一化,得到归一化矩阵

按照如下步骤进行：

a)从特征值Λ_R中选取前K个最大的特征值所对应的特征向量组成矩阵V，按照如下公式进行：

$V=\left[\frac{\stackrel{\‾}{A}}{{\stackrel{\‾}{B}}^T}\right]{\stackrel{\‾}{A}}^{-\frac{1}{2}}{\left(U_R\right)}_{:,1:K}{\left({{\mathrm{\Λ}}_R}^{-\frac{1}{2}}\right)}_{1:K,1:K}$

其中为相似矩阵A的共轭矩阵，

为相似矩阵B的共轭矩阵，

为矩阵

的转置矩阵,Λ_R为矩阵R的特征值，U_R为矩阵R的特征向量，K为选取的特征值个数；

b)对矩阵V进行行归一化处理，得到归一化矩阵

${\stackrel{\‾}{V}}_{\mathrm{ir}}=\frac{V_{\mathrm{ir}}}{\sqrt{\underset{r=1}{\overset{K}{\mathrm{\Σ}}}{V_{\mathrm{ir}}}^2}},i=1,...,n,r=1,...K$

其中n为像素个数，i和r为矩阵V中元素的坐标，K为选取的特征值个数。

Images