基于壓縮感知的混合光譜解析算法

伍 娟 妮

(日日升智能科技發(fā)展(山東)有限公司 山東 煙臺 264006)

0 引 言

壓縮感知(Compressed Sensing)[1-17]是一種近年發(fā)展起來的信號處理技術，它利用信號的稀疏性，通過求解欠定線性矩陣方程的解來有效地恢復重建信號。Donoho等科學家于2006首先提出壓縮感知的理論并給出了相關的證明。壓縮感知理論的出現(xiàn)引起了學術界的極大興趣，同時工業(yè)領域也表現(xiàn)出積極的熱情。壓縮感知理論研究的不斷推進，帶動了信息論、信號編碼理論和計算機科學發(fā)展，也為工程應用領域帶來新的技術支撐，進一步促進機器學習、模式識別、圖像處理、生物信息處理、醫(yī)學數(shù)據(jù)處理和高光譜遙感等領域的發(fā)展。

本文概要地介紹了范數(shù)、線性系統(tǒng)矩陣方程及常見求解算法[13，19]、稀疏度(Sparsity)[1]、受限等距性(Restricted Isometric Property,RIP)[1]和受限零空間性性質(Restricted Null space Property)[1]等壓縮感知問題的基本概念和數(shù)學模型。從經(jīng)典線性代數(shù)的角度來看，壓縮感知問題可以簡化為尋找如下欠定線性矩陣方程組的稀疏解的問題[1]：

y=Axx∈Rn×1

(1)

