張 猛,王華忠,任浩然,馮 波,隋志強(qiáng),王延光

(1.同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與反演成像研究組,上海200092；2.中國石油化工股份有限公司勝利油田分公司物探研究院,山東東營257022；3.浙江大學(xué)地球科學(xué)系,浙江杭州310027)

速度建模是地震偏移成像中的核心問題,高精度的速度建模是地震成像和屬性反演的基礎(chǔ)。全波形反演(Full waveform inversion,FWI)方法利用整個(gè)地震道集的全波形信息進(jìn)行非線性反演,理論上可以反演出寬波數(shù)分布的速度結(jié)構(gòu)。FWI的理論框架于20世紀(jì)80年代提出,Tarantola[1]詳細(xì)給出了FWI的實(shí)現(xiàn)框架,提出了基于廣義最小二乘反演理論的時(shí)間-空間域全波形反演方法。20世紀(jì)90年代,Pratt[2]又將其發(fā)展到頻率-空間域。

FWI雖然在理論上具有嚴(yán)謹(jǐn)?shù)臄?shù)學(xué)物理基礎(chǔ),然而卻一直難以得到大規(guī)模實(shí)際應(yīng)用,其原因是多方面的。首先,實(shí)際地震數(shù)據(jù)處理中難以提供一個(gè)滿足FWI要求的初始模型；缺少低頻和長偏移距數(shù)據(jù),FWI自身也不能給出可靠的初始速度模型。FWI具有強(qiáng)非線性,精度不高的初始速度模型很容易使得FWI陷入局部極小值。其次,實(shí)際地震數(shù)據(jù)的質(zhì)量也不能滿足FWI對數(shù)據(jù)的要求。噪聲和波動方程難以描述的復(fù)雜地震波現(xiàn)象極大地干擾了FWI反演,子波的未知和空變增加了反演解的非唯一性。再次,FWI的巨大計(jì)算量同樣制約著FWI的大規(guī)模實(shí)際應(yīng)用。

近年來,高密度寬方位角大偏移距地震采集方式的出現(xiàn)以及CPU/GPU異構(gòu)計(jì)算平臺的高速發(fā)展,為FWI提供了更好的數(shù)據(jù)準(zhǔn)備以及計(jì)算能力支撐。全波形反演或其變種方法有了在實(shí)際生產(chǎn)中應(yīng)用的可能性,成為勘探地球物理界的熱點(diǎn)研究課題。由于圖形處理器GPU(Graphic Processing Unit)在科學(xué)計(jì)算中相對CPU的計(jì)算效率優(yōu)勢,通過統(tǒng)一計(jì)算設(shè)備架構(gòu)平臺CUDA(Compute Unified Device Architecture)編寫程序可以讓GPU執(zhí)行計(jì)算任務(wù),近年來在地球物理領(lǐng)域得到了廣泛而成功的應(yīng)用。如Micikevicius[3]研究了基于GPU的有限差分實(shí)現(xiàn)問題,給出了快速實(shí)現(xiàn)有限差分的算法。2010年,趙磊等[4]和劉紅偉等[5]提出了基于CPU/GPU平臺的逆時(shí)偏移實(shí)現(xiàn)策略。Sun等[6]討論了TTI介質(zhì)逆時(shí)偏移計(jì)算對GPU所面臨的挑戰(zhàn)。

在數(shù)據(jù)質(zhì)量與計(jì)算能力得到更大保證的情況下,找到FWI在實(shí)際地震資料中的應(yīng)用方案成為了該理論方法是否實(shí)用的關(guān)鍵。近年來,國內(nèi)外學(xué)者分別探討了多種方案試圖使FWI走向?qū)嵱?。?dāng)前,比較系統(tǒng)的解決方案主要有:通過修改泛函格式和多尺度反演策略來降低反演解的非唯一性[7]；通過利用特征波場(反射波)來反演宏觀背景速度[8]；通過施加各種先驗(yàn)信息和約束條件來降低解空間的大小[9]。

