CN103530875B - 一种端元提取数据预处理方法 - Google Patents

Abstract

本发明提供了一种端元提取数据预处理方法，该方法通过建立基于模糊特征空间核空间引力模型，为高光谱数据像素点定义像元变异指数，实现了高光谱图像变异像素点的检测和移除。首先对高光谱数据集合进行模糊特征空间变换，得到模糊特征。利用像素的模糊特征，运用高斯径向基核函数优化的空间引力模型，计算3×3的空间邻域窗口内，邻域像素对中心像素的累加引力值，该引力值与像素变异指数成反比，对高变异指数像素进行移除。

Description

一种端元提取数据预处理方法

技术领域

本发明涉及数据处理领域，尤其涉及一种端元提取数据预处理方法。

背景技术

高光谱遥感技术已经在精细农业、环境检测、地质勘探、军事侦查等领域显现巨大的优势。然而，尽管高光谱遥感图像具有较高的光谱分辨率，却始终受空间分辨率的制约。较低的空间分辨率导致混合像元在高光谱数据中大量存在。近20年，混合像元分析成为该领域研究学者不容忽视的课题。

在线型混合模型假设下，若当前像素光谱曲线，有且仅由一种地物光谱构成，定义该像素点为对应于该地物的端元；反之，若当前像素光谱曲线由多种地物光谱线性叠加构成，定义该像素点为混合像素。传统的基于凸面单形体几何理论的端元提取方法，如：最大单形体体积(N-FINDR)、顶点成分分析(VCA)、正交子空间投影(OSP)及像素纯度指数(PPI)等诸多方法，旨在为每类地物寻找其对应的端元像素点。这种单一的端元提取方法——每种地物仅由某一条光谱曲线端元表示——忽略了像素光谱具有变异性的固有特性。

受光照条件、作物长势、含水量及阴影等条件的影响，像元光谱表现固有的变异性。C.A.Bateson，Gregory P.Asnerc，Roberts等学者针对光谱变异性展开了长期较深入的研究，以端元束集合替代单一端元集合，为待解混像素反演出最小均方误差准则下的丰度估计。这种做法尽管指出了光谱具有变异性的特点，却仍未从根本上解决变异光谱对真实端元光谱在光谱特征空间的干扰。

目前国内外尚无端元提取数据预处理方面的算法研究。

发明内容

本发明提供了一种端元提取数据预处理方法，用于检测并移除高光谱数据中具有较强光谱变异性的像素点。

一方面，提供了一种端元提取数据预处理方法，包括：对高光谱数据矩阵X进行模糊K-mean聚类，得到每个像元在模糊特征空间的模糊属性，以该模糊属性取代光谱特征；其次，进行3×3邻域窗口内核空间引力模型的高变异指数像素检测和移除；最后，设定判定阈值因子α，移除高变异指数像素，保留稳定像素点。模糊特征空间核空间引力模型的应用，实现了空间信息和光谱模糊特征信息的一体化，具有在高维非线性空间中从高光谱数据集中分离出高变异指数像素的能力。

优选的，对高光谱数据进行模糊K-mean聚类，以聚类后像素点的模糊特征向量替代原始光谱特征向量。其中，聚类数目等于由虚拟维数法估计出的端元数目。实现了光谱数据矩 阵由原始光谱特征空间到模糊特征空间的变换，将高维的光谱特征数据转换为较低维度的模糊特征数据。

可选的，运用高斯径向基核函数优化空间引力模型。模型的加权因子采用像素点距离加权因子，遍历计算全局像素点，为像素点定义变异指数。

优选的，设定判定阈值因子α，其中α非定值，因数据自身特性而变化，该因子需根据实际需要做适当调整。同时，该发明指出：判定阈值因子α可在均值引力误差附近范围取值。

本发明将高光谱数据的光谱和空间特性一体化，实现了对象元变异指数的判定，是一种高光谱端元提取数据预处理方法。该发明可作为现有的基于凸面单形体几何理论的端元提取算法的前期数据预处理方法，降低了变异像素对端元提取算法的干扰，提升了现有基于凸面单形体几何理论的端元提取算法的精度，具有较高的普适性。

