本發(fā)明涉及一種高光譜異常目標(biāo)檢測(cè)方法，特別是涉及一種基于結(jié)構(gòu)稀疏表示和內(nèi)部聚類濾波的高光譜異常目標(biāo)檢測(cè)方法。

背景技術(shù)：

高光譜異常目標(biāo)檢測(cè)技術(shù)是一種不需要提供待測(cè)目標(biāo)先驗(yàn)光譜信息的前提下，僅利用圖像像元之間的光譜差異性進(jìn)行目標(biāo)檢測(cè)的高光譜目標(biāo)檢測(cè)技術(shù)，在實(shí)際應(yīng)用中具有較強(qiáng)的實(shí)用性。

相對(duì)于整個(gè)高光譜圖像，異常目標(biāo)不僅出現(xiàn)的概率較低，而且所占的比例也較少，傳統(tǒng)的異常目標(biāo)檢測(cè)算法一般假設(shè)圖像背景服從高斯分布。在這種假設(shè)下可以利用全局或局部的統(tǒng)計(jì)特性來檢測(cè)異常目標(biāo)。然而，在實(shí)際應(yīng)用中由于受到空間分辨率限制，圖像背景具有較高的復(fù)雜性，僅通過假設(shè)背景服從高斯分布來進(jìn)行異常目標(biāo)檢測(cè)，不能取得很好的檢測(cè)效果。

近年來，基于稀疏表示的高光譜異常目標(biāo)檢測(cè)算法受到越來越多的關(guān)注。該類算法假設(shè)圖像中的像元可以通過背景和目標(biāo)的字典線性稀疏表示，在求得稀疏向量后，利用重建誤差的大小提取異常目標(biāo)。在實(shí)際應(yīng)用中具有很好的靈活性和可擴(kuò)展性。文獻(xiàn)“l(fā)ij,zhangh,zhangl,etal.hyperspectralanomalydetectionbytheuseofbackgroundjointsparserepresentation[j].selectedtopicsinappliedearthobservationsandremotesensing,ieeejournalof,2015,8(6):2523-2533.”公開了一種基于背景聯(lián)合稀疏表示的高光譜異常目標(biāo)檢測(cè)方法。該方法首先以全部圖像像元作為字典，通過稀疏表示模型求解稀疏向量，根據(jù)所求的非零稀疏向量位置確定局部背景像元，然后根據(jù)局部背景像元建立局部正交子空間提取異常目標(biāo)。文獻(xiàn)所述方法在選取局部背景像元時(shí)每次都使用了大規(guī)模的全局像元作為字典，給求解稀疏向量帶來較大的計(jì)算負(fù)擔(dān)；在求解稀疏向量時(shí)沒有考慮稀疏向量?jī)?nèi)部的結(jié)構(gòu)稀疏信息，造成求解稀疏向量不準(zhǔn)確；未有效考慮高光譜圖像中像元之間的空間約束，目標(biāo)檢測(cè)效率不高。

技術(shù)實(shí)現(xiàn)要素：

為了克服現(xiàn)有高光譜異常目標(biāo)檢測(cè)方法目標(biāo)檢測(cè)效率低的不足，本發(fā)明提供一種基于結(jié)構(gòu)稀疏表示和內(nèi)部聚類濾波的高光譜異常目標(biāo)檢測(cè)方法。該方法在稀疏表示理論框架下，首先基于局部rx算法選擇初始背景像元后，利用主成分分析字典學(xué)習(xí)法學(xué)習(xí)得到魯棒性的背景字典。稀疏向量求解和圖像重建過程中，引入重加權(quán)拉普拉斯先驗(yàn)，提高稀疏向量求解精度。計(jì)算原始圖像與重建圖像之間的誤差得到稀疏表示誤差。此外，利用內(nèi)部聚類濾波表示高光譜數(shù)據(jù)的空譜特性，通過計(jì)算待測(cè)像元與它同類其他像元線性表示結(jié)果的誤差，得到內(nèi)部聚類誤差，最后聯(lián)合稀疏表示誤差與內(nèi)部聚類誤差的線性加權(quán)和實(shí)現(xiàn)異常目標(biāo)的精確提取。在真實(shí)的高光譜衛(wèi)星圖像aviris、hydice和仿真的高光譜數(shù)據(jù)集上的試驗(yàn)結(jié)果表明，本發(fā)明獲得的檢測(cè)結(jié)果相對(duì)于已有的異常目標(biāo)檢測(cè)方法在恒虛警率的前提下檢測(cè)率提高了10％～15％。

