河流水质受到地形、降雨、景观、人类活动等因素的影响，普遍存在时空变异。河流水质的时空变异指在一定的区域内，不同时间不同河段的水质指标存在明显的差异性。河流水质时空变异研究是根据观测数据，分析反演水质指标参数的时空分化特征，探究参数自身及其影响因子的时空关系。如李跃飞等^[1]以秦淮河为研究对象，发现污染物浓度表现为旱季高于雨季、市区流域高于郊区流域的时空分布特征。张彦等^[2]比较了暴雨前后河南北部河流水质变异特征，发现暴雨前河流水质指标受到气象因子和河流水体内源的影响较大，暴雨后城市污染源对总氮的绝对贡献率最大。了解水质的时空变异并开展水质评价，对于污染发生的热点区域及时间识别、水污染高效处理、水资源科学管理具有重要意义^[3-11]。

由于受实际条件限制，水质观测在时间和空间上不可能连续，因此常需要结合时空特征分析方法全面探究水质时空变异规律。目前，常用的水质时空特征分析方法有箱线图法^[12]、主成分分析法^[12-13]、单因子指数法^[14-16]、内梅罗指数法^{[15, 17]}、灰色关联分析法^{[15-16, 18]}等，这些方法能较好地反映数据的离散分布情况、识别水质的主要影响指标、判别水质的基本变化趋势，但在实际应用中仍存在一系列问题。比如，如何量化水质变异时空次序？如何判别水质整体达标情况？如何更准确知道河流水质的背景值、最大值以及增长系(指)数？因此亟需探索一种新的简单、有效、准确的水质时空特征分析方法来解决以上问题。

水文经验频率法始于19世纪末，美国的科学家首先应用了频率曲线(当时称为历时曲线)开展水文研究。该方法分析某随机水文变量出现不同数值可能性的大小，在此基础上为相关水资源工作提供具有概率意义的水文序列数据^[19]。水文经验频率法因具有客观的分析过程，能将水文信息以统计规律自动反映出来，在操作上只有线型选择和参数估计两个步骤，操作简单可靠，结果客观清晰，因此被广泛地应用于防洪工作中^[20]。本文根据经验频率法原理，结合河流水质特征，将此方法引入水质的时空特征分析中，以期为河流水质时空变异分析提供简单可行的新思路。

为深入探讨“水文频率-水质”拟合曲线法的优势并科学识别水质变异规律，本文以秦淮河水质监测数据为例，利用箱线图法和“水文频率-水质”拟合曲线法对秦淮河水质时空分布变异特征进行对比分析，旨在为水质精准管理提供科学依据。具体来说，本文的主要目的如下：(1)对比箱线图法和“水文频率-水质”拟合曲线法的分析结果，并讨论“水文频率-水质”拟合曲线法的主要优势；(2)确定本研究中“水文频率-水质”拟合曲线在时间和空间上的最佳拟合形式及其深层次意义；(3)通过“水文频率-水质”拟合曲线法，表征秦淮河水质时空变异特征及不同等级的水质特征值状况。

1 材料与方法 1.1 研究区概况

秦淮河为江苏省南京市的“母亲河”，有南北两个源头，北源句容河发源于句容市宝华山南麓，南源溧水河发源于南京市溧水区东庐山，两河在南京市江宁区汇合成为秦淮河干流。秦淮河流域(31°35′~32°07′N，118°43′~119°18′E)从北源至三汊河全长110 km，全流域面积2631 km²，整体地势从南向北倾斜，大多地区高程低于50 m，上游坡度和扇面大，中下游坡度缓，从上游到下游由传统稻作农业区过渡到高发达城市区。全年雨季主要分布在6-9月，季风特征明显。

1.2 水质数据监测与划分

