時變子波提取及其在地震資料智能處理中的應(yīng)用

戴永壽, 張 鵬, 萬 勇, 孫偉峰, 張彧豪, 韓浩宇, 吳莎莎

(1.中國石油大學(xué)(華東)海洋與空間信息學(xué)院,山東青島 266580; 2.中國石油大學(xué)(華東)控制科學(xué)與工程學(xué)院,山東青島 266580)

近年來,中國在塔里木盆地庫車等地區(qū)發(fā)現(xiàn)了具有相當(dāng)規(guī)模的油氣藏資源,但主體目的層位于深層或超深層,且地表地形起伏大、地下構(gòu)造異常復(fù)雜,導(dǎo)致反射信號能量弱、地震資料品質(zhì)差、成像精度低[1]。常規(guī)地震數(shù)據(jù)處理方法已無法滿足新形勢下高分辨率地震勘探的需求[2-3],亟需融合智能信息處理新技術(shù)為深層油氣藏的勘探與開發(fā)提供新思路。地震子波是影響地震資料分辨率的關(guān)鍵因素,也是地震正反演、速度建模和偏移成像的重要基礎(chǔ)。早期研究通常假設(shè)地震子波時不變,但該假設(shè)僅適用于中淺層及簡單構(gòu)造資料的處理;而對于深層復(fù)雜構(gòu)造,地層的吸收與頻散作用導(dǎo)致子波在傳播過程中出現(xiàn)高頻衰減和相位畸變,因此具有時變性[4]。時變地震子波的準(zhǔn)確提取對于提高地震資料處理和解釋精度具有重要的意義。筆者將采用智能信息處理技術(shù)通過分別估計振幅譜與相位譜提取時變地震子波,并應(yīng)用于反褶積、偏移等地震數(shù)據(jù)處理主要環(huán)節(jié),提出基于稀疏表示的疊后非平穩(wěn)反褶積方法以及利用時變子波提高疊前偏移成像精度的方法。

1 時變子波提取方法

非平穩(wěn)褶積模型[5]可表示為

x(t)=w(t,τ)*r(t)+n(t).

(1)

式中,x(t)為非平穩(wěn)地震記錄;w(t,τ)為時變子波;r(t)為反射系數(shù)序列;n(t)為噪聲干擾,實際資料中既包含高斯噪聲,也包含非高斯噪聲。采用改進(jìn)的廣義S變換[6]對式(1)進(jìn)行時頻變換和濾波處理,則x(t)、w(t,τ)和r(t)在時頻域的關(guān)系可以表示為

X(t,f)=W(t,f)·R(t,f).

(2)

式中，X(t,f) 、W(t,f)和R(t,f)分別為地震記錄、時變子波和反射系數(shù)序列的時頻譜。

再對等式兩邊進(jìn)行對數(shù)變換,則地震記錄、子波和反射系數(shù)的對數(shù)譜之間的關(guān)系如下:

ln|X(t,f)|=ln|W(t,f)|+ln|R(t,f)|.

(3)

式中,ln|R(t,f)|為高頻信號,具有振蕩劇烈的特點;ln|W(t,f)|為低頻信號,具有光滑連續(xù)的特點。因此時變子波振幅譜的提取問題即可轉(zhuǎn)化為信號去噪問題?；パa集合經(jīng)驗?zāi)B(tài)分解(complementary ensemble empirical mode decomposition,CEEMD)是目前應(yīng)用效果最好的智能去噪方法[7],本文中利用CEEMD將不同時刻地震記錄的對數(shù)振幅譜分解為一組具有不同振蕩尺度的分量,濾除振蕩劇烈分量,重構(gòu)光滑連續(xù)分量,即實現(xiàn)了時變子波振幅譜的準(zhǔn)確提取。

在子波振幅譜提取的基礎(chǔ)上,研究采用智能尋優(yōu)方法提取時變子波相位譜。對子波振幅譜進(jìn)行希爾伯特變換即可估計時變最小相位子波[8]。應(yīng)用極零點模型描述最小相位子波Z域的系統(tǒng)函數(shù):

