郝金來，姚振興

中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

均勻彈性分層介質(zhì)模型中的同震位移、應(yīng)變以及應(yīng)力

郝金來，姚振興

中國科學(xué)院地質(zhì)與地球物理研究所，北京 100029

地震發(fā)生時產(chǎn)生的同震位移、應(yīng)變以及應(yīng)力變化，特別是同震應(yīng)力的變化，在地震觸發(fā)等問題的研究中有著重要的意義.本文進一步發(fā)展了基于均勻彈性水平層狀介質(zhì)，利用廣義反射透射系數(shù)矩陣和離散波數(shù)計算同震位移的方法，使之可以計算相應(yīng)的應(yīng)變、應(yīng)力以及同震庫侖應(yīng)力變化.可適用于多種情況，接收點可以位于地表以及地表以下，震源類型可以是剪切位錯源以及拉張位錯源.通過與半無限介質(zhì)的解析解相比較，結(jié)果一致，驗證了方法的可靠性以及計算精度，可以用于計算地震之后庫侖應(yīng)力變化，為判斷余震分布提供參考.在計算同震位移時，使用了梯形積分與Filon積分相結(jié)合的積分方式，即提高了同震位移計算的速度，又保證了計算精度，有利于反演問題研究.

廣義反射透射系數(shù)矩陣，離散波數(shù)，同震位移，庫侖應(yīng)力，均勻分層彈性介質(zhì)

1 引 言

在地震學(xué)研究中地震發(fā)生時產(chǎn)生的同震位移、應(yīng)變以及應(yīng)力變化有著非常重要的研究意義，同震位移可用于反演斷層破裂過程，應(yīng)力變化可用于地震成因以及地震觸發(fā)等研究［1－3］.地震發(fā)生后，周邊某些斷層上庫侖應(yīng)力增大，加大了斷層破裂的幾率，相反則減少了斷層上的負載［4－7］.在許多研究區(qū)域都找到了庫侖應(yīng)力與余震分布、地震觸發(fā)的對應(yīng)關(guān)系，例如1999年土耳其伊茲米特地震與杜賽地震等［8］.因此，計算地震發(fā)生后的同震位移、應(yīng)變、應(yīng)力變化是一項重要的基礎(chǔ)工作.

自從Steketee［9］首次將靜力學(xué)位錯理論引入地震學(xué)后，有關(guān)計算靜態(tài)位移場的理論研究迅速發(fā)展起來.根據(jù)介質(zhì)模型的不同，主要有以下幾類，首先是基于彈性半無限介質(zhì)模型的計算方法［10－11］，Okada［12－13］詳細地評述了有關(guān)的研究工作，給出了半無限介質(zhì)內(nèi)位錯點源和有限源靜態(tài)位移場的解析解.這種方法的優(yōu)點是給出的結(jié)果為解析解，計算速度快，缺點是模型過于簡單，特別是當?shù)乇泶嬖诘退賹訒r，半無限介質(zhì)模型給出的結(jié)果往往會存在較大誤差［14］.因此發(fā)展起來了基于均勻彈性層狀介質(zhì)的相關(guān)理論，Sato和Matsu′ura［14］，Ben－Menaham和Singh［15］，Singh［16］，Wang等［17］給出了水平層介質(zhì)中的計算方法.Israel等［18］，Pollitz［19］給出了球體分層模型中的計算方法.層狀模型可以較好地模擬地球介質(zhì).但模擬震后較長時間間隔內(nèi)的位移、應(yīng)變、應(yīng)力變化時，黏彈性松弛效應(yīng)的影響則不可忽略［20－22］，彈性模型無法考慮這一點，而黏彈性層狀介質(zhì)模型可以更好地擬合這種震后形變.隨著計算能力的增強，利用有限元等方法，基于更復(fù)雜的橫向不均勻介質(zhì)模型計算應(yīng)力變化，在近些年逐漸得到更廣泛地應(yīng)用［23－24］.這些不同的模型有其自身的適用范圍，可以根據(jù)不同的研究目的選擇合適方法加以利用.本文主要考慮均勻彈性水平層狀介質(zhì)模型下，靜力學(xué)位移、應(yīng)變、應(yīng)力的計算.

