基于GPU的天線組陣信號時延補償方法

毛飛龍, 焦義文, 馬 宏,*, 韓久江, 高澤夫, 李 超, 李 冬

(1. 航天工程大學(xué)電子與光學(xué)工程系, 北京 101400; 2. 國防科技大學(xué)電子科學(xué)學(xué)院, 湖南 長沙 410073)

0 引 言

隨著航天技術(shù)與深空探測技術(shù)的迅猛發(fā)展,各國對空間的探測向著更深、更遠(yuǎn)的范圍拓展。探測距離將從目前月球探測的40×104km拓展到火星探測的4×108km,以及木星探測的10×108km,這將帶來嚴(yán)重的傳輸損耗,進(jìn)而導(dǎo)致單天線接收信噪比(signal to noise ratio, SNR)越來越低。天線組陣[1]信號合成是解決單天線接收SNR低的一種重要方法,是指用多個天線構(gòu)成天線陣列對同一信源的信號進(jìn)行接收,之后對不同天線接收到的信號進(jìn)行相干合成,達(dá)到提高接收信噪比的目的。天線組陣合成技術(shù)的本質(zhì)是在消除不同天線信號間時延差和相位差后進(jìn)行相干相加,由于各天線間的噪聲(例如大氣衰減、雨衰)是隨機不相關(guān)的,理論上使用N個天線構(gòu)成的陣列合成信號的SNR是單天線接收信號的N倍[2]。

每個天線接收到的航天器信號存在由空間幾何路徑、大氣傳輸介質(zhì)和設(shè)備鏈路造成的時延差異。由空間幾何路徑產(chǎn)生的幾何時延和由大氣傳輸介質(zhì)產(chǎn)生的傳播時延可以根據(jù)先驗知識進(jìn)行標(biāo)校[3]。設(shè)備鏈路是變頻系統(tǒng),會產(chǎn)生設(shè)備時延和相位,需要通過設(shè)備鏈路標(biāo)校技術(shù)進(jìn)行標(biāo)校。受大氣湍流和設(shè)備溫漂等因素影響,采用上述方法標(biāo)校后的殘余時延和相位仍然存在較大的波動[4],因此信號合成系統(tǒng)還需要以自適應(yīng)的方式,通過相關(guān)處理來估計并實時修正殘余時延和相位。

天線組陣合成系統(tǒng)是一個數(shù)據(jù)密集且高度并行、邏輯復(fù)雜、對實時性要求高的系統(tǒng)。隨著大規(guī)模天線組陣合成系統(tǒng)的天線數(shù)量越來越多、接收信號的帶寬越來越大,使用單獨的圖形處理器(graphic processing unit, GPU)已不能滿足運算要求。因此,需要尋求能夠進(jìn)行大規(guī)模數(shù)據(jù)運算以及高效并行的體系架構(gòu)和硬件平臺。

近年來,隨著高性能計算的發(fā)展,GPU從專用于圖像領(lǐng)域的處理器逐漸向通用并行計算平臺轉(zhuǎn)變。GPU已發(fā)展成為一種高度并行化、多線程、多核的通用計算設(shè)備,具有杰出的計算能力和極高的存儲器帶寬,并被越來越多地應(yīng)用于圖形處理之外的計算領(lǐng)域。GPU的運算核心數(shù)遠(yuǎn)多于CPU,更適合于數(shù)據(jù)密集型計算的并行加速處理。NVIDIA于2007年推出了計算統(tǒng)一設(shè)備架構(gòu)(compute unified device architecture, CUDA),簡化了GPU系統(tǒng)的開發(fā)流程,使得GPU通用計算技術(shù)在信號處理領(lǐng)域得到了更為廣泛的應(yīng)用。目前,基于GPU的信號處理系統(tǒng)具有強大的并行處理能力,能夠在短時間內(nèi)通過并行處理完成大量數(shù)據(jù)的運算,在GPU資源得到充分利用時,可以實現(xiàn)對信號的實時處理。

因此,基于GPU的信號處理技術(shù)成為眾多領(lǐng)域的熱點,如射電天文[5-6]、雷達(dá)[7-8]、無線通信[9-10]、人工智能[11-12]等。本文主要研究天線組陣信號時延補償方法在GPU平臺下的實現(xiàn)與并行優(yōu)化。

1 天線組陣信號時延補償方法研究

1.1 天線組陣技術(shù)概括

天線組陣技術(shù)經(jīng)過幾十年的發(fā)展,逐漸形成了五種合成方案[13]:全頻譜合成(full-spectrum combining, FSC)、基帶合成(baseband combining, BC)、符號流合成(symbol-stream combining, SSC)、載波組陣合成(carrier arraying combining, CAC)、復(fù)符號合成(complex-symbol combining, CSC)。除FSC方法,上述方案都屬于基于載波跟蹤技術(shù)的合成方法,需要在陣元載波鎖定后合成,不適用于當(dāng)前大規(guī)模的深空天線組陣信號合成系統(tǒng)。

在天線組陣信號合成系統(tǒng)中,由于天線間下行鏈路的差異性,不同天線接收到同一信源的信號存在時延差,進(jìn)而導(dǎo)致合成的信號SNR達(dá)不到最低解調(diào)門限。為了保證能夠正確接收和解調(diào)信源的信號,就需要在合成信號之前對時延差進(jìn)行補償。根據(jù)時延差和系統(tǒng)采樣率的關(guān)系,時延補償可分為整數(shù)時延補償和小數(shù)時延補償。天線組陣中最常用的合成方法就是全頻譜合成,該方法直接在中頻進(jìn)行合成,通過相位差估計[14]算法進(jìn)行自適應(yīng)時延與相位的估計并進(jìn)行補償。圖1所示為天線組陣信號全頻譜合成示意圖[15]。

