王 偉

(上海市政工程設(shè)計(jì)研究總院〈集團(tuán)〉有限公司，上海 200092)

1 研究背景

某地區(qū)發(fā)生突發(fā)性水污染事件，如何在最短時(shí)間內(nèi)確定污染源位置、污染物釋放時(shí)間及強(qiáng)度等信息，并采取相應(yīng)措施進(jìn)行治理已是環(huán)境水力學(xué)中的主要研究課題。引起水污染事故的物質(zhì)分很多種，由于可溶性污染物質(zhì)在溶于水后會(huì)隨水運(yùn)動(dòng)并發(fā)生擴(kuò)散等過程，事故特征不易被發(fā)現(xiàn)，污染物濃度超標(biāo)情況一般只能通過水質(zhì)監(jiān)測(cè)站測(cè)出。因此為了減小發(fā)生該類事故后產(chǎn)生的影響，本文針對(duì)定點(diǎn)排放的可溶性污染物引起的水污染事故展開溯源研究。

早前對(duì)污染源識(shí)別的研究更多的是圍繞地下水開展[1-3]，而目前國內(nèi)外關(guān)于對(duì)污染物在水體中的遷移轉(zhuǎn)化規(guī)律以及相應(yīng)的數(shù)學(xué)模型研究已經(jīng)較為成熟，對(duì)突發(fā)性地表水污染事件溯源問題的研究也有了一些進(jìn)展：如辛小康等[4]采用遺傳算法作為優(yōu)化搜索工具，運(yùn)用一維非恒定水質(zhì)方程解析公式構(gòu)建適應(yīng)度函數(shù)，分別對(duì)單點(diǎn)源和多點(diǎn)源突發(fā)水污染事故進(jìn)行污染源識(shí)別；Cheng等[5]基于伴隨分析法，利用后向位置概率密度函數(shù)和前向概率密度函數(shù)耦合關(guān)系，提出了一個(gè)有效的地表水點(diǎn)源識(shí)別方法；朱嵩等[6]針對(duì)水動(dòng)力-水質(zhì)耦合模型，基于貝葉斯推理，利用馬爾科夫鏈蒙特卡羅抽樣對(duì)污染源強(qiáng)度和質(zhì)量進(jìn)行估計(jì)；楊海東等[7]運(yùn)用了微分蒙特卡羅方法，在Bayesian-MCMC方法的基礎(chǔ)上，結(jié)合微分進(jìn)化的繁殖思想，構(gòu)造Markov鏈識(shí)別污染源信息；汪守東等[8]采用拉格朗日追蹤法模擬油膜的輸移擴(kuò)散，實(shí)現(xiàn)了渤海海峽突發(fā)性溢油事故的模擬和預(yù)報(bào)。

近年來，由于沿海城市工業(yè)生產(chǎn)以及海上運(yùn)輸業(yè)的發(fā)展，河口濱海地區(qū)受到的水污染威脅越來越大[9]。海岸河口因其特殊的地理位置，受到潮汐作用影響，無論來自灣內(nèi)或外海水體中的污染物都會(huì)影響到該河口及其附近區(qū)域的水質(zhì)。雖然現(xiàn)在對(duì)于河口地區(qū)污染物溯源問題的研究并沒有河道污染溯源成熟，但在河口地區(qū)污染物溯源時(shí)可以運(yùn)用河道污染溯源的一些方法。文章結(jié)合了河道可溶性污染物溯源方法，分析河口地區(qū)潮流運(yùn)動(dòng)特點(diǎn)，以ELCIRC水動(dòng)力計(jì)算模塊及拉格朗日追蹤法為基礎(chǔ)，構(gòu)建了基于遺傳智能算法的河口地區(qū)突發(fā)性水污染溯源模型，通過構(gòu)造優(yōu)化模型分別回溯污染源位置及質(zhì)量。

2 二維恒定流水動(dòng)力基本方程及參數(shù)優(yōu)化模型

2.1 二維非恒定流水動(dòng)力模型

對(duì)于河口水污染問題，忽略垂直方向上的物質(zhì)濃度變化及密度流影響，主要考慮河口的橫向和縱向質(zhì)量濃度變化，等深條件下可溶性物質(zhì)的對(duì)流擴(kuò)散輸運(yùn)方程如式(1)[10]。

(1)

其中：H—水深,m；

C—等深的物質(zhì)濃度,g/L；

t—時(shí)間,s；

u，v—x、y方向上速度,m/s；

