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

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

隨著地震勘探技術(shù)的發(fā)展和勘探規(guī)模的擴(kuò)大，起伏地表?xiàng)l件下的地震勘探越來越普遍[1-2],通過有效的近地表速度建模方法解決靜校正問題成為地震資料處理的重點(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)：試射法通過迭代修改初始入射角獲得了準(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ù)來描述起伏地表以及地下界面,然后給定射線的初始點(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í)現(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ò)展，即從震源開始計(jì)算整個(gè)速度場(chǎng)的所有節(jié)點(diǎn)的最小走時(shí);②反向追蹤路徑，即從接收點(diǎn)開始向震源反向追蹤射線路徑。其核心思想是：首先利用Dijkstra[9]算法計(jì)算某節(jié)點(diǎn)旅行時(shí)，然后在計(jì)算過程中搜索滿足拋物插值的節(jié)點(diǎn)并進(jìn)行插值，將最小旅行時(shí)作為該計(jì)算節(jié)點(diǎn)的最終時(shí)間。射線追蹤過程是利用拋物旅行時(shí)插值從接收點(diǎn)網(wǎng)格開始一直追蹤到炮點(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è)過程，向前計(jì)算每個(gè)節(jié)點(diǎn)的旅行時(shí)[6]，向后追蹤射線路徑。將起伏地表網(wǎng)格化，如圖1所示，地表上速度為0，地表面以上速度區(qū)域?yàn)闊o效區(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)所屬的無效節(jié)點(diǎn)數(shù)(即A節(jié)點(diǎn)及其以上采樣點(diǎn)對(duì)應(yīng)的邊界節(jié)點(diǎn)和角節(jié)點(diǎn)無效節(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)橫向位置，則為無效節(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í)，為無效節(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)為無效節(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í)，說明此計(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è)為無效節(jié)點(diǎn);③有兩個(gè)節(jié)點(diǎn)為無效節(jié)點(diǎn)，僅有一個(gè)節(jié)點(diǎn)時(shí)間已知。因此，在起伏地表射線追蹤過程中，需要對(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)開始向炮點(diǎn)位置反向追蹤的過程，當(dāng)追蹤到炮點(diǎn)所處的網(wǎng)格時(shí)，需要判斷當(dāng)前接收點(diǎn)和炮點(diǎn)是否在同一個(gè)網(wǎng)格范圍內(nèi)。由于起伏地表的存在，射線路徑在同一網(wǎng)格內(nèi)經(jīng)過的網(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í)線，可見起伏地表連續(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)比較可見，近地表低頻信息得到了很好的反演，基本上反映了低、降速帶內(nèi)速度變化的趨勢(shì)。圖9為真實(shí)速度和層析速度的等時(shí)線和射線路徑，可見兩者基本吻合。經(jīng)過模型試算，證明起伏地表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可見，通過對(duì)初始模型(圖10a)進(jìn)行迭代處理，得到的層析建模速度場(chǎng)(圖10b)從南到北高程變化劇烈，呈現(xiàn)南低北高的趨勢(shì)，北側(cè)低、降速帶范圍埋深更淺，符合實(shí)際地形特征。利用層析建模得到的速度場(chǎng)計(jì)算靜校正量并對(duì)原始炮集進(jìn)行靜校正[11]，圖11為某單炮記錄初至層析靜校正前、后的結(jié)果。從圖11可見，靜校正前(圖11a)同相軸偏離雙曲規(guī)律，抖動(dòng)比較劇烈，初至同相軸不連續(xù);經(jīng)過靜校正后(圖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é)束語

射線追蹤是走時(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] 劉玉柱,程玖兵,董良國.面向起伏地表偏移成像的表層靜校正方法[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] 岳玉波，孫建國，楊昊，等.起伏地表下初值射線追蹤的實(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