起伏地表下基于拋物插值的最短路徑射線追蹤

桑運(yùn)云，孫軍曉，焦淑萍，金艷萍，陳院生

(1.中國(guó)石油天然氣集團(tuán)公司東方地球物理勘探有限責(zé)任公司研究院華北分院，河北任丘062552;2.中國(guó)石油天然氣集團(tuán)公司東方地球物理勘探有限責(zé)任公司研究院處理中心,河北涿州072751)

隨著地震勘探技術(shù)的發(fā)展和勘探規(guī)模的擴(kuò)大，起伏地表?xiàng)l件下的地震勘探越來(lái)越普遍[1-2],通過(guò)有效的近地表速度建模方法解決靜校正問(wèn)題成為地震資料處理的重點(diǎn)。初至波層析速度建模由于對(duì)初始模型的依賴性小，可以直接由初至波旅行時(shí)差反演近地表速度，精確地反演速度場(chǎng)的長(zhǎng)波長(zhǎng)分量，因而是常用的近地表速度建模方法。

初至波層析的關(guān)鍵是初至波旅行時(shí)計(jì)算和射線路徑的追蹤，關(guān)系到微商矩陣的準(zhǔn)確性和走時(shí)殘差誤差的大小，因此需要相對(duì)精確的初至波旅行時(shí)射線追蹤方法。傳統(tǒng)的射線追蹤方法各有優(yōu)缺點(diǎn)：試射法通過(guò)迭代修改初始入射角獲得了準(zhǔn)確的射線路徑，但計(jì)算比較耗時(shí);彎曲法利用最小走時(shí)條件描述射線走時(shí)方程，然后迭代修改射線得出正確的射線路徑，但常常會(huì)陷入局部收斂;波前法將地下介質(zhì)網(wǎng)格化，利用網(wǎng)格節(jié)點(diǎn)連線近似射線路徑，但是射線路徑常常為折線。起伏地表?xiàng)l件下射線追蹤的關(guān)鍵技術(shù)是起伏地表的網(wǎng)格化。岳玉波等[3]研究的初值射線追蹤方法首先根據(jù)三次樣條函數(shù)及線性函數(shù)來(lái)描述起伏地表以及地下界面,然后給定射線的初始點(diǎn)和初始方向,逐步計(jì)算射線路徑,因此計(jì)算精度在很大程度上依賴于逐步迭代追蹤的步長(zhǎng)。Cassell[4]采用網(wǎng)格剖分方法將地下介質(zhì)劃分為很多小塊，每塊速度為常數(shù)，按照一定初始條件從震源處逐個(gè)網(wǎng)格進(jìn)行射線追蹤，得到的射線路徑在每個(gè)網(wǎng)格內(nèi)都為直線段。Langan等[5]在此基礎(chǔ)上將網(wǎng)格內(nèi)的速度用梯度表示，用曲線表示射線路徑，提高了效率和精度。但是以上方法都不能精確地實(shí)現(xiàn)炮點(diǎn)和檢波點(diǎn)的點(diǎn)對(duì)點(diǎn)射線路徑追蹤，只能近似表示兩者之間的射線路徑。Moser[6]提出的最短路徑射線追蹤方法采用節(jié)點(diǎn)樹(shù)的形式實(shí)現(xiàn)了初至波射線追蹤，但得到的射線路徑往往成“之”字形。張建中等[7]利用線性插值方法改進(jìn)了最短路徑射線追蹤方法，克服了射線路徑只能走網(wǎng)格節(jié)點(diǎn)的缺陷，實(shí)現(xiàn)了動(dòng)態(tài)網(wǎng)格射線追蹤，但是在近地表?xiàng)l件下走時(shí)變化往往不滿足線性規(guī)律。

