統(tǒng)一計(jì)算設(shè)備架構(gòu)下的F-X域預(yù)測濾波并行算法

楊先鳳，貴紅軍*，傅春常

（1.西南石油大學(xué)計(jì)算機(jī)科學(xué)學(xué)院，成都 610500；2.西南民族大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，成都 610041）

（*通信作者電子郵箱gui18328655853@163.com）

0 引言

提高地震資料信噪比是進(jìn)行高分辨率地震資料處理的關(guān)鍵，而恰當(dāng)?shù)厝コS機(jī)噪聲是提高地震資料信噪比的重要環(huán)節(jié)。目前在地震資料處理流程中，F(xiàn)-X 域預(yù)測濾波是去除隨機(jī)噪聲的有效方法。該方法由Canals Luis 于1984 年提出，經(jīng)Gulunay［1］、HornBostel［2］等改進(jìn)后逐漸投入生產(chǎn)使用。國內(nèi)學(xué)者國九英等在文獻(xiàn)［3］與文獻(xiàn)［4］中分別提出了F-X 域預(yù)測與二維濾波并用以及F-X 域預(yù)測在三維場景下的使用，使F-X域預(yù)測技術(shù)進(jìn)一步趨向成熟。

近年來，為了實(shí)現(xiàn)對(duì)地下復(fù)雜結(jié)構(gòu)的精準(zhǔn)勘探，地震資料隨機(jī)噪聲去除技術(shù)進(jìn)一步向著高分辨率方向發(fā)展。主要有Oropeza 等［5］提出的多道奇異譜分析，以及Bekara 等［6］把F-X域經(jīng)驗(yàn)?zāi)B(tài)分解應(yīng)用于地震資料隨機(jī)噪聲的去除。以上兩者處理非穩(wěn)態(tài)信號(hào)的能力都強(qiáng)于F-X 域預(yù)測濾波技術(shù)，但計(jì)算復(fù)雜度更高。

隨著現(xiàn)代探測設(shè)備與新興處理技術(shù)的發(fā)展，地震資料處理精度已經(jīng)能滿足對(duì)地下復(fù)雜結(jié)構(gòu)進(jìn)行處理的需要，但是高分辨率技術(shù)往往意味著大量的計(jì)算時(shí)間。近十年來，統(tǒng)一計(jì)算設(shè)備架構(gòu)（Compute Unified Device Architecture，CUDA）在高性能計(jì)算領(lǐng)域發(fā)展迅速，并且在地震資料處理領(lǐng)域中應(yīng)用廣泛［7-9］。在F-X域去噪領(lǐng)域，目前關(guān)于算法提速的研究較少，其中邴萍萍等［10］提出了改進(jìn)的F-X 域經(jīng)驗(yàn)?zāi)B(tài)分解技術(shù)及其分布式并行實(shí)現(xiàn)，但是并行的實(shí)現(xiàn)基于Matlab平臺(tái)，在實(shí)際應(yīng)用中具有的意義有限。本文在簡述了F-X 預(yù)測濾波原理后，基于CUDA 模型對(duì)算法瓶頸進(jìn)行分析，設(shè)計(jì)并實(shí)現(xiàn)了基于CUDA的并行加速方案。

1 F-X域預(yù)測去噪原理

F-X 預(yù)測濾波技術(shù)的基本原理是有效波在F-X 域是連續(xù)的、可預(yù)測的，而隨機(jī)噪聲是不相干的、不可預(yù)測的［11］。即當(dāng)?shù)卣鹦盘?hào)的每一個(gè)線性同相軸的所有頻率分量都具有相同時(shí)差時(shí)，理論上其每一個(gè)頻率分量在空間x方向上是可以預(yù)測的，通過信噪比高的頻段求頻率分量的預(yù)測算子，然后用對(duì)應(yīng)的預(yù)測算子分別對(duì)x方向的復(fù)數(shù)序列進(jìn)行預(yù)測濾波，從而壓制隨機(jī)噪聲。

1.1 預(yù)測算子推導(dǎo)及求取方法

