基于CUDA的地震信號頻率補償并行算法

張 全 張杰明 雷 芩 彭 博*② 劉書妍②

(①西南石油大學計算機科學學院，四川成都 610500； ②油氣藏地質及開發工程國家重點實驗室(西南石油大學)，四川成都 610500； ③電子科技大學信息與通信工程學院，四川成都611731)

0 引言

隨著油氣勘探的不斷發展，低分辨率的地震數據無法揭示復雜構造的細節信息，難以滿足當今油氣勘探需求。提高地震資料的分辨率是地球物理領域的研究熱點之一。

眾所周知，地震資料的低頻成分可反映巖性變化，決定了基本的速度結構[1]。豐富的低頻信息可使儲層反演結果更清晰、可靠[2]，有利于巖性油氣藏的識別[3]。地震信號的頻帶寬度影響地震資料的分辨率[1]，但在實際采集過程中，地震波受地層吸收衰減、檢波器畸變等因素的影響，往往導致高、低頻信息的缺失。低頻信息的缺失會導致子波旁瓣增大，地震記錄出現假同相軸等[4-5]，此外，它會影響特殊構造的成像，造成深層構造不清晰[6]； 高頻信息的缺失會導致資料分辨率下降，地層細節顯示不清。針對這些問題，許多學者在頻率補償方面進行了研究。頻率補償傳統方法有反褶積、反Q濾波和譜均衡等，但都有各自的局限性。反褶積法可以恢復低頻成分，但同時會增強有效頻段外的低頻噪聲[7]。反Q濾波法存在精確估計Q值較困難的問題，Tonn[8]對常見的7種Q值計算方法進行了比較，發現沒有哪一種方法可適用于普遍情況，其效果依賴于地震資料的質量。近年來，一些學者運用壓縮感知[9-11]稀疏重構的思想重構地震信號，再用重構的信號恢復原始信號的頻率信息[12-13]。韓立國等[14]、張瑩[15]根據壓縮感知和稀疏約束的原理對地震信號進行全頻帶補償，全面拓寬了地震頻帶。丁燕等[16]用改進的寬帶俞氏低通濾波器，對地震信號進行了低頻補償，解決了頻譜疊加信號的失真問題。

雖然對地震信號頻率補償算法進行了改進，但面對龐大的地震數據，該算法計算效率仍然不高。大規模的地震數據處理屬于計算密集型任務，隨著高性能計算的發展，運用多核CPU或者運用圖形處理器(Graphics Prossessing Unit, GPU)進行地震資料處理已經成為時代的趨勢。張軍華等[17]分析了高性能計算在地震勘探行業的發展現狀； 趙虎等[18]運用OpenACC編程模型對逆時偏移算法進行并行加速； 張全等[19]在CPU/GPU異構平臺上對拋物線拉東變換與地震相干體算法進行并行加速； 張全等[20]、吳吉忠等[21]將高性能計算運用在地球物理勘探領域，有效提高了地震資料處理的效率。

本文基于CUDA(Compute Unified Device Ar-chitecture)平臺，對基于壓縮感知的頻率補償方法進行了分析和并行化實現。

1 頻率補償算法原理

根據地震褶積模型，地震信號s可描述為地震子波w和地下反射系數r的卷積與隨機噪聲n的和

s=w*r+n

(1)

式中“*”表示卷積。其頻率域表示形式為

S=WR+N=WFr+N

(2)

式中：S、W、R、N分別為地震信號、地震子波、反射系數、噪聲的傅里葉變換； F為傅里葉變換算子[22]。

由于地震子波在傳播途中受到大地濾波的作用，其頻帶是有限帶寬的。地震子波在與反射系數卷積的過程中，必然使反射系數的頻譜成分受到損失，最后得到的地震記錄在頻域也是有限帶寬的，在時域表現為分辨率不高。為得到更高分辨率的記錄，需估計出地震子波和反射系數，并對地震子波頻帶進行擴展。最后用寬頻子波與地層反射系數卷積得到寬頻的、高分辨率的地震記錄。

該算法處理對象為動校正后的疊前道集數據，整個算法處理流程主要包括如下六個步驟。

