龙里某塌陷时序InSAR变形监测的PS修正

doi:10.6046/zrzyyg.2021334

龙里某塌陷时序InSAR变形监测的PS修正

罗雪玮^,¹, 向喜琼^,¹, 吕亚东¹

1.贵州大学喀斯特地质资源与环境教育部重点实验室,贵阳 550025

2.贵州大学资源与环境工程学院,贵阳 550025

PS correction of InSAR time series deformation monitoring for a certain collapse in Longli County

LUO Xuewei^,¹, XIANG Xiqiong^,¹, LYU Yadong¹

1. Key Laboratory of Karst Georesources and Environment, Ministry of Education, Guizhou University, Guiyang 550025, China

2. College of Resources and Environmental Engineering, Guizhou University, Guiyang 550025, China

通讯作者: 向喜琼(1975-),男,博士,副教授,主要从事地质灾害评价与防治教学研究。Email:tujia@126.com。

责任编辑: 李瑜

收稿日期: 2021-10-12 修回日期: 2022-03-14

基金资助:

中央引导地方科技发展资金项目(黔科中引地[2021]4001)

Received: 2021-10-12 Revised: 2022-03-14

作者简介 About authors

罗雪玮(1997-),女,硕士研究生,主要从事地球探测与信息技术研究。Email: lxwei0608@163.com。

摘要

基于永久散射体的信号在很长时间范围内都能保持较高的干涉相干性的优点,为了解决小基线集合成孔径雷达干涉测量技术(small baseline subset interferometric synthetic aperture Radar,SBAS-InSAR)轨道精炼步骤时人工选择地面控制点可能会影响到监测结果这一问题,首先,该文将永久散射体与SBAS-InSAR结合,通过设置相干系数的阈值、振幅离差指数的阈值以及地表形变速率的阈值选出稳定的永久散射体,并将这些点作为SBAS-InSAR轨道精炼中的地面控制点,从而修正监测结果的准确度; 然后,选用2019年9月1日—2021年4月11日20景覆盖贵州省龙里县洗马镇的Sentinel-1A双极化影像为主要数据源,进行地表形变监测; 最后,将该方法所得结果、人工选择地面控制点的方法所得结果与北斗位移监测数据进行对比分析,可知该文方法比人工选择地面控制点的方法更精准。

关键词： 永久散射体; 多重阈值; SBAS-InSAR; 地面控制点

Abstract

The signals of permanent scatterers can maintain high interferometric coherence for a long time. To solve the problem that manually selecting ground control points may affect the monitoring results during the orbit refinement of the SBAS-InSAR, this study combined permanent scatterers with SBAS-InSAR. Firstly, by setting the thresholds of the coherence coefficient, the amplitude dispersion index, and the surface deformation rate, this study selected robust permanent scatterers as the ground control points in the orbit refinement of the SBAS-InSAR in order to correct the accuracy of the monitoring results. Then, this study selected 20 scenes of Sentinel-1A dual-polarization images that covered Xima Town, Longli County, Guizhou Province from September 1, 2019 to April 11, 2021 as the main data source for surface deformation monitoring. Finally, this study compared the results obtained using the proposed method and those obtained through manually selecting ground control points with the displacement monitoring data of the Beidou satellite, concluding that the data obtained using the method proposed in this study were more accurate.

Keywords： permanent scatterer; multi-threshold; SBAS-InSAR; ground control points

PDF (3674KB) 元数据多维度评价相关文章导出 EndNote| Ris| Bibtex 收藏本文

本文引用格式

罗雪玮, 向喜琼, 吕亚东. 龙里某塌陷时序InSAR变形监测的PS修正[J]. 自然资源遥感, 2022, 34(3): 82-87 doi:10.6046/zrzyyg.2021334

LUO Xuewei, XIANG Xiqiong, LYU Yadong. PS correction of InSAR time series deformation monitoring for a certain collapse in Longli County[J]. Remote Sensing for Land & Resources, 2022, 34(3): 82-87 doi:10.6046/zrzyyg.2021334

0 引言

