CN105957112A - 基于快速uncls的高光谱亚像素探测方法 - Google Patents

Abstract

本发明提出了一种基于快速UNCLS高光谱亚像素探测方法，用于解决存在的计算性能低的问题。包括如下步骤：输入高光谱图像，得到该图像的像素矢量矩阵；初始化相关参数，计算首个目标特征矢量；用光谱特征矩阵存储目标特征矢量；计算相关矩阵，并利用矩阵升序求逆方法对其求逆；利用矩阵降序求逆方法对引导矩阵求逆，得到修正的NLCS；利用修正的NLCS和相关矩阵的逆矩阵计算丰度矩阵；判断丰度矩阵的拟合优度是否小于设定阈值，满足则结束，否则计算新的目标特征矢量，执行步骤3。通过采用矩阵升序求逆和矩阵降序求逆，对现有基于UNLCS高光谱亚像素的探测方法中UNCLS进行修正，有效地提升了高光谱亚像素探测的计算性能。

Description

基于快速UNCLS的高光谱亚像素探测方法

技术领域

本发明属于图像处理技术领域，涉及一种基于快速UNCLS的高光谱亚像素探测方法，可用于对高光谱图像亚像素探测的计算性能进行有效提升。

背景技术

从20世纪60年代开始，随着信息技术、传感器技术和空间科学的进步，遥感技术得到了迅速的发展。遥感是从远处探测自然界物体的一种技术，自然中的一切物体都具有反射或发射电磁波的特性，物体不同，所处于的环境不同，物体发射或反射的电磁波也不相同。被物体发射或反射的电磁波中包含着物体的特征信息，能够作为感知不同物体的工具。基于此原理，遥感技术通过物理设备收集物体发射或反射的电磁波来感知和判别物体。通过遥感获得的图像就是遥感图像。

国际遥感界的共识是光谱分辨率在λ/10数量级范围的称为多光谱，这样的遥感器在可见光和近红外光谱区只有几个波段；而光谱分辨率在λ/100的遥感信息称之为高光谱遥感；随着遥感光谱分辨率的进一步提高，在达到λ/1000时，遥感即进入超高光谱阶段。

多光谱图像是我们平常接触较多的遥感影像,其在电磁波谱中的几个可见光波段上同时获取信息。高光谱遥感图像在电磁波谱的紫外、可见光、近红外和中红外区域中的数十至数百个非常窄的且连续的光谱段上同时获取信息。相比于传统的遥感图像如多光谱图像来说，高光谱图像可以提供丰富的信息而且有很高的光谱分辨率。

高光谱遥感信号是通过成像光谱仪获取的，由于通常高光谱图像的空间分辨率较低以及受到目标尺寸和地物复杂特性的影响，使得感兴趣目标在高光谱图像中往往处于亚像元级或弱信息状态。所谓的亚像元级目标是指尺寸小于遥感器空间分辨率的目标，而弱信息目标也叫小存在概率目标，即在图像中分布很少的目标。因此，传统的基于空间形态的目视解译方法无法实现对这类目标的探测识别。与传统的基于高空间分辨率遥感影像的目标探测原理不同，高光谱遥感目标探测主要依据目标与地物在光谱特征上存在的差异进行检测识别。虽然高光谱数据空间分辨率不高，但其丰富的光谱信息弥补了这一不足，而且在足够高的光谱分辨率下，感兴趣目标会表现出诊断性光谱特征或者会在背景地物中显示为一种“数据异常”，利用这些信息完全可以将亚像元级目标或小目标提取出来。

成像光谱仪所记录的地表物质的反射信号是以像元为单位标记的。地表物质的光谱信号是由多个像元组成的。如果像素仅包含一种类型的地表物质，则称此像素为纯像元或者端元；若该像素包含两种以上的地表物质，则称混合像素。高光谱图像中混合像素的存在，使得传统的基于像素级高光谱图像的分类、压缩和目标检测等应用受到限制。所以需要解决混合像素的分解问题。高光谱模型分为线性混合和非线性混合。在分析高光谱含有的混合像素时，一般应用线性光谱混合模型进行分析。基于线性混合模型的高光谱解混主要包括两部分：端元提取和混合像素分解。

