郯廬斷裂帶中南段前兆綜合信息量地震預測方法研究*

耿 杰，王 梅，單連君

(1.山東省地震局，山東濟南250014;2.河北省地震局，河北石家莊050021)

0 引言

郯廬斷裂帶作為中國大陸東部最顯著的深大斷裂，其未來的震情趨勢一直引起我們的高度關注。根據(jù)歷史地震活動與現(xiàn)階段中小震活動特點，該區(qū)存在發(fā)生中強地震的背景。本文研究在對該帶發(fā)生過的中強以上地震前出現(xiàn)的前兆異常進行綜合分析的基礎上，較完整地收集了2007年以來與研究區(qū)相關的省份的數(shù)字化前兆觀測資料及相關臺站和觀測信息 (干擾識別)資料，采用時頻分析方法對研究區(qū)前兆觀測資料進行了異常提取，“邊研究邊應用”地開展了研究區(qū)前兆異常綜合信息量提取方法及異常時間和空間綜合演化特征分析。

1 技術思路和技術途徑

地震預報的研究與實踐表明 (平建軍等，1999，2003a;王煒等，1996)，單項前兆手段的異常時空分布是極其復雜的。沒有一種單項前兆手段的異常在所有地震前都出現(xiàn)，亦沒有在任何一種單項前兆手段的異常變化出現(xiàn)后都有相應的地震發(fā)生，即單項前兆手段與地震的關系具有不確定性。盡管每一個前兆異常都不能反映地震孕育發(fā)生的全過程，但它卻可能反映地震孕育過程的某些信息，這些信息組合起來在一定的時間、空間內構成一個綜合的地震前兆信息場，這個信息場的時、空演變將有可能反映地殼介質破裂的震源物理過程。因此，深入開展地震前兆場的綜合信息量及其時、空動態(tài)演變研究，對于探尋地震孕育發(fā)生的普遍規(guī)律，提取地震預報的判據(jù)和指標，提高地震預報水平具有至關重要的意義。

目前我國的地震監(jiān)測有地形變、地磁、地電、水位、水溫、水化學、地應力等多種學科，各學科甚至同一學科不同觀測方法所測的物理量都不相同，要想將這些量綱各異的各種地震監(jiān)測資料有機地組合起來，進行地震預報綜合分析，首先要設法將每一種地震觀測異常統(tǒng)一轉化為同一信息量值，羅蘭格 (2002)對此進行了研究，提出了地震異常前兆信息量綜合數(shù)學表達式法。經對實際地震觀測資料進行計算，認為該表達式不僅可計算提取各種地震觀測資科異常持續(xù)變化的前兆信息量，還可定量地描述異常結束后其前兆信息的延續(xù)性，從而使其地震前兆信息更加突出合理。這就為我們進一步對比研究地震異常群體前兆信息特征、開展地震前兆綜合信息量的地震預報奠定了基礎。

前兆觀測資料中，往往存在著趨勢上升、下降、周期等非平穩(wěn)變化，這就需要選擇平滑濾波、傅里葉變換、變化率等數(shù)學方法，事先進行消周期、消趨勢項的數(shù)據(jù)預處理工作。地震前兆信息可能存在于不同的頻帶中，既可能存在于高頻段，也可能存在于低頻段，還可能兼而有之。在現(xiàn)行的數(shù)據(jù)預處理工作中，對前兆觀測數(shù)據(jù)或去低頻，或消高頻，通常會導致遺漏一些有用的異常信息，使異常信息不能完全表征出來。近幾年越來越被重視且得到廣泛應用的小波變換和經驗模態(tài)(EMD)分解技術 (程正興，武鐵敦，1994;武安緒等，2006)，為解決這一問題帶來了便利。通過對原始前兆觀測資料進行小波變換或EMD分解，可有效地把各種頻率成份從中分離出來，且分離出來的時頻序列資料數(shù)據(jù)變化平穩(wěn)，非常便于應用地震前兆信息量綜合數(shù)學表達式對其進行無量綱化震兆信息提取。

