常见混合像元分解方法简介二

端元就相当于一个像素里的亚像元，只包含一种地物的光谱信息，根据多光谱或高光谱的高光谱分辨率可以提取出来。

端元只包含一种地物信息，一般的像元都为混合像元，包括多种地物，在进行混合像元分解的时候，可以对一个像元中包括的几种端元进行定量描述，求得每个像元中几种端元在这个像元中的面积百分比，即端元的丰度。

混合像元分解
(2011-06-10 14:46:57)
转载▼
分类：ENVI/IDL学习
标签：
杂谈
混合像元是指在一个像元内存在有不同类型的地物，主要出现在地类的边界处。

混合像元的存在是影响识别分类精度的主要因素之一,特别是对线状地类和细小地物的分类识别影响较为突出，在土地利用遥感动态监测工作中,经常遇到混合像元的难题，解决这一问题的关键在于通过一定方法找出组成混合像元的各种典型地物的比例。

线性混合像元分解
由于线性模型是应用最广泛，也是研究最多的算法，下面重点介绍基于线性模型的混合像元分解算法。

一般而言，混合像元分解算法包括数据降维、端元选取和反演三个步骤。

１.数据降维
尽管数据降维不是混合像元分解算法的一个必需步骤，但由于大多数算法都将其作为一个流程，我们也将其当作一个步骤。

常用的降维算法有主成分分析（Principle Component Analysis，PCA）、最大噪声比变换（Maximum Noise Fraction，MNF）和奇异值分解（Singular Value Decomposition，SVD）。

(1) 主成分分析：遥感图像各波段之间经常是高度相关的，因此所有的波段参加分析是不必要的。

PCA就是一种去除波段之间相关性的变换。

PCA通过对原数据进行线性变换，获得新的一组变量，即主成分。

其中前几个主成分包含了原数据主要方差，同时各个主成分之间是不相关的。

(2) 最大噪声比变换：最大噪声比变换（Maximum Noise Fraction，MNF）[24]由Green等（1989）提出，该变换通过引入噪声协方差矩阵以实现对噪声比率的估计。

首先，通过一定方式（比如对图像进行高通滤波）获取噪声的协方差矩阵，然后将噪声协方差矩阵对角化和标准化，即可获得对图像的变换矩阵，该变换实现了噪声的去相关和标准化，即变换后的图像包含的噪声在各个波段上方差都为1，并且互不相关。

最后对变换后的图像再做主成分变换，从而实现了MNF变换，此时得到的图像的主成分的解释方差量对应于该主成分的信噪比大小。

(3) 奇异值分解：奇异值分解（Singular Value Decomposition, SVD）也是遥感图像处理中常用的变换，与PCA类似，SVD能够找出包含原始数据大部分方差的特征方向，不同的是，SVD特别适合于波段间高度相关的数据，而PCA在这种情况下很有可能会失败[25]。

2.端元选取
选取合适的端元是成功的混合像元分解的关键[26, 27]。

端元选取包括确定端元数量以及端元的光谱。

理论上，只要端元数量m小于等于L+1（L表示波段数），线性方程组就可以求解。

然而实际上由于端元波段间的相关性，选取过多的端元会导致分解结果更大的误差，尽管此时残差会减少[28]。

在能够描述一个场景内光谱的大部分方差的前提下，越少的端元数量是越好的选择[29]。

对于城市地区，最常用的端元选取方式是由Ridd等（1995）[30]
提出的植被-不透水层-土壤端元模型（Vegetation - Impervious surface–Soil, V-I-S），V-I-S
模型在很多研究中得到应用[2,31-35]。

而在非城市地区，一般采用植被-土壤-阴影（或干植被）端元模型。

端元的数量和类型确定后，下一步是确定端元的光谱。

端元光谱的确定有两种方式：(1) 使用光谱仪在地面或实验室测量到的“参考端元”；(2) 在遥感图像上得到的“图像端元”。

参考端元虽然可以精确测量，但由于各种因素（包括不同传感器、大气影响、辐射条件及物候等）造成的噪声，可能会导致其与图像上像元光谱的不一致。

