王蕾，張芳，王成鋒，徐振飛

（四川省地質(zhì)工程勘察院集團(tuán)有限公司，四川 成都 610000）

0 引言

場(chǎng)地地下水污染具有隱蔽性、復(fù)雜性和突發(fā)性等特點(diǎn)，導(dǎo)致污染物的來源識(shí)別存在較大困難，快速識(shí)別污染源發(fā)生的位置對(duì)于場(chǎng)地修復(fù)治理、防控措施的制定具有重要的研究?jī)r(jià)值（張辰等，2018）。要展開地下水污染治理工作，首要步驟是根據(jù)有限的觀測(cè)數(shù)據(jù)結(jié)合數(shù)據(jù)模型方法，找出污染源的位置及污染物的運(yùn)移路徑（楊海東等，2014）。目前地下水污染物溯源的方法，主要分為兩大類：地球化學(xué)足跡法和數(shù)學(xué)模擬法（龍玉橋等，2017）。由于污染物的擴(kuò)散和運(yùn)移是一種隨著時(shí)間演變的過程，故單獨(dú)使用地球化學(xué)足跡法不能完全解決污染物的溯源問題。而數(shù)學(xué)模擬法的研究多集中于簡(jiǎn)單的理想算例（趙勇勝，2012），實(shí)際工程應(yīng)用中，需要考慮研究場(chǎng)地的含水層參數(shù)、污染物類型等因素的影響（陳正俠等，2017），導(dǎo)致數(shù)學(xué)模型的建立對(duì)使用者的專業(yè)素質(zhì)要求較高，這限制了它在普通研究人員和用戶中的推廣使用。

本文提供一種場(chǎng)地地下水污染物自動(dòng)溯源算法，該算法利用地球化學(xué)足跡法和數(shù)學(xué)模型相結(jié)合的方式，通過彌散實(shí)驗(yàn)獲取研究場(chǎng)地水文地質(zhì)參數(shù)（沈婷婷，2015），將場(chǎng)地巡檢操作和污染物排查經(jīng)驗(yàn)量化為數(shù)學(xué)模型，該模型在研究場(chǎng)地的應(yīng)用將推進(jìn)類似場(chǎng)地地下水監(jiān)測(cè)信息服務(wù)的智能化，提高地下水污染源分析研判的效率，縮短污染事故的應(yīng)急處置時(shí)間，有效減少地下水污染帶來的環(huán)境破壞，在一定程度上減少因事件而造成的經(jīng)濟(jì)損失和社會(huì)影響。

1 研究區(qū)概況

研究區(qū)地下水資源豐富，區(qū)內(nèi)淺層地下水屬孔隙潛水，受自然地質(zhì)環(huán)境條件的限制，天然防護(hù)條件差，防污性能差，地下水環(huán)境脆弱。當(dāng)生產(chǎn)裝置發(fā)生滲漏時(shí)，如不及時(shí)監(jiān)測(cè)、預(yù)警和采取相應(yīng)處理措施，可能造成嚴(yán)重的地下水污染及快速向下徑流，影響到廠區(qū)下游城鎮(zhèn)及分散農(nóng)戶供水安全（姜鳳成等，2017；陳書客等，2019）。廠區(qū)面積約4.6 km2，廠外排水管線長(zhǎng)約77 km，廠內(nèi)合計(jì)監(jiān)測(cè)及應(yīng)急抽水井168 口，廠外排水管線及沿線監(jiān)測(cè)井25 口，共計(jì)193 口監(jiān)測(cè)及應(yīng)急抽水井。針對(duì)研究區(qū)的地下水監(jiān)測(cè)，已建成地下水監(jiān)測(cè)原始數(shù)據(jù)庫(kù)和預(yù)警數(shù)據(jù)庫(kù)，開發(fā)了地下水監(jiān)測(cè)信息管理系統(tǒng)，但隨著地下水監(jiān)測(cè)工作的深入，廠區(qū)地下水污染防控已經(jīng)由單純的監(jiān)測(cè)預(yù)警發(fā)展到污染物溯源、溶質(zhì)運(yùn)移等綜合分析方面，如何實(shí)現(xiàn)地下水異常時(shí)的及時(shí)溯源是目前工作區(qū)關(guān)注的重點(diǎn)之一，也是目前地下水監(jiān)測(cè)領(lǐng)域研究的難點(diǎn)之一。