本文結合地震前兆信息量綜合數(shù)學表達式和小波變換、EMD分解技術，分析處理了研究區(qū)及其附近地區(qū)67個臺點、91個測項的前兆觀測資料，以期為多頻域剖析前兆觀測資料的變化特征、客觀自然地多層次挖掘地震異常、最大程度地展示前兆群體異常在前兆綜合信息量時序和空間上的演化特征，為有效的進行震情預測提供可信度較高的依據(jù)。

2 地震綜合前兆信息量提取方法

對工作區(qū)域每一前兆觀測資料 (地殼形變、地電、地磁、水化、水位、水溫、應力、應變等)首先判別其是否存在非平穩(wěn)變化。若存在，則使用小波變換或EMD分解技術，將其分解為一組多階時頻序列，分析各時頻序列映震能力，并挑出映震效能最優(yōu)的時頻序列，應用地震前兆信息量綜合數(shù)學表達式，提取其無量綱前兆信息量;如不存在，則可根據(jù)地震前兆信息量綜合數(shù)學表達式，直接提取計算其無量綱前兆信息量。地震前兆信息量綜合數(shù)學表達式及其式中各參數(shù)的含義，請參見有關文獻 (平建軍等，1999，2013b，c;賈炯等2010;李廣鑫等，1992)，這里不再贅述。

2.1 單項前兆信息量的提取

數(shù)字化前兆資料多為典型的觀測信號，為非平穩(wěn)過程，具有不穩(wěn)定性、變化快等時頻特點。在多種天然與人為因素 (不包括觀測因素)影響下，實際前兆觀測資料往往存在趨勢和多種周期成份非平穩(wěn)的動態(tài)變化。其中某些周期性變化，可與具有一定周期性的構造活動以及地震孕育發(fā)生過程產生力學耦合作用，這類周期性變化，即為可能的地震異常。這就是說，前兆數(shù)據(jù)序列中僅是某些周期成分變化才與地震的孕育和發(fā)生有關。小波變換和EMD分解技術的特點，正好適合從多種周期變化中提取出與地震有關的周期性異常信息。

2.1.1 小波變換基礎

對于離散序列信號f(x)，在小波函數(shù)ψ(t)∈L2(R)中，尺度因子 (伸縮因子)a和平移因子b(a，b∈R)，也需要離散化，應用離散小波變換(DWT)作為不同頻率的信息識別基礎，即

在計算中，采用a=2k。隨著k的增加，信號從最高頻向低頻分解。當k=0時，信號為采樣頻率，k=1時，將頻率二等分，依此類推。

對于數(shù)字信號可以近似地表示為

2.1.2 小波選擇與頻率分析

在地球物理領域，Morlet小波應用比較廣泛。因此，我們在分處理地震活動與前兆信號時，均采用Morlet小波 (Morlet et al.，1882)。常用的是復值Morlet小波。

2.1.3 不同頻率的信息識別思路

由于數(shù)字化前兆的觀測精度提高，在干擾因素排除后，對趨勢異常與短期異常的識別與排除也是一大問題，而通過小波變換方法，對不同頻率范圍內的信息 (或高頻與低頻信息)進行識別與分解。根據(jù)式 (2)以及小波分解的近似部分f(x)、細節(jié)部分f(x)與頻率的關系，對觀測資料進行近似部分 (低頻)與細節(jié)部分 (高頻)信息分離。根據(jù)觀測量的物理含義 (如水動態(tài)觀測資料中反映固體潮汐的日波、半日波、半月波、月波等信息)，在分析干擾因素的基礎上，從近似部分中確定出趨勢變化信息;從細節(jié)部分識別出短期異常。

2.1.4 小波基的選取

