osp分析方法

① 混合像元分解

为了实现在无先验信息情况下的高光谱数据的端元提取，研究过程中提出了基于高光谱数据高阶统计量的改进独立成分分析（Independent Component Analysis，ICA）端元提取方法和基于扩展形态学与OSP的自动端元提取算法。

在端元提取的基础上，提出了基于信息散度的光谱混合分析方法，与基于混合调制匹配滤波的丰度估计方法相比，能够实现矿物组分的更精确估计。

6.3.1 改进ICA 模型端元提取

对于光谱解混问题，若端元光谱已知，则问题变得十分简单，可以通过极大似然、光谱信号/特征匹配、SAM、子空间投影等方法求解；但是在大多数情况下，物质的数据和它们的反射率是非已知的，光谱解混问题就转换成了盲源信号处理问题。因此，若不考虑噪声影响，ICA模型能够较好地实现。

（1）改进ICA算法的实现

ICA用于高光谱数据端元提取时的两个前提假设是：各种端元线性混合成观测到的信号及源信号（丰度）统计上独立。对于第一个假设，是比较容易满足的，因为目前线性混合模型就是在这个前提假设下建立的，并且通过前面的线性混合模型的物理机理分析也容易看出，实际问题中大多数情况下是满足的。但是对于第二个假设是很难满足的。因为从混合的物理意义上讲，每一种地物的丰度都应当大于0并且一个像元内各端元丰度之和等于1。

因此，如何将标准ICA模型改进，使其适合于高光谱数据混合像元分解是目前该领域的热点与难点。

首先分析了应用于混合像元分解的ICA模型与标准ICA模型的不同：标准ICA模型中的分离矩阵不必要是满秩且正交的。为了实现改进ICA的自动端元提取算法，在设计模型估计的学习算法时基于丰度的独立性考虑而不是针对分离矩阵考虑。在考虑丰度非负性且总和为1的约束的同时，采用扩展的信息最大化学习方法进行独立成分和混合矩阵的估计。

在非高斯性度量中，采用峰度的绝对值作为非高斯性的度量指标，即

高光谱遥感技术原理及矿产与能源勘查应用

为了极大化峰度的绝对值，可以从某个向量w开始，依据可观测的信号x₁ ，x₂ ，…，x_m（假设数据已经中心化和白化），计算出使y＝b^Tx峰度绝对值增大最快的方向，然后将向量b更新，转到该方向上，并标准化w，使‖b‖² ＝1。该操作可以利用梯度法及其扩展操作来实现，则有

高光谱遥感技术原理及矿产与能源勘查应用

则得到下面的梯度算法：

Δb∝E｛x（b^Tx）³｝-3b （6.26）

令式中峰度的梯度与b相等，则有

b∝E｛x（b^Tx）³｝-3b （6.27）

写成矩阵形式为

B←E｛X（B^TX）³｝-3B （6.28）

为了满足矩阵B中每一个分量‖b‖² ＝1，则需要对矩阵B进行去相关和标准化，实现方法如下：

B←（BB^T）^-1/2B （6.29）

因此，分离矩阵W＝B′× whitening_matrix可以计算得到，其中whitening_matrix是可观测信号X的白化矩阵，从而得到各个独立成分。

（2）算法应用与实验分析

利用本方法对图6.20所示的东天山局部区域反射率数据（基于改进辐射传输模型的大气校正方法获得）进行端元提取，提取结果如图6.21所示。通过分析该地区地质条件、主要矿物等，可以看出，本方法提取得到了该地区主要的矿物，包括：绿泥石、绿帘石、白云母、蛇纹石、方解石等九种蚀变矿物，并且本方法提取的是图像中存在的矿物光谱，端元光谱正确率较高。

图6.20 东天山局部区域反射率数据

图6.21 基于ICA方法的端元提取结果

利用混合比例20%～100% 混合矿物的解混结果（丰度估计结果）进行端元提取正确性的验证，如表6.21所示，表明丰度估计的均方根误差为0.019。

表6.21 目标丰度估计结果 单位：%

6.3.2 基于扩展形态学与OPS 的自动端元提取

