李琳，劉韜，胡天躍＊

1北京大學(xué)地球與空間科學(xué)學(xué)院，北京 100871

2中國(guó)石化石油勘探開發(fā)研究院，北京 100083

1 引言

有限差分方法（FDM）和有限元方法（FEM）是地震正演模擬中非常重要的兩種方法，其中有限差分方法因?yàn)槠涑绦驅(qū)崿F(xiàn)的便捷性和計(jì)算的高效性在生產(chǎn)上應(yīng)用廣泛，但是傳統(tǒng)有限差分方法算子是基于均勻介質(zhì)中規(guī)則網(wǎng)格條件下的泰勒展式推導(dǎo)，因而在模擬非規(guī)則的邊界問題當(dāng)中存在困難.此外，直接將差分算子應(yīng)用于非均質(zhì)問題中實(shí)際上也會(huì)引入相當(dāng)?shù)挠?jì)算誤差，這一現(xiàn)象在速度間斷面上尤為明顯，不少學(xué)者也對(duì)基于速度間斷面的差分算子進(jìn)行特殊處理（Kummer etal.，1987；Zahradnik etal.，1993；Lisitsa etal.，2010）.相比而言，有限元方法因?yàn)槠渚W(wǎng)格剖分的靈活性以及邊界條件的適應(yīng)性而日漸受工業(yè)生產(chǎn)的青睞，通過近幾十年的發(fā)展，已經(jīng)由最初的一維模型逐漸拓展到三維彈性波模擬（鄧玉瓊和張之力，1990；張劍鋒，1994；張美根等，2002；Zhang etal.，2002；Min etal.，2003；薛東川等，2007）.然而有限元面臨的一個(gè)主要問題還在于它的計(jì)算效率：在有限元計(jì)算過程當(dāng)中涉及到對(duì)質(zhì)量矩陣求逆計(jì)算.通常而言，有限元構(gòu)架下的質(zhì)量矩陣為一個(gè)對(duì)稱的窄帶寬的矩陣，使得質(zhì)量矩陣的求逆非常耗費(fèi)資源且復(fù)雜.為了解決這一問題，有些學(xué)者引入了集中質(zhì)量矩陣技術(shù)將質(zhì)量矩陣強(qiáng)行對(duì)角化以提高計(jì)算效率，但是在這一過程中會(huì)影響整體的計(jì)算精度（De Basabe，2009）.

為了解決有限元的計(jì)算效率以及精度問題，譜元法（SEM）被引入地震領(lǐng)域（Komatitsch and Vilotte，1998；Komatitsch and Tromp，1999；Komatitsch etal.，2001）.譜元法的核心思想是：在單元內(nèi)部求取質(zhì)量矩陣過程中，通過引入一個(gè)高精度的積分公式，并將單元內(nèi)節(jié)點(diǎn)配置于相應(yīng)的積分點(diǎn)之上，結(jié)合Lagrange插值形函數(shù)計(jì)算，最終獲得一個(gè)對(duì)角化的質(zhì)量矩陣.這種方法和集中質(zhì)量矩陣技術(shù)相比其優(yōu)勢(shì)在于，可以使得質(zhì)量矩陣對(duì)角化并且具備更高的計(jì)算精度.譜元方法在計(jì)算效率上的改進(jìn)，以及網(wǎng)格剖分中對(duì)于復(fù)雜結(jié)構(gòu)描述的適應(yīng)性，使得譜元法在短短的十幾年間得到了迅猛的發(fā)展，并應(yīng)用于正演模擬以及偏移成像當(dāng)中.值得注意的是，大部分的譜元方法研究是基于四邊形網(wǎng)格，其原因在于四邊形網(wǎng)格中Lagrange形函數(shù)可以通過內(nèi)積的方法求?。–ohen，2002），即可以方便地從一維Lagrange表達(dá)式直接推導(dǎo).此外更重要的一點(diǎn)是，四邊形網(wǎng)格中應(yīng)用的GLL（Gauss-Lobatto-Legedre）積分公式被證明是具有很高的計(jì)算精度，并且隨著積分點(diǎn)的增加其計(jì)算精度也在增加（Cohen，2002；De Basabe and Sen，2007）.

