與實(shí)際問題不同,假想例子沒有實(shí)際監(jiān)測(cè)數(shù)據(jù).所以將真實(shí)的污染源特征,輸入污染物質(zhì)運(yùn)移模擬模型,正向運(yùn)行模擬模型得到觀測(cè)井的污染物質(zhì)濃度,作為反向識(shí)別污染源特征時(shí)的實(shí)測(cè)數(shù)據(jù).水流和污染物質(zhì)在第1800和3600d的空間分布情況分別見圖3和圖4.觀測(cè)井的污染物質(zhì)濃度見圖5.

圖5 觀測(cè)井的污染物質(zhì)濃度實(shí)測(cè)數(shù)據(jù)

Fig.5 Measured data of contaminant concentration in observation wells

2.3 靈敏度分析

參數(shù)不確定性會(huì)影響污染物質(zhì)運(yùn)移模擬模型的輸出結(jié)果.考慮模型所有參數(shù)對(duì)輸出結(jié)果的影響,又會(huì)導(dǎo)致模型過于復(fù)雜.因此,應(yīng)用局部靈敏度分析方法[21],篩選出對(duì)模擬模型輸出結(jié)果影響最大的參數(shù),作為輸入模擬模型的隨機(jī)變量,其它參數(shù)作為模擬模型的確定性變量.

篩選隨機(jī)變量時(shí),首先,將輸入模擬模型的參數(shù)取均值(表3)輸入模型,運(yùn)行模型,得到觀測(cè)井的污染物質(zhì)濃度.然后,將輸入模型的某個(gè)參數(shù)加減10%,20%,其它參數(shù)保持不變.再次運(yùn)行模型,得到某一參數(shù)變化后,對(duì)應(yīng)的觀測(cè)井污染物質(zhì)濃度.利用式(2)計(jì)算出各參數(shù)的靈敏度系數(shù),選擇靈敏度系數(shù)最大的參數(shù),作為輸入模型的隨機(jī)變量.靈敏度分析結(jié)果見圖6.

表3 參數(shù)服從的概率分布及取值范圍

Table 3 Probability distribution and value range of parameters

由圖6可以看出,靈敏度系數(shù)最大的參數(shù)是滲透系數(shù).因此,將滲透系數(shù)作為輸入模型的隨機(jī)變量,其它參數(shù)作為確定性變量,取定值輸入模型.根據(jù)前人經(jīng)驗(yàn)[24],滲透系數(shù)服從的概率分布見表3.利用拉丁超立方抽樣方法[22-23]對(duì)滲透系數(shù)在給定范圍內(nèi)抽樣600組,每組參數(shù)樣本具有5個(gè)值,分別對(duì)應(yīng)研究區(qū)的5個(gè)參數(shù)分區(qū).同時(shí),利用該方法對(duì)污染物質(zhì)釋放強(qiáng)度抽樣600組,每組樣本具8個(gè)值,分別是2個(gè)污染源在4個(gè)釋放時(shí)段的污染物質(zhì)釋放強(qiáng)度.使?jié)B透系數(shù)和釋放強(qiáng)度隨機(jī)組合,得到600組由滲透系數(shù)和釋放強(qiáng)度組成的輸入樣本.

圖6 參數(shù)靈敏度分析

Fig.6 Parameter sensitivity analysis graph

應(yīng)用模擬-優(yōu)化方法識(shí)別地下水污染源時(shí),污染物質(zhì)運(yùn)移模擬模型會(huì)作為等式約束連接到優(yōu)化模型中,求解優(yōu)化模型時(shí),成百上千次的調(diào)用模擬模型,會(huì)產(chǎn)生大量的計(jì)算負(fù)荷,浪費(fèi)大量的計(jì)算時(shí)間.建立模擬模型的替代模型,能有效解決這一問題.

2.4 建立替代模型

應(yīng)用克里格方法建立污染物質(zhì)運(yùn)移模擬模型的替代模型.參考克里格方法的原理,利用MATLAB編寫克里格替代模型的訓(xùn)練程序.然后,建立污染物質(zhì)運(yùn)移模擬模型的替代模型.建立替代模型的步驟如下:

1.將600組輸入樣本依次輸入污染物質(zhì)運(yùn)移模擬模型,運(yùn)行模擬模型,可以得到與600組輸入樣本一一對(duì)應(yīng)的600組觀測(cè)井的污染物質(zhì)濃度.