光谱和空间信息同时利用进行高光谱数据建模能够提供高光谱数据处理的精度、可靠性和稳定性。因此，利用数学形态学理论将腐蚀和膨胀操作扩展到高光谱数据端元提取处理技术中，目前利用数学形态学理论进行端元提取的方法主要是Chang研究小组A.Plaza提出的自动数学形态学端元提取算法AMEE，该算法能够较好地得到图像中的端元，但是不同端元很难区分。因此，该部分主要针对这一问题引入了正交子空间的概念，对算法进行了改进。

6.3.2.1 算法原理

（1）正交子空间投影

正交子空间投影（Orthogonal Subspace Projection，OSP）主要是针对混合像元分解问题提出来的。在利用线性解混合模型进行混合像元分解时，假设像元向量r的光谱特征是m₁ ，m₂ ，…，m_p 的线性混合，则

r ＝ Mα ＋ n （6.30）

式中：端元光谱信号M ＝［m₁，m₂，…，m_p］；丰度矩阵α＝［α₁，α₂，…，α_p］；噪声或模型误差为n。

可以将上式继续分解为

r ＝ dα_p ＋ Uγ ＋ n （6.31）

式中：d＝m_p是感兴趣端元的光谱信号；U＝［ m₁ ，m₂ ，…，m_p-1 ］是非感兴趣端元的光谱信号。丰度矩阵则相应的分解，上式为OSP的基础。

通过设计正交子空间投影算子

实现非感兴趣端元或者干扰影响的抑制，由最小二乘估计得到投影算子

-UU^＃，其中，U^＃＝（U^TU）^-1U^T。

将

应用到上式，则得到OSP模型：

高光谱遥感技术原理及矿产与能源勘查应用

这里非感兴趣端元U已经被抑制，并且原来的噪声n也被压缩为

。

然后引入一个滤波算子w^T，为l × L向量，通过将输出信号信噪比最大化实现端元的提取。将w^T应用于OSP模型得到

高光谱遥感技术原理及矿产与能源勘查应用

由上式信噪比SNR可以得到

高光谱遥感技术原理及矿产与能源勘查应用

式中：σ_n是噪声的标准差；

是上式得一个最优解，则λ_max必是w^*其中之一。通过构造匹

配滤波算子得到OSP检测算子：

高光谱遥感技术原理及矿产与能源勘查应用

（2）扩展数学形态学操作

高光谱图像处理中，为了确定根据目标与背景差异的多维向量的排序关系，引入一个多维向量的度量算子，该度量算子由结构元素内各个像素累加距离计算得到，定义如下：

高光谱遥感技术原理及矿产与能源勘查应用

式中：dist是测量N维向量的逐点线性距离。为了有效地利用多/高光谱数据提供的光谱和空间信息，OPD方法源于正交子空间投影的概念和原理，计算得到的是正交投影后的残余。考虑到SAM只是从光谱波形形状出发，并且对噪声敏感，研究可采用OPD计算该距离，考虑两个N维光谱信号s _i ＝［s_i1，s_i2，…，s_iN］^T，s_j＝［s_j1，s_j2，…，s_jN］^T，则N维光谱信号s_i 和s_j之间的OPD表示为

高光谱遥感技术原理及矿产与能源勘查应用

其中，

＝I_N×N-s_k

，k＝i，j，并且I_N×N是N×N维的单位矩阵。

因此，累加距离D能够根据像元纯度的差异大小排序结构元素中的向量。根据以上定义和叙述，多/高光谱数据中腐蚀和膨胀操作分别定义如下：

高光谱遥感技术原理及矿产与能源勘查应用

式中arg_Max，arg_Min分别表示使得累加距离D达到最大和最小的像素向量。通过以上的分析表明，扩展到多/高光谱图像的膨胀结果得到的是在结构元素内纯度最大的像元，腐蚀结果得到的是在结构元素内混合度最大的像元。

（3）改进的端元自动提取模型

通过从原理上分析，AMEE方法能够提取得到纯端元，利用Ostu自动阈值分割方法很难区分不同类型的端元，尤其是光谱特征相似的端元，存在着较大的分离误差。因此，针对这一问题，在分析正交子空间投影OSP原理的基础上，引入向端元子空间投影的概念和理论，实现光谱和空间信息的端元提取。基于光谱和空间信息的端元提取改进模型实现的过程如图6.22所示。

