随着社会经济的发展，水资源已成为全球多地区发展的限制因子之一^[1]。作为生态系统的重要服务功能，区域水资源的产生与供给对确保灌溉农业、工业及旅游业的发展，促进居民生活水平的提高具有重要意义^[2]；因此，流域生态系统产水服务及其影响机制的研究逐渐成为水文学、流域规划与管理等领域的研究热点^[3]。已有研究表明：在长时间尺度上，气候是降水、产水等水资源变化的主导因素；但在较短时间尺度上，土地利用/覆被变化(land use and land cover change，LUCC)才是水文效应的重要影响因素^[4]。揭示和预测LUCC的水文效应，对指导区域土地利用规划，优化区域水土资源配置，以及流域的可持续发展至关重要^[5]。近年来，关于土地利用与水文过程关系的研究正逐渐增加，特征参数^[6]、卫星遥感^[7]、水文模型^[8]等方法技术被综合应用，植被变化、农田增减、城镇化等用地变化被认为可能影响水文过程^[9-11]。目前，大部分研究都关注于土地利用比例对水文过程的影响，而对景观格局的水文效应研究较少^[12]；而且，在探讨土地利用水文效应时，较多研究采用相关分析和多元线性回归分析的方法，没有考虑到空间因子，忽略了空间依赖性和空间非平稳性对水文效应的影响。

自20世纪60年代起，随着计算机技术的发展，水文效应的研究方法由传统的统计分析向水文模型转变。InVEST模型便是目前应用广泛的一种。InVEST(the integrated valuation of ecosystem services and tradeoffs tool)模型，全称“生态系统服务功能综合估价和权衡得失评估模型”，由美国斯坦福大学、大自然保护协会(TNC)和世界自然基金会(WWF)联合开发而成^[8]。使用该模型能够快捷地对包括产水功能在内的生态系统多种服务功能进行量化，从而很好地为政府及相关部门提供决策依据。

基于此，本文选择水资源供需矛盾突出的黑河中游为研究区，建立InVEST模型计算流域产水：一方面，将空间因子考虑在内，通过普通最小二乘回归、空间回归以及地理加权回归等多种方法，分析土地利用类型与流域产水的关系；另一方面，结合主成分分析筛选景观指数，进一步揭示景观格局与区域产水的关系，从而为区域土地利用规划与决策提供参考，促进水土资源的持续利用。

1 研究区概况

黑河是中国第2大内陆河，发源于祁连山中段，流经青海、甘肃、内蒙古3省自治区^[13]。黑河莺落峡至正义峡之间为中游，东靠武威和金昌，西至嘉峪关和酒泉，南与青海省接壤，北和内蒙古毗邻，具体范围介于E 97°20′~102°12′，N 37°28′~39°57′。中游干流全长185 km，总面积约1.96万km²。黑河中游地处青藏高原向内蒙古高原的过渡地带，海拔在1 234~4 886 m之间，地势东南高西北低，其北部山地和走廊平原属温带大陆性干旱气候，南部祁连山区为高寒半干旱气候。年降水量仅54.9~436.2 mm，且年内分布不均，多集中在6—9月；但全年蒸发量高达1 700 mm^[14]。黑河中游是典型的灌溉农业绿洲，也是整个流域人类活动最为密集、绿洲最为集中、经济最为发达的地区^[15]。流域人均可利用水资源量为1 250 m³(仅为全国平均水平的57%)，是典型的资源型缺水地区。

本研究所需数据包括基础地理信息数据、土地利用数据、气象数据和土壤数据等。基础地理信息数据来自中国科学院寒区旱区科学数据中心的黑河计划数据管理中心(www.heihedata.org)。土地利用数据包括黑河中游2000、2010和2014年土地利用栅格数据。其中，2000年土地利用数据来源于中国科学院“八五”重大应用项目^[16]，2010、2014年土地利用数据是通过黑河中游Landsat影像进行人机交互式目视解译得到，解译精度均＞88%。研究中，将土地利用类型划分为耕地、草地、林地、水域、建筑用地和未利用地6类，数据的空间分辨率为30 m。

气象数据来自中国气象数据网(cdc.nmic.cn/home.do)。对研究区及其周边的11个气象站点求得多年平均降水量，再通过普通克里金插值得到黑河中游多年平均降水量。使用“改进的哈格里夫斯”(modified Hargreaves)公式计算并插值得到多年平均参考蒸散量^[17]。土壤深度、植物可利用含水量数据均来自世界土壤数据库(HWSD)，空间分辨率为1 km。

