混合像元分解研究综述——线性混合像元分解算法

本文转载来自：http://hi.baidu.com/hilbertspace/blog/item/aae6a20f2cb564e4ab6457c5.html

--------------------------------------------------------------------------------------------------------------------------------------------------------

由于线性模型是应用最广泛，也是研究最多的算法，下面重点介绍基于线性模型的混合像元分解算法。一般而言，混合像元分解算法包括数据降维、端元选取和反演三个步骤。

1.数据降维

尽管数据降维不是混合像元分解算法的一个必需步骤，但由于大多数算法都将其作为一个流程，我们也将其当作一个步骤。常用的降维算法有主成分分析（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。考虑到端元比例的物理意义，需要满足约束条件Σf_i=1，则可通过拉格朗日乘子法求解部分约束最小二乘法（abundance sum-to-one constraint，ASC），也有解析解。Heinz等(2001)^[46] 提出全约束最小二乘（fully constrained least squares，FCLS）的求解方法。FCLS在原方程组里添加方程Σf_i=1，并给予其很大的权重，再用ANC对新的方程组求解，则求得的解近似满足两个约束条件。由于该方法能够同时满足两个约束条件以及高效的运算效率，因而得到广泛认可。但是也有学者认为，准确的混合像元模型的分解结果应当自动满足两个约束条件，强加的约束条件可能会导致新的误差^[47,48]。