CN109727280A - 一种基于正交基的高光谱图像丰度估计方法 - Google Patents

Abstract

一种基于正交基的高光谱图像丰度估计方法，涉及高光谱图像处理技术。以解决现有方法存在运算复杂度的问题。该算法利用Gram‑Schmidt方法计算端元向量组得到对应的正交基组，求解解混方程组，解得每个正交基的特征向量。将待解混光谱向量向特征向量投影，计算该投影向量的长度与正交基的长度比，得到正交基所代表端元的丰度估计。经过不同算法的比较分析，该算法只需进行向量内积运算，减小了运算复杂度，缩短高光谱图像丰度估计时间，提高了高光谱图像丰度估计效率。通过模拟数据和实际图像数据实验，验证了算法的有效性。

Description

一种基于正交基的高光谱图像丰度估计方法

技术领域

本发明涉及一种高光谱图像的丰度估计方法，涉及高光谱图像处理技术。

背景技术

高光谱图像包含丰富的空间信息和光谱信息，高光谱图像处理技术应用于地理、地质、农业、林业等领域。地物分布往往是复杂多样的，光谱成像仪的空间分辨率具有有限性，导致高光谱图像存在大量含有多种地物类型的混合像元。混合像元分解是高光谱图像处理的重要问题。混合像元分解包括端元提取和丰度估计。

像元混合模型主要包括2大类：线性光谱混合模型和非线性光谱混合模型。线性光谱混合模型是高光谱图像处理中应用最广泛的模型。根据设计思路解混算法可以分为凸面几何学算法和统计学算法。根据满足丰度约束条件的程度解混算法可分为无约束光谱解混算法、“和为1”约束光谱解混算法和“非负”约束光谱解混算法、全约束光谱解混算法。其中，无约束线性光谱解混方法是解混算法的根本，很多方法是在它的基础上发展起来的。普遍使用的无约束线性解混算法主要有：无约束最小二乘算法(unconstrained least squareerror,UCLS)、正交子空间投影算法(orthogonal subspace projection,OSP)、单形体体积算法(simplex volume,SV)、正交向量投影算法(orthogonal vector projection,OVP)。最小二乘法通过计算误差向量2范数的平方的最小值求得丰度估计矩阵，运算速度较快，是典型的丰度估计算法。该方法考虑全部端元的作用，不能单独分析感兴趣的端元。OSP算法将感兴趣的信号从背景信号中分离出来，逐个计算感兴趣端元的丰度估计。OSP算法涉及矩阵求逆运算，计算量大于UCLS算法。SV算法用待解混像元替换某个端元所得单形体与原单形体的体积比计算像元中端元的丰度，逐个得到端元的丰度。该算法物理意义明确，但需要进行行列式运算。OVP算法基于OSP思想，避免矩阵计算，减小计算量。

UCLS算法和OSP算法涉及矩阵求逆运算，SV算法需要行列式运算。OVP算法克服了UCLS算法、OSP算法和SV算法在计算复杂度上的不足。

CN103413292A的现有技术公开了一种基于约束最小二乘的高光谱图像非线性丰度估计方法。该发明通过在目标函数中引入丰度的非负及和为一约束以及非线性参数的有界约束，将高光谱图像非线性解混问题转化为求解丰度矢量和非线性参数的约束非线性最小二乘问题。进而，该发明采用一种交替迭代优化算法求解该问题。该发明从高光谱观测像素的混合模型出发，结合模型中丰度和非线性物理意义，有效地克服了线性解混的不足，同时具有良好的抗噪声性能，可以作为一种解决高光谱遥感图像非线性解混的有效手段。文献号为CN105976404A的现有技术公开了一种高光谱遥感图像的线性全约束丰度估计方法，构建原光谱解混问题的对偶问题，通过全约束最小二乘法对丰度进行初始化，结合偏离度ρ因素利用原始对偶内点法优化丰度结果，具体包括以下步骤：S1:建立图像的线性光谱混合模型：S2:对线性光谱混合模型施加约束条件；S3:采用原始对偶内点法和全约束最小二乘丰度估计算法对线性光谱混合模型进行混合像元的丰度估计。该发明有效地克服了全约束最小二乘法精度不高的问题，同时也解决了原始对偶内点法速度慢的问题。但现有技术中针对运算复杂度的问题没有提及。

发明内容

本发明的目的是提供一种基于正交基的高光谱图像丰度估计方法，以解决现有方法存在运算复杂度的问题。

本发明为解决上述技术问题采取的技术方案是：

一种基于正交基的高光谱图像丰度估计方法，所述方法的实现过程为：