式中：A∈Rm×n(m

(2)

s.t.y=As

此L0范數(shù)優(yōu)化問題是一個求解組合優(yōu)化的NP難問題，通常求解很困難。因此，要考慮此問題的L1松弛：

(3)

s.t.y=As

相關理論證明了當且僅當矩陣A滿足受限零空間性質時，L0范數(shù)問題與L1范數(shù)問題等價[1]。

針對L0范數(shù)和L1范數(shù)優(yōu)化問題，本文重點討論了三種恢復稀疏信號的求解算法。其中正交匹配追蹤(Orthogonal Matching Pursuit,OMP)[5]算法是一個近似求解L0問題的貪婪算法，其復雜度低且容易簡單。著名的最小角回歸(Least Angle Regression,LARS)[5-6]算法求解L1范數(shù)問題。分段正交匹配追蹤(Stage-wise Orthogonal Matching Pursuit,StOMP)[20]算法與迭代收縮算法存在廣泛的相似之處，并且結合了OMP和LARS的優(yōu)點。

數(shù)值實驗部分介紹了高光譜圖像處理和生物信息學以及化學信息學中的混合光譜模型。混合光譜解析[21-25]是近十年迅速發(fā)展的領域，也是高光譜分析領域和生物信息學光譜處理的重點和難點?；趬嚎s感知的混合光譜分解算法利用已知的光譜庫作為基礎，代替了端元集合用于高光譜混合像元分解的方法。近年來，隨著稀疏表示和壓縮感知理論成功地應用于機器學習、圖像處理和信號處理領域，研究人員對信號的稀疏性有了更深刻的認識。在高光譜分析領域，Iordache等[26]首次嘗試將稀疏性約束加入到高光譜混合像元分解模型中，自此之后，基于稀疏性的高光譜混合像元分解受到國內外學者的高度關注，并成為新的研究熱點。

通過兩種不同光譜庫的具體數(shù)值實驗，分析了幾種稀疏恢復算法在混合光譜解析時的性能和存在的問題，并給出了相應的優(yōu)化建議。

1 基礎概念

1.1 向量范數(shù)

向量x=(x1,x2,…,xn)的Lp范數(shù)定義為：

(4)

常見的有L0范數(shù)、L1范數(shù)、L2范數(shù)和無窮(極大)范數(shù):

(5)

(6)

(7)

(8)

1.2 線性模型

線性參數(shù)估計廣泛應用于科學和工程技術領域。線性模型可以表示為如下的線性矩陣方程形式：

Ax=y

(9)

式中：數(shù)據(jù)矩陣A∈Rm×n，假設y∈Rm×1，則x∈Rn×1。數(shù)據(jù)向量y和數(shù)據(jù)矩陣A的不同，矩陣方程可以有如下四種主要的類型：

(1)m=n。即方程個數(shù)等于未知量個數(shù)，矩陣A和向量y均已知且A非奇異，則方程有唯一解x=A-1y。

(2)m>n。稱為超定方程(又稱矛盾方程)，并且矩陣A和向量y均已知，其中之一或者二者可能存在誤差或干擾，此時，方程組沒有任何“嚴格”的解。

(3)m

(4) 盲矩陣方程。數(shù)據(jù)矩陣A未知，向量x未知，僅數(shù)據(jù)向量y已知。

矩陣方程的高效求解算法一直是數(shù)學和計算機算法領域研究的一個核心，本文概述上述矩陣方程的常見求解方法并重點說明欠定矩陣方程的適用條件和常用解法。

1.2.1超定方程m>n

超定矩陣方程中方程個數(shù)大于未知量個數(shù),最常用的求解方法是最小二乘法(Least Square,LS)。根據(jù)Gauss-Markov定理，如果假設量測誤差為高斯正態(tài)分布且具有零均等方差，則最小二乘估計量是唯一的一個具有最小方差的線性無偏估計量[19]。

針對矩陣A和向量y的不同情況，分別有普通最小二乘法、數(shù)據(jù)最小二乘法、Tikhonov正則化方法、總體最小二乘法和約束最小二乘法等多種算法[19]。

最小二乘法的核心思想就是尋找解向量x使得矩陣方程兩邊的誤差平方和最小。矩陣方程Ax=y的普通最小二乘解為：

(10)

若矩陣A列滿秩，即rank(A)=n,則ATA非奇異，矩陣方程有唯一解：

xLS=(ATA)-1ATy

(11)

若rank(A)

xLS=(ATA)-1ATy

(12)

式中:(ATA)-1為ATA的Moor-Penrose逆矩陣。

在工程應用中，數(shù)據(jù)矩陣A常為秩虧損矩陣，即rank(A)

(13)

式中:λ為正則化系數(shù)。其解為：

xLS=(ATA+λI)-1ATy

(14)

在信號處理領域該方法稱為松弛法。在機器學習和統(tǒng)計學中,稱該L2約束問題為嶺回歸。

Tikhonov正則化方法實質上是一種改良的最小二乘估計法，通過在矩陣ATA上加入一個懲罰項λI使得矩陣ATA+λI非奇異。它放棄了最小二乘法的無偏性，以損失部分信息、降低精度為代價提高模型的穩(wěn)健性，獲得的回歸系數(shù)更符合實際、更可靠，對病態(tài)數(shù)據(jù)的擬合要強于最小二乘法。

如果考慮矩陣A和向量y都存在觀測誤差，并且假設誤差是具有相同方差且滿足獨立同分布的隨機變量，則常采用總體最小二乘(Total Least Squares,TLS)法。如果噪聲統(tǒng)計相關或具有不同的方差，則考慮使用約束總體最小二乘法，具體內容參閱文獻[19]。

1.2.2盲矩陣方程

考慮盲矩陣方程:

X=AS

(15)

式中:X∈Rm×n為觀測矩陣;矩陣A∈Rm×d、S∈Rd×n均未知。若假設A列滿秩、S行滿秩，則可在已知X的情況下，求解未知矩陣A和S。求解盲矩陣方程的常用方法有：子空間法、非負矩陣分解法和獨立成分分析等，具體內容參閱文獻[19]及其參考文獻。

1.2.3欠定矩陣方程m

欠定矩陣方程中方程的個數(shù)小于未知量個數(shù)，適用于樣本數(shù)小于特征數(shù)的情況。此時需要考慮三個問題：該方程是否存在稀疏解，該稀疏解是否唯一，如何求解。

2 壓縮感知

2.1 概 述

壓縮感知理論認為：通過信號恢復算法，可以利用信號的稀疏性從比Shannon-Nyquist采樣定理所需要的信號少得多的樣本中恢復信號。壓縮感知理論中的兩個最關鍵的概念就是稀疏度和不相關性。稀疏度要求信號在某個域中是稀疏的，不相關性通過等距特性度量。一般而言，相關性越小，恢復算法的性能越好。

(2) 受限零空間性質。給定矩陣A∈Rm×n，其零空間為：

null(A)={x∈Rn|Ax=0}

(16)

對于給定的子集S∈{1,2,…,n}，若有：

null(A)∩C(S)={0}

則稱矩陣A∈Rm×n在S上滿足受限零空間性質，記為RN(S)。

(3) 受限等距性(Restricted Isometric Property,RIP)。最簡單度量方式是讓矩陣A中的各列做內積，一個更復雜的不相關概念與受限等距性(RIP)[5，8，27]有關。

對于線性系統(tǒng)Am×nx=y,m>n，若存在常數(shù)δk∈(0,1)，使得對于任意向量s∈Rk和A的子矩陣Ak∈Rm×k有:

(17)

則稱矩陣A滿足k限定等距性(kRIP)。此時可以通過下面的優(yōu)化問題近似完美地從y中恢復稀疏信號s，進而恢復出信號x：

(18)

s.t.y=As

此L0范數(shù)最小化問題是一個求解組合優(yōu)化的NP難問題，通常求解很困難。因此，要考慮此問題的L1松弛：

(19)

s.t.y=As

當且僅當矩陣Am×n滿足受限零空間性質時，L0與L1問題等價[5]。

2.2 數(shù)學模型

壓縮感知理論同時在信號處理和統(tǒng)計學習領域獨立發(fā)展，最后得到近似的概念[5]。從經(jīng)典線性代數(shù)的角度來看，壓縮感知問題可以簡化為尋找如下欠定線性系統(tǒng)的稀疏解的問題：

y=Axx∈Rn×1

(20)

式中：A∈Rm×n(m

圖1 CS模型

圖2 Lp模型

(1) 當0

(2) 當p=1時，Lp球是菱狀球，當球與直線Ax=y相交于坐標軸上時得到一個稀疏解。

(3) 當p>1時，Lp球是外凸球，直線Ax=y與外凸球切點一定不在坐標軸上，因而，此時的解是不稀疏解。

3 算 法

欠定稀疏矩陣方程是稀疏表示和壓縮感知領域需要求解的問題，其求解方法有貪婪追蹤、凸松弛方法、迭代收縮算法、貝葉斯框架和消息傳遞/置信傳播算法[1，17]等。

1) 貪婪算法是稀疏信號處理領域中廣泛使用的一種信號復原法算法，此類算法在某些情況下優(yōu)于凸集選擇算法[3]。貪婪算法一般分為貪婪追蹤算法和閾值算法。貪婪追蹤算法通常以零值為初始值，在每一步迭代中通過最優(yōu)化方法得到確定的非零值，經(jīng)過多輪迭代，逐步建立信號的估計。這類方法通常能非?？焖俚靥幚磔^大的數(shù)據(jù)，但是其理論性能比常規(guī)的凸集方法弱。常用的貪婪追蹤算法有：匹配追蹤(MP)、正交匹配追蹤(OMP)、梯度追蹤(GP)、共軛梯度追蹤(CGP)、分布正交匹配追蹤(StOMP)、分布弱梯度追蹤(StWGP)和順序遞歸匹配追蹤(ORMP)等。閾值類貪婪算法在迭代選擇非零元素的同時剔除非零元素。此類算法通常簡單且運行速度快而且有與凸集算法相當?shù)睦碚撔阅?，甚至能夠用來恢復非稀疏信號，但是計算精度不高。常用的閾值類算法有：迭代硬閾?IHT)、壓縮采樣匹配追蹤(CoSaMP)、子空間追蹤(SP)等。

