CN105550693A - 一种基于非负独立成分分析的布谷鸟搜索高光谱解混方法 - Google Patents

Abstract

本发明公开了一种结合非负独立成分分析方法和布谷鸟搜索优化方法的高光谱图像解混方法。本发明具体包括以下步骤：1、对观测数据进行白化预处理，去除二阶相关性，对于高光谱图像解混；2.在非负独立成分分析方法的基础上，引入丰度和为一约束项构建目标函数；3.利用步骤2中得到的参数最优值计算解混矩阵，估计丰度矩阵；4、利用非负最小二乘方法估计出端元光谱。仿真数据和真实高光谱数据的实验结果表明，本发明能有效地克服传统独立成分分析方法在解决高光谱解混问题时的局限性，获取良好的精度。

Description

一种基于非负独立成分分析的布谷鸟搜索高光谱解混方法

技术领域

本发明涉及高光谱图像的预处理手段领域，更具体的说，是涉及一种结合非负独立成分分析方法和布谷鸟搜索优化方法的无监督的高光谱图像解混方法。

背景技术

光谱解混作为高光谱图像的一种预处理手段,不仅是实现地物精确分类和识别的重要前提,而且是遥感技术向定量化深入发展的重要条件,它对发挥高光谱图像，在社会经济和国防建设方面的效益有着重要的现实意义。高光谱图像是由高光谱遥感成像仪捕获的包含空间和光谱信息的图像，能够帮助研究者们更好地提取地物信息。尽管高光谱图像的光谱分辨率很高，但因自然界地物复杂多样以及遥感成像仪空间分辨率有限，单像元点往往包含多种不同类别的地物，意味着“混合像元”广泛存在，使得传统的基于像元级光谱统计特性的分类方法无法直接使用，这不仅制约着地物识别的精度，而且成为高光谱遥感向定量化方向深入发展的主要障碍。因此，如何有效地解决混合像元问题，使遥感应用由像元级达到亚像元级，实现“光谱解混”，是高光谱图像处理技术面临的难题之一。“光谱解混”过程就是将混合像元分解为一些主要成分(端元)和每个成分所占比例(丰度)的过程。

线性混合模型和非线性混合模型是高光谱解混方法中的两个基本模型，相比于非线性混合模型，线性混合模型因其模型简单且物理意义明确而应用地更为广泛。高光谱线性混合模型忽略不同端元间的多次散射，把每个像元近似为不同端元光谱以其丰度为权重的线型组合。考虑到，高光谱数据为L个波段上的N个不同的像元，则某单像元点的观测值xR^L×1可以表示为

x＝As+e(1)

式中，x＝[x₁,x₂,L,x_L]R^L×1，L是波段数；A＝[a₁,a₂,L,a_P]R^L×P是端元光谱矩阵，a_iR^L×1代表第i个端元光谱，P是端元数；s＝[s₁,s₂,L,s_P]^TR^P×1是每个像元点的丰度向量；eR^L×1是误差项。其中，丰度向量s对应于每个端元所占比例，丰度必须满足丰度非负约束与丰度和为一约束，即

对于高光谱图像的N个不同像元所有观测向量组成观测矩阵XR^L×P，可以表示为

X＝AS+E(3)

式中，SR^P×N是N个像元的丰度向量组成的丰度矩阵，E对应得误差矩阵。现存的很多方法解混方法都是基于线性混合模型的，主要有三大类：基于几何原理的方法、基于稀疏回归的方法和基于统计方法的方法。基于几何方法的解混方法，如N-FINDR，顶点成分分析(vertexcomponentanalysis，VCA)，最小体积单形体分析(MinimumVolumeSimplexAnalysis，MVSA)等，以观测数据中存在“纯端元”为前提，把高光谱数据集构成单形体的顶点作为端元，往往实际地物分布并不能保证所有端元都有“纯”像元存在,故处理实际遥感数据时效果不理想。

基于稀疏回归的方法，假设观测像元能够用一组已知的“纯”的端元光谱线性表示(光谱库作为已知端元)，与整个光谱库的维数先比，每个像元点中的端元数很小，因此可以通过寻找该光谱库中能够表示混合像元的最佳光谱子集实现解混。然而，该解混结果很大程度上依赖合适光谱库的选取，很难获取与观测数据集相同条件下的光谱库，这是一个耗时且高代价的过程。

