徐 勇, 赵 纯, 窦世卿, 郝文强, 郑志威, 靖娟利

(桂林理工大学 测绘地理信息学院, 广西 桂林 541006)

1 研究区概况

环渤海地区以“C”字形环抱渤海，包括三省两市，即辽宁省、山东省、河北省、天津市以及北京市，位于34°22′—43°28′N，113°26′—125°52′E。气候类型为温带季风气候，主要植被带为暖温带落叶阔叶林带。地形主要以平原、丘陵和山地为主，整体上地势呈东南低西北高的分布格局。环渤海地区是中国最为发达的经济区之一，基础设施完善，自然资源丰沛，工业基础扎实，科技实力雄厚。2019年环渤海地区年末总人口约2.6亿(2020年统计年鉴)，约占全国总人口的18.38%，为中国人口最为稠密的地区之一，尤其是北京市和天津市，不仅是中国政治、经济和文化中心，2020年人口密度分别高达1 334和1 158人/km2，分别居全国第2和第3位，而山东省和河北省为中国人口大省，人口总量分居全国第2和第6(第7次人口普查数据)。

2 数据来源与研究方法

2.1 数据来源

本文所使用 MODIS NDVI数据来源于NASA MODIS陆地产品组发布的植被指数产品MOD13A3，该产品时间分辨率为每月，空间分辨率为1 km，时间跨度从2000年2月至2020年12月，该数据已经经过大气校正，去除了水、云、重气溶胶和云影的影响。本文通过对得到的每期MOD13A3数据集进行影像镶嵌、重采样、投影转换、数据格式转换、最大值合成等预处理，生成2000—2020年覆盖中国全境的年最大值植被NDVI时间序列，然后通过裁剪得到环渤海地区植被年最大合成NDVI时间序列。2000年1月植被NDVI数据缺失。由于1月为冬季，植被NDVI较低，故本文使用最大值合成法得到年最大合成植被NDVI，缺失月份不会对本文研究结果产生影响。人口密度数据来源于WorldPop发布的Population Density数据，时间分辨率为每年，空间分辨率为1 km，时间跨度为2000—2020年。人口密度数据需要经过镶嵌、投影转换、图像裁剪和重采样等预处理，裁剪得到2000—2020年环渤海地区人口密度时间序列。

2.2 研究方法

2.2.1 Theil-Sen Median趋势分析与Mann-Kendall显著性检验 本文采用Theil-Sen Median趋势分析和Mann-Kendall显著性检验[15]，分析环渤海地区2000—2010年、2010—2020年和2000—2020年植被NDVI时空演变特征及其变化的显著性。Theil-Sen Median趋势分析是一种非参数统计的趋势分析方法，可以有效去除离散值和异常值对趋势结果的影响，常被用于探究长时间序列数据的变化趋势。

Theil-Sen median趋势分析计算公式为：

(1)

式中：NDVIi，NDVIj分别表示第i和第j年的植被NDVI值；βNDVI代表植被NDVI的变化斜率,当βNDVI>0时,表示在研究时间段内植被NDVI呈上升趋势,植被覆盖有所增加; 当βNDVI=0时,表示在研究时间段内植被NDVI基本保持不变,植被覆盖较为稳定;当βNDVI<0时,表示在研究时间段内植被NDVI呈下降趋势,植被覆盖有所减少。

Mann-Kendall显著性检验是稳健的非参数统计方法，不需要方差正态性假设，在一定程度上避免异常值对分析结果的影响。本文利用Mann-Kendall显著性检验来判断植被NDVI时间序列变化趋势的显著性。具体公式为：

给定NDVI时间序列：

{NDVIi},i=2000，2001，2002，…,2020

(2)

计算检验统计量S：

(3)

符号函数sgn：

(4)

计算方差：

(5)

定义标准化统计量Z：

(6)

按照Mann-Kendall显著性检验的检验统计量Z值，将植被NDVI变化趋势分成4个等级：显著下降、不显著下降、不显著上升和显著上升。当|z|>1.96即置信水平为95%[16]下，为显著变化，当0<|z|<1.96，为不显著变化。

(7)

(8)

2.2.3 格兰杰因果分析 为检验X组时间序列是否为Y组时间序列的原因，采用格兰杰因果关系检验法验证其因果关系[18]。本文设Y为年际植被NDVI即Y={yt}，X为年际人口密度即X={xt}，假设两者之间是否具有因果解释能力，并设定滞后阶，进行格兰杰检验分析人口密度能否影响植被NDVI。其回归模型为：

(9)

