二維VTI介質(zhì)qP波斜率層析方法

谷巍巍,楊 鍇,王 瀟

(同濟(jì)大學(xué)海洋與地球科學(xué)學(xué)院反射地震學(xué)科組,上海200092)

立體層析是將一個(gè)局部相干同相軸在炮道集與檢波點(diǎn)道集內(nèi)的射線參數(shù)水平分量(psx、prx,以下簡(jiǎn)稱地表觀測(cè)px參數(shù))納入到數(shù)據(jù)空間之中,使得數(shù)據(jù)空間的信息相對(duì)傳統(tǒng)反射層析更為豐富。除了旅行時(shí)t外,地表觀測(cè)px參數(shù)與炮檢點(diǎn)坐標(biāo)也都納入到立體層析數(shù)據(jù)空間中,使得立體層析成為層析成像方法中唯一一種可以同時(shí)反演速度、反射點(diǎn)位置、反射張角與構(gòu)造傾角的方法[1-4]。

斜率層析屬于立體層析的一種特殊情況,在利用運(yùn)動(dòng)學(xué)反偏移重建數(shù)據(jù)空間的過程中,會(huì)出現(xiàn)炮檢點(diǎn)位置和旅行時(shí)信息已經(jīng)匹配而地表px參數(shù)無法匹配的情況,這時(shí)可利用px殘差完成對(duì)速度模型、反射點(diǎn)坐標(biāo)和散射角的更新[5-6]。斜率層析算法也因此而得名。本文將上述情況下的各向同性介質(zhì)斜率層析方法推廣到VTI介質(zhì)。與前人研究工作的不同之處在于,當(dāng)VTI介質(zhì)中的運(yùn)動(dòng)學(xué)反偏移完成之后,我們利用px殘差僅反演三個(gè)各向異性參數(shù),而不反演反射點(diǎn)坐標(biāo)和散射角這些與地下成像點(diǎn)有關(guān)的局部運(yùn)動(dòng)學(xué)信息。這是由于各向異性參數(shù)更新之后,會(huì)立即執(zhí)行VTI介質(zhì)中的運(yùn)動(dòng)學(xué)射線追蹤,以找到更新后介質(zhì)中的成像點(diǎn)位置(即運(yùn)動(dòng)學(xué)偏移),運(yùn)動(dòng)學(xué)偏移完成后相當(dāng)于更新了成像點(diǎn)的局部運(yùn)動(dòng)學(xué)信息。這樣做的優(yōu)點(diǎn)是緩解了算法實(shí)際應(yīng)用中需要聯(lián)合反演多種物理參數(shù)的壓力。二維典型理論數(shù)據(jù)試驗(yàn)結(jié)果表明,僅利用地表觀測(cè)的px殘差確實(shí)可以完成VTI介質(zhì)中三個(gè)各向異性參數(shù)的估計(jì),將其與VTI介質(zhì)中的運(yùn)動(dòng)學(xué)偏移/反偏移相結(jié)合,是一種非常實(shí)用的各向異性參數(shù)估算策略。本文考慮到qP波依然是大部分實(shí)際地震采集中需要面對(duì)的主要波動(dòng)現(xiàn)象,基于文獻(xiàn)[7]提出的VTI介質(zhì)擬聲波程函方程的漢密爾頓形式,推導(dǎo)了qP波斜率層析所需的Fréchet偏導(dǎo)數(shù)公式,給出了VTI介質(zhì)中的運(yùn)動(dòng)學(xué)反偏移與VTI介質(zhì)qP波斜率層析相結(jié)合的處理流程。需要強(qiáng)調(diào)的是,本文所討論的斜率層析只用到了px參數(shù)關(guān)于各向異性三參數(shù)的偏導(dǎo)數(shù),因此地表觀測(cè)位置(sx,rx)和旅行時(shí)t關(guān)于地下各向異性介質(zhì)的偏導(dǎo)數(shù)將被略去。

1 二維VTI介質(zhì)qP波斜率層析偏導(dǎo)數(shù)推導(dǎo)

任何一種層析反演算法都必須建立數(shù)據(jù)空間殘差與模型空間殘差之間的線性關(guān)系,即:

Δd=FΔm

(1)

基于方程(1),利用觀測(cè)到的數(shù)據(jù)殘差就可以計(jì)算出模型空間的更新量,達(dá)到更新初始模型的目的。其關(guān)鍵步驟是建立方程(1)所示的Fréchet偏導(dǎo)數(shù)矩陣(以下簡(jiǎn)稱F矩陣)。注意,在二維VTI介質(zhì)qP波立體層析中,數(shù)據(jù)空間為d=[sxrxtpsxprx],其中sx、rx為炮檢點(diǎn)坐標(biāo);t為走時(shí);psx、prx為炮檢點(diǎn)處的射線參數(shù)水平分量。地下模型空間為m=[x0z0θsθrvvεδ],其中x0、z0為地下反射點(diǎn)x的坐標(biāo);θs、θr分別代表從炮點(diǎn)一側(cè)與檢波點(diǎn)一側(cè)出射的散射相角;vv、ε、δ代表VTI介質(zhì)中的三個(gè)各向異性參數(shù)。

