汪陳躍，雷旭峰，李澤民，楊紹明，何 燕

（昆明物理研究所，云南 昆明 650223）

0 引言

空間變化離焦模糊是紅外系統(tǒng)普遍存在的一個(gè)問題，產(chǎn)生這一問題的根本原因是受紅外光學(xué)原理限制，視場(chǎng)中物距差別大的物體無法同時(shí)清晰成像。空間離焦模糊紅外圖像復(fù)原算法的目的就是通過圖像處理來減弱離焦導(dǎo)致的圖像模糊，提升圖像質(zhì)量。

理論研究認(rèn)為模糊圖像是由清晰圖像與一個(gè)點(diǎn)擴(kuò)散函數(shù)（point spread function，PSF）卷積并疊加上噪聲的結(jié)果。對(duì)于實(shí)際模糊圖像復(fù)原問題，模糊圖像已知，但模糊圖像的點(diǎn)擴(kuò)散函數(shù)是未知的，且點(diǎn)擴(kuò)散函數(shù)在空間分布上是不一致的，這種模糊稱為空間變化離焦模糊。為了復(fù)原模糊圖像，需估計(jì)圖像的點(diǎn)擴(kuò)散函數(shù)，為了簡(jiǎn)化估計(jì)點(diǎn)擴(kuò)散函數(shù)問題，可采用理論模型來近似實(shí)際點(diǎn)擴(kuò)散函數(shù)，離焦模糊的點(diǎn)擴(kuò)散函數(shù)的數(shù)學(xué)模型一般可以采用圓盤模型或高斯模型[1]，使用近似模型能夠?qū)Ⅻc(diǎn)擴(kuò)散函數(shù)估計(jì)簡(jiǎn)化為估計(jì)一個(gè)參數(shù)（圓盤半徑或高斯方差），高斯模型更接近實(shí)際圖像中的點(diǎn)擴(kuò)散函數(shù)，所以本文研究中選用高斯模型。

空間變化離焦模糊的點(diǎn)擴(kuò)散函數(shù)在空間域上不一致，無法直接使用某一個(gè)點(diǎn)擴(kuò)散函數(shù)來復(fù)原圖像，因而目前空間變化離焦模糊圖像復(fù)原算法主要有兩種方法，第一種是將圖像按模糊程度進(jìn)行分割，對(duì)分割后每個(gè)區(qū)域的點(diǎn)擴(kuò)散函數(shù)進(jìn)行精確估計(jì)，然后使用估計(jì)的點(diǎn)擴(kuò)散函數(shù)對(duì)對(duì)應(yīng)區(qū)域的圖像進(jìn)行復(fù)原處理，最后將復(fù)原后的區(qū)域圖像組合為完整的復(fù)原圖像[2-4]；第二種是對(duì)圖像的每一個(gè)像素的點(diǎn)擴(kuò)散函數(shù)進(jìn)行估計(jì)，并根據(jù)估計(jì)的點(diǎn)擴(kuò)散函數(shù)對(duì)每個(gè)像素進(jìn)行復(fù)原并組合為最終的復(fù)原圖像[5]。

由于點(diǎn)擴(kuò)散函數(shù)估計(jì)要達(dá)到較高的精度需要的計(jì)算量大，難以實(shí)現(xiàn)快速算法。為了實(shí)現(xiàn)空間變化離焦模糊快速算法且滿足實(shí)時(shí)性要求，本文提出使用固定的點(diǎn)擴(kuò)散函數(shù)組對(duì)圖像進(jìn)行復(fù)原，然后對(duì)復(fù)原后的圖像組進(jìn)行像素級(jí)圖像質(zhì)量評(píng)價(jià)，利用評(píng)價(jià)參數(shù)將多幅復(fù)原圖像中的最佳像素組合為最終復(fù)原圖像。模糊圖像復(fù)原算法選用截?cái)嗉s束最小二乘法（truncated constrained least-squares，TCLS）復(fù)原算法，復(fù)原圖像的圖像質(zhì)量評(píng)價(jià)采用基于復(fù)原圖像空間二階導(dǎo)（拉普拉斯算子）的特點(diǎn)提出新的評(píng)價(jià)方法。

