随着对全球变化的深入研究, 人们发现土地利用/覆盖变化(Land Use/Cover Change, LUCC)是影响全球环境变化的重要原因^[1-2]。LUCC不仅影响气候变化, 还进一步影响了区域生态系统服务价值(Ecosystem Service Value, ESV)的空间分布和景观格局变化的生态过程^[3-5]。ESV指人们从生态系统中获得的明显或潜在的利益^[6-7], 其主要指标包括供给服务、调节服务、支持服务和文化服务^[8-9]。目前, 经济价值评估的使用主要遵循由Costanza等人^[10]建立的ESV评估模型。根据我国实际情况, 谢高地等^[11-12]研究出了基于该模型的单位面积生态服务价值表, 并得到了广泛应用^[13-14]。景观格局通常指景观空间结构特征, 可以反映潜在的人类活动行为, 已逐渐成为土地利用变化研究的主要分析手段之一^[15]。目前, 景观格局的演变主要依靠景观指数分析和空间统计特征方法^[16-17], 如何清晰地理解和衡量土地利用、景观和生态系统服务价值的模式和过程的变化强度, 以及它们之间错综复杂的关系, 是一个亟待解决的问题。

强度分析是一种定量的数学理论框架, 相较于其他常见的LUCC分析模型, 强度分析可以系统地定量地分析土地利用动态变化过程中的内部转换状况^[18-19], 并检验在时间间隔、地类和转换水平上的强度变化^[20]。Huang^[21]等人第一次将强度分析方法应用在中国的东南部九龙河流域, 探究该地区的LUCC过程和模式。杨建新等^[22]人对武汉市土地利用强度分析时, 使用了一种转换模式交叉列联表展示各地类转换类型的稳定性与系统性特征。但大部分使用强度分析的LUCC分析中, 同时考虑景观格局与区域ESV之间相互影响关系的研究较少, 因此, 利用强度分析的优势弥补其他LUCC分析模型的不足, 进而深入挖掘其与景观格局和ESV的响应关系, 具有重要的现实意义。

为了进一步实现可持续发展的目标, 从景观格局变化的角度理解当今经济发展所造成的生态效益是研究热点之一^[23-24]。但是很少有研究分析土地利用变化的强度如何影响生态系统服务。因此, 本研究在强度分析中引入改进后的土地转换模式交叉列联表, 研究伊犁河流域(中国部分)1980—2015年LUCC、景观格局和ESV的变化情况；同时, 探究LUCC、景观格局与ESV之间的关系。以期为伊犁河流域未来土地资源持续利用和经济的可持续发展提供理论和方法支持。

1 研究区概况和数据来源 1.1 研究区概况

伊犁河流域(中国区域)地处新疆西部的伊犁哈萨克自治州, 位于东经80°09′42″—84°56′50″, 北纬42°14′16″—44°53′30″之间, 流域面积57550.80km²。流域内多种地貌类型并存, 地势东高西低, 东窄西宽, 形似朝西开口的喇叭状, 概括为“两山夹一盆”, 最高海拔5741m, 最低海拔532m。受地貌影响, 流域内降水丰富, 形成了西北干旱半干旱区独特的“湿岛”气候。年平均气温在2.9—10.4℃左右, 年平均降水量变化幅度在200—1000mm之间。由于流域整体较为湿润, 植被牧草发育良好, 是我国畜牧基地的典型代表^[25]。20世纪80年代以来, 新疆伊犁河流域土地开发活动趋于活跃, 三北防护林工程的实施、28个自然保护区的成立, 旅游业的发展等为流域带来了巨大的生态效益。

1.2 数据来源

