起伏地表彈性波傳播的間斷Galerkin有限元數(shù)值模擬方法

薛昭，董良國(guó)，李曉波，劉玉柱

1同濟(jì)大學(xué)海洋地質(zhì)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092

2中國(guó)石化地球物理重點(diǎn)實(shí)驗(yàn)室，南京 210000

1 引言

地表起伏是當(dāng)今山前帶地震勘探面臨的一個(gè)難題，它除了引起觀測(cè)系統(tǒng)設(shè)計(jì)、靜校正、地震成像等方面的問(wèn)題外，還會(huì)導(dǎo)致地震資料的低信噪比問(wèn)題，成為制約目前山前帶地震成像質(zhì)量的主要瓶頸之一.通過(guò)起伏地表?xiàng)l件下地震波傳播數(shù)值模擬，可以為山前帶地震勘探的數(shù)據(jù)采集、資料處理和解釋提供理論指導(dǎo)，這是目前解決山地地震勘探中的復(fù)雜觀測(cè)系統(tǒng)設(shè)計(jì)、低信噪比、靜校正以及地震成像等問(wèn)題的有效途徑.

目前，地震波傳播數(shù)值模擬有多種方法，其中，高階有限差分法因其高效便捷的優(yōu)點(diǎn)應(yīng)用最為廣泛（董良國(guó)等，2000a，2000b；Dablain，1986；Levander，1988）.但這些模擬方法多數(shù)所采用的規(guī)則網(wǎng)格在地下介質(zhì)橫向變速劇烈或地形起伏較大時(shí)，必然會(huì)導(dǎo)致模型的階梯狀離散，從而引起數(shù)值散射等問(wèn)題.另外，在起伏地表?xiàng)l件下，這些方法（裴正林，2004；董良國(guó)，2005；Hestholm and Ruud，1994；Tessmer etal.，1992）的自由邊界條件的實(shí)現(xiàn)精度較低，地形起伏劇烈時(shí)，一些方法也極易產(chǎn)生數(shù)值計(jì)算的不穩(wěn)定.而傳統(tǒng)的有限元法和有限體法能夠采用非結(jié)構(gòu)網(wǎng)格對(duì)復(fù)雜模型進(jìn)行貼體剖分，特別是采用三角形（或四面體）網(wǎng)格能夠?qū)θ我鈴?fù)雜的起伏地表模型進(jìn)行有效剖分.有限元法很早就被應(yīng)用于地震波傳播的數(shù)值模擬（Marfurt，1984；Serón etal.，1990；Padovani etal.，1994；邵秀民等，1995；周輝等，1997；Zhang and Verschuur，2002；楊頂輝，2002；黃自萍等，2004），但在地震勘探工業(yè)界并未得到廣泛應(yīng)用，主要是由于傳統(tǒng)有限元法存在以下幾方面的問(wèn)題：（1）計(jì)算復(fù)雜，且需要對(duì)大型質(zhì)量矩陣求逆，極大增加了計(jì)算量；（2）一般使用低階空間離散格式，存在較強(qiáng)的數(shù)值頻散；（3）使用高階空間格式離散時(shí)容易產(chǎn)生偽波等.而傳統(tǒng)的有限體法也因?yàn)榫炔桓?，在地震波傳播模擬中應(yīng)用不多.譜元法（SEM）（Komatitsch and Vilotte，1998）是一種高階有限元法，但SEM一般使用四邊形和六面體網(wǎng)格，而目前生成適應(yīng)復(fù)雜介質(zhì)分布的四邊形和六面體網(wǎng)格的技術(shù)還不是很成熟，難以適應(yīng)地表起伏劇烈、地下構(gòu)造復(fù)雜情況下地震波數(shù)值模擬的需要.

間斷Galerkin有限元法（DG-FEM）是近年來(lái)發(fā)展較快的一種改進(jìn)有限元法，它最早應(yīng)用于求解原子傳輸方程（Reed and Hill，1973）.Cockburn和Shu等（1989，2001）（Cockburn etal.，2004）將DGFEM與Runge-Kutta時(shí)間離散格式結(jié)合，發(fā)展了RKDG方法，并對(duì)該方法的發(fā)展做出了很大貢獻(xiàn)，之后，該方法被廣泛地運(yùn)用于電磁學(xué)（Cockburn etal.，2004）、氣動(dòng)聲學(xué)（Toulopoulos and Ekaterinaris，2004）等學(xué)科中.

基于數(shù)值流通量理論的DG-FEM本質(zhì)上是有限元法和有限體法的結(jié)合，在單元內(nèi)部使用有限元法處理，在單元邊界上采用了有限體法中數(shù)值流通量的處理思想.DG-FEM繼承了有限元法和有限體法的優(yōu)點(diǎn)，同時(shí)也克服了這兩個(gè)方法的一些缺陷，能夠?qū)崿F(xiàn)高精度、低頻散、有效的數(shù)值模擬.它可以使用非結(jié)構(gòu)網(wǎng)格單元（包括三角形或四面體網(wǎng)格），能夠根據(jù)介質(zhì)的分布特征設(shè)計(jì)出最優(yōu)的網(wǎng)格，因此，能夠適應(yīng)起伏地表及復(fù)雜構(gòu)造條件下地震波傳播數(shù)值模擬.它繼承了有限體法良好的局部特性，可以逐單元地求解彈性波方程，避免傳統(tǒng)有限元法的大型矩陣求逆過(guò)程，且其局部特性也非常有利于算法的并行化.同時(shí)，該特性也允許DG-FEM在單元內(nèi)部使用高階有限元，從而擁有高階精度.Kaser和Dumbser（2006a）將任意高階導(dǎo)數(shù)（ADEG）時(shí)間離散格式結(jié)合DG-FEM運(yùn)用于求解彈性波動(dòng)力學(xué)方程，ADEG格式采用空間導(dǎo)數(shù)替換時(shí)間導(dǎo)數(shù)的思想，因此該方法在時(shí)間和空間上均可以達(dá)到任意高階精度.之后，他們又將該方法拓展到三維（Dumbser and K?ser，2006）、黏彈（Kaser and Dumbser，2006）、各向異性（De la Puente etal.，2007）等介質(zhì)中彈性波傳播數(shù)值模擬中，并發(fā)展了局部時(shí)間步長(zhǎng)等（Dumbser etal.，2007；Hermann etal.，2011）求解策略.Etienne等（2010）基于中心數(shù)值流，采用低階（≤2）的DGFEM求解彈性波方程.因?yàn)樵摲椒ɑ谥行臄?shù)值流，速度場(chǎng)的更新只與應(yīng)力場(chǎng)有關(guān)，而應(yīng)力場(chǎng)的更新只與速度場(chǎng)有關(guān)，所以非常適用于蛙跳（leap-frog）時(shí)間離散格式.

