生态系统服务种类的多样性、空间分布的不均衡性及人类使用的选择性使得各种服务之间产生相互作用、相互联系、相互交织的动态关系，具体表现为此消彼长的权衡关系和同减同增（相互增益）的协同关系^［1］。权衡指某些类型生态系统服务的供给，由于其他类型生态系统服务使用的增加而减少的状况^［2］；协同指两种或多种生态系统服务同时增强的情形。近年来，全球和区域生态系统服务之间的权衡协同关系已成为研究的重要课题^［3-4］，主要研究内容包括生态系统服务时空分布特征^［5］、生态系统服务提供和受益地区的空间联系^［6］、多尺度效应的综合分析^［7］及权衡协同关系的驱动因素^［8］。目前评估生态系统服务间权衡协同关系的常用方法有统计方法^［9］、空间分析制图法^［10］、情景模拟法^［11］、生态系统服务流分析方法^［12］等。此外，生态系统服务间的权衡协同关系受到时间变化和空间尺度的影响。多种生态系统服务间的权衡协同关系可随着时间的推移而发生改变^［13］，Hou等^［13］发现2000—2010年黄土高原的产水量与蒸散量之间的权衡关系逐渐减弱，Han等^［14］发现2000—2012年三江源区净初级生产力（NPP）、产水量和水土保持的协同关系呈增强的趋势。除时间变化外，权衡协同关系还具有尺度依赖性，Hou等^［13］的研究结果表明粮食生产与产水量、粮食生产与NPP之间的相关性在不同的空间尺度上表现结果不同；Felipe-Lucia等^［15］在彼德拉河进行的研究表明，在斑块尺度上，食品生产与调节服务间没有显著性相关性，但当扩大到城市和景观水平尺度下，显著正相关；余玉洋等^［16］从秦巴山区全域、综合区域、样点间隔和典型样区等尺度上分析了NPP、粮食生产、土壤保持、水资源供给和生境质量5项生态系统服务间权衡协同关系，结果表明各服务间的权衡协同关系随尺度而发生变化，具有尺度依赖性。上述研究多集中整体尺度或者行政区划尺度分析区域生态系统服务及其权衡协同关系，而且多采用统计关系的数量分析权衡协同关系，缺乏流域内部权衡协同关系时空差异表达。

黄河流域既是中国重要的生态屏障同时是重要的经济地带，黄河流域是连接青藏高原、黄土高原、华北平原的生态廊道，拥有三江源等多个国家公园和国家重点生态功能区。上游是主要的产水区，因受大气环境影响，降水少且蒸散发能力强，导致产水系数低；中游土质疏松、降水强度大且坡度大导致水土流失严重；下游近年来由于城镇化进程加快和经济发展，生物多样性降低。此外，黄河流域布有大量的草地和林地，是重要的碳汇。作为生态环境脆弱的典型代表区域，黄河流域多生态系统服务间的权衡协同关系的时空特征研究，不仅可为区域生态修复和管理提供参考，也是区域内上、中、下游生态补偿和生物多样性保护的科学依据。因此从全域和二级流域尺度上全面探讨权衡协同关系是当前黄河流域生态系统服务与可持续发展研究亟需解决的问题；但是，当前对黄河流域生态系统服务权衡协同关系及其尺度效应研究相对较少。鉴于此，以黄河流域为研究对象，综合土地利用/覆被数据、气象数据、土壤数据和NDVI等多源数据，采用InVEST模型和CASA模型评估黄河流域1990—2018年产水量、土壤保持量、碳储量、生境质量和NPP 5项生态系统服务功能，采用Spearman相关性分析生态系统服务间的权衡协同关系，从不同尺度探讨生态系统服务权衡协同关系，有助于更加科学地管理生态系统，促进流域生态环境可持续发展。

