CN105513097A - 一种基于最小体积与优化约束条件的高光谱解混方法 - Google Patents

Abstract

本发明公开了一种基于最小体积与优化约束条件的高光谱解混方法，其步骤如下：步骤一：数据加载与预处理。步骤二：筛选图像采样点，构造优化的约束条件，寻找满足初始条件的数据。步骤三：将非负非线性规划问题转换为线性规划问题，对线性规划问题的目标函数结合优化后的约束条件进行求解，计算中间变量矩阵H_new、g_new。步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回步骤三，继续更新中间变量矩阵H_new、g_new。步骤五：解出满足非负性要求的端元矩阵，并计算丰度系数，完成图像的解混。本方法解决基于最小体积的高光谱解混算法中，约束条件过多所导致的数据存储空间大，运算时间长，精确性难以提高的问题。

Description

一种基于最小体积与优化约束条件的高光谱解混方法

技术领域

本发明涉及一种高光谱数据解混方法，具体涉及一种基于几何学最小体积算法解混的方法。

背景技术

高光谱遥感能利用很窄的电磁波波段从感兴趣的物体上获取有关数据。高光谱遥感是在电磁波谱的可见光、近红外、中红外和热红外波段范围内，获取许多非常窄的连续光谱数据的技术，目前先进的成像光谱仪可以收集到成百上千个非常窄的光谱波段信息。

利用高光谱解混技术能够确定高光谱图像中的端元数、端元光谱曲线以及每个像素点的丰度值，进而能够确定某种物质在图像观察区域的分布。在解混方法中，高光谱的混合模型分为两类，分别为线性混合模型以及非线性混合模型。线性模型假设物质成份在空间布局上是离散混合，忽略在不同类型物质之间的多重散射量，即每个像元点的光谱幅度仅是这个像元中存在的各物质成份光谱信号，以相对贡献量多少的线性叠加。对于高光谱线性混合模型来说，解混方法可以分为三大类，分别为：1)基于凸面几何学方法的线性解混模型；2)基于数据统计的方法，代表算法有ICA等，其核心思想在于假设原信号在体积上是独立分布的；3)基于稀疏模型的解混方法，代表算法有BP、BPDN等。

基于几何学的解混算法又能够被分为两类：一类是基于纯像元假设的解混算法，另一类是基于最小体积法的解混算法。基于纯像元假设的解混算法需要保证图像中的每个不同的端元至少完全占据一个像素，否则会出现端元错误估计的情况。而基于最小体积的解混算法即使在缺失纯像元的情况下也能够较高地估计出端元，但缺点是其优化目标函数是非线性且非凸的，很多情况下，只能得到近似的局部最优解，同时其约束条件数量非常多，使其算法效率没有基于纯像元假设的解混算法高。

发明内容

为了解决基于最小体积的高光谱解混算法中，约束条件过多所导致的数据存储空间大，运算时间长，精确性难以提高的问题，本发明提供了一种基于最小体积与优化约束条件的高光谱解混方法，通过优化约束条件改进基于最小体积法的高光谱解混方法，得出更准确的解混效果。

本发明的目的是通过以下技术方案实现的：

一种基于最小体积与优化约束条件的高光谱解混方法，对高光谱图像进行线性解混以求取端元矩阵时，考虑初值对求解端元矩阵的影响，并将求解体积最小问题的非线性问题转化为线性规划问题，近似求解得到非线性优化的局部最小解；并考虑求解端元的非负物理特性，从而提高求取端元的准确率。具体包括如下步骤：

步骤一：数据加载与预处理。

步骤二：筛选图像采样点，构造优化的约束条件，寻找满足初始条件的数据。

步骤三：将非负非线性规划问题转换为线性规划问题，对线性规划问题的目标函数结合优化后的约束条件进行求解，计算中间变量矩阵H_new、g_new。

步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回步骤三，继续更新中间变量矩阵H_new、g_new。

