基于相關(guān)系數(shù)優(yōu)化法的河流突發(fā)污染源項(xiàng)識(shí)別

陳媛華,王 鵬,姜繼平,郭 亮 (哈爾濱工業(yè)大學(xué)市政環(huán)境工程學(xué)院,黑龍江 哈爾濱 150090)

基于相關(guān)系數(shù)優(yōu)化法的河流突發(fā)污染源項(xiàng)識(shí)別

陳媛華,王 鵬*,姜繼平,郭 亮 (哈爾濱工業(yè)大學(xué)市政環(huán)境工程學(xué)院,黑龍江 哈爾濱 150090)

基于相關(guān)系數(shù)優(yōu)化法,結(jié)合地表水環(huán)境特征和污染物水質(zhì)過程特征,推導(dǎo)出一維河道中單點(diǎn)源瞬時(shí)排放的源項(xiàng)反演算法,得到了污染源排放特征與河流環(huán)境特征參數(shù)的反演公式.采用假想算例進(jìn)行數(shù)值試驗(yàn),綜合分析了流速信息、污染物衰減、監(jiān)測距離、監(jiān)測數(shù)據(jù)誤差及中間參數(shù)?T選取等因素對反演結(jié)果的影響,確定了該方法的適用條件和最優(yōu)條件的尋找方式.間隔10min進(jìn)行兩次監(jiān)測采樣,若監(jiān)測誤差小于5%,反演結(jié)果的相關(guān)系數(shù)達(dá)到-0.97,污染源位置和排放量反演結(jié)果的相對誤差均小于4%,綜合相對誤差在2%以內(nèi).并且方法具有監(jiān)測布點(diǎn)簡單高效,數(shù)據(jù)需求低,編程簡單等優(yōu)點(diǎn),值得在環(huán)境應(yīng)急管理中進(jìn)行實(shí)際應(yīng)用.

點(diǎn)源污染源；河流突發(fā)污染；源項(xiàng)反演；參數(shù)識(shí)別；相關(guān)系數(shù)優(yōu)化

世界各地都受到化學(xué)品泄漏事故的危害[1],其中地表水環(huán)境污染所占比重及對社會(huì)影響和損失最為顯著.一些突發(fā)性水污染事故很有可能在事故初期不知道污染源的位置以及排放特征,甚至有些污染事故在污染團(tuán)流出管轄范圍后管理者才知道污染的發(fā)生,這對流域水環(huán)境管理十分不利. 因此需要找到有效的手段減少事故危害.上述問題可以描述為,已知污染物時(shí)空分布求解模型參數(shù)、結(jié)構(gòu)和邊界條件等信息,在數(shù)學(xué)上可以歸納為反問題[2-5].反問題廣泛存在于環(huán)境水力學(xué)的各個(gè)研究領(lǐng)域,如廢水排放污染源控制問題,水環(huán)境容量計(jì)算問題,污染物總量控制及分配問題等[6-8].

目前,地表水污染源反演研究越來越受到關(guān)注.如 Boano等[9]應(yīng)用地學(xué)統(tǒng)計(jì)的方法識(shí)別河流中的污染源,張玉超等[10]采用支持向量機(jī)方法反演太湖葉綠素a的濃度分布,Cheng等[11]使用反向位置概率密度函數(shù)法和CCHE2D模型程序?qū)恿魑廴驹催M(jìn)行了識(shí)別,朱嵩等[12-13]采用貝葉斯方法進(jìn)行了污染源識(shí)別和模型參數(shù)反演研究,牟行洋[14]利用微分進(jìn)化算法研究了污染物源項(xiàng)識(shí)別反問題.