(4)

式中,a、b為待估計的參數(shù);m、n為模型的階次。針對式(4)所示的線性差分系統(tǒng),采用遞推最小二乘法辨識模型中的參數(shù),然后令分子和分母分別等于零,能夠估計得到最小相位子波Z域的極零點;再根據(jù)不同相位特征的子波極零點分布的規(guī)律及差異,對最小相位子波的極零點關(guān)于單位圓進(jìn)行對稱變換,并在局部相似度準(zhǔn)則的約束下確定最優(yōu)的組合,實現(xiàn)時變混合相位子波的準(zhǔn)確提取[9]。具體實現(xiàn)流程如下:

(1) 計算地震記錄的包絡(luò)。

(2) 利用不同極零點組合所對應(yīng)的子波相位譜對地震記錄進(jìn)行相位校正。

(3) 分別計算相位校正后地震記錄的局部包絡(luò)h(t)與原地震記錄的局部包絡(luò)s(t)的相似系數(shù)。計算公式為

(5)

式中，c為相似系數(shù)。

(4) 當(dāng)相似系數(shù)的倒數(shù)取最大時,對應(yīng)的極零點分布為最優(yōu)組合,對應(yīng)的子波即為最優(yōu)的時變子波。

2 疊后非平穩(wěn)反褶積

反褶積是壓縮地震子波、提高地震資料分辨率最為直接有效的處理技術(shù)。疊后水平層狀介質(zhì)中的反射系數(shù)序列具有稀疏性,而人工智能領(lǐng)域的稀疏表示理論與方法正是解決稀疏反演問題的有效手段。

式(1)的矩陣形式可以表示為

x=Wr+n.

(6)

式中，x為地震記錄向量;W為時變子波矩陣;r為反射系數(shù)序列向量;n為噪聲向量。

針對現(xiàn)有方法難以在高斯與非高斯噪聲共同干擾下準(zhǔn)確反演反射系數(shù)序列的問題,利用L2范數(shù)能夠壓制高斯噪聲、L1范數(shù)對非高斯噪聲具有較好抑制能力的特性,建立混合范數(shù)約束的反演目標(biāo)函數(shù)[10]:

(7)

式中,μ為正則化參數(shù);λ1和λ2為權(quán)重系數(shù),用來調(diào)節(jié)L1和L2范數(shù)約束的擬合誤差項的權(quán)值,取值范圍為[0, 1],且λ1+λ2=1。分裂Bregman迭代是對包含L1范數(shù)約束項的目標(biāo)函數(shù)進(jìn)行最小化求解最為快速有效的算法,其主要思想是利用中間變量將L1范數(shù)約束中的表達(dá)式分離出來,從而將單變量單目標(biāo)函數(shù)分裂為多變量多目標(biāo)函數(shù)的交替迭代,使其更容易求解。令d1=Wr-x,d2=r,反演目標(biāo)函數(shù)[11]可以分裂為

(8)

(9)

(10)

(11)

(12)

對式(8)直接求導(dǎo),再令導(dǎo)數(shù)等于0,即可估計每次迭代中的反射系數(shù)序列

(13)

式(9)和(10)均滿足鄰近算子的定義式[11-12],可分別改寫為

(14)

(15)

式中,sign(y)為符號函數(shù),y>0時取值為1,反之取值為-1。綜上所述,利用式(13)～(15)、(11)和(12)的交替迭代即可在高斯與非高斯噪聲共同干擾下實現(xiàn)反射系數(shù)序列的準(zhǔn)確反演。

3 疊前偏移成像方法

疊前資料的偏移成像方法由于包含豐富的波形信息近年來得到廣泛關(guān)注,而震源子波的估計精度不佳以及地震記錄頻帶較窄是影響偏移成像效果的主要因素[13]。理想的震源子波為尖脈沖,但實際上震源子波為能量集中在前端、具有一定延續(xù)長度的波形,符合最小相位特征[11]。對于疊前炮數(shù)據(jù),零偏移距地震道是最先接收到地震波的地震記錄,則該道地震記錄中初始時刻的子波最接近震源子波。因此本文中首先從零偏移距地震記錄中提取時變子波,然后對初始時刻的子波振幅譜進(jìn)行希爾伯特變換即可實現(xiàn)具有最小相位特征的震源子波估計。