在計算理論地震圖時，Kennett［25－26］發(fā)展了反射和透射系數(shù)矩陣方法，Bouchon［27］發(fā)展了一種可以有效地計算波數(shù)積分的離散波數(shù)方法，Yao和Harkrider［28］將兩者相結(jié)合，利用廣義反射透射系數(shù)矩陣和離散波數(shù)方法計算近場理論地震圖.謝小碧和姚振興［29］把這種處理波動問題方法加以推廣，用以處理靜態(tài)問題，給出了均勻彈性分層介質(zhì)中剪切位錯點源在地面產(chǎn)生的靜態(tài)位移場的解答，He等［30］給出了拉張源所產(chǎn)生的靜態(tài)位移場的解答.但是他們均未給出相應(yīng)的應(yīng)變、應(yīng)力的計算方法與結(jié)果.在他們工作的基礎(chǔ)上，我們進一步發(fā)展該方法，使之可以基于均勻分層彈性介質(zhì)模型，根據(jù)廣義反射透射系數(shù)矩陣和離散波數(shù)方法，計算同震應(yīng)變、應(yīng)力變化等.

謝小碧和姚振興［29］，He等［30］計算同震位移時，均采用了簡單的梯形積分方式.積分間隔的選取將直接影響到計算精度與速度，較小的積分間隔才能得到較為精確的數(shù)值解，但是這樣會增加計算時間，對于利用同震位移資料進行反演問題研究十分不利.理論地震圖計算中廣泛采用的Filon積分方法［31－33］，可以使得在選取較大的積分間隔時，計算精度仍然很高.在本文中，計算同震位移時，使用梯形積分與Filon積分相結(jié)合的方法，即提高了計算速度，又使得計算結(jié)果具有較高精度.

在本文主要考慮的均勻彈性水平層狀介質(zhì)模型中，前人研究主要存在以下一些問題：某些研究中震源類型并不全面，特別是拉張源的計算；某些研究中僅給出了位移的計算結(jié)果；某些研究中水平層的數(shù)量受到制約等［17］.Wang等［17］給出了一種全面的計算方法，使用了Thomson－Haskell矩陣進行計算，矩陣元素中存在導(dǎo)致溢出問題的指數(shù)項.Wang等［17］通過一些方法，如正交歸一化等避免這些問題的出現(xiàn).而廣義反射透射系數(shù)矩陣摒棄了造成數(shù)值不穩(wěn)定的指數(shù)項，解決了數(shù)值溢出等問題，因此我們利用廣義反射透射系數(shù)矩陣方法計算同震位移、應(yīng)變、應(yīng)力變化.

2 方 法

2.1 同震應(yīng)變、應(yīng)力以及庫侖應(yīng)力的計算

圖1 水平層狀介質(zhì)模型水平層狀介質(zhì)模型，1，2，3等表示水平層狀介質(zhì)，Z1，Z2，Z3等表示界面，r表示震中距，五角星標示了震源的位置，即Zs層.Fig.1 Layered medium modelLayered medium model，1，2，3represent the layered medium，Z1，Z2，Z3represent the interface，r represents the epicenter distance，and the star shows the location of the source，in the interface of Zs.

根據(jù)先前的研究結(jié)果［29－30］，在如圖1所示的水平層狀介質(zhì)中，三分量的靜力學(xué)位移解答可以表示成如下形式：其中wz，qr，vθ分別表示垂向、徑向和切向的靜力學(xué)位移，Σ為斷層面元的面積，Δui，i＝1，2，3分別對應(yīng)走滑、傾滑和拉張斷層的靜力學(xué)位錯的錯距，Aim為方向因子，qm，wm，vm分別為面諧矢量坐標系下的位移分量，Jm（kr），J′m（kr）分別為m階貝塞爾函數(shù)及其導(dǎo)數(shù).

根據(jù)位移解答求應(yīng)變與應(yīng)力的關(guān)鍵步驟為求位移分量的偏導(dǎo)數(shù)，即求解.公式（1）中與θ有關(guān)的量僅為方向因子Aim，與r有關(guān)的量僅為貝塞爾函數(shù)及其導(dǎo)數(shù)與z有關(guān)的量僅為qm，wm，vm.

對θ的偏導(dǎo)數(shù)表示為

根據(jù)方向性因子的表達式可得到相應(yīng)的偏導(dǎo)數(shù)：

對r的偏導(dǎo)數(shù)表示為

對z的偏導(dǎo)數(shù)表示為

當接收點位于地面或介質(zhì)內(nèi)部時，qmwm和vm表達式各不相同.當接收點位于地表時有

震源系數(shù)分別為

此處Pm，SVm和SHm分別表示P＋m，P－m，SV＋m，SV－m和SH＋m，SH－m，而

當接收點位于介質(zhì)內(nèi)部時，分接收點的深度zR小于震源深度和大于震源深度兩種情況.當接收點位于震源以上時，即當zR＜zs時，有

