田海俊，徐 洋，陳學(xué)雷，李長華，吳峰泉，汪群雄，劉 勇

(1.三峽大學(xué)，湖北 宜昌 443002；2.中國科學(xué)院國家天文臺，北京 100012)

射電干涉陣GPU相關(guān)器的研究初探

田海俊1，2，徐 洋2，陳學(xué)雷2，李長華2，吳峰泉2，汪群雄1，劉 勇1

(1.三峽大學(xué)，湖北 宜昌 443002；2.中國科學(xué)院國家天文臺，北京 100012)

射電干涉儀陣列規(guī)模的不斷擴大使其觀測能力越來越強，但與之俱來的密集型數(shù)據(jù)的實時處理對傳統(tǒng)解決方案的性能和成本等帶來巨大的挑戰(zhàn)，針對該挑戰(zhàn)設(shè)計并驗證了GPU解決方案。首先著重分析了射電干涉儀的相關(guān)器在運算和傳輸?shù)确矫娴幕拘枨?，然后根?jù)射電干涉陣列信號的特征，嘗試了多種GPU相關(guān)器模型，通過將相關(guān)器的計算任務(wù)有效映射為GPU線程模型，顯著提高了GPU的利用率，使單GPU的實際計算性能達到了理論峰值性能的77%。最后以我國正研制的天籟計劃為依托，開展了多種關(guān)鍵技術(shù)的先導(dǎo)實驗，為進一步在GPU集群環(huán)境下針對大型射電干涉儀的需求研發(fā)一套同時兼?zhèn)涑杀镜?、性能高、功耗少等?yōu)勢的相關(guān)器打下了堅實的基礎(chǔ)。

射電干涉儀；FX相關(guān)器；GPU；實時處理系統(tǒng)

CN53－1189/P ISSN1672－7673

典型的干涉陣列的信號處理包括:時間序列信號的模數(shù)轉(zhuǎn)化(Digitization)、干涉顯示度的計算(Complex Visibility)和校準、噪音的消除、天圖的傅里葉重構(gòu)。其中干涉顯示度的計算又主要包括各路信號的傅里葉變換(F-engine)和交叉互關(guān)聯(lián)(X-engine)兩部分，它是射電干涉儀數(shù)據(jù)實時處理過程中最關(guān)鍵也是計算量最大的部分，它的計算量按陣元數(shù)目的平方急劇增長，即O(N2)。

目前國際上干涉儀的陣元個數(shù)正向數(shù)百乃至數(shù)千發(fā)展。世界上大多數(shù)已建成的天文干涉儀陣列，如美國的甚大陣列(Very Large Array，VLA)，印度的巨型米波射電望遠鏡(Giant Meter wave Radio Telescope，GMRT)等都由幾十個單元組成。但隨著干涉儀技術(shù)的不斷成熟，人們已經(jīng)開始籌劃或正在建越來越大的陣列，例如，我國正研制的“天籟計劃”，以及國際上正在籌建的平方千米陣列望遠鏡(Square Kilometer Array，SKA)預(yù)期都由數(shù)百乃至數(shù)千個單元組成，歐洲即將建成的低頻射電干涉陣列(Low Frequency Array，LOFAR)由20000個天線構(gòu)成的48個基站組成。如此大規(guī)模的天線陣列對數(shù)據(jù)采集、傳輸、處理等有極高的技術(shù)要求，特別是隨著陣列單元數(shù)的增加，各天線間信號交叉相關(guān)的實時計算是一個極其困難的任務(wù)。以低頻射電干涉陣列為例，它采用了IBM BlueGene/P超級計算機集群專門用于數(shù)據(jù)的在線處理。它通過各基站分別計算，合成一些波束后再進行干涉的辦法來降低所需的數(shù)據(jù)傳輸率和計算速度，即使這樣也無法實現(xiàn)所有單元同時進行干涉。如何應(yīng)對這些挑戰(zhàn)，尤其是如何以低廉的成本應(yīng)對目前萬億次每秒甚至億億次每秒的實時處理需求是國際上關(guān)注的一個難題。

圖形處理器(Graphics Processing Unit，GPU)的不斷發(fā)展為這類問題的解決帶來了新的希望。圖形處理器的概念自NVIDIA公司于1999年首次以圖像加速器提出以來，如今已經(jīng)發(fā)展成為一種高度并行化的多線程、眾核處理器，不僅具備高質(zhì)量和高性能的圖形處理能力，而且具有杰出的浮點計算能力和存儲帶寬，這使其能作為中央處理器的協(xié)處理器有效地加速科學(xué)運算。

