王 戈， 于 强*， Yang Di， 牛 腾， 龙芊芊

1. 北京林业大学精准林业北京重点实验室， 北京 100083

2. Geographic Information Center, University of Wyoming, Laramie 82070, USA

引 言

与城市化相关的人类活动产生了一系列生态环境问题。 城市空气污染已成为严重的环境问题之一， 对人体健康造成了严重影响， 可引发心脏病、 肺癌、 哮喘和急性下呼吸道感染等疾病。 世界卫生组织估计， 每年有430万人因空气污染死亡； 如果可吸入颗粒物(pm2.5和pm10)从年平均70 μg·m-3降至年平均值20 μg·m-3， 预计将减少15%的死亡人数[1]。 减少大气颗粒物和尘埃粒子对于保障城市居民身体健康和提升城市生态环境质量至关重要。 研究表明[2]， 植被具有滞留大气颗粒物和尘埃粒子的作用， 一方面植物叶片表面绒毛， 油脂及粘性物质可以吸附一部分大气颗粒物和尘埃粒子， 另一方面植物群落可以降低附近风速， 增加空气湿度， 促使大气颗粒物和尘埃粒子降落。 因此以自然植被和人工植被为主要组分的绿色空间， 在吸滞大气颗粒物和尘埃粒子方面可以发挥一定作用。 在叶片尺度的滞尘研究上， 王会霞等[3]采用模拟降尘法测定植物叶片最大滞尘能力， 并探究了不同叶片性状对滞尘能力的影响； 贾彦等[4]利用电镜观测叶表结构， 从微观角度研究了叶片结构特征对滞尘能力的影响； 高君亮等[5]采用水洗过滤法测定单位叶面积滞尘量， 对比灌木与乔木滞尘效果， 对比了植物在枝叶和树冠等方面特征差异和植物滞尘量的相关性， 得到结论为： 灌木滞尘能力普遍大于乔木树种， 合理提高滞尘能力强树种比例可有效提高绿地滞尘能力。 罗娜娜等[6]以叶片光谱反射率曲线和偏最小二乘回归模型为基础， 分析了叶片尺度上尘埃量对植物光谱特征的影响。 谢英赞等[7]采用相关性分析法， 研究了同一树种在不同绿地类型中滞尘能力差异。 朱济友等[8]探究了大叶黄杨在不同尘埃浓度影响下, 叶面光谱及其叶片功能性状的响应机制及其权衡策略, 建立了滞尘量预测模型。 在城市或区域尺度上， 构建模型评估绿色空间的滞尘功能逐渐成为研究热点； 周志翔等[9]应用景观生态学原理和对比分析方法研究了绿地景观类型的空间结构及滞尘效应； 刘文平等[10]提出了由污染物浓度、 污染物沉降到叶表面的速率、 叶面积、 滞尘返还率、 滞尘时间共同决定的滞尘模型， 为量化绿色空间滞尘服务； Daniela等[11]利用城市森林效应评估法(UFORE)评估了城市植被吸收大气PM2.5效率为2.5%； Antoine等[12]基于空气动力学弥散效应模型评估了绿色空间对大气颗粒物的吸滞能力。 目前， 针对区域尺度绿色空间滞尘研究较少。

北京作为中国政治中心和京津冀城市群特大型城市， 近40年来城市化进程迅速， 大气颗粒物和尘埃粒子污染问题严重， 积极发挥北京市绿色空间滞尘功能具有现实意义。 研究时间段内落叶植被尘埃吸附能力较弱， 针叶植物滞留尘埃能力优于常绿阔叶植被， 但针叶植被在空间上分布不均。 据统计， 北京灌木种植面积是常绿针叶植被的3.23倍， 其中大叶黄杨种植面积占常绿灌木面积的78.03%[14]。 以北京市常见绿化树种大叶黄杨为研究对象， 利用地面采集叶片光谱数据和同一时期卫星遥感数据， 构建滞尘分布预测模型， 反演了北京市区绿色空间滞尘分布。 比较叶片滞尘前后的光谱曲线特征， 筛选对滞尘量高度相应的卫星波段， 对比所构建的十种指数与滞尘量间相关性， 构建滞尘反演模型。 针对地面采集的窄光谱发射率数据， 采用光谱响应函数将其转化为Sentinel-2影像的宽波光谱反射率数据。 利用滞尘反演模型获取绿色空间滞尘分布， 并获取尘埃污染程度分布， 尝试构建了一种快速估算绿色空间滞尘量方法， 为北京市绿色空间建设和城市环境治理提供参考。

