城市化影响下北京市植被物候变化特征

Spatiotemporal Dynamics of Vegetation Phenology in Beijing

摘要: 植被物候是环境变化的敏感指示器，植被物候的时空演化及其对城市化的响应研究有助于优化城市空间与绿地格局。基于MOD13Q1 EVI数据，本文采用S-G滤波进行时间序列重建，并借助动态阈值法提取物候参数，解析2002~2021年北京市物候的时空变化，同时在城–郊–村分区尺度上评估其对城市化强度的响应。结果表明：20年来全市植被物候呈现SOS提前(0.72 d∙a⁻¹)、EOS推迟(0.55 d∙a⁻¹)、GSL延长(1.26 d∙a⁻¹)的趋势；城市SOS最早(97.58 d)、EOS最晚(308.06 d)、GSL最长(210.02 d)。物候年际变化幅度呈梯度特征：城区最小(0.30 d∙a⁻¹)、郊区次之(1.07 d∙a⁻¹)、农村(1.41 d∙a⁻¹)更大。不透水面(ISP)增加普遍导致SOS提前、EOS推迟与GSL延长；在城市，郊区，农村三个分区中，ΔISP每增加10%，SOS分别提前0.96 d、0.97 d、0.56 d；EOS分别推迟0.28 d、0.31 d、0.26 d；GSL分别延长1.24 d、1.30 d、0.82 d，物候对城市化的响应具有明显的分区特征；将城郊村分区与不透水面占比结合，能够观察到物候对城市化敏感性的差异，从而提升物候对城市化响应的准确性。

Abstract: Vegetation phenology is a sensitive indicator of environmental change; elucidating its spatiotemporal evolution and responses to urbanization can inform urban spatial planning and green-space configuration. Using MOD13Q1 EVI, we reconstructed time series with a Savitzky-Golay filter and extracted phenological metrics via a dynamic-threshold method to characterize phenological changes in Beijing from 2002~2021, and assessed responses to urbanization intensity across urban-suburban-rural partitions. Citywide phenology exhibited earlier start of season (SOS; 0.72 d∙a⁻¹), later end of season (EOS; 0.55 d∙a⁻¹), and longer growing-season length (GSL; 1.26 d∙a⁻¹). In urban areas, SOS occurred earliest (DOY 97.58), EOS latest (DOY 308.06), and GSL was longest (210.02d). Interannual variability showed a clear gradient—smallest in urban zones (0.30 d∙a⁻¹), intermediate in suburban zones (1.07 d∙a⁻¹), and largest in rural zones (1.41 d∙a⁻¹). Increasing impervious surface percentage (ISP) generally advanced SOS, delayed EOS, and lengthened GSL. For each 10% increment in ΔISP, SOS advanced by 0.96 d, 0.97 d, and 0.56 d; EOS was delayed by 0.28 d, 0.31 d, and 0.26 d; and GSL extended by 1.24 d, 1.30 d, and 0.82 d in urban, suburban, and rural areas, respectively. These results demonstrate marked partition-dependent sensitivities of phenology to urbanization; combining the urban-suburban-rural framework with ISP effectively reveals heterogeneous sensitivities and improves attribution of phenological responses to urbanization.

1. 引言

植被物候泛指植物在环境条件驱动下出现的关键生长现象，包括萌芽、展叶、开花、结果以及叶片凋落等[1]，主要包括生长季开始期(start of season, SOS)、结束期(end of season, EOS)和生长季长度(growing season length, GSL)，是表征植被动态与气候变动的核心生物学指标[2]-[5]。植被物候通过调节碳循环，水分循环等生态过程，对城市的生态环境质量和可持续发展产生重要影响[3] [5]。同时，城市化进程的加快导致不透水面增加，改变城市地表的显热与潜热通量[6]，导致城市植被的SOS提前，EOS推迟[7] [8]。因此，在加速城市化背景下，监测城市植被物候的时空格局并解析其对城市化的响应机理，对推进城市生态保护与可持续发展具有重要理论意义[9]。

