李 璐, 殷樂宜, 牛浩博, 劉偉江, 陳 堅*

1.生態(tài)環(huán)境部環(huán)境規(guī)劃院, 長江經(jīng)濟帶生態(tài)環(huán)境聯(lián)合研究中心，北京 100012 2.生態(tài)環(huán)境部土壤與農(nóng)業(yè)農(nóng)村生態(tài)環(huán)境監(jiān)管技術中心，北京 100012

我國地下水環(huán)境風險源點多面廣，且地下水污染具有一定的隱蔽性，一旦發(fā)生污染，修復治理難度大且成本一般較高，但風險源周邊地下水監(jiān)測水平較低，因此，如何準確識別地下水污染來源對于地下水污染修復治理責任認定具有重要意義. 地下水污染源識別一般是基于一定的調(diào)查或監(jiān)測數(shù)據(jù)，推斷污染源的位置或污染路徑的過程[1]. 如何通過較少的監(jiān)測數(shù)據(jù)來識別污染來源，已經(jīng)成為近年來國內(nèi)外學者重點研究的領域.

目前，地下水污染源識別的主要方法包括地統(tǒng)計學法[2]、地球物理探測法[3-5]、同位素溯源法[6-7]和模型反演法等[8]. 其中，地統(tǒng)計學法要求多個位置甚至是長時間序列的監(jiān)測采樣結果，識別成本較高，一般用于已知污染源的地下水污染程度分析；地球物理探測法是基于地下水污染前后場地物性差異導致的物理場觀測值的變化進行污染分布情況分析的方法，但一般僅適用于與巖層存在較大物性差異以及30 m深度以內(nèi)的污染探測，且隨著精度要求越高，可探測范圍越淺；同位素溯源法目前已被廣泛應用于環(huán)境污染事件的仲裁、環(huán)境污染物的來源分析與示蹤中，但該方法只能對存在一定程度的特定無機物或有機物污染信息的污染源進行識別，且需要大量同位素和水化學監(jiān)測數(shù)據(jù)；模型反演方法中應用比較廣泛的求解方法主要包括模擬-優(yōu)化法[9-11]及不確定性分析方法等，其中，不確定性分析方法因?qū)ΡO(jiān)測數(shù)據(jù)需求量較小且可獲得相對可靠的反演結果，被越來越多的專家學者用于環(huán)境污染問題的反演識別，主要包括隨機游走粒子方法[12]、最小相對熵方法[13]、伴隨狀態(tài)方法[14]、貝葉斯統(tǒng)計方法[15-18]等.

貝葉斯統(tǒng)計方法以從監(jiān)測值中獲取參數(shù)信息為目標，將參數(shù)先驗概率密度函數(shù)與樣本似然函數(shù)結合在一起，形成一套非常靈活、直觀、易于理解的統(tǒng)計方法；此外，通過與其他分析方法聯(lián)用，還可以解決污染源貢獻率計算、反演精度提升等問題[19-20]. 張雙圣等[21-22]將貝葉斯模型與地下水二維水質(zhì)對流-擴散方程相耦合，建立了地下水污染源參數(shù)反演模型，實現(xiàn)了基于8個監(jiān)測點對1個污染源排放過程的反演. 然而，現(xiàn)實條件下很大概率上存在監(jiān)測數(shù)據(jù)的嚴重缺失，甚至只有單個監(jiān)測點一期數(shù)據(jù)的情況，為解決在監(jiān)測數(shù)據(jù)較少但存在多個地下水風險源的情景下，如何判定監(jiān)測數(shù)據(jù)異常的來源及風險源污染概率如何計算的問題，該研究提出了基于貝葉斯公式的地下水風險源污染概率估計研究方法，并開展案例分析，以期提供一種地下水污染來源判別方法的有效途徑，對地下水污染風險防控具有重要科學意義.

1 研究方法

1.1 貝葉斯概率計算公式

該研究利用貝葉斯公式，計算造成監(jiān)測點指標異常是由某一個風險源造成的概率，并將計算得到的概率最大值作為最優(yōu)解，從而形成地下水風險源污染概率估計方法.

貝葉斯模型可以表達為

P(A|B)=P(A)×P(B|A)/P(B)

(1)

式中：P(A|B)為后驗概率，代表在B監(jiān)測點m指標異常是由A風險源造成的概率；P(B|A)為似然度，代表A風險源造成B監(jiān)測點m指標異常的概率；P(A)為A風險源的先驗概率，代表A風險源發(fā)生m指標異常的概率；P(B)為標準化常量，代表B監(jiān)測點觀測到m指標出現(xiàn)異常的概率.

