礦業密集區地表熱環境時空異質性驅動機理

侯春華 ,李富平 ,3,4*,谷海紅 ,3,4,何寶杰 ,馬朋坤 ,宋 文 (1.華北理工大學礦業工程學院,河北 唐山063210；2.唐山市資源與環境遙感重點實驗室,河北 唐山 063210；3.河北省礦業開發與安全技術重點實驗室,河北 唐山063210；4.河北省礦區生態修復產業技術研究院,河北 唐山 063210；5.重慶大學建筑城規學院,重慶 400405；6.重慶大學,山地城鎮建設與新技術教育部重點實驗室,重慶 400405；7.澳大利亞新南威爾士大學建筑環境學院,澳大利亞 悉尼 2052；.北京城市氣象研究院城市邊界層與大氣環境研究所,北京 10009；9.中國科學院地理科學與資源研究所,陸地表層格局與模擬重點實驗室,北京100101)

礦業開發給國民經濟和社會發展作出貢獻的同時,會影響礦區及周邊環境和微氣候,但是目前人們對礦業密集區地表熱環境生物地球物理效應研究較少.礦業開采過程中大量剝離巖石和土壤,使得地表原有植被遭到破壞,壓占、堆積和采掘等擾動過程影響地表下墊面水分蒸發蒸騰循環[1],進而影響地表熱環境.地表溫度(Land Surface Temperature,LST)是區域與全球尺度地表物理過程分析與模擬的重要物理參量,在地表與大氣能量交換的過程中扮演著重要角色[2].當前主要集中在全球城市區域尺度范圍內,利用LST 開展城市空間熱環境及其氣候變化和環境整治研究,且多集中在熱島效應方面[3-4].

對于礦業密集區,無論是礦山企業數量還是礦山企業所占用土地數量,均相對集中[5].礦業開采改變了原有土地利用方式,必將對礦區及周邊區域地表熱環境造成一定影響,因此有必要分析礦業開發對地表熱環境的影響.區域地表熱環境定量評價是深入揭示地表熱環境形成機制的主要手段[6].在這方面衛星遙感數據具有定量估算時空尺度上地表生物地球物理參量的優勢[7].當前,通常利用遙感定量反演技術和GIS 空間分析方法,分析礦區地表熱環境和土地利用及其覆蓋變化(LUCC)、植被覆蓋度(VFC)等單一因素之間的定量關系,較少對人類活動強烈的礦業密集區地表熱環境異質性驅動因素進行綜合分析.如劉英等[8]以神東礦區為研究對象,基于1989～2013 年10 月 Landsat TM5/OLI 8 影像為數據源,采用單窗算法定量反演礦區地表溫度,并在礦區和礦井尺度上深入分析采區和非采區地表溫度的分布特征.朱星等[9]以永定礦區為研究區,選擇山地植被指數(NDMVI),利用像元二分模型方法,對礦區山地植被覆蓋度進行遙感估算,并利用熱紅外波段反演地表溫度來研究礦區的熱環境狀況變化.

礦業密集區內部景觀具有高度異質性,地表熱特性與眾多的生物物理特性密切相關,植被的生長狀況、蓋度和類型決定了地表反照率和水分蒸發蒸騰,從而影響了地表能量分配[10-11];地表濕度狀況決定了地表比熱特性,因而對LST 時空差異具有較好的揭示作用[12];礦業開采出現的大面積裸露地表由于可以改變地表層的潛熱和顯熱通量,因而也成為影響地表熱環境變化的關鍵因素[13-14].以上因素對礦業密集區地表能量交換和地表熱流均有影響,故單獨分析單一生物物理參數對地表熱環境的驅動機制是不夠的.

基于此,為揭示礦業密集區地表熱環境時空異質性特征及其驅動機制,本文利用2000～2018 年5 期Landsat 中分辨率遙感影像的熱紅外波段,采用大氣校正法反演研究區地表水熱因子中的LST,利用可見光波段結合研究區特點,提取地表下墊面擾動類型以及4 個表征地表熱環境時空異質性的生物物理參數;分別從地表下墊面擾動角度和生物物理角度分析地表熱環境時空異質性特征,借助統計學方法定量化分析不同生物物理參數對地表熱環境分異效應的驅動力,研究礦業密集區地表熱環境異質性的驅動機理及對不同生物物理參數的響應規律.研究結果對于準確了解礦業開發對區域地表熱環境的擾動規律、促進區域可持續發展等具有重要的理論和現實意義.

