利用概率密度分布提取流體觀測資料中的高頻異常信息
——以2008年汶川8.0級地震為例

孫小龍， 王廣才， 晏銳

1 中國地質(zhì)大學(北京)， 北京　100083 2 中國地震局地殼應力研究所， 北京　100085 3 中國地震臺網(wǎng)中心，北京　100045

利用概率密度分布提取流體觀測資料中的高頻異常信息

——以2008年汶川8.0級地震為例

孫小龍1,2， 王廣才1， 晏銳3*

1 中國地質(zhì)大學(北京)， 北京100083 2 中國地震局地殼應力研究所， 北京100085 3 中國地震臺網(wǎng)中心，北京100045

摘要隨著地下流體觀測技術(shù)的提高，尤其是地下流體觀測數(shù)字化改造以后，觀測資料的采樣頻率明顯提高，這些高頻采樣的觀測資料中蘊含著豐富的構(gòu)造信息，如何從分鐘值甚至更高頻率采樣的觀測資料中提取有用的異常信息，是目前從事地震地下流體資料分析人員最為關(guān)注的科學問題之一.本文引入概率密度分布法，分析了2008年汶川8.0級地震前南北地震帶及其附近區(qū)域72個測點的數(shù)字化水位和水溫分鐘值采樣高頻觀測資料，結(jié)果顯示：汶川8.0級地震前有16個測點水位和14個測點水溫出現(xiàn)高頻異常信息，出現(xiàn)高頻信息異常的觀測點多集中在滇西南構(gòu)造帶，異常出現(xiàn)的時間呈現(xiàn)出由南向北推移的特征.據(jù)此認為概率密度分布法在流體資料的高頻異常信息提取方面具有一定的可靠性與適用性，可為數(shù)字化高頻觀測資料異常提取提供借鑒.

關(guān)鍵詞流體資料； 概率密度分布； 高頻異常信息； 汶川8.0級地震

1引言

21世紀初，針對觀測井孔和觀測儀器老化、觀測技術(shù)水平跟不上技術(shù)現(xiàn)代化發(fā)展要求的現(xiàn)狀，中國地震局在“九五”、“十五”期間實施了前兆臺站觀測技術(shù)改造工程，對全國前兆觀測井進行了數(shù)字化、網(wǎng)絡化改造.經(jīng)歷了近10年的實踐與資料積累，數(shù)字化觀測的優(yōu)勢逐漸被行業(yè)內(nèi)同行所認可，與之前的模擬觀測相比，數(shù)字化觀測具有數(shù)據(jù)傳輸速度快、數(shù)據(jù)信息量豐富、人為誤差少等優(yōu)點，從而使地下水位的高頻、短周期信息大量增加，為捕捉地震孕育及發(fā)生過程中的前兆異常信息提供了有利條件.流體觀測技術(shù)經(jīng)歷了數(shù)字化改造后，觀測設備采樣率明顯提高，除了之前常采用的月均值、旬均值和日均值外，還出現(xiàn)了整點值、分鐘值、秒值以及更高頻的數(shù)據(jù).

但是，在日常的數(shù)據(jù)處理與震情跟蹤工作中，業(yè)內(nèi)普遍沿用先前針對模擬觀測所提出的各類分析方法(張素欣等，2002；楊明波等，2009)，數(shù)據(jù)類型的選取多采用月均值、日均值和整點值，數(shù)字化觀測資料中的分鐘值、秒值以及更高頻的數(shù)據(jù)信息，僅限于應用在流體對遠場地震的同震響應分析(劉耀煒等，2005；孫小龍和劉耀煒，2008；曹玲玲和高安泰，2010； Sun and Liu,2012)，更普遍的現(xiàn)象是在某些時候，工作人員為了異常分析需要，還須將數(shù)字化觀測資料中的分鐘值、整點值等轉(zhuǎn)化成日均值、月均值進行異常提取，變相增加了分析人員的工作量.出現(xiàn)以上現(xiàn)象的原因，并非是數(shù)字化觀測技術(shù)或觀測資料本身存在某種缺陷，而是由于在數(shù)字化改造實施后，尚缺乏與之相適應的分析處理方法.因此，在數(shù)字化觀測技術(shù)改造完成后，擺在地震監(jiān)測與預報人員面前的首要問題是盡快提出與數(shù)字化觀測資料相匹配的數(shù)據(jù)處理方法，即數(shù)字化觀測資料中高頻信息的異常提取方法研究.

