佘雅文 付廣裕 韋 進(jìn)

1 中國(guó)地震局地震預(yù)測(cè)重點(diǎn)實(shí)驗(yàn)室（中國(guó)地震局地震預(yù)測(cè)研究所），北京市復(fù)興路63號(hào)，100036

2 中國(guó)地震局地震研究所（地震大地測(cè)量實(shí)驗(yàn)室），武漢市洪山側(cè)路40號(hào)，430071

地下水變化與降雨是影響重力觀測(cè)的重要因素之一。近年來水文影響對(duì)連續(xù)重力觀測(cè)的研究已經(jīng)從經(jīng)驗(yàn)性轉(zhuǎn)變?yōu)槲锢砟Ｐ停?］。Tanaka使用兩臺(tái)gPhone獲取的重力殘差時(shí)間序列較好地反映了重力變化和降雨的響應(yīng)關(guān)系，并用Kazama建立的水文模型較好地模擬了重力變化［2－3］。利用降雨和靜水位的變化來確定長(zhǎng)期定點(diǎn)重力觀測(cè)的水文響應(yīng)模型，對(duì)去除干擾、獲取當(dāng)?shù)氐貧?nèi)部物質(zhì)遷移以及深部構(gòu)造活動(dòng)具有重要意義［4－5］。

本文利用十三陵地震臺(tái)的兩臺(tái)gPhone重力儀（序號(hào)109、118）2013－04～08 連續(xù)觀測(cè)數(shù)據(jù)進(jìn)行潮汐分析和提取重力殘差處理，從而對(duì)儀器的性能進(jìn)行判定，并在此基礎(chǔ)上分析gPhone連續(xù)重力觀測(cè)的水文響應(yīng)特征。

1 連續(xù)重力觀測(cè)

十三陵地震臺(tái)位于山區(qū)向平原過渡的斜坡地帶。兩臺(tái)gPhone重力儀并排放置在近似水平的平整瓷磚上，儀器旁放置數(shù)據(jù)采集和網(wǎng)絡(luò)通信設(shè)備，儀器和地面的接觸較為平整。利用兩臺(tái)重力儀在同一觀測(cè)條件下得到的數(shù)據(jù)來驗(yàn)證數(shù)據(jù)的可靠性，排除儀器自身原因造成的數(shù)據(jù)失真。由于斷電等原因，2013－06－19～07－04數(shù)據(jù)缺失。

2 數(shù)據(jù)處理

儀器原始數(shù)據(jù)為秒采樣，利用加漢寧窗的數(shù)字低通濾波器將秒采樣數(shù)據(jù)濾波為分鐘采樣，使用Tsoft軟件以480為窗口數(shù)、12cpd為采樣率將分鐘采樣降為小時(shí)采樣。采用Tsoft預(yù)處理軟件［6］對(duì)觀測(cè)曲線中的尖峰、地震和掉格等干擾信號(hào)進(jìn)行手工修正，對(duì)較短的記錄中斷用三次多項(xiàng)式進(jìn)行插值填補(bǔ)，保留長(zhǎng)時(shí)間中斷，并把秒采樣數(shù)據(jù)轉(zhuǎn)換成1min至1h采樣，對(duì)小時(shí)值數(shù)據(jù)進(jìn)行極移校正和水平校正。

利用BAYTAP－G 軟件進(jìn)行數(shù)據(jù)處理分析。BAYTAP－G 軟件可以將重力數(shù)據(jù)分離為隨機(jī)噪聲、固體潮汐、響應(yīng)數(shù)據(jù)項(xiàng)和漂移項(xiàng)4個(gè)部分，同時(shí)給出對(duì)數(shù)據(jù)潮汐因子的估計(jì)值［7］。使用輔助氣壓、傳感器溫度、水平校正值和相對(duì)重力數(shù)據(jù)一起作為輸入，其中，輔助氣壓、傳感器溫度和水平校正值為輔助數(shù)據(jù)。氣壓和傳感器溫度是重力觀測(cè)的重要影響因素，氣壓的影響大于傳感器溫度的影響。傳感器溫度雖然是自相關(guān)的，但其突變會(huì)引起觀測(cè)數(shù)據(jù)μGal級(jí)的變化，即每0.001 ℃可能會(huì)引起2.6μGal的重力變化［3］。選取水平校正作為輔助項(xiàng)是因?yàn)?，雖然儀器本身的預(yù)處理中包括水平校正，但是儀器的傾斜角度是隨時(shí)間不斷變化的，由儀器記錄的水平校正數(shù)據(jù)來看，影響在幾十μGal，且變化無規(guī)律。BAYTAP－G 在分離數(shù)據(jù)的過程中并不會(huì)因?yàn)檩o助數(shù)據(jù)在原始數(shù)據(jù)中不存在而產(chǎn)生錯(cuò)誤結(jié)果，因而忽略輔助數(shù)據(jù)的影響［7］。因?yàn)間Phone的漂移不是線性的［3］，所以選擇二階最小二乘法擬合觀測(cè)數(shù)據(jù)。對(duì)分離出的漂移信號(hào)作二階最小二乘法擬合，用漂移項(xiàng)數(shù)據(jù)減去擬合數(shù)據(jù)即可得到重力殘差［1］。