前兆模擬觀測資料的觀測周期主要為日觀測和整點值觀測，數(shù)字化觀測采樣周期主要為分鐘值觀測和整點值觀測。觀測資料中除含有地殼應力應變的信息外，還含有諸如降雨、溫度、開采和儀器故障等不同頻率的干擾成分，但大多數(shù)資料都呈現(xiàn)出規(guī)律性較好的年變化形態(tài)，因此，在小波分析時，選用正則性較好的離散正交小波(Daubechies小波)進行分析。采用db4小波對數(shù)據(jù)進行處理，首先考慮到db4小波在地球物理學研究范圍被廣泛采用;其次考慮到db4小波是緊支撐正交小波，可以使Mallat算法更快捷。它的光滑性也可以更高精度地模擬和分析信號。db4小波隨著分解層級的增加其正則性也增加，它抑制了該多項式信號在零階和一階的部分信號，而僅對該信號的二階部分及噪聲進行分解。因此，分解的細節(jié)信號部分db4中包含了噪聲信號的不規(guī)則性，其余各層細節(jié)中的信號周期性隨層級的增加而增大。另外也考慮到db4小波在時域和頻域局部化方面的強勁性。

本研究對于db4小波分析，將經預處理后的前兆觀測值分解成10個細節(jié) (10階)時頻序列，多頻域分析前兆觀測資料的變化特征和各時頻序列映震能力，客觀自然地多層次挖掘地震異常，在此基礎上選擇最優(yōu)時頻序列，應用地震前兆信息量綜合數(shù)學表達式，提取其無量綱前兆信息量，以最大程度地展示前兆群體異常映震效能，開展震情預測。

對研究區(qū)域前兆觀測資料，統(tǒng)一按上述方法，提取單項前兆信息量后，即可開展如下時間和空間綜合分析。

2.2 地震前兆信息量綜合分析

2.2.1 地震前兆信息量綜合時序分析

由于實際收集到的研究區(qū)域數(shù)字前兆觀測資料的起止時間有長有短，為抑制不同時段因存在的數(shù)字前兆觀測資料數(shù)目的多少而對綜合結果所產生的影響，在對研究區(qū)域進行地震綜合前兆信息量時序分析時，根據(jù)每一數(shù)字前兆觀測資料地震信息量的計算結果，逐時刻地截取并計算研究區(qū)域地震前兆綜合信息量時序值:

式中，S時序為某時刻工作區(qū)域地震綜合前兆信息量時序值;Si為該時刻某前兆觀測資料的前兆信息量值;N為該時刻所有前兆觀測資料總數(shù)。

2.2.2 地震前兆信息量綜合空間分析

根據(jù)研究區(qū)域每一數(shù)字前兆觀測資料地震信息量的計算結果及其臺站的空間經、緯度，再考慮同一觀測臺站不同觀測項目地震信息量間的差異性，以及前兆觀測臺站在空間分布上的不均勻性對地震信息場所產生的影響，按式 (6)逐時刻地截取并繪制研究區(qū)域地震綜合前兆信息量空間分布圖，通過研究中強震前區(qū)域震兆信息空間演化特征，從而建立其中強震空間發(fā)震地點判據(jù)指標。

式中，S窗為某時刻工作區(qū)域地震綜合前兆信息量時序值;Si為該時刻某前兆觀測資料的前兆信息量值;N為該時刻該經緯度節(jié)點掃描范圍內前兆觀測資料項目總數(shù)。

圖1 研究區(qū)前兆信息量分析臺 (點)分布圖Fig.1 Distribution map of the stations(points)for precursory information analysis in the study area

3 研究區(qū)地震前兆綜合信息量的計算及其映震效能分析

3.1 研究區(qū)前兆觀測資料基本情況

本文研究區(qū)域涉及遼寧、河北、北京、天津、山東、河南、江蘇、安徽8個省 (市)(圖1)，在對上述區(qū)域各臺站測項觀測資料進行充分收集、分析和調研基礎上，結合2012年度全國前兆觀測資料預效能評估結果，系統(tǒng)收集處理了各省 (市、區(qū))預報效能較好、干擾因素相對較少、觀測數(shù)據(jù)可信度較高且資料完整性相對較好的67個觀測臺 (點)合計91個測項2007-01-01～2014-06-10資料進行分析計算。

