CN103413292A - 基于约束最小二乘的高光谱图像非线性丰度估计方法 - Google Patents

Abstract

本发明属于遥感图像处理技术领域，具体为一种基于约束最小二乘的高光谱图像非线性丰度估计方法。该发明通过在目标函数中引入丰度的非负及和为一约束以及非线性参数的有界约束，将高光谱图像非线性解混问题转化为求解丰度矢量和非线性参数的约束非线性最小二乘问题。进而，该发明采用一种交替迭代优化算法求解该问题。该发明从高光谱观测像素的混合模型出发，结合模型中丰度和非线性物理意义，有效地克服了线性解混的不足，同时具有良好的抗噪声性能，可以作为一种解决高光谱遥感图像非线性解混的有效手段。在基于高光谱遥感图像的高精度的解混以及地面目标的检测和识别方面具有重要的应用价值。

Description

基于约束最小二乘的高光谱图像非线性丰度估计方法

技术领域

本发明属于遥感图像处理技术领域，具体涉及一种高光谱遥感图像非线性丰度估计方法。

背景技术

遥感技术是本世纪六十年代发展起来的新兴综合技术，与空间、电子光学、计算机、地理学等科学技术紧密相关，是研究地球资源环境的最有力的技术手段之一。高光谱遥感是将成像技术与光谱技术相结合的多维信息获取技术。高光谱成像仪在电磁波谱的数十至数百个非常窄且连续的光谱段上同时探测目标的二维几何空间与一维光谱信息。高光谱图像中，每一个观测像素都可以提取出一条完整连续的光谱曲线，为地物信息的提取和分析提供了极其丰富的信息，有助于更加精细的地物分类和目标识别。由于受到空间分辨率、多次散射效应以及地物分布异质性的限制，高光谱传感器获得的像元光谱通常是由多种地物光谱混合而成的。这种包含不止一种地物覆盖类型的像元被称为混合像元。混合像元的存在使得传统的像元级高光谱图像的分类、压缩和目标检测等应用受到限制。因此，以提取基本组成成分(端元)的特征光谱和估计它们在每个观测像素中所占的比例(丰度)为主要内容的高光谱解混，成为高光谱图像分析中的关键技术之一。它可以大大提高高光谱遥感对亚像元级目标的探测能力，被广泛应用于地质科学、水文科学、精准农业以及军事领域等[1]，[2]。

通常，混合像元所记录的光谱信号被认为是端元光谱的线性组合，其中的组合系数被认为是它们相应的丰度。满足这种假设的模型被称为线性混合模型(Linear Mixing Model，LMM)。LMM已经被广泛应用于高光谱解混的研究中。但是，对于一些包含砂石、矿物、植被以及水域等地物分布的高光谱图像来说，由于非线性混合现象的存在，LMM则会得到不精确地结果。为此，研究者根据高光谱数据自身特性以及相关的辐射传输理论(Radiative Transfer Theory，RTT)[3]，提出了一些非线性混合模型及其对应的解混算法。这些模型可以有效地克服LMM的局限性，提高解混的精度，在实际应用中有重要的现实意义[3]-[8]。常见的有Fan模型(FanModel，FM)[6]、广义双线性模型(Generalized Bilinear Model，GBM)[7]，多项式后非线性混合模型(Polynomial Postnonlinear Mixing Model，PPNMM)[8]等。通过在线性混合模型的基础上增加非线性项来描述非线性混合，这些模型可以有效的表示高光谱图像中的非线性混合现象。

下面介绍与本发明相关的一些概念：

线性和非线性光谱混合模型

对于高光谱数据集X=(x₁，x₂，·，x_N)∈·^L×N来说，每一个观测像素x∈·^L×1都可以用端元矩阵M=(m₁，m₂，·，m_p)∈·^L×p及其相应的丰度矢量s=(s₁，s₂，·，s_p)^T∈·^p×1表示如下

x=f(M，s，γ)+n          (10)

其中，L、N和p分别表示波段数、像素数和端元数。f(·)是一个含有固定参数γ的未知线性或非线性函数。n∈·^L×N则表示了噪声和可能的模型误差。同时，丰度矢量应满足丰度和为一约束(Abundance Sum-to-one Constraint，ASC)和丰度非负约束(Abundance NonnegativeConstraint，ANC)：

$\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{s}_i=1,s_i\≥0,\∀i\∈\{\mathrm{1,2},\·,p\}.---\left(11\right)$

对于LMM来说，

f(M，s)=Ms              (12)

即f(·)为线性函数。而根据GBM，非线性函数f(·)可以表述为

