楊　浩，柳炳利，張二喜，郭　科

（成都理工大學(xué) 數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室，四川 成都　610059）

基于奇異值分解的巖心高光譜數(shù)據(jù)降噪研究

楊浩，柳炳利，張二喜，郭科

（成都理工大學(xué) 數(shù)學(xué)地質(zhì)四川省重點(diǎn)實(shí)驗(yàn)室，四川 成都610059）

為了更有效地消除巖心高光譜數(shù)據(jù)中的噪聲，提出基于奇異值分解的巖心高光譜數(shù)據(jù)降噪方法，引入奇異值下降率的概念，利用奇異值下降率單調(diào)性的突變點(diǎn)來確定表征信號(hào)有用奇異值的個(gè)數(shù)。用該方法對(duì)地物光譜儀ASD Field?Spec?4實(shí)地采集到的巖心高光譜數(shù)據(jù)進(jìn)行降噪處理，并與依據(jù)奇異值相對(duì)強(qiáng)度確定奇異值突變點(diǎn)的降噪方法進(jìn)行對(duì)比，利用均方根誤差（RMSE）、信噪比（SNR）兩項(xiàng)指標(biāo)對(duì)降噪效果進(jìn)行評(píng)價(jià)。實(shí)驗(yàn)結(jié)果表明，該方法更能提高信噪比，降低均方根誤差，更能有效保持原始巖心高光譜曲線的吸收特征，消除高光譜曲線上的毛噪現(xiàn)象。

奇異值分解；奇異值下降率；巖心高光譜數(shù)據(jù)；降噪

0　引　言

高光譜數(shù)據(jù)蘊(yùn)含了豐富的巖石、礦物的反射光譜信息，但在實(shí)際野外采集光譜的過程中，周圍環(huán)境因素和一系列的人為因素都會(huì)使采集到的光譜數(shù)據(jù)含有大量噪聲，例如：光照條件的變化，周圍環(huán)境濕度、溫度等自然條件。這嚴(yán)重影響了地物反射光譜中的吸收特征，大大降低了數(shù)據(jù)的分析精度。因此，對(duì)高光譜數(shù)據(jù)的降噪研究是極其必要的。

國(guó)內(nèi)外不少學(xué)者對(duì)高光譜數(shù)據(jù)降噪進(jìn)行了研究。Boardman和Kruse提出了用最小噪聲分離技術(shù)（MNF）隔離高光譜數(shù)據(jù)中的噪聲［1］；徐冬等提出了一種基于多元線性回歸的高光譜遙感數(shù)據(jù)去噪方法［2］；印佳等使用主成分分析進(jìn)行高光譜數(shù)據(jù)去噪［3］；陳志剛等提出了基于經(jīng)驗(yàn)?zāi)B(tài)分解高光譜圖像光譜域去噪方法［4］；周丹等應(yīng)用小波閾值對(duì)高光譜進(jìn)行小波分析去除光譜噪聲［5］；路威等實(shí)現(xiàn)了高光譜遙感數(shù)據(jù)的三次光滑樣條去噪［6］。

以上方法大都基于光譜曲線的平滑處理，容易造成光譜特征丟失，而小波分析處理高光譜噪聲一般需要估計(jì)噪聲的統(tǒng)計(jì)特征來確定小波閾值，對(duì)噪聲水平有一定依賴性。為了尋求能夠適應(yīng)高光譜降噪要求的降噪方法，本文利用Hankel矩陣和SVD（奇異值分解）對(duì)野外實(shí)地采集到的巖心高光譜數(shù)據(jù)進(jìn)行降噪研究，引入了奇異值下降率的概念，以確定奇異值個(gè)數(shù)。

1　Hankel矩陣和SVD

Hankel矩陣是指每一條副對(duì)角線上的元素都相等的方陣。設(shè)方陣 A=［ai，j］m×n∈Cm×n，如果 ai，j=ai-1，j+1i-1，則稱矩陣 A為Hankel矩陣。SVD（Singular Value De?composition）即矩陣的奇異值分解，是現(xiàn)代數(shù)值分析最基本和最重要的工具之一。SVD在最優(yōu)化問題、統(tǒng)計(jì)學(xué)、信號(hào)處理以及工程技術(shù)等方面都有重要作用。

