邸海濱 許力生

（中國(guó)北京100081中國(guó)地震局地球物理研究所）

層狀介質(zhì)中有限震源引起的地面運(yùn)動(dòng)的計(jì)算對(duì)點(diǎn)源情形的拓展＊>

邸海濱 許力生

（中國(guó)北京100081中國(guó)地震局地球物理研究所）

從分層均勻介質(zhì)中地震波傳播的基本理論出發(fā)，參考已有的利用廣義反射透射系數(shù)矩陣方法計(jì)算分層均勻介質(zhì)中點(diǎn)源引起的地面運(yùn)動(dòng)的方法和程序，考慮了任意幾何特征的有限震源及有限震源的震源機(jī)制隨時(shí)間和空間發(fā)生變化的可能性，并考慮了并行計(jì)算的發(fā)展方向，本文對(duì)點(diǎn)源情況下的計(jì)算流程進(jìn)行了改進(jìn)，重新設(shè)計(jì)編寫(xiě)了計(jì)算程序GRTMatSyn，使其不但適應(yīng)于點(diǎn)源的情形，也適應(yīng)于任意幾何形狀的、有限的、震源機(jī)制隨時(shí)間和空間變化的有限震源的情形.為了驗(yàn)證該程序的可靠性和計(jì)算效率，設(shè)計(jì)了必要的有限震源模型和觀測(cè)點(diǎn)分布，計(jì)算了地面運(yùn)動(dòng)，并與已有的被廣泛認(rèn)可的計(jì)算有限平面斷層引起的地面運(yùn)動(dòng)的程序CompSyn的計(jì)算結(jié)果進(jìn)行了對(duì)比分析.結(jié)果表明，GRTMatSyn的計(jì)算結(jié)果可靠、計(jì)算效率較高.

廣義反透射系數(shù)矩陣 層狀介質(zhì) 有限震源 地面運(yùn)動(dòng)

Abstract：From fundamental theory of seismic wave propagation in stratified media，with reference to the techniques and applications already developed for calculating point－source－caused ground motion using generalized reflection and transmission matrices method，and with consideration of arbitrary geometry and tempo－spatially variable focal mechanisms of finite sources，as well as development trend of parallel computation，we improved computing－flow in the pointsource case and developed a new program，GRTMatSyn，to be suitable for both point sources and finite sources with arbitrary geometries and tempo－spatially variable mechanisms.Numerical tests were carried out with necessary models offinite sources and observation points for reliability and computing efficiency of the application.The results were compared with those from CompSyn，a widely recognized program for calculating ground motion caused by a finite plane fault.The comparison indicated that the application developed in this study is reliable and rather more efficient.

Key words：generalized reflection and transmission matrices method；stratified media；finite source；ground motion

引言

地震災(zāi)害直接與地震引起的地面運(yùn)動(dòng)有關(guān).及時(shí)掌握震中區(qū)地面運(yùn)動(dòng)的時(shí)空分布，無(wú)疑對(duì)了解災(zāi)情，進(jìn)而采取救援措施至關(guān)重要.理想的作法是利用高密度觀測(cè)站與高速暢通的通訊，實(shí)時(shí)監(jiān)測(cè)和傳遞地震現(xiàn)場(chǎng)的地面運(yùn)動(dòng)圖像.但實(shí)際上卻存在難以克服的困難，高密度的臺(tái)站分布幾乎不可能實(shí)現(xiàn)，通訊高速暢通也不可能保證，尤其在地震發(fā)生的時(shí)候.所以，權(quán)宜之計(jì)是利用遙遠(yuǎn)而稀疏的臺(tái)站觀測(cè)資料獲得地震震源的時(shí)空過(guò)程，然后利用震源的時(shí)空過(guò)程再現(xiàn)地震現(xiàn)場(chǎng)的地面運(yùn)動(dòng).

經(jīng)過(guò)近30年的努力，利用遠(yuǎn)場(chǎng)觀測(cè)資料獲取地震震源的時(shí)空破裂過(guò)程已獲得了長(zhǎng)足的進(jìn)展，震后5小時(shí)左右獲取震源的破裂圖像已成為可能（張勇，2008；張勇等，2008，2010）.震源的時(shí)空破裂圖像在地震應(yīng)急響應(yīng)中發(fā)揮的重要作用也日益得到體現(xiàn).例如，2008年四川汶川MS8.0地震（張勇等，2008），2010年青海玉樹(shù)MS7.1地震（張勇等，2010）等.然而，地震現(xiàn)場(chǎng)的地面運(yùn)動(dòng)圖像，包括地面位移圖像、地面速度圖像或地面加速度圖像，對(duì)于應(yīng)急決策更直觀.因此，利用地震模型快速計(jì)算近源的地面運(yùn)動(dòng)已成為現(xiàn)今社會(huì)發(fā)展的必然要求，也符合地震學(xué)學(xué)科發(fā)展的方向.