圖1 天線組陣信號FSC示意圖Fig.1 Schematic diagram of FSC of antenna array signals

為便于分析,假設(shè)天線信號為正弦信號。天線1為參考天線,天線2為待修正天線。天線1與天線2接收到的信號為

(1)

(2)

理想情況下,θ21=θ2(t)-θ1(t),即θ21為兩天線間的相位差值。之后利用ejθ21對天線2進(jìn)行補償:

s2(t)C=Asin[ωct+θ2(t)]·ejθ21

(3)

此時,兩天線的相位已對齊,將s1(t)與s2(t)C相干相加可得

s(t)=Asin[ωct+θ1(t)]+Asin[ωct+θ2(t)]·ejθ21

(4)

理論上來說,合成的信號s(t)的SNR為s1(t)的2倍,但實際傳輸信道中存在著各種噪聲干擾,使得合成信號的SNR小于s(t)SNR的50%[15]。

由圖1可知,天線組陣信號FSC方法中的兩個關(guān)鍵技術(shù)分別為時延估計和時延補償。時延估計算法主要通過對天線信號進(jìn)行互相關(guān)、對互相關(guān)譜求反正切等運算實現(xiàn)。典型的時延估計算法有Simple、Sumple[16-17]、Eigen[18]等算法。在準(zhǔn)確地估計出時延差和相位差后,就要進(jìn)行高精度時延補償,其主要目的是消除兩路信號的時延差,將兩路信號的時延與相位對齊,進(jìn)而進(jìn)行相干相加。因此,在估計性能達(dá)標(biāo)的情況下,合成性能優(yōu)劣取決于能否將估計所得的時延差和相位差高精度地補償?shù)礁鱾€天線。

1.2 典型時延補償方法

文獻(xiàn)[19-20]借鑒Sumple算法的思路,提出了一種以準(zhǔn)合成輸出信號作為虛擬參考信號的時延校準(zhǔn)方法。該方法將待修正天線以外的其他所有天線的信號作為虛擬參考天線,提高了虛擬參考信號的SNR,從而有效提高了時延估計性能。然而,時域方法將信號以整體形式(全帶寬)進(jìn)行時延估計和補償,存在兩方面缺點:① 相關(guān)時,需要循環(huán)卷積,運算量大;② 補償時,采用數(shù)字延遲線,只能進(jìn)行整數(shù)時延補償,補償精度受采樣率限制。美國航空航天局(National Aeronautics and Space Administration, NASA)深空網(wǎng)在其全頻譜處理陣(full spectrum processor array, FSPA)設(shè)備中,使用了一種基于上下邊帶的殘余時延和相位差聯(lián)合估計方法[3,14,21]。該方法首先通過信號濾波器提取各天線信號的上邊帶和下邊帶,然后利用相位估計算法分別估計待校準(zhǔn)天線與參考天線上下邊帶的相位差,最終兩相位點構(gòu)成的直線的斜率即為殘余時延,截距即為載波相位差。該方法的優(yōu)勢在于運算量小,而且將時延估計問題轉(zhuǎn)化為相位估計問題,在低SNR條件下,各邊帶可以采用高性能相位估計算法來提高估計精度。然而,該方法存在兩個問題:① 由于相位估計算法只能提取初相差,當(dāng)上下邊帶頻率差較大或殘余時延較大時,易出現(xiàn)模糊問題[22-23];② 時延補償精度受采樣率限制[24]。

由上述分析可知,時域合成方法均存在補償精度受采樣率限制的問題。即使假設(shè)時延估計誤差為零,補償后的小數(shù)時延也將導(dǎo)致頻帶內(nèi)相位滑動,從而惡化合成性能[3,25-26]。尤其當(dāng)帶寬較大時,上述方法不再適用。因此,引入更合理的小數(shù)時延補償方法是信號合成的重點和難點問題之一。

1.3 整數(shù)時延補償

整數(shù)時延通常而言就是信號的采樣點超前或滯后,現(xiàn)有的整數(shù)時延補償方法通過移動信號采樣點來實現(xiàn)。在現(xiàn)場可編程門陣列(field programmable gate array, FPGA)平臺下,整數(shù)時延補償可通過移位寄存器實現(xiàn)。在Matlab讀數(shù)據(jù)處理的方式下,可通過移動讀文件的起始地址實現(xiàn)整數(shù)時延補償。

在相位差估計模塊,需要在頻域進(jìn)行信號的共軛相乘、累加平均、求反正切等運算。其中,累加平均就是將M組512點(假設(shè)快速傅里葉變換(fast Fourier transform, FFT)點數(shù)為512)的頻域信號進(jìn)行能量累積(FFT積分),以提高估計的準(zhǔn)確性。在FPGA平臺下和Matlab仿真中,都能保證每次移動采樣點后,有足夠多的信號(大于等于積分長度M×512)以進(jìn)行后續(xù)的處理。

