基于對流擴(kuò)散方程的水質(zhì)預(yù)測模型研究

張晨璐， 宋金玲， 康燕， 張經(jīng)武， 樊劉炎

（河北科技師范學(xué)院 數(shù)學(xué)與信息科技學(xué)院 河北省農(nóng)業(yè)數(shù)據(jù)智能感知與應(yīng)用技術(shù)創(chuàng)新中心，河北 秦皇島 066004）

水是人類賴以生存的生命源泉， 隨著工農(nóng)業(yè)的快速發(fā)展， 水污染問題日益突出， 水污染對社會造成了一系列的危害， 不僅對農(nóng)作物和水產(chǎn)造成了巨大影響， 對人類的身體健康也危害頗深。 因此， 為了實現(xiàn)水污染的及時監(jiān)督治理， 對水質(zhì)進(jìn)行預(yù)測預(yù)警研究意義重大。

近年來， 隨著計算機(jī)以及預(yù)測技術(shù)的發(fā)展， 對于水質(zhì)預(yù)測的研究日益多元化， 目前研究成果主要包括基于統(tǒng)計回歸的水質(zhì)預(yù)測模型［1］、 基于ARIMA 的水質(zhì)預(yù)測模型［2－3］、 基于支持向量機(jī)的水質(zhì)預(yù)測模型［4］、 基于神經(jīng)網(wǎng)絡(luò)的水質(zhì)預(yù)測模型［5－7］。但是， 這些研究都是基于特定點(diǎn)位的歷史水質(zhì)數(shù)據(jù)預(yù)測該點(diǎn)位未來一段時間的水質(zhì)變化， 即在時間維度上進(jìn)行水質(zhì)預(yù)測［8－9］。 目前的河流水質(zhì)監(jiān)測中，限于經(jīng)費(fèi)問題， 監(jiān)測設(shè)備的布設(shè)一般比較稀疏， 未布設(shè)監(jiān)測設(shè)備位置的水質(zhì)情況則無法及時掌握， 不利于水污染事件的及時發(fā)現(xiàn)和治理。 因此， 為了實現(xiàn)河流水質(zhì)的細(xì)粒度監(jiān)測， 基于已知的監(jiān)測點(diǎn)位數(shù)據(jù)， 預(yù)測出河流中其他位置的水質(zhì)情況具有重要研究意義。 本研究基于流體力學(xué)的對流擴(kuò)散方程， 對河流上無監(jiān)測設(shè)備位置的水質(zhì)變化情況進(jìn)行空間維度預(yù)測， 以期實現(xiàn)河流水環(huán)境的細(xì)粒度監(jiān)測， 為環(huán)保部門和相關(guān)工作人員的水污染治理提供決策支持。

1 水質(zhì)對流擴(kuò)散方程

根據(jù)水中污染物的時間推進(jìn)和空間上的輸移、擴(kuò)散規(guī)律， 形成了不同的水質(zhì)模型， 并應(yīng)用于水體污染物遷移擴(kuò)散預(yù)測分析［10－12］。 按照系統(tǒng)內(nèi)參數(shù)的空間分布特性， 水質(zhì)模型可以分為一維、 二維和三維模型。 按照水質(zhì)參數(shù)的轉(zhuǎn)移特性， 水質(zhì)模型又分為隨流模型、 擴(kuò)散模型和隨流擴(kuò)散模型（或?qū)α鲾U(kuò)散模型）。 對于河流來說， 如果河流的深度和寬度遠(yuǎn)小于長度， 其在垂向和橫向的污染物濃度梯度可忽略不計， 這種情況下對河流不同點(diǎn)位的污染物濃度計算可以簡化為一維水質(zhì)模型。 且由于污染物在水體的遷移與轉(zhuǎn)化過程中產(chǎn)生了分子的擴(kuò)散， 在這種情況下一般考慮使用一維對流擴(kuò)散方程來模擬污染物的遷移擴(kuò)散， 如式（1）所示為一維對流擴(kuò)散方程。 水質(zhì)預(yù)測模型的建立過程， 是在選定的污染物輸移數(shù)學(xué)模型基礎(chǔ)上， 結(jié)合已知的參數(shù)和初始邊界條件， 從而來預(yù)測水體環(huán)境中污染物時空分布狀況的。

式中： u 為流速， km·h－1； E 為擴(kuò)散系數(shù)，km2·h－1； K 為污染物的衰減系數(shù)， h－1； c 為污染物質(zhì)量濃度， mg·L－1。

