基于波動方程重建震源子波的三維全波形反演

梁展源 吳國忱 王玉梅

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②中國石化勝利油田物探研究院， 山東東營 257022)

·偏移成像·

基于波動方程重建震源子波的三維全波形反演

梁展源*①吳國忱①王玉梅②

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②中國石化勝利油田物探研究院， 山東東營 257022)

提出一種基于波動方程重建震源子波的三維全波形反演方法，通過提取疊前炮集近炮檢距的直達波作為波動方程求解的邊界條件，重新求解波動方程，并記錄震源點處的時間序列作為地震子波。該方法利用震源波場重建原理模擬震源爆炸沿地表傳播的逆過程，在逆推過程中提取近炮檢距直達波，避免了折射波、反射波以及潛水波等干擾。該方法只使用模型表層速度，可以有效降低對初始模型的依賴。對三維SEG/EAGE推覆體模型以及實際工區數據的測試結果表明，基于波動方程震源波場重建子波與輸入子波只存在相位差異，經過相位調整后，將重建子波用于三維全波形反演，全波形反演速度接近真實速度，特別是在河道處的全波形反演結果更加可靠，觀測記錄與模擬記錄在波形、相位等方面均有較好的對應，驗證了算法的可行性與適用性。

地震子波 全波形反演 震源波場重建 波動方程 SEG/EAGE推覆體模型

1 引言

全波形反演(Full Waveform Inversion，簡稱FWI)可進行高分辨率地震速度成像，近年來逐漸發展為勘探地球物理研究的熱點。FWI基于波動方程全波場數值模擬，是一個局部尋優的非線性迭代反演問題[1-3]。FWI以觀測記錄和合成記錄的最小二乘函數作為目標函數，利用地震波的走時、振幅、相位等全波形信息反演地下介質參數，具有最高分辨率的特性[4-6]。同時，FWI也存在固有缺陷：一是效率問題，主要體現在反演伴隨著巨大計算量與存儲量；二是局部極小值問題，初始速度不準確、子波不匹配等因素都會使迭代陷入局部極小值。

隨著計算機技術的發展，特別是GPU并行程序的應用，FWI的計算效率明顯提高[7-11]，但是巨大的I/O吞吐問題依然存在。根據波動方程的求解原理，常采用震源波場重建方法解決存儲問題，該方法首先由Gauthier等[12]提出并用于逆時偏移過程，在震源波場正傳過程中存儲邊界波場，存儲的邊界波場的層數與使用的差分階數有關，在檢波點逆時傳播的同時進行源波場重建。該方法數值精度高，但是存儲量正比于有限差分算子的階數。為了進一步減少存儲量，Feng等[13]提出一種改進方案，只存儲一層網格的邊界波場，采用從外至內邊界波場從二階差分逐漸升到高階差分的方法進行震源波場重建。

FWI需要匹配正演模擬數據與實際觀測數據，若兩者的子波不匹配，會嚴重影響反演精度[14，15]。在FWI中，地震子波的處理策略主要分為以下兩大類。一類是不依賴子波的FWI方法，通過引入參考道消除地震子波對FWI的影響。Lee等[16]在頻率域將地震數據與參考道相除，得到歸一化波場用于消除地震子波的影響，但是當地震數據存在噪聲時會明顯降低反演精度[17，18]。Zhou 等[19]提出一種用多道平均作為參考道的策略，通過疊加可以消除噪聲影響，但在參考道中仍存在零值使反演不穩定等問題。Choi 等[20]在時間域提出不依賴子波的FWI，并證實該方法的抗噪性更好，但是仍需提供精確的初始速度。另一類是在迭代過程中更新地震子波的方法，首先由Song等[21]提出。在反演迭代過程中將地震子波作為一個參數，同時反演速度與地震子波，稱為震源信號的迭代反演方法(Iterative Estimation of the Source Signature，簡稱IES)。該方法的缺點在于若初始速度不準確，則每次迭代得到的地震子波也會存在巨大誤差，最終導致反演失敗。

