基于Matlab求解穩(wěn)定流的水文地質(zhì)參數(shù)

吳家杰，史麗娜

（江西省煤田地質(zhì)局一九五地質(zhì)隊(duì)，江西 南昌 330001）

1 概述

導(dǎo)水系數(shù)（T）和滲透系數(shù)（K）是進(jìn)行地下資源論證、地下水資源評(píng)價(jià)和地下水?dāng)?shù)值模擬的水文地質(zhì)參數(shù)。抽水試驗(yàn)法是確定含水層參數(shù)的主要方法。

MATLAB是矩陣實(shí)驗(yàn)室（Matrix Laboratory）的簡(jiǎn)稱，是美國(guó)MathWorks公司出品的商業(yè)數(shù)學(xué)軟件，可用于算法開發(fā)、數(shù)據(jù)可視化、數(shù)據(jù)分析以及數(shù)值計(jì)算。采用MATLAB軟件，不僅可以快速進(jìn)行大量數(shù)據(jù)的科學(xué)計(jì)算，還可以克服人工配線過程中產(chǎn)生的人工誤差。與其它軟件相比，MATLAB操作比較容易，編程相對(duì)較簡(jiǎn)單。本文依據(jù)梅山礦區(qū)抽水試驗(yàn)資料，采用MATLAB 程序來計(jì)算承壓含水層導(dǎo)水系數(shù)T 和滲透系數(shù)K，取得了較好的結(jié)果。

2 抽水試驗(yàn)

2.1 水文地質(zhì)概況

梅山礦區(qū)（南區(qū)）范圍內(nèi)分布的主要地層由老至新為：二疊系中統(tǒng)茅口組（P2m）、三疊系上統(tǒng)安源組紫家沖段（T3a1）、三疊系上統(tǒng)安源組龍王寨段（T3a2）、三疊系上統(tǒng)安源組高家段（T3a3）、侏羅系下統(tǒng)水北組（J1sh）、白堊系南雄組（K2n）、第四系（Q）。主要含水層為白堊系礫巖裂隙巖溶含水層、安源組紫家沖段底部礫巖含水層，茅口組灰?guī)r巖溶含水層，其它可視為隔水層或相對(duì)隔水層。

2.2 抽水試驗(yàn)情況

梅山礦區(qū)（南區(qū)）在詳查階段共進(jìn)行了二孔二層次的穩(wěn)定流抽水試驗(yàn)工作，兩個(gè)水文地質(zhì)鉆孔施工過程中均嚴(yán)格按照水文地質(zhì)鉆探施工要求進(jìn)行施工，對(duì)含水層采用清水鉆進(jìn)，并采用套管隔水，止水效果好，鉆探質(zhì)量等級(jí)為甲級(jí)；抽水試驗(yàn)前統(tǒng)一進(jìn)行含水層鹽化測(cè)井，測(cè)井成果準(zhǔn)確可靠，測(cè)井質(zhì)量等級(jí)為甲級(jí)；抽水試驗(yàn)嚴(yán)格遵照《煤礦床水文地質(zhì)、工程地質(zhì)及環(huán)境地質(zhì)勘查評(píng)價(jià)標(biāo)準(zhǔn)》進(jìn)行，抽水質(zhì)量均達(dá)合格以上，鉆孔質(zhì)量等級(jí)綜合評(píng)定均為甲級(jí)。

水1孔：該孔于2010年3月26日開孔，2010年6月3日終孔，對(duì)茅口組灰?guī)r進(jìn)行抽水。茅口組灰?guī)r抽水試驗(yàn)深度為469.68～530.44m，抽水延續(xù)時(shí)間107h（其中穩(wěn)定時(shí)間31h），共進(jìn)行了3個(gè)降深的試驗(yàn)工作，水位降深誤差分別為0.05%、0.02%、0.03%，流量誤差為0，流量、水位關(guān)系曲線為指數(shù)型，綜合評(píng)定質(zhì)量為優(yōu)質(zhì)。

水2孔：該孔于2010年4月18日開孔，2010年9月23日終孔，對(duì)茅口組灰?guī)r進(jìn)行抽水。茅口組灰?guī)r抽水試驗(yàn)深度為572.98～635.53m，抽水延續(xù)時(shí)間62.5h（其中穩(wěn)定時(shí)間28h），共進(jìn)行了3個(gè)降深的試驗(yàn)工作，水位降深誤差分別為0.1%、0.08%、0.18%，流量誤差為0，流量、水位關(guān)系曲線為直線型，綜合評(píng)定質(zhì)量為合格。

2個(gè)孔兩層抽水試驗(yàn)成果見表1。

表1 篩分和破碎設(shè)備

3 實(shí)際應(yīng)用

3.1 數(shù)據(jù)準(zhǔn)備