2.利用540組輸入樣本作為替代模型的輸入,污染物質(zhì)的濃度作為替代模型的輸出,訓(xùn)練替代模型.

3.將余下的60組輸入樣本,依次輸入替代模型,得到替代模型輸出的60組觀測(cè)井的污染物質(zhì)濃度.利用60組替代模型的輸出和模擬模型的輸出,檢驗(yàn)替代模型的精度.

應(yīng)用確定性系數(shù)和平均相對(duì)誤差檢驗(yàn)替代模型的精度,它們的計(jì)算方法見式14和15:

式中:yi是模擬模型輸出的第i個(gè)時(shí)段末的污染物質(zhì)濃度,mg/L;是替代模型輸出的第i個(gè)時(shí)段末的污染物質(zhì)濃度,mg/L;是模擬模型輸出的各時(shí)段末的污染物質(zhì)濃度均值,mg/L;n是觀測(cè)樣本的總數(shù)量.

2.5 建立優(yōu)化模型

建立替代模型后,需要建立識(shí)別地下水污染源特征的非線性優(yōu)化模型.非線性優(yōu)化模型包括目標(biāo)函數(shù)、決策變量和約束條件3部分.將觀測(cè)井的污染物質(zhì)實(shí)測(cè)濃度與模擬計(jì)算濃度擬合(替代模型的輸出濃度)最小平方誤差和,作為目標(biāo)函數(shù);將污染物質(zhì)的釋放強(qiáng)度,作為優(yōu)化模型的決策變量;建立識(shí)別地下水污染源特征的非線性優(yōu)化模型.優(yōu)化模型的約束條件包括污染源釋放強(qiáng)度不等式約束、5個(gè)分區(qū)滲透系數(shù)和地下水溶質(zhì)運(yùn)移規(guī)律等式約束;優(yōu)化模型表示如下:

式中:Ki是某個(gè)參數(shù)分區(qū)的滲透系數(shù),m/d;qm是污染物質(zhì)釋放強(qiáng)度,mg/d;是在觀測(cè)井處污染物質(zhì)的模擬計(jì)算濃度,mg/L;是在觀測(cè)井處的污染物質(zhì)實(shí)測(cè)濃度,mg/L.是替代模型.

利用遺傳算法[25]求解某滲透系數(shù)組約束下的優(yōu)化模型,只能得到與該組滲透系數(shù)影響下的污染源識(shí)別結(jié)果.而研究滲透系數(shù)不確定性對(duì)污染源識(shí)別結(jié)果的影響,需要統(tǒng)計(jì)多組滲透系數(shù)樣本影響下的污染源識(shí)別結(jié)果.因此,應(yīng)用蒙特卡羅方法進(jìn)行考慮滲透系數(shù)不確定性的污染源識(shí)別.

2.6 建立推算模型

應(yīng)用蒙特卡羅的思想,將優(yōu)化模型作為“試驗(yàn)”發(fā)生器,將某滲透系數(shù)樣本影響下的污染源識(shí)別結(jié)果,作為某次隨機(jī)“試驗(yàn)”的結(jié)果.通過對(duì)成千上萬(wàn)組隨機(jī)“試驗(yàn)”結(jié)果的統(tǒng)計(jì)與分析,考慮滲透系數(shù)不確定性對(duì)污染源識(shí)別結(jié)果的影響.

建立優(yōu)化模型后,再次應(yīng)用拉丁超立方方法,對(duì)滲透系數(shù)抽樣90組.將第1組滲透系數(shù),賦值到優(yōu)化模型中,然后求解優(yōu)化模型(整個(gè)求解過程中,滲透系數(shù)一直保持不變),得到與第1組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度(包括S1的4個(gè)污染物質(zhì)釋放強(qiáng)度值和S2的4個(gè)污染物質(zhì)釋放強(qiáng)度值).對(duì)第2組滲透系數(shù)做與第1組滲透系數(shù)相同處理,直到第90組滲透系數(shù).這樣就依次求解了90個(gè)優(yōu)化模型,獲得了90組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度.

