基于BFAST算法的神东矿区地表温度突变检测及影响分析

李园园，王藏姣，雷少刚，郭洋楠

(1.中国矿业大学 国土资源研究所，江苏 徐州 221116;2.矿山生态修复教育部工程研究中心，江苏 徐州 221116;3.神华神东煤炭集团有限责任公司 煤炭技术研究院，陕西 神木 719315)

0 引 言

地表温度是地球表面与大气间进行能量交互的结果，同时也是区域乃至全球陆面过程不可或缺的重要参数。[1-3]IPCC第五次评估报告明确指出，1880—2012年全球平均地表温度升高0.85 ℃[4]。地表温度变化对自然生态系统和人类社会经济系统均能产生重大影响，因此，人们对地表温度变化的关注日益增多[5-6]。

学者们已在地表温度突变时空分布特征、地表温度和土地覆被变化间的关系等方面作了大量研究[7-9]。如刘明春等[10]指出1961—2007年整个石羊河流域地表温度突变时间为1996年，但未能定量描述其突变幅度；王雪姣等[9]通过Mann-Kendall突变检验，指出1961—2015年间新疆地表温度呈显著波动上升趋势，地表温度突变存在明显的季节差异，但未对地表温度发生突变原因进行探究；M.Bokaie等[11]基于Landsat TM影像研究发现,地表温度与土地覆被呈负相关关系，但未能定量描述地表温度与土地覆被变化间的关系；T.D.Mushore等[12]研究指出，在津巴布韦土地覆被变化对城市温度升高的贡献大于0.5 ℃，但未对各种土地覆被转变与地表温度变化之间的关系进行解释。上述研究主要是针对地表温度突变时空分布特征、地表温度与土地覆被变化间的定性关系进行分析，具有较大的学术价值。相对于其他自然地理区域，在大规模高强度的煤炭开采及生态重建等多种因素共同作用下，矿区内地表水热特性存在巨大差异，地表温度变化也更为复杂。针对矿区大面积自然下垫面，开展地表温度突变及其与土地覆被变化间的定量关系研究具有一定学术价值。

基于可靠的地表温度数据和检测方法，准确检测出矿区地表温度突变是开展以上研究的前提。地表温度的获取分为传统定点观测和遥感监测两种，相对于传统定点观测范围小、受地理条件限制等不足，热红外遥感技术凭借其大范围且瞬时成像的优势，其时间序列已成功用于量化地表温度动态变化[13]。在中等分辨率的卫星传感器(如Landsat)不易获取足够多的长时间序列无云图像情况下，中分辨率成像光谱仪(MODIS)因其高时间分辨率、适中的空间分辨率成为更优选择[14]。目前用于地表温度突变检测的方法主要包括Mann-Kendall法、滑动t检验法、BFAST算法等[15-17]。Mann-Kendall法通常用于地表温度年数据检测，且易受季节性影响[18]；滑动t检验法检测出的突变点位置会随子序列长度的变化而发生飘移[16]；BFAST算法能够克服以上两种方法的缺陷，基于季节-趋势分解理论，分解时间序列，通过最小二乘移动求和(OLS-MOSUM)，判定时间序列中突变点位置，利用贝叶斯信息准则确定最佳突变点数量，可有效处理季节性数据，广泛应用于植被、气候、水文等领域时间序列的突变检测[19]。

基于以上研究，本文采用BFAST算法，以神东矿区为研究区，选取2000—2018年共867景MODIS地表温度数据MOD11A2构建地表温度时间序列，拟合其时序变化趋势，检测地表温度突变，通过空间叠加统计分析地表温度突变与土地覆被变化间的关系。具体研究内容为：(1)获取2000—2018年神东矿区地表温度长期变化趋势；(2)提取矿区地表温度突变时空分布特征和变化幅度；(3)分析半干旱矿区地表温度突变与土地覆被变化之间的关系，以期为矿区生态环境治理提供有效的决策信息。

1 材料与方法

1.1 研究区概况