3 數(shù)據(jù)分析

3.1 兩臺(tái)gPhone重力儀觀測(cè)數(shù)據(jù)的對(duì)比分析

圖1給出了兩臺(tái)gPhone重力儀在2013－04－20～06－19觀測(cè)到的預(yù)處理數(shù)據(jù)、輔助觀測(cè)值和殘差時(shí)間序列，其中1（a）為gPhone重力儀觀測(cè)到的原始數(shù)據(jù)下采樣后的小時(shí)值數(shù)據(jù)，1（b）為儀器傳感器溫度，1（c）為gPhone重力儀水平校正值，1（d）為儀器圍壓值，1（e）為重力殘差值。水平校正值和傳感器溫度作為BAYTAP－G的輔助項(xiàng)得到的重力漂移再減去其二階最小二乘法擬合值，得到重力殘差時(shí)間序列（05－09 23：00～05－11 12：00之間gPhone109無測(cè)量數(shù)據(jù)）。從圖1（b）和圖1（c）可以看出，gPhone109和gPhone118雖然放在同一個(gè)觀測(cè)室內(nèi)，但兩臺(tái)儀器的傳感器溫度和水平校正項(xiàng)有較大的差異。前者在2013－04－20～06－20時(shí)段內(nèi)傳感器溫度僅在53.463 ℃～53.465 ℃變化，且大部分時(shí)間維持恒溫，而后者傳感器溫度在53.094 ℃～53.103 ℃之間呈階梯狀變化，可見118號(hào)儀器的傳感器溫度變化沒有109號(hào)儀器穩(wěn)定。109和118號(hào)儀器的水平校正項(xiàng)變化幅度分別為10μGal和35μGal左右，水平穩(wěn)定性上109 號(hào)儀器也強(qiáng)于118 號(hào)。使用BAYTAP－G 分離出的109號(hào)和118號(hào)儀器漂移項(xiàng)通過一階最小二乘法擬合后計(jì)算的漂移速率分別約為5.2μGal／d和5.9μGal／d，109號(hào)儀器的漂移速率比118號(hào)低，造成這種情況的原因可能是109號(hào)儀器的水平性能和傳感器溫度變化都比118號(hào)穩(wěn)定?？傮w而言，109號(hào)gPhone重力儀的性能優(yōu)于118號(hào)重力儀。

圖2為兩臺(tái)重力儀在2013－04－20～06－20的觀測(cè)數(shù)據(jù)處理得到的重力殘差變化曲線，其中圖2（a）為BAYTAP－G 處理兩臺(tái)重力儀數(shù)據(jù)得到漂移項(xiàng)減去其二階最小二乘法擬合得到的重力殘差，圖2（b）為兩臺(tái)儀器二階重力殘差的三次樣條插值。從圖2（b）可以看出，重力變化的趨勢(shì)總體是一致的，局部略有差異。由于兩臺(tái)儀器外部環(huán)境相同，所以差異來源于儀器本身。從圖2（a）中可以看到，109號(hào)重力儀比118 號(hào)重力儀的重力殘差變化要平緩，且基本看不到周期性變化，而118重力殘差中可看到明顯的周期變化成分。

圖1 gPhone109和118號(hào)重力儀2013－04－20～06－20預(yù)處理數(shù)據(jù)、輔助觀測(cè)值和殘差時(shí)間序列Fig.1 Time series of preprocessing data，auxiliary data and gravity residual of gPhone meters（serial number109and 118）from 2013－04－20～06－20

圖2 2013－04－20～06－20兩臺(tái)重力儀殘差時(shí)間序列比較Fig.2 Comparison of time series of the gravity residuals of two gPhone gravimeters from 2013－04－20～06－20

3.2 兩臺(tái)重力儀記錄數(shù)據(jù)的調(diào)和分析

選取109和118號(hào)gPhone重力儀2013－04－20～06－20的連續(xù)觀測(cè)記錄進(jìn)行處理和分析。表1是109和118號(hào)gPhone重力儀記錄的重力數(shù)據(jù)通過下采樣為小時(shí)值后，使用BAYTAP－G［8］以及ETERNA33［9］進(jìn)行潮汐分析的結(jié)果。

表1 BAYTAP－G、ETERNA33和DEHANT理論潮汐因子對(duì)比Tab.1 Comparison of tidal factors of BAYTAP－G，ETERNA33and DEHANT

