地下水位中潮汐與氣壓效應(yīng)分析

劉 陽 何案華 趙 剛 張 帆 孫召華

1 海南省地震局,海口市美苑路49號，570203 2 中國地震局地殼應(yīng)力研究所(地殼動力學(xué)重點實驗室), 北京市安寧莊路1號，100085

3 河南省地震局, 鄭州市正光路10號，450016

?

地下水位中潮汐與氣壓效應(yīng)分析

劉 陽1何案華2趙 剛2張 帆1孫召華3

1 海南省地震局,海口市美苑路49號，570203 2 中國地震局地殼應(yīng)力研究所(地殼動力學(xué)重點實驗室), 北京市安寧莊路1號，100085

3 河南省地震局, 鄭州市正光路10號，450016

從井水位中潮汐與氣壓效應(yīng)的原理分析入手，著重解決輸入量間的多重共線性問題。采用偏最小二乘法回歸模型對水位數(shù)據(jù)中的潮汐與氣壓效應(yīng)進行計算發(fā)現(xiàn)，該方法得到的潮汐與氣壓響應(yīng)系數(shù)可正確反映潮汐引力與氣壓波動的物理過程；潮汐系數(shù)、氣壓系數(shù)、氣壓作用的滯后與延時常數(shù)等與含水層參數(shù)(儲水率、孔隙度、滲透率、體積模量等)以及井孔區(qū)域的地殼應(yīng)力應(yīng)變狀態(tài)密切相關(guān)。通過川03井連續(xù)60 d井水位(步長為9 d)的計算，可從響應(yīng)系數(shù)里識別明顯的地震同震響應(yīng)、震后調(diào)整過程和清晰的地震前兆異常。

潮汐效應(yīng)；氣壓效應(yīng)；偏最小二乘法；地震前兆

Matsumoto等[1-4]在地下水位觀測值中剔除氣壓效應(yīng)、潮汐效應(yīng)以及降雨效應(yīng)的主要思路是將地下水位觀測值按下式分解：

(1)

式中，yn為水位觀測值，xn為水位殘差，pn為氣壓效應(yīng)，Tn為潮汐效應(yīng)，Rn為降雨效應(yīng)；εn為測量噪聲，其均值為0、方差為σ2。

分析氣壓、潮汐與降雨各效應(yīng)的機理后[5]發(fā)現(xiàn)，響應(yīng)系數(shù)應(yīng)該是其響應(yīng)過程的體現(xiàn)，即通過系數(shù)的時間分布和數(shù)值大小可清晰反演潮汐引力、氣壓波動與降水荷載對井水位的作用過程。這些系數(shù)的時間分布和大小，跟含水層參數(shù)(儲水率、厚度、孔隙度、滲透率、體積模量等)、圍巖特性(骨架體積模量、覆蓋層厚度、巖性等)以及區(qū)域應(yīng)力應(yīng)變狀態(tài)息息相關(guān)[6-8]，但Matsumoto方法只存在純粹的數(shù)學(xué)模型，而得不到實際觀測中想研究的更多信息。本文就此進行改進。

1 氣壓效應(yīng)、潮汐效應(yīng)機理分析

1.1 氣壓效應(yīng)