1 实验部分

1.1 研究区

北京市位于北纬39°54′20″， 东经116°25′29″， 总面积约为16 410.54 km2， 年均气温为11～13 ℃， 年降水量为500～600 nm之间。 北京市季风特征明显， 全年降水量季节分布不均， 夏季降水量约占全年3/4， 其他季节较为干燥。 截至2020年末， 北京市常住人口2 189.3万人， 常住人口密度达1 289人·km-2。 全年客运量为36 255.8万人， 全年货运量为26 346.2万吨， 全市机动车保有量657万辆。 北京森林覆盖率为44.4%， 人均公共绿地面积为16.4 m2。 当前， 北京市绿色空间总量不足且分布不均匀、 人均公园绿地面积不足、 城市绿化覆盖率不高、 制约了城市环境改善[13]。

1.2 方法

1.2.1 数据采集与处理

选取北京市区绿色空间常见植被大叶黄杨(Euonymusjaponicus)作为试验对象， 大叶黄杨是多年生常绿灌木， 黄杨科黄杨属植物， 为中国特有植物， 具有耐阴抗旱、 枝叶茂密、 成活率高、 吸附有害物质能力强等特点。 基于广泛性、 均匀性、 规律性等采样点选取原则， 并参考先前研究成果， 分别设置实验采样点与验证采样点数量为50和25。 采样点均位于北京市区不同的绿色空间内， 采样点空间分布如图1。 采样时选择健康、 成熟与无损害的成熟叶片， 在每个采样点采集30个叶片。 在2021年4月7日上午8—10点进行采样， 采用透明密封袋盛装样本， 并置于恒温储藏箱中(温度6 ℃， 湿度80%)， 以保证叶片生物理化活性。

图1 采样点分布

按照上述原则所采集的叶片， 采用电子分析天平(精度0.000 1 g)称重， 叶片除尘前质量为W1， 除尘后叶片质量为W2， 叶片滞尘质量为W1和W2之差， 使用超纯净水冲洗叶片表面尘埃并采用吸水纸进行干燥处理。 将叶片经过扫描仪扫描后， 采用LA-S植物图像分析仪测定其叶面积S。 对不同叶片面积、 除尘前后质量进行多次测定求平均值， 从而获取单位面积滞尘量P(g·m-2)， 计算公式为

(1)

采用04多功能光谱仪进行叶片反射率光谱测定， ADS光谱仪波段范围为350～2 500 nm， 约有2 100个波段， 采样间隔为1 nm。 在室内进行叶片光谱测定， 光源为1 200 W卤素灯， 将光源置于距叶片表面70 cm位置， 与水平面夹角为70°， 将光谱仪镜头放置于叶片上方15 cm位置， 且与水平面保持垂直， 镜头视场角为25°。 对ADS光谱仪进行30 min预热， 用标准白板进行标定后开始测定叶片光谱， 读数稳定后保存数值， 每隔15 min重新优化校正。 每个叶片采集5条光谱曲线， 利用ASD光谱处理软件ViewSpecPro查看并剔除明显错误的光谱曲线。 对除尘前光谱反射率RDust和除尘后光谱反射率Pclean进行多次测定求平均值， 除尘前后光谱反射率比值γ， 计算公式为

(2)

1.2.2 遥感数据与处理

