基于FGG NUFFT的穿墻成像雷達(dá)快速BP算法

李家強(qiáng)，陳焱博，徐才秀，陳金立，劉 然，朱艷萍

(1. 南京信息工程大學(xué)氣象災(zāi)害預(yù)報預(yù)警與評估協(xié)同創(chuàng)新中心，江蘇南京 210044;2. 南京信息工程大學(xué)電子與信息工程學(xué)院，江蘇南京 210044;3. 江蘇北斗衛(wèi)星導(dǎo)航檢測中心有限公司，江蘇南京 210032)

0 引言

穿墻成像雷達(dá)(Through-Wall Imaging Radar, TWIR)利用電磁波良好的穿透特性，能夠?qū)ㄖ?、堡壘、草叢等隱蔽物后目標(biāo)進(jìn)行探測，其硬件成本低、質(zhì)量輕、體積小、便攜性好，并且能夠進(jìn)行實時目標(biāo)定位、分類與識別，因此被廣泛用于軍用與民用領(lǐng)域，諸如執(zhí)法行動、反恐斗爭、刑事偵查、搶險救援等[1-5]。

穿墻雷達(dá)常用的成像方法有壓縮感知(Compression Perception，CS)成像算法[6]、逆散射層析(Diffraction Tomography，DT)算法[7]、邊界估計(Shape Estimation Algorithm Based on BST and Extraction of Directly Scattered，SEABED)算法[8]、后向投影(Back Projection，BP)算法[9]等。CS成像對稀疏或可壓縮的回波信號，通過求解帶約束的l1范數(shù)最優(yōu)化問題，利用遠(yuǎn)低于Nyquist采樣定理所需采樣數(shù)重構(gòu)原始回波信號成像，雖然CS成像應(yīng)用前景廣闊，但仍存在加性噪聲和模型誤差干擾求解過程以及圖像結(jié)果等問題。DT算法能夠?qū)Τ上駡鼍暗膱D像準(zhǔn)確重構(gòu)，但其需要多次迭代運算，計算量巨大，在工程上難以應(yīng)用。利用目標(biāo)邊界形狀與接收脈沖延時之間存在逆邊界散射變化的SEABED算法，可以對目標(biāo)邊界進(jìn)行清晰成像，但存在障礙物時電磁波的傳播路徑和時延會發(fā)生變化，該算法不再滿足逆邊界散射變換，從而不能對目標(biāo)進(jìn)行精確成像。近年來，BP算法通過將雷達(dá)回波數(shù)據(jù)投影到成像區(qū)域的各個像素點，計算雷達(dá)回波在雷達(dá)天線和圖像像素之間距離的時延，在時域上進(jìn)行相干累加來實現(xiàn)高分辨率成像。由于BP算法計算過程簡單易于工程實現(xiàn)，且能精確補(bǔ)償墻體對信號產(chǎn)生的影響，因此被廣泛應(yīng)用于穿墻成像雷達(dá)中，但該算法運算數(shù)據(jù)量較大，存在冗余現(xiàn)象，從而導(dǎo)致計算復(fù)雜度較高，并且對計算機(jī)內(nèi)存需求大，為此，國內(nèi)外學(xué)者提出了一些相關(guān)加速BP算法的方法。文獻(xiàn)[10]提出了一種快速因子分解后向投影算法，該方法將子孔徑信號投影到局部極坐標(biāo)下的成像網(wǎng)格上，通過距離維數(shù)的偏移和角度維數(shù)的旋轉(zhuǎn)來校正得到粗糙的子圖像，再將子圖像進(jìn)行融合得到最后的成像圖，隨后Moon等人提出了一種新的分解BP成像方法[11]，該方法把圖像分割成單獨處理的列，將子圖像連貫地相加。該算法容易實現(xiàn)并行化，由于每一列都可以獨立于其他列形成，因此能夠降低計算復(fù)雜度。以上方法都是通過分割子圖像后進(jìn)行融合，雖然在計算復(fù)雜度上有所降低，但同時也降低了成像精度。文獻(xiàn)[12]提出了一種迭代子圖像的BP成像算法，該方法迭代子圖像進(jìn)行圖像重構(gòu)降低復(fù)雜度，但不適用于大規(guī)模的雷達(dá)成像。文獻(xiàn)[13]提出了一種快速BP算法，其原理是將成像區(qū)域分塊劃分，通過分級相干累加的方式減小BP算法的運算量，該算法雖然有效降低了一定的運算量，但是難以滿足現(xiàn)階段實時成像的需求。基于非均勻快速傅里葉變換的BP成像算法[14]利用BP算法中像素點幅值表達(dá)式滿足非均勻傅里葉變換表達(dá)式，采用快速傅里葉變換進(jìn)行計算，降低計算復(fù)雜度且保證成像質(zhì)量，適用于大規(guī)模實時成像，但該算法需要重復(fù)計算幅值表達(dá)式中的系數(shù)，對計算機(jī)內(nèi)存需求大。