1 離焦模糊圖像快速?gòu)?fù)原算法

1.1 離焦模糊圖像理論模型

模糊圖像的數(shù)學(xué)模型可用g＝f?k＋n表示。式中：g是模糊圖像；f是清晰圖像；k是點(diǎn)擴(kuò)散函數(shù)；n是噪聲；?是卷積符號(hào)。離焦模糊點(diǎn)擴(kuò)散函數(shù)的數(shù)學(xué)模型一般采用圓盤分布或高斯分布[1]。

圓盤分布來源于幾何光學(xué)中的模型，根據(jù)高斯成像公式和幾何相似性，可以推導(dǎo)出點(diǎn)光源成像為一個(gè)圓盤，其圓盤直徑由焦距、物距等共同決定[1]。

圓盤離焦模型的點(diǎn)擴(kuò)散函數(shù)可以表示為[1]：

式中：R是模糊半徑。

高斯模型來源于算法的實(shí)際實(shí)踐，暫未有明確的光學(xué)理論原理解釋，它的點(diǎn)擴(kuò)散函數(shù)數(shù)學(xué)表達(dá)如下[1]：

式中：σ是模糊參數(shù)，即高斯模糊標(biāo)準(zhǔn)差。

高斯模型更接近實(shí)際圖像中的點(diǎn)擴(kuò)散函數(shù)，所以研究中使用高斯模型。

1.2 TCLS復(fù)原算法

本文提出的算法需使用多個(gè)不同的點(diǎn)擴(kuò)散函數(shù)分別復(fù)原紅外圖像，因此需選擇能夠快速實(shí)現(xiàn)的復(fù)原算法以減少計(jì)算量，由于一般的復(fù)原算法都需轉(zhuǎn)換到變換域進(jìn)行復(fù)原（如維納濾波復(fù)原算法）或需要通過多次迭代（如Lucy-Richardson 圖像復(fù)原算法）來實(shí)現(xiàn)，不利于算法的硬件實(shí)現(xiàn)，本文選擇TCLS算法作為復(fù)原算法，TCLS算法的復(fù)原過程為空間域卷積計(jì)算，便于硬件實(shí)現(xiàn)快速處理。

TCLS算法由Kim 提出，TCLS算法將約束最小二乘法（constrained least-squares，CLS）在變換域的操作轉(zhuǎn)化到空間域，該算法執(zhí)行速度快[6]。

CLS的復(fù)原公式為[6]：

式中：H(u,v)是點(diǎn)擴(kuò)散函數(shù)h(x,y)的傅里葉變換；R(u,v)是空間域復(fù)原濾波器r(x,y)的傅里葉變換；P(u,v)是的傅里葉變換；λ是與圖像噪聲相關(guān)的系數(shù)，簡(jiǎn)稱噪聲系數(shù)，取值范圍(0,1)，圖像噪聲越小，噪聲系數(shù)取值越小。

CLS復(fù)原濾波器是頻域操作，使用傅里葉逆變換可以將它轉(zhuǎn)換到空間域，如下式[6]：

式中：F－1是傅里葉逆變換；r(x,y)是空間域復(fù)原濾波器。由于一個(gè)超大尺寸的濾波器計(jì)算是不現(xiàn)實(shí)的，我們需要對(duì)它進(jìn)行截?cái)?，Kim 提出使用升余弦函數(shù)進(jìn)行截?cái)?，相比高斯函?shù)，它能更好地保留高頻細(xì)節(jié)。最終的復(fù)原操作可以表達(dá)為[6]：

式中：rT(x,y)是截?cái)嗪蟮腃LS復(fù)原濾波器；g(x,y)是輸入的模糊圖像；fT(x,y)是初步復(fù)原的圖像。

TCLS 將復(fù)原操作轉(zhuǎn)化為卷積計(jì)算，減少了傅里葉變換過程，在實(shí)時(shí)成像系統(tǒng)中可以實(shí)現(xiàn)快速處理。

