引言

对承压含水层井孔而言，井水位对体应变潮汐的响应是井水位观测中普遍存在的一种现象．井水位潮汐中隐含有丰富的井－含水层系统信息，深入分析井水位潮汐对推算含水层的各种参数，衡量井水位对应力应变的响应量，了解含水层动态，解释井水位微动态，认识井水位潮汐的正常动态，评价井水位观测质量以及捕捉地震前兆异常等，都具有重要的科学意义和现实必要性．自１８７９年捷克斯洛伐克的一涌水煤矿里发现水位潮汐波动现象后，关于潮汐现象的报道不断增多，但直到２０世纪６０年代，关于井水位潮汐现象的深入研究才真正开始．Ｂｒｅｄｅｈｏｅｆｔ（１９６７）以体应变固体潮理论为基础，结合水文地质学的考虑，从理论上分析了井－含水层系统对体应变潮汐的响应，并提出一种利用井水位潮汐信息来确定含水层的孔隙度、骨架体变模量和储水系数的方法．之后，不同研究者对利用潮汐信息确定含水层参数的方法进行了发展（Ｃａｒｒ，ＶａｎＤｅｒＫａｍｐ，１９６９；Ｈｓｉｅｈｅｔａｌ，１９８７；Ｒｈｏａｄｓ，Ｒｏｂｉｎｓｏｎ，１９７９）．其中，Ｈｓｉｅｈ等推导了井水位对含水层孔隙压力响应的振幅比和相位差数学表达式，并给出了利用井水位潮汐相位和振幅比来确定井水位渗透系数和储水系数的方法．张昭栋等（２００２）在总结了水位对体应变潮汐、大气潮汐与荷载潮异同的基础上，给出了井水位潮汐的统一数学物理方程和解析解表达式，该解析解与Ｈｓｉｅｈ等（１９８７）给出的解一致．Ｅｌｋｈｏｕｒｙ等（２００６）利用Ｈｓｉｅｈ等（１９８７）提供的方法，分析了加利福利亚地区两口井水位观测资料的潮汐振幅比和相位差，认为大震发生后的潮汐振幅比和相位变化是由于地震波增加含水层渗透系数所致．同时，井水位是反演应力－应变变化的重要手段（Ｂｏｄｖａｒｓｓｏｎ，１９７０；Ｂｒｅｄｅｈｏｅｆｔｅｔａｌ，１９８６；Ｗａｋｉｔａ，１９７５；王六桥，李善因，２００１；张昭栋等，１９９４）．Ｂｏｄｖａｒｓｓｏｎ（１９７０）推导了应变引起的井底压力波动与井孔水位波动之间的表达式，认为井－含水层系统在有利的条件下可以作为非常灵敏的应变计使用；张昭栋等（１９９４）分析了井水位潮汐系数与含水层应力之间的定量关系，并以日本秋田７．７级地震引起的井水位变化为例，估计了中国大陆应力场的调整．此外，井水位潮汐是评价水位观测质量的一个重要指标，对井水位潮汐振幅和相位的动态跟踪还可能捕捉到地震前兆信息（张昭栋，高玉斌，１９９４），但对井水位潮汐振幅比和相位变化的可能原因、影响因素、正常动态及远距离大震对含水层影响的讨论较少，而且在利用井水位潮汐评价井水位观测质量时，对计算误差、影响因素及正常动态的考虑略显不足．本文以井－含水层系统对体应变潮汐的响应理论为基础，利用Ｂａｙｔａｐ－Ｇ提供的潮汐分析方法（Ｉｓｈｉｇｕｒｏ，Ｔａｍｕｒａ，１９８５），计算了云南曲靖井水位５个潮汐波的振幅和相位值，深入分析了２次抽水和几次中远场大震后潮汐振幅和相位的变化特征，讨论了变化的可能原因，证实了地震波能改变含水层渗透系数并导致孔隙压力重新分布的理论，为解释远距离地震触发和水文地质响应现象提供了实例．同时，分析井水位潮汐的动态变化，为分析含水层参数的正常动态特征、井水位微动态和井水位观测资料的动态评价等提供了参考依据．

１观测井孔基本情况与水位变化基本特征

