基于GPU的線頻調(diào)變標(biāo)算法并行實現(xiàn)?

唐扶光 劉 婭 鐘何平

（1.武漢輕工大學(xué)電氣與電子工程學(xué)院 武漢 430023）（2.海軍工程大學(xué)海軍水聲技術(shù)研究所 武漢 430033）

1 引言

隨著合成孔徑聲納（SAS）技術(shù)的不斷發(fā)展，成像分辨率不斷提高，測繪帶寬度不斷加大，它們共同作用導(dǎo)致了用于成像的原始數(shù)據(jù)量顯著增加［1～3］，嚴重影響著合成孔徑聲納系統(tǒng)的實時成像。傳統(tǒng)的距離多普勒成像算法［4］需要進行大量的插值運算，效率低下，而線性調(diào)頻變標(biāo)（CS）算法［5～7］中只用了FFT運算和復(fù)數(shù)乘/加運算，極大提高了成像效率，并且許多學(xué)者已經(jīng)對其并行化方法進行了研究［8］。目前解決SAS的實時成像問題在硬件選擇上通常有兩種方案:1）專用硬件方案，其特點是計算速度快，價格昂貴，開發(fā)周期長，軟件拓展性差；2）通用計算方案，其特點是編程靈活，開發(fā)周期短，軟件拓展性好。目前的合成孔徑成像通用計算方案主要采用的是計算機集群［9～11］，并且在合成孔徑成像算法的并行化方面取得了較好效果。但集群運算的缺點是計算設(shè)備體積龐大，價格昂貴，功耗大。近年來出現(xiàn)的圖形處理器（GPU）具有強大的浮點計算能力，并且NVIDIA公司提供多種CUDA環(huán)境下的常用函數(shù)庫，大大降低了編程難度，使得GPU成為了一種通用的高密度計算設(shè)備。目前GPU已廣泛應(yīng)用于各種科學(xué)計算中，并且都取得了很好的加速比，這為合成孔徑聲納CS算法的實時成像提供了新的途徑［12～14］。

本文提出了一種GPU環(huán)境下的多子陣合成孔徑聲納實時CS成像算法。首先分析了多子陣合成孔徑聲納的CS算法成像基本流程，然后對算法中的多子陣回波向單子陣回波轉(zhuǎn)化，方位向和距離向FFT變換和點乘運算的實現(xiàn)過程進行了GPU環(huán)境下的并行化設(shè)計，實現(xiàn)了GPU環(huán)境下的并行CS成像算法。最后通過實際SAS回波的成像試驗，驗證了所提并行CS成像算法的性能。

2 多子陣合成孔徑聲納CS成像算法

多子陣合成孔徑聲納CS成像算法的基礎(chǔ)是CS算法，該算法思想是基于線性調(diào)頻信號尺度變換的特點，然后采用CS操作消除距離徙動隨距離的空變特性，使得不同距離上的點目標(biāo)距離徙動軌跡相同，再通過在二維頻域上乘一個線性相位因子完成所有的距離徙動校正，從而避免插值操作［15～16］。CS算法優(yōu)點是不需要耗時的插值運算，只需要快速傅立葉變換和復(fù)數(shù)相乘就能完成距離遷徙的校正，使得成像效率大大提高，而且便于并行運算，相位保真度高。存在不足的是，CS算法是基于線性調(diào)頻信號發(fā)展起來的，因此它僅限于此種信號形式?；贑S算法的多子陣SAS成像算法的基本流程如圖1所示，主要分為以下幾個步驟:

1）根據(jù)等效相位中心近似條件PRI*V=D*M/2，確定有效子陣個數(shù)M，并進行陣元信號的重排。式中V表示聲基陣載體速度，PRI表示脈沖重復(fù)間隔，D表示陣元尺寸。

2）對多子陣SAS由單發(fā)多收和非停走停模式引入的固定相位誤差采用等效相位中心近似進行補償。

3）在距離時域、方位頻域采用CS算法原理將不同距離上的目標(biāo)徙動軌跡校正到與參考距離上的目標(biāo)徙動相同。

4）將距離向變?yōu)轭l域，在二維頻域中同時完成距離壓縮和距離徙動校正。

5）在距離時域、方位頻域進行方位脈沖壓縮和剩余相位補償。

多子陣合成孔徑聲納CS成像算法從計算角度看主要包含復(fù)數(shù)點乘、FFT和IFFT運算。但隨著合成孔徑聲納技術(shù)的發(fā)展，用于成像的原始數(shù)據(jù)維數(shù)越來越大，這三類基本運算的速度嚴重影響著CS算法的實時性。GPU中的大量并行計算單元為復(fù)數(shù)點乘、FFT和IFFT快速計算提供了基礎(chǔ)。此外，NVIDIA公司提供的CUDA編程環(huán)境中包含高效的FFT和IFFT函數(shù)庫，在加速FFT和IFFT計算的同時，極大簡化了GPU應(yīng)用程序開發(fā)。

