高階統計量在形變固體潮汐異常信息提取中的應用①

李希亮，劉希強，張 玲，董曉娜，盧雙苓，王 強

（山東省地震局，山東 濟南 250014）

0 引言

固體潮汐研究是研究地震成因、地震前兆與地震預測的重要內容之一。近年來通過對全球地震資料的分析，已明確了固體潮汐對地震的發生有著不可忽視的影響。地球固體潮汐是在日、月、地球等天體運動過程中所產生的潮汐形變，通過研究日、月引潮力引起的傾斜、重力值的變化，有可能了解到地殼局部巖石彈性性質的變化。固體潮汐觀測到的潮汐畸變、階變、脈沖等異常變化是震前地殼應力－應變場變化的反映，實際上與地殼介質的應力－應變狀態是密不可分的。在“八五”、“九五”、“十五”科技研究期間，利用固體潮觀測資料已經對震前地殼形變的中短期異常信息特征有所了解，已取得了一些重要研究成果［1］。

固體潮汐觀測數據可以很好地反映地球固體潮汐的變化特征及地震孕育過程的中短期異常特征，而數字化形變資料則提供了該方面研究的可能性。由于采樣率的大幅度提高，豐富了潮汐觀測的震前變化信息，如地震前短周期的脈動、突跳等短臨信息及全球大地震的同震變形等［2］。地球固體潮是一種地球整體受日月引潮力作用所產生的近似粘彈性形變的現象，是一種非平穩過程［3］。對于這種級聯影響的非線性非平穩過程，很難從平穩過程的角度直接用固體潮振幅異常（相位異常）進行地震預測。同時大量的固體潮觀測實踐與研究結果表明，在具有較高靈敏度及精度的潮汐觀測中，地震前兆信息對潮汐因子的影響是很微弱的，因此傳統的潮汐因子計算方法受到很大的限制。在系統分析和研究了地震預測領域常用的分析方法之后，本文用泰安地震臺數字化固體潮汐觀測數據，對濮陽ML4．6地震前后異常時間、空間上的特征進行分析，應用高階統計量方法來獲取異常背景信息。有可能對未來發震的時間及地點給出預估計。

1 理論模型與計算方法

1．1 高階統計量理論［4］

對非高斯信號來說，二階統計量只是其中一種信息，它不包含相位信息，因此對非最小相位系統的辨識便顯得無能為力。隨著信號處理應用領域的不斷擴大，非高斯、非最小相位、非因果、非平穩信號的處理問題已經成為現代信號處理研究的熱點，但基于二階統計量（如相關函數和功率譜等）的信號分析方法不能從信號中提取非最小相位信息及由于非高斯和非線性帶來的信息［5］。

高階統計量提供了傳統的二階統計量沒有的大量豐富信息，能適應于復雜信號或噪聲環境。當分析信號為非高斯信號，或者在時間序列分析中對非高斯的AR，MA，ARMA模型進行辨識時，高階統計量可以有效地抑制高斯或非高斯的有色噪聲，并抽取不同于高斯信號的多種信號特征等。理論上，凡是使用功率譜或相關函數等二階統計量進行過分析和處理，而又未得到滿意結果的問題都可以應用高階統計量方法來進行。高階統計量方法在地球物理資料數字信號處理和屬性參數提取等方面有較強的理論性和實用性。

高階統計量主要包括高階矩、高階累積量、高階矩譜、高階累積量譜等四種。

在實際應用中，常用的是隨機過程的高階累積量。設 ｛x（n）｝為平穩隨機過程，x1＝x（n），x2＝x（n＋τ1），…，xk＝x（n＋τk－1），則其中k階累積量和k階矩為

在k≥3時，k階累積量和k階矩統稱為高階統計量。

k階矩譜定義為k階矩的k－1維Fourier變換，即

k階累積譜定義為k階累積量的k－1維Fourier變換，即

1．2 TOC－AIC模型與計算方法

AIC信息準則由日本學者赤池弘治于1973年提出，其進行震相識別的實質就是求解背景噪聲與信號最佳劃分點的過程，換句話說就是找出信號時間序列的不穩定點，此點與AIC曲線極小點相對應，即為震相的到時點。AIC的一般計算式為

