盧昕婷, 韓立國, 張盼, 孫宏宇

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

?

基于全變分原理的多震源混合數(shù)據(jù)直接偏移方法

盧昕婷, 韓立國*, 張盼, 孫宏宇

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院， 長春 130026

多震源混合地震采集技術(shù)，即將多個震源以一定編碼方式連續(xù)地激發(fā)，得到多炮混合的地震數(shù)據(jù).該技術(shù)能減少地震采集時間，節(jié)約采集成本，但是混合數(shù)據(jù)的直接偏移會在成像剖面中引入嚴重的串?dāng)_噪聲，影響成像效果.從數(shù)學(xué)上看，地震成像屬于典型的數(shù)學(xué)物理反問題，可以采用線性反演方法求解一個正則化約束的最小二乘(LS)優(yōu)化問題，獲得更高質(zhì)量的成像結(jié)果.全變分(TV)正則化方法是圖像去噪和復(fù)原領(lǐng)域中廣泛應(yīng)用的熱點技術(shù),其能在去除噪聲的過程中保留圖像的邊緣信息和不連續(xù)性.在對TV圖像去噪復(fù)原方法原理分析的基礎(chǔ)上，本文將多震源混合數(shù)據(jù)直接偏移成像問題轉(zhuǎn)換成圖像復(fù)原的極小化能量泛函問題，用TV正則化代替?zhèn)鹘y(tǒng)最小二乘偏移(LSM)中的L2范數(shù)正則化，提出基于全變分原理的混合數(shù)據(jù)直接偏移方法.該方法使用基于梯度的快速迭代收縮閾值與快速梯度投影組合算法——FISTA/FGP求解最優(yōu)化問題，能有效壓制串?dāng)_噪聲，增強同相軸連續(xù)性，提高成像分辨率.理論模型測試結(jié)果表明：將本方法應(yīng)用于混合數(shù)據(jù)，無論是去噪效果還是成像精度都得到顯著改善.

多震源混合采集； 混合數(shù)據(jù)； 直接成像； 串?dāng)_噪聲； 全變分

1 引言

多震源混合采集(blended acquisition)是近年來興起的一種新型地震采集技術(shù)，該技術(shù)由Berkhout(2008)和Berkhout等(2009)在多震源同時采集技術(shù) (simultaneous acquisition)(Silverman,1979; Beasley et al., 1998; Bagaini, 2006; Beasley, 2008)的基礎(chǔ)上發(fā)展而來.在混合采集系統(tǒng)中，不同空間位置的多個震源以一定的編碼方式連續(xù)激發(fā)，得到的地震記錄為波場混疊的多炮混合數(shù)據(jù).由于震源之間的激發(fā)時間間隔很小，這種采集方式極大地提高了數(shù)據(jù)的獲取量和采集效率.

目前，針對混合數(shù)據(jù)的偏移成像方法可分為兩類：一是直接對混合數(shù)據(jù)進行偏移成像，不做任何預(yù)處理.該方法具有較高的處理效率，但是由于波場數(shù)據(jù)相互混疊干擾，直接偏移的成像剖面上會出現(xiàn)串?dāng)_噪聲(crosstalk)，嚴重降低成像質(zhì)量；二是先將多炮混合記錄進行炮分離處理(deblending)，對得到的單炮數(shù)據(jù)進行常規(guī)偏移成像.該方法處理效率較低，炮分離的效果將影響偏移成像結(jié)果.綜合比較，直接偏移方法將是未來混合數(shù)據(jù)成像的主要發(fā)展方向，其研究重點在于如何有效地壓制串?dāng)_噪聲.

Tang和Biondi(2009)在其論文中首次提出使用最小二乘偏移(LSM)來壓制混合數(shù)據(jù)直接成像的串?dāng)_噪聲.LSM可以使用多種偏移算子，具有較高靈活性，同時能很好地適應(yīng)混合采集技術(shù)的處理要求.以上性質(zhì)使得該方法成為現(xiàn)階段混合數(shù)據(jù)直接偏移的主流方法，近年來也陸續(xù)涌現(xiàn)出成功的應(yīng)用實例(Schuster and Dai, 2010; Verschuur and Berkhout,2011; Dai et al., 2011, 2012; Huang and Schuster, 2012).

