青藏高原東北緣地震活動性廣義帕累托模型的全域敏感性分析＊

任夢依劉哲

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

2） 中國河北唐山 063000 唐山地震監(jiān)測中心站

引言

地震危險性估計中，通常用地震活動性模型描述潛源區(qū)地震的時空分布、強(qiáng)度分布、頻度分布等地震活動特征（蔣溥，戴麗思，1993）.已有的分析方法可分為確定性方法和概率性方法兩個類別，概率性方法是目前構(gòu)建地震活動性模型中經(jīng)常使用的分析方法，該方法主要基于研究區(qū)的歷史地震活動性構(gòu)建模型，其中最常用的是利用經(jīng)典的古登堡-里克特的震級頻度關(guān)系（G-R關(guān)系）中對數(shù)線性關(guān)系的外推.然而對于震級上限而言，傳統(tǒng)的G-R關(guān)系模型因為線性外推震級上限無限大，與實(shí)際情況不符.任雪梅等（2012a）對我國分區(qū)域進(jìn)行了G-R關(guān)系的修正和震級極限值的確定，結(jié)果表明我國大陸存在非線性的震級-頻度關(guān)系，修正的G-R關(guān)系模型運(yùn)用于我國中強(qiáng)以上地震震級-頻度關(guān)系的擬合更為適合.

利用極值理論分析地震活動性也屬于概率性方法中常用的一種.近年來，一些學(xué)者探討了依據(jù)廣義極值理論基于廣義帕累托分布（generalized Pareto distribution，縮寫為GPD）構(gòu)建地震活動性模型的方法，用以研究震級分布的尾部特征.Pisarenko和Sornette （2003）分析了哈佛CMT目錄中1977——2000年12個地震帶的淺層地震，發(fā)現(xiàn)廣義帕累托分布可以較好地擬合震級分布的尾部特征；錢小仕等（2013a，b）利用廣義帕累托分布分析了云南地區(qū)歷史地震資料以及中國大陸各活動地塊邊界帶的強(qiáng)震震級分布特征，提出了一種利用強(qiáng)震數(shù)據(jù)推斷最大震級分布的可能途徑；田建偉等（2017）利用廣義帕累托分布估計了馬尼拉海溝俯沖帶潛在地震海嘯源區(qū)的地震危險性；任夢依（2018）利用地震活動性廣義帕累托模型估計了龍門山地區(qū)的震級上限；Dutfoy （2019）將廣義帕累托分布與泊松分布相結(jié)合，分析了法國南部地區(qū)地震震級分布的尾部特征并計算了年最大地震震級.

構(gòu)建地震活動性廣義帕累托模型，利用的是歷史地震記錄中大震級區(qū)段的記錄數(shù)據(jù)，假設(shè)大震級區(qū)段的震級分布符合廣義帕累托分布.這種方法適合歷史地震記錄時間長、中小地震記錄缺失的地區(qū).然而，構(gòu)建地震活動性廣義帕累托模型，需要人為選取歷史地震記錄起始時間和震級閾值，這些模型輸入?yún)⒘窟x取的不確定性可能導(dǎo)致地震危險性估計結(jié)果的不確定性.

