崔叢越 張獻(xiàn)兵 王彥賓

（中國北京100871北京大學(xué)地球與空間科學(xué)學(xué)院）

引言

地震波傳播數(shù)值計(jì)算是地震學(xué)中的重要研究領(lǐng)域，在理解復(fù)雜地球介質(zhì)中的地震波傳播過程、地震波反演以及強(qiáng)震地面運(yùn)動(dòng)計(jì)算等方面起著重要作用.隨著地震學(xué)研究的不斷發(fā)展，對地球內(nèi)部橫向非均勻性、復(fù)雜界面形態(tài)分布等復(fù)雜結(jié)構(gòu)的認(rèn)識(shí)不斷精細(xì)化，精確高效的數(shù)值方法是求解復(fù)雜介質(zhì)中地震波傳播頗為有效的手段之一.

常用的地震波傳播數(shù)值模擬方法包括有限差分法、偽譜法、有限元法等.對于大規(guī)模二維模型或三維模型，隨著對地震波場模擬的精度要求不斷提高，通常需要大量的計(jì)算機(jī)內(nèi)存和計(jì)算時(shí)間，為了提高計(jì)算效率，需要利用并行計(jì)算環(huán)境來實(shí)現(xiàn).常見的并行計(jì)算機(jī)架構(gòu)有對稱多處理結(jié)構(gòu)（symmetric multi-processing，簡寫為SMP）、大規(guī)模并行處理計(jì)算機(jī)（massive parallel processing，簡寫為MPP）、機(jī)群（cluster）等，軟件層面多采用信息傳遞接口（message passing interface，簡寫為MPI）.很多研究人員充分利用有限差分算子易于并行的特性，開發(fā)了基于機(jī)群的地震波場并行有限差分模擬方法（王德利等，2007；張明財(cái)?shù)龋?013），同時(shí)適宜于機(jī)群的有限元、偽譜法等并行技術(shù)也得以開發(fā)（Furumura et al，1998；王月英，孫成禹，2006；嚴(yán)九鵬，王彥賓，2008）.

相比于CPU，基于圖形處理器（graphic processing unit，簡寫為GPU）的通用計(jì)算技術(shù)（general purpose GPU，簡寫為GPGPU）起步較晚.早期的GPGPU直接使用圖形開發(fā)接口，將計(jì)算任務(wù)映射為紋理渲染.該方法編程困難，應(yīng)用范圍狹窄，實(shí)用價(jià)值有限.2006年，英偉達(dá)（NVIDIA）公司推出了統(tǒng)一計(jì)算設(shè)備架構(gòu)（compute unified device architecture，簡寫為CUDA）平臺(tái)，極大地降低了編程難度，拓寬了應(yīng)用范圍，使GPGPU被廣泛接受.由于GPU天然的并行性，CUDA在金融、影像處理、天文學(xué)、地球物理學(xué)等領(lǐng)域都有顯著的優(yōu)勢，近年來在地震波場模擬領(lǐng)域的應(yīng)用越來越廣泛.龍桂華等（2010）實(shí)現(xiàn)了基于CUDA的傅里葉偽譜微分矩陣算子的加速，獲得約200倍的加速比；Cai等（2012）利用CUDA實(shí)現(xiàn)了3D交錯(cuò)網(wǎng)格有限差分法加速，取得60—80倍的加速比；單蕊和朱偉（2013）嘗試了基于GPU并行的譜元法地震波場模擬；鄭亮等（2013）實(shí)現(xiàn)了基于GPU的三維有限單元地震波傳播模擬；Okamoto等（2010）和Aoi等（2012）實(shí)現(xiàn)了基于多個(gè)GPU的三維有限差分領(lǐng)域的分解算法，Okamoto等（2012）將其應(yīng)用于2011年日本東北地震的波場傳播計(jì)算.

在地震波傳播數(shù)值計(jì)算方法中，由于空間微分算子的局部特性，有限差分方法并行程度高；而傅里葉偽譜微分算子精度高，但是由于其全局特性，并行程度較低.魏星等（2010）給出了基于有限差分與偽譜的混合方法，即水平方向使用交錯(cuò)網(wǎng)格傅里葉偽譜法，垂直方向使用交錯(cuò)網(wǎng)格4階有限差分法，有效利用了兩種方法各自的優(yōu)點(diǎn).秦艷芳和王彥賓（2012）將其推廣至三維模型并實(shí)現(xiàn)了基于機(jī)群的并行計(jì)算.

