基于雙曲Radon變換的混采地震數(shù)據(jù)重建

王晗，韓立國

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

0 引言

在常規(guī)地震勘探采集過程中,作業(yè)周期長,勘探成本較高，隨著多源地震混合采集技術(shù)概念的提出,勘探成本有效降低。但由于混合震源采集地震數(shù)據(jù)容易受到地形限制、施工條件、儀器故障等影響，采集的地震數(shù)據(jù)會出現(xiàn)缺失及不規(guī)則的情況。在實際數(shù)據(jù)處理過程中，數(shù)據(jù)缺失情況會對后續(xù)處理產(chǎn)生不良影響，特別是基于多道處理的算法流程[1-4]。而混采地震數(shù)據(jù)重建比常規(guī)地震數(shù)據(jù)重建更為復(fù)雜，需要對混采數(shù)據(jù)進行炮分離，對分離的單炮數(shù)據(jù)進行重建。混采數(shù)據(jù)在分離和重建的過程中對精度要求較高，這就需要采用高精度分離方法和準(zhǔn)確的的重建方法。

在炮分離過程中，在偽分離后的共檢波點域使用多級中值濾波法，相比于傳統(tǒng)中值濾波法有更高的分離精度，能保留更多的有效信息，為分離后的重建提供條件。但若濾波窗口選擇過大，則會丟失一些有效信息，影響濾波結(jié)果。國外學(xué)者較早對地震混合采集技術(shù)進行了研究。Sliverman et al.[5]提出同時激發(fā)震源的新思路；Beasley et al.[6-8]使用脈沖型震源進行研究；Bagaini et al.[9]對各種方法進行了橫向比較，并將其分類，多源地震混合采集技術(shù)(blended acquisition)便是由此成型，國內(nèi)劉財將多級中值濾波技術(shù)應(yīng)用到地球物理領(lǐng)域[10]，Huo et al.[11]在混采數(shù)據(jù)的分離工作中使用了中值濾波方法。

針對分離后的缺失地震數(shù)據(jù)，目前有多種重建方法，基于壓縮感知的地震數(shù)據(jù)重建方法和基于Radon變換重建方法是應(yīng)用較廣泛的方法,其中壓縮感知方法中的稀疏變換、迭代算法和閾值模型等的選取將影響最終地震數(shù)據(jù)重建的效果和計算效率[12]。Radon變換在對分離的單炮數(shù)據(jù)進行重建后，會進一步消除數(shù)據(jù)的噪聲干擾，基于此，本文采取Radon變換的方法來重建炮分離后的數(shù)據(jù)。

Hampmson[13]提出拋物Radon變換，Kabiret al.[14]采用拋物Radon變換重建地震數(shù)據(jù)取得了良好效果。Sacchiet al.[15]嘗試用高分辨率Radon變換重建地震道，在犧牲計算效率的情況下提高分辨率。國內(nèi)王維紅等[16]提出采用加權(quán)拋物線Radon變換法重建地震數(shù)據(jù)，一定程度兼顧了重建效率和效果。

在所有的Radon變換方法中，稀疏約束的雙曲Radon變換分辨率最高，采用稀疏約束控制改變正則化參數(shù)λ，來平衡反演誤差與模型稀疏化程度，能夠在數(shù)據(jù)重建中獲得更高的精度[17]，符合對混采數(shù)據(jù)重建精度的追求。然而雙曲Radon變換重建地震數(shù)據(jù)面臨的主要問題是雙曲Radon變換算子只能在時間域表示，一般用共軛梯度法求解，計算量較大，求解效率低[18-20]。為此，筆者引入了Fista[21]算法來替換共軛梯度法，顯著地提高了計算效率。模型試算結(jié)果表明，針對混采數(shù)據(jù)，本文方法在保證計算效率的情況下，取得了很好的重建效果。

1 方法原理

1.1 多級中值濾波的原理

在混采地震數(shù)據(jù)重建中，第一步進行炮分離，中值濾波作為處理工具，經(jīng)過多年的發(fā)展，此技術(shù)已經(jīng)趨于成熟。中值濾波的原理是設(shè)一個濾波長度N(通常為奇數(shù))，以選中的樣點為中心取出N個樣點，將取出的N個按從小到大順序依次排列，重排序列中心位置的數(shù)值即為輸出值。重復(fù)此過程即可實現(xiàn)中值濾波[22]。

而二維中值濾波原理為：

(1)

式(1)中，W1(n1,n2)表示以a(n1,n2)為中心的數(shù)值序列:{a(n1-N,n2),…,a(n1,n2),…,a(n1+N,n2)}，同理可以得出W2(n1,n2)，W3(n1,n2)，W4(n1,n2)所表示的數(shù)值序列。但中值濾波往往存在一些問題，當(dāng)信號比濾波長度小時，這些較小的細節(jié)信號會丟失，導(dǎo)致圖像中很多包含重要信息的細節(jié)結(jié)構(gòu)受到破壞，這種破壞對數(shù)據(jù)的影響比噪聲本身更為嚴(yán)重。