井水位的氣壓效應(yīng)是由氣壓波動引起的井水位微動態(tài)[9]。圖1中，V1、V2分別為井孔內(nèi)和含水層內(nèi)水的體積，S1、S2分別為井孔及井孔-含水層交界面積。當氣壓(p0)同時作用在井區(qū)大地表面與觀測井水面上時，由于水體剛性很大，p0會直接傳遞到井-含水層界面(p2≈p0)；而作用在大地表面上的壓力因含水層以上的巖土變形而使其隨深度變小，作用到含水層頂板上的壓力明顯變小(p1

圖1 氣壓效應(yīng)原理圖Fig.1 Schematic diagram of atmospheric pressure effect mechanism of well level

1.2 潮汐效應(yīng)

Bredehoeft[8]指出，井孔的潮汐波動中蘊含了含水層力學(xué)和地殼應(yīng)變數(shù)據(jù)。井水位的地球固體潮效應(yīng)是指在太陽與月球等天體引力作用下含水層變形引起的井水位微動態(tài)。當?shù)厍蛳鄬ε蛎洉r含水層受張，孔隙壓力變小，引起井水回流到含水層，水頭下降；當?shù)厍蛳鄬嚎s時含水層受壓，孔隙壓力增大，引起含水層中地下水流入井中，水頭上升。太陽與月球的規(guī)律性運動，使地球與日月面距離、地球受到的膨脹與壓縮變形等都規(guī)律性地發(fā)生，井水位隨潮汐也規(guī)律性地變化。通過其動態(tài)過程可見，當潮汐理論值(體應(yīng)變值)增大時，水頭下降，水位觀測值增大；當潮汐理論值減小時，水頭上升，水位觀測值減小，即井水位觀測值與潮汐理論值之間呈正相關(guān)關(guān)系，因此潮汐響應(yīng)系數(shù)必須為正數(shù)。

通過上述分析可見，科學(xué)利用井水位觀測數(shù)據(jù)進行地震監(jiān)測與預(yù)測，不應(yīng)該只是簡單地剔除井水位中氣壓干擾與潮汐干擾，更重要的是從氣壓響應(yīng)系數(shù)與潮汐響應(yīng)系數(shù)里提取出隱含的含水層參數(shù)以及區(qū)域應(yīng)力應(yīng)變狀態(tài)。

2 狀態(tài)空間模型的局限性

對于氣壓效應(yīng)與潮汐效應(yīng)的剔除，利用狀態(tài)空間模型計算時，輸入量為當前到過去第n個氣壓觀測值、當前到過去第m個潮汐理論值。由于這些輸入量為當前值的滯后量，相互之間高度相關(guān)，即存在著嚴重的多重共線性[10]。

處理輸入量多重共線性問題的方法有[11]增加樣本容量、剔除法、嶺回歸法、偏最小二乘法、主成分分析法。通過比較發(fā)現(xiàn)，對于水位觀測值中的潮汐效應(yīng)與氣壓效應(yīng)分析，偏最小二乘法具有最小AIC值，這意味著偏最小二乘法可得到最接近實際過程的響應(yīng)模型。

3 偏最小二乘法計算潮汐響應(yīng)與氣壓響應(yīng)系數(shù)

Wold等[12-13]針對輸入量間的多重共線性問題，提出偏最小二乘法(PLS)。利用偏最小二乘法對2008-04-28～05-13川03井水位觀測值數(shù)據(jù)進行計算，如表1所示。結(jié)果顯示，潮汐分量取t0-t0、氣壓分量取p8-p8時，模型具有最小的AIC值。這意味著該時段井水位受潮汐引力作用與氣壓波動作用既沒有滯后、也沒有延時；模型SSE為0.195 546，決定系數(shù)為0.762 662 8，AIC值為-1 815.94。最優(yōu)模型如下：

yn=0.001 576 34tn+0.004 486 05pn-

3.147 360 89

(2)

注：tn-tm表示過去第n小時到過去第m小時的固體潮汐理論值，pn-pm表示過去n小時到過去第m小時的氣壓觀測值。表2同。

利用狀態(tài)空間模型對上述時間段的數(shù)據(jù)進行計算，并將Matsumoto原有模型中氣壓的起始量從i=0改為i=u,u∈(0,1,…,l)，即考慮氣壓作用的滯后因素。經(jīng)計算，修改后能得到AIC值最小的模型。

利用式(1)對2008-04-28～05-13川03井水位觀測值進行計算，其最優(yōu)模型如表2所示。結(jié)果顯示，當潮汐分量取t0-t15、氣壓分量取p7-p7時，有最小的AIC值，模型參數(shù)如表3所示。

