伍 娟 妮

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

0 引 言

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

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

y=Axx∈Rn×1

(1)

式中：A∈Rm×n(m

(2)

s.t.y=As

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

(3)

s.t.y=As

相關(guān)理論證明了當(dāng)且僅當(dāng)矩陣A滿足受限零空間性質(zhì)時(shí)，L0范數(shù)問(wèn)題與L1范數(shù)問(wèn)題等價(jià)[1]。

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

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

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

1 基礎(chǔ)概念

1.1 向量范數(shù)

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

(4)

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

(5)

(6)

(7)

(8)

1.2 線性模型

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

Ax=y

(9)

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

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

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

(3)m

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

矩陣方程的高效求解算法一直是數(shù)學(xué)和計(jì)算機(jī)算法領(lǐng)域研究的一個(gè)核心，本文概述上述矩陣方程的常見(jiàn)求解方法并重點(diǎn)說(shuō)明欠定矩陣方程的適用條件和常用解法。

1.2.1超定方程m>n

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

針對(duì)矩陣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逆矩陣。

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

(13)

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

xLS=(ATA+λI)-1ATy

(14)

在信號(hào)處理領(lǐng)域該方法稱為松弛法。在機(jī)器學(xué)習(xí)和統(tǒng)計(jì)學(xué)中,稱該L2約束問(wèn)題為嶺回歸。

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

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

1.2.2盲矩陣方程

考慮盲矩陣方程:

X=AS

(15)

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

1.2.3欠定矩陣方程m

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

2 壓縮感知

2.1 概 述

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

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

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

(16)

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

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

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

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

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

(17)

則稱矩陣A滿足k限定等距性(kRIP)。此時(shí)可以通過(guò)下面的優(yōu)化問(wèn)題近似完美地從y中恢復(fù)稀疏信號(hào)s，進(jìn)而恢復(fù)出信號(hào)x：

(18)

s.t.y=As

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

(19)

s.t.y=As

當(dāng)且僅當(dāng)矩陣Am×n滿足受限零空間性質(zhì)時(shí)，L0與L1問(wèn)題等價(jià)[5]。

2.2 數(shù)學(xué)模型

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

y=Axx∈Rn×1

(20)

式中：A∈Rm×n(m

圖1 CS模型

圖2 Lp模型

(1) 當(dāng)0

(2) 當(dāng)p=1時(shí)，Lp球是菱狀球，當(dāng)球與直線Ax=y相交于坐標(biāo)軸上時(shí)得到一個(gè)稀疏解。

(3) 當(dāng)p>1時(shí)，Lp球是外凸球，直線Ax=y與外凸球切點(diǎn)一定不在坐標(biāo)軸上，因而，此時(shí)的解是不稀疏解。

3 算 法

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

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

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

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

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

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

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

3.1 OMP算法

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

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

算法1OMP算法

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

輸出：向量x∈Rn。

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

步驟2辨識(shí)：求矩陣A中與殘差向量r最強(qiáng)相關(guān)的列，并將其加入標(biāo)簽集Ω。

Ωk=Ωk-1∪{jk}

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

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

步驟5更新矩陣A，設(shè)第jk列為零或刪除此列，因?yàn)榇肆幸呀?jīng)被選中。

步驟7輸出稀疏向量。

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

3.2 LARS算法

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

最小角回歸(LARS)是一種逐步回歸的同倫算法，是求解LASSO問(wèn)題[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具有最大相關(guān)性的矩陣A的列ai，即λ0=maxi{}。

S=S∪{i}

單變量矩陣As

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

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

(3) 關(guān)注λ=||,l?S，設(shè)第j個(gè)變量讓?duì)诉_(dá)到最大，即λ=||，則將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保證當(dāng)前殘差與已經(jīng)入選的變量之間的相關(guān)系數(shù)相等，選擇當(dāng)前殘差在已經(jīng)入選的變量構(gòu)成的子空間上的投影為求解路徑，繼續(xù)搜索，尋找新的滿足條件的變量，加入之前的子空間，并更新尋找路徑。

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

3.3 StOMP算法

與OMP算法相比，StOMP的每個(gè)階段都有多個(gè)系數(shù)可以進(jìn)入模型，而OMP的每個(gè)階段只能選擇一個(gè)系數(shù)。對(duì)于每步選擇多個(gè)元素的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殘差向后投影：

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

步驟3閾值選擇：

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

考慮一般情況，每次迭代可更新多個(gè)項(xiàng)。

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

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

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

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

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

步驟8返回稀疏向量x。

4 實(shí) 驗(yàn)

4.1 實(shí)驗(yàn)環(huán)境

本文所述實(shí)驗(yàn)用到的算法程序?yàn)樽跃帉?xiě)軟件，運(yùn)行環(huán)境:Windows操作系統(tǒng)；開(kāi)發(fā)環(huán)境:Visual Studio S2015；算法:C++語(yǔ)言和Armadillo庫(kù)。

4.2 光譜混合模型

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

4.3 實(shí)驗(yàn)設(shè)置與結(jié)果分析

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

4.3.1實(shí)驗(yàn)一

量測(cè)矩陣(譜庫(kù))A20×3 683由20個(gè)長(zhǎng)度為3 683的紅外純物質(zhì)光譜構(gòu)成，測(cè)試用混合光譜分別有帶噪聲和不帶噪聲兩個(gè)目標(biāo)光譜，均通過(guò)模擬計(jì)算得到，其擬合公式為：

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

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

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

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

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

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

表1 混合光譜分解1

4.3.2實(shí)驗(yàn)二

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

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

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

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

表2 混合光譜分解2

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

表3 Cosine距離

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

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

4.3.3結(jié)果分析

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

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

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

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

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

表4 第718號(hào)屬性與其他屬相之間的相關(guān)系數(shù)

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

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

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

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

5 結(jié) 語(yǔ)

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

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