崔慶輝， 芮擁軍， 秦寧， 劉鵬翔， 隆文韜

(中國石油化工股份有限公司勝利油田分公司物探研究院， 東營 257022)

隨著社會經(jīng)濟(jì)的快速發(fā)展，當(dāng)前在東部老油區(qū)城鎮(zhèn)、工廠、村莊、道路等障礙物密集分布，地震勘探的地表條件發(fā)生了巨大變化，造成基于規(guī)則采樣理論即Shannon采樣定理[1]設(shè)計的觀測系統(tǒng)受密集障礙物影響，在實際布設(shè)時點位發(fā)生嚴(yán)重偏離，特別是炮點位置變化更大，呈現(xiàn)出隨機分布特點，導(dǎo)致成像結(jié)果偏離預(yù)期。另一方面，一些老油田勘探開發(fā)程度較高，需要解決的地質(zhì)目標(biāo)表現(xiàn)出“深、小、薄、碎”等特征，在常規(guī)規(guī)則采樣理論下只能不斷提高采樣密度即炮道密度，這無疑會造成采集成本的急劇上升。

壓縮感知理論是一種新的采樣理論，基于信號的稀疏特性，在遠(yuǎn)小于Nyquist采樣率的條件下，通過隨機采樣獲取信號的離散樣點值，然后進(jìn)行信號的重建[2-3]。壓縮感知理論一經(jīng)提出，就引起學(xué)術(shù)界和工業(yè)界的廣泛關(guān)注，在各個領(lǐng)域得到迅速發(fā)展。

國外較早將壓縮感知理論應(yīng)用于海上地震勘探，取得了較好的應(yīng)用效果[4-5]。馬堅偉[6-7]較早研究了壓縮感知理論在地震勘探領(lǐng)域應(yīng)用的可行性，周松等[8]和呂公河等[9]在西部沙漠地區(qū)開展了基于壓縮感知的地震采集試驗，為中國陸上基于壓縮感知的地震勘探實際應(yīng)用奠定了基礎(chǔ)。黃小剛[10-11]提出了一套海上壓縮感知采集設(shè)計及數(shù)據(jù)重構(gòu)方法，并用模型數(shù)據(jù)展示了海上壓縮感知采集的效果?？傊诋?dāng)前低油價背景下，壓縮感知理論在地震勘探領(lǐng)域的研究得到了廣泛的重視[12]。

利用壓縮感知理論實現(xiàn)隨機地震采樣并重構(gòu)規(guī)則數(shù)據(jù)，首先要分析地震數(shù)據(jù)的稀疏性，地震信號(或波場)的時空復(fù)雜度(與地下及地表速度結(jié)構(gòu)、噪音特征等相關(guān))會顯著影響地震信號的稀疏性，進(jìn)而導(dǎo)致很多情形下隨機采集的實際數(shù)據(jù)的重建效果無法達(dá)到規(guī)則采樣數(shù)據(jù)的質(zhì)量[13]。中國陸上油田地表噪聲發(fā)育，地下波場復(fù)雜，基于壓縮感知理論進(jìn)行隨機地震數(shù)據(jù)采集及規(guī)則化重建具有更大難度。

首先研究基于壓縮感知理論的隨機地震觀測系統(tǒng)設(shè)計方法及基于地震數(shù)據(jù)高維空間低秩約束的隨機地震數(shù)據(jù)規(guī)則化重建方法，并通過三維模型數(shù)值模擬證明了以上方法在同等采樣密度下比常規(guī)規(guī)則采樣具有更好的成像效果。

1 隨機地震觀測系統(tǒng)設(shè)計基本理論

隨機地震采集需要滿足三個條件：地震數(shù)據(jù)具有稀疏性、合理的隨機采樣方法和有效的數(shù)據(jù)規(guī)則化重建算法。關(guān)于地震數(shù)據(jù)的稀疏性，很多學(xué)者的研究表明實際地震數(shù)據(jù)在時間域或其他域具有稀疏特征，即地震數(shù)據(jù)可稀疏表達(dá)[13-15]。現(xiàn)重點研究隨機地震觀測系統(tǒng)設(shè)計及數(shù)據(jù)規(guī)則化方法。

1.1 隨機地震采樣矩陣

野外地震數(shù)據(jù)觀測實際上是對信號的一種采樣，地震數(shù)據(jù)采集過程可以表示為

d=Df

(1)

式(1)中：f∈Rm為原始地震數(shù)據(jù)；d∈Rn為經(jīng)過采樣得到的地震數(shù)據(jù)；D∈Rn×m為采樣矩陣；通常有n

若f本身不稀疏，則可以使用f在變換域中的稀疏性來得到一種等價的表示。這里用C表示f的稀疏變換，則式(1)就可以轉(zhuǎn)換為