國際上目前在積極利用圖形處理器為各大射電望遠鏡系統(tǒng)尋求一種成本更低、性能更好、功耗更少的數(shù)據(jù)處理解決方案。2008年10月，文[1]作者利用模擬數(shù)據(jù)測試了圖形處理器在做傅里葉變換和信號交叉相關(guān)時的性能，他們著重考察了不同的圖形處理器并行方案對計算性能的影響，雖然該實驗比較簡單，但圖形處理器的高性能、低功耗的潛能已經(jīng)充分表現(xiàn)出來。2009年8月，文[2]作者利用圖形處理器專門針對澳大利亞的默奇森寬視場射電望遠鏡陣列(Murchison Widefield Array，MWA)的原型系統(tǒng)(32個天線陣元)進行FX加速實驗，該實驗基于單個圖形處理器(型號為Nvidia C1060)實現(xiàn)了32路信號交叉相關(guān)的實時計算，結(jié)果表明圖形處理器運算性能高于單線程中央處理器程序68倍，但是由于PCI-E(Host與Device之間的傳輸協(xié)議)帶寬的限制，單圖形處理器無法勝任較大的天線陣元數(shù)目。同年，文[3]作者對多核中央處理器和眾核架構(gòu)下(比如IBM Cell，IBM BG/P，Intel Core i7，ATI Radeon 4870，Nvidia C1060)射電干涉儀相關(guān)器的性能做了全面的測試，結(jié)果IBM的兩個平臺取得了很高的性能。事實上，上述3個小組的實驗在很大程度上受到數(shù)據(jù)傳輸帶寬的束縛，雖然他們都使用了基本的內(nèi)存管理策略，比如共享內(nèi)存的窗口配置策略和寄存器優(yōu)化策略，但系統(tǒng)整體性能仍然只發(fā)揮了圖形處理器理論性能的10%到30%，這表明圖形處理器更多潛能需要進一步去挖掘。在這種情況下，2011年7月，文[4]作者利用新一代的Fermi架構(gòu)圖形處理器卡GTX480，綜合了前面3個小組的經(jīng)驗，采用多層次的內(nèi)存優(yōu)化策略，將信號交叉相關(guān)運算的整體性能提高至理論性能的79%，大大增加了圖形處理器在干涉儀相關(guān)器中的競爭優(yōu)勢。

1 射電干涉儀在運算和傳輸?shù)确矫娴男枨蠓治?/h2>

射電干涉儀在做信號的實時處理時，采取不同的方式，計算的復(fù)雜度有很大的差別，因此首先簡要介紹兩種基本的相關(guān)器。

1.1 相關(guān)器

射電干涉儀在做數(shù)據(jù)的實時處理時，可以采用XF方式，即先做實域的交叉互關(guān)聯(lián)(X-engine)，再做傅里葉變換(F-engine)；也可以使用相反的處理方式，即FX方式。同樣的信號數(shù)據(jù)，XF方式要比FX方式計算復(fù)雜度較大，因為先做傅里葉變換，再做交叉相關(guān)時，只需要做同頻率的相關(guān)，不同頻率間的相關(guān)為0。因此，隨著天線陣元的增加，國際上更傾向FX方式。

在實際情況中，可能還會采用FFX算法，即先做一個時間序列信號的快速傅里葉變換，得到時間頻域數(shù)據(jù)，再做一個空間序列信號的快速傅里葉變換，得到空間頻域數(shù)據(jù)，最后再做相關(guān)，這時只需計算時域頻率、空域頻率都相同的信號，這樣可以進一步降低運算量。上述空域快速傅里葉變換，可以采用二維的，即對分布在一個平面上的陣列，沿兩個方向均做快速傅里葉變換[5]。比如國家天文臺陳學(xué)雷研究員正研制的“天籟計劃”將采用柱面望遠鏡陣列，沿南北方向的陣列單元較長，而沿東西方向則只有有限的幾個，使用快速傅里葉變換并不能減少計算量，因此只需沿同一柱面做一維快速傅里葉變換，不同柱面間則直接計算相關(guān)。