线性混合模型可表示为矩阵形式:

R＝Ma+n (1)

其中，M＝[m₁,m₂,...,m_p]是光谱特征矩阵，a＝[a₁,a₂,..,a_p]^T是每一个端元向量占有的丰度矩阵。n为一个L维的误差项，p个端元向量和对应的丰度矩阵都是未知量，根据实际的高光谱图像端元光谱和丰度矩阵分布情况，对于光谱和它对应的丰度矩阵有两个约束条件：

端元光谱及其丰度矩阵是非负的，即

m_i≥0,0≤a_i≤1 (2)

各个端元对应的丰度矩阵总和是1，即

$\underset{i=1}{\overset{p}{\Σ}}{a}_i=1,i=1,...,p---\left(3\right)$

基于算法实现的端元提取有：纯像元指数(Pixel Purity Index,PPI)是一种使用广泛成熟的方法，是由Boardman等研究人员在1995年基于凸单兴体理论提出并发展成熟的；N-FINDR利用高光谱图像在特征空间中呈凸面单形体的特性，通过搜索体积最大的单形体来寻找图像中的所有端元；迭代误差分析法(IEA)不需要对原始的高光谱图像进行降维或去冗余处理，可以直接对图像端元提取；顶点成分分析(Vertex Component Analysis,VCA)是在假设高光谱图像中存在纯像元的情况下基于凸面几何理论提出的。

基于算法实现的混合像素分解有以下典型方法：

1)基于最小二乘法(LS)实现混合像素分解，是基于线性混合模型在无约束条件下均方误差意义的最佳解，是目前应用最广泛的算法。但是该解的表达式多次用到矩阵的乘积及求逆，特别是在端元数目较多的情况下，使得计算复杂度比较大。另外，由于没有附加任何约束条件，混合像素分解效果不理想，会出现端元成分小于0或者大于1的情况。

结合非负和丰度矩阵和为一的约束条件提出了基于和一约束最小二乘方法(SCLS)、非负约束最小二乘方法(NLCS)、全约束最小二乘方法(FLCS)实现的混合像素分解。

2)基于投影寻迹实现的混合像素分解，是一种专门处理高维数据的降维方法。其基本思想是把高维数据投影到低维空间上,用一个设计的投影指标(PI),作为度量来寻找能反映原高维数据结构或特征的投影称为“令人感兴趣”的投影,然后通过分析和研究投影数据以达到了解原数据的目的。

3)基于正交子空间投影OSP实现的混合像素分解,认为所有像元均由感兴趣的目标和背景组成,在消除背景地物的同时,得到目标地物的最佳匹配效果。OSP的特点在于通过信号光谱的逐步分离来提取出感兴趣的信号。

4)还有基于最大似然估计，光谱信号匹配，光谱角映射，单形体体积法等算法实现的混合像素分解。

近年来，随着盲信号处理技术的兴起和发展，非监督的光谱解混称为研究热点。非监督的光谱解混是指在端元信息，如图像中端元的个数，端元光谱等。完全未知的情况下，将混合像素分解为端元向量及其在像素中所占的比例。对于非监督的光谱解混问题来说，端元向量及其在混合像素中所占的比例(丰度矩阵)都是未知的。传统的方法认为端元的丰度矩阵是未知的非随机参数，用统计的方法来估算，如最小二乘估计，学者Chang等提出基于UNCLS的高光谱亚像素探测方法。该方法是实现高光谱亚像素探测的常用算法，该方法中利用线性混合模型描述高光谱图像的像素矢量，利用线性回归模型实现丰度矩阵估计，利用拉格朗日乘子实现丰度矢量的非负约束，实现高光谱亚像素的探测。基于UNCLS的高光谱亚像素探测方法的优点是无论在没有先验的基础上执行多少次，都可以自动提取目标特征生成一个稳定的解。然而该方法依赖于传统的矩阵求逆方法，矩阵求逆的时间成本高，使得该方法的计算性能低。

发明内容

本发明的目的在于克服上述现有技术存在的缺陷，提出了一种基于快速UNCLS高光谱亚像素探测方法，通过对现有基于UNCLS的高光谱亚像素探测方法，分别采用矩阵升序求逆和矩阵降序求逆，实现对无监督非负最小二乘法UNCLS进行修正，用于解决现有技术存在的计算性能低的技术问题。