d=Ak

(2)

A=DCH

(3)

式中：H代表反變換矩陣；k∈Rn為f在稀疏域中的系數(shù)；n為變換域的維度(若C為正交變換，則n=m；若C為冗余變換，n>m)，A∈Rn×n。

要滿足A為隨機矩陣的要求，在變換矩陣C確定的情況下，采樣矩陣R也必須是隨機分布矩陣，那么就可以利用采樣數(shù)據(jù)高概率地恢復(fù)完整地震數(shù)據(jù)。

1.2 隨機地震觀測系統(tǒng)設(shè)計方法

通過高斯隨機函數(shù)進(jìn)行完全的隨機采樣，被基于壓縮感知的技術(shù)廣泛采用。但是，這種方法容易造成采樣點高度聚集或者高度分散的情況，導(dǎo)致對信號的某些部分采樣冗余或者采樣不足，不僅不符合高精度地震數(shù)據(jù)采集的要求，而且也會給后續(xù)數(shù)據(jù)恢復(fù)處理帶來困難，難以達(dá)到理想的重建效果。

計算機圖像學(xué)領(lǐng)域應(yīng)用的泊松碟隨機采樣方法[16]、抖動采樣方法(Jittered)[17]以及分段采樣方法不僅能保證測點的隨機性，還能保證測點的空間均勻分布，并具有明顯的“藍(lán)色噪聲”頻譜特征[14]?！八{(lán)色噪聲”頻譜特征是描述隨機類型特征的一個統(tǒng)計模型，其特點是在主要頻率周圍分布的其他頻率值都非常小，所以真值或者主頻率就很容易被檢測到。

結(jié)合實際地震數(shù)據(jù)特點，基于Jittered采樣原理提出了以下隨機地震觀測系統(tǒng)設(shè)計技術(shù)流程。

步驟1按照常規(guī)的三維地震觀測系統(tǒng)設(shè)計方法確定理想的觀測參數(shù)，在正交觀測系統(tǒng)下可以用炮點和檢波點的點距、線距、點數(shù)和線數(shù)表示，即{x0,y0,Δx,Δy,nx,ny}，其中(x0,y0)為采樣點的起始坐標(biāo)，Δx和Δy分別為采樣點之間的橫向和縱向間隔，nx和ny分別為橫向和縱向采樣點數(shù)。基于以上參數(shù)可得到規(guī)則的采樣點坐標(biāo){x,y}。

步驟2根據(jù)工區(qū)地表及地下復(fù)雜程度設(shè)置縱橫向壓縮采樣比例{rx,ry}，基于該壓縮采樣比例計算隨機采樣點坐標(biāo)的期望值{xe,ye}。

步驟3設(shè)置隨機采樣點坐標(biāo)相對于期望值的擾動量，即高斯隨機分布的標(biāo)準(zhǔn)差(σx,σy)，一般取(σx,σy)=(3Δx,3Δy)。

最后，根據(jù)步驟2中期望值，以及步驟3中標(biāo)準(zhǔn)差利用隨機函數(shù)生成隨機觀測網(wǎng)格{xrand,yrand}。

2 高維空間低秩約束的隨機地震數(shù)據(jù)重建

地震數(shù)據(jù)的偏移對數(shù)據(jù)規(guī)則性要求很高，實際采集的地震數(shù)據(jù)受觀測系統(tǒng)、地表變觀等影響需要進(jìn)行規(guī)則化處理，因此產(chǎn)生了很多的地震數(shù)據(jù)插值算法[18-19]，主要包括基于函數(shù)變換的重建、基于波場延拓的重建、低秩約束的地震數(shù)據(jù)重建。低秩約束的地震數(shù)據(jù)重建算法應(yīng)用最為廣泛，其基本原理為：完整地震數(shù)據(jù)在頻率域可以被一個低秩矩陣或者數(shù)組表示，不規(guī)則缺道數(shù)據(jù)被看作隨機噪聲，而隨機噪聲會增加矩陣或者高維數(shù)組的秩。地震數(shù)據(jù)重建就可以歸結(jié)為一個高秩矩陣或者數(shù)組的降秩問題，通過降秩實現(xiàn)重建。

地震數(shù)據(jù)中的反射波同相軸在局部表現(xiàn)為線性特征，利用數(shù)據(jù)的局部相似性，抽取局部五維數(shù)據(jù)，在高維空間進(jìn)行“稀疏-低秩聯(lián)合反演重建”，從而在壓縮感知理論下保證了重建數(shù)據(jù)的精度，算法主要流程如圖1所示。

圖1 低秩約束的高維空間地震數(shù)據(jù)重建算法流程Fig.1 Low-rank constrained high-dimensional spatial seismic data reconstruction algorithm flow

步驟1定義一個理想的規(guī)則工區(qū)網(wǎng)格G，按照規(guī)則觀測系統(tǒng)設(shè)計方法[20]在網(wǎng)格G上定義規(guī)則觀測系統(tǒng)，記為Ωreg。

步驟2獲取原始輸入數(shù)據(jù)集Dirr的非規(guī)則地震數(shù)據(jù)的觀測系統(tǒng)。以炮集為單位獲取炮檢點坐標(biāo)，取整到臨近網(wǎng)格，得到Dirr對應(yīng)的非規(guī)則觀測系統(tǒng)Ωirr。

步驟3以窗口滑動方式從原始數(shù)據(jù)集取得一個局部高維數(shù)據(jù)子集。通過在炮端檢端分別取局部窗口，抽取原始道集中落在窗口中的所有地震道，從而獲得一個不完備的五維局部數(shù)據(jù)體。通過滑動中心炮檢點，可以循環(huán)獲取不同位置的局部高維數(shù)據(jù)。

步驟4在步驟3獲取的不完備五維局部數(shù)據(jù)體上，在頻率域進(jìn)行規(guī)則數(shù)據(jù)重建，輸出一個完備且規(guī)則的局部五維局部數(shù)據(jù)體。計算過程即求解優(yōu)化問題，即

(4)

該優(yōu)化問題實現(xiàn)了對重建結(jié)果的三種約束。

(1)要求在觀測點位上的重建結(jié)果與原始觀測相吻合(L2項)。

(2)要求重建數(shù)據(jù)在四維波數(shù)具有稀疏性(L1項)。

(3)要求重建數(shù)據(jù)在進(jìn)行Hankel排列后具有低秩性(核范數(shù)項)。

式(4)所示優(yōu)化問題的求解，可通過以下迭代算法實現(xiàn)。

(1)關(guān)于L2項的梯度更新為

(5)

(6)

(7)

(2)關(guān)于L1項的軟閾值收縮為

(8)

(3)關(guān)于核范數(shù)項的奇異值收縮：

(9)

式(9)中：H(·)和H-1(·)分別是Hankel矩陣化正、反變換算子；Sμ(x)表示對x進(jìn)行奇異值軟閾值收縮，閾值取μ。

步驟5對步驟4重建的局部五維規(guī)則數(shù)據(jù)體做余弦窗函數(shù)加權(quán)，然后以累加方式更新到規(guī)則數(shù)據(jù)集Dreg中，并記錄輸出數(shù)據(jù)集中各個地震道的累加總權(quán)重Csum(xs,xr)。

更新規(guī)則數(shù)據(jù)集之前需要對重建數(shù)據(jù)的余弦加權(quán)操作，即每一道都需要乘以一個標(biāo)量系數(shù)，對于震源位于xs，檢波器位于xr的地震道d(xs,ys,xr,yr,t) ，它對應(yīng)的加權(quán)系數(shù)C計算公式為

(10)

(11)

同時，對式(10)進(jìn)行更新，每一次更新之后都需要將各個地震道的加權(quán)系數(shù)累加到總權(quán)重數(shù)組Csum中，以便后續(xù)消除加權(quán)幅值不均勻?qū)е碌幕儭８鹿綖?/p>