基于统计方法的解混方法是一种典型的无监督的方法，在不需要先验知识和人工干预的条件下，实现端元光谱的自动提取。实际地物的分布情况往往非常复杂，导致端元的无监督自动提取成为当前研究的一个难点。值得一提的是，随着盲信号分离技术的兴起和发展，为如何用无监督的方法更快更精确地实现高光谱盲解混提供了另外一条可能途径，己经成为当今遥感领域里的一个研究热点。独立分分析方法作为一种典型性的盲分离算法，未用到光谱的先验信息，对观测像元的纯度不做要求，具有很高的研究价值。

发明内容

本发明的目的在于提供基一种基于非负独立成分分析的布谷鸟搜索高光谱解混方法。

本发明一种结合非负独立成分分析方法和布谷鸟搜索优化方法的高光谱图像解混方法，具体包括以下步骤：

步骤一、对观测数据进行白化预处理，去除二阶相关性，对于高光谱图像解混，采用利用主成分分析方法(PrincipalComponentsAnalysis，PCA)对观测数据进行白化降维预处理，具体步骤如下：

(1)对原始观测数据不进行去均值处理，保证估计出的丰度源是非负的；(2)先求原始数据的自相关矩阵Corr＝XX^T/N，确保丰度向量间使相关的；(3)采用PCA计算出Corr的前P个较大特征值组成的对角矩阵D_P和对应的特征向量矩阵E_P，利用公式求出白化矩阵X'＝VXR^P×N，从而使得白化后的数据同时得到降维；

步骤二：

在非负独立成分分析方法的基础上，引入丰度和为一约束项构建目标函数为

式中，为权重；为非负独立成分分析方法的目标函数，其中，Y_＝min(Y,0)是矩阵Y的负的部分；为丰度和为一约束的目标函数，其中，λ是一个较大的正数，1R^P×1的全一列矢量；

利用基于Gives矩阵的QR分解理论，降低解混矩阵维数，将解混矩阵W转化为Gives矩阵乘积的形式如下：

W＝T_P1LΤ₂T₁

(2)

T₁＝T_1PT_1,P1LT₁₂,T₂＝T_2PT_2,P1LT₂₃,T_P1＝T_P1,P

式中，T_ij是P阶旋转矩阵，旋转角度θ[0,2π]；

采用布谷鸟搜索算法进行优化的迭代公式为

式中，M是参数维数，表示第k个巢在第t代的位置，表示点对点乘法，α₀取值0.01，是全局最优巢位置，Levy(step)～|step|^1β,0<β2表示随机步长step服从Levy分布；

采用布谷鸟搜索算法对目标函数进行优化具体步骤如下：

(1)初始化种群，设置巢个数为M，搜索空间维数为随机初始化巢位置计算每个巢的目标函数的适度值，记录当前的最优适度值f_min和最优位置best{1,2,L,M}；

(2)保留前一代最优巢位置其中t[1,T_max]，并按位置更新公式(10)对鸟巢位置进行更新，得到当前巢位置与上一代巢的位置进行比较，将适应值较好的巢位置作为当前最好位置，

(3)每个巢分配一个服从均匀分布的0-1的随机数rd与卵被宿主发现的概率p_a进行比较。若满足rdp_a，则将当前最好位置g_t保存；如果rd<p_a成立，则随机改变满足该条件的部分巢位置得到一组新巢，并与当前最优位置g_t比较，将适应值较好的位置取代较差的位置，得到一组新的最优位置

(4)计算(3)中得到的gg_t中最优巢位置和相应的最优适度值f_min，并进行判断。若满足迭代终止条件，则输出全局最优值和对应的全局最优位置，反之，则返回步骤二继续循环更新；一般情况下，迭代终止条件为：达到指定的最大迭代次数或规定的精度要求；本发明中采前者作为终止条件；

步骤三、利用步骤二中得到的参数最优值计算解混矩阵，估计丰度矩阵；

步骤四、利用非负最小二乘方法估计出端元光谱。

