6級左右歷史地震震源參數(shù)的模擬估計

譚毅培，陳棋福

1 中國地震局地球物理研究所，北京 100081

2 中國地震局地震預測研究所，北京 100036

3 天津市地震局，天津 300201

4 中國科學院地球深部研究重點實驗室，中國科學院地質(zhì)與地球物理研究所，北京 100029

1 引 言

我國歷史資料中有著大量歷史地震破壞程度的記載.與7級以上的大震或巨震相比，對6級左右中等強度歷史地震的震源參數(shù)研究相對并不充分.利用歷史資料定量估計得到較可靠的6級左右歷史地震震源參數(shù)，可有效彌補現(xiàn)代儀器記錄的地震資料相對短暫之不足，為強震危險性分析提供重要的基礎資料.

由歷史資料記載的地震宏觀破壞記錄確定相應烈度分布，并依據(jù)高烈度區(qū)的幾何中心、極震區(qū)和有感區(qū)面積及烈度分布形狀來估計歷史地震的震源參數(shù)（多為震中經(jīng)緯度），是普遍使用的方法［1－2］.但因當時人文、經(jīng)濟和地理環(huán)境差異導致史料記載的震害詳簡程度不一，由殘缺不全的歷史記錄推斷的6級左右地震烈度數(shù)據(jù)點相對稀少，且空間分布圖像不大完整，由此勾畫的烈度等震線分布不免帶入研究者的主觀假設，所給出的震中參數(shù)存在相當大的主觀性和不確定性.

為克服利用烈度數(shù)據(jù)估計震中位置和震級主觀性強以及不確定性無法定量估計的困難，Bakun和Wentworth［3］以美國加州地區(qū)為例，用已知震中位置和震級的現(xiàn)代儀器記錄地震作為訓練集，確定震級－震中距－烈度值衰減關(guān)系式，使用格點搜索試算法估計震中和震級，并通過對訓練集數(shù)據(jù)反復抽樣給出震中和震級估計的置信區(qū)間.Bakun等［4］近來又對方法做了改進，不通過訓練集而直接對地震烈度數(shù)據(jù)重采樣，使用Bootstrap方法估計震中和震級的置信區(qū)間.Bakun－Wentworth方法直接對烈度數(shù)據(jù)點擬合而放棄烈度分布圖，具有一定的可重復性和客觀性，因而得到了廣泛的應用，但其局限性在于必須依靠一定數(shù)量的現(xiàn)代中強地震來確定衰減關(guān)系，對如我國華北那樣歷史地震資料豐富，而有現(xiàn)代儀器記錄之破壞性地震相對缺乏的地區(qū)，Bakun－Wentworth方法應用較為困難.

通過理論計算不同震源參數(shù)在地表產(chǎn)生的烈度分布，并與烈度等震線圖相比較，也是估計歷史地震震源參數(shù)的一個途徑.陳培善［5］早在20世紀70年代就對利用位移表示定理的中遠場近似計算近場值進行了探索.Sirovich等［6］通過枚舉可能取值空間的設定震源參數(shù)，利用S波遠場輻射花樣計算各種設定震源產(chǎn)生的烈度分布圖，與歷史資料比較確定震源參數(shù)的方法（稱為KF方法），具有可估計多個震源參數(shù)，且不依賴現(xiàn)代中強地震記錄的優(yōu)勢.但KF方法僅使用S波遠場輻射花樣，沒有考慮地震波中體波低頻部分和面波的影響，計算出的雙峰型烈度等震線與一般高烈度區(qū)近橢圓烈度分布形態(tài)符合程度不佳，并且在擬合Ⅵ度以下的烈度等震線時情況較差［6］.

隨著計算機技術(shù)的不斷發(fā)展，利用地震波場模擬計算估計歷史地震震源參數(shù)的應用研究越來越多［7－10］.本研究將 Bakun－Wentworth方法直接擬合烈度數(shù)據(jù)點和KF方法參數(shù)枚舉的優(yōu)勢結(jié)合起來，引入近年來在震源破裂過程反演［11］和地震危險性分析［12］得到應用的模型選擇（model selection）和多模型統(tǒng)計推斷（multimodel inference）方法［13］，發(fā)展一種利用地震波場數(shù)值模擬和烈度數(shù)據(jù)估計多個震源參數(shù)及其不確定性的方法，克服KF方法僅使用S波高頻近似的問題和Bakun－Wentworth方法難以在現(xiàn)代中強地震活動記錄較少地區(qū)的使用問題.

2 模型選擇方法和算法流程

本方法將所有可能的震源參數(shù)組合進行枚舉，每個組合視作一個震源模型，通過模型選擇方法對所有可能的震源模型做出評估，給出各震源參數(shù)的權(quán)重值分布.

2.1 模型選擇的數(shù)學基礎

模型選擇方法通過計算各個模型與觀測數(shù)據(jù)之間的相對熵（relative entropy），對模型擬合數(shù)據(jù)的程度給出定量化估計.這里簡要介紹本研究所需的統(tǒng)計量及其表達式，對模型選擇理論詳細系統(tǒng)的論述見文獻［13］.

首先計算出觀測或推斷的數(shù)據(jù)點烈度值與理論模擬計算得到的烈度值間的方差

其中ISi和IBi分別為第i個點的模擬計算烈度值和烈度記錄值，N為烈度數(shù)據(jù)點總數(shù).接著采用適用于小樣本情況的赤池信息準則［13－14］來計算每個設定震源模型的AICc值，即各震源模型與烈度記錄間的相對熵

其中K為模型參數(shù)個數(shù)，例如本研究考慮震中經(jīng)度、緯度、震源深度和滑動角的情況下，K＝4.在所有設定震源模型的AICc值中找出最小的AICcmin，進而計算在給定某些參數(shù)下模型子空間中模型的AICc差值Δi

根據(jù)AICc差值可得到相應模型權(quán)重wi

