一种探测露头残煤自燃区域的方法与流程

本发明涉及矿山生态环境监测的技术领域，尤其是涉及一种探测露头残煤自燃区域的方法。

背景技术：

露头残煤主要是采掘揭露、煤层天然露头产生的，暴露在空气中，引起煤自燃现象。露头残煤自燃不但造成了大量煤炭资源的浪费，而且对生态环境造成了巨大的破坏，对矿区生活健康和生产安全造成了威胁。对煤自燃区域生态环境进行监测及修复，需将煤自燃区域提取出来。

遥感技术具有覆盖面广、成本低、获取便捷便利快捷的特点，因此被广泛应用于煤自燃方面的探测。如，《一种基于地表裂隙和温度的煤火勘探方法》，发明公开专利号为cn102176070b；《--一种煤田自燃的遥感监测方法》，发明公开专利号为cn104598874a；《一种基于卫星遥感的煤矿火灾监测方法》，发明公开专利号为cn102074093a；《热红外与雷达遥感联合探测煤火的方法》，发明公开专利号为cn108983309a。这些发明均是基于煤自燃的温度异常原理，适合大型的天然煤田火区。相对煤田火区，露头残煤自燃区规模小、分布零散，且伴随采掘扰动，地表有显著的地表覆盖变化特征，这给传统的温度异常探测带来了不利影响。对于露头残煤自燃区域的精准探测，目前仍缺乏有效方法。

技术实现要素：

为满足上述需求，本发明提供了一种探测露头残煤自燃区域的方法，为解决上述问题提供了一种可行性方案，为进行矿区生态环境修复，保证矿区的生活健康和生产安全提供了有力的技术支撑。

本发明采用的技术方案是一种探测露头残煤自燃区域的方法，包括以下步骤：

步骤1)获取研究区的煤田地质图，对煤田地质图进行矢量化，从矢量化的煤田地质图中选取为煤层的区域，该区域为研究区的煤层分布区，记为a；

步骤2)获取研究区的多期遥感影像，对影像进行预处理，得到预处理后影像；利用决策树方法对预处理后影像进行分类，得到分类后影像；再利用分类后变化检测方法提取露头煤和裸地增加的区域，记为b；其过程主要包括步骤2-1)至步骤2-3)；

步骤2-1)对遥感影像进行预处理，包括几何校正、辐射校正、图像融合及图像增强；

步骤2-2)利用决策树方法对预处理后影像进行分类；所述决策树方法的分类指标包括高程、坡度、归一化植被指数及光谱特征；利用决策树方法将预处理后影像分为林地、草地、耕地、建筑用地、裸地和露头煤六种地表覆盖类型；

步骤2-3)利用分类后变化检测方法提取露头煤和裸地增加的区域；所述分类后变化检测方法的具体步骤是：首先选取两个相邻时期的分类后影像，记为f和g，然后分别从f和g中提取出地表覆盖类型为露头煤和裸地的区域，记为f1和g1，其他区域为f2和g2；其次，叠置f1和g1，从g1中扣除g1与f1重合区域，得到两个相邻时期露头煤和裸地增加的区域，记为b；

步骤3)利用温度反演算法对预处理后遥感影像进行温度反演，得到地表温度影像；再利用自然间隔点分割法将地表温度影像划分为低、较低、中、较高、高五个级别；提取级别为高的区域作为高温异常区，记为c，其过程主要包括步骤3-1)至步骤3-4)；

步骤3-1)将预处理影像的像元dn值转转换为辐射亮度值，转换公式为：

式中lλ为辐射亮度值，单位为w·m^-2·sr^-1·μm^-1，dn为像元灰度值，lmax，lmin分别为相应波段灰度值分别为255和1时所对应的光谱辐射率，单位为w·m^-2·sr^-1·μm^-1；

步骤3-2)将辐射亮度值转换为亮度温度值，转换公式为：

trad＝k2/ln(k1/lλ+1)(公式2)

式中trad为亮度温度值，k1为常数，单位为w·m^-2·sr^-1·μm^-1；k²为常数，单位为k；

步骤3-3)求取地表比辐射率；自燃区域，即地表覆盖类型为裸地、露头煤类型的区域，地表比辐射率等于常数0.970；非自燃区域，即地表覆盖类型为林地、草地、耕地、建筑用地的区域，地表比辐射率的计算公式为：

式中εb为地表覆盖类型为建筑用地的地表比辐射率；εn为自然表面，即地表覆盖为林地、草地、耕地的地表比辐射率；pv为植被覆盖度，计算公式为：

pv＝[(ndvi-ndvimin)/(ndvimax-ndvimin)]²(公式4)

