PM_2.5污染在全球气候系统的驱动及交互影响下可引起巨大的社会经济损失, 在公共卫生领域承担着较大疾病负担贡献率^{[1, 2]}.为缓解PM_2.5污染及健康风险, 世界卫生组织^[3]、美国环保署^[4]和中国国家环境保护部^[5]等环境组织和国家制定了一系列干预措施与法规, 污染改善效果明显^[6].然而, 长三角地区由于人口密集、产业发展迅速、城市急剧扩张和人类活动活跃仍然为中国东部地区重要雾-霾区^[7].因此, 深入探讨长三角地区PM_2.5时空演变特征, 识别关键影响因素及其响应过程具有重要意义.然而, 当前多利用多时间尺度、时间序列数据及多空间尺度、空间截面数据进行统计分析, 而获取时空动态分异特征^{[8, 9]}, 忽略了地理过程的时空不可分离性, 缺乏对时空变化特征抽象规律的提取, 而聚类是时空模式提取的重要手段^[10].

PM_2.5浓度变化受自然因素与人为活动综合影响, 鉴于对PM_2.5源解析及间接传输扩散条件的物理作用的认识, 多数研究者从工业生产、交通排放和地面扬尘等直接排放因素^[11]、扩散沉降条件因素(地形、气象条件、地表条件)^{[12, 13]}和人类活动强度(城市化过程、城市形态)^{[14, 15]}等间接因素方面探讨PM_2.5驱动机制.而土地利用类型和景观格局可囊括自然/人工源与人类活动强度^{[16, 17]}而作为自然与人为复合因素, 其对PM_2.5的影响机制尚待深入.由于不同土地类型作为地表下垫面存在强烈的物理性质差异, 与PM_2.5浓度表现为显著的空间响应关系^[18], 其土地利用空间分布及景观结构的变化一方面会改变水热条件并与大气活动交互作用影响颗粒物的动态传播^[19], 另一方面会改变地表源-汇景观系统的空间格局^[20], 均会对PM_2.5浓度产生直接或间接的影响.目前针对土地利用与PM_2.5关系的探讨多采用多元线性回归^[20]、广义可加模型^[21]和空间计量模型^{[21, 22]}等.而这些传统统计分析模型易受限于变量选择对结果的敏感性, 如何选择合适的变量组合是关键.在识别关键影响因子中, 随机森林模型因具备运行高效、抗干扰与泛化能力强和无多重共线性, 在土壤属性^[23]和PM_2.5浓度^[24]等方向得到了较好地模拟效果.此外, 目前对PM_2.5浓度与土地利用空间关系的探索较少考虑到空间非平稳性及作用范围的差异性的影响^{[16, 17, 19]}, 较难把握PM_2.5的空间异质性与区域空间多尺度影响过程.因此, 本研究以我国长三角地区为例, 结合随机森林模型、空间计量模型和多尺度地理加权回归模型多角度尝试阐明土地利用/覆盖类型及转移对PM_2.5浓度的空间响应过程.其结果可精确分析土地利用/覆盖因素对PM_2.5的空间响应, 从区域空间联合角度为联防联控的政策提供精细理论支持, 期望助力共建人与自然生命共同体.

长三角地区工业化与城市化迅速, 人类活动密集, 人为污染源排放强度大, 土地覆盖变化明显, 是我国PM_2.5污染最为严重的区域之一, 在空间上呈现东南高西北低的格局(图 1).我国自2010年开展了区域性大气污染联防联控治理, 长三角地区在国家战略高度上于2013年成立区域大气污染防治协作机制, 以排放控制为主实施大气污染防治行动, 效果显著, 2013~2016年期间PM_2.5浓度下降了27%^[25], 但仍然有超过一半的区域高出年均浓度限值, 可见长三角地区PM_2.5污染并未得到根本性解决.长三角地区地形起伏、气候条件差异不明显, 但土地利用变化剧烈, 明确高强度人类活动影响下土地利用/覆盖转换对PM_2.5浓度的影响机制将为PM_2.5污染防控提供重要参考.

图 1

Fig. 1

图 1 2000~2018年研究区PM2.5浓度变化及土地利用转换特征 Fig. 1 PM2.5 concentration change and land use conversion characteristics in the study area from 2000-2018

考虑到地面监测站空间分布的不均匀性, 选取大气成分分析组织的中国区域地面PM_2.5空间分布数据集(V4.CH.03), 空间分辨率为0.01°×0.01°, 该数据集是通过化学传输模型(GEOS-Chem)基于NASA MODIS、MISR和SeaWIFS等卫星产品的AOD与PM_2.5之间的地球物理关系进行估算浓度值并利用地面测量数据校准^[26].为验证该数据集在研究区域的可靠性, 本文利用生态环境部发布的2015~2018年长三角区域地面监测站的PM_2.5指标计算站点年平均PM_2.5浓度, 通过空间匹配提取出对应的PM_2.5浓度值, 结果表明二者间的相关系数为0.732~0.927, 均方根误差为3.576~6.197 μg·m^-3, 信息损失不大, 能够捕捉到PM_2.5浓度变化信息, 可用于长三角区域PM_2.5浓度空间分析.