本文針對(duì)研究區(qū)的水文地質(zhì)條件，研究一套自動(dòng)溯源模型并將其設(shè)計(jì)轉(zhuǎn)化為計(jì)算機(jī)算法，該套算法根據(jù)地下水流場(chǎng)、各裝置特征污染物、各監(jiān)測(cè)井水質(zhì)情況進(jìn)行分析判斷，初步實(shí)現(xiàn)污染源追溯功能。算法將適用于具有類似地下水污染風(fēng)險(xiǎn)、布設(shè)有大量監(jiān)測(cè)孔的孔隙型場(chǎng)地地下水污染防控工作中。旨在減少由于外部因素對(duì)監(jiān)測(cè)結(jié)果造成的延判、誤判，為地下水防控提供更加合理、科學(xué)的決策依據(jù)（孫江虎等，2016；孫杰等，2019）。

2 模型方法

2.1 彌散試驗(yàn)應(yīng)用

通過彌散試驗(yàn)得到污染物從試驗(yàn)主孔到各觀測(cè)孔的最大和最小時(shí)間，從而利用時(shí)間位移公式計(jì)算出示蹤劑到達(dá)觀測(cè)孔的最小和最大速度。進(jìn)一步的，利用觀測(cè)到的水位數(shù)據(jù)和達(dá)西公式，計(jì)算監(jiān)測(cè)孔對(duì)應(yīng)的水力梯度和滲透系數(shù)。最后利用試驗(yàn)主孔和發(fā)現(xiàn)污染監(jiān)測(cè)孔之間水力梯度、最大滲透系數(shù)，計(jì)算出污染物通過污染監(jiān)測(cè)孔的可能最大速度，從而計(jì)算出污染物向下游可能運(yùn)移到的最大位移（馬佳等，2020；潘云鶴，2001）。

（1）采用多孔試驗(yàn)法開展彌散試驗(yàn)，試驗(yàn)主孔記為ZK1，主流向各觀測(cè)孔記為ZGC1、ZGC2……ZGCn，側(cè)向各觀測(cè)孔記為CGC1、CGC2……CGCn。

（2）彌散試驗(yàn)中主流向觀測(cè)孔ZGC1、ZGC2……ZGCn第一次捕捉到示蹤劑的時(shí)間分別記為TZGC1、TZGC2、……TZGCn，各主流向觀測(cè)孔至ZK1的距離分別記為L(zhǎng)ZGC1、LZGC2、……LZGCn，計(jì)算示蹤劑到達(dá)觀測(cè)孔的最大速度Vmax。

（3）彌散試驗(yàn)中側(cè)向（無限接近垂直）觀測(cè)孔CGC1、CGC2……CGCn第一次捕捉到示蹤劑的時(shí)間分別記為TCGC1、TCGC2、……TCGCn，各側(cè)向觀測(cè)孔至ZK1 的距離分別記為L(zhǎng)CGC1、LCGC2、……LCGCn，計(jì)算示蹤劑到達(dá)觀測(cè)孔的最小速度Vmin。

（4）記Vmax和Vmin的對(duì)應(yīng)觀測(cè)點(diǎn)為GCvmax、GCvmin，計(jì)算機(jī)程序利用觀測(cè)點(diǎn)的水位數(shù)據(jù)通過等值線算法繪制穿過GCvmax、GCvmin的等水位線，得到GCvmax、GCvmin的水位值Hvmax、Hvmin，進(jìn)而計(jì)算出GCvmax、GCvmin點(diǎn)位和試驗(yàn)主孔（ZKGC）之間對(duì)應(yīng)的水力梯度Ivmax、Ivmin。

其中，L為試驗(yàn)主孔到觀測(cè)孔的距離，HZ KGC為試驗(yàn)主孔的水位值。