對地震危險性估計的不確定性分析，一直為國內(nèi)外地震學(xué)家高度重視，也開展了一些相關(guān)研究.McGuire和 Shedlock （1981）及 McGuire （1987）是較早對地震危險性分析中模型輸入?yún)?shù)和模型輸出的不確定性進(jìn)行研究的學(xué)者之一，在對美國中部及東部地區(qū)進(jìn)行概率地震危險性評價時，就知識不完備導(dǎo)致的認(rèn)識不確定性，采用了邏輯樹的方法進(jìn)行表達(dá)和處理；對于地震危險性分析模型的輸入?yún)?shù)的不確定性滿足連續(xù)分布（如正態(tài)分布）的情況，Lawrence Livermore National Laboratory （1989）及 Cramer等（1996）提出了一種蒙特卡洛抽樣結(jié)合邏輯樹的不確定性表示方法；胡聿賢和鹿林（1990）在關(guān)于華北地區(qū)地震危險性研究中，分析了潛源區(qū)劃分的不確定性及地震活動性參數(shù)不確定性對地震危險性分析結(jié)果的影響；王健和高孟潭（1996）強(qiáng)調(diào)，除了要研究單個參數(shù)的敏感性外，還必須弄清參數(shù)間的交互作用以及多個參數(shù)同時變化時結(jié)果的范圍及其分布，他們利用因子設(shè)計中的正交設(shè)計法，分析了地震年平均發(fā)生率、b值、地震帶震級上限、潛在震源區(qū)震級上限與地震空間分布函數(shù)等幾個參數(shù)之間的交互作用，以及各參數(shù)取不同權(quán)重時對烈度概率分布的影響.美國公布的國家地震危險性圖，充分考慮了震級上限的不確定性，并將其列為2008年國家地震危險性圖編制工作的重大改進(jìn)（Petersenet al，2008）.我國第五代地震區(qū)劃圖在構(gòu)建地震活動性模型時，也強(qiáng)調(diào)了需充分認(rèn)識地震活動性的不確定性，并加以綜合分析和處理（潘華，李金臣，2016）.對于近些年應(yīng)用的廣義帕累托分布模型，Dutfoy （2019）利用自舉法（bootstrap）得到了地震活動性廣義帕累托模型參數(shù)的分布特征以及年度最大震級分布特征，分析討論了形狀參數(shù)概率密度分布的雙峰特征導(dǎo)致地震危險性估計結(jié)果的不確定性.

然而，已有的相關(guān)研究大都僅應(yīng)用局域敏感性分析方法（local sensitivity analysis method），只能分析當(dāng)單一輸入?yún)⒘孔兓溆鄥⒘勘３植蛔儠r，基于地震活動性模型所得到的地震危險性估計結(jié)果的不確定性，并未系統(tǒng)地分析輸入?yún)?shù)同時變化以及輸入?yún)?shù)之間存在相互作用時，地震危險性估計結(jié)果的不確定性，即缺少全域敏感性分析（global sensitivity analysis method）的研究.

因此，本文提出基于全域敏感性分析的地震危險性估計的不確定性分析的流程和方法.以青藏高原東北緣為示例研究區(qū)，應(yīng)用廣義帕累托模型，估計研究區(qū)的強(qiáng)震重現(xiàn)水平、震級上限及其相應(yīng)的置信區(qū)間，并將全域敏感性分析方法引入地震危險性估計的不確定性分析，利用拓展傅里葉幅度敏感性檢驗法（extended fourier amplitude sensitivity test，縮寫為 E-FAST）開展全域敏感性的定量分析.

1 方法

1.1 地震活動性廣義帕累托模型構(gòu)建

設(shè)X1，X2， ···，Xn是一個獨(dú)立且同分布的隨機(jī)變量序列，具有邊際分布函數(shù)，令Mn＝max{X1，···，Xn}，若存在{an＞0，bn∈R}和非退化分布函數(shù)G（x），使

則稱G為極值分布.Fisher和Tippett （1928）證明了極值分布的三大類型定理，指出G必屬于耿貝爾（Gumbel）分布、弗雷歇（Fréchet）分布、威布爾（Weibull）分布三種類型之一.廣義極值理論將這三種分布類型統(tǒng)一為一個分布類型，由位置參數(shù)μ，尺度參數(shù)σ，形狀參數(shù)ξ三個參數(shù)來表示（Coles，2001；史道濟(jì)，2006），即

對于足夠大的u，當(dāng)X＞u，y＝（X－u）的條件分布函數(shù)可以近似服從廣義帕累托分布，即

在構(gòu)建地震活動性廣義帕累托模型時，首先需確定震級閾值u.閾值u的選取通過兩種方法：① 依據(jù)震級平均超出量分布函數(shù)圖，震級平均超出量與震級閾值應(yīng)滿足線性關(guān)系；② 隨著閾值不斷增大，形狀參數(shù)和修飾的尺度參數(shù)估計值應(yīng)基本保持穩(wěn)定（Coles，2001）.

N年重現(xiàn)水平為預(yù)期每N年至少會發(fā)生一次震級大于xN的地震事件，如果每年地震發(fā)生次數(shù)為ny，當(dāng)ξ≠0時，N年重現(xiàn)水平為

式中，u為震級閾值，ζu為超過震級閾值的地震樣本數(shù)與地震樣本總數(shù)的比值.

