吳杰, 侯秦龍, 孫甲慶, 蘇勤,3, 王曉凱*, 寇龍江

1 中國(guó)石油天然氣股份有限公司勘探開(kāi)發(fā)研究院西北分院, 蘭州 7300202 西安交通大學(xué)信息與通信工程學(xué)院, 西安 7100493 電子科技大學(xué)環(huán)境與資源學(xué)院, 成都 611731

0 引言

隨著油氣勘探技術(shù)的不斷革新,勘探區(qū)域轉(zhuǎn)向更深層和地質(zhì)環(huán)境更加復(fù)雜區(qū)域.這些客觀因素不但對(duì)地震勘探帶來(lái)更高的挑戰(zhàn),而且使得采集到的地震資料包含更加復(fù)雜的噪聲,這對(duì)后續(xù)的地震數(shù)據(jù)處理解釋造成直接地影響.因此在地震勘探中,壓制噪聲是地震資料處理中的關(guān)鍵環(huán)節(jié),獲得高信噪比、高保真度和高分辨率地震資料,以便能夠更加有效地利用地震數(shù)據(jù),進(jìn)一步滿(mǎn)足復(fù)雜地區(qū)、復(fù)雜油氣藏勘探的需求.

傳統(tǒng)的傳感器具有易受電磁干擾、易腐蝕、在高溫高壓下效果差等缺點(diǎn).相較于傳統(tǒng)地震傳感器,光纖具有成本低廉、抗電磁干擾、耐腐蝕、耐高溫、耐高壓等優(yōu)點(diǎn).此外,光纖傳感能夠?qū)崿F(xiàn)信號(hào)獲取與傳輸一體化、空間采樣間隔小并且能夠準(zhǔn)確對(duì)光纖上振動(dòng)、溫度等信息進(jìn)行采集與監(jiān)測(cè)(馬國(guó)旗等,2020).因此,在鐵路、公路和橋梁的應(yīng)力檢測(cè)及溫度檢測(cè)方面,光纖傳感技術(shù)獲得了大量使用,并且基于此也開(kāi)發(fā)出很多基于光纖的儀器,比如溫度傳感器、光纖陀螺儀、加速度計(jì)和聲波傳感器等.光纖傳感技術(shù)在地震勘探、油氣開(kāi)發(fā)、天然地震學(xué)(張麗娜等,2020;曹衛(wèi)平等,2021;雷宇航等,2021;張衡等,2023;劉輝等,2022;林融冰等,2020,2022;江同文等,2022;張輝等,2023)及近地表探測(cè)等領(lǐng)域(Shao et al., 2022a,b;宋政宏等,2020)已引起了廣泛關(guān)注.在油氣勘探的VSP地震數(shù)據(jù)采集過(guò)程中,井中環(huán)境復(fù)雜且具有高溫高壓的特點(diǎn),而光纖具有耐腐蝕、耐高溫、耐高壓等優(yōu)點(diǎn),因此光纖在VSP地震數(shù)據(jù)采集及處理等過(guò)程中得到了廣泛地應(yīng)用(付建偉等,2004;邵婕等,2022;趙霏等,2022;武紹江等,2022;蔡志東等,2022).

Mestayer等(2011)在2009年開(kāi)始使用分布式傳感系統(tǒng)(Distribute Acoustic Sensing, DAS)進(jìn)行試驗(yàn)采集地震數(shù)據(jù),并于2010年得到采用DAS采集到的垂直地震剖面(Vertical Seismic Profile,VSP)地震數(shù)據(jù),已經(jīng)能夠達(dá)到良好成像和屬性提取質(zhì)量.相較于傳統(tǒng)采集方式,DAS系統(tǒng)獲取到的VSP地震數(shù)據(jù)質(zhì)量較高,因此其在工業(yè)上已經(jīng)獲得了廣泛地應(yīng)用(Mateeva et al.,2012;董新桐等, 2021).然而,在采用DAS系統(tǒng)進(jìn)行VSP數(shù)據(jù)采集的施工過(guò)程中,除了常規(guī)的干擾波及噪聲(呂公河等,2022;邵婕等,2022),還存在纜繩與測(cè)井耦合不好,容易受到振動(dòng)事件而擺動(dòng),因此在采用DAS采集VSP地震數(shù)據(jù)時(shí)會(huì)受到強(qiáng)烈的光纜耦合噪聲影響(Yu et al.,2016).光纜耦合噪聲嚴(yán)重影響地震資料的信噪比,對(duì)VSP地震資料后續(xù)的波場(chǎng)分離、成像及反演造成極大障礙(Bakku et al.,2014).圖1所示為DAS系統(tǒng)實(shí)際采的集VSP地震資料,可以看到黑色箭頭所指的呈現(xiàn)“Z”字形的強(qiáng)相干噪聲就是典型的光纜耦合噪聲.強(qiáng)烈的光纜耦合噪聲降低了該地震數(shù)據(jù)的信噪比,嚴(yán)重淹沒(méi)了上行波和下行波的同相軸,對(duì)于后續(xù)的波場(chǎng)分離、數(shù)據(jù)成像、屬性分析等應(yīng)用產(chǎn)生了巨大的影響.因此,業(yè)界不但需要能夠有效壓制DAS系統(tǒng)光纜耦合噪聲的方法,而且還要求該方法對(duì)有效信號(hào)具有較高的保真性.

圖1 DAS系統(tǒng)采集到的實(shí)際VSP地震數(shù)據(jù)

