基于GPU并行加速的VSP數據逆時偏移

郭雪豹，王建民，王維紅，王云專，柯 璇，劉詩竹

（1.東北石油大學 地球科學學院，黑龍江 大慶 163318； 2.大慶油田有限責任公司 勘探開發(fā)研究院，黑龍江 大慶163712）

0 引言

隨著我國各大油田勘探的不斷深入，探區(qū)內有利儲藏大多已勘探完畢，而剩余的油氣藏多具有油水分布復雜、薄儲層變化快等特點.對于日趨復雜多樣化的油氣藏類型，常規(guī)的偏移方法已無法對地下深層次的復雜構造準確成像.基于波動理論的逆時偏移與同樣基于波動理論的克?；舴蚱坪蛦纬滩ㄆ葡啾?，波動方程不存在任何假設，既不會受到速度橫向劇烈變化的影響，也不存在成像傾角的限制，理論上能夠對多次波、回轉波等準確成像，是目前成像精度最高的方法.

逆時偏移方法最早是在1983年被提出的，Baysal E等［1］提出不同的逆時偏移概念；Whitmore N D、Mc Mechan G A、Loeuenthal D等［2-4］將它應用于疊后偏移.隨著油氣資源的復雜多樣化，常規(guī)的地面地震逆時偏移方法已經無法滿足精細化勘探的需求，尤其對地下深層一些微幅構造及井旁的小斷層等信息.VSP技術是與地面地震觀測方法相對應的，它將檢波器置于井中垂向布置.因此，比地面檢波器更加接近目的層，減少地下淺層對地震信號的干擾，直接接收來自目的層的反射，同樣還能夠接收到一些地面檢波器所接收不到的陡傾角信息，具有比地面地震記錄更高的分辨率.郭建［5］分析VSP技術的應用現狀.朱金明等［6］通過按時間逆的次序求解雙程無反射波動方程，實現VSP逆時偏移.孫文博等［7］采用偽譜法實現逆時偏移方法，分別對分離波場偏移處理，得到較好的效果.Alok K S等［8］采用包含多次波在內的全波場成像，提高照明范圍.Sun Wenbo等［9］在角度域實現VSP逆時偏移.Mark E W等［10］利用walkaway VSP資料在鹽丘側翼上取得較好的成像效果.VSP資料在對地下深層構造成像方面具有優(yōu)勢，并且其內部豐富的波場信息也有助于更高精度成像.

采用VSP數據進行逆時偏移，利用VSP資料的高分辨率、高信噪比的特點，實現地下微幅構造的準確成像.根據設計的VSP觀測系統，筆者正演VSP地震記錄，利用VSP資料的豐富波場信息進行逆時偏移，以彌補常規(guī)地面地震在深層構造及井旁小構造方面的成像不足；逆時偏移算法中波場計算部分由GPU執(zhí)行，其余邏輯判斷部分由CPU執(zhí)行，充分挖掘兩者的計算特性來獲得最高加速比.算法經GPU加速后，計算效率得到了顯著提升.

1 VSP觀測系統

VSP觀測系統是一種由地表震源激發(fā)，井中檢波器接收的觀測方式（見圖1）.首先，與傳統地面檢波器接收到的反射信號路徑相比，VSP檢波器接收的反射波路徑明顯更短，避免了反射波回到地表時所要經過的淺層部分，減少地表低速帶對反射波信息中高頻成分的吸收，振幅畸變相對更小.其次，由于檢波器布置在地下，與地表的檢波器相比更加貼近目的層，因此能夠直接接收到來自目的層的反射信息，有助于研究來自目的層的單一反射信息；在地下存在陡傾角構造時，地表檢波器范圍有限，可能損失部分反射信息，而VSP的檢波器則能接收到.根據波場信息，VSP波場較地面地震波場更為豐富，不僅存在上行波（圖中以上行一次反射波為例），還有下行波（圖中以下行直達波為例），將有利于逆時偏移方法對VSP波場信息的充分利用.因此，與地面地震記錄相比，VSP數據的分辨率和信噪比更高，將給地下深層成像中帶來明顯優(yōu)勢.

2 逆時偏移原理