在本文的第2部分，介紹了利用任意高階間斷Galerkin有限元法模擬彈性波傳播的基本原理，并在附錄部分基于配點(diǎn)法對(duì)其完成細(xì)節(jié)作了詳細(xì)推導(dǎo).在時(shí)間離散方面，我們提出首先將空間離散后形成的非齊次線性常微分方程系統(tǒng)齊次化，然后再采用針對(duì)齊次問(wèn)題的強(qiáng)穩(wěn)定性保持龍格庫(kù)塔（SSP Runge-kutta）時(shí)間格式，最終獲得的算法在時(shí)間和空間上均可以達(dá)到任意高階.在第3部分，基于近最佳匹配層（NPML）的思想和復(fù)頻移（CFS）拉伸坐標(biāo)變換，我們推導(dǎo)了一種新的PML吸收邊界條件——CFS-NPML，該技術(shù)能夠有效地應(yīng)用到DG-FEM中.通過(guò)第4部分的數(shù)值試驗(yàn)來(lái)說(shuō)明發(fā)展的模擬方法的有效性.

另外，DG-FEM方法需要對(duì)介質(zhì)模型進(jìn)行貼體剖分，為此我們開(kāi)發(fā)了一種模型的自適應(yīng)三角剖分方法，該方法可以根據(jù)介質(zhì)參數(shù)的變化程度進(jìn)行自適應(yīng)貼體剖分.這種自適應(yīng)三角剖分方法我們將另文發(fā)表，這里不多贅述.

2 方法原理

二維速度-應(yīng)力彈性波方程可以表示為下面的一階變系數(shù)雙曲型偏微分方程系統(tǒng)：

2.1 彈性波方程的DG-FEM公式

設(shè)計(jì)算區(qū)域?yàn)棣福瑢?duì)其進(jìn)行有限單元離散，劃分為不重疊的子區(qū)域.間斷Galerkin有限元法采用分片（逐單元）光滑函數(shù)逼近波場(chǎng)，即不要求波場(chǎng)值在單元邊界上保持連續(xù).對(duì)方程（1）兩邊同時(shí)乘以標(biāo)量試驗(yàn)函數(shù)，并在第k個(gè)單元上進(jìn)行空間積分，有

其中，n為單元邊界?Ωk的外法向矢量.

由于在單元邊界上不要求波場(chǎng)連續(xù)，因此需要對(duì)單元邊界上波場(chǎng)值給予定義，這里暫用uk＊表示.對(duì)上式右端第一項(xiàng)可以再次使用高斯散度定理推論，得

假設(shè)單元內(nèi)物性參數(shù)為常數(shù)，則最終可以將上式整理并簡(jiǎn)寫(xiě)為

間斷Galerkin有限元法通過(guò)使用有限體法中的數(shù)值流通量來(lái)定義單元邊界上的波場(chǎng)值，這里，我們采用局部Lax-Friedrichs數(shù)值流通量.當(dāng)采用高階DG-FEM時(shí)，空間格式的精度主要由單元內(nèi)部自由度驅(qū)動(dòng)，使用局部Lax-Friedrichs數(shù)值流通量是比較合適的選擇.下面給出局部Lax-Friedrichs數(shù)值流通量的定義：

其中，cp為本單元的縱波速度.當(dāng)設(shè)（6）式中cp為零時(shí)，即為中心數(shù)值流通量.uk-，uk+分別為單元邊界處法線內(nèi)外側(cè)的波場(chǎng)（可以理解為左右極限），uk-由本單元控制，uk+由相鄰單元控制；Q定義為數(shù)值流通量系數(shù)矩陣.

圖1 物理單元到參考單元映射示意圖Fig.1 Mapping between physical and canonical element

為了計(jì)算方便，通過(guò)坐標(biāo)變換將物理單元Ωk映射到標(biāo)準(zhǔn)單元Ωc上（圖1）.當(dāng)采用直邊多邊形網(wǎng)格單元時(shí)，物理單元Ωk與參考單元Ωc之間以及單元對(duì)應(yīng)邊之間的Jacobi映射關(guān)系為常數(shù)，則（5）式可進(jìn)一步簡(jiǎn)化為

其中，Jk為單元體之間的Jacobi系數(shù)，為單元第e個(gè)直邊之間的Jacobi系數(shù).

2.2 空間離散方程

