基于Matlab龍江鎘污染系統(tǒng)模型的設(shè)計(jì)與實(shí)現(xiàn)

劉慶龍 金斌全 趙偉 李一卓

摘要：本文以廣西龍江河鎘污染事件的數(shù)據(jù)為依據(jù)，從重金屬鎘傳播特性角度出發(fā)，結(jié)合系統(tǒng)動(dòng)力學(xué)理論，率定了重金屬鎘污染擴(kuò)散系數(shù)，建立了鎘污染一維傳遞模型，探討了在不同流速下鎘在各個(gè)監(jiān)測(cè)點(diǎn)之間的傳播過(guò)程。利用模型預(yù)測(cè)了洛東電站鎘污染物濃度，與實(shí)際污染數(shù)據(jù)對(duì)比其[R2]=0.867，并進(jìn)行可視化展示。

關(guān)鍵詞：水體污染;遷移過(guò)程;數(shù)據(jù)分析;一維水質(zhì)模型;龍江河

中圖分類號(hào)：TP311 ? ? ?文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1009-3044（2021）32-0011-05

Design and Implementation of Cadmium Pollution System Model in Longjiang Based on MATLAB

LIU Qing-long， JIN Bin-quan， ZHAO Wei， LI Yi-zhuo

（College of Engineering and Design， Hunan Normal University， Changsha 410081， China）

Abstract： Based on the data of cadmium pollution event in Longjiang river of Guangxi， the diffusion coefficient of heavy metal cadmium pollution was calibrated from the perspective of transmission characteristics of heavy metal cadmium， combined with the theory of system dynamics. The one-dimensional transfer model of cadmium pollution was established， and the transmission process of cadmium between monitoring points at different flow rates was discussed. The model was used to predict the concentration of cadmium pollutants in Luodong power station， and compared with the actual pollution data， the R2 = 0.867， and the visualization was carried out.

Key words： water pollution; migration process; data analysis; one dimensional water quality model; Longjiang river

水資源是人類賴以生存的基礎(chǔ)。然而隨著經(jīng)濟(jì)的發(fā)展，近年我國(guó)突發(fā)水污染事故頻率增加，其中流域突發(fā)水污染事故約占突發(fā)環(huán)境事故的94.7%[1，2]。根據(jù)中華人民共和國(guó)生態(tài)環(huán)境部在2019年5月份發(fā)布的《2018中國(guó)生態(tài)環(huán)境狀況公報(bào)》顯示，2018年，全國(guó)地表水檢測(cè)的1935個(gè)水質(zhì)斷面（點(diǎn)位）中，Ⅰ～Ⅲ類比例為74.3%;劣Ⅴ類比例為6.9%[3]，如圖1所示。重金屬鎘在自然界中多以化合物的形態(tài)微量存在，我國(guó)規(guī)定居民用水中鎘的含量不得超過(guò)0.02mg/L[4]。由于鎘元素難以被生物體內(nèi)分解，鎘元素濃度高于標(biāo)準(zhǔn)的生活用水會(huì)導(dǎo)致鎘元素在生物體內(nèi)聚集，長(zhǎng)期使用被污染水體、食物等會(huì)嚴(yán)重影響居民的身體健康。水體中重金屬污染物的污染程度取決于重金屬的類型、含量和水體的形態(tài)。隨著重金屬?gòu)U水的排放，進(jìn)入水體的重金屬離子會(huì)發(fā)生吸附、擴(kuò)散、沉降等作用，但重金屬元素本身的衰減極其微弱。重金屬離子能夠通過(guò)水體自身的遷移作用進(jìn)行傳播，并在傳播過(guò)程中不斷被水體中的懸浮物以及河床中的沉積物吸附，從而達(dá)到新的平衡狀態(tài)[5]。相對(duì)于有機(jī)物污染，重金屬鎘的衰減性很弱，在研究中常稱之為保守型污染物質(zhì)或惰性污染物。根據(jù)重金屬污染物在水體中遷移的方向不同，可以分為水平、垂直方向的擴(kuò)散和沿水流方向的遷移[6]。