設(shè) A∈Cm×n，且r=rank（A），因?yàn)?AHA是半正定矩陣，且rank（AHA）=rank（A）=r，所以AHA的特征值可以排列 為 ：λ1≥λ2≥…≥λr＞0=λr+1=λr+2=…=λn則 稱為A的奇異值，λi為 AHA的特征值。進(jìn)一步，存在m階酉矩陣U和n階酉矩陣V，使得：

成立。其中Σ=diag（σ1，σ2，…，σr），而σi（i=1，2，…，r）稱為矩陣A的正奇異值，此時(shí)式（1）稱為矩陣A的奇異值分解式［7］。將式（1）寫成求和形式有：

式中：ui為AAH的第i個(gè)單位特征向量；vi為AHA的第i個(gè)單位特征向量。

2　SVD降噪原理

基于SVD分解的濾波方法是一種非線性濾波，它從矩陣的角度出發(fā)，將包含信號(hào)特征的矩陣分解到一系列奇異值和奇異值矢量對(duì)應(yīng)的子空間中［8］，本文將實(shí)測(cè)巖心高光譜數(shù)據(jù)構(gòu)造成Hankel矩陣，然后對(duì)Hankel矩陣進(jìn)行奇異值分解，進(jìn)行降噪處理，最后重構(gòu)巖心高光譜信號(hào)，得到降噪后的巖心高光譜信號(hào)。

野外實(shí)測(cè)巖心高光譜反射率序列x為：

式（4）中：s（i）代表去噪后的巖心高光譜反射率信號(hào)；n（i）代表噪聲信號(hào)，i=1，2，…，N，N為光譜反射率數(shù)據(jù)序列長(zhǎng)度。由H（i，j）=x（i+j-1）構(gòu)造m×n階Hankel矩陣［9］：

式中［10］，N=m+n-1。

將式（4）代入式（5）有：

式中：Hs為去噪后的巖心高光譜反射率信號(hào) s（i）的m×n階Hankel矩陣Hs（i，j）=s（i+j-1）；Hn為噪聲n（i）的m×n階Hankel矩陣 Hn（i，j）=n（i+j-1）；對(duì) H進(jìn)行奇異值分解，有：

參照式（2），將式（7）寫成：

式（7）、式（8）中，酉矩陣U∈Cm×m，酉矩陣V∈Cn×n，且U=［u1，u2，…，um］，V=［v1，v2，…，vn］，σi（i=1，2，…，r）為矩陣H的正奇異值，r=min（m，n），Σ=diag（σ1，σ2，…，σr），ui為 HHT的第i個(gè)單位特征向量，vi為 HTH的第i個(gè)單位特征向量。選取前L個(gè)奇異值及其對(duì)應(yīng)的特征向量重構(gòu)H的逼近矩陣：

經(jīng)過SVD分解并依據(jù)信號(hào)和噪聲各自的特點(diǎn)，即真實(shí)信號(hào)s（i）與噪聲信號(hào)n（i）之間的不相關(guān)性［11］，以及真實(shí)信號(hào)能量比較集中，而噪聲信號(hào)能量比較分散的特點(diǎn)，可以將由含噪的測(cè)量信號(hào)所構(gòu)成的Hankel矩陣分解成兩個(gè)互不相關(guān)的空間：真實(shí)信號(hào)空間和噪聲空間，從Hankel矩陣H中去除噪聲空間，就能夠得到經(jīng)過濾波后的信號(hào)空間，進(jìn)而得到降噪后的信號(hào)。一般而言，不再是一個(gè)Hankel矩陣，必須構(gòu)造一個(gè)Hankel矩陣如式（9）所示，來逼近，按照式（9）來構(gòu)造，對(duì)矩陣Hs的反對(duì)角線求平均值即可得到真實(shí)信號(hào)在每一時(shí)刻的值［9］s（i），i=1，2，…，N。同理，矩陣的反對(duì)角線求平均值，從而得到經(jīng)過降噪后，真實(shí)信號(hào)在每一時(shí)刻的估計(jì)，即SVD降噪后的信號(hào)。

3　奇異值下降率

SVD降噪的關(guān)鍵在于奇異值個(gè)數(shù)L的選取。Hanli Qiao提出以下規(guī)則［12］確定L：將H的奇異值從大到小排列，依次為σi（i=1，2，…，r），選取前L個(gè)奇異值之和與所有奇異值之和之比剛好超過90%的L作為SVD降噪的奇異值個(gè)數(shù)。張波使用奇異值相對(duì)強(qiáng)度曲線［8］的突變點(diǎn)確定降噪奇異值個(gè)數(shù)L，趙學(xué)智提出奇異值差分譜理論［13］，bi=σi-σi+1，i=1，2，…，r-1，以差分譜bi的第二個(gè)最大峰值所在坐標(biāo)i作為降噪奇異值個(gè)數(shù)L。