本發(fā)明解決其技術(shù)問題所采用的技術(shù)方案是：一種基于結(jié)構(gòu)稀疏表示和內(nèi)部聚類濾波的高光譜異常目標(biāo)檢測(cè)方法，其特點(diǎn)是包括以下步驟：

步驟一、輸入一幅高光譜圖像，選用局部雙窗rx算法獲取背景像元集合。局部rx算法計(jì)算公式如下：

dlocalrx(xi)＝(xi-μlocal(i))^t(∑local(i))^-1(xi-μlocal(i))(1)

式中，xi表示第i個(gè)像元，μlocal(i)表示xi像元所對(duì)應(yīng)的局部背景的均值，(∑local(i))^-1表示xi像元所對(duì)應(yīng)的局部背景的協(xié)方差矩陣的逆矩陣，利用公式(1)遍歷所有高光譜像元后選取閾值進(jìn)行背景分割，小于閾值的對(duì)應(yīng)像元為背景像元，得到背景像元矩陣zlocal，大小為nb×nz，nb表示高光譜數(shù)據(jù)波段個(gè)數(shù)，nz表示局部背景像元的個(gè)數(shù)。所述閾值取0.007～0.01；

步驟二、根據(jù)主成分分析學(xué)習(xí)算法，對(duì)步驟(一)得到的背景像元矩陣zlocal轉(zhuǎn)置后求其協(xié)方差矩陣，定義為covzlocal，大小為nb×nb。然后求協(xié)方差矩陣covzlocal的特征值v和特征向量p，將特征值按照從大到小排列后，相應(yīng)的特征向量也依次排列，得到的新的特征向量矩陣即為最終學(xué)習(xí)到的背景字典d。

步驟三、在稀疏表示框架下，高光譜圖像像元能被學(xué)習(xí)到的字典d稀疏表示，即

x＝dy+n(2)

式中，x為高光譜圖像，d為步驟(二)學(xué)習(xí)到的背景字典，y＝[y1,…,yn]為稀疏表示矩陣，n為稀疏表示誤差和圖像噪聲。將重加權(quán)拉普拉斯稀疏先驗(yàn)引入y，得到基于重加權(quán)拉普拉斯稀疏先驗(yàn)的結(jié)構(gòu)稀疏表示優(yōu)化模型。

步驟四、采用坐標(biāo)下降法求解優(yōu)化模型得到最優(yōu)稀疏信號(hào)yrec，然后通過公式得到重建后的高光譜圖像xrec。根據(jù)原圖像與重建圖像得到重建誤差，重建誤差計(jì)算公式如下：

rsprase(x)＝||x-xrec||2(3)

其中，x、xrec分別表示原圖像和重建后圖像，rsprase(x)表示基于稀疏表示的重建誤差，||·||2表示向量的2范數(shù)。

步驟五、采用k-means++算法對(duì)輸入的高光譜數(shù)據(jù)進(jìn)行聚類。假設(shè)輸入圖像x被聚為m個(gè)類，表示第m個(gè)聚類包含的像元集合，nb代表圖像波段數(shù)量，nm代表第m個(gè)聚類有n個(gè)像元。

步驟六、計(jì)算類內(nèi)像元之間的相似權(quán)值，計(jì)算公式如下：

式中，表示第m個(gè)聚類里第i、j個(gè)像元，h是尺度因子，為歸一化因子。

步驟七、針對(duì)每一個(gè)圖像像元，用與之同類型的其他像元在相似權(quán)值的基礎(chǔ)之上進(jìn)行線性表示，表示公式如下：