計(jì)算理論地震圖的經(jīng)典方法主要有基于射線理論的廣義射線方法（Helmberger，1968；Gilbert，Helmberger，1972；Chapman，1974）和基于波動(dòng)理論的反射率方法（Kennett，Kerry，1979；Bouchon，1981；Kennett，1980，1983；Chen，1999）.廣義射線方法需要指定感興趣的射線路徑，而且隨著射線路徑的復(fù)雜化，計(jì)算成本也逐漸增加，因此多用于遠(yuǎn)場(chǎng)較低頻的地震波計(jì)算.由于需要指定具體的射線路徑，所以，通常的計(jì)算結(jié)果并不是介質(zhì)的完全響應(yīng).比較而言，反射率方法卻不需要指定具體的傳播路徑，計(jì)算結(jié)果自動(dòng)包含介質(zhì)的完全響應(yīng)，因此更適合計(jì)算近斷層的地面運(yùn)動(dòng).

本文從分層均勻介質(zhì)中地震波傳播的基本理論出發(fā)（Kennett，1983），參考已有的利用廣義反射透射系數(shù)矩陣方法計(jì)算分層均勻介質(zhì)中點(diǎn)源引起的地面運(yùn)動(dòng)的方法和程序（李旭，1993；李旭，陳運(yùn)泰，1996；劉超，2011），針對(duì)任意幾何特征的有限震源，考慮震源機(jī)制在震源過(guò)程中的時(shí)空變化，順應(yīng)并行計(jì)算的發(fā)展方向，對(duì)計(jì)算流程進(jìn)行優(yōu)化，重新設(shè)計(jì)編寫(xiě)計(jì)算程序GRTMatSyn，并利用數(shù)值試驗(yàn)對(duì)該程序的可靠性和計(jì)算效率給出評(píng)價(jià).

1 分層均勻介質(zhì)中點(diǎn)源引起的地面運(yùn)動(dòng)

1.1 地面運(yùn)動(dòng)的表示

關(guān)于分層均勻介質(zhì)中點(diǎn)源引起的地震波傳播理論已經(jīng)相當(dāng)成熟（Kennett，Kerry，1979；Kennett，1980，1983；Yao，Harkrider，1983；Chen，1999）.下面只簡(jiǎn)要給出與計(jì)算地面運(yùn)動(dòng)直接有關(guān)的核心公式.

如圖1所示，在分層均勻介質(zhì)中，震源輻射的地震波總是可以分成兩部分：上行波（自震源向上傳播的波）和下行波（自震源向下傳播的波）.為了描述問(wèn)題的方便，在震源深度位置人為設(shè)置一個(gè)界面（這個(gè)界面其實(shí)不存在），使其平行于其它介質(zhì)界面，并以經(jīng)過(guò)震源垂直分層界面向下的方向?yàn)橹鴺?biāo)的中心軸建立柱坐標(biāo)系.

圖1 分層均勻介質(zhì)及坐標(biāo)構(gòu)架（Kennett，1983）z0表示自由表面，zL表示介質(zhì)底界面，zS表示震源處的虛界面.U和D分別表示上行波場(chǎng)和下行波場(chǎng)Fig.1 Stratified media and coordinate frame（Kennett，1983）z0denotes free surface，zL，bottom interface，and zS，imaginary interface of source.U and D represent up－going and down－going wave－fields，respectively

在柱坐標(biāo)系中，分層均勻介質(zhì)中點(diǎn)源引起的地面運(yùn)動(dòng)的面諧矢量位移解w0可以表示為（Kennett，1983；李旭，1993；李旭，陳運(yùn)泰，1996）

圖2 廣義反射透射系數(shù)矩陣（，，和）幾何解釋z0表示自由表面，zL表示介質(zhì)底界面，zS表示震源界面Fig.2 Geometrical illustration of the generalized reflection and transmission matrices，，andz0denotes free surface，zL，bottom interface，and zS，imaginary interface of source

由式（1）可以看出，面諧矢量位移解w0依賴(lài)于自由表面的接收系數(shù)矩陣、反射系數(shù)矩陣、分層介質(zhì)中上行波和下行波的廣義反透射系數(shù)，以及上行波與下行波地震波矢量間斷.如果通過(guò)這些量得到面諧矢量位移解w0，那么很容易通過(guò)Fourier－Hankel反變換得到時(shí)空域點(diǎn)源引起的地面運(yùn)動(dòng)解.可見(jiàn)，計(jì)算地面運(yùn)動(dòng)的關(guān)鍵在于計(jì)算面諧矢量位移解w0，而計(jì)算面諧矢量位移解w0的關(guān)鍵則在于計(jì)算廣義反透射系數(shù)矩陣和表征震源作用的地震波矢量間斷.因此，提高計(jì)算有限震源引起的地面運(yùn)動(dòng)效率的途徑之一是優(yōu)化廣義反透射系數(shù)矩陣和表示震源作用的地震波矢量間斷的計(jì)算流程.

