基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波方法

趙 巖，毛寧波，陳 旭

（1.油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室（長江大學(xué)），武漢 430100；2.長江大學(xué)地球物理與石油資源學(xué)院，武漢 430100；3.長江大學(xué)錄井技術(shù)與工程研究院，湖北荊州 434023）

0 引言

因地層吸收作用，地震波在地下介質(zhì)傳播中會發(fā)生振幅衰減和相位畸變，從而主頻降低，頻帶變窄，地震資料的分辨率降低［1-4］。反Q濾波［5-7］、Gabor反褶積［8-9］和Q偏移［10-11］等方法通常被用來進(jìn)行地震波的衰減補(bǔ)償，以提高地震資料的分辨率。就反Q濾波而言，由于硬件存儲精度的限制，為解決反Q濾波振幅補(bǔ)償會出現(xiàn)的數(shù)值穩(wěn)定性問題［5］，學(xué)者們提出了一些穩(wěn)定的反Q濾波方法，如穩(wěn)定因子法［5，12-13］、迭代法［14］、正則化法［15-16］、反演方法［3］、頻域方法［17］等。這些方法的振幅補(bǔ)償函數(shù)是固定的，對于任何地震記錄，其振幅補(bǔ)償強(qiáng)度和頻帶補(bǔ)償寬度均相同，但地震記錄是非平穩(wěn)信號，其振幅能量分布和有效頻帶范圍都是隨時(shí)間或深度變化而變化的。當(dāng)?shù)卣鹩涗浐性肼晻r(shí)，常規(guī)穩(wěn)定的反Q濾波方法無法結(jié)合地震記錄有效頻帶范圍對其進(jìn)行差別補(bǔ)償，導(dǎo)致其在某些時(shí)刻會同時(shí)增強(qiáng)有效振幅和噪聲能量，影響反Q濾波后地震記錄的信噪比。Zhao等［18］結(jié)合地震記錄的非平穩(wěn)性，提出一種時(shí)頻域信噪比的估算方法，并將其用于穩(wěn)定因子法反Q濾波，該方法基于時(shí)頻域信噪比估算的有效頻帶范圍，所以只在有效頻帶范圍內(nèi)進(jìn)行穩(wěn)定的反Q濾波補(bǔ)償，在一定程度上壓制了噪聲能量，但該方法振幅補(bǔ)償函數(shù)的增益限未與有效頻帶范圍相結(jié)合，對某些時(shí)刻的振幅能量也存在一定的壓制。張固瀾等［19］提出一種穩(wěn)定的自適應(yīng)增益限反Q濾波方法，其考慮了地震記錄的有效頻帶寬度，并使增益限和穩(wěn)定因子自適應(yīng)于有效頻帶的截止頻率，同時(shí)指出可通過分析地震資料的頻譜來獲取有效頻帶的截止頻率。

在自適應(yīng)增益限反Q濾波的基礎(chǔ)上，結(jié)合時(shí)頻域信噪比的概念，提出一種基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波方法。首先，從地震記錄相鄰地震道中提取信號和噪聲統(tǒng)計(jì)量，以多道平均的方式估算地震記錄的時(shí)頻域信噪比；然后，根據(jù)信噪比閾值確定有效頻帶截止頻率；最后，結(jié)合自適應(yīng)于截止頻率的增益限反Q濾波方法對地震記錄進(jìn)行時(shí)變增益振幅補(bǔ)償。

1 基本原理

1.1 穩(wěn)定的反Q 濾波方法

反Q濾波以一維波動方程［5］為基礎(chǔ)

式（1）中：U(z,ω) 為傳播距離為z時(shí)角頻率為ω的平面波；k為波數(shù)，?？紤]深度增量與走時(shí)增量的關(guān)系，則有

式（2）中：τ為傳播時(shí)間；Δτ為時(shí)間間隔；U(τ+Δτ,ω)和U(τ,ω)分別是τ+Δτ和τ時(shí)刻的波場為虛數(shù)單位；ωh為地震頻帶內(nèi)與最高頻率有關(guān)的調(diào)諧頻率；Q為品質(zhì)因子；ω為角頻率；γ=1/(πQ)πQ。因此，反Q濾波的延拓方程［15］可寫為

由于地震波衰減到一定程度時(shí)，環(huán)境噪聲的能量會高于信號能量，因此，反Q濾波會對環(huán)境噪聲進(jìn)行放大，造成數(shù)值不穩(wěn)定。這種不穩(wěn)定主要是由式（3）中的振幅校正算子所導(dǎo)致，而相位校正算子則是穩(wěn)定的。針對反Q濾波振幅校正存在的不穩(wěn)定性問題，通過在分子和分母中引入穩(wěn)定因子，可實(shí)現(xiàn)穩(wěn)定的反Q濾波算法［5］