LSM屬于數(shù)學(xué)上的病態(tài)問題，需要引入正則化約束來確保求解過程的穩(wěn)定.傳統(tǒng)的Tikhonov正則化(Tikhonov and Arsenin, 1977)采用光滑性約束如L2范數(shù)約束，具有良好的算法穩(wěn)定性，但它對目標(biāo)模型有平滑濾波效應(yīng)，容易導(dǎo)致模型過度光滑，邊緣模糊，不連續(xù)信息不能被精確重構(gòu).為了減弱這種平滑效應(yīng)，提高成像質(zhì)量，本文提出基于全變分(TV)(total variation)原理的TV范數(shù)約束偏移算法.

全變分(Rudin et al., 1992)是圖像去噪和復(fù)原領(lǐng)域的研究熱點技術(shù)之一，其基本原理是將圖像去噪/復(fù)原問題轉(zhuǎn)換成圖像模型的能量泛函最小化問題.全變分原理在地震圖像去噪(屈勇等，2011；公成敏，2014)、地震反演(Askan et al., 2007; Burstedde and Ghattas, 2009; 毛衡等, 2012)以及地震偏移速度模型重建(Anagaw, 2009; Anagaw and Sacchi, 2012)方面已經(jīng)得到應(yīng)用.地下介質(zhì)大多是不規(guī)則的，并且具有明顯的成層性，全變分模型采用非光滑性約束，在去噪時能夠不平滑目標(biāo)邊界保留圖像的不連續(xù)性特征.因此，為了在壓制串?dāng)_噪聲的同時最大限度地還原地下介質(zhì)的真實信息，本文將全變分應(yīng)用到多震源地震數(shù)據(jù)直接偏移成像中，建立全變分地震偏移成像模型，提出基于TV范數(shù)約束的混合數(shù)據(jù)直接偏移方法.考慮計算成本，本方法采用快速迭代收縮閾值與快速梯度投影組合算法——FISTA/FGP(Fast Iterative Shrinkage-Thresholding Algorithm/Fast Gradient Projection)(Beck and Teboulle, 2009)求解目標(biāo)函數(shù)的最優(yōu)化問題.理論模型的測試結(jié)果表明：TV范數(shù)約束的混合數(shù)據(jù)直接偏移能有效壓制噪聲，保護地震剖面的邊緣信息，增強剖面的結(jié)構(gòu)特征，提高成像精度.

2 混合采集與直接偏移成像

在波動方程波恩近似的條件下，地震數(shù)據(jù)d和地下反射系數(shù)m有如下線性關(guān)系：

d=Lm，

(1)

其中，L為正演模擬算子.

混合采集數(shù)據(jù)是來自多個震源的炮集疊加記錄，即超級炮，這種混合過程可表示如下：

(2)

將式(1)代入式(2)中有：

(3)

(4)

同常規(guī)單炮偏移算子一致，超級炮偏移算子也可以用正演算子的伴隨矩陣表示，即

(5)

(6)

將式(3)代入式(6)有：

(7)

由(7)式可以看出偏移成像相當(dāng)于地下真實反射系數(shù)的模糊表示.

3 全變分多震源地震數(shù)據(jù)偏移成像

3.1 全變分偏移成像數(shù)學(xué)模型

二維地震圖像的TV范數(shù)有兩種表示方法:

(1)基于L2范數(shù)的各向同性擴散模型(Chambolle, 2004):

(8)

式中，n1和n2為m的維數(shù).

(2)基于L1范數(shù)的各向異性擴散模型，即經(jīng)典的ROF模型:

(9)

