王鵬 鄭建常

1)山東省地震局，濟(jì)南 250102 2)東營市地震監(jiān)測中心，山東東營 257000

0 引言

2013年以來華北地區(qū)震群活動頻發(fā)，其中山東乳山震群尤為突出。2013年10月1日乳山地區(qū)地震叢集活動持續(xù)至2019年，共記錄到小震1.3萬余次，其中ML≥3.0地震31次，最大震級為2015年5月22日ML5.0，期間伴隨著多次地震活動增強(qiáng)—減弱的過程。

乳山震群所在的膠東半島地區(qū)，位于華北地塊東緣與蘇魯造山帶山東段的交匯位置，在大地構(gòu)造劃分上屬于華北克拉通與揚(yáng)子克拉通的結(jié)合部位(晁洪太等，2001；周瑤琪等，2015)。受蘇魯造山帶與古太平洋板塊俯沖的雙重影響，該區(qū)域巖漿活動十分頻繁，自太古代到新生代均有發(fā)現(xiàn)巖漿活動，其中以燕山晚期巖漿活動最為劇烈，同期還發(fā)生了大規(guī)模的巖漿侵入作用。在火山活動區(qū)域，高熱流區(qū)會導(dǎo)致較小的低粘度累積彈性應(yīng)變，相對于主余型序列而言，往往會產(chǎn)生更多的震群(Ben-Zion et al，2006；Yang et al，2009)。根據(jù)歷史地震記錄，膠東半島地區(qū)歷史上大震較少，卻經(jīng)常發(fā)生震群活動(李冬梅等，2011)，但這些小震序列(震群)持續(xù)時間通常很短，一般持續(xù)幾天，最多不會超過1個月。而乳山震群持續(xù)6年，且頻度和強(qiáng)度均較大，這在大陸東部地區(qū)甚為罕見。對乳山震群震源區(qū)活動特征的深入研究可為區(qū)域地震危險性分析和地震趨勢判斷提供有效的依據(jù)。

火山或地?zé)岬貐^(qū)的震群發(fā)震機(jī)制通常可以用流體運(yùn)移模型來解釋，在地殼中運(yùn)移的流體包括水和揮發(fā)性氣體(CO2)等，其通常被認(rèn)為是由中地殼到上地幔這一深度的巖漿源提供的(Takahata et al，2003；Ohba et al，2011)。在前人對乳山震群發(fā)震機(jī)制的研究中，也猜測流體可能是震群的觸發(fā)因素(Zheng et al，2017；鄭建常等，2019)，但始終缺乏直接的觀測證據(jù)，而在震群發(fā)震構(gòu)造(曲均浩等，2014、2015、2019；鄭建常等，2015a、2015b；李鉑等，2017；張斌等，2017)、震源參數(shù)(苗慶杰等，2016；王鵬等，2016；鄭建常等，2016；劉方斌等，2018)等問題的討論上也存在不同見解。究其原因?yàn)?，乳山地區(qū)斷裂錯綜復(fù)雜，次級斷裂(近SN向)、隱伏斷裂(NE、NW向)較為發(fā)育，斷裂格局總體呈現(xiàn)羽狀或雁列式形態(tài)(劉善寶，2005)，震中區(qū)并無活動斷層，距離最近的乳山斷裂15km。

統(tǒng)計(jì)學(xué)方法可以從復(fù)雜的現(xiàn)象中探究本質(zhì)規(guī)律，自O(shè)gata(1988)在“大森-宇津”公式(Omori，1894；Utsu，1961)的基礎(chǔ)上提出了“傳染型余震序列”(Epidemic-type Aftershock Sequence，ETAS)模型后，該模型已受到廣泛的關(guān)注，其在研究地震活動變化規(guī)律(Ogata，1992、2001；Zhuang，2000；蔣海昆等，2007；蔣長勝等，2013a、2013b、2013c、2015、2017)、探測應(yīng)力變化(雷興林等，2008；Lei et al，2017)、識別與流體作用有關(guān)的誘發(fā)地震(Hainzl et al，2005；Lei et al，2008；龍鋒等，2010；蔣海昆等，2011、2012)等方面得到了長足發(fā)展。

