一种基于边界投影最优梯度的高光谱非线性解混方法

摘要

本发明涉及一种基于边界投影最优梯度的高光谱非线性解混方法，边界投影最优梯度方法通过选择特殊的搜索点，步长由李普希兹常数决定的方式，大大加速了边界约束下的最优化的收敛速度，达到最优收敛速度

说明书

技术领域

本发明涉及高光谱图像解混领域，具体是涉及一种基于边界投影最优梯度 的高光谱非线性解混方法。

背景技术

在过去几十年中，高光谱成像已经是遥感应用的热点研究领域，如目标探 测、光谱解混以及目标匹配和分类。由于高光谱成像传感器和地表变化等原因， 混合像素广泛存在于高光谱成像中。在这种情况下，高光谱解混对于高光谱数 据后续的量化分析是非常必要的，光谱解混包括将混合像素分解成一系列的纯 净光谱特征，称为端元，以及每个像素中纯净端元所占的比例，称为丰度。光 谱解混的混合模型可以是线性的也可以是非线性的，这取决于待研究的高光谱 图像。

由于线性混合模型(LMM)相对简单且易于解释，LMM以及广泛应用在地球 科学和遥感处理领域。但是，LMM在很多情况下可能不适应，例如，当存在多散 射效应或紧密的相互作用时，非线性混合模型(NLMMs)提高了一种选择来克服 LMM的内在局限性。NLMMs以及在高光谱图像处理领域中提出，且可以分为两大 类。NLMMs的第一类是基于环境的自然特性的，它们包括基于双向反射的模型和 OlivierEches等人在《IEEEGEOSCIENCEANDREMOTESENSINGLETTERS》2014 年第11卷第4期《ABilinear–BilinearNonnegativeMatrixFactorization MethodforHyperspectralUnmixing》中提出的双线性混合模型(BMM)。NLMMs 的第二类为其它基于物理逼近的提供了更灵活的模型，它们包括Giorgio Licciardi等人在《IEEETRANSACTIONSONGEOSCIENCEANDREMOTESENSING》 2011年第49卷第11期《Pixelunmixinginhyperspectraldatabymeansof neuralnetworks》中提出的神经网络模型、YanfengGu等人在《IEEE TRANSACTIONSONGEOSCIENCEANDREMOTESENSING》2013年第51卷第7期 《Spectralunmixinginmultiple-kernelhilbertspaceforhyperspectral imagery》中提出的核模型和YoannAltmann等人在《IEEETRANSACTIONSONIMAGE PROCESSING》2014年第23卷第6期《Unsupervisedpost-nonlinearunmixing ofhyperspectralimagesusingahamiltonianmontecarloalgorithm》中 提出的后非线性模型等。BMM将地面和树冠之间的二次散射效应考虑进去，但是 BMM中的双线性互相作用通常显示出很强的相关性，这使得解混过程对噪声十分 敏感，AbderrahimHalimi等人在《IEEETRANSACTIONSONGEOSCIENCEANDREMOTE SENSING》2011年第49卷第11期《Nonlinearunmixingofhyperspectralimages usingageneralizedbilinearmodel》中提出的广义双线性模型(GBM)采用 一种有效的方式来克服BMM中的潜在假设。

基于GBM的高光谱图像非线性解混通常采用两阶段方法，主要包括两个步 骤。第一步称为端元提取，主要是从高光谱图像中提取纯净端元，这样二次交 互端元也可以通过提取的纯净端元得到。第二步称为丰度估计，主要是分别估 计纯净端元和二次交互端元对应的丰度。一些方法已经提出用来估计GBM的丰 度，如AbderrahimHalimi等人在《IEEETRANSACTIONSONGEOSCIENCEANDREMOTE SENSING》2011年第49卷第11期《Nonlinearunmixingofhyperspectralimages usingageneralizedbilinearmodel》中提出的贝叶斯方法、AbderrahimHalimi 等人在《IGARSS》2011年《Unmixinghyperspectralimagesusingthe generalizedbilinearmodel》中提出的梯度下降方法(GDA)和NaotoYokoya 等人在《IEEETRANSACTIONSONGEOSCIENCEANDREMOTESENSING》2014年第 52卷第2期《Nonlinearunmixingofhyperspectraldatausing semi-nonnegativematrixfactorization》中提出的半非负矩阵分解 (semi-NMF)等。但是，贝叶斯方法的代价是计算量大，semi-NMF方法容易收 敛到局部极值，且对初始值敏感，GDA是逐像素进行解混的方法，这阻碍我们应 用到大的高光谱图像解混中。