上面推導(dǎo)中一階雙曲型系統(tǒng)方程是一般形式，聲波方程、各向同性和各向異性彈性波方程、黏彈方程等均可以表達(dá)為一階雙曲方程系統(tǒng)，因此，它們都可以通過(guò)（5）式的DG-FEM求解.雖然DG-FEM并不限定網(wǎng)格單元類(lèi)型，但這里我們主要考慮直邊三角形網(wǎng)格單元.

按照有限元的思路，在每個(gè)單元內(nèi)，用多項(xiàng)式基函數(shù)逼近波場(chǎng)，當(dāng)采用的基函數(shù)與試驗(yàn)函數(shù)相同時(shí)，即為Galerkin方法.采用在三角網(wǎng)格單元上具有正交特性的Koornwinder-Dubiner多項(xiàng)式（附錄C），波場(chǎng)分量ukm（x，t）展開(kāi)表示為

假設(shè)單元Ωk內(nèi)介質(zhì)物性參數(shù)不變，即A、B為常系數(shù)矩陣.將（8）式代入（7）式，并選擇合適的數(shù)值流，經(jīng)過(guò)數(shù)學(xué)推導(dǎo)，最終在單元Ωk上的半離散方程為

獲得空間半離散方程后，再加上相應(yīng)邊界條件以及合適的時(shí)間離散格式，我們便可以對(duì)波場(chǎng)進(jìn)行時(shí)間方向上的更新.從（9）式中可以看到，第k個(gè)單元上波場(chǎng)的更新只與本單元以及相鄰三個(gè)單元有關(guān)，即該空間格式是緊致顯式的，我們可以逐單元更新波場(chǎng)，而不需要形成整體的有限元方程.

2.3 時(shí)間積分格式

作為概念上的表示，我們也可以將（9）式的空間離散方程整合成下面形式的線性方程系統(tǒng)：

其中，u（t）表示所有待更新的系數(shù)變量（全局自由度），設(shè)維度為n×1，L為n×n的全局空間離散算子，f為震源積分向量.

（10）式的解析解可以表示為

其中，t0為波場(chǎng)更新的初始時(shí)刻，u0為初始波場(chǎng)值.

Al-Mohy和Higham（2011）通過(guò)對(duì)f（t）進(jìn)行p階Taylor級(jí)數(shù)展開(kāi)，最終將（11）式的求解轉(zhuǎn)化為一個(gè)擴(kuò)展矩陣的指數(shù)函數(shù)的求解：

其中，In為n×n的單位矩，L為 （n+p）×（n+p）的擴(kuò)展矩陣，u0為 （n+p）×1的列向量，W＝為n×p的矩陣，I為（p-1）×（p-1）的單位矩陣，n，η有利于數(shù)值計(jì)算穩(wěn)定，ep為第p個(gè)元素不為零的單位向量.

Runge-Kutta方法是最早與DG-FEM結(jié)合使用的一種時(shí)間格式.目前，該方法從存儲(chǔ)以及精度上都已獲得了相當(dāng)大發(fā)展.其中，強(qiáng)穩(wěn)定保持Rugge-Kutta格式是一類(lèi)具有較強(qiáng)穩(wěn)定性的RK格式.盡管方法設(shè)計(jì)的出發(fā)點(diǎn)主要是針對(duì)非線性問(wèn)題，但事實(shí)證明，它對(duì)于常系數(shù)的線性問(wèn)題同樣有效，且能達(dá)到高階精度.對(duì)（13）式的齊次問(wèn)題，Gottlieb（2005）給出了下面的m級(jí)m-1階強(qiáng)穩(wěn)定性保持Rugge-Kutta時(shí)間格式（SSP-RK（m，m-1）），其表達(dá)式為

其中，αm，k為常系數(shù)，滿(mǎn)足下列關(guān)系式：

表1 SSP-RK（m，m-1）系數(shù)表Table 1 Coefficientsαm，jof linear SSP-RK（m，m-1）

需要指出的是矩陣-向量乘的過(guò)程就是實(shí)施DG-FEM有限元空間離散的過(guò)程，因此并不需要顯式的給出大型矩陣~L.

3 震源及邊界條件

3.1 震源

震源有不同的形式，一般需要模擬方向震源或爆炸震源的激發(fā).這里討論爆炸震源，設(shè)xs為震源位置，且在單元Ωks內(nèi).施加爆炸震源時(shí)，將子波作用在正應(yīng)力上比較方便，設(shè)正應(yīng)力項(xiàng)上加載震源為f（x，t）＝s（t）δ（x-xs），其中，s（t）為震動(dòng)函數(shù).〈φksi，f（x，t）〉Ωks＝s（t）〈φksi，δ（x-xs）〉Ωks＝s（t）φksi（xs）.因此，在模擬過(guò)程中，只需要在震源所在單元的正應(yīng)力對(duì)應(yīng)項(xiàng)加上向量

當(dāng)采用低階的DG法或者震源所在的網(wǎng)格單元較大時(shí)，這種施加方法產(chǎn)生的縱波源不純，殘留一些橫波的痕跡.為了壓制這種現(xiàn)象，可以在以震源為中心的高斯分布區(qū)域內(nèi)施加.在物理空間完成方式中，可以在所有節(jié)點(diǎn)（附錄B、C）上施加震源，具體可以直接將子波乘以節(jié)點(diǎn)對(duì)應(yīng)高斯分布系數(shù)，然后加到相應(yīng)正應(yīng)力項(xiàng)上.

3.2 自由表面條件