DEM数据来自美国马里兰大学全球土地覆被数据库，分辨率为30 m。基于DEM数据，在ArcGIS中使用Hydrology工具将研究区划分为20个子流域(图 1)。

2 研究方法 2.1 基于InVEST模型的产水量模拟

InVEST模型产水模块的计算原理，主要是通过降水量减去蒸发蒸腾水量来得到产水量。该模型忽略地表水与地下水的交互作用，不区分地表径流、地下径流以及基流。

产水模块是基于Budyko曲线和年平均降水量进行计算的。年产水量的计算公式为

$Y\left( x \right) = \left( {1 - \frac{{{\rm{AET}}\left( x \right)}}{{P\left( x \right)}}} \right)P\left( x \right)。$

(1)

式中：Y(x)为x类景观像素的年产水量，mm；AET(x)为年实际蒸散发量，mm；P(x)为年降水量，mm。

2.2 回归分析模型

1) 普通最小二乘回归(ordinary least squares，简称OLS)是最简单、最基本且应用最为广泛的模型之一。模型的一般形式为

$Y = \sum\limits_{i = 1}^n {{X_\mathit{i}}{\mathit{\beta }_\mathit{i}}} + \mathit{\varepsilon }\mathit{。} $

(2)

式中：Y为因变量，mm；X_i为第i个自变量，mm；β_i为回归系数；ε为随机误差项，mm。

2) 空间误差模型(spatial error model，简称SEM)通过空间误差项来引入空间依赖性。其表达式为

$Y = X\mathit{\beta + \varepsilon , } $

(3)

$\mathit{\varepsilon = \lambda W\varepsilon + \mu }\mathit{.} $

(4)

式中：λ为自回归参数，衡量样本观测值中的空间依赖作用；μ为空间残差项，mm；W为空间权重矩阵。

3) 地理加权回归(geographical weighted regression，简称GWR)是一种对空间变化关系建模的线性回归的局部形式，是对普通回归模型的扩展。该模型假定了其拟合系数是变量空间位置的函数，将变量的地理位置引入到回归参数中，再利用局部加权最小二乘法进行逐点的参数估计^[18]。其中，权是回归点所在空间位置到其他各采样点空间位置的距离函数。其表达式一般为

${Y_i} = {\mathit{\beta }_0}\left( {{u_i}, {v_\mathit{i}}} \right) + \sum\limits_{k = 1}^p {{\mathit{\beta }_k}\left( {{u_\mathit{i}}, {v_\mathit{i}}} \right)} {X_{ik}} + {\mathit{\varepsilon }_\mathit{i}}。$

(5)

式中：Y_i为第i个采样点的因变量，即待估测值，mm；X_ik为第i个采样点上的第k个解释变量，mm；(u_i, v_i)为第i个采样点的坐标；β_k(u_i, v_i)为第i个采样点上的第k个回归参数，是地理位置的函数。

在回归分析中，分别以土地利用类型和景观指数指标作为自变量。其中，土地利用类型指标以每个子流域的各类用地面积占该子流域总面积的比值来表示，计算公式如下：

${X_{\mathit{ij}}} = \frac{{{S_{ij}}}}{{{S_\mathit{i}}}}。$

(6)

式中：X_ij为第i个子流域中第j类用地面积指标；S_ij为第i个子流域内第j类用地的面积，m²；S_i为第i个子流域的总面积，m²。

3 结果分析 3.1 黑河中游产水量时空格局

2000、2010和2014年的产水量分别为1.73亿、1.63亿和1.38亿m³，呈现出逐年减少趋势，尤其是2010年至2014年，产水量减少幅度更大。

图 2中，虚线是整个黑河中游当年的平均产水深度，对比可知子流域间产水能力差异极大。三期中，产水集中于15~20号这6个子流域，15、16、18和19号子流域的产水深度远高于平均产水深度，而1~14号子流域产水深度则远低于平均产水深度，部分子流域甚至为0。

由图 3可见，三期产水量均呈现出东南多、西北少的空间分布。对比产水量与降水量分布格局，发现尽管产水量和降水量都有着东南多、西北少的大格局特征，但在小格局上的分布规律并不完全一致。说明产水量受到降水量的控制，但子流域间的局部差异还受到下垫面的影响。

3.2 黑河中游土地利用的水文效应 3.2.1 土地利用水文效应的普通最小二乘回归

在土地利用水文效应的回归建模中，6类指标间具有线性关系，若同时将其作为自变量，会造成变量冗余；因此，先去掉草地(16.29%)、未利用地(57.38%)这2类用地(两者在研究区所占面积最大，会对其它类型用地面积指标值影响较大)，只将耕地、林地、水体、建设用地4类用地指标作为自变量，将极差标准化后的产水指标作为因变量，进行普通最小二乘回归(表 1)。