本发明针对高光谱图像中“混合像元”广泛存在、线性混合模型中丰度非负约束和丰度和为一约束的特点，在传统的非负独立成分分析方法的基础上，加入光谱线性混合模型中的丰度和为一约束项，构成有约束的非负独立成分分析方法，引入用布谷鸟群智能搜索优化方法，利用其全局随机搜索性对目标函数进行优化，避免陷入局部极值，且对初始化选择没有严格要求；在优化过程中，为了减少参数维数并缩小布谷鸟搜索方法群体搜索范围，利用矩阵的QR分解理论，将对分离矩阵的搜索转化为对一系列Gives矩阵的识别，从而减小了计算量。仿真数据和真实高光谱数据的实验结果表明，本发明能有效地克服传统独立成分分析方法在解决高光谱解混问题时的局限性，获取良好的精度。

附图说明

图1是本发明的基于非负独立成分分析的布谷鸟搜索高光谱解混方法的流程图。

具体实施方式

下面结合附图对本发明做进一步说明。

有约束的非负独立成分分析方法

独立成分分析方法是一个无监督的盲源分离方法，其目标是在仅给出观测数据X的情况下，根据独立性测度，通过寻找一种线性变换分离出源，即

Y＝WX＝US(4)

式中，Y＝[y₁,y₂,L,y_L]^TR^L×N是源S的估计，W是通过最大化Y成分间的独立性得到的正交分离矩阵。通常情况下，独立性测度采用峭度、负熵和互信息等。

非负独立成分分析方法将传统独立成分分析方法对独立性要求与源的非负性结合起来进行分离，使分离结果尽可能的独立且是非负，通过解相关和非负旋转结合就可以实现，为了保证估计信号非负，通过最小化如下目标函数实现：

式中，Y_＝min(Y,0)是矩阵Y的负的部分，若Y是非负的，J_NICA(Y)接近于零；否则，J_NICA(Y)是一个较大的正数。

线性混合模型中的丰度和为一约束条件要求，对于每一个观测像元，其各个端元丰度之和必须为一，其目标函数可以定义为

式中，λ是一个较大的正数，1R^P×1的全一列矢量，当丰度矩阵中每个像元的丰度都满足和为一约束时,J_ASC(Y)接近于零，否则，J_ASC(Y)是一个较大的正数。

结合式(5)和式(6)，可以得到本发明的目标函数是

从而，得到丰度矩阵的估计为

式中，X'为白化后的观测数据。式(7)中，第一项能够使分离出的丰度变量满足非负约束，第二项可以实现丰度和为一，这两者结合可以确保丰度包含在特定的凸集中，来适应与真实丰度的分布情况。

布谷鸟搜索方法

布谷鸟搜索(Cuckoosearch，CS)是一种近年来发展起来的自然启发方法，属于全局随机搜索方法，由剑桥大学Xin-SheYang和SuashDeb在2009年提出。CS方法通过对某些布谷鸟物种的寄生育雏(BroodParasitism)行为进行模拟，利用相关的Levy飞行搜索机制，来求解优化问题。CS方法参数少、简单，与粒子群方法(ParticleSwarmOptimization，PSO)、遗传方法(GeneticAlgorithm，GA)和差分进化方法(DifferentiaEvolutionary，DE)等其他群体优化方法相比，CS方法优化求解能力更佳，搜索精度更高。

实际上，每个布谷鸟随机的或是类似随机寻地找适合自己产卵的巢，为了简单模拟布谷鸟寻巢产卵进行搜索的方式，我们使用以下三个理想化的规则：

(1)每个布谷鸟每次产生一个卵，堆放在一个随机选择的巢中；

(2)每次随机选择的巢中具有最高品质适应度值的巢将转到下一代；

(3)可用的巢的数量是固定的，布谷鸟的卵被宿主发现的概率为p_a[0,1]，可以进一步近似，即被新巢(新的解决方案)替换的概率。

基于以上三个规则，布谷鸟搜索路径和位置更新公式如下：

式中，M是参数维数，表示第k个巢在第t代的位置，表示点对点乘法，α>0，是步长控制系数，与感兴趣问题的尺度规模Scale有关，一般情况下取α＝O(Scale/100)。在为了适应巢质量间的差异，提出α的取值为式中，α₀取值0.01，是全局最优巢位置，意味着当前巢与最优巢的差异直接影响着生成新巢的可能性。Levy(step)～|step|^1β,0<β2表示随机步长step服从Levy分布。将步长step简化为式中，v～N(0,1)，代入式(9)，可得到新的位置更新公式为

