本發明屬于地震勘探技術領域，尤其涉及一種基于方向矢量的地震射線追蹤方法。

背景技術：

射線追蹤作為一種重要的正演方法在(微)地震定位、地震層析成像中有著廣泛的應用。它的運算效率及計算精度不僅影響著地震定位的速度和精度還影響反演的效率和結果的可靠性，因而提高射線追蹤的效率和精度將具有重要的意義。試射法(wesson,etal.,1971；chander,etal.,1977；langan,etal.,1985；sunetal.,1993)和基于邊值的彎曲法(jacob,etal.,1970；sorrels,etal.,1971；sleep,etal.,1973)奠定了射線追蹤的基本思路。而后發展的射線追蹤方法在計算效率和計算精度上有較大的提高，如發明人早期所提出的逐段迭代射線追蹤方法(高爾根等，1996；高爾根等，1998；高爾根等，2002；高爾根等，2002；gaoe,etal.,2002；gaoe,etal.,2008)。射線路徑與模型的交點稱為射線的中間點。逐段迭代射線追蹤方法的思想是在滿足射線參數相同的基礎上不斷調整中間點位置以滿足范數判斷條件從而得到射線的傳播路徑，該方法在水平層狀介質、二維速度隨機分布模型以及任意復雜界面下可以快速且準確地追蹤到射線路徑。但隨著應用的深入，對具有特定功能的射線追蹤算法的需求越發強烈，逐段迭代射線追蹤算法難以計算已知射線初始位置及方向矢量快速求取射線路徑及走時的要求。中國專利(cn103698813a)、中國專利(cn105068133a)及中國專利(cn1712991a)提出的射線追蹤算法皆未能給出透射射線、反射射線的顯式計算公式，計算中需依賴震源點、接收點位置作為初始邊界條件，并根據時間最小原理進行迭代。多次迭代的求解方式使得在計算單一射線的路徑及走時的效率低。不能計算以給定方向矢量出射，且在射線旅行時間t約束下的射線路徑及射線所代表的波陣面的位置，因而難以滿足實際需求。

綜上所述，現有技術存在的問題是：

(1)未能給出透射波、反射波的計算顯式公式，依賴于迭代算法進行求解，射線追蹤效率低。

(2)射線追蹤算法需依賴震源點、接收點的初始位置作為邊界條件。實際情況中，震源點位置難以確定，邊界條件退化為單一接收點空間位置及射線矢量，依賴震源點及接收點位置作為邊界條件的迭代式射線追蹤算法并不能快速地追蹤射線的路徑及走時。

技術實現要素：

針對現有技術存在的問題，本發明提供了一種基于方向矢量的地震射線追蹤方法。

本發明是這樣實現的，一種基于方向矢量的地震射線追蹤方法，所述基于方向矢量的地震射線追蹤方法包括以下步驟：

步驟一，根據射線起始點p1(x1,y1,z1)以方向矢量r1(rx1,ry1,rz1)及三維非均勻起伏地質界面的深度范圍，判斷與射線相交且射線走時最短的界面；

步驟二，根據射線上的點p1(x1,y1,z1)及射線的方向矢量r1(rx1,ry1,rz1)，構建射線的直線方程l(x,y,z)；根據地質界面的深度范圍，構造界面σ的最小逼近平面方程zmin及最大平面方程zmax；

步驟三，求解射線的直線方程l(x,y,z)與最小逼近平面zmin的交點dmin，與最大逼近平面zmax的交點dmax；判斷交點是否超出模型邊界，若超出模型區域則沿著射線移至模型區域內；將新得到點設置為初始點dmin及dmax；

步驟四，根據初始點dmin及dmax的位置，構造射線與三維非均勻起伏地質界面的求交迭代方程，計算射線與界面的交點d；

步驟五，基于fermat原理推導的三維空間中的透射波方程及反射波方程，并利用該方程求解射線在非均勻起伏界面處的反射rf波矢量及透射波射線矢量rt；

步驟六，根據所選定的目標層l或設定的射線旅行總時長t，判斷射線追蹤的終止條件；

步驟七，按照應用需求輸出不同的數據組合。

進一步，所述步驟一根據三維非均勻起伏地質模型的深度范圍判斷與射線相交且走時最短的界面，是通過對比射線的起始點p1(x1,y1，z1)與該點在各個地質界面上的投影點p1p(x1p,y1p,z1p)的深度來實施。

進一步，所述步驟二根據地質界面的深度范圍構造界面σ的最小逼近平面方程zmin及最大逼近平面方程zmax；其所構建的逼近平面為水平層且構成三維非均勻起伏地質界面的上下界。

進一步，所述步驟三判斷交點是否超出模型邊界，并將超出界面的交點移至模型區域內部；判斷射線與逼近平面的交點是否超出模型范圍是通過對比交點坐標與模型范圍來實施的；若交點超出模型的范圍，則將交點沿射線移至模型范圍內部；若交點在模型區域內，則將交點設置為求交算法的初始點。