Kxx，Kxy，Kyx，Kyy—二維擴(kuò)散系數(shù)張量,m2/s。

對(duì)于式(1)，一般采用有限體積法進(jìn)行差分求解，Kong等[11]利用算子分裂方法，將式(1)分成對(duì)流項(xiàng)和擴(kuò)散項(xiàng)兩部分進(jìn)行分步計(jì)算，建立水質(zhì)污染模型求解出非恒定流條件下的物質(zhì)濃度。

在求解河口污染物溯源問題時(shí)，未知參數(shù)的先驗(yàn)信息一般很難確定。文章基于近岸半日潮河口地區(qū)往復(fù)流特性，迅速判斷污染物釋放時(shí)間范圍。利用觀測(cè)濃度數(shù)據(jù)和水質(zhì)污染模型計(jì)算所得數(shù)據(jù)的相關(guān)性，采用改進(jìn)遺傳算法構(gòu)造參數(shù)優(yōu)化模型，快速有效地對(duì)河口地區(qū)污染源參數(shù)進(jìn)行識(shí)別。

2.2 河口地區(qū)半日潮雙峰特性

在河口地區(qū)，由于受到潮汐潮流作用，近岸灣內(nèi)海域以往復(fù)流為主，因此可以通過觀測(cè)點(diǎn)位置污染物濃度隨時(shí)間變化過程線大致判斷出該污染源排放至觀測(cè)點(diǎn)位置時(shí)間。圖1為污染物在投放點(diǎn)及觀測(cè)點(diǎn)之間運(yùn)動(dòng)示意圖，假設(shè)污染物從A點(diǎn)投入，觀測(cè)點(diǎn)為C點(diǎn)，B點(diǎn)為污染物到達(dá)的最遠(yuǎn)位置。如圖1所示，由于受潮汐潮流漲落影響，污染物從投入點(diǎn)隨落潮流運(yùn)動(dòng)到B點(diǎn)后會(huì)隨漲潮流重新回到A點(diǎn)，時(shí)間t1為污染物從A點(diǎn)投入后隨落潮流運(yùn)動(dòng)到C點(diǎn)的時(shí)間，t2為污染物從A點(diǎn)隨落潮流運(yùn)動(dòng)到B點(diǎn)后又隨漲潮流回到C點(diǎn)的時(shí)間。

圖1 污染物運(yùn)動(dòng)示意圖Fig.1 Diagrammatic Sketch of Contaminant Movement

圖2為污染物排放后在觀測(cè)點(diǎn)C觀測(cè)得到的污染物濃度與時(shí)間的關(guān)系，觀測(cè)時(shí)間為24 h，每隔0.5 h觀測(cè)一次，從污染物開始排放時(shí)即開始觀測(cè)。

圖2 觀測(cè)點(diǎn)位置污染物濃度隨時(shí)間變化過程線Fig.2 Contaminant Concentration Variation with Time at Observation Station

由圖1結(jié)合正規(guī)半日潮流性質(zhì)可得式(2)。

(2)

一般情況下，t1、t2的具體時(shí)間是未知的。根據(jù)所測(cè)得的數(shù)據(jù)，結(jié)合圖2，假設(shè)第一個(gè)峰值時(shí)間為t1，第二個(gè)峰值時(shí)間為t2，則t2-t1是可知的，約為6 h。由式(2)可知，t1的具體數(shù)值約為3 h，其與圖2中第一個(gè)峰值所在時(shí)間點(diǎn)一致，即從污染開始排放到第一次運(yùn)動(dòng)至觀測(cè)點(diǎn)位置時(shí)的時(shí)間為3 h，這說明，可以通過觀測(cè)點(diǎn)一系列的觀測(cè)數(shù)據(jù)以及半日潮雙峰特性，大致判別出污染物由源位置首次輸運(yùn)到觀測(cè)點(diǎn)的時(shí)間，即污染物釋放時(shí)間。

2.3 建立二維非恒定流河口溯源優(yōu)化模型

(3)

F(x,t)=min(abs(1-R))

(4)

(5)

通過式(5)構(gòu)成的優(yōu)化模型，可以得到污染源的排放位置x0，y0。

(6)

根據(jù)式(6)可以構(gòu)建優(yōu)化模型，進(jìn)一步推算污染源質(zhì)量。

系數(shù)ωi=1/(ci+1.0)2。

2.4 求解二維非恒定流溯源優(yōu)化模型

