王燦進(jìn)，孫 濤，王挺峰，郭 勁，劉玉龍

（1.中國(guó)科學(xué)院長(zhǎng)春光學(xué)精密機(jī)械與物理研究所 激光與物質(zhì)相互作用國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)春130033；2.吉林省煙草專賣局 信息中心，長(zhǎng)春130033）

0 引 言

激光主動(dòng)成像［1］系統(tǒng)具有分辨率高、抗干擾能力強(qiáng)、不受環(huán)境光限制等優(yōu)點(diǎn)，適用于微光、夜視以及遠(yuǎn)距離暗目標(biāo)的探測(cè)成像，目前已廣泛應(yīng)用于遠(yuǎn)距離監(jiān)視、自動(dòng)目標(biāo)識(shí)別等領(lǐng)域［2－3］。在主動(dòng)成像過(guò)程中，由于相干光經(jīng)過(guò)目標(biāo)表面反射或者通過(guò)無(wú)規(guī)則漲落的折射媒介傳播，會(huì)形成一種信號(hào)相關(guān)的乘性噪聲，即散斑噪聲，造成圖像質(zhì)量嚴(yán)重下降。文獻(xiàn)［4］通過(guò)計(jì)算灰度概率分布曲線證明散斑噪聲是激光主動(dòng)成像的主要噪聲，因此研究抑制散斑噪聲的算法是大多數(shù)激光主動(dòng)成像系統(tǒng)都必須考慮的內(nèi)容。

國(guó)內(nèi)外學(xué)者提出諸多抑制散斑噪聲的經(jīng)典算法，包括 均 值 濾 波［5］、Lee濾 波［6］、Frost濾 波［7］、Kuan濾波［8］及其改進(jìn)算法［9］。這些濾波器使用一定大小的濾波窗口在圖像上滑動(dòng)，使用不同算法計(jì)算出當(dāng)前像素的替代值，其中Lee濾波由于在濾除噪聲和保留信號(hào)完整性之間能夠達(dá)到較好的平衡而倍受關(guān)注。這類窗口算法雖然在同類區(qū)域去除散斑噪聲效果較好，但運(yùn)算時(shí)間較長(zhǎng)，對(duì)于圖像邊緣和細(xì)節(jié)信息的保留能力也不理想。變換域降噪算法中使用較多的是小波軟閾值濾波［10－11］，這類算法缺陷在于閾值的選取較為困難，閾值過(guò)大容易丟失圖像細(xì)節(jié)信息，閾值過(guò)小則不能很好地去除噪聲。近年來(lái)出現(xiàn)了不少針對(duì)散斑噪聲的新算法，如Jarabo－Amores等［12］提出一種基于“mean－shift”濾波的去噪方法，該方法不必事先假設(shè)噪聲的統(tǒng)計(jì)模型而能夠達(dá)到與統(tǒng)計(jì)去噪相當(dāng)?shù)男Ч?。Hyenkyun 等［13］提出基于交互技術(shù)（Shifting technology）的降噪方法，能夠快速有效地去除散斑噪聲。國(guó)內(nèi)，文獻(xiàn)［14］使用SURE（Stein′s unbiased risk estimate）獲得近似最優(yōu)的小波系數(shù)的權(quán)值，該方法優(yōu)勢(shì)在于計(jì)算量很小。

信號(hào)子空間方法的基本思想是將原始信號(hào)通過(guò)K－L變換分解為信號(hào)子空間和噪聲子空間，通過(guò)濾除噪聲子空間噪聲分量的同時(shí)保留信號(hào)子空間的信號(hào)部分，能夠恢復(fù)出較為純凈的信號(hào)，該方法已經(jīng)成功運(yùn)用于語(yǔ)音和雷達(dá)信號(hào)降噪領(lǐng)域［15－16］?；谄鋬?yōu)異的信噪分離效果，本文將該方法從一維推廣到二維空間，提出一種基于時(shí)域約束（TDC）的激光主動(dòng)成像散斑降噪方法。將本文的TDC 降噪方法同經(jīng)典的Lee、Frost和Kuan濾波作比較，結(jié)果表明本文算法具有出色的噪聲抑制效果，同時(shí)算法也具有很好的實(shí)時(shí)性。

