透射波旅行時Beam層析成像方法

羅 飛,王華忠,馮 波,彭仁艷,蔣先勇

(1.波現象與智能反演成像研究組(WPI),同濟大學海洋與地球科學學院,上海200092;2.中石化石油工程地球物理公司國際業務發展中心,江蘇南京210000;3.中國石油化工股份有限公司中原油田分公司物探研究院,河南濮陽457000)

在地震數據信噪比低的情況下,建立準確的背景速度很困難,而透射波信噪比高,受道間時差的影響小,只要觀測方式合適,利用透射波旅行時層析可以建立比較準確的表層(0～200m)、淺層(200～500m)甚至中層(500～1500m)的速度模型。對于面向高分辨和高保真的精確成像,即使淺層速度只存在較小的相對誤差,也會導致中深層成像達不到同相位疊加,從而嚴重影響高分辨率和高保真成像。

如今,地震數據采集技術進步明顯,“兩寬一高”(寬頻帶、寬方位、高密度)地震數據采集逐漸成為主流,研究與采集技術相匹配,能夠充分挖掘地下介質彈性參數信息的地震數據分析方法越來越重要。特別是,在“兩寬一高”地震數據采集中,偏移距的增大使透射波信息更加豐富,除了來自淺層的透射波,中深層的透射波信息也能被采集到,從而使透射波具備了反演中深層彈性參數的能力。寬方位采集方式還能優化透射波的照明范圍,降低了反問題的病態性。而高密度的地震采集則提供了更多的數據,配合適當的數據預條件方法,可以很好地降低層析核函數的稀疏性,使層析反演結果更加穩定。同時,相比反射波層析成像算法,透射波層析更容易實現自動化,對工業生產更具現實意義。因此非常有必要發展一套高精度的透射波層析反演方法與技術流程。

利用透射波信息進行層析成像,屬于數據域反演方法,根據利用的數據屬性,透射波層析成像技術分為走時層析、相位層析和波形層析。與波形層析相比,走時擾動對速度擾動的非線性程度較弱,對初始模型的要求較低,走時層析在實際生產中得到廣泛應用。

按正問題的不同,透射波層析則可分為射線層析、波束層析和波動方程層析。射線層析基于高頻射線理論,對初始模型精度的要求較低,但反演精度不高[1]。由于正問題基于射線理論,受制于焦散、多路徑和陰影區等射線理論的固有缺陷,且層析矩陣十分稀疏,病態性強,求解時收斂較慢,反演結果受射線照明的影響嚴重[2]。波動方程層析反演方法中,由于采用波動方程作為正演算子,因而可以同時預測地震數據的運動學和動力學信息。經典的數據域波動方程層析方法為全波形反演方法[3-4]。但是,正是因為采用波動方程計算子波的正傳播波場與殘差的反傳播波場,即兩次正演計算,所以計算量較大,加之海量的采集數據、復雜地區實際資料的信噪比低導致波形匹配困難等原因使得波動方程層析基本處于理論研究階段,很少能夠得到實際應用。

可見,射線層析與波動層析是層析反演的兩個極端,前者效率高但其效果受射線相關假設的制約,后者理論上完善但效率及初始模型限制了方法的實際應用。射線層析明顯與波傳播規律不符,可以從直觀上改進射線層析,引入波束的概念[5-8],采用Beam束模擬波傳播路徑,將走時殘差投影到Beam束路徑上,這樣做在不增加太多走時層析計算量的基礎上可以大幅度降低射線“盲區”數量及范圍。更重要的是,Beam束層析擺脫了高頻近似的射線理論,考慮了波的頻率信息,使得走時層析更加穩定,分辨率與精度也得到了顯著提高。因此,本文主要論述如何利用Beam束層析成像方法從透射波旅行時中挖掘與彈性參數相關的信息。

