基于分段函數(shù)的潛水井降深曲線在工程實(shí)踐中的應(yīng)用

朱悅璐，熊俊飛，鄧 煒

（南昌工程學(xué)院 水利與生態(tài)工程學(xué)院，江西 南昌 330099）

1 引 言

地下水向井運(yùn)動(dòng)在礦山排水、灌區(qū)農(nóng)業(yè)灌溉、地下水資源評(píng)價(jià)中經(jīng)常會(huì)遇到［1-3］。 自Dupuit［4］推導(dǎo)出潛水向井穩(wěn)定流公式以來，計(jì)算降深曲線、預(yù)測(cè)地下水位變動(dòng)等，一直是該領(lǐng)域的熱點(diǎn)［5-7］。 中國(guó)地質(zhì)大學(xué)陳崇希教授團(tuán)隊(duì)與南京大學(xué)薛禹群院士團(tuán)隊(duì)在2010—2014 年就穩(wěn)定流Dupuit 公式滲流場(chǎng)影響半徑問題進(jìn)行多次商榷［8-10］。 常云霞［11］利用Visual Modflow 軟件，模擬分析了礦井周邊區(qū)域地下水降深曲線變化趨勢(shì)。蘭盈盈［12］采用GMS（Groundwater Model System）模擬軟件，對(duì)2015—2028 年大坳灌區(qū)地下水位進(jìn)行動(dòng)態(tài)預(yù)測(cè)，結(jié)果表明新建灌區(qū)對(duì)地下水位影響較小，不會(huì)造成土壤鹽漬化問題。 閆峭［13］利用FlowHeat 1.0，對(duì)西安市某小區(qū)熱力井在各種布井方案和抽灌條件下的滲流場(chǎng)、溫度場(chǎng)進(jìn)行了模擬，并優(yōu)化了監(jiān)測(cè)井網(wǎng)的布局方案。 Pulat 等［14］、Michopoulos 等［15］模擬了抽水井熱源泵系統(tǒng)在不同降深下的水熱平衡問題。

地下水向井運(yùn)動(dòng)的研究成果已較為豐富，但仍存在局限性，如潛水井空間徑向流場(chǎng)的復(fù)雜分布，很難求得解析解。 為解決該問題，研究者大多假設(shè)向井流動(dòng)的潛水流近似水平，忽略垂向分速度，這一近似處理手段導(dǎo)致地下水計(jì)算軟件的底層框架大多演變?yōu)閭鹘y(tǒng)的Dupuit 假設(shè)，將潛水運(yùn)動(dòng)方程簡(jiǎn)化為一維或二維Boussinesq 方程。 在該假設(shè)下經(jīng)驗(yàn)公式存在的問題：①靠近井壁處，降深曲線計(jì)算誤差較大，不能用于實(shí)際工程；②當(dāng)r ＞1.5H0（r為與抽水井的距離，H0為潛水初始水頭）時(shí)，理論公式才與實(shí)際相符。 對(duì)于①，有學(xué)者提出井壁區(qū)域滲流場(chǎng)計(jì)算的相關(guān)問題，例如應(yīng)用有限元三角形相似原理［16］，在穩(wěn)定滲流場(chǎng)中確定一個(gè)能量損失率最小點(diǎn)，取井壁處的固定水頭點(diǎn)為溢出點(diǎn)，求解r ＜1.5H0這一區(qū)域的降深曲線，該方案迭代速度快，但仍有計(jì)算精度不足等缺陷；對(duì)于②，當(dāng)r ＞1.5H0時(shí)，實(shí)際潛水降深與理論計(jì)算值也有很大差距，Dupuit公式仍不可直接應(yīng)用，但目前對(duì)此區(qū)間計(jì)算精度進(jìn)行修正的研究較少。

基于上述問題，筆者利用16 口觀測(cè)井實(shí)測(cè)數(shù)據(jù)，在穩(wěn)定抽水條件下擬合方程，生成計(jì)算曲線，并采用試算法求得Dupuit 理論公式在實(shí)際工程中應(yīng)用的精度分界點(diǎn)，以該點(diǎn)為區(qū)間臨界點(diǎn)，將計(jì)算曲線與理論曲線組成分段函數(shù)，最終得出經(jīng)過修正后的潛水滲流場(chǎng)。該方案可有效減小傳統(tǒng)降深曲線在全區(qū)間的計(jì)算誤差，為同類場(chǎng)地條件下確定計(jì)算分界點(diǎn)、優(yōu)化觀測(cè)井布置等提供參考。