1 散斑噪聲分析

被激光照射的目標(biāo)表面相對(duì)于激光波長(zhǎng)而言凹凸不平，相干光被不同位置反射后在像面處發(fā)生干涉，即形成散斑圖樣。散斑噪聲的數(shù)學(xué)模型如式（1）（2）所示：

式中：I表示含噪圖像；I0表示無(wú)噪圖像；Isn表示散斑噪聲作用系數(shù)，服從均值為1，方差為1／L 的伽馬分布，L 是常數(shù)。

由式（1）可以看出散斑噪聲是一種乘性噪聲。為了將信號(hào)和噪聲分離，使用同態(tài)變換將其變?yōu)榧有栽肼暎?/p>

記上式為：

根據(jù)文獻(xiàn)［17］，同態(tài)變換后散斑噪聲的均值和方差的關(guān)系分別如式（5）（6）所示：

2 信號(hào)子空間濾波及TDC濾波方法

2.1 信號(hào)子空間

信號(hào)子空間技術(shù)是一種類似于傅里葉變換［18］、小波變換［19］的子空間變換技術(shù)。其原理是：帶噪信號(hào)的空間由信號(hào)子空間與噪聲子空間構(gòu)成，噪聲能量均勻分布于整個(gè)空間中，而信號(hào)能量則只分布于信號(hào)子空間中。將帶噪信號(hào)分解至信號(hào)子空間和噪聲子空間，通過(guò)消去噪聲子空間的分量，就可以恢復(fù)出近似純凈的信號(hào)。該方法對(duì)于語(yǔ)音信號(hào)具有良好的去相關(guān)性，因而廣泛用于語(yǔ)音增強(qiáng)領(lǐng)域。本文接下來(lái)將這種方法推廣到二維圖像空間。

設(shè)無(wú)噪圖像的估計(jì)為：

則誤差信號(hào)表示為：

定義失真信號(hào)能量和殘余噪聲能量分別為：

式中：tr（）表示取對(duì)角線上的元素。

為使失真信號(hào)能量達(dá)到最小，設(shè)置約束條件為：

式中：σ2是常數(shù)。對(duì)上述約束最小值問(wèn)題，可以采用Kuhn－Tucker條件得到Lagrange多項(xiàng)式：

根據(jù)條件

可以計(jì)算出濾波矩陣H 的估計(jì)值：

式中：RX＝E（XXT）是無(wú)噪圖像的協(xié)方差矩陣，RN＝E（NNT）是 噪 聲 的 協(xié) 方 差 矩 陣，μ 是Lagrange算子。將RX進(jìn)行特征值分解，表示成信號(hào)子空間和噪聲子空間兩部分：

式中：U1、U2分別是m×r和m×（n－r）的特征值矩陣，r是信號(hào)子空間的維數(shù)，n－r是噪聲子空間的維數(shù)；ΛX1、ΛX2是特征值對(duì)角陣；代表的是RX在信號(hào)子空間的分量；代表的是RX在噪聲子空間的分量。

取RX的信號(hào)子空間的分量代入式（15），可得濾波估計(jì)矩陣：

由式（16）（17）可知，為了得到HTDC，需要確定以下四個(gè)參數(shù)：信號(hào)子空間的維數(shù)r、無(wú)噪圖像的協(xié)方差矩陣RX、噪聲的協(xié)方差矩陣RN以及常數(shù)μ。

2.2 確定信號(hào)子空間維數(shù)r

估計(jì)信號(hào)子空間的維數(shù)，可以用奇異值分解實(shí)現(xiàn)。信號(hào)Y 的奇異值分解表示如下：

