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

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

1 研究背景

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

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

2 方法原理

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

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

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

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

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

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

圖1 理想時(shí)頻譜與不同控制參數(shù)的廣義S變換的時(shí)頻譜a.理想時(shí)頻；b.p=0.7，λ=1 對應(yīng)的時(shí)頻譜；c.p=0.9, λ=1 對應(yīng)的時(shí)頻譜；d.p=1，λ=1 標(biāo)準(zhǔn) S 變換時(shí)頻譜；e.p=1.2，λ=1 對應(yīng)的時(shí)頻譜；f.p=1，λ=1.2 對應(yīng)的時(shí)頻譜；g.p=1，λ=0.9 對應(yīng)的時(shí)頻譜；h.p=1，λ=0.6 對應(yīng)的時(shí)頻譜。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值由小變大時(shí)，低頻端與高頻端聚焦程度響應(yīng)不同。高頻端信號需要小的p值，而低頻端信號需要大的p值。為實(shí)現(xiàn)廣義S變換在時(shí)頻域的寬頻信號掃描，將參數(shù)p改為頻率的線性函數(shù)，即

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

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

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

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

圖3 含隨機(jī)噪聲的地震信號時(shí)頻譜圖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域時(shí)頻濾波

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖4 標(biāo)準(zhǔn) S 逆變換時(shí)頻濾波效果圖a.合成信號，b.合成信號時(shí)頻域?yàn)V波，c.時(shí)間域?yàn)V波效果。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é)果和不同濾波方法計(jì)算的差值a.準(zhǔn)確濾波結(jié)果，b.標(biāo)準(zhǔn)逆變換計(jì)算差值（紅色）和改進(jìn)逆變換（藍(lán)色）計(jì)算差值。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逆變換，可以將隨機(jī)噪音壓制概括為3步：① 對含噪信號進(jìn)行改進(jìn)廣義S變換得到信號的S譜；② 在S譜的基礎(chǔ)上利用譜比法[20]計(jì)算信號的瞬時(shí)信噪比；③ 根據(jù)瞬時(shí)信噪比的不同，采用指數(shù)衰減和閾值濾波兩種策略壓制隨機(jī)噪音。具體處理步驟如下：

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

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

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

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

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

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

3 合成信號試驗(yàn)

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

為了檢驗(yàn)本文改進(jìn)的廣義S變換方法在實(shí)際資料中的應(yīng)用效果，選擇高分辨率海洋二維地震疊前剖面進(jìn)行去噪處理，并將其與傳統(tǒng)廣義S變換濾波去噪法進(jìn)行對比。本文首先通過單炮初至前的噪聲信號確定去噪閾值，再利用改進(jìn)的廣義S變換進(jìn)行時(shí)頻濾波。圖9為地震疊前單炮剖面。傳統(tǒng)廣義 S 變換參數(shù)經(jīng)過對比后確定為（p=0.6,λ=1），改進(jìn)后的廣義 S 變換應(yīng)用參數(shù)（p=0.6→1,λ=1）。圖10是兩種方法壓制噪聲效果圖。從整體上可以看出，兩種廣義S變換時(shí)頻域去噪均能達(dá)到一定效果，隨機(jī)噪聲得到合理去除。但改進(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)的思路和計(jì)算過程，并將其應(yīng)用于隨機(jī)噪音壓制。通過對合成信號和實(shí)際資料的應(yīng)用，得到以下認(rèn)識：

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

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

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