研究区土地利用/覆盖数据来自于中国科学院资源环境科学数据中心(http://www.resdc.cn/), 空间分辨率为1 km, 时间上选取2000年与2018年.根据研究目的, 选取一级土地分类, 计算10 km×10 km格网单元内2000~2018年土地利用/覆盖转换类型的面积比例及格网单元内PM_2.5浓度平均值, 共获得3462个样本.

K-means聚类算法经多项研究检验可产生稳定的聚类边界, 考虑到时间变化特性, 采用余弦距离作为相似性度量, 是研究地理要素时空分异特征的重要手段之一^{[10, 27]}.另外, 聚类效果仍需体现PM_2.5浓度的空间分异性特征, 地理探测器则是一种广泛使用的空间分层异质性探测的空间统计方法, 通过q统计量量化空间分异与分层性来评估聚类内部与聚类之间空间上的相似性与差异性^{[28, 29]}.本研究将时间序列的PM_2.5数据中具有相似时间变化的格网分组在一起, 以每年原始PM_2.5浓度值为因变量, 聚类结果为自变量计算q统计量, 以获取聚类结果的平均q值, q值计算过程见王劲峰的研究^[28].越趋近于1的q值表示聚类结果有效地区分了长三角地区PM_2.5浓度的区域分异.

随机森林(RF)模型是由Breiman提出的一种基于决策树的Bagging典型算法^[30], 该模型的回归过程是通过Bagging随机且有放回的采样方法将原数据集按一定比例进行多轮抽取形成多组训练集, 对每个训练集分别进行随机选择特征进行递归构建相应决策树, 所有决策树结果的均值即为回归结果^{[24, 30]}.随机森林回归模型基于本身机器学习的算法具有运行高效、对异常值不敏感和无需特征选择等多个独特优势, 同时其变量重要性排序可增强模型解释能力^[24].本文以3462个样本所对应的2000~2018年PM_2.5浓度变化量与各土地利用/覆盖转移类型的面积比例, 分别作为因变量和解释变量构建随机森林模型, 通过变量重要性函数识别影响PM_2.5浓度变化的关键变量.

PM_2.5浓度变化受大气传输及区域聚集性人类活动的影响, 其空间分布存在较强的空间自相关性, 不同的土地利用/覆盖转换类型因其下垫面的性质变化及空间影响对PM_2.5浓度存在较强的空间溢出效应, 相邻地区的土地利用/覆盖转换及相邻地区PM_2.5浓度的变化均能影响本地的变化趋势^{[21, 22]}.因此, 通过引入空间回归模型量化不同土地转换类型对相邻地区的空间溢出效应信息, 包括空间误差模型(SEM)、空间滞后模型(SLM)和空间杜宾模型(SDM)^[31].本研究应用空间相关性来估计约束模型^{[22, 31]}, 利用拉格朗日检验各土地利用/覆盖转换类型对PM_2.5浓度变化的空间滞后和空间误差影响, 检验结果显示在0.01的置信水平下拒绝了空间滞后和空间误差影响的假设, 因此选择纳入解释变量和被解释变量的空间滞后变量构建空间杜宾模型.空间杜宾模型表达式为:

(1)

式中, Y为PM_2.5浓度变化值, β为自变量的回归估计系数, 系数值表示土地利用/覆盖转换方式对PM_2.5浓度变化的影响方向及程度, ρ为因变量Y空间滞后变量的空间自回归系数, WX为空间权重矩阵, W为自变量空间加权后对应的空间滞后自变量, γ为空间滞后自变量的回归系数, 其系数值表示相邻单元内土地利用/覆盖转换方式对本单元内PM_2.5浓度变化的影响, ε为随机误差项.

地理加权回归(GWR)是基于最小二乘线性回归(OLS)进行局部优化并使用空间关系加入运算的建模工具, 可以通过局部参数量化自变量和因变量之间的空间异质性^{[17, 32]}.考虑到各因子与因变量之间的空间关系在不同规模空间尺度下起不同作用与关系的尺度依赖性, Fotheringham等^[32]的研究在GWR基础上进行扩展, 放宽了模型中所有空间变化过程都在同一空间尺度上运行的假设(即所有解释变量均在相同空间尺度下作用于PM_2.5浓度), 为每个变量指定专用带宽以适应模型校准其多尺度效应, 该模型称之为多尺度地理加权回归模型(MGWR), MGWR模型可以表示为:

(2)

式中, (u_i, v_i)为样本点i的坐标, y_i和x_ik分别为因变量和第k个解释变量, 表示样本点i的第k个解释变量的局部参数估计, β₀(u_i, v_i)为截距, ε_i为误差项, bwk为第k个解释变量与自变量的条件关系的校准中使用的特定最佳带宽, 代表模型中解释变量与因变量之间的关系允许在空间上发生变化.

2000~2018年间, 长三角地区PM_2.5浓度变化分为4类模式时其轮廓系数(0.45)最大, 具备良好的簇内紧凑度和簇间分离度, 4类时变模式及空间分布见图 2.为表征聚类结果能否体现研究期间内年PM_2.5浓度的空间分异程度, 采用地理探测器q值评估聚类之间的空间差异, q值范围区间在[0.40, 0.73], 平均q值达0.59, 显著高于Chen等^[29]对中国区域各季节PM_2.5浓度时空聚类结果的平均q值(0.23~0.31).从时间变化趋势上看, 4类时变模式分别表现为长时间平稳后波动下降(Cluster 1); 剧烈波动上升后波动下降(Cluster 2); 非剧烈波动上升形成短暂平台后波动下降(Cluster 3); 非剧烈到剧烈波动上升形成稳定波动后波动下降(Cluster 4), 在浓度值上呈现分层模式.同时时变模式在空间上存在空间聚集性, 形成长三角地区南部, 中部, 西北部和东北部空间分异形态.在本聚类算法中并没有将空间属性纳入计算, 其结果的空间连续性及聚集性再一次体现了PM_2.5浓度变化的空间效应, 距离较近的浓度变化趋势是相似的, 而空间距离较远的浓度变化呈现较大的差异, 表明PM_2.5当地的污染改善就不仅仅取决于当地污染情况, 更取决于区域性的同步变化.

图 2

Fig. 2

图 2 长三角地区PM2.5浓度时空变化模式 Fig. 2 Spatial and temporal patterns of PM2.5 concentrations in the Yangtze River Delta region

人为排放及气象条件的变化对PM_2.5浓度存在的贡献不可忽略^[33].因此, 为剔除土地利用/覆盖转换之外因素的影响, 以长三角地区县域各地类未发生变化的PM_2.5浓度变化值为标准进行修正, 计算土地利用/覆盖转换的PM_2.5浓度变化值的相对量.由于未利用地占比较小, 仅考虑耕地、林地、草地、水体和建设用地的转换对PM_2.5浓度的相对影响.未发生土地利用/覆盖转换的PM_2.5浓度变化及空间计算修正后的不同土地利用/覆盖转换方式下的PM_2.5浓度的变化见表 1.

表 1 (Table 1)

表 1 不同土地利用/覆盖转换方式下的PM2.5浓度的相对变化量/μg·m-3 Table 1 Relative change in PM2.5 concentration under different land use conversions/μg·m-3 土地利用/覆盖转换类型 转换前土地类型的平均PM2.5变化量 修正变化量 耕地→林地 0.059 -0.286 耕地→水体 0.059 -0.160 耕地→建设用地 0.059 -0.044 林地→耕地 -13.386 0.052 林地→草地 -13.386 0.628 林地→建设用地 -13.386 0.136 草地→耕地 -8.147 1.563 草地→林地 -8.147 -0.016 水体→耕地 0.167 0.071 水体→建设用地 0.167 -0.027 建设用地→耕地 1.119 -0.081 建设用地→水体 1.119 0.680

表 1 不同土地利用/覆盖转换方式下的PM2.5浓度的相对变化量/μg·m-3 Table 1 Relative change in PM2.5 concentration under different land use conversions/μg·m-3

土地利用未发生转换时, PM_2.5浓度在耕地、建设用地、水体呈增幅状态, 林地与草地呈减幅状态, 林地与草地对PM_2.5浓度的缓解作用程度明显, 与多项土地利用/覆盖类型与PM_2.5关系研究结果一致^{[19, 22]}.土地转换对PM_2.5浓度的相对影响表现复杂.耕地向其他用地类型转换会引起PM_2.5浓度下降, 林地的转出则会引起PM_2.5浓度的上升, 耕地与林地的源-汇效应明显, 其他转换类型无一致性转换方向.PM_2.5浓度变化增幅最大的是草地转耕地及建设用地转水体, 放大了耕地的源作用和水体正向作用.水体由于吸湿作用对二次颗粒物形成、累积、扩散等过程影响对PM_2.5浓度的作用存在多向性^[34], 但水体的转入对PM_2.5浓度变化作用显著, 表现出倾向性的颗粒物累积作用.耕地、建设用地等人工性质下垫面的转出反映出向林地、水体等自然性质下垫面的转换对PM_2.5浓度的作用显著强于人工用地之间的相互转换.

土地利用/覆盖对PM_2.5浓度的影响过程并不是局部独立进行的, 地表的物理性质(尤其是热力性质)差异与大气传输的相互作用使得土地利用/覆盖对PM_2.5浓度的影响存在较强的空间溢出效应^[23], 区域性同步变化再次强调了PM_2.5浓度变化存在对周围环境变化的响应信息.因此, 通过对单一像元和周围聚集性土地利用/覆盖转换的PM_2.5浓度变化进行对比来明确周围的土地利用/覆盖转换类型对PM_2.5浓度的影响(表 2).结果表明, 周围聚集性未发生土地利用/覆盖转换区域耕地、林地、水体和建设用地比单一像元时对PM_2.5浓度变化的作用表现更强烈, 而草地的聚集性作用出现了弱化的情况, 但草地的变化程度最大.而周围聚集性发生土地转换区域草地向耕地的转换出现弱化, 林地向耕地出现反向情况(降低→升高), 其他转换类型均比单一像元时作用更为显著, 尤其以建设用地向耕地的转换与水体向草地的转换为甚, 周围聚集影响的PM_2.5浓度变化量是单一像元的近4倍.

表 2 (Table 2)

表 2 周围聚集性土地利用/覆盖转换对PM2.5浓度变化量的影响/μg·m-3 Table 2 Effect of peripheral aggregated land use conversion on the change in PM2.5 concentration/μg·m-3 转换类型 PM2.5浓度变化量 自身影响 周围聚集影响 耕地→耕地 0.059 1.223 耕地→水体 0.477 3.049 耕地→建设用地 -1.089 -1.220 林地→耕地 -3.386 1.828 林地→林地 -13.386 -13.969 林地→建设用地 -11.026 -13.861 草地→耕地 11.386 10.403 草地→草地 -8.147 -6.087 草地→水体 11.677 15.448 水体→耕地 2.775 8.310 水体→草地 3.153 12.628 水体→水体 0.167 1.748 水体→建设用地 2.232 6.219 建设用地→耕地 4.685 18.451 建设用地→水体 4.211 7.817 建设用地→建设用地 1.119 1.733

表 2 周围聚集性土地利用/覆盖转换对PM2.5浓度变化量的影响/μg·m-3 Table 2 Effect of peripheral aggregated land use conversion on the change in PM2.5 concentration/μg·m-3

相关系数表明PM_2.5浓度变化与大部分土地利用/覆盖转换类型均具备显著关系, 但强度与方向略有差异(表 3).其中, 在未发生土地转换的类型中, PM_2.5浓度变化与林地、草地呈显著负相关, 而与水体、建设用地和耕地呈显著正相关, 与多项研究证实的单一土地利用类型对PM_2.5的作用方向一致^{[19, 21]}, 同时耕地的正效应比建设用地更强, 表明农业活动比城市活动对PM_2.5浓度贡献更大.发生土地转换的类型中, PM_2.5浓度变化与林地、草地发生转换的类型多呈显著负相关, 而PM_2.5浓度变化与耕地、建设用地、水体之间的转换类型呈显著的正相关.另外, 值得注意的是林地和耕地之间的转换与PM_2.5浓度变化之间的关系强于林地和草地之间的转换, 接近于林地未发生转换的作用强度, 表明林地与耕地的相互转换中耕地的农作物植被效应发挥了较大的作用.

表 3 (Table 3)

表 3 PM2.5浓度变化量与土地利用/覆盖转换的相关系数1) Table 3 Correlation coefficient between the amount of change in PM2.5 concentration and land use conversion 类型 相关系数 类型 相关系数 类型 相关系数 类型 相关系数 类型 相关系数 V11 0.650** V21 -0.581** V31 -0.008 V41 0.251** V51 0.519** V12 -0.602** V22 -0.726** V32 -0.344** V42 -0.183** V52 -0.119** V13 -0.069** V23 -0.336** V33 -0.100** V43 0.082** V53 0.060** V14 0.271** V24 -0.215** V34 0.134** V44 0.146** V54 0.134** V15 0.424** V25 -0.257** V35 -0.007 V45 0.145** V55 0.250** 1) **表示相关性在0.01水平显著(双尾检验); 类型编号中首位数字表示转换前地类代码, 末位数字表示转换后地类代码, 1为耕地, 2为林地, 3为草地, 4为水体, 5为建设用地, 下划线表示强度排序前10的变量

