基于時間相依地震活動性模型的地震時間序列模擬

徐偉進， 吳健， 李雪婧， 高孟潭

1 中國地震局地球物理研究所， 北京 100081 2 中國地震災害防御中心， 北京 100029

0 引言

傳統(tǒng)概率地震危險性分析方法主要是基于震源模型、地震活動性模型、地震動預測模型以及場地條件等來計算地震動影響場，結果所表示的是所有震源對場址地震危險性的綜合影響.近些年來，隨著計算能力的提升，在地震危險性分析、地震災害損失預測以及地震巨災保險模型開發(fā)等領域，科學家和工程師們也希望采用一系列單個地震事件計算目標場址或區(qū)域的地震動影響(Ebel and Kafka, 1999; Musson, 2000; Hong and Goda, 2006; Assatourians and Atkinson, 2013; Faenza et al., 2017; Xu et al., 2021).這就需要將傳統(tǒng)潛在震源轉換成一系列的地震事件，而地震活動性模型是生成地震事件的重要基礎.按照地震發(fā)生概率(指未來某一確定時間段內的地震發(fā)生概率)是否隨時間變化，地震活動性模型分為兩類，即時間獨立(time-independent)模型和時間相依(time-dependent)模型.時間獨立的地震活動性模型假設地震發(fā)生概率不隨時間變化，地震的發(fā)生在時間上服從指數分布，也稱為泊松模型.在傳統(tǒng)概率地震危險性分析中，即假設地震發(fā)生在時間上是一個泊松過程(Cornell,1968; 高孟潭，1996，2015；潘華等，2013；Petersen et al., 2014).Xu等(2021)基于我國最新的地震活動性模型模擬生成了符合我國時空強分布特征的隨機地震事件集，該事件集在時間上符合泊松模型.

然而，近幾十年來，也有許多研究表明，一些斷層上的地震活動(地震發(fā)生概率)隨時間有明顯的變化，呈現時間相依特性，如果使用泊松模型，可能會高估或低估某一個時間段內的地震危險性水平.研究者們使用不同地區(qū)的地震目錄對地震的時間特性進行了實證研究，發(fā)現在許多情況下地震并不符合泊松模型，而時間相依模型能夠更好的描述地震的時間分布特征(Utsu，1984；Nishenko and Buland，1987; Ogata and Abe，1991；Tripathi，2006).研究者們提出了諸如采用伽馬(Gamma)模型、對數正態(tài)(Lognormal)模型、威布爾(Weibull)模型以及布朗過程時間(Brownian Passage Time，BPT)模型來描述地震的時間分布特征(Utsu，1984；Matthews et al.，2002；Tripathi，2006；Field and Jordan, 2015；Pasari and Dikshit，2015，2018).這些模型為地震預測和地震危險性分析提供了重要的理論支撐.國內外許多學者和機構開展了基于時間相依地震活動性模型的地震發(fā)生概率和地震危險性分析(WGCEP, 1988, 1990,1995, 2003, 2007；聞學澤,1998; 楊明和劉百篪，2000；冉洪流，2006; Field et al.，2015).

由于大地震時間相依的地震活動性特征，地震發(fā)生概率是隨時間變化的，未來某段時間內的地震等效發(fā)生率也會發(fā)生變化，若仍采用泊松模型進行隨機地震事件的模擬，那么在足夠長的時間內，得到的地震事件數量會多于或少于實際的地震事件數.由于大地震的危險性更大，這將對地震危險性分析結果產生顯著影響.因此，在具有時間相依地震活動特征的斷層(震源)上，應采用時間相依的地震活動性模型進行地震事件集的模擬，使模擬地震的時間序列具有時間相依特性.

具有時間相依特征的地震事件的模擬，可通過以下兩種方案進行：

(1)基于地震復發(fā)間隔、復發(fā)間隔的不確定性、最近地震的離逝時間等地震活動性參數，采用時間相依的地震活動性模型計算未來某段時間內的地震發(fā)生概率，然后轉換成等效地震年發(fā)生率.基于等效年發(fā)生率，仍采用泊松模型模擬地震的時間序列.該方案是將地震發(fā)生率進行了重新計算，但是地震事件的時間序列仍近似為泊松分布.目前時間相依的地震危險性分析大都是基于這一方案進行的(Petersen et al., 2007; Hebden and Stein, 2009; Field et al.，2015).