考慮到GPU在實時處理時以信號數(shù)據(jù)塊的方式進(jìn)行傳輸,下面分析典型的整數(shù)時延補償方法(移動信號采樣點),以及以這種方式實現(xiàn)的可行性,及其存在的問題。在GPU平臺下,要滿足信號合成的實時性,就必須在特定時間內(nèi)處理完成一定的數(shù)據(jù)量(例如在1 s之內(nèi)完成與采樣率數(shù)值大小一致的數(shù)據(jù)量),這就需要對數(shù)據(jù)流進(jìn)行分塊處理。圖2所示為GPU數(shù)據(jù)流分塊處理示意圖。

圖2 GPU數(shù)據(jù)流分塊處理示意圖Fig.2 Schematic diagram of GPU data flow block processing

每個數(shù)據(jù)塊的大小固定為M×512,每次迭代通過cudaMemcpy分別將天線1和參考天線的一個數(shù)據(jù)塊傳入GPU進(jìn)行運算。第1次迭代時,時延補償初始化為0,因此整數(shù)時延與小數(shù)時延均為0,不需要對天線1和天線2進(jìn)行采樣點偏移。第2次迭代時,利用第1次迭代估計的時延值進(jìn)行補償,首先將時延取整得到整數(shù)時延。若天線1相比參考天線時延超前D個采樣點,則天線1需要向前移動D個采樣點以實現(xiàn)整數(shù)時延補償,移動之后天線1只剩下(M×512-D)個有效數(shù)據(jù)。之后進(jìn)行時延差估計,需要對天線1和參考天線做FFT運算。此時,參考天線有M×512個有效數(shù)據(jù),而天線1只有(M×512-D)個有效數(shù)據(jù)。若FFT運算時點數(shù)不夠,則需自動進(jìn)行補零操作。因此,參與FFT運算的天線1信號并非完整的M×512數(shù)據(jù)點數(shù),而是末尾補了D個0的數(shù)據(jù)塊。同理,若天線2相比參考天線時延滯后D個采樣點,則天線2需要向后移動D個采樣點以實現(xiàn)整數(shù)時延補償,移動之后天線2只剩(M×512-D)個有效數(shù)據(jù)。在FFT運算時同樣需要進(jìn)行補零操作,這就會導(dǎo)致時延差估計時無法利用所有數(shù)據(jù)點的信息,估計的時延差誤差較大,進(jìn)而無法實現(xiàn)高精度時延補償。

1.4 小數(shù)時延補償

寬帶信號時域合成方法都采用整數(shù)時鐘移位的時延補償方法,存在時延補償精度低的問題。帶寬越大,時延補償精度要求越高。為保證合成損失小于0.1 dB,殘余時延引起的整帶內(nèi)的相位差必須小于20°。若帶寬為500 MHz,則殘余時延必須小于0.11 ns。設(shè)采樣時鐘為1 280 MHz,該補償方法的補償精度為0.39 ns(半個時鐘周期),已不能滿足要求,必須采用更高精度的小數(shù)時鐘時延補償方法。

小數(shù)時延就是時延數(shù)值中的小數(shù)部分,而小數(shù)時延無法通過移動信號的采樣點實現(xiàn),因此需要尋找其他方法進(jìn)行小數(shù)時延補償[27]。常見的小數(shù)時延補償主要分為兩類,分別是時域補償和頻域補償。二者的主要區(qū)別在于是否在頻域進(jìn)行小數(shù)時延補償。在FPGA平臺下,時域補償法主要通過各種最小誤差準(zhǔn)則逼近理想系統(tǒng)獲得有限沖擊響應(yīng),主要包括基于最小均方誤差(minimum mean square error, MMSE)準(zhǔn)則濾波法、拉格朗日插值法和基于Farrow結(jié)構(gòu)的濾波器組方法等[27]。而分?jǐn)?shù)時延數(shù)字濾波器大大增加了硬件開銷和實現(xiàn)的復(fù)雜度。

1.5 小 結(jié)

現(xiàn)有時延補償方法的不足主要體現(xiàn)在以下幾個方面。

(1) 在GPU平臺下,為了實現(xiàn)實時處理,每次處理的數(shù)據(jù)塊的大小都是固定的。而整數(shù)時延補償通過移動信號的起始采樣點來實現(xiàn),這將導(dǎo)致數(shù)據(jù)塊的有效點數(shù)減少,在之后的FFT運算時會自動進(jìn)行補零操作,從而使得估計時延差時沒有充分利用所有的數(shù)據(jù),導(dǎo)致時延差估計的精度降低,無法實現(xiàn)高精度時延補償。

(2) 時域小數(shù)時延補償通常采用分?jǐn)?shù)時延數(shù)字濾波器,而這種方法的實現(xiàn)復(fù)雜度高且增加了硬件開銷,實現(xiàn)信號的實時處理難度較大。

2 基于GPU的天線組陣信號時延補償方法

2.1 總體方案

本文對基于GPU的時延補償方法展開研究,圖3為基于GPU的天線組陣信號時延補償方法的總體方案框圖。N個天線信號經(jīng)采集卡輸出為N路并行的數(shù)字信號,其中某個天線固定為參考天線。首先,進(jìn)行整數(shù)時延補償,之后對所有天線信號進(jìn)行子帶拆分,此時信號已轉(zhuǎn)為頻域,之后進(jìn)行小數(shù)時延補償和相位補償。補償完成后,將信號分別傳輸至?xí)r延差估計模塊和子帶合成模塊。時延差估計模塊賦值計算出天線間的時延差和相位差,并將待補償?shù)闹祩鬏斨料乱淮蔚小Ｗ訋Ш铣赡K負(fù)責(zé)將天線信號合成并輸出至子帶重構(gòu)模塊,進(jìn)而輸出合成信號。