相比而言，使用三角單元的譜元方法研究和應(yīng)用較少，因?yàn)樵谌切尉W(wǎng)格當(dāng)中，Lagrange形函數(shù)的構(gòu)建比較復(fù)雜，不能直接通過內(nèi)積的方式獲得；更為關(guān)鍵的是，四邊形網(wǎng)格中效果很好的GLL積分公式在三角網(wǎng)格中不再適用（Cohen，2002）.然而這些問題的出現(xiàn)并不代表譜元法不能應(yīng)用在三角形網(wǎng)格中，相反因?yàn)槿蔷W(wǎng)格在模擬復(fù)雜的構(gòu)造問題中具有更高的靈活性，一直以來也有不少學(xué)者致力于三角譜元法的研究，并提出了一系列新的積分點(diǎn)去替代GLL積分公式（Blyth and Pozrikidis，2006；Chen and Babuska，1995；Hesthaven，1998；Taylor etal.，2000；Cohen etal.，2001）.這當(dāng)中最為常用的點(diǎn)包括Fekete點(diǎn)和Cohen點(diǎn)，其中Fekete點(diǎn)的位置是根據(jù)最優(yōu)性插值而提出來的（Taylor etal.，2000），這類節(jié)點(diǎn)在三角網(wǎng)格的邊界上，其配置方式與GLL點(diǎn)是一致的，因此可以很容易實(shí)現(xiàn)三角譜元法和四邊譜元法的融合（Komatitsch etal.，2001）.而Cohen點(diǎn)的配置是根據(jù)數(shù)值積分最優(yōu)的原理進(jìn)行求解.研究表明，對(duì)于N階的Lagrange形函數(shù)，采用Cohen點(diǎn)進(jìn)行積分計(jì)算具有2 N-1階的計(jì)算精度，已經(jīng)基本接近于GLL點(diǎn)的計(jì)算精度（Liu etal.，2012）.因此從計(jì)算質(zhì)量矩陣的精度而言，Cohen積分點(diǎn)的方式將更為精確.Liu等（2012）的文章中通過頻散分析詳細(xì)比較了基于Cohen點(diǎn)的和Fekete點(diǎn)的TSEM的精度，進(jìn)一步論證了Cohen點(diǎn)在精度上的優(yōu)越性.

本文主要研究采用2階Cohen點(diǎn)的三角譜元方法（TSEM）.首先是對(duì)三角譜元法的基本原理進(jìn)行闡述，包括Cohen積分點(diǎn)的位置，Lagrange形函數(shù)的求取，以及質(zhì)量矩陣的對(duì)角化處理.在這之后，本文將對(duì)2階三角譜元法的計(jì)算精度和穩(wěn)定性條件進(jìn)行分析.為了進(jìn)一步描述三角譜元法的計(jì)算特性，我們還將引入3階精度下的三角有限元（TFEM）算子進(jìn)行對(duì)比.分析表明，2階TSEM方法相對(duì)于3階TFEM方法具有更高的精度和穩(wěn)定性條件，并且需要更少數(shù)目的節(jié)點(diǎn)進(jìn)行刻畫，而數(shù)值算例也進(jìn)一步證明了分析結(jié)果的正確性.本文的最后部分將2階TSEM方法應(yīng)用于兩個(gè)具有孔洞結(jié)構(gòu)特征的地質(zhì)模型正演當(dāng)中，并探討孔洞在地震波場(chǎng)當(dāng)中的響應(yīng)特征.其波場(chǎng)的計(jì)算結(jié)果展示了三角譜元方法在地震正演模擬計(jì)算中的有效性及其廣闊的應(yīng)用前景.

2 三角譜元法的基本原理

考慮二維（x-z）的聲波方程，如果不考慮密度在空間上的變化，在進(jìn)行有限元離散之后滿足：

其中U代表波場(chǎng)函數(shù)，t代表時(shí)間，M和K分別代表質(zhì)量矩陣和剛度矩陣，具體的計(jì)算公式為：

其中Me和Ke對(duì)應(yīng)于單元內(nèi)的相應(yīng)矩陣，Se為從單元到全局的映射矩陣，λ為插值形函數(shù)，p和q對(duì)應(yīng)著單元內(nèi)矩陣的角標(biāo)，（l，m）對(duì)應(yīng)著單元內(nèi)部節(jié)點(diǎn)的位置，Ne為單元內(nèi)部節(jié)點(diǎn)的數(shù)目，而c代表單元內(nèi)介質(zhì)的聲波速度.

對(duì)于公式（1），在時(shí)間上采用2階精度的有限差分格式求解，可以得到：

可以看出，在時(shí)間域采用顯式差分格式之后需要對(duì)質(zhì)量矩陣M進(jìn)行求逆處理.在譜元法中，通過將單元內(nèi)節(jié)點(diǎn)配置于積分點(diǎn)之上，并引入Lagrange插值形函數(shù)和積分公式進(jìn)行求解，獲得對(duì)角化的質(zhì)量矩陣：

公式（4）的質(zhì)量矩陣在積分過程中被求和替代，其中ω表示積分求和中的加權(quán)因子；同樣的，p與（l，m）之間存在著線性映射關(guān)系.當(dāng)插值形函數(shù)λ（x，z）為L(zhǎng)agrange函數(shù)時(shí)滿足：

