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

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

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

基于相關(guān)系數(shù)優(yōu)化法,結(jié)合地表水環(huán)境特征和污染物水質(zhì)過程特征,推導(dǎo)出一維河道中單點源瞬時排放的源項反演算法,得到了污染源排放特征與河流環(huán)境特征參數(shù)的反演公式.采用假想算例進(jìn)行數(shù)值試驗,綜合分析了流速信息、污染物衰減、監(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)測布點簡單高效,數(shù)據(jù)需求低,編程簡單等優(yōu)點,值得在環(huán)境應(yīng)急管理中進(jìn)行實際應(yīng)用.

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

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

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

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

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

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

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

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

基于上述概化的水質(zhì)模型和解析解,對以監(jiā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未知,無法確定橫軸上各點的坐標(biāo), 因此尚不能確定此方程的斜率與截距.

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

n

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

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

求導(dǎo),可得:

其中,

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

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

可得,

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

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

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

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

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

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

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

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

2.1 算例應(yīng)用

2.1.1 算例來源 某點源污染源非常規(guī)瞬時排放進(jìn)河流中一定數(shù)量的化學(xué)物質(zhì) α.事發(fā)一段時間后,預(yù)警系統(tǒng)進(jìn)行響應(yīng),沿水流方向設(shè)置一系列監(jiān)測點反演其排放歷史,其位置見表 1.將解析解(式1)計算得到的結(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ù)據(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時刻瞬時排放點源的預(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)化法的反演計算公式,估算出未知參數(shù)見表 2.同時表2列出了將上述算例在不考慮污染物降解條件下進(jìn)行反演的計算結(jié)果.估計參數(shù)與真實值接近,反演效果十分理想.圖1為第一次采樣時刻(t=0),實測數(shù)據(jù)、反演數(shù)據(jù)與理論數(shù)據(jù)的曲線.

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

2.2 反演計算影響因素分析

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

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

表3 流速對瞬時點源算例計算結(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ì)量隨距離的增大而減少,因而在同一時刻同一位置其濃度小于非降解污染物;對于保守物質(zhì)(k=0),則濃度計算結(jié)果偏高.但從反演結(jié)果(表2)可以看出,兩者的計算誤差均在5%之內(nèi),其中保守物質(zhì)的計算結(jié)果較優(yōu)于非保守物質(zhì),而考慮可降解物質(zhì)的算法更適合應(yīng)用于實際情況.

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

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

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

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

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

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

3 結(jié)語

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

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

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

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

[3] 王彥飛.反演問題的計算方法及其應(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é)報, 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] 姜繼平.水污染危險區(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] 朱 嵩,劉國華,王立忠,等.水動力-水質(zhì)耦合模型污染源識別的貝葉斯方法 [J]. 四川大學(xué)學(xué)報, 2009,41(5):30-35.

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

[14] 牟行洋.基于微分進(jìn)化算法的污染物源項識別反問題研究 [J].水動力學(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] 郭建青,李 彥,王洪勝,等.突然排污時計算河流水質(zhì)參數(shù)的直線解析法 [J]. 水力發(fā)電學(xué)報, 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é)基金資助項目(50821002)

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

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

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