四面體網格剖分下速度與反射界面的同時反演

何雷宇 嚴 星 白超英*③

(①長安大學地質工程與測繪學院地球物理系，陜西西安 710054； ②新疆財經大學網絡與實驗教學中心，新疆烏魯木齊 830011； ③長安大學計算地球物理研究所，陜西西安 710054)

1 引言

層析成像始于醫(yī)學中的CT技術，Aki等[1]首次將其應用于地震學研究中。經過四十多年的發(fā)展和完善，地震層析成像已成為研究地下地質結構行之有效的主要技術途徑之一，并且涌現出多種成像的方法技術[2]。

傳統(tǒng)射線類走時成像方法一般采用規(guī)則網格單元(2D中的矩形網格單元、3D中的立方體網格單元)模型參數化方式，這種簡單模型參數化可以很容易地通過網格節(jié)點的編號確定該節(jié)點的空間坐標，從而節(jié)約計算機內存空間。然而這種簡單的模型參數化形式也存在諸多問題。例如：若要對復雜區(qū)域進行較為精確的描述，則必須進行網格單元的精細剖分，從而增加計算時間。更為重要的是規(guī)則網格單元參數化對起伏地表及地下不規(guī)則波阻抗界面無法進行準確的擬合，從而造成正演中不必要的計算誤差，進而導致反演中出現假象。

另一方面，從反演角度來說，規(guī)則網格的精細剖分，意味著更多的未知模型參數參與反演，從理論上就需要更多的射線進行覆蓋，故而求解的方程組的維數也將大大增加，數據對模型的約束也就越差。此外，更多的未知參數意味著參數矢量中欠定的和位于零空間的分量數目更多，性態(tài)大大變壞； 層析反演難度的增加，通常需加入正則化約束以改變方程組的病態(tài)，而這又增加了存儲量和計算量[3]。

上述問題的有效解決方案之一是采用不規(guī)則網格單元(包括三角網格或四面體網格)進行模型參數化。例如可根據模型中不同區(qū)域物性的復雜程度選擇不規(guī)則網格單元的大小，同時不規(guī)則網格單元剖分對起伏地表和地下波阻抗界面(含不規(guī)則異常體)的描述更加精確，從而確保了正演所需的計算精度。反演中由于未知模型參數相對較少，則所需射線也相對較少，故而待求解的方程組維數也大大減小，從總體上減少了反演的存儲量[3]。

三維情況下使用四面體網格單元的另一個優(yōu)勢在于：實際中地震臺站的分布是不均勻的，如北京圈附近的臺站密度與我國西部地區(qū)的臺站密度、陸地的臺站密度與海洋區(qū)域的臺站密度是不同的，加之天然地震的發(fā)生本身就具有不均勻性，兩者結合造成了地震數據覆蓋密度的極度差異性，所以有必要對數據密度大的區(qū)域進行更細致的剖分以增加數據的利用效率。

不規(guī)則網格單元在地震數值模擬中有較為廣泛的應用。正演方面，如在波場模擬中使用的不規(guī)則網格[4，5]、變網格[6]和貼體網格[7]等。而在地震射線追蹤研究方面，李波濤等[8]將三角剖分與波前重建法相結合，實現了初至波射線路徑和走時的快速計算； Yu等[9]引入單元密度控制函數，實現了三角網格下初至波射線追蹤的全局算法； Lan等[10]借助于流體力學中的貼體網格和坐標變換，引入Lax-Friedrichs快速掃描法求解與地形有關的程函方程以獲取地震波初至走時場。

上述幾位學者的研究主要是解決了2D模型中不規(guī)則網格單元模型參數化下初至波走時的計算與射線追蹤問題。Vinje等[11]在開放模型的波前構造法中將三角網格應用于反射界面的描述和波前面的描述，并實現了反射波射線追蹤。趙瑞等[12]將梯形網格與矩形網格結合，實現了多震相射線追蹤； Bai等[13]提出一種基于三角(2D)或四面體(3D)單元剖分下的最短路徑射線追蹤算法(Triangular Shortest-path Method，TSPM)，并實現了多震相地震射線路徑的追蹤與走時計算； 李曉玲等[14]在趙瑞等[12]研究的基礎上，將混合網格下的射線追蹤算法推廣至各向異性TI介質中。但混合網格僅可模擬起伏地表，無法對模型內部進行差異化剖分，優(yōu)勢有限； 隨后，李興旺等[15]實現了四面體網格參數化下TI介質中的多震相射線追蹤。