在秦淮河河段布置采样点(图 1)，采样时间为2010年6月至2013年4月，每月采样1次，14个采样点(z1~z14)均位于秦淮河干流的范围内，用于监测秦淮河沿岸氮、磷浓度的变化。采样点东至句容市东北方向的句容水库入口(z1)，西至南京市江宁区的外秦淮河(z14)。采样后将水样放于冷藏保温箱中储存，带回实验室先通过无菌的0.45 μm玻璃纤维滤膜过滤，然后测量总氮(TN)、总磷(TP)、硝态氮(NO₃^--N)和氨氮(NH₃-N)浓度。NO₃^--N和NH₃-N浓度采用流动分析仪测定(Skalar，荷兰)；TN和TP浓度先用过硫酸盐消解，然后用分光光度法测定硝酸盐和可溶性活性磷，TN和TP消解的回收率分别为98.4%和99.7%^[21]。

由于不同地区和不同时间影响河流环境的因素不一，导致河流水质浓度在时空上表现为不同的量值。本研究结合秦淮河流域实际情况，以6-9月为丰水期(简称“F”)，其他月为枯水期(简称“K”)来进行时间上的划分。同时，根据前期研究，基于人类活动强度，将z1~z3、z5~z8划为低活动区(简称“L”)采样点，其主要受到传统农业和农村居民区污染的影响；z4、z9~z14划为高活动区(简称“H”)采样点，其受到高集约农业和高发展水平城镇的影响^[21]。

1.3 箱线图法

箱线图是统计学中一种常见的统计数图，可以不受异常值的影响，以一种相对稳定的方式表现出数据的离散分布情况^[22]。箱线图主要展示6个关键信息：最大值、上四分位数(P₇₅)、中位数(P₅₀)、下四分位数(P₂₅)、最小值和平均值。箱体部分由2个四分位数构成，最大值和最小值则构成箱体上下2条端线。

为探索高/低活动区各月份水质变异特征和丰/枯水期各采样点水质变异特征，本文通过箱线图法进行初步展示，该分析用Origin 2021软件完成。

1.4 “水文频率-水质”拟合曲线法

经验频率法是将变量按顺序排列，根据数学期望理论建立累积经验频率，从而计算各变量值出现的频率，其计算公式如下：

$ p=\frac{m}{n+1} \times 100 \% $

(1)

式中，p为经验频率，m为样本序号，n为样本总数。式中经验频率的估算主要在于对样本序列中的每一项估算其对应的频率。首先，将总体随机地分为许多个样本，每个样本中都有相同个数的数据，并将其排序；然后，在每一个样本中取出同序项来研究，并计算这些同序项出现概率的均值；最后，当样本数为无穷大时，通过计算即可得到公式(1)。

此方法操作简单，理论可靠，因此在我国水文计算上广泛运用^[23]。本文以氮磷浓度为研究对象，将此方法的思路引入水质的时空评价中，其主要是根据水质样本经验点和最佳曲线形式求算拟合曲线的参数和方程式，是通过水质样本的经验分布探求总体分布的一种新方法。

“水文频率-水质”拟合曲线法的具体步骤如下：

1) 将水质数据从小到大排序，并给予每个水质点1，2，3，……，n的样本序号；

2) 用公式(1)求算每个采样点的经验频率p；

3) 以p为横坐标，水质浓度为纵坐标做散点图，并进行曲线拟合；

4) 通过对拟合曲线形式进行筛选，选择简单且拟合度最佳的曲线形式开展分析；

5) 按照拟合的方程式及其参数，对水质基本特征、水质等级特征值(背景值、最大值及水质环境质量标准特征值等)和水质达标情况进行评判。

通过在Origin 2021软件中不断尝试，发现线性曲线和指数曲线分别为水文频率与水质在时间和空间上的最佳拟合曲线形式，即：

在时间上分别计算高/低活动区各月份氮磷浓度平均值，通过以上方法，选取拟合度最佳的线性曲线形式进行拟合：

$ y=k x+b $

(2)

在空间上分别计算丰/枯水期各采样点氮磷浓度平均值，通过以上方法，选取拟合度最佳的指数曲线形式进行拟合：

$ y=h \mathrm{e}^{a x} $

(3)