云南曲靖井位于云南省曲靖市区的沉积盆地内，井孔口标高１８６４ｍ．该井地处南北向的曲靖断裂中部，位于南北向的小江深大断裂以东、北东向的寻甸—来宾断裂以南、北东向师宗断裂以北，距小江断裂约５０ｋｍ．水系属南盘江上游．观测井深２５９．２２ｍ，共下套管２５１ｍ，以下为裸孔．其中，变径无缝钢管２２１．６ｍ，打眼花管２９．４ｍ，过水段直径为１０８ｍｍ（图１ｂ）．含水层顶板埋深２２．３５ｍ，含水层为泥盆系灰岩岩溶裂隙承压水，水化学类型为Ｃａ２＋－ＣＯ３２－．成井时，水位埋深为１８．００ｍ，经钻探及抽水试验检查，止水效果良好．自１９９８年８月７日起正式开始观测，观测项目有静水位和水温．其中，１９９８年８月７日—１９９８年１２月１７日采用ＳＷ－４０型水位仪观测，１９９８年１２月１８日—２００２年１月１５日采用ＳＷ４０－１型水位仪观测，２００２年１月１５日至今则采用ＬＮ－３型数字水位仪进行水位观测，始测水位实际埋深为２２．７８ｍ，探头放置于距井孔２４．７８ｍ（即水下２ｍ）处．该井水位自观测以来，资料连续性较好，观测资料除２００１和２００３年受抽水影响外，其它时段观测环境无明显干扰．该井地下水补给来源于大气降水，井水位观测资料总体上呈现出雨季上升、旱季下降的年变化特征．图２为该井水位与降雨量的观测曲线．由图２可见，２００４年后，观测资料变化相对稳定．２０１０年井水位曾出现过大幅度下降，累积下降幅度达３．５ｍ．该变化是否与云南地区２０１０年上半年的大范围干旱有关尚需进一步研究．该井水位微动态特征明显，能记录到清晰的固体潮和同震信息，固体潮峰谷日变幅在１０—８０ｍｍ之间，为该井潮汐动态的研究提供了较好的基础资料．而且，该井自观测以来其附近１００ｋｍ范围内无５级以上地震发生，２００ｋｍ范围内无６级以上地震发生，这对研究远距离大震对该井－含水层系统的影响、水位潮汐振幅和相位正常动态及影响因素提供了有利条件．

２资料的预处理与分析方法

