CN102938072B - 一种基于分块低秩张量分析的高光谱图像降维和分类方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于分块的低秩张量分析的高光谱图像降维和分类方法。本发明根据高光谱图像的三维数据结构及其强烈的光谱特性和局部空间相关特性，将分块思想引入基于低秩张量分析的高光谱图像降维方法中，克服了图像的整体空间相关性较弱以及降维子空间维度的设定对降维效果的负面影响，最终得到一种能够大幅提高图像总体分类精度的新型降维方法——分块低秩张量分析法。该算法对各种不同的高光谱数据（包括仿真数据和实际数据集）都表现出良好的适用性。在基于高光谱遥感图像的高精度的地物分类方面具有重要的应用价值。

Description

一种基于分块低秩张量分析的高光谱图像降维和分类方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种基于分块的张量分析理论，解决高光谱遥感数据降维和分类问题的方法。

背景技术

遥感是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。近年来，随着高光谱成像技术的发展，高光谱遥感已经成为遥感领域一个快速发展的分支。作为一种多维信息获取技术，它将成像技术和光谱技术相结合，在电磁波谱的数十至数百个非常窄而且连续的光谱区间内同时获取信息，从而得到高光谱分辨率的连续、窄波段的图像数据，且每一个图像像元都可以提取一条完整连续的光谱曲线，大大扩展了图像解释地物的能力[1]。然而，由于高光谱图像（HyperspectralImagery,HSI）的相邻波段之间的波长差异往往只在纳米数量级，它们之间的相关系数一般都非常高。也就是说，在邻近的波段之间存在大量的冗余信息，很明显这会造成存储和处理能力的浪费。为了减少运算量、降低计算复杂度以及改善分类效率、提高分类精度，对HSI图像进行降维处理显得十分必要[2]。

当前应用广泛的主成分分析方法(PrincipalComponentAnalysis,PCA)，虽然易于实现，但必须将HSI每一个波段的二维图像化成向量，从而破坏了原图的空间分布特性，忽略了空间和光谱间的相互联系。该缺陷普遍存在于只能处理向量或者一维样本排列成的矩阵的降维方法。

要从根本上解决上述问题，就应该使降维方法契合HSI的数据结构，采用基于张量分析的降维方法。首先，用一个三阶张量存储高光谱数据，其中两个维度表示像素空间分布，余下的维度代表光谱分布，保存了HSI的空间分布特性，维护了空间和光谱间的关联性。其次，一个用于构建降维算法、与降维目的相适应的张量，不仅能挖掘物理意义明确的分量，更能充分利用空间和光谱间的关联性，将张量的价值最大化，从而提升降维效果。

基于张量分析的降维方法中，最具代表性的是Renard等人[2]的低秩张量分析法(LowerRankTensorAnalysis,LRTA)。LRTA是PCA在张量形式上的推广，建立在Tucker张量分解模型[3]的基础之上，将HSI投影到一个低秩的子空间进行降维，维护了原图的空间分布特性，在一定程度上利用了空间和光谱间的相互联系。该方法虽然考虑了空间与光谱间的相互作用，但它的局限性也很明显。首先，LRTA对于空间相关程度较弱的HSI并不奏效，只有强烈的空间相关性才能凸显出LRTA相对于PCA的优越性。而且，LRTA降维后的分类精度受子空间的影响非常大，现有的子空间估计方法（如Hysime、AIC、MDL[4-10]）却不可靠，往往不能找出使LRTA的性能达到最优的最佳子空间。

本发明是一种基于分块的低秩张量分析法(SubtensorbasedLowerRankTensorAnalysis,S-LRTA)，弥补了LRTA的缺陷。S-LRTA保持原来的光谱数据不变，将空间维划分成许多相同大小的子区域，每一子区域对应一个子张量，对每个子张量分别进行LRTA降维、分类，最终提高图像的整体分类精确度。分块的思想基于两点事实：在小范围的空间区域内，HSI对应的地物类别单一，像素相似，空间相关性尤为明显；而一般的HSI波段数目众多、波段宽度狭窄，光谱相关性非常明显。所以，只要将HSI张量的空间区域分块而无需改变原张量的光谱维，就可以最大限度地利用区域空间相关性和光谱间相关性。经S-LRTA降维后，每一子张量的子空间维度对分类精度的影响力度大为下降，从而无须在子空间的估计上耗费大量精力。即使是对整体空间相关性较弱的HSI，S-LRTA也能得到比LRTA更高的分类精度。