(2) 考慮到L0與L1問題的受限等價性，在求解欠定矩陣方程的稀疏解時往往使用更容易實現(xiàn)的L1凸松弛算法。Lp優(yōu)化問題在最優(yōu)化領域有著廣泛而深入的研究和大量的算法，比較常用的最小化L1算法有：Homotropy算法(最小角回歸算法LARS)、Chabollet-Pock算法[3]，F(xiàn)OCUSS(Focal Undetermined System Solver)算法(又稱為迭重加權最小二乘法(Iterative Reweighed Least Squares))[11]等。更深入的研究請參閱相關專著。

4) 稀疏貝葉斯學習(Sparse Bayesian Learning,SBL)的相關理論最早出現(xiàn)在Tipping關于相關向量機(Relevance Vector Machine,RVM)的文獻中，之后RVM模型被引入到壓縮感知信號恢復中[10]。文獻[16]中總結到：SBL算法與廣泛使用的L1懲罰算法(如LASSO、BP等)相比優(yōu)勢明顯，具體體現(xiàn)在以下幾個方面：(1) 在無噪情況下,如果不能滿足某些嚴格的條件，L1懲罰算法的全局最小點并不是真正的最稀疏解，而此時SBL可以更好地得到最稀疏的真實解。(2) 當感知矩陣A的列與列之間相關性很強時，大部分已知算法的性能都比較差，而SBL性能良好。(3) 學者已經(jīng)證明SBL算法等價于一種迭代加權L1算法，更易獲得真正的稀疏解。(4) 當信號不是特別稀疏甚至根本不稀疏時，SBL算法也表現(xiàn)出很好的優(yōu)勢?？傊?，Bayesian方法對先驗知識的使用使得SBL算法可以得到更好的結果。但是相較于其他算法，SBL算法計算速度偏慢。目前該方法已經(jīng)得到了越來越多的研究和發(fā)展。