步骤五：解出满足非负性要求的端元矩阵，并计算丰度系数，完成图像的解混。

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

1)通过筛选数据采样点，同时利用某一采样点必然能包含在由其他顶点构成的单形体内的性质，有效缩减了约束条件的个数使得计算开销变小，同时解混效率有较大提高。

2)不论是否存在纯像元，本方法都能够经过解算，得到高光谱图像的端元，进而得到丰度系数，完成高光谱图像的解混。

附图说明

图1为本发明的流程图；

图2为数据点集、真实端元集、估计端元集在2(p-1)维上的投影显示；

图3为真实端元光谱图与估计端元光谱图在224个波段上的显示。

具体实施方式

下面结合附图对本发明的技术方案作进一步的说明，但并不局限于此，凡是对本发明技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围，均应涵盖在本发明的保护范围中。

具体实施方式一：本实施方式提供了一种基于最小体积与优化约束条件的高光谱解混方法，假设光谱图像数据L是波段数，N是图像的像素点数量。若在光谱图像中有p个端元，首先对光谱信息进行降维处理，将数据从L维映射到p-1维；接着构造一个单形体，该单形体能够包含所有降维后的数据，通过循环求解线性规划子问题，不断更新单形体顶点矩阵，最终使单形体的体积收敛至稳定，最后通过单形体的顶点矩阵得到图像的端元矩阵，进而求出图像的丰度系数，完成高光谱图像的解混任务。

如图1所示，共分为五个步骤，具体步骤如下：

步骤一：数据加载与预处理。

1)加载数据矩阵其中L是波段数，N是图像的像素点数量；

2)对数据矩阵进行降维。已知光谱图像中，端元个数为p，将数据矩阵Y降维至p-1维数据矩阵具体降维方法为：

其中，U＝[y₁-d,y₂-d,...,y_N-d]，数据矩阵Y表示为列向量的形式：Y＝[y₁,y₂,...,y_N]；Q＝[q₁(UU^T),q₂(UU^T),...,q_p-1(UU^T)]，q_i(UU^T),i＝1,2,...,p-1表示矩阵UU^T第i大的特征值所对应的特征向量。

步骤二：采用优化后的约束条件，寻找满足初始条件的数据。

1)在数据矩阵中得到p个列向量：

$M_0=\[{\tilde{y}}_1^{\′},{\tilde{y}}_2^{\′},...,{\tilde{y}}_p^{\′}\]---\left(2\right).$

式中，向量选择满足：由这p个列向量作为顶点构成的单形体体积比其他任意p个列向量作为顶点构成的单形体体积大的条件。

2)在数据集中寻找数据点集，数据集中任意一点满足不在由M₀中的列向量作为顶点构成的单形体中的条件，最后得到满足要求的数据点集

具体寻找方法如下，首先令：

$B=\[{\tilde{y}}_1^{\′}-{\tilde{y}}_p^{\′},...,{\tilde{y}}_{p-1}^{\′}-{\tilde{y}}_p^{\′}\]---\left(3\right);$

$s=B^{-1}{\tilde{y}}_i-B^{-1}{\tilde{y}}_p^{\′}---\left(4\right);$

$\tilde{s}=\left[\begin{array}{c}s\\ 1-\underset{i=1}{\overset{p-1}{\Σ}}{s}_i\end{array}\right]---\left(5\right).$

其中， ${\tilde{y}}_i\∈\tilde{Y},$ s_i为s中的第i行元素。

对所有(N为像素点个数)求取其对应的向量若向量不满足的条件，则将加入数据集合中。假设最终有m(显然p≤m≤N，且一般情况下m＜＜N)个不满足上述条件的则最后得到的数据集合

3)优化的线性不等式约束条件为式(6)与式(7)：

${\mathrm{Hy}}_b\[i\]-g\≥0,\∀y_b\[i\]\∈{\tilde{Y}}_{block}---\left(6\right);$

