基于井水位潮汐效應(yīng)反演含水層參數(shù)的程序改進(jìn)

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

（1、防災(zāi)科技學(xué)院，河北 三河 065201 2、河北省地震動力重點(diǎn)實驗室，河北 三河 065201）

1 概述

含水層參數(shù)表征了含水層介質(zhì)的特征，其在水文地質(zhì)領(lǐng)域是非常重要的研究課題，同時也是研究各類水文地質(zhì)問題的基礎(chǔ)[1]。而含水層參數(shù)的計算除了傳統(tǒng)的抽水試驗和微水試驗以外，Hsieh等（1987）利用井水位與含水層孔壓變化之間的相位差推導(dǎo)了反演含水層導(dǎo)水系數(shù)的方法[2]?；诖死碚摶A(chǔ)，日本學(xué)者M(jìn).Ishiguro，Y.Tamura，T.Sato和M.Ooe利用Fortran77語言改進(jìn)開發(fā)了Baytap-G程序[3]。該程序利用貝葉斯建模程序來分析包含潮汐和其他變化的時間序列，以此可進(jìn)行：（1）估算振幅因子和相位差；（2）確定趨勢項和計算其頻譜；（3）對缺失數(shù)據(jù)進(jìn)行插值和對階梯變化量進(jìn)行估計；（4）粗略查找、識別異常數(shù)據(jù)；（5）計算ABIC值（Akaike's Bayesian Information Criterion）以此分析模型好壞[3]。將得到的相位差等參數(shù)手動輸入至Matlab[4]中即可計算出導(dǎo)水系數(shù)。此方法理論簡單易懂，但實際操作過程中需要手動輸入數(shù)據(jù)占用了一定的工作量。本次研究針對傳統(tǒng)計算方法步驟繁瑣，不利于研究者使用的問題對原有程序進(jìn)行了改進(jìn)和嫁接，使改進(jìn)程序?qū)崿F(xiàn)以Excel表格配置簡單參數(shù)，利用Matlab一鍵運(yùn)行即可得到所求參數(shù)，簡化操作步驟，提高研究者工作效率。

2 理論依據(jù)及原有計算方法

開放式承壓含水井水位響應(yīng)受地球固體潮擴(kuò)張引起的壓力水頭擾動，導(dǎo)致地下水的流動，當(dāng)含水層滲透率高時，那么井中的振蕩相位與受迫潮汐應(yīng)變幾乎協(xié)調(diào)一致；當(dāng)含水層滲透率低時，井中水位振蕩就較大的滯后于潮汐應(yīng)變。而實際情況位于兩者之間，因此在含水層潮汐應(yīng)變和井水位響應(yīng)之間存在相位差（時間滯后）[5]。（Elkhoury et al.，2007）分別計算出井水位和體應(yīng)變的相位差，則井水位與體應(yīng)變的相位差之差等于潮汐徑向流所產(chǎn)生的相位差；當(dāng)井水位潮汐振幅和相位差主要受潮汐徑向流影響，井水位潮汐相位差等于徑向流相位差。當(dāng)徑向流相位差已知時，可計算不同貯水系數(shù)情況下徑向?qū)禂?shù)。[6]

通過重力固體潮計算含水層參數(shù)，科研人員首先需要對原始水位數(shù)據(jù)進(jìn)行預(yù)處理包括水位的格式處理、異常值的處理、去除趨勢項以及分析數(shù)據(jù)并手工提取Baytap-G計算所需的參數(shù)，之后將處理好的井水位數(shù)據(jù)以及所需參數(shù)按照Baytap-G要求輸入軟件計算得到井水位振幅與相位差，最后篩選提取Baytap-G計算結(jié)果中的相位數(shù)據(jù)，轉(zhuǎn)化為絕對值數(shù)據(jù)之后輸入Matlab腳本中計算導(dǎo)水系數(shù)。整個方法需要研究者先對水位數(shù)據(jù)進(jìn)行預(yù)處理，之后再進(jìn)行井水位振幅、相位差以及導(dǎo)水系數(shù)的計算，過程繁瑣且易出現(xiàn)人為錯誤。

3 嫁接程序主要步驟

鑒于原有計算方法較繁瑣，本次研究對原有程序進(jìn)行了嫁接改進(jìn)，應(yīng)用Matlab程序?qū)Ω鞴δ艿膶崿F(xiàn)及對Baytap-G進(jìn)行調(diào)用的封裝，并且實現(xiàn)自動處理輸入數(shù)據(jù)格式和自動提取輸出的結(jié)果數(shù)據(jù)，一鍵運(yùn)行即得到所需要的計算結(jié)果。流程圖如圖1。

圖1 改進(jìn)程序操作流程圖

