地表反照率和植被覆盖度对矿区热环境的影响

侯春华，李富平，何宝杰，马朋坤，宋 文,6，赵菁菁

(1.华北理工大学 矿业工程学院，河北 唐山 063210;2.唐山师范学院 计算机科学系，河北 唐山 063000；3.重庆大学 建筑城规学院,重庆 400405;4.澳大利亚新南威尔士大学 建筑环境学院，澳大利亚 悉尼 2052;5.北京城市气象研究院城市边界层与大气环境研究所，北京 100089;6.中国科学院地理科学与资源研究所，北京100101)

地表温度(Land surface temperature，LST)是表征大气底层和陆地表层之间热量状况的重要指标，其变化主要受地表反照率(Albedo)和下垫面性质等因素的影响[1]。地表反照率是描述地表辐射特征的量纲一的参数[2]，时空变化受到自然过程以及人类活动的影响，例如工业化、城市化以及季节变化等都会对地表反照率造成影响，进而影响地表温度[3-4]。Hansen等[5]研究得出全球土地利用变化会导致地表反照率发生变化，从而引起全球净辐射减少，使得气温下降。矿产资源开采活动改变了原有地表下垫面性质，大量植被遭到破坏，使得地表反照率发生变化，进而影响地表辐射平衡和热量状况[6]。

以往研究只考虑植被覆盖的动态变化对地表热环境的驱动作用[7-9]。研究表明，植被的降温作用一定程度上是因为光合植被(Photosynthetic vegetation，PV)叶面的水分含量，如Amiri等[10]在伊朗的研究表明，高植被覆盖区因植被含水量变化导致地温在低温和高温特征之间相互转变。因而提取光合植被覆盖度(fractional cover of photosynthetic vegetation，fPV)信息，分析其动态变化引起地表反照率发生改变进而对地表热环境的驱动作用具有重要意义。近年来，地表反照率的时空特征研究一直备受关注，这些研究大多都是大、中尺度区域地表反照率时空分布变化特征[11-13]，对于矿业开采导致矿区内部光合植被覆盖改变引起地表反照率时空分布及其特征发生变化，进而加剧地表热环境异质性的分析研究相对薄弱。

基于此，本文基于2000年、2008年和2018年3期Landsat遥感影像的可见光波段，反演光合植被覆盖度与地表反照率，基于热红外波段反演地表温度。利用叠加分析、相关分析和回归分析法，定量化和可视化分析了fPV和Albedo与地表温度的相关性，进而探究由于矿产开发引起矿区光合植被覆盖时空变化导致地表反照率变化之间的潜在关系，并揭示矿区地表热环境异质性对光合植被覆盖度和地表反照率变化的响应规律。

1 研究方法

1.1 研究区概况及数据源

1.1.1 研究区概况

唐山首钢马兰庄铁矿位于河北省唐山市东北部的迁安市马兰庄镇境内(如图1所示)。迁安市以山地和丘陵为主，地势西北高，东南低，属暖温带、半湿润大陆性季风气候，全年气温变化较大。马兰庄铁矿矿区面积2.9 km2，于1970年开始建矿，1972年正式投产。矿山目前采用露天开采方式，主要包括白马山和沙河山2个采场，白马山采场于2007年末闭坑。矿区东临滦河、北依燕山、西与迁西接壤、南与首钢矿区相连，矿藏丰富。矿山于2008年、2009年和2010年分别在尾矿库、柳河峪排土场和白马山排土场持续实施植被恢复工程，主要种植刺槐、火炬树、沙棘等适宜于矿区自然条件的先锋植物和适宜树种[14]。

图1 研究区位置与三维示意图

1.1.2 数据源及预处理

1)遥感影像数据。所用数据分别为2000年9月6日、2008年9月12日和2018年9月8日，来源于美国地质调查局(USGS)Landsat官方网站(http://glovis.usgs.gov/)，Level 1T级中分辨率系列卫星产品(轨道号122/32)，分辨率为30 m，云量均小于2%，经查验 3 期影像过境时间均为上午10:00～11:00之间，卫星过境全天均无降水，天气状况较好，适宜进行地表反照率和LST反演。来源于地理空间数据云(http://www.gscloud.cn/)的，分辨率为250 m的与Landsat数据获取日期一致的MODIS地表反射率产品MOD09GQ，用于对Landsat反演的地表反照率进行验证。

预处理工作依托ENVI5.3软件平台，对各期Landsat影像进行辐射定标，大气校正和裁剪。各期MOD09GQ产品依托ArcGIS10.2软件进行投影变换使之与Landsat投影保持一致，最后拼接和裁剪。

2)实地调查数据。于2018年9月9日派出具有野外实地工作经验，并具有一定地理学理论基础的实测人员5人，采用实地测量法定量验证植被覆盖度反演精度。在矿区内选取10个有代表性的30 m×30 m的网格，网格中心点周围的植被覆盖度均存在一定差异，尽量包含0～100%的植被覆盖范围，分别记录下各网格单元中心点的经、纬度坐标。每个网格由5位观测员分别轮流站在网格单元的中心点进行目视观测，并分别将观测结果记录下来，观测结束后对每个网格的植被覆盖度观测结果求取均值。

