基于隨機陣列的NUFFT BP三維快速成像

向 高，張曉玲，吳宗亮

（電子科技大學(xué)電子工程學(xué)院，四川成都611731）

基于隨機陣列的NUFFT BP三維快速成像

向 高，張曉玲，吳宗亮

（電子科技大學(xué)電子工程學(xué)院，四川成都611731）

傳統(tǒng)的二維實陣列用于三維微波成像時，需要相對較高的脈沖重復(fù)頻率才能切換發(fā)射陣元。針對這種情況，利用隨機陣列來實現(xiàn)三維成像，目的是有效降低脈沖重復(fù)頻率。進(jìn)一步提出利用等效相位中心原理來實現(xiàn)任意隨機陣列布局的可行性方案。然后，基于隨機陣列，發(fā)現(xiàn)使用非均勻快速傅里葉變換（nonuniform fast Fourier transform，NUFFT）來實現(xiàn)后向投影（back projection，BP）算法的插值過程，能夠用較少的運算量實現(xiàn)BP算法的高精度插值。最后，結(jié)合并行處理技術(shù)來實現(xiàn)NUFFT BP算法，結(jié)果使得BP算法的執(zhí)行效率得到顯著提高。

隨機陣列；非均勻快速傅里葉變換；并行處理；后向投影

0 引 言

傳統(tǒng)的二維實陣列用于三維微波成像時，主要有均勻滿陣元方式和天線單元曲線運動合成面陣［1］兩種方式。對于滿陣元方式，除了成本高以外，當(dāng)要求所有陣元都要收／發(fā)信號的時候，需要很高的脈沖重復(fù)頻率（pulse repetition frequency，PRF）來實現(xiàn)切換陣元。同時，運動合成面陣容易受到陣元機械運動的影響，例如：存在平臺晃動，需要較長的運動合成時間等。

為了解決上述問題，本文提出使用隨機二維陣列進(jìn)行三維微波成像。盡管文獻(xiàn)［2］等已經(jīng)提出了隨機陣列成像的概念，但是并未給出合理的隨機陣列布局的生成方案。文獻(xiàn)［3- 4］雖然都利用等效相位中心（phase center approximation，PCA）原理來等效均勻陣列，然而該等效陣列僅僅只是一維的。文獻(xiàn)［5］將PCA原理進(jìn)一步推廣到了二維陣列，實現(xiàn)了多輸入多輸出（multiple input multiple output，MIMO）陣列的等效均勻二維陣列。本文正是以這種等效均勻二維陣列為基礎(chǔ)，實現(xiàn)了任意構(gòu)型的隨機二維陣列。本文以“十”字形布陣為基礎(chǔ)，分析了隨機陣列的隨機陣元數(shù)目的選擇。由于隨機陣列其陣元位置的不確定性，所以無法使用距離 多普勒［67］（range-Doppler，RD）等傳統(tǒng)算法進(jìn)行聚焦，因此提出使用Type-Ⅱ型非均勻快速傅里葉變換（nonuniform fast Fourier transform，NUFFT）［8］插值的后向投影（back projection，BP）算法進(jìn)行聚焦，并利用并行計算架構(gòu)［9－10］（compute unified device architecture，CUDA）實現(xiàn)了距離歷史和NUFFT插值的并行運算。

1 隨機二維陣列

隨機二維陣列其陣元在陣列平面內(nèi)的均勻網(wǎng)格上隨機分布，如圖1所示。當(dāng)這些隨機陣元不共線時，即可以產(chǎn)生空間二維角分辨率［2］，再結(jié)合脈沖壓縮技術(shù)，則能夠?qū)崿F(xiàn)三維微波成像。如何合理地選擇隨機二維陣列的構(gòu)成方案，這是主要考慮的問題。因此，本節(jié)提出了一種合理可行的隨機二維陣列的構(gòu)成方案，并給出了合理的隨機陣元數(shù)目的選擇方法。

圖1 隨機二維陣列成像的幾何構(gòu)型

1.1 隨機二維陣列的布局設(shè)計

如果采用直接隨機布置陣元的方式來實現(xiàn)隨機二維陣列，不僅成本高，難于實現(xiàn)，而且無法靈活改變陣元位置的分布，結(jié)果使得其實用性大大降低。因此，可以考慮利用虛擬陣元的方式來實現(xiàn)隨機二維陣列。圖2展示了一種基于二維MIMO陣列來實現(xiàn)虛擬隨機二維陣列的布局方案［5］，該方案基于PCA原理來實現(xiàn)虛擬二維均勻陣列，此時只要激活合適的收發(fā)陣元，就可以很容易地獲得任意構(gòu)型的隨機二維陣列。簡要敘述PCA原理：可以用一個虛擬的天線單元，等效出位置上分置的發(fā)射天線和接收天線兩者的功能，該虛擬天線單元同時具有收／發(fā)功能，并且位于分置的收發(fā)天線兩者連線的中心。

