杨 洁，曾 强

(1.新疆大学 资源与环境科学学院，新疆 乌鲁木齐 830046； 2.新疆大学 干旱生态环境研究所，新疆 乌鲁木齐 830046； 3.绿洲生态教育部重点实验室，新疆 乌鲁木齐 830046)

煤炭开采过程中，煤易发生自燃，进而导致地下煤火的发生。地下煤火的燃烧会影响煤矿及火区周边地区的社会和谐和可持续发展，地下煤火已成为煤矿开采区的主要灾害之一[1-4]。干旱炎热气候条件、易自燃煤炭属性及落后的开采方法，使新疆地区煤火不断发生[5-6]。燃烧火区沉降会导致地表产生裂隙，土壤结构退化，不断加剧水分蒸发，并对植物根系造成严重的机械损伤，使本已受到破坏的生态环境变得更加脆弱[7-10]，煤火造成的裂隙既是煤火燃烧的供氧通道，也是烟气逸散的通道。火区地表沉降与裂隙变化在一定程度上可反映火区的动态变化，以及火区空气进入和烟气逸散的特征。有学者研究表明：通过分析研究区干涉图，可以得到地下煤火燃烧导致的新沉降区及地下煤火的燃烧和地表面的沉降相互影响规律。其中合成雷达差分干涉(Differometric Synthetic Aperture Radar，D-InSAR)技术[11-12]是得到干涉图的有效方法之一。美国学者Gabriel等[13]使用D-InSAR技术监测地表形变，获取SEASET L波段影像数据经差分干涉处理后得到了研究地区的形变图，验证了该技术在灌溉区沉降监测的可行性；Donald等[14]利用机载干涉合成孔径雷达系统的数字地形DEM数据，对断层崖进行分析并获得地形地貌；李旺[15]、吴立新[16]等分别用 D-InSAR 技术对矿区、火山地区进行沉降监测，得到了研究区形变图，这些研究成果很好地指导了实践。

笔者采用双轨差分干涉测量方法对准南煤田水西沟火区的Sentinel数据进行处理，提取火区的地面沉降变化信息，尝试确定火区空气进入和烟气逸散区域，为监测火区的动态变化提供一种可行的途径和方法。

1 D-InSAR技术原理

D-InSAR技术是以雷达干涉测量理论为基础，通过复数据提取的相位信息为信息源获取地表三维信息和变化信息的一项技术，其基本原理如图1所示。

图1 D-InSAR基本原理

图1中：S1和S2分别为主、辅图像传感器；B为基线距离；α为基线与水平方向的倾角；θ为主图像入射角；H为主图像传感器相对地面高度；R1和R2分别为主、辅图像斜距；ΔR为主辅图像斜距变化量；P为地面目标点；h为高程；B∥和B⊥分别为基线B在斜距向平行和垂直方向上的投影分量。

雷达干涉相位φint的组成如下[17]：

φint=φtop+φdef+φflat+φatm+φnoi+2kπ

(1)

式中：φtop为地形起伏引起的地形相位；φdef为地表引起的形变相位；φflat为参考椭球面引起的参考相位，即平地相位；φatm为大气延迟引起的延迟相位；φnoi为处理误差及热噪声引起的噪声相位；k为整周模糊度；π为周期。

各相位分量见式(2)～(5)：

(2)

(3)

(4)

(5)

式中：λ为波长；Δh为目标物两次成像的形变量；Δr为沿雷达视线向的地面形变量；R为卫星到地面的斜距；r为外部DEM高程；δr为两个成像时刻在斜距方向上的大气延迟。

由式(1)～(5)可知，将其余相位去除可得到地表形变信息。

2 实验数据选择

研究区域位于准南煤田水西沟火区，位于东经88°55′～88°58′，北纬43°55′～43°57′之间。水西沟火区北部表现为20～30 m深的塌陷状态，而火区南部地表反潮且露头区有芒硝状结晶；火区燃烧煤层为6#、7#煤层，每年燃煤损失量为12.96万t[8]。采用NASA发布的90 m SRTM数据作为外部DEM数据用以消除干涉相位图中的地形因素影响，并获取了14景Sentinel卫星Level-1产品中单视复数图像(Single Look Complex，SLC)数据(https://earthdata.nasa.gov/)进行差分干涉处理。Sentinel-1(哨兵1号)卫星系统及数据主要参数见表1。

表1 Sentinel-1卫星数据主要参数

3 数据处理流程

使用双轨D-InSAR技术对2014年10月7日至2015年10月3日共14景单视复数图像(SLC)数据，进行差分干涉处理得到水西沟火区形变场，并提取分析了水西沟火区地表变化信息。数据处理流程如下：

1)获取干涉图：经配准、干涉、干涉图滤波、精密基线估计等步骤实施后得到干涉图。

2)模拟地形相位：将SRTM DEM数据与主图像进行配准，转换到雷达坐标系下，从而模拟得到地形相位。

3)获取形变量：选择主、辅图像，利用参数文件进行估算计算基线。

4)滤波处理：选用Goldstein自适应性滤波方法进行滤波处理。