附图说明

图1是高变异指数像素在光谱特征空间对真实端元造成淹没现象，单形体体积被错误估计示意图；

图2是高变异指数像素在模糊特征空间中的分布示意图；

图3是高斯径向基核函数优化的空间引力模型检测窗口示意图；

图4(a)是三种玉米种子图片；

图4(b)是玉米类型1、玉米类型2、玉米类型3对应的光谱特征曲线图；

图4(c)是叶片图像；

图4(d)是叶片光谱特性曲线图；

图5(a)是模拟数据玉米类型1丰度分布图；

图5(b)是模拟数据玉米类型2丰度分布图；

图5(c)是模拟数据玉米类型3丰度分布图；

图6(a)是高光谱AVIRIS印第安纳州农林数据第10波段图像；

图6(b)是印第安纳州农林数据16种地物类别的分布参照图；

图6(c)是印第安纳州农林数据16种地物对应灰度级的地物名称；

图7(a)是印第安纳州农林数据A(corn-notill)、B(corn-mino)、C(corn)的空间位置图；

图7(b)是corn-notill、corn-mino、corn的光谱特性曲线图；

图7(c)是N-FINDR算法在原始数据中端元提取结果图；

图7(d)是N-FINDR算法在FFS-KSAM预处理后的端元提取结果图；

图8是根据本发明实施例的端元提取数据预处理方法流程图。

具体实施方式

下面对本发明实施例涉及的定义进行说明：

像素点r_i的变异性大小由变异指数SI_i衡量，若某像素点具有较强的光谱变异性，并在光谱特征空间淹没了真实的端元，该象元被定义为具有较高变异指数的像素。

本发明所述的方法可以被定义为：模糊特征空间核空间引力模型（Fuzzy Feature Space Kernel Spatial Attraction Model,FFS-KSAM）。本发明首先对高光谱数据进行模糊特征空间变换，在模糊特征空间运用高斯径向基核函数优化的空间引力模型，在3×3的空间邻域窗口内对中心像素进行检测并判定该像元是否属于高变异指数像素。

线性混合模型假设：高光谱数据矩阵用X表示，其大小为L×N，其中L为波段总数，N为像元总数。第i个混合像素r_i的光谱特性为端元的线性叠加，用(1)式表达如下：

$r_i=\underset{j=1}{\overset{T}{\mathrm{\Σ}}}{a}_{i,j}{e}_j+n_i---\left(1\right)$

式中，e_j表示第j个端元对应的光谱向量，该向量大小为L×1，a_i,j为混合像素r_i中端元e_j的丰度系数，其中a_i,j满足约束条件：(a)a_i,j≥0（“非负性”约束）；(b)（和为“1”约束），端元总数为T，n_i为当前像素的噪声向量。

本发明实施例所述的端元提取数据预处理方法包括：

步骤一、对数据矩阵X进行模糊K-mean聚类，寻求代价函数J(Q,C)的极值，

$J(Q,C)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}\underset{j=1}{\overset{c}{\mathrm{\Σ}}}{\left(q_{i,j}\right)}^{\mathrm{\τ}}{\left({\mathrm{\ρ}}_{i,j}\right)}^2---\left(2\right)$

式中Q=[q_i,j](i=1,2,…,N;j=1,2,…,c)为模糊聚类的概率矩阵，N为像素总个数，c为聚类的类别数目，c值与端元总数相对应，本发明中，c值可以根据虚拟维数法进行估计，q_i,j满足约束条件：(a)0≤q_i,j≤1；(b)第i类的类别中心标记为c_i，聚类中心矩阵为C=[c_i]((i=1,2,…,c))，ρ_i,j根据(3)式计算，

(ρ_i,j)²=||r_i-c_j||²=(r_i-c_j)^T(r_i-c_j) (3)

经模糊聚类后，用其模糊概率特征替换原始像素向量，位于光谱空间的光谱特征向量r_i在模糊概率特征空间对应于模糊特征向量R_i，表示为(4)式，

