水平層狀VTI介質高效射線追蹤計算方法

2022-08-31 12:50:22蔣星達張偉陳曉非

地球物理學報 2022年9期

蔣星達，張偉，陳曉非

1 南方海洋科學與工程廣東省實驗室(廣州)，廣州 511458 2 南方科技大學深圳市深遠海油氣勘探技術重點實驗室，廣東深圳 518055 3 南方科技大學地球與空間科學系，廣東深圳 518055

0 引言

地下巖石受構造應力以及礦物定向排列等因素影響形成的周期性沉積、薄互層、定向裂縫等，可以表現出明顯的各向異性特征(Johnston and Christensen，1995；Helbig and Thomsen，2005).隨著油氣勘探對成像的精度要求提高，地下介質的各向異性性質越來越引起關注(Tsvankin et al.，2010；Thomsen，2014).以垂直方向為對稱軸的橫向各向同性層狀VTI介質已經被廣泛應用于地震疊前時間偏移成像(Douma and De Hoop，2006；黃中玉等，2009；Jenner，2011)和非常規油氣微地震監測(Li et al.，2013，2014；Perez et al.，2016；Michel and Tsvankin，2017；Gajek and Malinowski，2021；Kim et al.，2021)中.以上成像處理過程中均需要準確計算VTI介質中的地震波走時和射線路徑.各向異性介質中地震波能量以群速度方向(群角)傳播，但是程函方程卻是以相速度角度(相角)為自變量，導致各向異性介質中發展的各種射線追蹤算法需要搜索群角對應的相角，影響計算效率(趙愛華等，2006；Sadri and Riahi，2010；Bai et al.，2013；Wang，2014；趙后越和張美根，2014；楊穎航等，2018；Li et al.，2020).探索適用于疊前時間偏移和微地震監測的層狀VTI介質高效射線追蹤技術對提高油氣勘探數據處理效率和微地震實時監測具有重要意義.

基于初值問題的打靶法是經典的射線追蹤算法.打靶法通過二分法修改射線出射角度，直到能夠準確獲得震源到檢波點的射線路徑.打靶法原理簡單，易于實現，在各向同性介質走時計算中獲得廣泛應用(徐濤等，2004).打靶法的一個難點在于長偏移距時出射角度的微小變化可能導致射線路徑出現巨大擾動.田玥和陳曉非(2005)以射線在高速層位中的水平慢度代替射線出射角度有效的解決了這個問題，Fang和Chen(2019)進一步將該方法擴展到具有速度梯度的一維速度模型，但目前該方法還未見到在強各向異性介質中的應用(Perez et al.，2016).最短路徑法基于圖論，利用網格節點表示速度模型，通過計算鄰近節點之間的走時信息選擇連接震源點到檢波點的最小走時節點連線作為射線路徑(Nakanishi and Yamaguchi，1986；Moser，1991).最短路徑法適合于復雜介質，可以同時追蹤多個檢波點的初至波、反射波、折射波等各種震相走時，在各向同性和各向異性介質中均獲得大量的應用(Fischer and Lees，1993；劉洪等，1995；張美根等，2006a；唐小平和白超英，2009；Bai et al.，2009，2010；趙后越和張美根，2014；Li et al.，2020).最短路徑法的計算結果精度依賴于網格密度，節點連線形成的入射角和出射角不能連續變化，當網格稀疏時，會產生不滿足Snell定律的“之”字型射線路徑.提高網格節點的密度可以提升最短路徑方法的計算精度，但會顯著增加計算時間(Dijkstra，1959；Moser，1991；劉洪等，1995；尤建軍等，2006；張美根等，2006a).張美根等(2006b)和黃國嬌等(2020)提出僅在速度分界面處設置網格節點的界面源法，減少了網格節點的個數，提高了計算效率.張建中等(2004)提出動態網格射線追蹤方法，網格節點的走時和射線路徑的計算不受固定節點位置的約束，在不增加節點個數的情況下顯著提升了射線追蹤算法的精度(李培明等，2013；桑運云等，2013，2014).在各向異性介質中，最短路徑法計算節點之間的走時需要獲得兩點連線方向的群速度，計算給定群角的群速度顯著增加了計算時間(白海軍等，2011)，而前人發展出的各種群速度與群角近似公式(Byun et al.，1989；Sena，1991；趙愛華和丁志峰，2005)只適合于弱各向異性介質.

