滇池经过20多年的治理，已经遏制住了水质恶化的趋势，目前外海水质在V类~劣V类之间波动(高伟等，2013；Wang et al., 2012；王圣瑞等，2014).但从单个水质指标来看，BOD₅和COD基本维持在III类以下，主要超标因子为TN和TP(刘永等，2012；高伟等，2013)，而降低营养盐浓度的关键是改善入湖河道水质(邹锐等，2011；Baresel et al., 2007).就当前滇池流域管理来看，单纯的浓度监测已不能满足流域综合管理的要求，入湖污染负荷总量控制及其预报预警尚未引起重视.因此，有必要构建一套滇池流域TN、TP入湖负荷预报预警系统，建立河道TN、TP入湖负荷与湖体水质响应关系，形成流域-湖泊一体预报预警系统，从而有效监管入湖TN和TP污染排放情况.

目前，流域和河流的预报预警系统多集中在水质在线监测、水质超标和突发事故预警等方面，早期的如1997年德国与其他欧洲国家投入使用的多瑙河水污染事故应急预警系统(Pintér，1999；Jansky et al., 2004)，以及大庆市2004年建立的黎明河在线水质监测和管理系统(Yang et al., 2008).国内还基于GIS平台开发了一些水质预警系统，如汉江水质预警系统(窦明等，2002)、辽河流域水安全预警系统(王耕等，2007)、钱塘江水质预警预报系统(李佳等，2008)，这些预警系统侧重于整体功能的设计和信息发布，对于预警规则和体系的阐述较少.迟晓德等(2012)基于污染物总量控制思想构建了流域预警体系，特点在于利用层次分析法确定了流域预警指标，提出综合预警系统，并将预警等级划分为4个级别，但预警方式较为单一.且这些预警系统多以行政区划作为预警单元，对于流域(污染负荷产生单元)和纳污水体(容纳污染负荷的响应单元)整体预警较少.

因此，本文以子流域为基本单元建立滇池流域TN、TP入湖负荷预报预警系统，采用流域模型HSPF获取60个子流域入湖流量，利用贝叶斯递归回归树算法(BRRT)拟合各子流域入湖负荷计算方程，获取TN、TP负荷，并以此为计算核心提出预警体系.基于1999—2010年TN、TP负荷确定预警规则，考虑流域内降水和面源负荷的时间分配特征，建立两套预警体系，选取流域内以农业面源为主的柴河子流域以并“十二五”控制目标作为排放限值应用此系统对其2011年不同月份预警.

滇池流域(如图 1所示)位于昆明市境内，地处长江、红河、珠江三大水系分水岭地带(24°29′~25°28′N，102°29′~103°01′ E)，涉及五华、盘龙、官渡、西山、呈贡、晋宁、嵩明7个县(区)2920 km²的区域，滇池面积约309.5 km²(高伟等，2013).滇池是中国第六大淡水湖泊，流域面积不到云南省面积的1%，却集中了云南全省经济总量的1/3，由此可见滇池的重要性.而滇池目前污染严重，从“九五”开始，连续4个5年规划都将滇池治理纳入国家“三河三湖”治理规划的重点，为了加强滇池的保护和管理，防治水污染，改善流域生态环境，云南省自2013年1月1日起实施了《云南省滇池保护条例》.

图 1(Fig. 1)

图 1 滇池流域分布图 Fig. 1 Distribution of Lake Dianchi Watershed

滇池流域有7种主要污染源，高伟等(2013)就此建立了2008年的TN、TP排放清单(表 1).TN和TP排放量在空间分布上差异显著，在占滇池流域总面积16.2%的10个子流域(多分布于流域北部城区和东北部种植区)，却贡献了70%以上的TN排放量.TP排放量较为分散，累积贡献70%的子流域为21个，面积占37.9%.TN和TP排放量最高的子流域位置十分类似，排放量最高的前10名中有7个为相同的子流域.本文在此基础上建立滇池流域TN、TP预报预警系统，对入湖污染排放监控与管理给出一定科学依据.

