城市化背景下淮海经济区生态质量时空变化及驱动因素

Spatiotemporal changes and driving factors of ecological quality in Huaihai Economic Zone under urbanization

摘要

目的 探明2000—2020年淮海经济区生态质量的时空演变规律，揭示生态变化的主要驱动因素及其交互作用机制，为区域生态保护与土地利用优化提供科学依据。 方法 基于Google Earth Engine平台，利用2000—2020年5个时点的MODIS影像构建遥感生态指数（RSEI），并结合单因子与交互因子探测模型，定量分析生态质量变化特征、驱动机制及土地利用转化效应。 结果 （1） 近二十年间淮海经济区生态质量总体改善，75.89%的区域RSEI显著提升，均值由0.470升至0.577，空间分布呈显著正相关，60.43%的区域仍存在退化风险；（2） 土地利用类型、降水和温度是主要驱动因子，人为因素亦有显著影响，其中土地利用与降水的交互作用解释力最高；（3） 耕地与建设用地的相互转化主导土地利用变化，快速城市化导致局部生态质量下降。 结论 淮海经济区生态质量总体向好但区域差异明显，土地利用及其与气候因子的交互作用是生态变化的关键机制，应加强土地集约利用与生态保护，实现区域生态系统的持续提升。

Abstract

Objective This study aims to clarify the spatiotemporal evolution patterns of ecological quality in the Huaihai Economic Zone （2000—2020）， reveal the main driving factors of ecological changes and their interaction mechanisms， and provide a scientific basis for regional ecological protection and land-use optimization. Methods Based on the Google Earth Engine （GEE） platform， the remote sensing ecological index （RSEI） was constructed using moderate resolution imaging spectroradiometer （MODIS） images at five time points during 2000—2020. Single-factor and interaction-factor detector models were applied to quantitatively analyze the characteristics of ecological quality changes， driving mechanisms， and the effects of land-use conversions. Results （1） Over the past two decades， ecological quality in the Huaihai Economic Zone has generally improved， with 75.89% of the area showing a significant increase in RSEI. The mean RSEI rose from 0.470 to 0.577， and its spatial distribution showed a significant positive correlation. However， 60.43% of the area still faced degradation risks. （2） Land use type， precipitation， and temperature were the main driving factors， and anthropogenic factors also had significant influences. Among them， the interactions between land use and precipitation had the highest explanatory power. （3） The mutual conversions between cropland and construction land dominated land-use changes， and rapid urbanization led to declines in local ecological quality. Conclusion The ecological quality of the Huaihai Economic Zone has generally improved， but regional differences are obvious. Land use changes and their interactions with climatic factors are the key mechanisms of ecological changes. It is necessary to strengthen intensive land use and ecological protection to realize the sustained improvement of regional ecosystems.

Graphical abstract

关键词

生态环境质量 / 遥感生态指数（RSEI） / 土地利用变化 / 淮海经济区 / 谷歌地球引擎（GEE）

Key words

ecological and environmental quality / remote sensing ecological index （RSEI） / land use change / Huaihai Economic Zone / Google Earth Engine （GEE）

健康、稳定的生态系统是人类生存与发展的良好基础^［1］。迈入21世纪后，地球正面临日益严重的生态问题，人类行为则是生态环境产生剧烈变化的主导因素^［2］。随着城市化进程不断加快，自然资源过度利用、植被退化、水污染、全球气候变暖以及生物多样性减少等问题日趋严重，因此，对区域生态环境进行客观、准确的评价尤为重要。

遥感技术以其全天候、大范围实时监测的优势^［3］而被广泛应用到生态环境领域^［4］。早期的生态质量评价大多采用单一遥感指标，例如反映农作物长势和营养信息的归一化植被指数（Normalized Difference Vegetation Index， NDVI），用于检测热环境的变化，尤其是在城市热岛效应和全球变暖研究中有重要作用的地表温度产品（Land Surface Temperature， LST），用于监测大气中悬浮颗粒物浓度的气溶胶光学厚度（Aerosol Optical Depth， AOD）等。然而，生态质量评价需要综合考虑多方面因素的影响，单一指标虽然可以很好地反映生态质量的某一方面，却很难全面反映区域整体生态质量^［5］。基于此便提出了多种基于遥感指数的方法，例如中国政府于2006年发布的《生态环境状况评价技术规范》提出了一种具有启发性的利用遥感技术的生态环境状况指数（Ecological Index， EI）。然而，在实际应用中，EI往往面临指标获取困难、数据更新周期长等^［3］问题。2013年，徐涵秋^［6］在EI的基础上提出了一种由主成分变换来集成各个指标并且权重不是人为确定的遥感生态指数（Remote Sensing Ecological Index， RSEI）。RSEI作为一种能够客观、有效地衡量生态系统抵御干扰和维持平衡能力并且反映生态系统质量的指标，自提出以来已被广泛应用在矿区开采监测，城市生态评价^［7-8］以及流域生态监测等^［3，9］研究领域。RSEI自提出以来也被国外学者广泛接受，Naseri等^［10］对伊朗卡拉季地区引入景观格局指数分析生态质量与绿地空间结构的关系，反映城市扩张下生态空间格局的显著变化；Boon等^［11］基于RSEI对俄罗斯萨马拉地区2010—2020年生态状况的时空演变进行了定量分析，发现生态变化与人类干扰程度密切相关。2019年，宋美杰等^［12］提出的改进遥感生态指数模型（Modified Remote Sensing Ecological Index， MRSEI）在RSEI基础上，参考了《生态环境状况技术评价规范》（HJ192-2015），增加了计算生物丰度的功能。但是徐涵秋等^［13］于2022年通过解释主成分分析的机理以及雄安新区的应用实例，指出MRESI的构建会导致各主成分分量相互抵消从而使结果失真。