式中，x为自变量，即水文频率；y为因变量，即水质浓度；k、b、h和a为系数，e为自然常数。

此方法在水文经验频率法的基础上进行改进，由于环境领域对水质的重要性主要体现在较短时间(如一年、一季度等)的最大或平均值上，同时水质浓度存在周期性、可消减性、空间分布变化性强等特点，这不同于水文中最关注的多年径流量的理念，因此“水文频率-水质”拟合曲线的适用条件以及频率的意义也有变动：

1)“水文频率-水质”拟合曲线无需大量长时间基础数据，但需保证采样点、采样时间和采样频率足够典型，以便能反映河流在时间和空间上的真实水质浓度；同时需对原始数据进行简单处理(如平均值计算、去除异常点等)，从而保证数据的可用性和客观性。而在水文频率中需统计大量年径流量数据，其重点在于年径流量总值。

2)“水文频率-水质”拟合曲线用于水质时间分布特征分析时，其频率的意义为：河流水质共有频率对应百分比的时间小于频率对应浓度值，曲线上的每一个点都代表河流水质时间分布情况。如拟合曲线上某点对应的频率和TN浓度分别为0.8和5 mg/L，则说明该水体有80%的时间TN浓度小于5 mg/L。

3)“水文频率-水质”拟合曲线用于水质空间分布特征分析时，其频率的意义为：河流水质共有频率对应百分比的水体断面小于频率对应浓度值，曲线上的每一个点都代表河流水质空间分布情况。如拟合曲线上某点对应的频率和TN浓度分别为0.8和5 mg/L，则说明该水体有80%的水体断面TN浓度小于5 mg/L。

由于“水文频率-水质”拟合曲线具有上述特点，因此可根据行业水质标准或科研需要对水质关键界限点进行量化，确定背景值、最大值、不同等级标准下水质的特征值等信息，最终实现河流水质的科学管理。

2 结果与分析 2.1 高/低活动区水质时间变化箱线图

由于降雨、温度等气象因素的影响，在不同月份河流氮磷浓度呈现差异。由图 2可知，氮磷浓度在8月前后较低，1月前后较高，但NO₃^--N在各个月份的变异不如其他指标明显。从浓度大小来看，高活动区浓度整体高于低活动区。在高活动区，TN浓度大部分时期都高于Ⅴ类水标准(2 mg/L)，而在低活动区，高于Ⅴ类水的频次虽仍较多，但相比高活动区有所改善；对于NH₃-N，在高活动区大部分时期高于Ⅲ类水标准(1 mg/L)，而低活动区则大部分时期低于Ⅲ类水标准；在高活动区TP浓度在Ⅲ类水标准(0.2 mg/L)上下波动，而低活动区则大多时期低于Ⅲ类水标准。

2.2 丰/枯水期水质空间变化箱线图

由于人为活动强度等因素的影响，在不同区域河流氮磷浓度呈现不同差异。由图 3可知，氮磷浓度在传统农业区z1~z3较低，城市区z4点突增，之后逐渐下降直到进入城市区后又逐渐升高。从浓度大小来看，枯水期浓度整体高于丰水期，但在z1~z3处并不明显。TN不论在丰水期还是枯水期，大部分采样点浓度高于Ⅴ类水标准(2 mg/L)；NH₃-N在丰水期浓度大多位于Ⅲ类水标准(1 mg/L)之下，枯水期则大多位于Ⅲ类水之上；TP则不管在丰水期还是枯水期，大部分采样点浓度都低于Ⅲ类水标准(0.2 mg/L)。

2.3 基于“水文频率-水质”拟合曲线的河流水质时间特征

根据“水文频率-水质”拟合曲线，得到各氮磷指标在高/低活动区的时间特征。如图 4所示，所有R²均在0.82以上，说明拟合方程可信度高，拟合曲线可靠。