1.2 廣義反透射系數(shù)矩陣的表示

如圖3所示，用界面A和界面B表示相鄰的兩個(gè)界面，用界面B和C表示任意多層介質(zhì)之間的界面，用RD（z，z）與TD（z，z）分別表示地震波矢量通過(guò)界面A的下行波的反射系數(shù)和透射系數(shù)矩陣，用RU（z，z）與TU（z，z）分別表示地震波矢量通過(guò)界面A的上行波的反射系數(shù)和透射系數(shù)矩陣.具體某界面的反透射系數(shù)的計(jì)算比較簡(jiǎn)單（詳見(jiàn)附錄B）.相應(yīng)地，界面B與C之間的廣義反射透射系數(shù)矩陣用，，和來(lái)表示，那么z與z之間的廣義反射透射系數(shù)矩陣，，和由下式確定.

圖3 地震波矢量在分層介質(zhì)中傳播的幾何解釋以水平實(shí)線表示文中考慮的3個(gè)分界面zA，zB和zC；虛線表示任意界面；上標(biāo)“－”與“＋”分別表示分界面的上下邊界；帶箭頭的實(shí)線表示地震波矢量Fig.3 Geometrical illustration of the wave vectors in stratified media Solid horizontal lines represent the three interfaces zA，zB，and zC，involved in the text.Discontinued lines denote other interfaces，and superscript“－”and“＋”mean upper and lower boundaries of interface.Solid line with arrow represents wave vector

式中

式中，EABD＝diag｛exp［iωqα（zB－zA）］exp［iωqβ（zB－zA）］｝表示地震波在zA＜z＜zB之間發(fā)生的相位變化.其中，α與β分別表示P波與S波的波速；ω表示圓頻率；p表示慢度，qα＝（α－2－p2）1/2，qβ＝（β－2－p2）1/2，且Imωqα≥0，Imωqβ≥0.

1.3 地震波矢量間斷的表示

在式（1）中，震源作用以地震波矢量間斷∑U（zS）和∑D（zS）來(lái)描述（Kennett，1983），下列各式把描述一般點(diǎn)源的地震矩張量和地震波矢量間斷聯(lián)系起來(lái).

對(duì)于P－SV波

對(duì)于SH波

1）當(dāng)m＝0時(shí)

2）當(dāng)m＝±1時(shí)

3）當(dāng)m＝±2時(shí)

式中

其中，i表示虛數(shù)單位，ω表示圓頻率，p表示慢度，α與β分別表示P波與S波的波速，qα＝（α－2－p2）1/2，qβ＝（β－2－p2）1/2，且Imωqα≥0，Imωqβ≥0，εα＝（2ρqα）－1/2，εβ＝（2ρqβ）－1/2，m為貝塞爾函數(shù)的階數(shù)，ρ表示介質(zhì)密度.

2 有限震源引起的地面運(yùn)動(dòng)的計(jì)算

有限震源引起的地面運(yùn)動(dòng)，總是可以通過(guò)多個(gè)點(diǎn)源引起的地面運(yùn)動(dòng)的適當(dāng)疊加得到（姚振興，鄭天愉，1985）.在計(jì)算有限震源的地面運(yùn)動(dòng)時(shí)，我們必須考慮有限震源的一般性.因?yàn)椋瑢?shí)際地震的震源幾何形狀可能是多樣的，地震過(guò)程中震源機(jī)制可能隨時(shí)間和空間發(fā)生變化.考慮到幾何形狀的一般性，我們并不把震源局限在一個(gè)平面內(nèi)；考慮到震源機(jī)制的變化，我們使用地震矩張量表示震源.

2.1 考慮不同震源深度的廣義反透射系數(shù)矩陣的計(jì)算

由式（2）和式（3）可以看出，介質(zhì)內(nèi)部任意層之間的廣義反透射系數(shù)依賴(lài)于其介質(zhì)結(jié)構(gòu).如果它們的幾何分層和各層的物理屬性不發(fā)生變化，那么對(duì)應(yīng)的廣義反透射系數(shù)矩陣不會(huì)發(fā)生變化.對(duì)于有限震源情形，通常只有震源層界面發(fā)生變化，上行波和下行波的廣義反透射系數(shù)矩陣的變化僅發(fā)生在震源層與其鄰層之間.所以，對(duì)于多個(gè)震源界面的情況，其它層之間的廣義反透射系數(shù)矩陣無(wú)需重復(fù)計(jì)算.

