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

本文提出了一種GPU環(huán)境下的多子陣合成孔徑聲納實時CS成像算法。首先分析了多子陣合成孔徑聲納的CS算法成像基本流程，然后對算法中的多子陣回波向單子陣回波轉(zhuǎn)化，方位向和距離向FFT變換和點乘運算的實現(xiàn)過程進(jì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，并進(jìn)行陣元信號的重排。式中V表示聲基陣載體速度，PRI表示脈沖重復(fù)間隔，D表示陣元尺寸。

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

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

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

5）在距離時域、方位頻域進(jìn)行方位脈沖壓縮和剩余相位補(bǔ)償。

多子陣合成孔徑聲納CS成像算法從計算角度看主要包含復(fù)數(shù)點乘、FFT和IFFT運算。但隨著合成孔徑聲納技術(shù)的發(fā)展，用于成像的原始數(shù)據(jù)維數(shù)越來越大，這三類基本運算的速度嚴(yán)重影響著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ù)，然后采用單子陣合成孔徑成像算法進(jìn)行處理。多子陣SAS回波向單子陣回波轉(zhuǎn)化需要經(jīng)過兩個步驟:1）有效數(shù)據(jù)截?。?）固定相位補(bǔ)償。由于實際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ù)進(jìn)行有效截取可以直接采用高效內(nèi)存拷貝函數(shù)memcpy完成。固定相位補(bǔ)償分塊示意如圖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放置，因此固定相位補(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)系后，只需在時域中與相位補(bǔ)償因子相乘即可，其中 f0為中心頻率，c為聲速，r為距離，v載體速度，Δhi為第i個接收陣中心與發(fā)射陣中心的距離。

圖2 固定相位補(bǔ)償分塊示意圖

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變換要求，必須先進(jìn)行矩陣轉(zhuǎn)置操作，將方位向數(shù)據(jù)變?yōu)閮?nèi)存中的連續(xù)存儲形式。在采用CUDA進(jìn)行矩陣轉(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ù)來進(jìn)行方位向FFT正/逆變換，這時一次FFT變換的點數(shù)為方位向點數(shù)Na，連續(xù)變換次數(shù)為距離向點數(shù)Nr。最后再進(jìn)行一次矩陣轉(zhuǎn)置，完成方位向FFT變換。

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

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

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

為了驗證本文所提CS并行成像算法性能，在如下計算環(huán)境中進(jì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)境下分別進(jìn)行了單核CS算法成像、多核CS算法成像和GPU環(huán)境下的CS算法成像。

成像試驗中選用的原始數(shù)據(jù)是在千島湖進(jìn)行干涉合成孔徑聲納海試樣機(jī)試驗中獲取的，對應(yīng)的幅度信息如圖3（a）所示，最近采樣距離為51m，最遠(yuǎn)采樣距離231m，其距離向和方位向點數(shù)分別為9600和2880。從圖3（a）中可以看出遠(yuǎn)距離處信號回波較弱。試驗系統(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中可以看出，采用單核和多核計算時，不需要進(jìn)行數(shù)據(jù)傳輸。對于原始測試回波，采用CS單核成像算法，總的計算時間為27827ms。在共享內(nèi)存環(huán)境下，采用OpenMP并行后，總的計算時間降為7111ms，加速比為3.91。采用GPU計算，雖然需要進(jìn)行數(shù)據(jù)的上傳和下載操作，但因數(shù)據(jù)傳輸帶寬大，對于測試數(shù)據(jù)僅為202ms。由于采用了CUDA中高效的FFT變換函數(shù)庫，成像部分的FFT運算時間大幅度減小，此外算法中的點乘運算效率也大幅度提高。最后基于GPU的并行算法成像時間降為494.7ms，與單核計算相比，加速比達(dá)到56.25，滿足SAS實時成像需要。

5 結(jié)語

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