基于线性混合模型，利用OSP概念进行混合像元分解时，仍可以继续将端元矩阵分解，则混合模型可以表示为

r ＝ dα_d ＋ Uα_U ＋ n （6.40）

式中：d为第一个端元光谱信号，通过MEI图像中提取的端元数据，利用SAM方法得到其中一个端元光谱信号d＝e₁，然后利用向其正交子空间投影得到

来消除端元数据中的信号d＝e₁，其中

＝I-d（d^Td）^-1d^T。将

代入基于OSP的混合模型中得到

，可以看出，不仅消除了d＝e₁的影响而且有效地压制了噪声

。然后把

视为新的端元数据集，继续提取端元，经过m-1次投影之后得到d＝ ｛e₁ ，e₂ ，…，e_m｝，投影循环结束的条件为得到指定个数的端元或是由分解误差RMSE决定，从而得到了m个端元。

图6.22 光谱与空间信息结合的自动端元提取方法

该方法有效地克服了AMEE算法从数学机理上无法将端元区分的问题，提高了端元提取精度，实现了自动的端元提取。

6.3.2.2 应用实验与结果分析

应用美国内华达州Cuprite矿区的AVIRIS高光谱数据进行基于扩展数学形态学和正交子空间投影方法的自动端元提取研究。提取的四种主要蚀变矿物端元光谱如图6.23所示。

可以看出，图6.23得到了该矿区内四种典型的矿物端元光谱，通过与USGS光谱库数据比较，图6.23给出的提取的矿物光谱与参考光谱比较一致，说明了该方法的有效性和正确性。

该方法充分利用了高光谱数据提供的空间和光谱信息，并通过扩展数学形态学的理论将二者有效的结合，并综合利用，同时利用OSP的原理有效的区分不同类型的端元光谱，试验证明得到了较好的结果。

图6.23 基于扩展数学形态学和正交子空间投影方法得到的端元光谱

6.3.3 基于信息散度的光谱混合分析

该方法从总的端元组出发，每进行一次迭代循环，计算一个均方根误差，并去掉丰度最小的那个端元，进入到下一个循环，直到端元组中只剩下一个端元，停止循环；然后根据得到的均方根误差曲线，根据rms变化率准则来判定最优的端元子集。并且在实现过程中算法加入了端元的初选和二次选择，利用光谱信息散度（SID）作为最优端元组判定的准则，能够达到很好的端元选择效果和丰度估计精度。实现流程图如图6.24所示。

SID-SMA主要由三部分组成：端元的初选、端元的二次选择和最终的丰度估计。端元的初选是利用某些规则先去掉一些端元，然后进入迭代循环，再去掉一些端元，最后留下的那些端元我们就认为是某个像元真实存在的物质，最后利用这些端元对该像元进行丰度估计，实现光谱解混。

端元的初选有两种规则，根据数据信噪比的不同而不同。当数据的信噪比较高时，采用线性逆卷积的规则，也就是对原始的端元进行最小二乘估计，去掉丰度小于0的端元，再利用剩下的端元重新进行最小二乘估计，这样反复循环，直到最后最小二乘估计的丰度值没有负数为止。这种规则能够在较少的循环内去掉较多的无用端元，并且在信噪比较高时，端元选择的精度较高，这样既能提高算法的速度，又能取得好的精度。但当信噪比较低时，利用线性逆卷积的规则容易造成正确端元的遗漏，所以当信噪比降低到一定的值时，我们采用全限制最小二乘法来进行端元的初选。全限制最小二乘法首先计算初始端元的最小二乘估计

，然后计算端元的全限制最小二乘估计

-λs，其中λ＝

（1^Ts），s＝（M^TM）^-11，M是端元组矩阵，1为单位列向量。观察α_FCLS是否含有负值，如果有负值，则计算|α_FCLS，j/s_j|，去掉绝对值最大的对应的端元，然后循环，直到