1.2 FX相關(guān)器的數(shù)據(jù)傳輸和運算復(fù)雜度

本小節(jié)基于FX相關(guān)器，在給定干涉儀各種參數(shù)指標情況下，詳細討論數(shù)據(jù)傳輸壓力與運算復(fù)雜度。

假設(shè)干涉儀的天線數(shù)目為Na，每個天線的極化數(shù)為Np，采集的每個數(shù)據(jù)點通過模擬數(shù)字轉(zhuǎn)換為Nb位整型數(shù)據(jù)，干涉儀的頻率帶寬為Δν(根據(jù)Nyquist原理，干涉儀的采樣頻率Ns至少為2Δν)。相關(guān)器在做數(shù)據(jù)實時處理時，傅里葉變換的長度取Nf，相關(guān)運算需要考慮N2p個極化對和Na(Na＋1)/2個天線對(即基線的個數(shù)NB，其中包括了Na(Na－1)/2個交叉相關(guān)和Na個自相關(guān))，每一個相關(guān)運算的輸出結(jié)果是時間tI的積累。干涉儀的最大積分時間與它的基線長度Lb和地球的自轉(zhuǎn)速度ωr有關(guān)，基線的長度越長，允許的積分時間越短。因此，干涉儀的積分時間通常不能超過最長基線對應(yīng)的積分時間，即

這里ωr＝7.272×10－5rad/s；D為望遠鏡的口徑；f為過采樣的倍率。這里積分時間的限制是滿足干涉的基本要求，實際積分時間的設(shè)置還與科學(xué)目標有關(guān)，比如探尋快速射電暴需要毫秒級的積分時間。

1.2.1 相關(guān)器的數(shù)據(jù)傳輸率

圖1展示了FX相關(guān)器的數(shù)據(jù)傳輸示意圖，前端數(shù)據(jù)采集和數(shù)據(jù)的模擬數(shù)字轉(zhuǎn)換等部分已經(jīng)略去。該圖以模擬數(shù)字的輸出作為FX相關(guān)器的輸入，模擬數(shù)字的輸出為整型數(shù)據(jù)。

圖1 FX相關(guān)器的數(shù)據(jù)傳輸示意圖Fig.1 Illustration of the data transmission of the FX correlator

根據(jù)干涉儀的參數(shù)，可以估算每秒內(nèi)所有天線的數(shù)據(jù)產(chǎn)出，即模擬數(shù)字轉(zhuǎn)換的數(shù)據(jù)輸出為:

為了防止傅里葉變換或多向濾波器(Polyphase Filter，PPF)數(shù)據(jù)溢出，需要將輸入的整型數(shù)據(jù)擴展為高位的浮點型數(shù)據(jù)，這里因子2表示將數(shù)據(jù)的位數(shù)擴展了2倍。

在傅里葉變換階段，每一個天線采集的數(shù)據(jù)傅里葉變換接口的輸入或輸出數(shù)據(jù)率為(傅里葉變換不改變數(shù)據(jù)的存儲大小):

那么在傅里葉變換階段，總的數(shù)據(jù)輸入或輸出率為:

在交叉互關(guān)聯(lián)階段，即圖1所示的復(fù)數(shù)的乘累加(Complex-Multiplication-and-Accumulation，CMAC)階段，每個頻率通道向交叉互關(guān)聯(lián)的輸入數(shù)據(jù)率為:

交叉互關(guān)聯(lián)的輸出數(shù)據(jù)率較為復(fù)雜，因為它和積分時間有關(guān)，積分時間與基線長度、望遠鏡口徑、地球自轉(zhuǎn)速度等因素有關(guān)，積分時間的上限由(1)式給出。那么，交叉互關(guān)聯(lián)總的數(shù)據(jù)輸出率為:

因子2是指復(fù)數(shù)由實部和虛部兩部分實數(shù)組成，因子32意味著每個實數(shù)的字節(jié)數(shù)為32位，這是由望遠鏡設(shè)計者根據(jù)數(shù)據(jù)精度的要求確定。式中1/tI的意義為每秒數(shù)據(jù)輸出的次數(shù)，但實際數(shù)據(jù)的輸出頻率可能小于1/tI，比如積分時間為0.1 s，而數(shù)據(jù)輸出的頻率可能設(shè)置為每秒一次，這樣可以大大降低數(shù)據(jù)的輸出量。

