邢利英,張國珍

(1.蘭州交通大學(xué)環(huán)境與市政工程學(xué)院,甘肅蘭州 730070;2.南陽師范學(xué)院土木與建筑工程學(xué)院,河南南陽 473061)

基于改進(jìn)的共軛梯度算法重構(gòu)地下水污染源項

邢利英1,2,張國珍1

(1.蘭州交通大學(xué)環(huán)境與市政工程學(xué)院,甘肅蘭州 730070;2.南陽師范學(xué)院土木與建筑工程學(xué)院,河南南陽 473061)

針對地下水污染源項的識別問題,為重構(gòu)非恒定地下水污染源項,采用有限差分格式離散方程,同時附加終端觀測值,利用改進(jìn)的共軛梯度算法反演地下水的污染源項。結(jié)果表明:對于連續(xù)以及不連續(xù)不可導(dǎo)的污染源項均能夠被重構(gòu);此外,考慮測量誤差的影響,數(shù)值結(jié)果顯示在有測量誤差的情況下,該算法也能較好地反演污染源項,充分說明改進(jìn)的共軛梯度算法反演地下水污染源項是有效的。

地下水;污染源項;重構(gòu);反問題;共軛梯度法

隨著現(xiàn)代化的快速發(fā)展,大量工農(nóng)業(yè)污染物的排放導(dǎo)致地下水污染嚴(yán)重。地下水中污染物的遷移過程作用時間長,影響范圍大,在遷移過程中還涉及水文地理條件、物理反應(yīng)、化學(xué)反應(yīng)等諸多因素,使得污染物遷移過程變得極為復(fù)雜,并表現(xiàn)出非線性行為。由于地下水系統(tǒng)本身的隱藏性和復(fù)雜性,地下水污染并沒有引起業(yè)界足夠的重視。最近的一份調(diào)查報告顯示,作為主要淡水資源之一的地下水已經(jīng)被氮、碳?xì)浠衔铩⒂卸居泻Φ奈⒘吭貒?yán)重污染,并且呈現(xiàn)出由點(diǎn)到面、由城鎮(zhèn)向農(nóng)村發(fā)展的趨勢,所以此類問題變得越來越棘手[1-4]。

事實上,在地下水污染物的實際測量中,有很多信息難以準(zhǔn)確測量,甚至無法測量,然而可以通過一些已知信息來確定未知的信息,所以在地下水的反問題研究中包含參數(shù)識別、邊界確定、源項識別等幾類問題,而污染源項識別問題在地下水污染中具有特殊的意義。目前已有大量的文獻(xiàn)研究地下水污染源的識別問題,Li等[5]應(yīng)用全域多元二次曲面方法反演污染源項;Hazart等[6]設(shè)計基于馬爾科夫鏈的貝葉斯參數(shù)識別法重構(gòu)點(diǎn)源的位置、釋放時間和源強(qiáng)的大小;王建平等[7]基于馬爾科夫鏈的貝葉斯方法進(jìn)行了水質(zhì)參數(shù)的不確定性研究;李子等[8]采用時空全域徑向基函數(shù)配點(diǎn)法構(gòu)造懲罰函數(shù),通過搜尋形狀因子和比例因子來確定污染源項的位置及強(qiáng)度;Zeng等[9]采用稀疏網(wǎng)格差值的貝葉斯方法確定污染源項。關(guān)于污染源項反演的其他方法請參考文獻(xiàn)[10-14]。然而,之前的研究中很少采用改進(jìn)的共軛梯度法反演污染源項。

共軛梯度法在數(shù)學(xué)理論分析中應(yīng)用較多,Deng等[15]采用共軛梯度方法重構(gòu)了傳熱方程的輻射項系數(shù),數(shù)值結(jié)果表明該方法有效地反演了輻射項系數(shù);Yang等[16]應(yīng)用共軛梯度方法重構(gòu)傳熱方程的源項。本文設(shè)計一種改進(jìn)的共軛梯度方法反演地下水污染源項。