$f(M,s,\mathrm{\γ})=\mathrm{Ms}+\underset{i=1}{\overset{p-1}{\mathrm{\Σ}}}\underset{j=i+1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\γ}}_{i,j}{s}_i{s}_j{m}_i\·m_j---\left(13\right)$

其中，·表示Hadamard乘操作，γ=(γ_1,2,γ_1,3，·，γ_p-1，p)^T∈·^{(p(p-1)/2)×1}被称为非线性参数。γ满足约束条件

$0\≤{\mathrm{\γ}}_{i,j}\≤1,\∀i\∈\{\mathrm{1,2},\·,p-1\},\∀j\∈\{i+1,\·,p\}---\left(14\right)$

很显然，LMM是GBM的特例。当 ${\mathrm{\γ}}_{i,j}=1,\∀i\∈\{\mathrm{1,2},\·,p-1\},\∀j\∈\{i+1,\·,p\}$ 时，GBM则变成FM。在PPNMM下，函数f(·)可以表述为

f(M，s，b)=Ms+b(Ms)·(Ms)           (15)

其中，非线性参数b∈·。当b=0时，PPNMM转变为LMM。

发明内容

本发明的目的在于提出一种解混精度高且运算复杂度低的高光谱遥感图像非线性丰度估计方法。

本发明首先提出一种可以有效概括现有的线性和非线性混合模型的高光谱图非线性解混框架。该框架将待估计的参数分为引入线性混合的丰度矢量和引入非线性混合的表征非线性形式或强弱的非线性参数，并根据它们所满足的约束条件将非线性解混问题转化为一种约束的非线性最小二乘问题。在此基础上，本发明提出一种基于约束非线性最小二乘的交替迭代优化算法实现非线性解混。所提出的算法根据丰度矢量和非线性参数所遵循的约束条件设计不同的目标函数，进而采用不同的优化算法确定最优的估计结果。与其他优秀的同类方法相比，本发明具有更好的解混精度和抗噪声性能以及较低的运算复杂度，而且对于线性和非线性模型都具有较好的适应性。

本发明提出一种基于约束最小二乘的高光谱图像非线性丰度估计方法，具体内容介绍如下。

一、基于约束最小二乘的非线性解混框架的建立

不失一般性，式(9)可以表述为

$x=\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{\stackrel{\·}{m}}_i{s}_i+n=\stackrel{\·}{M}\left(\mathrm{\θ}\right)s+n,\mathrm{\θ}={(s^T,{\mathrm{\γ}}^T)}^T---\left(16\right)$

其中s=(s₁，s₂，·，s_p)^T∈·^p×1，γ=(γ₁，γ₂，·，γ_q)^T∈·^q×1，q表示非线性参数的个数。对于式(16)来说，丰度矢量s引入了线性，而参数θ引入了非线性。它们所要满足的约束条件为

$\underset{i=1}{\overset{p}{\mathrm{\Σ}}}{s}_i=1,s_i\≥0,\∀i\∈\{\mathrm{1,2},\·,p\}---\left(17\right)$

${\mathrm{\γ}}_i\∈{\·}_{{\mathrm{\γ}}_i},\∀i\∈\{\mathrm{1,2},\·,q\}$

·γ_i即为参数γ_i所属的实数区间。端元光谱矩阵M和引入非线性后的特征光谱矩阵

之间存在函数关系

g(·)：·^L×p→·^L×p

$\stackrel{\·}{M}=({\stackrel{\·}{m}}_1,{\stackrel{\·}{m}}_2,\·,{\stackrel{\·}{m}}_p)\·M=(m_1,m_2,\·,m_p)---\left(18\right)$

对于GBM来说

${\stackrel{\·}{m}}_i=m_i+\underset{j=i+1}{\overset{p}{\mathrm{\Σ}}}{\mathrm{\γ}}_{i,j}{s}_j{m}_i\·m_j,\∀i\∈\{\mathrm{1,2},\·,p\}---\left(19\right)$

而对于PPNMM来说

${\stackrel{\·}{m}}_i=m_i+b\underset{j=1}{\overset{p}{\mathrm{\Σ}}}{s}_j{m}_i\·m_j,\∀i\∈\{\mathrm{1,2},\·,p\}---\left(20\right)$

式(16)揭示了一种新的解混框架。从物理意义上讲，新框架将到达传感器的光谱仍然看作特征光谱的线性组合，而且组合系数仍然是各个端元所对应的丰度。所不同的是，此时的特征光谱不再是原始意义上的端元光谱，而是在原始端元光谱基础上通过某种非线性映射得到的衍生光谱。这种非线性映射的强弱是由参数θ决定的。它包括两部分：决定混合过程中发生的非线性映射的形式及强弱等因素的固有参数γ和地物的分布参数s。

