矿业密集区地表热环境时空异质性驱动机理

侯春华 ,李富平 ,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,说明地表温度的上升与裸土的增加关系最明显.