基于非参数回归改进的太湖总磷压力响应模型

华祖林，韩爱秋

(1．河海大学浅水湖泊综合治理与资源开发教育部重点实验室，江苏南京 210098;2．河海大学水资源高效利用与工程安全国家工程研究中心，江苏南京 210098;3．河海大学环境学院，江苏南京 210098)

湖泊参照状态指湖泊受到人类活动影响最小或环境系统可达到的最佳状态[1]，它是制定水质标准的科学依据，也是评估、预防、控制和管理湖泊富营养化的基础，更是水污染治理和生态修复期望达到的“终极”状态。压力-响应模型是最适用于确定浅水湖泊参照状态的模型之一。Lamon等[2]探索性地以贝叶斯多层次线性回归模型确立了美国不同类型湖泊压力变量和响应变量的响应关系;US EPA[3]介绍了利用线性回归、多元线性回归和非参数拐点分析等建立压力-响应模型的方法;Haggard等[4]利用线性回归和分类回归树相结合的方法确定了美国红河流域的营养物基准值;Qian等[5]构建了连续变量贝叶斯网络建模框架，该框架结合了贝叶斯网络模型和传统经验统计模型，并确定了美国俄亥俄州溪流的营养物基准值;在中东部平原湖泊生态区、云贵高原湖泊生态区等8个湖泊营养物一级生态区的基础上[6-7]，Huo等[8-11]引入压力-响应模型来制订中国受人类活动影响较严重湖泊的参照状态，利用线性回归压力-响应模型制订了东部平原湖泊生态区和云贵湖泊生态区的营养物基准，并在考虑季节影响水质的情况下采用分类回归树与拐点分析相结合的方法建立压力-响应模型深入研究湖泊生态分区参照状态的确定;吴超等[12]采用分类回归树和线性回归分别建立了压力-响应模型，并通过相互验证确定了太湖TN、TP等的参照状态;张亚丽等[13]利用贝叶斯拐点分析等3种非线性方法确定了中国9个典型湖泊的参照状态;霍守亮等[14-15]利用分类回归树、拐点分析法、贝叶斯层次线性回归建立了压力-响应模型，针对不同湖泊生态分区进行参照状态的确定和案例研究，并进行了归纳总结。

以上研究成果有值得商榷和改进之处，如线性回归模型要求响应关系是线性的、因变量的误差呈正态分布、采集的数据样本独立等，而实际上，光照、水深等环境因素会影响变量间的响应关系，导致其难以呈现出线性;此外，湖泊实际监测数据来源复杂，且为时间序列，很难满足线性回归分析设定的正态性和独立性等假设理论。贝叶斯层次回归模型虽然能有效缓解数据测量误差的问题，但需要在调查时以明显分层的方式采集原始数据，以便分层建立模型。分类回归树压力-响应模型适用于识别对模型结果变化有显著贡献的压力变量，但在数据量较少时缺乏稳健性，压力变量很小的变化可能引起模型结果较大的变化，使结果的准确性得不到保障;拐点分析法需要进一步研究低于营养物拐点对应的响应变量数值是否能够支持水体的指定用途[14-15]。这些限制性条件会导致推断的参照状态浓度过高或过低。压力-响应模型方法的改进一直是学者们孜孜不倦探索的课题。

本文针对原有压力-响应模型需要满足设定条件等特点，考虑到响应关系应是非线性的，以及异常点会降低模型稳健性等情况，尝试性地采用非参数LOWESS稳健回归来改进浅水湖泊的压力-响应模型。改进模型能消除异常点对估计结果的影响，对于非线性、非齐次问题有较好的效果。太湖作为东部平原湖泊生态区浅水湖泊的典型代表，参照状态的确定尤为重要[16-17]，以太湖为研究对象进行湖泊参照状态的确定计算。

1 非参数LOWESS稳健回归模型的构建

非参数LOWESS稳健回归是拟合难以用具体函数描述响应关系数据的有效方法，其基本思想是先用局部多项式进行拟合，然后定义稳健的权数并进行光滑处理，重复运算数次后就可以得到自变量与因变量稳健的拟合曲线[18]。具体步骤如下：

步骤1：将压力变量TP和响应变量Chl-a关系曲线表示为m(x)，利用局部多项式进行拟合，得到拟合值 ^m(xi)：

