王少偉， 徐 進， 楊偉濤

(煙臺大學 土木工程學院，煙臺 264005)

1 引 言

地下水是寶貴的自然資源，也是地質(zhì)環(huán)境的基本要素[1]。過量開采地下水會造成區(qū)域性降落漏斗，導致地面沉降和海水入侵等環(huán)境地質(zhì)問題。實現(xiàn)地下水流問題的高效求解對上述問題的預測和防治具有重要意義。

地下水流問題的常用求解方法包括解析法和數(shù)值法等。解析法[2,3]是利用積分變換等分析手段直接求解地下水流控制方程，得到問題的精確顯式解。Craig[4]提出了二維非均質(zhì)穩(wěn)定流問題的一般解，解的形式包括拉普拉斯方程的復多項式以及由多項式系數(shù)確定的附加非全純項，可以通過選擇多項式系數(shù)滿足不同邊界的流動條件。Zha等[5]定義了無限邊界、等效均質(zhì)含水層中抽水引起的二維和三維水流的近似穩(wěn)態(tài)解及穩(wěn)態(tài)條件，討論了非均質(zhì)含水層的穩(wěn)態(tài)解條件。然而，對于考慮復雜非均質(zhì)性和各向異性情形的三維地下水流問題，目前解析法仍難以應用。

數(shù)值方法(有限差分法和有限單元法等)是目前求解地下水流模型的主流方法，可以反映復雜水文地質(zhì)條件下的水流形態(tài)，得到廣泛研究和應用。潘樹來等[6]提出了一種無需對自由面作近似處理的精確算法，為實現(xiàn)有自由面非穩(wěn)定滲流問題的精細分析提供技術保證。Xie等[7]提出了一種新的有限體積多尺度有限元模型，該方法可以簡化仿真過程，有效模擬多孔介質(zhì)中的地下水流動。但是，在對區(qū)域性地下水流三維非穩(wěn)定流的精細化分析時，目前的數(shù)值方法仍然存在計算工作量大和計算時間長等問題[8,9]。除了傳統(tǒng)解析法和數(shù)值法，有限層方法[10]是一種混合算法，將解析手段和數(shù)值方法相結合，使得三維問題簡化為一維問題，可以實現(xiàn)地下水流模型的高效求解[11-14]。

值得注意的是，隨著地下水流模型的不斷完善，理論計算與數(shù)值模擬中涉及到的參數(shù)需要事先確定。如何通過實驗或現(xiàn)場資料有效獲取各類水文地質(zhì)參數(shù)非常重要[15-19]。然而，反演分析中需要多次調(diào)用正演程序，這加劇了數(shù)值方法本身存在的計算困難。另外，反演計算本身存在初值選擇問題，如果初值選取不合適，太過偏離真值，往往會出現(xiàn)收斂速度慢、局部收斂甚至發(fā)散等問題。本文主要結合同倫算法對參數(shù)反演方法進行改進，Keller等[20]首次將同倫法用于反演方法中，使得同倫算法的研究得到了迅速發(fā)展。文獻[21,22]利用同倫方法對不同模型參數(shù)進行反演計算，結果表明，同倫參數(shù)反演方法具有良好的穩(wěn)定性、收斂性和較快的計算速度。

本文針對地下水流有限層求解格式的解耦性，實現(xiàn)并行化，進一步提高其處理地下水三維流問題的計算效率。同時，結合大范圍收斂的非線性同倫方法[23,24]進行反演計算，避免初值選取問題。在此基礎上，利用MATLAB編制并行化計算程序和同倫反演程序，并通過算例對比驗證方法和程序的有效性和高效性。

2 地下水數(shù)學模型

根據(jù)質(zhì)量守恒定律和達西定律，在直角坐標系下，地下水三維、非穩(wěn)定流的數(shù)學模型用降深可表示為

(1)

式中S(x,y,z,t)為降深(L),Kx，Ky和Kz分別為x，y和z方向的滲透系數(shù)(LT-1)，q(x,y,z,t)為匯源項(T-1)，Ss為貯水率(L-1)。Sy為給水度(無量綱)，反映了潛水面單位降深重力排水的補給強度，當考慮承壓水流時，可以令Sy=0。

式(1)是地下水三維流的正演數(shù)學模型。對于參數(shù)反演問題，主要依據(jù)實測水位或者降深結果來反求含水層的未知水文地質(zhì)參數(shù)，使得

S(xi,yi,zi,tj)=Sc(xi,yi,zi,tj)