第一種確定奇異值個(gè)數(shù)的方法，參數(shù)90%的確定太主觀；第二種方法對(duì)于奇異值相對(duì)強(qiáng)度曲線突變點(diǎn)的界定僅靠人為觀察，會(huì)造成較大誤差；第三種方法當(dāng)兩相鄰奇異值差別不大時(shí)，方法失效，不能確定奇異值個(gè)數(shù)。

本文從奇異值突變角度出發(fā)，引入奇異值下降率的概念，以便科學(xué)地確定SVD降噪奇異值個(gè)數(shù)。設(shè)奇異值從大到小排序后的序列為（σ1，σ2，…，σr），為描述奇異值序列的突變特性，定義：

為第i+1個(gè)奇異值的奇異值下降率。由于奇異值是按降序排列，它描述了奇異值序列的下降快慢程度，前面的值較大的奇異值描述的是信號(hào)的信息。ci大表示第i+1個(gè)奇異值較第i個(gè)奇異值下降得快，ci小表示第i+1個(gè)奇異值較第i個(gè)奇異值下降得慢，ci單調(diào)遞減時(shí)奇異值表征了信號(hào)，當(dāng)ci在i點(diǎn)單調(diào)性發(fā)生變化，即此位置前后奇異值在性質(zhì)上出現(xiàn)了突變，說明i點(diǎn)以后的奇異值表征了噪聲的信息。取前i個(gè)奇異值進(jìn)行SVD降噪處理。

4　巖心高光譜數(shù)據(jù)SVD降噪處理

本次研究的巖心高光譜數(shù)據(jù)是利用美國(guó)ASD公司生產(chǎn)的地物光譜儀ASD FieldSpec?4在湖北省雞冠咀銅金礦巖心庫(kù)實(shí)地采集的。ASD FieldSpec?4是為針對(duì)地物環(huán)境遙感研究而專門設(shè)計(jì)的，能夠有效地獲取可見和近紅外光譜（Visible and Near?Infrared，VNIR），短波紅外光譜（Short?Wave Infrared，SWIR）。該儀器采集光譜范圍：350～2 500 nm。實(shí)測(cè)含噪巖心高光譜如圖1所示。由于ASD地物光譜儀采用了三個(gè)探測(cè)器：350～1 000 nm是512陣元PDA探測(cè)，1 000～1 800 nm 及1 800～2 500 nm兩個(gè)獨(dú)立的InGaAs探測(cè)器。導(dǎo)致采集到的巖心高光譜數(shù)據(jù)在波長(zhǎng)1 000 nm和1 800 nm處的反射率出現(xiàn)了階躍變化。需要對(duì)實(shí)測(cè)光譜進(jìn)行階躍處理，消除階躍，消除階躍后的巖心高光譜曲線如圖2所示。

由于Hankel矩的秩r=min（m，n）；r太大會(huì)導(dǎo)致Hankel矩陣非零奇異值個(gè)數(shù)過多，難以用少數(shù)幾個(gè)奇異值表征信號(hào)，r太小，會(huì)使噪聲和信號(hào)難以分離；取Han?kel矩陣的維數(shù)為2 080×72，由圖2中的巖心高光譜反射率數(shù)據(jù)構(gòu)造Hankel矩陣G，對(duì)矩陣G進(jìn)行奇異值分解。奇異值分布曲線如圖3所示；奇異值相對(duì)強(qiáng)度曲線如圖4所示；奇異值差分譜如圖5所示；奇異值下降率曲線如圖6所示。

圖1　實(shí)測(cè)含噪巖心高光譜曲線

圖2　處理階躍后含噪巖心高光譜曲線

圖3　奇異值分布曲線

圖4　奇異值相對(duì)強(qiáng)度曲線

圖5　奇異值差分譜曲線

圖6　奇異值下降率曲線

