礦業(yè)密集區(qū)地表熱環(huán)境時(shí)空演變特征及驅(qū)動(dòng)力分析

袁運(yùn)許，朱小龍，侯春華，劉碩，宋建偉，李小光

(1.河北省地質(zhì)環(huán)境監(jiān)測(cè)院，石家莊 050000；2.河北省地質(zhì)資源環(huán)境監(jiān)測(cè)與保護(hù)重點(diǎn)實(shí)驗(yàn)室，石家莊 050000；3.唐山師范學(xué)院計(jì)算機(jī)科學(xué)系，唐山 063000；4.唐山市資源與環(huán)境遙感重點(diǎn)實(shí)驗(yàn)室，唐山 063210；5.華北理工大學(xué)礦業(yè)工程學(xué)院，唐山 063210；6.河北省礦業(yè)開發(fā)與安全技術(shù)重點(diǎn)實(shí)驗(yàn)室，唐山 063210；7.河北省礦區(qū)生態(tài)修復(fù)產(chǎn)業(yè)技術(shù)研究院，唐山 063210)

中國(guó)是全球主要的礦業(yè)大國(guó)，礦產(chǎn)資源的開發(fā)在促進(jìn)地方文化和經(jīng)濟(jì)等快速發(fā)展的同時(shí)，也帶來一些生態(tài)環(huán)境問題。由于礦產(chǎn)資源的不當(dāng)開挖，導(dǎo)致原有地表遭到破壞，改變了城市地表下墊面的熱特性，進(jìn)而使得植被蒸散能力下降，不透水材料對(duì)太陽(yáng)熱輻射的吸收大大增加[1-2]。地表溫度是區(qū)域與全球尺度地表物理過程分析與模擬的重要物理參量[3]。當(dāng)前主要是在全球城市區(qū)域尺度范圍內(nèi)，利用遙感技術(shù)反演地表溫度，以探究城市尺度的熱島效應(yīng)為主要內(nèi)容[4-6]，較少對(duì)鎮(zhèn)域尺度礦業(yè)開采活動(dòng)強(qiáng)烈的礦業(yè)密集區(qū)，地表熱環(huán)境時(shí)空演變特征及驅(qū)動(dòng)力進(jìn)行詳細(xì)分析。米曉楠等[7]利用多源衛(wèi)星資料計(jì)算太原市2003—2019年城市熱島強(qiáng)度和城市熱島比例指數(shù)，分析太原市城市熱島效應(yīng)時(shí)空變化特征。王美雅等[8]基于多時(shí)相Landsat遙感影像，對(duì)中外6個(gè)典型超大城市(北京、上海、廣州、倫敦、紐約和東京)進(jìn)行城市熱環(huán)境變化對(duì)比及成因分析。黃曦濤等[9]利用多期Landsat遙感影像分析了西安市熱島效應(yīng)的分布與變化規(guī)律。

鄉(xiāng)鎮(zhèn)是中國(guó)基層鄉(xiāng)村發(fā)展管理和行政區(qū)劃的最基本單元[10]。對(duì)于鎮(zhèn)域尺度礦業(yè)密集區(qū)，無論是礦山企業(yè)數(shù)量，還是礦山企業(yè)所占用土地?cái)?shù)量來說都是相對(duì)集中的區(qū)域[11-12]。礦區(qū)及周邊區(qū)域的地表熱環(huán)境，會(huì)因?yàn)榈V業(yè)開采導(dǎo)致的原有地表覆被類型改變而發(fā)生變化[13-14]。但大多基于礦區(qū)尺度進(jìn)行研究。畢銀麗等[15]基于Landsat TM/OLI與 MODIS 衛(wèi)星遙感數(shù)據(jù)，分析了準(zhǔn)格爾礦區(qū)1988—2020年的植被覆蓋與地表溫度變化情況，結(jié)果表明：33a來準(zhǔn)格爾礦區(qū)外圍未擾動(dòng)區(qū)域植被覆蓋度呈整體上升趨勢(shì)，礦區(qū)地表溫度呈下降趨勢(shì)。韓凱旭[16]采用1990—2019年Landsat-8 OLI/TIRS與Landsat-5 TM系列遙感影像結(jié)合MODIS-MOD11 L2數(shù)據(jù)，結(jié)合土地覆被時(shí)空間格局變化規(guī)律和地表溫度空間分布，探究加拿大Alberta省油砂礦礦區(qū)油砂礦開采對(duì)地表溫度變化的影響作用。