1.2 自適應(yīng)增益限反Q 濾波方法

自適應(yīng)增益限反Q濾波方法［19］由以下公式實(shí)現(xiàn)

式（5）中：t為旅行時(shí)；ωd(t)為有效頻范圍的截止角頻率；A(t,ω)是振幅補(bǔ)償函數(shù)；c(t)是振幅補(bǔ)償函數(shù)U(t,ω)的時(shí)變增益限，dB；其與地震波有效頻范圍的截止角頻率相對應(yīng)；[c2(t)-2c(t)]-1是穩(wěn)定因子，主要用于壓制高頻噪聲和吉布斯效應(yīng)，以解決反Q濾波的非穩(wěn)定性問題；～ (t,ω) 和(t,ω) 分別為反Q濾波前后t時(shí)刻的地震記錄波場。

為了驗(yàn)證自適應(yīng)增益限反Q濾波的有效性，分別利用常規(guī)穩(wěn)定的反Q濾波方法和自適應(yīng)增益限反Q濾波方法對含噪聲合成地震記錄進(jìn)行振幅補(bǔ)償。圖1（a）為不含噪聲的合成地震記錄，子波主頻為35 Hz，采樣間隔為2 ms。圖1（b）為對圖1（a）中的記錄進(jìn)行Q=100 的衰減模擬結(jié)果，可以看到，隨著旅行時(shí)的增加，子波振幅逐漸衰減，相位畸變，波形產(chǎn)生變化。圖1（c）為均值為0 的隨機(jī)噪聲，其標(biāo)準(zhǔn)差為圖1（a）中記錄的最大幅度的0.1%。將圖1（c）中的隨機(jī)噪聲添加到圖1（b）中的衰減地震記錄中，得到含有噪聲的衰減合成地震記錄，分別在0.3 s，0.6 s，0.9 s，1.2 s，1.5 s 和1.8 s 處提取子波，計(jì)算其歸一化對數(shù)振幅譜［圖1（d）］，可以看到，隨著旅行時(shí)的增加，子波振幅譜的高頻成分逐漸減少，有效振幅頻帶逐漸“變窄”。

圖1 合成地震記錄Fig.1 Synthetic seismic records

圖2（a）為t=0.3 s，0.6 s，0.9 s，1.2 s，1.5 s，1.8 s時(shí)穩(wěn)定的反Q濾波振幅補(bǔ)償曲線，增益限為20 dB，可以看到，雖然不同時(shí)刻的振幅補(bǔ)償頻帶范圍不同，但是其振幅補(bǔ)償曲線的最大值是相同的。圖2（b）為t=0.3 s，0.6 s，0.9 s，1.2 s，1.5 s，1.8 s 時(shí)自適應(yīng)增益限的反Q濾波振幅補(bǔ)償曲線，其增益限隨著有效頻帶截止頻率的變化而變化。從圖1（d）中可以讀出上述6 個(gè)時(shí)刻的截止角頻率ωd分別為200π，180 π，164 π，156π，146 π 和134 π，由圖2（b）可看到，不同時(shí)刻的振幅補(bǔ)償頻帶范圍不同，而且其振幅補(bǔ)償曲線的最大值也是不同的，這體現(xiàn)了自適應(yīng)增益限反Q濾波的優(yōu)勢所在。分別用穩(wěn)定的反Q濾波和自適應(yīng)增益限反Q濾波方法對圖1 中的含噪聲衰減記錄進(jìn)行反Q濾波，其結(jié)果分別如圖2（c）和圖2（d）所示?？梢钥吹?，穩(wěn)定的反Q濾波方法可補(bǔ)償子波的振幅能量，但也對噪聲能量進(jìn)行了放大，而自適應(yīng)增益限反Q濾波方法既補(bǔ)償了子波的振幅能量，同時(shí)也壓制了大部分隨機(jī)噪聲能量。

圖2 圖1 中衰減合成地震記錄的反Q 濾波結(jié)果Fig.2 Inverse Q filtering results of the attenuated synthetic seismic records in Fig.1

1.3 基于時(shí)頻域信噪比的自適應(yīng)增益限反Q 濾波方法

由式（5）可知，自適應(yīng)增益限反Q濾波的振幅補(bǔ)償結(jié)果依賴于地震記錄的有效頻帶范圍，即有效頻帶的截止角頻率ωd(t) 。因此，本次研究將時(shí)頻域信噪比的概念和自適應(yīng)增益限反Q濾波方法結(jié)合起來，提出一種基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波方法。該方法首先在相關(guān)理論的基礎(chǔ)上從地震記錄相鄰地震道中提取信號和噪聲的平均瞬時(shí)功率譜，并采用多道平均的方式估算地震記錄的時(shí)頻域信噪比［18］