利用Gram-Schmidt方法计算高光谱图像的端元向量组得到对应的正交基组，求解解混线性方程组(正交基组与解混方程组相关，正交基是解混方程组的组成部分，求解解混线性方程组为了得到特征向量)，解得每个正交基的特征向量；将待解混光谱向量向特征向量投影得到投影向量，计算该投影向量的长度与正交基的长度比，得到正交基所代表端元的丰度估计；高光谱图像中的多个端元向量组成所述端元向量组。

进一步地，所述方法的具体实现过程为：

将端元向量组S₁,S₂,…,S_P利用Gram-Schmidt方法进行正交化，得到正交基β₁,β₂,…,β_P

其中，i为端元序数，j为正交基序数；S_i为端元向量；

经推导，

X为待解混光谱向量(像元矢量)，a₁,a₂,…,a_P为端元向量组S₁,S₂,…,S_P对应的丰度，即a_i为S_i的丰度；

由于β₁,β₂,…,β_P两两正交，可得

联立P个等式，可得解混线性方程组为

令d_i为正交基β_i的特征向量，即d_i为正交基β_i代表端元S_i在正交基空间的投影向量，则

则

表示正交基所代表端元的丰度估计。

进一步地，所述正交基组与解混方程组相关，正交基是解混方程组的组成部分，求解解混线性方程组为了得到特征向量。

进一步地，正交基所代表端元的丰度估计是将待解混光谱向量向正交基空间进行投影，由投影点在正交基上的解混系数获得待解混光谱向量的每个端元的对应丰度估计。

进一步地，对端元的个数稍作增减重复进行丰度估计时，采用正交基所代表端元的丰度估计(采用OB算法)中，式(10)和式(12)中仅需要计算新添加的项或者去掉多余项。

本发明的有益效果是：

本发明在OVP算法基础上，提出基于正交基的丰度估计算法(orthogonal bases,OB)，进一步减小了运算复杂度。本发明所述方法相较于现有方法(OVP,SV,OSP,UCLS算法)减少了运算量，大大缩短了高光谱图像丰度估计时间，提高了高光谱图像丰度估计效率。

大多数基于线性混合模型的混合像元分解方法在进行丰度估计时，涉及矩阵运算(矩阵求逆或方阵行列式的计算)，导致计算复杂度高。当端元数增加时，算法的计算量也会呈指数级增长。本发明避免矩阵运算的快速丰度估计算法计算复杂度低，可以提高运算效率。

丰度估计是高光谱图像混合像元分解的重要步骤，对定量分析地物成分和后续图像分类起到重要作用。针对大多数线性光谱混合模型丰度估计算法需要进行行列式运算和矩阵求逆运算，导致算法计算复杂度高的问题，提出基于正交基的丰度估计算法。该算法利用Gram-Schmidt方法计算端元向量组得到对应的正交基组，求解解混方程组，解得每个正交基的特征向量。将待解混光谱向量向特征向量投影，计算该投影向量的长度与正交基的长度比，得到正交基所代表端元的丰度估计。经过不同算法的比较分析，该算法只需进行向量内积运算，明显减小了运算复杂度。通过模拟数据和实际图像数据实验，验证了算法的有效性。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1是Cuprite Nevada地区的AVIRIS图像(193波段)，Cuprite Nevada表示内华达州Cuprite矿区；AVIRIS为专有名词，表示机载可见光/红外成像光谱仪；

图2是选择端元光谱图；

图3是模拟数据图像示意图；图4为Hydice图像9个端元示意图，图中：(a)为4种类型地物图，(b)为端元选择区域图。

具体实施方式一：如图1至4所示，本实施方式对所述的基于正交基的高光谱图像丰度估计方法的实现及与现有方法的对比进行如下阐述：

1、UCLS算法和OSP算法涉及矩阵求逆运算，SV算法需要行列式运算。OVP算法克服了UCLS算法、OSP算法和SV算法在计算复杂度上的不足。

2、无约束线性解混算法

2.1线性混合模型

目前研究较多的线性光谱混合模型可以表示为

X＝SA+N (1)

其中，像元矢量X为L行，端元矩阵S＝[S₁,S₂,…,S_P]为L×P维，丰度矢量A＝[a₁,a₂,…,a_P]^T为P行，T为矩阵转置运算，噪声N为L行。

2.2UCLS算法

UCLS算法是基于线性光谱混合模型最常用的丰度估计算法。根据LSMM的假设，线性方程组(1)中方程的数量L大于未知数的数量P，因此属于超定方程组。可以找到一个A，使得SA与X的差距尽可能的小，这个差距通常用向量2范数的平方表示，则最优解可按式(2)获得：

UCLS算法的丰度估计为

A^UCLS＝(S^TS)^-1S^TX (3)

2.3OSP算法

像元X包括目标端元d和背景端元U，则

X＝da_i+Ub+N (4)

