基于稀疏表示分步重構(gòu)算法的高光譜目標(biāo)檢測

廖佳俊，劉志剛，姜江軍，路志勇

?

基于稀疏表示分步重構(gòu)算法的高光譜目標(biāo)檢測

廖佳俊1，劉志剛1，姜江軍1，路志勇2

（1. 火箭軍工程大學(xué)，陜西 西安 710025；2. 96831 部隊(duì)，北京 100015）

針對(duì)傳統(tǒng)稀疏表示重構(gòu)算法在高光譜目標(biāo)檢測中表現(xiàn)出運(yùn)算速度慢的問題，提出了分步重構(gòu)算法（Two Steps Reconstruction，TSR）。該方法先求得個(gè)與待測像元最相似的字典原子，然后用這些原子線性表示待測像元以求解稀疏向量，舍棄了傳統(tǒng)重構(gòu)算法的迭代求解的方式，直接通過求解逆矩陣，簡化了運(yùn)算過程，使運(yùn)算速度大幅提高。本文給出了方法的具體過程并將其與傳統(tǒng)方法及其改進(jìn)方法進(jìn)行比較。實(shí)驗(yàn)結(jié)果表明，TSR在保證檢測精度不下降的同時(shí)能夠大幅提升運(yùn)算速度。

高光譜目標(biāo)檢測；重構(gòu)算法；稀疏表示

0 引言

目標(biāo)檢測是高光譜遙感技術(shù)應(yīng)用的重要方向之一，它涵蓋了環(huán)境監(jiān)測、城市調(diào)查、礦物填圖和軍事偵察等諸多領(lǐng)域。傳統(tǒng)的目標(biāo)檢測方法有光譜角填圖（Spectral Angle Mapping，SAM）、匹配子空間（Matched Subspace Detector，MSD）、正交子空間投影（Orthogonal Subspace Projection，OSP）等[1]。這些方法在高光譜目標(biāo)檢測方面有著廣泛的運(yùn)用，但也存在檢測結(jié)果精度不高、對(duì)分布模型假設(shè)準(zhǔn)確與否關(guān)系大等不足之處。今年來，稀疏表示理論在圖像處理及目標(biāo)檢測中取得了較好的發(fā)展[2-3]。2010年，Yi Chen等人將稀疏表示成功應(yīng)用于高光譜圖像分類和目標(biāo)檢測[4-5]，取得了不錯(cuò)的效果。越來越多的學(xué)者開始關(guān)注將稀疏表示應(yīng)用于高光譜圖像，稀疏表示為高光譜目標(biāo)檢測提供了一條新的研究思路?；谙∈璞硎镜姆椒m然能克服傳統(tǒng)方法精度不高以及需要假設(shè)分布模型等問題，但由于高光譜圖像光譜分辨率高帶來的海量數(shù)據(jù)，使得目前稀疏表示中應(yīng)用比較廣泛的正交匹配追蹤算法（Orthogonal matching pursuit，OMP）[6]在計(jì)算時(shí)由于迭代次數(shù)過多而造成運(yùn)算量大、速度慢的問題，特別是當(dāng)字典中原子數(shù)目過多或者稀疏度過高時(shí)會(huì)使得計(jì)算復(fù)雜度成指數(shù)增長。在壓縮感知領(lǐng)域，不少學(xué)者提出了許多OMP算法的改進(jìn)形式，包括gOMP、StOMP、SWOMP和CoSaMP[7-10]等；趙春暉等人將StOMP算法應(yīng)用到高光譜目標(biāo)檢測之中，取得了一定的效果[11]。我們將上述改進(jìn)方法應(yīng)用在高光譜目標(biāo)檢測中發(fā)現(xiàn)，由于光譜的復(fù)雜性使得高光譜待測像元并不能由字典中的原子精確重構(gòu)，有些算法在運(yùn)行時(shí)可能并不是由重構(gòu)誤差達(dá)到精度而完成重構(gòu)而是由于達(dá)到迭代次數(shù)而完成重構(gòu)，所以目標(biāo)檢測效果并不理想。針對(duì)這個(gè)問題，本文提出了一種分步重構(gòu)算法，放棄了原來迭代求解的思路，而通過直接求解逆矩陣實(shí)現(xiàn)重構(gòu)，使運(yùn)算速度得到大幅提高。

1 基于稀疏表示的高光譜目標(biāo)檢測及OMP算法概述

稀疏表示的思想就是將原始信號(hào)表示為字典與稀疏向量相乘的形式：

＝(1)

