信昆侖，盛希夫，陶 濤，項寧銀

（同濟大學(xué) 環(huán)境科學(xué)與工程學(xué)院，上海 200092）

近幾年，城市供水管網(wǎng)的突發(fā)污染事故屢有發(fā)生，除水源地污染引起管網(wǎng)飲用水問題之外，如管道錯接等供水管網(wǎng)內(nèi)部的原因也會造成水質(zhì)的局部乃至全局惡化.因此，如何快速準(zhǔn)確地定位污染源是飲用水安全保障急需解決的課題.關(guān)于污染源的定位研究，Laird等［1］提出了用非線性規(guī)劃方法來求解這一問題.首先用源頭追蹤算法替代EPANET水質(zhì)模擬計算引擎從而產(chǎn)生必要的數(shù)據(jù)，然后通過帶正則項的非線性規(guī)劃方法求出污染源侵入的時間和地點.針對這一方法解的非唯一性問題，Laird等［2］又提出了先用混合整數(shù)二次規(guī)劃的方法確定一個候選污染源節(jié)點集，然后在此污染源節(jié)點集的基礎(chǔ)上利用前一方法篩選出最有可能是污染源的節(jié)點.但這個方法對于大中型管網(wǎng)的適用性和此模型與實際管網(wǎng)相符程度仍然有待研究.Guan等［3］利用EPANET作為內(nèi)嵌的水質(zhì)模擬器再結(jié)合負梯度法提出了模擬優(yōu)化算法.這一方法在已知管網(wǎng)模型的運行模式和水質(zhì)監(jiān)測數(shù)據(jù)充足可靠的前提下能夠較好地解出污染源侵入的時間和地點.但只考慮了水質(zhì)反應(yīng)為線性的情況，也沒有給出確定候選污染源節(jié)點的方法.Pries等［4］針對此問題提出了關(guān)系樹-線性規(guī)劃算法.通過構(gòu)造上下游節(jié)點的污染物濃度關(guān)系樹并形成線性規(guī)則，再利用關(guān)系樹的線性規(guī)則反向求解線性規(guī)劃問題，得到污染源節(jié)點位置及污染注入屬性.Cristo等［5］提出了離散最優(yōu)化的方法來求解污染源的識別問題.主要通過構(gòu)造水質(zhì)污染矩陣來確定可能的污染源注入點，再通過最小化模擬值和監(jiān)測值的差來從可能的污染源點中確定污染源的位置.另外，Liu等［6］同樣利用EPANET作為內(nèi)嵌的水質(zhì)模擬計算引擎用遺傳算法隨機搜索最優(yōu)解.由于以上方法均需要管網(wǎng)的水質(zhì)監(jiān)測點可以報告實時的污染物濃度，對于水質(zhì)監(jiān)測點僅能判斷是否出現(xiàn)污染物的情況，Yang等［7］提出用貝葉斯方法計算已檢出污染物的水質(zhì)監(jiān)測節(jié)點其上游節(jié)點為污染源的概率，然后根據(jù)概率的大小判斷出可能的污染源.

本質(zhì)上來說，供水管網(wǎng)污染源的識別問題是管網(wǎng)水質(zhì)模擬的反演問題，從數(shù)學(xué)的角度來看也是一個優(yōu)化問題.本文基于Guan等［3］利用負梯度法求解供水管網(wǎng)污染源識別問題，結(jié)合Cristo等［5］提出的候選污染源節(jié)點集合方法，進一步考慮供水管網(wǎng)水質(zhì)反應(yīng)的復(fù)雜性，用差分步長法代替權(quán)重系數(shù)法求解梯度方向，采用EPANET作為內(nèi)嵌的水質(zhì)模擬計算引擎結(jié)合負梯度法來求解污染源侵入的時間、地點和侵入過程，并結(jié)合實例管網(wǎng)進行算法求解性能的分析.

1 模型構(gòu)建

