基于牛頓迭代法求解群速度的地震波射線追蹤模擬

張定文,李衛(wèi)東,段金龍,張學(xué)海

(1.河南工業(yè)大學(xué)信息科學(xué)與工程學(xué)院,河南 鄭州 450001;2.應(yīng)急管理部國(guó)家自然災(zāi)害防治研究院,北京 100085)

0 引言

研究表明,各向異性在地下介質(zhì)中是普遍存在的[1-2]。傳統(tǒng)的地下構(gòu)造研究多是在各向同性或完全彈性的介質(zhì)中進(jìn)行,忽略了由各向異性導(dǎo)致的群速度計(jì)算誤差而無(wú)法真實(shí)表達(dá)地下構(gòu)造。隨著計(jì)算機(jī)技術(shù),射線追蹤技術(shù)等發(fā)展,基于復(fù)雜各向異性地質(zhì)模型的地下構(gòu)造研究成為熱點(diǎn)[3-7]。

地震波射線追蹤作為一種快速有效的地震波場(chǎng)數(shù)值模擬方法,能夠直觀展現(xiàn)地震波的幾何傳播路徑,清晰表達(dá)地下內(nèi)部構(gòu)造不均勻性以及速度結(jié)構(gòu)各向異性,因此廣泛應(yīng)用于地震正演模擬,層析成像,偏移成像等領(lǐng)域。目前地震波射線追蹤方法眾多,主流的有兩點(diǎn)射線追蹤法、最短路徑法、有限差分法、走時(shí)插值法和波前構(gòu)建法等方法[8-9]。其中,最短路徑射線追蹤法最早由Nakanishi等[10]在1986年引入到地震波走時(shí)及射線路徑計(jì)算領(lǐng)域,之后國(guó)內(nèi)外學(xué)者相繼對(duì)算法進(jìn)行改進(jìn)和系統(tǒng)化[11-15]。通過(guò)對(duì)地震波進(jìn)行射線追蹤模擬,可以獲取地下速度結(jié)構(gòu)物性信息,獲取地下介質(zhì)中地震波走時(shí)分布,實(shí)現(xiàn)對(duì)地下介質(zhì)各向異性的直觀表達(dá),有效解讀地質(zhì)結(jié)構(gòu)。

在地震波的走時(shí)計(jì)算和射線追蹤模擬中,前人多是針對(duì)弱各向異性介質(zhì)采用群速度近似的表示方法來(lái)進(jìn)行計(jì)算[16-19],但當(dāng)各向異性強(qiáng)度較大時(shí),這些方法會(huì)造成很大誤差。為此,本文擬對(duì)群速度計(jì)算公式進(jìn)行推導(dǎo),并利用牛頓迭代法快速求解群速度,以有效進(jìn)行地震波在復(fù)雜三維地質(zhì)中的射線追蹤模擬。

1 群速度的計(jì)算

在各向異性介質(zhì)中,由于群速度和相速度發(fā)生分離,波前面不再是以震源為中心的標(biāo)準(zhǔn)球形[20]。在已知彈性參數(shù)的橫向各向同性介質(zhì)中,P波相速度與相角θ的關(guān)系式如下:

(1)

式中:ρ是介質(zhì)的密度;Cij(i,j=1,2,3,4,5,6)是介質(zhì)的彈性參數(shù);D(θ)的表達(dá)式如下:

D2(θ)=(C33-C44)2+2[2(C13+C44)2-

(C33-C44)(C11+C33-2C44)]sin2θ+

[(C11+C33-2C44)2-4(C13+C44)2]sin4θ

(2)

射角?和相角θ的關(guān)系如下:

?=?(θ)=θ+f(θ)

(3)

(4)

其中

C33-2C44)]sinθcosθ+[(C11+C33-

2C44)2-4(C13+C44)2]sin3θcosθ}

(5)

群速度關(guān)于相速度的表達(dá)式為:

[v(θ)secf(θ)]2

(6)

牛頓迭代法是一種在實(shí)數(shù)域和復(fù)數(shù)域上近似求解方程的方法,通過(guò)切線來(lái)逼近零點(diǎn),收斂速度很快。根據(jù)牛頓迭代法基本原理,在VP(θ)、ρ、Cij(i,j=1,2,3,4,5,6)已知的情況下,對(duì)式(1)和式(2)進(jìn)行聯(lián)立變形得:

g(θ)=4A1sin3θcosθ+2A2sinθcosθ

(7)

式中Ai(i=1,2,3)為常數(shù),由VP(θ)、ρ、Cij(i,j=1,2,3,4,5,6)確定。

結(jié)合牛頓迭代公式可得式(8),根據(jù)式(8)可快速求出相角θ的精確近似值,從而得到群速度。

(8)

2 射線追蹤

由于最短路徑法計(jì)算效率快,且在網(wǎng)格密度十分密集的情況下,最短路徑法得到的射線路徑可近似為真實(shí)的理論路徑,因此本文在牛頓迭代法求解群速度的基礎(chǔ)上,通過(guò)最短路徑法來(lái)實(shí)現(xiàn)地震波的射線追蹤模擬?；诘貙拥钠鸱闆r,采用Delaunay三角剖分法對(duì)三維地質(zhì)結(jié)構(gòu)進(jìn)行剖分,其單元結(jié)構(gòu)為不規(guī)則四面體。這樣的不規(guī)則四面體結(jié)構(gòu)不僅能夠更加真實(shí)的展現(xiàn)地質(zhì)結(jié)構(gòu),如突起、凹陷、裂縫等,還能方便地考慮各向異性的特點(diǎn)。

2.1 質(zhì)點(diǎn)設(shè)置

三維地質(zhì)速度結(jié)構(gòu)模型中每一個(gè)四面體單元可以看作為一個(gè)質(zhì)點(diǎn),在這個(gè)單元體中各物性參數(shù),各向異性強(qiáng)度可視為常數(shù),即在該單元體內(nèi)波速值大小一定。質(zhì)點(diǎn)的位置位于四面體的中心,其計(jì)算公式如下:

x=(x1+x2+x3+x4)/4

(9)

y=(y1+y2+y3+y4)/4

(10)

z=(z1+z2+z3+z4)/4

(11)

式中(x,y,z)為質(zhì)點(diǎn)坐標(biāo),(xi,yi,zi)i=(1,2,3,4)為四面體的四個(gè)頂點(diǎn)坐標(biāo)。

每個(gè)質(zhì)點(diǎn)的波速采用反距離加權(quán)法進(jìn)行求解。反距離加權(quán)法是一種應(yīng)用較廣泛的加權(quán)平均插值法,在四面體的四個(gè)頂點(diǎn)中,距離質(zhì)點(diǎn)較近的頂點(diǎn)波速所占權(quán)重較大,距離質(zhì)點(diǎn)較遠(yuǎn)的頂點(diǎn)波速所占的權(quán)重較小。根據(jù)式(12)可求出質(zhì)點(diǎn)的波速。

v=w1v1+w2v2+w3v3+w4v4

(12)

wi=di/D

(13)

式中:v為質(zhì)點(diǎn)速度;vi為四面體四個(gè)頂點(diǎn)的速度值;wi為各頂點(diǎn)所占的權(quán)重;di為各頂點(diǎn)到質(zhì)點(diǎn)的距離;D為四個(gè)頂點(diǎn)到質(zhì)點(diǎn)的距離之和。

2.2 算法原理

最短路徑算法實(shí)際上是一種局部尋找最優(yōu)解法的貪心算法。地震波從震源點(diǎn)出發(fā),先考慮與之相鄰的點(diǎn),在這些節(jié)點(diǎn)中選出走時(shí)最少的一點(diǎn)作為子震源,波會(huì)繼續(xù)傳向與子震源相鄰的節(jié)點(diǎn),在這些節(jié)點(diǎn)中再選出下一個(gè)子震源。按照這樣的思路,重復(fù)循環(huán)就可以快速求出各個(gè)節(jié)點(diǎn)的走時(shí)。在進(jìn)行射線路徑的追蹤時(shí),從接收點(diǎn)開始依次尋找上一級(jí)震源直到震源點(diǎn),記錄所經(jīng)過(guò)的節(jié)點(diǎn)就可得到射線追蹤軌跡。具體步驟如下:

(1)初始化兩個(gè)集合N和V。其中集合N保存所有質(zhì)點(diǎn)中未訪問(wèn)過(guò)的質(zhì)點(diǎn),集合V記錄已經(jīng)訪問(wèn)過(guò)的質(zhì)點(diǎn)。

(2)從震源點(diǎn)開始,尋找與其相連通且未被訪問(wèn)過(guò)的質(zhì)點(diǎn),將這些質(zhì)點(diǎn)存入集合N中。

(3)在N中找出從震源點(diǎn)到該質(zhì)點(diǎn)走時(shí)最小的點(diǎn),將該點(diǎn)視為子震源點(diǎn),存入集合V中,并找出該點(diǎn)的所有子震源點(diǎn)。

(4)遍歷訪問(wèn)子震源點(diǎn)的所有子質(zhì)點(diǎn),并求出從起始震源點(diǎn)到這些子質(zhì)點(diǎn)的走時(shí),把這些子質(zhì)點(diǎn)存入集合N中。

(5)循環(huán)步驟(3)和步驟(4),直到集合N為空。

在最短路徑的循環(huán)遍歷中,需不斷選取走時(shí)最小的點(diǎn),這就需要花費(fèi)時(shí)間,導(dǎo)致算法運(yùn)行時(shí)間較長(zhǎng)。對(duì)于本研究數(shù)據(jù)為無(wú)序序列,可采取快速排序算法,從而提高算法的效率,節(jié)省算法運(yùn)行時(shí)間。

3 實(shí)例應(yīng)用

本文選取的研究區(qū)為華北克拉通山西斷陷帶北部局部區(qū)域,數(shù)據(jù)引自華北地區(qū)地殼-上地幔地震波速度結(jié)構(gòu)模型v2.0,以P波波速進(jìn)行群速度的計(jì)算和射線追蹤模擬。華北克拉通是地球上最古老的克拉通之一,由于復(fù)雜的地質(zhì)結(jié)構(gòu),導(dǎo)致該地區(qū)地質(zhì)活動(dòng)頻繁,長(zhǎng)期的地質(zhì)演化、應(yīng)力作用、板塊運(yùn)動(dòng)及應(yīng)變等都會(huì)導(dǎo)致該地區(qū)地下速度結(jié)構(gòu)不均勻分布,從而該地區(qū)的三維速度結(jié)構(gòu)中留存著長(zhǎng)期的地質(zhì)構(gòu)造演化的信息[21]。因此,研究該強(qiáng)震區(qū)的地質(zhì)結(jié)構(gòu)具有十分重要的意義。

圖1中紅色矩形框內(nèi)為研究區(qū)域,范圍112°～113°E,40°～42°N,數(shù)據(jù)的分辨率在深度7 km以內(nèi)為0.2°×0.2°×0.2 km;深度為7～50 km時(shí)為0.2°×0.2°×1 km。

圖1 研究區(qū)范圍Fig.1 Study area

基于Paraview平臺(tái),利用其提供的開源腳本,使用Python語(yǔ)言進(jìn)行腳本開發(fā),進(jìn)行該研究區(qū)三維地質(zhì)模型的構(gòu)建。若原始數(shù)據(jù)中已有投影坐標(biāo),則可直接讀取投影坐標(biāo)進(jìn)行模型的構(gòu)建;若無(wú)投影坐標(biāo),則需先導(dǎo)入Proj4投影坐標(biāo)庫(kù),編寫投影坐標(biāo)函數(shù),將地理坐標(biāo)轉(zhuǎn)為投影坐標(biāo)。

對(duì)處理過(guò)的數(shù)據(jù)分別使用vtkDelaunay2D()類和vtkDelaunay3D()類進(jìn)行二維地層速度界面和三維速度結(jié)構(gòu)模型的構(gòu)建,并將二維地層速度界面和三維速度結(jié)構(gòu)模型相結(jié)合,得到三維地殼速度結(jié)構(gòu)模型(圖2)。

圖2 研究區(qū)三維地殼速度結(jié)構(gòu)模型Fig.2 Three dimensional crustal velocity structure model in the study area

