基于分布式模型的流域地表水質(zhì)模擬和污染來源識別方法研究

杜承陽　呂森偉　葉振宇　林奕影　章思捷

摘要：針對地表水控制斷面污染物濃度預(yù)測和污染來源識別問題，以河道控制斷面匯流區(qū)為基本單元，發(fā)展了直排源和非直排源負(fù)荷計(jì)算方法，基于大數(shù)據(jù)分析方法識別基本單元的污染物遷移轉(zhuǎn)化主控過程劃分了參數(shù)區(qū)，提出了分布式流域地表水質(zhì)模擬方法和污染來源識別方法。分別采用甌江一級支流好溪流域2018年-2022年水文流量、氨氮（NH3）和高錳酸鹽指數(shù)（CODMn）濃度監(jiān)測數(shù)據(jù)進(jìn)行了參數(shù)率定和模型驗(yàn)證。結(jié)果表明：在89.448%的因子覆蓋率的情況下，將好溪流域分為4類參數(shù)區(qū)，模型能夠有效的模擬流域水質(zhì)變化以及識別主要污染原來區(qū)域。

關(guān)鍵詞：地表水水質(zhì)；大數(shù)據(jù)分析；物理模型；好溪流域

中圖分類號：X830.2 文獻(xiàn)標(biāo)志碼：B

1 研究背景

動(dòng)態(tài)掌控流域地表水中污染物的來源及其成因，是實(shí)現(xiàn)流域污染精準(zhǔn)治理的基礎(chǔ)。由于污染來源多樣、驅(qū)動(dòng)一遷移轉(zhuǎn)化過程復(fù)雜，以及下墊面因子高度時(shí)空變異性等原因，難以直接確定污染負(fù)荷。

目前對流域地表水污染預(yù)測以及成因解析的方法主要可以分為二類，第一類基于分布式流域水文和污染負(fù)荷模型（如SWAT模型，HSPF模型），基于物理方法模擬流域內(nèi)各種污染物的遷移轉(zhuǎn)化過程，確定在斷面的污染物濃度，需要詳細(xì)的水文、氣象、水質(zhì)監(jiān)測等參數(shù)，并且由于河道污染負(fù)荷是多個(gè)過程疊加的結(jié)果，通常難以實(shí)現(xiàn)對地表水污染的來源的精確追溯。第二類模型通過了解流域內(nèi)污染負(fù)荷與各種因子之間的統(tǒng)計(jì)關(guān)系，基于因子的變化預(yù)測斷面水質(zhì)變化，包括LSTM模型，BP模型等。結(jié)合兩種方法的優(yōu)點(diǎn)，實(shí)現(xiàn)污染物遷移轉(zhuǎn)化驅(qū)動(dòng)因子諸多、遷移轉(zhuǎn)化過程復(fù)雜條件下的流域水質(zhì)預(yù)測和快速成因解析，對提高流域污染管控能力具有重要意義。基于以上需求，文章發(fā)展了一種流域地表水水質(zhì)預(yù)測和污染來源快速解析方法，在麗水市甌江一級支流好溪流域進(jìn)行了驗(yàn)證分析。

2 數(shù)據(jù)與方法

2.1 研究區(qū)概況

好溪為甌江一級支流，干流總長為45km。流域面積1305km2。流域地形高程為36-1025m。流域年多年平均降雨量1571mm，土壤主要包括紅壤、黃壤、紫色土、粗骨土和水稻土，分別占流域面積的36.93%、24.37%、1.40%、24.34%和11.91%-流域內(nèi)的主要污染負(fù)荷來源（農(nóng)村生活、城鎮(zhèn)、農(nóng)田等）如圖1（a）所示。好溪流域水系和水文、水質(zhì)監(jiān)測斷面以及匯流區(qū)劃分如圖1（b）所示。流域內(nèi)設(shè)有秋塘水文站，水東大橋等河道控制斷面水質(zhì)一動(dòng)監(jiān)測站。每4hr對氨氮（NH4）、高錳酸鹽指數(shù)（CODMn）等水質(zhì)指標(biāo)進(jìn)行一次監(jiān)測。

2.2 控制斷面污染負(fù)荷模型構(gòu)建

河道中污染負(fù)荷的變化主要受到直排源（工業(yè)點(diǎn)源、城鎮(zhèn)污水處理廠以及其他直接排入河道的污染源強(qiáng)）和非直排源（降雨沖刷形成的面源以及地下水出流本底源強(qiáng)）的影響。

2.2.1 直排源負(fù)荷

基本單元i內(nèi)直排源達(dá)到控制斷面的污染負(fù)荷如式（1）：

