楊光亮 周莉娟申重陽(yáng)韋進(jìn)

1)中國(guó)科學(xué)院研究生院地球科學(xué)學(xué)院地球動(dòng)力學(xué)實(shí)驗(yàn)室，北京1000392)中國(guó)地震局地震研究所，武漢4300713)地殼運(yùn)動(dòng)實(shí)驗(yàn)室，武漢430071

連續(xù)重力信號(hào)的非潮汐信息自動(dòng)提取*

楊光亮1，2，3)周莉娟2，3)申重陽(yáng)2，3)韋進(jìn)2，3)

1)中國(guó)科學(xué)院研究生院地球科學(xué)學(xué)院地球動(dòng)力學(xué)實(shí)驗(yàn)室，北京100039

2)中國(guó)地震局地震研究所，武漢430071
3)地殼運(yùn)動(dòng)實(shí)驗(yàn)室，武漢430071

采用基于EMD自動(dòng)提取算法及地震觸發(fā)判斷算法實(shí)現(xiàn)對(duì)臺(tái)站連續(xù)重力數(shù)據(jù)的潮汐、尖峰、臺(tái)階、儀器漂移和趨勢(shì)性變化的消除，達(dá)到連續(xù)重力非潮汐信息的自動(dòng)分離提取。研究實(shí)例表明，該算法可實(shí)現(xiàn)對(duì)重力臺(tái)站資料的自動(dòng)處理，適于對(duì)大量臺(tái)站連續(xù)重力觀測(cè)資料進(jìn)行批量化處理，快捷地提取非潮汐變化信息。

連續(xù)重力;潮汐;非潮汐;自動(dòng)提取;批量化處理

AbstractThe effect of tidal，peak，step，instrument drift and trend from the continuous tidal gravity data was removed by using the automatic extraction algorithm based on EMD and earthquake triggering algorithm，and achieved successfully the automatic extraction of tidal and non－tidal gravitational information.The examples show that the algorithm greatly reduces the manual work and it may be used for batch processing and quickly extracting the information of non－tidal changes of a large number of stations.

Key words:continuous gravity observation;tidal;non－tidal;automatic extraction;batch processing

1 引言

重力固體潮是由引潮力沿垂直方向的分量引起的地球重力場(chǎng)變化，其主要成分為長(zhǎng)周期潮、月潮、半月潮、日潮、半日潮、1/3日潮等等。連續(xù)重力臺(tái)站數(shù)據(jù)處理技術(shù)主要集中在潮汐部分的研究。PRETERNA和TSOFT可去除連續(xù)重力數(shù)據(jù)中的尖峰(突跳)、臺(tái)階和間斷等數(shù)據(jù)異常。各種調(diào)和分析方法雖在算法過(guò)程上略有差異，但功能類似，均能實(shí)現(xiàn)長(zhǎng)周期漂移與各種潮汐波信號(hào)的分離提取。

早期，連續(xù)重力觀測(cè)主要研究觀測(cè)數(shù)據(jù)的潮汐部分，如天體作用等，非潮汐部分通常作為干擾信息。隨著重力臺(tái)站連續(xù)重力數(shù)據(jù)由分鐘采樣到秒采樣過(guò)渡，重力固體潮汐不僅記錄了包括地球自身、日月及其他天體、地球外部大氣層等的引力作用，而且也記錄了地震發(fā)生前的巖石微破裂、錯(cuò)動(dòng)等地球內(nèi)部物理構(gòu)造變動(dòng)造成的密度分布變化及其他小區(qū)域地質(zhì)作用過(guò)程。因此，重力固體潮觀測(cè)不僅是研究地球內(nèi)部物質(zhì)分布、大氣負(fù)荷、地球液核動(dòng)力學(xué)、地殼運(yùn)動(dòng)、地球自轉(zhuǎn)等地球動(dòng)力學(xué)過(guò)程的利器，同時(shí)也可捕捉到反映地殼局部變化的高頻信息［1］。因而，連續(xù)重力觀測(cè)研究將在潮汐研究的基礎(chǔ)上，逐步向非潮汐研究延伸和擴(kuò)展。在充分利用現(xiàn)有觀測(cè)資料的前提下，進(jìn)一步發(fā)展新的連續(xù)重力非潮汐數(shù)據(jù)處理技術(shù)將是連續(xù)重力未來(lái)研究工作的重點(diǎn)。本文將對(duì)連續(xù)重力觀測(cè)的非潮汐信息快速、自動(dòng)提取技術(shù)進(jìn)行探討。