传统的研究方法在进行长时间序列分析时面临着数据难以获取、数据量大、数据处理复杂的挑战。Google Earth Engine（GEE）平台克服了数据采集的困难和局限性，其强大的计算功能可以快速高效地处理图像数据集。通过GEE平台处理RSEI具有范围广、时间长、信息量完整等优点^［14-15］。同时，GEE平台的部分数据集，如Landsat和中等分辨率成像光谱仪（Moderate Resolution Imaging Spectroradiometer， MODIS）系列，已经经过预处理，将原始数据转换为表观反射率甚至地表反射率，使其适合进一步分析，而不需要专门的辐射和大气校正软件进行额外处理^［16］。

淮海经济区位于中国东部，是一个跨省区域合作的经济区，涵盖江苏、山东、安徽和河南四省的部分地区。该区域以徐州为中心，具有重要的区位优势，连接中国东部沿海和中西部内陆地区。淮海经济区在中国的区域发展战略中占据重要地位，既是东中部经济联系的枢纽，又在“一带一路”倡议中发挥着桥梁作用，促进了区域经济一体化和协调发展。近年来国务院批准并正式印发《淮河生态经济带发展规划》，使得淮海经济区在推进资源枯竭城市、老工业基地转型升级的同时，既要巩固提升全国重要粮食生产基地的地位，又要重点保护区域内生态环境。淮海经济区生态系统类型多样，既包括农田、林地、草地等自然生态单元，也包含城市、工业园区等强人类干扰区域，具有显著的生态异质性和土地利用变迁频繁的特点，生态系统面临着退化与恢复交织共存、空间分布差异显著等复杂问题。因此，采用RSEI这种基于遥感的综合指标，能在大尺度、长时间序列上高效反映区域生态质量的时空演变特征，具有良好的适配性与应用前景。

本研究利用GEE平台计算了淮海经济区的RSEI，并在GEE平台中集成实现Sen斜率计算与Mann-Kendall显著性检验流程，高效完成生态质量时序变化趋势与显著性分析，探讨生态环境质量空间分布格局背后的驱动机制；利用GEE平台的算力，通过像元级分析建立RSEI变化与土地利用转化的精细对应关系，揭示城市扩张、退耕还林等动态过程对生态系统的影响机制；为淮海经济区生态保护分区管控、土地集约利用等政策制定提供数据支持与决策依据。

1 材料与方法

1.1 研究区概况

淮海经济区西起114°48′E，东至119°48′E，南自33°6′N，北到36°12′N。地处中国东部，包括江苏省的徐州、连云港、宿迁，山东省的菏泽、济宁、临沂、枣庄，河南省的商丘以及安徽省的淮北、宿州，共10个地级市（图1）。区域总面积约9.6万km²，气候为温带季风气候，年平均气温大致在14 ℃~16 ℃，年降水量约为800 mm到1 200 mm，四季分明。该区域的地形主要以平原为主，区域北部的临沂和枣庄的部分地区存在丘陵和低山地形，沂蒙山脉在临沂市有分布，区域内部有多条重要河流，包括淮河和黄河支流等，同时也有洪泽湖等一些湖泊。产业结构以农业为基础，工业为主导，服务业和物流业快速发展，形成了较为多元化的产业体系，支撑着区域经济的增长。

1.2 数据来源及预处理

MODIS数据提供丰富的波段组合，包含36个光谱波段，覆盖从可见光到热红外的宽光谱范围，同时，其监测周期短、覆盖范围广，在大规模生态研究中得到广泛应用。其中，绿度指标提取自MOD13 A1.061数据集（https：∥developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD13 A1）。该数据集空间分辨率为500 m，时间分辨率为16 d，直接提供每个像素位置的NDVI。热度指标提取自MOD11 A2.061数据集（https：∥developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD11 A2）。该数据集的空间分辨率为1 000 m，时间分辨率为8 d，其每个像素值是该8 d周期内所有对应MOD11 A1-LST像素的简单平均值。湿度指标与干度指标提取自MOD09 A1.061数据集（https：∥developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD09 A1）。该数据集的空间分辨率为500 m，时间分辨率为8 d，该产品包含7个反射率波段、两个质量层和4个观测波段。在进行RSEI相关分析时，本研究选取植被生长季节即每年3月到11月的影像数据^［17］。由于数据集的空间分辨率有所不同，为保证分析结果的准确性，使用GEE平台对MOD11 A2.061进行重采样，确保各数据集在空间分辨率上的一致性，为后续分析提供基础。

土地利用数据采用GLC_FCS30D数据集^［18］。该数据集以30 m的空间分辨率全面揭示从1985年到2022年期间的精细土地覆盖类别动态产品，其中包含35个土地覆被子类别和26个时间步长，2000年之前每五年更新一次，之后每年更新一次。

为满足研究需要，结合淮海经济区土地覆盖类型特征，本研究将这35种类别重新划分为耕地、森林、草地、建设用地、湿地以及稀疏植被与裸露区域。由于研究区域内并没有永久冰雪覆盖地，故不对此类地物进行分类以及讨论。分类结果如表1所示。

经计算，稀疏植被与裸露区域面积不足研究区域面积的1%，故在本研究后续内容中不对此类型进行研究与分析。本研究对MODIS数据和土地利用数据进行分辨率统一，使用GEE平台的重采样功能将土地利用数据重新采样至500 m分辨率，从而保证数据分析的准确性和一致性。

