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

劉軍勝

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

0 引言

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

1 方法原理

1.1 基本思路

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

1.2 射線路徑計(jì)算方法

1.2.1 單一路徑的計(jì)算

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

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

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

1)S-C路徑計(jì)算

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

(8)

(9)

2)C-R路徑計(jì)算

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

1.2.2 全路徑的計(jì)算

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

c=m*n*k2*i

(10)

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

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

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

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

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

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

2.1 模擬驗(yàn)證

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

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

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

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

圖5 新射線追蹤方法路徑計(jì)算結(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)

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

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

2.2 應(yīng)用分析

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

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

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

1)正演模擬

正演模擬主要計(jì)算在理論設(shè)計(jì)方案下,主要目的層的有效覆蓋次數(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ā)點(diǎn)優(yōu)化改善。

2)反演模擬

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

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

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

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

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

4)效果分析

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

3 結(jié)論

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

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

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

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