随着我国经济社会快速地发展，城镇化快速推进和工农业污染排放总量激增，东部发达地区浅水湖泊的水环境形势更加严峻，湖泊水体近年来黑臭现象频繁发生^[1-3].特别是太湖、巢湖等水体中黑臭问题尤为严重，成为当地政府与民众关注的热点.

湖泊的黑臭主要由有机物大量快速聚集、分解和水体严重缺氧引起^[4-5]：大量暴发、聚集的藻类死亡后其有机残体在适宜的气象、水文条件下发生厌氧分解消耗水体中的溶解氧（DO），使水体和沉积物快速进入厌氧状态；这导致水-沉积物界面附近发生还原反应，大量高价态的Fe、Mn和S被还原成低价态的Fe²⁺、Mn²⁺和S^2-并向表层水体释放，此过程中形成的FeS是主要的致黑物质^[6-8]，挥发性有机硫化物则是主要的致臭物质^[9-11]，Fe²⁺和S^2-浓度的变化对黑臭形成和发展至关重要^[12-14].此前黑臭水体定量关系研究多以氮磷等营养元素以及化学需氧量、生化需氧量等化学计量的变化为主，研究方法则主要以普通的线性回归为主^[15-16]，这显然不足以定量揭示黑臭现象的变化规律.此外，这些研究完全忽略了湖水流动特性对黑臭的影响，相关方面的研究较少.

水体流动对黑臭现象可能造成的影响分为3种：一是相同空间不同地点的DO、Fe²⁺和S^2-等物质通过水体流动而发生联系，使得这些物质的变化存在空间相关性；二是水体流动对污染物的分散作用，可能改变Fe²⁺、S^2-等物质在空间上的具体分布；三是流速会影响水中DO的复氧过程，进而影响水中Fe²⁺、S^2-等黑臭相关物质的浓度.总之，水体流动对黑臭现象的影响是重要而复杂的.

全面系统研究水体流动、DO、Fe²⁺、S^2-等对黑臭现象的影响，除实验观测等方法以外，环境生态水力学模型是运用较多的数学模型.这类模型虽然能够较准确地考虑水体流动、湖底地形、环境因子等要素的空间和时间变化对黑臭现象的影响，但建立这种模型需要多年累积的水文、气象和环境的观测数据，以率定和验证模型中涉及的几十甚至上百个参数^[17-19]；因此，在太湖^[18-21]、长江^[22-23]等有长期观测数据的水域中，较多使用这类模型研究营养盐循环及富营养化问题.此外，环境生态水力学模型是基于机理的模型，运用此类模型需要明晰黑臭水团物质与各种环境因子相互作用的机制，并将这些机制准确地定量化；而且求解模型需要专门的大型程序配合，这种程序的编制工作量非常巨大.相比之下巢湖的观测数据略有不足，巢湖黑臭现象发生的环境机制研究也不够成熟，难以运用环境生态水力学模型研究巢湖的黑臭现象.目前，巢湖的数学模型多是生态^[24-25]及水华暴发预测模型^[26-27]的研究，与水动力结合的环境数学模型研究较少，仅有的环境生态水力学案例大都集中于研究“引江济巢”工程对巢湖^{[17, 28]}的影响.相对而言，空间计量模型能在一定程度上考虑湖水流动对环境变量的影响，这种方法需要的数据较少，求解过程也没有环境水力学模型那么复杂，非常适用于对黑臭现象空间分布的初步探索性研究^[29-30]，在环境污染的研究领域也有成功的应用^[31].

本文使用空间计量模型，以巢湖南淝河口汇流湾区黑臭暴发时的情况为例，研究流速对该区域DO、Fe²⁺、S^2-浓度和分布的影响，结果可以为黑臭治理提供理论参考.

1 空间计量模型

相较于假设变量相互独立的统计学方法，空间计量模型能够考虑变量之间在空间分布等方面的关系；适合于探索环境变量的空间分布规律.

1.1 空间相关性度量

变量之间空间相关性有多种指数可以度量，如Moran's I指数、Geary's比率等方法.其中，Moran's I指数由于其物理清楚，计算相对简单而得到广泛应用.若随机变量X的具体观测值为x₁、x₂、……、x_n，则其自相关Moran's I指数计算方法由公式(1) 计算：

$\begin{split} I=\frac{n∑\limits^n_{i=1}∑\limits^n_{j=1}w_{ij}(x_i－\bar x_)(x_j－\bar x_)}{∑\limits^n_{i=1}∑\limits^n_{j=1}w_{ij}∑\limits^n_{i=1}{(x_i－\bar x_)}^2} \end{split}$