本研究使用的DEM数据为SRTMDEMUTM 90 M分辨率数字高程数据产品，来源于中国科学院计算机网络信息中心地理空间数据云站点（http：∥www.gscloud.cn）。文中图片使用的投影均为WGS_1984_UTM_zone_49N（WKID：32649，权限：EPSG）。

1.3 遥感生态指数构建

RSEI由4个指标通过主成分分析（Principal Com-ponent Analysis， PCA）进行集成后产生，这4个指标分别为绿度指标、干度指标、湿度指标和热度指标^［19］。

1.3.1 绿度指标

绿度指标选用NDVI来表示。NDVI与植被生长状况、植物生物量、叶面积指数密切相关。其计算方法如下：

\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}=V_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}\times {\mathrm{10}}^{-\mathrm{4}}（1）

式中：V_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}为MOD13 A1.061数据集的NDVI波段值。

1.3.2 干度指标

干度指标（Normalized Difference Building and Soil Index， NDBSI）由建筑指数IBI和土壤指数SI构建而成。其计算方法如下：

\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{S}\mathrm{I}=\frac{\mathrm{I}\mathrm{B}\mathrm{I}+\mathrm{S}\mathrm{I}}{\mathrm{2}}（2）

其中，

\mathrm{I}\mathrm{B}\mathrm{I}=\frac{\frac{\mathrm{2}{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}{{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}-\left(\frac{{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}{{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+{\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}}}+\frac{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}}{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}+{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}\right)}{\frac{\mathrm{2}{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}{{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}+\left(\frac{{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}{{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+{\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}}}+\frac{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}}{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}+{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}\right)}（3）

\mathrm{S}\mathrm{I}=\frac{\left({\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+{\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}}\right)-\left({\rho }_{\mathrm{b}\mathrm{l}\mathrm{u}\mathrm{e}}+{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}\right)}{\left({\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+{\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}}\right)+\left({\rho }_{\mathrm{b}\mathrm{l}\mathrm{u}\mathrm{e}}+{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}\right)}（4）

式中：{\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}},{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}},{\rho }_{\mathrm{b}\mathrm{l}\mathrm{u}\mathrm{e}},{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}},{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}分别代表MOD09 A1.061影像的波段1，波段2，波段3，波段4，波段6。

1.3.3 湿度指标

湿度指标选用缨帽变换的湿度分量（Tasseled Cap Wetness Index， WET），能够较好地描述地表环境状况。其计算方法如下：

\begin{array}{l}\mathrm{W}\mathrm{E}\mathrm{T}=a_{\mathrm{1}}{\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}}+a_{\mathrm{2}}{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+a_{\mathrm{3}}{\rho }_{\mathrm{b}\mathrm{l}\mathrm{u}\mathrm{e}}+a_{\mathrm{4}}{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}+\\ a_{\mathrm{5}}{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{2}}}+a_{\mathrm{6}}{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}+a_{\mathrm{7}}{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{2}}}\end{array}（5）

式中： {\rho }_{\mathrm{r}\mathrm{e}\mathrm{d}},{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}},{\rho }_{\mathrm{b}\mathrm{l}\mathrm{u}\mathrm{e}},{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}},{\rho }_{\mathrm{N}\mathrm{I}{\mathrm{R}}_{\mathrm{2}}},{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}},{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{2}}}分别代表MOD09 A1.061影像的波段1，波段2，波段3，波段4，波段5，波段6，波段7。对于MODIS影像，a_{\mathrm{1}}\mathrm{到}{a}_{\mathrm{7}}的系数分别为0.114， 0.248， 0.240， 0.313，-0.312，-0.641，-0.508^［20］。

1.3.4 热度指标

热度指标选用LST。热量由热传感器数据表示，而MODIS LST产品数据可以直接表示。LST是经常被用于研究气候变化以及研究干旱、蒸散发、植被密度和地表能量平衡的重要指标^［21］，其计算方法如下：

\mathrm{L}\mathrm{S}\mathrm{T}=V_{\mathrm{L}\mathrm{S}\mathrm{T}}\times \mathrm{0.02}-\mathrm{273.15}（6）

式中：LST为地表温度，单位为℃；V_{\mathrm{L}\mathrm{S}\mathrm{T}}为MOD11 A2.061数据集的LST_Day_1 km波段值。

1.3.5 主成分分析

在计算完成4个指标后，通过主成分分析的方法对4个指标进行集成。使用主成分分析时各指标的权重通过其对主成分的贡献自动计算得出而非由人为设定，因此这种方法可以避免由于人为因素或方法差异所导致的权重设定误差，从而提高结果的客观性^［19］。由于这4个指标的范围与单位不同，故在进行主成分分析之前，需要将这4个指标的值归一化，将四者量纲统一到［0，1］。归一化的计算方法如下：

\mathrm{N}{\mathrm{I}}_i=\frac{I_i-I_{\mathrm{m}\mathrm{i}\mathrm{n}}}{I_{\mathrm{m}\mathrm{a}\mathrm{x}}-I_{\mathrm{m}\mathrm{i}\mathrm{n}}}（7）

式中：\mathrm{N}{\mathrm{I}}_i为归一化后的某一指标值；I_i为该指标在像元i处的值；I_{\mathrm{m}\mathrm{a}\mathrm{x}}为该指标的最大值；I_{\mathrm{m}\mathrm{i}\mathrm{n}}为该指标的最小值。

归一化以后可以进行主成分分析以计算第一主成分（PC1）。为使 PC1较大的值表示更好的生态条件，可用1减去PC1获得初始生态指数：