基于不規(guī)則網格單元剖分下的層析成像研究最早可追溯到20世紀80年代。Tarantola等[16]意識到在地震層析成像中使用固定的均勻網格單元的一些局限性并提出了block-less模型參數化方法； Sambridge等[17]率先將Delaunay四面體和Voronoi多面體應用于地震層析成像領域； Gudmundsson等[18]利用Delaunay四面體和Voronoi多面體進行模型參數化，使用初至P波層析成像改進了區(qū)域上地幔地球參考模型； Curtis等[19]通過對給定的一組點采用Delaunay三角化來描述模型的速度結構，進而研究了2D情況下的井間反演問題； B?hm等[20]用Delaunay三角形和Voronoi多面體在3D反射層析成像中發(fā)展了一種自適應反演方法； Sambridge等[21]使用自適應Delaunay四面體網格單元進行了全球P波層析成像研究。

中國國內對不規(guī)則網格層析成像的研究較少。成谷等[3]對比了三角網格層析成像與矩形網格層析成像的優(yōu)缺點； 于師建等[22]采用三角網格模型參數化進行走時層析成像研究。他們所做的研究主要是不規(guī)則網格參數化下初至波走時層析成像問題，但并未引進反射波。鑒于此，本文研究了2D三角網格模型參數化的多震相聯(lián)合反演及速度與反射界面同時反演問題[23]。為了更好地解決實際地震工作的需要，還將該算法進一步推廣至3D四面體網格。在此過程中，要解決的主要問題是如何進行模型參數化、用參數化后的網格進行射線追蹤及如何求解Jacob矩陣元素并進行合理分配。顯然，3D情況下算法實現的難度總體上大于2D情況。

文中基于Multistage TSPM正演算法[13]，結合共軛梯度求解帶約束的阻尼最小二乘解反演算法，同時得到了四面體網格單元下走時關于速度、走時關于反射界面的偏導計算公式； 討論并實現了四面體網格單元模型參數化下3D復雜模型中地震多震相走時同時反演成像技術。本文主要研究內容包括： ①反射界面已知時地震多震相走時聯(lián)合成像技術； ②多震相走時聯(lián)合同時反演速度模型和反射界面的技術。數值模擬實驗結果表明，四面體網格單元模型參數化下的多波走時成像及同時反演成像結果較為可信，可作為起伏復雜介質中速度場和反射界面重建的一種有效方法。

2 多震相射線追蹤算法

2.1 模型參數化

本文采用的網格由開源的Tetgen生成[24]。利用Multistage TSPM進行射線追蹤之前，首先要對模型進行四面體網格單元參數化(圖1)。采用四面體單元進行模型參數化時，定義四面體單元的四個角點為主節(jié)點。為了增加射線出射角的覆蓋率(即保證正演的計算精度)，在四面體單元的四條邊上等間距插入次級節(jié)點，四條邊插值結束后再對四個面進行插值。四面體網格單元內部沒有節(jié)點，但炮點和檢波點可置于其中(圖1)。模型中速度的采樣是在主節(jié)點上進行，而次級節(jié)點上的速度則是通過速度插值得到(詳見下文)。

圖1 速度插值示意圖

2.2 最小走時的計算和速度插值

從炮點開始計算，節(jié)點走時的計算公式為(兩點間距離除以平均速度)

(1)

式中：i是波前點中當前的次級震源；j為待計算任意節(jié)點；Nj為待計算節(jié)點的集合；ti為炮點到達i節(jié)點射線的最小走時；D(xi,xj)是節(jié)點i與j之間的距離；v(xi)、v(xj)分別是i和j節(jié)點上的速度。當i節(jié)點不是主節(jié)點時，可通過下式求取該節(jié)點速度

(2)

式中：vj表示所在單元各主節(jié)點的速度，在三維模型中k=4；uj為點在四面體中的體積坐標，例如圖1中P點在四面體T1T2T3T4內的體積坐標uj可表示為

(3)

式中V(PT2T3T4)、V(T1T2T3T4)等表示各四面體的體積。顯然，四點共面時體積為零。

2.3 分區(qū)多步計算原理

