地下水流并行有限層方法及同倫反演研究

王少偉， 徐 進(jìn)， 楊偉濤

(煙臺(tái)大學(xué) 土木工程學(xué)院，煙臺(tái) 264005)

1 引 言

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

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

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

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

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

2 地下水?dāng)?shù)學(xué)模型

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

(1)

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

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

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

(2)

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

3 有限層求解格式及其并行特點(diǎn)

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

3.1 有限層求解格式

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

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

(3)

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

3.2 求解格式的特點(diǎn)

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

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

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

4 求解格式并行化方法

結(jié)合地下水流有限層求解格式的上述特點(diǎn)，利用MATLAB并行計(jì)算庫(kù)，對(duì)求解格式(3)進(jìn)行并行化處理。

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

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

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

5 地下水流降深反演計(jì)算

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

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

(4)

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

(5)

(6)

圖2 降深正演并行計(jì)算程序流程

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

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

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

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

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

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

圖3 降深反演計(jì)算并行程序流程

利用本文計(jì)算結(jié)果，圖5分別給出了第一承壓含水層和開(kāi)采含水層頂板處的降深隨距離的分布曲線(xiàn)。由于問(wèn)題的復(fù)雜性，不存在方便對(duì)比的解析解，本文采用有限差分計(jì)算結(jié)果進(jìn)行了對(duì)比?？梢钥闯觯疚牟⑿薪獯鹋c有限差分結(jié)果有較好的擬合度，計(jì)算結(jié)果很好地反映了多層結(jié)構(gòu)含水層的地下水三維越流特征。

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

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

圖4 承壓含水層完整井引起的水位降深曲線(xiàn)

圖5 三層結(jié)構(gòu)承壓含水層水位降深

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

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

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

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

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

6.5 并行計(jì)算效率分析

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

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

表1 選取不同初值情況下的反演計(jì)算結(jié)果

表2 計(jì)算機(jī)配置

圖8 不同函數(shù)級(jí)數(shù)項(xiàng)下的并行計(jì)算效率對(duì)比

7 結(jié) 論

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

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

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