在實際應用中最直接的獲得地震子波的方法是提取直達波作為震源信號。然而，由于實際介質中傳播路徑的影響及透射波、反射波以及散射波等的存在，直達波并不能作為真正的地震子波[22]。本文提出一種基于波動方程重建震源子波的三維FWI方法，通過提取疊前炮集近炮檢距的直達波作為波動方程求解的邊界條件，重新求解波動方程，并記錄震源點處的時間序列作為地震子波。該方法利用震源波場重建原理模擬震源爆炸沿地表傳播的逆過程，在逆推過程中提取近炮檢距直達波，避免了折射波、反射波以及潛水波等干擾。該方法只使用模型表層速度，可以有效降低對初始模型的依賴。將重建的地震子波用于三維模型及實際工區的FWI中，驗證了本文算法求取地震子波的可行性與適用性。

2 方法原理

2.1 FWI方法

在數學上FWI可以定義為一種局部尋優的最優化方法，通過建立最小二乘目標函數尋找目標函數的全局最小值。定義目標函數E，表示觀測數據d與模擬數據的L2范數[2]，即

(1)

式中F(m)為正演算子，其中m為模型參數。

在時間域，正演模擬數據一般通過雙程波波動方程獲得，實際上應該考慮波傳播過程中的衰減、各向異性等因素[23]。然而，由于計算資源的限制，通常采用減少各向異性參數的個數[24]、甚至選用各向同性介質的方法，本文的正演方程為以下常密度聲波方程

(2)

式中：P(x,t)為壓力波場；v(x)為速度參數，對應式(1)中的m，其中x為模型空間某點的坐標；f(x,t)為震源項。

由于FWI的高度非線性，通常采用局部線性最優化方法進行迭代反演。根據牛頓迭代法，每次迭代的參數更新可表示為[25]

(3)

(4)

式中：J為Jacobi矩陣；T為模擬地震數據中的最大時間；λ(x,t)表示殘差反傳波場。

2.2 波動方程源波場重建子波

式(2)為常密度聲波方程，在時間域可以用格林函數與震源表示該式的解

(5)

式中：ξ為震源的空間位置；τ為時間；V為體積分域。

求位移的過程可以看作震源項與格林函數在時間域的褶積。在聲介質假設下，式(2)與式(5)用于模擬獲得觀測數據的過程，如震源激發后，波傳播至檢波器時可接收到直達波(圖1)。

圖1 波傳播示意圖

設直達波為dΩ，則重新構建的常密度聲波方程為

(6)

式(6)為常密度聲波傳播方程的定解問題，與式(2)相比，式(6)中不含震源項。以地表接收到的直達波為方程求解的邊界條件，并假設初始條件為零。此過程模擬震源爆炸沿地表傳播的逆過程以求取地震子波，若正過程為震源項與格林函數在時間域的褶積，則求取子波的過程可表示為位移與格林函數在時間域的反褶積

(7)

式(7)為直達波逆時傳播求取地震子波的積分表示式，由于整個過程是數值求解，不需要求得地震子波的解析解。在重新求解時，取近炮檢距的直達波為方程求解的邊界條件，避免了折射波、反射波以及潛水波等干擾。由于利用直達波為方程求解的邊界條件，故式(6)中的v(x)為地表速度，只要求地表處速度接近真實速度，降低了對初始速度的依賴。

數值求解采用時間2階、空間N階的有限差分方法，通過保存不同層數的邊界值，可用不同的震源波場重建方法求解。文中保存一層邊界的直達波數據，使用Feng等[13]提出的震源波場重建方法，其原理如圖2所示。

圖2 波場重建原理圖

以模型地表處右邊界為例，黑色實心圓為模型邊界，此邊界在吸收邊界內部，d為記錄的直達波邊界值。由于只記錄一層直達波值，因此計算P2用2階差分，計算P3用4階差分，以此類推，…，計算PN/2用N階差分。在N/2層內，用以上方法計算，超出N/2層后，計算PK用N階差分

圖3 FWI流程圖

