李園園，王藏姣，雷少剛，郭洋楠

(1.中國(guó)礦業(yè)大學(xué) 國(guó)土資源研究所，江蘇 徐州 221116;2.礦山生態(tài)修復(fù)教育部工程研究中心，江蘇 徐州 221116;3.神華神東煤炭集團(tuán)有限責(zé)任公司 煤炭技術(shù)研究院，陜西 神木 719315)

0 引 言

地表溫度是地球表面與大氣間進(jìn)行能量交互的結(jié)果，同時(shí)也是區(qū)域乃至全球陸面過(guò)程不可或缺的重要參數(shù)。[1-3]IPCC第五次評(píng)估報(bào)告明確指出，1880—2012年全球平均地表溫度升高0.85 ℃[4]。地表溫度變化對(duì)自然生態(tài)系統(tǒng)和人類(lèi)社會(huì)經(jīng)濟(jì)系統(tǒng)均能產(chǎn)生重大影響，因此，人們對(duì)地表溫度變化的關(guān)注日益增多[5-6]。

學(xué)者們已在地表溫度突變時(shí)空分布特征、地表溫度和土地覆被變化間的關(guān)系等方面作了大量研究[7-9]。如劉明春等[10]指出1961—2007年整個(gè)石羊河流域地表溫度突變時(shí)間為1996年，但未能定量描述其突變幅度；王雪姣等[9]通過(guò)Mann-Kendall突變檢驗(yàn)，指出1961—2015年間新疆地表溫度呈顯著波動(dòng)上升趨勢(shì)，地表溫度突變存在明顯的季節(jié)差異，但未對(duì)地表溫度發(fā)生突變?cè)蜻M(jìn)行探究；M.Bokaie等[11]基于Landsat TM影像研究發(fā)現(xiàn),地表溫度與土地覆被呈負(fù)相關(guān)關(guān)系，但未能定量描述地表溫度與土地覆被變化間的關(guān)系；T.D.Mushore等[12]研究指出，在津巴布韋土地覆被變化對(duì)城市溫度升高的貢獻(xiàn)大于0.5 ℃，但未對(duì)各種土地覆被轉(zhuǎn)變與地表溫度變化之間的關(guān)系進(jìn)行解釋。上述研究主要是針對(duì)地表溫度突變時(shí)空分布特征、地表溫度與土地覆被變化間的定性關(guān)系進(jìn)行分析，具有較大的學(xué)術(shù)價(jià)值。相對(duì)于其他自然地理區(qū)域，在大規(guī)模高強(qiáng)度的煤炭開(kāi)采及生態(tài)重建等多種因素共同作用下，礦區(qū)內(nèi)地表水熱特性存在巨大差異，地表溫度變化也更為復(fù)雜。針對(duì)礦區(qū)大面積自然下墊面，開(kāi)展地表溫度突變及其與土地覆被變化間的定量關(guān)系研究具有一定學(xué)術(shù)價(jià)值。

基于可靠的地表溫度數(shù)據(jù)和檢測(cè)方法，準(zhǔn)確檢測(cè)出礦區(qū)地表溫度突變是開(kāi)展以上研究的前提。地表溫度的獲取分為傳統(tǒng)定點(diǎn)觀測(cè)和遙感監(jiān)測(cè)兩種，相對(duì)于傳統(tǒng)定點(diǎn)觀測(cè)范圍小、受地理?xiàng)l件限制等不足，熱紅外遙感技術(shù)憑借其大范圍且瞬時(shí)成像的優(yōu)勢(shì)，其時(shí)間序列已成功用于量化地表溫度動(dòng)態(tài)變化[13]。在中等分辨率的衛(wèi)星傳感器(如Landsat)不易獲取足夠多的長(zhǎng)時(shí)間序列無(wú)云圖像情況下，中分辨率成像光譜儀(MODIS)因其高時(shí)間分辨率、適中的空間分辨率成為更優(yōu)選擇[14]。目前用于地表溫度突變檢測(cè)的方法主要包括Mann-Kendall法、滑動(dòng)t檢驗(yàn)法、BFAST算法等[15-17]。Mann-Kendall法通常用于地表溫度年數(shù)據(jù)檢測(cè)，且易受季節(jié)性影響[18]；滑動(dòng)t檢驗(yàn)法檢測(cè)出的突變點(diǎn)位置會(huì)隨子序列長(zhǎng)度的變化而發(fā)生飄移[16]；BFAST算法能夠克服以上兩種方法的缺陷，基于季節(jié)-趨勢(shì)分解理論，分解時(shí)間序列，通過(guò)最小二乘移動(dòng)求和(OLS-MOSUM)，判定時(shí)間序列中突變點(diǎn)位置，利用貝葉斯信息準(zhǔn)則確定最佳突變點(diǎn)數(shù)量，可有效處理季節(jié)性數(shù)據(jù)，廣泛應(yīng)用于植被、氣候、水文等領(lǐng)域時(shí)間序列的突變檢測(cè)[19]。