(2)基于地震復發(fā)間隔和復發(fā)間隔的標準差以及地震離逝時間，采用基于布朗隨機過程的地震加卸載過程，生成具有準周期性的地震時間序列.

模擬具有時間相依特征的地震時間序列需要用到時間相依的地震活動性模型及其相關的發(fā)生概率計算，以及描述地震準周期發(fā)生的運動概率模型，將在下文中詳細介紹.

1 時間相依的地震活動性模型

特征地震的發(fā)生具有準周期性，但地震發(fā)生時刻受諸多不確定性因素的影響，不會按照某個固定時間間隔精確地發(fā)生.實際情況中，特征地震的復發(fā)間隔往往具有很大的不確定性.對特征地震的描述主要集中在震級-頻度分布模型和時間分布模型上.由于特征地震的震級比較固定，震級范圍在0.5個震級單位內.因此，建立合適的數學模型描述特征地震的復發(fā)間隔是時間相依地震危險性分析的關鍵環(huán)節(jié).目前使用較多的描述地震復發(fā)間隔的模型主要有對數正態(tài)分布模型、正態(tài)分布模型、布朗過程時間模型、伽馬模型以及威布爾模型等.美國加州地震預報工作組(Working Groups on California Earthquake Probabilities，WGCEP)在1988、1990和1995年的報告中(WGCEP, 1988, 1990,1995)以及Petersen等(2007)均采用了對數正態(tài)分布模型(lognormal distribution).WGCEP在2003和2007年的報告中(WGCEP, 2003, 2007)以及Field等(2015)采用了布朗過程時間(BPT)模型.下面對BPT模型和對數正態(tài)分布模型做簡要描述.

1.1 布朗過程時間(BPT)模型

Matthews等(2002)基于Reid(1910)的彈性回跳理論，并結合多個其他物理參數，提出了地震復發(fā)的BPT模型.該模型假設發(fā)震斷層以恒定的加載速率進行加載，采用標準布朗運動來描述累積加載過程.在實際使用中，模型的均值和標準偏差可直接采用特征地震序列計算得到.BPT模型在理論上反映了地震孕育和發(fā)生的內在物理機制，是近年來地震學家們廣泛使用的模型.BPT模型的概率密度函數可表示為：

(1)

式中，μ為復發(fā)間隔的均值，α=σ/μ為復發(fā)間隔的變異系數，σ為復發(fā)間隔的標準差.

BPT模型的累積概率函數為：

=Φ[u1(t)]+e2/α2Φ[-u2(t)],

(2)

其中

為標準正態(tài)分布函數的累積概率函數.圖1是μ為1，α分別為0.25、0.5、1、2時BPT模型的概率密度函數和累積概率函數，可以看出α對概率密度函數的形狀影響顯著，在α較小時，概率密度函數幾乎以均值對稱分布，更接近于正態(tài)分布.在α較大時，概率密度函數曲線峰值顯著向左移，與指數分布相似.

圖1 BPT模型的概率密度函數(a)和累積概率函數(b)，BPT(1，α)，α=0.25, 0.5, 1, 2Fig.1 Example probability distributions for the BPT model with a mean recurrence interval (μ) of 1 and aperiodicity (α) values of 0.25, 0.5, 1 and 2(a) The probability density function (PDF) for recurrence interval; (b) The associated cumulative distributions.

1.2 對數正態(tài)分布模型

對數正態(tài)分布模型是將變量的對數值看成是正態(tài)分布，其概率密度函數為：

(3)

其中μ,σ分別為變量對數的均值和標準差.在實際計算中，我們已知的是變量的期望值E(T)和方差var(T)，那么變量對數的均值μ和方差σ2可分別用(4)(5)式計算:

(4)

(5)

圖2中紅色虛線為對數正態(tài)分布的概率密度函數和累積分布函數.通過與BPT的分布(圖中藍線)比較，可以看出對數正態(tài)分布和BPT分布非常相似.