為提高DAS系統(tǒng)采集VSP地震資料的信噪比,諸多學(xué)者對(duì)此展開(kāi)了研究.針對(duì)DAS系統(tǒng)采集VSP地震資料中的隨機(jī)噪聲,近年來(lái)邵婕等學(xué)者利用深度學(xué)習(xí)等方法進(jìn)行了壓制(邵婕等,2022;Zhao et al.,2022),然而針對(duì)光纜耦合噪聲這種相干噪聲的壓制,仍需要進(jìn)一步的研究.Bakku等(2014)提出使用中值濾波來(lái)抑制光纜耦合噪聲干擾,但這種方法會(huì)損壞有效信號(hào).蔡志東等(2015)提出了一種基于模擬噪聲減去法來(lái)壓制光纜耦合噪聲.這種方法首先根據(jù)采集到的地震數(shù)據(jù),確定可能會(huì)受到光纜耦合噪聲干擾的區(qū)域,然后根據(jù)在確定的區(qū)域里所讀取出噪聲周期數(shù)以及相應(yīng)的時(shí)間長(zhǎng)度,得到時(shí)距曲線繼而得到振幅衰減曲線,然后通過(guò)褶積運(yùn)算得到擬合的光纜耦合噪聲數(shù)據(jù).得到噪聲數(shù)據(jù)后,通過(guò)原始數(shù)據(jù)減去光纜耦合噪聲,得到壓制光纜耦合噪聲后的地震數(shù)據(jù).然而,這種方法需要識(shí)別含有光纜耦合噪聲的區(qū)域,對(duì)檢測(cè)噪聲區(qū)域要求較高,有時(shí)會(huì)檢測(cè)不到受到較弱噪聲影響的區(qū)域.此外,VSP地震數(shù)據(jù)受到光纜耦合噪聲的影響是復(fù)雜的,采集到的大規(guī)模數(shù)據(jù)中噪聲分布與強(qiáng)弱不可預(yù)測(cè),因此會(huì)導(dǎo)致噪聲壓制不徹底或者對(duì)有效信號(hào)造成損傷.Huang等(2017)提出了一種擬合不同深度的DAS耦合噪聲并減去耦合噪聲的方法.由于DAS耦合噪聲比較復(fù)雜,這種方法對(duì)噪聲的抑制不是很徹底.Chen等(2019)提出了一種全局稀疏優(yōu)化的方法來(lái)抑制DAS耦合噪聲.此方法以形態(tài)成分分析(Morphological Component Analysis,MCA)為基礎(chǔ),通過(guò)研究DAS系統(tǒng)采集VSP地震數(shù)據(jù)中有效信號(hào)與光纜耦合噪聲在頻域與時(shí)頻域形態(tài)特征的差異,分別選擇能夠稀疏表示有效信號(hào)與光纜耦合噪聲的變換字典以構(gòu)成超完備字典,然后利用分塊坐標(biāo)松弛法(Block Coordinate Relaxation,BCR)迭代分離有效信號(hào)與光纜耦合噪聲(陳建友,2018;Chen et al.,2019).雖然全局稀疏優(yōu)化方法在處理實(shí)際數(shù)據(jù)時(shí)取得了應(yīng)用效果,但該方法對(duì)初至波處的DAS耦合噪聲的抑制效果不佳,此外,由于全局稀疏優(yōu)化方法所采用的全局參數(shù)未考慮有效信號(hào)和光纜耦合噪聲的時(shí)變特性,進(jìn)一步影響了光纜耦合噪聲的壓制效果.

本文對(duì)DAS系統(tǒng)中光纜耦合噪聲的特點(diǎn)進(jìn)行了研究,分析了有效信號(hào)和光纜耦合噪聲的時(shí)變特性,提出了一種VSP光纜耦合噪聲的局部稀疏優(yōu)化壓制方法.該方法采用分窗策略,利用局部形態(tài)成分分析方法來(lái)優(yōu)化得到最優(yōu)窗長(zhǎng)度,在局部數(shù)據(jù)上進(jìn)行分離,最終實(shí)現(xiàn)基于局部形態(tài)成分分析的光纜耦合噪聲局部稀疏優(yōu)化壓制方法.

1 方法

1.1 形態(tài)成分分析

對(duì)圖像或者地震信號(hào)的多種分量進(jìn)行分離的一種有效方案是基于信號(hào)稀疏表示理論的形態(tài)成分分析(MCA)方法.MCA理論假設(shè)對(duì)某待分析復(fù)雜信號(hào)s∈RN由多種具有不同形態(tài)特征的分量和噪聲相加得到.針對(duì)一般情況,通常將待分析的復(fù)雜信號(hào)建模成兩種具有不同形態(tài)特征的分量s1、s2與隨機(jī)噪聲n之和:

s=s1+s2+n,

(1)

在實(shí)際使用中,隨機(jī)噪聲常被忽略,因此式(1)可以進(jìn)一步簡(jiǎn)化如下簡(jiǎn)單模型:

s=s1+s2,

(2)

由于兩個(gè)分量s1及s2具有不同的形態(tài)特征,因此MCA理論假設(shè)對(duì)兩種信號(hào)分量s1、s2各存在稀疏表示字典Φ1及Φ2并做如下假設(shè):信號(hào)分量s1可以由字典Φ1稀疏表示,信號(hào)分量s2可以由字典Φ2稀疏表示,而信號(hào)分量s1無(wú)法由字典Φ2稀疏表示,同樣信號(hào)分量s2無(wú)法由字典Φ1稀疏表示.由于對(duì)于信號(hào)的每個(gè)分量都有對(duì)應(yīng)的字典稀疏表示,但是該稀疏字典對(duì)其他信號(hào)分量均不能進(jìn)行有效的稀疏表示.因此在形態(tài)成分分析理論中,選擇合適的字典對(duì)于信號(hào)分量的分離起到重要作用.

