徐 凱,孫贊東

(1.中國石油大學(xué)(北京),地質(zhì)地球物理綜合研究中心,北京102249;2.中國石油化工股份有限公司石油勘探開發(fā)研究院,北京100083)

隨著地震勘探開發(fā)的深入發(fā)展,復(fù)雜地質(zhì)條件下的地震勘探已經(jīng)成為目前學(xué)術(shù)界與工業(yè)界的研究重點與熱點。地球介質(zhì)廣泛存在粘滯性,地震波在地下傳播時會發(fā)生衰減,深層能量呈指數(shù)形式下降,如若忽略掉粘彈介質(zhì)的吸收衰減作用,會造成成像結(jié)果偏差大,頻散嚴(yán)重,深層能量不足且保幅性差,進(jìn)而嚴(yán)重影響后續(xù)的反演解釋結(jié)果。

針對粘介質(zhì)的吸收衰減作用,人們首先采用反Q濾波方法來補償高頻能量,但是該方法只適用于簡單介質(zhì)。隨著偏移技術(shù)的發(fā)展,利用Q偏移方法來恢復(fù)因吸收衰減作用而損失的能量逐漸成為研究的熱點。如RIBODETTI等[1]提出了基于射線理論的Q偏移方法,主要是利用射線理論做高頻近似。隨著逆時偏移技術(shù)的發(fā)展,基于波動方程理論的Q偏移也取得了長足的進(jìn)步。CAUSSE等[2]將粘聲介質(zhì)引入逆時偏移中,利用時間域波動方程實現(xiàn)了簡單的粘聲逆時偏移。ZHU等[3]在成像之前先進(jìn)行相位校正,再做基于粘聲介質(zhì)的逆時偏移,較好地解決了相位異常的問題。ZHANG等[4]推導(dǎo)了基于粘聲介質(zhì)常Q模型的微分方程,有效地補償了深層能量衰減與相位頻散影響。同時,針對高頻不穩(wěn)定的情況引入了規(guī)則化算子,有效解決了反向延拓過程中的不穩(wěn)定問題。與單程波偏移方法相比,逆時偏移有保幅方面的優(yōu)勢。

傳統(tǒng)的逆時偏移成像只是簡單的數(shù)據(jù)共軛轉(zhuǎn)置,并不是正演的逆,而最小二乘偏移是在非線性問題比較弱的情況下,利用Born近似對非線性問題局部尋優(yōu)迭代求解的一種反演方法。但是最小二乘偏移的計算量較大,需要大量的存儲空間,因此其發(fā)展較為緩慢。很多學(xué)者在反演思想的框架下,將逆時偏移與最小二乘偏移思想相結(jié)合。TARANTOLA[5]最早提出了最小二乘偏移思想,LEBRAS等[6]、LAMBARE等[7]在其基礎(chǔ)上進(jìn)行了完善,使用最速下降法反演相對于背景速度的速度擾動。NEMETH等[8]首先將最小二乘算法引入地震偏移以求取地下反射系數(shù),然后將Kirchhoff最小二乘偏移方法應(yīng)用于不完整的反射地震數(shù)據(jù)偏移。目前,許多學(xué)者將單程波方程和逆時偏移與最小二乘理論相結(jié)合,實現(xiàn)了基于波動方程的最小二乘偏移[9-16]。黃建平等[17]將最小二乘偏移方法應(yīng)用于碳酸鹽巖裂縫型儲層預(yù)測,詳細(xì)分析了最小二乘偏移理論在碳酸鹽巖成像中的優(yōu)勢,并將其應(yīng)用于地下裂縫模型測試。劉玉金等[18]通過引入預(yù)條件算子實現(xiàn)了擴展最小二乘逆時偏移,進(jìn)一步提高了最小二乘逆時偏移的精度和確定性。QU等[19]實現(xiàn)了彈性波最小二乘逆時偏移,更加真實地模擬了地下實際情況,取得了較好的成像效果。