表 3 PM2.5浓度变化量与土地利用/覆盖转换的相关系数1) Table 3 Correlation coefficient between the amount of change in PM2.5 concentration and land use conversion

通过随机森林模型获得各土地利用/覆盖转换类型对PM_2.5浓度变化的重要性排序(图 3).未发生转换的耕地(29.65%)、林地(26.98%)与建设用地(18.00%)以及发生转换的耕地-林地(17.84%)、耕地-建设用地(16.34%)、建设用地-耕地(22.58%)均对PM_2.5浓度变化有重要影响.从未发生转换的土地类型上看, 耕地、林地的重要性均强于建设用地, 从发生转换的土地变化类型上看, 建设用地-耕地的重要性均强于耕地-建设用地.另外, 随机森林模型的重要性及相关系数的强度排序结果相比, 耕地-耕地(29.65%及0.650)、林地-林地(26.98%及0.726)、建设用地-耕地(22.57%及0.519)、耕地-林地(17.84%及0.602)和耕地-建设用地(16.34%及0.424)与PM_2.5浓度变化之间表现出一致性的较强关系, 表明耕地的转换对PM_2.5浓度变化的影响可能是污染控制的关键.

图 3

Fig. 3

类型编号中首位数字表示转换前地类代码, 末位数字表示转换后地类代码, 1为耕地, 2为林地, 3为草地, 4为水体, 5为建设用地; IncMSE为精度减少平均值, IncNodePurity为节点不纯度平均减少值, 分别从模型估算误差和决策树节点两方面判断变量重要性程度 图 3 随机森林模型识别出的重要变量及重要性排序 Fig. 3 Significant variables identified by the random forest model and ranking of importance