3)其他数据。中国县级行政矢量图的研究区矢量边界数据；2018年9月28日空间分辨率为全色2 m、多光谱8 m的GF-1卫星遥感影像；谷歌地球高分辨率影像，用于辅助目视解译矿区地表扰动类型。

1.1.3 地表下垫面扰动类型提取

在参考中国土地利用现状调查分析体系和国土资源部土地利用分类体系的基础上，经过实地调研结合谷歌地球高分辨率影像，考虑迁安市首钢马兰庄铁矿地表下垫面覆被状况，利用目视解译法提取境内沙河山采场、尾矿库、柳河峪排土场、白马山排土场和水体等几类对地表热环境异质性影响较大的典型地物边界(如图2所示)。

图2 2018年研究区主要典型地物边界图

1.2 地表温度反演及等级划分

1.2.1 地表温度反演

利用辐射传输方程法(Radiative transfer equation，RTE)反演地表温度，基本原理为估计大气对地表热辐射的影响，将这一部分从卫星传感器所观测到的热辐射总量中减去，得到地表热辐射强度，然后借助大气辐射传输方程，将卫星所观测到地表热辐射强度转化为相应地表温度[15]。

卫星传感器接收到的热红外辐射亮度值的辐射传输方程为

Lλ=[εB(T)+(1-ε)Ld]τ+Lu

(1)

式中：Lλ为卫星热红外波段的辐射亮度，ε为地表比辐射率，T为地表真实温度，B(T)为黑体热辐射亮度)，τ为大气在热红外波段的透过率，Lu、Ld分别为大气上行和下行辐射亮度[16]。

假设大气、地表对热辐射具有朗伯体性质，则温度为T的黑体在热红外波段的辐射亮度B(T)为

B(T)=[Lλ-Lu-τ(1-ε)Ld]/τε

(2)

地表比辐射率计算公式为

ε=0.004Pv+0.986

(3)

通过覃志豪等[17]提出的方法计算地表比辐射率ε。地面真实温度T根据普朗克公式为

T=K2/ln(K1/B(T)+1)

(4)