對(duì)于Hnakel矩陣G的奇異值的突變點(diǎn)的界定，圖3和圖4的方法都指示第2個(gè)奇異值為突變點(diǎn)，而圖5差分譜并沒有出現(xiàn)波峰波谷，方法失效。取前2個(gè)奇異值進(jìn)行SVD降噪，降噪后光譜曲線如圖7所示。本文提出奇異值下降率曲線確定奇異值個(gè)數(shù)，奇異值下降率曲線如圖6所示，第9個(gè)點(diǎn)以前的奇異值下降率程依次遞減關(guān)系，在第10個(gè)點(diǎn)處的奇異值下降率大于第9個(gè)點(diǎn)的奇異值下降率，說明奇異值發(fā)生了突變，有噪聲摻雜到信號(hào)中了，取突變以前的9個(gè)奇異值進(jìn)行SVD降噪處理。降噪后光譜曲線如圖8所示。以兩種方式確定突變點(diǎn)的SVD降噪方法去噪指標(biāo)對(duì)比如表1所示。

圖7　奇異值相對(duì)強(qiáng)度確定突變點(diǎn)降噪后光譜曲線

由表1可以看出，與以奇異值相對(duì)強(qiáng)度確定奇異值突變點(diǎn)的SVD降噪方法相比，本文提出的以奇異值下降率確定SVD突變點(diǎn)的降噪方法具有高信噪比，低均方根誤差的特點(diǎn)。

圖8　奇異值下降率確定突變點(diǎn)降噪后光譜曲線

表1　兩種方法確定突變點(diǎn)SVD降噪指標(biāo)對(duì)比

對(duì)比圖2，圖7，由奇異值相對(duì)強(qiáng)度確定奇異值突變點(diǎn)的SVD降噪后巖心高光譜曲線在波長(zhǎng)約1 850 nm左右丟失了吸收峰，圖2中降噪前巖心高光譜曲線在波長(zhǎng)2 250 nm以后的位置隱約有吸收峰存在，但在圖7中，2 250 nm以后的吸收特征被去除了，降噪效果存在很大問題。

對(duì)比圖2，圖8，本文提出的奇異值下降率確定奇異值突變點(diǎn)的SVD降噪后光譜曲線保持了降噪前光譜曲線在1 850 nm波長(zhǎng)處的吸收特征，并去除了2 250 nm以后的巖心高光譜曲線毛噪現(xiàn)象，凸顯了2 300 nm以后的2個(gè)吸收峰位。研究表明［14］綠泥石化礦物主要的吸收峰位置為2 360 nm，碳酸鹽化礦物的主要吸收峰為2 342 nm，本文的降噪方法準(zhǔn)確揭示了2 300 nm以后的巖心高光譜數(shù)據(jù)的吸收峰，為后續(xù)巖心高光譜礦化蝕變信息的提取奠定了基礎(chǔ)。

5　結(jié)　論

本文提出了基于奇異值下降率確定奇異值突變點(diǎn)的SVD巖心高光譜數(shù)據(jù)降噪方法，并應(yīng)用此方法對(duì)巖心高光譜數(shù)據(jù)進(jìn)行降噪研究，研究結(jié)果表明，奇異值差分譜方法不能確定巖心高光譜數(shù)據(jù)降噪的SVD奇異值個(gè)數(shù)；奇異值相對(duì)強(qiáng)度確定突變點(diǎn)的巖心高光譜數(shù)據(jù)SVD降噪光譜曲線較原始高光譜曲線，丟失很多吸收峰；奇異值下降率能自動(dòng)識(shí)別巖心高光譜數(shù)據(jù)SVD降噪過程中的奇異值突變點(diǎn)，無需經(jīng)驗(yàn)判斷，解決了SVD降噪的關(guān)鍵問題；基于奇異值下降率的SVD巖心高光譜數(shù)據(jù)降噪方法有效去除了巖心高光譜數(shù)據(jù)中噪聲，消除了巖心高光譜曲線的毛噪現(xiàn)象，準(zhǔn)確提取了2 300 nm以后的巖心高光譜數(shù)據(jù)的吸收峰，為后續(xù)巖心高光譜礦化蝕變信息的提取奠定了基礎(chǔ)。

注：本文通訊作者為柳炳利。

［1］BOARDMAN J W，KRUSE F A.Automated spectral analysis：a geological example using AVIRIS data［C］//Proceedings of the Tenth Thematic Conference on Geologic Remote Sensing. North Grapevine Mountains，Nevada：［s.n.］，1994：407?418.

［2］徐冬，孫蕾，羅建書.基于多元線性回歸的高光譜遙感圖像小波去噪［J］.遙感信息，2013（6）：78?81.