其中接收函數(shù)矩陣和接收函數(shù)分別為

當接收點位于震源以下時，即當zR＞zs，有

其中接收函數(shù)矩陣和接收函數(shù)分別為

根據(jù)位移的偏導(dǎo)數(shù)，可得到應(yīng)變、應(yīng)力解答分別為［34－35］

其中err，eθθ，ezz，erθ，erz，eθz為應(yīng)變分量，τrr，τθθ，τzz，τrθ，τrz，τθz為應(yīng)力分量，urr，uθθ，uzz，urθ，urz，uθz為位移的偏導(dǎo)數(shù)，λ，μ為彈性常數(shù).根據(jù)應(yīng)力張量可以計算斷層面上的庫侖應(yīng)力變化［4，36－37］，如下式所示：

其中Δσ，Δτs，Δσn，μ′分別表示庫侖應(yīng)力，斷層面上的剪應(yīng)力，斷層面上的正應(yīng)力，以及等效摩擦系數(shù).

2.2 同震位移計算中積分方式的改進

同震位移計算中的關(guān)鍵為計算以下積分：

為了提高計算速度，需要選取較大的積分間隔，如果此時仍使用梯形積分，將會帶來較大的計算誤差.當kr較大時，Jm（kr）存在漸進展開式，可以采用Filon積分方式，選取較大的積分間隔，仍會有較高的計算精度［31－33］.因此在積分的計算過程中，采用分段積分的方式.當kr較小時，選取較小的積分間隔，采用梯形積分.當kr較大時，選取較大的積分間隔，采用紀晨和姚振興改進的Filon積分方法［32］，僅取貝塞爾函數(shù)Jm（kr）漸進展開式的零次項

這樣既保證了計算精度，又提高了計算速度.

3 數(shù)值檢驗

為了檢驗公式的正確性及數(shù)值方法的精度，我們做了一系列的數(shù)值檢驗.分別驗證均勻彈性半無限介質(zhì)模型中接收點位于不同深度時，以及分層介質(zhì)模型中接收點位于地表時，點源產(chǎn)生的位移以及應(yīng)力解的正確性.在笛卡爾坐標系中，設(shè)定沿著斷層走向為X軸正方向，逆時針旋轉(zhuǎn)90°為Y軸正方向，垂直向上為Z軸正方向，斷層左上角對應(yīng)于地表的投影點為坐標原點，如圖2所示.

圖2 坐標系示意圖Fig.2 The coordinate system

3.1 接收點位于地表時的點源檢驗

均勻彈性半無限介質(zhì)參數(shù)設(shè)定為P波速度4.0km／s，S波速度2.7km／s，密度2.5g／cm3.使用10層，每層厚度為3km，具有相同P波速度，S波速度以及密度的均勻彈性介質(zhì)層，之下設(shè)置為半無限，用來模擬均勻彈性半空間.震源參數(shù)設(shè)定如下：傾角為70°，斷層面積0.01km×0.01km，震源深度5.0km.接收點參數(shù)設(shè)定如下：接收點位于地表，在笛卡爾坐標系中，接收點橫坐標3.0km，縱坐標從－10.0km起始，到10.0km截止，共取41個點.分別計算了走滑斷層、傾滑斷層以及拉張斷層所產(chǎn)生的位移以及應(yīng)力結(jié)果，并將之與Okada給出的解析解［12－13］進行了對比.圖3a，圖3b分別給出了位移，應(yīng)力結(jié)果的對比圖，由于接收點位于地表，相應(yīng)的應(yīng)力分量只有徑向、切向正應(yīng)力以及徑向與切向之間的剪應(yīng)力.與解析解相比較，兩者一致.

3.2 接收點位于地下，深度小于震源深度時的點源檢驗

采用與3.1相同的介質(zhì)設(shè)定，震源參數(shù)設(shè)定如下：傾角為70°，斷層面積0.01km×0.01km，震源深度20.0km，接收點設(shè)定如下：接收點位于地下，深度16.0km，接收點橫坐標3.0km，縱坐標從－10.0km起始，到10.0km截止，共取41個點.圖4a，圖4b分別給出了位移，應(yīng)力結(jié)果與Okada的解析解［12－13］的對比圖，兩者一致.

3.3 接收點位于地下，深度大于震源深度時的點源檢驗