相關(guān)系數(shù)優(yōu)化算法在地下水污染源識(shí)別問題研究中曾有應(yīng)用[15],但地下水系統(tǒng)與地表水存在較大差異,此算法用于識(shí)別河流污染源的可行性有待研究.2007年郭建青等[16]嘗試將其引入地表水參數(shù)識(shí)別體系,但尚未考慮污染物質(zhì)在水體中的遷移轉(zhuǎn)化規(guī)律(水質(zhì)過程),并且算例的濃度數(shù)據(jù)直接來源于解析解的計(jì)算結(jié)果.因此極大地限制了算法的有效性和普適性.本研究在此基礎(chǔ)上,結(jié)合河流的環(huán)境特征和污染物水質(zhì)過程,考慮水力-水質(zhì)耦合過程改進(jìn)算法,得到了適用于可降解污染物的一維概化的地表水污染源識(shí)別方法.此算法可以反演出污染源特征(污染物排放量、污染源位置及初始排放時(shí)刻),同時(shí)計(jì)算出擴(kuò)散系數(shù)和水流速,并考察濃度監(jiān)測數(shù)據(jù)存在誤差條件下算法的準(zhǔn)確性和有效性.但本算法只適用于穩(wěn)態(tài)流動(dòng)條件下地表水污染物的遷移轉(zhuǎn)化過程.

1 相關(guān)系數(shù)優(yōu)化反演算法及修正

一維均一地表水瞬時(shí)點(diǎn)源排放流場的解析解[17]在許多教科書中可以找到:

式中:C(x,t)為瞬時(shí)脈沖排污引起的污染物濃度在水體的時(shí)空分布,mg/L;M為瞬時(shí)投放的污染物質(zhì)量,kg;A為河流斷面積,m2;Dx為縱向彌散系數(shù),m2/min; k為污染物一級反應(yīng)動(dòng)力學(xué)衰減速率系數(shù),min-1; u為河流平均速度, m/min; x0為污染源位置,m;t0為污染源排放初始時(shí)刻,min.

同一時(shí)刻監(jiān)測水體中不同斷面 xi污染物的濃度Ci,為簡化計(jì)算將首次監(jiān)測時(shí)刻設(shè)定為0時(shí)刻(t=0),污染物初始排放時(shí)刻先于監(jiān)測時(shí)刻,故為負(fù)值.根據(jù)已知的一系列濃度數(shù)據(jù)求得5個(gè)參數(shù), M、Dx、u、x0和t0(假定水的濃度ρ為1).

基于上述概化的水質(zhì)模型和解析解,對以監(jiān)測點(diǎn)與污染源位置的距離為自變量,污染物濃度對數(shù)值為變量的線性函數(shù)的相關(guān)系數(shù)取極值.求得污染源位置,然后選用另一組監(jiān)測數(shù)據(jù),求得其他參數(shù).具體過程如下:

首先,對式(1)取對數(shù),結(jié)果如下:

其中:

根據(jù)式(2)在直角坐標(biāo)系中繪制濃度對數(shù)值關(guān)于(x–xm)2的曲線,可得到一條以mx為斜率、以bx為截距的直線.但由于污染源位置xm未知,無法確定橫軸上各點(diǎn)的坐標(biāo), 因此尚不能確定此方程的斜率與截距.

對于任意的 xm值,可對式(2)進(jìn)行線性回歸,并且計(jì)算回歸的相關(guān)系數(shù).目標(biāo)就是選取一個(gè)適當(dāng)?shù)膞m值使相關(guān)系數(shù)達(dá)到最小,即最接近-1.式(2)的相關(guān)系數(shù)定義為:

n

式中:nx為監(jiān)測點(diǎn)數(shù);為平均距離,其表達(dá)式為:

關(guān)于xm對相關(guān)系數(shù)最小化:

求導(dǎo),可得:

其中,

可見,根據(jù)式(10)～(15)可以直接估算出污染源的位置(xm).實(shí)現(xiàn)了污染源識(shí)別的第一步.

確定xm后,通過式(4)可以計(jì)算得到Xi,然后,可對式(2)進(jìn)行最小二乘法線性擬合.

可得,

從式(10)、式(16) 和式(17)可得xm,mx和bx.通過更進(jìn)一步的檢驗(yàn)發(fā)現(xiàn)只有 2個(gè)可以獨(dú)立的用于求解的未知參數(shù).因此,其余參數(shù)需要通過附加的條件確定.由于未知參數(shù)不能通過同一時(shí)刻的污染物濃度數(shù)據(jù)確定,其與監(jiān)測點(diǎn)數(shù)目及濃度的時(shí)空分布無關(guān).這些信息一定通過不同時(shí)間監(jiān)測的污染物濃度提供.

為得到污染源的信息,需要估算全部的未知參數(shù).首先,在取樣的時(shí)刻(t=0)改寫式(1):

