基于噪聲可見(jiàn)性函數(shù)的SAR圖像增強(qiáng)快速算法

朱逸飛,楊 國(guó),王 強(qiáng),吳 文

(南京理工大學(xué) 近程高速目標(biāo)探測(cè)技術(shù)國(guó)防重點(diǎn)學(xué)科實(shí)驗(yàn)室，南京 210094)

0 引言

合成孔徑雷達(dá)SAR(synthetic aperture radar)是一種高分辨率成像雷達(dá)，具有全天時(shí)、全天候、穿透性強(qiáng)等特點(diǎn)，廣泛應(yīng)用于軍事和民用領(lǐng)域。但由于戰(zhàn)場(chǎng)環(huán)境復(fù)雜、大氣衰減、硬件設(shè)備干擾等原因，SAR成像結(jié)果會(huì)存在整體灰度分布低且集中、圖像對(duì)比度低、含有噪聲等問(wèn)題，造成目標(biāo)不清晰、圖像模糊，影響圖像識(shí)別效果，因此需要對(duì)SAR圖像進(jìn)行增強(qiáng)[1]。目前國(guó)內(nèi)外利用數(shù)據(jù)處理的方法提高SAR圖像質(zhì)量主要分為兩類[2]：一是在成像階段利用系統(tǒng)回波數(shù)據(jù)，通過(guò)超分辨成像算法獲得高分辨率的SAR圖像，即SAR的超分辨成像方法；二是在成像后的圖像數(shù)據(jù)基礎(chǔ)上進(jìn)行圖像質(zhì)量提高的處理。其中后者是本文的研究?jī)?nèi)容，典型算法如Kuan算法、Lee算法、雙邊濾波算法、偏微分算法等，這些算法能夠?yàn)V除部分噪聲并保留細(xì)節(jié)邊緣和目標(biāo)特征[3]。但是，典型的SAR圖像增強(qiáng)算法復(fù)雜度較高，面對(duì)高分辨率圖像時(shí)，需要的計(jì)算時(shí)間較長(zhǎng)，在目標(biāo)識(shí)別、跟蹤及災(zāi)難評(píng)估等對(duì)實(shí)時(shí)性要求較高的領(lǐng)域，無(wú)法滿足系統(tǒng)要求。

與中央處理器CPU(central processing unit)不同，圖形處理器GPU(graphic processing unit)專為密集型、高度并行化的計(jì)算而設(shè)計(jì)，其中用于數(shù)據(jù)處理的晶體管數(shù)量遠(yuǎn)遠(yuǎn)大于緩存和邏輯控制部分，這使得GPU更適合處理無(wú)邏輯關(guān)系數(shù)據(jù)的并行運(yùn)算[4]。2007年，NVIDIA公司推出通用并行計(jì)算架構(gòu)CUDA(compute unified device architecture)，目的是將GPU作為并行計(jì)算設(shè)備，進(jìn)行通用并行計(jì)算[5]。CUDA為開(kāi)發(fā)者提供了硬件的直接訪問(wèn)接口，可以方便地寫(xiě)出在GPU上執(zhí)行的程序，而不用像過(guò)去的GPGPU架構(gòu)，如OpenGL等，將計(jì)算映射到圖形API中，大大降低了開(kāi)發(fā)門(mén)檻和難度。自推出后，CUDA被廣泛應(yīng)用于計(jì)算機(jī)可視化、音視頻編解碼、流體力學(xué)模擬等領(lǐng)域，都能獲得較高的加速比，有著良好的應(yīng)用前景[6]。