其中，d＝S_i，S_i的丰度为a_i，1≤i≤P，U＝(S₁,…,S_i-1,S_i+1,…,S_P)，b＝(a₁,…,a_i-1,a_i+1,…,a_P)为U的丰度。

设U的正交补空间为U^⊥，将U^⊥作用于式(4)，则

U^⊥X＝U^⊥da_i (5)

其中，U^⊥＝I-UU^#，U^#为U的伪逆矩阵，U^#＝(U^TU)^-1U^T，(U^TU)^-1为(U^TU)的逆矩阵，单位矩阵I为P×P维。

将d^T作用于式(5)，则OSP算法的丰度估计为

2.4SV算法

根据线性光谱混合模型的几何学描述，可以用待解混像元替换某个端元所得单形体与原单形体的体积比计算像元中端元的丰度。

端元矩阵S＝[S₁,S₂,…,S_P]在L维特征空间中构成的P-1维单形体体积可根据式(7)计算

其中，S′＝[S₁-S_P,S₂-S_P,…,S_P-1-S_P]，将S中的S_i(1≤i≤P)替换为像元X，则X中第i个端元的丰度估计为

式中，当i＝P时，则令S_i′＝[S₂-S₁,S₃-S₁,…,X-S₁]。

2.5OVP算法

OVP算法基于OSP思想，利用Gram-Schmidt正交化计算端元的丰度估计。对端元集合S＝[S₁S₂…S_P]进行Gram-Schmidt正交化，得到正交基组按照不同端元顺序进行正交化，得到不同的正交基。将某个端元向量d＝S_P作为最后向量，得到的对应正交向量是唯一的。表示d正交于其他向量张成空间的分量。X中端元S_P的丰度估计a_P为X在上的投影长度与长度的比值。则

OVP算法仅仅需要向量内积运算，不需要矩阵乘积和矩阵求逆运算，运算复杂度较低。

2.6OB算法

将端元向量组进行Gram-Schmidt正交化得到正交基组，计算解混方程组，解得正交基对应的特征向量。将待解混光谱向量向特征向量投影，计算该投影向量的长度与正交基的长度比，得到正交基代表端元的丰度估计。

将S₁,S₂,…,S_PGram-Schmidt正交化，得到正交基β₁,β₂,…,β_P

经推导，

由于β₁,β₂,…,β_P两两正交，可得

联立P个等式，可得解混线性方程组为

OB算法的思路是：将待解混光谱向量向正交基空间进行投影，由投影点在正交基上的解混系数获得待解混光谱向量的每个端元的对应丰度估计。

令d_i为正交基β_i的特征向量，即d_i为正交基β_i代表端元S_i在正交基空间的投影向量，则

则

在丰度估计运算中，有时需要对端元的个数稍作增减重复进行丰度估计，采用OB算法，式(10)和式(12)中绝大多数项不需要重新运算，仅需要计算新添加的项或者去掉多余项。

3、OB,OVP,SV,OSP,UCLS算法之间关系及算法复杂度分析

3.1OB,OVP,SV,OSP,UCLS算法一致性

OVP,SV,OSP,UCLS算法具有一致性。

OB算法丰度估计结果与UCLS算法的一致性证明：

公式(2)可写为min[(X-SA)^T(X-SA)]，其中S＝[S₁,S₂,…,S_P]。求解该问题等价于寻找一个使得满足从而使得(X-SA)^T(X-SA)达到最小。

式(3)为X在[S₁,S₂,…,S_P]上的最小二乘估计。由最小二乘估计和正交子空间投影的等价关系，为向量X在子空间[S₁,S₂,…,S_P]上的正交投影点。

β₁,β₂,…,β_P是子空间S₁,S₂,…,S_P的一组正交基，所以式(3)为在S₁,S₂,…,S_P上的坐标。

对于OB算法，解混方程组式(11)可写成

将全部P个等式整理后左右分别相加可得到

因此，OB算法的解混结果同样也是在S₁,S₂,…,S_P上的坐标。所以，OB算法的解混结果是混合像元的最小二乘逼近。

3.2算法的计算复杂度

假设2个标量间的1次相乘记为1次基本运算，则5种算法的计算复杂度如表1所示。其中，[]表示取整函数。由表1可以看出，算法的复杂度主要受波段数量和端元数量2个因素影响。一般情况下，高光谱图像中波段数量N很大，在端元数量P较小的时候，算法复杂度主要受波段数量影响。随着端元数量P的增加，在P达到并超过一定限值后，算法的复杂度主要受端元数量P影响。由表1可以看出，OB算法的计算复杂度显著低于OVP，SV，OSP，UCLS算法。

表1 5种算法的计算复杂度比较

4、实验分析

4.1模拟数据实验

