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

李家強，陳焱博，徐才秀，陳金立，劉 然，朱艷萍

(1. 南京信息工程大學氣象災害預報預警與評估協同創新中心，江蘇南京 210044;2. 南京信息工程大學電子與信息工程學院，江蘇南京 210044;3. 江蘇北斗衛星導航檢測中心有限公司，江蘇南京 210032)

0 引言

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

穿墻雷達常用的成像方法有壓縮感知(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范數最優化問題，利用遠低于Nyquist采樣定理所需采樣數重構原始回波信號成像，雖然CS成像應用前景廣闊，但仍存在加性噪聲和模型誤差干擾求解過程以及圖像結果等問題。DT算法能夠對成像場景的圖像準確重構，但其需要多次迭代運算，計算量巨大，在工程上難以應用。利用目標邊界形狀與接收脈沖延時之間存在逆邊界散射變化的SEABED算法，可以對目標邊界進行清晰成像，但存在障礙物時電磁波的傳播路徑和時延會發生變化，該算法不再滿足逆邊界散射變換，從而不能對目標進行精確成像。近年來，BP算法通過將雷達回波數據投影到成像區域的各個像素點，計算雷達回波在雷達天線和圖像像素之間距離的時延，在時域上進行相干累加來實現高分辨率成像。由于BP算法計算過程簡單易于工程實現，且能精確補償墻體對信號產生的影響，因此被廣泛應用于穿墻成像雷達中，但該算法運算數據量較大，存在冗余現象，從而導致計算復雜度較高，并且對計算機內存需求大，為此，國內外學者提出了一些相關加速BP算法的方法。文獻[10]提出了一種快速因子分解后向投影算法，該方法將子孔徑信號投影到局部極坐標下的成像網格上，通過距離維數的偏移和角度維數的旋轉來校正得到粗糙的子圖像，再將子圖像進行融合得到最后的成像圖，隨后Moon等人提出了一種新的分解BP成像方法[11]，該方法把圖像分割成單獨處理的列，將子圖像連貫地相加。該算法容易實現并行化，由于每一列都可以獨立于其他列形成，因此能夠降低計算復雜度。以上方法都是通過分割子圖像后進行融合，雖然在計算復雜度上有所降低，但同時也降低了成像精度。文獻[12]提出了一種迭代子圖像的BP成像算法，該方法迭代子圖像進行圖像重構降低復雜度，但不適用于大規模的雷達成像。文獻[13]提出了一種快速BP算法，其原理是將成像區域分塊劃分，通過分級相干累加的方式減小BP算法的運算量，該算法雖然有效降低了一定的運算量，但是難以滿足現階段實時成像的需求。基于非均勻快速傅里葉變換的BP成像算法[14]利用BP算法中像素點幅值表達式滿足非均勻傅里葉變換表達式，采用快速傅里葉變換進行計算，降低計算復雜度且保證成像質量，適用于大規模實時成像，但該算法需要重復計算幅值表達式中的系數，對計算機內存需求大。

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

1 穿墻雷達成像模型

圖1 穿墻雷達探測場景模型

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

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

(1)

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

2 基于FGG NUFFT的穿墻雷達快速BP成像

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

2.1 聯合熵值法抑制墻體雜波

本文采用兩組天線陣元分別平行于墻體進行全向掃描得到兩組回波信號。對N個天線陣元回波信號分別進行采樣，記采樣次數為M，兩組回波信號數據分別可組成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)

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

(9)

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

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

(10)

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

(11)

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

2.2 基于FGG NUFFT的快速BP成像

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

(12)

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

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

(13)

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

(14)

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

(15)

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

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

(16)

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

(17)

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

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

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

(18)

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

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

(19)

(20)

E2·E3σ′·E4(σ′)

(21)

式中，-Ksp<σ′

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

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

(22)

其中，

(23)

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

3 實驗結果與分析

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

3.1 穿墻雷達探測場景成像

(a) 基于BP算法成像

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

圖5 成像點歸一化幅值

3.2 計算復雜度分析

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

表1 各類算法的計算復雜度比較

圖6 計算復雜度

3.3 存儲需求

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

表2 各類算法的內存需求比較

圖7 內存需求比較

4 結束語

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