馬昭陽(yáng)，鄧桂林，曾 旸，鄧 彬，楊 琪，王宏強(qiáng)

(國(guó)防科技大學(xué) 電子科學(xué)學(xué)院, 湖南 長(zhǎng)沙 410073)

毫米波是指波長(zhǎng)為1～10 mm、頻率為30～300 GHz的電磁波。與頻率較低的微波相比較，毫米波可以獲得更高的分辨率和更強(qiáng)的穿透性；與X射線相比較，毫米波波段的光子能量低，不會(huì)損害探測(cè)目標(biāo)。因此，三維毫米波全息成像技術(shù)在無(wú)損檢測(cè)中極具優(yōu)勢(shì)[1-2]。

毫米波全息成像系統(tǒng)主要包括基于平面陣列的成像系統(tǒng)和基于柱面陣列的成像系統(tǒng)。其中基于平面陣列的成像系統(tǒng)發(fā)展較成熟，已經(jīng)被成功應(yīng)用于各種領(lǐng)域，例如安全檢查、無(wú)損檢測(cè)、穿墻成像等[2-3]。由于收發(fā)天線均勻分布在平面中的二維單發(fā)單收(Two Dimensional Single Input Single Output, 2D-SISO)系統(tǒng)造價(jià)昂貴、難以普及，將其中一維用掃描方式替代實(shí)現(xiàn)的一維單發(fā)單收一維掃描(One Dimensional Single Input Single Output One Dimensional Scanning, 1D-SISO-1D-scanning)系統(tǒng)得到了廣泛的研究和應(yīng)用[4-6]。

另外，研究者們開始探索通過(guò)設(shè)計(jì)非均勻陣列提高成像效率[7-9]。在這種情況下，回波信號(hào)的非均勻性使得傳統(tǒng)的基于傅里葉變換的成像算法不再直接適用。因此，研究出一種適用于非均勻平面陣列的無(wú)損檢測(cè)三維成像算法是非常必要的。后向投影(Back Projection, BP)算法是雷達(dá)成像的標(biāo)準(zhǔn)算法，可以適用于各種成像場(chǎng)景。但是BP算法極高的運(yùn)算復(fù)雜度使其在陣元數(shù)目較多的情況下需要耗費(fèi)巨大的計(jì)算資源，難以滿足成像的時(shí)效性需求。距離徙動(dòng)算法(Range Migration Algorithm, RMA)是一種能兼顧成像質(zhì)量和速度的雷達(dá)成像算法，也是目前在實(shí)際場(chǎng)景中被應(yīng)用最廣泛的一種算法。然而RMA中處理非均勻數(shù)據(jù)的插值步驟往往會(huì)帶來(lái)運(yùn)算量過(guò)大或者精度不夠的問(wèn)題。非均勻快速傅里葉變換(NonUniform Fast Fourier Transform, NUFFT)可以在非均勻空間實(shí)現(xiàn)快速傅里葉變換，因此自被提出就展現(xiàn)出了其在雷達(dá)成像算法方面的巨大潛力[10-12]。

綜上所述，本文針對(duì)毫米波非均勻平面陣列，結(jié)合RMA和NUFFT的思想，提出了一種適用于無(wú)損檢測(cè)的三維成像算法。通過(guò)實(shí)驗(yàn)室搭建的毫米波非均勻1D-SISO-1D-scanning成像系統(tǒng)驗(yàn)證了所提算法的有效性。

1 成像模型和算法

毫米波非均勻平面陣列體制下的無(wú)損檢測(cè)成像幾何關(guān)系如圖1所示。陣面上的天線陣元位置為rR=(xR,yR,zR)，目標(biāo)中的缺陷位置為r=(x,y,z)。定義空氣部分為區(qū)域A，介電常數(shù)為εA；目標(biāo)部分為區(qū)域B，介電常數(shù)為εB；目標(biāo)中的缺陷部分為區(qū)域C，介電常數(shù)為εC。

圖1 成像幾何模型Fig.1 Imaging geometric model

缺陷區(qū)域的散射場(chǎng)可以近似寫為：

(1)

