一維地下水溶質(zhì)運(yùn)移模型多參數(shù)反演

程 林，韓龍喜，劉曉華，王 琰，陳奇良

(1.河海大學(xué)環(huán)境學(xué)院，江蘇南京 210098;2.南通天虹環(huán)境科學(xué)研究所有限公司，江蘇南通 226600)

一維地下水溶質(zhì)運(yùn)移模型多參數(shù)反演

程 林1，韓龍喜1，劉曉華2，王 琰1，陳奇良1

(1.河海大學(xué)環(huán)境學(xué)院，江蘇南京 210098;2.南通天虹環(huán)境科學(xué)研究所有限公司，江蘇南通 226600)

以一維地下水溶質(zhì)運(yùn)移模型為基礎(chǔ)，利用微分進(jìn)化算法重點(diǎn)探討監(jiān)測(cè)頻次、測(cè)量誤差等因素對(duì)縱向彌散系數(shù)DL、有效孔隙度n以及水流速度u等參數(shù)聯(lián)合反演結(jié)果的影響，并構(gòu)造瞬時(shí)源和連續(xù)源兩種工況的典型案例來驗(yàn)證該方法的可靠性。結(jié)果表明:微分進(jìn)化算法具有收斂速度快、精度高等特點(diǎn)，能夠運(yùn)用于瞬時(shí)源和連續(xù)源地下水溶質(zhì)運(yùn)移模型水力、水質(zhì)參數(shù)的反演問題。

地下水;溶質(zhì)運(yùn)移;微分進(jìn)化算法;多參數(shù)反演;瞬時(shí)源;連續(xù)源

地下水水力、水質(zhì)參數(shù)識(shí)別是環(huán)境水力學(xué)反問題中一類重要的問題。確定地下水水質(zhì)模型的相關(guān)參數(shù)，對(duì)進(jìn)一步的水質(zhì)預(yù)報(bào)預(yù)警起到了至關(guān)重要的作用，具有重要的研究價(jià)值和應(yīng)用背景。反問題基本原理是根據(jù)地下水水質(zhì)模型的控制方程、初始條件、邊界條件以及環(huán)境水力學(xué)的測(cè)量數(shù)據(jù)(濃度數(shù)據(jù))來識(shí)別方程中的部分或全部參數(shù)。Wagner［1］首次采用非線性最大似然估計(jì)法同時(shí)反演模型參數(shù)與污染源特征;王錦國等［2］采用實(shí)數(shù)編碼遺傳算法優(yōu)化反演地下水污染物運(yùn)移參數(shù);韓一龍等［3］利用模擬退火遺傳算法，結(jié)合二維非均質(zhì)各向異性承壓含水層模型，對(duì)導(dǎo)水系數(shù)、貯水系數(shù)進(jìn)行反演計(jì)算;魏連偉等［4］基于遺傳算法這一全局優(yōu)化技術(shù)，以及地下水的有限元模型，應(yīng)用到北京市應(yīng)急水源地水文地質(zhì)參數(shù)的反演中;江思珉等［5］根據(jù)已有的監(jiān)測(cè)數(shù)據(jù)(水位和濃度等)，利用和聲搜索算法對(duì)地下水污染源和未知含水層參數(shù)進(jìn)行同步反演。

總體而言，關(guān)于地下水相關(guān)參數(shù)反問題主要集中于水位、貯水系數(shù)、滲透系數(shù)等水力參數(shù)的反演或者源強(qiáng)與上述參數(shù)的聯(lián)合反演，對(duì)于水力、水質(zhì)參數(shù)的聯(lián)合反演研究較少;且地下水反問題中對(duì)參數(shù)反演的方法自身研究較多，而關(guān)于地下水多參數(shù)聯(lián)合反演的文獻(xiàn)未見反演精度的影響因素研究。針對(duì)以上問題，本文利用微分進(jìn)化算法，結(jié)合一維地下水溶質(zhì)運(yùn)移模型，構(gòu)造瞬時(shí)源和連續(xù)源兩種工況下的典型案例，探討監(jiān)測(cè)頻次、測(cè)量精度等因素對(duì)縱向彌散系數(shù)DL、有效孔隙度n以及水流速度u等多參數(shù)聯(lián)合反演結(jié)果的影響。

1 地下水溶質(zhì)運(yùn)移微分進(jìn)化反演

1.1 微分進(jìn)化算法概述