式中：S 是Y 的 奇 異 值 矩 陣，Sm×m＝diag（s1，s2，…，sm），s1≥s2≥…≥sm，即對(duì)角陣上的值采用降序排列。

信號(hào)某維度的奇異值大小代表該維度的重要程度，奇異值越大，對(duì)應(yīng)該維度越重要。取S的前r個(gè)分量，則Y 可以近似表示為：

文中選取r的條件為：

式中：Tr是一個(gè)接近于1的閾值，具體大小由實(shí)驗(yàn)決定。

由于像素值為0的像素點(diǎn)經(jīng)過(guò)log變換會(huì)出現(xiàn)無(wú)窮小，使Y 無(wú)法進(jìn)行奇異值分解，因此同態(tài)變換前先將每個(gè)像素亮度增加10，待取指數(shù)恢復(fù)之后再減去10。

2.3 確定無(wú)噪圖像的協(xié)方差矩陣RX、噪聲的協(xié)方差矩陣RN

經(jīng)過(guò)同態(tài)變換之后乘性噪聲變?yōu)榧有栽肼暎缡剑?）所示，則可以得到

式中：RY＝E（YYT）可以直接由含噪圖像求得。對(duì)于RN，式（6）給出了它和常數(shù)L的關(guān)系，但是L難以估計(jì)。本文從小波分析的角度出發(fā)：在小波域上，圖像信號(hào)集中在低頻分量，而噪聲信號(hào)均勻分布在整個(gè)小波空間［16］，因此可以利用高頻子帶的小波系數(shù)估計(jì)噪聲的強(qiáng)度：

2.4 確定常數(shù)μ

μ 的取值會(huì)影響濾波后圖像的質(zhì)量。μ 值大，可以較大程度地濾除散斑噪聲，卻會(huì)導(dǎo)致較嚴(yán)重的信號(hào)失真；μ 值小，能夠保留較多的細(xì)節(jié)信號(hào)，但會(huì)殘余更多的噪聲。選擇μ的值需要在信號(hào)完整性和濾波性能之間權(quán)衡。

根據(jù)文獻(xiàn)［16］，采用以下公式確定μ：式中：SNR 是圖像的信噪比；μ0 和s是常數(shù)，由實(shí)驗(yàn)確定。

式（23）表示的是：對(duì)于高信噪比的圖像，噪聲較小，μ取較小的值，以保護(hù)信號(hào)完整性為主要目標(biāo)；對(duì)于低信噪比的圖像，噪聲較大，μ 取較大的值，以濾波為主要目標(biāo)。

2.5 基于信號(hào)子空間的TDC濾波方法

本文基于信號(hào)子空間的散斑去噪方法如下：

（1）對(duì)增加亮度后的圖像進(jìn)行同態(tài)變換，將乘性噪聲變?yōu)榧有栽肼暋?/p>

（2）將圖像進(jìn)行小波變換，并利用高頻分量估計(jì)噪聲方差RN。

（3）對(duì)Y 進(jìn)行奇異值分解，估計(jì)信號(hào)子空間的維數(shù)r。

（4）計(jì)算無(wú)噪圖像的協(xié)方差矩陣RX＝RYRN，并根據(jù)式（16）對(duì)其進(jìn)行特征值分解。

（5）根據(jù)式（17）計(jì)算變換矩陣HTDC，得到無(wú)噪圖像的估計(jì)值

（6）同態(tài)逆變換，得到最終的去噪圖像。

3 實(shí)驗(yàn)驗(yàn)證