近幾年來，科研人員針對數(shù)字化觀測資料的特征與其所面臨的科學問題，探索性地提出了一些適合數(shù)字化資料的數(shù)據(jù)處理與異常分析方法.例如，楊叢杰等(2005)運用小波分析方法對1990 年常熟、1995 年蒼山和1996 年南黃海三個中強地震前江蘇地區(qū)井水位固體潮的變化特征進行了研究，發(fā)現(xiàn)井水位M2波潮汐因子在地震前幾個月幾乎都出現(xiàn)一個幅度較大、周期為半月或一個月左右的異常信號，表明小波分析方法在處理和分析井水位潮汐資料方面可能是一種有效的方法.邱澤華等(2010，2012)應用小波分解-超限率法處理并分析了汶川8.0級地震前姑咱臺和寧陜臺鉆孔應變高頻數(shù)字化資料在震前出現(xiàn)短周期“毛刺”的異?，F(xiàn)象，分析認為這種高頻的異常是震前構(gòu)造運動的表現(xiàn).劉琦和張晶(2011)引入時頻分析中的S變換法，對姑咱臺四分量鉆孔應變觀測曲線中出現(xiàn)的“壓脈沖”和“潮汐畸變”異常信號進行了分析，結(jié)果表明在背景信號的基礎(chǔ)上，周期約為10～60 min的信號在汶川地震前開始大量出現(xiàn)，震后逐漸衰減.晏銳等(2011)基于臨界慢化概念，運用表征臨界慢化現(xiàn)象的自相關(guān)系數(shù)和方差法，分析了汶川8.0級地震前氡濃度觀測資料，結(jié)果表明不同臺站的水氡濃度震前都存在明顯的臨界慢化現(xiàn)象.Manshour等(2009，2010)曾用概率密度分布法成功提取了全球12次強震發(fā)生前高頻垂直速度觀測數(shù)據(jù)中的異?，F(xiàn)象，結(jié)果表明該方法在提取高頻異常信息時效果明顯.

由此可見，數(shù)字化觀測資料的高頻成份中隱含著大量的前兆異常信息，運用適當?shù)臄?shù)據(jù)處理方法，是可以識別出這些前兆異常的.

2概率密度分布法識別異常信息

隨著地震觀測數(shù)據(jù)采樣率的不斷提高，其數(shù)據(jù)量也急劇增大.如何從大量的觀測數(shù)據(jù)中提取出分析所需的有用信息，是當前數(shù)字化觀測資料分析研究的重點.目前，在提取數(shù)字化觀測資料中的高頻信息時，普遍采用高通濾波法(邱澤華等，2010，2012)、小波分析法(Daubechies,1988；晏銳等，2011；邱澤華等，2012)、多項式擬合法(Manshour et al.,2009，2010)、經(jīng)驗模態(tài)分解法(Huang et al., 1998；孫小龍等，2011)等，這些方法既能有效地提取出觀測數(shù)據(jù)的低頻趨勢變化項，也能提取出高頻的短周期信息.

本文在地震流體數(shù)字化觀測資料的高頻信息成份提取過程中，采用了三階多項式法，該方法在日常的數(shù)據(jù)處理中多用來提取流體觀測資料中的長趨勢變化，而本文主要利用原始數(shù)據(jù)在消除趨勢變化后的剩余信息，即短周期的高頻信息成份.三階多項式法主要是利用多項式擬合法對原始數(shù)據(jù)進行趨勢擬合，得到一條最接近原始曲線趨勢變化的三階多項式變化曲線，利用原始曲線數(shù)據(jù)減去趨勢變化曲線數(shù)據(jù)，得到原始數(shù)據(jù)中的高頻信息成份.

對一組時序數(shù)據(jù)D(t)，以s為窗長進行三階多項式擬合，得到趨勢項擬合曲線為

T(t)=at3+bt2+ct+d，

(1)

其中，a、b、c和d均為擬合系數(shù)，t為時間.高頻信息成份Z(t)=D(t)-T(t)，給定窗長s后依次滑動進行擬合求參數(shù)，并最后得到原始數(shù)據(jù)中的高頻信息Z(t).

圖1為四川南溪臺水位分鐘值曲線及其利用三階多項式(窗長600 min)提取的高頻信息，其中圖1a和圖1b分別為一個窗長的原始觀測曲線和剔除趨勢后的高頻信息曲線(a中虛線為依據(jù)原始觀測值利用三階多項式擬合得到的趨勢變化線)，圖1c和圖1d分別為2007年9月1日至2008年5月11日的原始觀測曲線和剔除趨勢后的高頻信息曲線.從圖可以看出：

(1) 在一定窗長范圍內(nèi)，三階多項式方法基本上能擬合得到觀測值的趨勢變化線(圖1a)，進而可去除趨勢得到觀測曲線的高頻信息.如果觀測曲線中有明顯的潮汐或氣壓變化形態(tài)，如部分水位觀測曲線，該方法一定程度上能消除掉觀測曲線中固體潮或氣壓的半日潮變化(圖1b)，若想達到較好的擬合效果，窗長應小于半天.

(2) 從年尺度和月尺度的曲線對比來看，三階多項式方法能很好地剔除掉原始觀測曲線(圖1c)中的長期趨勢變化，得到原始觀測曲線中的高頻信息(圖1d).

(3) 圖1d中剔除掉趨勢變化后的高頻信息曲線，可看出類似于“水位固體潮半月潮”的波動變化.從圖1a可看出，雖然三階多項式法可很好地擬合得到觀測值的趨勢變化，但從局部來看，也有不吻合的部分(圖1b中3∶00—5∶00時段)，導致其去趨勢殘差(高頻信息)值較大.另外，水位半月潮是由于水位固體潮響應幅度和相位的規(guī)律性變化引起，在利用三階多項式進行擬合去趨勢時，這種規(guī)律性的變化也會直接影響到擬合殘差的波動變化.因此，結(jié)合圖1a和圖1b，這種類似于“水位固體潮半月潮”的信息，并非真正的水位半月潮信息，而是受其影響的擬合殘差的波動信息.

圖1　四川南溪臺水位分鐘值觀測曲線及其高頻信息Fig.1　Curves of water level recorded (a—b) and high-frequency information (c—d) at Nanxi station, Sichuan

