地下水位中潮汐與氣壓效應分析

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

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

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

?

地下水位中潮汐與氣壓效應分析

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

1 海南省地震局,?？谑忻涝仿?9號，570203 2 中國地震局地殼應力研究所(地殼動力學重點實驗室), 北京市安寧莊路1號，100085

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

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

潮汐效應；氣壓效應；偏最小二乘法；地震前兆

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

(1)

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

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

1 氣壓效應、潮汐效應機理分析

1.1 氣壓效應

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

圖1 氣壓效應原理圖Fig.1 Schematic diagram of atmospheric pressure effect mechanism of well level

1.2 潮汐效應

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

通過上述分析可見，科學利用井水位觀測數據進行地震監測與預測，不應該只是簡單地剔除井水位中氣壓干擾與潮汐干擾，更重要的是從氣壓響應系數與潮汐響應系數里提取出隱含的含水層參數以及區域應力應變狀態。

2 狀態空間模型的局限性

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

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

3 偏最小二乘法計算潮汐響應與氣壓響應系數

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

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

3.147 360 89

(2)

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

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

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

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

4 認識與討論

4.1 兩種模型的對比

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

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

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

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

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

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

4.2 潮汐響應沒有滯后與延時、氣壓響應既有滯后也有延時

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

4.3 潮汐系數、氣壓系數與當地應力應變狀態關系密切

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

4.4 潮汐響應系數中的地震前兆信息

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

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

5 結 語

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

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

[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].華中農業大學學報:社會科學版,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

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

劉陽，工程師，主要從事地下流體監測研究，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.