王忠慧，貢 力,2，康春濤，王 鴻，楊軼群

(1.蘭州交通大學(xué)土木工程學(xué)院，甘肅 蘭州 730070;2.蘭州交通大學(xué)調(diào)水工程及輸水安全研究所，甘肅 蘭州 730070)

突發(fā)水污染事件會(huì)使整個(gè)水生態(tài)環(huán)境遭到破壞，從而造成巨大的損失，甚至?xí)l(fā)社會(huì)動(dòng)蕩，需采取必要的應(yīng)急措施來應(yīng)對(duì)[1-2]。在水污染事件發(fā)生的第一時(shí)間判斷出污染源的位置、掌握其污染物的源強(qiáng)、獲知污染事件發(fā)生的時(shí)間[3]，是處理這種不確定性問題的關(guān)鍵所在，而能否快速有效地找到污染源的位置決定著能否最大限度減小突發(fā)水污染事件影響的范圍[4]。

突發(fā)水污染事情追蹤溯源就是了解水污染發(fā)生以及發(fā)展的全過程，而源項(xiàng)信息的識(shí)別尤其關(guān)鍵，國(guó)內(nèi)外在源項(xiàng)信息識(shí)別方面已取得了很多研究成果，但還沒有形成一個(gè)完整的解決突發(fā)水污染事情溯源問題的體系。源項(xiàng)信息識(shí)別常用的方法有確定性方法、概率統(tǒng)計(jì)方法和耦合的概率密度分析方法3類[5]。在確定性方法研究方面，辛小康等[6]通過遺傳算法與數(shù)學(xué)模型相結(jié)合的方法，對(duì)一維單點(diǎn)源單變量和多點(diǎn)源多變量問題分別進(jìn)行了研究；Boano等[7]利用地質(zhì)統(tǒng)計(jì)法有效恢復(fù)了任意分布源和多個(gè)獨(dú)立點(diǎn)源的水污染釋放歷史；曹宏桂等[8]將水污染問題利用有限元法進(jìn)行正演，用PSO-DE混合優(yōu)化算法結(jié)合移動(dòng)監(jiān)測(cè)平臺(tái)進(jìn)行反演，證明了該優(yōu)化算法在二維水環(huán)境下的適用性。確定性方法最大的特點(diǎn)是利用一組最優(yōu)的污染物濃度信息求解污染源項(xiàng)；缺點(diǎn)是計(jì)算量較大且沒有充分考慮數(shù)據(jù)的不確定性，而且當(dāng)信息不準(zhǔn)確時(shí)往往結(jié)果誤差很大。在概率統(tǒng)計(jì)方法研究方面，陳海洋等[9]在考慮有限寬度河流瞬時(shí)岸邊污染泄漏的情況下，建立了水體污染識(shí)別數(shù)學(xué)模型，并以典型Metropolis算法構(gòu)建馬爾科夫鏈取得后驗(yàn)概率分布；姜繼平等[10]通過Adaptive Metropolis算法對(duì)后驗(yàn)概率密度進(jìn)行采樣，得到了操作參數(shù)推薦值，為貝葉斯推理技術(shù)的應(yīng)用作出了重要貢獻(xiàn)；Wei等[11]通過分析正向模型的不確定性，結(jié)合AM算法對(duì)假想的河道排放污染物源項(xiàng)進(jìn)行了反演。概率統(tǒng)計(jì)方法的應(yīng)用主要是考慮到了突發(fā)水污染事件的不確定性，雖然能得到污染源項(xiàng)的“可能解”，但是計(jì)算時(shí)抽樣過程極為耗時(shí)且結(jié)果對(duì)隨機(jī)變量的分布信息過于依賴。在耦合的概率密度分析方法研究方面，早期的逆向概率密度函數(shù)來自地下水污染的研究，Neupauer等[12]研究表明，在一維和二維環(huán)境下的地下水，可以用逆向概率密度函數(shù)求得源項(xiàng)的位置和時(shí)間信息；程偉平等[13]將水污染釋放過程重構(gòu)，通過對(duì)比傳統(tǒng)的優(yōu)化方法與逆向概率密度的優(yōu)化方法，結(jié)果表明逆向概率密度優(yōu)化方法減少了計(jì)算量，提高了計(jì)算效率；王家彪等[14]提出結(jié)合概率密度分析的微分進(jìn)化法對(duì)突發(fā)水污染事情進(jìn)行溯源研究。耦合的概率密度分析方法通過實(shí)現(xiàn)源項(xiàng)信息之間的解耦，確保在計(jì)算過程中計(jì)算量不會(huì)過大，也能考慮到突發(fā)水污染事件的不確定性，是確定性方法和概率方法的結(jié)合。耦合的概率密度分析方法不僅簡(jiǎn)化了水污染溯源的過程，而且提高了溯源精度，是新一代的溯源研究方法。

