本發(fā)明屬于地震監(jiān)測(cè)技術(shù)領(lǐng)域，涉及一種高精度地震波走時(shí)射線追蹤方法。

背景技術(shù)：

傳統(tǒng)的射線追蹤方法包括試射法和彎線法，試射法是從震源出發(fā)，以一定的角度射出一條射線，射線根據(jù)snell定律進(jìn)行傳播，然后根據(jù)射線在地面上的落點(diǎn)和實(shí)際接收點(diǎn)的情況不斷調(diào)整初射角度，最后由最靠近接收點(diǎn)的兩條射線走時(shí)差值求出接收點(diǎn)的走時(shí)和射線路徑。彎線法是先假設(shè)震源和接收點(diǎn)之間存在一條初始的射線路徑，然后根據(jù)費(fèi)馬原理對(duì)射線路徑進(jìn)行改動(dòng)，從而求出接收點(diǎn)的最小走時(shí)和射線路徑。這兩種方法時(shí)常面臨陷于局部解的困擾，有時(shí)候不能求得全局最小解。

LTI算法是日本科學(xué)家Asakawa于1993年提出的一種射線追蹤方法，該方法分為向前-向后計(jì)算2個(gè)基本步驟，并以線性插值為基本假設(shè)。向前計(jì)算的目的是為了得到各個(gè)網(wǎng)格節(jié)點(diǎn)上的最小走時(shí)，向后計(jì)算是在向前計(jì)算的基礎(chǔ)上，由檢波點(diǎn)根據(jù)走時(shí)最小的原理逐步追蹤到炮點(diǎn)。

射線追蹤的基本理論是，在高頻近似的條件下，地震波的主要能量沿著射線路徑傳播。而射線追蹤的目的就是在給定速度模型(m行n列的網(wǎng)格，每個(gè)網(wǎng)格內(nèi)的速度為常數(shù)，不同網(wǎng)格間的速度可能不一樣)、震源點(diǎn)位置和接收點(diǎn)位置時(shí)，求出一條射線的傳播路徑，使得從震源點(diǎn)到接收點(diǎn)的傳播時(shí)間(走時(shí))最小。

旅行時(shí)線性插值(Linear Traveltime Interpolation，LTI)法射線追蹤具有較高的計(jì)算精度和效率，常用于二維速度模型的射線追蹤。然而，該方法在考慮射線的來向時(shí)只考慮了一個(gè)半平面，即炮點(diǎn)相對(duì)于網(wǎng)格節(jié)點(diǎn)所在的半平面，沒有考慮出現(xiàn)回折波的可能；對(duì)于地質(zhì)情況復(fù)雜時(shí)，不能適應(yīng)復(fù)雜速度模型的射線追蹤。

下面通過建模來說明，給出一種速度模型，如圖1所示，模型為1000m＊1000m的正方形，網(wǎng)格大小為10m＊10m，其中0～400m和600m～1000m深度上為1000m/s的低速體，中間400m～600m為1500m/s的高速體，該模型最先由Asakawa建立，用來驗(yàn)證射線追蹤的準(zhǔn)確性，本文稱之為Asakawa模型。現(xiàn)在以(0，200)的位置為震源，(1000，0)、(1000,100)、……、(1000，900)的位置為接收點(diǎn)，得到的射線路徑如圖1中的曲線所示。需要注意的是射線3和射線4兩條射線是經(jīng)過高速體的折射后到達(dá)的?，F(xiàn)在把模型順時(shí)針旋轉(zhuǎn)90°，震源和接收點(diǎn)相對(duì)于模型的位置不變，仍然用傳統(tǒng)LTI算法進(jìn)行射線追蹤，結(jié)果如圖2所示。對(duì)比圖1和圖2可以發(fā)現(xiàn)，射線3和射線4兩條射線的路徑不一樣，圖1中的兩條射線是經(jīng)過高速體的折射到達(dá)接收點(diǎn)的，而圖2中的兩條射線是直接由震源點(diǎn)到達(dá)接收點(diǎn)的。表1給出了圖1和圖2中射線走時(shí)的對(duì)比，可以看出，Asakawa模型旋轉(zhuǎn)后射線3和射線4兩條射線走時(shí)變大，說明傳統(tǒng)LTI算法在面對(duì)復(fù)雜模型時(shí)計(jì)算精度不夠。