地震流體觀測資料，其長周期的趨勢變化受構(gòu)造環(huán)境和水文地質(zhì)條件的影響較大(郭增建等，1974；楊明波等，2009)，而其短周期的高頻信息在無其他干擾因素的條件下，多表現(xiàn)為一種服從正態(tài)分布的隨機信號.從統(tǒng)計的角度分析，在去除趨勢變化后，且無其他影響因素存在的條件下，觀測資料變化值服從正態(tài)分布X(t)～N(0,σ).但是，如果存在某種干擾(如構(gòu)造活動或應力場變化)，σ值就會發(fā)生波動，并且σ值的波動會依干擾強度的變化而發(fā)生變化.從能量傳遞的角度考慮，其波動值服從對數(shù)正態(tài)分布lnσ～N(0,λ)(Castaing et al., 1990)，Z(t)值的概率密度分布可表示為

(2)

式中，σ和λ分別為Z(t)和lnσ的標準差，σ0為σ的數(shù)學期望值.

如圖2所示，分別服從正態(tài)分布和對數(shù)正態(tài)分布的兩組信號X(圖2a)和Y(圖2b)，其概率密度分布分別為圖2d和圖2e，由二者乘積合成的信號Z=XY(圖2c)的概率密度分布為圖2f.由圖可以看出，單純服從正態(tài)分布的原始信號X，在加入干擾信號Y后，其合成信號Z的概率密度分布P(Z)與X的概率密度分布P(X)形態(tài)有了明顯的不同，并且這種形態(tài)上的不同會依據(jù)干擾信號的強度而發(fā)生變化.如圖3所示，同一個原始信號X，在不同的干擾強度(干擾強度決定了λ值的大小)下，其概率密度分布表現(xiàn)出了明顯的不同，干擾強度越大，其合成信號的概率密度分布P(Z)與正態(tài)分布P(X)在形態(tài)上的偏離程度越大.

因此，在理想狀態(tài)下觀測資料中含有的高頻信息其概率密度分布應符合正態(tài)分布N(0,σ)，高頻信息的標準差σ相對穩(wěn)定，干擾信息的λ值接近于零.

圖2　隨機信號及其概率密度分布(a) 服從正態(tài)分布的隨機信號； (b) 服從對數(shù)正態(tài)分布的隨機信號； (c) 合成信號； (d) 正態(tài)分布的概率密度； (e) 對數(shù)正態(tài)分布的概率密度； (f) 合成信號的概率密度.Fig.2　Random signals and their probability density distribution(a) Normal distribution of random signals； (b) Logarithm normal distribution of signals； (c) Synthetic signals； (d) Normal distribution of probability density； (e) Logarithm normal distribution of probability density； (f) Probability density of synthetic signals.

圖3　隨機信號的疊加及其概率密度分布Fig.3　Stack of random signals and their probability density distributions

如果存在干擾因素，如地震孕育期間構(gòu)造應力狀態(tài)的變化，σ值就會有發(fā)生變化，那么觀測資料中含有的高頻信息其概率密度分布就會偏離正態(tài)分布，干擾信息的λ值就會大于零.對于一組類似于圖1中所提取的高頻信號Z(t)，由式(2)可知， Castaing等人(1990)將其概率密度分布的數(shù)學表達式表示為

(3)

式中，Z表示觀測數(shù)據(jù)，P為該數(shù)據(jù)的概率密度分布，F(xiàn)為符合正態(tài)分布的隨機信號，其標準差為σ，而G為未知的干擾信號，其標準差為λ.F和G的數(shù)學表達式分別為

(4)

(5)

實際應用時，可先統(tǒng)計得出原始數(shù)據(jù)的概率密度分布，然后依式(3)求得干擾G的標準差λ，一般用最大似然法或最小二乘法擬合求得.為得到最佳擬合值，用表達式為

(6)

式中，P(z)表示觀測數(shù)據(jù)的概率密度分布，其標準差為σ，PCastaing為利用式(3)計算得到的概率密度分布，σCastaing為其相應的標準差.

原始數(shù)據(jù)的概率密度分布和式(3)計算得到的概率密度分布，依據(jù)式(6)可求得二者之間的相似程度，當式(6)中的2為最小時，表示二者相似程度最高，進而可得到干擾項G的標準差λ.圖4為利用上述方法統(tǒng)計并計算的四川南溪臺水位分鐘值不同時段的高頻信息概率密度分布，可以看出水位高頻信息的概率密度分布值，在幾乎不存在干擾信號的時段(2007-09-12—2007-10-12)能很好地吻全正態(tài)分布(即λ值接近于0，λ=0.03)，而在臨近地震發(fā)生時段(2008-04-12—2008-05-12)，明顯偏離正態(tài)分布(即λ=0.32)，說明此時段存在某種干擾信號，這種干擾有可能與周邊構(gòu)造活動或地震孕育有關(guān).

圖5為利用以上方法處理的2007年9月12日—2008年5月11日四川南溪臺水位分鐘值高頻信息干擾信號強度時序值曲線(30天窗長，2日滑動，短橫線表示擬合誤差)，處理結(jié)果顯示在2008年1月上旬前，該臺水位受干擾強度很小(λ值接近于0)，而在2008年1月上旬后干擾強度逐漸增強，直到2008年5月12日汶川8.0級地震發(fā)生，截止地震發(fā)生前λ2值超過0.6.這種干擾信號的增強，預示著觀測點周邊構(gòu)造活動或應力積累水平的增加.由此可見，流體觀測資料高頻信息中隱含著與構(gòu)造活動相關(guān)的前兆信息，利用概率密度分析方法在一定程度上可以識別出這種前兆信息.

