引言

地下水广泛赋存于地壳岩体的空隙之中，能在各种驱动力的作用下自由地流动． 这一特性，决定了地下水对地壳运动，尤其是对地震的孕育与发生会有较灵敏的响应(孙丽娟等，2004)，这也是地下流体各测项长期作为地震前兆观测手段之一的主要原因，但并非所有的地下水异常变化均反映的是地震前兆信息． 地下水动态变化包括宏观动态和微观动态. 宏观动态是由含水层水量增减变化而引起的地下水动态随时间的变化过程，而微观动态则是由含水岩体应力应变状态改变而引起的地下水动态随时间的变化过程(贾化周等，1995)． 因此，地下水动态的异常变化既可能是区域构造活动的作用结果，也可能是水文因素的干扰所致. 如何有效地识别水位异常的形成原因，是当前利用地震地下流体资料进行地震预测预报的关键科学问题之一．

水位观测是地下流体观测中一项重要的观测手段，地下水位观测值是在多种因素影响下产生的综合物理量． 在没有构造应力引起水位变化的条件下，水位的变化常与地下水补给、 泄流状态密切相关，如降雨和地下水开采等(车用太等，1993； 晏锐，2008； 杨明波等，2009). 如何排除这些干扰因素的影响，是确认一项前兆异常之前的首要工作． 经过多年的实践，学者们提出并尝试了不同的分析方法，如在分析降雨对水位的影响时，提出了响应函数法( Besbes，De Marsily，1984； Park Parker，2008)、 卷积滤波法(王尤培等，1996)及组合水箱模型(王旭升等，2010)等方法. 在分析地下水开采对水位的影响时，常采用经典的Dupuit模型、 Thiem模型和Theis模型(陈崇希，1966； 薛禹群等，1997； 薛禹群,谢春红，2007)等，这些方法模型在实际应用中均取得了一定的实用效果．

从地震预报角度来讲，当前我国地震监测预报实践面临的突出问题之一，是如何更加有效地排除干扰异常并提取出有效的前兆异常信息(车用太等，2011)． 以当前的研究水平，尚无法直接准确地识别出地震前兆异常信息，只能用“排除法”，即在尽可能全面并有效地排除各类干扰因素的前提下，来确定某项异常与地震前兆之间的相关性． 因此，当出现地下水动态异常变化时，异常干扰分析工作显得尤为重要，它是准确识别地震前兆异常信息的必要前提． 云南姚安井水位于2009年10月后出现了持续下降的异常变化，截止2011年底，最大下降幅度约3 m，这种大幅度的水位异常下降，是环境干扰所致还是地震前兆异常信息？ 对此，我们需要做深入地分析研究． 本文拟通过调查分析水位与降水量、 地下水开采之间的相关性，来探讨该井水位大幅度下降的形成原因．

1.   观测井概况及水位异常变化

云南姚安井位于楚雄市姚安县内，其南侧约20 km处有南华—楚雄断裂经过(图1a). 井孔海拔1 894 m，井深85.94 m，0—85 m用钢套管和塑管屏闭，其中0—32 m是直径为146 mm的钢套管，32—68 m是直径为127 mm的钢套管，68—85 m是直径为108 mm的硬塑料套管，85 m至井底为裸孔． 观测水位呈弱承压性质，为第三系褐煤层间承压水． 该井采用“十五”数字化水位观测设备，始测于2007年6月．

图  1  云南姚安水位观测井地理位置(a)与水位观测曲线(b)图

Figure  1.  Geographic location (a) and water level observation (b) of Yaoan well

下载: 全尺寸图片 幻灯片

姚安井水位于2009年10月开始快速下降，截止2011年底，其最大下降幅度达3 m(图1b)． 一般而言，引起水位下降的原因主要有两大类： 一是环境干扰，如降水量减小或地下水开采等会引起含水层内水量的减小，进而引起水位的下降变化； 另外，含水层固体骨架在区域构造应力作用下发生变形，进而引起的孔隙压减小也会导致观测井内水位的下降变化(赵利飞等，2002)． 因此，水位异常具体的形成原因，需收集相关辅助资料做深入地分析研究． 针对姚安井水位的异常下降变化，本文收集了姚安井附近区域降水量和地下水开采量资料，结合水位的长期变化动态，来分析其水位下降的形成原因．

