基于測井約束的地震全波形反演方法

杜澤源 吳國忱 王玉梅

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266580； ③中國石化勝利油田物探研究院，山東東營 257022)

·偏移成像·

基于測井約束的地震全波形反演方法

杜澤源*①吳國忱①②王玉梅③

(①中國石油大學(華東)地球科學與技術學院，山東青島 266580； ②海洋國家實驗室海洋礦產資源評價與探測技術功能實驗室，山東青島 266580； ③中國石化勝利油田物探研究院，山東東營 257022)

在實際應用中全波形反演依賴于初始速度模型與低頻信息，若初始模型精確度低常會引起周波跳躍現象。測井資料具有垂向分辨率高的特點，含有豐富的超低頻與高頻成分，測井信息的引入可以在一定程度上彌補地震頻帶缺失、反射數據精度低、信噪比低的不足，故從兩個方面將先驗信息引入全波形反演。一是充分利用可靠的測井數據低頻信息建立較精確的初始速度模型，以彌補地震數據缺少低頻信息的問題；二是在全波形反演目標泛函中引入先驗模型約束項，利用測井數據約束反演過程，驅使反演沿給定的方向進行，以提高反演過程的穩定性，從而獲得較高精度的速度模型。通過EAGE推覆體模型對所提方法進行試算，相對不加約束全波形反演結果得到改善；實際地震資料應用結果表明，使用測井約束的反演結果能更清楚地展示砂巖層，具有較高的分辨率，反演效果更好。

全波形反演 低頻信息 測井約束 初始模型 目標泛函

1 引言

全波形反演是利用觀測地震數據反演地球介質地球物理參數的方法，其關鍵是基于波動方程有限差分的正演模擬技術，由于此方法基于雙程波動方程的全解[1]，因此稱為全波形反演。Tarantola[2]從波動方程出發，通過波場殘差逆時傳播計算梯度方向，提出基于最小二乘理論的時間域全波形反演理念，奠定了全波形反演的基礎。全波形反演最早在時間域實現，Petrick等[3]給出了較為完整的全波形反演定義，隨后，Pratt等[4,5]將其發展到頻率域，Shin等[6]將其推廣到Laplace域。全波形反演技術已被證明是建立高精度速度模型的有效方法[7,8]，得到的速度模型可用于改善深度偏移的質量，進而提供非常有價值的成像結果。全波形反演以減小模擬波場與實際觀測波場數據之間的殘差為目標，通過最小化目標泛函反演獲得高精度速度模型[2,9]。在迭代求解過程中，通過在模型空間尋找梯度和局部下降方向使數據空間的最小二乘目標泛函最小，從而獲得每次迭代的速度模型更新量[10]，并將更新速度模型作為下一次迭代的輸入模型。同時，全波形反演效果受諸多因素影響，可靠的大炮檢距數據、低頻信息、初始速度模型等是獲得較好反演結果的基本保障[11]。

在實際應用中全波形反演仍面臨諸多問題[12]，其中一個重要問題是全波形反演依賴于初始速度模型與低頻信息[13]。若初始模型精度低常會引起周波跳躍現象，Weibull等[14]認為全波形反演初始速度模型應該盡量接近真實模型以避免周期跳躍，建議采用波動方程代替射線層析方法獲得較精確的光滑初始模型。敖瑞德等[15]提出了基于包絡的全波形反演初始模型構建方法，可以得到較高精度的初始模型； Gong等[16]提出了時間域逆時偏移的三維全波形反演算法，可以減少反演結果對初始模型的依賴。然而，為了避免對初始模型的依賴，需要盡可能采集含有豐富低頻分量的地震數據，并且要求信噪比和保真度足夠高[17]，這在實際地震采集中很難達到，因此構建較準確的初始速度模型很有必要。段超然等[18]聯合波場重構反演和全波形反演以確保全波形反演在缺少低頻信息時仍能穩定收斂；王毓瑋等[19]提出了一種降低非線性程度的分步反演策略，在重建背景模型的基礎上恢復模型的中、高波數成分。