本文在反演思想的框架下,引入最小二乘理論,給出了基于廣義標(biāo)準(zhǔn)線性固體模型(GSLS)的粘聲介質(zhì)波動方程,并推導(dǎo)了對應(yīng)的反偏移算子與伴隨算子,將粘聲介質(zhì)逆時偏移與最小二乘理論結(jié)合,實現(xiàn)了基于粘聲衰減補償?shù)淖钚《四鏁r偏移(LSRTM)。通過層狀模型與Marmousi模型試算驗證了該方法的正確性及有效性。

1 方法原理

1.1 粘聲介質(zhì)波動方程

基于GSLS模型的波動方程可以寫成如下形式:

(1)

采用(1)式的交錯網(wǎng)格有限差分形式可以進(jìn)行波動方程模擬?；谝浑A速度-應(yīng)力交錯網(wǎng)格的有限差分模擬是目前使用較多的交錯網(wǎng)格模擬方式,交錯網(wǎng)格可以有效代替普通矩形網(wǎng)格,更好地模擬地下真實情況,提高正演模擬的精度。

1.2 粘聲介質(zhì)逆時偏移所對應(yīng)的反偏移算子

最小二乘逆時偏移的關(guān)鍵步驟是反偏移,主要是通過背景波場與慢度擾動的互相關(guān)產(chǎn)生所需要的反射波。反偏移并不像一般的正演過程那樣需要界面信息,可以在光滑的背景速度場中獲取反射波。

根據(jù)波場疊加原理,有:

(2)

這里,p0為背景波場,ps為擾動波場。背景波場的傳播方程與普通波場傳播方程一致,但是傳播速度存在較大區(qū)別。背景波場的速度場是平滑速度場,方程表示為:

(3)

定義慢度擾動為:

(4)

式中:v為全波場速度,v0為背景波場速度。擾動波場的傳播方程也與普通方程類似,而且速度場也是平滑后的速度場。擾動波場傳播與背景波場傳播最大的不同是震源不同,擾動波場的震源不再是傳統(tǒng)的地震子波震源,而是由背景波場與慢度擾動的乘積代替:

(5)

為了簡化上述方程,對該方程應(yīng)用Born近似,得到基于Born近似的粘聲介質(zhì)最小二乘逆時偏移反偏移算子:

(6)

擾動波場ps(x,t)由兩個過程計算得到:一是震源f(t)和背景速度v0產(chǎn)生p0(x,t);二是利用m(x)[?2p0(x,t)]/?2t和背景速度v0產(chǎn)生ps(x,t)。從公式(6)可以看出,背景場中沒有擾動場;擾動場中只存在一次擾動場。但Born近似依賴于地下尺度大小,當(dāng)存在小尺度體時,Born近似效果較好。這是當(dāng)前最小二乘疊前深度偏移成像的基礎(chǔ)。

1.3 粘聲介質(zhì)逆時偏移對應(yīng)的伴隨算子

粘聲介質(zhì)逆時偏移包括兩個過程:波的正向傳播與波的反向傳播。隨著時間的增加,粘聲介質(zhì)地震波在地下的正向傳播是穩(wěn)定的。在波的反推過程中,隨著時間的減少往往會出現(xiàn)不穩(wěn)定問題。解決該問題的方法主要有兩種:一種是在反傳過程中引入規(guī)則化算子,壓制反傳過程中的高頻震蕩。另一種是通過正演波動方程的伴隨方程實現(xiàn)波場的逆時延拓,有效地避免反推過程中的不穩(wěn)定問題。本文采用第二種方法解決不穩(wěn)定問題。

去掉震源項,正演算子T波動方程可以寫成:

(7)

令正演算子T的伴隨算子為T*,兩者的關(guān)系滿足=,則:

(8)

因此,正演的伴隨算子T*也是一個二階微分方程,表示為:

(9)

這里,

(10)

對公式(10)兩邊同時對時間求導(dǎo),則:

(11)

(12)

公式(9)和公式(12)構(gòu)成了粘聲介質(zhì)反向傳播方程。為了便于粘聲介質(zhì)微分方程進(jìn)行交錯網(wǎng)格數(shù)值模擬,引入ux和uz,將其空間二階微分方程降為空間一階微分方程,對空間微分進(jìn)行交錯網(wǎng)格求解,再建立復(fù)頻移完美匹配層(PML)邊界方程,進(jìn)行數(shù)值模擬。

