基于改進(jìn)的廣義S變換的海洋地震資料隨機噪音壓制

尉佳，岳龍，楊睿，徐清風(fēng)，李志強，劉云

1.自然資源部天然氣水合物重點實驗室，青島海洋地質(zhì)研究所，青島 266071 2.山東省地震局青島地震臺，青島 266001

1 研究背景

S變換[1]是由美國的地球物理學(xué)家Stockwell等[2-3]在短時傅里葉變換和小波變換的基礎(chǔ)上提出的。與短時傅里葉變換和小波變換相比，S變換既解決了短時傅里葉變換時窗固定的問題，同時具有小波變換的多分辨率分析。由于S變換所具有的獨特優(yōu)勢，其在信號分析和處理領(lǐng)域受到廣泛關(guān)注。由于高斯時窗和頻率的函數(shù)關(guān)系較為單一，降低了S變換在實際應(yīng)用中的靈活性，因此很多學(xué)者從不同角度對S變換進(jìn)行了推廣和改進(jìn)得到廣義S變換[4-8]，使得廣義S變換更加適用于地震資料處理、油氣檢測[9-11]等方面的信號處理。

廣義S變換作為時頻域的濾波工具，能夠壓制噪音從而提高地震資料信噪比[12]。趙淑紅等[13]首先在頻域和時域確定干擾波范圍，再利用廣義S變換處理VSP數(shù)據(jù)；李雪英等[14]考慮到地震信號有效頻帶隨時間逐漸變窄，在時頻域采用交互式的方法確定高頻噪音切除邊界；王云專等采用時頻譜疊加的方法計算共中心點道集中的濾波器進(jìn)而壓制隨機干擾[15]。陳學(xué)華等引入λ和p兩個參數(shù)控制高斯窗函數(shù)的變化趨勢，在應(yīng)用過程中通過試驗確定最優(yōu)參數(shù)[5]。本文對陳學(xué)華[5]的廣義S變換方法進(jìn)行了改進(jìn)，使其兼顧對高低頻端信號的分析能力，能夠更加適用于隨機噪音壓制；其次由于標(biāo)準(zhǔn)S逆變換存在能量泄露，通過修正S逆變換[16-19]，以消除濾波噪音。在壓制噪音時，首先計算信號瞬時信噪比，根據(jù)瞬時信噪比的高低結(jié)合不同處理策略實現(xiàn)噪音壓制。通過處理合成數(shù)據(jù)和實際資料，驗證改進(jìn)廣義S正反變換在隨機噪音壓制方面的有效性。

2 方法原理

2.1 改進(jìn)廣義S變換

信號u(t)的標(biāo)準(zhǔn)S正變換的表達(dá)式為

式中，S(τ,f)是信號u(t)的時頻譜，w(t,f)是高斯窗函數(shù)，其具體表達(dá)式為

廣義S變換就是通過調(diào)節(jié)尺度因子大小控制對應(yīng)的高斯窗函數(shù)的有效寬度。針對標(biāo)準(zhǔn)S變換高斯窗形態(tài)隨頻率變化趨勢單一進(jìn)而導(dǎo)致S變換時頻譜聚焦度低的問題，陳學(xué)華等[5]在標(biāo)準(zhǔn)S變換的基礎(chǔ)上引入了控制參數(shù)λ和p，目的為豐富尺度因子隨頻率的變化趨勢，針對不同的資料提供更多的參數(shù)組合。其改進(jìn)的尺度因子為

不同的控制參數(shù)λ和p影響尺度因子變化的快慢（圖1）。不同的λ值對尺度因子的主要影響集中在低頻端，即隨著λ值的增大，同一頻率對應(yīng)的尺度因子幅值增大。不同的p值對尺度因子的主要影響在高頻端，即隨著p值的增大，同一頻率對應(yīng)的尺度因子幅值減小。為了進(jìn)一步說明兩個參數(shù)對時變信號的聚焦作用，選取理論信號進(jìn)行效果分析。該信號是由兩個不同的寬頻帶調(diào)頻信號組成，即使得高斯窗尺度因子隨頻率的變化范圍比較大，可以有效檢驗方法的適用性。