下面介绍与本发明相关的一些概念：

1.张量

张量即多维数组，可记为，其中称为维度或模式，表示阶数，R表示实数。向量与矩阵分别为一阶和二阶张量，乃张量的特例。

2.张量的模式n展开

对张量进行模式n展开（mode-nunfolding）即把张量矩阵化。已知一个三阶张量，做模式n展开得到展开矩阵，元素的位置坐标由变为，其中行坐标保持不变，列坐标的计算规则如下：

模式1：(1)

模式2：(2)

模式3：(3)。

3.张量和矩阵的模式d乘法

一个张量和一个矩阵在模式d上相乘（mode-ntensormatrixproduct），，它们的积定义为：

(4)

具体地，乘积中每个元素的计算规则如下：

(5)

其中，。通俗地说，求一个张量和矩阵的n阶乘积，就是先把这个张量展开成模式n，再与矩阵相乘，然后将所得矩阵转换回张量的形式。

4.Tucker张量分解模型

张量分解模型分为PARAFAC和Tucker两大类。本发明采用的是Tucker模型。已知一个张量，对进行Tucker模型分解，如式(6)所示：

(6)

其中是核张量；矩阵可称为压缩矩阵、因子矩阵、模式矩阵、投影矩阵等；表示近似误差。Tucker分解模型的特点在于它充分地利用了各模式分量的相互作用[3]。

发明内容

本发明的目的在于提出一种能够最大限度地利用区域空间相关性和光谱间相关性、不依赖降维子空间维度、与图像整体空间相关程度无关的高光谱图像降维方法。

本发明提出的高光谱图像降维方法，是一种分块低秩张量分析法(SubtensorbasedLowerRankTensorAnalysis,记为S-LRTA)，即S-LRTA的重点在于分块处理，具体发明内容是：将空间维划分成若干相同大小的子区域，每一子区域对应一个子张量，对每一子张量依次进行LRTA降维和分类。

已知原高光谱图像张量，其中，第一、二维表示空间，第三维表示光谱；已知分块尺寸为，其中，；已知每块子张量的降维子空间维度为，其中在张量分析过程中用到，而不在降维结果中体现[2]。已知地物类别的编号为。具体步骤如下：

步骤1、保持的第三维不变，将其第一、二维的边界分别作镜像扩展，得到，使得，，其中表示取余运算。

步骤2、在空间上划分，总共得到块大小相同的子张量，每块子张量，其中，表示整除运算，下标表示子张量的编号。分块的方法如图1所示。

步骤3、对每一子张量进行基于低秩张量分析的降维，得到各子张量降维结果。

步骤4、采用有监督分类器对进行分类时，具体操作如下：

（a）设定一个样本个数的下限阈值，找出训练样本个数大于既定阈值的子张量进行分类，其中表示已分类的子张量的编号的集合，表示已分类的子张量的编号，表示在所有已分类的子张量中的序列号，表示已分类的子张量的总数。

（b）对于原训练样本数量不足的子张量，

其中表示未分类的子张量的编号的集合，表示未分类的子张量的编号，表示在所有未分类的子张量中的序列号，表示未分类的子张量的总数。按照，每次取出一个未分类的子张量，依次考察其上、下、左、右四个方位上的子张量。如果某方位上存在分类完毕的子张量，则将待分类的子张量与之合并后形成新的张量（上下两个子张量合并）或（左右两个子张量合并）再做低秩张量分析降维，以分类完毕的像元作为训练样本对待分类的子张量进行分类。如图1右侧所示，深灰色的块代表将要分类的子张量（尺寸为），带星号的块表示与之相邻的子张量，如图1左侧所示，带框的两个子张量合并后形成新的子张量（尺寸变为）。

