面向地表形变高精度监测的GNSS-InSAR融合方法

doi:10.11873/j.issn.1004-0323.2023.5.1028

面向地表形变高精度监测的GNSS-InSAR融合方法

柯福阳^,¹, 胡祥祥¹, 明璐璐¹, 刘学武², 尹继鑫², 刘宇航²

1.南京信息工程大学遥感与测绘工程学院，江苏南京 210440

2.西宁市测绘院，青海西宁 810000

GNSS-INSAR Fusion Method for High Precision Monitoring of Surface Deformation

KE Fuyang^,¹, HU Xiangxiang¹, MING Lulu¹, LIU Xuewu², YIN Jixin², LIU Yuhang²

1.School of Remote Sensing and Geomatics Engineering，Nanjing University of Information Science and Technology，Nanjing 210440，China

2.Xining Surveying and Mapping Institute，Xining 810000，China

收稿日期: 2022-05-30 修回日期: 2023-09-08

基金资助:

2022年度第六期“333人才”培养支持资助项目.  BRA2022042
江苏省“六大人才高峰”高层次人才项目.  XYDDX⁃045
西宁市科技计划项目.  2019⁃Y⁃12

Received: 2022-05-30 Revised: 2023-09-08

作者简介 About authors

柯福阳（1981-），男，福建惠安人，博士，教授，主要从事地质与气象灾害监测预警研究E⁃mail：kfy_0829@163.com , E-mail：kfy_0829@163.com

摘要

GNSS-InSAR数据融合进行监测地表形变是目前地表形变监测领域研究的热点问题，传统GNSS-InSAR数据融合方法融合简单、不能动态地反映地表形变特点，导致数据使用不充分、形变特征精度低等后果。提出了一种新的基于InSAR校正值和卡尔曼滤波的GNSS-InSAR融合方法。根据时间序列的GNSS观测值和InSAR校正观测值的时空相关性，通过卡尔曼滤波对两种数据进行融合，得到更精确的地表三维形变结果。利用2018年11月15日至2022年6月3日103景Sentinel-1A数据和同期13个GNSS点位数据进行处理，实验结果表明：校正后的InSAR观测值与GNSS观测值经卡尔曼滤波融合结果比未校正的InSAR观测值与GNSS观测值融合结果精度高45%，比InSAR观测值精度高57%。因此，基于InSAR校正值和卡尔曼滤波的GNSS-InSAR融合模型提高了InSAR变形监测的精度，拓展提升InSAR应用范围的广度和深度。

关键词： GNSS-InSAR融合 ; 拟合推估法 ; 卡尔曼方法 ; 地表三维形变测量 ; 地表形变

Abstract

Surface deformation is a geological phenomenon caused by natural or artificial factors， and its disaster-causing process is slow and irreversible. It is also a geological disaster with destructive solid power. Therefore， real-time and high-precision surface deformation monitoring is one of the most critical tasks in maintaining urban safety. However， due to the complex causes， long duration， wide range， and many triggering factors of surface deformation， there are many difficulties in monitoring surface deformation using single technology such as leveling， GNSS， INSAR， and optical remote sensing. Considering the characteristics and complementarities of InSAR and GNSS， the combination of InSAR and BeiDou/GNSS can improve the surface deformation monitoring capability in space and time at the same time. Unluckily， the traditional GNSS-InSAR data fusion method is simple to fuse and cannot dynamically reflect surface deformation characteristics， leading to insufficient data use and low accuracy of deformation features. A new fusion method is proposed based on the Kalman filter algorithm GNSS-InSAR correction values. The method mainly consists of two sequential processes， i.e.， the a priori processing of GNSS and INSAR data and the fusion process of GNSS-InSAR correction values based on the Kalman filter algorithm. The a priori processing of GNSS and INSAR data is to obtain the a priori deformation results using the fitted estimation model to correct the systematic errors in the InSAR observations. The fusion process of GNSS-InSAR correction values based on the Kalman filtering algorithm is to fuse the two data through Kalman filtering based on the spatial and temporal correlation between the time-series GNSS observations and the InSAR correction observations. The experiment was processed using 103 views of sentinel-1A data from November 15， 2018， to June 3， 2022， and 13 GNSS point data during the same period. The experimental results show that the fusion result of the corrected InSAR observations and GNSS observations by the Kalman filter is 45% more accurate than the fusion result of the uncorrected InSAR observations and the GNSS observations， which is 45% 57% higher than the accuracy of InSAR observations. Therefore， the fusion method model based on the Kalman filter algorithm of GNSS-InSAR corrected values proposed in this paper improves the accuracy of InSAR deformation monitoring and expands the breadth and depth of InSAR applications.

