CN105809105A - 基于端元约束非负矩阵分解的高光谱图像解混方法 - Google Patents

Abstract

基于端元约束非负矩阵分解的高光谱图像解混方法，属于高光谱图像分解方法领域。经典非负矩阵分解的目标函数具有非凸性，影响最优解的获得。一种基于端元约束非负矩阵分解的高光谱图像解混方法，采用光谱之间两两相关系数的绝对值之和作为相关性函数来衡量端元光谱相关性大小。端元光谱差别约束中引入自然对数函数减缓矩阵的迹运算的突变。采用投影梯度进行非负矩阵分解，综合图像分解误差和端元光谱的影响，得出目标函数。通过模拟数据实验和真实数据实验验证了算法的有效性。

Description

基于端元约束非负矩阵分解的高光谱图像解混方法

技术领域

本发明涉及一种基于端元约束非负矩阵分解的高光谱图像解混方法。

背景技术

由于光谱成像仪的空间分辨率限制和地物的复杂多样性，高光谱图像的某些像元中往往包含多种物质(即为端元)，这些包含其他端元的像元被称为混合像元。为了提高对真实地表覆盖的描述准确性，需要对混合像元进行分解，计算一种地物类型(端元)在该像元中所占的比例(即为丰度)。混合像元解混是高光谱图像定量分析的重要研究课题。1999年Lee和Seung在Nature杂志上提出了一种乘法迭代的非负矩阵分解方法，引起广泛关注。NMF算法具有强大的信息处理和问题求解能力，非负的限制符合许多实际问题的要求，该方法在工程领域得到广泛应用。经典NMF算法的目标函数具有明显的非凸性，存在局部极小值，影响最优解的获得。为了能使NMF算法应用于不同领域，通常需要根据不同应用场合的物理特性增加相应的约束条件。目前，已经提出平滑性限制(SC)、最小体积限制(MVC)、稀疏性限制等约束条件。

发明内容

本发明的目的是为了解决现有的经典非负矩阵分解的目标函数具有非凸性，影响最优解的获得，根据混合像元的特点，而提出一种以端元之间的相关性和不同端元之间的差别为约束条件，结合非负矩阵分解来进行混合像元分解，称为基于端元约束非负矩阵分解的高光谱图像解混方法(endmemberconstraintnonnegativematrixfactorization)，简称为EC-NMF。

一种基于端元约束非负矩阵分解的高光谱图像解混方法，所述方法通过以下步骤实现：

步骤一、设像元光谱矢量X，端元光谱矩阵S，N维向量的丰度矩阵A，随机噪声N，建立线性光谱混合模型：

X＝SA+N(1)

其中，所述端元是指光谱成像仪所呈的高光谱图像的像元中包含的多种物质，这些包含端元的像元被称为混合像元；

步骤二、将光谱之间两两相关系数的绝对值之和作为相关性函数，来衡量端元光谱相关性大小；

步骤三、增加引入自然对数函数的端元光谱差别约束，使端元光谱之间的差别达到最大；

步骤四、综合衡量图像分解误差和步骤二、步骤三获得的端元光谱约束的影响，建立目标函数。

本发明的有益效果为：

本发明提出的端元光谱相关性最小化约束和端元光谱差别最大化约束的非负矩阵分解算法，采用光谱之间两两相关系数的绝对值之和作为相关性函数来衡量端元光谱相关性大小。初次对非负矩阵分解进行约束后，还不能获得理想解混精度，还需要增加端元光谱差别约束，进一步提高解混精度，由于端元光谱应具有较大的独立性，良好的端元光谱应当使得端元光谱之间的差别达到最大，所以通过端元光谱自相关矩阵的迹的倒数表达端元光谱之间的差别，这种在端元光谱差别约束中引入自然对数函数的方法，能减缓矩阵的迹运算的突变。最后，采用投影梯度进行非负矩阵分解，综合图像分解误差和端元光谱的影响，得出目标函数。

通过模拟高光谱图像和真实高光谱图像验证了算法的有效性，通过值、值以及光谱角进行对比时，本发明的高光谱图像解混方法比其他方法图像解混精度提高10-15％；