由于傅立葉微分算子的全局性，以往通過MPI實(shí)現(xiàn)的并行混合方法只能對垂直方向的有限差分法進(jìn)行并行計(jì)算，性能提升較有限.因此，本文嘗試使用GPU實(shí)現(xiàn)這一過程，

利用CUDA平臺(tái)的cuFFT庫使得水平方向也能利用偽譜法進(jìn)行并行計(jì)算.

1 基本原理

1.1 交錯(cuò)網(wǎng)格混合方法

二維各向同性彈性介質(zhì)中P-SV波的波動(dòng)方程為

本構(gòu)關(guān)系為

式中，ux和uz分別為水平和垂直分量位移；fx和fz分別為x和z方向的體力；σxx，σzz和σxz為應(yīng)力分量；ρ為密度，λ為拉梅常數(shù)，μ為剪切模量.

按照圖1所示進(jìn)行網(wǎng)格劃分.對于實(shí)際地球模型介質(zhì)而言，其參數(shù)垂直方向變化較劇烈，因此在垂直方向采用4階精度交錯(cuò)網(wǎng)格有限差分法，水平方向采用交錯(cuò)網(wǎng)格偽譜法.對于有限差分法，交錯(cuò)網(wǎng)格可在計(jì)算量相當(dāng)?shù)那闆r下提高其計(jì)算精度；而對于偽譜法，交錯(cuò)網(wǎng)格可有效地減小由于傅里葉變換和介質(zhì)性質(zhì)急劇變化所造成的噪聲（魏星等，2010）.

交錯(cuò)網(wǎng)格偽譜法的微分公式為（秦艷芳，王彥賓，2012）

圖1 混合方法網(wǎng)格劃分Fig.1 Arrangement of the staggered grids

其中

4階精度交錯(cuò)網(wǎng)格有限差分法的內(nèi)部節(jié)點(diǎn)計(jì)算公式為（秦艷芳，王彥賓，2012）

自由地表?xiàng)l件通過滿足牽引力為零的條件來實(shí)現(xiàn)，模型的兩側(cè)和下部的人工邊界采用以下吸收邊界條件（Cerjan et al，1985）：

根據(jù)離散模型中的最小網(wǎng)格間距Δmin和最大波速vmax，垂直方向的有限差分法網(wǎng)格劃分要求（Daudt et al，1989）

水平方向的傅里葉偽譜法網(wǎng)格劃分要求（Kosloff，Baysal，1982）

且Δt只有小于最低限定條件才能保證解的穩(wěn)定.

1.2 CUDA編程模型

本文所用的程序均采用CUDA C開發(fā).CUDA C是基于標(biāo)準(zhǔn)C語言的擴(kuò)展，使用NVCC編譯器進(jìn)行編譯，兼容大部分NVIDIA顯卡.在GPU上運(yùn)行的函數(shù)稱為內(nèi)核函數(shù)（kernel），一個(gè)內(nèi)核函數(shù)以線程塊（block）為執(zhí)行單位，不同線程塊間并發(fā)執(zhí)行，但不能通信.每個(gè)線程塊中有一塊共享存儲(chǔ)器，可被線程塊中的所有線程訪問.同一線程塊中的線程可通過同步函數(shù)保證執(zhí)行順序.圖2給出了一個(gè)CUDA程序中線程塊與線程的層級關(guān)系.

圖2 CUDA程序中的線程結(jié)構(gòu)Fig.2 Arrangement of threads in CUDA

內(nèi)核函數(shù)運(yùn)行的過程中，線程塊被分為更小的線程束（warp）來執(zhí)行，線程束包含的線程數(shù)量由硬件決定，目前均為32.線程束的前一半線程或后一半線程組成半線程束（halfwarp）.共享存儲(chǔ)器被分為不同的存儲(chǔ)模塊（bank），每個(gè)存儲(chǔ)器有16或32個(gè)存儲(chǔ)模塊，每個(gè)存儲(chǔ)模塊無法在同一時(shí)間響應(yīng)多個(gè)請求.若位于同一半線程束的線程同時(shí)請求訪問一個(gè)存儲(chǔ)模塊中的不同地址，則會(huì)產(chǎn)生存儲(chǔ)體沖突（bank-conflict），對不同地址的請求只能被串行完成，嚴(yán)重降低訪問速度.因此，需要保證同一半線程束的線程在同一時(shí)間訪問的地址相同，或位于不同的存儲(chǔ)模塊.