在城市供水管網(wǎng)擁有相當(dāng)精度的水力水質(zhì)模型且管網(wǎng)上安裝了若各個水質(zhì)監(jiān)測點的情況下，污染源的識別問題可以轉(zhuǎn)化為一個數(shù)學(xué)上的最優(yōu)化問題.本文以模擬時段內(nèi)各個水質(zhì)監(jiān)測點的實際污染物濃度值與模擬值差的平方和為目標(biāo)函數(shù)，通過分析水質(zhì)監(jiān)測點不同時刻模擬濃度值與實測濃度值差來修正計算候選污染源在特定時刻的節(jié)點濃度.待模擬時段內(nèi)所有的候選污染源節(jié)點濃度更新后，再次進行水質(zhì)模擬.當(dāng)模擬結(jié)果達到預(yù)設(shè)的條件后，輸出候選污染源節(jié)點濃度矩陣.由于求解所得的候選污染源污染物累計注入量與其作為真實污染源的可能性成正比，因此可以此作為確定污染源節(jié)點的依據(jù).

首先定義

式中：xi（t）表示候選污染源節(jié)點i在t時刻的污染物濃度值，對于k個時間段，則xi是一個行向量，即xi＝［xi（t1），xi（t2），…，xi（tk）］；yj（t）表示水質(zhì)監(jiān)測點j在t時刻的模擬污染物濃度值；n表示候選污染源節(jié)點總數(shù)；m表示布置的水質(zhì)監(jiān)測點總數(shù)；t表示模擬時刻.

監(jiān)測點模擬濃度可以表示為

結(jié)合監(jiān)測點的實際監(jiān)測結(jié)果，污染源識別問題的目標(biāo)函數(shù)可以表示為［3］

式中：X＊為污染源識別問題的最優(yōu)解，y′i（t）為水質(zhì)監(jiān)測點i在t時刻的實際濃度值.為了便于求解，這里需作幾點假設(shè)：①污染物的流量相對于管網(wǎng)水量而言可忽略不計，不影響管網(wǎng)原來的水力工況；②污染物在管網(wǎng)節(jié)點處投加，投加點可以是多個，投加量可以隨時間變化，但在每一水質(zhì)時間步長之內(nèi)是恒定不變的；③污染物侵入管網(wǎng)中每個節(jié)點的概率是相同的，也就是說污染物是隨機侵入的；④管網(wǎng)中的水質(zhì)監(jiān)測點能夠監(jiān)測到任意濃度的污染物.

2 模型求解

由于式（3）表示的最優(yōu)化問題目標(biāo)函數(shù)和約束條件均為非線性函數(shù)，且約束條件是由供水管網(wǎng)恒定流方程、連續(xù)性方程、管道水質(zhì)反應(yīng)動力學(xué)方程共同決定的復(fù)雜非線性等式約束，因此該問題為高度復(fù)雜的非線性規(guī)劃問題.考慮到約束條件可采用成熟的EPANET水力水質(zhì)分析模塊進行隱式求解，故基于負梯度法，通過差分法求目標(biāo)函數(shù)梯度向量對此非線性規(guī)劃問題進行求解.求解過程如下：

步驟1 給定初始近似點X（0）及精度ε＞0，若‖▽f（X（0））‖2≤ε，則X（0）即為近似極小點.需要說明的是梯度向量的每一個元素均為目標(biāo)函數(shù)對X＝［x1，x2，…，xn］相應(yīng)元素的偏導(dǎo)數(shù).該偏導(dǎo)數(shù)通過差分法近似求解

式中：ΔX是第i個元素為δ且其余元素為0的列向量；δ為差分步長，其取值大小應(yīng)該合理確定；f即目標(biāo)函數(shù)（見式（3））.當(dāng)水質(zhì)反應(yīng)為線性關(guān)系時，水質(zhì)監(jiān)測點所測得的污染物濃度值可以表達成污染物侵入節(jié)點的污染物濃度加權(quán)值［3］，即

偏導(dǎo)數(shù)的求法可以簡化為

式中：aj，i為水質(zhì)監(jiān)測點j對候選污染源節(jié)點i的權(quán)重系數(shù)，即當(dāng)候選污染源節(jié)點i有一個單位的增量時，在水質(zhì)監(jiān)測點j有相應(yīng)的污染物濃度增量.權(quán)重系數(shù)可以通過EPANET水質(zhì)分析中的源頭追蹤功能得到.

