融合主被动遥感的乌海矿区土地损伤测度

郭山川, 汤 傲, 李效顺,2①, 卞正富, 侯湖平, 倪 衡

(1.中国矿业大学环境与测绘学院, 江苏 徐州 221116; 2.中国矿业大学江苏省资源环境信息工程重点实验室, 江苏 徐州 221116)

煤炭资源开发利用对矿区土地造成持续性扰动而引发的植被退化、地表形态变化、土壤肥力衰退与水土流失等土地损伤在土地资源承载力脆弱的干旱半干旱矿区尤为严重。与土地损毁相比,土地损伤着重于扰动效应引起的土地功能性与结构性衰变,其土地综合整治方法包含土地复垦、生态修复和自然恢复等[1-3]。土地损伤机理主要有:开采沉陷过程导致的土壤竖向迁移运动,驱使沉陷盆地内土壤的空间分布特征重建;开采、运输、沉陷等过程中重金属、铵态氮、速效磷等物质发生转移从而改变土壤化学特征;土壤微生物多样性、数量和种群结构在采矿胁迫下也呈现出显著响应特征。由此可见,矿区土地在采矿持续扰动下遭受长期功能性与结构性损伤,导致其植被覆盖与形态结构发生趋势性变化[3-6]。

传统土地损伤监测手段不仅费时、费力,且监测点呈离散分布,难以进行矿区全覆盖监测,亦不能揭示矿区土地损伤空间分布规律。近年来,国内外学者利用3S技术和新兴的合成孔径雷达干涉测量技术(InSAR)对矿区土地利用类型、景观格局、植被覆盖和地表沉陷等指标进行了全方位监测[7-11]。部分学者基于上述指标分析了采矿扰动下矿区土地生态系统的时空演变规律,进而讨论了矿区生态演替驱动力、恢复力以及协同机制等[12-17]。然而对矿区土地损伤的测度研究还相对滞后,尤其是融合主被动遥感提取矿区土地损伤因子,建立损伤理论模型以定量分析西北矿区土地损伤的研究较为少见。同时矿区土地损伤的空间分布特征、时序演替规律以及土地修复范式的科学匹配等问题均需在土地损伤有效测度的基础上开展进一步研究。

鉴于此,以典型西北干旱半干旱地区乌海矿区为例,从采矿扰动过程分析识别关键土地损伤因子,利用主、被遥感提取土地覆被类型、地表形变以及植被覆盖度等土地生态参量,构建以植被损伤指数和地表形变指数为主控参量的土地损伤测度模型,实现矿区土地损伤的空间全覆盖监测,定量揭示其土地损伤空间分布规律,有助于西北矿区土地修复范式的精准化研究。

1 研究区概况与数据介绍

1.1 研究区概况

乌海市是新兴煤炭资源型城市,辖区内矿产资源分布广泛,因此研究区乌海矿区范围以乌海市行政边界为参考。乌海市位于内蒙古自治区西南部(39.15°～39.52° N,106.36°～107.05° E),下辖海勃湾区、乌达区和海南区3个市辖区(图1)。地势呈东西高、中部低,3山成南北走向平行排列,中间形成2条平坦的谷地。地貌为构造侵蚀中低山地、剥蚀丘陵区、山前堆积冲洪积扇区和黄河冲积堆积阶地。乌海市属于典型干旱半干旱地区,年均温9.6 ℃,年均降水量159.8 mm,年均蒸发量3 289 mm。区域内累计查明煤炭资源储量为28.8亿t,是国家重要焦煤基地,持续大规模的煤炭资源开发给国家和地区带来巨大的经济效益,同时也对区域生态环境带来了负面影响。

图1 乌海市示意图Fig.1 The location of Wuhai City

1.2 数据准备

采用数据主要包括遥感数据、合成孔径雷达(SAR)数据、数字高程模型和行政区划数据等(表1)。

表1主要数据来源

Table1Datasources

