湖泊富营养化导致的藻类暴发一直是我国最为突出的水环境问题之一. 藻类过度生长是多种因素共同作用的结果，既包括水温、 光照、 风速等自然因素，也包括氮(N)、 磷(P)、 铁(Fe)等营养物过量排放的人类活动因素. 在诸多因素中，全球气候变暖背景下的水温变化与高强度人类活动所引起的N、 P排放增加被认为是导致湖泊富营养化最关键因素，因此，同时考虑水温、 N、 P因子变化的湖泊富营养化相关研究在逐渐增多，但温度与营养物对湖泊藻类生物量的交互作用等还需要深入研究[7]，比如水温、 N、 P促进藻类生物量增长的相对重要度的长期变化规律和季节性特征.

　　富营养化湖泊的藻类生长是自然界中一个非期望或非平均的现象，藻类生物量数据异质性很强，水华期间的藻类数据会呈“高峰厚尾”的分布，或存在显著的异方差等情况. 近年来在环境科学和生态学领域受到重视的分位数回归(quantile regression)方法特别适合处理这种波动性大、 异质性很强的环境数据. 该方法可针对回归变量任何一个分位点进行回归分析，且在存在极端值或重尾情况时仍能保持较好的稳健性，适宜处理应变量对自变量的极端响应，而不只是平均水平的响应，因此能更加全面地反映藻类生物量对水温、 N、 P等环境指标的响应特征. 本研究基于云南洱海长时间尺度(1990-2013年)的水质观测数据，运用分位数回归方法，按不同年份区间和不同季节分别分析洱海藻类生物量[以叶绿素a(Chl-a)表征]对N、 P、 水温的定量响应关系，探讨营养物因子与水温因子相对重要性的长时间尺度演变规律和季节性变化规律，对制定洱海富营养化控制策略提供科学依据.

　　1 材料与方法

　　1.1 研究区域

　　洱海是云南省第二大高原淡水湖泊，为滇西最大的断陷湖，跨洱源、 大理两县市，处于东经100°06′-100°17′，北纬25°36′-25°55′之间. 水面面积249.80 km2，汇水面积2 565.0 km2，最大水深21.0 m，平均水深10.5 m，库容28.8亿m3(图 1). 洱海是沿湖人民生活、 灌溉、 工业用水的主要水源地，是整个流域社会经济可持续发展的基础[14]. 洱海湖面多年平均海拔1 965.8 m，光照充足，辐射强，气温温和，为浮游藻类的大量繁殖提供了有利条件. 区年均气温15℃左右，年均降雨量1 055 mm，年均蒸发量1 970 mm. 流域水系发达，入湖河流大小共 117条.

图 1 洱海流域及水质监测点分布示意

　　1.2 数据来源

　　从云南省环境监测中心站、 中国大理洱海湖泊研究中心等环境监测和研究部门收集了洱海1990-2013年的水环境常规监测数据，全湖水质监测点为13个(图 1)，经纬度依次为： 25°51′36″N 100°10′12″E,25°51′18″N 100°11′24″ E,25°51′00″N 100°12′36″E,25°41′42″N 100°12′54″E,25°42′00″N 100°13′48″E,25°42′18″N 100°15′00″E,25°37′48″N 100°13′48″E,25°36′0″N 100°15′0″E,25°37′60″N 100°15′50″E,25°55′30″N 100° 6′54″E,25°54′54″N 100°8′42″E,25°54′36″N 100°10′48″E,25°47′29″N 100°11′43″E. 水质指标主要包括总氮(TN)、 总磷(TP)、 水温、 藻类叶绿素a(Chl-a)，采样频次为一年6次或12次. 各水质指标的主要测试方法如下：TN为过硫酸钾氧化紫外分光光度法，TP为钼锑抗分光光度法，水温为温度计现场实时测定，Chl-a为四波段分光光度法[15]. 数据主要统计特征见表 1，TP的变异系数最大(239.7%)，其次是Chl-a和TN，这可能与各指标的年份间差异和季节性差异有关. 各指标数据的Kolmogorov-Smirnov正态分布检验，除了水温数据呈正态分布以外，所有指标的数据都呈非正态分布.

　　表 1 洱海水温、 总氮、 总磷及藻类叶绿素a的统计值(n=1 419)

　　1.3 分位数回归方法