各向同性擴散模型在重構(gòu)時會增加多余的光滑作用，模糊圖像邊緣；各向異性擴散模型具有良好的邊緣保持特性，但是該模型可能會將噪聲處理成假邊緣，在平滑區(qū)出現(xiàn)階梯效應(yīng).由于地震圖像具有不連續(xù)性，通常表現(xiàn)為分段光滑的曲線或者直線奇異性，各向異性擴散模型在不同圖像特征區(qū)域內(nèi)擴散特性不一樣，有利于保持地震圖像的不連續(xù)性特征.特別說明一點，各向異性擴散模型在某種程度上可以看作特殊的L1范數(shù)稀疏約束，其本質(zhì)就是圖像梯度的L1范數(shù)，對于稀疏度比較高的模型，它能更好地保留解的內(nèi)在稀疏性.

本文將式(9)定義的TV范數(shù)約束應(yīng)用到混合數(shù)據(jù)直接偏移成像中以達到減弱串?dāng)_噪聲，提高成像精度的目的.

偏移成像的TV范數(shù)約束最優(yōu)化問題如下所示：

min‖m‖TV, subject toLm=d.

(10)

根據(jù)正則化理論，將上述優(yōu)化問題轉(zhuǎn)化為最小化成本函數(shù)J(m):

(11)

式(11)右邊由兩項組成：第一項為數(shù)據(jù)項，又叫保真項，反映了近似解與測量數(shù)據(jù)之間的誤差，其作用是在去噪過程中保留原始圖像特性，確保解的真實性；第二項為平滑項，即正則項，λ為正則化參數(shù)，用來調(diào)節(jié)數(shù)據(jù)項與平滑項的相對貢獻量.在全變分模型中，λ的取值對計算結(jié)果起著十分重要的影響：λ取值越大，TV范數(shù)占的權(quán)重越大，圖像去噪效果更好，但解可能會和測量數(shù)據(jù)有較大差異；若λ取值過小，會導(dǎo)致圖像邊緣比較模糊，得不到理想的去噪效果.在實際應(yīng)用中，λ的選取和圖像噪聲水平有關(guān).

3.2 全變分快速組合算法

本文采用Beck和Teboulle提出的快速組合算法FISTA/FGP求解式(11).該算法結(jié)合FISTA和FGP二者優(yōu)點，實現(xiàn)簡單，穩(wěn)定性好，有較高的計算效率.具體操作如下.

首先考慮L=I時的特殊情況，即

(12)

在FISTA/FGP中，F(xiàn)GP算法求解TV約束最優(yōu)化問題的基本思路是先構(gòu)造式(12)的對偶問題，然后基于梯度類方法求解.該對偶問題解的形式如式(13)所示：

m=PBl,u(d-λY(p,q)).

(13)

所需參數(shù)定義如下：

(1)Φ為矩陣對(p,q)的集合，其中p∈R(n1-1)×n2,q∈Rn1×(n2-1)且滿足：

(14)

(2)Y為線性運算符,其將矩陣對(p,q)映射為一個n1×n2的矩陣：

(Y(p,q))i,j=pi,j-pi-1,j+qi,j-qi,j-1,

i=1,…,n1-1,j=1,…,n2-1，

(15)

其中p0,j≡pn1,j≡qi,0≡qi,n2≡0.

(3)PBl,u為正交投影算子：

(16)

FGP算法實現(xiàn)流程如下：

初始賦值：

(r1,s1)=(p0,q0)=(O(n1-1)×n2,On1×(n2-1))，t1=1. 第k步迭代：

(pk,qk)=

(17a)

(17c)

(17d)

經(jīng)典的梯度方法中，關(guān)于式(11)的第k階迭代優(yōu)化公式為:

(18)

式中，α為迭代步長，控制收斂速度.FISTA/FGP中，確保迭代收斂需要保證α≥max eig(LHL)，即α的值必須大于LHL的最大特征值.

如果定義：

(19)

(20)

對比式(20)與式(12)，我們可以使用FISTA和FGP組合算法求解式(20).

FISTA的計算過程如下所示：

初始賦值：y1=m0，t1=1.

第k步迭代：

mk=DL(yk),

(21a)

(21b)

(21c)

其中，DL是算法的核心，以FGP替換DL，就成為FISTA/FGP算法.

4 理論模型測試