完成模型參數化和速度插值之后，為了減少計算量，避免對射線未經過區(qū)域的計算，可根據反射界面的起伏形態(tài)將模型分成若干(層狀)區(qū)域。從震源所在區(qū)域開始計算，根據給定的震相種類對特定射線進行追蹤并計算節(jié)點的最小走時。以圖2模型為例，該模型共有上下兩個界面，將整個模型分為3個區(qū)域，假設震源和臺站均位于地表，現在需追蹤第一個反射界面的反射轉換波P1S1(數字1表示計算區(qū)域編號，上標代表上行波，下標代表下行波)。計算步驟為： ①選擇第一區(qū)為計算區(qū)域，調用P波速度模型從震源開始計算直至一區(qū)所有網格全部計算完畢； ②保留一界面上的走時和射線路徑，同時選擇一區(qū)為計算區(qū)域，從界面一上搜索出走時最小點并將該點作為次級震源點，調用S波的速度模型，計算該區(qū)域的走時和射線路徑。按照上述原理可實現多震相地震射線的追蹤計算。

圖2 模型示意圖

3 反演算法

3.1 反演目標函數

在速度與反射界面同時反演中，走時對界面深度的偏導數遠小于走時對速度的偏導數，導致同時反演中界面基本不會更新。因此本文采用不同參數歸一化反演方法。多震相走時同時反演問題可歸結為帶約束的阻尼最小二乘最優(yōu)化問題，其目標函數為

(4)

(μCm+ATCdA)Zmm=ATCdd

(5)

式中Cm、Cd為模型和數據空間的協(xié)方差矩陣。該式是對非線性問題局部線性化的基本反演公式，其解具有局域解的特征。為了得到具有實際物理意義的解，可采用Zm、Cm及Cd的先驗信息。因此，選擇不同組合的(Cm，Cd)可得到不同形式的反演解。若μ=0,Cd=I則式(5)變?yōu)榻浀涞淖钚《藛栴}；若Cd=Cm=I，Zm=W，A=GW，這里W是表征射線寬度的帶寬矩陣，則式(5)變?yōu)镸eyerholtz等[26]提出的卷積壓制法。本文主要采用Tarantola等[27]提出的廣義帶約束的阻尼最小二乘反演解，即Cm、Cd分別為模型和數據空間的協(xié)方差矩陣的逆矩陣，其解為

(6)

該解的約束條件是

式(4)可用迭代的共軛梯度法進行求解[28，29]，其關鍵是如何求解四面體單元中具有偏導性質的Jacobi矩陣中的元素。

3.2 Jacobi矩陣的求解

根據射線理論可知，三維情況下沿某一射線路徑Rj的走時tj可用下式表示

(7)

式中：vc(x,y,z)為四面體單元內速度函數； ds為射線的微分元； 模型參數化后，射線路徑Rj上的走時可寫成求和形式

(8)

式中：n為該射線所穿過的四面體單元總數；Rj,k表示穿過第k個四面體單元的射線長度；vc,k(x,y,z)表示該四面體單元的速度分布。

速度與界面同時反演時的Jacobi矩陣包括兩部分：走時對速度變化的偏導數和走時對反射點深度變化的偏導數

(9)

式中：vk是第k個四面體單元的速度分布；tj是第j條射線的走時；zk是第k個反射點的深度值。

射線穿過某個四面體單元時，其走時對速度的一階偏導可用通過該四面體單元射線兩端點的平均值來代替[30]

(10)

其中

式中：ki和kj是穿過第k個四面體單元內射線的兩個端點；vb(k)是第k個四面體單元內四個主節(jié)點的速度值。

(11)

圖3 穿過某一個四面體網格的射線示意圖

式(10)中，界面p點處走時關于界面深度的導數為[20]

(12)

式中：γ1為入射向量與界面法向量的夾角；γ2為出射向量與界面法向量的夾角；v1,v2分別為入射向量側和出射向量側的速度。實際計算時，只需記錄入射點，反射點和出射點便可得到入射向量和出射向量，再求取在反射界面上反射點處的法向量即可得到走時對界面深度的偏導。

本文采用式(11)計算走時對速度的偏導，用式(12)計算走時對界面深度的偏導。在速度與界面的同時反演中，為了避免因反射界面更新引起某些區(qū)域的過度(或欠)更新，筆者對兩種不同類型的Jacobi矩陣元素進行了歸一化處理。

4 數值模擬分析

4.1 計算效率分析

