時(shí)空變自適應(yīng)反Q濾波

郭志偉，曹思遠(yuǎn)*，袁殿

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249

2 中國(guó)石油大學(xué)(北京)地球物理學(xué)院，北京 102249

3 北京億華通科技股份有限公司，北京 100084

1 引言

目前，全球油氣勘探呈現(xiàn)著明顯的全面(常規(guī)和非常規(guī))、多樣(各種類(lèi)型油氣資源)、更深(深層和超深層)的特點(diǎn)[1]。隨著地震勘探要求的日益提高，高分辨率的地震資料處理對(duì)后期的油藏描述和開(kāi)發(fā)顯得至關(guān)重要。地震波在地層介質(zhì)傳播過(guò)程中的吸收衰減是影響地震資料分辨率的主要原因，主要表現(xiàn)為振幅衰減、相位畸變、頻率降低等特征[2]。地層介質(zhì)的吸收衰減作用是造成地震波能量衰減的一個(gè)重要因素，頻率越高，衰減越快，嚴(yán)重影響了地震資料的中深層成像能力和儲(chǔ)層預(yù)測(cè)的精度[3]。反Q濾波是對(duì)吸收衰減進(jìn)行有效補(bǔ)償?shù)姆椒?，能較好地補(bǔ)償振幅和校正相位，具有較準(zhǔn)確的物理機(jī)制[4]。在反Q濾波的實(shí)現(xiàn)算法方面，不同學(xué)者進(jìn)行了大量的研究[5-14]。Robinson[15-16]提出了頻率域內(nèi)插值的相位反Q濾波及其改進(jìn)算法。Hale[17-18]根據(jù)Futterman模型提出了用泰勒級(jí)數(shù)展開(kāi)近似高頻補(bǔ)償?shù)姆碤濾波。Bickel等[19]提出了基于復(fù)函數(shù)平面波的反Q濾波算法。McCarley[20]提出了基于自回歸模型的常Q濾波方法。Varela等[21]提出了Hale反Q濾波算法的改進(jìn)方法。裴江云等[22]通過(guò)麥克勞林級(jí)數(shù)展式推導(dǎo)了基于Kjartansson模型的高階近似反Q濾波公式，其一階近似式與Hale的反Q濾波式一致。Bano[23]將常Q模型相位反Q濾波拓展到層Q值模型。在補(bǔ)償過(guò)程中，由于低信噪比和機(jī)器誤差等因素，高頻成分出現(xiàn)補(bǔ)償溢出的現(xiàn)象。針對(duì)該缺點(diǎn)，人們提出了增益補(bǔ)償函數(shù)，其作用是實(shí)現(xiàn)各頻率振幅合理有效的補(bǔ)償，通常在某個(gè)頻率點(diǎn)之后逐漸減小補(bǔ)償系數(shù)(即遵循低頻充分補(bǔ)償、中頻欠補(bǔ)償、高頻不補(bǔ)償?shù)脑瓌t)。Wang[24-28]基于波場(chǎng)向下延拓提出了穩(wěn)定高效的Gabor域反Q濾波，采用的增益函數(shù)只需一個(gè)控制參數(shù)(穩(wěn)定因子法)，即可實(shí)現(xiàn)自適應(yīng)的穩(wěn)定補(bǔ)償處理。李合群等[29]在對(duì)比不同增益控制方法的基礎(chǔ)上，提出了選擇性更靈活的補(bǔ)償算子，取得了較好的補(bǔ)償效果。張固瀾等[30]針對(duì)固定增益函數(shù)容易造成深層資料高頻補(bǔ)償不足的缺點(diǎn)，提出了時(shí)變?cè)鲆婧瘮?shù)，較好地恢復(fù)地震有效頻帶內(nèi)的能量。吳吉忠等[31]提出了一種基于等效Q值的反Q濾波算法，進(jìn)行補(bǔ)償高頻能量耗散和相位畸變校正。SHI等[32]將增益限與穩(wěn)定因子設(shè)為隨時(shí)間和地層Q值連續(xù)變化的參數(shù)，提出變?cè)鲆嫦薜姆碤濾波方法，該方法能夠補(bǔ)償深層能量衰減，提高分辨率。趙巖[33]提出了一種變穩(wěn)定因子的反Q濾波方法，其中穩(wěn)定因子是地震波的傳播距離、地震記錄的頻率以及品質(zhì)因子Q的函數(shù)。WEI等[34]引入整形正則矩陣來(lái)實(shí)現(xiàn)振幅補(bǔ)償?shù)姆€(wěn)定計(jì)算，利用增益控制因子和噪聲壓制系數(shù)來(lái)壓制噪聲，改進(jìn)了穩(wěn)定化反Q濾波方法。XUE等[35]提出一種自適應(yīng)穩(wěn)定的反Q濾波方法，該方法利用同步壓縮變換對(duì)地震信號(hào)主頻段有效成分進(jìn)行重構(gòu)，同時(shí)壓制主頻段外環(huán)境噪聲，能有效地補(bǔ)償振幅衰減，抑制環(huán)境噪聲，提高地震資料分辨率。

