張 帥，楊潤海，張 演，張惠菊，姚 遠(yuǎn)，沈婭宏，李 平

(1.云南省地震局，云南 昆明 650224；2.開遠(yuǎn)市地震局，云南 開遠(yuǎn) 661699)

0 引言

利用地震波研究地球內(nèi)部介質(zhì)變化依賴于重復(fù)性震源，使用天然重復(fù)震源進(jìn)行介質(zhì)變化測量的時(shí)間、空間分辨率及精度都受到一定程度的限制，而大部分人工震源監(jiān)測研究主要集中于淺層地殼(王寶善等，2016)。近年來，利用氣槍重復(fù)震源研究地球深淺部結(jié)構(gòu)和介質(zhì)變化成為一種有效手段。大容量氣槍震源與其它震源相比，有著激發(fā)性能高、低頻成分豐富、可重復(fù)性好、綠色環(huán)保、安全效能高等特點(diǎn)(Wang，2016，2020；Chen，2017；陳颙等，2007；王偉濤等，2017)，應(yīng)用前景十分廣闊。

在氣槍源探測過程中，地震波速變化測量本質(zhì)是要解決數(shù)字信號處理的時(shí)延估計(jì)問題，其中信噪比是影響信號波速測量精度的關(guān)鍵(楊潤海等，2020；蔣生淼等，2017)。而實(shí)際地震記錄信號由于存在隨機(jī)干擾和能量衰減，具有一定的非平穩(wěn)性，同時(shí)氣槍震源能量有限，傳播距離較遠(yuǎn)時(shí)信號信噪比低，甚至難以辨認(rèn)(向涯等，2019；李孝賓等，2017；徐逸鶴等，2016)。在波速變化的計(jì)算過程中，為了提高信噪比，通常會對信號進(jìn)行濾波和疊加等處理(劉自鳳等，2015；葉泵等，2017；譚俊卿等，2020；Xiang，2019)，但由于部分臺站距離太遠(yuǎn)，有效信號往往湮滅在噪聲之中，利用傳統(tǒng)降噪方法處理后的信號信噪比仍然較低，導(dǎo)致波速變化的計(jì)算結(jié)果誤差較大，甚至不能反映真實(shí)的介質(zhì)變化情況。因此，運(yùn)用新的降噪方法壓制噪聲，仍是提高地震資料信噪比和增強(qiáng)信號測量精度的主要途徑。

曲波變換是一種能夠?qū)π盘栠M(jìn)行最優(yōu)稀疏化變換的方法，它是通過基函數(shù)與信號的內(nèi)積來實(shí)現(xiàn)信號的稀疏表示(Candes，Donoho，1999；Candes，2005)，具有多分辨率特性、時(shí)頻局部性、多方向性等優(yōu)點(diǎn)，克服了小波變換表示邊緣、輪廓等高維奇異時(shí)的局限(Herrmann，Hennenfent，2008；仝中飛等，2008)。由于有效信號通過曲波變換后的曲波系數(shù)是稀疏的，而噪聲信號是隨機(jī)信號，在曲波域中的系數(shù)表現(xiàn)為非稀疏性。有效信號受到噪聲影響時(shí)，在曲波域中的稀疏性降低。傳統(tǒng)的曲波降噪方法是通過設(shè)定軟硬閾值進(jìn)行濾波，大于設(shè)定閾值的系數(shù)保留，小于閾值的去除，但如果閾值設(shè)置過大，在連續(xù)處容易造成不光滑的偽影或者偽吉布斯現(xiàn)象，而閾值設(shè)置過小，噪聲壓制能力減弱，將會嚴(yán)重影響信號的降噪效果(劉磊等，2011，王文波等，2006)，因此在曲波域中引入全變差技術(shù)通過調(diào)整域中大小系數(shù)，再進(jìn)行閾值濾波，可以使原始信號的主要特征得到保持(唐剛，楊慧朱，2010)。