数据 采集时间数据来源备注 Landsat TM2000年8月22日美国地质勘探局(USGS)云量0.34%,空间分辨率30 m Landsat TM2005年6月1日USGS云量0.12%,空间分辨率30 m Landsat TM2009年6月28日USGS云量0,空间分辨率30 m Landsat OLI2015年9月1日USGS云量0.08%,空间分辨率30 m 行政区划图2012年地方机构比例尺1∶10万 地形图2012年地方机构比例尺1∶10万 矿区分布图2012年地方机构比例尺1∶10万 矿山环境治理规划评估2012年地方机构比例尺1∶10万 矿山地质环境影响评价图2012年地方机构比例尺1∶10万 Sentinel-1A雷达数据2015-01-13—2015-12-27(共9景)欧洲航天局(ESA)空间分辨率3.7×13.9 m SRTM数字高程模型2003年美国航空航天局(NASA)空间分辨率30 m

USGS、ESA、NASA数据来源网站分别为https:∥glovis.usgs.gov/、https:∥scihub.copernicus.eu/和http:∥srtm.csi.cgiar.org/。

由于季节及时间累积效应差异等原因易使遥感影像信息产生伪变化,因此应尽量采用同时相、同时间间隔的遥感影像作为数据源进行处理分析。通过辐射校正、投影变化、几何校正、数据融合和影像剪裁等方法对研究区2000年8月、2005年6月、2009年6月(2010年夏季遥感影像缺失)的Landsat TM和2015年9月的Landsat OLI遥感影像进行预处理。SAR影像为2015年1月13日至2015年12月27日共9景3.7×13.9 m分辨率、干涉宽幅升轨模式(IW)的Sentinel-1A数据。

2 测度模型及方法

采矿活动对矿区土地的扰动形式主要有挖损、压占、沉陷、地裂缝和滑坡等，这些扰动过程的持续挤压式作用导致土地功能与结构损伤,反映在地表植被覆盖退化、土地覆被改变和地表形变等方面。因此,通过3S技术提取植被覆盖与土地覆被信息,InSAR技术反演地表形变,并构建损伤测度模型对矿区土地损伤进行定量测度和评估[18-19]。

2.1 土地损伤测度模型

以地表形变指数和植被损伤指数为主控变量,同时引入土地覆被类型变量,从空间像元级尺度构建西北矿区土地损伤测度模型,空间像元p土地损伤DL,p的测度模型为

DL,p=g(IND,p,TLC,p)+f(INVD,p,TLC,p)。

(1)

式(1)中,g(IND，p，TLC，p)为由像元p地表形变指数(normalized deformation index,NDI，IND)和土地覆被类型(land cover type,LCT，TLC)决定的土地结构模块损伤函数;f(INVD，p，TLC，p)为由像元p植被损伤指数(normalized vegetation damage index,NVDI，INVD)和LCT决定的土地功能模块损伤函数,IND，p，INVD，p∈[0，1]。

矿区土地损伤在持续扰动过程中存在累积效应和恢复力作用,在两者的耦合作用下,损伤测度边际函数呈先增大后减小特征,利用该特性建立形变变量和植被损伤变量与土地损伤的微分方程组,对微分方程组积分并进行欧拉变换得:

(2)

(3)

式(2)～(3)中,α、β为与TLC，p相关的权重系数,表示不同土地覆被类型对矿区土地损伤的敏感程度。结合乌海矿区实地调研资料和TD/T 1010—2015《土地利用动态遥感监测规程》[20],将研究区分为林草地、耕地、水域、裸地、建设用地和采煤用地6种土地覆被类型,利用德尔菲法对不同土地覆被类型的土地损伤敏感程度进行排序分析(表2)。

表2不同土地覆被类型的土地损伤敏感性分析

Table2Thesensitivityoflandusetypestolanddamageindex

土地覆被类型 土地损伤结构损伤函数系数(α)功能损伤函数系数(β)林草地DB耕地CA水域FF裸地EC建设用地AD采煤用地BE

A～F的分值赋为5～0。

如图2所示,测度模型输入参量为表达土地结构损伤的形变参量,表达土地功能损伤的植被损伤参量,以及由土地覆被类型决定的系数α、β。矿区损伤测度模型由结构损伤测度和功能损伤测度函数构成,输出结果为矿区内像元p的土地损伤值。

图2 乌海矿区土地损伤测度框架

2.2 损伤参量提取方法