我們將初至波路徑射線追蹤方法——基于拋物旅行時(shí)插值(PTI)的最短路徑射線追蹤方法[8]擴(kuò)展到起伏地表?xiàng)l件下，將起伏地表及地下界面精確地投影到剖分網(wǎng)格節(jié)點(diǎn)上，實(shí)現(xiàn)了起伏地表?xiàng)l件下的初至波路徑的射線追蹤。

1 基于PTI的最短路徑射線追蹤方法原理

基于PTI的最短路徑射線追蹤方法(以下簡(jiǎn)稱PTISPR法)分兩步進(jìn)行：①正向波前擴(kuò)展，即從震源開(kāi)始計(jì)算整個(gè)速度場(chǎng)的所有節(jié)點(diǎn)的最小走時(shí);②反向追蹤路徑，即從接收點(diǎn)開(kāi)始向震源反向追蹤射線路徑。其核心思想是：首先利用Dijkstra[9]算法計(jì)算某節(jié)點(diǎn)旅行時(shí)，然后在計(jì)算過(guò)程中搜索滿足拋物插值的節(jié)點(diǎn)并進(jìn)行插值，將最小旅行時(shí)作為該計(jì)算節(jié)點(diǎn)的最終時(shí)間。射線追蹤過(guò)程是利用拋物旅行時(shí)插值從接收點(diǎn)網(wǎng)格開(kāi)始一直追蹤到炮點(diǎn)網(wǎng)格。若A,B,C為同一邊界上已知走時(shí)的3個(gè)節(jié)點(diǎn)，其坐標(biāo)分別為(xA,zA)，(xB,zB)，(xC,zC)，對(duì)應(yīng)的走時(shí)為tA，tB，tC，則它們構(gòu)成了節(jié)點(diǎn)D的插值段，求該邊界上從R點(diǎn)(xR,zR)到達(dá)D點(diǎn)(xD,zD)的最小旅行時(shí)tD和射線路徑所利用的插值公式是[8]

(1)

其中，LAB=xA-xB，LAC=xA-xC，LBC=xB-xC;s為慢度。

2 起伏地表?xiàng)l件下PTISPR法射線追蹤的實(shí)現(xiàn)

2.1 待求旅行時(shí)的節(jié)點(diǎn)數(shù)量

PTISPR法射線追蹤分向前和向后兩個(gè)過(guò)程，向前計(jì)算每個(gè)節(jié)點(diǎn)的旅行時(shí)[6]，向后追蹤射線路徑。將起伏地表網(wǎng)格化，如圖1所示，地表上速度為0，地表面以上速度區(qū)域?yàn)闊o(wú)效區(qū)域，不需要計(jì)算節(jié)點(diǎn)時(shí)間。因此，首先需要找出速度場(chǎng)的有效區(qū)域以及要計(jì)算旅行時(shí)的節(jié)點(diǎn)數(shù)。如圖2所示，節(jié)點(diǎn)A和節(jié)點(diǎn)B為相鄰兩個(gè)橫向采樣點(diǎn)，求取兩者所在的網(wǎng)格范圍內(nèi)A節(jié)點(diǎn)所屬的無(wú)效節(jié)點(diǎn)數(shù)(即A節(jié)點(diǎn)及其以上采樣點(diǎn)對(duì)應(yīng)的邊界節(jié)點(diǎn)和角節(jié)點(diǎn)無(wú)效節(jié)點(diǎn)數(shù)目)。

圖1 起伏地表圖示

1) 節(jié)點(diǎn)A和節(jié)點(diǎn)B縱向采樣值相等(圖2a)。圖2a中①表示A和B都為角點(diǎn)節(jié)點(diǎn);②表示A為角點(diǎn)節(jié)點(diǎn)，B為下邊界節(jié)點(diǎn);③表示A為下邊界節(jié)點(diǎn)，B為角點(diǎn)節(jié)點(diǎn)或下邊界節(jié)點(diǎn)。