R_i=[q_i,1,q_i,2,q_i,3,…,q_i,j,…,q_i,c]^T (4)

q_i,j表示r_i属于类别j的模糊概率隶属度。

步骤二、根据步骤一的计算结果，在模糊特征空间，利用高斯径向基核函数空间引力模型，遍历计算3×3窗口内，中心像素点的变异指数。具体理论如下：

设检测窗口大小为m×m。（注：考虑理论上的简易性，本发明设计检测窗口为方形检测窗）模糊特征空间，当前检测窗口的中心像素点向量记为R_i，检测窗口内中心像素点的邻域像素向量记为R_l(l=1,2,…,i-1,i+1,…,m²-1)。根据(5)式，计算R_i和R_l的空间引力，该引力大小记为h_i,j，其物理意义为邻域像素点R_l对中心像素点R_i的空间引力，

上式中，V_n(r_i)和V_n(r_l)为原始光谱空间像素r_i及r_l属于类别n的概率，为相应的加权因子。需要指明的是：在本发明的理论分析中，若未特殊说明，原始空间的像素点r_i和r_l分别对应模糊特征空间向量R_i和R_l。

用如下公式(6)的高斯径向基核函数替代(5)式中的内积计算，对低维模糊概率特征空间进行非线性的高维特征空间映射，并按照公式(7)计算当前检测窗口内邻域像素对中心像素R_i的累加核空间引力，

K(R_i,R_l)=exp(-||R_i-R_l||²/δ²)=<R_i,R_l> (6)

R_i的累加核空间引力大小记为H_i，用像素r_i和r_l的空间距离加权因子替代加权因子的计算，t为距离指数因子，本发明取t=2，该指数可做适当调整，应满足t≥0的原则。高斯径向基核函数的参数δ根据公式(8)定义，

$\mathrm{\&delta;}=\frac{1}{c\&times;N}\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{q}_{i,j}|r_i-c_i|---\left(8\right)$

考虑极限情况下，检测窗口内像素点光谱特性相同，累加核空间引力H_i达到极值，该极值记为H_max，

步骤三、遍历计算像素点r_i的变异指数SI_i，SI_i定义如下，

${\mathrm{SI}}_{\left(i\right)}=H_{\max }-H_i=\left\{\begin{array}{cc}1& {\mathrm{SI}}_{\left(i\right)}\&GreaterEqual;\mathrm{\&alpha;}\\ 0& {\mathrm{SI}}_{\left(i\right)}<\mathrm{\&alpha;}\end{array}\right.---\left(10\right)$

α为判定门限因子，该判定因子α随待处理高光谱数据不同而略有变化，一般情况下可在均值引力误差H_mean附近范围取值，计算方法如(11)式所示，

$\mathrm{\&alpha;}\&ap;H_{\mathrm{mean}}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(H_{\max }-H_i)---\left(11\right)$

若经计算，变异指数SI_i为1，则该像素点视为高变异指数像素点，被移除，不参加端元提取算法的计算；反之，若变异指数SI_i为0，则该像素点视为稳定像素点，被保留，参加后续端元提取算法计算。所有像素点判定完毕，预处理过程结束。

本发明实施例提出的FFS-KSAM方法，属高光谱图像变异像素点检测和移除的预处理技术，是高光谱遥感图像混合像元分析领域研究前沿技术，国内外尚无该方面研究。该发明提供了服务于端元提取的前期数据预处理方案，可与任意现有的基于凸体几何理论的端元提取技术相结合，抑制了变异像素点对真实端元的淹没现象。该方法通过建立模糊特征空间核空间引力模型，实现了高光谱图像变异像素点检测和移除，为高光谱图像混合像元分析领域传统的基于凸面单形体几何理论的端元提取算法提供了海量数据的前期预处理方案，提升了端元提取的准确性。

本发明实施例的优点包括：

1)基于FFS-KSAM的端元提取数据预处理方法，可与任意现有的基于凸体几何理论的端元提取技术相结合，具有较高的可移植性。