式（1）的一維對流擴(kuò)散方程屬于偏微分方程， 目前對一維對流擴(kuò)散方程最常用的數(shù)值解法是有限差分法［13－14］， 有限差分法具有精度高、 涉及網(wǎng)格點(diǎn)少的優(yōu)點(diǎn)。 有限差分法又分為顯式差分法和隱式差分法， 顯式差分根據(jù)差分形式不同又細(xì)分為前向差分、后向差分和中心差分。 而對于整個有限差分法的求解過程則是利用離散的網(wǎng)格節(jié)點(diǎn)來取代連續(xù)的求解區(qū)域， 通過在每個離散的網(wǎng)格節(jié)點(diǎn)上進(jìn)行泰勒展開，將方程中的各階導(dǎo)數(shù)用差分表達(dá)式來近似表示， 再求解差分方程組便可得到偏微分方程的最終解。

差分的精度是以泰勒展開式近似值中誤差項的Δx 階數(shù)為準(zhǔn)， 通過對各種差分的泰勒展開式對比可知， 前向和后向差分均為一階精度、 中心差分是二階精度， 說明中心差分比前向和后向差分的精度更高。 另外， 雖然顯式差分法的形式比較簡單， 但是在很多狀況下計算時間相對較長， 現(xiàn)有文獻(xiàn)顯示隱式差分解比顯式差分解更加穩(wěn)定、 精確、 逼近真解［15］。 綜合上述分析， 本次試驗選用顯式中心差分法和隱式差分法分別對污染物濃度進(jìn)行預(yù)測， 下面對2 種方法的求解精度進(jìn)行驗證。

2 對流擴(kuò)散方程求解及精度驗證

現(xiàn)以某河流的水質(zhì)數(shù)據(jù)對對流擴(kuò)散方程的不同求解方法進(jìn)行驗證， 已知該河流在x ＝0 處有一個排放1 h 的點(diǎn)源， 該河流的u ＝5 km／h， E ＝2 km2／h， K ＝0.151 h－1， 令時間步長Δt ＝0.1 h， 空間步長Δx ＝0.5 km， 且初始時刻的上游邊界x ＝0 處的污染物質(zhì)量濃度為10 mg／L（監(jiān)測站點(diǎn)監(jiān)測的數(shù)據(jù)），下游邊界x ＝10 km 處的污染物質(zhì)量濃度為0， 通過以上條件對每個時間步各節(jié)點(diǎn)的污染物濃度進(jìn)行模擬。 本次試驗將監(jiān)測站點(diǎn)檢測到的污染源排放結(jié)束1.5 h 后部分節(jié)點(diǎn)的污染物（總氮）監(jiān)測值， 和對流擴(kuò)散方程各種解的模擬值進(jìn)行對比， 各節(jié)點(diǎn)的實測值具體如表1 所示。

表1 各監(jiān)測點(diǎn)實測值Tab.1 Measured value of each monitoring point mg·L－1

2.1 顯式差分法求解

根據(jù)求解偏微分方程的近似計算中利用差商代替微商的原理， 式（1）可以近似地用差分表示， 其中對x 的一階偏微分可用前向差分表示為：對x 的二階偏微分可用中心差分表示為：用前向差分代替和， 用中心差分代替后， 將同類項進(jìn)行合并， 式（1）可近似地表示為式（2）的差分方程：

由式（3）可知， 節(jié)點(diǎn)i 的下一時間步的濃度值需由當(dāng)前時間步的i－1， i， i＋1 三個節(jié)點(diǎn)的污染物濃度求得， 因此， 若想求n＋1 時間步N 個節(jié)點(diǎn)的濃度值， 則需要從第1 個時間步開始， 分別利用N個式（3）的方程求出下一時間步的各節(jié)點(diǎn)污染物濃度， 逐步迭代求出（） 后才能求得第n 個時間步各節(jié)點(diǎn)的污染物濃度值。 因此， 為了更方便地對各個時間步下各節(jié)點(diǎn)的污染物濃度值進(jìn)行求解， 可將式（3）表達(dá)成下式的矩陣形式。 根據(jù)式（4）和各已知條件， 按照時間步依次迭代可以求得每個時間步下各節(jié)點(diǎn)的污染物濃度， 從而得到原問題在求解域上的近似值。