開展鎮(zhèn)域尺度礦業(yè)密集區(qū)地表熱環(huán)境時(shí)空演變特征及驅(qū)動(dòng)力分析，是分析礦業(yè)開采擾動(dòng)規(guī)律和制定區(qū)域環(huán)境治理政策的關(guān)鍵。遙感技術(shù)在宏觀上開展礦區(qū)生態(tài)環(huán)境監(jiān)測(cè)具有不可比擬的優(yōu)勢(shì)[16]，為此，基于Landsat遙感影像數(shù)據(jù)的熱紅外波段，利用輻射傳輸方程法(radiative transfer equation，RTE)反演地表溫度，并進(jìn)行定性和定量精度驗(yàn)證?；贚andsat遙感影像數(shù)據(jù)的可見光波段，提取礦業(yè)密集區(qū)地表下墊面覆被類型及礦業(yè)用地時(shí)空動(dòng)態(tài)變化信息。在此基礎(chǔ)上，從地表下墊面覆被類型角度，借助“黃藍(lán)—差值影像法”、疊加分析和剖面分析等方法，對(duì)2000—2018年礦業(yè)密集區(qū)地表熱環(huán)境時(shí)空演變特征及驅(qū)動(dòng)力進(jìn)行分析。以期深化礦業(yè)密集區(qū)地表熱環(huán)境時(shí)空演變特征規(guī)律的認(rèn)識(shí)，為制定緩解高溫聚集效應(yīng)策略提供定量參考依據(jù)。

1 研究方法

1.1 研究區(qū)概況及數(shù)據(jù)源

1.1.1 研究區(qū)概況

馬蘭莊鎮(zhèn)位于河北省遷安市西北部，距市區(qū) 15 km，地處京、津、唐、秦腹地。東臨灤河，與遷安鎮(zhèn)、閻家店鄉(xiāng)、大崔莊鎮(zhèn)三個(gè)鎮(zhèn)鄉(xiāng)僅一河之隔，南連蔡園鎮(zhèn)，西、北與遷西縣接壤(圖1)。西部群山環(huán)抱，三面環(huán)山、一面鄰水，為半山區(qū)燕山隆起帶余脈南麓，總面積49.16 km2，耕地面積530.83 hm2[12]。

馬蘭莊鎮(zhèn)鎮(zhèn)內(nèi)礦產(chǎn)資源豐富，采礦場(chǎng)分布密集，是鐵礦石資源豐富的鋼鐵重鎮(zhèn)。全鎮(zhèn)現(xiàn)有礦山企業(yè)超過50家。2001年遷安市境內(nèi)19個(gè)鄉(xiāng)鎮(zhèn)，礦業(yè)用地密集程度的聚類分析統(tǒng)計(jì)結(jié)果表明，馬蘭莊鎮(zhèn)境內(nèi)的礦業(yè)用地占比為34.80%，聚類分析結(jié)果表明該鄉(xiāng)鎮(zhèn)礦業(yè)開發(fā)密集程度為1級(jí)，選取馬蘭莊鎮(zhèn)為重點(diǎn)分析的鎮(zhèn)域尺度礦業(yè)密集區(qū)[17]。馬蘭莊鎮(zhèn)境內(nèi)的礦業(yè)大多屬于露天開采方式，開采擾動(dòng)的具體表現(xiàn)形式是破壞地表原有植被，導(dǎo)致大量巖石和土壤出露。

1.1.2 數(shù)據(jù)源

(1)遙感影像數(shù)據(jù)。選取與研究區(qū)礦業(yè)開發(fā)歷史變化緊密相關(guān)的5個(gè)典型時(shí)期的影像(表1)。分別代表礦業(yè)開發(fā)不景氣時(shí)期、礦業(yè)開始大規(guī)模開采初期、礦業(yè)破壞最嚴(yán)重時(shí)期同時(shí)也是礦區(qū)生態(tài)修復(fù)起始時(shí)期、礦區(qū)生態(tài)修復(fù)效果初顯時(shí)期、礦區(qū)生態(tài)修復(fù)效果較顯著時(shí)期。

表1 Landsat遙感影像數(shù)據(jù)

由于夏季(6—8月)太陽(yáng)高度角大，太陽(yáng)輻射較強(qiáng)，地表溫度較高。并且不同下墊面覆被類型的比熱不同，導(dǎo)致地表熱環(huán)境異質(zhì)性特征尤為明顯，進(jìn)而能夠較為準(zhǔn)確地反映礦業(yè)密集區(qū)地表下墊面覆被和地表溫度的差異，因此選用數(shù)據(jù)時(shí)優(yōu)先考慮夏季數(shù)據(jù)。但因?yàn)檫w安市自7月初氣溫持續(xù)上升，7月初―9月初處于高溫期，因此選用4期8月初―9月初數(shù)據(jù)，由于2013年7月初—9月初數(shù)據(jù)云量較多，影響地表溫度反演效果，故2013年選取9月底影像作為數(shù)據(jù)源。

Landsat 8攜帶有OLI陸地成像儀和TIRS熱紅外傳感器，其中OLI陸地成像儀包括9個(gè)波段，TIRS包括2個(gè)單獨(dú)的熱紅外波段。

(2)輔助數(shù)據(jù)。谷歌地球高分辨率影像，用于對(duì)地表覆被類型進(jìn)行目視解譯的精度驗(yàn)證。

1.2 地表溫度反演及精度驗(yàn)證

1.2.1 基于輻射傳輸方程法反演地表溫度