黄河发源于青藏高原巴颜喀拉山北麓的约古宗列盆地，流经青海、四川、甘肃、宁夏、内蒙古、山西、陕西、河南、山东九省区，在山东省垦利县注入渤海，干流全长5 464 km，落差4 480 m。黄河流域位于32°—42°N、96°—119°E，东西长约1 900 km，南北宽约1 100 km。流域面积79.5万km²（包括内流区面积4.2万km²）。河口镇以上为黄河上游，河道长3 472 km，流域面积42.8万km²；河口镇至桃花峪为中游，河道长1 206 km，流域面积34.4万km²；桃花峪以下为下游，河道长786 km，流域面积只有2.3万km²。黄河流域幅员辽阔，山脉众多，东西高差悬殊，各区地貌差异也很大，流域内不同地区气候的差异显著。光照充足，太阳辐射较强，日照条件在全国范围内属于充足的区域，全年日照时数2 000—3 300 h；降水地区差别大，流域大部分地区年降水量为200—650 mm，中上游南部和下游地区多于650 mm，受地形影响较大的南界秦岭山脉北坡降水量可达700—1 000 mm，从西北到东南呈逐渐增加的特点，降水量分布不均，南北降雨量之比大于5。根据全国水资源分区，可将黄河流域分为8个二级流域，分别为龙羊峡以上、龙羊峡至兰州、兰州至河口镇、河口镇至龙门、龙门至三门峡、三门峡至花园口、花园口以下、内流区。

本研究采用多源数据集，包括土地利用/覆被数据集、卫星影像数据集、气象数据集、土壤数据集及相关辅助数据集。土地利用/覆被数据来源于地理国情监测云平台（http：//www.dsac.cn/），空间分辨率为1 km×1 km，土地利用/覆被类型分为耕地、林地、高覆盖度草地、中覆盖度草地、低覆盖度草地、水域、建设用地和未利用地8类；月度降水量数据来源于国家青藏高原数据中心（https：//data.tpdc.ac.cn）中国逐月降水量数据集doi： 10.5281/zenodo.3185722，空间分辨率为1 km×1 km，使用月度数据通过栅格求和得到年度降水数据；土壤沙粒、粉粒、黏粒、有机碳含量、土壤最大根系深度等数据来源于国家冰川冻土沙漠科学数据中心（http：//www.ncdc.ac.cn/）1∶100万土壤数据库；太阳大气顶层辐射数据来源于全球潜在蒸散和全球干旱指数数据集。NDVI数据来源于美国地质调查局（USGS，https：//www.usgs.gov/），采用最大值合成（MVC）方法计算每个像元的最大NDVI值。道路、河流、区县行政边界等基础地理数据等相关辅助数据集来源于国家地理信息中心。

（1） 植被净初级生产力（NPP）。采用CASA模型ENVI遥感估算模块进行评估，该方法将植被分类对评估结果的影响考虑在内，同时由于模型的参数的可获取性和操作的简单性已被广泛应用^［17］：

式中：NPP、APAP、ε分别表示植被净初级生产力（g·m^-2，以C计量）、吸收的光合有效辐射（MJ·m^-2）和实际光利用率（g·MJ^-1，以C计量）；SOL、FPAR分别为太阳总辐射量（MJ·m^-2）和植被对光合有效辐射的吸收比率；常数0.5表示植被利用的太阳辐射率；SR_min取值为 1.08，SR_max与植被类型有关；NDVI表示植被覆盖度；T_ε1、T_ε2和W_ε分别为低温、高温和水分条件对光利用率的胁影响。

（2） 产水量（WY）。采用InVEST模型产水模块评估产水量，该模型计算年产水量将年降水量、年潜在蒸散量以及植被类型、土壤深度和植物有效含水量作为数据输入^［18］，各栅格的降雨量减去实际蒸散量为研究区每个栅格单元x的年产水量Y_x^［19］：

