以共反射點為目標的二維射線追蹤方法研究

劉軍勝

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

0 引言

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

1 方法原理

1.1 基本思路

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

1.2 射線路徑計算方法

1.2.1 單一路徑的計算

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

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

圖1 炮點經CRP點到達地表的路徑示意圖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自帶的多元非線性方程組求解函數,每次僅能計算出初始賦值附近的一個根。求得方程組根以后,需要確保P1(x1,y1)在L11-L12之間,P2(x2,y2)在L21-L22之間,以及入射路徑和折射路徑位于法線的兩側。針對射線路徑的判斷,運用向量的夾角計算公式進行推導計算,方程根需滿足式8和式9。只有滿足以上所有條件的根才算方程組的有效解。

(8)

(9)

2)C-R路徑計算

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

1.2.2 全路徑的計算

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

c=m*n*k2*i

(10)

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

1.2.3 關鍵參數的設定

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

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

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

2 模擬驗證及應用分析

2.1 模擬驗證

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

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

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

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

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

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

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

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

2.2 應用分析

以共反射點為目標的二維射線追蹤方法可以計算目的層的有效覆蓋次數,同時可以計算基于共反射點的道集、炮集記錄。實際應用中可以通過道集、炮集分析優化觀測系統或激發點,實現目的層的較均勻采樣。特別是針對復雜構造形態,合理的優化激發點,能夠在保障目的層資料品質的前提下,最大限度地降低生產成本。

筆者將重點以利用該方法進行二維激發點優化來分析方法的應用效果。基本實現思路為通過S-C-R的路徑計算反射點覆蓋次數(正演),查找出問題區域、利用R-C-S的路徑計算實現問題區域較均勻覆蓋的炮集(反演)、優化激發點然后進行正演驗證激發點優化后正演驗證,具體實現流程見圖8,總體技術路線和實現共分10個步驟、多次循環。

圖8 激發點優化技術路線和流程Fig.8 Shots optimization technical route and flow

1)正演模擬

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

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

2)反演模擬

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

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

3)優化激發點

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

圖10 激發點優化后不同目的層有效覆蓋次數圖Fig.10 Effective folds diagram of different target layers before and after shots optimization

4)效果分析

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

3 結論

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

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

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

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