圖3 基于GPU的天線組陣信號時延補償方法總體方案框圖Fig.3 Block diagram of the overall scheme of time delay compensation method of antenna array signal based on GPU

2.2 整數(shù)時延補償

2.2.1 天線滯后參考天線

在整數(shù)時延補償方面,依舊采用移動信號采樣點的方法實現(xiàn)。針對數(shù)據(jù)塊的有效數(shù)據(jù)點數(shù)不足和進(jìn)行FFT時自動補零的問題,采用數(shù)據(jù)塊重疊保留的思想進(jìn)行改進(jìn)。圖4為天線滯后參考天線時的整數(shù)時延補償示意圖,具體實現(xiàn)步驟如下。

圖4 天線滯后參考天線時整數(shù)時延補償示意圖Fig.4 Schematic diagram of integer time delay compensation when the antenna lags behind the reference antenna

步驟 1第1次迭代開始,參考天線第1個數(shù)據(jù)塊A1(M×512個數(shù)據(jù))和天線X的第1個數(shù)據(jù)塊B1,通過PCIE總線由CPU內(nèi)存?zhèn)鬏數(shù)紾PU片上內(nèi)存。在這兩個數(shù)據(jù)塊前分別添加1個長度為N、數(shù)值為0的數(shù)據(jù)塊,并將后續(xù)處理的地址指向數(shù)據(jù)塊的最前端(補0數(shù)據(jù)塊的首地址),然后開始時延補償。由于初始化時延補償值為0,因此整數(shù)時延與小數(shù)時延均為0,此次迭代不需要進(jìn)行整數(shù)時延和小數(shù)時延補償,FFT運算的指針不需要偏移,仍為整個數(shù)據(jù)塊的最前端。FFT需要的數(shù)據(jù)點數(shù)為M×512,因此天線1和參考天線的數(shù)據(jù)塊的末尾N個數(shù)據(jù)未參與本次計算。之后進(jìn)行時延估計,假設(shè)天線X(B1數(shù)據(jù)塊)滯后參考天線(A1數(shù)據(jù)塊)D1個采樣點,則時延估計后的整數(shù)時延為D1。將參考天線數(shù)據(jù)塊末尾N個數(shù)據(jù)組成數(shù)據(jù)塊A1_tail,將天線X數(shù)據(jù)塊末尾(D1+N)個數(shù)據(jù)組成數(shù)據(jù)塊B1_tail,將A1_tail、B1_tail數(shù)據(jù)塊拷貝至下一個積分區(qū)間。

步驟 2第2次迭代開始,參考天線第2個數(shù)據(jù)塊A2(M×512個數(shù)據(jù))和天線X的第2個數(shù)據(jù)塊B2通過符合總線和接口標(biāo)準(zhǔn)(peripheral component interface express, PCIE)總線由CPU內(nèi)存?zhèn)鬏數(shù)紾PU片上內(nèi)存。將上一個積分區(qū)間拷貝而來的A1_tail添加至A2的前端構(gòu)成A2_new,將B1_tail添加至B2的前端構(gòu)成B2_new,并將后續(xù)處理的地址指向A2_new、B2_new數(shù)據(jù)塊的最前端。A1_tail為A1末尾的N個數(shù)據(jù),B1_tail為B1數(shù)據(jù)塊末尾的(D1+N)個數(shù)據(jù),此時B2_new(天線X)超前A2_new(參考天線)D1個采樣點,實現(xiàn)了天線X的整數(shù)時延補償。同時,保證了數(shù)據(jù)的連續(xù)性,在之后的FFT運算時數(shù)據(jù)點數(shù)大于M×512,不需要進(jìn)行自動補零。該方法在后續(xù)估計時延差時充分利用了所有的數(shù)據(jù),提高了估計的精度。假設(shè)天線X(B2_new數(shù)據(jù)塊)滯后參考天線(A2_new數(shù)據(jù)塊)D2個采樣點,則時延估計后的整數(shù)時延為D2。將參考天線數(shù)據(jù)塊(A2_new)末尾的N個數(shù)據(jù)組成數(shù)據(jù)塊A2_tail,將天線X數(shù)據(jù)塊(B2_new)末尾的(D2+D1+N)個數(shù)據(jù)組成數(shù)據(jù)塊B2_tail,將A2_tail、B2_tail數(shù)據(jù)塊拷貝至下一個積分區(qū)間。

步驟 3第3次迭代,之后每次迭代與第2次迭代的處理過程相同。

2.2.2 天線超前參考天線

天線超前參考天線時的處理方法與滯后參考天線時的處理方法類似。圖5為天線超前參考天線時的整數(shù)時延補償示意圖。

圖5 天線超前參考天線時整數(shù)時延補償示意圖Fig.5 Schematic diagram of integer time delay compensation when the antenna is ahead of the reference antenna