(12)

步驟6重復(fù)步驟3～步驟5，直到完成遍歷整個數(shù)據(jù)集。通過滑動局部窗遍歷數(shù)據(jù)的過程中，保持?jǐn)?shù)據(jù)窗尺寸不變，各方向的滑動步長取相應(yīng)方向的窗口尺寸的一半。

步驟7通過除以各個地震道的累加總權(quán)重Csum(xs,xr)，將計算得到的規(guī)則數(shù)據(jù)集Dreg權(quán)重歸一化。將重建的規(guī)則數(shù)據(jù)的各個地震道Dreg(xs,xr,t)除以其相應(yīng)的累加總幅值Csum(xs,xr)，以消除加權(quán)幅值不均勻?qū)е碌幕儭?/p>

3 模型數(shù)值模擬

選取東部某油田一個三維模型(圖2)，通過炮集數(shù)據(jù)的數(shù)值模擬及疊前深度偏移成像，驗證文中隨機地震觀測系統(tǒng)設(shè)計及數(shù)據(jù)重建算法的有效性。

圖2 東部油田某工區(qū)三維地質(zhì)模型Fig.2 Three-dimensional geological model of a work area in the eastern oilfield

首先設(shè)計如圖3(a)所示的三維規(guī)則觀測系統(tǒng)，炮點橫向間隔50 m?？v向間隔400 m，檢波點橫向間隔400 m，縱向間隔50 m?？紤]到目前實際生產(chǎn)中單方向隨機采集可行性最高，因此本次隨機觀測系統(tǒng)設(shè)計對炮點只在橫向上按照壓縮率50%進(jìn)行隨機采集設(shè)計[圖3(b)]，對檢波點只在縱向上按照壓縮率50%進(jìn)行隨機采集設(shè)計[圖3(c)]。