而多級中值濾波器可以在濾波過程中保護細節(jié)信息。其輸出定義為:

Ymin(n1,n2)=median[Ymax(n1,n2),Ymin(n1,n2),a(n1,n2)]

(2)

式(2)中，

(3)

(4)

Zi(n1,n2)=median[a(,)∈Wi(n1,n2)],

i=1,2,3,4

(5)

式中，median表示一般中級濾波[23]。

多級中值濾波具有自適應(yīng)特性，可以保護細節(jié)，稱之為“自適應(yīng)程度”，但這種“自適應(yīng)程度”具有局限性，僅限于Wi(i=1,2,3,4)中選擇兩個具有極大(極小)中值的基本子窗口，基本子窗口的中值與離散二維信號a共同確定輸出Ymin。濾波長度對去噪結(jié)果影響較大，降低濾波長度能夠保護有效信號，但去噪效果較差；提高濾波長度去噪效果較好，但會破壞有效信息[24]。

1.2 雙曲Radon變換原理

混采數(shù)據(jù)重建的第二步是進行地震道重建，對分離后的單炮數(shù)據(jù)進行Radon變換，經(jīng)過Radon正反變換后，缺失的的地震到得到恢復(fù)。Radon變換重建原理是基于水平層狀介質(zhì)模型反射波走時方程，走時方程可由級數(shù)展開，當(dāng)模型表現(xiàn)為各向同性時可以用簡單的雙曲線公式來描述縱波的同相軸，常規(guī)雙曲Radon[25]正變換為：

(6)

反變換為

(7)

式(7)中，τ是截距時間，p是慢度(速度倒數(shù))，x是炮檢距，d(t,x)表示道集數(shù)據(jù)，m(τ,p)是Radon變換域內(nèi)數(shù)據(jù)。在具體實現(xiàn)時，采用離散矩陣代替函數(shù)進行程序?qū)崿F(xiàn)，通常在最小二乘約束下的條件下求取正變換，矩陣形式的Radon反變換公式為:

d=Lm

(8)

式(8)中，d和m表示離散原始數(shù)據(jù)與Radon域數(shù)據(jù)的矩陣形式，定義L和LT分別為Radon變換算子和伴隨變換算子。建立最小平方意義下的線性化反演誤差目標(biāo)函數(shù)，正變換的最小平方解為[26]：

m=[(LTL)-1LTd]

(9)

1.3 稀疏約束反演

時間域內(nèi)稀疏約束的方法可以提高Radon域內(nèi)的分辨率，將Radon域數(shù)據(jù)作為稀疏條件進行約束反演，其中對稀疏約束模型取l1范數(shù)，對反演誤差取l2范數(shù)，變?yōu)榫€性反演問題的l1-l2混合范數(shù)求解，建立的目標(biāo)函數(shù)為：

(10)

式(10)中，λ是正則化參數(shù)，是平衡模型稀疏化與反演誤差的折中參數(shù)，λ的值越大，Radon域越稀疏，選取合適的λ值可以提高重建準(zhǔn)確度[27]。

1.4 參數(shù)選取

在重建時為避免產(chǎn)生假頻[28]，需要選擇合適的參數(shù)，為得到高分辨率的采樣結(jié)果，給出參數(shù)的選取準(zhǔn)則，參數(shù)的采樣率為

Δp=pk-pk-1

(11)

(12)

令Δv=vk-1-vk，得

(13)

當(dāng)vk和vk-1比較接近時，vk≈vk-1，式(13)變?yōu)?/p>

(14)

設(shè)原始數(shù)據(jù)中最大有效反射速度為vmax,要使式(14)對所有的有效波場均成立,則式中vk應(yīng)替換為vmax;同時,當(dāng)在CSP或CMP道集上做變換時, 選取的范圍通常為50～100, 得到臨界值關(guān)系為：

(15)

1.5 FISTA算法原理

時間域求解Radon變換時往往面臨計算量大的問題，其中，時間域算子L和LT是大型的稀疏矩陣，每次迭代都需要計算大型稀疏矩陣L和LT與向量的乘法[26]。求解‖d-Lm‖2時通常采用共軛梯度法，對欠定題給出最小范數(shù)解，對超定問題給出最小平方解，因此共軛梯度算法的遞推步驟較多，需要求解N*N階線性方程組，在實際運算中收斂速度也較慢。

快速迭代軟閾值算法(FISTA)可以有效地提高計算效率，加快反演的收斂速度，其迭代更新公式為：

m0=[(LTL)-1LTd]

(16a)

x0=m0

(16b)

t0=1

(16c)

當(dāng)k≥0時：

(17a)

(17b)

(17c)