圖4　四川南溪臺水位分鐘值高頻信息概率密度分布(不同時段)Fig.4　Probability density distributions of high-frequency information for water-level minute values at Nanxi station of Sichuan (in different time intervals)

圖5　四川南溪臺水位分鐘值高頻信息及其概率密度Fig.5　High-frequency information from water-level minute values at Nanxi station of Sichuan and its probability density

圖6　汶川8.0級地震前流體高頻信息異常曲線Fig.6　Curves of high-frequency fluid anomaly information before Wenchuan MS8.0 earthquake

3汶川8.0級地震前流體異常提取

本文利用上述概率密度分析方法，處理了中國大陸地區(qū)南北地震帶上72個流體觀測點的水位、水溫觀測數(shù)據(jù)(共計122條數(shù)據(jù))，數(shù)據(jù)時間段為2008年汶川8.0級震前一年(2007年9月12日—2008年5月11日).處理結(jié)果顯示，共有30條數(shù)據(jù)在地震前出現(xiàn)了異常變化，即干擾信號的λ值明顯增強.異常測項中水位出現(xiàn)異常的有16項，水溫有14項，個別觀測點水位和水溫均出現(xiàn)了概率密度分布異常(如云南開遠觀測點).圖6為其中部分異常變化的λ2值曲線，可以看出水位、水溫的高頻信息異常多出現(xiàn)在地震前3個月內(nèi)，并且這種異常多以震前上升和高值為主.另外，各觀測點異常出現(xiàn)的時間和λ2最高值不盡相同，在一定程度上反映了觀測點周邊應力水平受強震孕育影響的程度也各不相同.部分測點的λ2值在地震前出現(xiàn)了起伏變化(如圖6中的云南姚安水位)，這也反映了局部構(gòu)造活動或區(qū)域應力場在受強震孕育影響后的反復調(diào)整過程.

圖7a所示為本文所選南北地震帶流體觀測點中出現(xiàn)高頻信息異常變化和未出現(xiàn)異常變化測點的空間分布圖，由圖可以看出，出現(xiàn)異常變化的測點多分布在滇西南構(gòu)造帶，震源區(qū)以北出現(xiàn)異常的測點較少.圖7b所示為出現(xiàn)異常測點的λ2高值的大小及其出現(xiàn)時間的空間分布圖，雖然這些異常點的λ2高值在空間分布上無明顯規(guī)律，但其出現(xiàn)時間有明顯的由南向北推移的特征：如圖7b所示，較早出現(xiàn)的異常多集中在南北地震帶南段的云南南部區(qū)域，其異常出現(xiàn)時間多在地震前90天以上(藍色虛線以下，截止2008年5月12日12時)，之后出現(xiàn)的異常多集中在云南北部地區(qū)，其異常出現(xiàn)時間多為60—90天(藍色虛線與紅色虛線之間)，最后出現(xiàn)的異常點多集中在靠近震中的四川地區(qū)，其異常出現(xiàn)時間在震前30天左右(紅色虛線以上).

利用概率密度分布方法提取的汶川8.0級地震前流體高頻信息異常變化特征顯示： (1) 地震前出現(xiàn)高頻信息異常的觀測點多集中在滇西南構(gòu)造帶； (2) 異常出現(xiàn)的時間有由南向北推移的特征； (3) 異常的λ2值在空間分布上無明顯的變化規(guī)律.由概率密度分布法的分析原理可知，流體觀測資料高頻信息中的λ值的變化在一定程度上反映了觀測點受周邊構(gòu)造應力調(diào)整的作用，結(jié)合以上異常變化特征，分析認為在汶川8.0級孕育過程中，滇西南構(gòu)造帶的構(gòu)造應力作用或構(gòu)造活動受孕震過程的影響程度較高，并且這種影響有由南向北推進的過程，這種影響是否于大區(qū)域的構(gòu)造動力學背景有關(guān)，值得進一步深入分析.另外，異常的λ2值大小無明顯空間分布的特征表明，λ2值的大小可能與觀測點本身對構(gòu)造應力作用的敏感程度有關(guān)，如不同的構(gòu)造位置或不同的水文地質(zhì)條件等.

4結(jié)論與討論

作為一種統(tǒng)計學方法，概率密度分布法曾多次用在地震前兆異常識別(Manshour et al., 2009，2010)、太陽風擾動分析(Carbone et al., 1996； Burlaga and Lazarus, 2000；Sorriso et al., 2001)、外匯市場分析(Ghashghaie et al., 1996)以及人體心率異常檢測(Kiyono et al., 2005)等領(lǐng)域，本文將其引入地震地下流體數(shù)字化資料分析中，期望能為資料分析和異常識別提供一種新的手段.