本文在子波計算與梯度求取時均利用了震源波場重建方法(圖3)，前者利用直達波逆時回推得到地震子波，后者以計算代替存儲減少I/O的消耗。

3 數值試驗

3.1 模型測試

采用三維SEG/EAGE推覆體速度模型(圖4)，其中初始速度(圖4b)由真實速度(圖4a)經過一定平滑后得到。正演模擬所用子波為雷克子波(見圖7a)，正演模擬采用時間2階、空間10階有限差分算法。首先，輸入真實速度體(圖4a)，通過三維正演模擬得到疊前炮集。正演模擬僅為了得到沿地表傳播的直達波，因此并不考慮覆蓋次數以及炮檢距大小。為了避免折射波、潛水波等干擾，通常只考慮近炮檢距直達波。圖5為震源點位于(0，0，0)處的正演模擬記錄。

截取近炮檢距信噪比較高的直達波進行計算，以記錄所有炮點得到的震源子波(圖6)。圖7為重建子波與雷克子波。由圖可見，雷克子波與重建子波的頻譜幾乎一致(圖7b)，但兩者在相位上存在差異(圖7a)。因此對重建子波做一步相位調整，以使重建子波與雷克子波的相位一致，得到相位校正后的子波(圖8)。可見重建子波與輸入的雷克子波幾乎只存在相位差別，經過相位校正后，兩者形態幾乎一致。

圖4 三維SEG/EAGE推覆體模型

以圖左上角自定義三維坐標系為依據，x、y方向網格點數均為401，z方向網格點數為121，x、y、z方向網格尺寸均為10m。震源和檢波器都分布在地表處，共121炮，第一炮位于(0，0，0)處，第一個檢波點也位于(0，0，0)處。x、y方向的炮間距均為400m，x、y方向的檢波點間距均為100m，炮檢距為 1000m。地震記錄的時間采樣間隔為0.6ms，時間采樣點數為1600

圖5 震源點位于(0，0，0)處的正演模擬記錄

圖6 重建的三維地震子波

以重建的三維地震子波作為正演模擬子波進行三維FWI測試，測試模型為三維SEG/EAGE推覆體模型(圖4)。觀測地震數據由校正后的雷克子波(圖8中黑線所示)作為震源子波通過正演模擬獲得，模擬地震數據由校正后的重建三維地震子波(圖8中紅線所示)通過正演模擬獲得。 震源和檢波器都分布在地表，共629炮，第一炮位于(0，0，0)處，第一個檢波點也位于(0，0，0)處。x、y方向的炮間距均為160m，x、y方向的檢波點間距均為10m，炮檢距為1000m。地震記錄的時間采樣間隔為0.6ms，時間采樣點數為1600。數值模擬采用時間2階、空間10階有限差分算法，三個方向的空間采樣網格尺寸均為10m。

圖7 重建子波與雷克子波及其頻譜

圖8 相位校正后的子波黑線為雷克子波，紅線為重建子波

圖9為不同迭代次數的三維SEG/EAGE推覆體模型的FWI結果。由圖9可見，隨著迭代次數的增大，構造細節成像越來越清晰，如明顯恢復了z=0.6km處切片的河道(藍色橢圓內)細節。圖10為實際速度與FWI速度對比，圖11為河道處單道速度對比。由圖可見，FWI速度接近真實速度，特別在河道處(黃色虛線)的FWI結果更加可靠。

圖9 不同迭代次數的三維SEG/EAGE推覆體模型的FWI結果

圖10 真實速度(a)與FWI速度(b)對比

圖11 河道處單道速度對比

4 實際應用

選取中國東部三維M區的疊前地震資料進行試算。該區地表起伏平緩，地震資料信噪比較高，所用的初始速度場由偏移速度分析所得，且初始速度場的表層速度接近真實速度場的表層速度，滿足本文算法的假設。兼顧計算機與三維疊前FWI計算量與數據量的需求，選取較小塊的、觀測系統相對規范的、有井位驗證的三維區塊。所選區塊共882炮數據，炮檢距為 3000m，所獲得的疊前資料做過地表一致性校正等處理。截取近炮檢距處信噪比較高的直達波進行計算，記錄所有炮點得到的震源子波(圖12)，可見擬合得到的最終震源子波接近最小相位子波。