2 傳統(tǒng)Dupuit 公式及其不足

潛水井向井運(yùn)動(dòng)是潛水含水層中地下水運(yùn)動(dòng)的一個(gè)特例，因此其動(dòng)力學(xué)特性滿足潛水含水層中地下水運(yùn)動(dòng)的基本微分方程及假設(shè)條件，即Dupuit 假設(shè)：對(duì)于垂直的二維平面上任意一點(diǎn)（x，z），當(dāng)潛水流垂向坡角θ很?。撍蛔兓苄。r(shí)，可以忽略垂向分速度vz而僅考慮水平分速度vx（見圖1）。

圖1 Dupuit 假設(shè)

潛水井地下水向井運(yùn)動(dòng)的Dupuit 公式為

式中：H為潛水水頭；h為潛水含水層厚度；x、z分別為水平向、垂向坐標(biāo)。

式（1）為直角坐標(biāo)系下的標(biāo)準(zhǔn)微分方程，為方便與抽水井實(shí)際工況匹配，常用的方法是將式（1）轉(zhuǎn)化為柱坐標(biāo)：

大多數(shù)研究者認(rèn)為Dupuit 假設(shè)沒有考慮水流在入井時(shí)的能量損失，因此在抽水井壁附近的理論降深和實(shí)測(cè)降深差異較大［17-19］。 當(dāng)邊界條件為：r ＝rw、h ＝hw和r ＝Rw、h ＝H0（rw為抽水井半徑，Rw為抽水井影響半徑，hw為抽水井內(nèi)水頭）時(shí)，按文獻(xiàn)［17］對(duì)式（2）進(jìn)行積分運(yùn)算，得到潛水井Dupuit 理論曲線表達(dá)式：

潛水含水層厚度h可表示為觀測(cè)井到抽水井距離r的函數(shù)h ＝f（r） ，當(dāng)潛水初始水頭H0已知時(shí)，由Dupuit 公式得出的潛水井附近的水位降深計(jì)算公式為

式中：s為水位降深。

很多實(shí)際工程算例表明，直接使用式（4）會(huì)導(dǎo)致計(jì)算誤差偏大。 筆者在充分考慮前人研究成果的基礎(chǔ)上，以抽水井實(shí)測(cè)數(shù)據(jù)為基礎(chǔ)，擬合計(jì)算降深曲線，并通過相對(duì)誤差分析確定計(jì)算精度分界點(diǎn)，在分界點(diǎn)以內(nèi)采用擬合曲線，在分界點(diǎn)以外采用理論曲線，二者組成分段函數(shù)，建立全區(qū)間符合實(shí)際工況的潛水滲流場(chǎng)。

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

3.1 工程概況

某灌區(qū)試驗(yàn)場(chǎng)抽水井與周邊16 口觀測(cè)井W1、W2、…、W16的距離分別為43、60、90、120、135、140、170、220、270、480、510、590、650、700、730、780 m。 布井場(chǎng)地范圍內(nèi)地質(zhì)條件差異不大，可按各向同性考慮。抽水井是一個(gè)淺層開放性地下水系統(tǒng)，可接收大氣降水、地表水、灌溉回歸水等垂直入滲補(bǔ)給，通過淺層潛水蒸發(fā)、人工開采側(cè)向徑流等輸出，地下水水力特性屬于潛水-微承壓水。

3.2 實(shí)測(cè)降深曲線數(shù)學(xué)模型

由于各井實(shí)際工況不同，因此判斷抽水井在長(zhǎng)時(shí)間定流量抽水后是否在井附近形成穩(wěn)定降落漏斗，目前尚無統(tǒng)一標(biāo)準(zhǔn)。 一般方法是選取一個(gè)抽水時(shí)段Δt，如果在該時(shí)段觀測(cè)井水位基本無變化，則認(rèn)為達(dá)到穩(wěn)定狀態(tài)。 本研究選取Δt ＝200 min，監(jiān)測(cè)抽水10、210、1 600、1 800 min 時(shí)各觀測(cè)井的水位降深（見表1）。