表1：圖1和圖2中射線走時(shí)對(duì)比

技術(shù)實(shí)現(xiàn)要素：

為實(shí)現(xiàn)上述目的，本發(fā)明提供一種高精度地震波走時(shí)射線追蹤方法，解決了現(xiàn)有技術(shù)中由于在考慮射線來向的問題上存在缺陷，算法的適應(yīng)能力不強(qiáng)的問題。

本發(fā)明所采用的技術(shù)方案是，高精度地震波走時(shí)射線追蹤方法，按照以下步驟進(jìn)行：

步驟1，計(jì)算得到整個(gè)區(qū)域網(wǎng)格節(jié)點(diǎn)上的走時(shí)；

步驟2，除震源所在網(wǎng)格外，重新計(jì)算震源所在行網(wǎng)格上節(jié)點(diǎn)的走時(shí)，若比原來的小就替代；

步驟3，對(duì)于震源所在行之上的各行，先以網(wǎng)格下邊界為插值線段，重新計(jì)算網(wǎng)格上邊界上的2個(gè)節(jié)點(diǎn)的走時(shí)，若果比原來的小就替代；

步驟4，再以網(wǎng)格左邊界為插值線段，重新計(jì)算網(wǎng)格右邊界上2個(gè)節(jié)點(diǎn)的走時(shí)，若比原來的小就替代；

步驟5，再以網(wǎng)格右邊界為插值線段，重新計(jì)算網(wǎng)格左邊界上2個(gè)節(jié)點(diǎn)的走時(shí)，若比原來的小就替代；

步驟6，重復(fù)步驟3-5，直至到達(dá)空間內(nèi)的第一行。

進(jìn)一步的，所述步驟1中，計(jì)算得到整個(gè)區(qū)域網(wǎng)格節(jié)點(diǎn)上的走時(shí)具體按照以下步驟進(jìn)行：

步驟a，計(jì)算震源所在網(wǎng)格上4個(gè)節(jié)點(diǎn)的走時(shí)；

步驟b，計(jì)算震源所在列網(wǎng)格上各節(jié)點(diǎn)的最小走時(shí)，包括以網(wǎng)格下邊界為插值線段計(jì)算上方2個(gè)節(jié)點(diǎn)的最小走時(shí)，以及以網(wǎng)格上邊界為插值線段計(jì)算下方2個(gè)節(jié)點(diǎn)的最小走時(shí)，以此類推，直到計(jì)算到第一行或最后一行；

步驟c，計(jì)算炮點(diǎn)所在列右側(cè)一列網(wǎng)格節(jié)點(diǎn)走時(shí)，以網(wǎng)格左邊界為插值線段計(jì)算網(wǎng)格右方2個(gè)節(jié)點(diǎn)的最小走時(shí)，對(duì)于重復(fù)計(jì)算的節(jié)點(diǎn)，若計(jì)算所得值比原來的小就替代；

步驟d，計(jì)算炮點(diǎn)所在列右側(cè)一列網(wǎng)格節(jié)點(diǎn)走時(shí)，由下往上，以網(wǎng)格下邊邊界為插值線段，計(jì)算網(wǎng)格上方2個(gè)節(jié)點(diǎn)的最小走時(shí)，如果比原來的值小就替代；

步驟e，計(jì)算炮點(diǎn)所在列右側(cè)一列網(wǎng)格節(jié)點(diǎn)走時(shí)，由上往下，以網(wǎng)格上邊界為插值線段計(jì)算網(wǎng)格下邊界2個(gè)節(jié)點(diǎn)的最小走時(shí)，如果比原來的值小就替代；