自由邊界條件要求作用在自由表面上的應(yīng)力為零.在DG-FEM法中，數(shù)值流通量控制單元邊界上的波場(chǎng)值，因此，可以通過(guò)設(shè)置（6）式中的數(shù)值流通量來(lái)實(shí)施自由邊界條件.本文采用下面的方式來(lái)簡(jiǎn)化自由邊界條件的實(shí)施：在自由表面處，采用中心數(shù)值流通量，然后采用鏡像原理，將自由表面兩側(cè)應(yīng)力分量設(shè)置為反對(duì)稱(chēng)，而速度分量設(shè)置對(duì)稱(chēng)即可.而在模型的內(nèi)部單元邊界上我們使用了更精確的Lax-Friedrichs數(shù)值流通量理論.

3.3 CFS-NPML吸收邊界條件

地震波數(shù)值模擬一般計(jì)算有限區(qū)域內(nèi)地震波的傳播，因此需要引入人工吸收邊界條件來(lái)處理人為邊界反射問(wèn)題.在DG-FEM中，單元邊界處的數(shù)值流通量項(xiàng)包括兩部分：流出單元部分和進(jìn)入單元部分，進(jìn)入單元部分由相鄰單元控制.因此，在計(jì)算區(qū)域邊界處可以通過(guò)將數(shù)值流通量的進(jìn)入單元部分設(shè)置為零來(lái)控制邊界反射.該吸收邊界條件被稱(chēng)為開(kāi)放邊界條件.但開(kāi)放邊界條件對(duì)面波、倏逝波等的吸收效果并不理想，我們需要對(duì)吸收邊界條件作進(jìn)一步改善或采用最佳匹配層（PML）吸收邊界條件.

對(duì)于PML技術(shù)，拉伸坐標(biāo)變換是一個(gè)關(guān)鍵的手段，它通過(guò)一個(gè)拉伸函數(shù)將空間坐標(biāo)系映射到復(fù)空間中去，然后對(duì)入射匹配層的平面波進(jìn)行阻尼吸收.理論上，常規(guī)PML技術(shù)可以對(duì)進(jìn)入匹配層的入射波完全吸收，但經(jīng)過(guò)空間離散后，其對(duì)近平行入射到界面上的波以及低頻波、倏逝波等吸收效果仍然欠佳.Kuzuoglu等（1996）引入了復(fù)頻移（CFS）拉伸函數(shù)，并采用了不分裂式PML.基于復(fù)頻移拉伸函數(shù)的PML（一般簡(jiǎn)稱(chēng)為CFS-PML），對(duì)于近面波、倏逝波等都有良好的吸收效果.不分裂PML不需要分裂波場(chǎng)，但在時(shí)間域需要做大量褶積運(yùn)算，計(jì)算量較大.Komatitsch和Martin（2007）發(fā)展了一種無(wú)需顯式褶積計(jì)算、基于迭代格式的不分裂PML技術(shù)——卷積型完美匹配層CPML，大大降低了不分裂PML的時(shí)間域計(jì)算代價(jià).

Cummer（2003）介紹了一種新的不分裂PML技術(shù)——NPML（Nearly Perfectly Matched Layer）.該P(yáng)ML技術(shù)在推導(dǎo)過(guò)程中的一個(gè)假設(shè)前提是拉伸函數(shù)是空間不變的.事實(shí)上，從數(shù)學(xué)上也可以證明，在該假設(shè)條件下獲得的吸收方程仍然是一種PML技術(shù).本小節(jié)基于NPML的思想，同時(shí)采用復(fù)頻移拉伸（CFS）函數(shù)，通過(guò)引入新的輔助變量推導(dǎo)了一種新的PML技術(shù)（簡(jiǎn)稱(chēng)為CFS-NPML）.

定義復(fù)頻移（CFS）拉伸函數(shù)：

其中，θ∈x，｛｝z.復(fù)拉伸函數(shù)的極點(diǎn)移動(dòng)到復(fù)平面的虛軸上對(duì)常規(guī)PML進(jìn)行了改善，可以通過(guò)調(diào)節(jié)參數(shù)更好地吸收面波以及倏逝波的人為邊界反射.我們一般將參數(shù)kθ的值取為1.參數(shù)βθ＝其中，δθ為單元中心到PML邊界的深度，為PML的厚度，為縱波速度.參數(shù)其中為子波的主頻.（16）對(duì)應(yīng)的空間導(dǎo)數(shù)變換關(guān)系為

可以將（16）式寫(xiě)為

將彈性波方程（1）變換到頻率域，并按（17）、（18）式進(jìn)行坐標(biāo)變換，有

（19）式就可以表示為

將（20）式代入（19）式，并運(yùn)用NPML的思想作變量代換，有

（21）和（22）式反變換到時(shí)間域，就得到我們的CFSNPML吸收邊界條件：

這里必須指出，表面上CFS-NPML需要10個(gè)輔助變量，但實(shí)際上和在計(jì)算過(guò)程中是不需要的，即輔助變量個(gè)數(shù)也為8個(gè)（與CPML相同）.CFS-NPML可以寫(xiě)為下面的一階雙曲系統(tǒng)形式：

其中，

對(duì)于（24）式的一階雙曲系統(tǒng)，我們同樣可以采用前面介紹的DG-FEM方法對(duì)其進(jìn)行離散求解.

4 數(shù)值試驗(yàn)

下面，我們通過(guò)幾個(gè)數(shù)值試驗(yàn)來(lái)驗(yàn)證發(fā)展的DG-FEM方法在復(fù)雜構(gòu)造以及起伏地表情況下模擬彈性波傳播的有效性.

4.1 Lamb問(wèn)題