本文提出一種基于噪聲可見(jiàn)性函數(shù)的SAR圖像增強(qiáng)快速算法及其基于GPU的并行化實(shí)現(xiàn)方法。該算法運(yùn)用圖像分層理論，分離SAR圖像細(xì)節(jié)與噪聲，在對(duì)圖像細(xì)節(jié)層進(jìn)行處理時(shí)，結(jié)合人眼視覺(jué)特性，引入噪聲可見(jiàn)性函數(shù)，控制細(xì)節(jié)層增益系數(shù)，實(shí)現(xiàn)圖像增強(qiáng)。在提高圖像質(zhì)量、增強(qiáng)圖像細(xì)節(jié)的同時(shí)，充分運(yùn)用GPU的并行計(jì)算特點(diǎn)，具有較高的實(shí)時(shí)性。

1 基于圖像分層的增強(qiáng)算法

為了抑制圖像噪聲、增強(qiáng)圖像細(xì)節(jié)、提高圖像識(shí)別效果，Branchitta F等人提出了圖像分層理論[7]。其基本思想為：首先，選擇一種分層濾波器對(duì)原始輸入圖像進(jìn)行處理，抹平圖像細(xì)節(jié)和微小波動(dòng)，得到圖像基本層和細(xì)節(jié)層；其次，分別對(duì)圖像基本層和細(xì)節(jié)層進(jìn)行處理，因?yàn)榛緦又胁话瑘D像細(xì)節(jié)信息，所以處理過(guò)程中不必考慮細(xì)節(jié)損失問(wèn)題，重點(diǎn)在于圖像對(duì)比度的調(diào)整；由于細(xì)節(jié)層不僅包含原始圖像的細(xì)節(jié)信息，還有噪聲的存在，因此需要對(duì)細(xì)節(jié)層進(jìn)行噪聲抑制和細(xì)節(jié)增強(qiáng)；最后，將處理后的基本層和細(xì)節(jié)層融合并調(diào)整顯示范圍，得到輸出圖像。

研究基于圖像分層的增強(qiáng)算法，其重點(diǎn)在于分層濾波器、基本層圖像處理方法和細(xì)節(jié)層圖像增強(qiáng)方法的選擇。本文采用雙邊濾波作為分層濾波器，對(duì)SAR圖像進(jìn)行雙邊濾波，得到抹平圖像細(xì)節(jié)和微小波動(dòng)后的基本層，并與原圖像相減，得到包含細(xì)節(jié)信息和大量噪聲的細(xì)節(jié)層；對(duì)基本層采用γ<1的伽馬變換；通過(guò)噪聲可見(jiàn)性函數(shù)控制細(xì)節(jié)層增益系數(shù)。算法框架如圖1所示。

圖1 基于圖像分層的增強(qiáng)算法框架

雙邊濾波(Bilateralfilter)結(jié)合空域信息和灰度相似性，對(duì)圖像進(jìn)行基于空間分布的高斯濾波，是一種非線性的濾波方法。鄰域內(nèi)離中心像素較遠(yuǎn)的像素占有的權(quán)重較小，不會(huì)對(duì)中心像素的值造成太大影響，這樣有利于圖像細(xì)節(jié)邊緣像素值的保存，達(dá)到保邊去噪的目的。

設(shè)f(x,y)為原圖像，g(x,y)為原圖像經(jīng)雙邊濾波后的結(jié)果，即基本層，其表達(dá)式為：

(1)

式中，Sx,y表示以像素(x,y)為中心，半寬為N的鄰域，ω(i,j)為兩部分系數(shù)的乘積：

ω(i,j)=ωs(i,j)ωr(i,j)

(2)

式中，ωs(i,j)為空域?yàn)V波器系數(shù)，ωr(i,j)為值域?yàn)V波器系數(shù)：

(3)

(4)

對(duì)基本層g(x,y)采用γ<1的伽馬變換，拉伸圖像灰度范圍，提高對(duì)比度。伽馬變換，也叫做冪律變換，是一種非線性的灰度拉伸變換。設(shè)變換后的基本層為gγ(x,y)，變換形式為：

gγ(x,y)=C*g(x,y)γ

(5)

式中，γ為變換系數(shù)，C為控制像素變化范圍的常數(shù)。