当前监测植被物候的方法主要包括地面监测、近地面监测以及卫星监测[10]。其中，卫星监测凭借其成本小，多尺度的优势，被广泛应用于植被物候的监测与分析[11] [12]；卫星遥感通过MODIS [12]、Landsat [13]和Sentinel-2 [14]等平台，为植被物候的监测提供了可靠支持。归一化植被指数(NDVI) [15]、增强型植被指数(EVI) [14]和日光诱导叶绿素荧光(SIF) [16]是植被物候监测中广泛使用的植被指数；EVI因其在高植被覆盖区表现更稳定，被认为更适用于城市植被动态监[17]。然而，植被指数通常会因为大气扰动，传感器噪声等导致数据缺失或异常值，需要进行重构拟合以平滑数据并去除噪声[18]，常采用双Logistic函数法[11]、Savitzky-Golay (S-G)滤波法[19]等方法重构植被指数的时间序列；其中S-G滤波法能够保持时间序列原始结构特征，适合中高频噪声的去除[20]。在物候参数提取过程中，传统的导数法[21]和固定阈值法[22]在动态背景下存在不稳定性，而动态阈值法能够根据年内植被指数波动特征自适应确定物候拐点，具备更高的准确性[20]。因此，结合EVI指数，使用S-G滤波法和动态阈值法提取植被物候，已广泛应用于城市植被物候监测研究中[23] [24]。

随着城市化进程的不断推进，城市化对植被物候的影响逐渐成为物候研究的关注重点[25]。传统缓冲区方法常以城市边界为起点向外划分固定距离的缓冲区，通过对比缓冲区间物候均值差异分析植被物候对城市化的响应[7] [11]。已有研究基于缓冲区方法发现北京市城市地区的植被物候相比农村地区年均SOS提前约12.2 d，年均EOS推迟约18.9 d [11] [26]。然而，缓冲区内不同区域在开发程度和生态环境等方面存在差异，区域之间物候差异显著[27]，在精细尺度上研究植被物候差异时，缓冲区方法难以满足研究需求。为弥补上述不足，城市化强度法被提出[28]。该方法基于规则网格内不透水面百分比(Impervious Surface Percentage, ISP)将城市空间划分成0%~100%的100个连续梯度，依据ISP大小划分城乡梯度，以描述不同区域植被物候差异[12] [24]。目前尚无研究采用城市化强度法分析北京市植被物候对城市化的响应特征。然而，依据ISP大小划分城乡梯度容易产生误差，研究表明，许多ISP较低的网格包含受城市化影响的农村过渡带等区域，并非真实的乡村区域[29]，这会导致城乡物候差异被弱化。因此，本文提出在城市化强度法的基础上引入城乡空间分区策略，以期减小城乡判别偏差、并更清晰刻画植被物候对城市化变化的响应。

基于此，本研究选用MODIS的MOD13Q1 EVI数据集作为主数据源；利用S-G滤波和动态阈值法提取北京市植被物候参数SOS、EOS和GSL，分析植被物候时空变化特征采用S-G滤波与动态阈值法识别北京市植被物候参数(SOS、EOS、GSL)，并刻画其时空变化特征；利用“城–郊–村”范围数据集(Urban-Suburban-Rural Triad Structure Dataset, USR)进行城乡分区，在各分区中利用全球不透水面数据集(GAIA)量化城市化强度，探究北京市不同城乡分区上植被物候的时序变化特征及其对城市化的响应，为北京市生态环境保护和可持续发展战略制定等方面提供科学依据。

2. 研究区概况