\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}}=\mathrm{1}-\left\{\mathrm{P}\mathrm{C}\mathrm{1}〔f\left(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I},\mathrm{W}\mathrm{E}\mathrm{T},\mathrm{L}\mathrm{S}\mathrm{T},\mathrm{N}\mathrm{D}\mathrm{B}\mathrm{S}\mathrm{I}\right)〕\right\}（8）

式中：\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}}为初始生态指数；PC1为第一主成分。

对初始生态指数进行归一化处理，将其值统一到［0，1］之间，得到RSEI：

\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}=\frac{\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}}-\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}\_\mathrm{m}\mathrm{i}\mathrm{n}}}{\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}\_\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}\_\mathrm{m}\mathrm{i}\mathrm{n}}}（9）

式中：分别\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}\_\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_{\mathrm{0}\_\mathrm{m}\mathrm{i}\mathrm{n}}为初始生态指数的最大值与最小值；RSEI即遥感生态指数，其值越接近1代表生态质量越好。

1.3.6 水体掩膜

为避免研究区域内的大范围水体对研究结果产生影响^［13］，本研究采用改进的归一化差异水体指数（Modified Normalized Difference Water Index， MNDWI）对研究区域内的大片水体进行掩膜处理^［22］。MNDWI的计算方法如下：

\mathrm{M}\mathrm{N}\mathrm{D}\mathrm{W}\mathrm{I}=\frac{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}-{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}{{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}+{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}}（10）

式中：{\rho }_{\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}},{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}{\mathrm{R}}_{\mathrm{1}}}分别代表MOD09 A1.061影像的波段4，波段6。本研究将阈值设置为0，即认定MNDWI值大于0为水体，需进行掩膜处理。

1.4 空间聚类分析

聚类分析能够对具有高或低值的元素进行空间聚类，并可以识别具有统计意义的热点和冷点。莫兰指数（Moran′s I）是衡量空间单元间相关性和异质性的一种常用空间自相关指标。本研究计算全局和局部莫兰指数，以评估淮海经济区生境质量空间分布的自相关性，并探究其空间集聚特征。计算方法如下：

I_G=\frac{N{\displaystyle \sum _{i=\mathrm{1}}^N}{\displaystyle \sum _{j=\mathrm{1}}^N}{w}_{ij}\left(x_i-\overline{x}\right)\left(x_j-\overline{x}\right)}{\left({\displaystyle \sum _{i=\mathrm{1}}^N}{\left(x_i-\overline{x}\right)}^{\mathrm{2}}\right){\displaystyle \sum _{i=\mathrm{1}}^N}{\displaystyle \sum _{j=\mathrm{1}}^N}{w}_{ij}}（11）

I_L=\frac{{\displaystyle \sum _{j=\mathrm{1}}^N}{w}_{ij}\left(x_j-\overline{x}\right)}{{\displaystyle \sum _{j=\mathrm{1}}^N}{\left(x_j-\overline{x}\right)}^{\mathrm{2}}}（12）

式中：I_G和I_L分别表示全局莫兰指数与局部莫兰指数；N为像元数；W_{ij}为空间权重矩阵；x_i和x_j分别为第i个区域和第j个区域的RSEI值；\overline{x}为研究区域RSEI的均值。全局莫兰指数取值范围是［-1，1］。局部空间自相关性可分为不显著、高—高聚集、低—低聚集、高—低聚集和低—高聚集这5种类型。

1.5 Sen斜率 & Mann-Kendall显著性检验的趋势分析

Sen斜率估计是一种稳健的非参数统计的趋势计算方法。该方法计算效率高，对于测量误差和离群数据不敏感，常被用于长时间序列数据的趋势分析中^［23］。Mann-Kendall是一种非参数统计检验方法，其优点是不需要测量值服从正态分布，也不要求趋势是线性的，并且不受缺失值和异常值的影响，在长时间序列数据的趋势显著检验中有着十分广泛的应用^［24］。Sen斜率计算方法如下：

{\beta }_{\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}=\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left(\frac{\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_j-\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_i}{j-i}\right)（13）

式中：{\beta }_{\mathrm{R}\mathrm{S}\mathrm{E}\mathrm{I}}为RSEI时间序列趋势的Sen斜率估计量；其值为正表示RSEI时间序列呈增加趋势，其值为负表示RSEI时间序列呈减少趋势；i和j为时间序列数，n为时间序列长度，有0<i<j<n；\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_i和\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_j分别为第i时刻与第j时刻的RSEI值。

Mann-Kendall非参数统计，检验方法为：

S={\displaystyle \sum _{i=\mathrm{1}}^{n-\mathrm{1}}}{\displaystyle \sum _{j=i+\mathrm{1}}^n}\mathrm{s}\mathrm{g}\mathrm{n}\left(\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_j-\mathrm{R}\mathrm{S}\mathrm{E}{\mathrm{I}}_i\right)（14）

\mathrm{v}\mathrm{a}\mathrm{r}\left(S\right)=\frac{n\left(n-\mathrm{1}\right)\left(\mathrm{2}n+\mathrm{5}\right)}{\mathrm{18}}（15）

（16）

式中：S为检验统计量；var（S）为S的方差；Z为标准化检验统计量。显著性水平对应的阈值分别是±1.65（对应90%的置信水平）、±1.96（对应95%的置信水平）和±2.58（对应99%的置信水平）。本研究选择95%的置信水平进行显著性检验。

1.6 Hurst指数

基于重标极差法计算的Hurst指数可以定量描述淮海经济区生态环境质量变化的持续性特征^［25］， Hurst指数取值范围为［0，1］，其中值为0.5表示时间序列是随机的，无长期依赖性；高于0.5表明正相关性，意味着过去的趋势在某种程度上可以预测未来的走势；低于0.5则表示负相关性，表明过去趋势与未来走势相反。

