基于光谱解混的城市地物分类研究

黄作维 胡光伟 谢世雄

(湖南工业大学农牧业废弃物资源化综合利用湖南省重点实验室， 株洲 412000)

0 引言

现代城市是人工环境与自然环境的综合体，地物种类复杂多变，下垫面的组成成分复杂多样。传统的遥感数据分析方法对城市地物光谱特征和区分城市地物所需的光谱信息缺乏，难以反映复杂多样的城市地物覆盖类型,不利于精细分类[1]。高光谱遥感技术弥补了这一缺陷，极大地提升了对复杂地物的识别能力,从而能够区分光谱特征非常相似的城市地物。高光谱遥感能获取地面目标丰富的空间、辐射和光谱3方面的信息，提供了各类地物完整的连续光谱信息，光谱分辨率优于10 nm，使得准确、精细地解译地物成为可能[2]。由于高光谱图像探测距离远，大气传输过程中存在混合效应，成像空间范围场景复杂，空间分辨率低以及遥感仪器本身的混合效应这些因素导致混合像元广泛存在[3-4]。混合像元分解问题一直以来是当代遥感应用中的一个技术难点。高光谱混合像元分解模型一般分为两类，即线性光谱混合模型和非线性光谱混合模型，线性模型是迄今为止使用最多的一种模型，其突出优点是算法简单，物理含义明确。高光谱端元提取是分析和处理高光谱数据的前提条件，BOARDMAN[5]提出了利用凸面几何学分析的纯像元指数提取端元算法，文献[6-7]提出了N-FINDR算法，NASCIMENTO等[8]提出了顶点成分分析算法，文献[9-10]提出正交子空间投影方法。目前地物分类中利用的大部分是光谱信息，而对其空间领域信息利用较少，充分利用光谱信息和空域信息的结合来提高分类精度是一个亟待解决的问题[11-13]。本文利用相邻像素在空间上的相关性，提出一种优化的候选端元判断方法，运用改进的算法，在保证算法效率的同时可有效提高光谱解混的精度。

1 改进的混合像元分解方法

遥感传感器在瞬时视场IFOV内搜集目标，对地物的探测是利用光子探测器或热探测器检测地物所对应的地表物质光谱信号的综合，线性光谱解混的分解流程如图1所示,主要有：数据预处理、端元提取、精度评价、丰度估计。线性光谱解混首先进行端元提取，然后进行混合像元分解[14]。端元是组成混合像元的最基本单位,物理意义是代表某种具有相对固定光谱特征的地物类型，端元提取是混合像元分解的难点[15-17]。通过对高光谱遥感数据分析表明：成像外界因素(噪声)会引起地物光谱特征发生变化，且大多服从多元正态分布，关键图像中的像素是通过光谱和空间特征来共同描述的,高光谱遥感图像在空间分布上具有一定的连续性，在光谱特征上混合光谱(含有噪声)与端元光谱具有相关性(在一定的近邻范围内端元与混合光谱具有相似性)[18-19]。实验表明从物理和数学两个方面的像元空间信息相关参数进行候选端元判断更具有可靠性。

图1 混合像元分解Fig.1 Processing of mixed pixel decomposition

1.1 端元优化方法

端元在物理意义上必须是最或然端元，一定位置的端元与其邻域内的端元属于相同地物类别的概率比较大，而噪声则与其相邻像元光谱不具相似性。在数学上端元矢量矩阵必须是一个非病态矩阵[20-21]。可以用光谱角来描述像元光谱矢量间的相似性。对混合像元解混中存在的“噪声”，利用目标点邻域空间信息的进一步处理可以对一些误分结果起到修正的作用，有效提高分类精度[22-23]。

(1)空间领域方法：假设用M表示高光谱端元矩阵，图像在空间坐标(i,j)处地物位置用M(i,j)来表示。假设以像元为中心，邻域窗口大小为9×9，光谱波段编号为i，光谱波段数为n，对结果矩阵向四周扩张，它的邻域像元表示为

(1)

(2)

则以(i，j)为中心的像元点邻域信息表达为

得到邻域像元的端元集合为

Gi,j={Mi-1,j-1,Mi-1,j,Mi-1,j+1,Mi,j-1,

Mi,j+1,Mi+1,j-1,Mi+1,j,Mi+1,j+1}

(3)

高光谱像元相似性评价：一般可以用光谱角余弦、图像处理的窗口尺寸以及相似光谱百分数等参数来衡量，在地学统计分析方法中，可以用变差函数来描述某一波段光谱亮度的空间自相关性，变差函数为

