基于GPU 并行加速的黏聲最小二乘逆時偏移及應用

張 猛

（中國石化勝利油田分公司物探研究院，山東東營 257022）

0 引言

傳統的逆時偏移方法（RTM）具有無傾角限制的優點，能夠處理縱、橫向強變速問題，是目前業界公認的解決高陡復雜構造成像的較好方法［1-2］。然而，逆時偏移在淺層成像中存在較強的低頻噪音，難以實現保幅成像。最小二乘偏移在反演框架下進行偏移成像，能很好地克服上述問題的影響［3-5］。最小二乘偏移技術經歷了射線、單程波和雙程波動理論等發展階段［6-8］。其中，Kirchhoff 積分法最小二乘偏移計算效率高，但是精度較低；基于單程波算子的最小二乘疊前深度偏移雖然解決了多路徑、焦散等問題，且方便在頻率空間域進行Q值補償，但在復雜地質構造情況下，其計算精度依然無法滿足實際生產需要。當前，基于雙程波波動方程的最小二乘逆時偏移［9］（LSRTM）已成為地球物理界的研究熱點。傳統的偏移算法大都是基于聲波介質推導的，沒有考慮地球介質的黏滯性，存在深部儲層成像幅值弱的問題?；陴ぢ暯橘|的偏移成像方法在計算過程中對道集進行校正處理，實現Q值補償，可以達到保幅處理和提高成像分辨率的目的［10］。相對于聲波最小二乘逆時偏移方法，黏聲介質中的最小二乘逆時偏移技術在成像過程中考慮了地層的吸收衰減效應［11-13］，在理論上具有更加明顯的優勢。前人在黏聲最小二乘偏移技術方面做了大量的研究工作。Mulder 等［14］和Hak 等［15］在頻率空間域分析了黏聲介質中地震波偏移成像的多解性，并應用于線性及非線性反演；王雪君等［16］基于點擴散函數的分解與補償，實現了黏聲介質反演成像；李振春等［17］提出了黏聲介質時間域最小平方逆時偏移，在對黏滯性吸收進行校正的同時，還克服了頻率域方法不能應用于大規模模型的限制；還有學者通過low-rank 近似將算法應用到Q補償的RTM和LSRTM 中，并指出了基于Q補償的LSRTM 方法的成像優勢。黏聲最小二乘偏移由于計算能力的限制，早期發展較緩慢。近年來，計算機技術飛速發展，特別是高性能CPU 和GPU 等機器設備的誕生，極大地推動了各種先進地球物理方法技術的實用化進程［18］。Micikevicius［19］研究了基于GPU 的有限差分實現問題，給出了快速有限差分實現算法。Zhang 等［20］實現了基于CPU/GPU 異構平臺的全波形反演，并在陸上數據進行了初步應用。黏聲最小二乘逆時偏移技術在理論上具有明顯優勢，并已成功應用于海上數據，但在陸地數據應用的實例較少，在三維陸地資料的應用幾乎沒有成功案例。

本文基于廣義標準線性固體（GSLS）模型，建立黏聲介質波動方程，制定基于GPU 加速的黏聲介質中最小二乘逆時偏移實現流程，對理論模型進行測試，并將該方法應用于勝利油田某陸上探區三維地震資料處理，以期提升深層儲層成像精度，為巖性油氣藏的勘探開發提供依據。

1 基于黏聲介質的最小二乘逆時偏移理論

1.1 黏聲最小二乘逆時偏移基本原理

最小二乘反演成像是一個線性化反問題，用于估計地下介質的高波數擾動量，目前主要是指速度擾動量。其核心是建立一個線性化的正問題，用數值模擬結果（炮道集）去逼近實測的波場擾動。在L2 范數意義下，當逼近誤差最小時，高波數擾動解就是所要的解。

基于廣義標準線性固體（GSLS）模型的波動方程［21-23］表示為

式中：MR為松弛模量；τεi與τσi為松弛時間，ms；H(t)為單位階躍函數；I為標準線性體的個數；?為梯度算子；?·為散度算子；*為時間上的卷積算子；ρ為密度，kg/m3；f′為震源項。

首先，當只有速度為變化量時的黏聲介質波動方程為

式中：f=；p為波場；v為速度，m/s；f為震源項。

然后，引入上述方程的Born 近似線性化，其波動方程為

式中：ps為擾動波場；p0為背景波場；m(x) 為慢度擾動量。