式中：βt为时间序列X的趋势项；γt为时间序列Y的趋势项;n为滞后阶数;α为常数项;εt为残差项。

假定原假设H0：人口密度不是引起植被NDVI的格兰杰原因(H0:∀t∈(1,2,…,n),γt=0)；其余假设H1：人口密度是引起植被NDVI的格兰杰原因〔H1:∃t∈(1,2,…,n),γt≠0〕。当不考虑X的情况下(即假设βt=0)，设置回归模型的残差平方和RSR;当考虑X即人口密度的情况下(即假设βt≠0)，设置残差平方和RSU。构建F统计值：

(10)

选取置信水平为95%，若F≥F0.05(n,N-2n-1)，即拒绝原假设H0，认为人口密度是引起植被NDVI的格兰杰原因，反之则保持原假设。

式中：F为统计值； RSR为没有X约束下的残差平方和； RSU为在X约束下的残差平方和；n为待估参数的个数；N为样本容量。

2.2.4 重心转移模型和标准差椭圆 为分析2000—2020年环渤海地区植被NDVI和人口密度迁移特征，本次研究运用属性重心迁移模型[19]进行植被NDVI和人口密度演化模式的定量化表达，具体公式为:

(11)

(12)

标准差椭圆从数据范围、数据集聚性、变化方向和形状等多个角度，呈现地理要素的空间格局分布特征。方向角是正北方向与顺时针旋转的椭圆长轴形成的夹角，代表数据空间分布的主趋势方向，长半轴表示数据在主趋势方向上集聚性，短半轴反映数据在此趋势方向上集聚性。主要计算公式为：

tanθ=

(13)

(14)

(15)

3 结果与分析

3.1 植被NDVI的时空演变特征

3.1.1 NDVI的时间变化趋势 由图1可知，整体上，2000—2020年环渤海地区植被NDVI呈波动上升趋势，上升速率为0.022/10 a。2000—2020年环渤海地区植被NDVI均值为0.728，最小值出现在2000年，为0.666，最大值出现在2011年，为0.748，且2009和2014年植被NDVI值明显低于相邻年份。2000—2010年环渤海地区植被NDVI呈上升趋势，上升速率为0.057/10 a， 2000—2010年环渤海地区植被NDVI均值为0.716，略低于2000—2020年环渤海地区植被NDVI平均值。相较于2000—2010年和2010—2020年环渤海地区植被NDVI上升趋势较为平缓，上升速率仅为0.000 1/10 a。2010—2020年环渤海地区平均植被NDVI值最高，为0.738。综上所述，整体上，2000—2020年环渤海地区植被覆盖呈改善态势，改善趋势明显，且2000—2010年植被NDVI改善程度强于2010—2020年，这与徐勇等[15]和张鹏骞等[20]的研究结果一致。自2000年以来，针对环渤海地区采取了林草地植被保护、小流域生态治理和高原地区水土保持等一系列林业生态工程，且取得了良好的效果，加之为筹备2008年奥运会，北京市开展了一系列绿化工程，短期内明显提升了区域植被NDVI值。2009年植被NDVI值明显低于2008和2010年，这与李卓等[21]的研究结果一致，是由于2009—2010年环渤海地区遭受了历史极严重的旱灾，导致整体植被NDVI值较低。

图1 环渤海地区2000—2020年植被NDVI变化趋势