因此求和公式可以使得質(zhì)量矩陣對(duì)角化，其對(duì)角元素的取值取決于加權(quán)因子ω.在四邊形網(wǎng)格當(dāng)中，加權(quán)因子由GLL積分公式確定.但是對(duì)于三角網(wǎng)格的情況，GLL的積分點(diǎn)配置情況和積分公式不再適用；此外，三角網(wǎng)格中Lagrange形函數(shù)不能簡(jiǎn)單地通過內(nèi)積的方式進(jìn)行求取，因而無法直接尋找到一個(gè)顯式的表達(dá)式.為了解決以上問題，本文采用Cohen積分點(diǎn)對(duì)GLL積分點(diǎn)進(jìn)行替代：研究表明使用N階Cohen點(diǎn)計(jì)算質(zhì)量矩陣具有2 N-1階的計(jì)算精度，保證了質(zhì)量矩陣計(jì)算的精確性（Liu etal.，2012）.

圖1展示了2階Cohen積分點(diǎn)的配置方式：每個(gè)單元內(nèi)部具有7個(gè)節(jié)點(diǎn)，其中6個(gè)節(jié)點(diǎn)是均勻配置在三角單元的邊界之上，還有1個(gè)節(jié)點(diǎn)置于三角單元的中心.相比而言，傳統(tǒng)的2階TFEM使用的是6個(gè)節(jié)點(diǎn)離散，而3階TFEM使用的是10個(gè)節(jié)點(diǎn)離散.傳統(tǒng)有限元當(dāng)中的節(jié)點(diǎn)都是均勻配置于三角單元中.對(duì)于高階情況的Cohen積分點(diǎn)（如圖2所示），Cohen點(diǎn)不再均勻地配置在單元內(nèi)部.

三角單元中Lagrange形函數(shù)的求取可以通過基函數(shù)展開進(jìn)行：

圖1 不同三角單元的節(jié)點(diǎn)配置方式（a）2階三角有限元節(jié)點(diǎn)配置，1個(gè)單元內(nèi)部配置6個(gè)節(jié)點(diǎn)，均勻分布；（b）2階三角譜元法節(jié)點(diǎn)配置，單元內(nèi)部具有7個(gè)節(jié)點(diǎn)；（c）3階三角有限元節(jié)點(diǎn)配置，單元內(nèi)共有10個(gè)節(jié)點(diǎn)信息.Fig.1 The different node configurations in the triangular mesh（a）For the 2nd-order FEM scheme，there are 6nodes placed at the vertices；（b）For the 2nd-order Cohen case，there are 7nodes in the element；（c）For the 3rd-order FEM case，there are 10nodes inside the element.

圖2 高階Cohen積分點(diǎn)的配置情況（a）3階Cohen點(diǎn)配置，單元內(nèi)部一共有12個(gè)節(jié)點(diǎn)（Cohen etal.，1995）；（b）4階Cohen點(diǎn)配置，單元內(nèi)一共有18個(gè)節(jié)點(diǎn)（Mulder，1996）；（c）5階Cohen點(diǎn)配置，單元內(nèi)一共有30個(gè)節(jié)點(diǎn)（Chin-Joe-Kong etal.，1999）.Fig.2 The grid configuration for higher order Cohen nodes（a）The 3rd-order Cohen node，there are 12nodes inside the element（Cohen etal.，1995）；（b）The 4th-order Cohen nodes，there are 18nodes for each element（Mulder，1996）；（c）The 5th-order Cohen nodes（Chin-Joe-Kong etal.，1999），there are 30nodes inside.

其中φ（x，z）為基函數(shù)，cij代表φj（x，z）相應(yīng)的系數(shù)，n為單元內(nèi)部節(jié)點(diǎn)的個(gè)數(shù).根據(jù)Lagrange形函數(shù)的定義，以及選擇冪函數(shù)作為基函數(shù)，我們可以構(gòu)建出2階Cohen點(diǎn)下Lagrange形函數(shù)的表達(dá)式：

同樣地，在構(gòu)建過程中利用Vandermonde矩陣可以計(jì)算出權(quán)因子ω（Mercerat etal.，2006）：

在獲取三角網(wǎng)格中Lagrange插值形函數(shù)和相應(yīng)的Cohen積分點(diǎn)和加權(quán)因子之后，可以根據(jù)公式（2）計(jì)算出相應(yīng)的質(zhì)量矩陣和剛度矩陣進(jìn)行正演計(jì)算.

關(guān)于邊界的實(shí)施，在實(shí)際模型的地震聲波模擬中，最為常用的邊界條件是吸收邊界.其中，完全匹配層（PML）吸收邊界在有限差分模擬當(dāng)中應(yīng)用最為廣泛.盡管PML方法大部分情況下是基于有限差分法求解速度應(yīng)力的一階方程，同樣的，該方法可以應(yīng)用到二階情況下的有限元和譜元法，只需在剖分網(wǎng)格完畢之后，對(duì)在邊界上的單元進(jìn)行PML的處理（Komatitsch and Tromp，2003）.

3 頻散分析

