基于隨機森林的引江濟太入湖磷通量計算

陳年浩，陸 昊，饒文昕，劉 彤，錢 新

(南京大學環境學院污染控制與資源化研究國家重點實驗室，南京 210023)

前 言

跨流域調水作為控制湖泊富營養化問題的有力措施，在水質改善和富營養化控制過程中發揮積極作用[1-2]的同時，會挾帶引水區及流經通道的營養鹽和污染物等直接進入受水湖泊[3]，進而對流域水生態環境產生復雜影響。望虞河引江濟太工程是太湖流域綜合治理骨干工程之一，通過引長江水入太湖優化流域水資源配置。當前關于引江濟太工程對受水區貢湖及太湖全湖的水質影響有較多研究[4～6]，調水是否造成太湖總磷濃度反彈尚有爭議[7]，但調水過程中引入了一定量的磷進入太湖已被證實。磷是導致水體富營養化的重要元素之一[8]，相比于顆粒態磷，溶解態磷更易被生物利用[9]，耿雪等[10]指出應以溶解態總磷和溶解態無機磷來考量太湖表層水體磷的生物有效性?，F有研究多關注總磷，鮮有關注不同形態的磷的遷移，當前太湖藍藻暴發的風險依然較高，在未來規劃引水水量增加的情況下，摸清引江濟太不同形態的入湖磷通量具有重要現實意義。

一直以來，自動站駐測和人工巡測均只監測總磷而不區分磷形態，導致不同形態的磷通量無從計算。如何利用現有監測項目的歷史數據反演不同形態的磷成為關鍵，傳統的統計模型難以充分解釋磷形態與總磷及其他環境因子的復雜非線性關系，隨著大數據科學和人工智能的迅猛發展，機器學習的理論和方法在環境科學領域已得到廣泛應用，其中隨機森林(Random forest，RF)簡單易實現、計算開銷小，在很多現實任務中展現出強大的性能[11]，尤其適用于小樣本、高維數據集。在水環境領域中，隨機森林已成功應用于湖庫的葉綠素a模擬與預測[12]及地下水硝酸鹽污染的預測[13]、流域懸浮泥沙的評估[14]等。

本研究以引江濟太工程的望虞河-太湖為研究對象，利用現有監測項目，包括水質自動監測站的常規監測項目和氣象站的觀測項目的歷史數據，通過隨機森林模型建立一種溶解態磷反演方法，并根據模型結果反演多年來望虞河引江濟太的入湖溶解態磷和顆粒態磷通量，以期填補歷史監測數據空缺，為制定兼顧太湖水質控制和生態改善的引江濟太優化方案提供參考和支撐。

1 材料與方法

1.1 研究區域概況

望虞河位于江蘇省無錫市和蘇州市境內(31°27′N ～ 31°47′N，120°25′E ～ 120°51′E)，南起貢湖灣沙墩口，在耿涇口入長江，總長60.8 km，是連接太湖和長江距離最短的流域性河道，引江濟太工程在望虞河入江和入湖處分別通過常熟樞紐和望亭樞紐工程進行水利調度。望虞河流域屬典型平原河網區，兩岸支流口門眾多，東岸支流已全部閘控，西岸仍有口門開敞。西岸地區人口密集、經濟發達，污染負荷排放壓力較大，對引水水質有不利影響[15]。

1.2 樣品采集與實驗分析

在望虞河-太湖布設17個采樣點(圖1)，其中WY1和WY9分別與常熟樞紐和望亭立交樞紐重合，并位于樞紐閘內，引水入湖時WY9可代表望虞河入湖水質，TH2位于錫東水源地。分別于2021年引水期、排水期和關閘期的不同工況下開展現場環境監測，其中引水期2次，排水期1次，關閘期3次，共計6次。

圖1 研究區域及采樣點分布示意Fig.1 Study area and sampling sites

為保證模型在模擬與反演溶解態磷濃度時具有相同的輸入變量，每個采樣點通過JFEAAQ-177多參數水質儀原位監測pH、水溫(WT)、溶解氧(DO)、電導率(EC)、濁度(FTU)，通過1.5 L采水器采集水面下0.3 ～ 0.5 m水樣裝于潤洗過的500 mL PET瓶中，水樣4℃冷藏保存，運回實驗室。采用高錳酸鉀氧化分光光度法測定高錳酸鹽指數(CODMn)，納氏試劑分光光度法測定氨氮(NH3-N)，堿性過硫酸鉀紫外分光光度法測定總氮(TN)，鉬酸銨分光光度法測定總磷(TP)和溶解態總磷(DTP)，顆粒態磷(PP)為TP與DTP之差。同步收集采樣期間的氣象數據，包括大氣壓(AP)、風向(WD)、風速(WS)和相對濕度(RH)。

1.3 隨機森林算法