为实现上述目的，本发明采取的技术方案，包括如下步骤：

步骤1，读入原始高光谱图像，得到该图像的像素矢量矩阵；

步骤2，设定迭代变量k并初始化为1,遍历所述图像的像素矢量矩阵，并将遍历得到的第k个目标特征矢量进行存储，得到当前光谱特征矩阵其中是像素矢量矩阵中的元素，是像素矢量矩阵中元素的拟合优度；

步骤3，对得到的光谱特征矩阵进行矩阵运算，得到相关矩阵

步骤4，利用矩阵升序求逆公式，对得到的相关矩阵求逆，得到相关矩阵的逆矩阵R_k ^-1，具体实现步骤如下：

步骤4a,设定得到的相关矩阵的光谱特征表达式 $R_k=\left[\begin{array}{cc}{M}_{k-1}^T{M}_{k-1}& M_{k-1}^T{t}_k\\ t_k^T{M}_{k-1}& t_k^T{t}_k\end{array}\right],$

其中，M_k-1是上级光谱特征矩阵

步骤4b，利用矩阵升序求逆公式，对设定得到的相关矩阵的光谱特征表达式中的各分块矩阵进行替代，得到相关矩阵的逆矩阵R_k ^-1的矩阵表达式

其中，

$R_{11}={\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}+r_{22}\·{\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}{M}_{k-1}^T{t}_k{t}_k^T{M}_{k-1}^T{\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}$

$r_{12}=-r_{22}\·{\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}{M}_{k-1}^T{t}_k$

$r_{21}=r_{12}^T$

步骤5，设定相关矩阵的矩阵表达式

步骤6，采用矩阵降序求逆公式，对设定的相关矩阵的矩阵表达式进行运算，实现非负约束最小二乘法NCLS的修正；

步骤7，利用修正后的NCLS对光谱特征矩阵中目标特征矢量的丰度矩阵进行计算，得到

步骤8，设定所述图像的像素矢量矩阵中元素的拟合优度误差阈ε；

步骤9，利用得到的目标丰度矩阵计算所述图像的像素矢量矩阵中所有元素的拟合优度计算公式表示如下：

步骤10，判断得到的拟合优度和设定的拟合优度误差阈ε是否满足若是则算法结束，否则执行步骤11；

步骤11，对设定的迭代变量k加1，并利用得到的拟合优度计算第k个目标特征矢量得到光谱特征矩阵执行步骤3。

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

本发明由于在提取目标特征矢量时，采用了矩阵升序求逆的方法；在生成丰度矩阵时，采用了矩阵降序求逆的方法，实现了对UNCLS的修正，与现有技术中仅采用的矩阵求逆方法相比，缩短了计算时间，有效地提升了高光谱亚像素探测的计算性能。

附图说明

图1为本发明的实现流程框图；

图2为本发明与现有基于UNCLS高光谱亚像素探测的时间成本仿真对比图。

具体实施方式：

以下结合附图和具体实施例，对本发明作进一步详细说明：

参照附图1，本发明包括如下步骤：

步骤1，读入原始高光谱图像，得到该图像的像素矢量矩阵。读入的高光谱图像是通过AVIRIS(Airborne Visible/Infrared Imaging Spectrometer,机载可见光/红外成像光谱仪)采集得到。高光谱图像是数据立方体，其中包括二维空间的方向维和光谱维，分别携带高光谱图像的空间和光谱信息。该高光谱图像的像素矢量矩阵反映了图像中像元的光谱特征。

步骤2，设定迭代变量k并初始化为1,遍历所述图像的像素矢量矩阵，并将遍历得到的第k个目标特征矢量进行存储，得到当前光谱特征矩阵其中是像素矢量矩阵中的元素，是像素矢量矩阵中元素的拟合优度。

目标特征矢量的赋值公式是一个分段函数，在迭代变量k＝1，即第一次迭代时遍历所述图像利用分段函数中对应的公式k＝1得到第一个目标特征矢量，对其进行存储得到第一级光谱特征矩阵，此时光谱特征矩阵中只含有一个目标特征矢量；在迭代的过程中，通过将每次迭代产生的新目标特征矢量k＞1加入光谱特征矩阵实现对光谱特征矩阵的更新。