基于以上假設(shè)條件,將稀疏表示分量s1的字典Φ1和稀疏表示分量s2的字典Φ2組合構(gòu)成超完備字典并用于稀疏表示s.基于超完備字典,可通過(guò)求解如下優(yōu)化問(wèn)題得到信號(hào)s的稀疏表示:

s=Φ1x1+Φ2x2,

(3)

對(duì)于式(3)所示優(yōu)化問(wèn)題進(jìn)行求解過(guò)程中,由于存在L0范數(shù)項(xiàng),式(3)的優(yōu)化問(wèn)題為典型的非凸問(wèn)題,求解困難;其次信號(hào)受到噪聲影響,對(duì)其字典的匹配更加困難.因此,對(duì)于式(3)所描述的稀疏優(yōu)化問(wèn)題,為適當(dāng)松弛等式約束條件可將L0范數(shù)轉(zhuǎn)化為L(zhǎng)1范數(shù),并將式(3)由有約束優(yōu)化轉(zhuǎn)化為無(wú)約束條件,這樣可使求解困難的非凸問(wèn)題變得可解:

+λ(‖x1‖1+‖x2‖1).

(4)

為了給出式(4)中所示優(yōu)化問(wèn)題的求解形式,Bruce等(1998)給出分塊坐標(biāo)松弛(Block Coordinate Relaxation, BCR)算法來(lái)求解.BCR算法的核心思想是對(duì)于x1和x2兩個(gè)部分交替迭代閾值來(lái)進(jìn)行最優(yōu)化.式(4)中最優(yōu)化問(wèn)題的代價(jià)函數(shù)可以根據(jù)x1和x2兩個(gè)部分重新描述如下:

(5)

(6)

(7)

(8)

從而實(shí)現(xiàn)從混合信號(hào)s中分離出兩種信號(hào)分量s1和s2的目標(biāo).由于波場(chǎng)分離等問(wèn)題是地震信號(hào)處理與解釋中的關(guān)鍵問(wèn)題,因此MCA在面波噪聲壓制(Wang et al., 2012;陳文超等,2013; Chen et al.,2017)、工業(yè)電干擾壓制(Xu et al.,2013)、高鐵震源信號(hào)提取(王曉凱等,2019)等方面均取得了成功應(yīng)用.

對(duì)于DAS系統(tǒng)采集到的VSP地震數(shù)據(jù),假設(shè)采集得到地震數(shù)據(jù)包含有效信號(hào)、光纜耦合噪聲和隨機(jī)噪聲的疊加:

s=s1+s2+sn,

(9)

式中,s為VSP地震數(shù)據(jù),s1為有效信號(hào),s2為需要壓制的光纜耦合噪聲,sn為隨機(jī)噪聲.針對(duì)地震信號(hào)中有效信號(hào)和噪聲表現(xiàn)出不同的形態(tài)特征,可選取不同的字典分別稀疏表示,進(jìn)而構(gòu)成一個(gè)超完備字典,實(shí)現(xiàn)對(duì)地震信號(hào)的稀疏表示.在實(shí)際求解中可通過(guò)求解優(yōu)化式(10)來(lái)實(shí)現(xiàn):

(10)

式中,x1為有效信號(hào)系數(shù)矩陣,x2為需要壓制的噪聲系數(shù)矩陣,Φ1是對(duì)有效信號(hào)s1稀疏表示效果好、但對(duì)需要壓制的噪聲s2稀疏表示效果差的字典,而Φ2是對(duì)需要壓制的噪聲s2稀疏表示效果好、但對(duì)有效信號(hào)s1稀疏表示效果差的字典,ε是信號(hào)的重建門(mén)限.

1.2 VSP光纜耦合噪聲特征分析

圖1為DAS采集到的實(shí)際VSP地震數(shù)據(jù),共2035道,采樣點(diǎn)長(zhǎng)度為2000,采樣間隔為1 ms.可以看出,該VSP地震數(shù)據(jù)中所含的光纜耦合噪聲嚴(yán)重影響上下行波.圖中箭頭所指的區(qū)域的光纜耦合噪聲十分明顯,掩蓋了有效信號(hào)信息.為進(jìn)一步分析光纜耦合噪聲與有效信號(hào)的特征,我們對(duì)實(shí)際數(shù)據(jù)中光纜耦合噪聲占優(yōu)的地震道以及有效信號(hào)占優(yōu)的地震道進(jìn)行分析.圖1數(shù)據(jù)中第887道光纜耦合噪聲比較明顯,而第440道數(shù)據(jù)幾乎不含有光纜耦合噪聲,因此分別展示這兩道數(shù)據(jù)的振幅譜與時(shí)頻譜如圖2與圖3所示.可以看出,有效信號(hào)地震道(第440道數(shù)據(jù))的振幅譜表現(xiàn)為明顯的寬頻特征,而光纜耦合噪聲(第887道數(shù)據(jù))的振幅譜上表現(xiàn)為若干個(gè)窄帶譜疊加.此外,時(shí)頻譜上能觀察到光纜耦合噪聲窄帶能量隨時(shí)間的變化.上述分析結(jié)果表明,有效信號(hào)和光纜耦合噪聲的形態(tài)特征存在明顯差異.