式中：是一個(gè)過完備字典；稀疏向量是一個(gè)只有少數(shù)元素不為零的向量。

在高光譜目標(biāo)檢測中，一個(gè)待測像元可以通過大量端元中部分原子的線性組合表示。對(duì)應(yīng)于稀疏表示，可將大量端元構(gòu)成一個(gè)超完備字典，由目標(biāo)子字典t和背景子字典b級(jí)聯(lián)構(gòu)成，那么待測像元可以由下式表示：

式中：t和b分別對(duì)應(yīng)目標(biāo)子字典t和背景子字典b的稀疏向量。通過求解稀疏向量并判斷t和b是零向量還是稀疏向量即可判斷該待測像元是目標(biāo)還是背景。算法基本原理如圖1所示。

該方法關(guān)鍵性步驟之一就是求解稀疏向量的重構(gòu)算法，其中最為常見的就是OMP。其主要思路可以簡單總結(jié)如下：

Step 1：初始化殘差0＝，索引矩陣0＝?，字典支撐集0＝?，迭代次數(shù)＝1。

Step 4：更新殘差r＝－。

2 分步重構(gòu)算法

重構(gòu)算法的目標(biāo)是找到個(gè)合適的字典原子以表示待測像元。OMP算法采用逐個(gè)求解的思路，由于反復(fù)迭代，運(yùn)算時(shí)間很長。為了提高效率，可否一次性求解個(gè)合適的字典原子呢？根據(jù)這一思路，本文提出了一種分步重構(gòu)的算法，舍棄了迭代求解的方式，試圖使運(yùn)算過程得以簡化。該方法流程圖如圖2所示，具體步驟如下：

第一步：從字典中選取個(gè)與待測像元最相似的原子。

在這一步中本文提出兩種求解方法：線性表示法和直接表示法。相對(duì)而言，線性表示法在像素點(diǎn)多，稀疏度高的情況下速度更快；直接表示法在稀疏度稍低的情況下速度更快。其方法具體原理如下：

1）線性表示法（Two Steps Reconstruct with Linear Representation，TSR-L）

假設(shè)待測像元可以表示為如下等式：

＝11＋…＋da＝(3)

式中：表示待測像元；1, …,d表示字典中的各個(gè)原子；1, …,a表示系數(shù)。

由于通常是奇異矩陣，可以通過式(4)求解：

＝(D＋I)－1x(4)

式中：是一個(gè)極小的正則化項(xiàng)；是單位矩陣。

通過式(3)可以發(fā)現(xiàn)每一個(gè)原子都對(duì)待測像元的表示都做出了貢獻(xiàn)。計(jì)算每一個(gè)原子與待測像元的殘差：

其中ei可以表示每一個(gè)字典原子與待測像元之間的差別，較小的ei對(duì)應(yīng)字典原子對(duì)待測像元貢獻(xiàn)度越大。那么通過比較ei就能夠確定K個(gè)與待測像元最相似的字典原子，記為d1,…,dK，該方法通過逆矩陣一次求出每個(gè)原子對(duì)應(yīng)的貢獻(xiàn)度，對(duì)于稀疏度越高的情況效果越明顯。

Fig.1 Basic principle of Sparse representation

圖2 算法流程圖

Fig.2 Flow chart of algorithm

2）直接表示法（Two Steps Reconstruct with Direct Representation TSR-D）

第一步主要問題是求出個(gè)與待測像元相似的字典原子，由于待測信號(hào)與字典原子都進(jìn)行了歸一化處理，假設(shè)字典原子與待測像元的相似度可由下式表示：

L＝xd,＝1, …,(6)

式中：L越大表示相似度越高，通過排序后，取出前個(gè)值，即可求得與待測像元最相似的個(gè)字典原子，記為1,…,d，該方法由于計(jì)算簡單，在大多數(shù)情況下效果突出。如果待測像元是目標(biāo)像元，那么1,…,d中全部或者大部分都應(yīng)該屬于目標(biāo)子字典，否則應(yīng)全部或大部分屬于背景子字典。

第二步：用選出的個(gè)字典原子的線性組合進(jìn)一步表示待測像元，計(jì)算稀疏向量。

假設(shè)待測像元滿足如下等式：

＝11＋…＋bd＝bD(7)

式中：1,…,d表示與待測像元最相似的字典原子；1,…,b表示系數(shù)。由于通常也是奇異矩陣，那么可以通過式(7)求解：

式中：是一個(gè)極小的正則化項(xiàng)；是單位矩陣。