但是,僅僅考慮90組滲透系數(shù)對(duì)污染源識(shí)別結(jié)果的影響,很難具有代表性.而求解成千上萬(wàn)個(gè)優(yōu)化模型,又會(huì)產(chǎn)生海量的計(jì)算負(fù)荷,花費(fèi)大量的計(jì)算時(shí)間.因此提出了應(yīng)用克里格方法,建立能夠描述滲透系數(shù)與污染物質(zhì)釋放強(qiáng)度之間關(guān)系的推算模型,來(lái)推算不同滲透系數(shù)取值影響下的污染物質(zhì)釋放強(qiáng)度.將70組滲透系數(shù)和與之影響下的污染物質(zhì)釋放強(qiáng)度,分別作為推算模型的輸入與輸出,訓(xùn)練推算模型.利用余下的20組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度作為檢驗(yàn)數(shù)據(jù)(每組包括S1的4個(gè)污染物質(zhì)釋放強(qiáng)度值和S2的4個(gè)污染物質(zhì)釋放強(qiáng)度值),檢驗(yàn)推算模型的精度.推算模型與替代模型的精度評(píng)價(jià)指標(biāo)相同.

3 結(jié)果與討論

3.1 替代模型的精度評(píng)價(jià)

表4 替代模型的R2和MRE

Table 4R2and MRE of surrogate model

圖7 7口觀測(cè)井的污染物質(zhì)濃度擬合

Fig.7 Fitting graph of contaminant concentration of seven observation wells

(a)~(g)分別表示第1到第7口觀測(cè)井

對(duì)基于克里格方法替代模型的精度進(jìn)行了評(píng)價(jià).通常情況下,如果替代模型的R2高于0.9,MRE在10%以內(nèi),則認(rèn)為替代模型精度滿足研究需要.結(jié)合表4和圖7,圖8可以看出,替代模型與模擬模型輸出濃度的平均相對(duì)誤差均小于2%,確定性系數(shù)都大于0.99,接近1,替代模型與模擬模型的輸出結(jié)果擬合程度非常好,替代模型的精度很高,可以代替模擬模型連接到優(yōu)化模型中.

由圖8可知,利用7口觀測(cè)井濃度數(shù)據(jù)進(jìn)行替代模型的精度檢驗(yàn)時(shí),出現(xiàn)了個(gè)別的異常值,但是異常值的相對(duì)誤差均未超過4%,在10%以內(nèi),可見即使對(duì)于異常值,替代模型對(duì)模擬模型的擬合程度仍然很高.出現(xiàn)異常值的可能原因是,建立替代模型的訓(xùn)練數(shù)據(jù),在個(gè)別的輸入樣本取值及取值領(lǐng)域內(nèi)覆蓋的不夠全面,導(dǎo)致替代模型在該輸入樣本取值處泛化性能欠佳.因此,在獲取訓(xùn)練數(shù)據(jù)時(shí),盡量去覆蓋輸入樣本的取值范圍是很有必要的.

圖8 7口觀測(cè)井的污染物質(zhì)濃度相對(duì)誤差箱線圖

Fig.8 The relative error box plot of the contaminant concentration of seven observation wells

3.2 識(shí)別結(jié)果的不確定性分析

圖9 推算模型擬合曲線

Fig.9 Fitting graph of the inference model

應(yīng)用克里格方法,建立了描述滲透系數(shù)與污染物質(zhì)釋放強(qiáng)度之間關(guān)系的推算模型.由圖9可知, 推算模型的確定性系數(shù)達(dá)到了0.9895,與1非常接近;平均相對(duì)誤差為4.51%,在10%以內(nèi);推算模型的精度很高,可以用來(lái)推算任意滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度.

考慮滲透系數(shù)不確定性的地下水污染源識(shí)別,需要統(tǒng)計(jì)成百上千組滲透系數(shù)影響下的污染源識(shí)別結(jié)果.具體做法是:再次對(duì)滲透系數(shù)抽樣8000組,將其全部輸入推算模型,可以得到各組滲透系數(shù)影響下的污染物質(zhì)釋放強(qiáng)度.對(duì)8000組污染物質(zhì)釋放強(qiáng)度進(jìn)行統(tǒng)計(jì)與分析.分析因?yàn)闈B透系數(shù)不確定性,導(dǎo)致的污染源識(shí)別結(jié)果的不確定性.

由圖10可以看出,受滲透系數(shù)不確定性影響,污染源的識(shí)別結(jié)果也具有不確定性.應(yīng)用切比雪夫不等式[22]估計(jì)不同置信水平對(duì)應(yīng)的污染源識(shí)別結(jié)果置信區(qū)間.由表5可知,不同置信程度對(duì)應(yīng)的置信區(qū)間,均包含2個(gè)污染源的實(shí)際釋放強(qiáng)度.隨著置信程度增大,置信區(qū)間變大,反之則反.決策者可以根據(jù)不同的管理目標(biāo),選擇相信不同置信程度對(duì)應(yīng)的置信區(qū)間,做為污染源特征的識(shí)別結(jié)果.也可以選擇概率密度最大的釋放強(qiáng)度作為污染源識(shí)別結(jié)果(識(shí)別結(jié)果見表6).