结合相关分析和随机森林模型筛选出对PM_2.5浓度变化具有重要影响的变量构建空间杜宾模型, 结果见表 4.PM_2.5变化的自回归系数(ρ)为正且显著, 表明PM_2.5浓度变化存在显著的空间依赖性, 周围网格单元的PM_2.5浓度的快速增加会显著加剧本单元的污染.另外林地-林地、林地-耕地、林地-草地、林地-建设用地、建设用地-耕地这5个自变量的空间滞后变量在SDM模型中均通过95%的显著性检验, 而且其系数比其自变量系数数值更大, 说明以上5个自变量对于PM_2.5浓度变化的影响不仅存在显著的空间溢出效应, 而且其空间溢出效应比本身的影响强度更强.此外, 将相关系数与SDM模型中自变量系数相比, 发现林地-林地、林地-草地的关系方向出现差异, 但相关系数与滞后变量系数方向一致, 表明林地-林地与林地-草地对PM_2.5的作用是由空间滞后效应所导致的, 而其他变量的作用是由类型本身所决定的.

表 4 (Table 4)

表 4 多元线性回归和空间计量模型回归结果1) Table 4 Results of multiple linear regression and spatial econometric model regression 变量名 OLS SDM 变量名 OLS SDM ρ / 0.987** W_V21 / 0.017* Intercept 3.16E-07 0.001 W_V22 / -0.058** V11 0.085** 0.012** W_V25 / 0.019** V12 -0.169** -0.014** W_V51 / -0.020** V21 -0.102** -0.007* W_V23 / -0.008* V22 -0.472** 0.043** LM (lag) 11 387.345 / V23 -0.047** 0.007* LM (error) 16 587.258 4 / V25 -0.078** -0.005* Robust LM (lag) 841.247 4 / V51 0.035** 0.010** Robust LM (error) 6 041.160 8 / 1)**表示P < 0.01, *表示P < 0.05; ρ表示PM2.5浓度变化量的空间自回归系数; 变量V11~V51中首位数字表示转换前地类代码, 末位数字表示转换后地类代码, 1为耕地, 2为林地, 3为草地, 4为水体, 5为建设用地; W_V21~W_V23表示自变量对应的空间滞后自变量; “/”表示无相关数据