f1，i=Qp，icp，iexp（-k1R-k2li/v） 式（1）

其中，Qp，i為單元i內(nèi)全部直排污水排放量，cp，i為直排源污染物的平均濃度（由麗水市縉云縣和蓮都區(qū)排污口監(jiān)測數(shù)據(jù)確定），li為基本單元i到控制斷面的距離，R為降雨量，k2為河道中基礎(chǔ)流量（90%枯水年最枯月平均流量）對應(yīng)的污染物自凈系數(shù)，v為基礎(chǔ)流量對應(yīng)的流速，k1為反映河道中自凈能力隨河道流量增加的系數(shù)。

2.2.2 非直排源對河道污染負(fù)荷

根據(jù)降雨徑流曲線法，降雨形成的污染負(fù)荷表示為式（2）：

控制斷面匯流區(qū)的污染負(fù)荷為內(nèi)基本單元（數(shù)量為n）負(fù)荷的累加：

2.3 基于大數(shù)據(jù)解析的參數(shù)分區(qū)

基于30m精度DEM，將好溪流域劃分為87個(gè)基本單元圖1（b），參數(shù)的初始值基于scs法取值，計(jì)算各單元的污染負(fù)荷以及疊加后在控制斷面形成的通量。各基本單元之間的協(xié)方差矩陣為[r]ij，其中：

式（4）中，xi和xj分別為基本單元i和基本單元j的污染物濃度的時(shí)間序列，Cov（·）為協(xié)方差，σ（·）為標(biāo)準(zhǔn)差。基于矩陣正交求解協(xié)方差矩陣的特征值，特征值大小表征矩陣正交之后所對應(yīng)特征向量對整個(gè)矩陣的貢獻(xiàn)程度，即為主成分影響力度大小的指標(biāo)。若特征值小于1，說明該主成分的貢獻(xiàn)力度不如直接引入原變量的力度大，特征值大于l的情況下，作為獨(dú)立因子。累積貢獻(xiàn)率即因子對原始變量的整體貢獻(xiàn)程度。

3 結(jié)果與討論

3.1 參數(shù)分區(qū)與率定

流域87個(gè)基本單元污染負(fù)荷協(xié)方差矩陣特征值分析結(jié)果如表1所示，結(jié)果表明，特征值大于1.0的條件下，可將流域分為4類匯流區(qū)的情況下，全部因子對監(jiān)測信息的累積貢獻(xiàn)為89.448%。

在此基礎(chǔ)上，將流域分為4類參數(shù)區(qū)，采用控制斷面2018年-2021年水質(zhì)監(jiān)測數(shù)據(jù)參數(shù)率定。用Nash-Sutcliffe系數(shù)，相對均方根誤差對模擬誤差進(jìn)行評價(jià)。參數(shù)率定結(jié)果如表2所示。采用率定參數(shù)，以實(shí)測氣象和直排污染源為驅(qū)動(dòng)因子，對2022年1月-10月的流域斷面的污染物濃度變化和貢獻(xiàn)率進(jìn)行了模擬。如圖2（a）、圖2（b）所示模型驗(yàn)證期的降雨量（縉云縣），好溪出口斷面（水東大橋）NH3和CODMn模擬和實(shí)測結(jié)果的比較。采用率定參數(shù)，模型驗(yàn)證期NH3和CODMn模擬和實(shí)測均值的偏差分別為9.23%和7.20%，模擬值和實(shí)測值標(biāo)準(zhǔn)差的平均偏差分比為15.35%和2.34%，表明模擬值的變化區(qū)間與實(shí)測值一致。NH3和CODMn實(shí)測值的均值±標(biāo)準(zhǔn)差分別為0.114+0.070mg/l和2.226+0.727mg/l，模擬值的均值±標(biāo)準(zhǔn)差分別為0.124±0.046mg/l和1.934±0.727mg/l，N_S效率系數(shù)介于0.523-0.622，相對均方根誤差介于0.054-0.168之間。表明所采用的方法模擬地表水質(zhì)具有良好的精度。

NH3的模擬濃度與實(shí)測值的偏差小于CODMn的偏差，以溶解態(tài)形式為主的NH3的模擬精度更好，是由于所采用方法對于各種主控過程均采用了集總式的方法進(jìn)行計(jì)算，在CODMn包括了更多形態(tài)和遷移方式的情況下，集總式的方法不可能避免的產(chǎn)生更大的誤差。

3.2 主要負(fù)荷貢獻(xiàn)區(qū)模擬分析