其中：k=2πf/c為空間波數(shù)，f為發(fā)射信號(hào)的頻率，c為光速；Et(r,rR,k)=jkη0G(r,rR,k)為入射場(chǎng)，η0為自由空間的波阻抗；o(r)=εC(r)-εB(r)為反射率函數(shù)；G(rR,r,k)為半空間格林函數(shù)。

根據(jù)格林函數(shù)的性質(zhì)，可將式(1)改寫為：

(2)

其中，半空間格林函數(shù)的表達(dá)式為：

jkz(zR-z)]F(kx,kz,k)dkxdkz

(3)

(4)

近場(chǎng)雷達(dá)成像主要受相位項(xiàng)的影響，因此將式(3)代入式(2)，并忽略幅度項(xiàng)可得：

j(kAyyR-kByy)+jkz(zR-z)]dkxdkzdr

(5)

對(duì)式(5)進(jìn)行關(guān)于xR和zR的傅里葉變換可得：

j(kAyyR-kByy)-jkzz]dr

(6)

對(duì)式(6)進(jìn)行相位補(bǔ)償可得：

exp(-jkByy-jkxx-jkzz)dr

(7)

則待求解的目標(biāo)函數(shù)可以寫為：

exp(jkByy+jkxx+jkzz)dkxdkBydkz

(8)

在實(shí)際成像過(guò)程中，首先考慮到平面陣列是非均勻分布的，故引入NUFFT實(shí)現(xiàn)非均勻數(shù)據(jù)的傅里葉變換：

Er(kx,kz,k)=NUFFTxR,zR[Er(rR,k)]exp(jkAyyR)

(9)

其中，NUFFTxR,zR[·]表示進(jìn)行關(guān)于xR和zR的非均勻快速傅里葉變換。

另外考慮到kx和kz通常均勻分布，其逆傅里葉變換可以直接通過(guò)逆快速傅里葉變換(Inverse Fast Fourier Transform, IFFT)實(shí)現(xiàn)。 但kBy為非均勻頻域數(shù)據(jù)，RMA處理時(shí)會(huì)進(jìn)行插值操作，所提算法通過(guò)引入非均勻逆快速傅里葉變換(NonUniform Inverse Fast Fourier Transform, NUIFFT)取代插值和IFFT進(jìn)行處理：

exp(jkxx+jkzz)dkxdkz]

(10)

其中，NUIFFTkBy[·]表示進(jìn)行關(guān)于kBy的非均勻逆快速傅里葉變換。

因此所提算法的成像公式可總結(jié)為：

exp(jkAyyR)]]

(11)

其中，IFFTkx,kz[·]表示進(jìn)行關(guān)于kx和kz的逆快速傅里葉變換。

2 成像結(jié)果分析

介質(zhì)目標(biāo)設(shè)置為相對(duì)介電常數(shù)約為2.008的聚四氟乙烯。設(shè)置的非均勻平面陣列排布方式如圖2所示，在陣元數(shù)為101×101=10 201的均勻平面陣列中通過(guò)隨機(jī)抽樣的方式隨機(jī)抽取掉2 624個(gè)陣元，形成共7 577個(gè)陣元的非均勻平面陣列。發(fā)射信號(hào)的中心頻率為37.5 GHz，帶寬為12 GHz。

圖2 非均勻平面陣列構(gòu)型Fig.2 Nonuniform planar array topology

2.1 仿真結(jié)果

首先通過(guò)仿真驗(yàn)證所提算法的有效性并進(jìn)行成像性能的評(píng)估。目標(biāo)設(shè)置為一個(gè)坐標(biāo)為(0,0.05,0)m的理想真空點(diǎn)缺陷。分別利用BP算法和所提算法進(jìn)行成像，得到的結(jié)果如圖3所示。

由圖3可以看到，兩種算法均可以對(duì)介質(zhì)中的點(diǎn)目標(biāo)缺陷進(jìn)行準(zhǔn)確成像。為了定量分析兩種算法的性能，以沖激響應(yīng)寬度(Impulse Response Width, IRW)和峰值旁瓣比(Peak Side Lobe Ratio, PSLR)為精度評(píng)價(jià)指標(biāo)、運(yùn)行時(shí)間為速度評(píng)價(jià)指標(biāo)進(jìn)行對(duì)比，在相同計(jì)算機(jī)運(yùn)行環(huán)境下得到的結(jié)果如表1所示。

