Analysis of spatiotemporal evolution characteristics and driving factors of land desertification in Beijing from 1990 to 2020

研究区位于华北大平原的西北边缘，39°28′~41°05′N， 115°25′~117°35′E，总面积为1.64×10⁴ km²，地势西北高、东南低，西部、北部和东北部三面环山，东南部是一片缓缓向渤海倾斜的平原。境内流经的主要河流有永定河、潮白河、拒马河等。北京市的气候为暖温带半湿润半干旱季风气候，特点为春秋短，夏季高温多雨，冬季寒冷干燥。

北京市风沙灾害多发于3—5月。谢宜嘉等^[15]的研究结果表明：北京市植被的春季物候平均起始于第117天(即4月下旬)， 5月植被已处于生长阶段。为便于计算NDVI，基于Google Earth Engine云计算平台，本研究选用1990—2020年5月的Landsat系列遥感影像数据，分辨率为30 m×30 m。用ENVI 5.3软件对获取的遥感影像进行几何精校正、辐射校正、影像裁剪等预处理，得到北京市土地沙化遥感影像。此外，本研究为确保所有数据(表1)与遥感影像分辨率保持一致，利用ArcGIS 10.8软件对所有数据进行预处理，通过投影坐标、按掩膜提取、创建渔网、多值提取至点以及普通克里金插值等多项技术操作，将数据分辨率均调整为30 m×30 m，以便后续分析。

为定量表征地表植被覆盖与能量反射特征，本研究采用归一化植被指数(NDVI)与地表反照率(Albedo)等2个关键指标，构建Albedo-NDVI二维特征空间。其中，NDVI反映植被绿度与生物量状况^[18]，Albedo表征地表对太阳辐射的反射能力^[19]。同时为保障2个指标量纲统一并增强可比性，分别对NDVI和Albedo进行归一化处理，得到N_norm与A_norm。NDVI、Albedo、N_norm、A_norm参考文献[20]计算。

在Albedo-NDVI特征空间中，Albedo与NDVI之间存在显著的线性关系。本研究以A_norm为因变量、N_norm为自变量构建线性回归模型。基于该模型，可进一步定义荒漠化差值指数(I_DDI)，用以反映研究区内沙化状况。线性回归模型与I_DDI计算参考文献[21]。

在Albedo-NDVI特征空间中，不同地表类型沿特定轨迹分布，其中沙化土地的变化路径与特征空间中的某一边界(图1中AC边)接近^[21]。基于此，可通过垂直于该边的直线实现对沙化与非沙化土地的区分。本研究依据北京市第6次荒漠化与沙化监测结果及土壤类型分布，利用自然断点法^[20]对I_DDI进行分级，将其划分为沙化土地与未沙化土地2个类别，完成土地沙化状况的信息提取。

Figure 1.  Schematic illustrating the Albedo-NDVI spatial relationship

验证分类精度时，可信度高的方法是混淆矩阵与Kappa系数协同验证^[22]。这种组合方法不仅可以精细分类模型的识别能力，还稳健地提升了总体评价可信度。混淆矩阵采用构建参考数据与分类结果交叉对比矩阵的方式(采用总体精度、用户精度、生产者精度3个指标)来量化分类模型的预测性能。该矩阵采用n×n结构，其中行表示参考数据的真实类别，列表示模型的预测结果。Kappa系数可通过消除偶然一致性的影响客观反映模型分类结果与真实情况的吻合程度。总体精度、用户精度、生产者精度、Kappa系数参考文献[18]计算。

土地利用转移矩阵可直观反映研究期内不同土地利用类型的沙化加重与改善情况，其计算参考文献[23]。

单一动态度是一种用于估算沙化土地面积变化速度的方法，可定量分析研究期内各土地利用类型中沙化土地的面积变化强度，其计算参考文献[24]。

地理探测器是用来研究空间的分异规律及其影响因子的统计方法^[25]。本研究运用其因子探测器与交互探测器模块，探究了北京市土地沙化空间异质性及各因素驱动强度。①因子探测器。选取植被覆盖、气象因素及人类活动3类共7种驱动因素(包括植被指数、年降水量、年平均气温、年平均风速、地表蒸散发、人口密度与地区生产总值)，定量分析其对沙化土地空间分布的影响程度。各因素对沙化土地影响程度的解释力用统计量q表示(计算参考文献[25])，q值越大，影响越大。②交互探测器。探究7种驱动因素之间的交互作用，明确其对土地沙化的影响效应。交互类型见表2。③最优参数选择。基于地理探测器输入数据为类别型变量，采用R4.2.1软件的GD包^[26]，针对每个连续型驱动因子，分别测试等间隔、自然断点、分位数等5种离散化方法在3~10个分类区间下的表现^[27−28]，并选取能产生最大q值的组合作为该因子的最优离散方案，以确保探测结果的可靠性。

本研究在各期影像中随机选取200个样点，以NDVI为自变量(X)，Albedo为因变量(Y)构建特征空间散点图并得到其拟合方程。由图2可知：北京市1990—2020年Albedo-NDVI呈较强负相关，各期决定系数(R²)均高于0.68。