式中ndvi为归一化植被指数，ndvimax是ndvi的最大值，ndvimin是ndvi的最小值ndvi；当像元的ndvi大于ndvimax时，pv取值为1；当ndvi小于ndvimin时，pv取值为0；

步骤3-4)求取地表温度，计算公式为：

式中，ts为地表温度，λ为热红外波段的中心波长，c为光速，h为普朗克常数；σ为波尔兹曼常数；ε为地表比辐射率；

然后，利用自然间隔点分割法确定地表温度的自然间断点，利用自然间断点将地表温度分为低、较低、中、较高、高五个级别；提取地表温度级别为高的区域作为高温异常区，记为c；

步骤4)最后将煤层分布区a、露头煤和裸地增加区b、高温异常区c三者进行空间叠置分析，提取a、b、c三者重叠的区域作为露头残煤自燃区域。

本发明的技术效果在于：

(1)利用露头残煤自燃区域地质构造、地表覆盖变化、热异常特点，实现对露头残煤自燃区域的精准探测，排除了地面采掘扰动、露头残煤规模小、露头残煤分布零散等不利因素，保证了露头残煤自燃区域提取结果的可靠性。

(2)利用遥感影像实现露头残煤自燃区域的精准探测，不但成本较低，耗费时间较短，而且覆盖范围较广，可节约大量人力、物力。

(3)本发明最终提供露头残煤燃烧中心的地图，可以为煤自燃的环境监测、露头煤的开发和管理、生态精准修复提供直接支撑。

附图说明

图1为本发明的技术流程图；

图2为实施例的2002年和2006年的露头煤和裸地增加区图；

图3为实施例的2006年和2010年的露头煤和裸地增加区图；

图4为实施例的2002年份地表温度分级结果图；

图5为实施例的2006年份地表温度分级结果图；

图6为实施例的2010年份地表温度分级结果图；

图7为实施例的将2002-2006和2006-2010裸地和露头煤增加区域，2002年、2006年、2010年高温异常区叠加在一起的图。

具体实施方式

本发明中研究区为有露头残煤燃烧区域，遥感影像数据来源较广泛，多光谱及高光谱数据均可，例如：landsat、quickbird、sentinel、worldview、spot等数据。在本实施例中，以landsat影像为例，并结合图1进行详细描述，具体步骤如下：步骤1：选定车道沟、炭窑坪、石炭沟、寺儿沟四个采矿区作为研究区，四个采矿区的相对位置如图2和图3所示。获取研究区的煤田地质图，利用arcgis软件对煤田地质图进行矢量化，从矢量化的煤田地质图中选取为煤层的区域，并重采样至与遥感影像同一尺度，该区域为研究区的煤层分布区，记为a。

步骤2：获取研究区的多期遥感影像，对影像进行预处理，得到预处理后影像。利用决策树方法对预处理后影像进行分类，得到分类后影像。再利用分类后变化检测方法提取露头煤和裸地增加的区域，记为b，其过程综合使用envi与arcgis软件完成，主要包括步骤2-1至步骤2-3。

步骤2-1，利用envi软件对遥感影像进行预处理，包括几何校正、辐射校正、图像融合及图像增强。

步骤2-2，利用决策树方法对预处理后影像进行分类。所述决策树方法的分类指标包括高程、坡度、归一化植被指数及光谱特征。利用决策树方法将预处理后影像分为林地、草地、耕地、建筑用地、裸地和露头煤六种地表覆盖类型。决策树方法的分类过程是：设定高程阈值h1，归一化植被指数阈值ndvi1、ndvi2、ndvi3，坡度阈值s1，第五波段光谱阈值t1。对所有遥感影像的像元进行判断：当坡度小于h1、且归一化植被指数小于ndvi1时，该像元为耕地；当坡度小于h1、且归一化植被指数大于ndvi1时，该像元为建筑用地；当坡度大于h1、归一化植被指数小于t2、且第五波段光谱值小于t1时，该像元为露头煤；当坡度大于h1、归一化植被指数小于t2、且第五波段光谱值大于t1时，该像元为裸地；当坡度大于h1、归一化植被指数大于t2、且归一化植被指数大于t3时，该像元为林地；当坡度大于h1、归一化植被指数大于t2、归一化植被指数小于t3时、且坡度小于s1时，该像元为耕地；当坡度大于h1、归一化植被指数大于t2、归一化植被指数小于t3时、且坡度大于s1时，该像元为草地。