北京市(39˚26'~41˚03' N, 115˚25'~117˚35' E)位于华北平原北部边缘，地形西北高、东南低，西北部为燕山、太行山，东南部为平原区，是北京主要的城市发展空间(图1(a))。属暖温带半湿润季风气候，夏热多雨、冬冷干燥。北京市林地主要分布在西部和北部山区；耕地集中于东南部平原地带。北京市作为中国经济中心，城市化进程快、城市热岛效应显著，对城市植被物候产生影响，因此研究北京市植被物候具有重要意义。

注：该图基于自然资源部标准底图服务网站下载的审图号为GS(2024)0650号的标准地图制作，底图无修改。

Figure 1. Spatial distribution of elevation (a) and land-use/land-cover (LULC) types (b) in Beijing

图1. 北京市高程分布(a)及土地利用类型(b)

3. 数据和方法

3.1. 研究数据

本研究采用TIMESAT 3.3软件[30]提取北京市植被物候参数鉴于物候识别需以“三年滑动窗口”获取中间年度的结果，选取2001~2022年MOD13Q1数据集中的EVI作为输入。利用MODIS Reprojection Tool (MRT)对覆盖区影像开展批处理拼接、投影转换与范围裁切，生成2002~2021年北京地区的EVI时序，用于后续物候指标提取。GAIA为全球不透水面数据集，提供全球范围内每年人工不透水面的空间分布和变化情况，可用于2002~2021年北京市城市化强度的表征[31]。USR数据集选择Xiong等编制的2002~2021年1 km的“城–郊–村”范围数据集(USR) [32]，用于划分北京市2002~2021年每年的城市，郊区，农村分区(图2)，该数据通过夜间灯光分位数突变点检测划分每年的城市，郊区，农村区域，具有全国覆盖、长期一致性和高可靠性，适用于城市化、碳排放等环境研究。为消除不同数据源的分辨率差异，所有栅格数据在ArcGIS中统一重采样至250 m空间分辨率，用于后续分析。

Figure 2. Spatial distribution of urban, suburban, and rural areas in Beijing in 2021

图2. 北京市2021年城市郊区农村分布图

本研究使用的数据包括来自中分辨率成像光谱仪(MODIS)的 MOD13Q1、MCD12Q2产品数据集，以及不透水面(GAIA)、“城–郊–村”范围数据集(USR)，具体数据来源见表1。

Table 1. Data sources and applications

表1. 数据来源及用途

数据集	空间分辨率	时间	来源	用途
MOD13Q1 EVI数据	250 m	2001~2022	https://search.earthdata.nasa.gov/search	提取植被物候关键参数
MCD12Q2	500 m	2002~2021	https://search.earthdata.nasa.gov/search	植被物候参数验证
“城–郊–村”范围数据集(USR)	1 km	2002~2021	https://ieeexplore.ieee.org/ielx7/4609443/4609444/10352603	城市空间分区
全球不透水面数据集(GAIA)	30 m	2002-2021	https://data.ess.tsinghua.edu.cn/	构建城市化强度梯度

3.2. EVI时间序列重建

本研究选用Savizky-Golay (S-G)滤波法对北京市的EVI时序数据进行拟合重建。在原始植被特征保真方面表现良好，对人类活动密集区的适用性更高[33]。其公式为：

$Y_j^{\ast }=\frac{{\displaystyle \sum {}_{i=-n}^n}{C}_i{Y}_{j+i}}{N}$ (1)

式中：$Y_j^{\ast }$ 为时刻j上经过滤波平滑后的EVI值，$Y_{j+i}$ 为原始EVI时间序列，以j为中心，左右第i个位置的观测值，C_i为第i个位置的滤波器系数，n为滤波窗口宽度的1/2；N为滤波器的长度，N = 2n + 1，表示滑动窗口的总宽度。考虑到EVI的季节周期性，设置S-G滤波器参数为：n = 4，N = 9。

3.3. 植被物候参数提取

本研究使用动态阈值法提取物候参数，其原理为根据植被的生理特征，以植被指数(VI)曲线在上升段与下降段达到振幅某一固定占比的时刻，分别界定为SOS与EOS。使用较为灵活，其公式为：