我們從全波形反演理論基礎(chǔ)出發(fā),研究并提出了基于CPU/GPU異構(gòu)平臺的時(shí)空域聲波方程全波形反演算法實(shí)現(xiàn)流程。通過理論模型測試證明方法的高精度和高計(jì)算效率；同時(shí)給出實(shí)際地震資料全波形速度反演試處理的初步應(yīng)用效果。在此基礎(chǔ)上進(jìn)一步討論分析FWI對實(shí)際地震資料質(zhì)量的要求,針對陸上地震資料的生產(chǎn)性應(yīng)用提出FWI的實(shí)用化策略。

1 全波形反演基本原理

全波形反演方法利用疊前地震波場的運(yùn)動學(xué)和動力學(xué)信息重建地下速度結(jié)構(gòu),具有揭示復(fù)雜地質(zhì)背景下構(gòu)造與巖性細(xì)節(jié)信息的潛力。根據(jù)實(shí)際需要,全波形反演可以在時(shí)間、頻率、Laplace等域來實(shí)現(xiàn),從計(jì)算方法上來看,全波形反演的實(shí)現(xiàn)方式主要有梯度導(dǎo)引類方法和牛頓類方法(如高斯-牛頓法)。

以二維聲波方程為例,在時(shí)間-空間域,聲波方程表達(dá)為

(1)

利用全波形反演獲取地震波速度,一般將全波形反演問題定義為求解目標(biāo)函數(shù)取得極小值時(shí)對應(yīng)的速度,即定義誤差泛函:

(2)

式中:p(rg,t;rs)obs和p(rg,t;rs)cal分別為野外觀測與數(shù)值模擬的地震記錄；rs和rg分別表示炮點(diǎn)和檢波點(diǎn)坐標(biāo)；t為時(shí)間。全波形反演是通過模擬地震記錄與觀測地震記錄之間的最佳匹配來更新介質(zhì)速度(或其它參數(shù))的。

時(shí)空域波形反演中,通過將剩余波場的反向傳播與正向波場進(jìn)行互相關(guān)來計(jì)算反問題的負(fù)梯度方向,以實(shí)現(xiàn)速度更新過程。

2 基于CPU/GPU異構(gòu)平臺的全波形反演算法實(shí)現(xiàn)

CPU/GPU異構(gòu)并行是指采用CPU與GPU協(xié)同計(jì)算的模式。與CPU集群適合大粗粒度并行不同,GPU計(jì)算核心多對細(xì)粒度并行的問題更有優(yōu)勢。通常情況CPU負(fù)責(zé)數(shù)據(jù)輸入與輸出以及邏輯判斷,并將這些數(shù)據(jù)拷貝到GPU顯存,在GPU進(jìn)行計(jì)算核心計(jì)算；然后再將計(jì)算結(jié)果從顯存拷貝回CPU內(nèi)存,通過CPU輸出結(jié)果。

全波形反演計(jì)算量巨大,應(yīng)用傳統(tǒng)的CPU計(jì)算耗時(shí)長。綜合考慮CPU和GPU計(jì)算特點(diǎn),我們提出了基于CPU/GPU平臺的時(shí)空域聲波方程FWI實(shí)現(xiàn)流程(圖1)。

圖1 基于CPU/GPU平臺的時(shí)空域聲波方程FWI實(shí)現(xiàn)流程

分析圖1所示實(shí)現(xiàn)流程中的計(jì)算步驟,其中波動方程正演模擬和梯度計(jì)算耗時(shí)較長。在每次迭代過程中,求取梯度的部分需要計(jì)算兩次全波場模擬,而基于單炮的策略容易實(shí)現(xiàn)炮的計(jì)算機(jī)并行。首先,由CPU讀取地震子波、速度場和單炮數(shù)據(jù),然后,分別將以上數(shù)據(jù)信息拷貝到GPU顯存中,通過輸入的子波和速度場進(jìn)行地震波正演模擬計(jì)算。GPU具有超強(qiáng)的計(jì)算性能,應(yīng)用GPU計(jì)算全波形反演中的地震波正演模擬和梯度,數(shù)據(jù)的輸入、輸出和其它計(jì)算由CPU完成。計(jì)算結(jié)束后,將單炮梯度拷貝回CPU內(nèi)存,在內(nèi)存中進(jìn)行累加,得到總的梯度,最后生成目標(biāo)梯度函數(shù)。