由于疊前地震資料具有低信噪比、強衰減的特點,地震資料高低頻缺失嚴(yán)重、頻帶較窄,直接對疊前資料進(jìn)行偏移會導(dǎo)致虛假成像和模糊成像,不能準(zhǔn)確反映地下構(gòu)造真實的位置和形態(tài),進(jìn)而影響后續(xù)處理和解釋精度。本文中采用反Q濾波結(jié)合維納反褶積的方法對疊前地震資料進(jìn)行拓頻處理。地層品質(zhì)因子Q是描述介質(zhì)的非完全彈性特征,表征地層吸收衰減的重要參數(shù)。研究者們先后提出了多種Q值估計方法,其中頻率域的譜比法是最常用且效果最好的方法[14]:

(16)

式中，A(t0,f)和A(t1,f)為不同時刻提取的子波振幅譜。式(16)是關(guān)于f的線性函數(shù),其斜率為-2π(t1-t0)/Q。對不同時刻提取的子波振幅譜A(t0,f)和A(t1,f)作差并在有效頻帶內(nèi)進(jìn)行線性擬合求得斜率后,即可估計得到品質(zhì)因子Q[15]。在Q值估計的基礎(chǔ)上,利用下列公式可實現(xiàn)對地震資料的反Q濾波處理[16]:

XQ(t,f)=X(t,f)Λ(t,f),

(17)

Λ(t,f)=[β(t,f)+σ2]/[β2(t,f)+σ2],

(18)

(19)

式中,X(t,f)和XQ(t,f)分別為反Q濾波處理前后地震記錄的時頻譜;β(t,f)為濾波因子;σ2為穩(wěn)定因子;fr為參考頻率;Λ(t,f)為振幅補償因子。反Q濾波能夠?qū)ΟB前資料中的高頻衰減進(jìn)行有效補償,使得深層部分的能量得到提升。完成反Q濾波處理后,繼續(xù)對地震資料進(jìn)行維納反褶積[17]:

(20)

4 數(shù)值仿真驗證

為了驗證時變子波提取方法的有效性,本文中采用ARMA模型描述原始子波,并利用下式構(gòu)造時變子波[21]:

w(t,τ)=

(21)

圖1 合成地震記錄Fig.1 Synthetic seismogram

采用振幅譜與相位譜分別估計的時變子波提取方法從圖1(c)中估計不同時刻的子波,結(jié)果如圖2所示,圖2(a)～(d)分別為111、381、619及892 ms的提取結(jié)果。

圖2 時變地震子波提取結(jié)果Fig.2 Time-varying wavelet extraction results

從圖2中可以看出本文方法能夠準(zhǔn)確提取不同時刻的子波,且能真實反映動態(tài)衰減特征。利用提取的子波對圖1(c)進(jìn)行反褶積處理,應(yīng)用基于稀疏表示的非平穩(wěn)反褶積方法反演反射系數(shù)序列,結(jié)果如圖3所示,在一定強度的高斯與非高斯噪聲共同干擾下,本文方法也能有效壓縮地震子波,較為準(zhǔn)確地估計反射系數(shù)序列。

圖3 反射系數(shù)序列反演結(jié)果Fig.3 Results of reflection coefficient inversion

震源子波估計不準(zhǔn)及地震資料頻帶較窄是造成疊前偏移成像效果不佳的主要原因。本文中通過正演模擬合成40炮記錄,每炮記錄包括901道,采樣點3 750個,以第20炮為例(圖5(a))。首先采用本文方法從零偏移距地震道(圖5第401道)中估計不同時刻的子波,然后對初始時刻子波的振幅譜進(jìn)行希爾伯特變換,實現(xiàn)震源子波估計?；谑?16)進(jìn)行品質(zhì)因子Q值估計,并利用式(17)～(19)對地震資料進(jìn)行反Q濾波處理,以補償?shù)貙铀p。再采用式(20)對地震資料進(jìn)行維納反褶積,以壓縮子波、拓寬地震記錄頻帶。從圖4可以看出利用本文方法能夠準(zhǔn)確提取震源子波,有效解決震源子波估計的難題。