利用遥感方法计算矿区归一化差分植被指数(normalized difference vegetation index, NDVI,INDV)以获取土地损伤测度模型中的植被参量,NDVI是遥感影像中近红外波段的反射值和红光波段的反射值之差与上述两者之和的比值[21],其有效值域为[0,1],对地面覆盖为云、水、雪等负NDVI值区域进行零值处理。矿区下沉盆地土地植被覆盖变化能一定程度上反映其土地生产功能损伤情况[22-23],因此植被损伤指数NVDI与NDVI函数关系可表达为INVD=1-INDV。

矿区地表形变参量信息提取采用小基线集合成孔径雷达干涉测量(small baseline subset interferometric synthetic aperture radar, SBAS-InSAR)技术,该技术以全天候、全天时、高分辨率和连续空间覆盖的优势在地表形变测量中得到广泛关注[24-25]。SBAS-InSAR作为重要时序InSAR技术之一,能够有效减小一般InSAR易受失相干、大气延迟效应等因素引起的误差,在矿区连续、剧烈形变场反演中具有独一无二的优势[26-27]。对SBAS-InSAR提取的地表形变速率(ground deformation velocity, GDV,VGD)进行无量纲归一化处理,得到地表形变指数(IND):

(4)

另外,采用支持向量机算法(support vector machine, SVM)对乌海矿区遥感影像的光谱和纹理融合特征进行监督分类学习,提取乌海矿区土地覆被类型。

3 结果与分析

3.1 植被提取结果

乌海矿区2000、2005、2009和2015年4个时相的植被覆盖分布如图3所示。从空间分布上看,整体上低植被覆盖水平反映了乌海矿区突出的生态脆弱性,而且东南区域植被覆盖水平显著低于其他区域。

图3 2000—2015年乌海矿区NDVI的空间分布

为进一步分析乌海矿区植被覆盖信息提取结果的时序变化特征,计算乌海矿区范围内每个地理单元的平均NDVI值,结果见图4。可以发现,各期NDVI值均小于0.05,说明乌海矿区植被覆盖整体水平低;2000—2009年乌海矿区植被覆盖呈现逐期衰退趋势,期间植被覆盖水平下降22.9%,标准差呈递减趋势;而2009—2015年乌海矿区植被覆盖出现恢复反弹,反弹力度达10.8%,主要原因是地方政府和煤炭企业落实了生态文明建设政策方针[28-30],开展了系列土地整治、植树造林和环境治理等工程,同期标准差值激增印证了研究区植被覆盖恢复是土地综合整治工程造成的空间离散式恢复。

统计乌海矿区不同植被覆盖水平面积(表3),可以发现,2000—2005年,NDVI水平为>0.6～1的区域面积减小45.8 km2,>0.3～0.6的区域面积减小12.5 km2,>0.1～0.3的区域面积增加86.2 km2;2005—2009年,NDVI水平为>0.3～1的区域面积略有下降,共减小1.5 km2;2009—2015年,NDVI水平为>0.6～1的区域面积增加52.6 km2,>0.3～0.6的区域增加21.1 km2,生态治理效果明显。

3.2 沉陷反演结果

以时相为2015年9月22日和2015年11月9日的2景SAR图像差分干涉图为例(图5)展示SBAS-InSAR干涉结果,可以发现,干涉图纹清晰说明大气延迟、斑点噪声等误差相位通过空间低通和时间高通滤波得到有效削减,差分干涉结果质量较好;干涉图中呈现出系列沉陷漏斗,沉陷漏斗的位置、几何形状与矿区下沉盆地特征相似,说明区域在观测期间的地表沉陷主要由采矿活动引起。

图4 2000—2015年乌海矿区NDVI平均值变化

表32000—2015年乌海矿区不同NDVI级别的面积与比例

Table3TheareaandproportionofdifferentNDVIlevelsinthestudyareafrom2000to2015

NDVI值2000年2005年2009年2015年比例/%面积/km2比例/%面积/km2比例/%面积/km2比例/%面积/km2 0～0.166.31 163.364.71 135.447.0824.361.11 071.9 >0.1～0.324.3426.929.3513.147.1825.828.8504.4 >0.3～0.66.3110.55.698.05.799.46.9120.5 >0.6～13.053.30.47.40.34.63.357.1 总计1001 7541001 7541001 7541001 754