发明内容

为克服相应技术缺陷，本发明提出了一种基于边界投影最优梯度的高光谱 非线性解混方法方案。

本发明技术方案如下：

一种基于边界投影最优梯度的高光谱非线性解混方法，其特征在于，基于 高光谱图像解混的数学模型，并将高光谱非线性解混的广义双线性模型转化为 获取如下约束条件的最佳结果：

$\underset{A,B}{min}f(A,B)=\left|\right|Y-EA-FB|{|}_F^2,$

$s.t.A\ge 0,\underset{i=1}{\overset{M}{\Sigma }}{A}_{i,j}=1,\forall j=1,...,P,0\le B\le C,$

也就是在满足约束条件条件的情况下，获的目标函数 关于A和B的最小值，其中，A是高光谱图像纯净端元的丰 度矩阵；B是二次交互丰度矩阵；

其中，f(A,B)表示关于A和B的函数f，||.||_F表示矩阵的弗罗贝尼乌斯 Frobenius范数，s.t.表示使得条件满足，表示i从1到M求和，表示对 于任意的；Y∈R^D×P表示高光谱混合像素构成的矩阵，D和P分别表示高光谱图 像的光谱维的波段总数和空间维的像素总数，E＝[e₁,...,e_M]∈R^D×M表示高光谱图像 的纯净端元构成的矩阵，e_i(i＝1,...,M)表示第i个纯净端元，M表示高光谱图 像的纯净端元的个数，A＝[a₁,...,a_P]∈R^M×P表示丰度矩阵，a_i(i＝1,...,P)表示第i 个像素的丰度向量，表 示二次交互端元矩阵，表示阿达玛积(即点乘操作)，B∈R^M(M-1)/2×P表示二次交 互丰度矩阵，C_(i,j),k＝A_i,kA_j,k(k∈{1,...,P})，∈表示属于；

具体获取A和B的最小值包括如下步骤：

步骤1.1：初始化参数：随机初始化A^l和B^l，令l＝0，ε＝10^-6,A^l和B^l分别 表示A和B的第l次的迭代结果；

步骤1.2：更新A^l+1；通过边界投影最优梯度方法求解如下优化问题得到A^l+1：

$A^{l+1}=\underset{A}{min}f(A,B^l)=\left|\right|\left(\begin{array}{c}Y-F{B}^l\\ \delta 1_P^T\end{array}\right)-\left(\begin{array}{c}E\\ \delta 1_M^T\end{array}\right)A|{|}_F^2,$

s.t.A≥0,

其中δ用于控制丰度的和接近1的程度，表示元素全为1的P维行向量， 表示元素全为1的M维行向量；

步骤1.3：更新B^l+1；通过边界投影最优梯度方法求解如下优化问题得到B^l+1：

$B^{l+1}=\underset{B}{\min }f(A^{l+1},B)=\left|\right|Y-{\mathrm{EA}}^{l+1}-FB|{|}_F^2,$

s.t.0≤B≤C,

步骤1.4：判别收敛条件：若本流程结束，得到 高光谱图像纯净端元的丰度矩阵A和二次交互丰度矩阵B，ε是全局收敛门限， 若则l＝l+1，回转执行步骤2.2。

在上述的一种基于边界投影最优梯度的高光谱非线性解混方法，所述步骤 1.2和步骤1.3中，边界投影最优梯度方法基于以下步骤进行优化：

具体包括：

步骤1.1：初始化；V^k＝X^k,α^k＝1,k＝0,L＝||W^TW||₂，其中Z、W、H和U表示 已知的矩阵，V表示特殊选择的收缩点，β表示组合系数，L表示李普希兹常数， ||.||₂表示矩阵的2范数，k表示第k次迭代；

步骤1.2：更新X^k；

X^k＝P[Q],

$Q=V^k-\frac{1}{L}▿g\left(V^k\right)$

其中▽是梯度算子，P[Q]表示对矩阵Q每个元素Q_i,j进行如下边界投影操作：

步骤1.3：更新β^k+1；${\beta }^{k+1}=\frac{1+\sqrt{4{\left({\beta }^k\right)}^2+1}}{2}$

步骤1.4：更新V^k+1；$V^{k+1}=V^k+\frac{{\beta }^k-1}{{\beta }^{k+1}}(X^k-X^{k-1})$

