王本鋒,陸文凱,陳小宏,王志凱

(1.同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院,上海200092;2.清華大學(xué)自動(dòng)化系,智能技術(shù)與系統(tǒng)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京100084;3.中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京102249)

地震數(shù)據(jù)采集中,受障礙物、禁采區(qū)、海上拖纜羽狀漂移以及經(jīng)濟(jì)因素的制約,觀測(cè)的地震數(shù)據(jù)往往在空間方向具有不規(guī)則性,影響了后續(xù)自由表面多次波壓制(SRME)、全波形反演以及逆時(shí)偏移成像精度。另外,處理階段對(duì)廢炮、廢道的剔除也是造成數(shù)據(jù)不規(guī)則的重要因素之一。為了提高地震數(shù)據(jù)的橫向連續(xù)性,進(jìn)而提高后續(xù)多道處理以及精細(xì)油氣藏描述的精度,有必要進(jìn)行地震數(shù)據(jù)的插值重建[1-7]。

地震數(shù)據(jù)插值重建方法根據(jù)實(shí)現(xiàn)方式可以分為:基于數(shù)學(xué)變換的插值方法[4-5,8-10]、基于預(yù)測(cè)濾波的插值方法[11]、基于波動(dòng)方程的插值方法[12-13]以及近幾年來(lái)發(fā)展起來(lái)的基于降秩的插值方法[14-17]。其中基于數(shù)學(xué)變換的插值方法相對(duì)簡(jiǎn)單且容易實(shí)現(xiàn),隨著壓縮感知理論的發(fā)展得到了廣泛應(yīng)用。常用的稀疏變換包括Fourier變換、Curvelet變換、Dreamlet變換以及Seislet變換等[18-22]。其中Curvelet變換是一種新興的稀疏變換,其相對(duì)于小波變換增加了角度的信息,更加適用于曲線奇異性的表征,在地震數(shù)據(jù)插值重建以及去噪處理中得到了廣泛的應(yīng)用?；谙∈枳儞Q建立目標(biāo)泛函,借助優(yōu)化方法求解泛函,最終得到插值重建結(jié)果。在眾多優(yōu)化求解方法中,最流行的是凸集投影方法(projection onto convex sets,POCS),因?yàn)樗哂腥菀桌斫狻⒎奖銓?shí)現(xiàn)的優(yōu)點(diǎn)。POCS方法實(shí)際上是一種兩步法:首先利用迭代硬閾值方法得到變換域的解,再將其投影到觀測(cè)平面上[9-10,23]。由于Curvelet變換的冗余度較大(當(dāng)最細(xì)尺度為Curvelet時(shí),2D Curvelet變換的冗余度為8～10,3D Curvelet變換的冗余度為24～32)[22],計(jì)算效率較低,因此目前僅在2D地震資料處理中得到應(yīng)用。張華等[24]及王本鋒等[1]采用2D Curvelet變換對(duì)每一時(shí)間切片或頻率切片進(jìn)行插值,最終實(shí)現(xiàn)了3D地震數(shù)據(jù)的插值重建,但是該方法沒(méi)有利用時(shí)間或頻率的連續(xù)性約束。基于3D Curvelet變換的插值重建方法可以充分利用3D數(shù)據(jù)各個(gè)方向的連續(xù)性約束,但是計(jì)算量較大,必須有效提高計(jì)算效率。基于地震信號(hào)頻譜的共軛對(duì)稱(chēng)性,僅對(duì)有效頻率數(shù)據(jù)體進(jìn)行插值重建,計(jì)算效率可提高一倍以上[5]。實(shí)際上,影響插值重建計(jì)算效率的因素很多,例如稀疏變換種類(lèi)[4,25-26]、插值重建方法[27]、迭代誤差設(shè)計(jì)[23]以及數(shù)據(jù)體的合理表征[5]等。本文主要關(guān)注數(shù)據(jù)體的合理表征,采用規(guī)模較小的數(shù)據(jù)體表征原數(shù)據(jù),通過(guò)對(duì)規(guī)模較小的數(shù)據(jù)體進(jìn)行插值重建,有效提高插值重建的計(jì)算效率。