要将二者匹配起来需要进行复杂的校正，而且参考端元的微弱噪声就可能引起最后计算得到的端元比例有很大误差。

相对而言，直接从图像上寻找端元更加直接方便，因而得到广泛研究。

图像端元选取的方法大致可以分为两类：交互式端元提取和自动端元提取。

2.1交互式端元提取
从图像上选取端元的各种方法，都基于这样一种思想：在特征空间中，所有的混合像元都存在于由端元连接而成的多边形（或多面体）内。

这样的混合像元才能满足端元面积比例为正值并且总和为1的条件。

因此，最简单的交互式提取方法就是在特征空间中（通常是前两个或三个主成分构成的特征空间）目视寻找多边形的顶点作为端元。

为了减少目视选取的主观性，一些定量化方法被引进作为选取端元时的参考，这些方法的特点是仍然需要人工参与，所以被称为交互式提取。

(1)PPI指数[36]：像元纯度指数（Pixel Purity Index）是最成功的方法之一。

首先对图像进行MNF变换以实现数据降维，接着在由MNF的前几个主要成分组成的特征空间中，随机生成穿过数据云的测试向量，然后将数据点投影到测试向量上。

投影在测试向量两端的数据点有较大的可能性属于端元，用一个阈值选出在这个测试向量两端的极值点。

继续生成新的随机向量，重复上述步骤，记录图像中每个像元作为极值点的频度，即为PPI指数。

PPI指数越高意味着像元的纯度也越高。

(2)MEST算法[37]：MEST（Manual Endmember Selection Tool）算法通过主成分分析来确定混合物中端元的数目。

对于三维及三维以内的数据而言，可以直接通过目视选取多边形的顶点来确定端元，但对于更高维数的数据，显示上存在困难。

MEST提供了一种在高维空间中寻找端元光谱的方法。

(3)CAR和EAR指数[38]：由于端元内的光谱差异，一类端元往往对应很多条光谱，这两个指数用于解决从多条光谱中选择最具代表性的光谱问题。

CAR（Class Average RMSE）计算A类端元光谱用B类端元光谱来分解产生的残差,显然残差越小，A、B两类端元的混淆性越大。

EAR（Endmember average RMSE）计算A端元内某一条光谱用A端元内其他光谱来分解产生的残差。

显然EAR越低表明这条光谱的代表性越好，如果很高则证明这条光谱可能是离群点，没有代表性。

因此，利用CAR和EAR指数可以为端元的合并或分离提供依据，进而指导端元选取。

2.2自动端元提取
自动端元提取指不需要人工参与的端元提取方法。

相对于交互式提取方法，自动方法的提取
结果不受个人的主观性影响，但是由于自动提取方法一般采用纯数学判据，有时可能产生不具有物理意义的端元。

常见的自动提取方法有：N-FINDR、IEA、CCA和ORASIS。

(1) N-FINDR[39]：N-FINDR寻找这样一组像元，它们在特征空间中构成的单形体具有最大的体积。

首先对图像作MNF变换，随机选择m个像元作为端元，计算这些像元构成的单形体的体积。

然后依次将每个像元光谱代入各个端元位置，计算体积，若体积增大，则将该像元光谱替换为端元。

不断循环直至体积不再增大，此时得到的端元即最终的端元。

(2) IEA[40]：IEA（Iterative Error Analysis）寻找端元使得最小二乘法的拟合残差最小。

IEA 无需对数据进行降维，它首先选择图像均值作为初始向量，根据这个向量开始进行第一次分解，得到一幅残差图像。

选择一定数量的残差最大的像元形成像元集，对该集合中光谱角较小的像元光谱求平均作为第一个端元。

然后，以第一个端元为初始向量开始第二次分解，重复以上步骤得到了第二个端元。

以此类推，最后直到得到了规定数量的端元或残差图像小于某个阈值，停止运算，所得到的端元即为最终的端元。

