姚振岸 孫成禹 喻志超 馬 振

(①東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點學(xué)科實驗室，江西南昌 330013；②中國石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580；③北京大學(xué)地球與空間科學(xué)學(xué)院石油與天然氣研究中心，北京 100871)

0 引言

儲層孔隙中含有油氣水等流體，會引起地震波傳播過程中的強振幅衰減和相位畸變[1]，從而導(dǎo)致傳統(tǒng)的聲波或者彈性波偏移方法分辨率降低或者失效。關(guān)于地震波在黏彈性介質(zhì)中的傳播機理，主要有Maxwell模型、Kelvin-Vogit模型、標(biāo)準(zhǔn)線性固體、廣義標(biāo)準(zhǔn)線性體、分數(shù)階常Q模型等[2-6]。早期的地震資料衰減補償主要是使用反Q濾波在數(shù)據(jù)域?qū)崿F(xiàn)[7-10]，如Blanch等[11-12]提出用τ方法模擬常Q，并基于伴隨狀態(tài)法在τ-p域中用線性反演方法恢復(fù)黏聲介質(zhì)的短波長分量。考慮到地震波衰減是在傳播過程中發(fā)生的，融合衰減補償?shù)钠品椒ㄔ絹碓绞艿街匾?。Xin等[13]和Xie等[14]利用射線追蹤的方法在疊前深度偏移過程中進行衰減補償；白敏等[15]和代福材等[16]在高斯束疊前深度偏移過程中補償振幅衰減和校正相位畸變；Dai等[17]、Yu等[18]、Wang[19]和Valenciano等[20]在頻率域用單程波動方程偏移進行衰減補償；陳志德等[21]引入等效Q值實現(xiàn)頻率域波場延拓，基于單程波方程和穩(wěn)定相點原理實現(xiàn)黏聲疊前時間偏移，并提高了補償算法的穩(wěn)定性。朱峰等[22]推導(dǎo)了黏聲VTI介質(zhì)高階傅里葉有限差分波場延拓算子，實現(xiàn)了適應(yīng)高角度入射波場的黏聲疊前深度偏移方法。

不同學(xué)者提出了多種黏聲波動方程，在逆時偏移過程中分別對振幅損失和相位畸變進行補償和校正[23-27]。李振春等[28]基于廣義標(biāo)準(zhǔn)線性體模型對黏聲介質(zhì)波動方程進行線性化，并借助伴隨狀態(tài)法實現(xiàn)了黏聲介質(zhì)最小二乘逆時偏移(Q-LSRTM)的迭代算法，獲得了較高的成像分辨率。Dutta等[29]、Dai等[30]基于一階松弛機制時間域黏聲波動方程，用迭代的Q-LSRTM有效補償和校正了強衰減層導(dǎo)致的振幅損失和相位畸變。徐凱等[31]在反演框架下，構(gòu)建了Q-LSRTM多約束正則化目標(biāo)函數(shù)，并采用變分方法求解，消除了震源波場與接收點波場互相關(guān)引起的噪聲，提高了復(fù)雜地質(zhì)體成像的精度、增強了保幅性。相比于傳統(tǒng)聲波LSRTM，Q-LSRTM能實現(xiàn)反射體的精確歸位且振幅更均衡，收斂速率更快[32]。但Q-LSRTM也存在一定問題： 其一，由于在傳統(tǒng)Q-LSRTM中用于反傳數(shù)據(jù)殘差的伴隨Q傳播算子也是衰減的，導(dǎo)致成像結(jié)果分辨率降低；其二，基于時間域黏聲波動的LSRTM計算成本較高，尤其是需要很多次迭代才能達到預(yù)期成像效果。在完全彈性情況下，為了提高傳統(tǒng)LSRTM的收斂速度，李闖等[33]發(fā)展了一種預(yù)條件偏移算法，提高了鹽下構(gòu)造成像的分辨率和保幅性。