根據(jù)裘布依自流完整井公式，滲透系數(shù)K 和影響半徑R 可按照2個(gè)觀測(cè)孔及一個(gè)觀測(cè)孔、單孔抽水（其中R=10SwK0.5，K 的單位m／d）等三種形式求解，本次均為單孔抽水試驗(yàn)，未布設(shè)觀測(cè)孔，所以采用單孔抽水求解滲透系數(shù)和影響半徑。

利用裘布依自流完整井公式（K=0.366×Q×lg（R／r）／S×M），可將2個(gè)抽水孔資料分別整理（見表2）。

表2 抽水試驗(yàn)資料

3.2 MATLAB解算

傳統(tǒng)的求解滲透系數(shù)方法為迭代法，此法較為繁瑣，先假設(shè)一數(shù)值，多次迭代，直至相近為止。而MATLAB軟件中有一函數(shù)單變量函數(shù)求零值fzero，函數(shù)描述為若果x0是標(biāo)量，x=fzero（fun，x0）試圖求解函數(shù)fun在x0附近的值。由fzero返回的值x就是接近函數(shù)fun改變符號(hào)處的值，若果搜索失敗則返回NaN。

在MATLAB命令窗口中，依據(jù)表2中6個(gè)公式分別輸入：

x=fzero（＇x-1.4741×log10（x）-10.1833＇，12）；

x=fzero（＇x-1.1639×log10（x）-8.4411＇，12）；

x=fzero（＇x-0.9657×log10（x）-7.2323＇，12）；

x=fzero（＇x-21.6052×log10（x）-161.4252＇，177）；

x=fzero（＇x-23.1101×log10（x）-170.7527＇，177）；

x=fzero（＇x-24.4072×log10（x）-186.1250＇，177）；

x=fzero（＇x-1.4741×log10（x）-10.1833＇，12）。

解算結(jié)果見圖1，分別求出水1孔的3次降深的滲透系數(shù)K1、K2、K3分別為11.7613m／d、9.5835m／d、8.1101 m／d；水2 孔的3 次降深的滲透系數(shù)K1、K2、K3 分別為211.6716m／d、225.1171m／d、244.4124m／d。

圖1 MATLAB軟件滲透系數(shù)解算結(jié)果

MATLAB 求解滲透系數(shù)后，依據(jù)公式R=10 SwK0.5，求出各孔各次降深對(duì)應(yīng)的影響半徑R，水1孔的3次降深的影響半徑R1、R2、R3分別為537 m、720 m、870 m；水2孔的3次降深的R1、R2、R3分別為3722 m、3065 m、2653m。

2個(gè)抽水鉆孔各個(gè)水文地質(zhì)參數(shù)見表3。

表3 某礦區(qū)抽水試驗(yàn)水文地質(zhì)參數(shù)一覽

3.3 與傳統(tǒng)迭代法對(duì)比

MATLAB函數(shù)求解水文地質(zhì)參數(shù)與傳統(tǒng)迭代法求解水文地質(zhì)參數(shù)結(jié)果見表4，從表中可以看出，水1孔兩種方法求解水文地質(zhì)參數(shù)數(shù)值結(jié)果誤差不大，而水2孔兩種方法求解水文地質(zhì)參數(shù)數(shù)值結(jié)果產(chǎn)生一定誤差，可能導(dǎo)致估算涌水量產(chǎn)生較大偏差。

表4 MATLAB函數(shù)與傳統(tǒng)迭代法求解結(jié)果

4 結(jié)語

以往傳統(tǒng)求解抽水試驗(yàn)水文地質(zhì)參數(shù)（K、R），一般用迭代法，在實(shí)際過程應(yīng)用偏繁瑣，而通過裘布依自流完整井公式和影響半徑經(jīng)驗(yàn)公式，化為一元方程，利用MATLAB軟件fzero函數(shù)求解，優(yōu)點(diǎn)在于計(jì)算過程簡(jiǎn)單，產(chǎn)生誤差的概率小，精度高。而傳統(tǒng)迭代法先任意取值，一直迭代，不斷重復(fù)計(jì)算，直至數(shù)值接近目標(biāo)值，實(shí)際生產(chǎn)過程容易出差，且難以發(fā)現(xiàn)。因而利用MATLAB 函數(shù)求解水文地質(zhì)參數(shù)在實(shí)際生產(chǎn)過程運(yùn)用比傳統(tǒng)迭代法更為方便、可靠。

〔1〕李俊庭，王愈吉.地下水動(dòng)力學(xué)〔M〕.北京：地質(zhì)出版社，1986.

〔2〕石志遠(yuǎn)，許 模.抽水試驗(yàn)資料的遺傳算法反演水文地質(zhì)參數(shù)〔J〕.四川地質(zhì)學(xué)報(bào)，2011，31（4）：466-469，480.

〔3〕劉衛(wèi)國(guó).MATLAB升序設(shè)計(jì)與應(yīng)用（第二版）〔M〕.北京：高等教育出版社，2006.

〔4〕徐東艷，孟曉剛.MATLAB函數(shù)庫(kù)查詢辭典〔M〕.北京：中國(guó)鐵道出版社，2005.