式中：Y_x为栅格x的年产水量，由于年实际蒸散量无法直接测量获取，因此可以采用曲线对AET_x/P_x近似计算；R_x为栅格x的干燥度指数，无量纲；wx是一个经验参数，用来描述气候-土壤属性；AWCx是植被可利用含水量，由于土壤质地和有效土壤深度决定，用来确定土壤为植物生长储存和提供的总水量；Z为Zhang系数^［19］，为常数；ET0x是栅格x内的潜在蒸散量，数据难以获得，通常采用温度法、辐射法和综合法对其进行计算；PAWC表示植物可利用含水率；Kx表示植被或者作物的蒸散系数；MSDx（Max Soil Depthx）为最大土壤深度；RD_x为根系深度；

（3） 土壤保持量（SC）。根据InVEST模型泥沙淤积模块计算土壤保持量，该模型采用土壤流通方程^［19］预测流域每个像元的实际土壤侵蚀：

式中：USLE是年均土壤侵蚀量；R是降雨侵蚀因子；K是土壤侵蚀性因子；LS是坡度坡长因子；C是植被和管理因子；P是支持保护措施因子。该模型不考虑C、P因素计算原始土地覆盖的潜在侵蚀量，并假设景观是裸露的，得出土壤潜在侵蚀量^［20］：RKLS=R×K×LS，基于两种土壤侵蚀量，得出土壤保持量（SC）为：

（4） 碳储量（CS）。采用InVEST模型碳储量模块将生态系统的碳储量划分为4个基本碳库：地上生物碳（所有活植物材料）、地下生物碳（活的植物根系）、土壤碳（分布在有机土壤和矿质土壤中的有机碳）、死亡有机碳（凋落物、死亡树木中的碳）。该模型的Carbon模块以各土地利用/覆被类型为评估单元，以4种碳库的平均密度乘以各评估单元面积来计算研究区域的碳储量。生态系统的碳储量计算公式如下^［21］：

式中：C_tot为总碳储量；C_above、C_below、C_soil、C_dead分别为地上生物碳储量、地下生物碳储量、土壤碳储量、死亡有机碳储量。

基于各地类的碳密度和土地利用/覆被数据，流域内每种土地利用/覆被类型的碳储量计算公式如下：

式中：Ctoti为第i种土地利用/覆被类型的碳储量；Cabovei、Cbelowi、Csoili、Cdeadi分别为第i种土地利用/覆被类型地上、地下、土壤、死亡有机碳储量平均密度；Ai为该土地利用类型的面积。

（5） 生境质量（HQ）。采用InVEST模型中的生境质量模块进行评估，通过建立生境数据与威胁因子之间的联系，在综合考虑威胁因子影响距离和强度等因素，计算威胁因子对生境的负面影响，得到生境的退化程度，然后通过退化程度和生境适宜性计算生境质量，以此来反映区域的生物多样性^［22］。计算公式：

式中：Qxj为生境类型j中栅格x的生境质量指数；H_j为生境类型j的生境适宜度，取值范围在0—1；Dxj为生境退化指数；k为半饱和常数，一般为生境退化度最大值的1/2；z为归一化常量。

相关性分析法可有效反映两变量之间变化趋势的方向和程度^［23］，常用的相关系数有Pearson 系数、 Spearman 系数和 Kendall 系数，其中Speaman 秩相关系数主要利用变量的秩次大小进行线性相关分析，属于非参数统计方法，对原始数据的要求较低^［24］，对输入变量的总体分布形态或样本容量大小不做要求，只要两变量的观测值相互成对即可，甚至可对由连续变量转化而来的等级数据进行分析，适用范围很广^［25］。

式中：PX,Y为变量X和Y的Pearson相关系数；cov(X,Y)为两个变量的协方差；σX和σY为两变量的标准差；PrgX,rgY为应用于原始变量秩次的Spearman相关系数。

首先利用ArcGIS软件创建随机点，然后利用“值提取至点”工具获取各点的5项生态系统服务值；其次对获取样本数据进行z-score标准化处理，消除量级影响；再用KolmogorovSmirnov检测是否存在非正态分布；最后基于R语言的相关性分析函数包探究NPP、产水量、碳储量、土壤保持量和生境质量5种生态系统服务之间的权衡与协同关系。