設(shè)震源子波為w(t)，道間時(shí)差為Δt，地震記錄第一道反射記錄為w(t-t1)（t1為第一道記錄到達(dá)時(shí)間），經(jīng)傅氏變換后為W(f)ei2πft1，對(duì)于某頻率f0，沿空間x方向采樣形成的復(fù)數(shù)序列為：

將式（1）轉(zhuǎn)換為Z變換表示形式：

預(yù)測算子即ei2πf0Δt。一般地，對(duì)于剖面上m組不同時(shí)差的線性同相軸，其Z變換形式為：

式中：Q(Z)為關(guān)于Z的m次多項(xiàng)式，Q(0)=1；P(Z)為關(guān)于Z的m-1次多項(xiàng)式。為了進(jìn)一步說明，將式（3）寫為：

當(dāng)n≥m時(shí)，S(Z)Q(Z)中Zn的系數(shù)為0。若將Z視為一個(gè)預(yù)測誤差濾波器的Z變換，用它對(duì)S（Z）的系數(shù)序列進(jìn)行濾波，從輸出的第m個(gè)點(diǎn)開始，其預(yù)測誤差為0。即對(duì)于剖面上m組有不同時(shí)差的線性反射波同相軸，若預(yù)測算子長度L≥m，就可以將這些反射波預(yù)測出來。

然而在實(shí)際地震數(shù)據(jù)集中，并不能確定具有不同時(shí)差的線性同相軸的個(gè)數(shù)以及每組不同時(shí)差的信號(hào)的道間時(shí)差，因此F-X 預(yù)測濾波的預(yù)測算子由最小平方原理求取。設(shè)某一時(shí)間、空間窗口內(nèi)地震記錄數(shù)據(jù)為d(t，x)，變 換 到f?x域 為S(f，x)，對(duì)于某頻率f0，令預(yù)測算子為g(f0，l)，其預(yù)測誤差能量定義如下：

式中：n=1，2，…，N，N為地震道道數(shù)；l=1，2，…，L，L為預(yù)測算子長度。根據(jù)最小平方原理求解式（5）得到使E(f0)最小的g(f0，l)即對(duì)應(yīng)頻率的預(yù)測算子。

1.2 F-X域預(yù)測濾波流程

F-X 域預(yù)測濾波算法處理三維工區(qū)時(shí)，先進(jìn)行傅氏變換將T-X 域數(shù)據(jù)轉(zhuǎn)換到F-X 域，同時(shí)為提高數(shù)據(jù)讀取效率，將數(shù)據(jù)由F-X 排列轉(zhuǎn)變?yōu)閄-F 排列，即將每個(gè)頻率成分的數(shù)據(jù)按照空間位置順序存儲(chǔ)，如圖1（a）所示；之后每次處理一個(gè)頻率成分?jǐn)?shù)據(jù)，處理過程分為兩步，如圖1（b）所示。

1）在Crossline方向取Windowsize×Inline大小窗口；

2）窗口內(nèi)沿著Inline 方向每次取子窗口數(shù)據(jù)進(jìn)行一次濾波。

圖1（b）中算法對(duì)各個(gè)子窗口的計(jì)算是獨(dú)立的，因此對(duì)算法并行優(yōu)化時(shí)將各個(gè)子窗口的計(jì)算在GPU 上并行執(zhí)行，利用GPU多核多線程的特點(diǎn)提升算法效率。

圖1 F-X濾波流程Fig.1 Flowchart of F-X filtering

2 算法性能瓶頸分析

2.1 計(jì)算瓶頸

為分析算法的計(jì)算瓶頸，在本文4.1 節(jié)介紹的CPU 計(jì)算平臺(tái)上，對(duì)采樣點(diǎn)數(shù)量為1 001、大小為250 × 250 的工區(qū)進(jìn)行測試，對(duì)算法的快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）、濾波、逆快速傅里葉變換（Inverse Fast Fourier Transform，IFFT）三部分進(jìn)行耗時(shí)分析，結(jié)果如表1 所示?？梢钥闯鯢FT與濾波部分是算法較為耗時(shí)的部分。算法逐道對(duì)地震數(shù)據(jù)作短時(shí)傅里葉變換（Short-Time Fourier Transform，STFT），處理的時(shí)窗長為150，循環(huán)取時(shí)窗進(jìn)行計(jì)算的過程造成效率瓶頸。濾波部分以子窗口數(shù)據(jù)為處理單位，串行循環(huán)處理是算法執(zhí)行效率的主要瓶頸所在。