從表1可見，兩臺(tái)儀器的數(shù)據(jù)使用不同處理軟件得到的潮汐振幅因子及其中誤差較為接近，證明了本文數(shù)據(jù)處理的正確性。為了提供對(duì)比參考，將BAYTAP－G 和ETERNA33的振幅因子與DEHANT 理論模型振幅因子［10］進(jìn)行比較?？梢钥闯觯珺AYTAP－G 處理得到的潮汐振幅因子整體較為接近理論值，其原因可能是BAYTAP－G可使用輔助數(shù)據(jù)抑制數(shù)據(jù)校正不充分在潮汐分析中的干擾。此外，109號(hào)儀器與118號(hào)儀器的潮汐振幅因子相比整體更接近于理論值。

綜上，109號(hào)gPhone儀器整體性能優(yōu)于118號(hào)儀器。由于這個(gè)原因，之后有關(guān)水文活動(dòng)造成重力變化的討論都將圍繞109號(hào)gPhone重力儀的觀測(cè)數(shù)據(jù)進(jìn)行。

3.3 降雨造成的重力變化分析

圖3（a）和圖3（b）分別給出了2013－05－25～06－19和2013－07－05～08－04的109號(hào)gPhone重力儀重力殘差變化和降雨的時(shí)間序列圖。從圖中可以看到，大部分重力變化與降雨存在較好的對(duì)應(yīng)關(guān)系。大部分情況下，它們之間的對(duì)應(yīng)關(guān)系滿足降雨時(shí)重力變小、降雨后重力逐漸變大的規(guī)律。但是也有少數(shù)情況不符合這種模式，這可能是由于其他未知因素的重力變化引起的，其原因尚待進(jìn)一步研究。

為解釋這種重力隨降雨變化的特征，建立簡(jiǎn)單的物理模型（圖4（a））。圖4（b）中箭頭所指位置為109和118號(hào)gPhone重力儀的觀測(cè)室，觀測(cè)室西側(cè)緊鄰海拔為161m 的龍山，觀測(cè)室的海拔為101m。圖4（a）中，降雨開始時(shí)水層覆蓋在地表，由于龍山海拔高于觀測(cè)室，在龍山上的雨水對(duì)重力儀產(chǎn)生引力傾斜向上，與重力儀處于同一高度的雨水對(duì)重力儀的引力相互抵消，這時(shí)重力值變小。降雨結(jié)束后，由于山體的坡度，雨水更多地流向地勢(shì)低的地方，滲透到低于重力儀海拔的地下，這時(shí)重力值增大。Naujoks使用超導(dǎo)重力儀也觀測(cè)到了類似現(xiàn)象［4］。

圖3 gPhone109二階最小二乘法擬合殘差與降雨對(duì)應(yīng)的時(shí)間序列Fig.3 Correlation of time series of 2nd－order curve fitting residual gravity of gPhone meter（serial number 109）and rainfall

圖4 重力隨降雨變化的物理模型示意圖和儀器位置與周邊環(huán)境Fig.4 Response model of gravity to rainfall and the position and surrounding environment of gPhone meter

3.4 重力變化與地下水靜水位變化的對(duì)應(yīng)關(guān)系

圖5是119號(hào)gPhone重力儀觀測(cè)結(jié)果的二階最小二乘法擬合殘差和沙河地震臺(tái)靜水位關(guān)系圖。沙河地震臺(tái)與十三陵地震臺(tái)的距離約為15 km。由圖5可知，連續(xù)重力觀測(cè)殘差曲線與降雨量之間存在較好的相關(guān)性，水位下降時(shí)重力減小，水位上升時(shí)重力增大?？梢娫趯?duì)連續(xù)重力觀測(cè)進(jìn)行分析時(shí)，地下水的影響不可忽視，因此建議在全國(guó)的連續(xù)重力觀測(cè)網(wǎng)中配備地下水位觀測(cè)。

4月左右的沙河靜水位下降在重力觀測(cè)信號(hào)中無顯示，可能是由于gPhone重力儀在03－24才啟動(dòng)觀測(cè)，尚未完全穩(wěn)定所致。由圖5 可知，gPhone重力觀測(cè)大約需要1～2個(gè)月時(shí)間來穩(wěn)定信號(hào)，之后才能進(jìn)行較好的分析。109號(hào)gPhone重力儀觀測(cè)到的重力變化和沙河地震臺(tái)的靜水位存在4d左右的延遲對(duì)應(yīng)關(guān)系。延遲的原因可能是地下水流動(dòng)造成的，由于地下結(jié)構(gòu)比較復(fù)雜，地下水的活動(dòng)受很多地質(zhì)環(huán)境的影響，如地下水流動(dòng)、斷層的阻斷等，都會(huì)造成重力變化對(duì)水文活動(dòng)的延遲［4，11］。