1 研究方法

1.1 研究區概況及數據源

1.1.1 研究區概況 研究區位于河北省遷安市西北部的馬蘭莊鎮,距市區15km.東臨灤河,與遷安鎮、閆家店鄉、大崔莊鎮三個鎮鄉僅一河之隔,南連蔡園鎮,西、北與遷西縣接壤(圖1).西部群山環抱,三面環山、一面鄰水,為半山區燕山隆起帶余脈南麓,總面積49.16km2.屬溫帶大陸性季風型氣候,四季分明,全年平均溫度10.1°C,常年降雨量735mm 左右,年均日照1675h,年均無霜期172d.

鎮內礦產資源豐富,采礦場分布密集,是鐵礦石資源豐富的鋼鐵重鎮.馬蘭莊鎮為遷安市礦業開發1級密集區[15],基于模糊綜合聚類分析方法得出馬蘭莊鎮境內的礦區占比和聚集程度居遷安市19 個鄉鎮之首,因此確定馬蘭莊鎮為本文重點分析的區域.

圖1 研究區位置示意Fig.1 Location of study area

1.1.2 數據源1)遙感影像數據與預處理 選取覆蓋研究區5 個時段的Landsat 影像,來源于美國地質調查局(USGS)官方網站(http://glovis.usgs.gov/)的Level 1T 級中分辨率系列衛星產品(軌道號122/32).數據獲取日期分別為:2000 年9 月6 日(TM)、2003年8 月14 日(TM)和2008 年9 月12 日(TM),2013年9 月26 日(OLI/TIRS)和2018 年9 月8 日(OLI/TIRS).經查驗 5 期影像過境時間,均為上午10:00～11:00,晴朗無云,空氣擾動較小,可以真實反映LST 狀況,適宜進行LST 反演.

影像預處理主要基于ENVI 5.3 軟件,對各時期遙感影像的可見光和近紅外波段進行輻射定標,將影像灰度值轉換為傳感器的反射率(鑒于Landsat 8TIRS 11 波段的定標參數偏差較大,因此對TIRS 10波段單獨進行輻射定標用于LST 反演[16].)使用FLAASH 大氣校正工具對各期影像的可見光波段進行大氣校正,消除大氣和光照等因素對地物反射的影響;最后利用研究區矢量邊界對預處理結果進行裁剪.

2)MODIS 溫度產品

利用NASA 官網提供的全球每日1000m 分辨率的MODIS 地溫產品MOD11A1 白天的溫度波段,對LST 反演結果進行對比驗證.

3)其他數據

谷歌地球高分辨率影像,用于對地表擾動類型劃分進行精度驗證;中國縣級行政矢量圖的研究區矢量邊界數據.

1.2 地表溫度反演

1.2.1 采用大氣校正法反演地表溫度 目前,利用熱紅外波段反演LST 的方法有很多,如大氣校正法(也稱為輻射傳輸方程:Radiative Transfer Equation,RTE)[17]、 單 窗 算 法(Mono-Window Algorithm,MW)[18]、劈窗算法(Split Window Algorithm,SW)[19-21]等,大量研究證明了這些算法的可行性[22-23].本文在參照相關研究的基礎上,采用大氣校正法反演LST,其基本原理為去除大氣影響前提下,借助大氣輻射傳輸方程,將衛星所觀測到地表熱輻射強度轉化為相應LST[24].詳細步驟見文獻[25].

1.2.2 地表溫度驗證 利用相對驗證法對遷安市LST 反演結果與較為成熟的被廣泛應用于各研究領域,并且經過大量學者驗證的MODIS 地溫產品MOD11A1 數據進行比較[26],由圖2 可知,二者在5期數據上的總體波動趨勢基本一致.5 個年份Landsat 反演LST 結果與MOD11A1 地溫產品的差值分別為:-0.54,-1.68,-2.57,3.26,1.19.二者差值較小,可滿足研究需要,并能夠反映研究區地表溫度變化情況,反演結果可用于下一步研究.

