1. 引言
辽河平原地域辽阔,面积约3万多平方公里,长期的油气资源、矿产资源(煤炭)、地下水和地热资源的持续开采,使得辽河平原区域内地质环境问题逐步显现,地面沉降等地质灾害现象呈加剧发展趋势。地面沉降由于具备影响范围大,持续时间长,形成缓慢,不易察觉的特性,一旦形成,就会对地质环境系统造成不可逆的损害,破坏地面稳定性,造成建筑物倾斜、开裂等,严重威胁城市健康安全。传统沉降监测手段很难满足辽河平原监测需求 [1] 。合成孔径雷达干涉测量(InSAR)技术逐渐成熟并得到了工程化应用,已经成为地表形变监测的主要技术手段,在全球及区域性地形测图、大尺度地表形变监测中得到广泛应用,如地震前后的地表形变监测,由于地下资源过度开采等因素造成的城市地面沉降。因此本项目利用SBAS-InSAR技术,采用两种SAR数据源,对辽河平原整体地区进行沉降监测 [2] [3] [4] [5] 。
2. SBAS-InSAR技术原理
SBAS-InSAR方法可以通过限制空间基线和时间基线对满足条件的SAR影像进行干涉组合,能够有效降低因时间间隔较长而引起的差分干涉图去相干现象的影响。该方法除了将高质量的点目标作为地表形变信息的表达载体外,还能够充分利用SAR影像场景中具有一定相干性的分布式目标参与解译,进一步提高离散点位的监测密度,更加完整得反映形变区域的空间形态和演化规律。假设存在N + 1幅经配准的时序SAR影像,经基线组合选取后得到M个差分干涉对,且M满足条件
,第
个差分干涉对的第
个像素的相位模型为 [6] [7] [8] :
(1)
其中,
为雷达视向上的形变相位;
为DEM误差造成的残余地形相位;
为大气影响相位;
为轨道信息不精确引入的相位;
为其它噪声相位。对每个干涉图的模型相位进行相位解缠,通过奇异值分解(SVD)便能获取到地表形变信息 [9] [10] [11] 。
3. 数据处理
3.1. SAR数据源选择
由于辽河平原地区面积大,且整个区域分为重点监测区域与非重点监测区域,因此为节约成本并同时满足监测需求,本项目采用RADARSAT-2及Sentinel-1两种SAR数据,其中RADARSAT-2属于商业采购部分,应用于重点监测区域,其他区域采用Sentinel-1数据,其是免费的。
3.2. RADARSAT-2数据简介及处理方法
3.2.1. RADARSAT-2数据简介
RADARSAT-2卫星超宽精细(Extra Fine, XF)模式影像空间分辨率为5 m,幅宽125 km,采用C波段成像(波长5.6 cm),兼具高分辨率、大幅宽以及复杂地物环境下的地面沉降监测能力。根据项目监测需求查询监测区的卫星影像覆盖情况,综合考虑研究区大小、数据空间覆盖、时间序列完整性、空间基线等因素,本次处理选择2021年10月影像至2022年10月期间两组共16景SAR数据。
3.2.2. RADARSAT-2数据处理
针对RADARSAT-2数据,采用GAMMA结合StaMPS软件进行处理。数据预处理部分由GAMMA处理,时序分析部分由StaMPS处理 [12] [13] [14] [15] 。处理流程见图1。
Figure 1. RADARSAT-2 data processing flow chart
图1. RADARSAT-2数据处理流程图
成像处理:导入RADARSAT-2数据,将原始影像转换成软件可以处理的文件类型,获取每景影像对应的slc及slc.par参数文件。
影像配准及外部DEM配准:根据时间基线、空间基线选取配准主影像,将其他影像与主影像进行配准。SAR影像配准精度要求小于1/8像元。同时,将主影像与外部DEM进行配准,获取雷达坐标系下高程文件,为地理编码和差分干涉图做准备。
干涉对组合与差分干涉图生成:设置一定时间与空间基线阈值,将配准后的SLC影像进行干涉对组合,然后进行差分干涉处理,获取短基线集差分干涉图。
高相干目标点提取:将生成的差分干涉图转换格式,导入StaMPS中进行后续处理。迭代计算每个差分干涉图中每个候选点的噪声相位,结合SLC影像幅度离差对高相干目标点进行初选。计算初选点的噪声标准差,剔除噪点及质量不好的点,获取高相干目标点。
相位解缠:利用三维相位解缠方法对高相干目标点进行相位解缠。初步获取形变相位,此时形变相位包含轨道误差、DEM残余误差和大气相位等,需要分步进行剔除。
轨道、大气相位去除:对解缠后的形变相位进行最小二乘估计,估算DEM残余误差和主影像大气分量,利用曲面拟合方法估计轨道残余相位,然后通过时空滤波方法去除大气影响,最终获取地表形变信息。
3.3. Sentinel-1数据处理
Figure 2. Sentinel-1 data processing flow chart
图2. Sentinel-1数据处理流程图
数据导入与连接图生成:将原始卫星影像和精密轨道星历数据导入进行TOPSAR处理。生成SLC数据以及索引文件、SLC影像强度图、谷歌地球坐标文件、SARscape的索引识别文件、图形矢量文件。选择3D解缠模式,平台会根据导入的数据的相对关系自动选取最佳的超级主影像,并根据影像数据选取合适的时间基线和空间基线阈值进行配对,其他辅影像会以主影像作为参考影像,以此来生成辅影像与主影像的时间基线、空间基线信息,生成SAR的数据连接图。流程如图2所示。
干涉工作流及连接图编辑:因为使用的是哨兵数据,所以多视比设为4:1。SBAS-InSAR干涉流参数选取很大程度上决定了干涉结果的好坏,设置参数要根据影像数量及测区地表环境相关。利用覆盖研究区的DEM数据对影像进行辅助配准。本项目测区地质条件复杂,相干性较低,因此选择Goldstein算法进行干涉图滤波。由于研究区域内含有多处积水区,同时草木植被覆盖丰富,因此选取Delaunay三角网方法进行解缠。为了提高相位解缠质量,大幅降低相位解缠期间的错误率,可将分解等级调整为2。
轨道精炼和重去平:这步工作主要是去除上一步成果中的轨道、相位坡道等误差。这一步最重要的步骤就是选取地面控制点(GCP)。选取GCP点的依据:要远离形变区,在相位平稳处选点;在相位连续处选点,要远离高频相位;要避开个别解缠错误的相位点。
SBAS-InSAR第一次反演:基于线性模型计算出所有像对的形变(日期、速度、加速度和加速度变化)和高程(校正值和新的DEM)。
SBAS-InSAR第二次反演:计算时间序列上的位移,得到研究区地表的时序形变量。在第一步得到的形变速率基础上,进行大气滤波,估算和去除大气相位,得到最终位移结果。
地理编码:以上成果是在斜距坐标系下,此步可将成果的斜距坐标转换为地理坐标。本项目选择WGS-84坐标系,并得到LOS方向的形变结果。输出的结果还包括对应时期的累积形变量、最终的平均形变速率以及平均形变速率精度、平均形变加速度等。
4. 成果分析
本项目共监测出沉降区域十三处,沉降最严重地区沉降速率可达到90 mm/年。辽河平原整体沉降速率图见图3。部分沉降严重区域沉降速率图见图4~6。从图4中可看出,该区域大部分地区地表稳定,沉降值在10 mm以内,但中心区域的村庄发生了较大的沉降,最大可达到60 mm;图5中大部分地区沉降值在10 mm以内,但有多个村庄沉降值达到40 mm;图6中中心区域沉降达到60 mm,且周围有多处沉降达到30。正常地表沉降在10 mm以内我们都认为是正常现象,20 mm以上是我们的关注区。通过外业踏勘我们发现,沉降达到60 mm以上的地区,基本都在地下资源开采区。因此,原沉降成因主要受三方面影响,一是地下煤炭资源开采。比如小康矿、大平矿的煤炭开采在三台子水库附近形成大面积的沉降区;调兵山市多个煤矿在全市大范围内引发整体沉降;西马煤矿在西马峰镇的西大堡村和野老村等地形成了严重的沉降区。二是地下油气资源开采。辽河油田有大量采油厂,多年以来的持续开采,致使周边区域地下油气储层等液体资源储层大面积采空,且开采范围逐步扩大,地面沉降范围和沉降量值也将逐渐随时间扩大。同时由于油气储层埋藏深,采空区面积巨大,因而一般以缓慢、累进性的整体沉陷为主。三是地下水资源开采。根据资料实地调查取证,沉降区域内大面积种植水稻、玉米等农作物以及草莓等经济作物,在农作物生长周期内需水量巨大。调查区内水位埋深浅,平均在0.5~5 m左右,因其灌溉水源和饮用水源往往来自于深层地下承压水的开采,且其开采之后并无回灌,造成沉降区承压水位持续下降,土体有效应力增加产生压缩变形,从而表现为以开采中心为原点的地面沉降漏斗。
Figure 3. Overall settlement rate of Liaohe Plain
图3. 辽河平原整体沉降速率图
Figure 4. Monitoring results of settlement area I
图4. 沉降区域一监测成果
Figure 5. Monitoring results of settlement area II
图5. 沉降区域二监测成果
Figure 6. Monitoring results of settlement area III
图6. 沉降区域三监测成果
5. 结束语
通过本次项目作业,得出以下结论:
1) SBAS-InSAR技术对于地表沉降变化的监测是有效的,其结果的精度基本可以达到毫米级。
2) SBAS-InSAR技术特别适合对面积大的区域进行沉降监测。
3) SBAS-InSAR技术对于地表发生沉降变化尤其是长时间的缓慢变化非常敏感。
4) SBAS-InSAR技术成本低,目前国际上有哨兵免费的SAR数据源。
5) 本次作业发现辽河平原地区确实存在多处沉降区,其中沉降最大的位置在辽河油田,其次是煤矿区。对于地下资源开采的同时,采取科学的方法减少对我们生活环境的破坏是必要的。
参考文献