式中，k是當(dāng)前迭代次數(shù)，soft是取軟閾值算子，α是Lipschtiz constant常數(shù)，α≥max(eig(LTL))。僅經(jīng)過數(shù)次迭代即可獲得高分辨率的Radon變換結(jié)果[30]。

2 模型試算及實際資料處理

2.1 理論模型

為檢驗本文方法的效果，模擬了一個混采模型數(shù)據(jù)，該地震記錄共100道，道間距50 m，每道500個采樣點，采樣間隔4 ms，具體模型如圖1a所示。圖1b為模擬該模型缺失地震距數(shù)據(jù)的情況，使原始數(shù)據(jù)第40至50道共10道數(shù)據(jù)缺失，第1道和第20道間隔采樣。普通中值濾波法進行炮分離得到圖1c，采用多級中值濾波法進行炮分離得到圖1d。由圖像可見，多級中值濾波很好地完成了炮分離的任務(wù)。

對多級中值濾波后分離圖1d，使用高分辨率雙曲Radon變換進行數(shù)據(jù)重建，用Fista算法迭代50次，得到高分辨率雙曲變換Radon域內(nèi)結(jié)果圖2a，其重建后的結(jié)果為圖2b。比較圖1d和圖2b可以看出，重建結(jié)果十分接近原始數(shù)據(jù),并且可以將分離殘留部分進一步壓制。為進一步驗證重建效果，分別對缺失地震道數(shù)據(jù)圖1d和雙曲Radon變換重建后數(shù)據(jù)圖2b做譜得到圖2c和2d。對比二譜發(fā)現(xiàn)，不規(guī)則缺失數(shù)據(jù)出現(xiàn)嚴(yán)重假頻，雙曲Radon變換重建結(jié)果假頻明顯消除。說明雙曲Radon變換有很好的重建效果，能有效地消除假頻。

在求解雙曲Radon變換時，分別試算了共軛梯度算法和Fista算法，本文的模型測試是在Intel(R)Core(TM)i5-2300 CPU@2.80GHz，16GB DDR3 1333MHz內(nèi)存的環(huán)境下運行。比較兩種求解方法發(fā)現(xiàn)，在迭代50次的情況下Fista算法耗時123.8 s。而共軛梯度算法迭代150次則耗時324.4 s?？梢奆ista算法計算效率遠高于共軛梯度法。

a.高分辨率雙曲變換Radon域內(nèi)結(jié)果； b.高分辨率雙曲Radon變換重建結(jié)果； c.缺失地震道數(shù)據(jù)譜； d.雙曲Radon變換重建數(shù)據(jù)譜.圖2 Radon變換重建效果對比Fig.2 Contrast of Radon transform results

2.2 實際數(shù)據(jù)處理

圖3a為中國某地區(qū)實際地震資料的CMP道集，該數(shù)據(jù)，每道626個采樣點，采樣間隔4ms,采樣時長2.5 s，共126道，道間距24 m。其中有效反射波的速度范圍在1 000～2 000 m/s之間。圖3b是對實際數(shù)據(jù)前30到35道缺失后的數(shù)據(jù)。對實際數(shù)據(jù)進行普通中值濾波法炮分離得到圖3c，采用多級中值濾波法進行炮分離得到圖3d。對比圖3c和3d可見，多級中值濾波炮分離效果明顯更好。

a.原始數(shù)據(jù);b.不規(guī)則缺失地震道數(shù)據(jù); c.普通濾波炮分離后數(shù)據(jù); d.多級中值濾波炮分離后的數(shù)據(jù).圖3 實際數(shù)據(jù)炮分離后對比Fig.3 Contrast of separate field blended data

經(jīng)過雙曲Radon變換后Radon域結(jié)果如圖4a所示，對比重建結(jié)果圖4b,可以看出實際數(shù)據(jù)重建效果較好，缺失數(shù)據(jù)道的位置得到了有效的恢復(fù)。實際數(shù)據(jù)證明了本文方法有較高的實用性。

3 結(jié)論

(1)針對混采的地震數(shù)據(jù)，多級中值濾波與普通濾波法相比，能更好地完成炮分離的任務(wù)，確保數(shù)據(jù)的重建精度。

a.雙曲Radon正變換后Radon域內(nèi)結(jié)果；b.雙曲Radon變換后重建結(jié)果.圖4 實際數(shù)據(jù)Radon變換重建結(jié)果Fig.4 Radon transform reconstruction results of field data

(2)針對炮分離后的數(shù)據(jù)，雙曲Radon變換可以很好地恢復(fù)缺失地震數(shù)據(jù)，得到準(zhǔn)確的地震數(shù)據(jù)重建結(jié)果。

(3)Fista算法在時間域求解的過程中，優(yōu)化了算法的迭代流程，和常規(guī)的共軛梯度法相比，明顯提高了計算效率，一定程度上解決了雙曲Radon變換計算量大的問題。模型數(shù)據(jù)和實際數(shù)據(jù)的重建結(jié)果表明，本文方法重建效果好，計算效率高。