礦業密集區地表熱環境時空演變特征及驅動力分析

袁運許，朱小龍，侯春華，劉碩，宋建偉，李小光

(1.河北省地質環境監測院，石家莊 050000；2.河北省地質資源環境監測與保護重點實驗室，石家莊 050000；3.唐山師范學院計算機科學系，唐山 063000；4.唐山市資源與環境遙感重點實驗室，唐山 063210；5.華北理工大學礦業工程學院，唐山 063210；6.河北省礦業開發與安全技術重點實驗室，唐山 063210；7.河北省礦區生態修復產業技術研究院，唐山 063210)

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

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

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

1 研究方法

1.1 研究區概況及數據源

1.1.1 研究區概況

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

馬蘭莊鎮鎮內礦產資源豐富，采礦場分布密集，是鐵礦石資源豐富的鋼鐵重鎮。全鎮現有礦山企業超過50家。2001年遷安市境內19個鄉鎮，礦業用地密集程度的聚類分析統計結果表明，馬蘭莊鎮境內的礦業用地占比為34.80%，聚類分析結果表明該鄉鎮礦業開發密集程度為1級，選取馬蘭莊鎮為重點分析的鎮域尺度礦業密集區[17]。馬蘭莊鎮境內的礦業大多屬于露天開采方式，開采擾動的具體表現形式是破壞地表原有植被，導致大量巖石和土壤出露。

1.1.2 數據源

(1)遙感影像數據。選取與研究區礦業開發歷史變化緊密相關的5個典型時期的影像(表1)。分別代表礦業開發不景氣時期、礦業開始大規模開采初期、礦業破壞最嚴重時期同時也是礦區生態修復起始時期、礦區生態修復效果初顯時期、礦區生態修復效果較顯著時期。

表1 Landsat遙感影像數據

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

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

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

1.2 地表溫度反演及精度驗證

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

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

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

1.2.2 精度驗證

利用氣象站實測地溫進行精度驗證。由于馬蘭莊鎮沒有氣象站，因此利用遷安市氣象站數據進行驗證。由于遷安市氣象站缺少2000年、2003年和2008年實測地溫數據，因此只對2013年和2018年LST反演結果進行精度驗證。

2013年9月26日遷安市境內54439號氣象站，地溫實測點的經緯度分別為118°43′12″和 40°1′12″。2018年9月8日遷安市境內54439號氣象站，地溫實測點的經緯度分別為118°42′36″和 40°1′12″。兩個年份的實測點均位于Landsat影像的同一像元內。

經查驗該氣象站在2013年9月26日10：00點，0 cm和5 cm實測地溫分別為26 ℃和19.2 ℃。11：00，0 cm和5 cm實測地溫分別為29.2 ℃和22.3 ℃。2013年9月26日Landsat衛星在遷安市過境時間為10：49。對照谷歌高分影像發現該時刻氣象站地溫實測點位于稀疏草地像元內。因為地表溫度通常被定義為地球表面的溫度，在植被覆蓋度較低的區域，通常利用植被冠層以及土壤的平均溫度表征；在植被覆蓋度較高的區域，通常利用植被冠層的溫度表征；在裸地區域，通常利用土壤的表面溫度表征[20]。因此用0 cm(裸地地表)和5 cm(草地冠層)實測地溫均值來表示稀疏草地內的地溫值。

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

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

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

可知2013年和2018年地溫實測值與反演值之間誤差較小，反演結果精度較高。

1.3 礦業用地信息提取

1.3.1 影像分割和合并

基于ENVI5.3軟件平臺面向對象特征提取模塊 feature extraction的基于規則的面向對象影像分類算法，對馬蘭莊鎮境內礦業用地進行信息提取。

參考基于樣本的面向對象影像分類算法中分割和合并實驗方法，在實驗過程中通過預覽分割和合并效果確定最優分割尺度(圖2)。當分割尺度設置為10時，地物邊界清晰，但存在過分割現象，分割效果較破碎。當分割尺度設置為30時，地物邊界混雜，存在欠分割現象，不同地物可能會合并到同一地物類型中。最終經過多次實驗，當分割尺度為20，合并尺度為50時，各類地物邊緣特征能夠得到有效識別，因此作為最優閾值參與后續分類實驗。

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

1.3.2 特征分析

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

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

1.3.3 精度驗證