(1)統計子波。假設反射系數為白噪聲，子波的自相關近似等于地震記錄的自相關[23-24]。根據地震記錄的自相關可估計疊加道子波和隨炮檢距變化的子波集SW。

(2)計算疊加道。疊加道集中的所有道并求平均值，在一定程度上消除噪聲的影響。

(3)求解疊加道反射系數位置。為求解疊加道地層反射系數r，構建如下目標函數

(3)

(4)構建子波矩陣。根據反射系數位置可構建隨炮檢距變化的子波矩陣。每一道地震記錄對應的子波應當從自相關統計的子波集SW中選取。假設一個地震道集第m道地震數據為dm，選取SW中第m個子波，根據步驟(3)獲得的反射系數位置序列P，在強反射位置處為子波矩陣賦值，構建子波矩陣M(圖1)。

圖1 構造隨炮檢距變化子波矩陣示意圖

(5)計算每一道的反射系數。根據地震數據褶積模型，地震子波與反射系數卷積為線性變換過程，可以表示為子波矩陣與反射系數相乘的形式s=Mr(M可看作Toeplitz矩陣的特定列；r為反射系數除去零值后的序列，將卷積轉化為矩陣相乘的形式便于計算，且利用反射系數的稀疏性節省了計算量)，利用共軛梯度法[26-28]迭代求解可得到每一道的反射系數序列r，其步驟如下。

1)任取r(0)∈Rn計算f(0)=s-Mr(0)，取p(0)=f(0)。

2)對于k=1,2,…，計算

(4)

r(k+1)=r(k)+αkp(k)

(5)

f(k+1)=f(k)-αkMp(k)

(6)

(7)

p(k+1)=f(k+1)+βkp(k)

(8)

式中：k為迭代次數；f(k)為第k次迭代過程中的梯度向量；p(k)為第k次迭代過程中算法搜索方向；αk、βk分別為第k次迭代過程中的迭代梯度方向步長和搜索方向步長。

3)當f(k)=0或小于誤差限時停止計算，r(k)為反射系數的解。

(6)將反射系數與寬頻子波褶積得到拓頻后的道集。

2 頻率補償算法的并行實現

根據上述算法的計算過程，首先對算法的每一個步驟進行分段計時，分析算法的計算密集部分，然后針對該部分設計對應的并行方案。

2.1 計算瓶頸分析

本文以道集為單位對CPU代碼進行分段計時測試。測試環境：CPU為Intel Core(TM) i5900H； 主頻為2.40GHz； 內存為8GB； 操作系統為64位Windows10。測試數據為三維工區CRP道集，主測線范圍為2201～2401，聯絡測線范圍為1500～1700，每個道集有51道，每道有376個樣點，采樣間隔為4ms。測試結果為100次計時結果的平均值。主函數平均用時為183ms，步驟(1)～步驟(6)的用時和用時占比如表1所示。

由表1可知，步驟(4)～步驟(5)為算法相對耗時的部分，為并行改進的重點。

表1 各步驟耗時及占比

2.2 并行算法流程

步驟(4)、步驟(5)為算法最耗時的部分，也是本文算法改進的重點。在步驟(4)構造子波矩陣時，每一道的子波矩陣互相獨立，涉及到大量賦值操作； 在步驟(5)計算反射系數時，每一道反射系數的計算相互獨立，運用共軛梯度法求解需要大量矩陣運算，矩陣運算天然地適用于并行。

由于CPU和GPU不能直接互相訪問各自的存儲器，當GPU執行計算任務時，要將計算所需數據從CPU內存通過PCIE總線拷貝到GPU顯存。這種數據拷貝目前無法避免，會增加額外的存取開銷。為盡量減少額外開銷，本文將關聯性較強的子波矩陣和計算反射系數過程合并為一個處理模塊，卷積操作為一個處理模塊，兩個模塊在GPU端運行。將GPU端計算完成后的計算結果拷貝到CPU內存進行后續處理。根據分步驟的計時結果，統計子波、計算疊加道、確定反射系數位置這幾個模塊耗時較少，且不屬于計算密集過程，需要在CPU端單獨處理。并行改進后的算法流程如圖2所示。