为克服差分合成孔径雷达干涉测量技术(differential interferometric synthetic aperture Radar,D-InSAR)时空失相干、大气效应及噪声等因素的影响,一系列的时序合成孔径雷达干涉测量技术(interferometric synthetic aperture Radar,InSAR)技术应运而生。具有代表性的有Stacking-InSAR、永久散射体(permanent scatter,PS)PS-InSAR、小基线集(small baseline subset,SBAS)SBAS-InSAR以及相干永久目标分析(interferometric persistent target analysis)IPTA-InSAR等技术^[1]。Colombo等^[2]利用ESA ERS1/2卫星采集的合成孔径雷达(synthetic aperture Radar,SAR)数据,采用D-InSAR和PS-InSAR进行了干涉分析,在监测的工业区内发现了不同的地面沉降模式; Wu等^[3]采用SBAS-InSAR方法对延安新区2015—2019年的地表变形进行了调查,发现区内地基沉降主要来源于填方和城市荷载的作用,而建筑物的应力释放是造成地面隆起的主要因素; 刘一霖等^[4]基于SBAS-InSAR技术,引入偏移量追踪法、FFT过采样方法、多窗口迭代自适应滤波技术对矿区进行了地表时序监测,分析了研究区开采工作面内地表的时空变化规律; 葛大庆等^[5]融合了PS-InSAR和SBAS-InSAR两种方法,以短基线差分干涉生成相位图后利用点目标识别技术提取出相干目标,分析了德州地区地面沉降-回弹的动态变化过程; 李金超等^[6]将基于灰色模型和支持向量回归的组合模型(gray model and support vector regression,MG-SVR)与SBAS-InSAR技术相结合,对淮南潘集矿区发生沉降的区域进行了监测及预警。

在使用SBAS方法对地表形变进行监测时,需要在轨道精炼步骤中引入稳定的地面控制点(ground control point,GCP)估算和去除残余的恒定相位和解缠后还存在的相位坡道,从而使形变监测结果更为准确。但在没有精确GCP的情形下,如果研究区内存在较大高差、或者影像数据噪声过大,无法进行消除,人工来选择GCP可能会引起较大误差。

针对此问题,本文提出了一种基于PS的SBAS-InSAR方法。该方法通过设置相干系数的阈值、振幅离差指数的阈值以及地表形变速率的阈值选出稳定的PS点,并将这些点作为SBAS-InSAR轨道精炼中的地面控制点,以期能得到准确度更高的地表形变监测结果。最后,以位于贵州省龙里县洗马镇上石坎村的一处地面塌陷为实例验证了本文方法的有效性。

1 研究方法

1.1 相干系数法

相干系数法的基本思路为: 对于某一像元选定一定大小的移动窗口(m,n),利用该移动窗口内的像素复数信息来估计相干系数值。计算完毕后,移动该窗口,重复计算^[7-8]。相关系数 $γ$ 计算公式为:

(1)

\begin{matrix} γ = \frac{\overset{m}{\sum_{i = 1}} \overset{n}{\sum_{j = 1}} M (i, j) S^{*} (i, j)}{\sqrt{\overset{n}{\sum_{j = 1}} {|M (i, j)|}^{2} \overset{m}{\sum_{i = 1}} \overset{n}{\sum_{j = 1}} {|S^{*} (i, j)|}^{2}}} \end{matrix}

式中: M为主影像; S为辅影像; ^*为共轭相乘。

计算得到每个像元的 $γ$ 值后,再求得每个像元时序上的平均值 $\bar{γ}$ ,公式为:

(2)

\begin{matrix} \bar{γ} = \frac{1}{N} \overset{N}{\sum_{i = 1}} γ_{i} \end{matrix}

。

一般情况下, $γ$ 值越大,干涉相位的噪声越低,设定阈值 $T_{γ}$ ,当某一像素点的 $\bar{γ}$ 值大于 $T_{γ}$ 时,将该点判定为有效PS点; 否则视为无效PS点^[7-8]。

1.2 振幅离差法