步驟 1第1次迭代開始,參考天線第1個數(shù)據(jù)塊A1(M×512個數(shù)據(jù))和天線X的第1個數(shù)據(jù)塊B1通過PCIE總線由CPU內(nèi)存?zhèn)鬏數(shù)紾PU片上內(nèi)存。在這兩個數(shù)據(jù)塊前分別添加1個長度為N、數(shù)值為0的數(shù)據(jù)塊,并將后續(xù)處理的地址指向數(shù)據(jù)塊的最前端(補0數(shù)據(jù)塊的首地址),然后開始時延補償。由于初始化時延補償值為0,因此整數(shù)時延與小數(shù)時延均為0,此次迭代不需要進(jìn)行整數(shù)時延和小數(shù)時延補償,FFT運算的指針不需要偏移,仍為整個數(shù)據(jù)塊的最前端。FFT需要的數(shù)據(jù)點數(shù)為M×512,因此天線1和參考天線的數(shù)據(jù)塊的末尾N個數(shù)據(jù)未參與本次計算。之后進(jìn)行時延估計,假設(shè)天線X(B1數(shù)據(jù)塊)超前參考天線(A1數(shù)據(jù)塊)D1個采樣點,則時延估計后的整數(shù)時延為-D1。將參考天線數(shù)據(jù)塊末尾N個數(shù)據(jù)組成數(shù)據(jù)塊A1_tail,將天線X數(shù)據(jù)塊末尾(-D1+N)個數(shù)據(jù)組成數(shù)據(jù)塊B1_tail,將A1_tail、B1_tail數(shù)據(jù)塊拷貝至下一個積分區(qū)間。

步驟 2第2次迭代開始,參考天線第2個數(shù)據(jù)塊A2(M×512個數(shù)據(jù))和天線X的第2個數(shù)據(jù)塊B2通過PCIE總線由CPU內(nèi)存?zhèn)鬏數(shù)紾PU片上內(nèi)存。將上一個積分區(qū)間拷貝而來的A1_tail添加至A2的前端構(gòu)成A2_new,將B1_tail添加至B2的前端構(gòu)成B2_new,并將后續(xù)處理的地址指向A2_new、B2_new數(shù)據(jù)塊的最前端。A1_tail為A1末尾的N個數(shù)據(jù),B1_tail為B1數(shù)據(jù)塊末尾的(-D1+N)個數(shù)據(jù),此時B2_new(天線X)滯后A2_new(參考天線)D1個采樣點,實現(xiàn)了天線X的整數(shù)時延補償。同時保證了數(shù)據(jù)的連續(xù)性,在之后的FFT運算時數(shù)據(jù)點數(shù)大于M×512,不需要進(jìn)行自動補零。該方法在后續(xù)估計時延差時充分利用了所有的數(shù)據(jù),提高了估計的精度。假設(shè)天線X(B2_new數(shù)據(jù)塊)超前參考天線(A2_new數(shù)據(jù)塊)D2個采樣點,則時延估計后的整數(shù)時延為-D2。將參考天線數(shù)據(jù)塊的(A2_new)末尾N個數(shù)據(jù)組成數(shù)據(jù)塊A2_tail,將天線X數(shù)據(jù)塊(B2_new)的末尾(-D2-D1+N)個數(shù)據(jù)組成數(shù)據(jù)塊B2_tail,將A2_tail、B2_tail數(shù)據(jù)塊拷貝至下一個積分區(qū)間。

步驟 3第3次迭代,之后每次迭代與第2次迭代的處理過程相同。

2.3 小數(shù)時延與殘余相位補償

小數(shù)時延補償若采用可調(diào)分?jǐn)?shù)時延數(shù)字濾波器,則會增加硬件開銷和實現(xiàn)復(fù)雜度,信號實時處理的實現(xiàn)難度較大。而將信號轉(zhuǎn)為頻域,通過頻域進(jìn)行補償[28-30]可以有效解決上述問題。

首先,將寬帶信號拆分成多個子帶[31-32],子帶間隔決定了最大無模糊時延。例如,當(dāng)采樣率為1 280 MHz,子帶數(shù)為256時,子帶間隔為5 MHz,則最大無模糊時延為200 ns,大大減小了系統(tǒng)標(biāo)校的精度要求。其次,頻域合成方法[33-34]殘余的小數(shù)時延以及相位在頻域進(jìn)行補償,通過對各子帶進(jìn)行線性相位偏移來實現(xiàn)。這樣,由殘余小數(shù)時延引起的整帶內(nèi)相位滑動被限制在子帶范圍內(nèi),大大減小了小數(shù)時延引起的合成損失。仍以上個例子進(jìn)行說明,殘余小數(shù)時延最大值為0.39 ns(半個時鐘周期),在500 MHz帶寬內(nèi)引起的相位差約為70°,而在5 MHz子帶內(nèi)引起的相位差僅為0.7°,其引起的合成損失基本可以忽略。此外,頻域合成方法[35-36]通過對所有子帶相位差進(jìn)行線性最小二乘擬合,可以利用所有子帶的相位差信息,因此其殘余時延和相位的估計精度也更高。

因此,在小數(shù)時延補償方面,本文采用頻域補償?shù)姆椒?。在頻域?qū)ψ訋Х謩e進(jìn)行補償主要有以下優(yōu)勢:① 相比在全帶寬進(jìn)行小數(shù)時延補償,通過將各子帶中心頻率獨立進(jìn)行相位對齊,可將由小數(shù)時延引起的整帶內(nèi)相位滑動限制在子帶范圍內(nèi),大大減小了小數(shù)時延引起的合成損失;② 相比時域小數(shù)時延補償,具有結(jié)構(gòu)簡單、運算量小的優(yōu)勢。在頻域進(jìn)行小數(shù)時延補償,采用線性相位偏移實現(xiàn),對各個子帶獨立進(jìn)行相位偏移,可使用GPU并發(fā)進(jìn)行加速,大大提高了小數(shù)時延補償?shù)膶崟r性。

3 基于GPU的天線組陣信號時延補償方法

3.1 實現(xiàn)流程