2.   水位与降水量关系分析

引起含水层水量增减变化的主要因素为大气降水和地下水开采． 实际工作中，受施工抽水、 农田灌溉以及居民生活用水等因素的影响，要准确而全面地收集某区域的地下水开采量是非常困难的． 但是，收集局部区域大气降水量资料却是切实可行的，所以研究人员在进行水位宏观动态变化研究时，首先分析降水量对水位变化的影响，通过分析剔除了降水影响后的水位变化来进一步探讨其异常变化特征．

大气降水渗入补给是影响井水位变化的常见和主要因素，多数观测井的水位都有雨季上升、 旱季下降的年变规律. 降雨对水位观测的影响有两种最常见的形式: 一种是直接补给，即雨季降水时，井水位同步升高； 另一种是滞后补给，即井水位升高比大气降水滞后一段时间(白宝荣，付虹，2006). 云南姚安井即属于滞后补给型． 王旭升等(2010)针对水位滞后降水变化型观测井，提出了降雨-水位动态的组合水箱模型. 此模型可利用降水量变化来反演水位的动态变化，通过模拟值与实测值的对比分析剔除降水量对水位动态变化的影响． 该模型利用基于Gamma分布密度函数Γ(β)建立的单位脉冲响应函数γ(τ)来处理地下水补给的滞后延迟效果. 其数学表达式如下：

式中，t₀为一个时间因子，β为一个无量纲常数，γ(τ)表示单次降水在τ时刻形成的补给强度． 连续或间断的降水补给，可用上述单位脉冲响应函数的叠加来表示,即

式中，P_e(t)是t时刻产生的能够形成补给强度的有效降水量，参数T_jq反映井-含水层系统获得降水直接渗入补给的能力，T_iq反映井-含水层系统获得侧向补给的能力． 模型利用式(2)计算出某一时刻及之前降水所形成的累积补给强度R后，可预测下一时刻的水位值：

式中，H为水位值； 下标n表示当前时刻，n+1表示下一时刻； R_n是n时刻及之前降水量产生的累积补给强度； 参数z₀为排泄基准面的高度； T_h反映井-含水层系统的贮水能力．

该模型在利用降水量资料预测水位变化时，首先要利用正常动态变化时期的降水与水位资料来拟合求取z₀，β，t₀，T_h，T_iq和T_jq等6个参数，各参数均具有一定的物理含义，模型详细的推导见王旭升等(2010)． 在实际应用时，如果降水量和水位资料均为日值，则式(1)和(2)中的τ，t与n均以“天”为单位； 如果为月值，则以“月”为单位．

本文运用以上模型，模拟了姚安井水位的月均值变化曲线，以便更加客观地讨论降水量与水位变化之间的相关性特征(图2)． 模拟过程中6个参数是利用姚安井水位正常变化期间(2007—2008年)降水量与水位月均值拟合求得，其取值分别为z₀=0.8，β=2.2，t₀=2.1，T_h=0.1，T_iq=80，T_jq=13.6． 由图2可看出，在水位变化的正常动态时期模拟值与实测值能较好地吻合，模拟水位与实测水位的平均绝对误差为2.1 cm．

图  2  姚安井水位观测值与模拟值对比

Figure  2.  Comparison between observed and simulated water level in Yaoan well

下载: 全尺寸图片 幻灯片

图2中2009—2011年的水位模拟值是基于以上6个参数计算所得． 从模拟结果看，姚安井水位2009年6—9月的模拟值明显低于实测值. 这种水位的异常升高可能是与区域构造活动有关，是否与2009年7月9日姚安M_S6.0地震有关，尚有待进一步的研究． 2009年10月后实测水位急速下降，明显低于模拟水位，说明剔除降水量减小影响后的水位仍存在大幅下降的异常变化． 因此，仅从降水量补给角度看，姚安井2009年10月后出现的水位大幅下降与降水量减小并无直接相关性．