由于各反透射系數(shù)矩陣的計(jì)算都類(lèi)同，所以，我們借助于圖4和圖5僅闡述的計(jì)算過(guò)程.為了闡述問(wèn)題的方便，如圖4所示，我們假設(shè)僅有2個(gè)不同深度的震源，震源界面分別為zS1和zS2，兩者之間可以存在多個(gè)介質(zhì)層.

2.2 考慮震源機(jī)制可變的地震波矢量間斷的計(jì)算

式（4）和式（5）建立了地震波矢量間斷與描述一般點(diǎn)源的地震矩張量之間的聯(lián)系.我們知道，任意矩張量解或者任何一種震源機(jī)制解都可以表示為6個(gè)基本矩張量的加權(quán)和，即

式中，Mj（j＝1，6）表示6個(gè)基本矩張量，具體表示如下：

若以∑Uj（j＝1，6）表示震源處6個(gè)基本矩張量引起的地震波矢量間斷的上行波場(chǎng)，則

同理

因此，我們可以預(yù)先計(jì)算6個(gè)基本矩張量對(duì)應(yīng)的地震波矢量間斷（∑U（zS）和∑D（zS）），任意機(jī)制的震源引起的地震波矢量間斷都可以由此基本的地震波矢量間斷合成.很容易想到，如果有限震源的震源機(jī)制少于6種，采取這樣的技術(shù)路線是不劃算的；如果多于6種，則是節(jié)省時(shí)間的.然而實(shí)際情況是，一次破壞性地震的有限震源往往由遠(yuǎn)大于6的點(diǎn)源構(gòu)成.所以，在一般情況下，采用上面的技術(shù)路線是經(jīng)濟(jì)實(shí)用的.

2.3 有限震源引起的地面運(yùn)動(dòng)的合成

從上面討論可以看出，6個(gè)基本矩張量對(duì)應(yīng)6組地震波矢量間斷（∑Uj（zS）和∑Dj（zS），j＝1，6），每個(gè)震源深度對(duì)應(yīng)一套廣義反射透射系數(shù)矩陣R0SD，T0SU，RSLD和R0SU.將其代入式（1），即可得到6種基本震源引起的面諧矢量位移解w0j（j＝1，6），則任意矩張量引起的w0可以表示為

有限的震源可以看作多個(gè)點(diǎn)源的組合，每個(gè)點(diǎn)源有它自己的深度和震源機(jī)制.通過(guò)上述技術(shù)途徑可以很快得到每個(gè)點(diǎn)源的面諧矢量位移解，然后將這些面諧矢量位移解經(jīng)Fourier－Hankel反變換得到時(shí)間－空間域的位移，即地震位移，最后將所有點(diǎn)源在觀測(cè)點(diǎn)的位移疊加即可得到有限震源引起的地面運(yùn)動(dòng).概而言之，計(jì)算有限震源的地面運(yùn)動(dòng)的計(jì)算流程如下：

1）計(jì)算層狀介質(zhì)中各界面的反射透射系數(shù)矩陣、自由表面的反射系數(shù)矩陣和接收系數(shù)矩陣.

2）計(jì)算6個(gè)基本矩張量引起的地震波矢量間斷.

3）計(jì)算不同震源深度點(diǎn)源對(duì)應(yīng)的廣義反射透射系數(shù)矩陣.

4）將地震波矢量間斷、廣義反射透射系數(shù)矩陣、自由表面的反射系數(shù)矩陣和接收系數(shù)矩陣組合，得到不同震源深度點(diǎn)源的6個(gè)基本矩張量引起的面諧矢量位移解.

5）根據(jù)給定的有限震源的點(diǎn)源的深度和震源機(jī)制，調(diào)用4）的計(jì)算結(jié)果，生成此點(diǎn)源引起的面諧矢量位移解，并對(duì)其進(jìn)行Fourier－Hankel反變換，得到時(shí)－空域的介質(zhì)響應(yīng)或格林函數(shù).

6）將5）生成的介質(zhì)響應(yīng)與給定的震源時(shí)間函數(shù)褶積，生成考慮時(shí)間過(guò)程的點(diǎn)源引起的地面運(yùn)動(dòng).

7）疊加所有點(diǎn)源在給定觀測(cè)點(diǎn)引起的地面運(yùn)動(dòng)，生成有限震源引起的地面運(yùn)動(dòng).

2.4 并行計(jì)算的考慮