統(tǒng)計(jì)地震學(xué)方法很大程度上依賴于地震目錄，地震目錄的不完備勢必會影響到基于地震目錄的相關(guān)統(tǒng)計(jì)學(xué)結(jié)果的可靠性(李智超等，2014)。地震目錄的完備性主要依賴于研究區(qū)域地震臺網(wǎng)的監(jiān)測能力，而乳山震群周圍臺站分布格局不夠理想，地震序列短時間的叢集發(fā)生通常會造成地震的遺漏。通過遺漏地震檢測來使地震目錄更加完整(Schaff，2008)，可以提高分析結(jié)論的可靠性。用于微震檢測的模板匹配濾波技術(shù)(Gibbons et al，2006)提出后，該技術(shù)發(fā)展迅速并得到了廣泛的應(yīng)用，在檢測地脈動(Shelly et al，2007)、研究余震序列變化(Peng et al，2009；侯金欣等，2017；王鵬等，2016)、確定發(fā)震構(gòu)造(Yang et al，2009；譚毅培等，2016)、分析震源區(qū)應(yīng)力變化(Meng et al，2013)、識別重復(fù)地震(馬騰飛，2015)、監(jiān)視核爆(Zhang et al，2015)等諸多方面發(fā)揮著重要的作用。

本文利用基于GPU并行計(jì)算的模板匹配方法檢測乳山震群的遺漏地震，基于完備地震目錄對該震群ETAS模型的統(tǒng)計(jì)參數(shù)進(jìn)行分析，并初步推測乳山震群的發(fā)震機(jī)制。

1 理論方法

本文所用模板匹配方法的處理流程參照Peng等(2009)，并與侯金欣等(2017)的處理步驟類似，具體方法和處理流程可見Wang等(2019)，在此不再贅述，僅對ETAS模型做簡單介紹。

ETAS模型條件強(qiáng)度函數(shù)λ(t)(Ogata，1988；蔣海昆等，2012)為

(1)

式中，右側(cè)第二項(xiàng)為第i次地震對地震發(fā)生率的貢獻(xiàn)，i遍歷地震序列中的所有地震，Mz為參考震級，其必須大于等于最小完整性震級Mc；p代表地震序列衰減的快慢，p大則衰減快，反之則衰減慢；α表示觸發(fā)次級余震的能力，α越大表示越強(qiáng)的高階余震的激發(fā)能力；μ(t)最初的物理含義為剔除地震叢集后單位時間內(nèi)的理論地震數(shù)目，也表示由于外力作用觸發(fā)微震活動的強(qiáng)弱，與地震間的自激發(fā)過程無關(guān)(Hainzl等，2005)；c為數(shù)值極小的正值，其保證了右側(cè)第二項(xiàng)的分母不為0，與主震后短時間內(nèi)數(shù)據(jù)的完備程度(蔣海昆等，2007)和震源區(qū)的應(yīng)力狀態(tài)有關(guān)(Narteau et al，2009)；K為與主震震級及b、p、α、c、μ參數(shù)有關(guān)的地震活動水平，其中b為G-R關(guān)系中的比例系數(shù)(宋金等，2009)。

對于同一組觀測數(shù)據(jù)，數(shù)學(xué)上可采用多個統(tǒng)計(jì)模型對其進(jìn)行擬合，并可得到多組符合條件的參數(shù)。Ogata(1988)利用最大似然法估計(jì)ETAS模型參數(shù)，并利用最小化模型的信息量準(zhǔn)則(Akaike information criterion，AIC)進(jìn)行最佳模型參數(shù)的判別(Akaike，1974)。AIC準(zhǔn)則既考慮模型對觀測數(shù)據(jù)擬合的優(yōu)劣，也會懲罰為了提高擬合程度而無限制地增加模型參數(shù)的行為。ETAS模型的對數(shù)似然值為

(2)

AIC=-lgL+2k

(3)

其中，tb和te為用于AIC計(jì)算的起止時間；NAIC為計(jì)算時段內(nèi)的地震數(shù)目，為了保證計(jì)算結(jié)果的可靠性，計(jì)算時段內(nèi)的地震目錄必須保證基本完備；k為待擬合參數(shù)個數(shù)，本文中k=5，取AIC最小時對應(yīng)的模型參數(shù)作為最佳模型。