用實際資料重建的三維地震子波(圖12)對該區882炮數據進行三維FWI測試，圖13為迭代17次得到的速度連井剖面， 該剖面為反演三維速度場中抽取的一條過井剖面。由圖可見，反演速度與井曲線吻合較好。圖14為觀測地震記錄與模擬地震記錄對比。由圖可見，觀測記錄(圖14a)與模擬記錄(圖14b)在波形、相位等方面均有較好的對應。

圖12 實際資料重建的三維地震子波

圖13 迭代17次得到的速度連井剖面

圖14 觀測地震記錄(a)與模擬地震記錄(b)對比

5 結束語

本文提出一種基于波動方程重建震源子波的三維FWI方法。應用波動方程的求解理論，將地表接收到的直達波作為波動方程求解的邊界條件，以地表速度作為波動方程的速度參數，通過重新求解波動方程得到地震子波。該方法模擬震源爆炸沿地表傳播的逆過程，利用疊前地震資料近炮檢距直達波信息，避免了折射波、反射波以及潛水波等干擾。并且該方法只假定地表處速度接近真實速度，降低了對初始速度的依賴。

利用三維SEG/EAGE推覆體模型進行測試，對每個炮點都可以求取對應的地震子波，多個炮點擬合后得到的地震子波與初始輸入的雷克子波存在相位差異，但振幅譜無明顯差異，經過相位校正后即可進行三維FWI，并可得到較為準確的速度參數。利用實際資料求得的地震子波接近零相位子波，可直接用于FWI反演，通過與連井曲線對比，證明了本文算法的可行性與適用性。尚需指出，本文的算法研究是在假定地表速度接近真實速度的基礎上開展的，并且只考慮了速度因素，而影響地震波傳播的因素很多，如衰減、各向異性等，需考慮這些因素才能更好地模擬震源爆炸沿地表傳播的逆過程，從而求得更準確的地震子波。因此文中算法只適用于地表起伏平緩的地區，對于地表起伏劇烈的地區不一定適用。

[1] Lailly P.The seismic inverse problem as a sequence of before stack migrations.Conference on Inverse Scattering:Theory and Application.Society for Industrial and Applied Mathematics,Philadelphia,PA,1983,206-220.

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

[3] 王毓瑋,董良國,黃超等.降低彈性波全波形反演強烈非線性的分步反演策略.石油地球物理勘探,2016,51(2):288-294.

Wang Yuwei,Dong Liangguo,Huang Chao et al.A multi-step strategy for mitigating severe nonlinearity in elastic full-waveform inversion.OGP,2016,51(2):288-294.

[4] 李志曄,李振春,劉玉金等.基于波場分離的層析波形反演方法.石油地球物理勘探,2016,51(2):295-300.

Li Zhiye,Li Zhenchun,Liu Yujin et al.Tomographic waveform inversion with wave-field decomposition.OGP,2016,51(2):295-300.

[5] Virieux J,Operto S.An overview of full-waveform inversion in exploration geophysics.Geophysics,2009,74(6):WCC1-WCC26.

[6] 付繼有,李振春,楊國權等.聲介質波動方程反射旅行時反演方法.石油地球物理勘探,2015,50(6):1134-1140.

Fu Jiyou,Li Zhenchun,Yang Guoquan et al.Wave equation reflection travel-time inversion in acoustic media.OGP,2015,50(6):1134-1140.

[7] Maurer H,Greenhalgh S,Latzel S.Frequency and spatial sampling strategies for crosshole seismic waveform spectral inversion experiments.Geophysics,2009,74(6):C79-C89.

[8] Krebs J R,Anderson J E,Hinkley D et al.Fast full-wavefield seismic inversion using encoded sources.Geophysics,2009,74(6):C177-C188.

[9] Habashy T,Abubakar A,Pan G et al.Source-receiver compression scheme for full-waveform seismic inversion.Geophysics,2011,76(4):R95-R108.

