我要投稿
投稿指南
RSS订阅中分辨率成像光谱仪(Moderate Resolution Imaging Spectroradiometer, MODIS)是EOS系统中主要传感器之一。它具有较高的空间、光谱和时间分辨率[1]，应用前景广阔。但MODIS影像的最高空间分辨率为250m，混合像元广泛存在。混合像元是由于遥感器的成像原理所导致的遥感影像上的一种必然现象[2]。混合像元的存在是传统像元级遥感分类和面积测量精度难以达到使用要求的主要原因[3]。针对混合像元问题，学者们提出了光谱分解模型，通过选取恰当的地表覆盖类型组合，就能在一定程度上提高最终结果的精度[4]。混合像元分解包括两方面主要内容：(1)获取端元类型及其光谱特征值；(2)获取各像元内包含的端元类型及其丰度，即端元在像元内所占的面积百分比。选取合适的端元是成功进行混合像元分解的关键[5][6]。
端元选取包括确定端元数量以及端元光谱值。现有从遥感影像中获取端元的方法中，代表性的算法有像元纯度指数PPI算法[7]，N-FINDR算法[8][9][10]，序列最大角度凸锥SMACC[11]，顶点成分分析VCA[12][[13]，单体增长算法SGA[14]，正交子空间投影算法OSP[15]，序列投影算法SPA[16]，极小体积变换MVT[17]，凸锥分析算法CCA[18]，迭代误差分析IEA[19][20]等。其中PPI算法是最成功的方法之一。它通过将各像元投影到随机测试向量上，获取投影线段端点的像元作为端元。国内外研究者对利用PPI方法进行端元提取做了广泛研究，如Uttam Kumar et al[21]，Singh,et al[22]，Kaiwen Zhong[23]，郑有飞等[24]，樊风雷[25]，裴欢等[26]，李霞等[27]。PPI算法是在主成分变换(PCA)或最小噪声变换(MNF)降维之后进行的，Boardman等[7]指出，降维可能会“忽视”一些小目标，有时候这些小目标正是我们需要的。多光谱MODIS数据可用500m分辨率波段仅有7个波段，降维造成可用波段信息减少，同时使算法产生偏差，不降维情况下噪声增加会使端元提取产生误差。另一方面，陈水森[28]指出，PPI算法在提取纯像元时对端元选择存在不确定性，可能由此导致所提取地物端元纯度降低，或把同类稍有光谱偏差地物，即相同地物类型的不同亚类，当作不同端元对待，使提取的真实地物端元数减少。这种不确定性将加大图像基于端元的后续应用处理的误差。
为了充分利用MODIS七个波段的影像特征，同时避免噪声对PPI算法结果的影响，有效提取给定数目的端元光谱，本文提出一种基于MODIS影像的改进的PPI端元提取算法，在算法处理前不对影像进行降维处理，而是对提取到的备选端元利用光谱角的概念去除噪声像元，同时利用光谱角对同种地物类型的亚类进行标记，避免对同类端元的重复提取。论文最终将改进算法用于MODIS影像的端元提取中，并对利用提取到的端元进行混合像元分解的地物丰度图的精度进行了评价。
像元纯度指数(Pixel Purity Index, PPI)[7]算法认为图像数据集在N维空间中形成一个凸集，像素向量是N维空间中的矢量。首先通过最小噪声变换（MNF）对数据进行降维，接着生成大量穿过数据集合内部的随机测试向量，然后将光谱点分别往各个测试向量上投影，投影结果用一个确定门限选出在这个方向上的极值点。随着向量方向的不断变化，记录下图像中每个像素作为极值点的次数，最后认定出现频率最高的点就是要找的纯点。如图1所示的一条测试向量上，像元e3和e2的投影为极大和极小值，这两个像元将作为备选端元被考虑。
PPI算法有如下缺点[7]：
①PPI算法不是迭代算法，由于每次的投影向量都是随机生成的，获取端元的准确性无法保证，有可能出现同样的图像多次运行PPI算法而得到的端元却不同，或将同类地物像元当作不同端元提取出来的情况；
②PPI算法对图像噪声很敏感，因此必须要用MNF变换去除噪声后才能使用，但是降维会损失信息量；
③使用MNF变换对数据进行降维时，降维后的维度和波段的挑选没有合适的规则；
④端元数量和最大循环次数完全由用户自由指定，没有合适的规则；
针对上述缺点中的①，②和③，本文提出一种改进的PPI端元提取算法。该算法不对影像进行降维处理，利用光谱角来避免相同地物类型像元的选入，同时排除噪声像元的影响。
光谱角[29](Spectral Angle, SA)是指具有同样波长范围的两个光谱向量在光谱空间上所形成的夹角。它的基本原理是：把光谱作为矢量投影到N维空间上，其中N是选取的波段数。在N维空间中，各光谱曲线被看作有方向有长度的矢量，各光谱之间形成的夹角叫做光谱角。图2为二维空间的光谱矢量和光谱角的描述。
在N维空间，光谱角的数学表达为：假设
和
是两个光谱向量信号，SA通过比较光谱向量
和
的角度来测量光谱特性。即：
光谱角的变化范围是
，同类像元之间的光谱角很小，接近于0，因此其余弦值较大，接近于1，而不同类像元之间的光谱角则比较大，相应地其余弦值接近于0。光谱角考虑的是光谱矢量的方向而不是光谱矢量的长度，在计算影像光谱角时，仅考虑光谱矢量的方向而不用关注影像本身的亮度，所以同类稍有光谱偏差的地物会被赋予较高的相似性。光谱角在衡量像元光谱相似性方面具有明显的优越性，以其为度量基准的光谱角制图模型(SAM)在高光谱遥感影像分类、聚类、典型信息提取等方面都得到了广泛应用。
原始PPI算法将影像数据降维后提取端元，导致不能充分利用影像的光谱信息；另一方面PPI算法会把同类稍有光谱差异的像元作为不同端元来提取。为了充分利用MODIS数据7个波段的光谱信息，同时避免同类端元的重复提取，本文提出一种改进的PPI算法——Imp_PPI。针对原始PPI算法存在的问题，Imp_PPI算法所做的改进有两点，分别为：
改进一：对原始PPI算法提取到的N个端元（本文称作备选端元），两两计算光谱角，并将光谱角与给定的阈值
比较大小，对于光谱角小于
的两个备选端元，认为它们属于同种地物，并将其中的一个顺
次用后面较大PPI的像元替代。具体算法步骤如图3中红框所示。这一改进针对上述PPI算法缺点中的①，光谱角能够有效的衡量像元光谱的相似性，剔除备选端元中同种地物类型的像元。
改进二：由于本文的算法在提取端元前没有进行去噪处理，因而端元易受噪声的影响，为了避免结果端元中含有噪声光谱，本文的第二个改进为：对每个剔除了同类像元的备选端元，以它为中心开一个W*W的窗口，计算窗口中各邻域像元与中心像元的光谱角，并与给定阈值
比较大小，统计小于
的光谱角个数，并将其与给定阈值X相比较，若小于
的光谱角个数大于X，则说明该中心像元与其邻域中的大多数像元光谱相似，根据遥感图像上光谱值得连续性，可判断该中心像元不是噪声，并将它作为最终端元输出。这一改进针对PPI算法缺点中的②和③而提出。尤其对于MODIS数据，500m分辨率的7个波段经过MNF变换后，会造成严重的光谱信息损失，影响到各地物在影像上的光谱特征，使计算纯像元指数PPI的结果受到影响，最终导致提取到端元光谱的偏差。Imp_PPI算法的具体流程图如图3所示。
5 实验与结果分析
根据原始PPI和本文提出的Imp_PPI算法，作者在matlab R2010a下，用MODIS数据试验了两种算法提取端元光谱和混合像元分解的效果；即用两种算法提取了给定数目的端元光谱，并实现了混合像元分解。具体分如下两步进行，首先用PPI和Imp_PPI算法提取出端元光谱，然后再基于限制性模型求解各像元所占的比例。
本文选定的研究区位于东经84°50'~86°09'，北纬43°07'~43°50'，海拔1400~5218 m之间的天山中段。实验所用的MODIS L1B数据获取于2007年5月15日。该时段影像上主要的地物类型有积雪、植被、岩石等。数据经过去条带、去bow-tie、几何校正和辐射定标等预处理之后，得到覆盖研究区的反射率影像。为了评价两种端元提取算法对于混合像元分解结果的精度，本文用30m分辨率的TM影像作为参考“真值”，对分解结果进行定量评价。获取的同一时期的Landsat 5 TM数据来源于美国地质调查局USGS（Geological Survey of the United States），经辐射定标、拼接和裁剪，生成与MODIS数据同一坐标系下的研究区第1～7波段反射率影像。图4(a)是TM影像7,5,1波段的RGB合成图，图上对主要的地物类型进行了标示：图中S表示积雪；V代表植被；R是土壤；M代表阴影。
图3 Imp_PPI算法流程图
(a) (b)
(c) (d)
图4 研究区地物覆盖图
5.1 端元提取及结果分析
通过图4(a)的目视解译结果，指定影像上提取的端元数目N=4，分别为积雪、植被、岩石和阴影，并且通过多次试验以及对MODIS影像光谱值得分析，指定Imp_PPI算法中各阈值分别为：迭代次数iter_n=5000；光谱角阈值
=10°；光谱角阈值
=20°；X=10；窗口W=5；将以上各值代入PPI和Imp_PPI算法中，得到4个端元的光谱值，并根据地物光谱曲线，确定每个端元所代表的地物类别，得到图5(a)(b)所示的端元光谱曲线。
(a)端元数目N=4时，PPI算法提取端元光谱
(b)端元数目N=4时，Imp_PPI算法提取端元光谱
图5 Imp_PPI和PPI算法提取的四个端元光谱曲线
由图5(a)可知，PPI算法提取到的4个端元中，每两个属于同种类型，它们的光谱值在每个波段稍有差异，PPI算法却将其作为两种地物端元提取出来，而Imp_PPI算法能够有效识别相似端元，剔除同种地物类型的端元，从而提取到给定数目的端元光谱，如图5(b)所示。
为了利用PPI算法提取到MODIS影像上各地物的光谱，本文将算法的端元数目N设置为8，重新计算，算法最终输出表1所示的8个像元光谱。
表1 PPI算法提取的8个端元光谱值
|
反射率 |
端元1 |
端元2 |
端元3 |
端元4 |
端元5 |
端元6 |
端元7 |
端元8 |
|
|
0.0343 |
0.0343 |
0.4463 |
0.8327 |
0.8327 |
0.0348 |
0.6657 |
0.7712 |
|
Band2 |
0.3295 |
0.3295 |
0.4097 |
0.7781 |
0.7823 |
0.3210 |
0.6232 |
0.7152 |
|
Band3 |
0.0226 |
0.0254 |
0.4905 |
0.8816 |
0.8843 |
0.0326 |
0.7114 |
0.7704 |
|
Band4 |
0.0508 |
0.0535 |
0.4594 |
0.8407 |
0.8407 |
0.0544 |
0.6688 |
0.7648 |
|
Band5 |
0.3146 |
0.3146 |
0.1908 |
0.3926 |
0.4014 |
0.3493 |
0.3489 |
0.2695 |
|
Band6 |
0.2053 |
0.2053 |
0.0434 |
0.1241 |
0.1241 |
0.2051 |
0.0919 |
0.0593 |
|
Band7 |
0.0814 |
0.0814 |
0.0393 |
0.0393 |
0.0393 |
0.0813 |
0.0394 |
0.0392 |
结合地物光谱曲线分析表1可知，在PPI算法提取出的8个端元中，端元1，2和6是植被，端元4和5是积雪，端元3，7和8是积雪的亚类，仍没有提取出土壤和阴影端元。继续增加端元数目N=20，土壤能被提取出来，同时积雪、植被、土壤三种地物的亚类被提取出来，但是仍不能提取出阴影端元。
综合分析图5，表1和上述利用两种方法提取端元的过程可知，传统PPI算法在提取端元光谱时，采用的MNF变换在去除噪声的同时，造成的信息损失会导致某些端元无法被有效的识别，如本文实验中的阴影端元。由于阴影在MODIS影像各波段上均表现为较暗的特征，在MNF变换后的影像上不占据主要特征，所以无法通过计算纯净像元指数获取到。另一方面，传统PPI方法由于缺少对同种类型不同亚类像元的剔除过程，会导致提取出的给定数目的端元中包含的地物类型数目减少而亚类数目增多，如表1中的端元1，2和6，端元3，4，5，7和8。改进后的算法Imp_PPI通过省略传统算法中前期的MNF变换来保持影像的光谱特征，从而充分利用不同端元在影像波段上的特性将其提取出来，通过计算光谱角的方法来避免噪声端元和同种地物不同亚类被选入最终端元，能够有效的提取出满足数量要求且准确的端元光谱.
5.2 混合像元分解及结果分析
利用上述两种方法提取到的端元光谱值，采用限制性混合像元分解模型对MODIS影像进行混合像元分解，得到各地物的丰度图。具体做法如文献[24]所述。根据上述两种算法提取的端元，本文从TM影像中选取训练样本，利用最大似然分类法并结合目视解译得到了与端元对应的分类结果图。由于研究区影像上积雪为主要覆盖物，且不同地物后续精度分析原理相同，所以本文以积雪分类图为代表展开后续精度评估。如图4(b)(c)(d)所示，其中，图4(b)是TM影像的积雪分类图，图4(c)是以PPI算法结果作为端元光谱的混合像元分解积雪丰度图，图4(d)是以Imp_PPI算法结果作为端元光谱的混合像元分解积雪丰度图。
由图4(c)和(d)可知，线性混合模型能有效的分解混合端元，得到地物覆盖的丰度图。为了定量评价两种端元提取方法所得到的地物丰度图的精度，本研究采用数量精度作为指标来进行混合像元分类结果的精度评价，即通过数量比较来确定像元的分类精度。
对随机样本来说，其某一地类的数量精度
的计算公式为[30]：
(2)
其中，
为MODIS影像上第i个样方第j个类型的比例；
为TM影像上第i个样方第j个类型的比例。由于像元分解结果得到的像元精度值是连续的，依据连续尺度的样本大小计算公式，完成影像整体精度评价所需要的样本数n为：
(3)
式中，S为试探性样本的标准差；
为试探性样本精度值的均值；
为试探性样本的变异系数；Z为对应置信水平，从正态分布Z值的概率表上所查出的值，D为误差允许范围。如果抽样比n/N<0.05（N为研究区像元总数），就取n，否则对n进行修正，取修正值
。
利用抽样的方法，直接在MODIS地物丰度图上进行随机抽样。在对抽样像元进行检验时，为了避免或减少配准误差所带来的计算误差，本研究在MODIS积雪丰度图中抽取10个6×6的像元集作为试探性样本，与抽取的同位置的TM积雪分类图43×43像元集进行分类精度估算，得到各丰度图的样方精度值
，并统计样方精度值的均值
和S方差，如表2所示。
表2 以图5作为端元光谱的混合像元分解积雪丰度图的10个样方数量精度及其均值和标准差
|
|
||||||||||||
|
|
||||||||||||
通过表2的样方精度可知，在分解方法相同的情况下，Imp_PPI端元提取方法得到的积雪丰度图上的样方精度均值高于PPI方法，十个样方的平均精度提高了约2.7%。
当取置信水平为95%时，Z=1.96。取误差允许范围d±2.5%，依(3)式计算可得，本研究中完成各端元提取方法分类结果精度评价所需的样本总数分别为：
，即nPPI =102，nImp_PPI=96。
对MODIS积雪丰度图而言，研究区总像元数为N=42823，则n/N=0.0024和0.0022，均满足n/N<0.05的要求。在积雪丰度图上添加102和96个样方数的工作量较大，本文采用简略方法，再分别选取90个样本点，样本点分布尽量均匀。添加样本点后，积雪丰度图上平均样方精度值和样方精度值间的方差趋于稳定。计算可得，两种端元提取方法的分解精度分别为：QPPI=77.08%，QImp_PPI=79.27%。
通过以上对分解结果精度的评估可知，利用改进后的PPI算法Imp_PPI提取到的端元光谱进行混合像元分解，得到的地物丰度图的精度高于传统PPI算法的分解结果，整体精度提高了2.19%。这表明Imp_PPI提取的端元的个数和光谱值与影像的地物类型和光谱特征的符合度高于原始PPI算法，从而在利用这些端元光谱进行混合像元分解时，能更加准确的实现像元光谱的解混。另一方面，通过对不同位置上样方精度的分析可知，在大片地物内部，两种方法精度都较高，Imp_PPI分解精度改善不明显，如样方2和样方8，这是因为在影像中大片地物内部区域混合像元较少，主要是纯净像元，在求解线性光谱模型时，对这些像元基本不做分解处理，因为端元光谱对其影响较少。而在不同地物的边界处，两种方法中精度总体低于地物内部，且Imp_PPI方法分解精度有明显提高，如样方1,3,4,5,6,7,9和10。主要原因是混合像元主要出现在不同地物的边界处，端元光谱以及类型的差异都会影响到线性分解模型的分解结果，原始PPI算法不能有效的识别阴影端元，模型对各混合像元的分解没有考虑阴影的混合效果，因而会影像到对其他地物类型的分解精度。
总体而言，改进后的PPI方法能够有效利用MODIS影像的波段信息，高效提取端元光谱；利用改进后的PPI方法提取到的端元光谱进行混合像元分解，得到的地物丰度图的精度高于以原始PPI方法提取的端元光谱进行分解的结果，能够提供更准确的地表覆盖类型。
6 结论
本文针对MODIS数据的波段特征，对传统PPI算法进行了改进，弥补了传统PPI算法需要对数据进行降维处理带来影像信息损失，以及对同类端元的亚类重复提取的缺点，能够基于原始影像信息，提取端元光谱，并且有效的避免噪声像元的影响，也能避免将同类地物像元当作不同端元提取出来。本文利用改进后的PPI端元提取方法，结合线性光谱分解模型，对MODIS影像进行了线性光谱解混，利用数量精度进行了精度估计，实验的结果令人满意。
参考文献
[1] 曹云刚,刘闯.从AVHRR到MODIS的雪盖制图研究进展[J].地理与地理信息科学,2005,21.
[2] 胡茂桂,王劲峰.遥感影像混合像元分解及超分辨率重建研究进展[J].地理科学进展,2010,29(6):747-759.
[3] 赵英时.遥感应用分析原理与方法[M].北京:科学出版社.2003:328.
[4] 梁天刚,吴彩霞,等.北疆牧区积雪图像分类与雪深反演模型的研究[J].冰川冻土,2004,26(2):160-165.
[5] Tompkins S., Mustard J.F., et al.. Optimization of endmembers for spectral mixture analysis[J]. Remote Sensing of Environment, 1997.59(3): 472-489.
[6] Elmore A.J., Mustard J.F., et al.. Quantifying vegetation change in semiarid environments: Precision and accuracy of spectral mixture analysis and the Normalized Difference Vegetation Index[J]. Remote Sensing of Environment, 2000. 73(1): 87-102.
[7] Boardman J W, Kruse F A et al.. Mapping target signatures via partial unmixing of AVIRIS data: in Summaries[C]. Fifth JPL Airborne Earth Science Workshop. Pasadena: JPL Publication:23–26,1995.
[8] Winter M E. N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data[C]. Imaging Spectrometry V, SPIE Proceedings. Denver CO, USA, 3753: 266–275,1999a.
[9] Winter M E. Fast autonomous spectral endmember determination In hyperspectral data[C]. Proceedings of the Thirteenth International Conference on Applied Geologic Remote Sensing. Vancouver, B C, Canada: 337–344,1999b.
[10] Winter M E. A proof of the N-FINDR algorithm for the automated detection of end-members in a hysperspectral image Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery X[C], SPIE Proceedings. Orlando, USA: 12–15. 2004.
[11] Gruninger J H, Ratkowski A J, et al.. The sequential maximum angle convex cone (SMACC) endmember model. Algorithms for Multispectral, Hyperspectral and Ultraspectral Imagery X[C], SPIE Proceedings. Orlando, USA, 5425: 1–14. 2004.
[12] Nascimento J M P and Dias J M B. Vertex component analysis: a fast algorithm to unmix hyperspectral data[J]. IEEE Transactions on Geoscience and Remote Sensing, 2005. 43(4):898–910.
[13] Nascimento J M P. Supervised Hyperspectral Unmixing. Lisboa: Universidade Técnico de Lisboa: 2006.47–72.
[14] Chang C I, Wu C C, et la.. A new growing method for simplex-based endmember extraction algorithm[J]. IEEE Transactions on Geoscience and Remote Sensing, 2006.44(10): 2804–2819.
[15] Harsanyi J C, Chang C I. Hyperspectral image classification and dimensionality reduction: an orthogonal subspace projection approach [J]. IEEE Transactions on Geoscience and Remote Sensing, 1994.32(4): 779–785.
[16] Zhang J K, Rivard B, et al.. The successive projection algorithm (SPA), an algorithm with a spatial constraint for the automatic search of endmembers in hyperspectral data. Sensors, 2008,8:1321–1342.
[17] Craig M D.. Minimum-volume transforms for remotely sensed data[J]. IEEE Transactions on Geoscience and Remote Sensing, 1994,32(3): 542–552.
[18] Ifarraguerri A and Chang C I.. Multispectral and hyperspectral image anlysis with convex cones[J]. IEEE Transactions on Geoscience and Remote Sensing, 1999, 37(2): 756–770.
[19] Berman M, Kiiveri H, et al.. ICE: an automated statistical approach to identifying endmembers in hyperspectral images[J]. IEEE International Geoscience and Remote Sensing Symposium, France: Toulouse,2003, 1:279–283.
[20] Berman M, Kiiveri H, et al.. ICE: a statistical approach to identifying endmembers in hyperspectral images. IEEE Transactions on Geoscience and Remote Sensing[J]. 2004, 42(10): 2085–2095.
[21] Kumar, U.. Constrained Linear Spectral Unmixing Technique for Regional Land Cover Mapping Using MODIS Data[C]. International Joint Conferences on Computer, Information, and Systems Sciences, and Engineering, 2007.
[22] Singh, Dharmendta. A Pixel Purity Index and Curyelet based approach for the fusion of ASTER and MODIS data for land cover classficaiton[C].37th COSPAR Scientific Assembly.2008.
[23] Kaiwen Zhong, Xulong Liu, et al. Endmembers extraction using time series of MODIS and TM samples[J], Proc. SPIE 6790, 67900R (2007).
[24] 郑有飞,范旻昊,等.基于MODIS遥感数据的混合像元分解技术研究和应用[J].南京气象学院学报,2008,31(2):145-150.
[25] 樊风雷.基于线性光谱混合模型(LSMM)的两种不同端元值选取方法应用与评价——以广州市为例[J].遥感技术与应用,2008,23(03):272-277.
[26] 裴欢,房世峰,等.基于遥感的新疆北疆积雪盖度及雪深监测[J].自然灾害学报,2008.(05):52-57.
[27] 李霞,王飞,等.基于混合像元分解提取大豆种植面积的应用探讨[J].农业工程学报.2008,24(1):213-217.
[28] 陈水森.基于波谱库的作物纯像元识别与种植面积遥感估算[D].北京:中国科学院遥感应用研究所,2005.
[29] 唐宏,杜培军,等.光谱角制图模型的误差源分析与改进算法[J].光谱学与光谱分析,2005,25(8):1180-1183.
[30]吴健平,杨卫星.遥感数据分类结果的精度分析[J].遥感技术与应用,1995,10(1):17-24.