（5）利用達(dá)西公式，計(jì)算含污染物的流體在GCvmax、GCvmin點(diǎn)位上的滲透系數(shù)K主、K側(cè)。

（6）針對(duì)發(fā)現(xiàn)污染的監(jiān)測(cè)點(diǎn)Ji（i=1、2、3...n），根據(jù)第（4）步繪制的等水位線，記Ji 的水位數(shù)據(jù)為HJi，進(jìn)而計(jì)算出Ji水力梯度IJi。

其中，LJ i為試驗(yàn)主孔到Ji的距離。

（7）利用達(dá)西公式，計(jì)算污染物通過Ji的可能最大速度VJimax。

（8）計(jì)算場(chǎng)地裝置發(fā)生泄漏時(shí)，在監(jiān)測(cè)時(shí)間間隔T內(nèi)，污染物沿地下水流向運(yùn)移到監(jiān)測(cè)點(diǎn)Ji的可能最大距離LJimax。

2.2 半空間模型

在一個(gè)空間內(nèi)，半空間模型將所有樣本分成兩類，對(duì)于任意一個(gè)樣本xi∈X（X是一個(gè)含有m個(gè)樣本的樣本集合，i=1，2，…，m）只會(huì)有兩個(gè)狀態(tài)，即f（xi）≥0 或f（xi）＜0（王玉琨等，2008；鄭德濤等，1996）。在研究區(qū)的場(chǎng)地上進(jìn)行地下水污染物的溯源模擬，溯源對(duì)象（裝置）是否產(chǎn)生污染的狀態(tài)數(shù)據(jù)正好滿足半空間模型的樣本條件（Willem F.Bronsvoort et al.，1993），即裝置泄露和裝置未泄露。為了將裝置泄露的狀態(tài)數(shù)據(jù)利用數(shù)學(xué)表達(dá)式將其中的連續(xù)變量表達(dá)出來（方德洲和李淼，2007），建立研究區(qū)適用的半空間模型，劃分半空間Ai、Ai’，i=1、2、3…n，確定污染物溯源的有效搜索半空間是最終溯源的關(guān)鍵。具體實(shí)現(xiàn)方法如下：

（1）以發(fā)現(xiàn)污染監(jiān)測(cè)點(diǎn)Ji為圓心，LJimax為半徑畫圓，該圓圈定的范圍即污染源初始范圍OA。

（2）以Ji監(jiān)測(cè)點(diǎn)等水位線（水位值HJi）為分割法線，從空間上將污染源初始范圍OA分割為兩個(gè)半空間A1、A1’，A1內(nèi)的所有監(jiān)測(cè)點(diǎn)的水位值小于或等于HJi，即位于Ji下游；A1’內(nèi)的所有監(jiān)測(cè)點(diǎn)的水位值大于HJi，即位于Ji上游。

（3）標(biāo)定兩個(gè)半空間的屬性，其中半空間A1置屬性flag（A1）=1、半空間置A1’ 置屬性flag（A1）=0。

（4）針對(duì)屬性flag（A1）=1 的半空間A1內(nèi)所有監(jiān)測(cè)點(diǎn)進(jìn)行搜索，如果存在污染點(diǎn)，則以水位值最小污染點(diǎn)作為圓心，重復(fù)上述（1）～（3）步驟的操作，直至區(qū)域Ai無其他污染點(diǎn)存在。此時(shí)，半空間Ai’即為有效搜索半空間。

2.3 鎖定污染源

若有效搜索半空間Ai’內(nèi)只有一個(gè)裝置，則該裝置為污染源泄露裝置。若有多個(gè)裝置，則在Ai’內(nèi)搜索遍歷所有裝置，并利用半空間模型縮小污染源所在范圍。具體步驟如下：

（1）利用現(xiàn)有場(chǎng)地的等水位線內(nèi)插得到穿過裝置中心點(diǎn)的等水位線，并將裝置的水位記為HZ j，進(jìn)而根據(jù)觀測(cè)主孔和裝置之間的水頭差和距離計(jì)算出裝置對(duì)應(yīng)的水力梯度IZ j。

（2）利用達(dá)西公式和2.1 節(jié)中計(jì)算得出的K側(cè)，得到通過裝置Zj的可能最小流速VZjmin。