已知反演的耗時中正演的占比超過90%，所以同時反演的計算效率主要取決于選用的射線追蹤方法。為此對比了立方體模型參數化和四面體模型參數化下單炮正演的CPU計算時間，所選立方體模型的尺寸為80km×80km×48km。將主節(jié)點間距固定為8.0km，將次級節(jié)點間距從4.0km到0.25km逐步減半縮小，并記錄不同模型參數化下的計算時間(所用計算機：Dell，Intel(R) Core(TM)2 Duo CPU E7200，主頻2.52GHz)。

一般來說，無論采用何種模型參數化方式，當主節(jié)點的間距固定時，其正演的計算效率主要取決于次級節(jié)點的間距，但是隨著次級節(jié)點間距的減小，兩種模型參數化下的CPU時間差將大幅減小(表1，CPU時間差從次級節(jié)點間距為4.0km時的近200倍降到了次級節(jié)點間距為0.25km時的不到2倍)。從總節(jié)點數和CPU時間的關系圖(圖4)中也可得出類似結論。兩種模型參數化下的CPU時間都與次級節(jié)點間距成近似指數型的關系，但是立方體模型曲線的斜率更大。隨著總節(jié)點數目的增加，兩條曲線最后完全發(fā)散。

從圖4和表1可見，在需要高精度計算的情況(即更多的次級節(jié)點)下，上述兩種模型參數化的計算效率處于同一水平。由于CPU時間中90%以上都被正演所消耗，因此這個結論也同樣適用于后面的同時反演部分。

圖4 CPU時間隨節(jié)點總數的變化關系曲線

次級節(jié)點間距km立方體模型節(jié)點總數/個四面體模型節(jié)點總數/個立方體模型計算時間/s四面體模型計算時間/s4.00 5135 23848 0.03 5.422.0025827864740.116.031.001156913503751.4812.030.50489339136527524.7079.970.2520123155471683390.33723.03

4.2 多震相走時聯(lián)合反演及同時反演

為了驗證基于四面體網格單元下的反演算法的有效性，我們選擇了一個比較復雜的速度模型(圖2，模型尺度為80km×80km×48km)。地表起伏約為2km，共計3757個主節(jié)點。次級節(jié)點間距為1.0km，總共有613928個節(jié)點。共剖分出18432個四面體網格單元。在z=-30km附近有一個用正弦函數生成起伏界面，最大起伏為2km； 在z=-46km附近有一個上拱的弧形界面。30個震源(圖2中的灰色圓球)隨機分布在10～20km的深度范圍內，289個檢波器(圖2中的倒三角形)均勻分布在地表(注：以下所有數值模擬實驗中均采用此炮—檢排列)。圖5給出了單炮情況下的射線路徑分布圖。真實速度模型如圖6所示，可見水平面及垂直剖面上速度分布特征。下面討論初至波反演、初至波與反射波聯(lián)合反演、速度與反射界面同時反演。

4.2.1 多震相走時聯(lián)合反演成像

聯(lián)合反演中初始速度模型選用線性增加模型，反射界面為真實界面。反演結果均以與圖6相同的形式給出(圖7為水平切片，圖8為垂直切片)。其中圖7b和圖8b為僅用初至P波的反演結果，圖7c和圖8c為初至P波和P1P反射波聯(lián)合反演的結果，圖7d和圖8d為初至P波、P1P反射波和P2P反射波聯(lián)合反演的結果，圖9 (水平切片)和圖10(垂直切片)分別為反演結果與真實模型的相對誤差。從圖7～圖10可以看出，僅用初至P波進行成像在淺層可大體恢復速度場的起伏形狀，但其成像深度卻明顯不及聯(lián)合反演。在加入P1P反射波之后，由于增加了射線交錯概率，同時又提高了射線密度，從水平切片結果(圖7c)來看，聯(lián)合反演已基本能準確地反演速度場的異常圖像和異常幅度；從垂直切片的結果(圖8c)來看，其成像深度明顯增加；從圖9c來看，模型邊緣的速度場起伏的恢復效果也比初至P波的(圖9b)有明顯提升。加入P2P反射波之后，從兩個切片的結果(圖7d和圖8d)來看，速度場的恢復與真實模型已基本相差無幾，且從圖9d和圖10d可見，三種波聯(lián)合反演的精度最高。此數值模擬實驗表明，基于四面體網格的多震相走時聯(lián)合反演可更為準確地反演出異常體形狀和異常幅度，即效果明顯優(yōu)于單一震相的反演。

圖5 射線路徑示意圖

圖6 三維速度模型

4.2.2 多震相走時速度和界面同時反演成像