由式(1)可知，ETAS模型所表征的地震活動率包括兩部分，其中右側(cè)第二項(xiàng)是由于“余震激發(fā)余震”所導(dǎo)致的地震叢集，右側(cè)第一項(xiàng)是與余震激發(fā)無關(guān)、因外力作用而產(chǎn)生的地震活動。因而，通過對ETAS模型兩部分地震活動進(jìn)行分離，將外力觸發(fā)地震比例Rb用式(4)表示，而自激發(fā)地震所占比例Rt=1-Rb，可以根據(jù)各自所占的比重來討論外因觸發(fā)與地震活動的關(guān)系(蔣海昆等，2011)

(4)

2 數(shù)據(jù)處理

乳山震群周圍100km內(nèi)共有9個固定地震臺站，均位于震群的北部、西部和東部，而南側(cè)無固定臺站。乳山流動臺陣于2014年5月架設(shè)完成，共計(jì)18套寬頻帶地震儀，采樣率為100Hz(圖1(b))，平均臺間距約為0.5km，距離震群最大距離約10km。乳山臺陣對震群的包絡(luò)性較好，在震群南側(cè)架設(shè)的宮家島臺(GJDT)彌補(bǔ)了震群南側(cè)海島無固定臺站的缺陷。乳山臺陣自2014年5月投入使用，乳山震群持續(xù)時間較長，震群最大地震發(fā)生在2015年5月，地震活動主要集中在2015年6月，后續(xù)地震活動開始逐漸衰減。因此，選擇2014年5月至2015年6月期間固定地震臺站和流動臺陣的連續(xù)波形進(jìn)行遺漏地震檢測。

圖 1 乳山震群和地震臺站空間分布(a)中黑色三角形代表固定臺站，藍(lán)色三角形代表乳山臺陣，灰色圓圈為臺網(wǎng)定位結(jié)果，紅色圓圈為重定位結(jié)果，紅色方框?yàn)檠芯繀^(qū)域，紅色五角星為乳山震群；(b)為圖(a)中黑色方框區(qū)域放大

首先，對連續(xù)波形數(shù)據(jù)去均值、去線性趨勢，并采用零相移4階Butterworth帶通濾波器對2～10Hz進(jìn)行帶通濾波。選取臺網(wǎng)目錄(中國地震臺網(wǎng)中心)中ML≥2.0的地震作為模板地震，根據(jù)臺網(wǎng)目錄中發(fā)震時刻截取乳山臺陣中距離乳山震群最近且分布較好的RSH、GJDT、DFS、NHNC、WHUT和XJZT等6個臺站的三分量記錄作為模板地震波形，再依據(jù)震相報告讀取模板地震的P和S波到時，計(jì)算并選取信噪比大于5的109個地震作為模板地震，最大震級為ML4.5。然后，按照模板匹配方法的數(shù)據(jù)處理流程(Wang et al，2019)，利用模板地震與連續(xù)波形互相關(guān)計(jì)算后做疊加得到互相關(guān)系數(shù)波形，為了減少誤檢測事件，本文選取12倍的絕對離差中位數(shù)(Median Absolute Deviation，MAD)作為判定新檢測地震的閾值，并設(shè)定5s時間范圍內(nèi)僅檢測1次地震事件，選取互相關(guān)系數(shù)最大的遺漏事件，最后確定檢測事件的震級和發(fā)震時刻。