式中，表示第m個(gè)聚類中第i個(gè)元素被同一聚類中其他像元線性加權(quán)表示后的結(jié)果。

步驟八、計(jì)算原來像元與線性表示后的像元之間的表示誤差。如果該像元是背景像元，則能夠被該類其他背景像元線性表示，表示誤差較??；否則，若該像元為異常像元，則表示誤差較大。據(jù)此能夠提取異常目標(biāo)。表示誤差的計(jì)算公式如下：

式中，x、分別是原圖像和聚類表示后的圖像，rcluster(x)表示基于內(nèi)部聚類的重建誤差，||·||2表示向量的2范數(shù)。

步驟九、加權(quán)求和稀疏重建誤差和聚類誤差，最終得到異常目標(biāo)檢測(cè)結(jié)果。計(jì)算公式如下：

rfinal＝0.5×rsparse+0.5×rcluster(7)

在兩個(gè)重建誤差加權(quán)求和所得的灰度圖基礎(chǔ)之上，選取分割閾值δ，圖像中大于δ的像素標(biāo)記為1表示目標(biāo)，小于等于δ的像素標(biāo)記0表示背景，δ的取值范圍為0.1～0.3。最終得到只有0和1的二值結(jié)果圖，完成異常目標(biāo)檢測(cè)。

本發(fā)明的有益效果是：該方法在稀疏表示理論框架下，首先基于局部rx算法選擇初始背景像元后，利用主成分分析字典學(xué)習(xí)法學(xué)習(xí)得到魯棒性的背景字典。稀疏向量求解和圖像重建過程中，引入重加權(quán)拉普拉斯先驗(yàn)，提高稀疏向量求解精度。計(jì)算原始圖像與重建圖像之間的誤差得到稀疏表示誤差。此外，利用內(nèi)部聚類濾波表示高光譜數(shù)據(jù)的空譜特性，通過計(jì)算待測(cè)像元與它同類其他像元線性表示結(jié)果的誤差，得到內(nèi)部聚類誤差，最后聯(lián)合稀疏表示誤差與內(nèi)部聚類誤差的線性加權(quán)和實(shí)現(xiàn)異常目標(biāo)的精確提取。在真實(shí)的高光譜衛(wèi)星圖像aviris、hydice和仿真的高光譜數(shù)據(jù)集上的試驗(yàn)結(jié)果表明，本發(fā)明獲得的檢測(cè)結(jié)果相對(duì)于已有的異常目標(biāo)檢測(cè)方法在恒虛警率的前提下檢測(cè)率提高了10％～15％。

下面結(jié)合具體實(shí)施方式對(duì)本發(fā)明作詳細(xì)說明。

具體實(shí)施方式

本發(fā)明基于結(jié)構(gòu)稀疏表示和內(nèi)部聚類濾波的高光譜異常目標(biāo)檢測(cè)方法具體步驟如下：

假設(shè)輸入的高光譜圖像是一個(gè)三維數(shù)據(jù)立方體，包含nb個(gè)波段，每個(gè)波段是一幅nrow行和ncol列大小的圖像。為了計(jì)算方便，將每個(gè)波段拉伸成為一個(gè)行向量，所有行向量組成一個(gè)二維矩陣x，其中，x的每一列表示每一個(gè)像素對(duì)應(yīng)的光譜，該方向?yàn)楣庾V維；x的每一行對(duì)應(yīng)一個(gè)波段的所有像素值(即np＝nrow×ncol)，該方向?yàn)榭臻g維。具體如下：

1、基于主成分分析的魯棒性背景字典學(xué)習(xí)。

(1)采用雙窗局部rx算法獲取背景像素集合。

