基于Landsat 8遥感影像的峁底矿区土壤含水量反演研究

2023-02-22王丹萍康建荣练靖文崔腾飞

有色金属（矿山部分） 2023年1期

王丹萍，康建荣，练靖文，崔腾飞

(江苏师范大学地理测绘与城乡规划学院，江苏徐州 221116)

煤炭开采引起了一系列的环境问题，其中包括对土壤含水量的影响。煤炭资源的开采往往会造成地面沉陷，从而导致地下含水层遭到破坏，使得地下水位不断下降，因此造成土壤水分的减少[1]。由于采煤沉陷造成土壤含水量的降低，若使土壤含水量尽可能恢复，则需要一个较长的周期[2]。对于矿区植被的生长及受损后的恢复，土壤含水量是关键的环境因子[3]。同时，土壤含水量对农业生产有着不可替代的作用，其直接影响作物产量。所以土壤含水量的监测及预估对于研究区内生态环境、农业等领域的分析和治理具有很大的作用。

峁底煤矿自开采以来，周边地形发生了较大的变化，改变了土壤中水分含量及分布状况，影响到植被的生长，使得许多土地类型演变为裸土。由于研究区地形主要为山地，在降水较多时，因为缺少植被的保护，采空沉陷形成的坡面容易引起水土流失，甚至引发泥石流等自然灾害的发生。该矿区处于半干旱地区，地表水资源极其宝贵。土壤水分是制约半干旱地区植被生长和生态环境恢复的关键性因子[4]。

针对以上问题，为探索研究区内土壤含水量时空分布变化特征，本文采用植被供水指数法对研究区内的土壤含水量进行反演。以峁底煤矿所在乡镇为研究区，以长时间序列的Landsat 8遥感影像为数据源，反演研究区2013—2021年的土壤含水量，并对其整体的时空分布特征进行分析，进而揭示研究区土壤含水量的变化情况，最终结合煤矿开采工作对土壤含水量的影响，为矿区内水资源问题研究提供一定理论和数据帮助。

1 研究区与数据

1.1 研究区概况

奥家湾乡，隶属于山西省吕梁市兴县，地处兴县东部，位于晋西北丘陵地区，地势整体东高西低，境内群山蜿蜒，纵横交错，丘陵密布，地形支离破碎。该地气候为温带大陆性气候，四季分明，年平均气温较低，一般夏季多雨，冬季少雪，但近几年春夏季气候偏旱；境内河道属黄河流域，主要河道有蔚汾河、岚漪河；境内有丰富的矿产资源、牧草资源等。

峁底煤矿，隶属于奥家湾乡，位于38°25′54″～38°27′39″N，111°09′47″～111°11′23″E。矿区地处芦芽山东麓，此处地形较复杂，侵蚀冲刷严重，形成了接近南北走向的山脊及沟谷，整体地势南高北低(见图1)。井田中部为最高处，海拔为1 234.9 m；井田北部为最低处，海拔为1 028.0 m，最大相对高差为206.9 m，属于中低山区。矿区内的河流属于黄河流域蔚汾河水系，矿区内无常年性河流，各沟谷一般干涸无水，到雨季时作为泄洪通道。该地区一月份最低平均气温为-7.8～11.2 ℃，极端最低为-29.3 ℃；七月份最高平均气温为22.3～25 ℃，极端最高为38.4 ℃。年平均降水量为231.4～688.9 mm，60%以上的降水主要集中在夏季及秋初；年平均降水蒸发量为2 090.8 mm，达到降水量的4倍。矿区及周边地区所覆盖的植被主要以灌木丛为主，分布较稀疏，同时有少量农田分布。由于地处黄土高原地区，气候干旱，经气候变化及煤矿开采的影响，境内山地地表类型逐渐演变为裸土，植被日益减少。

图1 峁底煤矿西南部及周边环境的八月份无人机影像数据Fig.1 UAV image data of the southwest and surrounding environment of Maodi coal mine in August

1.2 研究数据

