黄河流域城镇化与生态系统服务的时空互动关系

邵雅静， 杨 悦， 员学锋

(1.长安大学土地工程学院，西安 710054；2.陕西省土地整治重点实验室，西安 710054)

城镇化是社会经济发展、土地利用模式转变，生态环境变化和生活方式转型的多维发展过程。改革开放以来，中国的城镇化发展取得了显著成果，城镇化率从1978年的17.9%增长到了2019年的60.6%，城镇化率整体提高了42.7%(《中国城市建设统计年鉴》)。在城镇化过程中，生态破坏、环境污染、资源浪费等问题逐渐凸显。人类活动逐渐加速了生态系统服务的退化，进一步扩大了生态系统服务的失衡，威胁人类福祉和全球的可持续发展。城镇化是社会经济发展不可避免的过程，如何平衡城镇化发展和生态系统保护是区域可持续发展的关键问题。

生态系统服务的概念萌芽于20世纪60年代。生态系统服务是指人类可以直接或者间接的从自然生态系统获得的，可以用于生产、生活的福祉和惠益。城镇化是影响生态系统结构、功能以及服务之间权衡协同关系的重要驱动因素。社会经济发展提高了人们对生态系统服务的需求水平。人口集聚导致了自然资源的再分配，同时也增加了粮食、水资源以及能源的消耗。城市建设用地扩张对耕地、林地和草地等生态用地的侵占，加剧了生态系统服务供给类型和供给数量的减少。理清城镇化和生态系统服务之间的作用关系成为生态系统服务领域的研究热点。近年来，研究多关注于快速城镇化地区生态系统服务的时空演变特征，且侧重于量化生态系统服务的经济价值，探析城镇化过程对生态系统服务价值的影响，少有研究关注城镇化指标如人口城镇化、土地城镇化和经济城镇化对生态系统服务供给的影响。目前针对黄河流域的研究较少，已有研究主要包括黄河流域生态系统服务之间的权衡关系量化，生态系统服务的供给和需求测算，土地利用驱动下的生态系统服务价值变化，以及社会经济发展对生态系统服务的影响，也有研究者开始关注黄河流域城镇化与生态系统服务之间的耦合协调关系。然而，少有研究揭示城镇化系统与生态系统服务供给之间的空间互动关系。因此，亟需从时空角度出发，深入探讨生态系统服务与城镇化之间的空间作用关系，为区域的生态管理工作提供一定支撑。

黄河流域是中国重要的生态屏障区和经济发展区。2018年黄河流域9个省区总人口4.2亿人，占中国总人口的30.3%，地区生产总值23.9亿元，占全国的26.5%。黄河流域生态环境脆弱，自然资源短板明显，区域经济发展不平衡问题凸出。因此在合理制定区域生态保护策略的需求下，如何实现城镇化和生态系统服务的协调发展是亟待解决的问题。基于此，本文分析了黄河流域2000-2018年人口城镇化、经济城镇化和土地城镇化过程，并且借助InVEST模型，RUSLE模型，CASA等生物物理模型评估了粮食生产、产水、土壤保持和固碳服务的时空变化，采用Pearson相关系数法和双变量空间自相关方法分析了城镇化与生态系统服务在县域层面上的相互作用关系，以期为黄河流域生态保护与可持续发展提供借鉴和参考。

1 研究区概况

黄河流域位于95°50′29″-119°6′26″E，32°6′53″-41°48′18″N，流域总面积达79.5万km。流经四川、甘肃、宁夏、内蒙、青海、陕西、山西、河南和山东9个省，总面积359.76万km。黄河流域横跨了中国的东中西三大区域，是我国重要的能源基地和重工业基地。本研究将黄河干、支流所流经的县域作为本文的黄河(沿黄)流域范围，共涉及9个省级行政单元，86个市级行政单元和400个县级行政单元(图1)。

图1 研究区地理位置和高程

2 数据来源及研究方法

2.1 数据来源