表 4 多元线性回归和空间计量模型回归结果1) Table 4 Results of multiple linear regression and spatial econometric model regression

考虑到空间关系的非平稳性特性, 为获取PM_2.5浓度变化对土地利用/覆盖转换的影响范围差异, 以消除模型内生性问题为目标, 选择了耕地-耕地、耕地-林地、耕地-建设用地、林地-耕地, 林地-林地、林地-草地、林地-建设用地、建设用地-耕地和建设用地-建设用地这9个变量进行多尺度地理加权回归(MGWR)建模, 并通过调整后的R²、赤池信息量准则(AIC)、修正赤池信息量准则(AICc)、残差平方和(RSS)评估MGWR模型的拟合优度, 从残差的Moran'I指数评估模型是否具备解决空间自相关问题的能力^[16](见表 5).通过比较发现MGWR调整后的R²、AIC远远高于OLS和RF, 但略低于SDM模型, RSS值远远低于OLS、RF和SDM模型, 表明MGWR模型在空间要素的拟合优度方面仅对传统回归分析有较大的优势.这主要是由于SDM模型中因变量滞后变量的自回归大大提高了模型的拟合度^[31].但是, MGWR的残差表现为非常低的空间自相关, 与线性回归残差假设相符合, 而SDM模型残差仍表现出较强的空间自相关, 同时与Tu等^[16]的研究(应用GWR模型)进行比较, MGWR模型残差随机性更强, 模型结果更为可靠.这表示土地利用/覆盖转换与PM_2.5浓度变化量之间的关系存在强烈的空间非平稳性特征, MGWR模型在量化空间关系中具备较大的模型优势.