随着遥感领域研究的不断发展，目前已经实现了大范围区域内土壤含水量的实时监测和动态预测，通过遥感影像监测土壤含水量能够较好地反映出土壤含水量在一定时间和空间上土壤水分的变化特点[5]。Landsat 8卫星携带的陆地成像仪(OLI)和热红外传感器(TIRS)的波段划分精细，对植被、温度以及土壤含水量都有极佳的监测效果[6]。Landsat 8影像数据量大且种类丰富，可以完全覆盖监测区域，因此，本文所采用的影像资料是来自于地理空间数据云(https：//www.gscloud.cn/)2013—2021年的Landsat 8 OLI_TIRS影像数据，云量均在10%以下，轨道号为126/33。

地表温度反演精度评定数据是来源于地理空间数据云MODIS陆地标准产品中的2015年12月17日和2016年4月7日的MOD11A1 1 KM地表温度/发射率产品。

1.3 数据预处理

在ENVI5.3中对Landsat 8遥感影像数据进行预处理，包括裁剪、辐射定标、大气校正等操作，以改善遥感影像的质量，为后续实验作好准备。同样在ENVI5.3中对MODIS地表温度产品数据进行预处理，包括格式转换、重采样、重投影、裁剪、温度换算等，为后续精度评定做准备。

2 研究方法

2.1 研究方法及原理

植被供水指数(VSWI)是一种综合植被状态和地表温度信息，反演有植被覆盖地区表层土壤含水量的监测方法。在植物供水状态正常的条件下，一定生长期内的植被指数和植被冠层温度都能维持在一定范围内，但在植被处于缺水状态时，植被的生长情况会遭受影响，植被指数会随之降低，而植被冠层温度则会随之上升[7]。VSWI值可表示植被相对干旱程度，VSWI值越小，植被的干旱程度越严重，土壤含水量越低，反之土壤含水量越高。公式为：

(1)

式中，NDVI为归一化植被指数；LST为地表温度。

2.2 地表参量确定

2.2.1 归一化植被指数(NDVI)

归一化植被指数一般与植被的繁盛程度有关，因此，可以利用归一化植被指数对植被覆盖度等进行研究分析。一般情况下，当作物缺水时，其生长情况会受到影响，植被指数也会随之降低。公式为：

NDVI=(ρNIR-ρRED)/(ρNIR+ρRED)

(2)

式中，ρNIR为近红外波段反射率；ρRED为红光波段反射率。本研究使用Landsat 8 OLI传感器的第四波段红光波段和第五波段近红外波段来计算NDVI值。

2.2.2 地表温度(LST)

地表温度是一个反映地表能量平衡的客观指标。对于裸地，地表温度被称为土壤温度；对于植被，地表温度被称为冠层温度[8]。由于土壤温度的变化能够间接反映土壤水分的变化，当冠层温度升高时，植被的缺水程度也逐渐升高[9]。本文中采用了大气校正法、单窗算法、基于影像的反演算法(Image-based Method，简称IB算法)来反演地表温度，并将三种方法的反演结果与MODIS地表温度产品进行精度验证，选出相关性最高的反演结果进行VSWI值的提取。

1)大气校正法。利用大气校正法反演地表温度的计算公式为：

(3)

式中，B(Ts)为黑体辐射亮度；K1、K2为指定波段的热转换参数，分别为774.885 3 W·(m2·sr·μm)-1和1 321.078 9 K。

黑体辐射亮度公式为：

(4)

式中，Lλ为影像热红外波段进行辐射校正后得到的星上辐射亮度；ε为地表比辐射率；τ、Lup、Ldown分别为大气透过率、大气向上辐射亮度、大气向下辐射亮度，三者可从NASA官网(https：//atmcorr.gsfc.nasa.gov/)获取。

地表比辐射率的计算可通过基于NDVI阈值的方法[10]进行估算，公式为：

εwater=0.995 (NDVI<0)

(5)

(0≤NDVI≤0.7)

(6)

(NDVI>0.7)

(7)

式中，PV为植被覆盖度，其计算公式为：

(8)

式中，NDVIsoil和NDVIveg分别为裸地和植被完全覆盖的像元NDVI值，本研究中选取NDVIsoil=0.05，NDVIveg=0.7。

2)单窗算法。本文中的单窗算法使用的是WANG和QIN等改进及验证后的单窗算法[11]，其计算公式为：

LST=

(9)

C=ετ

(10)

D=(1-τ)[1+(1-ε)τ]

(11)

式中，ε为地表比辐射率；T6为地面亮温值；a、b为常数，分别为-67.355 351和0.458 606；Ta为大气平均作用温度；τ为大气透过率，可通过NASA官网获取。

大气平均作用温度公式为：

Ta=16.0110+0.92621T0(中纬度夏季)

(12)

Tb=19.2704+0.91118T0(中纬度冬季)

(13)

T0=t+273.15

(14)