采用與3.1節(jié)相同的介質(zhì)設(shè)定，震源參數(shù)設(shè)定如下：傾角為70°，斷層面積0.01km×0.01km，震源深度10.0km，接收點設(shè)定如下：接收點位于地下，深度14.0km，接收點橫坐標3.0km，縱坐標從－10.0 km起始，到10.0km截止，共取41個點.圖5a，圖5b分別給出了位移，應(yīng)力結(jié)果與Okada的解析解［12－13］的對比圖，兩者一致.通過以上數(shù)值檢驗可以看出，基于均勻彈性分層介質(zhì)模型，根據(jù)廣義反射透射矩陣和離散波數(shù)方法計算出的位移與應(yīng)力結(jié)果與解析解一致，驗證了公式的正確性以及計算的精度.

3.4 分層介質(zhì)中接收點位于地表的點源檢驗

選取如表1所示的均勻彈性水平層狀介質(zhì)模型，將接收點設(shè)置于地表，計算點源產(chǎn)生的位移以及應(yīng)力解，通過與Wang等（2003）的計算結(jié)果比較，驗證在分層介質(zhì)中的數(shù)值計算精度.震源參數(shù)設(shè)定如下：傾角為70°，斷層面積0.01km×0.01km，震源深度5.0km.接收點參數(shù)設(shè)定如下：接收點位于地表，在笛卡爾坐標系中，接收點橫坐標3.0km，縱坐標從－10.0km起始，到10.0km截止，共取41個點.計算了走滑斷層與傾滑斷層所產(chǎn)生的位移以及應(yīng)力結(jié)果，圖6給出了位移，應(yīng)力結(jié)果的對比圖，由于接收點位于地表，相應(yīng)的應(yīng)力分量只有徑向、切向正應(yīng)力以及徑向與切向之間的剪應(yīng)力.由圖可見，兩者一致.

表1 分層介質(zhì)模型Table 1 Layered velocity model

4 庫侖應(yīng)力的計算

根據(jù)計算出的同震應(yīng)力張量，即可得出斷層面上的靜態(tài)庫侖應(yīng)力變化，用于地震觸發(fā)研究［4］.考慮一個斷層，走向0°，傾角70°，滑動量1m，滑動角45°，斷層長度20km，斷層寬度10km，計算庫侖應(yīng)力需要有斷層面以及滑動方向，假設(shè)接收點處的斷層面以及滑動方向與破裂斷層一致，計算深度為地下5.0km，接收點橫坐標從－20.0km起始，到40.0km截止，共取13個點，縱坐標從20.0km起始，到40.0km截止，共取5個點，首先基于半無限介質(zhì)模型，介質(zhì)參數(shù)為P波速度6.2km／s，S波速度3.5km／s，密度2.9g／cm3，分別根據(jù)Okada［12－13］和本文的方法計算庫侖應(yīng)力變化.圖7a分別給出了相應(yīng)的計算結(jié)果，以及絕對誤差，由圖可見兩者誤差很小，再次驗證了數(shù)值計算的精度.考慮分層介質(zhì)模型（表2），計算得到庫侖應(yīng)力分布，圖7b給出其與半無限結(jié)果的比較圖，兩者在庫侖應(yīng)力大小上存在著顯著差異，因此在計算應(yīng)力觸發(fā)時，考慮分層彈性半無限介質(zhì)模型是有必要的.

表2 速度結(jié)構(gòu)模型Table 2 The velocity model

圖3 均勻彈性半無限介質(zhì)中，接收點位于地表時，位移、應(yīng)力結(jié)果比較（a）位移結(jié)果比較，橫坐標為接收點位置，縱坐標為位移大小，左、中、右分別對應(yīng)于X、Y和Z方向位移，上、中、下分別對應(yīng)于走滑、傾滑和拉張斷層，其中方框為Okada給出的解析解結(jié)果，圓點為本文程序給出的計算結(jié)果；（b）應(yīng)力結(jié)果比較，橫坐標為接收點位置，縱坐標為應(yīng)力張量各分量大小，左、中、右分別對應(yīng)于走滑、傾滑和拉張斷層，上對應(yīng)于正應(yīng)力，其中曲線為Okada給出的解析解結(jié)果，圓點為本文程序給出的計算結(jié)果.黑色為X向正應(yīng)力，灰色為Y向正應(yīng)力，下為兩者之間的剪切應(yīng)力.Fig.3 The displacement and stress in the half－space model，the receiver is in the surface（a）The comparison of displacement results，the X－axis shows the location of the receiver，Y－axis shows the displacement，the left，middle，right columns represent the displacements of X，Y，Zcomponent respectively，the upper，middle，lower lines represent the strike－slip，dipslip，tensile sources respectively，the squares show the analytical solutions，the points show the results of our method；（b）The comparison of stress results，the X－axis shows the location of the receiver，Y－axis shows the stress，the left，middle，right columns represent the strike－slip，dip－slip，tensile sources respectively，the upper lines show the analytical solutions，the points show the results of our method，black，gray represent the normal stress ofX，Ycomponent respectively and the lower lines represent the shear stress.