有限元方法或者譜元方法作為一種數(shù)值模擬方法，在正演計(jì)算過程中存在著精度的問題.通常而言，對(duì)于固定的網(wǎng)格配置，地震頻率越高的情況下計(jì)算誤差也越大，這種誤差也稱為數(shù)值頻散.常用的頻散分析方法是通過一個(gè)平面波的近似解替換到控制方程中進(jìn)行分析，這一方法在有限差分算子當(dāng)中得到了廣泛的應(yīng)用（Alford etal.，974；Marfurt，1984；Moczo etal.，2000；劉紅偉等，2010）.但是有限元算法的頻散分析要復(fù)雜得多，考慮到有限元的求解算子在不同的網(wǎng)格點(diǎn)上是變化的，需要通過求解特征方程組的形式綜合考慮不同算子的影響（Carcione etal.，2002）.Mullen和Belytschko（1982）最早給出了一階有限元的頻散分析結(jié)果，他們考慮了矩形網(wǎng)格和不同形狀的三角網(wǎng)格，最終的結(jié)論是矩形網(wǎng)格比三角網(wǎng)格具有更高的精度；之后Mulder（1999）分析了一維情況下譜元法的頻散條件，得出在一維情況下譜元法的頻散特性要優(yōu)于常規(guī)有限元法；Cohen等（2001）接著對(duì)低階的四邊形譜元法進(jìn)行了頻散分析，論證了譜元法在使用GLL積分點(diǎn)的高精度性.De Basabe和Sen（2007）通過數(shù)值求解的方法對(duì)高階四邊形譜元法進(jìn)行了分析，并且證明高階譜元法能夠提高計(jì)算精度.Seriani和Oliveira（2008）在他們的研究中證明在彈性波問題中，高階譜元法也將提高計(jì)算精度.在本節(jié)中，我們對(duì)2階三角譜元法進(jìn)行頻散分析，并將分析結(jié)果和三角有限元方法進(jìn)行對(duì)比.

圖3 ‘X’形狀TSEM和TFEM的計(jì)算節(jié)點(diǎn)配置情況（a）2階TSEM方法，全局計(jì)算采用12個(gè)計(jì)算算子；（b）3階TFEM方法，全局計(jì)算采用18個(gè)算子.Fig.3 The grid configuration of the‘X’type triangular mesh（a）The grid configuration for the 2nd-order Cohen nodes.12classes of operators are considered in this case，which have been signed in the panel；（b）The grid configuration for the 3rd-order FEM，where 18classes of operators are considered in this case.

為了分析方便，我們選擇一種‘X’形狀的三角網(wǎng)格，如圖3所示.2階TSEM的單元內(nèi)部一共有7個(gè)節(jié)點(diǎn)，但是組集到全空間之后，全局有12種計(jì)算算子（如圖中所標(biāo)注）.對(duì)于3階TFEM格式，單元內(nèi)部一共有10個(gè)節(jié)點(diǎn)，相應(yīng)的全局算子為18種.

假設(shè)離散后的每個(gè)節(jié)點(diǎn)滿足平面域解析解：

其中kx＝kcosθ，kz＝ksinθ，k為波數(shù)，A為振幅，θ代表入射角，ω為圓頻率.

將公式（9）代入到求解方程（1）中，得到：

其中M為全局質(zhì)量矩陣，K為剛度矩陣，p（n）對(duì)應(yīng)著位置在（m，n）節(jié)點(diǎn)所對(duì)應(yīng)的算子類型.m的取值范圍為1≤m≤Nc，Nc對(duì)應(yīng)著不同類型算子的數(shù)目，對(duì)于2階TSEM，Nc＝12；對(duì)于3階TFEM，Nc＝18.

將公式（10）擴(kuò)展成一個(gè)矩陣方程：

其中Aα是Aj組成的一個(gè)向量，，p對(duì)應(yīng)著第j個(gè)計(jì)算點(diǎn)相鄰的離散節(jié)點(diǎn).假設(shè)求解公式（11）計(jì)算出的角頻率可以表達(dá)成：

其中C0為介質(zhì)速度，h為三角網(wǎng)格的平均長(zhǎng)度，Λi（k，θ）可以通過公式（12）計(jì)算出.在求取出角頻率ω之后，可以計(jì)算出相應(yīng)的相速度當(dāng)算子足夠精確的時(shí)候，我們認(rèn)為計(jì)算出的相速度Cp應(yīng)該接近于介質(zhì)速度C0，因此Cp／C0可以用于判斷數(shù)值頻散的大小.

