随着社会经济迅速发展并取得巨大进展[1-2]，人们对物质、环境等方面的需求不断增加[3]。生态系统在为人类社会提供食物、原材料和其他生态系统产品，以及气候调节和水土保持等生态系统服务的同时[4]，承受着社会经济发展带来的巨大压力。生态系统服务价值评估可以为区域可持续发展提供科学参考。通过评估生态系统产品和服务价值，定量核算其价值，可以加强对生态系统服务的了解，帮助了解生态系统对人类的贡献以及人类社会从生态系统中获得的利益[5]。
土地利用/覆盖变化（Land use and land cover change，LUCC）对生态系统服务的供给潜力有较大影响[6-7]。因此，研究土地利用变化的历史和未来变化轨迹及其对生态系统服务的影响，有利于对土地开发利用与保护做出科学决策，实现区域可持续发展。关于土地利用变化及其对生态系统服务影响的研究，主要集中在评估不同的土地利用类型在生态系统服务供给中的作用[8-9]，通过整合自然和经济因素，分析识别以往影响生态系统服务供给的驱动因素及机制[10]，或预测未来不同土地开发利用策略下土地利用变化引起的生态系统服务的演变趋势[11-12]。这些研究有助于了解生态系统服务的时间、空间异质性特征及其未来发展趋势，但其结果在预测和缓解未来问题方面指导性有待进一步增强。
本研究以黄河三角洲为研究区，以5年为时间间隔计算其土地利用面积变化、土地利用动态度以及土地开发利用强度，分析1985—2019年该区域土地利用变化的时空特征；进而计算其生态系统服务价值，分析土地利用变化背景下生态系统服务价值的响应变化特征；利用土地利用强度与生态系统服务价值构建耦合协调模型，分析其土地利用强度与生态系统服务协调程度的演变，为促进该区域的高质量发展提供参考。
1 材料与方法
1.1 研究区基本情况
研究区位于黄河入海口地带，地处中纬度，37°40′～38°10′ N，118°41′～119°16′ E，靠陆朝海，属于暖温带半湿润大陆性季风气候区。四季温差较为明显，年最高气温41.9 ℃，年最低气温-23.3 ℃，年平均气温11.7～12.6 ℃。土壤质地以轻壤土居多，新淤地较肥沃，新淤地淡质土层之下为老盐土和高矿化度的潜水，开采不当容易引起土壤的盐渍化。内陆地区土壤以盐碱地居多，沿岸淤地较为肥沃，水资源丰富，是重要的候鸟栖息地之一，具有生物多样性保护价值，也蕴含巨大的生态服务价值。
1.2 数据来源
本研究遥感数据来源于中国地理空间数据云（http://www.gscloud.cn），为1985、1995、2000、2005、2010、2015和2019年共7个时间段的Landsat影像。按照现有的土地利用分类标准以及土地利用特点，将研究区土地利用类型主要划分为耕地、草地、裸地、水体和建设用地5种。利用ENV I5.6软件对Landsat遥感数据进行FLAASH大气校正，通过监督分类得到各年份土地利用空间数据。生态系统服务价值换算涉及的社会经济统计数据来源于当地统计年鉴。
1.3 研究方法
1.3.1 土地利用动态变化及利用强度量化分析
（1）土地利用动态度评价。土地利用动态度能够量化某一土地利用类型的变化速率，可定量描述区域土地利用变化的速度，计算方法如式（1） 。
（2）土地利用强度评估。将土地利用强度定义为陆地生态系统经历由人类主导活动的影响强度。采用的土地利用强度计算公式参考文献[13]。该公式基于不同土地利用类型的面积以及单位面积土地使用能源和材料数量进行计算，已被证实是一种用于评估人类沿着连续梯度的累积效应的有效衡量指标。计算公式为每种土地利用类别的面积加权，如式（2） 所示。
表1 研究区各土地利用类型的利用强度系数 |
| 土地利用类型 | 土地利用强度系数（LDIi ） |
|---|---|
| 裸地 | 1.00 |
| 耕地 | 7.00 |
| 草地 | 3.41 |
| 建设用地 | 8.43 |
| 水体 | 1.00 |
1.3.2 生态系统服务价值估算及标准化处理
根据研究区实际修正生态系统服务价值的数值，利用统计年鉴数据计算2019年该区域单位面积粮食生产（小麦、水稻、玉米、小米、高粱、大豆和土豆）的年经济价值（1 605.83元/hm2），并以此作为研究区生态系统服务的单位经济价值。利用谢高地等[14]提出的“中国土地生态系统单位面积服务价值等值表”，估算研究区1985—2019年生态系统服务价值量。具体计算如式（3 ）～（4 ）。
式（3 ）～（4 ）中， 为生态系统服务价值，元；i为土地利用类型；j为生态系统服务类型； 为第i类土地利用类型的面积，hm2； 为第i类土地利用类型单位面积的生态系统服务价值，元/hm2； 为某类土地利用类型第j项生态系统服务价值当量（表2）；k为生态系统服务类型数量； 为单位生态系统服务的经济价值。
表2 生态系统服务价值分级标准 (元) |
| 生态系统服务 | 生态系统类型 | |||||
|---|---|---|---|---|---|---|
| 耕地 | 草地 | 建设用地 | 水体 | 裸地 | ||
| 合计 | 4.01 | 5.07 | 0 | 125.61 | 52.02 | |
| 供给服务 | 食物生产 | 0.85 | 0.10 | 0 | 0.80 | 0.51 |
| 原料生产 | 0.40 | 0.14 | 0 | 0.23 | 0.50 | |
| 水资源供给 | 0.02 | 0.08 | 0 | 8.29 | 2.59 | |
| 调节服务 | 气体调节 | 0.67 | 0.51 | 0 | 0.77 | 1.90 |
| 气候调节 | 0.36 | 1.34 | 0 | 2.29 | 3.60 | |
| 净化环境 | 0.10 | 0.44 | 0 | 5.55 | 3.60 | |
| 水文调节 | 0.27 | 0.98 | 0 | 102.24 | 24.23 | |
| 支持服务 | 土壤保持 | 1.03 | 0.62 | 0 | 0.93 | 2.31 |
| 维持养分循环 | 0.12 | 0.05 | 0 | 0.07 | 0.18 | |
| 生物多样性 | 0.13 | 0.56 | 0 | 2.55 | 7.87 | |
| 文化服务 | 美学景观 | 0.06 | 0.25 | 0 | 1.89 | 4.73 |
为统一比较，采用极差标准化法对数据进行消除量纲处理。标准化后，生态系统服务价值的取值范围均在0～1。使用ArcGIS的创建渔网工具，设置每个网格的尺度3.5 km×3.5 km，通过计算获得每个网格中土地的生态系统服务价值总和，使用极差标准化法计算出每个格子对应的归一化值，将归一化值的等级按照表3顺序分为5个等级。
表3 生态系统服务价值分级标准 |
| 生态系统服务价值标准化值 | 等级 |
|---|---|
| 0.6<ESV≤1.0 | 高 |
| 0.4<ESV≤0.6 | 较高 |
| 0.2<ESV≤0.4 | 中等 |
| 0.1<ESV≤0.2 | 较低 |
| 0<ESV≤0.1 | 低 |
1.4 耦合协调性分析
耦合概念源于物理学，意为两个或两个以上的系统通过相互的作用而互相影响，耦合协调度表示系统之间的相互作用程度[15]。为更好地反映土地利用综合集约度和生态系统服务价值，通过式（5） 计算。
$\begin{array}{c} S=\sqrt{U \times T} \\ U=2 \sqrt{U_{1} \times U_{2} / U_{1}+U_{2}} \\ T=0.5 U_{1}+0.5 U_{2} \end{array}$
表4 耦合协调度划分标准 |
| 类别 | 耦合协调度 | 等级 |
|---|---|---|
| 均衡发展 | 0.8<S≤1.0 | 高协调发展 |
| 0.6<S≤0.8 | 温和协调发展 | |
| 过度开发 | 0.5<S≤0.6 | 勉强协调发展 |
| 0.4<S≤0.5 | 轻微不平衡发展 | |
| 不平衡发展 | 0.2<S≤0.4 | 适度不平衡发展 |
| 0<S≤0.2 | 严重不均衡发展 |
2 结果与分析
2.1 土地利用变化及动态度分析
基于土地利用数据，统计1985—2019年研究区土地利用面积变化及动态度，如表5所示。从面积变化来看，1985—2019年该区域耕地、草地和裸地面积减少，水域和建设用地面积增加。从年际变化趋势来看，除建设用地一直呈增加趋势外，其他各类土地利用类型面积年际变化均呈现一定的波动特征，其中耕地、草地面积整体呈减少趋势，某一年际时段内面积增加，波动较小；裸地和水体面积变化年际波动较大，且均于1995—2000和2015—2019年两个时段出现波动，其中裸地在面积整体减少的背景下，面积呈现减少—增加—减少—增加的变化特征；水体面积呈现增加—减少—增加—减少的变化趋势。
表5 1985—2019年研究区土地利用面积年际变化及动态度 |
| 土地利用类型 | 1985年 | 1985—1995年 | 1995—2000年 | 2000—2005年 | 2005—2010年 | 2010—2015年 | 2015—2019年 | 2019年 | 1985—2019年 | |
|---|---|---|---|---|---|---|---|---|---|---|
| 耕 地 | 面积/hm2 | 467 007.66（63.55%） | -3 126.60 | 3 609.90 | -8 874.27 | -21 290.40 | -21 974.94 | -6 925.86 | 408 425.49 （55.58%） | -58 582.20 |
| 动态度/% | -0.07 | 0.16 | -0.38 | -0.93 | -1.00 | -0.42 | ||||
| 草 地 | 面积/hm2 | 20 956.86 （2.85%） | -8 019.00 | -3502.08 | -2 343.42 | 1 665.72 | -3 207.87 | -847.08 | 4 703.13 （0.64%） | -16 253.70 |
| 动态度/% | -3.83 | -5.41 | -4.97 | 4.70 | -7.33 | -3.82 | ||||
| 水 体 | 面积/hm2 | 41 646.15 （5.67%） | 5 232.87 | -7 221.15 | 18 277.83 | 23 428.17 | 38 929.05 | -9 783.99 | 110 508.93 （15.04%） | 68 862.78 |
| 动态度/% | 1.26 | -3.08 | 9.22 | 39.96 | 9.57 | -2.03 | ||||
| 裸 地 | 面积/hm2 | 157 743.09 （21.47%） | -9 384.75 | 1 517.04 | -8 142.48 | -11 795.40 | -28 928.88 | 9 469.08 | 110 477.70 （15.03%） | -47 265.40 |
| 动态度/% | -0.60 | 0.20 | -1.09 | -1.66 | -4.45 | 2.34 | ||||
| 建 设 用 地 | 面积/hm2 | 47 524.14 （6.47%） | 15 297.48 | 5 596.29 | 1 082.34 | 7 991.91 | 15 182.64 | 8 087.85 | 100 762.65 （13.71%） | 53 238.51 |
| 动态度/% | 3.22 | 1.78 | 0.32 | 2.30 | 3.92 | 2.18 | ||||
土地利用动态度变化趋势与土地利用面积变化基本一致。此外由动态度结果可知，不同土地利用类型变化速率的差异，耕地整体动态度数值较小，年际变化小；水体的动态度数值较大，年际变化大；其他土地利用类型的动态度数值居中，年际变化中偶有年份变异较大。
由1985—2019年研究区土地利用分布（图1）可知，内陆主要以耕地和建设用地为主，沿海地区主要为水域、裸地，草地集中分布于黄河入海口附近。研究区内陆地区的土地利用变化主要表现为建设用地面积增加。此外，值得注意1995—2000年耕地的面积有所增加，黄河入海口处发育有广阔的冲积平原，冲积平原形成后首先形成裸地，随后被开发利用为耕地。沿海地区土地利用变化主要表现为水域面积扩大，裸地面积减小，且西北部海岸新增水域面积大于东南部海岸，结合实地调查发现，水体面积的增加与沿岸海水养殖面积扩大有关。此外，集中分布于东北部海岸黄河入海口区域的草地面积减少，部分草地逐渐转变为裸地和耕地。
2.2 土地利用强度演变分析
基于土地利用数据计算得到网格尺度上土地利用强度，由图2可知，研究区由沿海向内陆土地利用强度逐渐增加。内陆地区除黄河开发强度较低外，整体开发强度较高；研究区外围沿海地区的开发强度整体较低，且低强度区域集中分布于西北部海岸和东南部海岸，这些区域的土地利用类型以裸地和水域为主，而东北部海岸耕地分布广泛，故土地利用强度较高。1985—2019年研究区土地利用强度发生变化的区域主要位于沿海地区，表现为西北部海岸土地利用强度降低，东北部海岸土地利用强度增加。
2.3 生态系统服务价值变化分析
基于土地利用数据换算生态系统服务价值，进而统计1985—2019年研究区生态系统服务价值年际变化。由表6可知，该时段研究区生态系统服务价值总体呈增加趋势，服务价值由1985年的247.55亿元增加为341.88亿元，增长了38.11%。从生态系统服务类型来看，供给服务和调节服务价值增加，其中调节服务增加比例最大（55.89%），供给服务价值增加比例较小（25.19%）；支持服务和文化服务降低，且从数值来看服务价值减少较小。此外，从年际变化特征来看，生态系统服务价值变化呈现波动性变化特征，总生态系统服务价值表现为增加—减少—增加—减少的波动特征。不同生态系统服务类型中，供给和调节服务年际变化波动较大，与总服务价值波动特征一致，支持服务和文化服务波动较小，呈减少—增加变化特征。
表6 1985—2019年研究区生态系统服务价值年际变化 (亿元) |
| 服务类型 | 1985年 | 1985—1995年 | 1995—2000年 | 2000—2005年 | 2005—2010年 | 2010—2015年 | 2015—2019年 | 2019年 | 1985—2019年 |
|---|---|---|---|---|---|---|---|---|---|
| 供给服务 | 24.98 | 0.14 | -0.94 | 2.07 | 2.40 | 3.69 | -1.06 | 31.28 | 6.30 |
| 调节服务 | 170.16 | 3.80 | -12.14 | 27.86 | 35.00 | 53.15 | -12.55 | 265.27 | 95.11 |
| 支持服务 | 38.63 | -1.49 | -0.15 | -0.54 | -1.03 | -3.11 | 0.86 | 33.17 | -5.46 |
| 文化服务 | 13.78 | -0.59 | -0.11 | -0.08 | -0.20 | -1.05 | 0.41 | 12.16 | -1.62 |
| 总价值 | 247.55 | 1.86 | -13.35 | 29.31 | 36.17 | 52.68 | -12.34 | 341.88 | 94.32 |
由图3可知，研究区生态系统服务价值内陆低、沿海高，呈现由内陆向沿海逐渐增加的空间分布特征。从时间演变来看，沿海地区生态系统服务价值增加，内陆地区基本保持不变。分区域来看，研究区西北部海岸的生态系统服务价值增加明显，1985—2019年由中等服务价值转变为高生态系统服务价值。1985—2019年部分区域由低生态系统服务转变为较高和高生态系统服务价值；东南部海岸大部分为较高等级服务价值；部分区域由高生态系统服务转变为较高生态系统服务价值。内陆地区除黄河生态系统服务价值稍高外，生态系统服务价值均较低，且1985—2019年服务等级变化不大。
2.4 耦合协调性分析
3 结论与讨论
本文以5年为时间间隔，揭示研究区土地利用与生态系统服务动态演变特征，结果表明，1985—2019年研究区土地利用变化与生态系统服务演变呈现较强的时空异质性。时间尺度上，土地利用变化及生态系统服务演变过程复杂、频繁，波动特征明显。空间分布上，研究区内陆地区以耕地和建设用地为主，整体土地利用强度大，土地利用变化过程表现为建设用地面积增加，生态系统服务价值相对较低；沿海地区土地利用类型以水体、裸地和草地等自然景观为主，土地利用强度较低，土地利用变化特征表现为水体面积增加，裸地和草地面积减少，其中西北部海岸水域面积增加明显，生态系统服务价值明显提升。结合土地利用变化分析，研究区耕地生态系统服务价值降低，水体生态系统服务价值增加，这与Wang等[3]和Zhang等[16]在黄河三角洲区域开展的生态系统服务演变研究结果基本一致。Zhang等[16]研究指出，随着该地区沿海水产养殖业的不断发展，越来越多的沿海滩涂被开发利用，未来该区域内的水域面积还会不断增加。此外，区域内草地具有水土保持、降尘减尘以及生物多样性维护等多种重要生态功能，草地开发为耕地会导致影响生态系统功能[17]。因此，从自然资源保护与生态系统服务持续提升的角度出发，需要特别关注区域内未利用土地的合理开发、滩涂陆海生态连通性的提升和高效生态种养模式的应用推广[18]。
1985—2019年研究区土地利用强度与生态系统服务价值耦合协调程度从以温和协调发展和勉强协调发展为主转变为以温和协调发展和高协调发展为主，表明研究区土地资源开发利用与生态环境供给的关系更加协调向好。类似的研究结果也出现在黄河三角洲城市化与生态系统服务价值耦合协调关系研究中。卢辞等[19]研究发现，2000年以来黄河三角洲城市化与生态系统服务价值耦合协调度呈上升趋势，且沿海地区耦合协调指数增幅较大。土地利用变化和土地利用强度变化较明显的区域，生态系统服务价值演变较剧烈的区域，以及耦合协调性转变明显的区域，均位于研究区沿海地区以及沿海向内陆的过渡地带。这些区域受自然条件和人为因素影响，环境较敏感，稳定性有待进一步提高，需要高度重视该区域资源开发利用与环境保护之间的关系[20]。此外，研究区内陆地区的黄河沿岸作为生态系统服务的热点区域之一，其环境保护问题同样需要特别关注。
在土地利用变化背景下，土地利用强度与生态系统服务价值的耦合协调程度可以反映区域土地资源开发与区域资源环境供给的关系。土地利用变化时，土地利用强度及生态系统服务供给均会做出响应，同时土地资源开发与区域资源环境供给的关系也会发生转变。生态环境保护成为高质量发展的主要任务之一，为此，有必要协调资源开发与资源环境的关系，增加土地利用与生态系统服务供给的耦合协调程度，实现区域可持续发展。
本文分析了1985—2019年研究区土地利用变化背景下，土地利用强度、生态系统服务价值的时空变化特征，计算两者的耦合协调性，主要研究结论如下。
（1）1985—2019年研究区土地利用变化较剧烈，表现为建设用地和水体面积增加，裸地、草地和耕地的面积呈减少趋势，其中耕地整体动态度数值较小，年际变化小，水体的动态度数值较大，年际变化大。该区域土地利用强度表现为由沿海向内陆逐渐增加的分布特征。
（2）研究区土地利用变化明显影响其生态系统服务时空分异，表现为内陆地区以耕地和建设用地为主，1985—2019年土地利用变化表现为建设用地面积扩张，生态系统服务价值相对较低；沿海地区土地利用类型主要为水体、裸地和草地，1985—2019年土地利用变化特征为水体面积扩大，裸地和草地面积减少，其中研究区西北部海岸水域面积增加明显，生态系统服务价值明显提升，东南部海岸黄河入海口区域草地面积减少。
（3）1985—2019年研究区生态系统服务价值与土地利用强度间的耦合协调性总体呈上升趋势。黄河及周边区域一直表现为高协调发展，新增高协调发展区域主要位于研究区内陆向海岸过渡的区域。