3.1.2 植被NDVI的空间变化特征 本文利用Theil-Sen Median趋势分析法得到2000—2010年、2010—2020年以及2000—2020年环渤海地区植被NDVI变化趋势，并根据Mann-Kendall显著性检验，检验其变化趋势的显著性。由表1和图2可知，2000—2010年和2010—2020年环渤海地区植被NDVI变化呈明显空间异质性，变化斜率分别在-0.076～0.072/a和-0.087～0.070/a之间。2000—2010年植被NDVI呈上升趋势的区域面积占83.83%，其中，呈显著上升的面积占30.26%，均匀分布在研究区内，无明显集现象；呈显著下降的面积仅占1.35%，零星分布在各大中城市中心及其周边地区。而2010—2020年植被NDVI呈上升趋势的面积为50%，其中，呈显著上升的面积仅占9.85%，主要分布在北京和河北省西北部以及辽宁省东部；呈显著下降的面积占5.71%，零星分布在辽宁省中部、天津以及河北和山东部分地区。由上可知，2000—2010年环渤海地区植被NDVI以上升为主，上升区域面积明显高于下降区域面积，而2010—2020年环渤海地区植被NDVI上升趋势明显减缓，植被NDVI呈上升趋势面积等于呈下降趋势面积。整体上，2000—2020年植被NDVI变化斜率在-0.034～0.037/a之间，植被NDVI整体呈上升趋势。植被NDVI上升区域面积远大于下降区域面积。植被NDVI呈上升趋势的面积占76.37%，下降趋势面积占23.68%。植被NDVI显著上升区域面积占总面积的45.18%，主要分布在辽宁省东部和西部、北京和河北西北部、山东省中部以及河北省和山东省交界处。这与Zheng等[5]和Jiang等[22]和研究结果一致，以上地区植被NDVI上升，主要得益于环渤海地区多年来实施退耕还林，京津风沙源治理工程和三北防护林工程等项目，使得地表植被覆盖程度得到提升，植被生态稳定性以及抵御风险的能力得到加强[6,8-9]。植被NDVI呈显著下降趋势区域占总面积的6.42%，主要分布在辽宁省中部、北京市城市中心和南部、天津市南部、河北省西南部至东北部沿海地区、山东省少部分地区，呈斑块状零星分布，这与Zheng等[5]、Jiang等[22]、孙涛等[23]和赵安周等[24]的研究结果一致，由于北京、天津、石家庄等城市周边的工业用地急剧扩张，耕地、有林地等生态空间大幅度缩减，人类活动对植被NDVI的变化产生负向影响[11,25-26]。

表1 环渤海地区植被NDVI变化趋势面积比例统计

3.2 植被NDVI与人口密度变化的相关性

图2 环渤海地区2000—2020年植被NDVI变化趋势及显著性检验

表2 环渤海地区人口密度冷热点区域面积比例统计 %

3.2.2 植被NDVI与人口密度变化的相关性 为检验环渤海地区人口密度和植被NDVI的因果关系，本文利用格兰杰因果关系检验法对两者的相关关系进行验证。该方法可判断两个时间序列之间是否存在因果性[30]。采用E-G两步法[31]对人口密度和植被NDVI两者的时间序列进行平稳性验证，由协整检验一般规律可得，两者的时间序列较平稳，存在协整关系进行格兰杰因果检验。滞后期取1～6，选择最优解作为结果进行分析。当p值<0.05时，拒绝原假设，即认定两者之间存在单向格兰杰原因[32]。根据表3可知，滞后期为2 a时，环渤海地区人口密度是植被NDVI的格兰杰原因，但植被NDVI不是人口密度的格兰杰原因，山东省人口密度是植被NDVI的格兰杰原因，但植被NDVI不是人口密度的格兰杰原因。滞后期为1 a时，北京市、河北省和辽宁省的人口密度是植被NDVI的格兰杰原因，但植被NDVI不是人口密度的格兰杰原因。滞后期为4 a时，天津市人口密度是植被NDVI的格兰杰原因，但植被NDVI不是人口密度的格兰杰原因。综上所述，环渤海地区人口密度对植被NDVI存在单向格兰杰原因，人口密度的改变是引起植被NDVI变化的格兰杰原因，而植被NDVI的变化不是造成人口密度改变的格兰杰原因。为探究植被NDVI变化与人口密度变化的空间耦合关系，本文将植被NDVI变化趋势结果与人口密度冷热点区进行叠置分析，得到植被NDVI与人口密度的空间耦合特征和数量关系，结果如图4所示。2000—2010年各冷热点区植被NDVI呈上升趋势的面积远大于呈下降趋势的面积。植被NDVI呈显著上升的面积占比随人口密度变化热点区到冷点区依次上升，其在冷点区面积占比最大，为37.11%，在热点区面积占比最小，为13.83%。而植被NDVI呈显著下降的面积占比在各冷点区面积占比与显著上升面积占比成相反，显著下降的面积占比在热点区最大，为31.16%，在冷点区最小，为9.17%。由些可见，2000—2010年环渤海地区的植被覆盖整体呈改善趋势，其中冷点区改善程度最为明显，次冷点区其次。

图3 环渤海地区人口密度冷热点区域分布