2 基于NVF的細(xì)節(jié)層增強(qiáng)

在研究水印估計(jì)問(wèn)題時(shí)，S.Voloshnovskiy等人提出了噪聲可見(jiàn)性函數(shù)NVF(noise visibility function)[8-9]：

(6)

歸一化處理后，式(6)變?yōu)椋?/p>

(7)

式中，θ為常數(shù)，m(x,y)為噪聲掩膜函數(shù)，與圖像的局部方差有關(guān)，可用于衡量圖像細(xì)節(jié)。

由式(7)可以看出，NVF(x,y)可以表示圖像中各像素點(diǎn)對(duì)噪聲的敏感程度，與圖像局部能量成反比。越是平坦的區(qū)域，m(x,y)越小，NVF(x,y)越接近1，圖像在該像素處較為敏感，噪聲對(duì)人眼刺激較大，視覺(jué)效果降低；越是變化劇烈的區(qū)域，m(x,y)越大，NVF(x,y)越接近0，圖像在該像素處允許較大的噪聲，其對(duì)圖像視覺(jué)效果的影響不明顯，噪聲表現(xiàn)出低可見(jiàn)性，不易被察覺(jué)。

將原圖像f(x,y)與基本層g(x,y)相減得到細(xì)節(jié)層圖像h(x,y)。h(x,y)中包含了原圖像中的細(xì)節(jié)信息，必須進(jìn)行保留或增強(qiáng)，但細(xì)節(jié)層中又不可避免地存在原圖像中的大部分噪聲。因此，在對(duì)細(xì)節(jié)層進(jìn)行處理時(shí)，不能直接簡(jiǎn)單地放大，這會(huì)造成平坦區(qū)域噪聲的過(guò)度放大。

增強(qiáng)細(xì)節(jié)層圖像時(shí)，通過(guò)NVF區(qū)分圖像的平坦與變化劇烈的區(qū)域，對(duì)不同的區(qū)域采用不同的增益系數(shù)。設(shè)細(xì)節(jié)層像素點(diǎn)(x,y)的增益系數(shù)為G(x,y)，對(duì)于平坦區(qū)域的像素點(diǎn)，NVF趨近于1，G(x,y)需要取較小值，避免噪聲被增強(qiáng)；對(duì)于變化劇烈的區(qū)域，NVF趨近于0，G(x,y)可以取較大值，在增強(qiáng)細(xì)節(jié)時(shí)不用考慮噪聲的問(wèn)題。設(shè)Gmax為G(x,y)的最大值，Gmin為G(x,y)的最小值，構(gòu)造增益系數(shù)公式如下：

G(x,y)=Gmin+(1-NVF(x,y))(Gmax-Gmin)

(8)

(9)

且有：

m(x,y)=k(x,y)-1-1

(10)

因?yàn)閗(x,y)∈[0，1]，無(wú)需調(diào)整范圍，所以取θ=1。由式(7)、式(8)和式(10)可得：

G(x,y)=Gmin+(1-k(x,y))(Gmax-Gmin)

(11)

增強(qiáng)后的細(xì)節(jié)層hNVF(x,y)可表示為：

hNVF(x,y)=h(x,y)·G(x,y)

(12)

3 基于CUDA的圖像增強(qiáng)算法實(shí)現(xiàn)

3.1 CUDA編程模型

GPU擁有強(qiáng)大的并行計(jì)算能力，但早期的GPU編程存在很多困難，如圖形API、交互接口復(fù)雜等，這些都限制了GPU在通用計(jì)算上的發(fā)展[10]。CUDA的出現(xiàn)解決了這一問(wèn)題，其包含一整套軟硬件體系，支持通用并行計(jì)算。由于CUDA采用類C語(yǔ)言進(jìn)行開(kāi)發(fā)，具有高性能計(jì)算指令開(kāi)發(fā)能力，開(kāi)發(fā)者能夠在GPU強(qiáng)大的并行計(jì)算能力基礎(chǔ)上建立更高效率的密集數(shù)據(jù)計(jì)算方案。