對(duì)于有限元類的方法，求解公式（12）將得到Nc個(gè)不同的角頻率，每個(gè)計(jì)算出的角頻率對(duì)應(yīng)著一個(gè)頻散因子，但只有一個(gè)頻散因子具有物理意義.根據(jù)頻散曲線的收斂性，我們可以確定相應(yīng)的頻散特性.圖4a為2階TSEM的頻散曲線，橫軸對(duì)應(yīng)著空間步長(zhǎng)和波長(zhǎng)的比值，縱軸為頻散誤差的大小，不同的曲線代表著不同的入射角度θ.分析結(jié)果表明，當(dāng)入射角度為30°時(shí)，頻散誤差最小.圖4b為3階TFEM的頻散曲線，此處有限元采用的是集中質(zhì)量矩陣進(jìn)行處理，分析表明，當(dāng)入射角為45°時(shí)，頻散最小.當(dāng)誤差閥值設(shè)為0.1%時(shí)，頻散曲線表明2階TSEM需要空間采樣率滿足每個(gè)波長(zhǎng)內(nèi)平均有10個(gè)計(jì)算節(jié)點(diǎn)；而TFEM需要每個(gè)波長(zhǎng)配置16個(gè)節(jié)點(diǎn).結(jié)果顯示，對(duì)于‘X’類型的網(wǎng)格，2階TSEM相比于3階TFEM具有更高的計(jì)算精度.

4 穩(wěn)定性分析

在得到頻散因子之后可以很容易分析其穩(wěn)定性條件.假設(shè)求解波動(dòng)方程時(shí)，時(shí)間因子通過2階差分進(jìn)行求解，則公式（11）變?yōu)椋?/p>

其中τ為時(shí)間步長(zhǎng).同樣地可以求解出角頻率ω：

根據(jù)公式（14），可以得到以下關(guān)系式：

5 實(shí)例分析

5.1 數(shù)值算例

為了比較2階TSEM和3階TFEM的數(shù)值精度，考慮求解以下波動(dòng)問題：

圖4 TSEM和TFEM的頻散曲線（a）2階TSEM頻散曲線，橫軸為空間采樣率，縱軸為頻散誤差.可以看出，采樣率越低，頻散誤差越小.當(dāng)誤差閥值取為0.1%時(shí)，對(duì)應(yīng)著空間采樣率應(yīng)當(dāng)大于每個(gè)波長(zhǎng)10個(gè)采樣點(diǎn)；（b）3階TFEM頻散曲線，當(dāng)頻散誤差閥值為0.1%時(shí)要求空間采樣率大于每個(gè)波長(zhǎng)內(nèi)配置16個(gè)計(jì)算節(jié)點(diǎn).Fig.4 Dispersion curves for TSEM and TFEM（a）The dispersive result for the 2nd-order Cohen nodes.It suggests that 10nodes per wavelength are sufficient to constrain the dispersive error into 0.1%；（b）The dispersive curve for 3rd-order FEM in triangular mesh，where the curves indicate that 16nodes per wavelength is sufficient for obtaining an accurate result.

圖5 穩(wěn)定性因子分析（a）2階TSEM的穩(wěn)定性條件；（b）3階TFEM的穩(wěn)定性條件.Fig.5 The computed stability condition（a）The stability condition curve for the 2nd-order Cohen nodes，the minimum stability factor is 0.43；（b）The stability condition curve for 3rd-order FEM with triangular mesh，the minimum stability factor is 0.35.

計(jì)算結(jié)果顯示，波場(chǎng)的主頻率大概在0.5Hz，如果假設(shè)最高頻率為主頻的1.5倍，那么相當(dāng)于最短波長(zhǎng)約為1.5m，即空間采樣率為一個(gè)波長(zhǎng)15個(gè)格點(diǎn)，屬于高采樣率范疇.從計(jì)算的誤差絕對(duì)值對(duì)比上看，3階TFEM的誤差絕對(duì)值要大于2階TSEM方法，即2階TSEM法具有更高的計(jì)算精度，這也與我們之前的頻散分析結(jié)果一致.

圖6 2階TSEM（虛線）和3階TFEM（實(shí)線）計(jì)算結(jié)果的誤差絕對(duì)值對(duì)比（a）記錄點(diǎn)放置于（1.0，1.0），其中2階TSEM的誤差絕對(duì)值要小于3階TFEM；（b）記錄點(diǎn)放置于（1.5，1.5），同樣地，2階TSEM的誤差絕對(duì)值要小于3階TFEM.對(duì)比顯示，2階TSEM精度要優(yōu)于3階TFEM.Fig.6 The comparison of computational errors for TSEM and TFEM，where two mesh systems are employed to assure the same sampling ratio（a）The computed results recorded at the position of（1.0，1.0），which shows that the 2nd-order TSEM has less computational errors than the 3rd-order TFEM；（b）The comparison results recorded at the position of（1.5，1.5），the computational error also indicates that the 2nd-order TSEM scheme generates more accurate result than the 3rd-order TFEM.

5.2 孔洞模型I

我們引入一個(gè)中國(guó)西部地區(qū)比較典型的孔洞模型，如圖7所示.這類孔洞介質(zhì)經(jīng)常在地震剖面上顯示出一種串珠狀的信號(hào)（圖8）.模型設(shè)計(jì)的深度為1km，寬2km，并考慮幾個(gè)不同大小形狀各異的孔洞（尺寸大小約在20m左右）：其中第一個(gè)孔洞的橫向?qū)挾茸畲螅谌齻€(gè)寬度最小，三種孔洞的速度都設(shè)置為1500m·s-1.我們通過使用三角網(wǎng)格對(duì)孔洞進(jìn)行刻畫，使用2階TSEM方法進(jìn)行地震正演計(jì)算，并分析不同大小孔洞對(duì)應(yīng)于地震波場(chǎng)的響應(yīng)特征以及串珠狀信號(hào)的形成機(jī)理.