圖10 不確定性分析

Fig.10 Uncertainty analysis diagram

(a)~(d)是S1各個(gè)釋放時(shí)段的污染物質(zhì)釋放強(qiáng)度;(e)~(h)是S2各個(gè)釋放時(shí)段的污染物質(zhì)釋放強(qiáng)度

表5 不同置信水平對(duì)應(yīng)的污染源識(shí)別結(jié)果置信區(qū)間

Table 5 Confidence interval of contamination source identification results corresponding to different confidence levels

由表6可知,識(shí)別得到的S1和S2的釋放強(qiáng)度與真實(shí)的污染源特征雖然有一定的誤差,但是整體上對(duì)污染源的真實(shí)特征逼近程度較好.其中,S1第4時(shí)段的釋放強(qiáng)度,與其真實(shí)釋放強(qiáng)度差距最大,導(dǎo)致這樣結(jié)果的原因可能有以下2方面:(1)替代模型和推算模型均具有誤差,會(huì)在一定程度上影響識(shí)別結(jié)果的準(zhǔn)確性.(2)抽樣的過程具有隨機(jī)性,抽樣的樣本不一定會(huì)完美覆蓋到真實(shí)的含水層滲透系數(shù).含水層滲透系數(shù)偏離真值,會(huì)導(dǎo)致污染源的識(shí)別結(jié)果偏離污染源特征的真實(shí)結(jié)果.

求解一個(gè)優(yōu)化模型,進(jìn)行1000次迭代計(jì)算,需要的時(shí)間約為0.45h,求解8000個(gè)優(yōu)化模型需要3600h.而建立推算模型需要約40.5h,應(yīng)用推算模型推算8000組滲透系數(shù)影響下的污染源特征,僅需要不到1s(瞬間輸出).應(yīng)用推算模型可以節(jié)省約99%的計(jì)算時(shí)間和計(jì)算負(fù)荷.

表6 概率密度最大的污染源識(shí)別結(jié)果

Table 6 Identification results of contamination sources with the highest probability density

4 結(jié)論

4.1 基于克里格方法建立的替代模型具有較高的精度,確定性系數(shù)高于0.99,平均相對(duì)誤差小于1.2%.將克里格替代模型做為等式約束條件連接到優(yōu)化模型中,供求解優(yōu)化模型調(diào)用,在保證一定精度的前提下,可以減少大量的負(fù)荷和時(shí)間.

4.2 應(yīng)用克里格方法,建立了描述滲透系數(shù)與污染物質(zhì)釋放強(qiáng)度之間關(guān)系的推算模型.建立的推算模型具有較高的精度,確定性系數(shù)和平均相對(duì)誤差分別為0.9895和4.51%.運(yùn)用克里格推算模型計(jì)算了8000滲透系數(shù)影響下的污染源特征,節(jié)省了約99%的計(jì)算負(fù)荷和時(shí)間.

4.3 對(duì)8000組污染源特征進(jìn)行定量的統(tǒng)計(jì)與分析,統(tǒng)計(jì)得到了置信水平為80%、60%、40%和20%對(duì)應(yīng)的污染源識(shí)別結(jié)果置信區(qū)間;分析出S1在4個(gè)釋放時(shí)段概率密度最大的釋放強(qiáng)度分別是74.22×105、53.67×105、40.35×105和16.36 ×105mg/d,S2在4個(gè)釋放時(shí)段概率密度最大的釋放強(qiáng)度分別是55.08× 105、50.50×105、39.41×105和21.01×105mg/d,將概率密度最大的釋放強(qiáng)度做為污染源識(shí)別結(jié)果,對(duì)污染源的真實(shí)特征有較好的逼近程度.

4.4 將模擬-優(yōu)化方法、靈敏度分析方法、蒙特卡羅方法和克里格方法結(jié)合,進(jìn)行考慮滲透系數(shù)不確定性的污染源識(shí)別研究,更加符合實(shí)際情況,而且能夠得到更加豐富的污染源識(shí)別結(jié)果,為決策者提供更多的參考依據(jù).