1 方法原理

1.1 3D Curvelet變換

3D Curvelet變換的數(shù)學(xué)表達(dá)式如下[22]:

(1)

1.2 3D地震數(shù)據(jù)插值重建方法

不規(guī)則數(shù)據(jù)體與規(guī)則完整數(shù)據(jù)體的關(guān)系可由公式(2)近似表征:

(2)

式中:dobs為觀測(cè)的3D不規(guī)則數(shù)據(jù)體,R為設(shè)計(jì)的采樣算子,d0為規(guī)則的完整3D數(shù)據(jù)體。地震數(shù)據(jù)插值重建的目的是從觀測(cè)數(shù)據(jù)以及采樣算子中恢復(fù)規(guī)則完整數(shù)據(jù)體。受觀測(cè)數(shù)據(jù)頻帶有限等因素的影響,方程(2)的求解是不適定的?？紤]到完整數(shù)據(jù)體d0在Curvelet變換域的稀疏性,基于壓縮感知理論,采用稀疏促進(jìn)策略,構(gòu)建目標(biāo)泛函Φ(·)如下:

(3)

式中:C為3D Curvelet變換,CT為3D Curvelet逆變換;λ為正則化因子,用以權(quán)衡擬合殘差以及Curvelet系數(shù)稀疏性。公式(3)可由POCS優(yōu)化方法[4,5,8-10]迭代求解:

(4)

式中:dk是第k步的更新解;I為與采樣算子R同維數(shù)的單位矩陣;Tλk是硬閾值函數(shù)[28-29],λk是閾值,可由改進(jìn)的指數(shù)閾值模型[1]確定。

(5)

式中:λmax,λmin為設(shè)置的最大、最小閾值;N為最大迭代次數(shù);k為當(dāng)前迭代次數(shù)。基于迭代公式(4),迭代次數(shù)達(dá)到設(shè)定最大迭代次數(shù)后,便可以得到插值重建后的地震數(shù)據(jù)。

由于觀測(cè)的地震信號(hào)為時(shí)間信號(hào),對(duì)時(shí)間信號(hào)進(jìn)行一維Fourier變換得到的頻譜具有共軛對(duì)稱(chēng)性質(zhì),因此負(fù)頻率成分可以由正頻率成分進(jìn)行精確估計(jì)。此外,由于觀測(cè)信號(hào)的頻帶有限,利用有效頻帶內(nèi)的頻率成分即可對(duì)原始信號(hào)進(jìn)行高精度的表征。該策略已在2D地震數(shù)據(jù)的插值重建中得到應(yīng)用[5],本文將該策略推廣應(yīng)用到3D地震數(shù)據(jù)的插值重建中。

(6)

(7)

2 數(shù)值算例

2.1 模擬數(shù)據(jù)

用模擬得到的3D數(shù)據(jù)對(duì)上述方法進(jìn)行了測(cè)試,該數(shù)據(jù)共有128炮,每炮128道,每道128個(gè)時(shí)間采樣點(diǎn);炮間距以及道間距均為24m,時(shí)間采樣率為8ms。圖1a為某一固定排列的三維數(shù)據(jù)體的平面展示(圖中①為等時(shí)間切片,②為共炮點(diǎn)道集,③為共檢

波點(diǎn)道集),其隨機(jī)缺失50%的地震道集如圖1b所示,可以看出,地震數(shù)據(jù)的橫向連續(xù)性遭到破壞,這將影響后續(xù)多道處理的精度。