上式對應一個線性化的正問題，可寫為

式中：dcal為計算得到的地震數據；m是地球物理參數矢量，如速度、密度等；L描述了依賴于m的地震波場正傳播過程。

最后，定義最小二乘偏移的目標函數為

式中：dobs為觀測到的地震數據；E()m為誤差泛函。應用迭代算法求解上式，迭代公式為

式中：μ為迭代步長；k為迭代次數。

迭代算法的關鍵是計算(Δδm)k+1，該項實質上是目標泛函的梯度項。梯度項的計算是正傳波場和逆時反傳波場的相關。

1.2 基于GPU 加速的黏聲最小二乘逆時偏移流程

GPU（圖形處理器）是專為執行復雜的數學和幾何計算而設計的處理器。在工程領域，通常采用CPU 與GPU 等2 種計算設備協同運算的模式對算法進行加速：由CPU 負責數據輸入、輸出以及邏輯判斷，并將這些數據拷貝到GPU 顯存，在GPU 進行核心算法運算，計算完成后，將計算結果從顯存拷貝回CPU 內存，最后通過CPU 輸出。

三維黏聲最小二乘逆時偏移計算量巨大，應用傳統的CPU 設備計算耗時長。分析計算流程可知，其中Born 正演和梯度計算這2 步計算量大?？紤]到GPU 具有超強的計算性能，應用GPU 計算黏聲LSRTM 中的Born 正演模擬和梯度，其他的數據輸入、輸出及簡單計算由CPU 完成。基于GPU 加速的黏聲LSRTM 的實現流程如圖1所示。

圖1 基于GPU 加速的黏聲最小二乘逆時偏移流程Fig.1 Process of least-squares reverse time migration in visco-acoustic medium based on GPU acceleration

2 理論模型測試

理論測試模型采用國際標準的SEG/EAGE 鹽丘模型。該模型是由來自全球50 多個石油公司、科研單位的地球物理學家、地質學家和計算機專家協同合作建立的，目前主要應用于疊前深度偏移成像算法、波動方程正演模擬、速度分析方法、觀測系統設計、計算軟硬件平臺等的模型驗證。本文數據的觀測參數為：總炮數為4 781；50 束線，炮線間距為160 m，炮點間距為80 m；檢波線間距為40 m，道間距為40 m；每束線96 炮，每炮65 道；記錄長度為5 s，時間采樣間隔為8 ms；空間采樣間隔為20 m；表層最小速度為1 500 m/s，鹽丘最大速度為4 450 m/s。

地層Q值參數反演的方法有很多，有野外測量法、經驗公式法、信號估算法和基于層析理論的估計方法等［24］。李慶忠［25］提出的經驗公式所計算的Q值在勝利油田東部探區具有較好的對應關系，并在油田勘探開發中得到了廣泛應用。其表達式為

式中：Qp為品質因子；vp為縱波速度，km/s。

為了驗證本文方法的有效性，首先進行正演模擬，獲取用于偏移算法驗證的理論數據。同時輸入速度模型［圖2（a）］和 用式（7）計算的Q值模型［圖2（b）］進行黏聲波動方程正演模擬，生成模擬炮集記錄，再分別進行聲波最小二乘逆時偏移和黏聲最小二乘逆時偏移計算，得到偏移剖面。結果顯示，由于考慮了地震波傳播過程中的吸收衰減效應，與聲波LSRTM［圖2（c）］相比，黏聲LSRTM［圖2（d）］偏移剖面地層層位更加清晰、斷層刻畫更加精細、準確，中、深層能量得到了有效補償，成像質量明顯改進。從這2 個偏移剖面中紅線處同一位置波形曲線［圖2（e）］可以看出，黏聲LSRTM 通過在偏移過程中對地震記錄進行了Q補償，其波形［圖2（e）中紅線］相對于沒有考慮吸收衰減的LSRTM的波形［圖2（e）中綠線］，主瓣更窄，分辨率更高，達到了Q補償提高分辨率的目的。同時，與對聲波正演資料進行聲波LSRTM 偏移的結果［圖2（e）中黑線］相比，波形擬合程度很高，從另一個角度驗證了本文方法的有效性和正確性。

圖2 理論測試中的速度模型、Q 值模型和最小二乘逆時偏移結果Fig.2 Velocity model，Q model and least-squares reverse time migration results in theoretical test