２．１井水位潮汐分析方法井水位潮汐分析包括实际观测潮汐与理论潮汐的对比．潮汐振幅比和相位差是分析井水位潮汐的有效参数，潮汐振幅比是指实际观测潮汐振幅与理论潮汐振幅之比，相位差是指实际潮汐相位与理论潮汐相位之差．根据Ｈｓｉｅｈ等（１９８７）和张昭栋等（２００２）给出的理想井－承压含水层系统潮汐方程解析解，对于饱和连续的各项同性均匀孔隙弹性含水层介质，井水位潮汐振幅比和相位变化体现了含水层储水系数和导水系数等参数的变化特征．振幅比δ和相位差Δφ表达式分别为式中，ｈ０和Ｈ０分别是频系数为ω的井水位和含水层压力水头的潮汐分波振幅值；Ｅ≈１－ωｒ２ｃＫｅｉ（α）／２Ｔ，Ｆ≈ωｒ２ｃＫｅｒ（α）／２Ｔ，α＝ｒｗ（ωＳ／槡Ｔ），Ｋｅｒ和Ｋｅｉ分别为开尔文函数的实部和虚部，ｒｃ和ｒｗ分别为井孔套管半径和含水层处井孔半径，Ｓ和Ｔ分别为含水层的储水系数和导水系数．潮汐振幅比和相位差的计算主要包括实际潮汐和理论潮汐的估计．目前理论潮汐的计算方法已经较为成熟，而且有不同的软件包可以直接使用，计算中只需要提供观测站的经度、纬度、高程、重力加速度等信息即可（Ａｇｎｅｗ，１９９７；Ｄｏａｎｅｔａｌ，２００６）．而实际观测资料不同频率成分的潮汐振幅和相位的估计，通常有时间域和频率域两种处理方法．时间域方法主要是使用潮汐调和分析，即将潮汐信号分解为不同频率谐波信号的和；频率域方法主要是对潮汐观测资料进行频谱分析，确定最靠近潮汐频率的潮汐信号．频率域分析方法除了要考虑谱遗漏的问题外，计算精度与观测资料的采样率有很大关系．由于离散傅里叶变化的频率要尽可能接近理论潮汐的频率，故观测资料的采样率必须要高（Ｅｌｋｈｏｕｒｙｅｔａｌ，２００６）．潮汐分析的程序较多，目前主要的潮汐分析程序有：Ｂａｙｔａｐ－Ｇ、ＥＴＥＲＮＡ、ＰＩＡＳＤ、ＳＰＯＴＬ、ＥＲＴＩＤ、ＴＳＯＦＴ等．这些程序很难说那个程序好，各程序有各自的侧重点和各自的优缺点．本文中我们使用了国际上通用的Ｂａｙｔａｐ－Ｇ程序提供的方法进行计算．该方法的好处是采用迭代法估计潮汐参数、漂移和气象等时间序列的回归系数，并且采用贝叶斯信息原理来判断最优参数估计，从而提高了计算精度（Ｉｓｈｉｇｕｒｏ，Ｔａｍｕｒａ，１９８５；Ｔａｍｕｒａｅｔａｌ，１９９１）．其基本原理是将井水位时间序列进行如下分解：式中，Ｍ为波群数；αｍ和βｍ分别为第ｍ个波群在潮汐分波表中的始末位置；Ａｍｎ，ωｍｎ，φｍｎ分别为第ｍ波群中第ｎ个潮汐分量的理论振幅、角频率和初始相位；δｍ和Δφｍ为第ｍ波群待估算的潮汐振幅比和相位滞后；Ｄｒ（ｔ）为漂移（或长趋势）项；Ｒ（ｔ）为气象等因素引起的干扰项；ｈ（ｔ）为人为或仪器引起的阶变项；ε（ｔ）为观测噪声．

２．２观测资料的预处理井水位观测资料不仅包含潮汐信息，而且还包含一些其它的非潮汐信息．为提高潮汐计算的精度，在进行井水位潮汐观测资料的分析前，对观测数据进行适当合理的预处理是非常必要的．因为观测资料中的阶变、脉冲和缺数会影响潮汐振幅和相位的计算精度．为了消除人为或仪器因素引起的计算误差，计算前，我们对观测资料中人为引起的阶变和脉冲信号进行人工修整，对观测中由于仪器等问题引起的短时间缺数，采用Ｍａｔｌａｂ提供的３次样条内插方法进行内插处理．３潮汐特征分析为提取云南曲靖井水位的潮汐和相位动态变化信息，我们利用Ｂａｙｔａｐ－Ｇ所提供的方法，取计算窗长为７２０ｈ（３０天观测资料）、滑动步长为３６０ｈ，分别计算了Ｍ２、Ｓ２、Ｏ１、Ｓ１Ｋ１波和Ｍ３波的潮汐振幅、相位及计算误差．计算结果表明，Ｍ２波的潮汐振幅和相位计算误差小，变化相对平稳，而Ｓ２、Ｏ１、Ｓ１Ｋ１波和Ｍ３波的潮汐振幅和相位计算误差相对较大，稳定性也相对较差．由Ｍ２波的动态变化特征可以看出，Ｍ２波的潮汐振幅和相位并非不变的常数，而是随时间动态变化的．自１９９８年观测至２００１年底，Ｍ２波的潮汐振幅和相位变化基本是稳定的，潮汐振幅比主要维持在０．３４—０．４ＧＰａ之间．２００１年年底潮汐振幅比和相位均出现了明显的下降变化；之后于２００３年９月、２００５年１月、２００７年６月、２００８年６月、２００９年７月和２０１１年３月均出现了明显的上升和缓慢衰减过程．同样的变化过程在Ｓ２和Ｏ１波的相位和潮汐振幅动态变化曲线中也能看到，但是Ｍ