當觀測資料中包含的地震異常比較顯著時，可從原始觀測曲線中直接識別(何案華等，2012)，但也有一些異常是隱含在觀測資料中的，需要用一些數(shù)學處理方法來識別(邱澤華等，2010,2012；劉琦和張晶，2011；晏銳等，2011).概率密度分布法識別地震異常前，需先進行觀測數(shù)據(jù)的高頻信息提取，提取方法也有多種，如三階多項式擬合去趨勢法、小波分析法等、經(jīng)驗模態(tài)法等，提取效果均可滿足實際工作需要.圖8為分別利用三階多項式和小波分析方法提取的四川南溪臺水位分鐘值觀測數(shù)據(jù)中的高頻信息及其概率密度分布曲線，從圖可以看出，盡管兩種方法的處理結(jié)果稍有不同，但利用其處理出的高頻信息均可有效地識別出地震前的異常信號，即地震前概率密度分布曲線的λ2值出現(xiàn)了明顯的高值變化.

本文利用概率密度分布方法提取的汶川8.0級地震前流體高頻信息異常變化結(jié)果表明，地震前存在一批流體高頻信息異常，這些異常多集中在滇西南構(gòu)造帶，并且異常出現(xiàn)的時間有由南向北推移的特征，這些異常有可能與大區(qū)域的構(gòu)造動力學背景有關(guān).

圖7　汶川8.0級地震前流體高頻信息異?？臻g分布Fig.7　Spatial distributions of high-frequency fluid anomaly information before Wenchuan MS8.0 earthquake

圖8　三階多項式與小波分解法對比Fig.8　Comparison of three-order polynomial and wavelet decomposition methods

從板塊運動的角度分析，汶川大地震是印度板塊推擠歐亞板塊的結(jié)果，其最根本動力來源是青藏高原和華南地塊之間的相對運動在斷裂帶上產(chǎn)生的能量積累和釋放(Wang et al., 2001；Zhang et al., 2004；張培震等，2008).這種大尺度范圍的構(gòu)造動力學作用，也會影響到地殼內(nèi)部巖體孔隙介質(zhì)的變形，進而引起地下水流狀態(tài)(如流速、流向、水溫等)的變化.結(jié)合概率密度分布法自身的原理，當流體敏感觀測點所在區(qū)域存在明顯的構(gòu)造動力作用時，原本服從正態(tài)分布規(guī)律的高頻信息，就會偏離正態(tài)分布，指示干擾強度的λ值就會出現(xiàn)高值異常變化.

從異常的空間分布來看，雖然震源區(qū)以北的西北地區(qū)和以南的西南地區(qū)背景觀測點數(shù)量相當，但絕大多數(shù)異常集中分布于震源區(qū)以南的滇西南構(gòu)造帶，而在震源區(qū)以北的區(qū)域分布很少.張培震等的研究表明，汶川8.0級特大地震孕育和發(fā)生是不同力學性質(zhì)、不同變形方式的多個地質(zhì)單元相互作用的結(jié)果，涉及3個地質(zhì)單元：川西高原由于地殼結(jié)構(gòu)的軟弱而發(fā)生強烈震前變形，是孕震的變形單元；龍門山斷裂帶的組成物質(zhì)強度很高、斷裂產(chǎn)狀不利于滑動，震前的變形很緩慢，是孕震的閉鎖單元；四川盆地由于剛性好不易變形而對青藏高原的向東擴展起著阻擋作用，是孕震的支撐單元.震前變形主要發(fā)生在變形單元(川西高原)，閉鎖單元(龍門山斷裂帶)只發(fā)生緩慢的變形和運動，支撐單元由于起著阻擋作用而與閉鎖單元之間幾乎不發(fā)生相對運動或變形(張培震等，2009).同時，川滇地區(qū)GPS現(xiàn)代地殼運動速度場和應變率場分析結(jié)果也表明，汶川8.0級地震震源區(qū)以南的滇西南構(gòu)造帶的地殼運動速度和應變率明顯高于震源以北區(qū)域(張培震等，2009； Wu et al., 2011).可見，流體高頻信息異?？臻g分布與周邊地殼運動場具有很好的一致性，因此認為二者之間有一定的關(guān)聯(lián)性.

從異常出現(xiàn)的時間來看，異常最先出現(xiàn)在距離震源區(qū)較遠的滇南地區(qū)，之后是川滇交界及其附近區(qū)域，最后向震源區(qū)推進.參考邱澤華等(2010)對異常形成機理的分析與Lockner等(1991)、Lockner(1996)的實驗結(jié)果，在斷層成核之前，聲發(fā)射源在一個相當廣大的區(qū)域上發(fā)展，伴隨成核過程開始，聲發(fā)射源向斷層部位集中，周圍遠處的聲發(fā)射源趨于消失.這種觀點可以解釋高頻信息的異常演化特征：汶川8.0級地震發(fā)生前，也即斷層破裂前，距離震中較遠的區(qū)域異常開始出現(xiàn)，之后逐步向震中位置集中.與本文研究結(jié)果類似，晏銳等的結(jié)果表明汶川震前水氡濃度都存在明顯的臨界慢化現(xiàn)象(晏銳等，2011).