為進(jìn)一步提高程序的計(jì)算效率，我們對(duì)計(jì)算過(guò)程中的兩個(gè)環(huán)節(jié)采用了并行計(jì)算設(shè)計(jì).第一個(gè)環(huán)節(jié)是廣義反透射系數(shù)矩陣，，和的計(jì)算.該環(huán)節(jié)包括自由表面與震源之間的廣義反透射系數(shù)矩陣，和的計(jì)算及震源與介質(zhì)底界面之間廣義反射系數(shù)矩陣的計(jì)算兩部分.二者并不發(fā)生聯(lián)系，因此，采用了并行計(jì)算.第二個(gè)環(huán)節(jié)是有限震源的地面運(yùn)動(dòng)計(jì)算.在這個(gè)環(huán)節(jié)上有兩種可行的并行方案：第一種是把有限的震源分成多組而觀測(cè)點(diǎn)為一組；第二種是把觀測(cè)點(diǎn)分為多組而震源為一組.考慮到在通常情況下觀測(cè)點(diǎn)數(shù)目遠(yuǎn)多于震源的數(shù)目，我們采用了后者，即并行計(jì)算同一有限震源在不同觀測(cè)點(diǎn)的地面運(yùn)動(dòng).

3 可靠性與效率的測(cè)試

為了測(cè)試程序的可靠性并檢驗(yàn)其計(jì)算效率，我們選擇相同的介質(zhì)模型和震源模型，分別利用GRTMatSyn和被廣泛使用的程序CompSyn（Spudich，Xu，2003）計(jì)算給定點(diǎn)的地面運(yùn)動(dòng)，對(duì)二者的計(jì)算結(jié)果和計(jì)算時(shí)間進(jìn)行對(duì)比分析.CompSyn采用Olson等（1984）的離散波數(shù)有限元方法計(jì)算格林函數(shù)，并采用Spudich和Archuleta（1987）提出的數(shù)字技術(shù)實(shí)現(xiàn)斷層面上的位移積分.CompSyn的特點(diǎn)是計(jì)算的格林函數(shù)反映了介質(zhì)的完全響應(yīng)（包括P波、S波、面波以及近場(chǎng)項(xiàng)），并適用于有限的平面斷層.然而，由于采用了有限元算法，所以計(jì)算量隨著觀測(cè)點(diǎn)離震源的距離增加而增加；由于時(shí)代的局限性，并沒(méi)有考慮并行計(jì)算；另外，由于局限于平面斷層，所以組成有限斷層的點(diǎn)源的走向和傾角不能變化.我們之所以使用CompSyn與GRTMatSyn進(jìn)行比較，是因?yàn)槎哂?jì)算的格林函數(shù)都包含了介質(zhì)的完全響應(yīng)，且都可以計(jì)算有限平面斷層引起的地面運(yùn)動(dòng).

由于CompSyn采用震源坐標(biāo)，所以我們不防定義Z軸垂直向下為正，N軸向北為正，E軸向東為正，坐標(biāo)系遵循右手法則.測(cè)試斷層如圖6所示，走向?yàn)?°，傾角為90°，滑動(dòng)角為0°；斷層埋深5km；斷層走向方向長(zhǎng)6km，傾向方向?qū)?km；滑動(dòng)量均為0.75 m，震源時(shí)間函數(shù)為持續(xù)時(shí)間1 s的箱型函數(shù).介質(zhì)模型選用兩層半介質(zhì)模型，具體參數(shù)如表1所示.另外，采用的頻率范圍為0—0.5Hz，慢度區(qū)間0.001—1 s/km，采樣間隔為0.25 s，波形時(shí)間長(zhǎng)度為128 s.我們計(jì)算了如圖6和表2所示的6個(gè)觀測(cè)點(diǎn)的地面運(yùn)動(dòng).

圖6 有限斷層模型與觀測(cè)臺(tái)站分布Fig.6 The model of finite source and observation stations

表1 介質(zhì)模型參數(shù)Table 1 Parameters of the stratified media

表2 6個(gè)觀測(cè)臺(tái)站的N－E坐標(biāo)Table 2 Coordinates of observation stations

為了同時(shí)測(cè)試程序的可靠性和計(jì)算效率，如圖7所示，我們?cè)O(shè)置了4種情形，分別利用1個(gè)點(diǎn)源、3個(gè)點(diǎn)源、6個(gè)點(diǎn)源和9個(gè)點(diǎn)源來(lái)表征有限斷層模型.為了定量比較兩組結(jié)果之間的一致性或者差異，計(jì)算和評(píng)價(jià)了對(duì)應(yīng)波形的相關(guān)系數(shù)和最大振幅差.最大振幅差由E＝｜P－P0｜/P0×100%計(jì)算得到，這里P0表示CompSyn計(jì)算結(jié)果的最大振幅，P表示GRTMatSyn計(jì)算結(jié)果的最大振幅.另外，為了節(jié)省篇幅，在此我們只展示南北分向的波形比較.更多的測(cè)試參見(jiàn)邸海濱（2011）文章.

