朱新峰，陸小霞，黎仁剛，王國(guó)海，王彥凌

(1.揚(yáng)州大學(xué)，江蘇 揚(yáng)州 225001；2.中國(guó)船舶重工集團(tuán)公司第七二三研究所，江蘇 揚(yáng)州 225101；3.中電科航空電子有限公司，四川 成都 611731;4.浙江理工大學(xué)，浙江 杭州 310018)

0 引 言

空間譜搜索[1]是陣列信號(hào)處理中的一個(gè)重要研究方向，在雷達(dá)、通信、聲吶等眾多領(lǐng)域有著極為廣闊的應(yīng)用前景。目前對(duì)于空間譜搜索主要用信號(hào)子空間類方法[2]來(lái)處理。信號(hào)子空間類方法主要由2步組成:一是對(duì)接收數(shù)據(jù)協(xié)方差矩陣進(jìn)行正交分解來(lái)估計(jì)信號(hào)(或噪聲)子空間;二是進(jìn)行譜峰搜索來(lái)提取波達(dá)方向估值。無(wú)論是正交分解還是譜值計(jì)算都面臨著巨大的運(yùn)算量。在信息技術(shù)的不斷更新中，傳統(tǒng)測(cè)向技術(shù)[3]已經(jīng)無(wú)法滿足當(dāng)下的需求，因此以MUSIC算法為代表的超分辨測(cè)向算法應(yīng)運(yùn)而生[4]。然而MUSIC算法中接收矩陣的子空間分解部分計(jì)算量龐大，影響了計(jì)算效率。為了高效實(shí)現(xiàn)MUSIC算法，因此考慮將MUSIC算法放在GPU下去實(shí)現(xiàn)，GPU可以提供數(shù)十倍乃至于上百倍于CPU的性能，能夠明顯提升運(yùn)行的速度。

1 MUSIC原理和運(yùn)算流程

1.1 MUSIC原理

在經(jīng)典窄帶MUSIC算法中，通常要求入射信號(hào)源數(shù)目小于陣元數(shù)，即P

設(shè)快拍數(shù)為N，根據(jù)式(1)所示接收數(shù)據(jù)模型，可求得陣列接收數(shù)據(jù)的協(xié)方差矩陣Rxx為：

Rxx= E[x(t)xH(t)]=

AE[s(t)sH(t)]AH+E[n(t)nH(t)]=

ARssAH+σ2I

(1)

式中：E[*]表示求算數(shù)平均；Rss表示原始信號(hào)的協(xié)方差矩陣，當(dāng)信號(hào)間相互獨(dú)立時(shí)，Rss為P維對(duì)角矩陣，對(duì)角線上值為對(duì)應(yīng)位置信號(hào)的平均功率；σ2為噪聲功率；I為M維單位矩陣。

因?yàn)橄嗷オ?dú)立的假設(shè)前提，不同陣元噪聲間、信號(hào)與噪聲間的協(xié)方差值接近0，所以接收數(shù)據(jù)的協(xié)方差矩陣可近似為信號(hào)矩陣的主對(duì)角線上加上噪聲功率。對(duì)Rxx進(jìn)行奇異值分解，得到特征向量矩陣U、V和對(duì)應(yīng)特征值組成的對(duì)角矩陣Σ：

Rxx=UΣV

(2)

式中：Σ對(duì)角線呈遞減排列時(shí)，可根據(jù)值大小將特征向量矩陣分為信號(hào)子空間Us與噪聲子空間Un。

(3)

U=[UsUn]

(4)

奇異值分解所得的特征向量?jī)蓛烧?，所得信?hào)子空間與噪聲子空間必然正交，且信號(hào)子空間與入射信號(hào)陣列流形張成的空間一致，所以陣列流形中任意導(dǎo)向矢量與噪聲子空間正交。

aiUn=0，i=1,2,…,P

(5)