圖2 低頻補償GPU并行算法流程

2.3 并行優化策略

根據并行算法的計算過程，本文對算法最后三個步驟設計了不同的并行方案。在構造子波矩陣過程和求解反射系數過程中，首先分析了計算任務之間的獨立性，然后將可并行的任務進行拆解，最后利用GPU大量輕量級的線程對任務進行細粒度的并行實現。兩信號的卷積部分利用卷積定理，先將兩個時域信號通過傅里葉變換轉換到頻域，再將兩個頻域信號點乘，再通過傅里葉逆變換將信號轉換回時間域，并通過cufft庫進行優化。此外，本文計算之前對地震信號數據進行了重新排列，并對數據訪存進行優化。

2.3.1 構造子波矩陣

CUDA的線程(thread)以計算網格(grid)的形式進行組織。一個grid可劃分成多個塊(block)，一個block可包含多個thread，每一個block和thread都有自己的索引編號。當核函數(kernel)執行時，處于激活狀態的thread將對kernel函數編譯后的指令同時調用、執行，根據thread編號實現對不同顯存區的操作。

對于每一個子波矩陣，需要對n個不同的位置賦值子波，彼此之間相互獨立，因此分配n個block對應n個計算位置。block0對應0號采樣點位置，block1對應1號采樣點位置，依此類推。由于線程數一般以2的整數次冪進行分配，因此對每一個block分配256個thread，每一個thread對應一個子波的采樣點。如此分配可以達到子波采樣點數量級別的并行度。

如果子波長度超過256采樣點，每一個block將進行多次計算，直到將子波的所有采樣點被覆蓋為止。對于計算能力強的GPU，可將每個block的線程數增加，比如分配512個thread。這些設置可根據計算任務的不同進行調整。

2.3.2 計算反射系數

用共軛梯度法進行反射系數求解時用到大量矩陣相乘運算。子波矩陣在邏輯上可看作一個二維矩陣，多個子波矩陣數據排列可看作一個三維的數據體。地震記錄可看作一維向量，多道地震記錄排列可以看作一個二維數據體。因此，本文對道集數據進行處理時，先將地震道集數據拷貝到顯存中統一以一維數組的形式進行存儲。每一道數據緊接著上一道數據進行排列，訪問數據時通過首地址加地址偏移量的方式尋址。道集對應的子波矩陣也采用同樣的排列方式存儲在一維數組中。這樣的內存存放有利于適應thread的訪存機制。對于同樣大小的數據量，CUDA的thread連續尋址比跳躍尋址更快。在計算矩陣相乘的時候，相鄰的thread訪問相鄰的內存區域，可大大提高存取效率。計算時，每一個子波矩陣與其對應的地震道進行迭代求解。在一個grid中分配多個block，每一個block負責一個子波矩陣和對應地震道的迭代求解。每一個block內，分配256個thread。

如圖3所示，以block0的計算為例：block0分配256個thread，每一個thread對應一個不同的采樣點，如果采樣點數n大于256則進行多輪計算。地震信號長度為n，反射系數位置個數為m，子波矩陣大小為n行、m列。將地震子波與反射系數卷積過程變為矩陣相乘的形式，r(0)為反射系數位置序列，計算矩陣與向量相乘時，線程編號與采樣點編號對應，一個線程負責m次相乘累加操作。r初始值可選取步驟(3)疊加道反射系數位置的值。計算α、β、f、p變量可根據式(4)～式(8)拆分成n×m維矩陣與長度為m的向量乘法操作、長度為n的向量內積、長度為n的向量加法操作。首先，調用kernel0計算f(0)=s-Mr(0)； 其次，調用kernel1計算α； 再次，調用kernel2更新r和f，調用kernel3再一次計算β； 最后，調用kernel4更新變量p。由于kernel函數包含隱式同步，可以保證上一步向量更新完畢再進行下一步計算，符合該算法的數據依賴關系。經過多次迭代直到f小于誤差限時跳出循環，得到反射系數序列r。在計算完成后，根據反射系數的位置索引，將r的值賦給初始化為零的長度為n的數組中，得到最終的反射系數序列。由于共軛梯度法具有二次終止性，理論上通過n次迭代可得到精確解。實際工程中可設置固定的迭代次數以滿足精度要求。本文最終選定的迭代次數為5。