上文將河口污染物溯源問題轉(zhuǎn)化為優(yōu)化問題，對(duì)于求解上述優(yōu)化模型，在傳統(tǒng)遺傳算法[15]基礎(chǔ)上，改進(jìn)遺傳算法使其加快收斂速度，提高計(jì)算效率。本文采用“加權(quán)法”結(jié)合遺傳算法中特定步驟的改進(jìn)，重新生成個(gè)體，具體計(jì)算步驟如下。

(1) 初始種群生成、變異、交叉這3個(gè)步驟均與傳統(tǒng)遺傳算法相同。

(7)

(8)

3 算例模擬結(jié)果分析

3.1 改進(jìn)遺傳算法驗(yàn)證

本節(jié)中用于驗(yàn)證的試驗(yàn)渠道信息及算例參數(shù)采集自參考文獻(xiàn)13中的實(shí)例，該實(shí)例為2014年3月22日南水北調(diào)中線的“水質(zhì)突發(fā)污染事件應(yīng)急示范工程”，針對(duì)該實(shí)例已有參數(shù)及觀測(cè)數(shù)據(jù)建立污染物溯源優(yōu)化模型。

本文對(duì)遺傳算法最大的改進(jìn)在于加入了權(quán)重系數(shù)，系數(shù)的大小影響改進(jìn)后遺傳算法的收斂效率。理論上，α越小，收斂速度越快。由表1可知，α在[1/10，1/2)取值時(shí)，收斂效率相差不大，且適應(yīng)度函數(shù)值變化不大。為了保證種群的多樣性，α取值區(qū)間為[1/10，1/3]。該驗(yàn)證算例中選取α值為1/10。

表1 α與n、F(x)的關(guān)系表Tab.1 Relationship of α, n and F(x)

如圖3所示，采用式(9)的均方差來評(píng)價(jià)收斂速度的快慢，均方差系數(shù)越大，種群散度相對(duì)較大。

(9)

其中：Fmin為x0，t0；m0為真實(shí)值時(shí)的適應(yīng)度函數(shù)，此時(shí)Fmin=0。由模擬結(jié)果可知，當(dāng)?shù)螖?shù)小于10時(shí)，改進(jìn)后算法中，適應(yīng)度函數(shù)的方差總是小于傳統(tǒng)算法適應(yīng)度函數(shù)的方差，這表明，相比傳統(tǒng)遺傳算法，改進(jìn)后算法能更快更穩(wěn)定地收斂于真實(shí)值附近，有助于快速確定污染源的相關(guān)參數(shù)。

圖3 適應(yīng)度函數(shù)的方差隨迭代次數(shù)變化情況Fig.3 Relationship of Variance of Fitness Function and Iteration Times

(2)改進(jìn)后的微分進(jìn)化算法收斂速度更快

假設(shè)7∶30時(shí)，t=0，那么9∶00排放污染物，t0的真實(shí)值為5 400 s。由圖4可知：傳統(tǒng)遺傳算法在本例中迭代20次，x0、t0趨于穩(wěn)定；而改進(jìn)后的遺傳算法只需迭代10次，x0、t0便趨于穩(wěn)定，即已收斂。

圖4 迭代過程曲線Fig.4 Iteration Process Curves of Different Algorithms

在實(shí)際問題中，有時(shí)并不需要精確確定污染源的相關(guān)參數(shù)，只需要確定參數(shù)范圍。在這種情況下，改進(jìn)后的遺傳算法更凸顯出優(yōu)勢(shì)。傳統(tǒng)遺傳算法迭代十幾次才能確定的范圍，改進(jìn)后的遺傳算法只需迭代2、3次便能做到。當(dāng)給定的邊界范圍增大時(shí)，改進(jìn)后算法能大大減小收斂時(shí)迭代次數(shù)，準(zhǔn)確地搜索到污染源的位置、時(shí)間和質(zhì)量，提高搜索效率。

3.2 二維河口非恒定流溯源算例