(3) ORASIS[42]：在ORASIS (Optical Real-time Adaptative Spectral Identification System)系统中，有一个称之为“样本选择”的算法来选择端元。

首先选择一组样本像元，再与图像的其他数据作光谱角匹配，去除光谱角相近的像元，只保留光谱角超过某个阈值的像元作为新的样本光谱。

接着利用修正的Gram–Schmidt正交化过程得到数目较少的一组基来构成子空间，再将样本光谱投影到子空间中，利用最小体积算法找到子空间中包含所有数据点的外接单形体，单形体的角点即为端元。

由于外接单形体的角点并不一定是图像上的数据点，因此即使图像上没有包含纯像元，ORASIS也可能找到端元光谱。

(4) 结合空间信息的方法：上述几种方法都只是利用了像元的光谱信息，近年来一些学者提出了结合空间信息的方法。

Plaza等（2002）[43]提出了一种基于数学形态学的AMEE算法(Automated Morphological Endmember Extraction)。

它首先将建立在二值图像上的腐蚀和膨胀算子拓展到多光谱图像中，Plaza等证明了腐蚀算子能够寻找到结构元里的混合最严重的像元，而膨胀算子则可寻找到其中的最纯像元，并提出MEI指数用于表示结构元中纯像元的纯度。

将结构元依次通过图像中的每个像元，可以在结构元定义的邻域中找出最纯的像元及其MEI指数。

不断增加结构元的大小，重复计算像元的MEI，并将不同大小结构元的MEI 作平均得到最终的MEI图像。

MEI指数大于一定阈值即为纯像元。

另外，Rogge等（2007）[44]提出的SSEE算法（Spatial-Spectral Endmember Extraction）通过结合先验的空间知识来获取更准确的端元光谱。

该算法的主要思想是：在图像的特定区域，混合像元包含特定的几个端元组分，进而可以将特征空间中的数据按照真实空间位置划分成不同的组，再做端元提取，以获取更高的精度。

Plaza等（2004）[45]比较了各种端元选取方法的优劣，结果表明结合空间信息的AMEE比其他单纯利用光谱信息的方法更具优势。

4.3反演
最简单的线性混合模型反演过程是无约束的最小二乘法，该方法假设观测误差w为高斯白噪声，此时估计的比例有解析表达式：F=(X’X)-1X’r。

考虑到端元比例的物理意义，需要满足约束条件Σfi=1，则可通过拉格朗日乘子法求解部分约束最小二乘法（abundance sum-to-one constraint，ASC），也有解析解。

Heinz等(2001)[46] 提出全约束最小二乘（fully constrained least squares，FCLS）的求解方法。

FCLS在原方程组里添加方程Σfi=1，并给予其很大的权重，再用ANC对新的方程组求解，则求得的解近似满足两个约束条件。

由于该方法能够同时满足两个约束条件以及高效的运算效率，因而得到广泛认可。

但是也有学者认为，准确的混合像元模型的分解结果应当自动满足两个约束条件，强加的约束条件可能会导致新的误差[47,48]。

光谱混合模型
在遥感图像中，一个像元往往覆盖几米甚至上千平方米的地表范围，其中可能包含着多种地物类型，这就形成了混合像元。

混合像元的存在主要有以下两个原因[1] ：一是传感器的空间分辨率较低，不同的地物可能存在于一个像元内，这种情况一般发生在遥感平台处于比较高的位置或者拥有宽视角；二是不同的地物组合形成同质均一化的地表类型，这种情况的发生不依赖于传感器的空间分辨率。

混合像元分解技术假设：在一个给定的场景里，地表由少数的几种地物（端元）组成，并且这些地物具有相对稳定的光谱特征，因此，遥感图像的像元反射率可以表示为端元的光谱特征和它们的面积比例（丰度或盖度）的函数。

混合像元分解途径一般通过建立光谱的混合模型实现。

由于假设地物（端元）具有相对稳定的光谱特征, 场景中不同像元间光谱的差异主要是端元比例变化的结果。