一個(gè)完整的CUDA程序包含在CPU(主機(jī)端，host)中執(zhí)行的部分和在GPU(設(shè)備端，device)中執(zhí)行的部分。如圖2所示，主機(jī)端負(fù)責(zé)數(shù)據(jù)的輸入輸出、內(nèi)存分配、主機(jī)與設(shè)備間的數(shù)據(jù)傳輸；而在設(shè)備端中運(yùn)行的核函數(shù)(kernel)承擔(dān)著并行計(jì)算的任務(wù)[11]。由于CPU對(duì)顯存中的數(shù)據(jù)進(jìn)行存取時(shí)只能通過(guò)PCIExpress接口，其理論帶寬為雙向各4GB/s，速度較慢，若頻繁地進(jìn)行存取顯存的操作會(huì)降低GPU計(jì)算效率[12]。

圖2 CUDA編程模型

在并行線程總數(shù)確定的情況下，為了有效利用GPU的計(jì)算資源，需要對(duì)核函數(shù)進(jìn)行合理的線程配置，根據(jù)GPU規(guī)格和需要計(jì)算的數(shù)據(jù)規(guī)模，合理地設(shè)定線程塊Block的數(shù)量和每個(gè)線程塊中線程Thread的數(shù)量。每個(gè)線程塊中線程數(shù)量應(yīng)為該GPU的warp大小的整數(shù)倍[13]。warp為最小單指令多線程執(zhí)行單元，以32個(gè)線程為一組，并行創(chuàng)建、調(diào)度和執(zhí)行[14]，考慮存儲(chǔ)器大小的限制，本文設(shè)定每個(gè)線程塊中線程數(shù)量為256，并根據(jù)不同的圖像大小劃分線程塊數(shù)量。

3.2 并行化圖像增強(qiáng)算法實(shí)現(xiàn)

上文中圖像增強(qiáng)算法復(fù)雜度較高，對(duì)于高分辨率圖像而言，串行實(shí)現(xiàn)方法實(shí)時(shí)性較低。該算法中像素點(diǎn)間的運(yùn)算相對(duì)獨(dú)立，數(shù)據(jù)依賴程度較低，具有較高的內(nèi)在并行性，適合GPU并行實(shí)現(xiàn)。在并行化雙邊濾波算法的基礎(chǔ)上[15]進(jìn)行線程內(nèi)的分層處理、基本層處理和基于NVF的細(xì)節(jié)層增強(qiáng)，并行完成各像素的全部計(jì)算過(guò)程，實(shí)現(xiàn)算法并行化設(shè)計(jì)。

具體流程如下：

1)獲取圖像數(shù)據(jù)，初始化變量，申請(qǐng)GPU顯存空間：cudaMalloc((void**)&Rdata,(width*height*sizeof(double)))；

2)完成CPU與GPU間的數(shù)據(jù)通信；

3)定義線程配置：

Dim3grid(blockx,blocky,1);

Dim3block(threadx,thready,1);

4)啟動(dòng)kernel函數(shù)，進(jìn)行并行化雙邊濾波和分層，根據(jù)式(11)計(jì)算細(xì)節(jié)層增益系數(shù)，完成細(xì)節(jié)層和基本層處理以及圖像融合；

5)實(shí)現(xiàn)圖像輸出，釋放內(nèi)存和顯存空間。

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

本文將從算法效果和加速性能兩方面進(jìn)行實(shí)驗(yàn)。硬件平臺(tái)：CPU為IntelCorei5 6400，4核，主頻2.70GHz；GPU為NVIDIAGeForceGTX960，1024核，運(yùn)行頻率1.22GHz；編譯環(huán)境為VisualStudio2012 和CUDA7.5。

4.1 算法效果