圖2 2000～2018 年反演地表溫度與當日MODIS 溫度產品對比Fig.2 Comparison of land surface temperature inversion from 2000 to 2018 with MODIS temperature products of the same day

1.3 地表熱環境異質性分析

為定量和可視化分析礦業密集區地表熱環境時空異質性特征,利用均值-標準差法對LST 進行分級,并分別從下墊面擾動類型角度和生物物理角度,對LST 和LUCC 以及4 個生物物理參數進行疊加分析和相關性分析,并分析午間LST 對下墊面擾動類型和生物物理參數時空變化的響應.

1.3.1 下墊面擾動類型分類 利用隨機森林分類法,將研究區下墊面地表擾動類型劃分為林地、耕地、工礦用地、建設用地和水域5 種類型,并對各地表擾動類型占比情況進行統計.精度驗證采用谷歌地球高分辨率影像和GF-1 影像,5 期影像總體分類精度和卡帕系數均高于90%和0.90,分類精度較高,效果較好.

1.3.2 地表溫度分級 為反映研究區LST 值數量及分布特征,選取合適的分級方法將反演所得溫度數據進行等級劃分.均值-標準差法可用于顯示要素屬性值與平均值之間的差異,使用與標準差成比例的等值范圍創建分類間隔,從而顯示不同溫度等級劃分情況與空間分布格局,能更好反映不同下墊面LST 差異[27-29].由于LST 會受到氣候等大背景環境因素年代間的不確定性和變化性的影響,直接利用LST 的絕對值進行地表熱環境的年代變化對比研究是不合理的,為了使不同年份的LST 數據能在統一的量綱下進行對比,參考前人研究[30],首先對LST 進行歸一化處理(NLST)[31],其計算公式如下:

式中:NLST 為歸一化后的地表溫度值;LST 為原始地表溫度值;LSTmax和LSTmin分別是LST 影像中的最小值和最大值.

根據NLST 的取值范圍,依托ArcGIS10.2 軟件,進一步將NLST 劃分為5 個熱力等級,分別為高溫、次高溫、中溫、次低溫和低溫(表1),更加系統地反映研究區LST 等級分布狀況[32].

表1 地表溫度等級劃分標準Table 1 Classification standard of land surface temperature

1.3.3 疊加分析 利用ArcGIS10.2 軟件的疊加分析功能,對研究區LST 和LUCC 進行疊加分析,統計不同下墊面擾動類型的LST 均值,據此從下墊面擾動角度探究研究區地表熱環境異質性驅動機制.

1.3.4 空間相關性分析 基于可見光波段分別提取4 個生物物理參數,從生物物理角度對研究區地表下墊面4 個生物物理參數和NLST 進行空間相關性分析,明晰研究區地表熱環境空間異質性對下墊面生物物理參數的響應規律.

1.4 生物物理參數

地表綠化程度、土壤濕度、地表裸露程度和不透水面密集程度等生物物理指標均能夠反映地表生態狀態,并且對LST 具有一定的正面或負面影響[33-35].

1.4.1 光合植被覆蓋度(fPV) 以往研究多集中于利用僅考慮裸土(Bare Soil,BS)和植被2 種組分的像元二分模型(Dimidiate Pixel Model,DPM)[36],反演區域VFC 探討其與地表熱環境之間的關系.植被的降溫作用一定程度上是因為光合植被(Photosynthetic Vegetation,PV)葉面的水分含量,如Amiri 等[37].在伊朗的研究表明,植被含水量變化導致高植被覆蓋區從低溫特征向高溫特征轉變,因而如何掌握光合植被覆蓋信息,對緩解區域地表熱環境高溫聚集效應具有重要意義[37].基于NDVI-DFI 特征空間的像元三分模型反演光合植被覆蓋度(fPV),定量分析其對LST 的驅動作用[38],具體計算步驟見文獻[39].

1.4.2 土壤濕度監測指數(SMMI) 地表水分狀況作為陸面過程的重要因子,是地表干旱信息最重要的表征參量,同時也是調控地-氣反饋的最重要參量之一[40],因而地表水分含量的變化對區域地表熱環境調控尤為重要.基于土壤濕度監測指數SMMI 表征研究區地表水分含量,該指數基于Nir－Red 光譜特征空間建立,具有不依賴土壤線而變化的優點[41],表達式為:

式中:Rred和Rnir分別為紅光波段和近紅外波段的光譜反射率.為使反演結果中像元值大的區域代表地表水分含量高,用1減去SMMI的反演值得到研究區地表水分含量影像.

1.4.3 增強型裸土指數(EBSI) 采礦活動的實施直接剝離掉大量礦區采場表土和巖石,由于裸土和巖石具有相似的反射率,出現的大面積裸露地表在遙感影像上與裸土的光譜特征相一致.采用增強型裸土指數EBSI 來量化研究區地表裸露程度[42].計算公式:

式中:Rswir1、Rred、Rblue和Rnir分別代表短波紅外1波段,紅波段,藍波段和近紅外波段的光譜反射率.

1.4.4 歸一化差值不透水面指數(NDISI) 地表非滲透表面具有較高的LST,一方面是因為其自身的物理特性,另一方面是因為不透水表面水分含量低,缺乏有效的蒸騰作用,因此有必要將不透水面也作為影響地表熱環境分異效應的一個重要因素進行研究[43].采用歸一化差值不透水面指數NDISI 量化研究區不透水面信息.公式為:

式中:Rmir、Rnir和Rswir2分別代表中紅外波段、近紅外波段和短波紅外2 波段的光譜反射率.MNDWI 為歸一化水體指數,公式為:

式中:Rgreen和Rmir分別代表綠波段和中紅外波段的光譜反射率.

2 結果與分析

2.1 地表熱環境時空異質性特征分析

2.1.1 下墊面擾動類型角度疊加分析 由于下墊面擾動類型發生變化會引起地表反射率的變化,從而導致LST 發生變化,為定量化、可視化的分析馬蘭莊鎮地表擾動類型的溫度分異特征,以及對研究區增溫和降溫的貢獻程度,從研究區下墊面擾動角度對5 期影像LST 和LUCC 進行疊加分析,對研究區進行下墊面擾動類型分類,并統計各年份地表擾動類型占比(圖3).5 期影像林地主要集中于中西部山區;耕地分散于境內中東部地區;工礦用地在研究區北部和中南部較為密集;居民地分散在工礦用地周邊;水域主要是位于各礦區尾礦庫內的少量尾水.2000 年礦業處于不景氣時期,境內的礦區較少,對生態破壞不大,綠地占比較大,工礦用地占比僅為26.56%;2000 年以后隨著礦業的大規模開采,大量綠地和耕地被占用,生態環境破壞程度加劇,工礦用地占比逐年增加,到2008 年綠地占比減少到11.54%,工礦用地占比增加到55.54%;2008 年以后,由于政府開始重視環境保護工作,礦區生態恢復開始逐步實施,工礦用地占比逐年縮小,但是由于生態復墾效果顯現時效滯后性的特點,到2018 年綠地占比僅比2008 年提高了2.13%.

為探究馬蘭莊鎮地表熱環境分異效應驅動機制,從地表下墊面擾動角度對LUCC 和LST 反演結果進行疊加分析,5 個年份境內不同下墊面擾動類型LST,均表現出工礦用地＞居民地＞耕地＞林地＞水域的分異特征(表2).

為進一步分析地表下墊面擾動對LST 的驅動規律,對NLST 進行均值-標準差分級,結果表明5 個年份境內地表熱環境均呈現顯著的空間異質性和規律性,低溫區和次低溫區主要分布于林地和水域,中溫區主要位于耕地,次高溫區主要分布于居民地,高溫區主要聚集在位于研究區北部和中南部的露 天礦區采場境內(圖4).

圖3 2000～2018 年馬蘭莊鎮下墊面擾動類型分類及占比(%)Fig.3 Land use classification and proportion of land disturbed types of Malanzhuang town from 2000 to 2018

表2 2000～2018 年馬蘭莊鎮各地物地表溫度均值(°C)Table 2 The mean value of LST from 2000 to 2018