根據(jù)上面式（4）的矩陣， 可以編程迭代求解出每個時間步各節(jié)點(diǎn)的污染物濃度， 排污結(jié)束1.5 h后各節(jié)點(diǎn)污染物濃度模擬值如表2 所示。 通過與監(jiān)測站點(diǎn)監(jiān)測到的實測值進(jìn)行對比， 模擬值與實測值的擬合度為0.510 2， 對比發(fā)現(xiàn)從第6 個節(jié)點(diǎn)開始模擬值和實測值相差較大， 原因可能是顯式中心差分法適用于順直天然河流， 而該河流并不是順直河流， 存在地勢落差和彎曲等情況， 導(dǎo)致顯示差分法的求解精度較低。

表2 顯式差分模擬值與實測值對比Tab.2 Comparison of explicit differential analog values and measured valuesmg·L－1

2.2 隱式差分法求解

隱式差分法和顯式差分法中求解偏微分方程的近似計算原理一致， 式（1）可以近似地用差分表示，其中對x 的一階偏微分可用前向差分表示為：＝； 對x 的二階偏微分可用中心差分表示為：用前向差分代替和， 用中心差分代替后， 將同類項進(jìn)行合并，合并同類項最終得到的隱式差分方程如式（5）所示：

式（7）的三角矩陣一般采用Thomas 方法進(jìn)行求解， 它是基于高斯消元法的算法， 算法分為2 個階段： 向前消元和回代。 Thomas 方法可以通過編程實現(xiàn)， 從而計算得到隱式差分法下空間維度各節(jié)點(diǎn)的水質(zhì)數(shù)據(jù)。 排污結(jié)束1.5 h 后各節(jié)點(diǎn)的污染物濃度模擬值如表3 所示， 模擬值與實測值之間的擬合度達(dá)到了0.748 1， 雖然各個監(jiān)測點(diǎn)的模擬值與實測值也存在著一定的誤差， 但是試驗結(jié)果表明， 隨著時間步的減小， 實測值與模擬值之間的誤差也在逐漸縮減， 說明隱式差分法的模擬值更加接近實測值。

表3 隱式差分模擬值與實測值對比Tab.3 Comparison of implicit differential analog values and measured valuesmg·L－1

2.3 解析法求解

解析法通常用來求解非線性方程的解， 能有效地提高準(zhǔn)確度和計算速度。 解析法是從對流擴(kuò)散方程解析出給定條件下物質(zhì)分布的解， 利用解析法可以求解任意初值問題， 方程的解中包含位置和時間。 解析解的求解方法有2 種： 一種是采用不定積分法求解； 另一種是采用漸進(jìn)分析法求解。 不定積分法是個求解微分方程的過程， 而漸進(jìn)分析法則是劃分為幾個漸進(jìn)步驟， 每一步都可以求得近似值，依次不斷迭代， 最終求得近似解。 漸進(jìn)分析法的對流擴(kuò)散方程解析解如式（8）所示。

其中， 初始條件為c（x，0）＝Mδ（x）， 邊界條件為c（±∞，t）＝0， ?c（±∞，t）／?x ＝0。 t 和x 分別表示時間和空間位置， u 表示流速， D 表示擴(kuò)散系數(shù)，M 表示河流排放污染物的總質(zhì)量， 本試驗中排放污染物的總質(zhì)量M ＝15.77 mg（通過已知的各節(jié)點(diǎn)污染物濃度擬合計算求得）。

通過編程計算式（8）求得解析解， 同樣將排污結(jié)束1.5 h 后各節(jié)點(diǎn)的污染物濃度模擬值與實測值進(jìn)行對比， 如表4 所示。 通過對比分析， 實測值與模擬值之間的擬合度達(dá)到了0.918 5， 且相關(guān)系數(shù)為1.000 0， 說明解析解的模擬值與實測值基本上比較接近， 因此有的文獻(xiàn)也將污染物擴(kuò)散的解析解作為實測值使用［16］。

表4 解析解模擬值與實測值對比Tab.4 Comparison of analytical solution analog values and measured valuesmg·L－1

2.4 求解方法對比及分析