其中構(gòu)造底數(shù)La＝1.5.如此得到的wi在相應震源模型參數(shù)空間的和為1，且每一個wi都不小于零.因而wi滿足完備性、非負性和可列可加性，由概率分布的公理化定義可將wi理解為一組概率分布.由此可對具體參數(shù)進行多模型統(tǒng)計推斷，即依據(jù)模型權(quán)重計算參數(shù)的權(quán)重值分布，估計相應的震源參數(shù).其計算方法類似于對多維離散隨機變量的全概率做邊緣化（marginalization），求某一參數(shù)維度的邊緣概率分布.

2.2 方法計算流程

本文考慮到發(fā)震斷層走向和傾角可充分借鑒地質(zhì)構(gòu)造調(diào)查、人工地震探測和烈度幾何分布等先驗估計，故僅關(guān)注震中位置、震源深度和滑動角構(gòu)成的四維震源參數(shù)模型空間.所采用方法首先計算震中位置搜索格點AICc值分布，并據(jù)此選擇可能的震中位置，再根據(jù)每個已選震中位置估計震源深度和滑動角.方法計算流程如下：

①震源深度、滑動角和發(fā)震斷層破裂方式在可能取值范圍以一定間隔進行枚舉，組合成為備選震源模型.

②計算每個備選震源模型在地表的理論地震波場，依據(jù)水平向峰值加速度（PGA）或峰值速度（PGV）與烈度的換算關(guān)系（如表1所列的中國地震烈度表GB／T 17742－2008或相應研究區(qū)的關(guān)系式）得到烈度分布.考慮到本研究使用的時域偽譜法（PSTD）［15］理論地震圖計算直接給出的是速度記錄及速度與地震波能量的相關(guān)性［7］，故本文采用將PGV轉(zhuǎn)化為烈度值的方式，而不采用將理論計算的速度記錄轉(zhuǎn)化為加速度再利用經(jīng)驗關(guān)系轉(zhuǎn)為烈度值的方式.

③首先利用式（1）和（2）計算每個震源模型的烈度方差及其AICc值，然后取每個震中位置相應模型的最小AICc值作為此震中格點的AICc值，得到震中搜索范圍內(nèi)的AICc值分布.

④根據(jù)震中AICc值分布，選擇AICc值較小的格點作為震中位置模型選擇初步結(jié)果，并采用式（3）和（4）計算震中格點的權(quán)重值.

⑤在步驟③得到的震源模型AICc值基礎上，對初步選擇的每個震中位置結(jié)果，依據(jù)相應的AICc值利用式（4）分別求得震源深度和滑動角的權(quán)重值分布，通過多模型統(tǒng)計推斷震源模型的震源深度和滑動角估計結(jié)果.

⑥結(jié)合地質(zhì)構(gòu)造、地震學和大地測量等研究成果，對以上震源參數(shù)模型選擇結(jié)果的合理性進行分析取舍，最后得到歷史地震震源參數(shù)估計的可能結(jié)果.

3 方法數(shù)值檢測

為檢測方法復原已知震源參數(shù)的能力，本研究基于表2所示的一維速度模型，將已知的震源參數(shù)作為雙力偶點源，采用時域偽譜法（PSTD）［15］來計算理論地震圖，并將理論波形的水平峰值速度（PGV）轉(zhuǎn)化成烈度值.所使用的轉(zhuǎn)化關(guān)系參考中國地震烈度表（表1），利用中值回歸得到烈度值與PGV的關(guān)系式為I＝3.26logPGV＋3.45，其中PGV單位是cm／s.對轉(zhuǎn)換所得的烈度值加入一定標準差σ生成的隨機噪聲擾動值后（取四舍五入后的整數(shù)），作為烈度數(shù)據(jù)點來估計震源參數(shù).

表1 中國地震烈度表（GB／T 17742－2008）Table 1 China earthquake intensity scale given by GB／T 17742－2008

表2 數(shù)值檢測使用的一維速度結(jié)構(gòu)模型Table 2 1－D structure used in this numerical test

鑒于我國6級左右歷史地震的烈度數(shù)據(jù)點記錄多為幾十個，故不失一般性以生成隨機分布的30個烈度數(shù)據(jù)點來進行數(shù)值檢測.所檢測震源參數(shù)的震中位置搜索格點間隔為2km，震源深度在3～20km范圍以1km為間隔取值，滑動角在0°～350°范圍以10°為間隔取值（滑動角取值規(guī)則參見文獻［16］）.

對簡單的一維速度結(jié)構(gòu)，所檢測震源參數(shù)模板的理論震中可取在正中間的位置，而理論震源深度和滑動角則分別參照1998年張北6.2級地震和2006年文安5.1級地震震源參數(shù)來設定，即分別取深度10km和滑動角130°［17］與深度9km和滑動角210°［18］.將參照的文安地震震源深度10km改為深度9km，是為了檢測不同震源深度模板的情況.這兩個設定的震源參數(shù)模板，一個是以逆沖為主兼有右旋走滑，一個以右旋走滑為主兼有正斷分量.檢測中在理論計算轉(zhuǎn)換的烈度值中加入均值為0、標準差σ分別為0.5和1的高斯分布隨機數(shù)作為觀測值，來模擬因地下三維速度結(jié)構(gòu)和場地效應及人為因素等導致的烈度數(shù)據(jù)的不確定性.

圖1 震源深度10km和滑動角130°的震源模板1000次測試的統(tǒng)計分析結(jié)果（a）加入不同標準差的隨機擾動后的烈度偏差統(tǒng)計直方圖；（b）只選取最大權(quán)重值推斷震中位置格點的頻次疊加圖；（c）震中估計結(jié)果中最大權(quán)重格點與理論震中位置的距離統(tǒng)計；（d）取權(quán)重值大于最大權(quán)重的30%的格點作為震中估計結(jié)果的頻次分布圖；（e）震源深度的最大權(quán)重估計值與理論值的偏差統(tǒng)計柱狀圖及其權(quán)重值分布；（f）滑動角最大權(quán)重估計值與理論值的偏差統(tǒng)計柱狀圖及其權(quán)重值分布.Fig.1 Statistical results of the numerical test with input parameters of focal depth 10km and rake 130°（a）Statistical distribution of intensity shift；（b）Statistical results of epicenter dislocation distance between maximum weight grid and the input location；（c）and（d）are frequency distributions of epicenter grids from model selection by the maximum weight and large than 30%of maximum weight respectively；（e）and（f）are weight distributions and misfit histograms between maximum weight solution and input value of focal depth and input value of rake angle respectively.