(1)

式中，n为观测值的总数， $\bar x$ 为流速、DO、Fe²⁺及S^2-浓度的平均值，x为相应的观测值， $w_{ij}$ 为空间权重矩阵W的第i行、j列的元素.

空间权重矩阵W的存在是Moran's I指数与普通相关系数的重大区别；空间权重矩阵用于衡量同一变量的采样点分布在纯粹地理空间上的相关性.生成空间权重矩阵需要先生成空间邻接矩阵用于表达不同采样点之间的邻居关系.有很多不同的邻接矩阵生成方式, 如传统的“后”式或“车”式邻接方法生成；也有较为先进的K-Nearest方法生成.生成邻接矩阵后只需要指定权重方式，即可生成空间权重矩阵，具体过程可参考文献[29].在此考虑观测点分布并不是完全规则，所以用K-Nearest方法生成空间邻接矩阵，邻接点个数为4，经过实验邻接点个数对空间权重矩阵的结果影响不大.

Moran's I指数与普通的相关系数一样，取值范围为[-1, 1]，但并非越接近于0其相关性就越小.正确的解读方式：若相邻单元的环境变量相关性越差，则计算的Moran's I指数越接近于随机分布的Moran's I指数.环境变量的n个观测值若为空间随机分布，则Moran's I指数表达式为：

$\begin{split} I_{\text{random}}=－\frac{1}{n－1} \end{split}$

(2)

若两个环境变量X、Y对应的观测值分别为x₁、x₂、…、x_n和y₁、y₂、…、y_n，则它们之间的互Moran's I指数由公式(3) 计算：

$\begin{split} I_{ij}=\frac{n∑\limits^n_{i=1}∑\limits^n_{j=1}w_{ij}(x_i－\bar x_)(y_j－\bar y_)}{∑\limits^n_{i=1}∑\limits^n_{j=1}w_{ij}∑\limits^n_{i=1}{(x_i－\bar x_)}^2} \end{split}$

(3)

Moran's I指数的统计显著性不能直接计算，必须用自助法(Bootstrap)抽样计算得到.

1.2 空间回归模型

考虑变量空间相关性之后，一般的二阶空间回归模型可以表达为:

$\begin{split} Y=ρW_1Y+Xβ+U\\U=λW_2U+ε\\ε\sim N(0,σ^2) \end{split}$

(4)

式中， $ρ、β、λ$ 为待定系数；W₁和W₂分别为因变量Y和残差U的空间权重矩阵，它们可以相同也可以不同.根据情况的不同，这一模型又包含了3大类空间模型，即一阶空间自回归模型：

$\begin{split} Y=ρW_1Y+ε\\ε\sim N(0,σ^2) \end{split}$

(5)

空间滞后模型：

$\begin{split} Y=ρW_1Y+βX+ε\\ε\sim N(0,σ^2) \end{split}$

(6)

空间误差模型：

$\begin{split} Y=Xβ+U\\U=λW_2U+ε\\ε\sim N(0,σ^2) \end{split}$

(7)

对于具体环境变量的观测值，可以通过拉格朗日算子法计算最适合，且统计上也显著成立的空间回归模型.对于具体的空间回归模型，需要使用极大似然法估计其中的系数^[29-30].

以上2种空间计量模型在研究中的作用有所不同.相关性度量衡量的是环境因素之间总体的空间相关程度，空间回归模型则侧重于环境因素之间局部的相关关系.

2 数据来源

数据来源于2013年5月河海大学和中国科学院南京地理与湖泊研究所在巢湖南淝河口发生黑臭现象时，联合采样所得数据及水样分析结果.采样点以南淝河口为中心，呈扇形分布（图 1）.由于南淝河是该区域主要的污染源，在表观上河口区黑臭较严重，所以在南淝河口附近对采样点进行加密.

在每个监测点同时采集水样, 同于测定流速、溶解氧以及Fe²⁺和S^2-指标.流速采用美国YSI公司的YR-1型声学多普勒流速仪(Flow Tracker)测定；DO浓度采用日本HORIBA公司U-53型多参数水质分析仪监测；S^2-浓度在实验室采用对氨基二甲基苯胺光度法测定^[32]；Fe²⁺浓度按Tamura等^[33]改进的邻菲啰啉光度法测定；共得到62组有效数据；采样结果的主要统计特征见表 1.