圖2 對數正態(tài)分布、BPT分布和指數分布的概率密度函數(a)和累積分布曲線(b)，各分布的均值為1，標準偏差為0.5(指數分布除外)Fig.2 Probability density (a) and cumulative distribution (b) functions of exponential, BPT, and lognormal models. All distributions have mean 1 and standard deviation 0.5 (except the exponential distribution)

2 時間相依模型的地震發(fā)生概率

2.1 時間相依的地震發(fā)生概率

根據上述介紹的時間相依的地震活動性模型，在已知大地震的離逝時間和大地震復發(fā)間隔及其不確定性時，可計算未來一段時間內地震發(fā)生的概率.設Te(escape time)為最近一次地震發(fā)生距今的時間，即所謂的離逝時間，那么未來ΔT時間內發(fā)生至少一次地震的條件概率為(Matthews et al., 2002)：

(6)

圖3為BPT模型中，地震復發(fā)間隔為μ=3000 a，α分別為0.3、0.5和0.7時，未來50年發(fā)生地震的條件概率隨離逝時間的變化曲線.可以看出隨著離逝時間的增加，地震發(fā)生的概率先是緩慢增大，然后急劇升高，最后歸于平緩.α對地震發(fā)生的條件概率影響顯著，在離逝時間較長時，α越小概率越大，相反，在離逝時間較短時，α越大則概率越大.

圖3 未來50年地震發(fā)生概率隨離逝時間變化的函數曲線，BPT(3000a，α)，α=0.3, 0.5, 0.7Fig.3 50-year probability as a function of elapsed time for BPT model with a mean recurrence interval (μ) of 3000 a and aperiodicity (α) values of 0.3, 0.5 and 0.7

圖4為地震復發(fā)間隔為μ=3000 a，離逝時間為1000 a時，采用BPT模型計算的在未來某一時間段內發(fā)生至少一次地震的條件概率.

圖4 未來一段時間內地震發(fā)生的BPT條件概率，BPT(3000 a，α)，α=0.3, 0.5, 0.7，T=1000 aFig.4 Probability of one or more events for the case in which the elapsed time is 1000 a for BPT model with a mean recurrence interval (μ) of 3000 a and aperiodicity (α) values of 0.3, 0.5 and 0.7

在地震學研究中，研究者們還關心某一斷層上自上次發(fā)生地震以來還無地震發(fā)生的年發(fā)生概率，稱為危險率(hazard rate):

(7)

圖5為采用指數模型、BPT模型和對數正態(tài)模型計算的危險率隨時間的變化曲線，各模型變量的均值為1，標準偏差為0.5.從圖中可以看出對于指數模型，危險率不隨時間變化.時間相依的危險率隨離逝時間的增加有起伏變化，與泊松模型的危險率顯著不同.還可看出BPT模型和對數正態(tài)模型的危險率曲線非常相似，僅在離逝時間相對很大時才有較為顯著的差異.考慮到BPT模型的物理意義及其使用的廣泛性，下文僅討論基于BPT模型的地震發(fā)生概率.

圖5 指數模型、對數正態(tài)分布模型和BPT模型的危險率隨時間變化函數，各分布模型均值為1，標準差為0.5(指數分布除外)Fig.5 Hazard-rate functions for exponential, lognormal, and BPT models. All distributions have mean 1 and standard deviation 0.5 (except the exponential distribution)

2.2 離逝時間未知的時間相依地震發(fā)生概率計算

離逝時間是計算時間相依地震發(fā)生概率的重要參數.許多情況下，一條斷層上最近一次地震的發(fā)生時間無法確定，因此計算不出具體的離逝時間.這時候需要考慮離逝時間的所有可能性.根據Field和Jordan(2015)的研究結果，離逝時間的概率分布可以用(8)式表示：

(8)

從式中可以看出離逝時間的概率密度函數為1減去地震復發(fā)間隔的累積分布函數除以平均復發(fā)間隔.

已知離逝時間的概率密度函數，結合地震發(fā)生條件概率的計算公式，基于全概率理論，對離逝時間的所有可能性進行積分，積分函數為離逝時間的概率分布乘以條件概率計算公式，具體計算公式為：

