基于濃度與流量突變的河流總磷通量估算

關鍵詞：突變分析；濃度與流量關系；總磷；LOADEST模型；潮河

河流水質在氣候和人類活動的驅動下時刻發生著變化，關系著河流生態系統的健康發展。農田徑流、畜禽養殖、生活污水排放、土地利用變化等人類活動導致河流水質的惡化，特別是自然水體中廣泛存在的營養物質元素（如氮和磷），如果過量不僅會對水生生物有害，而且會降低飲用水質量，損害人類健康。了解河流水文及水質特征變化及準確估算物質通量，對流域生態系統保護、水資源管理及可持續發展至關重要。

流域水文水質數據通常會由于洪水、干旱、人為活動或環境中發生的其他變化而出現突變點，探究流域流量及水質濃度的長期變化與突變特征，對估算河流物質通量及評估流域管理效果至關重要。突變分析是分析水資源長期變化的重要工具，通常選取某一水文或水質指標，通過相應的統計分析方法確定其時間序列變化發生的時間節點，常見的方法有：Mann-Kendall（M-K）突變檢驗法、Pettitt突變檢驗法、分段回歸法、Wild二值分割法等。這些方法通常應用于單因子時間序列的突變點分析，較難識別因子間（如濃度與流量）的關系變化。而研究表明河流氮磷濃度不僅受流量影響，且與流量間存在復雜的非線性關系，同時存在閾值效應，還可能存在時間變化特征。為了深入理解河流水文水質特征的變化，需要進一步識別濃度與流量關系的變化特征。

當前研究因子間關系變化的方法主要為貝葉斯模型法，可將時間視為獨立變量來檢測營養物質元素與流量、葉綠素含量等其他因子關系變化的時間。研究者建立了包含徑流量和多種水質因子的貝葉斯統計模型，更好地認識濃度與流量之間潛在關系的變化，理解河流水量、水質對人類活動變化的響應。比如，Alameddine等提出了一種結合溫度變化的貝葉斯突變點一閾值模型，該研究量化了總有機氮濃度與流量關系隨時間的變化，判斷美國Neuse河在實施一系列管理措施后流量變化發生的年份。Pirani等結合其他水質指標（水溫、溶解氧和電導率）的影響，識別河流水質組分時間動態及其對流量閾值效應的響應和潛在的季節性變化。這類模型對于認識水環境內部的復雜過程具有重要的理論意義和現實價值，可為分析濃度與流量關系的突變與閾值特征提供有力工具。

河流物質通量通常由徑流量和相對離散的濃度數據估算，LOADEST模型是估算河流物質通量的最常用方法之一。一些研究表明，LOADEST模型準確性受水質組分類型、采樣頻率、土地利用和流域面積等多個因素影響。而模型模擬結果產生偏差的根本原因可能來源于對濃度與流量關系的擬合性差。一方面，LOADEST模型在模擬通量時，低流量數據會對高流量數據的模擬造成影響，使得高流量數據的模擬偏差較大，特別是磷這類受到徑流沖刷影響較大的水質組分。另一方面，LOADEST模型假設濃度與流量之間的關系在不同時間和不同流量狀態下是固定不變的，并未考慮濃度與流量之間復雜的關系變化，采用統一的模型結構估算通量可能會產生較大的偏差。尤其是在長時間序列的模擬中，當河流受到明顯人為活動干擾，如點源的削減、土地利用改變或土地管理措施等，都可能會導致濃度與流量關系發生實質變化，從而使得物質通量評估效果并不理想，降低了模型的擬合效果和預測能力。

已有一些研究表明，在水質建模及通量模擬時考慮突變點，可以減少模型的時間不變性帶來的偏差，但不同方法、不同參數得到的突變點存在差異，當前仍需要加深對單因子（濃度或流量）突變點以及濃度一流量關系突變點對通量模擬影響的理解。為進一步加深對流量與水質濃度突變特征理解，探究濃度及流量突變點對通量模擬的作用，本研究以潮河流域為研究區，基于M-K突變檢驗法以及貝葉斯突變點模型對流量、總磷（TP）濃度以及二者之間的關系進行突變分析，并結合LOADEST模型估算其通量。潮河流域位于密云水庫上游，畜禽養殖、農業施肥、水土流失是影響其水質的主要原因。21世紀以來，該流域先后實施了多項水土保持和水華防治工程，有效減少了水土和氮磷的流失。長期的人類活動使流域水文水質產生的變化使其更適合于突變點的研究。本研究主要目標：分析流域長期徑流量和TP濃度的變化趨勢和突變特征；構建流域TP濃度與流量的非線性突變響應關系，識別二者關系發生突變的時間，分析濃度一流量關系變化特征；基于時間突變點構建LOADEST模型，評估流域的TP通量，比較不同突變點對通量模擬的影響。