表1 不同抽水時(shí)間觀測(cè)井水位降深實(shí)測(cè)值

將上述4 個(gè)抽水時(shí)間的水位降深曲線繪制在同一張圖上進(jìn)行對(duì)比（見圖2），可以看出：在抽水開始t＝10 min 和t＝210 min 時(shí)同一觀測(cè)井的水位降深差異較大，在距抽水井較近的范圍（ 0 ≤r≤300 m）差異較為顯著；當(dāng)t＝1 600 min和t＝1 800 min 時(shí)，同一觀測(cè)井水位降深差異很小，表明此時(shí)抽水范圍內(nèi)已形成相對(duì)穩(wěn)定的滲流場(chǎng)，降落漏斗形狀和水位線保持穩(wěn)定。

由圖2 可知，t＝1 800 min 時(shí)，各觀測(cè)井實(shí)測(cè)降深為穩(wěn)定流數(shù)據(jù)。 對(duì)該組數(shù)據(jù)應(yīng)用Matlab 軟件進(jìn)行擬合，根據(jù)曲線形態(tài)和工程經(jīng)驗(yàn)，選取冪函數(shù)、指數(shù)函數(shù)、有理數(shù)逼近以及傅里葉函數(shù)等4 種函數(shù)進(jìn)行擬合，結(jié)果見圖3。

圖2 不同抽水時(shí)間各觀測(cè)井水位降深曲線

圖3 不同函數(shù)擬合效果對(duì)比

以冪函數(shù)擬合的降深曲線表達(dá)式為

以指數(shù)函數(shù)擬合的降深曲線表達(dá)式為

以有理數(shù)逼近擬合的降深曲線表達(dá)式為

以傅里葉函數(shù)擬合的降深曲線表達(dá)式為

其中：a0＝4.39×107，a1＝-4.39×107，b1＝1.593×104，w ＝-5.602×10-7。

采用誤差平方和（SSE）、決定系數(shù)（R2）、校正決定系數(shù)（Adjusted R2）、均方根誤差（RMSE）對(duì)上述4 種函數(shù)擬合程度的優(yōu)劣進(jìn)行判斷，結(jié)果見表2。

表2 擬合程度指標(biāo)

SSE越接近0，說明模型選擇越合理、擬合越好、數(shù)據(jù)預(yù)測(cè)越成功；R2越接近1，說明擬合效果越好；Adjusted R2判斷標(biāo)準(zhǔn)與決定系數(shù)相同；RMSE越接近0擬合效果越好。 根據(jù)上述判斷標(biāo)準(zhǔn)，由表2 可知，通過有理數(shù)逼近擬合的水位降深曲線在4 種函數(shù)中效果最好，基本可以反映穩(wěn)定抽水時(shí)真實(shí)的水位降深狀態(tài)，其函數(shù)表達(dá)式為

3.3 潛水滲流場(chǎng)及其母線方程

以抽水井底部中心為坐標(biāo)原點(diǎn)，井軸線為z軸，建立直角坐標(biāo)系。 以式（9）為母線繞井軸線旋轉(zhuǎn)，此時(shí)在潛水含水層中形成一個(gè)旋轉(zhuǎn)曲面，該旋轉(zhuǎn)曲面即為潛水含水層中的實(shí)際降落漏斗，其表達(dá)式為

式（10）即為研究工況下地下水潛水位分布的空間方程，該式所表達(dá)的滲流場(chǎng)基于工程實(shí)測(cè)資料生成，因此在靠近井壁區(qū)域，即當(dāng)r≤1.5H0時(shí)，同樣符合實(shí)際情況。

傳統(tǒng)的此類研究，由于井壁附近邊界條件復(fù)雜，其概念化的微分方程甚至不能寫出顯性表達(dá)式，更難于求出解析解，因此工程意義有限；而本研究采用實(shí)測(cè)資料擬合方案所推導(dǎo)出的方程，過程清晰、物理概念明確，降落漏斗空間形態(tài)是由真實(shí)降深曲線為母線旋轉(zhuǎn)而來，且具有數(shù)值解，相比微分方程（例如二維潛水運(yùn)動(dòng)的Boussinesq 方程），更方便實(shí)際應(yīng)用。