然后，將中原子和根據(jù)式(8)求出來的對(duì)應(yīng)的系數(shù)的乘積之和分別求目標(biāo)和背景的重構(gòu)誤差t和b，如式(9)所示：

式中：t表示目標(biāo)子字典原子及其對(duì)應(yīng)系數(shù)相乘之和；b表示背景子字典原子及其對(duì)應(yīng)系數(shù)相乘之和；若待測像元為目標(biāo)，那么t應(yīng)當(dāng)越小，b應(yīng)當(dāng)越大，反之亦然。最后設(shè)置一個(gè)閾值，通過求t和b的差值()，如式(10)所示：

根據(jù)上述步驟，即可實(shí)現(xiàn)對(duì)高光譜圖像的目標(biāo)檢測處理?？梢钥闯觯痉椒ㄉ釛壛说蠼庀∈柘蛄康男问剑ㄟ^直接求逆矩陣的形式求解稀疏向量，簡化了運(yùn)算。

3 實(shí)驗(yàn)結(jié)果及分析

在本章中，采用3組實(shí)測高光譜圖像數(shù)據(jù)對(duì)本文提出的算法進(jìn)行驗(yàn)證，以接收機(jī)特性（Receiver Operating Characteristic，ROC）曲線下的面積（Area Under the Curve，AUC）[12]為檢測結(jié)果的效果指標(biāo)，將檢測效果和時(shí)間與OMP、StOMP、gOMP算法進(jìn)行對(duì)比。本文字典的選取采用文獻(xiàn)[3]提出的方法，其中目標(biāo)子字典原子個(gè)數(shù)依次分別為52，38，38；背景子字典窗口大小為內(nèi)窗7×7，外窗11×11，共計(jì)72個(gè)原子；稀疏度＝16。實(shí)驗(yàn)環(huán)境為Dell Precision T1650圖形工作站；8核處理器，主頻3.30GHz；內(nèi)存16G；操作系統(tǒng)：Microsoft Windows7；程序運(yùn)行平臺(tái)：MATLAB 2014a。

3.1 MSHyperSIS水杯數(shù)據(jù)

第一組實(shí)驗(yàn)數(shù)據(jù)是用地面MSHyperSIS成像光譜儀獲取的一幅高光譜圖像。圖像原始大小為260×336，共256個(gè)波段。本實(shí)驗(yàn)中所用的圖像數(shù)據(jù)是從原始圖像中裁取出來的，其大小為80×80，剔除其中的無效波段，保留有效波段181個(gè)。實(shí)驗(yàn)數(shù)據(jù)如圖3所示，圖像區(qū)域?yàn)橐黄嗟?，水泥地中放置了一個(gè)水杯為待檢測的目標(biāo)。其中，圖3(a)是第69個(gè)波段的灰度圖，圖3(b）是水杯的空間分布圖。

圖3 實(shí)驗(yàn)數(shù)據(jù)圖像

Fig.3 Experimental image

圖4(a)～(e)給出了OMP、StOMP、gOMP、TSR-L、TSR-D的目標(biāo)檢測效果圖。其目標(biāo)檢測所對(duì)應(yīng)的AUC值以及所消耗的時(shí)間如表1所示。

圖4 圖3的檢測結(jié)果對(duì)比

表1 圖3用不同重構(gòu)方法所得AUC值及時(shí)間

從圖4和表1中可以看出各種算法檢測結(jié)果相差不大，OMP算法運(yùn)算速度相對(duì)較慢；基于OMP算法改進(jìn)的StOMP算法速度有所提升，但檢測精度略有下降；gOMP算法運(yùn)算速度不升反降；本文提出的分步重構(gòu)算法對(duì)比于OMP算法在檢測精度略有提高的情況下，TSR-D運(yùn)算速度大幅提升，降低到OMP的15.2%，而TSR-L運(yùn)算時(shí)間則略有增加。

3.2 AVIRIS數(shù)據(jù)

第2組數(shù)據(jù)為機(jī)載可見光及紅外成像光譜儀（Airborne Visible/Infrared Imaging Spectrometer, AVIRIS）拍攝的美國圣地亞哥海軍機(jī)場的高光譜圖像。該高光譜圖像波長范圍為370～2500nm，空間分辨率為3.5m，大小為400×400個(gè)像元。截取其中100×100的子圖像，除去水汽吸收波段和信噪比較低的波段后，保留189個(gè)波段。圖5(a)是第204波段的灰度圖，圖5(b)是停靠飛機(jī)的空間分布圖。

圖6(a)～(e)給出了OMP、StOMP、gOMP、TSR-L、TSR-D的目標(biāo)檢測效果圖。其目標(biāo)檢測所對(duì)應(yīng)的ROC曲線下面積值以及所消耗的時(shí)間如表2所示。

圖5 實(shí)驗(yàn)數(shù)據(jù)圖像

從圖6和表2中可以看出OMP算法運(yùn)算速度相對(duì)較慢；基于OMP算法改進(jìn)的StOMP算法速度有所提升，檢測精度略有下降；gOMP算法運(yùn)算速度不升反降；本文提出的分步重構(gòu)算法對(duì)比于OMP算法在檢測精度略有提高的情況下，運(yùn)算速度大幅提升，其中TSR-L運(yùn)算時(shí)間降低到OMP的64.2%，TSR-D降低到OMP的29.3%。

圖6 圖5的檢測結(jié)果對(duì)比

表2 圖5用不同重構(gòu)方法所得AUC值及時(shí)間

3.3 MSHyperSIS迷彩服數(shù)據(jù)

第3組實(shí)驗(yàn)數(shù)據(jù)是用地面MSHyperSIS成像光譜儀獲取的一片以樹林為背景的高光譜圖像。圖像原始大小為500×336，共256個(gè)波段。本實(shí)驗(yàn)中所用的圖像數(shù)據(jù)是從原始圖像中裁取出來的，其大小為200×200，剔除其中的無效波段，保留有效波段181個(gè)。實(shí)驗(yàn)數(shù)據(jù)如圖7所示，圖像區(qū)域?yàn)橐黄瑯淞?，樹林中放置?件迷彩服為待檢測的目標(biāo)。其中，圖7(a)是第69波段的灰度圖，圖7(b)是迷彩服的空間分布圖。