由图4可见，2010—2020年植被NDVI呈显著上升的面积占比随人口集聚程度的下降呈现先减少后升高的趋势，其中次热点区面积占比最小，为3.97%，冷点区面积占比最大，为14.02%。同样，植被NDVI呈显著下降的面积占比呈现先减少后升高的趋势，其中次冷点区面积占比最小，为4.81%，热点区面积占比最大，为7.02%。相较于2000—2010年，2010—2020年植被NDVI变化趋势的面积占比变化规律不显著，且植被NDVI显著上升的面积占比有所减少，植被NDVI显著下降的面积占比增加。2000—2020年冷点区、次冷点区和次热点区的植被NDVI呈上升趋势变化的面积大于呈下降趋势面积。植被NDVI呈显著上升区域面积占比排序为：冷点区>次冷点区>次热点区>热点区，其中在冷点区的面积占比为65.14%，热点区的面积占比为19.94%。由此可知，植被NDVI呈显著上升的面积占比随人口集聚程度的下降呈上升趋势。植被NDVI显著下降变化趋势的面积占比排序为：热点区>次热点区>次冷点区>冷点区，其中热点区的面积占比为18.50%，冷点区的面积占比为1.88%，植被NDVI呈显著下降的面积占比随人口集聚程度的下降呈现降低的趋势。2000—2020年环渤海地区植被NDVI在次热点区、次冷点区和冷点区呈改善趋势，其中冷点区改善程度最为明显，但热点区的植被NDVI呈退化趋势。综上所述，环渤海地区人口密度变化是引起植被NDVI改变的格兰杰原因，且对区域植被生长具有双重作用。环渤海地区植被NDVI呈下降趋势的面积主要分布在热点区和次热点区，原因是热点区和次热点区主要位于各大中型城市及其周边地区，以上地区人口基数大，人口增长快，其中北京和天津的人口增长最为显著，其他省会城市次之。另外，以上地区城市化程度高，如2019年北京、天津、辽宁和山东的城市化率在分别为86.6%，83.5%，68.1%和61.5%，高于全国平均水平，相较于2000年，城市化率分别增长了9.1%，12.3%，13.9%和23.5%(中国统计年鉴2001—2020)；其次，热点区和次热点区城市的第二、三产业对劳动力需求较大，对流动人口具有较强的吸引力，导致热点区和次热点区城市人口集聚程度较高，城市扩张，人类活动较为频繁，对植被覆盖造成破坏，从而抑制植被生长。环渤海地区植被NDVI呈上升趋势的部分主要集中在冷点区，冷点区的地理位置缺乏优越区位条件，基础设施较为落后，导致人口向北京市、天津市两个直辖市以及石家庄、邯郸、济南、青岛和沈阳等经济核心城市流动。冷点区的经济开发活动较低，对自然环境的负向干扰程度较弱，为了保障环渤海地区经济的可持续发展，冷点区主要推行“三北”防护林工程、京津风沙源治理工程以及退耕还林还草工程等[6,8-9]，对植被生长起促进作用。

表3 人口密度和植被NDVI两个时间序列的Granger检验结果

图4 环渤海地区人口密度冷热点区域植被NDVI变化显著性面积比例统计

3.3 植被NDVI与人口密度重心迁移

3.3.1 植被NDVI重心迁移 由图5可知，2000—2020年环渤海地区植被NDVI标准差椭圆的重心经度在118°39′54″—118°46′13″之间，重心纬度在39°03′40″—39°10′52″之间，表明环渤海地区植被NDVI呈“东北—西南”走向，且在“东北—西南”空间格局进一步扩大。长半轴标准差由4.890 km下降至4.846 km，而短半轴标准差由2.628 km上升至2.668 km，标准差椭圆变化较小，表明植被NDVI年际变化较为平稳。总体上，2000—2020年环渤海地区植被NDVI重心往西北方向迁移。研究区内植被NDVI的改变以及重心迁移受人类活动因子的影响[33]。主要有两个方面原因：一方面，与环渤海地区，尤其是沿海地区，经济发展和城市扩张紧密相关。随着改革开放，以京津冀为经济中心带动辐射周围城市经济发展，城市用地的急剧扩张，人口密度显著增加，导致耕地、草地和林地等转为城市建设用地，人类活动频繁，环渤海地区沿海地区的植被覆盖率下降显著，植被NDVI重心向西北方向迁移。另一方面，三北防护林工程和环渤海地区农业生态补偿建设，实施封山育林、退耕还林等山区保护措施，退耕还林等举措初见成效，增加了环渤海地区北部的森林覆盖率，提高水土涵养能力，尤其是环渤海地区东北和西北部地区，植被覆盖得到明显的改善[34-35]，因此，环渤海地区植被NDVI重心向西北方向迁移。