本次實(shí)驗(yàn)中，考慮視覺(jué)效果和計(jì)算復(fù)雜度，取雙邊濾波參數(shù)σs=40、σr=20，鄰域半寬N=5，基本層伽馬變換系數(shù)γ=0.5，細(xì)節(jié)層增益系數(shù)最大值Gmax=1.5，最小值Gmin=1。算法實(shí)驗(yàn)效果如圖3所示。

圖3 算法效果

圖3(a)為1024×1024的原始SAR圖像，對(duì)圖像添加噪聲以模擬SAR成像結(jié)果，結(jié)果如圖3(b)所示。圖3(c)為經(jīng)雙邊濾波分層后得到的未處理的基本層，圖像的細(xì)節(jié)和微小波動(dòng)被抹平；圖3(d)為經(jīng)過(guò)伽馬變換后的基本層，圖像對(duì)比度得到增強(qiáng)；圖3(e)為未處理的細(xì)節(jié)層，圖像細(xì)節(jié)和噪聲同時(shí)存在；圖3(f)為經(jīng)過(guò)增強(qiáng)后的細(xì)節(jié)層，可以看出，在細(xì)節(jié)增強(qiáng)的同時(shí)，部分噪聲得到了抑制；圖3(g)和圖3(h)分別為本文所提算法的串行實(shí)現(xiàn)計(jì)算結(jié)果和并行化實(shí)現(xiàn)方法的計(jì)算結(jié)果，可以看到，圖像對(duì)比度提高的同時(shí)，噪聲得到了抑制，農(nóng)田部分的紋理變得更加清晰，房屋、河流等輪廓也變得十分明顯。

本文通過(guò)峰值信噪比PSNR來(lái)客觀評(píng)價(jià)圖像質(zhì)量[16]，如式(13)所示。其中m、n為圖像維度，f0(i,j)表示原圖像，f1(i,j)表示待評(píng)價(jià)圖像。PSNR值越大，說(shuō)明失真越少，圖像質(zhì)量越高。結(jié)果如表1所示。

(13)

表1 圖像PSNR值比較

從表1可以看出，本文所提算法能夠有效增強(qiáng)圖像細(xì)節(jié)，提高目標(biāo)識(shí)別能力，有較好的增強(qiáng)效果，并且對(duì)圖像質(zhì)量的改善能力優(yōu)于單獨(dú)的雙邊濾波；串行實(shí)現(xiàn)和并行實(shí)現(xiàn)不會(huì)影響算法處理效果，兩者得到的圖像質(zhì)量是一致的。

4.2 加速性能

分別對(duì)不同尺寸的SAR圖像進(jìn)行加速性能測(cè)試，記錄算法串行實(shí)現(xiàn)和并行實(shí)現(xiàn)的計(jì)算時(shí)間。為了減小誤差，提高準(zhǔn)確程度，取10次運(yùn)算結(jié)果的平均值，計(jì)算加速比。實(shí)驗(yàn)結(jié)果如表2所示。

表2 計(jì)算時(shí)間對(duì)比

從表2可以看出，本文所提算法在GPU上的執(zhí)行時(shí)間明顯少于在CPU上的執(zhí)行時(shí)間，基于CUDA的并行化實(shí)現(xiàn)方法具有顯著的加速效果，1024×1024的圖像的加速比達(dá)到75.756，這是因?yàn)楸疚乃崴惴ň哂休^好的并行特性，能夠充分利用GPU的架構(gòu)和存儲(chǔ)結(jié)構(gòu)特點(diǎn)，并行計(jì)算各像素的處理過(guò)程；隨著圖像尺寸的增大，加速比逐漸提高，這是因?yàn)橄到y(tǒng)調(diào)度和數(shù)據(jù)傳輸需要一定的時(shí)間開(kāi)銷，數(shù)據(jù)規(guī)模較大時(shí)，雖然調(diào)用的線程較多，數(shù)據(jù)傳輸時(shí)間較長(zhǎng)，但其在整體運(yùn)行時(shí)間中所占比重較小，因此加速比較大。

