山東地區(qū)b值的統(tǒng)計分析

薛靖楠,李 想,張 波,王永旺

(山東大學數(shù)學與系統(tǒng)科學學院,濟南 250014)

山東地區(qū)b值的統(tǒng)計分析

薛靖楠,李 想,張 波,王永旺

(山東大學數(shù)學與系統(tǒng)科學學院,濟南 250014)

利用山東省內(nèi)及其鄰近海域(34°～39°N,114°～125°E),1996年1月至2009年5月的地震資料,通過最小二乘法及極大似然法對山東地區(qū)近年來的b值進行了求解及分析,給出了本區(qū)b值的計算結(jié)果,討論了兩種方法得到結(jié)果的差異,并對近年來b值時序變化的一些特點及其所反映的問題進行了討論。

b值;最小二乘法;極大似然法;地震活動;統(tǒng)計分析

0 引言

反映地震頻度與震級關(guān)系的古登堡-里克特公式:lgN=a-bM是地震學中非常重要的經(jīng)驗公式。這個關(guān)系表明,大小地震是按一定比例發(fā)生的。其中,M是震級,N表示一定時間段中所統(tǒng)計地區(qū)內(nèi)發(fā)生的震級大于或等于M的地震個數(shù)。a、b是常數(shù),在短時間內(nèi)同一地區(qū)內(nèi)可以近似認為是恒定的[1]。

b值(即頻度-震級關(guān)系式中的斜率)具有明確的物理意義,反映大小地震間的比例關(guān)系。實驗表明,其主要與區(qū)域應(yīng)力水平和地殼結(jié)構(gòu)的均勻程度有關(guān)[2]。當區(qū)域應(yīng)力基本不變時,隨著地殼結(jié)構(gòu)的均勻程度增強或減弱,b值相應(yīng)減小或增強;當介質(zhì)結(jié)構(gòu)基本不變時,隨著區(qū)域應(yīng)力水平的增強或減弱,b值相應(yīng)減小或增大。在一般情況下,b值相對穩(wěn)定,在1附近波動;其中,大地震前 b值會發(fā)生異常變化,通常表現(xiàn)為減小。所以,根據(jù)b值的動態(tài)變化我們可以預(yù)測大地震的發(fā)生[3]。其次,根據(jù)古登堡公式可以推測出未來一段時間內(nèi)該地區(qū)發(fā)生最大地震的震級(即橫坐標的截距)[4]。因而,我們可以說b值在整個地震學研究中具有舉足輕重的地位。

1 方法與資料介紹

現(xiàn)在地震研究中,用來計算 b值的方法通常有2種:最小二乘法與極大似然法[5]。

1.1 最小二乘法

若我們假設(shè) x和y之間的關(guān)系表示成y=a+bx,那么,在 x取 x1,x2,…,xn時,y的取值分別為a+bx1,a+bx2,…,a+bxn,而實際觀測到的數(shù)據(jù)是

y1,y2,…,yn,因此提出目標量L(a,b)=(yi-(a+bxi))2,L(a,b)實際上就是 n對點的偏差平方和。我們要找 a,b,使L(a,b)達到最小。由:

1.2 極大似然法

設(shè)總體ζ服從參數(shù)為λ的指數(shù)分布,密度函數(shù)為f(x;λ)=λe-λx,λ >0,x>0。試求參數(shù)λ的極大似然估計量。

我們知道,λ的似然函數(shù)為

1.3 相關(guān)系數(shù)檢驗

在最小二乘法擬合中,最簡單的一種檢驗是相關(guān)系數(shù)檢驗,相關(guān)系數(shù)的定義是

其中 r2越接近于1,說明吻合程度越好

1.4 研究區(qū)域簡介