圖2 第887道數(shù)據(jù)的頻譜分析(a) 振幅譜; (b) 時(shí)頻譜.

圖3 第440道數(shù)據(jù)的頻譜分析(a) 振幅譜; (b) 時(shí)頻譜.

根據(jù)有效信號(hào)和光纜耦合噪聲的形態(tài)特征差異進(jìn)行選擇稀疏表示有效信號(hào)與光纜耦合噪聲的字典.從對(duì)圖2與圖3的分析可得知,DAS系統(tǒng)采集地震數(shù)據(jù)中光纜耦合噪聲表現(xiàn)為若干個(gè)窄帶譜疊加,而有效信號(hào)具有明顯的寬頻特征(尤其是淺層直達(dá)波,由于衰減較弱等原因,其帶寬可達(dá)200 Hz).這種形態(tài)特征差異提供了選擇稀疏字典的依據(jù).因此,選取連續(xù)小波變換稀疏表示有效信號(hào),選取離散余弦變換稀疏表示光纜耦合噪聲.

對(duì)于待分析地震數(shù)據(jù)道x(t),其連續(xù)小波變換定義為:

(11)

式中WTx(a,τ)為地震數(shù)據(jù)的連續(xù)小波變換系數(shù),a為尺度因子控制小波的伸縮,ψ(t)為母小波.可用式(12)計(jì)算連續(xù)小波變換反變換以重構(gòu)信號(hào):

(12)

離散余弦變換的定義為:

k=1,2,…,N,

(13)

式中,x[n]表示待分析地震數(shù)據(jù),CTx(k)表示離散余弦變換系數(shù),N表示采樣長(zhǎng)度.另外ω(k)的定義為:

(14)

離散余弦變換的反變換為:

n=1,2,…,N.

(15)

此外,離散余弦變換是針對(duì)實(shí)信號(hào)所定義的一種變換,因此實(shí)信號(hào)經(jīng)過(guò)離散余弦變換以后得到的還是一個(gè)實(shí)系數(shù).同時(shí)離散余弦變換還具有運(yùn)算量小、數(shù)據(jù)處理速度快等優(yōu)勢(shì).

圖4a、b分別為連續(xù)小波變換和離散余弦變換的基本原子示意圖,圖5a、b分別為有效信號(hào)和光纜耦合噪聲的波形示意圖.由圖4a與圖5a可知,連續(xù)小波變換的原子基函數(shù)與有效信號(hào)的波形更為相似;而由圖4b與圖5b可知,離散余弦變換原子基函數(shù)與光纜耦合噪聲的波形更加匹配.對(duì)于稀疏變換,本文采用如下稀疏度來(lái)衡量各字典稀疏表示信號(hào)的能力:

圖4 連續(xù)小波變換與離散余弦變換原子示意圖(a) 連續(xù)小波變換原子基函數(shù); (b) 離散余弦變換原子基函數(shù).

圖5 有效信號(hào)與光纜耦合噪聲時(shí)域波形(a) 有效信號(hào); (b) 光纜耦合噪聲.

(16)

式中,nx為系數(shù)矩陣x的長(zhǎng)度,‖x‖0為系數(shù)矩陣x非零元素的個(gè)數(shù).上述定義中,稀疏度越小表示系數(shù)矩陣中零系數(shù)越多而非零系數(shù)越少,即說(shuō)明信號(hào)的稀疏表示程度越高.

對(duì)圖1所示DAS采集的VSP地震數(shù)據(jù)中,選取幾乎只含有有效信號(hào)(第440道)及幾乎只含光纜耦合噪聲(第887道)的兩段數(shù)據(jù)分別進(jìn)行連續(xù)小波變換和離散余弦變換,兩段數(shù)據(jù)的采樣點(diǎn)長(zhǎng)度為2000,采樣間隔為1 ms,然后利用式(16)分別計(jì)算其稀疏度,稀疏度計(jì)算結(jié)果如表1所示.可以看到,連續(xù)小波變換相比離散余弦變換在表示有效信號(hào)時(shí)更加稀疏,而離散余弦變換相比于連續(xù)小波變換在表示光纜耦合噪聲時(shí)稀疏度更高.

表1 典型信號(hào)的連續(xù)小波變換和離散余弦變換的稀疏度

當(dāng)確定分別稀疏表示有效信號(hào)與光纜耦合噪聲的變換字典后,可構(gòu)建超完備字典,進(jìn)而基于MCA的相關(guān)理論并使用BCR算法可分離有效信號(hào)與光纜耦合噪聲,以實(shí)現(xiàn)對(duì)光纜耦合噪聲的壓制.全局稀疏優(yōu)化方法(Chen et al.,2019)根據(jù)有效信號(hào)與光纜耦合噪聲頻域與時(shí)頻域形態(tài)特征差異,依照表1所示的典型信號(hào)連續(xù)小波變換和離散余弦變換的稀疏度,選取不同的字典在整個(gè)時(shí)間段上分別表示有效信號(hào)和光纜耦合噪聲(即選取連續(xù)小波變換作為有效信號(hào)的稀疏表示字典,選擇離散余弦變換作為光纜耦合噪聲的稀疏表示字典),并將這兩種字典組成超完備字典,可稀疏表示DAS采集到的VSP地震數(shù)據(jù),進(jìn)而采用該超完備字典并結(jié)合BCR算法可以在整個(gè)時(shí)間持續(xù)期上壓制光纜耦合噪聲.上述全局方法并不需要光纜耦合噪聲的各窄帶譜峰的頻率位置,因此在壓制光纜耦合噪聲方面表現(xiàn)出很大的改進(jìn).然而,全局形態(tài)成分分析方法在壓制光纜耦合噪聲時(shí)未考慮有效信號(hào)和光纜耦合噪聲的時(shí)變特性,在數(shù)據(jù)的整個(gè)時(shí)間持續(xù)期上進(jìn)行優(yōu)化,影響了噪聲壓制效果,尤其在初至波附近的區(qū)域影響尤為明顯.

