同時(shí)震源數(shù)據(jù)的直接反演分離

王坤喜 毛偉建 張慶臣 李武群 詹 毅 孫鄖松

(①中國科學(xué)院測量與地球物理研究所計(jì)算與勘探地球物理研究中心，湖北武漢 430077； ②大地測量與地球動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，湖北武漢 430077； ③中國科學(xué)院大學(xué)，北京 100049； ④東方地球物理公司物探技術(shù)研究中心，河北涿州 072751)

0 引言

地震數(shù)據(jù)采集需要考慮成本與質(zhì)量之間的平衡。傳統(tǒng)的地震勘探野外數(shù)據(jù)采集過程中，相鄰炮的激發(fā)時(shí)間間隔較大，采集效率低，成本高。同時(shí)震源(simultaneous source)采集技術(shù)，也稱多震源地震數(shù)據(jù)采集技術(shù)，改變了傳統(tǒng)的采集方式。當(dāng)震源數(shù)量一定時(shí)，同時(shí)震源采集能夠縮短采集時(shí)間，增大覆蓋范圍；通過增加震源數(shù)量，提高覆蓋次數(shù)，改善對地下構(gòu)造的照明質(zhì)量[1-3]。

同時(shí)震源數(shù)據(jù)采集能夠提高采集效率和改善照明質(zhì)量，但是相鄰震源之間會(huì)產(chǎn)生“串?dāng)_噪聲”，因此同時(shí)震源數(shù)據(jù)不適用于傳統(tǒng)的處理流程。目前此類數(shù)據(jù)的處理方法主要有兩大類： ①直接成像，即直接用同時(shí)震源數(shù)據(jù)進(jìn)行偏移成像，這類方法受到串?dāng)_噪聲的干擾，成像質(zhì)量會(huì)降低[4-10]； ②分離法，即先對混合數(shù)據(jù)進(jìn)行分離，得到類似單個(gè)震源激發(fā)的數(shù)據(jù)，再用常規(guī)地震方法進(jìn)行處理。

分離法進(jìn)一步劃分為濾波類和反演類。從共炮域變換到其他域，在偽分離記錄中主震源信號(hào)是連續(xù)的，混合噪聲不再具有連續(xù)相干性。據(jù)此，Moore等[11]和Akerberg等[12]在共偏移距域進(jìn)行濾波； Huo等[13]利用多方向矢量中值濾波法對地層傾角進(jìn)行掃描，求得最佳傾角后，再利用中值濾波去除共中心點(diǎn)域的混合噪聲。Yang等[14]提出先利用PWD(plane-wave destruction)求取地層的傾角信息，再結(jié)合矢量中值濾波對串?dāng)_噪聲進(jìn)行壓制。反演類方法的思路是施加約束條件，將分離問題變?yōu)榉囱輪栴}。Doulgeris等[15]設(shè)計(jì)了迭代算子，在共檢波點(diǎn)域壓制噪聲； Mahdad等[16]在頻率—波數(shù)域分選偽分離數(shù)據(jù)，設(shè)計(jì)迭代算法，分離混合波場；韓立國等[17]結(jié)合多級(jí)中值濾波和Curvelet閾值迭代去噪算法，實(shí)現(xiàn)了混合數(shù)據(jù)的分離； Chen等[18]、祖紹環(huán)等[19]、Xue等[20]和宋家文等[21]分別將地震數(shù)據(jù)轉(zhuǎn)換到Seislet域、Curvelet域、Radon域和頻率—波數(shù)—波數(shù)域(FKK)，然后在稀疏域?qū)旌系卣饠?shù)據(jù)進(jìn)行正則化稀疏約束迭代，實(shí)現(xiàn)了混合數(shù)據(jù)的分離； Zu等[22]引入兩個(gè)卷積算子構(gòu)造逆問題，先求得地震信號(hào)的傾角，然后通過共軛梯度算法求解該逆問題，實(shí)現(xiàn)混疊數(shù)據(jù)的分離； Zhou等[23]在共檢波點(diǎn)域利用秩縮減濾波器和相干濾波器，通過迭代壓制非相干信號(hào)，也取得了不錯(cuò)的效果。