考慮到潛在震源區(qū)的地震事件的震級必為有限值，因此只有當(dāng)形狀參數(shù)ξ＜0，廣義帕累托分布具有有限上界，才符合震級存在上限的條件，震級上限估計值可表示為

1.2 定量的全域敏感性分析方法——E-FAST法

E-FAST 法是 Saltelli等（1999）結(jié)合 Sobol法（Sobol，1993）的優(yōu)點(diǎn)，對 FAST （Fourier amplitude sensitivity test）法（Cukieret al，1973，1978）進(jìn)行改進(jìn)后得到的一種基于方差的全域敏感性定量分析方法（洪明理等，2014；宋明丹等，2014）.該方法將模型的敏感性概括為單個輸入?yún)?shù)的敏感性及參數(shù)間交互作用的敏感性.模型輸出y的總方差V（y）可以分解為

式中，Vi＝V［E（y|xi） ］ ，Vi，j＝V［E（y|xi，xj） ］ －Vi－Vj，Vi，j，k~V1，2，···，n，以此類推.

定義xi的一階敏感性指標(biāo)或主效應(yīng)指標(biāo)為

式中，Sxi表示單個參數(shù)獨(dú)立作用的敏感性，依據(jù)Sxi的大小可以對輸入?yún)?shù)獨(dú)立作用的敏感性進(jìn)行排序.

定義xi，xj的二階敏感性指標(biāo)為

式中，Sxi，xj表示輸入?yún)?shù)xi與xj之間的交互作用對模型輸出結(jié)果不確定性的影響，更高階的敏感性指標(biāo)的含義以此類推.

定義xi的全效應(yīng)指標(biāo)為

E-FAST方法的優(yōu)點(diǎn)是：① 適用于各類模型，對于模型本身的性質(zhì)沒有特殊的限定，可用于分析多個輸入?yún)?shù)同時變化時，模型輸出結(jié)果的不確定性；② 可用于定量地分析各個輸入?yún)⒘康牟淮_定性，以及各個輸入?yún)⒘坎淮_定性的相互作用對模型輸出結(jié)果不確定性的影響.

2 案例研究區(qū)——青藏高原東北緣

2.1 地震地質(zhì)背景

青藏高原東北緣位于青藏高原與華北克拉通的交會地帶，西南一側(cè)是不斷向北東推擠并不斷增厚的青藏高原，北東一側(cè)是相對穩(wěn)定的鄂爾多斯地塊，北西一側(cè)是相對穩(wěn)定的阿拉善地塊.在印度板塊的推動作用下，青藏高原整體向北東一側(cè)運(yùn)動，東北緣地區(qū)就處在北東向擴(kuò)展的最前沿.

青藏高原東北緣的構(gòu)造多形成于印度板塊與歐亞大陸碰撞引起的SW-NE向擠壓，該擠壓始于 50 Ma前（Yin，Harrison，2000；Taylor，Yin，2009）.全球定位系統(tǒng)（GPS）的測量結(jié)果顯示，青藏高原東北緣約以每年10 mm的速率沿NE-SW向或NNE-SSW向水平縮短（Ganet al，2007）.長期以來，青藏高原東北緣構(gòu)造變形活動強(qiáng)烈、大震多發(fā)（Liet al，2017）.該地區(qū)分布著幾組具有不同活動性質(zhì)的斷裂帶，包括左旋走滑斷裂帶，右旋走滑斷裂帶，端部逆沖斷裂帶（徐化超，2018）.這些大型活動斷裂對該區(qū)地殼穩(wěn)定性起著重要的調(diào)節(jié)作用，而在彼此的共同作用下，該區(qū)新構(gòu)造運(yùn)動強(qiáng)烈，一直以來地震活動十分活躍.