考慮到抽水井實(shí)際的坐標(biāo)系模型，將式（10）通過坐標(biāo)變換：x ＝rcosφ，y ＝rsinφ，z ＝z，以柱坐標(biāo)形式表示，其方程可轉(zhuǎn)化為式（9）。

4 結(jié)果與討論

4.1 計(jì)算曲線與理論曲線對(duì)比

將本研究方案計(jì)算生成的降深曲線與按式（4）理論降深曲線在同一坐標(biāo)系下對(duì)比，見圖4（圖4（b）中理論降深曲線僅為一般經(jīng)驗(yàn)形態(tài)，無具體坐標(biāo)數(shù)值）?？梢钥闯?，計(jì)算曲線和理論曲線在趨勢(shì)上一致，整體上呈現(xiàn)離抽水井越遠(yuǎn)降深越小的趨勢(shì)，當(dāng)r大于一定值時(shí)，降深幾乎不再變化；但在距抽水井較近位置時(shí)，二者形態(tài)有較大差異。 本研究方案降深曲線與實(shí)測(cè)數(shù)據(jù)吻合較好，而理論降深曲線在井附近變動(dòng)趨勢(shì)更為劇烈，存在明顯的“陡降”。 一般認(rèn)為，造成該現(xiàn)象的主要原因有二：一是由流速造成地下水入井水頭損失，二是地下水觸井后由近似水平運(yùn)動(dòng)轉(zhuǎn)為垂直運(yùn)動(dòng)造成能量損失。

圖4 計(jì)算生成的降深曲線與理論降深曲線對(duì)比

4.2 理論曲線分界點(diǎn)計(jì)算結(jié)果

利用式（4）計(jì)算得到的降深曲線見圖5。 可以看出，理論計(jì)算結(jié)果與實(shí)測(cè)結(jié)果偏差較大。 傳統(tǒng)理論認(rèn)為“當(dāng)r≥1.5H0時(shí)，理論降深曲線與實(shí)際降深曲線一致”，但該結(jié)論與本研究實(shí)際情況不符。 研究區(qū)初始潛水含水層厚度為23.25 m，按照傳統(tǒng)理論，當(dāng)與觀測(cè)井距離r≥1.5H0（為方便計(jì)算，不妨取40 m）時(shí)，理論曲線與實(shí)際曲線應(yīng)當(dāng)一致，而事實(shí)上直至r＝300 m 時(shí)理論值與實(shí)測(cè)值才較為接近（抽水井與最近的觀測(cè)井距離設(shè)計(jì)為43 m）。

圖5 理論曲線與實(shí)測(cè)值對(duì)比

為定量分析不同區(qū)間理論降深和實(shí)測(cè)值的差異，按照40＜r ＜120 m、120 m＜r ＜400 m、400 m＜r ＜800 m 將觀測(cè)井分為3 個(gè)區(qū)間進(jìn)行試算，通過考察不同區(qū)間理論降深和實(shí)測(cè)值的相對(duì)誤差，確定符合實(shí)際工況的計(jì)算分界點(diǎn)。 相對(duì)誤差計(jì)算公式為

Er＝（ss- sm）／sm（11）

式中：Er為相對(duì)誤差；ss為理論降深；sm為實(shí)測(cè)降深。

40 m＜r＜120 m 區(qū)間理論降深和實(shí)測(cè)值的相對(duì)誤差見表3。 可以看出，該區(qū)間4 組數(shù)據(jù)平均相對(duì)誤差為63.75%，最大相對(duì)誤差為65%；應(yīng)用同樣方法可得，120 m＜r ＜400 m 區(qū)間4 組數(shù)據(jù)平均相對(duì)誤差為17%，最大相對(duì)誤差為21%。 上述說明理論曲線在本研究40 m＜r ＜120 m 和120 m＜r ＜400 m 這兩個(gè)區(qū)間都不適用，盡管該距離已完全滿足r≥1.5H0。 當(dāng)觀測(cè)井處于400 m＜r ＜800 m 這一區(qū)間時(shí)，平均相對(duì)誤差為7%、最大相對(duì)誤差為9%，這表明直到觀測(cè)半徑r＞400 m 時(shí)，理論曲線精度才可滿足工程需求，因此Dupuit 公式在實(shí)際應(yīng)用中需十分慎重。