通过分析，得到高/低活动区水质在时间上的不同等级特征值与特征参数：拟合曲线的截距b为经验频率x等于0时y的值，可视为河流在12个月份内氮磷浓度背景值；k为增长系数，其表示因变量(河流氮磷浓度)随自变量(频率)增加而增加的快慢程度，k值越大，河流氮磷浓度增加越快；当x等于1时y的值，可解释为一年内河流氮磷浓度最大值。以TN为例，TN_H和TN_L拟合曲线的b值分别为2.22和1.49 mg/L，表明在时间尺度上高活动区TN背景浓度高于低活动区。k值分别为5.57和2.29，表明高活动区TN浓度增长趋势快于低活动区。TN浓度最大值分别为7.78和4.39 mg/L，说明TN在高活动区最大浓度明显高于低活动区。

由图 4可知，NH₃-N在低活动区的背景值大于高活动区，TP在低活动区的增长系数略大于高活动区(相差很小)，其他情况氮磷指标的背景值、最大值和增长系数均为高活动区大于低活动区。

2.4 基于“水文频率-水质”拟合曲线的河流水质空间特征

根据“水文频率-水质”拟合曲线，得到各氮磷指标在丰、枯水期的空间特征，且所有氮磷指标R²均在0.87以上(图 5)，拟合曲线可靠。

由拟合曲线方程式可得，丰/枯水期水质在空间上的不同等级特征值与特征参数如下：拟合曲线的h为河流在某时期氮磷浓度背景值，a为增长指数，x等于1时y值为河流氮磷浓度最大值。同理，以TN为例，TN_F和TN_K背景值分别为1.04和1.62 mg/L，a值分别为1.78和1.83，TN浓度最大值分别为6.12和10.1 mg/L。由图 5可知，除了NH₃-N的增长指数在丰水期略大于枯水期之外，其他氮磷指标的背景值、最大值和增长指数均为枯水期大于丰水期。

3 讨论 3.1 高/低活动区水质时间特征分析

在高/低活动区以各月份氮磷浓度平均值为对象进行拟合，得到其最佳拟合方式为线性曲线形式。由线性曲线的性质可知，其因变量随自变量的变化而呈匀速变化趋势，这说明氮磷指标在时间上的变异是均匀的，一般不会发生水质的突变。通过各氮磷指标拟合曲线(图 4)可知，NO₃^--N在2个活动区的增长系数大致相同，其差异主要体现在高活动区背景值大于低活动区，因此在治理工作中相关部门要着重开展对背景值的控制。TP和NO₃^--N的变异特征较为相似，区别在于TP在低活动区的增长系数略大于高活动区，但相差很小。对于NH₃-N，其在低活动区的背景值略大于高活动区，可能是因为在NH₃-N排放最少的时期内，以传统农业为主的低活动区土壤中的NH₃-N遗留量多于以城市为主的高活动区土壤中的NH₃-N遗留量，因此传统农业土壤中NH₃-N的缓慢释放引发了以上现象；同时，由于高活动区增长系数明显大于低活动区，高活动区NH₃-N浓度在大部分时期内大于低活动区，因此需重点调查高活动区NH₃-N随时间快速增长的原因。而对于TN，高活动区的背景值和增长系数均大于低活动区，这就要求管理者在背景值和增长量上同时给予关注。以上结果表明，一方面高活动区水体氮磷浓度在时间上的分布与低活动区存在较大差异，另一方面不同氮磷指标的时间分布特征也各不相同。

同时，根据《地表水环境质量标准》，运用“水文频率-水质”拟合曲线方程，可解决箱线图法不能对水质关键界限点进行量化的问题。分析发现TN在高活动区全部时段内均未达到Ⅴ类水标准(2 mg/L)；在低活动区有17.61%的时段达到Ⅴ类水标准。NH₃-N在高活动区有24.39%的时段达到Ⅲ类水标准(1 mg/L)；在低活动区有75.38%的时段达到Ⅲ类水标准。TP在高活动区有45.63%的时段达到Ⅲ类水标准(0.2 mg/L)；在低活动区有90.93%的时段达到Ⅲ类水标准。