Ferretti等^[9]提出可以采用像素中振幅信息的时间序列来度量干涉相位的稳定性。R和I分别为影像的实部和虚部,若存在标准差为 $σ_{n}$ 的高斯噪声,则振幅值A服从Rice分布,公式为:

(3)

\begin{matrix} f (A) = \frac{A}{σ_{n}^{2}} I_{0} (\frac{A g}{σ_{n}^{2}}) e^{- \frac{A^{2} + g^{2}}{2 σ_{n}^{2}}} A > 0 \end{matrix}

式中: g为目标反射量,且g>0; I₀为Bessel函数^[6]。当信噪比 $\frac{g}{σ_{n}}$ 较小时,Rice分布演化为Rayleigh分布; 而在高信噪比目标( $\frac{g}{σ_{n}}$ >4)时,该分布趋于高斯分布。因此,当 $σ_{n} < < g$ 时,相位离差指数可以用振幅离差指数来估计,公式为:

(4)

\begin{matrix} σ_{v} ≅ \frac{σ_{n I}}{g} ≅ \frac{σ_{A}}{m_{A}} ≜ D_{A} \end{matrix}

式中: v为干涉相位; $σ_{v}$ 为相位离差指数; $σ_{n I}$ 为虚部的标准差; $m_{A}$ 和 $σ_{A}$ 为时序振幅的均值和标准差; $D_{A}$ 为振幅离差指数。

1.3 相干性、振幅离差指数以及形变速率结合的选点方法

首先采用相干系数阈值法初步识别出具有高相干系数的像元的PS点目标; 然后设定振幅离差值进一步筛选具有稳定强散射性的PS点作为目标; 最后设定一个形变速率区间来进行最终的PS点选取。如图1所示。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1 方法流程

Fig.1 Method and procedure

1.4 SBAS-InSAR技术

SBAS-InSAR技术可以限制时空基线失相干和大气延迟的影响^[10]。将N+1景覆盖同一区域的SAR影像序列,影像获取时间为 $(t_{0}, \dots, t_{N})$ ,根据空间基线和时间基线的条件组合干涉对得到M幅差分干涉图,M的数量与影像数量满足以下关系^[11]:

(5)

\begin{matrix} \frac{N + 1}{2} \leq M \leq N (\frac{N + 1}{2}) \end{matrix}

。

假设干涉图j由t_A和t_B时刻获取的影像干涉生成,在去除平地效应与地形相位影响后,干涉图j中距离向为r与方位向为x的某一像素的干涉相位表示为:

(6)

\begin{array}{l} δ ϕ_{j} (x, r) = ϕ (t_{B}, x, r) - ϕ (t_{A}, x, r) \\ \approx ϕ_{d e f, j} (x, r) + ϕ_{t o p p o, j} (x, r) \\ + ϕ_{a t m, j} (x, r) + ϕ_{n o i s e, j} (x, r) \end{array}

式中: $ϕ (t_{B}, x, r)$ 和 $ϕ (t_{A}, x, r)$ 分别为t_B和t_A时刻SAR影像上的相位值; $ϕ_{d e f, j} (x, r)$ 为t_A~t_B时刻内卫星视线向(los)的形变相位; $ϕ_{t o p p o, j} (x, r)$ 为地形相位误差; $ϕ_{a t m, j} (x, r)$ 为大气相位误差; $ϕ_{n o i s e, j} (x, r)$ 为噪声相位。

对所有干涉图构建差分干涉相位与影像获取时刻形变量之间的线性方程组,其矩阵形式为:

(7)

\begin{matrix} A ϕ = δ ϕ \end{matrix}

式中: $A$ 为M×N的系数矩阵; $ϕ$ 为参数矩阵,由未知像元(x,y)在N个时刻对应的未知形变相位值构成; $δ ϕ$ 为由M个式(6)所示的δϕ_j组成的观测方程组; M为差分干涉图的数量。

为求解研究区域各高相干点的形变速率,可用两景影像间的平均相位速率代替相位值,则式(7)变为:

(8)

\begin{matrix} B v = δ ϕ \end{matrix}