模型SSE為0.160 008 91，AIC值為-1 861.32，決定系數(shù)為 0.797 304 78。單從AIC與SSE值來看，多元線性回歸(MLR)較PLS有一定優(yōu)勢，但從模型的F-統(tǒng)計[14]，t-統(tǒng)計[15]和p值[16]中可明顯看出其缺陷，如存在多個p值>0.5；模型計算出的響應(yīng)系數(shù)也存在與實際相違背的地方，潮汐響應(yīng)系數(shù)既有正數(shù)又有負數(shù)。導(dǎo)致這一結(jié)果的原因是多元線性回歸模型并未考慮輸入量間的共線性問題。

4 認識與討論

4.1 兩種模型的對比

圖2為分別采用狀態(tài)空間模型和PLS得到的井水位觀測值剔除潮汐及氣壓干擾的效果圖。從曲線的光滑程度來看，狀態(tài)空間模型的計算結(jié)果優(yōu)于PLS。通過潮汐與氣壓干擾的剔除，兩種方法都可以清晰地看到汶川地震的同震響應(yīng)以及震后調(diào)整過程[17]，同震響應(yīng)均為水位階升0.035 m左右，且震后水位呈現(xiàn)趨勢性上升過程。

表2 多元線性回歸分析AIC計算結(jié)果匯總

表3 多元線性回歸分析最佳模型參數(shù)

樣本數(shù): 384;自由度: 366;均方根誤差: 0.020 9;可決系數(shù): 0.806,決定系數(shù): 0.797。

圖2 狀態(tài)空間模型與PLS對比Fig.2 Comparison between the state space model and PLS

單純從濾波、干擾剔除的角度，狀態(tài)空間模型明顯優(yōu)于PLS。但由于狀態(tài)空間模型未考慮解釋變量間的多重共線性，導(dǎo)致其結(jié)果只有數(shù)學(xué)意義，而不是物理過程的還原，這對于利用井水位進行地震監(jiān)測與預(yù)測來說遠遠不夠。

4.2 潮汐響應(yīng)沒有滯后與延時、氣壓響應(yīng)既有滯后也有延時

Melchior[18]認為，井水位固體潮在相位上與起潮力W2非常接近(φe=0±3°)。Rhoads等[19]認為含水層的體膨脹是固體潮汐擾動的直接結(jié)果，海潮和氣壓潮通過改變覆蓋層重量造成次生含水層體膨脹。對2008-04-28～06-08(步長為9 d)川03井的數(shù)據(jù)計算表明，潮汐對含水層的作用既沒有滯后、也沒有延時，而是與固體潮汐同步；氣壓作用則普遍滯后6～8 h，持續(xù)作用時間1～2 h。

4.3 潮汐系數(shù)、氣壓系數(shù)與當?shù)貞?yīng)力應(yīng)變狀態(tài)關(guān)系密切

利用最小二乘法計算2008-04-01～06-08(步長為9 d)川03井潮汐與氣壓系數(shù)(圖3)?？梢钥闯?，沒有地震干擾時，潮汐響應(yīng)系數(shù)為0.001 7 m/10-9左右，氣壓響應(yīng)系數(shù)為0.006 m/hPa左右，且潮汐響應(yīng)系數(shù)與氣壓響應(yīng)系數(shù)具有同步變化的趨勢，但氣壓響應(yīng)系數(shù)曲線更光滑。地震發(fā)生時，潮汐與氣壓響應(yīng)系數(shù)都表現(xiàn)出明顯的同震過程。潮汐響應(yīng)系數(shù)震后變化頻率明顯高于氣壓響應(yīng)系數(shù)，可能是由于井-含水系統(tǒng)對氣壓響應(yīng)過程有“濾波”作用，不像引潮力直接作用于含水層上，氣壓則是以流體與覆蓋層為載體而體現(xiàn)在水位變化上的。

4.4 潮汐響應(yīng)系數(shù)中的地震前兆信息