（c）综合（b)中当前待分类的子张量的每一个像素的所有分类结果（其中，表示像素的编号，表示与相邻的四个方位上所有分类完毕的子张量的编号，若存在个分类完毕的子张量，像素就有种分类结果），地物类别出现的频率最高就将该像素划归为此类；若有几种类别的频率相同，则参考与该子张量相邻的训练样本，训练样本中类别出现的频率最高就将该像素划归为此类；若仍不能判决，则随机选取训练样本中出现频率最高的地物类别中的一种作为该像素的类别；分类判决后，令，将该子张量的编号从中去除，并添加至。

（d）重复执行(a)，(b)，(c)，直至所有子张量都已分类完毕，即。

步骤5、将每一子张量的分类结果按其在空间上对应的位置合并，形成整幅图像的分类结果。

本发明的优点

本发明是一种新的基于张量分析和分块思想的降维和分类算法，不仅吸收了原张量降维方法的所有优点：如保持高光谱数据的空间分布特性、挖掘空间和光谱的联系，还省去了子空间维度的估计流程，并大幅提高了分类精度。理由如下：由于子张量涵盖的空间范围小、所包含的像素数少、像素的种类单一，则空间相关程度高；而且一般的HSI波段数众多、波段宽度狭窄，光谱相关性明显，无须划分原有的光谱，也能使每一子张量兼具强烈的空间相关性和谱间相关性。因此，在每一块子张量中，LRTA都能取得良好的降维效果，使得分类精度大幅提升，也就无需HSI具有明显的整体空间相关性。并且，每一子张量的子空间维度对分类精度的影响也大为下降，也就没有必要在子空间的估计上耗费精力，更何况现有的子空间估计算法本身并不可靠。

仿真和实地数据实验表明，我们的算法与LRTA算法相比明显具有更高的分类精度，无论HSI有无明显的整体空间相关性，即使任意设置子空间维度，S-LRTA的分类精度依然很高。因此，该发明具有重要的实际意义。对于实际高光谱遥感数据实验，该方法也得到了理想的结果，进一步证实了方法的有效性和对于各种不同数据的适用性。

附图说明

图1分块的方式及子张量的合并。

图2回形仿真HSI中，(a)9种地物的光谱曲线，（b）地物分布情况。

图3回形仿真HSI中，LRTA的降维效果与整体空间相关性的关系。(a)SVM,(b)KNN。

图4回形仿真HSI中，S-LRTA的设为(1,1,1)，LRTA的为Hysime，AIC或MDL算法的估计结果，比较S-LRTA和LRTA的降维效果。(a)SVM,(b)KNN。

图5带宽为8个像素的回形仿真HSI中，S-LRTA和LRTA的子空间设置相同时，比较两者的降维效果。

图6Indiana数据。(a)伪彩色图；(b)地物真实分布。

图7Indiana数据中，比较S-LRTA、LRTA以及PCA的降维效果。(a)SVM,(b)KNN。

图8Indiana数据中，当为(6,6)时，S-LRTA的降维效果与子空间维度的关系。

图9Indiana数据经过S-LRTA降维、SVM分类得到的地物分布图。及相应的总体分类精度为：(a)(2,2),100%,(b)(6,6),98.82%,(c)(10,10),92.37%。

图10Indiana数据中，设为(1,1,1)时，S-LRTA的降维效果与分块尺寸的关系。

具体实施方式

通过仿真数据和真实数据实验，将本发明提议的S-LRTA与现有的LRTA和PCA降维方法进行比较，证明S-LRTA的优越性。

为了量化HSI的“空间相关性”和“光谱间相关性”，本发明采用“平均相关系数”衡量相关程度，即对HSI张量在各模式上的展开矩阵分别求相关系数矩阵的平均值。这种方法简单客观。对于降维效果，则采用总体分类精度(OverallAccuracy,OA)[1]作为评判依据。已知地物真实，OA表示所有类别分类正确的样本点数量的总和的平均值，计算公式如(7)所示：

(7)

其中，共有类地物，表示类别的编号，表示本来为类别的测试样本被正确分到类别的总数，表示测试样本的总数。

