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

張定文,李衛東,段金龍,張學海

(1.河南工業大學信息科學與工程學院,河南 鄭州 450001;2.應急管理部國家自然災害防治研究院,北京 100085)

0 引言

研究表明,各向異性在地下介質中是普遍存在的[1-2]。傳統的地下構造研究多是在各向同性或完全彈性的介質中進行,忽略了由各向異性導致的群速度計算誤差而無法真實表達地下構造。隨著計算機技術,射線追蹤技術等發展,基于復雜各向異性地質模型的地下構造研究成為熱點[3-7]。

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

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

1 群速度的計算

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

(1)

式中:ρ是介質的密度;Cij(i,j=1,2,3,4,5,6)是介質的彈性參數;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)

射角?和相角θ的關系如下:

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

(3)

(4)

其中

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

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

(5)

群速度關于相速度的表達式為:

[v(θ)secf(θ)]2

(6)

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

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

(7)

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

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

(8)

2 射線追蹤

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

2.1 質點設置

三維地質速度結構模型中每一個四面體單元可以看作為一個質點,在這個單元體中各物性參數,各向異性強度可視為常數,即在該單元體內波速值大小一定。質點的位置位于四面體的中心,其計算公式如下:

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

(9)

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

(10)

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

(11)

式中(x,y,z)為質點坐標,(xi,yi,zi)i=(1,2,3,4)為四面體的四個頂點坐標。

每個質點的波速采用反距離加權法進行求解。反距離加權法是一種應用較廣泛的加權平均插值法,在四面體的四個頂點中,距離質點較近的頂點波速所占權重較大,距離質點較遠的頂點波速所占的權重較小。根據式(12)可求出質點的波速。

v=w1v1+w2v2+w3v3+w4v4

(12)

wi=di/D

(13)

式中:v為質點速度;vi為四面體四個頂點的速度值;wi為各頂點所占的權重;di為各頂點到質點的距離;D為四個頂點到質點的距離之和。

2.2 算法原理

最短路徑算法實際上是一種局部尋找最優解法的貪心算法。地震波從震源點出發,先考慮與之相鄰的點,在這些節點中選出走時最少的一點作為子震源,波會繼續傳向與子震源相鄰的節點,在這些節點中再選出下一個子震源。按照這樣的思路,重復循環就可以快速求出各個節點的走時。在進行射線路徑的追蹤時,從接收點開始依次尋找上一級震源直到震源點,記錄所經過的節點就可得到射線追蹤軌跡。具體步驟如下:

(1)初始化兩個集合N和V。其中集合N保存所有質點中未訪問過的質點,集合V記錄已經訪問過的質點。

(2)從震源點開始,尋找與其相連通且未被訪問過的質點,將這些質點存入集合N中。

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

(4)遍歷訪問子震源點的所有子質點,并求出從起始震源點到這些子質點的走時,把這些子質點存入集合N中。

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

在最短路徑的循環遍歷中,需不斷選取走時最小的點,這就需要花費時間,導致算法運行時間較長。對于本研究數據為無序序列,可采取快速排序算法,從而提高算法的效率,節省算法運行時間。

3 實例應用

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

圖1中紅色矩形框內為研究區域,范圍112°～113°E,40°～42°N,數據的分辨率在深度7 km以內為0.2°×0.2°×0.2 km;深度為7～50 km時為0.2°×0.2°×1 km。

圖1 研究區范圍Fig.1 Study area

基于Paraview平臺,利用其提供的開源腳本,使用Python語言進行腳本開發,進行該研究區三維地質模型的構建。若原始數據中已有投影坐標,則可直接讀取投影坐標進行模型的構建;若無投影坐標,則需先導入Proj4投影坐標庫,編寫投影坐標函數,將地理坐標轉為投影坐標。

對處理過的數據分別使用vtkDelaunay2D()類和vtkDelaunay3D()類進行二維地層速度界面和三維速度結構模型的構建,并將二維地層速度界面和三維速度結構模型相結合,得到三維地殼速度結構模型(圖2)。

圖2 研究區三維地殼速度結構模型Fig.2 Three dimensional crustal velocity structure model in the study area

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

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

在進行最短路徑射線追蹤時,需改進連通性的判斷算法。設震源點在最低層的一個單元四面體中,對模型中所有的四面體進行連通性判斷。在進行兩個四面體連通性判斷時,可通過如下算法進行判斷:

(1)在Paraview的腳本中通過GetCell()方法得到對應單元四面體的索引。

(2)在Paraview的腳本中通過GetPointId()方法得到對應單元四面體的四個頂點坐標。

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

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

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

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

4 結論

本文基于Paraview平臺自動化構建三維地質模型并進行射線追蹤模擬,為三維地質模型的構建和地震波射線追蹤模擬及可視化提供了一種新思路。根據地下普遍存在各向異性的事實和地震波基本傳播規律,利用牛頓迭代法高效求解群速度,并以華北克拉通山西斷陷帶北部局部區域為例,基于研究區三維地質模型采用最短路徑法進行地震波射線追蹤模擬及可視化。結果表明,該方法減少了由各向異性對地震波傳播帶來的影響,清晰表達了研究區地質結構和各向異性特點。

本文方法克服了目前基于二維平面、三維剖面或簡單三維空間的各向異性速度結構研究難點,實現地震波在真實地下復雜三維空間傳播規律的模擬,充分表達地下速度不均勻性和各向異性,更好地幫助研究人員解讀研究區地質結構。

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