圖2 參數(shù)為震源深度9km和滑動角210°的模板1000次測試統(tǒng)計分析結(jié)果（a）—（f）同圖1.Fig.2 Statistical results of the numerical test with input parameters of focal depth 9km and rake 210°Others are same as Fig.1.

Bakun等［4］利用Bootstrap重采樣1000次計算分析歷史地震震源參數(shù)估計方法的穩(wěn)定性和誤差.本研究同樣對兩組震源參數(shù)模板分別進行1000次參數(shù)估算，圖1和圖2分別是兩組模板參數(shù)測試的統(tǒng)計分析結(jié)果.其中圖1a和圖2a以直方圖形式給出了理論轉(zhuǎn)換烈度加入不同標準差的隨機干擾后的烈度值與理論值的偏離情況，可見考慮了烈度偏差達到1度甚至2～3度的情形.圖1b和圖2b給出了以最大權(quán)重作為震中估計結(jié)果的頻次疊加圖，圖1c和圖2c則給出了以最大權(quán)重估計的震中格點與理論震中的偏離統(tǒng)計（統(tǒng)計間隔為震中搜索間隔的2 km）；由此二圖可見重復多次以最大權(quán)重估計的震中結(jié)果都聚集在理論中心點周圍，與理論中心點的距離越大（大于4km）而被作為估值的頻率越低.考慮到歷史地震資料和震中估算精度的不確定性，圖1d和圖2d統(tǒng)計給出了以權(quán)重值超過30%最大權(quán)重的格點作為震中估算結(jié)果的頻次分布，即表示格點在1000次模型選擇測試中被選中的次數(shù)；由圖可見，即使在擾動標準差σ＝1時，理論震中位置格點被選中的頻次分別達912次和920次，雖然可選的震中格點增多且聚集范圍有所擴大.1000次隨機生成烈度點并添加一定干擾值的震中測試，模型選擇的結(jié)果大量聚集在理論震中位置及其周圍，說明本文所采用的分析方法對震中位置估計在統(tǒng)計上具有一致性和統(tǒng)計無偏性，以5km估算精度統(tǒng)計的震中選擇結(jié)果的可信度超過了90%.

圖1e和圖1f與圖2e和圖2f分別是震源深度和滑動角的相應統(tǒng)計結(jié)果，中間的累積柱狀圖分別表示以最大權(quán)重估計的震源深度和滑動角與理論值的偏差數(shù)，不同顏色表示不同的偏差值；左右二圖展示的震源深度和滑動角權(quán)重值分布顯示：理論震源深度和理論滑動角的估計值權(quán)重最大，說明本文所采用的分析方法對震源深度和滑動角的估計在統(tǒng)計上也具有一致性和一定的無偏性.

4 方法檢測實例：2004年美國Parkfield6.0級地震

在數(shù)值檢測結(jié)果的基礎上，我們進一步以獲得較多強震記錄的2004年美國加州Parkfield6.0級地震為例來檢測該分析方法的可行性.2004年P(guān)arkfield 6.0級地震發(fā)生在右旋走滑的San Andreas斷層Parkfield段，地震破裂長度約為40km，由起始破裂點向北西和南東方向不對稱雙向破裂，向北西破裂約30km，向南東破裂約10km，破裂速度約3km／s［19－20］.Berkeley大學的快速矩心矩 張 量 解 （http：／／quake.geo.berkeley.edu／mt／nc51147892.mt）給出的該地震震源深度為8.0km和滑動角為185°（－175°）.全球 GCMT解（http：／／www.globalcmt.org）給出震源深度12km 和滑動角188°（－178°）.利用強震波形和GPS數(shù)據(jù)聯(lián)合反演震源破裂過程［19，21］得出的震源深度8.1km.

圖3 Parkfield地震強震觀測點與震源搜索模型示意圖（a）黑色三角為強震臺站位置，也是本文使用的烈度數(shù)據(jù)點位置，虛線網(wǎng)格為震中搜索格點.（b）三種破裂方式示意圖，分別是對稱雙向破裂和向兩側(cè)的不對稱雙向破裂.Fig.3 Intensity data distribution and epicenter search grid of the Parkfield earthquake（a）Black triangles are strong motion stations and also intensity data points locations used in the test，epicenter search grids are given by dash lines.（b）Three kinds of rupture styles with symmetric and asymmetric bilateral rupture in two directions.

我們從 COSMOS強震數(shù)據(jù)中心（http：／／db.cosmos－eq.org）收集到如圖3a和表3所示的Parkfield地震的39個強震臺站記錄波形，并將強震波形記錄的水平向峰值速度采用表4的對應關(guān)系轉(zhuǎn)化為烈度值，檢測本文方法的參數(shù)估計結(jié)果與利用波形反演得到的震源參數(shù)是否一致.

表3 本文使用的強震臺站名稱及經(jīng)緯度Table 3 Strong motion stations used in the Parkfield earthquake test

表4 PGV轉(zhuǎn)化為烈度值對應表1）Table 4 Relationship between instrumental intensity and PGV from USGS1）

本研究的理論地震圖計算以San Andreas斷層為界采用不同的一維地下速度結(jié)構(gòu)模型（表5），即已用于Parkfield地震震源破裂過程反演［19，21］的速度模型.