圖1 CS算法流程圖

3 CS成像算法的并行化

3.1 多子陣SAS回波向單子陣SAS回波并行轉(zhuǎn)化

多子陣SAS成像的最簡單方式就是將多子陣回波數(shù)據(jù)先等效為單子陣回波數(shù)據(jù)，然后采用單子陣合成孔徑成像算法進行處理。多子陣SAS回波向單子陣回波轉(zhuǎn)化需要經(jīng)過兩個步驟:1）有效數(shù)據(jù)截??；2）固定相位補償。由于實際SAS系統(tǒng)中子陣個數(shù)M與陣元長度D、平臺速度V、脈沖重復(fù)間隔PRI不能恰好滿足臨界條件:PRI*V=D*M/2，一般作業(yè)過程中要求實際平臺運動速度V0＜D*M/2/PRI，此時成像需要的有效陣元數(shù)M0=2V0*PRI/D。由于采集原始信號是先按陣元后按脈沖順序排列的，因此對連續(xù)脈沖數(shù)據(jù)進行有效截取可以直接采用高效內(nèi)存拷貝函數(shù)memcpy完成。固定相位補償分塊示意如圖2所示，M0表示有效子陣數(shù)，Nr表示距離向點數(shù)，K表示當(dāng)前處理數(shù)據(jù)塊的脈沖個數(shù)。完成有效數(shù)據(jù)截取后，原始數(shù)據(jù)是以二維形式存儲的，二維數(shù)據(jù)的行數(shù)為M0*K，列數(shù)為Nr。如果把原始信號抽象為三維Nr*M0*K形式，每一個脈沖內(nèi)數(shù)據(jù)一層，然后按照脈沖順序?qū)⒒夭▽盈B放置，因此固定相位補償線程分塊適合采用三維形式。假設(shè)單個三維線程塊的維數(shù)為BI*BJ*BK，則總的三維線程分塊數(shù)為TBI*TBJ*TBK，其中TBI=（Nr+BI-1）/BI，TBJ=（M0+BJ-1）/BJ，TBK=（K+BK-1）/BK。現(xiàn)假設(shè)當(dāng)前線程塊的索引為（Si，Sj，Sk），線程塊內(nèi)線程索引為（Ti，Tj，Tk），則當(dāng)前線程處理的數(shù)據(jù)行坐標(biāo)TNr=BI*Si+Ti，列坐標(biāo)TNa=（BJ*Sj+Tj）+M0*（BK*Sk+Tk）。建立線程索引與原始回波數(shù)據(jù)的對應(yīng)關(guān)系后，只需在時域中與相位補償因子相乘即可，其中 f0為中心頻率，c為聲速，r為距離，v載體速度，Δhi為第i個接收陣中心與發(fā)射陣中心的距離。

圖2 固定相位補償分塊示意圖

3.2 距離向和方位向FFT變換

高效FFT變換是合成孔徑聲納CS算法快速成像的基礎(chǔ)，CUDA中提供高效的FFT函數(shù)庫極大簡化了CS算法的實現(xiàn)過程。在CS算法實現(xiàn)過程中，距離向和方位向FFT變換都是一維變換，這一過程的實現(xiàn)是通過調(diào)用CUDA中cufftExecC2C（cufftHandle plan， cufftComplex*idata， cufftComplex*odata，int direction）函數(shù)實現(xiàn)的，其中plan用于指定FFT變換基本參數(shù)，包括數(shù)據(jù)類型、一次FFT變換點數(shù)和FFT連續(xù)變換次數(shù)，idata表示輸出數(shù)據(jù)地址，odata表示輸出結(jié)果地址，direction表示FFT正變換/反變換。對于距離向的FFT正/逆變換可以直接通過調(diào)用cufftExecC2C函數(shù)來完成，在構(gòu)造FFT變換句柄plan時，指定一次FFT變換的點數(shù)為距離向點數(shù)Nr，連續(xù)變換次數(shù)為方位向點數(shù)Na。

方位向FFT變換不能直接調(diào)用cufftExecC2C函數(shù)，因為在內(nèi)存中方位向數(shù)據(jù)不是連續(xù)存儲的。為了滿足cufftExecC2C函數(shù)FFT變換要求，必須先進行矩陣轉(zhuǎn)置操作，將方位向數(shù)據(jù)變?yōu)閮?nèi)存中的連續(xù)存儲形式。在采用CUDA進行矩陣轉(zhuǎn)置時，線程塊先采用合并方式讀入對應(yīng)的原始數(shù)據(jù)，在線程塊內(nèi)部完成局部數(shù)據(jù)塊的轉(zhuǎn)置操作，然后再采用合并方式將轉(zhuǎn)置后的結(jié)果寫入顯存。采用合并訪問方式可以有效避免數(shù)據(jù)訪問沖突，提高數(shù)據(jù)訪問速度。完成矩陣轉(zhuǎn)置后，調(diào)用cufftExecC2C函數(shù)來進行方位向FFT正/逆變換，這時一次FFT變換的點數(shù)為方位向點數(shù)Na，連續(xù)變換次數(shù)為距離向點數(shù)Nr。最后再進行一次矩陣轉(zhuǎn)置，完成方位向FFT變換。