本文通過層狀介質(zhì)模型，SEG/EAGE斷層模型與Marmousi模型測試方法的準(zhǔn)確性和適應(yīng)性.其中，層狀模型主要測試方法的有效性，SEG/EAGE斷層模型主要測試方法的抗噪性和正則化參數(shù)選擇對偏移結(jié)果的影響，最后通過Marmousi模型測試方法在復(fù)雜介質(zhì)條件下的適應(yīng)性.

4.1 層狀介質(zhì)模型偏移

層狀介質(zhì)模型如圖1所示.模型縱橫向網(wǎng)格點數(shù)均為100，縱橫向網(wǎng)格間距均為10 m，背景速度為2500 m·s-1.檢波器和單震源個數(shù)均為100個，時間采樣間隔為1 ms，震源子波為主頻30 Hz的雷克子波.通過正演模擬得到100個單炮記錄，使用隨機時間間隔編碼算子，每5個單炮合成一個混合炮(圖2).數(shù)值試驗中所有偏移均采用克?；舴蚱扑阕樱鋵τ^測系統(tǒng)適應(yīng)性強，對速度模型敏感度比較低.

圖1 層狀模型(反射系數(shù))Fig.1 Layered model (reflectivity)

圖2 層狀模型混合炮記錄Fig.2 A blended shot gather of layered model

圖3a為混合數(shù)據(jù)常規(guī)克?；舴蚱瞥上衿拭?，圖3b為混合數(shù)據(jù)L2范數(shù)約束LSM成像剖面，圖3c為混合數(shù)據(jù)TV范數(shù)約束偏移成像剖面.從圖3可以看出，和常規(guī)偏移成像相比，LSM能有效去除一部分串?dāng)_噪聲(圖3矩形框內(nèi)所示)和偏移畫弧(圖3箭頭所示)，但在最終成像剖面上還有明顯的噪聲殘留，反射能量聚集也不夠集中.而TV范數(shù)約束偏移對串?dāng)_噪聲和偏移畫弧的壓制效果非常顯著，成像剖面變得干凈，并且偏移能量聚焦更佳，同相軸變清晰，更加接近真實反射系數(shù).由此驗證了TV范數(shù)約束在分層模型成像上的優(yōu)勢.

4.2 斷層模型偏移

圖3 混合數(shù)據(jù)(a)克希霍夫偏移結(jié)果，(b) LSM結(jié)果，(c) TV范數(shù)約束偏移結(jié)果Fig.3 Imaging results from (a) Kirchhoff migration, (b) LSM, (c) TV norm constrained migration for blended data

圖4 SEG/EAGE斷層模型(反射系數(shù))Fig.4 SEG/EAGE fault model (reflectivity)

由于實際地震數(shù)據(jù)中通常包含有一定程度的噪聲，本文接下來使用SEG/EAGE斷層模型的含噪數(shù)據(jù)進行數(shù)值試驗.為了減輕計算負擔(dān)，在不降低計算精確度和模型復(fù)雜性的前提下，將原有速度模型抽稀(圖4).重采樣之后模型的橫向和縱向網(wǎng)格點數(shù)均為100，網(wǎng)格間距均為10 m，背景速度為2500 m·s-1.檢波器和單震源個數(shù)均為100個，從模型網(wǎng)格的起始點依次排列.時間采樣間隔為2 ms，震源子波為主頻30 Hz的雷克子波，每5個單炮合成一個混合炮(圖5a).單炮混合之后，在數(shù)據(jù)中加入15%的高斯白噪聲(圖5b)，使用不同取值的λ進行混合數(shù)據(jù)TV范數(shù)約束偏移(圖6).

圖6a為混合數(shù)據(jù)克希霍夫偏移成像結(jié)果，圖6b，6c，6d分別為λ=0.00001，0.0001，0.001時的TV范數(shù)約束偏移成像結(jié)果.如圖所示，混合數(shù)據(jù)常規(guī)偏移成像剖面上，噪聲影響嚴重，成像質(zhì)量差.當(dāng)λ=0.0001時，TV范數(shù)約束偏移能較為理想的完成對串?dāng)_噪聲、偏移噪聲和隨機噪聲的衰減處理，同時成像結(jié)果貼近真實模型.當(dāng)λ=0.00001時，方法去噪力度不夠.隨著λ取值增大，去噪能力增強，但是當(dāng)λ=0.001時，深層的一部分真實信息也被去除.因此參數(shù)λ的取值，需要綜合考慮偏移成像的去噪效果和數(shù)據(jù)保真度.