图1为用于实验的Cuprite地区图像，选取图像中的明矾石Alunite(A)、水铵长石Buddingtonite(B)、玉髓Chalcedony(C)、高岭石Kaolinite(K)、蒙脱石Montmorillonite(M)共5种端元，并在图像中选取一块区域的光谱均值作为背景光谱b。图2为选择端元光谱。

TI为模拟图像，将b作为整幅高光谱图像的背景。5行像元块从上向下依次是A，B，C，K和M区域。5列像元块从左向右分别由5种矿物以不同丰度混合而成。即3×3纯像元，2×2纯像元，2×2主像元与其他4种像元各占50％混合，1×1主像元与背景各占50％混合，1×1主像元30％与背景70％混合，模拟数据图像如图3所示。TI1～TI3分别是无噪声、背景有噪声、像元块与背景像元都有噪声解混图像。模拟图像中加入的噪声为信噪比(SNR，signalto noise ratio)是20dB的Gauss白噪声。

对模拟数据图像进行OB，OVP，OSP和UCLS解混，解混结果如表2所示。可以看出，OB的解混结果与OVP，OSP，UCLS相同，验证了3.1的推导结论。

表2 4种算法的解混结果

4.2Hydice数据

去除水气吸收波段，Hydice数据保留169个波段，光谱空间分辨率为1.5m。选取4个不同类别地物光谱曲线(interferer,grass,tree,road)作为端元，如图4(a)所示。图4(b)中黄色区域的红色部分是目标区域，对应地面目标的15个像元，利用其中的5个目标光谱向量作为端元，组成9个端元。用这9个端元进行无约束线性丰度运算，UCLS，OSP，OVP和OB的计算结果完全相同。计算所需时间如表3所示。正如3.2分析，在端元数量较少的情况下，波段数量是运算时间主要的影响因素，所以4种算法运算时间相差不大。但是，当把已知端元增加到利用ATGP算法提取的144个端元时，则端元数量成为解混时间的主要影响因素，4种算法的解混时间急剧增加，如表4所示。

表3端元数量较小情况下4种算法的计算时间(S)

表4端元数量急剧增加时4种算法的计算时间(S)

从表3和表4可见，OB的运算时间最短，这与3.2的复杂度分析完全相符。

5、结论

UCLS，OSP，SV和OVP是常用的解混算法。UCLS，OSP和SV需要矩阵求逆或行列式运算，计算度复杂，且计算量比较大。OVP在运算复杂度上有所改进，但计算时间还是较长。本文在OVP基础上提出的OB算法，进一步减小了计算复杂度。理论上证明，该算法与UCLS，OSP，OVP解混结果一致。通过分析UCLS，OSP，OVP和OB算法，证明了OB算法的计算复杂度最小。同时，OB算法也可以逐个端元地进行丰度估计，以便更好地获取目标信息。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。

Claims (5)

1.一种基于正交基的高光谱图像丰度估计方法，其特征在于，所述方法的实现过程为：

利用Gram-Schmidt方法计算高光谱图像的端元向量组得到对应的正交基组，求解解混线性方程组，解得每个正交基的特征向量；将待解混光谱向量向特征向量投影得到投影向量，计算该投影向量的长度与正交基的长度比，得到正交基所代表端元的丰度估计；高光谱图像中的多个端元向量组成所述端元向量组。

2.根据权利要求1所述的一种基于正交基的高光谱图像丰度估计方法，其特征在于，所述方法的具体实现过程为：

将端元向量组S₁,S₂,…,S_P利用Gram-Schmidt方法进行正交化，得到正交基β₁,β₂,…,β_P

其中，i为端元序数，j为正交基序数；S_i为端元向量；

经推导，

X为待解混光谱向量，a₁,a₂,…,a_P为端元向量组S₁,S₂,…,S_P对应的丰度，即a_i为S_i的丰度；

由于β₁,β₂,…,β_P两两正交，可得

联立P个等式，可得解混线性方程组为

令d_i为正交基β_i的特征向量，即d_i为正交基β_i代表端元S_i在正交基空间的投影向量，则

则

表示正交基所代表端元的丰度估计。

3.根据权利要求1或2所述的一种基于正交基的高光谱图像丰度估计方法，其特征在于，所述正交基组与解混方程组相关，正交基是解混方程组的组成部分，求解解混线性方程组为了得到特征向量。

4.根据权利要求3所述的一种基于正交基的高光谱图像丰度估计方法，其特征在于，正交基所代表端元的丰度估计是将待解混光谱向量向正交基空间进行投影，由投影点在正交基上的解混系数获得待解混光谱向量的每个端元的对应丰度估计。

5.根据权利要求4所述的一种基于正交基的高光谱图像丰度估计方法，其特征在于，对端元的个数稍作增减重复进行丰度估计时，采用正交基所代表端元的丰度估计中，式(10)和式(12)中仅需要计算新添加的项或者去掉多余项。