通過(guò)Lamb問(wèn)題的求解，可以測(cè)試自由表面邊界條件實(shí)施的有效性.圖2是一個(gè)自由地表傾角為10°的均勻介質(zhì)模型，縱波速度為3200m·s-1，橫波速度為1847.5m·s-1，密度為2200kg·m-3.在自由表面（1720m，401.72m）處垂直地表施加震源（主頻為14.5Hz的雷克子波），在水平坐標(biāo)為2694.95m和3400.08m的地表上設(shè)置兩個(gè)檢波器.模擬采用空間3階，時(shí)間3階，時(shí)間步長(zhǎng)為0.2ms.

圖2 傾斜地表模型Fig.2 Geometry of the tilted surface model

圖3 為600ms時(shí)的水平和垂直速度分量波場(chǎng)快照，圖中P波、S波、Rayleigh面波以及Schmidt波等震相均清晰可見(jiàn).為了驗(yàn)證模擬方法的有效性，我們將地表記錄與解析解比較（圖4）.從圖中可以看到，無(wú)論是水平速度分量還是垂直速度分量，DGFEM模擬結(jié)果和解析解都吻合得非常好，說(shuō)明了DG-FEM方法可以獲得高精度的彈性波模擬結(jié)果.

4.2 CFS-NPML邊界吸收效果驗(yàn)證

本試驗(yàn)主要目的是驗(yàn)證本文推導(dǎo)的CFS-NPML技術(shù)對(duì)人為邊界反射的吸收效果.鑒于CPML是近年來(lái)提出的邊界吸收效果很好且得到廣泛應(yīng)用的吸收邊界條件，下面將本文發(fā)展的CFS-NPML與CPML進(jìn)行了對(duì)比分析.

試驗(yàn)?zāi)Ｐ蜑橐粋€(gè)簡(jiǎn)單的水平自由表面均勻模型，縱波速度為3200m·s-1，橫波速度為1847.5m·s-1，密度為2200kg·m-3，大小為2000m×2000m.在自由表面（50m，10m）處施加爆炸震源（主頻為25Hz的雷克子波）.采用P2階空間離散.模擬的時(shí)間步長(zhǎng)為0.25ms.PML的寬度為160m，7個(gè)單元的寬度.CPML和CFS-NPML的吸收參數(shù)相同，R＝10-7.

圖3 600ms時(shí)的波場(chǎng)快照（a）水平速度分量；（b）垂直速度分量.Fig.3 Snapshots of velocity components at 600ms for Lamb′s problem（a）Horizontal velocity component；（b）Vertical velocity component.

圖5 為采用上述兩種不同PML吸收邊界條件時(shí)1s時(shí)刻速度水平分量的波場(chǎng)快照（左圖為CFSNPML，右圖為CPML），為了顯示效果，對(duì)快照作了適當(dāng)增益.可以看到，不管是CPML，還是本文發(fā)展的CFS-PML，對(duì)面波等人為邊界反射的吸收都非常干凈.

圖4 DG-FEM（P3階）模擬結(jié)果與解析解的對(duì)比第1道：第1個(gè)檢波器接收記錄；第2道：第2個(gè)檢波器接收記錄.（a）水平分量；（b）垂直分量.Fig.4 Comparison between DG-FEM（P3）numerical results with analytical solution（a）Horizontal and（b）vertical components of the velocity.The 1st line seismogram at the first receiver，and the 2ndline at another receiver.

為了進(jìn)一步比較CFS-NPML和CPML的吸收效果，在地表（100m，0m）處設(shè)置檢波器，圖6a為水平速度分量記錄.直達(dá)波和面波的延續(xù)時(shí)間之后，記錄上未看到任何邊界反射，并且從圖6b的放大曲線可以看出，由于數(shù)值計(jì)算、邊界反射等造成的曲線震蕩峰值不到記錄面波峰值的0.267%，而且CPML和CFS-NPML對(duì)應(yīng)曲線吻合得很好，這進(jìn)一步說(shuō)明了CFS-NPML與CPML都具有非常理想的人為邊界反射吸收效果.

圖5 不同PML吸收邊界條件下，1s時(shí)刻水平速度分量波場(chǎng)快照，震源位置為（50m，10m）（a）CFS-NPML；（b）CPML.Fig.5 The snapshots of horizontal velocity component at 1scomputed by P2DG-FEM with CFS-NPML（a）and CPML（b）.The location of source is（50m，10m）

圖6 不同PML吸收邊界條件下，地表（100m，0m）處檢波器速度水平分量記錄（a）及其放大圖（b）對(duì)比Fig.6 Comparison of the horizontal velocity seismograms（a）and their local zoom（b）under CPML and CFS-NPML

表2列出了CFS-NPML與CPML計(jì)算效率對(duì)比.本實(shí)驗(yàn)的模擬總時(shí)長(zhǎng)2s，為了顯示方便，上圖只顯示了0～1s記錄.從表中可以看到，CFS-NPML較CPML具有更大的計(jì)算效率，本次試驗(yàn)節(jié)約計(jì)算時(shí)間約27%.

表2 CPML與CFS-NPML計(jì)算效率對(duì)比Table 2 The computational efficiency comparison of CPML and CFS-NPML

4.3 高陡構(gòu)造模型