針對上述問題，本文提出了一種基于快速高斯網(wǎng)格化的非均勻的快速傅里葉變換加速BP算法。該算法將BP算法中像素點幅值表達(dá)式與高斯核函數(shù)反卷積消除高斯平滑的影響，其次對均勻數(shù)據(jù)進(jìn)行快速傅里葉變換，最后對得到的數(shù)據(jù)進(jìn)行卷積運算平滑輸出數(shù)據(jù)。仿真結(jié)果表明，所提算法能夠避免重復(fù)計算和存儲BP算法中像素點幅值表達(dá)式中的系數(shù)，從而快速計算像素點幅值矩陣，并且在不改變成像質(zhì)量的同時，大大降低計算復(fù)雜度和內(nèi)存需求，適用于大規(guī)模成像和實時處理。

1 穿墻雷達(dá)成像模型

圖1 穿墻雷達(dá)探測場景模型

穿墻成像雷達(dá)模型如圖1所示：前墻與后墻的墻體厚度為d，相對介電常數(shù)為εw，探測目標(biāo)為半徑r的理想電導(dǎo)體，置于兩面墻之間，圓心與前墻的垂直距離為δ。設(shè)置平行于前墻體的兩組間距極小的收發(fā)共置天線，其距離墻體分別為h1和h2，兩組天線分別先后發(fā)射電磁波。電磁波首先穿透前墻體，經(jīng)過墻體另一側(cè)空氣傳播，到達(dá)目標(biāo)物體后按原路散射回天線，兩組天線具有相同的回波信號模型：

e(t)=et(t)+ew(t)+ea(t)

(1)

式中，et(t)為所需探測的目標(biāo)信號，ew(t)為前墻體與后墻體的反射信號，假定天線為理想天線，那么噪聲信號ea(t)=0。

2 基于FGG NUFFT的穿墻雷達(dá)快速BP成像

穿墻雷達(dá)探測過程中受墻體引起的強(qiáng)烈雜波信號干擾，往往目標(biāo)信號會被淹沒，無法進(jìn)行目標(biāo)檢測及成像，因此成像之前必須通過算法消除墻體雜波信號的影響。

2.1 聯(lián)合熵值法抑制墻體雜波

本文采用兩組天線陣元分別平行于墻體進(jìn)行全向掃描得到兩組回波信號。對N個天線陣元回波信號分別進(jìn)行采樣，記采樣次數(shù)為M，兩組回波信號數(shù)據(jù)分別可組成M×N維矩陣e1，e2：

e1=[y1,y2,…,yN]

(2)

(3)

將式(2)與式(3)相加，能夠得到

e=e1+e2=[z1,z2,…，zN]

(4)

計算矩陣e1，e2中每個采樣點值的概率：

n=1,2,…，N，m=1,2,…，M

(5)

n=1,2,…，N，m=1,2,…，M

(6)

(7)

(8)