圖2 MIMO二維面陣

簡單分析圖2中是如何實現(xiàn)虛擬二維均勻陣列。圖2使用了端發(fā)“密疏密”布陣模式［4］，能夠很好地利用陣列長度。在陣列平面的Ⅰ、Ⅱ、Ⅲ和Ⅳ區(qū)域布置發(fā)射陣元，用符號T表示，發(fā)射陣元沿x軸間距為dxT，沿y軸間距為dyT，發(fā)射陣元沿x軸和y軸數(shù)目記為MaT，要求其為偶數(shù)，那么總的發(fā)射陣元數(shù)目為MT＝M2aT；在陣列平面中心區(qū)域布置接收陣元，用符號R表示，接收陣元沿x軸間距為dxR，沿y軸間距為dyR，接收陣元沿x軸和y軸數(shù)目記為MaR，那么總的接收陣元數(shù)目為MR＝M2aR；最鄰近的發(fā)射陣元和接收陣元沿x軸和y軸的間距分別為dxTR和dyTR。如果上述陣元間距dxT、dyT、dxR、dyR、dxTR和dyTR滿足約束關(guān)系：dxT＝2da、dyT＝2da、dxR＝daMa／2、dyR＝daMa／2、dxTR＝da和dyTR＝da，那么就可以獲得一個虛擬的二維均勻陣列，此時該虛擬陣列沿x軸和y軸的間距都為da；虛擬陣元沿x軸和y軸的數(shù)目都為Ma。最終，可以用MT＋MR個陣元獲得MaTMaR×MaTMaR的虛擬的二維均勻陣列，該虛擬陣列邊長為La＝（n－1）da，且n＝MaTMaR。

在三維成像時，可以采用多發(fā)多收模式：被選中的發(fā)射單元同時發(fā)射不同中心頻率的寬帶信號，被選中的接收單元接收所有回波，然后在信號處理端通過數(shù)字濾波來獲得相應(yīng)的基帶信號。然而，此時需要采集端有較大的帶寬才能同時接收所有發(fā)射陣元的回波，否則，只能采用單發(fā)多收模式。相對于多發(fā)多收模式，單發(fā)多收模式可以靈活地控制和激活所需要的收發(fā)陣元，并能降低系統(tǒng)功耗。根據(jù)PCA原理，在單發(fā)多收模式下，一個發(fā)射陣元可以獲得與之對應(yīng)的多個虛擬陣元，并且每個虛擬陣元同時具有收／發(fā)功能。因此，與傳統(tǒng)二維陣列相比，單發(fā)多收模式下隨機二維陣列所需要的PRF大大降低。

1.2 隨機陣元數(shù)目的選擇

已經(jīng)分析了如何實現(xiàn)虛擬隨機二維陣列，為了簡化分析，后文的討論都是建立在第1.1節(jié)等效后的虛擬隨機二維陣列的基礎(chǔ)之上。然而，如何確定合適的隨機陣元數(shù)目，才能夠保證有較高質(zhì)量的三維微波成像效果呢？

文獻(xiàn)［2］定義了隨機陣列沿指定方向的模糊函數(shù)，該模糊函數(shù)主要取決于隨機陣列沿指定方向上的投影。為了保證陣列具有三維分辨能力和能夠充分利用陣列的長度，必須要求這些隨機陣元能夠在一對相互正交的方向上投影出滿陣列，這樣就可以在這一對正交方向上獲得與等效均勻線陣相當(dāng)?shù)木劢鼓芰?。例如，有若干個隨機陣元在陣列平面中心沿著x軸和y軸分布，這樣就可以構(gòu)成一個“十”字形正交陣列。顯然，要獲得更好的三維成像能力，則隨機陣元的數(shù)目Na不能少于2n－1。從N＝（MaTMaR）2個陣元網(wǎng)格中隨機選取Na個位置來布置陣元，那么能夠滿足投影要求的排列的概率可以表示為

