基于SBAS-InSAR与GM(1,1)模型的张博线采空区地表形变监测

刘瑞祥, 陶秋香*, 刘晓朋, 韩宇, 李雪鹏, 肖怡欣

(1.山东科技大学测绘与空间信息学院, 青岛 266590; 2.中铁工程设计咨询集团有限公司济南设计院, 济南 250022; 3.山东省青岛市勘查测绘研究院, 青岛 266033)

当今经济社会的快速发展需要大量的煤炭资源作为支撑,但在开发地下资源的过程中,形成了大量地下采空区,使原有的土质结构遭到破坏,很容易造成地面塌陷等人为灾害[1]。而拟建的张博线铁路沿线(即淄博市张店区至博山区路段)恰好经过了大量采空区,若采空区上覆地表尚未稳定,存在较严重的地表形变,就会对铁路建设以及运营造成安全隐患,给人们的出行带来危险。

尽管水准测量、全球导航卫星系统(global navigation satellite system, GNSS)测量等技术手段都能够进行精确的地表形变监测,但却具有一定的局限性,例如,人力、物力投入较多,难以实现大规模的覆盖和普查[2]。近年来,随着空间对地观测技术的发展,合成孔径雷达差分干涉测量(differential interferometric synthetic aperture radar, D-InSAR)技术[3-4]和小基线集合成孔径雷达差分干涉(small baseline subset interferometric synthetic aperture radar, SBAS-InSAR)技术[5]的发展,以其高效率、低成本、大范围的优点,极大地弥补了传统测量的不足,为地表形变监测提供了新的测量方法,之后的研究员利用这些技术取得了一系列的研究成果[6-11]。GM(1,1)灰度预测模型被提出后,为研究人员对地表形变的预测提供了一种极为简便可行的方法,石晓宇等[12]、刘哲[13]、郭在洁等[14]分别对GM(1,1)在不同情形下的应用进行了研究。

现利用SBAS-InSAR技术对43景Sentinel-1A SAR影像数据进行处理,获取拟建的张博线铁路沿线的地表形变信息;以张博线铁路沿线为中心线,对其左右两侧1 km缓冲区内总体地表形变情况、重点形变区域及其里程点的地表形变情况等进行分析;利用GM(1,1)灰度预测模型对地表特征点监测结果进行模拟和预测,以研究GM(1,1)在采空区地表形变监测中预测的可行性。

1 研究区与数据源

1.1 研究区概况

张博线位于山东省淄博市境内,北起胶济线淄博站,沿张南路南行至南定镇,转向西南,跨越省道S703及国道G205后沿G205至淄川区,向西跨越孝妇河,继续沿G205直至设计终点博山站。线路经过区域有淄博市张店区、淄川区及博山区,线路全长39.233 km。沿线村镇密集。图1为张博线地理位置及其经过的主要采空区。

图1 研究区地理位置Fig.1 Geographical location of study area

1.2 数据源

Sentinel-1SAR卫星为载有C波段的合成孔径雷达,具备干涉宽带(interferometric wide swath,IW)等多种工作模式,是欧洲航天局哥白尼计划中的极其重要的地球观测卫星。该卫星具有水平发射水平接收(VV)、水平发射垂直接收(VH)等不同的极化方式,能提供连续的、全天候的SAR影像,有着高重访频率、高覆盖能力等无可比拟的的优点,以及极好的时效性和可靠性[15]。现选取覆盖张博线,时间跨度为2019年1月5日—2022年1月1日的43景C波段、VV极化的哨兵1A升轨数据。表1为43景SAR影像的获取日期。

表1 Sentinel-1A SAR影像数据情况Table 1 The data information of Sentinel-1A SAR images

2 原理与方法

2.1 SBAS-InSAR技术原理

Berardino等[16]在2002年提出的SBAS-InSAR技术克服了大量传统InSAR的诸多弊端,有效地提高了地表形变监测的精度。

