高效高精度地震波初至行走計(jì)算

高效高精度地震波初至行走計(jì)算

0走時(shí)計(jì)算方法計(jì)算各向異性偏移時(shí)間是基于全波形反演、色度成像和其他各向異性地震模擬和成像技術(shù)的基礎(chǔ)。走時(shí)計(jì)算方法的選擇需要綜合考慮算法的精度、效率及穩(wěn)定性等因素。層析成像是一種層狀介質(zhì)速度和各向異性參數(shù)建模方法,其結(jié)果往往用于為全波形反演提供初始模型,但與波形反演中的旅行時(shí)信息不同,其核心之一就是模型初至波走時(shí)計(jì)算。層析成像是一種迭代反演算法,需要通過(guò)正演結(jié)果與實(shí)際數(shù)據(jù)的殘差不斷調(diào)整初始模型。走時(shí)計(jì)算的精度、效率與穩(wěn)定性和層析成像計(jì)算的精度、效率與穩(wěn)定性息息相關(guān)。因此,高效高精度且穩(wěn)定的走時(shí)計(jì)算方法是各向異性介質(zhì)偏移和層析成像等建模技術(shù)實(shí)現(xiàn)的關(guān)鍵。傳統(tǒng)的走時(shí)計(jì)算方法是射線(xiàn)追蹤方法,包括試射法(shootingtrace)和彎曲法(bendingtrace)射線(xiàn)追蹤技術(shù)結(jié)合圖形理論和復(fù)雜網(wǎng)格技術(shù),在局部進(jìn)行射線(xiàn)追蹤。主要衍生出的方法有插值法和最短路徑射線(xiàn)追蹤算法等。插值法走時(shí)計(jì)算的另一種方法是波前面追蹤,其包括兩種方法。一種方法是基于運(yùn)動(dòng)學(xué)射線(xiàn)追蹤系統(tǒng)的波前構(gòu)建法;另一種方法是通過(guò)求解波動(dòng)方程WKB(Wentzel,Kramers和Brillouin)展開(kāi)的零階近似后來(lái)衍生出的方法根據(jù)波前擴(kuò)展方法的不同可分為迎風(fēng)差分法FMM是由Sethian國(guó)內(nèi)對(duì)于各向異性走時(shí)計(jì)算的方法研究起步較晚,其中以何樵登在各向異性介質(zhì)中,走時(shí)計(jì)算不只依賴(lài)速度場(chǎng),隨著各向異性參數(shù)的增多,構(gòu)建迎風(fēng)差分格式求取走時(shí)也越來(lái)越困難。本文主要利用擾動(dòng)理論及泰勒公式將VTI介質(zhì)程函方程展開(kāi),完成了對(duì)VTI介質(zhì)程函方程的降維處理。對(duì)降維后的VTI介質(zhì)程函方程構(gòu)建迎風(fēng)差分格式并分別進(jìn)行求解,完成單點(diǎn)走時(shí)計(jì)算。結(jié)合窄帶推進(jìn)理論,構(gòu)成了各向異性的波前面追蹤走時(shí)計(jì)算方法。首先對(duì)均勻弱VTI各向異性介質(zhì)進(jìn)行計(jì)算并與解析值相減評(píng)估其誤差;然后為驗(yàn)證其穩(wěn)定性和有效性,又計(jì)算了VTI各向異性層狀介質(zhì)模型的走時(shí);最后以鹽丘模型作為測(cè)試模型,進(jìn)行了走時(shí)正演試算。1算法的原理和實(shí)現(xiàn)1.1vti介質(zhì)程函方程離散方程在VTI介質(zhì)中,聲波近似的程函方程為以下格式式中:τ(x,y,z)是震源到坐標(biāo)點(diǎn)(x,y,z)的走時(shí);v此四次方程的非線(xiàn)性程度強(qiáng)于各向同性和橢圓各向異性的程函方程,因此VTI介質(zhì)程函方程需要更為復(fù)雜的有限差分近似。為解此方程,我們可以利用擾動(dòng)理論式中,τ方程(4)是橢圓各向異性介質(zhì)的程函方程,獨(dú)立參數(shù)η相同階數(shù)項(xiàng)的系數(shù)相等。首先計(jì)算一次η項(xiàng)的系數(shù),將一次η項(xiàng)系數(shù)相等的公式通過(guò)式(4)簡(jiǎn)化后得到同理,可得到η平方項(xiàng)的系數(shù)公式:1.2基于橢圓異性介質(zhì)模型的背景模型為求取VTI介質(zhì)走時(shí)的解析解,首先從均勻介質(zhì)背景速度模型出發(fā)。為求取式(4)的解析解,設(shè)η=0,此時(shí)背景模型是橢圓各向異性介質(zhì)。在這種情況下,震源位于x=0、z=0處的走時(shí)解析解如下:該結(jié)果滿(mǎn)足方程(4)。通過(guò)式(7)可求出同理,計(jì)算將式(7)、(8)、(9)代入式(3)中,即可求取該點(diǎn)走時(shí)的解析值。1.3非均勻介質(zhì)走時(shí)計(jì)算在均勻介質(zhì)中,上述公式的推導(dǎo)是基于VTI介質(zhì)聲波近似程函方程(1)進(jìn)行的,因此得到的近似結(jié)果描述的是介質(zhì)的群速度。而在非均勻介質(zhì)正方形網(wǎng)格條件下,由于無(wú)法確定指定網(wǎng)格中的震源位置,因此上述公式不能簡(jiǎn)單地應(yīng)用于非均勻介質(zhì)走時(shí)計(jì)算中。對(duì)于此種情況,我們引入FMM。FMM的基本原理是將程函方程的時(shí)間梯度項(xiàng)近似為滿(mǎn)足波前傳播“熵守恒”理論的數(shù)值離散形式,通過(guò)數(shù)值計(jì)算來(lái)求取黏滯解。這種離散形式為式中:D該方法對(duì)于波前傳播過(guò)程中尖角拐角和拓?fù)浣Y(jié)構(gòu)變化等問(wèn)題的有效性在文獻(xiàn)然后,將式(12)代入式(4)中,通過(guò)計(jì)算可以求出點(diǎn)(i,j)處的τ1.4窄帶時(shí)走時(shí)計(jì)算FMM的實(shí)現(xiàn)策略主要為“熵守恒”理論,該理論可以這樣理解:傳播的初至波波前可以看做在紙上燃燒的火焰。如果紙上的某個(gè)位置被火焰燃燒了,那么這一位置將一直保持著被燒過(guò)的狀態(tài),后續(xù)不會(huì)再被燃燒。類(lèi)比于初至波在空間中的傳播,即初至波波前不會(huì)兩次經(jīng)過(guò)同一個(gè)點(diǎn),每個(gè)點(diǎn)的走時(shí)是關(guān)于位置的單值函數(shù)。同時(shí)由于波是單向傳播的,總是由走時(shí)較小點(diǎn)傳播至走時(shí)較大點(diǎn),因此,1.3節(jié)中所述單點(diǎn)走時(shí)計(jì)算的循環(huán)起點(diǎn)應(yīng)位于走時(shí)最小值點(diǎn)。文中采用了窄帶技術(shù)作為波前模擬技術(shù),通過(guò)排序選擇窄帶內(nèi)走時(shí)最小的點(diǎn)作為延拓點(diǎn),再通過(guò)對(duì)網(wǎng)格節(jié)點(diǎn)屬性的設(shè)定來(lái)“儲(chǔ)存”已經(jīng)完成計(jì)算的網(wǎng)格節(jié)點(diǎn),最后通過(guò)窄帶內(nèi)網(wǎng)格節(jié)點(diǎn)的走時(shí)更新來(lái)模擬波前的演化。整體步驟可分為初始化和循環(huán)兩部分。12啟動(dòng)式中:l為網(wǎng)格節(jié)點(diǎn)的間隔長(zhǎng)度;vi2確定窄帶以確定深度循環(huán)過(guò)程分為以下4個(gè)步驟。(1)從窄帶內(nèi)選取最小走時(shí)點(diǎn)τ(i(2)將點(diǎn)i((3)對(duì)點(diǎn)i((4)判斷窄帶是否為空,如果窄帶為空則終止循環(huán),否則跳回(1)繼續(xù)計(jì)算。2算法的精度和穩(wěn)定性分析2.1走時(shí)計(jì)算與解析值對(duì)比由于復(fù)雜界面非均勻模型的走時(shí)解析值較難求取,為了便于分析本文方法的正確性及精度,我們首先用一個(gè)水平界面均勻弱VTI各向異性模型計(jì)算走時(shí),并用走時(shí)的解析值與之對(duì)比(由于后文復(fù)雜模型通過(guò)速度及各向異性參數(shù)δ和ε定義,因此本文中所有模型均用v水平界面均勻介質(zhì)模型的大小為4000m×1200m,介質(zhì)速度為4000m/s,各向異性參數(shù)δ=0.1,ε=0.1,網(wǎng)格間距為5m,震源位置為(2000m,600m)。使用本文方法的走時(shí)計(jì)算結(jié)果見(jiàn)圖1a。走時(shí)解析值通過(guò)1.2節(jié)中敘述的方法求取,將走時(shí)計(jì)算值(圖1a)與解析值相減后得到誤差,誤差占解析值的比例(相對(duì)誤差)見(jiàn)圖1b。由圖1b可知:在震源附近走時(shí)計(jì)算結(jié)果的相對(duì)誤差較為明顯,主要分布在1.5%～4.0%之間,其原因是差分網(wǎng)格長(zhǎng)度太大,對(duì)震源處的差分網(wǎng)格進(jìn)行進(jìn)一步細(xì)分,該現(xiàn)象會(huì)有所改善;在遠(yuǎn)偏處誤差趨于穩(wěn)定,相對(duì)誤差穩(wěn)定在0.5%以下。2.2算法的穩(wěn)定性為了驗(yàn)證本算法對(duì)復(fù)雜介質(zhì)模型的穩(wěn)定性和有效性,分別選取雙層、多層和鹽丘介質(zhì)模型進(jìn)行試算。2.2.1走時(shí)模型試算對(duì)雙層VTI各向異性介質(zhì)模型進(jìn)行數(shù)據(jù)試算。該模型的大小為4000m×1200m,震源點(diǎn)位于坐標(biāo)(0,0)處。模型參數(shù)見(jiàn)表1,模型結(jié)構(gòu)見(jiàn)圖2a,使用本文方法試算的走時(shí)結(jié)果見(jiàn)圖2b。地震觀(guān)測(cè)點(diǎn)最先接收到的波稱(chēng)為初至波。在層狀介質(zhì)中,隨著偏移距及入射角的增大,下層介質(zhì)的折射波漸漸取代直達(dá)波被率先接收,在某節(jié)點(diǎn)產(chǎn)生“毛刺現(xiàn)象”。通過(guò)觀(guān)察圖2b可以發(fā)現(xiàn),第一層近偏處與均勻介質(zhì)傳播規(guī)律一致,隨著偏移距增加,在坐標(biāo)(500,600)處產(chǎn)生了“毛刺現(xiàn)象”,折射波取代了直達(dá)波率先被接收。2.2.2多層模型測(cè)試結(jié)果對(duì)多層VTI各向異性介質(zhì)模型進(jìn)行數(shù)據(jù)試算。該模型的大小為4000m×1200m,震源點(diǎn)位于坐標(biāo)(0,0)處。模型參數(shù)見(jiàn)表2,模型結(jié)構(gòu)見(jiàn)圖3a,使用本文方法試算的走時(shí)結(jié)果見(jiàn)如圖3b。由于多層模型的速度變化與雙層模型相比較慢,因此其等時(shí)線(xiàn)變化較為平緩,近偏處對(duì)速度突變層的表現(xiàn)不明顯。但在偏移距1000～2500m處,我們可以較清晰地分辨出3個(gè)層的折射波對(duì)于初至走時(shí)的影響,偏移距2500m之后,第二、三層的折射波漸漸被第四層的折射波所取代,并體現(xiàn)于等時(shí)線(xiàn)上。2.2.3波場(chǎng)快照的正演模擬對(duì)鹽丘介質(zhì)模型進(jìn)行數(shù)據(jù)試算。該模型大小為21940m×9144m,模型速度及各向異性參數(shù)結(jié)構(gòu)如圖4所示。該模型x軸2000～10000m處為一較明顯的高速體,15000～19000m處為一斷層,兩個(gè)比較突出的結(jié)構(gòu)體之間是由起伏層狀介質(zhì)連接(圖4a)。在層狀地質(zhì)體中,x軸10000～14000m、深度約4000m處有一強(qiáng)各向異性體(圖4b、c)。我們?cè)谀Ｐ蛒軸方向0、3650、12200、19000m處分別設(shè)置炮點(diǎn),觀(guān)察其在特殊地形下的試算效果,并與有限差分正演的波場(chǎng)快照進(jìn)行對(duì)比,結(jié)果如圖5所示。從圖5a可以看出,高速體對(duì)走時(shí)場(chǎng)的影響十分明顯。通過(guò)與波場(chǎng)快照(圖5b)對(duì)比可以發(fā)現(xiàn),排除波場(chǎng)快照中界面之間反射波的影響,3000ms等值線(xiàn)與波場(chǎng)快照中的波前面完全吻合。由于圖5c中炮點(diǎn)位置位于高速體深度最小的地面投影處,因此選擇2000ms的波場(chǎng)快照(圖5d)進(jìn)行對(duì)比?？梢钥闯?波場(chǎng)快照0～3000m深度處與等時(shí)線(xiàn)吻合;5000m以下受到高速體的影響,波場(chǎng)快照中波前發(fā)生分裂,但是最高能量的波前面與等時(shí)線(xiàn)依然吻合。排除反射波的影響,從圖5e、f可以清晰地看到,不論是在深度0～3000m的層狀介質(zhì)處,還是在深度3000～5000m的高速體和斷層處,波場(chǎng)快照都與波前面清晰地吻合。由于圖5g中震源點(diǎn)的位置是在斷層正上方,為觀(guān)察到高速體的影響,因此將波場(chǎng)快照的時(shí)間設(shè)置為5500ms(圖5h)?？梢郧逦乜吹缴疃?～3000m的波前面與等時(shí)線(xiàn)吻合;而對(duì)于高速體,由于偏移距過(guò)大的原因,有限差分正演受到反射波、繞射波等因素的影響,波場(chǎng)快照中最高能量波前面較難分辨。綜上所述,經(jīng)過(guò)多個(gè)模型的測(cè)試,本文提出的方法在不同地下構(gòu)造及不同偏移距下的穩(wěn)定性和準(zhǔn)確性得到了較好的驗(yàn)證,其計(jì)算結(jié)果合乎地震波傳播規(guī)律,能夠準(zhǔn)確地反映出地下介質(zhì)的結(jié)構(gòu)。3vti介質(zhì)聲波走時(shí)計(jì)算方法的適用性1)本文將基于波前面追蹤的走時(shí)計(jì)算方法引入各向異性介質(zhì)走時(shí)計(jì)算當(dāng)中,實(shí)現(xiàn)了網(wǎng)格化非均勻VTI介質(zhì)走時(shí)計(jì)算。首先對(duì)程函方程中的走時(shí)參數(shù)τ進(jìn)行展開(kāi),成功將四次方程分解為一個(gè)二次方程和兩個(gè)一次方程;然后將各向同性快速匹配走時(shí)計(jì)算方法引入VTI介質(zhì)中,對(duì)展開(kāi)的VTI介質(zhì)程函方程組進(jìn)行離散化并求解,得到介質(zhì)的走時(shí)。2)大量模型試算檢驗(yàn)證明,本文提出的VTI介質(zhì)聲波走時(shí)計(jì)算方法是有效的,該方法的模型試算結(jié)果能夠較好地反映出地下復(fù)雜介質(zhì)中地震波的傳播規(guī)律。雖然本文中用于試算的模型不能涵蓋地下介質(zhì)的所有情況,但本文中的計(jì)算方