實(shí)驗(yàn)搭建一套基于距離選通技術(shù)的激光主動(dòng)照明系統(tǒng)：使用532nm 波長(zhǎng)、5mrad發(fā)散角的半導(dǎo)體泵浦固體脈沖激光器發(fā)射脈沖光束照射3.4 km 處的目標(biāo)，用選通ICCD 接收反射信號(hào)，通過(guò)同步控制器控制選通快門(mén)的開(kāi)閉。其中選通ICCD 由選通像增強(qiáng)器和CCD 組成，具有高增益、高靈敏度、成像清晰等優(yōu)點(diǎn)。同步控制器控制ICCD 在納秒級(jí)的時(shí)間內(nèi)進(jìn)行選通，可以抑制雜散光，保證采集到較為清晰的含散斑圖像。圖1是實(shí)驗(yàn)系統(tǒng)框圖，圖2是實(shí)驗(yàn)系統(tǒng)實(shí)物圖。

圖1 實(shí)驗(yàn)系統(tǒng)框圖Fig.1 Block diagram of experimental system

圖2 實(shí)驗(yàn)系統(tǒng)實(shí)物圖Fig.2 Photo of experimental system

截取ICCD 輸出圖像中被照亮的區(qū)域，最終得到的實(shí)驗(yàn)圖像大小為150×150像素。分別采用Lee、Frost、Kuan和本文的TDC 濾波方法進(jìn)行實(shí)驗(yàn)，Lee、Frost和Kuan的濾波窗口均選為5×5。窗口大小選為5×5可以在濾除噪聲和保留圖像細(xì)節(jié)之間得到較好的平衡。實(shí)驗(yàn)編程環(huán)境為：Windows XP、i7－2600 3.40 GHz CPU、2G內(nèi)存、VS2008＋opencv軟件平臺(tái)。

圖3為4種算法的濾波效果對(duì)比，ICCD 的積分時(shí)間為50ns。

圖3 不同濾波方法的濾波結(jié)果Fig.3 Results of different denoising methods

從圖3直觀上看，Lee、Frost和Kuan結(jié)果圖中殘余較多的散斑噪聲，而TDC 濾波后的圖像最清晰，濾波效果最好。觀察圖像下方中間的門(mén)框，可以發(fā)現(xiàn)TDC 濾波后的邊緣更加清晰，這是因?yàn)槠溆? 種算法均使用鄰域像素進(jìn)行加權(quán)操作，必將導(dǎo)致細(xì)節(jié)模糊。而TDC 算法將圖像矩陣當(dāng)做整體進(jìn)行運(yùn)算，是對(duì)圖像信號(hào)的稀疏表示，信號(hào)子空間的信號(hào)部分沒(méi)有損失，故邊緣保持能力更強(qiáng)。

本文引入峰值信噪比PSNR［20］以評(píng)價(jià)濾波算法的有效性。改變ICCD 的積分時(shí)間，得到不同噪聲強(qiáng)度的圖像，分別用4種濾波方法進(jìn)行實(shí)驗(yàn)，計(jì)算其PSNR 值，比較濾波效果，結(jié)果如圖4所示。

圖4 不同濾波算法的PSNR 曲線Fig.4 PSNR curves of different denoising methods

由圖4看出，本文TDC算法對(duì)于不同強(qiáng)度的噪聲能夠取得比經(jīng)典的Lee、Frost和Kuan算法更好的濾除效果。觀察PSNR 曲線，隨著曝光時(shí)間的增加，TDC濾波算法的曲線與其余曲線的差值越來(lái)越大，這表示隨著噪聲的增大，TDC 的濾波性能優(yōu)勢(shì)相對(duì)于其余算法越來(lái)越突出。這是因?yàn)閷?duì)于Lee、Frost和Kuan 等窗口濾波而言，使用有限大小的濾波窗口信息，并不能很好地估計(jì)噪聲的方差，隨著噪聲的增加，濾波能力并不會(huì)顯著提升，濾波后殘余的噪聲會(huì)越來(lái)越大。而TDC算法的約束條件是將殘余噪聲限制在σ2范圍內(nèi)，算法在整個(gè)圖像上進(jìn)行噪聲強(qiáng)度的估計(jì)，并且只保留信號(hào)子空間中的非常小的噪聲能量，因此濾波效果更好。

