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

邸海濱 許力生

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

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

邸海濱 許力生

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

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

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

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

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

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

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

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

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

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

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

圖1 分層均勻介質(zhì)及坐標(biāo)構(gòu)架（Kennett，1983）z0表示自由表面，zL表示介質(zhì)底界面，zS表示震源處的虛界面.U和D分別表示上行波場和下行波場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)動的面諧矢量位移解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依賴于自由表面的接收系數(shù)矩陣、反射系數(shù)矩陣、分層介質(zhì)中上行波和下行波的廣義反透射系數(shù)，以及上行波與下行波地震波矢量間斷.如果通過這些量得到面諧矢量位移解w0，那么很容易通過Fourier－Hankel反變換得到時空域點(diǎn)源引起的地面運(yùn)動解.可見，計(jì)算地面運(yùn)動的關(guān)鍵在于計(jì)算面諧矢量位移解w0，而計(jì)算面諧矢量位移解w0的關(guān)鍵則在于計(jì)算廣義反透射系數(shù)矩陣和表征震源作用的地震波矢量間斷.因此，提高計(jì)算有限震源引起的地面運(yùn)動效率的途徑之一是優(yōu)化廣義反透射系數(shù)矩陣和表示震源作用的地震波矢量間斷的計(jì)算流程.

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

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

圖3 地震波矢量在分層介質(zhì)中傳播的幾何解釋以水平實(shí)線表示文中考慮的3個分界面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）來描述（Kennett，1983），下列各式把描述一般點(diǎn)源的地震矩張量和地震波矢量間斷聯(lián)系起來.

對于P－SV波

對于SH波

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

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

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

式中

其中，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)動的計(jì)算

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

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

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

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

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

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

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

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

同理

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

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

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

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

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

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

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

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

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

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

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

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

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

3 可靠性與效率的測試

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

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

圖6 有限斷層模型與觀測臺站分布Fig.6 The model of finite source and observation stations

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

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

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

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

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

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

4 討論與結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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

張勇，許力生，陳運(yùn)泰.2010.2010年青海玉樹地震震源破裂過程［J］.中國科學(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）為對于P－SV波

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

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

已知地震波在傳播過程中，遇到地層分界面會產(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

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

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

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

邸海濱，許力生.2012.層狀介質(zhì)中有限震源引起的地面運(yùn)動的計(jì)算對點(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.