附图说明

图1为本发明的流程图；

图2为本发明涉及的不同噪声强度时算法性能的比较；

图3为本发明涉及的不同像元个数时算法性能的比较；

图4为本发明真实实验数据中涉及的Cuprite地区172波段的图像；

图5为本发明涉及的Cuprite区域的AVIRIS数据含有的矿物基于EC-NMF算法的解混结果；其中，图5(a)为明矾石、图5(b)为水铵长石、图5(c)为玉髓、图5(d)为黄钾铁矾、图5(e)为高岭石#1、图5(f)为高岭石#2、图5(g)为高岭石#3、图5(h)为蒙脱石、图5(i)为白云母、图5(j)为绿脱石、图5(k)为榍石。

具体实施方式

具体实施方式一：

本实施方式的基于端元约束非负矩阵分解的高光谱图像解混方法，结合图1所示，所述方法通过以下步骤实现：

步骤一、设像元光谱矢量X，端元光谱矩阵S，N维向量的丰度矩阵A，随机噪声N，建立线性光谱混合模型：

X＝SA+N(1)

其中，所述端元是指光谱成像仪所呈的高光谱图像的像元中包含的多种物质，这些包含端元的像元被称为混合像元；

步骤二、将光谱之间两两相关系数的绝对值之和作为相关性函数，来衡量端元光谱相关性大小；

步骤三、增加引入自然对数函数的端元光谱差别约束，使端元光谱之间的差别达到最大；

步骤四、综合衡量图像分解误差和步骤二、步骤三获得的端元光谱约束的影响，建立目标函数。

具体实施方式二：

与具体实施方式一不同的是，本实施方式的基于端元约束非负矩阵分解的高光谱图像解混方法，步骤一所述线性光谱混合模型X＝SA+N中，端元光谱矩阵S＝[s₁,s₂,…,s_N]，端元光谱矩阵S中的元素s_i表示端元向量，i∈[1,N]；N维向量的丰度矩阵A＝[a₁,a₂,…,a_N]^T，N维矢量的丰度矩阵A中各分量元素表示对应端元的丰度，且

a_i≥0(2)

$\underset{i=1}{\overset{N}{\Σ}}{a}_i=1,i\∈\[1,N\]---\left(3\right)$

其中，所述丰度是指一种端元在该像元中所占的比例。

具体实施方式三：

与具体实施方式一或二不同的是，本实施方式的基于端元约束非负矩阵分解的高光谱图像解混方法，步骤二所述将光谱之间两两相关系数的绝对值之和作为相关性函数，来衡量端元光谱相关性大小的过程为，

步骤二一、进行基于端元约束的非负矩阵分解：

利用NMF算法，通过最小化欧式距离目标函数，在已知X的情况下解得S和A的最优解，

$f(S,A)=\frac{1}{2}\left|\right|X-SA|{|}^2---\left(4\right)$

迭代公式为：

S←S-β₁(SA-X)A^T(5)

A←A-β₂S^T(SA-X)(6)

基于NMF的光谱解混算法不需要判定是否存在纯像元，在提取端元的同时获取对应端元的丰度；

式中，f表示X和SA之间的欧式距离；β₁表示权值，β₂表示权值；

步骤二二、端元光谱相关性约束：

将经步骤一中公式(3)约束后的端元光谱进行归一化处理，得两条光谱向量的相关系数函数为：

$\mathrm{\ρ}(S_i,S_j)=\underset{k=1}{\overset{L}{\Σ}}{S}_{ki}{S}_{kj}---\left(7\right)$

公式(7)中，ρ表示相关系数函数；S_i、S_j表示端元光谱向量；k表示波段序数；S_ki、S_kj表示反射率；i表示N个端元中的第i个；j表示N个端元中的第j个；

公式(7)中，当2个端元向量的变化趋势一致时，相关系数函数为正值；当2个端元向量的变化趋势相反时，则相关系数函数为负值；当2个向量不相关时，相关系数函数为0；

步骤二三、定义由每2条光谱相关系数函数的绝对值的和组成光谱相关性函数：