圖6為基于GPU的時延補償方法流程圖。初始化時延和相位補償值為0,每次迭代過程先對信號進(jìn)行補償,再對殘余時延和殘余相位進(jìn)行估計,最后對時延和相位補償值進(jìn)行更新。其中,delay_now和phase_now為當(dāng)前時刻的殘余時延和殘余相位的估計值,mu為調(diào)整步長。因此,在每次迭代過程中,時延和相位補償值都由上次迭代過程估計得出,通過閉環(huán)反饋支路進(jìn)行更新。也就是說,前一個數(shù)據(jù)塊估計出的殘余時延與殘余相位都為下一個數(shù)據(jù)塊的時延和相位的補償值進(jìn)行了預(yù)報。

圖6 基于GPU的天線阻陣時延補償方法流程圖Fig.6 Flowchart of time delay compensation method of antenna array based on GPU

首先利用上一次迭代估計出的殘余時延計算出整數(shù)時延和小數(shù)時延;然后使用基于數(shù)據(jù)塊重疊保留的方法進(jìn)行整數(shù)時延補償;之后利用小數(shù)時延和殘余相位,通過多項式求解函數(shù)ployval反推出各個子帶需要補償?shù)南辔恢怠Ｈ缓?通過線性相位偏移進(jìn)行小數(shù)時延和殘余相位補償。

基于GPU的天線組陣時延補償方法實現(xiàn)流程圖如圖7所示。

圖7 基于GPU的天線組陣時延補償方法實現(xiàn)流程Fig.7 Implementation process of GPU-based antenna array time delay compensation method

本文設(shè)計的基于GPU的天線組陣信號時延補償方法基于GPU+CPU異構(gòu)處理平臺開發(fā),以CPU為控制器,以GPU為協(xié)處理器。由于GPU采用并行處理的架構(gòu)設(shè)計,其在進(jìn)行大量并行數(shù)據(jù)處理時計算性能明顯優(yōu)于CPU[37]。

3.2 時延補償算法GPU實現(xiàn)

根據(jù)第2.2節(jié)對算法的具體描述,本文設(shè)計的整數(shù)時延補償?shù)腃UDA代碼如算法1所示,主要通過使用CUDA中的數(shù)據(jù)傳輸函數(shù)cudaMemcpyDeviceToDevice實現(xiàn)。delay為上一次迭代估計出的時延補償值,step為當(dāng)前已迭代次數(shù)。signal為本次迭代由CPU傳至GPU內(nèi)存中的初始數(shù)據(jù)塊,signal_last為上一次迭代參與運算的數(shù)據(jù)塊,signal_now為本次迭代待參與運算的數(shù)據(jù)塊。N為重疊保留數(shù)據(jù)塊初始化大小,參與運算數(shù)據(jù)塊的大小固定為nsize,其中參考天線數(shù)據(jù)塊命名均以ref結(jié)尾。

算法 1 基于數(shù)據(jù)塊重疊保留的整數(shù)時延補償輸入: delay、step、signal、signal_last、N、nsize輸出: signal_nowInteger_delay=int (round (delay [step-1]));cudaMemcpy(signal_now, signal_last+(nsize-N)-Integer_delay, (N+Integer_delay)×sizeof(float2), cudaMemcpyDeviceToDevice);cudaMemcpy(signal_now_ref, signal_last_ref+(nsize-N), N×sizeof(float2), cudaMemcpyDeviceToDevice);cudaMemcpy(signal_now+(N+Integer_delay), signal, (nsize-(N+Integer_delay))×sizeof(float2), cudaMem-cpyDeviceToDevice);cudaMemcpy(signal_now_ref+N, signal_ref, (nsize-N)×sizeof(float2), cudaMemcpyDeviceToDevice)

步驟 1使用取整函數(shù)計算整數(shù)時延;

步驟 2將上一次迭代參與運算的數(shù)據(jù)塊(signal_last)末尾的(N+Integer_delay)個數(shù)據(jù)拷貝至signal_now的首地址處;

步驟 3將signal拷貝至signal_now的(N+Integer_delay)地址處。

線性相位偏移函數(shù)如算法2所示,rec為待處理數(shù)據(jù),dst為已處理數(shù)據(jù)。phase為相位偏移值,處理數(shù)據(jù)長度為length。

算法 2 線性相位偏移函數(shù)輸入: rec、phase、length 輸出: dst__global__ void process(float2×rec, float×phase, float2×dst, int length){int tid=threadIdx.x+blockIdx.x×blockDim.x;if (tid>=length)return;float Phase=phase[tid% NFFT];dst[tid].x=rec[tid].x×cos(Phase)+rec[tid].y×sin(Phase);dst[tid].y=-rec[tid].x×sin(Phase)+rec[tid].y×cos(Phase);}

小數(shù)時延與殘余相位補償如算法3所示。delay、phase分別為上一次迭代估計出的時延補償值、相位補償值,采樣率為Fre_sample,當(dāng)前已迭代次數(shù)為step,子帶頻率范圍為f_span,FFT點數(shù)為NFFT,參與運算數(shù)據(jù)塊大小固定為nsize。signal_FFT為待補償?shù)念l域信號,補償完成的頻域信號為signal_process。

步驟 1計算小數(shù)時延fraction_delay;

步驟 2以fraction_delay和phase為輸入,利用ployval反推出各個子帶需要補償?shù)南辔恢祊hase_offset;

步驟 3調(diào)用process函數(shù)對轉(zhuǎn)至頻域的信號signal_FFT進(jìn)行線性相位偏移,即小數(shù)時延與殘余相位補償。