采用不同控制參數(shù)λ和p對上述信號進(jìn)行廣義S變換，分析對應(yīng)的時頻譜聚焦程度。當(dāng)p值由小變大時，信號時頻譜在低頻端與高頻端的聚焦程度均響應(yīng)明顯。其中低頻端的聚焦程度隨p值增大而變好，高頻端隨p值增大而變差，甚至當(dāng)p＞1時，高頻端出現(xiàn)嚴(yán)重失焦（圖1a-d）。當(dāng)p值由小變大時，信號時頻譜在低頻端的聚焦程度無響應(yīng)，而高頻端的聚焦程度也變化不大（圖1c、f、g、h）。廣義S變換中參數(shù)λ和p只能取固定值，其單一性的變化趨勢不能兼顧高低頻端的信號聚焦程度。綜合以上實驗結(jié)果可知，p值對時頻域聚焦程度影響更大。改進(jìn)廣義S變換主要是更改p值變化趨勢，而固定λ=1。

圖1 理想時頻譜與不同控制參數(shù)的廣義S變換的時頻譜a.理想時頻；b.p=0.7，λ=1 對應(yīng)的時頻譜；c.p=0.9, λ=1 對應(yīng)的時頻譜；d.p=1，λ=1 標(biāo)準(zhǔn) S 變換時頻譜；e.p=1.2，λ=1 對應(yīng)的時頻譜；f.p=1，λ=1.2 對應(yīng)的時頻譜；g.p=1，λ=0.9 對應(yīng)的時頻譜；h.p=1，λ=0.6 對應(yīng)的時頻譜。Fig.1 Ideal time-frequency spectrum and time-frequency spectrum of S transform using different control factorsa.ideal time-frequency spectrum, b.p=0.7, λ=1 corresponding to time-frequency spectrum, c.p=0.9, λ=1 corresponding to time-frequency spectrum, d.p=1,λ=1 corresponding to time-frequency spectrum, e.p=1.2, λ=1 corresponding to time-frequency spectrum, f.p=1, λ=1.2 corresponding to time-frequency spectrum, g.p=1, λ=0.9 corresponding to time-frequency spectrum, h.p=1, λ=0.6 corresponding to time-frequency spectrum.

p值由小變大時，低頻端與高頻端聚焦程度響應(yīng)不同。高頻端信號需要小的p值，而低頻端信號需要大的p值。為實現(xiàn)廣義S變換在時頻域的寬頻信號掃描，將參數(shù)p改為頻率的線性函數(shù)，即

其中，初始值a＞0，變化率b＜0。

在S變換計算過程中，計算的最小頻率為0 Hz，最大頻率為待分析信號的Nyquist頻率，在這個頻帶范圍內(nèi)規(guī)定p值線性遞減，同時調(diào)整p值的初始值a和變化率b，使p值能夠自動匹配信號，增加信號時頻分析的聚焦程度。

為了檢驗該方法的效果，以理論模擬信號為例進(jìn)行時頻分析對比。該信號是由兩個頻率變化趨勢相反的線性調(diào)頻信號和一個正弦調(diào)頻信號組成，見公式（6）。圖2是標(biāo)準(zhǔn)S變換、廣義S變換的優(yōu)選方案（p=0.8，λ=1）和改進(jìn)廣義S變換的時頻譜圖（p=0→1，λ=1）。由圖2可知，標(biāo)準(zhǔn) S 變換在頻率方向出現(xiàn)了條帶狀能量團(tuán)，模糊了真實時頻譜。廣義S變換改善了時頻譜圖，增強了時頻譜的聚焦度，而改進(jìn)的廣義S變換其高頻的線性調(diào)頻信號時頻譜聚焦度更高，正弦型調(diào)頻信號能量分布更均勻，分析效果要優(yōu)于廣義S變換。

圖2 合成信號時頻譜圖a.標(biāo)準(zhǔn) S 變換，b.廣義 S 變換（p=0.8，λ=1），c.改進(jìn)的廣義 S 變換（p=0→1，λ=1）。Fig.2 Synthetic signal time-frequency spectruma.the standard S transform, b.conventional generalized S transform (p=0.8, λ=1), c.improved generalized S transform (p=0→1, λ=1).