逆時偏移過程可以分為震源波場正傳和檢波點處波場反傳.首先進行震源波場正傳，在正傳過程中為削弱用計算機有限存儲模擬地下無限介質的波場傳播所帶來的邊界反射，需要添加邊界條件，削弱邊界反射對波場的影響，并且保存每一時刻的波場信息；然后進行檢波點處波場反傳，與對應時刻的震源波場應用成像條件得到逆時偏移結果，對逆時偏移結果去噪處理得到最后的成像剖面.

二維聲波方程：

式中：p為波場值；v為介質速度.

通過有限差分法對波場進行數值模擬，時間二階、空間2L階精度的有限差分格式［12］為

式中：Δx、Δz分別為x、z方向的網格大??；Δt為時間步長；i、j分別為x、z方向網格點的位置；n為時間點；al為有限差分系數.

3 GPU加速技術

盡管逆時偏移方法成像精度極高，但需解決計算量大和存儲需求大的問題.GPU即圖形處理單元，經NVIDIA公司發(fā)布CUDA后，使它可以被開發(fā)人員編寫指令程序，便于讓GPU執(zhí)行CPU所分布的計算任務.GPU具有比CPU更多的計算核心，更易處理大量的并行計算任務，對于逆時偏移過程中的大量波場值計算更加具有優(yōu)勢［11］，因此可以顯著的提高算法的計算效率.

在逆時偏移過程中，計算量主要體現在模擬波場傳播計算中.對于一個二維的剖面，將它按縱橫向進行離散化，得到用網格點表示的剖面.在偏移過程中，不斷利用式（2）計算每一個網格點上不同時刻的波場值.對于同一時刻的波場剖面，每個網格點上的波場值計算并無先后順序，即同一時刻的所有網格點的計算是并行的.首先在GPU上開辟同一時刻所需要的線程數，然后將每一個網格點上的波場值計算放到對應的線程中，GPU上的計算核心多（計算核心即流多處理器，其本身包含標量流處理器，每個線程都是被發(fā)射到一個標量流處理器上執(zhí)行），因此比CPU更加適合并行計算［13-15］.

由式（2）可知，在計算一點的波場值時需要多個點的波場值參與運算，因此需要在GPU上開辟部分存儲以存放波場值.在GPU的存儲器中，全局存儲器是存儲空間最大的，但訪問速度遠不如其內部的共享存儲器.由于共享存儲器存儲有限，首先將計算的數據先讀入全局存儲器；然后在每個線程運算時，將所需要的數據再從全局存儲器復制到共享存儲器中，在計算時可以直接訪問共享存儲器，有效加快計算式的訪問速度，從而進一步提高計算效率.GPU加速計算流程（見圖2）.

4 模型測試

4.1 鹽丘模型

采用鹽丘模型測試算法（見圖3），網格大小及時間步長根據穩(wěn)定性條件選取.模型橫、縱向網格點數分別為800、600，橫、縱向網格大小分別為5、5m.采用地面放炮、井中接收的觀測方式，震源從地表左端50m處開始，向右每隔50m放一炮，共80炮；在分別距模型左端0、4 000m處布置2口井，井中每隔5m布置一個檢波器，每口井放置600個；震源處雷克子波頻率為40Hz，時間步長為0.5ms，時間采樣點為7 500.

在數值模擬過程中，采用保存部分波場信息的完全匹配層邊界條件［16-19］，吸收由計算機有限存儲所帶來的邊界反射衰減，成像條件為互相關成像條件［20］，對逆時偏移結果采用拉普拉斯去噪［21-23］處理.鹽丘模型運算采用GPU加速節(jié)點運行，硬件配置：CPU：E5-2630（2.3GHz／15M／6核心）＊2；內存：64 GB；系統硬盤：1TB；GPU：Nvidia Tesla K10＊3（3 072核心／4.57T單精度／0.19T雙精度／8G顯存）.

VSP逆時偏移的成像剖面見圖4.由圖3和圖4可知，地下構造基本成像清晰，淺層部分成像略顯不足，井旁地層成像準確且清晰，鹽丘輪廓清楚，證實VSP逆時偏移對井旁高分辨率成像的優(yōu)勢，能夠有效利用來自地下深層的反射信息.

4.2 Marmous模型