(2)

方程(2)展示了二維VTI立體層析Fréchet偏導(dǎo)數(shù)矩陣,其中第一行、第二行、第三行分別為旅行時(shí)t、炮點(diǎn)坐標(biāo)sx、檢波點(diǎn)坐標(biāo)rx關(guān)于立體層析模型空間各分量的偏導(dǎo)數(shù),本文不予討論。考慮到本文討論的斜率層析僅限于反演各向異性參數(shù),這里僅關(guān)注炮點(diǎn)射線參數(shù)水平分量psx、檢波點(diǎn)射線參數(shù)水平分量prx關(guān)于三個(gè)各向異性參數(shù)的偏導(dǎo)數(shù),關(guān)于散射相角和反射點(diǎn)位置的偏導(dǎo)數(shù)將不予討論。因此,本文使用的VTI介質(zhì)qP波斜率層析Fréchet偏導(dǎo)數(shù)矩陣退化為:

(3)

注意,公式(2)、公式(3)中的κ為針對(duì)矩陣中每一行的數(shù)量級(jí)不同設(shè)置的均衡系數(shù),目的是保證每一行的數(shù)值范圍基本一致。

WANG等[8]基于重新參數(shù)化后的二維VTI介質(zhì)擬聲波程函方程和射線擾動(dòng)理論[9],推導(dǎo)了二維qP波立體層析Fréchet偏導(dǎo)數(shù)矩陣。本文采用類似的方式推導(dǎo)出適用于二維qP波的斜率層析Fréchet偏導(dǎo)數(shù)矩陣。二維VTI介質(zhì)擬聲波程函方程如下[7]:

(4)

式中,vv是介質(zhì)垂向速度;vnmo是NMO速度。根據(jù)文獻(xiàn)[7]可知各向異性參數(shù)之間有如下關(guān)系:η=(ε-δ)/(1+2δ)和vnmo=vv(1+2δ)1/2。為了保證變量之間的獨(dú)立性,不妨采用一種等價(jià)的程函方程,其形式如下:

(5)

因此,二維VTI介質(zhì)下的qP波漢密爾頓算子可以寫為:

(6)

相應(yīng)地,與射線參數(shù)水平分量和垂直分量有關(guān)的射線追蹤一階微分方程組為:

(7)

式中:σ為沿射線方向的滑動(dòng)積分變量。

對(duì)方程(7)做全微分,可得到由初始相空間擾動(dòng)(Δpx,Δpz)及各向異性參數(shù)擾動(dòng)(Δvv,Δε,Δδ)所引起的射線軌跡上的相空間擾動(dòng):

(8)

其中,

(9a)

根據(jù)射線擾動(dòng)理論[9],關(guān)于模型空間3個(gè)各向異性參數(shù)的擾動(dòng)漢密爾頓形式可以寫為:

(9b)

(8)式中后三項(xiàng)形式如下:

(9c)

結(jié)合(9a)、(9b)、(9c)式,利用射線傳播矩陣?yán)碚揫10]即可求得方程(8)的傳播矩陣解:

(10)

式中:σ0代表射線出射位置,σ1代表射線傳播的最終位置,Π(σ1,σ0)表示從射線出射位置到最終位置的傳播矩陣。由于本文僅利用px作為數(shù)據(jù)殘差反演各向異性參數(shù),而不反演出射位置處的散射角信息,因此(10)式可進(jìn)一步簡(jiǎn)化為:

(11)

利用(11)式即可獲得二維VTI各向異性介質(zhì)qP波斜率層析所需的數(shù)據(jù)分量與模型空間各分量之間的一階擾動(dòng)關(guān)系,也意味著斜率層析Fréchet偏導(dǎo)數(shù)矩陣F得到建立,后續(xù)通過求解線性方程組(1),就可以反演各向異性參數(shù)。

2 二維VTI介質(zhì)qP波斜率層析中運(yùn)動(dòng)學(xué)反偏移及實(shí)現(xiàn)流程

