基于拐点分析的汉江水华暴发突变与归因研究

程兵芬 ，夏瑞，张远*，张楠 ，张新飞 ，刘成建,

1.北京师范大学水科学研究院，北京 100875；2.中国环境科学研究院，北京 100012；3.西北大学城市与环境学院/陕西省地表系统与环境承载力重点实验室，陕西 西安 710127

随着全球水体富营养化问题的加剧，国内外多个大型河流出现了严峻水华问题，严重影响了当地河流的水体结构与水生态环境，对人类饮水安全造成了极大的威胁（Jung et al.，2011；窦明等，2002；张远等，2017）。通常富营养化容易发生在水体流速相对静止的湖泊和水库，然而受多种因素的综合影响，国内外诸多大型河流频繁发生水华现象（Whitehead et al.，2015；梁开学等，2012；夏瑞等，2020）。营养盐含量、水文条件和气候条件的综合作用是水华发生的主要因素（Carey et al.，2012）；如营养盐总氮和总磷作为藻类生长的必要条件（吴兴华等，2017）；流速与流量通过改变水体滞留时间影响藻类的生长，水位可以通过影响水体分层现象影响藻类密度的变化（Ji et al.，2017）；降水的时空变化则会影响藻类的频率与强度（Phlips et al.，2007）；水温则会影响不同藻类物种的耐受性，不同的温度条件下，优势藻类结构有很大差异等（Jung et al.，2009；Kim et al.，2014）；如蓝藻适宜温度在30—35 ℃，绿藻适宜温度均在20—35 ℃，蓝绿藻适宜温度均在25—35 ℃，硅藻在0—15 ℃都能够的较好的生长，硅藻相比较蓝绿藻更适合低温条件（Chen et al.，2014）。综合分析营养物质、温度和水文条件的是探索河流藻华机制的必要条件（Kim et al.，2018）。在当前气候变化和人类活动共同作用背景下，河流水华未来的演变趋势及其主要影响因素也是当前亟需探索的科学问题。

汉江自 1992年第一次硅藻水华暴发以来，已前后暴发 10次不同程度的水华事件，严重时覆盖河段达到河口以上500多千米，单次持续时间最长可达20天之久（梁开学等，2012；夏瑞等，2020）。针对近年汉江水华的连年爆发过程，较多国内外学者对水华的现状及成因做了相关基础研究，如梁开学等（2012）发现2000年以后，汉江中下游硅藻水华与1992—2000年相比，面积缩小，但频率增加；景朝霞等（2019）指出2004—2014年汉江中下游非汛期富营养化严重，非汛期总磷和总氮浓度呈降低和增加趋势；潘晓洁等（2014）发现汉江春季水华暴发期间，不同采样点，硅藻种类占浮游植物比例在44.13%—51.29%之间；夏瑞等（2020）研究了识别了汉江河流型水华多影响要素，指出近50年来，汉江流域气候逐渐呈现暖干化趋势，且藻类密度的变化比营养盐和水文因子更为敏感；吴兴华等（2017）指出汉江硅藻水华发生的原因为适宜的气候条件、高硅氮比和低流量；殷大聪等（2012）发现高营养盐负荷、低流速、小流量、晴好的气候条件是导致汉江水华主要因素。然而，受多种因素制约，这些研究缺乏在连续长时间序列尺度下，对人类和自然双重影响下的河流水华变化特征及多影响要素综合识别。

本研究以汉江武汉段为研究区域，基于STL趋势分解法（Seasonal Trend Decomposition using LOSS）、MK突变检验法（Mann-Kendall）、冗余分析（Redundancy Analysis，RDA）等方法，主要分析了 2004—2017年 2—4月汉江武汉段藻类密度及其影响要素的突变规律，识别了藻类密度突变的主要影响因子，为汉江水华防治和预警预报提供技术支撑。

1 研究区域与方法

1.1 研究区域与数据