(9)

圖6為離逝時間未知，采用BPT模型計算的未來一段時間地震發(fā)生的條件概率，地震復發(fā)間隔的均值μ=3000 a，α分別為0.3、0.5和0.7.α越大，采用BPT模型計算的數值越接近泊松模型，這是由于α越大，地震復發(fā)間隔的不確定性越強，越接近于泊松分布.

圖6 離逝時間未知時采用BPT模型計算的未來一段時間至少發(fā)生一次地震的條件概率Fig.6 BPT conditional probability of one or more events for the case in which the date of the last event is unknown

圖7為地震復發(fā)間隔的均值μ=3000 a，α為0.5時，采用BPT模型計算的離逝時間已知和未知情況下未來一段時間發(fā)生至少一次地震的條件概率.從圖中可以看出由于離逝時間未知情況下是考慮了離逝時間的所有情況，因此計算的條件概率要大于離逝時間較小時的值，而小于離逝時間較大時的值.

圖7 BPT模型在離逝時間已知和未知情況下至少發(fā)生一次地震的條件概率的比較Fig.7 Comparison of BPT conditional probability of one or more events for the case in which the date of the last event is known and unknown

上述分析了在離逝時間未知情況下地震發(fā)生條件概率的計算，考慮了離逝時間所有的可能性.然而，在一些情況下，地震學家雖然未知斷層上最近一次地震發(fā)生的具體時間，但卻可以確定過去某一時刻到現在無地震發(fā)生，稱之為歷史開放間隔(historic open interval).這就縮小了離逝時間的范圍.離逝時間的概率分布可由(10)式表示：

(10)

式中TH為歷史開放時間間隔.那么未來一段時間至少有一次地震發(fā)生的條件概率可寫為：

(11)

圖8為已知歷史開放時間間隔，采用BPT模型計算的條件概率相對于離逝時間完全未知的計算值的比值，稱為概率增益.可以看出ΔT越小，概率增益越大，這是由于相對于離逝時間完全未知的情況，已知歷史開放時間間隔時，斷層更接近于失效狀態(tài)，計算的條件概率相對于離逝時間未知情況下更大.

圖8 考慮歷史開放時間和離逝時間完全未知情況下BPT模型計算未來一段時間地震發(fā)生的條件概率的比值隨歷史開放時間的變化曲線圖中不同曲線代表其所標注預測時間的結果，(a),(b)和(c)分別對應變異系數為0.3,0.5和0.7.Fig.8 BPT-model conditional probability of one or more events for the case in which the historic open interval is accounted for, divided by the probabilities for when the elapsed time is unknown, plotted against the open intervalThe alternative curves represent results for different forecast durations (ΔT) as labeled. (a), (b) and (c) are for aperiodicity (α) values of 0.3, 0.5 and 0.7, respectively.

2.3 等效地震發(fā)生率的計算

為了進行地震危險性計算，還需要知道特征地震在未來ΔT內的年發(fā)生率r.根據上文中的公式(6)、(9)和(11)可計算出特征地震在ΔT時段內的發(fā)生概率P.r仍然可看作服從泊松分布(Petersen et al., 2007)，那么P=1-exp(-r·ΔT)，則等效發(fā)生率為：

r=-ln(1-P)/ΔT.

(12)

圖9為根據BPT模型計算的未來一段時間內的等效發(fā)生率.其中μ=3000 a，α=0.5，離逝時間分別為1000 a、2000 a和3000 a.可以看出在離逝時間較短，ΔT較小時，BPT模型的地震發(fā)生率小于泊松模型，這是由于在離逝時間較短時，斷層上構造應力累積較小，發(fā)生地震的概率也較小.根據BPT模型的物理意義，某一特定斷層上，未來一段時間發(fā)生地震的概率與上一次發(fā)生地震的離逝時間相關，離逝時間越久，則表示構造應力累積越大，發(fā)生地震的概率就越大，因此計算的等效發(fā)生率也越大.

