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

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

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

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

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

按正問題的不同,透射波層析則可分為射線層析、波束層析和波動方程層析。射線層析基于高頻射線理論,對初始模型精度的要求較低,但反演精度不高[1]。由于正問題基于射線理論,受制于焦散、多路徑和陰影區(qū)等射線理論的固有缺陷,且層析矩陣十分稀疏,病態(tài)性強(qiáng),求解時收斂較慢,反演結(jié)果受射線照明的影響嚴(yán)重[2]。波動方程層析反演方法中,由于采用波動方程作為正演算子,因而可以同時預(yù)測地震數(shù)據(jù)的運(yùn)動學(xué)和動力學(xué)信息。經(jīng)典的數(shù)據(jù)域波動方程層析方法為全波形反演方法[3-4]。但是,正是因?yàn)椴捎貌▌臃匠逃嬎阕硬ǖ恼齻鞑ゲ▓雠c殘差的反傳播波場,即兩次正演計算,所以計算量較大,加之海量的采集數(shù)據(jù)、復(fù)雜地區(qū)實(shí)際資料的信噪比低導(dǎo)致波形匹配困難等原因使得波動方程層析基本處于理論研究階段,很少能夠得到實(shí)際應(yīng)用。

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

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

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

基于地表觀測的地震波數(shù)據(jù)估計地下介質(zhì)的彈性參數(shù)進(jìn)而進(jìn)行精確的油藏描述是勘探地震的核心問題。這是一個病態(tài)的、強(qiáng)非線性的問題,求解這樣的問題需要有充分體現(xiàn)其本質(zhì)的高度抽象的理論框架。地震數(shù)據(jù)是隨機(jī)的,確定性的反(散)射波場(同相軸)分布在滿足不同統(tǒng)計特征的隨機(jī)噪聲中?；诓煌昝赖牡乇碛^測數(shù)據(jù)和不完美的地震波預(yù)測模型進(jìn)行Bayes框架下的地震波反演成像就成了必然的理論選擇。Bayes估計理論框架下透射波反演成像方法的理論基礎(chǔ)描述如下。

1.1 透射波層析反問題數(shù)學(xué)表達(dá)與求解

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

(1)

(1)式為經(jīng)典數(shù)據(jù)域FWI的目標(biāo)泛函。式中:K(m)是依賴模型參數(shù)m的非線性函數(shù)(正問題),是模型空間M向數(shù)據(jù)空間D映射的有界非線性算子;CD表示與數(shù)據(jù)相關(guān)的協(xié)方差矩陣;CM則是與模型相關(guān)的協(xié)方差矩陣;mprior為先驗(yàn)?zāi)Ｐ托畔?dobs為觀測數(shù)據(jù)。

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

(2)

(3)

(4)

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

1.2 透射波數(shù)據(jù)殘差與預(yù)條件

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

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

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

建立透射波旅行時擾動與速度擾動之間的關(guān)系是透射波速度層析反演的最關(guān)鍵環(huán)節(jié)。

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

(5)

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

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

(6)

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

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

(7)

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

(8)

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

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

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

(11)

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

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

(12)

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

(13)

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

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

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

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

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

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

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

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

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

2.1 相位延遲

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

(16)

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

(17)

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

(18)

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

2.2 群延遲及瞬時走時

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

(19)

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

(20)

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

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

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

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

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

2.3 基于瞬時走時拾取的迭代優(yōu)化策略

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

(24)

(25)

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

(26)

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

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

2.4 數(shù)值試驗(yàn)

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

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

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

3 透射波旅行時Beam層析成像技術(shù)

3.1 Beam波路徑模擬

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

圖4 梯度速度場

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

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

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

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

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

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

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

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

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

圖8 過高速異常體縱向速度抽道對比結(jié)果a 左側(cè)異常體; b 右側(cè)異常體

3.3 觀測系統(tǒng)與透射波層析成像精度的關(guān)系分析

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

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

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

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

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

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

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

圖13 觀測系統(tǒng)示意(a)和真實(shí)速度場(b)

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

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

4 實(shí)際資料測試

4.1 瞬時旅行時拾取結(jié)果

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

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

圖16 采用傳統(tǒng)方法(綠線)和本文迭代優(yōu)化策略(紅線)進(jìn)行瞬時旅行時拾取的結(jié)果

4.2 反演結(jié)果

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

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

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

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

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

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

該工區(qū)數(shù)據(jù)生成層析敏感度矩陣,即使是稀疏存儲,仍然會占用計算機(jī)約378GB內(nèi)存。顯而易見,直接采用傳統(tǒng)LSQR算法求解,需要單節(jié)點(diǎn)至少擁有400GB內(nèi)存以上,否則難以實(shí)現(xiàn)。因此,本文采用MPI+OpenMP策略,即以MPI模式將靈敏度矩陣存儲分?jǐn)偟蕉鄠€節(jié)點(diǎn)上,單節(jié)點(diǎn)則采用OpenMP模式加速計算,以此降低對單節(jié)點(diǎn)內(nèi)存的要求。我們使用了7個節(jié)點(diǎn),每個節(jié)點(diǎn)上分?jǐn)偞蠹s54GB內(nèi)存。采用共軛梯度法求解,共耗時約17.5h,主要是數(shù)據(jù)IO以及多節(jié)點(diǎn)歸約等待。

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

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

5 結(jié)論

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

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

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

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