根據(jù)聯(lián)合熵值定理[15]，統(tǒng)計獨立信源的聯(lián)合熵值等于各信源熵之和，這里兩組回波信號前后分別獲得，可以視作相互獨立，因此兩者的聯(lián)合熵值表示如下：

(9)

由于天線陣元不同掃描位置的變化，目標(biāo)回波信號的變化強(qiáng)度較大，而墻體雜波的變化強(qiáng)度較小，因此目標(biāo)回波信號熵值較小，墻體雜波熵值較大[16]。根據(jù)聯(lián)合熵值定理[15]，統(tǒng)計獨立信源的聯(lián)合熵值等于它們分別熵值之和，且大于其任何一個信源的熵值，即墻體雜波的聯(lián)合熵值大于目標(biāo)回波的聯(lián)合熵值，可通過設(shè)置適當(dāng)?shù)拈T限值Z*能夠濾除墻體雜波的影響，對回波數(shù)據(jù)e中每個元素進(jìn)行處理：

Zn(m)=Z*·zn(m)

(10)

式中，Zn(m)為去除墻體雜波后所需的目標(biāo)回波信號，而Z*定義為

(11)

式中，log(N)為聯(lián)合熵值H的最大值，β為門限可調(diào)節(jié)因子，其調(diào)節(jié)范圍為(0,2)。根據(jù)最大離散熵值定理[15]，在離散信源等概率出現(xiàn)的情況下，信源的平均不確定性為最大，熵值將達(dá)到最大，即聯(lián)合熵值取得最大值。

2.2 基于FGG NUFFT的快速BP成像

圖2為穿墻雷達(dá)成像區(qū)域像素點網(wǎng)格化分，將整個探測場景劃分為P個像素網(wǎng)格點，l1，l2，l3分別表示發(fā)射電磁波天線陣元與前墻體的斜距、電磁波在墻體傳播的距離以及電磁波穿過墻體后到達(dá)目標(biāo)的距離，則穿墻雷達(dá)BP成像可表示為

(12)

式中，I(xp,yp)為第p(p=1,2,…,P-1)個像素網(wǎng)格點的復(fù)幅度值，fm=f0+mΔf為第m(m=1,2,…，M)個工作頻點，f0為發(fā)射信號的起始頻率，Δf為頻率間隔，τp為成像像素網(wǎng)格點p和第n(n=1,2,…，N)根天線之間的雙程時延：

τp=2(l1+l3)/c+2l2/v

(13)

NUFFT能夠直接對非均勻數(shù)據(jù)進(jìn)行快速處理，即可將均勻采樣數(shù)據(jù)變換到非均勻采樣數(shù)據(jù)[17]。對于均勻采樣序列F(k)∈C，k=0,1,…,K-1,K∈，則非均勻采樣序列的離散傅里葉變換定義為

(14)

式中xj∈[0,2π]。gτ(x)是[0,2π]上的一維周期高斯核函數(shù)，其表達(dá)式如下：

(15)

式中，τ為高斯核函數(shù)參數(shù)，決定高斯核函數(shù)的指數(shù)衰減速率。

對高斯核函數(shù)gτ(x)進(jìn)行傅里葉變換得到Gτ(k)，令Gτ(k)與均勻數(shù)據(jù)F(k)反卷積消除高斯平滑的影響，得到輔助函數(shù)F-τ(k)：

(16)

進(jìn)一步對式(16)進(jìn)行離散傅里葉反變換得到f-τ(x)：

(17)

式中，Kr=R×K，Kr為過采樣網(wǎng)格數(shù)，R為過采樣系數(shù)。由于f-τ(x)滿足均勻性，分布在均勻的網(wǎng)格點上，因此，可以用快速傅里葉變換對其進(jìn)行計算。

對f-τ(x)進(jìn)行卷積實現(xiàn)平滑輸出，得到期望值f(xj)：

f(xj)=f-τ*gτ(xj)=

(18)