乌海矿区SBAS-InSAR地表形变监测结果如图6所示。从采煤沉陷地地表形变场空间分布来看,地表沉降场与矿区地理位置有强空间相关性,沉降场主要分布在煤炭资源开采区,如乌海东部和中部山区,以及黄河以西的乌达矿区,而其余地方一般没有沉陷;采煤工作面地理位置多分布于地形起伏多变、地质条件复杂的山区,导致开采下沉盆地等沉线表现形式为系列同心不规则椭圆群。

利用GIS统计乌海矿区采煤沉陷监测结果,发现乌海矿区开采沉陷地面积为230.03 km2,占乌海矿区总面积的13.78%;年沉降速率为>200～800 mm·a-1的区域面积为44.98 km2,占沉陷地面积的19.55%;年沉降速率为>50～200 mm·a-1的区域面积为101.33 km2,占沉陷地面积的44.05%;年沉降速率为10～50 mm·a-1的区域面积为83.72 km2,占沉陷地面积的36.40%。

3.3 土地损伤测度

由图6可见,乌海矿区土地损伤空间分布呈现出以高强度开采工作面为中心向四周辐射减弱的特征:土地损伤值为>4.5～8.0的高损伤区地理空间位置主要集中在东部山区、中部山区以及黄河以西等煤炭资源开采强度大的区域;土地损伤值为>1.0～4.5的中度损伤区围绕在高损伤区周边;土地损伤值为0～1.0的非损伤区由于扩散衰减显著,其地理分布与矿区的空间相关性较弱。研究区共有586.57 km2土地在采矿扰动下遭到破坏,其中土地损伤值为>4.5～8.0的高损伤区面积为42.36 km2,占总面积的2.68%;中损伤区面积为544.21 km2,占比34.43%;研究区土地非损伤区面积为994.02 km2,占比62.89%。

4 讨论与结论

笔者对采矿扰动形式、矿区土地响应特征进行分析,同时针对土地损伤过程中其结构与功能的变化机制,构建西北矿区土地损伤测度模型,融合主被动遥感提取模型关键参量,以乌海矿区为例获取了该区域的土地损伤值,分析了其空间分布规律。以乌海市矿山环境治理规划评估图、矿山地质环境影响评价图、地面实测数据以及各矿区提供的矿山地质环境治理数据为验证数据对比分析笔者提出的土地损伤测度结果,发现测度结果与实际资料在空间上具有较强相关性,而且在时序上的关联性也得到了很好的检验。

(1)采矿扰动下,乌海矿区在2000—2009年间植被覆盖逐渐萎缩,2009—2015年期间植被覆盖水平提高是由于地方进行了局部区域的土地整治与修复工程,其标准差激增说明植被恢复过程是空间不连续的离散式事件。

图5 乌海矿区局部差分干涉图

(2)采矿扰动下,乌海矿区土地结构破坏表现在区域内出现大范围地表沉陷场,其沉陷速率最大超800 mm·a-1,监测到的沉陷漏斗与煤炭资源开采区有较强的空间相关性,其下沉等值线形式特征为系列不规则椭圆群。

(3)通过西北矿区土地损伤模型测算,发现土地损伤在空间分布上出现以开采工作面为中心向周边扩散衰减的现象,说明采矿扰动是导致乌海矿区土地功能性和结构性损伤的主因之一,另一方面也反映出乌海矿区土地修复工作刻不容缓。

(4)乌海矿区土地损伤测度结果与实际调研资料相比,具有较一致的空间相关性和时序关联性,印证了该模型的有效性。土地损伤测度模型有助于西北矿区长时序、全覆盖土地损伤监测体系的搭建、土地损伤精准修复策略的研究以及损伤-修复动态反馈机制的建立等,有利于“青山”与“金山”协同发展。

图6 乌海矿区地表沉降速率、土地覆被和土地损伤值示意

笔者利用提出的土地损伤测度模型能够有效应用于乌海矿区土地损伤调查,对采煤塌陷地治理工作有一定的现实意义,但也存在不足之处:土地损伤提取结果缺乏大量实测数据进行定量验证;提供SAR影像的Sentinel-1A卫星于2014年发射运行,形变监测时相与植被监测时相无法向前吻合,导致乌海矿区土地损伤监测时相较为单一。后续的研究重点是针对上述2点问题展开进一步分析,以更有力地指导矿区土地整治与生态修复工作。