陳學忠，李艷娥，陳麗娟

1 中國地震局地球物理研究所，北京 100081 2 重慶市地震局，重慶 401147

0 引言

1976年7月28日在河北唐山發(fā)生的MS7.8巨大地震給國家和人民造成了巨大災難.地震預報，特別是短臨地震預報，可以顯著減輕地震造成的人員傷亡和財產損失.但是，由于科技水平的局限性、對地震產生的成因以及震前觀測到的一些異常現象的認識，還有待深入.至今，地震預報仍然難以實現.對已發(fā)生的強震進行解剖，深入研究強震發(fā)生前的異?，F象，將有助于深化對地震孕育發(fā)生過程的認識，同時對推進地震預報深入研究也有裨益.唐山MS7.8地震發(fā)生之后，不少研究者對唐山地震前出現的地震學方面的異常進行了分析研究(朱傳鎮(zhèn)等，1977,1981；李全林等，1978；林懷存等，1987；錢兆霞，1988；薛艷和梅世蓉，1999；隗永剛等，2018；李赫等，2020)，研究內容主要涉及區(qū)域地震活動圖像、震群活動、b值、微震震源參數以及地震活動與地球自轉的關系等.關于唐山地震前b值的研究，李全林等(1978)研究了唐山震區(qū)及其鄰區(qū)b值隨時間的變化，發(fā)現震前b值存在下降過程，同時還對京津唐渤張區(qū)域進行了b值空間掃描分析，發(fā)現低b值區(qū)由西向東遷移的過程.目前，在b值分析時，研究者們主要在兩個方面開展工作：一個是特定時段的空間分布分析，另一個是選定空間范圍內的時間變化分析.b值的空間分布可以反映各個地區(qū)b值的差異性，僅僅根據b值的空間分布難以辨別出某一地區(qū)b值隨時間的變化情況.在b值的空間分布圖上，常常出現b值高低區(qū)域交錯分布，無規(guī)律可循，難以將其與地震孕育相聯(lián)系.同樣，在選定某空間范圍進行b值隨時間變化分析時，選取區(qū)域時沒有合理的標準，具有任意性，難以保證所選取的區(qū)域完全包含了b值異常區(qū)，因而不能反映全部b值異常區(qū)的b值隨時間的變化情況.假設地震發(fā)生前在震源附近的一個空間區(qū)域內存在b值異常變化，如果選取的空間范圍位于這個空間區(qū)域以外，與地震孕育過程無關，這時，b值隨時間的變化與地震的孕育之間難以建立有效的聯(lián)系；如果選取的空間范圍只是包含了這個空間區(qū)域的一部分，此時b值隨時間的變化也不能全面反映地震的孕育過程.一些震例研究表明，地震發(fā)生前b值存在數月-數年的下降變化過程(Imoto,1991; Jaumé et al.,1999; Enescu and Ito,2001; Nuannin et al.,2005; Chan et al.,2012; Nanjo et al.,2012; 史海霞等,2018).我們著眼于b值下降變化過程，通過格點搜索把b值下降變化的區(qū)域搜索出來，然后再分析這個區(qū)域內b值隨時間的變化特征，就可以全面地揭示b值異常的空間分布和時間變化特征.本文將根據上述思路對唐山地震前的b值異常特征進行研究.

1 資料和研究區(qū)域

圖1a給出了1976年唐山MS7.8地震及其余震的震中位置和研究區(qū)域，取圖中所示的經緯度范圍(37.5°N—42°N，113.5°E—121°E)為格點搜索范圍，即研究區(qū).這個區(qū)域顯然遠遠大于唐山MS7.8地震余震分布區(qū)域.對研究區(qū)內于1970年1月—1976年6月期間發(fā)生的地震，根據G-R關系得到的地震目錄完整性震級可達ML1.5.因此，在本文分析中，選取研究區(qū)內1970年1月—1976年6月ML≥1.5的地震目錄作為研究資料，共2214次地震.地震目錄來自于中國臺網中心匯編的全國小地震目錄庫.圖1b給出了這些地震的空間分布，其震級隨時間的變化見圖2.

圖1 (a)唐山MS7.8地震震中與余震分布,灰色圓圈為余震分布；(b)選擇的研究區(qū)域范圍與1970年1月—1976年6月期間發(fā)生的地震(ML≥1.5)震中分布,灰色圓圈為地震震中.紅色六角星為1976 唐山MS7.8地震.震源機制解來自于GCMT目錄