对于P＞0.05的自变量，其系数不具有显著性，尚不能认为有统计学意义，在此不予讨论^[19]。由表 1可知，在2000年的回归模型中，林地和水体的系数为正，表明林地和水体促进产水。在InVEST产水模型中，产水量与产流量的意义有所不同，模型计算的产水量包含在之后可能转化为的土壤含水量、凋落物持水量以及冠层截留量；因此，林地因其有树荫遮蔽，整体温度较低，蒸发速度慢，且其叶片、根系等具有持水功能，整体上表现为促进产水。而水体的蒸发形式是水面蒸发，主要受控于气象条件，由于水体比热容较大，其温度会低于陆面，使蒸发量往往更小，从而促进产水^[20]。2010年，耕地的回归系数为正，表明耕地也有利于产水。建设用地的系数为负，表明建设用地对产水有抑制作用，这可能由于城镇的建筑用地具有热岛效应，致使水分的蒸散发量较大，从而不利于产水。2014年的回归结果与前2期相似，耕地和林地促进产水，建设用地抑制产水。

同样地，将未利用地和草地这2类用地指标作为自变量，将极差标准化后的产水指标作为因变量，进行普通最小二乘回归(表 2)。

由表 2可见，草地的系数均不显著，2000、2010和2014年中未利用地的系数均为负，表明未利用地抑制产水。黑河中游的未利用地主要分布在北部地区，其地面多表现为荒漠，土壤类型主要为黏重盐沙，蓄水能力差，蒸散量大，故不容易产水。

3.2.2 土地利用水文效应的空间回归

空间依赖性检验能对模型选择做出判断。通过查看各普通最小二乘回归模型的空间依赖性检验结果，发现2010、2014年的回归模型中拉格朗日乘数均不显著(P>0.05)，表明对其保留OLS回归结果即可，而以2000年草地和未利用地为自变量的回归模型中，LM-lag不显著(P=0.28)而LM-error显著(P=0.03)，表明有必要使用空间误差回归模型进行改进。

于是，以2000年草地和未利用地指标为自变量，以产水深度为因变量，建立空间误差模型(表 3)。

由表 3可知，空间回归与普通最小二乘回归结果相同，未利用地系数为负，抑制产水。且空间误差模型的R²为0.71，而普通最小二乘回归的R²为0.59，可见模型性能得到的了较大提升。为检验空间回归模型效果，分别计算OLS和SEM的回归残差的全局Moran’s Ⅰ，从而检验回归模型的残差是否具有空间随机性。经过计算得到，OLS残差的Moran’s Ⅰ为-0.43(P=0.03)，表明OLS的残差具空间负相关性。而SEM残差的Moran’s Ⅰ为0.12(P=0.11)，与OLS相比，指数的绝对值明显减小，结合指数对应的P值，可以认为空间误差模型消除了残差的空间负相关性。

3.2.3 土地利用水文效应的地理加权回归

根据地理学第一定律，OLS中的平稳性假设并非恒成立，而地理加权回归能够容许一些不平稳的数据被模拟。经计算，2000年黑河中游产水深度的Moran’s Ⅰ为0.32(P=0.02)，表明因变量确实具有空间正相关性，适合采用地理加权回归模型进行拟合。

以2000年草地和未利用地指标为自变量，以产水深度为因变量，建立地理加权回归模型。得到的R²为0.90，远高于OLS的R²，证明地理加权回归明显地提升模型性能。地理加权回归模型对局部进行了回归，因此每个子流域都有1个回归系数，未利用地在各子流域的回归系数如图 4所示。

由图 4可见：与普通最小二乘回归和空间回归结果相一致，未利用地的地理加权回归系数整体上为负，但在各子流域中存在差异。在最南端，回归系数为正，表明未利用地在此处促进产水；在北部地区，系数为负，表明抑制产水；而中偏南部地区，回归系数达到最小值-1.17，表明未利用地在此处对产水的抑制作用最强。这种分布可能与黑河中游土壤质地及局部地形有关，南部土壤松软，容易蓄水持水，蒸发量相对减少，因此容易产水；由南向北土壤变得黏重，蓄水能力差，再加上中部主要为平原地形，接受太阳辐射多，蒸发量相对更大，因此不利于产水；而北部地区虽然土壤也属于黏重型，但因山地地形的存在，使当地有山地阴影、坡地小气候等对降水、蒸发产生影响，从而使北部未利用地对产水的抑制作用不及中部偏南地区强烈。