1材料與方法

1.1研究區概況

密云水庫是我國首都居民的主要地表飲用水源。潮河位于密云水庫上游，橫跨北京和河北兩省市，如圖1所示。潮河流域面積4888km2，約占密云水庫上游集水面積的40%。地貌以山地丘陵為主，占總面積80%左右，屬于大陸性季風氣候，多年平均氣溫為8.30℃，多年平均降水量為511mm，汛期（6-9月）占年降水量的75%以上。研究區共有人口約40.5萬，農業人口約占總人口的77%，城鎮人口約占23%。21世紀前該流域大部分地區經濟基礎較差，以農業為主，工業不發達，農田管理方式粗放。畜禽養殖、農業施肥等人類活動以及水土流失是影響潮河水質的主要原因。為了確保北京供水安全以及水庫的運行，1989-2005年在該流域陸續開展了多項水利工程建設并采取了水土保持措施。

1.2分析方法與數據來源

數據分析主要分為三部分：（1）采用線性回歸、局部加權回歸分析法（LOWESS）和M-K突變檢驗法分析流量與濃度的時間趨勢與突變特征。（2）通過改進的貝葉斯突變點模型識別濃度一流量關系變化；運用WinBUGS 14專業軟件包基于馬爾科夫鏈蒙特卡洛（Markov Chain Monte Carlo，MCMC）算法對TP濃度一流量關系突變識別模型進行主要參數估計；此部分將在1.3中詳細介紹。（3）基于以上突變點運用LOADEST模型進行分段擬合，并比較模型改進前后的擬合精度。

線性回歸、LOWESS法是常見的趨勢分析方法，LOWESS法相較于線性回歸更能夠處理與描述數據的局部波動特征。M-K突變檢驗法是一種非參數統計檢驗方法，被廣泛應用于水文氣象變化研究中；該方法可通過正向序列和反向序列的趨勢變化統計量UF和UB的交點判斷是否存在突變點；若UF和UB在95.00%顯著性（P=0.05）的置信區間[-1.96，1.96]內存在交點，則該交點為時間序列突變點。LOADEST模型運用中，選取觀測日期當日的流量和濃度監測數據作為平均流量和平均濃度，平均流量與平均濃度的乘積作為當日通量的觀測值，并與模擬結果進行比較驗證。最終模擬結果通過均方根誤差（Root MeanSquared Error. RMSE）和納什效率系數（Nash-Sut-cliffe Efficiency Coefficient．NSE）進行評價。主要數據分析基于MATLAB程序腳本實現，后續所有數據整理匯總與圖形繪制工作均用Origin軟件完成。

研究中潮河流域監測站（下會站）的流量、水質數據來自密云水庫管理處，其中徑流量數據為逐日尺度，水質數據包括：TP濃度、水溫（℃）、溶解氧（%）和電導率，為逐月尺度，篩除缺失值共262組水質數據。以上實際監測數據時間跨度為1992-2014年，共計23a。

為了探討流量和TP濃度變化成因，對潛在的影響因素進行分析，包括：降水、凈人為磷輸入（Net An-thropogenic Phosphorus Input，NAPI）和歸一化植被指數（Normalized Difference Vegetation Index， NDVI），降水數據來自中國氣象數據網（http：//data．cma．cn/），站點為密云站，時間跨度為1992-2014年；NAPI根據研究區各區縣統計年鑒提供的社會經濟數據計算，時間跨度為1995-2014年；NDVI數據來自于資源環境科學數據注冊與出版系統（http：//www．resdc.cn/DOI），時間跨度為1998-2014年。以上數據均以年為尺度，個別年份缺失數據根據已有數據的線性趨勢估算。

1.3貝葉斯突變點模型

