基于異構計算平臺的背景噪聲預處理并行算法*

吳 超,衛 謙,周俊偉,李會民,孫廣中

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

1 引言

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

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

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

本文基于GPU異構計算平臺設計了地震噪聲數據預處理的并行優化方法。

2 相關工作

2.1 NCF計算方法

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

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

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

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

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

2.2 NCF并行優化現狀

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

2.3 CUDA并行計算

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

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

3 地震噪聲數據預處理串行算法

3.1 算法流程

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

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

多臺站多天的預處理計算以單臺站單天的處理為基礎,分別對臺站和時間進行兩重循環實現對所有文件的遍歷處理,如算法1所述。

算法1多臺站多天的串行預處理計算

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

輸出：多臺站多天的頻譜數據S[]。

/* 遍歷臺站 */

1:fori=0 toNsta-1do

/* 遍歷天數 */

2:forj=0 toNday-1do

/* 調用單臺站單天預處理 */

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

4:endfor

5:endfor

3.2 熱點計算函數

3.2.1 去趨勢處理

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

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

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

(1)

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

(2)

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

算法2樸素串行求和計算

輸入：原始波形數據D[],波形點數Npts。

輸出：波形數據之和sum。

1:sum=0;

/* 遍歷波形數組求和 */

2:fori= 0 toNpts-1do

3:sum=sum+D[i];

4:endfor

求和計算屬于歸約算法。樸素的串行求和計算通過循環實現,每次迭代將一項元素累加進求和項。前Npts個元素之和依賴前Npts-1個元素之和,實現存在數據依賴關系。

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

3.2.2 時域歸一化處理

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

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

(3)

(4)

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

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

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

輸入：原始波形數據D[],波形點數Npts,窗長K。

輸出：時域歸一化波形數據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]。其中的主要計算是快速傅里葉變換?？焖俑道锶~變換是快速計算序列的離散傅里葉變換或其逆變換的方法。作為信號處理領域的核心算法,已有眾多文獻討論FFT算法及實現,FFTW庫[19]是其中影響較為廣泛的一個實現。本文的串行預處理計算調用FFTW庫執行快速傅里葉變換。

3.3 性能剖析

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

Figure 2 Perf profiling result of serial preprocessing

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

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

4 地震噪聲數據預處理算法的并行優化

4.1 分批處理

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

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

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

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

4.2 多流處理

在4.1節所述的分批處理流程中,CPU與GPU通過主機內存和設備內存之間的內存拷貝函數實現隱式同步。首先,CPU將GPU需要并行處理的批量波形數據通過cudaMemcpy函數從主機內存傳輸到設備內存;之后,GPU針對波形數據進行并行處理,得到批量頻譜數據;最后,GPU將批量頻譜數據通過cudaMemcpy函數從設備內存傳回主機內存,完成一輪GPU計算。

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

Figure 4 Workflow of multi-stream processing

4.3 并行計算框架

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

Figure 5 Structure of 3D Grid

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

算法4基于CUDA的預處理并行算法

輸入：單臺站單天原始波形數據D[]。

輸出：單臺站單天頻域數據S[]。

Step1H2D(Host2Device):波形數據D從內存拷貝到顯存。

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

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

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

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

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

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

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

Step6D2H(Device2Host):頻譜數據S從顯存拷貝到內存。

4.4 熱點函數的CUDA實現

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

算法5單Block求和CUDA核函數

輸入：原始波形數據d_data[]和波形點數Npts。

輸出：波形數據之和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

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

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

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

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

5.1 實驗平臺

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

5.2 實驗方法及數據

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

5.3 實驗結果

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

Figure 6 Runtime of FFT module

Figure 7 Runtime of normalization module

Figure 8 Runtime of detrend module

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

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

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

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

Table 1 Total time costs of preprocessing algorithms

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

Figure 9 Speedup of parallel preprocessing algorithms

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

6 結束語

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

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

致謝

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