1.3 基于稀疏優(yōu)化的光纜耦合噪聲壓制方法

全局稀疏優(yōu)化方法(Chen et al.,2019)的閾值在每道數(shù)據(jù)中是固定的,而光纜耦合噪聲與信號(hào)的強(qiáng)度之比會(huì)隨時(shí)間而變化,因此光纜耦合噪聲的壓制效果隨時(shí)間變化:遠(yuǎn)離初至波的光纜耦合噪聲較弱,壓制較為徹底;靠近初至波的光纜耦合噪聲較強(qiáng),壓制不徹底.本文提出將每道數(shù)據(jù)進(jìn)行分段,然后基于稀疏度選擇每個(gè)時(shí)窗的最優(yōu)長(zhǎng)度,在此基礎(chǔ)之上對(duì)每個(gè)時(shí)窗內(nèi)信號(hào)進(jìn)行稀疏優(yōu)化并實(shí)現(xiàn)有效信號(hào)和光纜耦合噪聲的分離,最終形成一種VSP光纜耦合噪聲的局部稀疏優(yōu)化壓制方法,以實(shí)現(xiàn)對(duì)光纜耦合噪聲的自適應(yīng)壓制.

假設(shè)地震記錄的一道數(shù)據(jù)L由分段l1,l2,…,lm組成,第i段li∈L的數(shù)據(jù)si表示為:

si=si_s+si_d+ni,

(17)

其中si_s表示有效信號(hào),si_d表示光纜耦合噪聲,ni為隨機(jī)噪聲.在每一個(gè)分段中,使用連續(xù)小波變換稀疏表示有效信號(hào),離散余弦變換稀疏表示光纜耦合噪聲,然后在分段內(nèi)求解優(yōu)化問(wèn)題:

(18)

按照上面的方法,對(duì)每一段數(shù)據(jù)s1,s2,…,sm分離可以得到有效信號(hào)s1_s,s2_s,…,sm_s和光纜耦合噪聲s1_d,s2_d,…,sm_d.考慮到分窗處理的邊界效應(yīng),對(duì)相鄰段重疊部分的有效信號(hào)進(jìn)行平均進(jìn)而得到整道的有效信號(hào).

與全局優(yōu)化方式不同,上述方法中一個(gè)最為重要的步驟即為窗的劃分方法.合適的時(shí)窗分割方式對(duì)于局部MCA分離有效信號(hào)和光纜耦合噪聲分離非常重要.本文將VSP地震數(shù)據(jù)初至波作為開(kāi)始截取時(shí)窗長(zhǎng)度為n1的一段作為第一段l1,而分割后的第一段地震數(shù)據(jù)s1可建模為有效信號(hào)s1_s和光纜耦合噪聲s1_d及隨機(jī)噪聲的疊加,進(jìn)而可采用BCR方法求解(18)式來(lái)實(shí)現(xiàn)分離.

稀疏系數(shù)矩陣x11是有效信號(hào)s1_s通過(guò)連續(xù)小波變換Φ1有效稀疏表示獲得,稀疏系數(shù)矩陣x22是光纜耦合噪聲s1_d通過(guò)離散余弦變換Φ2有效稀疏表示獲得,因此稀疏系數(shù)矩陣x11和x22的稀疏度高:

(19)

(20)

(21)

(22)

各稀疏表示字典均能夠稀疏表示對(duì)應(yīng)的成分但不能稀疏表示其他成分.因此,系數(shù)矩陣x11和x22較為稀疏,其稀疏矩陣非零項(xiàng)少(即‖x11‖0和‖x22‖0較小);系數(shù)矩陣x12和x21不稀疏,其稀疏矩陣非零項(xiàng)較多(即‖x12‖0和‖x21‖0較大).為綜合衡量超完備字典(由連續(xù)小波變換和離散余弦變換構(gòu)成)對(duì)分段數(shù)據(jù)l1的稀疏度,本文提出將分段l1內(nèi)的稀疏表示矩陣進(jìn)行綜合,提出了如下稀疏度測(cè)量公式:

(23)

式中n1為第一個(gè)分段的長(zhǎng)度.稀疏度測(cè)量公式(23)綜合反映了連續(xù)小波變換Φ1和離散余弦變換字典Φ2稀疏表示有效信號(hào)和光纜耦合噪聲的程度.Ψ(n1)越小,表明連續(xù)小波變換Φ1更加有效稀疏表示有效信號(hào)且離散余弦變換字典Φ2更能稀疏表示光纜耦合噪聲,因此而分離效果越好.由于分段l1的長(zhǎng)度n1在一定范圍內(nèi)變化,針對(duì)分段l1的每個(gè)可能長(zhǎng)度均按照(23)式計(jì)算稀疏度測(cè)量Ψ(n1),然后選擇最小的Ψ對(duì)應(yīng)于時(shí)窗長(zhǎng)度n1作為l1最佳的時(shí)窗長(zhǎng)度,并以此最優(yōu)時(shí)窗長(zhǎng)度和初至波的到達(dá)時(shí)刻來(lái)截取信號(hào),并在分段內(nèi)采用稀疏優(yōu)化方法實(shí)現(xiàn)光纜耦合噪聲的壓制.對(duì)于后續(xù)的數(shù)據(jù),可以用相同的最佳時(shí)窗長(zhǎng)度優(yōu)化方法進(jìn)行.