得到梯度和步長以后,速度更新過程按照(3)式執(zhí)行:

(3)

式中:g為梯度；v為速度；vopt為更新后的速度；αopt為步長。根據(jù)(3)式可以實(shí)現(xiàn)迭代的速度更新,使梯度導(dǎo)引的反演方法得以實(shí)現(xiàn)。當(dāng)然,還需要設(shè)定一個(gè)迭代終止的條件,可以為

(4)

(4)式規(guī)定了一個(gè)迭代收斂的規(guī)則,即第n次更新的量小于第n-1次迭代的一定比例,如ε=0.001,則認(rèn)為迭代收斂,反演可以終止；如果這個(gè)條件不滿足,則將更新結(jié)果作為輸入,再進(jìn)行下一次迭代。

3 理論模型測試與實(shí)際資料試處理

3.1 理論模型測試

理論模型測試采用Marmousi速度模型(圖2a)。其橫向有497個(gè)網(wǎng)格,縱向250個(gè)網(wǎng)格；橫向與縱向網(wǎng)格間距均為12.5m。正演計(jì)算采用PML吸收邊界。震源子波為10Hz的Ricker子波；炮間距為25m,共200炮；第1炮橫坐標(biāo)為625m。檢波器位置固定,每炮437個(gè)檢波器接收；檢波器間距為12.5m；第1炮檢波器橫坐標(biāo)為375m。時(shí)間步長為1ms；記錄長度為4.096s。

用Marmousi模型平滑后的速度模型(圖2b)作為反演的初始模型,用梯度方法對該模型進(jìn)行迭代,不同迭代次數(shù)的結(jié)果如圖3所示。圖3給出了迭代次數(shù)分別為1次、100次、200次和300次反演得到的速度模型,可以看出經(jīng)過300次迭代后,反演結(jié)果的邊界清晰,且速度值收斂到了真實(shí)值附近。

圖4給出了橫坐標(biāo)分別為3100,4200,5700和6800m處的速度反演結(jié)果?？梢婋S著迭代次數(shù)的增加,反演的速度曲線與真實(shí)模型越來越接近,最終反演結(jié)果真實(shí)可靠。

圖2 Marmousi速度模型(a)及其平滑后的模型(b)

圖3 迭代次數(shù)為1次(a)、100次(b)、200次(c)和300次(d)的全波形速度反演結(jié)果

圖4 橫坐標(biāo)為3100m(a),4200m(b),5700m(c)和6800m(d)處的速度反演結(jié)果垂向?qū)Ρ?/p>

進(jìn)一步對基于CPU程序和基于CPU/GPU異構(gòu)并行程序進(jìn)行效率對比測試。針對Marmousi速度模型(點(diǎn)數(shù)737,深度采樣250,炮數(shù)300),基于CPU/GPU平臺異構(gòu)并行計(jì)算(1塊Fermi2050圖形卡)完成一次反演迭代,每次反演相當(dāng)于3次正演的計(jì)算量,300炮迭代1次總共耗時(shí)約28min；在CPU平臺(雙路Xeon X5650,2.67Hz CPU),用MPI實(shí)現(xiàn)并行計(jì)算,采用主從并行模式,300炮迭代1次總共耗時(shí)164min；X5650 CPU 1個(gè)核單獨(dú)進(jìn)行串行計(jì)算時(shí)間為1784min。即單塊Fermi2050卡相對于單顆X5650 CPU的加速比約為164/28=5.8,相對于X5650 CPU其中一個(gè)核的加速比約為1784/28=63.8。

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

為了檢驗(yàn)全波形速度反演解決實(shí)際問題的能力,選取勝利油田某工區(qū)實(shí)際三維地震資料抽取的一條二維線數(shù)據(jù)進(jìn)行試處理。該工區(qū)陡坡帶中西段含油層系多,油藏類型豐富,主要以砂礫巖體油藏為主,油氣勘探潛力巨大。

