三角譜元法及其在地震正演模擬中的應(yīng)用

李琳，劉韜，胡天躍＊

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

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

1 引言

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

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

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

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

2 三角譜元法的基本原理

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

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

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

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

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

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

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

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

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

圖1 不同三角單元的節(jié)點配置方式（a）2階三角有限元節(jié)點配置，1個單元內(nèi)部配置6個節(jié)點，均勻分布；（b）2階三角譜元法節(jié)點配置，單元內(nèi)部具有7個節(jié)點；（c）3階三角有限元節(jié)點配置，單元內(nèi)共有10個節(jié)點信息.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積分點的配置情況（a）3階Cohen點配置，單元內(nèi)部一共有12個節(jié)點（Cohen etal.，1995）；（b）4階Cohen點配置，單元內(nèi)一共有18個節(jié)點（Mulder，1996）；（c）5階Cohen點配置，單元內(nèi)一共有30個節(jié)點（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é)點的個數(shù).根據(jù)Lagrange形函數(shù)的定義，以及選擇冪函數(shù)作為基函數(shù)，我們可以構(gòu)建出2階Cohen點下Lagrange形函數(shù)的表達式：

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

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

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

3 頻散分析

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

圖3 ‘X’形狀TSEM和TFEM的計算節(jié)點配置情況（a）2階TSEM方法，全局計算采用12個計算算子；（b）3階TFEM方法，全局計算采用18個算子.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個節(jié)點，但是組集到全空間之后，全局有12種計算算子（如圖中所標(biāo)注）.對于3階TFEM格式，單元內(nèi)部一共有10個節(jié)點，相應(yīng)的全局算子為18種.

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

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

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

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

將公式（10）擴展成一個矩陣方程：

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

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

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

4 穩(wěn)定性分析

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

其中τ為時間步長.同樣地可以求解出角頻率ω：

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

5 實例分析

5.1 數(shù)值算例

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

圖4 TSEM和TFEM的頻散曲線（a）2階TSEM頻散曲線，橫軸為空間采樣率，縱軸為頻散誤差.可以看出，采樣率越低，頻散誤差越小.當(dāng)誤差閥值取為0.1%時，對應(yīng)著空間采樣率應(yīng)當(dāng)大于每個波長10個采樣點；（b）3階TFEM頻散曲線，當(dāng)頻散誤差閥值為0.1%時要求空間采樣率大于每個波長內(nèi)配置16個計算節(jié)點.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.

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

圖6 2階TSEM（虛線）和3階TFEM（實線）計算結(jié)果的誤差絕對值對比（a）記錄點放置于（1.0，1.0），其中2階TSEM的誤差絕對值要小于3階TFEM；（b）記錄點放置于（1.5，1.5），同樣地，2階TSEM的誤差絕對值要小于3階TFEM.對比顯示，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

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

圖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正演模擬計算快照結(jié)果（a）360ms波場快照；（b）420ms波場快照；（c）480ms波場快照；（d）600ms波場快照.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.

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

5.3 孔洞模型II

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

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

6 結(jié)論

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

圖11 孔洞組合模型波場快照結(jié)果（a）250ms波場快照；（b）270ms波場快照；（c）300ms波場快照；（d）320ms波場快照.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.

附中文參考文獻

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

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

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

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

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