圖7 孔洞模型I示意圖，空洞速度填充為1500m·s-1（a）孔洞模型；（b）孔洞的具體形狀.Fig.7 The geological model with three caves of different shapes，where the velocity in the cave is set to be 1500m·s-1（a）The geological model；（b）The details for the caves.

圖9 2階TSEM正演模擬計(jì)算快照結(jié)果（a）360ms波場(chǎng)快照；（b）420ms波場(chǎng)快照；（c）480ms波場(chǎng)快照；（d）600ms波場(chǎng)快照.Fig.9 The computed results generated by the second-order TSEM with Cohen nodes（a）The snapshot for 360ms；（b）The snapshot for 420ms；（c）The snapshot for 480ms；（d）The snapshot for 600ms.

在正演計(jì)算中，震源采用主頻30Hz的Ricker子波，置于模型地表中間.不同時(shí)刻的波場(chǎng)快照結(jié)果如圖9所示：波場(chǎng)傳播到第二個(gè)孔洞的時(shí)候（420ms）產(chǎn)生了明顯的繞射波場(chǎng).這可以用惠更斯原理進(jìn)行解釋：波場(chǎng)傳播到一個(gè)介質(zhì)間斷面之后，相當(dāng)于在間斷面產(chǎn)生了一個(gè)新的波源進(jìn)行傳播.波場(chǎng)在480ms的時(shí)候，其他兩個(gè)孔洞也發(fā)生了繞射，從繞射能量而言，第一個(gè)孔洞產(chǎn)生的繞射波場(chǎng)最強(qiáng).到600ms的時(shí)候，第一個(gè)孔洞又產(chǎn)生了一個(gè)繞射波，這說明在孔洞介質(zhì)當(dāng)中，一次繞射波產(chǎn)生之后，波場(chǎng)又傳播到孔洞進(jìn)行了二次繞射，我們解釋這種現(xiàn)象為孔洞中產(chǎn)生的層間多次波.因?yàn)樵诳锥粗兴俣仍O(shè)置為水的速度，相對(duì)于圍巖是一個(gè)低速帶，這種情況很容易產(chǎn)生多次波.此外，在波場(chǎng)快照中其他幾個(gè)孔洞也相應(yīng)地產(chǎn)生了多次繞射波，但是能量相比于第一個(gè)孔洞，要小得多.地震正演模擬的結(jié)果展現(xiàn)出繞射波的能量與孔洞寬度相關(guān).寬度最大的第一個(gè)孔洞產(chǎn)生的繞射能量最強(qiáng)，而寬度最小的第三個(gè)孔洞產(chǎn)生的繞射能量最弱.這里可以用橫向分辨率去進(jìn)行解釋，在間斷面小于橫向分辨極限的情況下，所有在間斷面上的繞射波場(chǎng)傳播到檢波器之后被當(dāng)成一種信號(hào)接收.因此，當(dāng)孔洞的寬度越大，同時(shí)被接收到的信號(hào)越多，相應(yīng)的能量則越強(qiáng).

5.3 孔洞模型II

我們研究不同的孔洞組合對(duì)于地震波場(chǎng)的響應(yīng)特征.模型的大小為2km×1km，模型中考慮兩種孔洞的組合模式，分別為平行式孔洞組合（圖10b）以及三角式孔洞組合（圖10c）.同樣地使用三角網(wǎng)格對(duì)孔洞進(jìn)行刻畫，網(wǎng)格大小大約為10m，采用2階TSEM方法進(jìn)行正演模擬，地震主頻取為40Hz，位置放于模型地表中央.

正演計(jì)算之后得到的波場(chǎng)快照如圖11所示.250ms的時(shí)候，如圖11a，來自于地下起伏界面的反射波可以清晰地觀察到；當(dāng)波場(chǎng)傳播到270ms的時(shí)候，第一組孔洞組分別所產(chǎn)生的繞射波已經(jīng)出現(xiàn)在快照中；30ms之后，第一組孔洞組分別產(chǎn)生的多次波可以被觀察到，而此時(shí)第二組孔洞的繞射波已經(jīng)形成，但是各溶洞產(chǎn)生的繞射波場(chǎng)疊合在一起，比較難區(qū)分；在320ms的時(shí)候，第一組孔洞組中各孔洞產(chǎn)生的多次波場(chǎng)和繞射波場(chǎng)各自繼續(xù)傳播，并且具有較強(qiáng)的振幅；而第二組孔洞組合在波場(chǎng)上基本上體現(xiàn)為一個(gè)繞射波場(chǎng)和多次波場(chǎng)的傳播形式，且能量較第一組孔洞組更弱一些.以上快照結(jié)果顯示，平行式排布的孔洞組合在地震上具有更強(qiáng)的能量，而且各孔洞產(chǎn)生的多次波場(chǎng)比較明顯，在偏移結(jié)果上應(yīng)該會(huì)出現(xiàn)多組串珠現(xiàn)象.而三角式排布的孔洞組合在地震上的繞射能量比較弱，且各個(gè)孔洞產(chǎn)生的多次波場(chǎng)難以辨別，分析其原因是因?yàn)槿桥挪挤绞饺菀资沟貌▓?chǎng)干涉，因此預(yù)測(cè)該類孔洞組合在偏移剖面上的串珠現(xiàn)象會(huì)比較模糊，加上橫向分辨率的因素所產(chǎn)生的繞射波場(chǎng)不容易歸位形成明顯的串珠組現(xiàn)象.