圖3 用于數(shù)值模擬的觀測系統(tǒng)Fig.3 Geometry used in numerical simulation

利用三維波動方程正演得到模擬的規(guī)則炮集和隨機炮集數(shù)據(jù)[圖4(a)]，分別加入高斯噪聲與椒鹽噪聲[圖5(a)和圖6(a)]，利用本文提出的高維空間低秩約束的隨機地震數(shù)據(jù)重建算法得到重建后的規(guī)則地震數(shù)據(jù)，對比重建地震數(shù)據(jù)的波數(shù)譜。在波數(shù)譜中可以看出，模型正演地震記錄的波數(shù)譜展布區(qū)域較大[圖5(a)]，表明了地震信號的復(fù)雜度較高。此外，無噪聲及加入不同噪聲情況下，隨機地震數(shù)據(jù)重建結(jié)果較好且波數(shù)譜也得到了較好的恢復(fù)，并未在波數(shù)譜中引入明顯的空間假頻(圖4～圖6)，充分展示了高維空間低秩約束的隨機地震數(shù)據(jù)重建的性能。

將規(guī)則的地震數(shù)據(jù)按照地震道均勻抽取25%，得到稀疏的規(guī)則地震數(shù)據(jù)，基于完全相同的處理流程和速度模型進(jìn)行疊前深度偏移，對比在同樣采集密度下稀疏規(guī)則采樣和隨機采樣得到的數(shù)據(jù)成像效果(圖7、圖8)。

圖4 無噪聲情況下的隨機地震數(shù)據(jù)重建Fig.4 Reconstruction of random seismic data without noise

圖5 加入高斯噪聲的隨機地震數(shù)據(jù)重建Fig.5 Reconstruction of random seismic data with gause noise

圖6 加入椒鹽噪聲的隨機地震數(shù)據(jù)重建效果Fig.6 Reconstruction of random seismic data with salt and pepper noise

稀疏規(guī)則觀測系統(tǒng)造成了成像質(zhì)量明顯降低[圖 7(a)]，而隨機采集重建數(shù)據(jù)的成像結(jié)果與理想數(shù)據(jù)的成像效果基本相當(dāng)，但成像噪聲稍強[圖7(b)]。

為對比目標(biāo)層成像結(jié)果(大致位于X=2.4～4.5 km,Y=1.2～2.2 km的區(qū)域)，進(jìn)一步放大目標(biāo)層成像結(jié)果，如圖8所示，可以看出，相對于稀疏規(guī)則數(shù)據(jù)成像結(jié)果，隨機采集數(shù)據(jù)經(jīng)過重建后成像結(jié)果的信噪比更高，且對薄層的分辨能力更好。

圖7 不同數(shù)據(jù)的成像對比Fig.7 Imaging comparison of different data

圖8 不同觀測系統(tǒng)數(shù)據(jù)的局部成像對比Fig.8 Zoomed in imaging comparison of different data

4 結(jié)論

(1)針對具體的地質(zhì)任務(wù)，先進(jìn)行常規(guī)的規(guī)則地震觀測系統(tǒng)設(shè)計，在此基礎(chǔ)上結(jié)合工區(qū)地表及地下復(fù)雜程度，對規(guī)則的采樣點數(shù)進(jìn)行一定比例的壓縮，基于Jittered采樣原理進(jìn)行隨機觀測系統(tǒng)設(shè)計，是目前較為可行的方法。

(2)在進(jìn)行隨機觀測系統(tǒng)設(shè)計時，炮點和檢波點的隨機方式可以靈活選擇，在復(fù)雜地區(qū)可對炮點和檢波點只進(jìn)行線內(nèi)隨機；對于簡單的地區(qū)，可以進(jìn)行完全隨機。

(3)本文提出的基于低秩約束的高維地震數(shù)據(jù)重建算法，充分利用地震數(shù)據(jù)在局部具有的線性特征，在不同背景噪聲下均能夠完成數(shù)據(jù)的無假頻重建，可進(jìn)一步應(yīng)用于實際地震數(shù)據(jù)重建。

(4)模型數(shù)值模擬結(jié)果證明，在同樣采集密度下，經(jīng)過規(guī)則化數(shù)據(jù)重建后，隨機地震數(shù)據(jù)比稀疏的規(guī)則采樣地震數(shù)據(jù)具有更好的成像效果，為進(jìn)一步在實際地震勘探中開展隨機采集試驗奠定了基礎(chǔ)。