利用常規(guī)的規(guī)則網(wǎng)格對(duì)復(fù)雜構(gòu)造模型進(jìn)行剖分時(shí)會(huì)形成邊界的階梯狀離散，從而造成地震波傳播的數(shù)值散射問(wèn)題.為了說(shuō)明DG-FEM方法在模擬復(fù)雜構(gòu)造模型中地震波傳播時(shí)在剖分方面的優(yōu)勢(shì)，我們?cè)O(shè)計(jì)了一個(gè)含有高陡構(gòu)造的模型（圖7），上層介質(zhì)的縱橫波速度及密度分別為3464m·s-1、2000m·s-1和2000kg·m-3，下層介質(zhì)的縱橫波速度及密度分別為4000m·s-1、2300m·s-1和2300kg·m-3.從圖8的網(wǎng)格剖分局部放大圖中可以看到，DG-FEM可以沿著高陡構(gòu)造的邊界進(jìn)行貼體的剖分，而有限差分法（FD）采用的矩形網(wǎng)格（dx＝10m）離散后高陡構(gòu)造的邊界呈現(xiàn)階梯狀.圖9對(duì)比了有限差分法（空間4階，時(shí)間2階）和DG-FEM（空間2階，時(shí)間2階）的地表速度分量記錄，從地表記錄中方框內(nèi)可以看到FD模擬結(jié)果存在明顯的數(shù)值散射.在圖10的地表1500m處的檢波器記錄上，F(xiàn)D模擬方法因?yàn)殡A梯離散造成的這種數(shù)值散射更為明顯.

圖7 高陡傾角模型Fig.7 Geometry of the high steeply-dipping model

4.4 單峰模型

本模型（圖11）為一個(gè)單峰的單層均勻介質(zhì)起伏地表模型，地表起伏呈高斯指數(shù)函數(shù)分布.模型橫向展布4500m，單峰的中心軸大約在橫向坐標(biāo)2700m處，最高峰高程為400m.縱波速度為3464m·s-1，橫波速度為2000m·s-1，密度為2000kg·m-3.在自由表面（2000m，500m）處施加爆炸震源（主頻為25Hz的雷克子波）.檢波器排列在地表1000～4000m范圍內(nèi)，151個(gè)檢波器，相鄰檢波器橫向間距20m.模擬采用空間2階，時(shí)間2階，時(shí)間步長(zhǎng)為0.25ms.

圖12為用DG-FEM方法與基于縱向坐標(biāo)變換有限差分方法分別模擬的地表速度垂直分量對(duì)比.可以發(fā)現(xiàn)，相比有限差分方法，DG-FEM方法能夠更加精細(xì)地模擬地表起伏所引起的前向和后向散射波.

圖8 網(wǎng)格剖分局部放大圖（a）DG-FEM三角網(wǎng)格剖分；（b）FD矩形網(wǎng)格剖分.Fig.8 Local zoom of the discrete description for model（a）DG-FEM triangular mesh；（b）FD regular mesh.

圖9 地表單炮記錄（上：水平分量，下：垂直分量）（a，c）有限差分結(jié)果；（b，d）DG-FEM結(jié)果.Fig.9 Surface single-shot records（above：horizontal components，below：vertical components）（a，c）FD results；（b，d）DG-FEM results.

圖10 地表1500m處檢波器記錄（a）水平分量 ；（b）垂直分量.Fig.10 Seismograms at surface location 1500m（a）Horizontal component；（b）Vertical component.

圖11 單峰模型Fig.11 Geometry of the hill model

圖13 為單峰和凹陷模型下，DG-FEM模擬的地表速度垂直分量記錄.我們可以清晰地看到，當(dāng)?shù)卣鸩ㄏ蛴覀鞑サ降匦纹鸱帟r(shí)，會(huì)出現(xiàn)波場(chǎng)的分離和轉(zhuǎn)換等現(xiàn)象，例如P波和Schmidt波的分離、P波和面波遇到界面時(shí)的轉(zhuǎn)換及再發(fā)育，這些都導(dǎo)致了散射波發(fā)育豐富.從單峰和凹陷模型對(duì)比發(fā)現(xiàn)，地形隆起時(shí)，前、后向散射波發(fā)育都比較豐富，而地形凹陷時(shí)，前向散射更加豐富.

4.5 新疆模型

圖12 單峰模型地表垂直速度分量（a）DG-FEM結(jié)果；（b）FD結(jié)果.Fig.12 Seismograms of vertical velocity component computed for the hill model（a）DG-FEM results；（b）FD results.

圖13 單峰和凹陷模型下，DG-FEM模擬的地表垂直速度分量記錄（a）單峰模型；（b）凹陷模型.Fig.13 Seismograms of vertical velocity components simulated by DG-FEM for（a）the hill model and（b）the sag model

圖14 是中石油東方地球物理公司根據(jù)準(zhǔn)噶爾盆地一個(gè)工區(qū)的地表和地下構(gòu)造特點(diǎn)設(shè)計(jì)的一個(gè)起伏地表模型（圖中顯示的是縱波速度），地下構(gòu)造比較簡(jiǎn)單，但地形局部起伏比較劇烈，并且存在低速帶（縱波速度500m·s-1）和降速帶（縱波速度800m·s-1）.在自由表面（X＝4600m）處施加爆炸震源（主頻為25Hz的雷克子波），251個(gè)檢波器排列在地表1000～6000m范圍內(nèi)，相鄰檢波器間距20m，模擬采用空間2階，時(shí)間2階，時(shí)間步長(zhǎng)為0.2ms，模擬時(shí)長(zhǎng)2s.

從圖15的模擬記錄看出，由于地表起伏以及地表覆蓋的低降速帶的存在，近地表處地震波場(chǎng)的傳播非常復(fù)雜，特別是面波散射十分嚴(yán)重，一次反射波信號(hào)被嚴(yán)重污染，大幅度降低了記錄的信噪比.因此，起伏地表所引起的地震波場(chǎng)的散射是山地地震資料信噪比低的主要原因之一.

5 結(jié)論