歷史上青藏高原東北緣發(fā)生過M≥8.0地震5次，分別為1654年天水M8.0地震、1739年平羅M8.0地震、1879年武都M8.0地震、1920年海原M8.5地震和1927年古浪M8.0地震（圖1）.本文所用地震目錄主要源自國家地震科學(xué)數(shù)據(jù)共享中心的中國歷史地震目錄（顧功敘，1983）、國家地震臺網(wǎng)地震目錄、中國臺網(wǎng)正式地震目錄以及國家地震局震害防御司編制的中國歷史強(qiáng)震目錄（1995），共收集了青藏高原東北緣（95°E——107°E，32°N——40°N）1880——2020年M≥4.0歷史地震1 222次，其中M≥8.0地震3次，7.0≤M＜8.0地震13次，6.0≤M＜7.0地震54次，5.0≤M＜6.0地震228次.黃瑋瓊等（1994）對中國大陸大部分地區(qū)的歷史地震資料作了較詳細(xì)地分析與研究，顯示青藏高原東北緣M≥5.0地震資料在1885年之后基本完整.

圖1 青藏高原東北緣主要活動斷裂及M≥4.0地震分布圖Fig.1 Map of main active faults and M≥4.0 earthquakes distribution in northeastern Tibetan Plateau

為了使震級樣本盡量滿足構(gòu)建廣義帕累托模型時獨(dú)立同分布的條件，在進(jìn)行地震危險性估計之前需要盡量刪除地震目錄中的余震.目前國內(nèi)外地震危險性分析研究中，考慮到計余震后，共收集到青藏高原東北緣1880——2020年M≥4.0歷史地震864次，M≥5.0地震的M-t圖如圖2所示.

圖2 青藏高原東北緣 1880——2020 年M≥5.0地震的M-t圖Fig.2 M-t map of the M≥5.0 earthquakes from 1880 to 2020 in northeastern Tibetan Plateau

2.2 地震危險性估計

本文分析了1885——2020年青藏高原東北緣的地震目錄，利用廣義帕累托模型估算了青藏高原東北緣震級重現(xiàn)水平和震級上限.圖3為青藏高原東北緣震級數(shù)據(jù)的平均超出量分布圖，置信區(qū)間為95%.在平均超出量函數(shù)具有線性特征的區(qū)段選取對應(yīng)的閾值.當(dāng)u大于5.5時，平均超出量函數(shù)近似呈線性.另外，還要綜合考慮樣本數(shù)據(jù)的實(shí)際情算的簡便和工程應(yīng)用上的可操作性，普遍采用 Gardner和 Knopoff （1974）提出的時空窗法來刪除余震（徐偉進(jìn)，2012）.空間窗口的確定方法為（Knopoff，Gardner，1972）

圖3 震級平均超出量分布函數(shù)圖Fig.3 Diagram of mean excess function of magnitude

式中，R為窗口半徑，a，b為固定參數(shù)值，在我國a，b的取值分別為0.5和?1.78 （劉杰等，1996）.

最終確定的時間窗口如表1所示.刪除況，通過觀察選取不同閾值時，形狀參數(shù)和修飾的尺度參數(shù)估計值的變化情況（圖4）來進(jìn)行取舍.應(yīng)選取形狀參數(shù)和修飾的尺度參數(shù)較為穩(wěn)定呈近似水平直線的區(qū)段所對應(yīng)的閾值.由圖4容易看出，當(dāng)閾值大于5.5時，兩個參數(shù)的估計值隨閾值的增大相對比較穩(wěn)定.另外需要考慮的是，閾值取得過大，超閾值的樣本量較少，會導(dǎo)致估計量的方差過大，反之閾值取得太小，使得超出量分布與廣義帕累托分布出現(xiàn)較大偏差，會導(dǎo)致估計量的標(biāo)準(zhǔn)差過大.綜合考慮超出量分布函數(shù)的線性要求和分布參數(shù)估計的穩(wěn)定性要求，選取閾值u＝5.5.

表1 不同震級地震的余震時間窗（Gardner，Knopoff，1974）Table 1 Aftershock time windows of different magnitudes（Gardner，Knopoff，1974）

圖4 修飾的尺度參數(shù)（a）和形狀參數(shù)（b）估計值隨震級閾值選取的變化Fig.4 Variation diagram of estimators for modified scale parameter （a） and shape parameter （b）according to the magnitude threshold

地震活動性廣義帕累托模型，可以充分利用歷史地震記錄中大震級區(qū)段的信息，為研究震級分布的尾部特征提供統(tǒng)計分析方法.本文模型設(shè)置震級閾值為5.5，因而，如果所采用的地震目錄小于M5.5的地震記錄缺失，將不會對模型結(jié)果產(chǎn)生影響，該模型的特點(diǎn)也正是適合歷史地震記錄時間長、中小地震記錄缺失的地區(qū).