1.7 地理探测器

地理探测器是一套用来探测地理要素空间异质性及其驱动力的统计学方法^［26］。本研究利用因子探测器和交互探测器探测各因子对淮海经济区生态环境质量空间异质性的解释力以及各因子之间的关系。各因子对生态环境质量空间异质性的解释力以及因子间的关系可以用q值来度量。计算公式如下：

q=\mathrm{1}-{\displaystyle \sum _{h=\mathrm{1}}^L}{N}_h{\sigma }_h^{\mathrm{2}}/N{\sigma }^{\mathrm{2}}（17）

式中：N为整个区域的样本单元数量；N_h为子区域的样本单元数量；L为子区域数量；{\sigma }^{\mathrm{2}}为整个区域的方差；{\sigma }_h^{\mathrm{2}}为子区域的方差；q的取值范围为［0，1］， q值为1时表示因子完全控制生态环境质量的空间分布，q值为0时表示因子与生态环境质量空间分布无关。两个因子之间的交互作用有以下5种情况：（1），表明非线性减弱；（2），表明单因子非线性减弱；（3）q\left(X_{\mathrm{1}}\bigcap X_{\mathrm{2}}\right)>\mathrm{m}\mathrm{a}\mathrm{x}〔q\left(X_{\mathrm{1}}\right),q\left(X_{\mathrm{2}}\right)〕，表明双因子增强；（4）q\left(X_{\mathrm{1}}\bigcap X_{\mathrm{2}}\right)=q\left(X_{\mathrm{1}}\right)+q\left(X_{\mathrm{2}}\right)，表明两个因子相互独立；（5） q\left(X_{\mathrm{1}}\bigcap X_{\mathrm{2}}\right)>q\left(X_{\mathrm{1}}\right)+q\left(X_{\mathrm{2}}\right)，表明非线性增强。

2 结果与分析

2.1 淮海经济区RSEI时空变化

由于以五年为时间间隔可以较好地反映出长期变化趋势，并且不会忽略其中的细微变化，故本研究选取淮海经济区2000年、2005年、2010年、2015年、2020年的影像进行试验分析。

表2展示了各年份淮海经济区RSEI变化趋势，可以看出淮海经济区在2000—2020年期间的平均RSEI为0.519，其生态环境质量处于中等水平，同时总体呈现上升趋势，这表明区域生态环境在此期间有所改善。若以十年为周期，在2000—2010年这十年间，RSEI从0.470增加到0.533，年均增幅为0.063，增幅为13.4%，这一阶段的生态环境改善幅度较为显著；而在2010—2020年这十年间，RSEI从0.533增加到0.577，年均增幅为0.044，增幅为8.3%。虽然RSEI仍在增长，但增长速度较上一阶段有所放缓。

2.1.1 淮海经济区生态环境质量等级评价

根据现有研究^［27］，本研究将RSEI按等间隔分为5个等级，即：差 （0~0.2），较差 （0.2~0.4），中等 （0.4~0.6），良好 （0.4~0.6）和优秀 （0.8~1）（表3）。

从表3可以看出，从2000—2020年淮海经济区生态环境质量的总体趋势向好的方向发展。虽然在2005年和2010年期间有一定的波动，特别是“差”和“较差”等级的比例有所增加，但随着时间的推移，“良好”和“优秀”等级的比例显著提升。具体而言，在2000年时，淮海经济区有71.15%的区域生态环境质量处于“中等”水平，表现为生态环境较为单一，但改善潜力巨大；同时，仅有6.13%的区域达到 “良好”的标准，而“优秀”的比例更是微乎其微，仅为0.28%。2005年，虽然“差”和“中等”类别的面积比例分别上升至2.92%和59.86%，显示出部分地区的生态环境质量有所恶化，但“良好”与“优秀”类别的比例也有所上升，分别达到15.36%和0.38%。2010年，区域生态环境质量有显著改善，“良好”比例上升至32.73%，几乎是2005年的两倍，而“优秀”的比例也有所提升，达到1.03%。此时，“差”类别的区域大幅减少，比例降至1.32%。2015年，淮海经济区的生态环境质量进一步改善，“差”的比例下降至1.20%，“良好”的比例略有回落至30.59%，但“优秀”类别的比例上升至1.40%。到2020年，生态环境质量的改善更加明显，“差”和“较差”类别的区域比例分别下降至0.05%和8.49%，“良好”与“优秀”类别的区域比例大幅上升，分别达到40.53%和4.11%，这标志着大部分区域已进入“良好”及以上的生态环境等级。这表明，区域生态质量“由低向中、由中向高”的改善趋势明显。

2.1.2 淮海经济区生态环境质量变化趋势分析

2000—2020年淮海经济区RESI的Sen斜率 & Mann-Kendall显著性检验的趋势分析结果如表4和图2所示。近20年来，淮海经济区约80%的区域生态环境得到改善，其中显著改善的区域面积占比几乎达到区域面积的一半。生态质量下降的区域主要分布在济宁市西南部、徐州市中部、临沂市中部、商丘市西部。这一结果表明，淮海经济区在这20年间的生态环境质量总体向好发展，尤其是显著改善的区域分布广泛，反映出区域内的生态修复和环境治理成效显著。区域西南部的商丘市、中部的徐州市、东北部的临沂市未来生态质量变化趋势会像过去20年间一样持续改善，但通过GEE平台的计算，仅39.57%的区域未来生态质量变化趋势会与过去20年相同，不足四成，因此未来仍需注重对区域的生态保护工作。