本文結(jié)合針對(duì)齊次問(wèn)題的強(qiáng)穩(wěn)定性保持龍格庫(kù)塔（SSP Runge-kutta）算法，將DG-FEM方法推廣至?xí)r間任意高階精度，用于求解彈性波方程.Lamb問(wèn)題等一系列數(shù)值模擬計(jì)算實(shí)例表明，DG-FEM方法是一個(gè)高階高精度的數(shù)值模擬方法.相比有限差分方法，DG-FEM方法能夠高精度地實(shí)現(xiàn)起伏地表自由邊界條件.同時(shí)，因?yàn)椴捎觅N體的網(wǎng)格剖分，它能夠避免有限差分等方法因階梯狀離散造成的數(shù)值散射等問(wèn)題，而三角形網(wǎng)格使DG-FEM方法能

圖14 新疆模型（部分）Fig.14 Geometry of the Xijiang model（partly）

圖15 地表單炮記錄（a）水平分量；（b）垂直分量.Fig.15 Surface single-shot records（a）Horizontal component；（b）Vertical component.

夠有效模擬任意復(fù)雜起伏地表?xiàng)l件下彈性波的傳播.

基于NPML的思想和復(fù)頻移拉伸坐標(biāo)變換，本文推導(dǎo)的CFS-NPML吸收邊界條件能夠與DGFEM方法有效結(jié)合，并且對(duì)面波的人為邊界反射也具有良好的吸收效果.通過(guò)數(shù)值試驗(yàn)對(duì)比CFSNPML和CPML在DG-FEM模擬彈性波傳播中的表現(xiàn)可以發(fā)現(xiàn)，二者都能有效吸收近平行邊界入射的包括面波等震相的人為邊界反射.在同樣吸收參數(shù)條件下，兩者吸收效果相當(dāng)，但本文發(fā)展的CFSNPML技術(shù)具有更高的計(jì)算效率.

本文發(fā)展的高精度的數(shù)值模擬方法可以為我們分析起伏地表?xiàng)l件下彈性波的傳播規(guī)律以及波場(chǎng)響應(yīng)特征提供一種有效的彈性波傳播數(shù)值模擬工具.

附錄A 高斯散度定理及其推論

高斯散度定理：

高斯散度定理推論：

在文中，我們作如下統(tǒng)一表示：小寫(xiě)黑體字母表示矢量；大寫(xiě)黑體字母表示矩陣；帶箭頭的小寫(xiě)黑體字母表示元素為矢量的矢量；帶箭頭的大寫(xiě)黑體字母表示元素為矩陣的矢量.

附錄B DG-FEM空間離散算子推導(dǎo)

（1）?？臻g（Modal Space）

函數(shù)f（x）定義在單元Ωk上，該式可以看作是函數(shù)在多項(xiàng)式空間的展開(kāi)，是每一個(gè)方向上的模長(zhǎng).記?？臻g：

（2）物理空間（Physics Space）

在單元Ωk上選一組離散節(jié)點(diǎn)（節(jié)點(diǎn)的選擇見(jiàn)C部分），物理空間：

（3）空間關(guān)系

?？臻g和物理空間存在映射關(guān)系，定義廣義范德蒙特矩陣：

由（B1）式及廣義范德蒙特矩陣，建立（B2）及（B3）的映射關(guān)系為

（4）函數(shù)內(nèi)積

函數(shù)g（x）、f（x）定義在單元Ωk上，則它們之間的內(nèi)積表示為

采用正交基函數(shù)，則

（5）偏導(dǎo)數(shù)內(nèi)積

函數(shù)g（x）、f（x）定義在單元Ωk上，定義下面形式的內(nèi)積：

同理，由（B4）有

將偏導(dǎo)數(shù)在?？臻g直接展開(kāi)：

同理，

則偏導(dǎo)內(nèi)積表示為

（6）函數(shù)邊界內(nèi)積

由?？臻g和物理空間存在映射關(guān)系，有

另外從區(qū)域上看，由?？臻g和物理空間存在映射關(guān)系：

最終，有：

附錄C 基函數(shù)和節(jié)點(diǎn)選擇

我們一般先將物理單元映射到標(biāo)準(zhǔn)單元上，然后在參考單元上選擇一組權(quán)函數(shù)用于上面的計(jì)算需求.在三角形單元上，我們一般采用Koornwinder-Dubiner基函數(shù)：

對(duì)于物理空間上的節(jié)點(diǎn)，我們選擇使廣義范德蒙特矩陣有更好求逆性質(zhì)的Fekete點(diǎn)集（數(shù)據(jù)來(lái)至Alvise Sommariva的個(gè)人網(wǎng)站）.

Al-M ohy A H，Higham N J.2011.Computing the Action of the Matrix Exponential，with an Application to Exponential Integrators.SIAM J.Sci.Comput.，33（2）：488-511.

Cockburn B，Li F Y，Shu C W.2004.Locally divergence-free discontinuous Galerkin methods for the Maxwell equations.Journal of Computational Physics，194（2）：588-610.

Cockburn B，Shu C W.1989.Tvb Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws ii：general framework.Mathematics of Computation，52（186）：411-435.

Cockburn B，Shu C W.2001.Runge-Kutta discontinuous Galerkin methods for convection-dominated problems.Journal of Scientific Computing，16（3）：173-261.

Cummer S A.2003.A simple，nearly perfectly matched layer for gneral electromagnetic media.IEEE Microwave and Wireless Components Letters，13（2）：128-130.

Dablain M A.1986.The application of high-order differencing to the scalar wave equation.Geophysics，51（1）：54-66.

De la Puente J，K?ser M，Dumbser M，etal.2007.An arbitrary high-order discontinuous galerkin method for elastic waves on unstructured meshes IV：anisotropy.Geophysical Journal International，169（3）：1210-1228.