為了重建曲波域中由于閾值的設(shè)定導(dǎo)致系數(shù)過度缺失，造成信號關(guān)鍵細(xì)節(jié)隨機(jī)丟失的問題，本文引入壓縮感知理論。壓縮感知理論是由Candes等(2006)正式提出(Hennenfent，Herrmann，2008)，隨著近幾年的不斷發(fā)展，已經(jīng)成功應(yīng)用于多個(gè)地震學(xué)和地球物理反演問題中(張帥等，2021；唐剛，2010；姚華建，2013；Yao，2011)。壓縮感知數(shù)據(jù)重建過程涉及稀疏變換、采樣方法(測量矩陣)和重構(gòu)算法，大部分研究都是基于這3個(gè)方面開展(Yin，2014)。傳統(tǒng)地震數(shù)據(jù)重建方法大都受 Nyquist采樣理論的限制，對于不滿足采樣定理的超稀疏地震數(shù)據(jù)，重建效果不佳(郭念民等，2014)。而壓縮感知理論表明：即使是隨機(jī)缺失的數(shù)據(jù)，也可以恢復(fù)出滿足一定精度要求的完整數(shù)據(jù)(唐剛，楊慧朱，2010；孔麗云等，2012，白蘭淑等，2010)。因此利用壓縮感知信號重建方法，可以解決因閾值選取不當(dāng)導(dǎo)致信號曲波域中重要系數(shù)缺失的問題，最大程度地重構(gòu)出有效信號的系數(shù)。

如上所述，綜合兩種方法的優(yōu)勢，本文構(gòu)建了一種基于離散曲波變換的壓縮感知方法，將曲波域中的閾值濾波和全變差技術(shù)與壓縮感知正交匹配重構(gòu)算法(OMP)的稀疏促進(jìn)策略充分結(jié)合起來，對有效信號和噪聲信號進(jìn)行有效分離，從而實(shí)現(xiàn)信號的高效降噪目的。

1 研究方法原理

1.1 曲波變換

信號的稀疏表示至關(guān)重要，離散余弦變換、小波變換、曲波變換等都能稀疏表征信號，因此常常被用作壓縮感知的稀疏變換基。曲波變換具有各向異性和多尺度、多方向的諸多優(yōu)點(diǎn)，較其它2種方法稀疏效果更好(Herrmann，Hennenfent，2008)，是地震數(shù)據(jù)稀疏表達(dá)最有效的方法之一，用離散曲波變換來稀疏表征地震信號的方法原理如下：

假定輸入信號為[,]，其中≥0,<，<，基于wrapping算法核心思想：圍繞原點(diǎn)包裹，對任意區(qū)域映射到原點(diǎn)的仿射區(qū)域(楊連剛等，2019)，其具體變換流程如下：

①將[,]通過二維的傅立葉變換和FFT域采樣得到表達(dá)式為：

(1)

②每個(gè)尺度和角度與窗函數(shù)作用后得到：

(2)

③對于每一個(gè)大小為×的窗數(shù)據(jù)，圍繞其原點(diǎn)，可得：

(3)

(4)

⑤取軟硬閾值如下：

(5)

式中：是閾值。

1.2 壓縮感知理論

壓縮感知理論核心思想是如果信號是可壓縮的或在某個(gè)變換域是稀疏的，用一個(gè)與變換基不相關(guān)的觀測矩陣將高維信號投影到低維空間上，然后通過重構(gòu)算法就可以從低維信號重構(gòu)出原始高維信號(Candes，2005)。

假設(shè)氣槍信號為()，其重構(gòu)問題可以通過以下的模型進(jìn)行表示：

()=[()]

(6)

(7)

式中：是的等價(jià)或者稀疏表示。

根據(jù)缺失的位置設(shè)計(jì)一個(gè)平穩(wěn)的、變換基不相關(guān)的×維的測量矩陣，這里選用一個(gè)服從(0,1)、正態(tài)分布、大小為×高斯隨機(jī)矩陣來滿足壓縮感知理論對于隨機(jī)投影矩陣的要求。對進(jìn)行觀測得到觀測數(shù)據(jù)為：

==

(8)

式(8)為欠定方程組，有多個(gè)方程解，但具有稀疏性，可以通過不斷促進(jìn)的稀疏性進(jìn)行反演求解。這就為實(shí)現(xiàn)缺失數(shù)據(jù)的精確重構(gòu)成為可能。其中的稀疏促進(jìn)策略(Herrmann，Hennenfent，2008)如下：

(9)

