翟浩然，么嘉棋，王光辉，唐新明
(1.自然资源部国土卫星遥感应用中心，北京 100048；
2.首都师范大学 资源环境与旅游学院，北京 100048；
3.天津师范大学 京津冀生态文明发展研究院，天津 300380)
颗粒物(particulate matter，PM)是大气中存在的各种固态或液态颗粒状物质的总称，是大气的重要组成部分，各种颗粒状物质均匀地分散在大气中构成一个相对稳定的庞大的悬浮体系，即气溶胶体系[1]。近年来，在京津冀区域，快速的城市化进程使人类活动对颗粒物排放的贡献程度不断增加，PM2.5、PM10等颗粒物浓度严重超标，远远超过了世界卫生组织提供的指导值[2-3]。实现快速降低、削减大气颗粒物浓度，对城市规划、可持续健康发展具有重要意义。
植被是陆地生态系统的主体，在调节大气、减缓温室气体浓度的上升、维持整个生态系统的各个方面都发挥着重要作用。植被可以有效降低颗粒物浓度，这主要是由于以下两方面原因：一是植被的叶片、茎秆等表面分布着气孔，可以直接捕获空气中的颗粒物，起到滞尘作用；
二是植被使地面扬尘大大减少，并且由于地表粗糙度的增加、风速的降低，进一步促进颗粒物沉降[4-6]。韦晶等[7]发现林地与无污染和轻微污染天气均表现为正相关，而与高污染天气表现为负相关。陈文波等[8]发现南昌市中心城区植被覆盖度与PM2.5浓度显著负相关，但在不同季节、不同空间尺度上，二者的相关性均存在差异。莫莉等[9]发现北京市PM2.5、PM10等颗粒物浓度与林木覆盖率发展趋势相反，呈显著负相关。Escobedo等[10]发现智利圣地亚哥城市森林一年去除污染物约3 500 t，约占总量的60%，其中PM10的去除率最大。植被对降低颗粒物浓度、净化大气环境发挥着至关重要的作用，故此得到越来越多的关注。通过科学实验探究植被覆盖与颗粒物污染之间的关系，具有重要的现实意义。
本文基于卫星遥感反演的植被指数数据以及地面站空气质量监测数据，分析京津冀区域植被覆盖的时空变化规律，并从时间和空间两个角度分别探究植被覆盖对大气颗粒物浓度的影响。
1.1 研究区概况
京津冀地区位于113°27′E～119°51′E、36°05′N～42°40′N之间，由西北向东南可大致划分为坝上高原、燕山—太行山山脉以及海河平原3个部分，属温带半湿润大陆性季风气候，冬夏分明、晨昏各异。2020年，京津冀地区国土面积约占全国的2.3%，而地区生产总值和总人口则分别约占全国的8.5%和7.8%。作为中国经济最具活力、开放程度最高、创新能力最强、吸纳外来人口最多的地区之一，在经济快速发展的同时，京津冀地区空气污染现象频发，严重威胁着生态环境及人类健康。
1.2 数据
1) NDVI数据。归一化植被指数(normalized difference vegetation index，NDVI)是应用最广泛的植被指数之一，可以象征植被的覆盖程度。实验中用到的NDVI数据为MODIS传感器3级产品MOD13A3的C6版本，时间分辨率为1个月，空间分辨率为1 km[11]。实验共处理数据864景，处理后数据的有效取值范围为[-0.2，1]。为清晰直观地描述NDVI的空间分布和变化趋势，按数值大小对NDVI值做以下等级划分：[-0.2，0)为1级、[0，0.1)为2级、[0.1，0.2)为3级、[0.2，0.3)为4级、[0.3，0.4)为5级、[0.4，0.5)为6级、[0.5，0.6)为7级、[0.6，0.7)为8级、[0.7，0.8)为9级、[0.8，0.9)为10级、[0.9，1]为11级，且规定称1～4级为低覆盖度，5～7级为中覆盖度，8～11级为高覆盖度。
2)空气质量数据。实验中用到的空气质量数据为中国环境监测总站通过全国城市空气质量实时发布平台(http：//106.37.208.233：20035/)发布的全国各监测点位逐小时空气质量自动监测结果。实验选取分布于京津冀地区及周边城市的215个站点，2014年12月1日0时至2018年12月31日23时的逐小时PM2.5、PM10浓度数据。首先对数据进行计算整合，对于某一时间范围存在少量缺失的数据则以该月份其他日期的该小时均值进行代替；
随后使用克里金插值方法，获取研究区范围内各时间单位颗粒物浓度的空间分布。
1.3 方法
1) 空间重心计算。重心的物理学概念是指物体各部分所受重力之合力的作用点。地理学空间重心一般借鉴物理学上的概念，是指区域内某地理要素在空间平面上力矩达到平衡的点，通过比较分析不同时间重心的空间位置迁移轨迹，可以判断该地理要素的分布和变化规律[12-14]。将各时间求取的重心结果按照时间顺序连接，得到的折线即为重心迁移轨迹，它是要素属性的时间变化在空间位置上的直观反映。
2) 曲线估计。当两个变量间存在显著相关性时，一般通过建立回归模型的方式描述两个变量间的关系。为描述植被覆盖情况对颗粒物浓度的影响，实验中借助IBM SPSS 25.0软件平台进行曲线估计。用于曲线估计的回归模型包括线性、对数、逆、二次、三次、幂、复合、S曲线、Logistic、增长和指数共11种，通过拟合优度R2的大小判断各模型拟合结果的优劣。
2.1 植被覆盖的时空特征
分别统计京津冀地区NDVI月均值(图1)。总体来看，研究区整体月均NDVI表现出显著的周期性规律。1、2、3月除冬小麦外，其他主要作物和树木均未进入生长期，低覆盖度区域覆盖了研究区90%以上的范围，各月NDVI月均值在0.2左右，其中2月为全年最低值0.194。3月后，光照逐渐增强，温度迅速回升，大量作物开始播种，树木发芽展叶，各类植物纷纷进入生长期，植被覆盖度逐渐升高，低覆盖度范围下降，高覆盖度面积快速上升，至8月许多植被达到全年最茂盛的状态，NDVI月均值达到全年最高值0.698，高覆盖度区域覆盖了研究区82%以上的范围。8月后，天气逐渐转凉，光照强度下降，树木叶片掉落，草木枯萎，除冬小麦外的多数作物收获后不再播种，植被覆盖度迅速下降，低值重新取代高值占据了研究区大部分区域，至12月NDVI月均值下降至0.238，低覆盖度区域覆盖范围超过80%。
图1 NDVI月均值的时间变化
NDVI月均值空间分布的变化规律与各区域植被类型的空间差异有显著关系，海河平原大部区域以冬小麦-夏玉米为主，廊坊东部、坝上地区多种植单季玉米，邢台、邯郸、衡水东部是棉花的主要产区，东部沿海的天津、唐山、秦皇岛则多以单季水稻为主，林果种植的核心区域面积较小，分散在石家庄、衡水、沧州等地，蔬菜则主要分布在各大城市的郊区[15-17]。分别统计京津冀地区NDVI各月均值的空间分布(图2)，并分别计算相邻两月间的NDVI变化(图3)。结果显示，1—2月，除海河平原南部处于越冬期的冬小麦、燕山山脉部分常绿针叶林外，其他区域大多缺少植被覆盖，区域整体NDVI月均值较低，绝大多数区域为3级或4级低覆盖度。3—4月，海河平原冬小麦进入返青、拔节期，NDVI值有一定增长，达到6～8级中、高覆盖度，同时树木进入展叶期，山区的NDVI值有小幅增加，达到5～7级中覆盖度。进入5月，树木叶片生长加速，东北部燕山和西南部太行山地区NDVI值显著升高，部分区域达到9～10级高覆盖度，同时海河平原上棉花、水稻、单季玉米等作物陆续开始播种，使这些作物产区的植被覆盖度有所升高。6月，冬小麦进入收获期，夏玉米开始播种，因此造成海河平原大部分区域NDVI值有明显回落；
同时棉花、水稻、玉米等开始迅速生长，这些作物产区NDVI值大幅增加，达到6～8级中、高覆盖度。7—8月，各类树木和作物生长愈发茂盛，各地区NDVI值均有大幅增长，大部分区域达到9～10级高覆盖度。9—10月，随着玉米、棉花、水稻等各类作物的陆续收获和草木的枯萎，西北部坝上草原地区和东南部海河平原地区NDVI值显著下降至5～6级中覆盖度，同时海河平原冬小麦开始播种。11—12月，树木进入落叶期，东北部燕山和西南部太行山地区NDVI值显著降低，冬小麦开始出苗和分蘖，海河平原部分区域NDVI值小幅回升，全区NDVI值多在4～5级中、低覆盖度。
图2 NDVI月均值的空间分布
图3 NDVI月均值的各月间变化
分别对各月NDVI月均值求取空间重心(图4)。结果表明，京津冀地区NDVI各月月均值重心均集中在北京、保定、廊坊3市的交界地区，且迁移轨迹表现出一定规律性。1月，重心位于中部；
2月，由于冬小麦种植区NDVI值略有下降，重心向北微移；
3—4月，随着东南部冬小麦进入返青、拔节期，重心显著向南迁移；
5月，北部林地内树木进入展叶阶段，使重心大幅北移；
6月，北部NDVI继续增加，同时南部冬小麦收割、夏玉米播种使NDVI降低，导致重心继续北移；
7月，海河平原各类作物生长，重心向南移动；
8月，由于东北部林地NDVI稳定在较高水平，其他区域各植物继续生长，重心向西南微移；
9月，各地区NDVI值降低，其中西北部的坝上高原NDVI降速稍快，使重心向东南微移；
10月，随着东南部大面积农田进入收获期使NDVI下降较快，重心向北微移；
11—12月，东南部的冬小麦出苗、分蘖,使NDVI略有回升，重心向南微移；
下一年1月，气温快速下降，各处NDVI值继续降低，而东北部林地的降速较慢，重心向北微移。
图4 NDVI月均值的重心迁移轨迹
2.2 植被覆盖对大气颗粒物浓度的影响
研究区的NDVI月均值均以年为单位存在显著周期性，且在一年中表现为单峰单谷规律，这与颗粒物浓度月均值的变化规律显著相关。夏季，植被覆盖度较高，树木、草地和作物通过干沉降等方式对颗粒物浓度起到较强的削弱和抑制作用；
冬季，由于植被覆盖度较低，使得地表颗粒物暴露于空气中，增加了颗粒物浓度。为准确量化植被覆盖与颗粒物浓度之间的相关性，分别计算2015—2018年研究区NDVI与PM2.5、PM10浓度月均值时间序列的皮尔逊相关系数。结果表明，研究区NDVI与PM2.5、PM10浓度月均值表现分别为负相关，相关系数分别为-0.683和-0.737，二者均在0.01水平上显著相关。分别使用11种曲线模型对NDVI值与颗粒物浓度进行回归分析。各拟合模型中，NDVI与PM2.5浓度的拟合中，最佳拟合模型为幂函数，拟合优度R2达到0.580；
NDVI与PM10浓度的拟合中，最佳拟合模型为指数函数，拟合优度R2达到0.601。结果表明,PM10与NDVI的相关程度略强于PM2.5。模型方程如式(1)、式(2)所示。
PM2.5=31.488·NDVI-0.606
(1)
PM10=183.358·e-1.301·NDVI
(2)
以NDVI月均值为横坐标，PM2.5、PM10浓度月均值为纵坐标，绘制散点及模型曲线图(图5)。
图5 PM2.5、PM10浓度与NDVI散点及模型曲线
植被覆盖对颗粒物浓度存在较强的影响，除对浓度时间变化的影响外，植被对于颗粒物浓度的影响存在明显的空间效应。已有研究表明，土地利用和土地覆被类型对于颗粒物浓度的影响范围一般在5 km以内[18-21]，因此从0.5 km到5.0 km设定10个不同的缓冲区等级，各等级缓冲区半径分别相差0.5 km，分别统计各空气质量监测站周边区域内的NDVI均值，并分别在各缓冲区尺度上对NDVI和颗粒物浓度年均值、季均值进行相关分析(图6)。结果表明，除春季NDVI季均值与PM2.5浓度在1 km以下相关性较低外，其他季节在不同缓冲区范围内均表现为负相关，结果均通过显著性检验，再次印证了植被指数对颗粒物浓度的削减和抑制作用。随着缓冲区半径的增加，NDVI值与颗粒物年均值和各季节均值的相关性均表现出了先升后降的趋势，且不同季节中相关性最强的空间尺度有所差异。其中NDVI值与PM2.5和PM10浓度年均值分别在3 km和2.5 km尺度上相关性最强，冬季均在2 km尺度上相关性最强，春季分别在3.5 km和2.5 km尺度上相关性最强，夏季分别在3.5 km和3 km尺度上相关性最强，秋季分别在3 km和2.5 km尺度上相关性最强。各季节相关性最强的空间尺度中，夏季最大、冬季最小，春秋两季介于冬夏两季之间，这是由于夏季植被覆盖度最大，茂盛的植被对于颗粒物浓度的影响作用较强、范围较大，冬季植被覆盖度很低，对颗粒物的影响作用较弱、范围较小。另外，各季节NDVI对PM2.5的最佳影响尺度均稍大于PM10，表明植被覆盖在各季节对PM2.5的作用范围均略大于PM10。
图6 各尺度NDVI与PM2.5、PM10浓度的相关性
2.3 典型站点的比较
为进一步探究植被覆盖对颗粒物浓度的影响，分别选择位于植被区和非植被区的典型空气质量监测站点，对不同站点的颗粒物浓度变化规律进行对比。分别选择了位于北京市的昌平镇和定陵两个站点。昌平镇站位于昌平区城区内，是北京市的23个城市环境评价点之一，站点周边土地利用类型以建设用地为主，植被覆盖程度相对较低；
定陵站位于明十三陵景区内，是北京市唯一一个城市清洁对照点，周边土地利用类型以林地为主，植被覆盖程度相对较高。昌平镇站和定陵站之间的直线距离约9 km，二者地理位置接近、气候条件相似，但植被覆盖程度存在显著差异，可以作为对照进行对比实验。
前文实验表明，NDVI与PM2.5和PM10约在2.5～3 km范围内的相关性最高，据此分别统计两个典型监测站点周边3 km范围内NDVI月均值(图7)。结果表明，2015—2018年定陵站NDVI均值为0.523，其中冬春夏秋四季的NDVI均值分别为0.296、0.470、0.733、0.591；
昌平镇站NDVI均值为0.324，冬春夏秋四季的NDVI均值分别为0.183、0.281、0.457、0.374；
定陵站各月的NDVI值均显著高于昌平镇站，NDVI值为昌平镇站的1.6倍左右。NDVI结果证明了所选取站点植被覆盖程度的显著差异，可以分别代表植被区和非植被区的典型站点。
图7 典型站点NDVI月均值的时间变化
分别统计昌平镇站和定陵站不同季节的颗粒物浓度时均值(图8)。总体来看，植被区PM2.5和PM10浓度均显著低于非植被区，表明植被覆盖对颗粒物浓度有显著的削减作用。各季节中，植被区的PM2.5和PM10浓度时均值一般表现为双峰双谷变化规律，峰值分别出现在上午10至11时和晚上20至21时，谷值分别出现在早晨4至6时和下午15至17时，双峰规律较明显。与植被区不同，非植被区各季节PM2.5和PM10浓度时均值多呈单峰单谷或不明显的双峰双谷变化规律，具体表现为凌晨到早晨，随着人类活动强度的逐渐降低，在清晨出现第一个谷值，这段时间非植被区的浓度变化规律和植被区基本一致；
之后，随着人类活动的增强，颗粒物浓度有所上升，但植被区浓度上升幅度较小，且在9至10时便达到第一个峰值后浓度不再大幅上升，这比非植被区第一个峰值的出现提前约1个小时；
之后，随着气温升高，边界层高度的上升，大气湍流活动增强，非植被区的浓度有较大下降，但这在植被区表现不明显，仅表现为小幅下降或保持平稳，因此从第一个波峰到第二个波谷的变化不明显；
之后，随着气温下降和晚高峰的到来，非植被区颗粒物浓度有所上升，但植被区的上升幅度较小且持续时间较短，一般在18至19时便达到第二个峰值，这比非植被区提前约2个小时；
随后，植被区浓度快速下降，而此时非植被区浓度仍有上升趋势，因此植被区、非植被区的浓度之差一般在此时段达到最大，这与夜晚植被区内的人类活动强度显著下降有关。总体来看，非植被区一日内浓度的上下波动幅度较大，人类活动是造成两次浓度升高的主要原因；
植被区一日内浓度的上下波动幅度显著小于非植被区，植被覆盖一方面有效抑制了一日内两次颗粒物浓度升高的幅度，另一方面显著缩短了颗粒物浓度上升的持续时间，结果证明了植被覆盖对颗粒物浓度的降低具有重要作用。
图8 典型站点大气颗粒物浓度时均值
本文以京津冀地区作为实验区，以多种时间、空间尺度分析植被覆盖对大气颗粒物浓度的影响，相关结论如下。
1)NDVI月均值呈显著的周期性变化规律，3至8月稳定上升，8至11月则快速下降，11月至次年3月稳定在较低覆盖度。植被覆盖时间变化对大气颗粒物浓度有显著影响。分别使用不同曲线模型对NDVI值与颗粒物浓度进行回归分析，发现指数函数或幂函数最能反映二者关系，研究区整体NDVI与PM2.5、PM10浓度月均值回归方程的R2分别达到0.580、0.601。
2)植被对于颗粒物浓度的影响作用同样存在明显的空间效应，分别统计0.5～5 km范围内10种尺度半径缓冲区内的NDVI均值，并与颗粒物浓度年均值、季均值进行相关分析。结果表明，植被覆盖对PM2.5和PM10浓度年均值的最佳影响尺度分别为3 km和2.5 km，各季节的最佳影响尺度中夏季最大、冬季最小，表明茂盛的植被对于颗粒物浓度的影响作用较强、范围较大。
3)植被区颗粒物浓度在各季节均显著低于非植被区，非植被区四季颗粒物浓度均在一日内表现出明显的双峰双谷规律，而植被区则表现为单峰单谷或不明显的双峰双谷规律。植被覆盖既有效抑制了颗粒物浓度一日内两次升高的幅度，又显著缩短了浓度上升的持续时间，使颗粒物浓度的变化比非植被区更加平缓。
本文分别从时间、空间角度入手，分析了NDVI和PM2.5、PM10的时空特征，为探究植被覆盖对大气颗粒物浓度的影响提供了思路和方法。但由于数据源或实验方法的限制，实验可能存在一定不足：实验数据的时间为2015—2018年，时间跨度有限；
颗粒物浓度空间分布是基于实测数据，通过克里金方法空间插值获得，由于站点数量有限、空间分布不均，插值结果可能存在一定误差；
实验主要聚焦于区域宏观尺度，而对于颗粒物的来源、植被削减颗粒物作用机理等微观过程未作深入分析；
植被覆盖对颗粒物浓度的最佳影响尺度的科学意义仍有待进一步挖掘。以上均是后续研究中重点关注的方向。