圖 2為檢測到的一次ML2.7遺漏地震的示例。其中，圖2(a)為互相關(guān)系數(shù)和MAD的倍數(shù)，可見在1s的時間內(nèi)互相關(guān)系數(shù)出現(xiàn)多次高值，且均高于12倍MAD的閾值，但根據(jù)假定只取互相關(guān)系數(shù)最大的點(diǎn)，即圖中紅色圓圈所對應(yīng)的時間；圖2(b)表示了遺漏地震所在的波形序列的位置，可見在遺漏地震之前連續(xù)發(fā)生了2次較大地震，分別是2014年9月16日14時42分59秒的ML4.1和14時43分33秒的ML3.9地震，而該遺漏地震幾乎淹沒在前面較大地震的尾波里，幅值是前面較大地震的1/10，不易分辨，這也是地震漏檢測的原因；圖2(c)、2(d)為RSH臺和DFS臺記錄的模板地震在連續(xù)波形上的掃描結(jié)果，可以看到模板地震與檢測出的遺漏事件相關(guān)性較好，相關(guān)系數(shù)為0.61。同時，模板地震的S波的幅值與遺漏事件相當(dāng)，因此根據(jù)S波振幅比計(jì)算的遺漏地震的震級也為ML2.7。依據(jù)檢測事件與所對應(yīng)的模板事件具有相同走時的假設(shè)，確定檢測地震事件的發(fā)震時刻為2014年9月16日14時44分03秒。

圖 2 模板匹配方法檢測遺漏地震示例(以ML2.7地震為例)(a)為經(jīng)過偏移、疊加和取平均值后的互相關(guān)系數(shù)和MAD的倍數(shù)圖，其中相關(guān)系數(shù)大小用豎線表示，MAD的倍數(shù)用橙色點(diǎn)表示；(b)中紅色點(diǎn)線框?yàn)檫z漏地震的位置，即圖(a)所示的時間范圍；(c)、(d)為DFS臺和RSH臺記錄的模板地震和遺漏地震事件的波形對比，其中紅色為2014年11月19日3時24分30秒的ML2.7模板地震，藍(lán)色為該時段內(nèi)的連續(xù)波形，檢測地震事件為2014年9月16日14時44分3秒的ML2.7地震

3 結(jié)果分析

3.1 遺漏地震檢測結(jié)果

2014年5月5日至2015年6月30日期間，共檢測到7500個地震事件，其中109個為模板事件的自檢測，7391個為檢測地震事件，檢測地震事件也包括了臺網(wǎng)已記錄但未用作模板的地震事件。山東地震臺網(wǎng)記錄到地震1937個，檢測地震事件是臺網(wǎng)目錄的近4倍，由模板匹配方法獲得的地震目錄簡稱為檢測目錄。根據(jù)檢測目錄和臺網(wǎng)目錄的頻度對比(圖3)可見，與臺網(wǎng)目錄相比，檢測目錄中的遺漏地震主要集中在ML≤1.0的震級范圍，1.02.5時幾乎沒有遺漏的地震，這與震群所在區(qū)域臺網(wǎng)的監(jiān)測能力有一定關(guān)系。

圖 3 檢測目錄和臺網(wǎng)目錄頻度和累積頻度對比

利用最大曲率法對比分析了臺網(wǎng)目錄和檢測目錄的完備震級，該方法選取震級-頻率曲線(Frequency Magnitude Curve)中斜率最大值所對應(yīng)的震級作為Mc(Wiemer et al，2000)。在實(shí)際應(yīng)用中，Mc震級通常對應(yīng)非累積頻度-震級分布(Non Cumulative Frequency Magnitude Distribution)中具有最大地震頻度的震級(黃亦磊等，2016)。

由圖3 可以看出，在非累積頻度-震級分布中，檢測目錄和臺網(wǎng)目錄的最大值對應(yīng)的完備震級分別為0.1級和0.9級，由于最大曲率方法得到的完備震級通常偏小，利用最大似然法進(jìn)行擬合，發(fā)現(xiàn)檢測目錄和臺網(wǎng)目錄分別在震級0.2和1.0級以上b值的擬合較為穩(wěn)定，最小完整性震級分別取0.2級和1.0級，計(jì)算的b值分別為 0.95±0.01和 0.89±0.03。基于模板匹配方法的遺漏地震檢測將研究時段內(nèi)乳山震群的完備震級從1.0級降至0.2級，彌補(bǔ)了臺網(wǎng)目錄中1.0級以下地震的缺失，使地震目錄更加完整。

3.2 利用ETAS模型分析地震活動特征