5)相位解缠：采用最小二乘法，将图像划分为若干个小正方形格网，对图像中的所有区域所有像元进行处理，对相干性过低且小于阈值的部分掩膜进行处理，将干涉相位和干涉图的成像几何关系联系起来从而获得地表高程。

2014年10月7日至2014年10月31日经过差分干涉处理的干涉图、差分干涉图、去平地效应图和相位解缠效果图见图2。

(a)干涉图

(b)差分干涉图

(c)去平地效应图

(d)相位解缠效果图

2014年10月7日至2014年10月31日解压后的干涉图如图2(a)所示；差分处理后得到的差分干涉图如图2(b)所示，图中可看到上部有明显的干涉条纹，在其他部分表现不是很明显；经过去平地效应之后去除一部分噪声得到图2(c)，图中条纹更加清晰可见，线条更加完整；相位解缠处理得到解缠效果图如图2(d)所示，整幅图像均出现了条纹，干涉信息得到进一步处理。

4 差分干涉处理和结果分析

差分干涉对时间间隔和垂直基线、入射角及中心经纬度参数如表2所示。

表2 干涉对参数

由表2可知，除2015年4月18日至2015年 6月 5日干涉对时间基线为48 d外，其余均为24 d，而2014年10月31日至2014年11月24日空间基线最长为213.27 m。基线长度是相干的一个重要限制因素，垂直基线过长可能导致失相干，过短效果不明显。每一幅图像中心经纬度相差不大，均覆盖研究区域。

按照上述数据处理流程共得到11组干涉对，通过ArcGIS软件对得到的沉降量图像进一步分析处理，获取该研究区域在各时间段内的具体地面沉降分布，其变化如图3所示。

图3 地面沉降分布变化

由图3可知，2014年和2015年Sentinel数据经过差分干涉处理生成的地面形变图中沉降范围大体一致，说明了双轨法可以从整体上确定矿区的沉降分布情况。引入外部DEM数据去除相关误差，得到具体沉降范围如表3所示。根据表3数据进一步绘制干涉对沉降变化折线图，如图4所示。

表3 2014—2015年干涉对沉降分布范围

图4 不同干涉对沉降变化折线图

由表3和图4可知，不同干涉对的时间间隔比较均匀(14 d或14×2 d)，并且沉降变化明显。由图4 可知，随着影像干涉对的时间推移，研究区沉降高低值呈现相似变化规律。2014年10月7日至2015年4月18日期间，沉降高、低值呈波动变化；随后于2015年6月29日分别减少至-21.25 mm(沉降高值)、-84.19 mm(沉降低值)；2015年6月29至2015年10月3日总体呈现增长趋势。

5 融合遥感综合研究

通过单窗算法和植被模型反演2014年遥感影像数据得到地表温度和植被覆盖图像。利用ArcGIS软件重分类计算得到高温区域和不同级别植被覆盖区域(将植被盖度五等分，其中[0,20%)为一级植被覆盖区域，[20%,40%)为二级植被覆盖区域，[40%,60%)为三级植被覆盖区域，[60%,80%)为四级植被覆盖区域，[80%,100%]为五级植被覆盖区域)影像图，叠加分析高温区域和一级植被覆盖区域初步得到火区裂隙可能存在的位置，火区裂隙圈定如图5所示。

图5 火区裂隙圈定图

结合地表温度和植被覆盖场去接受或者拒绝一个潜在的热异常和植被异常像素，从而圈定出A、B、C、D、E 5个火区裂隙范围。通过shp图像提取A、B、C、D、E 5个范围2014年的沉降变化，具体信息见表4。

表4 2014年沉降变化信息

由表4数据进一步绘制所圈定的5个火区位置沉降变化折线图，地表沉降最大值变化折线如图6所示，沉降最小值变化折线如图7所示。

图6 地表沉降最大值变化折线图

图7 地表沉降最小值变化折线图

由图6可知，位置A、B、C地表沉降最大值随着时间变化一直减小；位置D稍有增加后骤减；位置E先明显增大后减小。

由图7可知，位置A地表沉降最小值变化不明显；位置B、C均随时间变化先明显减小后增大，位置C表现为沉降，无隆升量；位置D先明显减小后略减小；位置E先增大后显著减小。

已有Sentinel数据处理后表明2014年研究区裂隙可疑位置沉降趋势增加明显：2014年10月 7日至2014年10月31日干涉对仅位置C出现沉降，其他位置地表均为隆升；2014年12月7日至2014年12月31日干涉对5个位置沉降最小值及均值均为负值，表示发生沉降。C位置沉降幅度变化最小，相对较稳定，但其沉降结果为负值表示沉降。

6 结语

1)采用双轨法差分干涉处理方法测量水西沟火区Sentinel数据，验证了差分干涉测量(D-InSAR)技术用于煤田火区地面沉降监测的可行性。

2)遥感反演处理得到水西沟火区温度异常区及一级植被覆盖叠加图，可确定火区裂隙可能存在A、B、C、D、E共5个区域，为远程宏观监测火区裂隙区变化提供了实用方法，为进一步分析火区燃烧系统动态演化提供了可行途径。