模型空間限制在x的一個(gè)L-ball 中。在(x,x)平面上可以畫出目標(biāo)函數(shù)直線,而約束條件則成為平面上大小不斷變化的一個(gè)norm ball 。直線與 norm ball 首次相交的地方就是最優(yōu)解。與L范數(shù)相比,L范數(shù)更有可能得到值為0的解(圖1),使得x更具稀疏性(姚華建,2013)。

圖1 L1范數(shù)(a)和L2范數(shù)(b)求解示意圖(p表示范數(shù))Fig.1 Schematic diagram of solving L1 norm(a) and L2 norm(b)

1.3 基于曲波域全變差的壓縮感知降噪策略

氣槍信號可以表征為()=+，其中，代表有效信號，表示為噪聲信號?；谇ㄓ虻拈撝禐V波方法(CSCT)，域中小系數(shù)容易被過度去除，缺失信號關(guān)鍵細(xì)節(jié)信息，并且不連續(xù)的區(qū)域常會出現(xiàn)震蕩的偽影現(xiàn)象(唐剛，楊慧朱，2010)，影響數(shù)據(jù)的降噪資料。因此，需要首先利用壓縮感知重建策略對因閾值濾波導(dǎo)致缺失的系數(shù)進(jìn)行重建，解決信號稀疏特征的連續(xù)性問題，再通過全變差最小化技術(shù)對重建后的系數(shù)進(jìn)行調(diào)整(Newark，1998；沈維燕，2006)，達(dá)到聯(lián)合降噪的目的。

利用全變差最小技術(shù)對重建的曲波系數(shù)進(jìn)行調(diào)整主要有兩種方式：一種是調(diào)整大于設(shè)定閾值的系數(shù)：

+1=-(())

(10)

另一種是調(diào)整小于閾值設(shè)定的系數(shù)：

+1=-(())

(11)

式中：表示在小于設(shè)定閾值的曲波系數(shù)空間上的投影；表示在大于設(shè)定閾值的曲波系數(shù)空間上的投影；設(shè)置最大迭代次數(shù)和初始值，=1。()通過以下過程求解。

假設(shè)是中的有限子集，曲波域中系數(shù)的全變差定義為：

(12)

其中

(13)

對上式進(jìn)行分部積分得到系數(shù)函數(shù)的全變差為：

(14)

式中：為圖像的支撐區(qū)間；∈為圖像的坐標(biāo)向量。

基于全變差的降噪方法可通過最小化()來實(shí)現(xiàn)：

(15)

(16)

假設(shè)其滿足Neumann邊界條件，可通過投影梯度算法進(jìn)行求解：

+1=-(())

(17)

(18)

式中：>0為一個(gè)很小的正值，以解決全變差函數(shù)在某些點(diǎn)的不可微性，取步長

(19)

式中：=+1，判斷是否接近，如果是迭代結(jié)束，否則繼續(xù)求解次梯度。

基于曲波城的壓縮感知降噪方法具體的流程如圖2所示。

圖2 基于曲波域的壓縮感知降噪方法流程圖Fig.2 The compress sensing denoising method based on curve wave domain

2 數(shù)值模擬

為了驗(yàn)證方法的有效性，利用信噪比值及均方根誤差來衡量氣槍信號的降噪效果：

(20)

(21)

2.1 氣槍信號數(shù)值模擬

基于曲波域全變差的壓縮感知降噪方法對低信噪比的氣槍信號實(shí)現(xiàn)了高效穩(wěn)定的降噪效果。為了驗(yàn)證該方法應(yīng)用于氣槍信號降噪的處理效果，其中，已知的理論一維信號長度為1 400個(gè)采樣點(diǎn)，信號帶寬2～8 Hz，采樣率100 Hz，且具有兩個(gè)正弦周期(向涯等，2019)。在每道信號上隨機(jī)加入6 dB高斯白噪聲，再從100道中隨機(jī)選擇其中的9道按照線性遞增的方式再次加入不同比例的隨機(jī)噪聲，噪聲最高信號信噪比達(dá)到-15 dB，如圖3a所示，從圖中可以看出，加入噪聲的9道信號相對原始信號已經(jīng)基本看不出有效信號。