步驟2 若‖▽f（X（0））‖2＞ε，求步長λ0，并計算

求最佳步長本文采用黃金分割法，最佳步長即為使得目標(biāo)函數(shù)最小的步長.若‖▽f（X（1））‖2≤ε，則輸出X（1），否則轉(zhuǎn)入步驟3.

步驟3 一般地說，若‖▽f（X（k））‖2≤ε，則X（k）即為所求的近似解；若‖▽f（X（k））‖2＞ε，則求步長λk，并確定下一個近似點

如此迭代下去，直至達到所設(shè)定精度要求為止.

3 算例研究

以供水管網(wǎng)實例Net3為例進行污染源的識別問題求解.Net3管網(wǎng)模型含節(jié)點數(shù)92個，水庫2座，網(wǎng)中水塔3個，管段數(shù)117個，水泵2臺，如圖1所示.圖中虛線范圍內(nèi)的節(jié)點為候選污染源節(jié)點.

圖1 管網(wǎng)拓撲結(jié)構(gòu)圖Fig.1 Topology layout of the case study network

假定在節(jié)點101和119同時注入污染物，注入方式為設(shè)置點注入［8］，注入時間段為3∶00am～10∶00am，污染物注入質(zhì)量濃度變化模式見圖2.管網(wǎng)水質(zhì)模擬時間總長24h，水力時間步長1h，水質(zhì)時間步長5min.在管網(wǎng)中隨機選取6個節(jié)點作為水質(zhì)監(jiān)測點（分別為連接節(jié)點121，151，111，267，205，213），如圖1所示.通過水質(zhì)模擬，以這些節(jié)點模擬污染物質(zhì)量濃度數(shù)據(jù)作為監(jiān)測結(jié)果數(shù)據(jù)（見表1）.

表1 監(jiān)測點模擬污染物質(zhì)量濃度Tab.1 Simulated concentrations at monitoring sites mg·L-1

3.1 確定候選節(jié)點集

候選節(jié)點集是指根據(jù)水質(zhì)監(jiān)測點對污染物的監(jiān)測情況所確定的可能是管網(wǎng)污染事故污染源發(fā)生位置的節(jié)點集合.目前關(guān)于管網(wǎng)突發(fā)污染事故污染源追蹤的研究，仍然以通過經(jīng)驗法確定候選節(jié)點集合為主［3］，候選節(jié)點集合確定的范圍和準(zhǔn)確程度都會影響到污染源識別問題最優(yōu)化算法的性能.為此，本研究根據(jù)文獻［5，9］提出的最小覆蓋集方法進行候選節(jié)點集合的確定.首先求出節(jié)點的監(jiān)測時間矩陣，以24h為它們的服務(wù)水平，可以確定出每個監(jiān)測點所覆蓋的范圍.通過表1可知，本次污染事故假定模式下，實際上僅水質(zhì)監(jiān)測點151，111，267，205，213能監(jiān)測到污染物.因此可以把5個監(jiān)測點共同覆蓋的節(jié)點作為候選節(jié)點集（如圖1所示），即

3.2 結(jié)果分析

考慮到按照負梯度方向搜索得到的是局部最優(yōu)解，且絕大多數(shù)候選污染源節(jié)點均不是污染源注入點，可將初始解向量的每個元素均設(shè)置為0.求解過程中，差分步長δ取0.0100，最高迭代次數(shù)取360.圖2中顯示采用差分法計算梯度向量求解得到的結(jié)果，由于求解的多數(shù)候選點的質(zhì)量濃度值均很小，圖中僅顯示30個候選污染源節(jié)點中所求解的注入質(zhì)量濃度較大的4個節(jié)點的注入質(zhì)量濃度隨時間變化情況.候選污染源節(jié)點101是各時刻累計投加質(zhì)量濃度值最大的節(jié)點，初步選定為污染源節(jié)點.同時，候選污染源節(jié)點10從3∶00am到10∶00am基本上以100mg·L-1的質(zhì)量濃度持續(xù)注入.圖1表明，節(jié)點10和節(jié)點101是同一個管段的起止節(jié)點.該現(xiàn)象表明有一個污染源在以節(jié)點10和節(jié)點101為起止節(jié)點的管道中.事實上，候選污染源節(jié)點101的注入質(zhì)量濃度變化曲線基本與實際污染源節(jié)點101的質(zhì)量濃度變化曲線相吻合.此外，候選污染源節(jié)點119的投加質(zhì)量濃度呈拋物線型變化也與實際的污染源119投加質(zhì)量濃度變化相吻合，可以確定這里也是一個污染源注入點.另外，候選污染源節(jié)點151作為污染源節(jié)點119的鄰近節(jié)點，也有一定質(zhì)量濃度的污染物注入，這就更表明了另一個污染源就是節(jié)點119或者在其附近.