设研究区域有N+1幅雷达影像,它们获取的时间分别为t0、t1、…、tn,现在选择其中一幅影像作为公共主影像,任何影像都能与其他至少一幅影像进行干涉,可以得到M幅干涉图,M满足如下关系式。

(1)

利用在tA和tB时刻获得的两景已经去除地形相位的SAR影像(tA

δφi(x,r)=φ(tB,r,x)-φ(tA,r,x)

(2)

式(2)中:φ(tA,r,x)和φ(tB,r,x)分别为像元(r,x)在tA与tB时刻相对于参考时刻t0的雷达视线向的累计形变相位;λ为雷达的波长;d(tA,r,x)和d(tB,r,x)分别为tA与tB时刻相对于t0的雷达视线向的累计形变量,d(t0,r,x)≡0。有

(3)

假设IE=[ IE1IE2… IEM]、IS=[ IS1IS2… ISM]分别为干涉数据处理时按照时间顺序排列的主、辅影像序列,并满足条件为

IEj>ISj,j=1,2,…,M

(4)

那么所有差分干涉图相位就可以组成的观测方程为

δφi=φ(tIEi)-φ(tISi)

(5)

将式(5)转换为线性方程可得

Aφ=δφ

(6)

式(6)中:A为M×N维矩阵,当所有的SAR影像都被分在一组,且M≥N,矩阵A的秩是N,通过对式(6)采用最小二乘法求解,可得

(7)

2.2 GM(1,1)模型

GM(1,1)模型是一种用于研究数据量与信息较少的不确定性问题的理论方法,是邓聚龙教授于1982年提出的一种灰色系统理论,利用该理论可以利用较少的数据来获得较好的预测效果[17]。原理如下。

设有原始观测数列为

X(0)=[X(0)(1),X(0)(2),…,X(0)(n)]

(8)

将原始观测数列进行一次累加,得到新的数列为

X(1)=[X(1)(1),X(1)(2),…,X(1)(n)]

(9)

对数列X(1)建立一阶微分方程可得

(10)

式(10)中:a为发展灰数,反映X(1)及X(0)的发展趋势;u为内生控制灰数,反映了数据间变化的关系。

k=2,3,…,n

(11)

令

Z(1)(k)=μX(1)(k)+(1-μ)X(1)(k-1),

k=2,3,…,n

(12)

式(12)中:Z(1)(k)被称为式(10)的背景值;μ∈[0,1],μ为权重系数。

假定μ取值为0.5,则有

(13)

将式(10)离散化,则有

X(1)(k)-X(1)(k-1)+aZ(1)(k)=u,

k=2,3,…,n

(14)

利用最小二乘法求解式(14)可得

(15)

求得a和u后,继续求解式(10),得

(16)

将式(16)离散化得

(17)

(18)

代入式(17)得

k=0,1,…,n-1

(19)

对式(19)做累减还原便可得到X(0)的灰色预测模型为

k=2,3,…,n-1

(20)

2.3 数据处理流程

数据处理采用的平台是世界著名的瑞士 GAMMA 遥感公司开发的专门用于处理SAR数据的商业软件:GAMMA,采用的方法是SBAS-InSAR技术,充分利用现有43景sentinel-1A数据,得到最优结果。首先,利用连续的常规差分方案,通过设定的时空基线指标使43景数据进行自由配对组合的方式,并选取其中2020年1月24日获取的SAR影像作为公共主影像,然后从剩下的42个单主影像组成的干涉对中选取时空基线指标较好的像对,联合外部参照 DEM 处理为雷达坐标系下的模拟解缠地形相位,之后再进行常规差分处理, 获取区域内形变概略位置,以及形变区的演变过程,据此提取强度阈值、相位相干性阈值等经验参数;然后进行SBAS-InSAR 处理,通过时空基线指标进行自由配对组合,以时间基线小于96 d,空间基线小于 194.9 m 为时空临界基线,选取146个组合对,进行差分干涉处理,得到差分干涉图;对差分干涉图进行滤波处理,去除大气延迟误差、轨道误差、残余DEM误差等;对滤波后的差分干涉图进行相位解缠;利用解缠后的干涉图进行精密基线估计;然后利用估计的精密基线再次进行“差分干涉处理→滤波处理→差分干涉图解缠”处理;对第二次获取的解缠后的差分干涉图进行浏览,剔除依然存在轨道误差、大气延迟误差的差分干涉图,形成最优组合的解缠后的差分干涉图序列;最后进行地表形变反演和归一化处理,得到研究区的地表形变时间序列;对求取的形变结果进行克吕格插值方法分析,基于参考点进行整体修正,得到最终形变结果。具体的数据处理流程如图2所示。

图2 数据处理流程Fig.2 Data processing flow

3 实验结果分析

3.1 铁路沿线沉降分析

利用图2所示的数据处理流程和方法对表1的43景Sentine-1A SAR数据进行处理,获取2019年1月—2022年1月张博线铁路沿线的地表形变信息,对张博线铁路沿线左右500 m缓冲区内的地表形变进行分析,得到如图3所示的铁路沿线形变速率图和如图4所示的1 km缓冲区内最终累计形变图。

图3 铁路沿线1 km缓冲区内形变速率Fig.3 The deformation rate in the 1 km buffer zone along the railway

分析图3和图4可以看出,整个铁路沿线在监测周期中的地表形变以沉降为主,最大沉降速率与累计沉降量分别为-7.9 mm/a与-23.8 mm;小部分区域存在微小抬升,最大抬升速率与累计抬升量为2.7 mm/a与7.7 mm。其大部分区域年平均形变速率在-2.0～2.0 mm/a,累积形变量集中在 -5.0～5.0 mm,整个区域相对稳定;结合矿区范围来看,各矿区域范围内都存在沉降区,且沉降范围与沉降量各不相同。除此之外,在不存在采空区的十里铺村西南侧存在范围较大的持续沉降区域,应进行实地调查以查明原因。

为了更好地分析铁路沿线的形变区域进行分析,将形变范围分成了3个等级,大于等于0 mm、小于等于0 mm、小于等于-5 mm共3段,如图5所示。

图5 监测周期内3个区间内累计地表形变区域分布Fig.5 The regional distribution of cumulative surface deformation in 3 intervals during the monitoring period

分析图5,存在沉降的区域主要位于GDK1～GDK11、GDK16～GDK21、GDK28～GDK30、GDK31～GDK34、GDK37～GDK39,沉降区域的面积为 18.95 km2;沉降较严重(沉降量多于-5.0 mm)的区域主要位于GDK16～GDK19区域,在 GDK23～GDK24、GDK8～GDK9附近也有少许,通过进行面积统计发现,沉降小于等于-5.0 mm的区域的面积为3.88 km2,小于整个监测区域的1/10,整体相对稳定;GDK3～GDK5、GDK11～GDK16、GDK21～GDK24、GDK25～GDK28、GDK34～GDK37等附近区域未沉降或略有抬升的区域面积为21.22 km2。

分别提取铁路中线以及左右两侧各500 m在监测周期内的累计形变量,制作剖面图如图6所示。

图6 铁路沿线中心及其两侧各500 m沿线形变速率剖面图Fig.6 The deformation rate profile along the center of the railway and 500 m on both sides

分析图6,在2019年1月—2022年1月期间,铁路沿线的形变速率整体都较小,最大不超过 -4.0 mm/a,从中线可以看出,形变速率主要集中在-1.0～1.0 mm/a,其中在16～18 km,沉降速率超过了-1.0 mm/a,最大达到-1.7 mm/a;从形变稳定性上来看,中线整体波动较小,但是在11～23 km 处,形变出现了较大波动,形变速率由 -0.2 mm/a 变化到-1.7 mm/a,再由-1.7 mm/a变化到-0.2 mm/a,变化幅度相对于其他路段较大;对于左侧500 m,相对于中线来说,铁路沿线的形变速率和形变波动幅度,都有了一定程度的增大,在6～9 km处出现了最大形变速率,为-2.6 mm/a,在16～18 km,沉降速率均超过了-2.0 mm/a;在形变稳定性上,在6～10 km和23～26 km处,形变呈现出W形起伏,在11～23 km处,形变呈现出V形起伏,且起伏程度相对较大,但考虑到该区域存在河道,应进行实地勘察,排除不稳定沉降对铁路线的威胁;对于右侧500 m,从整体来看,大部分地区形变速率相对较为稳定,形变速率在-2～1.3 mm/a, 在16.4 km处出现了最大沉降速率,为-4.0 mm/a;在15～21 km处(光正实业西北部河道内),形变呈现出W形起伏,且起伏状态较大,稳定较差,须进行实地排查。

3.2 GM(1,1)模型在地表形变监测中的应用

选择沉降较为严重的区域选择P1、P2、P3三点,以及出现抬升的P4点,如图7所示。提取此4个点在43个时期内的累计沉降值,以前40期数据为初始数据建立GM(1,1)模型,并用该模型对后三期进行预测,与实际数据相比较。

图7 选取测试点的位置Fig.7 The location of the selected test point

预测结果如图8～图11所示。其中,P1点的发展灰数a=-0.037 8;P2点的发展灰数a=-0.046 9;P3点的发展灰数a=-0.050 1;P4点的发展灰数a=-0.036 2。

图8 P1点GM(1,1)预测模型与监测值对比Fig.8 The GM (1,1) prediction model of P1 point is compared with the monitoring value

图9 P2点GM(1,1)预测模型与监测值对比Fig.9 The GM (1,1) prediction model of P2 point is compared with the monitoring value

图10 P3点GM(1,1)预测模型与监测值对比Fig.10 The GM (1,1) prediction model of P3 point is compared with the monitoring value

由文献[18]可知:

(1)当-a<0.3时,GM(1,1)可用于中长期预测。

(2)当0.3<-a<0.5时,GM(1,1)可用于短期预测,中长期预测慎用。

(3)当0.5<-a<0.8时,用GM(1,1)作短期预测应十分谨慎。

(4)当0.8 <-a<1时,应采用参差修正GM(1,1)模型。

(5)当-a>1时,不宜采用GM(1,1)模型。

即GM(1,1)模型在本次地表形变监测中可以用于中长期预测。

4 结论

利用43景Sentinel-1A数据,基于SBAS-InSAR技术获取了2019年1月5日—2022年1月1日淄博市拟建张博线铁路沿线地表形变信息,对其分析可得出以下主要结论。

(1)2019年1月—2022年1月期间,沉降区域的面积在逐年增加,统计发现三年来最终呈现沉降的区域面积为22.83 km2,占总面积的51.8%,但沉降量级不大,大部分区域年平均形变速率在-2.0～2.0 mm/a,整个区域比较稳定;沉降最严重的区域位于GDK23+550左侧232.2 m,出现了一个明显的沉降中心,最大沉降速率达到-7.9 mm/a, 最大沉降量达到-23.8 mm,且呈现出线性沉降趋势,对其沉降情况需要进一步关注。

(2)大部分里程点的形变较为稳定,形变速率在-1.0 ～ 1.0 mm/a,累计形变量在-3.0～3.0 mm。形变严重的点为GDK16、GDK17、GDK18,沉降值分别达到-9.0、-11.0、-11.0 mm,需要进一步关注。

(3)铁路中线大部分地区处于稳定状态,铁路左侧500 m和右侧500 m沿线形变速率和波动幅度要略大于铁路中线,需要进一步关注,为铁路的安全建设与运行提供保障。

(4)GM(1,1)模型在地表形变监测中拟合度较高,表现较好,可用于中长期的预测。