基于波場上傳試射的三維射線追蹤方法

辛全圣，賈 雨，舒 引

(1.成都理工大學核技術與自動化工程學院，成都 610059;2.電子科技大學通信與信息工程學院，成都 611731)

引 言

射線追蹤方法在地震勘探中有廣泛應用，它在地震定位、地球物理層析成像、地震勘探資料的偏移處理、地震模擬及其它地震數(shù)字處理與反演中都起著極重要的作用。傳統(tǒng)的射線追蹤方法有試射法［1］和彎曲法［2］，后來發(fā)展起來的有，有限差分法［3-4］、旅行時線性插值法［5-7］、最短路徑法［8-9］、波前重建方法［10］等。與其他方法相比，試射法精確性高，在全面搜索和適應復雜地質(zhì)模型等方面有很大優(yōu)勢，是兩點追蹤問題的一種重要方法［11］。

試射法，又稱打靶法，是基于斯奈爾(Snell)定理的計算方法。試射法是最早提出和使用的射線追蹤方法，在數(shù)學上屬于初值問題。試射法根據(jù)由震源出發(fā)的一束射線到達接收點的情況對射線出射角及其密度進行調(diào)整，最后由最靠近接收點的兩條射線走時內(nèi)插求出接收點走時［12］，從而實現(xiàn)兩點間的射線追蹤計算。試射法雖然是一種經(jīng)典的射線追蹤方法，但不能誤認為它是低級的、行將淘汰的算法，相反，一直在不斷的發(fā)展。試射法在部分問題上仍需進一步改進:其一是在計算精度方面存在難以解決每一個炮檢對存在多條射線的問題［其二是在計算效率方面存在每一次觀測系統(tǒng)變化都需要進行試射和迭代的完整流程而導致計算效率較慢的問題。為此，本文對CRP統(tǒng)計數(shù)據(jù)的準確度和三維射線追蹤的計算效率進行改進。

1 方法原理

通過試射獲得可能存在射線的區(qū)域，通過迭代計算準確的射線軌跡。首先以目標反射層為基準，計算反射層的每一個控制點到地表的全方位試射射線，并作為迭代第一步的初始數(shù)據(jù)，只要模型不變化，則試射結果一直有效;然后給定觀測系統(tǒng)，對控制點到每一個炮點的射線進行迭代，并記錄有效射線作為迭代第二步的初始數(shù)據(jù)，只要觀測系統(tǒng)的炮位置不變化，則炮點到目標層控制點計算結果一直有效;最后遍歷每個目標層三角形，迭代出該反射點在該三角形內(nèi)部的所有射線。

1.1 試射

傳統(tǒng)的試射方法主要采用從炮點進行試射，水平角度的范圍為0到360度，垂直角度范圍為0到180度，往反射層試射。到達反射層之前通過和子面求交點，然后發(fā)生折射，如果折射不成功，則這個角度試射失敗，用下一角度繼續(xù)進行試射。當射線到達反射層時，根據(jù)反射定理，利用當前方向和反射點的法向求得反射方向，從反射層往地表計算射線，簡單的二維示意圖如圖1所示。本文采用的試射方法對傳統(tǒng)的試射方法做了相應改進，本文采用的試射方法從反射層的反射點進行試射，水平角度的范圍為0到180度，垂直角度范圍為0到90度。簡單的二維示意圖如圖2所示。

圖1 從炮點進行試射二維示意圖

圖2 從反射點進行試射二維示意圖

1.2 迭代

迭代部分是把試射部分的結果作為原始數(shù)據(jù)，通過三角形面積迭代法最終找出炮點到檢波點的有效射線。先根據(jù)試射文件求出每一個炮點到反射層每一個頂點經(jīng)反射后再到接收層的射線，記錄文件，如圖3所示。然后根據(jù)上一步求出的射線找出炮點經(jīng)反射層反射到其關聯(lián)檢波器的射線，即為射線追蹤最終射線結果，如圖4所示。

1.3 三角形面積迭代法

圖3 迭代第一步示意圖

圖4 迭代第二步示意圖