的元素都大于0为止。利用全限制最小二乘法进行端元初选每次循环只能去掉一个端元，效率比线性逆卷积要低，但是它在数据信噪比比较低时选择端元的正确率较高，遗漏的端元较少，因此，我们在信噪比低时使用它作为端元的初选。

初选之后，进入二次选择过程。二次选择时先利用剩下的端元组对混合光谱进行最小二乘估计， 利用估计的丰度值及相应的端元重新建模得到建模光谱，计算建模光谱与原光谱的信息散度（SID），然后去掉丰度最小的那个端元，得到较小的端元组，进入下一个循环，这样，直到端元组里只剩下一个端元为止，停止循环。到循环结束时会得到一条SID曲线，根据这条曲线能够判定算法进行到第几个循环后留下的是正确的端元。最后利用这些正确的端元对混合光谱进行最小二乘估计，得到最终的光谱解混结果。

图6.24 SID-SMA流程图

这里的光谱信息散度衡量的是两条光谱之间的信息差异。假设两条光谱分别是x＝（x₁，x₂，…，x_l）^T和y＝（y₁，y₂，…y_l）^T，可以得到两条光谱的概率向量分别是p＝（p₁，p₂，…，p_l）^T 和q＝（q₁， q₂，…，q₁）^T，其中p_i＝x_i/

，q_i＝y_i/

。根据信息理论，我们可以得到x和y的自信息为

I_i（x）＝-logp_i和I_i（y）＝-logq_i

通过上式，可以得到y关于x的相对熵：

高光谱遥感技术原理及矿产与能源勘查应用

同理可得x关于y的相对熵：

高光谱遥感技术原理及矿产与能源勘查应用

而x和y的光谱信息差异为

SID（x，y）＝ D（x‖y）＋ D（y‖x） （6.43）

SID是利用光谱信息造成的相对熵对两条光谱进行相似度测评的度量，它的效果要好于光谱角度调制（SAM）。另外在二次选择的时候，我们使用的SID判定准则是两个循环SID的变换量：

ΔSID ＝ SID_it-SID_it-1 （6.44）

当变化量大于某个阈值时，我们认为已经有正确的端元被排除出端元组，而这个循环之后剩下的端元都是正确的端元。在实际应用时，由于正确的端元个数相对整个端元组总是比较少的，所以我们判定最佳端元组时，往往从最后一个循环往前推，当SID的变化量小于某个阈值时，则认为在这个循环以后正确端元开始被排除出端元组。

为了验证SID-SMA方法的性能，利用USGS光谱库数据进行不同信噪比、不同端元数对算法性能的影响，丰度估计的误差由式（6.45）给出：

高光谱遥感技术原理及矿产与能源勘查应用

式中：n为端元组中端元的个数，试验中选取了USGS中的29条光谱，故设为29；m为混合光谱的条数；a_ij为预先设定的第i条混合光谱第j个端元的丰度值；e_ij为计算得到的第i条混合光谱第j个端元的丰度值。

在进行算法验证时，针对混合光谱由3～10个端元合成时的情况，每种情况都有1000条混合光谱，每次试验时从29个端元中随机选取。并且利用全限制最小二乘（FCLS）、线性卷积（LD）、迭代光谱混合分析（ISMA）及SID-SMA方法进行性能分析与验证。

表6.22是四种算法端元选择精度的比较，其中三个比较参数分别是用于最终解混选择的端元个数（实际的个数为6个）、端元选择的正确率、正确端元遗漏的个数。从中可以看出，当噪声比较小的时候，LD算法的端元选择性能比较好，正确端元遗漏比较少，而且端元选择的正确率比FCLS要高，但是当信噪比下降到25和12时，LD遗漏的端元就要比FCLS要多了，这也会影响到LD丰度估计的精度。而SID-SMA与ISMA相比，虽然在端元选择的正确率上比不上ISMA，但是SID-SMA遗漏的端元要比ISMA少得多，特别是信噪比为12 的时候，ISMA遗漏的端元达到了3.43 个，也就是说已经漏掉了实际端元的一半多，而它的较高的端元选择正确率也是建立在此之上，这大大影响了ISMA在低信噪比时丰度估计性能。

表6.22 四种算法的精度比较