式中，F(xiàn)Na∈R＋表示重復(fù)度，并且是關(guān)于Na的非線性函數(shù)。從物理意義上分析，當(dāng)Na增大時，隨機陣元能夠覆蓋的“十”字形陣列的可能性就越大，即p（Na）是關(guān)于參數(shù)Na的單調(diào)不減函數(shù)。FNa在實際中很難確定，如果令其取1，則可以獲得Na＜N時函數(shù)p（Na）的上界；而當(dāng)Na＝N時，p（Na）≡1。由此，可以得到

由式（2）可以估計出使p（Na）滿足給定概率門限的Na的值。一般取Na為n的整數(shù)倍，如果成像效果不理想，可以逐步提高倍數(shù)。實驗發(fā)現(xiàn)當(dāng)Na約占20%的網(wǎng)格時，就能取得較好的成像效果。

當(dāng)N較小時，可以通過仿真實驗來確定合適的Na值。如圖3所示，通過2萬次蒙特卡羅仿真試驗，可以繪制出大概的p（Na）曲線，其步驟如下：

步驟1 生成隨機陣元位置序列，然后向x和y方向進(jìn)行投影；

步驟2 統(tǒng)計存在空行、空列以及同一陣元在x和y方向都被投影的情況。然后就可以繪制出p（Na）的曲線。

圖3 不同陣元間隔時p（Na）曲線

2 隨機二維陣列的三維微波成像

2.1 信號模型

假定隨機陣列工作在單發(fā)多收模式下，等效后的隨機陣列如圖1所示。盡管在一個PRF內(nèi)可以獲得多個陣元的數(shù)據(jù)，但是為了描述方便，假設(shè)在每個PRF內(nèi)，只有第i（i＝0，1，…，Na－1）個陣元Γi收／發(fā)信號，用符號t表示快時間，那么，對于觀測場景Ω中任一散射目標(biāo)P∈Ω其回波為

式中，I（P）表示散射目標(biāo)P的散射系數(shù)；wa（Γi）表示陣列對Γi的加權(quán)，可假設(shè)為1；pr［·］表示發(fā)射的線性調(diào)頻信號的波形；R（P，Γi）＝2‖P－Γi‖2表示雙程距離歷史；c表示電磁波傳播速度；λ為發(fā)射信號波長。那么，對于整個觀測場景Ω，陣元Γi∈Ξ（Ξ＝，觀測到的回波可以表示為

2.2 BP算法簡介

由于隨機二維陣列的陣元是隨機選擇的，所以其不再適合利用傳統(tǒng)的RD算法［6-7］、線頻調(diào)變標(biāo)（chirp scaling， CS）算法［1112］等進(jìn)行微波成像。由于BP算法［1315］適合于任意的天線相位中心軌跡，因此，BP算法成為最合適的選擇。首先對回波信號進(jìn)行距離壓縮，那么壓縮后的時域信號可以表示為BP算法的實質(zhì)是時域相關(guān)積累。于是，只需要對陣元Γi的回波進(jìn)行目標(biāo)Q∈Ω的相位補償，便可以得到Q的單陣元圖像（Q；Γi）為

式中，tQ（Γi）＝R（Q，Γi）／c，表示目標(biāo)Q和陣元Γi之間的回波延遲。最終，將所有單陣元的圖像進(jìn)行累加，就可以得到目標(biāo)Q的BP成像結(jié)果為

在這里需要強調(diào)的是，ssc（tQ（Γi）；Γi）需要利用插值計算來完成。在傳統(tǒng)的BP成像算法中，插值計算使用辛克函數(shù)插值［15］。在后文中，將使用NUFFT［8，16-18］來實現(xiàn)該插值過程。

3 NUFFT BP及加速實現(xiàn)

由于R（Q，Γi）相對于目標(biāo)Q而言是非均勻的，所以無法直接利用傳統(tǒng)的快速傅里葉變換（fast Fourier transform，F(xiàn)FT）來快速求解ssc（tQ（Γi）；Γi）。然而，Type-ⅡNUFFT恰好能夠?qū)崿F(xiàn)從均勻頻域數(shù)據(jù)插值得到高精度的非均勻時域數(shù)據(jù)。因此，本節(jié)利用快速高斯柵格的Type-ⅡNUFFT［8］來完成ssc（tQ（Γi）；Γi）的高精度插值。

3.1 NUFFT BP算法的提出

對于選定的陣元Γi，利用一維Type-ⅡNUFFT可以同時實現(xiàn)對成像空間中所有的時間節(jié)點tQ（Γi）進(jìn)行插值。根據(jù)離散時間信號理論，ssc（tQ（Γi）；Γi）可用離散傅里葉變換進(jìn)行求解，即