采样时该区域肉眼即可观察到明显的黑色水团，并伴有强烈恶臭气味；说明该地区正处于黑臭现象较为严重的时期；因此，所得数据可用于研究巢湖湖湾区黑臭水体的水质状态.

3 结果与分析

Fe²⁺、S^2-、DO和流速的空间相关Moran's I指数及其显著性检验结果见表 2.其空间权重矩阵以距离倒数为空间权重，显著性检验结果则由自助法计算得到.

若S^2-、Fe²⁺、DO和流速这4个量的观测值在空间上完全随机，即无空间相关性，则公式(2) 中的n取62，可得对应的Moran's I指数为-0.016. 4个量的自相关系数远远高于该值，其显著性检验结果也表明在95%的概率上这4个量的分布在空间上自相关性显著. S^2-、Fe²⁺、DO和流速在某采样点的值与其周围采样点的值有较强的正相关性，这说明由水体流动引起物质对流扩散作用，导致水质等物理化学量的空间相关显著.因此，在研究黑臭现象及黑臭水体的运动规律时，不能孤立研究暴发黑臭的样点处水体变化规律，必须将这种空间相关性纳入考虑，同时需要关注黑臭水体区域整体的水体水质.

与空间自相关Moran's I指数的结果不同，Fe²⁺、DO和流速的二元Moran's I指数结果表明，Fe²⁺与DO和流速均存在着显著负相关（表 2），即Fe²⁺浓度不仅与采样点本身的DO和流速有关，而且与其附近的DO和流动情况有关.这是由于水体的流动不断改变黑水团的位置，而DO与Fe²⁺之间的化学反应则相对较慢，所以湖湾区某点Fe²⁺浓度情况一定与黑水团从南淝河口至采样点的整个漂移路径上DO的分布有关，不仅只受采样点本地DO浓度的影响. Moran's I指数的结果也表明，Fe²⁺浓度受其周围DO浓度的影响要高于周围流速对其的影响. Fe²⁺与S^2-浓度呈显著正相关，这是由于2种离子均是引起黑臭的主要原因，因此在黑臭水体中Fe²⁺和S^2-中一种离子的浓度较高可能伴随着另一种离子浓度也很高.

S^2-浓度的相关性与Fe²⁺浓度的情况有所不同：S^2-浓度与流速也呈显著负相关，但其与DO浓度之间不存在相关性.这是由于水体流动对水中的离子主要是物理作用，流动对S^2-浓度的作用与其对Fe²⁺浓度的作用类似；而生成S^2-和Fe²⁺所需的缺氧还原条件不同，因此造成它们与DO浓度的相关性显著不同.增加DO浓度对去除Fe²⁺效果比去除S^2-的效果好.

通过流速和DO浓度之间的显著性检验，在95%置信水平上不能发现这两者之间存在关系.虽然理论上流速会影响DO的复氧速率，但由于湖湾区域污染较严重，水中耗氧物质的浓度对DO浓度的影响占绝对主导地位；再加上采样时天气情况良好，湖流流速绝对值较小.湖流的最大速度仅为10.9 cm/s，绝大多数采样点的流速小于10 cm/s（表 1）；因此实际上湖流的复氧作用不强，其对DO浓度的影响较小.这也进一步说明湖湾区的湖流对Fe²⁺浓度分布的影响主要是通过分散作用，而并非是通过对DO来影响Fe²⁺浓度.

为更好地分析各变量之间的空间分布关系，可利用观测值分位点对采样点的DO、Fe²⁺、S^2-浓度和流速的情况进行分类.分位点虽然不能表示采样点环境变量的绝对状况，但是能反映采样点环境变量的相对情况；便于同一采样条件下不同样点情况的比较.考虑到采样点的数目不多，以及更清楚地表现DO、Fe²⁺、S^2-浓度和流速在汇流湾区分布规律的要求，采用三分法表示，即利用观测值1/3分位点和2/3分位点将结果分隔为高、中和低3种情况；分位点使用经验分布函数法确定^[34]，计算结果见表 1. 图 2是Fe²⁺、S^2-、DO浓度和流速的分布情况，其多边形用泰森多边形方法生成.

南淝河口湖湾区的Fe²⁺和S^2-浓度分布表现为东高西低；DO浓度的分布则表现为西高东低.而流速在河口很高，当远离南淝河时其入河动量迅速消散，且受地形较高和风吹程不足的影响，流速迅速减小；当离河口较远，风吹程逐渐增大，流速又随之增大(图 2).