除了算法優(yōu)化，CUDA程序的主要性能瓶頸出現(xiàn)在顯存和內(nèi)存間的帶寬上.GPU并行程度足夠高的程序才能彌補(bǔ)兩者數(shù)據(jù)傳輸所造成的延遲，而一個(gè)內(nèi)核函數(shù)中調(diào)用的線程塊和線程的數(shù)量均有一定的限制.本文的測試環(huán)境中，線程塊的數(shù)量上限遠(yuǎn)大于所需值，可忽略線程塊的數(shù)量限制，每個(gè)線程塊最多可有1 024個(gè)線程.

1.3 統(tǒng)一內(nèi)存

統(tǒng)一內(nèi)存是cuda 6.0引入的新特性，它可以通過cudaMallocManaged函數(shù)開辟受管理的內(nèi)存，通過單一指針接受來自CPU與GPU的訪問.由于內(nèi)存和顯存彼此獨(dú)立，GPU程序必須分別在CPU和GPU開辟存儲(chǔ)空間，而且在調(diào)用GPU函數(shù)前需將內(nèi)存數(shù)據(jù)復(fù)制到顯存，并在GPU執(zhí)行完成后將數(shù)據(jù)復(fù)制回來.統(tǒng)一內(nèi)存通過一個(gè)單一指針避免了這個(gè)繁瑣的過程（圖3）.目前的顯卡硬件中，統(tǒng)一內(nèi)存本質(zhì)上是將手動(dòng)復(fù)制數(shù)據(jù)的操作自動(dòng)化，而下一代的麥克斯韋（Maxwell）架構(gòu)顯卡將通過在硬件層面上引入統(tǒng)一虛擬內(nèi)存，從而進(jìn)一步提升性能.

與手動(dòng)內(nèi)存管理相比，統(tǒng)一內(nèi)存主要有以下兩大優(yōu)點(diǎn)：

圖3 統(tǒng)一內(nèi)存示意圖Fig.3 Skematic diagram of unified memory

1）降低編程難度，方便移植現(xiàn)有程序.統(tǒng)一內(nèi)存使得GPU程序不需要引入額外變量，從而可以使原有C程序的代碼基本保持不變；而手動(dòng)管理內(nèi)存至少要引用一個(gè)額外的變量，在多維數(shù)據(jù)等情況下還要引入更多，這使得函數(shù)的參數(shù)形式和寫法均需進(jìn)行較大調(diào)整.特別是在使用鏈表等復(fù)雜數(shù)據(jù)結(jié)構(gòu)的情況下，內(nèi)存與顯存之間的復(fù)制會(huì)變得非常麻煩，而統(tǒng)一內(nèi)存可以避免該問題，減小工作量和降低出錯(cuò)概率.

2）提升數(shù)據(jù)傳輸性能.統(tǒng)一內(nèi)存會(huì)按需在內(nèi)存與顯存間轉(zhuǎn)移數(shù)據(jù)，通過一系列復(fù)雜的算法盡量多地利用帶寬，通常情況下效率會(huì)高于未經(jīng)優(yōu)化的手動(dòng)管理內(nèi)存的程序.同時(shí)，統(tǒng)一內(nèi)存仍然支持直接開辟顯存與異步內(nèi)存?zhèn)鬏敚瑥亩鴿M足對性能要求更高的程序，在降低編程難度的同時(shí)并未犧牲靈活性.

1.4 cuFFT庫

混合方法在水平方向采用傅里葉偽譜法計(jì)算空間微分，因此需要用到快速傅里葉變換.cuFFT庫是CUDA自帶的快速傅里葉變換函數(shù)庫，與串行傅里葉變換相比加速比可達(dá)10倍（NVIDIA，2014）.cuFFT庫包括cuFFT和cuFFTW兩個(gè)子庫，其中cuFFT子庫提供了較高的計(jì)算效率，而cuFFTW子庫為移植使用C語言中快速傅里葉變換庫的現(xiàn)有程序提供了方便.cuFFT庫為不同規(guī)模的輸入提供高度優(yōu)化的算法，可以選擇單精度或雙精度，支持批量變換.