３波的变化并不明显．其中，２００１年和２００３年年底该井有抽水的影响，潮汐振幅和相位的下降变化可能与抽水有关，而之后的几次小幅度上升和缓慢的恢复过程则可能与中强地震的发生有一定的关系．这几次中强地震在图３中已用竖线标出，地震的基本情况见表１．由图３和表１可以看出，在几次能量密度较大的地震发生后，潮汐振幅比和相位均发生了明显的变化，而２００９年７月９日姚安ＭＳ６．０地震和２０１０年４月１４日青海玉树ＭＳ７．１地震的能量密度相对较小，其引起的变化也不是很显著；２００１年１１月１４日昆仑山口西ＭＳ８．１地震时该井正好受抽水影响，给分析带来一定困难．误差分析是评价计算参数可靠性的有效手段．不同潮汐波的均方根误差（ｒｏｏｔｍｅａｎｓｑｕａｒｅｅｒｒｏｒ，简写为ＲＭＳＥ）分析表明（图４），Ｍ２潮汐波的计算误差最小，Ｍ３潮汐波的计算误差最大，而Ｏ１波、Ｓ２波和Ｓ１Ｋ１波计算误差则位于Ｍ２与Ｍ３波之间，其中Ｍ２波潮汐振幅误差主要分布在０—０．０２ＧＰａ之间，相位误差主要分布在０°—１．５°之间；Ｓ２波的潮汐振幅误差主要分布在０—０．０６ＧＰａ之间，相位误差主要分布在１°—６°之间；Ｏ１波的潮汐振幅误差主要分布在０—０．０６ＧＰａ之间，相位误差主要分布在１°—１０°之间；Ｓ１Ｋ１波的潮汐振幅误差主要分布在０—０．２ＧＰａ之间，相位误差主要分布在５°—２５°之间；Ｍ３波潮汐振幅误差主要分布在０．１—０．２ＧＰａ之间，相位误差主要分布在１５°—５０°之间．计算误差的大小不仅与各潮汐波振幅有关，而且还受其它干扰因素的影响．Ｍ２波的计算误差之所以最小，主要是由于该潮汐波是井水位潮汐观测中振幅最大、影响因素最少的一个谐波；而Ｏ１波和Ｍ３波观测幅度相对较小，故计算误差有所增大；Ｓ２波和Ｓ１Ｋ１波除观测振幅比Ｍ２波小外，还受到气温和气压扰动的影响，故计算误差也有所增大．

４讨论