首先介紹文獻(xiàn)[5]和文獻(xiàn)[6]提出的運(yùn)動(dòng)學(xué)反偏移概念以及由此導(dǎo)出的各向同性介質(zhì)中的斜率層析方法。如圖1所示,v是地下成像點(diǎn)(x,z)處的速度,ts、tr分別是地下成像點(diǎn)(x,z)到炮點(diǎn)s、檢波點(diǎn)r的單程走時(shí),h是半偏移距,m是炮點(diǎn)s與檢波點(diǎn)r的中點(diǎn)。圖1a表示地下模型信息與地表數(shù)據(jù)信息的幾何關(guān)系;圖1b為共偏移距偏移成像剖面,在該剖面上拾取構(gòu)造傾角ξ;圖1c為偏移距域共成像點(diǎn)道集(CIG),在該道集上拾取剩余曲率(RMO)的斜率信息tanφ。基于偏移后共偏移距成像剖面上搜索到的構(gòu)造傾角ξ,向地表發(fā)射線找到對(duì)應(yīng)的炮檢點(diǎn)坐標(biāo),對(duì)應(yīng)的射線水平參數(shù)分量設(shè)為psx和prx。

圖1 運(yùn)動(dòng)學(xué)反偏移原理(據(jù)文獻(xiàn)[5])a 地下模型信息與地表數(shù)據(jù)信息的幾何關(guān)系; b 偏移后共偏移距剖面; c 偏移后共成像點(diǎn)道集

(12a)

(12b)

式中,α=2vcosβcosξ,β是反射角,它與兩個(gè)散射角θs、θr以及構(gòu)造傾角ξ的關(guān)系如下:

cosθs+cosθr=2cosβcosξ

(12c)

將上述各向同性介質(zhì)中的運(yùn)動(dòng)學(xué)反偏移推廣到VTI介質(zhì)非常簡(jiǎn)單,只需要將各向同性介質(zhì)的運(yùn)動(dòng)學(xué)射線追蹤更換為VTI介質(zhì)的運(yùn)動(dòng)學(xué)射線追蹤,將各向同性克希霍夫疊前深度偏移更換為VTI介質(zhì)克希霍夫疊前深度偏移即可。

圖2 VTI介質(zhì)的運(yùn)動(dòng)學(xué)反偏移及斜率層析流程

3 理論數(shù)據(jù)測(cè)試

利用二維VTI多層背斜模型數(shù)據(jù),驗(yàn)證了本文計(jì)算公式的正確性和實(shí)現(xiàn)流程的可靠性。圖3為6層背斜VTI參數(shù)模型,基于該模型進(jìn)行VTI波動(dòng)方程正演模擬。正演數(shù)據(jù)共601炮,每炮801道,炮間距20m,道間距10m,最大偏移距4000m,記錄時(shí)間為4s。圖4展示了正演數(shù)據(jù)的3個(gè)單炮記錄。大量實(shí)際資料分析表明,δ參數(shù)通常是地震響應(yīng)較弱的一個(gè)參數(shù),目前工業(yè)界主要通過利用井資料對(duì)各向同性深度偏移成像剖面進(jìn)行標(biāo)定后,再利用真實(shí)深度與成像深度的差轉(zhuǎn)換得到δ剖面,這是一個(gè)比較成熟的技術(shù)[13]。因此,本文重點(diǎn)討論如何利用提取出的地表px參數(shù)殘差反演垂直速度vv和ε,將δ剖面設(shè)置為正確值,不對(duì)其進(jìn)行反演。

圖5展示了初始vv剖面和初始ε剖面(常梯度)。為了檢驗(yàn)本文算法聯(lián)合反演這兩個(gè)參數(shù)的能力,這里使用的初始vv各向同性模型為經(jīng)過傾角時(shí)差校正(DMO)速度分析后轉(zhuǎn)到深度域、并做了大尺度平滑后的速度模型。

基于上述兩個(gè)剖面以及正確的δ剖面,實(shí)施第1輪VTI介質(zhì)克希霍夫深度偏移,獲得初始小偏移距疊加成像剖面(偏移距范圍0～500m)和若干初始共成像點(diǎn)道集(圖6)。基于圖6a所示的初始成像剖面,拾取了700多個(gè)數(shù)據(jù)點(diǎn)(圖6b),然后從每個(gè)數(shù)據(jù)點(diǎn)出發(fā),利用結(jié)構(gòu)張量算法從零偏移距到最大偏移距每隔5個(gè)偏移距向地表發(fā)出一對(duì)射線,這樣得到的數(shù)據(jù)空間一共有1.4×104個(gè)射線對(duì)。圖6d 展示了1500m偏移距的偏移成像剖面局部放大結(jié)果,以及CDP1430處的CIG,圖中水平方向紅線對(duì)應(yīng)深度為820m。兩條紅線的交點(diǎn)為A,利用結(jié)構(gòu)張量算法可以快速地在1500m偏移距成像剖面內(nèi)搜索出關(guān)于A點(diǎn)的構(gòu)造傾角ξ,同時(shí)在CDP1430處的CIG內(nèi)搜索出RMO曲線的斜率信息tanφ。