5) 置信傳播(Belief Propagation)[10]是一種在圖模型上進行推斷的消息傳遞算法。其主要思想是：對于馬爾可夫隨機場中的每一個節(jié)點，通過消息傳播，把該節(jié)點的概率分布狀態(tài)傳遞給相鄰的節(jié)點，從而影響相鄰節(jié)點的概率分布狀態(tài)，經(jīng)過一定次數(shù)的迭代，每個節(jié)點的概率分布將收斂于一個穩(wěn)態(tài)。在求解欠定矩陣方程時，常見的消息傳遞算法有最小和算法(Min-Sum Algorithm)[3]和近似消息傳遞(Approximate Message Passing，AMP)算法[3]等。最小和算法可以準確計算LASSO估計的高維極限，而AMP算法允許沿著不同大小的實例序列進行漸近精確分析，特別是在大系統(tǒng)中，AMP算法以指數(shù)速度收斂于LASSO最優(yōu)。

欠定矩陣方程稀疏求解算法的研究仍在不斷的發(fā)展和更新中，當為某一個具體的應用選擇算法時，通常需要權衡不同算法的特性，如：計算速度、存儲量、是否易于實現(xiàn)、靈活性、還原效率等因素。本文重點研究以下三個算法：

(3) OMP和LARS算法在遇到高維問題時，性能較差。分步正交匹配算法(StOMP)[12]的設計者將OMP與LARS相結合，并且與迭代收縮算法存在廣泛的相似之處。

3.1 OMP算法

學者已經(jīng)證明了只要采樣矩陣A滿足RIP參數(shù)，OMP算法就可以使用k個步長精確恢復任意k個稀疏信號[5]。

OMP算法的每一步迭代都基于最佳匹配原則，找出稀疏解的一個非零元素，然后用偽逆修正此非零元素，循環(huán)迭代，直到找到所有的非零元素。算法要求已知非零元素的個數(shù)即稀疏度k。OMP算法如算法1所示。

算法1OMP算法

輸入：給定矩陣A∈Rm×n，向量y∈Rm,誤差閾值ε。

輸出：向量x∈Rn。

步驟1初始化：令殘差r=y,標簽集Ω=?。

步驟2辨識：求矩陣A中與殘差向量r最強相關的列，并將其加入標簽集Ω。

Ωk=Ωk-1∪{jk}

其中AΩk=[aω1,aω1,…,aω1],ω1∈Ωk

步驟4更新殘差：rk=y-AΩkxk。

步驟5更新矩陣A，設第jk列為零或刪除此列，因為此列已經(jīng)被選中。

步驟7輸出稀疏向量。

OMP算法時間復雜度為O(kmn)，它的變化形式在精度和時間復雜度都有所改進，而且得到了廣泛的應用。OMP算法在統(tǒng)計建模中被稱為向前逐步回歸(Forward Stepwise Regress,FSW)[1]，常見的改進還有：純匹配追蹤(Pure MP)[1]、松弛匹配追蹤(Relaxed MP)[1]和弱匹配追蹤(Weak MP)[1]等。

3.2 LARS算法

在統(tǒng)計學習中,L1范數(shù)約束問題被稱為LASSO(Least Absolute Shrinkage and Selection Operator)。LASSO由斯坦福大學統(tǒng)計系的Hastie等提出[5]，是一種不等式約束的普通最小二乘法，具有收縮和選擇兩種基本功能，常用于回歸問題。LASSO方法和嶺回歸(Ridge Regression,RR)方法一樣，都有助于降低過擬合的風險，并且LASSO更易獲得“稀疏”解。