表1 整體耗時(shí)分析Tab.1 Overall time consumption analysis

2.2 數(shù)據(jù)冗余讀取

為方便說明，以原算法中參數(shù)為說明對(duì)象。如圖1 所示，每次處理時(shí)，CPU 從硬盤讀取一個(gè)頻率切片的數(shù)據(jù)到內(nèi)存。CPU 每次取20 × 20 鄰域數(shù)據(jù)進(jìn)行濾波，之后沿Inline 方向偏移17 道繼續(xù)處理；Inline 方向一次處理完成后，沿Crossline 方向的偏移為17 道。此處存在兩次數(shù)據(jù)冗余讀取，以一個(gè)子窗口出發(fā)說明。第一處是在Inline方向上，相鄰子窗口存在20 ×3 的數(shù)據(jù)冗余讀??；第二處是在Crossline 方向上，相鄰子窗口存在20 × 3的數(shù)據(jù)冗余讀取。

3 并行優(yōu)化策略

針對(duì)第2 章介紹的算法性能瓶頸，主要從兩方面對(duì)算法進(jìn)行效率優(yōu)化：一方面是基于CUDA 對(duì)FFT 與濾波部分設(shè)計(jì)并行優(yōu)化策略；另一方面是GPU 端基于共享內(nèi)存與CUDA 流技術(shù)優(yōu)化數(shù)據(jù)冗余讀取。

3.1 FFT并行優(yōu)化策略

本文算法對(duì)數(shù)據(jù)進(jìn)行傅里葉變換時(shí)，以數(shù)據(jù)道為單位，每次取道上150 個(gè)數(shù)據(jù)，補(bǔ)零成256 個(gè)點(diǎn)進(jìn)行一次傅里葉變換（如圖2（a）），數(shù)據(jù)體較龐大時(shí)，這部分有大量的循環(huán)計(jì)算造成耗時(shí)，針對(duì)這種情況，利用并行FFT來提高效率。

對(duì)FFT進(jìn)行并行化時(shí)，為發(fā)揮多線程計(jì)算優(yōu)勢，將一維數(shù)據(jù)排列成二維數(shù)據(jù)矩陣，在行方向、列方向各進(jìn)行一次傅里葉變換。具體操作是將一個(gè)時(shí)窗的256 點(diǎn)排列為16 × 16 矩陣，Block 中64 個(gè)線程，16 個(gè)線程為一組完成一個(gè)時(shí)窗的傅里葉變換。先進(jìn)行行方向變換，各個(gè)線程按編號(hào)從顯存讀取對(duì)應(yīng)行的數(shù)據(jù)，此時(shí)的讀取滿足合并內(nèi)存訪問，不足的數(shù)據(jù)自動(dòng)在線程內(nèi)初始化為0，這樣同時(shí)也優(yōu)化了CPU 中補(bǔ)零過程的耗時(shí)，各線程完成處理后將數(shù)據(jù)存入共享內(nèi)存；再進(jìn)行列方向變換，讀取對(duì)應(yīng)列數(shù)據(jù)時(shí)進(jìn)行合并線程訪問，即線程束的線程訪問相鄰數(shù)據(jù)，且讀取的開始位置是對(duì)齊的，則一組線程的讀取會(huì)一次命中。本文數(shù)據(jù)是浮點(diǎn)型復(fù)數(shù)，每個(gè)數(shù)據(jù)占8 B，一行數(shù)據(jù)占128 B，本文顯卡基于開普勒架構(gòu)，支持32 B、64 B、128 B 三種緩存行模式，因此16 線程按相鄰位置讀取時(shí)滿足合并線程訪問，保證了數(shù)據(jù)讀取效率。