式中，sSc（fk；Γi）＝F［ssc（t；Γi）］，表示ssc（t；Γi）的傅里葉變換；fk為快時間頻率；t′Q（Γi）＝2πFstQ（Γi）／Nr，表示將tQ（Γi）歸一化到［0，2π］，Nr為距離向采樣點數(shù)；Fs為回波信號的采樣頻率。

定義輔助函數(shù)ssc，－τ（t；Γi）為

式中，F(xiàn)－1［·］表示逆傅里葉變換；u＝0，1，…，Mτ－1，Mτ＝2Nr為輔助函數(shù)的點數(shù)。然后，利用高斯核函數(shù)（x∈［0，2π］）對輔助函數(shù)ssc，－τ（t；Γi）進(jìn)行插值，從而就可以得到ssc（tQ（Γi）；Γi），即

式中，符號*表示卷積。

通常τ非常小，所以式（10）中的高斯核函數(shù)為類脈沖函數(shù)，實際中其截斷長度選擇為2Msp即可，其中Msp表示截斷的高斯核函數(shù)的一半的點數(shù)。當(dāng)Msp＝6時，NUFFT可以獲得6位數(shù)值精度，而Msp＝12時，則可以獲得12位數(shù)值精度［8］。具體在操作時，參數(shù)τ的取值通常選擇為Msp／N2r。在本節(jié)中，由于BP算法的具體實現(xiàn)過程最終是由NUFFT插值來完成的，所以稱這種BP算法為NUFFT BP算法。

3.2 CUDA加速實現(xiàn)

雖然，上述的NUFFT算法快速實現(xiàn)了ssc（tQ（Γi）；Γi）的插值。但是，當(dāng)成像空間的點數(shù)增多時，求解R（Q，Γi）的計算量仍然很大。例如，對60×60×60點（21.6萬點）的區(qū)域進(jìn)行成像時，對于陣元Γi需要計算21.6萬次R（Q，Γi）。利用圖形處理單元（graphics processing unit，GPU）并行計算技術(shù)，可以將求解R（Q，Γi）并行化，同時NUFFT插值也可以進(jìn)行并行處理，從而能夠使得BP成像處理效率得到進(jìn)一步提高。CUDA加速的NUFFT BP成像過程可以用如圖4所示的偽代碼進(jìn)行描述。

圖4 CUDA加速的NUFFT BP

圖4 中，核函數(shù)cuda_Range HistoryQGamma用來求解Γi和場景中所有像素點的距離歷史；核函數(shù)cuda_Nu FFTInterp完成距離歷史歸一化插值時間節(jié)點和NUFFT插值過程；核函數(shù)cuda_Phase Compensation And Accumulation完成相位補償和積累操作。

需要強調(diào)的是，用CUDA來加速NUFFT時，能夠?qū)崿F(xiàn)的點數(shù)受CUDA能夠?qū)崿F(xiàn)的FFT的最多的點數(shù)（約226復(fù)數(shù)點）和GPU全局內(nèi)存大小的限制［19］。由于成像過程中不止NUFFT占用GPU內(nèi)存，所以通常假設(shè)成像時NUFFT處理的點數(shù)不超過800萬點。當(dāng)成像區(qū)域像素點數(shù)大于這一限制時，需要對成像區(qū)域進(jìn)行適當(dāng)?shù)貏澐衷偬幚?。把用CUDA加速實現(xiàn)的NUFFT BP成像過程簡稱為CUDA NUFFT BP。

3.3 運算量分析

本小節(jié)對傳統(tǒng)BP算法（使用辛克函數(shù)插值），NUFFT BP算法和CUDA NUFFT BP算法的運算量進(jìn)行簡要分析。由于乘法運算和開方運算對執(zhí)行效率影響較大，故分析時主要考慮乘法運算和開方運算的運算量。

如圖4所示，NQ表示整個成像空間的像素點數(shù)，sX、sY和sZ分別表示成像空間的沿x軸、y軸和z軸的像素點數(shù)，于是NQ＝sX×sY×sZ。辛克函數(shù)插值核點數(shù)用MBP表示，NQ點數(shù)據(jù)開方操作的運算量為Osq（NQ）。

傳統(tǒng)BP算法的運算主要由距離壓縮（包括FFT、頻域相乘和逆快速傅里葉變換（inverse fast Fourier transform，IFFT））、距離歷史計算（包括乘法和開方運算）和辛克函數(shù)插值（查表實現(xiàn)）構(gòu)成。于是，傳統(tǒng)BP算法的運算量大約為

如果采用直接計算辛克函數(shù)的值，那么插值操作的運算量大約為10MBPNQNa。