圖7 實(shí)驗(yàn)數(shù)據(jù)圖像

圖8(a)～(e)給出了OMP、StOMP、gOMP、TSR-L、TSR-D的目標(biāo)檢測效果圖。其目標(biāo)檢測所對(duì)應(yīng)的ROC曲線下面積值以及所消耗的時(shí)間如表3所示。

圖8 圖7的檢測結(jié)果對(duì)比

表3 圖7用不同重構(gòu)方法所得AUC值及時(shí)間

從圖8和表3中可以看出OMP算法運(yùn)算速度也相對(duì)較慢；基于OMP算法改進(jìn)的StOMP算法速度有所提升，檢測精度略有下降；gOMP算法運(yùn)算速度不升反降；本文提出的分步重構(gòu)算法對(duì)比于OMP算法在檢測精度略有提高的情況下，運(yùn)算速度大幅提升，其中TSR-L運(yùn)算時(shí)間降低到OMP的62.5%，TSR-D降低OMP的28.9%。

3.4 實(shí)驗(yàn)結(jié)果分析

對(duì)比上述3組實(shí)驗(yàn)結(jié)果可以看到，StOMP算法相比于OMP算法在運(yùn)算速度上有一定的提升，但精度有少許降低，這可能是因?yàn)槊看蔚x擇多個(gè)原子不如選擇單個(gè)原子精確度高。gOMP算法在運(yùn)算速度上不升反降，這可能是因?yàn)檫\(yùn)行時(shí)不是由重構(gòu)誤差達(dá)到精度而完成重構(gòu)而是由于達(dá)到迭代次數(shù)而完成重構(gòu)。TSR-D在所有數(shù)據(jù)中均有最好的表現(xiàn)，不僅精度高而且運(yùn)算速度比傳統(tǒng)的OMP算法提升超過70%；TSR-L在像素點(diǎn)比較多的情況下能保持比較高的精度，同時(shí)速度也比OMP和gOMP算法快，但相比TSR-D速度上還有一定差距。本文進(jìn)一步研究發(fā)現(xiàn)，當(dāng)稀疏度越高時(shí)，TSR-L相比TSR-D有更快的速度，如圖9所示，這可能是因?yàn)門SR-L在第一步時(shí)也是通過逆矩陣求解，而相比之下求解逆矩陣對(duì)稀疏度的增加更不敏感。進(jìn)一步對(duì)比本文算法以及OMP類算法的原理，我們發(fā)現(xiàn)在進(jìn)行字典原子的選擇時(shí)，由于原理不同各自選出的原子可能并不相同，本文算法選出的原子本身和待測像元更相似，所以該方法是合理的。綜上，本文提出的兩種方法在檢測精度，特別是運(yùn)算速度上相比傳統(tǒng)算法有明顯提高，并且在稀疏度高低不同的情況下有各自的優(yōu)勢。