6 結(jié)論

圖8 中國(guó)西部地區(qū)典型的串珠狀現(xiàn)象Fig.8 One seismic section obtained from north-western China.Some typical phenomenal of“serials reflections”can be observed in the section after employing time migration，as signed in the graph

圖10 孔洞組合模型，其中孔洞大小約為50m（a）孔洞模型II，考慮兩種不同的孔洞組合模型；（b）平行排列孔洞組合模式；（c）三角排列孔洞組合模式.Fig.10 The geological model with different types of caves（a）The velocity model，where two types of caves group are considered；（b）The triangular meshes characterizing the first group caves，where three caves are listed in a line；（c）The triangular mesh discretized for the second group of caves，where the caves are placed in a shape of triangle.

譜元法通過將單元內(nèi)部節(jié)點(diǎn)布置于特殊的積分點(diǎn)之上，并結(jié)合Lagrange形函數(shù)可以得到對(duì)角化的質(zhì)量矩陣，從而在保證精度和網(wǎng)格靈活性的同時(shí)增加了其計(jì)算效率.常規(guī)的譜元法都是應(yīng)用在四邊形網(wǎng)格之上，而本文介紹了在三角網(wǎng)格中譜元法計(jì)算框架的構(gòu)建，通過采用Cohen積分點(diǎn)去替代四邊形網(wǎng)格中的GLL積分點(diǎn)，同樣地可以保留質(zhì)量矩陣對(duì)角化的優(yōu)勢(shì).頻散分析和穩(wěn)定性分析表明，采用2階Cohen積分點(diǎn)的三角譜元方法，相對(duì)于使用集中質(zhì)量矩陣的3階三角有限元方法具有更高的計(jì)算精度以及更寬松的穩(wěn)定性條件.在本文的最后引入了兩個(gè)中國(guó)西部地區(qū)典型的孔洞地質(zhì)模型，并使用2階三角譜元法進(jìn)行孔洞的刻畫和正演模擬，其計(jì)算結(jié)果證明了該方法在模擬復(fù)雜構(gòu)造情況下的有效性.

圖11 孔洞組合模型波場(chǎng)快照結(jié)果（a）250ms波場(chǎng)快照；（b）270ms波場(chǎng)快照；（c）300ms波場(chǎng)快照；（d）320ms波場(chǎng)快照.Fig.11 The snapshot results for model with caves combination（a）The snapshot for 250ms；（b）The snapshot for 270ms；（c）The snapshot for 300ms；（d）The snapshot for 320ms.

Alford R M，Kelly K R，Boore D M.1974.Accuracy of finitedifference modeling of the acoustic wave equation.Geophysics，39（6）：834-842.

Blyth M G，Pozrikidis C.2006.A Lobatto interpolation grid over the triangle.IMA J.Appl.Math，71（1）：153-169.

Carcione J M，Herman G C，Ten Kroode A P E.2002.Seismic modeling.Geophysics，67（4）：1304-1325.

Chen Q，Babuska I.1995.Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle.Comput.Methods Appl.Mech.Eng.，128（3-4）：405-417.

Chin-joe-kong M J S，Mulder W A，Van Veldhuizen M.1999.Higher-order triangular and tetrahedral finite elements with mass lumping for solving the wave equation.J.Eng.Math.，35（4）：405-426.

Cohen G，Joly P，Tordjman N.1995.Higher order triangular finite elements with mass lumping for the wave equation.／／Third international conference on Mathematical and Numerical Aspects of wave propagation.Philadelphia，PA：SIAM，270-279.

Cohen G，Joly P，Roberts J E，etal.2001.Higher order triangular finite elements with mass lumping for the wave equation.SIAM J.Numer.Anal.，38（6）：2047-2078.

Cohen G.2002.Higher-Order Numerical Methods for Transient Wave Equations：Scientific Computation.New York：Springer-Verlag.

De Basabe J D.2009.High-order finite element methods for seismic wave propagation［Ph.D.thesis］.Texas，University of Texas at Austin.

De Basabe J D，Sen M K.2007.Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations.Geophysics，72（6）：81-95.