針對(duì)TCLS算法復(fù)原圖像存在的振鈴和噪聲問題，Kim 采用了空間自適應(yīng)噪聲平滑（spatially adaptive noise smoothing,SANS）算法進(jìn)行處理，它利用輸入模糊圖像的局部方差將輸入模糊圖像和復(fù)原后的圖像進(jìn)行融合，融合方程和依據(jù)局部方差計(jì)算的比例參數(shù)見下式[6]：

TCLS算法的優(yōu)點(diǎn)是復(fù)原過程為空間域的卷積計(jì)算，不需要對(duì)輸入圖像進(jìn)行傅里葉變換或小波變換，硬件上可以采用有限脈沖響應(yīng)濾波器實(shí)現(xiàn)，所以TCLS復(fù)原算法的硬件實(shí)現(xiàn)難度較小。TCLS算法的缺點(diǎn)是隨著需要復(fù)原的模糊的標(biāo)準(zhǔn)差增大，空間域復(fù)原濾波器的尺寸也增大，濾波器尺寸過大會(huì)導(dǎo)致計(jì)算量少的優(yōu)勢(shì)消失，所以TCLS算法適合于模糊程度較小時(shí)，有一定的局限性。

圖1是TCLS 處理的效果，其中(a)是調(diào)焦清晰圖，(b)是同場(chǎng)景調(diào)焦模糊圖，(c)是(b)經(jīng)過TCLS復(fù)原后的圖像（復(fù)原濾波器對(duì)應(yīng)標(biāo)準(zhǔn)差為1.8，噪聲系數(shù)0.0001），復(fù)原圖像細(xì)節(jié)得到增強(qiáng)，但噪聲和振鈴明顯，(d)是(c)經(jīng)SANS 處理后的圖像，SANS 實(shí)現(xiàn)了降噪和去振鈴，去振鈴后的圖像能夠使用簡(jiǎn)單的圖像質(zhì)量評(píng)價(jià)方法進(jìn)行評(píng)價(jià)。

圖1 TCLS算法復(fù)原效果對(duì)比Fig.1 Comparison of image of blurred and restored by TCLS

2 空間變化離焦模糊圖像復(fù)原算法

2.1 提出的空間變化離焦模糊圖像復(fù)原方法

Cheong 利用模糊圖像的局部梯度對(duì)每個(gè)像素的點(diǎn)擴(kuò)散函數(shù)進(jìn)行估計(jì)，然后根據(jù)估計(jì)的點(diǎn)擴(kuò)散函數(shù)參數(shù)從預(yù)設(shè)的多個(gè)TCLS 濾波器中選擇并計(jì)算最終使用的TCLS 濾波器，實(shí)現(xiàn)了不需迭代的快速?gòu)?fù)原[5]。Cheong的算法的局限性是該算法的點(diǎn)擴(kuò)散函數(shù)估計(jì)推導(dǎo)是基于孤立邊的，當(dāng)圖像復(fù)雜時(shí)，點(diǎn)擴(kuò)散函數(shù)的估計(jì)精度下降明顯[5]。為了去掉點(diǎn)擴(kuò)散函數(shù)估計(jì)的流程，本文提出針對(duì)不同參數(shù)TCLS復(fù)原的圖像進(jìn)行區(qū)域圖像質(zhì)量評(píng)價(jià)的方法，將不同參數(shù)復(fù)原圖像中的最佳像素組合為最終的輸出圖像。