表 1(Table 1)

表 1 2008年滇池流域主要污染源TN和TP贡献率 Table 1 Emissions of TN and TP from different sources in Lake Dianchi Watershed in 2008

表 1 2008年滇池流域主要污染源TN和TP贡献率 Table 1 Emissions of TN and TP from different sources in Lake Dianchi Watershed in 2008 污染物类型 TN贡献率 TP贡献率 点源 城镇生活 72.70% 42.80% 畜禽养殖 3.60% 13.70% 第三产业 3.70% 9.80% 工业 1.90% 4.30% 小计 81.80% 70.60% 非点源 种植业 11.70% 11.90% 城市径流 5.80% 9.60% 农村生活 0.60% 7.90% 小计 18.20% 29.40% 合计 100% 100%

采用HSPF(Hydrological Simulation Program-Fortran，Bicknell et al., 2001; 董延军等，2009)作为水文模拟工具获取各入湖河流的流量.HSPF是一个半分布式流域模型，能模拟流域内长时间连续的水文和水力过程，水文模块是在SWM(Stanford Watershed Model)斯坦福模型的基础上发展而来的，应用广泛，研究比较成熟(董延军等，2009；李兆富等，2012).HSPF模型在流域水文模拟中表现出了较高的精度和可靠性，在长时间水文过程模拟中误差较小(Albek et al., 2004)，与SWM模拟精度相似(Johnson et al., 2003)，且比SWAT(Soil and Water Assessment Tool)精度高(Xie et al., 2013).

滇池流域水文模拟的HSPF模型已由北京大学滇池水专项课题组完成(郭怀成等，2013；高伟等，2014).需要的数据包括 DEM(1∶50000)、流域边界及河道尺寸、土地利用(30 m×30 m)、土壤类型、18个站点的1999—2010年气象和水文数据(日最高/最低气温(℃)、日照时间(h · d^-1)、小时蒸发量(mm · h^-1)和日降雨量(mm · d^-1)).基于建好的流域HSPF水文模型，输入气象和水文数据获取各子流域入湖流量.

采用BRRT(Bayesian Recursive Regression Tree，周丰等，2010；Zhou et al., 2014)拟合各子流域污染负荷排放的分段函数以获取TN、TP入湖负荷.HSPF模型的水质模块对输入数据要求较高，而且输入文件需要第三方软件操作，对于预报预警系统的稳定运行是一个较大限制.因此，选择BRRT替代模型在保证预测准确性的情况下简化计算过程.BRRT模型在水库水质模拟(周丰等，2010)、大尺度水文模拟(周丰，2014)、农田氮淋溶(徐鹏等，2013)和农田N₂O排放模拟(Zhou et al., 2014)中得到应用.类似利用高级统计模型替代非线性响应机理模型实现高效计算已经在流域系统优化中有很多研究(Qin et al., 2007；He et al., 2008；郭怀成等，2012).BRRT原理是通过随机搜索和局部贪婪搜索，以似然函数值最大为目标，建立因变量和自变量之间的时间线性分段函数.

首先，针对TN、TP负荷分别建立流域入湖负荷计算方程：

式中，0.011574为单位转换系数；i为入湖流域编号；t为某天；C_(i，t)为流域入湖口日平均浓度(mg · L^-1)；Q_(i，t)为流域入湖口日平均流量(m³ · s^-1)；Load_PS_(i，t)为点源排放负荷(kg · d^-1)；Load_NPS_(i，t)为面源排放负荷(kg · d^-1)；k为点源折减系数，由于点源排放相对规律，因此，以只考虑点源排放情况的HSPF模型模拟的入湖口负荷与点源排放负荷的比值表示.

式(1)中，k×Load_PS_(i，t)点源负荷比较容易获取，重点在于Load_NPS_(i，t)面源负荷的计算.采用BRRT基于1999—2010年滇池流域HSPF水质模型的数据样本，建立面源负荷回归的线性分段函数，因变量为面源污染负荷(Load_NPS_(i，t))，分类变量为化肥施用量(FERT_(i，t))和流域入湖流量(Q_(i，t))，回归方程中自变量则只包括Q_(i，t)，函数形式如下：