我們同樣采用 Ma等［19］的發(fā)震斷層垂直于地表和斷層走向140°的簡化假設，設斷層破裂長度為40km，并假設在整個斷層面上地震滑動角一致，斷層面上應力降分布為從起始破裂點向周圍衰減.備選震源模型構(gòu)建考慮三個因素：震源深度、滑動角和斷層破裂方式.震源深度取值范圍3km到15km，以1km為間隔枚舉，滑動角取值范圍0°到350°，以10°為間隔枚舉.如圖3b斷層破裂方式分為對稱雙向破裂和分別向兩側(cè)傳播的不對稱雙向破裂，起始破裂點分別取在斷層北西端的1／4、1／2和3／4處，稱為L方式、M方式和R方式.圖3a中網(wǎng)格為震中估計格點搜索范圍示意圖，搜索格點間隔為5km，震源球所指位置是GCMT解給出的震中位置.

圖4給出了Parkfield地震震源參數(shù)估計過程的主要計算結(jié)果，其中圖4a為震中搜索格點的AICc值分布；圖4b為震中模型權(quán)重值分布，權(quán)重最大的震中格點為圖中模型1所示的震中位置，其權(quán)重值為0.367；以不低于最大權(quán)重值的30%來選擇，則有如圖所示的4個震中位置估計模型，即模型1至模型4（權(quán)重值相應為0.367、0.16、0.108和0.108）.其中模型1權(quán)重值明顯大于其它三個，且與圖3a中GCMT解給出的震中位置一致.

表5 Parkfield地震地震動模擬計算所用的一維速度結(jié)構(gòu)模型Table 5 Crustal velocity structure used in the Parkfield earthquake test

圖4 Parkfield地震震源參數(shù)估計結(jié)果（a）震中格點搜索AICc值分布圖；（b）震中模型選擇結(jié)果，圓圈表示選中的4個模型震中位置；（c）以模型1為震中估計震源深度和滑動角的AICc值分布圖；（d）震源深度的權(quán)重值分布；（e）滑動角的權(quán)重值分布.Fig.4 Source parameters estimation result of the Parkfield earthquake（a）AICcdistribution of epicenter search grids；（b）Epicenter model selection results with 4circles；（c）AICcdistribution with different focal depths and rake angles for the epicenter Model 1；（d）and（e）are weight distributions with focal depth and rake angle，respectively.

以模型1的震中位置為依據(jù)，進一步估計Parkfield地震震源深度和滑動角.圖4c是震中位置取模型1位置時，不同震源深度和滑動角模型的AICc值分布，圖4d和圖4e分別是震源深度和滑動角的權(quán)重值分布.結(jié)果顯示震源深度在7km處取得最大權(quán)重值0.29，8km深度權(quán)重值略小為0.27，這一結(jié)果與Berkeley快速矩心矩張量解震源深度8.0km 和震源破裂過程反演［19，21］得出的8.1km震源深度結(jié)果基本一致，但與GCMT解給出的12km深度結(jié)果有差距.滑動角權(quán)重值分布曲線有兩個峰值區(qū)，一個代表的震源機制是左旋走滑兼有逆沖分量，另一個的震源機制則是右旋走滑兼有正斷分量，需要根據(jù)地質(zhì)構(gòu)造等其它信息選擇其中一個.考慮到San Andreas的右旋走滑斷層屬性，及1966年同一斷層段發(fā)生的6級地震的震源機制為右旋走滑并帶有較小的傾滑分量［22－23］，因此可以判定2004年P(guān)arkfield地震震源機制應取右旋走滑兼有正斷分量的峰值，含左旋走滑分量的峰值被舍去.滑動角在200°取得最大權(quán)重值0.23（舍去一個峰值則權(quán)重值比圖4e上的讀數(shù)增加一倍），在170°、190°和210°權(quán)重值相近約為0.11，220°和230°權(quán)重值也超過了0.05，因而滑動角估計結(jié)果為［170°，230°］，震源機制以右旋走滑為主，兼有少量正斷分量.此結(jié)果與Berkeley快速矩心矩張量解滑動角185°（－175°）和 GCMT解滑動角188°（－178°）的結(jié)果基本一致.

取模型選擇結(jié)果中權(quán)重值最大的震中位置、震源深度8km和滑動角200°，計算了39個臺站的理論地震圖（圖5）與實際記錄波形進行對比，兩個水平方向分別為平行斷層和正交斷層方向.波形長度取15秒，經(jīng)過0.2～1Hz的4階零相位Butterworth濾波器帶通濾波.模擬波形與實際記錄波形的符合程度是可以接受的，說明震源模型選擇結(jié)果與實際發(fā)震斷層性質(zhì)相近.

需要說明的是，由于本研究是以Parkfield 6.0級地震進行應用檢測，震中位置已有不少研究給出較確切的結(jié)果，故僅以模型1的斷層位置進一步估計震源深度和滑動角.而在應用于具體的6級左右歷史地震震源參數(shù)估計時，則需對出現(xiàn)的多個震中位置估計結(jié)果，分別推斷震源深度和滑動角，得到多組可能的估值.最后根據(jù)區(qū)域地質(zhì)構(gòu)造和大地測量等研究結(jié)果，綜合判斷這幾組估值的合理性.

5 方法應用實例：1882年河北深縣6級地震

圖5 Parkfield地震水平向模擬波形和記錄波形對比圖紅色為模擬波形，黑色為記錄波形.兩個方向分別為平行斷層和正交斷層，波形經(jīng)過0.2～1Hz的4階零相位Butterworth濾波器帶通濾波.Fig.5 Waveforms of horizontal components from observed ground motions（black）and synthetic calculations（red）with the best selected source model Left and right seismograms of each station are components parallel to and perpendicular to the fault.All the waveforms are band－pass filtered between 0.2to 1Hz with a zero－phase forth order Butterworth filter.