Deng Y Q，Zhang Z L.1990.Three-dimensional finite element modelling of elastic waves.Chinese J.Geophys.（in Chinese），33（1）：41-53.

Gerald R.1989.Applied Numerical Analysis.Addison-Wesley Publishing Co.

Hesthaven J S.1998.From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex.SIAM J.Numer.Anal.，35（2）：655-676.

Kummer B，Behle A，Dorau F.1987.Hybrid modeling of elastic wave propagation in two dimensional laterally inhomogeneous media.Geophysics，52（6）：765-771.

Lisitsa V，Podgornova O，Tcheverda V.2010.On the interface error analysis for finite difference wave simulation.Computational Geosciences，14（4）：769-776.

Komatitsch D，Vilotte J P.1998.The spectral-element method：an efficient tool to simulate the seismic response of 2Dand 3D geological structures.Bull.Seism.Soc.Am.，88（2）：368-392.

Komatitsch D，Tromp J.1999.Introduction to the spectral element method for three dimensional seismic wave propagation.Geophys.J.Int.，139（3）：806-822.

Komatitsch D，Martin R，Tromp J，etal.2001.Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles.J.Compu.Acoust.，9（2）：703-718.

Komatitsch D，Tromp J.2003.A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation.Geophys.J.Int.，154（1）：146-153.

Liu H W，Li B，Liu H，etal.2010.The algorithm of high order finite difference pre-stack reverse time migration and GPU implementation.ChineseJ.Geophys.（in Chinese），53（7）：1725-1733.

Liu T，Sen K M，Hu T Y，etal.2012.Dispersion analysis of the spectral element method using a triangular mesh.Wave Motion，49（4）：474-483.

Marfurt K J.1984.Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations.Geophysics，49（5）：533-549.

Mercerat E D，Vilotte J P，Sánchez-Sesma F J.2006.Triangular spectral element simulation of two dimensional elastic wave propagation using unstructured triangular grids.Geophys.J.Int.，166（2）：679-698.

Min D，Shin C，Pratt R G，etal.2003.Weighted-averaging Finiteelement method for 2DElastic wave equations in the frequency domain.Bull.Seism.Soc.Am.，93（2）：904-921.

Moczo P，Kristek J，Halada L.2000.3Dfourth-order staggered grid finite difference schemes：stability and grid dispersion.Bull.Seism.Soc.Am.，90（3）：587-603.

Mulder W A.1996.A comparison between higher-order finite elements and finite differences for solving the wave equation.Proceedings of the Second ECCOMAS Conference，1-7.

Mulder W A.1999.Spurious modes in finite-element discretizations of the wave equation may not be all that bad.Appl.Numer.Math.，30（4）：425-445.

Mullen R，Belytschko T.1982.Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation.Int.J.Numer.Methods Eng.，18（1）：11-29.

Seriani G，Oliveira S P.2008.Dispersion analysis of spectral element methods for elastic wave propagation.Wave Motion，45（6）：729-744.

Taylor M A，Wingate B A，Vincent R E.2000.An algorithm for computing Fekete points in the triangle.SIAM J.Numer.Anal.，38（5）：1708-1720.

Xue D C，Wang S X，Jiao S J.2007.Wave equation finite-element modeling including rugged topography and complicated medium.Progress in Geophysics（in Chinese），22（2）：522-529.

Zahradnik J，Moczo P，Horn F.1993.Testing four elastic finite difference schemes for behavior at discontinuities.Bull.Seism.Soc.Am.，83（1）：107-129.

Zhang J，Verschuur D J.2002.Elastic wave propagation in heterogeneous anisotropic media using the lumped finite-element method.Geophysics，67（2）：625-638.

Zhang J F.1994.An accurate method for finite-element elastic wave field modeling.Oil Geophysical Prospecting（in Chinese），29（1）：29-36.

Zhang M G，Wang M Y，Li X F，etal.2002.Finite element forward modeling of anisotropic elastic waves.Progress in Geophysics（in Chinese），17（3）：384-389.

附中文參考文獻(xiàn)

鄧玉瓊，張之力.1990.彈性波的三維有限元模擬.地球物理學(xué)報(bào)，33（1）：41-53.

劉紅偉，李博，劉洪等.2010.地震疊前逆時(shí)偏移高階有限差分算法及GPU實(shí)現(xiàn).地球物理學(xué)報(bào)，53（7）：1725-1733.

薛東川，王尚旭，焦淑靜.2007.起伏地表復(fù)雜介質(zhì)波動(dòng)方程有限元數(shù)值模擬方法.地球物理學(xué)進(jìn)展，22（2）：522-529.

張劍鋒.1994.一種高精度有限元彈性波場(chǎng)正演方法.石油地球物理勘探，29（1）：29-36.

張美根，王妙月，李小凡等.2002.各向異性彈性波場(chǎng)的有限元數(shù)值模擬.地球物理學(xué)進(jìn)展，17（3）：384-389.