1.2 先驗概率修正

對于某個風險源來說，其是否發(fā)生污染取決于風險源是否排放目標污染物，而與B監(jiān)測點指標是否異常無關，在相關參數(shù)未知的情況下，風險源排放目標污染物和不排放目標污染物的概率相等，即可將初始先驗概率〔p0(Si)〕設為0.5. 然而，實際上各污染源性質(zhì)并不相同，若確定風險源排放目標污染物時，可將p0(Si)設置為1.0；若確定風險源不排放目標污染物時，可將p0(Si)設置為0；對于不同目標污染物來說，p0(Si)也各不相同. 除此之外，還需要利用風險源的原材料、生產(chǎn)工藝、產(chǎn)品、存在時間、廢水排放特征、環(huán)境保護設施運行情況等相關軟數(shù)據(jù)信息進行先驗概率修正. 對于Si風險源來說，修正后的先驗概率〔p(Si)〕與初始先驗概率〔p0(Si)〕、污染釋放可能性系數(shù)(Li)和污染物可能釋放量系數(shù)(Qi)有關，即

p(Si)=p0(Si)×Li×Qi

(2)

式(2)中先驗概率修正時所需數(shù)據(jù)通常來源于場地調(diào)查、環(huán)境統(tǒng)計、專家建議和文獻資料等. 該研究基于《地下水污染防治分區(qū)劃分工作指南》[23]中所述污染源荷載評價辦法，提出了利用風險源存在時間和廢水排放量等進行Li和Qi估算的方法，其中，Li與風險源存在時間有關，存在時間為(0,5]、(5,10]、(10,20]、(20,30]、>30 a時，Li可分別按照0.1、0.2、0.5、0.8、1.0取值；Qi與廢水排放量有關，廢水排放量為(0,1×104]、(1×104,1×105]、(1×105,5×105]、(5×105,1×106]、>1×106m3/a時，Qi可分別按照0.2、0.4、0.6、0.8、1.0取值.

若風險源為廢棄地塊，其廢水排放量可能為0，無法按照上述計算方法進行先驗概率修正，因而提出了基于《地下水污染源強評價、分類與防控技術研究》[24]中點源源強計算方法以及風險源特征的Li和Qi的計算方法，其中，Li與風險源的防滲措施運行情況有關. 若風險源已開展防滲，防滲措施運行(0,1]、(1,5]和>5 a時，Li可分別按照0.2、0.6、0.8取值，若未開展防滲的可將Li設置為1.0；Qi與滲水面積有關，滲水面積為(0,1×103]、(1×103,1×104]、(1×104,1×105]、(1×105,1×106]、>1×106m2時，Qi可分別按照0.2、0.4、0.6、0.8、1.0取值.

1.3 似然度計算

似然度的計算一般涉及風險源與監(jiān)測點的向量距離、水頭、流向、孔隙度、滲透系數(shù)、彌散度、滲漏濃度、監(jiān)測點觀測值等參數(shù). 似然度也是一種對計算值和觀測值之間接近程度的度量，計算值和觀測值越接近，則似然度的值就越大，反之亦然. 若已知多點或多次觀測值，可通過極大似然思想和溶質(zhì)運移方程求解似然函數(shù)[25-26]，但受監(jiān)測條件限制而無法準確求解溶質(zhì)運移方程時，也可通過對流彌散方程對似然度進行簡化求解.

假設研究區(qū)地下水流可概化為一個等厚的各向同性均質(zhì)含水層穩(wěn)定流，根據(jù)流場特征可知，似然度p(Si,m)可以表示為

p(Si,m)∝cosα×Δhi/Li2

(3)

也可以表示為

p(Si,m)/W=cosα×Δhi/Li2

(4)

式中：α為Si風險源和監(jiān)測點連線與地下水流向的夾角；Δhi為水頭差，m；Li為Si風險源與監(jiān)測點的距離，m；W為常數(shù).

1.4 標準化常量及后驗概率計算

標準化常量為歸一化的積分常數(shù)，可以表示為

(5)

也可以表示為

(6)

則后驗概率公式可以表示為

(7)

根據(jù)式(7)計算得到的后驗概率最優(yōu)解則為可能造成目標污染物超標的污染源[1].

2 案例分析

2.1 研究區(qū)概況