這是由于在接收同等太陽輻射能量的條件下,裸土和不透水面比熱容小于水體,導致裸土和不透水面升溫更快.居民區和公共設施用地雖然為不透水面等建設用地,但溫度未及工礦用地,原因之一是居民區和公共設施用地內部植被較之豐富,土壤水分含量及蒸散量大于工礦用地,另一方面工礦用地因其裸露地表分布廣、含水量低以及地表反照率高使其在接收同等太陽輻射量的前提下,能吸收更多的熱量而具有極高的LST,形成局部熱島中心.可據此從地表下墊面擾動類型角度初步判斷馬蘭莊鎮地表熱環境分異效應的驅動機制.

2.1.2 生物物理角度空間相關性分析 為從多角度剖析研究區地表熱環境異質性特征及驅動機制,從生物物理角度分別對研究區下墊面4 個生物物理參數和NLST 進行空間相關性分析.為了使不同年份的4 個生物物理參數和NLST 能在統一的量綱下進行對比,對各年份的值進行歸一化處理,使其值域分別在[0,1]范圍內. 以2018 年為例,結合下墊面擾動類型分類及占比圖(圖3e)和地表熱力等級空間分布圖(4e),分析研究區地表熱環境空間異質性對下墊面生物物理參數的響應規律(圖5).

圖5a 和圖5b 的高值區均對應于林地、水域等低溫區,低值區對應于工礦用地和不透水面等高溫區;圖5c 和圖5d 的高值區均對應于工礦用地和不透水面等高溫區,低值區對應于林地、水域等低溫區.上述研究表明,地表生態狀況對礦業密集區地表熱環境分布的影響較為明顯,4 個生物物理參數和LST均存在著較強的相關性.這是由于不同地表覆被對于直接太陽輻射和間接太陽輻射的反射與散射導致地表熱環境的異質性.

圖4 馬蘭莊鎮2000～2018 年地表熱力等級空間分布Fig.4 Spatial distribution of NLST in study area from 2000 to 2018

圖5 2018 年馬蘭莊鎮4 個生物物理參數Fig.5 The four factors’ patterns over the study area

為進一步探究4 個生物物理參數與LST 之間的關系,以2018 年為例,利用ENVI5.3 對NLST 和4 個生態參數分別作散點圖分析(圖6).由散點圖可知,fPV和SMMI 與LST 均表現出明顯的負相關關系,EBSI 和NDISI 與LST 均表現出明顯的正相 關關系.

圖6 2018 年馬蘭莊鎮LST 與4 個生物物理參數散點圖Fig.6 The scatter plot of LST and four factors in study area

2.2 四個生物物理參數對NLST 的驅動規律

為進一步了解研究區地表熱環境異質性對以上4 個生物物理指標時空變化的響應規律,分別探討2000～2008 和2008～2018 年兩個時段4 個生物物理指標對地表熱環境分異效應的驅動規律(表3).

表3 2000～2018 年4 個生物物理指標統計值Table 3 The mean value of four factors from 2000 to 2018

由表3 可知,2000～2018 年,與LST 表現出明顯負相關關系的fPV和SMMI 分別從2000 年的0.4951 和0.7340 下降到2018 年的0.3542 和0.7263,與LST 表現出明顯正相關關系的EBSI 和NDISI 分別從2000年的0.6132 和0.4769 上升到2018 年的0.7093 和0.6125,并且2000 年、2008 年和2018 年的NLST 均值分別為0.4131、0.4303 和0.4861,由NLST 均值統計結果表明研究區LST 在19a 間整體呈現上升趨勢.這是由于2000 年到2008 年,fPV和SMMI 總體呈現明顯下降趨勢, NDISI 和EBSI 總體呈現顯著上升態勢,在這期間礦業大規模擴張,改變了原有的下墊面景觀格局,植被遭到破壞,原有土地資源銳減,這些變化降低了植被覆蓋和地表蒸散,增加了不透水材料對太陽熱輻射的吸收,使得地表輻射能量平衡發生顯著變化. 2008 年以后大力開展了綠色礦山建設,礦區生態恢復工作開始實施,2008 年到2018 年,fPV和SMMI 呈現上升趨勢,植被覆蓋度提高加大了地表蒸散,但是由于生態復墾效果顯現時效滯后性的特點,fPV和SMMI 較2008 年表現出微弱增長態勢.