4 實驗驗證與分析

4.1 系統(tǒng)及核心參數(shù)設(shè)計

利用相位校正(phase calibration, PCAL)信號進(jìn)行性能分析,PCAL信號的頻率特性具有梳狀頻譜特性和線性相位特性。通過在接收機鏈路輸出端提取PCAL信號的自相位譜,可以得到PCAL信號經(jīng)過整個傳輸鏈路的附加時延、相位以及非線性相位失真。

選用表1和表2所示的GPU和CPU仿真平臺和具體仿真環(huán)境參數(shù),仿真參數(shù)設(shè)置如表3所示。

表1 GPU參數(shù)

表2 CPU參數(shù)

表3 仿真實驗參數(shù)

本文中所有幅度均為仿真信號的歸一化幅度。圖8所示為時延補償結(jié)果對比圖。圖8(a)為天線1與參考天線初始信號的時域波形細(xì)節(jié)圖(為方便比較,選取其中一段數(shù)據(jù))。從圖8(a)可以看出,天線1與參考天線的時延為2個采樣點,信號未對齊。圖8(b)為現(xiàn)有方法(即FFT時點數(shù)不夠的自動補零方法)迭代30次時的補償細(xì)節(jié)圖,圖8(c)為本文提出的基于數(shù)據(jù)塊重疊保留的整數(shù)時延補償?shù)?0次時的補償細(xì)節(jié)圖。從圖8(b)和圖8(c)可以看出,兩種方法的天線1與參考天線都實現(xiàn)了對齊,且合成信號的幅度約為參考天線的2倍。從時域波形無法看出二者之間的區(qū)別。

圖8 時延補償結(jié)果對比圖Fig.8 Comparison diagram of time delay compensation results

圖9為時延補償誤差對比圖,即迭代第30次時兩種方法的補償誤差(天線1補償后信號幅度與參考天線之差)對比圖,現(xiàn)有方法的補償誤差在10-1量級,而本文方法的補償誤差在10-7量級。

圖9 時延補償誤差對比圖Fig.9 Time delay compensation error comparison chart

圖10為兩種方法的時延補償值隨迭代次數(shù)變化的曲線。delayi=delayi-1+mu×delaynow,其中delaynow為本次迭代計算出的時延,delayi-1為上一次迭代時的時延補償值,delayi為本次更新后的時延補償值。在時延真值為2個采樣點的情況下,20次迭代后,兩種方法的時延補償值都收斂于某一值。本文方法的時延補償值收斂于真值2.000,現(xiàn)有方法的時延補償值收斂于2.003。這是由于FFT自動補零導(dǎo)致參與運算的有效數(shù)據(jù)點數(shù)減少,導(dǎo)致時延差估計的精度降低,無法實現(xiàn)高精度時延估計與補償。

圖10 時延補償值對比圖Fig.10 Time delay compensation value comparison

在得到所有子帶的相位差后,需要利用線性最小二乘擬合方法來估計寬帶信號的殘余時延和相位。擬合直線的斜率即為殘余時延,截距即為零頻處的殘余相位。因此,擬合的精度直接影響估計的殘余時延和相位的精度,進(jìn)而影響補償對齊的精度。下面分析、對比兩種方法在最小二乘擬合時的殘差。

圖11為第1次迭代擬合殘差對比圖,圖11 (a)為現(xiàn)有方法的第1次迭代擬合殘差,擬合殘差為10-12量級。由于初始時延補償值為0,因此不需要進(jìn)行整數(shù)時延補償,在FFT時也不需要進(jìn)行補零操作,充分利用了所有信號數(shù)據(jù),因此其擬合殘差很小,擬合性能好。本文提出的方法如圖11(b)所示。在第1次迭代時,擬合殘差達(dá)10-5量級。由于初始化進(jìn)行了信號補零,導(dǎo)致第一次迭代時沒有充分利用所有信號數(shù)據(jù),擬合殘差較大。

圖11 第1次迭代擬合殘差對比圖Fig.11 Comparison of the first iteration fitting residual

圖12為第30次迭代(收斂時)的擬合殘差對比圖。現(xiàn)有方法除第1次迭代,擬合殘差均達(dá)到了10-5量級。由于未進(jìn)行數(shù)據(jù)的重疊保留處理,在FFT時需要進(jìn)行補零操作,因此其擬合殘差較大,擬合質(zhì)量差。而本文提出的方法在第1次迭代之后,擬合殘差為10-10量級。由于擬合前的數(shù)據(jù)連續(xù)且充分,擬合殘差小,擬合質(zhì)量高。

圖12 第30次迭代擬合殘差對比圖Fig.12 Comparison of the 30th iteration fitting residual

4.2 加速效果驗證與分析

利用NVIDIA提供的分析工具NVIDIA Nsight Systems軟件[37]詳細(xì)分析了算法運行過程中各模塊的耗時情況,算法運行時序圖如圖13所示。

圖13 算法運行時序圖Fig.13 Algorithm runtime diagram

之后對1 000次迭代的耗時進(jìn)行記錄,并取平均值,得到算法在GPU端運算與CPU端運算兩部分的具體耗時,如表4所示。其中,GPU端運算占55%,CPU端運算占45%,GPU端耗時略多于CPU端。之后對GPU端總耗時中數(shù)據(jù)傳輸與核函數(shù)實際執(zhí)行時間的占比進(jìn)行分析,數(shù)據(jù)傳輸包括三個方向:CPU到GPU;GPU到CPU;GPU到GPU。也即MemcpyHtoD、MemcpyDtoH、MemcpyDtoD,分析結(jié)果如表5和圖14所示。