但在實(shí)際情況中，由于噪聲的存在，以及入射信號(hào)復(fù)包絡(luò)近似不變的假設(shè)，導(dǎo)致估計(jì)結(jié)果存在一定偏差，導(dǎo)向矢量與噪聲子空間不可能完全正交。

因此，通過(guò)掃描不同入射角構(gòu)造各來(lái)波方向上的陣列流形，求取式(3)模值平方，取得極小值的點(diǎn)即為來(lái)波方向，如式(6):

(6)

式中：argmin{*}表示求*最小值時(shí)掃描變量的值。

為了更加直觀地表示來(lái)波方向，且使接近的2個(gè)峰值更易于分辨，繪制空間譜圖時(shí)通常取倒數(shù)，如式(7):

(7)

式中：Pmusic表示空間譜。

根據(jù)MUSIC算法的基本原理，繪出MUSIC算法的一般流程，如圖1所示。

圖1 MUSIC算法流程圖

1.2 算法運(yùn)行時(shí)間資源分配

以五元非均勻陣列為目標(biāo)陣形，入射方向分別為10°，30°，60°，信噪比為10 dB，高斯白噪聲為500的線性調(diào)頻信號(hào)作為入射信號(hào)，MUSIC算法計(jì)算出的空間譜如圖2所示。

圖2 5元陣空間譜

圖2中，3個(gè)最大譜峰處的方向角坐標(biāo)值即為來(lái)波方向，算法能夠準(zhǔn)確估算出3個(gè)入射信號(hào)的方向，具有優(yōu)越的多信號(hào)同時(shí)測(cè)向能力。

在時(shí)間上，算法中陣元數(shù)為9，使用信噪比是60，中心頻率是18e9，帶寬是20e6，采樣頻率是40e6，根據(jù)上面數(shù)據(jù)計(jì)算出采樣時(shí)間、快拍持續(xù)時(shí)間、頻率變化率，并在Matlab平臺(tái)下進(jìn)行仿真，如圖3所示。

圖3 九元陣方向角空間譜

使用的總時(shí)間為2.148 926 s，并且單獨(dú)測(cè)試了計(jì)算占用時(shí)間、接收矩陣、協(xié)方差計(jì)算、奇異值分解(SVD)，譜峰搜索分別穩(wěn)定在0.001 462 s，0.002 166 s，0.003 893 s，1.963 615 s，各占大約0.07%,0.14%,0.21%,0.58%?？梢钥闯銎渲蠸VD和譜峰搜索占了大部分的計(jì)算時(shí)間。本文著重進(jìn)行SVD和譜峰搜索部分的優(yōu)化，以提高計(jì)算效率。

2 使用GPU提升MUSIC算法運(yùn)算效率

2.1 MUSIC算法分解

2.1.1 接收數(shù)據(jù)

寬帶信號(hào)與窄帶信號(hào)[5]是一個(gè)相對(duì)概念，不同系統(tǒng)中的寬帶信號(hào)與窄帶信號(hào)具有不同的區(qū)分標(biāo)準(zhǔn)，且寬帶信號(hào)與窄帶信號(hào)的處理方法也有較大差距。

窄帶信號(hào)是指信號(hào)帶寬遠(yuǎn)遠(yuǎn)小于信號(hào)波前跨越陣列最大口徑所需要的時(shí)間的倒數(shù)，即：

(8)

式中：L為陣列最大口徑；F和λ為信號(hào)中心頻率和該頻率對(duì)應(yīng)的波長(zhǎng)。

2.1.2 奇異值分解

基于單邊Jacobi[6]旋轉(zhuǎn)的SVD算法,相對(duì)于雙邊、單邊的計(jì)算量小，并且容易并行實(shí)現(xiàn)。單邊Jacobi方法直接對(duì)原矩陣A進(jìn)行單邊正交旋轉(zhuǎn)，A可以是任意矩陣。AJ1J2…JK=B，同樣每一輪的迭代都要使A的任意兩列都正交一次，迭代退出的條件是B的任意兩列都正交。單邊Jacobi旋轉(zhuǎn)有這么一個(gè)性質(zhì)：旋轉(zhuǎn)前若‖ap‖>‖aq‖，則旋轉(zhuǎn)后依然是‖ap‖>‖aq‖；反之亦然。把矩陣B每列的模長(zhǎng)提取出來(lái)，B=U·Σ，把J1J2J3…Jk記為V，則A=V·Σ·VT。