該研究選取石家莊市某工業(yè)集聚區(qū)作為研究對象，該區(qū)域位于滹沱河沖洪積扇平原地帶，地勢平坦，屬華北暖溫帶半濕潤地區(qū)，受大陸性季風氣候影響，年均氣溫12.2 ℃，年均降水量474.0 mm. 該工業(yè)集聚區(qū)于2000年底批復建設的裝備制造業(yè)及重化工產(chǎn)業(yè)基地，入駐了包括無機鹽制造、化學農(nóng)藥制造、其他合成材料制造、其他基礎化學原料制造、有機化學原料制造等行業(yè)類別在內(nèi)的多家企業(yè)，目前部分企業(yè)已停產(chǎn).

通過2019年4月現(xiàn)場調(diào)查采樣發(fā)現(xiàn)，該區(qū)域東南部某一農(nóng)灌井中，Cr6+含量超過GB 5749—2006《生活飲用水衛(wèi)生標準》指標限值及GB/T 14848—2017《地下水質(zhì)量標準》Ⅲ類標準限值3倍多，超過GB 5084—2005《農(nóng)田灌溉水質(zhì)標準》指標限值1倍多；CHCl3有檢出，但其含量未超過GB/T 14848—2017《地下水質(zhì)量標準》Ⅲ類標準限值. 為查找該監(jiān)測點Cr6+和CHCl3含量異常的原因，擬根據(jù)筆者所建立的風險源污染概率估計方法進行分析.

基于收集到的該工業(yè)集聚區(qū)內(nèi)企業(yè)的建成時間、廢水排放量、原材料、產(chǎn)品等相關軟數(shù)據(jù)，得到研究區(qū)域內(nèi)8個可能造成污染事件的風險源(記為S1～S8)，風險源及監(jiān)測點(記為O1)分布位置見圖1. 研究區(qū)含水層巖性主要為中砂、黏土互層，滲透系數(shù)為15～20 m/d，含水層平均厚度為18 m，地下水由西北流向東南，與區(qū)域多年水位方向一致.

2.2 風險源先驗概率計算

通過對風險源相關軟數(shù)據(jù)進行分析，與Cr6+含量和CHCl3含量具有極強相關性的風險源為S6和S3，其中，S6風險源的原材料包括鉻礦，鉻礦及其礦渣中所含Cr6+極易通過雨水淋濾作用進入地下水，而S3風險源的原材料中包括CHCl4，在進入土壤或地下水中極易降解為CHCl3等. 8個風險源中，由于S6風險源廢水排放量為0，需按廢棄地塊進行先驗概率計算. 8個風險源的行業(yè)類別、存在時間、廢水排放量等軟數(shù)據(jù)信息及先驗概率計算結果見表1.

2.3 污染事件似然度計算

根據(jù)研究區(qū)水文地質(zhì)條件，筆者將地下水流概化為均質(zhì)各向同性的潛水平面二維穩(wěn)定流. 根據(jù)風險源與監(jiān)測點的相對位置、水頭差等數(shù)據(jù)，利用式(4)計算污染事件發(fā)生的似然度與標準化常量的比值〔p(Si,m)/W〕，結果見表2.

圖1 地下水風險源及監(jiān)測點位置及地下水位等值線示意Fig.1 Sketch of the groundwater risk sources location, the monitoring point location and the groundwater table

表1 研究區(qū)域內(nèi)風險源相關信息及先驗概率值

2.4 后驗概率計算結果分析

基于2.2節(jié)和2.3節(jié)計算得到的先驗概率和似然度，按照式(7)計算得到不同指標異常時不同風險源的后驗概率，計算結果如表3所示.

由表3可知，O1監(jiān)測點Cr6+含量異常由S6風險源造成的后驗概率最大，即O1監(jiān)測點Cr6+含量異常由S6風險源造成的概率為76.2%，由此表明，監(jiān)測點Cr6+含量異常最有可能由某無機鹽制造業(yè)污染源造成；O1監(jiān)測點CHCl3含量異常由S1和S3風險源造藥制造業(yè)污染源造成.

表2 似然度計算結果

成的后驗概率分別為32.7%和23.6%，由此表明，監(jiān)測點CHCl3含量異常最有可能由一個或兩個化學農(nóng)

表3 后驗概率計算結果