$1_{p-1}^{T}H-1_{p-1}^{T}gq\≤q---\left(7\right);$

$q=1_m^T{\tilde{Y}}_{block}^T{\left({\tilde{Y}}_{block}{\tilde{Y}}_{block}^T\right)}^{-1}---\left(8\right).$

其中，q∈R^1×(p-1)，H与g为要优化的变量矩阵。式(6)与式(7)将不等式约束个数降为了(p-1)×(m+1)个，不等式约束的个数被大大减小，从而大大提高了求解的速度。

本步骤通过求解线性规划问题的方法，求解出一个满足线性不等式约束式(6)与式(7)关于变量矩阵H与g的一组可行解矩阵H₀与g₀，作为后续步骤的初始值。

本步骤中的＜、≤、≥分别表示所有元素都＜、≤、≥。

步骤三：对线性规划问题的目标函数结合优化后的约束条件进行求解，计算中间变量矩阵H_new、g_new。迭代计算中，变量矩阵H与g的阶段性最优变量值存储在矩阵H_new、g_new中。

构造优化函数，通过代数余子式的展开，将非线性的优化问题改造成为线性规划子问题，通过周期循环迭代的方式，更新得到新的矩阵H_new、g_new，直到满足停止迭代的条件。

1)构造非线性目标函数：

$f(H,g)=\underset{H,g}{max}\left|\det \left(H\right)\right|---\left(9\right).$

2)将非线性目标函数转化为线性的目标函数：

首先将H的行列式det(H)通过代数余子式的方式展开：

$\det \left(H\right)=\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right)---\left(10\right).$

式中h_ij是矩阵H的第i行第j列的元素，H_ij是将矩阵H移除第i行元素与第j列元素后得到的子阵。进而目标函数转化为：

$f(H,g)=\underset{k_i^T,g_i}{max}\left|\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right)\right|---\left(11\right).$

由于目标优化函数|det(H)|仍然是非线性的，因此将其分解为下面两个线性目标函数：

$p^{*}(h_i,g_i)=\underset{h_i^T,g_i}{max}\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right)---\left(12\right);$

$q^{*}(h_i,g_i)=\underset{h_i^T,g_i}{min}\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+1}{h}_{ij}\det \left(H_{ij}\right)---\left(13\right).$

3)将线性目标函数与优化后的线性约束条件结合，构造如下双线性规划子问题：

$p^{*}=\underset{h_i^T,g_i}{\max }\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{k}_{ij}\det \left(H_{ij}\right)---\left(14\right);$

$\begin{array}{cc}s.t.& h_i^T{y}_b\[i\]-g_i\≥0,\∀y_b\[i\]\∈{\tilde{Y}}_{block},1_{p-1}^{T}H-1_{p-1}^{T}gq\≤q\end{array}$

$q^{*}=\underset{h_i^T,g_i}{min}\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+1}{k}_{ij}\det \left(H_{ij}\right)---\left(15\right).$

$\begin{array}{cc}s.t.& h_i^T{y}_b\[i\]-g_i\≥0,\∀y_b\[i\]\∈{\tilde{Y}}_{block},1_{p-1}^{T}H-1_{p-1}^{T}gq\≤q\end{array}$

其中，h_i为矩阵H第i行的行向量，g_i为向量g的第i个元素。

4)求解上述两个线性规划子问题，根据p^*与q^*的值，选择性地更新对应的矩阵元素h_i与g_i，直到将矩阵H与g中所有的元素都更新一遍。

具体更新h_i与g_i的规则为：若求解线性规划问题式(14)得到的解为求解线性规划问题式(15)得到的解为h _i、g _i，若|p^*|＞|q^*|，则将h_i与g_i更新为与否则，将h_i与g_i更新为h _i与g _i。

在这一步骤中，逐行更新H与g的值：令i＝1,...,p-1，依次得到h_i与g_i的值，最终得到矩阵H与g中所有元素都更新后的新矩阵H_new与g_new。