1 地下水污染控制方程

地下水污染溶質(zhì)的輸運(yùn)和遷移過程是一個相當(dāng)復(fù)雜的物理化學(xué)過程,主要包括對流、擴(kuò)散、吸附、交換和其他的物理化學(xué)過程。簡單起見,這里只考慮沿主流方向,結(jié)合初邊界條件,污染物遷移模型可用一維非穩(wěn)定對流擴(kuò)散方程表示:

式中:c(x,t)為測點(diǎn)(x,t)的污染物濃度,mg/L; q(x)為污染源項,僅考慮為空間變量的函數(shù),mg/(L·s); aL為污染物沿x方向的擴(kuò)散度,m;υ為地下水實際的流速,m/s,通過水動力模型計算給出;λ為污染物的衰減系數(shù),s-1;ne為有效地孔隙率,無量綱; Q為計算區(qū)域;T為溶質(zhì)輸運(yùn)時間,s;L為主流方向長度,m;f1(t)和f2(t)為邊界條件;φ(x)為初始條件。

簡言之,如果方程組(1)中的參數(shù)aL、υ、λ、ne、q(x)以及初邊界條件均已知,則方程組(1)構(gòu)成一個正問題,可以通過實驗和數(shù)值模擬的方法預(yù)測污染物濃度c(x,t)。然而,大多數(shù)情況下,地下水的污染源項并不能準(zhǔn)確地獲得,因此可以通過已知的參數(shù)和初邊界條件推測未知的污染源項,方程組(1)則轉(zhuǎn)變?yōu)橐粋€反問題。

解決上述反問題,通常可以附加一組容易獲得的終端觀測值,一般應(yīng)用如下的形式:

式中,cT(x)為終端觀測值。

由此,式(1)和(2)就構(gòu)成了一個關(guān)于源項識別的反問題,可以應(yīng)用迭代的方法求解。

2 差分方程

假設(shè)計算區(qū)域為Q=[0,L]×[0,T],用有限差分法求解偏微分方程的源項識別問題,首先對計算區(qū)域Q進(jìn)行網(wǎng)格剖分。將區(qū)域Q劃分M×N等分,則空間步長和時間步長分別為測點(diǎn)(xj,tn)是網(wǎng)格節(jié)點(diǎn)。

式中:h,τ分別為空間步長和時間步長;M、N為兩個整數(shù)。則方程組(1)中的第一個式子可轉(zhuǎn)化如下簡單的形式:

式中:ct、cxx、cx分別為污染物濃度c關(guān)于時間的一階導(dǎo)數(shù)、空間的二階導(dǎo)數(shù)、空間的一階導(dǎo)數(shù)。四點(diǎn)隱格式的差分方程為

由于隱格式的差分方程是無條件穩(wěn)定的,因此可應(yīng)用任意步長求解線性方程(9)。下面將詳細(xì)地介紹改進(jìn)的共軛梯度算法。

3 改進(jìn)的共軛梯度算法

共軛梯度算法由Fletcher等[17]首次提出的,作為一種簡單有效的迭代技術(shù),廣泛應(yīng)用于地理學(xué)和傳熱學(xué)方程的源項及擴(kuò)散項系數(shù)反演。在大中型線性和非線性的無約束優(yōu)化問題中,共軛梯度算法尤其有效。共軛梯度算法的基本思路為結(jié)合共軛性和最速下降法,構(gòu)建一組共軛負(fù)梯度方向,沿著下降的方向?qū)ふ夷繕?biāo)函數(shù)的最小值,合適終止條件的共軛梯度法實際上是一類迭代正則化方法。鑒于共軛梯度的特點(diǎn),本文擬進(jìn)一步改進(jìn)共軛梯度算法重構(gòu)地下水污染源項。

設(shè)非線性優(yōu)化問題的控制方程可以表示為

其中D(F):X→Y為一有界線性算子,X和Y屬于Hilbert空間。對于x∈X,通常設(shè)F為一連續(xù)性的算子,所以上述非線性優(yōu)化問題可轉(zhuǎn)化為