(4)

式中h——像元之间的欧氏距离

N(h)——距离h之间的像元对个数

位置为xi的像元在j波段的光谱角为dnj(·)，位置为xi+h的像元在k波段的光谱角为dnk(·)。

(2)为充分利用混合光谱的近邻特征，区分端元光谱和噪声的敏感度，用K1表示像元相似性阈值，根据M中每个窗口的像元与中心像元的光谱角余弦进行排序，然后对每个窗口中光谱角余弦取中值，最后对窗口所有中值求平均值，作为像元相似性阈值

(5)

式中k——窗口中光谱编号

C——平均值D——中值

S——排序ri——光谱角余弦

优化准则: 根据高光谱像元分解后误差最大(空间差异性最大)的像元如果是非噪声像元，则其是候选端元。在端元集合Gi，j中，假设候选端元与M中各个端元的夹角余弦最小值为Kmin，利用光谱角阈值判断最大误差的像元矢量是否为候选端元方法是：根据阈值K1,若K1

1.2 改进算法

混合像元分解过程中不可避免地涉及大量矩阵运算和迭代过程，采用单机串行方法处理消耗大量时间，限制了像元分解的效率。如何快速地进行海量数据运算是目前遥感技术发展中亟待解决的问题。随着网络传输性能的提高和图形处理器(Graphic processing unit，GPU)的快速发展，针对高光谱数据计算密集，数据处理过程复杂的特点，GPU能够充分利用NVIDIA的图形处理器的并行计算引擎，由于GPU有强大的浮点计算能力，采用粗粒度的线程块间的并行计算和细粒度的线程间的双层并行计算，随着CPU技术的快速发展，目前对称多处理器(Symmetrical multi-processing，SMP)已经成为最流行的并行计算系统之一。并行计算中的GPU 函数称为内核函数，它可以分为两个方面的并行，就是每个Grid 里有多个块(block)并行，每个block又有多个线程(thread)并行。

在SMP集群的并行结构中，具有分布式内存和共享内存两种并行计算的优势，基于SMP集群下的并行优势，提出了一种改进的端元提取算法。在SMP处理器共享内存模型中，为了提高ICA计算过程中每个并行处理单元的执行效率，算法中把每个处理单元称为线程单元(Thread element，TE)。在实际集群运算中，并行能力能够延伸到各计算节点之间(计算节点产生进程进行节点间的并行计算)，这种实体被称为处理单元(Process element，PE)，为保证运算速度，采用多个TE进行并行计算。

改进的端元提取算法如下：

(1)从高光谱图像中随机抽取像元组成初始端元矩阵，利用PCA变换对图像进行去相关处理，得到新图像Y，并进行数据中心化与白化处理，处理的变换矩阵如下

(6)

其中D=diag(d1,d2,…,dn)E=(e1,e2,…,en)

式中D——协方差矩阵特征值作为对角元素的对角矩阵

E——协方差矩阵单位特征向量为列的矩阵

(2)对矩阵W进行初始化分离，假设一共有a个TE，基于划分原则把图像划分为{Y1,Y2, …,Ya}，第j个TE负责对Yj(j=1,2,…,a)的划分，第j个TE计算公式为

W(j)(k+1)=[E{Yjg(WTX)}-E{g′(WTYj)}W(k)]

(7)

(3)利用式(7)对每一个TE进行计算，当所有TE完成运算后，对W(k+1)进行规约为

(8)

(4)根据以上优化的候选端元判断方法，在端元集合Gi，j中选取候选端元，完成后得到第j个TE(j=1,2，…,a)的丰度估计值为

(9)

(5)对SMP集群的并行结构中的同构的计算节点(PE)进行统计，对主PE图像进行粗划分，假设一共有b个PE，每个PE含有a个TE，对每个PE进行图像划分Y={Y1，Y2，…，Yb}，然后执行步骤(2)、(3)，把第p(p=1, 2，…,b)个从PE细划分为Yp={Yp,1,Yp, 2, …,Yp, a}，运算后得到W(k+1)，并把结果传输给主PE。

(6) 对主PE进行规约得到W(k+1)，进行归一化和对称正交化后传送到各从PE，迭代完成后，从PE获得分离矩阵，根据各从PE丰度最大绝对值进行最大值规约，同时得到所对应的端元ei, j，其中i=1,2，…,m;j=1,2，…,a。