步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回步骤三，继续更新中间变量矩阵H_new，g_new。

1)设置变化率检测终止条件如下：

{|det(H_new)|-|det(H₀)|}/|det(H)|＜ε(16)。

式中，H₀代表在步骤三开始迭代之前，矩阵H的初值。

2)判断中间变量矩阵H_new是否满足条件式(16)，如果满足终止条件，则令H^*＝H_new，g^*＝g_new，H^*与g^*作为矩阵H与g的最终更新结果。否则，令H＝H_new，H₀＝H_new，g＝g_new，转入步骤三的第4)步。

步骤五：还原端元矩阵，计算丰度系数，完成图像的解混。

根据H^*与g^*，还原端元矩阵，并计算丰度系数，完成图像的解混，端元矩阵应满足非负性的要求。

1)令α_p＝(H^*)^-1g^*，α＝[α₁,α₂,...,α_p]，则端元矩阵A为：

A＝Qα+d(17)。

2)若矩阵A是非负矩阵，计算丰度系数矩阵：

$X={\left[\begin{array}{cc}{(H^{*}\tilde{Y}-g^{*}\×I_N^T)}^T& 1_N-{(H^{*}\tilde{Y}-g^{*}\×1_N^T)}^T{1}_{p-1}\end{array}\right]}^T---\left(18\right).$

式中为行向量，且向量中每个元素值都为1，向量中元素个数为N；1_N为列向量，且向量中每个元素值都为1，向量中元素个数为N；1_p-1为列向量，且向量中每个元素值都为1，向量中元素个数为(p-1)；x[i]为第i个像素点丰度系数向量，已由步骤一计算得出。

具体实施方式二：本实施方式中，用于解混实验的高光谱图像遥感数据集由标准光谱库数据混合得到，混合光谱数据中端元个数p＝3，每个像素的丰度值满足和为1并且非负的条件。数据的信噪比为40dB。数据中的端元光谱曲线在USGS光谱库中选择。数据大小为224×5000，224为波段数，5000为像素点数。该仿真数据集中不存在纯像元，像素的丰度系数均小于0.8。

执行步骤一：数据加载与预处理。

加载存储原始数据的矩阵已知数据中端元数p＝3。

对数据矩阵image0进行降维处理，由于端元数为3，因此，通过式(1)得到降至2维的数据矩阵，保存在image1中，image1代表降维后的数据矩阵。

执行步骤二：寻找满足初始条件数据。

将降维后的数据矩阵image1作为输入矩阵，根据具体实施方式一种步骤二第1)步的方法，在数据矩阵image1中选择3个列向量，将结果存入矩阵M₀，

对降维后的数据矩阵image1中所有的列向量进行如下判断：首先，根据式(3)、式(4)和式(5)得到矩阵将结果存入中间变量矩阵S_0中；然后，再判断S_0中所有数值是否都大于0，如果不是，则将存入数据筛选矩阵image2中；最后，得到的数据矩筛选阵m≤5000，并且在很多情况下m＜＜5000，对于本实施方式m＝2114。

通过式(8)得到数据矩阵q。以中间变量矩阵H与g的元素作为变量，根据线性不等式(6)和式(7)，得到一组可行解，将可行解存入中间变量初始化矩阵H_0与g_0中，作为后续迭代初始值。

执行步骤三：循环求解双线性规划子问题，选择性更新最优值。

若是由步骤二转入步骤三，则中间变量迭代初值矩阵H_origin与g_origin有H_origin＝H_0，g_origin＝g_0。

令中间变量矩阵H＝H_origin，g＝g_origin，同时采用如下方法更新数据矩阵H与g的第i行的所有元素：