中外專家學(xué)者根據(jù)不同衛(wèi)星傳感器熱紅外波段數(shù)據(jù)的差異，提出了不同的溫度反演算法。其中基于Landsat數(shù)據(jù)的熱紅外波段進(jìn)行地表溫度的反演算法主要有輻射傳輸方程法、單窗算法(mono-window algorithm，MWA)和單通道算法(single-channel method, SC)。其中輻射傳輸方程法計(jì)算過程較為復(fù)雜，但反演精度較為準(zhǔn)確，精度達(dá) 0.6 ℃。單窗算法計(jì)算過程比較簡(jiǎn)單，但是較難獲取大氣平均溫度和大氣透射率等參數(shù)[18]。單通道算法操作過程簡(jiǎn)單，但受大氣水汽含量影響較大，由于遷安市夏季大氣水汽含量較大，因此不適合用該模型來反演地表溫度。

結(jié)合以上3種算法的特點(diǎn)，選取輻射傳輸方程法(又稱為大氣校正法)反演地表溫度，詳細(xì)步驟參見文獻(xiàn)[19]。

1.2.2 精度驗(yàn)證

利用氣象站實(shí)測(cè)地溫進(jìn)行精度驗(yàn)證。由于馬蘭莊鎮(zhèn)沒有氣象站，因此利用遷安市氣象站數(shù)據(jù)進(jìn)行驗(yàn)證。由于遷安市氣象站缺少2000年、2003年和2008年實(shí)測(cè)地溫?cái)?shù)據(jù)，因此只對(duì)2013年和2018年LST反演結(jié)果進(jìn)行精度驗(yàn)證。

2013年9月26日遷安市境內(nèi)54439號(hào)氣象站，地溫實(shí)測(cè)點(diǎn)的經(jīng)緯度分別為118°43′12″和 40°1′12″。2018年9月8日遷安市境內(nèi)54439號(hào)氣象站，地溫實(shí)測(cè)點(diǎn)的經(jīng)緯度分別為118°42′36″和 40°1′12″。兩個(gè)年份的實(shí)測(cè)點(diǎn)均位于Landsat影像的同一像元內(nèi)。

經(jīng)查驗(yàn)該氣象站在2013年9月26日10：00點(diǎn)，0 cm和5 cm實(shí)測(cè)地溫分別為26 ℃和19.2 ℃。11：00，0 cm和5 cm實(shí)測(cè)地溫分別為29.2 ℃和22.3 ℃。2013年9月26日Landsat衛(wèi)星在遷安市過境時(shí)間為10：49。對(duì)照谷歌高分影像發(fā)現(xiàn)該時(shí)刻氣象站地溫實(shí)測(cè)點(diǎn)位于稀疏草地像元內(nèi)。因?yàn)榈乇頊囟韧ǔ１欢x為地球表面的溫度，在植被覆蓋度較低的區(qū)域，通常利用植被冠層以及土壤的平均溫度表征；在植被覆蓋度較高的區(qū)域，通常利用植被冠層的溫度表征；在裸地區(qū)域，通常利用土壤的表面溫度表征[20]。因此用0 cm(裸地地表)和5 cm(草地冠層)實(shí)測(cè)地溫均值來表示稀疏草地內(nèi)的地溫值。

由于實(shí)測(cè)地溫為地表冠層的整點(diǎn)溫度，因此首先利用線性插值法[21]，內(nèi)插得到與衛(wèi)星過境時(shí)刻相對(duì)應(yīng)的0 cm(裸地地表)和5 cm(草地冠層)的實(shí)測(cè)地溫分別為28.62 ℃和21.74 ℃，計(jì)算二者均值得到稀疏草地實(shí)測(cè)地溫值為25.18 ℃。經(jīng)查驗(yàn)Landsat衛(wèi)星過境時(shí)刻與該氣象站實(shí)測(cè)地溫同經(jīng)緯度的像元內(nèi)地溫反演值為25.09 ℃。計(jì)算得出2013年9月26日相同時(shí)刻相同經(jīng)緯度，地溫實(shí)測(cè)值和反演值之間誤差為0.09 ℃。

經(jīng)查驗(yàn)該氣象站在2018年9月8日10：00，0 cm 和5 cm實(shí)測(cè)地溫分別為23.3 ℃和26.5 ℃。11：00，0 cm和5 cm實(shí)測(cè)地溫分別為24.2 ℃和29.5 ℃。2018年9月8日Landsat衛(wèi)星在遷安市過境時(shí)間為10：46。對(duì)照谷歌高分影像發(fā)現(xiàn)該時(shí)刻氣象站地溫實(shí)測(cè)點(diǎn)位于濃密草地像元內(nèi)，因此用5 cm(草地冠層)實(shí)測(cè)地溫值來表示濃密草地內(nèi)的地溫值。

由于實(shí)測(cè)地溫為地表冠層的整點(diǎn)溫度，因此利用線性插值法[16]，內(nèi)插得到與衛(wèi)星過境時(shí)刻相對(duì)應(yīng)的5 cm(草地冠層)的實(shí)測(cè)地溫，即得到濃密草地實(shí)測(cè)地溫值為28.81 ℃。經(jīng)查驗(yàn)Landsat衛(wèi)星過境時(shí)刻與該氣象站同經(jīng)緯度的像元內(nèi)地溫反演值為28.68 ℃。計(jì)算得出2018年9月8日相同時(shí)刻相同經(jīng)緯度，地溫實(shí)測(cè)值與反演值之間的誤差為0.13 ℃。