圖2 FFT并行化流程Fig.2 FFT Parallel flow diagram

3.2 濾波部分并行優(yōu)化策略

本文算法以20 × 20鄰域數(shù)據(jù)為單位進(jìn)行濾波，為方便說明，Inline、Crossline 方向分別稱為x、y方向。對(duì)于每一個(gè)數(shù)據(jù)點(diǎn)，考慮其與x、y方向上鄰近數(shù)據(jù)的相關(guān)性，線性預(yù)測算子演變?yōu)榫匦晤A(yù)測算子。預(yù)測算子的求?。?2］如下所示：

其中：R(t)、R(t+a)是地震記錄的不同延遲時(shí)刻的自相關(guān)；h(t)為濾波因子。將式（6）表示成矩陣形式為：

其中，A稱為Toeplitz矩陣。

濾波過程分為求Toeplitz 矩陣A、求濾波因子、反褶積濾波三步，對(duì)濾波部分的并行設(shè)計(jì)由這三部分展開。

3.2.1 求Toeplitz矩陣并行優(yōu)化策略

求矩陣A時(shí)，將數(shù)據(jù)每道做自相關(guān)運(yùn)算并將不同時(shí)刻的值相加求和，就得到A的第一行，同時(shí)也不難得到其余n行。本文預(yù)測算子長度為7× 7，表示為h(-3：3，-3：3)，分別對(duì)應(yīng)待預(yù)測值上下左右三鄰域數(shù)據(jù)點(diǎn)權(quán)值。預(yù)測算子長度為49，由式（7），A列數(shù)為49，若以R(20，20)表示x方向上1～20個(gè)點(diǎn)，y方向上1～20 道每道做自相關(guān)的值之和，則矩陣A第一行對(duì)應(yīng) 的 值 為R(20 -3：20+3，20 -3：20+3)，其 中，R(20 -3，20)表示數(shù)據(jù)矩陣沿著x方向左移三位對(duì)應(yīng)鄰域內(nèi)各道做自相關(guān)運(yùn)算的值求和。

進(jìn)行并行優(yōu)化時(shí)，將對(duì)應(yīng)子窗口數(shù)據(jù)拷貝到Block共享內(nèi)存中，64 個(gè)線程一組先進(jìn)行矩陣A、b的計(jì)算，每個(gè)線程計(jì)算A第一行中對(duì)應(yīng)時(shí)刻的值。由于數(shù)據(jù)矩陣在x、y方向移動(dòng)時(shí)，中心14 × 14區(qū)域始終參與計(jì)算，本文先對(duì)這部分?jǐn)?shù)據(jù)進(jìn)行并行求和歸約如圖3，將20 × 20 數(shù)據(jù)處理成7× 7 數(shù)據(jù)后，各線程讀取相應(yīng)值完成A、b的計(jì)算。

圖3 并行歸約流程Fig.3 Flowchart of paralleling and reduction

3.2.2 求預(yù)測算子并行優(yōu)化策略

已知A、b后，采用雅可比迭代根據(jù)式（7）對(duì)預(yù)測算子h進(jìn)行求解。將式（7）每行改寫為hi的表達(dá)式形式，記迭代次數(shù)為k，hik表示第k次迭代hi的值。hik值由下式計(jì)算：

表示通過hik-1計(jì)算hik，不斷迭代使hik趨近解。本文對(duì)雅可比迭代進(jìn)行并行優(yōu)化來提高計(jì)算效率，初始化hi0為0.01，收斂閾值ε為0.000 1，雅可比迭代的并行計(jì)算過程如圖4所示。其中：Flag 是對(duì)所有線程可見的共享變量；hThreadIDk與hThreadIDk-1是線程中存儲(chǔ)的對(duì)應(yīng)下標(biāo)兩次迭代的值。每次更新h時(shí)，不同下標(biāo)的hi互不影響，因此讓一個(gè)線程進(jìn)行式（7）中一個(gè)hi方程的迭代。Flag 初始化為false，若當(dāng)前線程兩次迭代值誤差超過閾值，則置Flag為true，繼續(xù)迭代。

