應(yīng)用反射理論的最小二乘逆時偏移成像

段心標 王華忠 鄧光校

(①同濟大學(xué)海洋與地球科學(xué)學(xué)院波現(xiàn)象與智能反演成像研究組，上海 200092；②中國石化石油物探技術(shù)研究院，江蘇南京 211103； ③中國石化西北油田分公司，新疆烏魯木齊 830011)

0 引言

從地震數(shù)據(jù)中準確提取并利用復(fù)雜的地下巖性與儲層信息，是油氣勘探的關(guān)鍵環(huán)節(jié)。地震波成像是連接地震觀測數(shù)據(jù)與儲層解釋的橋梁，常規(guī)偏移成像技術(shù)[1-5]主要目的是獲得能精細描述地質(zhì)構(gòu)造的地震成果資料，用于部署和落實井位。

隨著油氣勘探程度的日趨深入和地震成像技術(shù)的不斷進步，直接估算地下彈性參數(shù)的反演成像方法受到越來越多的關(guān)注。早在20世紀60年代，Backus等[6-7]從數(shù)學(xué)角度研究了地球物理反問題，奠定了早期地球物理反演的理論基礎(chǔ)。Tarantola[8-9]基于Bayes估計理論，建立了一套地震波全波形反演的理論方法體系，推動了地震反演成像技術(shù)的快速發(fā)展。然而，全波形反演具有很強的非線性特征，難以由觀測數(shù)據(jù)直接反演得到較準確的介質(zhì)彈性參數(shù)，無法滿足油藏描述的精度需求。實際中應(yīng)用的依然是經(jīng)典的地震波成像處理和儲層描述流程，而如何獲得更高偏移成像質(zhì)量是其中重要環(huán)節(jié)。

最小二乘偏移[10-18]是常規(guī)偏移的后續(xù)發(fā)展，其核心思想是在線性反演理論下求解模型空間的精確解，相當于在常規(guī)偏移構(gòu)造成像基礎(chǔ)上，考慮了振幅補償和波形校正，從而得到更保幅和更高分辨率的成像結(jié)果。

經(jīng)典最小二乘偏移以散射波線性化表達為基礎(chǔ)，認為模型參數(shù)是速度(或慢度)的高波數(shù)擾動，其成像實質(zhì)是估計介質(zhì)的散射強度。Zhang等[19-20]提出基于逆時偏移疊加剖面的反偏移方法，并形成相應(yīng)的最小二乘偏移技術(shù)，但并未給出反偏移方法的詳細理論依據(jù)。陳生昌等[21-22]提出地震反射波數(shù)據(jù)偏移成像方法和相應(yīng)的反射波動方程最小二乘偏移方法，并引入“反射率”概念，即最小二乘偏移的實質(zhì)是估計反射面上的反射率。本文認為：因地下介質(zhì)以層狀為主，獲取高保真反射系數(shù)是地震成像走向儲層刻畫的關(guān)鍵，估算不同角度的反射系數(shù)是最小二乘偏移的目標。然而，考慮到反演不同角度反射系數(shù)的計算量巨大，且高質(zhì)量偏移疊加成像依然是當今實際應(yīng)用的主要目標，因此以最小二乘偏移獲得忽略角度信息的平均反射系數(shù)具有現(xiàn)實意義。

本文首先簡介線性最小二乘偏移方法基礎(chǔ)理論，重點討論散射理論線性化表達方法和反射理論線性化表達方法，并分別建立散射強度、反射系數(shù)與地震數(shù)據(jù)的線性關(guān)系式； 然后從偏移成像實際應(yīng)用需求出發(fā)，選擇基于反射理論并忽略角度信息的Kirchhoff正演反偏移表達方法，實現(xiàn)與實際應(yīng)用需求相一致的最小二乘逆時偏移技術(shù)； 最后，通過模型數(shù)據(jù)和實際資料驗證本文方法的效果。

1 方法原理

1.1 線性最小二乘偏移方法

基于Bayes估計理論，地震反演成像可表示為后驗概率密度最大化問題； 在地震反演通常假設(shè)條件下，進而可轉(zhuǎn)化為如下目標函數(shù)最小化問題

(1)

式中：m為地下介質(zhì)參數(shù)，如速度、散射強度、反射系數(shù)或波阻抗等；dobs為觀測數(shù)據(jù)；L(·)表示波場傳播算子。該式非線性反演問題求解難度大，需對非線性反演問題逐步進行線性化，依次反演參數(shù)的低、中、高波數(shù)信息。最小二乘偏移是一種線性化的地震反演成像方法，用于估計介質(zhì)參數(shù)的高波數(shù)成分，其目標函數(shù)可表示為

(2)

式中L為矩陣形式的線性化波場傳播算子。

令目標函數(shù)梯度為零，則有如下的法方程

LTLm=LTdobs

(3)

最小二乘偏移成像結(jié)果表示為

(4)

