毛飛龍,焦義文,馬宏,*,張宇翔,聶欣林,高澤夫

1.航天工程大學(xué) 電子與光學(xué)工程系,北京 101416 2.國防科技大學(xué) 電子科學(xué)學(xué)院,長沙 710072 3.吉林大學(xué)白求恩第一醫(yī)院,長春 130000

1 引言

隨著航天技術(shù)與深空探測技術(shù)的迅猛發(fā)展,各國對空間的探測向著更深更遠(yuǎn)的范圍拓展。探測目標(biāo)向火星以及更遙遠(yuǎn)的太陽系大天體延伸,探測距離將從目前月球探測的4×105km拓展到火星探測的4×108km和木星探測的109km。距離越遠(yuǎn),信號到達(dá)地面的增益越低,這將對地面接收設(shè)備帶來巨大挑戰(zhàn)。通過天線組陣信號合成方法獲得更高的深空通信增益是值得大力發(fā)展的方法。天線組陣中的一項關(guān)鍵技術(shù)是信號間相位差估計,兩天線接收到同一信源的信號,互相關(guān)運(yùn)算得到天線間的互相關(guān)譜,之后經(jīng)過反正切相位鑒別器求得天線間的相位差。由于反正切運(yùn)算的存在,獲得的相位值被限制在[-π,π]之間。這些處于[-π,π]之間的相位就被稱作卷繞相位,為了得到真實(shí)的相位值信息,就必須將卷繞相位恢復(fù)為真實(shí)值,該過程就是相位解卷繞。

在行星探測任務(wù)中,探測器的跟蹤及精密測定軌在任務(wù)中占據(jù)重要地位,是完成工程任務(wù)和科學(xué)探測的基礎(chǔ)。中國嫦娥四號任務(wù)相時延處理軟件中,探測器卷繞相關(guān)相位的范圍是[-π,π],也需要進(jìn)行解卷繞處理。相位解卷繞對于探測器的定位和定軌,特別在探測器變軌、捕獲及下降著陸等關(guān)鍵弧段至關(guān)重要。因此,如何準(zhǔn)確且快速地進(jìn)行相位解卷繞對深空天線組陣信號合成、行星探測任務(wù)的變軌、捕獲等具有重要的工程意義。此外,相位解卷繞還是陣列信號相位差估計[1]、無線電干涉測量[2-5]、光學(xué)干涉儀[6-10]、核磁共振成像[11-14]等技術(shù)中的一個關(guān)鍵且基礎(chǔ)的算法。

在干擾較小時,只要從相位的第一個值開始逐點(diǎn)向后判斷并通過加減2π還原真值就可完成解卷繞。但信號處理、光學(xué)成像等系統(tǒng)對實(shí)時性的要求較高,使用上述的串行算法對大量數(shù)據(jù)進(jìn)行解卷繞將導(dǎo)致實(shí)時性急劇惡化。而目前對解卷繞的研究主要集中在提升其抗干擾性能[15-18],關(guān)于相位解卷繞并行實(shí)現(xiàn)的研究較少。因此,迫切地需要尋求一種準(zhǔn)確性好、實(shí)時性強(qiáng)的解卷繞方法。

近年來,隨著高性能計算的發(fā)展,圖形處理器(graphic processing unit,GPU)從專用于圖像領(lǐng)域的處理器逐漸向著通用并行計算平臺轉(zhuǎn)變。GPU已發(fā)展成為一種高度并行化、多線程、多核的通用計算設(shè)備,具有杰出的計算能力和極高的存儲器帶寬,并被越來越多地應(yīng)用于圖形處理之外的計算領(lǐng)域。GPU的運(yùn)算核心數(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)具有強(qiáng)大的并行處理能力,能夠在短時間內(nèi)通過并行處理完成大量數(shù)據(jù)的運(yùn)算,在GPU資源得到充分利用時,可以實(shí)現(xiàn)對信號的實(shí)時處理要求。因此,基于GPU的信號處理技術(shù)成為眾多領(lǐng)域的熱點(diǎn),如射電天文[19-20]、雷達(dá)[21-31]、無線通信[32-33]、人工智能[34-35]等。