表3 40 m＜r＜120 m 區(qū)間理論降深和實(shí)測(cè)值的相對(duì)誤差

4.3 計(jì)算曲線-理論曲線組合方程

在實(shí)際工程中一般無法大量設(shè)置觀測(cè)井，因此判斷精度分界點(diǎn)的位置顯得尤為重要。 由上述理論曲線計(jì)算結(jié)果可知，r＝400 m 可作為本研究計(jì)算誤差分界點(diǎn)，在該范圍之外潛水降深理論值與實(shí)測(cè)值一致，可不設(shè)置觀測(cè)井；在該范圍之內(nèi)理論值與實(shí)測(cè)值誤差較大，需設(shè)置觀測(cè)井，不可用理論公式直接計(jì)算。 顯然對(duì)于具體工程，采用本文的計(jì)算分界點(diǎn)方案，比直接以r ＝1.5H0作為分界點(diǎn)更為精確且更具有實(shí)際意義。

將式（4）與式（9）聯(lián)立，在實(shí)際分界點(diǎn)之內(nèi)采用本研究構(gòu)建的有理數(shù)逼近方程，在實(shí)際分界點(diǎn)之外采用Dupuit 方程推導(dǎo)出的理論降深公式，這樣即可建立適用于本研究對(duì)象的自井附近至無限遠(yuǎn)處全區(qū)間范圍降深分段函數(shù)：

4.4 分段函數(shù)工程驗(yàn)證

將觀測(cè)井與抽水井的距離及其他相關(guān)數(shù)據(jù)代入式（12）進(jìn)行計(jì)算，得出分段降深曲線（見圖6）。 以該曲線為母線繞井軸旋轉(zhuǎn)所形成的降落漏斗即在本研究條件下的滲流場(chǎng)。 為驗(yàn)證式（12）在實(shí)際工程中的適用性，停止抽水一段時(shí)間，待試驗(yàn)場(chǎng)地潛水分布回歸穩(wěn)定狀態(tài)后重新進(jìn)行獨(dú)立抽水試驗(yàn)。 重新抽水后，取t＝2 000 min，對(duì)比各觀測(cè)井實(shí)測(cè)降深與計(jì)算降深。

圖6 分段降深曲線

以理論分界點(diǎn)r＝1.5H0（取40 m）和計(jì)算分界點(diǎn)r＝400 m 將試驗(yàn)場(chǎng)地分為r＜40 m、40 m＜r＜400 m、r＞400 m 3 個(gè)區(qū)域，將各區(qū)間降深實(shí)測(cè)值、分段函數(shù)計(jì)算值、Dupuit 理論曲線值三者對(duì)比，以判斷計(jì)算曲線和理論曲線在實(shí)際應(yīng)用中的差異。

在r＜40 m 區(qū)間，選擇井壁附近5、20、30 m 處3 個(gè)觀測(cè)孔，分別對(duì)比實(shí)測(cè)值、計(jì)算值、理論值，見表4。 可以看出，本文所構(gòu)建的分段函數(shù)計(jì)算值皆優(yōu)于Dupuit理論曲線計(jì)算值，尤其在靠近井壁區(qū)域。 當(dāng)觀測(cè)井距抽水井為5 m 時(shí)，Dupuit 理論曲線值與實(shí)測(cè)值相差2.67 m，相對(duì)誤差為56%，理論值與實(shí)測(cè)值偏差較大，Dupuit 公式顯然無法應(yīng)用于實(shí)際工程；相同位置處，本文構(gòu)建的函數(shù)曲線計(jì)算結(jié)果與實(shí)測(cè)值絕對(duì)誤差僅為0.16 m，相對(duì)誤差僅為3%，其計(jì)算精度遠(yuǎn)遠(yuǎn)高于理論曲線的。

表4 井壁附近觀測(cè)降深與計(jì)算降深對(duì)比