汉江发源于陕西省宁强县潘冢山，流经陕西、湖北两省，于武汉市进入长江，全长1577 km，落差1964 km，流域面积1.59×105km2。汉江流域属北亚热季风气候区，气候温暖湿润，夏季平均气温26 ℃，冬季平均气温在 0 ℃间，年平均降水量约为897 mm，降水集中在5—10月，多年平均径流量为5.91×1010m3。汉江武汉段位于平原地区，全长50 km，是武汉市居民用水和工业用水的主要水源（殷大聪等，2017）。汉江武汉段径流缓慢，既受到上游南水北调取水工程影响下的来水来沙条件额度制约，又受到长江顶托作用的影响，具有十分复杂的水文情势，水质及水华问题关注度高（谢平等，2004）。

环境与生物监测数据来源于国家城市供水水质监测网武汉监测站提供的宗关、白鹤嘴与琴断口3个站点的水质监测数据，主要监测指标为总氮（TN，mg·L−1）、总磷（TP，mg·L−1）、藻类密度（107cell·L−1），监测时段为 2004—2017 年 2—4 月每个自然月的上中下旬各10 d左右；这个时段汉江处于枯水期同时水华现象较为频繁（梁开学等，2012；夏瑞等，2020）。在每个采样点，采集1.0  L的表层水加鲁哥氏液固定，静置48 h后采用虹吸方法吸取上清液，留100 mL备用。藻类细胞的数量在24  h内分析测定，吸取0.1 mL滴入计数框，用视野法计数，每个样计数3次，取平均值；采用分光光度法对叶绿素含量进行测定；采用钼酸铵分光光度法和碱性钾消解紫外分光光度法，对1  L水样进行原位过滤并带回实验室进行营养成分分析，包括 TP和TN 等（Xia et al.，2019；国家环境保护总局，2002）。水文监测数据来源于长江水利委员会水文局（http://www.cjh.com.cn/），均为同时期仙桃水文站的流量（Q）、流速（v）、水位（WL1）以及水温（t）。同时采用了汉口水文站的水位数据（WL2），仙桃水文站与汉口站水位差 ΔH。研究区域及观测站点位置、分类见图1。

图1 研究区域及采样点位分布Fig.1 Geographical location and distribution of monitoring station

1.2 研究方法

研究联合使用STL趋势分解法、MK检验法、以及 RDA冗余分析法，旨在综合识别与河流水华暴发密切关联的藻类密度变化趋势和突变点，辨析水华年和非水华年期间藻类密度主要驱动因子和演变规律，其中STL趋势分解法主要识别藻类密度的时间周期变化趋势及突变点，MK检验法则检验变化趋势置信度水平，RDA冗余分析则识别影响藻类密度变化的主要驱动因子及贡献。具体研究方法如下图2。

图2 河流水华拐点突变分析方法Fig.2 Mutation Test Method for the evolution and mutation attribution of water bloom

1.2.1 STL趋势分解法

STL（A Seasonal and Trend Decomposition Using Loess）最早由Cleveland（1979）提出，它是以鲁棒局部加权回归作为平滑方法的时间序列分解方法。STL方法定义T(t)时刻原始数据Yt经STL分解为趋势项Tt、季节项St和残差项Rt：

式中，趋势项Tt被认为是相对稳定的数据变化项，代表了低频数据的变化趋势，是对污染源和气象条件的综合反应；季节性成分St为高频率变化的周期性稳定的数据扰动项，主要是由气象条件的季节变化引起的；残差Rt是不规则的随机扰动引起的。STL方法最早应用于经济学领域特别是对美国失业人口数据分析上，近年在环境科学、生态学、水质科学等学科中有着广泛的应用（Silawan et al.，2008）。

1.2.2 MK突变检验法

Mann-Kendall（MK）统计检验方法由Mann and Kendall两人在20世纪中叶发明并完善，它不需要样本遵从一定的分布，且受异常值影响较小（Mann et al.，1945；Kendall et al.，1975）。MK 统计检验方法在降水、水文、气候变化等时间趋势研究上广泛应用（Evans et al.，2000）。

设时间序列数据（x1,x2,…xn）是n个独立的、随机变量同分布的样本，定义标准的正态统计变量S：

式中，sign()为符号函数。当Xi−Xj＜0、=或＞0 时，sign(Xi−Xj)为−1；当Xi−Xj=0 时，sign(Xi−Xj)为 0；当Xi−Xj＞0 时，sign(Xi−Xj)为 1。