微分進(jìn)化算法是Rainer Storn和Kenneth Price提出的一種更簡(jiǎn)單和更有效的基于種群進(jìn)化的多點(diǎn)搜索算法，同樣應(yīng)用了“優(yōu)勝劣汰，適者生存”的自然進(jìn)化法則，故應(yīng)當(dāng)歸屬于演化算法。筆者將一維地下水溶質(zhì)運(yùn)移模型的參數(shù)識(shí)別反問題轉(zhuǎn)化為非線性優(yōu)化問題，即將參數(shù)的反演問題轉(zhuǎn)化為求解式(1)待求優(yōu)化變量S來獲取最優(yōu)參數(shù)的過程。

式中:ρ(x，ti)為x處不同時(shí)刻溶質(zhì)質(zhì)量濃度計(jì)算值;ρobs(x，ti)為x處不同時(shí)刻溶質(zhì)質(zhì)量濃度觀測(cè)值。

1.2 約束條件

約束條件包括縱向彌散系數(shù)DL、有效孔隙度n與水流速度u的上下界約束以及式(2)的約束。

式中:DL為縱向彌散系數(shù)，cm2/min;αL為彌散度，cm;u為水流速度，cm/min。

彌散度是1個(gè)與實(shí)驗(yàn)尺度有關(guān)的參數(shù)，根據(jù)大量的彌散度土柱實(shí)驗(yàn)結(jié)果［6］，土柱實(shí)驗(yàn)測(cè)得的彌散度αL取值為［0.01，1］;有效孔隙度n取值范圍為(0，1.0);鑒于地下水流速土柱實(shí)驗(yàn)的數(shù)據(jù)缺少一些綜述性的資料，筆者采用前人的一些實(shí)驗(yàn)成果［7］，并在此基礎(chǔ)上以足夠大的數(shù)據(jù)作為水流速度u的一個(gè)上界約束，因此，水流速度u的取值范圍可為(0，10.0］。根據(jù)式(1)的約束關(guān)系，縱向彌散系數(shù)DL取值范圍為(0，10.0］。

1.3 反演步驟

a.產(chǎn)生初始種群。采用隨機(jī)方法產(chǎn)生初始種群:

式中:(0)為初始種群，Si(0)為種群中的第i個(gè)個(gè)體，i=1，2，…，N。

每個(gè)個(gè)體中有n個(gè)待求變量，那么Si(0)可表示為

其中，各個(gè)變量的初始值用隨機(jī)方法產(chǎn)生:

式中:x(i)j(0)為種群ˉX(0)中第i個(gè)個(gè)體Si(0)第j個(gè)待求變量的初始值;x(i)j，min(0)，x(i)j，max(0)分別為種群ˉX(0)中第i個(gè)個(gè)體Si(0)第j個(gè)待求變量的最小可能值和最大可能值。

b.計(jì)算適應(yīng)值。計(jì)算種群中各個(gè)個(gè)體Si(t)的適應(yīng)值f(Si(t))，在本問題中，適應(yīng)值為監(jiān)測(cè)點(diǎn)不同采樣時(shí)間上的溶質(zhì)濃度計(jì)算值和給定監(jiān)測(cè)值差值的平方和，然后取其平均值:

式中:W為監(jiān)測(cè)數(shù)據(jù)總個(gè)數(shù)。

c.繁殖。微分進(jìn)化算法采用特殊尋優(yōu)原則進(jìn)行繁殖，對(duì)種群中的每個(gè)個(gè)體Si(t)，通過隨機(jī)數(shù)方法分別生成3個(gè)兩兩不相等的隨機(jī)整數(shù)，r1，r2，r3∈{1，2，…，N}和隨機(jī)整數(shù)jrand∈{1，2，…n}，通過式(8)計(jì)算得到1個(gè)中間個(gè)體:

d.選擇進(jìn)化。當(dāng)且僅當(dāng)式(7)中產(chǎn)生的中間個(gè)體適應(yīng)值更好時(shí)，該個(gè)體才作為進(jìn)化的個(gè)體進(jìn)入下一代種群中，否則，仍然保留原來的個(gè)體，作為下一代生成中間個(gè)體的父?jìng)€(gè)體。即，

e.終止檢驗(yàn)。當(dāng)種群代數(shù)到達(dá)設(shè)定的最大代數(shù)時(shí)，算法結(jié)束，將適應(yīng)值最小的個(gè)體作為最優(yōu)解輸出。

2 瞬時(shí)源條件下一維地下水溶質(zhì)運(yùn)移模型的多參數(shù)識(shí)別

2.1 瞬時(shí)源一維地下水溶質(zhì)運(yùn)移模型正問題

不妨假定某含水層中(-∞ ＜x＜+∞)存在一維均勻流場(chǎng)，流動(dòng)方向?yàn)閤軸的正向，當(dāng)對(duì)該含水層瞬時(shí)投放示蹤劑進(jìn)行試驗(yàn)時(shí)，污染源可以概化為瞬時(shí)源。當(dāng)瞬時(shí)排放質(zhì)量為M的污染物質(zhì)時(shí)，排泄區(qū)下游污染物濃度變化規(guī)律的一維溶質(zhì)運(yùn)移模型控制方程為