本文綜合打靶法和最短路徑法的優點提出了一種層狀VTI介質高效射線追蹤方法.該方法不需要尋找特定群角對應的群速度，而是以相角出發，從震源點出射一系列射線，基于打靶法計算射線同各層速度邊界的交點位置和走時.將射線交點作為網格節點，基于節點信息計算檢波點走時信息和射線路徑.當檢波點周圍的網格密度不足時，通過自適應加密增加檢波點周圍的網格節點個數，最終利用插值和動態網格方法(張建中等，2004)獲得檢波點的走時和射線路徑.與傳統VTI介質最短路徑方法相比，本方法網格節點設置靈活，射線路徑不受節點位置限制，在內存和CPU占用資源少的基礎上可以快速同時獲得多個檢波點的射線路徑和走時信息.該方法適用于任意強度層狀VTI介質，計算效率高，在需要計算直達波、反射波的微地震監測、垂直地震剖面(Vertical Seismic Profiling, VSP)成像、井間地震成像、疊前偏移等油氣勘探領域具有廣泛的應用價值.

1 方法描述

本文提出了一種根據稀疏網格節點計算水平層狀VTI介質中震源到檢波點射線路徑和走時的高效射線追蹤方法.該方法流程如圖1所示，主要包括網格模型建立和檢波點射線追蹤計算兩個步驟.在建立網格模型過程中，我們根據從震源出射的稀疏相角精確計算VTI介質中的群速度和群角，根據群角方向計算出射射線與各層位的交點作為稀疏網格節點，并獲得各個節點的精確走時.若檢波點周圍網格節點走時差距較大，則在檢波點周圍局部空間內加密網格節點，精確計算各節點的走時信息.在檢波點射線追蹤計算過程中，通過三次樣條插值算法獲得檢波點的走時信息，并利用動態網格法反傳快速計算檢波點到震源的射線路徑信息.

1.1 VTI介質中的群速度和群角

地震波在各向異性介質中傳播時會分裂為qP、qSV、qSH三種波形.VTI介質中相速度與相角具有如下關系(Thomsen, 1986)：

圖1 本文高效射線追蹤方法流程圖Fig.1 The flowchart of our high-efficient raytracing method

(1)

其中[vP,vSV,vSH]表示qP、qSV、qSH波相速度值(圖2)，[α0,β0,ε,δ,γ]是表征VTI介質性質的Thomsen參數.α0為qP波垂向傳播相速度；β0為qS波垂向傳播相速度；ε為qP波水平和垂直方向速度變化率大小，用來度量qP波各向異性強度；γ為qSH波水平和垂直方向速度變化率大小，用來度量qSH波各向異性強度;δ為qP波垂向各向異性變化的快慢程度(Tsvankin, 1997).

圖2 各向異性介質波前示意圖(Thomsen，1986)其中φ為群角，θ為相角.Fig.2 Schematic diagram of the wavefront in anisotropic media (From Thomsen, 1986)φ is group angle, θ is phase angle.

相角和相速度是地震波沿波前面垂直方向變化所展現的表觀性質，實際上的地震波能量是以群速度大小沿著群角方向傳播.VTI介質中，群速度/群角和相速度/相角具有如下關系：

(2)

其中[V,φ]為群速度與群角.將相速度公式(1)代入(2)式，即可獲得qP、qSV、qSH波給定相角時的群速度和群角表達式(具體表達式見附錄A).

由公式(2)可知，強各向異性VTI介質中群速度和群角之間沒有直接表達式，而是通過相角連接.給定相角時可以通過公式(A1)和(A2)計算出群速度和群角，但是如果計算特定群角方向的群速度，則需要先搜索確定相應的相角，才能通過公式(A1)計算出群速度.盡管可以在走時計算前通過掃描相角建立群角和群速度查找表，提高傳統最短路徑方法節點之間走時計算的效率，但建立高分辨率的查找表仍然對計算效率具有明顯影響.

1.2 網格模型建立

與傳統最短路徑法設置固定網格節點不同，本文以從震源按不同相角出射的一系列射線與界面的交點作為網格節點(圖3).