在實際勘探過程中可獲得豐富的先驗信息，如聲波測井、地質信息等，傳統全波形反演一般不考慮先驗信息。當觀測地震資料缺少大炮檢距數據和低頻信息時，在大多數情況下反演結果不理想。因此，可以在全波形反演中利用上述先驗信息。測井資料具有垂向分辨率高的特點，含有豐富的超低頻與高頻成分，測井信息的引入可以在一定程度上彌補地震頻帶缺失、反射數據精度低、信噪比低的不足，故可將測井信息引入初始速度模型建立。全波形反演需要一個平滑的初始模型以避免在迭代初始就陷入局部極小，而測井數據的分辨率較高，具有高精度的速度及可靠的低頻信息，利用測井數據的低頻信息對初始模型進行約束，可以得到較精確的初始速度模型。另一方面，可以從目標泛函約束的角度引入先驗模型信息，但該過程的操作難度較大，常用于AVO反演。吳國忱等[20]在AVO反演過程中加入尺度約束項，以提高反演的穩定性及可靠性；馮國峰等[21]針對二維波動方程，將大范圍收斂的同倫方法引入算子參數識別過程，并利用測井資料進行約束反演；Hu等[22]建議在倍增的正則化項中使用先驗模型；王守東[23]針對井間地震資料走時層析成像分辨率低、穩定性差的問題，提出了一種井間地震資料成像新方法——全變差正則化波形反演方法(TVWI)；Han等[24]、傅紅筍等[25]研究了測井約束下地震波形反演的同倫攝動求解方法；Asnaashari等[26]提出通過測井、地質解釋建立先驗模型，并使用正則化項和先驗模型項擴展反演目標函數，將先驗模型用于反演優化過程，約束最小化目標泛函；成景旺等[27]使用基于柯西分布的頻率域目標函數，將低頻反演結果作為高頻反演的初始模型以減少解的非唯一性。模型測試結果表明，上述方法可獲得較為可靠的反演速度模型，但很少投入實際應用。

本文從兩個方面將先驗信息引入全波形反演。一是充分利用可靠的測井數據低頻信息建立較精確的初始速度模型，以彌補地震數據缺少低頻信息的不足；二是在全波形反演目標泛函中引入先驗模型約束項，利用測井數據約束反演過程，驅使反演沿給定的方向進行，以提高反演過程的穩定性，從而獲得較高精度的速度模型。EAGE推覆體模型數據測試及實際地震資料應用結果均證明了文中方法的可行性及有效性。

2 理論方法

2.1 測井約束低頻模型構建

鑒于全波形反演對初始模型和低頻信息的依賴性，而地震資料由于采集條件的限制，缺少必要的低頻數據。通過測井數據約束常規初始速度模型，在常規背景速度場中加入可靠的測井數據低頻信息，構建更加精確的初始速度模型，從而避免周波跳躍現象，提高反演精度。

2.1.1 測井資料預處理

測井曲線是地層物理性質隨深度變化的記錄，包含了高分辨率的地層信息，可以通過測井縱波速度數據提取可靠的低頻速度信息。首先要對測井數據進行相應的分析與整理，根據工區的地質背景信息，分析所選測井縱波速度數據，去除奇異值，保證測井數據準確、可靠，然后再進行井震匹配標定，得到準確的時深關系。

實際測井數據一般是深度域的，要想獲得測井數據的低頻信息，需要把井震標定后的測井曲線進行深時轉換，由深度域轉到時間域。然后再分析測井數據的頻譜特征，并與地震數據的頻譜進行對比，綜合考慮地質概況，通過適當的濾波方法對測井數據進行濾波處理，得到低頻測井曲線。若預先獲得的常規初始速度模型是時間域的，則可以在時間域將測井曲線與初始速度模型進行匹配，通過低頻測井數據約束得到新的初始速度模型，然后將模型轉換到深度域；若常規初始速度模型是深度域的，則需要把時間域的低頻測井數據曲線經過時深轉換變換到深度域，然后對常規初始速度模型進行深度域約束，得到可直接使用的初始速度模型。圖1為模型數據實際測井速度曲線與低頻測井曲線。

