基于GPU 并行加速的黏聲最小二乘逆時(shí)偏移及應(yīng)用

張 猛

（中國(guó)石化勝利油田分公司物探研究院，山東東營(yíng) 257022）

0 引言

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

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

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

1.1 黏聲最小二乘逆時(shí)偏移基本原理

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

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

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

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

式中：f=；p為波場(chǎng)；v為速度，m/s；f為震源項(xiàng)。

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

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

上式對(duì)應(yīng)一個(gè)線性化的正問(wèn)題，可寫為

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

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

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

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

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

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

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

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

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

2 理論模型測(cè)試

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

地層Q值參數(shù)反演的方法有很多，有野外測(cè)量法、經(jīng)驗(yàn)公式法、信號(hào)估算法和基于層析理論的估計(jì)方法等［24］。李慶忠［25］提出的經(jīng)驗(yàn)公式所計(jì)算的Q值在勝利油田東部探區(qū)具有較好的對(duì)應(yīng)關(guān)系，并在油田勘探開(kāi)發(fā)中得到了廣泛應(yīng)用。其表達(dá)式為

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

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

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

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

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

3 實(shí)際資料試處理

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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