2.1.3 淮海经济区生态环境质量空间分异特征

本研究对淮海经济区RSEI空间自相关性进行了全局和局部莫兰指数分析。2000—2020年淮海经济区全局莫兰指数均大于0，且均通过显著性检验，表明RSEI呈正向空间自相关性，并在空间上显著聚集。全局莫兰指数呈先升后降趋势，从2000年的0.764持续上升至2010年的0.773，再逐步下降至2020年的0.701。这一结果表明淮海经济区RSEI的空间聚集性呈先升后降趋势。

图3是通过ArcMap 10.8.2对矢量影像通过栅格转点后进行的局部莫兰指数分析结果，可以看出2000—2020年聚类分布变化不明显。高—高聚集与低—低聚集分布最为广泛，高生态质量的区域与低生态质量的区域空间相关性显著。高—高聚集区域主要分布在商丘市、徐州市、连云港市；低—低聚集区域主要分布在临沂市、枣庄市以及济宁市东北部；高—低聚集与低—高聚集区域在淮海经济区零星分布。

2.2 淮海经济区生态环境质量驱动因素分析

生态环境质量的空间分布特征是多种因素共同作用的结果。本研究选取了土地利用、海拔、坡度、温度、降水、GDP、人口密度7个因素，分析了它们对淮海经济区环境质量的解释力及主导因子。

因子探测的结果如表5所示，除GDP与人口密度外，所有因子对RSEI的空间异质性都有较强的解释力，其中，土地利用是RSEI空间异质性的主导因素，解释了RSEI过半的空间异质性。

图4是交互探测结果热力图，各因子间的相互作用增强了RSEI空间异质性的解释力，表明RSEI的空间异质性受多因子协同作用的影响。因子组合相互作用解释力最高的是土地利用与降水（0.625），其次是土地利用与温度（0.577），其中温度与降水的相互作用较为明显。总体而言，土地利用、温度、降水是影响淮海经济区RSEI空间分布的主要因素。

2.3 淮海经济区生态环境质量变化与土地利用变化的关系

根据表5结果，土地利用是淮海经济区生态环境质量变化的主导因素，本节对淮海经济区生态环境质量变化与土地利用变化的关系进行探究。考虑到RSEI不适合用于水体评价^［13］，本研究仅量化了耕地、林地、草地、建设用地的RSEI变化，结果如图5所示。可以看出，同一种土地利用类型与RSEI具有较好一致性，即建筑用地RSEI值偏低，森林、草地RSEI值居中，耕地RSEI值略高。

2.3.1 淮海经济区土地利用变化特征

近20年来淮海经济区土地利用分布如图6所示，耕地是淮海经济区的主要土地利用类型，占区域总面积的70%以上，建设用地总体上分散在不同的区域，并随时间推移不断从中心向周边扩张。草地与林地变化不明显。

2000—2020年淮海经济区土地利用变化情况见表6所示。20年间耕地向建设用地的转化面积最大，为6 260.09 km²。从这一现象不难看出城市化进程对耕地的长期压力。湿地和草地也受到较大的影响，湿地转化为耕地和建设用地的面积分别为216.04 km²和173.58 km²；草地向耕地和建设用地的转化分别为588.02 km²和178.11 km²。这表明除耕地外，其他自然土地类型也面临显著的土地利用变化压力。森林的变化则主要体现在向耕地的转化，共42.25 km²，这表明在城市化和农业扩展的过程中，森林资源仍然受到一定程度的侵占。

耕地的大量减少表明城市化进程对农业用地的压力在过去20年中持续存在。湿地和草地亦然，这两个土地类型在此期间也经历了显著的转化，主要是向耕地和建设用地的转化，反映了农业扩展和城市化对生态系统的压力。

2.3.2 不同土地利用类型转换导致的RSEI均值的变化

本研究量化了不同土地利用类型转换引起的RSEI值变化，结果如表7所示。在2000—2020年的整个时间段内，耕地、林地、草地向建设用地的转化均会在不同程度上引起RSEI的下降。其中，耕地转化为建设用地使RSEI下降9.6×10^-2，而草地转化为建设用地则使RSEI减少5.9×10^-3。与此相对，草地向耕地的转化在这期间带来较为显著的RSEI增加，达到0.156 3，林地向耕地的转化也导致0.157 7的增加。这表明，尽管在一些情况下土地利用转化对RSEI产生积极影响，但建设用地的扩张仍是导致RSEI下降的主要因素。

3 讨 论

淮海经济区作为老牌煤炭工业基地与主要粮食生产区^［28］，近20年来愈发注重生态环境保护工作。2024年年底，区域内各市生态环境局共同签署《淮海经济区生态资源共建共享协议》与《淮海经济区生态环境联合执法协作协议》更加彰显区域生态保护的协同趋势和重要程度，体现了政府部门于政策层面对生态环境问题的重视与应对意愿。本研究对了解淮海经济区生态质量变化具有一定参考价值，为区域可持续发展提供基础数据与决策依据。

本研究发现，近20年来，淮海经济区生态质量呈上升趋势，这一结果与过去的研究结论相符。尽管不同研究采用的遥感数据来源、评估模型等有所差异，但这些研究在趋势判断上表现出高度一致性，进一步印证RSEI在大尺度生态监测中的适用性。本研究则从整个淮海经济区的尺度进行系统分析，涵盖城镇、农业区、工业区等多种典型区域，在空间分布与驱动因素方面进行了更细致的刻画，表明区域生态保护措施和政策的实施可能对生态质量产生积极的影响。未来的研究可进一步探讨区域内不同城市生态质量变化的驱动因素及其机制，尤其是在政策实施、土地利用变化和气候变化等多重因素的综合影响下，生态质量变化的空间异质性及其内在机理仍有待深入分析。