圖6 《中國歷史強震目錄》給出的1882年深縣6級地震烈度等震線圖（a）和深縣凹陷區(qū)斷裂體系分布圖（b，改自梁富康等［24］）Fig.6 Isoseismal map of the Shenxian earthquake given by China historical earthquake catalog（a）and regional faults systems in the Shenxian depression（b，modified from Liang et al.［24］）.

本研究進一步以有著較好歷史地震破壞記錄的1882年河北深縣6級地震為例，采用同樣的途徑對深縣地震的震源參數(shù)進行估計.《中國歷史強震目錄》［1］中給出的深縣地震（圖6a）震中位于38.1°N、115.5°E（震中精度屬于2類，誤差≤25km），極震區(qū)在深縣北部，較高烈度的Ⅷ度區(qū)和Ⅶ度區(qū)位于深縣凹陷區(qū)內(nèi).梁富康等［24］據(jù)反射勘探剖面編制的深縣凹陷區(qū)斷裂體系分布圖（圖6b）為我們的應用研究提供了構(gòu)造參考依據(jù).梁富康等［24］將深縣凹陷區(qū)的斷層規(guī)模分成三級，即2條一級斷層（衡水斷層和虎北斷層），4條二級斷層，從北到南依次為舊城北斷層、何莊斷層、深西斷層和深南斷層，及廣泛發(fā)育的三級或級別更小的斷層.梁富康等［24］提供的通過深縣凹陷主體的構(gòu)造剖面顯示，深縣凹陷區(qū)主要斷層為正斷性質(zhì)，區(qū)內(nèi)沉積層厚度在5km左右.衡水斷層為走向北西、傾向北東、具有左旋走滑分量的正斷層，斷層傾角較陡，一般55°左右，有些地段斷面陡直［24－25］.舊城北斷裂走向近東西，傾向南，何莊斷層也為東西走向，向南傾斜［24］.虎北斷層為北東東走向、北北西傾向，深西斷層為北東東走向、南南東傾向，深南斷層也是北東東走向，而其傾向為北北西［24］.

龍峰等［26］以系統(tǒng)收集整理獲得的1965年以來華北地區(qū)發(fā)生的34次中強地震的可靠破裂尺度參數(shù)，回歸得到華北地區(qū)6級左右地震斷層破裂長度在15km左右.據(jù)此判斷，圖6中的三級以下斷層規(guī)模無法成為深縣6級地震的發(fā)震斷層，4條二級斷層和2條一級斷層有可能是發(fā)生深縣6級地震的斷裂.據(jù)深縣地震的烈度等震線形態(tài)判斷，發(fā)震構(gòu)造似乎應為北西—南東走向的斷層，但北西向的衡水斷層與極震區(qū)距離在20km以上，且處于Ⅶ度區(qū)以外，似乎不具備成為深縣地震發(fā)震斷層的條件.因而1882年深縣地震的發(fā)震構(gòu)造是一個值得研究的問題.

5.1 烈度資料

依據(jù)前人整理的《河北省地震資料匯編》［27］和《中國歷史強震目錄》［1］，深縣地震給出有破壞程度記載的地點名稱42個，深縣、深澤等較大城鎮(zhèn)在烈度圖（圖6a）上有明確標識，可以以此為依據(jù)確定烈度值；高烈度區(qū)及周邊的村莊則根據(jù)周圍城鎮(zhèn)的烈度值和文獻記載中關(guān)于“輕”和“重”的描述來判定烈度值.深縣處于Ⅷ度和Ⅶ度區(qū)之間，則其中“最重”和“較重”的村判為Ⅷ度和Ⅷ度弱（郭家莊、北安、深州城關(guān)四廂、西八弓），“較輕”和“又輕”的村判為Ⅶ度（西馬莊、南官村、北官村）.束鹿處于Ⅶ度和Ⅵ度區(qū)之間，其中“最重”和“次重”的村判為Ⅶ度（雙井、北龐營、南龐營、北呂彩），“較輕”的判為Ⅵ度（沙河）.安平的烈度為Ⅵ＋，郭家店（郭屯）是周圍村莊中較重的，因此判為Ⅶ度.對于文字資料沒有單獨明確描述的，參考烈度等震線圖判斷烈度.“城北之郭家莊等村，城西之杜家莊等村，西北與安平、束鹿毗連之張村等村，同時震塌房屋，間有傷斃人口等情”.根據(jù)烈度等震線圖，張村在深縣西北，與圖中給出的Ⅷ度區(qū)很接近，與郭家莊一致判為Ⅷ度，杜家莊在深縣正西，與Ⅷ度區(qū)相距較遠，判為Ⅶ度.破壞程度記錄不明確和無法查明確切經(jīng)緯度的烈度點舍去.最終使用的深縣地震39個烈度數(shù)據(jù)點位置及烈度值見圖7和表6.

圖7 深縣地震烈度數(shù)據(jù)點分布圖Fig.7 Intensity data points distribution of the Shenxian earthquake

表6 1882年深縣地震烈度數(shù)據(jù)點位置分布Table 6 Intensity data points of the 1882Shenxian earthquake

5.2 震源參數(shù)估計

根據(jù)深縣地震的烈度數(shù)據(jù)點及其等震線和區(qū)域構(gòu)造分布（圖7和圖6），震源參數(shù)估計以三個可能的斷層走向來分析，分別為與衡水斷層相近的北西走向（走向設為135°）、與舊城北斷層和何莊斷層相近的近東西走向（設為90°）、及與深南斷層、深西斷層和虎北斷層相近的北東東走向（設為60°）.理論地震圖計算采用表2所示的速度模型，參考中國地震烈度表（表1）將理論波形水平峰值速度（PGV）轉(zhuǎn)化成烈度值，作為烈度數(shù)據(jù)點來估計震源參數(shù).