圖9 不同稀疏度下TSR-D和TSR-L運(yùn)行時(shí)間對(duì)比圖

4 結(jié)束語

本文針對(duì)傳統(tǒng)稀疏表示重構(gòu)算法在高光譜目標(biāo)檢測中表現(xiàn)出的運(yùn)算速度較慢的問題，提出了分步重構(gòu)算法。該算法舍棄了傳統(tǒng)方法使用迭代求解稀疏向量的方法，通過直接求解逆矩陣使得運(yùn)算簡化，從而達(dá)到了大幅提升算法運(yùn)算速度的效果，并且兩種方法在稀疏度高低不同的情況下有各自的優(yōu)勢。本文提出的方法簡單易行，與傳統(tǒng)算法及其改進(jìn)算法進(jìn)行目標(biāo)檢測實(shí)驗(yàn)對(duì)比的結(jié)果表明，該算法在保證檢測精度不下降的同時(shí)大幅提升了運(yùn)算速度。

[1] 張兵, 高連如. 高光譜圖像分類與目標(biāo)探測[M]. 北京: 科學(xué)出版社, 2011: 223, 248.

ZHANG B, GAO L R.[M]. Beijing:,2011: 223, 248.

[2] 楊春偉, 王仕成, 廖守億,等. 基于核稀疏編碼的紅外目標(biāo)識(shí)別方法[J]. 紅外技術(shù), 2016, 38(3):230-235.

YANG C W, WANG S C, LIAO S Y.An infrared target recognition method based on Kernel Sparse coding[J]., 2016, 38(3): 230-235.

[3] 孫君頂, 趙慧慧. 圖像稀疏表示及其在圖像處理中的應(yīng)用[J]. 紅外技術(shù), 2014, 36(7): 533-537.

SUN J X, ZHAO H H.Sparse representation and applications in image processing[J]., 2014, 36(7):533-537.

[4] CHEN Yi, Nasrabadi N M, Tran T D. Sparsity-based classification of hyperspectral imagery[C]//2010(IGARSS), 2010: 2796-2799.

[5] CHEN Yi, Nasrabadi N M, Tran T D. Sparse representation for target detection in hyperspectral imagery[J]., 2011, 5(3): 629-640.

[6] Pati Y C, Rezaiifar R, Krishnaprasad P S. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition[C]//, 1993, 1: 1-3.

[7] Donoho D L, Tsaig Y, Drori I, et al. Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit[J]., 2012, 58(2): 1094-1121.

[8] WANG J, Kwon S, Shim B. Generalized Orthogonal Matching Pursuit[J]., 2011, 60(12): 6202-6216.

[9] Blumensath, T, Davies Mike E. Stagewise weak gradient pursuits[J]., 2009, 57(11): 4333-4346.

[10] Needell D, Tropp J A. CoSaMP: Iterative signal recovery from incomplete and inaccurate samples[J]., 2009, 26(3): 301-321.

[11] 趙春暉, 靖曉昊, 李威. 基于StOMP稀疏方法的高光譜圖像目標(biāo)檢測[J]. 哈爾濱工程大學(xué)學(xué)報(bào), 2015, 36(7): 992-996.

ZHAO C, JING X, LI W. Hyperspectral imagery target detection algorithm based on StOMP sparse representation[J]., 2015, 36(7): 992-996.

[12] WANG Y, HUANG S, LIU D, et al. A novel band selection method based on curve area and genetic theory[J]., 2014, 43(3): 193-202.

Target Detection in Hyperspectral Image Using Two Steps Reconstruction Based on Sparse Representation

LIAO Jiajun1，LIU Zhigang1，JIANG Jiangjun1，LU Zhiyong2

(1.,710025,; 2.96831,,100015,)

Aiming at the efficiency problem of traditional reconstruction algorithm for target detection in hyperspectral image, a two steps reconstruction algorithm was proposed.nearest atoms with text pixel were calculated for representing the text pixel to calculate sparse vectors. This method solves the problem by calculating inverse matrix instead of iteration, which simplifies the process. The specific process of the method is given and the method is compared with the traditional and improved ones. Experimental results show that this method can improve the computation speed without decreasing the accuracy of detection.

hyperspectral target detection，reconstruction algorithm，sparse representation

TP75

A

1001-8891(2016)08-0699-06

2016-04-08；

2016-05-17.

廖佳?。?992-），男，湖南湘潭人，碩士研究生，主要從事遙感圖像處理研究。

國家自然科學(xué)基金項(xiàng)目（41574008）。