光谱混合模型
1. 线性模型
线性混合模型（Linear Mixing Model, LMM）[8,9]是最简单，也是应用最广泛的光谱混合模型。

在线性模型中，混合光谱等于端元光谱与端元面积比例的线性组合。

该模型基于以下假设[10] ：到达遥感传感器的光子与唯一地物发生作用（即不同地物间没有多次散射）。

线性模型的数学表达式如下：r=Σfixi+w。

为使LMM具有物理意义，需要受到两个约束条件限制：一是端元面积比例之和为1，即；二是所有的端元比例都为非负的。

2.非线性模型
(1) 高次多项式模型
高次多项式模型通过考虑端元之间的交叉项来描述光谱混合的非线性效应。

二次多项式的混合模型的数学表达式为：r=Σfixi+Σfijxixj。

Ray等（1996）[12]首次提出该模型，并将其应用于沙漠植被盖度反演，Zhang等(1998)[13]将该模型用于潘阳湖地区的土壤、植被盖度研究。

该模型尽管可以很好地模拟混合光谱中多次散射的贡献，但是所计算得到的端元比例却不能直接对应于地物的盖度，需要将多次散射项的虚拟比例（fij）合理分配到各端元的盖度上。

一般认为，多次散射主要发生在植被内部，故将多次散射项的虚拟比例全部分配给植被端元[14]。

Ben Sorners等(2009)[15]通过地面实验发现将交叉项的虚拟比例按计算得到的一次散射项的比例分配给各个端元能够取得更好的精度。

(2) 几何光学模型
几何光学模型[16]适用于植被覆盖地区，它将植被简化为一系列规则的几何体（如圆锥、椭球体等），进而假设地表由四种组分（端元）构成：植被光照面（C）、植被阴影面（T）、光照背景面（G）、背景阴影面（Z）。

这些组分所占的面积比例与树冠、树高、树密度、太阳入射角、观测角都有关。

因此，在实际应用中，往往要对地表特征进行简化，假设树冠具有相同的规则几何形状，假设树木的高度，位置满足一定的分布。

几何光学模型基于植被景观的三维几何特征，需要关于植被的详细几何参数。

(3) 其他非线性模型
概率模型基于Marsh等人（1980）[17]提出的近似最大似然法。

该模型将端元比例求取问题转化为分类概率问题。

通过各端元光谱与混合光谱的相似性判别分析，产生一个判别值，根据判别值决定端元所占的比例。

该模型只考虑两种地物混合的情况。

另一个有代表性的工作是由Ju等（2003）[18]提出的MDA（Mixture Discriminant Analysis）方法。

MDA方法将混合像元光谱分布用多个高斯分布的端元光谱曲线叠加来表示，其中分布的参数估计由期望最
大算法(Expectation Maximization, EM)算法完成，进而根据贝叶斯原理得到像元归属于各类端元的后验概率。

由于后验概率满足非负且总和为1的两个约束条件，因此可以直接用后验概率来表示混合像元的端元比例。

模糊模型由Wang（1990）[19]提出，其基本原理与概率模型相似，即将端元比例求取问题转化为模糊分类的隶属度问题。

该模型根据模糊数学的运算法则确定各端元光谱的模糊均值和模糊协方差矩阵，进而得到混合像元对各端元的隶属度，隶属度即可对应于相应的端元比例。

神经网络本质上是一种机器学习算法。

与基于神经网络的硬分类方法不同点在于：用于混合像元分解神经网络的输入层为混合像元光谱，输出层为像元归属于各类端元的判别值，即相应的端元比例。

而用于混合像元分解的神经网络的训练样本不再只归属于某一类，而是将训练样本的端元比例作为输出层的值进行学习。

由于神经网络是非线性加权模型，适当的神经网络能够模拟出非线性混合过程。

但是神经网络的表现效果与其自身的结构有很大的关系，常见的用于混合像元分解的神经网络结构有：BP神经网络（back-propagation neural network）[20]、ARTMAP[21]、ART-MMAP[22]等。

评论
线性模型建立在像元内相同地物都有相同的光谱特征以及光谱线性可加性基础上，同时忽略了多次散射过程。