2) 節(jié)點(diǎn)A的縱向采樣值大于節(jié)點(diǎn)B的縱向采樣值(圖2b)。根據(jù)A和B的坐標(biāo)位置，可以得到連接A,B的直線段，然后依次求取與A及A以上節(jié)點(diǎn)對(duì)應(yīng)的右邊界的交點(diǎn)，若節(jié)點(diǎn)的橫向位置小于交點(diǎn)橫向位置，則為無(wú)效節(jié)點(diǎn)。對(duì)于A節(jié)點(diǎn)對(duì)應(yīng)的下邊界節(jié)點(diǎn)，則根據(jù)A的縱向坐標(biāo)計(jì)算。

3) 節(jié)點(diǎn)A的縱向采樣值小于節(jié)點(diǎn)B的縱向采樣值(圖2c)。計(jì)算原理同2)，僅僅當(dāng)節(jié)點(diǎn)的橫向位置大于交點(diǎn)橫向位置時(shí)，為無(wú)效節(jié)點(diǎn)。

圖2 起伏地表在一個(gè)網(wǎng)格內(nèi)的不同分布情況a 節(jié)點(diǎn)A,B縱向采樣值相等; b 節(jié)點(diǎn)A縱向采樣值大于節(jié)點(diǎn)B; c 節(jié)點(diǎn)A縱向采樣值小于節(jié)點(diǎn)B

2.2 速度場(chǎng)的處理

受網(wǎng)格化和起伏地表的影響，在起伏地表范圍內(nèi)常常會(huì)出現(xiàn)角節(jié)點(diǎn)為無(wú)效節(jié)點(diǎn)而邊界節(jié)點(diǎn)為有效節(jié)點(diǎn)的情況，在求取網(wǎng)格或路徑速度時(shí)就會(huì)出現(xiàn)某一節(jié)點(diǎn)采樣值對(duì)應(yīng)的速度為0的結(jié)果，進(jìn)而影響網(wǎng)格內(nèi)速度和某一射線路徑速度的計(jì)算。解決方法是：當(dāng)某一節(jié)點(diǎn)采樣值對(duì)應(yīng)的速度為0時(shí)，說(shuō)明此計(jì)算節(jié)點(diǎn)為邊界節(jié)點(diǎn)，因此首先判斷其屬于右邊界節(jié)點(diǎn)還是下邊界節(jié)點(diǎn)，若為右邊界節(jié)點(diǎn)則將采樣值橫向加1，縱向不變，得到的采樣值對(duì)應(yīng)的速度作為該節(jié)點(diǎn)的速度；若為下邊界節(jié)點(diǎn)則將采樣值縱向加1，橫向不變，得到的采樣值對(duì)應(yīng)的速度作為該節(jié)點(diǎn)的速度。

2.3 插值

拋物插值是根據(jù)某一網(wǎng)格內(nèi)已知的相鄰3個(gè)節(jié)點(diǎn)時(shí)間求其它節(jié)點(diǎn)的時(shí)間。在起伏地表所處的網(wǎng)格邊界，相鄰3個(gè)節(jié)點(diǎn)會(huì)出現(xiàn)3種情況：①相鄰3個(gè)節(jié)點(diǎn)時(shí)間全部已知;②相鄰兩個(gè)節(jié)點(diǎn)時(shí)間已知，另一個(gè)為無(wú)效節(jié)點(diǎn);③有兩個(gè)節(jié)點(diǎn)為無(wú)效節(jié)點(diǎn)，僅有一個(gè)節(jié)點(diǎn)時(shí)間已知。因此，在起伏地表射線追蹤過(guò)程中，需要對(duì)網(wǎng)格邊界上的節(jié)點(diǎn)進(jìn)行判斷，對(duì)應(yīng)上述3種情況分別進(jìn)行拋物插值、線性插值和直接計(jì)算兩節(jié)點(diǎn)的直線時(shí)間。

2.4 起伏地表所在網(wǎng)格內(nèi)射線路徑的確定