式中：LTL為Hessian矩陣，可用H表示；LTdobs是對觀測數(shù)據(jù)進行常規(guī)偏移成像。最小二乘偏移是Hessian矩陣的逆H-1作用于常規(guī)偏移成像，能夠消除Hessian矩陣導(dǎo)致的成像模糊化效應(yīng)，提高成像分辨率和保幅性，減少成像假象。

由于“目標函數(shù)梯度為零”僅是一種理想假設(shè)，且Hessian矩陣非常龐大，其計算和存儲都很困難，Hessian矩陣的求逆更難以完成，因而無法直接求解法方程以實現(xiàn)最小二乘偏移。為此，通常采用迭代法求解最小二乘偏移問題，假設(shè)目標函數(shù)在初始值附近滿足二次型，從初始解逐步逼近真實解。

1.2 地震波場的線性化表達

1.2.1 散射理論的Born正演線性化表達

一階Born近似假設(shè)下，散射波場可寫成

(5)

(6)

進而，散射波場可表示為

G0(x′,xs,ω)dx′

(7)

式(7)定義了一個Fredholm第一類積分方程，它描述了一個線性問題。其積分核可寫成

(8)

在散射假設(shè)下，可認為地表接收到的攜帶地下信息的信號就是散射波場記錄，則可將式(7)寫成如下線性算子表達式

Lm=dobs

(9)

此式中的m設(shè)定為慢度平方擾動，則該式表達了散射強度與地震數(shù)據(jù)之間的線性關(guān)系。

1.2.2 反射理論的Kirchoff正演線性化表達

反射波理論假設(shè)下，地下界面Ω上入射波場引起的二次場可表示為

(10)

地下界面Ω上波場的法向?qū)?shù)為

(11)

利用WKBJ近似方法[23]，有

(12)

由Green定理，反射面上方地面上任意一點x的反射波場可寫成如下積分式

(13)

式中Gr(x,x′,ω)是反射波傳播的Green函數(shù)(圖1)。

圖1 波場傳播示意圖

將式(10)～式(12)代入式(13)，可得

(14)

根據(jù)程函方程關(guān)系，式(14)可進一步寫成

Gs(x′,xs,ω)Gr(x,x′,ω)dΩ

(15)

式中s(x′)為慢度矢量。

式(15)建立了地震數(shù)據(jù)與不同方向反射系數(shù)之間的聯(lián)系[24-25]。然而，在緩變速度結(jié)構(gòu)中，不同角度的反射波能量較均衡； 在速度劇變結(jié)構(gòu)中，不同角度的反射系數(shù)難以估計。同時，考慮到直接基于式(15)進行波場模擬的計算量巨大，因此建立地震數(shù)據(jù)與平均反射系數(shù)之間的關(guān)系就具有現(xiàn)實意義。若不考慮背景速度對振幅的影響，式(15)就退化為

(16)

式中：D表示散射體；R′表示平均反射系數(shù)。

在反射假設(shè)下，認為地表接收到的攜帶地下信息的信號就是反射波場記錄。式(16)同樣可寫成線性算子(式(9))表達形式，此時m表示忽略地下角度信息的平均反射系數(shù)，其積分核為

(17)

1.3 最小二乘逆時偏移實現(xiàn)方法

逆時偏移技術(shù)現(xiàn)已在實際處理中得到廣泛應(yīng)用，其成像基礎(chǔ)是Claerbout[26]給出的成像方法：波場外推+成像條件。成像條件為：成像點處，入射波的到達時等于反射波的出發(fā)時。顯然，常規(guī)逆時偏移是估計地下介質(zhì)的反射系數(shù)。由于地下介質(zhì)以層狀為主，以定位反射界面的位置或產(chǎn)生角度反射系數(shù)道集為目標的反射系數(shù)成像是當前地震成像的主流需求。

如前所述，由于角度反射系數(shù)反演的計算量巨大，因而實際應(yīng)用中，最小二乘偏移的目的是在常規(guī)逆時偏移基礎(chǔ)上進一步提高偏移疊加成像質(zhì)量，使成像結(jié)果分辨率更高，成像假象更少。

在反演成像方法中，反偏移采用Kirchhoff正演思路，具體迭代反演步驟如下。

(1)以設(shè)定的偏移速度場對觀測數(shù)據(jù)dobs作逆時偏移成像，作為最小二乘偏移的初始值m0。

(2)對第i次迭代后偏移成像值mi(i=0,1,2,…)做反偏移，正演模擬預(yù)測數(shù)據(jù)dcal，反偏移公式為

(18)

(3)計算觀測數(shù)據(jù)dobs與預(yù)測數(shù)據(jù)dcal的殘差波場δd，當殘差滿足設(shè)定條件時迭代停止。

(4)對δd進行逆時偏移成像，得到目標泛函的梯度gi。

(5)計算共軛梯度方向

(19)

(6)計算迭代步長

(20)

(7)更新偏移成像值

mk+1=mk-μk+1hk

(21)

(8)轉(zhuǎn)到步驟(2)，展開新一輪迭代。

由于迭代算法計算量巨大，三維最小二乘逆時偏移中采用了GPU加速計算技術(shù)，在此不再贅述。