根據(jù)輸入高光譜圖像分辨率選擇窗口尺寸都為奇數(shù)的矩形外窗outer和內(nèi)窗inner(outer＝nouter×nouter,inner＝ninner×ninner,nouter>ninnner,nouter、ninner分別表示外窗和內(nèi)窗的尺寸，都為奇數(shù))。以每一個(gè)輸入像素xi∈x,i＝1,…,np為中心，根據(jù)外窗和內(nèi)窗尺寸獲取局部背景區(qū)域nlocal＝outer-inner，nlocal為nb×nlocal大小的局部背景矩陣，然后計(jì)算nlocal的均值和協(xié)方差矩陣分別得到μlocal＝[μ1,…,μn]^t和∑local。局部rx計(jì)算公式如下：

dlocalrx(xi)＝(xi-μlocal(i))^t(∑local(i))^-1(xi-μlocal(i))(1)

式中，μlocal(i)表示第i個(gè)輸入像元所對(duì)應(yīng)的局部背景的均值，(∑local(i))^-1表示第i輸入像元所對(duì)應(yīng)的局部背景的協(xié)方差矩陣的逆矩陣，利用公式(1)遍歷所有高光譜像元后選取閾值進(jìn)行背景分割，小于閾值的對(duì)應(yīng)像元為背景像元，得到背景像元矩陣zlocal，大小為nb×nz，nz表示局部背景像元的個(gè)數(shù)。為了得到純凈的背景像元，選擇較小的閾值較為合適，建議取值為0.007～0.01。

(2)基于主成分分析的魯棒性背景字典學(xué)習(xí)。

將背景像元矩陣zlocal轉(zhuǎn)置后求其協(xié)方差矩陣，定義為covzlocal，大小為nb×nb。然后求協(xié)方差矩陣covzlocal的特征值v和特征向量p(每一個(gè)特征值對(duì)應(yīng)一個(gè)特征向量)，將特征值按照從大到小排列后，相應(yīng)的特征向量也依次排列，得到的新的特征向量矩陣即為最終學(xué)習(xí)到的背景字典d。

2、建立基于重加權(quán)拉普拉斯稀疏先驗(yàn)的結(jié)構(gòu)稀疏表示模型。

在稀疏表示框架下，高光譜圖像像元可以被學(xué)習(xí)到的字典d稀疏表示，即

x＝dy+n(2)

d為已學(xué)習(xí)到的背景字典，y＝[y1,…,yn]為稀疏表示矩陣，n為稀疏表示誤差和圖像噪聲。

假設(shè)n服從的矩陣正態(tài)分布，則對(duì)應(yīng)的似然函數(shù)為

其中，σn＝diag(λ)是以λ元素為對(duì)角線元素的對(duì)角矩陣，用于表示誤差和噪聲的強(qiáng)度。表示q矩陣的加權(quán)跡范數(shù)。

為了表示稀疏向量?jī)?nèi)部的結(jié)構(gòu)稀疏性，本發(fā)明將重加權(quán)拉普拉斯稀疏先驗(yàn)引入y。首先假設(shè)y服從如下分布

其中，σy＝diag(γ)表示以γ的元素為對(duì)角線元素的對(duì)角矩陣，控制y中每一行的稀疏度，γi＝0表示y的第i行為0。假設(shè)則其中任一列yi服從高斯分布。假設(shè)超參數(shù)γ服從以下的伽馬分布，

以上兩級(jí)先驗(yàn)，等價(jià)于重加權(quán)拉普拉斯分布，因?yàn)閷?duì)于yi有

其中，為保證求解便利，本發(fā)明中不直接使用重加權(quán)拉普拉斯先驗(yàn)，而是采用級(jí)聯(lián)先驗(yàn)，而且λ，γ和κ均為待估計(jì)參數(shù)。

由于λ，γ和κ未知，無法采用最大后驗(yàn)估計(jì)(maximumaposteriorestimation，map)直接對(duì)稀疏信號(hào)y進(jìn)行估計(jì)。因此，本發(fā)明根據(jù)經(jīng)驗(yàn)貝葉斯框架，先基于輸入的高光譜數(shù)據(jù)x利用map估計(jì)未知參數(shù)λ，γ和κ，如下

其中，λopt,γopt,κopt分別表示最優(yōu)的λ，γ，κ。通過積分，并引入-2log運(yùn)算，容易得知式子(7)等價(jià)于最小化如下的式子