(a) BP算法x-y平面投影(a) x-y plane projection utilizing BP algorithm

(b) 所提算法x-y平面投影(b) x-y plane projection utilizing the proposed algorithm

(c) BP算法x-z平面投影(c) x-z plane projection utilizing BP algorithm

(d) 所提算法x-z平面投影(d) x-z plane projection utilizing the proposed algorithm

(e) BP算法y-z平面投影(e) y-z plane projection utilizing BP algorithm

(f) 所提算法y-z平面投影(f) y-z plane projection utilizing the proposed algorithm圖3 仿真成像結(jié)果Fig.3 Simulation imaging results

通過(guò)表1的對(duì)比可以看到，在三個(gè)坐標(biāo)軸方向上BP算法的IRW都略微小于所提算法的IRW，兩者相差不大，說(shuō)明兩種算法的分辨率幾乎一致。而在三個(gè)坐標(biāo)軸方向上所提算法的PSLR都小于BP算法的PSLR，說(shuō)明所提算法的旁瓣抑制效果更好?？偟膩?lái)說(shuō)，所提算法的成像精度能夠與BP算法基本相當(dāng)。在成像速度方面，所提算法遠(yuǎn)遠(yuǎn)快于BP算法。仿真結(jié)果有效證明了所提算法的成像效率優(yōu)勢(shì)。

2.2 實(shí)驗(yàn)結(jié)果

在實(shí)驗(yàn)室搭建的毫米波非均勻1D-SISO-1D-scanning成像系統(tǒng)中，金屬“檸檬片”夾在前后厚度分別為2 cm和0.5 cm的聚四氟乙烯塊狀材料之間。目標(biāo)通過(guò)泡沫塊固定在毫米波非均勻陣列前方。

分別利用BP算法和所提算法進(jìn)行成像，得到的結(jié)果如圖4所示。在實(shí)驗(yàn)的三維成像結(jié)果中，介質(zhì)目標(biāo)中的金屬“檸檬片”清晰可見，且可以看到聚四氟乙烯塊狀材料的輪廓。兩種算法的成像質(zhì)量基本相同，但是BP算法得到的成像結(jié)果中旁瓣更多，這也與前文對(duì)仿真結(jié)果的分析一致。

由表1中的運(yùn)行時(shí)間可知，此時(shí)BP算法耗費(fèi)的計(jì)算資源仍然是巨大的，而所提算法依舊能夠在較短時(shí)間內(nèi)完成成像。

表1 性能比較

(a) BP算法x-y平面投影(a) x-y plane projection utilizing BP algorithm

(b) 所提算法x-y平面投影(b) x-y plane projection utilizing the proposed algorithm

(c) BP算法x-z平面投影(c) x-z plane projection utilizing BP algorithm

(d) 所提算法x-z平面投影(d) x-z plane projection utilizing the proposed algorithm

(e) BP算法y-z平面投影(e) y-z plane projection utilizing BP algorithm

(f) 所提算法y-z平面投影(f) y-z plane projection utilizing the proposed algorithm圖4 實(shí)驗(yàn)成像結(jié)果Fig.4 Experimental imaging results

3 結(jié)論

本文提出了一種毫米波非均勻平面陣列無(wú)損檢測(cè)三維成像算法，解決了傳統(tǒng)基于傅里葉變換的成像算法不能適用于非均勻信號(hào)和算法中插值操作難以兼顧精度與速度的問(wèn)題。仿真中，利用理想點(diǎn)缺陷目標(biāo)定量分析了兩種算法的成像性能。同時(shí)利用搭建的毫米波非均勻1D-SISO-1D-scanning成像系統(tǒng)進(jìn)行了實(shí)驗(yàn)。結(jié)果表明，相較于BP算法，所提算法能夠?qū)崿F(xiàn)相當(dāng)?shù)某上窬龋瑫r(shí)在成像速度方面具有明顯優(yōu)勢(shì)。本文工作可為未來(lái)毫米波成像技術(shù)在無(wú)損檢測(cè)領(lǐng)域的發(fā)展提供思路和支持。