步骤1.5：判别收敛条件：||T[▽g(X^k)]||_F≤max(10^-3,ε)||T[▽g(X⁰)]||_F，其中 max(a,b)表示取a和b的较大值，T[▽g(X^k)]表示对梯度矩阵▽g(X^k)中的每个元 素▽g(X^k)_i,j进行如下操作：

其中min(a,b)表示取a和b的较小值；

若||T[▽g(X^k)]||_F＞max(10^-3,ε)||T[▽g(X⁰)]||_F，则k＝k+1，回转执行步骤1.2；若 ||T[▽g(X^k)]||_F≤max(10^-3,ε)||T[▽g(X⁰)]||_F，本流程结束，求解得到X。

因此，本发明边界投影最优梯度方法通过选择特殊的搜索点，步长由李普 希兹常数决定的方式，大大加速了边界约束下的最优化的收敛速度，达到最优 收敛速度此外，边界投影最优梯度方法可以有效地应用于基于GBM的高 光谱非线性解混中，它具有收敛速度快、对初始值选择不敏感等优点。

附图说明

图1是本发明实施例的流程图。

图2是本发明实施例的MoffettField数据解混后植被的丰度图。

图3是本发明实施例的MoffettField数据解混后水的丰度图。

图4是本发明实施例的MoffettField数据解混后土壤的丰度图。

图5是本发明实施例的MoffettField数据解混后植被-水的丰度图。

图6是本发明实施例的MoffettField数据解混后植被-土壤的丰度图。

图7是本发明实施例的MoffettField数据解混后水-土壤的丰度图。

具体实施方式

下面结合附图和实施例对本发明进行进一步的说明。

参照附图1，本发明主要由2个步骤组成：建立高光谱图像解混的数学模型， 求解基于GBM的非线性解混模型。实施例选取的真实数据是MoffettField高 光谱数据集，它之前已经用于基于GBM的解混，我们选取50×50的图像来验证 实验效果，除去水汽吸收的波段剩下203个波段，它主要包含三种端元：植被、 水和土壤。为了验证提出方法的有效性，我们采用如下指标：重构误差(RE) 和平均光谱角距离(SMAD)，其定义如下：

$RE(Y,Y)=\sqrt{\frac{1}{LP}\underset{i=1}{\overset{P}{\Sigma }}\left|\right|y_i-y_i|{|}^2},$

$SMAD(Y,Y)=\underset{i=1}{\overset{P}{\Sigma }}{\cos }^{-1}\left(\frac{y_i^T{y}_i}{\sqrt{y_i^T{y}_i}\sqrt{{y_i}^T{y}_i}}\right),$

其中y_i和y_i分别表示重构像素和参考像素，cos^-1表示反余弦。

具体实施时，本发明技术方案可采用计算机软件技术实现自动运行流程。 实施例执行步骤如下：

(1)建立高光谱图像解混的数学模型，将高光谱非线性解混的广义双线性 模型(GBM)转化为如下优化问题：

$\underset{A,B}{min}f(A,B)=\left|\right|Y-EA-FB|{|}_F^2,$

$s.t.A\ge 0,\underset{i=1}{\overset{M}{\Sigma }}{A}_{i,j}=1,\forall j=1,...,P,0\le B\le C,$

上述式子表示求目标函数关于A和B的最小值，并满足条 件$A\ge 0,\underset{i=1}{\overset{M}{\Sigma }}{A}_{i,j}=1,\forall j=1,...,P;$

其中，min是最小化算子，f(A,B)表示关于A和B的函数f，||.||_F表示矩阵的 弗罗贝尼乌斯Frobenius范数，s.t.表示使得条件满足，表示i从1到M求 和，表示对于任意的。Y∈R^D×P表示高光谱混合像素构成的矩阵，D和P分别表 示高光谱图像的光谱维的波段总数和空间维的像素总数，E＝[e₁,...,e_M]∈R^D×M表示 高光谱图像的纯净端元构成的矩阵，e_i(i＝1,...,M)表示第i个纯净端元，M表 示高光谱图像的纯净端元的个数，A＝[a₁,...,a_P]∈R^M×P表示丰度矩阵，a_i(i＝1,...,P) 表示第i个像素的丰度向量，表示二次交互端元矩阵，表示阿达玛积(即点乘操作)， B∈R^M(M-1)/2×P表示二次交互丰度矩阵，C_(i,j),k＝A_i,kA_j,k(k∈{1,...,P})，∈表示属于；