本文提出的方法的流程圖見圖2。首先使用m個(gè)不同標(biāo)準(zhǔn)差σi的點(diǎn)擴(kuò)散函數(shù)對(duì)應(yīng)的TCLS復(fù)原算法復(fù)原模糊圖像并使用SANS算法去振鈴，其中標(biāo)準(zhǔn)差σi(i＝1～m)為遞增序列，而輸入圖像可以認(rèn)為是標(biāo)準(zhǔn)差為0的TCLS＋SANS復(fù)原算法復(fù)原的圖像，將這m＋1個(gè)圖像按標(biāo)準(zhǔn)差遞增的關(guān)系編號(hào)，編號(hào)分別為1到m＋1，然后對(duì)每幅復(fù)原圖像按像素進(jìn)行區(qū)域圖像質(zhì)量評(píng)價(jià)獲得m＋1幅圖像質(zhì)量圖，將每個(gè)像素位置的圖像質(zhì)量評(píng)價(jià)指標(biāo)最大的圖像編號(hào)輸出，生成最佳像素的復(fù)原圖像編號(hào)圖，由于圖像質(zhì)量評(píng)價(jià)僅使用了區(qū)域信息，會(huì)存在錯(cuò)誤，使用中值濾波處理圖像編號(hào)圖，去除錯(cuò)誤值，最后使用處理后的圖像編號(hào)圖將復(fù)原圖像組合為復(fù)原圖像。本文提出的方法無需對(duì)點(diǎn)擴(kuò)散函數(shù)進(jìn)行估計(jì)，而且利用了TCLS算法可以并行實(shí)現(xiàn)的特點(diǎn)，可以實(shí)現(xiàn)空間變化離焦模糊圖像的快速?gòu)?fù)原。

2.2 空間變化離焦模糊復(fù)原算法實(shí)現(xiàn)

模糊圖像與清晰圖像的區(qū)別主要集中在邊附近，模糊圖像中邊的寬度大，邊不夠清晰，而復(fù)原圖像中邊較為清晰。為了分析圖像邊的情況，使用拉普拉斯算子與圖像卷積獲得拉普拉斯處理的圖像，拉普拉斯算子強(qiáng)調(diào)的是圖像中灰度的突變，而灰度突變一般對(duì)應(yīng)于圖像中清晰邊。拉普拉斯濾波圖像中的“零交叉”位置是圖像的邊，一個(gè)像素為“零交叉”像素時(shí)，在其四鄰域內(nèi)同時(shí)存在正負(fù)像素。隨著復(fù)原濾波器的標(biāo)準(zhǔn)差不斷增大，復(fù)原圖像中的邊逐漸變清晰，邊的邊緣出現(xiàn)明顯的“零交叉”，如果繼續(xù)增大，圖像中出現(xiàn)超出表達(dá)范圍的結(jié)果（壓減到能表達(dá)的最大值），圖像反而質(zhì)量下降。式(8)為拉普拉斯算子與圖像卷積計(jì)算：

式中：fi(x,y)是按照復(fù)原濾波器的標(biāo)準(zhǔn)差由小到大對(duì)復(fù)原后圖像按1～m＋1 進(jìn)行編號(hào)得到復(fù)原圖像組；Li(x,y)為對(duì)應(yīng)的拉普拉斯濾波后的圖像。

檢測(cè)圖像中的“零交叉”可以使用統(tǒng)計(jì)的方法，計(jì)算一個(gè)像素的四鄰域的像素灰度值的絕對(duì)值之和與四鄰域的像素灰度值和的絕對(duì)值的差值，當(dāng)該像素為“零交叉”時(shí)，差值大于0。

Ldi(x,y)為計(jì)算的表征邊強(qiáng)度的指標(biāo)，會(huì)受圖像區(qū)域像素分布影響，為了能夠比較不同標(biāo)準(zhǔn)差的復(fù)原濾波器的復(fù)原結(jié)果，用Ldi(x,y)除以Li(x,y)的絕對(duì)值來去除像素周邊分布的影響，計(jì)算過程為式(10)，其中ε是為了避免除零運(yùn)算而設(shè)的大于零的值，ε的取值較小即可。

圖2 提出的算法流程Fig.2 The proposed algorithm flow

Ldi(x,y)最大的像素認(rèn)為是圖像質(zhì)量最好的像素，使用對(duì)應(yīng)的圖像編號(hào)i作為該像素位置的最優(yōu)像素所在的圖像編號(hào)，生成最優(yōu)像素圖像編號(hào)圖Id(x,y)。