根据式(16)，高光谱解混可以通过估计θ=(s^T，γ^T)^T实现，需要解决一个约束的非线性最小二乘问题

$\mathrm{\θ}=\underset{\mathrm{\θ}}{\arg \min }J\left(\mathrm{\θ}\right)=\underset{\mathrm{\θ}}{\arg \min }\left\{\frac{1}{2}{\left|\right|x-\stackrel{\·}{M}\left(\mathrm{\θ}\right)s\left|\right|}^2\right\},s.t.\left(17\right)---\left(21\right)$

其中，x表示输入观测像素，θ表示非线性系统的输出。

本文中采用交替迭代优化的算法求解最优的θ。

二、约束条件的添加

根据丰度的物理意义，任意观测像素中，各个端元的丰度之和必须为1。本发明中ASC可以通过最小化目标函数实现

$f_{\mathrm{ASC}}\left(s\right)=\frac{1}{2}{\left|\right|(1^{T}s-1)\left|\right|}^2---\left(22\right)$

其中1∈·^p×1是元素全为一的列矢量。当丰度矢量s满足ASC时，f_ASC(s)=0，否则f_ASC(s)是一个较大的正数。

由式(17)的定义可知，任意观测像素对应的丰度和函数f(·)所固有的非线性参数(如果该参数存在的话)都可以认为是满足一定的有界约束。例如，丰度矢量s中任意元素s_i应满足s_i∈[0，1]，GBM中任一非线性参数均在区间[0，1]上。这样，ANC也可以被满足。文中设计了一种新的目标函数以强制这些参数在特定区间上。

以非线性参数γ为例，假设γ=(γ₁，γ₂，·，γ_q)^T中元素均在区间[a，b]上，则可以通过最小化目标函数实现这一有界约束：

$f_B\left(\mathrm{\γ}\right)=\underset{i=1}{\overset{q}{\mathrm{\Σ}}}\mathrm{\φ}\left({\mathrm{\γ}}_i\right),\mathrm{\φ}\left({\mathrm{\γ}}_i\right)=\left\{\begin{array}{cc}0,& {\mathrm{\γ}}_i\∈[a,b]\\ +\∞,& {\mathrm{\γ}}_i\∉[a,b]\end{array}\right.---\left(23\right)$

很显然上述φ(γ_i)使用起来较为困难。基于此，我们采用与文献[9]中类似的幂函数来代替式(23)中的φ(·)

$\mathrm{\&phi;}\left(x\right)=\left\{\begin{array}{cc}{(x-a)}^{2n}/2n,& x<a\\ 0,& a\&le;x\&le;b\\ {(x-b)}^{2n}/2n,& x>b\end{array}\right.---\left(24\right)$

其中n为正整数。当γ中所有元素均在区间[a，b]上时，f_B(γ)=0；否则，f_B(γ)＞0。该函数的导函数为

$\frac{\&PartialD;f_B\left(\mathrm{\&gamma;}\right)}{\&PartialD;\mathrm{\&gamma;}}={(\frac{\&PartialD;f_B\left(\mathrm{\&gamma;}\right)}{\&PartialD;{\mathrm{\&gamma;}}_1},\frac{\&PartialD;f_B\left(\mathrm{\&gamma;}\right)}{\&PartialD;{\mathrm{\&gamma;}}_2},\&CenterDot;,\frac{\&PartialD;f_B\left(\mathrm{\&gamma;}\right)}{\&PartialD;{\mathrm{\&gamma;}}_q})}^T---\left(25\right)$

其中

$\frac{\&PartialD;f_B\left(\mathrm{\&gamma;}\right)}{\&PartialD;{\mathrm{\&gamma;}}_i}=\frac{\&PartialD;\mathrm{\&phi;}\left({\mathrm{\&gamma;}}_i\right)}{\&PartialD;{\mathrm{\&gamma;}}_i}=\left\{\begin{array}{cc}{({\mathrm{\&gamma;}}_i-a)}^{2n-1},& {\mathrm{\&gamma;}}_i<a\\ 0,& a\&le;{\mathrm{\&gamma;}}_i\&le;b,\\ {({\mathrm{\&gamma;}}_i-b)}^{2n-1},& {\mathrm{\&gamma;}}_i>b\end{array}\right.i=\mathrm{1,2},\&CenterDot;,q.---\left(26\right)$

文中取n=1。新的目标函数连续可导，可以保证迅速收敛到预定空间。

三、约束的非线性最小二乘算法

综合考虑目标函数(20)和上述两个约束条件的添加方式，解混问题转化为两个约束优化问题如下

${\mathrm{\&gamma;}}^{(t+1)}=\underset{\mathrm{\&gamma;}}{\arg \min }{J}_1\left(\mathrm{\&gamma;}\right)=\underset{\mathrm{\&gamma;}}{\arg \min }\{\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}(s^{\left(t\right)},\mathrm{\&gamma;})s^{\left(t\right)}\left|\right|}^2+{\mathrm{\&lambda;}}_1{f}_B\left(\mathrm{\&gamma;}\right)\}---\left(27\right)$

$s^{(t+1)}=\underset{s}{\arg \min }{J}_2\left(s\right)=\underset{s}{\arg \min }\{\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}\left({\mathrm{\&theta;}}^{(t+1)}\right)s\left|\right|}^2+{\mathrm{\&lambda;}}_1{f}_{\mathrm{ASC}}\left(s\right)+{\mathrm{\&lambda;}}_2{f}_B\left(s\right)\}---\left(28\right)$