对于统计值Z来说，Z＞0时，上升趋势；Z＜0时，下降趋势。|Z|＞1.28、|Z|＞1.64、|Z|＞2.32 分别表示时间变化趋势通过了置信度α=90%、α=95%和α=99%的显著性检验。

当采用非参数 Mann-Kendall法进行突变检测时，检验统计量与上述Z有所不同。设有一时间序列如下：x1,x2, …,xn，构造秩序列ri，ri含义为xi＞xj（1≤j≤i）的累计数。定义：

Sk均值E(Sk)与方差Var(Sk)的定义如下：

假定一组时间序列是随机独立的变量，则对这组时间序列的统计量UFk作如下定义：

式中，UFk为标准正态分布；给定显著水平α，查询正态分布表，得到临界值Uα。通过统计序列UFk和反序列UBk交叉点对序列x变化趋势进一步进行分析，同时明确突变时间，并指出突变区域。

1.2.3 RDA冗余分析法

20世纪50年代开始生态学家开始基于排序的方法研究植被群落的连续分布。冗余分析（Redundancy Analysis，RDA）是一种使用物种和环境因子组成数据的排序技术，是当前生态学领域较为常用的直接梯度变量约束排序分析方法（Ordination Analysis。Vonwehrden et al.，2009）。冗余分析最早由Vandenwollenberg A L发明，它的主要原理简要如下：

式中，Z1代表独立的内生变量n×r矩阵；Z2代表自变量或外生变量为n×t矩阵；W为t×d权重矩阵；A为d×r载荷矩阵；E为n×r残差矩阵；式（9）为降秩回归模型。

RDA计算结果可以将环境因子和物种因子等样点投射到排序轴构成的二维平面上，并评价各因子之间的相关性和贡献度（赖江山，2013）。

2 结果与讨论

2.1 藻类密度时间变化及突变分析

按照水华发生的统计标准（窦明等，2002），汉江武汉段在2008—2011年连续4年暴发水华事件，藻类密度均达到 1.0×107cell·L−1以上。由图 3 可知，2004—2017年2—4月汉江武汉段宗关、白鹤嘴与琴断口 3个观测站藻类密度均值为 0.37×107cell·L−1；2008年水华事件中藻类密度达到最高值，即 5.3×107cell·L−1；2016 年藻类密度反弹上升，出现一次水华事件。从不同点位水华期间藻类密度比较上看，2008、2011、2016年水华暴发期间3个站点藻类密度峰值基本相近，2009年与2010年水华暴发期间 3个站点藻类密度峰值则呈现出宗关＞琴断口＞白鹤嘴的趋势。

图3 2004—2017年汉江藻类密度（D）时间变化趋势Fig.3 Variations of algae density (D) in Hanjiang River in 2004‒2017

从STL分解的趋势项变化上看（图4、表1），2004—2017年，汉江3点位（宗关、琴断口、白鹤嘴）藻类密度时间变化趋势基本一致，2004—2008年呈现明显的上升趋势，2008—2017年整体为波动下降的趋势。MK检验结果显示，2004—2017年，藻类密度上升趋势不显著，2004—2008年藻类密度呈显著上升趋势（α=0.05），2008—2014年呈显著下降趋势（α≤0.1），2008年前后则是藻类密度变化的主要突变时段。

图4 基于STL趋势分解的2004—2017年汉江藻类密度趋势项时间变化Fig.4 Trends of algae density in Hanjiang River in 2004‒2017 based on STL method

表1 2004—2017年汉江藻类密度MK趋势检验及拐点分析Table 1 Trend test and breakpoint analysis of algae density in Hanjiang River in spring 2004‒2017 based on MK method

从季节变化上看（图3，表2），3个站点的藻类密度变化具有明显的季节性变化，2月中旬至 3月上旬为汉江武汉段水华易发期，且藻类密度呈现双峰现象，其峰值分别在出现在2月中旬与3月上旬，说明这两个时间较2月下旬比，更易发生水华；2008、2010、2016年峰值主要出现在 3月上旬，2009、2011年峰值出现在2月中旬；2004—2008年藻类密度最高值出现在3月上旬，2009—2017年其最高值出现在2月中旬，说明水华暴发时间有前移的趋势。查看汉江中下游历年水华发生情况统计数据，汉江水华主要发生在每年2—3月，主要优势种均为适宜温度较低的硅藻。2月上旬至3月下旬，汉江武汉段水体温度均在0—15 ℃，均值温度分别为6.8、9.0、7.7、9.8、13.0 ℃，且各旬温度的最大值均超过 10 ℃，在营养盐较充足的情况下，硅藻对低温条件下也能有很好的适应性，这与辛小康等（2019）与夏瑞等（2020）研究结果一致，这一定程度上导致水华峰值呈前移趋势。