式中：K1、K2为热红外波段的定标参数，对于Landsat5 TM影像第6波段，K1=607.76 W/(m2·μm·sr)，K2=1 260.56 K；对于Landsat8 TIRS影像第10波段，K1=774.89 W/(m2·μm·sr)，K2=1 321.08 K[18]。大气校正参数τ、Lu和Ld在美国国家航空航天局NASA官网(http://atmcorr.gsfc.nasa.gov/)中输入影像中心经纬度、成像时间和其他参数信息获取。

1.2.2 地表温度分级

为反映研究区地表温度值等级及分布特征，利用密度分割法将反演所得地表温度数据进行等级划分。为使不同年份的地表温度数据能在统一的量纲下进行对比，并消除气候等大背景环境因素年代间的不确定性和变化性的影响，参考已有研究[19]，首先对地表温度进行正规化处理(Normalized land surface temperature,NLST)[20]，其计算公式如下：

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

(5)

式中：NLST为正规化后的地表温度值，LST为原始地表温度值，LSTmax、LSTmin分别为LST影像中的最大值和最小值。

根据NLST的取值范围，依托ArcGIS10.2软件，以0.2为间隔将NLST划分为5个热力等级，分别为低温、次低温、中温、次高温和高温(见表1)，更加系统地反映研究区地表温度等级分布状况。

表1 地表温度等级划分标准

1.3 光合植被覆盖度反演及验证

1.3.1 光合植被覆盖度反演

通过构建NDVI-DFI特征空间提取纯净端元特征值，利用像元三分模型对研究区混合像元进行分解，分解为光合植被覆盖度、非光合植被覆盖度(fractional cover of NPV，fNPV)和裸土覆盖度，单独保存各波段以此来获取光合植被覆盖度信息。经实地采样验证，像元三分模型在研究区fPV的估算精度较高。具体计算步骤详见文献[21]。为了对不同年份的fPV在统一的量纲下进行对比，利用以下公式对fPV进行正规化处理，使其值域在[0，1]之间：

NfPV=(fPV-(fPV)min)/((fPV)max-(fPV)min)

(6)

式中：NfPV为正规化后的光合植被覆盖度,fPV为原始光合植被覆盖度,(fPV)max、(fPV)min分别为原始光合植被覆盖度的最大值和最小值。

1.3.2 验 证

分别利用定性和定量方法进行验证。首先利用GF-1假彩色影像对植被覆盖度(Vegetation fractional coverage，VFC)反演结果进行了定性验证。绘制了研究区2018年9月8日NDVI-DFI特征空间散点图，并截取同期研究区内典型地物fPV、fNPV和fBS的RGB空间分布图，与近同期2018年9月28日GF-1假彩色影像上的真实地物做对比(如图3所示)。图3(a)为非光合植被反演结果；图3(b)为对应GF-1影像中的非光合植被；图3(c)为裸土反演结果；图3(d)为对应GF-1影像中的裸土；图3(i)为光合植被反演结果；图3(j)为对应GF-1影像中的光合植被；图3(k)、3(g)为非光合植被伴生于光合植被中的反演结果；图3(l)、3(h)为对应的GF-1影像；图3(e)为混合光合/非光合植被和裸土的反演结果；图3(f)为对应的GF-1影像。结果表明和实际地物吻合度较高，反演效果较好。

图3 利用GF-1影像对基于NDVI-DFI像元三分模型反演VFC的结果作定性验证

然后按照经、纬度一致原则，提取实地采样点相应位置像元fPV、fNPV和fBS的均值，对实测值和反演值进行线性拟合进行定量验证(如图4所示)。拟合结果表明，fPV、fNPV和fBS的实测值和反演值的决定系数R2分别为0.78、0.70和0.76，在p<0.01水平上显著相关，拟合度较好，反演精度较高。

图4 fPV、fNPV和fBS的反演值和实测值

1.4 地表反照率(Albedo)反演及验证

1.4.1 地表反照率(Albedo)反演

地表反照率反演算法参考Liang[22]于2011年建立的，适用于不同大气和下垫面类型的公式：

Albedo=0.356×Rblue+0.130×Rred+

0.373×Rnir+0.085×Rswir1+

0.072×Rswir2-0.001 8

(7)

式中：Rblue、Rred、Rnir、Rswir1和Rswir2分别为传感器对应的蓝、红、近红外、短波红外1和短波红外2波段的光谱反射率，Albedo为地表反照率。

为了对不同年份的Albedo在统一的量纲下进行对比，利用以下公式对Albedo行正规化处理，使其值域在[0，1]之间：

NAlbedo=(Albedo-(Albedo)min)/

((Albedo)max-(Albedo)min)

(8)

式中：NAlbedo为正规化后的地表反照率,Albedo为原始地表反照率,(Albedo)max、(Albedo)min分别为原始地表反照率的最大值和最小值。

1.4.2 地表反照率验证

分别利用定性和定量方法进行验证。首先，以2000年分辨率为250 m的的MODIS的MOD09GQ陆地反射率产品为数据源，对以同年分辨率为30 m的基于Landsat影像反演的Albedo结果的空间一致性进行相对验证。由于MOD09GQ空间分辨率较低，矿区尺度较小，以矿区范围进行对比不能很好的观察二者空间对应情况。因此以整个迁安市2000年影像为例，对反演结果进行相对验证，可以看出二者的高值区和低值区的对应效果较好，并且Landsat反演结果的空间分辨率更高，更适合对小尺度矿区境内的地表反照率进行分析(如图5所示)。

然后利用定量验证方法，从研究区境内随机选取20个采样点，对基于Landsat数据反演的Albedo结果与MOD09GQ地表反射率产品进行线性拟合(如图6所示)。可以看出二者具有较高的相关性，决定系数R2达0.86，说明研究区适宜利用该方法反演Albedo，反演结果精度较高。

图5 2000年Landsat反演地表反照率与MOD09GQ反射率产品

图6 基于Landsat的Albedo反演值和MOD09GQ反射率产品

2 结果与讨论

2.1 基于fPV和Albedo的地表热环境异质性特征

基于3期Landsat影像反演的地表温度结果，均呈现明显的空间异质性和规律性，说明Landsat数据的空间分辨率能达到该研究区尺度下的地表温度反演要求，满足地表热环境的空间分布研究(如图7所示)。对照图1可知，高温区主要分布于北部沙河山采场、中部柳河峪排土场和南部未复垦的白马山排土场，中温区主要位于高温区边缘以及其他裸露地表，低温区主要位于北部已复垦尾矿库和已复垦排土场以及境内绿地和水体。

为进一步明晰矿区地表热环境异质性特征，利用ArcGIS10.2的叠加分析功能，分别计算矿区各地物地温均值(见表2)。结果表明，2000年地温排序为：白马山排土场>沙河山采场>柳河峪排土场>尾矿库>水域，2008年地温排序为：白马山排土场>柳河峪排土场>沙河山采场>尾矿库>水域；2018年地温均值排序为：沙河山采场>柳河峪排土场>白马山排土场>尾矿库>水域。19年间白马山排土场的地温排序变动情况可知，2000—2008年矿区处于大规模开采期，白马山排土场地温在这9年间始终居高不下；2008年开始随着生态复垦工作的实施，白马山排土场开始栽植大量灌木和乔木等耐旱耐高温植被，使得2018年该区域地温均值低于沙河山采场和柳河峪排土场，说明2008—2018年间，白马山排土场的生态修复工作效果较好，随着植被覆盖和地表湿度的明显上升，地表反照率下降，地表温度明显下降，此外，还应考虑到不同地表覆被对于直接太阳辐射和间接太阳辐射的反射与散射问题导致的地表热环境异质性问题。

图7 2000—2018年NLST空间分布

表2 2000—2018年马兰庄矿各地物地表温度均值

2.2 fPV和Albedo与地表温度相关性分析

为探究fPV和Albedo与地表温度的相关关系，利用ENVI5.3软件的密度分割工具，分别对fPV、Albedo和NLST进行等级划分，划分规则为每隔0.2为一个等级(如图8所示)，并利用ArcGIS10.2软件的统计工具计算fPV、Albedo和NLST各级别面积占比(见表3～5)。

图8 2000—2018年fPV、Albedo和NLST等级占比

表3 2000—2018年fPV各等级占比

表4 2000—2018年Albedo各等级占比

表5 2000—2018年NLST各等级占比

2000—2008年间，由于矿业大规模开采，植被遭到破坏，fPV的高、次高和中等级别面积总占比由2000年的49.19%，下降到2008年的23.67%，下降了25.52个百分点。随着fPV的下降，采场裸岩逐渐出露，和矿区采场、未复垦排土场和未复垦尾矿库等裸露地表相重叠的高温、次高温和中温区，面积总占比由2000年的78.8%上升到2008年的82.51%，增加了3.71个百分点。2008—2018年间，随着矿区生态修复工作着手实施，尾矿库和排土场等区域栽植了大量植被和作物，fPV由2008年的23.67%，上升到2018年的43.81%，上升了20.14个百分点。随着fPV的上升，高温、次高温和中温区，面积总占比由2008年的82.51%下降到2018年的68.03%，下降了14.48个百分点。2000—2018年19年间，高温、次高温和中温区面积占比下降了10.77个百分点，Albedo等级高、次高和中级面积占比下降了25.52个百分点，说明生态修复工作成效显著。

为进一步探究fPV和Albedo与地表温度之间的关系，以2018年为例，依托ENVI5.3软件，分别利用fPV和Albedo与地表温度做散点图(如图9所示)，结果表明fPV与LST呈现明显的负相关关系，Albedo与LST呈现明显的正相关关系。可据此初步认为，fPV的增加有助于缓解地表热强度，Albedo的增大有助于增强地表热强度。

2.3 光合植被覆盖度对地表热环境异质性的驱动作用

2000—2018年fPV空间分布具有明显的异质性和规律性(如图10所示)，高值区分布于已复垦排土场和尾矿库以及其他植被覆盖区域，中低值区位于矿区采场裸岩和排土场裸土等裸露区域和水体。利用ArcGIS10.2的叠加分析功能，计算3期影像不同下垫面LST和fPV均值并排序：矿区采场裸岩LST均值35.29℃>排土场矿渣LST均值35.02℃>复垦植被LST均值25.25 ℃>水体LST均值21.77℃；复垦植被fPV均值0.73>排土场矿渣fPV均值0.09>矿区采场裸岩fPV均值0.03>水体fPV均值0.00。

图9 2018年马兰庄矿地表温度与fPV和Albedo散点图

图10 2000—2018年光合植被覆盖度空间分布

为定量分析地表热环境异质性对fPV变化的响应规律。依托ArcGIS10.2软件提取随机点工具，以及SPSS22.0软件的相关分析功能对3期影像的fPV和LST反演值进行回归分析(如图11所示)。结果表明，2000、2008、2018年决定系数R2分别为0.63、0.55和0.65，fPV与LST在0.01水平(双侧)呈线性负相关关系(p<0.01)，说明fPV的增加对地表热环境具有降温效应，由回归系数可知，2000、2008、2018年，fPV每增加10%，会使LST相应下降0.66、0.74、1.09 ℃。

图11 2000—2018年fPV与LST回归分析

2.4 地表反照率对地表热环境异质性的驱动作用

2000—2018年Albedo空间分布具有明显的异质性和规律性(如图12所示)，高值区分布于矿区采场裸岩和排土场裸土等裸露区域，中低值区位于植被和水体。利用ArcGIS10.2的叠加分析功能，计算3期影像不同下垫面LST和Albedo均值及排序：矿区采场裸岩LST均值35.29 ℃>排土场矿渣LST均值35.02 ℃>复垦植被LST均值25.25 ℃>水体LST均值21.77 ℃；矿区采场裸岩Albedo均值0.79>排土场矿渣Albedo均值0.68>复垦植被Albedo均值0.46>水体Albedo均值0.12。

为定量分析地表热环境异质性对Albedo变化的响应规律。依托ArcGIS10.2软件提取随机点工具以及SPSS22.0软件的相关分析功能对3期影像的Albedo和LST反演值进行回归分析(如图13所示)。结果表明，2000、2008、2018年决定系数R2分别为0.35、0.40、0.48，Albedo与LST在0.01水平(双侧)呈线性正相关关系(p<0.01)，说明Albedo的增加对地表热环境具有增温效应，但是相较于fPV，Albedo和地表温度的相关性较小。由回归系数可知，2000、2008、2018年，Albedo每增加10%，会使LST相应上升1.09、1.36、1.76 ℃。此外，不同地表覆被对于直接太阳辐射和间接太阳辐射的反射与散射可能对该结果造成一定的干扰。

图12 2000—2018年地表反照率空间分布

图13 2000—2018年Albedo与LST回归分析

3 结 论

1)研究区地表温度空间分布呈现明显的异质性和规律性，高温区主要分布于北部沙河山采场、中部柳河峪排土场和南部未复垦的白马山排土场，中温区主要位于高温区边缘以及其他裸露地表，低温区主要位于北部已复垦尾矿库和已复垦排土场以及境内绿地和水体。