圖9 根據BPT模型和指數模型計算的未來一段時間的等效發(fā)生率隨著ΔT的變化曲線Fig.9 Seismic rate as a function of forecast durations (ΔT) for BPT model and Poisson model

3 基于布朗隨機過程的地震加卸載過程

統(tǒng)計模型廣泛應用于科學研究中，通過對事件時間序列的統(tǒng)計分析，可獲知事件發(fā)生的時間特征.在地震學研究中，研究者們采用指數分布、對數正態(tài)分布、威布爾分布以及伽馬分布等描述地震事件的時間序列過程.通過對這些模型進行經驗擬合，研究者們可獲知地震的時間分布特征.然而這些統(tǒng)計模型不能對地震發(fā)生的時間過程進行深刻的描述，存在物理意義不明確的問題(Matthews et al., 2002).Matthews等(2002)基于Reid(1910)的彈性回跳理論，并考慮構造應力、介質物性等物理量，提出了采用布朗運動來描述地震震源的加載過程.

設Y(t)為t時刻的斷層震源加載狀態(tài)，震源的加卸載過程可用(13)式表示：

Y(t)=x0+X(t)-X[R(t)],

(13)

其中x0為震源加載的初始狀態(tài)，即Y(0)=x0；X(t)為從時間0到時間t的累積加載，X(t)=λt，λ為加載率；X[R(t)]為最近一次斷層失效(地震發(fā)生)時的震源的累積加載，在斷層失效后，震源立即回到初始加載狀態(tài).

將布朗隨機擾動加入累積加載便可得到加載狀態(tài)的隨機過程：

X(t)=λt+σW(t),t≥0,

(14)

其中W為標準布朗運動，σ為隨機擾動的幅度參數.標準布朗運動可通過對平穩(wěn)增量進行積分得到，該平穩(wěn)增量的分布為均值為零方差為常數的正態(tài)分布.

斷層失效狀態(tài)可寫為：

xf=x0+δ,δ>0,

(15)

其中δ為可接受的最高累積加載.實際情況下失效狀態(tài)下的加載是已知的.

根據上述描述，布朗更新過程可用平均加載率λ、隨機擾動率σ、加載初始狀態(tài)x0以及失效狀態(tài)xf等四個參數來確定.布朗運動和布朗加卸載過程可分別用BM(x0,λ,σ)和BRO(x0,λ,σ,xf)來表示.圖10顯示了三種布朗加卸載過程：BRO(0,1,σ,1)，σ取值為0.25、0.5和1.0.可以看出在σ較小時，地震復發(fā)呈現明顯的周期性.當σ較大時，地震發(fā)生的隨機性顯著增強.

圖10 λ=xf= 1時布朗加卸載過程路徑狀態(tài)曲線(a) σ=0.25, (b) σ=0.5, (c) σ=1.0.Fig.10 Load-state paths for a Brownian relaxation oscillator with λ=xf=1

地震的發(fā)生率(復發(fā)間隔的倒數)即是加載過程的加載率，復發(fā)間隔的不確定性是加載過程的隨機擾動.離逝時間是加載時間.由此可預測未來一段時間斷層應力的加載狀態(tài).上述中失效狀態(tài)的時刻即是地震發(fā)生的時間，記錄這些時間便得到符合時間相依分布的地震時間序列.由此可見，根據第2節(jié)計算的地震發(fā)生率(復發(fā)間隔)、復發(fā)間隔的不確定性以及地震離逝時間等參數，可采用基于布朗隨機過程的地震加卸載過程來模擬地震的時間序列，這是基于時間相依地震事件模擬的最重要環(huán)節(jié).

4 基于BPT模型的特征地震時間序列模擬實例

上文中系統(tǒng)介紹了時間相依的地震活動性模型、時間相依地震發(fā)生概率的計算方法以及基于布朗隨機過程的地震加載過程理論，基于這些理論、模型和方法，我們模擬了符合時間相依準周期分布特征的地震序列.假設某一斷層上地震復發(fā)間隔的均值為μ=3000 a，標準差α=0.25，以此參數模擬了時間相依的地震序列，并與具有相同μ值的泊松分布的地震序列進行了比較(圖11).可以發(fā)現基于BPT模型模擬的地震序列呈現準周期分布的特征，而泊松地震序列則呈現出隨機分布特征.