震源參數(shù)模型的地震震中搜索圍繞39個烈度數(shù)據(jù)點中心位置（115.677°E、38.179°N）的100km×100km范圍進行，搜索間隔5km，格點的排列方向與可能的發(fā)震斷層走向一致；震源深度則在3～20km之間以1km為間隔枚舉，斷層傾角在30°～90°之間以15°為間隔枚舉，滑動角在0°～360°間以10°間隔枚舉；斷層破裂長度設為20km，破裂方式分為從斷層兩側(cè)開始的不對稱雙向破裂和從斷層中間開始的對稱雙向破裂三種破裂模式；考慮到歷史地震震級估計的不確定性，以及淺部沉積層對地震波的放大作用，震級在5.5～6.5級間以0.05級為間隔枚舉.

圖8—圖10分別給出了以北西走向、東西走向和北東東走向斷層為可能的發(fā)震斷層的震源參數(shù)估計結(jié)果.發(fā)震斷層為北西走向的震中格點搜索結(jié)果（圖8）中，有12個格點的AICc為最小值－77.0，另有7個格點AICc值為次小的－68.3，這19個格點作為震中位置模型選擇的初步結(jié)果，如圖8b中圓圈位置.考慮到深縣凹陷區(qū)的沉積層厚度及地震較少發(fā)生在沉積層內(nèi)，因此舍去對應震源深度最大權(quán)重值小于5km的模型.若不同的震中格點對應模型為相同傾角和破裂模式，則只選擇其中一個AICc值最小的格點為此傾角和破裂模式對應的震中位置.則有4個格點成為北西走向發(fā)震斷層震中模型選擇的最終結(jié)果，稱作模型1—模型4，模型2和模型3震源機制為正斷為主兼有走滑分量，模型1和模型4以左旋走滑為主.

圖8 北西走向斷層震源參數(shù)估計結(jié)果（a）震中格點搜索AICc值分布；（b）震中模型權(quán)重值分布和震中模型選擇初步結(jié)果，圖中圓圈是初步選擇的震中位置；（c）震中模型選擇的可能結(jié)果，不同顏色標出了可能選擇模型的震中位置與發(fā)震斷層位置和破裂情況，對不同模型有所重疊的斷層稍作偏移以清晰展示結(jié)果；（d）—（g）相應震中選擇模型的震源深度和滑動角估計結(jié)果.Fig.8 Source parameters estimation results with NW direction strike focal structure（a）AICcdistribution of epicenter searching grid；（b）Model selection result of preliminary epicenters，the circles demonstrate preliminary model selection epicenter locations；（c）Model selection result of preferred epicenters named Model 1—4with relative fault and its rupture style；（d）—（g）Multimodel inference results of focal depth and rake angle for four models.

對東西走向發(fā)震斷層設定的震中格點搜索AICc值分布（圖9a），共有6個格點AICc取得最小值－68.3，13個格點取得 AICc次小值－61.2，這19個格點作為震中位置模型選擇的初步結(jié)果，如格點權(quán)重值分布圖9b中所示圓圈位置.對初步選擇的每個震中位置，使用多模型統(tǒng)計推斷估計震源深度和滑動角（圖9d—9g）.同樣對震中初步結(jié)果做多模型統(tǒng)計推斷，得到4個震中位置作為東西走向發(fā)震斷層震中模型選擇的最終結(jié)果（圖9c），稱作模型5—模型8，震中位置比較集中于深縣以北15km范圍內(nèi).考慮到研究區(qū)內(nèi)主要斷層的正斷屬性，因而在滑動角兩個峰值中取含有正斷分量的部分（180°～360°），4個模型都以正斷為主兼有少量走滑分量.

圖9 東西走向斷層震源參數(shù)估計結(jié)果（a）—（g）意義同圖8.Fig.9 Source parameters estimation results with EW direction strike focal structure Others are same as Fig.8.

發(fā)震斷層為北東東走向震中位置估計初步結(jié)果如圖10a，有2個格點的AICc為最小值－77.0，另有9個格點AICc值為次小的－68.3，這11個格點成為震中位置模型選擇的初步結(jié)果，5個模型成為最終的模型選擇結(jié)果，稱作模型9—模型13.震中位置分布比較集中，其中模型10、模型11、模型12和模型13震源機制以正斷為主，模型9兼有正斷和左旋走滑分量.

5.3 模型選擇結(jié)果討論

表7 深縣地震模型選擇結(jié)果匯總Table 7 Model selection results of the Shenxian earthquake

圖10 北東東走向斷層震源參數(shù)估計結(jié)果（a）—（k）意義同圖8.Fig.10 Source parameters estimation results with NEE direction strike focal structure Others are same as Fig.8.

圖11a為斷層走向北西的模型震中和斷層位置，模型1和模型4震源機制估計結(jié)果為走滑兼有少量正斷分量，與區(qū)域地質(zhì)構(gòu)造［24－25］北西走向衡水斷層正斷層屬性不符.模型2和模型3與衡水斷裂地表出露位置相距較遠，且其震源深度估計值不大于10km，因而很難將模型與衡水斷層聯(lián)系起來.

圖11b表示斷層走向為東西向的四個模型震中與斷層和烈度數(shù)據(jù)點的相對位置，四個模型震源機制估計結(jié)果都為正斷型，與深縣凹陷區(qū)內(nèi)區(qū)域地質(zhì)構(gòu)造環(huán)境相符.其中模型7和模型8震中位置與舊城北斷裂和何莊斷裂位置接近，模型5和模型6震中位置在舊城北斷層以北，而兩條斷層傾向為南南東，因而模型5和模型6震中位置似乎與東西向斷層不符.

圖11c中發(fā)震斷層為北東東走向的5個可能模型，震中位置都在三條斷層以北.其中模型11和模型12距離深西斷層較近，震源機制估計結(jié)果為正斷型，也符合區(qū)域地質(zhì)構(gòu)造環(huán)境.

圖11 深縣地震模型選擇結(jié)果綜合分析圖虛線圓圈表示震中模型，內(nèi)部數(shù)字表示模型號碼.深縣凹陷區(qū)一級二級斷層位置、烈度數(shù)據(jù)點位置和三個走向發(fā)震斷層震中模型選擇結(jié)果等表示同圖8至圖10.Fig.11 Source parameters model selection results of the Shenxian earthquake Dash circle demonstrate epicenter location with model ID.Others are same as Figs.8to 10.