2.1.2 加權約束模型構建

通過測井數據空間加權約束的方法建立新的初始速度模型。首先，根據工區概況，確定需要參與模型構建的井數及井位信息。將常規方法獲得的初始速度模型看成一個三維矩陣，記為Vinitial(x,y,z)，其中x、y、z為模型空間點位坐標。將低頻測井數據看成縱向矢量，記為Vwell(xk,yk,z)，其中k=1、2、…、N，N為井的總數，xk、yk為井的水平坐標，圖2為測井約束控制示意圖。通過測井數據空間加權約束初始速度模型，得到高精度速度模型Vcon(x,y,z)，構建方法通過下式實現

(1)

圖1 模型數據實際測井速度曲線與低頻測井曲線

測井約束項的權重是服從Logistic函數形態分布的，則測井約束空間權重因子可以表示為

(2)

2.2 測井約束全波形反演方法

通常地球物理反演問題都基于Bayes框架[1]。地震波的正傳播過程可抽象表示為

d=G(m)

(3)

式中：m為地球物理參數矢量，如地下地層介質速度、密度或其他彈性參數；d為地震數據；G(m)為基于地層介質模型m的地震波傳播過程，描述了地震波場的傳播。研究地層介質模型m中地震波傳播現象，得到觀測地震數據d的過程稱為地震波正演，由地震數據d推測地層介質模型參數m的過程稱為地震反演。

地震反演過程可以表示為

m=G-1(d)

(4)

全波形反演的誤差泛函一般表示為[28]

(5)

式中dobs為觀測數據。式(5)為最小二乘目標函數。假設觀測數據與模型數據服從高斯分布，在目標函數中加入先驗模型項，定義含先驗模型信息的目標泛函為

(6)

式中：mprior為測井約束的先驗模型；λ為權重參數；CD和CM分別為觀測數據和模型的協方差。

本文采用梯度引導類方法[29]求解反演問題，梯度引導類方法通過求取目標泛函對模型的導數尋找迭代更新方向。基于聲波方程的假設，模型參數為縱波速度場m，觀測數據為采集的地震炮集波場dobs，由地震波正演模擬得到的波場為u=G(m)，G表示地震波場正傳過程。則誤差泛函為

(7)

式中：t為波場傳播時間采樣點；xr為與炮點xs相對應的檢波點坐標。

基于線性化假設，利用變分原理推導泛函梯度，得到泛函對速度模型的梯度為

(8)

式中：G*表示波場殘差的反傳過程；u(t,xs,xr)為炮點xs對應的檢波點xr位置的波場。

在反演過程中，權重參數λ的選取是動態的，誤差泛函中數據誤差項會隨著迭代次數的增加而不斷減小。因此，為了保證先驗模型項與數據誤差項對誤差泛函貢獻的平衡，要求λ隨著迭代次數的增加而減小。圖3為測井約束全波形反演流程圖。

圖3 測井約束全波形反演流程圖

3 模型試算

本文采用EAGE推覆體模型進行模型試算，抽取600、1400m處垂向單道模型數據作為測井數據(圖4a)。

首先，對測井數據進行分析與處理，將測井數據由深度域轉到時間域，并分析測井數據與地震記錄數據的頻譜特征，對測井數據進行濾波處理，得到低頻信息，進而將濾波后的時間域測井曲線轉換到深度域(圖5)。然后根據本文提出的方法，用得到的低頻測井數據對常規初始速度模型(圖4b)進行加權約束，建立新的基于測井約束的高精度初始速度模型(圖4c)。 圖6為常規初始速度模型曲線與測井約束模型曲線。由圖可見，測井約束速度模型依然是一個平滑的模型，其精度比常規初始速度模型高。

圖4 真實速度模型、反演初始速度模型及先驗模型

圖5 實際測井曲線與濾波后低頻測井曲線

圖6 常規初始速度模型曲線與測井約束模型曲線