本文采用耦合的概率密度分析方法，將突發(fā)水污染溯源的模型簡(jiǎn)化，利用天牛須搜索(beetle antennae search，BAS)算法[14]識(shí)別一維河渠污染物源項(xiàng)，并通過Matlab進(jìn)行模型仿真試驗(yàn)來驗(yàn)證BAS算法的可行性及應(yīng)用前景。

1 數(shù)學(xué)模型

水污染溯源問題就是通過已檢測(cè)到的固定時(shí)刻的質(zhì)量濃度信息確定污染源位置x、排放時(shí)刻t及源強(qiáng)m03個(gè)污染物源項(xiàng)[15]。本文通過水力學(xué)計(jì)算方法，利用正向質(zhì)量濃度概率密度函數(shù)(污染物從污染源傳播到下游某一斷面的概率)與逆向位置概率密度函數(shù)(從觀測(cè)者角度出發(fā)，由觀測(cè)斷面判斷污染物可能來自任意位置的概率)之間的關(guān)系，實(shí)現(xiàn)源項(xiàng)3個(gè)未知量之間的解耦，進(jìn)而通過質(zhì)量濃度與位置概率密度函數(shù)之間的相關(guān)關(guān)系以及質(zhì)量濃度與正向質(zhì)量濃度概率密度函數(shù)之間的關(guān)系，構(gòu)建水力學(xué)模型。

假設(shè)河渠在長(zhǎng)度方向上遠(yuǎn)遠(yuǎn)大于寬度和深度方向，則污染物排入水體后能在短時(shí)間內(nèi)與水體混合均勻，并且水流流速與污染物質(zhì)量濃度均勻分布。在污染物只隨流程方向變化的情況下，污染物輸移數(shù)學(xué)模型可以簡(jiǎn)化為一維的水流水質(zhì)耦合模型[16]，計(jì)算公式為

(1)

式中：ρ(x,t)為污染物在x斷面t時(shí)刻的平均質(zhì)量濃度，mg/L；u為河渠斷面平均流速，m/s；E為渠道縱向離散系數(shù)，m2/s；K為污染物降解系數(shù)，s-1。

若水污染的排放形式為瞬時(shí)排放，則污染源以下河渠中的污染物質(zhì)量濃度隨河長(zhǎng)的變化規(guī)律為

ρ(x,t)=

(2)

式中：M為瞬時(shí)排放的污染物面源強(qiáng)度，即源強(qiáng)，g/m2；x0為污染物排放位置，m；t0為污染物排放時(shí)間，s。

任何形式的水污染在河渠輸移過程中都可以用兩種概率密度函數(shù)來表示：一種表示污染物在河渠中的質(zhì)量濃度分布，為正向質(zhì)量濃度概率密度函數(shù)，即污染物從源頭傳播到下游某一位置的質(zhì)量濃度分布概率；另一種表示從觀測(cè)斷面判斷污染物可能在上游的某個(gè)斷面，為逆向位置概率密度函數(shù)，反映污染物在不同位置的可能性大小。而耦合的概率密度分析方法就是通過對(duì)比發(fā)現(xiàn)正向質(zhì)量濃度概率密度函數(shù)與逆向位置概率密度函數(shù)之間的關(guān)系，實(shí)現(xiàn)污染物排放位置、時(shí)間和源強(qiáng)之間的解耦。

該方法起初用于對(duì)地下水的研究，Neupauer等[17]得到了一維地下水域條件下正向質(zhì)量濃度概率密度與逆向位置概率密度之間的關(guān)系：

(3)

式中：P為t′時(shí)刻從斷面x2推斷污染物的源頭位置x1處的逆向位置概率密度，具有m-1的量綱；t′為逆向計(jì)算時(shí)間點(diǎn)；P′為t時(shí)刻污染物由x1斷面輸移到x2斷面的正向質(zhì)量濃度概率密度，也具有m-1的量綱。

因?yàn)镻與P′具有相同的量綱，分別表示逆向位置概率密度輸送過程與正向質(zhì)量濃度概率密度輸送過程，兩者之間是相互伴隨的過程，所以兩者之間具有高度的耦合性，即下游給定斷面觀測(cè)到的質(zhì)量濃度過程與污染源在某一時(shí)刻出現(xiàn)在某一位置的概率相對(duì)應(yīng)。