根據(jù)式(18)可以計(jì)算得到3個(gè)不同組合的參數(shù):M,Dxt0和x0-ut0,其估計(jì)取值如下:

然后為了獨(dú)立地求解出各個(gè)參數(shù),污染物濃度數(shù)據(jù)必須取自一定時(shí)間后(?t),同樣按照式(10)和式(16)求得xm′和mx′,它們的定義是:

最后,根據(jù)上式能夠計(jì)算出其余4個(gè)參數(shù):

此外,根據(jù)式(21)可以計(jì)算 M,得到這些參數(shù)后,就可以實(shí)現(xiàn)反演計(jì)算過程.算法實(shí)際應(yīng)用過程見圖1.

圖1 算法實(shí)際應(yīng)用過程Fig.1 Process of applying the inversion algorithm

2 案例應(yīng)用及反演計(jì)算影響因素分析

2.1 算例應(yīng)用

2.1.1 算例來源 某點(diǎn)源污染源非常規(guī)瞬時(shí)排放進(jìn)河流中一定數(shù)量的化學(xué)物質(zhì) α.事發(fā)一段時(shí)間后,預(yù)警系統(tǒng)進(jìn)行響應(yīng),沿水流方向設(shè)置一系列監(jiān)測點(diǎn)反演其排放歷史,其位置見表 1.將解析解(式1)計(jì)算得到的結(jié)果Ci作為相應(yīng)斷面濃度數(shù)據(jù)的期望,加入高斯分布的白噪聲[11],此處為5%.詳細(xì)“觀測”數(shù)據(jù)見表1.其他參數(shù)如下: M=3000kg; u=18m/min;Dx=60m2/min;k=0.001min-1;t0=-85min (第一次取樣前 85min);x0=-1000m; A=20m2.為確定所有參數(shù),在10min后進(jìn)行第二次采樣.

表1 瞬時(shí)點(diǎn)源算例的“觀測”數(shù)據(jù)Table 1 Observed concentration data for the instantaneous source case

表2 算例污染源與參數(shù)反演結(jié)果Table 2 Estimated and real parameters for the conservative and decay instantaneous source case

圖2 t=0時(shí)刻瞬時(shí)排放點(diǎn)源的預(yù)測濃度數(shù)據(jù)和“觀測”濃度數(shù)據(jù)Fig.2 Predicted and observed concentration for a conservative and a decay solute at t=0

2.1.2 反演結(jié)果 基于修正的相關(guān)系數(shù)優(yōu)化法的反演計(jì)算公式,估算出未知參數(shù)見表 2.同時(shí)表2列出了將上述算例在不考慮污染物降解條件下進(jìn)行反演的計(jì)算結(jié)果.估計(jì)參數(shù)與真實(shí)值接近,反演效果十分理想.圖1為第一次采樣時(shí)刻(t=0),實(shí)測數(shù)據(jù)、反演數(shù)據(jù)與理論數(shù)據(jù)的曲線.

由表2和圖2可見,計(jì)算結(jié)果較接近假想算例中的真實(shí)值,相關(guān)系數(shù)絕對值達(dá)到了0.9747.

2.2 反演計(jì)算影響因素分析

2.2.1 河流流速的影響 實(shí)際研究中,河流的流速可方便獲得.因此,考察水流速度條件已知與否對反演結(jié)果的影響.

由表3可知,速度條件對相關(guān)系數(shù)R沒有影響.但其他參數(shù)的計(jì)算結(jié)果有所差別.本算例中流速已知的反演結(jié)果比未知時(shí)略差.

表3 流速對瞬時(shí)點(diǎn)源算例計(jì)算結(jié)果的影響Table 3 The impact of known velocity and unknown velocity for the instantaneous source case

2.2.2 污染物降解過程的影響 BOD是反映水體中污染物濃度的綜合指標(biāo)之一,其在河流中的降解系數(shù)一般為 0.039～20d-1[18].參照 BOD的降解系數(shù)范圍,考察污染物在河流中遷移轉(zhuǎn)化過程對反演結(jié)果的影響,因此選取化學(xué)物質(zhì) α的降解系數(shù)為0.001min-1,即1.44d-1.由圖3可見,對于可降解的污染物質(zhì)(k=0.001min-1),其在水體中的總質(zhì)量隨距離的增大而減少,因而在同一時(shí)刻同一位置其濃度小于非降解污染物;對于保守物質(zhì)(k=0),則濃度計(jì)算結(jié)果偏高.但從反演結(jié)果(表2)可以看出,兩者的計(jì)算誤差均在5%之內(nèi),其中保守物質(zhì)的計(jì)算結(jié)果較優(yōu)于非保守物質(zhì),而考慮可降解物質(zhì)的算法更適合應(yīng)用于實(shí)際情況.