3.2.3 反褶積濾波并行優(yōu)化策略

反褶積濾波是預(yù)測算子與數(shù)據(jù)點(diǎn)逐點(diǎn)卷積的過程。預(yù)測算子為7× 7，數(shù)據(jù)矩陣為20 × 20，遍歷數(shù)據(jù)矩陣每個(gè)值，以其為中心，取上下左右三鄰域內(nèi)7× 7 的數(shù)據(jù)（不足的數(shù)據(jù)為0）與預(yù)測算子卷積，結(jié)果表示當(dāng)前值的濾波結(jié)果。

為了提高讀寫帶寬，共享內(nèi)存一般分成32塊，即bank，一個(gè)bank 大小為4 B。線程對(duì)不同塊以及同塊同位置數(shù)據(jù)的訪問不會(huì)產(chǎn)生沖突，但訪問同塊不同位置時(shí)會(huì)產(chǎn)生沖突，此時(shí)讀取效率低下。對(duì)反褶積部分進(jìn)行并行化時(shí)，數(shù)據(jù)矩陣存儲(chǔ)在共享內(nèi)存中，Block中線程讀取對(duì)應(yīng)區(qū)域數(shù)據(jù)進(jìn)行一個(gè)數(shù)據(jù)的的卷積計(jì)算，此時(shí)20 × 20線程與20 × 20數(shù)據(jù)一一對(duì)應(yīng)，由于線程對(duì)共享內(nèi)存的數(shù)據(jù)訪問重疊性較大，很容易發(fā)生Bank conflict。

為了減少Bank conflict對(duì)效率的影響，記編號(hào)在同一行的線程為一組，線程按行分為20 組。同一組的線程從同一起始位開始讀取，不同行的線程組設(shè)置不同的起始位start，start 的取值為（0：21），此時(shí)20 組線程的讀取起始位剛好在矩陣次對(duì)角線上，如圖5，有16 組的讀取不會(huì)產(chǎn)生Bank conflict，之后組號(hào)和為20 的線程交換起始位，這樣各組線程則完成了一行數(shù)據(jù)的讀取。之后各組上下移動(dòng)三次，各個(gè)線程就完成了所需數(shù)據(jù)的讀取。

圖4 雅可比迭代流程Fig.4 Flowchart of Jacobi iteration

圖5 數(shù)據(jù)矩陣讀取示意圖Fig.5 Schematic diagram of data matrix reading

3.3 冗余讀取優(yōu)化策略

解決數(shù)據(jù)冗余讀取主要使用了CUDA 流。CUDA 模型中，CUDA 流表示GPU 的一個(gè)操作隊(duì)列，在支持異步拷貝引擎的Tesla 顯卡上，數(shù)據(jù)拷貝與計(jì)算可同時(shí)執(zhí)行，某一條流的數(shù)據(jù)拷貝時(shí)間可以被其他流的計(jì)算時(shí)間掩蓋，CUDA 流的應(yīng)用價(jià)值已經(jīng)得到驗(yàn)證［13-15］。將20個(gè)窗口求Toeplitz矩陣、求預(yù)測算子和反褶積濾波的過程視為一個(gè)計(jì)算任務(wù)，由一個(gè)流執(zhí)行一個(gè)任務(wù)的數(shù)據(jù)拷貝與計(jì)算。對(duì)多個(gè)任務(wù)的計(jì)算采用流水線方式如圖6，Stream2的數(shù)據(jù)拷貝與Stream1的數(shù)據(jù)計(jì)算相互重疊，因此可提高整體效率。

圖6 冗余讀取優(yōu)化示意圖Fig.6 Schematic diagram of redundant reading optimization

一個(gè)流中各個(gè)子窗口的計(jì)算相互獨(dú)立，可再細(xì)分到Block中并行處理。為降低Inline 方向數(shù)據(jù)的冗余讀取對(duì)效率的影響，將一個(gè)流中的8 個(gè)子窗口打包成一個(gè)數(shù)據(jù)段，由一個(gè)Block執(zhí)行一個(gè)數(shù)據(jù)段的計(jì)算，利用共享內(nèi)存提高子窗口間重疊數(shù)據(jù)的訪問效率。若工區(qū)一條測線上數(shù)據(jù)道數(shù)為2 000，則一個(gè)流中有300個(gè)數(shù)據(jù)段，即需要申請(qǐng)300個(gè)Block進(jìn)行計(jì)算。

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