式中:K為降解系數(shù)，1/d;w為橫截面面積，m2;M為污染物質(zhì)量，kg。

當(dāng)溶質(zhì)為不可降解物質(zhì)時(shí)，即K=0時(shí)，該方程通過變化得到如下解析解為

2.2 數(shù)值試驗(yàn)

算例1:在某多孔介質(zhì)柱體中投放10 g示蹤劑，柱體橫截面面積w為1m2，首先設(shè)定準(zhǔn)確的水力、水質(zhì)參數(shù)，其中有效孔隙度n為0.3，平均流速u為5.0 m/d，縱向彌散系數(shù)DL為0.5 m2/d。為確定參數(shù)，利用解析解模型計(jì)算得到的一組數(shù)據(jù)作為一組污染物濃度的監(jiān)測(cè)資料，其數(shù)學(xué)描述為

以監(jiān)測(cè)點(diǎn)位(x=1.0 m)前提下，分別假定污染物投放0.10 d后，每隔0.05 d監(jiān)測(cè)一次，計(jì)算結(jié)果見表1。采用本文所提出的微分進(jìn)化算法自動(dòng)識(shí)別參數(shù)DL、n和u，并分析監(jiān)測(cè)頻次、測(cè)量精度對(duì)反演結(jié)果的影響。

表1 瞬時(shí)源示蹤劑試驗(yàn)觀測(cè)數(shù)據(jù)

2.2.1 監(jiān)測(cè)頻次對(duì)參數(shù)反演結(jié)果的影響

為分析監(jiān)測(cè)頻次對(duì)參數(shù)反演精度的獨(dú)立影響，通過情景分析方法，以監(jiān)測(cè)頻次為1次、2次、3次、4次4種工況下進(jìn)行計(jì)算，結(jié)果見表2。

表2 不同監(jiān)測(cè)頻次瞬時(shí)源參數(shù)反演結(jié)果

2.2.2 測(cè)量誤差對(duì)參數(shù)反演結(jié)果的影響

不同的監(jiān)測(cè)儀器使得監(jiān)測(cè)的數(shù)據(jù)精度不同，為研究不同測(cè)量誤差情況下的微分進(jìn)化算法的反演精度情況，本次試驗(yàn)觀測(cè)數(shù)據(jù)誤差采用不同的小數(shù)位數(shù)來模擬不同儀器下的測(cè)量精度，為獨(dú)立反應(yīng)情況，4種工況都以監(jiān)測(cè)頻次3次為基礎(chǔ)，計(jì)算結(jié)果見表3。

由表2和表3的數(shù)值試驗(yàn)可以看出，監(jiān)測(cè)頻次和測(cè)量誤差對(duì)瞬時(shí)源工況下一維地下水溶質(zhì)運(yùn)移模型的縱向彌散系數(shù)DL、有效孔隙度n以及水流速度u等多參數(shù)的聯(lián)合反演結(jié)果有一定的影響，當(dāng)監(jiān)測(cè)數(shù)據(jù)較少，或者測(cè)量誤差較大時(shí)，反演結(jié)果出現(xiàn)不穩(wěn)定的情況，誤差較大;但當(dāng)監(jiān)測(cè)頻次大于2，且測(cè)量誤差多于1位時(shí)，反演結(jié)果具有較好的穩(wěn)定性，且能夠給出較高精度的反演結(jié)果。

表3 不同測(cè)量誤差瞬時(shí)源參數(shù)反演結(jié)果

3 連續(xù)源條件下一維地下水溶質(zhì)運(yùn)移模型的多參數(shù)識(shí)別

3.1 連續(xù)源一維地下水溶質(zhì)運(yùn)移模型正問題