cuFFT庫將傅里葉變換分為3大類：復(fù)數(shù)域→復(fù)數(shù)域（C2C），實(shí)數(shù)域→復(fù)數(shù)域（R2C），復(fù)數(shù)域→實(shí)數(shù)域（C2R）.當(dāng)輸入數(shù)據(jù)為實(shí)數(shù)域的采樣點(diǎn)時(shí)，輸出數(shù)據(jù)滿足Xk＝X＊N－K，選擇以R2C的方式進(jìn)行變換可利用此對稱性，將計(jì)算量降至最低.

調(diào)用cuFFT庫的一般步驟為：

1）利用cufftPlan1D／cufftPlan2D／cufftPlan3D創(chuàng)建一維、二維或三維變換計(jì)劃，或利用cufftPlayMany函數(shù)創(chuàng)建一個(gè)批量變換計(jì)劃；

2）根據(jù)變換類型，調(diào)用cufftExecC2C，cufftExecR2C或cufftExecC2R函數(shù)執(zhí)行變換；

3）調(diào)用cufftDestroy函數(shù)銷毀計(jì)劃，釋放資源.

2 計(jì)算流程實(shí)現(xiàn)

2.1 混合方法計(jì)算步驟

GPU程序的并行部分直接影響執(zhí)行效率，因此要盡可能將相互獨(dú)立的計(jì)算步驟放到GPU中執(zhí)行，同時(shí)將并行執(zhí)行效率低的部分交給CPU執(zhí)行.其中，微分、應(yīng)力及位移的計(jì)算部分由GPU完成，數(shù)據(jù)導(dǎo)入導(dǎo)出和循環(huán)控制部分由CPU完成.設(shè)水平方向網(wǎng)格數(shù)為m，垂直方向網(wǎng)格數(shù)為n，線程塊數(shù)與每個(gè)線程塊中線程數(shù)的比值為x，求解方程（1）、（2）的具體步驟如下：

1）開辟內(nèi)存、顯存，讀取數(shù)據(jù)并使之初始化；

2）利用傅里葉偽譜法求出?ux／?x和?uy／?x；

3）利用有限差分法求出?ux／?y和?uy／?y；

4）根據(jù)式（2）求出σxx，σxy和σyy；

5）利用傅里葉偽譜法求出?σxx／?x和?σxy／?x；

6）利用有限差分法求出?σxy／?y和?σyy／?y；

7）利用式（1）計(jì)算ux和uy.

2.2 邊界條件

式（7）中取α＝0.015，I0＝20.由于垂直方向采用有限差分法，無需對自由表面進(jìn)行處理.計(jì)算得到大小為m×n的系數(shù)數(shù)組，由于其被多次使用且在計(jì)算過程中保持不變，適合將其放入紋理存儲(chǔ)器.同樣適合放入紋理存儲(chǔ)器的數(shù)據(jù)還有密度等介質(zhì)參數(shù).

2.3 有限差分法

有限差分法對兩個(gè)方向求空間微分在同一內(nèi)核函數(shù)并行執(zhí)行，調(diào)用的線程塊數(shù)為2mx，線程數(shù)為n／x.根據(jù)實(shí)際測試，取x＝m／64時(shí)計(jì)算效率最高.實(shí)現(xiàn)方式與串行代碼完全類似，只需將循環(huán)改為分配給不同的線程，線程的編號(hào)與循環(huán)中的步數(shù)一一對應(yīng).運(yùn)行完成后需執(zhí)行一次cudaDeviceSynchronize函數(shù)，將顯存中的數(shù)據(jù)同步至內(nèi)存.

2.4 傅里葉偽譜法

計(jì)算中需要對2n個(gè)數(shù)組進(jìn)行傅里葉正變換和反變換，可用cufftPlanMany函數(shù)創(chuàng)建計(jì)劃批量進(jìn)行.由于輸入的數(shù)據(jù)均為實(shí)數(shù)，理論上正變換采用R2C類型計(jì)算效率較高，但R2C必須開辟一個(gè)實(shí)數(shù)和一個(gè)復(fù)數(shù)共兩個(gè)數(shù)組.由于cuFFT支持就地變換（變換后的結(jié)果直接保存在原數(shù)組），使用C2C類型只需開辟一個(gè)復(fù)數(shù)數(shù)組.因此，在顯存比較有限的情況下，使用C2C類型效率反而略高.