1.4 粘聲介質(zhì)最小二乘理論

最小二乘逆時偏移的主要思想是利用波動方程來模擬反射系數(shù)的響應(yīng),得到正演記錄,然后將其與輸入的實際地震數(shù)據(jù)對比,如果存在誤差,就采用最小二乘方法校正,直至誤差減小到一個合理的范圍內(nèi)。

定義最小二乘優(yōu)化的目標(biāo)函數(shù):

(13)

式中:m代表反射系數(shù)模型,即成像結(jié)果;L為Born近似下的反偏移算子;d為實際輸入數(shù)據(jù)(野外炮記錄);a為正則化參數(shù);H1和H2代表1范數(shù)和2范數(shù)。當(dāng)目標(biāo)函數(shù)f(m)最小化時,我們認(rèn)為此時求得的m就是最優(yōu)解。求解(13)式的首要問題是計算目標(biāo)函數(shù)關(guān)于參數(shù)變化的導(dǎo)數(shù),即目標(biāo)函數(shù)的梯度,本文采用以下變分方法。

對任意的介質(zhì)擾動δm和擾動步長t,總存在如下目標(biāo)函數(shù)變分:

(14)

對上式定義的變分求偏導(dǎo)數(shù):

(15)

m(k+1)=m(k)-mg(m)=m(k)-mL*

(Lm-dobs)=m(k)-mL*(Lm(k)-dobs)=

m(k)-mL*(dcal-dobs)

(16)

該最優(yōu)化問題也可以采用共軛梯度法來求解,實現(xiàn)思想與最速下降法類似,就是以模型m0為初始模型,建立迭代公式,求取梯度與步長,最后檢測誤差是否在允許范圍內(nèi)。在目標(biāo)函數(shù)對模型參數(shù)m求導(dǎo)時,誤差的梯度為:

(17)

其迭代公式為:

式中:α為迭代步長。

在最小二乘疊前深度偏移中,一般令梯度等于零(實際梯度不等于零)時,可以得到線性方程:

(20)

該方程求解的核心問題是求解Hessian矩陣的逆,但是Hessian矩陣的逆無法顯式求解,因此難以得到其精確解。由于迭代求解的算法不是直接顯式求解Hessian矩陣的逆,而是通過隱式求解的方法得到,因此最小二乘的方法被廣泛應(yīng)用于偏移與全波形反演中。

具體實現(xiàn)流程如圖1所示。

圖1 最小二乘逆時偏移實現(xiàn)流程

2 模型試算

2.1 層狀模型試算

采用如圖2所示水平層狀介質(zhì)模型(Q均值為101)進(jìn)行了測試。模型縱、橫向采樣點數(shù)為150×150,縱、橫向網(wǎng)格間距為10m。震源為主頻20Hz的雷克子波,共50炮,炮間距30m,炮點分布在網(wǎng)格點上。每炮150個檢波器,采用全接收排列,檢波器間隔10m,采用粘聲介質(zhì)交錯網(wǎng)格正演模擬得到炮記錄。圖3給出了聲波最小二乘逆時偏移第1次迭代結(jié)果、粘聲波最小二乘逆時偏移第1次迭代結(jié)果、聲波最小二乘逆時偏移第40次迭代結(jié)果以及粘聲波最小二乘逆時偏移第40次迭代結(jié)果。

圖2 水平層狀模型

從圖3可以看出,無論是逆時偏移還是最小二乘逆時偏移,采用粘聲波方程的效果明顯好于采用普通聲波方程的效果。由于采樣射線不足,常規(guī)聲波逆時偏移對邊界的刻畫存在一定不足(圖3a,圖3c),但粘聲波逆時偏移的結(jié)果分辨率較高,對淺層噪聲壓制效果較好(圖3b,圖3d)。由圖3c和圖3d可以看出,經(jīng)過40次迭代,聲波最小二乘逆時偏移與粘聲波最小二乘逆時偏移效果都有了較大提高,同相軸能量加強,能量分布更加均衡。但是聲波最小二乘逆時偏移存在較強的淺層噪聲,邊界處成像效果較差,深層能量不足,同相軸不夠清晰。而粘聲波最小二乘逆時偏移能明顯壓制淺層噪聲,成像分辨率更高。特別是在第二層(深層)界面,同相軸更加清晰,分辨率較高,能量更加均衡,有效補償了吸收衰減損失的能量。