表2 2004—2017年水华期间不同月旬汉江藻类密度峰值密度统计分布Table 2 Peaks of algae density in Hanjiang River during algae blooms in 2004‒2017 in different month 107 cell·L−1

2.2 水华暴发相关因子突变分析

营养盐的含量、温度的适宜性以及水动力条件是影响藻类生长繁殖的主要影响因素，为解释藻类时间变化趋势和突变点，研究进一步分析了营养要素（TN、TP、N/P）、气候要素（水温t）、水文要素（流量Q、流速v、水位WL1与水位WL2，水位差ΔH）的时间变化趋势与突变点分布。

从营养要素年变化趋势上看（图 5、表 3），2004—2017年2—4月TP处于 0.095—0.128 mg·L−1之间，与藻类密度的相关系数达到0.27（P＜0.05），相关系数较小且相关性不明显，这可能与TP变幅较小而藻类密度波动较大有关；但在所有的因子中TP仍与藻类密度的相关系数最大，且 2004—2017年TP时间变化趋势不显著。硅藻生长所需的TP为0.01 mg·L−1，TP＞0.01 mg·L−1时，硅藻大量生长（Potapova et al.，2007）；TP直接影响硅藻种群的稳定性，水华初期磷酸盐浓度为0.06 mg·L−1，水华中期为 0.05 mg·L−1，而水华末期下降至 0.04 mg·L−1（Soininen et al.，2004）。2004—2017年2—4月TN处于 1.578—2.143 mg·L−1之间，且 2004—2017 年TN 时间变化趋势不显著。TN＜0.2 mg·L−1时，low N硅藻大量生长，TN＞3.0 mg·L−1时，high N 硅藻大量生长，TN 在 0.2—3.0 mg·L−1则为 low N 和 high N硅藻共同生长（Chen et al.，2004）。2004—2017年，汉江武汉段2—4月TN在low N和high N硅藻均可生长的浓度范围，含量充足，这可能是造成藻类密度与TN含量没有相关性的原因。2004—2017年2—4月N/P处于16.344—23.168之间，与藻类密度的相关性系数为−0.28（P＜0.1），且 2004—2017年N/P无明显上升趋势（P＞0.1）。在水华易发生的2月中旬到3月上旬，N/P均大于16，且随着TP的变化出现相反的变化。研究表明汉江武汉段硅酸盐 2月与 4月的含量均在 3 mg·L−1左右（殷大聪等，2017），其含量远远大于韩国Lower HanRiver硅藻水华（优势种同样为冠盘藻）发生时河段中硅的含量（0—1.5 mg·L−1）（Jung et al.，2009）；硅负荷已不再是汉江水华的限制性因素；由此可见，与TN、N/P、硅酸盐浓度相比，汉江武汉段水华的发生更易受到TP含量变化的影响，这与Xia et al.（2020）研究结果一致。

图5 2004—2017年汉江藻类密度及相关因子时间变化趋势Fig.5 Variations of algal density and influencing factors from 2004 to 2017

表3 2004—2017年汉江藻类密度及相关因子变化趋势显著性检验Table 3 Trends of algal density and influencing factors in Hanjiang River from 2004 to 2017 and their significance test with algae density