全波形反演中地震子波的估計(jì)方法是根據(jù)地震直達(dá)波到達(dá)時(shí)與偏移距的關(guān)系,首先推算出淺層速度約1740m/s。根據(jù)實(shí)際地震資料波形,設(shè)定選用Ricker子波。通過實(shí)驗(yàn)對比確定子波參數(shù):主頻20Hz,時(shí)延40ms,振幅10000。

圖5為通過常規(guī)速度分析得到的初始速度模型。圖6為經(jīng)50次迭代后的模擬數(shù)據(jù)(圖6a)與實(shí)際觀測數(shù)據(jù)(圖6b)的對比,可以看出,模擬數(shù)據(jù)與觀測數(shù)據(jù)非常接近,基本滿足全波形反演需要。

圖5 常規(guī)速度分析獲得的初始速度模型

圖6 50次迭代后的模擬炮集(a)和實(shí)際觀測炮集(b)

圖7給出了經(jīng)過50次迭代反演后最終的速度模型。圖8給出了采用圖5所示初始速度模型逆時(shí)偏移結(jié)果(圖8a)和全波形反演最終速度模型逆時(shí)偏移結(jié)果(圖8b)。由圖8可見,在淺層和大斷面處,全波形速度反演速度模型的成像質(zhì)量較初始速度模型的成像質(zhì)量有一定的提升,淺層同相軸的連續(xù)性明顯增強(qiáng),分辨率有所提高；淺層斷層更加清晰,中深層斷面信息更加豐富,成像精度有所改進(jìn)。

但是,由于該二維線并非直接由采集而來,而是從三維數(shù)據(jù)中抽取得到的。雖然應(yīng)用全波形反演速度模型進(jìn)行逆時(shí)偏移成像結(jié)果的質(zhì)量有所改進(jìn),但這還不足以說明全波形速度反演技術(shù)在實(shí)際應(yīng)用中可能達(dá)到的應(yīng)有水平。

圖7 經(jīng)過50次迭代反演后的最終速度模型

圖8 采用圖5所示初始速度模型逆時(shí)偏移結(jié)果(a)和全波形反演速度模型逆時(shí)偏移結(jié)果(b)

4 FWI在陸上地震資料應(yīng)用中的問題分析

二維全波形速度反演要取得較好的結(jié)果,需要有高質(zhì)量數(shù)據(jù)體的支持,具體來說對地震資料有如下要求:測線盡量與構(gòu)造走向垂直；inline方向偏移距分布豐富,覆蓋次數(shù)高,含有大偏移距數(shù)據(jù)；數(shù)據(jù)空間采樣規(guī)則,地表高程變化已校正；資料信噪比較高,壓制了面波與隨機(jī)噪聲；激發(fā)子波能量一致性已校正等。此外,對數(shù)據(jù)進(jìn)行子波零相位化整形也是非常重要的。

除了資料因素外,還必須解決的技術(shù)問題包括:選擇盡可能逼真地描述地震波在地下傳播的正演算子(考慮彈性波、吸收衰減和各向異性等)；選擇凸性更好的目標(biāo)反演以及合適的約束條件；建立一個(gè)較為精確的初始地球物理模型；構(gòu)造適當(dāng)?shù)恼齽t化參數(shù),突出地質(zhì)和偏移成像的先驗(yàn)信息參與度；重建地震數(shù)據(jù)中的低頻信息。

但是,在陸上實(shí)際地震資料采集、處理環(huán)節(jié),無論是獲得高質(zhì)量的原始資料還是解決上述技術(shù)問題,均存在較大的困難。為了解決這一瓶頸問題,提出如下實(shí)現(xiàn)策略和流程:

1) 應(yīng)用低頻大偏移距數(shù)據(jù),通過透射波(折射波或回轉(zhuǎn)波)層析,建立宏觀背景速度模型,獲得低頻信息；

2) 應(yīng)用上述背景速度,通過基于最小二乘的疊前深度偏移獲取反射界面信息；

3) 在背景速度和反射界面信息的控制下,通過反射波FWI或者反射波波形層析,獲取整個(gè)工區(qū)的速度信息。