圖3 偏移結(jié)果對比a 聲波最小二乘逆時偏移第1次迭代結(jié)果; b 粘聲波最小二乘逆時偏移第1次迭代結(jié)果; c 聲波最小二乘逆時偏移第40次迭代結(jié)果; d 粘聲波最小二乘逆時偏移第40次迭代結(jié)果

圖4給出了最小二乘逆時偏移數(shù)據(jù)殘差隨迭代次數(shù)變化的曲線。可以看出,隨著迭代次數(shù)的增加,殘差穩(wěn)定下降,直至最后收斂。粘聲最小二乘逆時偏移由于對模型粘性吸收衰減進(jìn)行了補償,收斂得更快。圖5,圖6,圖7分別對比了x=250m,750m,1250m處粘聲波最小二乘逆時偏移第1次和第40次迭代結(jié)果、聲波最小二乘逆時偏移第1次和第40次迭代結(jié)果與真實反射系數(shù)。由于觀測系統(tǒng)分布及覆蓋次數(shù)等問題,不同位置能量可能存在一定差異,圖5和圖7能量較弱,圖6能量較強。對比圖5,圖6和圖7 可以看出,第1次迭代采用聲波和粘聲波方程偏移的結(jié)果均與真實反射系數(shù)差距較大,但是粘聲波逆時偏移明顯能量更強,效果更好,更加接近真實反射系數(shù)。經(jīng)過40次迭代以后,聲波和粘聲波偏移效果都有提高,但是很明顯,粘聲波最小二乘逆時偏移已經(jīng)接近真實反射系數(shù),能量更加均衡,尤其是第二層(深層)界面,粘聲波最小二乘逆時偏移能與真實反射系數(shù)接近,保幅性較好。說明該方法能很好地消除地層吸收等因素對地震波振幅及相位的影響,有效提高成像分辨率,實現(xiàn)保幅成像。

圖4 最小二乘逆時偏移數(shù)據(jù)殘差與迭代次數(shù)關(guān)系曲線

2.2 Marmousi模型試算

采用國際通用的Marmousi模型進(jìn)行了試算,圖8 是Marmousi模型的速度場和Q值模型,該模型含有較多的斷層、背斜等,在速度為2000～4000m/s范圍內(nèi)存在較強的衰減。模型縱、橫向采樣點數(shù)為737×380,縱、橫向網(wǎng)格間距為10m。震源為主頻20Hz的雷克子波,共184炮,炮間距為40m,炮點分布在網(wǎng)格點上。每炮737個接收點,接收點間距10m,采用粘聲介質(zhì)交錯網(wǎng)格正演模擬得到炮記錄,時間長度為6.0s,采樣間隔為1.0ms。

圖5 x=250m處成像結(jié)果單道對比a 聲波最小二乘逆時偏移第1次迭代結(jié)果與真實反射系數(shù); b 粘聲波最小二乘逆時偏移第1次迭代結(jié)果與真實反射系數(shù); c 聲波最小二乘逆時偏移第40次迭代結(jié)果與真實反射系數(shù); d 粘聲波最小二乘逆時偏移第40次迭代結(jié)果與真實反射系數(shù)

圖6 x=750m處成像結(jié)果單道對比a 聲波最小二乘逆時偏移第1次迭代結(jié)果與真實反射系數(shù); b 粘聲波最小二乘逆時偏移第1次迭代結(jié)果與真實反射系數(shù); c 聲波最小二乘逆時偏移第40次迭代結(jié)果與真實反射系數(shù); d 粘聲波最小二乘逆時偏移第40次迭代結(jié)果與真實反射系數(shù)

圖7 x=1250m處成像結(jié)果單道對比a 聲波最小二乘逆時偏移第1次迭代結(jié)果與真實反射系數(shù); b 粘聲波最小二乘逆時偏移第1次迭代結(jié)果與真實反射系數(shù); c 聲波最小二乘逆時偏移第40次迭代結(jié)果與真實反射系數(shù); d 粘聲波最小二乘逆時偏移第40次迭代結(jié)果與真實反射系數(shù)