其中，tr(·)表示跡范數(shù)，σby＝σn+dσyd^t，為代價(jià)函數(shù)。

通過變形式子(8)，得到稀疏信號(hào)y的非分離稀疏約束模型。首先，對(duì)式子(8)的第一份部分進(jìn)行變形

然后，將式子(9)帶入到式子(8)中，得

接著，引入新的代價(jià)方程如下

顯然，而且可以證明，最小化式子(8)再對(duì)稀疏信號(hào)y進(jìn)行map估計(jì)，與直接最小化式子(11)得到的λ，γ和κ相同，關(guān)于y的解僅相差一個(gè)常量因此，式子(11)可以看作是關(guān)于稀疏信號(hào)y的正則化回歸模型，其中為稀疏信號(hào)的非分離稀疏約束。該約束不能拆分成對(duì)于y中每一行的獨(dú)立約束，因此該約束能同時(shí)約束稀疏信號(hào)中非零元素，潛在地考量非零元素之間的相關(guān)性。此外，σby中包含了表征噪聲強(qiáng)度的λ，因此，得到的稀疏約束可以隨著估計(jì)的噪聲強(qiáng)度自適應(yīng)的變化，具有噪聲魯棒性。

最終得到了如下的基于重加權(quán)拉普拉斯稀疏先驗(yàn)的結(jié)構(gòu)稀疏表示優(yōu)化模型：

3、模型求解并重建高光譜圖像。

已知待檢測(cè)高光譜數(shù)據(jù)x，本發(fā)明采用坐標(biāo)下降法求解式子(12)，每次迭代中僅優(yōu)化一個(gè)變量而固定剩余的其他變量。λ⁰，γ⁰，κ⁰分別表示初始的數(shù)值，t記錄迭代次數(shù)，η表示更新閾值，λ^t，γ^t，κ^t，y^t分別表示迭代到第t次時(shí)的數(shù)值，λ^t+1，γ^t+1，κ^t+1，y^t+1分別表示迭代到第t+1次時(shí)的數(shù)值。

具體步驟如下：

①初始化，λ⁰，γ⁰，κ⁰均初始化為對(duì)應(yīng)長(zhǎng)度的全1向量，計(jì)數(shù)變量t＝0；

②更新中間變量σn＝diag(λ^t)，σy＝diag(γ^t)，σby＝σn+dσyd^t；

③固定λ^t，γ^t和κ^t，根據(jù)式子(11)得到關(guān)于y的優(yōu)化形式，如下

求解得到y(tǒng)的更新規(guī)則如下，

④固定y^t+1，λ^t和κ^t，得到關(guān)于γ的優(yōu)化形式，如下

求解得到如下的更新形式：

其中，為γ^t+1的第i個(gè)元素，代表v^t+y^t+1(y^t+1)^t的對(duì)角線元素組成的向量，zi為z的第i個(gè)元素；

⑤固定y^t+1，γ^t+1和κ^t，得到關(guān)于λ的優(yōu)化形式，如下

求解得到如下的更新形式：

其中，根號(hào)運(yùn)算表示向量每一個(gè)元素開方后組成的向量，./運(yùn)算代表兩個(gè)向量對(duì)應(yīng)元素相除后組成的向量，代表對(duì)角線元素組成的向量。

⑥固定y^t+1，γ^t+1和λ^t+1，得到關(guān)于κ的優(yōu)化形式，如下:

求解得到如下的更新形式：

上式子中的加法和除法運(yùn)算均作用在向量的每一個(gè)元素上，得到一個(gè)新向量，d＝10^-6的引入是為了確保γ中出現(xiàn)0時(shí)，式子(20)依然有意義。

⑦計(jì)算稀疏信號(hào)y更新前后的差異，如下

其中，表示對(duì)y^t+1內(nèi)的每一個(gè)元素乘以||·||f表示弗羅貝尼烏斯范數(shù)(frobeniusnorm),如果計(jì)數(shù)器t>5或者更新差異η<10^-3，則退出循環(huán)；否則t+1，循環(huán)執(zhí)行②至⑦。

