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

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

1 引言

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

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

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

近年來(lái),隨著高性能計(jì)算的發(fā)展,圖形處理器(graphic processing unit,GPU)從專用于圖像領(lǐng)域的處理器逐漸向著通用并行計(jì)算平臺(tái)轉(zhuǎn)變。GPU已發(fā)展成為一種高度并行化、多線程、多核的通用計(jì)算設(shè)備,具有杰出的計(jì)算能力和極高的存儲(chǔ)器帶寬,并被越來(lái)越多地應(yīng)用于圖形處理之外的計(jì)算領(lǐng)域。GPU的運(yùn)算核心數(shù)遠(yuǎn)多于CPU,更適合于數(shù)據(jù)密集型計(jì)算的并行加速處理。NVIDIA于2007年推出了計(jì)算統(tǒng)一設(shè)備架構(gòu)(compute unified device architecture,CUDA),簡(jiǎn)化了GPU系統(tǒng)的開發(fā)流程,使得GPU通用計(jì)算技術(shù)在信號(hào)處理領(lǐng)域得到更為廣泛的應(yīng)用。目前,基于GPU的信號(hào)處理系統(tǒng)具有強(qiáng)大的并行處理能力,能夠在短時(shí)間內(nèi)通過(guò)并行處理完成大量數(shù)據(jù)的運(yùn)算,在GPU資源得到充分利用時(shí),可以實(shí)現(xiàn)對(duì)信號(hào)的實(shí)時(shí)處理要求。因此,基于GPU的信號(hào)處理技術(shù)成為眾多領(lǐng)域的熱點(diǎn),如射電天文[19-20]、雷達(dá)[21-31]、無(wú)線通信[32-33]、人工智能[34-35]等。

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

2 相位解卷繞算法研究

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

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

θ=φ+2k·π

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

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

2.2 相位解卷繞方法

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

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

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

(1)

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

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

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

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

3.1 算法并行性分析

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

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

對(duì)解卷繞過(guò)程分析可知,對(duì)存在卷繞現(xiàn)象的某一值解卷繞時(shí)只需確定其與真實(shí)相位的偏差。利用該偏差值對(duì)卷繞相位進(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個(gè)相位發(fā)生了卷繞,需要進(jìn)行解卷繞處理。每組的補(bǔ)償值相同,只需計(jì)算出第一個(gè)卷繞相位的補(bǔ)償值即可。

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

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

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

圖3為曲線型相位解卷繞示意圖,下面對(duì)曲線型卷繞的補(bǔ)償值進(jìn)行分析,從第2組開始,每組的第1個(gè)相位為卷繞相位。第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個(gè)值。由于每組有且只有第1個(gè)相位為卷繞相位且第1組未發(fā)生卷繞現(xiàn)象。則B可化簡(jiǎn)為:

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

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

(2)

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

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

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

(3)

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

Aunwrap=A+offset

(4)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

并行解卷繞算法的第一個(gè)模塊就是獲取卷繞相位信息,卷繞信息包括卷繞點(diǎn)的位置、卷繞類型和卷繞點(diǎn)數(shù)量。定義結(jié)構(gòu)體Wrap,包含id和type數(shù)據(jù)項(xiàng),分別對(duì)應(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為該序列的長(zhǎng)度,difference為相鄰兩個(gè)相位之間的差,wrap為輸出的卷繞信息結(jié)構(gòu)體序列,num為卷繞點(diǎn)的數(shù)量。調(diào)用size個(gè)線程,每個(gè)線程分別獲取各點(diǎn)卷繞信息,達(dá)到并行運(yùn)算的目的。每個(gè)線程內(nèi)的具體操作為:首先計(jì)算出該點(diǎn)與前1相位的差值difference,通過(guò)判斷該值是否在[-π,π]范圍內(nèi),來(lái)判斷該點(diǎn)是否為卷繞點(diǎn)。若是卷繞點(diǎn),則將該點(diǎn)的位置和卷繞類型寫入卷繞信息序列wrap。

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

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

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

獲取卷繞模塊已經(jīng)得到了卷繞信息序列、相鄰點(diǎn)的差值序列difference和卷繞點(diǎn)數(shù)量num,接下來(lái)需要建立補(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個(gè)線程,對(duì)應(yīng)為各個(gè)分組,每個(gè)線程獨(dú)立地計(jì)算該組的補(bǔ)償值。根據(jù)公式(2)計(jì)算出每組的補(bǔ)償值,并將每組的補(bǔ)償值拓展到該組內(nèi)的所有位置處,使得同一分組內(nèi)的補(bǔ)償值都相同,進(jìn)而得到補(bǔ)償序列Offset。

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

并行解卷繞算法的最后一個(gè)模塊就是并行補(bǔ)償模塊,利用補(bǔ)償序列Offset對(duì)初始卷繞相位序列的所有值進(jìn)行補(bǔ)償。調(diào)用size個(gè)線程,size為原始相位序列的長(zhǎng)度。每個(gè)線程獨(dú)立地對(duì)每個(gè)值進(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 仿真驗(yàn)證與結(jié)果分析

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10 加速效果對(duì)比2Fig.10 Comparison 2 of acceleration ratio

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

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

表5 卷繞相位數(shù)量固定情況下GPU并行算法較GPU串行算法加速效果對(duì)比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串行算法加速效果對(duì)比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串行解卷繞的加速效果對(duì)比圖。圖12為GPU并行解卷繞相比CPU串行解卷繞的加速效果對(duì)比圖。

圖11 加速效果對(duì)比3Fig.11 Comparison 3 of acceleration ratio

圖12 加速效果對(duì)比4Fig.12 Comparison 4 of acceleration ratio

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

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

6 結(jié)論

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

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