由于反Q濾波處理過(guò)程中需要兼顧信噪比，一般情況下，地震資料的信噪比越高，則期望補(bǔ)償越充分。穩(wěn)定因子法中的增益函數(shù)只有一個(gè)控制參量，其值越大，對(duì)高頻的補(bǔ)償越充分?；谏鲜鎏攸c(diǎn)和處理目標(biāo)，通過(guò)建立增益參數(shù)和信噪比的映射關(guān)系，對(duì)資料不同信噪比的區(qū)域進(jìn)行不同程度的補(bǔ)償處理，即自適應(yīng)增益反Q濾波。

2 基本原理

2.1 反Q濾波

根據(jù)地層吸收衰減等價(jià)于偏移的思想，地表波場(chǎng)的向下延拓可以表示為：

其中，U(r,ω)是角頻率為ω傳播距離為r的平面波波場(chǎng)，k(ω)代表波數(shù)，i是虛數(shù)單位。地震剖面圖記錄了從震源到反射點(diǎn)然后再回到地表的反射波。在波數(shù)的定義中引入地層品質(zhì)因子Q：

v(ω)為與頻率有關(guān)的相速度，將復(fù)波數(shù)k代入式(1)得：

用旅行時(shí)增量Δτ=Δ(ωr)代替距離增量，得到大地Q濾波的表達(dá)式：

ν(ωr)參考頻率ωr對(duì)應(yīng)的相速度。根據(jù)Wang提出的改進(jìn)的Kolsky頻散模型，相速度ν(ω)近似為：

ωh為調(diào)諧頻率，一般取值為地震帶寬的最高頻率，可以在反Q濾波時(shí)提供相對(duì)準(zhǔn)確的相位校正。其中γ=(πQr)-1。注意：在頻散式(5)中只考慮正頻率，將之代入式(4)得：

式(6)為基于波場(chǎng)延拓的Q濾波的基本公式，根據(jù)該式得到反Q濾波式：

式中有兩個(gè)指數(shù)算子，前者為振幅補(bǔ)償算子，后者為相位補(bǔ)償算子。反Q濾波是在頻率域進(jìn)行的，最后對(duì)所有平面波求和得到補(bǔ)償后的時(shí)間域地震信號(hào)：

2.2 自適應(yīng)增益函數(shù)

直接對(duì)振幅進(jìn)行補(bǔ)償時(shí)會(huì)放大隨機(jī)噪聲與機(jī)器噪聲，產(chǎn)生極大的不穩(wěn)定性。為解決振幅補(bǔ)償值溢出和高頻噪音增大問(wèn)題，采用穩(wěn)定因子法構(gòu)建增益函數(shù)。

衰減函數(shù)為β(τ,ω)：

增益函數(shù)為Λ(τ,ω)：

其中，σ2為穩(wěn)定因子，Glim為增益參數(shù)(一般取10～100)，值越高，補(bǔ)償越充分。根據(jù)式(10)，當(dāng)ω非常大時(shí)，Λ(τ,ω)趨向于1，表明增益函數(shù)對(duì)該高頻段不補(bǔ)償也不壓制。

我們接收到的地震記錄一般包括有效地震反射波和各種各樣的噪音，不同時(shí)刻的信噪比是不同的，即記錄的信噪比是時(shí)間和空間的函數(shù)。以信噪比函數(shù)作為自適應(yīng)增益參數(shù)的選取準(zhǔn)則，建立Glim與信噪SNR(τ)之間的映射關(guān)系F，則有

式中，穩(wěn)定因子和增益參數(shù)與地震資料的信噪比有關(guān)，是時(shí)變參量，增益函數(shù)Λ~(τ,ω)具有自適應(yīng)特征。當(dāng)資料信噪比在時(shí)空域變化較大時(shí)，有必要采用式(11)來(lái)確定自適應(yīng)穩(wěn)定的增益參數(shù)進(jìn)行反Q濾波。