Keywords： GNSS-InSAR data fusion method ; Fitting estimation method ; Kalman filtering ; 3D surface deformation measurement ; The surface deformation

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

本文引用格式

柯福阳, 胡祥祥, 明璐璐, 刘学武, 尹继鑫, 刘宇航. 面向地表形变高精度监测的GNSS-InSAR融合方法. 遥感技术与应用[J], 2023, 38(5): 1028-1041 doi:10.11873/j.issn.1004-0323.2023.5.1028

KE Fuyang, HU Xiangxiang, MING Lulu, LIU Xuewu, YIN Jixin, LIU Yuhang. GNSS-INSAR Fusion Method for High Precision Monitoring of Surface Deformation. Remote Sensing Technology and Application[J], 2023, 38(5): 1028-1041 doi:10.11873/j.issn.1004-0323.2023.5.1028

1 引言

随着城市基础设施建设规模的逐步扩大，地质灾害频发。常见的诸如泥石流、滑坡、地面沉陷等地质灾害多因为地表形变引起的^［1-2］，由此对地表三维形变监测^［3］具有重大现实意义。目前，地表形变监测方法主要有全球导航卫星定位系统（GNSS）和卫星雷达干涉测量^［4］（InSAR）两种。其中，GNSS 技术的优势在于具有多方向、高精度、自动化观测等特点，但受限于监测设备，其GNSS点位稀疏，一些灾害点人工无法布设控制点^［5］。InSAR技术的优势在于其具有大范围、高精度、不需要地面架设设备等特点^［6］。人工无法布设控制点的地方，InSAR都可以覆盖到。但InSAR受SAR卫星观测周期的影响^［7］，同一个位置监测周期长，且只能获取该位置视线向（LOS）上的一维形变，制约了InSAR监测精度^［8］。如何利用两者不同优点进行融合，解决高时空分辨率、高精度的三维形变场问题成为地表形变监测领域研究的热点^［9-10］。Bock等^［11］率先提出GNSS和InSAR的融合思想，Ge等^［12］提出双插双估计（DIDP）方法，该方法利用GNSS数据来提高InSAR观测值的精度，再利用InSAR反过来加密GNSS观测结果。罗海滨等^［13］提出了直接分解法，该方法成功利用GNSS和InSAR的优点，算法简单，但是十分依赖GNSS的精度和GNSS插值的精度，若GNSS精度或GNSS插值精度较差，最后得到的结果会受到影响。胡俊等^［14］提出了赫尔默特方差估计法，该方法根据GNSS和InSAR观测值方差的估计值反复定权，以达到两者观测值的权合理，其结果较最小二乘法有提高。汪友军等^［15］利用应力应变模型结合方差分量估计算法，得到更高精度的三维地表形变结果。赵增鹏等^［16］分析了对流层对InSAR的影响，通过GNSS数据处理InSAR数据中存在的对流层延迟误差，提高了D-InSAR的监测精度。江克贵等^［17］提出了利用BK模型结合单视线D-InSAR的监测方法，但由于观测区域较远，可能导致了预计效果不佳。吕佳凝等^［18］在直接分解法中进一步提出拟合推估法^［19］。该方法同样利用了GNSS数据的特性，对InSAR影像中存在的系统误差进行校正，提高了融合精度。但是方法末尾结合了最小二乘法，当有较大粗差出现时，将不能得到较好的结果。熊陆云等^［20］提出了一种基于贝叶斯定理的最大后验估计准则的虚拟观测迭代最小二乘（ VOILS ）方法来校正GPS位移内插过程引起的误差。其对垂直方向的改善尤为明显。雷坤超等^［21］利用GNSS水平精度高和InSAR垂直精度高的原理，将GNSS水平位移测量与InSAR垂直观测结合，得到了北京平原区的地表三维形变场。但是，该方法降低了GNSS的时间分辨率和InSAR的空间分辨率。薛学明等^［22］基于无迹卡尔曼滤波，结合GNSS和InSAR时间序列观测模拟了火山变形，GNSS和InSAR数据拟合结果在2~4 km深度处拟合的较好。武帅莹等^［23］基于机器学习的GNSS大气建模方法对InSAR进行大气改正，提升InSAR监测精度。

