以共反射點為目標(biāo)的二維射線追蹤方法研究

劉軍勝

(中石化地球物理有限公司,北京 100000)

0 引言

地震資料采集的過程中,目的層的有效覆蓋次數(shù)是影響資料信噪比和成像效果的關(guān)鍵因素,資料信噪比與覆蓋次數(shù)的平方根成正比。技術(shù)人員一般利用照明分析,射線追蹤模擬方法解決地下復(fù)雜構(gòu)造形態(tài)下覆蓋次數(shù)均勻性問題,而其中射線追蹤方法是最普遍的方法[1-4]。射線追蹤自上世紀(jì)70 年代首次提出至今,以不同的數(shù)值計算方法為關(guān)鍵技術(shù),以計算描述地震波運動學(xué)特征的相關(guān)參數(shù)(走時、射線路徑)等為核心目標(biāo),相繼提出了十多種方法[5]。當(dāng)前最常用的射線追蹤的方法為試射法、迭代法、波前追蹤法,地震勘探采集設(shè)計軟件的射線追蹤模塊主要基于以上三種方法。KLSEIS軟件的二維射線追蹤模塊綜合運用了試射法和迭代法,MESA軟件運用了波前追蹤法,OMNI則基于試射法。三種方法計算精度和計算速度存在較大的差異[6-9]。例如試射法的計算精度主要取決于設(shè)定的入射角度范圍和增量,迭代法的計算精度則受計算網(wǎng)格大小的影響。復(fù)雜構(gòu)造形態(tài)條件下,若參數(shù)設(shè)置不合理,各種方法均存在遺漏路徑的問題。計算精度較高的迭代法在穿越相同炮檢點存在多路徑時,僅能計算一條路徑,同時無法限定反射和入射角范圍。由于漏解、精度的問題,常用的設(shè)計軟件射線追蹤模塊均缺少覆蓋次數(shù)計算功能。更為關(guān)鍵的是,常用的射線追蹤方法直接計算地震波自激發(fā)點到接收點的路徑,限制了研究人員查找和判斷影響目的層關(guān)鍵區(qū)域的炮集或道集。設(shè)計人員利用設(shè)計軟件射線追蹤模塊僅能查看目的層區(qū)域射線的稀疏和稠密程度,發(fā)現(xiàn)問題,無法給出明確的解決方案。設(shè)計人員需要借助其他方法或是通過不斷的調(diào)試來查找影響問題區(qū)域的炮檢點進(jìn)而解決問題。專業(yè)模型軟件Tesseral Pro在這一方面有所改善,可以查看目的層選定段的所有射線路徑。但Tesseral Pro僅能按照流程布設(shè)理論觀測形式,對野外生產(chǎn)中的特殊觀測系統(tǒng)設(shè)計優(yōu)化指導(dǎo)意義較小。筆者改變了一般的射線追蹤思路,將上覆地層界面進(jìn)行區(qū)間化,以目的層共反射點為目標(biāo),分段計算地震波傳播路徑,利用Matlab軟件編程實現(xiàn),克服了三種主流射線追蹤方法在路徑計算方面的缺點,準(zhǔn)確確定目的層覆蓋次數(shù)分布特征,查找影響覆蓋次數(shù)的關(guān)鍵炮集,優(yōu)化激發(fā)點,彌補常用設(shè)計軟件的不足,指導(dǎo)野外生產(chǎn),提升資料品質(zhì)。

1 方法原理

1.1 基本思路

假設(shè)地層為均勻連續(xù)層狀介質(zhì),則射線由起始點向下穿越不同的目的層經(jīng)過反射,折射后到達(dá)目的層固定反射點位置,然后再經(jīng)過不同地層反射、折射后回到地面的路徑是可以精確計算的。當(dāng)以激發(fā)點為起始時,在限定偏移距的條件下,可以計算出經(jīng)過每個反射點的路徑數(shù)量,得出目的層的有效覆蓋次數(shù)。當(dāng)以接收點為起始時,可以計算出實現(xiàn)目的層覆蓋次數(shù)均勻條件下的所有炮點分布特征。無需限定終止點的位置,射線與地面的交點若未處于整道位置,通過微調(diào)反射點的位置即能實現(xiàn)交點歸位,而且微調(diào)后的反射點與原反射點仍屬于同一線元。