Figure 2.  Albedo−NDVI feature space in 1990－2020

根据Albedo-NDVI拟合方程得到1990—2020年荒漠化差值指数(I_DDI)计算公式及经自然断点法分级后，I_DDI与沙化土地关系见表3。

本研究基于每期影像随机布设的300个样本点(均匀覆盖各土地类型)，结合Landsat真彩色影像与Google Earth高分辨率地图进行目视解译，完成荒漠化差值指数模型精度评价。由表4可知：模型总体精度为81.00%~90.67%，用户精度与生产者精度分别为71.14%~92.90%和77.40%~92.09%。尽管部分年份如1990和2000年未沙化土地用户精度、2000年沙化土地生产者精度较低，但整体分类精度良好。此外，7期Kappa系数为0.619~0.806，其中2005年系数值最高，较最低的2000年提升约30%。参照Kappa系数分级标准，所有系数值均高于0.600，模型分类结果一致性较高，满足研究要求。

研究结果表明：1990—2020年期间北京市沙化土地呈现西北—东南的地理分布格局，集中分布在平原与沿河两岸。在平原地区，沙化土地集中分布在延庆的康庄地区、昌平的南口地区，怀柔和密云的大沙河地区，在山前风沙区和房山拒马河地区也有零星分布。在沿河两岸，沙化土地主要分布于潮白河、永定河、温榆河沿岸。以行政区划来看，1990—2020年北京市沙化土地分布于延庆区、怀柔区、密云区、平谷区、昌平区、顺义区、朝阳区、丰台区、房山区、大兴区、通州区共11个区。

由图3可知：2010年前北京及各区沙化土地面积呈缓慢下降趋势(仅1995年略增)，总面积由1990年的5.54×10⁴ hm²降至2010年的5.24×10⁴ hm²，降幅约3%~4%。2010—2020年间则骤减，至2020年沙化土地面积仅为2.23×10⁴ hm²，其中2010—2015年降幅最大，达47.46%，2015—2020年降幅减缓至19.23%。从区域分布看，1990—2020年，大兴区沙化面积占比由32.14%减少至24.52%，位于第1位；延庆区占比由18.45%增至25.15%，成为第二大分布区；朝阳区、丰台区占比最小，均不足1.00%。整体来看，近30 a间沙化土地面积缩减明显，尤其在2010年后治理效果尤为突出。

Figure 3.  Summary of desertified land area in Beijing from 1990 to 2020

为深入解析北京市1990—2020年沙化土地的来源类型，更好地量化沙化土地中不同土地利用类型的面积变化强度，本研究利用土地利用转移矩阵和单一动态度2种方法计算得出不同土地利用类型中沙化土地的面积(表5)及动态度(表6)。

由表5可知：1990—2020年间北京市不同土地利用类型的沙化面积总体呈下降趋势。7期数据中，林地沙化面积占比均最高，但整体持续减少；未利用地占比最低，仅1995年与2020年有轻微沙化外，其余年份无沙化。

由表6可知：各土地利用类型在不同时期的沙化变化差异明显，沙化重心持续转移。其中1990—1995年水域动态度最高(15.14%)，1995—2000年未利用地动态度最明显(−20.00%)，2000—2005年建筑等用地动态度最大(3.65%)，2005—2010年耕地动态度最明显(−11.59%)。2010—2015年所有地类动态度均为负值，建筑等用地减少最快(−12.95%)；2015—2020年草地动态度转为最高(17.42%)。整体而言，沙化土地来源受土地利用类型变化影响，具有明显的时段差异与类型更替特征。

本研究选取植被指数(x₁)、年降水量(x₂)、年平均气温(x₃)、年平均风速(x₄)、地表蒸散发(x₅)、人口密度(x₆)与地区生产总值(x₇)共7个自然与人为因素，运用地理探测器中的因子探测与交互探测模型，分别评估各因子的独立解释力及交互作用对土地沙化的影响。

根据7 a数据的平均解释力(q)排序(表7)，各因子对沙化土地的影响权重从高到低依次为人口密度(x₆)、地表蒸散发(x₅)、地区生产总值(x₇)、年平均气温(x₃)、年降水量(x₂)、植被指数(x₁)、年平均风速(x₄)。人口密度是土地沙化的关键驱动因素，地表蒸散发、地区生产总值、年平均气温、年降水量也具有一定影响。植被指数和年平均风速的贡献较弱(q＜0.150)。

同时，各因素对土地沙化的解释力存在显著年际差异。人口密度在多数年份解释力（q）均高于0.350，2010年达到峰值0.428，为最关键驱动要素；地表蒸散发解释力次之，2020年最高，为0.443。地区生产总值、年平均气温与年降水量解释力相近，约0.190，均在2020年达到峰值。而植被指数与风速的解释力长期处于低位且变化平缓，表明其影响较为稳定且有限。