常规GNSS-InSAR融合方法多为简单融合，不能时序地反映地表三维形变特征。已有的卡尔曼滤波融合方法，都没有很好的解决数据本身存在的误差，融合后结果精度较低。本文提出一种新的基于InSAR校正值和卡尔曼滤波的GNSS-InSAR融合方法。该方法通过GNSS数据和InSAR数据存在时空相关性，利用GNSS数据对InSAR观测值中存在的误差进行校正。然后利用GNSS、InSAR数据之间的时空关联建立方程，基于卡尔曼滤波模型融合GNSS、InSAR观测数据。使用西宁市的GNSS和InSAR数据进行了实测实验，实测实验证明了该方法的有效性。

2 融合模型技术原理

2.1　拟合推估法函数模型

拟合推估法是一种推估方法。该方法是通过一定的拟合法则，对数据进行推估。由于SAR卫星存在雷达参数、相位观测量和地形数据等不确定性，其误差会传播到形变结果中去^［24］。若直接与GNSS数据融合，其结果也包含这些误差。但是，InSAR的系统误差与GNSS存在较高的空间相关性^［25］，可通过GNSS数据进行削弱^［26-27］。将高精度的GNSS观测值对InSAR观测值进行约束，得到同点位上GNSS和InSAR观测值的差值，通过拟合法则将差值拟合出来，将差值进行校正，最后得到精度较高的InSAR观测值^［28］。

为了对InSAR观测值进行校正，需要将GNSS观测值投影到LOS向：

L_{G P S}^{l o s} = S_{x} * E + S_{y} * N + S_{z} * U

（1）

其中： $L_{G P S}^{l o s}$ 为监测点GNSS视线向形变量； $E$ 、 $N$ 、 $U$ 为监测点东西向、南北向、垂直向上的形变量； $S_{x}$ 、 $S_{y}$ 、 $S_{z}$ 分别是InSAR观测值的单位投影矢量。其中 $α$ 是卫星航向方位角（以正北方向为起点顺时针旋转）， $θ$ 是雷达入射角。 $S_{x} = - s i n θ * s i n (α - \frac{3}{2} π)$ ， $S_{y} = - s i n θ * c o s (α - \frac{3}{2} π)$ ， $S_{z} = c o s θ$ 。

用GNSS观测值对InSAR观测值进行约束后得到两者的差值，如式（2），利用式（3）将差值进行拟合：

V = L_{G P S}^{l o s} - L_{I n S A R}^{l o s}

（2）

V = a_{0} + a_{1} * B + a_{2} * L + a_{3} * B^{2}

+ a_{4} * L^{2} + a_{5} * B * L + Δ

（3）

其中： $L_{I n S A R}^{l o s}$ 为InSAR观测值； $B$ 为监测点的纬度； $L$ 为监测点的经度； $V$ 为GNSS归算值与InSAR的差值； $a_{0}$ 、 $a_{1}$ 、 $a_{2}$ 、 $a_{3}$ 、 $a_{4}$ 、 $a_{5}$ 为二次拟合模型参数； $Δ$ 为观测噪声。

通过上述公式解方程组即可解算得到待求参数 $a_{0}$ 、 $a_{1}$ 、 $a_{2}$ 、 $a_{3}$ 、 $a_{4}$ 、 $a_{5}$ ，将任一点位的经纬度带入上述公式后，即可得到该点的差值，然后将InSAR观测值进行补偿得到改正后的InSAR视线向观测值 ${\bar{L}}_{I n S A R}^{l o s}$ 。

2.2　基于卡尔曼滤波的GNSS-InSAR融合模型

卡尔曼滤波是一种能对系统状态进行最优估计的算法。利用卡尔曼滤波对GNSS和InSAR数据进行融合，首先需要建立卡尔曼滤波观测方程和状态方程。其观测方程如下：

L_{k} = A * X_{k} + Δ

（4）