圖3 計算反射系數

在子波矩陣和地震道數據相乘時，地震道的數據是重復利用的，子波矩陣的每一行需與對應的反射系數相乘。因此，可將初始化反射系數事先存入共享內存進行加速。全局內存對所有的thread可見，共享內存只對同一個block的thread可見。但是共享內存帶寬大約為1.5TB·s-1，遠高于全局內存帶寬190GB·s-1。在子波矩陣與地震數據進行乘法計算時，從共享內存中讀取數據，大大提高了存取效率。

2.3.3 重構寬頻地震記錄

原始版本的寬頻地震記錄重構方法是將反射系數與寬頻子波卷積得到拓頻的地震記錄，假設反射系數與子波采樣點數均為N，對于每一個采樣點位置，需要進行N次相乘再相加的操作，時間復雜度為O(N2)。根據頻率域卷積定理，兩個函數卷積的傅里葉變換是函數傅里葉變換的乘積。

快速傅里葉變換時間復雜度為O(NlbN)，向量點乘時間復雜度為O(N)，改進后的算法需要兩次傅里葉正變換，一次傅里葉逆變換，一次長度為N的向量點乘，總體時間復雜度為3O(NlbN)+O(N)。

可見將時間域卷積計算方式改為頻率域點乘，再進行逆傅里葉變換，可有效降低算法的時間復雜度。在頻率域點乘改進的基礎上，將算法部署到GPU設備上進行，傅里葉變換部分采用cuFFT庫進行。在本次計算任務中，傅里葉正變換為從實向量到復數向量的變換，該變換后的向量為共軛對稱復向量。假設向量長度為J，運用cufftExecR2C句柄只用計算前J/2+1個值。向量點乘時由于共軛對稱復向量點乘計算結果仍然共軛對稱，所以只計算J/2+1長度的向量點乘結果。如此，相對于常規的復數形式的傅里葉變換可以節省約一半的計算量。計算后的向量仍然共軛對稱，運用cufftExec-C2R句柄，將長度為J/2+1的復向量逆變換到時域，結果為長度為J的實向量。計算共軛復向量點乘時，由于每一個采樣點的計算互相獨立，具有良好的并行性。核函數進行計算時，每一個thread負責一個采樣點的乘法計算并且互不影響，如此可達到采樣點數量級別的并行度。

3 應用

用實際地震數據對本文的GPU改進版進行精度和速度測試，計算環境與算法瓶頸分析平臺相同。精度測試中，比較了GPU改進版本與CPU版本對單一道集處理結果的誤差，展示了算法對單一道集的處理效果，以及該算法對疊加剖面成圖效果的影響； 速度測試中，比較了對單一道集和多道集數據兩個版本的計算效率。

3.1 CPU與GPU算法精度比較

測試數據為川西地區實際地震數據，Inline、Xline測線數均為1000，每一個道集包含51道，每一道包含2000個采樣點。圖4a為算法處理前的共反射點道集； 圖4b為CPU版本處理后的共反射點道集； 圖4c為GPU版本算法處理后的共反射點道集。對比圖4a和圖4b可以看出，低頻補償之前同相軸不連續，低頻補償后，強軸反射位置更加突出，一些隱蔽的同相軸得以顯現。比較CPU(圖4b)與GPU版本(圖4c)的計算結果，可見二者基本一致，相對誤差在10-6數量級內，跟浮點數的存儲形式和浮點數計算舍入誤差有關，在單精度浮點數計算誤差范圍內。

圖4 CPU與GPU版本道集處理結果對比

圖5a、圖5b分別為頻率補償前和本文的CPU版本處理后的疊加剖面?？梢?，頻率補償后疊加剖面層位特征更明顯，綠色箭頭處標志層位的細節更加豐富。對比頻率補償前、后的疊加剖面的頻譜可以發現：補償前中心頻率Qc為39.0Hz，絕對帶寬(Q2-Q1)為28.0Hz(圖5c)； 補償后中心頻率為48.5Hz，絕對帶寬為39.0Hz(圖5d)，主頻提升了9.0Hz，有效帶寬拓展了11.0Hz。