2.2.3 采樣位置和時(shí)間的影響 構(gòu)建一個(gè)評價(jià)指標(biāo),式(28)用向量的方法計(jì)算相對誤差,用于考察采樣位置和時(shí)間同污染源間的時(shí)空距離對算例中反演結(jié)果的綜合影響.相對誤差式中:為反演求得的參數(shù)向量;I為向量第i個(gè)參數(shù);為真實(shí)值的參數(shù)向量.

由圖4可見,在反演過程中,監(jiān)測位置與時(shí)間是一對耦合的參數(shù),共同影響反演結(jié)果.在二維的時(shí)空坐標(biāo)下,源排放后 65～85min內(nèi)在距污染源800～1000m的范圍內(nèi)進(jìn)行第一次采樣, 反演效果較好.

圖4 監(jiān)測位置和時(shí)間對反演結(jié)果的綜合影響Fig.4 The impact of monitoring location and time on the inversion results

2.2.4 監(jiān)測數(shù)據(jù)誤差及兩次采樣時(shí)間間隔對反演結(jié)果的影響 利用相對誤差考察采樣噪聲和時(shí)間間隔(?t)對算例中反演計(jì)算結(jié)果的綜合影響.

由圖5可見,當(dāng)采樣的濃度噪聲不超過5%,取樣時(shí)間間隔對反演結(jié)果的影響較小;當(dāng)濃度誤差超過 5%時(shí),隨著誤差增大,計(jì)算結(jié)果與真實(shí)值間的偏差越大.

圖5 采樣噪聲和時(shí)間間隔對反演結(jié)果的綜合影響Fig.5 The comprehensive impact of sampling noises and interval on the inversion results

3 結(jié)語

污染物在實(shí)際河流的遷移轉(zhuǎn)化規(guī)律十分復(fù)雜,這對水質(zhì)預(yù)測帶了很大的困難.本文考慮了污染物水質(zhì)特征,給出適于可降解污染物質(zhì)的修正的相關(guān)系數(shù)優(yōu)化反演算法,使本身運(yùn)用于地下水的方法向地表水的特定環(huán)境下開展實(shí)際運(yùn)用推進(jìn)了一大步.本方法更適合實(shí)際情景的應(yīng)用.

假想算例計(jì)算表明,本方法的計(jì)算結(jié)果接近真實(shí)情況,并探討和分析論證其在實(shí)際運(yùn)用時(shí)可能遇到的問題,具有一定的應(yīng)用價(jià)值.

[1] 楊 潔,畢 軍,張海燕,等.中國環(huán)境污染事故發(fā)生與經(jīng)濟(jì)發(fā)展的動(dòng)態(tài)關(guān)系 [J]. 中國環(huán)境科學(xué), 2010,30(4):571-576.

[2] 劉家琦.數(shù)學(xué)物理中的反問題及其數(shù)值解法 [J]. 科學(xué)探索1983,(3):75-83.

[3] 王彥飛.反演問題的計(jì)算方法及其應(yīng)用 [M]. 北京:高等教育出版社, 2005:3-13.

[4] 劉曉東,姚 琪,薛紅琴.環(huán)境水力學(xué)反問題研究進(jìn)展 [J]. 水科學(xué)進(jìn)展, 2009,20(6):885-893.

[5] Atmadja J, Bagtzoglou. A C. State of the art report on mathematical method for groundwater pollution source identification [J]. Environmental Forensics, 2001,(2):205-214.

[6] 閔 濤,周孝德.污染物一維非恒定擴(kuò)散逆過程反問題的數(shù)值求解 [J]. 西安理工大學(xué)學(xué)報(bào), 2003,19(1):1-5.

[7] Kitanidis P K, Vomvoris E G. A geostatistical approach to inverse modeling in ground water modeling and one-dimensional simulations [J]. Water Resource Research, 1983,19(3):677-690.