$\text{SOS}={\text{EVI}}_{\text{MIN1}}+\left({\text{EVI}}_{\text{MAX}}-{\text{EVI}}_{\text{MIN1}}\right)\times \alpha$ (2)

$\text{EOS}={\text{EVI}}_{\text{MIN2}}+\left({\text{EVI}}_{\text{MAX}}-{\text{EVI}}_{\text{MIN2}}\right)\times \beta$ (3)

$\text{GSL}=\text{EOS}-\text{SOS}$ (4)

式中，EVI_MAX表示EVI在全年时间序列最大值，EVI_MIN1表示曲线上升阶段EVI最小值；EVI(EOS)表示提取EOS的EVI阈值，EVI_MIN2表示曲线下降阶段EVI最小值。参数α和β为动态阈值，本文根据前人研究和实际情况[34]，将SOS的阈值α设置为0.2，EOS的阈值β设置为0.4。提取2002~2021年间北京市的SOS与EOS，在提取的SOS和EOS数据中，分别随机抽取了30个样本点，并在MOD12Q2物候产品中空间位置相对应的植被物候参数进行对比验证(图3)。SOS_MCD12Q2和SOS_MOD13Q1之间的相关性显著(R² = 0.852, p < 0.01, RMSE = 4.917 d)，EOS_MCD12Q2和EOS_MOD13Q1之间的相关性显著(R² = 0.747, p < 0.01, RMSE = 5.074 d)。该结果说明，本文提取的SOS与EOS与MODIS官方产品具有较高一致性，证明了提取结果的准确性。

Figure 3. Phenology parameter validation

图3. 物候参数验证

3.4. Theil-Sen斜率与Mann-Kendall检验

本研究采用Theil-Sen Median趋势分析法对北京市植被物候的变化趋势进行计算，分析其时空演变格局。Theil-Sen斜率(T-S)是一种基于成对点的所有线的斜率的中值的非参数趋势计算方法，其计算公式为：

$\beta =\text{Median}-\frac{\left(x_j-x_i\right)}{j-i},\forall j>i>1$ (5)

式中：β为斜率；Median为中值，x_j和x_i分别对应物候序列第i年和第j年的值。β的正负表示趋势方向；当β > 0时，物候呈推迟或者延长趋势；当β < 0时，则反映为提前或者缩短趋势。

采用Mann-Kendall趋势检验对植被物候年际变化进行显著性测试。Mann-Kendall显著性检验是一种非参数检验方法，样本无正态分布要求，并且对离群点不敏感，适用于长时间序列分析。

$S={\displaystyle \sum _{i=1}^{n-1}{\displaystyle {\sum }_{j=i+1}^{n}f\left(x_j-x_i\right)}}$ (6)