射線路徑的確定是從接收點(diǎn)開(kāi)始向炮點(diǎn)位置反向追蹤的過(guò)程，當(dāng)追蹤到炮點(diǎn)所處的網(wǎng)格時(shí)，需要判斷當(dāng)前接收點(diǎn)和炮點(diǎn)是否在同一個(gè)網(wǎng)格范圍內(nèi)。由于起伏地表的存在，射線路徑在同一網(wǎng)格內(nèi)經(jīng)過(guò)的網(wǎng)格點(diǎn)數(shù)往往大于2，使射線路徑出現(xiàn)折線形狀。處理方法是：當(dāng)接收點(diǎn)和炮點(diǎn)坐標(biāo)橫向距離不大于橫向采樣間隔時(shí)，若縱向距離不大于縱向采樣間隔的1/2，則直接將接收點(diǎn)修正為炮點(diǎn);同理，當(dāng)接收點(diǎn)和炮點(diǎn)坐標(biāo)縱向距離不大于縱向采樣間隔時(shí)，若橫向距離不大于橫向采樣間隔的1/2，則直接將接收點(diǎn)修正為炮點(diǎn)。

2.5 模型試算

2.5.1 起伏地表水平層狀模型

圖3是加入了起伏地表的4層水平模型速度場(chǎng)(600m×600m)。將模型網(wǎng)格化，橫向和深度采樣間隔都為4m，網(wǎng)格數(shù)為151×151。炮點(diǎn)位于起伏地表橫向300m處，模型底部和兩個(gè)邊界接收，共83道，道間距為20m。其中底部31道，第1道坐標(biāo)為(0,599m);兩個(gè)邊界各26道，第1個(gè)接收點(diǎn)位置坐標(biāo)分別為(0,80m)，(600m,80m)。圖4是采用PTISPR法由上述觀測(cè)系統(tǒng)得到的等時(shí)線和射線路徑，從等時(shí)線圖上可以明顯看到3個(gè)分界面處波前面發(fā)生的變化(圖4a)，從射線路徑圖上可以看到直達(dá)波和折射波兩種不同類型的初至波路徑(圖4b)。

圖3 起伏地表水平層狀模型

圖4 起伏地表水平層狀模型等時(shí)線(a)和射線路徑(b)

2.5.2 起伏地表連續(xù)介質(zhì)模型

圖5為起伏地表連續(xù)介質(zhì)速度模型(1500m×500m)，速度場(chǎng)從淺到深的變化范圍為300～2800m/s，網(wǎng)格節(jié)點(diǎn)數(shù)為101×101，橫向采樣間隔15m，深度采樣間隔5m。地面地震觀測(cè)系統(tǒng)：炮點(diǎn)和接收點(diǎn)都位于起伏地表，模型中間位置放炮，兩邊接收，第1道橫向位置0，道間隔30m，共51道。圖6顯示了應(yīng)用PTISPR法由上述觀測(cè)系統(tǒng)得到的回折波射線路徑和對(duì)應(yīng)的等時(shí)線，可見(jiàn)起伏地表連續(xù)介質(zhì)模型對(duì)應(yīng)的等時(shí)線是一系列半徑隨深度增加的圓弧，射線是從震源出發(fā)向下到達(dá)某一深度后又向上拐回地面到達(dá)觀測(cè)點(diǎn)的兩套圓弧射線，稱為“回折波”。

圖5 起伏地表連續(xù)介質(zhì)模型

圖6 起伏地表連續(xù)介質(zhì)模型等時(shí)線(a)和射線路徑(b)

3 起伏地表?xiàng)l件下PTISPR法初至波走時(shí)層析

對(duì)于某一模型速度場(chǎng)，初至波層析速度分析大體流程為：

1) 建立初始速度模型;

2) 分別對(duì)初始速度場(chǎng)和模型速度場(chǎng)進(jìn)行初至波路徑射線追蹤，由此得到初至波旅行時(shí)(將初始速度場(chǎng)和模型速度場(chǎng)的旅行時(shí)相減得到的走時(shí)殘差)和層析所需的微商矩陣(射線在每個(gè)網(wǎng)格內(nèi)的路徑長(zhǎng)度);