對(duì)不規(guī)則缺失數(shù)據(jù)體進(jìn)行頻譜分析后,選取最大有效頻率為奈奎斯特頻率,即fmax=fN=62.5Hz,利用本文提出的有效頻率域POCS方法以及常用的時(shí)間域POCS方法對(duì)不規(guī)則數(shù)據(jù)體進(jìn)行插值重建。設(shè)最大迭代次數(shù)為50,重建結(jié)果以及重構(gòu)誤差如圖2所示。由圖2a和圖2c可以看出,兩種方法均可以對(duì)缺失地震數(shù)據(jù)進(jìn)行較好的重建,且重構(gòu)誤差(圖2b,圖2d)在允許的誤差范圍內(nèi)。通過(guò)重構(gòu)信噪比進(jìn)一步比較了兩種方法的異同。圖3為兩種方法重構(gòu)信噪比曲線,可見(jiàn)隨著迭代次數(shù)的增加,兩種方法均趨于收斂,收斂時(shí)的信噪比分別為24.0dB和26.3dB,有效頻率域POCS方法精度略低。將重建后的有效頻率數(shù)據(jù)體轉(zhuǎn)換到時(shí)間域,最終的重構(gòu)信噪比為22.1dB。有效頻率域POCS方法的精度相對(duì)時(shí)間域POCS方法的精度稍低的原因之一是時(shí)間采樣率較大,當(dāng)時(shí)間采樣率較小時(shí),有效頻率域POCS方法與時(shí)間域POCS方法的精度相當(dāng),甚至略好[5],具體的影響因素及其量化分析留作后續(xù)研究。在計(jì)算效率上,時(shí)間域POCS方法所用時(shí)間為2918s,有效頻率域POCS方法所用時(shí)間為1512s,計(jì)算效率提高近一倍,驗(yàn)證了理論分析的正確性。

圖1 模擬數(shù)據(jù)a 完整數(shù)據(jù)體; b 不規(guī)則缺失數(shù)據(jù)體(缺失50%)

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

從實(shí)際海洋拖纜數(shù)據(jù)中截取一小塊3D數(shù)據(jù)體對(duì)本文方法的有效性進(jìn)行了測(cè)試。該數(shù)據(jù)體共120炮,每炮120道,每道151個(gè)時(shí)間采樣點(diǎn);炮間距及道間距均為25m,時(shí)間采樣率為8ms。圖4a為該數(shù)據(jù)體平面展示圖(圖中①為等時(shí)間切片,②為共炮點(diǎn)道集,③為共檢波點(diǎn)道集),其隨機(jī)缺失50%的地震道集如圖4b所示,可以看出,缺失后地震數(shù)據(jù)的橫向連續(xù)性變差。

圖2 模擬數(shù)據(jù)測(cè)試結(jié)果a,b 有效頻率域POCS方法重建結(jié)果及重構(gòu)誤差; c,d 時(shí)間域POCS方法重建結(jié)果及重構(gòu)誤差

圖3 模擬算例重構(gòu)信噪比曲線

對(duì)實(shí)際數(shù)據(jù)進(jìn)行頻譜分析,根據(jù)有效頻帶分布范圍,選擇fmax=fN=62.5Hz,利用本文提出的有效頻率域POCS方法及時(shí)間域POCS方法對(duì)圖4b中不規(guī)則數(shù)據(jù)體進(jìn)行插值重建,最大迭代次數(shù)根據(jù)經(jīng)驗(yàn)設(shè)置為50。圖5a為本文方法重建結(jié)果,圖5b為時(shí)間域POCS方法重建結(jié)果,兩種方法重建結(jié)果與完整數(shù)據(jù)體的一致性均較好。圖6展示了兩種方法重構(gòu)信噪比曲線,可見(jiàn)隨著迭代次數(shù)的增加,兩種方法均趨于收斂。將重建后的有效頻率數(shù)據(jù)體轉(zhuǎn)換到時(shí)間域,最終的重構(gòu)信噪比為17.0dB,與時(shí)間域方法的重構(gòu)信噪比17.3dB相近。但是,有效頻率域POCS方法所用時(shí)間為1719s,時(shí)間域POCS方法所用時(shí)間為3396s,即有效頻率域POCS方法的計(jì)算效率相對(duì)時(shí)間域POCS方法提高了近一倍,與理論分析結(jié)果一致,進(jìn)一步驗(yàn)證了本文方法的有效性。當(dāng)減小最大有效頻率時(shí),得到的有效頻率數(shù)據(jù)體規(guī)模將減小,可以提高插值重建的計(jì)算效率,但是以犧牲插值重建精度為代價(jià),具體的量化分析將留作后續(xù)研究。