由式(2)(3)可推導(dǎo)得出逆向的位置概率密度為

P(x1,x2,t′)=

(4)

式中t2-t′為污染物從x1斷面(源頭位置)到x2斷面所用的時(shí)間。

由式(4)可知在污染源位置和釋放時(shí)間信息已知的情況下，就可以計(jì)算得到逆向位置概率密度P。

2 優(yōu)化模型

2.1 優(yōu)化模型的構(gòu)建

由式(3)可知，逆向位置概率密度P與觀測(cè)的質(zhì)量濃度系列ρ之間線性相關(guān)，其相關(guān)系數(shù)r=1。相關(guān)系數(shù)的表達(dá)式為

(5)

優(yōu)化模型建立的目的就是尋找最優(yōu)解，既然逆向位置概率密度函數(shù)已經(jīng)實(shí)現(xiàn)了污染源排放位置、時(shí)間和源強(qiáng)之間的解耦，便可以先從水污染現(xiàn)場(chǎng)推測(cè)得到的一系列污染物排放位置、時(shí)間，即x′0、t′0中尋優(yōu)，確定最優(yōu)的x′0、t′0為x0和t0，再計(jì)算得到M。

根據(jù)觀測(cè)的質(zhì)量濃度系列ρi，假定的x′0、t′0以及對(duì)應(yīng)時(shí)刻的質(zhì)量濃度信息得到Pi的表達(dá)式，構(gòu)造目標(biāo)函數(shù)如式(6)所示。對(duì)于目標(biāo)函數(shù)式(6)，當(dāng)且僅當(dāng)x′0=x0、t′0=t0時(shí)，相關(guān)系數(shù)r=1，此時(shí)目標(biāo)函數(shù)達(dá)到最優(yōu)狀態(tài)。

min(1-r)

(6)

由式(6)得到污染物排放位置和時(shí)間，則源強(qiáng)M可由下式大致推算：

(7)

然后利用正向質(zhì)量濃度概率密度函數(shù)與質(zhì)量濃度信息之間的關(guān)系，進(jìn)一步構(gòu)建優(yōu)化模型計(jì)算最終源強(qiáng)，其目標(biāo)函數(shù)為

(8)

2.2 優(yōu)化模型的求解

本文使用BAS算法求解一維河渠穩(wěn)態(tài)點(diǎn)源擴(kuò)散模型，該算法是2017年提出的一種高效的智能優(yōu)化算法[18]，它通過模擬天牛搜索食物的過程進(jìn)行優(yōu)化計(jì)算。天牛是甲蟲一族，它最大的特征是擁有兩只比身體還要長(zhǎng)的觸角(須)，當(dāng)天牛覓食的時(shí)候，利用兩須搜索食物氣味，通過食物的氣味(也就是目標(biāo)函數(shù))控制移動(dòng)行為和轉(zhuǎn)向行為來搜索食物，當(dāng)天牛搜索到食物的時(shí)候也就是計(jì)算得到最終的優(yōu)化結(jié)果的時(shí)候。為了使模型的計(jì)算結(jié)果精確，本文將BAS算法中的步長(zhǎng)進(jìn)行改進(jìn)，采用變步長(zhǎng)的BAS算法求解模型，隨著步長(zhǎng)的減小，算法的計(jì)算速度由快變慢，完成從全局到局部尋優(yōu)的一個(gè)過程，可以有效縮短算法的收斂速度，提高算法的計(jì)算精度。具體計(jì)算步驟如下：

步驟1：確定優(yōu)化模型的目標(biāo)函數(shù)(式(6)(8))。

步驟2：以al表示左須位置，ar表示右須位置，o表示質(zhì)心，用d0表示兩須之間的距離。假設(shè)天牛頭朝向任意，因而從天牛右須指向左須的向量朝向也是任意的，由歸一化的隨機(jī)向量l表示：

(9)

式中:rands(m,1)為一隨機(jī)向量，與天牛須的指向有關(guān)；m為維數(shù)，即目標(biāo)函數(shù)未知數(shù)的個(gè)數(shù)。

步驟3：經(jīng)過G次迭代后，天牛左右兩須的位置可以根據(jù)l確定為

(10)

式中：aG為第G次迭代后質(zhì)心的位置；dG為第G次迭代時(shí)天牛兩須的距離。