1.2.2 相關(guān)器的運算復(fù)雜度

這里運算復(fù)雜度主要指相關(guān)器在執(zhí)行傅里葉變換和交叉互關(guān)聯(lián)操作時每秒的浮點運算次數(shù)。

Cooley-Tukey是一種常用的快速傅里葉變換算法，它的運算復(fù)雜度為(3Nf/2)log2Nf(其中包括(Nf/2)log2Nf次乘法運算和Nflog2Nf次加法運算)。由該復(fù)雜度公式可知，傅里葉變換的計算復(fù)雜度除了與天線的個數(shù)Na和極化數(shù)Np等因素有關(guān)之外，最主要是由傅里葉變換的長度Nf決定，即:

交叉互關(guān)聯(lián)是相關(guān)器計算量最大的一部分，也是相關(guān)器最核心的部分。對于每秒的輸入數(shù)據(jù)，相關(guān)運算都需要執(zhí)行NsNBN2p/Nf次復(fù)數(shù)的乘累加操作，每次復(fù)數(shù)的乘累加需要8次的浮點運算，所以每個頻率通道每秒的浮點運算次數(shù)為:

每秒浮點運算總量為:

1.2.3 “天籟計劃”的性能需求分析

“天籟計劃”是我國正在建設(shè)的一個大型射電干涉儀陣列，主要科學(xué)目標是探測宇宙中的暗能量[6]。

表1給出了“天籟計劃”一期和二期工程設(shè)備擬采用的參數(shù)設(shè)置。表2是根據(jù)表1的參數(shù)設(shè)置，估算了FX相關(guān)器在3個主要環(huán)節(jié)中數(shù)據(jù)傳輸和運算的需求。對于“天籟計劃”運算需求的簡單估算:一塊典型四核中央處理器(比如:Intel Core i7)的計算性能大概為70GFLOPS (0.54GFLOPS/Watt)，根據(jù)表2的估算結(jié)果，僅僅完成天籟一期工程交叉相關(guān)運算，至少需要430個Intel Core i7處理器完成交叉相關(guān)的計算；而一塊同時期的典型圖形處理器(比如Nvidia Geforce GTX 580)的計算性能至少為1.5TFLOPS (6.5GFLOPS/Watt)，完成天籟一期工程的實時數(shù)據(jù)處理一共只需30多個圖形處理器即可，無論從能耗還是從硬件成本上，圖形處理器方案都將遠遠低于中央處理器的方案。

表1 “天籟計劃”一期和二期工程擬采用的參數(shù)設(shè)置Table 1 The parameters planned for the Phase 1 and Phase 2 of the Tianlai project

表2 “天籟計劃”的數(shù)據(jù)傳輸和運算性能需求Table 2 The requirements on the data transmission and computing performances in the Tianlai project

2 FX相關(guān)器圖形處理器實現(xiàn)的優(yōu)化模型

如(7)式和(9)式所示，交叉互關(guān)聯(lián)的計算量將遠大于傅里葉變換的計算量，因此本文將重點討論交叉互關(guān)聯(lián)的性能優(yōu)化。

利用圖形處理器實現(xiàn)相關(guān)器，其性能主要由數(shù)據(jù)傳輸和運算兩部分決定，但數(shù)據(jù)傳輸通常是圖形處理器整體性能的瓶頸，受限于數(shù)據(jù)的傳輸速率，圖形處理器的實際性能一般只能達到理論性能的10%～30%。為了提升圖形處理器的實際性能，必須減少數(shù)據(jù)的反復(fù)讀寫，提高數(shù)據(jù)的利用率。圖形處理器的存儲器模型設(shè)計了多種存儲類型，包括全局存儲器、紋理存儲器、共享存儲器、寄存器等，根據(jù)不同存儲器的特點，合理使用多層次存儲器對數(shù)據(jù)進行緩存加速，可以有效減少圖形處理器對顯存的反復(fù)訪問，提升程序的性能。

FX相關(guān)器的設(shè)計，正是利用多層次存儲器，將線程模塊與數(shù)據(jù)模塊合理映射，實現(xiàn)了數(shù)據(jù)的有效調(diào)度，下面從基本的模型出發(fā)詳細討論線程模塊與數(shù)據(jù)模塊的配比。