$f\left(x_j-x_i\right)=\{\begin{array}{l}+1,x_j-x_i>0\hfill \\ 0,x_j-x_i=0\hfill \\ -1,x_j-x_i<0\hfill \end{array}$ (7)

$Var\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{18}$ (8)

$Z=\{\begin{array}{l}\frac{S-1}{\sqrt{Var\left(S\right)}},S>0\hfill \\ 0,\left(S=0\right)\hfill \\ \frac{S+1}{\sqrt{Var\left(S\right)}},S<0\hfill \end{array}$ (9)

式中：x_j和x_i分别表示物候序列中第j年和第i年的值，n为年份总数，S为检验统计量，Z为标准化检验统计量。显著性水平区p < 0.05，当$\left|Z\right|$ > 1.96时，则趋势通过显著性检验。确定植被物候趋势显著性的标准见表2。

Table 2. Phenological trend significance category determination

表2. 植被物候变化趋势显著性类别判断

变化趋势	β值	Z值
极显著减少	β < 0	$\left|Z\right|$ > 2.58
显著减少	β < 0	2.58 ≥ $\left|Z\right|$ ≥ 1.96
不显著减少	β < 0	1.96 ≥ $\left|Z\right|$ ≥ 0
不显著增加	β > 0	1.96 ≥ $\left|Z\right|$ ≥ 0
显著增加	β > 0	2.58 ≥ $\left|Z\right|$ ≥ 1.96
极显著增加	β > 0	$\left|Z\right|$ > 2.58

3.5. 植被物候稳定性分析

以变异系数(CV)衡量物候变化的稳定性。CV表征长时序数据的相对离散度；数值越小，序列越稳定；数值越大，说明波动越明显。公式如下：

$\text{CV}=\frac{\sqrt{\frac{1}{n-1}{\displaystyle {\sum }_{i=1}^n{\left(x_i-\overline{x}\right)}^2}}}{\overline{x}}$ (10)

式中：CV为变异系数，n为研究时段长度，i表示年份($i=1,2,\cdots ,n$ )，x_i是第i年的SOS (EOS或GSL)，$\overline{x}$ 为2002~2021年的SOS (EOS或GSL)平均值。

3.6. 城市分区及城市化强度梯度的设立

本研究依据Xiong等[32]编制的“城–郊–村”范围数据集(USR)对北京市进行空间分区，形成2002~2021年每年的城市、郊区、农村三类类空间分区。进行分区的目的是：减小城市化强度法中城乡判别的偏差；凸显城郊村不同分区中植被物候对城市化强度的敏感性差异，更真实地反映物候对城市化的响应。在每一类空间分区中，城市化强度定义为研究区内250 m × 250 m网格内的不透水面百分比(ISP) [28]。在每一个空间分区内部，划分城市化强度梯度从0%到100%进行排序，本研究将北京市2002~2021年各年度的城市化强度按1个百分点的等宽间隔划分为100个等级(0%~1%、1%~2%、……、99%~100%)，以表征三类空间分区中各自的城市化强度梯度。

3.7. 城市化强度对植被物候的影响

为分析城市化对植被物候的响应，先在各分区内以ISP = 0%的网格为基准，计算其物候平均值；随后将其他强度等级与该基准相减，得到每1%间隔的ΔISP及对应的ΔSOS、ΔEOS、ΔGSL，并以线性回归评估相对ISP = 0%的物候影响。如下所示：

$\Delta {\text{SOS}}_{ti}=a_1\Delta {\text{ISP}}_{ti}+b$ (11)

$\Delta {\text{EOS}}_{ti}=a_2\Delta {\text{ISP}}_{ti}+c$ (12)

$\Delta {\text{GSL}}_{ti}=a_3\Delta {\text{ISP}}_{ti}+d$ (13)

其中ΔSOS_ti、ΔEOS_ti和ΔGSL_ti分别为第t年城市化强度为i时的ΔSOS、ΔEOS和ΔGSL；ΔISP_ti为第t年城市化强度为i时与初始基准网格(ISP = 0%)的城市化强度之差；i的值在1~100之间。参数a表示城市化强度对植被物候的影响程度，b，c，d分别表示当ISP = 0%时，SOS，EOS，GSL的基准值。

4. 结果

4.1. 植被物候趋势分析及时空变化特征

北京市2002~2021年均SOS呈现提前趋势，年均EOS呈现推迟趋势，年均GSL呈现延长趋势(图4)。北京市20年SOS平均为107.80 d，平均每年提前0.72 d∙a⁻¹，其中，2017年最早的SOS为96.58 d；平均EOS为297.37 d，平均每年延迟0.55 d∙a⁻¹，最晚的的EOS是2012年的303.77天；平均GSL为189.56 d，平均每年延长1.26 d∙a⁻¹，2020年最长GSL为202.12 d。

Figure 4. Interannual trend in vegetation phenology in Beijing, 2002~2021

图4. 2002~2021北京市植被物候年际变化趋势

2002~2021年间北京市植被物候空间变化显著，且城市附近的植被物候变化较大。结果显示，SOS以提前为主：88.46%的像元呈提前，仅6.93%表现为推迟，其中17.00%达到显著提前(p < 0.05) (图5(a)、图5(b))。EOS则整体推迟：82.20%的像元为推迟，11.28%为提前，其中10.49%显著推迟(p < 0.05) (图5(c)、图5(d))。GSL普遍延长：90.10%的像元延长，6.33%缩短，其中23.52%显著延长(p < 0.05) (图5(e)、图5(f))。总体来看，2002~2021年间北京市整体植被物候呈现出SOS提前、EOS推迟、GSL延长的显著趋势。

Figure 5. Spatial trends in vegetation phenology in Beijing and statistically significant areas

图5. 北京市植被物候空间变化趋势及其显著性区域

2002~2021年北京市SOS主要集中在95~110 d。城市地区SOS最早，主要集中在70 d (图6(a))。整体分布特征为由西北向东南提前。EOS主要集中在290~305 d之间。其中北京市城市周围EOS最晚，主要集中在320 d (图6(b))。西北部EOS出现的时间相对较早，整体分布特征为由西北向东南推迟。GSL主要集中在180~95 d之间，城市周围GSL最长，主要集中在270 d，整体分布特征为由西北向东南延长(图6(c))。

SOS的变异系数主要集中在0.08~0.12之间，呈现由北向南逐渐扩大的趋势(图6(d))。中部和西北部的农田地区CV变异系数较高，集中在0.16~0.2之间；EOS的变异系数主要集中在0.02~0.03之间(图6(e))，证明EOS较为稳定。中部的城市地区和农田以及西北部的农田CV变异系数略高，在0.04~0.05之间；GSL的变异系数在0.1~0.14之间(图6(f))。整体趋势为由北向南逐渐扩大，南部和西北部的农田地区CV变异系数较大，在0.14~0.15之间。

Figure 6. (a)~(c) Mean vegetation phenology in Beijing, 2002~2021; (d)~(f) Coefficient of variation of vegetation phenology in Beijing, 2002~2021

图6. (a)~(c)北京市2002~2021年植被物候均值；(d)~(f)北京市2002~2021年植被物候变异系数

4.2. 各分区植被物候的变化趋势

4.2.1. 各分区植被物候分布

各分区中，城市区域SOS最早，EOS最晚，GSL最长。城市SOS均值为97.58 d，而郊区和农村分别为106.31 d、108.81 d (图7(a))；城市EOS均值308.06 d，郊区和农村分别为301.00 d、298.10 d (图7(b))；城市GSL均值210.02 d，郊区和农村分别为194.21 d、188.82 d (图7(c))。综上所述，随着区域距离城市越远，SOS逐渐推迟、EOS逐渐提前、GSL逐渐缩短，表明城市化促使SOS提前、EOS推迟，GSL延长。

Figure 7. Distribution of vegetation phenology across Urban-Suburban-Rural (USR) Zones in Beijing

图7. 北京市城郊村分区植被物候分布

4.2.2. 各分区上植被物候的年际变化

2002~2021年，各分区的SOS均呈提前趋势：城区0.23 d∙a⁻¹、郊区0.70 d∙a⁻¹、农村0.82 d∙a⁻¹，农村区域的提前幅度最大(图8(a))；EOS在四区均呈推迟趋势：城区0.07 d∙a⁻¹、郊区0.38 d∙a⁻¹、农村0.60 d∙a⁻¹、农村推迟最显著(图8(b))；GSL均呈延长趋势：城区0.30 d∙a⁻¹、郊区1.07 d∙a⁻¹、农村1.41 d∙a⁻¹，其中农村区域延长最显著(图8(c))。总体来看，物候年际变化幅度呈梯度分布，农村变化幅度最大，城区变化最缓。

Figure 8. Interannual trends of vegetation phenology across Urban-Suburban-Rural (USR) zones in Beijing

图8. 北京市城郊村分区植被物候年际变化趋势

4.3. 各分区植被物候对城市化的响应

2002~2021年间，北京市各分区的SOS与ISP显著负相关，EOS和GSL与ISP显著正相关(图9)。表明随着城市化强度增强，各分区均呈现SOS提前、EOS推迟、GSL延长的一致性趋势，当ISP从10%提升至100%时，城市，郊区，农村的SOS分别提前11.27 d，11.05 d，6.61 d；EOS分别推迟2.56 d，3.35 d，2.81 d；GSL分别延长13.83 d，14.40 d，9.42 d。总体来看，郊区和郊区对ISP响应最敏感，并且物候变化幅度相似，农村区域最弱。

Figure 9. Spatial distribution of vegetation phenology across urbanization-intensity levels in Beijing’s Urban-Suburban-Rural (USR) zones

图9. 北京市城郊村分区不同城市化强度下的植被物候分布

各分区对比ISP = 0%的网格，城市和郊区区域中，ΔISP每增加10%，城市和郊区区域的SOS分别提前0.96 d，0.97 d (p < 0.01)，EOS推迟0.28 d，0.31 d (p < 0.01)，GSL延长1.24 d，1.30 d (p < 0.01) (图10(a)、图10(b))；农村区域中ΔISP每增加10%，SOS提前0.56 d，(p < 0.01)，EOS推迟0.26 d (p < 0.01)，GSL延长0.82 d，相同的ΔISP变化幅度下，农村区域的物候变化幅度远小于城市和郊区(图10(c))。

Figure 10. Effects of urbanization intensity on vegetation phenology across different subregions of Beijing: (a) urban area; (b) suburban area; (c) rural area

图10. 北京市不同分区的城市化强度对植被物候的影响。(a) 城市区域；(b) 郊区区域；(c) 农村区域

5. 讨论

2002~2021年间，北京市SOS平均每年提前约0.72 d、EOS每年推迟约0.55 d、GSL每年延长约1.26 d；这一物候趋势与国内外多数城市的物候研究相一致[9] [11]。热岛效应导致春季积温的升高促进植被休眠解除提前，造成SOS的提前，秋季冷量累积不足导致EOS的推迟，最终导致GSL增加[9]。东南部的农田区域和城市周边变化更显著，并且这些地区的SOS和GSL波动相对更大(CV 0.10~0.14)，频繁的人为管理是物候波的重要原因。植被物候呈现出东南–西北方向的空间梯度也与北京“西北山地、东南平原”的温度–海拔梯度一致[26] [35]。

从各分区植被物候来看，城市区的SOS最早(97.60 d)、EOS最晚(308.10 d)、GSL最长(210.00 d)；城市，郊区，农村分区中，表现为SOS逐渐推迟、EOS逐渐提前、GSL逐渐缩短的梯度。与缓冲区的物候特征结论一致[11] [26]，证明分区的准确性。城郊村的物候差异主要是由城市热岛和人为管理造成春季物候提前和秋季物候推迟[3] [5]。年际变化幅度呈现城市核心区最小、其他分区更大的现象，这与一些城市的“核心区变化趋于平稳”相一致，这是因为植被对城市中心热环境产生了适应性，所以年际变化幅度较小[24] [36]；其次，城市区域树种配置更趋稳定且受管理约束，物候的变化趋势变小[15]；相对的，农村地区受年际间气候变化影响更大，因而物候变化幅度更大[33]。

沿不透水面比例(ISP)从低到高，各分区均出现“随着城市化强度的提升，SOS提前、EOS推迟、GSL延长”的一致响应，这与华北地区的研究结果一致[28]。但敏感度则呈现不同分区特征：在城市与郊区，ISP占比较大，热岛效应更强，因此城市和郊区的SOS和GSL对ISP的变化更加敏感[12] [24]；农村虽有农田灌溉，但管理节律受农事影响，因此物候对ISP的敏感度低于城市[36]。相对各分区ISP = 0%的地区，ΔISP对植被物候具有显著且同向的驱动效应，物候对城市化强度的敏感性为：郊区 > 城市 > 农村。可能在于城市与郊区更强的热岛效应响应[5] [7]；城市核心区因地类混合与热岛影响，且处于相对稳定的生长环境中，致使物候随城市化强度的变化幅度小于郊区[24] [33]。农村区因ISP较低、像元由森林、草地、农田主导，且夜间灯光等人为扰动对物候的影响相较温度较弱，因此对ISP变化表现出不敏感[3] [36]。值得注意的是，三类分区中的EOS对ΔISP敏感性较低，这主要是因为秋季EOS主要受温度和风速控制[35]，因此各分区植被物候对ΔISP的响应不敏感。

除热岛效应外，城市化还伴随一系列环境因子的叠加影响。例如，氮沉降在城市及郊区更加显著，可通过改变土壤养分状况与叶绿素含量促进SOS提前[36]；城市和郊区的大气CO₂浓度升高增强了光合速率与碳同化能力，从而延长GSL [24]。此外，光污染、空气污染物对辐射吸收与散射的调节作用也可能间接影响物候[16]。这些因素的交互效应值得未来研究中进一步量化与模型化。

本文在各分区内用0%~100%连续梯度量化城市化强度，既减少城乡区域的错分的情况，又能够区分不同分区对城市化强度的不同响应，因此该结合策略以减少城市化强度法在植被物候研究中的不确定性。

本研究综合使用了空间分辨率分别为30 m、250 m及1 km的多源遥感数据。尽管本文通过统一重采样与空间对齐减少误差，但不同分辨率数据在研究过程中仍可能会造成一些偏差。例如可能影响ISP与物候参数之间的拟合强度(R²)及显著性水平。因此，未来可考虑采用分辨率匹配或多尺度敏感性分析以进一步量化其不确定性。此外，部分分区中城市化强度与物候指标的回归结果(图7、图9)中R²值相对较低，说明年际变化中存在较强的非线性与多因子干扰。城市植被物候不仅受城市化强度影响，还受到气候波动、土地管理和物种组成等多重因素的交互作用。线性趋势模型虽能揭示总体方向，但难以刻画复杂的非线性响应特征[12]，这也是城市生态系统高复杂性的体现。未来研究可结合非线性拟合或机器学习方法(如GAM或RF模型)以更准确地解释城市化与物候之间的多维关系。

6. 结论

本研究基于2002~2021年北京市EVI数据和全球不透水面数据，量化北京市植被物候时空演变特征，利用城市化强度法和城郊村范围数据集，揭示植被物候对城市化的响应。结论如下：

2002~2021年间，北京市植被物候表现为“SOS提前、EOS推迟、GSL延长”，东南平原的农田地区与城市周边变化更为显著，西北山区相对较弱，表明城市化对植被物候产生影响。

北京市的物候在城市区域SOS最早(97.60 d)、EOS最晚(308.10 d)、GSL最长(210.00 d)；城市–郊区–农村之间依次表现为SOS推迟、EOS提前、GSL缩短。同时，北京市物候的年际变化幅度出现“城市区域变化幅度较小(0.30 d∙a⁻¹)，郊区(1.07 d∙a⁻¹)，农村区域变化幅度更大(1.41 d∙a⁻¹)”的特征。

随着不透水面比例(ISP)升高，城市–郊区–农村各分区普遍表现为SOS提前、EOS推迟、GSL延长，但敏感度显著分区化：ΔISP每增加10%，城市与郊区GSL延长1.24 d，1.30 d，农村GSL延长0.82 d。