王 羽 王 磊 李靖宇 江克贵 李世保 杨 科 刘洪江
(1.安徽理工大学空间信息与测绘工程学院,安徽 淮南 232001;2.安徽理工大学矿山采动灾害空天地协同监测与预警安徽普通高校重点实验室,安徽 淮南 232001;3.安徽理工大学矿区环境与灾害协同监测煤炭行业工程研究中心,安徽 淮南 232001;4.中国矿业大学(北京)地球科学与测绘工程学院,北京 100083;5.安徽理工大学矿业工程学院,安徽 淮南 232001)
地下矿产资源开采会使采空区上方岩层受到重力影响而产生形变,易导致位于采空区上方的建(筑)物、民房、耕地以及交通设施等遭到破坏,还会引发地质灾害,影响了区域居民生产生活,并对矿区安全生产造成了很大威胁[1-5]。通常情况下,矿区的地表形变主要表现为地面沉降。传统的矿区地表形变观测技术主要包括地面观测技术和三维激光扫描技术。其中,地面观测技术主要在地面布设观测站,采用全站仪、水准仪、GPS-RTK[6-7]等测量仪器获取观测站各点的形变信息,再根据监测点上的形变来反演整个矿区的地表形变。该类方法虽然观测精度很高,但局限性较明显,即无法动态地反映整个区域的下沉状态,并且人力物力成本较高,不利于对一个区域进行长时间的追踪观测。三维激光扫描技术可以获得实时动态的地表形变信息,且观测精度可达到毫米级别,但是三维激光点云数据处理耗时费力、监测效率低、观测成本较高,使得该技术难以大范围推广应用[8]。受到多方面的条件限制,传统的矿区地表形变观测技术已经无法适应现代矿区灾害预警体系中对于地表形变全面且高效监测的要求。
20世纪60年代以来,合成孔径雷达技术(Synthetic Aperture Radar,SAR)的发展使得地表形变监测的研究迎来重大技术性突破,由于雷达波束的特殊性,使得该技术可不受气候条件和昼夜的影响,甚至可以透过部分地表或者植被获取地面信息[9-10]。该类特性使其在地震、冰川和地表形变监测等领域具有非常广阔的应用场景,尤其是近些年来迅速发展的差分干涉测量技术(Differential Interferometry Synthetic Aperture Radar,D-InSAR)能够以较低成本获取大区域较高精度的地表形变信息[11]。国外学者较早应用该技术对地表形变进行研究,并取得了较好的成果。1989年,GABRIEL等[12]应用D-InSAR技术成功测量了加利福尼亚河谷的地表形变,并且监测精度可达到厘米级,该项研究第一次证实了D-InSAR技术可实现地表形变的高精度监测。但是,D-InSAR技术在当时并未引起太大的反响,直至1993年,CLAUDIE等[13]使用雷达数据提取了美国加利福尼亚州兰德斯地震的位移场,并且通过与实测数据比较证明了D-InSAR监测结果比较准确,至此该技术开始受到学术界的重视。国内在D-InSAR形变观测方面的研究虽然起步比较晚,但也取得了很好的成绩。在2000年前我国学者主要是针对合成孔径雷达干涉测量(Interferometry Synthetic Aperture Radar,InSAR)技术的一些关键算法和应用进行了展望和试验研究,具体的技术应用案例比较少。2001年,刘国祥等[14]使用卫星雷达数据获取了香港近海回填地基1 a内的地表沉降信息,其结果精度优于1 cm。2011年,吴宏安等[15]使用ENVI ASAR影像数据获取了太原市2003—2009年的地表形变,论证了差分干涉技术可以应用于城市地表形变监测。2016年张永红等[16]提出了多主影像相干目标小基线干涉技术,并使用该技术提取了京津冀地区22 a间的地面沉降信息,并从时间和空间角度分析了地表形变的变化特征。
研究发现:对于矿区形变梯度比较小的区域,DInSAR能够有效获取其真实的地表形变量,但是当矿区形变梯度过大而超过了D-InSAR的监测能力时,会导致该技术无法获取矿区真实的地表形变信息。这是因为当形变梯度过大时,D-InSAR在进行相位解缠时会使整周相位丢失,从而使得最终解算的形变量远小于实际形变量。MASSONNET等[17]认为雷达数据分辨率、波长等都影响D-InSAR技术对地面沉降的监测能力,当同一位置的像元在某一时段内的形变量过大,且当该处像元在进行干涉测量时,由于相干性较低甚至出现失相干的现象,会导致该处像元无法直接通过D-InSAR技术获取真实准确的形变值。蒋弥等[18]研究发现,最大形变梯度和形变量不仅与相干性有关,多视的视数对可监测最大形变梯度和最大形变量的能力也有影响。冷红伟[19]对比分析了TerraSAR-X雷达数据和Sentinel-1A雷达数据分别在多视和全分辨率情况下监测大梯度形变的能力。
对于上述存在的问题,不少学者从不同角度提出了解决方法,其中在1999年,MICHEL等[20]提出了像素偏移量估计(Offset-tracking)技术,该技术利用SAR影像强度信息,通过影像配准原理进行互相关计算从而提取了像元间的偏移量。像素偏移量估计技术对于雷达影像的相干性要求很低,可避免相位解缠错误的问题,在一定程度上可实现大梯度地表形变监测,但是该技术监测精度较低,无法满足实际应用需求,并且在积水区域雷达的强度信息减弱,这些因素都导致了像素偏移量估计技术在地表形变监测中的应用受到制约。蒋创等[21]结合开采沉陷地表动态移动变形机理,研究了基于动态概率积分法模型(Dynamic Probability Integral Method,DPIM)约束的InSAR相位解缠模型和地表三维变形提取方法,克服了矿区大梯度形变区域解缠失败的不足,恢复了密集区缠绕相位,在较大程度上实现了矿区开采沉陷盆地三维形变提取。
针对D-InSAR因矿区形变梯度过大而无法获取真实地表形变量的问题,实现矿区全尺度梯度地表形变监测,本研究提出了一种融合水准和D-InSAR的矿区地表形变监测方法(LEV-InSAR模型)。该模型结合D-InSAR地表形变监测的基本原理,通过水准测量技术弥补D-InSAR相位解缠失败的不足,获取矿区完整连续且准确的地表形变信息。为了验证该方法的监测精度和适用性,针对山西省长治市霍尔辛赫矿区3501工作面,选取9景Sentinel-1A SAR影像,采用双轨D-InSAR技术获取2021年12月23日—2022年3月29日的矿区地表形变监测结果,并结合同时期观测线水准观测数据,成功解算了工作面内准确的地表形变量,取得了良好的监测效果。本研究提出的监测方法对于矿区大梯度形变监测具有一定的实用价值。
1.1 双轨D-InSAR技术基本原理与流程
双轨D-InSAR技术基本原理是利用同一区域不同时相的两景雷达影像,依据雷达影像中包含的相位信息进行干涉处理,得到干涉相位图。这些干涉相位中包含多种相位信息(式(1)),通过各种方法解算出形变相位,再将形变相位按照一定的关系式转化成形变量。
式中,Φtopo为地形相位;Φflat为平地效应相位;Φdef为地表形变相位;Φatm为大气延迟相位;Φnoise为噪声相位。
得到真实地表形变相位Φdef后,进而计算雷达沿视线方向的地表形变ΔR,公式为
式中,λ为雷达波长,mm。
双轨D-InSAR基本流程如图1所示。
图1 双轨D-InSAR技术流程Fig.1 Flow of dual-track D-InSAR technique
1.2 LEV-InSAR方法基本原理与流程
在实际应用中,当影像相邻像元的形变梯度过大时,D-InSAR的整周模糊数相位解缠结果不正确,导致丢失若干整周相位(解算的相位差与真实相位差之间相差2π的整数倍),从而使最终的D-InSAR监测形变结果与真实量之间相差若干个整周相位所对应的形变量。理论上,通过D-InSAR得到的地表形变监测值中包含的不足整周部分相位对应的形变量是较为准确的。
针对D-InSAR监测方法的不足,提出了一种融合水准数据的D-InSAR矿区全尺度梯度地表形变监测方法(LEV-InSAR模型)。该模型从观测线水准数据出发,借鉴概率积分法思想,应用观测线拟合曲线方程得到工作面上方地表任意点的模拟形变值,由于模拟形变值与地表真实形变值的量级相当,实际上,根据模拟形变值可以提取出D-InSAR监测像元的实际整周模糊数,最后整合解算的整周相位形变量和去除整周相位后的D-InSAR监测值,即可获得开采工作面全尺度梯度地表形变监测信息。具体思路是:首先在矿区工作面上方地表以走向线和倾向线为基准构建XY工作面坐标系,并将水准监测点所在的大地坐标系转换到工作面坐标系内;然后将水准测量值沿工作面走向和倾向方向拟合得到两个拟合曲线方程,即倾向线拟合曲线方程和走向线拟合曲线方程,记为W(x)和W(y)。
根据概率积分法理论,对于工作面内任意一点的形变值可表示为
式中,k为系数,通过最小二乘原理拟合得到。
通过式(3)估算得到的任一点形变值W(x,y)并未进行相位缠绕,则可以通过该值计算出工作面任意一点真正的整周模糊数,再由整周模糊数计算得到形变量WZ(x,y):
式中, 一个整周相位相当于1/2的雷达波长量级的形变;int·()为取整函数。
通过D-InSAR获得的地表任意一点的形变值记为WD(x,y),当矿区形变梯度过大时,由于影像失相干导致在相位解缠时丢失了若干整周相位,使得形变值WD(x,y)与真实值之间相差若干个整周相位量级的形变量,然而形变值WD(x,y)中包含的不足整周部分相位对应的形变量较为准确。一般来说,WD(x,y)中也可能含有部分整周相位形变量,为了不与WZ(x,y)中部分形变量重复,应从WD(x,y)中去除整周相位形变量,只保留不足整周部分相位对应的形变值,记为WD′(x,y)。
其中,n为非负整数;WDi(x,y)为任意一点处第i期的D-InSAR监测值,mm;m为D-InSAR监测总期数;WDi′(x,y)为经过去除整周相位后的任意一点第i期D-InSAR监测值,mm;WD′(x,y)为经过去除整周相位后任意一点的D-InSAR累计监测值,mm;“±”号由WDi(x,y)值正负决定,当WDi(x,y)为正时,取“-”号,当WDi(x,y)为负时,取“+”号。
将通过水准数据模拟得到的整周相位形变值WZ(x,y)与从D-InSAR监测值中提取的不足整周部分相位对应的形变值WD′(x,y)求和,即可得到工作面内任意位置的地表形变值W′(x,y):
LEV-InSAR模型的具体流程如图2所示。
图2 LEV-InSAR模型的具体流程Fig.2 Detailed flow of LEV-InSAR model
2.1 研究区概况
选取山西省长治市长子县霍尔辛赫矿区3501工作面作为研究区域。该工作面位于长子县城区东北方向,距离长子县主城区约3.4 km,地理坐标为东经112°54′30″～112°55′10″,北纬36°7′40″～36°8′10″。工作面东西宽约300 m,南北长约800 m,占地面积约0.24 km2。3501工作面开采煤层平均厚度为5.6 m,煤层倾角为1°～12°,平均开采深度为576 m,工作面日推进速度为3～4刀,约为2.4 m/d。3501工作面于2022年2月初开始回采,2022年2月8日开始离层注浆,2022年6月15日该工作面基本停采,2022年9月15日注浆停止。工作面内地形起伏较小,地表高程为935～939 m。该工作面内EW向县道横穿工作面,道路南侧有一焦化厂,其余部分为大面积的农田。工作面内共设置了两条观测线,一条为EW向的倾向观测线,点号标记为Q01～Q22,另一条为SN向的走向观测线,点号标记为Z01～Z27,共49个观测站。工作面位置及观测线布置如图3所示。
图3 3501工作面位置及观测线布置Fig.3 3501 working face position and layout of observation lines
2.2 数据源
2.2.1 Sentinel-1A SAR影像数据
根据3501工作面的开采时间和进度,并且为了避免工作面上方夏秋季节农田植被对D-InSAR干涉相位造成较大影响,选定2021年12月—2022年3月为研究时段。在此基础上选取了覆盖研究区域的9景C波段Sentinel-1A 升轨数据,数据类型为单视复数影像(Single Look Complex,SLC),成像模式为干涉宽幅(Interferometric Wide swath,IW),极化方式为VV(垂直传送和垂直接收),成像日期为2021年12月23日—2022年3月29日,影像的具体参数取值见表1,数据覆盖范围如图4所示。DEM数据来源于美国航空航天局(National Aeronautics and Space Administration,NASA)和美国国家地理空间情报局(National Geospatial-Intelligence Agency,NGA)免费对外开放的航天飞机雷达地形测绘任务(Shuttle Radar Topography Mission,SRTM)数据,SRTM DEM数据的空间分辨率为30 m,绝对高程精度为±16 m,DEM数据将用于去除差分干涉测量过程中由地形起伏因素引起的地形相位。
图4 Sentinel-1A SAR影像覆盖范围Fig.4 Coverage of Sentinel-1A SAR image
表1 Sentinel-1A SAR影像参数Table 1 Sentinel-1A SAR image parameters
2.2.2 水准数据
本研究依据9景Sentinel-1A SAR影像成像时间,选取该工作面2021年12月22日—2022年3月29日走向线和倾向线上地表移动观测站的观测数据。该观测数据来源于17期三等水准测量的实测数据,第1期观测时间为2021年12月22日,第17期观测时间为2022年3月28日,平均每5 d观测1次,测量仪器为莱卡sprinter100数字水准仪,标称精度为每千米往返测高程中误差0.4 mm,采用平差软件进行水准网平差,平差后最弱点中误差为2.93 mm。原定观测方案中共布设了49个观测站,但在后期观测过程中,部分观测站因外部因素影响而遭到破坏。为保证水准数据的可靠性和有效性,本研究在选取观测站数据时,只保留了部分观测站水准数据,其中包括20个倾向线观测站和26个走向线观测站,共计46个观测站的水准观测数据。
2.3 结果与分析
2.3.1 D-InSAR地表形变监测结果分析
对获取的9景Sentinel-1A数据进行雷达差分干涉测量处理,按照时间序列构建8个干涉对,每个干涉对时间基线均为12 d。按照D-InSAR实用化数据处理流程对8个干涉对进行处理,得到相邻影像之间的形变信息,通过统计矿区各时间段的沉降信息,得到该工作面2021年12月23日—2022年3月29日的累计形变图以及观测线下沉曲线如图5所示。
图5 D-InSAR地表形变监测结果Fig.5 Monitoring results of ground deformation by D-InSAR
由图5可知:该矿区在倾向线以北发生的形变较为明显,沉降中心位于走向观测点Z06附近,沉降分布范围及形态类似于一个漏斗,在沉降中心区域沉降量最大,边缘部分沉降量较小。2021年12月23日—2022年3月29日,D-InSAR所监测的最大累计沉降量为172.5 mm,与水准实测数据相比,D-InSAR地表形变监测结果的均方根误差为112.0 mm。在沉降漏斗边缘及外部区域,D-InSAR监测到的累计沉降量小于20 mm。这些微量形变形成原因为:① 受到矿区开采影响而造成的微小形变;② 可能是卫星接收雷达信号时产生的误差以及卫星本身固有的误差;③ 数据处理过程中产生的误差。
为了从数值上反映D-InSAR地表形变监测效果,按照观测点位置提取了若干个观测点的D-InSAR地表形变监测值,并与水准实测数据进行了对比,结果见表2。
表2 D-InSAR地表形变监测误差统计Table 2 Statistics of D-InSAR surface deformation monitoring errors
由表2可知:在矿区沉降盆地边缘,地表形变量较小,D-InSAR监测结果与水准监测结果吻合程度较高,但随着沉降程度的增加,在沉降严重的水准点上,D-InSAR监测结果与水准监测结果之间的绝对差值明显增加。这一现象很好地说明了在地面形变梯度过大区域,差分干涉测量的相位解缠错误使得D-InSAR监测结果与真实地表形变信息之间有较大差异,并且这种差异在沉降盆地中心表现得更为明显。差分干涉测量技术虽然在矿区形变监测中具有显著优势,但从监测结果来看,其更适用于矿区周边缓慢微小形变监测,对于形变梯度大的沉降中心区域监测能力明显不足。
2.3.2 LEV-InSAR模型监测结果分析
基于本研究提出的融合D-InSAR和水准的地表形变监测方法获取3501工作面地表形变信息,按照LEV-InSAR模型的基本思路,首先需要构建LEVInSAR模型监测方程。根据观测站点位信息,构建基于剖面观测线的XOY标系(图6),赋予各观测站点横(x)纵(y)坐标。
图6 XOY坐标系示意Fig.6 Schematic of XOY of coordinate system
按照观测线水准观测数据,分别拟合得到走向和倾向观测线拟合曲线方程:
式中,x0,xc,ω1,A1,y0,yc,ω2,A2均为刻画观测站点坐标与其形变值之间函数关系的参数,无实际含义,当坐标信息和形变值确定时,这些参数均为实际值。
根据MATLAB非线性拟合得到的式(7)中各参数取值分别为:x0=-4.860 72,xc=543.543 81,ω1=160.062 82,A1=-136 062.176 7;y0=-32.115 81,yc=321.732 65,ω2=110.944 8,A2=-117.203 99。
由式(3)可得工作面内任意一点处的形变值估算方程为
式中,k取-230。
根据式(5),分别对每一期D-InSAR监测结果进行甄别,经过去除整周相位形变量后累加,得到WD′(x,y),再根据式(4)和式(6)得到本研究工程实例的LEV-InSAR模型监测方程为
其中,波长λ取56 mm。
按照式(9)模型,融合水准观测数据和D-InSAR形变量解算得到3501工作面的LEV-InSAR模型地表形变监测结果,如图7所示。
图7 LEV-InSAR模型地表形变监测结果Fig.7 Surface deformation monitoring results by LEV-InSAR model
为验证该方法的监测精度,分别提取了两条观测线的下沉曲线,如图8所示。此外,统计了各观测点的LEV-InSAR模型地表形变监测值相对于水准观测数据的绝对误差和相对误差,见表3。
图8 主断面下沉曲线Fig.8 Subsidence curves of main section
表3 LEV-InSAR模型监测误差统计Table 3 Statistics of monitoring errors LEV-InSAR model
由图8、表3可知:倾向线和走向线上观测点的监测结果与实际测量的下沉值基本吻合,倾向线和走向线的下沉曲线和水准曲线的形态特征趋近。与水准数据相比,倾向线上各观测点误差为-16.9～16.7 mm,最大监测误差为16.9 mm,约为最大沉降值的3.3%,最小监测误差为0.8 mm,倾向线均方根误差为10.2 mm,平均相对误差为14.4%;走向线上各观测点误差为-32.0～41.0 mm,且大部分观测点的监测误差小于20 mm,最大监测误差为41.0 mm,约为最大沉降值的8%,最小监测误差为0.7 mm,走向线均方根误差为15.8 mm,平均相对误差为26.6%。从单点误差来看,由于个别点实际下沉值较小,使得这些观测点的相对误差对于绝对误差比较敏感,使得该类观测点的相对误差和整体的平均相对误差较大。对于累计沉降量大的观测点而言,相对误差较小,由此可见,LEV-InSAR模型对于矿区大梯度地表形变监测具有更明显的优势。
针对D-InSAR无法实现矿区全尺度梯度地表形变监测的不足,提出了一种融合D-InSAR和水准的矿区全尺度梯度地表形变监测方法(LEV-InSAR模型)。研究表明:LEV-InSAR模型对于矿区大梯度地表形变监测具有较好的适用性和可靠性,可以提取到矿区空间连续的真实地表形变信息,具有较好的应用价值。
猜你喜欢观测线水准梯度一个带重启步的改进PRP型谱共轭梯度法数学物理学报(2022年1期)2022-03-16一个改进的WYL型三项共轭梯度法数学物理学报(2021年6期)2021-12-21一种改进的水准网条件平差算法北京测绘(2021年1期)2021-02-05煤矿开采地表移动变形观测与分析陕西煤炭(2020年6期)2020-11-23一种自适应Dai-Liao共轭梯度法应用数学(2020年2期)2020-06-24采煤工作面地表岩移观测方法设计陕西煤炭(2019年5期)2019-09-06亿欣煤业15#煤综采工作面地表岩移观测站的设计山西煤炭(2019年2期)2019-08-29一个具梯度项的p-Laplace 方程弱解的存在性华东师范大学学报（自然科学版）(2019年3期)2019-06-24媲美激光光源的成像水准Acer宏碁E8620C家庭影院技术(2018年11期)2019-01-21高速铁路轨道控制网(CPⅢ)高程网建立方法探讨铁道勘察(2011年1期)2011-06-07