$\mathrm{\μ}\left(S\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=i+1}{\overset{N}{\Σ}}\left|\mathrm{\ρ}(S_i,S_j)\right|---\left(8\right)$

来衡量N条端元光谱的整体相关性；

其中，|·|为绝对值函数；当每2条光谱的相关性减小时，函数|ρ(S_i,S_j)|的值减小并趋近于0，故μ(S)的值也随之减小，反之亦然；

步骤二四、采用公式(9)表示的光谱相关性函数作为约束条件，

$J_1\left(S\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=i+1}{\overset{N}{\Σ}}\left|\mathrm{\ρ}(S_i,S_j)\right|---\left(9\right)$

步骤二五、对公式(9)求导数将该约束最小化，整理得到公式(10)表示的最小化的光谱相关性函数约束条件，

${\▿}_S{J}_1\left(S\right)=\underset{i=1,i\≠j}{\overset{N}{\Σ}}sign\left(\mathrm{\ρ}\right(S_i,S_j\left)\right)S_j^T{S}_i---\left(10\right);$

式中，▽_S表示求导数计算；sign(·)表示符号函数。

具体实施方式四：

与具体实施方式三不同的是，本实施方式的基于端元约束非负矩阵分解的高光谱图像解混方法，步骤三所述增加引入自然对数函数的端元光谱差别约束，使端元光谱之间的差别达到最大的过程为，

步骤三一、增加端元光谱差别约束，进一步提高解混精度：

通过公式(11)表示的端元光谱自相关矩阵的迹的倒数表达端元光谱之间的差别：

J₂(S)＝(Tr(A^TA))^-1(11)

式中，Tr表示矩阵的迹运算；J₂(S)值越小，端元光谱之间的差别越大；A^TA为对称正定矩阵，A^TA的全部特征值λ_i≥0，则∑λ_i＞0，则Tr(A^TA)＝∑λ_i，则J₂(S)＞0；

步骤三二、由于公式(11)的函数值受矩阵变化影响较大，引入ln自然对数函数减缓Tr函数值的突变：

采用公式(12)表示的优化的函数项，在不改变Tr函数单调性的同时减弱公式(11)对目标函数稳定性的影响，

J₂(S)＝ln((Tr(A^TA))^-1)＝-ln(Tr(A^TA))(12)

步骤三三、对公式(12)求导数将该约束最小化，整理得到公式(13)表示的最小化的约束条件：

${\▿}_S{J}_2\left(S\right)=-\frac{A}{Tr\left(A^{T}A\right)}---\left(13\right).$

具体实施方式五：

与具体实施方式一、二或四不同的是，本实施方式的基于端元约束非负矩阵分解的高光谱图像解混方法，步骤四所述建立目标函数的过程为，

步骤四一、通过公式(1)的NMF线性混合模型，使公式(2)代表的端元的丰度的约束得到满足；采用公式(14)的形式，使式(3)的约束得到满足；即全约束最小二乘算法(FCLS)，

$\stackrel{\‾}{X}=\left[\begin{array}{c}X\\ \mathrm{\δ}{1}_N^T\end{array}\right],\stackrel{\‾}{S}=\left[\begin{array}{c}S\\ \mathrm{\δ}{1}_N^T\end{array}\right]---\left(14\right)$

式中，1表示其分量都为1的向量，即δ为权值；迭代计算时，用代替X，用代替S；

EC-NMF的投影梯度迭代公式为：

$\stackrel{\‾}{S}\←P_{\mathrm{\Ω}}\left(x\right)(\stackrel{\‾}{S}-{\mathrm{\β}}_1(\left(\stackrel{\‾}{S}A-\stackrel{\‾}{X}\right)A^T-{\mathrm{\α}}_1{\▿}_{\stackrel{\‾}{S}}{J}_1\left(\stackrel{\‾}{S}\right)-{\mathrm{\α}}_2{\▿}_{\stackrel{\‾}{S}}{J}_2\left(\stackrel{\‾}{S}\right)\left)\right)---\left(15\right)$