由于Id(x,y)的計(jì)算僅使用了像素周邊區(qū)域的圖像信息，當(dāng)周邊區(qū)域?yàn)槠交瑓^(qū)域，沒有邊的信息時(shí)，計(jì)算結(jié)果會(huì)出現(xiàn)錯(cuò)誤，而相鄰像素的離焦程度具有相關(guān)性，在很小的區(qū)域內(nèi)，最優(yōu)像素圖像編號(hào)不應(yīng)出現(xiàn)劇烈變化，所以對(duì)編號(hào)圖Id(x,y)按一定大小的區(qū)域（如7×7）取中值作為區(qū)域中心像素的編號(hào)。對(duì)于有邊的區(qū)域，計(jì)算結(jié)果正確的像素占多數(shù)，取中值處理能剔除錯(cuò)誤值；對(duì)于大的平滑區(qū)域，由于不同參數(shù)的復(fù)原圖像基本無變化，取中值并不會(huì)導(dǎo)致復(fù)原圖像出現(xiàn)錯(cuò)誤像素。

對(duì)Id(x,y)的每個(gè)像素按區(qū)域取中值獲得最終的最優(yōu)像素圖像編號(hào)圖Idf(x,y)，將Idf(x,y)中每個(gè)像素Idf(i,j)對(duì)應(yīng)的復(fù)原圖像中(i,j)位置的像素值作為最終復(fù)原圖像(i,j)位置的像素值，生成最終復(fù)原圖像。

圖3是使用提出的復(fù)原方法復(fù)原的圖像和過程圖像，使用提出的圖像質(zhì)量評(píng)價(jià)方法能獲得較好的效果。

3 空間變化離焦模糊圖像復(fù)原算法性能分析對(duì)比

3.1 空間變化離焦模糊圖像復(fù)原算法對(duì)比

為了測(cè)試本文提出的算法，使用本文算法、Zhang的算法[4]和Cheong的算法[5]分別對(duì)兩幅空間變化離焦模糊圖像進(jìn)行復(fù)原。Zhang的算法是根據(jù)輸入圖像估計(jì)出點(diǎn)擴(kuò)散函數(shù)分布圖，然后根據(jù)估計(jì)的點(diǎn)擴(kuò)散函數(shù)分割圖像，對(duì)每個(gè)區(qū)域分別使用估計(jì)的點(diǎn)擴(kuò)散函數(shù)進(jìn)行復(fù)原，然后組合為復(fù)原圖像，復(fù)原圖像經(jīng)降噪和去振鈴處理后輸出。Cheong的算法則是根據(jù)圖像的梯度估計(jì)點(diǎn)擴(kuò)散函數(shù)，然后依據(jù)每個(gè)像素估計(jì)的點(diǎn)擴(kuò)散函數(shù)使用TCLS算法復(fù)原像素，然后組合輸出最終圖像。

圖4是對(duì)圖像1 使用不同算法復(fù)原的圖像，圖5是對(duì)圖像2 使用不同算法復(fù)原的圖像，本文提出的算法復(fù)原的圖像最好，圖像的細(xì)節(jié)豐富且對(duì)比度大，Cheong的算法復(fù)原的圖像細(xì)節(jié)丟失多，而Zhang的算法復(fù)原的圖像效果不顯著。Cheong和Zhang的算法的復(fù)原效果差是因?yàn)辄c(diǎn)擴(kuò)散函數(shù)估計(jì)出現(xiàn)了較大誤差，高精度的點(diǎn)擴(kuò)散函數(shù)估計(jì)是很難實(shí)現(xiàn)的。

圖3 復(fù)原圖像和部分過程圖像Fig.3 Restored image and partial process images

表1為使用信息熵和標(biāo)準(zhǔn)差進(jìn)行量化分析復(fù)原圖像的結(jié)果，Zhang的算法的指標(biāo)與輸入圖像最接近，與Zhang的算法復(fù)原的圖像無顯著變化一致，Cheong的算法的指標(biāo)提升大，但實(shí)際圖像細(xì)節(jié)丟失嚴(yán)重，本文提出的算法兩個(gè)指標(biāo)均有提升。

圖4 圖像1 復(fù)原圖像對(duì)比Fig.4 Comparison of the restored image of image 1

圖5 圖像2 復(fù)原圖像對(duì)比Fig.5 Comparison of restored images of image 2