表 5 (Table 5)

表 5 模型效果评估1) Table 5 Model effectiveness evaluation 参数 OLS RF SDM MGWR Adjusted R2 0.601 0.685 0.985 0.978 AIC 6 594.786 / -4 377.77 -2 828.084 AICc 6 596.862 / / -2 634.646 RSS 1 354.048 34 631.330 942.141 65.866 残差Moran's I 0.783 / 0.721 0.007 1) OLS和MGWR为指筛选变量后所有样本量的模型计算结果; RF为按7:3划分的测试集的模型计算结果; “/”表示无相关数据

表 5 模型效果评估1) Table 5 Model effectiveness evaluation

MGWR模型中不同带宽可以揭示不同土地利用/覆盖转换方式作用于PM_2.5浓度变化的尺度效应.具体来说, 较大的带宽对应的土地利用/覆盖转换方式对PM_2.5浓度变化的作用方式及强度在较大的空间尺度范围内大致相同, 存在较小的空间梯度, 其空间关系距离衰减程度较弱, 空间异质性程度低, 空间关系趋近平稳, 对空间位置较不敏感.结果表明, 不同土地利用/覆盖转换方式对PM_2.5浓度变化影响的尺度效应有所不同(表 6), 呈现出: V22 < V11 < V12 < V15 < V21 < V25 < V23 < V55 < V51的顺序.未发生土地利用/覆盖转换的类型上, 耕地与林地影响空间尺度明显比建设用地小得多, 距离衰减程度较大, 空间异质性程度更强, 表明PM_2.5浓度变化量对耕地与林地的位置十分敏感.而发生土地类别变换的土地转换类型与未发生土地类别变换相比, 其与PM_2.5浓度变化量的空间关系在作用尺度上呈现扩张状态, 空间关系距离衰减梯度减缓.另外, 空间尺度呈现出耕地向其他用地的转换 < 林地向其他用地的转换 < 建设用地的转换的现象, 耕地转换的空间影响明显受距离所限制, 需要根据其作用方向考虑耕地转换的空间位置及密集度.发生类型转换与转换前类型的带宽相比较来看, 耕地转出类型带宽大小接近耕地-耕地, 建设用地转出类型带宽大小接近建设用地-建设用地, 而相反, 林地转出类型带宽大小与林地-林地差异较大, 那可以推测耕地与建设用地的转换受转换前土地类型的空间尺度效应所影响, 林地的转换对PM_2.5浓度变化量作用空间尺度变化较大, 其影响机制可能更为复杂.

表 6 (Table 6)

表 6 MGWR模型中各变量的空间带宽1) Table 6 Spatial bandwidth of each variable in the MGWR model 变量 V11 V12 V15 V21 V22 V23 V25 V51 V55 带宽/个 57 70 119 295 44 476 423 1 222 1 093 1) 带宽实际意义是纳入局部回归的样本个数, 类型编号中首位数字表示转换前地类代码, 末位数字表示转换后地类代码, 1为耕地, 2为林地, 3为草地, 4为水体, 5为建设用地

表 6 MGWR模型中各变量的空间带宽1) Table 6 Spatial bandwidth of each variable in the MGWR model

MGWR模型输出结果中以每个样本点的回归系数来揭示各个变量对PM_2.5浓度变化的影响的空间异质性特征(图 4).为寻找空间关系与PM_2.5浓度时空变化的关联性, 利用q统计量量化时空聚类分区能否表征空间关系的空间分异程度.结果显示耕地-建设用地(q=0.437)、林地-耕地(q=0.614)、林地-林地(q=0.533)和建设用地-耕地(q=0.833)对PM_2.5浓度变化的空间响应关系与PM_2.5浓度时空演变存在较强的关联性, 那么PM_2.5时空模式的地理分区同样适用于耕地-建设用地、林地-耕地、林地-林地和建设用地-耕地这4个土地转换方式进行地区差别管理.

图 4

Fig. 4

ma为绝对强度平均值; 各土地利用/覆盖转换类型通过95%置信检验的回归系数(即P < 0.05)占全样本量的49.6%~100%; 图中空白区域表明回归系数不具备显著性 图 4 不同土地转换类型与PM2.5浓度变化量回归系数空间分布 Fig. 4 Spatial distribution of regression coefficients between different land conversion types and the change in PM2.5 concentration