山東地區(qū)位于華北地震區(qū)東南部。本次研究的區(qū)域范圍為 114°～125°E,34°～39°N。其主要涉及山東內(nèi)陸及其附近海域。中國東部規(guī)模最大的北北東向郯廬斷裂帶縱貫山東省中部,北北東向聊考及滄東斷裂帶展布于魯豫冀交界處,北西向燕山—渤海—威海斷裂帶穿越渤海直達膠東半島北部,并與郯廬斷裂帶交匯在渤海中部。這些大型活動斷裂帶,也是地球物理異常帶及深部地殼變異帶。山東內(nèi)陸及近海還分布著其它許多規(guī)模不等、方向不一的次級活動斷裂(圖1)。這樣特定的地質(zhì)構(gòu)造背景,決定了山東內(nèi)陸及近海是一個多震地區(qū)。自有歷史地震記載以來,山東及附近海域共發(fā)生5級以上地震70余次(圖2),其中7級以上地震7次,約占華北7級以上地震次數(shù)的γ3。其中公元1668年的郯城8.5級地震更是中國東部強度最大的地震,造成山崩地裂、房屋夷平、死亡5萬余人,成為曠古奇災(zāi)。

圖1 山東地區(qū)構(gòu)造分布圖

圖2 山東地區(qū) MS≥4 3/4地震分布圖

2 通過最小二乘法和極大似然法計算b值及最大震級的估計

2.1 利用最小二乘法計算值

使用山東省內(nèi)及其鄰近海域(34°～39°N,114°～125°E)1996年1月至2009年5月期間的地震目錄數(shù)據(jù),通過最小二乘法進行了地震b值的計算,結(jié)果如圖3。

貨幣是指任何一種可以執(zhí)行交換媒介、價值尺度、延期支付標準和完全流動的財富儲藏手段等功能的商品，貨幣是商品交換發(fā)展到一定階段的產(chǎn)物。貨幣的本質(zhì)就是一般等價物。貨幣可以分為現(xiàn)金、存款、央票、短期銀行理財、貨幣基金等；

圖3 1996年1月至2009年5月山東地區(qū)震級頻度關(guān)系圖

首先,我們分別統(tǒng)計了上述范圍內(nèi)大于山東地震臺網(wǎng)可以監(jiān)控的震級 Mi(Mi=2.0,2.1,……,5.5)的地震次數(shù) Ni,然后以震級 M為橫坐標,lgN為縱坐標畫出散點圖。如圖3所示。通過圖像我們可以非常清晰地看出散點幾乎位于一條直線上,這表明了經(jīng)驗公式的準確性。

接著我們通過最小二乘法,綜合上述數(shù)據(jù),計算出了b值,結(jié)果為0.9386,這個b值可以視作山東省近年來b值的平均值或基準值,具體計算方法如下:

由1.1最小二乘法可直接推得古登堡公式lgN=a-bM中a,b的估計值

計算其相關(guān)系數(shù)的平方,約為0.9920,說明直線吻合得很好。

最后,通過延長直線到橫軸,我們得到了其在橫坐標上的截距,約為5.74。這即是未來山東地區(qū)所發(fā)生地震最大震級的估計值

2.2 計算b值的誤差分析及克服方法

保證b值準確的條件是:資料準確,樣本數(shù)足夠,計算方法合理。但是通常情況下,這三方面都存在一定的問題,因而導致計算結(jié)果誤差較大。為了能夠在最大限度內(nèi)減小誤差,增強b值的計算準確度,我們采取了拋投法。

由于b值的具體數(shù)值與起算震級的選取有關(guān)系,因而選取合適的起算震級可以極大地減小誤差。這里我們采取拋棄lgN值較小的點(即提高起算震級)來克服“掉頭”現(xiàn)象[6](即在 X值較小時,散點常常偏離了直線,接近于曲線。通常認為導致這一現(xiàn)象的主要原因是小震漏記),當然,起算震級也不能選的太高,以免丟失信息。我們根據(jù)前人的經(jīng)驗,起算震級從1.5級開始,每隔0.1級做一次 b值計算,并計算出相關(guān)系數(shù)的平方。結(jié)果發(fā)現(xiàn),相關(guān)系數(shù)從較小值逐漸增大,達到一極大值時又逐漸減小。其中,最大值集中在2.0附近,這與計算得到的臺網(wǎng)監(jiān)測能力基本一致,因而我們選取2.0作為起算震級。

2.3 利用極大似然法計算b值

由lgN = a-bM,我們可以推出 N =eln10(a-bM)(M>M0),從而震級M的概率密度:

通過觀察可知,此式即為1.2中密度函數(shù)當 x=M-M0,λ=bln 10的特例,因而通過解似然方程可得:

代入數(shù)據(jù)(依然取起算震級=2.0),求得 b=0.790,這個b值可以視作山東省近年來以極大似然法計算得到的b值的平均值或基準值。

3 近年來b值的長期變化規(guī)律及其意義

研究的空間和時間范圍仍與第一節(jié)相同。

取12個月為單位時間段,1個月為滑動步長,用每12個月的地震記錄通過最小二乘法(起算震級仍取2.0)計算出一個 b值,總共得到了150個 b值數(shù)據(jù),按時間順序連成曲線,即圖4。注意圖中標明的時間是計算b值的終止時間(縱坐標為 b值,橫坐標為時間t)。

為了更好地研究b值變化與較大地震的對應(yīng)關(guān)系,我們在曲線圖上標明了所有不小于4.5級地震的發(fā)生時間(地震參數(shù)見表1)。

表1 山東地區(qū)1996—2009年4.5級以上地震表

極大似然法的 b值計算結(jié)果見圖5,由曲線變化形態(tài)看,與圖4最小二乘法的結(jié)果存在一定差異。分析產(chǎn)生這種差異的原因,出現(xiàn)該差異的原因是由于所采用方法對離散數(shù)據(jù)所設(shè)置的偏離權(quán)重不同,最小二乘法以離散點到擬合直線的距離為權(quán)重,而極大似然法基于概率分布函數(shù),當數(shù)據(jù)的概率分布函數(shù)是正態(tài)分布時,極大似然估計與最小二乘估計相同,但在實際應(yīng)用中,常常會因為我們用于計算b值的數(shù)據(jù)并不滿足正態(tài)分布,因而會出現(xiàn)二者計算結(jié)果的差異。謝華章通過數(shù)值模擬對極大似然法計算b值的條件進行了詳細研究,認為震級范圍和分檔間隔以及參與統(tǒng)計的地震數(shù)都對計算結(jié)果存在影響,并且認為當參與統(tǒng)計的地震數(shù)目較少時,由于個別較大地震隨機性的影響,用最小二乘法計算b值可能會產(chǎn)生相當大的系統(tǒng)偏差[7]。

圖4 山東地區(qū)b值變化曲線圖(1996年1月至2009年5月,M≥2.0級地震)

圖5 山東地區(qū)極大似然法b值變化曲線圖

從兩種方法的b值曲線圖上我們可以看出一些初步結(jié)果:

(1)4.5級以上的地震幾乎全部集中在曲線的低b值區(qū)域,對應(yīng)關(guān)系十分明顯。如果按第1節(jié)中計算出的最小二乘法b值基準值0.9386,取參考值為0.9,那么可以從圖4中看出,所有ML4.5級地震發(fā)生時,b值均在參考值之下。同樣,在極大似然法的b值變化曲線中,該方法的 b值基準值0.790,取參考值為0.7,則從圖5中也同樣可以看出,幾乎所有的ML4.5級以上地震也都發(fā)生在接近該參考值或發(fā)生在參考值以下。巖石破裂實驗為根據(jù)b值動態(tài)變化預(yù)測地震提供了物理基礎(chǔ)。肖爾茨最先從巖石破裂實驗中得出,巖石處于高應(yīng)力狀態(tài)時,b值低[8],大量的震例研究中發(fā)現(xiàn)較大地震前 b值下降[9-10]。

實際上,這種現(xiàn)象并不難理解。我們知道,大地震發(fā)生之前一段時間內(nèi)巖石會發(fā)生一些裂隙或微破裂,為大地震即將發(fā)生的前兆。起初,巖石破裂程度較小,所以常常只會引發(fā)一些較小的地震,小地震的數(shù)量增多會導致震級-頻度關(guān)系圖中直線下降更陡,即斜率更大,因而此時 b值較大;而后(即在大地震發(fā)生前較短時間內(nèi)),巖石破裂程度增大,此時中等地震增多,這段時間內(nèi)直線下降趨勢勢必變緩,即斜率變小,因而此時b值較小。這就說明了為什么大地震常常發(fā)生在低b值之后。