（3）計(jì)算監(jiān)測(cè)間隔時(shí)間T內(nèi)，污染物沿地下水流向運(yùn)移到該點(diǎn)的可能最小距離LZjmin。

（4）以Zj為圓心，以該裝置發(fā)生污染后污染物在T內(nèi)的最小運(yùn)移距離LZjmin為半徑，繪制圓形OB。

（5）以穿過該裝置的等水位線HZj為分割法線，從空間上將OB分割為兩個(gè)空間Bj、Bj’，Bj內(nèi)的所有監(jiān)測(cè)點(diǎn)的水位值小于裝置點(diǎn)水位值，即位于裝置下游；Bj’內(nèi)的所有監(jiān)測(cè)點(diǎn)的水位值大于或等于裝置點(diǎn)水位值，即位于裝置上游。

（6）標(biāo)定兩個(gè)半空間的屬性，其中半空間Bj置flag（B）=1；半空間Bj’置flag（B）=0，然后過濾掉flag（B）=0 無效半空間，保留有效半空間Bj。

（7）搜索Ai與Bj的交集區(qū)域內(nèi)除了Ji以外的監(jiān)測(cè)井，若存在無污染的監(jiān)測(cè)井，則排除該裝置是污染源的可能；若無監(jiān)測(cè)井，或者有監(jiān)測(cè)井且監(jiān)測(cè)井全部污染，則裝置Zj 有可能是污染源，記為Wk，k=1、2、3…m，m為可能發(fā)生滲漏的污染裝置總數(shù)。

（8）Wk所在的范圍就是污染源所在范圍。

（9）將監(jiān)測(cè)到的污染物類型與Wk內(nèi)的特征污染物類型進(jìn)行比對(duì)，具有同類型的污染物的裝置即為泄露裝置，即鎖定污染源。

3 試驗(yàn)案例

3.1 試驗(yàn)數(shù)據(jù)

本課題選擇的研究區(qū)早在2013 年就建立了地下水監(jiān)測(cè)信息管理系統(tǒng)，系統(tǒng)包含自動(dòng)化監(jiān)測(cè)和人工監(jiān)測(cè)兩大部分，為研究區(qū)的地下水監(jiān)測(cè)數(shù)據(jù)的傳輸、查詢、分析、上報(bào)等工作提供了便捷的管理方式。但隨著業(yè)務(wù)的持續(xù)開展，研究區(qū)地下水污染防控的需求已經(jīng)由單純的監(jiān)測(cè)預(yù)警發(fā)展到污染物溯源、溶質(zhì)運(yùn)移等綜合分析方面，而本課題的開展正是為了解決地下水中發(fā)現(xiàn)污染后污染物溯源的難題，課題組將第2 部分的模型方法編寫為計(jì)算機(jī)算法，并形成單獨(dú)的功能模塊供系統(tǒng)調(diào)用，而基礎(chǔ)數(shù)據(jù)（如監(jiān)測(cè)井、裝置等）則調(diào)用現(xiàn)有的系統(tǒng)數(shù)據(jù)庫(kù)，課題的實(shí)驗(yàn)全部在該模塊中展開。

表1 試驗(yàn)數(shù)據(jù)

試驗(yàn)實(shí)施前，課題組通過彌散試驗(yàn)已得到該場(chǎng)地示蹤劑到達(dá)觀測(cè)孔的最小和最大速度，并隨機(jī)選擇了兩口污染監(jiān)測(cè)井J1，J2，在系統(tǒng)操作界面錄入污染監(jiān)測(cè)井編號(hào)后，系統(tǒng)依照本算法流程開始了如下計(jì)算。