基于以上研究，本文采用BFAST算法，以神東礦區(qū)為研究區(qū)，選取2000—2018年共867景MODIS地表溫度數(shù)據(jù)MOD11A2構(gòu)建地表溫度時(shí)間序列，擬合其時(shí)序變化趨勢(shì)，檢測(cè)地表溫度突變，通過(guò)空間疊加統(tǒng)計(jì)分析地表溫度突變與土地覆被變化間的關(guān)系。具體研究?jī)?nèi)容為：(1)獲取2000—2018年神東礦區(qū)地表溫度長(zhǎng)期變化趨勢(shì)；(2)提取礦區(qū)地表溫度突變時(shí)空分布特征和變化幅度；(3)分析半干旱礦區(qū)地表溫度突變與土地覆被變化之間的關(guān)系，以期為礦區(qū)生態(tài)環(huán)境治理提供有效的決策信息。

1 材料與方法

1.1 研究區(qū)概況

神東礦區(qū)位于晉陜蒙三省交界處，地理坐標(biāo)為經(jīng)度109°51′17″～110°49′05″ E，緯度38°52′32″～39°39′09″ N，地處毛烏素沙地的東南邊緣，總面積為3 537 km2(圖1)。研究區(qū)屬于典型的干旱、半干旱大陸性氣候，干燥少雨，且降雨主要集中在6—9月，年均降水量360 mm左右，年均日照時(shí)數(shù)2 875.9 h，年均氣溫7.29 ℃，年均蒸發(fā)量2 000 mm左右，年均濕度56%。土壤類(lèi)型多為風(fēng)沙土，結(jié)構(gòu)松散，透水性強(qiáng)，保水保肥能力差，土壤貧瘠。礦區(qū)地處干草原和森林草原的過(guò)渡地帶，主要植被類(lèi)型為干草原、落葉闊葉灌叢和沙生類(lèi)型植被。

圖1 研究區(qū)位置示意圖

1.2 數(shù)據(jù)源及處理