通過上面3 種對流擴(kuò)散方程求解方法的對比結(jié)果發(fā)現(xiàn)， 解析解的模擬值與實測值之間的誤差最小， 說明解析解的精度最高、 最接近實測值， 其次是隱式差分法， 最后是顯示差分法。 解析法得到的模擬值與實測值的相關(guān)系數(shù)也最高， 說明解析法比隱式差分法的穩(wěn)定性高， 隱式差分法又比顯式差分法穩(wěn)定性高。 另外， 用解析法求解過程中， 污染物總質(zhì)量M 是一個已知參數(shù)， 而實際情況下很難獲知污染源的排污總量， 因此使用解析法的條件比較苛刻， 不過在污染物總質(zhì)量M 能擬合求解的情況下則可以使用該方法。 隱式差分法的求解過程只需要知道幾個必須的參數(shù)便能進(jìn)行計算， 且計算快速、 計算精度比較高， 因此隱式差分法比較適合實際使用。 綜合3 種求解方法的特點(diǎn)， 在下面的水質(zhì)預(yù)測應(yīng)用中， 由于缺乏各時間步各節(jié)點(diǎn)的實測值，將解析解模擬值作為實測值使用， 采用隱式差分法對各節(jié)點(diǎn)的污染物濃度進(jìn)行預(yù)測。

3 對流擴(kuò)散方程應(yīng)用

在實際河流中應(yīng)用對流擴(kuò)散方程， 采用隱式差分法求解污染物濃度， 與解析法模擬的實測值進(jìn)行對比。 由于上述試驗已經(jīng)證明了隱式差分法比顯示差分法的精確度更高， 更接近實測值， 其可行性更高， 且文獻(xiàn)［16］也表明了隱式差分法的高效可行性， 而對于污染物總質(zhì)量M 的驗證在文獻(xiàn)［17］中已經(jīng)證明了其擬合求解方法的有效性， 利用擬合的污染物總質(zhì)量M 求解解析解可行， 而將解析解模擬值作為實測值使用也是目前一些文獻(xiàn)［16，18］的常用方法， 因此， 利用擬合污染物總質(zhì)量M 求得解析解將其作為實測值， 并與隱式差分法進(jìn)行對比的方法具有一定可靠性。

3.1 研究區(qū)概述

本次試驗以木蘭溪流域為研究區(qū)， 利用對流擴(kuò)散方程進(jìn)行水質(zhì)預(yù)測。 木蘭溪流域位于福建省沿海中部， 地勢自西北向東南沿海傾斜， 上游西北部龍巖峰為最高， 海拔1 267 m。 流域呈扇形， 西部和北部以山地為主， 低山、 峽谷、 盆地錯雜其間， 中部和東部為沖積平原和海積平原。 木蘭溪干流總長105 km， 自源頭至仙游溪口大橋長28.7 km， 溪口大橋至木蘭陂長50.5 km， 平均坡降為1.41‰。

由于整個木蘭溪流域由西北向東南存在高山、盆地、 平原等地勢地貌， 受流域地勢和各支流匯入的影響， 導(dǎo)致不同河段的流速、 污染物擴(kuò)散系數(shù)等參數(shù)會有所變化， 而且該河流的整體走向比較曲折不屬于順直河流， 如果對整個河流采用一個對流擴(kuò)散方程進(jìn)行預(yù)測， 可能導(dǎo)致模擬的污染物濃度出現(xiàn)較大誤差， 甚至不滿足順直河流的變化趨勢而失效。 綜合以上各種因素， 決定按照木蘭溪的地勢及支流匯入情況， 對河流分段建立對流擴(kuò)散模型， 對相應(yīng)河段的污染物濃度進(jìn)行預(yù)測， 以減小模擬值與真實值之間的誤差。

3.2 水質(zhì)預(yù)測模型構(gòu)建

通過分析整個木蘭溪流域， 發(fā)現(xiàn)其干流有7 個主要支流匯入， 各支流匯入點(diǎn)均布設(shè)了監(jiān)測設(shè)備（有水質(zhì)監(jiān)測數(shù)據(jù)）， 考慮到支流可能有污染物匯入干流， 而且河段地勢的不同會導(dǎo)致水文參數(shù)不一樣， 從而影響預(yù)測模型， 因此將支流匯入干流的交叉口作為分段點(diǎn)。 每個河段的起點(diǎn)可以看成一個連續(xù)穩(wěn)定的污染源， 對每個河段分別構(gòu)建對流擴(kuò)散水質(zhì)預(yù)測模型， 根據(jù)每個水質(zhì)預(yù)測模型的隱式差分解， 將上下2 個監(jiān)測點(diǎn)的污染物濃度作為邊界條件代入到差分方程中， 就可以求得每個河段內(nèi)某時間步不同空間步節(jié)點(diǎn)的污染物濃度預(yù)測值。 由于每個河段只是具體的流速、 擴(kuò)散系數(shù)和衰減系數(shù)不同，污染物濃度的模擬過程和計算過程及計算公式完全相同， 因此以某段為例， 利用對流擴(kuò)散方程對污染物的濃度變化進(jìn)行模擬。

