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

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

Fig.5 Measured data of contaminant concentration in observation wells

2.3 靈敏度分析

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

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

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

Table 3 Probability distribution and value range of parameters

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

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

Fig.6 Parameter sensitivity analysis graph

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

2.4 建立替代模型

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

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

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

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

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

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

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

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

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

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

2.6 建立推算模型

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

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

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

3 結果與討論

3.1 替代模型的精度評價

表4 替代模型的R2和MRE

Table 4R2and MRE of surrogate model

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

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

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

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

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

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

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

3.2 識別結果的不確定性分析

圖9 推算模型擬合曲線

Fig.9 Fitting graph of the inference model

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

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

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

圖10 不確定性分析

Fig.10 Uncertainty analysis diagram

(a)~(d)是S1各個釋放時段的污染物質(zhì)釋放強度;(e)~(h)是S2各個釋放時段的污染物質(zhì)釋放強度

表5 不同置信水平對應的污染源識別結果置信區(qū)間

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

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

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

表6 概率密度最大的污染源識別結果

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

4 結論

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

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

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

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