可知2013年和2018年地溫實(shí)測(cè)值與反演值之間誤差較小，反演結(jié)果精度較高。

1.3 礦業(yè)用地信息提取

1.3.1 影像分割和合并

基于ENVI5.3軟件平臺(tái)面向?qū)ο筇卣魈崛∧K feature extraction的基于規(guī)則的面向?qū)ο笥跋穹诸愃惴ǎ瑢?duì)馬蘭莊鎮(zhèn)境內(nèi)礦業(yè)用地進(jìn)行信息提取。

參考基于樣本的面向?qū)ο笥跋穹诸愃惴ㄖ蟹指詈秃喜?shí)驗(yàn)方法，在實(shí)驗(yàn)過程中通過預(yù)覽分割和合并效果確定最優(yōu)分割尺度(圖2)。當(dāng)分割尺度設(shè)置為10時(shí)，地物邊界清晰，但存在過分割現(xiàn)象，分割效果較破碎。當(dāng)分割尺度設(shè)置為30時(shí)，地物邊界混雜，存在欠分割現(xiàn)象，不同地物可能會(huì)合并到同一地物類型中。最終經(jīng)過多次實(shí)驗(yàn)，當(dāng)分割尺度為20，合并尺度為50時(shí)，各類地物邊緣特征能夠得到有效識(shí)別，因此作為最優(yōu)閾值參與后續(xù)分類實(shí)驗(yàn)。

圖2 2018年馬蘭莊鎮(zhèn)Landsat影像不同分割合并尺度下的效果

1.3.2 特征分析

以2018年影像為例，圖3展示了研究區(qū)不同地表下墊面覆被類型的平均光譜曲線。5類地物在海岸波段(第1波段)、藍(lán)波段(第2波段)、綠波段(第3波段)和紅波段(第4波段)的DN值相差不大，區(qū)分度不高。礦業(yè)用地僅在在近紅外波段(第5波段)與其他地物表現(xiàn)出了明顯差異，5類地物像元亮度值(digital number, DN)在該波段排序?yàn)榱值?耕地>居民地>礦業(yè)用地>水域。因此不能單獨(dú)利用光譜信息來提取礦業(yè)用地信息，需要結(jié)合紋理特征和空間特征對(duì)研究區(qū)礦業(yè)用地進(jìn)行信息提取。

圖3 2018年研究區(qū)地表下墊面覆被類型光譜曲線

1.3.3 精度驗(yàn)證

在谷歌地球高分辨率影像上，利用隨機(jī)選點(diǎn)的方式選擇礦業(yè)用地和其他地類各50個(gè)樣本點(diǎn)，對(duì)鄉(xiāng)鎮(zhèn)A境內(nèi)礦業(yè)用地信息提取結(jié)果進(jìn)行了精度驗(yàn)證。表2為礦業(yè)用地信息提取混淆矩陣，總體分類精度和Kappa系數(shù)分別為93.00%和0.86。礦業(yè)用地的制圖精度和用戶精度分別為92.00%和93.88%，說明礦業(yè)用地信息提取結(jié)果精度較高，可用于下一步分析。

表2 礦業(yè)用地信息提取混淆矩陣

2 結(jié)果與分析

2.1 礦業(yè)用地時(shí)空動(dòng)態(tài)演變分析

2.1.1 礦業(yè)用地空間變化分析

馬蘭莊鎮(zhèn)境內(nèi)礦業(yè)用地的演變過程以2008年為界，分為兩個(gè)階段。2000年以前礦業(yè)處于開發(fā)不景氣時(shí)期，境內(nèi)的礦區(qū)占地面積較少，主要集中于北部地區(qū)，東南部地區(qū)零星分布。2000年以后隨著礦業(yè)的大規(guī)模開采，采礦境界不斷擴(kuò)張，2003年在2000年基礎(chǔ)上沿原開采境界向外圍擴(kuò)展，2008年在2003年基礎(chǔ)上向外圍擴(kuò)張。2008年以后隨著礦區(qū)生態(tài)修復(fù)工作的逐步實(shí)施，2013年礦區(qū)占地在2008年基礎(chǔ)上逐步縮小，2018年在2013年基礎(chǔ)上逐步縮小，縮小的區(qū)域主要被植被覆蓋(圖4)。

圖4 5期影像馬蘭莊鎮(zhèn)境內(nèi)礦業(yè)用地分布

2.1.2 礦業(yè)用地時(shí)間演變統(tǒng)計(jì)