以上分離方法都要施加約束條件，比如地震信號(hào)的相干性和稀疏性。Wapenaar等[24-26]在地震多維反褶積干涉的基礎(chǔ)上，對同時(shí)震源數(shù)據(jù)的分離進(jìn)行了研究，認(rèn)為可以不施加這些額外的約束條件。在時(shí)間延遲算子的偽逆中，增加了基礎(chǔ)點(diǎn)擴(kuò)散矩陣，在滿足一定的分離條件下，將混合數(shù)據(jù)按延遲時(shí)間進(jìn)行歸位，直接得到分離結(jié)果。

本文通過推廣基礎(chǔ)點(diǎn)擴(kuò)散矩陣，將檢波點(diǎn)和炮點(diǎn)規(guī)則排列采集的同時(shí)震源數(shù)據(jù)分離方法拓展應(yīng)用于不規(guī)則排列采集的同時(shí)震源數(shù)據(jù)分離。

1 方法原理

1.1 數(shù)據(jù)混合的編碼方式

(1)

(2)

式中ω為角頻率。對每一個(gè)角頻率分量，在頻率域的混合炮數(shù)據(jù)可用矩陣的形式(圖1)表示[1]為

Psim=PΓ

(3)

圖1 混合編碼過程(每個(gè)頻率分量)

1.2 混合數(shù)據(jù)的偽分離

僅從方程的角度看，同時(shí)震源數(shù)據(jù)的分離就是通過式(3)求解原始數(shù)據(jù)矩陣P，所以形式上可以得到分離結(jié)果

(4)

(5)

(6)

式(6)表明偽分離法是先復(fù)制n份混合數(shù)據(jù)，再按主炮的延遲時(shí)間歸位[17]，但仍然保留其他炮的信號(hào)，形成串?dāng)_噪聲。

1.3 構(gòu)建基礎(chǔ)點(diǎn)擴(kuò)散矩陣

在偽分離中，干擾炮按照主炮的延遲時(shí)間處理，而非按照自己的延遲時(shí)間歸位，形成所謂的“串?dāng)_噪聲”。不同于常規(guī)加約束的迭代反演方法，這里不是簡單地將這些干擾炮的“串?dāng)_”看成噪聲，而是認(rèn)為這些“串?dāng)_噪聲”沒有得到正確歸位。通過討論同時(shí)震源數(shù)據(jù)的分離與多維反褶積地震干涉技術(shù)在原理上的相似性，引入基礎(chǔ)點(diǎn)擴(kuò)散矩陣對混合數(shù)據(jù)進(jìn)行處理[25-26]，將“串?dāng)_噪聲”按各自的延遲時(shí)間進(jìn)行歸位。

地震干涉技術(shù)可以將任意兩個(gè)檢波器接收到的數(shù)據(jù)合成為在若干檢波器之間傳播的波，相當(dāng)于將其中的一個(gè)檢波器看作虛擬震源(圖2)。

圖2 地震干涉的卷積模型

(7)

(8)

(9)

(10)

(11)

式中：Cseq(xB,xA,t)為相關(guān)函數(shù)；Rseq(x,xA,t)為點(diǎn)擴(kuò)散函數(shù)。其表達(dá)式分別為

(12)

(13)

(14)

(15)

在式(9)兩邊同時(shí)與s(i)(t-ti)褶積，并對超級(jí)炮σ(m)內(nèi)的所有單炮求和，得到類似于式(10)的表達(dá)式

(16)

對界面T上選取的一個(gè)固定參考檢波器xA，其記錄的波場為uin(xA,σ(m),t)，在式(16)兩端同時(shí)與uin(xA,σ(m),t)做互相關(guān)計(jì)算，并對所有超級(jí)炮σ(m)求和，則式(11)可以改寫為同時(shí)震源模型下的相關(guān)函數(shù)

(17)

式中：Csim(xB,xA,t)代表偽分離后的數(shù)據(jù)；Rsim(x,xA,t)為點(diǎn)擴(kuò)散函數(shù)。其表達(dá)式分別為