1.2 射線路徑計算方法

1.2.1 單一路徑的計算

以共反射點位目標(biāo)的射線追蹤方法將射線路徑分成由起始點到反射點和由反射點到終止點兩部分。將上覆地層細(xì)分成若區(qū)間,地震波的穿越不同地層每一段區(qū)間的傳播路徑間的相交,在地層界面上的反射和折射關(guān)系表示為多個非線性方程,計算反射、透射點的過程也就轉(zhuǎn)換成非線性方程組求根的過程。Matlab中的多元非線性方程組快速求解算法、微分方法、插值算法使射線路徑的計算得以快速實現(xiàn)。

圖1為簡易三層界面地震模型,v1、v2、v3為不同地層的層速度;L11-L12,L21-L22,L31-L32,L23-L24,L13-L14為不同界面上有限長度的線段區(qū)間。假設(shè)存在S到C然后到R的路徑,經(jīng)過L11-L12之間的P1,L21-L22之間的P2,L23-L24之間的P3,L13-L14之間的P4,那么射線路徑計算方法如下:

圖1 炮點經(jīng)CRP點到達(dá)地表的路徑示意圖Fig.1 Source-CRP-surface pathway diagram

1)S-C路徑計算

y1-y0=k11*(x1-x0)

(1)

y1-y12=k1*(x1-x12)

(2)

y2-y1=k22*(x2-x1)

(3)

y2-y22=k2*(x2-x22)

(4)

(5)

cos(tan-1k11-tan-1k1)*v2=

cos(tan-1k22-tan-1k1)*v1

(6)

cos(tan-1k22-tan-1k2)*v3=

cos(tan-1k33-tan-1k2)*v2

(7)

Matlab自帶的多元非線性方程組求解函數(shù),每次僅能計算出初始賦值附近的一個根。求得方程組根以后,需要確保P1(x1,y1)在L11-L12之間,P2(x2,y2)在L21-L22之間,以及入射路徑和折射路徑位于法線的兩側(cè)。針對射線路徑的判斷,運用向量的夾角計算公式進(jìn)行推導(dǎo)計算,方程根需滿足式8和式9。只有滿足以上所有條件的根才算方程組的有效解。

(8)

(9)

2)C-R路徑計算

在計算得到P2以后,可計算P2C的斜率,進(jìn)而得出P3C的斜率。利用線性插值和微分方法計算出L23、L24、L13、L14的坐標(biāo),L23L24,L13L14的斜率k5,k6。計算出由反射點C到P4的路徑需要建立基于P3P4斜率,P3,P4坐標(biāo)的5元非線性方程組,包含過P3、P4的4個線段的斜率,滿足P3處斯奈爾定律的5個方程,與S-C路徑的方法類似。滿足P3在L23-L24之間,P4在L13-L14之間,以及入射路徑和折射路徑位于法線兩側(cè)解方程的解為有效解。當(dāng)?shù)貙訑?shù)量為n時,方程組的元數(shù)為3n-1。得到交點P3,P4以后,進(jìn)而計算出R的位置。

1.2.2 全路徑的計算

多元非線性方程具有多解性的特點,但在一定的范圍內(nèi),解是唯一或是不存在的。對于已知的觀測系統(tǒng)和地質(zhì)模型,計算過目的層所有線元的射線,需要將上覆地層界面以固定的區(qū)間進(jìn)行劃分,利用S-C,C-R的路徑計算方法和循環(huán)運算,逐個判斷和查找穿越所有區(qū)間的路徑。假設(shè)地表物理點的個數(shù)為m,目的層CRP點為n,每層上覆地層界面數(shù)量為i,每層區(qū)間的數(shù)量為k,則循環(huán)運算的次數(shù)c滿足式5。