乳山震群自2013年10月1日開始，期間經(jīng)歷了多個增強(qiáng)—減弱的階段，根據(jù)日頻度和累積頻度隨時間的演化特征(圖4)，將整個震群發(fā)展過程分為5個活動時段：第一階段為2013年10月1日至2014年1月6日，2013年10月1日12時7分55秒震群初始地震ML3.8地震發(fā)生后，地震活動短時間內(nèi)集中發(fā)生，日頻次達(dá)到高值，并在接下來的2個月內(nèi)逐漸減弱至2014年1月6日；第二階段開始于2014年1月7日，乳山震群發(fā)生了ML4.7地震，隨后地震活動再一次增強(qiáng)，日活動頻度逐漸減弱后又出現(xiàn)了短時間的增強(qiáng)，至2014年4月3日又趨于減弱；第三階段則開始于2014年4月4日的ML4.6地震，與前一階段類似，地震活動在短時間的增強(qiáng)后迅速減弱；第四階段為震群最大地震2015年5月22日ML4.6地震發(fā)生后，日頻度達(dá)到最高值，隨后地震活動經(jīng)過半年多時間的調(diào)整，期間也伴隨著增強(qiáng)—減弱的過程，至2015年底基本趨于穩(wěn)定；第五階段為2016年以后，地震活動較弱，總體較為穩(wěn)定，期間僅有個別地震活動小幅度增強(qiáng)的現(xiàn)象。

圖 4 乳山震群日頻度和累積頻度隨時間的演化左側(cè)縱坐標(biāo)表示日頻度，用黑色點(diǎn)線表示；右側(cè)縱坐標(biāo)表示累積頻度，用紅色線表示；藍(lán)色點(diǎn)線將整個發(fā)展過程分為5個階段，黑色點(diǎn)線又將每個階段分為若干個子階段

圖 5 乳山震群震級、頻度-時間圖(a)為檢測目錄(黑點(diǎn))和臺網(wǎng)目錄(灰點(diǎn))的震級隨時間變化圖，其中紅色點(diǎn)線框?yàn)楸疚睦媚０迤ヅ浞椒z測的時間段； (b)為圖(a)紅色框的放大，并加入累積頻度的統(tǒng)計(jì)，其中藍(lán)點(diǎn)為檢測目錄，黑點(diǎn)為臺網(wǎng)目錄，紅色線為累積頻度曲線

根據(jù)震群發(fā)展過程，每個階段又分為若干個小的階段(圖4)。將檢測得到的完備震級加入到臺網(wǎng)目錄中，利用ETAS模型對這些不同的階段分別進(jìn)行反演，得到各模型參數(shù)的變化特征，從而得出外力作用和自激發(fā)作用對震群發(fā)展的影響。

不同的截止震級對ETAS模型參數(shù)的反演影響較大(蔣海昆等，2012；蔣長勝等，2013c)，要保證每個階段的ETAS模型參數(shù)反演中截止震級必須大于完備震級。因此首先利用震級-序號法考察完備震級在震群發(fā)展過程中的變化(圖6)，震級-序號法是按地震發(fā)生的次序排列，通過分析不同震級的地震概率密度分布來定性分析完備震級Mc的方法，其中地震概率密度最大的位置所對應(yīng)的震級即為完備震級Mc。之所以使用地震發(fā)生次序，而不是地震時間，主要是為了避免余震或震群中“叢集事件”的影響(蔣長勝等，2013a)。這里以500個地震作為窗長，以0.1級為間隔，計(jì)算窗內(nèi)各個震級地震的概率密度，可以看到在震群發(fā)生初期完備震級在0.9級左右，在檢測目錄加入后使得完備震級明顯降低，在檢測的2014年5月至2015年6月期間，完備震級基本維持在0.1級左右，在2015年6月之后完備震級又上升到0.9級左右。為了分析乳山震群整個發(fā)展過程中ETAS模型參數(shù)的趨勢變化，并能更好地與其他階段的參數(shù)變化進(jìn)行對比，因此統(tǒng)一將截止震級取為0.9級，按照圖4 中的階段劃分，反演了各個階段的ETAS模型參數(shù)，見表1，其中第一階段(Ⅰ)至第四階段(Ⅳ)以及每個階段的子階段(Ⅲ1、Ⅲ2，Ⅳ1、Ⅳ2)的擬合圖見圖7。