圖2中三維地質(zhì)結(jié)構(gòu)模型包括了地球內(nèi)部的主要間斷面(上地殼面、中地殼面和下地殼面),x、y、z軸分別代表三維空間的三個(gè)方向(緯度、經(jīng)度、深度)。由圖可知,研究區(qū)沉積層深度范圍大致在0～4.8 km左右,P波波速在6.1 km/s以下;研究區(qū)上地殼深度范圍在5～15 km左右,P波波速在6.1～6.25 km/s范圍之間;研究區(qū)中地殼深度范圍在15～40 km左右,P波波速在6.25～6.7 km/s范圍之間;研究區(qū)下地殼深度范圍在40 km以上,P波波速在6.7 km/s以上。

利用牛頓迭代法求取群速度,根據(jù)最短路徑法原理,基于華北克拉通山西斷陷帶北部局部區(qū)域三維地質(zhì)模型進(jìn)行射線追蹤。

在進(jìn)行最短路徑射線追蹤時(shí),需改進(jìn)連通性的判斷算法。設(shè)震源點(diǎn)在最低層的一個(gè)單元四面體中,對(duì)模型中所有的四面體進(jìn)行連通性判斷。在進(jìn)行兩個(gè)四面體連通性判斷時(shí),可通過(guò)如下算法進(jìn)行判斷:

(1)在Paraview的腳本中通過(guò)GetCell()方法得到對(duì)應(yīng)單元四面體的索引。

(2)在Paraview的腳本中通過(guò)GetPointId()方法得到對(duì)應(yīng)單元四面體的四個(gè)頂點(diǎn)坐標(biāo)。

(3)如果兩個(gè)四面體都處于同一層地質(zhì)介質(zhì)中,則只需判斷兩四面體的空間點(diǎn)坐標(biāo)是否有三個(gè)共同點(diǎn),若有三個(gè)共同的點(diǎn),則兩個(gè)四面體相連通,反之,則不連通。

(4)如果兩個(gè)四面體分別處于不同的地質(zhì)介質(zhì)中,但滿足條件(3),則需遵循snell定理判斷是否能夠發(fā)生折射,若未達(dá)到臨界角,則兩個(gè)四面體相連通,反之,則不連通。

采用最短路徑算法可求出所有與震源連通的節(jié)點(diǎn)走時(shí),并找出到檢波點(diǎn)用時(shí)最小的路徑。如圖3所示,圖中紫線選中的單元分別為震源單元和檢波單元,紅線即為震源點(diǎn)至檢波點(diǎn)的模擬路徑。

圖3 最短路徑射線追蹤效果圖Fig.3 Effect picture of shortest path ray tracing

4 結(jié)論

本文基于Paraview平臺(tái)自動(dòng)化構(gòu)建三維地質(zhì)模型并進(jìn)行射線追蹤模擬,為三維地質(zhì)模型的構(gòu)建和地震波射線追蹤模擬及可視化提供了一種新思路。根據(jù)地下普遍存在各向異性的事實(shí)和地震波基本傳播規(guī)律,利用牛頓迭代法高效求解群速度,并以華北克拉通山西斷陷帶北部局部區(qū)域?yàn)槔?基于研究區(qū)三維地質(zhì)模型采用最短路徑法進(jìn)行地震波射線追蹤模擬及可視化。結(jié)果表明,該方法減少了由各向異性對(duì)地震波傳播帶來(lái)的影響,清晰表達(dá)了研究區(qū)地質(zhì)結(jié)構(gòu)和各向異性特點(diǎn)。

本文方法克服了目前基于二維平面、三維剖面或簡(jiǎn)單三維空間的各向異性速度結(jié)構(gòu)研究難點(diǎn),實(shí)現(xiàn)地震波在真實(shí)地下復(fù)雜三維空間傳播規(guī)律的模擬,充分表達(dá)地下速度不均勻性和各向異性,更好地幫助研究人員解讀研究區(qū)地質(zhì)結(jié)構(gòu)。

致謝:本研究使用的華北克拉通中部山西斷陷帶區(qū)域速度結(jié)構(gòu)模型數(shù)據(jù)引自華北地區(qū)地殼-上地幔地震波速度結(jié)構(gòu)模型v2.0,(鄭天愉、段永紅、許衛(wèi)衛(wèi)、艾印雙、陳凌、趙亮、張耀陽(yáng)、徐小兵,2015),數(shù)據(jù)鏈接網(wǎng)址:http://www.craton.cn/data,在此表示感謝。