2 數(shù)值模型試驗及實際應(yīng)用

首先利用二維Sigsbee2a模型(圖2a)做數(shù)值試驗，驗證本文方法效果； 進而將三維最小二乘逆時偏移技術(shù)應(yīng)用于實際地震資料成像處理。

圖2 Sigsbee2a模型及其成像結(jié)果

2.1 Sigsbee2a模型數(shù)據(jù)

為了便于計算，把Sigsbee2a經(jīng)典模型速度值單位轉(zhuǎn)為m/s，并重新定義速度網(wǎng)格，新網(wǎng)格的尺寸為3200×1200，間隔為5m×5m。對該模型正演模擬120炮數(shù)據(jù)，炮間距為100m，檢波點間隔為20m，炮檢距范圍是-2000～2000m，記錄長度為7s。

利用本文所提最小二乘逆時偏移方法進行反演迭代，得到第20次迭代結(jié)果(圖2c)，并與常規(guī)逆時偏移剖面(圖2b)進行對比； 圖3是目標泛函迭代誤差分析曲線，可見反演逐步收斂。從整體剖面及鹽丘左下局部放大圖(圖4)的對比可見，相對于常規(guī)逆時偏移方法，本文方法的成像效果(圖4b)顯著改善，分辨率變高，深部成像更清晰。

圖3 迭代誤差曲線

圖4 Sigsbee2a模型的常規(guī)逆時偏移(a)和最小二乘逆時偏移(b)成像的局部放大

2.2 實際資料

本次試驗所用實際數(shù)據(jù)取自中國西北M油田，測試處理的目的是提高儲層的成像精度，力求對奧陶系縫洞體及小斷裂精確成像。所選數(shù)據(jù)滿覆蓋面積為401km2，共48000炮，成像面元是30m×30m。利用100GPU節(jié)點(每節(jié)點4塊K10卡)3次迭代，共用時59天。

圖5和圖6是兩條不同測線的常規(guī)逆時偏移與最小二乘逆時偏移成像剖面對比，可見最小二乘逆時偏移剖面中“串珠”成像收斂更好，對深層大斷裂、層間小斷層的刻畫更清晰，T78(奧陶系蓬萊壩組頂界面)以下成像能量更強，成像分辨率得到提高。對比奧陶系目的層的成像頻譜(圖7)，可見最小二乘逆時偏移成像拓寬了頻譜(圖7b，頻帶展寬5Hz)； T74(一間房組頂界面)下40～120ms均方根振幅，也證實本文方法成像結(jié)果(圖8b)對“串珠”刻畫更清晰。

以上結(jié)果均展示了最小二乘逆時偏移成像在提高成像分辨率、振幅均衡性及成像精度等方面的效果。

圖5 L1測線的常規(guī)逆時偏移(a)和最小二乘逆時偏移(b)成像剖面

圖6 L2測線的常規(guī)逆時偏移(a)和最小二乘逆時偏移(b)成像剖面

圖7 常規(guī)逆時偏移(藍線)和最小二乘逆

3 結(jié)論與討論

最小二乘偏移是一種線性化的地震反演成像方法，線性近似可以是“背景速度+速度擾動”模式，也可以是“背景速度+反射系數(shù)”模式。由于地下介質(zhì)以層狀為主，現(xiàn)今業(yè)界常規(guī)偏移成像的主要目標是反射系數(shù)。本文分析了兩種線性化模式的差異，選擇基于忽略角度信息的反射系數(shù)的反偏移正問題表達方法，提出一種估計平均反射系數(shù)的最小二乘逆時偏移方法。理論模型測試結(jié)果表明本方法能提高成像分辨率，深部成像更清晰。將研發(fā)的三維最小二乘逆時偏移技術(shù)應(yīng)用于中國西部探區(qū)資料處理，所得“串珠”成像比常規(guī)逆時偏移結(jié)果收斂更好，深層大斷裂、層間小斷層的刻畫更清晰，成像分辨率獲得提高，充分體現(xiàn)了最小二乘逆時偏移成像在提高成像分辨率和振幅均衡性及提高成像精度方面的優(yōu)勢。

圖8 常規(guī)逆時偏移(a)與最小二乘逆時偏移(b)的均方根振幅平面屬性對比

然而，目前最小二乘逆時偏移實際應(yīng)用效果仍不甚理想，主要原因包括反偏移得到的無噪聲一次波模擬數(shù)據(jù)與實際觀測記錄差異較大、實際地震子波未知、現(xiàn)有建模精度達不到最小二乘逆時偏移要求及三維技術(shù)計算量大等。為了降低最小二乘逆時偏移對速度模型和地震子波的依賴，近年也發(fā)展了動態(tài)拉伸時差校正方法和空變子波估計方法，一定程度上減弱了反演的病態(tài)性。但考慮到反問題的復(fù)雜性，這些應(yīng)用瓶頸問題并未得到有效解決，它們會持續(xù)伴隨最小二乘逆時偏移的發(fā)展過程，有待逐步克服和解決。