基于雙曲Radon變換的混采地震數據重建

王晗，韓立國

吉林大學 地球探測科學與技術學院，長春 130026

0 引言

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

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

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

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

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

1 方法原理

1.1 多級中值濾波的原理

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

而二維中值濾波原理為：

(1)

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

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

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]。

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

1.2 雙曲Radon變換原理

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

(6)

反變換為

(7)

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

d=Lm

(8)

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

m=[(LTL)-1LTd]

(9)

1.3 稀疏約束反演

時間域內稀疏約束的方法可以提高Radon域內的分辨率，將Radon域數據作為稀疏條件進行約束反演，其中對稀疏約束模型取l1范數，對反演誤差取l2范數，變為線性反演問題的l1-l2混合范數求解，建立的目標函數為：

(10)

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

1.4 參數選取

在重建時為避免產生假頻[28]，需要選擇合適的參數，為得到高分辨率的采樣結果，給出參數的選取準則，參數的采樣率為

Δp=pk-pk-1

(11)

(12)

令Δv=vk-1-vk，得

(13)

當vk和vk-1比較接近時，vk≈vk-1，式(13)變為

(14)

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

(15)

1.5 FISTA算法原理

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

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

m0=[(LTL)-1LTd]

(16a)

x0=m0

(16b)

t0=1

(16c)

當k≥0時：

(17a)

(17b)

(17c)

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

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

2.1 理論模型

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

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

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

a.高分辨率雙曲變換Radon域內結果； b.高分辨率雙曲Radon變換重建結果； c.缺失地震道數據譜； d.雙曲Radon變換重建數據譜.圖2 Radon變換重建效果對比Fig.2 Contrast of Radon transform results

2.2 實際數據處理

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

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

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

3 結論

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

a.雙曲Radon正變換后Radon域內結果；b.雙曲Radon變換后重建結果.圖4 實際數據Radon變換重建結果Fig.4 Radon transform reconstruction results of field data

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

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