圖4 均勻彈性半無限介質(zhì)中，接收點位于地下，深度小于震源深度時，位移、應(yīng)力結(jié)果比較（a）位移結(jié)果比較，橫坐標為接收點位置，縱坐標為位移大小，左、中、右分別對應(yīng)于X、Y和Z方向位移，上、中、下分別對應(yīng)于走滑、傾滑和拉張斷層，其中紅色方框為Okada給出的解析解結(jié)果，藍色圓點為本文程序給出的計算結(jié)果；（b）應(yīng)力結(jié)果比較，橫坐標為接收點位置，縱坐標為應(yīng)力張量各分量大小，左、中、右分別對應(yīng)于走滑、傾滑和拉張斷層，其中直線為Okada給出的解析解結(jié)果，圓點為本文程序給出的計算結(jié)果.上半部分為正應(yīng)力，黑、紅、藍色分別表示X、Y和Z方向的正應(yīng)力，下半部分為剪切應(yīng)力，黑、紅、藍色分別表示XY、XZ和YZ方向的剪切應(yīng)力.Fig.4 The displacement and stress in the half－space model，the receiver is underground and the depth of it is above that of the source（a）The comparison of displacement results，the X－axis shows the location of the receiver，Y－axis shows the displacement，the left，middle，right columns represent the displacements of X，Y，Zcomponent respectively，the upper，middle，lower lines represent the strike－slip，dip－slip，tensile sources respectively，the red squares show the analytical solutions，the blue points show the results of our method；（b）The comparison of stress results，the X－axis shows the location of the receiver，Y－axis shows the stress，the left，middle，right columns represent the strike－slip，dip－slip，tensile sources respectively，the upper，lower lines represent the normal，shear stress respectively，in the normal stress，black，red，blue represent the normal stress ofX，Y，Zcomponent respectively，in the shear stress，black，red，blue represent the shear stress of XY，XZ，YZcomponent respectively.

圖5 均勻彈性半無限介質(zhì)中，接收點位于地下，深度大于震源深度時，位移、應(yīng)力結(jié)果比較（a）位移結(jié)果比較，橫坐標為接收點位置，縱坐標為位移大小，左、中、右分別對應(yīng)于X、Y和Z方向位移，上、中、下分別對應(yīng)于走滑、傾滑和拉張斷層，其中紅色方框為Okada給出的解析解結(jié)果，藍色圓點為本文程序給出的計算結(jié)果；（b）應(yīng)力結(jié)果比較，橫坐標為接收點位置，縱坐標為應(yīng)力張量各分量大小，左、中、右分別對應(yīng)于走滑、傾滑和拉張斷層，其中直線為Okada給出的解析解結(jié)果，圓點為本文程序給出的計算結(jié)果.上半部分為正應(yīng)力，黑、紅、藍色分別表示X、Y和Z方向的正應(yīng)力，下半部分為剪切應(yīng)力，黑、紅、藍色分別表示XY、XZ和YZ方向的剪切應(yīng)力.Fig.5 The displacement and stress in the half－space model，the receiver is underground and the depth of it is blow that of the source（a）The comparison of displacement results，the X－axis shows the location of the receiver，Y－axis shows the displacement，the left，middle，right columns represent the displacements of X，Y，Zcomponent respectively，the upper，middle，lower lines represent the strike－slip，dip－slip，tensile sources respectively，the red squares show the analytical solutions，the blue points show the results of our method；（b）The comparison of stress results，the X－axis shows the location of the receiver，Y－axis shows the stress，the left，middle，right columns represent the strike－slip，dip－slip，tensile sources respectively，the upper，lower lines represent the normal，shear stress respectively，in the normal stress，black，red，blue represent the normal stress of X，Y，Zcomponent respectively，in the shear stress，black，red，blue represent the shear stress of XY，XZ，YZcomponent respectively.