式（6）中：Px(t,ω)，Ps(t,ω)和Pn(t,ω)分別是地震道、信號和噪聲的平均瞬時(shí)功率譜。

然后，根據(jù)地震記錄的信噪比特征及地震數(shù)據(jù)處理的實(shí)際需求，選擇閾值C0，并利用SNR(t,ω)≥C0計(jì)算地震記錄隨旅行時(shí)t變化的有效頻帶范圍

從式（7）可看出，該頻帶范圍會隨閾值C0的變化而變化，當(dāng)C0變大時(shí)，有效頻帶“變窄”，反之，有效頻帶“變寬”。最后，確定反Q濾波中有效頻帶的截止角頻率ωd(t)

并利用式（5）進(jìn)行自適應(yīng)增益限反Q濾波。

圖3 中A 為理論合成的不含衰減地震記錄，共512 個(gè)采樣點(diǎn)，采樣間隔為2 ms，在此作為參考道。B 為對A中的地震記錄進(jìn)行Q=50 的衰減模擬結(jié)果，同時(shí)添加隨機(jī)噪聲。C—E 分別為閾值C0=0.005，0.040 和0.150 情況下的基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波結(jié)果。由圖3 可以看到，不同的閾值C0，其反Q濾波結(jié)果具有明顯差異。結(jié)合圖3和圖4 可以看到，不同的閾值C0可得到不同的截止角頻率ωd(t)，進(jìn)而得到不同的反Q濾波結(jié)果。以參考道為標(biāo)準(zhǔn)，當(dāng)C0較小時(shí)，有效頻帶“過寬”，會出現(xiàn)過補(bǔ)償，放大了噪聲能量，如圖3 中C 所示；當(dāng)C0較大時(shí)，有效頻帶“過窄”，深層地震記錄振幅出現(xiàn)欠補(bǔ)償，如圖3 中E 所示；當(dāng)C0=0.04 時(shí)，其反Q濾波結(jié)果最接近參考道，如圖3 中D 所示，此時(shí)閾值C0得到的有效頻帶截止頻率較為合適。

圖3 含噪聲合成地震記錄及不同閾值情況下的基于時(shí)頻域信噪比的自適應(yīng)增益限反Q 濾波結(jié)果Fig.3 Noise-added synthetic seismic records and the results of self-adaptive gain-limit inverse Q filtering based on SNR in time-frequency domain under different thresholds

圖4 圖3 中衰減地震記錄B 的時(shí)頻域信噪比（a）及不同閾值情況下的有效頻帶截止頻率曲線（b）Fig.4 SNR in time-frequency domain of the attenuated seismic record B in Fig.3（a）and the cut-off frequency curves of effective band under different thresholds（b）

對于如何選取閾值C0，可以根據(jù)地震數(shù)據(jù)的時(shí)頻域特征采用計(jì)算時(shí)頻譜的方法定性分析，再結(jié)合圖3中的理論合成記錄舉例說明。計(jì)算圖3 中B—E 地震記錄的時(shí)頻譜（圖5）。由圖5（a）可看出，其淺層振幅能量較強(qiáng)，深層振幅能量較弱，噪聲能量分布在整個(gè)頻帶范圍內(nèi)；由圖5（b）可看出，其深層出現(xiàn)異常強(qiáng)振幅能量，表現(xiàn)為高頻噪聲能量的增強(qiáng)；由圖5（c）可看出，補(bǔ)償后的深層振幅能量與淺層基本一致，時(shí)頻譜振幅能量比較均衡；把圖5（d）與圖5（a）相比，其深層振幅能量有所增強(qiáng)，但與圖5（c）相比，其振幅能量只得到部分補(bǔ)償，表現(xiàn)為有效振幅能量相對較弱。由此可見，確定閾值C0=0.04比較合適，對于實(shí)際地震數(shù)據(jù)，亦可采用上述方法優(yōu)選合適的閾值。

圖5 圖3 中衰減地震記錄在不同閾值情況下反Q 濾波結(jié)果的時(shí)頻譜Fig.5 Time-frequency spectra of the inverse Q filtering results of the attenuated seismic records in Fig.3 under different thresholds

2 理論數(shù)據(jù)測試