將本文所提VSP光纜耦合噪聲的自適應(yīng)稀疏優(yōu)化方法總結(jié)如下:

(1)第一步,以VSP數(shù)據(jù)單道信號(hào)初至波為起點(diǎn)并按照設(shè)定的時(shí)窗長(zhǎng)度范圍截取一系列不同長(zhǎng)度的地震信號(hào)l1;

(2)第二步,針對(duì)不同長(zhǎng)度的地震信號(hào)l1,均通過(guò)求解(18)所示優(yōu)化問(wèn)題進(jìn)行光纜耦合噪聲和有效信號(hào)的分離;

(3)第三步,針對(duì)不同長(zhǎng)度的地震信號(hào)l1的分離結(jié)果,采用(23)式計(jì)算稀疏度測(cè)量公式;

(4)第四步,針對(duì)不同長(zhǎng)度的地震信號(hào)l1分離結(jié)果的稀疏度測(cè)量,選取稀疏度測(cè)量最小對(duì)應(yīng)的長(zhǎng)度為最佳分段長(zhǎng)度,并在此最佳分段內(nèi)進(jìn)行光纜耦合噪聲和有效信號(hào)的分離;

(5)第五步,以上個(gè)最佳分段的終點(diǎn)為起點(diǎn),并按照設(shè)定的時(shí)窗長(zhǎng)度范圍截取一系列不同長(zhǎng)度的地震信號(hào)lk;

(6)第六步,返回第二步對(duì)后續(xù)分段進(jìn)行分離,直到單道信號(hào)中的每一個(gè)分段均完成分離;

(7)第七步,將每個(gè)分段的分離結(jié)果進(jìn)行拼接,形成單道數(shù)據(jù)的有效信號(hào)和光纜耦合噪聲的分離,以達(dá)到光纜耦合噪聲的壓制.

在實(shí)際處理過(guò)程中,為保證離散余弦變換的頻率分辨率及連續(xù)小波變換小尺度小波的有效支撐域,時(shí)窗的大小不宜過(guò)小,通常情況下時(shí)窗長(zhǎng)度約為1 s左右;另外,為減少分窗處理的邊界效應(yīng),建議各分段之間設(shè)置少數(shù)重疊樣點(diǎn),在處理完畢后對(duì)分離結(jié)果進(jìn)行加斜坡函數(shù)并求和以減少各分段之間的邊界效應(yīng).

2 算例

2.1 合成地震數(shù)據(jù)算例

本文使用模型數(shù)據(jù)進(jìn)行實(shí)驗(yàn),驗(yàn)證本文所提方法的有效性.使用震源子波為50 Hz的Ricker子波,采樣點(diǎn)長(zhǎng)度為2000,時(shí)間采樣間隔為1 ms,按照Ganley(1981)所提的方法,使用圖6所示的地層模型可以得到圖7a所示的合成零偏移距VSP記錄.光纜耦合噪聲模型是通過(guò)式(24)給出的數(shù)據(jù)s1[t]、s2[t]分別與70 Hz、140 Hz和40 Hz、80 Hz 的正弦信號(hào)做卷積得到的(為了便于書(shū)寫(xiě),式(24)只展示了其中一段噪聲):

圖6 地層模型結(jié)構(gòu)

圖7 合成的有效信號(hào)與光纜耦合噪聲(a) 合成的有效信號(hào); (b) 合成的光纜耦合噪聲; (c) 合成的70 Hz、140 Hz光纜耦合噪聲; (d) 合成的40 Hz、80 Hz光纜耦合噪聲.

(24)

其中圖7c為s1[t]與70 Hz及140 Hz混合正弦信號(hào)卷積得到的合成光纜耦合噪聲,而圖7d為s2[t]與40 Hz及80 Hz混合正弦信號(hào)卷積得到的合成光纜耦合噪聲.將圖7c的合成光纜耦合噪聲和圖7d的合成光纜耦合噪聲模型進(jìn)行疊加得到圖7b的最終合成光纜耦合噪聲.

分別抽取合成的有效信號(hào)和合成光纜耦合噪聲第50道數(shù)據(jù)并展示時(shí)頻譜如圖8a及圖8b所示.可以看出,有效信號(hào)在時(shí)頻譜表現(xiàn)為明顯的寬頻特征,而光纜耦合噪聲的振幅譜上表現(xiàn)為若干個(gè)窄帶譜疊加,光纜耦合噪聲在70 Hz和140 Hz上沿時(shí)間方向逐漸衰減,而在40 Hz和80 Hz上持續(xù)到1 s.

圖8 第50道合成數(shù)據(jù)中有效信號(hào)和光耦耦合噪聲的時(shí)頻圖(a) 有效信號(hào)的時(shí)頻圖; (b) 光纜耦合噪聲的時(shí)頻圖.