如圖5所示，已知目標點P和所在的出射三角形T1T2T3，面積為S，P將T1T2T3分成三個三角形，面積分別為S1，S2，S3。三個頂點對應的射線的方向向量分別為點在三角形中的面積坐標加權求和得出新的方向分量以該方向分量試射并得到出射點T，如果出射點T和接收器的距離小于給定的精度，則試射迭代結束。否則，用出射點T把原來的三角形分割成三個新的三角形 TT1T2、TT2T3和TT1T3。找到P點所在三角形，并作為下次試射迭代的初始三角形。按上述方法繼續(xù)試射直到出射點和接收器的距離小于給定的精度為止。

圖5 三角形面積迭代法示意圖

2 基于波場上傳試射算法

2.1 試射

試射是射線追蹤核心算法的第一步，試射的結果以射線的格式存放文件，一條射線記錄了該射線的方向、反射點坐標、射線路徑中和每一個層面的交點的坐標和法向、以及接收層的入射點和出射點。試射的第一步是確定試射角度。理論上試射射線在接收層越密集，即角度間隔越小，射線追蹤的結果越準確，但是角度間隔越小，試射部分所花的時間就越多，可根據(jù)工區(qū)的實際情況選擇合適的角度間隔。

2.2 迭代第一步

(1)首先初始化觀測系統(tǒng)，讀取所有炮點和檢波點，獲得炮點和檢波點的x值和y值，通過向接收層投影，給炮點和檢波點賦z值。

(2)逐一處理每一個反射點，首先讀取反射點的射線文件，根據(jù)每條射線的出射點x，y坐標，將該射線投入網(wǎng)格中，利用一個二維的vector變量存放對應射線的索引。

(3)遍歷處理每一個炮點，首先根據(jù)炮點的x，y坐標找到炮點屬于的網(wǎng)格，然后根據(jù)二維的vector變量數(shù)據(jù)，取出該網(wǎng)格以及周圍一圈網(wǎng)格中的射線索引，通過這些射線索引獲取每一條射線，從這些射線中選出三條出射點包圍該炮點的射線，如果存在多條，則選出構成的三角形面積最小的三條射線。如圖3所示，三條紫色射線r1，r2，r3為選出的三條滿足要求的射線。

(4)根據(jù)找出的三條射線和炮點利用三角形面積迭代法找出炮點到反射點的射線r，如果迭代過程中不收斂，即某一個方向的試射射線的出射點不在迭代的三角形中，則返回失敗。

2.3 迭代第二步

(1)對每一個炮點逐一處理，讀取該炮點到反射點經(jīng)反射后再到接收層的所有射線。

(2)遍歷反射層每一個三角形，取出三角形三個反射點，判斷當前炮點到三個反射點是否都有射線，如果到某個點沒有射線，則跳過當前三角形，處理下一個三角形。

(3)如果三角形三個頂點都有射線，則取出三條射線。如圖4所示，r1，r2，r3為炮點到反射層當前三角形三個頂點的三條射線，三條射線在接收層的出射點構成的三角形為T1T2T3。遍歷當前炮點的關聯(lián)檢波點，如果檢波點在三角形內(nèi)部，則用三角形面積迭代法找出最終射線r。需要注意的是，這里三角形迭代法找出最終射線的比例比迭代第一步找出射線的比例低，因為第一步中是利用三角形面積迭代法找出反射點到炮點的射線，而這一步是找出從炮點到檢波點經(jīng)過反射層反射的射線，隨著折射反射越多，收斂比例越低。

3 數(shù)值模擬及分析