实验采用两种有监督的分类器：支撑向量机（SupportVectorMachine，SVM）[11]和最近邻分类法（K-NearestNeighborhood，KNN）[12]分别测试降维方法的效果。对每一组仿真数据和真实数据，都随机抽取了其中样本的20%作为SVM或KNN的训练样本。SVM的核都采用径向基函数，其gamma参数均设为8。在S-LRTA的步骤4中，将判断子张量是否具有足够的训练样本的下限设为1。作图时，降维子空间的维度统一用表示，空间区域分块的尺寸用表示。

下面，分别用仿真数据和实际高光谱图像数据为例说明本发明的具体的实施方式:

1.仿真数据

我们采用仿真数据揭露LRTA的缺陷以及测试S-LRTA的性能。构造仿真HSI的光谱数据来源于USGC的数据库librarysplib06a，从中随机选取了9种类型地物的186个波段，光谱特性曲线如图2（a）所示。实验一共用到了4组回形仿真HSI，每幅图像由9层宽度相同的回形带组成，最里面的地物在空间上成正方形分布，第2种地物包裹其外，在空间上成回形带状分布，第3种地物再成回形带状分布并包裹于第2种地物之外，以此类推，9种地物层层包裹形成所需的仿真数据。4组HSI的尺寸为85×85×186、102×102×186、119×119×186以及136×136×186，相应的回形带宽分别为5、6、7、8个像素，每条回形随机取一种物质。以带宽为8个像素的仿真HSI为例，其地物分布如图2（b）所示。将为了模拟真实情况，还给每组仿真HSI加上同样的随机噪声，该噪声在光谱上呈高斯分布、方差为1/300、信噪比为40dB。4组仿真HSI的空间相关程度不同而光谱相关程度相同。

实验1整体空间相关性与LRTA性能关系实验为了单独考察整体空间相关性对LRTA降维效果的影响，令4组回形仿真HSI的子空间维度相同。综合图3和表1可以看到，回形带越窄，即整体空间相关性越弱，则LRTA降维效果越差，分类精度越低，说明LRTA的降维效果对空间相关性的变化极为敏感。如果HSI的空间相关程度不够高，则LRTA无法有效地处理空间信息。

表1平均相关系数与仿真HSI的回形带宽的对应关系

。

实验2子空间维度与LRTA性能关系实验图4的LRTA结果说明Hysime、AIC和MDL三种算法对子空间维度的估计均不可靠，因为对于每4组回形仿真数据，结合Hysime、AIC、MDL的LRTA产生的SVM和KNN分类精度均低于70%。以带宽为8的回形仿真数据为例，采用LRTA时，图5中两种分类精度随子空间维度的变化都极其明显。图4、5说明LRTA的性能极其依赖于子空间的维度，但是使LRTA有效的子空间难以用现有的估计算法确定。

实验3S-LRTA与LRTA性能比较实验从图4的S-LRTA结果可见，采用S-LRTA后，对于任意设定的子空间维度，4组回形仿真数据的SVM和KNN分类精度都达到了100%，说明整体空间相关性对分类精度的影响不再重要。以带宽为8的回形仿真数据为例，将其分成若干8×8×186的子张量，如图5所示，无论每块子张量的子空间维度如何改变，分类精度都是100%，证明子空间维度的设置对S-LRTA几乎无影响。对比图4和图5中S-LRTA与LRTA的曲线，可知S-LRTA的降维效果普遍优于LRTA。

由仿真实验结果可知，本发明提出的S-LRTA能够克服LRTA的缺陷：消除整体空间相关性以及子空间维度对降维效果的影响。

2.实际数据

本发明使用实际的高光谱图像数据对所提出S-LRTA的性能进行测试。

采用的实际数据是一幅Indiana地区的AVIRIS高光谱遥感数据[13]。该数据成像于1992年7月，包含220个波段，波长范围为，光谱分辨率为10nm，空间分辨率为17m，图像大小为145×145(共21025pixels)。该区域内主要的地物覆盖类型有：各种农作物（包括大豆、玉米、小麦等）、植被（包括树林、草地等）以及各种人工建筑（高速公路、铁塔、房屋等）。同时，美国Purdue大学研究组给出一份该地区实地勘测结果可供参考[14]。