[8] 姜繼平.水污染危險(xiǎn)區(qū)域判定及其原型系統(tǒng)開發(fā) [D]. 哈爾濱:哈爾濱工業(yè)大學(xué), 2008.

[9] Boano F, Revelli R, Ridolfi L. Source identification in river pollution problems: A geostatistical approach [J]. Water Resource Research, 2005,41(7):W07023.doi:1029/2004WR003754.

[10] 張玉超,錢 新,錢 瑜,等.支持向量機(jī)在太湖葉綠素 a非線性反演中的應(yīng)用 [J]. 中國環(huán)境科學(xué), 2009,29(1):78-83.

[11] Cheng W P, Jia Y. Identification of contaminant point source in surface waters based on backward location probability density function method [J]. Advances in Water Resources, 2010,33(4): 397-410.

[12] 朱 嵩,劉國華,王立忠,等.水動(dòng)力-水質(zhì)耦合模型污染源識(shí)別的貝葉斯方法 [J]. 四川大學(xué)學(xué)報(bào), 2009,41(5):30-35.

[13] 朱 嵩,劉國華,毛欣煒,等.基于貝葉斯推理的標(biāo)準(zhǔn) k-ε湍流模型參數(shù)識(shí)別 [J]. 四川大學(xué)學(xué)報(bào), 2010,42(4):78-82.

[14] 牟行洋.基于微分進(jìn)化算法的污染物源項(xiàng)識(shí)別反問題研究 [J].水動(dòng)力學(xué)研究與進(jìn)展, 2011,26(1):24-30.

[15] Sidauruk P, Cheng A H D, Ouazar D. Groundwater contaminant source and transport parameter identification by correlation coefficient optimization [J]. Ground Water, 1998:36(2):208-214.

[16] 郭建青,李 彥,王洪勝,等.突然排污時(shí)計(jì)算河流水質(zhì)參數(shù)的直線解析法 [J]. 水力發(fā)電學(xué)報(bào), 2007,26(4):61-65.

[17] 鄭 彤.環(huán)境系統(tǒng)數(shù)學(xué)模型 [M]. 北京:化學(xué)工業(yè)出版社, 2003: 43-46.

[18] 謝永明.環(huán)境水質(zhì)模型概論 [M]. 北京:中國科學(xué)技術(shù)出版社,1996:124-130.

Contaminant point source identification of rivers chemical spills based on correlation coefficients optimization method.

CHEN Yuan-hua, WANG Peng*, JIANG Ji-ping, GUO Liang (School of Municipal and Environmental Engineering, Harbin Institute of Technology, Harbin 150090, China). China Environmental Science, 2011,31(11)：1802～1807

A novel inversion algorithm based on an optimization approach for river point pollution sources was developed. Mass transport and kinetics processes of the contaminants in surface waters were combined along with the discharge history. And other relative parameters were deduced under the scenario that singular source instantly discharges degradable and soluble chemicals into one-dimensional rivers. A series of numerical experiments were carried out based on the hypothetic cases to analyze inversion effects associated with ambient river flow rates, contaminant decay rates, monitoring sites setting, sampling data errors and time intervals between two groups of sampling. When the monitoring time interval was less than 10 minutes and sampling data errors were controlled fewer than 5% approximately, the relative errors of pollution source location, total released mass and synthetical relative error are under 4%, 4% and 2%, respectively. Results show that parameters calculated fit well with the real values. In addition, the algorithms had the advantages such as efficient sampling process, minimum data requirement as well as easy programming. It was worthwhile to utilize this method for emergency environmental management practices.

point pollution source；river chemical spill；pollution source inversion；parameter identification；correlation coefficient optimization

X703.1

A

1000-6923(2011)11-1802-06

2011-02-01

國家自然科學(xué)基金資助項(xiàng)目(50821002)

* 責(zé)任作者, 教授, pwang73@hit.edu.cn

致謝：文章的英文部分由美國弗吉尼亞大學(xué)土木與環(huán)境工程系龍梧生教授修改與指導(dǎo),在此表示感謝.

陳媛華(1986-),女,黑龍江哈爾濱人,哈爾濱工業(yè)大學(xué)市政環(huán)境工程學(xué)院碩士研究生,主要從事環(huán)境模擬與預(yù)警研究.發(fā)表論文1篇.