拟合优度是一个统计术语，是检验生成的模型对样本的拟合程度。在高光谱图像处理中的拟合优度是对混合像元解混后结果的评价标准。

步骤3，对得到的光谱特征矩阵进行矩阵运算，得到相关矩阵相关矩阵是光谱特征矩阵M_k的转置矩阵和本身的矩阵运算。

步骤4，利用矩阵升序求逆公式，对得到的相关矩阵求逆，得到相关矩阵的逆矩阵R_k ^-1，具体实现步骤如下：

步骤4a,设定得到的相关矩阵的光谱特征表达式其中，M_k-1是上级光谱特征矩阵在矩阵升序求逆公式中需要将矩阵表示为分块矩阵的形式进行运算，因此设定相关矩阵的光谱特征表达式。

步骤4b，利用矩阵升序求逆公式，对设定得到的相关矩阵的光谱特征表达式中的各分块矩阵进行替代，得到相关矩阵的逆矩阵R_k ^-1的矩阵表达式分块矩阵中r₂₂、R₁₁是方阵，r₂₁等于r₁₂的转置，因为相关矩阵的光谱特征表达式中的分块矩阵是由迭代中已知的变量组成，因此与传统的直接求逆方法相比，矩阵升序求逆方法可将求逆的过程大大简化。

步骤5，设定相关矩阵的矩阵表达式在矩阵降序求逆公式中，也需要将矩阵表示为分块矩阵的形式进行运算，因此设定相关矩阵的矩阵表达式。

步骤6，采用矩阵降序求逆公式，对设定的相关矩阵的矩阵表达式进行运算，实现非负约束最小二乘法NCLS的修正。

其中涉及到NCLS算法，NCLS是UNCLS中用于求取目标特征矢量丰度矩阵的算法，实现步骤如下：

步骤6a，初始化非负索引集P^(k)＝{1,2,...,p}，负索引集R^(k)＝φ，其中p表示目标特征矢量的数量，非负索引集用来存放非负丰度矩阵的索引值，负索引集用来存放负的丰度矩阵矢量的索引值，在迭代中按照丰度矩阵的正负将其索引值放置在对应的索引集中，实现对非负索引集和负索引集的更新。

步骤6b，计算线性回归模型并判断是否所有的元素都非负，满足算法终止，否则执行步骤6c，线性回归模型是用来估计丰度矩阵引入的变量。

步骤6c，将非负索引集P^(k)中与得到的线性回归模型中负元素的索引值对应的元素移动到负索引集R^(k)中。

步骤6d，将得到的线性回归模型与得到的负索引集R^(k)相对应的元素进行复制，得到负丰度矩阵

步骤6e，将相关矩阵中与非负索引集P^(k)元素对应的行和列删除，得到引导矩阵

步骤6f，利用矩阵降序求逆公式对得到的引导矩阵求逆，得到引导矩阵的逆矩阵

步骤6g，利用得到的引导矩阵的逆矩阵和得到的负丰度矩阵进行计算，得到乘子矢量因目标特征矢量丰度矩阵需满足非负约束，所以引入乘子矢量。

步骤6h，判断得到的乘子矢量中是否所有元素为负，满足执行步骤，否则执行步骤6i。

步骤6i，将得到的乘子矢量中对应的索引值从负索引集R^(k)移动到非负索引集P^(k)，其中表示λ^(k)的第j个元素。

步骤6j，删除得到的非负索引集P^(k)在得到的相关矩阵中对应的列，得到辅助矩阵Ψ^(k)。

步骤6k，利用得到的辅助矩阵Ψ^(k)、乘子矢量和线性回归模型进行计算，得到非负丰度矩阵

步骤6l,判断得到的非负丰度矩阵是否所有的元素都非负，满足算法终止，否则执行步骤6c。

步骤7，利用修正后的NCLS对光谱特征矩阵中目标特征矢量的丰度矩阵进行计算，得到计算得到的丰度矩阵满足非负约束：

步骤8，设定所述图像的像素矢量矩阵中元素的拟合优度误差阈ε，拟合优度阈值ε用来实现终止条件的设定，该值是通过大量实验和统计得到的值。

