基于Matlab龍江鎘污染系統模型的設計與實現

劉慶龍 金斌全 趙偉 李一卓

摘要：本文以廣西龍江河鎘污染事件的數據為依據，從重金屬鎘傳播特性角度出發，結合系統動力學理論，率定了重金屬鎘污染擴散系數，建立了鎘污染一維傳遞模型，探討了在不同流速下鎘在各個監測點之間的傳播過程。利用模型預測了洛東電站鎘污染物濃度，與實際污染數據對比其[R2]=0.867，并進行可視化展示。

關鍵詞：水體污染;遷移過程;數據分析;一維水質模型;龍江河

中圖分類號：TP311 ? ? ?文獻標識碼：A

文章編號：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

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

現在學界廣泛采用的水體重金屬污染傳播模型，大部分根據質量平衡原理構建。這類模型的優勢在于可以分析河床對重金屬的吸附，能夠有效提高模型的精度，但其弊端在于建模過程比簡單模型更復雜，模型建立所需的時間較長。突發型水污染事件往往在應急方面有很高的要求，模擬的精度方面可以通過后期的應急監測來不斷調整和完善，一般的水相污染模型可以滿足突發型水污染事件的要求，故本設計討論的污染模型不考慮底泥對重金屬的吸附和降解作用，主要考慮重金屬污染在水體中的水平輸移過程。

本文通過分析廣西龍江河基本水體特征，從重金屬鎘傳播特性角度出發[7]，對龍江河2012年水質數據進行處理，利用系統動力學理論建立一維龍江水質模型，體現出鎘在各個監測點之間的傳播過程，并進行可視化展示。

1水污染事故水質模擬案例分析

1.1研究背景及概況

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

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

根據龍江水域的基本水文數據，采樣點附近的水力學特征參數如表1所示。

1.2實驗數據篩選

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

對龍江水污染數據進行篩選和處理，經處理后的數據較處理前的數據相比，更能夠清晰、直觀地體現出龍江水污染各個點鎘濃度峰值及其變化曲線，處理后的數據也更加準確，便于進行后續龍江水污染模型的建立與分析。

2 系統模型簡介

2.1一維水質模型的建立

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

系統動力學以系統分析為基礎，強調動態的、復雜的（非線性、反饋回路、延遲和隨機性等）系統結構，認為系統的行為模式與特征主要根植于系統內部的反饋結構與機制，因此適合于研究復雜系統的結構、功能與動態行為之間的關系[10]。

根據龍江水污染事故的監測以及鎘元素傳播途徑特點分析，可以建立一維水質模型如式（1）所示;

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

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

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

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

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

2.2彌散系數求解

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

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

從圖5可以看出：

l當濃度[Δc]為負值時，求得彌散系數[Ex]也為負值，也就是點落在第三象限;同理，當濃度差[Δc]為正值時，求得彌散系數[Ex]也為正值，也就是點落在第一象限。彌散系數[Ex]的絕對值隨濃度差[Δc]絕對值的增大而增大，即彌散程度隨濃度差擴大而提高;

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

l當濃度差[Δc]為正值時，有的點會取到負值，即有點落在第二象限;當濃度差[Δc]為負值時，有的點會取到正值，即有點落在第四象限。簡而言之，若有點落在第二象限或者第四象限，原因是每一個位置到達峰值的時間不同，會導致公式中各個位置之間的濃度差不同，出現正負性不一致的情況。

采用二次擬合方法，得出的龍江水質模型的彌散系數：

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

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

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

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

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

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

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

表3 ? 不同流速下的彌散系數變化

[濃度差[Δ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 ]

將模型參數求解出后，為了驗證模型是否能在鎘污染的實際復雜情況下保持合理與準確性，結合數據，分段模擬污染物傳播預測濃度，將預測實現可視化展示，并與實際監測數據進行比對。

當采用擬合[Ex]的數據組時，污染物的實際數值相對于預測值偏低，考慮可能是因為河流底泥對于鎘的吸附作用。從整體預測效果來看，污染物濃度的變化趨勢、變化速率均能較好匹配，數據比對如圖8所示。但是，圖8中某些實際數據會出現0值，這是因為實際測量中存在一定測量誤差，在預測中模型能將這些誤差值計算出來進行補足。

9數據比對（去除0值）

將實際數據中0值點進行去除處理，然后再次進行數據比對，如圖9所示。計算實際數據及預測數據之間差距，其中均方差為：[7.53×10?4]，[R2]為：0.867。

3模型可視化

龍江水污染模型的預測系統整體設計滿足了系統要求，可視化系統所實現主要功能模塊包含了用戶登錄，讀取已知數據，設置參數等模塊。具體實現形式如流程圖10所示。

系統主要包含了用戶登錄界面和水污染數值預測界面兩大部分，其中用戶登錄界面如圖11所示，水污染數值預測界面如圖12所示。

水污染預測系統界面主要通過了GUI的文本框按鈕，圖形窗，編輯文本，以及靜態文本來構成。靜態文本主要用作展示各模塊名稱，編輯文本用于導入數據，設定參數。點擊函數回調按鈕會將設定點的預測結果在圖形axes窗口中畫出圖來，達到預測的目的。

4 結論

本文以龍江鎘污染的實測數據來建立基于系統動力學的鎘污染預測模型，求取了[Ex]，并對[ux]的變化進行定性分析，實現了對洛東電站下一時刻鎘污染濃度變化的預測，[R2]=0.867，為突發性的鎘污染提供理論預測模型。

參考文獻：

[1] 艾恒雨，劉同威.2000—2011年國內重大突發性水污染事件統計分析[J].安全與環境學報，2013，13（4）：284-288.

[2] 劉丹，黃俊.流域突發水污染事件應急能力體系建設[J].人民長江，2015，46（19）：71-74.

[3] 2018年《中國生態環境狀況公報》（摘錄二）[J].環境保護，2019，47（12）：50-55.

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

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

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

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

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

[9] 張永祥，蔡德所，吳攀高.廣西龍江河沉積物重金屬污染特征及生態風險評價[J].中國農村水利水電，2014（6）：50-53.

[10] 王其藩.系統動力學理論與方法的新進展[J].系統工程理論方法應用，1995（4）：6-12.

【通聯編輯：唐一東】