(7)当各TE搜索完候选端元后，根据各端元的丰度最大绝对值进行规约，主PE同步等待各从PE得到的端元提取结果，完毕后获得最终端元集合，算法结束。

2 实验

2.1 实验数据与处理

研究区位于四川省成都市郊区，高光谱数据采集时间是2014年5月，传感器采用上海技术物理研究所研制的PHI。在数据获取期间，对实验区进行了详细的地面调查，选择了具有代表性的1块数据进行实验，其中主要地物类型为道路、水体、建筑物、植被等，实验中选择了0.38～0.96 μm的112个波段，该数据的第3～6、23～34、66、77～89、99～108波段由于为水气体吸收波段和信噪比低而被剔除，其余的72个波段被用于进一步处理，图像大小为223 MB。经过图像配准校正和辐射校正后，选取第7、45、90个波段分别作为R、G、B进行数据合成，合成后的假彩色图像如图2所示。

图2 研究区影像图Fig.2 Image of study area

为了定量评价算法的性能，在分解之前要获得参考地物真实分布情况的各端元的参考光谱，首先从高光谱图像上用 ROI 方法挑选出一片纯净地物(代表一种地物的均匀的区域)，并求出所有像元的均值，然后以此均值为标准建立起来光谱库，该实验的参考光谱是根据图像上的 ROI 区域按照上述方法选取的。在生成过程中一共有92个较纯的像元被选取出来，其中对应植被、水体、道路、建筑的像元数分别为20、12、28、32，取其均值作为该端元的参考光谱。

为了进一步定量评价算法，采用经典算法N-FINDR、ICA 和 OSP等方法在同条件下进行分析。实验运行环境：计算机操作系统为 Windows 7(R)操作系统，CPU型号是intel core(TM)i7，内存4 GB，主频为2.80 GHz，GPU 程序开发环境为 Microsoft Visual Studio.net 2008，并行处理的编程软件采MPICH2以及Microsoft Visual C++2008 OpenMP并行程序设计，能有效降低并行计算的难度和复杂度。

2.2 端元提取实验

根据优化端元提取方法，利用优化后的改进算法(I-ME)，参数的选取与仿真数据相同，不同地物种类植被、水体、道路、建筑物分别用数字编号 1、2、3、4 表示。图3给出了在不同波段下4种端元的光谱曲线比较结果。从图3可以看出，端元各不相同, 另外还采用和参考光谱相比较的方法，提取结果和参考光谱之间的光谱角余弦都较大，各提取端元与参考光谱的比较如图4所示。

图3 端元提取结果Fig.3 Result of endmember extraction

为了定量评价所提取的端元的准确度，本实验利用光谱角距离(SAD)和光谱散度(SID)参数来衡量，表1给出了所提取端元的精度，结果表明，相比其他提取方法该算法(I-ME)有更好的精度，算法的精度由高到低依次为：I-ME、ICA、OSP、N-FINDR。

针对高光谱数据量大的特点，系统运用GPU高性能并行改进算法，随着计算节点的增加，各算法的运行时间逐步减少。图5为几种不同方法的加速比曲线。由图5可知，采用并行算法能获得更高的效率，改进算法的加速比高于其他的方法(ICA、N-FINDR)，比N-FINDR提高了8.5%，充分说明了改进算法的有效性。

2.3 地物分类实验及分析

参考数据是指真实地物图，为考察准确性而选择的像素，本次精度的检验数据是利用了野外调绘的真实数据，共包含4个真实地物类别，标定样本13 219个，随机选择标定样本的5%～20%作为训练样本，剩余作为测试样本，利用地面真实数据与分类结果图进行比较处理。当遥感图像有较多的混合像元时，这种分类方法所得结果的精度(特别是数量精度)会下降。混淆矩阵中数据来源有：已知的典型区域的地物类型图，类前选择的训练区和训练样本时确定的各个类别及其空间分布图和实地调查的分类结果图。本研究采用基I-ME算法进行混合像元分解后再分类，利用 ENVI 软件的Confusion Matrix 功能对上述的分类结果图进行精度评价，精度结果是可靠的。

图4 提取端元与参考光谱的比较Fig.4 Comparison result of extracted endmember and reference spectra

E1(SAD/SID)E2(SAD/SID)E3(SAD/SID)E4(SAD/SID)N-FINDER0.0921/0.03310.1364/0.05360.0862/0.00690.2313/0.0479ICA0.0903/0.03260.1371/0.04980.0812/0.00640.2325/0.0475OSP0.0910/0.03340.1358/0.05080.0834/0.00670.2337/0.0438I-ME0.0876/0.03010.1204/0.04310.0778/0.00680.2066/0.0416