神东矿区位于晋陕蒙三省交界处，地理坐标为经度109°51′17″～110°49′05″ E，纬度38°52′32″～39°39′09″ N，地处毛乌素沙地的东南边缘，总面积为3 537 km2(图1)。研究区属于典型的干旱、半干旱大陆性气候，干燥少雨，且降雨主要集中在6—9月，年均降水量360 mm左右，年均日照时数2 875.9 h，年均气温7.29 ℃，年均蒸发量2 000 mm左右，年均湿度56%。土壤类型多为风沙土，结构松散，透水性强，保水保肥能力差，土壤贫瘠。矿区地处干草原和森林草原的过渡地带，主要植被类型为干草原、落叶阔叶灌丛和沙生类型植被。

图1 研究区位置示意图

1.2 数据源及处理

为构建地表温度时间序列，获取 2000-01-01—2018-12-31共867景MODIS/Terra卫星MOD11A2地表温度数据。数据来源于美国航空航天局(https://www.nasa.gov)，该数据产品依据分裂窗算法计算地表温度，数据级别为L3级，投影方式为Sinusoidal投影，时间分辨率为8 d，空间分辨率为1 km。运用MRT(MODIS reprojection tool)将MOD11A2转换至UTM投影，WGS84-49N坐标，沿研究区边界进行裁剪。基于timesat3.0平台，选用Savitsky-Golay方式对裁剪后的影像进行滤波，从而剔除影像中因外强迫和仪器测量误差等引起的噪声和扰动，避免此类虚假信息造成误判突变点[16]。然后，通过最大值合成将其转换为16 d最大值合成数据，并将DN值转换为地表温度[20]。

为进一步分析地表温度突变与土地覆被变化的定量关系，对突变显著的2002-06-08的Landsat 7 ETM+和2003-08-17 的Landsat 5 TM 共2景影像利用随机森林分类算法进行土地利用分类，结合土地利用分类体系以及研究区特点，将研究区土地利用类型分为耕地、林地、草地、建设用地、沙地及水域6类。分类影像下载于美国航空航天局(https://www.nasa.gov)，该影像为L1T级数据，空间分辨率为30 m，影像质量良好。L1T级数据已经过系统几何校正，达亚像元精度，对其进行辐射校正、大气校正、裁剪等预处理工作[21]。两景影像分类总体精度均达到90%以上，Kappa系数达到0.86以上。

1.3 研究方法

1.3.1 趋势分析

为分析神东矿区地表温度时间序列长期变化趋势，运用最小二乘法逐像元拟合地表温度年际变化趋势，拟合斜率即为地表温度随时间变化的线性回归系数，公式为

(1)

式中：Slope为拟合斜率；n为时间序列长度；LSTi为第i年的LST值。

通过相关系数r判定变化趋势的显著程度，采用t检验法对相关系数r的显著性水平进行检验，公式为

(2)

(3)

式中，t的检验自由度为n-2。

1.3.2 BFAST算法

BFAST算法最早由J.Verbesselt等[22]于2010年提出，该算法通过迭代时间序列的方法将原始时间序列分解为季节组分、趋势组分和残差组分，能够有效检测半干旱地区遥感影像时间序列中发生的季节与趋势突变[23]。理论模型为

Yt′=Tt′+St′+et′，(t′=1,2,…,n),

(4)

式中：Yt′为观测数据；Tt′为趋势组分；St′为季节组分；et′为残差组分；t′为观测时间。

Tt′=αi+βit′,

(5)

式中：i=1，…，m,m为突变点的总数；αi和βi分别为线性模型的截距和斜率。

(6)

式中：j为突变点所在位置，取值从1到q，q为突变点总数；k为调和项数目；αj,k为振幅；δj,k为相位；f为频率(本研究采用的地表温度遥感影像合成间隔为16 d，f为23)。

季节组分和趋势组分中突变点的识别需要确定突变点的数量和突变点在时间序列中的位置。此处利用最小二乘移动求和从季节组分和趋势组分中检验判断是否存在突变点。季节组分和趋势组分中突变点数量和位置的计算：首先分别剔除季节组分(Yt′-St′)和趋势组分(Yt′-Tt′)，其次采用贝叶斯信息准则确定突变点的最优个数，最后通过最小二乘法估计突变点在时间序列中的位置。本文BFAST算法由R语言的bfast包实现(https://cran.r-project.org/web/packages/bfast/index.html)。

图2 技术路线图

2 研究结果

2.1 地表温度长期变化趋势

2000—2018年，神东矿区98.63%的区域地表温度呈下降趋势，1.37%的区域地表温度呈上升趋势(表1)，但变化趋势均未达到显著性水平(p>0.05,-0.456

表1 地表温度长期变化趋势表

2.2 地表温度突变时空分布特征

本文重点分析地表温度长期变化趋势中发生的突变，地表温度趋势组分突变的时空分布特征分为两个方面：(1)地表温度突变频率的时空分布特征。统计不同年份不同区域的地表温度突变分布。(2)地表温度最大突变发生的时间和幅度。提取地表温度发生一次或多次突变的区域中最大突变发生时间和幅度。

2.2.1 突变频率的时空分布特征

2000—2018年，神东矿区地表温度突变次数为0～4次，其中41.56%的区域地表温度发生了至少1次突变，发生1～4次突变的区域分别占研究区面积的32.88%，7.86%，0.79%和0.03%(图3(a))。井工采区内地表温度突变面积的比例远低于整体水平，仅19.02%的区域地表温度发生突变，其中13.41%的区域地表温度发生了1次突变，5.59%的区域地表温度发生了2次突变，0.01%的区域内地表温度发生了3次突变。露天采区内地表温度突变面积的比例高于井工采区，露天采区内34.66%的区域地表温度发生突变，其中地表温度发生了1次突变的区域面积占比达到32.16%，发生2次突变的区域面积占比为2.47%，发生了3次突变的区域面积占比为0.03%。

神东矿区地表温度突变时间主要集中在2001—2005年和2014年，期间地表温度突变总面积占神东矿区面积的44.81%，其中2003年地表温度突变面积占比最高，高达20.07%(图3(b))。2006—2018年，地表温度发生突变面积大幅减少，呈现稳定变化态势。井工采区地表温度在2002年发生突变面积最多，占井工采区面积的7.90%；露天采区地表温度2014年发生突变面积最多，占露天采区面积的16.25%。

图3 神东矿区地表温度突变时空分布图(2000—2018年)

2.2.2 最大突变发生时间和突变幅度

从地表温度发生最大突变时间看，神东矿区2000，2006，2015，2016，2017和2018年共6年地表温度未出现最大突变，最大突变出现的年份主要为2003和2014年，占发生突变区域面积的比例分别达到46.53%和13.38%(图4(a))。井工采区内地表温度发生最大突变最多的年份为2002年，面积达到了65 km2。露天采区内地表温度发生最大突变最多的年份为2014年，面积达到了16 km2。

从地表温度发生最大突变幅度看，神东矿区地表温度发生最大突变幅度为-6.82～7.16 ℃，其中20.43%区域地表温度负向突变，79.57%区域地表温度正向突变(图4(b))。井工采区地表温度发生最大突变幅度为-6.80～6.48 ℃，其中30.22%区域地表温度正向突变，69.78%区域地表温度负向突变。露天采区地表温度最大突变幅度为-6.55～5.10 ℃，其中80%区域地表温度正向突变，20%区域地表温度负向突变。

图4 神东矿区地表温度最大突变发生时间和幅度分布图(2000—2018年)

2.3 地表温度突变与土地覆被变化间关系

本研究选取2003年地表温度发生最大突变的区域分析地表温度突变与土地覆被变化间关系。基于ENVI变化检测方法获取2003年地表温度发生最大突变区域土地覆被变化信息,通过空间叠加分析，统计得出土地覆被变化与地表温度突变间的关系。从地表温度突变幅度看，地表温度突变均表现为正向突变，变化幅度在1.33～6.93 ℃间(图5(a))。从土地覆被变化看，建设用地、沙地以及林地面积均有所增加，分别增加了15，5，3 km2；草地、耕地和水域面积均减少，分别耕地转为沙地，地表温度突变幅度为1.94～2.49 ℃；草地转为沙地，地表温度突变幅度为1.92～2.42 ℃。水域转为草地，地表温度突变幅度为2.21 ℃；林地转为草地，地表温度突变幅度为1.92～2.62 ℃；耕地转为草地，地表温度突变幅度为1.74～2.76 ℃。草地转为林地，地表温度突变幅度为1.84～2.82 ℃。林地转为耕地，地表温度突变幅度为2.32～2.61 ℃；草地转为耕地，地表温度突变幅度为1.33～2.79 ℃。草地转为建设用地，地表温度突变幅度最大，达1.89～6.93 ℃。

图5 2003年地表温度突变与土地覆被变化间的关系

3 讨 论

植被在调节地表温度方面扮演着重要角色，如植物叶片的蒸腾作用，持续吸收热辐射，能够有效调节湿度，降低温度。[25]李宏宇等[26]研究表明，在中纬度地区，地表温度随植被覆盖度的增加而降低。近年来，神东矿区的植被覆盖状况得到了较大改善，矿区整体植被覆盖度由矿区开采初期的3%～11%提高到59%以上[27]，其中马家塔煤矿的植被覆盖率达到80%[28]。整体植被覆盖度的提高对地表产生的降温效果进一步证明了19年来神东矿区地表温度长期呈下降趋势的合理性。

此外，就不同开采方式而言，露天开采比井工开采更易发生地表温度突变。可归结为两种开采方式对地表扰动程度的不同。神东矿区井工开采除变形较大的区域外，采区与非采区并无明显差异[29,30]。比如，沉陷区域的中性区随着工作面推过，地表会逐渐稳定，开采过程中形成的临时性地裂缝能够自动愈合。[31]相对于井工开采，露天开采对地表的破坏更为明显和直接，露天开采过程如岩土剥离、采掘场开挖、排土场堆积等均对地表形态造成较大破坏，地表受到大幅度扰动后，露天矿区产生明显的增温效应[32]，本研究中露天采区地表温度发生最大突变的区域80%为正向突变。

地表温度突变方向取决于土地覆被变化类型和规模，梁保平等[33]指出建设用地和沙地面积增加会导致地表温度升高，原因为沙地和建设用地的地表潜热容量小，接收太阳辐射之后，地表温度迅速升高。草地、耕地面积减少同样能引起地表温度的升高，原因为植被能够吸收太阳辐射并进行蒸腾作用，可减少到达地表的热量。水域比热容大，其面积减少会直接造成地表温度升高。林地面积增加虽然能够有效降低地表温度，但林地面积增加很小时，冷却效应并不明显。[34]因此，2003年地表温度发生最大突变的区域受建设用地、沙地面积增加，耕地、草地和水域面积减少影响较大，而受林地面积增加影响较小，地表温度均表现为正向突变。同样地，地表温度突变幅度同样取决于土地覆被变化，不同土地覆被的地表温度存在较大差异[35]，因而不同土地覆被变化引起地表温度突变的幅度也不尽相同。李润林等[36]研究表明在甘肃省张掖市，耕地、草地、林地、水域4种土地覆被与沙地的地表温度差值为1.20～4.90 ℃。本研究中耕地、草地转为沙地地表温度突变幅度均为1.92～2.49 ℃。潘竟虎等[37]研究表明，甘肃省张掖市水域、耕地和建设用地3种土地覆被与草地的地表温度差值为1.66～7.42 ℃，本研究中水域、林地、耕地转为草地时地表温度突变幅度为1.74～2.76 ℃。本文研究区与张掖市地理位置、气候条件等具有可比性，因此认为本文研究结果合理。

4 结 论

(1)2000—2018年，在我国西北地区向暖湿化转型的背景下，神东矿区整体地表温度呈下降趋势，可归因于矿区整体植被覆盖的大幅改善，表明生态修复对区域气候改善具有重要意义。

(2)露天采区比井工采区更容易发生地表温度突变，露天采区内地表温度突变面积占比为34.66%，而井工采区内仅19.02%的区域地表温度发生突变，原因为露天开采对地表形态扰动程度远高于井工开采。露天采区地表温度发生最大突变的区域中80%为温度正向突变，表明露天开采具有明显的增温效应。

(3)土地覆被变化引起地表温度发生突变，但地表温度突变方向和幅度取决于土地覆被变化类型、规模等因素。建设用地、沙地面积增加，耕地、草地、水域面积减少可导致地表温度发生正向突变。在土地覆被转变类型中，草地转为建设用地，地表温度突变幅度最大，为1.89～6.93 ℃。本研究采用中分辨率地表温度数据产品，下一步可选用分辨率更高的遥感影像，从而对地表温度突变进行更精细的检测。此外，可将地表温度突变检测与土地覆被连续变化监测相结合，以研究地表温度突变与土地覆被变化驱动间的长期变化规律。