綜合模型選擇結(jié)果和地質(zhì)構(gòu)造資料分析認為：深縣地震最可能的震源參數(shù)估計有兩組：一是圍繞模型3和模型4震中位置（震中格點中心分別為115.55°E、38.02°N，115.61°E、38.02°N），走向東西向，震源深度較深（最大權(quán)重值深度為15km和16km），震源機制為以正斷為主；另一組是模型11和模型12震中位置（震中格點中心分別為115.51°E、38.01°N，115.57°E、38.04°N），走向北東東，震源深度較深（最大權(quán)重值深度為14km和17km），震源機制為以正斷為主.東西向的舊城北斷層和何莊斷層，以及北東東走向的深西斷層為1882年深縣地震發(fā)震構(gòu)造的可能性較大.這一結(jié)果不同于圖6a所示的主要以Ⅵ度烈度等震線形態(tài)估計的北西走向發(fā)震構(gòu)造可能，其原因很可能在于對地處石家莊—保定間的定縣這一烈度點值的判定和Ⅵ度等震線的勾畫（圖7）.對因局部構(gòu)造或場地因素或人為因素出現(xiàn)的單一烈度異常點，若將其視為強震動的平均衰減結(jié)果則可能導致烈度等震線圖勾畫的偏差.

6 結(jié)論和討論

本文借鑒Bakun－Wentworth方法［3］直接擬合烈度數(shù)據(jù)點和KF方法［6］利用參數(shù)枚舉設定震源的做法，提出了一種基于烈度資料和地震動波場模擬估計歷史地震震源參數(shù)的模型選擇方法.數(shù)值恢復測試與Parkfield地震實例測試表明，本文提出的分析方法所估計的震源參數(shù)具有統(tǒng)計一致性和一定的無偏性.該方法可通過計算AICc值和模型權(quán)重，根據(jù)歷史地震烈度數(shù)據(jù)估計震中位置、震源深度和滑動角等震源參數(shù)及其不確定性，克服了Bakun－Wentworth方法需要利用研究區(qū)域的現(xiàn)代地震資料確定衰減關(guān)系式的問題，及KF方法在擬合Ⅵ度以下烈度等震線情況較差的問題，使方法具有一定客觀性和廣泛的適用性.

所提出方法一個顯著特點在于將利用歷史地震烈度資料估計震源參數(shù)的不確定性表示為多種可能結(jié)果的形式.與現(xiàn)代地震的波形記錄和烈度調(diào)查資料相比，由6級左右歷史地震破壞資料推斷得到的烈度數(shù)據(jù)不確定性較強，所含信息量較少，因而利用歷史地震烈度資料估計震源參數(shù)的結(jié)果也具有較大的誤差.KF方法［6］計算設定震源模擬烈度場與烈度分布圖的方差，僅將取得最小方差的震源參數(shù)范圍作為估計值，忽略比最小方差稍大的震源參數(shù)估值則可能拋棄潛在的結(jié)果.本文方法在分析震中格點搜索結(jié)果時，不僅選擇AICc值最小的震中格點，同時也將AICc值比最小值稍大的格點也納入震中模型選擇初步結(jié)果，同時通過計算每個震中格點的權(quán)重值，嘗試對震中模型的可能性做出定量化估計.將震源參數(shù)估計結(jié)果以多個可能的方式給出，充分考慮到烈度數(shù)據(jù)的不確定性，一定程度上增加了結(jié)果的可信度.

將本文提出方法應用于1882年河北深縣地震的震源參數(shù)估計，基于歷史資料推斷的地震烈度數(shù)據(jù)點，分析認為深縣凹陷區(qū)的北西向、東西向和北東東走向斷層對烈度數(shù)據(jù)點擬合程度接近，三個走向的斷層都有可能成為發(fā)震斷層.根據(jù)模型選擇結(jié)果和地質(zhì)資料的綜合分析，東西向的舊城北斷層或何莊斷層，及北東東走向的深西斷層成為發(fā)震構(gòu)造的可能性較大.

不可否認，受計算條件和研究程度的限制，本研究中將震源模型的發(fā)震斷層設定在矩形平面，與實際地震的發(fā)震斷層形態(tài)和可能較為復雜的斷層破裂方式有所差距；地震動模擬計算高頻地震波輻射還存在一定困難，可能會使由理論地震波場得到的烈度分布，與近斷層的烈度資料相差較大；以及由歷史記載的地震破壞程度來估計烈度和PGV或PGA與烈度轉(zhuǎn)換的不確定性，對本文提出的方法會產(chǎn)生一定的影響.地下結(jié)構(gòu)的橫向不均勻性，尤其是淺部沉積構(gòu)造導致的地震動效應，亦會影響模型選擇的結(jié)果.隨著地下速度結(jié)構(gòu)探測研究的深入，使用高精度的三維結(jié)構(gòu)模型進行波場模擬，或采用場地修正方法（site correction）對烈度數(shù)據(jù)點進行合理的修正，應可提高本文方法的應用效果.

（References）

［1］國家地震局震害防御司.中國歷史強震目錄（公元前23世紀—公元1911年）.北京：地震出版社，1995.Disaster Prevention Department of China Earthquake Administration.China Historical Earthquake Catalog （23rd Century BC—1911AD）（in Chinese）.Beijing：Seismological Press，1995.

［2］Musson R M W，Jeménes M J.Macroseismic estimation of earthquake parameters.Sixth Framework Programme Report，2008.

［3］Bakun W H， Wentworth C M.Estimating earthquake location and magnitude from seismic intensity data.Bull.Seismol.Soc.Am.，1997，89（2）：557.

［4］Bakun W H，Cópera A G，Stucchi M.Epistemic uncertainty in the location and magnitude of earthquakes in Italy from macroseismic data.Bull.Seismol.Soc.Am.，2011，101（6）：2712－2725.