圖6 均勻彈性水平層狀介質(zhì)中，接收點位于地表時的位移、應(yīng)力結(jié)果比較（a）位移結(jié)果比較，橫坐標為接收點位置，縱坐標為位移大小，左、中、右分別對應(yīng)于X、Y和Z方向位移，上、下分別對應(yīng)于走滑、傾滑斷層，其中方框為Wang et al.（2003）［17］給出的結(jié)果，圓點為本文程序給出的計算結(jié)果；（b）應(yīng)力結(jié)果比較，橫坐標為接收點位置，縱坐標為應(yīng)力張量各分量大小，左、右分別對應(yīng)于走滑、傾滑斷層，上對應(yīng)于正應(yīng)力，其中曲線為Wang et al.（2003）［17］給出的結(jié)果，圓點為本文程序給出的計算結(jié)果.黑色為X向正應(yīng)力，灰色為Y向正應(yīng)力，下為兩者之間的剪切應(yīng)力.Fig.6 The displacement and stress in the layered velocity model，the receiver is in the surface（a）The comparison of displacement results，the X－axis shows the location of the receiver，Y－axis shows the displacement，the left，middle，right columns represent the displacements of X，Y，Zcomponent respectively，the upper，lower lines represent the strike－slip，dip－slip sources respectively，the squares show the solutions of Wang et al.（2003）［17］，the points show the results of our method；（b）The comparison of stress results，the X－axis shows the location of the receiver，Y－axis shows the stress，the left，right columns represent the strike－slip，dip－slip sources respectively，the upper lines represent the normal stress，the lines show the solutions of Wang et al.（2003）［17］，the points show the results of our method，black，gray represent the normal stress of X，Ycomponent respectively and the lower represent the shear stress the shear stress.

圖7 庫侖應(yīng)力比較（a）半無限介質(zhì)解析解與本文計算結(jié)果比較，左、中、右分別為解析解、本文結(jié)果以及誤差絕對值；（b）半無限介質(zhì)與分層介質(zhì)結(jié)果比較，左、中、右分別為半無限介質(zhì)結(jié)果、分層介質(zhì)結(jié)果以及誤差絕對值.Fig.7 The Coulomb stress（a）The results of analytical and our method in the half－space model，the left，middle，right columns represent the analytical solution，result of our method and absolute value of discrepancy between them respectively；（b）The comparison of the half－space model and layered model，the left，middle，right columns represent the solution in the half－space，result in the layered model and absolute value of discrepancy between them respectively.

表3 不同積分方式的計算結(jié)果比較Table 3 Comparison of different integration methods

5 同震位移計算中積分方式的比較

為了提高同震位移的計算速度，利于反演問題研究，我們改進了積分方式.通過不同的積分方法計算了半無限介質(zhì)模型中點源產(chǎn)生的同震位移，考察了改進方式的優(yōu)越性.接收點設(shè)定于地表，坐標為（20km，250km）.震源參數(shù)設(shè)定為：震源深度7.0km，斷層面積0.01km×0.01km，傾角為60°.

由計算結(jié)果（表3）可見，采用改進的積分方式（積分方式二）計算同震位移時，核函數(shù)計算次數(shù)減少了2／3左右，而誤差均控在3%以內(nèi).這樣既保證了計算精度，又提高了計算速度，有利于反演問題研究.

6 結(jié) 論

本文進一步發(fā)展了基于均勻彈性水平層狀介質(zhì)，利用廣義反射透射系數(shù)矩陣和離散波數(shù)計算同震位移的方法［29－30］，使之可以計算相應(yīng)的應(yīng)變、應(yīng)力以及同震庫侖應(yīng)力變化.接收點可以位于地表以及地下，可以適應(yīng)于剪切位錯源和拉張源等.根據(jù)應(yīng)力張量計算同震庫侖應(yīng)力的變化，可用于地震觸發(fā)研究，為判斷余震分布提供參考.謝小碧和姚振興［29］，He等［30］計算同震位移時，均采用了簡單的梯形積分方式，本文使用了梯形積分與Filon積分相結(jié)合的積分方式，既提高了同震位移計算的速度，又保證了計算精度，有利于反演問題研究.

（References）

［1］ Stein R S，King G C P，Lin J.Change in failure stress on the southern San Andreas fault system caused by the 1992 magnitude＝7.4Landers earthquake.Science，1992，258（5086）：1328－1332.

［2］ Bouchon M.The state of stress on some faults of the San Andreas system as inferred from near－field strong motion data.J.Geophys.Res.，1997，102（B6）：11731－11744.

［3］ 沈正康，萬永革，甘衛(wèi)軍等.東昆侖活動斷裂帶大地震之間的黏彈性應(yīng)力觸發(fā)研究.地球物理學(xué)報，2003，46（6）：786－795.Shen Z K，Wan Y G，Gan W J，et al.Viscoelastic triggering among large earthquakes along the east Kunlun fault system.Chinese J.Geophys.（in Chinese），2003，46（6）：786－795.