按照本文方法的處理流程，對加噪聲信號進(jìn)行降噪處理，在參數(shù)的選擇過程中，其對曲波域閾值的選擇綜合對比了3種閾值(硬閾值、軟閾值、軟硬閾值結(jié)合)的降噪效果，其中軟硬閾值結(jié)合的降噪效果較好，折衷系數(shù)選擇0.001，能夠最大限度保留信號的細(xì)節(jié)信息。在對曲波域系數(shù)進(jìn)行壓縮重建過程中，信號的缺失率不能低于30%，否則信號不能得到有效恢復(fù)。迭代次數(shù)選擇也尤為重要，次數(shù)較大會引入隨機(jī)噪聲，較小對信號的細(xì)節(jié)信號恢復(fù)效果較差，綜合嘗試后認(rèn)為選擇迭代6次較為合適。處理結(jié)果如圖3b所示，從圖中可以看到，即使信噪比很低的信號，隨機(jī)噪聲都能得到很好地壓制，與模擬信號的波形和振幅一致性仍較強(qiáng)，相位恢復(fù)精度較高，獲得較好的處理效果。

圖3 模擬信號加載不同比例噪聲(a)以及本文方法降噪結(jié)果(b)Fig.3 The original signal loaded with different proportions of noise(a)，and the noise reduction by the method in this paper(b)

為了進(jìn)一步驗(yàn)證該方法的可行性和優(yōu)勢，本文引入帶通濾波(BDF)、小波變換濾波(DWTF)、S變換模板濾波(STTF)，分別對加噪聲信號進(jìn)行處理，過程如下：用帶通濾波對信號進(jìn)行2～8 Hz帶通濾波。小波變換濾波利用小波函數(shù)對信號進(jìn)行小波3層分解后，其中閾值統(tǒng)一取1，進(jìn)行閾值濾波。S變換模板濾波對信號進(jìn)行疊加后，變換到時(shí)頻域中，通過設(shè)定時(shí)頻域區(qū)間的閾值窗口得到濾波模板，利用濾波模板再對信號進(jìn)行濾波。4種方法的處理結(jié)果如圖4所示，選擇其中一道進(jìn)行處理效果評價(jià)(圖5)。

圖4 不同濾波方法與本文方法處理效果對比Fig.4 Comparison of processing effects between various filtering methods and the method in this paper

圖5 不同濾波降噪方法處理結(jié)果的信噪比(a)和均方根誤差對比(b)Fig.5 Comparison of the results of SNR(a)and MSE(b)by different filtering and noise reduction methods

如圖5所示，通過對4種方法處理結(jié)果的信噪比和均方根誤差可以看出，隨著噪聲比例越高，傳統(tǒng)的幾種方法處理后的信號越低，均方根誤差越高。而本文的方法對不同比例的噪聲信號降噪效果明顯高于其它方法，即使信噪比很低也能達(dá)到很好的效果，并且穩(wěn)定性更強(qiáng)。從單一道的處理效果(表1)來看，BDF方法可以將高頻成分去除，但是與有效信號頻率相近的噪聲依然存在。用DWTF方法處理后雖然有效信號得到加強(qiáng)，但是高頻成分隨機(jī)噪聲依然存在，STTF方法雖然克服了前兩者的缺點(diǎn)，但是部分相位沒有得到有效恢復(fù)，造成細(xì)節(jié)缺失，而本文的方法不僅很好地壓制了隨機(jī)噪聲，而且恢復(fù)了對相位細(xì)節(jié)。

表1 不同濾波降噪方法的處理結(jié)果對比評價(jià)Tab.1 Comparison and evaluation of processing results of different filtering and noise reduction methods

2.2 波速測量精度評價(jià)

為了進(jìn)一步驗(yàn)證本文方法在波速測量中的高效穩(wěn)定性，由已知一維氣槍信號利用移動(dòng)窗口壓縮-拉伸法(Meier，2010；劉志坤，2010)生成波速變化率不超過±5%的100道二維模擬信號(圖6a)。在每道信號上隨機(jī)加入6 dB高斯白噪聲，再隨機(jī)從100道中選擇其中的9道按照線性遞增的方式再次加入不同比例的隨機(jī)噪聲(同2.1節(jié))，如圖6b所示。按照上述的處理流程分別對加噪聲信號進(jìn)行降噪處理，對含噪信號進(jìn)行2～8 Hz帶通濾波。利用小波函數(shù)對信號進(jìn)行小波3層分解后，閾值統(tǒng)一取2，進(jìn)行小波閾值濾波。對信號所有道進(jìn)行疊加后，變換到時(shí)頻域中，通過設(shè)定時(shí)頻域區(qū)間的閾值窗口，得到濾波模板，從而達(dá)到利用疊加信號進(jìn)行模板濾波的目的。