在實際野外地震勘探中，由于稀疏的震源和檢波器，記錄偏移孔徑有限、速度變化快及構(gòu)造復(fù)雜，以及不均勻的地震照明等原因，偏移剖面當(dāng)中往往存在多種假象[34]。為了消除偏移假象、提升成像分辨率，偏移反褶積被廣泛采用[35]。偏移成像本質(zhì)上就是真實反射率的一種模糊，模糊算子即為Hessian算子，而偏移反褶積處理就等效于一種近似的Hessian逆處理。Hu等[34-35]在空間—波數(shù)域設(shè)計了一個偏移反褶積算子壓制疊后偏移假象；Yu等[36]將偏移反褶積從疊后偏移拓展到疊前深度偏移。Guitton[37]直接在空間域構(gòu)建了一系列匹配濾波器近似Hessian矩陣的逆實現(xiàn)偏移反褶積，后來這些匹配濾波器被Aoki等[38]作為預(yù)條件化因子用于傳統(tǒng)的LSRTM，獲得了良好的成像效果，這些匹配濾波器也被稱為去模糊濾波器。針對多震源LSRTM，Dai等[39-40]用去模糊濾波器減少串?dāng)_噪聲并且加速多震源LSRTM的收斂。但目前偏移反褶積的原理和實施都是基于無耗散介質(zhì)假設(shè)，而地下儲層往往具有衰減性質(zhì)。

為了克服Q-LSRTM的缺陷，本文將匹配濾波偏移反褶積與Q-LSRTM相結(jié)合。首先基于點擴散函數(shù)(PSF)構(gòu)建黏聲去模糊濾波器，然后將其作為Q-LSRTM迭代過程中的預(yù)條件化因子，最終實現(xiàn)高分辨的衰減補償偏移。相比于Q-LSRTM，預(yù)條件Q-LSRTM成像結(jié)果分辨率更高，反射體振幅更均衡，收斂速度更快。就計算成本來說，盡管預(yù)條件Q-LSRTM需要預(yù)先計算黏聲去模糊濾波器，但其達到預(yù)期成像效果需要的迭代次數(shù)遠比Q-LSRTM少。

1 理論方法

1.1 點擴散函數(shù)與去模糊濾波器

在聲學(xué)介質(zhì)中，基于Born近似，觀測地震記錄d的形成過程可以表示為

d=Lm0

(1)

式中：L表示線性模擬算子，與觀測系統(tǒng)、震源子波和地下介質(zhì)模型參數(shù)有關(guān)；m0是地下反射率模型。對應(yīng)L，可以獲得其伴隨偏移算子LT，則偏移結(jié)果mmig可表示為

(2)

mmig其實是真實反射率模型m0的模糊版本，LTL就是一個模糊算子，即偏移格林函數(shù)[41]，也被稱為點擴散函數(shù)(PSF)。由于震源子波、觀測系統(tǒng)等多種因素的影響，偏移成像中往往存在較多的假象。為了獲得更清晰的成像結(jié)果，可以取LTL的逆并作用于成像結(jié)果

(3)

但在地震成像時，直接計算(LTL)-1是不現(xiàn)實的，為此采用迭代LSRTM方法予以求解，但其計算量往往比標(biāo)準(zhǔn)偏移成像大一個數(shù)量級。為了加速收斂，設(shè)計去模糊算子的近似逆算子(LTL)-1，并在LSRTM中作為預(yù)條件化算子[35-39]。

實際應(yīng)用中，去模糊濾波器基于背景模型即偏移速度場和參考反射率模型獲得，以均勻介質(zhì)速度背景場和均勻分布點散射體反射率模型為例闡述去模糊濾波器的獲取過程。圖1a為均勻分布點散射體模型，圖1b為對應(yīng)的疊前深度偏移成像結(jié)果，其中震源和檢波器均勻分布于模型上方。單個散射體的PSF響應(yīng)即為其偏移結(jié)果。記參考反射率模型為mref，則在Bron近似下利用式(1)可獲得正演模擬數(shù)據(jù)dref，繼而可以獲得標(biāo)準(zhǔn)的參考偏移結(jié)果

mmig-ref=LTLmref=LTdref

(4)

將點散射體模型分解成多個子區(qū)域(如圖1a中黑色矩形框)，使點散射體位于矩形窗口的中央，窗口大小為wx×wz，足以覆蓋其PSF響應(yīng)的有效能量(圖1b中矩形框)。在這個局部窗口內(nèi)，假設(shè)PSF保持不變。在wx×wz局部窗口區(qū)域內(nèi)，去模糊濾波器將參考偏移成像結(jié)果和參考反射率模型聯(lián)系起來[38，40]，即

mmig-ref(x0,z0)idV

(5)

式中：i為局部窗口序號；積分限Vi為整個局部窗口區(qū)域。式(5)寫成矩陣形式為

[mref]i=Fi?[mmig-ref]i

(6)