圖2 ML≥1.5地震震級隨時間變化

2 唐山地震前b值變化特征

2.1 分析方法

采用極大似然法計算b值(Aki，1965):

(1)

b值的95%置信度的標準差為

(2)

先將研究區(qū)域劃分成0.1°×0.1°的格點，以每個格點為中心點，選取2°×2°的空間子區(qū)域，計算每個子區(qū)內地震的b值隨時間的變化曲線.然后，挑選出震前b值呈明顯下降變化形態(tài)的子區(qū)，并以該子區(qū)的中心點作為一個異常點，可以把多個異常點集中的區(qū)域看成是b值異常區(qū).然后，選取這個異常區(qū)內的地震計算b值隨時間的變化，這樣就可以得到整個b值異常區(qū)的范圍及其b值隨時間變化的特征.

計算每個子區(qū)b值隨時間的變化時，使用的是等事件數的樣本窗口，將該窗口以相等事件數的增量進行滑動.因b值是一個統(tǒng)計量，一般會受樣本量的影響，選取等事件數量的窗口，目的是為了確保分析結果不受樣本大小差異的影響.

2.2 結果

對于每個子區(qū)，以100個地震計算b值，5個地震滑動，計算b值隨時間的變化曲線.共計算了2392個子區(qū)的b值隨時間的變化.有的子區(qū)b值個數較少，其時間變化曲線難以反映出可靠的b值變化信息.我們選取具有15個及以上b值數據的子區(qū)進行分析，共1747個子區(qū).圖3a中給出了上述2392個計算了b值隨時間變化的子區(qū)中心點位置，其中灰色圓點“”為b值數據小于15個的子區(qū)，紅色圓點“”為具有15個及以上b值數據的子區(qū).從1747個子區(qū)的b值隨時間的變化曲線中找出震前b值呈下降變化形態(tài)的子區(qū)，靠人工是很費事的.我們在b值時間變化曲線上取定1973年1月—1976年6月間的m個b值數據，對于相鄰兩個b值，計算出后1個b值與前1個b值之差，一共得到m-1個b值差值.再計算出m個b值中最后1個b值與第1個b值之間的相對變化量.挑選b值呈下降變化形態(tài)的子區(qū)的條件設定為：(1)m-1個b值差值中有55%為負值，即多數b值是下降的；(2)相對下降變化量≥10%.根據這兩個條件，可用程序自動挑選出滿足設定條件的子區(qū).圖3b給出了按設定條件挑選出來的子區(qū)的中心點的空間位置，以綠色圓點表示，共391個點.這些綠色圓點中有376個點集中分布在京—津—唐一帶地區(qū).我們從376個點中隨機抽取了10個點，圖4中給出了它們所在的子區(qū)的b值隨時間的變化，還給出了各子區(qū)的中心點的經緯度，豎直線“│”表示唐山MS7.8地震的發(fā)生時間，曲線末端的數字為10個子區(qū)的序號，這些序號所在的空間位置也用紅色數字標注在圖3b中.從圖4中可以看出，這10個子區(qū)的b值在震前確實呈下降變化形態(tài)，說明前述設定的選擇條件是可行的.選取圖3b中黑色粗實線圍成的區(qū)域，這個區(qū)域比376個綠色圓點集中區(qū)稍大一點，可以基本包含整個b值下降異常區(qū)，以100個地震計算b值，5個地震滑動，可以得到b值隨時間的變化曲線(圖5).可以發(fā)現，在唐山MS7.8地震發(fā)生前，這個區(qū)域內b值大約在1973年6月從1.1開始下降，到1974年10月達到最低值0.72，相對下降幅度約為35%.從1974年11月到1976年6月，b值維持在低位，沒有明顯變化.

圖3 (a)b值時間變化的子區(qū)中心點位置分布.灰色圓“”為b值數據小于15個的子區(qū)，紅色圓“”為具有15個及以上b值數據的子區(qū)；(b)b值時間變化呈下降變化形態(tài)的子區(qū)中心點位置分布(綠色圓點).六角星為唐山地震震中

圖4 隨機選取的10個格點為中心點，2°×2°的空間子區(qū)域內地震的b值隨時間的變化

圖5 b值(黑色粗實線)隨時間的變化