2)FFS-KSAM方法，将光谱特征信息与空间邻域信息一体化，强化了对变异像素点的检测、识别能力；

3)FFS-KSAM方法预处理后的数据，抑制了变异像素点对真实端元的淹没现象，提升了端元提取的精度；

4)FFS-KSAM方法在一定程度上对冗余的高光谱数据进行了精简，降低了端元提取的计算次数，缩短了计算时间。

下面结合附图举例对本发明做更详细地描述。

如图8所示，本发明实施例还提供了一种端元提取数据预处理方法，在该方法中，首先，读入高光谱数据，根据虚拟维数法对端元数目进行估计，估计出模糊特征后，进行空间变换，然后，进行核函数空间引力模型检测识别，判断是否检测出异常像素点，如果否，则保留数据，否则移除异常像素点。移除之后，判断检测识别步骤是否已经遍历所有像素，如果否，则继续进行核函数空间引力模型检测识别，否则，流程结束。

参照图1，是高变异指数像素在光谱特征空间对真实端元造成淹没现象，单形体体积被错误估计示意图。该图显示了：在二维光谱特征空间中，地物类别A存在较明显的高变异指数像素点C'，在光谱特征空间淹没了地物类别C的真实端元像素C，单形体S'被错误估计。该图显示了，高光谱数据中，高变异指数像素在常规基于凸体几何的端元提取算法中的干扰 现象。

参照图2，是模糊特征空间中高变异指数像素的分布示意图。该图显示了，经模糊特征空间变换后，数据存在于以地物类别概率为坐标轴的概率单形体T中，其中，具有高变异指数的像素位于概率单形体T的中心体积内部，成为端元概率较高的像素点分布于单形体T的顶点部分。该图展示了，经模糊特征空间变换，高变异指数像素点与稳定像素点的可分离现象。

参照图3，是高斯径向基核函数优化的空间引力模型检测窗口示意图。以3×3为检测窗口，像素R_i为当前检测窗口中心像素，R_l为邻域像素，为空间距离加权因子。该检测窗口利用高斯径向基核函数空间引力模型的高维特征非线性映射，根据模糊特征，检测并移除高变异指数像素，实现了光谱信息和空间信息的一体化。

参照图4(a)，是运用推扫式PIKA II型高光谱成像仪采集的玉米种子图像。每列为5粒且每列为同种遗传基因类型玉米种子，由于其遗传基因不同，从左至右分别定义为玉米类型1、玉米类型2、玉米类型3。图中总计15粒玉米种子，每列3粒，每列玉米种子具有相同的遗传基因，不同列玉米种子遗传基因不同，从左至右按遗传基因不同区分标记为玉米类型1、玉米类型2、玉米类型3。该图像光谱覆盖范围为405～907nm，可用波段数为160。

参照图4(b)，是三种转基因玉米种子——玉米类型1、玉米类型2、玉米类型3，对应的光谱特征曲线图；

参照图4(c)，是运用推扫式PIKA II型高光谱成像仪采集的叶片图像。（PIKA II型高光谱成像仪的参数同参照图4(a)）；

参照图4(d)，是叶片光谱特性曲线图；

利用参照图4的真实光谱特性曲线，制作大小100×100×160的高光谱模拟数据。端元光谱为玉米类型1、玉米类型2和玉米类型3。叶片光谱曲线不作为端元存在，以丰度0～0.02的大小线性混合在模拟高光谱数据中，旨在模拟高光谱数据像素光谱特征的变异性。

参照图5(a)，是摸拟数据中玉米类型1的丰度分布图。高亮矩形区域为高纯度光谱含量区，该区域空间大小为70×20像素，该区域内由玉米类型1的光谱，其丰度为0.98～1.00和叶片光谱，其丰度为0～0.02，线性混合而成。高亮度矩形区域外的丰度值大小与高亮矩形区域中心点距离成反比。