現(xiàn)在學(xué)界廣泛采用的水體重金屬污染傳播模型，大部分根據(jù)質(zhì)量平衡原理構(gòu)建。這類模型的優(yōu)勢(shì)在于可以分析河床對(duì)重金屬的吸附，能夠有效提高模型的精度，但其弊端在于建模過(guò)程比簡(jiǎn)單模型更復(fù)雜，模型建立所需的時(shí)間較長(zhǎng)。突發(fā)型水污染事件往往在應(yīng)急方面有很高的要求，模擬的精度方面可以通過(guò)后期的應(yīng)急監(jiān)測(cè)來(lái)不斷調(diào)整和完善，一般的水相污染模型可以滿足突發(fā)型水污染事件的要求，故本設(shè)計(jì)討論的污染模型不考慮底泥對(duì)重金屬的吸附和降解作用，主要考慮重金屬污染在水體中的水平輸移過(guò)程。

本文通過(guò)分析廣西龍江河基本水體特征，從重金屬鎘傳播特性角度出發(fā)[7]，對(duì)龍江河2012年水質(zhì)數(shù)據(jù)進(jìn)行處理，利用系統(tǒng)動(dòng)力學(xué)理論建立一維龍江水質(zhì)模型，體現(xiàn)出鎘在各個(gè)監(jiān)測(cè)點(diǎn)之間的傳播過(guò)程，并進(jìn)行可視化展示。

1水污染事故水質(zhì)模擬案例分析

1.1研究背景及概況

廣西龍江河地處黔桂珠江水系典型的亞熱帶喀斯特巖溶區(qū)域，具有獨(dú)特的水質(zhì)特征。廣西河池市作為中國(guó)著名的有色金屬產(chǎn)地，各類有色金屬礦產(chǎn)資源儲(chǔ)備量高居全國(guó)前列。龍江河發(fā)源于拉浪水庫(kù)，截止于紅花電站。龍江水域地圖如圖2所示。

龍江河總水域面積約68km2，多年平均流量為67.2m3/s，河床寬100～200m，龍江最大流量為4240m3/s，最大流速2.7m/s，最大水面寬172m。正常流量147m3/s，正常流速1.03m/s，正常水面寬105m?？菟诹髁?963年僅13.1m3/s[8]。

根據(jù)龍江水域的基本水文數(shù)據(jù)，采樣點(diǎn)附近的水力學(xué)特征參數(shù)如表1所示。

1.2實(shí)驗(yàn)數(shù)據(jù)篩選

由于鎘污染事發(fā)突然，受采樣站點(diǎn)的布置等因素的影響，原始采集數(shù)據(jù)中存在一定的異常值（負(fù)值、異常突變值等），所以需要對(duì)龍江水污染原始數(shù)據(jù)中的異常值進(jìn)行缺失值替換處理，預(yù)處理前后數(shù)據(jù)如圖3所示。

對(duì)龍江水污染數(shù)據(jù)進(jìn)行篩選和處理，經(jīng)處理后的數(shù)據(jù)較處理前的數(shù)據(jù)相比，更能夠清晰、直觀地體現(xiàn)出龍江水污染各個(gè)點(diǎn)鎘濃度峰值及其變化曲線，處理后的數(shù)據(jù)也更加準(zhǔn)確，便于進(jìn)行后續(xù)龍江水污染模型的建立與分析。

2 系統(tǒng)模型簡(jiǎn)介

2.1一維水質(zhì)模型的建立

對(duì)于突發(fā)的水污染事故而言，進(jìn)入水體的一般為守恒污染物，即在自然條件下不會(huì)進(jìn)行主動(dòng)降解的污染物（如重金屬、高分子有機(jī)物等），故一般不考慮污染物的衰減過(guò)程。而污染物的分散作用過(guò)程中，一般情況下河流污染物的彌散過(guò)程影響要遠(yuǎn)大于污染物的分子擴(kuò)散與湍流擴(kuò)散，因此對(duì)于污染物的分散作用，一般只考慮彌散過(guò)程。河流污染物的遷移原理圖如圖4所示。