映射關(guān)系F可根據(jù)信噪比SNR(τ)的實(shí)際特征設(shè)計(jì)，本文給出一種簡(jiǎn)單的線性定義：

式中，待定參量Gmin、Gmax分別為最小和最大自適應(yīng)增益參數(shù)，SNRmin、SNRmax分別為SNR(τ)的最小和最大值。

3 模型測(cè)試

3.1 無(wú)噪模型

設(shè)計(jì)5個(gè)反射系數(shù)界面的模型，反射系數(shù)值為1，分別位于 200 ms、350 ms、500 ms、650 ms和 800 ms，地面初始地震子波為50 Hz主頻的Ricker子波，地下介質(zhì)Q值為100。圖1a是基于波場(chǎng)延拓正演的無(wú)噪衰減記錄(不含幾何擴(kuò)散、反射/透射損失等固有衰減)，隨著傳播距離的增加，地震子波的延續(xù)長(zhǎng)度逐漸增加、振幅衰減、相位改變。圖1b是利用式(10)反Q濾波的結(jié)果，增益參數(shù)Glim取50?？梢钥吹?，在無(wú)噪情況下，中深層地震子波的波長(zhǎng)、振幅和相位都得到較好的恢復(fù)和校正。圖2是反Q濾波前(2a)、后(2b)時(shí)頻譜的對(duì)比，吸收衰減降低了中深層記錄的主頻，反Q濾波較好地恢復(fù)了中深層記錄的中高頻成分。

3.2 含噪模型

對(duì)無(wú)噪模型添加不同程度的隨機(jī)噪音，考慮在不同信噪比條件下，對(duì)比固定增益參數(shù)和自適應(yīng)增益參數(shù)的反Q濾波效果。

圖3a是含噪記錄一，在無(wú)噪模型的基礎(chǔ)上添加了0.75%的隨機(jī)噪音，圖3b是固定增益參數(shù)反Q濾波結(jié)果，圖3c是自適應(yīng)增益參數(shù)反Q濾波結(jié)果，其中，固定增益參數(shù)(見(jiàn)圖4a)取30，自適應(yīng)增益參數(shù)(見(jiàn)圖4b)最大值取30，最小值取5。對(duì)比圖3a和圖3b，采用固定增益參數(shù)的反Q濾波較好地恢復(fù)了深層的弱反射信息(如800 ms處的同相軸)，同時(shí)擴(kuò)大了隨機(jī)噪音，大幅降低了信噪比水平；對(duì)比圖3b和圖3c，采用自適應(yīng)增益參數(shù)的反Q濾波在恢復(fù)弱反射信息的同時(shí)，有效地控制了噪音的放大，深層的信噪比得到較好的改善。圖5是反Q濾波前后的時(shí)頻譜對(duì)比，圖5a是原始含噪記錄一的時(shí)頻譜，圖5b是固定增益參數(shù)反Q濾波記錄的時(shí)頻譜，圖5c是自適應(yīng)增益參數(shù)反Q濾波記錄的時(shí)頻譜。可以看到，固定增益參數(shù)的反Q濾波對(duì)噪音中高頻成分的放大效果較為明顯，自適應(yīng)增益參數(shù)的反Q濾波則較好地控制了低信噪比區(qū)域的噪音(見(jiàn)紅色圈中部分)，提高了時(shí)頻譜的整體信噪比水平。

圖1 無(wú)噪模型反Q濾波剖面對(duì)比(a)反Q濾波前；(b)反Q濾波后Fig. 1 Comparison of inverse Q filtering pro file of the noise-free model. (a)before inverse Q filtering; (b)after inverse Q filtering

圖2 無(wú)噪模型反Q濾波時(shí)頻譜對(duì)比(a)反Q濾波前；(b)反Q濾波后Fig. 2 Spectrum comparison of the noise-free model with inverse Q filtering. (a)before inverse Q filtering; (b)after inverse Q filtering