AIC（n）＝－2log（模型的最大似然函數）＋2n其中n為模型的獨立參數的個數。

在震相識別中常見的AIC方法有：自回歸AIC方法、基于神經網絡的AIC方法、基于小波變換的AIC方法等。AR－AIC 算法（M．LEONARD，1999），即自回歸的AIC方法，首先假設可以把地震記錄分為兩個平穩過程，在每個部分分別建立自回歸（AR）模型。這兩個AR模型是完全不同的，即一個只包括背景噪聲（噪聲AR模型）而另一個只包括地震信號（信號AR模型）。當AIC取最小值時，AR模型取最優的階數，此時AR模型與地震記錄這個時間序列最為相符，也就是噪聲與地震信號的最優劃分點，這一點就是震相的到時點。AIC計算式為

其中為噪聲模型的方差為信號模型的方差；mF為噪聲AR模型的系數的個數；mB為信號R模型的系數的個數。此做法是在假設地震信號是平穩的，而實際是非平穩的，所以對于低信噪比且初動成隱伏狀的地震信號，估計的初至誤差較大。

Maeta（1985）提出了不同與AR－AIC的算法，他的做法是直接從地震記錄圖中計算AIC的值，而沒有利用AR系數，其計算式為

在選定的時間窗內對k逐點搜索，AIC的值反映x這個時間序列以k為分界的前后兩段方差的最小值，AIC極小值的點即為P波到時。

AIC方法在震相識別中的應用已經很成熟，AR－AIC算法和VAR－AIC算法是比較常用的兩種方法，在快速、有效的處理地震記錄中起到了很大的作用，推進了地震預警的發展。

震相到時檢測技術就是假設震相到時前后的地震記錄是兩個不同的穩態過程，開始估計的到時把時間序列分成兩個局部穩態的時間序列，以檢測震相的到時。但是信噪比高低和震級大小對AIC值的影響起決定性作用。同樣的在地震前兆數字化觀測記錄中的時間序列中，如果存在異常信號的介入，將出現兩個不同的穩態過程。兩個局部穩態的時間序列的劃分也很關鍵，用以檢測地震前兆異常信號。

考慮到地震前兆信號的復雜性，不平穩性，以及低信噪比等特點，本文提出了一種直接從前兆數字化觀測資料中計算AIC函數的方法。具體思路是用兩個時間段數據的三階累積量代替VAR－AIC方法中的方差，稱為TOC－AIC方法（third－order cumulant）。AIC函數表示為

模擬實驗證明，當隨機高斯噪聲中存在異常信號時，此方法可以抑制噪聲，精確地提出異常信號的初至信息。在不影響觀測分析結果精度的前提下，要對觀測數據進行預處理，盡量消除固體潮汐的周期性變化對計算的影響，排除干擾，突出異常信息。如果有新的地震前兆信號到達時，地震前兆數字化記錄中就會出現突變和漸變。采用TOC－AIC方法，從背景噪聲和地震前兆異常信號中找出信號時間序列中的不穩定點，也就是AIC函數的最小值，來識別地震前兆異常的初至。

2 震例研究

2．1 臺址及資料情況

泰安臺位于泰山南麓，地處魯中隆起，泰山山前第四紀活動斷裂北側。臺基為太古代花崗片麻巖，巖體完整致密均勻，測量信噪比高，可靠性良好，為全國一類形變基準臺。儀器洞室處坐標是東經117°07′22″，北緯36°12″35′。主洞進深約76m，最大覆蓋厚度約29m，室溫年變幅約0．06℃。布設有水平擺傾斜儀、垂直擺傾斜儀，還平行布設數字化SS－Y伸縮儀、DSQ水管傾斜儀各一套，共計11個測項，均為“九五”數字化前兆臺網改造之后新上的數字化儀器，已經積累了近10年的資料，資料可靠性高（圖1）。數字化傾斜儀的精度可以達到高于模擬記錄的水平，水平擺傾斜儀精度、穩定性都高于垂直擺傾斜儀［6］。這就是在同一洞室里面水平擺的潮汐異常最明顯的原因所在。由于臺站及其周邊地區中強地震較少，相關的形變潮汐短臨異常出現的很少，因此濮陽ML4．6地震潮汐短臨異常的識別對今后的地震研究工作有一定的借鑒意義。