表1統(tǒng)計(jì)的是不同算法對(duì)于150×150圖像的平均運(yùn)算時(shí)間。可見(jiàn)本文算法比其余3種算法耗時(shí)要小得多。這是因?yàn)門(mén)DC 算法對(duì)整個(gè)圖像矩陣同時(shí)進(jìn)行計(jì)算，而其余3種算法均需要滑動(dòng)濾波窗口遍歷整個(gè)圖像區(qū)域。在本文的實(shí)驗(yàn)條件下（成像區(qū)域大小為150×150，幀頻為25 幀／s），TDC濾波算法能夠滿足實(shí)時(shí)處理的要求。

表1 不同算法的平均運(yùn)算時(shí)間Table 1 Average running time of different denoising methods

4 結(jié)束語(yǔ)

受語(yǔ)音增強(qiáng)領(lǐng)域的子空間分解方法的啟發(fā)，本文在分析散斑噪聲數(shù)學(xué)模型的基礎(chǔ)上，提出一種基于信號(hào)子空間的激光主動(dòng)成像去噪方法，并搭建了一套基于距離選通技術(shù)的激光主動(dòng)照明系統(tǒng)進(jìn)行實(shí)驗(yàn)驗(yàn)證。本文算法將含噪圖像分解為信號(hào)子空間和噪聲子空間，通過(guò)去除噪聲子空間并盡量保留信號(hào)子空間的能量，以取得良好的去噪效果。結(jié)果證明本文TDC 方法具有比經(jīng)典的Lee、Frost和Kuan算法更好的去噪性能，同時(shí)縮短了計(jì)算時(shí)間，滿足激光主動(dòng)成像系統(tǒng)實(shí)時(shí)性的要求。

［1］王燦進(jìn)，孫濤，石寧寧，等.基于雙隱含層BP算法的激光主動(dòng)成像識(shí)別系統(tǒng)［J］.光學(xué)精密工程，2014，22（6）：1639－1647.Wang Can－jin，Sun Tao，Shi Ning－ning，et al.Laser active imaging and recognition system based on double hidden layer BP algorithm［J］.Opt Precision Eng，2014，22（6）：1639－1647.

［2］錢(qián)方，孫濤，郭勁，等.無(wú)參考的特征點(diǎn)復(fù)雜度激光干擾圖像評(píng)估［J］.光學(xué)精密工程，2015，23（4）：1179－1186.Qian Fang，Sun Tao，Guo Jin，et al.No－reference laser－dazzling image quality assessment based on feature－point complexity［J］.Opt Precision Eng，2015，23（4）：1179－1186.

［3］趙建川，王弟男，陳長(zhǎng)青，等.紅外激光主動(dòng)成像和識(shí)別［J］.中國(guó)光學(xué)，2013，6（5）：795－802.Zhao Jian－chuan，Wang Di－nan，Chen Chang－qing，et al.Infrared laser active imaging and recognition technology［J］.Chinese Optics，2013，6（5）：795－802.

［4］李自勤，李琦，王騏.由統(tǒng)計(jì)特性分析激光主動(dòng)成像系統(tǒng)圖像的噪聲性質(zhì)［J］.中國(guó)激光，2004，31（9）：1081－1085.Li Zi－qin，Li Qi，Wang Qi.Noise characteristic in active laser imaging system by statistic analysis［J］.Chinese Journal of Lasers，2004，31（9）：1081－1085.

［5］Loupas T，McDicken W，Allan P.An adaptive weighted median filter for speckle suppression in medical ultrasonic images［J］.IEEE Trans Circuits System，1989，36（1）：129－135.

［6］Lee J S.Digital image enhancement and noise filtering by use of local statistics［J］.IEEE Trans Pattern Analysis and Machine Intell，1980，20：165－168.