4.3 復(fù)雜模型偏移

最后使用垂直射線和常密度假設(shè)下的局部Marmousi模型(圖7)考察本方法對復(fù)雜地質(zhì)構(gòu)造的成像效果.模型橫向和縱向網(wǎng)格點數(shù)分別為551和221，縱橫向網(wǎng)格間距均為10 m，背景速度為3000 m·s-1.檢波器和單震源個數(shù)均為276個，炮間距和道間距均為20 m.時間采樣間隔為1 ms，震源子波為主頻30 Hz的雷克子波.通過正演模擬得到276個單炮記錄，每4個單炮合成一個混合炮.

圖5 (a)混合炮記錄；(b)含信噪比15%高斯白噪的混合炮記錄Fig.5 (a) is the blended shot gather. (b) is the blended shot gather that is added with white Gaussian noise (SNR=15%).

圖6 混合數(shù)據(jù)(a)克?；舴蚱平Y(jié)果，(b) λ=0.00001時TV范數(shù)約束偏移結(jié)果，(c) λ=0.0001時TV范數(shù)約束偏移結(jié)果，(d) λ=0.001時TV范數(shù)約束偏移結(jié)果Fig.6 Imaging results for blended data from (a) Kirchhoff migration, (b) TV norm constrained migration (λ=0.00001), (c) TV norm constrained migration (λ=0.0001), (d) TV norm constrained migration (λ=0.001)

圖7 Marmousi模型(反射系數(shù))Fig.7 Marmousi model (reflectivity)

混合數(shù)據(jù)常規(guī)偏移成像如圖8a所示，圖8b和圖8c分別為常規(guī)單炮數(shù)據(jù)和混合數(shù)據(jù)TV范數(shù)約束偏移成像.由圖可以看出，混合數(shù)據(jù)常規(guī)偏移成像效果差，剖面上存在嚴重的串?dāng)_噪聲和偏移噪聲污染，成像精度低.TV范數(shù)約束偏移通過多次迭代對繞射點能量反復(fù)聚焦，使繞射波的歸位更準(zhǔn)確，串?dāng)_噪聲和偏移噪聲得到有效壓制.同時，整個成像剖面上粗糙的同相軸變得清晰，連續(xù)性增強，背斜、斷層、薄互層等復(fù)雜地質(zhì)構(gòu)造都達到了很好的成像效果與精度，成像分辨率顯著提高.

取圖7和圖8中方框1和方框2內(nèi)的圖像放大顯示(圖9和圖10)，能更明顯地看出TV范數(shù)約束偏移的成像效果.圖11為Marmousi模型偏移成像和反射系數(shù)任意兩道數(shù)據(jù)的波形對比圖.相比較之下，TV范數(shù)約束偏移的成像結(jié)果和真實反射系數(shù)有較高相似度，波峰波谷吻合度較好.

5 計算速率分析

對于多震源混合數(shù)據(jù)直接偏移成像，TV范數(shù)約束偏移每次迭代都需要進行1次正演和1次偏移.因此N次迭代，其計算成本大約是常規(guī)偏移的2N倍，計算量十分巨大.

圖10 方框2內(nèi)數(shù)據(jù)對比(a) 混合數(shù)據(jù)克?；舴蚱?； (b) 反射系數(shù)； (c) 常規(guī)單炮數(shù)據(jù)TV范數(shù)約束偏移； (d) 混合數(shù)據(jù)TV范數(shù)約束偏移.Fig.10 The data in rectangle 2(a) Kirchhoff migration for blended data; (b) Reflectivity; (c) TV norm constrained migration for conventional sources data; (d) TV norm constrained migration for blended data.