對聲波最小二乘逆時偏移方法和粘聲波最小二乘逆時偏移方法進(jìn)行了對比測試。圖9a為聲波最小二乘逆時偏移第1次迭代結(jié)果,圖9b為粘聲波最小二乘逆時偏移第1次迭代結(jié)果,可以看出,無論是聲波最小二乘逆時偏移,還是粘聲波最小二乘逆時偏移,第1次迭代結(jié)果都存在較強的淺層噪聲,成像結(jié)果模糊不清,分辨率較低,但是在模型底部區(qū)域,粘聲波最小二乘逆時偏移的能量比聲波最小二乘逆時偏移強,同相軸更加明顯,說明該方法具有一定的能量補償作用。圖9c與圖9d分別是聲波最小二乘逆時偏移與粘聲波最小二乘逆時偏移第50次迭代后的結(jié)果?？梢钥闯?隨著迭代次數(shù)的增加,兩種方法都能有效地壓制低頻噪聲,成像結(jié)果分辨率更高,同相軸更加清晰。但在高衰減值區(qū)域,粘聲波最小二乘逆時偏移結(jié)果的能量明顯強于聲波最小二乘偏移結(jié)果,層位更加清晰,同相軸更加明顯,成像位置準(zhǔn)確,構(gòu)造刻畫較好,能量更加均衡,保幅性更好,說明該方法具有較好的振幅衰減補償作用,能夠有效補償因為地下介質(zhì)吸收衰減而造成的能量損失,避免成像位置錯亂等問題。圖10對比了第301道數(shù)據(jù)偏移結(jié)果,可以看出,本文方法偏移的結(jié)果更加接近于真實反射系數(shù),中深層能量對比更加均衡。

圖8 Marmousi模型參數(shù)a 速度場; b 品質(zhì)因子模型

圖9 Marmousi模型偏移結(jié)果a 聲波最小二乘逆時偏移第1次迭代結(jié)果; b 粘聲波最小二乘逆時偏移第1次迭代結(jié)果; c 聲波最小二乘逆時偏移第50次迭代結(jié)果; d 粘聲波最小二乘逆時偏移第50次迭代結(jié)果

圖10 Marmousi模型第50次迭代單道對比(第301道)a 聲波/粘聲波最小二乘逆時偏移第1次迭代結(jié)果; b 聲波/粘聲波最小二乘逆時偏移第50次迭代結(jié)果

3 結(jié)束語

本文基于GSLS波動方程,在反演理論框架下,引入粘聲介質(zhì)最小二乘理論,推導(dǎo)了粘聲介質(zhì)最小二乘逆時偏移中的反偏移算子,針對偏移中的不穩(wěn)定現(xiàn)象,構(gòu)建了對應(yīng)的粘聲逆時偏移伴隨算子,建立了粘聲介質(zhì)最小二乘逆時偏移的目標(biāo)函數(shù),給出了最小二乘逆時偏移的實現(xiàn)流程與步驟；通過層狀模型和Marmousi模型試算,驗證了粘聲介質(zhì)最小二乘逆時偏移方法的正確性與有效性。研究結(jié)果表明：該方法可以合理消除地層吸收、透射及幾何擴散等因素對地震振幅、頻率及相位的影響,提高成像分辨率；消除常規(guī)逆時偏移方法由于觀測系統(tǒng)對地下一些地區(qū)照明弱導(dǎo)致的能量不均衡對成像效果的影響,以及震源波場和接收點波場互相關(guān)引起的噪聲,改善復(fù)雜介質(zhì)條件下地質(zhì)體成像的振幅屬性,實現(xiàn)保幅、保真、高精度成像。

但是,實際生產(chǎn)中地震資料質(zhì)量較差,噪聲較多,如何在低信噪比情況下取得良好的成像效果將是我們下一步研究的課題。而且實際資料中通常會存在傾角約束、井約束等先驗信息,如何利用先驗信息降低對背景速度的依賴也是我們下一步努力的方向。