從圖3(a)兩個圈里的潮汐響應(yīng)系數(shù)變化趨勢可以看出，汶川地震前潮汐響應(yīng)系數(shù)表現(xiàn)出明顯的異常。黃色圈里，在地震前21 d左右，潮汐響應(yīng)系數(shù)出現(xiàn)階降，體現(xiàn)了一種構(gòu)造活動。綠色圈內(nèi)，在地震前8 d，潮汐系數(shù)出現(xiàn)下降趨勢，但此時氣壓響應(yīng)系數(shù)仍基本保持穩(wěn)定，這可能是由于氣壓響應(yīng)系數(shù)經(jīng)過井-含水層的濾波，使其曲線趨于平緩導(dǎo)致。這種變化趨勢比較合理的解釋是：汶川地震前，由于斷層的蠕動導(dǎo)致川03井區(qū)域應(yīng)力處在張應(yīng)力不斷增強的過程[20]。

圖3 川03井潮汐與氣壓系數(shù)Fig.3 Tidal and atmospheric pressure coefficient of the Chuan 03# well

5 結(jié) 語

在解決輸入量多重共線性問題的基礎(chǔ)上，通過PLS模型計算水位中潮汐與氣壓效應(yīng)系數(shù)被認為是最接近實際作用過程的模型，其潮汐與氣壓響應(yīng)系數(shù)都表現(xiàn)出明顯的同震響應(yīng)過程。在曲線形態(tài)上，氣壓響應(yīng)系數(shù)較潮汐響應(yīng)系數(shù)要平滑，相當于經(jīng)過“低通的物理濾波器”。潮汐引力直接作用于含水層水體上，響應(yīng)過程沒有滯后與延時，而氣壓波動與之相反，體現(xiàn)出明顯的滯后與延時，且滯后與延時常數(shù)也體現(xiàn)出明顯的同震過程，潮汐響應(yīng)系數(shù)里隱含有明顯的震前異常。由于濾波作用，氣壓響應(yīng)系數(shù)里沒有體現(xiàn)異常信息。

致謝：本文水位觀測數(shù)據(jù)來源于中國地震臺網(wǎng)中心，氣象數(shù)據(jù)來源于美國國家氣象信息中心，觀測井資料由四川省地震局提供，在此一并表示感謝！

[1] Matsumoto N. Regression Analysis for Anomalous Changes of Groundwater Level due to Earthquakes[J]. Geophysical Research Letters,1992,19(12):1 193-1 196

[2] Matsumoto N, Takahashi M.Time Series Analysis for Detecting Changes of Groundwater Level due to Earthquakes at Hamaoka Observation Well, Shizuoka Prefecture, Central Japan[J]. J Seism Soc Japan,1993,45:407-415

[3] Matsumoto N, Kitagawa G, Roeloffs E A.Hydrological Response to Earthquakes in the Haibara Well, Central Japan-Ⅰ:Groundwater Level Changes Revealed Using State Space Decomposition of Atmospheric Pressure, Rainfall and Tidal Responses[J].Geophysical Journal International,2003,155(3):885-898

[4] Matsumoto N, Roeloffs E A.Hydrological Response to Earthquakes in the Haibara Well, Central Japan-Ⅱ. Possible Mechanism Inferred from Time-Varying Hydraulic Properties[J].Geophysical Journal International,2003,155(3):899-913

[5] Kamp G, Gale J E.Theory of Earth Tide and Barometric Effects in Porous Formations with Compressible Grains[J].Water Resources Research,1983,19(2):538-544

[6] Rexin E E, Oliver J, Prentiss D.Seismically-Induced Fluctuations of the Groundwater Level in the Nunn-Bush Well in Milwaukee[J].Bulletin of the Seismological Society of America,1962,52(1): 17-25

[7] Cooper H H, Bredehoeft J D, Papadopulos I S.The Response of Well-Aquifer Systems to Seismic Waves[J]. Journal of Geophysical Research,1965,70(16):3 915-3 926

[8] Bredehoeft J D.Response of Well-Aquifer Systems to Earth Tides[J].Journal of Geophysical Research, 1967, 72(12):3 075-3 087

[9] Jacob C E.On the Flow of Water in an Elastic Artesian Aquifer[J]. EOS, Transactions American Geophysical Union, 1940,21(2): 574-586

[11]馬雄威. 線性回歸方程中多重共線性診斷方法及其實證分析[J].華中農(nóng)業(yè)大學(xué)學(xué)報:社會科學(xué)版,2008,74 (2):78-81(Ma Xiongwei. Diagnosis and Empirical Analysis on Multicollinearity in Linear Regression Model[J]. Journal of Huazhong Agricultural University:Social Sciences Edition,2008, 74(2):78-81)