式中，a、b为回归系数；FERT_(i，t)为折纯后的化肥施用量(kg · d^-1).

从而建立了不同影响因素(FERT_(i，t)、Q_(i，t))水平下面源污染负荷Load_NPS_(i，t)的响应关系：

式中，p(T X，Y)为入湖面源负荷的响应决策树，X为影响因素，Y为面源负荷，Θ为模型参数，p(T)和p(Θ T)为决策树和参数的先验分布，p(Y X，Θ，T)为似然函数.

BRRT具体原理和应用请参考文献(周丰等，2010)，文中重点给出参数设置和计算过程.具体参数设置包括：10次重启(经验而定)，每次2万次迭代，决策树先验分布的参数α和β及参数c和λ分别为0.95、1、1和0.1173(为敏感性分析最佳组合).

最终，由式(1)和(2)可得，流域入湖污染负荷计算方程为：

HSPF模型的校准和验证见文献(郭怀成等，2013；高伟等，2014).面源负荷BRRT模型的校准和验证如下.模型中重点考虑农业面源负荷，选取流域内以农业面源负荷为主的柴河子流域说明校准和验证过程.样本来源为滇池流域HSPF模型1999—2010年逐日模拟结果，分别随机选取1100个有效数据进行校准与验证(1000个用来校准，100个用来验证)，TN、TP入湖负荷校准与验证结果分别如图 2、3所示.

图 2(Fig. 2)

图 2 TN入湖负荷BRRT模型校准与验证(a.校准期，b.验证期) Fig. 2 BRRT model calibration and validation of TN pollution loading(a.calibration; b.validation)

图 3(Fig. 3)

图 3 TP入湖负荷BRRT模型校准与验证(a.校准期，b.验证期) Fig. 3 BRRT model calibration and validation of TP pollution loading(a.calibration; b.validation)

由图 2和3可以看出，BRRT模拟TN入湖负荷校准和验证的可决系数R²分别为0.8383和0.8996，模拟TP入湖负荷校准和验证的可决系数R²分别为0.8207和0.8250，说明在柴河子流域BRRT拟合函数可以替代HSPF机理模型定量计算入湖污染负荷.当然对于全流域来说，当R²达到0.7以上就可以认为拟合良好.但也存在个别子流域R²小于0.5的情况，这些子流域集中在滇池北部和东部的城市区，可能是由于这几个子流域的城市面源贡献较大，而在BRRT模拟过程中无法量化城市面源负荷.从全流域整体模拟效果来看，BRRT拟合函数可以达到预报预警系统TN、TP入湖负荷的计算要求.

滇池流域污染负荷预警的目的是根据TN、TP的年控制目标，以历史污染负荷和模型预测为基础，对当年各月份实际负荷及可能的全年总负荷状况进行预报，评价不同的预警等级和行动方案.提出预警指数(EWI，Early Warning Index)的概念进行预警等级划分，预警指数EWI可以作为量化预警等级的指标，当EWI>100%时表示被评价数据大于标准值(负荷超过控制目标)，EWI的具体计算后续将给出.本文中预警等级共划分为5个水平(表 2).

表 2(Table 2)

表 2 预警等级 Table 2 Early warning grades

表 2 预警等级 Table 2 Early warning grades 预警水平 预警指数(EWI) 预警信号 行动措施 正常水平 EWI≤50% 绿色预警 正常监测与监管 启动水平 50%<EWI≤ 80% 黄色预警 持续监测与监管 警告水平 80%<EWI≤100% 橙色预警 增加流域监测频率，对重点污染源进行排查和不定期抽查 超标水平 100%<EWI≤150% 红色预警 严格实施监测，并实施重点污染源减排措施，控制入湖污染物负荷 严重超标水平 EWI>150% 黑色预警 严格控制河流沿线企业污染源的排放，对排入河流的污染物进行源追踪，对城市和生活污染源加大处理力度，并启动制定好的应急预备方案