本文從層析成像的理論框架出發,系統論述了波動方程層析、波束層析以及射線層析的異同,剖析了參數迭代公式中各項的含義,總結出利用透射波旅行時層析成像方法挖掘“兩寬一高”地震數據信息時應注意的若干核心問題(層析正問題、數據的預條件與模型的正則化),并重點闡述了加權走時的含義和自動化的獲取加權地震走時信息的方法。最后,論述了Beam層析算法理論以及對觀測系統的適應性,結合MPI+OpenMP策略在某探區實際資料上進行了應用,其處理結果說明了透射波旅行時Beam層析處理流程在“兩寬一高”海量地震數據處理中的可行性。

1 透射波層析成像的理論框架

基于地表觀測的地震波數據估計地下介質的彈性參數進而進行精確的油藏描述是勘探地震的核心問題。這是一個病態的、強非線性的問題,求解這樣的問題需要有充分體現其本質的高度抽象的理論框架。地震數據是隨機的,確定性的反(散)射波場(同相軸)分布在滿足不同統計特征的隨機噪聲中。基于不完美的地表觀測數據和不完美的地震波預測模型進行Bayes框架下的地震波反演成像就成了必然的理論選擇。Bayes估計理論框架下透射波反演成像方法的理論基礎描述如下。

1.1 透射波層析反問題數學表達與求解

基于Bayes反演理論[9-10],在L2范數下定義誤差泛函為:

(1)

(1)式為經典數據域FWI的目標泛函。式中:K(m)是依賴模型參數m的非線性函數(正問題),是模型空間M向數據空間D映射的有界非線性算子;CD表示與數據相關的協方差矩陣;CM則是與模型相關的協方差矩陣;mprior為先驗模型信息;dobs為觀測數據。

通常而言,上述泛函采用局部最優化方法求解。若給定初始模型mk-1,牛頓迭代算法寫為:

(2)

(3)

(4)

式中:μk-1為最速下降迭代步長。

1.2 透射波數據殘差與預條件

如公式(4)所述,我們根據不同的需求,靈活選擇不同的透射波屬性(δd)進行反傳,屬性可以是走時、相位、振幅或波形等。理論上,為了充分挖掘透射波攜帶的信息,應該從透射波走時反演出發,逐漸過渡到波形反演。結合前文分析,本文重點論述如何充分挖掘透射波走時(δt)中攜帶的地下彈性參數信息。

1.3 透射波旅行時層析正問題(K(m))

為了求解公式(4),必須研究透射波走時層析的正問題t=K(m)。最經典的透射波層析方法是波形層析方法,譬如初至波或早至波波形層析。低頻長偏移數據下的FWI本質上實現了波形層析,疊前數據Laplace變換下的FWI是一種典型的做法[13-14],通過調節Laplace變換中的系數控制參與層析反演的早至波成分的多少。如此進行透射波波形層析的優勢是基本不需要人工干預,但初至波及后續一段時間內的波現象非常復雜,成因不明,加上初始模型不準,數據中缺乏低頻成分,透射波波形層析很難收斂。因此,針對陸上數據,最好是進行波動理論的初至波旅行時層析,如果早至波到達時檢測問題解決得好,還可以進行波動理論的早至波旅行時層析。

建立透射波旅行時擾動與速度擾動之間的關系是透射波速度層析反演的最關鍵環節。

通過波場的Rytov近似可以導出旅行時擾動與速度擾動之間的線性關系,得到層析反演的線性矩陣方程,然后在最小二乘意義下進行層析反演。馮波基于波動方程的Rytov近似并采用模擬數據的歸一化能譜作為單頻走時擾動的加權函數,導出帶限地震波的走時敏感度核函數[15-16]:

(5)