图6.25是四种算法的丰度估计误差随信噪比变化而变化的情况。从中可以看出，无论信噪比为多少，丰度估计误差最小的都是SID-SMA，当信噪比为100时，SID-SMA的误差小于0.1。而LD算法在噪声比较小的时候性能比较好，但当信噪比下降到25时，丰度估计误差已经大于FCLS，说明LD比较适合在高信噪比的情况下使用。而ISMA正如我们上面所说的那样，当信噪比下降时，它由于遗漏太多的正确端元，从而导致估计的误差急剧增大，到信噪比下降到12时，估计误差已经达到0.9，明显高于其他三种算法。FCLS则在信噪比比较低的时候性能稍微差，但是它的抗噪性能比较强。

图6.26给出了SNR为100∶1情况时LD，FCLS，ISMA和SID-SMA的估计误差随参加混合端元数的变化而变化的情况。当参加混合的端元数增加时，算法的整体精度下降。而ISMA 的精度下降的最快。当实际端元数量2～4 个时，ISMA 的误差较小，基本与 SID-SMA 相仿，并且要明显小于FCLS，但是当端元数增加到10个时，ISMA的误差已经超过其他两种算法很多。而SID-SMA算法在端元数2～10个时，整体精度都很高，不过当实际端元数达到9～10个时，误差有加速增大的趋势。但是，从整体上来讲，SID-SMA的表现是最好的。

图6.25 四种算法的不同SNR情况下丰度估计整体误差

图6.26 四种算法的丰度估计误差随实际端元个数变化情况（SNR100∶1）

利用SID-SMA方法进行Cuprite矿区主要蚀变矿物明矾石、高岭石、热液硅石及布丁石的丰度估计结果与迭代SMA的结果进行比较如图6.27所示。从图6.27分析可知，ISMA得到的矿物丰度图中有大量的点的丰度小于0 ，四种矿物丰度小于0 的点的比例分别是1.55%，1.48%，4.86% 和5.68%，而SID-SMA得到的丰度值没有负数值。此外，虽然对于明矾石与高岭石矿物的丰度结果相似，但是对于布丁石、热液硅石等吸收特征较宽且不明显的矿物丰度估计结果SID-SMA方法明显优于ISMA方法。

利用图6.20所示的东天山局部区域反射率数据进行基于SID-SMA方法的矿物组分含量估计，结果如彩图6.1所示。

图6.27 由SID-SMA（每组左侧）和ISMA（每组右侧）得到的矿物丰度图

（a）高岭石；（b）明矾石；（c）布丁石；（d）热液硅石

6.3.4 基于混合调制匹配滤波的丰度估计

由于在自然环境中，线性混合模型将受到两个约束条件的限制，这两个约束条件制约着混合系数（即端元丰度）的大小，分别为非负性约束和归一化约束。它们的物理意义非常明显：光谱是能量的表现，不可能存在负值；混合能量的大小是存在限定的，不可能无限的大，从线性混合模型物理机理分析可以看出，混合能量满足归一化约束。为此，在端元提取后的丰度估计中，考虑在上述约束条件下进行端元的丰度估计。

匹配滤波通过最大化已知端元波谱信号，压制未知复合背景的响应信号实现端元波谱匹配，该方法能生成类似于波谱分离的影像，但由于不需要已知所有的端元光谱而使计算量显着降低。

混合调制匹配滤波技术是匹配滤波技术和线性混合分解理论的复合方法。该方法将上述描述的匹配滤波不需要已知其他背景端元波谱的优点与线性混合分解理论中的物理条件限制（给定像元的信号是包含在该像元中的单一物质成分的线性组合，同时各个组分的含量为正且和为1）结合起来，因而提高了矿物的检出限，能探测出其他方法不能检测出的岩石中微量的矿物成分。混合调制匹配滤波的结果为灰度在0～1.0之间的匹配滤波图像，反映了参考波谱的相对匹配程度，即相对丰度图像。其中，1.0代表完全匹配，即相对参考波谱丰度为1。

利用该方法进行混合比例40%～100% 混合矿物的丰度估计误差分析，得到结果如表6.23所示，得到丰度估计的均方根误差为0.12。

表6.23 丰度估计结果比较 单位：%