在圖形處理器線程模型中[7]，多個線程構(gòu)成一個線程塊(block)，線程塊為線程的基本組織單位，多個線程塊組成一個線程網(wǎng)格(grid)。每個線程都擁有自己的寄存器，各線程寄存器之間不能直接通信。每個線程塊中所有的線程共享一個共享存儲器(shared memory)，同一個線程塊中的線程可以自由讀寫共享存儲器。雖然寄存器速度最快，但是寄存器為線程私有，線程之間不能通過寄存器共享數(shù)據(jù)，因而共享存儲器為圖形處理器上數(shù)據(jù)緩存方式的首選。

前文講過，F(xiàn)X相關(guān)器只需要計算同頻率信號之間的交叉關(guān)聯(lián)。各天線在某頻率f1上的信號構(gòu)成一個大小為Na的向量，那么該頻率上的交叉關(guān)聯(lián)為一個Na×Na的矩陣，其有效部分為矩陣的下(或上)三角(對應(yīng)基線個數(shù)NB)，如圖2中左邊的三角形所示。為了提高數(shù)據(jù)的利用率，需要盡可能多地將數(shù)據(jù)緩存在共享存儲器中，但是圖形處理器上的每個共享存儲器的大小有限，只能存儲部分數(shù)據(jù)，因而需要將數(shù)據(jù)分塊，分散到不同的線程塊分別進行處理。假設(shè)處理時將Na分解為寬度為W的B個正方形數(shù)據(jù)塊分別進行處理(Na＝WB)，那么對于一個頻率為f1的交叉關(guān)聯(lián)運算，圖形處理器的數(shù)據(jù)加載總次數(shù)為[4]:

圖2(左)所示，每個不同顏色的區(qū)域代表一個W×W大小的線程塊，其總數(shù)量為B(B＋1)/2個，每個線程塊每次處理在X方向和Y方向上都需要W個數(shù)據(jù)點。如果將這B(B＋1)/2個不同顏色的格子按照顏色從淺到深的方向展開，即可形成右圖中網(wǎng)格在gx方向的布局，如果再將gy方向上表示為不同的頻率，那么一個圖形處理器的線程網(wǎng)格就可以完成干涉儀陣列中所有基線上某序列頻率段(f1～fn)的交叉相關(guān)任務(wù)的計算。從(12)式可知，當W越大時，圖形處理器的數(shù)據(jù)加載次數(shù)越小，當W＝Na時，也即一個線程塊一次處理Na個數(shù)據(jù)，此時，數(shù)據(jù)加載次數(shù)Rin最小，即2Na。但是實際上每個線程塊的共享存儲器大小有限，且線程個數(shù)(處理能力)有限，在Na很大的情況下，W不可能等于Na，且當W逐漸增大時，線程塊中的可用資源(寄存器、指令發(fā)射單元、執(zhí)行單元等)的利用率會逐漸趨向于飽和，雖然Rin在減小，但是總的計算時間可能會成倍地增長，因而需要根據(jù)干涉儀的規(guī)模，同時結(jié)合圖形處理器硬件配置等情況合理地選擇W的大小，將Na分解到多個線程塊中去。

圖2 交叉關(guān)聯(lián)部分計算任務(wù)與線程數(shù)目的對應(yīng)關(guān)系[4]Fig.2 The mapping of data blocks in cross-correlation computations to grid threads[4]

(10)式中Rin的大小直接影響圖形處理器線程模型的布局，圖2中，圖形處理器線程模型中包括線程網(wǎng)格布局、線程塊布局和數(shù)據(jù)塊3部分，其中數(shù)據(jù)塊代表每個線程每次處理的數(shù)據(jù)大小。下面將對線程塊布局和數(shù)據(jù)塊的關(guān)系進行量化分析:

假設(shè)線程塊的寬度和高度分別為Bx和By，線程塊中每個線程每次讀取k個浮點數(shù)據(jù)，每個線程處理的數(shù)據(jù)塊的寬度和高度分別為Dx和Dy，線程塊所需的數(shù)據(jù)即為線程塊中所有線程讀取的數(shù)據(jù)之和:

式中，BxDx和ByDy分別為計算模型在X軸和Y軸方向的單位數(shù)據(jù)長度，對應(yīng)圖2中的W；數(shù)字4代表雙極浮點數(shù)據(jù)的大小(字節(jié))。由上面的分析可知每個線程能夠讀取的數(shù)據(jù)k的大小、數(shù)據(jù)塊的寬度Dx和高度Dy都需要根據(jù)圖形處理器的硬件參數(shù)進行適當調(diào)整。實際程序?qū)崿F(xiàn)時，大部分圖形處理器的配置并不一樣，而且由于摩爾定律的存在，圖形處理器的計算性能每年都會翻倍增長(圖形處理器硬件架構(gòu)或各硬件參數(shù)都將隨著改變)，因而這些參數(shù)的最優(yōu)取值并不會一成不變，只有根據(jù)實際硬件配置進行適當調(diào)整，才能使系統(tǒng)整體的性能達到最大。例如對先導(dǎo)實驗所使用的Nvidia GTX480顯卡，當Dx＝Dy＝2，Bx＝By＝8時，實際性能可以達到理論性能的77%。

3 關(guān)鍵技術(shù)的先導(dǎo)實驗

經(jīng)過充分的調(diào)研之后，分別在兩種不同圖形處理器架構(gòu)(AMD和NVIDIA)下，基于兩種圖形處理器編程模型(OpenCL和Cuda C)，利用復(fù)正弦模擬信號做了傅里葉和多路信號交叉相關(guān)的原型實驗。實驗過程中使用圖形處理器相關(guān)模型對信號進行互關(guān)聯(lián)，同時結(jié)合各種優(yōu)化措施，比如紋理緩存、架構(gòu)(cuda)流并行技術(shù)等，在GTX480顯卡上將交叉互關(guān)聯(lián)的計算性能從最初的100多GFLOPS(約理論性能的10%)提升至1033GFLOPS(約理論性能的77%)，如圖3。但從數(shù)據(jù)的傳輸速率角度看，主要瓶頸在于主機到設(shè)備的傳輸(即PCI-E的傳輸速率)，如圖4。

從圖3可以看出GTX460的Kernel性能(即只包括內(nèi)核程序執(zhí)行而不包括數(shù)據(jù)傳輸?shù)男阅?最高達到401Gflops，約占理論性能的44%，總體性能最高達到335GFlops，約占理論性能的37%。而GTX480的Kernel性能最高達到1033Gflops，約占理論性能的77%，總體性能最高達到786GFlops，約占理論性能的58%，具體如表3?？梢姛o論是Kernel性能還是總體性能，GTX480都比GTX460高出至少20個百分點，這里最主要的區(qū)別在于GTX480與GTX460硬件架構(gòu)不同:GTX480的核心架構(gòu)為GF100，每個核心中有32個流處理器；而GTX460的核心架構(gòu)為GF104，每個核心中有48個流處理器。但是他們的指令發(fā)射單元都是兩個，每個指令發(fā)射單元每次發(fā)射16個指令，因此GTX480中每次發(fā)射的指令數(shù)剛好與流處理器個數(shù)相符，而GTX460中有一個指令發(fā)射單元需要連續(xù)發(fā)射兩次指令才能滿足48個流處理器的處理需求，這就導(dǎo)致GTX460中的指令發(fā)射密度高于GTX480，即GTX460架構(gòu)需要更高的指令級并行，這種指令級并行最終要對應(yīng)到代碼的優(yōu)化上，本先導(dǎo)實驗的代碼設(shè)計主要是針對GTX480，因此GTX460的實驗性能相對理論性能較低。

圖3 GTX460與GTX480在不同天線個數(shù)時的計算性能曲線圖Fig.3 Computing speeds of GTX460 and GTX480 for different antenna numbers

圖4 GTX460與GTX480在不同天線個數(shù)時的數(shù)據(jù)傳輸速率曲線圖Fig.4 Data transmission rates of GTX460 and GTX480 for different antenna numbers

表3 GTX460與GTX480的Kernel性能和總體性能與理論性能的百分比Table 3 The performances of the kernel and overall performances of GTX460 and GTX480 (expressed in percentages of theoretically possible values)

這里計算的性能只包含算法本身的計算量，不包含系統(tǒng)實際運行所消耗的計算量。對于GTX480的77%的Kernel性能，未包括系統(tǒng)實際運行所消耗的運算量(如指令調(diào)度、內(nèi)存尋址等)，因此，Kernel的實際運算性能可能達到理論值的90%以上。所以，在今后性能優(yōu)化方面，Kernel算法優(yōu)化空間不大，將著重考慮數(shù)據(jù)傳輸方面的優(yōu)化。