圖9為云南高大井水溫高頻信息概率密度自觀測以來至今λ2值時序曲線與周邊地震的對應圖，從圖中可以看出在周邊發(fā)生的多次中強地震前，該井水溫概率密度λ2值均出現(xiàn)了不同幅度的升高(λ2值高于0.02).對比異常信息表(表1)，認為異常幅度值與地震震級存在一定的正相關(guān)，即震級大越大異常幅度越高.同時，異常幅度值與震中距也有一定的關(guān)系，除了在周邊中強地震前λ2值出現(xiàn)異常升高外，在周邊發(fā)生較近區(qū)域個別中小地震發(fā)生前也出現(xiàn)了明顯的異常升高現(xiàn)象.畢竟概率密度分布法是一種統(tǒng)計學手段，雖然能有效地識別出部分地震前的異常信號，但在某些λ2值異常升高時段周邊并沒有發(fā)生地震，如2008年9—11月λ2值的異常升高并無相應地震對應，這可能與不同地震的孕震機理不同相關(guān).

表1　云南高大井周邊地震信息表

識別地震前兆信息是一項較困難的工作，如何界定某個信息是正常還是異常更加不易.在積累探索過程中，多以未發(fā)生地震前的信息為背景值，當信號值與該背景值之間出現(xiàn)明顯差異可認為是異常信號.但是，由于每個臺站自身的井孔條件、信息捕捉能力以及所處的構(gòu)造位置不同，其地震前出現(xiàn)異常的形態(tài)和幅值也不盡相同，如圖6中各臺點高頻信息的概率密度分布λ2值，在汶川8.0級地震前其異常幅度也各不相同.如何在地震前有效地識別出異常信息，這需要一個不斷積累的過程，須對每個臺站、每個測項的歷史資料震例進行總結(jié)，進而歸納提煉出其異常指標，方可在日后的地震預報應用中加以應用.從圖9可知，本文中云南高大井水溫的高頻信息的λ2峰值高于0.02即可認為是異常信號，且其異常峰值多出現(xiàn)在震前45天內(nèi).進一步的統(tǒng)計結(jié)果顯示，該井水溫震前的高頻信息異常除了與震級、震中距有關(guān)外，還與地震的震源機制密切相關(guān)，如圖10所示.圖10b中方塊大小表示有地震對應的高頻信息λ2峰值時間與發(fā)震時間的時間差(dT，單位為天，地震序號與表1相對應)，方塊顏色表示時間差的大小，實線為時間差約為10天的分界；圖10a為各地震的震源機制解及其空間分布，藍色震源球表示異常峰值時間與發(fā)震時間的差值dT小于10天，紅色震源球表示dT大于10天.從圖10可以看出，表1所列21個地震中，異常峰值時間與發(fā)震時間的差值大于10天的地震幾乎全為走滑型地震(2008年汶川8.0級地震EQ1例外)，小于10天的地震多為非走滑型地震.這種震源機制解與異常持續(xù)時間之間的相關(guān)性，一定程度上反映了地震在孕震、發(fā)震過程中的區(qū)域構(gòu)造應力作用或大范圍的構(gòu)造動力學背景.

圖9　云南高大井水溫高頻信息概率密度時間序列Fig.9　Temporal sequence of probability density of high-frequency information from water temperature data of Gaoda well in Yunnan

圖10　云南高大井水溫高頻信息異常對應的地震信息(a)地震震源機制解及其空間分布；(b)地震震級、震中距及異常持續(xù)時間.Fig.10　High-frequency anomalies in water temperature of Gaoda well in Yunnan and corresponding earthquake data(a) Focal mechanism solutions of earthquakes and its spatial distribution； (b) Magnitude, epicenter distance and duration of abnormities.

因此，與晏銳等(2011)的研究類似，本文的研究是回溯性的個案研究，雖然概率密度分布法可有效地識別出數(shù)字化流體觀測資料中的高頻信息異常，但僅是一種數(shù)學處理方法，對高頻信息異常的形成機理以及與地震孕育過程相互關(guān)系的分析僅限于初步討論，有待于積累更多的震例分析和進一步的深入研究.

致謝本項工作得到中國地震局監(jiān)測預報司預報管理處的大力支持，中國地震局預測研究所邵志剛博士在研究過程中給予了悉心指導和幫助，審稿專家的修改建議使作者對該項研究的內(nèi)容有了更深層次認識，在此表示衷心的感謝.

References

Burlaga L F, Lazarus A J. 2000. Lognormal distributions and spectra of solar wind plasma fluctuations: Wind 1995—1998.JournalofGeophysicalResearch, 105(A2): 2357-2364.

Cao L L, Gao A T. 2010. Coseismic response characteristics of digital records of water level and water temperature in Gansu caused by WenchuanMS8.0 earthquake.ActaSeismologicaSinica(in Chinese), 32(3): 290-299.

Carbone V, Veltri P, Bruno R. 1996. Solar wind low-frequency magnetohydrodynamic turbulence: extended self-similarity and scaling laws.NonlinearProcessesinGeophysics, 3: 247-261.

Castaing B, Gagne Y, Hopfinger E J. 1990. Velocity probability density functions of high Reynolds number turbulence.PhysicsD, 46(2): 177-200.

Daubechies I. 1988. Orthonormal bases of compactly supported wavelets.CommunicationsonPureandAppliedMathematics, 41(7): 909-996.

Ghashghaie S, Breymann W, Peinke J, et al. 1996. Turbulent cascades in foreign exchange markets.Nature, 381(6585): 767-770.

Guo Z J, Qin B Y, Feng X C. 1974. Discussion on the change of ground-water level preceding a large earthquake from an earthquake source model.ActaGeophysicaSinica(in Chinese), 17(2): 99-105.