从系数的空间分布上看, 不同土地利用/覆盖转换方式对PM_2.5影响方式存在显著差异.以因子系数绝对值的平均值(ma)判定各变量对PM_2.5浓度变化的影响强度, 发现各土地利用/覆盖转换类型对PM_2.5浓度变化的影响强度依次为: V22>V12>V11>V15>V21>V51>V25>V23>V55, 与相关系数强度排序相比可明显发现耕地-林地、耕地-建设用地和林地-建设用地的影响强度增强, 表明人工用地与自然用地的相互转换对PM_2.5的影响需要引起重视.另外从7种土地转换方式对PM_2.5浓度变化的作用信息来看(见表 7), 不同土地转换对PM_2.5浓度变化驱动方式不同, 耕地-耕地、耕地-林地与林地-林地在线性模型和非线性模型中均对PM_2.5变化的作用表现较强, 说明这3类土地转换方式在整个土地转换与PM_2.5浓度变化交互系统中是关键性的地类, 尤其要发挥林地-林地的强负效应及强负空间溢出效应.而林地-耕地与建设用地-耕地对PM_2.5浓度变化的多种作用方式下在全局作用上表现相对明显, 即使两者带宽差异明显, 但仍然要从全区域性宏观调控入手更为有效.

表 7 (Table 7)

表 7 7种土地转换方式对PM2.5浓度变化量的作用信息对比1) Table 7 Comparison of information on the effect of seven land conversions on the change in PM2.5 concentration 变量 MGWR模型 相关系数 随机森林/% SDM模型 V11 [-0.470, 0.210]; 0.111 0.650 29.65 0.012 V12 [-0.304, 0.184]; 0.138 -0.602 17.84 -0.014 V21 [-0.157, -0.021]; 0.078 -0.581 14.13 -0.007(0.017) V22 [-0.813, -0.089]; 0.457 -0.726 17.84 0.043(-0.058) V23 [-0.075, -0.016]; 0.04 -0.336 7.33 0.007(-0.008) V25 [-0.085, 0.064]; 0.051 -0.257 9.37 -0.005(0.019) V51 [0.022, 0.132]; 0.064 0.519 22.57 0.010(-0.02) 1) 类型编号中首位数字表示转换前地类代码, 末位数字表示转换后地类代码, 1为耕地, 2为林地, 3为草地, 4为水体, 5为建设用地; MGWR模型、相关系数、随机森林和SDM模型这4列分别表示4类不同模型结果中各变量与因变量的关系数值表示, 其中MGWR模型括号内表示系数范围, 括号外表示绝对强度平均值, SDM模型括号内为空间滞后变量系数

表 7 7种土地转换方式对PM2.5浓度变化量的作用信息对比1) Table 7 Comparison of information on the effect of seven land conversions on the change in PM2.5 concentration

耕地未发生转换或转出对PM_2.5的影响在区域范围内均存在正负向双重作用[图 4(a)~4(c)].苏沪沿海及浙南地区呈现出耕地本身的正效应[图 4(a)], 而耕地向林地的转换在长三角地区大范围呈现较强的负效应, 尤其在安徽中北部及江苏东部沿海地区[图 4(b)], 耕地向建设用地转换在江苏、安徽大部分地区的回归系数不显著, 一定程度上反映出建设用地的低解释性[图 4(c)].耕地向自然用地与人工用地的转换对PM_2.5浓度的影响出现了较大的差异, 对于其差异性仍然需要根据其效应关系区别差异化管理.多项研究表示耕地对PM_2.5的影响由于农作物的植被效应及秸秆燃烧、农药施肥等农业活动引发颗粒源贡献而存在双向性^{[21, 22]}, 如耕地变化过程中秸秆燃烧、农药施肥等农业活动的变化也可通过与大气的交互及大气对流引发的远距离传输等一系列反馈机制^{[19, 21, 22]}及空间溢出效应而影响区域PM_2.5浓度变化, 其空间影响距离是对耕地及其转换合理调控的重要考虑因素.MGWR模型结果显示耕地及耕地转换的带宽为57~119个, 空间作用半径尺度约为42~61 km, 而黑龙江农作物秸秆燃烧对PM_2.5的影响区域集中在50~80 km之间^[35].耕地上生物质燃烧、土壤翻种等农业活动被认为对PM_2.5具备重要贡献, 但由于农作物的植被效应的季节性变化使耕地成为了动态的源-汇景观^{[21, 22]}.长三角地区禁止秸秆燃烧的政策实施, 2016年后生物质燃烧火点下降60%^[36], 其颗粒源的远距离传输受到一定阻断.此外, 耕地向林地的转换使得植被类型的变化及土壤翻种等农业活动大大减少, 同时所引发的暴雨、致冷效应等间接气候效应使气象条件变化来影响其对PM_2.5浓度的空间影响距离^[37].因此, 土地转换情况下空间影响尺度的收缩意味着需尤为注意对PM_2.5浓度正效应发生区域的耕地转换中心.