2.3 NLST 與下墊面生物物理參數回歸分析

2.3.1 單因素回歸分析 為進一步定量分析研究區地表熱環境分異效應對下墊面4 個生物物理參數的響應規律.利用SPSS22.0 軟件分別將4 個生物物理參數作為自變量,NLST 作為因變量進行回歸分析(表4).結果表明,fPV和SMMI 與NLST 在0.01 水平(雙側)均呈顯著線性負相關關系(P＜0.01);EBSI 和NDISI 與NLST 在0.01 水平(雙側)均呈顯著線性正相關(P＜0.01).由回歸系數可知,19a 間,fPV和SMMI每增加10%,會使NLST 相應下降0.0520～0.1060 和0.0550～0.1130.EBSI 和NDISI 每增加10%,會使NLST 相應上升0.0500～0.1040 和0.0570～0.0960,因此研究區地表溫度的下降與光合植被覆蓋度和土壤水分含量的升高有明顯的關系,地表溫度的上升與裸土和不透水面的增加有明顯的關系.

表4 2000～2018 年4 個生物物理指標與NLST 單因素回歸分析Table 4 Single factor regression impacts in study area on the relationship between four biophysical parameters and NLST from 2000 to 2018

表5 2000～2018 年4 個生物物理指標與NLST 多因素回歸分析Table 5 Multiple regression results in study area based on four biophysical parameters and NLST throughout 2000 to 2018

2.3.2 多因素綜合分析 為進一步了解4 個生物物理參數對研究區地表熱環境的綜合驅動作用,對4個生物物理參數與NLST 進行多元線性回歸分析(表5).各年份多元回歸方程均通過了皮爾遜相關系數的顯著性P＜0.01 的顯著性檢驗.并且各年份多元回歸方程的決定系數R2均大于一元回歸方程,表明利用4 個生物物理參數綜合考量地表熱環境異質性,比單獨利用單一參數效果更好.另外,以2018 年fPV為例,在單因素回歸模型中,fPV每增加10%,會使NLST 相應下降0.0610,而在綜合多因素的回歸模型中,fPV每增加10%,會使NLST 相應下降0.0150,兩者相差4 倍.可見單因素分析并不能全面反映下墊面生物物理參數與LST 的關系,并且可能高估各個對象之間真實的回歸系數.從方程回歸系數來看,5 期影像與LST 呈負相關關系的SMMI 回歸系數均大于fPV,說明地表溫度的下降與土壤濕度的增加關系最明顯,與LST 呈正相關關系的EBSI 回歸系數均大于NDISI,說明地表溫度的上升與裸土的增加關系最明顯.

3 結論

3.1 利用Landsat 遙感數據基于大氣校正法反演LST,反演結果與地面實測氣溫數據變化趨勢基本一致,反演結果誤差較小,精度較高,表明利用遙感數據在無地面同步資料情況下進行溫度反演,能夠滿足地表熱環境異質性研究的需要.

3.2 從下墊面擾動角度分析得知,研究區地表熱環境的時空異質性特征非常明顯.5 期影像礦業用地LST 均居首位,屬于熱島效應聚集區,且5 期影像的LST 均表現出工礦用地＞居民地＞耕地＞林地＞水域的異質性特征.從生物物理角度分析得知,4 個生物物理參數和LST 存在著較強的相關性,并且fPV和SMMI 與LST 表現出明顯的負相關關系,EBSI 和NDISI 與LST 表現出明顯的正相關關系.

3.3 從單因素回歸分析結果來看,4 個生物物理參數與NLST 均通過了P＜0.01 的顯著性檢驗,4 個驅動因子對研究區LST 均存在顯著影響;其中fPV和SMMI 與NLST 在0.01 水平(雙側)均呈顯著線性負相關關系,EBSI和NDISI與NLST在0.01水平(雙側)均呈顯著線性正相關關系.

3.4 從多因素回歸分析結果來看,利用4 個生物物理參數綜合考量地表熱環境異質性,比單獨利用單一參數效果更好.從方程回歸系數來看,地表溫度的下降與土壤濕度的增加關系最明顯;與LST 呈正相關關系的EBSI 回歸系數均大于NDISI,說明地表溫度的上升與裸土的增加關系最明顯.