為構(gòu)建地表溫度時(shí)間序列，獲取 2000-01-01—2018-12-31共867景MODIS/Terra衛(wèi)星MOD11A2地表溫度數(shù)據(jù)。數(shù)據(jù)來(lái)源于美國(guó)航空航天局(https://www.nasa.gov)，該數(shù)據(jù)產(chǎn)品依據(jù)分裂窗算法計(jì)算地表溫度，數(shù)據(jù)級(jí)別為L(zhǎng)3級(jí)，投影方式為Sinusoidal投影，時(shí)間分辨率為8 d，空間分辨率為1 km。運(yùn)用MRT(MODIS reprojection tool)將MOD11A2轉(zhuǎn)換至UTM投影，WGS84-49N坐標(biāo)，沿研究區(qū)邊界進(jìn)行裁剪?；趖imesat3.0平臺(tái)，選用Savitsky-Golay方式對(duì)裁剪后的影像進(jìn)行濾波，從而剔除影像中因外強(qiáng)迫和儀器測(cè)量誤差等引起的噪聲和擾動(dòng)，避免此類(lèi)虛假信息造成誤判突變點(diǎn)[16]。然后，通過(guò)最大值合成將其轉(zhuǎn)換為16 d最大值合成數(shù)據(jù)，并將DN值轉(zhuǎn)換為地表溫度[20]。

為進(jìn)一步分析地表溫度突變與土地覆被變化的定量關(guān)系，對(duì)突變顯著的2002-06-08的Landsat 7 ETM+和2003-08-17 的Landsat 5 TM 共2景影像利用隨機(jī)森林分類(lèi)算法進(jìn)行土地利用分類(lèi)，結(jié)合土地利用分類(lèi)體系以及研究區(qū)特點(diǎn)，將研究區(qū)土地利用類(lèi)型分為耕地、林地、草地、建設(shè)用地、沙地及水域6類(lèi)。分類(lèi)影像下載于美國(guó)航空航天局(https://www.nasa.gov)，該影像為L(zhǎng)1T級(jí)數(shù)據(jù)，空間分辨率為30 m，影像質(zhì)量良好。L1T級(jí)數(shù)據(jù)已經(jīng)過(guò)系統(tǒng)幾何校正，達(dá)亞像元精度，對(duì)其進(jìn)行輻射校正、大氣校正、裁剪等預(yù)處理工作[21]。兩景影像分類(lèi)總體精度均達(dá)到90%以上，Kappa系數(shù)達(dá)到0.86以上。

1.3 研究方法

1.3.1 趨勢(shì)分析

為分析神東礦區(qū)地表溫度時(shí)間序列長(zhǎng)期變化趨勢(shì)，運(yùn)用最小二乘法逐像元擬合地表溫度年際變化趨勢(shì)，擬合斜率即為地表溫度隨時(shí)間變化的線性回歸系數(shù)，公式為

(1)

式中：Slope為擬合斜率；n為時(shí)間序列長(zhǎng)度；LSTi為第i年的LST值。

通過(guò)相關(guān)系數(shù)r判定變化趨勢(shì)的顯著程度，采用t檢驗(yàn)法對(duì)相關(guān)系數(shù)r的顯著性水平進(jìn)行檢驗(yàn)，公式為

(2)

(3)

式中，t的檢驗(yàn)自由度為n-2。

1.3.2 BFAST算法

BFAST算法最早由J.Verbesselt等[22]于2010年提出，該算法通過(guò)迭代時(shí)間序列的方法將原始時(shí)間序列分解為季節(jié)組分、趨勢(shì)組分和殘差組分，能夠有效檢測(cè)半干旱地區(qū)遙感影像時(shí)間序列中發(fā)生的季節(jié)與趨勢(shì)突變[23]。理論模型為

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

(4)

式中：Yt′為觀測(cè)數(shù)據(jù)；Tt′為趨勢(shì)組分；St′為季節(jié)組分；et′為殘差組分；t′為觀測(cè)時(shí)間。

Tt′=αi+βit′,

(5)

式中：i=1，…，m,m為突變點(diǎn)的總數(shù)；αi和βi分別為線性模型的截距和斜率。

(6)

式中：j為突變點(diǎn)所在位置，取值從1到q，q為突變點(diǎn)總數(shù)；k為調(diào)和項(xiàng)數(shù)目；αj,k為振幅；δj,k為相位；f為頻率(本研究采用的地表溫度遙感影像合成間隔為16 d，f為23)。

季節(jié)組分和趨勢(shì)組分中突變點(diǎn)的識(shí)別需要確定突變點(diǎn)的數(shù)量和突變點(diǎn)在時(shí)間序列中的位置。此處利用最小二乘移動(dòng)求和從季節(jié)組分和趨勢(shì)組分中檢驗(yàn)判斷是否存在突變點(diǎn)。季節(jié)組分和趨勢(shì)組分中突變點(diǎn)數(shù)量和位置的計(jì)算：首先分別剔除季節(jié)組分(Yt′-St′)和趨勢(shì)組分(Yt′-Tt′)，其次采用貝葉斯信息準(zhǔn)則確定突變點(diǎn)的最優(yōu)個(gè)數(shù)，最后通過(guò)最小二乘法估計(jì)突變點(diǎn)在時(shí)間序列中的位置。本文BFAST算法由R語(yǔ)言的bfast包實(shí)現(xiàn)(https://cran.r-project.org/web/packages/bfast/index.html)。

圖2 技術(shù)路線圖

2 研究結(jié)果

2.1 地表溫度長(zhǎng)期變化趨勢(shì)

2000—2018年，神東礦區(qū)98.63%的區(qū)域地表溫度呈下降趨勢(shì)，1.37%的區(qū)域地表溫度呈上升趨勢(shì)(表1)，但變化趨勢(shì)均未達(dá)到顯著性水平(p>0.05,-0.456

表1 地表溫度長(zhǎng)期變化趨勢(shì)表

2.2 地表溫度突變時(shí)空分布特征

本文重點(diǎn)分析地表溫度長(zhǎng)期變化趨勢(shì)中發(fā)生的突變，地表溫度趨勢(shì)組分突變的時(shí)空分布特征分為兩個(gè)方面：(1)地表溫度突變頻率的時(shí)空分布特征。統(tǒng)計(jì)不同年份不同區(qū)域的地表溫度突變分布。(2)地表溫度最大突變發(fā)生的時(shí)間和幅度。提取地表溫度發(fā)生一次或多次突變的區(qū)域中最大突變發(fā)生時(shí)間和幅度。

2.2.1 突變頻率的時(shí)空分布特征

2000—2018年，神東礦區(qū)地表溫度突變次數(shù)為0～4次，其中41.56%的區(qū)域地表溫度發(fā)生了至少1次突變，發(fā)生1～4次突變的區(qū)域分別占研究區(qū)面積的32.88%，7.86%，0.79%和0.03%(圖3(a))。井工采區(qū)內(nèi)地表溫度突變面積的比例遠(yuǎn)低于整體水平，僅19.02%的區(qū)域地表溫度發(fā)生突變，其中13.41%的區(qū)域地表溫度發(fā)生了1次突變，5.59%的區(qū)域地表溫度發(fā)生了2次突變，0.01%的區(qū)域內(nèi)地表溫度發(fā)生了3次突變。露天采區(qū)內(nèi)地表溫度突變面積的比例高于井工采區(qū)，露天采區(qū)內(nèi)34.66%的區(qū)域地表溫度發(fā)生突變，其中地表溫度發(fā)生了1次突變的區(qū)域面積占比達(dá)到32.16%，發(fā)生2次突變的區(qū)域面積占比為2.47%，發(fā)生了3次突變的區(qū)域面積占比為0.03%。

神東礦區(qū)地表溫度突變時(shí)間主要集中在2001—2005年和2014年，期間地表溫度突變總面積占神東礦區(qū)面積的44.81%，其中2003年地表溫度突變面積占比最高，高達(dá)20.07%(圖3(b))。2006—2018年，地表溫度發(fā)生突變面積大幅減少，呈現(xiàn)穩(wěn)定變化態(tài)勢(shì)。井工采區(qū)地表溫度在2002年發(fā)生突變面積最多，占井工采區(qū)面積的7.90%；露天采區(qū)地表溫度2014年發(fā)生突變面積最多，占露天采區(qū)面積的16.25%。

圖3 神東礦區(qū)地表溫度突變時(shí)空分布圖(2000—2018年)

2.2.2 最大突變發(fā)生時(shí)間和突變幅度

從地表溫度發(fā)生最大突變時(shí)間看，神東礦區(qū)2000，2006，2015，2016，2017和2018年共6年地表溫度未出現(xiàn)最大突變，最大突變出現(xiàn)的年份主要為2003和2014年，占發(fā)生突變區(qū)域面積的比例分別達(dá)到46.53%和13.38%(圖4(a))。井工采區(qū)內(nèi)地表溫度發(fā)生最大突變最多的年份為2002年，面積達(dá)到了65 km2。露天采區(qū)內(nèi)地表溫度發(fā)生最大突變最多的年份為2014年，面積達(dá)到了16 km2。

從地表溫度發(fā)生最大突變幅度看，神東礦區(qū)地表溫度發(fā)生最大突變幅度為-6.82～7.16 ℃，其中20.43%區(qū)域地表溫度負(fù)向突變，79.57%區(qū)域地表溫度正向突變(圖4(b))。井工采區(qū)地表溫度發(fā)生最大突變幅度為-6.80～6.48 ℃，其中30.22%區(qū)域地表溫度正向突變，69.78%區(qū)域地表溫度負(fù)向突變。露天采區(qū)地表溫度最大突變幅度為-6.55～5.10 ℃，其中80%區(qū)域地表溫度正向突變，20%區(qū)域地表溫度負(fù)向突變。

圖4 神東礦區(qū)地表溫度最大突變發(fā)生時(shí)間和幅度分布圖(2000—2018年)

2.3 地表溫度突變與土地覆被變化間關(guān)系

本研究選取2003年地表溫度發(fā)生最大突變的區(qū)域分析地表溫度突變與土地覆被變化間關(guān)系?；贓NVI變化檢測(cè)方法獲取2003年地表溫度發(fā)生最大突變區(qū)域土地覆被變化信息,通過(guò)空間疊加分析，統(tǒng)計(jì)得出土地覆被變化與地表溫度突變間的關(guān)系。從地表溫度突變幅度看，地表溫度突變均表現(xiàn)為正向突變，變化幅度在1.33～6.93 ℃間(圖5(a))。從土地覆被變化看，建設(shè)用地、沙地以及林地面積均有所增加，分別增加了15，5，3 km2；草地、耕地和水域面積均減少，分別耕地轉(zhuǎn)為沙地，地表溫度突變幅度為1.94～2.49 ℃；草地轉(zhuǎn)為沙地，地表溫度突變幅度為1.92～2.42 ℃。水域轉(zhuǎn)為草地，地表溫度突變幅度為2.21 ℃；林地轉(zhuǎn)為草地，地表溫度突變幅度為1.92～2.62 ℃；耕地轉(zhuǎn)為草地，地表溫度突變幅度為1.74～2.76 ℃。草地轉(zhuǎn)為林地，地表溫度突變幅度為1.84～2.82 ℃。林地轉(zhuǎn)為耕地，地表溫度突變幅度為2.32～2.61 ℃；草地轉(zhuǎn)為耕地，地表溫度突變幅度為1.33～2.79 ℃。草地轉(zhuǎn)為建設(shè)用地，地表溫度突變幅度最大，達(dá)1.89～6.93 ℃。

圖5 2003年地表溫度突變與土地覆被變化間的關(guān)系

3 討 論

植被在調(diào)節(jié)地表溫度方面扮演著重要角色，如植物葉片的蒸騰作用，持續(xù)吸收熱輻射，能夠有效調(diào)節(jié)濕度，降低溫度。[25]李宏宇等[26]研究表明，在中緯度地區(qū)，地表溫度隨植被覆蓋度的增加而降低。近年來(lái)，神東礦區(qū)的植被覆蓋狀況得到了較大改善，礦區(qū)整體植被覆蓋度由礦區(qū)開(kāi)采初期的3%～11%提高到59%以上[27]，其中馬家塔煤礦的植被覆蓋率達(dá)到80%[28]。整體植被覆蓋度的提高對(duì)地表產(chǎn)生的降溫效果進(jìn)一步證明了19年來(lái)神東礦區(qū)地表溫度長(zhǎng)期呈下降趨勢(shì)的合理性。

此外，就不同開(kāi)采方式而言，露天開(kāi)采比井工開(kāi)采更易發(fā)生地表溫度突變。可歸結(jié)為兩種開(kāi)采方式對(duì)地表擾動(dòng)程度的不同。神東礦區(qū)井工開(kāi)采除變形較大的區(qū)域外，采區(qū)與非采區(qū)并無(wú)明顯差異[29,30]。比如，沉陷區(qū)域的中性區(qū)隨著工作面推過(guò)，地表會(huì)逐漸穩(wěn)定，開(kāi)采過(guò)程中形成的臨時(shí)性地裂縫能夠自動(dòng)愈合。[31]相對(duì)于井工開(kāi)采，露天開(kāi)采對(duì)地表的破壞更為明顯和直接，露天開(kāi)采過(guò)程如巖土剝離、采掘場(chǎng)開(kāi)挖、排土場(chǎng)堆積等均對(duì)地表形態(tài)造成較大破壞，地表受到大幅度擾動(dòng)后，露天礦區(qū)產(chǎn)生明顯的增溫效應(yīng)[32]，本研究中露天采區(qū)地表溫度發(fā)生最大突變的區(qū)域80%為正向突變。

地表溫度突變方向取決于土地覆被變化類(lèi)型和規(guī)模，梁保平等[33]指出建設(shè)用地和沙地面積增加會(huì)導(dǎo)致地表溫度升高，原因?yàn)樯车睾徒ㄔO(shè)用地的地表潛熱容量小，接收太陽(yáng)輻射之后，地表溫度迅速升高。草地、耕地面積減少同樣能引起地表溫度的升高，原因?yàn)橹脖荒軌蛭仗?yáng)輻射并進(jìn)行蒸騰作用，可減少到達(dá)地表的熱量。水域比熱容大，其面積減少會(huì)直接造成地表溫度升高。林地面積增加雖然能夠有效降低地表溫度，但林地面積增加很小時(shí)，冷卻效應(yīng)并不明顯。[34]因此，2003年地表溫度發(fā)生最大突變的區(qū)域受建設(shè)用地、沙地面積增加，耕地、草地和水域面積減少影響較大，而受林地面積增加影響較小，地表溫度均表現(xiàn)為正向突變。同樣地，地表溫度突變幅度同樣取決于土地覆被變化，不同土地覆被的地表溫度存在較大差異[35]，因而不同土地覆被變化引起地表溫度突變的幅度也不盡相同。李潤(rùn)林等[36]研究表明在甘肅省張掖市，耕地、草地、林地、水域4種土地覆被與沙地的地表溫度差值為1.20～4.90 ℃。本研究中耕地、草地轉(zhuǎn)為沙地地表溫度突變幅度均為1.92～2.49 ℃。潘竟虎等[37]研究表明，甘肅省張掖市水域、耕地和建設(shè)用地3種土地覆被與草地的地表溫度差值為1.66～7.42 ℃，本研究中水域、林地、耕地轉(zhuǎn)為草地時(shí)地表溫度突變幅度為1.74～2.76 ℃。本文研究區(qū)與張掖市地理位置、氣候條件等具有可比性，因此認(rèn)為本文研究結(jié)果合理。

4 結(jié) 論

(1)2000—2018年，在我國(guó)西北地區(qū)向暖濕化轉(zhuǎn)型的背景下，神東礦區(qū)整體地表溫度呈下降趨勢(shì)，可歸因于礦區(qū)整體植被覆蓋的大幅改善，表明生態(tài)修復(fù)對(duì)區(qū)域氣候改善具有重要意義。

(2)露天采區(qū)比井工采區(qū)更容易發(fā)生地表溫度突變，露天采區(qū)內(nèi)地表溫度突變面積占比為34.66%，而井工采區(qū)內(nèi)僅19.02%的區(qū)域地表溫度發(fā)生突變，原因?yàn)槁短扉_(kāi)采對(duì)地表形態(tài)擾動(dòng)程度遠(yuǎn)高于井工開(kāi)采。露天采區(qū)地表溫度發(fā)生最大突變的區(qū)域中80%為溫度正向突變，表明露天開(kāi)采具有明顯的增溫效應(yīng)。

(3)土地覆被變化引起地表溫度發(fā)生突變，但地表溫度突變方向和幅度取決于土地覆被變化類(lèi)型、規(guī)模等因素。建設(shè)用地、沙地面積增加，耕地、草地、水域面積減少可導(dǎo)致地表溫度發(fā)生正向突變。在土地覆被轉(zhuǎn)變類(lèi)型中，草地轉(zhuǎn)為建設(shè)用地，地表溫度突變幅度最大，為1.89～6.93 ℃。本研究采用中分辨率地表溫度數(shù)據(jù)產(chǎn)品，下一步可選用分辨率更高的遙感影像，從而對(duì)地表溫度突變進(jìn)行更精細(xì)的檢測(cè)。此外，可將地表溫度突變檢測(cè)與土地覆被連續(xù)變化監(jiān)測(cè)相結(jié)合，以研究地表溫度突變與土地覆被變化驅(qū)動(dòng)間的長(zhǎng)期變化規(guī)律。