綜上所述，唐山地震前在京—津—唐一帶出現大面積b值下降異常，從1973年中到1974年10月，b值呈趨勢下降變化形態(tài).之后，直到臨近地震發(fā)生時，b值穩(wěn)定在低值，沒有明顯變化.在空間上異常面積較大，在時間上表現為先下降，然后維持在低值平穩(wěn)變化，這就是唐山地震前b值的變化特征.

3 結論與討論

著眼于b值隨時間變化過程中的下降變化形態(tài)，采用格點搜索法，將研究區(qū)域劃分成0.1°×0.1°的格點，以每個格點為中心點，選取2°×2°的空間子區(qū)域，計算每個子區(qū)b值隨時間的變化曲線.根據各子區(qū)b值隨時間變化曲線，得到了唐山MS7.8地震前b值下降異常分布區(qū)域.b值下降異常集中分布在京—津—唐一帶地區(qū)，位于(116°E—119°E，38.8°N—41°N)的經緯度范圍內，所占面積大約為35000 km2.

在b值下降異常區(qū)內，唐山MS7.8地震發(fā)生前b值大約從1973年6月到1974年10月經歷了17個月的下降過程，從1.1降至0.72，相對下降幅度約為35%.從1974年11月—1976年6月的20個月中，b值一直維持在低值，沒有明顯變化.一般把余震分布區(qū)作為地震的破裂區(qū)，唐山地震的破裂區(qū)約占b值異常區(qū)內的26%.震前破裂區(qū)內發(fā)生的ML≥1.5總共111個，顯然難以給出破裂區(qū)內可靠的b值隨時間變化曲線來.

為了考察計算b值的地震個數的影響，分別以150、120和80個地震為一組計算b值，以5個地震滑動得到b值隨時間的變化.以同樣的挑選b值呈下降變化形態(tài)的子區(qū)的設定條件進行挑選，得到的b值呈下降變化形態(tài)的子區(qū)的中心點分布如圖6.當計算b值的地震個數增加時，滿足挑選條件的子區(qū)數會減少，b值呈下降變化形態(tài)的子區(qū)的中心點也會有所減少；當計算b值的地震個數減少時，滿足挑選條件的子區(qū)數會增多，b值呈下降變化形態(tài)的子區(qū)的中心點略有增加.總體上來講，b值下降異?？臻g分布概貌沒有本質的變化.

圖6 b值時間變化呈下降變化形態(tài)的子區(qū)中心點位置分布(綠色圓點)

b值變化除了與應力有關外，還與介質的均勻度和熱梯度有關.應力增加，b值降低(Scholz，1968；Wyss，1973)；介質越均勻，b值越低(Mogi, 1962)；熱梯度降低，b值也降低(Warren and Latham, 1970).一般認為強震發(fā)生前存在區(qū)域應力增強過程，所以認為震前幾個月或幾年的b值下降變化與應力有關，這更容易理解.2011年日本本州MW9.1地震前，震區(qū)b值自2005年起下降，一直持續(xù)至臨近地震發(fā)生前(Nanjo et al.，2012).同期，視應力呈上升趨勢變化(李艷娥和陳學忠，2011).即日本本州MW9.1地震前，視應力與b值呈反向變化，即負相關.視應力上升，同時b值下降.視應力與b值之間的負相關關系是構造應力增強的結果.因此，盡管缺乏唐山地震前視應力變化結果，但仍然可以把b值的下降過程理解為構造應力的增強.

選取圖3b中黑色粗實線圍成的區(qū)域內的地震，以40天為時間窗口統(tǒng)計地震的應變釋放，然后將統(tǒng)計的應變釋放進行累加，可以得到累積應變釋放曲線，如圖7a.根據累積應變釋放數據，以25個數據為窗口，1個數據滑動，通過直線擬合可以得到累積應變釋放曲線的斜率K隨時間的變化，如圖7b.對于應變減速釋放，斜率K逐漸減??；對于應變加速釋放，斜率K逐漸增大.圖7b中顯示，1975年4月以前，K逐漸減小，應變釋放是減速的；從1975年4月直到臨近地震發(fā)生時，K逐漸增大，應變釋放是加速的.當K逐漸增大時，b值維持在低值，沒有明顯變化.因此，低b值條件下的應變釋放加速可能更具有預測意義.

圖7 (a)累積應變釋放曲線；(b)累積應變釋放曲線斜率K隨時間的變化.“↓”為唐山地震