式中：xi为TP的观测值，mg/L;β^(xi)为对局部参数β(xi)的最小二乘加权回归。

步骤2：定义并计算稳健权数δk：

其中 ek=m(xk)-m^(xk)

式中：B为二次核函数，定义为B(x)=(1-x2)2I(I为符号函数，当|x|＜1时I=1;当|x|≥1时I=0);ek为残差;xk为加权的xi;s为残差绝对值的中位数，s=median(e1， e2，…， en)。

步骤3：利用δkwk(xi)对m(x)进行局部加权最小二乘估计，得到新的模型残差ek。其中

式中：W为距离权函数;hi为光滑参数，hi= xi-xj，xj(j=1，2，…，n)为 xi相邻的 TP 观测值，hi为计算得出的r个数中最小的一个。

步骤4：重复第2步和第3步，最后第N次迭代得到的拟合值产生TP和Chl-a的响应关系曲线。

稳健权数可将异常值排除在外，并且初始残差大(小)的观测值在下一次加权最小二乘法中的权数就小(大)，重复几次后就可把异常值不断地排除在外，最终得到稳健的响应关系曲线。Cleveland[18]推荐迭代次数N=3。

非参数LOWESS稳健回归压力-响应模型构建流程见图1。

图1 压力-响应模型构建流程

2 模型数据与压力变量的确定

2.1 模型数据

利用LOWESS稳健回归方法建立压力-响应模型推断太湖参照状态，数据来源于文献[19]。文献[19]记录了图2太湖湖泊生态系统国家野外科学观测研究站(简称太湖站)全太湖32个野外站点中8个站点自1991—2006年的逐月观测数据。考虑到参照状态是受人类活动影响最小或环境系统可达到的最佳状态，所以本文采用位于湖心区域受人类活动影响相对较小的7号和8号站点的数据。华祖林等[20]采用季节分解的非参数法发现1999—2001年太湖湖心区域的TP与Chl-a的观测值比1995—2006年的观测值更适合用于推断太湖的参照状态。本文使用7、8号站点1999—2001年共3年108个观测值进行分析。

图2 太湖站监测点图

2.2 压力变量确定和数据处理

首先对模型数据进行分析。1999—2001年，太湖湖心TN、TP以及Chl-a质量浓度平均值、中位数、均方差等见表1。1999—2001年太湖湖心观测值稳定性按从大到小顺序为TP、TN、Chl-a。在峰度和偏度检验中，TN、TP和Chl-a均存在正偏态现象，偏重程度按从大到小顺序为Chl-a、TN、TP，且三者均偏离正态分布。综合考虑数据稳定性及峰度和偏度检验结果，在选择压力变量时应优先考虑TP。

在研究湖泊氮磷与Chl-a关系时，TN与TP质量浓度比值(ρ(TN)/ρ(TP))是一个关键参数，常根据ρ(TN)/ρ(TP)判别藻类受营养盐限制的类型。ρ(TN)/ρ(TP)比较大(＞17)时藻类受磷限制，ρ(TN)/ρ(TP)比较小( ＜10)时藻类受氮限制，ρ(TN)/ρ(TP)中等(10～17)时藻类受二者共同限制。表1显示太湖1999—2001年TN质量浓度平均值为1．69 mg/L，TP 质量浓度平均值为 0．08 mg/L，ρ(TN)/ρ(TP)=21．12，大于 17。 可见，1999—2001年太湖为磷限制型湖泊。张波等[21]研究表明，太湖水体存在一定的固氮作用，固氮速率为每小时1．53 ng/L，导致水体中TN与Chl-a响应关系不明显。因此分析太湖TP和Chl-a的响应关系更符合压力-响应模型的释义。

图3 TP与Chl-a幂变换数据Q-Q图

表1 1999—2001年TN、TP与Chl-a统计特征比较

由于原始逐月观测数据TP与Chl-a存在数量级差别，因此需要对原始数据进行处理。本文采用Qian[22]提出的幂变换对TP和Chl-a数据进行处理。图3是TP数据经过对数变换后和Chl-a 0．1次方变化后的Q Q图。从图3可以明显地看出，TP数据点与标准正态分布直线拟合效果较差，存在严重的偏离。结合表1峰度和偏度检验结果，说明1999—2001年太湖湖心TP观测数据和幂变换后均不满足正态性假设，所以基于上述假设建立的压力-响应模型具有不合理性，进一步说明建立非参数方法压力-响应模型的必要性。

