吳 超,衛(wèi) 謙,周俊偉,李會(huì)民,孫廣中

(1.中國(guó)科學(xué)技術(shù)大學(xué)網(wǎng)絡(luò)信息中心,安徽 合肥 230026;2.中國(guó)科學(xué)技術(shù)大學(xué)物理學(xué)院,安徽 合肥 230026;3.中國(guó)科學(xué)技術(shù)大學(xué)地球和空間科學(xué)學(xué)院,安徽 合肥 230026;4.中國(guó)科學(xué)技術(shù)大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,安徽 合肥 230027)

1 引言

傳統(tǒng)的地震學(xué)研究采集天然或人工震源產(chǎn)生的地震波信號(hào),利用地震波反演地質(zhì)結(jié)構(gòu)信息。地震臺(tái)站觀測(cè)的波形信號(hào)除了地震信號(hào),還包含大量的背景噪聲。通常認(rèn)為背景噪聲是無(wú)價(jià)值的信號(hào),在資料處理中應(yīng)予去除[1]。近年來(lái)隨著跨學(xué)科交叉研究的進(jìn)展,地球物理學(xué)家發(fā)現(xiàn)通過(guò)計(jì)算地震臺(tái)站間的背景噪聲互相關(guān)函數(shù)NCF(Noise Cross-correlation Function)可以提取出臺(tái)站間的格林函數(shù)[2],進(jìn)而研究地殼的速度、結(jié)構(gòu)等信息。21世紀(jì)以來(lái)基于背景噪聲的地震學(xué)研究得到蓬勃發(fā)展[3],相關(guān)成果已廣泛應(yīng)用于地球內(nèi)部結(jié)構(gòu)研究[4]、淺地表地質(zhì)調(diào)查[5]和油氣勘探等諸多領(lǐng)域。

然而,直接利用2個(gè)地震臺(tái)站的噪聲記錄進(jìn)行互相關(guān)計(jì)算,通常很難得到高信噪比的格林函數(shù)。這主要是由于噪聲的頻譜存在優(yōu)勢(shì)頻率,其較強(qiáng)的能量抑制了其它頻段的信號(hào)。為了減輕這種不利因素的影響,需要對(duì)原始背景噪聲記錄進(jìn)行預(yù)處理,通過(guò)時(shí)域和頻域信號(hào)處理,提升特定頻段的信噪比,之后進(jìn)行互相關(guān)計(jì)算獲得所需的信息[6]。

近年來(lái),世界主要大國(guó)紛紛加大密集臺(tái)陣的布設(shè)(如美國(guó)的USArray和中國(guó)地震科學(xué)臺(tái)陣探測(cè)計(jì)劃ChinArray)。這些觀測(cè)臺(tái)陣通常由成百上千個(gè)觀測(cè)臺(tái)站組成,對(duì)地震波形進(jìn)行長(zhǎng)年累月的觀測(cè)積累。密集的地震觀測(cè)臺(tái)站和長(zhǎng)程的時(shí)間積累記錄下海量的地震波形數(shù)據(jù),結(jié)合背景噪聲互相關(guān)計(jì)算,可以反演得到高精度的地下結(jié)構(gòu)圖像。與之對(duì)應(yīng)地,獲取密集臺(tái)站間NCF所需的計(jì)算耗時(shí)也快速增長(zhǎng)。為解決數(shù)據(jù)量激增帶來(lái)的計(jì)算壓力,研究人員嘗試數(shù)據(jù)并行路線,采用云計(jì)算加速地震數(shù)據(jù)處理,并取得了一定的成果[7]。隨著圖形處理器GPU(Graphics Processing Unit)、現(xiàn)場(chǎng)可編程邏輯門(mén)陣列FPGA(Field-Programmable Gate Array)等異構(gòu)計(jì)算設(shè)備的發(fā)展和日益成熟,異構(gòu)計(jì)算設(shè)備提供了遠(yuǎn)高于中央處理器CPU(Central Processing Unit)的數(shù)值計(jì)算能力,異構(gòu)并行計(jì)算技術(shù)為加速地震數(shù)據(jù)處理提供了新的技術(shù)選擇。

本文基于GPU異構(gòu)計(jì)算平臺(tái)設(shè)計(jì)了地震噪聲數(shù)據(jù)預(yù)處理的并行優(yōu)化方法。

2 相關(guān)工作

2.1 NCF計(jì)算方法

在背景噪聲地震學(xué)研究中,通常將臺(tái)站與臺(tái)站之間形成的“臺(tái)站對(duì)”稱為路徑。N個(gè)地震臺(tái)站兩兩之間形成的有效路徑(以臺(tái)站A和臺(tái)站B為例,臺(tái)站與自身的結(jié)對(duì)如路徑AA和BB去除,頭尾對(duì)換的2條路徑如AB和BA視為等同,只保留其一)數(shù)量為N(N-1)/2。NCF計(jì)算的主要過(guò)程分為3個(gè)步驟:

(1)預(yù)處理計(jì)算:對(duì)臺(tái)站單天的時(shí)域信號(hào)進(jìn)行去趨勢(shì)處理、帶通濾波、時(shí)域歸一化和譜白化等處理得到頻譜;

(2)互相關(guān)計(jì)算:對(duì)路徑在同一天的頻譜進(jìn)行共軛相乘并經(jīng)傅里葉反變換得到該路徑時(shí)間域的NCF;

(3)疊加計(jì)算:對(duì)路徑長(zhǎng)時(shí)間(根據(jù)尺度從小到大,范圍可以為周、月、年)的NCF進(jìn)行疊加求平均得到該路徑最終的互相關(guān)信息[8]。

3個(gè)步驟中互相關(guān)計(jì)算與預(yù)處理計(jì)算算法復(fù)雜,耗時(shí)比較高。預(yù)處理部分計(jì)算量正比于參與計(jì)算的臺(tái)站個(gè)數(shù)以及疊加的天數(shù),互相關(guān)部分計(jì)算量正比于參與計(jì)算的臺(tái)站個(gè)數(shù)的平方以及疊加的天數(shù)。

2.2 NCF并行優(yōu)化現(xiàn)狀

隨著密集地震臺(tái)陣的建設(shè)和地震波形數(shù)據(jù)長(zhǎng)時(shí)間的積累,計(jì)算耗時(shí)長(zhǎng)的問(wèn)題逐漸成為NCF應(yīng)用的一個(gè)主要阻礙。以ChinArray為例,據(jù)測(cè)算其臺(tái)陣674個(gè)臺(tái)站1年的記錄數(shù)據(jù)(降采樣到1 Hz)的NCF計(jì)算在1臺(tái)高性能服務(wù)器(Intel?雙路Xeon?E5-2620,64 GB內(nèi)存)上需要耗時(shí)3 840 h,接近6個(gè)月的時(shí)間[9]。為了提升數(shù)據(jù)處理的計(jì)算性能,研究人員嘗試采用不同的方法對(duì)NCF算法進(jìn)行優(yōu)化加速。文獻(xiàn)[10]使用云計(jì)算技術(shù)加速 NCF計(jì)算;文獻(xiàn)[11]基于分布式內(nèi)存并行計(jì)算機(jī)設(shè)計(jì)了NCF并行算法;文獻(xiàn)[12,13]分別實(shí)現(xiàn)了Julia語(yǔ)言和Python語(yǔ)言的NCF高性能計(jì)算方法;文獻(xiàn)[14]提出了CPU-GPU架構(gòu)超級(jí)計(jì)算機(jī)的NCF的高效計(jì)算方法;文獻(xiàn)[15]描述了在統(tǒng)一計(jì)算架構(gòu)CUDA(Compute Unified Device Architecture)上實(shí)現(xiàn)NCF計(jì)算的方法。上述研究都選擇NCF計(jì)算中計(jì)算復(fù)雜度最高的互相關(guān)計(jì)算部分進(jìn)行優(yōu)化,對(duì)預(yù)處理計(jì)算的加速優(yōu)化則沒(méi)有涉及。由于預(yù)處理算法的多樣性和可變性,在實(shí)際的工程實(shí)踐中預(yù)處理計(jì)算實(shí)施的次數(shù)通常遠(yuǎn)高于互相關(guān)計(jì)算;同時(shí),隨著高頻地震信號(hào)(如10 Hz～100 Hz)的測(cè)量和分析,預(yù)處理計(jì)算的數(shù)據(jù)規(guī)模也呈數(shù)量級(jí)增長(zhǎng)。因此,作為NCF計(jì)算的熱點(diǎn)問(wèn)題之一,預(yù)處理計(jì)算同樣具備性能優(yōu)化的必要性和價(jià)值。

2.3 CUDA并行計(jì)算

統(tǒng)一計(jì)算架構(gòu)CUDA是英偉達(dá)公司為推廣CPU+GPU異構(gòu)計(jì)算推出的軟件開(kāi)發(fā)包,支持C/C++、FORTRAN等多種編程語(yǔ)言。在CUDA開(kāi)發(fā)環(huán)境中,CPU和內(nèi)存屬于主機(jī)(host)端,GPU和顯存屬于設(shè)備(device)端。CUDA程序包含host端代碼和device端代碼,host端代碼以邏輯運(yùn)算為主運(yùn)行在CPU上,device端以數(shù)值計(jì)算為主運(yùn)行在GPU上。通過(guò)使用CUDA內(nèi)置的存儲(chǔ)API,實(shí)現(xiàn)數(shù)據(jù)在host和device之間的流動(dòng)。