其中，θ^(t+1)=((s^(t))^T，(γ^(t+1))^T)^T，λ₁和λ₂是表示所添加约束的强弱的正数。

问题(26)是一个约束的非线性优化问题。我们采用梯度下降法求解该问题。其迭代步骤为

${\mathrm{\&gamma;}}^{(t+1)}={\mathrm{\&gamma;}}^{\left(t\right)}-{\mathrm{\&lambda;}}^{\left(t\right)}\frac{\&PartialD;J_1\left(\mathrm{\&gamma;}\right)}{\&PartialD;\mathrm{\&gamma;}}---\left(29\right)$

其中λ^(t)可以采用固定步长也可以通过线性搜索的方法求得。文中我们采用固定步长λ^(t)=0.5。

问题(28)则是一个全约束的线性最小二乘问题，同样可以采用梯度法求解。文中我们采用FCLS算法[10]求解。

方法：基于约束非线性最小二乘的丰度估计算法(Constrained Nonlinear LeastSquares-based Abundance Estimation，CNLS-AE)

根据上述内容，本发明采用的算法中采用的具体步骤归纳如下：

1、建立基于约束最小二乘的非线性解混框架；

2、采用交替迭代优化算法，输入：高光谱数据集X∈·^L×N，端元光谱矩阵M∈·^L×p输出：估计的丰度矩阵S=(s₁，s₂，·，s_N)∈·^p×N，非线性参数矩阵Γ=(γ₁，γ₂，·，γ_N)∈·^q×N。交替迭代优化算法具体步骤如下，其中γ^(t)和s^(t)分别表示第t迭代后的结果。

步骤1.初始化S和Γ采用FCLS算法或者随机初始化得到初始的丰度矩阵

随机选取满足要求的初始非线性参数矩阵

${\mathrm{\&Gamma;}}^0=({\mathrm{\&gamma;}}_1^0,{\mathrm{\&gamma;}}_2^0,\&CenterDot;,{\mathrm{\&gamma;}}_N^0);$

其中，

为第n个丰度矢量的初始值，

为第n个非线性参数的初始值，n=1，2，·，N。

步骤2.对于每一个观测像素x_n(n=1，2，·，N)执行以下循环：

2a)设置迭代次数t=0，x=x_n，

2b)依次求解优化问题(26)和(27)得到γ^(t+1)和s^(t+1)；

2c)判断函数J(θ)是否已收敛：若收敛输出x_n的丰度矢量s_n=s^(t+1)和非线性矢量γ_n=γ^(t+1)，并处理下一个观测像素，否则返回2b)继续执行。

步骤3.输出结果S=(s₁，s₂，·，s_N)和Γ=(γ₁，γ₂，·，γ_N)。

我们采用每次迭代与前一次迭代之间J(θ)的相对变化

$\mathrm{\&xi;}=\frac{{J\left(\mathrm{\&theta;}\right)}^{t+1}-{J\left(\mathrm{\&theta;}\right)}^t}{{J\left(\mathrm{\&theta;}\right)}^t},$ $J\left(\mathrm{\&theta;}\right)=\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}\left(\mathrm{\&theta;}\right)s\left|\right|}^2+{\mathrm{\&lambda;}}_1{f}_B\left(\mathrm{\&gamma;}\right)+{\mathrm{\&lambda;}}_1{f}_{\mathrm{ASC}}\left(s\right)+{\mathrm{\&lambda;}}_2{f}_B\left(s\right)$

                                                                         (30)作为判断是否收敛的标准。当该值小于某个阈值η(我们取η=10^-6)时，可以认为目标函数已经收敛。

本发明的有益效果在于：其从高光谱观测像素的混合模型出发，结合模型中丰度和非线性物理意义，有效地克服了线性解混的不足，同时具有良好的抗噪声性能，可以作为一种解决高光谱遥感图像非线性解混的有效手段。在基于高光谱遥感图像的高精度的解混以及地面目标的检测和识别方面具有重要的应用价值。