式中，t为大气平均气温，本文通过天气网(http：//lishi.tianqi.com/)查询到当天最高和最低气温，根据遥感影像拍摄时间进行估计。

3)IB算法。利用IB算法反演地表温度的计算公式为：

(15)

(16)

式中，T为辐射亮温，单位为K；λ为热红外波段的中心波长值，对于Landsat 8影像，其值为10.9 μm；δ为Boltzmann常数，值为1.38×10-23J/K；h为Plank常数，值为6.626×10-34Js；c代表光速，值为2.998×108ms；ε为地表比辐射率。

3 结果与分析

3.1 实验分析

为了获取最佳的地表温度数据，将三种方法反演出的地表温度结果分别与MODIS温度产品进行相关性分析(图2)。大气校正法、单窗算法、IB算法与MODIS温度产品分析结果从整体上来看，都具有良好的线性关系，相关系数分别为0.677 7、0.677 6、0.677 6。

图2 三种不同算法反演的地表温度与MODIS温度产品线性拟合结果Fig.2 Linear fitting results of land surface temperature and MODIS temperature products by inversion of three different algorithms

根据相关性比较，本研究采用了大气校正法反演地表温度，由于无法精确地收集到研究区的大气平均气温，所以反演得到的地表温度精度会略低，但与MODIS温度产品数据仍具有较强的相关性。

通过采用植被供水指数法反演得出研究区土壤含水量的时空分布情况。从分布情况来看，颜色由浅到深，土壤含水量则由低到高，且土壤含水量程度主要以较旱为主。研究区春夏季气候较干旱，所以大部分地区土壤含水量较低，得到的结果与实际调查结果基本一致，此方法对于反演研究区土壤含水量具有一定的可行性。

3.2 土壤含水量时空分布特征

3.2.1 空间分布特征

根据研究区土壤含水量分布情况及相关资料，本文将土壤含水量等级分为5级：干旱(<20%)、轻旱(20%～40%)、适宜(40%～60%)、湿润(60%～80%)、潮湿(>80%)。从2013—2021年研究区土壤含水量的分布情况(图3)来看，土壤含水量等级主要以轻旱和干旱为主。研究区北部植被覆盖度高，主要以潮湿和湿润为主；中部主要以适宜和轻旱为主，境内有省道贯通，所以省道所处位置表现为干旱；南部植被覆盖度低，有大量裸土分布，主要以轻旱和干旱为主。

图3 2013—2021年研究区土壤含水量空间分布情况图Fig.3 Spatial distribution of soil moisture content in the study area from 2013 to 2021

2013年，土壤含水量整体偏高，总体上呈“北部和中部高，南部较低”的特点，北部土壤水分含量最高，以潮湿为主，而省道及附近村庄建筑土壤含水量则最低。2014年，土壤含水量呈现出“北部和东南部高，西部低”的特点，北部土壤水分含量最高，以潮湿为主，省道、附近乡村建筑及西部地区土壤水分含量较少。2015年、2017年和2020年土壤含水量分布特点均与2014年基本一致，土壤含水量总体上与2014年大致相同。2016年和2019年土壤含水量的分布特点与2013年比较接近，2016年土壤含水量总体上与2013年大致相同，但2019年整体土壤含水量低于二者。2018年，土壤含水量呈“北部和东部高，西部低”的分布特点，东北部土壤水分含量最高，以适宜为主，而省道、附近乡村建筑及西部地区则相对较少。2021年，土壤含水量总体偏低，呈现出“北高南低”的特点，北部土壤含水量偏高，以轻旱为主，省道、附近村庄建筑及中、南部大部分地区土壤水分含量都比较低。

3.2.2 时间分布特征

2013—2021年土壤含水量情况表现出整体波动降低的态势。总体上，2013—2015年，土壤水分含量呈递减趋势；2015—2018年，土壤水分含量也呈现出明显的降低趋势；而2018—2021年，土壤含水量出现了先增加后降低的变化。2013—2021年土壤含水量各类型(表1)所占面积及比例变化明显。

表1 2013—2021年研究区土壤含水量时间变化情况Table 1 Temporal change of soil moisture content in the study area from 2013 to 2021