利用Marmous模型進行測試，網格大小及時間采樣由模型發(fā)布方參數決定.模型縱、橫向網格點數分別為750、993，縱、橫向網格大小分別為6.25、4.00m（見圖5）；采用地面放炮、井中接收的觀測方式，震源從地表左端62.50m處開始，向右每隔62.50m放一炮，共99炮；在分別位于模型的62.50、3 125.00、6 187.50m處3口井，井中每隔4.00m布置一個檢波器，每口井放置750個；震源處雷克子波頻率為40Hz，時間步長為0.4ms，時間采樣點為10 000.

VSP逆時偏移的成像剖面見圖6，地面地震的偏移結果見圖7.由圖5可知，Marmous模型地下地層較多，地勢起伏大，給成像帶來一定難度.由圖7可知，地面地震逆時偏移可以較清晰地對它進行成像，尤其在淺層部分，地層清晰連續(xù)，但在2 000m以下的部位開始模糊.在圖6中2 000m以下的部位，構造形態(tài)基本正確，成像清晰，突出VSP逆時偏移對于深層構造的成像優(yōu)勢.在計算過程中，采用CPU偏移一炮需要2 844s，在相同參數下采用GPU偏移一炮僅需36 s，有效提高偏移的速度70倍以上.

5 結論

（1）逆時偏移方法能夠充分利用VSP數據的豐富波場信息及高分辨率、高信噪比的優(yōu)勢，在井旁及深層構造成像方面具有比地面地震逆時偏移更好的成像效果.但由于對淺層反射信號接收不足，在淺層部分成像效果不如地面地震逆時偏移成像清晰.

（2）逆時偏移算法對計算機的計算能力要求較高，充分利用GPU上的共享存儲器，使其達到顯著的加速比，增強算法的實用性.

（3）二維VSP逆時偏移在一定程度上彌補地面地震逆時偏移對深層成像不足的問題，但隨著勘探目標的日趨復雜，二維VSP逆時偏移也難以達到復雜構造精細勘探的需求，開展三維VSP逆時偏移技術是未來的研究方向.

［1］Baysal E，Kosloff D D，Sherwood W C.Reverse time migration［J］.Geophysics，1983，48（11）：1514-1524.

［2］Whitmore N D.Iterative depth migration by backward time propagation［C］∥SEG Technical Program Expanded Abstracts，1983：382-385.

［3］McMechan G A.Migration by extrapolation of time-dependent boundary values［J］.Geophysical Prospecting，1983，31（3）：413-420.

［4］Loewenthal D，Mulfti I R.Reversed time migration in spatial frequency domain［J］.Geophysics，1983，48（5）：627-635.

［5］郭建.VSP技術應用現狀及發(fā)展趨勢［J］.勘探地球物理進展，2004，27（1）：1-8.Guo Jian.The application status and development trends of VSP technology［J］.Progress in Exploration Geophysics，2004，27（1）：1-8.

［6］朱金明，顏俊華.VSP逆時偏移［J］.石油地球物理勘探，1991，26（5）：564-570.Zhu Jinming，Yan Junhua.VSP reverse-time migration［J］.Oil Geophysical Prospecting，1991，26（5）：564-570.

［7］孫文博，孫贊東.基于偽譜法的 VSP逆時偏移及其應用研究［J］.地球物理學報，2010，53（9）：2196-2203.Sun Wenbo，Sun Zandong.VSP reverse time migration based on the pseudo-spectral method and its applications［J］.Chinese J.Geophys，2010，53（9）：2196-2203.

［8］Alok K S，Xander S，Eric V.VSP imaging using all multiples：Full wavefield migration approach［C］∥SEG Technical Program Expanded Abstracts，2012：1-6.

［9］Sun Wenbo，Sun Zandong.Angle-domain VSP reverse-time migration and its application［C］∥SEG Technical Program Expanded Abstracts，2010：4339-4343.

［10］Mark E W，Lu Rongrong，Xander C，et al.A novel application of time-reversed acoustic：Salt-dome flank imaging using walkaway VSP survey［J］.Geophysics，2006，71（2）：A7-A11.

［11］劉紅偉，李博，劉洪.地震疊前逆時偏移高階有限差分算法及 GPU實現［J］.地球物理學報，2010，53（7）：1725-1733.Liu Hongwei，Li Bo，Liu Hong.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation［J］.Chinese J.Geophys，2010，53（7）：1725-1733.