系統(tǒng)動(dòng)力學(xué)以系統(tǒng)分析為基礎(chǔ)，強(qiáng)調(diào)動(dòng)態(tài)的、復(fù)雜的（非線性、反饋回路、延遲和隨機(jī)性等）系統(tǒng)結(jié)構(gòu)，認(rèn)為系統(tǒng)的行為模式與特征主要根植于系統(tǒng)內(nèi)部的反饋結(jié)構(gòu)與機(jī)制，因此適合于研究復(fù)雜系統(tǒng)的結(jié)構(gòu)、功能與動(dòng)態(tài)行為之間的關(guān)系[10]。

根據(jù)龍江水污染事故的監(jiān)測(cè)以及鎘元素傳播途徑特點(diǎn)分析，可以建立一維水質(zhì)模型如式（1）所示;

[?c?t=Ex?2c?x2?ux?c?x?kc] ? ? ? ? ? ? ? ? （1）

其中，假設(shè)在體積微元處[（x）t]時(shí)刻的污染物質(zhì)量濃度為[c]，[ux]為流速，[Ex]為彌散系數(shù)，污染物的衰減速率系數(shù)為[k]，i和j分別表示當(dāng)前位置和當(dāng)前時(shí)刻，i-1和i+1分別表示上一位置和下一位置，j-1和j+1分別表示上一時(shí)刻和下一時(shí)刻，[?x]表示中間站點(diǎn)到兩邊站點(diǎn)距離的平均值，因此式（1）變形為：

[cj+1i?cji?t=Eicji+1?2cji+cji?1?x2?uicji?cji?1?x?kicji]   （2）

當(dāng)污染物為鎘這類難溶性的守恒污染物時(shí)，設(shè)定式（2）中的[ki=0]，則整理可得彌散系數(shù)如式（3）所示：

[Ei=（cj+1i?cji?t+uicji?cji?1?x）/（cji+1?2cji+cji?1?x2）] ?（3）

2.2彌散系數(shù)求解

將龍江水域水力學(xué)特征數(shù)據(jù)帶入一維水質(zhì)模型。本文取糯米灘、西門崖、露塘斷面三個(gè)測(cè)量點(diǎn)的鎘濃度差值以及流速[ui]為0.28m/s代入公式（3），來(lái)求取各個(gè)點(diǎn)的彌散系數(shù)[Ex]值，如表2所示。

然后通過(guò)多組數(shù)據(jù)求出彌散系數(shù)[Ex]關(guān)于鎘濃度差[Δc]的函數(shù)關(guān)系，最終確定龍江一維水質(zhì)模型參數(shù)[Ex]的非線性擬合，R2=0.5153，如圖5所示。

從圖5可以看出：

l當(dāng)濃度[Δc]為負(fù)值時(shí)，求得彌散系數(shù)[Ex]也為負(fù)值，也就是點(diǎn)落在第三象限;同理，當(dāng)濃度差[Δc]為正值時(shí)，求得彌散系數(shù)[Ex]也為正值，也就是點(diǎn)落在第一象限。彌散系數(shù)[Ex]的絕對(duì)值隨濃度差[Δc]絕對(duì)值的增大而增大，即彌散程度隨濃度差擴(kuò)大而提高;

l當(dāng)濃度差[Δc]為0mg/L時(shí)，彌散系數(shù)[Ex]理論上也為0m2/s;當(dāng)濃度差Δc為0mg/L時(shí)，彌散系數(shù)[Ex]為216.4m2/s，是由于污染物濃度的數(shù)據(jù)測(cè)量存在一定誤差。