圖 6 乳山震群震級-序號圖色標(biāo)表示計(jì)算窗內(nèi)地震在某個震級段的概率密度分布

表 1乳山震群不同階段的ETAS模型參數(shù)

圖 7 乳山震群各階段ETAS模型擬合結(jié)果

對于第一階段，背景地震或外力作用的比例為35.5%，說明有35.5%的地震是由外力作用觸發(fā)的，剩余64.5%的地震由余震的自激發(fā)產(chǎn)生。α值為0.947，低于典型的構(gòu)造地震(2.0)，震群型地震的α值通常小于1(Ogata，2001)。日本地區(qū)震群的α值一般介于0.35～0.85之間，而一般的主余序列則在1.2～3.1之間(Ogata，1992)。α值表征余震的激發(fā)能力，這一階段的α值接近1，表明有較強(qiáng)的余震的激發(fā)能力。p值為余震衰減率，這里p值(0.99)與典型構(gòu)造地震(1.0左右)相當(dāng)，表明序列屬于正常衰減。同時，在這一階段模型理論值與觀測值擬合較好，表明參數(shù)估計(jì)較為準(zhǔn)確，外力作用在震群發(fā)生初期占一定比重，但主要還是受余震的自激發(fā)作用影響。

第二階段的地震活動伴隨著較大地震(2014年1月7日ML4.7地震)短時間內(nèi)集中發(fā)生，隨后逐漸衰減，期間發(fā)生過1次日頻度突然增強(qiáng)的現(xiàn)象。這一階段外力作用所占比例有所下降，從第一階段的35.5%降至24.3%，表明外力作用較弱，p值無變化，仍為0.99，表明序列仍處于正常衰減，但α值較高(1.258)，且高于第一階段，表明余震的激發(fā)能力增強(qiáng)。

第三階段與第二階段類似，伴隨著較大地震(2014年4月4日ML4.6地震)的發(fā)生，地震活動出現(xiàn)叢集，隨后趨于衰減，由于地震活動速率變化較大，將該階段分為2個子階段(Ⅲ1、Ⅲ2，圖4)，其中Ⅲ1階段外力作用所占比例有所上升(39.9%)，表明外力作用有增強(qiáng)的趨勢，這種短期叢集的地震受外力作用觸發(fā)的較多。α值明顯下降，由1.258降至0.362，表明余震的激發(fā)能力明顯減弱，與外力作用增強(qiáng)一致。同時p值也明顯增大，變?yōu)?.5，表明序列的衰減作用增強(qiáng)。在這一階段，模型理論值與觀測值擬合較差，實(shí)際地震發(fā)生率(圖7(c)藍(lán)色實(shí)線)明顯高于擬合的理論曲線(圖7(c)灰色點(diǎn)線)，表明大森定律型自觸發(fā)只有很小的促進(jìn)作用，外力作用引起的地震活動的比例可能被低估(Lei et al，2008)。在Ⅲ2階段，外力作用所占比例繼續(xù)增大(45.9)，說明將近一半的地震是由于外力作用觸發(fā)的，同時α值明顯升高(1.444)，表明這一活動時期地震活動的自激發(fā)作用也相應(yīng)增強(qiáng)，p值也再次回到1.01，序列呈現(xiàn)天然地震的正常衰減速率。

第四階段地震活動速率明顯增加，也伴隨著震群最大地震(2015年5月22日ML5.0地震)的發(fā)生，日頻次達(dá)到最高值，隨后快速下降。將這一活動時段也分為2個子階段(Ⅳ1、Ⅳ2，圖4)，其中Ⅳ1階段外力作用較上一階段有所降低(35.2%)，約有 1/3 的地震是由外力因素觸發(fā)的，而α值依然保持較高值(1.002)，表明這一階段更多的是受最大地震(ML5.0)的影響，激發(fā)了短期內(nèi)的余震發(fā)生，p值為1.06，表明仍在正常衰減。Ⅳ2階段在衰減(p值)變化不大的情況下，α值略微下降(0.969)，說明受最大地震(ML5.0)的影響減弱，同時受外力作用明顯增強(qiáng)(43.6%)，與Ⅲ1階段類似，該階段的實(shí)際地震發(fā)生率(圖7(f)藍(lán)色實(shí)線)也明顯高于擬合的理論曲線(圖7(f)灰色點(diǎn)線)，表明外力作用也可能被低估(Lei et al，2008)。外力作用在這一階段可能起主要作用，而地震的自激發(fā)作用次之。