在求解式(18)時，每次計算非均勻采樣點xj時都要遍歷所有的均勻網(wǎng)格點，因此計算量巨大，由于高斯函數(shù)的指數(shù)衰減特性，對遠(yuǎn)離xj的網(wǎng)格點忽略不計，因此可設(shè)置網(wǎng)格擴(kuò)散范圍，僅考慮xj附近的Ksp個點，Ksp為高斯核函數(shù)單邊可延伸覆蓋的網(wǎng)格點數(shù)量(Ksp=6表示單精度，Ksp=12表示雙精度)。

式(18)中f-τ、gτ可表示為

(19)

(20)

E2·E3σ′·E4(σ′)

(21)

式中，-Ksp<σ′

由于目標(biāo)回波數(shù)據(jù)Zn(m)是均勻數(shù)據(jù)，時延τp對于單個像素點的幅值I(xp,yp)是非均勻的，無法直接用快速傅里葉變換進(jìn)行計算，而本文所提FGG NUFFT算法適應(yīng)非均勻性，能夠?qū)Ψ蔷鶆驍?shù)據(jù)進(jìn)行快速傅里葉變換。該算法的主要思想就是將均勻樣本與高斯函數(shù)反卷積消除高斯平滑的影響，其次進(jìn)行快速傅里葉變換，最后將樣本卷積高斯函數(shù)實現(xiàn)平滑輸出。因此，運用FGG NUFFT算法對像素點的幅值I(xp,yp)進(jìn)行快速計算。

式(12)中BP成像算法可變形為

(22)

其中，

(23)

式(23)能夠通過FGG NUFFT來計算，式中每個像素點的幅值In(xp,yp)對應(yīng)式(14)中的非均勻采樣點f(xj)，均勻目標(biāo)回波數(shù)據(jù)Zn(m)對應(yīng)式(14)中的均勻采樣點F(k)，頻率間隔與雙程時延的乘積2πΔfτp對應(yīng)式(14)中的非均勻數(shù)據(jù)xj。

3 實驗結(jié)果與分析

為了驗證所提算法的有效性，本文利用基于FDTD的正演數(shù)值模擬軟件GprMax2D/3D模擬穿墻雷達(dá)場景。模型設(shè)置如圖1所示：天線陣元距離墻體分別為h1=0.05 m，h2=0.04 m；沿水平方向等間距掃描N=50次，掃描范圍為0.1～2.1 m；天線發(fā)射電磁波的時間窗為24 ns。墻體為均勻介質(zhì)，其厚度d=0.2 m，相對介電常數(shù)εw=6.4；理想電導(dǎo)體目標(biāo)的半徑r=0.1 m，其圓心距離墻體δ=1.0 m。圖3為經(jīng)過聯(lián)合熵值抑制墻體雜波后得到的目標(biāo)時域回波圖。

3.1 穿墻雷達(dá)探測場景成像

(a) 基于BP算法成像

分別運用BP成像算法，NUFFT BP成像算法[17]和本文所提算法對穿墻雷達(dá)探測區(qū)域進(jìn)行成像。選取電磁波頻段范圍為1～2 GHz，頻率間隔Δf=0.49 MHz，頻率點M=2 036。成像區(qū)域設(shè)置橫向距離2.2 m，縱向距離2.1 m，按照橫向的和縱向劃分279×293個像素點。圖4為3種算法的成像結(jié)果，其中圖4(a)為基于BP算法成像、圖4(b)為基于NUFFT BP算法成像、圖4(c)為基于FGG NUFFT BP算法成像。為了比較3種算法的成像分辨率，對像素點矩陣進(jìn)行歸一化幅值，圖5為3種算法的成像像素點歸一化幅值曲線圖。取成像最大幅值下降-3 dB處的像素點寬度，BP成像為3.8個像素點，NUFFT BP成像為5.3個像素點，本文算法成像為4.9個像素點，3種算法的寬度相差較小，均能對目標(biāo)區(qū)域進(jìn)行精確成像，且本文算法保持較好效果的同時優(yōu)于NUFFT BP成像。