隨機森林由Breiman[16]在2001年正式提出，在以CART(Classification and Regression Tree)決策樹為基學習器構建Bagging[17]集成的基礎上，融合Ho[18]的隨機子空間(random subspace)思想，進一步在決策樹的訓練過程中引入隨機屬性選擇，關于隨機森林的算法原理本文不再贅述。在Bagging的采樣過程中，未出現在采樣數據集的數據稱為袋外(Out-of-bag，OOB)數據，通過計算袋外數據的誤差(OOB error)可估算模型的誤差。研究表明，OOB error是RF模型泛化準確率的無偏估計[19]，在某些情況下比交叉驗證的效果更好[20]。

1.4 磷通量估算

根據望虞河的逐月出入湖水量(太湖流域管理局)和入湖斷面WY9(無錫312國道橋自動監測站)的TP以及反演的DTP，通過濃度與水量相乘即可計算望虞河出入湖的TP通量和DTP通量，PP通量為兩者之差。

2 模型構建

2.1 數據預處理

10個水質參數和4個氣象參數的描述性統計特征如表1所示，17個點位的DTP濃度為0.007 ～ 0.166 mg/L，平均值為0.042 mg/L，占TP的比例在8% ～ 98%之間，平均值為51%，可見DTP占TP的比例極差較大，與TP不存在簡單的線性關系。變量間的Pearson相關性特征如圖2所示，總磷、濁度和風向與DTP呈顯著相關，總磷還與溶解氧、總氮和濁度呈顯著相關，其他變量間存在不同程度的相關性。

表1 水質和氣象參數描述性統計Tab.1 Descriptive statistics of water quality and meteorological parameters

圖2 參數間的相關性特征Fig.2 Correlation characteristics of water quality and meteorological parameters

本研究數據集較小(102組)，按8∶2的比例隨機劃分訓練集和測試集，劃分后訓練集共有數據81組，測試集21組。

2.2 模型訓練與驗證

基于拆分的訓練集，輸出變量設置為DTP，其余13個參數作為輸入變量。隨機森林中有兩個顯著影響模型性能和運行效率的超參數，需要在開始學習過程之前人為設置：①ntree：指定隨機森林所包含的決策樹數量；②mtry：指定決策樹節點隨機選取的屬性個數。

首先優化mtry，將1到自變量(屬性)個數之間的所有整數賦值給mtry，觀察每個取值下模型的OOB error，OOB error最小時的mtry值為最佳值。然后優化ntree，觀察模型內誤差隨決策樹數量變化的情況，根據奧卡姆剃刀(Occam’s razor)原理，選擇使模型誤差穩定時的最小樹數量。

2.3 變量篩選

采用后向變量終止法進行變量篩選，即先將全部變量選入模型，每次擬合刪除重要性最低的預測變量，觀察模型性能是否有提升，若有提升，剔除該變量，余下變量重新擬合模型，重復上述步驟，直到模型性能下降或無明顯提升為止。變量重要性通過殘差平方和來度量，節點純度(node purity)增加等同于殘差平方和的減少，IncNodePurity即increase in node purity，代表了每個變量對分類樹每個節點上觀測值的異質性的影響，該值越大表示該變量的重要性越大。

2.4 模型評估

選用相關系數(R2)、均方根誤差(Root Mean Squared Error，RMSE)和納什效率系數(Nash-Sutcliffe Efficiency，NSE)評估模型性能和泛化能力。R2和NSE越接近1，RMSE越小表示模型性能越好。

本研究所有的數據分析與加工和模型構建工作在基于R4.1.2的RStudio軟件中實現。

3 結果與討論

3.1 模型參數

超參數mtry遍歷1到13(輸入變量個數)的模型OOB error和模型內誤差隨決策樹數量變化的結果如圖3所示，可確定最佳的超參數組合為mtry=8和ntree=700。按最佳超參數組合代入訓練集中的81組數據訓練模型。

圖3 RF模型的超參數(mtry和ntree)尋優Fig.3 The determination of optimal hyperparameters (mtry and ntree)for RF model

3.2 變量重要性分析

不同水質參數對預測DTP濃度的重要性如圖4所示，總磷(TP)、風向(WD)、濁度(FTU)、水溫(WT)、pH和溶解氧(DO)是重要性排序靠前的6個變量。

總磷涵蓋了各種形態的磷，水體中總磷主要以顆粒態磷的形式存在[21]，水體中的磷主要以懸浮顆粒物為媒介進行輸送，較高的濁度意味著水體中有更多的懸浮物等不溶性顆粒，進而影響磷的賦存形態及含量。水溫對沉積物中營養物質的釋放速率和釋放量及生物對磷的吸收和利用效率等均有重要影響[22]，pH同樣影響沉積物中營養物質的釋放[22]，合適的pH還有利于藻類生長[23]，因而水溫和pH是水體中磷循環的重要影響因素。溶解氧則是影響底泥沉積物磷素釋放[24]和微生物生長的重要因素。風場通過水面風擾動產生的剪切力對水體濁度產生影響，是泥沙再懸浮和遷移的關鍵[25]，風速重要性較低可能是因為輸入變量中已有濁度，通過部分依賴圖分析發現當風向為東風(90°附近)時對模型輸出結果影響最大，這可能和變量本身分布特征及流域盛行風向有關。