由表3可以看出，2000—2018年，礦區(qū)占地面積變化明顯，以2008年為界，分為兩個(gè)階段。第一個(gè)階段是2000—2008年，2000年以前礦業(yè)處于開發(fā)不景氣時(shí)期，境內(nèi)的礦區(qū)占地面積較少，對(duì)生態(tài)破壞不大，2000年礦業(yè)用地占比僅為25.78%。2000年以后隨著國(guó)內(nèi)市場(chǎng)對(duì)于礦產(chǎn)資源需求的加大，礦產(chǎn)資源開采規(guī)模急劇擴(kuò)張，大量植被和耕地遭到破壞，大量土地被壓占，生態(tài)環(huán)境破壞程度加劇，礦業(yè)用地占比逐年增加，2003年占比增加到31.91%。2008年是礦業(yè)破壞最嚴(yán)重時(shí)期，礦區(qū)面積占比高達(dá)38.78%。第二個(gè)階段是2008—2018年，由于礦業(yè)開采導(dǎo)致的生態(tài)環(huán)境問題日益突出，從2008 年開始，政府開始重視環(huán)境保護(hù)工作，礦區(qū)生態(tài)修復(fù)開始逐步實(shí)施，礦業(yè)用地占比逐年縮小，2013年礦區(qū)生態(tài)修復(fù)效果初顯，礦區(qū)面積占比下降到32.07%，礦區(qū)生態(tài)問題得到一定程度緩解。直至2018年礦區(qū)面積占比下降到27.94%，礦區(qū)生態(tài)修復(fù)效果初見成效。由此可知，2000—2003年和2003—2008年這兩個(gè)時(shí)期，是礦業(yè)用地占比的急劇上升時(shí)期，2008—2013年和2013—2018年這兩個(gè)時(shí)期，是礦業(yè)用地占比的緩慢下降時(shí)期。

表3 5期影像馬蘭莊鎮(zhèn)境內(nèi)礦業(yè)用地變化統(tǒng)計(jì)

為進(jìn)一步定量化和可視化的探究馬蘭莊鎮(zhèn)境內(nèi)，2000—2018年礦業(yè)用地時(shí)空動(dòng)態(tài)變化過程，利用遙感變化檢測(cè)技術(shù)黃藍(lán)-差值影像法，分別獲得了5個(gè)時(shí)段(2000—2003年、2003—2008年、2008—2013年、2013—2018年、2000-2018年)，礦業(yè)用地和非礦業(yè)用地的變化影像(圖5)，并統(tǒng)計(jì)了礦業(yè)用地和非礦業(yè)用地面積占比變動(dòng)情況。

黃色區(qū)域表示礦業(yè)用地增加區(qū)域(非礦業(yè)用地減少區(qū)域)；藍(lán)色區(qū)域表示礦業(yè)用地減少區(qū)域(非礦業(yè)用地增加區(qū)域)

2000—2008年，由于礦產(chǎn)資源的大規(guī)模開采，大量綠地被占用，致使礦業(yè)用地在2000—2003年和2003—2008年兩個(gè)時(shí)段內(nèi)大幅度增加。2003年礦業(yè)用地占比相較于2000年增加了6.71%，2008年礦業(yè)用地占比相較于2003年增加了6.87%。

2008—2018年，礦區(qū)生態(tài)環(huán)境質(zhì)量在深度治理中明顯改善，使得礦業(yè)用地在2008—2013年和2013—2018年這兩個(gè)時(shí)段內(nèi)大幅度減少。2013年礦業(yè)用地占比相較于2008年減少了8.39%，2018年礦業(yè)用地占比相較于2013年減少了4.13%。雖然2008年以后政府開始重視環(huán)境保護(hù)工作，但是由于生態(tài)修復(fù)效果顯現(xiàn)時(shí)效滯后性的特點(diǎn)，2018年礦業(yè)用地占比仍高于2000年2.16%?；凇包S藍(lán)-差值影像法”得到的鄉(xiāng)鎮(zhèn)A境內(nèi)5個(gè)時(shí)段的礦業(yè)用地變動(dòng)情況，與鄉(xiāng)鎮(zhèn)A境內(nèi)5個(gè)年份礦業(yè)用地的時(shí)空演變特征相吻合。

2.2 地表熱環(huán)境特征疊加分析

為了定量分析礦業(yè)開發(fā)對(duì)地表熱環(huán)境的擾動(dòng)規(guī)律，對(duì)歸一化地表溫度(normalized land surface temperature, NLST)進(jìn)行均值-標(biāo)準(zhǔn)差分級(jí)。根據(jù)NLST的取值范圍，依托ArcGIS10.2軟件，將NLST劃分為5個(gè)熱力等級(jí)(表4)，更加系統(tǒng)地反映研究區(qū)地表溫度等級(jí)分布狀況[22]。結(jié)果表明，5個(gè)年份境內(nèi)地表熱環(huán)境均呈現(xiàn)顯著的空間異質(zhì)性和規(guī)律性。林地和水域分別是低溫和次低溫區(qū)，耕地主要為中溫區(qū)，居民地是次高溫區(qū)，礦業(yè)用地為高溫區(qū)(圖6)。

圖6 5期影像鄉(xiāng)鎮(zhèn)A地表溫度等級(jí)分布

表4 NLST等級(jí)劃分標(biāo)準(zhǔn)