3 計(jì)算效率測試分析

為了測試上述流程的計(jì)算效率，本文采用二維均勻介質(zhì)模型，對比CPU與GPU版本的運(yùn)行時(shí)間.運(yùn)行時(shí)間以循環(huán)開始作為起點(diǎn)，不計(jì)入數(shù)據(jù)讀取、初始化及導(dǎo)出的時(shí)間.測試環(huán)境如表1所示.使用clock函數(shù)計(jì)量時(shí)間，波動(dòng)一般在1%以內(nèi).

CPU與GPU的運(yùn)行結(jié)果如表2所示，圖4給出了加速比與網(wǎng)格數(shù)的關(guān)系.通過測試可見，使用GPU可取得明顯的加速效果.隨著網(wǎng)格數(shù)的增加，并行程度提高，加速也有明顯的增加.但由于顯存的限制，增幅會(huì)逐漸減緩.

表1 測試環(huán)境Table 1 Test environment

表2 計(jì)算耗時(shí)表Table 2 Elapsed time

4 應(yīng)用實(shí)例

為了驗(yàn)證本文所給出的GPU并行方法的有效性，我們選取蘭州盆地的一個(gè)二維剖面，計(jì)算點(diǎn)源激發(fā)的地震波在其中的傳播過程.蘭州盆地模型引自嚴(yán)武建等（2013），沉積盆地上部包含一個(gè)深度達(dá)393m的第四紀(jì)沉積層和一個(gè)深度達(dá)3.0km的第三系沉積層，底部為上地殼.模型寬70km，深35km.震源為雙力偶點(diǎn)源，矩張量分量Mxx＝1.0×1012N·m，Mzz＝－1.0×1012N·m.網(wǎng)絡(luò)總數(shù)為2 048×1 024，水平和垂直方向的網(wǎng)格大小均為34.2m.震源時(shí)間函數(shù)的寬度為0.1s，時(shí)間間隔為0.001s，共計(jì)算3萬時(shí)間步.

圖4 不同網(wǎng)格數(shù)的加速比Fig.4 Speedup in various number of grids

圖5 蘭州盆地模型地表水平分量ux和垂直分量uz的合成地震圖Fig.5 Horizontal and vertical component of synthetic waveforms in Lanzhou basin model

圖5 給出了地表垂直分量和水平分量的地震圖.基巖地表的震相簡單，可以看到直達(dá)P、SV波和上地殼底部的反射波.當(dāng)直達(dá)波傳播進(jìn)入沉積盆地時(shí)，在30km距離上開始看到直達(dá)波及其在盆地邊界上產(chǎn)生的轉(zhuǎn)換波，同時(shí)可以看到由這些波相干疊加而形成的面波，其振幅大于直達(dá)波和轉(zhuǎn)換波.當(dāng)波傳播到低速的第四系沉積層時(shí)，在其邊界上又產(chǎn)生各種波的透射與轉(zhuǎn)換，使得低速層內(nèi)的波場比第三系沉積層更為復(fù)雜，最終形成振幅更大的面波.圖6給出了4，8，10，14，22和30s的波場快照，可以看出P波和SV波在基巖中的傳播、盆地邊界所發(fā)生的震相轉(zhuǎn)換以及面波的形成過程.

圖6 不同時(shí)刻時(shí)地震波傳播的波場快照.圖中紅色和綠色分別表示P和SV波Fig.6 Wavefield snapshots at t＝4，8，10，14，22，30s，where red and green represent P and SV waves，respectively

由上述結(jié)果可以看出，低速沉積盆地造成了地震波明顯的放大效應(yīng)；由于各種震相的相干疊加以及邊緣面波的作用，盆地內(nèi)地震波的振幅會(huì)增大，地震波傳播的持續(xù)時(shí)間會(huì)變長.這些模擬結(jié)果與嚴(yán)武建等（2013）的結(jié)果完全一致，從另一方面說明了本文的GPU并行模擬方法在實(shí)際應(yīng)用中的可行性.

5 討論與結(jié)論