如圖3（a）-圖3（c）所示水東大橋斷面日降雨量分別為0、19.4mm和30.4mm條件下的NH3主要來源區(qū)域的分布。在降雨量為0的情況下，主要污染負(fù)荷區(qū)與直排源強(qiáng)分布一致。隨著降雨量的增強(qiáng)，首先是Ⅲ類參數(shù)區(qū)的貢獻(xiàn)率增加，其實(shí)是IV類參數(shù)區(qū)的貢獻(xiàn)率增加，然而NH3負(fù)荷的分擔(dān)比例則隨著降雨量的增加而下降，Ⅰ類參數(shù)區(qū)，對NH3負(fù)荷的影響非常小。

圖3（d）為19.4mm降雨量條件下，CODMn高負(fù)荷區(qū)的比較，與圖3（c）比較可以看出，NH3和CODMn的主要污染分布區(qū)域顯然不同，CODMn主要來源于Ⅱ類匯流區(qū)，并且CODMn貢獻(xiàn)率分布的集中程度也明顯的超過NH3貢獻(xiàn)率分布的集中程度。

3.3 與SWAT模型對比分析

分布式流域水文和污染負(fù)荷模型，如SWAT模型，描述了陸面一河道完整水文和污染物遷移轉(zhuǎn)化過程，模型本身則包括了大量的參數(shù)。文章所提出的方法，采用大數(shù)據(jù)方法，通過主控因子識別，保留89.448%的污染成因的條件下，將流域部分分為4類水文區(qū)，每一類水文區(qū)包括7個(gè)參數(shù)，描述了污染物達(dá)到控制斷面整個(gè)遷移轉(zhuǎn)化過程，在保留主控物理機(jī)制的情況下，極大的削減了模型參數(shù)。相比SWAT等分布式模型，文章所采用提出的方法在快速識別主要污染區(qū)域方面具有更大的靈活性。

文章提出的模型與SWAT模型在物理機(jī)制上，均為均衡模型，以集總的方式考慮污染物遷移的動(dòng)態(tài)過程對濃度變化的影響。因此，在物理機(jī)制上，率定的參數(shù)是設(shè)定時(shí)間步長的最優(yōu)解。因而模型設(shè)定的時(shí)間步長發(fā)生變化，模擬精度將在一定程度降低。

SWAT模型等物理模型的誤差成因主要在于：模型各種物理過程的概化（如采用線性過程對非線性過程進(jìn)行近似）和誤差的傳遞性（如SWAT模型首先模擬流量和泥沙過程，水文模擬過程誤差的傳遞與水質(zhì)模擬本身模擬誤差的疊加造成了水質(zhì)模擬誤差顯著的超過了水文過程的模擬誤差）。在大部分流域，盡管水文流量的N＿S系數(shù)能夠達(dá)到0.7，但水質(zhì)指標(biāo)的模擬誤差通常很難超過0.6。文章所發(fā)展的模型，誤差產(chǎn)生的原因主要是兩個(gè)因素，影響控制斷面水質(zhì)各種因子總的覆蓋度為89.448%，顯著降低模型參數(shù)數(shù)量的情況下，同時(shí)也有10.5%的因子覆蓋度并沒有考慮，一定程度造成了計(jì)算偏差。此外，相比SWAT模型，文章對于各種污染物的物理化學(xué)和生物過程所采用的是集總式的方法，陸面水文過程和河道水文過程中污染物各種轉(zhuǎn)化過程的速度顯著不同，在時(shí)間步長較小的情況下，集總式的方法相對更為準(zhǔn)確的反映了各種物理、化學(xué)和生物過程疊加對污染物濃度產(chǎn)生的影響，文章中模型的時(shí)間步長為4hr，較好的模擬了污染物的動(dòng)態(tài)過程也表明了這一現(xiàn)象。

4 結(jié)論

發(fā)展了一種流域地表水質(zhì)預(yù)測和主要負(fù)荷來源區(qū)識別的方法，通過大數(shù)據(jù)進(jìn)行主控因子識別和參數(shù)分區(qū)，顯著的減少了分布式模型的參數(shù)數(shù)量。在浙江省麗水市甌江一級支流好溪流域應(yīng)用結(jié)果表明，該方法經(jīng)過參數(shù)率定后，有效的識別流域污染負(fù)荷及其來源。NH3和CODMn濃度模擬值的系統(tǒng)性偏差小于10%，N_S系數(shù)超過0.5。此方法直接建立了污染源一徑流條件一污染負(fù)荷的關(guān)系，從而實(shí)現(xiàn)了污染負(fù)荷貢獻(xiàn)率的直接核算。