基于D-InSAR的济宁市采煤沉陷监测与分析

2019-02-13姬宗皓薄怀志王东平高彦生宋炳忠

山东国土资源 2019年2期

姬宗皓，薄怀志，王东平，高彦生，宋炳忠

(1.济宁市采煤塌陷地治理中心，山东济宁 272000；2.自然资源部采煤沉陷区综合治理与生态修复工程技术创新中心，山东济宁 272100；3.山东省鲁南地质工程勘察院，山东济宁 272100；4.山东省地矿局采煤塌陷地综合治理工程技术研究中心，山东济宁 272100)

0 引言

济宁是鲁西煤炭基地的重要组成部分，煤田面积3920km2，主要由兖州煤田、济宁煤田、宁汶煤田、滕县煤田、金乡煤田、梁山煤田、巨野煤田组成，煤系地层主要为石炭-二叠纪月门沟群太原组和山西组[1]。经过几十年高强度开采，济宁出现了大面积的采煤沉陷。采煤沉陷造成严重的地质环境问题和社会问题，主要表现在：耕地保有量锐减、严重破坏矿区生态环境、引发诸多复杂社会问题、严重制约城市建设发展等[2-3]。查清济宁市采煤沉陷现状，对于采煤塌陷地监测监管意义重大[4]。与土地利用现状变更不同，开采沉陷主要表现为垂向空间变化，大范围监测难度较大。目前对于单个开采工作面的沉陷监测多使用传统的测量方法，即布设沉陷监测点，使用水准仪或GNSS-RTK定期监测沉陷情况[5]。该方法测量精度较高，但观测周期长、费用高、观测点离散分布难以反映大范围全面形变情况和连续形变规律。秦晓敏等[6]提出利用多时相DEM数据高程变化、结合高分辨率遥感影像、矿区开采资料提取开采沉陷范围的方法，但精度较低，难以监测沉陷量小于0.5m的塌陷地。

D-InSAR技术具有全天候、大范围、精度高等优点，近年来发展迅速，已经在地面沉降监测、地质灾害防治等领域得到了广泛的应用[7-9]。将D-InSAR技术应用于济宁市大范围采煤沉陷区动态监测，对开采沉陷监管、治理规划设计、政府决策等方面都具有重要意义。RADARSAT-2是一颗搭载C波段传感器的高分辨率商用雷达卫星，由加拿大太空署与MDA公司合作，于2007年12月14日在哈萨克斯坦拜科努尔基地发射升空。该文使用RADARSAT-2的数据产品SLC(单视复型数据)，研究了D-InSAR技术在济宁市大范围采煤沉陷区监测及其时序变化特征分析中的应用。

1 D-InSAR采煤沉陷监测的原理

D-InSAR技术是由InSAR技术发展而来的，它是以合成孔径雷达复数图像的相位信息获取地表变化信息的技术。根据成像时间分类，InSAR可以分为单次轨道和重复轨道2种模式。单次轨道干涉是指在同一机载或星载平台上装载两幅天线，其中一幅天线发射信号，两幅天线都接受地面回波信号，并利用获取的数据进行干涉处理。重复轨道干涉是指同一传感器或相似传感器按照平行轨道2次对地成像，分别发、收信号，利用得到的数据进行干涉处理。目前，常用的D-InSAR地表形变监测通常为星载单天线重复轨道模式。

1.1 InSAR测量原理

假设雷达传感器S1、S2在不同时间、两次相互平行地对同一区域重复观测，生成两幅影像(主影像、从影像)，S1、S2的距离为空间基线，按雷达信号入射方向，基线B可投影为平行分量和垂直分量，θ为雷达侧视角，α为空间基线水平倾角，传感器S1到地面分辨单元P的斜距为r，与S2到该点的斜距差为Δr，如图1所示。P在主(从)影像上的灰度、相位值可表示为：

主影像和从影像经配准后共轭相乘生成干涉图，该分辨单元的干涉相位值为：

ψ=arctan(u1·u2*)=ψ1-ψ2ψ∈[-π,π]

SAR系统观测值ψi(i=1,2)与雷达波至地面单元的距离有关，同时受地面分辨单元散射特性的影响。假定2次成像时的地表散射特性相同，则干涉相位可表示为：

图1 InSAR几何示意图

1.2 D-InSAR沉陷监测基本原理

如图2所示，两次成像期间地表像元P由于地下开采发生了沉陷ΔH,且沉陷量在雷达视线方向上的形变分量为Δd，为了获取形变相位，需将地形相位从干涉相位中去除，采用外部DEM差分之前先地理编码至雷达坐标系，与SAR影像精密配准，其模拟的地形相位[10-12]：

从干涉相位中去除模拟地形相位后，忽略大气延迟及噪声误差，得到形变相位：

当形变相位变化2π时，对应P的形变为：

d2π为形变模糊度，与SAR成像波长有关。获取P点在雷达视线方向上的形变分量Δd后，即可根据ΔH与Δd的几何关系计算得到P点的沉陷量ΔH。

图2 D-InSAR几何示意图

2 采煤沉陷区范围的提取