伊犁河流域1980、1990、1995、2000、2005、2010、2015年的土地利用数据, 均来自中国科学院资源环境科学与数据中心(http://www.resdc.cn/)。各年份生态系统服务价值计算所用到的粮食产量、粮食价格等统计数据来源于《伊犁统计年鉴》、《伊犁哈萨克自治州统计年鉴》、《全国农产品成本收益资料汇编》。考虑到建筑用地对生态系统的正面效应远小于负面效应, 故本文不纳入生态系统服务价值。

2 研究方法 2.1 强度分析

强度分析可以定量的表达不同时间内各土地利用类型的变化强度, 采用空间层次分析法, 通过三个层次：时间间隔层次、地物类型层次、转移层次进行分析^[26]。表 1给出了本文使用强度分析所用到的数学符号及其含义。

(1) 时间间隔层次：探究每个时间段的年平均变化强度S_t与整个研究期内的总变化强度U相比, 变化速度的快慢。如果LUCC在研究期内是绝对稳定变化的, 那么S_t=U；如果研究期内LUCC变化是快速的, 那么S_t>U；如果研究期内LUCC变化是缓慢的, 那么S_t < U。时间间隔层次的公式为：

(1)

(2)

(2) 地物类型层次：探究在某时间段内, 某一土地利用类型的变化强度与年平均变化强度相比相对活跃还是平稳。如果某时间段内地类的增加或减少强度是绝对稳定的, 那么L_ti=S_t或G_tj=S_t。如果某时间段内地类的增加或减少强度是活跃的, 那么L_ti>S_t或G_tj>S_t；如果某时间段内地类的增加或减少强度是平稳的, 那么L_ti < S_t或G_tj < S_t。地物类型层次的公式为：

(3)

(4)

(3) 转移层次：为了区分在某时间段内, 不同土地利用类型之间的转变是否强烈。当地类n的转移强度R_tin>平均增加强度W_tn时, 说明n地类的增加可能来源于对i地类的占用, 当地类n的转移强度R_tin＜平均增加强度W_tn时, 说明n地类的增加避免于对i地类的占用；当转移强度Q_tmj>平均增加强度V_tm时, 说明地类m的减少可能转变为j地类, 当转移强度Q_tmj＜平均增加强度V_tm时, 说明地类m的减少未转变为j地类。转移层次的公式为：

(5)

(6)

(7)

(8)

在进行转换层次的强度分析时, 为了直观的分析不同地类在较多的时间节点上的变化情况, 本文在杨建新^[22]等人的研究基础上, 设计了面积增长和面积减少两种土地转换模式交叉列联表。如图 1所示, 每一横排①—⑤代表 5个时间间隔, 颜色代表地类j是否趋向于或避免于转换为地类i, 数值越大代表越稳定地趋向于转换为地类i, 颜色也就越深, 数值越小代表稳定地避免于转换为地类i, 颜色也就越浅。横向上填充相同的颜色说明研究期间j地类向i地类的转换较稳定；纵向上填充的颜色相同说明某时间间隔内地类j向地类i的转换较稳定。

2.2 景观指数选取

景观格局指数的选择要满足可以代表景观特征的不同方面、指数不具有冗余性和以其他相关研究的指数选择作为参考三个标准^[27-28]。通过我们对研究的需求, 为了全面反映流域景观格局特征, 分析整合了先前研究中选择的指标, 景观水平上选取斑块数量(NP)、平均斑块面积(AREA_MN)、香农多样性指数(SHDI)、香农均匀度指数(SHEI)、蔓延度(CONTAG)、聚集度指数(AI)；在景观类型水平上选取斑块密度(PD)、平均斑块面积(AREA_MN)、最大斑块指数(LPI)、景观形状指数(LSI)、散布与并列指数(IJI)。以上指标可以反应景观破碎度、形状面积、多样性、聚散性。

2.3 生态系统服务价值评估

基于土地利用类型的分类结果, 参照Costanza等的ESV评估方法^[9]来计算伊犁河流域1980—2015年的ESV。此外, 本文还采用了谢高地等^{[11, 28]}提出的生态系统系数修正方法：1个生态服务价值当量因子的经济价值为单位面积粮食产量的1/7^[11]。为消除不同年份农作物价格波动对总价值的影响, 以新疆2015年3种主要农作物(稻谷、小麦、玉米)的播种面积、产量和平均价格作为基础数据, 由式(9)计算出伊犁河流域单位面积农田生态系统粮食作物的经济价值为1728.01元/hm², 估算公式为：

(9)

其中, VC₀为伊犁河流域单位面积生态服务当量系数的经济价值, o_i为第i种粮食作物面积(hm²)；p_i为第i种粮食作物单产(kg/hm²)；q_i为第i种粮食作物的平均价格(元/kg)；M为3种粮食作物总面积；测算得出伊犁河流域各景观类型单位面积生态系统服务价值(表 2)。

ESV计算公式为：

(10)

(11)

其中, ESV代表区域生态系统服务总价值(元)；A_i代表第i种景观类型的面积(hm²)；VC_i代表第i种景观类型的价值系数[(元hm^-2a^-1)]；ESV_f代表单项生态系统服务价值(元)；VC_if代表单项服务价值系数[(元hm^-2a^-1)]。

3 结果与分析 3.1 强度分析 3.1.1 时间间隔水平上的强度分析

图 2中, 左侧区域表示在每个时间间隔上实际发生变化的地类面积占比情况, 右侧区域表示每个时间间隔内地类的平均变化强度(S_t), 即年均变化面积在总面积的占比情况。右侧区域的虚线表示整个研究时间间隔内的总变化强度(U)。由左侧区域可以看出, 伊犁河流域各土地利用类型在1980—2000年的变化面积占比较大, 三个时间间隔的变化面积分别为22.6%、17.44%和19.61%, 2000—2015年的三个时间间隔中, 各土地利用类型的变化面积占比明显减少, 分别是1.55%、0.49%和1.23%；右侧区域1980—2000年间三个时间间隔内的S_t呈增大趋势, 都要大于U, 说明变化是快速的。2000—2015年的三个时间间隔中, S_t都要小于U, 表现为慢速变化。整体来看, 35年来, 伊犁河流域的S_t呈现先上升后下降趋势。

3.1.2 地物类型层次的强度分析

图 3显示了6个时间间隔的地物类型层次强度分析。在1980—2015年期间, 研究区的地类变化较为复杂, 图 3的0刻度线左侧区域表示各地类的年均增加面积和减少面积, 右侧区域表示各地类的年均增加强度(G_tj)和减少强度(L_ti), 虚线表示对应时间段内, 各地类的平均变化强度(S_t)。可以看出, 1980—2015年, 草地的净减少面积最大, 耕地的净增加面积是最大的, 城镇建设用地的变化最小, 除了草地的面积减少外, 其余地类的面积均有所增加。土地面积变化幅度较大的是耕地(增加了3.66%)、城镇建设用地(增加了0.35%), 林地变化最小(增加了0.07%)；草地虽然减少量较大, 但其本身的基数大, 面积占比依旧较稳定。

年均变化强度以2000年为节点表现出了明显的阶段性, 2000年以前, 各地类的增减强度都较大, 除了草地, 其余地类的G_tj都非常活跃, G_tj前三位的地类是城镇建设用地、林地和水域。2000年以后, 各地类的活跃度相对降低, 且G_tj明显高于L_ti, 耕地和城镇建设用地成为流域内变化强度最活跃的地类, 林地和未利用土地的G_tj和L_ti趋于平稳, 水域的G_tj活跃度增加。总体来看, 耕地和城镇建设用地的G_tj在6个时间段内都大于S_t, 最为活跃, 草地的变化强度较平稳, G_tj始终小于S_t和L_ti。土地利用类型由变化剧烈趋向于平衡状态, 研究区城镇建设的发展也从盲目逐步到理性发展过度。

3.1.3 转换层次的强度分析

图 4、图 5显示的是转入(地类面积增加)和转出(地类面积减少)两种模式的转换层次交叉列联表, 伊犁河流域土地利用变化呈现出多种转换模式。在整个研究期间, 展现出稳定地倾向性系统转换模式的有：耕地-城镇建设用地、草地-林地(1980—2000年)、草地-耕地(2000—2015年)、水域-未利用土地(1980—2000年)、水域-耕地(2000—2015年)、城镇建设用地-耕地；展现出稳定地避免性系统转换模式的有：耕地-林地、林地-水域、草地-水域、水域-草地、城镇建设用地-林地、城镇建设用地-草地。

从转换层次强度分析的结果中可以得到, 水域面积的增长很大一部分是对耕地和未利用土地的占用, 但耕地和未利用地转出为水域的强度却在大部分时间间隔中表现出避免性, 这表明耕地、未利用土地转出为水域面积较大的原因主要是耕地和未利用土地的分布面积相对较大, 并且二者与水域的空间分布格局关系较为密切, 当水域面积增加时从二者转入的概率较大。

3.2 景观格局时空演变分析

本文选取的景观指数可以从破碎化、多样性水平和聚散性三个方面来描述1980—2015年伊犁河流域的景观异质性特征。

3.2.1 景观水平分析

从图 6可以看出伊犁河流域斑块数量(NP)在1980—2015年间呈先下降后上升趋势, 平均斑块面积(AREA_MN)趋势与NP相反, 蔓延度(CONTAG)和聚集度指数(AI)表征不同景观类型的聚集状态, 两个指数都呈先上升后大幅下降趋势, 说明流域内景观斑块的聚集性减弱, 破碎化程度加大。香农多样性指数(SHDI)和香农均匀度指数(SHEI)可以表征景观的多样性, 将二者结合分析发现流域内景观多样性在研究期内有所增加, 各景观类型区域均衡化发展。

3.2.2 景观类型水平分析

通过分析不同地类斑块的景观指数(图 7), 斑块密度(PD)和平均斑块面积(AREA_MN)不仅反映景观的异质性, 还与破碎度有一定的正相关性, 伊犁河流域水域和林地的PD在1980—2000年间波动较大, 早期经济发展的需求较大, 对水域和林地的破坏导致斑块破碎化程度加大, 后期随着西天山国家级自然保护区的建立和三北防护林等政策的实施, 使得破碎化程度有所减缓。耕地的PD逐渐增加, 景观形状指数(LSI)呈现先下降后上升趋势。最大斑块指数(LPI)显示草地景观在流域内占有极大的优势度, 由于人类活动干扰度增加, 城镇化发展水平持续上升, 对草地的开采程度加大, 导致许多大斑块的景观被分割为小斑块, LSI的上升趋势说明了草地景观的板块外形逐渐复杂化。城镇建设用地在1980—2000年间PD波动较大, 后期基本处于稳定状态。1980—2015年, 流域内各景观类型的散布与并列指数(IJI)都有增加趋势, 说明各地类景观与其周围相临近景观类型的丰富度有所增加。

3.3 生态系统服务价值变化分析

土地利用的变化也引起了伊犁河流域ESV的波动, 总体ESV变化如表 3所示。1980—2015年, 流域内总ESV呈先减少后增加的变化趋势, 2005年后稳定在998亿元左右。1980—1990年流域内水域过于活跃的减少强度和草地约1228km²的减少面积导致流域内ESV大幅下降, 1995—2000年ESV的快速回升除了耕地面积对草地和城镇建设用地的占用外, 更主要的是水域活跃的增加强度。水域由于其单位生态服务价值大, 导致ESV占比较大, 有效的补充了生态系统服务价值的损失。总的来看, 流域ESV在2000—2015年的增长值(1.72亿元)大于1980—2000年的增长值(-1.23亿元)。

由表 4可知, 伊犁河流域单项ESV中, 占前三位的是支持服务中的废物处理、水文调节、保持土壤。1980—2015年, 废物处理、水文调节、食物生产、美学景观、原材料生产的价值有所增长, 分别增长约3.38亿元、2.76亿元、2.17亿元、0.86亿元和0.30亿元, 其它单项ESV均呈减少趋势, 土壤保持亏损的最多, 约4.17亿元, 这与城镇建设用地的扩张有着直接的关系。所有单项ESV中, 食物生产的增幅最大, 达到7%, 这与伊犁河流域的耕地面积增加有关。

3.4 土地利用变化强度与景观指数、生态系统服务价值之间的关系

不同景观指数与各一级分类的ESV的相关性如表 5所示, 可以看出, 总ESV与NP、LSI、LPI、SHEI、SHDI呈正相关, 与AREA_MN、AI、CONTAG呈负相关。选择与每个一级ESV高度相关的景观指标(P < 0.01, P < 0.05)以构建多元回归模型(表 6), 供给服务、支持服务、文化服务都与所选择的景观指标显著相关, 可以看出SHEI对供给服务和支持服务的影响要大于其他景观指数, 文化服务与NP、LSI显著正相关, 与AREA_MN、AI显著负相关, 其中调节服务与各景观指数相关性较强, 这意味着优势景观斑块面积的增加、形状多样、各斑块聚集化发展会促使流域生态系统服务价值的增加, 景观斑块数量的增加、分布越分散, 即随着破碎化的增加, 使生态系统服务价值逐渐减少。为了进一步了解土地利用变化强度对生态系统价值的直接影响, 对伊犁河流域6个时间间隔的年均变化强度和7个年份的生态系统服务价值进行线性回归分析(图 8), 结果表明, 土地利用变化强度与ESV之间存在较强的正相关关系。

4 讨论 4.1 景观格局时空演变分析

伊犁河流域1980—2015年林地、草地面积减少的原因一部分是耕地和城镇建设用地的增加, 过度放牧活动占用了天然绿洲^[26], 另一方面是水资源的开发利用导致流域径流减少, 河道地下水位降低而引起了林草地的生态萎缩。1994年以来, 国家和地方政府相继出台了一系列政策法规, 例如《基本农田保护条例》, 在一定程度上规范了农田的退耕^[29]。2000年以后“国家级自然保护区”“三北防护林”等生态工程的成立, 在一定程度上保护了伊犁河流域的林地和水域资源, 带来了显著的生态、经济和社会效益^[30]。

目前流域内土地利用向着多样化方向发展, 通过强度分析, 可以帮助研究人员深入的了解流域内土地利用变化模式与过程, 并且对那些表现出活跃和稳定的关键变化模式给予更多关注^[31]。此外, 强度分析可以计算额外的指标来分析未知或不精确的土地利用数据^[32]。与那些仅通过传统遥感方法分析沿海地区动态的土地使用变化研究相比, 强度分析可以提供更多的信息, 像本文中使用了分辨率较小的土地利用类型数据, 强度分析模型的多层次分析就可以很好地弥补这一缺点。

4.2 土地利用、景观格局和生态系统服务价值的相互关系

土地利用变化强度的程度导致的景观格局的改变影响着生态系统的物质循环和能量流动, 导致流域内生态系统服务价值量的变化^[33-34]。Ayinuer ^[35]等分析了艾比湖自然保护区的ESV与景观格局的关系, 发现ESV与NP、PD和SHDI负相关, 与ED、SIDI、CONTAG呈正相关。岑晓腾^[36]以杭州湾南岸为例, 研究ESV与景观格局的关系, 发现土地利用越丰富、景观连通性和多样性越高, 对区域ESV有提升作用。从长远来看, 需要了解LUCC对景观格局变化的影响并确定ESV, 即了解景观格局与ESV之间的相关性, 由此可以将景观格局的变化情况作为指导区域ESV的可持续发展的关键因素。在对多元回归模型进行分析的基础上, 探索了ESV与景观指标之间的关系, 以量化景观变化对ESV的影响(表 6), 结果表明, NP、AREA_MN、AI均与各一级ESV呈负相关, 景观格局上的碎片化和分散性可能会对生态系统服务价值产生负面效应^[37]。

目前生态经济发展总体处于不协调水平, 为了继续提升伊犁河流域生态系统服务功能, 需对流域内自然资源和人类活动进行科学合理的规划管理^[38]。开垦荒地、合理放牧, 规整农田和牧草地, 以避免斑块趋于破碎化；控制城镇建设用地的扩展；保护流域内水资源, 科学的进行旅游景区和商业林用地的开发, 只有当人类活动与生态系统各项功能发生交集时, 才可以提升其价值, 为人类提供服务^[36]。

本研究从时间间隔层次、地类层次、转移层次对土地利用变化指标进行整合, 能够客观反应地类间的转换优势；基于当量因子法对ESV进行修正, 通过相关分析和回归分析法揭示伊犁河流域内土地利用变化和景观格局与ESV之间的相关程度。由于影像的分辨率可能会影响强度分析的结果或景观格局、ESV的变化, 因此需要根据Aldwaik等^[39]提出的框架进一步分析每个层次实际变化强度与标准变化强度的差异, 用更精细的分辨率数据来分析景观变化的研究。因为生态系统服务价值的异质性, 修正后的价值未必是最精确的, 并且研究区ESV的变化也受其它自然和社会因素的影响, 因此今后将会对其它驱动因子进行进一步剖析研究。

5 结论

通过土地利用变化强度分析, 系统性地递进式分析了1980—2015年间伊犁河流域土地利用变化过程与模式, 较为清晰直观地展示了不同层次上的土地利用类型变化。研究结果表明：

(1) 伊犁河流域35年来土地利用变化强度先增强后逐渐减弱, 其中耕地、水域和城镇建设用地出稳定性增长的趋势。在转移层次强度分析方面, 耕地向城镇建设用地展现出稳定地倾向性转换模式, 草地和水域、城镇建设用地展现出稳定地避免性转换模式。

(2) 流域内土地利用变化导致景观格局在研究期间斑块聚集性减弱, 破碎化加剧, 景观多样性有所增加, 各地类景观均衡化发展。1980—2015年, 伊犁河流域总生态系统服务价值增加了0.49亿元, 其中水域和林地是提供生态系统服务价值的主要地类, 共占流域总生态系统服务价值的47.94%。

(3) 伊犁河流域优势景观斑块面积的增加、形状多样化、各斑块聚集化发展会促使流域生态系统服务价值的增加, 土地利用变化强度与ESV存在明显的正相关关系。