$A\←P_{\mathrm{\Ω}}\left(x\right)(A-{\mathrm{\β}}_2{\stackrel{\‾}{S}}^T(\stackrel{\‾}{S}A-\stackrel{\‾}{X}\left)\right)---\left(16\right)$

投影函数为：P_Ω(x)＝max(0,x)，P_Ω表示最大值函数，在进行式(15)的迭代时，x表示式在进行式(16)的迭代时，x表示式 $A-{\mathrm{\β}}_2{\stackrel{\‾}{S}}^T(\stackrel{\‾}{S}A-\stackrel{\‾}{X});$

步骤四二、建立EC-NMF的目标函数为：

J(S,A)＝f(S,A)-α₁J₁(S,A)-α₂J₂(S,A)(17)

式中，α₁、α₂表示权值。

仿真实验：

(一)性能指标

采用光谱角距离(SpectralAngleDistance，SAD)和均方根误差(RootMeanSquareError，RMSE)这两个普遍应用的指标来衡量解混效果。他们常被用来计算光谱和丰度解混估计值与真实值的近似程度。对于第i个端元，SAD定义为：

${\mathrm{SAD}}_i={\cos }^{-1}\left(\frac{S_i^T{\hat{S}}_i}{\left|\right|S_i^T\left|\right|\left|\right|{\hat{S}}_i\left|\right|}\right)---\left(18\right)$

其中，S为真实端元光谱，为其估计值。

RMSE定义为：

${\mathrm{RMSE}}_i={\left(\frac{1}{N}\underset{j=1}{\overset{N}{\Σ}}{\left(A_{ij}-{\hat{A}}_{ij}\right)}^2\right)}^{\frac{1}{2}}$

其中A_i为端元对应真实丰度，为其估计值。

(二)模拟数据实验

从USGS矿物光谱库中选取5种线性独立的端元光谱(明矾石、水铵长石、方解石、高岭石、白云母)，按Dirichlet分布混合，端元丰度之和进行归一化，加上不同的白噪声，形成模拟实验数据。

抗噪性能实验：对EC-NMF算法与MVCNMF、SCNMF算法不同噪声情况下解混性能进行比较。信噪比在∞(不含噪声)、30、20、10dB是进行实验，表示数学期望。

如图2所示的不同噪声强度时算法性能的比较，其中，图2(a)纵坐标表示数学期望，图2(b)纵坐标表示数学期望。

不同像元个数解混性能实验：像元个数不同时，EC-NMF算法与MVCNMF、SCNMF算法解混性能比较。像元数目分别为1600、3600、6400、10000，信噪比固定为30dB。

如图3所示的不同像元个数时算法性能的比较，图3(a)纵坐标表示数学期望，图3(b)纵坐标表示数学期望。

观察实验结果，EC-NMF算法解混性能要优于MVCNMF、SCNMF算法。

(三)真实数据实验

采用美国JPL实验室的AVIRIS获得的内华达州Cuprite高光谱图像，图4为第172波段图像。去除信噪比低和水蒸气吸收波段余下49个有效波段。

图5(a)到图5(k)为Cuprite区域的AVIRIS数据含有的矿物基于EC-NMF算法的解混结果，表1为三种算法的结果，整体来说，EC-NMF性能优于另两种方法。

表1Cuprite数据的光谱角比较

(四)总结

针对端元独立性的特点，提出一种基于端元约束的非负矩阵分解的混合像元分解方法。该方法适用于无纯像元的情况。通过模拟数据与真实数据实验分析，表明本文提出的方法能够提高混合像元分解的精度。投影梯度算法的收敛速度较慢，需采用更有效的优化算法来提高解混的效率。可以考虑采用其他形式的约束条件来提高混合像元的分解精度。

Claims (5)

1.一种基于端元约束非负矩阵分解的高光谱图像解混方法，其特征在于：所述方法通过以下步骤实现：

步骤一、设像元光谱向量X，端元光谱矩阵S，N维向量的丰度矩阵A，随机噪声N，建立线性光谱混合模型：

X＝SA+N(1)