有限差分／偽譜混合方法結(jié)合了有限差分和偽譜兩種方法的優(yōu)點(diǎn)，保持了較高的計(jì)算精度，是計(jì)算二維非均勻介質(zhì)地震波傳播的有效方法.同時(shí)，該方法能夠有效處理自由表面及垂直方向的不連續(xù)面.由于水平方向與垂直方向的計(jì)算相對獨(dú)立，混合方法可以較容易地推廣至三維，擴(kuò)展其應(yīng)用范圍.本文實(shí)現(xiàn)了該混合方法的GPU編程，顯著提高了計(jì)算速度.實(shí)際算例表明，在處理好顯存上限的前提下，基于GPU的混合方法程序可完全替代其CPU版本.隨著顯卡硬件的發(fā)展與開發(fā)工具的不斷完善，基于GPU的模擬程序?qū)⒃絹碓狡毡?

由于cuFFT庫的局限性，目前傅里葉變換只能從CPU中調(diào)用，因此與純有限差分或純偽譜法相比，混合方法增加了額外的內(nèi)存拷貝時(shí)間，也無法實(shí)現(xiàn)偽譜與有限差分并行執(zhí)行.不過NVIDIA公司已在開發(fā)支持從內(nèi)核函數(shù)調(diào)用的cuFFT版本，所以此問題有望解決.

在利用混合方法計(jì)算三維模型時(shí)，單個(gè)GPU的性能和顯存很難滿足計(jì)算需求.因此可考慮采用MPI／CUDA混合方法，根據(jù)顯存大小或線程塊數(shù)量限制，沿垂直方向?qū)⒛Ｐ头殖啥鄠€(gè)區(qū)域，在GPU集群上進(jìn)行計(jì)算.在MVAPICH等支持CUDA的MPI環(huán)境中，可以直接發(fā)送和接收GPU數(shù)據(jù)，從而方便高效地完成計(jì)算.在三維問題中，數(shù)據(jù)拷貝帶來的延遲影響將更加突出，因此優(yōu)化數(shù)據(jù)結(jié)構(gòu)，并且更精細(xì)地控制內(nèi)存與顯存之間的傳輸將是本文所用程序亟需改進(jìn)的地方.

龍桂華，李小凡，江東輝.2010.基于交錯(cuò)網(wǎng)格Fourier偽譜微分矩陣算子的地震波場模擬GPU加速方案［J］.地球物理學(xué)報(bào)，53（12）：2964-2971.

Long G H，Li X F，Jiang D H.2010.Accelerating seismic modeling with staggered-grid Fourier pseudo-spectral differentiation matrix operator method on graphics processing unit［J］.Chinese Journal of Geophysics，53（12）：2964-2971（in Chinese）.

秦艷芳，王彥賓.2012.地震波傳播的三維偽譜和高階有限差分混合方法并行模擬［J］.地震學(xué)報(bào)，34（2）：147-156.

Qin Y F，Wang Y B.2012.Three-dimensional parallel hybrid PSM／FDM simulation of seismic wave propagation［J］.Acta Seismologica Sinica，34（2）：147-156（in Chinese）.

單蕊，朱偉.2013.基于GPU并行的譜元法地震波數(shù)值模擬［J］.中國煤炭地質(zhì)，25（3）：52-54，62.

Shan R，Zhu W.2013.Parallel SEM seismic wave simulation based on GPU［J］.Coal Geology of China，25（3）：52-54，62（in Chinese）.

王德利，雍運(yùn)動(dòng)，韓立國，廉玉廣.2007.三維粘彈介質(zhì)地震波場有限差分并行模擬［J］.西北地震學(xué)報(bào)，29（1）：30-34.

Wang D L，Yong Y D，Han L G，Lian Y G.2007.Parallel simulation of finite difference for seismic wavefield modeling in 3-D viscoelastic media［J］.Northwestern Seismological Journal，29（1）：30-34（in Chinese）.

王月英，孫成禹.2006.彈性波動(dòng)方程數(shù)值解的有限元并行算法［J］.中國石油大學(xué)學(xué)報(bào)：自然科學(xué)版，30（5）：27-30.

Wang Y Y，Sun C Y.2006.Finite element parallel algorithm for numerical solution of elastic wave equation［J］.Journal of China University of Petroleum，30（5）：27-30（in Chinese）.