本文以串行解卷繞算法為基礎(chǔ),通過對算法的并行映射尋求對實(shí)時性能的提升。本文在GPU平臺下設(shè)計了一種通用的并行相位解卷繞算法,并驗證了算法的正確性與實(shí)時性。

2 相位解卷繞算法研究

2.1 相位卷繞現(xiàn)象分析

相位卷繞現(xiàn)象一般出現(xiàn)在反正切相位鑒別之后,相位的主值被限制在了[-π,π]之間。這些被限制在[-π,π]之間的相位,與真實(shí)的相位相差2k·π(k為整數(shù))[4]。假設(shè)真實(shí)的相位為θ,卷繞的相位為φ,則有:

θ=φ+2k·π

本文將常見的相位卷繞現(xiàn)象分為兩類,分別為直線型相位卷繞和曲線型相位卷繞。直線型相位卷繞主要出現(xiàn)在干涉測量中,表現(xiàn)為直線的截斷。曲線型相位卷繞主要是余弦波形的相位卷繞,如旋轉(zhuǎn)相位干涉儀的卷繞現(xiàn)象,就是典型的曲線型相位卷繞[36],表現(xiàn)為對余弦波形的截斷。圖1為相位卷繞示意圖,圖1(a)為直線型相位卷繞原始信號波形,圖1(b)為其卷繞后的相位圖像。圖1(c)為曲線型相位卷繞原始信號波形,圖1(d)為其卷繞后的相位圖像。圖中兩條虛線表示[-π,π]范圍,可見卷繞后的相位值被限制在了該范圍內(nèi)。

圖1 相位卷繞現(xiàn)象示意Fig.1 Schematic of phase wrapping

2.2 相位解卷繞方法

目前,相位解卷繞的方法有很多,包括連續(xù)逐點(diǎn)解卷繞方法、拉普拉斯相位解卷繞方法等。其中最常用的方法就是通過逐點(diǎn)判斷并還原真實(shí)相位值的解卷繞方法。在卷繞相位中,由于發(fā)生卷繞點(diǎn)的相位會出現(xiàn)正負(fù)2π的跳變現(xiàn)象,而未發(fā)生卷繞的區(qū)域相位是近似連續(xù)的。因此,可以通過相鄰相位值之間的差來判斷是否出現(xiàn)相位跳變現(xiàn)象,也即卷繞現(xiàn)象。之后將卷繞的相位值通過加減2π來得到連續(xù)的、非卷繞的相位曲線。在一些文獻(xiàn)中,這種方法也被稱作區(qū)域生長法。

具體步驟為:如果相位卷繞圖中相鄰相位發(fā)生了大于+π的跳變量,則從跳變點(diǎn)開始后面所有相位值全都減去2π,如果相位跳變小于-π,則從跳變點(diǎn)開始后面所有相位全都加上2π,相鄰相位之間的差值處于-π和+π之間則不做處理,逐點(diǎn)處理完所有相位即完成相位解卷繞。

旋轉(zhuǎn)單基線可以利用數(shù)字積分器進(jìn)行相位的累加處理以達(dá)到解卷繞的作用。一種典型的數(shù)字積分器計算公式如(1)所示[37]。

(1)

式中:φi是當(dāng)前時刻的相位差,φi-1是上一時刻的相位差,φ(i)是積分器當(dāng)前累加的相位差,φ(i-1)是積分器上一次的相位差。該算法也是通過逐點(diǎn)判斷兩點(diǎn)間的相位差并還原來實(shí)現(xiàn)解卷繞。

在連續(xù)逐點(diǎn)解卷繞方法中,由于是對相位信息的逐點(diǎn)解卷繞,得到的解卷繞的結(jié)果非常準(zhǔn)確,但是由于該方法屬于串行過程,導(dǎo)致解卷繞的速度非常緩慢,限制了其在實(shí)時性要求較高的場景使用。拉普拉斯方法是由嚴(yán)格的數(shù)學(xué)公式推導(dǎo)而來,是對純數(shù)學(xué)公式的應(yīng)用,因此其解卷繞的速度明顯快于連續(xù)逐點(diǎn)解卷繞方法,不足之處就是拉普拉斯方法解卷繞準(zhǔn)確性不如連續(xù)逐點(diǎn)解卷繞方法,高頻區(qū)域的處理結(jié)果過于平滑。