圖11 偏移數(shù)據(jù)與反射系數(shù)波形對比圖1號線為混合數(shù)據(jù)克希霍夫偏移，2號線為混合數(shù)據(jù)TV范數(shù)約束偏移，3號線為常規(guī)單炮數(shù)據(jù)TV范數(shù)約束偏移，4號線為真實反射系數(shù).(a)第300道的數(shù)據(jù)；(b)第450道的數(shù)據(jù).Fig.11 The waveform comparison of migration data and reflectivityThe first line is Kirchhoff migration for blended data, the second line is TV norm constrained migration for blended data, the third line is TV norm constrained migration for conventional sources data, the fourth line is reflectivity. (a) is the 300th traces, (b) is the 450th traces.

此類迭代偏移方法以計算成本換取成像精度，實際使用時，最終迭代次數(shù)需要權(quán)衡這兩方面的影響，兼顧方法的有效性和實用性.隨著計算機計算性能提高，GPU等快速計算設(shè)備的應(yīng)用以及并行編程技術(shù)的發(fā)展，有望解決計算效率難題，使本方法對二維乃至三維大規(guī)模多震源混合數(shù)據(jù)的處理成為可能.

6 結(jié)論

多震源混合地震采集技術(shù)是對傳統(tǒng)單炮地震采集技術(shù)的重大革新，而混合數(shù)據(jù)直接成像處理將是支撐該技術(shù)的發(fā)展方向之一.針對混合數(shù)據(jù)直接偏移所引起的串?dāng)_噪聲，本文提出基于全變分原理的混合數(shù)據(jù)TV范數(shù)約束偏移.使用TV范數(shù)作為正則化約束項可以測度圖像的跳躍性,不需要強制要求解的平滑性，圖像的邊緣紋理在成像過程中能夠得到保存.

本方法屬于混合數(shù)據(jù)直接偏移，無需炮分離處理，提高了數(shù)據(jù)處理效率.利用FISTA/FGP快速算法求解優(yōu)化問題，實現(xiàn)方便，有較高收斂速度.層狀介質(zhì)模型，SEG/EAGE斷層模型和Marmousi模型的TV范數(shù)約束偏移試算都獲得理想的成像結(jié)果：在最終成像剖面上，隨機噪聲、偏移弧狀噪聲和串?dāng)_噪聲都得到有效壓制；偏移能量正確歸位；各尺度地質(zhì)構(gòu)造都有較好的成像精度，成像分辨率明顯提高.TV范數(shù)約束偏移適用于常見的偏移算法(波動方程偏移，逆時偏移等)，計算效率和結(jié)果精度因所選偏移算子而有一定的差異.

最后需要注意，TV范數(shù)約束偏移是基于最小二乘優(yōu)化的迭代偏移算法，通過增加偏移迭代次數(shù)來達到去噪和提高成像分辨率的目的.本算法能改善傳統(tǒng)偏移成像方法的成像效果，但是計算耗時也相應(yīng)增加，因此迭代次數(shù)的選擇要綜合考慮計算成本和成像精度.隨著計算機運行速度的提升以及并行編程技術(shù)的發(fā)展，計算時間存在大量減少的可能，為本方法的實際應(yīng)用提供了有希望的前景.

Anagaw A Y. 2009. Total variation and adjoint state methods for seismic wavefield imaging [Master′s thesis]. Alberta: Physics Department of University of Alberta.

Anagaw A Y, Sacchi M D. 2012. Edge-preserving seismic imaging using the total variation method.JournalofGeophysicsandEngineering, 9(2): 138-146.

Askan A, Akcelik V, Bielak J, et al. 2007. Full waveform inversion for seismic velocity and anelastic losses in heterogeneous structures.BulletinoftheSeismologicalSocietyofAmerica, 97(6): 1990-2008.

Bagaini C. 2006. Overview of simultaneous vibroseis acquisition methods. ∥ 76th Annual International Meeting, SEG, Expanded Abstracts, 70-74. Beasley C J, Chambers R E, Jiang Z R. 1998. A new look at simultaneous sources. ∥ SEG Technical Program Expanded Abstracts, 133-135.