分析計算所使用原始數(shù)據(jù)來源于上述省份及國家臺網(wǎng)中心數(shù)字化前兆臺網(wǎng)Oracle原始和預處理數(shù)據(jù)庫，在此基礎上生成整點值或日值數(shù)據(jù)，作為計算的基礎數(shù)據(jù)。原始觀測數(shù)據(jù)中，由于人為、儀器原因 (標定、校測和故障)等因素，經常會產生觀測值的非正常變化，直接影響數(shù)據(jù)的分析，對于這些非正常變化、非客觀的資料動態(tài)，不同形態(tài)特征要進行具體分析，針對經常出現(xiàn)的非正常變化情況，在計算前分別對其進行了預處理，本次研究的主要預處理方法有:

缺值:由于儀器原因，在某些時間段內缺值，數(shù)據(jù)中以“999999”表示。對于缺值這種情況，依據(jù)一定的方法進行了插值。常用的插值法有線性插值法、拋物線插值法、拉格朗日插值法和多元線性回歸擬合外推插值法等。

突降:數(shù)據(jù)中某些點由于儀器等干擾原因造成突降。這種情況采取了根據(jù)突升前后值，將突降后數(shù)據(jù)接回突降以前的處理。

突跳:由于已知的干擾方面原因，數(shù)據(jù)出現(xiàn)幾個點高于或低于周圍點的突跳變化情況，根據(jù)突跳點兩端數(shù)據(jù)的平均值，插入突跳處，以消除突跳點。

3.2 單項前兆信息量的計算與分析

對于研究區(qū)域每一前兆觀測資料 (形變、電磁、流體等)，首先判別其是否存在非平穩(wěn)變化。若存在非平穩(wěn)變化，則使用小波變換或 EMD分解技術，將其分解為一組10階的時頻序列，分析各時頻序列映震能力，并選擇出映震效能最優(yōu)的時頻序列，應用地震前兆信息量綜合數(shù)學表達式，提取其無量綱前兆信息量。

對研究地區(qū)前兆觀測資料進行小波或經驗模態(tài) (EMD)分解，結果顯示該方法對非平穩(wěn)觀測數(shù)據(jù)的異常信息具有較強的挖掘和識別能力，部分原始數(shù)據(jù)曲線上異常特征不明顯的測項，經識別后，顯示有明顯的異常信息存在。

3.3 地震前兆信息量綜合時序和空間分析

對研究區(qū)域前兆觀測資料，統(tǒng)一按上述方法提取前兆信息量后，即可開展時、空綜合分析。

3.3.1 地震前兆信息量綜合時序

2007-01-01～2014-06-10期間，研究區(qū)及其附近地區(qū)發(fā)生ML≥5.0地震4次，顯著震群事件1次，我們對上述地震進行地震前信息量綜合時序回溯性檢驗。對經過小波變換或EMD分解后的單項前兆異常信息量提取結果按 (5)式進行綜合信息量綜合時序分析，結果見表1、圖2。

除安徽安慶地震前信息量綜合時序曲線未出現(xiàn)異常外，其余4次地震前均顯示出現(xiàn)明顯的高值異常出現(xiàn)，異常出現(xiàn)在震前1～105 d，異常幅度最大為0.158。上述研究區(qū)內發(fā)生的5次地震均發(fā)生在前兆綜合信息量時序曲線的高值異常出現(xiàn)后，由于研究時段內所發(fā)生的地震震級差異不大，在前兆綜合信息量時序曲線上顯示，異常幅度和異常持續(xù)時間差異不大。

表1 研究區(qū)及其附近地區(qū)M L≥5.0地震前兆綜合信息量時序異常特征 (S時序≥0.1)(2007-01-01～2014-06-10)Tab.1 The time series anomaly characteristics(S time≥0.1)of earthquake precursor information with M L≥5.0 earthquakes in the study area and its vicinity from Jan.1，2007 to June.10，2014

圖2 研究區(qū)前兆綜合信息量時序曲線(2012-01～2014-09)Fig.2 The time sequence curve of earthquake precursor comprehensive information in the study area from Jan.2012 to Sep.2014

圖3 1995年蒼山M S5.2地震前后前兆綜合信息量時序曲線Fig.3 The time sequence curve of earthquake precursor comprehensive information Before and after the Cangshan M S5.2 earthquake in 1995