地理探测器的因子探测结果中，坡度与海拔对RSEI的影响程度并不高，原因很可能是淮海经济区大部分处于平原地区，海拔与坡度总体变化不明显；而人为因素对RSEI影响程度低很可能是因为经济活动集中在城市核心区，而城市核心区分布如图6所示，集中在固定小部分区域，同时人口密度较高的地区可能既有高质量生态区域（如城市公园）也有低质量生态区域（如工业区）。

淮海经济区近20年来生态质量改善的主要机制可能包括以下几方面：一是大尺度的土地利用格局调整，特别是退耕还林、矿山修复等工程性措施，有效提高了区域生态系统的结构稳定性与生态服务功能；二是近年来区域降水增加可能在一定程度上为生态系统恢复提供了有利气候条件；三是以徐州为代表的城市在经济建设过程中强化绿地系统建设、构建生态廊道，使城市建成区内部生态质量得到提升。本研究同样发现，区域内生态质量在未来有可能会退化，未来需要严控建设用地无序扩张，特别是在耕地与林地集中区域，应优先利用存量建设用地，减少对自然生态空间的侵占。同时加强绿地建设，通过增设城市绿地、公园等手段，改善生态破碎格局，建立常态化生态监测机制，利用遥感平台每年更新生态质量评估结果，作为生态管控与土地规划的重要参考依据。

本研究还存在一些不足：（1） 地理探测器结果表明气温与降水对RSEI的解释力强，但未能进一步结合具体气候指标做分区分析，未来可引入气候分区方法或长期气象序列数据，细化气候对生态的影响机制。同时，淮海经济区由煤炭资源型区域向生态修复型区域转型的典型性决定了煤矿资源的空间分布与生态质量具有高度相关性，后续研究可以采取结合资源开采强度与修复投资数据的方法，探讨资源转型对区域生态系统的中长期影响。（2） 本研究在探究淮海经济区RSEI变化与土地利用变化的关系时，仅通过GEE平台进行逐像元的统计，没有采用统计建模方法进一步定量刻画两者的因果机制。未来可考虑引入多元回归、地理加权回归或机器学习方法构建RSEI预测模型，提高驱动因素识别的准确性与空间分辨率。 （3） 由于不同数据集的空间分辨率存在差异，重采样过程中不可避免地会带来一定的细节损失和数据同质化。未来的研究可以考虑使用相同空间分辨率的数据源，以减少重采样带来的误差，从而提高分析结果的精度和可靠性。

4 结 论

（1） 近20年间淮海经济区75.89%的区域生态质量得到改善。区域的平均RSEI值从2000年的0.470小幅上升至2020年的0.577。RSEI的空间分布呈正向空间自相关性，并在空间上显著聚集。60.43%的区域在未来存在生态环境质量退化的风险。

（2） 土地利用类型、降水、温度是影响淮海经济区生态环境质量的关键驱动因素，人为因素如GDP与人口密度也是重要影响因子。土地利用类型与降水的交互作用解释能力最强。

（3） 耕地与建设用地的相互转化是淮海经济区土地利用变化的主要形式。2000—2020年淮海经济区快速城市化导致大量耕地和林地转化为建设用地，这种转化导致RSEI下降。因此，推进土地集约利用，减少建设用地扩大对耕地和林地的占用，是保证区域可持续发展的关键。

参考文献

[1] Song X P， Hansen M C， Stehman S V， et al. Global land change from 1982 to 2016［J］. Nature， 2018，560（7720）：639-643.

[2] Baste I A， Watson R T. Tackling the climate， biodiversity and pollution emergencies by making peace with nature 50 years after the Stockholm Conference［J］. Global Environmental Change， 2022，73：102466.

[3] Xiong Y， Xu W， Lu N， et al. Assessment of spatial-temporal changes of ecological environment quality based on RSEI and GEE： a case study in Erhai Lake Basin， Yunnan Province， China［J］. Ecological Indicators， 2021，125：107518.

[4] Karbalaei Saleh S， Amoushahi S， Gholipour M. Spatiotemporal ecological quality assessment of metropolitan cities： a case study of central Iran［J］. Environmental Monitoring and Assessment， 2021，193（5）：305.

[5] 汪士为，吴伟.近31年嘉陵江流域生态环境质量时空演变及驱动因子探测［J］.水土保持研究，2024，31（1）：427-439.

[6] Wang S， Wu W. Spatiotemporal variation of ecological environmental quality and its response to different driving factors in Jialing River Basin in recent 31 years［J］. Research of Soil and Water Conservation， 2024，31（1）：427-439.

[7] 徐涵秋.区域生态环境变化的遥感评价指数［J］.中国环境科学，2013，33（5）：889-897.

[8] Xu H Q. A remote sensing index for assessment of regional ecological changes［J］. China Environmental Science， 2013，33（5）：889-897.

[9] Xu C， Li B， Kong F， et al. Spatial-temporal variation， driving mechanism and management zoning of ecological resilience based on RSEI in a coastal metropolitan area［J］. Ecological Indicators， 2024，158：111447.

[10] Huang H， Chen W， Zhang Y， et al. Analysis of ecological quality in Lhasa Metropolitan Area during 1990—2017 based on remote sensing and Google Earth Engine platform［J］. Journal of Geographical Sciences， 2021，31（2）：265-280.

