毛飛龍,焦義文,馬宏,*,張宇翔,聶欣林,高澤夫
1.航天工程大學 電子與光學工程系,北京 101416 2.國防科技大學 電子科學學院,長沙 710072 3.吉林大學白求恩第一醫院,長春 130000
隨著航天技術與深空探測技術的迅猛發展,各國對空間的探測向著更深更遠的范圍拓展。探測目標向火星以及更遙遠的太陽系大天體延伸,探測距離將從目前月球探測的4×105km拓展到火星探測的4×108km和木星探測的109km。距離越遠,信號到達地面的增益越低,這將對地面接收設備帶來巨大挑戰。通過天線組陣信號合成方法獲得更高的深空通信增益是值得大力發展的方法。天線組陣中的一項關鍵技術是信號間相位差估計,兩天線接收到同一信源的信號,互相關運算得到天線間的互相關譜,之后經過反正切相位鑒別器求得天線間的相位差。由于反正切運算的存在,獲得的相位值被限制在[-π,π]之間。這些處于[-π,π]之間的相位就被稱作卷繞相位,為了得到真實的相位值信息,就必須將卷繞相位恢復為真實值,該過程就是相位解卷繞。
在行星探測任務中,探測器的跟蹤及精密測定軌在任務中占據重要地位,是完成工程任務和科學探測的基礎。中國嫦娥四號任務相時延處理軟件中,探測器卷繞相關相位的范圍是[-π,π],也需要進行解卷繞處理。相位解卷繞對于探測器的定位和定軌,特別在探測器變軌、捕獲及下降著陸等關鍵弧段至關重要。因此,如何準確且快速地進行相位解卷繞對深空天線組陣信號合成、行星探測任務的變軌、捕獲等具有重要的工程意義。此外,相位解卷繞還是陣列信號相位差估計[1]、無線電干涉測量[2-5]、光學干涉儀[6-10]、核磁共振成像[11-14]等技術中的一個關鍵且基礎的算法。
在干擾較小時,只要從相位的第一個值開始逐點向后判斷并通過加減2π還原真值就可完成解卷繞。但信號處理、光學成像等系統對實時性的要求較高,使用上述的串行算法對大量數據進行解卷繞將導致實時性急劇惡化。而目前對解卷繞的研究主要集中在提升其抗干擾性能[15-18],關于相位解卷繞并行實現的研究較少。因此,迫切地需要尋求一種準確性好、實時性強的解卷繞方法。
近年來,隨著高性能計算的發展,圖形處理器(graphic processing unit,GPU)從專用于圖像領域的處理器逐漸向著通用并行計算平臺轉變。GPU已發展成為一種高度并行化、多線程、多核的通用計算設備,具有杰出的計算能力和極高的存儲器帶寬,并被越來越多地應用于圖形處理之外的計算領域。GPU的運算核心數遠多于CPU,更適合于數據密集型計算的并行加速處理。NVIDIA于2007年推出了計算統一設備架構(compute unified device architecture,CUDA),簡化了GPU系統的開發流程,使得GPU通用計算技術在信號處理領域得到更為廣泛的應用。目前,基于GPU的信號處理系統具有強大的并行處理能力,能夠在短時間內通過并行處理完成大量數據的運算,在GPU資源得到充分利用時,可以實現對信號的實時處理要求。因此,基于GPU的信號處理技術成為眾多領域的熱點,如射電天文[19-20]、雷達[21-31]、無線通信[32-33]、人工智能[34-35]等。
本文以串行解卷繞算法為基礎,通過對算法的并行映射尋求對實時性能的提升。本文在GPU平臺下設計了一種通用的并行相位解卷繞算法,并驗證了算法的正確性與實時性。
相位卷繞現象一般出現在反正切相位鑒別之后,相位的主值被限制在了[-π,π]之間。這些被限制在[-π,π]之間的相位,與真實的相位相差2k·π(k為整數)[4]。假設真實的相位為θ,卷繞的相位為φ,則有:
θ=φ+2k·π
本文將常見的相位卷繞現象分為兩類,分別為直線型相位卷繞和曲線型相位卷繞。直線型相位卷繞主要出現在干涉測量中,表現為直線的截斷。曲線型相位卷繞主要是余弦波形的相位卷繞,如旋轉相位干涉儀的卷繞現象,就是典型的曲線型相位卷繞[36],表現為對余弦波形的截斷。圖1為相位卷繞示意圖,圖1(a)為直線型相位卷繞原始信號波形,圖1(b)為其卷繞后的相位圖像。圖1(c)為曲線型相位卷繞原始信號波形,圖1(d)為其卷繞后的相位圖像。圖中兩條虛線表示[-π,π]范圍,可見卷繞后的相位值被限制在了該范圍內。
圖1 相位卷繞現象示意Fig.1 Schematic of phase wrapping
目前,相位解卷繞的方法有很多,包括連續逐點解卷繞方法、拉普拉斯相位解卷繞方法等。其中最常用的方法就是通過逐點判斷并還原真實相位值的解卷繞方法。在卷繞相位中,由于發生卷繞點的相位會出現正負2π的跳變現象,而未發生卷繞的區域相位是近似連續的。因此,可以通過相鄰相位值之間的差來判斷是否出現相位跳變現象,也即卷繞現象。之后將卷繞的相位值通過加減2π來得到連續的、非卷繞的相位曲線。在一些文獻中,這種方法也被稱作區域生長法。
具體步驟為:如果相位卷繞圖中相鄰相位發生了大于+π的跳變量,則從跳變點開始后面所有相位值全都減去2π,如果相位跳變小于-π,則從跳變點開始后面所有相位全都加上2π,相鄰相位之間的差值處于-π和+π之間則不做處理,逐點處理完所有相位即完成相位解卷繞。
旋轉單基線可以利用數字積分器進行相位的累加處理以達到解卷繞的作用。一種典型的數字積分器計算公式如(1)所示[37]。
(1)
式中:φi是當前時刻的相位差,φi-1是上一時刻的相位差,φ(i)是積分器當前累加的相位差,φ(i-1)是積分器上一次的相位差。該算法也是通過逐點判斷兩點間的相位差并還原來實現解卷繞。
在連續逐點解卷繞方法中,由于是對相位信息的逐點解卷繞,得到的解卷繞的結果非常準確,但是由于該方法屬于串行過程,導致解卷繞的速度非常緩慢,限制了其在實時性要求較高的場景使用。拉普拉斯方法是由嚴格的數學公式推導而來,是對純數學公式的應用,因此其解卷繞的速度明顯快于連續逐點解卷繞方法,不足之處就是拉普拉斯方法解卷繞準確性不如連續逐點解卷繞方法,高頻區域的處理結果過于平滑。
本文旨在尋求一種準確性好、實時性強的解卷繞方法。因此,本文嘗試在連續逐點串行解卷繞算法的基礎上,設計一種基于GPU的并行相位解卷繞方法。
在串行連續逐點解卷繞方法中,要確認某一相位值是否卷繞,就要判斷該值與上一相位值的差是否在[-π,π]之間。主要有三種情況,若差值大于π,則該點有卷繞現象,此時卷繞相位的大小為真實相位值加2π所得的數值;若差值小于-π,則該點也有卷繞現象,此時卷繞相位的大小為真實相位值減2π所得的數值;若差值在[-π,π]范圍內,則該點不存在卷繞現象,不需要進行解卷繞。串行連續逐點解卷繞方法對某一值解卷繞時需滿足其之前的所有值均為真實相位值(不存在卷繞現象),因此,該方法只能從頭到尾逐次進行解卷繞,其耗時的主要原因也就在于此。
針對傳統方法串行的處理方法,若能將串行的處理過程轉化為對所有點的并行解卷繞,即可很好地解決實時性差的問題。
對解卷繞過程分析可知,對存在卷繞現象的某一值解卷繞時只需確定其與真實相位的偏差。利用該偏差值對卷繞相位進行補償,即可完成該點的解卷繞。而在該卷繞相位之后至下一卷繞相位之間范圍內(即不存在卷繞現象的范圍)各點的補償值與該卷繞相位的補償值相同。如圖2所示,將所有相位值按卷繞相位出現的位置分為6組(圖中六種不同顏色表示)。其中第1組不存在卷繞現象,不需要解卷繞,第2、3、4、5、6組的第1個相位發生了卷繞,需要進行解卷繞處理。每組的補償值相同,只需計算出第一個卷繞相位的補償值即可。
圖2 直線型相位解卷繞示意Fig.2 Schematic of linear phase unwrapping
因此,要實現對各點的并行解卷繞,找出發生卷繞現象的位置并計算出各個卷繞相位的補償值是關鍵所在。
下面對直線型卷繞的補償值進行分析,從第2組開始,每組的第1個相位為卷繞相位。第1組不需要解卷繞,補償值為0;第2組相位值向下跳變了2π,因此該組的補償值為2π×1;同理,第3組的補償值為2π×2;第4組的補償值為2π×3。
圖3為曲線型相位解卷繞示意圖,下面對曲線型卷繞的補償值進行分析,從第2組開始,每組的第1個相位為卷繞相位。第1組不需要解卷繞,補償值為0;第2組相位值向下跳變了2π,因此該組的補償值為2π×1;同理,第3組的補償值為2π×2;第4組相位值向上跳變了2π,其補償值為上一組補償值的基礎上減2π,因此,第4組的補償值為2π×1;同理,第5組的補償值為0。
綜合兩種類型的卷繞現象可以得出,按卷繞點相位出現的位置進行分組后,每組的補償值由之前所有卷繞點(包括當前組的)的數量和卷繞類型決定。定義卷繞相位的卷繞類型:0為未發生卷繞,-1為向下跳變,1為向上跳變。設某一相位序列A為:
A={x1,x2,x3,…,xL}
定義序列A中各相位的卷繞類型構成的序列為B:
B={a11,a12,a13,…,a21,a22,a23,…,a31…}
其中,aij代表第i組的第j個值。由于每組有且只有第1個相位為卷繞相位且第1組未發生卷繞現象。則B可化簡為:
B={a11,0,0,…,0,a21,0,0,…,0,a31…}
其中,a11=0,則第n組的解卷繞補償值valuen為:
(2)
當卷繞相位為直線型時,所有分組的補償值相同,解卷繞補償值valuen可進一步化簡,其中k為直線斜率:
當卷繞相位為曲線型或任意卷繞類型時,解卷繞補償值通過式(2)計算得出。顯然,所有點的補償值序列offset為:
offset={value1,…,value2,…,valuen,…}
(3)
根據offset序列對所有點并行進行補償(將現有相位與補償值相加),即可實現對所有相位的并行解卷繞。
Aunwrap=A+offset
(4)
本文設計的基于GPU的相位解卷繞算法,相比于基于CPU的算法,主要改進如下:利用GPU眾核優勢,將解卷繞算法拆分為3個獨立的模塊,分別使每個模塊進行多線程并發。傳統算法按數據的順序依次進行獲取卷繞信息、建立補償值、進行補償三個模塊。本文中的并行算法不再按照數據點的順序進行運算,而是根據3.1節中公式推導結果對所有數據點同時獲取卷繞信息,在得到所有點的卷繞信息之后,根據卷繞點的數量和位置將數據進行分組,利用多個線程同時對不同的分組賦予不同的補償值。最后并發多個線程對所有點進行相位補償。
圖4為并行相位解卷繞實現流程,具體可分為以下步驟:
圖4 并行相位解卷繞實現流程Fig.4 Implementation of parallel phase unwrapping
步驟1:設原始相位序列的長度為L。調用L個線程,判斷所有相鄰相位值之間的差是否在[-π,π]范圍內從而獲取所有卷繞點的信息。包括卷繞點的位置、卷繞類型和卷繞點數量。
步驟2:按照卷繞點出現的位置對相位序列進行分組,n個卷繞相位可分為n+1組。
步驟3:調用n+1個線程,根據公式(3)計算出各組的解卷繞相位補償值value,并根據公式(3)建立所有點的補償值序列offset。
然而筆者認為,兩個罪名不僅調整范圍相差甚異,而且各自有行為的評價側重點。比較二者的客體不同點即可看出:
步驟4:調用L個線程,根據公式(4)對所有相位值進行補償。
以上步驟將并行解卷繞過程拆分為三個模塊:獲取卷繞模塊、建立補償模塊、并行補償模塊,分別命名為get_wrap、get_offset、parallel_unwrap。
本節主要討論如何基于GPU的編程架構對解卷繞算法進行并行設計優化,以進一步提高運算效率,滿足實際工程應用中的實時信號處理需求。
CUDA是以大量線程來實現高吞吐量數據的實時并行處理,線程間越獨立,加速的效果越明顯。根據上文對算法并行性的分析,每個模塊內每個線程的運算都是獨立的。因此,本文設計的并行解卷繞算法適合用GPU進行并行加速,可以極大地提高運算效率,為后續信號處理的實時處理打下基礎。
圖5為GPU并行優化示意圖,首先通過CPU分配好主機內存和設備內存,初始化各個變量。并利用cudaMemcpy將原始卷繞相位序列寫入GPU片上內存。在實際處理過程中,由于GPU的線程數目有限,根據原始相位序列的長度選取GPU的數量。在每個GPU中以單指令多數據流(single instruction multiple data,SIMD)并行架構進行各模塊的運算。所有模塊計算完成后,將解卷繞后的序列傳回CPU內存或繼續在GPU中進行后續的處理。
圖5 GPU并行優化示意Fig.5 Schematic of GPU parallel optimization
并行解卷繞算法的第一個模塊就是獲取卷繞相位信息,卷繞信息包括卷繞點的位置、卷繞類型和卷繞點數量。定義結構體Wrap,包含id和type數據項,分別對應卷繞的位置和類型。如算法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為輸出的卷繞信息結構體序列,num為卷繞點的數量。調用size個線程,每個線程分別獲取各點卷繞信息,達到并行運算的目的。每個線程內的具體操作為:首先計算出該點與前1相位的差值difference,通過判斷該值是否在[-π,π]范圍內,來判斷該點是否為卷繞點。若是卷繞點,則將該點的位置和卷繞類型寫入卷繞信息序列wrap。
由于線程并發時,存在多個線程同時訪問wrap內存的問題,即訪問沖突導致寫入數據出錯的問題。采用原子操作方法進行解決,原子操作是指不會被線程調度機制打斷的操作;這種操作一旦開始,就一直運行到結束,中間不會有任何context switch(切換到另一個線程)。在多進程(線程)訪問共享資源時,能夠確保所有其他的進程(線程)都不在同一時間內訪問相同的資源[38-42]。本文使用atomicAdd原子相加操作,確保各個線程不會同時訪問wrap,進而得到所有卷繞點的信息。
同時,由于各個線程運算結束的時間也各不相同,這就導致所有線程運算完成后,得到一個無序的卷繞信息序列。需要進一步對wrap序列按照其id數據項進行排列。
獲取卷繞模塊已經得到了卷繞信息序列、相鄰點的差值序列difference和卷繞點數量num,接下來需要建立補償序列Offset如算法2。
算法2:建立補償序列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;} 在本模塊,調用num個線程,對應為各個分組,每個線程獨立地計算該組的補償值。根據公式(2)計算出每組的補償值,并將每組的補償值拓展到該組內的所有位置處,使得同一分組內的補償值都相同,進而得到補償序列Offset。 并行解卷繞算法的最后一個模塊就是并行補償模塊,利用補償序列Offset對初始卷繞相位序列的所有值進行補償。調用size個線程,size為原始相位序列的長度。每個線程獨立地對每個值進行補償,得到解卷繞之后的數據unwrap_data。 算法3:并行補償 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];} 在并行解卷繞仿真實驗中,選用表1和表2所示的GPU和CPU仿真平臺和具體仿真環境參數。 表1 GPU參數Table 1 Parameters of GPU 表2 CPU參數Table 2 Parameters of CPU 首先對并行解卷繞算法的正確性進行驗證,仿真產生了1000點的卷繞相位序列,對序列進行加噪處理,并分別使用CPU串行解卷繞方法和GPU并行解卷繞方法進行運算,解卷繞結果如圖6。 圖6 解卷繞結果Fig.6 Result of phase unwrapping 從圖中可以看出,在信噪比為30dB和60dB兩種情況下,本文設計的GPU并行解卷繞算法與傳統的CPU串行解卷繞算法結果一致,驗證了GPU并行解卷繞方法的正確性。 之后利用NVIDIA提供的分析工具NVIDIA Nsight Systems軟件統計了GPU并行解卷繞時間中CPU和GPU間的數據傳輸時間,以及GPU實際執行各自模塊所占總執行時間的比例,具體時序圖如圖7。 圖7 解卷繞算法運行時序Fig.7 Timeline of phase unwrapping 本文設計的基于GPU的并行相位解卷繞算法通過三個模塊(獲取卷繞模塊、建立補償模塊、并行補償模塊)實現。統計了這些模塊的耗時占比,除此之外,還統計了CPU和GPU間的數據傳輸時間占比,即Memcpy HtoD和Memcpy DtoH。GPU并行相位解卷繞算法的各模塊耗時所占比例如圖8所示,其中,所占比例最大的是建立補償模塊,達到了49.58%,接下來是CPU到GPU數據傳輸和GPU到CPU數據傳輸,分別達到了25.39%、23.59%。占用時間最少的是獲取卷繞模塊和并行補償模塊,分別為0.05%、1.4%。 圖8 解卷繞算法各模塊耗時占比Fig.8 Time consumption ratio of each module of phase unwrapping 之后對并行解卷繞的加速比進行仿真分析,考慮到在信號處理過程中,數據在GPU與CPU之間的傳輸耗時占比不可忽略,若信號處理系統通過GPU進行運算,則傳回至CPU進行串行解卷繞需要耗費大量的傳輸時間,而在GPU中直接調用單核進行串行解卷繞是一種解決方法。因此本文仿真的對象分別為CPU串行解卷繞(CPU多線程)、GPU串行解卷繞(GPU單線程)、GPU并行解卷繞(GPU多線程)。 考慮到影響算法耗時的主要因素是原始相位數據量的大小和其中卷繞相位的數量。首先仿真產生6組原始相位數據,數據量固定為1000000,卷繞相位數量為12~12000不等。對GPU并行算法和GPU串行算法分別進行計時,得到表3所示中的耗時和加速比結果;對GPU并行算法和CPU串行算法分別進行計時,得到表4所示中的耗時和加速比結果。 表3 數據量固定情況下GPU并行算法較GPU串行算法加速效果對比Table 3 Comparison of acceleration effects between GPU parallel algorithm and GPU serial algorithm under fixed data volume 表4 數據量固定情況下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 可以看到,在數據量固定時,隨著卷繞相位數量增加,GPU并行解卷繞算法的耗時先升高后降低。分析原因如下:建立補償序列模塊為本算法中耗時占比最大的模塊,在本模塊,調用線程的數量為卷繞相位的數量,每個線程對應一個分組,每個線程獨立地計算該組的補償值。由此可知,本模塊內線程并行的數量與卷繞相位的數量一致。因此,在數據量固定且卷繞相位小于6912時,卷繞相位數量越多,線程并行的數量也就越多,也就能更加充分地利用GPU的運算資源。此時,本模塊總耗時與單個線程耗時的最大值相當,而線程并行的數量越多,單個線程處理的數據量就越少,耗時也就越短。而當卷繞相位大于6912時,線程數超過了設備的最大值,線程調度出現了擁塞,進而導致耗時增加。 下面考慮卷繞相位對算法耗時的影響,仿真產生3組原始相位數據,卷繞相位數量固定為600,每組的數據量為100000-1000000不等。對GPU并行算法和GPU串行算法分別進行計時,得到表5所示中的耗時和加速比結果;對GPU并行算法和CPU串行算法分別進行計時,得到表6所示中的耗時和加速比結果。 表5 卷繞相位數量固定情況下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 卷繞相位數量固定情況下GPU并行算法較CPU串行算法加速效果對比Table 6 Comparison of acceleration effects between GPU parallel algorithm and CPU serial algorithm when the number of winding phases is fixed 在卷繞相位數量固定的情況下,圖11為GPU并行解卷繞相比GPU串行解卷繞的加速效果對比圖。圖12為GPU并行解卷繞相比CPU串行解卷繞的加速效果對比圖。 圖11 加速效果對比3Fig.11 Comparison 3 of acceleration ratio 圖12 加速效果對比4Fig.12 Comparison 4 of acceleration ratio 隨著總數據量的增加,加速比呈現升高的趨勢。數據量越大,GPU并行的優勢越明顯,這也印證了GPU確實適合大規模數據的實時處理。GPU并行解卷繞相比GPU串行解卷繞有63倍的加速比,GPU并行解卷繞相比CPU串行解卷繞約有3.5倍的加速比。出現兩個數值不同的原因主要包括并發線程的數量和處理器的主頻兩方面。本文使用的GPU與CPU規格參數,GPU型號為NVIDIA TESLA A100,其時鐘頻率約為1.41GHz、最大核心數為6912。GPU并行算法與GPU串行算法都采用此GPU進行運算,二者唯一的區別就是線程并發的數量,GPU并行算法可同時并發成百上千個線程,而GPU串行算法則僅使用一個線程進行運算,因此加速比測試結果可達到文中的63倍。 GPU并行算法與CPU串行算法采用不同的處理器,CPU串行算法采用兩塊Intel(R)Xeon(R)Gold 6226R處理器,雖然其線程并發的最大數量為32個,但其主頻可達到2.9GHz。綜合線程并發數與主頻兩個因素,CPU串行算法的線程并發數與主頻都高于GPU串行算法,因此,GPU并行算法相比CPU串行算法無法達到63倍(GPU串行算法)的加速比,經過多組實驗測試只能達到3.5倍加速比。 針對大量數據串行相位解卷繞實時性較差的問題,分析了算法并行的可行性并設計了基于GPU的并行相位解卷繞算法。通過獲取卷繞模塊、建立補償模塊、并行補償模塊實現并行解卷繞。并利用線程并行、GPU并行、SIMD、原子操作等方法對算法進行優化。 本文中算法在不同信噪比條件下與傳統串行解卷繞方法結果一致。證明了本文中算法在GPU平臺下實現的正確性。之后對GPU并行解卷繞、GPU串行解卷繞、CPU串行解卷繞三種算法的實時處理能力進行了對比。實驗表明,基于GPU的并行解卷繞算法相比CPU串行解卷繞算法有約3.5倍的加速比,相比GPU串行解卷繞算法有63倍的加速比。因此,本文設計的基于GPU的并行相位解卷繞算法對深空天線組陣信號合成、行星探測任務的變軌、捕獲等航天任務具有重要的工程意義。4.4 并行補償模塊GPU實現
5 仿真驗證與結果分析
6 結論