⑧假設(shè)上述循環(huán)結(jié)束得到的最優(yōu)稀疏信號(hào)為yrec，則待重建的高光譜圖像xrec可以通過如下方式得到：

4、計(jì)算重建誤差。

顯然xrec是基于背景字典重建的結(jié)果，由此可知原高光譜圖像中背景像素可以被精確重建，而異常像素不能被精確重建。重建誤差計(jì)算公式如下：

rsprase(x)＝||x-xrec||2(23)

其中，x、xrec分別表示原圖像和重建后圖像，都為nb×np二維矩陣，||·||2表示向量的2范數(shù)。

5、基于內(nèi)部聚類濾波算法的高光譜圖像聚類。

(1)利用k-means++算法對(duì)圖像聚類。

針對(duì)輸入高光譜圖像x，設(shè)置聚類個(gè)數(shù)k的值(根據(jù)不同圖像數(shù)值不同，通常范圍為30～50)，然后進(jìn)行以下具體步驟：

①隨機(jī)選擇一個(gè)像元作為第一個(gè)聚類中心；

②對(duì)于圖像中的每一個(gè)像元x，計(jì)算它與最近聚類中心(指已選擇的聚類中心)的距離d(x)；

③選擇一個(gè)新的像元作為新的聚類中心，選擇原則是d(x)較大的點(diǎn)，被選取作為下一個(gè)聚類中心；

④重復(fù)步驟②和③直到k個(gè)聚類中心被選出來；

⑤遍歷每個(gè)像元計(jì)算它與k個(gè)聚類中心的距離，該像元與哪個(gè)聚類中心最近就將它劃分到該類，這樣將所有圖像像元聚集為k類；

⑥分別計(jì)算由⑤得到的k個(gè)聚類的質(zhì)心，得到k個(gè)新的聚類中心；

⑦重復(fù)步驟⑤和⑥，直到質(zhì)心不再變化，得到將原圖像x聚為k個(gè)類的結(jié)果。假設(shè)輸入圖像x被聚為m個(gè)類表示第m個(gè)聚類包含的像元集合，nb代表圖像波段數(shù)量，nm代表第m個(gè)聚類有n個(gè)像元。；

(2)計(jì)算類內(nèi)像元之間的相似權(quán)值，計(jì)算公式如下：

式中表示第m個(gè)聚類里第i、j個(gè)像元，h是尺度因子，為歸一化因子。

(3)針對(duì)每一個(gè)圖像像元，用與之同類型的其他像元在相似權(quán)值的基。

礎(chǔ)之上進(jìn)行線性表示，表示公式如下：

式中，表示第m個(gè)聚類中第i個(gè)元素被同一聚類中其他像元線性加權(quán)表示后的結(jié)果。

(4)計(jì)算原來像元與線性表示后的像元之間的表示誤差。

如果該像元是背景像元，則可以被該類其他背景像元線性表示，表示誤差較??；否則，若該像元為異常像元，則表示誤差較大。據(jù)此可以提取異常目標(biāo)。表示誤差的計(jì)算公式如下：

式中，x、分別是原圖像和聚類表示后的圖像。

6、加權(quán)求和稀疏重建誤差和聚類誤差，最終得到異常目標(biāo)檢測(cè)結(jié)果。通過上述計(jì)算，分別得到了結(jié)構(gòu)稀疏表示和內(nèi)部聚類誤差，通過加權(quán)

計(jì)算兩者誤差之和，最終分割出異常目標(biāo)，得到檢測(cè)結(jié)果。計(jì)算公式如下：

rfinal＝0.5×rsparse+0.5×rcluster(27)

根據(jù)灰度圖rfinal選取分割閾值δ，圖像中大于δ的像素標(biāo)記為1表示目標(biāo)，小于等于δ的像素標(biāo)記0表示背景，δ的建議取值范圍為0.1～0.3。最終得到只有0和1的二值結(jié)果圖，完成異常目標(biāo)檢測(cè)，結(jié)束。