首先剔除了坏波段1~4，78-82，103-115，148-166以及211-220，这些波段都是水吸收或者低SNR波段，然后将剩下的总共169个波段数据用于算法验证工作。取第9、17、117波段分别作为R、G、B分量合成伪彩色图，如图6(a)所示。根据端元数目确定算法的估计结果，并结合实地勘测结果，最终确定共取9类地物，真实分布情况如图6(b)所示，各类的名称及样本数如表2所示。

表2Indiana数据中的地物分类样本个数

类别	对应地物	样本个数
			C1	玉米未耕地(Corn-notill)	1256
C2	玉米疏耕地(Corn-min)	726
			C3	牧草(Grass/Pasture)	431
C4	林地(Grass/Tress)	626
			C5	干草(Hay-windrowed)	443
C6	大豆未耕地(Soybeans-notill)	828
			C7	大豆疏耕地(Soybeans-min)	2284
C8	大豆已耕地(Soybeans-clean)	503
			C9	树林(Woods)	1198

实验1整体空间相关性与S-LRTA性能关系实验由表3可知，Indiana数据的空间相关性不明显。而如图7所示，采用S-LRTA后，无论降维后的光谱维度怎样变化，SVM和KNN的分类精度都远远超过了LRTA和PCA使SVM和KNN达到的分类精度。可见，通过充分利用图像局部的空间相关性，S-LRTA不仅大幅提高了分类精度，还摆脱了对HSI的的整体空间相关性的依赖。

表3Indiana数据各模式上的平均相关系数

。

实验2子空间维度与S-LRTA性能关系实验如图8所示，无论子空间维度如何变化，S-LRTA都能使SVM和KNN的分类精度维持在98.75%以上，说明S-LRTA的降维效果受子空间维度的变化的影响甚微。所以，运用S-LRTA时，可以任意设定子空间维度，无需借助原本就不可靠子空间维度估计算法。

实验3分块尺寸与S-LRTA性能关系实验分块方法自身的瑕疵可能会带来负面效应，但不会影响全局的分类精度。如图9(c)所示，只有当分块的尺寸偏大时，分块负面效应才会显现，而且错分的像素只集中在两种或多种地物的交界处，绝大部分区域都分类正确。图10也表明，S-LRTA的降维效果良好并且对分块的大小并不敏感：只要分块的尺寸足够小，例如小于10，无论其怎样变化，总体分类精度都保持在95%以上。因此，可以将任意设置成一组小于(10,10)的值。但也不宜过小，以免过多增加程序耗时。

综上，无论是仿真数据还是实际数据的实验，都有力地验证了S-LRTA的优越性和实用性。S-LRTA的降维效果对分块尺寸、子空间维度的改变均不敏感，对HSI的空间相关程度要求不高，从而S-LRTA无需借助子空间维度的估计算法，能够适用于各种HSI。

参考文献

[1]ChangC-I.Hyperspectralimaging:techniquesforspectraldetectionandclassification[M].NewYork:Plenum,2003

[2]N.RenardandS.BourennaneDimensionalityreductionbasedontensormodelingforclassificationmethods[J].IEEETrans.onGeosci.RemoteSens.,47(4):`801-805April2009.

[3]A.Cichochi,R.Zdunek,A.H.PhanandS.AmariNonnegativematrixandtensorfactorizations[M].NewYork:Wiley,2009.

[4]M.WaxandT.KailathDetectionofsignalsbyinformationtheoreticcriteria[J]IEEETrans.onGeosci.RemoteSens.,33(2):387-392April1985

[5]J.Bioucas-DiasandJ.NascimentoHyperspectralsubspaceidentification[J]IEEETrans.onGeosci.RemoteSens.,46(8):2435-2445August2008

[6]C.ChangandQ.DuEstimationofnumberofspectrallydistinctsignalsourcesinhyperspectralimagery[J].IEEETrans.onGeosci.RemoteSens.,42(3):`608-619March2004.