2)研究区光合植被覆盖度在2000—2018年19年间，呈现先下降再上升的变化趋势，相应的高温、次高温和中温区呈现先上升后下降的变化趋势，说明光合植被覆盖度与地表温度呈现负相关关系，并且空间上呈现相反的变化特征。研究区地表反照率在2000—2018年19年间，呈现先微弱下降再大幅度下降趋势，总体呈现下降趋势，与地表温度19年间高温、次高温和中温区的总体呈现下降趋势的变化特征一致，说明地表反照率与地表温度呈现正相关关系，并且空间上整体变化趋势相一致。

3)光合植被覆盖度和地表温度的回归分析结果表明二者在0.01水平(双侧)呈线性负相关关系(p<0.01)，说明光合植被覆盖度的增加对地表热环境具有降温效应，决定系数R2分别为0.63、0.55和0.65，由回归系数可知，2000、2008、2018年，光合植被覆盖度每增加10%，会使地表温度分别相应下降0.66、0.74、1.09 ℃。地表反照率与地表温度的回归分析结果表明二者在0.01水平(双侧)呈线性正相关关系(p<0.01)，说明地表反照率的增加对地表热环境具有增温效应，决定系数R2分别为0.35、0.40和0.48，由决定系数来看，相较于光合植被覆盖度，地表反照率和地表温度的相关性较小，由回归系数可知，2000、2008、2018年，地表反照率每增加10%，会使地表温度分别相应上升1.09、1.36、1.76℃。

4)地表下垫面覆被类型对地表反照率的分布和变化特征具有较大影响，矿业开采活动使得原有地表覆被遭到破坏，大面积裸露地表出现，改变了地表下垫面性质，在接收同等太阳辐射的条件下，出露的岩石和裸土热容量小，升温更快。并且由于采场出露的裸岩和排土场矿渣具有的高反照率对其周围区域的净辐射损失有重要影响，是造成采场和排土场地表温度迅速升高的主要原因。对于矿业开发密集区，今后应采取边开采边复垦的土地整治措施，以达到降低裸露地表面积，提高光合植被覆盖度，进而改善地表热环境高温聚集效应的目的。全球气候变化导致的气温升高也会造成一定程度的地表温度升高，但是如何剥离全球气温上升造成的影响是很难做到的。本文研究弥补了矿业开采导致矿区内部光合植被覆盖改变引起地表反照率时空分布变化，进而加剧地表热环境异质性分析的不足，有助于人们理解光合植被覆盖度和地表反照率之间潜在的关系。