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 豐度估計結果比較 單位：%