確定震級閾值后，得到研究區(qū)震級的超出量樣本.依據(jù)式（3）進(jìn)行廣義帕累托分布的擬合，利用極大似然估計方法，得到修飾的尺度參數(shù)和形狀參數(shù)的估計值分別為1.12和?0.32.由于形狀參數(shù)的估計值為負(fù)值，因此廣義帕累托分布為威布爾分布，右端點(diǎn)為有限值，意味著研究區(qū)存在震級上限.將閾值、形狀參數(shù)和修飾的尺度參數(shù)估計值分別代入式（4）和式（5）進(jìn)行研究區(qū)地震危險性估計，結(jié)果列于表2.

表2 青藏高原東北緣地震危險性估計Table 2 Seismic hazard estimation for northeastern Tibetan Plateau

對廣義帕累托分布擬合的情況進(jìn)行診斷，通常根據(jù)實(shí)際數(shù)據(jù)與分布模型的一致性來判斷模型的準(zhǔn)確性.診斷結(jié)果由4個圖（概率圖、分位數(shù)圖、重現(xiàn)水平圖、密度曲線圖）表示（圖5）.概率圖根據(jù)樣本數(shù)據(jù)的累積概率與模型分布的累積概率之間的關(guān)系繪制，分位數(shù)圖依據(jù)樣本數(shù)據(jù)分布的分位數(shù)與模型分布的分位數(shù)之間的關(guān)系繪制.當(dāng)樣本數(shù)據(jù)符合理論假設(shè)的分布時，概率圖（圖5a）和分位數(shù)圖（圖5b）應(yīng)近似為直線.圖5c代表的是對應(yīng)于不同重現(xiàn)期的重現(xiàn)水平，樣本數(shù)據(jù)應(yīng)落在給定分布分位數(shù)估計置信區(qū)間內(nèi).由于ξ＜0，為威布爾分布，符合函數(shù)存在上限值的條件，因此重現(xiàn)水平曲線為凸曲線，且逐漸趨于有限值.圖5d中擬合的密度曲線與數(shù)據(jù)的直方圖也較為一致.綜上，從圖5可以看出，各散點(diǎn)數(shù)據(jù)基本緊密圍繞各參考線分布，表明擬合狀態(tài)良好，上述診斷結(jié)果不能否定利用廣義帕累托分布分析青藏高原東北緣強(qiáng)震震級分布特征的合理性和適用性.

圖5 廣義帕累托分布擬合診斷圖（a） 概率圖；（b） 分位數(shù)圖；（c） 重現(xiàn)水平圖；（d） 概率密度圖Fig.5 Diagnostic plots of the generalized Pareto distribution fitted to magnitude（a） Probability plot；（b） Quantile plot；（c） Return level plot；（d） Probability density plot

研究區(qū)震級上限估計值為8.95，置信度95%的置信區(qū)間為 ［ 8.03，9.87 ］ ，而研究區(qū)記錄到的歷史最大地震震級為8.5.

震級上限的含義，是指“在地震帶震級-頻度關(guān)系式中，累積頻度趨于零的震級極限值”（胡聿賢，1999）.地震帶的震級上限與帶內(nèi)歷史上最大地震震級、地震帶的地震活動水平、b值以及該帶最大地震的重現(xiàn)期有關(guān)，即

式中：MT為T年內(nèi)累積頻度為1的地震震級；A＝lg（v0T），v0為M≥4.0地震年平均發(fā)生率；T為所考慮時間段.從式中可以看到，大于1，則M＞MT，也就是說，地震帶的震級上限應(yīng)大于該帶歷史最大地震震級.當(dāng)然，由于構(gòu)造規(guī)模和巖石強(qiáng)度都有極限，實(shí)際上M8.6地震已為數(shù)極少（胡聿賢，1999）.

任雪梅等（2012b）基于修正的G-R關(guān)系模型，對祁連山——六盤山地震帶1970年以來ML3.0以上地震的震級-頻度關(guān)系進(jìn)行擬合，得到最大震級為9.66.錢小仕等（2013b）基于廣義帕累托分布的超閾值分布模型，對中國大陸活動地塊邊界帶強(qiáng)震震級分布特征開展研究，以1921年為起始年，震級閾值為5.3，估計了海源——祁連帶震級上限為9.95.