仿真和实际高光谱数据实验表明，与同类算法相比，本发明方法中采取的算法具有更好的丰度估计精度和抗噪声性能以及较低的运算复杂度，而且对于线性和非线性模型都具有较好的适应性，这对于处理复杂的高光谱数据来说具有重要的实际意义。

附图说明

图1USGS光谱库中得到的五种光谱：(a)Alunite GDS83Na63，(b)Nontronite GDS41，(c)Pyrope WS474，(d)Buddingtonite NHB2301，(e)Montmorillonite CM20。

图2LMM下算法性能的比较.(a)RMSE，(b)RE。

图3FM下算法性能的比较.(a)RMSE，(b)RE。

图4GBM下算法性能的比较.(a)RMSE，(b)RE。

图5PPNMM下算法性能的比较.(a)RMSE，(b)RE。

图6不同模型下四种算法的运行时间比较.(a)LMM，(b)FM，(c)GBM，(d)PPNMM。

图7不同端元数目下算法性能的比较.(a)RMSE，(b)RE。

图8Moffet Field数据集图像块：(a)波段50的灰度图，(b)由N-FINDR算法得到的三种端元光谱。

图9四种算法得到的RE图：(a)FCLS，(b)GBM-GDA，(c)PPNMM-GDA，(d)CNLS-AE。

具体实施方式

下面，分别用仿真数据和实际遥感图像数据为例说明本发明的具体的实施方式：

本发明中采用的基于约束最小二乘的非线性丰度估计方法(算法)用CNLS-AE表示。

1、仿真数据实验

在本节中，我们采用人工产生的仿真数据测试所提出的算法的性能。我们将本文提出的算法与文献[7]和[8]中提到的分别基于GBM和PPNMM的两种梯度法，即GBM-GDA和PPNMM-GDA算法作比较。此外，我们还要与性能较优的基于LMM的FCLS算法[10]进行比较。

我们使用均方根误差(Root Mean Square Error，RMSE)和重构误差(Reconstruction Error，RE)来衡量丰度估计算法的优劣。RMSE用于衡量丰度矩阵的估计结果与真实值的近似程度。假设算法得到的端元丰度为

真实丰度为S={s₁，s₂，·，s_N}∈·^p×N，则RMSE被定义为[11]

$\mathrm{RMSE}=\frac{1}{p}\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}\left[\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\stackrel{\&CenterDot;}{S}}_{\mathrm{ij}}-S_{\mathrm{ij}})}^2}\right],---\left(31\right)$

其中和S_ij分别为矩阵

和S的第i行第j列的元素。RE则衡量了算法对模型的适用程度，它被定义为[12]

$\mathrm{RE}=\frac{1}{L}\underset{i=1}{\overset{L}{\mathrm{\&Sigma;}}}\left[\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{\mathrm{\&Sigma;}}}{({\hat{X}}_{\mathrm{ij}}-X_{\mathrm{ij}})}^2}\right],---\left(32\right)$

其中X=(x₁，x₂，·，x_N)∈·^L×N和

分别是原始数据矩阵和从估计结果重构得到的数据矩阵，即 ${\stackrel{\&CenterDot;}{x}}_n=f(M,{\stackrel{\&CenterDot;}{s}}_n,\mathrm{\&xi;}),$ $\&ForAll;n\&Element;\{\mathrm{1,2},\&CenterDot;,N\}.$ 和X_ij分别为矩阵和x的第i行第j列的元素。

我们还测算各个算法的运行时间以定量衡量它们的复杂度。所有实验的软硬件环境为Intel(R)Xeon(R)E5504CPU2.00GHz，24GB内存和Windows7及Matlab7.0.

仿真试验中采用的模拟数据共有四组，分别依据LMM、FM、GBM和PPNMM产生。其中GBM的非线性参数γ_s为区间[0，1]上的随机实数，PPNMM中的非线性参数b∈[-0.3，0.3]。端元光谱是从USGS光谱库得到20种不同的矿物光谱数据，其中五种矿物光谱曲线如图1所示。它们共有224个波段，覆盖0.37～2.48μm的波长范围。而丰度矩阵则是由Dirichlet分布得到，它满足ANC和ASC。每组数据包含2500个像素。构造的图像上叠加由信噪比(Signal-to-Noise Ratio，SNR)衡量的不同强度的高斯白噪声，就得到了实验用的仿真数据。每组解混试验都在同等条件下执行40次，然后使用平均结果作为最终的输出结果，以避免单次实验造成的误差。