上式即为发明中最终采用的巢位置的更新公式。

总结

高光谱图像解混问题可以归结为目标函数式(7)所示的优化求解问题。本发明采用布谷鸟搜索方法对式(7)进行优化求解，得到解混矩阵W。对于P个端元，求解解混矩阵W为P×P维，求解参数维数大小是D＝P²，随着端元个数的增加，求解参数维数迅速增加，计算量也随之快速增大。为了减少参数维数并缩小布谷鸟搜索方法群体搜索范围，本发明利用矩阵QR分解理论，将对解混矩阵的搜索转化为对一系列Gives矩阵的识别，将参数维数降低至从而减小了计算量。

解混矩阵W转化为Gives矩阵乘积的形式如下：

W＝T_P1LΤ₂T₁

(11)

T₁＝T_1PT_1,P1LT₁₂,T₂＝T_2PT_2,P1LT₂₃,T_P1＝T_P1,P

式中，T_fq是P阶旋转矩阵，旋转角度为θ[0,2π]。式(11)表明，W为个Gives矩阵的乘积，每个Gives矩阵只有一个位置元素θ，实现将参数维数降低P² 维。

如图1所示，首先进行预处理：对观测数据进行白化处理，去除二阶相关性。但是，对于高光谱图像解混，数据维数比较大，大大增加了计算量，需要对数据进行降维。对传统白化过程从三方面进行改进：(1)对原始观测数据不进行去均值处理，保证估计出的丰度源是非负的；(2)先求原始数据的自相关矩阵Corr＝XX^T/N，确保丰度向量间使相关的；(3)利用主成分分析方法(PrincipalComponentsAnalysis，PCA)计算出Corr的前P个较大特征值组成的对角矩阵D_P和对应的特征向量矩阵E_P，利用公式求出白化矩阵，从而使得白化后的数据同时得到降维X'＝VXR^P×N，此时，式(4)可以变换为

Y＝WX'(12)

其次进行CS优化求解：对白化后的数据进行旋转，得到旋转矩阵W使所有的数据尽可能的落在正象限中，通过最小化非负约束来实现。本发明采用布谷鸟搜索方法对式(7)进行优化求解，得到解混矩阵W。CS方法的对式(7)进行优化的步骤可以概括为：

步骤1初始化种群，设置巢个数为M，搜索空间维数为随机初始化巢位置计算每个巢的目标函数的适度值，记录当前的最优适度值f_min和最优位置best{1,2,L,M}。

步骤2保留前一代最优巢位置其中t[1,T_max]，并按位置更新公式(10)对鸟巢位置进行更新，得到当前巢位置与上一代巢的位置进行比较，将适应值较好的巢位置作为当前最好位置，

步骤3每个巢分配一个服从均匀分布的0-1的随机数rd与卵被宿主发现的概率p_a进行比较。若满足rdp_a，则将当前最好位置g_t保存；如果rd<p_a成立，则随机改变满足该条件的部分巢位置得到一组新巢，并与当前最优位置g_t比较，将适应值较好的位置取代较差的位置，得到一组新的最优位置

步骤4计算步骤3中得到的gg_t中最优巢位置和相应的最优适度值f_min，并进行判断。若满足迭代终止条件，则输出全局最优值和对应的全局最优位置，反之，则返回步骤二继续循环更新。一般情况下，迭代终止条件为：达到指定的最大迭代次数或规定的精度要求。本发明中采前者作为终止条件。

最后是得到解混结果：将最优巢位置代入式(11)得到解混矩阵W，根据式(8)丰度估计

端元光谱可通过非负最小二乘方法求解。

Claims (1)

1.一种结合非负独立成分分析方法和布谷鸟搜索优化方法的高光谱图像解混方法，其特征是，具体包括以下步骤：

步骤一、对观测数据进行白化预处理，去除二阶相关性，对于高光谱图像解混，采用利用主成分分析方法对观测数据进行白化降维预处理，具体步骤如下：