對于廣義帕累托模型，由于當(dāng)ξ逐漸趨近于0時，模型估計的震級上限值可能會變得不穩(wěn)定，Pisarenko等（2008，2014）提出可以利用分位數(shù)來進(jìn)一步量化尾部特征.錢小仕等（2013a，b）也提到在工程地震危險性分析中，潛在最大震級的確定往往與特定建筑結(jié)構(gòu)的抗震設(shè)計年限相關(guān)，他們利用0.999 97高分位數(shù)作為最大震級估計，其值往往低于廣義帕累托模型公式直接計算的震級上限估計值.本文旨在提出基于全域敏感性分析的地震危險性估計不確定性分析流程和方法，震級上限估計值未采用分位數(shù)的形式.

基于現(xiàn)有的地震活動性模型，進(jìn)行地震危險性估計，都不可避免地存在不確定性.

2.3 不確定性分析

構(gòu)建地震活動性廣義帕累托模型，需要人為選取歷史地震記錄起始時間和震級閾值.這些輸入?yún)⒘窟x取的不確定性可能導(dǎo)致地震危險性估計結(jié)果的不確定性.本文采用E-FAST方法開展定量的全域敏感性分析.

地震活動性廣義帕累托模型的輸入?yún)?shù)為地震目錄起始時間ts和震級閾值u，輸出結(jié)果為不同時期的震級重現(xiàn)水平、震級上限及其對應(yīng)的置信區(qū)間.青藏高原東北緣地震目錄起始時間范圍的選取是基于對該地區(qū)地震目錄完整性的研究（黃瑋瓊等，1994）；震級閾值的選擇范圍通過上述1.1節(jié)中閾值的兩種選擇方法確定.設(shè)地震目錄起始時間范圍為1880——1890年，震級閾值的取值范圍介于5.5——5.9之間，并假設(shè)它們在其不確定性范圍內(nèi)服從均勻分布，利用E-FAST方法進(jìn)行蒙特卡洛抽樣，共采集輸入?yún)?shù)的抽樣樣本194組.

對青藏高原東北緣地震活動性廣義帕累托模型進(jìn)行參數(shù)敏感性分析.首先，根據(jù)輸入?yún)?shù)的樣本組合，計算相應(yīng)的震級上限，以及對應(yīng)一定重現(xiàn)期的震級重現(xiàn)水平；其次，將輸入?yún)?shù)和相應(yīng)的地震危險性分析結(jié)果文件導(dǎo)入敏感性分析軟件中，利用E-FAST方法，計算廣義帕累托模型地震危險性估計結(jié)果對各輸入?yún)?shù)的敏感性指標(biāo)，即ts和u的主效應(yīng)指標(biāo)Si和全效應(yīng)指標(biāo)（表3）.為了更直觀地顯示各輸入?yún)?shù)的敏感性指標(biāo)，繪制如圖6——8所示的直方圖.

圖6 地震危險性估計結(jié)果對地震目錄起始時間ts和震級閾值u的主效應(yīng)指標(biāo)Si直方圖橫軸數(shù)字代表的含義下同F(xiàn)ig.6 Histogram of the first-order effect Si of seismic hazard estimation on earthquake catalogue start time ts and magnitude threshold u The meanings of the number represent are the same below

表3 地震危險性估計結(jié)果對地震目錄起始時間ts和震級閾值u的主效應(yīng)指標(biāo)Si和全效應(yīng)指標(biāo)Table 3 Total and the first-order effects of the seismic hazard estimation on earthquake catalogue initial time ts and magnitude threshold u

表3 地震危險性估計結(jié)果對地震目錄起始時間ts和震級閾值u的主效應(yīng)指標(biāo)Si和全效應(yīng)指標(biāo)Table 3 Total and the first-order effects of the seismic hazard estimation on earthquake catalogue initial time ts and magnitude threshold u