圖5 gPhone二階最小二乘法擬合殘差與沙河地震臺(tái)靜水位Fig.5 Time series of 2nd－order curve fitting residual gravity and static water level of Shahe seismic station

4 結(jié) 語

本文對(duì)十三陵地震臺(tái)的兩臺(tái)gPhone重力儀（109和118號(hào)）在2013－04～08的連續(xù)觀測(cè)數(shù)據(jù)進(jìn)行潮汐分析和提取重力殘差處理，對(duì)儀器的性能進(jìn)行判定，并在此基礎(chǔ)上分析gPhone重力觀測(cè)的水文響應(yīng)。結(jié)果顯示：1）影響gPhone重力儀觀測(cè)精度的儀器自身因素主要來自傳感器溫度和儀器的水平穩(wěn)定性。2）利用BAYTAP－G 軟件和ETERNA33軟件獲得的潮汐分析結(jié)果大體一致，證實(shí)了兩套處理軟件的一致性。相對(duì)而言，BAYTAP－G 得到的主要振幅因子更接近DEHANT 給出的理論值。3）依據(jù)109號(hào)gPhone重力儀觀測(cè)結(jié)果給出的二階重力殘差變化和十三陵地震臺(tái)的降雨有μGal量級(jí)的對(duì)應(yīng)關(guān)系，并且降雨開始時(shí)重力值減小，降雨結(jié)束后重力值增大。該現(xiàn)象可根據(jù)觀測(cè)站周圍特定的地形條件定性解釋。4）109號(hào)gPhone重力儀獲取的二階重力殘差變化和距離十三陵地震臺(tái)15km 的沙河地震臺(tái)靜水位變化存在4d左右延遲的正相關(guān)對(duì)應(yīng)關(guān)系，該對(duì)應(yīng)關(guān)系顯示了地下水位觀測(cè)對(duì)重力觀測(cè)解釋的必要性，而4d左右的延遲則是由于重力觀測(cè)與地下水位觀測(cè)點(diǎn)存在15km 的差異所致。重力變化與地下水位變化之間的正相關(guān)對(duì)應(yīng)關(guān)系表明，在全國(guó)性連續(xù)重力觀測(cè)網(wǎng)中，有必要進(jìn)一步配備地下水位觀測(cè)，以提高連續(xù)重力觀測(cè)的數(shù)據(jù)處理精度。

［1］Tanaka T，Miyajima R，Asai H，et al.Hydrological Gravity Response Detection Using agPhone Below－and Aboveground［J］.Earth，Planets and Space，2013，65（2）：59－66

［2］Kazama T，Okubo S.Hydrological Modeling of Groundwater Disturbances to Observed Gravity：Theory and Application to Asama Volcano，Central Japan［J］.Journal of Geophysical Research：Solid Earth，2009，114（B8）

［3］Tanaka，T.Introduction／Adjustment of gPhone and Its Present Condition of the Data［J］.Rep Tono Res Inst Earthq Sci，2010，25：75－81

［4］Naujoks M，Kroner C，Weise A，et al.Evaluating Local Hydrological Modelling by Temporal Gravity Observations and a Gravimetric Three－Dimensional Model［J］.Geophysical Journal International，2010，182（1）：233－249

［5］韋進(jìn)，李輝，劉子維，等.武漢九峰地震臺(tái)超導(dǎo)重力儀觀測(cè)分析研究［J］.地球物理學(xué)報(bào)，2012（6）：1 894－1 902（Wei Jin，Li Hui，Liu Ziwei，et al.Observation of Superconducting Gravimeter at Jiufeng Seismic Station［J］.Chinese J Geophys，2012（6）：1 894－1 902）

［6］Camp M V，Vauterin P.Tsoft：Graphical and Interactive Software for the Analysis of Time Series and Earth Tides［J］.Computers ＆Geosciences，2005，31（5）：631－640

［7］Tamura Y，Sato T，Ooe M，et al.A Procedure for Tidal Analysis with a Bayesian Information Criterion［J］.Geophysical Journal International，1991，104（3）：507－516

［8］Akaike H.Likelihood and the Bayes Procedure［J］.Trabajos de Estadística Y de Investigación Operativa，1980，31（1）：143－166

［9］Wenzel H G.The Nanogal Software：Earth Tide Data Processing Package ETERNA 3.30［J］.Bull Inf Marées Terrestres，1996，124：9 425－9 439

［10］Dehant V.Tidal Parameters for an Inelastic Earth［J］.Physics of the Earth and Planetary Interiors，1987，49（1），97－116

［11］Pool D R，Schmidt W.Measurement of Ground－Water Storage Change and Specific Yield Using the Temporal－Gravity Method near Rillito Creek，Tucson，Arizona［R］.Tucson：US Geological Survey，1997