第一作者: 李 峰(1979-),男,博士后,副教授/工程师,主要从事测量与遥感方面的研究。Email:lif1223@aliyun.com。
为了评估内蒙古乌达煤田煤火治理的效果,采用2008年治理前、2011和2013年治理中以及2015年治理后期的Landsat 5/8热红外波段影像,利用单窗算法反演4个年度的乌达煤田地表温度; 采用自适应梯度阈值(self-adaptive gradient-based thresholding,SAGBT)法圈定该煤田煤火区,分析煤火空间分布格局演变动态,并对探测结果进行了实地验证。结果表明,所识别煤火区的准确率为75%; 煤火区面积由2008年的1.194 km2演变到2015年的0.873 km2,呈显著下降趋势。总的来说,煤火治理取得了初步成效,实现了煤火区面积26.88%的减量; 但仍有73.12%的存量,说明今后仍须加大煤火治理力度。
In order to assess the fire-fighting effects in the Wuda coal field of Inner Mongolia, the authors adopted two Landsat 5 images(
煤火是指一定规模的煤体在自然条件下与空气接触后, 发生氧化、自燃乃至剧烈燃烧而引发地下煤裂隙、煤层露头和煤矸石堆火情的现象[1]。煤火的诱发因子多、可控性差。煤火不仅会烧毁大量的煤炭资源、危及矿井生产安全, 而且释放的有毒气体和元素会污染环境, 直接危害居民健康, 所排放的温室气体也会加剧全球变暖; 同时, 煤火灾害还会引发地面沉降、塌陷、地裂缝等地质灾害, 威胁当地基础设施和居民生命财产安全。对煤火的探测与防治对促进当地社会经济发展和环境保护具有重要的现实意义。
目前煤火探测的方法主要基于煤火区(以下简称“ 火区” )的地磁场特征、电场特征、地球化学场特征和热场特征[2]进行探测。其中, 基于地磁场特征、电场特征、地球化学场特征的火区探测方法的成本高, 效率低, 不适于圈定大面积火区范围; 利用卫星遥感数据反演地表温度特征、建立火区温度场的方法较为成熟, 已在煤田煤火监测中发挥了重要作用。蒋卫国等[3]根据1989— 2005年间获取的ETM+/TM 6数据, 应用单通道算法分析了乌达煤田煤火面积变化情况。邱程锦等[4]使用辐射传导方程反演了2001— 2010年间获取的4景TM 6影像的地表温度, 从宏观上分析了乌达煤火的分布状况。Huo等[5]采用Weng法用TM 6影像反演了1999— 2006年间的地表温度, 用于预测煤火走势。除TM6的热红外波段外, ASTER数据也被用于煤火识别中[6]。虽然ASTER数据的分辨率和热敏度都较高, 但因部分区域存档数据少、温度反演算法较复杂等原因, 妨碍了对其进一步的研究与应用。随着对高分辨率热红外数据需求的不断增长, 航空热红外遥感技术也常被应用于煤火监测。孙玉林等[7]根据高分辨率的航空三通道(可见光、中红外和热红外)影像, 圈定了2005年乌达煤田18处共350万m2的火区范围。Wang等[8]引入手持式热红外成像仪和无人机航空影像, 识别出大同地区的煤火。这类数据的分辨率和探测精度虽高, 但操作复杂、成本很高, 难以推广应用。与ETM+/TM 6相比, Landsat 8的TIRS传感器采用了更先进的红外光子检测技术, 使传感器的热敏度更高。徐涵秋[9]认为, 与Landsat 7热红外数据相比, Landsat 8数据在地表温度反演中更具优势。Yu等[10]认为大气传导方程法的精度最高, 劈窗算法精度次之, 单通道算法精度最低。姬洪亮等[11]认为单窗算法和单通道算法精度接近, Weng算法精度较低。Wang等[12]的研究结果则表明单窗算法优于单通道算法。地表温度反演算法是煤火识别的基础, 而将煤火高温异常值从背景值中分离出来则是煤火识别的关键。当前煤火高温异常提取方法有地表温度均值加其1~2倍标准差法、密度分割算法、移动窗口直方图法和自适应梯度阈值(self-adaptive gradient-based thresholding, SAGBT)法等。采用移动窗口直方图法提取煤火热异常的精度可达到80%; 采用SAGBT法提取煤火热异常, 在ASTER传感器过境时同步进行了地面温度采样与验证, 精度则达到了84.6%[13, 14]。
综上所述, Landsat 8 TIR数据比ETM+/TM 6数据精度有所提高, 用单窗算法反演地表温度的精度较高, SAGBT法在煤火热异常提取方面更加可靠。乌达煤田煤火的起火时间已经很久, 但多数研究限于资金和时间等原因, 只能集中在煤火燃烧最猛烈的时期, 未能进行可持续监测研究。鉴于2010年乌达煤田启动了为期4 a的全面灭火工程的情况, 本文将基于2008年煤火治理前、2011和2013年煤火治理中、2015年煤火治理后期的Landsat 5/8热红外波段影像, 分别采用单窗算法和SAGBT法提取乌达煤火高温燃烧中心作为火区范围, 并量化分析3个时段煤火的变化规律, 以期对煤火治理效果给予初步评价。
乌达煤田位于内蒙古自治区乌海市境内, 地处贺兰山北端、乌兰布和沙漠的南部边缘, 东距黄河约10 km。该煤田为一不对称向斜盆地, 东部被乌达逆断层所切割, 煤系地层埋藏较浅, 为近距离缓倾斜煤层群, 含煤地层属于石炭二叠系。整个煤田主要由苏海图、黄白茨和五虎山3个矿区组成。火区遍布煤田内, 多分布在露头废弃的小煤窑处, 并由露头向深部蔓延。
本文所采用数据为Landsat 5 TM(获取时间为2008年11月16日和2011年9月22日)和Landsat 8 TIRS(获取时间为2013年10月29日和2015年3月25日)共4景影像的热红外波段, 其中TM热红外波段(B6)分辨率为120 m, TIRS热红外波段(B10)分辨率为100 m。Landsat遥感影像数据来源于美国地质调查局网站(http: //glovis.usgs.gov/), 已被重采样为30 m分辨率, 主要用于反演煤田地表温度。MODIS的MOD021KM数据来源于美国国家航空航天局(https: //ladsweb.nascom.nasa.gov/data/search.html), 用于计算大气水汽含量。反演地表温度所用到的当日气温数据来源于中国气象科学数据共享服务网(http: //cdc.nmic.cn/home.do)。
2.1.1 地表比辐射率计算
为了有效排除煤田周边沙漠、居民地等地物对地表温度的干扰, 使用含煤地层的煤田边界选取研究区域。乌达煤田的原始地表破坏严重, 2008―2015年间的煤田地表主要以裸岩、裸土、煤矸石和煤炭为主, 可直接将煤田地表类型归为自然地表类型。地表比辐射率主要受地表组分、地表粗糙度和波长的影响, 因TM6波段的光谱范围基本覆盖了Landsat 8 B10的光谱范围, 故乌达煤田地表比辐射率ε 的计算仍可沿用TM6的比辐射率计算方法[15, 16], 即
ε =0.962 5+0.061 4 fv-0.046 1
fv=(NDVI-NDVIs)/(NDVIv-NDVIs) , (2)
NDVI=(NIR-R)/(NIR+R) , (3)
式中: fv为植被覆盖度; NIR和R分别为近红外和红光波段的反射率; NDVI为归一化差值植被指数; NDVIs为土壤的植被指数; NDVIv为植被的植被指数。
2.1.2 大气水汽含量计算
利用MODIS L1B数据的B2和B19的比值计算大气水汽含量, 可部分消除地表反射率随波长变化对大气透过率的影响, 提高大气水汽含量的反演精度。用大气水汽含量w(g/cm2)进一步估算大气透射率, 其计算公式[17] 为
w={[α -ln(ρ 19/ρ 2)]/β }2 , (4)
式中: ρ 2和ρ 19分别为B2和B19的反射率; 常数α 和β 分别为0.02和0.651。
2.1.3 星上亮温计算
根据定标系数对Landsat 5/8影像的DN值进行辐射定标, 转换成星上辐亮度值。利用普朗克反函数, 求解星上亮温Tt, 即
Tt=K2/ln(1+K1/Lt) , (5)
式中: Tt为星上亮温, K; Lt为卫星热红外波段星上辐射亮度值, W· m-2· sr-1· μ m-1; 对于Landsat 5 B6, K1=607.76 W· m-2· sr-1· μ m-1, K2=1 260.56 K; 对于Landsat 8 TIRS B10, K1=774.89 W· m-2· sr-1· μ m-1, K2=1 321.08 K[18]。
2.1.4 单窗算法
基于Landsat影像的地表温度反演算法主要有普适性单通道算法、单窗算法、大气辐射传导方程法和Weng方法等。单窗算法对Landsat 5/8热红外数据均具有较好的稳定性, 故本文主要采用单窗算法反演地表温度, 其计算公式为
Ts={a(1-C-D)+[b(1-C-D)+C+D]Tt-DTa}/C , (6)
式中: Ts为反演的地表温度, K; 对于Landsat 5卫星, 回归系数a和b的值分别为-67.355 35和0.458 608[19]; 对于Landsat 8卫星, 回归系数a和b在0° ~30℃时分别为-59.139 1和0.421 3, 在0° ~40℃时分别为-60.919 6和0.427 6; C和D为中间变量, C=ε τ , D=(1-τ )[1+(1-ε )τ ], ε 和τ 分别为地表比辐射率和大气透射率。
Landsat 5与Landsat 8的大气透射率估算方程分别见表1和表2[20]。
根据式(4)计算出2008年11月16日、2011年9月22日、2013年10月29日和2015年3月25日 的水分含量分别为0.411 g/cm2, 1.291 g/cm2, 0.426 g/cm2和0.555 g/cm2。选用的大气透射率方程分别为表1中的式③、式①及表2中的式⑤。
地下煤火的燃烧会在地表一定范围内引起高温异常, 在反演得到煤田地表温度的基础上, 通过对地表温度取阈值可从低温背景温度中识别出相对高温的火区。煤田内的环境、煤的特性、地理位置与地形的差异决定了只根据温度均值加标准差的阈值很难精确识别火区。为了更加准确地提取火区, 本文选用Du等[21]基于ASTER数据提出的SAGBT法识别火区, 识别精度可高达84.6%。SAGBT算法基于火区边界过渡到低温背景区时温度出现快速衰减的假设, 首先, 对反演的地表温度与扩展的Sobel算子进行卷积运算, 生成温度梯度图; 根据经验值, 设置梯度下界为[(gm+0.5σ g), (gm+1.5σ g)](σ g为梯度的标准差, gm为梯度均值), 梯度上界为(gm+3.2σ g), 在[(gm+0.5σ g), (gm+3.2σ g)]到[(gm+1.5σ g), (gm+3.2σ g)]之间, 分别以0.1σ g进行递增, 形成11景潜在的高梯度缓冲区影像, 绝大多数的潜在高梯度缓冲区中心与火区中心重合; 然后, 使用经典的数学形态学细化算法, 将潜在高梯度缓冲区骨架化成仅有1个像元的高梯度极值线; 为了排除极少数位于低温区的高梯度极值线, 设置经验阈值tm+σ t(σ t为温度的标准差, tm为温度均值)将反演的温度图像分割出高温缓冲区; 叠置高梯度极值线与高温缓冲区, 统计两者重合部分的温度值并取其均值, 作为该景潜在高梯度缓冲区影像的火区温度阈值; 最后, 为进一步提高精度, 取11景影像的温度阈值的均值作为最终的分割火区的温度阈值。该算法程序在IDL8.3下实现。
以2008年11月16日获取的Landsat 5影像为例, 采用单窗算法反演的地表温度如图1(a)所示; 用SAGBT算法提取煤田火区的过程包括: ①根据单窗算法反演的地表温度生成温度梯度(图1(b)); ②以[(gm+0.5σ g), (gm+3.2σ g)]为下界、[(gm+1.5σ g), (gm+3.2σ g)]为上界, 形成11景潜在高梯度缓冲区(图1(c)); ③通过数学形态学细化算法将潜在高梯度缓冲区骨架化, 生成高梯度极值线(图1(d)); ④通过叠置高梯度极值线与蓝色高温缓冲区(图1(e))计算该景影像的平均温度值; ⑤最后计算11景影像的平均温度值作为最终的温度阈值, 划分出煤田(分为A区―F区共6个子区)内的火区(图1(f))。
为了验证上述探测火区方法的可靠性并估算其精度, 笔者采用了2种对比方式: ①卫星影像识别火区与航空影像识别火区的对比; ②卫星影像识别火区与野外调查火区位置的对比。
首先, 对使用2008年11月16日Landsat 5 TM热红外波段提取的火区与利用2009年3月航空三通道影像调查的火区进行叠置对比分析(如图1(f)所示), 发现2009年3月航空影像提取并经野外调查所圈定的2-1号、2-2号、4号、6-2号、15号、16-2号和17号共7个火区未能从2008年11月16日的卫星影像中识别出来, 未识别火区中的17号和15号火区稍大, 其余的规模较小; 其他14个火区均被正确识别, 被航空影像和地面调查所证实的火区占总火区的67%。通过分析认为, 7个火区未能探测出来的原因极有可能是因为本文所用Landsat卫星数据的过境时间为2008年11月16日, 比2009年3月航空调查数据获取时间早4个多月, 而煤田在2008年— 2010年8月正式实施灭火前, 火区面积呈持续递增趋势, 因而在获取时间较迟的影像中可以识别出更多的火区。笔者推测, 如果选取更接近航空影像拍摄日期的卫星影像(如使用2009年3月获取的Landsat 5数据)来识别火区, 估计通过卫星影像和航空影像所探测到的火区会有更高的重合率。
然后, 进行卫星热红外影像识别火区与野外实测火区的位置对比。将使用2015年3月25日Landsat 8 TIRS热红外波段提取的火区与2015年8月现场调查的部分明显烟点进行叠置对比分析(图2(c))后发现, 8个调查的烟点中仅有A区东南角的2个火点未与探测的火区重叠, 其余6个火点均在探测的火区范围内, 估算的火区探测准确率为75%。综上所述, 可认为本文所用的火区识别方法有大于70%的准确率, 有较高的可信度。
由于基于梯度的自适应阈值法对温度影像的反演精度依赖较低, 对火区的识别主要与影像热分布及温度快速衰减边缘有关, 因而采用该方法为不同时期火区的变化检测提供了可能。在方法可靠性得到了验证的前提下, 本文使用4景2008— 2015年间获取的热红外影像进行了乌达煤田内火区的动态变化分析。
为了便于统计4个年度的火区变化情况, 结合2009年3月调查的火区边界、苏海图矿界、黄白茨矿界以及五虎山矿界, 将煤田划分为A区、B区、C区、D区、E区和F区共6个子区(图2)。通过ArcGIS中的Summary Statistics功能统计每个年度各子区内的火区面积, 根据式(7)计算火区面积演变的相对变化率vS, 即
vS=
式中: Sf和Se分别为变化前和变化后的火区面积。
表3示出2008―2015年间乌达煤田火区演变状况。
由表3可知, 2008―2011年间, A区、B区、C区和F区的火区范围出现不同程度的缩小(但D区和E区的火区范围增大, 出现新火区, 火势增强), 整个煤田内火区总面积由2008年的1.951 km2减少到2011年的1.194 km2, 降幅达38.80%。这说明从2010年10月起实施的灭火工程取得了较好的初步成效。2011— 2013年间, C区和D区的火区出现明显萎缩, 火势减弱; A区、B区和E区的火区有少量增加; 而F区的火区扩张且离散地分布在整个F区, 整个煤田火区总面积下降到了1.092 km2, 降幅为8.54%, 相对于2008―2011年有所减缓。2013―2015年间, A区的火区面积稍有增长, 零星火区减少, 火区变得相对集中; F区的火区连接成片成为大火区; B区、C区、D区和E区的火区急剧萎缩, 火区数量减少且相对集中; 火区面积继续下降到0.873 km2, 降幅为20.05%, 火区面积仍在缩小。综上所述, 2008―2015年间, 乌达煤田内火区面积在持续缩减, 灭火工作初见成效; 但F区的火区有连通聚集成片的蔓延趋势, 应引起有关部门的重视。
1)本文基于Landsat系列卫星数据, 利用单窗算法与自适应梯度阈值(SAGBT)法提取了乌达煤田火区。通过将2008年11月16日卫星影像提取的火区与2009年3月航空手段圈定的火区叠置比对, 发现重合率达67%; 根据2015年8月份现场调查的部分烟点与2015年3月提取的火区叠置比对, 表明有75%的烟点位于所提取的火区内。通过验证认为, 该算法的有较高的准确率, 可满足火区变化检测与动态分析的要求。
2)在2008―2011年间, 发现除2个火区面积增加外, 整体火区面积降幅达38.8 %; 2011―2013年间, 除1个火区面积有明显增加外各火区火势减弱, 火区面积整体降幅为8.54%, 较上个时段有所减缓; 2013―2015年间, 各火区变得相对集中, 火区整体面积降到0.873 km2, 降幅为20.05%。纵观2008―2015年间火区演变过程, 火区由治理前的相对集中, 变为治理期间的离散零星分布, 再到治理后期的集中分布, 全程反映了火区治理的动态变化过程。火区面积由治理前的1.194 km2下降到0.873 km2, 表明内蒙古乌达煤火治理实现了26.88%的减量, 但仍有73.12%的存量需要治理, 任务仍然十分艰巨。
本文利用Landsat卫星影像提取的火区是指由煤火导致的地表热异常区(包括高温燃烧中心和烟气释放高温区), 而非温度异常不能到达地表的地下煤层缓慢自燃区。此外, 太阳高度角、地形起伏、土地利用类型、气象条件等环境因素也影响着Landsat卫星反演地表温度的精度, 进而会影响对火区的识别。在后续的研究中, 采用高光谱分辨率和高空间分辨率卫星数据联合地表物化探手段, 可有效提高火区识别的精度。未来随着无人机热红外遥感技术的发展和成熟, 实时、廉价的遥感热红外火区探测, 不仅会提升火区探测的精度, 而且也会为开展持续的煤火动态监测提供可能。
志谢: 感谢中国地质调查局发展研究中心杜晓敏博士给予本文的帮助!
The authors have declared that no competing interests exist.
[1] |
|
[2] |
|
[3] |
|
[4] |
|
[5] |
|
[6] |
|
[7] |
|
[8] |
|
[9] |
|
[10] |
|
[11] |
|
[12] |
|
[13] |
|
[14] |
|
[15] |
|
[16] |
|
[17] |
|
[18] |
|
[19] |
|
[20] |
|
[21] |
|