式中:f(x):Rn→R是一連續(xù)可微的函數(shù);‖·‖2為歐幾里得模數(shù);δ為誤差范圍。

尤其是當(dāng)n很大時,共軛梯度相當(dāng)有效,式(12)則可轉(zhuǎn)化為

詳細(xì)的改進(jìn)共軛梯度算法步驟[15,20]如下:

步驟1:設(shè)k=0,選擇一個迭代初始值q0(x)= 0,x∈[0,L];

步驟2:求解初邊值問題(式(1)),得到其解ck,此時q(x)=qk(x);

步驟3:確定殘差rk=gδ-ck,求解下面的共軛方程,其解為uk(x,t);

步驟4:計算sk=uk+αk-1sk-1,其中

并得到其解為c:=dk(x,t),令

其中wk-1=uk-uk-1,則有qk+1=qk+βksk。

步驟6:增加迭代次數(shù)k,執(zhí)行循環(huán)到步驟2。任取γ>1,當(dāng)k∈Ν,第一次滿足‖gk-g‖L2(0,L)≤γδ時終止循環(huán)。

4 算 例

本文設(shè)計了3個算例來檢驗上述算法的有效性。a.試驗1。采用一假設(shè)的數(shù)值算例驗證改進(jìn)共軛梯度法的有效性,控制方程如下:

當(dāng)式(18)取第一組參數(shù):a(x)=1,φ(x)=sin (2πx),q(x)=2π2sin(2πx),時間步長和空間步長分別取1/100s和1/100 m,應(yīng)用改進(jìn)的共軛梯度算法反演污染源項的結(jié)果見圖1(a)。從圖1(a)中可以看出,當(dāng)?shù)螖?shù)K=50時,反演結(jié)果已經(jīng)很好地與真實解吻合。

當(dāng)取第二組參數(shù):φ(x)=sin(x),q(x)=|(x -1)(x-3)|,x∈[0,4],則反演結(jié)果見圖1(b)。從圖1(b)中可以看出,除了2個不可導(dǎo)點(diǎn),源項被較好地反演。以上2組試驗很好地說明了共軛梯度法能有效地反演控制方程源項,可以應(yīng)用到實際問題源項的重構(gòu)中。

b.試驗2。對淄博灃水南部地區(qū)的地下水硫酸污染源強(qiáng)度進(jìn)行數(shù)值反演,相關(guān)數(shù)據(jù)參考文獻(xiàn)[21]。根據(jù)計算區(qū)域的水文地質(zhì)條件,取x=4 000 m,T= 11 a,并取αL=1 m,u=1 m/d,ne=0.25,空間步長h=100 m,時間步長τ=0.5 a。通過線性擬合,初邊值條件表達(dá)式為

附加終端觀測值為擬合后的線性函數(shù)c(x,11)=0.102 6x+152。

圖1 試驗1的反演結(jié)果

應(yīng)用改進(jìn)的共軛梯度算法進(jìn)行反演,附加數(shù)據(jù)和反演結(jié)果的對比見圖2。從圖2中可以看出當(dāng)?shù)螖?shù)K=30時,污染源項已經(jīng)較好地被重構(gòu),充分說明該方法收斂速度較快,反演結(jié)果穩(wěn)定。采用cδT(x)=cδ(x,T)=c(x,T)[1+δrandom(x)]計算測量誤差,當(dāng)測量誤差δ分別為2%和5%時,污染源項的真值和反演結(jié)果的計算誤差分別為3.04%和2.83%,比參考文獻(xiàn)[21]的誤差4.18%要小。

c.試驗3。為了進(jìn)一步驗證改進(jìn)共軛梯度算法的有效性,將污染源項函數(shù)變化為一分段函數(shù)(式(20)),空間距離增至為6 000 m,其他模型參數(shù)與試驗2相同。不同迭代次數(shù)K和不同測量誤差δ的反演結(jié)果見圖3。