1990—2018年，黄河流域产水量总体波动增加，2000年最低，2015年最高，总体呈先减少后增加的趋势；碳储量共增加了0.18亿t，总体表现出先减少再增加的变化特征；土壤保持量呈持续减少的趋势，共减少了1.560亿t，减幅为16.94%；平均生境质量指数为0.683，生境质量整体较好，1990—2018年全流域生境质量指数整体上升，增量为0.002；NPP呈先减少后增加的趋势，增幅为17.02%（表1）。

黄河流域5项生态系统服务表现出空间异质性（图1）。产水深度整体呈现出比较一致的规律性，西北地区产水深度低，西南及东南区域产水深度高，龙羊峡以上等黄河上游地区流域产水量最高，兰州至河口镇流域为产水量的低值区，因为该区域蒸发旺盛，且降水偏少。碳储量整体呈东南高、西北低的空间分布特征，具体表现为洛河流域、汾河流域以及渭河流域是碳储量的高值区域。全流域土壤保持量较高的地区主要分布在黄河流域上游以及黄土高原中部地区，而兰州至河口镇以及河套平原地区、黄河下游地区土壤保持量相对较低。生境质量高值区主要分布在上游兰州以上，低值区域主要分布在关中盆地、渭河谷地以及黄河下游平原区。NPP高值区域主要分布在黄河流域下游，包括渭河-关中盆地区域、汾河流域以及花园口以下区域，而低值区域主要分布在黄河流域上中游地区，特别是兰州至河口镇流域等，整体呈现出东高西低的地带性规律。

1990—2018年，全域各生态系统服务间的相关系数均小于0.5（图2），表明各项生态系统服务间相关程度不高，其中NPP与碳储量的相关系数最高，生境质量与碳储量的相关系数次之。产水量与土壤保持量呈协同关系，产水量与碳储量呈权衡关系；NPP与土壤保持量、碳储量、生境质量为协同关系；土壤保持量与碳储量、生境质量为协同关系；碳储量与生境质量为协同关系。产水量与NPP在1990、2000、2010年为协同关系，1995、2005、2018年为权衡关系；产水量与生境质量在1990、1995、2000、2010年为权衡关系，2005、2018年转为弱协同。其他服务间权衡协同的关系并未发生变化，但各权衡协同程度发生变化，产水量与土壤保持量、碳储量与土壤保持量的相关系数基本保持不变，其他各服务间的相关系数波动变化，产水量与生境质量由负相关逐渐变为正相关，相关系数逐年增加，其他各组间的相关系数呈减少的趋势，生境质量与碳储量间的相关系数减少尤为明显。

8个二级流域产水量与碳储量、产水量与生境质量、碳储量与土壤保持量、碳储量与生境质量、碳储量与NPP、土壤保持量与生境质量之间的关系与黄河流域全域均保持一致（表2）。产水量与土壤保持量、产水量与NPP、土壤保持量与NPP以及生境质量与NPP的关系在流域间出现了差异。三门峡至花园口、河口镇至龙门、龙羊峡以上、龙羊峡至兰州以及内流区产水量与土壤保持量与全域关系不同，为权衡关系，而在花园口以下、兰州至河口镇、龙门至三门峡表现为相互增益的协同关系。花园口以下，生境质量与NPP、土壤保持量与NPP的关系与其他流域以及全域不同，为权衡关系。

