姚振岸 孫成禹 喻志超 馬 振

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

0 引言

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

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

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

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

1 理論方法

1.1 點(diǎn)擴(kuò)散函數(shù)與去模糊濾波器

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

d=Lm0

(1)

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

(2)

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

(3)

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

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

mmig-ref=LTLmref=LTdref

(4)

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

mmig-ref(x0,z0)idV

(5)

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

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

(6)

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

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

由卷積的交換律可得

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

(7)

將上述卷積運(yùn)算寫(xiě)成矩陣運(yùn)算形式，參考偏移結(jié)果[mmin-ref]i變成一個(gè)(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個(gè)局部窗口內(nèi)近似等于(LTL)-1。在實(shí)際應(yīng)用中，為了保證計(jì)算效果，建議局部窗口長(zhǎng)度不超過(guò)1.5倍地震波長(zhǎng)。

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

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

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

(10)

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

(11)

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

(12)

(13)

(14)

(15)

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

(16)

1.3 黏聲最小二乘逆時(shí)偏移

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

(17)

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

(18)

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

(19)

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

δP(xr,t;xs)

(20)

式中“*”表示時(shí)間方向的褶積運(yùn)算。

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

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í)際觀測(cè)數(shù)據(jù)和正演數(shù)據(jù)之間的殘差。

成像擾動(dòng)δm與體積模量擾動(dòng)δK近似線性相關(guān)，而δK可由背景場(chǎng)(通過(guò)式(19)計(jì)算)和伴隨場(chǎng)(通過(guò)式(22)計(jì)算)之間的零時(shí)間延遲互相關(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最小二乘逆時(shí)偏移

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

(26)

(27)

(28)

式中α是更新步長(zhǎng)。

2 Marmousi Ⅱ模型測(cè)試

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

圖2 Marmousi Ⅱ速度及Q模型

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

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

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

圖5 圖3、圖4藍(lán)色矩形框局部放大顯示

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

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

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

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

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

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

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

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

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

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

3 結(jié)論與認(rèn)識(shí)

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

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

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

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

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