將經(jīng)典的FWI實(shí)現(xiàn)方法拆分為上述技術(shù)步驟來實(shí)現(xiàn),降低了FWI的非線性性,這將是FWI走向?qū)嵱没闹匾緩?。同時(shí),在具體實(shí)現(xiàn)過程中,必須要充分考慮初始模型的先驗(yàn)信息,將多層次、多尺度的多方面信息綜合[10],以人機(jī)結(jié)合的方式進(jìn)行融合實(shí)現(xiàn)。

5 結(jié)束語

我們研究并實(shí)現(xiàn)了基于CPU/GPU異構(gòu)平臺的全波形速度反演算法。理論數(shù)據(jù)的測試結(jié)果表明,基于Marmousi模型合成數(shù)據(jù)的全波形速度反演可以反演出高分辨率的速度場；通過GPU加速計(jì)算可大幅提高計(jì)算效率。進(jìn)一步對勝利油田某探區(qū)實(shí)際地震資料進(jìn)行試處理,取得了初步的應(yīng)用效果。在此基礎(chǔ)上,分析了全波形反演目前在實(shí)際應(yīng)用中存在的問題,提出了全波形反演實(shí)用化的策略。

為了使全波形速度反演實(shí)用化,需要應(yīng)用不同的波形信息分步驟進(jìn)行速度反演,而不是直接利用全波場信息進(jìn)行反演,這樣可以降低經(jīng)典全波形反演方法的非線性性。同時(shí)充分考慮初始模型的先驗(yàn)信息,在實(shí)現(xiàn)過程中綜合應(yīng)用多層次、多尺度的多方面信息。只有在滿足以上條件時(shí),全波形反演技術(shù)才有可能真正解決復(fù)雜的實(shí)際地質(zhì)問題。

致謝:感謝中國石油化工股份有限公司勝利油田分公司物探研究院的資助。

參 考 文 獻(xiàn)

[1] Tarantola A.Inversion of seismic reflection data in the acoustic approximation[J].Geophysics,1984,49(8):1259-1266

[2] Pratt R G.Seismic waveform inversion in the frequency domain,Part I:theory and verification in a physical scale model[J].Geophysics,1999,64(3):888-901

[3] Micikevicius P.3D finite difference computation on GPUs using CUDA[J].Expanded Abstracts of 78thAnnual Internat SEG Mtg,2008,79-84

[4] 趙磊,王華忠,劉守偉.逆時(shí)深度偏移成像方法及其在CPU/GPU 異構(gòu)平臺上的實(shí)現(xiàn)[J].巖性油氣藏,2010,F06:36-41

Zhao L,Wang H Z,Liu S W.Reverse time migration method and its implementation on CPU/GPU heterogeneous platform[J].Lithologic Reservoirs(in Chinese),2010,F06:36-41

[5] 劉紅偉,李博,劉洪,等.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn)[J].地球物理學(xué)報(bào),2010,53(7):1725-1733

Liu H W,Li B,Liu H,et al.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation[J].Chinese Journal Geophysics(in Chinese),2010,53(7):1725-1733

[6] Sun X,Suh S.Maximizing throughput for high performance TTI-RTM:from CPU-RTM to GPU-RTM[J].Expanded Abstracts of 81stAnnual Internat SEG Mtg,2011,3179-3183

[7] Shin C,Min D J.Waveform inversion using alogarithmic wavefield[J].Geophysics,2006,71(3):R31-R42

[8] Xu S,Wang D,Chen F,et al.Inversion on reflected seismic wave[J].Expanded Abstracts of 82ndAnnual Internat SEG Mtg,2012,1473-1476

[9] Ma Y,Hale D,Gong B,et al.Image-guided sparse-model full waveform inversion[J].Geophysics,2012,77(4):R189-R198

[10] 楊勤勇,胡光輝,王立歆.全波形反演研究現(xiàn)狀及發(fā)展趨勢[J].石油物探,2014,53(1):77-83

Yang Q Y,Hu G H,Wang L X.Research status and development trend of full waveform inversion[J].Geophysical Prospecting for Petroleum,2014,53(1):77-83