圖11 基于BPT模型和泊松模型的模擬隨機地震序列Fig.11 Simulation of random earthquake sequence based on BPT and Poisson models

基于BPT模型模擬的地震時間分布受標準差的影響，為了了解標準差對模擬地震序列時間分布特征的影響，我們模擬了不同標準差下(α=0.1，α=0.3，α=0.5，α=0.7)的BPT隨機地震序列(圖12)，可以看出，在標準差較小時，地震序列呈準周期分布，當標準差較大時，地震序列呈現隨機性，這與上述BPT模型概率密度函數的分析是一致的.

圖12 在平均復發(fā)間隔為3000 a，變異系數分別為0.1、0.3、0.5和0.7時，基于BPT模型的隨機模擬地震時間序列Fig.12 Simulated earthquake time series based on BPT model with a mean recurrence interval (μ) of 3000 a and aperiodicity (α) values of 0.1, 0.3, 0.5 and 0.7

上文完成了地震時間序列的實例模擬，模擬的地震時間序列能夠體現時間相依特征，與理論模型是吻合的.一套地震序列，除了發(fā)震時間外，還有地震震級、發(fā)震位置等參數.對于震級，主要根據震級-頻度關系(G-R關系)來模擬，也可根據歷史地震目錄中地震的震級隨機抽取.關于地震位置，可依據地震的空間分布概率來確定，對于單條斷層的特征地震，可根據斷層幾何形狀、產狀以及斷層面上的物理參數來確定.具體方法可參考Xu等(2021).通過上述方法，可模擬一系列地震的時間、地震位置以及震級等參數，生成地震事件集，可用于地震預測、地震危險性分析以及地震災害評估等領域.

5 結論與討論

為了模擬具有時間相依特性的地震時間序列，本文系統(tǒng)研究了描述時間相依地震活動特征的理論、模型與方法.根據時間相依地震活動性模型及其參數特征，給出了已知和未知大地震離逝時間的情況下地震發(fā)生概率的計算流程以及等效發(fā)生率的計算方法.基于BPT加卸載過程，建立了大地震準周期復發(fā)模擬方法.以特征地震模型為例，進行了特征地震時間序列的模擬.與泊松時間序列相比，基于BPT模型的時間序列具有準周期性.此外，模擬特征地震時間序列受地震復發(fā)間隔標準差的影響，標準差越小，序列越趨于周期性，反之，時間序列則呈現隨機性.

本文中，基于時間相依地震活動性模型模擬的地震事件與基于泊松模型的相比，主要有兩個方面的差異，一是地震發(fā)生率的區(qū)別，對于泊松模型，地震年發(fā)生率是不隨時間變化的，而時間相依地震活動性模型則相反，并且在地震離逝時間較長的情況下，基于時間相依的地震活動性模型計算的地震發(fā)生率要顯著高于泊松模型，這對地震危險性分析具有重要影響.二者的第二個差異是在時間軸上的表現不同，本研究中模擬的時間相依的特征地震時間序列呈現周期性，而泊松模型時間序列是完全隨機的.

本文模擬地震時間序列基于的是BPT模型，該模型的理論基礎是彈性回跳學說，彈性回跳理論認為，在斷層上發(fā)生一次特征地震后，構造應力急劇釋放，后續(xù)短期內發(fā)生大地震的概率急劇減小，隨著時間推移，構造應力重新累積，發(fā)震概率也隨之增大.該理論目前被廣泛用于地震預測和地震危險性分析.但是，實際情況下，由于一條斷層上特征地震記錄非常少，很難驗證該理論的正確性.以后隨著地震記錄的增多，在一條斷層上有多個特征地震發(fā)生，可對該理論進行驗證.

本文介紹的方法對于提高地震概率預測水平具有重要意義.在任一斷層上，通過本文介紹方法，可模擬一系列具有時間相依特征的地震事件集，并用于地震預測、地震危險性分析以及地震災害評估等領域.本文研究成果已應用于中國地震巨災保險模型的構建.