[10] Ben-Hadj-Ali H,Operto S,Virieux J.An efficient frequency-domain full waveform inversion method using simultaneous encoded sources.Geophysics,2011,76(4):R109-R124.

[11] Li X,Aravkin A,van Leeuwen T et al.Fast randomized full-waveform inversion with compressive sensing.Geophysics,2012,77(3):A13-A17.

[12] Gauthier O,Virieux J,Tarantola A.Two-dimensional nonlinear inversion of seismic waveforms:Numerical results.Geophysics,1986,51(7):1387-1403.

[13] Feng Bo,Wang Huazhong,Tian Lixin et al.A strategy for source wavefield reconstruction in reverse time migration.SEG Technical Program Expanded Abstracts,2011,20:3164-3168.

[14] 陳樂壽,王光鍔.大地電磁測深法.北京:地質出版社,1984.

[15] 楊培杰,印興耀.地震子波提取方法綜述.石油地球物理勘探,2008,43(1):123-128.

Yang Peijie,Yin Xingyao.Summary of seismic wavelet pick-up.OGP,2008,43(1):123-128.

[16] Lee K H,Kim H J.Source-independent full-waveform inversion of seismic data.Geophysics,2003,68(6):2010-2015.

[17] Xu K,Greenhalgh S A,Wang M Y.Comparison of source-independent methods of elastic waveform inversion.Geophysics,2006,71(6):R91-R100.

[18] Choi Y,Min D J.Source-independent elastic waveform inversion using a logarithmic wavefield.Journal of Applied Geophysics,2012,76:13-22.

[19] Zhou B,Greenhalgh S A.Crosshole seismic inversion with normalized full-waveform amplitude data.Geophysics,2003,68(4):1320-1330.

[20] Choi Y,Alkhalifah T.Source-independent time-domain waveform inversion using convolved wavefields:Application to the encoded multisource waveform inversion.Geophysics,2011,76(5):R125-R134.

[21] Song Z M,Williamson P R,Pratt R G.Frequency-domain acoustic-wave modeling and inversion of crosshole data:Part Ⅱ-Inversion method,synthetic experiments and real-data results.Geophysics,1995,60(3):796-809.

[22] 敖瑞德,董良國,遲本鑫.不依賴子波、基于包絡的FWI初始模型建立方法研究.地球物理學報,2015,58(6):1998-2010.

Ao Ruide,Dong Liangguo,Chi Benxin.Source-independent envelope-based FWI to build an initial model.Chinese Journal of Geophysics,2015,58(6):1998-2010.

[23] Zhou H,Amundsen L,Zhang G.Fundamental issues in full waveform inversion.SEG Technical Program Expanded Abstracts,2012,31.

[24] Warner M,Ratcliffe A,Nangoo T et al.Anisotropic 3D full-waveform inversion.Geophysics,2013,78(2):R59-R80.

[25] Ma Y,Hale D.Quasi-Newton full-waveform inversion with a projected Hessian matrix.Geophysics,2012,77(5):R207-R216.

*山東省青島市黃島區長江西路66號中國石油大學(華東)地球科學與技術學院，266580。Email:lzy09017216@163.com

本文于2016年11月28日收到，最終修改稿于2017年9月6日收到。

本項研究受國家“973”重大專項(2013CB228604)和國家自然科學基金—石油化工基金聯合重點項目(U1562215)聯合資助。

1000-7210(2017)06-1200-08

梁展源，吳國忱，王玉梅.基于波動方程重建震源子波的三維全波形反演.石油地球物理勘探，2017,52(6):1200-1207.

P631

A

10.13810/j.cnki.issn.1000-7210.2017.06.010

(本文編輯：劉勇)

梁展源 博士研究生，1989年生；分別于2013、2016年獲中國石油大學(華東)地球物理學專業學士學位、地質資源與地質工程專業碩士學位； 目前在中國石油大學(華東)地球科學與技術學院攻讀地質資源與地質工程專業博士學位，研究方向為應用地球物理正、反演方法。