2 連續(xù)重力信號(hào)特征

設(shè)t時(shí)刻的重力觀測(cè)值y(t)為:

其中s(t)為重力潮汐，d(t)為漂移，p(t)為非潮汐信息，ε(t)為觀測(cè)誤差。

連續(xù)重力的潮汐分析主要是提取各種潮波，這些潮波的周期從幾小時(shí)到幾千小時(shí)間變化［2－6］，在頻譜上表現(xiàn)為甚低頻特征［3］，連續(xù)重力原始記錄中的尖峰、間斷、臺(tái)階、儀器漂移及趨勢(shì)項(xiàng)也基本表現(xiàn)為低頻特征，在頻譜上與潮汐信息存在交叉，因而很難通過(guò)頻率濾波的方式去除，因此，連續(xù)重力數(shù)據(jù)中的尖峰、間斷、臺(tái)階通常是以人工操作的方式去除，趨勢(shì)項(xiàng)和儀器漂移等通過(guò)多項(xiàng)式擬合去除，潮波主要通過(guò)理論計(jì)算來(lái)提取。該過(guò)程通常比較繁瑣，也是連續(xù)重力數(shù)據(jù)處理過(guò)程中需要人工完成的部分。

隨著連續(xù)重力研究向局部、動(dòng)態(tài)地震地質(zhì)變化研究的延伸，在連續(xù)重力處理中通常認(rèn)為是干擾的高頻信號(hào)的提取和分析也成為連續(xù)重力分析的一個(gè)重要方向?，F(xiàn)有研究理論認(rèn)為，較大地震發(fā)生前都伴隨著巖石的微小破裂，因而連續(xù)重力信號(hào)可能包含有臨震或同震信息，可以通過(guò)對(duì)大量臺(tái)站的連續(xù)重力記錄進(jìn)行分析，提取出具有共同特征的信息，捕捉這種臨震或同震的信號(hào)。本文基于EMD的自動(dòng)提取算法［7］對(duì)潮汐信號(hào)進(jìn)行分離提取，以期對(duì)大量臺(tái)站的連續(xù)重力進(jìn)行快速處理。

3 非潮汐信號(hào)自動(dòng)提取

模態(tài)分解(EMD)技術(shù)是一種信號(hào)處理技術(shù)，它可以根據(jù)信號(hào)本身的尺度特征對(duì)連續(xù)時(shí)間序列進(jìn)行模態(tài)提?。?－12］。