圖4 震源子波估計結(jié)果Fig.4 Estimation of source wavelet

圖5 模擬單炮記錄Fig.5 Simulated single shot record

圖5表明利用提取的子波進(jìn)行拓頻處理后地震資料的深層能量得到了提升,邊界更加清晰。最后應(yīng)用最小二乘逆時偏移方法對地震資料進(jìn)行成像,結(jié)果如圖6所示。對比圖6(a)和(b)可以看出,在不經(jīng)過處理的情況下直接進(jìn)行偏移成像,地下構(gòu)造出現(xiàn)了虛假的雙重邊界,而利用本文方法估計震源子波并實現(xiàn)高頻延拓,再進(jìn)行偏移處理能夠有效提高成像精度。

圖6 偏移成像結(jié)果Fig.6 Migration imaging results

5 實際資料處理

為了驗證基于時變子波的地震資料智能處理技術(shù)在實際地震數(shù)據(jù)處理方面的有效性,首先選取勝利油田某區(qū)塊疊后地震剖面(圖7(a))并提取時變地震子波,然后利用分裂Bregman迭代法對地震剖面進(jìn)行反褶積處理,結(jié)果如圖7(b)所示。對比圖7(a)和(b)可以得出,經(jīng)本文方法處理后疊后地震剖面的分辨率更高,同相軸更細(xì),能夠為地震資料的解釋提供重要的參考。圖8為圖7反褶積結(jié)果的頻譜,從中可以看出經(jīng)處理后地震資料的頻帶寬度得到了明顯的提升。

圖7 實際疊后地震資料反褶積結(jié)果Fig.7 Deconvolution result of actual poststack seismic data

圖8 反褶積前后頻譜的對比Fig.8 Comparison of average normalized amplitude spectra of images shown

為了進(jìn)一步驗證本文方法在疊前偏移成像中的重要作用,從圖9(a)所示疊前炮數(shù)據(jù)中估計震源子波,再利用反Q濾波與維納反褶積相結(jié)合的方法對地震資料進(jìn)行拓頻處理,結(jié)果如圖9(b)所示。可以看出，深層能量得到了有效補償,同相軸更細(xì)。在此基礎(chǔ)上應(yīng)用最小二乘逆時偏移方法進(jìn)行成像,結(jié)果如圖10所示。

圖9 疊前炮數(shù)據(jù)的拓頻處理Fig.9 Frequency expanding of prestack seismic data

圖10 疊前地震資料的偏移成像結(jié)果Fig.10 Migration imaging results of prestack seismic data

從橢圓標(biāo)記處可以看出經(jīng)處理后的成像結(jié)果所受的干擾更少,對層位的反映更清晰。將虛線框標(biāo)記的區(qū)域進(jìn)行放大,表明本文方法能夠使同相軸連續(xù)性更好,振幅能量更均衡。

6 結(jié)束語

初步探究了時變子波提取及其在地震資料智能處理中的應(yīng)用。利用智能信息處理技術(shù)通過分別估計不同時刻子波的振幅譜和相位譜實現(xiàn)了時變地震子波提取;提出了基于稀疏表示的疊后非平穩(wěn)反褶積方法,在復(fù)雜噪聲干擾下反演反射系數(shù)序列,提高地震資料的分辨率;研究了時變子波在疊前偏移中的應(yīng)用方法,通過震源子波估計與地震記錄的拓頻處理,進(jìn)一步提高了成像精度。今后將繼續(xù)研究利用隨機共振方法的低信噪比地震資料無損降噪技術(shù)、基于卷積神經(jīng)網(wǎng)絡(luò)的疊前子波提取方法以及基于壓縮感知的地震數(shù)據(jù)低頻補償方法,最終有效應(yīng)用于塔里木深層地震資料的處理,實現(xiàn)高精度全波形反演與偏移成像。