當然,這只是一個普遍情況。實際上,有時候大地震也會發(fā)生在高b值之后(例如本次調(diào)查中2008年3月的4.8級地震發(fā)生時對應(yīng)的最小二乘法b值就處于相對較高的位置,但此時極大似然法的b值曲線位于極小值)。這種情況一般出現(xiàn)在與地震震級相比研究區(qū)范圍較大,或地震發(fā)生在研究區(qū)邊緣、外圍。2002年4月22日的隆堯5.4級和2008年3月10日的蘭考4.8級地震均屬于這種情況。

4 結(jié)束語

G-R關(guān)系是迄今地震活動性研究中普適性最好的統(tǒng)計關(guān)系之一,b值物理意義明確,簡單易求,現(xiàn)在被不少人用于地震預(yù)測研究,尤其在中國應(yīng)用十分廣泛。本文利用現(xiàn)有地震資料,通過最小二乘法及極大似然法對山東地區(qū)近年來的b值進行了求解及分析,確定了該區(qū)b值的基準值,結(jié)果表明山東地區(qū)1996年以來12次 ML4.5級以上地震都發(fā)生在b值基準值以下。該結(jié)果有助于在日常分析預(yù)報中確定選擇合適的低b值作為異常指標。當然,由于可能存在一些小震的疏漏及方法本身固有的局限性,因而本次b值計算結(jié)果多少存在一些誤差。同時由于資料的時間尺度短,和我們對震前b值變化機理的認識尚待深入,因此文中的分析和結(jié)論也有待進一步驗證。最后,衷心希望本文的研究結(jié)果能對山東地區(qū)的地震預(yù)報和地震安全性評價工作有一定的參考價值。

致謝:感謝周翠英研究員和鄭建常副研究員的指導和幫助。

[1] 李全林.地震頻度-震級關(guān)系的時空掃描[M].北京:地震出版社,1979:1-102.

[2] 傅征祥,呂曉健,邵輝成,等.中國大陸及其分區(qū)余震序列b值的統(tǒng)計特征分析[J].地震,2008,28(3):1-7.

[3] 黃德瑜,張宇霞.b值時空掃描與地震預(yù)報[C]//地震監(jiān)測與預(yù)報方法清理成果匯編(測震學分冊).北京:地震出版社,1989:215-223.

[4] 吳開統(tǒng),張智,焦遠碧,等.用b值橫截距法預(yù)報強余震震級的方法研究[C]//地震預(yù)報方法實用化研究文集(地震學專輯).北京:地震出版社,1989:183-193.

[5] 張建中,宋良玉.地震b值的估計方法及其標準誤差[J].地震學報,1981,(3):292-301.

[6] 陳培善,白彤霞,李保昆.b值和地震復(fù)發(fā)周期[J].地球物理學報,2003,46(4):510-519.

[7] 謝華章.b值數(shù)字模擬的再研究[J].華北地震科學,1991,9(1):28-34.

[8] 韓渭賓.b值在地震預(yù)測中的三類應(yīng)用及其物理基礎(chǔ)與須注意的問題[J].四川地震,2003,(1):1-5.

[9] 邊慶凱,李守忠.張北6.2級地震前地震活動性異常及序列特征研究[J].華北地震科學,1999,17(3):35-42.

[10] 陳紹緒,李淑蓮.2000年6月25日唐山4.4級地震預(yù)測實踐[J].華北地震科學,2000,18(4):53-57.

Statistical Analysis of B-values in Shandong Area

XUE Jing-nan,LI Xiang,ZHANGBo,WANG Yong-wang
(School of Mathematics&System Science Shandong University Jinan,Shandong 250014)

Based on small earthquakes data,this paper calculates and analyzes b-values from 1996Jan.to 2009May.in Shandong area using the least square method and maximum likelihood method and compares these two methods based on the results.Finally,we discuss the characteristics of b-values temporal variation in Shandong area and problems reflected by it.

b-value;least square method;maximum likelihood method;seismic activity;statistical analysis

P315.5

A

1003-1375(2011)01-0001-05

2010-09-01

薛靖楠(1990-),男(漢族),山東濟南人,山東大學數(shù)學與系統(tǒng)科學學院在讀學生,主要從事概率統(tǒng)計研究.E-mail:xjn1990@126.com.