其中，所述端元是指光谱成像仪所呈的高光谱图像的像元中包含的多种物质，这些包含端元的像元被称为混合像元；

步骤二、将光谱之间两两相关系数的绝对值之和作为相关性函数，来衡量端元光谱相关性大小；

步骤三、增加引入自然对数函数的端元光谱差别约束，使端元光谱之间的差别达到最大；

步骤四、综合衡量图像分解误差和步骤二、步骤三获得的端元光谱约束的影响，建立目标函数。

2.根据权利要求1所述的基于端元约束非负矩阵分解的高光谱图像解混方法，其特征在于：步骤一所述线性光谱混合模型X＝SA+N中，端元光谱矩阵S＝[s₁,s₂,‖,s_N]，端元光谱矩阵S中的元素s_i表示端元向量，i∈[1,N]；N维向量的丰度矩阵A＝[a₁,a₂,‖,a_N]^T，N维向量的丰度矩阵A中各分量元素表示对应端元的丰度，且

a_i≥0(2)

$\underset{i=1}{\overset{N}{\Σ}}{a}_i=1,i\∈\[1,N\]---\left(3\right)$

其中，所述丰度是指一种端元在该像元中所占的比例。

3.根据权利要求1或2所述的基于端元约束非负矩阵分解的高光谱图像解混方法，其特征在于：步骤二所述将光谱之间两两相关系数的绝对值之和作为相关性函数，来衡量端元光谱相关性大小的过程为，

步骤二一、进行基于端元约束的非负矩阵分解：

利用NMF算法，通过最小化欧式距离目标函数，在已知X的情况下解得S和A的最优解，

$f(S,A)=\frac{1}{2}\left|\right|X-SA|{|}^2---\left(4\right)$

迭代公式为：

S←S-β₁(SA-X)A^T(5)

A←A-β₂S^T(SA-X)(6)

基于NMF的光谱解混算法不需要判定是否存在纯像元，在提取端元的同时获取对应端元的丰度；

式中，f表示X和SA之间的欧式距离；β₁表示权值，β₂表示权值；

步骤二二、端元光谱相关性约束：

将经步骤一中公式(3)约束后的端元光谱进行归一化处理，得两条光谱向量的相关系数函数为：

$\mathrm{\ρ}(S_i,S_j)=\underset{k=1}{\overset{L}{\Σ}}{S}_{ki}{S}_{kj}---\left(7\right)$

公式(7)中，ρ表示相关系数函数；S_i、S_j表示端元光谱向量；k表示波段序数；S_ki、S_kj表示反射率；i表示N个端元中的第i个；j表示N个端元中的第j个；

公式(7)中，当2个端元向量的变化趋势一致时，相关系数函数为正值；当2个端元向量的变化趋势相反时，则相关系数函数为负值；当2个向量不相关时，相关系数函数为0；

步骤二三、定义由每2条光谱相关系数函数的绝对值的和组成光谱相关性函数：

$\mathrm{\μ}\left(S\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=i+1}{\overset{N}{\Σ}}\left|\mathrm{\ρ}(S_i,S_j)\right|---\left(8\right)$

来衡量N条端元光谱的整体相关性；

其中，|·|为绝对值函数；当每2条光谱的相关性减小时，函数|ρ(S_i,S_j)|的值减小并趋近于0，故μ(S)的值也随之减小，反之亦然；

步骤二四、采用公式(9)表示的光谱相关性函数作为约束条件，

$J_1\left(S\right)=\underset{i=1}{\overset{N}{\Σ}}\underset{j=i+1}{\overset{N}{\Σ}}\left|\mathrm{\ρ}(S_i,S_j)\right|---\left(9\right)$

步骤二五、对公式(9)求导数将该约束最小化，整理得到公式(10)表示的最小化的光谱相关性函数约束条件，

${\▿}_S{J}_1\left(S\right)=\underset{i=1,i\≠j}{\overset{N}{\Σ}}sign\left(\mathrm{\ρ}\right(S_i,S_j\left)\right)S_j^T{S}_i---\left(10\right);$