1995年9月20日山東蒼山MS5.2地震是發(fā)生在郯廬斷裂帶上的顯著地震事件之一，蒼山地震前兆綜合信息量時序值計算結果顯示 (圖3)，1995年8月10日和8月20日分別出現(xiàn) S時序為0.28和0.23的高值異常，震后時序曲線迅速恢復到正常背景值波動。

2015年度山東省地震趨勢會商分析，2012-01-01～2014-09-30對山東地區(qū)前兆測項觀測資料進行時序前兆綜合信息量計算結果顯示，2013年11月23日萊州ML5.0地震前，于2013-08-10～2013-09-10出現(xiàn)S時序≥0.2的高值異常時段，并于2013年9月10日出現(xiàn)最高值達0.301;2014-01-01～2014-09-30 S時序計算結果，僅在2014年7月31日出現(xiàn)0.216的異常值，8月恢復至正常波動，再無高值異常顯示，至2015年3月，山東地區(qū)無ML≥5.0地震發(fā)生。

3.3.2 地震前兆信息量綜合空間分析

對2007-01-01～2014-06-10時段研究區(qū)前兆測項資料經過小波或EMD分解后的單項前兆異常信息量提取結果按式 (6)進行綜合信息量空間分析，結果見表2。

對前述5次地震進行地震前信息量綜合空間回溯性檢驗表明，除安徽安慶地震前信息量綜合空間掃描未出現(xiàn)S窗≥0.2的異常區(qū)外，其余4次地震前均顯示有明顯的異常區(qū)出現(xiàn) (圖4)。

對2007-01-01～2014-06-10時段研究區(qū)綜合信息量空間掃描結果，2007-01-01～2014-05-31無S窗≥0.2的異常區(qū)出現(xiàn)。研究區(qū)在2012-06-10出現(xiàn)S窗≥0.2的異常區(qū)后，2012-06-30在該異常區(qū)西側又出現(xiàn)一個S窗≥0.2的異常區(qū)，之后相鄰的兩個異常區(qū)逐漸擴展，至2012-07-10達最大，之后逐漸收縮。2012年7月20日高郵ML5.3地震發(fā)生在最初出現(xiàn)S窗≥0.2異常區(qū)邊緣，S窗最大值0.98。

表2 研究區(qū)及其附近地區(qū)M L≥5.0地震前綜合信息量 (S窗≥0.2)空間演化特征 (2007-01-01～2014-06-10)Tab.2 Spatial evolution characteristics of earthquake precursor comprehensive information(S time≥0.20)before M L≥5.0 earthquakes in the study area and its vicinity from Jan 1.2007 to Jun.10，2014

圖4 2013年11月23日山東萊州地震前地震前兆綜合信息量空間演化特征(a)2013-06-30;(b)2013-07-10;(c)2013-07-31;(c)2013-08-10;(d)2013-09-10Fig.4 The spatial evolution map of earthquake precursor comprehensive information before Shandong Laizhou M L5.0 earthquake on Jan.1，2013 in Shandong Province

2013年1月23日燈塔ML5.1地震的回溯性檢驗顯示:2012-11-10研究區(qū)開始出現(xiàn)S窗≥0.2的異常區(qū)，異常區(qū)逐漸擴大，至2012-12-10達到最大后逐漸收縮，至2013-01-20后異常區(qū)域消失。2013年1月23日燈塔ML5.1地震發(fā)生在最初出現(xiàn)S窗≥0.2異常區(qū)邊緣，S窗最大值0.999 6。

2013-01-20～2013-06-20期間，研究區(qū)內無S窗≥0.2異常區(qū)出現(xiàn)。2013-06-30在沂沭斷裂帶南段西側首先出現(xiàn)S窗≥0.2的異常區(qū)，該異常區(qū)逐漸擴展，至2013-07-10達最大，并且在萊州ML5.0地震區(qū)域附近出現(xiàn)異常區(qū)，2013-07-31后萊州地震附近區(qū)域的異常區(qū)逐漸收縮，至2013-08-10完全消失，與此同時，前述沂沭斷裂帶南段西側的異常區(qū)也經歷了逐漸縮小的過程。2013年11月23日萊州ML5.0地震發(fā)生在最初出現(xiàn) S窗≥0.2異常區(qū)邊緣，S窗最大值0.997(圖4)。