以柴河子流域为例，制定预警规则从而计算预警指数，其它子流域制定方法一致.基于1999—2010年柴河子流域TN和TP入湖负荷计算每个月TN和TP的平均排放率(每个月排放的污染负荷占全年排放负荷的比)和月平均累积排放率(图 4和5).可以看出，6月底之前TN和TP累积排放率分别为34.0%和31.5%，6—11月入湖污染负荷增加，主要是由于降雨增加和农田施肥(滇池流域4—7、10—12月施肥)的原因.显然，对6月之前和之后采用同一套预警规则不合理，因此，制定两套预警体系(EWS，Early Warning System).

图 4(Fig. 4)

图 4 TN污染负荷月排放变化率与贡献率(a.TN累积排放率与月排放率变化图，b.TN月贡献率) Fig. 4 Monthly change rate and contribution rate of TN pollution loading(a.Cumulative emission rate and monthly emission rate of TN; b.Monthly contribution rate of TN)

图 5(Fig. 5)

图 5 TP污染负荷月排放变化率与贡献率(a. TP累积排放率与月排放率变化图，b. TP月贡献率) Fig. 5 Monthly change rate and contribution rate of TP pollution loading(a.Cumulative emission rate and monthly emission rate of TP; b.Monthly contribution rate of TP)

1)预警体系一(EWS-1)

针对6月以前预警：考虑到只用当前排放量与排放限值对比，则6月之前预警状态一直处于正常水平，不能达到实际预警效果.因此，采用累积负荷排放率与历史平均累积负荷排放率对比作为预警指数(EWI)，具体计算如下：

式中，i为1~6；R_i为i月累积负荷排放率；Load_i为i月累积负荷排放量(kg)；Load₀为流域负荷排放限值(kg).

式中，EWI为预警指数，R_0-i为i月历史平均累积负荷排放率.

2)预警体系二(EWS-2)

针对6月以后预警：①由于8月份流域TN和TP的平均累积排放率均超过了60%，知道现状负荷排放量占排放限值的比例很有意义，因此，采用现状负荷排放量占排放限值的比作为预警指数EWI₁；②采用累积负荷排放率与历史平均累积负荷排放率的比作为预警指数EWI₂；③现状排放量加预测负荷计算预警指数EWI₃，6月以后到12月的预测跨度相比6月以前要短，预测准确度将提高，叠加预测负荷的预警能更好的控制全年污染负荷排放.最终，给出3个预警指数EWI₁、EWI₂和EWI₃，并用最严格的作为综合预警指数EWI发布预警信号.具体计算如下：

对于现状排放量加预测负荷的预警考虑两种情况：①由于今年和去年流域治理力度相差不会很大，所以采用现状负荷排放量加去年该月到12月的负荷排放量作为综合负荷排放量Load_Z1；②利用皮尔逊III分布作为频率计算的线型，用韦布尔公式计算经验频率，通过适线法调整频率曲线的统计参数和设计值，对大观楼1951—2010年年降雨量进行频率分析，得出丰、平枯典型年降雨数据：特枯年，2009年；枯水年，2003年；平水年，2006年；丰水年，2001年；特丰年，1999年.将计算年份匹配到其中一种降雨典型年，利用该典型年的降雨数据并考虑计算年份的点源和化肥施用削减率基于HSPF-BRRT预测该月到12月的负荷排放量，加上截止该月的现状累积排放负荷作为综合负荷排放量Load_Z2.最终，取Load_Z1和Load_Z2中较大值与排放限值对比得出预警指数EWI₃.具体计算如下：

式中，j为i+1~12；Load_j为去年第j月负荷排放量(kg).

式中，k为i+1~12；Load_k为考虑点源和化肥施用削减率在匹配降雨典型年基于HSPF-BRRT预测的第k月负荷排放量(kg).