式中：“?”表示空間卷積；[mref]i和[mmig-ref]i分別為第i個局部窗口內(nèi)參考反射率模型和參考偏移成像結(jié)果。與PSF相同，在每個局部窗口內(nèi)去模糊濾波器保持不變，定義濾波器的空間尺寸為fx×fz。

圖1 點散射模型及其偏移成像轉(zhuǎn)換示意圖

由卷積的交換律可得

Fi?[mmig-ref]i=[mmig-ref]i?Fi

(7)

將上述卷積運算寫成矩陣運算形式，參考偏移結(jié)果[mmin-ref]i變成一個(N+M-1)×N階卷積矩陣，去模糊濾波器Fi變?yōu)镹×1階向量f，其中N=fx×fz且M=wx×wz，即為

[mmig-ref]ifi=[mref]i

(8)

(9)

繼而采用三角(LU)分解方法求解式(9)即可獲得去模糊濾波器，在第i個局部窗口內(nèi)近似等于(LTL)-1。在實際應(yīng)用中，為了保證計算效果，建議局部窗口長度不超過1.5倍地震波長。

1.2 黏聲介質(zhì)中的去模糊濾波器

將去模糊濾波器應(yīng)用到偏移成像的過程就是一種偏移反褶積處理。前人對偏移反褶積的研究都集中在無耗散介質(zhì)當(dāng)中，但強吸收會對地震波振幅和相位產(chǎn)生明顯的衰減和畸變[1]。為了補償衰減效應(yīng)，Dutta等[29]發(fā)展了時間域Q-LSRTM技術(shù)，成像結(jié)果歸位更精確，振幅更均衡。在傳統(tǒng)Q-LSRTM中，用于反傳數(shù)據(jù)殘差的伴隨Q傳播算子也是衰減的，因此成像分辨率有一定的損失。

這種分辨率的損失能通過黏聲介質(zhì)中的PSF即偏移格林函數(shù)予以解釋。給定速度為v0的均勻介質(zhì)和位于xs處角頻率為ω的簡諧點源，聲學(xué)格林函數(shù)可以表示為

(10)

如果介質(zhì)是耗散的，則將完全彈性介質(zhì)中的實速度v0替換成復(fù)速度[1]

(11)

從而得到黏聲介質(zhì)格林函數(shù)

(12)

(13)

(14)

(15)

如1.1所述，在偏移成像的不同子區(qū)域內(nèi)，通過匹配參考成像數(shù)據(jù)和參考反射率模型，得到局部去模糊濾波器，即

(16)

1.3 黏聲最小二乘逆時偏移

基于標(biāo)準(zhǔn)線性固體(SLS)黏彈性模型一階松弛機制，二維黏聲介質(zhì)波動程可寫為[11，42]

(17)

式中：ρ為密度；P為壓力；K為體積模量；S為震源子波；γp為記憶變量；v為質(zhì)點速度向量；τ與品質(zhì)因子Q有關(guān)，可表示為

(18)

其中τσ和τε分別表示應(yīng)力和應(yīng)變松弛時間。在Born近似的框架下，對體積模量K做一定擾動δK，則擾動場可以寫為

(19)

式中忽略了二階擾動變量。若用壓力和記憶變量格林函數(shù)GP(xr,t;x0,0)和Gγp(xr,t;x0,0)，式(19)也可表示為

δP(xr,t;xs)

(20)

式中“*”表示時間方向的褶積運算。

在Q-LSRTM的框架下，式(20)的解等效于矩陣—向量運算

dQ=LQm0

(21)

式中：dQ為Born模擬的黏聲地震數(shù)據(jù)；LQ為線性黏聲模擬算子；m0為地下介質(zhì)反射率。用伴隨狀態(tài)法可推導(dǎo)式(19)的伴隨方程[11，12，43]為

(22)

式中： (q，u，s)為(P，v，γp)的伴隨狀態(tài)變量；Δd表示每次迭代中實際觀測數(shù)據(jù)和正演數(shù)據(jù)之間的殘差。

成像擾動δm與體積模量擾動δK近似線性相關(guān)，而δK可由背景場(通過式(19)計算)和伴隨場(通過式(22)計算)之間的零時間延遲互相關(guān)得到，即

(23)

上式可用格林函數(shù)表示為

[GP(xr,-t;x,0)*ΔP(xr,t;xs)]+

[GτP(xr,-t;x,0)*ΔP(xr,t;xs)]}

(24)

用矩陣可以表示為

(25)

1.4 預(yù)條件Q最小二乘逆時偏移

