基于井水位潮汐效应反演含水层参数的程序改进

魏海滨 谷洪彪* 翟泽宇 朱柏正 徐邑荣

（1、防灾科技学院，河北 三河 065201 2、河北省地震动力重点实验室，河北 三河 065201）

1 概述

含水层参数表征了含水层介质的特征，其在水文地质领域是非常重要的研究课题，同时也是研究各类水文地质问题的基础[1]。而含水层参数的计算除了传统的抽水试验和微水试验以外，Hsieh等（1987）利用井水位与含水层孔压变化之间的相位差推导了反演含水层导水系数的方法[2]。基于此理论基础，日本学者M.Ishiguro，Y.Tamura，T.Sato和M.Ooe利用Fortran77语言改进开发了Baytap-G程序[3]。该程序利用贝叶斯建模程序来分析包含潮汐和其他变化的时间序列，以此可进行：（1）估算振幅因子和相位差；（2）确定趋势项和计算其频谱；（3）对缺失数据进行插值和对阶梯变化量进行估计；（4）粗略查找、识别异常数据；（5）计算ABIC值（Akaike's Bayesian Information Criterion）以此分析模型好坏[3]。将得到的相位差等参数手动输入至Matlab[4]中即可计算出导水系数。此方法理论简单易懂，但实际操作过程中需要手动输入数据占用了一定的工作量。本次研究针对传统计算方法步骤繁琐，不利于研究者使用的问题对原有程序进行了改进和嫁接，使改进程序实现以Excel表格配置简单参数，利用Matlab一键运行即可得到所求参数，简化操作步骤，提高研究者工作效率。

2 理论依据及原有计算方法

开放式承压含水井水位响应受地球固体潮扩张引起的压力水头扰动，导致地下水的流动，当含水层渗透率高时，那么井中的振荡相位与受迫潮汐应变几乎协调一致；当含水层渗透率低时，井中水位振荡就较大的滞后于潮汐应变。而实际情况位于两者之间，因此在含水层潮汐应变和井水位响应之间存在相位差（时间滞后）[5]。（Elkhoury et al.，2007）分别计算出井水位和体应变的相位差，则井水位与体应变的相位差之差等于潮汐径向流所产生的相位差；当井水位潮汐振幅和相位差主要受潮汐径向流影响，井水位潮汐相位差等于径向流相位差。当径向流相位差已知时，可计算不同贮水系数情况下径向导水系数。[6]

通过重力固体潮计算含水层参数，科研人员首先需要对原始水位数据进行预处理包括水位的格式处理、异常值的处理、去除趋势项以及分析数据并手工提取Baytap-G计算所需的参数，之后将处理好的井水位数据以及所需参数按照Baytap-G要求输入软件计算得到井水位振幅与相位差，最后筛选提取Baytap-G计算结果中的相位数据，转化为绝对值数据之后输入Matlab脚本中计算导水系数。整个方法需要研究者先对水位数据进行预处理，之后再进行井水位振幅、相位差以及导水系数的计算，过程繁琐且易出现人为错误。

3 嫁接程序主要步骤

鉴于原有计算方法较繁琐，本次研究对原有程序进行了嫁接改进，应用Matlab程序对各功能的实现及对Baytap-G进行调用的封装，并且实现自动处理输入数据格式和自动提取输出的结果数据，一键运行即得到所需要的计算结果。流程图如图1。

图1 改进程序操作流程图

本次程序嫁接主要在Matlab中实现以下功能：读取配置及数据文件、循环处理功能、临时数据清理、异常数据处理、趋势项消除、生成Baytap控制文件和数据文件、调用Baytap程序、Baytap运行结果解析、导水系数计算、输出最终结果。

4 案例测试

本次测试采用辛庄台孔隙介质的井进行分析。井孔参数如下：

辛庄井位于沧县隆起上的白塘口凹陷内。1968年开始人工观测水位，后转为仪器观测。2001年进行了数字化改造并与2002年投入使用。井口标高4m，井孔深度648.12 m，观测层岩性为砂质黏土夹薄层细沙、黏土、中砂，地下水埋藏类型为承压水。

4.1 辛庄井井水位数据处理结果

将辛庄井2008年全年静水位数据同样利用改进程序计算导水系数，计算过程与良乡井计算过程一致，计算结果如表1所示。

表1 改进程序导水系数结果表

同样的数据利用传统方法计算M2波相位数据见表2，将其转化为绝对值数据后写到导水系数计算的脚本中并运行计算，导水系数计算结果见表3。

表2 潮汐异常处理软件处理结果M2波相位数据表

表3 传统方法导水系数计算结果表

4.2 处理结果对比分析

通过辛庄井2008全年静水位数据，分别利用本文改进程序与传统方法计算导水系数，发现两者计算结果基本相同但仍有细微差别（见表4、图2）。

图2 辛庄井导水系数计算结果对比柱状图

表4 辛庄井导水系数计算结果对比表

本文改进程序与传统计算程序采用相同的理论实现，而计算结果产生细微差别，但误差很小在可接受范围内可忽略不计，基于科学的严谨性对于误差的产生我们做了细致的分析发现，误差产生的原因在于两者Baytap生成的数据文件中的井水位数据的精确度不同，在井水位单位统一为米时，本文改进程序生成数据精度为小数点后四位，传统计算程序限于软件限制生成数据精度为小数点后三位，本文计算数据精度更高。除此之外两个程序在理论计算方面再无其他差别，在许多学者多次运用并认可传统计算方法的基础上本文改进程序计算导水系数结果无误。

5 结论

本文依据科研人员在开展井水位的固体潮效应反演含水层参数研究中遇到的实际问题，对原有计算程序进行了改进，得到以下结论：（1）本文利用Matlab控制整个流程，调用各个功能模块对传统方法的不同业务进行处理。改进程序将原有计算方法繁琐的步骤简化，使研究者可以配置完简单参数后，在Matlab一个程序中即可一键获得所需参数。（2）通过对良乡台（裂隙介质井）、辛庄台（孔隙介质井）两口不同介质井的井水位使用传统计算方法与改进程序计算导水系数，对比发现由于两者Baytap生成的数据文件中的井水位数据的精确度不同使得计算结果有细微差别，误差在可接受范围内。（3）本次改进程序在不降低计算精度的前提下，极大的提高了计算效率，降低了计算时由于步骤繁琐而出现的人为误差，为研究人员进行相关方面科学研究提供了极大的便利。