比較圖9a所示的有效信號(hào),圖9b所示含有光纜耦合噪聲的合成VSP數(shù)據(jù)中有效信號(hào)受到很強(qiáng)的光纜耦合噪聲干擾,部分同相軸甚至被掩蓋.對(duì)圖9b中的合成地震數(shù)據(jù)中每道數(shù)據(jù)用全局稀疏優(yōu)化方法進(jìn)行分離,同時(shí)按照本文方法分為三段l1,l2和l3并采取局部稀疏優(yōu)化進(jìn)行分離(各分段時(shí)窗變化范圍為600 ms至1000 ms).將全局稀疏優(yōu)化方法分離結(jié)果也分為三段l1,l2和l3并計(jì)算每段的稀疏度測(cè)量,如表2中第二行所示;利用局部稀疏優(yōu)化方法在不同分段l1,l2和l3進(jìn)行分離,計(jì)算的每段稀疏度測(cè)量如表2中的第三行所示.可以看出與全局稀疏優(yōu)化方法相比,雖然本文所提方法需要在多個(gè)分段內(nèi)進(jìn)行優(yōu)化會(huì)導(dǎo)致計(jì)算效率低于全局方法,但每個(gè)分段的局部稀疏優(yōu)化方法Ψ值更小,即局部稀疏優(yōu)化方法能夠更加稀疏地表示有效信號(hào)和光纜耦合噪聲.使用全局稀疏優(yōu)化方法對(duì)該合成地震數(shù)據(jù)進(jìn)行光纜耦合噪聲壓制,獲得有效信號(hào)和光纜耦合噪聲如圖9c、d所示.使用本文的局部稀疏優(yōu)化方法進(jìn)行光纜耦合噪聲壓制,獲得有效信號(hào)和光纜耦合噪聲如圖9e、f所示.可以看到,兩種方法都能夠壓制光纜耦合噪聲:全局稀疏優(yōu)化方法雖然效率較高但壓噪后在靠近初至波處具有明顯的諧波噪聲殘留;局部稀疏優(yōu)化方法雖然效率較低但對(duì)光纜耦合噪聲的壓制結(jié)果更加徹底.

表2 合成信號(hào)全局稀疏優(yōu)化和本文局部稀疏優(yōu)化所得的Ψ

圖9 全局稀疏優(yōu)化和局部稀疏優(yōu)化針對(duì)合成地震記錄光纜耦合噪聲壓制對(duì)比(a) 未受光纜耦合噪聲干擾有效信號(hào); (b) 包含光纜耦合噪聲對(duì)的合成VSP數(shù)據(jù); (c) 全局稀疏優(yōu)化方法得到的有效信號(hào); (d) 全局稀疏優(yōu)化方法壓制的光纜耦合噪聲; (e) 本文方法分離的有效信號(hào); (f) 本文方法分離的光纜耦合噪聲.

為了進(jìn)一步驗(yàn)證,抽取第50道數(shù)據(jù)的分離結(jié)果并做時(shí)頻分析,結(jié)果如圖10所示.原始數(shù)據(jù)中時(shí)頻譜含有若干個(gè)窄帶譜(圖10a);全局稀疏優(yōu)化方法殘留光纜耦合噪聲比較明顯(圖10b),局部稀疏優(yōu)化壓噪方法能夠徹底壓光纜耦合噪聲(圖10d);此外,相較于全局稀疏優(yōu)化方法提取的光纜耦合噪聲時(shí)頻譜(圖10c),本文局部稀疏優(yōu)化方法提取噪聲的時(shí)頻譜中并沒(méi)有有效信號(hào)(圖10e),因此本文所提方法對(duì)有效信號(hào)具有高保真性.

圖10 全局稀疏優(yōu)化和局部稀疏優(yōu)化方法壓噪前后單道數(shù)據(jù)(第50道)時(shí)頻譜對(duì)比(a) 原始數(shù)據(jù)時(shí)頻譜; (b) 全局稀疏優(yōu)化方法得到的有效信號(hào)的時(shí)頻譜; (c) 全局稀疏優(yōu)化方法得到的光纜耦合噪聲時(shí)頻譜; (d) 本文稀疏優(yōu)化方法得到的有效信號(hào)的時(shí)頻譜; (e) 本文稀疏優(yōu)化方法得到的光纜耦合噪聲時(shí)頻譜.

2.2 實(shí)際地震數(shù)據(jù)算例

圖11所示為DAS系統(tǒng)采集到的實(shí)際地震數(shù)據(jù),共有1070道,采樣點(diǎn)長(zhǎng)度2000,采樣間隔1 ms,該實(shí)際數(shù)據(jù)受到嚴(yán)重的光纜耦合噪聲干擾,降低了該地震數(shù)據(jù)的信噪比,并嚴(yán)重淹沒(méi)了上行波和下行波的同相軸,對(duì)后續(xù)的地震資料分析以及解釋造成嚴(yán)重影響.將圖11中的實(shí)際地震資料每道地震數(shù)據(jù)按上述方法也分為三段l1,l2和l3(各分段時(shí)窗變化范圍為600 ms至1000 ms),表3給出了圖11實(shí)際地震數(shù)據(jù)不同段l1,l2和l3的全局稀疏優(yōu)化方法和局部稀疏優(yōu)化方法的稀疏度測(cè)量Ψ.可以看出與全局稀疏優(yōu)化方法相比,局部稀疏優(yōu)化方法的Ψ值更小,說(shuō)明局部稀疏優(yōu)化方法能夠更加稀疏的表示有效信號(hào)和光纜耦合噪聲.