［7］Frost V S，Stiles J A，Shanmugan K S，et al.A mode for radar images and its application to adaptive digital filtering of multiplicative noise［J］.IEEE Trans Pattern Analysis and Machine Intell，1982，4：157－165.

［8］Kuan D，Sawchuk A，Strand T，et al.Adaptive restoration of images with speckle［J］.IEEE Trans Acoust，Speech and Signal Process，1987，35（3）：373－383.

［9］Lu Y H，Tan S Y，Yeo T S，et al.Adaptive filtering algorithms for SAR speckle reduction［C］∥Geoscience and Remote Sensing Symposium，Lincoln，NE，USA，1996.

［10］Donoho D L.Denoising by soft－thresholding［J］.IEEE TransInform Theory，1995，41（3）：613－627.

［11］葉樹(shù)亮，張玉德，張煒.齒輪視覺(jué)檢測(cè)中的尺度與方向相關(guān)性聯(lián)合降噪［J］.光學(xué)精密工程，2014，22（6）：1622－1630.Ye Shu－liang，Zhang Yu－de，Zhang Wei.Scale and directional correlation combined denoiseing in gear visual inspection［J］.Opt Precision Eng，2014，22（6）：1622－1630.

［12］Jarabo－Amores P，Rosa－Zurera M，Mata－Moya D，et al.Mean－shift filtering to reduce speckle noise in SAR images［C］∥Instrumentation and Measurement Technology Conference，Singapore，2009.

［13］Hyenkyun W，Yun S.Alternating minimization algorithm for speckle reduction with a shifting technique［J］.Image Processing，2012，21（4）：1701－1714.

［14］李曉峰，徐軍，羅積軍，等.激光主動(dòng)成像圖像噪聲分析與抑制［J］.紅外與激光工程，2011，40（2）：332－337.Li Xiao－feng，Xu Jun，Luo Ji－jun，et al.Noise analyzing and denoising of intensity image for laser active imaging system［J］.Infrared and Laser Engineering，2011，40（2）：332－337.

［15］Ikeda S，Ohtsuki T，Tsuji H.Signal－subspace－partition event filtering for eigenvector－based security system using radio waves［C］∥Personal，Indoor and Mobile Radio Communications，Tokyo，Japan，2009.

［16］Gu J，Yang J，Zhang H，et al.Speckle filtering in polarimetric SAR data based on the subspace decomposition［J］.Geoscience and Remote Sensing，2004，42（8）：1635－1641.

［17］Donoho D L，Johnstone I M.Ideal spatial adaptation by wavelet shrinkage［J］.Biometrika，1994，81（3）：425－455.

［18］鄧建青，劉晶紅.基于Fourier－Mellin變換和Keren算法的改進(jìn)運(yùn)動(dòng)估計(jì)算法［J］.液晶與顯示，2011，26（3）：364－369.Deng Jian－qing，Liu Jing－h(huán)ong.Improved motion estimation algorithm based on Fourier－Mellin transform and Keren algorithm［J］.Chinese Journal of Liquid Crystals and Displays，2011，26（3）：364－369.

［19］任志英，高誠(chéng)輝，申丁，等.雙樹(shù)復(fù)小波穩(wěn)健濾波在工程表面粗糙度評(píng)定中的應(yīng)用［J］.光學(xué)精密工程，2014，22（7）：1820－1827.Ren Zhi－ying，Gao Cheng－h(huán)ui，Shen Ding，et al.Application of DT－CWT robust filtering to evaluation of engineering surface roughness［J］.Opt Precision Eng，2014，22（7）：1820－1827.

［20］任文琦，王元全.基于梯度矢量卷積場(chǎng)的四階各向異性擴(kuò)散及圖像去噪［J］.光學(xué)精密工程，2013，21（10）：2713－2719.Ren Wen－qi，Wang Yuan－quan.GVC－based fourthorder anisotropic diffusion for image denoising［J］.Opt Precision Eng，2013，21（10）：2713－2719.