除了降水补给强度减弱会引起井水位的下降变化外，地下水过量开采也同样会引起井水位的下降变化(杨明波等，2009)． 由图2中姚安井附近区域降水量数据可看出，该地区2009—2011年降水总量连年下降，表明此段时间姚安地区气候相对较干旱． 因此，虽然姚安井水位下降与降水量减小无直接相关性，但由于同期气候干旱、 降水量连年降低，为满足居民生活用水、 农田灌溉等需要，会引起地下水开采量的增加． 姚安观测井深度不到100 m，很容易受地下水开采量的增加而导致水位出现大幅下降． 所以，姚安井2009年10月后所出现的水位大幅下降可能与同期气候干旱导致的地下水开采量增加有一定关系，这需要收集相关区域的地下水开采量资料做进一步地深入研究．

3.   水位与地下水开采关系分析

地下水作为重要的生产和生活资源，随着各地经济与社会飞速发展，开采的地区和开采井数量越来越多，开采强度也越来越大，地下水开采对地震地下水动态的干扰也变得越来越普遍和严重(车用太等，2011)． 据相关报道 中国国土资源网“姚安： 岩溶红层区打成5口深井”.[2010-03-31]http://www.clr.cn/front/read/read.asp?ID=191750. ，2009年底—2010年初，云南姚安地区遭遇百年不遇的旱情，持续数月无有效降雨． 为此，姚安县国土资源局组织相关部门在旱情最重、 人畜饮水最困难、 受灾人口相对集中的地区，进行了打井取水工程，先后在栋川、 光禄、 大龙口、 官屯、 弥兴等乡镇打井若干口． 截止2010年3月底，姚安县5口深井实现日出水量260 m³，解决了邻近群众的生活用水问题． 这些抽水井的抽水层位多在50 m以下深度. 姚安井相关资料 《云南省地震监测志》第三章第二节“楚雄彝族自治州地震监测台网”. 显示，附近区域埋深50—86 m层位主要为粘土层和褐煤层岩溶裂隙含水层，而姚安井水位观测层在85 m处． 因此，一方面这种大面积的干旱，会直接影响到地下水补给量的减小，另一方面地下水开采层位与姚安井观测层位的地下水之间存在一定的水力联系，周边抽水活动有可能对姚安井水位动态变化产生一定的影响．

从图2中姚安井水位与降水资料的对比可看出，水位的大幅下降出现在降水量减少的时段，而水位恢复则出现在降水量充裕的时段． 这种水位动态表现为“抽水过程中水位下降，停止抽水时水位恢复”的特征，表明姚安井水位的动态变化一定程度上与周边降水减小引起的地下水开采有关． 图3为姚安井附近区域地质图及周边开采井的空间分布图. 可以看出，各开采井与观测井同处白垩系砂岩、 泥岩夹泥灰岩岩层，构造位置上均处于牟定断裂与南华—楚雄断裂之间． 这种岩性以砂岩为主，区域内断层分布于相对比较简单的地质环境，更有利于地下水的渗流与水力联系． 各开采点与观测井之间的距离分别约为： 栋川2 km、 官屯5 km、 大龙口6 km、 光禄9 km、 弥兴15 km． 陈崇希(1966，2010)曾指出，在理想情况下，井孔抽水引起的水位下降范围是随时间延续发展的，假如没有因抽水而得到补给量的增加或排泄量的减少，那么影响范围将无限制地发展，达到整个含水层． 但实际情况下，由于自然界中潜水位或承压水头并非水平，即存在天然径流-补给，只要地下水开采量小于补给量就可以形成稳定井流． 车用太等(2011)指出，在基岩裂隙含水层中，水文地质条件较复杂地区，其干扰源(抽水点)对观测井水位动态产生影响的最大距离约为10 km． 因此，从空间上看，除了弥兴井外，其它4个抽水井均有可能对姚安观测井的水位动态产生不同程度的干扰．