表3 實(shí)測(cè)信號(hào)全局稀疏優(yōu)化和局部稀疏優(yōu)化所得的Ψ

首先使用全局稀疏優(yōu)化方法對(duì)該實(shí)際地震數(shù)據(jù)進(jìn)行去噪處理,獲得有效信號(hào)和光纜耦合噪聲如圖11b、c所示.接下來(lái)使用本文的局部稀疏優(yōu)化方法進(jìn)行去噪處理,獲得的去噪結(jié)果和光纜耦合噪聲如圖11d、e所示.可以看到,兩種方法都能夠壓制光纜耦合噪聲,但全局稀疏優(yōu)化方法壓噪后在靠近初至波處具有明顯的諧波噪聲殘留,而本文局部稀疏優(yōu)化方法對(duì)光纜耦合噪聲的壓制結(jié)果更加徹底且上下行波形態(tài)清晰.

圖12a為圖11a數(shù)據(jù)的放大顯示,使用全局稀疏優(yōu)化方法獲得有效信號(hào)和光纜耦合噪聲如圖12b、c所示,使用本文的局部稀疏優(yōu)化方法獲得的去噪結(jié)果和DAS噪聲如圖12d、e所示.對(duì)比圖12黃色矩形框區(qū)域可知,全局稀疏優(yōu)化方法壓噪后具有明顯的光纜耦合噪聲殘留(具有明顯的“Z”字形的強(qiáng)相干周期噪聲),而局部稀疏優(yōu)化方法對(duì)光纜耦合噪聲的壓制結(jié)果更加徹底,幾乎沒(méi)有光纜耦合噪聲殘留.同時(shí)可以看到,光纜耦合噪聲剖面中幾乎不含有效信號(hào)能量,說(shuō)明本文方法在壓制光纜耦合噪聲的同時(shí),對(duì)有效信號(hào)具有高保真性.

圖12 全局稀疏優(yōu)化和局部稀疏優(yōu)化針對(duì)實(shí)際VSP數(shù)據(jù)光纜耦合噪聲壓制的放大對(duì)比(a) 實(shí)際VSP數(shù)據(jù); (b) 全局稀疏優(yōu)化方法得到的有效信號(hào); (c) 全局稀疏優(yōu)化方法壓制的光纜耦合噪聲; (d) 本文方法分離的有效信號(hào); (e) 本文方法分離的光纜耦合噪聲.

為了進(jìn)一步驗(yàn)證,抽取第25道數(shù)據(jù)分離結(jié)果進(jìn)行時(shí)頻分析.原始數(shù)據(jù)中時(shí)頻譜含有若干個(gè)窄帶譜(圖13a),說(shuō)明該道數(shù)據(jù)存在光纜耦合噪聲.對(duì)采用全局稀疏優(yōu)化和局部稀疏方法得到的有效信號(hào)做時(shí)頻分析,時(shí)頻譜分別如圖13b、d所示,全局稀疏優(yōu)化方法殘留光纜耦合噪聲比較明顯,而本文所提局部稀疏優(yōu)化方法能夠徹底壓光纜耦合噪聲.對(duì)采用全局稀疏優(yōu)化和局部稀疏方法壓制掉的光纜耦合噪聲做時(shí)頻分析,時(shí)頻譜分別如圖13c、e所示,本文所提局部稀疏方法提取的噪聲并沒(méi)有損傷有效信號(hào),因此對(duì)有效信號(hào)具有高保真性.

圖13 全局稀疏優(yōu)化和局部稀疏優(yōu)化方法壓噪前后單道數(shù)據(jù)(第25道實(shí)際數(shù)據(jù))時(shí)頻譜對(duì)比(a) 原始數(shù)據(jù)時(shí)頻譜; (b) 全局稀疏優(yōu)化方法得到的有效信號(hào)的時(shí)頻譜; (c) 全局稀疏優(yōu)化方法得到的光纜耦合噪聲時(shí)頻譜; (d) 本文稀疏優(yōu)化方法得到的有效信號(hào)的時(shí)頻譜; (e) 本文稀疏優(yōu)化方法得到的光纜耦合噪聲時(shí)頻譜.

3 結(jié)論

本文分析DAS采集的VSP數(shù)據(jù)中有效信號(hào)和光纜耦合噪聲的形態(tài)特征,為有效信號(hào)和光纜耦合噪聲分別選取稀疏表示字典;在此基礎(chǔ)之上,本文提出了一種綜合有效信號(hào)和光纜耦合噪聲的稀疏度測(cè)量,并以此獲得單道VSP數(shù)據(jù)的最佳分段選取策略,最終提出一種基于局部稀疏優(yōu)化的光纜耦合噪聲自適應(yīng)壓制方法.通過(guò)合成地震數(shù)據(jù)和實(shí)際地震數(shù)據(jù)驗(yàn)證了本文提出的基于局部稀疏優(yōu)化光纜耦合噪聲自適應(yīng)壓制方法,能夠有效地壓制光纜耦合噪聲,同時(shí)對(duì)有效信號(hào)具有高保真性.但本文所提局部稀疏優(yōu)化方法需要進(jìn)行多次局部稀疏優(yōu)化,會(huì)導(dǎo)致其計(jì)算量高于傳統(tǒng)的全局方法,因此后續(xù)在處理海量地震數(shù)據(jù)時(shí)需采用并行計(jì)算或者深度學(xué)習(xí)類(lèi)方法來(lái)減少計(jì)算時(shí)間.