[11] Li X， Mo X， Zhang C， et al. Unraveling the impacts of river network connectivity on ecological quality dynamics at a basin scale［J］. Remote Sensing， 2024，16（13）：2370.

[12] Naseri N， Mostafazadeh R. Spatial relationship of Remote Sensing Ecological Indicator（RSEI） and landscape metrics under urban development intensification［J］. Earth Science Informatics， 2023，16（4）：3797-3810.

[13] Boori M S， Choudhary K， Paringer R， et al. Spatiotemporal ecological vulnerability analysis with statistical correlation based on satellite remote sensing in Samara， Russia ［J］. Journal of Environmental Management， 2021，285：112138.

[14] 宋美杰，罗艳云，段利民.基于改进遥感生态指数模型的锡林郭勒草原生态环境评价［J］.干旱区研究，2019，36（6）：1521-1527.

[15] Song M， Luo Y， Duan L. Evaluation of ecological environment in the Xilin Gol Steppe based on modified remote sensing ecological index model［J］. Arid Zone Research， 2019，36（6）：1521-1527.

[16] 徐涵秋，邓文慧. MRSEI指数的合理性分析及其与RSEI指数的区别［J］.遥感技术与应用，2022，37（1）：1-7.

[17] Xu H Q， Deng W H. Rationality analysis of MRSEI and its difference with RSEI［J］. Remote Sensing Technology and Application， 2022，37（1）：1-7.

[18] Parastatidis D， Mitraka Z， Chrysoulakis N， et al. Online global land surface temperature estimation from Landsat［J］. Remote Sensing， 2017，9（12）：1208.

[19] Zhang J， Zhang P， Deng X， et al. Study on the spatial and temporal trends of ecological environment quality and influencing factors in Xinjiang Oasis［J］. Remote Sensing， 2024，16（11）：1980.

[20] Kumar L， Mutanga O. Google earth engine applications since inception： usage， trends， and potential［J］. Remote Sensing， 2018，10（10）：1509.

[21] 陈怀亮，刘玉洁，杜子璇，等.黄淮海地区植被生长季变化及其气候变化响应［J］.应用气象学报，2011，22（4）：437-444.

[22] Chen H L， Liu Y J， Du Z X， et al. The change of growing season of the vegetation in Huanghe-Huaihe-Haihe Region and its responses to climate changes［J］. Journal of Applied Meteorological Science， 2011，22（4）：437-444.

[23] Zhang X， Zhao T， Xu H， et al. GLC_FCS30D： The first global 30 m land-cover dynamics monitoring product with a fine classification system for the period from 1985 to 2022 generated using dense-time-series Landsat imagery and the continuous change-detection method［J］. Earth System Science Data， 2024，16（3）：1353-1381.

[24] 徐涵秋.城市遥感生态指数的创建及其应用［J］.生态学报，2013，33（24）：7853-7862.

[25] Xu H Q. A remote sensing urban ecological index and its application［J］. Acta Ecologica Sinica， 2013，33（24）：7853-7862.

[26] Lobser S E， Cohen W B. MODIS tasselled cap： land cover characteristics expressed through transformed MODIS data［J］. International Journal of Remote Sensing， 2007，28（22）：5079-5101.

[27] Xu H， Wang Y， Guan H， et al. Detecting ecological changes with a remote sensing based ecological index（RSEI） produced time series and change vector analysis［J］. Remote Sensing， 2019，11（20）：2345.

[28] 徐涵秋.利用改进的归一化差异水体指数（MNDWI）提取水体信息的研究［J］.遥感学报，2005，（5）：589-595.

[29] Xu H Q. A study on information extraction of water body with the modified normalized difference water index（MNDWI）［J］. Journal of Remote Sensing， 2005，（5）：589-595.

[30] 徐勇，黄雯婷，窦世卿，等.2000—2020年西南地区植被NDVI对气候变化和人类活动响应特征［J］.环境科学，2022，43（6）：3230-3240.

[31] Xu Y， Huang W T， Dou S Q， et al. Responding mechanism of vegetation cover to climate change and human activities in southwest China from 2000 to 2020［J］. Environmental Science， 2022，43（6）：3230-3240.

[32] 赖仕久，胡晋山，康建荣，等.煤炭资源型区域生态环境演变研究：以山西省为例［J］.自然资源遥感，2024，36 （4）：62-74.

[33] Lai S J， Hu J S， Kang J R， et al. Ecological evolution of coal resource-based regions： a case study of Shanxi Province［J］. Remote Sensing for Natural Resources， 2024，36（4）：62-74.

[34] Liu W M， Cheng Z Y， Li J， et al. Assessment of ecological asset quality and its drivers in Agro-pastoral Ecotone of China［J］. Ecological Indicators， 2025，170：113072.

[35] 王劲峰，徐成东.地理探测器：原理与展望［J］.地理学报，2017，72（1）：116-134.

[36] Wang J F， Xu C D. Geodetector： Principle and prospective［J］. Acta Geographica Sinica， 2017，72（1）：116-134.

[37] Zhou J B， Liu W Q. Monitoring and evaluation of eco-environment quality based on remote sensing-based ecological index（RSEI） in Taihu Lake Basin， China［J］. Sustainability， 2022，14（9）：5642.

[38] Tao J Y， Lu Y Q， Ge D Z， et al. The spatial pattern of agricultural ecosystem services from the production-living-ecology perspective： a case study of the Huaihai Economic Zone， China［J］. Land Use Policy， 2022，122：106355.

基金资助

国家自然科学基金面上项目(42471457)

在徐高校服务“343”产业发展项目(gx2023010)

江苏省研究生科研与实践创新计划(KYCX25_3181)