圖1 泰安臺形變儀器布設示意圖Fig．1 Schematic diagram of deformation instruments at Tai＇an seismic station．

2．2 正常背景場固體潮汐TOC－AIC方法研究

在正常背景場下，固體潮汐觀測數據的TOCAIC方法研究結果如圖2所示。正常背景場下的處理結果比較平滑，反映不出異常信息的介入。

圖2 2006年1月泰安臺垂直擺NS向的原始資料和TOC－AIC值Fig．2 The recorded data and TOC － AIC values of NS vertical pendulum at Tai＇an station in January 2006．

2．3 異常背景場固體潮TOC－AIC方法研究

2006年4月9日在與山東莘縣相鄰的河南濮陽境內發生了ML4．6地震，震源深度12km，震中距泰安地震臺165km。在地震前泰安臺的潮汐形變觀測不同程度的出現了短臨異常。將TOC－AIC方法分析引入潮汐形變資料處理中，明顯地提取到了我們想要得到的地震前兆信息。圖3、4顯示結果對異常信息的反映程度雖有不同，但都能明顯的反映出異常信息的介入。其中垂直擺的資料分析結果較好的反映了臨震信號的介入（圖4（a）、（b）），水平擺的分析結果則較好的反映出了異常信號的介入和臨震信號的介入（圖4（c）、（d））。

圖3 濮陽地震前后傾斜儀原始記錄Fig．3 The recorded data of pendulum before and after the Puyang earthquake．

3 異常機理分析及其干擾排除

3．1 區域應力場背景特征

山東地區整體平均主壓應力方向為近EW向。圖5為1970年以來初動符號得到的山東及鄰區中小地震的震源機制解得到的主壓應力水平投影，主壓應力軸優勢方位為80°左右，主張應力軸優勢方位為350°左右。殷海濤等在山東GPS觀測網的基礎上，對山東地區地殼運動特征進行了分析，認為山東地區的主應變場呈NEE－SWW向擠壓，這是由于魯西塊體主要受來自青藏塊體NEE向的擠壓［7］。

濮陽ML4．6地震發生在聊考斷裂帶的西側。聊考斷裂帶為一被第四系掩蓋的隱伏斷裂帶，南起河南蘭考，北至聊城以北與齊廣斷裂交會。全長達270多公里。根據物探和鉆探資料，斷裂走向NNE20°～40°，傾向 NW，傾角35°～60°，為一東升西落的正斷層。新第三紀以來，仍表現出強烈活動性，也是山東地震源地之一。濮陽ML4．6地震的震源機制解顯示主破裂面走向為312．00°，傾角82°，滑動角13°，主壓應力方向為85°，主張應力方向為356°，與該區正常情況下主壓應力場方向85°±15°、主張應力場方向350°±15°比較一致，同時符合山東地區地殼運動特征。初步推測此次地震的發震構造為聊考斷裂帶內規模較小的NW或NNW之次級隱伏斷裂，NNE向聊考斷裂為控震構造。

3．2 異常機理分析

泰安臺震前出現短、臨前兆異常，其中水平擺兩項短期異常最為顯著可靠。結合濮陽ML4．6地震震源機制解和泰安臺與震中的位置可見，地震發生前泰安臺附近主要受來自NNW方向的壓應力和NEE方向的張應力作用，與水平擺出現的NE向傾斜異常相一致。這說明泰安臺水平擺出現的異常很可能是震中區應力增強導致的近源區應力變化的結果。隨著斷層上應力的增強—釋放—減弱，異常的表現也隨之出現相應的變化，直到最后恢復到原來正常水平［8］。泰安臺與濮陽震中之間有多條斷層縱橫交錯，雖然許多斷層之問沒有直接相連，但它們之間互相交錯，使得彼此又具有內在的聯系。尤其是NNE向的聊考斷裂帶與EW、NNW向的多條斷裂相交，其活動對與之相關的次級斷裂將起到控制作用。因而，泰安臺出現的傾斜異常應歸為“場兆”異常，是大區域應力場及伴生的多條斷層綜合效應的結果。

