孫小龍， 王廣才， 晏銳

1 中國(guó)地質(zhì)大學(xué)(北京)， 北京　100083 2 中國(guó)地震局地殼應(yīng)力研究所， 北京　100085 3 中國(guó)地震臺(tái)網(wǎng)中心，北京　100045

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

——以2008年汶川8.0級(jí)地震為例

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

1 中國(guó)地質(zhì)大學(xué)(北京)， 北京100083 2 中國(guó)地震局地殼應(yīng)力研究所， 北京100085 3 中國(guó)地震臺(tái)網(wǎng)中心，北京100045

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

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

1引言

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

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

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

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

2概率密度分布法識(shí)別異常信息

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

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

對(duì)一組時(shí)序數(shù)據(jù)D(t)，以s為窗長(zhǎng)進(jìn)行三階多項(xiàng)式擬合，得到趨勢(shì)項(xiàng)擬合曲線為

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

(1)

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

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

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

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

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

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

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

(2)

式中，σ和λ分別為Z(t)和lnσ的標(biāo)準(zhǔn)差，σ0為σ的數(shù)學(xué)期望值.

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

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

圖2　隨機(jī)信號(hào)及其概率密度分布(a) 服從正態(tài)分布的隨機(jī)信號(hào)； (b) 服從對(duì)數(shù)正態(tài)分布的隨機(jī)信號(hào)； (c) 合成信號(hào)； (d) 正態(tài)分布的概率密度； (e) 對(duì)數(shù)正態(tài)分布的概率密度； (f) 合成信號(hào)的概率密度.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　隨機(jī)信號(hào)的疊加及其概率密度分布Fig.3　Stack of random signals and their probability density distributions

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

(3)

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

(4)

(5)

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

(6)

式中，P(z)表示觀測(cè)數(shù)據(jù)的概率密度分布，其標(biāo)準(zhǔn)差為σ，PCastaing為利用式(3)計(jì)算得到的概率密度分布，σCastaing為其相應(yīng)的標(biāo)準(zhǔn)差.

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

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

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

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

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

3汶川8.0級(jí)地震前流體異常提取

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

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

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

4結(jié)論與討論

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

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

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

圖7　汶川8.0級(jí)地震前流體高頻信息異常空間分布Fig.7　Spatial distributions of high-frequency fluid anomaly information before Wenchuan MS8.0 earthquake

圖8　三階多項(xiàng)式與小波分解法對(duì)比Fig.8　Comparison of three-order polynomial and wavelet decomposition methods

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

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

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

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

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

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

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

圖10　云南高大井水溫高頻信息異常對(duì)應(yīng)的地震信息(a)地震震源機(jī)制解及其空間分布；(b)地震震級(jí)、震中距及異常持續(xù)時(shí)間.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)的研究類似，本文的研究是回溯性的個(gè)案研究，雖然概率密度分布法可有效地識(shí)別出數(shù)字化流體觀測(cè)資料中的高頻信息異常，但僅是一種數(shù)學(xué)處理方法，對(duì)高頻信息異常的形成機(jī)理以及與地震孕育過程相互關(guān)系的分析僅限于初步討論，有待于積累更多的震例分析和進(jìn)一步的深入研究.

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

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.

附中文參考文獻(xiàn)

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

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

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

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

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

邱澤華, 張寶紅, 池順良等. 2010. 汶川地震前姑咱臺(tái)觀測(cè)的異常應(yīng)變變化. 中國(guó)科學(xué): 地球科學(xué), 40(8): 1031-1039.

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

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

孫小龍, 劉耀煒, 晏銳. 2011. 經(jīng)驗(yàn)?zāi)B(tài)分解法在地下水資料處理中的應(yīng)用. 大地測(cè)量與地球動(dòng)力學(xué), 31(2): 80-83.

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

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

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

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

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

張素欣, 張子廣, 劉俊明等. 2002. 數(shù)字化水位觀測(cè)資料的應(yīng)用研究. 地震, 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

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

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

doi:10.6038/cjg20160512 中圖分類號(hào)P315

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

孫小龍， 王廣才， 晏銳. 2016. 利用概率密度分布提取流體觀測(cè)資料中的高頻異常信息——以2008年汶川8.0級(jí)地震為例.地球物理學(xué)報(bào)，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.