3.3.2 人口密度重心迁移 2000—2020年环渤海地区人口密度标准差椭圆形成东北—西南的空间分布格局，标准差椭圆参数变化较为稳定，重心经度由118°05′38″逐步下降至117°58′30″，而重心纬度由38°12′00″逐步上升至38°17′21″，方位角由59.48°逐步扩大至61.32°。长半轴标准差从4.254 km缩小至4.061 km，表明环渤海地区人口密度在“东北—西南”方向上呈收敛趋势，空间集聚性轻微上升，而短半轴标准差由2.405 km扩大至2.419 km，表明2000—2020年的环渤海地区人口密度在“西北—东南”方向上呈扩张趋势，空间集聚性轻微下降。由图6可知，2000—2020年环渤海地区人口密度变化较为平稳，重心逐步向西北方向迁移。影响环渤海地区人口密度重心迁移的主要因素是本地人口基数大以及外部人口的迁入[36]，使得人口密度逐年上升。辽中南经济集聚区、京津唐工业区的城市经济发展良好，城市化和工业化水平高，带动了产业经济由东南沿海逐步向西北内陆扩张，而产业发展和城市建设需更多人口弥补劳动力不足，吸引外来人口向该地区迁移，导致环渤海地区的人口重心向西北内陆迁移。2000—2020年的环渤海地区植被NDVI重心与人口密度重心均往西北方向迁移，表明人口重心迁移方向和NDVI重心迁移方向呈现一定相关性，人口密度重心迁移轨迹比较均一，单方向往西北迁移，植被NDVI重心整体往西北迁移，但迁移轨迹较为复杂，无明显迁移模型，振荡幅度较大，不与人口密度重心迁移轨迹呈明显的线性关系。原因是除人口密度外，植被生长还与区域地形、气候条件以及林业生态工程等息息相关，已有研究[11]表明，环渤海地区植被与降水、气温以及相对湿度呈明显的正向效应，近年来，环渤海地区年际气温和降水呈波动上升趋势，温暖湿润的环境对区域植被生长具有促进作用。此外，环渤海地区积极开展退耕还林还草工程、京津风沙源治理工程和三北防护林工程，以上工程不仅抵消人类活动对植被生长带来的负向影响，还提高了环渤海地区植被覆盖程度，增加了区域植被NDVI[37]。综上所述，环渤海地区人口密度与植被NDVI重心变化具有一致性，但人口密度与植被NDVI变化没有明显的线性相关性，人口密度是影响环渤海地区植被NDVI的因素，但该地区的植被NDVI同时还受气候变化和生态林业工程等因素的影响[6,8-9,11]。

图5 环渤海地区2000—2020年植被NDVI重心转移空间变化

图6 环渤海地区2000—2020年人口密度重心转移空间变化

4 结 论

(1) 2000—2020年环渤海地区植被NDVI整体呈上升趋势，上升斜率为0.022/10 a。2000—2010年植被NDVI上升速率高于2010—2020年植被NDVI上升速率。2000—2020年植被NDVI呈上升趋势的面积大于呈下降趋势的面积，其中显著上升区域面积占总面积的45.18%，主要分布于研究区的西北部和东北部。

(2) 2000—2020年环渤海地区人口密度变化以冷点区为主，主要从辽宁省东北部延伸至河北省西北部，以及山东省东南地区。相较于2000—2010年，2010—2020年人口密度变化的冷点区、次冷点区以及热点区均有所减少，主要表现为次热点区增加。

(3) 环渤海地区人口密度变化是引起植被NDVI变化的格兰杰原因，而植被NDVI的变化不是造成人口密度改变的格兰杰原因。且人口密度变化对植被NDVI变化的负向效应要强于正向效应，植被NDVI呈显著下降趋势面积主要集中于热点区，呈显著上升面积主要位于冷点区。

(4) 2000—2020年的环渤海地区植被NDVI重心与人口密度重心均往西北方向迁移，表明人口重心迁移方向和NDVI重心迁移方向呈现一定相关性，人口密度重心迁移轨迹比较均一，但植被NDVI重心迁移轨迹较为复杂，不与人口密度重心迁移轨迹呈明显的线性关系。

本文的研究结果基本符合一般规律，表明利用人口密度分析研究对植被NDVI的研究方法切实可行，但仍存在不足之处。首先考虑的人为因素较少，人口密度作为单一因子进行分析导致的结论可能过于片面，后期的研究可以加入地区GDP、人口劳动力等人为因子[38]，从多角度分析人类活动对植被NDVI的影响。其次未采用最优滞后期模型计算人口密度与植被NDVI的最优滞后阶，导致格兰杰因果关系方法的检验结果比较粗略，在一定程度上影响评价结果的精度。因此，后续的研究应从适当增加研究因子、采用先进数学模型计算相关性等方面对人类活动与植被NDVI的相关性进行深入研究，使得研究结果与结论更加准确与科学。