3) 建立層析方程組[10]AΔs=Δt,其中A為微商矩陣，Δt為初始速度場(chǎng)和模型速度場(chǎng)的初至旅行時(shí)差，Δs為慢度更新量;

4) 求解層析方程組得到慢度更新量，更新初始模型，然后轉(zhuǎn)至步驟2);

5) 如此迭代，直到更新后的速度場(chǎng)滿足需要(兩次迭代旅行時(shí)變化很小，速度改變量很小)為止。

3.1 模型試算

圖7所示為含起伏地下界面起伏地表模型的真實(shí)速度場(chǎng)(1500m×1500m)，采樣點(diǎn)數(shù)為101×301，橫向采樣間隔15m，深度采樣間隔5m。地面地震觀測(cè)系統(tǒng)：炮點(diǎn)和檢波點(diǎn)都位于起伏地表，第1炮橫向位置0，炮間隔250m，共7炮;每炮

101道接收，初始道橫向位置0，道間隔15m。初始速度場(chǎng)深層速度與真實(shí)速度場(chǎng)一致，淺層是梯度模型；迭代18次后平均走時(shí)殘差由-160.0000ms下降到-0.0078ms，得到圖8所示的層析建模速度場(chǎng)。與真實(shí)速度場(chǎng)(圖7)比較可見(jiàn)，近地表低頻信息得到了很好的反演，基本上反映了低、降速帶內(nèi)速度變化的趨勢(shì)。圖9為真實(shí)速度和層析速度的等時(shí)線和射線路徑，可見(jiàn)兩者基本吻合。經(jīng)過(guò)模型試算，證明起伏地表PTISPR法可以用于近地表建模。

3.2 實(shí)際資料試算

某地形起伏地區(qū)實(shí)際地震資料的炮檢排列沿南北走向，每個(gè)排列50炮。沿南北向網(wǎng)格間隔100m，深度采樣間隔5m，網(wǎng)格化后采樣點(diǎn)數(shù)為321×121。根據(jù)微測(cè)井和小折射資料，構(gòu)建沿深度遞增的梯度初始模型，梯度變化量為25m/s。圖10a 為初始速度場(chǎng)，圖10b為層析建模速度場(chǎng)。由圖10可見(jiàn)，通過(guò)對(duì)初始模型(圖10a)進(jìn)行迭代處理，得到的層析建模速度場(chǎng)(圖10b)從南到北高程變化劇烈，呈現(xiàn)南低北高的趨勢(shì)，北側(cè)低、降速帶范圍埋深更淺，符合實(shí)際地形特征。利用層析建模得到的速度場(chǎng)計(jì)算靜校正量并對(duì)原始炮集進(jìn)行靜校正[11]，圖11為某單炮記錄初至層析靜校正前、后的結(jié)果。從圖11可見(jiàn)，靜校正前(圖11a)同相軸偏離雙曲規(guī)律，抖動(dòng)比較劇烈，初至同相軸不連續(xù);經(jīng)過(guò)靜校正后(圖11b)，抖動(dòng)現(xiàn)象消除，同相軸更加符合雙曲規(guī)律，初至同相軸連續(xù)性更好。實(shí)際資料試算結(jié)果證明本文提出的PTISPR法適用于起伏地表情況下地震資料的初至層析靜校正處理。

圖7 含起伏地下界面的起伏地表模型真實(shí)速度場(chǎng)

圖8 含起伏地下界面的起伏地表模型層析速度場(chǎng)

圖9 含起伏地下界面的起伏地表模型真實(shí)速度和層析速度對(duì)應(yīng)的等時(shí)線和射線路徑a 真實(shí)速度等時(shí)線; b 層析速度等時(shí)線; c 真實(shí)速度射線路徑; d 層析速度射線路徑