圖4 輸入變量的重要性排序Fig.4 The importance ranking of input variables

3.3 變量篩選與模擬結果

上述訓練好的模型命名為RF#1，在RF#1中剔除重要性最低的輸入變量RH，生成的模型為RF#2，再剔除WS的模型為RF#3，3個模型的性能如表2所示。

表2 RF模型性能比較Tab.2 Comparison of RF model performance

剔除RH后，模型在訓練集和測試集上的表現均有一定程度的提升，在測試集上的性能提升更為明顯。進一步剔除WS后，盡管模型在訓練集的表現有輕微提升，但在測試集上的性能幾乎沒有變化，訓練集和測試集的性能差異增大，會增加模型過擬合的風險。因此綜合比選后選擇RF#2作為最終模型，其模擬結果如圖5所示。

圖5 最終模型的模擬結果Fig.5 Simulation results of the final RF model

結合表2及圖5可見，模型RF#2的擬合效果良好，但當DTP濃度較高時，模擬值與實測值偏差較大。產生高值低估有如下原因：一是目前樣本量有限，DTP濃度高值(大于0.1 mg/L)較少，導致最終訓練出的模型傾向于低值端。二是隨機森林不能夠做出超越訓練集數據范圍的預測，這可能導致在某些噪音比較大的樣本集上隨機森林容易陷入過擬合，也是模型在測試集上的表現不如訓練集的原因。

3.4 磷通量計算結果

通過隨機森林模型反演的入湖磷通量計算結果如圖6所示。2010～2021年，望虞河入湖的DTP通量為8.44～73.14 t，占TP通量的比例為39.5%～57.6%，12年累計407.11 t，占TP通量的44.1%；PP通量為8.17～111.11 t，12年累計通量515.28 t，占TP通量的55.9%。望虞河入湖的磷以顆粒態磷為主，但近年來入湖DTP通量占TP通量的比例有升高的趨勢。

圖6 望虞河溶解態磷和顆粒態磷入湖通量Fig.6 DTP fluxes and PP fluxes of Wangyu River into Taihu Lake

根據水利部太湖流域管理局發布的《太湖健康狀況報告》(2010 ～ 2018年)，2010 ～ 2018年太湖全湖累計入湖磷通量為19260 t，同期望虞河入湖TP通量占比不足5%。結合貢湖年均蓄水量[7]和巡測的TP數據計算貢湖的蓄磷量，2010 ～ 2021年，貢湖的蓄磷量為13.86 ～ 23.17 t，平均值為20.76 t。望虞河年均入湖的TP通量是貢湖蓄磷量的3.9倍，最高可達9.1倍，其中易被生物利用的DTP年均入湖通量是貢湖蓄磷量的1.6倍，最高可達3.6倍，這將極大地促進灣內藍藻生長。無法被吸收的DTP則會和PP一起進入沉積物中，最高可達貢湖灣水體蓄磷量5.5倍，平均2.1倍的PP通量進入貢湖灣后，將大量沉積在灣內，成為灣內磷的重要內源和長期來源，在受到風浪等外力因素的擾動及藻類生物泵吸作用的影響時，向上覆水釋放蓄積的磷，有較高的釋放風險。

由此可見引江濟太帶來的磷通量對太湖全湖影響較小，但對貢湖灣形成復合的高負荷磷素沖擊。近年來貢湖灣有向藻型湖區轉化的趨勢，灣內有錫東、沙渚和金墅灣3個水源地，調水引致的磷輸入對太湖磷循環及水環境的影響需要引起足夠的重視并進一步采取措施保障水源地供水安全。

4 結 論

(1)本研究以望虞河-太湖為研究對象，基于隨機森林的機器學習模型，利用現有監測項目的歷史數據建立了反演DTP濃度的方法，最終模型的R2、RMSE、NSE分別為0.690、0.0110和0.651，具有良好的擬合和泛化性能，為歷史DTP濃度數據的反演提供了新思路。

(2)本研究10個水質參數和4個氣象參數間存在不同程度的相關性，其中總磷、風向、濁度、水溫、pH和溶解氧是本研究隨機森林模型中的重要預測因子，這些參數直接或間接地影響水體中磷的賦存形態、含量及循環過程。

(3)2010 ～ 2021年，望虞河入湖的DTP通量為8.44 ～ 73.14 t，PP通量為8.17 ～ 111.11 t，累計通量占TP通量的比例分別為44.1%和55.9%。望虞河年均入湖的TP通量是貢湖蓄磷量的3.9倍，DTP通量是1.6倍，對直接受水區形成復合的高負荷磷素沖擊，引江濟太調水引致的磷輸入對貢湖磷循環及水環境的影響需要引起足夠的重視。

(4)本研究通過隨機森林估算出不同形態的入湖磷通量，可為研究引江濟太工程對太湖磷循環及水質的長期影響提供關鍵數據，并為優化引江濟太引水方案提供技術支持。