分布式聲波傳感VSP數(shù)據(jù)中的光纜耦合噪聲局部稀疏優(yōu)化壓制方法

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

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

0 引言

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

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

Mestayer等(2011)在2009年開始使用分布式傳感系統(tǒng)(Distribute Acoustic Sensing, DAS)進(jì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ī)的干擾波及噪聲(呂公河等,2022;邵婕等,2022),還存在纜繩與測井耦合不好,容易受到振動事件而擺動,因此在采用DAS采集VSP地震數(shù)據(jù)時會受到強烈的光纜耦合噪聲影響(Yu et al.,2016).光纜耦合噪聲嚴(yán)重影響地震資料的信噪比,對VSP地震資料后續(xù)的波場分離、成像及反演造成極大障礙(Bakku et al.,2014).圖1所示為DAS系統(tǒng)實際采的集VSP地震資料,可以看到黑色箭頭所指的呈現(xiàn)“Z”字形的強相干噪聲就是典型的光纜耦合噪聲.強烈的光纜耦合噪聲降低了該地震數(shù)據(jù)的信噪比,嚴(yán)重淹沒了上行波和下行波的同相軸,對于后續(xù)的波場分離、數(shù)據(jù)成像、屬性分析等應(yīng)用產(chǎn)生了巨大的影響.因此,業(yè)界不但需要能夠有效壓制DAS系統(tǒng)光纜耦合噪聲的方法,而且還要求該方法對有效信號具有較高的保真性.

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

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

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

1 方法

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

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

s=s1+s2+n,

(1)

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

s=s1+s2,

(2)

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

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

s=Φ1x1+Φ2x2,

(3)

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

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

(4)

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

(5)

(6)

(7)

(8)

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

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

s=s1+s2+sn,

(9)

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

(10)

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

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

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

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

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

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

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

(11)

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

(12)

離散余弦變換的定義為:

k=1,2,…,N,

(13)

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

(14)

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

n=1,2,…,N.

(15)

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

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

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

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

(16)

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

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

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

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

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

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

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

si=si_s+si_d+ni,

(17)

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

(18)

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

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

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

(19)

(20)

(21)

(22)

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

(23)

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

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

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

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

(3)第三步,針對不同長度的地震信號l1的分離結(jié)果,采用(23)式計算稀疏度測量公式;

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

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

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

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

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

2 算例

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

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

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

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

(24)

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

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

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

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

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

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

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

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

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

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

表3 實測信號全局稀疏優(yōu)化和局部稀疏優(yōu)化所得的Ψ

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

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

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

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

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

3 結(jié)論

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