图  3  姚安井附近区域地质图

Figure  3.  Regional geological map around Yaoan well

下载: 全尺寸图片 幻灯片

由上可知，姚安井水位与周边地下水开采之间存在着成因与空间上的相关性，但这尚不能确定二者之间存在必然的因果关系，还需要分析其时间与强度上的相关性． 为便于分析，本文假定姚安井含水层系统为均质各向同性的承压含水层，其地下水非稳定运动的基本微分方程可表示为(薛愚群等，1997)

式中，H为水头，H₀为初始水头值，t表示时间，K表示含水层的渗透系数，μ_s为含水层贮水率； W为单位体积流量，用以表示从单位体积含水层中流入(为正)或流出(为负)的水量，Γ₂为第二类边界条件． 上述模型从实用角度出发，假设含水层中的水流服从Darcy定律，渗透系数K、 贮水率μ_s不受含水层孔隙度变化(骨架变形)的影响．

式(4)加上相应的初始条件和边界条件，便可构成一个描述地下水流动体系的数学模型. 该模型计算时需输入含水层渗透系数K和贮水率μ_s的值，二者均可通过抽水试验得知. 但在实际应用中由于受客观条件的限制，无法进行抽水试验来获取含水层参数． 为此，很多学者提出了利用观测水位的固体潮调和分析法来获取含水层的介质参数(Munk Macdonald，1960； Bredehoeft，1967； Marine，1975； Hsieh et al，1987； 唐九安，1999； 张昭栋等，2002)，即通过固体潮调和分析法可得知含水层的贮水率μ_s和导水系数T. 而导水系数T与渗透系数K之间存在定量关系，即T＝KM，M为含水层等效厚度． 本文利用姚安井水位的M₂波固体潮调和分析结果，求得该井观测含水层的贮水率μ_s和导水系数T的逐月变化值(图4)． 从图4中可以看出，2009年后该井含水层贮水率和导水系数的变化相对比较稳定，说明姚安井含水层介质参数在水位大幅下降的过程中并未发生明显的变化，依据式(4)所建立的模型能基本反映实际的情况．

图  4  姚安井水位M₂波调和分析结果与含水层参数

Figure  4.  Water level harmonic analysis results of M₂ wave and aquifer parameters of Yaoan well

下载: 全尺寸图片 幻灯片

式(4)表示的微分方程，一般很难求得其解析解． 因此，各种各样的数值法被用来求得该式的近似解(薛禹群,谢春红，2007). 本文采用有限差分法(向后差分法)来求解． 计算过程中，除了给定渗透系数K和贮水率μ_s外，还需要给定模型的边界条件和初始条件． 考虑到姚安井观测含水层为承压含水层，其含水层上、 下均存在隔水层，故模型中垂向边界条件设为隔水边界． 为了探讨姚安井附近地区5口抽水井抽水的影响范围，模型中水平方向的边界条件设定为自由边界，即在抽水期间，降深漏斗之外的区域可视为定水头边界． 由于本文的关注点主要为抽水引起的水位下降，所以，模型中初始水头设定为统一值，其水位降深初始值均为0 m，即在抽水前含水层中的水头值处处相等． 另外，查阅相关资料 《云南省地震监测志》第三章第二节“楚雄彝族自治州地震监测台网”. 得知研究区含水层厚度约为50 m，其平均渗透系数K约在4×10^-6—6×10^-6 m/s之间. 由图4结果可知，姚安井含水层贮水率μ_s约为7×10^-7—9×10^-7/m之间．

姚安地区地下水开采井是为抗旱而打成，所以其地下水开采量与降水量之间应该存在一定的相关性，即降水量充足时段其开采量减小，而降水量减小时段其开采量会相应增强． 另外，姚安井水位多年动态特征表明，其水位与降水之间存在约2个月的延迟(图2)，属于滞后补给型观测系统． 那么，降水量减小引起的干旱在时间上也会存在一定程度的滞后效果，其开采量与降水量之间也应存在时间延迟． 因此，本文在模拟过程中，模型中抽水量数据的选取是在参考“5口深井最大出水量260 m³”的基础上，依据降水量变化而定，二者之间呈线性相关．