首先采用内点对偶法求解线性规划问题式(14)，得到矩阵行向量hp_new，元素gp_new为此线性规划问题的最优可行解，以及在这组可行解下的目标函数值p_val；然后采用同样的方法求解线性规划问题式(15)，得到最优可行解矩阵向量hq_new，元素gq_new，以及在这组可行解下的目标函数值q_val。如果满足|p_val|＞|q_val|的条件，则用hp_new、gp_new分别代替数据矩阵H与g的第i行所有元素；否则用hq_new、gq_new分别代替数据矩阵H与g的第i行所有元素。

采用上述更新H与g的第i行所有元素的方式，将中间变量矩阵H与g从第1行到第p-1(即3-1＝2)行的所有元素都替换一遍。进而完成H与g的全局更新。

执行步骤四：根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回步骤三，继续更新中间变量矩阵H与g。

首先，设置比较阈值th＝ε＝0.05；其次，计算数值e＝(|det(H)|-|det(H_origin)|)/|det(H_origin)|；最后，若e＜ε，则令中间变量终值矩阵H_Final＝H，g_Final＝g；否则令中间变量迭代初值矩阵H_origin＝H，g_origin＝g，并转回执行步骤三。

执行步骤五：还原端元矩阵，计算丰度系数，完成图像的解混。

计算矩阵alpha_p＝(H_Final)^-1×g_Final；计算矩阵 $alpha\_p1=alpha\_p\×1_2^T+{(H\_Final)}^{-1};$ 计算矩阵alpha＝[alpha_p1alpha_p]；根据式(17)计算端元矩阵A；端元矩阵A应满足非负性条件，按照式(18)计算所有像素点的丰度系数，存入数据矩阵X中。得到丰度系数矩阵完成高光谱的解混。

图2为本实施例数据点集、真实端元集、估计端元集在p-1(即3-1＝2)维上的投影显示，图3为真实端元光谱图与估计端元光谱图在224个波段上的显示。由图2可以看出估计端元与真实端元的二维表示几乎重合，图3显示估计的端元曲线与真实的端元曲线也几乎重合，这些都充分显示出了该解混算法的有效性及高精度筛选特性。

Claims (6)

1.一种基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述方法步骤如下：

步骤一、数据加载与预处理：

1)加载数据矩阵其中L是波段数，N是图像的像素点数量；

2)对数据矩阵进行降维，将数据矩阵Y降维至p-1维数据矩阵p为端元个数；

步骤二、筛选图像采样点，构造优化的约束条件，寻找满足初始条件的数据；

步骤三、将非负非线性规划问题转换为线性规划问题，对线性规划问题的目标函数结合优化后的约束条件进行求解，计算中间变量矩阵H_new、g_new；

步骤四、根据变化率检测终止条件判断是否终止迭代计算，若不满足终止条件，则返回步骤三，继续更新中间变量矩阵H_new、g_new；

步骤五、解出满足非负性要求的端元矩阵，并计算丰度系数，完成图像的解混。

2.根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述步骤二的具体步骤如下：

1)在数据矩阵中得到p个列向量：

$M_0=\[{\tilde{y}}_1^{\′},{\tilde{y}}_2^{\′},...,{\tilde{y}}_p^{\′}\],$

式中，

2)在数据集中寻找数据点集，数据集中任意一点满足不在由M₀中的列向量作为顶点构成的单形体中的条件，最后得到满足要求的数据点集

3)优化的线性不等式约束式：

$1_{p-1}^{T}H-1_{p-1}^{T}gq\≤q;$

$q=1_m^T{\tilde{Y}}_{block}^T{\left({\tilde{Y}}_{block}{\tilde{Y}}_{block}^T\right)}^{-1};$

其中，q∈R^1×(p-1)，H与g为要优化的变量矩阵；

4)求解出满足线性不等式约束式的可行解矩阵H₀与g₀，作为后续步骤的初始值。

3.根据权利要求2所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述步骤2)中，具体寻找方法如下：

首先令：

$B=\[{\tilde{y}}_1^{\′}-{\tilde{y}}_p^{\′},...,{\tilde{y}}_{p-1}^{\′}-{\tilde{y}}_p^{\′}\];$