井水位潮汐不仅体现了体应变固体潮的变化，同时还受到气压潮的影响（Ｍａｔｓｕｍｏｔｏｅｔａｌ，２００３；Ｑｕｉｌｔｙ，Ｒｏｅｌｏｆｆｓ，１９９１；ｖａｎｄｅｒＫａｍｐ，Ｇａｌｅ，１９８３）．由于大气压往往受到温度的调制，因此与太阳有关的气压潮汐波Ｓ１Ｋ１和Ｓ２往往与热效应混合在一起，使得确定Ｓ１Ｋ１和Ｓ２波的潮汐信息变得复杂化．而月球引起的Ｏ１波和Ｍ２波受气压影响的程度较小，理论上源于月球的Ｍ２和Ｏ１波对于分析水位固体潮都是理想的，但由于Ｏ１波的幅度较小，信息往往掩埋在噪声中，而Ｍ２波幅度为水位潮汐信息中最大的一个谐波，因此Ｍ２波是分析井水位固体潮以及含水层特征最为有效的潮汐分量．井水位对固体潮的响应可以用渗流条件下的水流运动方程来表示．根据井水位潮汐现象统一数学物理方程的解析解，潮汐响应的幅度和相位是含水层导水系数、频率以及井孔半径的函数．对潮汐而言，其频率是固定的，如果在井孔没有变径的条件下，井水位潮汐幅度和相位的变化主要受含水层导水系数的影响（Ｂｒｅｄｅｈｏｅｆｔ，１９６７；Ｅｌｋｈｏｕｒｙｅｔａｌ，２００６；Ｈｓｉｅｈｅｔａｌ，１９８８；Ｒｏｊｓｔａｃｚｅｒ，Ｒｉｌｅｙ，１９９０；汪成民等，１９８８；张昭栋，王昌文，１９９１；张昭栋，郑金涵，１９８９；车用太等，２００６）．由图３可以看出，曲靖井水位Ｍ２波的潮汐振幅和相位出现了明显的变化，该变化体现了与井孔相连的含水层渗透系数的变化．地震能影响水文地质系统，从而引起水位、流量、化学量的变化，甚至引起泥火山和沙土液化等现象（Ｍａｎｇａ，Ｗａｎｇ，２００７）．然而很多类似的变化距离震中较远，在与断层相关的静态应力无法解释这些变化的情况下，很多研究者倾向于用地震波引起的动态应力变化来解释类似的变化．动态应力本身不能引起井水位的持续变化，但是动态应力能够清除裂隙中的障碍体，从而改变含水层的渗透系数并导致孔隙压力的重新分布（Ｂｒｏｄｓｋｙｅｔａｌ，２００３；Ｅｌｋｈｏｕｒｙｅｔａｌ，２００６；Ｒｏｅｌｏｆｆｓ，１９９８；Ｒｏｊｓｔａｃｚｅｒｅｔａｌ，１９９５；Ｒｏｊｓｔａｃｚｅｒ，Ｗｏｌｆ，１９９２；Ｗａｎｇ，２００７；Ｗａｎｇ，Ｃｈｉａ，２００８；Ｗａｎｇｅｔａｌ，２００９）．曲靖井水位在几次大震后出现的Ｍ２波潮汐振幅和相位变化，说明该井水位的同震效应与渗透系数的变化有关．地震波引起的含水层渗透系数变化的幅度与地震波峰值加速度呈正比例关系（Ｅｌｋ－ｈｏｕｒｙｅｔａｌ，２００６）．同时也有研究者注意到，井水位持续变化的幅度与地震波的峰值地面速度呈正比例关系（Ｒｏｅｌｏｆｆｓ，１９９８）．Ｗａｎｇ和Ｍａｎｇａ（２０１０）总结了地震波能量密度与地震引起的一些水文现象之间的关系，给出了地震波引起含水层渗透系数变化的能量密度下限约为１０－４Ｊ•ｍ－３．其中当能量密度大于１０－３Ｊ•ｍ－３时，地震波能量能更有效地清除裂隙中的障碍体，从而引起含水层渗透系数的变化更为显著．从云南曲靖井水位Ｍ２潮汐相位和振幅比的变化情况也可以清除地看出，在２００４年苏门答腊ＭＷ９．１，２００８年５月１２日四川汶川ＭＳ８．０，以及２０１１年３月１１日日本ＭＷ９．０地震引起的含水层渗透系数变化最为显著，但从能量密度与振幅比和相位差之间的关系来看，似乎能量密度与振幅比和相位差之间并没有明显的线性关系（图５ａ）．从同震响应引起振幅比和相位差变化的几次地震震级、震中距与能量密度之间的关系来看（图５ｂ），引起振幅比和相位差变化的能量密度下限为１０－３Ｊ•ｍ－３．这一结果与Ｗａｎｇ和Ｍａｎｇａ（２０１０）给出的结果是一致的．这可能说明，当地震的能量密度大于一定阀值时，地震波产生的能量能引起含水层参数发生变化，但这种变化似乎不是简单的线性关系．震后渗透系数的恢复过程体现了断层的愈合过程．１９９５年Ｈｙｏｇｏｋｅｎ－－Ｎａｎｂｕ地震发生后，为研究Ｎｏｊｉｍａ断层的愈合过程，Ｎｏｊｉｍａ项目每３年进行一次重复注水试验，以研究断层区渗透系数的变化情况．结果表明，震后３年内渗透系数至少下降了５０％（Ｋｉｔａｇａ－ｗａｅｔａｌ，２００２，２００７）．Ｍ２波相位的变化体现了井－含水层渗透系数的变化．Ｅｌｋｈｏｕｒｙ等（２００６）分析了地震引起的Ｍ２波相位变化情况，发现地震发生时，地震波引起渗透系数增加，之后逐渐恢复到一个稳定水平，说明了与含水层相连的裂隙的愈合过程．但从几次大震后振幅比和相位差的恢复情况看，并没有完全恢复到震前的水平，说明渗透系数在一定程度可能存在永久性的改变．同时，从２００１年抽水引起的振幅比和相位差变化情况可以看出，抽水引起的井水位下降在２００１年底就结束了，但潮汐因子和相位的低值一直持续到２００３年缅甸ＭＳ７．２地震的发生．这也说明，该井含水层参数改变后不能完全恢复到原来的水平．Ｌｉｕ和Ｍａｎｇａ（２００９）利用岩石物理实验模拟不排水条件下，在含孔隙和裂隙的砂岩中加入淤泥粒度的微小颗粒，测量渗透系数对重复震荡的响应情况．结果表明，重复震荡使裂隙中的微小颗粒固化，阻碍流体通道，从而降低渗透系数．Ｍｉｔｃｈｅｌｌ和Ｆａｕｌｋｎｅｒ（２００８）利用完整结晶岩石进行三轴压缩实验，以测量渗透系数和孔隙体积随差应力、应变、有效应力等因素的演化情况．结果表明，岩石在发生宏观破裂前，随差应力的增加渗透系数至少增加２个数量级；发生宏观破裂后的渗透系数比加压前增加了３个数量级；周期性加载情况下，随着差应力的增加渗透系数有所增加，随着差应力的减小渗透系数维持在相对较高的水平．据此认为，对能反映渗透系数变化的Ｍ２波的考察可能能捕捉到震前和震时的有用信息．抽水也可能引起含水层渗透系数的变化．含水层不仅受到上覆岩层的荷载作用，而且还受到来自水平方向的地应力作用．由于抽水导致储存在含水层中的水流出，引起含水层压力下降，使含水层发生压缩变形，孔隙空间随之减少，从而使其渗透系数也随之降低；流动阻力上升，流速降低，致使潮汐振幅和相位也随之变化．从图３可以看出，曲靖井Ｍ２波潮汐振幅和相位２００１年底出现了明显的下降，这可能与周边抽水有关．而２００３年底的振幅比和相位虽然总体呈缓慢上升变化趋势，但这种缓慢变化有别于几次大震引起的快速变化．这种缓慢上升可能是由于抽水引起的潮汐振幅比下降和缅甸地震引起的潮汐振幅比上升共同作用的结果．云南曲靖井水位自观测以来，曾出现过３次比较明显的大幅度下降．其中，２００１和２００３年的下降可能与附近地区的抽水有关，而２０１０年的下降没有找到干扰原因．统计３次井水位下降时段相应的水位潮汐振幅比和相位差变化（图６）可以看出，３次水位变化幅度与振幅比和相位差之间存在线性关系．这种关系说明，２０１０年井水位大幅度下降也引起了潮汐振幅比和相位差的变化．但由于样本数量较少，类似的关系是否具有普适性，将有待进一步证实．