最小角回歸(LARS)是一種逐步回歸的同倫算法，是求解LASSO問題[30-32]的有效算法。LARS算法如算法2所示。

算法2LARS算法

輸入：給定矩陣A∈Rm×n，向量y∈Rm。

輸出：向量x∈Rn。

步驟1矩陣A歸一化，使其樣本均值為零，L2范數(shù)為1。

步驟2初始化x=0，支撐集S=supp(x)=?,殘差r0=y-Ax=y。

步驟3尋找與殘差r具有最大相關性的矩陣A的列ai，即λ0=maxi{}。

S=S∪{i}

單變量矩陣As

步驟4fork=1 toK=min(m-1,n)

(2) 以向量ξ為方向，從xk-1開始，朝著它們的最小二乘解移動系數(shù)x，其中As:x(λ)=xk-1+(λk-1-λ)ξ,0≤λ≤λk-1。進一步得到新的殘差r(λ)=y-Ax(λ)。

(3) 關注λ=||,l?S，設第j個變量讓λ達到最大，即λ=||，則將j加入S，有：S=S∪{j}，并令λk=λ。

(4)xk=x(λk)=xk-1+(λk-1-λ)ξ。

(5)rk=y-Axk。

步驟5返回稀疏向量xK。

在LARS算法中，向量集Am×n的角平分線向量(Equiangular vector)μ的定義如下，它與Am×n中的向量的夾角都相等：

μ=Equi(A)=A(ATA)-1In

(21)

式中：In=(1,1,…,1)n；Equi表示求角平分線向量。

LARS保證當前殘差與已經(jīng)入選的變量之間的相關系數(shù)相等，選擇當前殘差在已經(jīng)入選的變量構成的子空間上的投影為求解路徑，繼續(xù)搜索，尋找新的滿足條件的變量，加入之前的子空間，并更新尋找路徑。

LARS算法利用LASSO系數(shù)路徑是分段線性的事實，通常情況下，LARS算法更加準確，計算復雜度更低，并且對于低維問題能夠得到一個更好的解。一般情況下，LARS的步驟等于A的維度，因此在高維問題上，其計算代價更高。

3.3 StOMP算法

與OMP算法相比，StOMP的每個階段都有多個系數(shù)可以進入模型，而OMP的每個階段只能選擇一個系數(shù)。對于每步選擇多個元素的StOMP而言，迭代多次(建議10次)可以得到近似線性系統(tǒng)y≈Ax所需的稀疏解。StOMP算法如算法3所示。

算法3StOMP算法

輸入：給定矩陣A∈Rm×n，向量y∈Rm，誤差閾值ε。

輸出：向量x∈Rn。

步驟1初始化x0=0，支撐集S=supp(x)=?,殘差r0=y-Ax0=y。

步驟2殘差向后投影：

計算ek=AT(y-Axk-1)=ATrk-1

步驟3閾值選擇：

如果選擇閾值使得|Jk|=1，即每次更新一項，則得到OMP。

考慮一般情況，每次迭代可更新多個項。

步驟4支撐集更新：Sk=Sk-1∪Jk。

步驟5帶約束最小二乘：給定支撐集，并假設|Sk|

式中：P∈Rn×|Sk|是一個選擇算子，從向量x中選擇允許為非零項的|Sk|的個數(shù)。

步驟6更新殘差：rk=y-Axk。

步驟7循環(huán)執(zhí)行步驟2到步驟6，直到停止迭代。

步驟8返回稀疏向量x。

4 實 驗

4.1 實驗環(huán)境

本文所述實驗用到的算法程序為自編寫軟件，運行環(huán)境:Windows操作系統(tǒng)；開發(fā)環(huán)境:Visual Studio S2015；算法:C++語言和Armadillo庫。

4.2 光譜混合模型