CUDA在device端對(duì)GPU設(shè)備的線程進(jìn)行多層抽象,經(jīng)由網(wǎng)格(Grid)-線程塊(Block)-線程(Thread)3層結(jié)構(gòu)組織線程。通過(guò)規(guī)劃程序在host端和device端的計(jì)算任務(wù),合理使用Grid-Block-Thread組織結(jié)構(gòu),用戶可以充分利用CPU+GPU執(zhí)行異構(gòu)計(jì)算,獲得遠(yuǎn)超CPU串行計(jì)算的運(yùn)算能力。目前,CUDA并行計(jì)算已成功應(yīng)用于人工智能、科學(xué)計(jì)算和大數(shù)據(jù)等諸多應(yīng)用領(lǐng)域。

3 地震噪聲數(shù)據(jù)預(yù)處理串行算法

3.1 算法流程

單臺(tái)站單天的預(yù)處理計(jì)算讀取臺(tái)站記錄的波形文件作為輸入。首先,需對(duì)輸入的波形數(shù)據(jù)進(jìn)行全局去趨勢(shì)(Detrend)處理;然后,再將波形數(shù)據(jù)進(jìn)行分段(Segmentation,分段之間允許有交疊)。對(duì)于每個(gè)分段的處理,首先是對(duì)分段波形(Segment Beam)執(zhí)行去趨勢(shì)處理、時(shí)域歸一化(Time Normalization)處理等時(shí)域操作;然后,對(duì)分段波形進(jìn)行快速傅里葉變換FFT(Fast Fourier Transform)得到分段頻譜,并對(duì)其執(zhí)行譜白化(Spectrum Whitening)處理;最后,將所有分段頻譜(Segment Spectrum)拼接(Concatenate)得到完整頻譜,寫(xiě)入頻譜文件。預(yù)處理計(jì)算流程如圖1所示。

Figure 1 Serial preprocessing of single-day-beam-file of single station

多臺(tái)站多天的預(yù)處理計(jì)算以單臺(tái)站單天的處理為基礎(chǔ),分別對(duì)臺(tái)站和時(shí)間進(jìn)行兩重循環(huán)實(shí)現(xiàn)對(duì)所有文件的遍歷處理,如算法1所述。

算法1多臺(tái)站多天的串行預(yù)處理計(jì)算

輸入：多臺(tái)站(Nsta)多天(Nday)的波形數(shù)據(jù)D[]。

輸出：多臺(tái)站多天的頻譜數(shù)據(jù)S[]。

/* 遍歷臺(tái)站 */

1:fori=0 toNsta-1do

/* 遍歷天數(shù) */

2:forj=0 toNday-1do

/* 調(diào)用單臺(tái)站單天預(yù)處理 */

3:S[i,j]=PreProcess(D[i,j]);

4:endfor

5:endfor

3.2 熱點(diǎn)計(jì)算函數(shù)

3.2.1 去趨勢(shì)處理

地震波形數(shù)據(jù)常存在趨勢(shì)性變化,這種變化可能來(lái)自儀器本身的漂移,也可能是地殼蠕變的結(jié)果。趨勢(shì)性變化容易掩蓋細(xì)微的異常信息,因此在信號(hào)預(yù)處理中需要進(jìn)行去趨勢(shì)處理[16]。

對(duì)于時(shí)間序列(如地震波形數(shù)據(jù))的趨勢(shì)計(jì)算,常用的有3種計(jì)算方法:線性方法(基于最小二乘擬合Least Squares Fitting)、非線性方法(基于經(jīng)驗(yàn)?zāi)B(tài)分解Empirical Mode Decomposition)和一階差分法FD(First Differencing)[17]。地震背景噪聲預(yù)處理采用的是基于最小二乘法擬合的線性方法計(jì)算波形趨勢(shì)。

假設(shè)Dn(n=1,2,…,N)是采樣間隔為T(mén)的均勻采樣的波形數(shù)據(jù),采樣時(shí)間t=T,2T,…,NT。Dn包含一階趨勢(shì)項(xiàng)a0+a1t,記為式(1):

(1)

基于最小二乘方法擬合出趨勢(shì)項(xiàng),如式(2)所示:

(2)

其中,C(N,T)為確定系數(shù),可由求和公式、平方和公式、二階矩陣求逆得出解析解。由式(1)和式(2)可知,最小二乘法主要的計(jì)算量集中在求和∑Di(以及類似計(jì)算∑iDi)。樸素的串行求和計(jì)算如算法2所示。

算法2樸素串行求和計(jì)算

輸入：原始波形數(shù)據(jù)D[],波形點(diǎn)數(shù)Npts。

輸出：波形數(shù)據(jù)之和sum。

1:sum=0;

/* 遍歷波形數(shù)組求和 */

2:fori= 0 toNpts-1do