(2)

式中S(xi,yi,zi,tj)為計算降深值，Sc(xi,yi,zi,tj) 為i位置在j時刻的實測降深值。在求解地下水流正演數(shù)學模型(1)的基礎上，結合一定的優(yōu)化算法可以實現(xiàn)對反問題(2)的分析。

3 有限層求解格式及其并行特點

根據(jù)含水層系統(tǒng)的層狀非均質(zhì)特點，將含水層系統(tǒng)沿z方向離散成L個不同厚度的層元，每一層可以具有不同的水文地質(zhì)參數(shù)以考慮這種層狀非均質(zhì)。

3.1 有限層求解格式

降深在平面x和y方向分別選用M和N項正交完備的解析函數(shù)系表示，z方向則采用常規(guī)有限元線性形函數(shù)來離散，最終得到解析數(shù)值混合形式的降深試探函數(shù)[12]。

將降深試探函數(shù)代入正演數(shù)學模型(1)，基于伽遼金原理或者泛函極值原理，可以得到正演數(shù)學模型的求解格式，用矩陣形式表示為[12]

(3)

式中[G]和[B]分別為整體滲透矩陣以及整體貯水矩陣，{Q}為水量向量，{Φ}為總待求系數(shù)向量。

3.2 求解格式的特點

(1) 求解格式(3)的[G]和[B]等整體系數(shù)矩陣和向量具有顯式的解析表達式，無需數(shù)值積分。與傳統(tǒng)的數(shù)值方法相比，極大減輕了計算工作量，節(jié)省了運算時間。矩陣元素具體解析式詳見文獻[12]。

(2) 利用所引入解析函數(shù)的正交性，可以實現(xiàn)求解格式(3)的整體解耦，將其拆分為M×N個獨立的L+1階子線性系統(tǒng)，拆分后的系數(shù)矩陣[G]和[B]如圖1所示。

(3) 整體矩陣拆解后的子矩陣互不耦合，相互獨立，可以分別求解，具有天然的并行性。利用這一特點進行并行計算時，各進程之間無需數(shù)據(jù)交換，不同進程單獨處理，節(jié)省運算時間。

4 求解格式并行化方法

結合地下水流有限層求解格式的上述特點，利用MATLAB并行計算庫，對求解格式(3)進行并行化處理。

針對求解格式關于級數(shù)項M和N的解耦性，采用SPMD(single program multiple data)的方式進行并行化處理[25]。通過SPMD和End命令定義一個SPMD并行結構代碼段，各進程使用相同主程序，具有不同數(shù)據(jù)。程序中使用parpool(K)或parpool指令打開MATLAB并行池，打開K個線程(worker)或默認線程。根據(jù)worker的數(shù)目，利用labindex命令將級數(shù)項M和N分成M1，M2，…，Mi和N1，N2，…，Nj組數(shù)據(jù)，分別計算降深值(i和j通過worker數(shù)目決定)。

根據(jù)求解格式特點，各進程單獨計算，無需通訊環(huán)節(jié)，待各進程計算結束，將結果返還給直接累加，可以得到最終降深值，程序流程如圖2所示。

圖1 整體系數(shù)矩陣的解耦性

5 地下水流降深反演計算

在并行正演方法的基礎上，使用同倫方法進行反演計算。同倫方法是一種非線性優(yōu)化算法，主要采用逐步逼近的計算方式對初值逐次優(yōu)化，具有大范圍收斂的特點。

同倫方法的基本思想是將一個簡單問題與所要求的復雜問題放到一個一般問題中，首先求出簡單問題的解，然后逐次連續(xù)過渡到復雜目標問題的解[23,24]。

(4)

根據(jù)同倫方法的思想，構造如下同倫函數(shù)

(5)

(6)

圖2 降深正演并行計算程序流程

6 數(shù)值算例及驗證

根據(jù)上述并行策略和反演算法，分別編寫MATLAB計算程序以實現(xiàn)地下水流正反演問題的分析。為了驗證本文方法和計算程序的正確性，選取了承壓水完整井、多層結構含水層系統(tǒng)非穩(wěn)定流問題進行了分析，并與已有解答進行對比。最后，通過數(shù)值算例探討并行程序的計算效率和反演程序的收斂性。

6.1 單層承壓含水層的地下水流分析