实验1抗噪声性能实验在这个实验中，我们比较四种算法的抗噪声性能。我们改变仿真数据的SNR，将其从∞(无噪声)、50dB降至10dB，每隔5dB为一档。所采用的端元光谱为如图1所示。图2～5分别给出了不同的混合模型下，四种算法的RMSE和RE关于SNR的变化曲线。我们可以看到随着噪声的增加，四种算法的性能都逐渐下降。在LMM下，FCLS算法总是得到最小的RMSE，而我们提出的CNLS-AE算法与之很接近。由此可见，CNLS-AE算法对线性模型的数据也有较好的性能。在非线性混合模型，即FM、GBM和PPNMM下，FCLS性能则较差。三种非线性丰度估计算法相比，CNLS-AE算法总是具有较好的性能，即具有相对小的RMSE和RE。此外，我们还可以看出，在FM和GBM下，PPNMM-GDA也可以得到优于GBM-GDA的结果。这是由于当实际混合模型与所假设的模型存在偏差时，PPNMM与FM和GBM相比具有较强的鲁棒性，这也吻合了文献[10]的结论。

为了测试算法的复杂度，在这部分实验中，我们还比较了四种算法在不同模型下的运行时间。图6(a)～(d)给出了四种模型下四种算法的运行时间变化情况。很显然，FCLS算法的运行时间总是最小的。三种非线性丰度估计算法相比，提出的CNLS-AE算法的运行时间则是最小的。这说明提出的算法具有较小的计算复杂度。

实验2端元数目与算法性能实验在这个实验，我们调查四种算法的性能随着端元数变化的情况。我们采用的端元数目从3增加到15。每次试验时，我们使用从20种光谱数据中随机选取所需数目的光谱产生模拟数据。模拟数据是依据GBM产生的，而且所含噪声的SNR为30dB。

图7给出了四种算法的RMSE和RE随端元数目增加的变化情况。在所选的端元数目范围内，提出的算法的RMSE总是小于GBM-GDA和FCLS算法，具有较好的性能。

2、实际数据实验

在本节中，我们使用实际的高光谱遥感图像数据集对所提出算法的性能进行测试。我们使用的数据集是由机载可见光及红外成像光谱仪(Airborne Visible/Infrared Imaging Spectrometer，AVIRIS)拍摄于1997年6月19日的Moffett Field数据集。该数据集包含512×614个像素，224个波段，波长范围是0.37-2.48μm，光谱分辨率为10nm。波段1～10，108～113，153～168和222-224因为信噪比太低或为水吸收波段而被移除，留下的189个波段被用于算法验证。我们选取大小为50×50的较为典型的子图像块用于实验。图8(a)显示了该图像块的波段50的灰度图。实地勘测可知，该地区包含三种地物，即土壤，水和植被。而且该地区存在明显的非线性现象。图8(b)给出了采用N-FINDR[13]算法得到的这三种矿物的光谱。

图9给出了四种算法得到RE图。可以看出，我们提出的算法具有普遍较小的RE。为了定量的比较四种算法的性能，我们在表1中给出了它们对应的RE和单像素的平均运算时间。由表1可知，针对该数据集，我们提出的算法可以得到最小的RE，而且与其他两种非线性丰度估计算法相比，运行时间也是最少的。

表1四种算法的平均RE和单像素的平均运行时间比较

综上可知，对于模拟和实际高光谱数据来说，我们提出的算法相对于其他类似算法而言，都具有较好的估计精度和较快的收敛速度，可以快速而精确地实现高光谱图像的解混。

参考文献：

[1]C.I.Chang.Hyperspectral imaging：techniques for spectral detection and classification.NewYork：Plenum，2003.

[2]N.Keshava and J.F.Mustard，“Spectral unmixing，”IEEE Signal Processing Magazine，vol.19，no.1，pp.44-57，Jan.2002.

[3]S.Chandrasekhar，Radiative Transfer：New York：Dover，1960.

[4]B.Hapke，“Bidirection reflectance spectroscopy.I.theory，”J.Geophys.Res.，vol.86，pp.3039-3054，1981.

[5]P.Kulbelka and F.Munk，“Reflection characteristics of paints，”Zeitschriftfur Technische Physik，vol.12，pp.593-601，1931.

[6]W.Fan，B.Hu，J.Miller，and M.Li，“Comparative study between a new nonlinear model andcommon linear model for analysing laboratory simulated-forest hyperspectral data，”Int.J.Remote Sens.，vol.30，no.11，pp.2951-2962，Jun.2009.

[7]A.Halimi，Y.Altmann，N.Dobigeon，and J.-Y.Tourneret，“Nonlinear unmixing of hyperspectralimages using a generalized bilinear model，”IEEE Trans.Geosci.Remote Sens.，vol.49，no.11，pp.4153-4162，Nov.2011.