圖3 直達波網格節點分布示意圖其中S代表震源，k代表層位序號.Fig.3 Schematic diagram of grid nodes in direct waveS is source, k is layer number.

1.2.1 確定射線出射相角分布

如果在震源處采用等相角分布出射射線，有可能在遠偏移距處界面節點間距過大.為此，本文采用了一種簡易的不均勻出射相角分布：以震源與相鄰層位上等間距分布的點連線構成的角度作為出射射線相角(圖4)，相應的相角角度為：

圖4 射線出射相角計算示意圖Fig.4 Schematic diagram of the calculation of ray emergence phase angle

(3)

h為震源與鄰近層位的垂直距離，Rn為層位上的某一個點距離震源的水平偏移距.Rn=nΔR，ΔR為震源相鄰層位節點間隔，n為出射相角個數.需要注意的是，以上公式只是為了給出一個非均勻的出射相角分布，真實的射線路徑按照相角對應的群角出射.

1.2.2 射線在層狀VTI介質中傳播

射線在不同層位傳播過程中相角遵循Snell定律，

(4)

其中下標k代表層位序號(圖3)，[θk+1,θk]為同一層位界面的入射和出射相角，[vk+1(θk+1),vk(θk)]為對應的入射和出射相速度，p為射線參數.當入射相角θk+1已知時，根據公式(1)和(4)可以確定相速度vk+1(θk+1)和射線參數p，出射相角θk可以推導獲得(具體推導過程請見附錄B).

對于qP波和qSV波,

(5)

對于qSH波，

(6)

以上公式中的介質參數為k層的值，這里為了表達簡潔沒有顯示下標.在qP波和qSV波中選擇符合0°≤θk≤90°的正實數作為最終的出射相角，進而可以計算該層位射線傳播的群角和群速度(附錄A).當遇到qSV波多值射線問題時，我們選擇最小的走時(即按照最大速度傳播)作為射線到達檢波點時間.

1.2.3 網格節點位置和走時

確定出射相角后可以計算對應的群角以及群速度，進而可以計算該層射線與相鄰層位的交點，確定相鄰層位節點的位置和走時(圖5).該位置和走時可以表示為：

圖5 射線傳播路徑和網格節點位置計算Fig.5 The computation of grid node position and ray propagation path

(7)

1.2.4 局部網格節點加密

在獲得各層網格節點位置后，由于檢波點與網格節點位置基本不會正好重合，因此需要根據網格節點走時和射線路徑計算檢波點的走時和射線路徑信息.當檢波點兩端網格節點走時差距較大時，我們選擇在兩個節點對應的兩條射線之間添加更多的射線，加密檢波點周圍節點的密度.具體步驟表示為：

(1)獲取檢波點兩側最近節點之間的走時差Δt，判斷Δt是否大于規定時間差ΔT；

(3)重新獲取檢波點兩側最近節點走時差Δt′，若Δt′>ΔT，則重復步驟(2)；若Δt′≤ΔT，則停止.

通過以上操作，可以保證檢波點最終的走時誤差在小于ΔT的精度范圍.

1.3 檢波點走時和射線路徑計算

檢波點周圍的網格節點如圖6所示.根據檢波點所在層位的網格節點時間，檢波點D的走時tr,D可以利用三次樣條函數插值獲得：

tr,D=aD+bD(xr,D-xr,A2)+cD(xr,D-xr,A2)2

+dD(xr,D-xr,A2)3，(8)

系數[aD,bD,cD,dD]由A2、B2界面所有節點位置和走時獲得.

圖6 檢波點D的射線路徑示意圖CD為需要計算的射線路徑；A1、A2、B1、B2為檢波點兩端網格節點；r為層位.Fig.6 Schematic diagram of the raypath of receiver DCD is the raypath to be calculated; A1, A2, B1, B2 are grid nodes on both sides of the receiver; r is the layer number.

傳統最短路徑法通常以檢波點同網格節點連接具有最小走時的路徑作為射線路徑，但是存在“之”字形路徑問題.張建中等(2004)提出的動態網絡射線追蹤方法不受固定節點限制，基于Fermat原理確定出不位于網格節點上的射線路徑，本文應用該方法實現VTI介質的更準確的射線路徑計算，其原理可表示為：

