統一計算設備架構下的F-X域預測濾波并行算法

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

（1.西南石油大學計算機科學學院，成都 610500；2.西南民族大學計算機科學與技術學院，成都 610041）

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

0 引言

提高地震資料信噪比是進行高分辨率地震資料處理的關鍵，而恰當地去除隨機噪聲是提高地震資料信噪比的重要環節。目前在地震資料處理流程中，F-X 域預測濾波是去除隨機噪聲的有效方法。該方法由Canals Luis 于1984 年提出，經Gulunay［1］、HornBostel［2］等改進后逐漸投入生產使用。國內學者國九英等在文獻［3］與文獻［4］中分別提出了F-X 域預測與二維濾波并用以及F-X 域預測在三維場景下的使用，使F-X域預測技術進一步趨向成熟。

近年來，為了實現對地下復雜結構的精準勘探，地震資料隨機噪聲去除技術進一步向著高分辨率方向發展。主要有Oropeza 等［5］提出的多道奇異譜分析，以及Bekara 等［6］把F-X域經驗模態分解應用于地震資料隨機噪聲的去除。以上兩者處理非穩態信號的能力都強于F-X 域預測濾波技術，但計算復雜度更高。

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

1 F-X域預測去噪原理

F-X 預測濾波技術的基本原理是有效波在F-X 域是連續的、可預測的，而隨機噪聲是不相干的、不可預測的［11］。即當地震信號的每一個線性同相軸的所有頻率分量都具有相同時差時，理論上其每一個頻率分量在空間x方向上是可以預測的，通過信噪比高的頻段求頻率分量的預測算子，然后用對應的預測算子分別對x方向的復數序列進行預測濾波，從而壓制隨機噪聲。

1.1 預測算子推導及求取方法

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

將式（1）轉換為Z變換表示形式：

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

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

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

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

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

1.2 F-X域預測濾波流程

F-X 域預測濾波算法處理三維工區時，先進行傅氏變換將T-X 域數據轉換到F-X 域，同時為提高數據讀取效率，將數據由F-X 排列轉變為X-F 排列，即將每個頻率成分的數據按照空間位置順序存儲，如圖1（a）所示；之后每次處理一個頻率成分數據，處理過程分為兩步，如圖1（b）所示。

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

2）窗口內沿著Inline 方向每次取子窗口數據進行一次濾波。

圖1（b）中算法對各個子窗口的計算是獨立的，因此對算法并行優化時將各個子窗口的計算在GPU 上并行執行，利用GPU多核多線程的特點提升算法效率。

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

2 算法性能瓶頸分析

2.1 計算瓶頸

為分析算法的計算瓶頸，在本文4.1 節介紹的CPU 計算平臺上，對采樣點數量為1 001、大小為250 × 250 的工區進行測試，對算法的快速傅里葉變換（Fast Fourier Transform，FFT）、濾波、逆快速傅里葉變換（Inverse Fast Fourier Transform，IFFT）三部分進行耗時分析，結果如表1 所示。可以看出FFT與濾波部分是算法較為耗時的部分。算法逐道對地震數據作短時傅里葉變換（Short-Time Fourier Transform，STFT），處理的時窗長為150，循環取時窗進行計算的過程造成效率瓶頸。濾波部分以子窗口數據為處理單位，串行循環處理是算法執行效率的主要瓶頸所在。

表1 整體耗時分析Tab.1 Overall time consumption analysis

2.2 數據冗余讀取

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

3 并行優化策略

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

3.1 FFT并行優化策略

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

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

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

3.2 濾波部分并行優化策略

本文算法以20 × 20鄰域數據為單位進行濾波，為方便說明，Inline、Crossline 方向分別稱為x、y方向。對于每一個數據點，考慮其與x、y方向上鄰近數據的相關性，線性預測算子演變為矩形預測算子。預測算子的求取［12］如下所示：

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

其中，A稱為Toeplitz矩陣。

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

3.2.1 求Toeplitz矩陣并行優化策略

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

進行并行優化時，將對應子窗口數據拷貝到Block共享內存中，64 個線程一組先進行矩陣A、b的計算，每個線程計算A第一行中對應時刻的值。由于數據矩陣在x、y方向移動時，中心14 × 14區域始終參與計算，本文先對這部分數據進行并行求和歸約如圖3，將20 × 20 數據處理成7× 7 數據后，各線程讀取相應值完成A、b的計算。

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

3.2.2 求預測算子并行優化策略

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

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

3.2.3 反褶積濾波并行優化策略

反褶積濾波是預測算子與數據點逐點卷積的過程。預測算子為7× 7，數據矩陣為20 × 20，遍歷數據矩陣每個值，以其為中心，取上下左右三鄰域內7× 7 的數據（不足的數據為0）與預測算子卷積，結果表示當前值的濾波結果。

為了提高讀寫帶寬，共享內存一般分成32塊，即bank，一個bank 大小為4 B。線程對不同塊以及同塊同位置數據的訪問不會產生沖突，但訪問同塊不同位置時會產生沖突，此時讀取效率低下。對反褶積部分進行并行化時，數據矩陣存儲在共享內存中，Block中線程讀取對應區域數據進行一個數據的的卷積計算，此時20 × 20線程與20 × 20數據一一對應，由于線程對共享內存的數據訪問重疊性較大，很容易發生Bank conflict。

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

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

圖5 數據矩陣讀取示意圖Fig.5 Schematic diagram of data matrix reading

3.3 冗余讀取優化策略

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

圖6 冗余讀取優化示意圖Fig.6 Schematic diagram of redundant reading optimization

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

4 實驗結果分析

4.1 實驗環境

本文實驗采用的硬件配置：Inter Corei5-8300H，主頻2.30 GHz；內存8 GB；NVIDIA Tesla K20c，顯存為4 733 MB。軟件環境：Windows 10（64 位）操作系統；VS2010 +QT 4.8.7+CUDA9.1 編程環境。測試數據為某工區疊后數據，具體參數如下：數據大小為30.9 GB，線號范圍為1 120～1 370，共計88 801道，每道采樣點數為1 001。

4.2 計算精度分析

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

4.3 效率對比

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

從整體來看，算法效率隨著數據規模的提升而提升。總體加速Ⅰ為優化冗余讀取前算法整體在CPU/GPU 平臺上的效率比較，由于CPU 處理部分的開銷，算法效率提升了6.4～8.1 倍；總體加速Ⅱ為優化冗余讀取后算法的效率對比情況，可以看出優化后并行算法效率是串行算法的8.9～11.9 倍，表明了這部分工作對本文算法的必要性。

5 結語

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