[8]Y.Altmann，A.Halimi，N.Dobigeon，and J.-Y.Tourneret，“Supervised nonlinear spectralunmixing using a postnonlinear mixing model for hyperspectral imagery，”IEEE Trans.ImageProcessing，vol.21，no.6，pp.3017-3025，2012.

[9]W.Xia，X.S.Liu，B.Wang and L.M.Zhang，“Independent Component Analysis for BlindUnmixing of Hyperspectral Imagery With Additional Constraints，”IEEE Trans.Geosci.RemoteSens.，Vol.49，No.6，pp.21652179，Jun.2011

[10]D.C.Heinz and C.I.Chang，“Fully constrained least square linear spectral unmixing analysismethod for material quantification in hyperspectral imagery，”IEEE Trans.Geosci.Remote Sens.，vol.39，no.3，pp.529-545，Mar.2001.

[11]A.Plaza，P.Martinez，R.Perez，et al.“A quantitative and comparative analysis of endmemberextraction algorithms from hyperspectral data”.IEEE Trans.Geosci.RemoteSens.，vol.42，no.3，pp.650-663，2004

[12]J.Plaza，E.M.T.Hendrix，I.García，et al.“On endmember identification in hyperspectralimages without pure pixels：A comparison of algorithms”.Journal of Mathematical Imaging andVision，vol.42，pp.163-175，2012

[13]M.EWinter，“N-FINDR：An algorithm for fast autonomous spectral endmember determinationin hyperspectral data，”inProc.of the SPIE conference on Imaging Spectrometry V，1999，vol.3753，pp.266-277。

Claims (6)

1.一种基于约束最小二乘的高光谱图像非线性丰度估计方法，其特征在于，在现有混合模型的基础上，将高光谱图像解混问题转换为一种约束的非线性最小二乘问题，并采用交替迭代优化算法求解该问题；具体步骤如下：

(1)建立基于约束最小二乘的非线性解混框架，将高光谱数据集X=(x₁，x₂，·，x_N)∈·^L×N中的每一个观测像素x∈·^L×1用端元光谱矩阵M=(m₁，m₂，·，m_p)∈·^L×p及其相应的丰度矢量s=(s₁，s₂，·，s_p)^T∈·^p×1表示：

$x=\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{\stackrel{\&CenterDot;}{m}}_i{s}_i+n=\stackrel{\&CenterDot;}{M}\left(\mathrm{\&theta;}\right)s+n,\mathrm{\&theta;}={(s^T,{\mathrm{\&gamma;}}^T)}^T---\left(1\right)$

其中，n表示噪声和可能的模型误差；L、N、p和q分别表示波段数、像素数、端元数和非线性参数的个数；θ为参数，决定了非线性映射的强弱；为由θ决定的引入非线性后的特征光谱矩阵；γ为非线性参数，表示混合过程中发生的非线性映射的形式及强弱等因素，γ=(γ₁，γ₂，·，γ_q)^T∈·^q×1；

同时，s、γ满足以下约束条件：

$\underset{i=1}{\overset{p}{\mathrm{\&Sigma;}}}{s}_i=1,s_i\&GreaterEqual;0,\&ForAll;i\&Element;\{\mathrm{1,2},\&CenterDot;,p\}---\left(2\right)$

${\mathrm{\&gamma;}}_i\&Element;{\&CenterDot;}_{{\mathrm{\&gamma;}}_i},\&ForAll;i\&Element;\{\mathrm{1,2},\&CenterDot;,q\}$

·γ_i为参数γ_i所属的实数区间；

端元光谱矩阵M和引入非线性后的特征光谱矩阵

之间存在函数关系g(·)，它被定义为：

g(·)：·^L×p→·^L×p

$\stackrel{\&CenterDot;}{M}=({\stackrel{\&CenterDot;}{m}}_1,{\stackrel{\&CenterDot;}{m}}_2,\&CenterDot;{,\stackrel{\&CenterDot;}{m}}_p)\&CenterDot;M=(m_1,m_2,\&CenterDot;,m_p)---\left(3\right)$

(2)采用交替迭代优化算法，输入高光谱数据集X∈·^L×N和端元光谱矩阵M∈·^L×p，求解式(4)所示约束的非线性最小二乘问题

$\mathrm{\&theta;}=\underset{\mathrm{\&theta;}}{\arg \min }J\left(\mathrm{\&theta;}\right)=\underset{\mathrm{\&theta;}}{\arg \min }\left\{\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}\left(\mathrm{\&theta;}\right)s\left|\right|}^2\right\},s.t.\left(2\right)---\left(4\right)$