图5为姚安井水位实测值与模拟值对比曲线，图5a,b中含水层参数K和μ_s分别取其平均值5×10^-6 m/s和8×10^-7/m． 结合图2中降水量资料，图5a中的抽水量是依据降水量计算得出，二者在数值上满足P_抽水=50-10×P_降水的线性关系. 该图中2010年初模拟水位下降值较实测值低，说明该时间段实际抽水量要比本文模型中所赋的值偏小; 而2011年初的模拟水位下降值较实测值高，说明该时间段实际抽水量较模型中所赋的抽水量值要大． 笔者分析，本文模拟过程中的抽水活动是自2009年11月开始的，假定5口观测井同时抽水，且各个抽水井的抽水量相等． 实际情况是，姚安县内各抽水井的形成并出水是有时间先后顺序的，部分抽水井在2010年3月才开始投入使用，所以实际的抽水量要比本文模型中所赋的抽水量小，其水位的下降幅度也会较模拟值小． 2010年后，地方政府依“找水行动计划”继续打井并投入使用，从而导致2011年枯水期(降水量减小时段)抽水量较2010年增加的现象，所以实际的抽水量要比本文模型中所赋抽水量要大，其水位的下降幅度也会较模拟值大． 为了更接近实际情况，图5b中设定2010年11月前的各月抽水量是图5a中抽水量的80%，而2010年11月后的各月抽水量是图5a中抽水量的160%． 从水位模拟结果来看，调整抽水量数据后的模拟水位更逼近实测水位值．

图  5  姚安井水位观测值与模拟值对比

Figure  5.  Comparison between observed and simulated water level in Yaoan well

下载: 全尺寸图片 幻灯片

图6为姚安观测井附近周边约50 km范围内自2009年底开始抽水，截止2010年7月(图5a)和2011年7月(图5b)的模拟水位降落漏斗分布图(抽水量与图5b相同，含水层参数K和μ_s分别取5×10^-6 m/s和8×10^-7/m)． 从图6可以看出，5口抽水井在持续抽水约半年后，可对其周边地下水水头产生明显的影响，但其影响范围相对有限，水位降深值也较小． 随着抽水量增大和抽水时间增长后，其影响范围也逐渐扩大，水位降深值也随之加大． 在井群抽水情况下，各个抽水井抽水引起的降落漏斗之间存在相互叠加效应，可以说在整个水位降落漏斗区域，水位的下降变化不同程度上受到了各个抽水井的综合影响． 但是，由于各抽水井之间相隔距离不同，降落漏斗区内不同位置的水位降落值也存在一定的差异． 就姚安井来说，其水位的下降主要受离其较近的栋川井和官屯井的影响，其次是受大龙口井的影响．

图  6  2010年7月(a)和2011年7月(b)抽水引起的水位降落漏斗分布图

Figure  6.  Distribution of groundwater depression cone induced by pumping. (a) July 2010; (b) July 2011

下载: 全尺寸图片 幻灯片

由于受客观条件的限制，本文无法得到姚安井水位观测含水层渗透系数K和贮水率μ_s的确切值，只得到了二者的取值范围． 所以，为了分析含水层参数变化对计算结果的影响，在其取值范围内，本文分别计算了取不同含水层参数时的水位变化，计算结果见图7． 图7a中取K为定值5×10^-6 m/s，而μ_s分别取7×10^-7/m，8×10^-7/m和9×10^-7/m； 图7b中取μ_s为定值7×10^-7/m，而K分别取4×10^-6 m/s，5×10^-6 m/s和6×10^-6 m/s． 对比图7a，b可看出，在本文模型所取参数变化范围内，贮水率μ_s的波动对水位计算结果的影响不是很大，而渗透系数K的波动则对水位计算结果的影响相对较大． 但总体来看，在所取参数范围内，水位的模拟值与实测值能较好地吻合．