3:sum=sum+D[i];

4:endfor

求和計(jì)算屬于歸約算法。樸素的串行求和計(jì)算通過(guò)循環(huán)實(shí)現(xiàn),每次迭代將一項(xiàng)元素累加進(jìn)求和項(xiàng)。前Npts個(gè)元素之和依賴前Npts-1個(gè)元素之和,實(shí)現(xiàn)存在數(shù)據(jù)依賴關(guān)系。

根據(jù)式(2)計(jì)算趨勢(shì)項(xiàng)a0+a1t的估計(jì)值,波形數(shù)據(jù)逐項(xiàng)減去趨勢(shì)項(xiàng)的估計(jì)值實(shí)現(xiàn)去趨勢(shì)處理。

3.2.2 時(shí)域歸一化處理

數(shù)據(jù)預(yù)處理中最重要的一步是時(shí)間域歸一化處理,其目的是為了去除地震信號(hào)、儀器故障引起的畸變信號(hào)以及地震臺(tái)站附近噪聲對(duì)互相關(guān)計(jì)算結(jié)果的影響。

目前,常用的時(shí)間域歸一化方法有2種:one-bit方法和滑動(dòng)絕對(duì)平均方法。one-bit方法遍歷波形原始記錄,將正值替換為+1,負(fù)值替換為-1?；瑒?dòng)絕對(duì)平均方法則是根據(jù)式(3)計(jì)算一定時(shí)窗內(nèi)波形數(shù)據(jù)絕對(duì)振幅的平均值,即該時(shí)窗中心點(diǎn)的權(quán)重,逐段移動(dòng)時(shí)窗,根據(jù)式(4)每點(diǎn)的幅值除以權(quán)重得到新的波形序列。

(3)

(4)

窗的長(zhǎng)度2k+1決定了有多少振幅信息可以保留。當(dāng)k=0時(shí),滑動(dòng)絕對(duì)平均方法等效于one-bit方法。為計(jì)算全部波形數(shù)據(jù)的時(shí)域歸一化,下標(biāo)n從0遍歷到Npts-1,在計(jì)算過(guò)程中當(dāng)j超過(guò)波形數(shù)組索引范圍時(shí)(如j<0或j>Npts-1),超出部分的dj值設(shè)置為0進(jìn)行計(jì)算。

本文預(yù)處理計(jì)算使用滑動(dòng)絕對(duì)平均時(shí)域歸一化方法,如算法3所示。

算法3滑動(dòng)絕對(duì)平均時(shí)域歸一化串行算法

輸入：原始波形數(shù)據(jù)D[],波形點(diǎn)數(shù)Npts,窗長(zhǎng)K。

輸出：時(shí)域歸一化波形數(shù)據(jù)D[]。

1:k=[K/2];

2:fori=0 toNpts-1do

3:temp=0;

4:forj=0 to 2kdo

5:ifi-k+j≥0 andi-k+j

6:temp=temp+fabs(D[i-k+j]);

7:endif

8:endfor

9:temp=temp/K;

10:D[i]=D[i]/temp;

11:endfor

3.2.3 譜白化處理

譜白化處理利用快速傅里葉變換和反變換,將時(shí)域波形的頻譜進(jìn)行分頻,對(duì)分頻波形進(jìn)行時(shí)變?cè)鲆?最后重新合成獲得白化信號(hào)[18]。其中的主要計(jì)算是快速傅里葉變換?？焖俑道锶~變換是快速計(jì)算序列的離散傅里葉變換或其逆變換的方法。作為信號(hào)處理領(lǐng)域的核心算法,已有眾多文獻(xiàn)討論FFT算法及實(shí)現(xiàn),FFTW庫(kù)[19]是其中影響較為廣泛的一個(gè)實(shí)現(xiàn)。本文的串行預(yù)處理計(jì)算調(diào)用FFTW庫(kù)執(zhí)行快速傅里葉變換。

3.3 性能剖析

本節(jié)使用Linux kernel中的系統(tǒng)性能剖析工具Perf對(duì)串行預(yù)處理計(jì)算程序典型計(jì)算場(chǎng)景(1 000個(gè)波形數(shù)據(jù)文件,1 Hz數(shù)據(jù)采樣率)的運(yùn)行時(shí)情況進(jìn)行性能剖析,各計(jì)算模塊耗時(shí)比例如圖2所示。

Figure 2 Perf profiling result of serial preprocessing

根據(jù)串行預(yù)處理程序的性能剖析,主要耗時(shí)的計(jì)算模塊為去趨勢(shì)處理、時(shí)域歸一化處理和快速傅里葉變換(譜白化處理)。其中,占比最高的部分為快速傅里葉變換(FFT),達(dá)到總耗時(shí)的84%;占比第2的是時(shí)域歸一化處理(Time Normalization),占總耗時(shí)的8%;占比第3的是去趨勢(shì)處理(Detrend),占比4%。3個(gè)計(jì)算模塊總耗時(shí)超過(guò)95%。I/O時(shí)間約占程序整體運(yùn)行時(shí)間的3%左右。