魏星，王彥賓，陳曉非.2010.模擬地震波場的偽譜和高階有限差分混合方法［J］.地震學(xué)報(bào)，32（4）：392-400.

Wei X，Wang Y B，Chen X F.2010.Hybrid PSM／FDM method for seismic wavefield simulation［J］.Acta Seismologica Sinica，32（4）：392-400（in Chinese）.

嚴(yán)九鵬，王彥賓.2008.重疊區(qū)域偽譜法計(jì)算非均勻地球介質(zhì)地震波傳播［J］.地震學(xué)報(bào)，30（1）：47-58.

Yan J P，Wang Y B.2008.Modeling seismic wave propagation in heterogeneous medium using overlap domain pseudospectral method［J］.Acta Seismologica Sinica，30（1）：47-58（in Chinese）.

嚴(yán)武建，王彥賓，石玉成.2013.基于偽譜和有限差分混合方法的蘭州盆地強(qiáng)地面運(yùn)動(dòng)二維數(shù)值模擬［J］.地震工程學(xué)報(bào)，35（2）：302-310.

Yan W J，Wang Y B，Shi Y C.2013.2Dsimulation of the strong ground motion in Lanzhou basin with hybrid PSM／FDM method［J］.China Earthquake Engineering Journal，35（2）：302-310（in Chinese）.

張明財(cái)，熊章強(qiáng)，張大洲.2013.基于 MPI的三維瑞雷面波有限差分并行模擬［J］.石油物探，52（4）：354-362.

Zhang M C，Xiong Z Q，Zhang D Z.2013.3Dfinite difference parallel simulation of Rayleigh wave based on message passing interface［J］.Geophysical Prospecting for Petroleum，52（4）：354-362（in Chinese）.

鄭亮，張懷，尹鳳玲，石耀霖.2013.基于GPU的大規(guī)模三維地震波傳播有限元程序?qū)崿F(xiàn)［C］∥中國地球物理2013：第七分會(huì)場論文集.昆明，2013年10月13日.

Zheng L，Zhang H，Yin F L，Shi Y L.2013.Finite element program implementation of large scale seismic wave propagation modeling based on GPU［C］∥Chinese Geophysics 2013：Proceedings of the 7th Branch.Kunming，October 13，2013.

Aoi S，Nishizawa N，Aoki T.2012.Large scale simulation of seismic wave propagation using GPGPU［C］∥15th World Conference on Earthquake Engineering.Lisbon，Portugal，September 24to 28，2012.

Cai C，Chen H Q，Deng Z，Chen D，Khan S U，Zeng K，Wu M X.2012.GPGPU-aided 3Dstaggered-grid finite-difference seismic wave modeling［C／OL］∥12th International Conference on Scalable Computing and Communications（ScalCom）.［2014-05-12］.http：∥newtheme.hgpu.org／？p＝8610.

Cerjan C，Kosloff D，Kosloff R，Reshef M.1985.A nonreflecting boundary condition for discrete acoustic and elastic wave equations［J］.Geophysics，50（4）：705-708.

Daudt C R，Braille L W，Nowack R L，Chiang C S.1989.A comparison of finite-difference and Fourier method calculations of synthetic seismograms［J］.Bull Seismol Soc Am，79（4）：1210-1230.

Furumura T，Kennett B L N，Takenaka H.1998.Parallel 3-D pseudospectral simulation of seismic wave propagation［J］.Geophysics，63（1）：279-288.

Kosloff D，Baysal E.1982.Forward modeling by a Fourier method［J］.Geophysics，47（10）：1402-1412.

NVIDIA.2014.CUFFTLibraryUser’s Guide［M／OL］.［2014-04-26］.http：∥docs.nvidia.com／cuda／cufft／.

Okamoto T，Takenaka H，Nakamura T，Aoki T.2010.Accelerating large-scale simulation of seismic wave propagation by multi-GPUs and three-dimensional domain decomposition［J］.Earth Planets Space，62（12）：939-942.

Okamoto T，Takenaka H，Nakamura T，Aoki T.2012.Large-scale simulation of seismic wave propagation of the 2011 Tohoku-Oki M9earthquake［C］∥Proceedings of the International Symposium on Engineering Lessons Learned from the 2011 Great East Japan Earthquake.Tokyo：Japan Association for Earthquake Engineering：349-360.