為了進(jìn)一步檢驗改進(jìn)的廣義S變換在地震信號分析中的適用性，選取一道添加高斯白噪聲的實際地震信號進(jìn)行分析。在高頻噪聲存在的情況下，利用上述三種方法得到時頻譜的能量分布（圖3）。由圖3中立體圖與平面圖對比可知，標(biāo)準(zhǔn)S變換對低頻部分的信號分析結(jié)果聚焦度比較好，但其高頻噪音的時頻譜的能量被過分放大，如果用這種方法在時頻域采用閾值濾波的方法進(jìn)行去噪，顯然會把部分信號當(dāng)成噪音同時也無法把高頻噪音去掉。廣義S變換的時頻譜改善了高頻噪音部分的聚焦度，但在低頻端能量團(tuán)幾乎連在一起，時間分辨率變差。改進(jìn)的廣義S變換和標(biāo)準(zhǔn)S變換在低頻端的能量分布相近，在高頻端使高頻噪音的能量分散聚焦，而且各個時頻點的能量比較穩(wěn)定，沒有被過分放大?？梢钥闯龈倪M(jìn)的廣義S變換對地震資料進(jìn)行分析或者時頻域濾波，效果顯然要優(yōu)。

圖3 含隨機噪聲的地震信號時頻譜圖a.標(biāo)準(zhǔn) S 變換，b.廣義 S 變換（p=0.8, λ=1），c.改進(jìn)的廣義 S 變換（p=0→1, λ=1）。Fig.3 Time-frequency spectrum of seismic signal with random noisea.the standard S transform, b.conventional generalized S transform (p=0.8, λ=1), c.improved generalized S transform (p=0→1, λ=1).

2.2 改進(jìn)S域時頻濾波

由于廣義S變換具有逆變換，因此，可以利用廣義S正反變換實現(xiàn)時頻濾波。其常規(guī)濾波思路是利用二維濾波函數(shù)F(τ，f)與S變換時頻譜相乘，對乘積結(jié)果沿時間軸積分，再反傅氏變換回時間域。濾波公式為：

式中，F(xiàn)(τ,f)是時頻域濾波函數(shù)。通過計算發(fā)現(xiàn)當(dāng)F(τ,f) =1 時，即不對時頻譜作任何處理可以完全重構(gòu)信號，但當(dāng)F(τ,f) ≠1 時，實際濾波結(jié)果和理想濾波結(jié)果存在差異，無法保證時頻濾波的準(zhǔn)確性[10]，如圖4所示，在時頻域?qū)⒑铣尚盘枮?00～600 ms的時頻譜濾掉，反變換回時間域就出現(xiàn)了能量泄露。

為了解決這個問題，需要修正變換公式得到精確濾波的改進(jìn)S域時頻濾波公式。

式中，u(t)是待分析信號，w(τ-t,f)是窗函數(shù)。

根據(jù)高斯窗函數(shù)的表達(dá)式，引入一個新的變量x(τ,t)。

對式（8）針對變量t做傅里葉變換：

對比式（7）和式（9）可得

對式（10）針對變量f作反傅里葉變換：

由式（8）可知，當(dāng)τ=t時，u(t)和x(τ,t)之間具有對應(yīng)關(guān)系：

因此，令τ=t代入式（11）可得：

對u(t)進(jìn)一步展開：

式中，?代表褶積，函數(shù)m(t)的表達(dá)式：

由式（14）可知，計算得到的結(jié)果是原始信號和信號m(t)的褶積，因此在時頻域濾波結(jié)果的基礎(chǔ)上需要再做一步濾波修正處理，最終的濾波表達(dá)式為：

式 中 ，F(xiàn)(t,f) 是 時 頻 域 濾 波 函 數(shù) ，Ufilt2(f)是ufilt2(t)的頻譜，M(f)是m(t)的頻譜，IFT表示反傅里葉變換。式（16）和式（17）是最終的濾波公式。

為了檢驗對S逆變換的改進(jìn)效果，分別用標(biāo)準(zhǔn)逆變換、改進(jìn)逆變換將圖4b中的時頻譜逆變換，并計算與準(zhǔn)確濾波結(jié)果的差值，如圖5所示。對比結(jié)果顯示，改進(jìn)逆變換消除了濾波區(qū)信號的畸變，避免產(chǎn)生濾波噪音。