圖6a是含噪記錄二，在無(wú)噪模型的基礎(chǔ)上添加了5%的隨機(jī)噪音，圖6b是固定增益參數(shù)反Q濾波結(jié)果，圖6c是自適應(yīng)增益參數(shù)反Q濾波結(jié)果，其中，固定增益參數(shù)(見(jiàn)圖7a)取20，自適應(yīng)增益參數(shù)(見(jiàn)圖7b)最大值取20，最小值取1。圖6a中的含噪記錄，隨著傳播距離的增加，信噪比越來(lái)越低，500 ms和650 ms處的有效反射波較弱，同相軸信息不明顯，800 ms處的反射信息淹沒(méi)在隨機(jī)噪音中，幾乎看不到同相軸。圖6b是采用固定增益參數(shù)的反Q濾波結(jié)果，與反Q濾波前剖面相比，信噪比進(jìn)一步降低，500 ms和650 ms處的反射信息淹沒(méi)在隨機(jī)噪音中，分辨能力還不如反Q濾波前的剖面。圖6c是采用自適應(yīng)增益參數(shù)的反Q濾波結(jié)果，500 ms和650 ms處的反射信息得到較好的增強(qiáng)，信噪比得到較好的控制，800 ms處的弱反射同相軸得到一定的恢復(fù)，分辨能力得到一定的改善。圖8是反Q濾波前后的時(shí)頻譜對(duì)比，圖8a是原始含噪記錄二的時(shí)頻譜，圖8b是固定增益參數(shù)反Q濾波記錄的時(shí)頻譜，圖8c是自適應(yīng)增益參數(shù)反Q濾波記錄的時(shí)頻譜。在原始含噪記錄的時(shí)頻譜中(圖8a)，650 ms以下的隨機(jī)噪音與反射波的能量在同一量級(jí)；如果采用固定增益參數(shù)進(jìn)行反Q濾波，對(duì)噪音的增強(qiáng)效果大于對(duì)有效反射波的恢復(fù)效果(圖8b中紅圈)，導(dǎo)致剖面中(圖6b)500 ms以下的同相軸淹沒(méi)在隨機(jī)噪音中，不可分辨；如果采用自適應(yīng)增益參數(shù)進(jìn)行反Q濾波，650 ms和800 ms處的有效反射振幅得到較好的恢復(fù)，有效反射的能量團(tuán)在時(shí)頻譜中仍然占主要成分(圖8c紅圈)，清晰可辨。

圖3 含噪模型一(0.75%隨機(jī)噪音)反Q濾波剖面對(duì)比(a)原始含噪記錄；(b)固定增益反Q濾波；(c)自適應(yīng)增益反Q濾波Fig. 3 Comparison of inverse Qfiltering profile in noise model 1 (0.75% random noise) (a)original noise recording; (b)fixed gain inverse Qfiltering; (c)adaptive gain inverse Qfiltering

圖4 模型一增益參數(shù)(a)固定增益參數(shù)；(b)自適應(yīng)增益參數(shù)Fig. 4 Model 1 gain parameter (a)fixed gain parameter; (b)adaptive gain parameter

圖5 含噪模型一(0.75%隨機(jī)噪音)反Q濾波時(shí)頻譜對(duì)比(a)原始含噪記錄；(b)固定增益反Q濾波；(c)自適應(yīng)增益反Q濾波Fig. 5 Spectrum comparison of inverse Qfiltering with noise model 1 (0.75% random noise) (a)original noise records;(b)fixed gain inverse Q filtering; (c)adaptive gain inverse Qfiltering

圖6 含噪模型二(5%隨機(jī)噪音)反Q濾波剖面對(duì)比(a)原始含噪記錄；(b)固定增益反Q濾波；(c)自適應(yīng)增益反Q濾波Fig. 6 Comparison of inverse Q filtering profile of noise model 2 (5% random noise) (a)original noise records; (b)fixed gain inverse Qfiltering; (c)adaptive gain inverse Qfiltering

圖7 模型二增益參數(shù)((a)固定增益參數(shù)；(b)自適應(yīng)增益參數(shù))Fig. 7 Model 2 gain parameter ((a)fixed gain parameter; (b)adaptive gain parameter)

圖8 含噪模型二(5%隨機(jī)噪音)反Q濾波時(shí)頻譜對(duì)比(a)原始含噪記錄；(b)固定增益反Q濾波；(c)自適應(yīng)增益反Q濾波Fig. 8 Spectrum comparison of inverse Qfiltering with noise model 2 (5% random noise) (a)original noise records; (b)fixed gain inverse Qfiltering; (c)adaptive gain inverse Qfiltering

4 資料處理