本次程序嫁接主要在Matlab中實現(xiàn)以下功能：讀取配置及數(shù)據(jù)文件、循環(huán)處理功能、臨時數(shù)據(jù)清理、異常數(shù)據(jù)處理、趨勢項消除、生成Baytap控制文件和數(shù)據(jù)文件、調(diào)用Baytap程序、Baytap運(yùn)行結(jié)果解析、導(dǎo)水系數(shù)計算、輸出最終結(jié)果。

4 案例測試

本次測試采用辛莊臺孔隙介質(zhì)的井進(jìn)行分析。井孔參數(shù)如下：

辛莊井位于滄縣隆起上的白塘口凹陷內(nèi)。1968年開始人工觀測水位，后轉(zhuǎn)為儀器觀測。2001年進(jìn)行了數(shù)字化改造并與2002年投入使用。井口標(biāo)高4m，井孔深度648.12 m，觀測層巖性為砂質(zhì)黏土夾薄層細(xì)沙、黏土、中砂，地下水埋藏類型為承壓水。

4.1 辛莊井井水位數(shù)據(jù)處理結(jié)果

將辛莊井2008年全年靜水位數(shù)據(jù)同樣利用改進(jìn)程序計算導(dǎo)水系數(shù)，計算過程與良鄉(xiāng)井計算過程一致，計算結(jié)果如表1所示。

表1 改進(jìn)程序?qū)禂?shù)結(jié)果表

同樣的數(shù)據(jù)利用傳統(tǒng)方法計算M2波相位數(shù)據(jù)見表2，將其轉(zhuǎn)化為絕對值數(shù)據(jù)后寫到導(dǎo)水系數(shù)計算的腳本中并運(yùn)行計算，導(dǎo)水系數(shù)計算結(jié)果見表3。

表2 潮汐異常處理軟件處理結(jié)果M2波相位數(shù)據(jù)表

表3 傳統(tǒng)方法導(dǎo)水系數(shù)計算結(jié)果表

4.2 處理結(jié)果對比分析

通過辛莊井2008全年靜水位數(shù)據(jù)，分別利用本文改進(jìn)程序與傳統(tǒng)方法計算導(dǎo)水系數(shù)，發(fā)現(xiàn)兩者計算結(jié)果基本相同但仍有細(xì)微差別（見表4、圖2）。

圖2 辛莊井導(dǎo)水系數(shù)計算結(jié)果對比柱狀圖

表4 辛莊井導(dǎo)水系數(shù)計算結(jié)果對比表

本文改進(jìn)程序與傳統(tǒng)計算程序采用相同的理論實現(xiàn)，而計算結(jié)果產(chǎn)生細(xì)微差別，但誤差很小在可接受范圍內(nèi)可忽略不計，基于科學(xué)的嚴(yán)謹(jǐn)性對于誤差的產(chǎn)生我們做了細(xì)致的分析發(fā)現(xiàn)，誤差產(chǎn)生的原因在于兩者Baytap生成的數(shù)據(jù)文件中的井水位數(shù)據(jù)的精確度不同，在井水位單位統(tǒng)一為米時，本文改進(jìn)程序生成數(shù)據(jù)精度為小數(shù)點(diǎn)后四位，傳統(tǒng)計算程序限于軟件限制生成數(shù)據(jù)精度為小數(shù)點(diǎn)后三位，本文計算數(shù)據(jù)精度更高。除此之外兩個程序在理論計算方面再無其他差別，在許多學(xué)者多次運(yùn)用并認(rèn)可傳統(tǒng)計算方法的基礎(chǔ)上本文改進(jìn)程序計算導(dǎo)水系數(shù)結(jié)果無誤。

5 結(jié)論

本文依據(jù)科研人員在開展井水位的固體潮效應(yīng)反演含水層參數(shù)研究中遇到的實際問題，對原有計算程序進(jìn)行了改進(jìn)，得到以下結(jié)論：（1）本文利用Matlab控制整個流程，調(diào)用各個功能模塊對傳統(tǒng)方法的不同業(yè)務(wù)進(jìn)行處理。改進(jìn)程序?qū)⒃杏嬎惴椒ǚ爆嵉牟襟E簡化，使研究者可以配置完簡單參數(shù)后，在Matlab一個程序中即可一鍵獲得所需參數(shù)。（2）通過對良鄉(xiāng)臺（裂隙介質(zhì)井）、辛莊臺（孔隙介質(zhì)井）兩口不同介質(zhì)井的井水位使用傳統(tǒng)計算方法與改進(jìn)程序計算導(dǎo)水系數(shù)，對比發(fā)現(xiàn)由于兩者Baytap生成的數(shù)據(jù)文件中的井水位數(shù)據(jù)的精確度不同使得計算結(jié)果有細(xì)微差別，誤差在可接受范圍內(nèi)。（3）本次改進(jìn)程序在不降低計算精度的前提下，極大的提高了計算效率，降低了計算時由于步驟繁瑣而出現(xiàn)的人為誤差，為研究人員進(jìn)行相關(guān)方面科學(xué)研究提供了極大的便利。