在40 m＜r＜400 m 區(qū)間，選取6 個(gè)觀測(cè)點(diǎn)，進(jìn)行理論值、計(jì)算值與實(shí)測(cè)值對(duì)比，見圖7。 可以看出，實(shí)測(cè)降深與理論值在該區(qū)間仍有較大偏差，絕對(duì)誤差最大值出現(xiàn)在r＝135 m 處，其值為0.31 m，最大相對(duì)誤差為15%；實(shí)測(cè)值與計(jì)算值趨勢(shì)則較為吻合，絕對(duì)誤差最大值出現(xiàn)在r＝265 m 處，其值為0.15 m，最大相對(duì)誤差為8%。

圖7 較遠(yuǎn)處理論值、計(jì)算值與實(shí)測(cè)值對(duì)比

進(jìn)一步，選用建模序列對(duì)比中常用的總量偏差BLAS、均方根納什效率系數(shù)SNSE兩個(gè)校驗(yàn)指標(biāo)進(jìn)行判斷，其中BLAS越小則效果越好，SNSE越接近1 效果越好，其數(shù)學(xué)意義見參考文獻(xiàn)［20］，二者計(jì)算公式分別為

式中：ssn為降深計(jì)算值；son為降深觀測(cè)值；N為觀測(cè)數(shù)據(jù)的數(shù)量。

理論值與實(shí)測(cè)值對(duì)比檢驗(yàn)結(jié)果為BLAS＝0.17、SNSE＝0.54，表明理論值不能滿足工程精度［20］要求；實(shí)測(cè)值與計(jì)算值對(duì)比檢驗(yàn)結(jié)果為BLAS＝0.08、SNSE＝0.86，表明計(jì)算值可以滿足工程精度要求。 本組對(duì)比結(jié)果再次表明，Dupuit 理論曲線即便在較遠(yuǎn)處與實(shí)測(cè)值仍有一定誤差，直接使用Dupuit 公式需謹(jǐn)慎，這一結(jié)論與大量前人研究結(jié)果相吻合［21-24］。

5 結(jié) 論

（1）冪函數(shù)、指數(shù)函數(shù)、有理數(shù)逼近、傅立葉函數(shù)等4 種擬合函數(shù)計(jì)算結(jié)果表明，對(duì)于本研究潛水井定流量抽水，有理數(shù)逼近擬合效果最好，擬合方程的決定系數(shù)R2值接近0.99。

（2）傳統(tǒng)方案理論曲線以r＝1.5H0作為計(jì)算精度分界點(diǎn)在實(shí)際工程中并不適用。 以本研究為例，應(yīng)用理論公式所計(jì)算的降深，不但在井壁附近與實(shí)際降深差別較大，在相當(dāng)大的范圍內(nèi)（r＜400 m）都存在較大誤差，該范圍遠(yuǎn)遠(yuǎn)大于本研究對(duì)應(yīng)的1.5H0≈40 m，因此應(yīng)用Dupuit 公式時(shí)應(yīng)慎重。

（3）以傳統(tǒng)精度分界點(diǎn)r＜40 m 為依據(jù)計(jì)算時(shí)，理論值與實(shí)測(cè)值最大誤差為62%，Dupuit 曲線在此區(qū)間不適用于本研究實(shí)例；計(jì)算值與實(shí)測(cè)值最大誤差為3%，計(jì)算曲線在此區(qū)間適用于本研究實(shí)例。

（4）40 m＜r＜400 m 區(qū)間降深理論值與實(shí)測(cè)值最大誤差為15%，在原本理論曲線認(rèn)為計(jì)算精度較高的區(qū)間，實(shí)際工程應(yīng)用仍有較大誤差，表明Dupuit 曲線在這一區(qū)間仍然適用性不強(qiáng)；計(jì)算值與實(shí)測(cè)值最大誤差為8%，表明分段函數(shù)曲線在該區(qū)間可用于本研究實(shí)例。

（5）當(dāng)r＞400 m 時(shí)，分段函數(shù)為原理論曲線表達(dá)式，理論值與實(shí)測(cè)值平均誤差為7%，最大誤差為9%，表明該區(qū)間可以采用Dupuit 公式。 因此，在無觀測(cè)條件時(shí)，要想采用Dupuit 公式計(jì)算潛水降深，為安全起見，需選擇相當(dāng)大的計(jì)算半徑。