在采用压力-响应模型推断参照状态时，需要给定一个响应变量基准值。在保证水体使用功能的前提下推断营养物基准时，Chl-a是联系营养物浓度的重要响应变量[8]。参考 USEPA[23]和霍守亮等[14，24-25]对Chl-a基准值的研究，本文在建立压力-响应模型时设定4．73 μg/L Chl-a为太湖响应变量基准值。

3 模型验证

采用LOWESS稳健回归建立1999—2001年太湖湖心区域TP与Chl-a的压力-响应模型，平滑参数f=0．2，迭代次数 N=3。图4为太湖湖心区域1999—2001年Chl-a与TP压力-响应模型及模型残差。从图4的散点可以看出，Chl-a与TP并不呈现出简单的线性响应关系，而是复杂的非线性响应关系。根据确定的响应变量Chl-a基准值4．73 μg/L，可以得出压力变量TP的参照状态质量浓度为0．018 mg/L。此外，结合置信区间和模型残差评估模型结果的准确性，利用自助法求得1999—2001年太湖TP参照状态质量浓度的95%置信区间为0．013～0．030 mg/L。 从图4(b)可以看出，模型残差均值为0．0mg/L左右，下25%分位点为-0．06 mg/L，上25%分位点为0．04mg/L，说明LOWESS稳健回归方法可用于建立压力-响应模型进行参照状态的推断。

图4 TP与Chl-a压力-响应模型及模型残差

表2显示的是不同地区TP参照状态质量浓度。郑丙辉等[26]应用频率分析法建立了太湖TP参照状态，并结合时间参照状态法和沉积物反演法进行验证，综合分析得出TP参照状态为0．03 mg/L;华祖林等[20]利用1999—2001年TP数据5%分位点的值作为太湖参照状态质量浓度，其结果为0．03 mg/L，并利用几何分布分块自助法确定其95%置信区间为0．025～0．046 mg/L。采用频率分析法在确立湖泊参照状态时需要确定频率分布分位点，通常都是选择上(下)25%或5%，由于25%或5%是人为确定的，因而结果受到人为影响大。吴超等[12]利用分类回归树和线性回归法建立的压力-响应模型相互验证，确定太湖流域TP质量浓度基准值为0．67mg/L。

表2 不同地区TP参照状态浓度

霍守亮等[14]利用非参数拐点分析和贝叶斯拐点分析，得到不同湖泊生态区压力变量TP突变点的质量浓度范围为0．015～0．222mg/L，利用简单线性回归模型得到中东部湖泊生态分区TP参照状态质量浓度范围为(0．022±0．007)mg/L;顾莉等[27]将1960年总碱度实测值代入其建立的改进MEI模型得到太湖TP质量浓度为0．018mg/L。在20世纪60年代对太湖的调查研究表明，太湖TP质量浓度为0．01～0．05 mg/L[28]。 太湖蓝藻暴发最近30年才出现，60年代太湖的状态应该比较接近其历史上受人类活动影响较小时的情况。本文采用季节分解的非参数局部线性回归模型推断得出更适合建立太湖参照状态的数据，而且非参数LOWESS稳健回归不需要设定分位点，计算过程受人为影响小，计算得到的太湖TP参照状态质量浓度为0．018 mg/L，在线性回归和拐点方法得出的参照状态质量浓度范围内。因此，本文得到的结果可以作为太湖参照状态质量浓度，并且是可信的。

4 结 语

以季节分解的非参数方法验证的适合推断太湖参照状态的数据为基础，结合峰度和偏度检验，并考虑湖泊营养盐的限制问题，以TP作为压力变量，构建了稳健性高、对非线性问题适应性强的LOWESS稳健回归的湖泊压力-响应模型。根据设定的Chl-a推断TP的参照状态质量浓度为0．018 mg/L;通过自助法得出其95%置信区间为0．013～0．030 mg/L，同时从多角度对所建模型进行了验证，表明结果较为合理可信。