圖4 泰安臺不同擺不同方向的TOC－AIC值Fig．4 TOC － AIC values of vertical and horizontal pendulum in different directions at Tai＇an station．

3．3 異常干擾排除

形變觀測受多種因素影響，其中包含了氣象變化、人為作用、儀器故障及未知的諸多情況，使得干擾排除成為識別異常的一項必不可少的工作。但在目前的認識水平下，干擾形態與地震前兆形態極為相似，給干擾與前兆的判定帶來很大困難［9－10］。經過與臺站工作人員的共同落實，排除了氣壓、雷電、洞室進人、儀器故障及附加載荷等干擾因素對數字化形變高頻記錄的影響，提高了濮陽ML4．6地震泰安臺形變潮汐短臨異常的可靠性。

4 結論與討論

地震孕育過程中固體潮潮汐伴隨著地震孕育產生一定的變化，這種變化在具有較高靈敏度及穩定性的潮汐觀測過程中，采用較高分辨率的高階統計量分析方法能較好的識別。處理結果中的低點就是異常信號的介入點，這對地震異常信息的提取有指導意義，此方法能較好的提取潮汐形變短臨異常。

圖5 山東及鄰區單震震源機制解主壓應力軸水平投影Fig．5 The horizontal projection of principal compression stress axis from single earthquake in Shandong province and its adjacent areas．

泰安臺水平擺潮汐短臨異常形態明顯，震后完全恢復到正常年變形態。異常初期表現為地面NEE向傾斜，正好背向震中，與多數傾斜異常統計結果一致。研究表明，在同一洞室水平擺和垂直擺對地震異常的反應程度有差異，雖然它們同是記錄傾斜潮汐變化的儀器，但由于原理與結構的不一致，對地震異常的反映程度不一樣，高階統計量方法分析結果中的異常初至也不一致。高靈敏度的水平擺傾斜儀器撲捉地震前兆的能力要好一些，對地殼表面運動反應更敏感，更有利于獲取短臨異常。泰安臺距震中約165km，進一步驗證了6．5級以下的中強地震，突變點集中分布在震源附近200km范圍內的結論［11］。

本文只是應用某些前兆資料來說明高階統計量分析方法在提取固體潮潮汐異常的可行性和有效性，對于前兆機理和前兆異常的分析也是比較膚淺的。本文采用高階統計量方法對地震前兆信號進行分析與特征提取，研究了信號偏離高斯分布的程度，為地震前兆的識別提供了一條新的思路。期望給地震預測提供一些有用的參考。

［1］ 劉序儼，洪星．論潮汐波與地震波［J］．大地測量與地球動力學，2003，23（2）：73－76．

［2］ 張燕，吳云，劉永啟，等．小波分析在地殼形變資料處理中的應用［J］．地震學報，2004，26（增）：103－109．

［3］ 李瑞浩．重力學引論［M］．北京：科學出版社，1988．

［4］ 王書明，朱培民，李宏偉，等．地球物理學中的高階統計量方法［M］．北京：科學出版社，2006．

［5］ 李希亮，劉希強，董曉娜，等．高階統計量方法在地球物理學中的應用與展望［J］．西北地震學報，2010，32（2）：201－205．

［6］ 王梅，宋治平，李峰，等．由泰安地震臺形變資料看固體潮數字化觀測運行［J］．地震研究，2004，27：49－56．

［7］ 殷海濤，李 杰，張 玲，等．基于GPS觀測網的山東地區地殼運動特征分析［J］．西北地震學報，2008，20（3）：276－281．

［8］ 李杰，李希亮，盧雙苓，等．濮陽ML4．6地震前泰安臺形變異常特征分析［J］．大地測量與地球動力學，2007，27（4）：100－104．

［9］ 王梅，李峰，孔向陽，等．數字化形變觀測干擾識別［J］．大地測量與地球動力學，2004，24（1）：94－98．

［10］ 王梅，宋治平，李峰，等．形變數字化資料綜合分析［J］．大地測量與地球動力學，2003，23（4）：60－64．

［11］ 牛安福，張晶，高福旺，等．地殼持續加速變形與地震關系研究［J］．大地測量與地球動力學，2002，22（1）：29－33．