4.1 實(shí)驗(yàn)環(huán)境

本文實(shí)驗(yàn)采用的硬件配置：Inter Corei5-8300H，主頻2.30 GHz；內(nèi)存8 GB；NVIDIA Tesla K20c，顯存為4 733 MB。軟件環(huán)境：Windows 10（64 位）操作系統(tǒng)；VS2010 +QT 4.8.7+CUDA9.1 編程環(huán)境。測試數(shù)據(jù)為某工區(qū)疊后數(shù)據(jù)，具體參數(shù)如下：數(shù)據(jù)大小為30.9 GB，線號(hào)范圍為1 120～1 370，共計(jì)88 801道，每道采樣點(diǎn)數(shù)為1 001。

4.2 計(jì)算精度分析

圖7 為1 130 測線在CPU 及GPU 平臺(tái)上的計(jì)算剖面效果對(duì)比，無論從整體趨勢還是細(xì)節(jié)來看，CPU 與GPU 平臺(tái)計(jì)算結(jié)果基本一致。從兩平臺(tái)的計(jì)算結(jié)果誤差圖（c）可以看出，本文并行算法的誤差絕對(duì)值達(dá)到10-9，能滿足實(shí)際工程處理對(duì)于精度的要求。

4.3 效率對(duì)比

為了驗(yàn)證并行算法在性能上的提升情況，選取采樣點(diǎn)均為1 001的不同工區(qū)，分別測試了優(yōu)化冗余讀取前后算法各部分在GPU 與CPU 平臺(tái)上的計(jì)算時(shí)間，結(jié)果如表2 所示。從局部來看，F(xiàn)FT 過程有較多從硬盤拷貝數(shù)據(jù)的過程，這個(gè)過程會(huì)隨著工區(qū)規(guī)模增加而明顯影響程序的性能，因而FFT 部分的加速比增幅略低。求Toeplitz矩陣過程中，各線程訪問共享內(nèi)存時(shí)存在難以優(yōu)化的Bank conflict，這是提高其優(yōu)化效率的主要瓶頸。利用雅可比迭代求預(yù)測算子時(shí)，每次需要對(duì)各線程迭代結(jié)果進(jìn)行同步，之后在單線程中進(jìn)行收斂判斷，這是算法流程的固有瓶頸。

從整體來看，算法效率隨著數(shù)據(jù)規(guī)模的提升而提升?？傮w加速Ⅰ為優(yōu)化冗余讀取前算法整體在CPU/GPU 平臺(tái)上的效率比較，由于CPU 處理部分的開銷，算法效率提升了6.4～8.1 倍；總體加速Ⅱ?yàn)閮?yōu)化冗余讀取后算法的效率對(duì)比情況，可以看出優(yōu)化后并行算法效率是串行算法的8.9～11.9 倍，表明了這部分工作對(duì)本文算法的必要性。

5 結(jié)語

本文在NVIDIA GPU 平臺(tái)上利用CUDA 并行編程語言，對(duì)F-X 域預(yù)測濾波算法進(jìn)行并行優(yōu)化使算法效率達(dá)到7.4～9.1的加速比；并在優(yōu)化CPU/GPU 數(shù)據(jù)拷貝過程后，使算法效率進(jìn)一步達(dá)到8.9～11.9 的加速比。本文算法還存在的瓶頸有兩方面：一是數(shù)據(jù)在FFT 與濾波處理后生成的中間數(shù)據(jù)較大，需要存儲(chǔ)在硬盤上而造成耗時(shí)；二是FFT后的數(shù)據(jù)轉(zhuǎn)存為X-F 順序數(shù)據(jù)時(shí)不滿足合并內(nèi)存訪問，數(shù)據(jù)體較大時(shí)，這部分耗時(shí)也較大。