由于污染源項的不連續(xù)性,當(dāng)?shù)螖?shù)增大到100次時,污染源項才被較好地重構(gòu)。類似地,考慮測量誤差δ的影響時,盡管污染源項的不連續(xù)性,而且測量誤差也會產(chǎn)生影響,但是模擬結(jié)果也比較正常,充分證明改進(jìn)共軛梯度算法的穩(wěn)定性和較大的容錯性。

5 結(jié) 語

圖2 實際問題的反演結(jié)果與附加數(shù)據(jù)的比較

圖3 試驗3的反演結(jié)果

采用有限差分法離散地下水控制方程,應(yīng)用改進(jìn)的共軛梯度算法重構(gòu)地下水污染源項。數(shù)值模擬結(jié)果表明,對于性質(zhì)友好的源項以及不連續(xù)不可導(dǎo)的源項均能夠被反演。考慮測量誤差的影響,結(jié)果顯示在有測量誤差的情況下也能得到穩(wěn)定解,充分說明改進(jìn)的共軛梯度算法是一種有效的反演方法。污染源項如果是連續(xù)的函數(shù),迭代次數(shù)大約為50次,反演結(jié)果趨于穩(wěn)定;污染源項如果是不連續(xù)的函數(shù),迭代次數(shù)大約為100次時,反演結(jié)果趨于真實解,從而說明改進(jìn)的共軛梯度算法是一種快速穩(wěn)定的反演方法。

[1]ZHOU H Y,GOMEZ-HEMANDEZ J J,LI L P. Inverse methods in hydrogeology:evolution and recent trends[J].Advances in Water Resources, 2014,63:22-37.

[2]孫才志,陳相濤,陳雪姣,等.地下水污染風(fēng)險評價研究進(jìn)展[J].水利水電科技進(jìn)展,2015,35(5):152-161. (SUN Caizhi,CHEN Xiangtao,CHEN Xuejiao,et al. Recent advances in groundwater contamination risk assessment[J].Advances in Science and Technology of Water Resources,2015,35(5):152-161.(in Chinese))

[3]YEH W W G.Review:optimization methods for groundwater modeling and management[J]. Hydrogeology Journal,2015,23(6):1051-1065.

[4]朱嵩.基于貝葉斯推理的環(huán)境水力學(xué)反問題研究[D].杭州:浙江大學(xué),2008.

[5]LI Zi,MAO Xiaokang,ZHOU Kang.Global multiquadric collocation method for groundwater contaminant source identification[J].Environmental Modelling&Software,2011,26(12):1611-1621.

[6]HAZART A,GIOVANNELLI J F,DUBOST S,et al.Inverse transport problem of estimating point-like source using a Bayesian parametric method with MCMC[J].Signal Processing,2014,96:346-361.

[7]王建平,程聲通,賈海峰.基于MCMC法的水質(zhì)模型參數(shù)不確定性研究[J].環(huán)境科學(xué),2006,27(1):24-30. (WANG Jianping,CHENG Shengtong,JIA Haifeng. Markov chain monte carlo scheme for parameter uncertainty analysis in water quality model[J]. Environmental Science,2006,27(1):24-30.(in Chinese)).

[8]李子,毛獻(xiàn)忠,周康.污染物非恒定輸運(yùn)逆過程反演模型研究[J].水力發(fā)電學(xué)報,2013,32(6):115-121. (LI Zi,MAO Xianzhong,ZHOU Kang.Inversion model for the inverse process of unsteady pollutant transport[J].Journal of Hydroelectric Engineering, 2013,32(6):115-121.(in Chinese))

[9]ZENG Lingzao,SHI Liangsheng,ZHANG Dongxiao, et al.A sparse grid based Bayesian method for contaminant source identification[J].Advances in Water Resources,2012,37:1-9.

[10]WILLIAMS B,CHRISTENSEN W F,REESE C S. Pollution source direction identification:embedding dispersion models to solve an inverse problem[J]. Environmetrics,2011,22(8):962-974.

