陳文超,陳 昕,王 偉,王曉凱,陳建友,3
(1.西安交通大學(xué),陜西西安710049;2.海洋石油勘探國家工程實(shí)驗(yàn)室,陜西西安710049;3.酒泉衛(wèi)星發(fā)射中心,甘肅酒泉732750)
對(duì)地震信號(hào)進(jìn)行稀疏表示可以有效揭示地震信號(hào)的本質(zhì)特征,有利于更為清晰和直觀地認(rèn)識(shí)和描述地震信號(hào)。在數(shù)字信號(hào)處理中,常常將信號(hào)在給定的函數(shù)(或矢量)集上進(jìn)行投影表示,把這種在變換域上表達(dá)原始信號(hào)的方式稱為信號(hào)表示或者信號(hào)變換。信號(hào)的稀疏表示就是在變換域上尋找盡量少的基函數(shù)重構(gòu)逼近原始信號(hào),從而獲得原始信號(hào)簡潔而有效的表達(dá)形式。目前常用的變換方法主要有傅里葉變換、離散余弦變換(Discrete Cosine Transform,DCT)、小波變換、曲波(Curvelet)變換等,以及聯(lián)合這些變換方法組成超完備字典的信號(hào)分解方法。而稀疏性是壓縮感知的前提,如果采集很少一部分?jǐn)?shù)據(jù)并且希望從這些少量數(shù)據(jù)中解壓出大量信息,就需要保證:①這些采集到的少量數(shù)據(jù)包含了原始信號(hào)的全局信息;②存在一種算法能夠從這些少量的數(shù)據(jù)中還原出原先的信息。因此,信號(hào)的稀疏表示與壓縮感知息息相關(guān)。
以上基于特定數(shù)學(xué)變換的信號(hào)表示方法,根據(jù)各自變換基的特性,可對(duì)相應(yīng)結(jié)構(gòu)模式的信號(hào)進(jìn)行有效表示。然而,考慮到實(shí)際地震數(shù)據(jù)的復(fù)雜性,如疊前炮集中的體波、面波、工業(yè)干擾、抽油機(jī)振動(dòng)噪聲、滑動(dòng)掃描接收方式中的諧波噪聲等,其波形結(jié)構(gòu)差異明顯。另外,疊后地震數(shù)據(jù)中穩(wěn)定沉積地層及與陸相儲(chǔ)層相關(guān)的地質(zhì)體(如點(diǎn)壩、河道等)的地震響應(yīng)信號(hào)波形結(jié)構(gòu)也存在明顯差異。這樣,在分析實(shí)際地震信號(hào)時(shí),僅選擇單一的變換方法,往往很難對(duì)地震信號(hào)實(shí)現(xiàn)稀疏分解。如果將以上多種信號(hào)表示方法聯(lián)合起來組成超完備冗余表示字典,就有可能實(shí)現(xiàn)對(duì)復(fù)雜信號(hào)的稀疏表示。基于近年來信號(hào)稀疏表示方面研究取得的進(jìn)展,相關(guān)學(xué)者發(fā)展了形態(tài)成分分析(Morphological Component Analysis,MCA)的信號(hào)分析方法[1-2],就是利用超完備的復(fù)合波形字典對(duì)信號(hào)進(jìn)行稀疏分解。這種理論認(rèn)為復(fù)雜信號(hào)可以表示成若干信號(hào)成分的線性組合,假設(shè)每種成分都能找到對(duì)其進(jìn)行稀疏表示的波形字典,而該字典對(duì)其它成分不能稀疏表示,那么復(fù)雜信號(hào)就能夠由這些波形字典聯(lián)合組成的超完備冗余字典進(jìn)行稀疏表示。另外,如果能夠根據(jù)地震數(shù)據(jù)自身的特點(diǎn),自適應(yīng)地構(gòu)造出適用于這種地震數(shù)據(jù)的變換基,也可以實(shí)現(xiàn)對(duì)地震數(shù)據(jù)的稀疏分解。OLSHAUSEN等[3]從待分析信號(hào)自身出發(fā),采用機(jī)器學(xué)習(xí)的方法自適應(yīng)地構(gòu)造出超完備字典,帶來信號(hào)表示的稀疏化。RUBINSTEIN等[4]結(jié)合數(shù)學(xué)變換,給出基于稀疏字典學(xué)習(xí)的信號(hào)稀疏表示方法,該方法提高了自適應(yīng)波形字典的普適性,降低了其對(duì)訓(xùn)練樣本數(shù)量的要求,并且使自適應(yīng)字典具有快速的前向和伴隨算子,從而使得自適應(yīng)字典方法適用于高維字典的構(gòu)造和高維信號(hào)的分析。因此,將復(fù)合型和學(xué)習(xí)型超完備波形字典的信號(hào)表示方法用于地震信號(hào)處理和解釋,有利于更好地理解地震信號(hào)包含的信息,是非常有意義的發(fā)展方向。
本文根據(jù)地震數(shù)據(jù)中波形形態(tài)特征的差異,構(gòu)造了兩種基于稀疏表示的超完備字典,根據(jù)其功能分為:體波與面波干擾分離以及強(qiáng)反射屏蔽效應(yīng)的消除,并且利用實(shí)際數(shù)據(jù)驗(yàn)證了方法的有效性。
本文根據(jù)地震信號(hào)有效波和干擾波的波形結(jié)構(gòu)特點(diǎn),構(gòu)造了用于壓制面波、工業(yè)噪聲、諧波噪聲及抽油機(jī)噪聲等干擾的一維炮集數(shù)據(jù)信噪分離模型。對(duì)于炮集數(shù)據(jù),每道信號(hào)可看作有效信號(hào)和噪聲兩種分量的線性疊加,地震記錄s可表示為:
(1)
式中:s1為有效信號(hào)分量;s2為噪聲分量;n為隨機(jī)噪聲分量。
(1)式引入了加性隨機(jī)噪聲n,從而提高了信號(hào)模型的普適性,使其不僅能夠描述受到噪聲干擾的地震數(shù)據(jù),也允許實(shí)際數(shù)據(jù)包含適當(dāng)?shù)臏y量誤差或其它雜波干擾。假設(shè)隨機(jī)噪聲分量n服從零均值的高斯白噪分布,其在每個(gè)位置的取值都是一個(gè)N(0,σ2)隨機(jī)過程。
(1)式表示在噪聲干擾下的一種多源混合信號(hào)恢復(fù)重建和分離問題,這里的兩種源信號(hào)分別指有效波信號(hào)s1和噪聲s2。MCA的目標(biāo)便是分別提取出s1,s2兩種成分。根據(jù)s1和s2具有不同的形態(tài)特征,MCA理論假設(shè)s1和s2可以分別由字典Φ1和Φ2有效地稀疏表示,但是用Φ2稀疏表示s1和用Φ1稀疏表示s2時(shí)稀疏性差,那么有效波信號(hào)和干擾噪聲可以在Φ1和Φ2聯(lián)合的超完備波形字典中得到稀疏表示,從而通過求解如下的最優(yōu)化問題實(shí)現(xiàn)兩種信號(hào)分量的恢復(fù)和分離:
(2)
式中:α1為有效波信號(hào)在字典Φ1中的向量表示形式;α2為噪聲在字典Φ2中的向量表示形式;ε為信號(hào)重建的誤差門限,使得該優(yōu)化問題在面對(duì)噪聲干擾及數(shù)據(jù)缺失時(shí)仍然能夠取得合理的稀疏解。
BRUCE等[5]提出了分塊坐標(biāo)松弛算法(Block-Coordinate-Relaxation,BCR)對(duì)(2)式進(jìn)行求解,其核心算法為通過每次迭代輪流對(duì)系數(shù)向量α1和α2進(jìn)行更新。將(2)式中的優(yōu)化問題表示為:
(3)
(4)
(5)
(6)
如此反復(fù)迭代,求解兩部分的稀疏系數(shù),最終從混合信號(hào)s中將具有不同波形特征的信號(hào)s1和s2分量進(jìn)行分離。
另外,根據(jù)疊后地震資料中穩(wěn)定沉積地層及與儲(chǔ)層相關(guān)的地質(zhì)體(河道、點(diǎn)壩等)對(duì)應(yīng)的地震響應(yīng)波形結(jié)構(gòu)特點(diǎn),我們同時(shí)構(gòu)造了用于提取高精度陸相沉積地震響應(yīng)信號(hào)的疊后地震資料模型。其數(shù)學(xué)表達(dá)形式與(1)式類似。
面波是陸上地震記錄中一種具有明顯的高振幅和頻散特性的規(guī)則干擾波,具有能量強(qiáng)、頻率低、視速度低等特點(diǎn),在單炮地震記錄中一般呈扇形分布,遮
蓋了絕大部分的體波信息。合理地壓制面波干擾,有利于提高地震資料的信噪比。常用的壓制面波干擾的方法包括帶通濾波、高通濾波和多通道濾波以及基于視速度特征差異的F-K濾波、沿時(shí)間-偏移距的加窗F-K濾波等[6-9]?;隗w波及面波在變換域分布區(qū)域的不同,小波變換、Curvelet變換、Radon變換等方法也被成功用于面波壓制。
考慮到實(shí)際地震資料非常復(fù)雜,僅選擇單一的變換方法往往很難對(duì)地震信號(hào)實(shí)現(xiàn)稀疏表示。如果把以上多種信號(hào)表示方法聯(lián)合起來組成超完備冗余表示字典,就有可能實(shí)現(xiàn)對(duì)復(fù)雜信號(hào)的稀疏表示。陳文超等[10]提出了聯(lián)合連續(xù)小波變換與局部余弦變換的超完備字典,成功實(shí)現(xiàn)了面波及體波分量的分離。TRAD等[11-12]提出將線性Radon變換和雙曲Radon變換聯(lián)合構(gòu)成一種混合Radon變換,該方法可對(duì)符合模型要求的地震資料取得較好的波場分離效果。YARHAM等[13]將兩種波場分量的初步估計(jì)用作Curvelet域中待求解真實(shí)波場分量的加權(quán),在Curvelet域建立稀疏優(yōu)化分離模型,得到較好的波場分離結(jié)果。WANG等[14]提出將一維平穩(wěn)小波變換與局部離散余弦變換聯(lián)合構(gòu)成超完備波形字典,初步探索了兩種波場分量的自適應(yīng)稀疏分解途徑。
本文根據(jù)體波和面波信號(hào)的波形特征差異,利用Q可調(diào)小波函數(shù)[15],分別構(gòu)造了能夠稀疏表示體波和面波分量的Q小波函數(shù)。針對(duì)低震蕩體波選擇低Q值小波變換作為其稀疏表示字典Φ1,針對(duì)高震蕩面波選擇高Q值小波變換作為其稀疏表示字典Φ2,建立基于雙波形字典的優(yōu)化問題,如公式(2)所示,求解該稀疏優(yōu)化問題,按照公式(5)和公式(6)反復(fù)迭代求解,得到最終的信噪分離結(jié)果。
圖1a和圖1b分別為能稀疏表示體波和面波的Q小波函數(shù)。由圖1可以看出,Q比較小的母小波函數(shù)具有寬帶特性,能稀疏表示頻帶較寬的體波信號(hào);而Q較大的小波函數(shù)能稀疏表示頻帶較窄的面波干擾。圖2為我國東部某油田實(shí)際炮集記錄分離結(jié)果。其中,圖2a為原始炮集記錄,該記錄受到了強(qiáng)面波噪聲的干擾;圖2b為采用本文方法分離出的體波信號(hào),可以看出,面波能量得到明顯壓制,中深部目的層段同相軸連續(xù)性較好;圖2c為分離出的面波干擾噪聲剖面,可以看出,面波噪聲中幾乎觀察不到有效信號(hào)殘余。也就是說,本文方法在明顯壓制面波干擾的同時(shí),對(duì)有效信號(hào)有良好的保真性能。
圖1 不同Q值的小波母函數(shù)a Q為1的小波母函數(shù); b Q為4的小波母函數(shù)
圖2 我國東部某油田實(shí)際炮集記錄分離結(jié)果a 原始炮集記錄; b 本文方法分離出的體波信號(hào); c 本文方法分離出的面波干擾
河道砂體的預(yù)測與識(shí)別是陸相儲(chǔ)層預(yù)測面臨的主要問題之一。近年來,地震資料處理、解釋方法在砂體預(yù)測與識(shí)別方面取得了顯著的效果。例如,相干體算法能清晰地刻畫河道、斷層等地質(zhì)體邊界;曲率屬性可以幫助解釋人員識(shí)別沉積地層形態(tài)特征,研究表明,曲率屬性對(duì)于描述斷層、裂縫走向及分布十分有效;頻譜成像方法利用地層調(diào)諧效應(yīng)可以增強(qiáng)對(duì)薄砂體的識(shí)別與預(yù)測能力。
上述地震處理、解釋方法應(yīng)用成功的前提條件是儲(chǔ)層結(jié)構(gòu)相對(duì)簡單,而且需要高信噪比地震數(shù)據(jù)作為基礎(chǔ)[16]。然而實(shí)際勘探開發(fā)所面臨的目標(biāo)復(fù)雜。以大慶油田的扶余油層為例,扶余油層與上部青山口組泥巖存在較明顯的波阻抗差,因而扶余油層頂面T2反射同相軸表現(xiàn)為連續(xù)的強(qiáng)反射特征。在T2強(qiáng)反射背景影響下,扶余油藏的F11油層組砂體地震響應(yīng)識(shí)別困難,儲(chǔ)層預(yù)測成功率低[17]。另外,煤層的強(qiáng)反射特征也是影響煤系儲(chǔ)層識(shí)別的主要因素[18]。目前多子波分解技術(shù)、匹配追蹤方法及EMD分解等方法被用于壓制強(qiáng)地震反射背景,提高了儲(chǔ)層識(shí)別率。
本文根據(jù)河道砂體等與儲(chǔ)層相關(guān)的地質(zhì)體空間展布范圍較小的特點(diǎn),考慮到強(qiáng)背景反射連續(xù)性好、穩(wěn)定等特征,構(gòu)造了能分別稀疏表示河道砂體等與儲(chǔ)層相關(guān)的地質(zhì)體及穩(wěn)定地層的地震響應(yīng)的超完備波形字典Φ1與Φ2,建立稀疏優(yōu)化問題,如公式(2)所示,通過公式(5)、公式(6)迭代求解,實(shí)現(xiàn)了強(qiáng)背景干擾分離,突出了砂體地震響應(yīng)信號(hào)。
圖3a為能稀疏表示穩(wěn)定沉積地層地震響應(yīng)的Curvelet變換原子,圖3b為能稀疏表示河道砂體等與儲(chǔ)層相關(guān)的地質(zhì)體地震響應(yīng)的二維平穩(wěn)小波變換原子。由圖3可以看出,Curvelet變換原子適宜于具有各向異性的曲線狀結(jié)構(gòu)信號(hào)的多尺度和多方向分析;二維平穩(wěn)小波變換原子適宜點(diǎn)狀結(jié)構(gòu)信號(hào)的多尺度分析。
圖4是我國東部某油田實(shí)際疊后地震剖面,其中黑、白線之間區(qū)域?yàn)槟繕?biāo)層。可以看出,目標(biāo)層內(nèi)砂體地震響應(yīng)完全被上覆泥巖強(qiáng)反射所掩蓋,根據(jù)原始數(shù)據(jù)很難刻畫儲(chǔ)層位置。圖5a為圖4中沿白線所示時(shí)間位置的沿層切片;圖5b為用本文方法分離出的河道砂體等與儲(chǔ)層相關(guān)的地震響應(yīng)信號(hào)沿層切片。對(duì)比圖5a和圖5b可以看出,經(jīng)過強(qiáng)屏蔽層地震響應(yīng)的壓制,河道砂體等地震響應(yīng)信號(hào)明顯凸顯,特別是圖中黑色橢圓所示位置,河道空間形態(tài)清晰。圖5c
圖3 Curvelet變換(a)與二維平穩(wěn)小波變換(b)原子對(duì)比
圖4 我國東部某油田實(shí)際地震剖面
圖5 某油田壓制強(qiáng)反射背景后的砂體地震響應(yīng)切片a 圖4所示地震剖面沿白線時(shí)間切片; b 本文方法分離出的河道砂體等與儲(chǔ)層相關(guān)的地震響應(yīng)信號(hào)沿層切片; c 本文方法分離出的屏蔽層地震響應(yīng)信號(hào)沿層切片
為本文方法分離出的屏蔽層地震響應(yīng)信號(hào)沿層切片,可以看出屏蔽層的空間大致變化情況。
本文根據(jù)地震數(shù)據(jù)中特征波形形態(tài)差異,構(gòu)造能分別稀疏表示數(shù)據(jù)組成成分的超完備字典,建立了基于超完備字典的地震信號(hào)表示方法,實(shí)現(xiàn)基于Q可調(diào)小波的面波干擾壓制和基于強(qiáng)反射背景壓制的砂體地震響應(yīng)信號(hào)提取技術(shù)。通過實(shí)際數(shù)據(jù)分析,檢驗(yàn)了方法的有效性。與目前常用的基于單變換信號(hào)表示方法相比,基于超完備字典的方法,充分考慮了地震資料的復(fù)雜性,結(jié)合了不同變換對(duì)地震數(shù)據(jù)中波形特征差異明顯的多種成分的稀疏表示性能,明顯增強(qiáng)了在變換域中解決實(shí)際問題的能力。
致謝:感謝大慶油田勘探開發(fā)研究院的陳樹民副院長、王建民副總工、王成和趙忠華專家,中國科學(xué)院地質(zhì)與地球物理研究所李幼銘研究員及海洋石油勘探國家工程實(shí)驗(yàn)室副主任高靜懷教授對(duì)本文工作的指導(dǎo)與支持。
[1] STARCK J L,ELAD M,DONOHO D L.Redundant multiscale transforms and their application for morphological component analysis[J].Advances in Imaging and Electron Physics,2004,132(4):287-348
[2] STARCK J L,ELAD M,DONOHO D L.Image decomposition via the combination of sparse representations and a variational approach[J].IEEE Transactions on Image Processing,2005,14(10):1570-1582
[3] OLSHAUSEN B A,FIELD D J.Emergence of simple-cell receptive field properties by learning a sparse code for natural images[J].Nature,1996,381(6583):607-609
[4] RUBINSTEIN R,ZIBULEVSKY M,ELAD M.Double sparsity:learning sparse dictionaries for sparse signal approximation[J].IEEE Transactions on Signal Processing,2010,58(3):1553-1564
[5] BRUCE A G,SARDY S,TSENG P.Block coordinate relaxation methods for nonparamatric signal denoising[C]∥Proceedings of International Society for Optical Engineering.USA:Society of Photo-optical Instrumentation Engineers,1998:75-86
[6] EMBREE P,BURG J P,BACKUS M M.Wide-band velocity filtering-the pie-slice process[J].Geophysics,1963,28(6):948-974
[7] GELISLI K,KARSLI H.F-K filtering using the Hartley transform[J].Journal of Seismic Exploration,1988,7(2):101-108
[8] TREITEL S,SHANKS J L,FRASIER C W.Some aspects of fan filtering[J].Geophysics,1967,32(5):789-800
[9] 劉法啟,張關(guān)泉.小波變換與F-K算法在濾波中的應(yīng)用[J].石油地球物理勘探,1996,31(6):782-791
LIU F Q,ZHANG G Q.Application of wavelet transform and F-K algorithm in filtering[J].Oil Geophysical Prospecting,1996,31(6):782-791
[10] 陳文超,王偉,高靜懷.基于地震信號(hào)波形形態(tài)差異的面波噪聲稀疏優(yōu)化分離方法[J].地球物理學(xué)報(bào),2013,56(8):2771-2782
CHEN W C,WANG W,GAO J H.Sparsity optimized separation of Ground-roll noise based on morphological diversity of seismic waveform components[J].Chinese Journal of Geophysics,2013,56(8):2771-2782
[11] TRAD D,SACCHI M,ULRYCH T J.A hybrid linear-hyperbolic Radon transform[J].Journal of Seismic Exploration,2001,9(4):303-318
[12] TRAD D,ULRYCH T J,SACCHI M.Latest views of the sparse radon transform[J].Geophysics,2003,68(2):386-399
[13] YARHAM C,BOENIGER U,HERRMANN F J.Curvelet-based ground roll removal[J].Expanded Abstracts of 76thAnnual Internat SEG Mtg,2006:2777-2781
[14] WANG W,GAO J H,CHEN W C,et al.Data adaptive ground-roll attenuation via sparsity promotion[J].Journal of Applied Geophysics,2012,83(1):19-28
[15] SELESNICK I W.Wavelet transform with tunable Q-factor[J].IEEE Transactions on Signal Processing,2011,59(8):3560-3575
[16] DE MATOS M C,YENUGU M,ANGELO S M,et al.Integrated seismic texture segmentation and cluster analysis applied to channel delineation and chert reservoir characterization[J].Geophysics,2011,76(5):11-21
[17] GUAN X W,XIE C L,FENG X Y,et al.Application of wavelet decomposition technology in thin sand identification under strong reflection[C].SPG/SEG Beijing 2016 International Geophysical Conference Technical Program Expanded Abstracts.Beijing:SPG&SEG,2016:417-420
[18] ZHANG Z J,ZHANG J H,GONG X N,et al.The method research of stripping strong reflection of coal seams[J].Expanded Abstracts of 77thEAGE Conference and Exhibition,2015:1-5