進一步調(diào)研發(fā)現(xiàn)，S6風險源于2014年5月停產(chǎn)，鉻渣貯存庫房頂棚、圍墻破損嚴重，庫房內(nèi)大量鉻渣堆存尚未處置，2016年5月已對其實施掛牌督辦. 但此次研究結果顯示，Cr6+污染羽已超過1 500 m，存在巨大安全隱患，亟需治理，進一步佐證了監(jiān)測點Cr6+污染最有可能由S6風險源造成. 但S1風險源和S3風險源從生產(chǎn)工藝過程中涉及的原輔料和產(chǎn)物來說，S3風險源主要以吡啶、四氯化碳、氰化鈉、百草谷二氯鹽原藥等生產(chǎn)有機農(nóng)藥，S1風險源主要以阿維菌素、吡蟲啉、粉劑、乳劑等生產(chǎn)農(nóng)藥制品，但是粉劑和乳劑的成分不清，需要進一步補充分析S1風險源是否釋放或產(chǎn)生CHCl3，因此，S1和S3風險源是否造成CHCl3含量異常仍需多個監(jiān)測點或多期監(jiān)測數(shù)據(jù)進行反演驗證.

2.5 討論

該研究建立的地下水污染來源識別方法通過將風險源與監(jiān)測點之間污染物遷移的相關關系進行了簡化和求解，在此過程中篩選了關鍵參數(shù)進行分析，主要目的是在有限數(shù)據(jù)中挖掘最可能解，在地下水污染發(fā)現(xiàn)初期鎖定最可能的污染來源，以期減少常規(guī)地下水污染調(diào)查的工作量，也可以使后期的地下水污染調(diào)查更具有針對性，為地下水風險管控提供了有效方法和手段.

但是，在研究過程中還發(fā)現(xiàn)一些可能影響判定結果的情景，需要下一步重點研究：①當存在行業(yè)相似、后驗概率結果無顯著差異的兩個或多個風險源時，需要視現(xiàn)場條件開展少量補充監(jiān)測，并利用監(jiān)測結果進行后驗概率反演，進一步利用不確定分析理論甄別污染來源，在監(jiān)測數(shù)據(jù)充足時還可實現(xiàn)對污染釋放過程的反演[27-28]. ②當存在水位波動較大等情況時，利用穩(wěn)定流進行污染來源判定可能造成準確性下降，需收集區(qū)域內(nèi)多年水位變化情況，分時段進行概率分析再進行加權概率疊加. ③當區(qū)域空間異質(zhì)性顯著或存在優(yōu)勢通道時，則無法將研究區(qū)水流簡單概化為均質(zhì)穩(wěn)定流進行風險源污染概率計算，需要通過大量水文地質(zhì)資料及水位監(jiān)測數(shù)據(jù)進行研究區(qū)概念模型更新，并利用溶質(zhì)運移方程得到濃度計算值與實測值的擬合結果，通過不斷迭代優(yōu)化似然函數(shù)[21,26]，理論上可以實現(xiàn)非均質(zhì)條件下的污染溯源，這也是目前和未來一段時間相關研究的重要方向，同時進一步提升不確定性分析方法在實際案例中的應用水平.

3 結論

a) 利用已知的多個風險源的建成時間、廢水排放量、原材料、產(chǎn)品等相關軟數(shù)據(jù)信息修正先驗概率，并基于對流彌散方程提出似然度的簡化求解方法，提出了基于貝葉斯模型的地下水風險源污染概率估計方法，用于計算得到監(jiān)測點指標異常來源于每個風險源的后驗概率，通過結果判定監(jiān)測值異常的最可能來源，為我國目前監(jiān)測水平較低現(xiàn)狀下污染來源識別問題的解決提供了有效途徑.

b) 對石家莊市某工業(yè)集聚區(qū)內(nèi)8個風險源下游一個農(nóng)灌井中Cr6+含量和CHCl3含量異常事件的分析發(fā)現(xiàn)，Cr6+含量異常來源于S6風險源的后驗概率為76.2%，即Cr6+含量異常最有可能由某無機鹽制造業(yè)污染源造成；CHCl3含量異常來源于S1和S3風險源的后驗概率分別為32.7%和23.6%，監(jiān)測點CHCl3含量異常最有可能由一個或兩個化學農(nóng)藥制造業(yè)污染源造成.

c) 該研究僅針對有限信息下的污染來源進行概率判斷，風險源是否造成地下水污染仍需多個監(jiān)測點數(shù)據(jù)進行驗證分析，并結合溶質(zhì)運移方程繼續(xù)開展特殊情景下的方法優(yōu)化研究，使其更加滿足地下水污染溯源分析和責任認定應用需求.