李 璐, 殷樂(lè)宜, 牛浩博, 劉偉江, 陳 堅(jiān)*

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

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

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

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

1 研究方法

1.1 貝葉斯概率計(jì)算公式

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

貝葉斯模型可以表達(dá)為

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

(1)

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

1.2 先驗(yàn)概率修正

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

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

(2)

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

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

1.3 似然度計(jì)算

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

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

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

(3)

也可以表示為

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

(4)

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

1.4 標(biāo)準(zhǔn)化常量及后驗(yàn)概率計(jì)算

標(biāo)準(zhǔn)化常量為歸一化的積分常數(shù)，可以表示為

(5)

也可以表示為

(6)

則后驗(yàn)概率公式可以表示為

(7)

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

2 案例分析

2.1 研究區(qū)概況

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

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

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

2.2 風(fēng)險(xiǎn)源先驗(yàn)概率計(jì)算

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

2.3 污染事件似然度計(jì)算

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

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

表1 研究區(qū)域內(nèi)風(fēng)險(xiǎn)源相關(guān)信息及先驗(yàn)概率值

2.4 后驗(yàn)概率計(jì)算結(jié)果分析

基于2.2節(jié)和2.3節(jié)計(jì)算得到的先驗(yàn)概率和似然度，按照式(7)計(jì)算得到不同指標(biāo)異常時(shí)不同風(fēng)險(xiǎn)源的后驗(yàn)概率，計(jì)算結(jié)果如表3所示.

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

表2 似然度計(jì)算結(jié)果

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

表3 后驗(yàn)概率計(jì)算結(jié)果

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

2.5 討論

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

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

3 結(jié)論

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

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

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