从气候要素（环境因子）年变化趋势上看，2004—2017年2—4月水温处于9.593—14.437 ℃之间，与藻类密度的相关性系数达到−0.25（P＜0.1），且2004—2017年无明显的时间变化趋势（P＞0.1）。水华（Density＞1.0×107cell·L−1）发生时，水体实测温度范围为10.011—13.031 ℃。汉江硅藻水华的主要优势种一般为冠盘藻，其在低温条件下也能有很好的适应性（夏瑞等，2020）。从水文要素年变化趋势上看，2004—2017年2—4月水位、流量呈现明显的下降趋势，且流量与藻类密度相关系数为−0.21（P＜0.1）。藻类密度与流速与水位均无明显相关性。2004—2017年2—4月水位差范围为7.714—10.056 m，与藻类密度的相关系数达到0.30（P＜0.05）；水位差对藻类密度的影响与汉江回水区的顶托对汉江下游河段的水动力条件变化密切相关（Qu et al.，2014）。由此可见，单独水文要素的变化，并不足以引起汉江水华的发生，流量减少、流速下降，加之长江水位对汉口的顶托作用，易导致汉江武汉段出现水华。

从营养盐含量、水文条件和气候条件的变化趋势与突变点时间分布上看（图 6），TP、TN的 UF和 UB曲线交点主要出现在 2007—2008年和2015—2016年，水温UF和UB曲线交点主要出现在2004—2007年和2014年，流量UF和UB曲线交点主要出现在2004年和2014年，水位差UF和UB曲线交点主要出现在2008—2009年和2014年。殷大聪等（2012）研究表明汉江武汉段TN和硅酸盐等营养盐负荷已经满足春季暴发硅藻水华所需，TP的UF和UB曲线交点出现时间与藻类密度变化拐点（2008年藻类密度最高，2016年水华反弹）一致，表明TP与TN和硅酸盐等营养盐负荷相比，TP限制影响更大。夏瑞等（2020）发现，2008—2011年，汉江中下游主要河段年平均流量和水位均较1961—1999年明显偏少，处于历史上相对较低的一段时期。由此可见，汉江下游武汉段藻类密度的变化趋势是多种因素共同作用的结果，营养负荷、气候要素与水文要素综合变化是汉江武汉段藻类密度2008年前后出现峰值及拐点，而2016年藻类密度反弹上升是出现一次水华事件的主要原因。

图6 2004—2017年汉江环境因子与水文因子变化趋势及突变分析Fig.6 Trends and mutation analysis of environmental and hydrological factors in Hanjiang River from 2004 to 2017

2.3 水华暴发关键驱动因子分析

为对比分析水华年和非水华年影响因子的区别及影响，研究进一步选取 WL1为仙桃水文站的水位数据，代表汉江武汉段的水文情势；WL2数据来源为汉口水文站水位数据，其位于汉口汇入长江口的下游，代表长江的水文情势；水位差ΔH为汉口水文站与仙桃水文站水位差值，表征长江对汉江的顶托作用。

从2004—2017年2—4月汉江下游水华年与非水华年气象水文要素统计结果上看（图7），与非水华年相比，叶绿素浓度偏高14.4%，TN、TP均偏高1.9%左右，水温偏低1 ℃，流速偏低0.05 m·s−1左右，仙桃、汉口平均水位分别偏低0.283、0.266 m，水位差偏小0.171 m，流量偏低15.0%。可见，水华发生流量小、氮磷负荷高且温度偏低的年份。营养盐含量已经不是河流藻类水华的主要限制因素（殷大聪等，2017）；当河流中的总磷含量超过 0.030 mg·L−1时，营养盐含量将不会成为主要限制因素（Sivapragasam et al.，2010）。Liu et al.（2012）研究表明，与湖泊富营养化相比，人类活动导致的河流水动力条件和水循环条件的变化对河流水华影响更大。Mitrovic et al.（2011）发现低流量提高了水的透明度，并创造了一个适当的光和温度条件，以促进藻类生长和繁殖。卢大远等（2000）指出快速水流会改变水平和垂直水交换率加速污染物的扩散和降解，最终降低河流中藻类繁殖率，从而抑制藻类生长。王红萍等（2004）发现长江回水区的顶托对汉江下游河段的水动力条件也有较大影响，使得汉江武汉段水体滞留时间增加，加大了汉江水华发生的几率；当长江水位较高时，假如汉江武汉段流量较小，汉江水位受到的顶托作用明显，小流量对应高水位，当汉江流量较大时，长江顶托作用相对较小。卢大远等（2000）指出，当汉江水流速度≤0.8 m·s−1、水流量≤500 m3·s−1、水温≥10 ℃、pH＞8.0、DO≥12 mg·L−1、色度≥15 度等是汉江硅藻水华指示因子的“预警值”；王红萍等（2004）以藻类密度500×104cell·L−1为汉江春季水华预警值，提出宗关段水华预警流速为0.225 m·s−1。本研究中流量已经明显大于500 m3·s−1；非水华年和水华年的流速均大于0.225 m·s−1；由此可见，低流速、小流量、高水位差等水动力条件易导致汉江硅藻生长和水华形成，近年汉江水华指示因子的“预警值”已经改变；汉江水华现象是多种因素综合作用的结果，单一的限制因素分析已经满足不了水华现象研究的需求。