根據圖6，已知節點A1、B1的二維坐標和走時分別為[xr+1,A1,zr+1,A1,tr+1,A1]，[xr+1,B1,zr+1,B1,tr+1,B1].在已知檢波點(或射線路徑與層位的交點)D的坐標和走時[xr,D,zr,D,tr,D]和其射線傳播路徑節點C的縱坐標zr+1,C條件下，C點的走時tr+1,C可以通過A1、B1界面節點走時三次樣條插值獲得：

tr+1,C=aC+bC(xr+1,C-xr+1,A1)+cC(xr+1,C

-xr+1,A1)2+dC(xr+1,C-xr+1,A1)3，(9)

系數[aC,bC,cC,dC]由A1、B1界面所有節點位置和走時獲得.

tr+1,C和tr,D關系可以表示為：

(10)

sCD為射線CD方向的群慢度，可根據A1、B1界面所有已知節點的群慢度插值獲得：

sCD=aCD+bCD(xr+1,C-xr+1,A1)+cCD(xr+1,C-xr+1,A1)2

+dCD(xr+1,C-xr+1,A1)3，(11)

系數[aCD,bCD,cCD,dCD]由A1、B1界面所有節點位置和群慢度獲得.

(1)確定檢波點所在的層位以及兩側網格節點，根據公式(9)獲得檢波點的走時信息，計算檢波點與相鄰層位的射線路徑交點；

(2)以該層位射線路徑交點為新的“檢波點”，重復以上步驟(1)獲得下一個層位的射線交點.根據新的射線交點繼續計算，直到到達震源所在層位；

(3)將以上各層獲得的射線交點與震源和檢波點依次連接，獲得該檢波點在VTI模型中的射線路徑信息.

反射波網格節點的確定與直達波相似(圖7)，除卻在反射層位依據反射定律射線沿反射角傳播外，各層節點位置以及走時信息均可以通過公式(7)獲得.進而根據插值方法計算獲得各個檢波點的走時信息和射線路徑.

圖7 反射波網格節點示意圖Fig.7 Schematic diagram of grid nodes in reflected wave

2 數值實驗

為了驗證本文算法的有效性以及實用性，我們設計了一個含高速夾層的強各向異性多層VTI介質模型(圖8).位于地面和井中的震源以及檢波器用來模擬疊前時間偏移和微地震井中監測的數據采集方式.兩個震源的位置分別為(500 m, 1200 m)和(500 m, 0 m).井中檢波器共20個，水平位置為1000 m，垂直方向沿290 m至670 m按20 m等間隔排列.地面共有40個檢波器，按照50 m等間隔均勻水平排列.

圖8 VTI介質模型(a) 震源點(500 m, 1200 m)與井中檢波器； (b) 震源點(500 m, 0 m)與地面檢波器； (c) 各層Thomsen參數大小.Fig.8 VTI medium model(a) Source (500 m, 1200 m) and the receivers in the downhole well; (b) Source (500 m, 0 m) and the receivers on the surface; (c) The Thomsen parameters in each layer.

2.1 單個檢波器走時和射線路徑計算

首先以震源(500 m, 1200 m)到地面單個檢波器(1500 m, 0 m)的走時和射線路徑為例展示本文方法的計算過程和結果(圖9).對震源上下相鄰的1100 m和1500 m層位，以50 m為間距計算出射射線，根據震源到每個射線點的連線計算出射相角，在震源兩側共計算得到74條射線(上下層位各37條).計算這些射線沿著各層傳播的相角，根據相角對應的群角和群速度獲得各層速度界面的網格節點位置和走時.由圖9a所示，射線密度在遠偏移距并未出現明顯擴散，動態網格節點分布較密集.檢波器兩側節點位置為(1471.6 m, 0 m)(1526.9 m, 0 m)，間隔超過50 m.節點走時為565.57 ms和572.90 ms，走時差為7.33 ms.檢波器并未與網格節點重合，走時無法直接獲得(圖9b).為了保持檢波器兩側最近節點之間的走時差小于1.00 ms，計算加密射線條數為8條(7.33/1.00≈8)，在檢波器兩側最近節點之間從震源出發添加8條射線,獲得添加射線的走時和路徑信息(圖9c和圖9d).此時，檢波器兩側節點位置分別為(1497.9 m, 0 m)(1503.8 m, 0 m)，兩側節點走時為569.03 ms和569.81 ms.因為兩側節點走時差為0.78 ms，小于預設精度1.00 ms，停止插入新的射線.通過三次樣條插值方法獲得檢波器的走時為569.30 ms，與真實值569.20 ms相對誤差小于0.1%.通過動態網格法逐層回傳計算獲得檢波器到震源的射線路徑，該射線路徑與檢波器真實射線路徑基本重合.射線路徑和走時的準確計算驗證了本文方法的可行性.