由前文分析知，高活动区水质明显劣于低活动区，这与很多学者的研究结果相同^[24-28]。在高活动区，工农产业集中、人口密集、交通发达，属于城镇区或城郊集约农业区。在高活动区往往伴随着污染物的强排放，且该区域建设用地密集，大量的不透水面导致了氮磷营养物快速迁移至河流，污染物来不及被截留或消纳就进入水体，因而易导致水质恶化^[29]。朱珍香等^[24]在厦门后溪流域的研究发现，城市的扩张以及污水处理设施的不足均会导致水体中污染物浓度增加；廖炜等^[30]以斧头湖流域为对象，发现在城镇化发展过程中，人畜、耕地及城镇生活源是面源污染的主要来源；张亚娟等^[31]以城镇化趋势明显的于桥水库流域为例，发现城镇居民用地与水质指标存在极显著的相关性，说明居民用地是流域水质污染重要的贡献源。而在低活动区，主要以传统稻作农业为主，人口相对稀少，污染排放量低，且污染物在排放至水体后会被消纳，因而相较于高活动区其水质较佳。

3.2 丰/枯水期水质空间特征分析

在丰/枯水期以各采样断面氮磷浓度平均值为对象进行拟合，得到其最佳拟合方式为指数曲线形式。由指数曲线的性质可知，其因变量随自变量的变化而可能呈现快速变化趋势，结合本研究结果易知，氮磷浓度随x的增加呈指数级增长趋势，因此在空间上可能发生水质的突变。通过各氮磷指标拟合曲线(图 5)可知，NO₃^--N和NH₃-N丰、枯水期的增长指数大致相同，其差异为枯水期背景值明显大于丰水期，因此相关部门需在枯水期加强对背景值的控制；而对于TN和TP，枯水期的背景值和增长指数均大于丰水期，这就要求管理者在背景值和增长量上同时给予关注。以上结果表明，一方面枯水期水体氮磷浓度在空间上的分布与丰水期存在较大差异，另一方面不同氮磷指标的空间分布特征也各不相同。

同时，根据《地表水环境质量标准》，结合“水文频率-水质”拟合曲线方程，分析发现TN在丰水期有37.06% 的河段达到Ⅴ类水标准(2 mg/L)；在枯水期有11.59%的河段达到Ⅴ类水标准。NH₃-N在丰水期有77.44% 的河段达到Ⅲ类水标准(1 mg/L)；在枯水期有37.90%的河段达到Ⅲ类水标准。TP在丰水期有79.89%的河段达到Ⅲ类水标准(0.2 mg/L)；在枯水期有61.99%的河段达到Ⅲ类水标准。

由以上分析知，河流的大部分河段丰水期水质优于枯水期，结合实际情况分析，虽然在丰水期经过降雨径流冲刷可将地表各类营养物带入水体，但由于水量增加，可能对河流氮磷浓度起到了稀释作用；另一方面，由于夏季温度较高，适合各类植物和微生物的生长，增强了水体内植物的吸收和微生物的反硝化作用^[32]等，从而进入水体的营养物可以被快速消纳，最终使得水体氮磷浓度降低，金鑫等^{[8, 33-34]}在研究中也发现了相似的现象。但项颂等^[35-37]研究却得出了相反的结果，其发现在丰水期水质浓度会明显升高，这可能是因为其研究区域内丰水期受面源污染影响严重，或研究的流域范围较小，污染物来不及消纳便进入水体从而导致河流水质浓度升高。因此，在开展水质的分布特征研究中，需结合实际下垫面、降雨及人类活动情况具体探讨^[38]。根据以上分析并结合实际情况，虽然蓝藻水华等恶性水体事件常发生在雨量充沛的丰水期，但枯水期的高氮磷浓度仍导致水体存在较高的蓝藻水华事件暴发风险，因此也必须给予密切关注。

3.3 “水文频率-水质”拟合曲线法优缺点分析

由分析可知，“水文频率-水质”拟合曲线法和箱线图法同样操作简便，结果清晰直观，适用于水质时空变异特征的初步探索。与此同时，“水文频率-水质”拟合曲线法还具有以下特有优势：