c=m*n*k2*i

(10)

此方法不但可以計算S-C-R的路徑,也可以將S替換成R計算R-C-S的路徑。當(dāng)所有的路徑計算完成以后,根據(jù)觀測系統(tǒng)設(shè)定限制偏移距和反射角大小后,進(jìn)而統(tǒng)計穿越不同線元的路徑數(shù)量可以得到目的層的有效覆蓋次數(shù)。若不限定偏移距范圍,通過S-C-R的路徑計算可以得出在炮點均勻布設(shè)條件下實現(xiàn)較均勻采樣的道集,進(jìn)而指導(dǎo)偏移距的選擇。在限定偏移距和反射角的條件下,也可以通過R-C-S逆向路徑的計算,得到檢波點均勻布設(shè)條件下,實現(xiàn)較均勻采樣的炮集,指導(dǎo)激發(fā)點優(yōu)化。

1.2.3 關(guān)鍵參數(shù)的設(shè)定

在射線路徑的計算中,穿越目的層固定線元的最終路徑的多少除受地層特征和觀測系統(tǒng)設(shè)定影響外,還受到反射角限定和上覆地層劃分區(qū)間大小的影響。

地層入射角方面,華凱等人[10]在入射角對地震動特性研究中指出隨著入射角增的增大,P波作用引起的水平位移動與S波作用引起的垂直位移通常會增大,地震動軌跡構(gòu)成的輪廓面積隨之增大,持續(xù)時間變小。當(dāng)入射角達(dá)到臨界角以后,進(jìn)而產(chǎn)生折射波等。區(qū)別于其他方法,本方法中研究人員可以根據(jù)研究需要自主限定反射角大小,滿足不同的研究條件。

地層區(qū)間劃分方面,類似于迭代法的網(wǎng)格設(shè)定,上覆地層區(qū)間設(shè)定是影響計算精度和效率的關(guān)鍵參數(shù)。地下地層受構(gòu)造變動影響往往為彎曲的界面,但地層界面在有限的范圍內(nèi)可以近似表示為直線,斜率可以準(zhǔn)確計算。范圍越小,利用多段線段刻畫出來的界面就越準(zhǔn)確。隨著范圍增大,計算精度逐漸發(fā)生變化。同時由于MATLAB每次僅能計算出方程組的一個解,當(dāng)同一區(qū)間存在多路徑時會出現(xiàn)漏解的問題。研究人員對比分析了將目的層的上覆地層以線元(1/2道距),1/2線元,2倍線元劃分,以及有限折點法(將地層按照轉(zhuǎn)折點分成少量的線段)等區(qū)間劃分方式,除有限折點法外,其余三種方法得出的結(jié)果一致,研究人員認(rèn)為在線元或者線元大小相近的尺度上劃分區(qū)間大小能夠確保路徑計算的完整性,避免漏解。

2 模擬驗證及應(yīng)用分析

2.1 模擬驗證

建立圖2所示的二維三層地質(zhì)模型,層速度分別為2 500 m/s,2 800 m/s,3 200 m/s,理論設(shè)計方案道距為50 m,炮點距為100 m,觀測系統(tǒng)為2 475-25-50-25-2 475。在2 475 m～7 575 m范圍內(nèi)共布設(shè)52炮。分別利用Tesseral Pro專業(yè)模型軟件和基于共反射點的射線追蹤方法進(jìn)行對比分析。圖3為Tesseral Pro 射線追蹤模擬結(jié)果,可以看出在凹陷的兩翼存在明顯的“反射盲區(qū)”,無有效反射路徑。圖4為利用Tesseral Pro進(jìn)行波動方程模擬放炮、疊加、偏移、時深轉(zhuǎn)換后的深度剖面,從結(jié)果中顯示在凹陷的兩翼存在弱反射能量,說明這些照明強度不足,同時也說明圖3顯示的“反射盲區(qū)”并非真正的盲區(qū)。利用新的射線追蹤方法(本次限定反射角小于臨界角的前提,實際可根據(jù)研究需要進(jìn)行調(diào)整),結(jié)果見圖5和圖6。圖5顯示給定的地震勘探模型無明顯的盲區(qū)。圖6為圖5的局部放大,顯示在凹陷兩翼的射線路徑明顯減少,且存在少量線元無射線路徑?；诠卜瓷潼c的射線追蹤方法結(jié)論與波動方程模擬結(jié)果相近。