通過計算測井數據得到測井約束全波形反演的先驗模型(圖4d)。在水平方向上對井數據進行線性距離加權內插和外推得到插值剖面，然后對其進行輕微的高斯平滑。以平滑后的插值速度作為先驗模型約束反演目標泛函，充分利用了測井數據的先驗信息，降低了反演陷入局部極小的可能。

為了驗證所提方法的反演效果，分別使用常規初始速度模型及本文方法獲得的測井約束速度模型進反演。圖7為常規初始速度模型反演結果與測井約束反演結果，圖8為常規反演與測井約束反演的誤差下降曲線。可見測井約束反演方法可以較快地收斂得到最優反演結果。為了更直觀地比較兩種方法的反演效果，抽取1400m位置處垂向單道速度反演結果進行對比(圖9)。可見通過測井約束反演得到的結果更接近真實速度模型，尤其在深層能得到較好的反演結果。

圖7 常規初始速度模型反演結果(a)與測井約束反演結果(b)

震源和檢波點都位于地表處，共有25炮，第一炮位于10m處，炮間距為80m，第一個檢波點位于0處，檢波點間隔為10m，地震記錄采樣間隔為0.0005s，采樣點數為1500，總記錄時間為0.75s，采用高階有限差分法進行二維聲波正演，理論子波是主頻為30Hz的雷克子波

圖8 常規反演與測井約束反演的誤差下降曲線

圖9 1400m位置處垂向單道速度反演結果

4 實際資料應用

將本方法應用于A工區二維資料。地震資料為一條炮線數據(圖10)，反演區域面積為12.45km×4km，初始速度為通過偏移速度分析獲得的背景速度場，網格尺寸為10m×10m，在炮線3.55km處有一口探井，測井縱波數據深度范圍為1890～3220m。由于該區地震數據帶寬有限，低頻分量較少，地震數據主頻為25Hz，低頻信息嚴重缺失(圖12a)。對測井數據進行處理，圖11為經過井震標定后的深度域測井數據。根據時深關系，將井曲線由深度域轉換到時間域，再對測井數據進行頻譜分析，可以看出測井數據既有高分辨率的地層介質速度信息又有豐富的低頻分量信息(圖12b)。對測井數據進行濾波處理，提取測井數據的低頻分量，并選取截止頻率為5Hz，在一般情況下截止頻率不能取太高，以免在反演一開始就陷入局部極值。圖13為時間域測井速度曲線以及濾波后的低頻測井曲線。通過時深轉換將時間域低頻測井曲線變換到深度域。

必須注意到，一般實際測井深度有限，測井數據并非從地表貫穿井底，測井曲線的范圍遠小于地震資料的時間長度，則低頻測井曲線不能完全覆蓋研究區域的頂、底界面。如果僅簡單地在縱波速度范圍對初始速度模型進行校正，則可能在真實井曲線數據范圍的頂、底處會產生速度突變，為了防止這一現象，需要對測井數據進行上、下延拓，并逐步在向外延伸區域以常規初始速度為準。圖14為測井約束曲線對比。由圖可見，低頻測井曲線相對于常規初始模型速度更精確，可利用延拓后的低頻測井曲線約束常規初始速度模型。

圖15為常規初始速度模型與測井約束初始速度模型。由圖可見，使用測井約束低頻模型構建方法，利用低頻測井數據對常規初始速度模型(圖15a)進行空間加權約束，可得到測井約束下的初始速度模型(圖15b)。

圖10 單炮地震記錄

圖11 經過井震標定后的深度域測井數據采樣間隔為0.125m

圖12 地震數據 (a)、測井數據(b)頻譜

圖13 時間域測井速度曲線及濾波后低頻測井曲線對比

圖14 測井約束曲線對比

