摘 要：后向投影（back projection， BP）算法是一種精確的時域成像算法，但BP算法的計(jì)算復(fù)雜度高，難以實(shí)現(xiàn)實(shí)時性成像，特別是在考慮三維成像時，BP算法的計(jì)算復(fù)雜度會進(jìn)一步增加。提出一種應(yīng)用在合成孔徑聲納（synthetic aperture sonar， SAS）上的三維快速分解BP（fast factorized BP， FFBP）成像算法，并利用圖形處理器（graphics processing unit， GPU）加速三維FFBP算法。經(jīng)過對點(diǎn)目標(biāo)的測試，計(jì)算時間從原本的263 s降低到了2.3 s，解決了SAS中的三維成像實(shí)時性問題。同時，驗(yàn)證了所提算法在非理想航跡下的成像效果。結(jié)果表明，在添加幅度不超過0.1 m（一個波長以內(nèi)）的正弦擾動時，所提算法對點(diǎn)目標(biāo)仍有良好的聚焦效果。

關(guān)鍵詞： 快速分解后向投影; 并行計(jì)算; 圖形處理器; 合成孔徑聲納; 三維成像

中圖分類號： TN 911.7 文獻(xiàn)標(biāo)志碼： A""" DOI：10.12305/j.issn.1001-506X.2024.10.01

Three-dimensional fast factorized back projection SAS imaging

algorithm accelerated by GPU

TAO Hongbo^1， ZHANG Dongsheng HUANG Yong³

（1. Laboratory of Ocean Acoustic Technology， Institute of Acoustics， Chinese Academy of Sciences，

Beijing 100190， China; 2. School of Electronic， Electrical and Communication Engineering，

University of Chinese Academy of Sciences， Beijing 100049， China; 3. Ningbo “Five-in-One” Campus

Education Development Center， Zhejiang University， Ningbo 315100， China）

Abstract： The back projection （BP） algorithm is an accurate time-domain imaging algorithm， but the high computational complexity of the BP algorithm makes it difficult to achieve real-time imaging. Especially when three-dimensional imaging is considered， the computational complexity of the BP algorithm will further increase. The paper proposes a three-dimensional fast factorized BP （FFBP） imaging algorithm applied to synthetic aperture sonar （SAS）， and accelerates the three-dimensional FFBP algorithm using a graphics processing unit （GPU）. After testing on point targets， the calculation time is reduced from the original 263 s to 2.3 s， solving the real-time imaging problem in SAS. At the same time， the imaging effect of the proposed algorithm under non-ideal trajectories is verified. The results show that when a sinusoidal disturbance with an amplitude not exceeding 0.1 m （within one wavelength） is added， there is still a good focusing effect on the point target.

Keywords： fast factorized back projection （FFBP）; parallel computing; graphics processing unit （GPU）; synthetic aperture sonar （SAS）; three-dimensional imaging

0 引 言

合成孔徑聲納（synthetic aperture sonar， SAS）通過收集一連串的腫沖相干組合來提高聲納圖像的沿航跡分辨率，經(jīng)過信號相干疊加后，SAS的航跡向分辨率不再與成像距離和工作頻率有關(guān)^［14］。如今，為了提高區(qū)域覆蓋率，大多采用多子陣SAS，即一個發(fā)射機(jī)與多個接收水聽器組成的陣列^［56］。后向投影（back projection， BP）算法是一種精確的時域成像算法^［710］，相比于頻域成像算法，BP算法有更好的保相性，適用于任意的聲納陣列構(gòu)型，可以更好地適應(yīng)多子陣SAS中不均勻采樣的問題。1983年，Newbold^［11］成功將計(jì)算機(jī)層析（computer tomography， CT）技術(shù)引入到聚束合成孔徑雷達(dá)（synthetic aperture radar， SAR）成像領(lǐng)域，這為后來的時域成像算法發(fā)展提供了方向。但由于時域成像算法的計(jì)算負(fù)載太大，當(dāng)時并沒有被應(yīng)用到實(shí)際設(shè)備中。在1999年的雷達(dá)會議上，Yegulalp^［12］正式提出了快速BP（fast BP， FBP）算法，F(xiàn)BP算法將BP算法的計(jì)算復(fù)雜度從O（N3）降低到了O（N^2.5），并指出在極坐標(biāo)系下可以實(shí)現(xiàn)較粗的分辨率成像而不損失圖像信息。2000年前后，瑞典國防研究局（Swedish Defence Research Agency）的Ulander等^［13］正式提出了快速分解BP（fast factorized BP， FFBP）算法，F(xiàn)FBP算法的計(jì)算復(fù)雜度為O（N2log₂N），F(xiàn)FBP算法極大地提高了BP算法的計(jì)算效率，使得時域成像算法應(yīng)用于實(shí)時性成像設(shè)備成為了可能。近年來，Moon等^［14］研究了一種基于雷達(dá)的四叉樹BP算法，這是在直角坐標(biāo)系下的二維FFBP算法，也有學(xué)者研究了BP算法與FFBP算法在利用圖形處理單元（graphics processing unit， GPU）與現(xiàn)場可編程門陣列（field programmable gate array， FPGA）實(shí)現(xiàn)并行計(jì)算時的加速效果^［1520］。Cao等^［21］研究了基于FPGA下的SAR實(shí)時BP算法的并行優(yōu)化與實(shí)現(xiàn)，Wielage等^［22］研究了FFBP算法在GPU與FPGA下高性能與低功耗的表現(xiàn)情況。

由于二維成像只能獲得目標(biāo)的航跡距離向信息，無法直接獲得目標(biāo)的深度信息，因此想要直接獲得目標(biāo)的三維特征就必須考慮三維成像。在SAS上，二維成像是通過布置一條線陣實(shí)現(xiàn)的，拓展到三維成像時，需要沿垂直航跡向航跡向布置一平面接收陣，因此在考慮三維成像時，算法的計(jì)算復(fù)雜度會急劇上升。本文提出了適用于SAS的三維FFBP算法，通過沿航跡向與垂直航跡向的子孔徑劃分降低算法的計(jì)算復(fù)雜度，并通過中央處理器（central processing unit， CPU）與GPU加速效果的對比，將本文所提出的三維FFBP算法利用GPU實(shí)現(xiàn)大規(guī)模并行計(jì)算，以滿足SAS三維成像實(shí)時性的需求。同時，對比了FFBP算法與BP算法同時采用GPU加速后的計(jì)算時間以及成像質(zhì)量。最后，在仿真軟件中驗(yàn)證了本文所提算法在非理想航跡下的成像效果。

1 FFBP算法的原理

FFBP算法主要由Ulander等^［13］提出，相比于BP算法，F(xiàn)FBP算法快速化的實(shí)質(zhì)就是由于初始階段子孔徑的角域帶寬很小，相應(yīng)的角域采樣率可以很低，應(yīng)用較少的角域采樣點(diǎn)數(shù)就能得到 “粗”角分辨率圖像，即沒有必要將每個孔徑對應(yīng)的數(shù)據(jù)BP到成像網(wǎng)格的每個像素點(diǎn)^［23］。如圖1所示，其中t為孔徑點(diǎn)位置，r為原點(diǎn)到真實(shí)像素點(diǎn)與投影像素點(diǎn)的位置，R（·）為孔徑點(diǎn)到真實(shí)像素點(diǎn)與投影像素點(diǎn)的位置。這種近似不可避免地會帶來距離誤差，Ulander等^［13］在論文中給出了距離誤差的限制范圍：

ΔR≤dDr_min（1）

式中：d為子孔徑沿航跡向長度;D為子圖像沿航跡向長度;r_min為接收陣到目標(biāo)的最短距離。

同時，F(xiàn)FBP算法對角度與距離的采樣要求為

Δcosv≤λ_min2d（2）

Δr≤speed2B（3）

式中：λ_min為最短波長;speed為聲速;B為信號帶寬。

FFBP算法大體可分為3個階段：初始階段、迭代階段、最終階段^［24］。如圖2所示，初始階段可以認(rèn)為是先將整個孔徑劃分成一定數(shù)量的初始子孔徑，并根據(jù)式（1）設(shè)置好子圖像尺寸，計(jì)算每個子孔徑上的脈沖點(diǎn)與每個子圖像上的像素點(diǎn)距離并相干疊加得到初始圖像。迭代階段是FFBP算法的核心，隨著子孔徑的合并，子圖像的尺寸也隨之變短，所需像素點(diǎn)也隨之增多，這時便從粗分辨率的初始圖像不斷融合成為高分辨的中間圖像。在最后階段，當(dāng)子孔徑融合成為整個孔徑時，這時生成的圖像便是最終圖像。

FFBP算法的計(jì)算效率與選取初始子孔徑長度有關(guān)，可以認(rèn)為初始子孔徑長度越短，其計(jì)算效率越高，因?yàn)榇藭r所需像素點(diǎn)較少，但是當(dāng)初始子孔徑長度選取得過于短時，可能會影響成像質(zhì)量，甚至導(dǎo)致無法成像，因此對于初始子孔徑長度的選取需要根據(jù)實(shí)際情況綜合考慮，F(xiàn)FBP算法是一種計(jì)算效率與圖像質(zhì)量折衷的算法。

2 GPU加速下的三維FFBP算法

2.1 三維FFBP算法原理

Moon等^［14］提出了一種基于雷達(dá)的四叉樹二維BP算法，這是一種二維FFBP算法，其采用直角坐標(biāo)系進(jìn)行成像，相比于極坐標(biāo)系，直角坐標(biāo)系更加直觀且簡潔。本文將這種四叉樹二維BP算法拓展到三維并應(yīng)用在了SAS上，下面將對算法進(jìn)行簡要介紹。

對于二維成像，SAS采集到的回波信號表達(dá)式為

s（τ，t_m）=sinspeed·［Bτ－2R（t_m）speed］·

exp［－j4πλR（t_m）］（4）

式中：j為虛數(shù)；sinspeed（τ）=［sin（πτ）］/（πτ）;τ為距離向時間;t_m為航跡向時間;λ為信號波長。

BP算法的表達(dá)式為

I（x，y）=∑tms（τ，t_m）exp［j4πλR_p（t_m）］（5）

式中：R_p（t_m）為每個脈沖點(diǎn)到每個像素點(diǎn)的距離。

可以發(fā)現(xiàn)，BP算法需要計(jì)算每一個脈沖點(diǎn)與像素點(diǎn)的距離，在考慮三維成像時其計(jì)算復(fù)雜度會急劇增加，因此本文提出了一種適用于SAS的三維FFBP算法，以降低算法的運(yùn)算量。對于二維FFBP算法，子孔徑劃分是通過一條虛擬線陣進(jìn)行的，拓展到三維FFBP算法，其子孔徑劃分則是通過虛擬面陣來實(shí)現(xiàn)的，需要同時考慮沿航跡向與垂直航跡向的子孔徑劃分，才能最大程度地提升算法的計(jì)算效率。

如圖3所示，這是本文所提的三維FFBP算法的原理圖，為了方便考慮，只選取了虛擬面陣的一列與圖像空間的一列，并將初始脈沖點(diǎn)設(shè)置為4個，圖像空間像素點(diǎn)為1個。

首先將整個孔徑劃分成4個子孔徑，每個脈沖點(diǎn)坐標(biāo)記為n^（0）_i（i=0，1，2，3，）其中上角標(biāo)0的意思為第一級，此時假設(shè)圖像空間僅為一個像素點(diǎn)，其坐標(biāo)記為p^（0）_k=（p^（0）_n，p^（0）_m，p^（0）_t），其中p^（0）_n、p^（0）_m、p^（0）_t分別為航跡向、垂直航跡向、距離向像素點(diǎn)坐標(biāo)，k，n，m，t=0。因此，按照BP算法計(jì)算出此時每個脈沖點(diǎn)與每個像素點(diǎn)的距離并相干疊加，具體表達(dá)式為

d［n^（0）_i，p^（0）_n，p^（0）_m］=n^（0）_i－p^（0）_k2（6）

得到的初始圖像表達(dá)式為

I（d［n^（0）_i，p^（0）_n，p^（0）_m］）=

∑s（τ，t_i，l_x）expj4πλd［n^（0）_i，p^（0）_n，p^（0）_m］（7）

式中：l_x為垂直航跡向第x個陣元的坐標(biāo);d［n^（0）_i，p^（0）_n，p^（0）_m］為每個脈沖點(diǎn)到圖像空間中每個像素點(diǎn)的距離。

將得到的初始圖像作為下一級的輸入圖像，本文所提三維FFBP算法可以認(rèn)為是通過對輸入圖像做相位補(bǔ)償來得到本級圖像。上一級的每個脈沖點(diǎn)到本級圖像空間中的每個像素點(diǎn)的距離為

d［n^（0）_j，p^（1）_n，p^（1）_m］=n^（0）_j－p^（1）_k2（8）

這里是算法關(guān)鍵的步驟，圖3（c）意味著在中間階段，本級的每個子孔徑包含上一級兩個脈沖點(diǎn)，也可以認(rèn)為是本級的每個子孔徑包含兩個父孔徑。同時，計(jì)算出本級每個脈沖點(diǎn)到本級圖像空間每個像素點(diǎn)的距離為

d［n^（1）_i，p^（1）_n，p^（1）_m］=n^（1）_i－p^（1）_k2（9）

可以得到距離誤差表達(dá)式為

Δr_j=d［n^（0）_j，p^（1）_n，p^（1）_m］－d［n^（1）_i，p^（1）_n，p^（1）_m］（10）

得到的本級圖像表達(dá)式為

I（d［n^（1）_i，p^（1）_n，p^（1）_m］）=

∑ni∈njI（d［n^（0）_j，p^（1）_n，p^（1）_m］）expj4πλΔr_j（11）

接下來是最核心的思想，本算法假設(shè)了一個近似條件為

I（d［n^（s－1）_j，p^（s）_n，p^（s）_m］）≈

I（d［n^（s－1）_j，p^（s－1）_n/2，p^（s－1）_m/4）（12）

對于二維FFBP算法來說，其只考慮了沿航跡向的近似，這個假設(shè)意味著利用上一級生成的圖像表示下一級圖像的輸入值，假設(shè)沿航跡向每一次增加了兩個像素點(diǎn)，即假設(shè)第一級圖像有兩個像素點(diǎn)，而第二級圖像有4個像素點(diǎn)。根據(jù)式（12），將第一級生成圖像中的第1個像素值分配給第二級輸入圖像中的第1個和第2個像素值。同理，將第一級生成圖像中的第2個像素值分配給第二級輸入圖像中的第3個和第4個像素值，后續(xù)迭代也以此類推下去。對于三維FFBP算法來說，本文將垂直航跡向每一次產(chǎn)生的4個像素點(diǎn)，即沿垂直航跡向的上一級圖像的第1個像素值分配給下一級輸入圖像的前4個像素值并作為初始值，以此類推。

式（12）中這種近似不存在誤差的條件為

d［n^（s-1）_j，p^（s）_n，p^（s）_m］=d［n^（s－1）_j，p^（s－1）_n/2，p^（s－1）_m/4］ （13）

當(dāng)然，當(dāng)不滿足式（13）時，式（12）中的近似關(guān)系就會出現(xiàn)誤差，具體在SAS上如何設(shè)置這種距離誤差，本文認(rèn)為需要在實(shí)驗(yàn)中不斷調(diào)試以達(dá)到比較好的成像效果，這種距離誤差越小越好。

通過上面的不斷迭代，當(dāng)所有子孔徑最終融合成一整個孔徑時便恢復(fù)了圖像，最后一個階段與前面的迭代階段有所不同，最終生成的圖像表達(dá)式為

I（d［n^（s）_c，p^（s）_n，p^（s）_m］）=

∑nc∈njI（d［n^（s－1）_j，p^（s）_n，p^（s）_m］）·

expj4πλd［n^（s）_c，p^（s）_n，p^（s）_m］（14）

最終的距離誤差表達(dá)式為

d［n^（s）_c，p^（s）_n，p^（s）_m］=n^（s）_c－p^（s）_k2（15）

式中：n^（s）_c為整個完整孔徑中心坐標(biāo);d［n^（s）_c，p^（s）_n，p^（s）_m］為完整孔徑中心到圖像空間每個像素點(diǎn)的距離。

以上便是完整的三維FFBP算法流程，通過子孔徑劃分，使得原始三維BP算法的計(jì)算復(fù)雜度由O（N⁵）降低到了O（N3log₂Nlog₂N）。為了便于分析，這里假定垂直航跡向陣元數(shù)、航跡向脈沖數(shù)、垂直航跡向、航跡向、距離向像素點(diǎn)數(shù)均為N，每次迭代合并因子為2。這只是在算法層面上的加速，下面本文將利用GPU實(shí)現(xiàn)并行計(jì)算，相比于串行計(jì)算，并行計(jì)算可以更加快速地處理三維成像中數(shù)據(jù)量大的問題。

2.2 利用GPU加速三維FFBP算法

GPU與CPU設(shè)備的架構(gòu)是迥異的，CPU的設(shè)計(jì)是用來運(yùn)行少量比較復(fù)雜的任務(wù)，GPU的設(shè)計(jì)則是用來運(yùn)行大量比較簡單的任務(wù)，CPU的設(shè)計(jì)主要針對大量離散而不相關(guān)任務(wù)的系統(tǒng)。而GPU的設(shè)計(jì)主要針對解決那些可以分解成成千上萬個小塊并可獨(dú)立運(yùn)行的問題^［25］。如圖4所示，圖中N代表第N個線程塊與線程。一般來說，一塊GPU中一般包含多個流多處理器（stream multiprocessors， SM），每個SM上有許多個塊，每個塊中又有許多線程（一般至少為128個），正是因?yàn)镚PU可以同時調(diào)度如此多的線程數(shù)量，所以GPU更適合進(jìn)行大規(guī)模并行計(jì)算。需要注意的是，GPU設(shè)置塊與線程數(shù)量時不是隨意選取的，對于不同版本的GPU，所支持的SM數(shù)量與每個SM中擁有的塊數(shù)量也是不同的，同時每個塊中最多擁有1 024個線程。在處理大規(guī)模并行計(jì)算問題時，需要綜合考慮硬件性能，塊的數(shù)量與每個塊中的線程數(shù)量并不一定是越多越好，具體參數(shù)如何設(shè)置需要根據(jù)設(shè)備性能與算法需求，通過塊與線程數(shù)量的不斷修改與調(diào)試，來達(dá)到最優(yōu)的計(jì)算效率。一般來講，塊的數(shù)量要遠(yuǎn)多于SM的數(shù)量，這使得設(shè)備能夠得到充分并行，同時要保證每個塊中線程的數(shù)量為32的整數(shù)倍。

英偉達(dá)為GPU增加了一個易用的編程接口，也就是所謂的統(tǒng)一計(jì)算架構(gòu)（compute unified device architecture， CUDA），本文使用了CUDA對所提三維FFBP算法進(jìn)行硬件層面的加速，并通過仿真實(shí)驗(yàn)驗(yàn)證了GPU在實(shí)現(xiàn)并行計(jì)算時的巨大潛力與優(yōu)勢。

3 非理想航跡下的誤差分析

在第2節(jié)中，假設(shè)了聲納基陣做勻速直線運(yùn)動。然而在真實(shí)情況下，SAS基陣載體的運(yùn)動情況往往會受到許多因素的干擾（如風(fēng)浪、潮汐等）從而偏離航跡。在150 kHz的工作頻率下，聲波波長約為1 cm。超過1 mm的擺動誤差會導(dǎo)致SAS圖像散焦^［26］。一般來講，為了實(shí)現(xiàn)圖像聚焦而不出現(xiàn)明顯模糊，運(yùn)動擾動必須保持在波長的1/8甚至1/16以內(nèi)^［2731］。超過這個范圍可能會導(dǎo)致圖像退化，影響成像效果。下面將進(jìn)行必要的理論推導(dǎo)，聲納基陣載體的非理想航跡的距離表達(dá)式為

d₁=（x+Δu－t_x）²+（y－t_y）²+（z－t_z）²（16）

聲納基陣載體的理想航跡距離表達(dá)式為

d₂=（x－t_x）²+（y－t_y）²+（z－t_z）²（17）

式中：x、y、z為陣元的垂直航跡向、航跡向、距離向坐標(biāo);Δu為垂直航跡向偏移距離;t_x、t_y、t_z為目標(biāo)垂直航跡向、航跡向、距離向坐標(biāo)。

因此，運(yùn)動偏差會導(dǎo)致延時的變化，附加的延時表達(dá)式為

Δτ=2（d₁－d₂）speed（18）

當(dāng)這種附加延時Δτ超出一定范圍時，如果不進(jìn)行補(bǔ)償，便會影響成像質(zhì)量。

本文模擬了一種搖擺情況，即在垂直航跡方向添加了一個周期內(nèi)的正弦擾動，這種正弦擾動的表達(dá)式為

Δu=Asin θ（19）

式中：A為擾動幅度;θ為擾動角度。

一般來講，偏移距離的限制條件為

Δu≤λ8（20）

本文算法對點(diǎn)目標(biāo)成像進(jìn)行了一定的加窗處理，因此在偏移量不是很大（1個波長范圍內(nèi)）時，可以有效地抑制由基陣載體的不規(guī)則運(yùn)動帶來的干擾。圖5為1個波長范圍內(nèi)的擺動誤差。當(dāng)運(yùn)動誤差較大時，就必須考慮運(yùn)動補(bǔ)償，通過對回波信號的補(bǔ)償來提高圖像的質(zhì)量以及相位信息的準(zhǔn)確性。

4 仿真實(shí)驗(yàn)結(jié)果與分析

本文仿真了平面陣下的三維FFBP算法，然后利用C++對點(diǎn)目標(biāo)進(jìn)行了驗(yàn)證，最后考慮了非理想航跡的情況，并在仿真軟件中利用本文算法對點(diǎn)目標(biāo)的峰值旁瓣比指標(biāo)進(jìn)行了對比分析。

4.1 FFBP算法距離誤差分析與計(jì)算效率對比

表1為三維成像聲納的設(shè)備參數(shù)，本文采用虛擬面陣的陣元數(shù)為16×64（垂直航跡向×航跡向），圖像空間設(shè)置的像素點(diǎn)為128×128×128（垂直航跡向×航跡向×距離向），圖像空間的場景大小為-1.6～1.6 m（垂直航跡向），-0.1～0.1 m（航跡向），19.68～20.32 m（距離向）。

表2給出了FFBP算法初始時的距離誤差參數(shù)設(shè)置，其中垂直航跡向距離誤差ΔR_x、航跡向距離誤差ΔR_y由式（1）計(jì)算得到，ΔR_x為λ/50，ΔR_y為λ/256，父、子孔徑間的平均距離誤差Δr_mean、父、子孔徑間的最大距離誤差Δr_max由式（10）計(jì)算得到，可以看到誤差均在1個波長范圍內(nèi)。

在考慮計(jì)算效率時，首先在CUDA中對GPU進(jìn)行了塊與線程數(shù)量的配置，通過不斷調(diào)整塊與線程數(shù)量以便得到最優(yōu)的計(jì)算效率。然后，分別在不同仿真環(huán)境下利用OpenMP（open multi-processing）、在vs2019中利用CUDA進(jìn)行了運(yùn)行時間的對比。所計(jì)算的運(yùn)行時間包括采集回波后的脈沖壓縮、距離向升采樣以及成像算法，筆記本的配置為Intel（R） Core（TM） i7-10870H CPU@ 2.21 GHz，顯卡型號為GTX1650TI。

表3給出了GPU在對脈沖壓縮、距離向升采樣進(jìn)行處理時的參數(shù)配置，可以發(fā)現(xiàn)，通過調(diào)整塊與線程數(shù)量，其計(jì)算效率也有所不同。表3證明了在設(shè)置GPU參數(shù)時，并不一定是選取的塊數(shù)量越多，其計(jì)算效率就越快，當(dāng)然GPU中選取的塊數(shù)量也不能太少，否則設(shè)備將無法得到充分并行，表3的最后一行中出現(xiàn)了錯誤，這是因?yàn)榭偟木€程數(shù)量超過1 024個，GPU中每個塊的線程數(shù)量最多為1 024個。維度方向?yàn)榇怪焙桔E向×航跡向×距離向。

表4是GPU對三維FFBP算法處理時的參數(shù)配置，當(dāng)選取4個塊、每個塊中有32個線程時其計(jì)算效率最好，而當(dāng)只有1個塊時，三維FFBP算法的計(jì)算效率明顯下降。正如前文所言，選取的塊數(shù)量太少導(dǎo)致設(shè)備沒有得到充分并行。需要說明的是，GPU在對三維FFBP算法進(jìn)行處理時，只能沿距離向進(jìn)行并行計(jì)算，這是因?yàn)槿SFFBP是一種迭代算法，可并行計(jì)算的前提是循環(huán)之間互不影響，所以三維FFBP算法沿航跡向與垂直航跡向是通過子孔徑劃分來提高計(jì)算效率的。其中，沿垂直航跡向迭代次數(shù)為4，每次迭代合并2個子孔徑，生成4個像素點(diǎn);沿航跡向迭代次數(shù)為6，每次迭代合并2個子孔徑，生成2個像素點(diǎn)，維度方向?yàn)榫嚯x向。

表5分別對比了在仿真中運(yùn)行三維BP算法、三維FFBP算法以及在vs2019中運(yùn)行三維FFBP算法、三維FFBP算法利用CPU加速（即采用OpenMP、三維FFBP算法）利用GPU加速（即采用CUDA）的計(jì)算時間，其中在CUDA中針對脈沖壓縮與距離向升采樣最終選取的塊數(shù)量為16×32×512，線程數(shù)量為1×2×64，對于三維FFBP算法最終選取的塊數(shù)量為4，線程數(shù)量為32。

可以看到在合成孔徑長度只有20 m時，三維BP算法就已經(jīng)需要計(jì)算接近270 s，難以實(shí)現(xiàn)實(shí)時性成像。利用本文所提出的三維FFBP算法將計(jì)算時間縮短為38 s，相比于三維BP算法，計(jì)算時間縮短至約為此前的1/7。將仿真代碼改寫成C++代碼后計(jì)算時間又得到了加速，這是因?yàn)镃++語言更加高效，在處理循環(huán)結(jié)構(gòu)時表現(xiàn)更好。接著，本文利用CPU對三維FFBP算法進(jìn)行了加速，在vs2019中利用OpenMP對三維FFBP算法進(jìn)行CPU層面的并行計(jì)算，vs2019中編譯器支持OpenMP2.0版本，可以看到計(jì)算時間縮短到了26 s。最后，本文利用GPU加速三維FFBP算法，對三維FFBP算法進(jìn)行GPU層面的并行計(jì)算，利用CUDA實(shí)現(xiàn)三維FFBP算法后，最終計(jì)算時間僅需2 s左右，這表明利用GPU實(shí)現(xiàn)的大規(guī)模并行計(jì)算可以很好地解決三維成像中的大數(shù)據(jù)量問題。

4.2 三維FFBP算法成像結(jié)果

圖6是三維FFBP算法對點(diǎn)目標(biāo)在3個維度上的成像效果，可以看到本文所提算法對點(diǎn)目標(biāo)在3個維度上均得到了很好的聚焦效果。

圖7是各個方向的剖面圖，本文所提算法在各個方向均進(jìn)行加窗處理，所用窗函數(shù)為Kaiser窗。從圖7可以看到3個方向的峰值旁瓣比均滿足成像需求。雖然加窗會使得主瓣展寬，降低分辨率，但通過主瓣的適當(dāng)展寬來換取旁瓣的降低是可以接受的。

如圖8所示，經(jīng)過三維FFBP算法后點(diǎn)目標(biāo)位置在（0.012 5，0.000 62，20.022 5）m，由此可以驗(yàn)證本文所提算法精度很高，對目標(biāo)的定位很準(zhǔn)確。

4.3 BP算法與FFBP算法對比

本節(jié)將從BP算法與FFBP算法的計(jì)算效率與成像質(zhì)量兩方面進(jìn)行對比分析。

4.3.1 計(jì)算效率對比

相比FFBP算法，BP算法在考慮脈沖壓縮與距離向升采樣時沒有任何區(qū)別。區(qū)別在于BP算法由于不涉及迭代過程，因此采用GPU時可以實(shí)現(xiàn)3個維度的并行計(jì)算。同樣地，表6給出了GPU運(yùn)行三維BP成像算法時的參數(shù)配置。維度方向?yàn)榇怪焙桔E向×航跡向×距離向。

通過表6可知，在進(jìn)行三維BP成像算法時選取128×128×4個塊以及1×1×32個線程。表7給出了各個階段的計(jì)算時間對比以及總時間對比。

通過表7可以看到，在同樣使用GPU時，F(xiàn)FBP算法要比BP算法快了1.08 s，這對于本就只有2～3 s的成像時間而言，計(jì)算效率的提升是很可觀的。

4.3.2 成像質(zhì)量對比

FFBP算法是一種以犧牲一定成像質(zhì)量來換取計(jì)算效率提升的算法，因此本節(jié)將會對比FFBP算法與BP算法沿3個維度方向的峰值旁瓣比。圖9是BP算法沿3個方向的剖面圖?？梢钥吹?，與圖7中的FFBP算法相比，各個方向的峰值旁瓣比均有所降低。

圖10是BP算法在3個維度上對點(diǎn)目標(biāo)的成像效果，可以看到聚焦效果很好。

4.4 非理想航跡下的成像指標(biāo)分析

聲納基陣平臺在實(shí)際運(yùn)動時會受到風(fēng)浪、潮汐等因素影響導(dǎo)致基陣載體偏離航跡。因此，本文模擬了一種搖擺情況，即在垂直航跡方向上添加一個周期內(nèi)的不同幅度的正弦擾動。表8是在仿真環(huán)境中不同正弦擾動幅度下的3個維度峰值旁瓣比指標(biāo)，可以看到，當(dāng)正弦擾動幅度為0.2 m（2個波長）時，距離向峰值旁瓣比略高于-13 dB;當(dāng)擾動幅度到0.4 m（4個波長）時，距離向峰值旁瓣比已經(jīng)遠(yuǎn)高于-13 dB，此時已經(jīng)無法滿足成像需求，必須考慮采用運(yùn)動補(bǔ)償來保證成像質(zhì)量。

5 結(jié)束語

由于BP算法的計(jì)算復(fù)雜度較高，在考慮三維成像時，難以實(shí)現(xiàn)實(shí)時性成像。因此，本文提出了一種三維FFBP算法，算法計(jì)算時間從263 s降低到了38 s，計(jì)算效率相比于三維BP算法提升了近7倍。將仿真語言改寫成C++語言后，計(jì)算效率有小幅度提升。同時，在vs2019中對比了CPU與GPU加速效果，發(fā)現(xiàn)CPU加速需要26 s，而GPU加速僅需2 s左右。這是因?yàn)镚PU更適合進(jìn)行大規(guī)模并行計(jì)算，在考慮三維成像時，可以充分發(fā)揮GPU的優(yōu)勢。同時，在CUDA中設(shè)置塊與線程數(shù)量時，需要綜合考慮硬件性能，才能達(dá)到最優(yōu)的計(jì)算效率。通過子孔徑劃分與GPU加速，經(jīng)過對點(diǎn)目標(biāo)的驗(yàn)證，本文提出的三維FFBP算法在成像質(zhì)量與精度上均有良好效果。接著，對比了三維FFBP算法與三維BP算法同時采用GPU加速后的計(jì)算效率以及成像質(zhì)量。三維FFBP算法在計(jì)算效率上提高了1.08 s，效率提升對于僅有2 s左右的成像時間來說是顯著的。與此同時，三維FFBP算法沿3個維度方向上的峰值旁瓣比雖然有所升高，但均在較理想范圍內(nèi)。最后，考慮了非理想航跡的情況，即沿垂直航跡方向添加了不同程度的正弦擾動。在幅度不超過0.1 m（1個波長范圍內(nèi)）時，本文算法仍然適用，但當(dāng)這種擾動超出一定范圍時，就必須考慮采用運(yùn)動補(bǔ)償來保證一定的成像質(zhì)量。

參考文獻(xiàn)

［1］ 李寶奇， 黃海寧， 劉紀(jì)元， 等. 基于改進(jìn)CycleGAN的光學(xué)圖像遷移生成水下小目標(biāo)合成孔徑聲納圖像算法研究［J］. 電子學(xué)報， 202 49（9）： 17461753.

LI B Q， HUANG H N， LIU J Y， et al. Research on the algorithm of generating underwater small target synthetic aperture sonar images based on improved CycleGAN［J］. Acta Electronica Sinica， 202 49（9）： 17461753.

［2］ TINH N D， KHANH T D. A new back-projection algorithm in frequency domain for multi-receiver synthetic aperture sonar［C］∥Proc.of the 8th NAFOSTED Conference on Information and Computer Science， 2021： 3944.

［3］ THOMAS B， HUNTER A， DUGELAY S. Phase wrap error correction by random sample consensus with application to synthetic aperture sonar micronavigation［J］. IEEE Journal of Oceanic Engineering， 2020， 46（1）： 221235.

［4］ GERG I D， MONGA V. Deep multi-look sequence processing for synthetic aperture sonar image segmentation［J］. IEEE Trans.on Geoscience and Remote Sensing， 2023， 61： 4200915.

［5］ TINH N D， KHANH T D. A new imaging geometry model for multi-receiver synthetic aperture sonar considering variation of the speed of sound in seawater［J］. IEEE Trans.on Smart Processing and Computing， 202 10（4）： 302308.

［6］ ZHANG X B， YANG P X. Back projection algorithm for multi-receiver synthetic aperture sonar based on two interpolators［J］. Journal of Marine Science and Engineering， 202 10（6）： 718.

［7］ FOCSA A， ANGHEL A， DATCU M. Accelerated back-projection SAR processor on arbitrary elliptical imaging grid with azimuth spectrum unfolding［J］. IEEE Geoscience and Remote Sensing Letters， 2023， 20： 4008005.

［8］ 李彥君， 劉佳， 徐秋鋒. 基于切比雪夫擬合的BP自聚焦算法［J］. 系統(tǒng)工程與電子技術(shù)， 202 44（10）： 30203028.

LI Y J， LIU J， XU Q F. BP self-focusing algorithm based on Chebyshev fitting［J］. Systems Engineering and Electronics， 202 44（10）： 30203028.

［9］ SOMMER A， OSTERMANN J. Backprojection subimage autofocus of moving ships for synthetic aperture radar［J］. IEEE Trans.on Geoscience and Remote Sensing， 2019， 57（11）： 83838393.

［10］ ZHANG B J， XU G， ZHOU R， et al. Multi-channel back-projection algorithm for mmWave automotive MIMO SAR imaging with Doppler-division multiplexing［J］. IEEE Journal of Selected Topics in Signal Processing， 202 17（2）： 445457.

［11］ NEWBOLD H G. Method and apparatus for measuring x-or radiation absorption or transmission at plural angles and analyzing the data［P］. U.S.Patent： 3778614， 19731211.

［12］ YEGULALP A F. Fast backprojection algorithm for synthetic aperture radar［C］∥Proc.of the IEEE Radar Conference， 1999： 6065.

［13］ ULANDER L M H， HELLSTEN H， STENSTROM G. Synthetic-aperture radar processing using fast factorized back-projection［J］. IEEE Trans.on Aerospace and Electronic Systems， 2003， 39（3）： 760776.

［14］ MOON K， LONG D G. A new factorized backprojection algorithm for stripmap synthetic aperture radar［J］. Scientific Research， 2013， 4（1）： 4256.

［15］ MOTA D， CRUZ H， MIRANDA P R， et al. Onboard processing of synthetic aperture radar backprojection algorithm in FPGA［J］. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 202 15： 36003611.

［16］ GONG H， LIU Y， CHEN X Y， et al. Scene optimization of GPU-based back-projection algorithm［J］. The Journal of Supercomputing， 2023， 79（4）： 41924214.

［17］ XU Y， ZHANG Z， CHEN L Y， et al. The adaptive streaming SAR back-projection algorithm based on half-precision in GPU［J］. Electronics， 202 11（18）： 2807.

［18］ 韋順軍， 蒲羚， 張曉玲， 等. 復(fù)雜軌跡合成孔徑雷達(dá)后向投影算法圖像流GPU成像［J］. 電訊技術(shù)， 2016， 56（8）： 879886.

WEI S J， PU L， ZHANG X L， et al. GPU imaging of complex trajectory synthesis aperture radar backprojection algorithm image flow［J］. Telecommunications Technology， 2016， 56（8）： 879886.

［19］ PARK J， TANG P T P， SMELYANSKIY M， et al. Efficient backprojection-based synthetic aperture radar computation with many-core processors［J］. Scientific Programming， 2013， 21（34）： 165179.

［20］ HETTIARACHCHI D L N， BALSTER E J. Fixed-point processing of the SAR back-projection algorithm on FPGA［J］. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 202 14： 1088910902.

［21］ CAO Y， GUO S C， JIANG S， et al. Parallel optimisation and implementation of a real-time back projection （BP） algorithm for SAR based on FPGA［J］. Sensors， 202 22（6）： 2292.

［22］ WIELAGE M， CHOLEWA F， FAHNEMANN C， et al. High performance and low power architectures： GPU vs. FPGA for fast factorized backprojection［C］∥Proc.of the 5th International Symposium on Computing and Networking， 2017： 351357.

［23］ 李浩林， 張磊， 楊磊， 等. 一種改進(jìn)的快速分解后向投影 SAR 成像算法［J］. 電子與信息學(xué)報， 2013， 35（6）： 14351441.

LI H L， ZHANG L， YANG L， et al. An improved fast factorized back-projection SAR imaging algorithm［J］. Journal of Electronics amp; Information Technology， 2013， 35（6）： 14351441.

［24］ BANK S M. Studies in high resolution synthetic aperture sonar［D］. London： University College London， 2003.

［25］ SHANE C. CUDA并行程序設(shè)計(jì)： GPU編程指南［M］. 蘇統(tǒng)華， 李東， 李松澤， 等， 譯. 北京： 機(jī)械工業(yè)出版社， 2014.

SHANE C. CUDA programming： a developer’s guide to parallel computing with GPUs［M］. SU T H， LI D， LI S Z， et al， Trans. Beijing： China Machine Press， 2014.

［26］ SHIPPEY G， JONSSON M， PIHL J "N B. Position correction using echoes from a navigation fix for synthetic aperture sonar imaging［J］. IEEE Journal of Oceanic Engineering， 2009， 34（3）： 294306.

［27］ BELLETTIN A， PINTO M A. Theoretical accuracy of synthetic aperture sonar micronavigation using a displaced phase-center antenna［J］. IEEE Journal of Oceanic Engineering， 200 27（4）： 780789.

［28］ BONEFANT W W. Interferometic synthetic aperture sonar processing［M］. Atlanta： Georgia Institute of Technology， 1999.

［29］ THOMAS B. Phase preserving 3D micro-navigation for interferometric synthetic aperture sonar［D］. Bath： University of Bath， 2020.

［30］ CALLOW H J. Signal processing for synthetic aperture sonar image enhancement［D］. Christchurch： University of Canterbury， 2003.

［31］ HANSEN R E， SAEBO T O， GADE K， et al. Signal processing for UAV based interferometric synthetic aperture sonar［C］∥Proc.of the Oceans， 2003： 24382444.

作者簡介

陶鴻博（1998—），男，助理工程師，碩士，主要研究方向?yàn)樗曅盘柼幚?、高分辨聲成像算法?/p>

張東升（1980—），男，研究員，博士，主要研究方向?yàn)樗曅盘柼幚怼⒏叻直鏈y深側(cè)掃聲納。

黃 勇（1965—），男，教授，博士，主要研究方向?yàn)樗曅盘柼幚?、三維高分辨成像。