進一步，所述步驟四構造射線與三維非均勻起伏地質界面的求交迭代方程；方程是以射線的方向矢量為約束，以射線與三維非均勻起伏界面的距離為判斷條件，根據射線與三維非均勻起伏地質界面距離差構建目標函數：

f(x,y)＝l(x,y)-f(x,y)；

其中z＝l(x,y,z)為射線的直線方程，z＝f(x,y)為界面的一般性方程。在交點p0(x0,y0，z0)處滿足f(p0)＝0。

進一步，以步驟三的結果(點p1(x1,y1，z1)及點p2(x2,y2，z2))為起始點，且這兩點分布在界面兩側的點，可構造三維空間中的求根的迭代算法：

(1)計算f(p1)和f(p2)的值，并比較其絕對值；

(2)若|f(p1)|＞|f(p2)|則利用：

計算出p的坐標作為新的p1點；

(3)若|f(p1)|＜|f(p2)|則利用：

計算出p的坐標作為新的p2點；

(4)若p1、p2滿足：

||p1p2||2＜ε；

則可得交點p0的坐標，迭代結束；若不滿足重復迭代步驟(1)-步驟(4)。

進一步，所述步驟五中計算射線在三維非均勻起伏地質界面處的透射波方程及反射波方程為：

(1)透射波方程：

其中，r1x，r1y,r1z:表示初始射線矢量r的分量；

v1,v2:表示界面兩側的速度；

f’x,f’y,f’z,:表示交點處界面的法向矢量f’分量；

r2x，r2y,r2z:表示透射波射線矢量r2的分量；

(2)反射波方程：

r2x＝r1x+r1zf′x-r2zf′x；

r2y＝r1y+r1zf′y-r2zf′y；

其中，r1x,r1y,r1z:表示初始射線矢量r1的分量；

v1,v2:表示界面兩側的速度；

f’x,f’y,f’z,:表示交點處界面的法向矢量f’分量；

r2x,r2y,r2z:表示反射波射線矢量r2的分量。

進一步，所述步驟六中終止條件包括：

(1)根據所選定的目標層l；

(2)根據給定的射線旅行時長度t；

(3)根據設定的射線方向矢量v及射線旅行時長度t；

若以選定的目標地層l為終止條件，則可追蹤獲得射線由震源出射，到達目標層l的射線的路徑及走時；若以時間長度t為約束條件，則可獲取由震源出射，走時為t的波陣面的位置；若以射線方向矢量v及射線旅行時長度t為約束條件，則可獲得由震源出射，走時為t，法向為v的波陣面的位置；

所述步驟七中輸出數據包括：射線的路徑、射線的走時、走時為t，法向為v的波陣面的位置。

本發明的另一目的在于提供一種應用所述基于方向矢量的地震射線追蹤方法的地震定儀。

本發明的另一目的在于提供一種應用所述基于方向矢量的地震射線追蹤方法的地震層析成像儀。

本發明的優點及積極效果為：基于空間平面約束的射線矢量與三維任意起伏界面的求交迭代算法，可計算三維空間中任意起伏界面上的射線追蹤問題；基于fermat原理推導的三維任意起伏速度模型上的透射波及反射波的顯式計算公式，可避免利用迭代求解射線中間點；基于不同需求而設計的射線追蹤算法終止條件，可根據旅行時計算波前位置；基于不同需求的算法輸入及輸出設計；能根據射線的矢量信息計算射線的路徑及走時，也能基于方向約束及走時約束獲取波陣面的位置；射線追蹤算法在地震偏移、(微)地震定位及地震層析成像等領域中有廣泛的應用前景。本發明為了計算在給定地震方向矢量的前提下，求取地震射線的路徑及波前的位置。本發明亦能應用于常規地震射線的追蹤問題的求解，因而具有普適性。

附圖說明

圖1是本發明實施例提供的基于方向矢量的地震射線追蹤方法流程圖。

圖2是本發明實施例提供的基于方向矢量的射線追蹤方法的流程示意圖。

圖3是本發明實施例提供的空間平面約束的射線矢量與三維任意起伏界面的求交迭代算法示意圖。

圖4是本發明實施例提供的透射波射線路徑示意圖。

圖5是本發明實施例提供的反射波射線路徑示意圖。

圖6是本發明實施例提供的地質模型及正演射線示意圖。

具體實施方式

為了使本發明的目的、技術方案及優點更加清楚明白，以下結合實施例，對本發明進行進一步詳細說明。應當理解，此處所描述的具體實施例僅僅用以解釋本發明，并不用于限定本發明。

下面結合附圖對本發明的應用原理作詳細的描述。

如圖1所示，本發明實施例提供的基于方向矢量的地震射線追蹤方法包括以下步驟：

s101：根據射線起始點p1(x1,y1,z1)以方向矢量r1(rx1,ry1,rz1)及三維非均勻起伏地質界面的深度范圍，判斷與射線相交且射線走時最短的界面；