总体上，Fe²⁺和DO浓度的分布相反：若DO浓度高，则Fe²⁺浓度较低，反之亦然.这一关系在DO浓度最高的区域表现的尤其明显；而在DO浓度相对较低的区域表现的相对并不明显.一方面是因为DO浓度对Fe²⁺的氧化速率受水体中其它化学因素，如pH值的影响非常大^[11].另一方面则是因为湖体流动对Fe²⁺浓度分布有明显的影响：从图中可以看到，南淝河作为污染的主要来源之一，其DO浓度很低，但河口并非是Fe²⁺浓度最高的区域；这是因为河口流速高，使得南淝河排出含高浓度Fe²⁺的污水团被迅速地推离河口；由于受地形等因素影响，南淝河口东部区域的流速很慢，使得污染水团堆积于东部区域，并进一步造成该区域DO浓度的减少；远离南淝河口的地方在DO浓度不高的情况下，其大部分区域Fe²⁺浓度也较低，这也可以部分归结为该区域相对较高的流速，使得各种污染物分散，从而浓度迅速减小；同时水体的复氧过程则需要一个相对较长的时间才能使DO浓度明显上升，这造成该区域DO浓度并未随着Fe²⁺浓度的减小而立刻上升.这表明，虽然在总体上DO浓度对Fe²⁺浓度的影响要高于流速对Fe²⁺浓度的影响，但在局部区域湖体流速对Fe²⁺的浓度分布有着重要的影响. S^2-浓度的分布和Fe²⁺浓度分布非常类似，这也与S^2-和Fe²⁺的强正相关结果符合.

通过拉格朗日算子法的检验，在空间滞后模型、空间误差模型和空间通用模型这3种空间关系模型中，只有空间滞后模型的显著性系数小于0.05；其余两种模型的显著性系数均大于0.05，因此只需要用空间滞后模型进行回归分析.以Fe²⁺和S^2-为因变量，以DO为自变量时，线性回归模型与空间滞后模型的回归系数都不能通过统计检验，即无论是Fe²⁺或是S^2-浓度都与DO浓度之间不符合这2种模型（表 3）.这一结果反映了DO对Fe²⁺、S^2-的影响复杂，并不是线性关系.

当流速为自变量时，对Fe²⁺和S^2-的线性结果则相反，其线性回归模型和空间滞后模型的回归系数均能通过显著性检验（表 3）.线性回归模型与空间滞后模型相比较，空间滞后模型的显著性更好；且其相应的空间滞后项系数也通过显著性检验.这进一步说明周边水体中Fe²⁺和S^2-浓度分布对本地离子的浓度影响很大.空间滞后模型与线性回归模型的结果比较表明，流速对Fe²⁺和S^2-浓度的影响主要是依靠其对污染物质的分散作用，即水体流动减小Fe²⁺和S^2-浓度；而水体流动对DO复氧产生影响，但影响Fe²⁺和S^2-浓度的作用较小.较之空间滞后模型，线性回归模型低估了水体流动对减小Fe²⁺和S^2-浓度的作用, 低估程度分别为7%和12%.

Moran's I指数侧重的是Fe²⁺、S^2-与流速、DO在空间整体的相关性，而空间回归则反映了三者在局部空间上的关系.空间计量模型结果清楚地表明，在观测期间DO的作用是减小Fe²⁺的总体浓度，而对S^2-浓度总体影响不大.流速主要改变Fe²⁺、S^2-在空间上的局部分布.由于南淝河口湖湾区流速较小，因此该区域水体的流动不足以影响Fe²⁺和S^2-总体浓度.

4 结论

使用空间相关分析和空间回归模型分析巢湖南淝河口区黑臭水体中流速和DO对Fe²⁺、S^2-浓度的分布.结果表明：

1) 该区域流速、DO和Fe²⁺、S^2-浓度的分布都有较强的自空间相关性.互空间相关性研究表明，DO和流速对Fe²⁺浓度在空间上都有较强的互相关性；流速与S^2-也有较强的空间互相关性，但DO与S^2-之间无显著相关性.

2) 增加DO可以在总体上减小Fe²⁺浓度，但对S^2-浓度效果较弱.水体流动使不同地点的DO、Fe²⁺及S^2-产生空间相关性，并改变Fe²⁺和S^2-的局部分布.水体流动通过增强DO复氧，间接影响Fe²⁺浓度的作用较小.

3) 流速与Fe²⁺和S^2-之间均满足线性空间滞后模型.空间滞后模型结果显示，线性回归模型中水流分散Fe²⁺浓度的作用低估7%左右，对S^2-浓度的作用则低估12%左右.