uin(xA,σ(m),-t)]

(18)

uin(xA,σ(m),-t)]

(19)

將式(9)變換到頻率域，每個(gè)頻率分量都有如下離散化后的矩陣表達(dá)式

(20)

式(14)和式(15)表示同時(shí)震源模型中的波場響應(yīng)，其離散矩陣為

Uin=GinΓ

(21)

Uout=GoutΓ

(22)

在式(20)兩端同乘以Γ，得到式(16)的離散矩陣表達(dá)式

(23)

(24)

根據(jù)式(19)，Rsim在頻率域的矩陣表達(dá)式為

Rsim=(Uin)?Uin

(25)

將式(21)代入式(25)，得到

Rsim=Γ?(Gin)?GinΓ=Γ?RΓ

(26)

式中(Gin)?Gin被定義為基礎(chǔ)點(diǎn)擴(kuò)散矩陣(the basic point-spread function)，用符號(hào)R表示

R=(Gin)?Gin

(27)

(28)

式中：xi為第i個(gè)震源到第1個(gè)震源的距離；γ(xi,ω)與二維偏移中的分辨率函數(shù)類似[24,27]，與sinc插值函數(shù)成正比[25,28]。在信號(hào)重建中，重建函數(shù)[29]

(29)

表示原始信號(hào)g(x)是待重建信號(hào)gd(x)與sinc(x)的卷積，由此將sinc(x)離散化。同理將γ(xi,ω)離散化并表示為矩陣R

(30)

式中γi=γ(xi,ω)=γ(iΔs,ω)，Δs為炮間距。

(31)

將式(24)代入上式，并用Gin替代W+[26]，得到

(32)

其中W-Uout為混合后的矩陣Psim，分別將式(21)和式(27)代入上式，得到同時(shí)震源的直接反演分離公式

=Psim[(GinΓ)?GinΓ]-1(GinΓ)?Gin

=Psim[Γ?(Gin)?GinΓ]-1Γ?(Gin)?Gin

=Psim(Γ?RΓ)-1Γ?R

(33)

(34)

式中b1=exp(-iωt1),b2=exp(-iωt2)。將式(30)和式(34)代入式Rsim=Γ?RΓ中，得到

(35)

式中Δt=t2-t1。式(35)中右邊第一項(xiàng)是式(30)中基礎(chǔ)點(diǎn)擴(kuò)散矩陣R的重采樣變形式(采樣因子為2)，第二和第三項(xiàng)也是基礎(chǔ)點(diǎn)擴(kuò)散矩陣的重采樣變形式，但是對距離的改變量為±Δs，在空間上對串?dāng)_噪聲進(jìn)行了歸位；因子exp(?iωΔt)表示時(shí)間改變量±Δt，對串?dāng)_噪聲在時(shí)間上進(jìn)行了歸位。式中的第二和第三項(xiàng)解釋了對同時(shí)震源數(shù)據(jù)中的串?dāng)_噪聲的處理方式。在對上式的分析中，假設(shè)一個(gè)超級(jí)炮是由相鄰兩個(gè)單炮混合組成的。如果一個(gè)超級(jí)炮由n個(gè)相鄰震源組成，也能夠得到類似于式(35)的表達(dá)式，只是基礎(chǔ)點(diǎn)擴(kuò)散矩陣R的重采樣變形式在空間和時(shí)間的改變量會(huì)隨之發(fā)生變化[26]。

直接反演分離的具體流程如圖3所示。

1.4 分離適用條件

直接反演分離方法利用地震數(shù)據(jù)帶寬有限的特點(diǎn)，使用與sinc插值函數(shù)類似的γ(xi,ω)時(shí)，要求直接反演方法滿足空間采樣定理，不能產(chǎn)生空間假頻，這意味著混合度n不能夠無限制地增大。

根據(jù)遠(yuǎn)場理論，在共炮域，一個(gè)頻率為f的平面簡諧波入射到地面檢波器上，空間假頻出現(xiàn)的截止頻率為

(36)

圖3 同時(shí)震源數(shù)據(jù)的直接反演分離流程圖