式中: $B$ 为M×N的系数矩阵; $v$ 为形变速率向量, $v^{Τ}$ 可表示为:

(9)

\begin{matrix} v^{Τ} = [\begin{array}{l} v_{1} = \frac{φ_{1}}{t_{1} - t_{0}}, \\ \dots, v_{N} = \frac{φ_{N} - φ_{N - 1}}{t_{N} - t_{N - 1}} \end{array}] \end{matrix}

。

当系数矩阵 $B$ 为满秩(即M≥N)时,可用最小二乘法求解形变速率; 当M<N时,矩阵 $B$ 出现秩亏,可用奇异值分解法获取研究区域的形变速率。获取形变速率后,依据研究区SAR影像间的时间区间可以求得相应时间段内的形变量。

2 实验数据和处理

2.1 数据

本文以龙里县洗马镇的一处地面塌陷为实例(图2),选取20景由欧洲航天局提供的C波段Sentinel-1A SAR影像数据,时间范围为2019年9月1日—2021年4月11日。选用POD精密定轨星历数据去除轨道相位,采用SRTM-DEM去除地形相位并进行坐标系转换。表1列举了Sentinel-1A影像数据的基本参数。

图2

新窗口打开| 下载原图ZIP| 生成PPT

图2 Sentinel-1A影像位置

Fig.2 Sentinel-1A image position

表1 Sentinel-1A数据参数

Tab.1 Sentinel-1A parameter

参数	值
数据类型成像模式重访周期/d 极化方式距离分辨率/m 方位分辨率/m 入射角/(°) 轨道方向	SLC IW 12 VV 5 20 44 升轨

新窗口打开| 下载CSV

2.2 PS点的选取

首先对20景SAR数据进行配准、干涉、去平处理并计算出数据的时序振幅标准差与平均值。先采用相干系数法识别出具有高相干系数的像元作为PS点选取目标,考虑到研究区为城镇,建筑物较多,故设置相干性阈值为0.75; 然后使用振幅离差指数进一步筛选出具有稳定相位的PS点,根据Ferretti等^[9]的研究,通常振幅离差指数只考虑小于0.25的像素点,故设置振幅离差指数阈值为0.25; 最后基于形变速率来确定PS点,因为本文选取的PS点需要作为SBAS-InSAR轨道精炼时的GCP点,而GCP点的基本选取原则为远离形变区域,故本文将形变速率区间设置为[-0.1,0.1] mm/a。基于以上3个限制条件,共选出相位稳定的PS点36个,如图3所示,且将PS点导入地图对比查看后发现大部分PS点都位于稳定的建筑物(道路、屋顶等)上。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3 基于3个阈值得到的PS点

Fig.3 PS points based on the three thresholds

2.3 SBAS-InSAR处理

将得到的PS点由地理坐标转为雷达斜距坐标,然后将其作为GCP点与干涉处理后的像对进行轨道精炼,估算和去除残余的恒定相位和解缠后还存在的相位坡道; 接着通过2次反演估算形变速率和残余地形,并去除大气相位; 最后获取了研究区地表形变速率(图4)和时序形变结果(图5)。可以看出: 区内平均地表形变速率在-6.5~5.7 mm/a之间,年平均最大沉降速率为-6.52 mm,说明该地变形程度轻微,且监测到的形变主要集中在塌陷区边界及道路周围,与实地地质调查的结果相符合。除开植被覆盖区域出现失相干外,塌陷区内大部分农田没有监测到相关形变是由于春耕时节居民对农田进行翻土种地等操作,农田里长出了新作物,所以也导致了该区域的部分失相干情况。以2019年9月1日为第一时相,其余各期影像相对于第一时相的形变如图5所示。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4 SBAS-InSAR形变速率图

Fig.4 Deformation rate diagram of SBAS-InSAR

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5 研究区地表形变时间序列

Fig.5 Time series of surface deformation in the study area

3 结果评述