在速度與界面的同時反演中采用的初始速度模型為向下線性增加模型，實際中初始反射界面可通過一些先驗信息獲得。這里初始反射界面我們選為水平界面(上界面為z=-30km，下界面為z=-47km)。對兩種不同的Jacobi偏導元素進行了最大值歸一化處理。圖11a和圖12a為初至P波和P1P反射波聯(lián)合同時反演的速度結果，圖11b和圖12b為初至P波、P1P和P2P聯(lián)合同時反演的速度結果，圖11c、圖11d、圖12c、圖12d為對應的相對誤差。圖13b為初至波和P1P反射波聯(lián)合同時反演的上界面結果，圖13c、圖13e分別為初至波、P1P和P2P聯(lián)合同時反演的上、下界面的反演結果。

從速度反演結果的水平切片來看，同時反演已大體恢復速度場的起伏形態(tài)，由于模型邊界附近相比模型中央區(qū)域數據覆蓋密度要小，而速度與界面的耦合關系又決定了界面的恢復情況必將影響速度的恢復，所以在水平切片(圖11)中可看出模型中央速度恢復較好，但邊緣區(qū)域恢復稍差，有些許失真，這種失真現象在雙界面的同時反演(圖11b、圖11d)中由于較多的數據覆蓋而使之得到了一些改善。

圖7 不同震相走時聯(lián)合反演成像的水平切片(z=-16km)

圖8 不同震相走時聯(lián)合反演成像的垂直剖面(y=40km)

從垂直剖面(圖12)中可看出，單界面同時反演時-30km以上的速度場的起伏已經基本恢復，但在模型的邊緣有些許失真(圖12a、圖12c)。雙界面的同時反演(圖12b、圖12d)已經基本恢復速度場的起伏形態(tài)，相比于單界面的情況，模型邊緣及下方速度場的失真情況也有了明顯改善。

從圖13中可見，無論是初至P波與P1P波同時反演，還是初至P波與P1P、P2P波同時反演，上界面的形態(tài)都已基本恢復。由于在雙界面的同時反演中模型邊界附近的數據密度大于單一界面同時反演情況下的數據密度，而大數據密度意味著對速度和界面更好的約束，也意味著更好的反演結果，這一點可從圖13b與圖13c的對比中看出，無論是幅值的恢復情況，還是界面形態(tài)的相似度，圖13c的結果都要比圖13b好。從雙界面的同時反演結果中已基本可見下界面近似弧形的形態(tài)特征，但其幅值的恢復不如上界面，這是因為在上界面之上有直達波信息，可幫助更加準確地反演出速度分布，進而更好地反演上界面的形態(tài)，而處在上界面和下界面之間的區(qū)域則沒有這個優(yōu)勢，所以在雙界面的同時反演中，上界面的反演結果優(yōu)于下界面。從這個數值模擬實驗中可得出以下結論：基于四面體網格的同時反演得到的速度和界面結果較為可信。

圖9 反演水平切片與真實模型的相對誤差(z=-16km)

圖10 反演垂直剖面和真實模型的相對誤差(y=40km)

圖11 不同震相走時聯(lián)合同時反演成像水平切片及相對誤差(z=-16km)

5 結果分析與討論

本文基于四面體模型參數化下的多震相射線追蹤技術，結合共軛梯度法求解帶約束的阻尼最小二乘問題反演算法，討論并實現了四面體模型參數化下多震相走時聯(lián)合及同時反演技術。此外還對四面體模型參數化下的計算精度和CPU時間做了詳細討論，可知在總節(jié)點數目大致相同的情況下，其CPU時間遠小于立方體模型下的計算時間。在同時反演中為了平衡速度和界面對走時的影響，本文對兩種不同類型的Jacobi元素進行了最大值歸一化處理。數值模擬實驗表明，四面體模型參數化下的多震相走時聯(lián)合成像技術具有很好的反演能力，能較準確地反演速度場的起伏特征和分布情況；對于速度與界面的同時反演，可較好地反演速度場的起伏特征，并可基本恢復反射界面的幾何形態(tài)。在現有的地震臺站分布極不均勻的情況形下，四面體網格單元可對數據密度不同的區(qū)域進行不同尺度的剖分，進而在走時反演中更顯著地提高數據的利用效率，且四面體網格單元對起伏地表和異常體的擬合精度明顯優(yōu)于立方體單元網格，故該算法具有廣泛的應用價值。