[12]Wold S, Albano C,Dunn W J. Pattern Recognition Finding and Using Regularities in Multivariate Data[J]. Food Research and Data Analysis,1983,3:147-188

[13]Wold S, Martens H, Wold H.The Multivariate Calibration Problem in Chemistry Solved by the PLS Method. In Matrix Pencils[M]. Berlin Heidelberg: Springer,1983

[14]Weir B S.EstimatingF-Statistics:A Historical View[J].Philosophy of Science,2012,79(5):637-643

[15]Mikusheva A. Second Order Expansion of thet-Statistic in AR (1) Models[J].Econometric Thery,2015,31(3): 426-448

[16]Nuzzo R.Statistical Errors[J]. Nature, 2014, 506(7 487):150-152

[17]Yan R, Woith H, Wang R J.Ground Water Level Changes Induced by the 2011 Tohoku Earthquake in China Mainland[J].Geophysical Journal International, 2014,199(1):533-548

[18]Melchior P.Diurnal Earth Tides and the Earth’s Liquid Core[J].Geophysical Journal International, 1966, 12(1): 15-21

[19]Rhoads G H, Robinson E S.Determination of Aquifer Parameters from Well Tides[J].Journal of Geophysical Research:Solid Earth,1979,84(B11):6 071-6 082

[20]Lippincott D K, Bredehoeft J D, Moyle W R.Recent Movement on the Garlock Fault as Suggested by Water Level Fluctuations in a Well in Fremont Valley[J].Journal of Geophysical Research:Solid Earth,1985,90(B2):1 911-1 924

Analysis of Tidal and Atmospheric Pressure Effects in the Water Level

LIUYang1HEAnhua2ZHAOGang2ZHANGFan1SUNZhaohua3

1 Earthquake Administration of Hainan Province, 49 Meiyuan Road, Haikou 570203, China 2 Key Laboratory of Crustal Dynamics, Institute of Crustal Dynamics, CEA, 1 Anningzhuang Road, Beijing 100085, China 3 Earthquake Administration of Henan Province,10 Zhengguang Road,Zhengzhou 450016, China

In this paper, we analyze the principle of tidal and atmospheric pressure effects of the water level, focusing on solving the problem of multicollinearity between the input variables. The Partial Least Squares Regression (PLS) model is used to eliminate the tidal and atmospheric pressure effects, and the results show that the response coefficients can correctly reflect the physical processes of tidal force and atmospheric pressure fluctuations. Both the magnitude of the tidal and atmospheric pressure coefficients, and the hysteresis and delay constants of atmospheric pressure, are closely related to the parameters of the aquifer (i.e., storage, porosity, permeability, bulk modulus, etc.),the crustal stress, and strain state around the borehole. Through continuous 60 days (Step size is 9 days) calculation of Chuan 03# well, the seismic and co-seismic responses, the adjustment processes after earthquakes, and even precursor anomalies are clearly identified in the tidal and atmospheric pressure response coefficients.

tidal effect; atmospheric effect; partial least square method; earthquake precursor

Special Fund for Basic Scientific Research of Central Public Research Institutes, No.ZDJ2014-04; National Natural Science Foundation of China,No.41104051.

HE Anhua，associate researcher，majors in method and theory of underground fluid, E-mail:dqs_hah@163.com.

2016-01-02

項目來源：中央級公益性科研院所基本科研業(yè)務(wù)費專項(ZDJ2014-04)；國家自然科學(xué)基金(41104051)。

劉陽，工程師，主要從事地下流體監(jiān)測研究，E-mail：liuyang_07@163.com。

何案華，副研究員，主要從事地下流體方法與理論研究，E-mail：dqs_hah@163.com。

10.14075/j.jgg.2016.12.018

1671-5942(2016)012-1112-05

P315

A

About the first author:LIU Yang, engineer, majors in research of underground fluid, E-mail:liuyang_07@163.com.