图7 2004—2017年2—4月汉江下游水华年与非水华年藻类密度及相关因子变化特征Fig.7 Characteristics of algal density and influencing factors in algal bloom year and non-algal bloom year in spring of 2004‒2017 in Hanjiang River

为进一步厘清主要影响因子对藻类密度的影响和贡献，先对水华年与非水华年汉江藻类密度与影响要素数据进行lg(x+1)转换，然后基于RDA方法对其关系进行分析。由图8可知，水华年前两个排序轴对藻类密度变化的解释率为40.9%，其中对藻类密度影响较大的影响要素为总磷与水位差，对叶绿素影响较大的影响要素为总磷，蒙特卡罗置换检验结果发现对汉江藻类变化趋势的影响较为显著的影响要素为水位差（F=7.4，P=0.010）与总磷（F=4.2，P=0.044），其他影响要素不显著；非水华年前两个排序轴对藻类密度变化的解释率为35.1%，其中对藻类密度影响较大的影响要素为总磷、流量与水位，对叶绿素影响较大的影响要素为水位差与总氮，蒙特卡罗置换检验结果发现对汉江藻类变化趋势的影响较为显著的影响要素为流量（F=3.6，P=0.06）与总磷（F=2.8，P=0.088）。由此可见，磷是水华年与非水华年汉江藻类生长的主要限制营养元素，水华年水位差即长江与汉江水文交互作用对藻类密度生长影响显著（Xia et al.，2020），非水华年汉江的流量对藻类密度影响显著。

图8 2004—2017年2—4月汉江下游水华年与非水华年RDA分析排序图Fig.8 Results of RDA during algal and non-algal bloom years in Hanjiang River in 2004‒2017

3 结论

汉江武汉段是武汉市人民的饮用水水源地，水华的发生对饮用水安全造成极大的威胁，研究建立了基于拐点分析的河流水华诊断与归因分析方法。综合采用STL趋势分解法、MK突变检验法、以及冗余分析等多源统计诊断与归因分析方法，重点分析了2004—2017年2—4月汉江武汉段藻类密度的时间变化趋势和突变点，识别了水华年和非水华年期间藻类密度主要影响因子。

（1）汉江武汉段水华年主要发生在2008—2011年和2016年，2月中旬至3月上旬藻类密度最大，均值为 0.79×107cell·L−1，2008 年水华事件为近 20年最严重。2004—2017年2—4月汉江武汉段宗关、白鹤嘴与琴断口 3个观测站藻类密度均值为0.37×107cell·L−1，峰值密度为 5.3×107cell·L−1；三站点藻类密度差异较小且时间变化趋势一致，2004—2008年呈显著上升趋势，2008—2014年呈显著下降趋势，2008年前后为藻类密度突变点。

（2）营养负荷、气候要素与水文要素综合变化是导致 2004—2017年春季汉江水华现象的主要原因。近年汉江水华指示因子的“预警值”已经改变；总磷是水华年与非水华年汉江藻类生长的主要限制营养元素；低流速、小流量、低水位差等水动力条件易导致汉江硅藻生长和水华形成；与非水华年相比，总氮、总磷均偏高1.9%左右，水温偏低1 ℃，流速偏低0.05 m·s−1左右，流量偏低15.0%。

（3）针对汉江水华的治理，建议除深入分析水华特征和影响因素外，应开展汉江总磷源解析，明确总磷来源，制定总磷污染防控对策；监控汉江武汉段的流量与汉口的水位，加强中上游丹江口水库生态水量调度，预防“顶托”的形成。