重現(xiàn)期/a 參數(shù)M Si images/BZ_126_783_2200_820_2242.png20 u 0.742 4 0.937 8 ts 0.059 5 0.184 0 50 u 0.749 5 0.868 0 ts 0.117 8 0.178 2 100 u 0.374 1 0.678 8 ts 0.228 6 0.522 3 200 u 0.487 6 0.936 6 ts 0.024 9 0.405 1 500 u 0.568 6 0.975 0 ts 0.004 8 0.329 7∞u 0.460 0 0.958 8 ts 0.015 2 0.473 8置信度95%的置信區(qū)間上端點(diǎn)置信度95%的置信區(qū)間下端點(diǎn)Si images/BZ_126_1265_2200_1302_2242.png0.779 1 0.968 4 0.030 2 0.158 9 0.740 6 0.982 8 0.015 3 0.182 4 0.709 5 0.987 7 0.008 8 0.204 5 0.685 8 0.988 9 0.005 6 0.223 1 0.661 6 0.987 6 0.004 0 0.244 3 0.345 0 0.929 2 0.029 5 0.623 5 Si images/BZ_126_1795_2200_1832_2242.png0.681 3 0.880 5 0.112 5 0.229 4 0.402 3 0.797 9 0.138 9 0.502 4 0.570 5 0.964 7 0.011 8 0.329 7 0.600 7 0.978 9 0.004 3 0.300 2 0.609 4 0.982 2 0.003 3 0.292 1 0.386 0 0.940 2 0.024 2 0.573 0

圖6為廣義帕累托模型的地震危險性估計結(jié)果對兩個輸入?yún)?shù)的主效應(yīng)指標(biāo)直方圖.從圖上可以看出，廣義帕累托模型估計結(jié)果明顯對u更為敏感.對于ts，除了100年震級重現(xiàn)水平對應(yīng)的主效應(yīng)指標(biāo)略高于0.2，其余的主效應(yīng)指標(biāo)均不到0.2.u的主效應(yīng)指標(biāo)遠(yuǎn)高于ts的主效應(yīng)指標(biāo)，因此在建立模型時，首先要考慮震級閾值u選取的合理性和準(zhǔn)確性.

結(jié)合表3綜合分析，對應(yīng)不同重現(xiàn)期的u的主效應(yīng)指標(biāo)有所差別，表明震級重現(xiàn)水平及震級上限對u的敏感性程度也有所不同.除了100年震級重現(xiàn)水平對應(yīng)的u的主效應(yīng)指標(biāo)為0.374 1，其余重現(xiàn)期為20，50，200和500年的震級重現(xiàn)水平以及震級上限對應(yīng)的u的主效應(yīng)指標(biāo)均高于 0.45，分別為 0.742 4，0.749 5，0.487 6，0.568 6 和 0.460 0.

圖7為廣義帕累托模型的地震危險性估計結(jié)果對兩個輸入?yún)?shù)的全效應(yīng)指標(biāo)直方圖.主效應(yīng)與全效應(yīng)的大小差異，即交互效應(yīng)，體現(xiàn)了輸入?yún)?shù)之間對輸出結(jié)果不確定性的影響存在非線性效應(yīng)與否.圖8主要體現(xiàn)了廣義帕累托模型的地震危險性估計結(jié)果對兩個輸入?yún)?shù)的交互效應(yīng)指標(biāo).圖中柱狀長度對應(yīng)的數(shù)值代表了輸入?yún)?shù)的全效應(yīng)指標(biāo)大小，其中，藍(lán)色部分代表主效應(yīng)指標(biāo)在全效應(yīng)指標(biāo)中的占比，橘色部分代表交互效應(yīng)指標(biāo)在全效應(yīng)指標(biāo)中的占比.

圖7 地震危險性估計結(jié)果對地震目錄起始時間ts和震級閾值u的全效應(yīng)指標(biāo)直方圖Fig.7 Histogram of the total effect of seismic hazard estimation on earthquake catalogue start time ts and magnitude threshold u

圖8 地震危險性估計結(jié)果對地震目錄起始時間ts和震級閾值u的主效應(yīng)指標(biāo)和交互效應(yīng)指標(biāo)直方圖Fig.8 Histogram of the first-order effect and interactions on earthquake catalogue start time ts and magnitude threshold u

