基于随机森林的引江济太入湖磷通量计算

陈年浩，陆 昊，饶文昕，刘 彤，钱 新

(南京大学环境学院污染控制与资源化研究国家重点实验室，南京 210023)

前 言

跨流域调水作为控制湖泊富营养化问题的有力措施，在水质改善和富营养化控制过程中发挥积极作用[1-2]的同时，会挟带引水区及流经通道的营养盐和污染物等直接进入受水湖泊[3]，进而对流域水生态环境产生复杂影响。望虞河引江济太工程是太湖流域综合治理骨干工程之一，通过引长江水入太湖优化流域水资源配置。当前关于引江济太工程对受水区贡湖及太湖全湖的水质影响有较多研究[4～6]，调水是否造成太湖总磷浓度反弹尚有争议[7]，但调水过程中引入了一定量的磷进入太湖已被证实。磷是导致水体富营养化的重要元素之一[8]，相比于颗粒态磷，溶解态磷更易被生物利用[9]，耿雪等[10]指出应以溶解态总磷和溶解态无机磷来考量太湖表层水体磷的生物有效性。现有研究多关注总磷，鲜有关注不同形态的磷的迁移，当前太湖蓝藻暴发的风险依然较高，在未来规划引水水量增加的情况下，摸清引江济太不同形态的入湖磷通量具有重要现实意义。

一直以来，自动站驻测和人工巡测均只监测总磷而不区分磷形态，导致不同形态的磷通量无从计算。如何利用现有监测项目的历史数据反演不同形态的磷成为关键，传统的统计模型难以充分解释磷形态与总磷及其他环境因子的复杂非线性关系，随着大数据科学和人工智能的迅猛发展，机器学习的理论和方法在环境科学领域已得到广泛应用，其中随机森林(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)本研究通过随机森林估算出不同形态的入湖磷通量，可为研究引江济太工程对太湖磷循环及水质的长期影响提供关键数据，并为优化引江济太引水方案提供技术支持。