表4 算法耗時對比表

表5 GPU端耗時對比表

圖14 GPU端耗時對比餅圖Fig.14 Time consumed comparison pie chart of GPU

由表5可以看出,在GPU端總耗時中,核函數(shù)實際執(zhí)行時間占總耗時的45.13%,數(shù)據(jù)傳輸時間占總耗時的54.87%?？梢?數(shù)據(jù)傳輸時間占比稍多于核函數(shù)實際執(zhí)行時間占比,因此算法還可通過GPU異步流并發(fā)的方式進(jìn)行進(jìn)一步優(yōu)化,隱藏部分?jǐn)?shù)據(jù)傳輸耗時,使核函數(shù)實際執(zhí)行的時間占比進(jìn)一步提高。

為了進(jìn)一步研究時延補償模塊及整個系統(tǒng)的耗時情況,針對1 s的數(shù)據(jù)進(jìn)行加速比分析和實時性分析。圖15為時延補償加速對比圖,其中,本文設(shè)計的基于GPU的時延補償方法模塊處理完成1 s的數(shù)據(jù)量共耗時約為6.2 ms。

圖15 時延補償加速對比圖Fig.15 Comparison chart of time delay compensation acceleration

CPU串行時延補償與并行算法的區(qū)別主要有兩點:其一在于整數(shù)時延補償?shù)姆绞?串行方法采用memcpy函數(shù)在CPU端進(jìn)行,并行算法采用cudaMemcpyDtoD函數(shù)進(jìn)行傳輸;其二在于串行方法的小數(shù)時延與相位補償沒有進(jìn)行并行優(yōu)化,需要從頭至尾遍歷每一個值并乘以補償?shù)南辔弧Ｔ谔幚硐嗤瑪?shù)據(jù)量時耗時114 ms,加速比約為18。

圖16為信號合成實時性對比圖。在信號合成系統(tǒng)中,分別采用基于GPU的時延補償方法和CPU串行時延補償方法。除時延補償,其余模塊的處理方式相同。對比兩種方法在合成1 s數(shù)據(jù)的耗時可以看出,采用GPU的時延補償方法耗時約783 ms,采用CPU串行時延補償方法耗時約9 105 ms。實時性的要求是,處理1 s的數(shù)據(jù)耗時少于1 s (1 000 ms)。因此,本文基于GPU的時延補償方法在被應(yīng)用于合成系統(tǒng)中時達(dá)到了實時處理的能力,而在被應(yīng)用于現(xiàn)有的CPU串行時延補償方法的合成系統(tǒng)時遠(yuǎn)未達(dá)到實時性,只能應(yīng)用先采集數(shù)據(jù)、之后進(jìn)行非實時處理的方式。

圖16 信號合成實時性對比圖Fig.16 Comparison chart of signal synthesis real-time property

最后,對訪存帶寬需求進(jìn)行了分析,利用NVIDIA-smi命令監(jiān)視算法對GPU訪存的需求。從圖17可以看出,GPU的內(nèi)存總共約為80 GB,在文中表1的各項參數(shù)下,GPU已占用約為1.2 GB,該數(shù)值對相關(guān)的仿真實驗有一定的參考價值。

圖17 NVIDIA-smi命令監(jiān)視窗口Fig.17 NVIDIA-smi command watching window

5 結(jié) 論

針對大規(guī)模天線組陣合成系統(tǒng)天線數(shù)量越來越多、接收信號的帶寬越來越大,單獨使用CPU已不能滿足運算需求這一問題,本文研究了基于GPU的天線組陣信號時延補償方法,分析了算法并行的可行性,并設(shè)計了基于GPU的天線組陣信號時延補償方法,通過仿真對比實驗驗證了本文提出的并行算法的優(yōu)越性。本文所做主要工作如下。

(1) 針對GPU平臺以數(shù)據(jù)塊為單位進(jìn)行處理的方式,本文首先分析了典型的整數(shù)時延補償方法在GPU平臺實現(xiàn)的可行性,發(fā)現(xiàn)將現(xiàn)有方法直接應(yīng)用于GPU無法保證正確性。在此基礎(chǔ)上,設(shè)計了基于數(shù)據(jù)塊重疊保留的整數(shù)時延補償方法。

(2) 分析了傳統(tǒng)小數(shù)時延方法的優(yōu)劣,并設(shè)計了適合于GPU并行優(yōu)化的頻域小數(shù)時延補償方法。

(3) 通過實測數(shù)據(jù)對基于GPU的天線組陣信號時延補償方法進(jìn)行了實驗驗證。首先驗證了算法的正確性,證實了本文方法的補償誤差遠(yuǎn)低于現(xiàn)有的補償方法,在無噪聲的情況下達(dá)到10-7量級,證明了本文方法在GPU平臺上實現(xiàn)的正確性。之后,對基于GPU的天線組陣信號時延補償方法和CPU串行時延補償方法的實時處理性能進(jìn)行了對比。實驗表明,本文方法相比CPU串行方法有約18倍的加速比,利用本文方法可以實時地進(jìn)行信號合成。

本文的并行算法可以運用到未來的大規(guī)模天線組陣信號合成系統(tǒng)中,實現(xiàn)高速、準(zhǔn)確、實時的信號合成,有力支持未來的一系列深空探測任務(wù)。