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

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

1 引言

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

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

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

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

2 相關(guān)工作

2.1 NCF計算方法

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

(1)預(yù)處理計算:對臺站單天的時域信號進行去趨勢處理、帶通濾波、時域歸一化和譜白化等處理得到頻譜;

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

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

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

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

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

2.3 CUDA并行計算

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

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

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

3.1 算法流程

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

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

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

算法1多臺站多天的串行預(yù)處理計算

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

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

/* 遍歷臺站 */

1:fori=0 toNsta-1do

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

2:forj=0 toNday-1do

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

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

4:endfor

5:endfor

3.2 熱點計算函數(shù)

3.2.1 去趨勢處理

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

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

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

(1)

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

(2)

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

算法2樸素串行求和計算

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

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

1:sum=0;

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

2:fori= 0 toNpts-1do

3:sum=sum+D[i];

4:endfor

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

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

3.2.2 時域歸一化處理

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

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

(3)

(4)

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

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

算法3滑動絕對平均時域歸一化串行算法

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

輸出：時域歸一化波形數(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 譜白化處理

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

3.3 性能剖析

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

Figure 2 Perf profiling result of serial preprocessing

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

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

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

4.1 分批處理

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

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

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

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

4.2 多流處理

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

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

Figure 4 Workflow of multi-stream processing

4.3 并行計算框架

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

Figure 5 Structure of 3D Grid

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

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

輸入：單臺站單天原始波形數(shù)據(jù)D[]。

輸出：單臺站單天頻域數(shù)據(jù)S[]。

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

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

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

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

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

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

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

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

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

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

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

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

輸入：原始波形數(shù)據(jù)d_data[]和波形點數(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

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

算法6滑動絕對平均時域歸一化核函數(shù)

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

輸出：時域歸一化分段波形數(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 實驗及結(jié)果分析

5.1 實驗平臺

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

5.2 實驗方法及數(shù)據(jù)

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

5.3 實驗結(jié)果

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

Figure 6 Runtime of FFT module

Figure 7 Runtime of normalization module

Figure 8 Runtime of detrend module

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

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

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

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

Table 1 Total time costs of preprocessing algorithms

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

Figure 9 Speedup of parallel preprocessing algorithms

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

6 結(jié)束語

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

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

致謝

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