圖2 二維模型Fig.2 2D model

圖3 Tesseral Pro 射線追蹤結(jié)果Fig.3 Ray tracing result of Tesseral Pro

圖4 Tesseral Pro波動方程模擬深度剖面Fig.4 Depth section of wave equation simulation with Tesseral Pro

圖5 新射線追蹤方法路徑計算結(jié)果Fig.5 Ray tracing result of new method

圖6 新射線追蹤方法路徑結(jié)果(右側(cè)凹陷左翼放大)Fig.6 Ray tracing result of new method(left wing of right sunken zooming in)

利用新方法對射線路徑進(jìn)行統(tǒng)計,得出目的層的有效覆蓋次數(shù),見圖7。通過覆蓋次數(shù)顯示可以看出第二層左凹陷的兩翼和右凹陷的左翼存在覆蓋次數(shù)明顯低于周邊地層的問題,特別是右側(cè)凹陷的左翼在7 287.5處覆蓋次數(shù)為零,該運算結(jié)果與波動方程模擬結(jié)果基本一致,進(jìn)一步驗證了方法的準(zhǔn)確性。

圖7 目的層有效覆蓋次數(shù)Fig.7 Effective folds of target layer(藍(lán)線:第一層 紅線:第二層)

2.2 應(yīng)用分析

以共反射點為目標(biāo)的二維射線追蹤方法可以計算目的層的有效覆蓋次數(shù),同時可以計算基于共反射點的道集、炮集記錄。實際應(yīng)用中可以通過道集、炮集分析優(yōu)化觀測系統(tǒng)或激發(fā)點,實現(xiàn)目的層的較均勻采樣。特別是針對復(fù)雜構(gòu)造形態(tài),合理的優(yōu)化激發(fā)點,能夠在保障目的層資料品質(zhì)的前提下,最大限度地降低生產(chǎn)成本。

筆者將重點以利用該方法進(jìn)行二維激發(fā)點優(yōu)化來分析方法的應(yīng)用效果?；緦崿F(xiàn)思路為通過S-C-R的路徑計算反射點覆蓋次數(shù)(正演),查找出問題區(qū)域、利用R-C-S的路徑計算實現(xiàn)問題區(qū)域較均勻覆蓋的炮集(反演)、優(yōu)化激發(fā)點然后進(jìn)行正演驗證激發(fā)點優(yōu)化后正演驗證,具體實現(xiàn)流程見圖8,總體技術(shù)路線和實現(xiàn)共分10個步驟、多次循環(huán)。

圖8 激發(fā)點優(yōu)化技術(shù)路線和流程Fig.8 Shots optimization technical route and flow

1)正演模擬

正演模擬主要計算在理論設(shè)計方案下,主要目的層的有效覆蓋次數(shù),從而界定需要改善的區(qū)域,改善區(qū)域主要分為三類:①低覆蓋區(qū)域,需要進(jìn)行加密激發(fā);②高覆蓋區(qū)域,需要進(jìn)行稀疏激發(fā);③覆蓋適中區(qū)域,可根據(jù)需要加密或稀疏激發(fā)。

根據(jù)圖7可以清晰看出滿覆蓋區(qū)域目的層反射界面的有效覆蓋次數(shù)在25次左右。形態(tài)復(fù)雜的第二界面線元4 562.5～4 937.5,6 812.5～7 287.5處屬于低覆蓋區(qū)域;線元7 812.5～8 637.5正常應(yīng)處于覆蓋次數(shù)漸減帶,但覆蓋次數(shù)較高,兩段需要通過激發(fā)點優(yōu)化改善。