第五階段隨著地震活動的減弱，外力作用比例減弱，α值和p值也有所下降，說明這一階段外力作用和自激發(fā)作用均在減弱，同時衰減也相應(yīng)減弱，與這一時期持續(xù)較長時間的弱地震活動水平有關(guān)。

我們還考察了乳山震群的ETAS模型參數(shù)Rb、α和p隨時間的演化規(guī)律(圖8)，可以看到除個別點(diǎn)外，Rb、α的變化趨勢較為一致，總體呈現(xiàn)先上升后下降的趨勢，而p值除在個別點(diǎn)外，其他時段變化較小，基本維持在1.0左右，表明乳山震群一直保持著天然地震正常的平均衰減速率。Rb值在整個過程中基本維持在20%～50%之間，表明外力作用在震群的發(fā)展過程起到一定的作用，其在某些階段甚至接近50%，但總體來看仍是震群前期地震的自激發(fā)起主導(dǎo)作用。

圖 8 ETAS模型參數(shù)隨時間的變化

4 討論

4.1 不同截止震級對ETAS模型反演結(jié)果的影響

通過對2014年5月至2015年6月期間的遺漏地震進(jìn)行檢測，使得完備震級大幅減小，在完備震級隨時間的演化中也能明顯看到這種特征。2014年5月至2015年6月這一時間段對應(yīng)Ⅲ2、Ⅳ1階段(圖4)，因此，測試這2個階段在截止震級分別取0.1級和0.9級時ETAS模型的擬合結(jié)果，比較不同參數(shù)的結(jié)果(圖9)。從圖9 可以看到，外力觸發(fā)地震的比例Rb差別較大，截止震級取0.1級時2個階段的Rb值明顯低于截止震級取0.9級的情況。由于任何一個地震均能以一定概率觸發(fā)后續(xù)地震，低于Mc的大量微震事件將可能顯著改變背景地震概率(Rb值)(蔣長勝等，2013b)。在III2階段，隨著截止震級的減小，α值和μ值明顯減小，表明序列觸發(fā)微震的能力和外力作用觸發(fā)微震活動的能力均有所減弱。在Ⅳ1階段，隨著截止震級的減小，α值和μ值則有所增加，表明自激發(fā)作用和外力作用均有所增強(qiáng)。這種現(xiàn)象可能與2個階段的地震活動特征不同有關(guān)，III2階段地震活動趨勢較為平穩(wěn)，而Ⅳ1階段由于震群最大地震的發(fā)生，地震活動頻度迅速增加，日頻度達(dá)到最大值(圖5)，震級越大，激發(fā)高階余震的能力越強(qiáng)，對后續(xù)地震的影響也就越大。盡管α值和μ值均有所增加，受震群最大地震的影響，外力觸發(fā)地震的比例卻很低(3.7%)，地震活動中仍是自激發(fā)作用的貢獻(xiàn)更大。另一方面由于震群中最大地震發(fā)生后，短期內(nèi)ETAS模型參數(shù)隨時間變化較為明顯且波動較大，ETAS模型參數(shù)與地震序列持續(xù)時間也有一定關(guān)系(蔣長勝等，2013c)，2個階段的持續(xù)時間相差較大，可能也是一方面因素。此外，震源區(qū)的破裂特征、斷層的愈合速度、應(yīng)力加載、周邊構(gòu)造場的調(diào)整等諸多因素均會使地震序列的衰減特征和激發(fā)余震能力出現(xiàn)明顯差異。

由于ETAS模型假定每一次地震均可觸發(fā)其后發(fā)生的任何震級的地震，因此，截止震級在確保地震目錄完整性的同時，也可能切斷了截止震級以下地震與以上地震的觸發(fā)關(guān)系，從而造成“鏈接缺失”現(xiàn)象(Wang et al，2010)。通過遺漏地震檢測，使得完備震級減小，而截止震級的減小盡可能地減少了這種現(xiàn)象的出現(xiàn)，因此地震目錄的完整性對ETAS模型參數(shù)反演至關(guān)重要。