圖4 實(shí)際資料a 完整數(shù)據(jù)體; b 不規(guī)則缺失數(shù)據(jù)體(缺失50%)

圖5 實(shí)際資料重建結(jié)果a 有效頻率域POCS方法; b 時(shí)間域POCS方法

圖6 實(shí)際資料重構(gòu)信噪比曲線

3 結(jié)論與認(rèn)識(shí)

本文基于3D Curvelet稀疏變換,利用頻率域地震數(shù)據(jù)的共軛對(duì)稱(chēng)性質(zhì),得到規(guī)模減半的有效頻率數(shù)據(jù)體,應(yīng)用凸集投影(POCS)方法,研究了高效的3D地震數(shù)據(jù)插值重建方法。研究結(jié)果表明:有效頻率域POCS方法是一種有效的插值重建手段;對(duì)規(guī)模減半的有效頻率數(shù)據(jù)體進(jìn)行插值重建,能在保證插值精度的同時(shí),提高約一倍的計(jì)算效率。

本文僅僅圍繞地震數(shù)據(jù)插值重建問(wèn)題展開(kāi)研究,關(guān)于插值去噪一體化以及考慮衰減信息的插值重建,將是我們下一步的研究計(jì)劃。由于Curvelet變換的計(jì)算效率較低,因此高效的稀疏變換以及稀疏字典的研究是我們另一個(gè)研究方向。

[1] 王本鋒,陳小宏,李景葉,等.POCS聯(lián)合改進(jìn)的Jitter采樣理論曲波域地震數(shù)據(jù)重建[J].石油地球物理勘探,2015,50(1):20-28

WANG B F,CHEN X H,LI J Y,et al.Seismic data reconstruction based on POCS and improved Jittered sampling in the Curvelet domain[J].Oil Geophysical Prospecting,2015,50(1):20-28

[2] 王本鋒.基于反演理論的地震信號(hào)處理方法研究[D].北京:中國(guó)石油大學(xué)(北京),2015

WANG B F.Research on seismic data processing technology by inverse theory[D].Beijing:China University of Petroleum,2015

[3] 曹靜杰,王本鋒.基于一種改進(jìn)凸集投影方法的地震數(shù)據(jù)同時(shí)插值和去噪[J].地球物理學(xué)報(bào),2015,58(8):2935-2947

CAO J J,WANG B F.An improved projection onto convex sets method for simultaneous interpolation and denoising[J].Chinese Journal of Geophysics,2015,58(8):2935-2947

[4] WANG B F,WU R S,GENG Y,et al.Dreamlet-based interpolation using POCS method[J].Journal of Applied Geophysics,2014,109(10):256-265

[5] WANG B F.An efficient POCS interpolation method in the frequency-space domain[J].IEEE Geoscience and Remote Sensing Letters,2016,13(9):1384-1387

[6] HERRMANN F J,HENNENFENT G.Non-parametric seismic data recovery with Curvelet frames[J].Geophysical Journal International,2008,173(1):233-248

[7] HENNENFENT G,HERRMANN F J.Simply denoise:wavefield reconstruction via Jittered undersampling[J].Geophysics,2008,73(3):V19-V28