統(tǒng)計(jì)5期影像鄉(xiāng)鎮(zhèn)A境內(nèi)地表溫度等級(jí)面積占比(表5)，2000年境內(nèi)高溫區(qū)面積占比為25.11%，隨著礦業(yè)用地的擴(kuò)張，2003年和2008年高溫區(qū)面積占比分別達(dá)到31.05%和39.93%；隨著2008年之后礦區(qū)生態(tài)修復(fù)工作的大規(guī)模實(shí)施，礦業(yè)用地占比逐年縮小，2013年境內(nèi)高溫區(qū)面積占比下降到30.75%，隨著生態(tài)修復(fù)工作的推進(jìn)，直至2018年境內(nèi)高溫區(qū)面積占比下降到28.71%，礦區(qū)生態(tài)修復(fù)效果初見成效。由此可知，2000—2003年和2003—2008年境內(nèi)高溫區(qū)面積占比處于急劇增加時(shí)期，2008—2013年和2013—2018年，境內(nèi)高溫區(qū)面積占比處于緩慢下降時(shí)期。高溫區(qū)面積占比與境內(nèi)礦業(yè)用地面積占比的波動(dòng)情況基本一致，并且經(jīng)過疊加分析得知境內(nèi)高溫區(qū)范圍與礦業(yè)用地占地范圍基本重合。

表5 5期影像鄉(xiāng)鎮(zhèn)A境內(nèi)NLST各等級(jí)面積占比

2.3 地表熱環(huán)境特征剖面分析

在地表熱環(huán)境研究中，剖面分析方法被廣泛應(yīng)用于表征不同空間上的溫度差異[23]。貫穿礦業(yè)用地的溫度剖面線具有較大波動(dòng)，呈現(xiàn)出多個(gè)波峰和波谷特征，能較好地從空間上清晰地展現(xiàn)礦業(yè)開發(fā)密集區(qū)不同地表下墊面覆被類型的地表熱環(huán)境差異。

以2000年、2008年和2018年3期地表溫度數(shù)據(jù)為基礎(chǔ)，在ArcGIS10.2軟件的支持下，選取具有代表性的貫穿大石河礦區(qū)、馬蘭莊礦區(qū)和水廠礦區(qū)等大型礦區(qū)的4條剖面線A-A、B-B、C-C和D-D，對(duì)馬蘭莊鎮(zhèn)境內(nèi)不同地表下墊面覆被類型對(duì)地表熱環(huán)境的影響進(jìn)行剖面分析(圖7)。

圖7 鄉(xiāng)鎮(zhèn)A境內(nèi)4條剖面線位置

為了使不同年份的地表溫度數(shù)據(jù)能在統(tǒng)一的量綱下進(jìn)行對(duì)比，參考文獻(xiàn)[24]，首先對(duì)地表溫度進(jìn)行正規(guī)化處理NLST，其計(jì)算公式[25]為

NLST=(LST-LSTmin)/(LSTmax-LSTmin)

(1)

式(1)中：NLST為正規(guī)化后的地表溫度值；LST為原始地表溫度值；LSTmax為L(zhǎng)ST影像中的最大值；LSTmin為L(zhǎng)ST影像中的最小值。

2.3.1 剖面線A-A溫度與地表下墊面覆被類型關(guān)系

剖面線A-A長(zhǎng)度約為7 300 m，橫穿1個(gè)礦區(qū)采場(chǎng)、1個(gè)排土場(chǎng)和1尾礦庫(kù)(圖8)。

黃色溫度剖面線表示2000年地表溫度；橙色溫度剖面線表示2008年地表溫度；紅色溫度剖面線表示2018年地表溫度

整個(gè)剖面線范圍內(nèi)，2008年和2018年溫度剖面線總體變化趨勢(shì)基本一致。2000年地表溫度普遍低于2008和2018年，這是由于2000年之后礦業(yè)大規(guī)模開發(fā)占用了大量綠地，出現(xiàn)的裸露地表具有熱傳導(dǎo)率大，熱慣量和熱容量小的熱學(xué)特性，使得地表溫度在太陽(yáng)直接輻射作用下迅速升高。最為明顯的是2 000～3 000 m范圍，2000年該范圍內(nèi)覆蓋著大量地表原有植被，2000年之后隨著礦產(chǎn)資源開采工作的實(shí)施，原有地表覆被逐漸被剝離，大量植被遭到破壞，形成了礦區(qū)采場(chǎng)等裸露地表，因此2008年和2018年的溫度剖面線明顯高于2000年。4 900～5 800 m的礦區(qū)排土場(chǎng)范圍，由于2008年之后對(duì)該排土場(chǎng)進(jìn)行的生態(tài)修復(fù)工作成效顯著，大量綠色植被覆蓋其上，使得到達(dá)地表的一部分太陽(yáng)輻射被阻擋和吸收，蒸騰量增加，與裸地相比地溫上升較慢，因此2018年的溫度剖面線明顯低于2000年和2008年。

整個(gè)剖面線內(nèi)的低溫區(qū)(波谷)位于綠地和水域。高溫區(qū)(波峰)主要分布于1 400～3 000 m的礦區(qū)采場(chǎng)裸巖，4 800～5 800 m未修復(fù)排土場(chǎng)礦渣和5 900～6 800 m的尾礦庫(kù)尾砂。整個(gè)剖面線范圍內(nèi)的高溫和低溫覆被類型交替出現(xiàn)，導(dǎo)致溫度上升和下降幅度較大，溫度剖面線“波峰”和“波谷”轉(zhuǎn)化較快，振幅較大。

