基于井水位潮汐效應反演含水層參數的程序改進

魏海濱 谷洪彪* 翟澤宇 朱柏正 徐邑榮

（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）本次改進程序在不降低計算精度的前提下，極大的提高了計算效率，降低了計算時由于步驟繁瑣而出現的人為誤差，為研究人員進行相關方面科學研究提供了極大的便利。