圖5 成像點歸一化幅值

3.2 計算復(fù)雜度分析

分別對BP成像、NUFFT BP成像和本文提出算法成像的計算復(fù)雜度進(jìn)行簡要的分析。假設(shè)忽略天線根數(shù)的影響，取N=1，由于天線的采樣點遠(yuǎn)遠(yuǎn)小于像素點個數(shù)，即M?P。BP算法的計算復(fù)雜度為C0=MP，NUFFT BP的計算復(fù)雜度為C1=MlogM+|logε|P≈|logε|P，其中ε為理想計算精度[18-19]。本文設(shè)置高斯核函數(shù)單邊可延伸覆蓋6個網(wǎng)格點，即運算精度Ksp=6；過采樣系數(shù)R=2，取值細(xì)節(jié)見文獻(xiàn)[20]。本文所提的算法中E，E1，E2，E3，E4能預(yù)先被計算和存儲，不需要重復(fù)計算，因此本文算法的計算復(fù)雜度為C2=(Ksp+1)M+(KrlogKr)+P≈P。圖6為3種算法的計算復(fù)雜度與像素點網(wǎng)格點變化的曲線圖，通過比較可以發(fā)現(xiàn)在同等條件下，對于同樣的像素點本文算法計算復(fù)雜度要遠(yuǎn)低于其他兩種算法計算復(fù)雜度，各類算法復(fù)雜度見表1。

表1 各類算法的計算復(fù)雜度比較

圖6 計算復(fù)雜度

3.3 存儲需求

分別對BP成像、NUFFT BP成像和本文提出算法成像的內(nèi)存需求進(jìn)行簡要的分析。實部和虛部都用浮點數(shù)表示，每個實部的內(nèi)存需求為4 Byte，每個虛部的內(nèi)存需求為8 Byte。BP的像素點幅值矩陣I∈CM×P，因此直接計算BP成像的內(nèi)存需求為S0=8MP。 NUFFT BP成像的內(nèi)存需求為S1=8(MlogM+|logε|P)≈8|logε|P。本文算法成像只需要存儲系數(shù)E∈RP，E1∈C2Ksp×M，E2∈R2Ksp×P，E3∈R2Ksp×P，E4∈C2Ksp，則S2=4P(1+4Ksp)+8Ksp(1+2M)≈100P。圖7為3種算法的內(nèi)存需求與像素點網(wǎng)格點變化的曲線圖，通過比較可以發(fā)現(xiàn)在同等條件下，對于同樣的像素點本文算法內(nèi)存需求要低于其他兩種算法內(nèi)存需求，各類算法內(nèi)存需求見表2。

表2 各類算法的內(nèi)存需求比較

圖7 內(nèi)存需求比較

4 結(jié)束語

傳統(tǒng)的BP成像算法需要計算每個回波數(shù)據(jù)與圖像像素之間雙程并進(jìn)行相干累加，因此存在計算復(fù)雜度高，內(nèi)存需求量巨大，需要重復(fù)計算等問題。本文在對墻體雜波進(jìn)行聯(lián)合熵值法抑制的基礎(chǔ)上，將FGG NUFFT算法運用到傳統(tǒng)的BP成像中，首先將BP算法中像素點幅值表達(dá)式與高斯核函數(shù)反卷積消除高斯平滑的影響，其次對均勻數(shù)據(jù)進(jìn)行快速傅里葉變換，最后對得到的數(shù)據(jù)進(jìn)行卷積運算實現(xiàn)對數(shù)據(jù)均勻平滑輸出。仿真實驗結(jié)果表明，本文的方法在保證成像質(zhì)量的情況下，解決了計算復(fù)雜度和內(nèi)存需求隨著工作頻點、天線數(shù)量、像素點的增加而增加的問題，為穿墻雷達(dá)大規(guī)模成像與實時成像提供了新的思路與解決途徑。