［4］ King G C P，Stein R S，Lin J.Static stress changes and the triggering of earthquakes.Bull.Seism.Soc.Am.，1994，84（3）：935－953.

［5］ 石耀霖，曹建玲.庫侖應(yīng)力計算及應(yīng)用過程中若干問題的討論—以漢川地震為例.地球物理學(xué)報，2010，53（1）：102－110.Shi Y L，Cao J L.Some aspects in static stress change calculation—case study on Wenchuan earthquake.Chinese J.Geophys.（in Chinese），2010，53（1）：102－110.

［6］ Toda S，Stein R S，Reasenberg P A，et al.Stress transferred by the 1995 Mw＝6.9Kobe，Japan，shock：effect on aftershocks and future earthquake probabilities.J.Geophys.Res.，1998，103：24543－24565.

［7］ 萬永革.“地震靜態(tài)應(yīng)力觸發(fā)”問題的研究.北京：中國地震局地球物理研究所，2001.Wan Y G.The study about“Static stress triggering”（in Chinese）.Beijing：Institute of Geophysics China Earthquake Administration，2001.

［8］ Barka A.The 17August 1999Izmit earthquake.Science，1999，258（5435）：1858－1859.

［9］ Steketee J A.On volterra′s dislocations in a semi－infinite elastic medium.Can.J.Phys.，1958，36（2）：192－205.

［10］ Maruyama T.Static elastic dislocation in an infinite and semiinfinite medium.Bull.Earthquake Res.Inst.Tokyo Univ.，1964，42：289－368.

［11］ Press F.Displacements，strains，and tilts at teleseismic distances.J.Geophys.Res.，1965，70（10）：2395－2412.

［12］ Okada Y.Surface deformation due to shear and tensile faults in a half－space.Bull.Seism.Soc.Am.，1985，75（4）：1135－1154.

［13］ Okada Y.Internal deformation due to shear and tensile faults in a half－space.Bull.Seism.Soc.Am.，1992，82：1018－1040.

［14］ Sato R，Matsu′ura M.Static deformations due to the fault spreading over several layers in a multilayered medium，I，displacement.J.Phys.Earth，1973，21（3）：227－249.

［15］ Ben－Menahem A，Singh S J.Multipolar elastic fields in a layered half space.Bull.Seism.Soc.Am.，1968，58（5）：1519－1572.

［16］ Singh S J.Static deformation of a multilayered half－space by internal sources.J.Geophys.Res.，1970，75（17）：3257－3263.

［17］ Wang R，Martin F L，Roth F.Computation of deformation induced by earthquakes in a multi－layered elastic crust：FORTRAN programs EDGRN／EDCMP.Comp.Geosci.，2003，29（2）：195－207.

［18］ Israel M，Ben－Menahem A，Singh S J.Residual deformation of real Earth models with application to the Chandler wobble.Geophysical Journal of the Royal Astronomical Society，1973，32（2）：219－247.

［19］ Pollitz F F.Coseismic deformation from earthquake faulting on a layered spherical earth.Geophys.J.Int.，1996，125（1）：1－14.

［20］ Nur A，Mavko G.Postseismic viscoelastic rebound.Science，1974，183（4121）：204－206.

［21］ Wang R，Martin F L，Roth F.PSGRN／PSCMP—a new code for calculating co－and post－seismic deformation，geoid and gravity changes based on the viscoelastic－gravitational dislocation theory.Comp.Geosci.，2006，32（4）：527－541.

［22］ Pollitz F F.Postseismic relaxation theory on the spherical earth.Bull.Seism.Soc.Am.，1992，82（1）：422－453.

［23］ Masterlark T，Wang H F.Transient stress coupling between the 1992Landers and 1999Hector Mine，California，earthquakes.Bull.Seism.Soc.Am.，2002，92（2）：1470－1480.

［24］ 張懷，吳忠良，張東寧等.虛擬川滇—基于千萬網(wǎng)格并行有限元計算的區(qū)域強震演化過程數(shù)值模型設(shè)計和構(gòu)建.中國科學(xué)D輯，2009，39（3）：260－270.Zhang H，Wu Z L，Zhang D N，et al.Quasi Sichuan－Yunnan－Based on thousands of grid parallel finite element computations strong earthquake evolution process numerical model design and construction.Science Chian（D）（in Chinese），2009，39（3）：260－270.

［25］ Kennett B L N.Reflections，rays，and reverberations.Bull.Seism.Soc.Am.，1974，64（6）：1685－1696.

［26］ Kennett B L N.Seismic wave propagation in stratified media.New York：Cambridge U.Press，1983.