(1)对原始观测数据不进行去均值处理，保证估计出的丰度源是非负的；(2)先求原始数据的自相关矩阵Corr＝XX^T/N，确保丰度向量间使相关的；(3)采用PCA计算出Corr的前P个较大特征值组成的对角矩阵D_P和对应的特征向量矩阵E_P，利用公式求出白化矩阵X'＝VXR^P×N，从而使得白化后的数据同时得到降维；

步骤二：

在非负独立成分分析方法的基础上，引入丰度和为一约束项构建目标函数为

式中，为权重；为非负独立成分分析方法的目标函数，其中，Y_＝min(Y,0)是矩阵Y的负的部分；为丰度和为一约束的目标函数，其中，λ是一个较大的正数，1R^P×1的全一列矢量；

利用基于Gives矩阵的QR分解理论，降低解混矩阵维数，将解混矩阵W转化为Gives矩阵乘积的形式如下：

W＝T_P1LΤ₂T₁

(2)

T₁＝T_1PT_1,P1LT₁₂,T₂＝T_2PT_2,P1LT₂₃,T_P1＝T_P1,P

式中，T_ij是P阶旋转矩阵，旋转角度θ[0,2π]；

采用布谷鸟搜索算法进行优化的迭代公式为

${\mathrm{pos}}_k^{\left(t\right)}={\mathrm{pos}}_k^{(t-1)}+{\mathrm{\&alpha;}}_0\left(\begin{array}{cc}{\mathrm{pos}}_k^{(t-1)}& {\mathrm{pos}}_b^{(t-1)}\end{array}\right)*Levy\left(step\right),k=1,2,3,L,M,b\{1,2,L,M\}---\left(3\right)$

式中，M是参数维数，表示第k个巢在第t代的位置，表示点对点乘法，α₀取值0.01，是全局最优巢位置，Levy(step)～|step|^1β,0<β2表示随机步长step服从Levy分布；

采用布谷鸟搜索算法对目标函数进行优化具体步骤如下：

(1)初始化种群，设置巢个数为M，搜索空间维数为随机初始化巢位置计算每个巢的目标函数的适度值，记录当前的最优适度值f_min和最优位置best{1,2,L,M}；

(2)保留前一代最优巢位置其中t[1,T_max]，并按位置更新公式(10)对鸟巢位置进行更新，得到当前巢位置 ${\mathrm{pos}}_t={\&lsqb;{\mathrm{\&theta;}}_1^{\left(t\right)},{\mathrm{\&theta;}}_2^{\left(t\right)},L,{\mathrm{\&theta;}}_M^{\left(t\right)}\&rsqb;}^T,$ 与上一代巢的位置 ${\mathrm{pos}}_{t1}={\&lsqb;{\mathrm{\&theta;}}_1^{\left(t1\right)},{\mathrm{\&theta;}}_2^{\left(t1\right)},L,{\mathrm{\&theta;}}_M^{\left(t1\right)}\&rsqb;}^T$ 进行比较，将适应值较好的巢位置作为当前最好位置， $g_t={\&lsqb;{\mathrm{\&theta;}}_1^{\left(t\right)},{\mathrm{\&theta;}}_2^{\left(t\right)},L,{\mathrm{\&theta;}}_M^{\left(t\right)}\&rsqb;}^T;$

(3)每个巢分配一个服从均匀分布的0-1的随机数rd与卵被宿主发现的概率p_a进行比较。若满足rdp_a，则将当前最好位置g_t保存；如果rd<p_a成立，则随机改变满足该条件的部分巢位置得到一组新巢，并与当前最优位置g_t比较，将适应值较好的位置取代较差的位置，得到一组新的最优位置 ${\mathrm{pos}}_t={\&lsqb;{\mathrm{\&theta;}}_1^{\left(t\right)},{\mathrm{\&theta;}}_2^{\left(t\right)},L,{\mathrm{\&theta;}}_M^{\left(t\right)}\&rsqb;}^T;$

(4)计算(3)中得到的gg_t中最优巢位置和相应的最优适度值f_min，并进行判断。若满足迭代终止条件，则输出全局最优值和对应的全局最优位置，反之，则返回步骤二继续循环更新；迭代终止条件为：达到指定的最大迭代次数；

步骤三、利用步骤二中得到的参数最优值计算解混矩阵，估计丰度矩阵；

步骤四、利用非负最小二乘方法估计出端元光谱。