一般认为， 当降水量达15 mm或风速达17 m·s-1可以去除叶片表面尘埃， 可以模拟实验室人工叶片除尘。 6月2日白天持续降雨并伴随大风天气， 6月3日为晴天， 因此将6月3日影像作为无尘植被影像， 选择4月7日作为植被滞尘影像。 从Copernicus Open Access Hub(哥白尼开放访问中心https://scihub.copernicus.eu/)获取本研究所需的Sentinel-2遥感数据， Sentinel-2影像有13个波段， 其中2， 3， 4和8波段的空间分辨率为10 m， 5， 6， 7， 8 A波段的空间分辨率为20 m， 1， 9和10波段的空间分辨率为60 m。 Sentinel-2数据已经过几何校正与辐射定标， 为避免大气吸收与散射所造成的辐射量误差， 采用SNAP软件对下载数据进行大气校正。

图2 室内叶片滞尘量及光谱测定过程(a)和

1.2.3 绿色空间识别

城市绿色空间是由园林绿地、 城市森林、 立体空间绿化、 都市农田和水域湿地等构成的绿色系统[15]。 将北京市区绿色空间界定为林地、 草地和农田为组分的绿色开敞空间， 采用国产GF-2遥感数据进行北京市区绿色空间识别， 三景影像数据采集时间是2020年9月2—3日， 数据来源于中国资源卫星应用中心-数据服务平台。 利用ENVI软件进行影像预处理， 主要包括： 正射校正、 辐射定标、 图像融合和大气校正。 使用2016年World-view Ⅱ影像为配准参照， 采用Gram-Schmidt Pan Sharpening模块进行GF-2全色波段和多光谱影像融合， 采用人工目视解译法， 将土地覆被类型分为农田、 林地、 草地、 湿地与建设区， 生成北京市土地覆被分布(图3) 。

图3 2020年北京市区土地覆被及绿色空间分布图

1.2.4 植被指数

遥感植被指数可以有效表征植被状况， 选取10种遥感植被指数， 计算滞尘前后植被指数比， 本研究采用植被指数比值反演滞尘含量， 具体表达式如表1。

表1 植被指数公式

1.2.5 空间自相关模型

利用ArcGIS10.4 Create Fishnet工具将研究区划分为150 m×150 m正方形网格单元， 共计38 415个， 计算每个网格单元中单位面积滞尘量平均值。 采用全局空间自相关和局部空间自相关来刻画空气污染程度的空间聚集特征[16]。

全局自相关用于反映空气污染程度空间关联程度的总体特征， 其值介于-1和1之间， 表达式如式(2)

(2)

局部自相关用于反映每个正方形网格单元之间的空间关联程度， 通过莫兰散点图将其划分为高-高(H-H)、 高-低(H-L)、 低-高(L-H)和低-低(L-L)四种类型。 其中， H-H(L-L)聚集型表示相邻网格单元间存在正的空间自相关， 即存在高(低)空气污染程度的空间聚集效应， H-L(L-H)表示相邻格网单元存在负的空间自相关， 高(低)空气污染程度被低(高)空气污染程度网格单元包围。

(3)

2 结果与讨论

2.1 滞尘对叶片光谱的影响

叶面滞尘会导致植物叶片光谱曲线发生变化。 在暗室中利用ASD光谱仪测定“除尘前”和“除尘后”叶片光谱曲线， 分别计算除尘前叶片与除尘后叶片的不同波段光谱反射率均值， 结果如图4(a)所示。 除尘前后叶片光谱反射率曲线存在一定差异， 但是光谱曲线变化趋势基本相同。 除尘前后叶片光谱曲线在可见光波(380～780 nm)和近红外波段(780～1 100 nm)存在两个吸收峰和三个吸收谷。 其中， 反射峰位于553和1 088 nm， 是植被典型的光谱反射带。 吸收谷分别位于390～500， 630～690和975～980 nm， 是植被典型的光谱吸收带。 在670～760 nm之间， 光谱曲线斜率及光谱反射率变化较大， 出现了“红边效应”。 在760～1 425 nm之间出现一个较高的“反射平台”， 这是由于叶片海绵组织对近红外波段具有较强反射率。 在350～700, 740～1 870和2 000～2 400 nm波段, 除尘后叶片光谱曲线均高于除尘前， 可能由于尘埃对光子的吸收与散射作用较强有关。

图4 滞尘对光谱特征的影响(a)和除尘前后

光谱曲线的红边、 黄边、 蓝边是反映植被生长状况的重要参数。 三边参数相关含义见表2， 统计除尘前后叶片光谱反射率一阶导数。 如图4(b)所示， 滞尘对红边、 黄边、 蓝边位置没有明显影响， 蓝边为521 nm， 黄边为628 nm， 红边为704 nm， 叶面滞尘对三边位置干扰不敏感。 如表2叶片滞尘对蓝边与黄边幅值影响不明显， 对红边幅值影响明显。 叶片滞尘会导致三边面积降低， 但是红边面积对滞尘敏感性大于蓝边和黄边面积， 这是由于叶面滞尘大多为固体颗粒， 这些物质会大幅降低绿色植物对近红波段敏感程度。

表2 滞尘对光谱蓝、 黄、 红边的影响

2.2 叶片滞尘量与大叶黄杨叶片光谱数据的相关分析

2.2.1 叶片滞尘量与单波段反射率相关分析

计算50个叶片光谱反射率比值(除尘前/除尘后)和对应叶面滞尘量的相关度， 结果如图5所示。 由图5可知， 光谱反射率比值和滞尘量在354～705 nm波段为正相关， 最大正相关系数位于653 nm波段， 其值为0.420， 在510～590 nm处明显凹形， 该波段位于绿光波段内， 表明绿光波段对尘埃的敏感性低于其他波段。 706～1 110 nm波段为负相关， 最大负相关系数位于801 nm波段， 其值为-0.729。 由相关性曲线分析可知， 红波段和近红波段对滞尘高度敏感。

图5 滞尘量与光谱反射率比值(除尘前后)之间的关系

2.2.2 卫星波段光谱反射率与滞尘量的相关分析

针对地面采集的窄光谱发射率数据， 采用光谱响应函数将其转化为Sentinel-2， Landsat8， ZY-3， Spot6和MODIS卫星的宽波段光谱反射率数据。 选取350～1 110 nm波段内5种卫星的对应波段， 计算除尘前后卫星波段的植被指数比值， 将其与滞尘量进行相关分析， 结果如表3所示。 由表3可知， 不同卫星植被指数比值与单位面积滞尘量相关性变化具有相同趋势， 除PVI指数外， 其余9个植被指数均与滞尘量为负相关。 Sentinel-2的EVI指数相关性最大， 其相关系数为-0.809。 因此， 选择Sentinel-2作为遥感数据来源。

表3 五种卫星的植被指数比值与滞尘量之间的关系

2.3 绿色空间滞尘量估测模型构建

由于Sentinel-2的EVI指数与滞尘量相关性最高， 建立EVI指数滞尘前后比值与滞尘量的回归模型如图6所示。 如表4所示， 线性拟合模型与二次拟合模型的决定系数(R2)分别为0.705和0.751， 均方根差(RMSE)分别为1.69和1.298， 二次回归模型精度较高， 利用EVI指数建立的回归模型可以较为准确反演北京市区绿色空间滞尘分布。 因此， 采用EVI指数二次回归模型作为滞尘反演模型。

图6 基于EVI指数的反演模型

表4 滞尘反演模型

2.4 植被滞尘反演与验证

利用EVI指数滞尘反演模型进行反演研究。 采用ENVI5.3软件， 将2021年4月7日EVI指数与2021年6月3日EVI指数比值作为自变量， 采用EVI指数滞尘反演模型进行叶面滞尘量反演， 并结合北京市区绿色空间分布， 获取绿色空间滞尘空间分布。 4月7日北京市区绿色空间滞尘量分布如图7(a)， 在整体上滞尘量分布趋势表现为： 北部滞尘量比南部高， 东部比西部高， 中心城区高于郊区。 道路两旁的绿色空间滞尘量较高， 绿色空间单位面积与滞尘量距道路距离成反比。 同一区域内， 绿色空间内部受尘埃影响程度受绿色空间斑块面积影响较大， 其斑块面积与内部滞尘量成反比。 市区中部的绿色空间滞尘量较高， 因为附近分布有北京西站、 北京站与北京南站， 人流与交通流量较大。 据统计， 2016年北京西站、 北京站和北京南站日均客流量约为46万人次， 聚集的人流与交通流量会导致颗粒物污染。 望京地区附近绿色空间滞尘量较高， 该区域人员与商业区聚集。

图7 北京市绿色空间的滞尘分布(a)和滞尘反演精度验证(b)

将25组滞尘反演量和对应实测值进行相关分析， 利用决定系数(R2)与均方根差(RMSE)评价反演模型精度。 结果如图7(b)所示， 采用滞尘量反演预测值与实测数据具有很高拟合度，R2为0.807， RMSE为6.287， 表明EVI指数反演模型可以较为准确反映北京市区绿色空间滞尘量空间分布。

2.5 北京市空气污染程度空间分析

植被滞尘分布可以作为表征空气污染的指标， 植被滞尘量一定程度可以反映尘埃污染程度。 在绿色空间范围内选取1 500个采样点， 利用ArcGIS Extract Values to Points工具， 获取采样点位的滞尘量数值， 采用克里金(Kriging)插值法获取北京市区尘埃污染分布。 如图8(a)所示， 4月7日尘埃污染集中分布在北京市区中部、 北部与东部， 在市区西部与南部， 尘埃污染分布较为分散， 表明北京市区中部、 北部和东部易产生尘埃污染， 污染区域大多分布在人口和商业密集区域附近。

图8 北京市区大气尘埃污染分布和空间自相关分析

北京市区尘埃污染分布的全局Moran’I值0.912，p值小于0.05， 说明研究区域尘埃污染分布在整体上呈现显著的正向空间自相关关系， 具有明显聚集性， 空间分布并非完全随机。 由图8(b)可知， 植被滞尘 LISA分布图(p<0.05)直观反映植被滞尘在空间聚集与分异的位置分布特征， 在北京市区中部， 东部和东北部出现了大面积的尘埃污染高-高聚集(H-H聚集)区域， 尘埃污染低-低聚集(L-L聚集)主要分布在北京市区西部和南部。

3 结 论

以北京市区为研究区域， 以大叶黄杨为研究对象， 经过室外采样与室内试验获取了叶片样本滞尘量、 光谱反射率和叶面积等数据， 分析了滞尘程度对光谱反射率影响， 确定了滞尘响应的高度敏感波段， 建立了对应卫星波段植被指数比值与滞尘量的回归模型， 并结合绿色空间分布数据， 采用Sentinel-2影像获取北京市区绿色空间滞尘空间分布， 并获取了尘埃污染分布， 采用空间自相关模型检验其空间聚集特性。

(1) 在740～1 870 nm波段, 滞尘后叶片光谱反射率明显低于滞尘前叶片。 滞尘对红边、 黄边、 蓝边位置没有明显影响， 对“红边幅值”和“红边面积”影响较明显。 叶片光谱除尘前后反射率比值和滞尘量的相关分析表明： 在354～705 nm波段为正相关， 706～1 110 nm波段为负相关， 最大正相关系数和最大负相关系数分别位于红光波段和近红光波段范围内， 红和近红波段对滞尘高度敏感。

(2) 进行宽波段滞尘前后植被指数比值和单位面积滞尘量相关分析： Sentinel-2， Landsat8， ZY-3， Spot6和MODIS的相关性相差较小。 其中， Sentinel-2计算的EVI指数比值与单位面积滞尘量相关系数最大， 构建对应的滞尘反演模型， 所构建的线性与二次回归模型决定系数分别为0.705和0.751。 因此， 选取EVI指数二次回归模型作为滞尘反演模型。

(3) 利用2021年4月7日和6月3日的Sentinel-2影像， 采用EVI指数滞尘反演模型获取了北京市区绿色空间滞尘分布， 并利用实测数据对滞尘反演结果的合理性进行验证， 绿色空间的滞尘量分布表现为中心高于郊区， 东部高于西部， 北部高于南部。 决定系数(R2)与均方根差 (RMSE)分别为0.807和6.287， 反演模型可以较为准确反映北京市区绿色空间滞尘量空间分布。 北京市区中部、 北部和东部易产生尘埃污染， 且污染分布具有明显聚集性， 并非完全随机。