摘要：煤矿高强度开采会导致严重的地面形变及次生地质灾害。时序合成孔径雷达干涉(synthetic aperture Radar interferometry, InSAR)具有较强的形变监测能力,但在开采核心及周边低相干区域时无法监测到足够的目标点。该文尝试将分布式目标(distributed target, DT)和缓慢去相关滤波相位目标(slowly-decorrelating filtered phase target, SDFPT)进行联合,以提高矿区形变监测点的密度和覆盖度。分别采用快速同质点选取(fast statistically homogenous pixel selection, FaSHPS)法和振幅离差指数法选取DT和SDFPT候选点,分别对2类点进行相位优化和稳定性分析,筛选出符合条件的DT和SDFPT形成融合点集,并对其进行三维相位解缠、恢复相位时间序列和时空滤波,最终得到融合点集的形变时间序列和年均形变速率。选取2018年4月—2020年4月获取的覆盖布尔台煤矿的60景Sentinel-1影像进行形变监测,结果表明,融合DT和SDFPT后形变点密度和覆盖度显著提升,可监测最大形变量级也随之增加。实验区域内存在5处形变漏斗,最大累积形变量达到-309.76 mm; 形变影响范围和不同年份时序形变量的差异与矿区开采活动密切相关。
引言
煤炭是我国主要的能源和化工资源之一。统计资料表明,煤炭资源在我国能源结构中占比接近60%。煤炭资源的开采极大地带动了国民经济的发展,但长期的开采破坏了矿区及周边的底层结构和生态环境,诱发了地面形变、滑坡等地质环境问题,严重危害到矿区及周边人民的生命财产安全。开展矿区地表形变监测及分析是确保煤矿资源稳定开采、减少人力财力损失的必要举措。
以永久散射体(persistent scatterer, PS)干涉和小基线集(small baseline subset, SBAS)干涉为代表的时序合成孔径雷达干涉(synthetic aperture Radar interferometry, InSAR)技术在地表形变监测领域中得到了广泛的应用。其主要思想是从时序合成孔径雷达(synthetic aperture Radar, SAR)影像中选取相位稳定的相干目标点,并提取差分干涉相位,从中分离出地形残余相位、大气延迟和轨道误差相位、噪声相位以及所需的形变相位,形变监测精度最高可达毫米级。PS-InSAR和SBAS-InSAR显著提高了地面形变监测精度,但对于植被覆盖率较高、人工目标较少的区域,获取的地物目标数量有限,从而无法完全反映研究区域的整体形变。
为了提升此类区域的监测点密度和覆盖度,国内外学者将研究重点转移到自然界中散射特性中等且集群式分布的分布式目标(distributed target, DT)。联合DT与PS的新一代分布式目标InSAR (DT-InSAR,或分布式散射体InSAR (distributed scatterer InSAR, DS-InSAR))方法应运而生。该方法通过同质检验选取DT并对其相位进行优化,然后将DT和PS进行联合处理。研究表明,DT在数量、密度和覆盖度方面远高于PS,显著提高了形变监测效果。DT-InSAR逐渐在铁路、滑坡、矿区等多个领域的形变监测中得到较多应用。在传统的DT探测思路基础上,蒋弥等提出一种快速同质点选取(fast statistically homogenous pixel selection, FaSHPS)方法,相较于其他同质点选取方法计算效率更高,且在影像数量较少的情况下也能有效地对同质点进行识别。
相比于一般的监测区域,煤矿开采区的形变梯度较大,加重了区域失相干的程度。因此,PS-InSAR和SBAS-InSAR在矿区形变监测中很难取得令人满意的效果。DT-InSAR在矿区形变监测中的应用在很大程度上弥补了PS-InSAR和SBAS-InSAR的不足。现有研究均是将DT和传统的PS目标进行融合,然后采用单一主影像或者SBAS干涉组合策略进行相位提取和形变建模。根据现有研究可知,PS点密度在非城市区域非常低,尤其是在矿区形变漏斗区域,基本无法有效识别出PS目标。因此,对于矿区形变监测,DT-InSAR监测点的密度有待进一步提高。
2008年,Hooper在其所提出的StaMPS方法的基础上,提出了一种联合PS和SBAS的形变建模和解算思路,并提出了缓慢去相关滤波相位(slowly-decorrelating filtered phase, SDFP)的概念和处理思路,显著提高了形变监测点的数量和密度。传统的PS探测采用较为严格的振幅离差指数(amplitude deviation index, ADI)阈值,而Hooper则首先采用较宽松的ADI阈值获取目标候选点,然后对低质量候选点进行相位滤波提高相位质量,与高质量候选点共同参与后续的形变建模与解算。研究表明,该方法比传统的PS,StaMPS和SBAS干涉具有更好的形变监测效果,尤其是在失相干区域具有明显的应用优势。
本文以全球单井开采量最大的井工煤矿——中国布尔台煤矿局部及周边区域为实验区域,尝试在时序InSAR中将DT和SDFP目标(SDFP target, SDFPT)进行融合,分别通过FaSHPS和ADI法选取DT和SDFPT候选点并进行相位优化,以期进一步提高形变监测点的密度和空间覆盖度,提升矿区形变监测效果。选取覆盖实验区域的60景Sentinel-1影像数据,利用融合多类型相干目标的时序InSAR方法进行数据处理和分析,并与StaMPS-SBAS监测结果进行对比,验证该方法的可靠性及监测效果的提升程度。基于监测结果,对布尔台矿区的地表形变进行分析和解释,为矿区形变区域的灾害防治提供数据支撑。
1 研究区概况与数据源
1.1 研究区概况
研究区域位于内蒙古鄂尔多斯市内伊金霍洛镇东部及乌兰木伦河西南部,地处神府东胜(神东)煤矿群中的布尔台煤矿、寸草塔二矿和霍洛湾煤矿一带。具体的监测范围如图1所示。该区域地貌主要以黄土丘陵为主,地势总体呈现出高态势。气候为大陆性季风气候,降雨量较少,主要集中在每年的7,8,9这3个月份。由于受干旱和人类活动等因素的长期影响,植被覆盖率相较于其他地区不高,主要以农作物和杂草为主。由于近些年来开采活动加剧,导致地表形变影响面积扩大,对矿区周围环境有很大的安全隐患。
近年来,一些学者采用传统的差分InSAR (differential InSAR, DInSAR),SBAS-InSAR或DS-InSAR获取了神东煤矿不同矿区2012—2013年及2017—2019年间的形变分布。研究表明,该区域存在显著的形变漏斗,常规的DInSAR和SBAS-InSAR在矿区形变监测中效果并不理想,DS-InSAR具有较好的监测效果。但是,该区域内缺少后续的持续跟踪研究。并且,目前采用时序InSAR的研究没有覆盖本文所选实验区域(尤其是布尔台煤矿)。布尔台煤矿是目前世界上单井开采量最大的井工煤矿,煤层呈现出近水平和浅埋藏的特点,工作面的采煤平均厚度为2.5 m,采深为265 m,具有较高的监测必要性和研究价值。
1.2 数据源
SAR影像数据来源于Sentinel-1卫星。该卫星是欧洲航天局哥白尼计划的地球观测卫星,主要由Sentinel-1A和Sentinel-1B这2颗卫星组成,载有C波段SAR传感器,波长约为5.6 cm,卫星采用近极地太阳同步轨道,轨道高度约700 km,双星协同工作下,同一地区的重复观测周期为6 d。本文选取实验区域2018年4月—2020年4月的60景干涉宽幅模式SAR影像为数据源。
采用由日本宇宙航空研究开发机构(Japan Aerospace Exploration Agency,JAXA)提供的具有30 m空间分辨率、5 m高程精度的AW3D30 DSM数据作为外部地形数据。在本文中,结合StaMPS小基线干涉方法及矿区形变的特点,干涉时空基线阈值的设置需要考虑2方面的因素: 一是矿区的形变量级较大,当时空基线较大时,干涉对的失相干较为严重,且形变相位梯度也较大,因而进一步增加恢复相位模糊度的难度,导致难以准确解算研究区的形变; 二是StaMPS小基线干涉方法采用最小二乘恢复相位时间序列,需要干涉图集的时空基线连续分布,以避免恢复相位时间序列时出现秩亏,而过小的时空基线限制会导致干涉图集不连续,进而导致秩亏问题。为了降低时空失相干、相位梯度对监测结果的影响,并保持干涉图集的连续性,本文进行了多次实验,当设置空间基线阈值为80 m,时间基线阈值为48 d时,可满足上述2个条件,取得较好的形变监测结果。干涉图集的时空基线如图2所示。
2 数据处理方法
2.1 DT探测及同质滤波
DT点的选取主要分为2个步骤,同质点选取和时序相位优化。同质点的选取的基本思想是判别SAR影像数据集中一定空间范围内的2个像元是否属于同一地物,算法的主要原理是抽取同一像元在时间维度上的强度信息,度量2个样本的相似程度,从而判断是否为同质点(DT候选点)。
本文对于同质点选取使用的是FaSHPS算法,该算法不同于传统的假设检验方法,其核心思想是转换为置信区间估计,通过简单的逻辑运算来识别同质点。将三维数据在时间维求得一个平均值,目标像素的平均值作为真值,而目标像素的邻近像素值作为待估值。
2.2 SDFPT探测
基于SBAS-InSAR的技术原理,Hooper在2008年的基础上进行了改进,提出了StaMPS-SBAS方法,对相干点的选取方法进行了优化,目标像素点被称为SDFP点。
2.3 DT与SDFPT融合及形变解算
前已述及,SDFPT候选点的探测采用了ADI阈值法。其中,
ADI
≤0.25的像元被标记为PS点,
ADI
>0.25且小于某个阈值(文献[26]中建议采用0.6)的像元被标记为SDFPT候选点。为保证高质量PS点相位不受影响,不做任何处理,直接参与后续计算。
2.4 数据处理流程
数据处理流程如图3所示。对Sentinel-1A数据和DSM数据进行地理编码、主辅影像配准和差分干涉等预处理后得到差分干涉图和强度图,采用FaSHPS方法选取得到DT候选点集合,对DT候选点的差分干涉图进行同质滤波后,提取出候选点的像素坐标和相位,再进行时序相位优化后输出像素坐标和相位矩阵。同时,根据振幅离差指数选取出SDFPT候选点,迭代分析候选点的相位稳定性后设定时间相干性阈值后得到SDFPT点。
图3 数据处理流程
将探测得到的DT点和SDFPT点融合得到目标点集合,后续处理以StaMPS-SBAS方法的流程为基础,对目标点集合进行三维相位解缠,对解缠结果进行最小二乘解算,计算出空间不相关侧视角误差,并去除大气延迟、轨道误差和噪声等,从而得到LOS向形变时间序列和形变速率。
3 结果与分析
3.1 DT点质量分析
后验相干性及时序相干性(相干矩阵)是衡量DT目标质量的重要指标,二者的取值范围均为[0,1]。如图4所示为本文中研究区域的后验相干系数和时序相干系数。后验相干性的空间分布表明,本文研究区域内大部分自然地表具有较高的相干性,该区域内目标点的整体质量较好。其中,后验相干性超过0.5的目标占80.28%。对时序相干性矩阵进行统计表明,任意干涉对的平均相干性均大于0.25,表明经过相位优化后干涉数据集的整体干涉质量较好。为保证所筛选的DT具有较好的质量,建议选取后验相干性值大于后验相干性取值范围一半的目标作为DT。同时,为保证具有足够多的DT参与后续形变建模和解算,提高监测点密度和覆盖度,后验相干性阈值不宜设置过大。本文中,以0.6作为DT筛选的后验相干性阈值。
图4 数据集质量评价
图5
展示了20190407—20190501原始差分干涉相位和优化后相位的对比。可以看出,在低相干区域,采用时序相位优化策略可利用时间维信息有效实现噪声区域信噪比的提升,可以使空间相位更加连贯,进而使整个干涉图集的相位质量得到提升。综上所述,本文中所筛选出的DT具有较可靠的相位质量。
图5 优化前后差分干涉相位对比
3.2 融合方法效果分析
分别使用常规StaMPS-SBAS和联合DT和SDFPT的时序分析方法对覆盖实验区域的60景Sentinel-1A影像进行处理,得到了实验区域2018年4月24日—2020年4月25日期间的LOS向地表时序形变信息,实验区域年平均形变速率如图6所示。常规StaMPS-SBAS时序方法在实验区域共监测到135 955个目标点,选点密度约为591个/km
2
,大量的点位分布在村庄内的建筑物等附近,而联合DT和SDFPT的方法共监测到2 089 323个目标点,选点密度约为9 091个/km
2
,在非村庄区域也分布有大量的点目标,监测的点密度约为常规StaMPS-SBAS的15.4倍,大大增加了实验区域内点目标的空间分布盖度与密度。尤其值得注意的是,在矿区(图中黑色椭圆及方框所示)位置监测点的密度和覆盖度提升非常明显,这对矿区形变监测十分有利。这表明融合DT和SDFPT的时序InSAR方法在煤矿形变监测中具有非常好的效果。
图6 研究区LOS年平均形变速率
如图6 (a)和6 (b)所示,对比2种方法监测得到的实验区域形变分布结果可知,二者在空间上形变趋势的分布具有良好的一致性,形变量级也基本相仿。DT和SDFPT融合方法在低相干区域也能获取到充足的点目标。从所示形变速率图可以明显看出实验区域内有5个明显的形变区,分别属于寸草塔二矿、布尔台煤矿和霍洛湾煤矿的范围。形变最为严重的区域为寸草塔二矿,常规StaMPS-SBAS监测得到的最大形变速率为-163.2 mm/a,但是监测点非常稀少,不能准确反映矿区形变分布特征。联合DT和SDFPT的时序InSAR方法得到的最大形变速率为-172.74 mm/a,在形变区域内监测到足够数量且空间分布均匀的目标点,且获取到了更大形变量的目标点,提供了更加详实的煤矿区域地面形变信息,可以有效克服StaMPS-SBAS等常规时序监测方法在低相干区域(尤其是矿区)获取的点目标不足等问题。
3.3 形变结果检验
为了验证联合DT和SDFPT的时序方法监测结果的可靠性,本文通过对比分析2种方法得到的LOS向形变速率进行交叉验证,将常规StaMPS-SBAS方法监测得到的目标点作为参考值,选取与融合方法监测结果最邻近范围内的目标点为同名点(100 83对),根据同名点之间的形变速率值及其差异,得到了同名点形变速率相关性和差异分布直方图,如图7所示。如图7 (a)所示,采用线性拟合可得到
y
=0.956
x
-6.739的线性方程,Pearson相关系数为0.94,
R
2
为0.89,StaMPS-SBAS和联合DT和SDFPT的时序InSAR这2种方法监测得到的LOS向地面形变速率之间显示了较高的相关性。如图7 (b)形变速率差异直方图所示,2种时序监测方法得到的同名点地面形变速率之间的差异较小,标准差为4.056 mm/a。其中,差异区间为-10~10 mm/a的同名点占比最高,约为84.97%。形变速率差异绝对值在10~20 mm/a的同名点占比约为13.56%。
图7 基于SBAS-InSAR的形变精度验证
为了进一步验证融合方法的可靠性,本文采用与SBAS干涉策略不同的PS-InSAR对研究区域的数据进行处理,将得到的形变结果与融合方法进行对比。由于矿区的形变量级较大,而PS-InSAR方法未对时空基线进行限制,长基线干涉对失相干较为严重,且矿区形变相位梯度较大,无法有效获取沉降中心的形变结果,因此仅能选取2种方法沉降区边缘和非沉降区的同名点进行交叉验证。根据2组形变速率,得到其相关性和差异分布,结果如
图8
所示。同名点形变速率相关性分析得到
y
=0.795
x
-3.475的线性方程,Pearson相关系数为0.77,
R
2
为0.60,显示了较高的相关性。2种方法的同名点形变速率差异区间主要在-10~10 mm/a,标准差为4.286 mm/a。但由于PS-InSAR时空失相干更为严重,因此相关性低于SBAS与融合方法的相关性。综上所述表明,联合DT和SDFPT的时序InSAR方法获得的地表形变结果具有较好的可靠性。
图8 基于PS-InSAR的形变精度验证
3.4 实验区域地表形变分析
如图9所示,实验区域内共有5个明显的形变区域,从上到下依次命名A,B,C,D,E区。这5个形变区分别位于寸草塔二矿、布尔台煤矿和霍洛湾煤矿的范围之内。其余地区的形变量级较小,实验区域LOS向年均形变速率分布区间为-172.74~30.24 mm/a。在5个煤矿形变区中,区域A和区域C表现出更加严重的地面形变。区域A位于寸草塔二矿范围内,是形变速率最大的一个区域,中心形变速率约为-172.74 mm/a,形变范围面积约为2.81 km
2
。区域B,C,D均位于布尔台煤矿范围内,区域C是形变范围面积最大的一个区域,中心形变速率约为-170.84 mm/a,形变范围面积约为3.81 km
2
,区域B中心形变速率约为-158.56 mm/a,形变范围面积约为3.43 km
2
,区域D中心形变速率约为-146.29 mm/a,形变范围面积约为2.71 km
2
。区域E位于霍洛湾煤矿范围内,中心形变速率约为-126.45 mm /a,形变范围较小,约为0.79 km
2
。
为分析煤矿形变区的时序形变规律,在区域A,B,C,D,E内的形变中心和边缘分别选取a,b,c,d,e和a',b',c',d',e'共10个点作为特征点进行时序形变量分析,特征点的分布见
图9
。在形变结果中提取10个特征点2018年4月24日—2020年4月25日内的时序累积形变图,如
图10
所示。特征点的时序形变图在一定程度上反映了研究区这2 a内的形变趋势,在各个区域中心和附近选取的特征点的时序形变趋势基本一致。但由
图7
可知,5个形变区域的累积形变量大小和形变趋势都各不相同,在每个时段内累积形变量变化的速度也不相同。区域A在2018年4月—2019年4月间,矿区内地表形变量持续增大,地面形变变化速率较快,到2019年10月之前,地面形变变化有放缓的趋势,而后又继续保持,地表累积形变量为-307.96 mm。区域B和区域E的时序形变趋势大致相同,在监测时段内地表形变量变化较为平缓,但一直在持续发生。区域B的最大地表累积形变量为-186.76 mm,区域E特征点的最大地表累积形变量为-102.34 mm。区域C和区域D的时序形变趋势大致相同,但相较于其他区域较为复杂,从监测时段开始到2018年10月之前,2个区域的地表形变量变化缓慢,而后形变量变化速度持续增大直到2019年4月放缓。区域D直到监测时段结束形变趋势都较为平缓,地表累积形变量最大为-153.98 mm。不同于区域D,区域C从2019年10月—2020年4月地表形变量持续增大,地表累积形变量最大为-190.64 mm。
图10 特征点时间序列累积形变量
从特征点时序形变图可知,5个区域均有不同程度的显著形变,表明这些矿区在监测时间段内存在着煤矿资源的开采活动。从整体上分析这5个区域的时序地表形变规律,在春夏和冬季地面形变变化速率相对较快。该实验区域是大陆性季风气候,降雨量较少,主要集中在7—9月份,调查发现研究区2018年和2019年7—9月份均有不同程度的降雨,因此煤矿开采程度较其他时间段有所减弱,导致研究区地表形变量变化较为平缓。
3.5 典型矿区监测结果分析
为了进一步分析矿区的地面形变特征,选取形变区中形变速率较大的区域A、区域B和区域C作为典型矿区分析形变区的空间分布特征。分别沿形变区的形变方向在3个区域内绘制了3条剖面线,分别为F-F'(图11 (a))、G-G'(图12 (a))、H-H'(图13 (a)),并沿着剖面线提取地表累积形变量。为了更好地研究矿区内时序形变趋势,在监测时间段内以6个月左右为一期,共提取2018年10月2号、2019年4月7日、2019年10月4号和2020年4月25日4期形变剖面数据,分别如图11 (b)、图12 (b)和图13 (b)所示。从图11 (a)剖面图可知,区域A形变范围为椭圆形,沿西北-东南方向的地表形变影响距离较长,沿西南-东北方向的地表形变影响距离较短。结合图11 (b)剖面线形变图可知,区域A在2018年4月24日—2018年10月2日期间内的累积形变量级相较于其他时间段内较大,最大累积形变量达到了-179.68 mm,在2019年4月7日—2019年10月4日之间累积形变量较小。寸草塔二矿的形变剖面呈现明显漏斗状,说明该矿区只存在单个形变中心,且剖面线上的累积形变量随时间持续增长,剖面线上的形变量折线呈近似对称状,形变的严重程度自剖面线的两端向形变中心逐渐加剧。从图12 (a)剖面线示意图可知,区域B沿西南-东北方向的形变影响范围较小,沿西北-东南方向的形变影响距离较大,结合剖面线时序形变图分析,区域B在2019年4月7日—2019年10月4日期间的形变量较大,增幅约为95 mm,而后到监测时段结束时形变量的增幅都较小。剖面线形变量折线图只有一个形变中心,形变严重程度自两端向形变中心加剧,且形变中心离边缘的形变梯度非常大。从图13 (a)剖面线示意图可知,区域C地面形变涉及范围较大,沿南北方向形变影响距离较长,沿东西方向形变影响距离较小。结合图13 (b)可知,在2018年10月2日之前区域C的地表形变量较小,最大形变量为-35.5 mm。而在2019年10月4日—2020年4月25日期间地表形变量相对于其他时间段增加较多,增加量级约69.21 mm。剖面线上提取的形变量折线只有一个形变中心,且形变量随着时间持续增加,形变中心与边缘的差异较区域A更小,表明形变沿形变中心自东南向西北方向扩散。
图11 区域A剖面线时序形变图
图12 区域B剖面线时序形变图
图13 区域C剖面线时序形变图
4 结论
针对常规时序InSAR以及PS联合DT方法在矿区形变监测中的问题,本文尝试将DT和SDFPT进行融合。以Sentinel-1 SAR影像为数据基础,分别基于常规StaMPS-SBAS方法和联合DT和SDFPT的时序InSAR方法获取了布尔台矿区附近2018年4月—2020年4月时段内的LOS向地表形变结果,并分析了研究区的地表形变特征,结论如下:
1) 联合DT和SDFPT的时序InSAR方法在实验区域内监测得到的相干点密度相比于常规StaMPS-SBAS方法提高了15.4倍,点位分布均匀且在形变中心监测到更大形变速率的目标点,能够更好地反映实验区域的整体形变情况。
2) 对比分析常规SBAS-InSAR方法和常规PS-InSAR方法与融合方法得到的监测结果,2种方法均与融合方法具有较高相关性,Pearson相关系数分别为0.94和0.77,形变速率差异标准差分别为4.056 mm/a和4.286 mm/a,形变速率差异绝对值小于10 mm/a的同名点占比均大于84%。
3) 实验区域内存在5个明显的形变区域,形变影响范围面积最大的为3.8 km
2