3.3 CS算法中的點乘運算并行

將多子陣回波數(shù)據(jù)等效為單子陣回波以后，單陣CS算法在成像過程中共需三次進行點成運算:1）在距離時域-方位頻域（tr-fa）將不同距離的徙動曲線校正成一樣；2）在距離頻域-方位頻域（fr-fa）對不同距離的回波進行統(tǒng)一的距離徙動校正；3）在距離時域-方位頻域（tr-fa）進行剩余相位補償。在這三部分點乘運算中，1）和3）是在距離向上依次對每列數(shù)據(jù)進行點乘運算，2）是在方位向上依次對每行進行點乘運算，它們的共同點是對每行或每列進行點乘運算時不存在先后順序。對于CS算法中點乘運算采用一維線程組完成，假設(shè)線程組的維數(shù)為L×1，原始回波數(shù)據(jù)距離向和方位向維數(shù)分別為Nr×Na，在進行計算任務(wù)分配時，行坐標(biāo)為Tnr的列采用線程Tnr mod L計算點乘，列坐標(biāo)為Tna的行采用線程Tna mod L計算點乘，這樣可以最大限度地均衡使用計算資源。在進行點乘運算前，計算過程經(jīng)常使用的常量，包括距離向距離和頻率坐標(biāo)、方位向距離和頻率坐標(biāo)、CS因子等先在主機中完成計算，再上傳至顯存，這樣進行點乘計算時可以直接使用，避免重復(fù)計算，節(jié)省計算時間。

4 試驗結(jié)果與分析

為了驗證本文所提CS并行成像算法性能，在如下計算環(huán)境中進行了成像試驗:處理器Intel（R）Xeon（R）CPU X5650 2.67G（2處理器12核）；內(nèi)存48G；顯卡Tesla C2050；操作系統(tǒng)Windows 7專業(yè)版；軟件環(huán)境VS2008。為了充分驗證所提算法性能，在上述環(huán)境下分別進行了單核CS算法成像、多核CS算法成像和GPU環(huán)境下的CS算法成像。

成像試驗中選用的原始數(shù)據(jù)是在千島湖進行干涉合成孔徑聲納海試樣機試驗中獲取的，對應(yīng)的幅度信息如圖3（a）所示，最近采樣距離為51m，最遠采樣距離231m，其距離向和方位向點數(shù)分別為9600和2880。從圖3（a）中可以看出遠距離處信號回波較弱。試驗系統(tǒng)基本參數(shù)如表1所示。

表1 試驗系統(tǒng)基本參數(shù)

三種計算方式下的CS算法成像結(jié)果分別如圖3（b）～（d）所示，可以看出三種成像結(jié)果是一致的，成像結(jié)果非常清晰，滿足水下高分辨率成像需求，為水下小目標(biāo)探測提供了保證。

圖3 CS算法成像結(jié)果

表2 CS算法成像效率比較（ms）

三種計算方式下的CS算法效率比較如表2所示。從表2中可以看出，采用單核和多核計算時，不需要進行數(shù)據(jù)傳輸。對于原始測試回波，采用CS單核成像算法，總的計算時間為27827ms。在共享內(nèi)存環(huán)境下，采用OpenMP并行后，總的計算時間降為7111ms，加速比為3.91。采用GPU計算，雖然需要進行數(shù)據(jù)的上傳和下載操作，但因數(shù)據(jù)傳輸帶寬大，對于測試數(shù)據(jù)僅為202ms。由于采用了CUDA中高效的FFT變換函數(shù)庫，成像部分的FFT運算時間大幅度減小，此外算法中的點乘運算效率也大幅度提高。最后基于GPU的并行算法成像時間降為494.7ms，與單核計算相比，加速比達到56.25，滿足SAS實時成像需要。

5 結(jié)語

本文提出了一種基于GPU的實時合成孔徑聲納CS成像算法。給出了多子陣合成孔徑聲納CS成像算法的基本流程，單子陣SAS回波向單子陣回波轉(zhuǎn)換、方位向和距離向FFT變換和點乘運算的GPU并行實現(xiàn)方法。最后在GPU并行平臺上進行了CS算法成像試驗，驗證了所提并行算法性能，滿足SAS成像的實時性需求。