已知該段的長度大約為4 km， 該段起點(diǎn)的污染物（總氮）質(zhì)量濃度為10 mg／L（監(jiān)測站點(diǎn)監(jiān)測的數(shù)據(jù)）， 持續(xù)時間為1 h， 因此可將該起點(diǎn)看作排放1 h 的污染源， 該段終點(diǎn)的污染物質(zhì)量濃度為6.8 mg／L（為1 h 后監(jiān)測值）。 由人工測量得到的該段水文參數(shù)分別為： u ＝6 km／h， E ＝3 km2／h、 K ＝0.03 h－1， 取時間步長Δt ＝0.05 h， 空間步長Δx ＝0.5 km。

由于本試驗中的污染源數(shù)據(jù)只有排放濃度， 為了使用解析法求解各時間步各個節(jié)點(diǎn)的污染物濃度， 首先需要擬合計算污染物排放總質(zhì)量M， 已知在污染源排放0.75 h 后該段的各個節(jié)點(diǎn)污染物濃度已經(jīng)監(jiān)測（如圖1 所示）， 將各節(jié)點(diǎn)的監(jiān)測值擬合成一個污染源的排放量函數(shù)， 進(jìn)一步對排放量函數(shù)進(jìn)行定積分求解得到M ＝35.26 mg。 將M ＝35.26 mg和其他已知條件代入式（8）， 可以求得各時間步各節(jié)點(diǎn)的污染物濃度， 并在下文作為實測值使用， 其中將表5 中排污結(jié)束1.5 h 后的各節(jié)點(diǎn)污染物濃度值作為對比數(shù)據(jù)。

圖1 污染源排放量擬合函數(shù)Fig.1 Pollution source emissions fitting function

表5 利用解析解求得的污染物濃度Tab.5 Pollutant concentration obtained by analytical solutionmg·L－1

根據(jù)該段起、 終監(jiān)測站點(diǎn)的污染物濃度， 以及該段的流速、 擴(kuò)散系數(shù)、 衰減系數(shù)、 時間步長、 空間步長等數(shù)據(jù)， 利用隱式差分法計算出該段排污結(jié)束1.5 h 后各節(jié)點(diǎn)的污染物濃度， 模擬結(jié)果如表6所示。

表6 隱式差分法模擬值Tab.6 Analog values of implicit difference decomposition methodmg·L－1

通過對比分析表6 的數(shù)據(jù)發(fā)現(xiàn)， 隱式差分法計算得到的模擬值比較接近實測值， 擬合度較高、 誤差較小， 由于解析解中的污染源排污總質(zhì)量為擬合求得， 與實際排污量本身存在誤差， 而且實際河流中污染物的擴(kuò)散與支流匯入以及河道坡度的升降均有明顯關(guān)系， 解析解模擬值和實測值存在一定誤差屬于正?，F(xiàn)象， 因此利用隱式差分法預(yù)測的各節(jié)點(diǎn)污染物濃度基本符合自然河流的水質(zhì)情況。

4 結(jié)論

本次試驗利用對流擴(kuò)散方程進(jìn)行了空間維度的水質(zhì)預(yù)測， 并將其在木蘭溪流域進(jìn)行了應(yīng)用。

（1） 通過試驗將對流擴(kuò)散方程的顯式差分法、隱式差分法和解析法的求解精度進(jìn)行了對比分析，驗證結(jié)果表明解析法高于隱式差分法、 隱式差分法高于顯式差分法。

（2） 在木蘭溪流域的實際應(yīng)用中， 對木蘭溪干流分段建立了對流擴(kuò)散預(yù)測模型， 并以解析法計算的各時間步各節(jié)點(diǎn)的污染物濃度值作為實測值， 與隱式差分法計算的模擬值進(jìn)行了對比。 通過應(yīng)用可知采用隱式差分法進(jìn)行空間維度預(yù)測具有一定可行性， 但是還需要進(jìn)一步提高預(yù)測精度。

（3） 利用解析法預(yù)測在實際應(yīng)用中雖然存在一定難度， 但是在污染源排污量已知或者可以擬合獲知的情況下， 具有較好的應(yīng)用效果。