假定半無限含水層中(0≤x＜+∞)存在一維流場(chǎng)，流動(dòng)方向?yàn)閤軸正向，在起始端(x=0 m)連續(xù)釋放示蹤劑，若不考慮混合過程段，則在該斷面處質(zhì)量濃度為ρ0，污染物的遷移擴(kuò)散規(guī)律可以概化為如下公式:

當(dāng)污染物為不可降解物質(zhì)時(shí)，即K=0時(shí)，通過拉普拉斯變換可得這一問題的解析解為

3.2 數(shù)值試驗(yàn)

算例2:本算例2的數(shù)據(jù)引自文獻(xiàn)［8］，原為采用erfc(x)近似公式求解地下水彌散系數(shù)的算例。在某多孔介質(zhì)柱體中連續(xù)投放濃度為ρ0=3 333.33 g/m3的示蹤劑，u=5.0 m/d，DL=0.5 m2/d，在x=1.0m處的不同時(shí)刻溶質(zhì)質(zhì)量濃度與初始質(zhì)量濃度的比值見表4。同時(shí)采用微分進(jìn)化算法反演參數(shù)DL、u，并分析監(jiān)測(cè)頻次、測(cè)量精度對(duì)反演結(jié)果的影響。

表4 不同時(shí)刻溶質(zhì)濃度與初始濃度的比值

3.2.1 監(jiān)測(cè)頻次對(duì)參數(shù)反演結(jié)果的影響

為分析監(jiān)測(cè)頻次對(duì)連續(xù)源工況下地下水相關(guān)參數(shù)反演精度的獨(dú)立影響，通過情景分析方法，以監(jiān)測(cè)頻次分別為1次、2次、3次、4次這4種工況下進(jìn)行計(jì)算，結(jié)果見表5。

表5 不同監(jiān)測(cè)頻次連續(xù)源參數(shù)反演結(jié)果

3.2.2 測(cè)量誤差對(duì)參數(shù)反演結(jié)果的影響

不同的監(jiān)測(cè)儀器使得監(jiān)測(cè)的數(shù)據(jù)精度不同，為研究不同測(cè)量誤差情況下的微分進(jìn)化算法的反演精度情況，本次試驗(yàn)觀測(cè)數(shù)據(jù)誤差采用不同的小數(shù)位數(shù)來模擬不同儀器下的測(cè)量精度，為獨(dú)立反應(yīng)情況，4種工況下都以監(jiān)測(cè)頻次3次為基礎(chǔ)，計(jì)算結(jié)果見表6。

表6 不同測(cè)量誤差連續(xù)源參數(shù)反演結(jié)果

由表5和表6的數(shù)值試驗(yàn)可以看出，監(jiān)測(cè)頻次和測(cè)量誤差對(duì)連續(xù)源一維地下水溶質(zhì)運(yùn)移模型的縱向彌散系數(shù)DL和水流速度u等參數(shù)聯(lián)合反演結(jié)果有一定的影響，當(dāng)監(jiān)測(cè)數(shù)據(jù)較少時(shí)，反演結(jié)果出現(xiàn)不穩(wěn)定的情況，誤差較大;而測(cè)量誤差對(duì)結(jié)果反演影響較小，但當(dāng)監(jiān)測(cè)頻次大于1時(shí)，反演結(jié)果具有較好的穩(wěn)定性，且能夠給出較高精度的反演結(jié)果。

4結(jié)語

a.從數(shù)值試驗(yàn)結(jié)果可知，微分進(jìn)化算法對(duì)瞬時(shí)源和連續(xù)源兩種工況下一維地下水溶質(zhì)運(yùn)移模型的水力、水質(zhì)參數(shù)的聯(lián)合反演能夠給出較高的精度。但監(jiān)測(cè)頻次和測(cè)量誤差對(duì)參數(shù)的反演精度有一定的影響，因此，利用該算法進(jìn)行反演時(shí)，監(jiān)測(cè)數(shù)據(jù)需要達(dá)到一定的要求。

b.在利用微分進(jìn)化算法的過程中發(fā)現(xiàn)，該算法具有收斂速度較快的優(yōu)點(diǎn)。

［1］ WAGNER B J.Simultaneously parameter estimation and contaminant source characterization for coupled groundwater flow and contaminant transport modeling［J］.Journal of Hydrology，1992，135:275-303.