其中： $L_{k} = [L_{G P S}^{E} L_{G P S}^{N} L_{G P S}^{U} {L_{I n S A R}]}^{T}$ 为观测量； $X_{k} = {[d_{e} d_{n} d_{u} v_{e} v_{n} v_{u}]}^{T}$ 为状态矢量； $Δ$ 为观测噪声。其中 $d$ 为监测点的形变量， $v$ 为监测点的变形速率，3个状态向量通过最小二乘法求得^［29］。 $A = {[\begin{matrix} I_{3} & 0 * I_{3} \\ S & 0 \end{matrix}]}_{6 * 6}$ 为观测方程设计矩阵， $I_{3}$ 为一个3×3的单位阵。

其状态方程如下所示：

X_{k} = φ_{k � k - 1} * X_{k - 1} + Γ_{k - 1} * W_{k - 1}

(5)

其中： $Φ_{k � k - 1} = {[\begin{matrix} I_{3} & Δ t * I_{3} \\ 0 & I_{3} \end{matrix}]}_{6 * 6}$ 为状态转移矩阵； $Γ_{k - 1} = {[\begin{matrix} \frac{1}{2} * Δ t^{2} * I_{3}^{} \\ Δ t * I_{3} \end{matrix}]}_{6 * 3}$ 为系统动态过程噪声分布矩阵， $W_{k - 1}$ 为系统动态过程噪声向量； $W_{k - 1} = {[1 1 1]}^{T}$ 或者取3个方向加速度的平均值；经验证 $W_{k - 1}$ 的改变对结果影响不大。

在状态方程和观测方程建好后，需给定系统初始状态。取初始时刻的形变量及方差为0，初始时刻的 $X_{0}$ 及对应的方差矩阵 $Σ_{x_{0}}$ 为GNSS和InSAR数据平差后的形变速率和方差。观测噪声的方差阵^［30］ $Σ_{w} = 4 Δ t_{k}^{- 4} \sum_{0}$ 。

2.3　融合模型流程

如图1所示，先利用拟合推估模型对InSAR观测值中的系统误差进行校正，将校正后的InSAR观测值与对应的GNSS观测值进行基于卡尔曼滤波的融合，最后得到地表三维形变结果。

图1

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

图1 融合流程图

Fig.1 Flow chart of fusion algorithm

3 研究区域及研究数据

3.1　研究区域

西宁位于中国西北，四面环山，东西呈条带状。湟水两岸（分别为南山和北山）沟壑纵横，地形切割强烈，人类工程活动强烈，生态地质环境脆弱，崩塌、滑坡、泥石流等突发性地质灾害发育，特别是滑坡、崩塌高发，对坡脚居民的生命财产构成严重威胁^［31-32］。研究区和GNSS站点分布如图2所示。

图2

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

图2 研究区及GNSS站点示意图

Fig.2 Schematic diagram of the research area and GNSS site

3.2　SAR数据

研究采用的数据是由欧空局（ESA）发布的2018年11月15日至2022年6月3日103景Sentinel-1A卫星降轨数据以及美国国家航空航天局发布的分辨率为90 m的SrtmV4DEM数据，相邻影像时间间隔为12 d。影像覆盖范围如图2黄色方框。InSAR卫星传感器相关数据见表1。

表1 InSAR卫星相关数据

Table 1 InSAR satellite related data

入射角	方位角	极化方式	轨道	时间间隔
34.039	193.223	VV	降轨	20181115~20220603

新窗口打开| 下载CSV

3.3　GNSS数据

研究使用的GNSS数据来自课题组承担的西宁市南山滑坡体自动化监测项目。为了缓解城区交通压力，西宁市规划设计了凤凰山快速路，快速路穿越南山山体。为保证南山隧道周边建筑和人民生命财产安全，在南山西坡布置了GNSS监测站，对南山滑坡体进行自动化变形监测。GNSS站点如图2所示。GNSS数据平面精度 $\pm (2.5 + 1 \times 10^{- 6} \times D)$ mm，高程精度 $\pm (5 + 1 \times 10^{- 6} \times D)$ mm，其中D为基线距离（单位：km）。基于开源rtklib开发的高精度GNSS变形监测软件，在剔除一些粗差较大和数据缺失的点位后，得到了与InSAR同步观测的13个GNSS变形结果。

4 数据处理与分析