圖7 利用不同數(shù)量點(diǎn)源表征的有限震源模型（a）1個(gè)點(diǎn)源；（b）3個(gè)點(diǎn)源；（c）6個(gè)點(diǎn)源；（d）9個(gè)點(diǎn)源Fig.7 Four cases of finite source approximated with point－sources（a）1 point－source；（b）3 point－sources；（c）6 point－sources；（d）9 point－sources

圖8 4種情況下CompSyn（粗線）與本文程序（細(xì)線）計(jì)算南北向波形的比較（a）采用1個(gè)點(diǎn)源表示有限震源模型；（b）采用3個(gè)點(diǎn)源表示有限震源模型；（c）采用6個(gè)點(diǎn)源表示有限震源模型；（d）采用9個(gè)點(diǎn)源表示有限震源模型Fig.8 Comparison of the waveforms computed with CompSyn（bold line）and GRTMatSyn（light line）in four cases（a）With 1 point－source；（b）With 3 point－sources；（c）With 6 point－sources；（d）With 9 point－sources

圖8展示了4種情況下GRTMatSyn與CompSyn計(jì)算結(jié)果的比較.在每個(gè)子圖中，粗線表示CompSyn計(jì)算的地面運(yùn)動(dòng)，細(xì)線表示GRTMatSyn計(jì)算的地面運(yùn)動(dòng).子圖左邊的數(shù)字從上而下依次為：CompSyn計(jì)算結(jié)果的最大振幅、兩個(gè)結(jié)果之間的相關(guān)系數(shù)和GRTMatSyn計(jì)算結(jié)果的最大振幅；右邊的字符自上而下依次為觀測(cè)臺(tái)站編號(hào)與南北向分量；最下方給出波形對(duì)應(yīng)的時(shí)間刻度.從圖8可以看出，僅用一個(gè)點(diǎn)源代替有限斷層時(shí)，距震源較近的觀測(cè)點(diǎn)的波形的一致性較差，相關(guān)系數(shù)較小，最大振幅差較大.例如，最近的觀測(cè)點(diǎn)（No.1）的相關(guān)系數(shù)僅為0.72，最大振幅差達(dá)40%.隨著震中距的增大，二者之間的相關(guān)系數(shù)增大，最大振幅差減小.例如，最遠(yuǎn)的觀測(cè)點(diǎn)（No.6）的相關(guān)系數(shù)達(dá)到0.92，振幅差減小至24.6%.不過(guò)，隨著點(diǎn)源數(shù)目的增加（從1到9），兩組結(jié)果之間的差異減小，即使在最近的觀測(cè)點(diǎn)，兩個(gè)波形之間的相關(guān)系數(shù)也可以達(dá)到0.96，最大振幅差不超過(guò)10%.我們把點(diǎn)源數(shù)目增加到更多（例如15），但相關(guān)系數(shù)和最大振幅差也基本不變.我們認(rèn)為，兩組結(jié)果之間的這種差異很可能與采用的數(shù)值算法有關(guān).比如，在計(jì)算貝塞爾函數(shù)時(shí)，GRTMat－Syn采用精確計(jì)算的方法，而CompSyn采用查表并插值的方法.

為了解程序的計(jì)算效率，我們統(tǒng)計(jì)了相同情況下，GRTMatSyn與CompSyn的計(jì)算時(shí)間.我們使用的計(jì)算機(jī)主要參數(shù)為四核CPU、4G內(nèi)存、500G硬盤(pán)、512M獨(dú)立顯卡.表3僅給出了3個(gè)點(diǎn)源情況下兩個(gè)程序的計(jì)算耗時(shí).不難看出，計(jì)算同一觀測(cè)點(diǎn)的地面運(yùn)動(dòng)，GRTMatSyn所用時(shí)間較少，而CompSyn所用的時(shí)間較多；計(jì)算不同的觀測(cè)點(diǎn)的地面運(yùn)動(dòng)，GRTMatSyn所用時(shí)間不隨震中距的變化而變化，而CompSyn所用的時(shí)間隨著震中距的增大而增大.

表3 本文程序（GRTMatSyn）與CompSyn計(jì)算耗時(shí)的比較Table 3 Comparison of the computing time between our program and CompSyn

4 討論與結(jié)論

在程序的適用性方面，GRTMatSyn的設(shè)計(jì)考慮了任意幾何特征的有限震源情況以及有限震源的震源機(jī)制隨時(shí)間和空間發(fā)生變化的情況.考慮到大地震斷層可能不是一個(gè)平面，所以，GRTMatSyn的設(shè)計(jì)沒(méi)有局限于平面斷層；考慮到震源機(jī)制在震源過(guò)程中可能會(huì)發(fā)生變化，所以，GRTMatSyn從底層設(shè)計(jì)時(shí)就考慮了任意震源機(jī)制的需求.即便是一組和多組爆炸源或其它類(lèi)型的震源引起的地面運(yùn)動(dòng)，GRTMatSyn也能勝任.