式中：v為近地表地震波的傳播速度；Δx為檢波器間距；θ為地震波的傳播角度。根據(jù)波場互換原理，把共檢波點(diǎn)道集中的檢波點(diǎn)視為炮點(diǎn)，炮點(diǎn)視為檢波點(diǎn)，則允許截止頻率為

(37)

在同時(shí)震源數(shù)據(jù)采集中，相鄰n個(gè)單炮幾乎同時(shí)激發(fā)，形成一個(gè)超級(jí)炮，超級(jí)炮的間距為

ΔS=nΔs

(38)

由全部超級(jí)炮形成的混合波場中，允許的截止頻率為

(39)

實(shí)際震源激發(fā)的頻率f應(yīng)不大于截止頻率

(40)

可得

(41)

當(dāng)滿足式(41)時(shí)，空間假頻就不會(huì)出現(xiàn)，才能使用γ(xi,ω)，混合震源才會(huì)被正確解混為單個(gè)震源。

2 檢波點(diǎn)與炮點(diǎn)的不規(guī)則排列

在前面原理部分，待分離的同時(shí)震源數(shù)據(jù)體是通過固定炮間距和固定檢波點(diǎn)間距的觀測系統(tǒng)混合采集得到。在實(shí)際生產(chǎn)中，遇到無法避開的障礙物時(shí)，不可避免地要改變檢波點(diǎn)和炮點(diǎn)間距，形成不規(guī)則排列。為了使直接分離算法得到廣泛應(yīng)用，對檢波點(diǎn)和炮點(diǎn)都不規(guī)則排列的采集條件進(jìn)行研究。

2.1 檢波點(diǎn)不規(guī)則排列

(42)

2.2 炮點(diǎn)不規(guī)則排列

(43)

在炮點(diǎn)不規(guī)則排列的觀測系統(tǒng)中，防止出現(xiàn)空間假頻的分離適用條件(式(41))變?yōu)?/p>

(44)

式中Δsi為第i炮與第i+1炮之間的距離， Δs1+Δs2+Δs3+…+Δsn代表一個(gè)超級(jí)炮在地面的覆蓋范圍。

3 算列

3.1 規(guī)則排列下的模擬數(shù)據(jù)

圖4為傳播速度具有橫向和垂直變化的Marmousi模型。模型的最低速度為1500m/s，正演的震源截止頻率為25Hz。在地表創(chuàng)建765個(gè)炮點(diǎn)，炮間距為10m，以3個(gè)相鄰單炮組成一個(gè)超級(jí)炮，即混合度n=3，因此共有255個(gè)超級(jí)炮，每個(gè)單炮的延遲時(shí)間是隨機(jī)的，取值為0～2s； 用767個(gè)檢波器接收信號(hào)，道間距為10m，采樣間隔為4ms。模擬計(jì)算在Intel Xeon E7-4807、內(nèi)存為128G的IBM服務(wù)器上運(yùn)行，使用的編程語言為MATLAB。混合數(shù)據(jù)體(其中部分道集如圖5所示)經(jīng)過一次處理可以得到全部解混后的數(shù)據(jù)。

圖4 Marmousi速度模型

如果僅采用最小二乘法進(jìn)行分離，只對主炮進(jìn)行了歸位。圖5e為在共檢波點(diǎn)道集中用最小二乘法求得的分離結(jié)果。對比圖5d，可以看出主炮已經(jīng)歸位，但是其他干擾炮形成的“串?dāng)_”特征依然存在，這些“串?dāng)_”數(shù)據(jù)攜帶其他干擾炮的信息，并不能簡單地理解為噪聲。所以在每個(gè)頻率分量中，首先構(gòu)建基礎(chǔ)點(diǎn)擴(kuò)散矩陣R(圖6)，然后計(jì)算點(diǎn)擴(kuò)散矩陣Rsim。在共檢波點(diǎn)域，矩陣Rsim把這些“串?dāng)_”在時(shí)間上按各炮的延遲時(shí)間進(jìn)行歸位，在空間上按炮間距進(jìn)行歸位，因此混疊的同時(shí)震源數(shù)據(jù)得到了分離。在實(shí)際計(jì)算中，為了提高矩陣Rsim求逆的穩(wěn)定性，一般在矩陣Rsim中增加正則化參數(shù)β，即Rsim=Γ?RΓ+βI，試驗(yàn)證明β=1.0×10-6ε是一個(gè)合理的選擇，其中ε是矩陣Rsim中的最大元素。