［5］陳培善.震源機制與烈度分布的關(guān)系（二）.地球物理學報，1977，20（1）：9－19.Chen P S.Relation between the seismic source mechanism and the intensity distribution （continued）.Chinese J.Geophys.（Acta Geophysica Sinica）（in Chinese），1977，20（1）：9－19.

［6］Sirovich L. A simple algorithm for tracing synthetic isoseismals.Bull.Seismol.Soc.Am.，1996，86（4）：1019－1027.

［7］Fitzko F，Suhadolc P，Audia A，et al.Constraints on the location and mechanism of the 1511Western－Slovenia earthquake from active tectonics and modeling of macroseismic data.Tectonophysics，2005，404（1－2）：77－90.

［8］Nunziata C，Sacco C，Panza G F.Modeling of ground motion at Napoli for the 1688scenario earthquake.Pure.Appl.Geophys.，2010，doi：10.1007／s00024－010－0113－1.

［9］周青云，周仕勇，陳棋福等.正演推算1730年北京西郊地震的發(fā)震斷層和滑動角.地震研究，2008，31（4）：369－376.Zhou Q Y，Zhou S Y，Chen Q F，et al.Forward inference of the causative fault and the slip angles of the faults of the 1730 western Beijing earthquake.J.Seism.Res.（in Chinese），2008，31（4）：369－376.

［10］Ding Z，Romanelli F，Chen Y T，et al.Realistic modeling of seismic wave ground motion in Beijing city.Pure.Appl.Geophys.，2004，161（5－6）：1093－1106，doi：10.1007／s00024－003－2498－6.

［11］Emolo A，Zollo A.Kinematic source parameters for the 1989 Loma Prieta earthquake from the nonlinear inversion of accelerograms.Bull.Seismol.Soc.Am.，2005，95（3）：981－994.

［12］Scherbaum F，Delavaud E，Riggelsen C.Model selection in seismic hazard analysis：An information－theoretic perspective.Bull.Seismol.Soc.Am.，2009，99（6）：3234－3247.

［13］Burnham K P，Anderson D R.Model Selection and Multimodel Inference，A Practical Information－Theoretic Approach，2nd ed.New York：Springer－Verlag，2002.

［14］Akaike H.A new look at the statistical model identification.IEEE Trans.Autom.Control，1974，19（6）：716－723.

［15］Liu L B，Arcone S A.Propagation of radar pulses from a horizontal dipole in variable dielectric ground：A numerical approach.Subsurface Sensing Technologies and Applications，2005，6（1）：15－24.

［16］Aki K，Richards P G.Quantitative Seismology，2nd ed.Sausalito，California：University Science Books，2002.

［17］中國地震局監(jiān)測預報司.一九九八年張北地震.北京：地震出版社，1998：2.Monitoring and Forecasting Department of China Earthquake Administration.1998Zhangbei Earthquake （in Chinese）.Beijing：Seismological Press，1998：2.

［18］黃建平，倪四道，傅容珊等.綜合近震及遠震波形反演2006文安地震（Mw5.1）的震源機制解.地球物理學報，2009，52（1）：120－130.Huang J P，Ni S D，F(xiàn)u R S，et al.Source mechanism of the 2006 Mw5.1Wen′an Earthquake determined from a joint inversion of local and teleseismic broadband waveform data.Chinese J.Geophys.（in Chinese），2009，52（1）：120－130.

［19］Ma S，Custódio S，Archuleta R J，et al.Dynamic modeling of the 2004 Mw6.0Parkfield，California，earthquake.J.Geophys.Res.，2008，113：B02301，doi：10.1029／2007JB005216.

［20］Mendoza C，Hartzell S.Finite－fault analysis of the 2004 Parkfield，California，earthquake using pnl waveforms.Bull.Seismol.Soc.Am.，2008，98（6）：2746－2755.

［21］Custódio S，Liu P C，Archuleta R J.The 2004 Mw6.0 Parkfield，California，earthquake：Inversion of near－source ground motion using multiple data sets.Geophys.Res.Lett.，2005，32（23）：L23312，doi：10.1029／2005GL024417.

［22］McEvilly T V，Bakun W H，Casaday K B.The Parkfield，California，earthquakes of 1966.Bull.Seismol.Soc.Am.，1967，57（6）：1221－1244.

［23］Wu F T.Parkfield earthquake of June 28，1966：magnitude and source mechanism.Bull.Seismol.Soc.Am.，1968，58（2）：689－709.

［24］梁富康，于興河，李先平等.冀中坳陷深縣凹陷的生長斷層特點及其對沉積的控制作用.中國地質(zhì)，2011，8（2）：263－270.Liang F K，Yu X H，Li X P，et al.Growth faults in Shenxian depression and their control over the sedimentation.Geology in China （in Chinese），2011，8（2）：263－270.

［25］楊旭升，劉池陽，楊斌誼等.冀中坳陷衡水轉(zhuǎn)換斷裂帶特征及演化.煤田地質(zhì)與勘探，2004，32（3）：5－8.Yang X S，Liu C Y，Yang B Y，et al.Characteristics and evolution of Hengshui transfer fault zones in Jizhong depression.Coal Geology ＆ Exploration （in Chinese），2004，32（3）：5－8.

［26］龍峰，聞學澤，徐錫偉.華北地區(qū)地震活斷層的震級－破裂長度、破裂面積的經(jīng)驗關(guān)系.地震地質(zhì)，2006，28（4）：511－535.Long F，Wen X Z，Xu X W.Empirical relationships between magnitude and rupture length，and rupture area，for seismogenic active faults in North China.Seismology and Geology （in Chinese），2006，28（4）：511－535.

［27］張秀梅.河北省地震資料匯編.北京：地震出版社，1990：324－327.Zhang X M.Hebei Earthquake Material Collection （in Chinese）.Beijing：Seismological Press，1990：324－327.