2.1.3 譜峰搜索

應(yīng)用于MUSIC算法的譜峰搜索方法多種多樣，主要有柵格化多細(xì)粒度搜索，應(yīng)用于特殊陣列的免搜索或快速搜索算法，以及最優(yōu)化求解的智能優(yōu)化算法。根據(jù)系統(tǒng)資源與實(shí)現(xiàn)需求的不同，譜峰搜索的選擇也不盡相同[5]。

由于本文系統(tǒng)的實(shí)時(shí)性要求較高，選取了三級(jí)細(xì)粒度直接搜索的方法。分別采用5°，1°，0.25°對(duì)全方位，一級(jí)點(diǎn)[-10°,10°]，二級(jí)點(diǎn)[-2°,2°]作為柵格大小，如圖4所示。

圖4 常用譜峰搜索方法

2.2 算法的GPU高效實(shí)現(xiàn)

為了高效實(shí)現(xiàn)MUSIC算法，將MUSIC算法中部分復(fù)雜運(yùn)算放在GPU下去實(shí)現(xiàn)，GPU可以提供數(shù)十倍乃至于上百倍CPU的性能，將接收信號(hào)的模型、協(xié)方差的計(jì)算、SVD放入GPU里去并行計(jì)算，能夠明顯提升運(yùn)行的速度。

3 算法仿真

3.1 改進(jìn)SVD

為了提高復(fù)數(shù)SVD的計(jì)算效率，將獲取列范數(shù)，歸并排序，矩陣排序，單邊Jacobi以及奇異值的排序放在GPU里去進(jìn)行并行計(jì)算。CPU與GPU上的運(yùn)行時(shí)間對(duì)比如圖5所示。

圖5 CPU與GPU上的運(yùn)行時(shí)間對(duì)比

由圖5可以看出，在GPU里的運(yùn)行時(shí)間最終穩(wěn)定在0.00 255 s，與原先在CPU里的運(yùn)行時(shí)間0.003 893 s相比速度明顯加快。

下面在譜峰搜索部分GPU的實(shí)現(xiàn)中進(jìn)行了使用共享內(nèi)存、分配鎖頁(yè)內(nèi)存和增加CUDA流操作3個(gè)方面的優(yōu)化。

3.2 使用共享內(nèi)存

首先將信噪比、中心頻率、帶寬、采樣頻率、采樣時(shí)間、快拍持續(xù)時(shí)間、頻率變化率設(shè)置為全局變量，減少大量的內(nèi)存獲取請(qǐng)求，從而減少時(shí)間的交互。

然后再將譜峰搜索中的計(jì)數(shù)器設(shè)置為共享內(nèi)存，在調(diào)用過(guò)程中便不需要進(jìn)行繁瑣的交換數(shù)據(jù)。并且將調(diào)用__syncthread()函數(shù)來(lái)完成線程同步，這個(gè)函數(shù)的調(diào)用將確保線程塊中的每個(gè)線程都執(zhí)行完__syncthreads()前面的語(yǔ)句，才會(huì)執(zhí)行下一條語(yǔ)句。這樣可以確保對(duì)共享內(nèi)存進(jìn)行讀取之前，想要寫入的操作已經(jīng)完成。優(yōu)化前后的時(shí)間對(duì)比如圖6所示。

圖6 GPU上使用共享內(nèi)存前后的時(shí)間對(duì)比

從圖6可以看出，沒(méi)有使用共享內(nèi)存時(shí)，GPU的時(shí)間穩(wěn)定在350 ms左右，使用共享內(nèi)存后時(shí)間穩(wěn)定在315 ms左右，有了執(zhí)行時(shí)間降低了35 ms左右，性能提升10%左右。