利用本文計算程序?qū)?jīng)典的單層承壓含水層完整井引起的地下水非穩(wěn)定流問題進行計算。在算例中,含水層厚度c=100 m，貯水率Ss=2×10-5/m，滲透系數(shù)Kx=Ky=Kz=5 m/d，抽水量Q=2000 m3/d。利用計算結果，圖4給出了距離井中心x=y=10 m處不同時刻的完整井抽水降深曲線，可以看出本文解與文獻[26]的解析解曲線非常吻合，證明了并行化正演程序的正確性和計算精度。

6.2 多層含水層系統(tǒng)的地下水三維流分析

為了驗證并行正演程序處理復雜地下水流問題的能力，選取三層結構含水層系統(tǒng)三維流問題進行了分析，如圖5所示，整個含水層由厚為20 m的第一承壓含水層、中間厚為10 m的弱透水層和厚為30 m的第二承壓含水層組成，各層的水文地質(zhì)參數(shù)見文獻[12]。

圖3 降深反演計算并行程序流程

利用本文計算結果，圖5分別給出了第一承壓含水層和開采含水層頂板處的降深隨距離的分布曲線。由于問題的復雜性，不存在方便對比的解析解，本文采用有限差分計算結果進行了對比?？梢钥闯?，本文并行解答與有限差分結果有較好的擬合度，計算結果很好地反映了多層結構含水層的地下水三維越流特征。

6.3 單層含水層地下水流模型參數(shù)反演分析

為了驗證本文地下水流同倫反演程序的正確性，設置一單層承壓含水層抽水降深算例，利用Theis解，代入算例6.1計算參數(shù)生成水位降深時間序列，記錄時間段為t=0 d～2 d，時間間隔為Δt= 0.1 d，測點坐標為x=y=10 m。將該降深時間序列作為反演計算中的實測值，假定K為未知參數(shù)進行反演。

圖4 承壓含水層完整井引起的水位降深曲線

圖5 三層結構承壓含水層水位降深

6.4 多層結構含水層地下水流模型參數(shù)反演分析

利用算例6.2的計算參數(shù)，記第一和第二兩含水層的滲透系數(shù)分別為K1和K2，通過正演程序計算生成含水層頂和底板距井點10 m處的降深時間序列 (記錄時間段為t=0 d～2 d，時間間隔為 Δt=0.1 d)，以此作為資料代入同倫程序來反演滲透系數(shù)K1和K2。

將初始值選為K0=(10,10) m/d，同倫曲線如圖7所示，其中

圖6 參數(shù)反演同倫曲線

圖7 參數(shù)反演同倫曲線

6.5 并行計算效率分析

利用有限層求解格式存在的解耦性進行并行化處理，從而達到提高計算效率的目的。區(qū)別于有限元等純數(shù)值方法，所采用的解析函數(shù)項M和N對有限層方法計算結果的精度有重要影響。為了說明并行求解效率，針對不同解析級數(shù)項數(shù)，對比反演計算中不同線程的運算時間，計算機硬件配置列入表2。

計算中將計算時間固定，生成一組距井點不同位置的降深值作為反演實測值。反演計算中采用四線程并行程序，觀察在不同級數(shù)項數(shù)下的串行和并行運算時間。圖8將計算時間固定為1 d，由曲線可以看出，當級數(shù)項數(shù)較小時，并行運算時間要高于串行計算；計算循環(huán)量大于60時，并行優(yōu)化的高效性開始體現(xiàn)；在100項處，計算效率提升了近 1倍。這是由于并行反演計算程序較為復雜，程序間存在數(shù)據(jù)傳輸問題，占用一定運行時間。

表1 選取不同初值情況下的反演計算結果

表2 計算機配置

圖8 不同函數(shù)級數(shù)項下的并行計算效率對比

7 結 論

(1) 基于地下水流三維有限層求解格式，建立了該方法的并行化算法，通過MATLAB編制相應程序?qū)崿F(xiàn)了地下水流有限層方法的并行化，利用解析解以及有限差分解驗證了并行方法與計算程序的正確性。

(2) 利用上述并行算法作為正演程序，基于非線性同倫原理，進一步開發(fā)了地下水參數(shù)反演程序，利用數(shù)值算例驗證了程序的正確性；同時，研究表明同倫反演算法不依賴于參數(shù)的初值選取，具有大范圍收斂的特點，適用于地下水流反演計算問題的數(shù)值求解。

(3) 通過數(shù)值算例說明了并行化程序的計算效率。利用本文并行算法的解耦性，在并行計算中可以節(jié)省計算成本，便于編程實現(xiàn)。