参照图5(b)，是摸拟数据中玉米类型2丰度分布图。高亮矩形区域为高纯度光谱含量区，该区域空间大小为70×20像素，该区域内由玉米类型2的光谱，其丰度为0.98～1.00和叶片光谱，其丰度为0～0.02，线性混合而成。高亮度矩形区域外的丰度值大小与高亮矩形区域中心点距离成反比。

参照图5(c)，是摸拟数据中玉米类型3丰度分布图。高亮矩形区域为高纯度光谱含量区， 该区域空间大小为70×20像素，该区域内由玉米类型3的光谱，其丰度为0.98～1.00和叶片光谱，其丰度为0～0.02，线性混合而成。高亮度矩形区域外的丰度值大小与高亮矩形区域中心点距离成反比。

需指明的是，在模拟数据中，高亮矩形区域分别为对应玉米类型光谱丰度为0.98～1.00，该部分区域模拟地物高纯度区域；亮矩形区域外部为三种类型玉米种子的线性混合，该部分区域模拟不同地物类别线性混合区域。全局叶片光谱的添加，丰度为0～0.02，旨在模拟光谱变异性。

表1是模拟数据端元提取算法得到的平均光谱角及丰度反演均方根误差，参照表1，是模拟数据端元提取算法得到的平均光谱角及丰度反演均方根误差；

表1

参照图6(a)是真实的高光谱AVIRIS数据图像，该区域为印第安纳州农林区数据第10波段图像，该数据于1992年6月采集，光谱覆盖范围为0.4～2.5um，去除低信噪比后，我们采用该数据中的100个波段。

参照图6(b)是印第安纳州农林数据16种地物类别的分布参照图，以不同的灰度级标识不同的地物类别；

参照图6(c)是印第安纳州农林数据16种地物对应灰度级的地物名称；

参照图7(a)是印第安纳州农林数据A(corn-notill)、B(corn-mino)、C(corn)的空间位置图，此三种地物空间分布近邻，地物属性相似，较易产生明显的光谱变异现象；

参照图7(b)是corn-notill、corn-mino、corn的光谱特性曲线图；

参照图7(c)是N-FINDR算法在原始数据中端元提取结果图，该图显示了未经预处理的高光谱数据，在主成分空间因受高变异指数像素的影响，导致了端元误提取现象；

参照图7(d)是N-FINDR算法在FFS-KSAM方法预处理后的高光谱数据进行端元提取结果图，该发明提出的FFS-KSAM方法有效地抑制了高变异像素点；

表2是印第安纳州农林数据16种地物光谱重建均方根误差及统计平均均方根误差，参照表2是印第安纳州农林数据光谱重建均方根误差（经FFS-KSAM方法预处理后，运用VCA、N-FINDR、OSP算法，判定门限因子分别为α=0.5×10^-4、α=0.7×10^-4及α=1.0×10^-4，16种地物对应的光谱重构均方根误差及平均光谱重构均方根误差）；

表2

表3是印第安纳州农林数据经FFS-KSAM预处理后，端元提取算法（VCA、N-FINDR、OSP）的运行时间，参照表3是印第安纳州农林数据经FFS-KSAM预处理后，端元提取算法（VCA、N-FINDR、OSP算法）的运行时间；

表3

上述为本发明特举之实例，并非用以限定本发明。本发明提供的“一种新型的端元提取数据预处理方法”，同样适用于对其他的高光谱图像进行前期预处理。在不脱离本发明的实质和范围内，可做些许的调整和优化，以本发明的保护范围以权利要求为准。

Claims (5)

1.一种端元提取数据预处理方法,其特征在于，包括：

获取高光谱图像；

对所述高光谱图像对应的高光谱数据矩阵进行模糊聚类，得到每个像元在模糊特征空间的模糊属性，以该模糊属性特征向量取代所述矩阵中的光谱特征向量；

对所述矩阵进行3×3邻域窗口内核空间引力模型的高变异指数像素检测和移除，移除高变异指数像素，保留稳定像素点；

所述的端元提取数据预处理方法包括：

(1)对数据矩阵X进行模糊K-mean聚类，寻求代价函数J(Q,C)的极值，

