最小二乘疊前時間偏移在地震數(shù)據(jù)規(guī)則化中的應(yīng)用

吳 丹 龔仁彬 王從鑌 胥小馬 吳海莉

(中國石油勘探開發(fā)研究院西北分院計算機(jī)所，北京 100191)

0 引言

在三維地震勘探中，受地表障礙物、觀測設(shè)備、施工條件等因素的限制，采集到的數(shù)據(jù)往往出現(xiàn)不規(guī)則、缺道、壞道等現(xiàn)象。地震數(shù)據(jù)的這種不規(guī)則性會嚴(yán)重影響后續(xù)的多次波壓制[1]、偏移成像[2]、四維地震監(jiān)測[3]等環(huán)節(jié)。為了消除地震數(shù)據(jù)不規(guī)則的影響，一般采用地震數(shù)據(jù)規(guī)則化方法重建缺失數(shù)據(jù)。但是，由于數(shù)據(jù)規(guī)則化本質(zhì)上是一個欠定的反演問題，因此一直以來都是地震數(shù)據(jù)處理領(lǐng)域十分具有挑戰(zhàn)性的研究課題。

目前，地震數(shù)據(jù)規(guī)則化方法主要分為三類：第一類是在信號變換域引入稀疏約束恢復(fù)地震數(shù)據(jù)，如Radon變換[4]、傅里葉變換[5-8]、小波變換[9]、曲波變換[10-11]等；第二類是基于褶積算子的預(yù)測濾波器方法，根據(jù)已知數(shù)據(jù)提取局部地震屬性(如局部傾角)，再擴(kuò)展該信息以恢復(fù)缺失數(shù)據(jù)[12-18]；第三類是模型驅(qū)動方法，利用已知的地下模型信息，采用積分延拓算子重建缺失數(shù)據(jù)，常見的延拓算子包括炮檢距延拓[19]、炮點延拓[20]、方位角延拓[21]等。其中前兩類方法屬于數(shù)據(jù)驅(qū)動的規(guī)則化方法，