2.3.2 剖面線B-B溫度與地表下墊面覆被類型關(guān)系

剖面線B-B長(zhǎng)度約為9 700 m，橫穿2個(gè)礦區(qū)采場(chǎng)、1個(gè)排土場(chǎng)和1個(gè)尾礦庫(kù)(圖9)。

黃色溫度剖面線表示2000年地表溫度；橙色溫度剖面線表示2008年地表溫度；紅色溫度剖面線表示2018年地表溫度

可以看出，300～1 400 m范圍，雖然3個(gè)年份在該區(qū)域內(nèi)都是礦區(qū)采場(chǎng)裸巖，但是由于2008年和2018年礦區(qū)采場(chǎng)裸巖內(nèi)存在部分陰影，陰影區(qū)域由于太陽(yáng)短波輻射被遮蔽，導(dǎo)致地溫升溫較慢，因此2008年和2018年的溫度剖面線普遍低于2000年。3 000～3 900 m 和5 400～6 200 m范圍，2000年的溫度剖面線普遍低于2008年和2018年，這是由于2000年的大量綠地區(qū)域，被轉(zhuǎn)變?yōu)?008年和2018年的礦業(yè)用地和居民地等裸露地表和不透水面，裸地和不透水面的熱傳導(dǎo)系數(shù)大、熱慣量和熱容量小，吸收太陽(yáng)輻射后升溫較快。3 900～5 000 m范圍，2000年溫度剖面線普遍高于2008年和2018年，這是因?yàn)?008年和2018年對(duì)該區(qū)域的尾礦庫(kù)進(jìn)行了生態(tài)修復(fù)工作，植物的光合作用和蒸騰作用會(huì)消耗較多能量，使得地表溫度升溫較慢。從 7 000 m 開始，2008年和2018年的溫度剖面線與2000年呈現(xiàn)相反的變化趨勢(shì)，主要原因是2000年的綠地和水域，轉(zhuǎn)化為了2008年和2018年的居民地，植被和水體的減少使得蒸發(fā)蒸騰作用大幅度降低，導(dǎo)致地溫急劇上升。

整個(gè)剖面線內(nèi)的低溫區(qū)(波谷)位于綠地和水域。高溫區(qū)(波峰)分布于300～1 400 m的礦區(qū)采場(chǎng)裸巖，1 500～3 000 m的排土場(chǎng)礦渣，3 900～4 900 m 的尾礦庫(kù)尾砂，以及5 400～6 200 m的礦區(qū)采場(chǎng)裸巖，由于這4個(gè)區(qū)域內(nèi)的高溫和低溫覆被類型交替出現(xiàn)，因此溫度剖面線“波峰”和“波谷”交替突變，振幅較大。

2.3.3 剖面線C-C溫度與地表下墊面覆被類型關(guān)系

剖面線C-C長(zhǎng)度約為7 100 m，橫穿1個(gè)尾礦庫(kù)和2個(gè)礦區(qū)采場(chǎng)(圖10)。

黃色溫度剖面線表示2000年地表溫度；橙色溫度剖面線表示2008年地表溫度；紅色溫度剖面線表示2018年地表溫度

可以看出，4 500 m以內(nèi)的區(qū)域，2000年、2008年和2018年的溫度剖面線出現(xiàn)交替變化，主要原因是2000年大部分原有綠地轉(zhuǎn)化為了2008年和2018年堆積尾砂的尾礦庫(kù)。植被具有遮蔽太陽(yáng)輻射和蒸騰吸熱作用，而尾礦庫(kù)內(nèi)尾砂較之植被能夠吸收和發(fā)射更多的長(zhǎng)波輻射，對(duì)地-氣間的顯熱交換更強(qiáng)烈，導(dǎo)致地表升溫較快。4 500～7 100 m范圍，2000年的溫度剖面線普遍高于2008年和2018年，這是由于2008年和2018年溫度剖面線貫穿礦區(qū)采場(chǎng)陰影區(qū)域，陰影區(qū)域由于太陽(yáng)短波輻射被遮蔽，因而地表溫度普遍較低。

整個(gè)剖面線范圍內(nèi)的低溫區(qū)(波谷)主要位于綠地和水域。高溫區(qū)(波峰)分布于500～2 600 m的未修復(fù)尾礦庫(kù)，2 700～4 100 m以及4 300～7 000 m的礦區(qū)采場(chǎng)裸巖。

2.3.4 剖面線D-D溫度與地表下墊面覆被類型關(guān)系

剖面線D-D長(zhǎng)度約為6 400 m，橫穿2個(gè)排土場(chǎng)和1個(gè)礦區(qū)采場(chǎng)(圖11)。

黃色溫度剖面線表示2000年地表溫度；橙色溫度剖面線表示2008年地表溫度；紅色溫度剖面線表示2018年地表溫度