步驟f，重復(fù)步驟c-e，一直到最右邊的一列網(wǎng)格；

對(duì)于炮點(diǎn)左側(cè)網(wǎng)格節(jié)點(diǎn)最小走時(shí)的計(jì)算，采用和右側(cè)相對(duì)應(yīng)的方式得到。

本發(fā)明的有益效果是：在傳統(tǒng)LTI算法的基礎(chǔ)上進(jìn)行了改進(jìn)，考慮了網(wǎng)格節(jié)點(diǎn)射線從各個(gè)方向傳來的可能，實(shí)驗(yàn)結(jié)果表明：改進(jìn)后的LTI方法具有較強(qiáng)的適應(yīng)性，網(wǎng)格節(jié)點(diǎn)走時(shí)計(jì)算精度得到提高，相應(yīng)的射線路徑也更加準(zhǔn)確。

附圖說明

為了更清楚地說明本發(fā)明實(shí)施例或現(xiàn)有技術(shù)中的技術(shù)方案，下面將對(duì)實(shí)施例或現(xiàn)有技術(shù)描述中所需要使用的附圖作簡(jiǎn)單地介紹，顯而易見地，下面描述中的附圖僅僅是本發(fā)明的一些實(shí)施例，對(duì)于本領(lǐng)域普通技術(shù)人員來講，在不付出創(chuàng)造性勞動(dòng)的前提下，還可以根據(jù)這些附圖獲得其他的附圖。

圖1是Asakawa模型由傳統(tǒng)LTI算法追蹤出的射線路徑圖。

圖2是Asakawa模型旋轉(zhuǎn)以后用傳統(tǒng)LTI算法追蹤出的射線路徑圖。

圖3是傳統(tǒng)的旅行時(shí)線性插值(LTI)原理示意圖。

圖4是向前計(jì)算確定各節(jié)點(diǎn)最小走時(shí)的實(shí)現(xiàn)過程圖。

圖5是向后計(jì)算確定射線路徑的過程圖。

圖6是本發(fā)明實(shí)施例算法的實(shí)現(xiàn)過程圖。

圖7是Asakawa模型旋轉(zhuǎn)以后用本發(fā)明方法追蹤出的射線路徑。

具體實(shí)施方式

下面將結(jié)合本發(fā)明實(shí)施例中，對(duì)本發(fā)明實(shí)施例中的技術(shù)方案進(jìn)行清楚、完整地描述，顯然，所描述的實(shí)施例僅僅是本發(fā)明一部分實(shí)施例，而不是全部的實(shí)施例?；诒景l(fā)明中的實(shí)施例，本領(lǐng)域普通技術(shù)人員在沒有做出創(chuàng)造性勞動(dòng)前提下所獲得的所有其他實(shí)施例，都屬于本發(fā)明保護(hù)的范圍。

先詳細(xì)介紹傳統(tǒng)的LTI算法的基本原理：

如圖3所示，求某點(diǎn)P最小走時(shí)的問題可以作如下描述：假設(shè)A和B兩點(diǎn)的坐標(biāo)和走時(shí)已知，求經(jīng)過線段AB上的某點(diǎn)C到達(dá)點(diǎn)P的最小走時(shí)。設(shè)A點(diǎn)為坐標(biāo)原點(diǎn)，AB長(zhǎng)度為L(zhǎng)，AC長(zhǎng)度為r，P點(diǎn)坐標(biāo)為(x，y)，根據(jù)LTI的基本假設(shè)，C點(diǎn)的走時(shí)T_C可以由A和B兩點(diǎn)的走時(shí)T_A和T_B線性插值得到，即：

則P點(diǎn)的走時(shí)為C點(diǎn)的走時(shí)加上CP線段所用的時(shí)間：

式中：ΔT＝T_B-T_A，s為空間內(nèi)的慢度。對(duì)T_P求關(guān)于r的導(dǎo)數(shù)，并令其等于0

可得：

當(dāng)時(shí)，