(2)，求解基于GBM的非线性解混模型，获得高光谱图像纯净端元的丰度 矩阵A和二次交互丰度矩阵B，其具体求解流程如下：

(2.1)：初始化参数：

随机初始化A^l和B^l，令l＝0，ε＝10^-6,A^l和B^l分别表示A和B的第l次的迭 代结果；

(2.2)：更新A^l+1；

通过边界投影最优梯度方法求解如下优化问题得到A^l+1：

$A^{l+1}=\underset{A}{min}f(A,B^l)=\left|\right|\left(\begin{array}{c}Y-F{B}^l\\ \delta 1_P^T\end{array}\right)-\left(\begin{array}{c}E\\ \delta 1_M^T\end{array}\right)A|{|}_F^2,$

s.t.A≥0,

其中δ用于控制丰度的和接近1的程度，表示元素全为1的P维行向量， 表示元素全为1的M维行向量；

(2.3)：更新B^l+1；

通过边界投影最优梯度方法求解如下优化问题得到B^l+1：

$B^{l+1}=\underset{B}{\min }f(A^{l+1},B)=\left|\right|Y-{\mathrm{EA}}^{l+1}-FB|{|}_F^2,$

s.t.0≤B≤C,

(2.4)：判别收敛条件：

若本流程结束，得到高光谱图像纯净端元的丰 度矩阵A和二次交互丰度矩阵B，ε是全局收敛门限，若 则l＝l+1，回转执行步骤2.2。

步骤(2)中所述的边界投影最优梯度方法，它是用于求解如下边界约束的 最优化问题：

$\underset{X}{\min }g\left(X\right)=\left|\right|Z-WX|{|}_F^2,$

s.t.H≤X≤U,

其求解过程包括以下子步骤：

步骤1)：初始化；

V^k＝X^k,α^k＝1,k＝0,L＝||W^TW||₂，其中Z、W、H和U表示已知的矩阵，V表示 特殊选择的收缩点，β表示组合系数，L表示李普希兹常数，||.||₂表示矩阵的 2范数，k表示第k次迭代；

步骤2)：更新X^k；

X^k＝P[Q],

$Q=V^k-\frac{1}{L}▿g\left(V^k\right)$

其中▽是梯度算子，P[Q]表示对矩阵Q每个元素Q_i,j进行如下边界投影操作：

步骤3)：更新β^k+1；

${\beta }^{k+1}=\frac{1+\sqrt{4{\left({\beta }^k\right)}^2+1}}{2}$

步骤4)：更新V^k+1；

$V^{k+1}=V^k+\frac{{\beta }^k-1}{{\beta }^{k+1}}(X^k-X^{k-1})$

步骤5)：判别收敛条件：

||T[▽g(X^k)]||_F≤max(10^-3,ε)||T[▽g(X⁰)]||_F，其中max(a,b)表示取a和b的较大 值，T[▽g(X^k)]表示对梯度矩阵▽g(X^k)中的每个元素▽g(X^k)_i,j进行如下操作：

其中min(a,b)表示取a和b的较小值；

若||T[▽g(X^k)]||_F＞max(10^-3,ε)||T[▽g(X⁰)]||_F，则k＝k+1，回转执行步骤1.2；若 ||T[▽g(X^k)]||_F≤max(10^-3,ε)||T[▽g(X⁰)]||_F，本流程结束，求解得到X。边界投影最优 梯度方法收敛速度为

实施例中，D＝203，P＝2500，M＝3，ε＝10^-6。对MoffettField高光谱数据集 进行非线性解混后，得到的植被、水、土壤、植被-水、植被-土壤和水-土壤的 丰度图分别如附图2-7所示。通过对比算法GDA和semi-NMF，RE和SMAD如表1 所示。

表1比较REs(10^-2)和SMADs(10^-2)

应当理解的是，本说明书未详细阐述的部分均属于现有技术。

应当理解的是，上述针对实施例的描述较为详细，并不能因此而认为是对 本发明专利保护范围的限制，本领域的普通技术人员在本发明的启示下，在不 脱离本发明权利要求所保护的范围情况下，还可以做出替换或变形，均落入本 发明的保护范围之内，本发明的请求保护范围应以所附权利要求为准。

总体来说，提出了一种基于边界投影最优梯度的高光谱非线性解混方法， 边界投影最优梯度方法收敛速度为可以有效地应用于基于GBM的高光谱 非线性解混中，它具有收敛速度快、对初始值选择不敏感等优点。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属 技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采 用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定 义的范围。