地震信號(hào)與其相應(yīng)的干擾背景在頻譜與能量上比較接近，而從連續(xù)重力主要潮波的周期上看，其潮汐與非潮汐信號(hào)在頻譜和能量上相差較大，理論上更容易實(shí)現(xiàn)分離提取。一般，24小時(shí)或48小時(shí)的連續(xù)重力包含了主要能量的潮汐波，因而連續(xù)重力數(shù)據(jù)的分析長(zhǎng)度通常是24或48小時(shí)。該時(shí)間長(zhǎng)度的秒采樣數(shù)據(jù)量較大，小波分析由于涉及大量傅里葉變換技術(shù)及要求數(shù)據(jù)長(zhǎng)度是2的n次方，而受到諸多限制。EMD算法不涉及譜域的計(jì)算，對(duì)數(shù)據(jù)長(zhǎng)度也不作限制，在計(jì)算效率上優(yōu)勢(shì)明顯。采用HHT算法對(duì)北京地球物理觀象臺(tái)24小時(shí)連續(xù)重力數(shù)據(jù)進(jìn)行進(jìn)行了分離提取(圖1)，共提取出自動(dòng)提取出12個(gè)模態(tài)，其中模態(tài)1是儀器漂移、長(zhǎng)周期潮汐及趨勢(shì)項(xiàng)信息，模態(tài)2是主要的潮汐信息，其余各模態(tài)是高頻的非潮汐信息。理論計(jì)算的潮汐與分離出的主要潮汐(圖1(e))得到的潮汐因子與通過(guò)潮汐分析軟件(Tsoft)得出的潮汐因子結(jié)果相近。通過(guò)該方法分離出了潮汐與非潮汐信號(hào)(圖1(f))。從分離結(jié)果看，分離出的非潮汐信號(hào)基本不再含有潮汐信息及長(zhǎng)周期變化，整個(gè)處理過(guò)程均可自動(dòng)完成，無(wú)需人工干預(yù)。因此，該方法適合對(duì)大量臺(tái)站連續(xù)重力數(shù)據(jù)進(jìn)行批量、自動(dòng)處理，從海量數(shù)據(jù)中初步提取非潮汐信息。

采用文獻(xiàn)［7］的算法，分別對(duì)成都臺(tái)連續(xù)重力進(jìn)行了自動(dòng)潮汐與非潮汐分離計(jì)算，較好地分離出潮汐與非潮汐信息(圖2)。

4 非潮汐信號(hào)尖峰消除

STA/LTA觸發(fā)算法是檢測(cè)瞬時(shí)突變信號(hào)的經(jīng)典算法，對(duì)文獻(xiàn)［7］中設(shè)計(jì)的地震觸發(fā)算法作適當(dāng)調(diào)整，可自動(dòng)判斷非潮汐時(shí)間序列中是否存在突變點(diǎn)，從而判斷尖峰位置。計(jì)算表明，上述比值設(shè)置為4時(shí)，可以基本識(shí)別出非潮汐信號(hào)的尖峰(圖3(a)、(b))。通過(guò)這種方式每次可去掉一部分尖峰，再對(duì)剩余時(shí)間序列重復(fù)上述步驟(1)和(2)，直到完全達(dá)到去除非潮汐時(shí)間序列中的尖峰為止。

圖3(c)和圖3(d)分別為迭代循環(huán)兩次后的結(jié)果，該非潮汐時(shí)間序列的峰值從600退化到100。通過(guò)人機(jī)交互運(yùn)算，尖峰可逐漸消除。

由此可見(jiàn)，使用該方法可使這一過(guò)程變得簡(jiǎn)便。

5 結(jié)論

采用基于EMD模態(tài)分解的自動(dòng)分離提取算法對(duì)連續(xù)重力信號(hào)進(jìn)行處理，能自動(dòng)或交互完成連續(xù)記錄中的潮汐(固體潮、海潮)、尖峰、臺(tái)階、漂移及趨勢(shì)性變化等信息，成功提取到重力非潮汐變化信息。通過(guò)實(shí)例處理結(jié)果表明，該方法適合對(duì)大量連續(xù)重力數(shù)據(jù)進(jìn)行批量處理。

但本工作只是從海量連續(xù)重力數(shù)據(jù)中初步自動(dòng)分離出了非潮汐信息，目前其作用主要體現(xiàn)在對(duì)海量連續(xù)重力數(shù)據(jù)的特征信息進(jìn)行初步篩選，更詳細(xì)、精確的分析還應(yīng)結(jié)合Tsoft、Baytap－G、gotic2等專業(yè)潮汐分析處理軟件。本文是實(shí)現(xiàn)自動(dòng)非潮汐信號(hào)提取的嘗試，后續(xù)還需對(duì)非潮汐信息進(jìn)行更詳細(xì)的分解，分析其各頻帶所包含的物理意義，特別是對(duì)多臺(tái)站地震前后的信息進(jìn)行提取分析，以期能捕捉到一些有意義的地震前兆觀測(cè)結(jié)果。

