張正帥,陳時軍,周 晨,趙 瑞

(1.山東省地震局，濟南 250014；2.黑龍江省地震局，哈爾濱 150000)

0 引言

地震活動時空分布往往表現(xiàn)出明顯的非均勻性[1]。盡管大量的地震事件表現(xiàn)出較好的統(tǒng)計規(guī)律，用以刻畫地震活動的總體規(guī)律，如描述余震衰減時間分布的Omori模型[2]以及隨后發(fā)展的ETAS模型[3]，描述震級分布的Gutenberg-Richter震級頻度關系[4]等。然而，地震事件時空分布往往表現(xiàn)出強烈的非獨立、非隨機分布的特性，使得這些統(tǒng)計規(guī)律的普適性受到挑戰(zhàn)。進一步研究發(fā)現(xiàn)，不同地震事件之間存在著大時空范圍的相互作用或相互關聯(lián)[5]，一次強震的震源斷裂錯動所造成的應力場變化可能擾動其后的區(qū)域地震活動[6]。例如，美國蘭德斯地震后誘發(fā)的距震中遠達1 250 km的14個區(qū)域出現(xiàn)地震活動增強現(xiàn)象[1]等，表明地震活動時空分布具有明顯的復雜性，而這種復雜性與地殼內(nèi)存在的復雜構造體系以及復雜的動力學系統(tǒng)有關。近年來發(fā)展起來的復雜網(wǎng)絡模型，大量用于研究自然科學、社會科學、工程技術等領域存在相互關聯(lián)的復雜動力系統(tǒng)。復雜網(wǎng)絡是對真實系統(tǒng)的抽象，它可以從整體的角度研究復雜系統(tǒng)的結構和功能，探究復雜系統(tǒng)內(nèi)在的相互聯(lián)系。復雜網(wǎng)絡研究的開創(chuàng)性工作來源于Watts等[7]發(fā)現(xiàn)的真實網(wǎng)絡具有小世界性質及Barabási等[8]發(fā)現(xiàn)的真實網(wǎng)絡的無標度特性。復雜網(wǎng)絡由此作為一種新理論受到科學界的大量關注，已經(jīng)得到應用的領域涉及生物、電力、交通和社交網(wǎng)絡等。近些年由Abe等[9-14]將復雜網(wǎng)絡的概念引入到地震學研究之中，并做了大量工作，之后陸續(xù)有成果發(fā)表：2011年，謝周敏[15]利用加權復雜網(wǎng)絡模型研究了地震網(wǎng)絡的拓撲結構和動力學行為，2014年，何璇等[16-17]提出基于時空影響域的地震網(wǎng)絡構造方法，2015年，趙海等[18]研究了美國南加州地區(qū)地震網(wǎng)絡的規(guī)模和熵演化；2015年，李光光等[19]基于地震網(wǎng)絡的k-核解析了地震活動分布特征。2017年，Denisse Pastén等[20]對智利南北兩個地區(qū)的地震活動不同震級閾值，不同網(wǎng)格劃分等情況的復雜網(wǎng)絡分析，發(fā)現(xiàn)介數(shù)中心性的非普遍性；2018年，Soghra Rezaei等[21]對動態(tài)地震網(wǎng)絡的偏好連接性進行了研究，尋找最容易被影響的地震節(jié)點。本文采用Abe等[12]提出的地震網(wǎng)絡構造方法，將地震目錄轉換成基于多圖的復雜網(wǎng)絡，由此構造的網(wǎng)絡表示地震活動的動態(tài)信息，通過對該復雜網(wǎng)絡量化分析，探究地震活動的復雜性。

1 地震網(wǎng)絡構造方法及數(shù)據(jù)

1.1 地震網(wǎng)絡構造方法