其优点是物理意义明确、构模简单。

但其基本假设光谱线性可加性和对多次散射过程的忽略处理是对现实情况的粗略近似，可能直接影响到像元分解的精度。

非线性模型中的高次多项式模型和几何光学模型考虑了多次散射过程。

高次多项式模型通过考虑端元之间的交叉项来描述光谱混合的非线性效应，几何光学模型则把地面当成三维物体考虑多次散射。

而其他非线性模型则主要将将端元比例求取问题转化为软分类问题。

这些模型虽然简单方便，但缺乏清晰的物理意义。

在不同的情况下使用不同的混合模型是非常重要的。

有研究表明[23]，将线性模型用于非线性光谱混合的数据，反演的端元比例误差会达到30%。

然而，已有研究尚未给出各种混合模型的适用条件及其定量判据，模型的选择只能根据使用者的经验进行判断。

如何根据地物端元特征选择混合模型的研究仍需要大力加强。

同时，尽管在特定的表面使用特定的非线性模型有明显的优点，但非线性模型并没有得到广泛的应用，主要原因是地物的详细散射参数通常很难获取，而这些参数往往会对最后的反演结果有着相当大的影响[1]。

端元内光谱差异问题或叫做端元变异，端元不稳定(Endmember variation)。

一般的混合像元分解算法假设相同地物都有相同的光谱特征，因而对整幅图像采用相同的端元光谱。

但由于同物异谱现象的存在，端元的光谱并非恒定的值，这就是端元内光谱差异现象。

这种现象的存在常常会导致分解结果的误差。

目前，解决该问题的方法可以分为四类：(1) 多端元方法
多端元方法指对每一类地物选取多个端元光谱参与混合像元分解。

其中最典型的方法是由Roberts等(1998)[49]提出的MESMA(Multiple Endmember Spectral Mixture Analysis)方法。

该方法首先为每类地物选取多条光谱，并以此生成多个端元组合（每个端元组合由不同地物中的某一条光谱组成），接着对每个像元寻找最小二乘法误差最小的端元组合，进而求出每个像元的端元比例。

该方法在很多研究中被证实是十分有效的[50-54]。

Bateson等(2000)[55]
提出了一种端元束的方法，该方法对每类地物生成端元束（一个端元束由许多同一类地物的
光谱组成），将所有端元束的光谱作为端元进行混合像元分解。

因为端元数目超过光谱波段数，方程组欠定，所以只能求解出每一类地物（也就是一个端元束内所有光谱的比例之和）的最小值和最大值，再对其作平均得到每类地物的比例。

该方法的优点在于可以得到每类地物比例的误差范围。

多端元方法机制明确，但计算复杂，耗时过长。

(2) 光谱变换
在很多情况下，同类地物的光谱的差别来自绝对值的变化，而光谱形状是相似的。

因此通过对光谱进行一定的变换可以减少端元的光谱差异。

Wu（2003）[56]提出将光谱除以各个波段的均值，再作混合像元分解，并应用于城市监测；Garcia-Haro等（2005）[57]将光谱作标准化后再作混合像元分解；Asner等(2003)[58]将光谱作微分后再作混合像元分解。

Juan Pablo Guerschman等(2009)[59]利用原始光谱计算出归一化差分植被指数（Normalized Difference Vegetation Index, NDVI）和纤维素吸收指数（Cellulose Absorption Index，CAI），假设两个指数也满足线性混合模型，利用两个指数求得光合植被、非光合植被及裸土的比例。

由于这两个指数只利用了光合植被及非光合植被的特征波段，因此端元光谱差异也能得到一定的压缩。

变换光谱的方法虽然计算方便，但在光谱变换中线性混合模型假设被破坏，使得混合机理不清楚，由此可能带来新的误差。

(3) 基于概率的方法
Song(2005)[60]提出了一种贝叶斯混合光谱分析方法（Bayesian Spectral Mixture Analysis, BSMA）, BSMA不再将端元光谱用一个值表示，而是用概率密度函数来描述端元光谱的分布，进而根据贝叶斯理论推导出最似然的端元比例。