NUFFT BP算法的運算主要包括計算頻域距離壓縮數(shù)據(jù)，距離歷史計算以及NUFFT插值。其中，直接計算頻域距離壓縮數(shù)據(jù)，可以省掉IFFT的運算量；同時，逐陣元方式還可以有效降低距離歷史計算的運算量。于是，NUFFT BP算法的運算量大約為

CUDA NUFFT BP算法為了便于并行化處理，其距離歷史計算和傳統(tǒng)BP算法一樣，其余運算量同NUFFT BP算法。那么，CUDA NUFFT BP算法的運算量大約為

觀察式（11）～式（13）可以發(fā)現(xiàn)，當(dāng)成像的像素點數(shù)NQ較大時，3種算法的運算量主要取決于插值運算和距離歷史計算。其中，NUFFT BP算法的運算量和傳統(tǒng)BP算法用查表實現(xiàn)時基本相當(dāng)，且低于直接計算辛克函數(shù)的值的方式；同時，NUFFT插值精度也遠(yuǎn)優(yōu)于辛克函數(shù)插值。

4 仿真實驗

本節(jié)利用計算機仿真實驗來驗證隨機二維陣列實現(xiàn)三維微波成像的可行性。同時，本節(jié)也進(jìn)行了CUDA加速的NUFFT BP的性能實驗。實驗中，設(shè)置9個點目標(biāo)，它們的位置分別為（0 m，0 m，0 m）、（30 m，－30 m，10 m）、（－30 m，30 m，10 m）、（－30 m，－30 m，10 m）、（30 m，30 m，10 m）、（30 m，－30 m，－10 m）、（－30 m，30 m，－10 m）、（－30 m，－30 m，－10 m）和（30 m，30 m，－10 m）。其他主要參數(shù)如表1所示。

表1 陣列參數(shù)和波長

4.1 陣列成像對比實驗

圖5 不同陣元分布及其三維成像結(jié)果

圖5 展示了4種不同的陣列布局方案，以及它們對應(yīng)的9個點目標(biāo)的三維成像結(jié)果（10 dB動態(tài)范圍）。從圖5（a）和圖5（e）可以看出，“十”字形陣列布局盡管可以實現(xiàn)目標(biāo)三維聚焦，但是對于非中心點目標(biāo)，旁瓣非常明顯。當(dāng)陣元數(shù)目增加構(gòu)成“米”字形陣列時，旁瓣效應(yīng)明顯減弱，但是在鄰近目標(biāo)的幾個分辨單元內(nèi)旁瓣依然較強，如圖5（b）和圖5（f）所示。利用和“米”字形陣列數(shù)目（4n－3個）相近的隨機陣元（4n個）成像，如圖5（c）和圖5（g）所示，其成像結(jié)果雖然也會散布一些非目標(biāo)強散射單元，但相比“米”字形陣列成像結(jié)果，它們在成像空間里更分散，強度更弱。而當(dāng)使用滿陣元進(jìn)行成像時，就得到非常理想的成像結(jié)果，如圖5（d）和圖5（h）所示。

進(jìn)一步對點目標(biāo)成像性能進(jìn)行分析。分別取占全部陣元20%、40%和100%的陣元數(shù)目對目標(biāo)點（0，0，0）進(jìn)行成像，并沿著y軸作該點目標(biāo)的一維剖面圖，結(jié)果如圖6所示。由圖6可以發(fā)現(xiàn)，在上述3種情況下，剖面圖的主瓣寬度幾乎一致，因此其分辨率基本相同（約29．2 m）。這是因為足夠數(shù)目的隨機陣元在陣列平面內(nèi)均勻隨機分布，能夠在x軸和y軸方向都投影出滿陣元線性陣列，所以使得3種情況的有效陣列長度均相同，從而它們的分辨率基本一致。由圖6還可發(fā)現(xiàn)，上述3種情況的峰值旁瓣比也基本相同（約－13．3 dB）。

圖6 不同數(shù)目隨機陣元時點目標(biāo)剖面圖

需要注意，圖6并不能反映出積分旁瓣比的性能，因為陣列天線產(chǎn)生的是角分辨率，這造成點目標(biāo)的旁瓣會向陣列方向發(fā)生彎曲，如圖5（e）所示。然而綜合圖5和圖6可知，積分旁瓣比隨著所選用的陣元數(shù)目的增加可以得到不斷改善，最終達(dá)到和滿陣元時一致。