當(dāng)時(shí)，

對(duì)一階導(dǎo)數(shù)繼續(xù)求導(dǎo)，可以得到T_P關(guān)于r的二階導(dǎo)數(shù)，可以看出二階導(dǎo)數(shù)在這兩個(gè)解處大于0，可見這兩個(gè)解都是函數(shù)T_P的極小值，而T_P又是關(guān)于r的初等函數(shù)，在區(qū)間0≤r≤L內(nèi)連續(xù)，對(duì)比這2個(gè)解以及AB兩端點(diǎn)的值，就可以得出T_P在0≤r≤L區(qū)間內(nèi)的最小值。

下面介紹LTI算法實(shí)現(xiàn)步驟

LTI分為向前計(jì)算和向后計(jì)算兩個(gè)基本步驟：向前計(jì)算是從炮點(diǎn)開始，由公式(4)求出空間內(nèi)所有網(wǎng)格節(jié)點(diǎn)的最小走時(shí)；向后計(jì)算是從接收點(diǎn)開始，由公式(3)確定射線路徑上與各個(gè)網(wǎng)格邊界的交點(diǎn)。

向前計(jì)算的步驟為：

步驟1，計(jì)算震源所在網(wǎng)格上4個(gè)節(jié)點(diǎn)的走時(shí)，如圖4(a)所示；

步驟2，計(jì)算震源所在列網(wǎng)格上各節(jié)點(diǎn)的最小走時(shí)，包括以網(wǎng)格下邊界為插值線段計(jì)算上方2個(gè)節(jié)點(diǎn)的最小走時(shí)，以及以網(wǎng)格上邊界為插值線段計(jì)算下方2個(gè)節(jié)點(diǎn)的最小走時(shí)，以此類推，直到計(jì)算到第一行或最后一行，如圖4(b)所示；

步驟3，計(jì)算炮點(diǎn)所在列右側(cè)一列網(wǎng)格節(jié)點(diǎn)走時(shí)，以網(wǎng)格左邊界為插值線段計(jì)算網(wǎng)格右方2個(gè)節(jié)點(diǎn)的最小走時(shí)，對(duì)于重復(fù)計(jì)算的節(jié)點(diǎn)，若計(jì)算所得值比原來的小就替代，如圖4(c)所示；

步驟4，計(jì)算炮點(diǎn)所在列右側(cè)一列網(wǎng)格節(jié)點(diǎn)走時(shí)，由下往上，以網(wǎng)格下邊邊界為插值線段，計(jì)算網(wǎng)格上方2個(gè)節(jié)點(diǎn)的最小走時(shí)，如果比原來的值小就替代，如圖4(d)所示；

步驟5，計(jì)算炮點(diǎn)所在列右側(cè)一列網(wǎng)格節(jié)點(diǎn)走時(shí)，由上往下，以網(wǎng)格上邊界為插值線段計(jì)算網(wǎng)格下邊界2個(gè)節(jié)點(diǎn)的最小走時(shí)，如果比原來的值小就替代，如圖4(e)所示；

步驟6，重復(fù)步驟3-5，一直到最右邊的一列網(wǎng)格，如圖4(f)。

對(duì)于炮點(diǎn)左側(cè)網(wǎng)格節(jié)點(diǎn)最小走時(shí)的計(jì)算，可采用和右側(cè)相似的方式得到。

向后計(jì)算的過程為：

步驟1，根據(jù)公式t＝t_向前+sd計(jì)算經(jīng)接收點(diǎn)所在網(wǎng)格上4個(gè)節(jié)點(diǎn)到達(dá)接收點(diǎn)的走時(shí)，找出其中最小的如圖5(a)所示，其中t_向前為向前計(jì)算得到的網(wǎng)格上節(jié)點(diǎn)的走時(shí)，s為網(wǎng)格的慢度，d為接收點(diǎn)到節(jié)點(diǎn)距離；