l當(dāng)濃度差[Δc]為正值時(shí)，有的點(diǎn)會(huì)取到負(fù)值，即有點(diǎn)落在第二象限;當(dāng)濃度差[Δc]為負(fù)值時(shí)，有的點(diǎn)會(huì)取到正值，即有點(diǎn)落在第四象限。簡(jiǎn)而言之，若有點(diǎn)落在第二象限或者第四象限，原因是每一個(gè)位置到達(dá)峰值的時(shí)間不同，會(huì)導(dǎo)致公式中各個(gè)位置之間的濃度差不同，出現(xiàn)正負(fù)性不一致的情況。

采用二次擬合方法，得出的龍江水質(zhì)模型的彌散系數(shù)：

[Ex=8313x2+8120x+216.4] ? ? ? ? ? ? ? （4）

由于模型受到地理環(huán)境等外界不可抗因素的影響，會(huì)導(dǎo)致河流流速發(fā)生改變，從而影響彌散系數(shù)的大小。根據(jù)實(shí)際的理論分析，可以得出當(dāng)流速[ui]增大時(shí)，會(huì)使得彌散系數(shù)[Ex']增大;當(dāng)流速[ui]減小時(shí)，會(huì)使得彌散系數(shù)[Ex'']減小，如圖6和圖7所示。

采用二次擬合方法，得出流速增大到0.36m/s時(shí)龍江水質(zhì)模型的彌散系數(shù)：

[Ex'=1.88×105x2+8.595×104x+2133] ? （5）

采用二次擬合方法，得出流速減小到0.2m/s時(shí)龍江水質(zhì)模型的彌散系數(shù)：

[Ex''=3850x2+6404x+15.93] ? ? ? ? ? ? （6）

取濃度差[Δc]的值為0.0341、0.0083、-0.0191、-0.0285、-0.1581，通過(guò)其相對(duì)的[Ex]、[Ex']和[Ex'']可以看出彌散系數(shù)的變化，和[Ex]的絕對(duì)值相比，[Ex']的絕對(duì)值是增大的，而[Ex'']的絕對(duì)值是減小的，如表3所示。

表3 ? 不同流速下的彌散系數(shù)變化

[濃度差[Δc]（mg/L） [Ex]（m2/s） [Ex']（m2/s） [Ex'']（m2/s） 0.0341 405.60 794.03 234.24 0.0083 1553.97 1735.30 1372.64 -0.0191 1263.78 1844.07 683.49 -0.0285 1493.81 2537.57 450.06 -0.1581 -4925.68 -7752.26 -2099.10 ]

將模型參數(shù)求解出后，為了驗(yàn)證模型是否能在鎘污染的實(shí)際復(fù)雜情況下保持合理與準(zhǔn)確性，結(jié)合數(shù)據(jù)，分段模擬污染物傳播預(yù)測(cè)濃度，將預(yù)測(cè)實(shí)現(xiàn)可視化展示，并與實(shí)際監(jiān)測(cè)數(shù)據(jù)進(jìn)行比對(duì)。

當(dāng)采用擬合[Ex]的數(shù)據(jù)組時(shí)，污染物的實(shí)際數(shù)值相對(duì)于預(yù)測(cè)值偏低，考慮可能是因?yàn)楹恿鞯啄鄬?duì)于鎘的吸附作用。從整體預(yù)測(cè)效果來(lái)看，污染物濃度的變化趨勢(shì)、變化速率均能較好匹配，數(shù)據(jù)比對(duì)如圖8所示。但是，圖8中某些實(shí)際數(shù)據(jù)會(huì)出現(xiàn)0值，這是因?yàn)閷?shí)際測(cè)量中存在一定測(cè)量誤差，在預(yù)測(cè)中模型能將這些誤差值計(jì)算出來(lái)進(jìn)行補(bǔ)足。

9數(shù)據(jù)比對(duì)（去除0值）

將實(shí)際數(shù)據(jù)中0值點(diǎn)進(jìn)行去除處理，然后再次進(jìn)行數(shù)據(jù)比對(duì)，如圖9所示。計(jì)算實(shí)際數(shù)據(jù)及預(yù)測(cè)數(shù)據(jù)之間差距，其中均方差為：[7.53×10?4]，[R2]為：0.867。