為了比較算法運(yùn)行速度，在同一仿真平臺(tái)和軟件上進(jìn)行算法仿真并統(tǒng)計(jì)算法運(yùn)行時(shí)間。仿真平臺(tái)，CPU：I5-8400，內(nèi)存：16 GB，仿真軟件為Matlab 2017b。3個(gè)算法對(duì)測(cè)試圖像的復(fù)原運(yùn)行的時(shí)間見表2，其中Zhang的算法的復(fù)原運(yùn)行時(shí)間最長(zhǎng)，而本文提出的算法運(yùn)行時(shí)間最短。

表1 復(fù)原圖像質(zhì)量評(píng)價(jià)對(duì)比Table1 Comparison of image quality assessment of restoration algorithms

表2 復(fù)原算法運(yùn)行時(shí)間對(duì)比Table2 Comparison of running time of restoration algorithms s

3.2 本文提出的算法的硬件實(shí)現(xiàn)分析

對(duì)于M×N的輸入圖像，若算法預(yù)設(shè)的復(fù)原濾波器對(duì)應(yīng)的離焦模糊標(biāo)準(zhǔn)差的最大值為σ（σ≤3），則最大的復(fù)原濾波器尺寸為n×n(n＝6×σ＋1)，可計(jì)算得本文提出的算法的時(shí)間復(fù)雜度為O(M×N×n2)。

本文提出的算法復(fù)原過程為多個(gè)預(yù)設(shè)的復(fù)原濾波器分別進(jìn)行卷積計(jì)算，卷積計(jì)算在硬件上可以使用有限脈沖響應(yīng)濾波器實(shí)現(xiàn)，且多個(gè)預(yù)設(shè)的復(fù)原濾波器可以進(jìn)行并行計(jì)算；復(fù)原濾波器的尺寸限定在n×n，卷積計(jì)算硬件優(yōu)化后能夠在3n行之內(nèi)完成，后續(xù)的SANS和圖像質(zhì)量評(píng)價(jià)都使用的是簡(jiǎn)單算法，不會(huì)產(chǎn)生大的延遲，所以本文提出的算法能夠滿足實(shí)時(shí)性要求，但由于存在n×n的卷積計(jì)算，算法的計(jì)算量大，對(duì)硬件資源要求高，后續(xù)研究中可進(jìn)一步優(yōu)化復(fù)原的卷積過程來降低計(jì)算量。硬件實(shí)現(xiàn)時(shí)還需考慮算法對(duì)存儲(chǔ)和傳輸帶寬的要求；本文提出的算法中使用的復(fù)原濾波器是預(yù)設(shè)好的，在硬件實(shí)現(xiàn)中可以直接寫入程序，無需進(jìn)行大量讀寫操作，對(duì)傳輸帶寬要求低；本文提出的算法無需存儲(chǔ)整幅圖像，對(duì)存儲(chǔ)要求低。綜上所述，本文提出的算法能夠進(jìn)行硬件實(shí)現(xiàn)，并能達(dá)到實(shí)時(shí)性要求，但計(jì)算量大，對(duì)硬件資源要求高。

4 結(jié)論

空間變化離焦模糊圖像復(fù)原是提升紅外圖像的一個(gè)重要途徑，本文提出的基于圖像質(zhì)量評(píng)價(jià)的方法能夠快速較好地復(fù)原模糊圖像，算法達(dá)到很好的復(fù)原效果，與基于估計(jì)點(diǎn)擴(kuò)散函數(shù)的方法對(duì)比，該算法能大大降低計(jì)算量，同時(shí)不需要迭代的過程，便于硬件實(shí)現(xiàn)。本文復(fù)原后圖像仍然存在一定的振鈴，部分原因是采用的高斯模型與實(shí)際點(diǎn)擴(kuò)散函數(shù)有誤差，后續(xù)研究可針對(duì)實(shí)際的光學(xué)系統(tǒng)離焦模糊進(jìn)行標(biāo)定，獲得更精準(zhǔn)的模糊模型，進(jìn)而提升復(fù)原圖像的效果。