劉 陽(yáng) 何案華 趙 剛 張 帆 孫召華

1 海南省地震局,?？谑忻涝仿?9號(hào)，570203 2 中國(guó)地震局地殼應(yīng)力研究所(地殼動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室), 北京市安寧莊路1號(hào)，100085

3 河南省地震局, 鄭州市正光路10號(hào)，450016

?

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

劉 陽(yáng)1何案華2趙 剛2張 帆1孫召華3

1 海南省地震局,?？谑忻涝仿?9號(hào)，570203 2 中國(guó)地震局地殼應(yīng)力研究所(地殼動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室), 北京市安寧莊路1號(hào)，100085

3 河南省地震局, 鄭州市正光路10號(hào)，450016

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

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

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

(1)

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

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

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

1.1 氣壓效應(yīng)

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

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

1.2 潮汐效應(yīng)

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

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

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

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

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

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

Wold等[12-13]針對(duì)輸入量間的多重共線性問(wèn)題，提出偏最小二乘法(PLS)。利用偏最小二乘法對(duì)2008-04-28～05-13川03井水位觀測(cè)值數(shù)據(jù)進(jìn)行計(jì)算，如表1所示。結(jié)果顯示，潮汐分量取t0-t0、氣壓分量取p8-p8時(shí)，模型具有最小的AIC值。這意味著該時(shí)段井水位受潮汐引力作用與氣壓波動(dòng)作用既沒(méi)有滯后、也沒(méi)有延時(shí)；模型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表示過(guò)去第n小時(shí)到過(guò)去第m小時(shí)的固體潮汐理論值，pn-pm表示過(guò)去n小時(shí)到過(guò)去第m小時(shí)的氣壓觀測(cè)值。表2同。

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

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

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

4 認(rèn)識(shí)與討論

4.1 兩種模型的對(duì)比

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

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

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 語(yǔ)

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

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

[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]馬雄威. 線性回歸方程中多重共線性診斷方法及其實(shí)證分析[J].華中農(nóng)業(yè)大學(xué)學(xué)報(bào):社會(huì)科學(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

項(xiàng)目來(lái)源：中央級(jí)公益性科研院所基本科研業(yè)務(wù)費(fèi)專(zhuān)項(xiàng)(ZDJ2014-04)；國(guó)家自然科學(xué)基金(41104051)。

劉陽(yáng)，工程師，主要從事地下流體監(jiān)測(cè)研究，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.