本文旨在尋求一種準(zhǔn)確性好、實(shí)時性強(qiáng)的解卷繞方法。因此,本文嘗試在連續(xù)逐點(diǎn)串行解卷繞算法的基礎(chǔ)上,設(shè)計一種基于GPU的并行相位解卷繞方法。

3 并行相位解卷繞算法設(shè)計

3.1 算法并行性分析

在串行連續(xù)逐點(diǎn)解卷繞方法中,要確認(rèn)某一相位值是否卷繞,就要判斷該值與上一相位值的差是否在[-π,π]之間。主要有三種情況,若差值大于π,則該點(diǎn)有卷繞現(xiàn)象,此時卷繞相位的大小為真實(shí)相位值加2π所得的數(shù)值;若差值小于-π,則該點(diǎn)也有卷繞現(xiàn)象,此時卷繞相位的大小為真實(shí)相位值減2π所得的數(shù)值;若差值在[-π,π]范圍內(nèi),則該點(diǎn)不存在卷繞現(xiàn)象,不需要進(jìn)行解卷繞。串行連續(xù)逐點(diǎn)解卷繞方法對某一值解卷繞時需滿足其之前的所有值均為真實(shí)相位值(不存在卷繞現(xiàn)象),因此,該方法只能從頭到尾逐次進(jìn)行解卷繞,其耗時的主要原因也就在于此。

針對傳統(tǒng)方法串行的處理方法,若能將串行的處理過程轉(zhuǎn)化為對所有點(diǎn)的并行解卷繞,即可很好地解決實(shí)時性差的問題。

對解卷繞過程分析可知,對存在卷繞現(xiàn)象的某一值解卷繞時只需確定其與真實(shí)相位的偏差。利用該偏差值對卷繞相位進(jìn)行補(bǔ)償,即可完成該點(diǎn)的解卷繞。而在該卷繞相位之后至下一卷繞相位之間范圍內(nèi)(即不存在卷繞現(xiàn)象的范圍)各點(diǎn)的補(bǔ)償值與該卷繞相位的補(bǔ)償值相同。如圖2所示,將所有相位值按卷繞相位出現(xiàn)的位置分為6組(圖中六種不同顏色表示)。其中第1組不存在卷繞現(xiàn)象,不需要解卷繞,第2、3、4、5、6組的第1個相位發(fā)生了卷繞,需要進(jìn)行解卷繞處理。每組的補(bǔ)償值相同,只需計算出第一個卷繞相位的補(bǔ)償值即可。

圖2 直線型相位解卷繞示意Fig.2 Schematic of linear phase unwrapping

因此,要實(shí)現(xiàn)對各點(diǎn)的并行解卷繞,找出發(fā)生卷繞現(xiàn)象的位置并計算出各個卷繞相位的補(bǔ)償值是關(guān)鍵所在。

下面對直線型卷繞的補(bǔ)償值進(jìn)行分析,從第2組開始,每組的第1個相位為卷繞相位。第1組不需要解卷繞,補(bǔ)償值為0;第2組相位值向下跳變了2π,因此該組的補(bǔ)償值為2π×1;同理,第3組的補(bǔ)償值為2π×2;第4組的補(bǔ)償值為2π×3。

圖3為曲線型相位解卷繞示意圖,下面對曲線型卷繞的補(bǔ)償值進(jìn)行分析,從第2組開始,每組的第1個相位為卷繞相位。第1組不需要解卷繞,補(bǔ)償值為0;第2組相位值向下跳變了2π,因此該組的補(bǔ)償值為2π×1;同理,第3組的補(bǔ)償值為2π×2;第4組相位值向上跳變了2π,其補(bǔ)償值為上一組補(bǔ)償值的基礎(chǔ)上減2π,因此,第4組的補(bǔ)償值為2π×1;同理,第5組的補(bǔ)償值為0。

綜合兩種類型的卷繞現(xiàn)象可以得出,按卷繞點(diǎn)相位出現(xiàn)的位置進(jìn)行分組后,每組的補(bǔ)償值由之前所有卷繞點(diǎn)(包括當(dāng)前組的)的數(shù)量和卷繞類型決定。定義卷繞相位的卷繞類型:0為未發(fā)生卷繞,-1為向下跳變,1為向上跳變。設(shè)某一相位序列A為:

A={x1,x2,x3,…,xL}

定義序列A中各相位的卷繞類型構(gòu)成的序列為B:

B={a11,a12,a13,…,a21,a22,a23,…,a31…}

其中,aij代表第i組的第j個值。由于每組有且只有第1個相位為卷繞相位且第1組未發(fā)生卷繞現(xiàn)象。則B可化簡為:

B={a11,0,0,…,0,a21,0,0,…,0,a31…}

其中,a11=0,則第n組的解卷繞補(bǔ)償值valuen為:

(2)

當(dāng)卷繞相位為直線型時,所有分組的補(bǔ)償值相同,解卷繞補(bǔ)償值valuen可進(jìn)一步化簡,其中k為直線斜率:

當(dāng)卷繞相位為曲線型或任意卷繞類型時,解卷繞補(bǔ)償值通過式(2)計算得出。顯然,所有點(diǎn)的補(bǔ)償值序列offset為:

offset={value1,…,value2,…,valuen,…}

(3)

根據(jù)offset序列對所有點(diǎn)并行進(jìn)行補(bǔ)償(將現(xiàn)有相位與補(bǔ)償值相加),即可實(shí)現(xiàn)對所有相位的并行解卷繞。

Aunwrap=A+offset

(4)

3.2 算法實(shí)現(xiàn)流程

本文設(shè)計的基于GPU的相位解卷繞算法,相比于基于CPU的算法,主要改進(jìn)如下:利用GPU眾核優(yōu)勢,將解卷繞算法拆分為3個獨(dú)立的模塊,分別使每個模塊進(jìn)行多線程并發(fā)。傳統(tǒng)算法按數(shù)據(jù)的順序依次進(jìn)行獲取卷繞信息、建立補(bǔ)償值、進(jìn)行補(bǔ)償三個模塊。本文中的并行算法不再按照數(shù)據(jù)點(diǎn)的順序進(jìn)行運(yùn)算,而是根據(jù)3.1節(jié)中公式推導(dǎo)結(jié)果對所有數(shù)據(jù)點(diǎn)同時獲取卷繞信息,在得到所有點(diǎn)的卷繞信息之后,根據(jù)卷繞點(diǎn)的數(shù)量和位置將數(shù)據(jù)進(jìn)行分組,利用多個線程同時對不同的分組賦予不同的補(bǔ)償值。最后并發(fā)多個線程對所有點(diǎn)進(jìn)行相位補(bǔ)償。

圖4為并行相位解卷繞實(shí)現(xiàn)流程,具體可分為以下步驟:

圖4 并行相位解卷繞實(shí)現(xiàn)流程Fig.4 Implementation of parallel phase unwrapping

步驟1:設(shè)原始相位序列的長度為L。調(diào)用L個線程,判斷所有相鄰相位值之間的差是否在[-π,π]范圍內(nèi)從而獲取所有卷繞點(diǎn)的信息。包括卷繞點(diǎn)的位置、卷繞類型和卷繞點(diǎn)數(shù)量。

步驟2:按照卷繞點(diǎn)出現(xiàn)的位置對相位序列進(jìn)行分組,n個卷繞相位可分為n+1組。

步驟3:調(diào)用n+1個線程,根據(jù)公式(3)計算出各組的解卷繞相位補(bǔ)償值value,并根據(jù)公式(3)建立所有點(diǎn)的補(bǔ)償值序列offset。

然而筆者認(rèn)為，兩個罪名不僅調(diào)整范圍相差甚異，而且各自有行為的評價側(cè)重點(diǎn)。比較二者的客體不同點(diǎn)即可看出:

步驟4:調(diào)用L個線程,根據(jù)公式(4)對所有相位值進(jìn)行補(bǔ)償。

以上步驟將并行解卷繞過程拆分為三個模塊:獲取卷繞模塊、建立補(bǔ)償模塊、并行補(bǔ)償模塊,分別命名為get_wrap、get_offset、parallel_unwrap。