式中，▽_S表示求导数计算；sign(·)表示符号函数。

4.根据权利要求3所述的基于端元约束非负矩阵分解的高光谱图像解混方法，其特征在于：步骤三所述增加引入自然对数函数的端元光谱差别约束，使端元光谱之间的差别达到最大的过程为，

步骤三一、增加端元光谱差别约束，进一步提高解混精度：

通过公式(11)表示的端元光谱自相关矩阵的迹的倒数表达端元光谱之间的差别：

J₂(S)＝(Tr(A^TA))^-1(11)

式中，Tr表示矩阵的迹运算；J₂(S)值越小，端元光谱之间的差别越大；A^TA为对称正定矩阵，A^TA的全部特征值λ_i≥0，则Σλ_i＞0，则Tr(A^TA)＝Σλ_i，则J₂(S)＞0；

步骤三二、引入ln自然对数函数减缓Tr函数值的突变：

采用公式(12)表示的优化的函数项，在不改变Tr函数单调性的同时减弱公式(11)对目标函数稳定性的影响，

J₂(S)＝ln((Tr(A^TA))^-1)＝-ln(Tr(A^TA))(12)

步骤三三、对公式(12)求导数将该约束最小化，整理得到公式(13)表示的最小化的约束条件：

${\▿}_S{J}_2\left(S\right)=-\frac{A}{Tr\left(A^{T}A\right)}---\left(13\right).$

5.根据权利要求1、2或4所述的基于端元约束非负矩阵分解的高光谱图像解混方法，其特征在于：步骤四所述建立目标函数的过程为，

步骤四一、通过公式(1)的NMF线性混合模型，使公式(2)代表的端元的丰度的约束得到满足；采用公式(14)的形式，使式(3)的约束得到满足，即全约束最小二乘算法，

$\stackrel{\‾}{X}=\left[\begin{array}{c}X\\ \mathrm{\δ}{1}_N^T\end{array}\right],\stackrel{\‾}{S}=\left[\begin{array}{c}S\\ \mathrm{\δ}{1}_N^T\end{array}\right]---\left(14\right)$

式中，1表示其分量都为1的向量，即δ为权值；迭代计算时，用代替X，用代替S；

EC-NMF的投影梯度迭代公式为：

$\stackrel{\‾}{S}\←P_{\mathrm{\Ω}}\left(x\right)(\stackrel{\‾}{S}-{\mathrm{\β}}_1(\left(\stackrel{\‾}{S}A-\stackrel{\‾}{X}\right)A^T-{\mathrm{\α}}_1{\▿}_{\stackrel{\‾}{S}}{J}_1\left(\stackrel{\‾}{S}\right)-{\mathrm{\α}}_2{\▿}_{\stackrel{\‾}{S}}{J}_2\left(\stackrel{\‾}{S}\right)\left)\right)---\left(15\right)$

$A\←P_{\mathrm{\Ω}}\left(x\right)(A-{\mathrm{\β}}_2{\stackrel{\‾}{S}}^T(\stackrel{\‾}{S}A-\stackrel{\‾}{X}\left)\right)---\left(16\right)$

投影函数为：P_Ω(x)＝max(0,x)，P_Ω表示最大值函数，在进行式(15)的迭代时，x表示式 $\stackrel{\‾}{S}-{\mathrm{\β}}_1\left(\right(\stackrel{\‾}{S}A-\stackrel{\‾}{X})A^T-{\mathrm{\α}}_1{\▿}_{\stackrel{\‾}{S}}{J}_1(\stackrel{\‾}{S})-{\mathrm{\α}}_2{\▿}_{\stackrel{\‾}{S}}{J}_2(\stackrel{\‾}{S}\left)\right),$ 在进行式(16)的迭代时，x表示式 $A-{\mathrm{\β}}_2{\stackrel{\‾}{S}}^T(\stackrel{\‾}{S}A-\stackrel{\‾}{X});$

步骤四二、建立EC-NMF的目标函数为：

J(S,A)＝f(S,A)-α₁J₁(S,A)-α₂J₂(S,A)(17)

式中，α₁、α₂表示权值。