s102：根據射線上的點p1(x1,y1,z1)及射線的方向矢量r1(rx1,ry1,rz1)，構建射線的直線方程l(x,y,z)；根據地質界面的深度范圍，構造界面σ的最小逼近平面方程zmin及最大平面方程zmax；

s103：求解射線的直線方程l(x,y,z)與最小逼近平面zmin的交點dmin，與最大逼近平面zmax的交點dmax；判斷交點是否超出模型邊界，若超出模型區域則沿著射線移至模型區域內；將新得到點設置為初始點dmin及dmax；

s104：根據初始點dmin及dmax的位置，構造射線與三維非均勻起伏地質界面的求交迭代方程，計算射線與界面的交點d；

s105：基于fermat原理推導的三維空間中的透射波方程及反射波方程，并利用該方程求解射線在非均勻起伏界面處的反射rf波矢量及透射波射線矢量rt；

s106：根據所選定的目標層l或設定的射線旅行總時長t，判斷射線追蹤的終止條件；

s107：按照應用需求輸出不同的數據組合。

下面結合附圖對本發明的應用原理作進一步的描述。

如圖2所示，本發明實施例提出的基于方向矢量的地震射線追蹤方法具體實現步驟包括：

步驟一，根據射線起始點p1(x1,y1,z1)以方向矢量r1(rx1,ry1,rz1)及三維非均勻起伏地質界面的深度范圍，判斷與射線相交且射線走時最短的界面。該步驟作為求解射線矢量與三維非均勻起伏地質界面的初始條件，將確定與射線最先相交的地質界面σ。

步驟二，根據射線上的點p1(x1,y1,z1)及射線的方向矢量r1(rx1,ry1,rz1)，構建射線的直線方程l(x,y,z)。根據地質界面的深度范圍，構造界面σ的最小逼近平面方程zmin及最大平面方程zmax。

步驟三，求解射線的直線方程l(x,y,z)與最小逼近平面zmin的交點dmin，與最大逼近平面zmax的交點dmax。判斷交點是否超出模型邊界，若超出模型區域則沿著射線移至模型區域內。將新得到點設置為初始點dmin及dmax。

步驟四，根據初始點dmin及dmax的位置，構造射線與三維非均勻起伏地質界面的求交迭代方程，計算射線與界面的交點d。

步驟五，基于fermat原理推導的三維空間中的透射波方程及反射波方程。并利用該方程求解射線在非均勻起伏界面處的反射rf波矢量及透射波射線矢量rt。

步驟六，根據所選定的目標層l或設定的射線旅行總時長t，判斷射線追蹤的終止條件。

步驟七，按照應用需求輸出不同的數據組合。

圖3為基于空間平面約束的射線矢量與三維任意起伏界面的求交迭代算法示意圖，求交迭代方程以射線的方向矢量為約束，以射線與三維非均勻起伏界面的距離為判斷條件，根據射線與三維非均勻起伏地質界面距離差構建目標函數：

f(x,y)＝l(x,y)-f(x,y)；

其中z＝l(x,y,z)為射線的直線方程，z＝f(x,y)為界面的一般性方程。在交點p0(x0,y0,z0)處滿足f(p0)＝0。以步驟103的結果(點p1(x1,y1，z1)及點p2(x2,y2，z2))為起始點，且這兩點分布在界面兩側的點，可構造三維空間中的求根的迭代算法：

(1)計算f(p1)和f(p2)的值，并比較其絕對值。

(2)若|f(p1)|＞|f(p2)|則利用：

計算出p的坐標作為新的p1點。

(3)若|f(p1)|＜|f(p2)|則利用：

計算出p的坐標作為新的p2點。

(4)若p1、p2滿足：

||p1p2||2＜ε；

則可得交點p0的坐標，迭代結束；若不滿足重復迭代步驟(1)-(4)。

圖4為透射波射線路徑示意圖，利用fermat原理計算射線在界面的處的透射波方程方程，其應用公式為：

其中，r1x,r1y,r1z:表示初始射線矢量r1的分量；

v1,v2:表示界面兩側的速度；

f’x,f’y,f’z,:表示交點處界面的法向矢量f’分量；

r2x，r2y,r2z:表示透射波射線矢量r2的分量；

圖5為透射波射線路徑示意圖，利用fermat原理計算射線在界面的處的反射波方程方程，其應用公式為：

r2x＝r1x+r1zf′x-r2zf′x；

r2y＝r1y+r1zf′y-r2zf′y；

其中，r1x，r1y,r1z:表示初始射線矢量r的分量；

v1,v2:表示界面兩側的速度；

f’x,f’y,f’z,:表示交點處界面的法向矢量f’分量；

r2x,r2y,r2z:表示反射波射線矢量r2的分量；

圖6為三維空間中多層起伏地質模型中射線追蹤的正演算例，可見射線路徑與真實的射線路徑相近，這也表明了該射線追蹤算法的可靠性。

以上所述僅為本發明的較佳實施例而已，并不用以限制本發明，凡在本發明的精神和原則之內所作的任何修改、等同替換和改進等，均應包含在本發明的保護范圍之內。