由圖5和圖6，容易發(fā)現(xiàn)，當(dāng)陣元均勻隨機分布時，只要選擇適當(dāng)數(shù)目的隨機陣元，就能使目標(biāo)的三維聚焦效果和滿陣元時相當(dāng)；當(dāng)陣元數(shù)目增加時，其聚焦效果也逐步逼近滿陣元的結(jié)果。在實驗過程中還可以得出如下結(jié)論：一般地，當(dāng)隨機陣元數(shù)目占總網(wǎng)格數(shù)的18%～20%時，已經(jīng)可以有很好的三維聚焦效果（10 d B動態(tài)范圍）。這意味著即使在使用滿陣元成像時，也可以先利用20%的陣元聚焦，然后提取10 d B動態(tài)范圍，再積累剩余陣元的數(shù)據(jù)，這樣操作大約可以提高3～5倍的成像效率。

4.2 NUFFT插值精度實驗

本小節(jié)比較NUFFT插值和辛克函數(shù)插值的插值精度。其中，NUFFT插值核長度為12；辛克函數(shù)插值核長度為16，并采用直接計算辛克函數(shù)的值的方式。選擇某一陣元數(shù)據(jù)，對1 331個插值節(jié)點使用這兩種方法分別插值。精確的插值結(jié)果則按式（8）用距離壓縮頻域數(shù)據(jù)進(jìn)行離散傅里葉變換得到。圖7展示了利用離散傅里葉變換精確插值的結(jié)果，以及辛克函數(shù)插值和NUFFT插值的誤差。

圖7 辛克函數(shù)插值和NUFFT插值的誤差

由圖7可以發(fā)現(xiàn)，辛克函數(shù)插值的誤差分布在0～1；而NUFFT插值的誤差則分布在0～4×10－6，這和前文提到的6位數(shù)值精度是一致的。上述實驗結(jié)果表明12點NUFFT插值的精度遠(yuǎn)遠(yuǎn)優(yōu)于16點辛克函數(shù)插值。根據(jù)前文運算量分析可知，12點NUFFT插值所需要的運算量也遠(yuǎn)低于16點辛克函數(shù)插值。

4.3 CUDA性能實驗

為了充分驗證算法的性能和有效性，本節(jié)選擇兩臺不同配置計算機進(jìn)行實驗。計算機I運行環(huán)境如下：操作系統(tǒng)為Windows 7 64bit，Matlab 2012b，Intel i5－3550 CPU（3．30GHz），內(nèi)存8G，Nvidia GeForce GTX 670以及CUDA 5．5。計算機Ⅱ的運行環(huán)境如下：Windows 7 64bit，Matlab 2012b，Intel i5－2320 CPU（3GHz），內(nèi)存12G，Nvidia GeForce GT 560Ti以及CUDA 5．5。

根據(jù)前文的分析可以知道，CUDA NUFFT BP按逐陣元方式批量處理成像區(qū)域的所有像素點，所以，要處理的像素點越多，就更能體現(xiàn)出CUDA的優(yōu)勢。設(shè)計兩組實驗，利用滿陣元方式成像，分別選擇61×61×61點（實驗Ⅰ）和121×121×121點（實驗Ⅱ）的成像區(qū)域。利用Matlab的GPU并行處理功能實現(xiàn)距離歷史計算的NUFFT BP稱為Matlab GPU NUFFT BP?？疾?種BP成像方法：傳統(tǒng)BP、NUFFT BP、Matlab GPU NUFFT BP和CUDA NUFFT BP。需要說明的是，前3種方法是利用Matlab實現(xiàn)的；傳統(tǒng)BP的插值為辛克函數(shù)插值；而其余方法都是利用NUFFT來實現(xiàn)BP的插值過程。另外，由于Matlab對于循環(huán)計算執(zhí)行效率很低，這直接影響了距離歷史計算和辛克函數(shù)插值的效率，不能真實反映算法本身的運算量。因此，本文將距離歷史計算和辛克函數(shù)插值都采用C語言實現(xiàn)，然后由Matlab進(jìn)行調(diào)用；并且辛克函數(shù)插值也采用批量計算插值方式，而不是循環(huán)逐點插值方式。

由表2可以看出，當(dāng)成像搜索區(qū)域有226 981點（實驗Ⅰ）時，與傳統(tǒng)BP相比，利用NUFFT插值的BP都在一定程度上節(jié)約了計算時間，而CUDA NUFFT BP的加速效果則更為顯著。對于CUDA NUFFT BP，計算機Ⅰ加速倍數(shù)約為31倍，計算機Ⅱ加速倍數(shù)約為22倍。兩臺計算機的CUDA NUFFT BP性能不同，主要是和它們的GPU并行計算能力有關(guān)。