步驟2，根據(jù)步驟1得到的網(wǎng)格節(jié)點(diǎn)，用公式(3)和公式(4)計(jì)算到達(dá)接收點(diǎn)的最小走時(shí)以及對(duì)應(yīng)的插值點(diǎn)；

步驟3，將步驟2中求得的最小值對(duì)應(yīng)的插值點(diǎn)作為新的接收點(diǎn)，重復(fù)步驟1和步驟2，得到下一個(gè)最小走時(shí)對(duì)應(yīng)的插值點(diǎn)，如圖5(b)所示；

步驟4，重復(fù)步驟3，直到求得的插值點(diǎn)到達(dá)炮點(diǎn)所在網(wǎng)格的邊界，如圖5(c)所示；

步驟4，從接收點(diǎn)開始，依次連接最小走時(shí)所對(duì)應(yīng)的插值點(diǎn)，最后到達(dá)震源點(diǎn)，這樣就得到了由震源點(diǎn)到接收點(diǎn)最小走時(shí)的路徑。

因此，結(jié)合背景技術(shù)所說，傳統(tǒng)LTI算法由于在考慮射線來向的問題上存在缺陷，導(dǎo)致算法的適應(yīng)能力不強(qiáng)，現(xiàn)在針對(duì)這一點(diǎn)進(jìn)行改進(jìn)，提出了高精度LTI算法，實(shí)現(xiàn)過程如圖6所示，具體按照以下步驟進(jìn)行：

步驟1，先按照傳統(tǒng)的方法計(jì)算得到整個(gè)區(qū)域網(wǎng)格節(jié)點(diǎn)上的走時(shí)；

步驟2，除震源所在網(wǎng)格外，重新計(jì)算震源所在行網(wǎng)格上節(jié)點(diǎn)的走時(shí)，若比原來的小就替代，如圖6(a)所示；

步驟3，對(duì)于震源所在行之上的各行，先以網(wǎng)格下邊界為插值線段，重新計(jì)算網(wǎng)格上邊界上的2個(gè)節(jié)點(diǎn)的走時(shí)，若果比原來的小就替代，如圖6(b)所示；

步驟4，再以網(wǎng)格左邊界為插值線段，重新計(jì)算網(wǎng)格右邊界上2個(gè)節(jié)點(diǎn)的走時(shí)，若比原來的小就替代，如圖6(c)所示；

步驟5，再以網(wǎng)格右邊界為插值線段，重新計(jì)算網(wǎng)格左邊界上2個(gè)節(jié)點(diǎn)的走時(shí)，若比原來的小就替代，如圖6(d)所示；

步驟6，重復(fù)步驟3-5，直至到達(dá)空間內(nèi)的第一行。

對(duì)于炮點(diǎn)所在網(wǎng)格下方各行節(jié)點(diǎn)走時(shí)的重新計(jì)算，可采用相似的方式。

當(dāng)采用本發(fā)明方法之后，重新按照?qǐng)D2的模型和震源-接收點(diǎn)關(guān)系進(jìn)行射線追蹤，得到的結(jié)果如圖7所示，由圖7可以看出，射線3和射線4兩條射線也是經(jīng)由高速層折射到達(dá)的接收點(diǎn)，這與圖1中的結(jié)果相同，同時(shí)對(duì)比表2中的走時(shí)信息，改進(jìn)后的走時(shí)比原來的更小，也與圖1中對(duì)應(yīng)接收點(diǎn)的走時(shí)相同，說明改進(jìn)后的LTI算法比傳統(tǒng)LTI算法得到的走時(shí)精度更高，追蹤出來的射線路徑也更準(zhǔn)確。

表2：圖2和圖7中射線走時(shí)對(duì)比

以上所述僅為本發(fā)明的較佳實(shí)施例而已，并非用于限定本發(fā)明的保護(hù)范圍。凡在本發(fā)明的精神和原則之內(nèi)所作的任何修改、等同替換、改進(jìn)等，均包含在本發(fā)明的保護(hù)范圍內(nèi)。