在程序的計(jì)算效率方面，GRTMatSyn的設(shè)計(jì)考慮了計(jì)算流程的優(yōu)化和并行計(jì)算.計(jì)算流程的優(yōu)化主要出現(xiàn)在通過(guò)一次循環(huán)計(jì)算所有不同深度的震源對(duì)應(yīng)的廣義反透射系數(shù)矩陣的環(huán)節(jié)上，避開(kāi)早期的針對(duì)不同震源深度的多次循環(huán).并行計(jì)算主要出現(xiàn)在兩個(gè)環(huán)節(jié)：第一個(gè)環(huán)節(jié)是廣義反透射系數(shù)矩陣的計(jì)算.廣義反透射系數(shù)矩陣的計(jì)算分上行波和下行波兩種，二者并不發(fā)生聯(lián)系，因此，設(shè)計(jì)了并行計(jì)算；第二個(gè)環(huán)節(jié)是有限震源的地面運(yùn)動(dòng)計(jì)算.在這個(gè)環(huán)節(jié)上，把觀測(cè)點(diǎn)分為多組進(jìn)行并行計(jì)算.

由于對(duì)適用性和計(jì)算效率的考慮，GRTMatSyn不但可以計(jì)算任意有限震源引起的地面運(yùn)動(dòng)，同樣也能夠計(jì)算點(diǎn)源引起的地面運(yùn)動(dòng).GRTMatSyn的計(jì)算結(jié)果與CompSyn的計(jì)算結(jié)果的對(duì)比分析表明，GRTMatSyn的計(jì)算結(jié)果是可靠的，計(jì)算效率也得到了提高.

除計(jì)算點(diǎn)源和有限震源引起的地面運(yùn)動(dòng)外，GRTMatSyn在基礎(chǔ)研究方面也將發(fā)揮積極的作用.例如，介質(zhì)結(jié)構(gòu)的波形反演問(wèn)題.我們知道，介質(zhì)結(jié)構(gòu)的波形反演問(wèn)題是一個(gè)非線性反演問(wèn)題，因此，反演效率直接取決于波形的計(jì)算效率.如果波形計(jì)算的效率高，介質(zhì)結(jié)構(gòu)的波形反演就可行且效率高；反之，即不可行或效率不高.又如，區(qū)域地震或地方地震的震源機(jī)制反演問(wèn)題，不同于遠(yuǎn)場(chǎng)地震震源機(jī)制的波形反演.區(qū)域地震或地方地震的震源機(jī)制反演，由于介質(zhì)響應(yīng)的復(fù)雜性和震源深度對(duì)波形的影響，也通常采用非線性反演，需要反復(fù)計(jì)算不同震源深度、基本震源機(jī)制的震源產(chǎn)生的波形.因此，反演效率也直接取決于波形的計(jì)算效率.

當(dāng)然，該程序也有其局限性，僅適用于計(jì)算分層均勻介質(zhì)中震源引起的地面運(yùn)動(dòng).在計(jì)算機(jī)技術(shù)高速發(fā)展的今天，三維介質(zhì)中地震波的計(jì)算已成為可能，或?qū)⒑芸鞈?yīng)用于實(shí)際地震的應(yīng)急響應(yīng).

邸海濱.2011.層狀介質(zhì)中有限震源引起的地面運(yùn)動(dòng)計(jì)算［D］.北京：中國(guó)地震局地球物理研究所：51－－90.

李旭.1993.用地震波形資料反演1990年青海共和地震的震源過(guò)程［D］.北京：國(guó)家地震局地球物理研究所：1－－45.

李旭，陳運(yùn)泰.1996.合成地震圖的廣義反射透射系數(shù)矩陣方法［J］.地震地磁觀測(cè)與研究，17（2）：1－－20.

劉超.2011.斷層厚度的地震效應(yīng)和非對(duì)稱(chēng)矩張量［D］.北京：中國(guó)地震局地球物理研究所：1－－110.

姚振興，鄭天愉.1985.對(duì)二灘水電站壩區(qū)場(chǎng)地地面運(yùn)動(dòng)的估計(jì)［J］.地球物理學(xué)報(bào)，28（2）：170－－180.

張勇.2008.震源破裂過(guò)程反演方法研究［D］.北京：北京大學(xué)：1－－16.

張勇，馮萬(wàn)鵬，許力生，周成虎，陳運(yùn)泰.2008.2008年汶川大地震的時(shí)空破裂過(guò)程［J］.中國(guó)科學(xué)：D輯，38（10）：1186－－1194.

張勇，許力生，陳運(yùn)泰.2010.2010年青海玉樹(shù)地震震源破裂過(guò)程［J］.中國(guó)科學(xué)：D輯，40（7）：819－－821.