Amdahl定律(Amdahl’s Law)指出:固定計(jì)算工作負(fù)載的前提下,并行系統(tǒng)隨著處理器數(shù)目無(wú)限增大所能達(dá)到的加速比上限為1/Ts,其中Ts為程序中串行分量所占比重[20]。由Amdahl定律可知,僅針對(duì)預(yù)處理計(jì)算耗時(shí)最多的快速傅里葉變換模塊進(jìn)行并行加速,并行系統(tǒng)的加速比不會(huì)超過(guò)8倍。而如果對(duì)耗時(shí)前三的計(jì)算模塊進(jìn)行并行改造,則系統(tǒng)的加速比上限可以達(dá)到或接近33倍。為了實(shí)現(xiàn)海量波形文件的準(zhǔn)實(shí)時(shí)預(yù)處理,盡可能提高系統(tǒng)的加速效果,本文在后續(xù)章節(jié)中將對(duì)預(yù)處理程序的所有計(jì)算模塊進(jìn)行基于CUDA的移植和優(yōu)化。

4 地震噪聲數(shù)據(jù)預(yù)處理算法的并行優(yōu)化

4.1 分批處理

多臺(tái)站多天的預(yù)處理計(jì)算是數(shù)據(jù)密集型任務(wù)。密集臺(tái)陣長(zhǎng)時(shí)間積累的波形數(shù)據(jù)容量往往遠(yuǎn)超單個(gè)GPU卡的顯存容量,也超出了單臺(tái)服務(wù)器的內(nèi)存容量。在實(shí)際生產(chǎn)環(huán)境中,針對(duì)計(jì)算服務(wù)器數(shù)量有限,需要使用單臺(tái)服務(wù)器完成密集臺(tái)陣波形數(shù)據(jù)的預(yù)處理任務(wù)的情況,通過(guò)采用分批處理的策略完成長(zhǎng)時(shí)間積累的波形數(shù)據(jù)在單服務(wù)器上的任務(wù)計(jì)算。分批處理流程如圖3所示。

Figure 3 Batch CUDA preprocessing of multi-day-beam-files of multiple stations

分批處理使用二重循環(huán)實(shí)現(xiàn)。外層循環(huán)實(shí)現(xiàn)的是硬盤(pán)與主機(jī)內(nèi)存的數(shù)據(jù)(波形數(shù)據(jù),算法輸入為時(shí)域波形,輸出為頻域波形)傳輸,在外層迭代計(jì)算之前和之后分別執(zhí)行文件的輸入和輸出I/O操作。內(nèi)存循環(huán)實(shí)現(xiàn)的是主機(jī)內(nèi)存與設(shè)備內(nèi)存的數(shù)據(jù)傳輸,在內(nèi)層迭代之前將數(shù)據(jù)從主機(jī)內(nèi)存?zhèn)鬏數(shù)皆O(shè)備內(nèi)存,在內(nèi)層迭代之后將數(shù)據(jù)從設(shè)備內(nèi)存?zhèn)鬏數(shù)街鳈C(jī)內(nèi)存。對(duì)于一次計(jì)算任務(wù),需要在主機(jī)內(nèi)存和設(shè)備內(nèi)存之間傳輸?shù)目倲?shù)據(jù)量是確定不變的,傳輸次數(shù)是影響數(shù)據(jù)傳輸開(kāi)銷的主要因素,傳輸次數(shù)越多則傳輸開(kāi)銷越大。為了減少CPU和GPU之間數(shù)據(jù)傳輸次數(shù),分批處理進(jìn)行迭代的依據(jù)是每批處理的波形文件批數(shù)以盡量占滿主機(jī)內(nèi)存(或設(shè)備內(nèi)存)容量為準(zhǔn),即外層循環(huán)的單次迭代處理批數(shù)由主機(jī)內(nèi)存可容納的最大波形文件數(shù)確定,內(nèi)存循環(huán)的單次迭代處理批數(shù)由設(shè)備內(nèi)存可容納的最大波形文件數(shù)確定。