以泉州灣河口污染物溯源為例，泉州灣流場(chǎng)為近岸非恒定往復(fù)流，忽略密度流影響，該算例借鑒文獻(xiàn)[11]中的泉州灣河口模型及計(jì)算方案，采用ELCIRC海洋模型模擬該地區(qū)流場(chǎng)，模擬出的流場(chǎng)應(yīng)用于水質(zhì)模型計(jì)算污染物濃度。對(duì)于泉州灣的復(fù)雜地形，ELCIRC模型采用不規(guī)則三角形網(wǎng)格，利用SMS生成模型所需不規(guī)則三角形網(wǎng)格，三角形網(wǎng)格數(shù)量為31 247個(gè)，節(jié)點(diǎn)數(shù)量為15 934個(gè)。模型開邊界最大潮位位于平均海平面上4.02 m，最小潮位位于平均海平面下2.46 m。該算例水動(dòng)力和水質(zhì)模塊是基于文獻(xiàn)[11]中泉州灣水動(dòng)力及水質(zhì)模型，水動(dòng)力模型計(jì)算出的觀測(cè)點(diǎn)流速與流向與該點(diǎn)實(shí)測(cè)值擬合情況良好，因此文章中對(duì)于模型的可靠性及所取參數(shù)不再贅述。

該地形如圖5所示，位置坐標(biāo)采用平面二維坐標(biāo)，污染源位置位于(9.796，18.717)，上午9點(diǎn)在該污染源位置處傾倒360 kg的污染物。在觀測(cè)點(diǎn)位置(14.405，16.967)，從上午9時(shí)開始監(jiān)測(cè)，每隔0.5 h觀測(cè)一次，連續(xù)觀測(cè)24 h，在上午11∶30開始首次觀測(cè)到污染濃度數(shù)據(jù)，濃度數(shù)據(jù)如圖6所示。污染源和觀測(cè)位置分別在圖5中以十字和空心圓標(biāo)出。

圖5 泉州灣地形及污染源、觀測(cè)點(diǎn)位置示意圖Fig.5 Diagrammatic Sketch of Quanzhou Bay Terrain and Locations of Source and Observation Station

圖6 觀測(cè)點(diǎn)觀測(cè)數(shù)據(jù)Fig.6 Observation Data at Observation Station

3.2.1 算例系數(shù)α率定

假設(shè)已知算例中真實(shí)釋放時(shí)間，對(duì)污染源位置進(jìn)行優(yōu)化，率定改進(jìn)遺傳算法系數(shù)α。權(quán)重系數(shù)α為[0.1,1]，步長(zhǎng)為0.1，分別對(duì)不同α值進(jìn)行溯源計(jì)算，比較達(dá)到相同目標(biāo)函數(shù)值時(shí)所迭代計(jì)算的次數(shù)，為了減小計(jì)算結(jié)果的隨機(jī)性，對(duì)每個(gè)α值重復(fù)溯源計(jì)算10次，取其迭代次數(shù)的平均值，如表2所示。

由表2可知，在α值為0.1時(shí)，其達(dá)到目標(biāo)函數(shù)設(shè)定值的迭代次數(shù)相比于α取其他值時(shí)小得多，說明當(dāng)α取0.1時(shí)可以更好地提高污染物位置溯源的效率，具有很廣的適用性。

3.2.2 算例結(jié)果分析

由圖6觀測(cè)數(shù)據(jù)及灣內(nèi)速度場(chǎng)，通過半日潮雙峰特性估計(jì)的污染釋放時(shí)間為165 min?？紤]到觀測(cè)數(shù)據(jù)峰值可能存在誤差及泉州灣正規(guī)半日潮潮汐，將污染釋放時(shí)間適當(dāng)放大為2～3 h。根據(jù)拉格朗日追蹤法，反算出大致污染源位置分別為(10.024，18.928)、(10.015，19.797)，在圖7中以十字標(biāo)出。根據(jù)該推算位置框定大致溯源的位置，在圖7中以矩形框示出，x為(7.223，11.637)，y為(17.029，20.287)，在該范圍內(nèi)給定初始化種群，個(gè)數(shù)為40，選取適應(yīng)度函數(shù)值較小的前20個(gè)種群，初始位置種群在圖7中以淺色圓圈示出。按照傳統(tǒng)遺傳算法和改進(jìn)遺傳算法思路迭代，分別迭代5次和10次，得到種群位置如圖8(a)～圖8(d)所示。為了更清晰地表示出種群位置，圖8僅截取圖7中矩形框示出，圖8(a)、8(b)最淺色叉形分別表示傳統(tǒng)算法迭代5次和10次后種群位置分布，圖8(c)、8(d)最淺色叉形分別表示改進(jìn)算法迭代5次和10次后種群位置分布。

表2 不同分配系數(shù)下目標(biāo)函數(shù)達(dá)到設(shè)定值(10-1以下)的迭代次數(shù)Tab.2 Iteration Steps of Different Distribution Coefficient When Objective Function Reaching Setting Value (below 10-1)