Q-LSRTM的誤差泛函[29]為

(26)

(27)

(28)

式中α是更新步長。

2 Marmousi Ⅱ模型測試

為了說明本文提出的預(yù)條件Q-LSRTM的有效性，應(yīng)用Marmousi Ⅱ模型進行測試。圖2a和圖2b分別為正演速度模型和Q值模型，其中Q只取兩個值，Q=20表示衰減層，Q=10000可視為彈性介質(zhì)。整個模型的網(wǎng)格點數(shù)為801×351，模型縱橫向采樣間隔均為10m。在合成黏聲地震記錄時，炮間隔為50m，共計150炮均勻分布在模型表面，道間距為10m，震源采用主頻為15Hz的Ricker子波。采集系統(tǒng)位于地表，且排列長度固定為整個模型的長度。計算中采用時間二階精度、空間八階精度的時間域交錯網(wǎng)格有限差分法及C-PML邊界條件。本文共用到兩套數(shù)據(jù)，一套是基于聲波方程的正演結(jié)果，另一套是基于一階松弛機制黏聲波動方程的正演結(jié)果。對真實速度模型和Q值模型進行二維平滑，得到偏移速度模型(圖2c)和偏移Q值模型(圖2d)。分別執(zhí)行聲波LSRTM、Q-LSRTM、預(yù)條件Q-LSRTM三種最小二乘逆時偏移，迭代次數(shù)均為20次，其中第1次迭代結(jié)果即為聲波逆時偏移(RTM)、Q-RTM和預(yù)條件Q-RTM的成像結(jié)果。圖3為聲波數(shù)據(jù)聲波RTM和LSRTM成像結(jié)果，并以此作為基準(zhǔn)對比數(shù)據(jù)。

圖2 Marmousi Ⅱ速度及Q模型

圖3 聲波數(shù)據(jù)聲波RTM(a)和LSRTM(b)成像結(jié)果

圖4 黏聲數(shù)據(jù)不同方法成像結(jié)果對比

圖5、圖6分別為聲波數(shù)據(jù)LSRTM(圖3b)、黏聲數(shù)據(jù)LSRTM(圖4b)、Q-LSRTM(圖4d)、預(yù)條件Q-LSRTM(圖4f)成像結(jié)果藍色、紅色矩形框局部放大顯示，可以看出：預(yù)條件Q-LSRTM提高了黏聲數(shù)據(jù)的成像精度和分辨率(藍色箭頭所指)；采用LSRTM對黏聲數(shù)據(jù)進行成像時，不僅存在成像分辨率低的問題，而且由于介質(zhì)黏彈性導(dǎo)致地震波振幅和相位的畸變，還會導(dǎo)致深層反射體歸位出現(xiàn)錯亂(圖6b紅色箭頭所指)。相對而言，盡管Q-LSRTM能有效歸位主要反射界面，但分辨率較低，且不能呈現(xiàn)構(gòu)造細節(jié)，預(yù)條件Q-LSRTM成像與基準(zhǔn)聲波LSRTM結(jié)果比較接近，分辨率較高。

圖7為聲波數(shù)據(jù)LSRTM(圖5a)、黏聲數(shù)據(jù)LSRTM(圖5b)、Q-LSRTM(圖5c)、預(yù)條件Q-LSRTM(圖5d)在x=5.0km處成像結(jié)果的中深層歸一化波數(shù)譜對比，可以看出，傳統(tǒng)LSRTM和Q-LSRTM的主波數(shù)都相對偏低，只有預(yù)條件Q-LSRTM的波數(shù)譜與基準(zhǔn)聲波LSRTM接近，也說明預(yù)條件Q-LSRTM提高了黏聲數(shù)據(jù)的成像分辨率。

圖5 圖3、圖4藍色矩形框局部放大顯示

圖6 圖3、圖4紅色矩形框局部放大顯示

圖8為聲波數(shù)據(jù)LSRTM(圖6a)、黏聲數(shù)據(jù)LSRTM(圖6b)、Q-LSRTM(圖6c)、預(yù)條件Q-LSRTM(圖6d)在x=2.6km處成像結(jié)果的深層歸一化波數(shù)譜對比，可以看出，黏聲數(shù)據(jù)的聲波LSRTM波數(shù)譜嚴重偏離基準(zhǔn)數(shù)據(jù)譜，Q-LSRTM有所改善，預(yù)條件Q-LSRTM能恢復(fù)出高波數(shù)細節(jié)信息，最接近于基準(zhǔn)數(shù)據(jù)。

