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

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

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

關鍵詞：地表水水質；大數據分析；物理模型；好溪流域

中圖分類號：X830.2 文獻標志碼：B

1 研究背景

動態掌控流域地表水中污染物的來源及其成因，是實現流域污染精準治理的基礎。由于污染來源多樣、驅動一遷移轉化過程復雜，以及下墊面因子高度時空變異性等原因，難以直接確定污染負荷。

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

2 數據與方法

2.1 研究區概況

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

2.2 控制斷面污染負荷模型構建

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

2.2.1 直排源負荷

基本單元i內直排源達到控制斷面的污染負荷如式（1）：

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

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

2.2.2 非直排源對河道污染負荷

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

控制斷面匯流區的污染負荷為內基本單元（數量為n）負荷的累加：

2.3 基于大數據解析的參數分區

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

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

3 結果與討論

3.1 參數分區與率定

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

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

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

3.2 主要負荷貢獻區模擬分析

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

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

3.3 與SWAT模型對比分析

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

文章提出的模型與SWAT模型在物理機制上，均為均衡模型，以集總的方式考慮污染物遷移的動態過程對濃度變化的影響。因此，在物理機制上，率定的參數是設定時間步長的最優解。因而模型設定的時間步長發生變化，模擬精度將在一定程度降低。

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

4 結論

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