2)反演模擬

利用程序計算出對應(yīng)的需要加密、稀疏或保持不變的炮集區(qū)域(計算的過程中限定了激發(fā)點的范圍為2 025～8 075)。結(jié)果見圖9,2 025～4 275段、6 875～8 075段的為低覆蓋CRP炮集段,通過該區(qū)域的激發(fā)點加密能夠提升有效覆蓋次數(shù)。5 175～8 075段則為高覆蓋CRP炮集段,通過對該區(qū)域的稀疏采樣能夠降低有效覆蓋次數(shù)。

圖9 實現(xiàn)不同目的層區(qū)域均勻覆蓋的炮集分布Fig.9 Shot gather to achieve the uniform coverage of different target layer section

3)優(yōu)化激發(fā)點

對反演模擬的結(jié)果進(jìn)行綜合分析,對所有激發(fā)點進(jìn)行重新規(guī)劃,在保證激發(fā)點總數(shù)變化較小的基礎(chǔ)上,確保目的層有效覆蓋次數(shù)均衡。針對圖3所示的模型,在綜合考慮目的層成像效果和勘探成本的基礎(chǔ)上,確定在2 425～4 175段采用3道2炮的加密激發(fā),4175～6875段采用3道1炮的稀疏激發(fā),6 875～7 925段進(jìn)行1道1炮加密激發(fā)的激發(fā)點優(yōu)化方案,總體設(shè)計58炮,較原方案僅增加6炮。但關(guān)鍵低覆蓋區(qū)域的覆蓋次數(shù)得到較大的改善,特別是線元7 287.5處的覆蓋次數(shù)由0次提升到了7次,見圖10。(a:二維模型;b:第一層有效覆蓋次數(shù)變化 c:第二層有效覆蓋次數(shù)變化)

圖10 激發(fā)點優(yōu)化后不同目的層有效覆蓋次數(shù)圖Fig.10 Effective folds diagram of different target layers before and after shots optimization

4)效果分析

以共反射點為目標(biāo)的射線追蹤思路,使該方法在提升目的層成像效果方面更有針對性,也更為便捷。基于本方法的整套Matlab程序能夠快速、直接地查找二維地震勘探過程中目地層的低覆蓋區(qū)域,明確影響覆蓋次數(shù)的關(guān)鍵炮點和檢波點,實現(xiàn)物理點的精準(zhǔn)優(yōu)化,以最小的投入提升資料的品質(zhì)。

3 結(jié)論

筆者改變了由激發(fā)點到接收點的射線追蹤思路,提出了由激發(fā)點或接收點經(jīng)過CRP點再到地表的射線路徑計算方法,并利用軟件編程實現(xiàn)自動化計算,用于指導(dǎo)和驗證二維激發(fā)點和觀測系統(tǒng)優(yōu)化及效果分析。但本方法主要存在兩方面的不足:

1)研究主要針對P波入射,P波反射和透射計算,未計算轉(zhuǎn)換波、折射波等衍生地震波,解決有效覆蓋的問題而無法進(jìn)行波組分析。

2)本文中的方法編程實現(xiàn)利用了MATLAB在矩陣運算,方程求解等方面的優(yōu)勢,降低了編程難度。但MATLAB由于其自身基于向量的運算方式,在循環(huán)運算方面存在效率低的短板。針對圖2所示的模型,完成圖8的全流程計算需要1個小時。

筆者所述的方法也可以通過其他編程語言實現(xiàn)或是在并行計算方面進(jìn)行改進(jìn),進(jìn)一步提升運算效率。此方法可作為傳統(tǒng)射線追蹤方法或是波動方程模擬方法的補充,快速解決重點區(qū)域構(gòu)造成像問題。在實施過程中為了設(shè)計工作者可以僅主要針對重點目的層段進(jìn)行分析和優(yōu)化,或者對形態(tài)單一的目的層的上覆地層進(jìn)行簡化處理,降低計算工作量,使運算更為便捷。