然而，模拟实验与真实数据的检验并未证实该方法相对于FCLS的优越性。

(4) 基于光谱匹配的方法
陈晋等(2009)[61]提出了基于光谱匹配的混合像元分解技术将混合光谱与端元按照其比例线性混合的光谱进行匹配，以匹配效果最好的端元比例作为分解的最终结果。

这种算法能够在一定程度上克服由于大气或地形因素影响引起的端元光谱变异。

精度评价
混合像元分解结果等同于软分类结果，目前一般采用软分类结果的评价方法。

最简单且应用最广泛的精度评价是均方根误差RMSE。

RMSE计算了某类端元的分解总体精度，但对误差来源的错分、漏分误差没有进行区别。

借鉴硬分类中基于混淆矩阵的精度评价方法，近年来一些学者试图利用混淆矩阵对混合像元分解结果进行精度评价。

混合像元分解结果的混淆矩阵如表1所示，由于混合像元分解不能直接提供端元比例在混合像元中空间分布信息，因此需要通过一些算子估计混淆矩阵中的元素。

比较有代表性的算子有：取小算子（MIN）[62]、相似度指数（Similarity Index, SI）[63]和乘法算子（Product Operator, PROD）[64]。

取小算子基于模糊集理论，得到最理想情况下的精度估计；相似度算子引入标准化过程，使得矩阵元素落在0到1之间；乘法算子基于概率的思想，得到概率意义下的期望精度估计。

表2
归纳了这三个算子的计算方式及物理意义。

Silvan-Cardenas 等(2008)[65] 指出用于混淆矩阵的算子需要满足一些约束条件，并提出SCM评价方法（Sub-pixel confusion–uncertainty matrix），该算法引入复合算子，分别在混淆矩阵的对角元和非对角元使用不同的算子，不仅给出错分、漏分的精度，同时给出误差不确定性范围。

但是，到目前为止，没有直接的证据表明哪一种精度评价方式更优越。

参考文献
[1]Keshava N. and Mustard J.F., Spectral unmixing[J]. IEEE Signal Processing Magazine, 2002. 19(1): 44-57.
[2]Phinn S., Stanford M., Scarth P., Murray A.T., Shyy P.T., Monitoring the composition of urban environments based on the vegetation-impervious surface-soil (VIS) model by subpixel analysis techniques[J]. International Journal of Remote Sensing, 2002. 23(20): 4131-4153.
[3]Small C., Multitemporal analysis of urban reflectance[J]. Remote Sensing of Environment, 2002. 81(2): 427-442.
[4]Kameyama S., Yamagata Y., Nakamura F., Kaneko M., Development of WTI and turbidity estimation model using SMA - application to Kushiro Mire, eastern Hokkaido, Japan[J]. Remote Sensing of Environment, 2001. 77(1): 1-9.
[5]Rogan J., Franklin J., Roberts D.A., A comparison of methods for monitoring multitemporal vegetation change using Thematic Mapper imagery[J]. Remote Sensing of Environment, 2002. 80(1): 143-156.
[6]Riano D., Chuvieco E., Ustin S., Zomer R., Dennison P., Assessment of vegetation regeneration after fire through multitemporal analysis of AVIRIS images in the Santa Monica Mountains[J]. Remote Sensing of Environment, 2002. 79(1): 60-71.
[7]Pinet P.C., Shevchenko V.V., Chevrel S.D., Daydou Y. Rosemberg C., Local and regional lunar regolith characteristics at Reiner Gamma Formation: Optical and spectroscopic properties from Clementine and Earth-based data[J]. Journal of Geophysical Research-Planets, 2000.
105(E4): 9457-9475.
[8]Adams J.B., Smith M.O., Gillespie A.R., Simple Models For Complex Natural Surfaces: A Strategy For The Hyperspectral Era Of Remote Sensing[C]. IGARSS'89. 12th Canadian Symposium on Remote Sensing. 1989. 1.。