为了验证该方法的优越性,本文使用相同时间段的SAR影像数据,基于选择GCP点的重要标准^[12]: ①没有残余地形条纹; ②没有形变条纹,远离形变区域; ③没有相位跃变; ④至少20~30个点。通过对比滤波后的干涉图和解缠后的干涉图,在没有相位跃变和残余地形相位的区域,人工手动选择了30个控制点,得到研究区的形变速率。并收集了研究区内设置的3处北斗位移监测站点的监测结果(选取时间分别为2020年11月6日,2020年12月12日,2021年1月17日,北斗监测结果的时间跨度为2020年10月28日—2021年11月25日,与本文InSAR时间序列相对一致),然后在2种SBAS形变结果中选取了监测站点周围不超过50 m的3个变形点,计算这3个变形点在同一位移监测时间内的相对垂直向形变量,并与位移监测数据进行比较,比较结果如表2,图中点号为北斗监测的站点名。

表2 2种方法与北斗监测点比较情况

Tab.2 Comparison between the two methods and the Beidou monitoring station(mm)

点号	位移监测数据	本文方法		传统SBAS方法
点号	位移监测数据	形变量	误差	形变量	误差
SC-07	-2.41	1.07	3.48	4.53	6.94
SC-06	-0.94	-0.63	0.31	-0.11	0.83
SC-05	-2.90	-5.79	-2.89	-9.28	-6.99

新窗口打开| 下载CSV

由表2看出,2种方法中,结合了PS点的SBAS-InSAR方法的监测结果比人工选取GCP点的传统SBAS-InSAR方法更接近位移监测数据,且2种监测结果与位移监测数据均有一定的误差。通过查阅文献再结合实地考察,分析造成误差的原因主要有以下3点: ①由于客观条件的限制,研究区内北斗位移监测站点的设置与所选择的InSAR形变点的位置并不是完全重合的,故选择了在监测站点50 m以内的3个邻近InSAR形变点来进行对比分析; ②InSAR形变结果的分析处理过程是通过像元而不是单个点来处理的,像元覆盖了一定的范围,所以得到的结果可能会受到周围临近像元及整体区域上的影响而造成误差; ③除了InSAR技术本身存在的误差源,在时序InSAR中还存在雷达影像对精度的影响、干涉相位对函数模型观测量与未知量之间函数关系的影响^[13]。InSAR的相位组成可通过函数表现为以下形式:

(10)

\begin{array}{l} φ = φ_{d e f o} + φ_{a t m o} + φ_{o r b i t} \\ + φ_{d e m} + ε = A x + ε \end{array}

式中: $ε$ 为失相干随机噪声, $E {φ} = 0$ 。基于时序InSAR技术,可以将多幅干涉图组成其函数模型观测量,增加多余观测,从而进行误差的解算与分离。

4 结论

1)在进行SBAS-InSAR处理时,为了使选择的GCP点不引起较大误差从而影响到形变结果,本文通过设置相干系数的阈值、振幅离差指数的阈值以及地表形变速率的阈值选取了稳定的PS点,然后将其作为GCP点进行轨道精炼及后续的反演处理得到地表形变信息。并将本文方法与采用人工选取GCP点的方法所得结果和北斗位移监测结果进行对比。由本文方法得到的形变结果比通过人工选择GCP的SBAS方法得到的形变结果更为精准,能更加有效地分析研究区地表时序形变。在实例研究中,研究区内平均地表形变速率在-6.52~5.66 mm/a之间,年平均最大形变速率为-6.52 mm,该地整体变形程度轻微,变形主要集中在塌陷区边界及道路周围。

2)由于本文方法所选取的PS点大部分都位于建筑物上,当建筑物没有发生竖直方向上的变形时,也会被判定为PS点。那么基于这样的PS点得到的时序形变结果也可能会出现较大误差。因此,后续的研究可以通过将升降轨数据进行融合得到三维的形变量,再进一步解算水平与竖直方向的形变,通过共同判断x和y方向上的位移,这样得到的永久散射体就会更为可靠。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

曹翠华

. InSAR技术监测西南喀斯特地表形变应用研究[D]. 西安: 西安科技大学, 2020.