圖6 使用不同濾波降噪方法對二維波速變化率的模擬信號處理后的降噪結(jié)果對比Fig.6 Comparison of noise reduction results of analog signals with two-dimensional wave velocity change rate by different filtering and noise reduction methods

由圖6c可以看出，本文方法不僅可以很好地壓制隨機(jī)噪聲，也去除了隨機(jī)的線性干擾；DWTF雖然恢復(fù)了有效信號的大部分能量，但是對高頻成分的噪聲壓制力較弱(圖6d)；STTF方法對隨機(jī)噪聲的壓制力較強(qiáng)，但是由于模板閾值的選取存在一定的誤差，導(dǎo)致對信號細(xì)節(jié)的恢復(fù)較差，并且對加入的較強(qiáng)隨機(jī)線性干擾的壓制能力較弱，一些線性干擾依然存在(圖6e)；BDF方法只能對固定頻段的噪聲進(jìn)行壓制，高頻噪聲和線性干擾依然存在(圖6f)。

為進(jìn)一步驗(yàn)證本文方法的穩(wěn)定性和高效性，對上述幾種降噪方法處理后的信號進(jìn)行波速測量精度評價(jià)，如圖7所示。

從圖7可以看到，利用本文方法降噪后測量的信號波速變化率與模擬信號(藍(lán)線)的量級一致，波速變化曲線基本一致，且擬合度較高。用DWTF方法處理后測量的信號波速變化率雖然與模擬信號的趨勢一致，但是依然存在一定的誤差，并且在較強(qiáng)的線性干擾處存在突變，用BDF方法處理后測量的信號波速變化率由于對隨機(jī)噪聲和線性干擾壓制較弱，導(dǎo)致波速測量的精度嚴(yán)重不足，并且在較強(qiáng)的線性干擾處存在突變量。

圖7 不同濾波降噪方法處理后波速測量結(jié)果Fig.7 Wave speed measurement results after processing by different filtering and noise reduction methods

3 實(shí)際氣槍信號處理

云南賓川氣槍主動(dòng)源區(qū)域(25.4°～ 26.4°N，99.8°～101.2°E)位于由紅河斷裂、劍川—麗江斷裂和程海斷裂圍成的三角形塊體內(nèi)，該區(qū)主體在大理—麗江活斷層系(王彬等，2016)，是由正斷和左旋走滑斷裂構(gòu)成的“Z”型張扭性復(fù)合變形帶，包括麗江—大理斷裂帶和程海—賓川斷裂帶，一系列近SN向與NW向的斷裂在此交匯疊加(張?jiān)迄i等，2020)。近年來，隨著賓川主動(dòng)源的建成，研究區(qū)內(nèi)不僅建有40多個(gè)主動(dòng)源流動(dòng)地震觀測臺陣，還有多個(gè)固定地震觀測臺站，臺站密度明顯提高，并積累了大量地震觀測資料，為開展相關(guān)科學(xué)研究提供了良好的數(shù)據(jù)保障和基礎(chǔ)(Wang，2016；Jiang，2019)。但是震中距大于20 km臺站的記錄信號信噪比普遍較低(陳蒙，2014；向涯等，2017)。雖然可以利用氣槍信號的高度重復(fù)性的特點(diǎn)通過信號疊加提高信號的信噪比，但是提高單次激發(fā)氣槍信號的信噪比，不僅可以提高氣槍激發(fā)信號的利用效率，而且對提高介質(zhì)波速時(shí)空變化測量精度具有重要意義。

本文選取距離氣槍發(fā)射源81 km的53251臺站的信號，先對信號去均值、去趨勢等處理，然后再利用本文的方法對信號進(jìn)行降噪處理。實(shí)際氣槍數(shù)據(jù)采樣率為 100 Hz，由于距離較遠(yuǎn)，隨機(jī)噪聲干擾較多，尤其是信號受到很多線性干擾的影響，使得數(shù)據(jù)信噪比較低(圖8a)。本文利用不同濾波方法對實(shí)際氣槍信號進(jìn)行降噪(圖8b～e)。