4 算法GPU實(shí)現(xiàn)

本節(jié)主要討論如何基于GPU的編程架構(gòu)對解卷繞算法進(jìn)行并行設(shè)計優(yōu)化,以進(jìn)一步提高運(yùn)算效率,滿足實(shí)際工程應(yīng)用中的實(shí)時信號處理需求。

4.1 并行優(yōu)化設(shè)計

CUDA是以大量線程來實(shí)現(xiàn)高吞吐量數(shù)據(jù)的實(shí)時并行處理,線程間越獨(dú)立,加速的效果越明顯。根據(jù)上文對算法并行性的分析,每個模塊內(nèi)每個線程的運(yùn)算都是獨(dú)立的。因此,本文設(shè)計的并行解卷繞算法適合用GPU進(jìn)行并行加速,可以極大地提高運(yùn)算效率,為后續(xù)信號處理的實(shí)時處理打下基礎(chǔ)。

圖5為GPU并行優(yōu)化示意圖,首先通過CPU分配好主機(jī)內(nèi)存和設(shè)備內(nèi)存,初始化各個變量。并利用cudaMemcpy將原始卷繞相位序列寫入GPU片上內(nèi)存。在實(shí)際處理過程中,由于GPU的線程數(shù)目有限,根據(jù)原始相位序列的長度選取GPU的數(shù)量。在每個GPU中以單指令多數(shù)據(jù)流(single instruction multiple data,SIMD)并行架構(gòu)進(jìn)行各模塊的運(yùn)算。所有模塊計算完成后,將解卷繞后的序列傳回CPU內(nèi)存或繼續(xù)在GPU中進(jìn)行后續(xù)的處理。

圖5 GPU并行優(yōu)化示意Fig.5 Schematic of GPU parallel optimization

4.2 獲取卷繞模塊GPU實(shí)現(xiàn)

并行解卷繞算法的第一個模塊就是獲取卷繞相位信息,卷繞信息包括卷繞點(diǎn)的位置、卷繞類型和卷繞點(diǎn)數(shù)量。定義結(jié)構(gòu)體Wrap,包含id和type數(shù)據(jù)項,分別對應(yīng)卷繞的位置和類型。如算法1。

算法1:獲取卷繞相位

Input:wrap_data、size

struct Wrap { int id;int type;};

__global__ void get_wrap(double *wrap_data,float *difference,Wrap *wrap,int size,int num)