研究使用ENVI SARScape软件，利用短基线集技术（Small Baseline Subset，SBAS）对所得103景Sentinel-1A卫星降轨数据和对应精轨数据进行处理。①设置时间基线、空间基线^［33］阈值等相关参数。②对像对进行干涉处理并利用DRTMV4作为外部DEM数据去除干涉对地形效应。③估算和去除残余的恒定相位和解缠^［34］后还存在的相位坡道^［35］。④若解缠结果不理想，会对干涉图做二次解缠来优化。其中，线性模型最为稳定。⑤大气延迟是InSAR中主要误差源。利用定制的高通滤波和低通滤波去除大气相位。⑥将所得SBAS结果进行编码并对位移速率图进行分析。图3为基于SBAS-InSAR技术的LOS方向2018年11月15日至2022年6月3日间地表平均形变速率图^［36］，黑色椭圆区域为GNSS布置区。图4为同期的地表位移时间序列图。

图3

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

图3 基于SBAS-InSAR技术的LOS方向时间序列形变速率场（2018年11月15日至2022年6月3日）

Fig.3 Deformation rate field of LOS direction time series based on SBAS-InSAR technology （November 15， 2018 to June 3， 2022）

图4

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

图4 基于SBAS-InSAR技术的LOS方向时间序列形变场（2018 年 11 月 15 日至 2022 年 6月3 日）

Fig.4 Deformation field of LOS direction time series based on SBAS-InSAR technology（November 15， 2018 to June 3， 2022）

如图3、图4所示，南山研究区山体陡峭，南山隧道正位于山脚之下，山体滑坡将会对南山隧道造成严重威胁。研究区内最大形变速率-41.6 mm/a，最大位移形变-75.6 mm，威胁着南山隧道周边建筑和人民生命财产安全。

4.1　拟合推估法分析

InSAR处理生成的矢量点带有经纬度坐标信息，根据GNSS点的经纬度坐标和envi定位工具，找到对应位置的矢量点。为比较GNSS和InSAR的监测结果，利用公式（1）将GNSS测得方向转换成LOS向，如图5所示。黑色曲线表示InSAR的LOS向累积形变值，红色曲线表示GNSS投影到LOS方向的累积形变值。表2给出了各点绝对误差和均方根误差统计结果。

图5

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

图5 2018 年 11 月 15 日至 2022 年 6月3 日监测点视线向位移变化曲线

Fig.5 Change curve of LOS direction displacement at monitoring points from November 15， 2018 to June 3， 2022

表2 InSAR监测结果精度统计/mm

Table 2 Accuracy statistics of InSAR monitoring results/mm

点号	最大绝对误差	最小绝对误差	平均绝对误差	均方根误差
Jc01	14.7	1.9	7.6	8.2
Jc02	10.0	0.1	3.4	5.0
Jc03	9.5	0.0	0.3	4.0
Jc05	19.6	0.2	11.1	12.4
Jc06	20.0	0.2	9.4	12.0
Jc08	8.4	0.1	1.3	4.4
Jc09	13.3	0.5	2.6	4.7
Jc10	10.3	0.0	2.6	4.1
Jc11	10.7	0.2	5.2	6.5
Jc12	9.5	0.2	5.5	6.3
Jc13	9.5	0.4	1.5	4.2
Jc16	11.0	0.0	7.3	8.3
Jc18	18.4	0.0	10.0	11.6

新窗口打开| 下载CSV

由图5可知，GNSS视线向形变量呈周期性变化，InSAR视线向形变量周期性较弱。InSAR和GNSS监测曲线的方向大致相同^［37］，但仍存在一些点位的InSAR和GNSS监测结果有明显差距，如5号、6号和18号点。这些点位可能是由于SAR卫星轨道误差、大气延迟等误差影响，导致与同期GNSS监测结果有差距。为了对InSAR观测值进行校正，在研究区同期观测的13个GNSS站中均匀选择7个GNSS点位作为先验点参与模型解算，剩余6个点位参与精度评定。如图6所示，黑色曲线表示InSAR视线向累积形变值，红色曲线表示投影到LOS方向的GNSS累积形变值，蓝色曲线表示经校正后InSAR累积形变值。

图6

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