步骤9，利用得到的目标特征矢量的丰度矩阵计算所述图像的像素矢量矩阵中所有元素的拟合优度计算公式表示如下：

步骤10，判断得到的拟合优度和设定的拟合优度误差阈ε是否满足若是则算法结束，否则执行步骤11。

步骤11，对设定的迭代变量k加1，并利用得到的拟合优度计算第k个目标特征矢量将得到第k个目标特征矢量加入光谱特征矩阵实现光谱特征矩阵的更新，然后执行步骤3，进行下一次迭代。

以下结合仿真实验，对本发明的技术效果作进一步的说明。

1、仿真条件和内容：

本实施例采用处理器为Inter Core i5-3210M，主频为2.5GHz，内存4GB的硬件测试平台，64位Windows 7操作系统的Matlab R2013a的软件平台，采用WTC的真实高光谱图像，对本发明与现有基于UNCLS高光谱亚像素探测的时间成本进行仿真对比，其结果如图2。

2、仿真结果分析：

参照图2，横坐标是目标特征矢量的个数，纵坐标从下至上是实现高光谱亚像素探测的运算时间和加速比。随着目标特征矢量的数目增加，本发明的时间成本呈稳定状态，对应速度保持稳定；现有基于UNCLS高光谱亚像素探测的时间成本呈现上升趋势，对应速度逐渐降低。

仿真结果表明，相对于现有基于UNCLS高光谱亚像素探测方法，本发明的时间成本明显降低，计算性能明显提升。

Claims (2)

1.一种基于快速UNCLS的高光谱亚像素探测方法，其特征在于，包括有如下步骤：

(1)读入原始高光谱图像，得到该图像的像素矢量矩阵；

(2)设定迭代变量k并初始化为1,遍历所述图像的像素矢量矩阵，并将遍历得到的第k个目标特征矢量进行存储，得到当前光谱特征矩阵其中是像素矢量矩阵中的元素，是像素矢量矩阵中元素的拟合优度；

(3)对得到的光谱特征矩阵进行矩阵运算，得到相关矩阵

(4)利用矩阵升序求逆公式，对得到的相关矩阵求逆，得到相关矩阵的逆矩阵R_k ^-1，具体实现步骤如下：

(4a)设定得到的相关矩阵的光谱特征表达式

其中，M_k-1是上级光谱特征矩阵

(4b)利用矩阵升序求逆公式，对设定得到的相关矩阵的光谱特征表达式中的各分块矩阵进行替代，得到相关矩阵的逆矩阵R_k ^-1的矩阵表达式

其中，

$R_{11}={\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}+r_{22}\·{\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}{M}_{k-1}^T{t}_k{t}_k^T{M}_{k-1}^T{\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}$

$r_{12}=-r_{22}\·{\left(M_{k-1}^T{M}_{k-1}\right)}^{-1}{M}_{k-1}^T{t}_k$

$r_{21}=r_{12}^T$

(5)设定相关矩阵的矩阵表达式

(6)采用矩阵降序求逆公式，对设定的相关矩阵的矩阵表达式进行运算，实现无监督非负约束最小二乘法UNCLS的修正；

(7)利用修正后的UNCLS对光谱特征矩阵中目标特征矢量的丰度值进行计算，得到

(8)设定所述图像的像素矢量矩阵中元素的拟合优度误差阈ε；

(9)利用得到的目标丰度值计算所述图像的像素矢量矩阵中所有元素的拟合优度计算公式表示如下：

(10)判断得到的拟合优度和设定的拟合优度误差阈ε是否满足若是则算法结束，否则执行步骤(11)；

(11)对设定的迭代变量k加1，并利用得到的拟合优度计算第k个目标特征矢量得到光谱特征矩阵执行步骤(3)。

2.根据权利要求1所述的一种基于快速UNCLS的高光谱亚像素的探测方法，其特征在于，步骤(6)中所述的运算，采用矩阵降序求逆公式(P)^-1＝(A-BD^-1C)^-1，对设定的相关矩阵的矩阵表达式进行求逆，得到无监督非负约束最小二乘法UNCLS的修正中引导矩阵的逆矩阵其中P表示降维矩阵。