為進(jìn)一步測試地震數(shù)據(jù)信噪比對研究方法的影響，采用3 個(gè)不同信噪比的理論合成地震記錄進(jìn)行常規(guī)穩(wěn)定的反Q濾波和基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波，并對其結(jié)果進(jìn)行對比分析（圖6）。圖6 中的A 為不含噪聲的合成地震記錄，共512 個(gè)采樣點(diǎn)，采樣間隔為2 ms，在此作為參考道。對A 中地震記錄進(jìn)行Q=100 的衰減模擬，并分別添加標(biāo)準(zhǔn)差為A 中地震記錄最大幅度1.5%，3.0% 和6.0% 的隨機(jī)噪聲，分別如圖6 中的B 所示，可以看到，隨著傳播時(shí)間的增大，地震記錄的有效振幅能量逐漸變?nèi)?，不同記錄的噪聲能量逐漸增強(qiáng)，信噪比逐漸降低。分別對3 個(gè)不同信噪比的衰減地震記錄進(jìn)行常規(guī)穩(wěn)定的反Q濾波，可以看到，隨著信噪比的降低，常規(guī)穩(wěn)定的反Q濾波對噪聲能量的放大越來越明顯，噪聲能量越來越強(qiáng)，如圖6中的C 所示。作為對比，分別對它們進(jìn)行基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波，如圖6 中的D所示。由圖6 可以看到，當(dāng)原始衰減記錄中的噪聲能量較弱時(shí)，本文方法可以有效補(bǔ)償振幅能量，并壓制噪聲能量，如圖6（a）和圖6（b）中的D 所示。在噪聲能量繼續(xù)增強(qiáng)，信噪比繼續(xù)降低的情況下，本文方法在補(bǔ)償有效振幅能量的同時(shí)，亦會部分放大噪聲能量，如圖6（c）中的D 所示，但與常規(guī)穩(wěn)定的反Q濾波結(jié)果相比，其對噪聲能量的壓制效果非常突出。

圖6 基于時(shí)頻域信噪比的自適應(yīng)增益限反Q 濾波Fig.6 Self-adaptive gain-limit inverse Q filtering based on SNR in time-frequency domain

3 實(shí)際資料應(yīng)用

將本次研究方法用于實(shí)際地震記錄，疊后地震記錄如圖7（a）所示，雖然地震記錄的整體信噪比較高，但仍有未被完全壓制的隨機(jī)干擾。圖7（b）為對圖7（a）中的地震剖面進(jìn)行穩(wěn)定的反Q濾波結(jié)果，其中Q值為150，增益限為10 dB?？梢钥吹?，深層地震記錄的振幅能量增強(qiáng)，且分辨率有所提高，但地震記錄也變得“粗糙”，這是因?yàn)榉碤濾波也增強(qiáng)了隨機(jī)噪聲的能量。由圖7（c）可以看到，其分辨率明顯得到改善，且地震剖面比圖7（b）“干凈”，這是因?yàn)楸敬窝芯糠椒▔褐屏瞬糠蛛S機(jī)噪聲的能量。對圖7 黃色虛線框內(nèi)的地震剖面進(jìn)行放大顯示，如圖8 所示。利用穩(wěn)定的反Q濾波和本次研究方法處理得到的地震記錄的分辨率都得到明顯提高，但是穩(wěn)定的反Q濾波方法處理后的剖面具有明顯的“毛刺”，隨機(jī)干擾能量較強(qiáng)；而利用本此研究方法處理得到的結(jié)果，明顯改善了地震記錄的分辨率，同時(shí)也壓制了噪聲能量。圖9 為圖8 中相應(yīng)地震記錄的歸一化對數(shù)振幅譜，其中藍(lán)線為原始地震記錄，紅線為穩(wěn)定的反Q濾波結(jié)果，黃線為本次研究方法的處理結(jié)果?？梢钥吹?，與穩(wěn)定的反Q濾波方法相比，基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波方法在補(bǔ)償高頻振幅能量的同時(shí)，壓制了部分高頻帶內(nèi)的隨機(jī)噪聲，有效提高了地震記錄的分辨率和信噪比。

圖7 實(shí)際地震記錄及不同反Q 濾波方法處理結(jié)果Fig.7 Real seismic records and processing results of different inverse Q filtering methods

圖8 圖7 中局部地震記錄的放大顯示Fig.8 Enlarged display of the local seismic records in Fig.7

圖9 圖8 中相應(yīng)地震記錄的歸一化振幅譜Fig.9 Normalized amplitude spectra of the corresponding seismic records in Fig.8

4 結(jié)論

（1）基于時(shí)頻域信噪比的自適應(yīng)增益限反Q濾波方法利用時(shí)頻域信噪比計(jì)算地震記錄有效頻帶及相應(yīng)的截止頻率，對地震記錄進(jìn)行時(shí)變增益的振幅補(bǔ)償，有效壓制了部分高頻噪聲能量，明顯改善了反Q濾波后地震記錄的信噪比和分辨率。

（2）本次研究方法的關(guān)鍵在于如何有效選取時(shí)頻域信噪比的閾值C0，進(jìn)而估算有效頻帶和截止頻率，其會直接影響自適應(yīng)增益限反Q濾波的結(jié)果。

（3）本次研究方法在一定的信噪比限度內(nèi)可取得較好的應(yīng)用效果，隨著信噪比的降低，其對噪聲的壓制效果逐漸變差。