圖4 標(biāo)準(zhǔn) S 逆變換時頻濾波效果圖a.合成信號，b.合成信號時頻域濾波，c.時間域濾波效果。Fig.4 Effect of time and frequency filtering of standard S inverse transforma.synthetic signal, b.time and frequency filtering of synthetic signal, c.filtering effect of time field.

圖5 準(zhǔn)確濾波結(jié)果和不同濾波方法計算的差值a.準(zhǔn)確濾波結(jié)果，b.標(biāo)準(zhǔn)逆變換計算差值（紅色）和改進(jìn)逆變換（藍(lán)色）計算差值。Fig.5 Accurate filtering result and difference calculated by different filtering methodsa.accurate filtering result, b.difference calculated by standard S inverse transform (red) and difference calculated by improved S inverse transform (blue).

基于上述改進(jìn)廣義S變換和修正的S逆變換，可以將隨機噪音壓制概括為3步：① 對含噪信號進(jìn)行改進(jìn)廣義S變換得到信號的S譜；② 在S譜的基礎(chǔ)上利用譜比法[20]計算信號的瞬時信噪比；③ 根據(jù)瞬時信噪比的不同，采用指數(shù)衰減和閾值濾波兩種策略壓制隨機噪音。具體處理步驟如下：

（1）計算信號的S變換時頻譜，利用該時頻譜計算瞬時信噪比和閾值。閾值的計算思路是：計算初至前兩個頻率段內(nèi)噪音的時頻譜的平均值，作為初始閾值，具體的頻率段可以分成有效頻帶內(nèi)和有效頻帶外兩段，計算兩個閾值。由于高頻隨機干擾的時頻譜幅值的均值隨時間變化不大，而相對淺層數(shù)據(jù)深層數(shù)據(jù)的高頻成分衰減快，能量弱，高頻噪音能量占優(yōu)，因此在無法利用初至前記錄計算閾值時也可利用深層數(shù)據(jù)較高頻段的時頻譜計算初始閾值。

式中，Emean(t,f)是閾值，S(t,f)是時頻譜，T是參與計算時間點個數(shù)，F(xiàn)是參與計算頻率點個數(shù)，t1到t2是參與計算的時間段，f1到f2是參與計算的頻率段。

（2）實際地震數(shù)據(jù)需要考慮地層的吸收衰減，且高頻衰減更快，因此深部有效頻帶變窄，有效信號能量弱[21]，此時如果將深部地震記錄按照閾值處理，會損失很多有效信號，因此為了保留一部分有效信號以供后續(xù)處理和進(jìn)一步使用，此時對瞬時信噪比低的數(shù)據(jù)段采用指數(shù)衰減的方法處理，即對信噪比小于1的時間點，對該時間點大于有效頻帶高頻端的S譜的幅值進(jìn)行一個指數(shù)衰減的加權(quán)，即將S(t,f)乘以一個以f為變量的指數(shù)衰減函數(shù)H(f)，如式（19）：

式中，μ是一個衰減系數(shù)，f是頻率，μ越大對高頻衰減越小，反之對高頻衰減越大。

（3）當(dāng)瞬時信噪比大于1時，有效信號能量占優(yōu)，且隨機噪音的時頻譜分布在整個頻帶上，對該時間點所有S譜采用閾值濾波法，當(dāng)某頻率點的幅值小于r倍的Emean(t,f)的時候，定義該點的S譜為噪音，將其充零，當(dāng)其幅值大于r倍的Emean(t,f)的時候，將其時頻譜的幅值減去r倍的Emean(t,f)，并將處理之后的時頻譜反變換回時間域。

在實際應(yīng)用中，應(yīng)根據(jù)實際地震資料的情況，進(jìn)行反復(fù)試驗，選取合適的μ和r值，既能壓制隨機噪音，又能保護(hù)高頻端的有效信號，為數(shù)據(jù)拓頻留有空間。

3 合成信號試驗

為了檢驗方法的有效性，首先利用零相位雷克子波模擬48道經(jīng)地層吸收衰減的地震記錄，加入隨機白噪聲后，利用改進(jìn)廣義S變換對地震記錄壓制隨機噪音（圖6）。抽取單道地震模擬數(shù)據(jù)對比，能夠清晰地看出改進(jìn)廣義S變換壓制隨機噪聲的效果，尤其深部低信噪比信號也得到有效的恢復(fù)（圖7）。分析時頻域中濾波效果，可看到有效信號時頻譜得以保留，隨機噪音時頻譜得到有效壓制（圖8）。