圖16為實際數據測井約束反演與常規反演的誤差下降曲線。由圖可見，利用測井約束全波形反演方法能較快地收斂得到最優結果。圖17為常規初始速度模型反演結果與測井約束反演結果。由圖可見：由于震源點主要集中于5～7km范圍，故模型兩側覆蓋次數較低，反演效果較覆蓋次數高的區域差；測井解釋結果表明，在井位置的3000m深度處存在一套砂巖層(圖中黑色箭頭所指處)，與常規不加約束的反演結果(圖17a)相比，使用測井約束的反演結果能更清楚地展示砂巖層，具有較高的分辨率，反演效果更好(圖17b)。

圖15 測井曲線附近目標層區域的常規初始速度模型(a)與測井約束初始速度模型(b)

圖16 實際數據測井約束反演與常規反演的誤差下降曲線

圖17 常規初始速度模型反演結果(a)與測井約束反演結果(b)

5 結論與討論

針對全波形反演對初始模型及低頻信息的依賴性，本文研究了基于測井約束的初始模型構建方法及全波形反演方法，在算法實現的基礎上，通過模型數據測試證明了文中方法的有效性，并將該方法成功應用于實際地震數據，主要得到以下認識。

(1)由于采集條件的限制，地震數據缺少低頻信息。通過提取測井數據中可靠的低頻信息可以校正常規初始速度模型中的低頻信息，提高初始速度模型的精度。同時，在全波形反演目標泛函中加入先驗模型約束項，將測井約束信息引入反演優化過程，可以使反演向著給定的方向進行，避免反演陷入局部極小，可提高反演的穩定性。

(2)測井約束模型構建能在一定程度上提高初始模型的精度，但當構造復雜或地層傾角較大時存在破壞原始構造信息的風險，需要進一步改進模型構建方法。其中權重參數需要隨著迭代次數增加動態變化，以防止先驗模型項對誤差反演貢獻過大而在梯度求取中產生新的局部極值，導致反演無法收斂。

(3)模型測試和實際資料應用得到了較好的反演結果，證明了本文方法的有效性。但在實際應用中存在許多亟待解決的問題，如子波的選取、地震資料信噪比等。

[1] Tarantola A.Inverse Problem Theory,Methods for Data Fitting and Model Parameter Estimation.Elsevier Science Publ Co Inc,Amsterdam and New York，1987.

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

[3] Petrick W R,Borup D T,Steven A J et al.Seismic borehole tomography using full waveform inversion.SEG Technical Program Expanded Abstracts,1988,7:1250-1252.

[4] Pratt R G,Shin C,Hick G J.Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion.Geophysical Journal International,1998,133(2):341-362.

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

[6] Shin C,Cha Y H.Waveform inversion in the Laplace domain.Geophysical Journal International,2008,173(3):922-931.

[7] Sirgue L,Barkved O I,Van Gestel J P et al.3D waveform inversion on Valhall wide-azimuth OBC.EAGE 71st Conference and Exhibition Extended Abstracts,2009.

[8] Warner M,Umpleby A,Stekl I et al.3D full-wavefield tomography:imaging beneath heterogeneous overburden.EAGE 72st Conference and Exhibition Extended Abstracts,2010.

[9] Pratt R G,Song Z M,Williamson P et al.Two-dimensional velocity models from wide-angle seismic data by waveform inversion.Geophysical Journal International,1996,124(2):323-340.

[10] Smithyman B,Pratt G,Hayles J et al.Detecting near-surface objects with seismic waveform tomography.Geophysics,2009,74(6):WCC119-WCC127.

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

[12] 楊午陽,王西文,雍學善等.地震全波形反演方法研究綜述.地球物理學進展,2013,28(2):766-776.

Yang Wuyang,Wang Xiwen,Yong Xueshan et al.The review of seismic full waveform inversion method.Progress in Geophysics,2013,28(2):766-776.

[13] Shah N,Warner M,Nangoo T et al.Quality assured full-waveform inversion:Ensuring starting model adequacy.SEG Technical Program Expanded Abstracts,2012,31.

[14] Weibull W,Arntsen B,Nilsen E.Initial velocity models for full waveform inversion.SEG Technical Program Expanded Abstracts,2012,31.

[15] 敖瑞德,董良國,遲本鑫.不依賴子波、基于包絡的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.

