融合點擴散函數(shù)的預條件黏聲最小二乘逆時偏移

姚振岸 孫成禹 喻志超 馬 振

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

0 引言

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

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

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

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

1 理論方法

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

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

d=Lm0

(1)

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

(2)

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

(3)

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

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

mmig-ref=LTLmref=LTdref

(4)

將點散射體模型分解成多個子區(qū)域(如圖1a中黑色矩形框)，使點散射體位于矩形窗口的中央，窗口大小為wx×wz，足以覆蓋其PSF響應的有效能量(圖1b中矩形框)。在這個局部窗口內，假設PSF保持不變。在wx×wz局部窗口區(qū)域內，去模糊濾波器將參考偏移成像結果和參考反射率模型聯(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個局部窗口內參考反射率模型和參考偏移成像結果。與PSF相同，在每個局部窗口內去模糊濾波器保持不變，定義濾波器的空間尺寸為fx×fz。

圖1 點散射模型及其偏移成像轉換示意圖

由卷積的交換律可得

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

(7)

將上述卷積運算寫成矩陣運算形式，參考偏移結果[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個局部窗口內近似等于(LTL)-1。在實際應用中，為了保證計算效果，建議局部窗口長度不超過1.5倍地震波長。

1.2 黏聲介質中的去模糊濾波器

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

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

(10)

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

(11)

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

(12)

(13)

(14)

(15)

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

(16)

1.3 黏聲最小二乘逆時偏移

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

(17)

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

(18)

其中τσ和τε分別表示應力和應變松弛時間。在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為地下介質反射率。用伴隨狀態(tài)法可推導式(19)的伴隨方程[11，12，43]為

(22)

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

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

(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 預條件Q最小二乘逆時偏移

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

(26)

(27)

(28)

式中α是更新步長。

2 Marmousi Ⅱ模型測試

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

圖2 Marmousi Ⅱ速度及Q模型

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

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

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

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

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

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

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

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

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

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

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

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

圖10 不同Q值偏移模型的預條件Q-LSRTM偏移結果

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

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

3 結論與認識

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

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

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

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

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