[8] WANG B F,CHEN X H,LI J Y,et al.An improved weighted projection onto convex sets method for seismic data interpolation and denoising[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2016,9(1):228-235

[9] WANG B F,LI J Y,CHEN X H.A novel method for simultaneous seismic data interpolation and noise removal based on the L0norm constraint[J].Journal of Seismic Exploration,2015,24(2):187-204

[10] WANG B F,WU R S,CHEN X H,et al.Simultaneous seismic data interpolation and denoising with a new adaptive method based on dreamlet transform[J].Geophysical Journal International,2015,201(2):1180-1192

[11] SPITZ S.Seismic trace interpolation in the F-X domain[J].Geophysics,1991,56(6):785-794

[12] RONEN J.Wave-equation trace interpolation[J].Geophysics,1987,52(7):973-984

[13] FOMEL S.Seismic reflection data interpolation with differential offset and shot continuation[J].Geophysics,2003,68(2):733-744

[14] GAO J J,Sacchi D M,Chen X H.A fast reduced-rank interpolation method for prestack seismic volumes that depend on four spatial dimensions[J].Geophysics,2013,78(1):V21-V30

[16] MA J W.Three-dimensional irregular seismic data reconstruction via low-rank matrix completion[J].Geophysics,2013,78(5):V181-V192

[17] 馬繼濤,王建花,劉國(guó)昌.基于頻率域奇異值分解的地震數(shù)據(jù)插值去噪方法研究[J].石油物探,2016,55(2):205-213

MA J T,WANG J H,LIU G C.Seismic data noise attenuation and interpolation using singular value decomposition in frequency domain[J].Geophysical Prospecting for Petroleum,2016,55(2):205-213

[18] 高建軍,陳小宏,李景葉,等.基于 POCS 方法指數(shù)閾值模型的不規(guī)則地震數(shù)據(jù)重建[J].應(yīng)用地球物理,2010,7(3):229-238

GAO J J,CHEN X H,LI J Y,et al.Irregular seismic data reconstruction based on exponential threshold model of POCS method[J].Applied Geophysics,2010,7(3):229-238

[19] GAN S W,WANG S D,CHEN Y K,et al.Compressive sensing for seismic data reconstruction via fast projection onto convex sets based on seislet transform[J].Journal of Applied Geophysics,2016,130(7),194-208

[20] GAN S W,WANG S D,CHEN Y K,et al.Dealiased seismic data interpolation using seislet transform with low-frequency constraint[J].IEEE Geoscience and Remote Sensing Letters,2015,12(10):2150-2154

[21] LIU W,CAO S Y,GAN S W,et al.One-step slope estimation for dealiased seismic data reconstruction via iterative seislet thresholding[J].IEEE Geoscience and Remote Sensing Letters,2016,13(10):1462-1466

[23] WANG B F,CHEN X H,WANG T.A novel reconstruction error definition for seismic data interpolation[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015,3875-3879

[24] 張華,陳小宏.基于 Jitter 采樣和曲波變換的三維地震數(shù)據(jù)重建[J].地球物理學(xué)報(bào),2013,56(5):1637-1649

ZHANG H,CHEN X H.Seismic data reconstruction based on Jittered sampling and Curvelet transform[J].Chinese Journal of Geophysics,2013,56(5):1637-1649

[25] CAO J J,ZHAO J T.Simultaneous seismic interpolation and denoising based on sparse inversion with a 3D low redundancy Curvelet transform[J].Exploration Geophysics,2017,48(4):422-429

[26] CAO J J,ZHAO J T.3D seismic interpolation with a low redundancy,fast Curvelet transform[J].Journal of Seismic Exploration,2015,24(2):121-134

[27] ZU S H,ZHOU H,CHEN Y K,et al.Interpolating big gaps using inversion with slope constraint[J].IEEE Geoscience and Remote Sensing Letters,2016,13(9):1369-1373

[28] BLUMENSATH T,DAVIES M E.Iterative hard thresholding for compressed sensing[J].Applied and Computational Harmonic Analysis,2009,27(3):265-274

[29] BLUMENSATH T,DAVIES M E.Iterative thresholding for sparse approximations[J].Journal of Fourier Analysis and Applications,2008,14(5):629-654