圖3 為黏聲LSRTM 算法迭代的收斂曲線。分析全程30 次迭代的收斂曲線可以看出，偏移算法在開始迭代時曲線下降較快，說明收斂速度較快，大約在10 次迭代后，曲線變緩，說明迭代殘差趨于穩定，達到了較好的收斂結果。

圖3 黏聲最小二乘逆時偏移收斂曲線Fig.3 Convergence curve of least-squares reverse time migration in visco-acoustic medium

3 實際資料試處理

本次研究選取了勝利油田某探區三維實際資料進行黏聲最小二乘逆時偏移試處理。工區位于東營凹陷東部，處于斷裂帶和中央隆起帶，地下構造復雜，斷裂發育，油氣資源豐富，油氣勘探潛力巨大，油藏類型以構造-巖性油藏為主。

在進行偏移前，首先對炮集數據進行了常規疊前預處理，包括靜校正、能量均衡、去噪、直達波切除、面波剔除等，以保證處理后的炮集數據與模型正演的數據能夠較好匹配。同時，為了獲得準確的速度模型，需要在時間域進行常規的偏移速度分析以獲得初始速度模型，然后應用DIX 公式將時間域速度轉換成深度域速度，最后通過射線層析反演，計算出精確的深度域速度模型。

該區的地層Q值是應用勝利經驗公式估算的。勝利淺層Q值計算經驗公式為

中深層Q值計算經驗公式為

圖4 為不同的Q值經驗公式計算結果對比。勝利經驗公式是對公式（7）的進一步細化和完善，由于該公式是針對勝利探區實際地質情況的經驗總結，更符合勝利探區的實際，由這2 種計算方法得出的Q值模型總體趨勢是一致的［圖5］。

圖4 不同的Q 值經驗公式計算結果對比Fig.4 Comparison of calculation results of empirical formulas with different Q values

圖5 勝利油田某探區速度模型和Q 值模型Fig.5 Velocity model and Q model of an exploration area in Shengli Oilfield

采用速度模型［圖5（a）］和應用勝利經驗公式計算得到Q值模型［圖5（c）］，分別進行常規聲波逆時偏移和黏聲最小二乘逆時偏移，其結果如圖6 所示。應用黏聲最小二乘逆時偏移對實際資料進行處理，共進行了10 次偏移迭代，其成像質量較常規聲波逆時偏移有一定的提升，在中、淺層同相軸的連續性增強，分辨率有所提升，斷面信息更加豐富，成像精度有所提高，如圖6（a）和（b）中紅色圈定部分所示。通過分析深度域水平切片，黏聲最小二乘逆時偏移成像剖面能有效地刻畫洼陷沉積區域砂礫巖邊界，該方法處理的地震數據反映的內幕信息更加豐富，如圖6（c）和（d）中紅色箭頭所示。

為了進一步分析其效果，將2 種方法處理的深度域地震數據轉換到時間域，然后分別求取了瞬時振幅，同時抽取瞬時振幅數據體的切片［圖6（e）—（f）］，分析可知，黏聲最小二乘逆時偏移比常規逆時偏移提取的瞬時振幅屬性更加清晰，分辨率更高，這一結果將有助于后續巖性油氣藏的有效識別。

為了分析計算效率，抽取實際資料中的某一單炮數據，分別應用基于CPU 的程序和基于GPU 并行加速的程序進行了黏聲LSRTM 處理（僅進行1次迭代），其計算效率如表1 所列。對比可見，單塊TaslaK20 X 相對于Xeon E5-2650 中一個核加速比約為78.5。

4 結論

（1）吸收衰減補償可以校正地震波的振幅和相位，對振幅的校正能夠使中深層成像更加清楚，對相位的校正能夠使成像分辨率得到提升，對分辨薄儲層有利，這項技術是保幅處理的關鍵步驟之一。

（2）相對于常規逆時偏移（RTM），黏聲最小二乘逆時偏移在衰減介質中進行黏性補償，不僅能有效壓制淺層噪音，而且能對地下弱照明區域進行照明補償，提升了地震資料的分辨率，實現了真振幅成像，這對于巖性油氣藏的勘探具有重要意義。

（3）基于GPU 加速的黏聲介質中最小二乘逆時偏移實現流程，相對于基于CPU 的算法，效率大幅提升，解決了制約黏聲最小二乘逆時偏移技術在實際生產中大規模實際應用的計算難題。