從圖8b可以看出BDF方法只能對固定頻段的噪聲進(jìn)行壓制，高頻噪聲和線性干擾依然存在。

圖8 不同濾波降噪方法對實(shí)際氣槍信號的處理效果對比Fig.8 Comparison of processing effects of different filtering and noise reduction methods on airgun source singnals

DWTF方法雖然恢復(fù)了有效信號的大部分能量，但是線性干擾依然存在，對高頻成分的噪聲壓制力較弱(圖8c)。STTF方法對隨機(jī)噪聲的壓制力較強(qiáng)，但是由于模板閾值選取的不穩(wěn)定，部分有效信號能量沒有恢復(fù)，并且對加入的較強(qiáng)隨機(jī)線性干擾的壓制能力較弱，一些線性干擾依然存在(圖8d)。本文濾波方法不僅可以很好的壓制隨機(jī)噪聲，也很好地去除了隨機(jī)的線性干擾，對同相軸的連續(xù)性和一致性恢復(fù)較好(圖8e)。

為了更加清晰地對比多種方法降噪效果，本文對上述實(shí)際氣槍信號進(jìn)行1 500次疊加，并進(jìn)行2～7 Hz帶通濾波，選擇其中一道氣槍信號進(jìn)行對比分析，并將處理后的信號與疊加信號進(jìn)行互相關(guān)分析，結(jié)果如圖9所示。

圖9 不同濾波降噪方法處理結(jié)果與疊加信號對比評價(jià)Fig.9 Comparison and evaluation of processing results of different filtering noise rodution methods and superimposed signals

由圖9可見，經(jīng)BDF方法處理后高頻成分被去除，但是與有效信號頻率相近的噪聲依然存在，互相關(guān)系數(shù)僅為0.880 7。經(jīng)DWTF方法處理后雖然有效信號得到加強(qiáng)，恢復(fù)了信號的大部分能量，但是高頻成分隨機(jī)噪聲依然存在，互相關(guān)系數(shù)為0.888 2，STTF方法雖然克服了前兩者的缺點(diǎn)，互相關(guān)系數(shù)達(dá)到0.941 7，但是部分相位沒有得到有效恢復(fù)，造成細(xì)節(jié)缺失。而本文的濾波方法處理后互相關(guān)系數(shù)達(dá)到0.986 0，不僅很好地壓制了隨機(jī)噪聲，而且恢復(fù)了相位細(xì)節(jié)，取得了較好的降噪效果。

4 結(jié)論

本文根據(jù)氣槍信號的特點(diǎn)，構(gòu)建了一種基于離散曲波變換的壓縮感知方法，將曲波域中的閾值濾波和全變差技術(shù)與壓縮感知正交匹配重構(gòu)算法的稀疏促進(jìn)策略充分結(jié)合起來，并詳細(xì)介紹了該方法的基本原理和處理流程，通過不斷驗(yàn)證，實(shí)現(xiàn)了對信噪比較低氣槍信號高效的降噪處理，主要得到以下結(jié)論：

(1)相較于傳統(tǒng)的幾種濾波方法，本文構(gòu)建的濾波方法具有明顯的降噪優(yōu)勢，即使信噪比很低的信號，隨機(jī)噪聲都能得到很好地壓制，穩(wěn)定性較好，也可以有效去除線性干擾，與模擬信號和疊加信號相比，其波形和振幅一致性恢復(fù)較好，相位恢復(fù)精度較高，獲得比較好的處理效果。

(2)信號的信噪比對波速測量的精度影響較大，相較于其它濾波方法降噪后的波速測量結(jié)果，利用本文方法降噪后測量的信號波速變化率，與模擬信號相比，兩者量級一致，波速變化曲線基本一致，且擬合度較高。

(3)在曲波域系數(shù)的重構(gòu)過程中，重構(gòu)算法和迭代次數(shù)對系數(shù)的重構(gòu)精度有一定的影響，嘗試不同的重構(gòu)算法有利于獲得更好的重構(gòu)結(jié)果，盡可能恢復(fù)信號的更多細(xì)節(jié)信息。