經(jīng)過遺漏地震檢測，在保證地震目錄完整性的基礎(chǔ)上，降低了截止震級，計(jì)算得到的ETAS模型參數(shù)也更加可靠。從理論曲線和實(shí)際地震活動的擬合來看，截止震級取0.1級時擬合較好，同時AIC值也更低(圖9)。但由于乳山臺陣于2014年5月才架設(shè)運(yùn)行，而乳山震群2013年10月開始集中活動，因此基于臺陣波形數(shù)據(jù)的遺漏地震檢測的時間段僅是震群發(fā)展過程中的一部分(圖5)。為了考察乳山震群整個發(fā)展過程中ETAS模型參數(shù)的趨勢變化(圖8)，同時能更好地與其他階段進(jìn)行對比(圖7)，對于整個震群發(fā)展過程取截止震級0.9級以上的完備地震目錄進(jìn)行反演計(jì)算。將截止震級設(shè)為高于完備震級(0.9級)的某個震級，利用ETAS模型進(jìn)行反演，盡管各參數(shù)的數(shù)值不完全相同，但反映出各參數(shù)隨時間的變化規(guī)律是一致的。

圖 9 不同截止震級下ETAS模型擬合結(jié)果對比(a)、(b)截止震級為0.1級時2個階段的擬合結(jié)果；(c)、(d)截止震級為0.9級時2個階段的擬合結(jié)果

4.2 發(fā)震機(jī)制討論

本文利用外力作用引起的地震與自激發(fā)地震的比例系數(shù)Rb來考察外力因素對震群活動的影響。除個別階段外，Rb和α的總體變化趨勢較為一致，呈現(xiàn)出先升高后降低的過程，而p值除個別階段外，總體變化較小，維持在天然地震正常的衰減速率。Rb值在整個過程中基本維持在20%～50%之間，表明外力作用在震群的發(fā)生、發(fā)展過程中均占不小的比例，尤其在震群初期或較大地震發(fā)生之前，地震可能是受外力因素觸發(fā)或間接引起的應(yīng)力調(diào)整引發(fā)的，但從整個過程的比例來看，余震的自激發(fā)作用所占比重更大，發(fā)揮著更大的作用。震群發(fā)展后期(V階段)外力作用和自激發(fā)作用均在減弱，地震活動呈現(xiàn)較長時間的弱活動水平。

前人在對乳山震群空間分布、震中遷移速率、ETAS模型、高頻衰減系數(shù)等方面的研究中推測乳山震群可能是由流體入侵觸發(fā)的(鄭建常等，2015b、2016；Zheng et al，2017)，但目前仍沒有確切的觀測證據(jù)。通過分析ETAS模型參數(shù)特征和其隨時間演化的規(guī)律，也可以看到外力因素在乳山震群的發(fā)生、發(fā)展過程中起著一定的作用，而流體入侵或運(yùn)移也可能是外力因素之一，今后還需結(jié)合其他依據(jù)做詳細(xì)的論證。

5 結(jié)論

本文利用基于GPU并行運(yùn)算的模板匹配濾波技術(shù)，對于乳山震群2014年5月5日至2015年6月30日期間的連續(xù)地震波形進(jìn)行檢測，識別出7500個地震事件，其中109個為模板事件的自檢測，7391個為檢測地震事件。檢測事件是臺網(wǎng)目錄的4倍，利用較為準(zhǔn)確的b值擬合方法得到的震群完備震級從1.0級降至0.2級。

乳山震群2013年10月開始集中活動，而乳山臺陣自2014年5月5日才開始運(yùn)行，因此，將檢測目錄補(bǔ)充到整個乳山震群的地震目錄中，并劃分整個震群發(fā)展過程為5個活動時段，借助ETAS模型分析了乳山震群的地震活動性以及在各個階段的變化特征。結(jié)果表明，乳山震群外力觸發(fā)的地震所占比例較大，結(jié)合前人研究，初步推測流體作用可能是外力因素之一，且不同時期外力對微震活動的觸發(fā)強(qiáng)度不同。