3.3 分配鎖頁(yè)內(nèi)存

PCI-E數(shù)據(jù)總線傳輸?shù)臄?shù)據(jù)始終是鎖頁(yè)內(nèi)存中的數(shù)據(jù)，若應(yīng)用程序沒(méi)有對(duì)鎖頁(yè)內(nèi)存進(jìn)行分配，則CUDA驅(qū)動(dòng)會(huì)在后臺(tái)為程序分配。這之間必然存在分頁(yè)內(nèi)存與鎖頁(yè)內(nèi)存中不必要的拷貝操作，所以，通過(guò)分配鎖頁(yè)內(nèi)存來(lái)減少這些拷貝操作，提高效率[8]。

在應(yīng)用程序中，將主機(jī)端分配內(nèi)存的代碼塊:

double*host_n=(double *)malloc(size_of_n_in_bytes);

替換為：

double *host;

CUDA_CALL(cudaMallocHost(&host_n,size_of_n_in_bytes));

并將最后釋放內(nèi)存的代碼：

Free(host_n);

替換為：

CUDA_CALL(host_n)。

測(cè)試運(yùn)行時(shí)間，仿真結(jié)果如圖7所示。

圖7 GPU上分頁(yè)內(nèi)存與鎖存內(nèi)存的時(shí)間對(duì)比

從圖7可以看出，使用分頁(yè)內(nèi)存時(shí)運(yùn)行時(shí)間穩(wěn)定在315 ms，使用鎖頁(yè)內(nèi)存時(shí)穩(wěn)定在310 ms左右，時(shí)間減少5 ms，性能上提升了1.6%，提升的空間并沒(méi)有前面使用的共享內(nèi)存優(yōu)化那么明顯。

3.4 增加CUDA流操作

CUDA流[8]表示一個(gè)GPU操作隊(duì)列，并且該隊(duì)列中的操作將以指定的順序執(zhí)行。我們可以在流中添加一些操作，如核函數(shù)啟動(dòng)、內(nèi)存復(fù)制等。將這些操作添加到流的順序也就是他們的執(zhí)行順序。將每個(gè)流視為GPU上的一個(gè)任務(wù)，并且這些任務(wù)可以并行執(zhí)行。

首先是創(chuàng)建一個(gè)流：

cudaStream_t aes_async_stream;

CUDA_CALL(cudaStreamCreate(&aes_async_stream));

相對(duì)應(yīng)地會(huì)有銷毀流：

CUDA_CALL(cudaSteamDestroy(aes_async_stream));

然后將需要的事件壓入流中：

CUDA_CALL(cudaEventRecord(start_time,aes_async_stream));

CUDA_CALL(cudaMemcpyAsync(dev_xt,host_xt,size_of_xt_in_bytes,cudaMemcpyHostToDevice,aes_async_stream));

運(yùn)行時(shí)間如圖8所示。

圖8 GPU上使用CUDA流前后的時(shí)間對(duì)比

從圖8可以看出，在沒(méi)有使用CUDA流之前，時(shí)間穩(wěn)定在310 ms左右，使用CUDA流之后，時(shí)間穩(wěn)定在284 ms左右，時(shí)間減少26 ms，性能上提升9.15%。

4 結(jié)束語(yǔ)

本文針對(duì)MUSIC算法的高效GPU實(shí)現(xiàn)，經(jīng)過(guò)復(fù)數(shù)SVD的GPU并行計(jì)算，在譜峰搜索部分GPU實(shí)現(xiàn)中使用共享內(nèi)存、分配鎖頁(yè)內(nèi)存、增加CUDA流操作等幾個(gè)方面的改進(jìn)，整個(gè)計(jì)算過(guò)程減少56 ms的時(shí)間?？偟膩?lái)說(shuō)，性能提升了20%的效率。實(shí)驗(yàn)結(jié)果表明，性能的提升能夠符合預(yù)想的結(jié)果。