地震網(wǎng)絡構造方法如下：首先對研究區(qū)的經(jīng)度、緯度和深度進行三維網(wǎng)格化，得到若干大小相等的地理立方單元，如圖1所示。在所研究的時空范圍內(nèi)，若該地理單元內(nèi)至少發(fā)生過1次地震(震中位于單元的地理范圍內(nèi))，則將該地理單元抽象為地震網(wǎng)絡中一個節(jié)點；若2次相繼發(fā)生的地震震中處于不同的地理單元內(nèi)，那么這兩個節(jié)點間產(chǎn)生1條邊；若2次相繼發(fā)生的地震震中處于同一個地理單元內(nèi)，那么該節(jié)點產(chǎn)生1個自環(huán)。網(wǎng)格大小的取值可依據(jù)地震定位精度以及網(wǎng)絡的穩(wěn)定性確定，謝周敏[15]認為，網(wǎng)格邊長在5km至10km范圍內(nèi)變化時，其所構造的地震網(wǎng)絡具備較好的穩(wěn)定性[22]。根據(jù)本文所選資料的地震定位精度，每個網(wǎng)格的緯度、經(jīng)度及深度大小分別選擇為5km×5km×5km和10km×10km×10km進行研究，并對結果進行對比分析。

地理單元尺寸為0.2°×0.2°圖1 研究區(qū)網(wǎng)格劃分示意圖Fig.1 A schematic description of meshing in the research region

度大于5圖2 2013年地震網(wǎng)絡拓撲結構Fig.2 Earthquake network topology in 2003

1.2 數(shù)據(jù)