表2 CUDA性能實驗Ⅰ結(jié)果

當(dāng)像素點數(shù)增加到1 771 561點（實驗Ⅱ）時，利用CUDA并行計算距離歷史的優(yōu)勢就更為明顯，如表3所示。由于Matlab利用GPU并行求解距離歷史涉及到頻繁的內(nèi)存拷貝，所以Matlab GPU NUFFT BP的效率提高地并不顯著。對于計算機Ⅰ，其CUDA NUFFT BP相對于Matlab GPU NUFFT BP的加速倍數(shù)達(dá)到了約20倍；而計算機Ⅱ也有約13倍的加速倍數(shù)。

表3 CUDA性能實驗Ⅱ結(jié)果

綜合表2和表3，可以得出如下結(jié)論：NUFFT插值能夠提高BP成像算法插值過程的效率；由于CUDA NUFFT BP并行實現(xiàn)了計算距離歷史以及并行實現(xiàn)了NUFFT，從而顯著地加速了BP成像算法的執(zhí)行效率。

5 結(jié) 論

本文首先利用隨機二維陣列進(jìn)行三維微波成像，降低了對PRF的要求和系統(tǒng)功耗。然后，基于等效相位中心原理，本文提出了一種利用二維MIMO陣列天線等效虛擬隨機二維陣列的方法，結(jié)果降低了陣列成本，同時因為該方法容易實現(xiàn)任意隨機分布的二維陣列，所以又增強了陣列設(shè)計的靈活性。研究結(jié)果表明，通過TypeⅡNUFFT實現(xiàn)BP的插值過程，可以使得BP的計算效率得到一定程度地提高，而利用CUDA并行運算又能夠進(jìn)一步地提高BP的計算效率。最后，仿真實驗證實了利用隨機二維陣列進(jìn)行三維微波成像的有效性和可行性，其中，NUFFT BP的插值精度要遠(yuǎn)優(yōu)于傳統(tǒng)BP的辛克函數(shù)插值的插值精度，并且在算法性能方面，CUDA NUFFT BP相對于傳統(tǒng)的BP算法也有極大地提升。另外，CUDA NUFFT BP不僅適用于隨機二維陣列的三維成像，而且還可以推廣到其他任意合成孔徑雷達(dá)的成像模式之中。

［1］Shi J，Zhang X L，Yang J Y，et al．APC trajectory design for“one-active”linear-array three-dimensional imaging SAR［J］. IEEE Trans.on Geoscience and Remote Sensing，2010，48（3）：1470- 1486.

［2］Shi J，Zhang X L，Xiang G，et al．Signal processing for microwave array imaging：TDC and sparse recovery［J］.IEEE Trans. on Geoscience and Remote Sensing，2012，50（11）：4584- 4598.

［3］Ender J H，Klare J．System architectures and algorithms for radar imaging by MIMO-SAR［C］∥Proc.of the IEEE Radar Conference，2009：1- 6.

［4］Xie W C，Zhang X L，Shi J．MIMO antenna array design for airborne down-looking 3D imaging SAR［C］∥Proc.of the IEEE International Conference on Signal Processing Systems，2010：452- 456.

［5］Zhang X L，Xiang G，Shi J．An implementation of sparse MIMO array for radar antenna［P］．China：201210000517，2013- 06-12．（張曉玲，向高，師君．一種稀疏MIMO平面陣列雷達(dá)天線構(gòu)建方法［P］．中國：201210000517，2013- 06- 12．）

［6］Wu Q S，Zhang Y D，Amin M G，et al．Focusing of tandem bistatic SAR data using range Doppler algorithm［C］∥Proc.of the IEEE Radar Conference，2014：927- 931.

［7］Zare A，Masnadi-Shirazi M A，Samadi S．Range-Doppler algorithm for processing bistatic SAR data based on the LBF in the constant-offset constellation［C］∥Proc.of the IEEE Radar Conference，2012：17- 21.

［8］Greengard L，Lee J Y．Accelerating the nonuniform fast Fourier transform［J］.Society for Industrial and Applied Mathematics Review，2004，46（3）：443- 454.

［9］Jin X X，Seok-Bum K．GPU-based parallel implementation of SAR imaging［C］∥Proc.of the International Symposium on Electronic System Design，2012：125- 129.

［10］Malanowski M，Krawczyk G，Samczynski P，et al．Real-time high-resolution SAR processor using CUDA technology［C］∥Proc.of the International Radar Symposium，2013：673- 678.