[7]D.Letexier,S.BourennaneandJ.Blanc-TalonNonorthogonaltensormatricizationforhyperspectralimagefiltering[J].IEEEGeosci.RemoteSens.Letters,5(4):`3-7January2008.

[8]N.Renard,S.BourennaneandJ.Blanc-TalonMultiwayfilteringappliedonhyperspectralimages[J].ACIVS2006,LNCS4179:127-137,2006

[9]D.LetexierandS.BourennaneEstimationofn-moderanksofhyperspectralimagesfortensordenoising[J].EUSIPCO2009,2594-2597

[10]A.HuckandM.GuillaumeEstimationofthehyperspectraltuckerranks[J]ICASSP2009,1281-1284

[11]MelganiF,andBruzzoneL.Classificationofhyperspectralremotesensingimageswithsupportvectormachines[J].IEEETrans.onGeosci.RemoteSens.,42(8):1778-1790,Aug.2004.

[12]CoverT,andHartP.Nearestneighborpatternclassification[J].IEEETrans.onInf.Theory,13(1):21-27,Jan.1967.

[13]http://cobweb.ecn.purdue.edu/~biehl/Multispec/documentation.html

[14]Landgrebe.MultispectralDataAnalysis:ASignalTheoryPerspective[R].SchoolofElectr.Comput.Eng.,PurdueUniv.,WestLafayette,IN,1998.。

Claims (1)

1.一种基于分块低秩张量分析的高光谱图像降维和分类方法，其特点在于用三阶张量描述高光谱图像，将空间维划分成若干相同大小的子区域，而保持原来的光谱维不变，每一子区域对应一个子张量，对每一子张量依次进行基于低秩张量分析的降维和分类；

已知原高光谱图像张量，其中，第一、二维表示空间，第三维表示光谱；已知分块尺寸为，其中，；已知每块子张量的降维子空间维度为，其中在张量分析过程中用到，已知地物类别的编号为；其中：

表示空间域的维度，即为行数、为列数，表示光谱域的维度，即波段数”；“表示子张量空间域的维度，即为行数、为列数”；“为子张量空间域降维后的维度，为光谱域降维后的维度；

具体步骤如下：

步骤1、保持的第三维不变，将其第一、二维的边界分别作镜像扩展，得到，使得，，其中表示取余运算；

步骤2、在空间上划分，总共得到块大小相同的子张量，每块子张量，其中，表示整除运算，下标表示子张量的编号；

步骤3、对每一子张量进行基于低秩张量分析的降维，得到各子张量降维结果；

步骤4、采用有监督分类器对进行分类时，具体操作如下：

(a)设定一个样本个数的下限阈值，找出训练样本个数大于既定阈值的子张量进行分类，其中表示已分类的子张量的编号的集合，表示已分类的子张量的编号，表示在中的编号，表示已分类的子张量的总数；

(b)对于原训练样本数量不足的子张量其中表示未分类的子张量的编号的集合，表示未分类的子张量的编号，表示在中的编号，表示未分类的子张量的总数；按照，每次取出一个未分类的子张量，依次考察其上、下、左、右四个方位上的子张量，如果某方位上存在分类完毕的子张量，则将待分类的子张量与之合并后形成新的张量或再做低秩张量分析降维，以分类完毕的像元作为训练样本对待分类的子张量进行分类；

(c)综合(b)中当前待分类的子张量的每一个像素的所有分类结果，其中，表示像素的编号，表示与相邻的四个方位上所有分类完毕的子张量的编号，若存在个分类完毕的子张量，像素就有种分类结果，地物类别出现的频率最高就将该像素划归为此类；若有几种类别的频率相同，则参考与该子张量相邻的训练样本，训练样本中类别出现的频率最高就将该像素划归为此类；若仍不能判决，则随机选取训练样本中出现频率最高的地物类别中的一种作为该像素的类别；分类判决后，令，将该子张量的编号从中去除，并添加至；

(d)重复执行(a),(b),(c)，直至所有子张量都已分类完毕，即；

步骤5、将每一子张量的分类结果按其在空间上对应的位置合并，形成整幅图像的分类结果。