He A H, Zhao G, Liu C L, et al. 2012. The anomaly characteristics before Wenchuan earthquake and Yushu earthquake in Qinghai Yushu and Delingha geothermal observation wells.ChineseJ.Geophys. (in Chinese), 55(4): 1261-1268, doi: 10.6038/j.issn.0001-5733.2012.04.021.

Huang N E, Shen Z, Long S R, et al. 1998. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis.Proc.R.Soc.LandA, 454(1971): 903-995.Kiyono K, Struzik Z R, Aoyagi N, et al. 2005. Phase transition in a healthy human heart rate.PhysicalReviewLetters, 95(5): 058101.

Liu Q, Zhang J. 2011. Application of S transform in analysis of strain changes before and after Wenchuan earthquake.JournalofGeodesyandGeodynamics(in Chinese), 31(4): 6-9.

Liu Y W, Yang X H, Liu Y M, et al. 2005. The characters of groundwater responses induced by Sumatra Ms 8.7 earthquake. In: Survey and Forecast Department of CEA ed. Effect in Chinese Mainland Induced by Sumatra 8.7 Earthquake in Indonesia, 2004 (in Chinese). Beijing: Earthquake Press, 131-258.

Lockner D A. 1996. Brittle fracture as an analog to earthquakes: Can acoustic emission be used to develop a viable prediction strategy.J.Acous.Emission., 14(3-4): S88-S101.

Lockner D A, Byerlee J D, Kuksenko V, et al. 1991. Quasi-static fault growth and shear fracture energy in granite.Nature, 350(6313): 39-42.

Manshour P, Ghasemi F, Matsumoto T, et al. 2010. Anomalous fluctuations of vertical velocity of Earth and their possible implications for earthquakes.PhysicalReviewE, 82(3): 036105.Manshour P, Saberi S, Sahimi M, et al. 2009. Turbulencelike behavior of seismic time series.PhysicalReviewLetters, 102(1): 014101.

Qiu Z H, Tang L, Zhang B H, et al. 2012. Extracting anomaly of the Wenchuan earthquake from the dilatometer recording at NSH by means of wavelet-Overrun Rate Analysis.ChineseJ.Geophys. (in Chinese), 55(2): 538-546, doi: 10.6038/j.issn.0001-5733.2012.02.016.

Qiu Z H, Zhang B H, Chi S L, et al. 2011. Abnormal strain changes observed at Guza before the Wenchuan earthquake.ScienceChinaEarthSciences, 54(2): 233-240.

Sorriso-Valvo L, Carbone V, Giuliani P, et al. 2001. Intermittency in plasma turbulence.PlanetaryandSpaceScience, 49(12): 1193-1200.

Sun X L, Liu Y W. 2008. Characteristics and mechanism of co-seismic responses in the Tayuan well.EarthquakeResearchinChina(in Chinese), 24(2): 105-115.

Sun X L, Liu Y W. 2012. Changes in groundwater level and temperature induced by distant earthquakes.GeosciencesJournal, 16(3): 327-337.Sun X L, Liu Y W, Yan R. 2011. Application of empirical mode decomposition method to groundwater data.JournalofGeodesyandGeodynamics(in Chinese), 31(2): 80-83.

Wang Q, Zhang P Z, Freymueller J T, et al. 2001. Present-day crustal deformation in China constrained by Global Positioning System measurements.Science, 294(5542): 574-577.

Wu Y Q, Jiang Z S, Yang G H, et al. 2011. Comparison of GPS strain rate computing methods and their reliability.GeophysicalJournalInternational, 185(2): 703-717, doi: 10.1111/j.1365-246X.2011.04976.x.

Yan R, Jiang C S, Zhang L P. 2011. Study on critical slowing down phenomenon of radon concentrations in water before the WenchuanMS8.0 earthquake.ChineseJ.Geophys. (in Chinese), 54(7): 1817-1826.Yang C J, Feng Z S, Song D W, et al. 2005. Preliminary application of wavelet analysis method to extract the characters of tidal factor in well water level before earthquakes.NorthwesternSeismologicalJournal(in Chinese), 27(2): 163-167. Yang M B, Kang Y H, Zhang Q, et al. 2009. Tendencious fall of groundwater table in Beijing region and recognition of earthquake precursor information.ActaSeismologicaSinica(in Chinese), 31(3): 282-289.

Zhang P Z, Shen Z K, Wang M, et al. 2004. Continuous deformation of the Tibetan Plateau from global positioning system data.Geology, 32(9): 809-812.

Zhang P Z, Wen X Z, Xu X W, et al. 2009. Tectonic model of the great Wenchuan earthquake of May 12, 2008, Sichuan, China.ChineseSci.Bull. (in Chinese), 54(7): 944-953.

Zhang P Z, Xu X W, Wen X Z, et al. 2008. Slip rates and recurrence intervals of the Longmen Shan active fault zone, and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan, China.ChineseJ.Geophys. (in

Chinese), 51(4): 1066-1073.Zhang S X, Zhang Z G, Liu J M, et al. 2002. Application and analysis of digitized observation data of underground water level in Hebei Province.Earthquake(in Chinese), 22(4): 89-93.

附中文參考文獻

曹玲玲, 高安泰. 2010. 汶川MS8.0地震引起的甘肅數(shù)字化水位、水溫同震響應特征分析. 地震學報, 32(3): 290-299.