Beasley C J. 2008. Simultaneous sources: A technology whose time has come. ∥ 78th Annual International Meeting, SEG, Expanded Abstracts, 2796-2800.Beck A, Teboulle M. 2009. Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems.IEEETransactionsonImageProcessing, 18(11): 2419-2434.

Berkhout A J. 2008. Changing the mindset in seismic data acquisition.TheLeadingEdge, 27(7): 924-938.

Berkhout A J, Blacquière G, Verschuur D J. 2009. The concept of double blending: Combining incoherent shooting with incoherent sensing.Geophysics, 74(4): A59-A62.Burstedde C, Ghattas O. 2009. Algorithmic strategies for full waveform inversion: 1D experiments.Geophysics, 74(6): WCC37-WCC46.

Chambolle A. 2004. An algorithm for total variation minimization and applications.JournalofMathematicalImagingandVision, 20(1): 89-97.

Dai W, Wang X, Schuster G T. 2011. Least-squares migration of multisource data with a deblurring filter.Geophysics, 76(5): R135-R146.

Dai W, Fowler P, Schuster G T. 2012. Multi-source least-squares reverse time migration.GeophysicalProspecting, 60(4): 681-695. Gong C M. 2014. Application of variational principle in seismic data denoising.Computer&DigitalEngineering(in Chinese), 42(7): 1271-1274.

Huang Y S, Schuster G T. 2012. Multisource least-squares migration of marine streamer and land data with frequency-division encoding.GeophysicalProspecting, 60(4): 663-680.

Mao H, Wang W, Han B. 2012. Seismic waveform inversion based on Bayesian-sparsity constraint regularization method.CommunicationonAppliedMathematicsandComputation(in Chinese), 26(3): 285-297.

Qu Y, Cao J X, Zhu H D, et al. 2011. An improved total variation technique for seismic image denoising.ActaPetroleiSinica(in Chinese), 32(5): 815-819.

Rudin L I, Osher S, Fatemi E. 1992. Nonlinear total variation based noise removal algorithms.PhysicaD:NonlinearPhenomena, 60(1-4): 259-268.Schuster G T, Dai W. 2010. Multi-source wave equation least-squares migration with a deblurring filter. ∥ 72nd EAGE Conference & Exhibition, Extended Abstracts, 276-279. Silverman D. 1979. Method of three dimensional seismic prospecting. U. S. Patent 4, 159, 463, 1979.Tang Y X, Biondi B. 2009. Least-squares migration/inversion of blended data. ∥ SEG Technical Program Expanded Abstracts, 28: 2859-2863.

Tikhonov A N, Arsenin V Y. 1977. Solutions of Ill-Posed Problems. Washington, DC: V. H. Winston and Sons.

Verschuur D J, Berkhout A J. 2011. Seismic migration of blended shot records with surface-related multiple scattering.Geophysics, 76(1): A7-A13.

附中文參考文獻

公成敏. 2014. 全變分原理在地震數(shù)據(jù)去噪中的應(yīng)用. 計算機與數(shù)字工程, 42(7): 1271-1274.

毛衡, 王薇, 韓波. 2012. 基于貝葉斯-稀疏約束正則化方法的地震波形反演. 應(yīng)用數(shù)學(xué)與計算數(shù)學(xué)學(xué)報, 26(3): 285-297.

屈勇, 曹俊興, 朱海東等. 2011. 一種改進的全變分地震圖像去噪技術(shù). 石油學(xué)報, 32(5): 815-819.

(本文編輯 何燕)

Direct migration method of multi-source blended data based on total variation

LU Xin-Ting, HAN Li-Guo*, ZHANG Pan, SUN Hong-Yu

CollegeofGeo-explorationSciencesandTechnology,JilinUniversity，Changchun130026,China