在對(duì)波形文件進(jìn)行批數(shù)劃分時(shí),采用針對(duì)時(shí)間維度的數(shù)據(jù)劃分,設(shè)計(jì)主機(jī)內(nèi)存、設(shè)備內(nèi)存分層的批處理方法。首先,估算所有臺(tái)站的單天數(shù)據(jù)預(yù)處理需要的主機(jī)內(nèi)存和設(shè)備內(nèi)存空間。根據(jù)服務(wù)器和GPU卡的主機(jī)內(nèi)存和設(shè)備內(nèi)存容量,估算主機(jī)內(nèi)存和設(shè)備內(nèi)存單次最多可容納計(jì)算的天數(shù)分別為T(mén)M和TV(一般有TM>TV)。當(dāng)遇到主機(jī)內(nèi)存容量小于設(shè)備內(nèi)存容量(TM

4.2 多流處理

在4.1節(jié)所述的分批處理流程中,CPU與GPU通過(guò)主機(jī)內(nèi)存和設(shè)備內(nèi)存之間的內(nèi)存拷貝函數(shù)實(shí)現(xiàn)隱式同步。首先,CPU將GPU需要并行處理的批量波形數(shù)據(jù)通過(guò)cudaMemcpy函數(shù)從主機(jī)內(nèi)存?zhèn)鬏數(shù)皆O(shè)備內(nèi)存;之后,GPU針對(duì)波形數(shù)據(jù)進(jìn)行并行處理,得到批量頻譜數(shù)據(jù);最后,GPU將批量頻譜數(shù)據(jù)通過(guò)cudaMemcpy函數(shù)從設(shè)備內(nèi)存?zhèn)骰刂鳈C(jī)內(nèi)存,完成一輪GPU計(jì)算。

CUDA流(Stream)也稱異步流,表示GPU操作隊(duì)列。同一個(gè)流內(nèi)操作的執(zhí)行順序是有序的,但每個(gè)操作的執(zhí)行開(kāi)始時(shí)間是不可控的,并且不同流的執(zhí)行順序也不能顯式控制[21]。為了提升預(yù)處理計(jì)算的執(zhí)行效率,通過(guò)利用GPU的多流處理能力,將批量處理的數(shù)據(jù)平均分配到GPU多個(gè)流上。不同流數(shù)據(jù)從主機(jī)內(nèi)存到設(shè)備內(nèi)存的傳輸、結(jié)果從設(shè)備內(nèi)存到主機(jī)內(nèi)存的傳輸與GPU端的核函數(shù)計(jì)算可進(jìn)行并行處理,實(shí)現(xiàn)傳輸延遲的隱藏,減少傳輸和計(jì)算的耗時(shí)。多流處理流程如圖4所示。

Figure 4 Workflow of multi-stream processing

4.3 并行計(jì)算框架

經(jīng)由串行程序的靜態(tài)分析可知,預(yù)處理計(jì)算的計(jì)算任務(wù)主要集中在分段波形數(shù)據(jù)的循環(huán)處理,循環(huán)中每次迭代處理一個(gè)分段的波形數(shù)據(jù)。由于循環(huán)內(nèi)每個(gè)分段的計(jì)算獨(dú)立,分段之間不存在數(shù)據(jù)依賴或數(shù)據(jù)競(jìng)爭(zhēng),所以單臺(tái)站單天預(yù)處理計(jì)算的整個(gè)循環(huán)可以使用CUDA的一維Grid結(jié)構(gòu)并行處理,每個(gè)分段波形數(shù)據(jù)分配1個(gè)Block完成處理。多臺(tái)站多天的預(yù)處理,則是在臺(tái)站維和時(shí)間維進(jìn)行拓展,使用三維Grid結(jié)構(gòu)完成計(jì)算,如圖5所示。Grid結(jié)構(gòu)X軸、Y軸、Z軸3個(gè)維度分別對(duì)應(yīng)預(yù)處理計(jì)算處理的分段、臺(tái)站及時(shí)間。

Figure 5 Structure of 3D Grid

單臺(tái)站單天預(yù)處理計(jì)算整體的算法如算法4所示。

算法4基于CUDA的預(yù)處理并行算法

輸入：?jiǎn)闻_(tái)站單天原始波形數(shù)據(jù)D[]。

輸出：?jiǎn)闻_(tái)站單天頻域數(shù)據(jù)S[]。

Step1H2D(Host2Device):波形數(shù)據(jù)D從內(nèi)存拷貝到顯存。

Step2分配1個(gè)Block完成全局波形數(shù)據(jù)的去趨勢(shì)處理。

Step3全局波形數(shù)據(jù)切分成L段,每個(gè)分段分配1個(gè)Block進(jìn)行處理。

Step4遍歷k從0到L-1,并行處理:

Step4.1Block(k)完成第k段波形的去趨勢(shì)處理、帶通濾波和時(shí)域歸一化處理;

Step4.2第k段波形執(zhí)行快速傅里葉變換,獲得分段頻譜;

Step4.3Block(k)完成第k段波形的譜白化處理。

Step5所有分段頻譜拼接成全局頻譜。

Step6D2H(Device2Host):頻譜數(shù)據(jù)S從顯存拷貝到內(nèi)存。

4.4 熱點(diǎn)函數(shù)的CUDA實(shí)現(xiàn)