圖5 疊加剖面計算結果

以上實例說明本文的GPU改進版本的有效性。

3.2 CPU與GPU版本算法速度比較

首先討論了共軛梯度法求解反射系數過程中，迭代次數對于算法效率及精度的影響。然后比較了在相同迭代次數條件下CPU版本與GPU版本的計算效率。

3.2.1 迭代次數對算法耗時和精度的影響

本文對不同精度的計算結果進行了測試。迭代次數為3時，相對誤差較大不滿足要求； 迭代次數為5次時，前、后兩次計算結果的相對誤差在10-2數量級內； 迭代次數為11時，計算結果相對誤差在10-3數量級； 迭代次數為21時，計算結果相對誤差在10-5數量級。迭代次數越多，算法的精度越高，計算越耗時。但是，迭代次數對CPU版本算法的效率和GPU版本算法的影響并不相同。本文對工區中200個道集的數據進行了不同迭代次數的處理，測試結果如表2所示。

表2 迭代次數對算法性能的影響

可見，迭代次數的增加對CPU版本的影響較大，對GPU版本則影響較小。該結果間接驗證了并行算法的有效性。如果實際工程要求更高精度的計算結果，需要人為增加迭代次數，GPU版本增加的時間成本遠小于CPU版本，為更高精度的頻率補償提供了可能。由于相對誤差在10-2數量級內已經滿足精度要求，后續測試的迭代次數均設為5。

3.2.2 相同迭代次數CPU與GPU版本的計算速度

首先對單一道集進行測試，用C++計時函數計時。由于改進后算法耗時較少，為保持測試結果相對準確，GPU版本計時采用1000次計算的平均值，CPU串行程序版本和GPU并行程序版本計時結果如圖6所示。

圖6 CPU與GPU版本各步驟耗時對比

為測試該算法模塊的總體加速效果，選用該工區多道地震數據進行測試。測試軟件為RevScope 3.5.4，計時工具為軟件自帶計時模塊。在保證結果正確的前提下，計算GPU版本的加速比

(9)

式中：Tg為GPU改進版本總體運行時間；Tc為CPU版本總體運行時間。測試結果如表3所示。

表3 GPU版本算法性能測試結果

單一道集測試結果表明，算法計算時間由原來的183ms下降到24ms，加速比可達7.6倍，改進版本加速效果明顯。步驟(4)+步驟(5)的計算時間由原來的135ms下降到3.2ms； 步驟(6)計算時間由30ms下降到2.8ms。改進后算法瓶頸在于步驟(3)，但是這一模塊的求解算法為快速迭代收縮閾值算法，此迭代算法求解過程中，數據依賴非常強，不適合并行化，因此本文沒有對該模塊并行優化。

實際軟件運行該處理模塊時，工區測試結果加速比較道集測試結果加速比偏低。這是由于除了本文算法部分，測試軟件端還有額外的開銷，比如額外的函數調用。此外在處理規模較大工區數據時，會先將道集數據按批次從硬盤拷貝至內存。這些額外的開銷會增加算法處理時間，導致加速比降低，屬于正?，F象。

4 結論與認識

本文對頻率補償算法進行分析，發現耗時的主要步驟是計算反射系數、子波與反射系數的卷積。

(1)本文對該算法耗時的部分進行重新設計，將求解反射系數部分設計成并行算法，通過CUDA將算法重寫并部署在GPU上進行。將計算地震子波與地層反射系數卷積部分，從時間域計算改為頻率域計算，并部署在GPU上進行。在PC端運算時，GPU版本相對于CPU版本算法加速比可達4倍以上，明顯降低了地震數據處理的時間成本。

(2)本文對比了CPU版本和GPU版本共軛梯度法求解反射系數過程中，迭代次數對計算精度與計算時間的影響。結果表明，迭代次數增加在提高計算精度的同時，也會增加計算時間，并且對CPU版本計算效率的影響較大。本文并行算法可有效降低迭代次數增加帶來的時間成本。

成都晶石石油軟件有限公司提供了GEO-SCOPE地質放大鏡軟件，在此深表感謝！