结合研究区2000—2014年土地利用变化趋势分析，总体来看，黑河中游产水能力的下降受到了未利用地面积增加以及林地、水体、耕地面积减少的影响。

3.3 黑河中游景观格局的水文效应 3.3.1 景观格局指数主成分分析

参考已有研究^[21-23]，选择目前常用、且具有代表性的景观指数共9个，包括最大斑块占景观面积比例(LPI，largest patch index)、斑块面积算术平均(AREA_MN，mean of patch area)、斑块面积标准差(AREA_SD，standard deviation of patch area)、形状指数加权平均(SHAPE_AM，area-weighted mean of shape index)、分形维数加权平均(FRAC_AM，area-weighted mean of fractal index)、斑块多度(PR，patch richness)、斑块多度密度(PRD，patch richness density)、修正Simpson多样性指数(MSIDI，modified Simpson’s diversity index)、修正Simpson均匀度指数(MSIEI，modified Simpson’s evenness index)。基于2000、2010及2014年土地利用数据，使用Fragstats 4.2^[24]计算出黑河中游的9种景观指数。

对黑河中游3期景观指数进行主成分分析，结果如表 4所示。

由表 4可见，3期数据的前3个特征根的累计贡献率均大于85%，即3个年份的景观格局均可使用前3个主成分描述。对比不同年份的主成分载荷可以发现，第1主成分主要由LPI、FRAC_AM、MSIDI、MSIEI来解释，主要代表景观丰富度、均匀度；第2主成分主要由AREA_SD解释，代表斑块面积特征；第3主成分主要由PRD解释，代表景观斑块相互间影响情况^{[23, 25]}。

3.3.2 景观格局水文效应的普通最小二乘回归

基于主成分分析，将极差标准化后的主成分得分作为自变量，将产水深度作为因变量，进行普通最小二乘回归(表 5)。

由表 5可见，3期数据中第1主成分的系数均显著，且均为正，表明景观丰富度、均匀度与产水成正相关。2000年中的第2主成分系数也显著，其系数为负，第2主成分主要由斑块面积标准差AREA_SD解释，表明斑块间面积差距越大即景观均匀程度越低，越不利于产水。综合推断，在黑河中游，景观类型越丰富多样，景观斑块分布越均匀，则越有利于产水。

总体来看，2000—2014年期间研究区景观多样性和景观均匀度的下降，也是黑河中游产水能力下降的原因之一。

4 结论与讨论

1) 黑河中游产水能力在2000—2014年期间呈现出减小趋势，与未利用地面积增加和林地面积减少、以及景观多样性和均匀度下降等景观格局变化有直接关系。2000、2010和2014年产水量分别为1.73亿、1.63亿和1.38亿m³。产水量空间分布不均匀，东南多、西北少，子流域间产水能力差异大，这与王蓓等^[26]的研究结论基本一致。

2) 黑河中游水文效应的尺度特征明显。在土地利用尺度上，林地、水体、耕地等有助于促进产水，而建设用地和未利用地则抑制产水；在景观尺度上，景观类型越丰富、斑块分布越均匀，越有利于产水，这与前人的研究结论相近^{[20, 27]}。因此，结合黑河中游生态安全建设^[28]，对于北部地区，加强未利用地的封育保护，增加荒漠植被群落的稳定性，防止沙漠化的进一步扩展；对于水土条件较为优越的中部偏南地区，在维持地类多样性的同时，控制农业绿洲无序扩展，协调各类用地比例，提升景观斑块的均匀度，合理配置生态用水；对于产水量较高的南部地区，加快林草植被的保护与恢复，保障流域源头的森林生态屏障。

3) 普通最小二乘回归、空间回归以及地理加权回归结果均能表明黑河中游未利用地整体上抑制产水。但与普通最小二乘回归模型相比，空间回归模型消除了回归残差的空间负相关性，使模型性能得到提升。地理加权回归使模型性能得到进一步提升，并揭示出未利用地回归系数的空间差异。

4) 笔者从土地利用和景观格局2个尺度揭示黑河中游水文效应的影响因素与影响机制，为区域生态系统管理提供参考。本研究关注揭示土地利用活动对水文效应的影响，没有考虑气候变化对产水的影响。今后在更长时间序列的研究中，可结合气候变化来分析对流域产水的综合影响。另外，研究中所选的景观格局指数虽依据主成分分析和已有研究，但仍难避免主观因素的影响。未来研究可结合研究区实际情况改进景观格局指数筛选方法，深入探讨景观格局的水文效应。