圖3 6層背斜VTI模型參數(shù)a 垂直速度vv剖面; b ε剖面; c δ剖面

圖4 二維VTI介質(zhì)波動(dòng)方程模擬的3個(gè)單炮記錄

圖5 初始垂直速度vv(a)與初始ε(b)剖面

圖6 VTI介質(zhì)中通過運(yùn)動(dòng)學(xué)反偏移獲取數(shù)據(jù)空間的實(shí)現(xiàn)過程a 初始小偏移距疊加成像剖面(偏移距范圍0～500m); b 基于圖6a拾取的數(shù)據(jù)點(diǎn)位置; c 利用圖5所示初始vv、初始ε以及正確的δ剖面進(jìn)行VTI克希霍夫疊前深度偏移獲得的共成像點(diǎn)道集; d 1500m共偏移距成像剖面局部放大結(jié)果及CDP1430處的CIG道集; e 坐標(biāo)為7100m的共炮點(diǎn)記錄; f 坐標(biāo)為8600m的共檢波點(diǎn)記錄

圖9、圖10分別為基于初始模型和最終反演模型的CIG道集和VTI疊前深度偏移結(jié)果。圖9a對(duì)比了5km處基于初始模型和最終反演模型得到的CIG道集,圖9b 對(duì)比了10km處基于初始模型和最終反演模型得到的CIG道集,可見基于最終反演模型獲得的CIG均得到了有效拉平。圖10對(duì)比了基于初始模型和最終模型的PSDM偏移剖面,可以看出中深部的成像品質(zhì)大為提升,說明本文方法的反演結(jié)果能夠滿足疊前深度偏移成像的要求。

圖7 聯(lián)合反演結(jié)果a vv剖面; b ε剖面; c 誤差泛函下降曲線

圖8 不同水平位置vv和ε抽線對(duì)比a 5km處vv抽線對(duì)比; b 10km處vv抽線對(duì)比; c 5km處ε抽線對(duì)比; d 10km處ε抽線對(duì)比

圖9 5km(a)和10km(b)處的參數(shù)模型更新前后CIG道集對(duì)比

圖10 基于初始模型(a)和最終模型(b)PSDM偏移剖面對(duì)比

4 結(jié)論和討論

本文將各向同性介質(zhì)中的斜率層析方法推廣到二維VTI介質(zhì)qP波的情況。基于射線擾動(dòng)理論推導(dǎo)了地表px參數(shù)關(guān)于地下三個(gè)各向異性參數(shù)的偏導(dǎo)數(shù),建立了適合于二維VTI介質(zhì)qP波斜率層析的線性方程組,并將各向同性介質(zhì)中的運(yùn)動(dòng)學(xué)偏移/反偏移策略推廣到VTI介質(zhì)中,最終實(shí)現(xiàn)了二維VTI介質(zhì)qP波斜率層析方法。研究結(jié)果表明,VTI介質(zhì)中的運(yùn)動(dòng)學(xué)偏移/反偏移與二維VTI介質(zhì)qP波斜率層析線性方程組相結(jié)合,可以很好地實(shí)現(xiàn)二維VTI介質(zhì)中的各向異性參數(shù)重建。該方法也可以方便地推廣到三維或者TTI介質(zhì)。

需要指出的是,無論是各向同性介質(zhì)還是各向異性介質(zhì)中的斜率層析方法,對(duì)初始模型都有一定要求:初始成像體的品質(zhì)能夠保證CIG上提取的RMO斜率tanφ和共偏移距成像剖面上提取的構(gòu)造傾角ξ都是可靠的。如果這個(gè)要求無法滿足,那么后續(xù)通過運(yùn)動(dòng)學(xué)反偏移獲取的地表px參數(shù)精度將難以保證。

其次,在斜率層析中忽略旅行時(shí)對(duì)反射點(diǎn)位置和散射角的導(dǎo)數(shù)項(xiàng)不會(huì)降低算法的適用性。這是因?yàn)?如果在立體層析反演中能夠?qū)?shù)據(jù)殘差通過合理的方式轉(zhuǎn)移到某一個(gè)數(shù)據(jù)分量上,則會(huì)使得數(shù)據(jù)空間得到合理的壓縮,這對(duì)提高反演的穩(wěn)定性是有益的。這一點(diǎn)已經(jīng)在前人發(fā)展的斜率層析和控制定向?qū)游龇椒ㄖ械玫搅俗C實(shí)。同時(shí),本文方法沒有將反射點(diǎn)位置和散射角納入到模型空間內(nèi),實(shí)際上又壓縮了傳統(tǒng)斜率層析的模型空間,進(jìn)一步提高了反演的穩(wěn)定性,有利于本文方法的推廣應(yīng)用。