［27］ Bouchon M.Discrete wave number representation of elastic wave fields in three－space dimensions.J.Geophys.Res.，1979，84（B7）：3609－3614.

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

［29］ 謝小碧，姚振興.計算分層介質(zhì)中位錯點源靜態(tài)位移場的廣義反射、透射系數(shù)矩陣和離散波數(shù)方法.地球物理學(xué)報，1989，32（3）：270－280.Xie X B，Yao Z X.A generalized reflection－transmition coefficient matrix method to calculate static displacement field of a stratified half－space by dislocation source.Chinese J.Geophys.（in Chinese），1989，32（3）：270－280.

［30］ He Y M，Wang W M，Yao Z X.Static deformation due to shear and tensile faults in a layered half－space.Bull.Seism.Soc.Am.，2003，93（5）：2253－2263.

［31］ Mallick S，F(xiàn)razer L N.Practical aspects of reflectivity modeling.Geophysics，1988，55（10）：1355－1364.

［32］ 紀晨，姚振興.區(qū)域地震范圍的寬頻帶理論地震圖算法研究.地球物理學(xué)報，1995，38（4）：460－468.Ji C，Yao Z X.The study of the method for broadband regional synthetic seismogram.Chinese J.Geophys.（in Chinese），1995，38（4）：460－468.

［33］ Chen X F，Zhang H M.An efficient method for computing Green′s functions for a layered half－space at large epicentral distances.Bull.Seism.Soc.Am.，2001，91（4）：858－869.

［34］ Cotton F，Coutant O.Dynamic stress variations due to shear faults in a plane－layered medium.Geophys.J.Int.，1997，128（3）：676－688.

［35］ Aki K，Richards P G.Quantitative Seismology：Theory and Methods.San Francisco：W.H.Freeman and Co.，1980.

［36］ 汪建軍.同震、震后和震間應(yīng)力觸發(fā).武漢：武漢大學(xué)，2010.Wang J J.Coseismic，Postseismic and Interseismic Stress Triggerings（in Chinese）.Wuhan：Wuhan University，2010.

［37］ 周宇明，單斌，熊熊.靜態(tài)應(yīng)力觸發(fā)中影響庫侖應(yīng)力變化的參數(shù)敏感性分析.大地測量與地球動力學(xué)，2008，28（5）：21－26.Zhou Y M，Shan B，Xiong X.Parameters sensitivity analysis of Coulomb stress change in static stress triggering.Journal of Geodesy and Geodynamics（in Chinese），2008，28（5）：21－26.

The coseismic displacement，strain and stress in the layered elastic model

HAO Jin－Lai，YAO Zhen－Xing
Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing100029，China

The coseismic displacement，strain and stress，especially the change of the coseismic stress，are significant in the earthquake triggering investigation.Based on the layered elastic model，we develop a method to calculate the coseismic dislocation by the generalized reflectiontransmission coefficient matrix and discrete wavenumber method.We can use it to calculate the coseismic strain，stress，and the change of Coulomb stress.It can be applied in the general situation.The receiver can be located on the surface or underground.The earthquake source can be pure shear dislocation source or tensile source.For the half－space model，the result of our method coincides with the analytical solution，so the precision of our method is reliable.The method can be used to calculate the Coulomb stress after an earthquake，and the change of Coulomb stress can give a reference to estimate the distribution of the aftershocks.We combine the trapezoidal integration method with the Filon integration method to calculate the coseismic displacement.It improves the speed of the calculation and there is no accuracy loss.It is useful for the inversion research based on the coseismic displacement.

Generalized reflection－transmission coefficient matrix，Discrete wavenumber，Coseismic displacement，Coulomb stress，Layered elastic model

10.6038／j.issn.0001－5733.2012.05.025

P315

2011－04－23，2012－03－28收修定稿

國家重點基礎(chǔ)研究發(fā)展計劃（2008CB425701）和國家自然科學(xué)基金（40974028；41030319）資助.

郝金來，男，1984年生，理學(xué)博士，主要從事震源機制與同震形變等方面的研究.E－mail：haojinlai＠gmail.com

郝金來，姚振興.均勻彈性分層介質(zhì)模型中的同震位移、應(yīng)變以及應(yīng)力.地球物理學(xué)報，2012，55（5）：1682－1694，

10.6038／j.issn.0001－5733.2012.05.025.

Hao J L，Yao Z X.The coseismic displacement，strain and stress in the layered elastic model.Chinese J.Geophys.（in Chinese），2012，55（5）：1682－1694，doi：10.6038／j.issn.0001－5733.2012.05.025.

（本文編輯 胡素芳）