模型定義時間t上觀測到的TP濃度為因變量，而徑流量和其他水質參數被視作為自變量；引入了術語“突變點”來表示時間上的變化點，以“閾值”表示流量范圍的變化。假設TP濃度對數服從正態分布：

模型各參數取值以及不確定范圍確定采用MCMC算法，該方法已廣泛應用于統計模型參數估計及參數不確定性分析等領域。在WinBUGS 14的MC-MC參數估計中，其中一共設置了3條馬爾科夫鏈，每條鏈采樣20000次。各項水質影響因素的懲罰樣條計算參照相關文獻。

2結果與分析

2.1流量與TP濃度趨勢及突變分析

圖2展示了潮河流域流量、TP濃度的趨勢與突變分析結果?？梢钥吹搅饔驈搅髁靠傮w呈下降趨勢（圖2a）。M-K突變分析中統計量UF和UB在95.00%顯著性臨界限的上限和下限區間，即[-1.96，+1.96]內存在的交點，以此判斷流域徑流量在1998年初出現突變（圖2b），突變前后的平均流量分別為7.92m3·s-1和2.86m3·s-1．突變前的平均流量明顯高于突變后的平均流量。

潮河流域TP濃度呈現明顯的下降趨勢。根據線性趨勢結果，TP濃度以平均每年1.38x10-3mg·L-1的速度降低（圖2c）。根據M-K突變分析結果，TP濃度在1993年初和1996年底出現突變（圖2d）。據此，TP濃度時間變化總結為三個階段：下降階段（1992-1993年）、速增階段（1994-1996年）和波動降低階段（1997-2014年）；三個時間段平均濃度分別為0.08、0.06mg·L-1和0.03mg·L-1。在1997-2014年，TP濃度于2004年底出現峰值，但未形成突變；2005年以后，TP濃度基本維持在0.03mg·L-1以下的低濃度水平。

2.2濃度一流量關系突變的識別

表1展示了濃度一流量關系的突變點一閾值模型相關參數。結果顯示，潮河流域TP濃度與流量的非線性關系在2004年12月前后發生突變（圖3）。可據此將TP濃度一流量關系劃分為前后獨立的兩個階段：1992-2004年和2005-2014年。前后兩個階段的流量閾值分別為2.36m3·s-1和9.08m3·s-1。模型中溫度、電導率對TP濃度的影響系數為正，產生正面影響；溶解氧則表現出負面影響。

河流物質濃度一流量關系可以分為兩種模式：稀釋主導模式和流量驅動模式。將大于或等于閾值的流量稱為高流量，小于閾值的流量稱為低流量。突變點與流量閾值將流域出口的流量與對應的TP濃度分為了4組（表1）：突變前一低流量；突變前一高流量；突變后一低流量；突變后一高流量。結果顯示，在突變點前，TP濃度與流量的關系是典型的流量驅動模式。在突變點后，濃度與流量關系發生轉變，在低流量情況下仍是流量驅動模式，但在高流量情況下轉變為稀釋主導模式。

2.3基于突變點的TP通量估算

依據流量、TP濃度以及TP濃度一流量關系的突變分析得到了三類突變點，據此將觀測時間序列劃分相應的子序列，利用LOADEST模型分別估算每個子序列對應的污染物通量方程，并與總時間序列（ 1992-2014年）進行比較（表2），并將不同模型的綜合模擬效果繪制散點圖（圖4），比較不同突變點設置對負荷模擬的影響。根據LOADEST模型統計指標，LOAD-EST模型所模擬的不同時段的TP通量與實測值的吻合程度較高，方程決定系數R2在0.66-0.93之間，模型優化系數均具有顯著性統計學意義（顯著性水平Plt;0.05），說明LOADEST模型能夠滿足潮河流域TP的通量評估。

整體上看，模擬結果表現出分段模型優于未分段模型。其中，基于TP濃度突變點的模型（模型3）表現最佳，使NES由未分段模型（模型1）的0.50提高到0.96；與之表現相近是基于流量突變點的模型（模型2，NES=0.95）；模型2和模型3的線性擬合方程斜率均接近于l，RMSE相比模型l有所降低（圖4b、圖4c和表2）?；赥P濃度一流量關系突變點分段模型的NES雖然沒有模型2和模型3高，但比模型l模擬TP的NES提高了0.20（圖4d），整體的RMSE也有所降低，TP通量的擬合曲線斜率由0.34提高到0.58，增加了0.24。