圖2 污染物注入質(zhì)量濃度變化曲線（差分法）Fig.2 Injection concentrations of contamination source（finite difference method）

3.3 影響求解結(jié)果的因素分析

（1）差分步長

對于圖1所顯示的小型算例管網(wǎng)，將迭代次數(shù)均設(shè)定為360次，差分步長δ依次取值為0.1000，0.0100，0.0010和0.0001時，計算結(jié)果的梯度向量模長值依次為16.52，0.24，0.15和0.33.為了定量地衡量計算結(jié)果與實際結(jié)果的擬合程度，計算污染源節(jié)點101和119各個時段假定注入濃度和解出注入濃度差的平方和總和，依次為：1564.09，1914.74，7703.17，64201.58，如表2所示.

由表2可知，差分步長取0.0100時對應(yīng)的計算結(jié)果（梯度向量模長為0.24，目標(biāo)函數(shù)值為1914.74，而梯度向量模長和目標(biāo)函數(shù)值是反映梯度法求解質(zhì)量的兩個主要指標(biāo)）最好，此時差分步長取值正好與EPANET水質(zhì)選項中設(shè)定的水質(zhì)公差值相同.其原因在于水質(zhì)公差表示管網(wǎng)水質(zhì)的最小變化值，也即一部分水體與另一部分水體水質(zhì)的最小差值［9］，這個值也表征了管網(wǎng)中水質(zhì)監(jiān)測點的污染物最低監(jiān)測限度，因此當(dāng)差分步長取水質(zhì)公差值時，計算結(jié)果最好.若差分步長繼續(xù)減小，由于梯度計算將逐步超出管網(wǎng)模型的計算精度范圍，計算結(jié)果反而變差.

表2 差分步長求解結(jié)果對比Tab.2 Comparison of results with differentδs

（2）梯度計算方法

針對本文研究的算例管網(wǎng)，分別分析了通過差分法（見式4）以及線性比例法（見式6）進行梯度計算的求解效果.線性比例法計算后得到的結(jié)果如圖3所示，其中迭代次數(shù)為400次，收斂時間為1min，梯度向量的模長為3230.83；通過差分法編碼計算后得到的結(jié)果如圖2所示，其中迭代次數(shù)為400次，收斂時間為5min，梯度向量的模長為0.25.

圖3 污染物注入質(zhì)量濃度變化曲線（線性比例法）Fig.3 Injection concentrations of contamination sources with linearly proportional method

圖3表明在候選污染源節(jié)點115，117，119和120處有污染物注入，所定位出的污染源位置在假定污染源注入點附近區(qū)域，但對于污染源節(jié)點污染物濃度的變化識別精度較差；而通過差分法計算（見圖2）則能夠直接確定污染源注入的時間、地點和侵入過程.此外，迭代次數(shù)相同時，采用線性比例法計算結(jié)果的梯度向量模長遠大于采用差分法計算的梯度向量模長.因而采用線性比例法計算的結(jié)果要比用差分步長編碼計算的結(jié)果差，但用線性方式編碼計算的收斂速度快于用差分法計算的速度.在采用線性比例法進行計算時，還需要為迭代變量設(shè)置一個可取值區(qū)間，以保證解的非負性和所能取的最大值，否則將會出現(xiàn)無法收斂的情況.而用差分法計算時，僅需要保證解的非負性，無需設(shè)置解的最大限值.分析其原因在于，線性比例法計算時，由于水質(zhì)監(jiān)測點對于某一候選污染源節(jié)點的權(quán)重系數(shù)，是取相應(yīng)時間步長時段內(nèi)的權(quán)重系數(shù)平均值，因此所求得的梯度向量并不能精確地給出目標(biāo)函數(shù)的收斂方向，導(dǎo)致計算結(jié)果較差.因此，建議在需要快速確定污染源點位時可采用線性比例法，當(dāng)需要更精確地分析污染物注入時間、濃度及變化模式時，采用基于管網(wǎng)水質(zhì)模擬的差分法進行計算.