[11]WANG Hui,JIN Xin.Characterization of groundwater contaminant source using Bayesian method[J]. Stochastic Environmental Research and Risk Assessment,2012,27(4):867-876.

[12]OU Yang,WANG Xiaoyan.Identification of critical source areas for non-point source pollution in Miyun Reservoir watershed near Beijing,China[J].Water Science&Technology,2008,58(11):2235-2241.

[13]劉蘊(yùn),方晟,李紅,等.基于四維變分資料同化的核事故源項反演[J].清華大學(xué)學(xué)報(自然科學(xué)版),2015,55 (1):98-104.(LIU Yun,FANG Sheng,LI Hong,et al.Source inversion in nuclear accidents based on 4D variational data assimilation[J].Journal of Tsinghua University(Science and Technology),2015,55(1): 98-104.(in Chinese))

[14]軒曉博,逄勇,李一平,等.金屬礦區(qū)重金屬遷移對水體影響的數(shù)值模擬[J].水資源保護(hù),2015,31(2):30-35.(XUAN Xiaobo,PANG Yong,LI Yiping,et al. Numerical simulation of influence of heavy metal migration on water in metallic mining areas[J]. Water Resources Protection,2015,31(2):30-35.(in Chinese))

[15]DENG Zuicha,QIAN Kun,RAO Xiaobo,et.al.An inverse problem of identifying the source coefficient in a degenerate heat equation[J].Inverse Problems in Science and Engineering,2014,23(3):498-517.

[16]YANG Liu,DENG Zuicha,YU Jianning,et al. Optimization method for the inverse problem of reconstructing the source term in a parabolic equation [J].Mathematics and Computers in Simulation, 2009,80(2):314-326.

[17]FLETCHER R,REEVES C M.Function minimization by conjugate gradients[J].The Computer Journal, 1964,7(2):149-154.

[18]肖庭延,于慎根,王彥飛.反問題的數(shù)值解法[M].北京:科學(xué)出版社,2003:115-120.

[19]姚勝偉.幾類共軛梯度算法的研究[D].上海:華東理工大學(xué),2014.

[20]YAO Shengwei,LU Xiwen,NING Liangshuo,et al. A class of one parameter conjugate gradient methods [J].Applied Mathematics and Computation,2015, 265:708-722.

[21]范小平,李功勝.確定地下水污染強(qiáng)度的一種改進(jìn)的遺傳算法[J].計算物理,2007,24(2):187-191.(FAN Xiaoping,LI Gongsheng.An improved genetic algorithm for groundwater pollution[J].Chinese Journal of Computational Physics,2007,24(2):187-191.(in Chinese))

Reconstruction of groundwater pollution source term with improved conjugate gradient algorithm

XING Liying1,2,ZHANG Guozhen1
(1.School of Environmental and Municipal Engineering,Lanzhou Jiaotong University,Lanzhou 730070,China; 2.School of Civil and Architecture Engineering,Nanyang Normal University,Nanyang 473061,China)

In order to reconstruct the unsteady pollution source term of groundwater for its identification, the finite difference scheme was used to discretize the equations,terminal observed data were appended, and an improved conjugate gradient algorithm was used to reconstruct the pollution source term.The results show that,whether the pollution source term is a continuous function or a discontinuous and nondifferentiable function,it can be reconstructed.In addition,with consideration of measurement errors,the proposed method can reconstruct the source term effectively,indicating that the improved conjugate gradient algorithm is effective in reconstructing the groundwater pollution source term.

groundwater;pollution source term;reconstruction;inverse problem;conjugate gradient algorithm

X523

A

1004-6933(2017)03- 0042- 05

2016 06-08 編輯:王 芳)

10.3880/ji.ssn.1004-6933.2017.03.009

國家自然科學(xué)基金(51468033)

邢利英(1980—),女,講師,博士研究生,主要從事環(huán)境水力學(xué)方面研究。E-mail:lyxing2011@163.com