图  7  姚安井水位观测值与模拟值对比

Figure  7.  Comparison between observed and simulated water level in Yaoan well

下载: 全尺寸图片 幻灯片

结合以上水位变化与地下水开采活动的分析，可以看出，云南姚安井水位动态异常变化与地下水开采活动之间在成因、 空间、 时间和强度等4个方面均存在一定的相关性，基本判定为周边抽水引起的干扰异常．

4.   讨论与结论

引起地下水位异常变化的主要因素为构造活动和环境干扰，要判断水位异常是否包含地震前兆信息，首要工作需判定其异常是否存在干扰因素． 本文主要从降水影响和地下水开采活动两个方面，分析了姚安井水位2009年10月后出现异常下降的可能形成原因． 环境干扰方面，依据姚安井-含水层系统水文地质特征，建立了相应的数学物理模拟，分析了姚安井水位与降水量和地下水开采量之间的相关性特征，得出以下几点认识：

1) 通过分析姚安井水位与其降水量之间的定量关系发现，仅从降水量补给角度来看，姚安井降水量2009，2010和2011年的连续减小与2009年10月后出现的水位大幅度下降并无直接相关性． 但考虑到水位异常变化期间云南地区大面积的气候干旱因素，分析认为干旱导致的地下水开采量增加有可能是该异常变化的主要原因．

2) 依据姚安县内各开采点抽水量资料，通过对比模拟水位与实测水位，分析认为区域地下水开采引起的水位下降与实测水位下降之间存在成因、 空间、 时间和强度上的相关性． 虽然模拟结果与实测值在局部时间段内有所差异，但考虑到本文所构建地下水流模型的理想化与各井抽水量资料的完整性，模拟结果与实测值存在误差也在情理之中． 总体来看，姚安井水位2009年后出现的持续下降异常变化，与观测井所处区域遭遇的旱灾存在一定的相关性，周边打井抽水是水位下降的主要干扰源．

3) 综合分析认为姚安井及其附近区域降水量明显减小导致出现区域气候干旱，进而引起该区域地下水开采活动明显增强，而开采层与姚安井观测层基本属于同一层位，所以地下水开采引起的区域地下水位下降，是姚安井2009年10月水位异常下降的主要原因．

多年的监测实践表明，地下水动态的异常变化既有可能是震前的前兆异常，也有可能是某种环境干扰所致，及时识别与排除干扰异常，对于正确把握震情与成功识别前兆异常和预测至关重要(车用太等，2011)． 地下水动态干扰因素中最为常见的是大气降水的渗入补给与地下水开采干扰，且干扰方式多种多样，这给干扰排除工作造成了极大的困难. 除了要定性地确定干扰源与地下水动态异常之间的相关性，更关键的是要确定二者在空间、 时间和强度上的定量相关性，即异常出现的时间与干扰源作用时间、 异常幅度与干扰源作用强度之间的相关性特征分析．

本文在分析姚安井水位异常变化与大气降水、 地下水开采的相关性时，建立了一定的数学物理模型，通过数值模拟结果与实测值的对比来定量化地探讨地下水动态与干扰因素之间的相关性特征． 异常识别与分析过程中所采用的研究方法，虽然具有一定的局限性，但对今后地下水动态异常识别及判定工作具有一定的借鉴意义． 另外，采用数值模拟手段来进行异常分析过程中，关键是确定数学物理模型中各参数(如边界条件、 初始条件、 井-含水层介质参数)的合理性，各参数的取值均要依据实际的井-含水层系统来确定，任何一个参数的设定有偏差，都会导致模拟结果的失真，从而导致模拟值与实测值的对比毫无意义，异常判定结果就会有误． 即便各参数的取值合理，由于模型建立过程中的概念化、 理想化及其适用范围的局限性，模拟结果往往与实测值之间也不会完全吻合，这就需要研究人员在结果分析时进行必要的深入剖析．

审稿专家对本文提出的建设性意见，使作者对姚安井水位异常变化及其成因有了更深入地认识，在此表示衷心地感谢．