最终，以EWI评价预警等级发布预警信号，同时给出EWI₁、EWI₂和EWI₃作为决策参考.预警体系选择、预警指数计算和预警等级评价如图 6所示.

图 6(Fig. 6)

图 6 预警指数计算流程 Fig. 6 Calculation process of early warning index

以柴河子流域为例，应用上述预警体系对2011年2、4、9和11月分别预警.《滇池流域水污染防治“十二五”规划编制大纲》(中国环境科学研究院等，2010)中指出，在“十二五”末期外海南岸控制单元的入外海通量要在2008年污染负荷的基础上减少30%以上，以“十二五”规划目标作为柴河子流域的污染负荷排放限值，计算可得TN排放限值为130.2 t · a^-1，TP排放限值为6.8 t · a^-1.2011年匹配到降雨典型年为2009年，应用上述预警体系计算2011年柴河子流域污染负荷、EWI和预警等级，具体见表 3和4.

表 3(Table 3)

表 3 2011年柴河子流域TN预警 Table 3 TN early-warning of Chaihe catchment in 2011

表 3 2011年柴河子流域TN预警 Table 3 TN early-warning of Chaihe catchment in 2011 月份 Load0/t R0-i Ri EWS-1 EWS-2 预警等级 EWI Loadz1/t Loadz2/t EWI1 EWI2 EWI3 EWI 2 130.2 9.97% 18.16% 182.1% - - - - - - 黑色预警 4 19.12% 34.6%4 181.2% - - - - - - 黑色预警 9 70.91% 88.76% - 217.86 150.08 88.7% 125.2% 167.3% 167.3% 黑色预警 11 93.50% 103.7% - 159.71 148.48 104.4% 110.9% 122.7% 122.7% 红色预警

表 4(Table 4)

表 4 2011年柴河子流域TP预警 Table 4 TP early-warning of Chaihe catchment in 2011

表 4 2011年柴河子流域TP预警 Table 4 TP early-warning of Chaihe catchment in 2011 月份 Load0/t R0-i Ri EWS-1 EWS-2 预警等级 EWI Loadz1/t Loadz2/t EWI1 EWI2 EWI3 EWI 2 6.8 9.13% 16.17% 177.1% - - - - - - 黑色预警 4 17.97% 32.06% 1784% - - - - - - 黑色预警 9 77.19% 72.79% - 9.61 6.62 72.8% 94.3% 141.3% 141.3% 红色预警 11 94.42% 88.53% - 7.31 6.57 88.5% 93.8% 107.5% 107.5% 红色预警

由表 3和4可以看出，2011年柴河子流域TN和TP污染负荷不同月份预警基本处于红色和黑色预警，在2011年化肥施用水平和污染处理力度下，必然会超过“十二五”规划的目标，因此，流域要对农业面源施肥合理管理，并加强点源处理和实施流域截污，从而减少入湖污染负荷.

1)BRRT替代模型在流域内以农业面源为主的柴河子流域拟合效果良好，TN、TP入湖负荷校准和验证的可决系数均大于0.8，拟合函数可以替代机理模型定量计算入湖污染负荷.

2)提出了两套以预警指数(EWI)评价预警等级的滇池流域TN、TP入湖负荷预警体系(EWS).预警体系一(EWS-1)适用于6月份之前，以月累积负荷排放率和月历史平均累积负荷排放率的比值表征EWI.预警体系二(EWS-2)适用于6月份之后，考虑现状排放量、历史平均排放量、排放限值和预测排放量之间的关系计算3个预警指数，以其最严格者作为综合预警指数EWI评价预警水平发布预警信号.

3)根据滇池流域水污染防治“十二五”规划目标，确定柴河子流域2011年TN排放限值为130.2 t · a^-1，TP排放限值为6.8 t · a^-1，并对2011年2、4、9和11月分别进行模拟预警，结果显示，预警等级基本处于红色和黑色预警，需要加强农业面源施肥管理、点源处理和流域截污等，该系统能够为控制入湖污染负荷提供预警支持.