2014年4月1日發(fā)生乳山ML4.6震群，該震群是2014年度山東地區(qū)最為顯著的地震事件，在華北地區(qū)也十分突出。研究區(qū)空間綜合信息量掃描結果顯示，2013-12-10開始出現(xiàn)S窗≥0.2的異常區(qū)，之后異常區(qū)逐漸擴展，至2014-01-20達最大后逐漸縮小 (2014-02-10)。2014-02-20在異常區(qū)西側出現(xiàn)異常區(qū)，2014-03-10乳山震群附近的異常區(qū)消失。乳山震群發(fā)生在最初出現(xiàn)S窗≥0.2異常區(qū)邊緣，S窗最大值0.997。

4 分析結論

(1)小波和EMD分解技術的應用，可有效地把各種頻率成份從原始前兆觀測資料中分離出來，多頻率域的充分展示前兆觀測資料的異常信息，增加了對地震異常的辨識和挖掘力度。將地震前兆信息量綜合數(shù)學表達式和小波、EMD分解技術結合在一起，分析處理前兆觀測資料，為多頻域剖析前兆觀測資料的變化特征，客觀自然地多層次挖掘地震異常、最大程度地展示前兆群體映震效能提供了新的途徑。

(2)分析研究區(qū)地震綜合前兆信息量時序曲線，在ML≥5.0地震前會出現(xiàn)不同程度的高值異常 (S時序≥0.10)，一般高值異常出現(xiàn)后1～6個月內研究區(qū)有發(fā)生ML≥5.0地震的可能。

(3)對研究區(qū)地震前兆信息量綜合空間掃描的結果顯示，未來的震中區(qū)多位于最早出現(xiàn)S窗≥0.2異常區(qū)的邊緣區(qū)域，但異常區(qū)域的遷移規(guī)律性特征不明顯，這可能與資料分析時段研究區(qū)內發(fā)生地震的震級較小及前兆臺點的空間分布不勻有一定的關系。前兆臺點的空間分布不勻某種程度上導致前兆異常分布不勻，也會導致等值線連線在一定程度上存在失真現(xiàn)象。

程正興，武鐵敦.1994.小波的發(fā)展與應用［J］.微機發(fā)展，(5):8-10.

賈炯，平建軍，刁桂苓，等.2010.華東地區(qū)地震綜合前兆信息量及其映震能力分析［J］.華北地震科學，28(2):48-52.

李廣鑫，程式，李振興.1992.綜合積分預報方法在四川地區(qū)的應用研究［J］.四川地震，(3):41-46.

羅蘭格.2002.我國地震綜合預報方法研究的回顧與展望［J］.華北地震科學，20(4):1-l8.

平建軍，王秀英.1999.地震前兆信息量的研究(二)地震前兆信息量的綜合數(shù)學表達式［J］.華北地震科學，17(3):8-l4.

平建軍，張永仙，單連君，等.2003a.華北地區(qū)地震短期異常分析及其前兆信息量的提取與計算［J］.華北地震科學，21(3):1-7.

平建軍，張永仙，單連君，等.2013b.震災防御技術，8(4):397-407.

平建軍，張永仙，單連君.2013c.華北地區(qū)地震短期異常分析及其前兆信息量的提取與計算［J］.華北地震科學，21(3):1-7.

王煒，吳耿鋒，黃冰樹，等.1996.基于模糊神經網(wǎng)絡和符號的地震預報專家系統(tǒng) NGESEP［J］.中國地震，12(4):339-346.

武安緒，林向東，穆會泳，等.2006.EMD新技術在數(shù)字波形預處理中的初步應用［J］.華南地震，26(1):l33-138.

Morlet J.G.，Arens G.，F(xiàn)ourgeau E..1982.Wave propagation and sampling theory:complex signal and scattering in muti 2 layered media［J］.Geophysics，47(2):203-211.