$J(Q,C)=\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}{\left(q_{i,j}\right)}^{\mathrm{\&tau;}}{\left({\mathrm{\&rho;}}_{i,j}\right)}^2$

式中Q＝[q_i,j]，i＝1,2,…,N；j＝1,2,…,c为模糊聚类的概率矩阵，N为像素总个数，c为聚类的类别数目，q_i,j满足约束条件：第i类的类别中心标记为c_i，聚类中心矩阵为C＝[c_i]，i＝1,2,…,c，

(ρ_i,j)²＝||r_i-c_j||²＝(r_i-c_j)^T(r_i-c_j)

经模糊聚类后，用其模糊概率特征替换原始像素向量，位于光谱空间的光谱特征向量r_i在模糊概率特征空间对应于模糊特征向量R_i，

R_i＝[q_i,1,q_i,2,q_i,3,…,q_i,j,…,q_i,c]^T

q_i,j表示r_i属于类别j的模糊概率隶属度；

(2)根据步骤(1)的计算结果，在模糊特征空间，利用高斯径向基核函数空间引力模型，遍历计算3×3窗口内，中心像素点的变异指数：

检测窗口大小为m×m，模糊特征空间，当前检测窗口的中心像素点向量记为R_i，检测窗口内中心像素点的邻域像素向量记为R_l，l＝1,2,…,i-1,i+1,…,m²-1，计算R_i和R_l的空间引力，引力大小记为h_i,j，其物理意义为邻域像素点R_l对中心像素点R_i的空间引力，

上式中，V_n(r_i)和V_n(r_l)为原始光谱空间像素r_i及r_l属于类别n的概率，为相应的加权因子，原始空间的像素点r_i和r_l分别对应模糊特征空间向量R_i和R_l；

用高斯径向基核函数替代内积计算，对低维模糊概率特征空间进行非线性的高维特征空间映射，计算当前检测窗口内邻域像素对中心像素R_i的累加核空间引力，

K(R_i,R_l)＝exp(-||R_i-R_l||²/δ²)＝＜R_i,R_l＞

R_i的累加核空间引力大小记为H_i，用像素r_i和r_l的空间距离加权因子替代加权因子的计算，t为距离指数因子，t＝2，t≥0，高斯径向基核函数的参数δ

$\mathrm{\&delta;}=\frac{1}{c\&times;N}\underset{j=1}{\overset{c}{\mathrm{\&Sigma;}}}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{q}_{i,j}\left|r_i-c_i\right|$

考虑极限情况下，检测窗口内像素点光谱特性相同，累加核空间引力H_i达到极值，该极值记为H_max，

(3)遍历计算像素点r_i的变异指数SI_i，SI_i定义如下，

${\mathrm{SI}}_{\left(i\right)}=H_{max}-H_i=\left\{\begin{array}{cc}1& {\mathrm{SI}}_{\left(i\right)}\&GreaterEqual;\mathrm{\&alpha;}\\ 0& {\mathrm{SI}}_{\left(i\right)}<\mathrm{\&alpha;}\end{array}\right.$

α为判定门限因子，

$\mathrm{\&alpha;}=H_{mean}=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}(H_{max}-H_i)$

若经计算，变异指数SI_i为1，则该像素点视为高变异指数像素点，被移除，不参加端元提取算法的计算；若变异指数SI_i为0，则该像素点视为稳定像素点，被保留，参加后续端元提取算法计算，所有像素点判定完毕，预处理过程结束。

2.根据权利要求1所述的方法，其特征在于，模糊聚类的数目等于由虚拟维数法估计出的端元数目。

3.根据权利要求1所述的方法，其特征在于：运用高斯径向基核函数优化空间引力模型,模型的加权因子采用像素点距离加权因子，遍历计算全局像素点，为像素点定义变异指数。

4.根据权利要求1所述的方法，其特征在于：设定判定阈值因子α，其中α根据数据自身特性进行变化。

5.根据权利要求1所述的方法，其特征在于：所述模糊聚类包括模糊K-mean聚类。