5 結(jié)論

本文提出一種基于噪聲可見(jiàn)性函數(shù)的SAR圖像增強(qiáng)快速算法及其基于GPU的并行化實(shí)現(xiàn)方法。實(shí)驗(yàn)結(jié)果表明，該算法能夠充分利用圖像分層的內(nèi)在并行性，實(shí)現(xiàn)并行加速，在有效提高圖像質(zhì)量，增強(qiáng)圖像細(xì)節(jié)的同時(shí)，減少計(jì)算時(shí)間，進(jìn)一步滿足SAR圖像增強(qiáng)系統(tǒng)實(shí)時(shí)性要求。

[1] 陳建飛,李躍華,王劍橋. 毫米波被動(dòng)合成孔徑成像系統(tǒng)[J].微波學(xué)報(bào),2012,S1:253-256.

[2] 尹奎英. SAR圖像處理及地面目標(biāo)識(shí)別技術(shù)研究[D]. 西安：西安電子科技大學(xué),2011.

[3] 陳 沫. SAR圖像軍事目標(biāo)識(shí)別方法研究[D]. 合肥：安徽大學(xué),2013.

[4] 張 舒. GPU高性能運(yùn)算之CUDA[M]. 北京：中國(guó)水利水電出版社,2009.

[5] 姚 旺,胡 欣,劉 飛,等. 基于GPU的高性能并行計(jì)算技術(shù)[J]. 計(jì)算機(jī)測(cè)量與控制,2014,22(12):4160-4162.

[6] 宋曉麗,王 慶. 基于GPGPU的數(shù)字圖像并行化預(yù)處理[J]. 計(jì)算機(jī)測(cè)量與控制,2009,17(6):1169-1171.

[7] Branchitta F,Diani M,Romagnoli M. New technique for the visualization of high dynamic range infrared images[J]. Optical Engineering,2009,48(9):6401.

[8] Voloshynovskiy S,Pereira S,Iquise V,et al. Attack modelling: towards a second generation watermarking benchmark[J]. Signal Processing,2001,81(6):1177-1214.

[9] 湯代佳,尚東方. 基于噪聲可見(jiàn)性函數(shù)的交通流量異常預(yù)警算法[J]. 上海電機(jī)學(xué)院學(xué)報(bào),2012,15(6):386-389.

[10] 班陽(yáng)陽(yáng),張勁東,陳家瑞,等. 后向投影成像算法的 GPU 優(yōu)化方法研究[J]. 雷達(dá)科學(xué)與技術(shù),2014(6):659-665.

[11] 張紹良,閆鈞華,劉 成,等. 基于CUDA的紅外圖像快速增強(qiáng)算法研究[J]. 電子設(shè)計(jì)工程,2012,20(17):153-157.

[12] Zhao C,You W,Wang Y,et al. GPU implementation for real-time hyperspectral anomaly detection[A]. IEEE International Conference on Digital Signal Processing[C]. IEEE,2015:940-943.

[13] Yiu B Y S,Tsang I K H,Yu A C H. Real-time GPU-based software beamformer designed for advanced imaging methods research[C]. Ultrasonics Symposium. IEEE,2011:1920-1923.

[14] 李大禹,胡立發(fā),穆全全,等. CUDA架構(gòu)下的液晶自適應(yīng)波面數(shù)值解析[J]. 光學(xué)精密工程,2010,18(4):848-854.

[15] 曾炫杰,陳強(qiáng),譚海鵬,等. 基于CUDA的加速雙邊濾波算法[J]. 計(jì)算機(jī)科學(xué),2015,42(s1).

[16] 段 群,吳粉俠,劉小豫. 基于CUDA的椒鹽噪聲自適應(yīng)中值濾波算法[J]. 航空計(jì)算技術(shù),2013,43(2):60-63.