4 結(jié)論

負梯度法可有效解決供水管網(wǎng)單點污染物注入模式和多點污染物注入模式下的污染源定位問題.相對于遺傳算法等隨機搜索算法，沿負梯度方向搜索具有更高的搜索效率，通過引入差分算法可有效降低目標(biāo)函數(shù)對于決策變量求導(dǎo)的難度.而通過將利用最小覆蓋集確定候選污染源節(jié)點集合的方法與梯度法相結(jié)合，可進一步提高對污染源識別問題的求解能力.此外，通過不同差分步長下的計算結(jié)果表明，當(dāng)其值取水質(zhì)模擬所采用的水質(zhì)公差值時，算法求解結(jié)果最好.而相對于用線性比例法求解梯度向量，差分法可獲得更高的求解質(zhì)量，不但可以定位污染源，并能以較高精度識別污染物注入的時間和濃度變化模式.

由于污染源識別問題是利用有限個水質(zhì)監(jiān)測點獲得的數(shù)據(jù)反演污染源的侵入過程，因此，水質(zhì)監(jiān)測點的數(shù)量及空間布置對于問題的求解有直接的影響.因此，在將來的研究中，本方法應(yīng)進一步與水質(zhì)監(jiān)測點的優(yōu)化布置結(jié)合，以更進一步提高算法的求解效率.此外，對于無法有效獲取管網(wǎng)中水質(zhì)的精確濃度數(shù)據(jù)的情況，如何對突發(fā)污染事件進行污染源的辨識，是未來需要進一步研究的問題.

［1］Laird C D，Biegler L T，van Bloemen Waanders B G，et al.Contamination source determination for water networks［J］.Water Resources Planning and Management，2005，131（2）：125.

［2］Laird C D，Biegler L T，van Bloemen Waanders B G.Mixedinteger approach for obtaining unique solutions in source inversion of water networks［J］.Water Resources Planning and Management，2006，132（4）：242.

［3］GUAN Jiabao，Aral M M，Maslia M L，et al.Identification of contaminant sources in water distribution systems using simulation optimization method：case study ［J］. Water Resources Planning and Management，2006，132（4）：252.

［4］Pries A，Ostfeld A.Contamination source identification in water systems：a hybrid model trees-linear programming scheme［J］.Water Resources Planning and Management，2006，132（4）：263.

［5］Cristo D C，Leopardi A.Pollution source identification of accidental contamination in water distribution networks［J］.Water Resources Planning and Management，2008，134（2）：197.

［6］LIU Li，Ranji Ranjithan S，Mahinthakumar G.Contamination source identification in water distribution systems using an adaptive dynamic optimization procedure［J］.Water Resources Planning and Management，2011，137（2）：183.

［7］YANG Xueyao，Boccelli D L，De Sanctis A E.A Bayesian approach for probabilistic contaminantion source identification［C］／／World Environmental and Water Resource Congress.［S.l.］：ASCE，2011.

［8］Rossman L A.EPANET2 users’manual［M］.Cincinnati：U.S.Environmental Protection Agency，2000.

［9］王鳳仙.供水管網(wǎng)水質(zhì)管理軟件應(yīng)用和預(yù)警技術(shù)研究［D］.上海：同濟大學(xué)環(huán)境科學(xué)與工程學(xué)院，2010.WANG Fengxian.Study on application of water quality management software and early-warning technology for water distribution systems［D］.Shanghai：College of Enviromental Science and Engineering of Tongji University，2010.