步驟4：由左右兩須的適應(yīng)值確定天牛下一時(shí)刻的位置坐標(biāo)；

aG=aG-1+slsign[f(al)-f(ar)]

(11)

其中

s=cd0

式中：s為步長(zhǎng)；c為常數(shù)；f(al)為al的適應(yīng)度值，表示向左走；f(ar)為ar的適應(yīng)度值，表示向右走；sign為符號(hào)函數(shù)，若f(al)f(ar)，此時(shí)天牛向右須方向移動(dòng)。每進(jìn)行一次迭代，搜索距離和步長(zhǎng)的變化如下：

dG=η1dG-1+d0

(12)

sG=η2sG-1

(13)

式中η1、η2分別是搜索距離更新系數(shù)與步長(zhǎng)衰減系數(shù)。

步驟5：將步驟4計(jì)算得到的aG代入目標(biāo)函數(shù)式中，直到找到目標(biāo)函數(shù)的最優(yōu)解或達(dá)到最大迭代次數(shù)，結(jié)束迭代。

BAS算法流程如圖1所示。

圖1 BAS算法流程Fig.1 BAS algorithm flow chart

3 模型仿真試驗(yàn)

綜合考慮了南水北調(diào)中線[19]、引大入秦輸水渠道等一些實(shí)際河渠的基本輸水狀況后，選定一段沒有支流匯入或突然涌泄水的、可視為一維恒定流進(jìn)行計(jì)算的矩形斷面的河渠進(jìn)行仿真試驗(yàn)。河渠為漿砌塊石護(hù)面，糙率系數(shù)為0.025，底寬為10 m，水深 4 m，坡降0.002 8%，流速0.36 m/s，河渠縱向離散系數(shù)為2.0 m2/s。假設(shè)上午9:00在A斷面處(假定此處為坐標(biāo)原點(diǎn))瞬時(shí)排入1 000 kg的污染物，為了便于對(duì)比分析污染物的質(zhì)量濃度變化過程以及驗(yàn)證BAS算法的計(jì)算結(jié)果，在下游河道設(shè)置了距A點(diǎn)800 m和1 200 m的A1和A2兩個(gè)監(jiān)測(cè)斷面。根據(jù)污染物擴(kuò)散模型(式(2))，在不考慮降解作用的情況下，分別計(jì)算兩個(gè)斷面的質(zhì)量濃度過程，并在得到的計(jì)算結(jié)果中人為加入10%的觀測(cè)誤差(也就是在污染物擴(kuò)散模型中得到的質(zhì)量濃度數(shù)據(jù)中加入10%以內(nèi)的變幅)，生成污染物質(zhì)量濃度變化過程如表1所示。

為便于說明BAS算法計(jì)算本文優(yōu)化模型的可行性，將得到的兩斷面質(zhì)量濃度信息分別進(jìn)行兩次水污染溯源計(jì)算。選擇Matlab軟件進(jìn)行模型仿真計(jì)算，取初始步長(zhǎng)為1 000，迭代200次，再進(jìn)行50次計(jì)算求取平均值后，得到計(jì)算結(jié)果如表2 所示。

從表2可看出，由A1斷面計(jì)算得到的污染物質(zhì)量為975.72 kg，與實(shí)際值相比誤差為2.43%，計(jì)算得到的排放位置差25.43 m，排放時(shí)間超出7 min；由A2斷面計(jì)算得到的污染物質(zhì)量為964.32 kg，與實(shí)際值的誤差為3.57%，計(jì)算得到的排放位置差56.22 m，排放時(shí)間超出9 min。通過不同觀測(cè)斷面的兩次計(jì)算結(jié)果發(fā)現(xiàn)，由A1斷面的濃度信息計(jì)算得到的源項(xiàng)信息的精度比A2斷面高。因?yàn)樵谟?jì)算過程中沒有考慮降解作用，所以不能從水力學(xué)的角度來分析誤差。考慮到BAS算法在計(jì)算時(shí)對(duì)初始步長(zhǎng)的選取具有一定的隨機(jī)性,并且通過多次計(jì)算發(fā)現(xiàn)當(dāng)初始步長(zhǎng)在大于自變量2/3的范圍內(nèi)取值時(shí)，計(jì)算結(jié)果的精度較高。由于兩個(gè)斷面所求的污染物排放位置和時(shí)間(即自變量)范圍不同，A2斷面的自變量范圍較A1斷面大，因此，初始步長(zhǎng)為 1 000時(shí)，對(duì)于A2斷面來說，得到的溯源精度就沒有A1斷面高。在實(shí)際發(fā)生突發(fā)水污染事情時(shí)，可以通過現(xiàn)場(chǎng)考察或是通過已有的資料分析得到自變量的大致范圍，再計(jì)算源項(xiàng)，便可得到更加精準(zhǔn)的污染物排放位置和時(shí)間。