光譜分析法是根據(jù)物質的光譜來鑒別物質及確定其化學組成和相對含量的方法。在高光譜圖像(Hyperspectral image)[33-41]分析領域，混合光譜常稱為混合像元，純物質稱為“端元”(Endmember)，各個端元所占的比例稱作“豐度”(Abundance)。高光譜解混(或稱混合光譜分解)(Hyperspectral Unmixing,HU)或稱光譜混合分析(Spectral Mixture Analysis)的目的就是分析出混合光譜內包含的不同純物質的光譜信號及其含量。高光譜混合模型主要有線性和非線性兩種。線性模型(Linear Spectral Mixing Model，LSMM)假設單一光子僅能影響一類地物并將反射結果線性累積到像元光譜中。在生物信息學和化學信息學中[42]，如何有效地分析復雜多組分體系中的組分及其含量一直是分析化學研究的難點。自然界的物質往往以混合物的形式存在，比如中草藥成分中的有效和有毒成分分析、痕量物質的檢測等。在生物信息學和化學信息學的光譜分析中存在類似的線性光譜混合模型，根據(jù)多組分體系的朗伯-比爾定律和加和定律，混合物的光譜可以表示為：y=Ax+e，其中:y為目標光譜向量；x為光譜濃度向量；A為已知光譜樣本；e為殘差?；旌瞎庾V解析也就是在已知光譜庫中找到某一光譜或者光譜組合可以最好地逼近目標光譜。

4.3 實驗設置與結果分析

本文描述如下兩個數(shù)據(jù)源不同的實驗：實驗一中的數(shù)據(jù)源光譜矩陣A20×3 683中的各個數(shù)據(jù)之間的距離比較分散，而實驗二中的數(shù)據(jù)源矩陣A91×1 867中包含非常近似的數(shù)據(jù)。因為MDS算法能夠保證原始空間中的距離在低維空間中得以保持，數(shù)據(jù)源均采用多維尺度分析(Multidimensional Scaling,MDS)[43-44]算法進行降維可視化以幫助分析高維數(shù)據(jù)之間的關系。在本實驗中，MDS算法選用光譜分析中常用的Cosine距離作為度量標準。

4.3.1實驗一

量測矩陣(譜庫)A20×3 683由20個長度為3 683的紅外純物質光譜構成，測試用混合光譜分別有帶噪聲和不帶噪聲兩個目標光譜，均通過模擬計算得到，其擬合公式為：

T=0.24S0+0.69S1+0.07S2(+Gaussian)

式中：S0、S1、S2表示光譜；Gaussian表示高斯噪音。

圖3、圖4和圖5分別為量測矩陣(譜庫)A20×3 683的MDS降維可視化、混合物光譜、純物質光譜。

圖4 混合物光譜(0.24S0+0.75S1+0.01S2+Gaussian)

圖5 純物質光譜1(從上到下依次為S0、S1、S2)

OMP、LARS和StOMP三個算法的計算結果如表1所示?？梢钥闯觯惴∣MP、LARS和StOMP在此測試集上均能很好地計算混合物所含純物質及其含量。

表1 混合光譜分解1

4.3.2實驗二

量測矩陣(譜庫)A91×1 867由91個長度為1 867的紅外純物質光譜構成，測試用混合光譜分別有帶噪聲和不帶噪聲兩個目標光譜，均通過模擬計算得到，其擬合公式為：T=0.24S0+0.69S1+0.07S2(+Gaussian)。圖6、圖7和圖8分別為矩陣量測矩陣(譜庫)A91×1 867的MDS降維可視化、混合物光譜純物質光譜，表2為本文算法的計算結果。

圖6 量測矩陣(譜庫)A91×1 867的MDS降維

圖7 混合物光譜(0.24S0+0.75S1+0.01S2+Gaussian)

圖8 純物質光譜2(從上到下依次為S0、S1、S2)

表2 混合光譜分解2

此實驗數(shù)據(jù)表明，算法OMP、LARS和StOMP在此測試集上均存在誤差，因此引入弱匹配追蹤(Weak Matching Pursuit,WMP)[1]和子空間追蹤(Subspace Pursuit,SP)[1]進行對比研究。為了分析算法的計算精度，表3給出了在數(shù)據(jù)源中與S0、S1和S2三種物質的Cosine距離最小的5個光譜。

表3 Cosine距離

當混合光譜不包含噪聲時，LARS和StOMP在此測試集上能很好地計算混合物所含純物質及其含量。但是OMP算法得到的結果有偏差，其中S1的含量最大，得到的結果正確，但是從含量數(shù)據(jù)分析，OMP算法錯誤地把S0和S11混淆，因為這兩個光譜的Cosine距離等于0.998 2；把S2和S83混淆，因為這兩個光譜的Cosine距離等于0.941 0。