研究地區(qū)位于青藏高原東緣，發(fā)育有多條活動斷裂帶。新生代以來，受青藏高原物質東南方向擠出作用影響而構造活動劇烈，GPS測量結果顯示該地區(qū)現(xiàn)今地殼活動強烈且構造變形機制復雜。該地區(qū)是中國大陸最顯著的地震活動區(qū)域之一，地震樣本量大，為本研究提供了良好的數(shù)據(jù)基礎。為分析該地區(qū)現(xiàn)代地震活動特征，本文選用研究區(qū)(29°N-34°N，100°E-106°E)2006年1月1日至2017年12月31日之間的地震目錄構建地震網(wǎng)絡，期間發(fā)生M0.0以上地震共計132 377次。由于地震網(wǎng)絡是對真實地震現(xiàn)象的抽象，小震也會對整個復雜網(wǎng)絡系統(tǒng)產(chǎn)生一定的影響，因此本研究不可忽略小地震樣本?？紤]到地震目錄的完備性和一致性，根據(jù)劉麗芳等[23]，龍鋒等[24]的研究，研究區(qū)內(nèi)的最小完整性震級Mc=1.5左右，因此本文構造地震網(wǎng)絡的震級下限取M1.5，震源深度取0km至80km，并剔除非天然地震事件，經(jīng)過篩選，研究區(qū)范圍10年內(nèi)共計發(fā)生地震49 826次。本文構造地震網(wǎng)絡采用的地震目錄來自于國家地震科學數(shù)據(jù)共享中心(http://data.earthquake.cn)。圖1為研究區(qū)網(wǎng)格劃分示意圖。作為示例，圖2給出了2013年地震數(shù)據(jù)構造的地震網(wǎng)絡，其中共有節(jié)點2 232個，邊數(shù)6 734條，為了更直觀展示地震網(wǎng)絡的拓撲結構，圖2中只給出了度大于5的節(jié)點予以顯示。

2 地震網(wǎng)絡結構特性

該部分對所構建的地震網(wǎng)絡的特性予以研究。首先，分析地震網(wǎng)絡的無標度特性；然后從地震網(wǎng)絡的平均路徑長度和聚類系數(shù)的考察中，分析了地震網(wǎng)絡的小世界特性。地震網(wǎng)絡本身是一種有向圖，方向性對研究地震事件周期性(揭示在發(fā)生了多少次后續(xù)地震后，地震事件返回到初始節(jié)點)是重要的[25]。然而，考慮到本文在分析地震網(wǎng)絡時，只關心節(jié)點間是否連通的靜態(tài)特性；另外，當研究地震網(wǎng)絡的小世界性質時，必須忽略方向性，并將路徑長度定義為連接一對頂點的邊數(shù)中可能最小的值[25]。本文將構造的地震網(wǎng)絡轉換為無向的簡單圖。首先移除網(wǎng)絡中的自環(huán)，其次將兩個節(jié)點之間存在的多邊轉換為單邊，并且忽略邊的指向。

2.1 無標度特性

無標度網(wǎng)絡中大部分節(jié)點只有少數(shù)的連接，而某些節(jié)點卻擁有與其他節(jié)點的大量連接，度分布滿足冪律性，即P(k)～k-γ。P(k)是度值為k的節(jié)點出現(xiàn)的概率，其中冪指數(shù)γ為正指數(shù)，范圍在1到2之間。本文通過最大似然估計方法計算得到冪指數(shù)。圖3和圖4給出了研究區(qū)2006年1月1日至2017年1月31日共10年尺度的地震網(wǎng)絡度分布。其中圖3、圖4分別是網(wǎng)格大小為5km×5km×5km和10km×10km×10km的結果。

圖3 節(jié)點度分布雙對數(shù)曲線(網(wǎng)格大小為5km×5km×5km)Fig.3 The log-log plot of degree distribution with the cell size 5km×5km×5km

圖4 節(jié)點度分布雙對數(shù)曲線(網(wǎng)格大小為10km×10km×10km)Fig.4 The log-log plot of degree distribution with the cell size 10km×10km×10km

從圖3、圖4中可見，兩種不同尺寸的網(wǎng)格所構造的地震網(wǎng)絡的度分布P(k)滿足冪律分布，顯示該地震網(wǎng)絡具有無標度特點。兩個地震網(wǎng)絡度分布都呈現(xiàn)冪律性，表明地震活動中震源之間的關聯(lián)程度具有較強的異質性，即震源與震源之間關聯(lián)程度很大與關聯(lián)程度很小的各種情況都存在[15]。冪指數(shù)γ是有所差異的，γ越小，說明網(wǎng)絡在度分布上的非均勻性越強，即某些中心節(jié)點的度越大；反之，度分布均勻性越強。相對來說以10km×10km×10km構造的地震網(wǎng)絡比以5km×5km×5km構造的地震網(wǎng)絡的度分布相對更均勻一點。兩種網(wǎng)格尺寸下地震網(wǎng)絡都是在大量地震事件統(tǒng)計意義下滿足冪律性，地震網(wǎng)絡都表現(xiàn)出較好的穩(wěn)定性。另外，這種冪律性特征也說明不同地震事件之間互相作用蘊含深層的物理意義。

2.2 小世界特性

對比而言，根據(jù)所構建的地震網(wǎng)絡節(jié)點數(shù)和邊數(shù)隨機產(chǎn)生的Erd?s-Rényi隨機網(wǎng)絡模型(以下簡稱ER模型)，其網(wǎng)絡中兩個節(jié)點之間不論是否具有共同的鄰居節(jié)點，其連接概率均為p，因此ER模型的聚類系數(shù)為：CER=〈k〉/(N-1)≈〈k〉/N=p≤1，其中〈k〉表示網(wǎng)絡的平均度。因此，相比ER模型，具有小世界特性的網(wǎng)絡具有更高的聚類系數(shù)。

具有小世界性質的網(wǎng)絡平均路徑長度L與ER模型一樣具有較小的平均路徑長度，由于ER模型的平均路徑長度LER∝lnN/ln〈k〉，lnN隨著N增長得很慢，使得即使規(guī)模很大的ER網(wǎng)絡也具有很小的平均路徑長度。表1給出了本文構建的地震網(wǎng)絡與ER網(wǎng)絡的參數(shù)對比。

表1 研究區(qū)地震網(wǎng)絡與ER網(wǎng)絡結構參數(shù)對比(2006年至2017年)Tab.1 The structure parameter comparison for earthquake network and ER network(2006—2017)

從表1可以看出，本文構建的地震網(wǎng)絡相比ER網(wǎng)絡具有明顯的高聚類系數(shù)，而且兩者之間都具有較短的平均路徑，說明該地震網(wǎng)絡具備小世界特性。

3 大震前后地震網(wǎng)絡演化分析

就地震活動而言，地震的不斷發(fā)生意味著地震網(wǎng)絡結構不斷發(fā)生改變，通過分析地震網(wǎng)絡結構特性的演變，對探究大地震發(fā)生前后中小地震活動特征從而指導大地震預測具有重要減災意義。結合Douglas等[27]提出的基于時間窗口的地震網(wǎng)絡分析方法，本文設置時間窗口T，并以一定步長dT對地震目錄進行滑動掃描，將每一個窗口所限定時間范圍內(nèi)的地震事件構造成一個新的地震網(wǎng)絡，其中每個時間窗的窗尾作為窗口標志時間點，從而計算隨窗口演化的地震網(wǎng)絡特征參數(shù)值。

本研究主要考察研究地區(qū)M7.0以上地震前后的地震網(wǎng)絡演化特征，包括2008年5月12日汶川M8.0地震、2013年4月20日蘆山M7.0地震、2017年8月8日九寨溝M7.0地震。本文網(wǎng)格單元選擇5km×5km×5km，窗口參數(shù)選擇為T=10d，dT=1d，d表示天數(shù)。通過這樣一組參數(shù)對地震前后6個月內(nèi)地震活動數(shù)據(jù)進行滑動掃描計算，進行演化分析。

3.1 地震網(wǎng)絡規(guī)模演化

地震網(wǎng)絡的節(jié)點數(shù)和邊數(shù)反映了網(wǎng)絡的規(guī)模(大小)，是網(wǎng)絡最直觀的特征參量，所以首先考察地震前后網(wǎng)絡規(guī)模的變化情況。圖5、6、7分別給出了汶川M8.0地震，蘆山M7.0地震和九寨溝M7.0地震的計算結果。其中，圖a表示地震活動的M-T圖(震級-時間圖)。

從圖5、6、7中的M-T圖可以看出，無論是汶川地震，還是蘆山地震與九寨溝，震前地震活動性比較均勻，難以發(fā)現(xiàn)明顯的變化情況。根據(jù)網(wǎng)絡規(guī)模演化分析，汶川M8.0地震前3個月，網(wǎng)絡節(jié)點數(shù)和邊數(shù)出現(xiàn)一定的增高—降低現(xiàn)象，而蘆山M7.0地震、九寨溝M7.0地震前網(wǎng)絡邊數(shù)沒有明顯變化，但節(jié)點數(shù)存在明顯的波動變化—降低的過程，表明網(wǎng)絡規(guī)模的穩(wěn)定性變差。然而，3個大地震之前的規(guī)模演化情況并不一致，所以僅僅通過網(wǎng)絡規(guī)模的變化情況，很難提取到明顯的異常信息，也就難以判斷是否會發(fā)生大地震。但是，地震發(fā)生后，網(wǎng)絡規(guī)模出現(xiàn)大尺度變化，且隨著時間流逝，節(jié)點數(shù)和邊數(shù)逐漸減少。這個現(xiàn)象與余震衰減規(guī)律是一致的，主要的物理機制在于大地震激發(fā)了彈性波能量的釋放，進而觸發(fā)大量余震，導致地震網(wǎng)絡規(guī)模的變化，這也從復雜網(wǎng)絡的角度對余震變化規(guī)律有了更深層次的理解。

3.2 k-核解析

Abe等[28]證明了地震網(wǎng)絡具有層次結構。度小的節(jié)點群以分層的方式組織起來。然而，我們不清楚哪些節(jié)點屬于哪個層，不同層之間的差異性還不夠清晰。研究表明，隨著網(wǎng)絡規(guī)模的增長，其節(jié)點k-核值的變化會趨于穩(wěn)定[29]，k-核值是一種更加穩(wěn)定、簡單的參數(shù)[30]。所以，為了獲得網(wǎng)絡層次的細節(jié)信息，本文利用k-core分解對網(wǎng)絡進行解析，通過聚焦于地震網(wǎng)絡最高核，研究最高核隨時間的演化特征。

圖5 汶川M8.0地震前后地震網(wǎng)絡規(guī)模演化Fig.5 Scale evolution of the earthquake network before and after the Wenchuan M8.0 earthquake

圖6 蘆山M7.0地震前后地震網(wǎng)絡規(guī)模演化Fig.6 Scale evolution of the earthquake network before and after the Lushan M7.0 earthquake

圖7 九寨溝M7.0地震前后地震網(wǎng)絡規(guī)模演化Fig.7 Scale evolution of the earthquake network before and after the Jiuzhaigou M7.0 earthquake

k-核是1983年由Seidman等[31]提出的一種簡化網(wǎng)絡拓撲結構的方法。k-核在復雜網(wǎng)絡中的應用是由Gaertler等[32]以及Gkantsidis等[32]提出的，其定義為：一個網(wǎng)絡圖中的k-核是指反復去掉度數(shù)小于和等于k的節(jié)點后剩余的子圖。節(jié)點的核數(shù)表示包含該節(jié)點的最深的核，即節(jié)點存在于k-核中，但是在(k+1)-核中被移除，則該節(jié)點的核數(shù)為k。因此可以通過k-核解析由外層至內(nèi)層，一層一層地解析網(wǎng)絡，從而揭示網(wǎng)絡的層次結構性質。節(jié)點核數(shù)的最大值為網(wǎng)絡的最大核數(shù)，也稱最高核數(shù)[34]。在地震網(wǎng)絡中，k-核的重要特點是節(jié)點的連通性，一個節(jié)點的核數(shù)越大，則通過該節(jié)點到達其他節(jié)點的路徑就越多，也就是說這個節(jié)點的影響域更大。根據(jù)3.1部分選取的網(wǎng)格大小以及窗口參數(shù)，對地震前后數(shù)月內(nèi)的地震網(wǎng)絡進行最高核數(shù)計算，結果如圖8所示。

圖8給出了汶川M8.0地震、蘆山M7.0地震、九寨溝M7.0地震期間，地震網(wǎng)絡最高核數(shù)的演化量。對于蘆山M7.0地震和九寨溝M7.0地震而言，震前3個月左右出現(xiàn)高值異常，且兩者具有一定的相似性。震后最高核的值會逐步下降到5以下，恢復較為穩(wěn)定的狀態(tài)。但汶川M8.0地震前地震網(wǎng)絡k-核并未出現(xiàn)明顯異常，表明地震網(wǎng)絡層次性結構可能未發(fā)生明顯改變，這可能與汶川地震孕震范圍大，亦或與汶川地震破裂方式、孕育環(huán)境等因素有關。由于網(wǎng)絡層次性表現(xiàn)為度很小的節(jié)點具有高聚類系數(shù)，且屬于高度連接的小模塊，相反，度很高的節(jié)點具有低聚類系數(shù)，其作用是將不同模塊連接起來[26]。大量實證研究表明，許多真實網(wǎng)絡中節(jié)點的聚類系數(shù)與度存在近似的倒數(shù)關系[35-36]，說明復雜網(wǎng)絡中聚類系數(shù)在一定程度上可以刻畫網(wǎng)絡的層次結構。為進一步描述這3次地震前后地震網(wǎng)絡結構的層次特征，本文計算了上述3個地震前后幾個月的地震網(wǎng)絡聚類系數(shù)，結果如圖9所示。

左箭頭表示發(fā)震時間，黑色虛線表示最高核數(shù)均值線圖8 地震網(wǎng)絡最大核數(shù)的演化Fig.8 The evolution of the highest layer of the earthquake network

左箭頭表示發(fā)震時間，黑色虛線表示平均聚類系數(shù)均值線圖9 地震網(wǎng)絡平均聚類系數(shù)演化Fig.9 The evolution of the average clustering coefficient of the earthquake network

從圖9中可見，汶川M8.0地震前地震網(wǎng)絡的平均聚類系數(shù)并沒有發(fā)生大的波動，說明網(wǎng)絡的層次結構未明顯改變。所以，汶川地震之前，最高核數(shù)并未出現(xiàn)明顯的異常特征。

4 結論與討論

本文采用基于時間序列構造地震網(wǎng)絡的方法，將研究地區(qū)的地震目錄數(shù)據(jù)映射為一種復雜網(wǎng)絡的形式，通過運用復雜網(wǎng)絡的方法研究地震數(shù)據(jù)的統(tǒng)計規(guī)律，根據(jù)本文的研究可得到以下結論：

1)通過計算地震網(wǎng)絡的聚類系數(shù)和最短路徑長度，通過與同規(guī)模的ER隨機網(wǎng)絡相比，發(fā)現(xiàn)地震網(wǎng)絡具有明顯的高聚類系數(shù)及較短的平均路徑的特點，說明地震網(wǎng)絡的小世界特性。