郭增建, 秦保燕, 馮學才. 1974. 從震源孕育模式討論大震前地下水的變化. 地球物理學報, 17(2): 99-105.

何案華, 趙剛, 劉成龍等. 2012. 青海玉樹與德令哈地熱觀測井在汶川與玉樹地震前的異常特征. 地球物理學報, 55(4): 1261-1268, doi: 10.6038/j.issn.0001-5733.2012.04.021.

劉琦, 張晶. 2011. S變換在汶川地震前后應變變化分析中的應用. 大地測量與地球動力學, 31(4): 6-9.

劉耀煒, 楊選輝, 劉永銘等. 2005. 地下流體對蘇門答臘8.7 級地震的響應特征.∥中國地震局監(jiān)測預報司編. 2004年印度尼西亞蘇門答臘8. 7 級大地震及其對中國大陸地區(qū)的影響. 北京: 地震出版社, 131-258.

邱澤華, 張寶紅, 池順良等. 2010. 汶川地震前姑咱臺觀測的異常應變變化. 中國科學: 地球科學, 40(8): 1031-1039.

邱澤華, 唐磊, 張寶紅等. 2012. 用小波-超限率分析提取寧陜臺汶川地震體應變異常. 地球物理學報, 55(2): 538-546, doi: 10.6038/j.issn.0001-5733.2012.02.016.

孫小龍, 劉耀煒. 2008. 塔院井水位和水溫的同震響應特征及其機理探討. 中國地震, 24(2): 105-115.

孫小龍, 劉耀煒, 晏銳. 2011. 經(jīng)驗模態(tài)分解法在地下水資料處理中的應用. 大地測量與地球動力學, 31(2): 80-83.

晏銳, 蔣長勝, 張浪平. 2011. 汶川8.0級地震前水氡濃度的臨界慢化現(xiàn)象研究. 地球物理學報, 54(7): 1817-1826.

楊叢杰, 馮志生, 宋德偉等. 2005. 小波分析方法在提取井水位潮汐因子震前變化特征的初步應用. 西北地震學報, 27(2): 163-167.

楊明波, 康躍虎, 張慶等. 2009. 北京地下水位趨勢下降動態(tài)及地震前兆信息識別. 地震學報, 31(3): 282-289.

張培震, 聞學澤, 徐錫偉等. 2009. 2008年汶川8.0級特大地震孕育和發(fā)生的多單元組合模式. 科學通報, 54(7): 944-953.

張培震, 徐錫偉, 聞學澤等. 2008. 2008年汶川8.0級地震發(fā)震斷裂的滑動速率、復發(fā)周期和構(gòu)造成因. 地球物理學報, 51(4): 1066-1073.

張素欣, 張子廣, 劉俊明等. 2002. 數(shù)字化水位觀測資料的應用研究. 地震, 22(4): 89-93.

(本文編輯張正峰)

Extracting high-frequency anomaly information from fluid observational data:a case study of the WenchuanMS8.0 earthquake of 2008

SUN Xiao-Long1,2, WANG Guang-Cai1, YAN Rui3*

1ChinaUniversityofGeosciences(Beijing),Beijing100083，China2InstituteofCrustalDynamics，Beijing100085,China3ChinaEarthquakeNetworksCenter,Beijing100045,China

AbstractWith improvement of underground fluid observation technology, the sampling frequency of data has been obviously improved, especially after the digital transformation. There is a great deal of tectonic information in high-frequency observational data. It is a key problem how to extract the useful high-frequency information from theses data. We introduced the PDF (probability density distribution) method, and analyzed the observed fluid data of water level and water temperature at 72 stations in the area of the north-south seismic zone before the Wenchuan MS8.0 earthquake on 12 May, 2008. The analysis results show that 16 water level and 14 water temperature data had the high-frequency anomalies before the MS8.0 earthquake, and the points with abnormal information were concentrated in the tectonic belts of southwest Yunnan. The high-frequency anomalies appeared from south to north over time. On the basis of these results we considered that the PDF method has certain reliability and applicability to extract the high-frequency anomaly information from observational data of fluid. It can provide a reference for information extraction of high-frequency digital observation.

KeywordsFluid observational data； Probability density distribution； High-frequency abnormal information； Wenchuan MS8.0 earthquake

基金項目地震行業(yè)科研專項經(jīng)費項目(201408019-03)和國家自然科學基金資助項目(41502239)共同資助.

作者簡介孫小龍,男,1981年生,副研究員，主要從事地震地下流體研究及地震預測研究.E-mail:xlsun04@163.com *通訊作者晏銳，男，1978年生，副研究員，主要從事地震地下流體研究與地震預測研究.E-mail:ray_2005@163.com

doi:10.6038/cjg20160512 中圖分類號P315

收稿日期2015-06-20，2016-01-12收修定稿

孫小龍， 王廣才， 晏銳. 2016. 利用概率密度分布提取流體觀測資料中的高頻異常信息——以2008年汶川8.0級地震為例.地球物理學報，59(5):1673-1684,doi:10.6038/cjg20160512.

Sun X L, Wang G C, Yan R. 2016. Extracting high-frequency anomaly information from fluid observational data: a case study of the WenchuanMS8.0 earthquake of 2008.ChineseJ.Geophys. (in Chinese),59(5):1673-1684,doi:10.6038/cjg20160512.