當混合光譜含有均值為零、方差為1的高斯噪聲時，對于含量最大的S1，五個算法能得到幾乎都正確的結果，但是WMP的計算精度最高。對于其他混淆的光譜S2與S19，其Cosine距離為0.998 0。

4.3.3結果分析

總體上來說，對于含量最大的光譜S1(含量0.75)，在含噪聲和無噪聲模型中，均能得到比較滿意的結果。對于含量為0.24的光譜S0，OMP算法容易把它與其相近的光譜混淆。但是對于含量只有0.01的光譜S3，因為數(shù)據(jù)源中有與其非常相似的光譜，所以只有LARS和StOMP在不含噪聲的模型中給出了準確的結果，在其他算法和模型中，該物質均沒有被正確識別。

混合物光譜解析的一個難點在于含噪聲系統(tǒng)中精確計算痕量物質的性質及其含量，可能的優(yōu)化方法有：

(1) 譜庫篩選：根據(jù)目標光譜篩選構成矩陣A的譜庫，在其滿足算法成立的數(shù)學條件情況下，盡可能少地包含相似數(shù)據(jù)。比較圖3和圖6，實驗一的數(shù)據(jù)源中的各個光譜的距離較遠(圖3)，而實驗二的數(shù)據(jù)源中有些光譜的距離非常近(圖6)。

(2) 特征選擇：從n個維度中選擇出能代表整體特性的k(k

以實驗一為例，數(shù)據(jù)矩陣A20×3 683的第718號屬性是一個特征峰位置，此特征峰與其他屬性的Pearson相關系數(shù)如表4所示，數(shù)據(jù)顯示，其中有22%的屬性與該屬性的相關系數(shù)大于0.7。

表4 第718號屬性與其他屬相之間的相關系數(shù)

需要說明的是，雖然相關系數(shù)接近于1的程度與數(shù)據(jù)維數(shù)n相關。當n較小時，相關系數(shù)的波動比較大，僅僅根據(jù)相關系數(shù)判定變量x與y之間的線性關系時誤差較大。但是在光譜數(shù)據(jù)中，特征數(shù)n一般都比較大，簡單起見，本文采用相關系數(shù)作為度量標準來說明特征之間的線性相關性。

(3) 濾波：采用合適的濾波技術和算法，盡可能地提高混合光譜的信噪比。

(4) 恢復算法的選擇和優(yōu)化：根據(jù)具體問題選擇合適的算法，提高計算精度和性能。

本文未分析各個算法的時間復雜度，特別是面對大數(shù)據(jù)量時，時間復雜度會是一個衡量算法性能的關鍵因素。除了發(fā)展速度更快的算法外，并行計算、GPU、通過小波變化壓縮數(shù)據(jù)等都是提高算法計算速度的有效方法。

5 結 語

壓縮感知自創(chuàng)始以來受到了學術界和工業(yè)界的廣泛關注，它利用了信號的稀疏特性，通過求解欠定線性系統(tǒng)的稀疏解來有效地恢復重建信號，改變了傳統(tǒng)信號處理中數(shù)據(jù)采集然后壓縮(特征提取)再傳輸(處理)的模式，在信息理論、信號/圖像處理、醫(yī)學成像、模式識別、生物信息、光學/雷達成像和無線通信等諸多領域有著重要的理論研究和應用價值。

壓縮感知相關的理論比較復雜，本文從欠定線性矩陣方程開始，介紹了壓縮感知的基礎知識。最后，通過兩種不同光譜庫的具體數(shù)值實驗，重點討論了OMP、LARS和StOMP稀疏恢復算法在混合光譜解析時的性能和存在的問題，并給出相應的優(yōu)化建議。通常認為求解壓縮感知問題的方法有貪婪追蹤、凸松弛方法、迭代收縮算法、貝葉斯框架、消息傳遞/置信傳播(Belief Propagation)算法等。通過壓縮感知算法精確恢復信號的方法基于光譜庫完備的假設，隨著壓縮感知和稀疏表示理論的蓬勃發(fā)展，學者們對稀疏性有了更為深刻的認識，稀疏化建模在光譜學等一些高維數(shù)據(jù)分析領域中都有著廣闊的應用前景和深入的研究空間。