地下水位中潮汐与气压效应分析

刘 阳 何案华 赵 刚 张 帆 孙召华

1 海南省地震局,海口市美苑路49号，570203 2 中国地震局地壳应力研究所(地壳动力学重点实验室), 北京市安宁庄路1号，100085

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



地下水位中潮汐与气压效应分析

刘 阳1何案华2赵 刚2张 帆1孙召华3

1 海南省地震局,海口市美苑路49号，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.