圖6 含噪地震記錄去噪效果圖a.48 道模擬地震記錄， b.含白噪聲的模擬地震記錄， c.去除噪聲的模擬地震記錄。Fig.6 Denoising effect for seismic records with noisea.48 channels simulated seismic record, b.simulated seismic record with white noise, c.simulated seismic record after denoising.

圖7 單道信號對比圖a.加白噪聲后與未加噪聲單道信號對比， b.去噪后與未加噪聲單道信號對比。Fig.7 Comparison of single channel signalsa.comparison of single channel signals with noise and those without noise, b.comparison of single channel signals without noise and those after noise reduction.

圖8 時頻域濾波效果對比a.不含白噪聲模擬地震道的時頻譜，b.含白噪聲模擬地震道的時頻譜，c.去除噪聲模擬地震道的時頻譜。Fig.8 Comparison of filtering effects in time-frequency domaina.time and frequency spectrum of simulated seismic channel without noise, b.time and frequency spectrum of simulated seismic channel with white noise,c.time and frequency spectrum of simulated seismic channel after noise reduction.

4 實際資料處理

為了檢驗本文改進(jìn)的廣義S變換方法在實際資料中的應(yīng)用效果，選擇高分辨率海洋二維地震疊前剖面進(jìn)行去噪處理，并將其與傳統(tǒng)廣義S變換濾波去噪法進(jìn)行對比。本文首先通過單炮初至前的噪聲信號確定去噪閾值，再利用改進(jìn)的廣義S變換進(jìn)行時頻濾波。圖9為地震疊前單炮剖面。傳統(tǒng)廣義 S 變換參數(shù)經(jīng)過對比后確定為（p=0.6,λ=1），改進(jìn)后的廣義 S 變換應(yīng)用參數(shù)（p=0.6→1,λ=1）。圖10是兩種方法壓制噪聲效果圖。從整體上可以看出，兩種廣義S變換時頻域去噪均能達(dá)到一定效果，隨機噪聲得到合理去除。但改進(jìn)的廣義S變換在遠(yuǎn)偏移距處的弱信號去噪效果更好，深部反射同相軸連續(xù)性更優(yōu)，有效信號的損失較小，對弱信號起到保護(hù)作用。

圖9 地震疊前單炮剖面Fig.9 Seismic pre-stack single shot profile

圖10 兩種方法壓制噪聲效果對比a.改進(jìn)廣義S變換壓噪后剖面，b.傳統(tǒng)廣義S變換壓噪后剖面，c.改進(jìn)廣義S變換去除的噪聲，d.傳統(tǒng)廣義S變換去除的噪聲。Fig.10 Comparison of noise suppression effects between two methodsa.denoised profile by improved generalized S transform, b.denoised profile by traditional generalized S transform, c.noise removal by improved generalized S transform, d.noise removal by generalized S transform.

5 結(jié)論

本文針對廣義S變換高斯窗參數(shù)λ和p固定、逆變換存在信號泄露的問題，提出了改進(jìn)的思路和計算過程，并將其應(yīng)用于隨機噪音壓制。通過對合成信號和實際資料的應(yīng)用，得到以下認(rèn)識：

（1）改進(jìn)廣義S變換時頻譜聚焦度要優(yōu)于廣義S變換，而且參數(shù)p與頻率f呈線性變化，以適應(yīng)不同的應(yīng)用目的；改進(jìn)后的逆變換可以實現(xiàn)精確濾波，消除了濾波噪音。

（2）利用單炮初至前或深部較高頻段時頻譜計算閾值，根據(jù)每個時刻時頻譜的幅值與閾值的差異壓制噪音，同時考慮到地震信號的有效頻帶和瞬時信噪比隨時間變化的情況，對瞬時信噪比較低的采用指數(shù)衰減的策略壓制噪音，可以在保護(hù)有效信號的基礎(chǔ)上實現(xiàn)去噪目的。

（3）高斯窗參數(shù)p的范圍選擇，應(yīng)當(dāng)根據(jù)高低頻端的聚焦程度而改變。通過試驗確定的參數(shù)p，既能夠有效壓制噪音又能避免弱有效信號的損失，為后續(xù)處理打好基礎(chǔ)。