為了方便比較，根據所有類型的突變點，將各模型的序列進行分段比較，共分為四段（表3，1997年與1998均存在突變點，由于時間較近，為方便統計統一為1997年）。在1992-1996年的兩個階段，分段模型的表現均優于未分段模型，其中以模型3的NES最高。在1997-2004年階段，模型2和模型3效果均較好，NESgt;0.80，但模型4表現不佳；而對于2005-2014階段，模型4的NES最高，對比模型1它使NES由-0.31提高到0.89。

綜合來看，以突變點對河流TP通量進行分段估算，一定程度上有助于提高LOADEST模型模擬河流物質通量的準確性，而不同的突變點類型對模型的優化效果不同。

3討論

3.1流量和TP濃度突變的成因分析

流量和TP濃度整體呈下降趨勢，時間序列突變發生在1992-1998年之間，雖然突變發生的時間節點不一致但是時間相對靠近。相關分析顯示流量及TP濃度Pearson相關系數為0.45 （Plt;0.05），一定程度上說明TP濃度變化與流量關系密切，TP濃度的突變可能與流量的突變存在關聯。研究結合了潛在的影響因素（降水、NAPI和NDVI）進行討論，這些因素的變化趨勢如圖5所示。NAPI作為磷的輸入來源，總體呈上升趨勢，無顯著突變點；與流量及TP濃度呈負相關，Pearson相關系數分別為-0.55 （Plt;0.01）和-0.42（Plt;0.05），說明存在其他因素控制了磷的輸出和變化。ND VI整體也呈上升趨勢，與流量呈反比（Pearson相關系數：-0.54，Plt;0.01），體現了植被覆蓋增加對于水土保持具有積極作用，但由于NDVI不存在突變點，與流量和TP濃度的突變缺少關聯性，其影響仍需進一步討論。

降水是徑流的主要驅動因素。在1992-1998年期間，降水出現了多個突變點（1992年、1993年和1996年），與流量以及TP濃度的突變點的時間較為一致。這期間多年為豐水年，降雨量較大，降雨徑流沖刷使得泥沙進入河道，增大了河流中的TP濃度；而后1999-2002年降水相對較少，多為平水年及枯水年，徑流和TP濃度也相應減少。雖然1989-1993年流域已實施了大量的水土保持措施，以及20世紀70-80年代建成的水利工程也開始運作，但由于1994年及1998年出現暴雨和山洪，水土保持措施的減水作用可能被暴雨所導致的強度產流作用所掩蓋。綜合以上討論，1992-1998年期間徑流量及TP濃度的時間序列突變的發生主要受降水變化的控制。

3.2TP濃度一流量關系突變的成因分析

潮河流域TP濃度一流量關系在2005年前后發生了轉變，突變前為典型的流量驅動模式，突變后在高流量情況下為稀釋主導模式。雖然濃度一流量關系是復雜的非線性關系，但利用線性趨勢依然可以看出，突變前在高、低流量區流量與TP濃度均呈弱正相關（圖6a），TP濃度隨流量的增加而增加；突變后，線性趨勢由弱正相關轉為弱負相關（圖6b），這可能表明超過流量閾值后，流域景觀的磷流失不再隨著流量的增加而增加。這與在英國Hampshire Avon河對磷的研究結果類似，在突變點之前高流量與磷酸鹽增加的關系并不顯著，突變點后較高水平的河流流量似乎與磷酸鹽稀釋有關，即稀釋主導模式。由于這項研究時間跨度較短，研究者認為此現象受到極端暴雨事件的影響。