表1 人工生成的觀測(cè)斷面污染物質(zhì)量濃度Table 1 Mass concentration information in observed section

表2 仿真溯源結(jié)果Table 2 Table of simulation traceability results

為了測(cè)試溯源結(jié)果的準(zhǔn)確性，把由A1斷面計(jì)算得到的源項(xiàng)信息代入污染物擴(kuò)散模型(式(2))中，得到污染物在A1斷面處的計(jì)算質(zhì)量濃度變化過程；將計(jì)算質(zhì)量濃度變化過程與觀測(cè)質(zhì)量濃度變化過程進(jìn)行比較(圖2)可知，在質(zhì)量濃度未達(dá)到峰值之前計(jì)算質(zhì)量濃度較觀測(cè)質(zhì)量濃度高，但質(zhì)量濃度達(dá)到峰值之后，觀測(cè)質(zhì)量濃度系列比計(jì)算質(zhì)量濃度系列大。總體上，兩條曲線的相關(guān)性較高，溯源的計(jì)算結(jié)果具有較高的精度，因此溯源模型重構(gòu)的突發(fā)水污染事件能夠反映實(shí)際水污染全過程。

為說明BAS算法在突發(fā)水污染溯源計(jì)算中的特點(diǎn)，在模型仿真試驗(yàn)中采用粒子群算法(partical swarm optimization，PSO)進(jìn)行了計(jì)算，結(jié)果表明，同樣取50次計(jì)算結(jié)果的平均值，PSO算法的精度比BAS算法高。但進(jìn)一步分析發(fā)現(xiàn)BAS算法在迭代了106次之后趨于穩(wěn)定，而PSO算法迭代了564次后才趨于穩(wěn)定；BAS算法計(jì)算50次后取平均值的計(jì)算精度與PSO算法計(jì)算42次后取平均值的計(jì)算精度相同，但BAS算法所用的總計(jì)算時(shí)間為30 s，而PSO算法為125 s?？梢夿AS算法的單次計(jì)算精度比PSO算法低，但是BAS算法收斂速度快，達(dá)到相同的計(jì)算精度BAS算法所用的計(jì)算時(shí)長(zhǎng)也明顯低于PSO算法。

圖2 A1斷面污物物質(zhì)量濃度過程Fig.2 Concentration process diagram of A1 section

總體來說，BAS算法計(jì)算結(jié)果精度高，收斂速度快，算法簡(jiǎn)單，在將來的基于污染物源項(xiàng)解耦的突發(fā)水污染溯源方面有較好的應(yīng)用前景。但BAS算法在計(jì)算過程中步長(zhǎng)一直在隨著迭代次數(shù)的增加而減小，無論能否進(jìn)一步找到最優(yōu)位置，步長(zhǎng)一直在衰減，缺少一種反饋機(jī)制，這也是BAS算法在水污染溯源研究中有待提高的地方。

4 結(jié) 語

本文利用耦合的概率密度分析方法將一維水力學(xué)條件下的正向的質(zhì)量濃度信息與逆向的位置信息進(jìn)行解耦，構(gòu)建了污染物排放時(shí)間、位置和源強(qiáng)解耦后的水力學(xué)模型。將優(yōu)化后的水力學(xué)模型利用BAS算法進(jìn)行模型仿真計(jì)算，發(fā)現(xiàn)BAS算法可以很好地實(shí)現(xiàn)源項(xiàng)信息的求解；通過不同斷面的溯源計(jì)算以及與PSO算法進(jìn)行比較，結(jié)果表明BAS算法收斂速度快，運(yùn)算量小，但是步長(zhǎng)具有衰減性，不管每一步計(jì)算得到的結(jié)果是否變得更優(yōu)，步長(zhǎng)總會(huì)衰減，但多次求取平均值后可以達(dá)到計(jì)算精度要求。本文論證了BAS算法在突發(fā)水污染溯源領(lǐng)域的可行性，拓寬了BAS算法的應(yīng)用領(lǐng)域，但是對(duì)于復(fù)雜的水力學(xué)模型以及更高維問題的研究還有待深入。