　　分位数回归模型(quantile regression)是依据因变量的条件分位数对自变量进行回归，得到所有分位水平的回归模型.

　　假设随机变量的分布函数为F(y)=Prob(Y≤y)，Y 的τ 分位数为满足F(y) ≥τ的最小y 值:

　　F(y)的τ 分位点Q(τ)由最小化关于ξ 的目标函数得到:

　　其中ρτ(u)称为检验函数:

　　假设因变量Y和自变量X 在τ 分位的线性函数关系为Y=X′β+ε，给定X=x时，Y的条件分布函数为FY (y | x)，则τ分位数为：

　　线性条件分位数通常表示为:

　　分位数回归能在不同的分位数τ得到不同的分位函数，残差计算方法不同于最小二乘法，具体如下：

　　1.4 数据处理与分析

　　考虑长时间尺度变化，将1990-2013年的监测数据按照相同年限(6 a)分成4组：1990-1995年、 1996-2001年、 2002-2007年、 2008-2013年，分别对这4组进行三元分位数回归分析. 为分析季节变化，将所有监测数据按春、 夏、 秋、 冬四季分成4组，按不同季节分别进行三元分位数回归分析，比较各因子的斜率值变化，并计算95%置信区间(CI)，CI覆盖0(包含0)，表示相关性“不显著”，而CI未覆盖0(不包含0)，说明相关性“显著”. 为便于比较3个因子的斜率值变化，在进行多元分位数回归之前，所有原始数据取对数(lg) 处理.

　　使用统计软件STATA V.12.0进行分位数回归分析和参数估计.

　　2 结果与讨论

　　2.1 水质指标的年际变化与季节变化

　　2.1.1 年际变化

　　从TN变化趋势图看，2001年之前，TN浓度呈缓慢上升趋势，但总体处于较好的Ⅱ类水质. 在2001-2003年期间有一次跳跃式增长，从Ⅱ类区跨入Ⅲ类区，近10年(2003-2013年)TN浓度维持在Ⅲ类水平，达到2001年之前浓度水平的近2倍[图 2(a)]. TP浓度年均值的变化趋势表现为有升有降，基本在Ⅱ类水平线上下浮动，浓度高峰出现在2002-2004年区间和2013年. 近10余年(2002-2013年)的TP浓度水平比前10年(1990-2001年)的浓度水平总体高出15%左右[图 2(b)]. 洱海的年平均水温在16-19℃之间周期性波动，未见明显的上升或下降趋势[图 2(c)].

图 2 洱海TN、 TP、 水温和藻类Chl-a的逐年变化(1990-2013年)

　　从藻类Chl-a变化趋势看，2000年之前洱海全湖Chl-a维持在低浓度水平，变化范围为1.0-3.5 mg·m-3. 2000-2002的3年期间，Chl-a浓度呈直线上升，藻类生物量从低水平跨越至较高水平. 2003年之后的Chl-a平均浓度达到2000年之前的10余倍[图 2(d)]. Pearson相关性分析显示藻类Chl-a与TN以及水温均呈显著正相关性(P<0.01)，而与TP的相关性不显著(P>0.05).

　　2.1.2 季节变化

　　4项水质指标具有明显的季节性特征，均表现为夏、 秋两季较高，且两季数据接近，而春、 冬季数值较低. 8-10月TN浓度最高，5月是一年中TN上升的起点，11月是TN开始下降的拐点[图 3(a)]. 一年中TP浓度在4-7月是洱海上升速率最快的时段[图 3(b)]. 水温在夏季和秋季最高，分别为23.3℃和20.2℃，很适宜藻类生长[图 3(c)]. 藻类Chl-a浓度在6-7月出现跳跃式上升，平均值由10 mg·m-3跃升至20 mg·m-3，在7-11月维持在最高水平[图 3(d)].

图 3 洱海TN、 TP、 水温和藻类Chl-a的季节变化

　　2.2 藻类对氮、 磷及水温的长时间尺度响应特征

　　通过三元分位数回归方法分析洱海N、 P和水温对藻类Chl-a的耦合效应随时间变化特征. 4个年份区间的lg(TN)、 lg(TP)和lg(水温)的斜率值随分位点的变化谱图见图 4，斜率值及其95%置信区间的统计特征见表 2.

　　表 2 不同年份区间的多元分位数回归的斜率参数统计

蒸发器，污水处理设备：www.gdgrsw.com