林地未发生转换或转出对PM_2.5的影响在区域范围内呈现出负向作用[图 4(d)~4(g)], 林地转换过程中植被类型的变化、植被配置结构以及林地的空间溢出效应也会由于一系列反馈机制使得对PM_2.5作用存在转换类型间的差异.多项研究表示不同林地配置结构因郁闭度、通透度的不同而对PM_2.5的调控能力存在较大差异^[38], 林地-林地对PM_2.5浓度变化量的负效应强度呈现由南向北递增的空间趋势[图 4(e)], 这与林地的空间布局完全相反, 表示绿地规模并不是对PM_2.5浓度呈单调作用, 植物生长也与PM_2.5排放周期的影响绿地滞留效果与需求之间存在错位^[39], 因此仅考虑林地规模对PM_2.5浓度进行调控是不可行的.从图 4(d)~4(e)比较来看, 江苏和安徽地区对比当地林地-林地的负效应, 林地向耕地的转换的负效应呈断崖式减弱, 而在林地广布耕地零散的浙江地区反而是林地-耕地的强负效应区域.林地转换过程中乔木等植被向农作物的植被更换使PM_2.5的调控能力减弱, 另外Cai等^[20]的研究指出表示绿色空间与非绿色空间的相互作用使得非绿色空间与PM_2.5之间的关系随着规模的增加大大降低, 浙江的林地-耕地与较大规模的林地相互作用而弱化了耕地的关系, 同时江苏和安徽地区林地-耕地周围的耕地的空间溢出效应所影响本地PM_2.5浓度的上升使负效应的大幅度减弱.因此, 相同转换类型对PM_2.5的影响会引发空间溢出效应而发生变化.林地向草地的转换则会弱化林地的负效应[图 4(f)], 更是突出植被类型对PM_2.5浓度的影响差异性.而林地向建设用地转换类型通过显著性检验的回归系数在空间上集中在浙江区域, 其强度呈东向西递减的空间趋势[图 4(g)], 呈现以杭州湾为中心的辐射作用.

建设用地未发生转换和向耕地转换对PM_2.5的影响在研究区内呈现出正向作用, 并从西南向东北方向递减, 且空间异质性程度较低, 加上其平均强度较弱, 表现出建设用地的大作用范围及弱正向作用[图 4(g)~4(i)].建设用地作为城市活动的表征, 加上其下垫面易升温的热力性质^[40], 建设用地对PM_2.5的重要贡献与城市内部较弱的扩散条件均在各级城市尺度研究中得到证实^[17~19].但是仅仅考虑土地利用/覆盖转换来看, 未发生转换的建设用地本身对PM_2.5变化贡献并不明显, 说明城市中心区域不是PM_2.5污染控制的关键区域.而城市中心外围所发生的建设用地的转入与转出的影响更强, 尤其是建设用地与耕地之间人类活动的相互触发.其中, 建设用地转向耕地表现出全区域的显著性作用, 作用空间尺度最大, 距离衰减梯度小, 城市联合调控较为合适.

(1) 2000~2018年长三角地区PM_2.5浓度出现了4类空间连续聚集的时空变化模式, 形成南-中-西北-东北空间分异形态, 区域性的同步变化十分强烈.PM_2.5时空模式的地理分区适用于耕地-建设用地、林地-耕地、林地-林地和建设用地-耕地这4个土地转换方式进行地区差别管理.

(2) 土地转换对PM_2.5浓度的相对影响表现复杂, 耕地与林地的源-汇效应明显.周围聚集性土地利用/覆盖转换普遍比单一像元时对PM_2.5浓度变化的作用表现更显著, 空间效应明显.

(3) PM_2.5浓度变化与林地、草地转换类型大多呈显著负相关, 而与耕地、建设用地和水体之间的转换类型呈显著的正相关, 耕地-耕地、耕地-林地、耕地-建设用地、林地-林地和建设用地-耕地之间转换对PM_2.5浓度变化贡献度较高.PM_2.5浓度变化存在显著的空间依赖性, 林地-林地、林地-耕地、林地-草地、林地-建设用地和建设用地-耕地存在较强的空间溢出效应.

(4) MGWR模型在拟合优度及解决空间自相关能力上具备优势, 并揭示了不同土地利用/覆盖转换类型对PM_2.5浓度变化的多尺度空间作用关系.其空间关系表现出强烈的转移类型差异, 而多模型结果显示不同土地转换对PM_2.5浓度变化的驱动方式不同, 因此需分类别分层制定针对性联合管理策略.