頻譜泄漏相關(guān)的測試。用離散傅里葉變換作譜分析時，由于信號長度的截斷，總是會伴隨頻譜泄漏，這會對信號處理帶來諸多負面影響，對于這個問題也做了一定的研究和實驗。減小頻譜泄漏的方法，往往通過使用窗口函數(shù)，在先導(dǎo)實驗中，實驗了較為復(fù)雜的多相濾波器方法[8]，圖5是多相濾波器在不同環(huán)境下的實現(xiàn)結(jié)果，圖中每個子圖都有4條使用多相濾波器濾波后的曲線，多相濾波器有一個參數(shù)taps，即對信號分段濾波并疊加的次數(shù)，如圖當taps＝4時，圖中的黑線，頻譜泄漏已非常小。對于大數(shù)據(jù)量的信號，中央處理器架構(gòu)下運行性能往往要低于圖形處理器，但它可以檢查圖形處理器結(jié)果的精確度和正確性。如圖5所示，左圖為在中央處理器下用Matlab計算的結(jié)果，中間圖為在中央處理器下用C語言實現(xiàn)濾波用fftw函數(shù)庫做傅里葉變換的結(jié)果，右圖為在圖形處理器下用C語言實現(xiàn)濾波并用架構(gòu)的FFT庫做傅里葉變換的結(jié)果。

4 總結(jié)與展望

本文從目前國際上射電干涉儀陣列數(shù)據(jù)實時處理的運算瓶頸入手，提出了圖形處理器相關(guān)器的解決方案，然后詳細分析了FX相關(guān)器進行數(shù)據(jù)的實時處理時在各主要環(huán)節(jié)上數(shù)據(jù)的傳輸率和運算復(fù)雜度等問題，并根據(jù)“天籟計劃”的具體參數(shù)設(shè)置估算了該項目在數(shù)據(jù)傳輸和運算等方面的需求，詳細討論了使用圖形處理器進行信號交叉關(guān)聯(lián)時計算模型與線程模型的對應(yīng)關(guān)系。最后，在NVIDIA的單塊圖形處理器架構(gòu)下，基于Cuda C，利用復(fù)正弦模擬信號完成了傅里葉和多路信號交叉相關(guān)的原型實驗，并通過圖形處理器相關(guān)模型結(jié)合各種優(yōu)化技術(shù)將計算性能從最初的100多GFLOPS(約理論性能的10%)提升至1033GFLOPS(約理論性能的77%)，該結(jié)果證實了圖形處理器杰出的浮點計算能力和很高的傳輸帶寬，使圖形處理器在射電干涉儀相關(guān)器的開發(fā)上具有極大的潛力。更重要的是，單圖形處理器相關(guān)器的研發(fā)和實驗可以為今后圖形處理器集群環(huán)境下相關(guān)器的實現(xiàn)打下了良好的基礎(chǔ)。

圖5 頻譜泄漏測試中央處理器與圖形處理器結(jié)果驗證比較Fig.5 Test results of spectral leakage for a CPU architecture and a GPU architecture

單圖形處理器只能處理數(shù)個天線的信號處理，面對由幾百乃至上千個天線組成的大型射電干涉陣，需要圖形處理器集群才能應(yīng)對。先導(dǎo)實驗表明，圖形處理器的計算性能足以應(yīng)對復(fù)雜的運算，雖然圖形處理器各節(jié)點可以通過OpenMP或MPI等標準的消息傳遞協(xié)議完成任務(wù)的調(diào)度，但是上千路信號之間的交叉互相關(guān)，勢必導(dǎo)致圖形處理器節(jié)點之間的數(shù)據(jù)傳輸壓力過大。為解決這個問題，可以嘗試分頻分布式信號并行處理技術(shù):將同頻率的數(shù)據(jù)直接匯聚到相同的計算單元上，每個計算單元僅僅負責(zé)部分頻率點的相關(guān)性運算。這樣做的好處是，不需要計算單元之間的相互通信，減少了數(shù)據(jù)交換的壓力，并且方便計算節(jié)點的擴充，可以適用于將來更大規(guī)模的接收機。