(1)土地利用数据来源于地理空间数据云平台(https://www.gscloud.cn/)，分辨率为30 m×30 m，采用中科院土地资源分类，主要分为6大类，耕地、林地、草地、水域、建设用地和未利用地。

(2)气象数据来源于国家气象科学数据中心(http://data.cma.cn/)，包括降雨、温度、日照时间和平均风速等；

(3)高程数据(DEM)和NDVI数据来源于地理空间数据云平台，空间分辨率分别为90，250 m；

(4)矢量数据(全国行政区划图，黄河流域矢量范围)来源于中国科学院资源环境科学与数据中心(http://www.resdc.cn/)和国家地球系统科学数据中心(http://www.geodata.cn/)；

(5)土壤数据(土壤类型、土壤质地、土壤有机质含量和根系深度等)来源于兰州寒旱所提供的世界土壤数据库(Harmonized World Soil Database v1.2)，比例尺为1∶100万；

(6)1 km网格人口密度数据来源于中国科学院资源环境科学与数据中心( http://www.resdc.cn/)；

(7)统计数据来源于青海、四川、甘肃、宁夏、内蒙、陕西、山西、山东、河南9省区2000-2018年的统计年鉴；

对于不同来源的数据，在ArcGIS中进行重采样和投影转换，最终都将投影坐标统一为Krasovsky_1940_Albers，分辨率统一为1 km。

2.2 研究方法

2.2.1 服务的量化

(1)粮食生产。参考李双成等的研究成果，根据不同土地利用类型的NDVI值将粮食、油料、水果和蔬菜按照农用地面积，肉类和奶类按照草地面积，水产品按照水域面积进行空间化处理，最终获得全部土地利用类型上的粮食供给空间分布特征，具体计算过程为：

式中：为第个地市第个栅格的粮食产量；为第个地市的粮食产量(kg/hm)；NDVI为第个地市第个栅格的NDVI；NDVI为第个地市的NDVI。

(2)产水服务。本研究使用InVEST模型中的产水模型(water yield)来量化黄河流域的产水量，该模型基于Budyko水热耦合平衡假设，将区域的产水量近似等于年均降水量和实际蒸散发的差值，计算公式为：

式中：()为栅格单元的年产水量(mm)；AET()为栅格单元的年实际蒸发量(mm)；()为栅格单元的年降水量(mm)。

(3)固碳服务。植被净初级生产力(net primary productivity，NPP)可以直接反映植被群落在自然环境条件下的生产力，表征固碳服务。CASA(carnegie ames-stanford approach)模型被广泛用于估算全球生态系统服务的净初级生产力，计算公式为：

NPP(,)=APAR(,)×(,)

式中：NPP(，)为像元在月的净初级生产力(g C/m)；APAR(，)为像元在月吸收的光合有效辐射(MJ/m),(,)为像元在月的实际光能利用率(g C/MJ)。

(4)土壤保持。土壤保持量(soil conservation)是潜在土壤侵蚀量和实际土壤侵蚀量之间的差值，修改的土壤流通方程RUSLE(revised universal soil loss equation)被广泛应用到土壤侵蚀量的计算过程中，计算公式为：

=RKLS-USLE=×××(1-×)

式中：为土壤保持量(t/(hm·a))；RKLS为在没有植被覆盖和土壤保持措施下土壤的潜在侵蚀量；USLE为区域土壤的实际侵蚀量；为降雨侵蚀性因子((MJ·mm)/(hm·h·a))，通过Wischmeier等提出的公式计算；为土壤可蚀性因子((t·hm·h)/(MJ·hm·mm))，利用Williams等提出的EPIC公式计算；和分别为坡长和坡度因子，坡长采用Wischmeier等提出的公式计算，坡度因子采用Mc Cool等提出的公式计算；为植被覆盖因子，参照蔡崇法等的研究成果计算；为土壤保持措施因子，参考Lufafa等方法计算。

2.2.2 城镇化指标 城镇化主要表现为区域经济增长、人口集聚和城市用地的扩张。根据黄河流域的城镇化现状和数据的可获性，本文主要从人口、土地和经济3个城镇化维度来分析黄河流域的城镇化过程，采用人口密度(人/km)、GDP密度(万元/km)和建设用地占比(%)来分别表示人口城镇化、经济城镇化和土地城镇化指标，用3个指标标准化处理之后的叠加结果表示综合城镇化指数(comprehensive urbanization index，CUI)。

2.2.3 两者关系量化

(1)Pearson相关系数。本研究运用相关系数法分析城镇化与生态系统服务之间关系的时序演变特征。若相关系数为正并且通过了显著性检验则为正相关，相关系数为负并且通过了显著性检验则为负相关。具体的计算过程为：

(2)双变量空间自相关分析。本研究采用双变量空间自相关(bivariate spatial autocorrelation analysis)分析城镇化与生态系统服务在空间上的相关特征及分布情况。双变量全局自相关和局部自相关方法能够反映2个变量在整个研究区的相关关系以及空间关系的局部集聚特征。Moran’s取值[-1，1]，Moran’s>0表示正相关，Moran’s<0表示负相关，Moran’s=0则表示属性值是随机分布的。借助GeoDa软件，在值检验(=0.05)的基础上绘制LISA聚类图(local Indicator of spatial association)来反映城镇化和生态系统服务在空间上的依赖性和关联性。输出结果包括：高值聚类(HH)、低值聚类(LL)、高低集聚(HL)、低高集聚(LH)和不相关(NN)。

3 结果与分析

3.1 2000-2018年黄河流域的城镇化发展过程

由图2可知，2000-2018年黄河流域综合城镇化水平普遍有所提升，综合城镇化率提高了102.60%，其中人口城镇化(人口密度)均值从2000年的269.56人/km增加到了2018年的322.70人/km，土地城镇化(建设用地占比)均值从2000年的5.12%增长到了2018年的8.06%，经济城镇化(GDP密度)均值从2000年的228.56万元/km增长到了2018年的1 792.07万元/km，变化幅度分别为19.71%，57.36%，684.06%。从变化速率来看，整体表现为经济城镇化速度>土地城镇化速度>人口城镇化速度。

图2 2000-2018年黄河流域城镇化空间分布格局及变化

从城镇化指标的空间分布特征来看，人口城镇化、土地城镇化和经济城镇化的区域不平衡特征显著。2000-2018年黄河流域78.25%的地区人口密度有所增加，同时21.75%的地区人口密度有所减少。2000-2018年黄河上游地区人口密度基本处于200人/km以下，人口密度>800人/km的地区基本分布在下游的山东省和河南省以及中游的陕西省和山西省。2000-2018年黄河流域94.25%地区建设用地的占比呈增长趋势，而5.75%的地区建设用地的占比有所减少，黄河中下游地区土地城镇化的集聚特征最为显著。2000-2018年黄河流域建设用地占比超过20%的地区面积占比从2000年的5%左右增长到了2018年24%。2000-2018年黄河流域GDP密度整体增加，经济城镇化水平普遍提升且增长的区域主要分布在省会城市及其周边地区。整体来看，黄河流域的城镇化过程具有显著的空间差异性，黄河中下游地区的经济发展优势和人口集聚程度显著高于上游地区，是黄河流域人口城镇化、经济城镇化和土地城镇化的核心地区。

3.2 2000-2018年黄河流域生态系统服务的空间分布特征

2000-2018年黄河流域的粮食生产、固碳和土壤保持服务均值分别增长65.32%，59.47%，46.27%，而产水服务均值减少21.15%(图3)。整体来看，黄河流域70%以上的地区粮食生产、固碳和土壤保持呈现增长趋势。2000-2018年黄河流域粮食生产的增长区集中分布在下游的山东省和河南省，这些地区位于华北平原，水热资源丰富，平均粮食单量高于4 000 kg/hm。2000-2018年固碳服务在黄河中下游地区的增长最为显著，2018年固碳服务的高值区主要集中在黄河中游陕西省南部的秦巴山区，固碳服务均值基本大于500 g C/m，这些地区植被覆盖度较高，是我国可持续发展的生物基因库，也是我国重要的生态屏障之一，固碳服务供给优势明显。2000-2018年土壤保持服务在整个黄河流域呈增长趋势，位于黄河上游的青海省东南地区增长最为显著，这些地区位于青藏高原，植被覆盖度较高，区域的土壤保持能力较强。2000-2018年黄河流域的产水服务均值由2000年93.40 mm减少到2018年的73.65 mm，产水服务在黄河上游青海省地区增长显著，而在中下游地区普遍减少，减少的区域占黄河流域总面积的55.50%，主要分布在山西省和河南省的交界地区。

图3 2000-2018年黄河流域生态系统服务的空间分布特征及变化

从综合服务的变化特征来看，2000-2018年黄河流域综合生态系统服务普遍增长，综合生态系统服务均值增长21.09%，增长的区域占比98.50%。

从分布特征来看，位于黄河流域中游的甘肃省庆阳市和平凉市增长最为显著，增长量均高于0.3，其次为中游地区陕西省的延安市和宝鸡市，增长量在0.2～0.3(图4)，黄河中下游地区的综合生态系统服务显著高于上游地区。

图4 2000-2018年黄河流域综合生态系统服务时空分布及变化趋势

3.3 黄河流域城镇化和生态系统服务的时空关系分析

(1)2000-2018年黄河流域城镇化和生态系统服务的相关关系。本文利用Pearson相关系数量化了黄河流域城镇化和生态系统服务之间关系的时序特征(图5)。2000-2018年黄河流域城镇化指标和生态系统服务之间均具有显著的相关性，其中土地城镇化、人口城镇化和经济城镇化对固碳和土壤保持服务有消极影响(<0)，而对产水服务和粮食生产有一定的积极影响(>0)。随着时间的变化，城镇化对生态系统服务的影响程度在不断的变化。土地城镇化、人口城镇化和经济城镇化对固碳服务和土壤保持的影响在不断的加强，而对粮食生产和产水服务的影响有所减弱，也就是说在时间尺度上城镇化发展对生态系统服务的积极影响正在减弱而消极影响正在增强。

注：*表示p<0.05；**表示p<0.01；***表示p<0.001。

(2)2000-2018年黄河流域城镇化和生态系统服务的空间关系。为了进一步揭示城镇化与生态系统服务在空间上的作用关系，本文借助GeoDa软件来分析黄河流域城镇化与生态系统服务之间的空间互动关系。在1%的显著性水平下，2000年、2010年和2018年城镇化和生态系统服务之间的Moran’指数分别为0.140 0，0.090 4，-0.045 7(图6)。从全局自相关结果来看，黄河流域城镇化与生态系统服务之间的空间关系从正相关关系逐渐转变成为负相关关系。

在LISA聚类结果(图6)中，HH代表城镇化和生态系统服务高值-高值相聚集的状态，即高城镇化水平对应了较高的生态系统服务水平；HL表示城镇化和生态系统服务的高值-低值聚集状态，即城镇化水平高而生态系统服务水平低；LH和HL则刚好相反；LL则代表城镇化和生态系统服务的低-低聚集状态，即区域的城镇化水平低且生态系统服务水平也较低。黄河流域城镇化与生态系统服务HH区分布较少且零星分布在黄河中下游地区，2000-2018年HH区分布面积减少5.5%，LL区分布面积减少1.75%。LL分布在黄河流域上游以及中游的北部地区，黄河上游是典型的城镇化发展滞后区，尤其是上游地区的青海、内蒙和甘肃等地，脆弱的生态环境和匮乏的水资源很大程度上限制了区域的土地资源开发和规模利用，制约了城镇化发展。2000-2018年HL区和LH区分别增加0.25%和7.75%，其中LH区主要分布在黄河流域南部，HL零星分布在黄河上游地区的内蒙地区。整体来看，黄河流域城镇化和生态系统服务之间发展不协调的地区逐渐增多。因此，基于城镇化和生态系统服务之间的空间关系因地制宜地制定区域生态保护措施显得尤为重要。

图6 黄河流域城镇化与生态系统服务之间的全局Moran’I指数(P=1%)和局部聚类LISA图(P=5%)

4 讨 论

理清黄河流域城镇化与生态系统服务的时空关系有助于实现区域社会经济发展和生态保护的“双赢”目标。城镇化发展会改变生态系统结构，从而影响生态系统服务的供给。黄河流域的城镇化发展迅速，经济城镇化和土地城镇化速度均高于人口城镇化速度，符合城镇化的发展趋势。黄河流域的生态系统服务具有显著的空间异质性，在黄河上游、中游和下游地区差异较大，这也与Zhang等的研究成果一致。2000-2018年黄河流域的综合生态系统服务有所增长，尤其在黄河中上游地区最为显著，这主要益于我国从1978年开始实施的“三北防护林工程”“退耕还林还草工程”等一系列的生态恢复工程，使得黄河中上游地区的生态环境得到了显著改善。

城镇化过程会改变生态系统结构，从而影响生态系统服务供给。城镇化要素的集聚特征也进一步加剧了生态系统服务的空间再分配。2000-2018年黄河流域城镇化和生态系统服务的空间作用关系逐渐从正相关转变为负相关，这与Tian等的研究成果相似。主要是因为在城镇化发展的初期阶段，政府制定了一系列的生态恢复政策来实现区域的生态保护，这一阶段城镇化和生态系统服务均处于稳步增长阶段，2010年以后，随着城镇化进程的加快，生态环境的压力骤增，加剧了城镇化对生态系统服务的负面效应。同时城镇化以资源的高消耗、生态系统服务的流失为代价的发展模式，以及资源的不合理开发和过度利用进一步加速了生态系统服务的衰减。尤其是大量的耕地、林地和草地被城市开发建设所占用，导致生态用地面积不断减少，生态空间受到挤压，最终影响了区域的生态服务供给。2000-2018年黄河流域城镇化与生态系统服务呈负相关且发展不协调的空间区域(HL区和LH区)有所增长，而城镇化和生态系统服务HH区和LL区逐渐减少。城镇化和生态系统服务HH区的分布面积较少且零星分布在黄河中下游地区，这些地区的城镇化水平和生态系统服务供给水平较高。这也印证了Pan等的研究，在城镇化发展水平较高的地区生态系统服务的可达性更高，而且城市发展的可能性更大。未来可以根据城镇化和生态系统服务关系的区域特征进行分类施策以提高区域生态保护效率。可以将城镇化和生态系统服务的HH区作为黄河流域高质量发展的重点生态保护区，LL区可以作为生态限制区，HL区和LH区分别作为生态恢复区和生态控制区，同时重视对HL区和LH区的监管。此外，根据土地城镇化、人口城镇化和经济城镇化对黄河流域不同地区生态系统服务的影响机制，进一步结合流域的实际情况，采取有效措施控制城镇化发展对区域生态系统服务的过渡干扰，加强城市绿色基础设施的建设，通过针对性的策略来实现黄河流域城镇化和生态系统服务的可持续发展。

5 结 论

(1)2000-2018年黄河流域的综合城镇化水平整体提升，其中经济城镇化实现了爆发性的增长，其次是土地城镇化和人口城镇化发展，变化率分别为684.06%,57.36%，19.71%。

(2)2000-2018年黄河流域的粮食生产、固碳和土壤保持服务均呈增长趋势，而产水服务逐渐减少，变化率分别为65.32%,59.47%，46.27%，-21.15%。整体来看，黄河上游地区土壤保持和产水服务明显增加，黄河中下游地区粮食生产和固碳服务显著增长而产水服务明显降低。

(3)在时间尺度上，黄河流域人口、土地和经济城镇化对固碳服务和土壤保持服务的负向影响在不断加强，而对粮食生产和产水服务的正向影响不断减弱。在空间尺度上，城镇化对生态系统服务的负向影响逐渐增强，并且城镇化发展和生态系统服务不协调区域不断增加。