[16] Gong B,Chen G,Yingst D et al.3D waveform inversion based on reverse time migration engine.SEG Technical Program Expanded Abstracts,2008,27:1900-1904.

[17] Baeten G,Maag J W,Plessix R E et al.The use of low frequencies in a full-waveform inversion and impedance inversion land seismic case study.Geophysical Prospecting,2013,61(4):701-711.

[18] 段超然,韓立國.主成分分析波場重構反演與全波形反演聯合速度重構.石油地球物理勘探,2016,51(6):1134-1140.

Duan Chaoran,Han Liguo.A joint velocity reconstruction method:principal component analysis based wavefield-reconstruction inversion combined with full waveform inversion.OGP,2016,51(6):1134-1140.

[19] 王毓瑋,董良國,黃超等.降低彈性波全波形反演強烈非線性的分步反演策略.石油地球物理勘探,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.

[20] 吳國忱,田軍.基于尺度約束項的AVO反演方法.石油地球物理勘探,2014，49(6):1157-1164.

Wu Guochen,Tian Jun.AVO inversion based on scale constraints.OGP,2014,49(6):1157-1164.

[21] 馮國峰,韓波,劉家琦.二維波動方程約束反演的大范圍收斂廣義脈沖譜方法.地球物理學報,2003,46(2):265-270.

Feng Guofeng,Han Bo,Liu Jiaqi.Widely convergent generalized pulse spectrum methods for 2D wave equation inversion.Chinese Journal of Geophysics,2003,46(2):265-270.

[22] Hu W,Abubakar A,Habashy T M.Simultaneous multifrequency inversion of full-waveform seismic data.Geo-physics,2009,74(2):R1-R14.

[23] 王守東.井間地震資料全變差正則化波形反演.石油物探,2011,50(3):234-240.

Wang Shoudong.Total variation regularization waveform inversion of crosswell seismic data.GPP,2011,50(3):234-240.

[24] Han B,Fu H S,Liu H.A homotopy method for well-log constraint waveform inversion.Geophysics,2007,72(1):R1-R7.

[25] 傅紅筍,曹莉,韓波.測井約束地震波形反演的同倫攝動法.地球物理學報,2012,55(9):3173-3179.

Fu Hongsun,Cao Li,Han Bo.A homotopy perturbation method for well log constrained seismic waveform inversion.Chinese Journal of Geophysics,2012,55(9):3173-3179.

[26] Asnaashari A,Brossier R,Garambois S et al.Regularized seismic full waveform inversion with prior model information.Geophysics,2013,78(2):R25-R36.

[27] 成景旺,呂曉春,顧漢明等.基于柯西分布的頻率域全波形反演.石油地球物理勘探,2014,49(5):940-945.

Cheng Jingwang,Lü Xiaochun,Gu Hanming et al.Full waveform inversion with Cauchy distribution in the frequency domain.OGP,2014,49(5):940-945.

[28] 任浩然.聲介質地震波反演成像方法研究[學位論文].上海:同濟大學,2011.

Ren Haoran.A Study of the Seismic Wave Inversion Method for Imaging in the Acoustic Media[D].Tongji University,Shanghai,2011.

[29] Mora P.Nonlinear two-dimensional elastic inversion of multioffset seismic data.Geophysics,1987,52(9):1211-1228.

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

本文于2016年10月28日收到，最終修改稿于2017年8月21日收到。

本項研究受國家重大科技專項(2016ZX05024-001-008)資助。

1000-7210(2017)06-1184-09

杜澤源，吳國忱，王玉梅.基于測井約束的地震全波形反演方法.石油地球物理勘探，2017,52(6):1184-1192.

P631

A

10.13810/j.cnki.issn.1000-7210.2017.06.008

(本文編輯：劉勇)

杜澤源 博士研究生，1988年生；2013年獲中國石油大學(華東)勘查技術與工程專業學士學位。現在該校地球科學與技術學院攻讀地質資源與地質工程專業博士學位，主要從事地震波正、反演方面的研究。