在谷歌地球高分辨率影像上，利用隨機選點的方式選擇礦業用地和其他地類各50個樣本點，對鄉鎮A境內礦業用地信息提取結果進行了精度驗證。表2為礦業用地信息提取混淆矩陣，總體分類精度和Kappa系數分別為93.00%和0.86。礦業用地的制圖精度和用戶精度分別為92.00%和93.88%，說明礦業用地信息提取結果精度較高，可用于下一步分析。

表2 礦業用地信息提取混淆矩陣

2 結果與分析

2.1 礦業用地時空動態演變分析

2.1.1 礦業用地空間變化分析

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

圖4 5期影像馬蘭莊鎮境內礦業用地分布

2.1.2 礦業用地時間演變統計

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

表3 5期影像馬蘭莊鎮境內礦業用地變化統計

為進一步定量化和可視化的探究馬蘭莊鎮境內，2000—2018年礦業用地時空動態變化過程，利用遙感變化檢測技術黃藍-差值影像法，分別獲得了5個時段(2000—2003年、2003—2008年、2008—2013年、2013—2018年、2000-2018年)，礦業用地和非礦業用地的變化影像(圖5)，并統計了礦業用地和非礦業用地面積占比變動情況。

黃色區域表示礦業用地增加區域(非礦業用地減少區域)；藍色區域表示礦業用地減少區域(非礦業用地增加區域)

2000—2008年，由于礦產資源的大規模開采，大量綠地被占用，致使礦業用地在2000—2003年和2003—2008年兩個時段內大幅度增加。2003年礦業用地占比相較于2000年增加了6.71%，2008年礦業用地占比相較于2003年增加了6.87%。

2008—2018年，礦區生態環境質量在深度治理中明顯改善，使得礦業用地在2008—2013年和2013—2018年這兩個時段內大幅度減少。2013年礦業用地占比相較于2008年減少了8.39%，2018年礦業用地占比相較于2013年減少了4.13%。雖然2008年以后政府開始重視環境保護工作，但是由于生態修復效果顯現時效滯后性的特點，2018年礦業用地占比仍高于2000年2.16%。基于“黃藍-差值影像法”得到的鄉鎮A境內5個時段的礦業用地變動情況，與鄉鎮A境內5個年份礦業用地的時空演變特征相吻合。

2.2 地表熱環境特征疊加分析

為了定量分析礦業開發對地表熱環境的擾動規律，對歸一化地表溫度(normalized land surface temperature, NLST)進行均值-標準差分級。根據NLST的取值范圍，依托ArcGIS10.2軟件，將NLST劃分為5個熱力等級(表4)，更加系統地反映研究區地表溫度等級分布狀況[22]。結果表明，5個年份境內地表熱環境均呈現顯著的空間異質性和規律性。林地和水域分別是低溫和次低溫區，耕地主要為中溫區，居民地是次高溫區，礦業用地為高溫區(圖6)。

圖6 5期影像鄉鎮A地表溫度等級分布

表4 NLST等級劃分標準

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

表5 5期影像鄉鎮A境內NLST各等級面積占比

2.3 地表熱環境特征剖面分析

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

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

圖7 鄉鎮A境內4條剖面線位置

為了使不同年份的地表溫度數據能在統一的量綱下進行對比，參考文獻[24]，首先對地表溫度進行正規化處理NLST，其計算公式[25]為

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

(1)

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

2.3.1 剖面線A-A溫度與地表下墊面覆被類型關系

剖面線A-A長度約為7 300 m，橫穿1個礦區采場、1個排土場和1尾礦庫(圖8)。

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

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

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

2.3.2 剖面線B-B溫度與地表下墊面覆被類型關系

剖面線B-B長度約為9 700 m，橫穿2個礦區采場、1個排土場和1個尾礦庫(圖9)。

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

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

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

2.3.3 剖面線C-C溫度與地表下墊面覆被類型關系

剖面線C-C長度約為7 100 m，橫穿1個尾礦庫和2個礦區采場(圖10)。

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

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

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

2.3.4 剖面線D-D溫度與地表下墊面覆被類型關系

剖面線D-D長度約為6 400 m，橫穿2個排土場和1個礦區采場(圖11)。

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

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

3 結論

在探究鎮域尺度礦業開發密集區，午間地表熱環境特征及驅動力過程中，基于輻射傳輸方程法(RTE)反演地表溫度，并進行定量精度驗證。基于“黃藍差值影像法”、疊加分析和剖面分析等方法，揭示礦業開采過程對地表熱環境動態變化的驅動力，彌補了礦業開采導致礦業密集區地表熱環境特征及驅動力分析的不足。得出以下主要結論。

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

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

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

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