龙羊峡以上流域产水量与土壤保持量、碳储量和生境质量呈权衡关系，其他各生态系统服务间均为协同关系，其中碳储量与生境质量的相关性尤为显著，平均相关系数为0.707，说明该流域碳储量与生境质量的相互增益关系相较于全域更为显著（图3），高于碳储量与生境质量间的相关系数0.413，NPP与碳储量相关系数较高为0.239，高于其他服务间相关系数，但略低于黄河流域两者间的平均相关系数。龙羊峡以上流域产水量与土壤保持量关系为权衡关系，这与全域不同。龙羊峡至兰州产水量与NPP、产水量与土壤保持量间的关系与全流域不同，由协同关系转为权衡，其他各关系均未发生变化。兰州至河口镇流域5项生态系统服务间的关系与全域完全一致，仅产水量与生境质量及碳储量呈权衡关系，其余生态服务功能之间均为协同关系，但协同程度低于兰州以上地区。河口镇至龙门流域各生态系统服务间的权衡协同关系与全域一致，碳储量与NPP、生境质量与NPP以及生境质量与碳储量均表现出较强的相关性，同时产水量与NPP、碳储量以及生境质量表现出较强权衡关系。相比全流域，该区域土地利用类型多为林地、植被覆盖程度相对较高，NPP、碳储量及生境质量较高，但该区域潜在蒸散发大而产水量相对较小，表现出了较强的权衡关系。内流区产水量与土壤保持量、产水量与NPP与全域有所差异外，其他生态服务权衡协同关系与全流域一致。龙门至三门峡流域产水量与NPP表现出权衡关系，与全域相比，该流域产水量较低，主要林地多分布于此，支持林地的生长需消耗大量的水资源，因此产水量低，同时林地较其他地类的碳密度值较高，具较高碳储量能力，因此二者的权衡关系较强。产水量与碳储量的相关系数变大，说明二者间的权衡关系显著化，与全域的权衡协同关系相比，产水量与碳储量间的权衡关系变强。三门峡至花园口流域产水量与土壤保持量、产水量与NPP的权衡协同关系与全域相反，表现出弱权衡特征，该区域产水量不高，但土壤保持量能力较强，NPP较高，表现出权衡关系。

各生态系统服务整体表现出了明显的流域差异且显示地域规律（图3）。较明显的是在黄河中游流域，生境质量与NPP、碳储量与NPP及土壤保持量与NPP表现出较高的相关性，而上游和下游地区其相关系数较低，说明中游地区协同程度高于上、下游地区，尤其河口镇至龙门以及龙门至三门峡流域，是黄河流域林地分布主要区域，且有大面积草地，生境质量与NPP总体较高。同时，该区域是黄土高原核心区域，由于土壤质地及地形原因，且植被覆盖度高，实际侵蚀量少，土壤保持功能强。碳储量与生境质量相关协同相关系数较高，尤其上游区域以及三门峡以下，协同系数均高于0.5，明显高于兰州至河口镇、河口至龙门镇以及龙门至三门峡流域，因为上游是高协同区，生境质量与碳储量服务功能均高，而下游地区是明显的低值区域，生境质量与碳储量服务功能均弱，两种情况均表现为强协同。

上游流域产水量与生境质量、产水量与碳储量之间的权衡系数小于下游地区，说明下游地区产水量与生境质量、产水量与碳储量的权衡程度更高。上游地区降水丰富，且蒸散发小，以草地为主要地类，产水能力较大，产水量相比下游地区较高，上游区域生境质量相对全域最好，下游地区生境质量较低，同时下游地区因为林地面积分布较集中且平均碳密度大，总体碳储量服务功能较强，因此上游地区的权衡程度小于下游地区，从全域整体看产水量与生境质量为此消彼长的权衡关系。

生态系统功能源于景观格局（斑块结构、廊道布局和景观多样性等）与生态过程（区域水土过程和全球气候变化）间的相互作用及作用的尺度约束，即生态系统服务功能具有空间异质性与尺度效应。这就决定了生态系统服务功能间的相互作用关系也具有了尺度效应，同时系统服务功能间直接相互作用或共同驱动因素的时空差异也会导致服务间的权衡协同关系发生变化^［25］。黄河流域范围广，东西、南北跨度范围大，是中国东西过渡地带，过渡性对区域的阻隔和分异作用产生影响，造成了地形及气候等多方面的差异性，这会对生态系统服务产生影响，促使生态系统服务关系的差异性。