步骤2-3，利用分类后变化检测方法提取露头煤和裸地增加的区域。所述分类后变化检测方法的具体步骤是：首先选取两个相邻时期的分类后影像，记为f和g。然后利用arcgis软件分别从f和g中提取出地表覆盖类型为露头煤和裸地的区域，记为f1和g1，其他区域为f2和g2。其次，叠置f1和g1，从g1中扣除g1与f1重合区域，得到两个相邻时期露头煤和裸地增加的区域，记为b。本实施例中考虑了两组相邻时期，分为为2002和2006年、2006和2010年。图2和图3中分别给出了2002年、2006年以及2006年、2010年的矿区内露头煤和裸地区域。2002年与2006年的露头煤和裸地叠加(如图2所示)、2006年与2010年的露头煤和裸地叠加(如图3所示)，叠加区域就是相邻时期露头煤和裸地增加的区域。

步骤3：利用温度反演算法对预处理后遥感影像进行温度反演，得到地表温度影像。再利用自然间隔点分割法将地表温度影像划分为低、较低、中、较高、高五个级别。提取级别为高的地表温度影像作为高温异常区，记为c，其过程主要包括步骤3-1至步骤3-4。

步骤3-1，利用经过几何校正的影像像元的dn值计算辐射亮度值。计算公式为：

式中lλ为辐射亮度值，单位为w·m^-2·sr^-1·μm^-1，dn为像元灰度值，lmax，lmin分别为相应波段灰度值分别为255和1时所对应的光谱辐射率，单位为w·m^-2·sr^-1·μm^-1。针对landsat影像，对于tm第6波段：lmax＝15.303w·m^-2·sr^-1·μm^-1，lmin＝1.238w·m^-2·sr^-1·μm^-1；对于etm+第6波段低增益波段lmax＝17.040w·m^-2·sr^-1·μm^-1，lmin＝0，高增益波段lmax＝12.650w·m^-2·sr^-1·μm^-1，lmin＝3.200w·m^-2·sr^-1·μm^-1。将其带入公式1中，计算tm/etm+第6波段辐射亮度值的公式为：

对于tm6，公式为：

lλ＝0.05537dn+1.1826(公式6)

对于etm+6，公式为：

lλl＝0.067dn-0.067

lλh＝0.0372dn+3.16(公式7)

式中lλh，lλl分别为低增益和高增益波段的辐射亮度。

步骤3-2，利用辐射亮度值计算亮度温度值，计算公式如下：

trad＝k2/ln(k1/lλ+1)(公式2)

式中lλ为上面求取的辐射亮度值，trad为量度温度值，k1为常数，单位为w·m^-2·sr^-1·μm^-1。k2为常数，单位为k。对于tm数据，k1＝607.76w·m^-2·sr^-1·μm^-1，k2＝1260.56k；对于etm+数据k1＝666.09w·m^-2·sr^-1·μm^-1，k2＝1282.71k。

步骤3-3，计算地表比辐射率。自燃区域，即地表覆盖类型为裸地、露头煤类型的区域，地表比辐射率等于常数0.970。非自燃区域，即地表覆盖类型为林地、草地、耕地、建筑用地的区域，地表比辐射率的计算公式为：

式中εb为地表覆盖类型为建筑用地的地表比辐射率；εn为自然表面，即地表覆盖为林地、草地、耕地的地表比辐射率；pv为植被覆盖度，计算公式为：

pv＝[(ndvi-ndvimin)/(ndvimax-ndvimin)]²(公式4)

式中ndvi为归一化植被指数，ndvimax是ndvi的最大值，ndvimin是ndvi的最小值ndvi。当像元的ndvi大于ndvimax时，pv取值为1；当ndvi小于ndvimin时，pv取值为0。

步骤3-4，求取地表温度，计算公式为：

式中，ts为地表温度，λ为热红外波段的中心波长，c为光速，h为普朗克常数；σ为波尔兹曼常数；ε为地表比辐射率。

然后，利用自然间隔点分割法确定地表温度的自然间断点，利用自然间断点将地表温度分为低、较低、中、较高、高五个级别。提取地表温度级别为高的区域作为高温异常区，记为c。图4、图5和图6为实施例的研究区域2002、2006、2010三个年份地表温度分级结果，图中1、2、3、4、5分别代表低、较低、中、较高、高五个级别。

步骤4：利用arcgis软件对煤层分布区a、露头煤和裸地增加区b、高温异常区c三者进行空间叠置分析，提取a、b、c三者重叠区域作为露头残煤自燃区域。图7中将2002-2006裸地和露头煤增加区域，2006-2010裸地和露头煤增加区域、2002年高温异常区、2006年高温异常区、2010年高温异常区叠加在一起，重叠区域即是露头残煤自燃区域，在露头残煤自燃区域标记了燃烧中心，如图7所示。该图件可以为露头残煤的开发、环境影响评价、土地复垦提供支撑。