3模型可視化

龍江水污染模型的預(yù)測(cè)系統(tǒng)整體設(shè)計(jì)滿足了系統(tǒng)要求，可視化系統(tǒng)所實(shí)現(xiàn)主要功能模塊包含了用戶登錄，讀取已知數(shù)據(jù)，設(shè)置參數(shù)等模塊。具體實(shí)現(xiàn)形式如流程圖10所示。

系統(tǒng)主要包含了用戶登錄界面和水污染數(shù)值預(yù)測(cè)界面兩大部分，其中用戶登錄界面如圖11所示，水污染數(shù)值預(yù)測(cè)界面如圖12所示。

水污染預(yù)測(cè)系統(tǒng)界面主要通過(guò)了GUI的文本框按鈕，圖形窗，編輯文本，以及靜態(tài)文本來(lái)構(gòu)成。靜態(tài)文本主要用作展示各模塊名稱，編輯文本用于導(dǎo)入數(shù)據(jù)，設(shè)定參數(shù)。點(diǎn)擊函數(shù)回調(diào)按鈕會(huì)將設(shè)定點(diǎn)的預(yù)測(cè)結(jié)果在圖形axes窗口中畫出圖來(lái)，達(dá)到預(yù)測(cè)的目的。

4 結(jié)論

本文以龍江鎘污染的實(shí)測(cè)數(shù)據(jù)來(lái)建立基于系統(tǒng)動(dòng)力學(xué)的鎘污染預(yù)測(cè)模型，求取了[Ex]，并對(duì)[ux]的變化進(jìn)行定性分析，實(shí)現(xiàn)了對(duì)洛東電站下一時(shí)刻鎘污染濃度變化的預(yù)測(cè)，[R2]=0.867，為突發(fā)性的鎘污染提供理論預(yù)測(cè)模型。

參考文獻(xiàn)：

[1] 艾恒雨，劉同威.2000—2011年國(guó)內(nèi)重大突發(fā)性水污染事件統(tǒng)計(jì)分析[J].安全與環(huán)境學(xué)報(bào)，2013，13（4）：284-288.

[2] 劉丹，黃俊.流域突發(fā)水污染事件應(yīng)急能力體系建設(shè)[J].人民長(zhǎng)江，2015，46（19）：71-74.

[3] 2018年《中國(guó)生態(tài)環(huán)境狀況公報(bào)》（摘錄二）[J].環(huán)境保護(hù)，2019，47（12）：50-55.

[4] 林立豐，潘尚霞，王鳳巖，等.基準(zhǔn)劑量法在突發(fā)性鎘污染飲用水應(yīng)急閾值研究中的應(yīng)用[J].環(huán)境與健康雜志，2018，35（4）：295-298.

[5] 宮凱悅.松花江哈爾濱段河流底泥重金屬污染及內(nèi)源釋放規(guī)律研究[D].哈爾濱：哈爾濱工業(yè)大學(xué)，2014.

[6] 趙華林.城市半封閉河道水體重金屬污染分配特征研究[J].四川環(huán)境，2019，38（6）：120-124.

[7] 劉潭秋，沈新平，王漢華.基于一種時(shí)間序列模型的河流重金屬污染濃度預(yù)測(cè)研究[J].計(jì)算技術(shù)與自動(dòng)化，2012，31（3）：29-33.

[8] 童娟.珠江流域概況及水文特性分析[J].水利科技與經(jīng)濟(jì)，2007，13（1）：31-33.

[9] 張永祥，蔡德所，吳攀高.廣西龍江河沉積物重金屬污染特征及生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)[J].中國(guó)農(nóng)村水利水電，2014（6）：50-53.

[10] 王其藩.系統(tǒng)動(dòng)力學(xué)理論與方法的新進(jìn)展[J].系統(tǒng)工程理論方法應(yīng)用，1995（4）：6-12.

【通聯(lián)編輯：唐一東】