張 猛

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

0 引言

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

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

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

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

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

基于廣義標(biāo)準(zhǔn)線性固體（GSLS）模型的波動方程［21-23］表示為

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

首先，當(dāng)只有速度為變化量時的黏聲介質(zhì)波動方程為

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

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

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

上式對應(yīng)一個線性化的正問題，可寫為

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

最后，定義最小二乘偏移的目標(biāo)函數(shù)為

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

式中：μ為迭代步長；k為迭代次數(shù)。

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

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

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

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

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

2 理論模型測試

理論測試模型采用國際標(biāo)準(zhǔn)的SEG/EAGE 鹽丘模型。該模型是由來自全球50 多個石油公司、科研單位的地球物理學(xué)家、地質(zhì)學(xué)家和計算機專家協(xié)同合作建立的，目前主要應(yīng)用于疊前深度偏移成像算法、波動方程正演模擬、速度分析方法、觀測系統(tǒng)設(shè)計、計算軟硬件平臺等的模型驗證。本文數(shù)據(jù)的觀測參數(shù)為：總炮數(shù)為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值參數(shù)反演的方法有很多，有野外測量法、經(jīng)驗公式法、信號估算法和基于層析理論的估計方法等［24］。李慶忠［25］提出的經(jīng)驗公式所計算的Q值在勝利油田東部探區(qū)具有較好的對應(yīng)關(guān)系，并在油田勘探開發(fā)中得到了廣泛應(yīng)用。其表達式為

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

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

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

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

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

3 實際資料試處理

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

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

該區(qū)的地層Q值是應(yīng)用勝利經(jīng)驗公式估算的。勝利淺層Q值計算經(jīng)驗公式為

中深層Q值計算經(jīng)驗公式為

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

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

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

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

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

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

4 結(jié)論

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

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

（3）基于GPU 加速的黏聲介質(zhì)中最小二乘逆時偏移實現(xiàn)流程，相對于基于CPU 的算法，效率大幅提升，解決了制約黏聲最小二乘逆時偏移技術(shù)在實際生產(chǎn)中大規(guī)模實際應(yīng)用的計算難題。