Dong L G，Ma Z T，Cao J Z，etal.2000a.A staggered-grid highorder difference method of one-order elastic wave equation.Chinese J.Geophys.（in Chinese），43（3）：411-419.

Dong L G，Ma Z T，Cao J Z.2000b.A study on stability of the staggered grid high order difference method of first order elastic wave equation.Chinese J.Geophys.（in Chinese），43（6）：856-864.

Dong L G.2005.Numerical simulation of seismic wave propagation under complex near surface conditions.Progress in Exploration Geophysics（in Chinese），28（3）：187-194.

Dumbser M，K?ser M.2006.An arbitrary high order discontinuous galerkin method for elastic waves on unstructuredmeshes II：the three-dimensional isotropic case.Geophysical Journal International，167（1）：319-336.

Dumbser M，Kaser M，Toro E.2007.An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes V：local time stepping and p-adaptivity.Geophysical Journal International，171（2）：695-717.

Etienne V，Chaljub E，Virieux J，etal.2010.An hp-adaptive discontinuous Galerkin finite-element method for 3Delastic wave modelling.Geophysical Journal International，183（2）：941-962.

Gottlieb S.2005.On high order strong stability preserving Runge-Kutta and multi step time discretizations.Journal of Scientific Computing，25（1）：105-128.

Hermann V，Kaser M，Castro C E.2011.Non-conforming hybrid meshes for efficient 2Dwave propagation using the Discontinuous Galerkin method.Geophysical Journal International，184（2）：746-758.

Hestholm S，Ruud B.1994.2Dfinite-difference elastic wave modelling including surface topography.Geophysical Prospecting，42（5）：371-390.

Huang Z P，Zhang M，Wu W Q，etal.2004.A domain decomposition method for numerical simulation of the elastic wave propagation.Chinese J.Geophys.（in Chinese），47（6）：1094-1100.

K?ser M，Dumbser M，de la Puente J，etal.2007.An arbitrary high order Discontinuous Galerkin method for elastic waves on unstructured meshes III：Viscoelastic attenuation.Geophysical Journal International，168（1）：224-242.

Kaser M，Dumbser M.2006.An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes I：the two-dimensional isotropic case with external source terms.Geophysical Journal International，166（2）：855-877.

Komatitsch D，Martin R.2007.An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation.Geophysics，72（5）：155-167.

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.

Kuzuoglu M，Mittra R.1996.Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers.IEEE Microwave and Guided Wave Letters，6（12）：447-449.

Levander A R.1988.Four-order finite-difference P-SV seismograms.Geophysics，53（11）：1425-1436.

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

Padovani E，Priolo E，Seriani G.1994.Low-and high-order finite element method：experience in seismic modeling.J.Comp.Acous.，2（4）：371-422.

Reed W H，Hill T R.1973.Triangular mesh methods for the neutron transport equation，Technical Report，LA-UR-73-479，Los Alamos Scientific Laboratory.

Serón F J，Sanz F J，Kindelán M，etal.1990.Finite-element method for elastic wave propagation.Comm.Appl.Numerical Methods，6（5）：359-368.

Shao X M，Lan Z L.1995.Numerical simulation of the seismic wave propagation in inhomogeneous isotropic elastic media.Chinese J.Geophys.（in Chinese），38（S1）：39-55.

Tessmer E，Kosloff D，Behle A.1992.Elastic wave propagation simulation in the presence of surface topography.Geophys.J.Int.，108（2）：621-632.

Toulopoulos I，Ekaterinaris J A.2004.High-order discontinuous Galerkin discretizations for computational aeroacoustics in complex domains.AIAA，44（3）：502-511.

Yang D H.2002.Finite element method of the elastic wave equation and wavefield simulation in two-phase anisotropic media.Chinese J.Geophys.（in Chinese），45（4）：575-583.

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

Zhou H，Xu S Z，Liu B，etal.1997.Modeling of wave equations in anisotropic media using finite element method and its stability condition.Chinese J.Geophys.（in Chinese），40（6）：833-841.

Pei Z L.2004.Numerical modeling using staggered-grid high-order finite-difference of elastic wave equation on arbitrary relief surface.Oil Geophysical Prospecting（in Chinese），39（6）：629-634.

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

董良國(guó)，馬在田，曹景忠等.2000a.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法.地球物理學(xué)報(bào)，43（3）：411-419.

董良國(guó)，馬在田，曹景忠.2000b.一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性研究.地球物理學(xué)報(bào)，43（6）：856-864.

董良國(guó).2005.復(fù)雜地表?xiàng)l件下地震波傳播數(shù)值模擬.勘探地球物理進(jìn)展，28（3）：187-194.

黃自萍，張銘，吳文青等.2004.彈性波傳播數(shù)值模擬的區(qū)域分裂法.地球物理學(xué)報(bào)，47（6）：1094-1100.

邵秀民，藍(lán)志凌.1995.非均勻各向同性彈性介質(zhì)中地震波傳播的數(shù)值模擬.地球物理學(xué)報(bào)，38（S1）：39-55.

楊頂輝.2002.雙相各向異性介質(zhì)中彈性波方程的有限元解法及波場(chǎng)模擬.地球物理學(xué)報(bào)，45（4）：575-583.

周輝，徐世浙，劉斌等.1997.各向異性介質(zhì)中波動(dòng)方程有限元法模擬及其穩(wěn)定性.地球物理學(xué)報(bào)，40（6）：833-841.

裴正林.2004.任意起伏地表彈性波方程交錯(cuò)網(wǎng)格高階有限差分法數(shù)值模擬.石油地球物理勘探，39（6）：629-634.