基于CUDA實(shí)現(xiàn)高效的求和計(jì)算是實(shí)現(xiàn)去趨勢(shì)處理的關(guān)鍵。傳統(tǒng)的基于相鄰元素的歸約樹(shù)并行實(shí)現(xiàn),在相鄰的線程之間導(dǎo)致分支發(fā)散(Branch Divergence),性能不佳[22]。通過(guò)基于可變間隔的規(guī)約計(jì)算,以及運(yùn)用共享內(nèi)存減少訪存開(kāi)銷,可實(shí)現(xiàn)線程塊內(nèi)高效的并行求和計(jì)算,如算法5所示。

算法5單Block求和CUDA核函數(shù)

輸入：原始波形數(shù)據(jù)d_data[]和波形點(diǎn)數(shù)Npts。

輸出：波形數(shù)據(jù)之和d_sum。

1:inttx=threadIdx.x;

2:inti,stride;

3:doublesum=0;

4:fori=tx;i

5:sum+=d_data[i];

6:endfor

7:extern_shared_doublepartial[];

8:partial[tx]=sum;

9:_syncthreads();

10:forstride=(blockDim.x+1)/2;stride≥1;stride=(stride+1)/2do

11:iftx+stride

12:partial[tx]+=partial[tx+stride];

13:partial[tx+stride]=0;

14:endif

15: _syncthreads();

16:endfor

17:iftx==0then

18: *d_sum=partial[0];

19:endif

基于滑動(dòng)絕對(duì)平均方法的時(shí)域歸一化處理中波形數(shù)據(jù)格點(diǎn)的更新是根據(jù)周?chē)顸c(diǎn)在時(shí)間上加權(quán)求和得出,本質(zhì)是一維Stencil計(jì)算。Stencil計(jì)算是訪存密集型計(jì)算,使用共享內(nèi)存緩存格點(diǎn)數(shù)據(jù)優(yōu)化程序的訪存方式?；诠蚕韮?nèi)存的時(shí)域歸一化處理核函數(shù)算法如算法6所示。

算法6滑動(dòng)絕對(duì)平均時(shí)域歸一化核函數(shù)

輸入：分段波形數(shù)據(jù)d_data[],分段波形點(diǎn)數(shù)Nsegpts,總波形點(diǎn)數(shù)Npts,窗長(zhǎng)K。

輸出：時(shí)域歸一化分段波形數(shù)據(jù)d_data[]。

1:inti,j,idx,start;

2:idx=blockIdx.x*blockDim.x+threadIdx.x;

3:tx=threadIdx.x;

4:blkstart=blockIdx.x*blockDim.x;

5:nextblkstart=(blockIdx.x+1)*blockDim.x;

6:_shared_doubled_share[Nsegpts];

7:d_share[tx]=d_data[idx];

8:_syncthreads();

9:doubletmp=0;

10:start=idx-K/2;

11:fori=0;i

12:j=start+i;