1) 结果直观有规律。从“水文频率-水质”拟合曲线中可直观地看出高活动区氮磷浓度大于低活动区(图 4)，枯水期氮磷浓度高于丰水期(图 5)。相较于箱线图法，其可以更好地将水质时空特征按次序有规律地展示出来，同时生成曲线方程，有助于科研和政府部门更加深入地认识河流水环境。

2) 可量化关键水质评判点。结合地表水环境质量标准，“水文频率-水质”拟合曲线法可将各类水质的界限值表现出来，从而得到河流水质整体达标情况。如在高/低活动区时间特征“水文频率-水质”拟合曲线中(图 4)，通过拟合方程的简单计算，就可得出“TN在高活动区全部时间内均未达到Ⅴ类水标准(2 mg/L)；在低活动区有17.61%的时间达到Ⅴ类水标准”的推论。而在对应的箱线图(图 2)中却无法实现这个目标。

3) 可展示水质特征值及特征参数，使水质时空变异的过程更为清晰。根据“水文频率-水质”拟合曲线法的性质，可从拟合曲线方程中推算河流水质的背景值、最大值以及增长系(指)数。如在丰/枯水期空间特征“水文频率-水质”拟合曲线中(图 5)，TN_F和TN_K背景值分别为1.04和1.62 mg/L，增长指数值分别为1.78和1.83，这说明枯水期水体的基础污染负荷强度高于丰水期，同时由于枯水期增长指数更大，故枯水期TN浓度也可能以更快的速度增长，相关部门应给予重视，其他指标同理。以上分析使得河流水质变异过程更加清晰。

虽然“水文频率-水质”拟合曲线法具有以上诸多优点，但得出的结论仍存在一定的不确定性。如使用“水文频率-水质”拟合曲线开展河流水质探讨时，可能由于水质采样点、采样时间和采样频率代表性不足，或未对原始数据进行预处理(求平均值、去除异常点等)，将导致在某些区间内拟合值和实际值发生偏离，从而增加了对河流水质真实值估计的不确定性。因此，在初期河流采样布点时，需全面考虑其采样点的典型性，即采样点是否布设于顺直河段、河床稳定、水流集中、无浅滩处，且避开死水区、回水区、排水口等情况。同时，还要充分关注采样时间与采样频率是否合适，如每次的采样时间是否具有代表性；采样前后是否发生特殊事件(大型降水事件、灌溉事件、污水偷排等)并做好记录；在水质浓度变化较大的月份是否增加采样频率等等。最后，对得到的原始数据要进行一定的预处理，如求各月份的平均值或求各采样点的平均值后再进行方法的调用，以防异常数据干扰拟合结果。此外，由于“水文频率-水质”拟合曲线是根据统计学概率将现实水质分布反映出来，现实中的小概率事件可能无法客观地体现在拟合曲线中，因此要理性辩证地认识拟合曲线的参数与水质变异过程。

4 结论

本研究采用箱线图法和“水文频率-水质”拟合曲线法对秦淮河干流水质变化趋势进行探索，根据流域特点划分高/低活动区和丰/枯水期开展研究分析，主要结论如下：

1) 相较于箱线图法，“水文频率-水质”拟合曲线能够更加直观有规律地将水质整体情况进行展示，同时也可以量化关键水质评判点与特征值，使水质时空变异的过程更为清晰。

2) 高/低活动区各月份氮磷指标“水文频率-水质”拟合曲线的最佳形式为线性曲线，丰/枯水期各采样点氮磷指标“水文频率-水质”拟合曲线的最佳形式为指数曲线。

3) 河流氮磷浓度在时间和空间上均受到背景值和增长系(指)数的影响，但在时间上一般不会发生突变，而在空间上则有较大可能发生水质的突增。

4) 氮磷指标在时间和空间上均表现出不同差异。高活动区水体氮磷浓度高于低活动区，是水质恶化的热点区域；枯水期水体氮磷浓度高于丰水期，是水质控制的关键时期。