$\tilde{s}=\left[\begin{array}{c}s\\ 1-\underset{i=1}{\overset{p-1}{\Σ}}{s}_i\end{array}\right];$

其中，s_i为s中的第i行元素；

对所有1≤i≤N求取其对应的向量若向量不满足的条件，则将加入数据集合中；假设最终有m个不满足上述条件的则最后得到的数据集合p≤m≤N，且m＜＜N。

4.根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述步骤三的具体步骤如下：

1)构造非线性目标函数：

$f(H,g)=\underset{H,g}{max}\left|\det \left(H\right)\right|;$

2)将非线性目标函数转化为线性的目标函数：

首先将H的行列式det(H)通过代数余子式的方式展开：

$\det \left(H\right)=\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right);$

式中h_ij是矩阵H的第i行第j列的元素，H_ij是将矩阵H移除第i行元素与第j列元素后得到的子阵，进而目标函数转化为：

$f(H,g)=\underset{h_i^T,g_i}{\max }\left|\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right)\right|;$

将上述目标函数分解为下面两个线性目标函数：

$p^{*}(h_i,g_i)=\underset{h_i^T,g_i}{\max }\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right),$

$q^{*}(h_i,g_i)=\underset{h_i^T,g_i}{\min }\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right);$

3)将线性目标函数与优化后的线性约束条件结合，构造如下双线性规划子问题：

$p^{*}=\underset{h_i^T,g_i}{\max }\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right),$

$s.t.h_i^T{y}_b\[i\]-g_i\≥0,\∀y_b\[i\]\∈{\tilde{Y}}_{block},1_{p-1}^{T}H-1_{p-1}^{T}gq\≤q$

$q^{*}=\underset{h_i^T,g_i}{\min }\underset{j=1}{\overset{N-1}{\Σ}}{(-1)}^{i+j}{h}_{ij}\det \left(H_{ij}\right),$

$s.t.h_i^T{y}_b\[i\]-g_i\≥0,\∀y_b\[i\]\∈{\tilde{Y}}_{block},1_{p-1}^{T}H-1_{p-1}^{T}gq\≤q$

其中，h_i为矩阵H第i行的行向量，g_i为向量g的第i个元素；

4)求解上述两个线性规划子问题，根据p^*与q^*的值，选择性地更新对应的矩阵元素h_i与g_i，直到将矩阵H与g中所有的元素都更新一遍，最终得到矩阵H与g中所有元素都更新后的新矩阵H_new与g_new。

5.根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述步骤四的具体步骤如下：

1)设置变化率检测终止条件如下：

{|det(H_new)|-|det(H₀)|}/|det(H)|＜ε，

式中，H₀代表在步骤三开始迭代之前，矩阵H的初值；

2)判断中间变量矩阵H_new是否满足终止条件，如果满足终止条件，则令H^*＝H_new，g^*＝g_new，H^*与g^*作为矩阵H与g的最终更新结果，否则，令H＝H_new，H₀＝H_new，g＝g_new，转入步骤三。

6.根据权利要求1所述的基于最小体积与优化约束条件的高光谱解混方法，其特征在于所述步骤五的具体步骤如下：

1)令α_p＝(H^*)^-1g^*，α＝[α₁,α₂,...,α_p]，则端元矩阵A为：

A＝Qα+d；

2)若矩阵A是非负矩阵，计算丰度系数矩阵为：

$X={\left[\begin{array}{cc}{(H^{*}\tilde{Y}-g^{*}\×1_N^T)}^T& 1_N-{(H^{*}\tilde{Y}-g^{*}\×1_N^T)}^T{1}_{p-1}\end{array}\right]}^T;$

式中，为行向量，且向量中每个元素值都为1，向量中元素个数为N；1_N为列向量，且向量中每个元素值都为1，向量中元素个数为N；1_p-1为列向量，且向量中每个元素值都为1，向量中元素个数为(p-1)。