圖5 規(guī)則排列下不同時(shí)域的原始采集數(shù)據(jù)與混合采集數(shù)據(jù)

圖6 頻率分量12Hz處的基礎(chǔ)點(diǎn)擴(kuò)散矩陣R

根據(jù)圖3完成后續(xù)步驟，得到全部765炮的分離數(shù)據(jù)。圖5中混合采集數(shù)據(jù)的分離結(jié)果見圖7。

定義如下信噪比公式定量描述數(shù)據(jù)分離的效果

(45)

本次實(shí)驗(yàn)中，共炮點(diǎn)數(shù)據(jù)(圖5a與圖7a)和共檢波點(diǎn)數(shù)據(jù)(圖5b與圖7b)的分離信噪比SNR分別為37.8dB和35.3dB。抽取第90道和第200道的原始數(shù)據(jù)(圖5b)與分離數(shù)據(jù)(圖7b)進(jìn)行對比，結(jié)果見圖8。

在前述的計(jì)算機(jī)環(huán)境下，一次性分離全部混合數(shù)據(jù)體(由255個(gè)超級(jí)炮，767個(gè)檢波點(diǎn)和每道1200個(gè)采樣點(diǎn)組成)耗時(shí)381s，用時(shí)較短。對比分離后數(shù)據(jù)(圖7a、圖7b)與原始數(shù)據(jù)(圖5a、圖5b)，可以看出，不僅淺層數(shù)據(jù)得到了很好的恢復(fù)，深部細(xì)節(jié)也得到保留，殘差較小，信噪比較高，且兩條波形曲線的吻合度很高，這證明了本文方法的準(zhǔn)確性。

3.2 規(guī)則排列下的實(shí)際數(shù)據(jù)

實(shí)際數(shù)據(jù)的炮間距為20m，采樣間隔為1ms，觀測時(shí)間為2.3s，混合度n=2?？紤]到截止頻率對本方法的影響較大，對實(shí)際數(shù)據(jù)濾波，去掉了高頻分量。分離結(jié)果及殘差見圖9c、圖9d?？梢钥闯?，雖然中間記錄道的強(qiáng)“串?dāng)_噪聲”沒能完全被壓制，但是分離后的深部反射層的同相軸信號(hào)得到了較好的恢復(fù)，最終的SNR為10.5dB。

圖7 規(guī)則排列下不同時(shí)域的分離數(shù)據(jù)與殘差

圖8 單道分離數(shù)據(jù)與原始采集數(shù)據(jù)的波形對比圖

3.3 不規(guī)則排列的模擬數(shù)據(jù)

不規(guī)則排列的正演速度模型仍為Marmousi模型(圖4)。炮點(diǎn)數(shù)為765，炮間距隨機(jī)，但是需要滿足分離適用條件(式(44))。檢波點(diǎn)數(shù)為150，分為3段，每段包含50個(gè)檢波點(diǎn)，分別布置在地面500～1500m、3000～4000m和5500～6500m，道間距也隨機(jī)取值。炮點(diǎn)與檢波點(diǎn)覆蓋范圍如圖10所示，混合度n=3，其他條件與前述規(guī)則排列正演相同。分別討論當(dāng)檢波點(diǎn)或炮點(diǎn)不規(guī)則排列下，同時(shí)震源數(shù)據(jù)的直接反演分離效果。