參 考 文 獻(xiàn)

[1] RIBODETTI A,VIRIEUX J.Asymptotic theory for imaging the attenuation factors QP and QS:lecture notes in physics[M].Berlin,Heidelberg:Springer,1997:334-353

[2] CAUSSE E,URSIN B.Viscoacoustic reverse-time migration[J].Journal of Seismic Exploration,2000,9(2):165-184

[3] ZHU T,HARRIS J M,BIONDI B.Q-compensated reverse-time migration[J].Geophysics,2014,79(3): S77-S87

[4] ZHANG Y,ZHANG P,ZHANG H.Compensating for visco-acoustic effects in reverse-time migration[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010:3160-3164

[5] TARANTOLA A.Linearized inversion of seismic reflection data[J].Geophysical Prospecting,1984,32(6):998-1015

[6] LEBRAS R,CLAYTON R W.An iterative inversion of back-scattered acoustic waves[J].Geophysics,1988,53(4):501-508

[7] LAMBRARE G,VIRIEUX J,MADARIAGA R,et al.Iterative asymptotic inversion in the acoustic approximation[J].Geophysics,1992,57(9):1138-1154

[8] NEMETH T,WU C,SCHUSTER G T.Least-squares migration of incomplete reflection data[J].Geophysics,1999,64(1):208-221

[9] LI Q Y,HUANG J P,LI Z C.Cross-correlation least-squares reverse time migration in the pseudo-time domain[J].Journal of Geophysics and Engineering,2017,14(4):841-851

[10] YAO G,JAKUBOWICZ H.Least-squares reverse-time migration in a matrix-based formulation[J].Geophysical Prospecting,2016,64(3):611-621

[11] QU Y M,LI Z C,HUANG J P,et al.Elastic full waveform inversion for surface topography[J].Geophysics,2017,82(5):R269-R285

[12] QU Y M,LI Z C,HUANG J P,et al.Pre-stack reverse time migration based on layered mapping method for irregular surface[J].Extended Abstract of 77thEAGE Conference & Exhibition,2015:234-237

[13] 郭書娟,馬方正,段心標(biāo),等.最小二乘逆時偏移成像方法的實現(xiàn)與應(yīng)用研究[J].石油物探,2015,54(3):301-308

GUO S J,MA F Z,DUAN X B,et al.Research of least-squares reverse-time migration imaging method and its application[J].Geophysical Prospecting for Petroleum,2015,54(3):301-308

[14] 李振春,周麗穎,黃建平,等.最小二乘逆時偏移在碳酸鹽巖縫洞型儲層成像中的應(yīng)用[J].地球物理學(xué)進(jìn)展,2017,32(2):664-671

LI Z C,ZHOU L Y,HUANG J P,et al,Application of least square reverse time migration in the imaging of fractured-type carbonate reservoirs[J].Progress in Geophysics,2017,32(2):664-671

[15] QU Y M,LI J L,HUANG J P.Elastic least-squares reverse time migration with velocities and density perturbation[J].Geophysical Journal International,2017,212(2):1033-1056

[16] QU Y M,HUANG J P,LI Z C,et al.Attenuation compensation in anisotropic least-squares reverse time migration[J].Geophysics,2017,82(6):S411-S423

[17] 黃建平,李振春,孔雪,等.碳酸鹽巖裂縫型儲層最小二乘偏移成像方法研究[J].地球物理學(xué)報,2013,56(5):1716-1725

HUANG J P,LI Z C,KONG X,et al.A study of least-squares migration imaging method for fractured-type carbonate reservoir[J].Chinese Journal of Geophysics,2013,56(5):1716-1725

[18] 劉玉金,李振春,吳丹,等.局部傾角約束最小二乘偏移方法研究[J].地球物理學(xué)報,2013,56(3):1003-1011

LIU Y J,LI Z C,WU D,et al.The research on local slope constrained least-squares migration[J].Chinese Journal of Geophysics,2013,56(3):1003-1011

[19] QU Y M,HUANG J P,LI Z C,et al,Viscoacoustic anisotropic full waveform inversion[J].Journal of Applied Geophysics,2017,136(1):484-497