本文通過對復雜地質(zhì)構造理論模型進行射線路徑的計算，驗證該算法的有效性。圖6為射線路徑圖，炮點為2000個，25條線，每條線80個炮點，每炮關聯(lián)檢波器為1152個，紅色排列片代表接收的檢波點，綠色的點代表未接收的檢波點。圖7為水平界面的理論CMP面元情況與實際射線追蹤統(tǒng)計的CRP面元情況對比圖。可以看出，當界面為水平界面時，CRP方位角和炮檢距統(tǒng)計圖與基于CMP理論計算的完全一致，且不同偏移距范圍的覆蓋次數(shù)呈規(guī)則分布。CRP位于炮檢點的中點，與理論CMP重合，射線路徑對稱，CRP覆蓋次數(shù)分布規(guī)律與理論計算的覆蓋次數(shù)分布規(guī)律完全一致?？梢则炞CCRP統(tǒng)計數(shù)據(jù)的準確度明顯提高。以馮家灣三維工區(qū)為例，試射單線程單機計算估算15天左右，迭代第一步計算量為試射的1/5左右，迭代第二步相對于前兩步計算量很小，估計少一個數(shù)量級。計算效率的提高取決于實際的工區(qū)和觀測系統(tǒng)，數(shù)據(jù)量越大、觀測系統(tǒng)變化的次數(shù)越多，計算效率提高越明顯。傳統(tǒng)的試射方法每一次觀測系統(tǒng)變化都需要進行試射和迭代的完整流程。本文的試射方法只要模型不變化，則試射結果一直有效;只要觀測系統(tǒng)的炮位置不變化，則炮點到目標層控制點計算結果一直有效，提高了三維射線追蹤的計算效率。

圖6 射線路徑圖

圖7 水平界面的理論CMP面元情況與實際射線追蹤統(tǒng)計的CRP面元情況對比圖

4 結束語

傳統(tǒng)的試射方法從炮點出發(fā)，通過調(diào)整入射角使得射線逐步逼近檢波點位置，有些檢波點無法接收到射線，存在盲區(qū)。針對傳統(tǒng)的試射方法難以解決每一個炮檢對存在多條射線以及每一次觀測系統(tǒng)變化都需要進行試射和迭代的完整流程而導致計算效率較慢的問題，本文提出了基于波場上傳試射的三維射線追蹤方法，從原理上解決了每一個炮檢對存在多條射線的問題。同時，只要模型不變化，則試射結果一直有效;只要觀測系統(tǒng)的炮位置不變化，則炮點到目標層控制點計算結果一直有效。數(shù)值模擬結果表明，該算法能夠提高CRP統(tǒng)計數(shù)據(jù)的準確度，根據(jù)實際工區(qū)和觀測系統(tǒng)的不同，計算效率有不同程度的提高，適于復雜地質(zhì)構造三維射線追蹤的工程應用。

［1］馬爭鳴，李衍達.二步法射線追蹤［J］.地球物理學報，1991，34(4):501-508.

［2］高爾根，徐果明，蔣先藝，等.三維結構下逐段迭代射線追蹤方法［J］.石油地球物理勘探，2002，37(1):11-16.

［3］ Vidale J E.Finite-difference calculation of traveltimes［J］.Bulletin of the seismological society of A-merica，1988，78(6):2067-2076.

［4］ Vidale J E.Finite-difference calculation of traveltimes in three dimension［J］.Geophysics，1990，55:521-526.

［5］張東，張婷婷，喬友鋒，等.三維旅行時場B樣條插值射線追蹤方法［J］.石油地球物理勘探，2013，48(4):559-566.

［6］ Zhang D，Zhang T T，Zhang X L.A new 3-D ray tracing method based on LTI using successive partitioning of cell interfacesand traveltime gradients［J］.Journalof Applied Geophysics，2013，92:20-29.

［7］葉佩，李慶春.旅行時線性插值射線追蹤提高計算精度和效率的改進方法［J］.吉林大學學報，2013，43(1):291-298.

［8］ Moser T J.Shortest path calculation of seismic rays［J］.Geophysics，1991，56(1):59-67.

［9］桑運云，李振春，張凱.拋物旅行時插值最短路徑射線追蹤［J］.石油地球物理勘探，2013，48(3):403-409.

［10］王慶林，白超英.相空間波前構造法地震多值射線追蹤的改進［J］.石油地球物理勘探，2013，48(1):15-21.

［11］徐 濤，徐果明，高爾根，等.三維復雜介質(zhì)的塊狀建模和試射射線追蹤［J］.地球物理學報，2004，47(6):1118-1126.

［12］張釙，劉 洪，李幼銘.射線追蹤方法的發(fā)展現(xiàn)狀［J］.地球物理學進展，2000，15(1):36-45.