圖7 x=5.0km處不同方法成像結(jié)果的中深層波數(shù)譜對比

圖8 x=2.6km處不同方法成像結(jié)果的深層波數(shù)譜對比

圖9為黏聲數(shù)據(jù)聲波LSRTM、Q-LSRTM、預(yù)條件Q-LSRTM(預(yù)條件20次)、預(yù)條件Q-LSRTM(預(yù)條件6次)四種情況下數(shù)據(jù)殘差隨著迭代次數(shù)的變化曲線，可以看出，四種偏移成像結(jié)果的精度是逐漸升高的。就收斂速率來看，預(yù)條件Q-LSRTM第4次迭代和第8次迭代時的數(shù)據(jù)殘差接近于Q-LSRTM第8次和第20次迭代的數(shù)據(jù)殘差，因此預(yù)條件Q-LSRTM能減少約50%的計算量，尤其在早期的幾次迭代當(dāng)中體現(xiàn)的更為明顯。預(yù)條件Q-LSRTM收斂速率提升的原因在于黏聲去模糊濾波器是Hessian算子逆的一種有效近似估計，盡管需要預(yù)先計算去模糊濾波器，但計算成本還是要比Q-LSRTM低得多。兼顧到計算精度和計算效率，參考反射率模型中點散射體的布置應(yīng)該相對稀疏，由此會引入一些全局噪聲，如果預(yù)條件在第6次迭代后終止，數(shù)據(jù)殘差會進一步減小，相應(yīng)地成像精度也有所提高。

為了測試預(yù)條件Q-LSRTM對Q值的敏感性，分別應(yīng)用Q=20、40、80、100四種模型進行偏移成像，結(jié)果如圖4f、圖10所示，圖11為對應(yīng)的參考反射率模型的黏聲PSF，可以看出：隨著Q值的增大，衰減層及其下部的反射體的成像振幅逐漸減弱，成像分辨率逐漸降低；相反地，當(dāng)Q值越大時，衰減越弱，參考反射率模型中的點散射體偏移響應(yīng)即PSF越強，從而用其構(gòu)建黏聲去模糊濾波器時，考慮到的衰減效應(yīng)越弱、成像精度越低。

圖9 數(shù)據(jù)殘差隨著迭代次數(shù)的變化曲線

圖10 不同Q值偏移模型的預(yù)條件Q-LSRTM偏移結(jié)果

因此，一個相對精確的偏移Q值模型是預(yù)條件Q-LSRTM取得良好成像效果的前提。

圖11 不同Q值偏移模型的參考反射率模型的黏聲PSF

3 結(jié)論與認識

(1)在傳統(tǒng)Q-LSRTM成像當(dāng)中，數(shù)據(jù)殘差的反傳過程也是衰減的，會導(dǎo)致成像分辨率降低；而基于PSF的預(yù)條件Q-LSRTM方法采用黏聲去模糊濾波器補償由于地下強衰減導(dǎo)致的振幅和分辨率降低，能得到較高的成像精度和更均衡的成像振幅。

(2)由于地下模型是未知的，構(gòu)建黏聲去模糊濾波器時一般選用均勻分布的點散射體模型?；谄扑俣葓龊蚎值模型，采用和實際一致的觀測系統(tǒng)和地震子波，通過正演模擬和偏移成像獲得黏聲PSF，最終采用匹配濾波的思想構(gòu)建黏聲去模糊濾波器。兼顧到計算精度和計算效率，參考反射率模型中點散射體密度和局部窗口長度需根據(jù)構(gòu)造復(fù)雜度確定。

(3)與Q-LSRTM相比，預(yù)條件Q-LSRTM的收斂效率得到大幅度提升，在前幾次迭代中更為明顯。為了獲得更高的成像精度和計算效率，預(yù)條件Q-LSRTM只在前幾次迭代中使用為宜。就計算成本而言，盡管預(yù)條件Q-LSRTM需要預(yù)先計算黏聲去模糊濾波器，但計算耗時仍然比Q-LSRTM小得多。

(4)與Q-LSRTM類似，為了獲得較好的成像效果，預(yù)條件Q-LSRTM也需要一個相對精確的Q值模型估計，同時也必須有一個相對精確的偏移速度模型。

(5)本文的方法可拓展到黏彈各向異性等更復(fù)雜介質(zhì)的偏移成像當(dāng)中，考慮到計算量，建議采用并行加速計算。