圖9 檢波器(1500 m, 0 m)射線追蹤計算(a) 研究介質中的網格節點； (b) 地面網格節點走時； (c) 檢波器走時計算； (d) 檢波器射線路徑計算. Node:動態網格節點. New:重新插入的射線. True:真實射線路徑和走時. Receiver:計算的檢波器射線路徑和走時.Fig.9 The raytracing calculation of receiver (1500 m, 0 m)(a) The grid nodes in the study medium; (b) The traveltime of the grid nodes on the surface; (c) The calculation of the traveltime of the receiver; (d) The calculation of the raypath of the receiver. Node:The nodes in dynamic networks. New: The new inserted rays. True: The true raypath and traveltime. Receiver: The calculated raypath and traveltime of the receiver.

圖11 不同方法反射波射線路徑和走時對比(a) qP波真實射線路徑與最短路徑法結果對比； (b) qP波真實射線路徑與本文高效射線追蹤法結果對比； (c) 走時對比. True:真實結果. Shortest:傳統最短路徑法結果. High-eff:本文高效射線追蹤方法結果.Fig.11 The comparison of raypath and traveltimes with different methods in reflected wave(a) The comparison of raypath between the true results and the shortest-path method in qP wave; (b) The comparison of raypath between the true results and our high-efficient method in qP wave; (c) The comparison between different traveltimes. True: True result. Shortest: The result of the traditional shortest-path method. High-eff: The result of our high-efficient method.

圖12 最短路徑法和本文射線追蹤方法對比(a) 不同網格間距兩種方法的計算時間和計算誤差對比； (b) 不同偏移距兩種方法的計算時間和計算誤差對比. Time:Shortest:最短路徑法消耗時間. Time:High-eff: 本文高效射線追蹤方法消耗時間. Misfit:Shortest: 最短路徑法計算誤差. Misfit: High-eff: 本文高效射線追蹤方法計算誤差.Fig.12 The comparison between the shortest-path method and our high-effective raytracing method(a) The comparison of calculation time and misfit in different grid spaces between the two methods; (b) The comparison of calculation time and misfit in different offsets between the two methods. Time:Shortest: The calculation time of the shortest-path method. Time:High-eff: The calculation time of our high-effective raytracing method. Misfit:Shortest: The calculation misfit of the shortest-path method. Misfit:High-eff: The calculation misfit of our high-effective raytracing method.

2.2 精度對比

我們通過同VTI介質最短路徑射線追蹤算法結果對比驗證本文算法精度.傳統最短路徑法沿層位界面設置固定節點，節點間隔為1.0 m，依據Dijkstra算法依次計算單個節點與相鄰層位所有節點的走時，選擇其中的最小走時節點作為該節點的地震波射線傳播路徑.由于每次計算節點之間走時需要獲得群角方向對應的群速度，本文按照1.0°等間隔計算了群角在[0°，90°]區間與群速度對應表，兩節點之間走時計算需要通過查表獲得相應的群速度.

2.2.1 井中監測

本測試模擬了微地震井中監測的數據采集方式.震源(500 m, 1200 m)與垂直檢波器的射線路徑和走時如圖10所示(射線路徑以qP波為例，qSV和qSH波類似).相比于最短路徑法，本文射線追蹤算法獲得的射線路徑與準確射線路徑重疊程度更高，在各個層位的交點更加接近.qP、qSV、qSH三種波的走時對比顯示本文射線追蹤算法走時曲線與真實曲線更加接近，計算精度更高.

2.2.2 地面監測

