CN105844619A - 一种高光谱遥感影像端元提取方法 - Google Patents

Abstract

本发明公开了一种高光谱遥感影像端元提取方法，针对纯净像元假设不一定成立的情况，该方法将线性混合模型与粒子群优化算法相结合，配合高维粒子编码方式，推导出新的有效的粒子进化方程。同时，实验证明该方法针对端元提取问题更加有效。

Description

一种高光谱遥感影像端元提取方法

技术领域

本发明属于高光谱遥感影像处理技术领域，涉及一种高光谱遥感影像端元提取方法，尤其涉及一种线性混合模型约束的粒子群优化的高光谱遥感影像端元提取方法。

背景技术

高光谱遥感影像具有光谱分辨率高，图谱合一的特性。每个像元可以提供几乎连续的地物光谱曲线，使利用高光谱遥感反演陆地细节成为可能。高光谱遥感技术已经成为当前遥感领域的前沿技术之一。然而，由于高光谱传感器对每一个像元提供很窄的成像波段，这一波段具有较低的能量，使得只有扩大瞬时视场面积，传感器才能获取到有效信息。因此，高光谱传感器一般具有较低空间分辨率，在进行对地成像时，受自然界地物分布复杂多样性的影响，像元很难由单一均匀的地表覆盖类型组成，因此高光谱遥感影像中存在大量的混合像元。

混合像元分解技术是解决混合像元问题的有效途径，该方法主要包括端元提取和丰度估计两个核心任务。其中，端元提取是丰度估计的前提条件。针对混合像元端元提取问题，国内外学者提出了很多方法，传统的端元提取算法基于纯净像元存在假设。该类方法认为高光谱图像中每类地物至少存在一个纯净的像元。经典的端元提取算法包括：纯像元指数、N-FINDR、正交子空间投影、顶点成分分析等。

然而，由于成像条件的影响以及地物的复杂性，纯净像元存在的假设往往不成立。那么传统的基于纯净像元假设的算法就失效了，现有的一些端元生成算法发展还不成熟。因此，需要发展一种有效的端元提取算法，适用于纯净像元不存在的情况。

发明内容

为了解决上述技术问题，本发明提供了一种线性混合模型约束的粒子群优化算法用于高光谱遥感影像端元提取方法。

本发明所采用的技术方案为：一种高光谱遥感影像端元提取方法，其特征在于：将高光谱遥感影像分别读入大小为l×n的矩阵R＝{r₁,...,r_j,...,r_n}中，矩阵中每列向量r_j＝(r_1j,...,r_ij,...,r_lj)^Τ为像素的光谱辐射值，l为高光谱遥感影像的波段数，n为遥感影像的像素数；对高光谱遥感影像矩阵R进行如下操作：

步骤1：初始化种群，X_i(t)＝(e_i,1(t),e_i,2(t),...,e_i,D(t))i＝1,2,...,M；其中，t表示迭代的次数，初始化过程中令t＝0；X_i(t)表示种群中的一个粒子；M表示该种群中粒子的数目；e_ij(0)∈R，表示初始化过程中第i个粒子中的第j个向量；基于纯净像元存在假设的顶点成分分析VCA(Vertex Component Analysis)算法提取的结果作为初始的种群；

步骤2：计算每初始粒子的个体最优位置P_i(0)及所有粒子的全局最优位置G(0)；

P_i(0)＝X_i(0)

$G\left(0\right)=\underset{P_i}{\mathrm{argmin}}f\left(P_i\right(0\left)\right)$

$f\left(X_i\right(t\left)\right)=\mathrm{\γ}*vol\left(X_i\right(t\left)\right)+mean\left(abs\right(\underset{j=1}{\overset{D}{\mathrm{\Σ}}}{a}_{i,j}\left(t\right)-1\left)\right)+mean\left(\underset{j=1}{\overset{D}{\mathrm{\Σ}}}\right|\min \left(a_{i,j}\right(t),0)\left|\right)$

其中，f(·)表示适应度函数；vol(X_i(t))表示求解X_i(t)所构成的单形体的体积；丰度a_i,j(t)满足γ表示两个约束之间的比例系数。

步骤3：粒子进化；粒子进化方程为：