{int tid = blockIdx.x * blockDim.x + threadIdx.x;

if(tid >= size - 1) return;

difference[tid + 1] = wrap_data[tid + 1] - wrap_data[tid];

if(difference[tid + 1] >CU_PI)

temp = atomicAdd(num,1);

wrap[temp].id = tid + 2; wrap[temp].type = 1;

if(difference[tid + 1] <-CU_PI)

temp = atomicAdd(num,1);

wrap[temp].id = tid + 2; wrap[temp].type = -1;

sort(wrap,wrap + num);

wrap_data為初始相位序列,size為該序列的長度,difference為相鄰兩個相位之間的差,wrap為輸出的卷繞信息結(jié)構(gòu)體序列,num為卷繞點(diǎn)的數(shù)量。調(diào)用size個線程,每個線程分別獲取各點(diǎn)卷繞信息,達(dá)到并行運(yùn)算的目的。每個線程內(nèi)的具體操作為:首先計算出該點(diǎn)與前1相位的差值difference,通過判斷該值是否在[-π,π]范圍內(nèi),來判斷該點(diǎn)是否為卷繞點(diǎn)。若是卷繞點(diǎn),則將該點(diǎn)的位置和卷繞類型寫入卷繞信息序列wrap。

由于線程并發(fā)時,存在多個線程同時訪問wrap內(nèi)存的問題,即訪問沖突導(dǎo)致寫入數(shù)據(jù)出錯的問題。采用原子操作方法進(jìn)行解決,原子操作是指不會被線程調(diào)度機(jī)制打斷的操作;這種操作一旦開始,就一直運(yùn)行到結(jié)束,中間不會有任何context switch(切換到另一個線程)。在多進(jìn)程(線程)訪問共享資源時,能夠確保所有其他的進(jìn)程(線程)都不在同一時間內(nèi)訪問相同的資源[38-42]。本文使用atomicAdd原子相加操作,確保各個線程不會同時訪問wrap,進(jìn)而得到所有卷繞點(diǎn)的信息。

同時,由于各個線程運(yùn)算結(jié)束的時間也各不相同,這就導(dǎo)致所有線程運(yùn)算完成后,得到一個無序的卷繞信息序列。需要進(jìn)一步對wrap序列按照其id數(shù)據(jù)項進(jìn)行排列。

4.3 建立補(bǔ)償模塊GPU實(shí)現(xiàn)

獲取卷繞模塊已經(jīng)得到了卷繞信息序列、相鄰點(diǎn)的差值序列difference和卷繞點(diǎn)數(shù)量num,接下來需要建立補(bǔ)償序列Offset如算法2。

算法2:建立補(bǔ)償序列Offset

Input:difference、wrap、num

Output:Offset

__global__ void get_offset(float *difference,Wrap *wrap,float *Offset,int num)

{int tid = blockIdx.x * blockDim.x + threadIdx.x;

if(tid >= num) return;

for(int i = 0;i <= tid - 1;i++)

{sum = wrap[i].type + sum;}

A =(-1)* sum *(2 * CU_PI);

seg_start = int(wrap[tid - 1].id);

seg_end = int(wrap[tid].id);

for(int i = seg_start - 1;i

Offset[i] = A;}

在本模塊,調(diào)用num個線程,對應(yīng)為各個分組,每個線程獨(dú)立地計算該組的補(bǔ)償值。根據(jù)公式(2)計算出每組的補(bǔ)償值,并將每組的補(bǔ)償值拓展到該組內(nèi)的所有位置處,使得同一分組內(nèi)的補(bǔ)償值都相同,進(jìn)而得到補(bǔ)償序列Offset。

4.4 并行補(bǔ)償模塊GPU實(shí)現(xiàn)

并行解卷繞算法的最后一個模塊就是并行補(bǔ)償模塊,利用補(bǔ)償序列Offset對初始卷繞相位序列的所有值進(jìn)行補(bǔ)償。調(diào)用size個線程,size為原始相位序列的長度。每個線程獨(dú)立地對每個值進(jìn)行補(bǔ)償,得到解卷繞之后的數(shù)據(jù)unwrap_data。

算法3:并行補(bǔ)償

Input:wrap_data、Offset

Output:unwrap_data、size

__global__ void parallel_unwrap(double *wrap_data,float *unwrap_data,float *Offset,int size)

{int tid = blockIdx.x * blockDim.x + threadIdx.x;

if(tid >= size) return;

unwrap_data[tid] = wrap_data[tid] + Offset[tid];}

5 仿真驗證與結(jié)果分析

在并行解卷繞仿真實(shí)驗中,選用表1和表2所示的GPU和CPU仿真平臺和具體仿真環(huán)境參數(shù)。

表1 GPU參數(shù)Table 1 Parameters of GPU

表2 CPU參數(shù)Table 2 Parameters of CPU

首先對并行解卷繞算法的正確性進(jìn)行驗證,仿真產(chǎn)生了1000點(diǎn)的卷繞相位序列,對序列進(jìn)行加噪處理,并分別使用CPU串行解卷繞方法和GPU并行解卷繞方法進(jìn)行運(yùn)算,解卷繞結(jié)果如圖6。

圖6 解卷繞結(jié)果Fig.6 Result of phase unwrapping

從圖中可以看出,在信噪比為30dB和60dB兩種情況下,本文設(shè)計的GPU并行解卷繞算法與傳統(tǒng)的CPU串行解卷繞算法結(jié)果一致,驗證了GPU并行解卷繞方法的正確性。

之后利用NVIDIA提供的分析工具NVIDIA Nsight Systems軟件統(tǒng)計了GPU并行解卷繞時間中CPU和GPU間的數(shù)據(jù)傳輸時間,以及GPU實(shí)際執(zhí)行各自模塊所占總執(zhí)行時間的比例,具體時序圖如圖7。

圖7 解卷繞算法運(yùn)行時序Fig.7 Timeline of phase unwrapping

本文設(shè)計的基于GPU的并行相位解卷繞算法通過三個模塊(獲取卷繞模塊、建立補(bǔ)償模塊、并行補(bǔ)償模塊)實(shí)現(xiàn)。統(tǒng)計了這些模塊的耗時占比,除此之外,還統(tǒng)計了CPU和GPU間的數(shù)據(jù)傳輸時間占比,即Memcpy HtoD和Memcpy DtoH。GPU并行相位解卷繞算法的各模塊耗時所占比例如圖8所示,其中,所占比例最大的是建立補(bǔ)償模塊,達(dá)到了49.58%,接下來是CPU到GPU數(shù)據(jù)傳輸和GPU到CPU數(shù)據(jù)傳輸,分別達(dá)到了25.39%、23.59%。占用時間最少的是獲取卷繞模塊和并行補(bǔ)償模塊,分別為0.05%、1.4%。

圖8 解卷繞算法各模塊耗時占比Fig.8 Time consumption ratio of each module of phase unwrapping

之后對并行解卷繞的加速比進(jìn)行仿真分析,考慮到在信號處理過程中,數(shù)據(jù)在GPU與CPU之間的傳輸耗時占比不可忽略,若信號處理系統(tǒng)通過GPU進(jìn)行運(yùn)算,則傳回至CPU進(jìn)行串行解卷繞需要耗費(fèi)大量的傳輸時間,而在GPU中直接調(diào)用單核進(jìn)行串行解卷繞是一種解決方法。因此本文仿真的對象分別為CPU串行解卷繞(CPU多線程)、GPU串行解卷繞(GPU單線程)、GPU并行解卷繞(GPU多線程)。

考慮到影響算法耗時的主要因素是原始相位數(shù)據(jù)量的大小和其中卷繞相位的數(shù)量。首先仿真產(chǎn)生6組原始相位數(shù)據(jù),數(shù)據(jù)量固定為1000000,卷繞相位數(shù)量為12～12000不等。對GPU并行算法和GPU串行算法分別進(jìn)行計時,得到表3所示中的耗時和加速比結(jié)果;對GPU并行算法和CPU串行算法分別進(jìn)行計時,得到表4所示中的耗時和加速比結(jié)果。

表3 數(shù)據(jù)量固定情況下GPU并行算法較GPU串行算法加速效果對比Table 3 Comparison of acceleration effects between GPU parallel algorithm and GPU serial algorithm under fixed data volume

表4 數(shù)據(jù)量固定情況下GPU并行算法較CPU串行算法加速效果對比Table 4 Comparison of acceleration effects between GPU parallel algorithm and CPU serial algorithm under fixed data volume

圖9為GPU并行解卷繞相比GPU串行解卷繞的加速效果對比圖。圖10為GPU并行解卷繞相比CPU串行解卷繞的加速效果對比圖。

圖10 加速效果對比2Fig.10 Comparison 2 of acceleration ratio

可以看到,在數(shù)據(jù)量固定時,隨著卷繞相位數(shù)量增加,GPU并行解卷繞算法的耗時先升高后降低。分析原因如下:建立補(bǔ)償序列模塊為本算法中耗時占比最大的模塊,在本模塊,調(diào)用線程的數(shù)量為卷繞相位的數(shù)量,每個線程對應(yīng)一個分組,每個線程獨(dú)立地計算該組的補(bǔ)償值。由此可知,本模塊內(nèi)線程并行的數(shù)量與卷繞相位的數(shù)量一致。因此,在數(shù)據(jù)量固定且卷繞相位小于6912時,卷繞相位數(shù)量越多,線程并行的數(shù)量也就越多,也就能更加充分地利用GPU的運(yùn)算資源。此時,本模塊總耗時與單個線程耗時的最大值相當(dāng),而線程并行的數(shù)量越多,單個線程處理的數(shù)據(jù)量就越少,耗時也就越短。而當(dāng)卷繞相位大于6912時,線程數(shù)超過了設(shè)備的最大值,線程調(diào)度出現(xiàn)了擁塞,進(jìn)而導(dǎo)致耗時增加。

下面考慮卷繞相位對算法耗時的影響,仿真產(chǎn)生3組原始相位數(shù)據(jù),卷繞相位數(shù)量固定為600,每組的數(shù)據(jù)量為100000-1000000不等。對GPU并行算法和GPU串行算法分別進(jìn)行計時,得到表5所示中的耗時和加速比結(jié)果;對GPU并行算法和CPU串行算法分別進(jìn)行計時,得到表6所示中的耗時和加速比結(jié)果。

表5 卷繞相位數(shù)量固定情況下GPU并行算法較GPU串行算法加速效果對比Table 5 Comparison of acceleration effects between GPU parallel algorithm and GPU serial algorithm when the numder of winding phases is fixed

表6 卷繞相位數(shù)量固定情況下GPU并行算法較CPU串行算法加速效果對比Table 6 Comparison of acceleration effects between GPU parallel algorithm and CPU serial algorithm when the number of winding phases is fixed

在卷繞相位數(shù)量固定的情況下,圖11為GPU并行解卷繞相比GPU串行解卷繞的加速效果對比圖。圖12為GPU并行解卷繞相比CPU串行解卷繞的加速效果對比圖。

圖11 加速效果對比3Fig.11 Comparison 3 of acceleration ratio

圖12 加速效果對比4Fig.12 Comparison 4 of acceleration ratio

隨著總數(shù)據(jù)量的增加,加速比呈現(xiàn)升高的趨勢。數(shù)據(jù)量越大,GPU并行的優(yōu)勢越明顯,這也印證了GPU確實(shí)適合大規(guī)模數(shù)據(jù)的實(shí)時處理。GPU并行解卷繞相比GPU串行解卷繞有63倍的加速比,GPU并行解卷繞相比CPU串行解卷繞約有3.5倍的加速比。出現(xiàn)兩個數(shù)值不同的原因主要包括并發(fā)線程的數(shù)量和處理器的主頻兩方面。本文使用的GPU與CPU規(guī)格參數(shù),GPU型號為NVIDIA TESLA A100,其時鐘頻率約為1.41GHz、最大核心數(shù)為6912。GPU并行算法與GPU串行算法都采用此GPU進(jìn)行運(yùn)算,二者唯一的區(qū)別就是線程并發(fā)的數(shù)量,GPU并行算法可同時并發(fā)成百上千個線程,而GPU串行算法則僅使用一個線程進(jìn)行運(yùn)算,因此加速比測試結(jié)果可達(dá)到文中的63倍。

GPU并行算法與CPU串行算法采用不同的處理器,CPU串行算法采用兩塊Intel(R)Xeon(R)Gold 6226R處理器,雖然其線程并發(fā)的最大數(shù)量為32個,但其主頻可達(dá)到2.9GHz。綜合線程并發(fā)數(shù)與主頻兩個因素,CPU串行算法的線程并發(fā)數(shù)與主頻都高于GPU串行算法,因此,GPU并行算法相比CPU串行算法無法達(dá)到63倍(GPU串行算法)的加速比,經(jīng)過多組實(shí)驗測試只能達(dá)到3.5倍加速比。

6 結(jié)論

針對大量數(shù)據(jù)串行相位解卷繞實(shí)時性較差的問題,分析了算法并行的可行性并設(shè)計了基于GPU的并行相位解卷繞算法。通過獲取卷繞模塊、建立補(bǔ)償模塊、并行補(bǔ)償模塊實(shí)現(xiàn)并行解卷繞。并利用線程并行、GPU并行、SIMD、原子操作等方法對算法進(jìn)行優(yōu)化。

本文中算法在不同信噪比條件下與傳統(tǒng)串行解卷繞方法結(jié)果一致。證明了本文中算法在GPU平臺下實(shí)現(xiàn)的正確性。之后對GPU并行解卷繞、GPU串行解卷繞、CPU串行解卷繞三種算法的實(shí)時處理能力進(jìn)行了對比。實(shí)驗表明,基于GPU的并行解卷繞算法相比CPU串行解卷繞算法有約3.5倍的加速比,相比GPU串行解卷繞算法有63倍的加速比。因此,本文設(shè)計的基于GPU的并行相位解卷繞算法對深空天線組陣信號合成、行星探測任務(wù)的變軌、捕獲等航天任務(wù)具有重要的工程意義。