５结论

通过本文的分析，得到以下结论和初步认识：１）井水位潮汐振幅和相位是反映含水层参数的有用信息．深入分析井水位潮汐变化特征，对了解含水层参数的正常动态特征，解释井水位微动态，评价井水位观测资料，以及捕捉地震前兆异常，无疑都具有重要的现实意义．２）云南曲靖井水位潮汐振幅和相位的动态变化可能与抽水和远距离强震等因素有关．其中，２００１和２００３年该井抽水期间，井水位潮汐振幅和相位均出现明显的下降变化．而２００４年苏门答腊ＭＷ９．１地震，２００８年四川汶川ＭＳ８．０地震，以及２０１１年日本ＭＷ９．０地震后，该井水位潮汐与相位振幅比和相位差均出现了明显的升高和缓慢的恢复过程．该变化证实了远距离大震地震波能引起含水层渗透系数增加的说法，其变化幅度可能与地震的能量密度有关，其恢复过程体现了裂隙的愈合过程．３）为评价井水位的观测质量和捕捉地震前兆信息，有必要对井水位潮汐的实时动态进行跟踪．同时为了评价计算结果的可靠性，对计算误差的考察也是必要的．尽管本文是对单一井孔的分析和讨论，但由于该井具有一定的代表性，其分析方法能为其它井孔潮汐的分析提供参考．

论文格式没有一个固定或都官方的样式，但是约定俗成，多看看不同的期刊的要求，结合自己的的需要再进行修改创作吧！