以往基于模型驅(qū)動的數(shù)據(jù)規(guī)則化方法雖然利用了已知的地下模型信息，但是處理精度往往過度依賴于模型的準(zhǔn)確性。當(dāng)?shù)叵履Ｐ痛嬖谝欢ㄕ`差時，規(guī)則化結(jié)果不理想。本文采用疊前時間偏移和反偏移算子作為積分延拓算子，在成像空間保留了一個維度的冗余信息(如地表炮檢距)，使映射到成像空間的過程中不丟失模型的誤差信息； 再通過引入更加穩(wěn)定的預(yù)條件算子，可以在成像空間更好地提取和應(yīng)用反射波同相軸特征，從而改善數(shù)據(jù)規(guī)則化質(zhì)量。最后通過理論模型和實際資料處理驗證本文方法能否夠有效改善數(shù)據(jù)規(guī)則化效果，改進(jìn)最終成像質(zhì)量。

1 方法原理

1.1 最小二乘偏移理論

基于一階Born近似理論，反射地震數(shù)據(jù)d可以通過地下反射系數(shù)模型m線性映射得到

d=Lm

(1)

式中L為線性映射算子，也稱為Born正演算子，可以是疊前/疊后反時間偏移算子，也可以是疊前/疊后反深度偏移算子。該算子與速度模型有關(guān)，由于規(guī)則化過程中速度模型固定，因此式(1)忽略了速度這個參數(shù)。最小二乘偏移最小化目標(biāo)函數(shù)為

(2)

式中：dobs為觀測數(shù)據(jù)；D為正則化算子；μ為阻尼因子。式(2)的顯式解為

m=(LTL+μDTD)-1LTdobs

(3)

定義一個新的向量p，滿足

m=Pp

(4)

則式(2)變?yōu)?/p>

(5)

式中P=D-1為預(yù)條件算子。上式的解為

p=(PTLTLP+μI)-1PTLTdobs

(6)

式中I為單位矩陣。將式(6)代入式(4)，可得

m=P(PTLTLP+μI)-1PTLTdobs

(7)

Clearbout[22]指出，合理的預(yù)條件算子比對應(yīng)的正則化算子可以更快地加速迭代反演的收斂速度。

由于Born正演(或反偏移)算子L一般無法用矩陣顯式表達(dá)，而且模型和數(shù)據(jù)的規(guī)模很大，因此一般采用共軛梯度法近似求解式(7)，得到分辨率更高、振幅保真度更好的反射系數(shù)模型。在相同的理論框架下，選擇不同的反偏移算子，可以構(gòu)建不同的最小二乘偏移算法，如最小二乘疊前時間偏移[23-24]、最小二乘逆時偏移[25-26]等。

1.2 地震數(shù)據(jù)規(guī)則化

當(dāng)觀測數(shù)據(jù)不規(guī)則時，式(5)可以修改為

(8)

式中K為Mask算子，該算子的作用是將缺失地震道置0，非缺失地震道保持不變。通過共軛梯度法最小化目標(biāo)函數(shù)式(8)得到最優(yōu)解popt，再應(yīng)用預(yù)條件算子和反偏移算子，則可以重建缺失地震數(shù)據(jù)

dreg=LPpopt

(9)

重建后的地震數(shù)據(jù)可以表示為

dreg=LP(PTLTKTKLP+μI)-1PTLTKTdobs

(10)

分三種情況分析式(10)中各項對數(shù)據(jù)規(guī)則化的影響。

第一種情況，忽略Hessian矩陣求逆過程(式(10)中的求逆項)，即直接采用偏移+反偏移進(jìn)行數(shù)據(jù)重建，式(10)變?yōu)?/p>

dreg=LLTKTdobs

(11)

這種做法雖然也能夠起到數(shù)據(jù)插值的作用，但是振幅信息不可靠。通過設(shè)計一個匹配濾波器可以在一定程度上解決這個問題，但是匹配濾波器的估計本身也是一個欠定的反問題，存在多解性。

第二種情況，忽略預(yù)條件項，即直接采用最小二乘偏移+反偏移進(jìn)行數(shù)據(jù)重建，式(10)變?yōu)?/p>

dreg=L(LTKTKL)-1LTKTdobs

(12)

這種做法可以解決偏移+反偏移中存在的振幅匹配問題，但是無法消除數(shù)據(jù)缺失帶來的成像噪聲，該噪聲反投影到數(shù)據(jù)域，會降低插值數(shù)據(jù)的信噪比。

第三種情況，完整地求解式(10)，即采用預(yù)條件最小二乘偏移+反偏移進(jìn)行數(shù)據(jù)重建，這種做法既考慮了振幅不匹配的問題，又可以壓制數(shù)據(jù)不規(guī)則帶來的噪聲。

實際上，預(yù)條件算子的選取對數(shù)據(jù)規(guī)則化效果的影響很大。本文采用因果積分算子作為預(yù)條件算子，該算子的作用是使振幅隨炮檢距平滑變化。當(dāng)速度比較準(zhǔn)確時，共成像點道集拉平，滿足振幅平滑變化的條件；當(dāng)速度誤差較大時，相鄰道之間的振幅差異也不會太大，因此也基本能滿足條件；當(dāng)速度誤差很大時，需要進(jìn)一步采用更加合適的預(yù)條件算子，以保留原始數(shù)據(jù)中的剩余旅行時信息。

1.3 疊前時間偏移/反偏移

綜合考慮算子的靈活性、計算效率以及預(yù)處理階段對規(guī)則化的需求，本文采用Kirchhoff疊前時間反偏移作為線性正演算子L，該算子可以表示為

(13)

式中：m(x，h，τ)為時間偏移共炮檢距道集；d(y，H，t)為疊前共炮檢距數(shù)據(jù)；x和h分別為成像域中心點坐標(biāo)和炮檢距；τ為成像域雙程旅行時；y和H分別為數(shù)據(jù)域中心點坐標(biāo)和炮檢距，且H∈[h-Δh,h+Δh],Δh為炮檢距間隔；w為加權(quán)函數(shù)，合適的加權(quán)函數(shù)可以加快反演的收斂速度。在各向同性介質(zhì)和直射線假設(shè)下，雙程旅行時滿足雙平方根方程

(14)

式中v為均方根速度。

正演算子式(13)的共軛轉(zhuǎn)置LT為Kirchhoff疊前時間偏移，其表達(dá)式為

(15)

1.4 速度的影響

疊前時間偏移共成像點道集屬于擴(kuò)展模型的一種，該模型除了和物理位置有關(guān)，還允許模型隨炮檢距變化。因此，不論背景速度是否正確，使用該模型均可以模擬出和觀測數(shù)據(jù)相同的反射波旅行時信息，都可以在不發(fā)生周波跳躍的前提下正確地匹配波形信息[27]，即不論速度正確與否，隨著迭代次數(shù)的增加，目標(biāo)函數(shù)式(2)中的第一項(數(shù)據(jù)擬合項)總是可以逐漸趨近于零。但是，如果速度模型不準(zhǔn)確，式(2)中的正則化算子需要適應(yīng)速度不準(zhǔn)確帶來的旅行匹配誤差。

2 數(shù)值試驗

2.1 Marmousi模型測試

圖1a為深度域Marmousi層速度模型，圖1b為Dix公式轉(zhuǎn)換后的時間域均方根速度模型。圖1c是公開的波動方程有限差分正演共炮點道集[28]經(jīng)過炮檢距分選后的(時間，中心點，炮檢距)域的疊前數(shù)據(jù)體。圖1d為圖1c數(shù)據(jù)隨機(jī)剔除50%地震道后的疊前數(shù)據(jù)體。

圖1 Marmousi模型及其正演數(shù)據(jù)體

使用均方根速度場(圖1b)和不同的疊前數(shù)據(jù)體進(jìn)行疊前時間偏移和最小二乘疊前時間偏移試驗。首先對比炮檢距為-0.325km的共炮檢距成像剖面。圖2a和圖2b為使用完整的地震數(shù)據(jù)(圖1c)分別進(jìn)行常規(guī)疊前時間偏移和最小二乘疊前時間偏移得到的共炮檢距成像剖面，可以看出：盡管Marmousi模型構(gòu)造復(fù)雜，速度橫向變化劇烈，但是時間偏移仍然可以大致準(zhǔn)確地對地下構(gòu)造形態(tài)進(jìn)行成像；經(jīng)過最小二乘反演之后的成像結(jié)果分辨率得到顯著提高。圖2c、圖2d和圖2e為使用不完整數(shù)據(jù)(圖1d)分別進(jìn)行常規(guī)疊前時間偏移、最小二乘疊前時間偏移和預(yù)條件最小二乘疊前時間偏移得到的共炮檢距成像剖面，可以看出，數(shù)據(jù)缺失給成像結(jié)果帶來許多噪聲，降低了信噪比，經(jīng)過最小二乘反演，雖然可以提高分辨率，但是并不能消除數(shù)據(jù)缺失的影響。通過加入合理的先驗信息(比如本文所采用的共成像點道集同相軸沿炮檢距平滑變化)，既可以提高成像分辨率，又可以壓制成像噪聲。

圖3a和圖3b為采用完整數(shù)據(jù)分別進(jìn)行常規(guī)疊前時間偏移和最小二乘疊前時間偏移的共成像點道集，可以看出：盡管模型復(fù)雜，速度橫向變化劇烈，但是采用Dix公式轉(zhuǎn)化而來的均方根速度場(圖1b)，共成像點道集基本可以被拉平，該信息可以作為先驗信息約束反演成像；另外，經(jīng)過最小二乘反演，成像分辨率有一定的提升，而且淺層大炮檢距的拉伸效應(yīng)也在一定程度上得到壓制。圖3c、圖3d、圖3e為采用不完整地震數(shù)據(jù)分別進(jìn)行常規(guī)疊前時間偏移、最小二乘疊前時間偏移和預(yù)條件最小二乘疊前時間偏移的共成像點道集，可以看出，和共炮檢距成像剖面類似，數(shù)據(jù)缺失在成像道集中也引入噪聲，通過最小二乘反演能夠提高分辨率，但是無法壓制這類噪聲，利用共成像點道集同相軸沿炮檢距平滑變化這個先驗信息，再進(jìn)行預(yù)條件最小二乘疊前時間偏移，不僅能夠提高分辨率，而且能夠有效壓制數(shù)據(jù)缺失帶來的噪聲。

(a)完整數(shù)據(jù)常規(guī)疊前時間偏移；

(b)完整數(shù)據(jù)最小二乘疊前時間偏移；

(c)不完整數(shù)據(jù)常規(guī)疊前時間偏移；

(d)不完整數(shù)據(jù)最小二乘疊前時間偏移；

(e)不完整數(shù)據(jù)預(yù)條件最小二乘疊前時間偏移

(a)完整數(shù)據(jù)常規(guī)疊前時間偏移；

(b)完整數(shù)據(jù)最小二乘疊前時間偏移；

(c)不完整數(shù)據(jù)常規(guī)疊前時間偏移；

(d)不完整數(shù)據(jù)最小二乘疊前時間偏移；

(e)不完整數(shù)據(jù)預(yù)條件最小二乘疊前時間偏移

對上述成像結(jié)果應(yīng)用疊前時間反偏移即可以進(jìn)行數(shù)據(jù)規(guī)則化。圖4為5.300km位置處的CMP道集規(guī)則化結(jié)果，圖5為規(guī)則化結(jié)果與原始道集誤差。

圖4 Marmousi模型不同方法數(shù)據(jù)規(guī)則化的5.300km處CMP道集對比

圖5 Marmousi模型不同方法數(shù)據(jù)規(guī)則化的5.300km處CMP道集誤差對比

由圖4和圖5可以看出，直接采用常規(guī)偏移方法也可以對數(shù)據(jù)進(jìn)行插值，但是，由于沒有考慮Hessian矩陣的影響，插值結(jié)果的振幅信息不準(zhǔn)確，需要考慮振幅匹配問題，而這個問題也是一個欠定問題，存在多解性。通過最小二乘反演可以解決振幅匹配問題，但是在數(shù)據(jù)缺失嚴(yán)重的區(qū)域，沒有先驗信息做約束，容易將成像結(jié)果中的噪聲反投影到數(shù)據(jù)空間，使插值結(jié)果信噪比較低。預(yù)條件最小二乘偏移規(guī)則化可以同時解決振幅不匹配和信噪比的問題，使得插值結(jié)果振幅屬性更加可靠、信噪比更高。

需要注意的是，對于本文方法，如果均方根速度比較準(zhǔn)確，可以采用沿著炮檢距方向的平滑算子作為預(yù)條件算子，但是如果均方根速度不準(zhǔn)確，則需要采用更加復(fù)雜的預(yù)條件算子，比如雙曲Radon變換、非穩(wěn)態(tài)平面波構(gòu)建濾波器[8]等。盡管如此，通過偏移算子的作用，成像空間依然會比數(shù)據(jù)空間更簡單，因此有助于簡化預(yù)條件算子的設(shè)計，這也是模型驅(qū)動類數(shù)據(jù)規(guī)則化方法的意義所在。

2.2 實際資料處理

為了進(jìn)一步驗證本文方法的有效性，對海上實際不規(guī)則地震數(shù)據(jù)進(jìn)行規(guī)則化處理。圖6a為海上實際觀測的疊前數(shù)據(jù)體，受拖纜漂移等因素的影響，觀測數(shù)據(jù)中存在明顯的地震道缺失；圖6b為疊加速度分析得到的速度場，本文使用該速度場進(jìn)行時間域成像，并將其作為模型先驗信息進(jìn)行數(shù)據(jù)規(guī)則化處理。

使用圖6b所示的速度場，采用常規(guī)疊前時間偏移、最小二乘疊前時間偏移、預(yù)條件最小二乘疊前時間偏移進(jìn)行成像。圖7是炮檢距為-0.425km的三種方法共炮檢距成像剖面，可以看出，受數(shù)據(jù)不規(guī)則的影響，常規(guī)偏移結(jié)果中存在明顯的成像噪聲。經(jīng)過最小二乘反演，成像分辨率得到顯著提高，但是信噪比依然較低。通過在共成像點道集引入平滑約束，進(jìn)行預(yù)條件最小二乘反演之后，不僅可以提高成像分辨率，而且可以有效壓制成像噪聲。

圖8為常規(guī)疊前時間偏移、最小二乘疊前時間偏移、預(yù)條件最小二乘疊前時間偏移后得到的位于2.360km處的炮檢距域共成像點道集，可以看出，共成像點道集基本拉平，說明圖6b所示的速度場比較準(zhǔn)確。另外，也能夠明顯的看出，受數(shù)據(jù)缺失的影響，常規(guī)疊前時間偏移成像道集的信噪比較低(圖8a)。

圖6 海上實際觀測疊前數(shù)據(jù)體(a)及其疊加速度場(b)

(a)常規(guī)疊前時間偏移；

(b)最小二乘疊前時間偏移；

(c)預(yù)條件最小二乘疊前時間偏移

圖7 海上實際數(shù)據(jù)炮檢距為-0.425m的不同方法共炮檢距偏移剖面

通過最小二乘反演進(jìn)行成像，雖然可以提高分辨率，但是無法壓制數(shù)據(jù)缺失帶來的噪聲(圖8b)。進(jìn)一步在成像道集引入同相軸沿炮檢距方向平滑變化這一先驗信息，可以顯著提高道集的信噪比，同時也能夠保持反演帶來的分辨率(圖8c)。

最后，應(yīng)用疊前時間反偏移對缺失數(shù)據(jù)進(jìn)行重建。圖9a為2.360km處的CMP道集，該道集存在明顯的數(shù)據(jù)缺失。圖9b～圖9d分別為采用常規(guī)疊前時間偏移、最小二乘疊前時間偏移和預(yù)條件最小二乘疊前時間偏移的規(guī)則化結(jié)果。從結(jié)果可以看出，對常規(guī)疊前時間偏移結(jié)果進(jìn)行反偏移，雖然也可以起到數(shù)據(jù)插值的作用，但是由于沒有考慮Hessian矩陣的影響，插值出來的振幅信息不可靠。經(jīng)過最小二乘反演，雖然可以解決振幅信息不匹配的問題，但是受不規(guī)則數(shù)據(jù)引入噪聲的影響，插值結(jié)果的信噪比較低。通過預(yù)條件化技術(shù)既可以恢復(fù)振幅信息，又能夠提高規(guī)則化后數(shù)據(jù)的信噪比。

(a)常規(guī)疊前時間偏移；

(b)最小二乘疊前時間偏移；

(c)預(yù)條件最小二乘疊前時間偏移

圖9 海上實際數(shù)據(jù)不同方法數(shù)據(jù)規(guī)則化的CMP道集

3 結(jié)論與討論

本文在最小二乘反演的理論框架下，采用最小二乘疊前時間偏移進(jìn)行數(shù)據(jù)規(guī)則化，利用合理的預(yù)條件算子改善數(shù)據(jù)規(guī)則化的處理效果。從理論模型和實際資料處理結(jié)果可以得到如下幾點結(jié)論及認(rèn)識：

(1)疊前時間偏移共成像點道集是一種擴(kuò)展模型，可以緩解模型驅(qū)動類數(shù)據(jù)規(guī)則化方法對模型的過度依賴；

(2)最小二乘疊前時間偏移可以提高成像分辨率，改善振幅保真度，是本文規(guī)則化方法的中間結(jié)果，可以用作質(zhì)量監(jiān)控和后續(xù)的成像結(jié)果分析；

(3)成像空間比數(shù)據(jù)空間的信息更加簡單，更有利于設(shè)計預(yù)條件算子；

(4)合理的預(yù)條件算子可以改善規(guī)則化的效果。當(dāng)然，在復(fù)雜地表、復(fù)雜地質(zhì)條件下，也可以采用更加復(fù)雜的成像算子，以簡化成像空間預(yù)條件算子的設(shè)計，從而進(jìn)一步改善規(guī)則化效果，但是相應(yīng)地計算成本會增加，對速度模型的依賴性也會變強(qiáng)。