黄河流域生态服务功能权衡协同关系存在尺度差异性，随着尺度的缩小，在二级流域尺度范围发现生态系统服务之间权衡协同关系的差异性，主要表现在产水量与土壤保持量、产水量与NPP、土壤保持量与NPP、生境质量与NPP在二级流域的权衡协同关系与全域不同，产水量与土壤保持量在三门峡至花园口、河口镇至龙门等5个二级流域表现为权衡关系，主要是因为龙羊峡以上、龙羊峡至兰州这两个二级流域位于黄河流域上游，降水量丰沛且蒸散发低，产水量较高，同时该区域坡度大，土地利用类型多为高寒草甸，草层低，截留降水能力差导致水土保持能力较差。而在河口镇至龙门、内流区和三门峡至花园口这些区域，降水不足且潜在蒸散发较高，而且分布有林地和高覆盖草地，坡度相对平缓，导致产水量低而土壤保持量较高。 产水量与NPP间的关系和产水量与NPP间关系相似，在上游龙羊峡至兰州，草层低且降水充足；三门峡至花园口、河口镇至龙门和内流区3个二级流域位于黄河流域内陆，是典型的干旱半干旱气候，降水量少且蒸散发大，是黄河流域 产水量低值区，与此同时，由于2000年后实施“退耕还林还草”政策致使绿色植被大面积增加，从而使NPP增加，而且NPP的增加可以通过增强林冠截留来降低产水量^［26］。在花园口以下，光热条件充沛，适宜绿色植物的生长，而且该区是人口密集和工业加工的重点区域，用水量非常大导致产水量低。土壤保持量与NPP、生境质量与NPP在二级流域的权衡协同关系一致，两者在花园口以下均表现为权衡关系，这主要与花园口是建设用地和经济发展之地密切相关，威胁生境的因子较多，地面多为硬化或者封闭的建设用地，因此土壤保持能力相对较弱。而在其余7个二级流域的关系与全域一致为协同关系。内流区深处黄河流域内陆，是典型的干旱半干旱气候，降水量少且蒸散发大，是黄河流域产水量低值区，而该区域地势相对较为平坦，不宜造成水土流失，土壤保持量较大，产水量与土壤保持量表现出了弱权衡。因此，以流域作为生态系统服务管理单元，能够降低生态调节措施风险，兼顾土壤保持服务和水源供给服务之间的权衡协同关系。但是也存在一定程度的不足，如本研究评估生态系统服务种类不足，无法全面系统地认识黄河流域的生态系统服务，与此同时，模型采用的部分参数是参考了前人相关研究相似区域的文献数据，造成误差，因此未来应结合观测数据与遥感数据，评估更多种类的生态系统服务，定量探究生态系统服务间的相互关系受自然环境、社会经济等因素的影响程度及差异，为生态系统服务管理筛选出适宜的环境条件，促进区域生态系统多种服务可持续发展。

1990—2018年，黄河流域产水量波动增加，碳储量整体减少，土壤保持量波动减少，生境质量和NPP均上升。各项生态系统服务在空间上表现较为显著的差异性，5项生态系统服务的高值区与低值区有较明显的差异。

黄河流域全域各生态系统服务间的相关系数均低于0.5，而且6期相关系数变化不大。产水量与土壤保持量呈协同关系，与碳储量呈权衡关系；NPP与土壤保持量、碳储量、生境质量为协同关系；土壤保持量与碳储量、生境质量为协同关系；碳储量与生境质量为协同关系。

8个二级流域产水量与碳储量、产水量与生境质量、碳储量与土壤保持量、碳储量与生境质量、碳储量与NPP、土壤保持量与生境质量之间的关系与黄河流域全域均保持一致。产水量与土壤保持量、产水量与NPP、土壤保持量与NPP以及生境质量与NPP的关系在流域间出现了差异。各项生态系统服务功能整体表现出了明显二级流域差异性。