［12］劉洋，李承楚，牟永光.任意偶數階精度的有限差分數值模擬方法［J］.石油地球物理勘探，1998，33（1）：1-10.Liu Yang，Li Chengchu，Mu Yongguo.Finite-difference numerical modeling of any even-order accuracy［J］.Oil Geophysical Prospecting，1998，33（1）：1-10.

［13］李博，劉紅偉，劉國峰，等.地震疊前逆時偏移算法的CPU／GPU實施對策［J］.地球物理學報，2010，53（12）：2938-2943.Li Bo，Liu Hongwei，Liu Guofeng，et al.Computational strategy of seismic pre-stack reverse time migration on CPU／GPU［J］.Chinese J.Geophys，2010，53（12）：2938-2943.

［14］石穎，陸加敏，柯璇，等.基于 GPU 并行加速的疊前逆時偏移方法研究［J］.東北石油大學學報，2012，36（4）：111-115.Shi Ying，Lu Jiamin，Ke Xuan，et al.Prestack reverse time migration based on GPU parallel accelerating algorithm［J］.Journal of Northeast Petroleum University，2012，36（4）：111-115.

［15］柯璇，石穎，劉詩竹.地震資料逆時偏移中的圖像處理器加速算法［J］.計算機系統應用，2013，22（11）：115-118.Ke Xuan，Shi Ying，Liu Shizhu.GPU acceleration algorithm of the reverse time migration with seismic data［J］.Computer Systems＆ Applications，2013，22（11）：115-118.

［16］郭念民，吳國忱.基于PML邊界的變網格高階有限差分聲波方程逆時偏移［J］.石油地球物理勘探，2012，47（2）：256-265.Guo Nianmin，Wu Guochen.Highorder finite difference method in reverse-time migration with variable grids based on PML boundary condition［J］.Oil Geophysical Prospecting，2012，47（2）：256-265.

［17］王守東.聲波方程完全匹配層吸收邊界［J］.石油地球物理勘探，2003，38（1）：31-34.Wang Shoudong.Absorbing boundary condition for acoustic wave equation by perfectly matched layer［J］.Oil Geophysical Prospecting，2003，38（1）：31-34.

［18］王永剛，邢文軍，謝萬學，等.完全匹配層吸收邊界條件的研究［J］.中國石油大學學報：自然科學版，2007，31（1）：19-24.Wang Yonggang，Xing Wenjun，Xie Wanxue，et al.Study of absorbing boundary condition by perfectly matched layer［J］.Journal of China University of Petroleum：Natural Science Edition，2007，31（1）：19-24.

［19］胡昊，劉伊克，常旭，等.逆時偏移計算中的邊界處理分析及應用［J］.地球物理學報，2013，56（6）：2033-2042.Hu Hao，Liu Yike，Chang Xu，et al.Analysis and application of boundary treatment for the computation of reverse-time migration［J］.Chinese J.Geophys，2013，56（6）：2033-2042.

［20］薛東川，張云鵬，朱振宇，等.幾種疊前逆時偏移成像條件的比較［C］.中國地球物理，2012：581.Xue Dongchuan，Zhang Yunpeng，Zhu Zhenyu，et al.Comparison of several imaging condition for prestack reverse-time migration［C］.The Chinese Geophysical，2012：581.

［21］劉紅偉，劉洪，鄒振，等.地震疊前逆時偏移中的去噪與存儲［J］.地球物理學報，2010，53（9）：2171-2180.Liu Hongwei，Liu Hong，Zou Zhen，et al.The problems of denoise storage in seismic reverse time migration［J］.Chinese J.Geophys.，2010，53（9）：2171-2180.

［22］Zhang Y，Sun J.Practical issues in reverse time migration：true amplitude gathers，noise removal and harmonic source encoding［J］.First Break，2009，26（1）：29-35.

［23］陳可洋.基于拉普拉斯算子的疊前逆時噪聲壓制方法［J］.巖性油氣藏，2011，23（5）：87-95.Chen Keyang.Pre-stack reverse-time noise suppressing method based on Laplacian operator［J］.Lithologic Reservoirs，2011，23（5）：87-95.