基于分段函數的潛水井降深曲線在工程實踐中的應用

朱悅璐，熊俊飛，鄧 煒

（南昌工程學院 水利與生態工程學院，江西 南昌 330099）

1 引 言

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

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

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

2 傳統Dupuit 公式及其不足

潛水井向井運動是潛水含水層中地下水運動的一個特例，因此其動力學特性滿足潛水含水層中地下水運動的基本微分方程及假設條件，即Dupuit 假設：對于垂直的二維平面上任意一點（x，z），當潛水流垂向坡角θ很小（潛水位變化很小）時，可以忽略垂向分速度vz而僅考慮水平分速度vx（見圖1）。

圖1 Dupuit 假設

潛水井地下水向井運動的Dupuit 公式為

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

式（1）為直角坐標系下的標準微分方程，為方便與抽水井實際工況匹配，常用的方法是將式（1）轉化為柱坐標：

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

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

式中：s為水位降深。

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

3 實例應用

3.1 工程概況

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

3.2 實測降深曲線數學模型

由于各井實際工況不同，因此判斷抽水井在長時間定流量抽水后是否在井附近形成穩定降落漏斗，目前尚無統一標準。 一般方法是選取一個抽水時段Δt，如果在該時段觀測井水位基本無變化，則認為達到穩定狀態。 本研究選取Δt ＝200 min，監測抽水10、210、1 600、1 800 min 時各觀測井的水位降深（見表1）。

表1 不同抽水時間觀測井水位降深實測值

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

由圖2 可知，t＝1 800 min 時，各觀測井實測降深為穩定流數據。 對該組數據應用Matlab 軟件進行擬合，根據曲線形態和工程經驗，選取冪函數、指數函數、有理數逼近以及傅里葉函數等4 種函數進行擬合，結果見圖3。

圖2 不同抽水時間各觀測井水位降深曲線

圖3 不同函數擬合效果對比

以冪函數擬合的降深曲線表達式為

以指數函數擬合的降深曲線表達式為

以有理數逼近擬合的降深曲線表達式為

以傅里葉函數擬合的降深曲線表達式為

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

采用誤差平方和（SSE）、決定系數（R2）、校正決定系數（Adjusted R2）、均方根誤差（RMSE）對上述4 種函數擬合程度的優劣進行判斷，結果見表2。

表2 擬合程度指標

SSE越接近0，說明模型選擇越合理、擬合越好、數據預測越成功；R2越接近1，說明擬合效果越好；Adjusted R2判斷標準與決定系數相同；RMSE越接近0擬合效果越好。 根據上述判斷標準，由表2 可知，通過有理數逼近擬合的水位降深曲線在4 種函數中效果最好，基本可以反映穩定抽水時真實的水位降深狀態，其函數表達式為

3.3 潛水滲流場及其母線方程

以抽水井底部中心為坐標原點，井軸線為z軸，建立直角坐標系。 以式（9）為母線繞井軸線旋轉，此時在潛水含水層中形成一個旋轉曲面，該旋轉曲面即為潛水含水層中的實際降落漏斗，其表達式為

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

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

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

4 結果與討論

4.1 計算曲線與理論曲線對比

將本研究方案計算生成的降深曲線與按式（4）理論降深曲線在同一坐標系下對比，見圖4（圖4（b）中理論降深曲線僅為一般經驗形態，無具體坐標數值）。可以看出，計算曲線和理論曲線在趨勢上一致，整體上呈現離抽水井越遠降深越小的趨勢，當r大于一定值時，降深幾乎不再變化；但在距抽水井較近位置時，二者形態有較大差異。 本研究方案降深曲線與實測數據吻合較好，而理論降深曲線在井附近變動趨勢更為劇烈，存在明顯的“陡降”。 一般認為，造成該現象的主要原因有二：一是由流速造成地下水入井水頭損失，二是地下水觸井后由近似水平運動轉為垂直運動造成能量損失。

圖4 計算生成的降深曲線與理論降深曲線對比

4.2 理論曲線分界點計算結果

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

圖5 理論曲線與實測值對比

為定量分析不同區間理論降深和實測值的差異，按照40＜r ＜120 m、120 m＜r ＜400 m、400 m＜r ＜800 m 將觀測井分為3 個區間進行試算，通過考察不同區間理論降深和實測值的相對誤差，確定符合實際工況的計算分界點。 相對誤差計算公式為

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

式中：Er為相對誤差；ss為理論降深；sm為實測降深。

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

表3 40 m＜r＜120 m 區間理論降深和實測值的相對誤差

4.3 計算曲線-理論曲線組合方程

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

將式（4）與式（9）聯立，在實際分界點之內采用本研究構建的有理數逼近方程，在實際分界點之外采用Dupuit 方程推導出的理論降深公式，這樣即可建立適用于本研究對象的自井附近至無限遠處全區間范圍降深分段函數：

4.4 分段函數工程驗證

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

圖6 分段降深曲線

以理論分界點r＝1.5H0（取40 m）和計算分界點r＝400 m 將試驗場地分為r＜40 m、40 m＜r＜400 m、r＞400 m 3 個區域，將各區間降深實測值、分段函數計算值、Dupuit 理論曲線值三者對比，以判斷計算曲線和理論曲線在實際應用中的差異。

在r＜40 m 區間，選擇井壁附近5、20、30 m 處3 個觀測孔，分別對比實測值、計算值、理論值，見表4。 可以看出，本文所構建的分段函數計算值皆優于Dupuit理論曲線計算值，尤其在靠近井壁區域。 當觀測井距抽水井為5 m 時，Dupuit 理論曲線值與實測值相差2.67 m，相對誤差為56%，理論值與實測值偏差較大，Dupuit 公式顯然無法應用于實際工程；相同位置處，本文構建的函數曲線計算結果與實測值絕對誤差僅為0.16 m，相對誤差僅為3%，其計算精度遠遠高于理論曲線的。

表4 井壁附近觀測降深與計算降深對比

在40 m＜r＜400 m 區間，選取6 個觀測點，進行理論值、計算值與實測值對比，見圖7。 可以看出，實測降深與理論值在該區間仍有較大偏差，絕對誤差最大值出現在r＝135 m 處，其值為0.31 m，最大相對誤差為15%；實測值與計算值趨勢則較為吻合，絕對誤差最大值出現在r＝265 m 處，其值為0.15 m，最大相對誤差為8%。

圖7 較遠處理論值、計算值與實測值對比

進一步，選用建模序列對比中常用的總量偏差BLAS、均方根納什效率系數SNSE兩個校驗指標進行判斷，其中BLAS越小則效果越好，SNSE越接近1 效果越好，其數學意義見參考文獻［20］，二者計算公式分別為

式中：ssn為降深計算值；son為降深觀測值；N為觀測數據的數量。

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

5 結 論

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

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

（3）以傳統精度分界點r＜40 m 為依據計算時，理論值與實測值最大誤差為62%，Dupuit 曲線在此區間不適用于本研究實例；計算值與實測值最大誤差為3%，計算曲線在此區間適用于本研究實例。

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

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