图5 并行算法的加速比Fig.5 Speedup ratio of parallel algorithm

针对传统的地物分类方法(MLC、SVM)，一个重要环节就是特征提取和选择，即为分类寻找最优特征，实现最大可分性。为了表明改进的混合像元解混方法对地物分类精度的影响，将原始图像进行降维采样后，进行对比实验分析。第1组实验是高光谱图像进行预处理后直接利用传统的方法(SVM、MLC)进行分类；第2组实验是先将原始图像进行预处理后利用本文的改进方法进行混合像元分解，再用SVM进行分类(I-ME)；第3组实验是将原始图像进行预处理后利用N-FINDR方法进行混合像元分解后，再利用分类器SVM进行分类。通过比较评价指标来间接评价解混的效果,采用了混淆矩阵、Kappa 系数、总体分类精度、识别精度、虚警率等指标对结果进行评价。

表2～4分别给出了基于SVM分类结果的混淆矩阵、基于N-FINDER算法的混合像元分解后SVM分类的混淆矩阵、基于改进算法的混合像元分解后SVM分类的混淆矩阵。表5～7给出了经过SVM分类后各类地物分类精度、基于N-FINDR混合像元分解后SVM分类的各类地物分类精度、基于改进算法混合像元分解后SVM的分类精度。可以看出，经过解混后(图像空间信息特征的引入)可以改善分类精度。表5给出了用传统方法(SVM)的分类结果，与解混后的表6、7对比，从Kappa系数、虚警率和识别精度来看，解混后的精度有一定提高。

图6 不同分类方法的对比Fig.6 Comparison result of different classification methods

建筑物水体植被道路总和建筑物3653405333731水体 1102382831622737植被 58815638355811道路 1168130668937总和 38322571585690113216

表3 基于N-FINDR光谱解混后SVM分类的混淆矩阵Tab.3 Confusion matrix after spectral unmixing based on N-FINDR

表4 基于改进算法光谱解混后SVM分类的混淆矩阵Tab.4 Confusion matrix based on improved spectral unmixing SVM classification

表5 SVM分类后各类地物分类精度Tab.5 Classification accuracy based on SVM classification %

基于改进算法分解后SVM分类结果与基于N-FINDER算法的混合像元分解后相比，除了道路的Kappa系数从79.34%降到了78.53%，基于改进算法比基于N-FINDER算法的分类指标有较大改善，特别是虚警率指标有明显降低。

表6 解混后各类地物分类精度(N-FINDR)Tab.6 Classification accuracy after spectral unmixing (N-FINDR) %

表7 解混后各类地物分类精度 (改进算法)Tab.7 Classification accuracy after spectral unmixing (improved method) %

实验定量分析如图6所示，图6a为经典SVM分类的结果图，图6b为N-FINDR混合像元分解后SVM分类的结果图，图6c为改进算法混合像元分解后SVM分类的结果。可以看出改进算法能有效减少端元“噪声”像元的数目，与基于N-FINDER算法的混合像元分解相比，改善了地物分类效果，可以看出绝大部分的地物(如建筑物、植被、道路、水体)都能被正确地分类，与实际情况大体相符，证明该分类方法的精度较高。但是由于部分地物之间复杂的空间拓扑关系以及其他不可预见的因素也导致了一些误分的情况，如在植被和裸土相连的部分，裸土上的植被非常少，各种特征不明显而被错误的分类为裸土，同时一部分裸土上有少量的植物被错分为植被，又如当两个中间有植被的建筑物相距很近，而且建筑物和植被高度非常接近时，两者很难分类，但总体的分类精度还是比较理想的。

3 结束语

高光谱图像分类是高光谱遥感应用领域的基础，也是遥感图像处理的核心之一。光谱解混是实现高光谱图像高精度分类的前提，准确地分解混合像元是高光谱遥感技术得以广泛应用的关键性问题。本文通过优化的候选端元判断方法提出了一种改进的混合像元分解方法，相关实验表明解混后得到的分类精度得到了明显改善，证明了改进算法对混合像元分解的有效性。今后还要进一步改进对各种特征因子的空间分析性能，把CPU+GPU异构系统并行优化方法应用于遥感图像处理应用中,提取更优的端元来对数据进行描述，从而获得更高的分类精度。