從(5)式可以看出,對于任意一對炮檢(xr,xs),均需計算一次波場正演(計算u0(x,t;xs))以及一次波場逆時傳播(計算p(x,t;xr),因此計算量遠高于射線追蹤或波束層析。

相應地,地震波走時層析正問題矩陣形式可表示為:

(6)

式中:δt為透射波的走時殘差;δv為速度參數擾動。

從程函方程的角度也可以建立速度擾動和旅行時擾動之間的關系進行射線理論的透射波旅行時層析反演。為此,首先給出波動方程的平面波解:

(7)

式中:p(xi,t)為空間某一點xi在t時刻的壓力場;P(xi)為振幅項;T(xi)為旅行時;ω=2πf為圓頻。(7)式也是波動方程WKBJ近似的零階解。將(7)式回代到聲波波動方程中,得到:

(8)

當圓頻ω→∞(高頻近似)時,(8)式中第三項可以忽略,同時為了使公式(8)成立,需要滿足如下兩個方程:

公式(9)和公式(10)便是波動方程在高頻近似下推導出來的程函方程和輸運方程。公式(9)是關于旅行時T(xi)的一階非線性偏微分方程,是射線旅行時理論的基礎;公式(10)則為振幅項P(xi)的一階線性偏微分方程。

求解程函方程可以進行射線追蹤,建立起如下的沿射線路徑的速度和旅行時之間的積分關系:

(11)

式中:s表示慢度場模型;s(xj)表示在點xj處的慢度值,與速度場的關系為s(xj)=v-1(xj);dli(xj,s)表示在慢度模型s中第i根射線在網格xj處的射線長度;Ri(s)代表在模型s中第i根射線的射線路徑;ti是第i根射線旅行時。(11)式屬于Fredholm第一類積分方程(Radon變換也屬于此類方程)。

在此基礎上,引入介質慢度擾動量δs,并對旅行時進行了Taylor展開并只保留一階項,假設介質的小擾動不會引起射線路徑太大的變化,即假設射線在光滑介質中穿過,可以進一步得到慢度擾動量和旅行時擾動量之間的近似線性關系:

(12)

公式(12)引入射線層析核函數矩陣KRay,則射線層析正問題的矩陣方程形式可寫為:

(13)

基于公式(13),當引入波束的概念時,傳統的射線類走時層析核函數矩陣KRay和正問題的矩陣形式可重寫為:

式中:KBeam為引入波束概念后的核函數,以射線交點為中心呈高斯衰減,表示Beam-ray計算的波場振幅可以作為加權因子,來提高透射波層析的反演精度;(l,q)是射線中心坐標系,l是沿射線的弧長,q代表垂直于射線的距離;σ是標準差,與頻率相關。這里需要指出的是,根據對KBeam的不同定義,在保證沿垂直射線的橢圓(三維)或線段(二維)的積分為1的前提下,可以形成不同形式的波束(Beam)。

結合公式(6)、公式(13)、公式(15)以及公式(4)可以看出,從波動方程出發,通過不同的假設近似,透射波旅行時層析反演從波動理論逐漸退化到射線理論,復雜程度逐步下降,挖掘透射波走時能力也在遞減,其計算耗時也呈指數形式下降。因此,折中選取研究波束走時層析,在精度和效率兩者間做個平衡,更具有現實意義。

1.4 透射波層析模型正則化(CM)

地震波反演問題是一個典型的病態問題,但求解這一病態問題具有現實意義。求解病態反問題的基本思想是施加各種正則化方法,使得反問題的解能解決實際問題。

從目前的實際情況來看,促進層析成像技術走向實用化的關鍵之一就是反演過程中模型正則化的施加。本質上,模型正則化是利用模型的先驗信息約束模型估計本身及其分布范圍,降低模型估計的多解性。通過模型參數協方差(逆)算子對反演解進行正則化約束是非常根本的方法。MANDELLI等[17]系統對比了不同模型正則化(傾角濾波、TV以及各向異性擴散)的差異,認為把地質構造約束定量化地引入到層析成像中,將增加反演的穩定性。KAZEI等[18]則從另一個角度出發,將模型正則化轉化為懲罰函數選擇問題,從模型的L2范數問題過渡到Lp范數,從而實現由低波數向高波數的過渡。ROMAHN等[19]利用已知井信息構建模型正則化算子CM,相當于構建了一個匹配濾波器,作用于反演速度擾動上,達到用井信息修正梯度的目的。

2 透射波層析成像中到達時測量問題

地震波走時的拾取是透射波層析成像中重要的一步。一般意義下,地震波走時拾取表達一個物理實體的到達,數學上抽象為一個質點的到達。雖然數學上一個質點的到達概念很明確,但是不具可測性,因此,實體的到達時必然引起不確定性。同理,子波的到達時測量也存在不確定性。理論意義下,在時間域中,測理論脈沖的到達時是準確的,測帶限波形的到達時存在不確定性;在頻率域中,用單頻波的相位同相軸測時差是可行的,而測單頻波的到達時存在困難。傳統的旅行時拾取算法大體可分為兩類:滑動時窗方法和相干類方法。在滑動時窗方法[20-21]中,地震信號序列的屬性在連續或重疊的移動窗口中計算。相干類方法[22]依賴于使用一些相似度測量技術,比如利用互相關來比較單個或多個波形,從而得到走時信息。

本文在相速度和群速度的基礎上,定義了基于相位延遲和群延遲[23]的兩種加權帶限子波到達時。并在文獻[24]和文獻[25]的基礎上,討論了基于群延遲的自動走時拾取的迭代優化策略。通過引入一個濾波算子,迭代拾取的策略能夠自動拾取走時,用于透射波層析反演算法。

2.1 相位延遲

在真空中,所有波長的電磁波均以相同的相速度——真空中的光速傳播。在色散介質中,由于介質的折射率與光的波長(頻率)有關,不同波長(頻率)的電磁波具有不同的相速度。程乾生[23]給出了相應“相位延遲”的概念,一般的濾波器頻譜S(ω)為:

(16)

式中:A(ω)=|S(ω)|為振幅譜;Φ(ω)為相位譜。濾波器的相位延遲Tp(ω)為:

(17)

相位延遲表示頻率為ω的單一正弦波經過濾波器后的延遲時間。如果將檢波器接收到的信號看成若干頻率的諧波經過地球濾波后被檢波器接收,則帶限子波走時可以看作相位延遲的一個加權效應。XIE等[26]指出帶限地震信號的走時擾動可以用單頻相位延遲擾動加權疊加表示:

(18)

式中:P(ω)為加權算子。馮波[15]詳細論述了公式(18)的物理含義,在此不再贅述。

2.2 群延遲及瞬時走時

實際存在的波不是單頻的,傳播介質對這個(或這些)波必然是色散的,那么,傳播中的波由于各不同頻率的成分運動快慢不一致,會出現擴散,但假若(注意這個假設)這個波由一系列頻率差別不大的簡諧波組成,這時在相當長的傳播過程中整個波組仍將維持為一個整體,以固定的速度傳播,這個特殊的波群稱為“波包”,這個速度稱為群速度,也稱包絡延遲。與相速度不同,群速度比波包的中心相速度小,并且二者的差值同中心相速度隨波長而變化的平均率成正比。群速度是波包的能量傳播速度,也是波包所表達信號的傳播速度。濾波器的群延遲Tg(ω)為:

(19)

群延遲(包絡延遲)反映某一頻率的包絡經過濾波器后的延遲時間,體現了在頻率ω的鄰域內局部延遲性質。如果我們將檢波器接收到的信號看成若干頻率的地震波經過地球濾波后被檢波器接收,則帶限子波走時可以看作群延遲在一定頻率內的平均效應(平均時間),可寫為:

(20)

LUO等[27-28]闡明了瞬時旅行時的具體意義,其等價于群延遲。在時-頻域內瞬時走時τ(t)的表達形式為:

式中:S(t,ω)為信號s(t)的時頻變換;Im{·}代表取虛部。(22)式中,計算導數dS(t,ω)/dω可能不穩定,因此采用FOMEL[29]提出的平滑除法以減小除法帶來的不穩定。

根據公式(21),計算帶限信號的瞬時走時可以分為互不干涉的兩步:首先,使用公式(22)將單道地震數據映射到時間-頻率域;然后,通過算子Mω(公式(23))從時頻域重新映射回時間域。我們通過圖1和圖2 進一步說明如何計算瞬時走時。圖1a展示的是由兩個地震事件(0.32s,1.28s)合成的某一地震信號。圖1b和圖1c分別為合成信號的時頻譜S(t,ω)和采用公式(22)計算得到的時頻圖τ(t,ω)。最終計算的結果如圖2所示。其中圖2a表示將τ(t,ω)通過映射算子作用,也就是在有效帶寬內對頻率平均后的結果,即τ(t),當τ(t)=t時,其對應的走時t則為瞬時旅行時。而旅行時的自動拾取結果t位于0.32s和1.28s處,如圖2b中脈沖位置所示(紅線處)。

圖1 合成地震信號(t=0.32,1.28s)(a)、合成信號的時頻譜(b)和τ(t,ω)時頻譜(c)

圖2 τ(t,ω)經過有效帶寬內對頻率平均后的結果τ(t)(a)和旅行時拾取結果(紅線)(b)

2.3 基于瞬時走時拾取的迭代優化策略

這里,我們引入一個阻尼濾波算子:

(24)

(25)

η(i)={t|t∈Method(STA/LTA,etc.)}

(26)

式中:w(i,m)為定義的迭代濾波算子,該算子作用于每一道;a為阻尼因子,決定了濾波器的寬度;T0和M分別代表地震信號的采樣時間和最大迭代次數;η(i)相當于一個粗略的初至走時。

這種方法可以分為以下幾個步驟:①通過公式(26)計算每一道的η(i);②為每一次迭代設計濾波器算子w(i,m)并作用于每一道上;③利用公式(21)和公式(23)計算每一道的瞬時旅行時;④重復步驟②和步驟③,直到迭代結束。

2.4 數值試驗

我們使用YILMAZ[30]提供的Data6作為實際測試數據。該地震數據為炸藥震源激發,48道接收,道間距為100m,采樣間隔為4ms。該實際數據存在多個接近初至波信號的來自更深層的折射波信號,兩次反射波信號約1.5s和2.7s。

我們分別采用Saragiotis的方法和基于瞬時走時拾取的迭代優化策略對該實際數據進行處理,結果如圖3所示。對比圖3a與圖3b可以看出,改進后的方法能夠更加準確地拾取初至以及反射走時信息,在此基礎上,配合一定的質量控制(QC)準則,可以逐層選取合適的透射旅行時信息用于層析成像。

圖3 采用Saragiotis的方法(a)和基于瞬時走時拾取的迭代優化策略(b)得到的旅行時拾取結果

3 透射波旅行時Beam層析成像技術

3.1 Beam波路徑模擬

當引入波束的概念時,可以用公式(11)和(12)模擬Beam波路徑。圖4顯示了三維梯度速度場的Beam波路徑模擬情況。該模型X方向3750m,共301個網格采樣點,網格大小為12.50m;Y方向2500m,

圖4 梯度速度場

共201個網格采樣點,網格大小為12.50m;深度h=1875m,共301個網格采樣點,網格大小為6.25m。梯度速度模型可表示為:v(h)=600+(h/6.25)·3,單位m/s。

將炮點放置在(100m,800m,25m)處,檢波點放置在(1625m,800m,25m)處。圖5為不同頻率的三維Beam波路徑,可以看出,沿著半徑方向能量呈高斯衰減,隨著頻率的增加Beam波半徑變小,當頻率趨近于無窮大時,Beam束近似為射線。

3.2 三維透射波旅行時Beam層析與射線層析對比分析

為了具體說明Beam束相對于射線層析的優勢,采用如圖6a所示的速度模型進行數值試驗。該模型在近地表附近存在兩個高速異常體。該模型X方向3750m,共301個網格采樣點,網格大小為12.50m;Y方向2500m,共201個網格采樣點,網格大小為12.50m;深度為1875m,共301個網格采樣點,網格大小為6.25m。背景梯度速度場可表示為:v(h)=600+(h/6.25)·3,單位m/s,高速異常體速度為1600m/s。每條炮線上有11炮,總共16條炮線,共計176炮,炮線距為250m,炮間距為250m,道間距125m。在此基礎上,只選擇滿覆蓋左側異常體的4條線數據進行反演。

圖5 正演模擬的Beam波路徑a 沿X方向的切片(20Hz); b 沿Y方向的切片(20Hz); c 沿X方向的切片(30Hz); d 沿Y方向的切片(30Hz)

以圖6b所示的背景梯度場為初始速度場,經過5次反演迭代的結果如圖7所示。圖7a為傳統射線類三維透射波層析反演結果(未加入任何約束)。圖7b 為三維Beam束層析反演結果。對比圖7a和圖7b可見,三維Beam束層析結果分辨率更高。為了進一步說明效果,抽取分別穿過兩個異常體的縱向速度曲線進行對比,結果如圖8所示,可以看出,經過5次迭代的Beam層析(綠線)結果明顯優于經過5次迭代的Ray層析(紅線)結果。

相對射線走時層析而言,Beam束走時層析在不增加太多計算量的基礎上可以大幅度降低射線“盲區”數量;同時,該方法考慮了波的頻率信息,使得走時層析更加穩定,分辨率與精度也得到了顯著提高。在同等觀測系統條件(不完備)下,Beam束參數反演縱橫分辨率以及反演深度都比射線類方法要好。一般,Beam束層析耗時高于射線層析,但遠低于傳統波動方程層析,因此可以在生產上,特別是處理“兩寬一高”地震數據時發揮一定作用。

圖6 三維高速異常體模型(a)和初始速度場(b)

圖7 三維透射波旅行時Beam層析迭代5次反演結果a 射線層析反演結果; b Beam束層析反演結果

圖8 過高速異常體縱向速度抽道對比結果a 左側異常體; b 右側異常體

3.3 觀測系統與透射波層析成像精度的關系分析

隨著“兩寬一高”采集的出現,地震波中包含的信息越來越豐富。ZHOU等[31]通過對實際資料的處理評價了透射波層析成像技術的能力,如圖9所示,采用不同偏移距數據進行透射波層析反演得到的速度模型精度存在一定差異。可見,在充分利用中、深層折射波以及回轉波信息(長/超長偏移距生成的初至/早至波)時,透射波層析反演的精度能夠得到一定程度的提升。

首先,截取BP模型的部分數據,采用不同道間距數據進行反演。圖10顯示了真實速度模型以及初始速度模型。圖11為不同道間距的反演結果。圖12 為橫向速度的抽道對比結果。通過分析可知,對道間距的基本要求是:反演目標越接近地表、異常體尺度越小,道間距越小。同時,考慮到成本以及效率的原因,在保證較好反演分辨率和消除道間距時差的基礎上,可以適當選取較大的道間距進行反演。

接著,設計如圖13a所示的寬方位觀測系統(藍線為炮線,紅線為檢波線),5條炮線,炮線距200m,炮線長度為16km,每條炮線上有639炮,炮間距25m,因此存在長偏移距信息。37條檢波線,檢波線間距為200m,檢波線長度為4800m,檢波點間距為50m,相當于1炮激發,3552道接收。并以圖13b所示真實速度場正演數據,使用梯度速度場作為初始速度。通過對比不同形式地震數據的反傳梯度的區別,簡要總結了寬方位觀測系統對Beam透射波層析成像方法技術的影響。

圖9 真實速度場和炮記錄(a)、5km偏移距數據反演結果以及正演記錄(b)以及全偏移距數據(14km)反演結果和正演記錄(c)[31]

圖10 真實速度場(a)和初始速度場(b)

圖14為二炮線與五炮線數據速度反傳梯度的對比結果,可以看出,隨著炮線的增加,有效方位增加,從而射線覆蓋更好,某些在二炮線時未反演出來的區域在五炮線情況下出現。圖15為不同偏移距情況下的反傳梯度。可見,隨著偏移距的逐漸增加,深部照明能量也在不斷增加。

結合前面的數值分析結果可知,透射波層析反演的精度與觀測系統密切相關。在各向同性介質假設下,寬方位的采集模式能夠優化照明范圍,提升反演精度。同時,需要采集長偏移距/超長偏移距數據,這樣可以加大反演深度,得到更深層的反演結果。在此基礎上,道間距并非越小越好,而是應該根據處理數據的規模、反演目標的尺度/深度以及生產成本等信息,權衡分辨率與成本,在高密度數據中抽取合適的道間距。總之,合適的高密度寬方位長偏移距觀測系統,能夠提高透射波層析成像的精度以及深度。

圖13 觀測系統示意(a)和真實速度場(b)

圖14 二炮線(a)與五炮線(b)數據反傳梯度的對比結果

圖15 最大偏移距為9600m(a)和13500m(b)情況下的反傳梯度

4 實際資料測試

4.1 瞬時旅行時拾取結果

海上某“兩寬一高”工區的實際地震記錄如圖16所示,共639道,道間距25m,最大偏移距為13625m,縱向采樣4100個,采樣間隔為1ms。該數據在紅框處存在嚴重的吸收衰減現象,信噪比很低;在綠框處由于偏移距較大導致深層折射波與初至波相鄰很近,嚴重影響了拾取精度。分別采用傳統方法和本文迭代優化策略進行瞬時旅行時拾取,其結果分別如圖16中的綠線和紅線所示,圖17a和圖17b分別為圖16中紅框處和綠框處的放大顯示。對比圖17a 和17b可以發現,采用本文方法能夠得到更加精確的初至旅行時。

按第3節中介紹的方法,半自動拾取了該工區共101309153道的瞬時旅行時。通過QC處理保留大約85961691道,有效道數占比為84.85%。需要指出的,QC剔除的道,主要集中在吸收衰減區域以及超長偏移距處(圖16),這是因為采用半自動化(加入一些人為QC策略)算法時同樣的參數可能不適合所有的數據,所以拾取質量有差異。

圖16 采用傳統方法(綠線)和本文迭代優化策略(紅線)進行瞬時旅行時拾取的結果

4.2 反演結果

獲取合適的透射波旅行時信息后,我們從圖18a所示初始速度場出發,采用透射波旅行時Beam層析算法進行了兩次迭代,其迭代更新量如圖18b所示,圖18c為更新后的速度場,共處理大約1.6TB數據。

圖17 圖16中紅框處(a)、綠框處(b)的放大顯示

分析迭代更新量可以發現,在速度場的東北方向存在一個高速異常界面,與地震數據上出現遠偏移距吸收衰減嚴重以及生成折射波現象相吻合。圖19為走時差對比結果(紅色箭頭指示為拾取走時與更新后速度場正演走時作差的結果,而綠色箭頭則為拾取走時與初始速度場正演走時作差的結果)。可見,采用更新后的速度正演模擬走時的結果更趨近于拾取結果。分別對更新前、后的速度場進行偏移成像,圖20為采用雙平方根疊前深度偏移得到的X方向的偏移剖面切片,對比圖20a和圖20b的綠色橢圓處可以看出,采用更新后速度得到的偏移結果同相軸更加連續聚焦。結合迭代更新量的速度異常分布、走時差對比以及偏移結果,說明透射波旅行時Beam層析速度反演有能力處理大規模的“兩寬一高”地震數據。

圖18 某工區處理結果a 初始速度場; b 迭代更新量; c 更新后的速度場

圖19 炮檢互逆后某炮所有檢波點的走時差 圖中,紅色箭頭指示拾取走時與更新后速度場正演走時差的結果;綠色箭頭指示為拾取走時與初始速度場正演走時差的結果。

4.3 大型稀疏矩陣的求解性能分析

該工區數據生成層析敏感度矩陣,即使是稀疏存儲,仍然會占用計算機約378GB內存。顯而易見,直接采用傳統LSQR算法求解,需要單節點至少擁有400GB內存以上,否則難以實現。因此,本文采用MPI+OpenMP策略,即以MPI模式將靈敏度矩陣存儲分攤到多個節點上,單節點則采用OpenMP模式加速計算,以此降低對單節點內存的要求。我們使用了7個節點,每個節點上分攤大約54GB內存。采用共軛梯度法求解,共耗時約17.5h,主要是數據IO以及多節點歸約等待。

可見,面對“兩寬一高”采集數據,特別是TB級規模的數據,層析核函數的稀疏存儲與求解問題仍然是工業生產需要面對的難題。一方面,需要研究無需顯式計算和存儲核函數的求解層析反問題的算法,比如隱式矩陣向量乘方法、伴隨狀態法[32]等。這類算法將層析反問題轉化為單道并行處理,降低對計算機內存的要求。另一方面,“兩寬一高”的地震數據為我們提供了更加充足完備的樣本,我們可以借鑒機器學習中的優化算法,引入隨機梯度的概念,將地震數據分成若干樣本子集,每次只使用部分數據進行反演。總之,如何在現有計算機條件下,完成“兩寬一高”地震數據的大規模層析成像處理,是值得我們關注且具有現實意義的。

圖20 采用雙平方根疊前深度偏移得到的X方向的偏移剖面切片a 更新前速度偏移結果; b 更新后速度偏移結果

5 結論

復雜介質情況下,“兩寬一高”采集的地震數據,需要有針對性的參數估計方法。高密度采集加上炮集中透射波具有較高的信噪比,利用透射波(初至波)進行層析成像對于提高淺中層速度(包括Q值)建模的精度具有重要意義。

透射特征波場的有效提取是利用透射波信息進行層析反演的關鍵環節。初至帶限子波的到達時測量(計算瞬時走時),應該分為四步:首先,提取數據某種特征;然后,定義走時與該特征的關系;再次,度量觀測數據和模擬數據該種特征的距離;最后,對拾取結果進行質量監控,剔除精度較差的結果。同時,隨著地震數據量的日益增大,孤立地討論單道的地震事件提取而忽視數據的橫向空間連續性,這樣沒有充分利用數據的特征信息,從而降低特征提取的精度。因此,在機器學習框架下,將地震數據拓展到高維空間,基于多域、多屬性特征的自動地震旅行時拾取技術,或許是一個很好的思路。

相對射線走時層析而言,Beam束走時層析在不增加太多計算量的基礎上可以大幅度降低射線“盲區”數量,并且考慮了波的頻率信息,使得走時層析更加穩定,分辨率與精度也得到了顯著提高。數值試驗結果表明,相同情況下Beam束層析反演縱橫分辨率以及反演深度都比射線類方法好。經過海上實際資料測試,透射波(初至和早至波)Beam層析反演的深度并不局限于淺層,隨著采集偏移距的增加,其反演深度將隨之增加(達到淺層甚至中層)。采用寬方位數據進行透射波層析反演,增加了波場的照明范圍,使反演解更加穩定。高密度采集則為我們提供了充足的數據信息,允許我們靈活設計合適的求解策略。因此,和透射波層析反演算法相匹配的觀測系統應該是寬方位高密度長偏移距觀測系統。同時,Beam束層析耗時遠低于傳統波動方程層析,但仍高于射線層析。為了更好地在“兩寬一高”采集數據中利用透射波信息(特別是TB級數據),透射波旅行時Beam層析反演的精度與效率問題仍需繼續探究。

致謝:感謝中國石油天然氣股份有限公司勘探開發研究院、西北分院,中海石油(中國)有限公司研究中心、湛江分公司以及中國石油化工股份有限公司物探技術研究院、勝利油田分公司對波現象與智能反演成像研究組(WPI)研究工作的資助與支持。