图6 2018 年 11 月 15 日至 2022 年 6月3 日6个站点InSAR观测值、InSAR校正值和GNSS视线向观测值值趋势对比图

Fig.6 Comparison Diagram of InSAR Observation Values， InSAR Correction Values and GNSS Line-of-sight Observation Values at Six Stations from November 15， 2018 to June 3， 2022

从图6可以看出，校正后的InSAR观测值与GNSS观测值差值明显减少，尤其是12号和16号，校正后的InSAR观测值曲线十分接近GNSS观测值曲线，在一定程度上也表现出明显的周期性。6号点校正结果稍差，但仍比未校正之前的结果要好，这是由于先验点的选择和数量不同，所求得的模型参数不同，无法适配所有点位，导致某些点的校正效果稍差。由表3可以看出，9号和12号校正后最大误差和均方根差不超过3 mm，校正效果最好。6号点校正后最大误差13.0 mm，均方根误差7.6 mm，也同比校正前减少了4.4 mm。6个点最大误差平均减少了3.4 mm，均方差根平均减少了3.5 mm。表4为校正前后InSAR整体监测结果精度对比表。

表3 校正后InSAR监测结果精度统计/mm

Table 3 Accuracy statistics of InSAR monitoring results after correction/mm

点号	最大绝对误差	最小绝对误差	平均绝对误差	均方根误差
Jc06	13.0	0.7	5.0	7.6
Jc09	4.6	0.0	1.6	2.5
Jc12	4.2	0.0	1.0	1.6
Jc13	4.6	0.1	2.8	3.5
Jc16	7.6	0.0	1.9	3.5
Jc18	11.0	0.7	5.7	7.3

新窗口打开| 下载CSV

表4 校正前后InSAR整体监测结果精度表/mm

Table 4 Precision table of InSAR overall monitoring results before and after correction/mm

改正前			改正后
	平均绝对误差	均方根误差	平均绝对误差	均方根误差
整体优化精度	6.4	7.8	3.0	4.3

新窗口打开| 下载CSV

表4精度检验显示，校正后的InSAR观测值更适合参与到GNSS-InSAR数据的融合解算中。校正后InSAR观测值的平均误差和均方根误差均比校正前误差小，其参与校正的6个点位，形变监测精度改善极为显著。整体平均误差校正后提高了约53%，均方根误差提高了约44%，说明使用拟合推估方法较好地抵制了InSAR 观测值中因变形不均所引起的系统误差，更有效的减弱观测值中的系统误差的影响。

4.2　卡尔曼滤波分析

为便于开展实验准确性的验证，首先将7个参与InSAR观测值校正的GNSS点除去，然后剩余6个GNSS点位的观测结果进行卡尔曼滤波的融合解算。根据采集数据的时间分布滤波计算了102期。图7所示是监测点在E、N、U方向上不同方法得到的融合观测值。其中，黑色曲线表示GNSS观测值，红色曲线表示未校正的InSAR观测值和GNSS观测值融合后的观测值，蓝色曲线表示校正后的InSAR观测值和GNSS观测值融合后的观测值。图8是本文方法在E、N、U方向拟合精度和拟合残差分布图。

图7

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

图7 E、N、U 3个方向上GNSS观测值、直接融合值和校正正后融合值对比图

Fig.7 Comparison of GNSS observed value， direct fusion value and fusion value after correction of the Remaining Stations in the Three Directions E， N and U

图8

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

图8 E、N、U方向拟合残差分布图

Fig.8 Residual distribution maps were fitted in E， N and U directions

从图7可以看出，经校正的InSAR观测值和GNSS观测值融合后在E、N、U 3个方向上的曲线更接近GNSS观测值，周期性也更明显。而未校正的InSAR观测值和GNSS观测值融合后与GNSS观测值仍有较大误差。说明拟合推估顾及到观测区内不符的整体变化趋势，可以高精度的拟合区域内的不符值从而削弱 InSAR 观测值中的系统误差的影响。同时，GNSS-InSAR 数据通过卡尔曼滤波融合后，充分地利用InSAR、GNSS观测数据，动态地反映地表形变特点，得到更精确的地表三维形变结果。