2.1 InSAR数据的选择

根据开采沉陷的特点，综合考虑存档数据的覆盖情况、时间序列完整性、空间基线等因素，结合济宁煤炭开采范围，选取15期RADARSAT-2 Wide模式数据进行开采沉陷监测，详细数据参数见表1[13-14]。

表1 数据获取参数

2.2 数据处理流程

根据以往InSAR数据处理的经验，采用D-InSAR多基线累积结果作为济宁采煤沉陷的形变监测结果。多基线干涉相位累积叠加方法以两轨差分为基础，选取相邻影像组合干涉对，即前一个干涉对的辅影像为后一个干涉对的主影像，组成首尾相连的干涉对组合，对每一个干涉对采用两轨法进行处理得到解缠相位，将各个干涉对的解缠相位连续累加并转化为形变量，由此求取的形变即是对应观测时段内的累积形变量。该方法利用一系列解缠的差分干涉图，通过累积的方法得到相应的形变量。处理时对独立干涉像对大气相位进行了平均，达到抑制大气噪声的目的。

D-InSAR数据处理主要包括主辅影像预处理、影像配准及重采样、干涉图生成、基线估计、地形相位差分、差分图滤波、相位解缠、地理编码等内容[15-16]。

由于采用SAR数据为C波段成像，时间去相干比较严重，为保证干涉的质量，对满足时间基线小于200天，空间垂直基线小于300m的34个干涉对进行干涉、地形相位差分、滤波、解缠等批处理。所选干涉对的时间、空间垂直基线见表2。

表2 干涉对的时间、空间垂直基线数据

为减小SAR成像系统噪声，在基于强度的SAR影像配准过程对距离向、方位向进行1∶5多视，以20140617为主影像进行影像配准，采用强度互相关的多级匹配方法，距离向及方位向配准精度均在0.05像元以内，满足干涉要求。主辅影像重采样后进行干涉处理，得到相干图。采用分辨率30m的SRTM1 DEM，结合卫星星历数据、成像几何参数等与主影像强度图进行精配准、重采样，配准精度优于0.1像元(距离向0.06，方位向0.08)，由精配准并编码至SAR坐标系的高程进行地形相位模拟，并从原始干涉图中去除地形相位，得到差分干涉图，如图3所示[17]。

图3 34个干涉对对应的差分干涉图

经解缠处理恢复真实的差分相位。解缠参考点位于稳定区域，采用相干性阈值0.3～0.4对低相干区域掩膜后使用最小费用流法整块解缠得到初始解缠相位图。解缠前只对少数几个差分图利用条纹频率法进行基线估计，可以看出大部分解缠相位中含有轨道残余相位，除此还存在几处解缠误差，由局部相干性过低造成；采用二次曲面拟合的方法去除轨道残余相位，并对地表水汽相位采用高程相关的线性模型进行拟合，最后采用最邻近插值法对掩膜区的解缠断裂进行空间插值，得到最终解缠相位，如图4所示。

图4 34个干涉对对应的最终解缠图

2.3 结果分析

采用多基线累积叠加结果作为该次矿区沉陷的InSAR监测结果。根据多基线累积叠加方法相关分析，利用未插值累积叠加结果用于提取超过InSAR形变监测能力的“空值区”，“空值区”包含在5cm范围线内；利用插值累积叠加结果提取塌陷5cm范围线。同时，结合现场调查、矿山采掘工程平面图、概率积分法沉陷预测确定塌陷地范围[18](图5)。图中黑色多边形区域为矿区范围线，在矿区范围内监测到多个采空区的塌陷活动，采空区塌陷在空间上表现为漏斗形状，与矿区开采引起的地表塌陷空间特征符合的较好。

图5 基于累积形变提取的塌陷5cm范围线

2.4 精度分析

为了验证D-InSAR开采沉陷监测的可靠性，选取了某煤矿岩移观测数据与D-InSAR监测数据进行对比。为了真实反映D-InSAR精度，需保证InSAR数据获取时间与岩移观测时间基本一致，共选取4期数据、150个岩移观测点进行对比，结果见表3、图6。

表3 D-InSAR精度检核

图6 精度检核较差范围图

3 采煤沉陷区时序变化特征

多基线累积叠加在得到矿区总体形变信息的同时，也可很好地反映矿区形变发生过程[19]。以花园煤矿为例，分析矿区地面沉陷时序变化情况。该处相干性较好，沉降特征明显，并且沉降中心较为连续。区域累积沉陷结果见图7。

图7 区域累积沉陷结果图

分析其沉陷过程，2013年4月—6月间矿区形变较为稳定；2013年6月—9月出现形变加剧的趋势；2013年9月—2014年6月间形变较为平稳；2014年6月—9月间形变进一步加剧；2014年9月—2015年5月间形变平稳；2015年5月—10月，形变趋势进一步加剧，沉陷呈现为非线性逐步增大的过程。结合地质采矿及工作面布设情况，经分析，A1点周期性间隔沉陷主要是受到临近工作面开采沉陷影响。A1点沉陷时序变化曲线见图8、区域沉陷时序变化见图9。