13:if(j≥0) and (j

14:if(j≥blkstart) and (j

15:tmp+=fabs(d_share[tx]);

16:else

17:tmp+=fabs(d_data[idx]);

18:endif

19:endif

20:endfor

21:d_data[idx]=tmp/K;

5 實(shí)驗(yàn)及結(jié)果分析

5.1 實(shí)驗(yàn)平臺(tái)

本文的實(shí)驗(yàn)平臺(tái)硬件配置如下:CPU為 Intel?CoreTMI5-7500,架構(gòu)為Kaby Lake,主頻為3.4 GHz,核心數(shù)為4,內(nèi)存大小為16 GB;GPU使用的是NVIDIA的GeForce?RTX 2080Ti,4 352個(gè)CUDA核心,雙精度浮點(diǎn)性能達(dá)到420.2 GFlops,顯存大小為11 GB,顯存帶寬可達(dá)616.0 GB/s。

5.2 實(shí)驗(yàn)方法及數(shù)據(jù)

實(shí)驗(yàn)方法是針對(duì)預(yù)處理計(jì)算的主要計(jì)算模塊耗時(shí)和程序整體耗時(shí),隨著處理波形文件數(shù)量的變換,呈現(xiàn)CPU與GPU之間的運(yùn)行時(shí)間和加速比。加速比為CPU運(yùn)行時(shí)間除以GPU運(yùn)行時(shí)間。測(cè)試數(shù)據(jù)為美國(guó)地震觀測(cè)臺(tái)陣USArray采集的部分波形數(shù)據(jù)。每個(gè)波形文件記錄了單個(gè)臺(tái)站單天采集到的地震信號(hào),數(shù)據(jù)采樣率為1 Hz,數(shù)據(jù)點(diǎn)數(shù)為86 400,共有10 000個(gè)波形文件。

5.3 實(shí)驗(yàn)結(jié)果

本文利用NVIDIA公司提供的性能剖析工具中的nvprof軟件[23],剖析了基于GPU的并行預(yù)處理計(jì)算各主要計(jì)算模塊的運(yùn)行時(shí)性能。圖6～圖8分別記錄了快速傅里葉變換、時(shí)域歸一化處理和去趨勢(shì)處理基于CPU的串行和基于GPU的并行版本針對(duì)不同數(shù)量波形文件的計(jì)算耗時(shí)。

Figure 6 Runtime of FFT module

Figure 7 Runtime of normalization module

Figure 8 Runtime of detrend module

快速傅里葉變換模塊是串行程序中耗時(shí)占比第一的模塊。針對(duì)數(shù)據(jù)規(guī)模從1 000到5 000個(gè)SAC文件,基于GPU的并行版本相比基于CPU的串行版本加速139～203倍。

時(shí)域歸一化處理模塊是串行程序中耗時(shí)占比第2的模塊。針對(duì)數(shù)據(jù)規(guī)模從1 000到5 000個(gè)SAC文件,基于GPU的并行版本相比基于CPU的串行版本加速25～30倍。

去趨勢(shì)處理模塊是串行程序中耗時(shí)占比第3的模塊。針對(duì)數(shù)據(jù)規(guī)模從1 000到5 000個(gè)SAC文件,基于GPU的并行版本相比基于CPU的串行版本加速58～72倍。

基于CPU的串行預(yù)處理計(jì)算、基于GPU的并行預(yù)處理計(jì)算(單流處理)的總處理時(shí)間實(shí)驗(yàn)結(jié)果如表1所示。表中數(shù)據(jù)是對(duì)10次實(shí)驗(yàn)取平均得到的結(jié)果。

Table 1 Total time costs of preprocessing algorithms

以表1中CPU串行預(yù)處理計(jì)算耗時(shí)為基準(zhǔn),計(jì)算預(yù)處理計(jì)算GPU單流并行和GPU多流并行的加速比,結(jié)果如圖9所示。

Figure 9 Speedup of parallel preprocessing algorithms

根據(jù)實(shí)驗(yàn)結(jié)果可知:(1) 基于GPU的并行預(yù)處理算法性能顯著優(yōu)于基于CPU的串行預(yù)處理算法。(2) 隨著波形文件數(shù)據(jù)量的增加,基于GPU的并行預(yù)處理算法的加速比呈線性增長(zhǎng),當(dāng)文件數(shù)量達(dá)到5 000時(shí),加速比達(dá)到最大,約為95.27倍;當(dāng)文件數(shù)量超過(guò)5 000時(shí),基于GPU的并行預(yù)處理算法的加速比略有下降。原因分析如下:文件數(shù)量達(dá)到5 000時(shí),GPU計(jì)算能力得到充分利用,GPU計(jì)算時(shí)間與文件I/O時(shí)間達(dá)到均衡,加速比達(dá)到峰值。之后再增加文件數(shù)量,由于I/O時(shí)間的影響,所以并行系統(tǒng)整體加速比下降。(3)使用多流處理技術(shù),可以隱藏部分?jǐn)?shù)據(jù)內(nèi)存?zhèn)鬏數(shù)难舆t,在單流GPU并行加速的基礎(chǔ)上可以進(jìn)一步獲得平均約8%的性能優(yōu)化效果。

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

對(duì)比傳統(tǒng)的基于CPU的串行地震噪聲預(yù)處理算法,本文提出的基于GPU的并行算法顯著提升了計(jì)算性能,且具有較好的并行度。通過(guò)實(shí)驗(yàn)仿真可知,基于GPU的并行預(yù)處理算法適用于處理地震波形文件數(shù)量較多的場(chǎng)景,為地震數(shù)據(jù)預(yù)處理的GPU加速提供了研究思路。

當(dāng)前,針對(duì)GPU并行模型在地震數(shù)據(jù)預(yù)處理算法中應(yīng)用的研究文獻(xiàn)較少。通過(guò)本文實(shí)驗(yàn)可知,并行模型對(duì)包含大量信號(hào)處理計(jì)算的地震預(yù)處理算法具有良好的加速效果。由于分批處理各批量之間無(wú)數(shù)據(jù)相關(guān)性,本文算法非常適于擴(kuò)展為分布式版本。下一步將研究MapReduce分布式計(jì)算框架與本文加速方法相結(jié)合的策略,以及研究單計(jì)算節(jié)點(diǎn)多GPU卡并行計(jì)算及多計(jì)算節(jié)點(diǎn)并行計(jì)算在地震噪聲預(yù)處理算法中的應(yīng)用。

致謝

感謝中國(guó)地震局地球物理研究所王偉濤研究員和王芳博士在本文工作完成過(guò)程中給予的指導(dǎo)和幫助。