X_i,j(t+1)＝(p_i,j(t)-α_j(t)·(p_i(t)*β_j(t)))/(1-α_j(t))

其中，α_j(t)为一个随机生成数，表示图像的混合程度；β_j(t)是一个随机的丰度向量，需要满足的条件是：该向量中的第j个元素是0，且满足和为一条件；表示的是一个0～1的随机生成数；

步骤4：更新粒子的个体最优位置及所有粒子的全局最优位置；

粒子个体最优位置：

$P_i(t+1)=\left\{\begin{array}{cc}{P}_i\left(t\right)& \begin{array}{cc}if& f\left(\right\{X_i\left(t\right),X_{i,j}\left(t+1\right)\left\}\right)\&GreaterEqual;f\left(P_i\right(t\left)\right)\end{array}\\ \{X_i\left(t\right),X_{i,j}(t+1)\}& \begin{array}{cc}if& f\left(\right\{X_i\left(t\right),X_{i,j}\left(t+1\right)\left\}\right)<f\left(P_i\right(t\left)\right)\end{array}\end{array}\right.;$

其中，{X_i(t),X_i,j(t+1)}表示的是X_i(t)粒子的第j个列向量已经更新；f(x)表示适应度函数；

所有粒子的全局最优位置：

$G(t+1)=\underset{P_i}{\mathrm{argmin}}f\left(P_i\right(t+1\left)\right);$

步骤5：判断是否达到循环终止条件；

如果否，重新执行第3步；

如果是，则循环终止，输出最终的结果G(t)。

作为优选，步骤1中所述基于纯净像元存在假设的VCA算法，其具体实现过程是：(1)选择适当的超平面，将影像中的像元数据投影到选择的超平面上得到投影后的数据，后续的步骤在投影后的像元数据上操作；(2)将像元数据投影到随机方向，并将具有最大投影的像元作为第一个端元；(3)通过迭代地将数据投影到与由已提取端元构成的子空间正交的方向上的方法来提取其余端元，被极限投影所对应的像元作为新的端元，直到所有的像元被提取出来。

作为优选，步骤3中所述α_j(t)选择α_j(t)～U(-0.01,0.01)。

作为优选，步骤5中所述循环终止条件是X(t)收敛或达到最大迭代次数。

作为优选，所述方法还包括对最后的输出结果进行校正，校正公式为：

S_final＝(S_VCA*A_VCA ^T)*inv(A_VCA*A_VCA ^T)

$A_{VCA}=f_{unmixing}(S_{VCA},S_{LMMC-PSO})=inv(S_{LMMC-PSO}^T*S_{LMMC-PSO})*(S_{LMMC-PSO}^T*S_{VCA})$

其中，S_VCA∈R^l*D为VCA方法提取的端元结果，认为是混合像元；S_LMMC-PSO∈R^l*D为线性混合模型约束的粒子群优化算法端元提取的结果，认为是纯净像元；f_unmixing是丰度估计的函数，选择无约束的最小二乘方法；A_VCA是采用纯净像元S_LMMC-PSO求解混合像元S_VCA得到的丰度；校正之后的端元光谱曲线S_final为最终的结果。

本发明具有以下优点和积极效果：

(1)创造性的提出线性混合模型约束的粒子群优化算法，该方法没有速度更新公式，根据端元提取问题，利用线性混合模型更新粒子位置。

(2)保持了高光谱图像光谱特征物理意义。

(3)本发明模型简单，参数少，易于实现。

附图说明

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

图2是本发明实施例的实验结果示意图。

具体实施方式

为了便于本领域普通技术人员理解和实施本发明，下面结合附图及实施例对本发明作进一步的详细描述，应当理解，此处所描述的实施示例仅用于说明和解释本发明，并不用于限定本发明。

本发明需要解决的问题包括：(1)在纯净像元不存在的情况下，粒子群优化算法不能直接应用于端元提取问题；(2)高光谱图像光谱特征物理意义的保持。

本发明关键发明点为创造性提出线性混合模型约束的粒子群优化算法，并由于解决纯净像元不存在情况下的端元提取问题。

请见图1，本发明提供的一种高光谱遥感影像端元提取方法，将高光谱遥感影像分别读入大小为l×n的矩阵R＝{r₁,...,r_j,...,r_n}中，矩阵中每列向量r_j＝(r_1j,...,r_ij,...,r_lj)^Τ为像素的光谱辐射值，l为高光谱遥感影像的波段数，n为遥感影像的像素数；对高光谱遥感影像矩阵R进行如下操作：

步骤1：初始化种群，X_i(t)＝(e_i,1(t),e_i,2(t),...,e_i,D(t))i＝1,2,...,M；其中，t表示迭代的次数，初始化过程中令t＝0；X_i(t)表示种群中的一个粒子；M表示该种群中粒子的数目；e_ij(0)∈R，表示初始化过程中第i个粒子中的第j个向量；基于纯净像元存在假设的顶点成分分析VCA(Vertex Component Analysis)算法提取的结果作为初始的种群；

VCA算法是比较成熟的、效果较好的端元提取算法，简单的操作步骤为：(1)选择适当的超平面，将影像中的像元数据投影到选择的超平面上得到投影后的数据，后续的步骤在投影后的像元数据上操作；(2)将像元数据投影到随机方向，并将具有最大投影的像元作为第一个端元；(3)通过迭代地将数据投影到与由已提取端元构成的子空间正交的方向上的方法来提取其余端元，被极限投影所对应的像元作为新的端元，直到所有的像元被提取出来。

步骤2：计算每初始粒子的个体最优位置P_i(0)及所有粒子的全局最优位置G(0)；

P_i(0)＝X_i(0)

$G\left(0\right)=\underset{P_i}{\mathrm{argmin}}f\left(P_i\right(0\left)\right)$

$f\left(X_i\right(t\left)\right)=\mathrm{\&gamma;}*vol\left(X_i\right(t\left)\right)+mean\left(abs\right(\underset{j=1}{\overset{D}{\mathrm{\&Sigma;}}}{a}_{i,j}\left(t\right)-1\left)\right)+mean\left(\underset{j=1}{\overset{D}{\mathrm{\&Sigma;}}}\right|\min \left(a_{i,j}\right(t),0)\left|\right)$

其中，f(·)表示适应度函数；vol(X_i(t))表示求解X_i(t)所构成的单形体的体积；丰度a_i,j(t)满足γ表示两个约束之间的比例系数。

步骤3：粒子进化；粒子进化方程为：

X_i,j(t+1)＝(p_i,j(t)-α_j(t)·(p_i(t)*β_j(t)))/(1-α_j(t))

其中，α_j(t)为一个随机生成数，表示图像的混合程度，在算法中选择α_j(t)～U(-0.01,0.01)；β_j(t)是一个随机的丰度向量，需要满足的条件是：该向量中的第j个元素是0，且满足和为一条件；表示的是一个0～1的随机生成数；

步骤4：更新粒子的个体最优位置及所有粒子的全局最优位置；

由于本发明所用的高维编码方式，粒子的每个列向量都是一个端元光谱曲线，每个列向量的进化都会引起适应度函数的变化。为了避免每个列向量同时变化带来的叠加结果的影响，在本发明中，采用每个粒子的每个列向量更新的同时更新粒子的个体最优以及所有粒子的全局最优位置。

粒子个体最优位置：

$P_i(t+1)=\left\{\begin{array}{cc}{P}_i\left(t\right)& \begin{array}{cc}if& f\left(\right\{X_i\left(t\right),X_{i,j}\left(t+1\right)\left\}\right)\&GreaterEqual;f\left(P_i\right(t\left)\right)\end{array}\\ \{X_i\left(t\right),X_{i,j}(t+1)\}& \begin{array}{cc}if& f\left(\right\{X_i\left(t\right),X_{i,j}\left(t+1\right)\left\}\right)<f\left(P_i\right(t\left)\right)\end{array}\end{array}\right.;$

其中，{X_i(t),X_i,j(t+1)}表示的是X_i(t)粒子的第j个列向量已经更新；f(x)表示适应度函数，步骤2已经给出定义；

所有粒子的全局最优位置：

$G(t+1)=\underset{P_i}{\mathrm{argmin}}f\left(P_i\right(t+1\left)\right);$

步骤5：判断X(t)是否收敛或达到最大迭代次数；

如果否，重新执行第3步；

如果是，则循环终止，输出最终的结果G(t)。

在上述过程中，无法保证高光谱图像光谱特征物理意义，因此会出现光谱曲线光滑性很差的现象。后处理主要是是根据高光谱成像的物理机制，也就是线性混合模型对上述的结果进行校正，校正公式为：

S_final＝(S_VCA*A_VCA ^T)*inv(A_VCA*A_VCA ^T)

$A_{VCA}=f_{unmixing}(S_{VCA},S_{LMMC-PSO})=inv(S_{LMMC-PSO}^T*S_{LMMC-PSO})*(S_{LMMC-PSO}^T*S_{VCA})$

其中，S_VCA∈R^l*D为VCA方法提取的端元结果，认为是混合像元；S_LMMC-PSO∈R^l*D为线性混合模型约束的粒子群优化算法端元提取的结果，认为是纯净像元；f_unmixing是丰度估计的函数，选择无约束的最小二乘方法；A_VCA是采用纯净像元S_LMMC-PSO求解混合像元S_VCA得到的丰度；校正之后的端元光谱曲线S_final为最终的结果。

校正之后的端元光谱曲线就是最终的结果。

基于矩阵R，对高光谱遥感影像进行端元提取，在具体实现的时候还有以下注意事项：

(1)粒子的定义：采用高维编码方式，每个粒子为一个矩阵而不是一个向量，粒子中的每个列向量代表一条待定的端元光谱曲线，这样的定义方式更符合问题的物理意义。

(2)自适应度函数f的选择。自适应度函数包括两个部分：一个部分是体积，负责粒子向着单形体的内内部移动，避免单形体的体积过大，另外一个部分是丰度约束，负责粒子向着单形体的外部移动，避免单形体的体积过小。

$f\left(X_i\right(t\left)\right)=\mathrm{\&gamma;}*vol\left(X_i\right(t\left)\right)+mean\left(abs\right(\underset{j=1}{\overset{D}{\mathrm{\&Sigma;}}}{a}_{i,j}\left(t\right)-1\left)\right)+mean\left(\underset{j=1}{\overset{D}{\mathrm{\&Sigma;}}}\right|\min \left(a_{i,j}\right(t),0)\left|\right);$

其中，a_i,j为解混丰度，γ为比例系数，两个约束条件的比例。

以下通过对比试验来验证本发明的有益效果。

本试验采用的数据为仿真数据，共224个波段，光谱覆盖范围为：0.38-2.5μm，影像尺寸64像素×64像素。分别采用VCA、MVC-NMF、MVSA、MVES、SISAL和本发明方法进行端元提取。

试验采用的定量评价指标为光谱角距离，比较所提取的端元光谱与真实参考光谱的光谱角距离，光谱角距离越小，说明该方法的效果越好。每个方法运行10次，取统计结果，实验结果见附图2。从试验结果上看，我们的方法效果明显优于对比算法。

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

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

Claims (5)

1.一种高光谱遥感影像端元提取方法，其特征在于：将高光谱遥感影像分别读入大小为l×n的矩阵R＝{r₁,...,r_j,...,r_n}中，矩阵中每列向量r_j＝(r_1j,...,r_ij,...,r_lj)^Τ为像素的光谱辐射值，l为高光谱遥感影像的波段数，n为遥感影像的像素数；对高光谱遥感影像矩阵R进行如下操作：

步骤1：初始化种群，X_i(t)＝(e_i,1(t),e_i,2(t),...,e_i,D(t))i＝1,2,...,M；其中，t表示迭代的次数，初始化过程中令t＝0；X_i(t)表示种群中的一个粒子；M表示该种群中粒子的数目；e_ij(0)∈R，表示初始化过程中第i个粒子中的第j个向量；基于纯净像元存在假设的顶点成分分析VCA算法提取的结果作为初始的种群；

步骤2：计算每初始粒子的个体最优位置P_i(0)及所有粒子的全局最优位置G(0)；

P_i(0)＝X_i(0)

$G\left(0\right)=\underset{P_i}{\mathrm{argmin}}f\left(P_i\right(0\left)\right)$

$f\left(X_i\right(t\left)\right)=\mathrm{\&gamma;}*vol\left(X_i\right(t\left)\right)+mean\left(abs\right(\underset{j=1}{\overset{D}{\&Sigma;}}{a}_{i,j}\left(t\right)-1\left)\right)+mean\left(\underset{j=1}{\overset{D}{\&Sigma;}}\right|min\left(a_{i,j}\left(t\right),0\right)\left|\right)$

其中，f(·)表示适应度函数；vol(X_i(t))表示求解X_i(t)所构成的单形体的体积；丰度a_i,j(t)满足γ表示两个约束之间的比例系数；

步骤3：粒子进化；粒子进化方程为：

X_i,j(t+1)＝(p_i,j(t)-α_j(t)·(p_i(t)*β_j(t)))/(1-α_j(t))

其中，α_j(t)为一个随机生成数，表示图像的混合程度；β_j(t)是一个随机的丰度向量，需要满足的条件是：该向量中的第j个元素是0，且满足和为一条件；表示的是一个0～1的随机生成数；

步骤4：更新粒子的个体最优位置及所有粒子的全局最优位置；

粒子个体最优位置：

$P_i(t+1)=\left\{\begin{array}{cc}{P}_i\left(t\right)& iff\left(\right\{X_i\left(t\right),X_{i,j}(t+1)\left\}\right)\&GreaterEqual;f\left(P_i\right(t))\\ \{X_i\left(t\right),X_{i,j}(t+1)\}& iff\left(\right\{X_i\left(t\right),X_{i,j}(t+1)\left\}\right)<f\left(P_i\right(t))\end{array}\right.;$ 其中，{X_i(t),X_i,j(t+1)}表示的是X_i(t)粒子的第j个列向量已经更新；f(x)表示适应度函数；

所有粒子的全局最优位置：

$G(t+1)=\underset{P_i}{\arg \min }f\left(P_i\right(t+1\left)\right);$

步骤5：判断是否达到循环终止条件；

如果否，重新执行第3步；

如果是，则循环终止，输出最终的结果G(t)。

2.根据权利要求1所述的高光谱遥感影像端元提取方法，其特征在于：步骤1中所述基于纯净像元存在假设的VCA算法，其具体实现过程是：(1)选择适当的超平面，将影像中的像元数据投影到选择的超平面上得到投影后的数据，后续的步骤在投影后的像元数据上操作；(2)将像元数据投影到随机方向，并将具有最大投影的像元作为第一个端元；(3)通过迭代地将数据投影到与由已提取端元构成的子空间正交的方向上的方法来提取其余端元，被极限投影所对应的像元作为新的端元，直到所有的像元被提取出来。

3.根据权利要求1所述的高光谱遥感影像端元提取方法，其特征在于：步骤3中所述α_j(t)选择α_j(t)～U(-0.01,0.01)。

4.根据权利要求1所述的高光谱遥感影像端元提取方法，其特征在于：步骤5中所述循环终止条件是X(t)收敛或达到最大迭代次数。

5.根据权利要求1所述的高光谱遥感影像端元提取方法，其特征在于：所述方法还包括对最后的输出结果进行校正，校正公式为：

S_final＝(S_VCA*A_VCA ^T)*inv(A_VCA*A_VCA ^T)

$A_{\mathrm{VCA}}=f_{\mathrm{unmixing}}(S_{\mathrm{VCA}},S_{\mathrm{LMMC}-\mathrm{PSO}})=\mathrm{inv}(S_{\mathrm{LMMC}-\mathrm{PSO}}^T*S_{\mathrm{LMMC}-\mathrm{PSO}})*(S_{\mathrm{LMMC}-\mathrm{PSO}}^T*S_{\mathrm{VCA}})$

其中，S_VCA∈R^l*D为VCA方法提取的端元结果，认为是混合像元；S_LMMC-PSO∈R^l*D为线性混合模型约束的粒子群优化算法端元提取的结果，认为是纯净像元；f_unmixing是丰度估计的函数，选择无约束的最小二乘方法；A_VCA是采用纯净像元S_LMMC-PSO求解混合像元S_VCA得到的丰度；校正之后的端元光谱曲线S_final为最终的结果。