Bouchon M.1981.A simple method to calculate Green′s functions for elastic layered media［J］.Bull Seism Soc Amer，71（4）：957－－971.

Chapman C H.1974.Generalized ray theory for an inhomogeneous medium［J］.Geophys J R astr Soc，36（3）：673－－704.

Chen X F.1999.Seismogram synthesis in multi－layered half－space partⅠ：Theoretical formulations［J］.Earthquake Research in China，13（2）：149－－174.

Gilbert F，Helmberger D V.1972.Generalized ray theory for a layered sphere［J］.Geophys J Inter，27（1）：57－－80.

Helmberger D V.1968.The crust－mantle transition in the Bering Sea［J］.Bull Seism Soc Amer，58（1）：179－－214.

Kennett B L N，Kerry N J.1979.Seismic waves in a stratified half space［J］.Geophys J R astr Soc，58：179－－214.

Kennett B L N.1980.Seismic waves in a stratified spaceⅡ：Theoretical seismograms［J］.Geophys J R astr Soc，61（1）：1－－10.

Kennett B L N.1983.Seismic Wave Propagation in a Stratified Media［M］.Cambridge：Cambridge University Press：1－－342.

Olson A H，Orcutt J A，F(xiàn)razier G A.1984.The discrete wave number/finite element method for synthetic seismograms［J］.Geophys J R astr Soc，77（2）：421－－460.

Spudich P，Archuleta R.1987.Techniques for earthquake ground motion calculation with applications to source parameterization of finite faults［G］∥Bolt B A ed.Seismic Strong Motion Synthetics.Orlando，F(xiàn)lorida：Academic Press：205－－265.

Spudich P，Xu L S.2003.Software for calculating earthquake ground motion from finite faults in vertically varying media［G］∥Lee W H K，Kanamori H，Jennings P C，Kisslinger C eds.International Handbook of Earthquake and Engineering Seismology：1633－－1634.

Yao Z X，Harkrider D G.1983.A generalized reflection transmission coefficient matrix and discrete wave number method for synthetic seismograms［J］.Bull Seism Soc Amer，73（6A）：1685－－1699.

附錄A 自由表面的反射系數(shù)和接收系數(shù)矩陣

自由表面的反射系數(shù)矩陣～R和～W接收系數(shù)矩陣表達(dá)式（Kennett，1983）為對(duì)于P－SV波

式中，i表示虛數(shù)單位，p表示慢度，α0與β0分別表示自由表面的P波與S波波速且表示自由表面處的介質(zhì)密度.

附錄B 反射透射系數(shù)矩陣的計(jì)算方法

已知地震波在傳播過(guò)程中，遇到地層分界面會(huì)產(chǎn)生反射與透射作用，反射與透射系數(shù)矩陣被用于描述這一作用.假設(shè)在深度z處存在一地層分界面，上下層面分別以z－和z＋表示.根據(jù)廣義反射透射系數(shù)矩陣方法的基本理論（Kennett，1983），以RD（z－，z＋）和TD（z－，z＋）表示下行波的反射系數(shù)和透射系數(shù)矩陣，RU（z－，z＋）和TU（z－，z＋）表示上行波的反射系數(shù)和透射系數(shù)矩陣，則

已知地層分界面處的傳遞矩陣Q（z－，z＋）具有如下表達(dá)式（Kennett，1983）：

式中，上標(biāo)“T”表示矩陣轉(zhuǎn)置，下標(biāo)“－”和“＋”表示分界面上下兩層介質(zhì).SV

式中，i表示虛數(shù)單位，α與β分別表示P波與S波的波速，p表示慢度，ρ表示介質(zhì)密度，且

將公式（B－3）代入（B－1）和（B－2），即可完成地層分界面處反射透射系數(shù)矩陣的計(jì)算.

Calculation of the ground motion generated by a finite source in stratified media：An extension of the point－source case

Di Haibin Xu Lisheng
（Institute of Geophysics，China Earthquake Administration，Beijing 100081，China）

10.3969/j.issn.0253－3782.2012.04.001

P315.3＋3

A

中國(guó)地震局地球物理研究所基本業(yè)務(wù)費(fèi)（DQJB11C12）和國(guó)家自然科學(xué)基金（40874026）資助.

2011－06－30收到初稿，2011－07－12決定采用修改稿.

e－mail：xuls＠cea－igp.ac.cn

邸海濱，許力生.2012.層狀介質(zhì)中有限震源引起的地面運(yùn)動(dòng)的計(jì)算對(duì)點(diǎn)源情形的拓展.地震學(xué)報(bào)，34（4）：425－－438.

Di Haibin，Xu Lisheng.2012.Calculation of the ground motion generated by a finite source in stratified media：An extension of the point－source case.Acta Seismologica Sinica，34（4）：425－－438.