［11］Li D，Liao G S，Wang W，et al．Extended azimuth nonlinear chirp scaling algorithm for bistatic SAR processing in high-resolution highly squinted mode［J］.IEEE Geoscience and Remote Sensing Letters，2014，11（6）：1134- 1138.

［12］Chen S C，Xing M D，Yang T L，et al．A nonlinear chirp scaling algorithm for tandem bistatic SAR［C］∥Proc.of the IEEE International Geoscience and Remote Sensing Symposium，2013：2485- 2488.

［13］Wang W，Ma Y H，Wang X P．High computation efficiency BP imaging algorithm for MIMO radar［J］.Systems Engineering and Electronics，2013，35（10）：2080- 2085．（王偉，馬躍華，王咸鵬．一種高運算效率的MIMO雷達(dá)BP成像算法［J］．系統(tǒng)工程與電子技術(shù)，2013，35（10）：2080- 2085．）

［14］Yu M，Zhang X L，Liu Z．Acceleration of fast factorized back projection algorithm for bistatic SAR［C］∥Proc.of the IEEE International Geoscience and Remote Sensing Symposium，2013：2493- 2496.

［15］Shi J．Researches on the principle and imaging technology of bistatic SAR and linear array SAR［D］．Chengdu：University of Electronic Science and Technology，2009．（師君．雙基地SAR與線陣SAR原理及成像技術(shù)研究［D］．成都：電子科技大學(xué)，2009．）

［16］Yang H G，Yuan Y，Wu J J，et al．NUFFT applied to motion compensation in the near-space SAR imaging［C］∥Proc.of the IEEE International Geoscience and Remote Sensing Symposium，2012：3883- 3886.

［17］Sun S L，Zhu G F．Accelerating CS radar imaging by NUFFT［C］∥Proc.of the IEEE International Conference on Signal Processing，Communication and Computing，2013：1- 4.

［18］Capozzoli A，Curcio C，Liseno A，et al．NUFFT-based SAR back projection on multiple GPUs［C］∥Proc.of the IEEE Tyrrhenian Workshop on Advances in Radar and Remote Sensing，2012：62- 68.

［19］NVIDIA C．CUFFT library［EB／OL］．［2014- 2- 13］．http：∥docs．nvidia．com／cuda／cufft／

Fast 3D imaging method for random array based on NUFFT BP

XIANG Gao，ZHANG Xiao-ling，WU Zong-liang
（School of Electronic Engineering，University of Electronic Science and Technology，Chengdu 611731，China）

Traditional two dimensional（2D）antenna array requires high pulse repetition frequency（PRF）to switch transmitters，when it is applied to 3D microwave imaging．Therefore，3D imaging is gotten by the random antenna array，as a result it reduces the requirement for PRF．Besides，a reasonable plan based on the phase center approximation（PCA）is proposed to implement the random antenna array of arbitrary geometry．The back projection（BP）algorithm interpolated by the nonuniform fast Fourier transform（NUFFT）is called NUFFT BP，and achieves high precision interpolation with relatively less computational complexity．Moreover，the NUFFT BP is further accelerated by the compute unified device architecture（CUDA）parallel computing technology，named CUDA NUFFT BP．It significantly improves the execution efficiency of the BPalgorithm．

random antenna array；nonuniform fast Fourier transform（NUFFT）；parallel computing；back projection（BP）

TN95

A

10．3969／j．issn．1001-506X．2015．04．01

向 高（1985-），男，博士研究生，主要研究方向為雷達(dá)信號處理技術(shù)、自適應(yīng)信號處理、稀疏信號處理技術(shù)。E-mail：moran422＠gmail．com

1001-506X（2015）04-0725-07

2014- 05- 09；

2014- 09- 28；網(wǎng)絡(luò)優(yōu)先出版日期：2014- 11- 21。

網(wǎng)絡(luò)優(yōu)先出版地址：http：／／w ww．cnki．net／kcms／detail／11．2422．TN．20141121．0940．010．html

國家自然科學(xué)基金（61101170）；高等學(xué)校博士學(xué)科點專項科研基金（20110185110001）資助課題

張曉玲（1964-），女，教授，博士研究生導(dǎo)師，博士，主要研究方向為雷達(dá)信號處理技術(shù)、分類與識別、導(dǎo)航定位技術(shù)。E-mail：xlzhang＠uestc．edu．cn

吳宗亮（1979-），男，博士研究生，主要研究方向為雷達(dá)信號處理技術(shù)、自適應(yīng)信號處理、雷達(dá)目標(biāo)識別。E-mail：nanyuanfeixue＠sina．com