谷巍巍,楊 鍇,王 瀟

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

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

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

1 二維VTI介質qP波斜率層析偏導數推導

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

Δd=FΔm

(1)

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

(2)

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

(3)

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

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

(4)

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

(5)

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

(6)

相應地,與射線參數水平分量和垂直分量有關的射線追蹤一階微分方程組為:

(7)

式中:σ為沿射線方向的滑動積分變量。

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

(8)

其中,

(9a)

根據射線擾動理論[9],關于模型空間3個各向異性參數的擾動漢密爾頓形式可以寫為:

(9b)

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

(9c)

結合(9a)、(9b)、(9c)式,利用射線傳播矩陣理論[10]即可求得方程(8)的傳播矩陣解:

(10)

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

(11)

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

2 二維VTI介質qP波斜率層析中運動學反偏移及實現流程

首先介紹文獻[5]和文獻[6]提出的運動學反偏移概念以及由此導出的各向同性介質中的斜率層析方法。如圖1所示,v是地下成像點(x,z)處的速度,ts、tr分別是地下成像點(x,z)到炮點s、檢波點r的單程走時,h是半偏移距,m是炮點s與檢波點r的中點。圖1a表示地下模型信息與地表數據信息的幾何關系;圖1b為共偏移距偏移成像剖面,在該剖面上拾取構造傾角ξ;圖1c為偏移距域共成像點道集(CIG),在該道集上拾取剩余曲率(RMO)的斜率信息tanφ?；谄坪蠊财凭喑上衿拭嫔纤阉鞯降臉嬙靸A角ξ,向地表發(fā)射線找到對應的炮檢點坐標,對應的射線水平參數分量設為psx和prx。

圖1 運動學反偏移原理(據文獻[5])a 地下模型信息與地表數據信息的幾何關系; b 偏移后共偏移距剖面; c 偏移后共成像點道集

(12a)

(12b)

式中,α=2vcosβcosξ,β是反射角,它與兩個散射角θs、θr以及構造傾角ξ的關系如下:

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

(12c)

將上述各向同性介質中的運動學反偏移推廣到VTI介質非常簡單,只需要將各向同性介質的運動學射線追蹤更換為VTI介質的運動學射線追蹤,將各向同性克希霍夫疊前深度偏移更換為VTI介質克?；舴虔B前深度偏移即可。

圖2 VTI介質的運動學反偏移及斜率層析流程

3 理論數據測試

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

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

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

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

圖4 二維VTI介質波動方程模擬的3個單炮記錄

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

圖6 VTI介質中通過運動學反偏移獲取數據空間的實現過程a 初始小偏移距疊加成像剖面(偏移距范圍0～500m); b 基于圖6a拾取的數據點位置; c 利用圖5所示初始vv、初始ε以及正確的δ剖面進行VTI克?；舴虔B前深度偏移獲得的共成像點道集; d 1500m共偏移距成像剖面局部放大結果及CDP1430處的CIG道集; e 坐標為7100m的共炮點記錄; f 坐標為8600m的共檢波點記錄

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

圖7 聯合反演結果a vv剖面; b ε剖面; c 誤差泛函下降曲線

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

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

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

4 結論和討論

本文將各向同性介質中的斜率層析方法推廣到二維VTI介質qP波的情況。基于射線擾動理論推導了地表px參數關于地下三個各向異性參數的偏導數,建立了適合于二維VTI介質qP波斜率層析的線性方程組,并將各向同性介質中的運動學偏移/反偏移策略推廣到VTI介質中,最終實現了二維VTI介質qP波斜率層析方法。研究結果表明,VTI介質中的運動學偏移/反偏移與二維VTI介質qP波斜率層析線性方程組相結合,可以很好地實現二維VTI介質中的各向異性參數重建。該方法也可以方便地推廣到三維或者TTI介質。

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

其次,在斜率層析中忽略旅行時對反射點位置和散射角的導數項不會降低算法的適用性。這是因為,如果在立體層析反演中能夠將數據殘差通過合理的方式轉移到某一個數據分量上,則會使得數據空間得到合理的壓縮,這對提高反演的穩(wěn)定性是有益的。這一點已經在前人發(fā)展的斜率層析和控制定向層析方法中得到了證實。同時,本文方法沒有將反射點位置和散射角納入到模型空間內,實際上又壓縮了傳統斜率層析的模型空間,進一步提高了反演的穩(wěn)定性,有利于本文方法的推廣應用。