从图8、图9可以看出，E、N、U 3个方向拟合残差最稳定且最小的是N方向，残差变化最大的是U方向上。E方向上的拟合精度初始值最大，但随着卡尔曼滤波期数的增加，E方向上的拟合残差逐渐减小，这是由于卡尔曼滤波的初始状态往往难以确定，利用经验给定，导致误差较大。但是随着滤波进行，之后各期拟合残差较之前变小。U、N方向上的拟合残差从开始到结束基本趋于稳定。

图9

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

图9 E、N、U方向拟合精度

Fig.9 E， N， U direction fitting accuracy

为说明在GNSS-InSAR融合前对InSAR进行校正的必要性，分别用已校正的InSAR观测值和未校正的InSAR观测值与GNSS观测值融合，图10为InSAR校正前后与GNSS融合的均方根误差统计结果。表5为两者均方根误差统计结果。

图10

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

图10 InSAR校正前后与GNSS融合在E、N、U方向上均方根误差比较（单位：mm）

Fig.10 Comparison of root mean square error between InSAR and GNSS fusion in E， N and U directions before and after calibration（Units：mm）

表5 整体均方根误差比较（LOS向）/mm

Table 3 Comparison of global root mean square error (LOS direction) /mm

点位	Klman滤波InSAR（未校正）-GNSS融合值	Klman滤波InSAR（校正后）-GNSS融合值
jc06	9.4	6.1
jc09	3.2	1.8
jc12	4.7	1.2
jc13	3.4	2.8
jc16	6.7	2.4
jc18	9.4	5.9

新窗口打开| 下载CSV

从图10可以看出，未校正的InSAR观测值和GNSS融合后，N和U方向均方根差很大，尤其是N方向，说明未校正InSAR数据存在着一系列的系统误差。将InSAR观测值校正后，N、U方向上的均方根差明显降低。从表5可以看出，误差较大的6号点在校正前后均方根差变化最为明显，其余点的均方根差也有明显变化，这说明InSAR观测值中系统误差的存在会大大降低GNSS-InSAR数据融合解算的精度，而InSAR观测值经过系统误差的校正后，融合结果的精度也会显著增加。

基于InSAR校正值和卡尔曼滤波的GNSS-InSAR融合方法得到的E、N、U方向上时序三维形变场如图11所示。

图11

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

图11 E、N、U时序三维形变场

Fig.11 3D deformation fields of E， N and U time series

从图11得到的西宁南山E方向时序三维形变场可以看出，南山西南部地区有向东运动趋势，从各期形变速率大小可以看出这种运动趋势相对稳定。从N方向时序三维形变场可以看出，南山西南部地区有明显向南运动趋势，且形变速率较快。从U方向时序三维形变场可以看出，南山中部区域有下沉趋势，从各期形变速率大小可以看出这种运动趋势相对稳定，南山东部有上升趋势。整体上看，南山西南区域呈东南运动趋势，中部地区下沉。

5 总结

常规GNSS、InSAR数据在联合解算过程中，无法充分利用GNSS、InSAR数据之间强互补性，降低了两者的时间分辨率。在融合算法中较多的考虑到算法本身，忽略了GNSS、InSAR数据本身的系统误差，导致精度结果无法达到最优。本文通过对拟合推估和卡尔曼滤波模型的研究，设计了拟合推估法来校正InSAR数据中的存在的系统误差。利用卡尔曼滤波将系统校正后的InSAR观测值和GNSS观测值融合解算，得到了精度更高的地表三维形变结果。

以西宁地区为例，利用西宁的103景sentinel-1A影像和同期GNSS监测数据，将GNSS观测数据和InSAR数据进行融合研究。实验结果表明：通过拟合推估法校正的InSAR观测值比未校正的InSAR观测值整体平均误差提高了大约53%，均方根误差提高了约44%。校正后的InSAR观测值与GNSS观测值经卡尔曼滤波融合结果比未校正的InSAR观测值与GNSS观测值融合结果精度高45%，比InSAR观测值精度高57%。本文方法先对InSAR数据进行了校正，利用卡尔曼滤波对两种数据进行融合，在充分利用数据的前提下，动态地反映地表形变特点，提高了形变精度，对国家防灾减灾工作具有现实意义。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

WEI

Gang

， YIN

Zhiqiang

， SHI

Liqun

，et al.

Geological characteristics and stability analysis of linjiaya landslide in Beishan area， Xining City

［J］. Geology and Resources，2015，24（2）：146-151.