雖然氣候變化或極端天氣的發生也可能導致濃度一流量關系突變，但從長時間尺度來看，潮河流域TP濃度一流量關系突變更可能與流域水土保持措施有關。潮河流域于1998-2000年啟動退耕還林還草和京津風沙源治理工程水土保持項目，至2005年均屬于實施階段。近20多年來潮河流域自然景觀的恢復使得林草地面積比重明顯增加，1998-2014年流域的NDVI由0.73波動上升至0.83（圖5）。植被覆蓋面積的增加很大程度上降低了流域的徑流量和土壤侵蝕強度，也減少了伴隨泥沙輸出的顆粒態磷。另一方面，城市化以及逐年增加的人口增大了流域內的用水需求，也可能導致流域產流減少。但是通過對比可以發現，除去已經討論的1992-1998年期間，2004年降水量也出現較高的峰值，雖然徑流量較低（年均2.46m3·s-1），但是出現了較高的TP濃度峰值（0.29mg·L-1），這可能是強降雨沖刷導致顆粒態磷的流失；關系突變點后，2010年流域也經歷了較大的降水，但對比1992-1998年間和2004年的數據，2010年流量年均為2.81m3·s-1．TP濃度峰值為0.04mg·L-1，這表明隨著強降雨被沖刷進入河流的TP減少了。在Zhang等研究中也得到了類似的結論，他們發現潮河流域氣候和土地利用變化對年輸沙量減少的貢獻率分別為4.12%和95.88%．土地利用變化對徑流和泥沙減少的貢獻大于氣候變化。潮河流域土地利用變化的主要驅動力為上述工程與管理措施的實施，它們對流域水文、水質特征產生了一定的影響。

3.3突變點對通量模擬的影響

研究比較了基于三類突變點分段對TP通量模擬的效果，并了解到基于突變點識別的分段建模有助于改善LOADEST模型的模擬效果，雖然基于TP濃度突變點的分段模型的效果最佳，但不同類型突變點側重點不同。由于TP濃度與流量呈正相關，基于M-K分析的流量突變點與TP濃度突變點都識別了高流量引起的突變，以此為分段結合LOADEST模型模擬，可以減少低流量數據對高流量的影響，提高高流量數據的模擬精確度。

運用突變點一閾值模型識別濃度一流量關系變化，可以更好地理解流域中的物質運輸過程及動態，也可以較客觀和明確地評價流域所采用的管理政策成功與否和起效時間，可為管理者評估措施實施效果以及制定未來管理計劃提供參考。而基于濃度一流量關系突變點的模型，將突變點一閾值模型與LOAD-EST模型結合，在一定程度上可以克服濃度一流量關系時間不變性這一缺點，提高了突變點后對濃度一流量關系的把握。同時也應該看到，當前突變點一閾值模型仍存在一定缺陷。首先目前的方法中只假設了一個突變點的存在，但研究結果表明考慮多個突變點更能把握濃度與流量的變化，因此未來研究中可以推廣到多個突變點。另一方面，突變點一閾值模型考慮了水溫、電導率以及溶解氧等水質參數對濃度一流量關系的影響，而對降水、NDVI等自然或人為的環境因素考慮較為欠缺。

此外，需要考慮分段模擬可能帶來新的不確定性。因為數據集是可能影響模擬準確度的原因之一，而數據的缺失對LOADEST模型的模擬結果會造成一定的影響。對時間序列進行分段會造成每個子序列的數據量減少，從而可能降低模擬效果。目前，僅將濃度與流量的突變識別與LOADEST模型結合，尚未應用于更多模型，因此對于其他模型的影響效果尚不清楚。因此，在未來的研究中，或許可以將突變點分析與更多模型結合，以及考慮開發涉及多個突變點的模型。

4結論

研究基于Mann-Kendall突變檢驗和貝葉斯突變點模型識別潮河流域流量、TP濃度和TP濃度一流量關系的變化點，并將突變點與LOADEST模型結合，以更準確地模擬河流TP的通量。研究發現：

（1）潮河流域流量、TP濃度整體均呈下降趨勢，流量突變點發生于1998年；TP濃度突變點發生于1993年和1996年。流量、TP濃度突變點的形成時間相近，可能受到降水的控制。

（2）TP濃度一流量的非線性關系在2004年12月前后發生突變。突變點前，TP濃度與流量的關系是典型的流量驅動模式；突變點后二者關系會在高流量情況下轉變為稀釋主導模式。TP濃度一流量關系的突變可能與流域內水土保持措施有關。

（3）基于突變點識別的分段建模有助于改善LOADEST模型的模擬效果，而不同類型突變點各有優勢?；赥P濃度突變點的分段模型的整體效果最佳，基于濃度一流量關系突變點的模型對關系突變后TP通量的模擬效果最佳。