图4为1990—2020年7个驱动因素对土地沙化协同作用的解释力。结果显示：任意两因子的交互作用均明显高于其单独作用，因子间普遍存在协同增强效应。1990、2010、2020年地区生产总值与其余因子交互作用，及2000年、2010年植被指数与其余因子交互作用，均呈双向增强；2015年地区生产总值、植被指数分别与其余因子呈非线性增强；其余因子交互作用则在双向增强与非线性增强间波动。

Figure 4.  Interactive detection results of driving factors of land desertification from 1990 to 2020

从不同年份看，驱动组合随时间发生演变。1990年与2000年，年降水量与人口密度的交互作用最强(q＞0.500)；1995—2010年，地表蒸散发与人口密度的解释力最高(q＞0.600)；2015年后，地表蒸散发与地区生产总值的交互作用跃居主导地位(q＞0.450)。从驱动因素看，1990—2020年地表蒸散发与人口密度的交互作用明显，与各自高贡献度一致。年平均风速除2020年外，与其他因子的交互作用均较弱，其中与植被指数的交互作用在部分年份仅约0.100。

土地沙化是多因子协同驱动的结果，不同阶段自然与人为因素对沙化过程反映出复合的影响机制。本研究基于1990—2020年Landsat遥感影像，构建Albedo-NDVI特征空间与荒漠化差值指数模型，实现了北京市沙化土地的动态识别。与张国桢^[29]、刘克等^[14]采用的专家打分法、遥感影像目视解译相比，本方法在提升提取精度的同时有效降低了人为误差，为城市沙化动态监测提供了可靠技术路径。

北京市沙化土地分布格局受自然基底与外部沙源的共同影响。陈广庭^[30]指出：北京市平原地区长期受潮白河、永定河等河流冲积作用，土壤以砂粒含量高的亚砂土、亚黏土和粗细砂为主，结构松散、黏结力弱，抗侵蚀能力差，为沙化提供了物质基础。同时，北京受周边沙源区(如内蒙古科尔沁沙地、河北宣化盆地)及主要风道(如延庆狼山风口、密云古北口)影响^[31]，尽管年平均风速常年低于2.6 m·s⁻¹，风蚀贡献有限，但外源沙尘输送仍是不可忽视的驱动因素。

自20世纪90年代以来，随着“三北”防护林、京津风沙源治理、平原造林^[32]等重点治沙工程持续推进，以及两轮百万亩造林绿化工程的实施，城市公园、口袋公园兴建，北京市完成退化草原治理2.61×10⁷ hm²，新增森林面积3.28×10⁶ hm²，植被覆盖度显著提升^[33]，沙化状况显著好转。从土地利用类型看，林地沙化土地占比最高，主要因其组成以固定及半固定沙地为主，经过长时间的发展各方面都趋于稳定^[29]，长期内不会发生显著改善。而耕地因春季植被覆盖度低，NDVI值波动小，导致其在驱动分析中解释力较弱。因此，生态修复措施对沙化逆转具有主导作用，但遥感时相与植被物候的匹配程度也影响驱动机制的准确识别。

全球气候变暖加剧了北京地区的环境压力，但2010年后实施的生态补水政策^[34]有效改善了区域水分条件，促使2015年地表蒸散发对沙化的影响程度达到峰值。与此同时，城市化快速发展导致人口集聚与建设扩张，大量占用耕地、林地且过度开采的行为加剧了沙化恶化趋势。2014年“非首都功能疏解”政策^[35]实施后，人口密度骤减，人口压力显著缓解。人口密度与地表蒸散发交互作用的解释力多年最高，表明北京市沙化过程受“人为−水分”协同机制调控，水分改善可在一定程度上缓解人为压力对沙化的驱动作用。该机制是2015年沙化面积降至最低的关键原因，也为生态调控与人口疏解相结合的治理策略提供了理论依据。

本研究在揭示北京市沙化时空演变及其驱动机制方面仍存在一定局限：①模型未引入土壤湿度等指标，后期可引入以提升沙化土地的识别精度；②受数据限制，未引入空间计量模型，难以量化沙化格局的稳定性；③NDVI在低覆盖区敏感性有限，后续可融合MODIS-EVI等数据提升动态监测能力。

本研究采用1990—2020年Landsat长时间序列遥感影像，基于Albedo−NDVI特征空间构建荒漠化差值指数模型，实现了北京市高精度土地沙化空间信息提取，再结合地理探测器定量分析了其关键驱动机制。主要结论如下：①30 a间北京市土地沙化集中分布于平原和沿河两岸，沙化面积显著减少3.3×10⁴ hm²，整体呈现显著逆转趋势；②驱动机制分析表明：人口密度是关键驱动因素，且年降水量、地表蒸散发量与人口密度的双因子交互作用具有最强的解释力，揭示了“人类活动-水分条件”的协同作用对城市沙化的主导影响；③未来治理应聚焦沿河植被带和季节性耕地，采取水土保持和植树造林等综合措施，并构建“卫星遥感-地面站点-无人机巡检”三位一体的土地沙化监测网络，以实现精准防治与土地资源可持续利用。