基于改進的共軛梯度算法重構地下水污染源項

邢利英,張國珍

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

基于改進的共軛梯度算法重構地下水污染源項

邢利英1,2,張國珍1

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

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

地下水;污染源項;重構;反問題;共軛梯度法

隨著現代化的快速發展,大量工農業污染物的排放導致地下水污染嚴重。地下水中污染物的遷移過程作用時間長,影響范圍大,在遷移過程中還涉及水文地理條件、物理反應、化學反應等諸多因素,使得污染物遷移過程變得極為復雜,并表現出非線性行為。由于地下水系統本身的隱藏性和復雜性,地下水污染并沒有引起業界足夠的重視。最近的一份調查報告顯示,作為主要淡水資源之一的地下水已經被氮、碳氫化合物、有毒有害的微量元素嚴重污染,并且呈現出由點到面、由城鎮向農村發展的趨勢,所以此類問題變得越來越棘手[1-4]。

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

共軛梯度法在數學理論分析中應用較多,Deng等[15]采用共軛梯度方法重構了傳熱方程的輻射項系數,數值結果表明該方法有效地反演了輻射項系數;Yang等[16]應用共軛梯度方法重構傳熱方程的源項。本文設計一種改進的共軛梯度方法反演地下水污染源項。

1 地下水污染控制方程

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

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

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

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

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

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

2 差分方程

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

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

式中:ct、cxx、cx分別為污染物濃度c關于時間的一階導數、空間的二階導數、空間的一階導數。四點隱格式的差分方程為

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

3 改進的共軛梯度算法

共軛梯度算法由Fletcher等[17]首次提出的,作為一種簡單有效的迭代技術,廣泛應用于地理學和傳熱學方程的源項及擴散項系數反演。在大中型線性和非線性的無約束優化問題中,共軛梯度算法尤其有效。共軛梯度算法的基本思路為結合共軛性和最速下降法,構建一組共軛負梯度方向,沿著下降的方向尋找目標函數的最小值,合適終止條件的共軛梯度法實際上是一類迭代正則化方法。鑒于共軛梯度的特點,本文擬進一步改進共軛梯度算法重構地下水污染源項。

設非線性優化問題的控制方程可以表示為

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

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

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

詳細的改進共軛梯度算法步驟[15,20]如下:

步驟1:設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:增加迭代次數k,執行循環到步驟2。任取γ>1,當k∈Ν,第一次滿足‖gk-g‖L2(0,L)≤γδ時終止循環。

4 算 例

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

當式(18)取第一組參數:a(x)=1,φ(x)=sin (2πx),q(x)=2π2sin(2πx),時間步長和空間步長分別取1/100s和1/100 m,應用改進的共軛梯度算法反演污染源項的結果見圖1(a)。從圖1(a)中可以看出,當迭代次數K=50時,反演結果已經很好地與真實解吻合。

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

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

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

圖1 試驗1的反演結果

應用改進的共軛梯度算法進行反演,附加數據和反演結果的對比見圖2。從圖2中可以看出當迭代次數K=30時,污染源項已經較好地被重構,充分說明該方法收斂速度較快,反演結果穩定。采用cδT(x)=cδ(x,T)=c(x,T)[1+δrandom(x)]計算測量誤差,當測量誤差δ分別為2%和5%時,污染源項的真值和反演結果的計算誤差分別為3.04%和2.83%,比參考文獻[21]的誤差4.18%要小。

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

由于污染源項的不連續性,當迭代次數增大到100次時,污染源項才被較好地重構。類似地,考慮測量誤差δ的影響時,盡管污染源項的不連續性,而且測量誤差也會產生影響,但是模擬結果也比較正常,充分證明改進共軛梯度算法的穩定性和較大的容錯性。

5 結 語

圖2 實際問題的反演結果與附加數據的比較

圖3 試驗3的反演結果

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

[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]孫才志,陳相濤,陳雪姣,等.地下水污染風險評價研究進展[J].水利水電科技進展,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]朱嵩.基于貝葉斯推理的環境水力學反問題研究[D].杭州:浙江大學,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法的水質模型參數不確定性研究[J].環境科學,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]李子,毛獻忠,周康.污染物非恒定輸運逆過程反演模型研究[J].水力發電學報,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]劉蘊,方晟,李紅,等.基于四維變分資料同化的核事故源項反演[J].清華大學學報(自然科學版),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]軒曉博,逄勇,李一平,等.金屬礦區重金屬遷移對水體影響的數值模擬[J].水資源保護,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]肖庭延,于慎根,王彥飛.反問題的數值解法[M].北京:科學出版社,2003:115-120.

[19]姚勝偉.幾類共軛梯度算法的研究[D].上海:華東理工大學,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]范小平,李功勝.確定地下水污染強度的一種改進的遺傳算法[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

國家自然科學基金(51468033)

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