圖7 溯源示意圖Fig.7 Diagrammatic Sketch of Tracing to the Source

圖8 傳統(tǒng)算法與改進(jìn)算法收斂速度對(duì)比Fig.8 Comparison of Convergence Rate between Traditional Algorithm and Improved Algorithm

由圖8(a)～圖8(d)可知，在傳統(tǒng)和改進(jìn)遺傳算法分別迭代到第5次時(shí)，改進(jìn)后算法中種群收斂速度明顯快于傳統(tǒng)算法，在兩種算法分別迭代到第10次時(shí)，改進(jìn)后算法已經(jīng)收斂，并且收斂位置與源位置基本重合，而傳統(tǒng)算法中種群雖有收斂趨勢(shì)，但較改進(jìn)算法速度慢得多。這說明上文提出的改進(jìn)遺傳算法在求解河口污染物溯源優(yōu)化模型時(shí)可以更快速、更準(zhǔn)確，極大地提高了計(jì)算效率。

得到污染源位置參數(shù)后，可以在此基礎(chǔ)上利用優(yōu)化模型進(jìn)一步推算污染源質(zhì)量，由于在優(yōu)化質(zhì)量參數(shù)時(shí)無需再計(jì)算水質(zhì)模型，計(jì)算時(shí)間較短，改進(jìn)遺傳算法優(yōu)勢(shì)不明顯。為了減小模型隨機(jī)性的影響，重復(fù)計(jì)算20次，計(jì)算結(jié)果如圖9所示。

圖9 污染源質(zhì)量回溯結(jié)果Fig.9 Result of Tracing to the Pollution Source Quality

由圖9可知，在迭代多次后，污染源質(zhì)量計(jì)算值為350～370 kg，和算例設(shè)定源強(qiáng)度360 kg基本一致，收斂穩(wěn)定且效果很好，說明在得知污染源位置參數(shù)后，利用優(yōu)化模型可以快速高效地進(jìn)行污染源質(zhì)量反演。本研究為了突出污染效果，將污染源質(zhì)量參數(shù)放大，實(shí)際情況中污染源質(zhì)量可能會(huì)略小。由于在計(jì)算污染源位置時(shí)利用相關(guān)系數(shù)法，假定的污染源質(zhì)量大小不會(huì)影響到源位置的判斷。

4 結(jié)論

本研究利用數(shù)據(jù)相關(guān)性，構(gòu)建了河口溯源模型，結(jié)合半日潮河口地區(qū)潮流特性，采用改進(jìn)遺傳算法對(duì)模型未知參數(shù)進(jìn)行優(yōu)化，實(shí)現(xiàn)了河口地區(qū)可溶性污染物快速溯源，得到了以下主要結(jié)論。

(1)根據(jù)觀測(cè)點(diǎn)的濃度觀測(cè)數(shù)據(jù)，結(jié)合半日潮河口地區(qū)潮流特性，可以單獨(dú)實(shí)現(xiàn)對(duì)污染物投放時(shí)間的快速估計(jì)。

(2)模型采用相關(guān)系數(shù)法，實(shí)現(xiàn)了污染源位置、質(zhì)量的解耦，解放出污染源質(zhì)量參數(shù)，減小了優(yōu)化模型的計(jì)算量，且污染源質(zhì)量大小對(duì)模型的適用性沒有影響。

(3)通過改進(jìn)遺傳算法選擇步驟，加速了算法的收斂，提高了相同迭代步數(shù)下的計(jì)算精度，這對(duì)于快速定位污染源位置及掌握污染源強(qiáng)度等信息并進(jìn)行應(yīng)急處理具有很現(xiàn)實(shí)的指導(dǎo)意義。

本研究使用的河口可溶性污染物水質(zhì)模型為二維模型，沒有考慮垂向濃度的變化以及密度流，適用于近海河口以往復(fù)流為主的地區(qū)，而外海多以旋轉(zhuǎn)流為主，對(duì)外海海域適用性并不好。另外在發(fā)生實(shí)際污染事故時(shí)，關(guān)于污染源未知參數(shù)的先驗(yàn)范圍可能并不全面，而在遺傳算法中對(duì)先驗(yàn)范圍估計(jì)得過大會(huì)導(dǎo)致尋優(yōu)效率的降低，因此可以進(jìn)一步探討對(duì)未知參數(shù)先驗(yàn)范圍的合理估計(jì)。