［2］王錦國，周志芳，黃勇，等.基于實(shí)碼遺傳算法的地下水污染物運(yùn)移參數(shù)反演［J］.水文，2002，22(5):9-11.(WANG Jinguo，ZHOU Zhifang，HUANG Yong，et al.Inverse analysis for parameters of contamination transport in groundwater based on real-coded genetic algorithm［J］.Journal of China Hydrology，2002，22(5):9-11.(in Chinese))

［3］韓一龍，單永明.運(yùn)用模擬退火遺傳算法估計(jì)地下水反演參數(shù)［J］.計(jì)算機(jī)工程與應(yīng)用，2012，48(12):224-228.(HAN Yilong，SHAN Yongming. Using improved simulated annealing genetic algorithm to estimate parameters in groundwater inverse problem［J］.Computer Engineering and Applications，2012，48(12):224-228.(in Chinese))

［4］魏連偉，邵景力，張建力，等.遺傳算法在水文地質(zhì)參數(shù)反演中的應(yīng)用［J］.工程勘察，2004(3):28-31.(WEI Lianwei，SHAO Jingli，ZHANG Jianli，et al.Application of generic algorithm to inversion of hydrogeological parameters［J］.Geotechnical Investigation ＆ Surveying，2004(3):28-31.(in Chinese))

［5］江思珉，蔡奕，王敏，等.基于和聲搜索算法的地下水污染源與未知含水層參數(shù)的同步反演研究［J］.水利學(xué)報(bào)，2012，43(12):1470-1477.(JIANG Simin，CAI Yi，WANG Min，et al. Simultaneous identification of groundwater contaminant source and aquifer parameters by harmony search algorithm［J］. Journal of Hydraulic Engineering，2012，43(12):1470-1477.(in Chinese))

［6］馬建良，陳喜，程勤波，等.一維變密度溶質(zhì)運(yùn)移實(shí)驗(yàn)及參數(shù)推求［J］.水資源保護(hù)，2008，24(3):8-11.(MA Jianliang，CHEN Xi，CHENG Qinbo，et al.Identification of hydrodynamic parameters based on one dimensional variable density and solute transport numerical model［J］.Water Resources Protection，2008，24(3):8-11.(in Chinese))

［7］肖先煊，許模，蔡國軍，等.基于潛水滲流模型的地下水實(shí)際流速［J］.實(shí)驗(yàn)室研究與探索，2013，32(4):11-14.(XIAO Xianxuan，XU Mo，CAI Guojun，et al.On actual flow velocity of groundwater based on diving percolation model tests［J］.Research and Exploration in Laboratory，2013，32(4):11-14.(in Chinese))

［8］王軍輝，周志芳.erfc(x)近似公式及其在求解地下水運(yùn)移參數(shù)中的應(yīng)用［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，2001，29(3):111-114.(WANG Junhui，ZHOU Zhifang.Approximation of erfc(x)and its application in determining the groundwater transport parameters［J］.Journal of Hohai University:Natural Sciences，2001，29(3):111-114.(in Chinese))

Multi-parameter inversion of one-dimensional groundwater solute transport model

CHENG Lin1，HAN Longxi1，LIU Xiaohua2，WANG Yan1，CHEN Qiliang1
(1.College of Environment，Hohai University，Nanjing 210098，China;2.Nantong Tianhong Environmental Science Research Institute Co.，Ltd.，Nantong 226600，China)

Based on a one-dimensional groundwater solute transport model，the influences of monitoring frequency and measurement precision on joint inversion results of the longitudinal dispersion coefficientDL，the effective porosityn，and the flow velocityuwere studied using the differential evolution algorithm.Two typical cases，the instantaneous source and continuous source conditions，were constructed to verify the reliability of the differential evolution algorithm.Numerical experimental results show that the differential evolution algorithm has the characteristics of fast convergence and high precision，and it can be applied to the inversion of hydraulic and water quality parameters of solute transport models under the conditions of an instantaneous source and a continuous source.

groundwater;solute transport;differentialevolution algorithm;multi-parameterinversion;instantaneous source;continuous source

P641.2

A

1004-6933(2014)03-0005-04

10.3969/j.issn.1004-6933.2014.03.002

程林(1987—)，男，碩士研究生，研究方向?yàn)榄h(huán)境水力學(xué)。E-mail:chenglin0891@163.com

韓龍喜，教授。E-mail:hanlongxi@sina.com

(收稿日期:2013-07-31 編輯:高渭文)