第一步：調(diào)取經(jīng)過彌散試驗(yàn)得出示蹤劑到達(dá)觀測(cè)孔的最大、最小速度Vmax、Vmin，并自動(dòng)定位對(duì)應(yīng)觀測(cè)孔，取出對(duì)應(yīng)觀測(cè)孔的水位進(jìn)而計(jì)算孔的水力梯度Ivmax、Ivmin，利用達(dá)西公式，計(jì)算出含污染物的流體在GCvmax、GCvmin點(diǎn)位上的滲透系數(shù)K主、K側(cè)。最后，針對(duì)本次發(fā)現(xiàn)的污染點(diǎn)J1，J2，計(jì)算其水力梯度IJ1、IJ2和污染物通過這兩個(gè)監(jiān)測(cè)孔的最大速度VJ1max、VJ2max。輸入J1，J2的記錄觀測(cè)時(shí)間間隔T，計(jì)算出污染物可能運(yùn)移到檢測(cè)孔的最大距離LJ1max、LJ2max。試驗(yàn)數(shù)據(jù)如表1 所示。

第二步：彌散試驗(yàn)應(yīng)用完成后，算法將繼續(xù)執(zhí)行，利用改進(jìn)后的半空間模型確定本次試驗(yàn)的有效半空間，分別為A1′和A2′，由于J2在J1的下游，最終確定A2′為有效搜索半空間。

第三步：結(jié)合裝置再次利用改進(jìn)后的半空間模型進(jìn)行污染源范圍的縮小。在半空間A2’內(nèi)搜索到裝置Z1、Z2，發(fā)現(xiàn)Z1內(nèi)對(duì)應(yīng)污染監(jiān)測(cè)井J1，Z2內(nèi)對(duì)應(yīng)無污染監(jiān)測(cè)井J3。即可判斷出Z1是疑似滲漏裝置，Z2不是泄露裝置（表2）。

表2 污染源鎖定

3.2 效果分析

通過將上述模型算法應(yīng)用于研究區(qū)現(xiàn)有的地下水監(jiān)測(cè)系統(tǒng)中，經(jīng)驗(yàn)證，該模型算法可行，起到了確定或縮小污染源范圍的作用。圖1 為算法驗(yàn)證流程圖。圖2 在監(jiān)測(cè)系統(tǒng)上模型運(yùn)行的結(jié)果，最終鎖定一個(gè)裝置為疑似污染對(duì)象，大大縮小了排查范圍，提高了現(xiàn)場(chǎng)工作人員的工作效率。

4 結(jié)論

（1）本文選用的研究場(chǎng)地是項(xiàng)目組成員多年開展地下水監(jiān)測(cè)、巡檢和維護(hù)項(xiàng)目的實(shí)施地，是布設(shè)有大量監(jiān)測(cè)孔的孔隙型場(chǎng)地。筆者對(duì)場(chǎng)地的地質(zhì)、水文地質(zhì)條件以及基本水文地質(zhì)參數(shù)等資料均已獲取，故利用彌散試驗(yàn)對(duì)污染物的最大、最小運(yùn)移半徑進(jìn)行計(jì)算是較好的選擇。

（2）將專家經(jīng)驗(yàn)量化為數(shù)學(xué)模型，可以在一定程度上彌補(bǔ)當(dāng)前數(shù)學(xué)模型理想算例脫離現(xiàn)實(shí)的情況，同時(shí)將其轉(zhuǎn)化為地下水監(jiān)測(cè)系統(tǒng)中的溯源算法并提供系統(tǒng)后臺(tái)服務(wù)，用戶通過填寫巡檢過程中的相關(guān)參數(shù)來運(yùn)行服務(wù)，對(duì)用戶的專業(yè)深入程度要求不高，可操作性較強(qiáng)。

（3）半空間模型結(jié)合場(chǎng)地裝置的泄露狀態(tài)數(shù)據(jù)，可以很好地發(fā)揮半空間模型的優(yōu)勢(shì)，同時(shí)也能提高算法的搜索的效率和準(zhǔn)確率。

圖1 算法驗(yàn)證流程圖

圖2 污染物溯源效果

（4）溯源解析結(jié)果和現(xiàn)場(chǎng)實(shí)際經(jīng)驗(yàn)研判基本吻合，因此，基于半空間模型和彌散試驗(yàn)的場(chǎng)地地下水污染物自動(dòng)溯源算法是可行的。該算法的應(yīng)用可為同類場(chǎng)地地下水污染事件的處理提供參考。