2)度分布滿足冪律分布，說明地震網(wǎng)絡具有無標度特性。

3)汶川M8.0地震、蘆山M7.0地震和九寨溝M7.0地震前3個月左右，網(wǎng)絡規(guī)模穩(wěn)定性有一定的擾動，其中汶川M8.0級地震前無論網(wǎng)絡節(jié)點數(shù)和邊數(shù)均出現(xiàn)小幅高值-降低的過程，而蘆山M7.0級地震和九寨溝M7.0級地震前3個月內(nèi)網(wǎng)絡邊數(shù)穩(wěn)定但節(jié)點數(shù)存在波動升高-降低的過程。但是，僅僅通過網(wǎng)絡規(guī)模的變化情況，很難提取到明顯的大震前的異常信息。

4)大地震之后，網(wǎng)絡規(guī)模有明顯增大，節(jié)點數(shù)和邊數(shù)在震后的數(shù)量顯著增多，其主要原因在于主震觸發(fā)了震區(qū)不同范圍內(nèi)的大量余震，對網(wǎng)絡規(guī)模造成了影響，這也從復雜網(wǎng)絡的角度說明了地震事件之間存在內(nèi)在的動力學行為。

5)通過對k-核解析分析發(fā)現(xiàn)，蘆山M7.0地震和九寨溝M7.0地震發(fā)生前地震網(wǎng)絡的最高核存在顯著的高值異常，但是，汶川地震之前未發(fā)現(xiàn)異常。通過對網(wǎng)絡聚類系數(shù)的分析發(fā)現(xiàn)，可能的原因在于地震網(wǎng)絡層次結構未發(fā)生變化。這與汶川地震的孕震條件、特定的地震動力學過程有一定相關性，值得進一步分析。

雖然本研究采用5km和10km尺度進行網(wǎng)格劃分得到的結果是穩(wěn)定的，但不同地區(qū)網(wǎng)格劃分尺度需要進一步斟酌。此外，汶川地震前最高核變化未出現(xiàn)明顯異常，而李光光等[19]在美國加州地區(qū)的計算結果也顯示并不是所有大地震都能引起地震網(wǎng)絡最高核數(shù)的異常變化。所以，考慮到單參數(shù)對地震網(wǎng)絡特性描述的片面性，難以綜合反映地震網(wǎng)絡的動力學演化規(guī)律，因此在后續(xù)研究中建議考慮多參數(shù)融合分析，從而客觀全面地探究震前地震活動異常特征。