得到估计的丰度矢量s和非线性参数γ；进而估计θ=(s^T，γ^T)^T，实现高光谱解混；其中所述交替迭代优化算法具体步骤如下，其中γ^(t)和s^(t)分别表示第t迭代后的结果；

步骤1.初始化S和Г采用FCLS算法或者随机初始化得到初始的丰度矩阵

随机选取满足要求的初始非线性参数矩阵其中，

为第n个丰度矢量的初始值，

为第n个非线性参数的初始值，n=1，2，·，N；

步骤2.对于每一个观测像素x_n，n=1，2，·，N，执行以下操作：

2a)设置迭代次数t=0，x=x_n，

2b)依次求解优化问题

${\mathrm{\&gamma;}}^{(t+1)}=\underset{\mathrm{\&gamma;}}{\arg \min }{J}_1\left(\mathrm{\&gamma;}\right)=\underset{\mathrm{\&gamma;}}{\arg \min }\{\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}(s^{\left(t\right)},\mathrm{\&gamma;})s^{\left(t\right)}\left|\right|}^2+{\mathrm{\&lambda;}}_1{f}_B\left(\mathrm{\&gamma;}\right)\}---\left(5\right)$

和

$s^{(t+1)}=\underset{s}{\arg \min }{J}_2\left(s\right)=\underset{s}{\arg \min }\{\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}\left({\mathrm{\&theta;}}^{(t+1)}\right)s\left|\right|}^2+{\mathrm{\&lambda;}}_1{f}_{\mathrm{ASC}}\left(s\right)+{\mathrm{\&lambda;}}_2{f}_B\left(s\right)\}---\left(6\right)$

得到γ^(t+1)和s^(t+1)，其中θ^(t+1)=((s^(t))^T，(γ^(t+1)^T)^T，

假设γ_i或s_i所在的区间为[a，b]，则函数φ(x)定义为

$\mathrm{\&phi;}\left(x\right)=\left\{\begin{array}{cc}{(x-a)}^{2n}/2n,& x<a\\ 0,& a\&le;x\&le;b\\ {(x-b)}^{2n}/2n,& x>b\end{array}\right.,n\&Element;{\&CenterDot;}^{+}.---\left(7\right)$

n为正整数，λ₁和λ₂分别是表示所添加约束的强弱的正买数；

2c)判断目标函数J(θ)是否已收敛：若收敛输出x_n的丰度矢量s_n=s^(t+1)和非线性矢量γ_n=γ^(t+1)，并处理下一个观测像素，否则返回2b)继续执行；

步骤3.输出结果S=(s₁，s₂，·，s_N)和Г=(γ₁，γ₂，·，γ_N)。

2.根据权利要求1所述的估计方法，其特征在于：步骤2b)中，采用梯度下降法求解γ^(t+1)，其迭代步骤为

${\mathrm{\&gamma;}}^{(t+1)}={\mathrm{\&gamma;}}^{\left(t\right)}-{\mathrm{\&lambda;}}^{\left(t\right)}\frac{\&PartialD;J_1\left(\mathrm{\&gamma;}\right)}{\&PartialD;\mathrm{\&gamma;}}---\left(8\right)$

其中λ^(t)采用固定步长或线性搜索的方法求得。

3.根据权利要求1所述的估计方法，其特征在于：步骤2b)中，采用梯度法求解s^(t+1)。

4.根据权利要求3所述的估计方法，其特征在于：步骤2b)中，采用FCLS算法求解s^(t+1)。

5.根据权利要求1所述的估计方法，其特征在于：步骤2c)中，采用每次迭代与前一次迭代之间J(θ)的相对变化ξ，判断目标函数J(θ)是否收敛；ξ和J(θ)的计算具体公式如下：

$\mathrm{\&xi;}=\frac{J{\left(\mathrm{\&theta;}\right)}^{t+1}-J{\left(\mathrm{\&theta;}\right)}^t}{J{\left(\mathrm{\&theta;}\right)}^t},J\left(\mathrm{\&theta;}\right)=\frac{1}{2}{\left|\right|x-\stackrel{\&CenterDot;}{M}\left(\mathrm{\&theta;}\right)s\left|\right|}^2+{\mathrm{\&lambda;}}_1{f}_B\left(\mathrm{\&gamma;}\right)+{\mathrm{\&lambda;}}_1{f}_{\mathrm{ASC}}\left(s\right)+{\mathrm{\&lambda;}}_2{f}_B\left(s\right)---\left(9\right)$

6.根据权利要求2所述的估计方法，其特征在于：所述相对变化ξ小于10^-6时，目标函数J(θ)收敛。

Images