可以看出，0～1 250 m范圍，2008年和2018年的溫度剖面線變化趨勢(shì)較一致，2000年溫度剖面線明顯低于2008年和2018年，這是由于2000年地表原有植被較多，而2008年和2018年大量綠地占用使得裸露地表出現(xiàn)，地表能夠吸收更多的太陽(yáng)輻射，并且裸地極為干燥，熱容量小，導(dǎo)致地溫急劇升高。1 500～2 500 m范圍，2000年和2008年溫度剖面線普遍高于2018年，主要原因是2008年之后對(duì)該區(qū)域的排土場(chǎng)進(jìn)行了生態(tài)修復(fù)，栽植了大量綠色植被，當(dāng)植被覆蓋密度較大時(shí)，通過蒸發(fā)蒸騰作用帶走較多熱量，因而2018年地表溫度明顯低于2000年和2008年。2 750～3 250 m范圍，2000年溫度剖面線明顯高于2008年和2018年，這是由于2008年和2018年該排土場(chǎng)內(nèi)水域面積增大，水體在各地表下墊面覆被類型中具有最大的熱容量和反射能力以及最強(qiáng)的蒸散作用，溫度上升緩慢，因此地表溫度最低。3 250～6 300 m范圍，3個(gè)年份地表溫度呈現(xiàn)交替相反的波段趨勢(shì)，這是由于該區(qū)域內(nèi)排土場(chǎng)和采場(chǎng)境內(nèi)的水域、綠地、陰影和裸巖、礦渣等地物交替出現(xiàn)導(dǎo)致的。

3 結(jié)論

在探究鎮(zhèn)域尺度礦業(yè)開發(fā)密集區(qū)，午間地表熱環(huán)境特征及驅(qū)動(dòng)力過程中，基于輻射傳輸方程法(RTE)反演地表溫度，并進(jìn)行定量精度驗(yàn)證?；凇包S藍(lán)差值影像法”、疊加分析和剖面分析等方法，揭示礦業(yè)開采過程對(duì)地表熱環(huán)境動(dòng)態(tài)變化的驅(qū)動(dòng)力，彌補(bǔ)了礦業(yè)開采導(dǎo)致礦業(yè)密集區(qū)地表熱環(huán)境特征及驅(qū)動(dòng)力分析的不足。得出以下主要結(jié)論。

(1)利用5期Landsat遙感影像熱紅外波段，基于RTE反演的地表溫度結(jié)果，與氣象站實(shí)測(cè)地溫?cái)?shù)據(jù)之間誤差約為0.1 ℃，表明利用Landsat影像熱紅外數(shù)據(jù)基于RTE，在無地面同步資料情況下進(jìn)行地表溫度反演精度較高，能夠滿足研究區(qū)地表熱環(huán)境異質(zhì)性研究的需要。

(2)3期影像各條溫度剖面線的波谷(低溫區(qū))均位于綠地和水域，波峰(高溫區(qū))均位于礦區(qū)采場(chǎng)裸巖、未進(jìn)行生態(tài)修復(fù)的排土場(chǎng)礦渣和尾礦庫(kù)內(nèi)尾砂。主要是由于礦業(yè)用地的裸地占比較大，裸地由于沒有植被覆蓋，因此很難通過植被的蒸騰作用達(dá)到降低地表溫度的作用，因而地表溫度較高。而綠地內(nèi)的植物和農(nóng)作物蒸騰蒸散較高，地表溫度較低。水域的比熱容較大，升溫較慢。

(3)馬蘭莊鎮(zhèn)的礦業(yè)在 2000 年以前處于開發(fā)不景氣時(shí)期，境內(nèi)的礦區(qū)占比較少僅為25.78%；2000年以后隨著礦業(yè)的大規(guī)模開發(fā)，采礦境界不斷擴(kuò)張，境內(nèi)礦區(qū)面積在2000年基礎(chǔ)上沿原開采境界不斷向外圍擴(kuò)展，致使2008年上升到38.78%；2008年以后隨著礦區(qū)生態(tài)修復(fù)工作的逐步實(shí)施，礦區(qū)占地面積逐步縮小，縮小的區(qū)域主要被植被覆蓋，使得2018年礦區(qū)面積下降到27.94%。

(4)馬蘭莊鎮(zhèn)境內(nèi)高溫區(qū)均位于礦業(yè)用地，并且高溫區(qū)面積占比從2000年的25.11%，上升到2008年的39.93%，再下降到2018年的28.71%，高溫區(qū)面積占比與礦業(yè)用地面積占比的變化趨勢(shì)基本一致，即均以2008年為拐點(diǎn)分為兩個(gè)階段，分別是2000—2008年的上升階段，和2008—2018年的下降階段。這是由于2000年開始礦業(yè)大規(guī)模開發(fā)，地表原有植被被破壞，變成了以裸地形式存在的礦業(yè)用地。裸地上由于沒有植被覆蓋，因此很難通過植被的蒸騰作用來達(dá)到降低地表溫度的目的，因此2008年地溫高溫區(qū)面積較之2000年明顯增加。2008年之后由于開始進(jìn)行生態(tài)修復(fù)工作，經(jīng)過人工栽植植被后，一部分裸地轉(zhuǎn)化為了植被覆蓋區(qū)，植被可以有效阻擋太陽(yáng)輻射直達(dá)地表，保持地表濕度，還可以通過蒸騰作用將吸收的太陽(yáng)輻射釋放出去，因而達(dá)到降低地表溫度的目的。