［3］印佳，杜戰(zhàn)戰(zhàn).基于主成分分析的高光譜遙感圖像非局部去噪［J］.現(xiàn)代電子技術(shù)，2015，38（11）：70?72.

［4］陳志剛，束炯.高光譜圖像光譜域噪聲去除的經(jīng)驗(yàn)?zāi)B(tài)分解方法［J］.紅外與毫米波學(xué)報(bào)，2008，27（5）：378?382.

［5］周丹，王欽軍，田慶久，等.小波分析及其在高光譜噪聲去除中的應(yīng)用［J］.光譜學(xué)與光譜分析，2009，29（7）：1941?1945.

［6］路威，余旭初，劉娟.高光譜遙感數(shù)據(jù)三次光滑樣條濾噪［J］.測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2005，22（1）：11?13.

［7］吳昌愨，魏洪增.矩陣?yán)碚撆c方法［M］.北京：電子工業(yè)出版社，2013：79?80.

［8］張波，李健君.基于Hankel矩陣與奇異值分解（SVD）的濾波方法以及在飛機(jī)顫振試驗(yàn)數(shù)據(jù)預(yù)處理中的應(yīng)用［J］.振動(dòng)與沖擊，2009，28（2）：162?166.

［9］DOCLO S，DOLOGLOU I，MOONEN M.A novel iterative sig?nal enhancement algorithm for noise reduction in speech［R/ OL］.［1999?07?18］.http：//www.researchgate.net/publication/ 28264.

［10］GOLAFSHAN Reza，SANLITURK Kenan Yuce.SVD and Hankel matrix based de?noising approach for ball bearing fault detection and its assessment using artificial faults［J］. Mechanical systems and signal processing，2016（70/71）：36?50.

［11］高建，張飛艷，謝偉，等.利用能量變分方法進(jìn)行高光譜數(shù)據(jù)去噪處理［J］.武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2012，37（3）：322?325.

［12］QIAO Hanli.New SVD based initialization strategy for non?negative matrix factorization［J］.Pattern recognition letters，2015，63：71?77.

［13］趙學(xué)智，葉邦彥，陳統(tǒng)堅(jiān).奇異值差分譜理論及其在車床主軸箱故障診斷中的應(yīng)用［J］.機(jī)械工程學(xué)報(bào)，2010，46（1）：100?108.

［14］張媛，張杰林，趙學(xué)勝，等.基于峰值權(quán)重的巖心高光譜礦化蝕變信息提取［J］.國(guó)土資源遙感，2015（2）：154?159.

A new method about core hyperspectral data denoising based on singular value decomposition

YANG Hao，LIU Bingli，ZHANG Erxi，GUO Ke
（Sichuan Province Key Laboratory for Mathematical Geology，Chengdu University of Technology，Chengdu 610059，China）

In order to eliminate the noise in the core hyperspectral data more effectively，a new method about core hyper?spectral data denoising based on singular value decomposition is proposed，in which the concept of singular value decline rate is brought.The number of useful singular value of the characterization signal is determined by the abrupt change point of the singu?lar value decline rate.This method is used to denoise the core hyperspectral data collected by ASD FieldSpec?4，and compared with the denoising method that relies on the singular value relative strength to determine singular value mutation point.Its denois?ing effect was evaluated with the root mean square error（RMSE）and the signal?to?noise ratio（SNR）.The results show that the method can improve SNR，reduce RMSE，keep the absorption characteristics of original core hyperspectral curve more effective?ly，and eliminate the frizz phenomenon of core hyperspectral curve.

singular value decomposition；singular value decrease rate；core hyperspectral data；denoising

TN957.54?34；TN911.74

A

1004?373X（2016）18?0004?05

10.16652/j.issn.1004?373x.2016.18.002

楊浩（1990—），男，彝族，四川喜德人，碩士研究生。主要從事高光譜數(shù)據(jù)處理的研究工作。柳炳利（1981—），男，講師，博士。主要從事數(shù)學(xué)地質(zhì)、遙感地質(zhì)方面的研究工作。

2016?01?27

國(guó)家自然科學(xué)基金資助項(xiàng)目（41272363）；中國(guó)地質(zhì)調(diào)查局地質(zhì)調(diào)查項(xiàng)目：基于巖心高光譜的原生暈地球化學(xué)部找礦預(yù)測(cè)研究（12120114002001）