首先，針對檢波點(diǎn)不規(guī)則排列得到的同時(shí)震源數(shù)據(jù)進(jìn)行分離測試(圖11)，分離后的信噪比SNR=39.5dB。然后，針對炮點(diǎn)不規(guī)則排列的情況進(jìn)行了測試(圖12)，分離的信噪比SNR=38.1dB。對比分離結(jié)果(圖11c、圖12c)與原始采集數(shù)據(jù)(圖11a、圖12a)，可以看出，混合波場從淺層到深層都得到了很好的分離，在分離的殘差(圖11d、圖12d)中幾乎看不到連續(xù)的有效反射信號(hào)。證明本文方法不僅適用于規(guī)則排列采集的同時(shí)震源數(shù)據(jù)的分離，通過拓展后同樣適用于不規(guī)則排列采集的同時(shí)震源數(shù)據(jù)的分離。

圖9 規(guī)則排列下實(shí)際共檢波點(diǎn)道集數(shù)據(jù)的直接反演分離

圖10 炮點(diǎn)與檢波點(diǎn)在地面的覆蓋范圍示意圖

圖11 檢波點(diǎn)不規(guī)則排列條件下共炮點(diǎn)道集數(shù)據(jù)的直接反演分離

圖12 炮點(diǎn)不規(guī)則排列下共檢波點(diǎn)道集數(shù)據(jù)的直接反演分離結(jié)果

4 討論

本文在求取矩陣Rsim=Γ?RΓ+βI的逆的過程中采用了正則化方法，該求逆過程會(huì)引入較多的假頻噪聲。通過對正則化參數(shù)的研究，選擇合適的參數(shù)β，可以將這種影響降到最低。圖13展示了取不同正則化參數(shù)β時(shí)，規(guī)則排列的Marmousi混合數(shù)據(jù)體中，第100個(gè)檢波器記錄的共檢波點(diǎn)道集的分離信噪比變化曲線?？梢钥闯霎?dāng)參數(shù)β取值較大或者較小時(shí)，信噪比都會(huì)降低，假頻噪聲增強(qiáng)，分離效果較差； 當(dāng)β取1.0×10-4ε～1.0×10-8ε時(shí)，可以得到較高的信噪比。

算例中，對于檢波點(diǎn)與炮點(diǎn)規(guī)則排列采集的Marmousi正演混合數(shù)據(jù)體，分離后可以得到765個(gè)共炮點(diǎn)道集和767個(gè)共檢波點(diǎn)道集，這些共炮點(diǎn)道集與共檢波點(diǎn)道集各自的分離信噪比見圖14?？梢钥闯?，靠近觀測系統(tǒng)中部采集的混合數(shù)據(jù)分離效果較好，處于兩端的分離效果相對較差，具體原因還需進(jìn)一步研究。

圖13 規(guī)則排列下的模擬數(shù)據(jù)算列中正則化參數(shù)β對分離SNR的影響

圖14 Marmousi模型檢波點(diǎn)與炮點(diǎn)規(guī)則排列下數(shù)據(jù)體的全部共炮點(diǎn)道集與共檢波點(diǎn)道集分離SNR曲線

5 結(jié)論

一般而言，進(jìn)行高密度震源采樣時(shí)，若震源間距較小、激發(fā)的延遲時(shí)間很短，激發(fā)波場容易產(chǎn)生強(qiáng)相干性，這不利于同時(shí)震源數(shù)據(jù)的分離。本文采用的直接反演方法則不受這些條件限制，對分離激發(fā)延遲時(shí)間短、炮間距較小的高密度同時(shí)震源數(shù)據(jù)具有一定的優(yōu)勢。該方法只需在頻率域進(jìn)行一次性分離，不需要迭代，效率高、準(zhǔn)確度高，對同時(shí)震源數(shù)據(jù)的分離有很好的應(yīng)用前景。

通過適當(dāng)改變基礎(chǔ)點(diǎn)擴(kuò)散矩陣，該方法被推廣到檢波點(diǎn)和炮點(diǎn)不規(guī)則排列的觀測系統(tǒng)中，不僅滿足陸地復(fù)雜地形的同時(shí)震源數(shù)據(jù)的分離，而且適用于海上拖纜的動(dòng)態(tài)不規(guī)則采集的混合數(shù)據(jù)的分離。

需要注意的是，直接反演分離算法有一定的適用范圍，要求混合度、炮間距和地震信號(hào)的頻率滿足本文中的分離適用條件。