In acquisition of multi-source blended data, multiple sources are continuously shot and geophones continuously receive seismic signals to acquire blended shot records overlapping in both the spatial and temporal domains. Relative to the traditional seismic acquisition, blended acquisition can decrease the effective survey duration and reduce the cost of surveys. However, this continuously shooting and recording strategy can lead to complex blended wave fields and create great difficulties to the following migration imaging. There are two migration imaging methods about the blended data at this stage: to perform the migration processing directly on the blended data without any pre-separation, or recovery the blended data to individual shot data, which is called “deblending”, then conduct standard migration processing to these deblended data. The first method has a high processing efficiency, but the multi-source direct migration is usually not satisfactory because of the crosstalk contamination. The second method has low efficiency and the source separation has serious influence on imaging quality. All things considered, the direct migration method has a better prospect of application, so this article studies a direct migration method of blended data which can effectively remove the crosstalk.Mathematically, the seismic imaging can be regarded as a typical mathematical and physical inverse problem and a better imaging result can be obtained by solving a regularized least-squares (LS) problem with a linear inversion method. The total variation (TV) regularization method is widely used in image denoising and restoration fields. It can produce a denoised image where edges and discontinuities are preserved. On the basis of analyzing the principles of TV image denoising and restoration method, in this paper, we formulate the direct migration of multi-source seismic data as a problem of minimizing the energy functional of image restoration, use the TV regularization to replace the L2 norm regularization in the traditional least-squares migration (LSM) and propose a TV norm constrained migration method. The proposed method employs a gradient-based algorithm—a combination of fast iterative shrinkage-thresholding algorithm and fast gradient projection method (FISTA/FGP) to solve the optimization problem. It combines the advantages of FISTA and FGP, which is simple in algorithm, convenient to realize and stable.To test the validity and adaptability of our migration method, numerical experimentations are carried out on a horizontally-layered medium model, SEG/EAGE fault model and 2D Marmousi model, respectively. First, we compute the prestack Kirchhoff migration, the L2 norm constrained LSM and TV norm constrained migration for the horizontally-layered medium model. The results show that the TV norm constrained migration can effectively suppress the migration artifacts and crosstalk and the image quality is better than the other two. For the case of data containing noise, this article adds 15% Gauss white noise into the synthetic blended data of the SEG/EAGE fault model and applies the TV norm constrained migration with different values of the optimization parameters. The results indicate that our algorithm has good anti-noise performance and point out that there is influence of the selection of regularization parameters on the imaging results. So considering the denoise effect and fidelity of migration imaging, an assigned range of these optimization parameters is proposed. For the complex medium case, we compute the standard migration and TV norm constrained migration for the blended data of the 2D Marmousi model and use the conventional data TV norm constrained migration as a reference. The imaging results and waveform comparison both indicate that our method can yield desired migration imaging of complex geological structures.For the crosstalk noise caused by direct imaging of blended data, we propose a direct migration method based on TV regularization. The results of theoretical model experiments show that our method can effectively suppress the random noise, migration artifacts and crosstalk noise, enhance the continuity of seismic events and improve the imaging resolution. TV norm constrained migration belongs to an iterative migration algorithm based on least-square optimization and has a prominent imaging effect. However, it is also a computing-expensive algorithm. With the development of the computer performance and the parallel programming technology, the computational cost can be reduced, which can enable our method to practical applications.

Multi-source blended acquisition; Blended data; Direct imaging; Crosstalk noise; Total variation

盧昕婷, 韓立國, 張盼等. 2015. 基于全變分原理的多震源混合數(shù)據(jù)直接偏移方法.地球物理學(xué)報，58(9):3335-3345,

10.6038/cjg20150926.

Lu X T, Han L G, Zhang P，et al. 2015. Direct migration method of multi-source blended data based on total variation.ChineseJ.Geophys. (in Chinese),58(9):3335-3345,doi:10.6038/cjg20150926.

10.6038/cjg20150926

P631

2014-12-31，2015-09-01收修定稿

國家自然科學(xué)基金項目(41374115)，國家高技術(shù)研究發(fā)展計劃(863計劃)重大項目課題(2014AA06A605)聯(lián)合資助.

盧昕婷，女，1988年生，博士生，主要從事混合震源數(shù)據(jù)處理.E-mail:luxt13@mails.jlu.edu.cn

*通訊作者 韓立國，男，1961年生，教授，博士生導(dǎo)師，主要從事地震數(shù)據(jù)處理解釋工作.E-mail:hanliguo@jlu.edu.cn