圖1 連續(xù)重力數(shù)據(jù)的HHT提取Fig.1HHT extraction of continuous gravity data

圖2 連續(xù)重力潮汐與非潮汐信息自動(dòng)提取Fig.2Automatic separation and extraction of tidal and nontidal information from continuous gravity data

圖3 連續(xù)重力非潮汐尖峰消除Fig.3Elimination of non－tidal peaks of gravity

1周江存，等.中國(guó)大陸精密重力潮汐改正模型.地球物理學(xué)報(bào)，2009，52(6):1 474－1 482.(Zhou Jiangcun，et al.Accurate correction models for tidal gravity in Chinese continent［J］.Chinese J.Geophys.，2009，52(6):1 474－1 482)

2周擎，等.固體潮的地震預(yù)測(cè)研究與地球動(dòng)力學(xué)研究之分析比較［J］.地球物理學(xué)進(jìn)展，2005，20(1):118－122.(Zhou Qing，et al.Comparing the earthquake forecast of the earth tide with the geodynamics of the earth tide［J］.Progress in Geophysics.，2005，20(1):118－122)

3蔣駿，等.固體潮潮汐因子的特征頻譜及其應(yīng)用［J］.地殼形變與地震，1994，(2):88－95.(Jiang Jun，et al.Characteristics spectrum of the Earth tidal factor and its application［J］.Crustal deformation and Earthquake，1994，(2):88－95)

4Doodson A T.The Harmonic development of the tide－generating potential［J］.Proceedings of the Royal Society of London.Series A，1921，100(704):305－329.

5Casotto and Biscani F.A fully analytical approach to the harmonic development of the tide－generating potential accounting for precession，nutation，and perturbations due to figure and planetary terms［J］.AAS Division on Dynamical Astronomy，2004，36(2):67.

6Cartwright D E.Tides:a scientific history［M］.Cambridge University Press，2001.

7楊光亮，等.基于HHT的地震信號(hào)自動(dòng)去噪算法［J］.大地測(cè)量與地球動(dòng)力學(xué)，2010，(3):39－42.(Yang Guangliang，et al.Automatic de－nosing algorithm of earthquake signal based on HHT decomposition［J］.Journal of Geodesy and Geodynamics，2010，(3):39－42)

8Norden E，et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and nonstationary time series analysis［J］.Proceedings of the Royal Society of London，1998，454:903－995.

9Norden E，et al.A new view of nonlinear water waves——the Hilbert spectrum［J］.Annual Review of Fluid Mechanics，1999，31:417－457.

10Flandrin P，Rilling G and Gon?alves P.Empirical mode decomposition as a filterbank［J］.IEEE Signal Proc Lett.，2004，11:112－114.

11于德介，程圣軍，楊宇.希爾伯特黃變換在齒輪故障診斷中的應(yīng)用［J］.機(jī)械工程學(xué)報(bào)，2005，41(6):102－107.(Yu Dejie，Cheng Shengjun and Yang Yu.Hilbert－Huang Transform in gear fault diagnosis［J］.Mechanical Engineering，2005，41(6):102－107.)

AUTOMATIC SEPARATION AND EXTRACTION OF NON-TIDAL INFORMATION FROM CONTINUOUS GRAVITY OBSERVATION

Yang Guangliang1，2，3)，Zhou Lijuan2，3)，Shen Chongyang2，3)and Wei Jin2，3)
1)Geodynamics Laboratory，College of Earth Sciences，Graduate University of CAS，Beijing100039
2)Institute of Seismology，China Earthquake Administration，Wuhan430071
3)Crustal Movement Laboratory，Wuhan 430071

P207

A

1671-5942(2011)03-0075-04

2011－02－13

中國(guó)地震局地震研究所重點(diǎn)基金(IS200916004)

楊光亮，男，博士，主要從事重力動(dòng)態(tài)變化與重力殼幔結(jié)構(gòu)反演研究工作.E－mail:vforyang@gmail.com