2014年与2013年相比，潮湿、湿润、适宜类型所占比重明显下降，适宜类型比重下降更为显著，面积从68.272 2 km2变为17.918 1 km2，而较旱、干旱类型所占比重则显著上升，面积分别从13.989 6、3.416 4 km2变为56.2032、39.088 8 km2。2016年与2015年相比，潮湿、湿润、适宜类型所占比重明显上升，适宜类型比重上升更为显著，面积从16.213 5 km2变为57.483 9 km2，较旱、干旱类型所占比重则明显下降，干旱类型比重下降更为显著，面积从54.349 2 km2变为2.184 3 km2。2018年与2017年相比，潮湿、湿润、适宜、干旱类型所占比重明显下降，湿润类型比重下降更为显著，面积从12.094 2 km2变为1.307 7 km2，而较旱类型所占比重则明显上升，面积从55.027 8 km2变为86.551 2 km2。2020年与2019年相比，湿润、适宜类型所占比重明显下降，适宜类型比重下降更为显著，面积从59.208 3 km2变为17.027 1 km2，而较旱、干旱类型所占比重则明显上升，干旱类型比重上升更为显著，面积从5.507 1 km2变为40.470 3 km2。2021年与2020年相比，潮湿、湿润、适宜、较旱类型所占比重明显下降，适宜类型比重下降更为显著，面积从17.027 1 km2变为3.554 1 km2，而干旱类型所占比重则明显上升，面积从40.470 3 km2变为72.568 8 km2。

从整体上来看，2014年与2015年、2017年、2020年相比，土壤含水量各类型面积差别较小，土壤含水量较低，集中分布在20%～40%，面积分别占研究区总面积的43.68%、36.96%、42.77%、43.40%；2013年与2016年相比，土壤含水量各类型面积差别也较小，而土壤含水量则相对较高，集中分布在40%～60%，面积分别占研究区总面积的53.06%、44.68%。总体上，2013年、2016年及2019年的土壤含水量分布状况要优于其他年份，从土壤含水量程度占比(图4)中可明显看出，这三个年份的土壤湿度相对较高。

图4 2013—2021年土壤含水量程度占比柱状图Fig.4 Bar chart of the percentage of soil moisture content from 2013 to 2021

3.2.3 土壤含水量动态分析

研究区各类型土壤水分的面积变化有明显的差异(表2)，其中变化最大的是干旱类型，土壤含水量<20%，减少了69.1524 km2，变化率为53.75%；面积变化最小的是潮湿类型，土壤含水量>80%，增加了11.3724 km2，变化率为8.84%。

表2 2013—2021年研究区土壤含水量面积变化差异统计Table 2 Difference statistics of soil moisture content area change in the study area from 2013 to 2021

将土壤含水量面积变化率划分为7个等级：剧烈减少区(变化率<-30%)，减少区(-30%≤变化率<-10%)，轻度减少区(-10%≤变化率<-5%)，稳定区(-5%≤变化率<5%)，轻度增加区(5%≤变化率<10%)，增加区(10%≤变化率<30%)和剧烈增加区(变化率≥30%)[5]。从土壤含水量类型面积变化情况来看，干旱类型面积变化最大，可归为剧烈增加区；适宜类型面积变化程度仅次于干旱类型，可归为剧烈减少区。总体来看，2013—2021年研究区土壤含水量呈下降趋势。

3.3 煤矿开采对土壤含水量的影响

2013—2021年，峁底煤矿均处于开采状态。采矿区与省际干线公路相邻，位于公路南侧，从图3中可以看出采矿区的土壤含水量属于较高的区域，由于矿井会进行一系列的排水工作使矿区的开采能正常进行，所以此处的土壤含水量会较其他地区略高，而变化趋势与整体变化趋势基本一致。矿区与其他区域的土壤含水量在研究年份前几年差异比较显著，但后面几年差异逐渐缩小，考虑到煤矿开采会导致土壤受损，且土壤的自主修复周期较长，会产生许多无法闭合的裂缝，而这些裂缝会加速土壤水分的流失和蒸发[4]，因此矿区土壤含水量的变化趋势是整体降低的。针对矿区土壤含水量减少的问题，应及时对受损土壤进行人工修复，及时填补裂缝，减少土壤水分的流失。

4 讨论

采用植被供水指数法对土壤含水量进行反演，研究结果与实际调查结果基本相同，但由于缺乏实测数据，因此VSWI反演的土壤含水量结果与实际土壤含水量相关性较弱，且由于影像拍摄时间对应的当天降雨量均为0 mm，所以并未将降雨量作为影响土壤含水量的因素进行分析，后续研究会将数据时空范围继续扩大，并将降雨量、大气温度、地形等因素考虑进来，便于更精确地监测土壤水分的变化。同时，通过长时间序列的反演对比，可清晰地得出该地土壤含水量时空分布特征，对改善当地的生态环境和提高农业产量具有重大的意义。