本測試模擬了疊前時間偏移的數據采集方式.地面檢波器接收地面震源(500 m, 0 m)在介質層位深度1000.0 m的反射波射線路徑和走時如圖11所示(以qP波為例).同直達波計算結果相同，本文方法在反射層位的反射點更加接近準確反射點的位置，射線路徑與準確路徑更加相似.本文方法相比于最短路徑法qP、qSV、qSH波的走時精度更高，走時誤差均小于0.1%，顯示本文射線追蹤方法在反射波計算方面同樣具有優越性.

2.3 效率對比

我們對比了最短路徑法和本文射線追蹤方法的計算效率.兩種方法消耗時間和對應的計算精度如圖12所示.圖12a顯示了震源(500 m, 1200 m)到檢波點(1500 m, 0 m)在不同網格間距下最短路徑法和本文射線追蹤方法計算直達波消耗時間和計算誤差對比.從圖中可以明顯看出，隨著節點間距的增大，最短路徑法效率逐漸提高，但是計算誤差也逐漸增大.而本文射線追蹤方法隨著初射間隔的增大(圖4中Rn)，消耗時間變化不明顯，計算誤差也一直保持小于1%.在節點間距為1 m時，最短路徑法與本文方法計算誤差相當，但是計算消耗時間是本文方法的幾十倍.節點間距的增大提高了最短路徑法的計算效率，但是相比本文射線追蹤方法計算誤差明顯增加，在節點間距為10 m時誤差已經超過了10%.圖12a顯示了震源(500 m, 0 m)以1000.0 m為反射層位在不同偏移距下最短路徑法節點間距為1 m時消耗時間和計算精度對比.在較小偏移距時(0～500 m)，最短路徑法網格范圍較小，計算消耗時間較少，同時誤差較低.隨著偏移距的增大，最短路徑法計算網格范圍較大，計算量增加，消耗時間增加，由于網格離散的特點計算走時誤差也逐漸積累.而本文射線追蹤方法隨著偏移距增大計算耗時未出現明顯變化，計算誤差也一直保持在0.1%左右的較高精度，證明本文方法受震源到檢波點的距離影響較小.

由上述兩個例子可以看出，最短路徑法按照1 m的間隔設置網格節點，提高了走時計算精度，但是在研究區域范圍較大時顯著增加了計算量，計算網格節點之間的走時信息需要查找群角和群速度表，降低了計算效率.本文射線追蹤方法得益于稀疏網格和局部節點加密的特點，計算精度受初射間距影響較少，同時不用計算群角和群速度對應關系，在較大網格間距和較遠的偏移距時仍可以通過插值獲得高精度的走時信息.由此可看出本文射線追蹤方法相比最短路徑法計算時間耗費少，效率高，非常適于層狀VTI介質快速射線追蹤正演計算.

3 結論

本文提出一種適用于強各向異性的水平層狀VTI介質高效射線追蹤方法.本方法利用從震源按照不同相角出射的一系列射線與速度界面的交點作為稀疏網格節點，精確計算網格節點的位置和走時信息.在檢波點周圍自適應增加網格節點密度.利用三次樣條插值法獲得檢波點的高精度走時信息，同時基于動態網格思想計算獲得射線路徑信息.本方法吸收了打靶法快速高效的特點，同時避免了長偏移距射線路徑巨大擾動的難題，在借鑒最短路徑法利用網格節點進行多個檢波點射線追蹤的基礎上，避免了高密度網格節點引起的計算效率降低問題.本方法不受網格節點位置限制，可以獲得不位于網格節點的射線路徑，克服了節點稀疏情況下的射線路徑“之”字形缺陷.基于嚴格計算群角與群速度關系，本方法適用于任意強度VTI介質qP、qSV、qSH波射線追蹤.本方法不用多次遍歷群角和群速度關系，在提高計算效率的基礎上可以高精度計算直達波、反射波等射線路徑信息.同傳統最短路徑方法的對比驗證了本文方法走時和射線路徑計算的高精度和高效性.本方法占用硬件資源少，耗時短，非常適合于需要多次快速正演計算的地震疊前時間偏移、井中微地震實時監測等情況.

附錄A 群速度和群角計算公式

將qP、qSV和qSH波的相速度公式(1)代入群速度和群角公式(2)，可以得到任意強度VTI介質群速度/群角和相角的關系：