圖10 起伏地表區(qū)實(shí)際地震資料層析初始速度場(chǎng)(a)和建模速度場(chǎng)(b)

圖11 起伏地表區(qū)實(shí)際地震資料單炮記錄靜校正前(a)和初至層析靜校正后(b)

4 結(jié)束語(yǔ)

射線追蹤是走時(shí)層析的核心技術(shù)。良好的初至波射線追蹤方法可以得到精確的初至走時(shí)和準(zhǔn)確的層析核函數(shù)，進(jìn)而保證層析方程組的準(zhǔn)確性，得到更好的層析反演速度場(chǎng)。PTISPR方法在起伏地表?xiàng)l件下可以正確地追蹤初至波射線路徑，是初至波走時(shí)層析速度建模的核心技術(shù)。理論模型和實(shí)際資料試算結(jié)果證明，本文提出的PTISPR法是正確有效的，適用于起伏地表?xiàng)l件下的初至走時(shí)層析建模和靜校正。

參 考 文 獻(xiàn)

[1] 劉玉柱,程玖兵,董良國(guó).面向起伏地表偏移成像的表層靜校正方法[J].石油物探,2012,51(6):584-589

Liu Y Z,Cheng J B,Dong L G.A new static correction method for the migration from rugged topography[J].Geophysical Prospecting for Petroleum,2012,51(6):584-589

[2] 宋桂橋，于世煥.山前帶地震勘探技術(shù)進(jìn)展與對(duì)策研究[J].石油物探,2012,51(6):539-547

Song G Q,Yu S H.Progress and strategy of the seismic exploration in foothill area[J].Geophysical Prospecting for Petroleum,2012,51(6):539-547

[3] 岳玉波，孫建國(guó)，楊昊，等.起伏地表下初值射線追蹤的實(shí)現(xiàn)[J].勘探地球物理進(jìn)展，2007,30(5):388-391

Yue Y B,Sun J G,Yang H，et al.Realization of initial value ray tracing for rugged topography[J].Progress in Exploration Geophysics,2007,30(5):388-391

[4] Cassel B R.A method for calculating synthetic seismograms in laterally varying media[J].Geophysical Journal International,1982,69(2),339-354

[5] Langan R T,Lerche I,Cutler R T.Tracing of rays through heterogeneous media:an accurate and efficient procedure[J].Geophysics,1985,50(9),1456-

1465

[6] Moser T J.Shortest path calculation of seismic rays[J].Geophysics,1991,56(1):59-67

[7] 張建中，陳世軍，許初偉.動(dòng)態(tài)網(wǎng)絡(luò)最短路徑射線追蹤[J].地球物理學(xué)報(bào)，2004,47(5)：889-904

Zhang J Z，Chen S J，Xu C W.A method of shortest path raytracing with dynamic networks[J].Chinese J Geophys(in chinese),2004,47(5):889-904

[8] 桑運(yùn)云，李振春，張凱.基于拋物旅行時(shí)插值的最短路徑射線追蹤[J].石油地球物理勘探,2013,48(3):403-409

Sang Y Y,Li Z C,Zhang K.Shortest path raytracing based on parabolic traveltime interpolation[J].Oil Geophysical Prospecting,2013,48(3):403-409

[9] Dijkstra E W.A note on two problems in connection with graphs [J].Numerische Mathematik,1959,1(1):269-271

[10] Stork C,Clayton R W.Linear aspects of tomographic velocity analysis [J].Geophysics，1991，56(4):483-495

[11] 韓曉麗，楊長(zhǎng)春，麻三懷.復(fù)雜山區(qū)初至波層析反演靜校正[J].地球物理學(xué)進(jìn)展,2008,23(2)：475-483

Han X L,Yang C C,Ma S H.Static of tomographic inversion by first break in complex areas [J].Progress in Geophysics,2008,23(2):475-483