在后面的研究工作中，將進一步嘗試多種優(yōu)化策略，比如使用多流水線的技巧，使計算和數(shù)據(jù)傳輸同時進行，以避免因傳輸?shù)钠款i而導(dǎo)致計算資源的等待，甚至還可以嘗試去掉從主機到設(shè)備的傳輸環(huán)節(jié)，可以將數(shù)據(jù)從望遠鏡的模擬數(shù)字轉(zhuǎn)換器直接送到圖形處理器的設(shè)備存儲器中。從GTX460和GTX480的整體對比情況來看，圖形處理器硬件的更新對FX相關(guān)器的影響也非常大，目前圖形處理器正處于快速發(fā)展時期，NVIDIA公司已經(jīng)推出新一代開普勒(Kepler)架構(gòu)圖形處理器，其性能比Fermi(GTX480即為Fermi架構(gòu))強3到4倍，而功耗更低，下一步將對該架構(gòu)圖形處理器進行測試。

總之，經(jīng)過不斷的測試，尋求各種最佳方案，利用圖形處理器逐步實現(xiàn)一套成本更低、性能更高、功耗更少的相關(guān)器，更好地滿足天籟計劃的需求。

參考文獻：

[1] Harris C，Haines K，Staveley-Smith L.GPU accelerated radio astronomy signal convolution[J]. Experimental Astronomy，2008，22(1-2):129－141.

[2] Wayth R B，Greenhill L J，Briggs F H.A GPU-based real-time software correlation system for the murchison wide-field array prototype[J].Publications of the Astronomical Society of the Pacific，2009，121:857－865.

[3] van Nieuwpoort R V，Romein J W.Using many-core hardware to correlate radio astronomy signals [C]//the 23rd ACM International Conference.2009:440－449.

[4] Clark M A，La Plante P C，Greenhill L J.Accelerating radio astronomy cross-correlation with graphics processing units[J].International Journal of High Performance Computing Applications，2013，27(2):178－192.

[5] Tegmark M，Zaldarriaga M.Omniscopes:Large area telescope arrays with only N logN computational cost[J].Physical Review D:Particles，F(xiàn)ields，Gravitation and Cosmology，2010，82:10pp.

[6] 陳學(xué)雷.暗能量的射電探測——天籟計劃簡介[J].中國科學(xué):物理學(xué)力學(xué)天文學(xué)，2011，41 (12):1358－1366.

Chen Xuelei.Radio detection of dark energy-the Tianlai project[J].Scientia Sinica:Physica，Mechanica＆Astronomica，2011，41(12):1358－1366.

[7] 張舒，褚艷利，趙開勇，等.GPU高性能運算之CUDA[M].北京:中國水利水電出版社，2009.

[8] Harris C J.A parallel model for the heterogeneous computation of radio astronomy signal correlation [M].Australia:The University of Western Australia，2009:33－38.

A Preliminary Study on GPU-based Correlators for a Radio Interferometer Array

Tian Haijun1，2，Xu Yang2，Chen Xuelei2，Li Changhua2，Wu Fengquan2，Wang Qunxiong1，Liu Yong1
(1.China Three Gorges University，Yichang 443002，China，Email:hjtian＠lamost.org；2.National Astronomical Observatories，Chinese Academy of Sciences，Beijing 100012，China)

As scales of radio interferometer telescope arrays increase，it becomes a tremendous challenge for the traditional solution to deal with intense real-time radio-astronomy data with high performances and low costs.In this paper we propose a creative and effective GPU-based solution to fit the needs of processing data from new arrays.We first give a brief description of the requirements on an effective solution，i.e.a highspeed low-cost procedure to transmit data and compute for real-time processing of data from a large-scale radio interferometer.We then present our GPU-based model，which improves practical computing performances of a single GPU card through accurately mapping data blocks in cross-correlation computations to GPU grid threads. The model achieves performances of up to 77%of those theoretically possible.Finally，several tests are carried by tailoring the model parameters to the requirements of the Tianlai project.The tests pave the way for realizing an on-line signal processing system with even much higher performances and a lower cost in an environment of GPU clusters.Such a system is to work for future more complicated and much larger radio interferometers.

Radio interferometer；FX Correlator；GPU；Real-time processing system

TP302

A

1672－7673(2014)03－0209－09

2013－10－21；

2013－11－14

田?？?，男，博士后.研究方向:天文信息學(xué).Email:hjtian＠lamost.org