實(shí)際資料來(lái)自于沙漠工區(qū)，由于近地表吸收衰減嚴(yán)重，導(dǎo)致記錄的主頻較低(約為15～20 Hz)，嚴(yán)重影響了中深層的成像精度。圖9a是原始疊加剖面，目的層4 s附近，圖9b是采用固定增益參數(shù)的反Q濾波結(jié)果，圖9c是采用自適應(yīng)增益參數(shù)的反Q濾波結(jié)果，圖9d是圖9b與圖9c所示的剖面之差。圖10為固定增益參數(shù)與自適應(yīng)增益參數(shù)對(duì)比，其中，固定增益參數(shù)為40(圖10a)；自適應(yīng)增益參數(shù)(圖10b)最小值為10，最大值為60，主要分布范圍在25～45。對(duì)比圖9a、圖9b和圖9c，反Q濾波后剖面的同相軸變細(xì)，垂向分辨率得到提高(見(jiàn)黑色箭頭)；對(duì)比圖9b和圖9c，兩個(gè)剖面的分辨率相近，但圖9c的信噪比高于圖9b(如黑圈部分)，圖9d顯示了兩者之差，主要是凌亂的反射和隨機(jī)噪音，如果把這部分信息加到圖9c中，則只會(huì)降低資料的信噪比，而對(duì)分辨率的提升卻有限，這表明采用自適應(yīng)增益參數(shù)的反Q濾波能較好地控制低信噪比區(qū)域的噪音擴(kuò)大，有效改善反Q濾波剖面的信噪比。

圖11是疊加記錄反Q濾波前后的時(shí)頻譜對(duì)比，圖11a是原始疊加記錄的時(shí)頻譜，圖11b是固定增益參數(shù)反Q濾波記錄的時(shí)頻譜，圖11c是自適應(yīng)增益參數(shù)反Q濾波記錄的時(shí)頻譜。對(duì)比發(fā)現(xiàn)，反Q濾波前疊加記錄的時(shí)頻譜(圖11a)能量主要分布在10～30 Hz，經(jīng)反Q濾波處理后，時(shí)頻譜(圖11b和圖11c)能量主要集中在15～40 Hz之間，中高頻成分得到較好的恢復(fù)；圖11c的高頻成分比圖11b弱，這是由于該段記錄的信噪比相對(duì)較低，采用了較小的增益參數(shù)(小于固定增益參數(shù)40)進(jìn)行反Q濾波。值得注意的是，在圖11b中，4000 ms處的主頻(約為40 Hz)明顯高于4500 ms處的主頻(約為32 Hz)，在圖11c中，4000 ms和4500 ms處的主頻基本一致(約為30 Hz)，這說(shuō)明自適應(yīng)增益參數(shù)能較好地調(diào)節(jié)不同信噪比區(qū)域的補(bǔ)償效果。

圖9 疊后資料反Q濾波剖面對(duì)比(a)原始疊加記錄；(b)固定增益反Q濾波；(c)自適應(yīng)增益反Q濾波；(d)剖面b與剖面c之差Fig. 9 Comparison of inverse Qfiltering profile of post-stack data (a)original overlay records; (b)fixed gain inverse Qfiltering;(c)adaptive gain inverse Qfiltering; (d)the difference between section b and section c

5 結(jié)論

穩(wěn)定的反Q濾波處理為了兼顧分辨率和信噪比，需要不斷調(diào)試增益參數(shù)，將原始資料的信噪比作為先驗(yàn)信息，用于約束增益參數(shù)的選取，可實(shí)現(xiàn)地震資料自適應(yīng)的反Q濾波處理。采用自適應(yīng)增益參數(shù)的反Q濾波，根據(jù)資料不同區(qū)域的信噪比水平，能自動(dòng)調(diào)節(jié)補(bǔ)償效果，有效控制低信噪比區(qū)域的補(bǔ)償，避免噪音被過(guò)度放大。需要注意的是，為避免反Q濾波剖面的突變，自適應(yīng)增益參數(shù)應(yīng)進(jìn)行適當(dāng)?shù)亩S平滑處理。

圖10 實(shí)際資料增益參數(shù)。(a)固定增益；(b)自適應(yīng)增益Fig. 10 Actual data gain parameter. (a)fixed gain; (b)adaptive gain

圖11 疊加記錄反Q濾波時(shí)頻譜對(duì)比。(a)原始疊加記錄；(b)固定增益反Q濾波；(c)自適應(yīng)增益反Q濾波Fig. 11 Spectrum comparison inverse Qfiltering is recorded in superposition records. (a)original superposition records;(b)fixed gain inverse Q filtering; (c)adaptive gain inverse Qfiltering