從圖7和8可以看出，兩個參數(shù)對地震危險性估計結(jié)果不確定性的影響均存在非線性效應(yīng)，即兩個參數(shù)之間存在相互作用.對應(yīng)不同的重現(xiàn)期，在影響地震危險性估計結(jié)果不確定性上，兩個參數(shù)之間的相互作用程度不同.對比不同重現(xiàn)期的震級重現(xiàn)水平，輸入?yún)?shù)u的全效應(yīng)指標(biāo)中，主效應(yīng)占主要部分，交互效應(yīng)的比重隨重現(xiàn)期的增大而逐漸增大；而另一輸入?yún)?shù)ts的全效應(yīng)指標(biāo)中，交互效應(yīng)占主要部分.對于震級上限而言，兩個輸入?yún)?shù)的全效應(yīng)指標(biāo)中，交互效應(yīng)占主要部分，表明這兩個參數(shù)主要是以交互作用的形式對震級上限的不確定性產(chǎn)生影響.

3 討論與結(jié)論

本文選取青藏高原東北緣為案例研究區(qū)，提出了基于全域敏感性分析的地震危險性估計不確定性分析流程和方法.

基于廣義帕累托分布，構(gòu)建了研究區(qū)地震活動性模型，進(jìn)行了地震危險性估計（包括震級重現(xiàn)水平和震級上限的估計）.鑒于構(gòu)建模型需要人為選取ts和u這兩個輸入?yún)?shù)，其取值的不確定性，可能導(dǎo)致地震危險性估計結(jié)果的不確定性，因此，選取具有全域敏感性分析功能的E-FAST方法，定量分析了上述兩個輸入?yún)?shù)的不確定性對地震危險性估計結(jié)果不確定性的影響.結(jié)果表明：地震危險性估計結(jié)果，對兩個輸入?yún)?shù)中的u更為敏感；對應(yīng)不同重現(xiàn)期，地震危險性估計結(jié)果對u的敏感性程度有所不同；兩個輸入?yún)?shù)對地震危險性估計結(jié)果不確定性的影響，存在非線性效應(yīng)；對不同的重現(xiàn)期而言，在影響地震危險性估計結(jié)果不確定性上，兩個輸入?yún)?shù)之間的相互作用程度不同.

目前，確定某一構(gòu)造或構(gòu)造區(qū)的最大震級，常用的方法有概率統(tǒng)計方法和確定性方法.前者有基于G-R關(guān)系統(tǒng)計分析的方法或利用震級-斷層破裂尺度之間經(jīng)驗關(guān)系統(tǒng)計分析的方法；后者有基于研究區(qū)歷史最大地震估計的方法和構(gòu)造類比法等.新近的研究，開始傾向于以形變觀測等資料作為約束，結(jié)合構(gòu)造或區(qū)域的歷史地震信息，估計該區(qū)域的最大震級.

本文選用基于廣義帕累托分布構(gòu)建的地震活動性模型，進(jìn)行地震危險性估計，所用的方法屬于概率統(tǒng)計方法.這一方法，利用歷史地震記錄中大震級區(qū)段的信息，無需先驗地選定震級下限，較為適合歷史地震記錄時間長但低震級地震記錄缺失的地區(qū)，便于構(gòu)建研究區(qū)強(qiáng)震活動性模型，和其它概率統(tǒng)計方法一樣，也有其局限性，即未能考慮地震發(fā)生的物理機(jī)制和物理過程，因為僅利用歷史地震記錄中大震級區(qū)段的信息，有時統(tǒng)計所用的歷史地震記錄樣本偏少.

必須指出，不論通過哪一種方法構(gòu)建地震活動性模型進(jìn)行地震危險性估計，都難以避免不確定性，因此有必要開展相應(yīng)的不確定性和敏感性分析研究，分析地震活動性模型不確定性的來源和特征，即分析輸入?yún)?shù)的不確定性如何影響和制約地震危險性估計的不確定性.這些分析結(jié)果，可以指導(dǎo)地震活動性模型構(gòu)建時輸入?yún)?shù)的選取和調(diào)整，從而使地震活動性模型構(gòu)建得更為完善.

本文提出的定量全域敏感性分析方法，對地震活動性模型的性質(zhì)沒有特殊的限定，便于在地震危險性估計不確定性分析中推廣應(yīng)用.