CN101770027B

CN101770027B - Surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion

Info

Publication number: CN101770027B
Application number: CN2010101067941A
Authority: CN
Inventors: 何秀凤; 何敏; 罗海滨; 黄其欢
Original assignee: Hohai University HHU
Current assignee: Hohai University HHU
Priority date: 2010-02-05
Filing date: 2010-02-05
Publication date: 2012-05-16
Anticipated expiration: 2030-02-05
Also published as: CN101770027A

Abstract

The invention discloses a surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion, which comprises the following steps: 1. laying GPS points; 2. collecting GPS signals and SAR images; 3. resolving a GPS point; 4. inhibiting speckle noise of the SAR image; 5. establishing a coordinate conversion relation between GPS and SAR data; 6. correcting SAR satellite orbit errors by utilizing GPS data; 7. forming an InSAR differential interference pattern; 8. inverting atmospheric delay by using a GPS point; 9. correcting InSAR atmospheric delay errors; 10. unwrapping the InSAR differential interferogram; 11. converting the unwrapping phase into deformation information; InSAR deformation information geocoding; 13. fusing InSAR and GPS deformation information to obtain a high-precision surface three-dimensional deformation result; 14. and estimating and obtaining a surface three-dimensional deformation result with high space-time resolution by using Kalman filtering. The invention adopts the earth surface three-dimensional deformation monitoring technology of the InSAR and GPS data fusion, thereby breaking through the application limitation of a single monitoring technology and greatly improving the space-time resolution of the three-dimensional deformation monitoring result. The monitoring precision is high, especially the precision in the vertical direction, and meanwhile the overall reliability of the system is enhanced.

Description

3D Surface Deformation Monitoring Method Based on InSAR and GPS Data Fusion

技术领域 technical field

本发明涉及一种空间对地监测地表三维形变的方法，具体说是一种基于InSAR与GPS数据融合的地表三维形变监测方法，适用于城市环境、滑坡和地震等区域高时空分辨率、高精度的地表三维形变监测。The invention relates to a method for space-to-ground monitoring of three-dimensional deformation of the surface, specifically a method for monitoring three-dimensional deformation of the surface based on the fusion of InSAR and GPS data, which is suitable for urban environments, landslides, earthquakes and other areas with high temporal and spatial resolution and high precision 3D deformation monitoring of the ground surface.

背景技术 Background technique

近年来，随着卫星定位理论的日臻成熟以及接收设备性能的提高，全球定位系统(GPS)已广泛应用于形变监测。GPS可以提供高时间分辨率、高精度的三维形变测量值，而且不受天气的影响、点位间毋需通视、布站灵活、对地表状况适应能力强。但由于接收机数量和布网阵列限制，GPS空间分辨率较低、覆盖区域小，可能遗漏一些重大安全隐患。具体说来，现有GPS形变监测方法主要存在以下不足：(1)采用布网观测的方法，劳动强度大、监测系统成本高；(2)采用点观测的方法，因而监测数据的空间分辨率低，难以满足大区域变形分析和地质灾害预测的要求；(3)在深山峡谷或城市建筑物密集区，GPS接收机天线受到遮挡，这使得接收到的GPS卫星数减少，导致GPS定位精度大大下降，达不到形变观测精度要求。In recent years, with the maturation of satellite positioning theory and the improvement of the performance of receiving equipment, Global Positioning System (GPS) has been widely used in deformation monitoring. GPS can provide high-time resolution, high-precision three-dimensional deformation measurement value, and it is not affected by the weather, there is no need for direct vision between points, flexible station layout, and strong adaptability to surface conditions. However, due to the limitation of the number of receivers and network arrays, the spatial resolution of GPS is low and the coverage area is small, which may miss some major security risks. Specifically, the existing GPS deformation monitoring methods mainly have the following deficiencies: (1) The method of network deployment observation is labor-intensive and the cost of the monitoring system is high; (2) The method of point observation is used, so the spatial resolution of monitoring data is limited. It is difficult to meet the requirements of deformation analysis and geological disaster prediction in large areas; (3) In deep mountains and canyons or densely built urban areas, the GPS receiver antenna is blocked, which reduces the number of received GPS satellites, resulting in a significant increase in GPS positioning accuracy Decrease, can not meet the deformation observation accuracy requirements.

合成孔径雷达干涉测量(InSAR)是近年来迅速发展起来的地表探测新技术，已有越来越多的国家和地区利用InSAR的差分技术来探测诸如由矿山开采、地震、火山运动等引起的地面形变现象。与GPS技术相比，InSAR技术具有垂直监测精度高、监测范围大、空间近连续性等优点。但InSAR技术时间分辨率低且对大气延迟误差、卫星轨道误差、地表状况以及时态不相干等因素非常敏感，这造成了InSAR技术在地表形变探测应用中的困难。Synthetic Aperture Radar Interferometry (InSAR) is a new technology for surface detection that has developed rapidly in recent years. More and more countries and regions have used InSAR differential technology to detect ground surface areas caused by mining, earthquakes, volcanic movements, etc. Deformation phenomenon. Compared with GPS technology, InSAR technology has the advantages of high vertical monitoring accuracy, large monitoring range, and near spatial continuity. However, the time resolution of InSAR technology is low and it is very sensitive to factors such as atmospheric delay error, satellite orbit error, surface conditions, and temporal irrelevance, which makes it difficult for InSAR technology to be used in surface deformation detection applications.

发明内容 Contents of the invention

本发明的目在于，克服现有技术存在的缺陷，综合利用InSAR与GPS技术，突破单一技术的应用局限，提供一种基于InSAR与GPS数据融合的地表三维形变监测方法。The purpose of the present invention is to overcome the defects of the existing technology, comprehensively utilize InSAR and GPS technology, break through the application limitation of a single technology, and provide a three-dimensional surface deformation monitoring method based on InSAR and GPS data fusion.

本发明的基本原理为：Basic principle of the present invention is:

InSAR技术是近年来迅速发展起来的地表探测新技术，其原理是通过处理接收到的不同时相地表回波来监测地表的形变。与GPS技术相比，InSAR技术具有垂直监测精度高、监测范围大、空间近连续性等优点。但InSAR技术时间分辨率低且对大气延迟误差、卫星轨道误差、地表状况以及时态不相干等因素非常敏感，这造成了InSAR技术在地表形变探测应用中的困难。InSAR technology is a new surface detection technology developed rapidly in recent years. Its principle is to monitor the deformation of the surface by processing the received surface echoes in different phases. Compared with GPS technology, InSAR technology has the advantages of high vertical monitoring accuracy, large monitoring range, and near spatial continuity. However, the time resolution of InSAR technology is low and it is very sensitive to factors such as atmospheric delay error, satellite orbit error, surface conditions, and temporal irrelevance, which makes it difficult for InSAR technology to be used in surface deformation detection applications.

本发明考虑将GPS与InSAR数据融合起来，突破单一技术应用的局限，发挥其各自优势，极大地改善空间域和时间域的分辨能力。但是将GPS与InSAR的数据融合存在以下需要解决关键技术：(1)利用GPS改正InSAR卫星轨道误差；(2)利用GPS改正InSAR数据的大气延迟影响；(3)GPS辅助InSAR相位解缠；(4)InSAR与GPS数据在时空域融合。The present invention considers the fusion of GPS and InSAR data, breaks through the limitations of a single technology application, exerts their respective advantages, and greatly improves the resolution capabilities of the space domain and the time domain. However, the following key technologies need to be solved in the fusion of GPS and InSAR data: (1) use GPS to correct InSAR satellite orbit error; (2) use GPS to correct the influence of atmospheric delay of InSAR data; (3) GPS-assisted InSAR phase unwrapping; ( 4) InSAR and GPS data fusion in time and space domain.

考虑到要发明一种具有高精度的监测方法，就要对观测数据的误差进行改正，为此根据InSAR数据处理的流程，将GPS数据融合其中，利用GPS改正的InSAR误差，然后采用基于改进马尔科夫随机场的GPS和InSAR融合方法，实现InSAR与GPS数据在时空域的融合，获取高时空分辨率的地表三维形变信息。Considering that in order to invent a high-precision monitoring method, it is necessary to correct the error of the observation data. Therefore, according to the flow of InSAR data processing, the GPS data is integrated into it, and the InSAR error corrected by GPS is used. The GPS and InSAR fusion method of the Cove random field realizes the fusion of InSAR and GPS data in the space-time domain, and obtains the three-dimensional deformation information of the surface with high space-time resolution.

1、GPS改正InSAR数据的轨道误差1. GPS corrects the orbit error of InSAR data

在InSAR图像配准、相位到高程的转换和地理编码等步骤都要用到SAR卫星的轨道信息。可见，获取精确的轨道信息对提高InSAR测量的精度至关重要。由于InSAR数据自身带的轨道信息精度不是很高，如，ERS的基本轨道数据的精度在飞行方向2-4m，垂直轨迹方向1-2m，径向5m，不能满足测量精度的要求，为此，本发明引入GPS观测数据，利用GPS精确测得InSAR图像上特征地物(如河流和道路的交叉口，水塔等)或人工角反射器的地面坐标，再利用地理编码的逆过程求得SAR卫星精确轨道坐标，达到改正InSAR数据的轨道误差的目的。The orbital information of SAR satellites is used in the steps of InSAR image registration, phase-to-elevation conversion, and geocoding. It can be seen that obtaining accurate orbital information is crucial to improving the accuracy of InSAR measurement. Because the accuracy of the orbit information carried by InSAR data itself is not very high, for example, the accuracy of the basic orbit data of ERS is 2-4m in the flight direction, 1-2m in the vertical trajectory direction, and 5m in the radial direction, which cannot meet the measurement accuracy requirements. Therefore, The present invention introduces GPS observation data, utilizes GPS to accurately measure the ground coordinates of characteristic objects (such as intersections of rivers and roads, water towers, etc.) or artificial corner reflectors on InSAR images, and then uses the inverse process of geographic coding to obtain SAR satellites Accurate orbit coordinates to achieve the purpose of correcting the orbit error of InSAR data.

2、GPS反演大气延迟改正InSAR大气延迟误差2. GPS inversion atmospheric delay corrects InSAR atmospheric delay error

SAR传感器发射的电磁波在穿过大气时，传播速度将发生变化，传播路径也将发生弯曲，从而产生延迟现象。同一区域不同时间观测得到的SAR影像的大气延迟不完全相同，因而不可避免的对InSAR干涉相位图产生污染。为此，本发明利用GPS反演大气延迟，首先采用随机过程法获取GPS的对流层估计，其次采用站间时域双差法估计大气延迟改正，考虑到大气延迟与高程具有较强的相关性，根据地理空间相关性理论，提出采用融合GPS高程信息的改进反距离加权内插法对其加密，然后将加密后的大气延迟改正从InSAR干涉相位图中消除，从而改正InSAR大气延迟误差。When the electromagnetic wave emitted by the SAR sensor passes through the atmosphere, the propagation speed will change, and the propagation path will also bend, resulting in delay phenomenon. The atmospheric delays of SAR images observed at different times in the same area are not exactly the same, which inevitably pollutes the InSAR interferogram. For this reason, the present invention utilizes GPS inversion atmospheric delay, first adopts stochastic process method to obtain the tropospheric estimate of GPS, secondly adopts inter-station time domain double-difference method to estimate atmospheric delay correction, considering that atmospheric delay has a strong correlation with elevation, According to the geospatial correlation theory, an improved inverse distance weighted interpolation method fused with GPS elevation information is proposed to encrypt it, and then the encrypted atmospheric delay correction is eliminated from the InSAR interferogram to correct the InSAR atmospheric delay error.

3、GPS辅助InSAR相位解缠3. GPS-assisted InSAR phase unwrapping

为了将相位转换为高程或形变信息，必须对干涉图中主值区间内的缠绕相位进行恢复得到其绝对相位，即要进行相位的解缠。残差点的存在使得解缠相位不连续甚至出现无法解缠的“孤岛”。为此，本发明引入少量的GPS控制点，以这些GPS控制点作为新的解缠起算点，消除“孤岛”并使整个区域的解缠相位连续。In order to convert the phase into elevation or deformation information, it is necessary to restore the wrapped phase in the principal value interval of the interferogram to obtain its absolute phase, that is, to unwrap the phase. The existence of residual points makes the unwrapped phases discontinuous and even "islands" that cannot be unwrapped appear. For this reason, the present invention introduces a small number of GPS control points, and uses these GPS control points as new unwrapping starting points to eliminate "islands" and make the unwrapping phases of the whole area continuous.

本发明根据解缠算法误差传播的特点，提出了两种不同的GPS辅助InSAR相位解缠的方法，即，GPS辅助InSAR Goldstein枝切线解缠算法和引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法。According to the characteristics of unwrapping algorithm error propagation, the present invention proposes two different GPS-assisted InSAR phase unwrapping methods, that is, GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm and GPS-assisted phase unwrapping algorithm based on Markov random field that introduces residual points. InSAR unwrapping algorithm.

4、基于改进马尔科夫随机场的InSAR和GPS数据融合方法4. InSAR and GPS data fusion method based on improved Markov random field

对于像滑坡、地震这样的地质灾害，会在短时间内引起大范围较大的形变，InSAR一个月左右的重复观测周期是不够的，需要在时间域内插值估计。另外，GPS是基于有限点的观测，点的空间密度不能满足监测的需求，需要在空间域内插值估计。本发明采用基于改进马尔科夫随机场的InSAR和GPS数据融合方法，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息，获得高空间分辨率的三维形变场；然后对高时间频率GPS数据(其所采集的地面沉降量)在时间域上对上述高空间分辨率的三维形变场进行克里格插值和加密，从而将准GPS结果在时间域内加密成准GPS时间序列，获取高时间分辨率的三维形变场；最后利用卡尔曼滤波对格网中的所有点进行估计获得高时空分辨率的地表三维形变结果，实现地表形变高时空分辨率的监测。For geological disasters such as landslides and earthquakes, which will cause large-scale and large deformations in a short period of time, InSAR's repeated observation cycle of about one month is not enough, and interpolation estimation in the time domain is required. In addition, GPS is based on the observation of limited points, and the spatial density of points cannot meet the needs of monitoring, so it needs to be estimated by interpolation in the spatial domain. The present invention adopts the InSAR and GPS data fusion method based on the improved Markov random field, decomposes the one-dimensional deformation result of the InSAR line of sight into three-dimensional deformation information through GPS constraints, and obtains a three-dimensional deformation field with high spatial resolution; The frequency GPS data (the land subsidence collected by it) performs Kriging interpolation and encryption on the above-mentioned three-dimensional deformation field with high spatial resolution in the time domain, so that the quasi-GPS results are encrypted into a quasi-GPS time series in the time domain, and the obtained Three-dimensional deformation field with high temporal resolution; finally, Kalman filter is used to estimate all points in the grid to obtain three-dimensional surface deformation results with high temporal and spatial resolution, so as to realize the monitoring of surface deformation with high temporal and spatial resolution.

本发明基于InSAR与GPS数据融合的地表三维形变监测方法，其步骤是：The present invention is based on the surface three-dimensional deformation monitoring method of InSAR and GPS data fusion, and its steps are:

1、布设GPS点1. Set up GPS points

布设的GPS点，既要作为控制点，用于GPS和InSAR数据坐标转化和改正InSAR误差改正，又要作为监测点，对测区的重点部位进行监测。因此，布设的GPS点分为两类：第一类是监测区内稳定的特征点(称为第一类GPS点)(稳定的特征点的稳定是指不受形变影响)，如河流、道路的交叉处，这类点主要用于将GPS数据和InSAR数据统一到同一参考坐标系下、改正SAR卫星轨道误差以及反演大气延迟；第二类是监测区的重点监测部位(称为第二类GPS点)。The deployed GPS points should not only be used as control points for GPS and InSAR data coordinate transformation and correction of InSAR error correction, but also as monitoring points to monitor key parts of the survey area. Therefore, the deployed GPS points are divided into two categories: the first type is the stable feature points in the monitoring area (called the first type of GPS points) (the stability of stable feature points means that they are not affected by deformation), such as rivers, roads, etc. This type of point is mainly used to unify GPS data and InSAR data into the same reference coordinate system, correct SAR satellite orbit error and invert atmospheric delay; the second type is the key monitoring part of the monitoring area (called the second like GPS points).

在布设的第一类GPS点时，要注意点的分布，尽量使其均匀分布在测区内。When laying out the first type of GPS points, pay attention to the distribution of points and try to make them evenly distributed in the survey area.

2、采集并记录GPS信号和收集SAR图像2. Collect and record GPS signals and collect SAR images

因为要用GPS数据改正InSAR数据的轨道误差和大气延迟误差，采集的GPS信号要与SAR图像获取的时间同步，GPS信息的采集时间段最好能向前和向后顺延10分钟左右。Because GPS data is used to correct the orbit error and atmospheric delay error of InSAR data, the collected GPS signal must be synchronized with the time of SAR image acquisition, and the GPS information collection time period should preferably be extended forward and backward by about 10 minutes.

3、依据GPS信号解算出两类GPS点的三维坐标信息3. Calculate the three-dimensional coordinate information of two types of GPS points according to the GPS signal

两类GPS点分开解算，解算的第一类GPS点的三维坐标信息包括：平面位置信息和高程信息，第二类GPS点的三维坐标信息包括：平面位移和垂直沉降。The two types of GPS points are calculated separately. The three-dimensional coordinate information of the first type of GPS point includes: plane position information and elevation information, and the three-dimensional coordinate information of the second type of GPS point includes: plane displacement and vertical settlement.

4、SAR图像的斑点噪声抑制4. Speckle noise suppression for SAR images

有效地抑制SAR图像的斑点噪声，有利于提高InSAR监测结果精度和确保监测结果的可靠性。本发明采用改进的Lee算法(是现有的算法)。Effective suppression of speckle noise in SAR images is beneficial to improving the accuracy of InSAR monitoring results and ensuring the reliability of monitoring results. The present invention adopts the improved Lee algorithm (an existing algorithm).

5、搜索SAR图像上的特征点，并建立搜索到的特征点与第一类GPS点之间的坐标转换关系5. Search for the feature points on the SAR image, and establish the coordinate transformation relationship between the searched feature points and the first type of GPS points

为了将InSAR和GPS的数据进行融合，必须将它们的数据统一到同一参考坐标系下。搜索InSAR图像上特征地物，并用GPS测定其地面坐标。假设特征点目标的影像坐标为(I_row，I_col)，其地面坐标为(G_lat，G_lan)。在统一投影关系下，两套坐标系统的变换可表达为：In order to fuse InSAR and GPS data, it is necessary to unify their data into the same reference coordinate system. Search for feature objects on InSAR images, and use GPS to measure their ground coordinates. Assume that the image coordinates of the feature point target are (I _row , I _col ), and their ground coordinates are (G _lat , G _lan ). Under the unified projection relationship, the transformation of the two sets of coordinate systems can be expressed as:

$[\begin{matrix} {G G}_{lat lat} & {G G}_{lon the lon} \end{matrix}] = = [\begin{matrix} {I I}_{row row} & {I I}_{col col} & 11 \end{matrix}] [\begin{matrix} {a a}_{11} & {a a}_{22} \\ {b b}_{11} & {b b}_{22} \\ {c c}_{11} & {c c}_{22} \end{matrix}] - - - - - - ((11))$

或or

L＝BX (2)L=BX (2)

其中，L为特征点目标的GPS坐标矩阵，B为其对应影像坐标矩阵，x是3×2阶的转换矩阵。因此，至少需要3个特征点数据的数据才可以确定转换矩阵中的6个参数。当有更多的特征点数据时，可用最小二乘平差求得6个转换参数。假设观测误差为

，则观测方程为：Among them, L is the GPS coordinate matrix of the feature point target, B is the corresponding image coordinate matrix, and x is a transformation matrix of order 3×2. Therefore, at least 3 feature point data are required to determine the 6 parameters in the transformation matrix. When there are more feature point data, the least squares adjustment can be used to obtain 6 conversion parameters. Suppose the observation error is

, then the observation equation is:

L＝BX+□ (3)L＝BX+□ (3)

转换矩阵的估计值为：The estimated value of the transformation matrix is:

$\overset{^^}{X x} = = {(({B B}^{T T} PB PB))}^{- - 11} {B B}^{T T} PL PL - - - - - - ((44))$

其中，，D_□是

的方差矩阵。最后，利用转换矩阵将GPS求得的大气延迟改正、轨道改正等数据平移、旋转转换为InSAR影像坐标系下的数据。in, , _D is

The variance matrix of . Finally, the transformation matrix is used to translate and rotate the atmospheric delay correction, orbit correction and other data obtained by GPS into the data in the InSAR image coordinate system.

6、利用GPS数据改正SAR卫星的轨道误差6. Using GPS data to correct the orbit error of SAR satellites

采用测图方程法和最小二乘估计相结合的方法，具体是：A method combining surveying equation method and least squares estimation is adopted, specifically:

在轨道误差的改正中，噪声抑制后的SAR图像作为输入数据，搜索InSAR图像上稳定的特征地物，并用GPS测定其地面坐标。根据特征地物的雷达坐标(即特征地物到InSAR图像的斜距和视角)与其地面坐标的转换关系，构建地形测图方程，即In the correction of orbit error, the noise-suppressed SAR image is used as input data, and the stable feature objects on the InSAR image are searched, and their ground coordinates are measured by GPS. According to the conversion relationship between the radar coordinates of the feature features (that is, the slant distance and viewing angle from the feature features to the InSAR image) and its ground coordinates, the topographic mapping equation is constructed, namely

$[\begin{matrix} X x \\ Y Y \\ Z Z \end{matrix}] = = {λ λ}_{00} [\begin{matrix} {a a}_{11} & {a a}_{22} & {a a}_{33} \\ {b b}_{11} & {b b}_{22} & {b b}_{33} \\ {c c}_{11} & {c c}_{22} & {c c}_{33} \end{matrix}] [\begin{matrix} 00 \\ r r sin sin θ θ \\ r r cos cos θ θ \end{matrix}] + + [\begin{matrix} {X x}_{s the s} \\ {Y Y}_{s the s} \\ {Z Z}_{s the s} \end{matrix}] - - - - - - ((55))$

其中，X，Y，Z为特征点的地面坐标，a₁，a₂，…，c₃为SAR卫星姿态角

构成的旋转矩阵元素，λ₀是比例因子，r为征地物到SAR图像的斜距，θ为视角，X_s，Y_s，Z_s为SAR卫星的轨道位置。Among them, X, Y, Z are the ground coordinates of the feature points, a ₁ , a ₂ ,..., c ₃ are the attitude angles of the SAR satellite

λ ₀ is the scale factor, r is the slant distance from the land requisition object to the SAR image, θ is the angle of view, X _s , Y _s , and Z _s are the orbital positions of the SAR satellite.

然后，将SAR卫星的轨道参数和基线参数表示成为时间的线性函数，即：Then, the orbital parameters and baseline parameters of the SAR satellite are expressed as a linear function of time, namely:

X_s＝X_s0+X_s(t-t₀)X _s ＝X _s0 +X _s (tt ₀ )

. (6).

B＝B₀+B(t-t₀)B＝B ₀ +B(tt ₀ )

对地形测图方程线性化，可得观测误差方程：By linearizing the topographic mapping equation, the observation error equation can be obtained:

V＝AX-L (7)V＝AX-L (7)

式中，In the formula,

V＝[V_X，V_Y，V_Z]^T V=[V _X , V _Y , V _Z ] ^T

$A A = = [\begin{matrix} {a a}_{1111} & {a a}_{1212} & . . . . . . & {a a}_{117117} \\ {a a}_{21 twenty one} & {a a}_{22 twenty two} & . . . . . . & {a a}_{217217} \\ {a a}_{3131} & {a a}_{3232} & . . . . . . & {a a}_{317317} \end{matrix}]$

L＝[L_X，L_Y，L_Z]L=[L _X , L _Y , L _Z ]

其中，ΔX_s0，ΔY_s0，ΔZ_s0和

分别为t₀时刻的轨道位置和姿态改正数；

和

分别为相应的变化率改正数；Δλ₀，ΔB₀分别为t₀时刻的比例尺和基线改正数；

分别为相应的变化率改正数；Δφ₀为相位常数改正数；a₁₁，a₁₂，…a₃₁₇为各未知参数的偏导数；V_X，V_Y，V_Z分别为X，Y，Z的已知值与其近似值的差。Among them, ΔX _s0 , ΔY _s0 , ΔZ _s0 and

are the orbital position and attitude correction number at time t ₀ , respectively;

and

are the corresponding rate-of-change corrections; Δλ ₀ and ΔB ₀ are the scale and baseline corrections at t ₀ , respectively;

are the correction numbers of the corresponding rate of change; Δφ ₀ is the correction number of the phase constant; a ₁₁ , a ₁₂ ,...a ₃₁₇ are the partial derivatives of each unknown parameter; V _X , V _Y , V _Z are X, Y, Z respectively The difference between a known value and its approximate value.

由于误差方程含有17个未知数，因此至少需要6个以上的第一类GPS点解算这些未知参数。本发明采用最小二乘法估计观测误差方程中的17个未知参数，由此得到SAR卫星轨道数据误差的改正值。Since the error equation contains 17 unknowns, at least 6 GPS points of the first type are needed to solve these unknown parameters. The invention adopts the least square method to estimate 17 unknown parameters in the observation error equation, thereby obtaining the correction value of the SAR satellite orbit data error.

7、形成InSAR差分干涉图7. Form InSAR differential interferogram

首先对主副SAR图像进行精配准，配准精度要达到1/8像元；其次将副图像重采样到主图像的坐标系中；再次最后将配准后的主副图像作复共轭相乘，生成干涉图；然后消除平地效应相位；最后采用两轨法或三轨法进行差分干涉处理，获得差分干涉图。Firstly, the primary and secondary SAR images are finely registered, and the registration accuracy should reach 1/8 pixel; secondly, the secondary image is resampled to the coordinate system of the main image; and finally, the registered primary and secondary images are complex conjugated Multiply to generate an interferogram; then eliminate the flat-earth effect phase; finally use the two-track or three-track method for differential interferometry to obtain a differential interferogram.

8、利用第一类GPS点反演大气延迟，并将其统一到InSAR的坐标系统中8. Use the first type of GPS points to invert the atmospheric delay and unify it into the InSAR coordinate system

利用GPS反演大气延迟，首先采用随机过程法获取GPS的对流层估计，其次采用站间时域双差法估计大气延迟改正，并将其统一到InSAR的坐标系统中。其中，站间时域双差法为：Using GPS to invert the atmospheric delay, the stochastic process method is used to obtain the GPS tropospheric estimate first, and then the inter-station time-domain double-difference method is used to estimate the atmospheric delay correction, which is unified into the InSAR coordinate system. Among them, the inter-station time domain double difference method is:

假设A点是SAR图像上的一个GPS点，B是同一幅SAR图像上的另外一GPS点，如果将GPS解算得到的A、B两点在历元j的对流层延迟分别记为D_A ^j、D_B ^j，则站间对流层延迟的差值为Assuming that point A is a GPS point on the SAR image, and B is another GPS point on the same SAR image, if the tropospheric delays of points A and B at epoch j obtained by GPS solution are recorded as D _A ^j , D _B ^j , then the difference in tropospheric delay between stations is

${D D.}_{AB AB}^{j j} = = {D D.}_{B B}^{j j} - - {D D.}_{A A}^{j j} - - - - - - ((88))$

以A为参考站，利用式(8)同理可推得其它GPS点的对流层延迟站间单差。Taking A as the reference station, the tropospheric delay inter-station single difference of other GPS points can be deduced similarly by using formula (8).

假设两个站点A、B，在两个历元j(主图像获取时刻)、k(副图像获取时刻)由式(8)建立的两个单差分别为：Assuming two sites A and B, the two single differences established by formula (8) in two epochs j (primary image acquisition time) and k (secondary image acquisition time) are respectively:

${D D.}_{AB AB}^{j j} = = {D D.}_{B B}^{j j} - - {D D.}_{A A}^{j j} - - - - - - ((99))$

${D D.}_{AB AB}^{k k} = = {D D.}_{B B}^{k k} - - {D D.}_{A A}^{k k}$

由式(9)的两个单差可以得到历元双差为：From the two single differences in formula (9), the epoch double difference can be obtained as:

${D D.}_{AB AB}^{jk jk} = = {D D.}_{AB AB}^{k k} - - {D D.}_{AB AB}^{j j} = = (({D D.}_{B B}^{k k} - - {D D.}_{A A}^{k k})) - - (({D D.}_{B B}^{j j} - - {D D.}_{A A}^{j j})) - - - - - - ((1010))$

$= = (({D D.}_{B B}^{k k} - - {D D.}_{B B}^{j j})) - - (({D D.}_{A A}^{k k} - - {D D.}_{A A}^{j j}))$

由式(10)可以看出，估计大气改正的双差有两种形式：一种是站间时域双差，即先站间差然后再历元间差，另一种是时域站间双差，即先历元间差再站间差。因为站间差仅与单幅SAR图像有关，从而可以更自由、灵活的与其它具有站间差的SAR图像组合形成历元间差，所以，本发明采用更有利于实际应用的站间时域双差。It can be seen from formula (10) that there are two forms of double difference in estimated atmospheric correction: one is inter-station time-domain double-difference, that is, inter-station difference first and then inter-epoch difference, and the other is time-domain inter-station difference Double difference, that is, first the difference between epochs and then the difference between stations. Because the inter-station difference is only related to a single SAR image, it can be more freely and flexibly combined with other SAR images with inter-station differences to form an inter-epoch difference. Therefore, the present invention uses an inter-station time domain that is more conducive to practical applications. double difference.

9、利用GPS反演的大气延迟改正InSAR差分干涉图大气延迟误差9. Using the atmospheric delay retrieved by GPS to correct the atmospheric delay error of the InSAR differential interferogram

考虑到大气延迟与高程具有较强的相关性，根据地理空间相关性理论，引入高程影响因子，构建空间异相关模型，提出采用融合GPS高程信息的改进反距离加权内插法对其加密。对GPS反演的大气延迟进行加密，然后将加密后的大气延迟从InSAR干涉相位图中消除。其中，融合GPS高程信息的改进反距离加权内插法的基本原理如下：Considering the strong correlation between atmospheric delay and elevation, according to the geospatial correlation theory, the elevation influence factor is introduced to build a spatial heterogeneous correlation model, and the improved inverse distance weighted interpolation method combined with GPS elevation information is proposed to encrypt it. The atmospheric delay retrieved by GPS is encrypted, and then the encrypted atmospheric delay is eliminated from the InSAR interferogram. Among them, the basic principle of the improved inverse distance weighted interpolation method fused with GPS elevation information is as follows:

反距离加权内插法是一种常用的空间插值方法，它假设两个事物的相似程度随着彼此间距离的缩短而增加。因此，反距离加权(IDW)内插法以估值点和观测点间的距离为权重进行加权平均，距离估值点越近的观测点赋予的权重越大。反距离加权(IDW)内插法可表示为：Inverse distance weighted interpolation is a commonly used spatial interpolation method, which assumes that the similarity between two things increases as the distance between them decreases. Therefore, the inverse distance weighted (IDW) interpolation method uses the distance between the evaluation point and the observation point as the weight to carry out weighted average, and the closer the observation point is to the evaluation point, the greater the weight is given. The inverse distance weighted (IDW) interpolation method can be expressed as:

$\{\begin{matrix} Z Z (({x x}_{00},, {y the y}_{00},, {z z}_{00})) = = {Σ Σ}_{i i = = 11}^{n no} {λ λ}_{i i} Z Z (({x x}_{i i},, {y the y}_{i i},, {z z}_{i i})) \\ {λ λ}_{i i} = = \frac{{d d}_{i i 00}^{- - p p}}{{Σ Σ}_{i i = = 11}^{N N} {d d}_{i i 00}^{- - p p}} \\ {d d}_{i i 00} = = \sqrt{{(({x x}_{i i} - - {x x}_{00}))}^{22} + + {(({y the y}_{i i} - - {y the y}_{00}))}^{22} + + {(({z z}_{i i} - - {z z}_{00}))}^{22}},, ((i i = = 11 . . . . . . n no)) \end{matrix} - - - - - - ((1111))$

Z(x₀，y₀，z₀)是点(x₀，y₀，z₀)处的估计值，Z(x_i，y_i，z_i)是点(x_i，y_i，z_i)处的观测值，n是用于估计的观测值的个数，λ_i是点(x_i，y_i，z_i)处观测值对应的权，d_i0(i＝1…N)是观测点i到估值点的距离，p(p＞0)为距离的幂，且有

实际应用中，通常p＝2。Z(x ₀ , y ₀ , z ₀ ) is the estimated value at point (x ₀ , y ₀ , z ₀ ), Z(x _i , y _i , z _i ) is the point (xi , y _i _, z _i ), n is the number of observations used for estimation, λ _i is the weight corresponding to the observation at point (xi _, y _i , z _i ), d _i0 (i=1…N) is the observation The distance from point i to the evaluation point, p(p>0) is the power of the distance, and

In practical application, usually p=2.

由式(11)可以看出，反距离加权内插法将估值点和观测点间的水平距离和高程差对估值点的影响看成是等同的。但实际并非如此，如大气延迟的插值问题，估值点和观测点间的水平距离和高程差对估值点的影响并不相同。为了体现估值点和观测点间的水平距离和高程差对估值影响的不同，根据地理空间相关性理论，引入高程影响因子α，构建空间异相关模型，提出融合高程信息的改进反距离加权内插法，具体的插值公式如下：It can be seen from formula (11) that the inverse distance weighted interpolation method regards the influence of the horizontal distance and elevation difference between the evaluation point and the observation point on the evaluation point as equal. But this is not the case in reality, such as the interpolation problem of atmospheric delay, the horizontal distance and elevation difference between the estimation point and the observation point have different influences on the estimation point. In order to reflect the difference in the impact of the horizontal distance and elevation difference between the evaluation point and the observation point on the evaluation, according to the geospatial correlation theory, the elevation influence factor α is introduced, a spatial heterogeneous correlation model is constructed, and an improved inverse distance weighted method that integrates elevation information is proposed. Interpolation method, the specific interpolation formula is as follows:

$\{\begin{matrix} Z Z (({x x}_{00},, {y the y}_{00},, {z z}_{00})) = = {Σ Σ}_{i i = = 11}^{n no} {λ λ}_{i i} [[Z Z (({x x}_{i i},, {y the y}_{i i},, {z z}_{i i})) - - {m m}_{i i}]] + + {m m}_{00} \\ {m m}_{i i} = = {b b}_{00} + + {b b}_{i i} {h h}_{i i} \\ {λ λ}_{i i} = = \frac{{d d}_{i i 00}^{- - p p}}{{Σ Σ}_{i i = = 11}^{N N} {d d}_{i i 00}^{- - p p}} \\ {d d}_{i i 00} = = \frac{\sqrt{{((α α {d d}_{i i 00}^{s the s}))}^{22} + + {((((11 - - α α)) {d d}_{i i 00}^{h h}))}^{22}}}{α α} \\ {d d}_{i i 00}^{s the s} = = \sqrt{{(({x x}_{i i} - - {x x}_{00}))}^{22} + + {(({y the y}_{i i} - - {y the y}_{00}))}^{22}} \\ {d d}_{i i 00}^{h h} = = {h h}_{i i} - - {h h}_{00} ((i i = = 11,, . . . . . .,, n no)) \\ α α = = \frac{abs abs ((k k 11))}{abs abs ((k k 11)) + + abs abs ((k k 22))} \end{matrix} - - - - - - ((1212))$

其中，Z(x₀，y₀，z₀)是点(x₀，y₀，z₀)处的估计值，Z(x_i，y_i，z_i)是点(x_i，y_i，z_i)处的观测值，n是用于估计的观测值的个数，λ_i是点(x_i，y_i，z_i)处观测值对应的权，h_i为点(x_i，y_i，z_i)处的高程，h₀为点(x₀，y₀，z₀)处的高程，b₀、b₁为待求拟合系数，d_i0 ^s为观测点i到估值点的水平距离，d_i0 ^h为观测点i与估值点的高程差，k₁、k₂分别为大气延迟的差值随水平距离和高程差变化关系的线性拟合曲线的斜率，abs(·)代表取绝对值。Among them, Z(x ₀ , y ₀ , z ₀ ) is the estimated value at point (x ₀ , y ₀ , z ₀ ), Z( _xi , y _i , z _i ) is the point (xi , y _i _, z _i ), n is the number of observations used for estimation, λ _i is the weight corresponding to the observation at point (xi _, y _i , z _i ), h _i is the point ( _xi , y _i , z _i ), h ₀ is the height at point (x ₀ , y ₀ , z ₀ ), b ₀ , b ₁ are the fitting coefficients to be found, and d _i0 ^s is the distance from observation point i to evaluation point d _i0 ^h is the elevation difference between the observation point i and the evaluation point, k ₁ and k ₂ are the slopes of the linear fitting curves of the relationship between the atmospheric delay difference and the horizontal distance and elevation difference, respectively, abs(· ) means to take the absolute value.

融合高程信息的改进反距离加权内插法是根据地理空间相关性理论，在空间自相关假设基础上，进一步假设多种地理要素的空间分布具有内在关联性，由一种要素的空间分布推测其它要素空间分布的空间异相关模型。因此，融合高程信息的改进反距离加权内插法更有利于当待估站点高程较周围采样站点高程有较大变化时插值精度的提高。The improved inverse distance weighted interpolation method for fusion of elevation information is based on the geospatial correlation theory, on the basis of the spatial autocorrelation assumption, and further assumes that the spatial distribution of various geographical elements is inherently relevant, and the spatial distribution of one element is used to infer the other. A spatial heterocorrelation model for the spatial distribution of features. Therefore, the improved inverse distance weighted interpolation method fused with elevation information is more conducive to the improvement of interpolation accuracy when the elevation of the station to be estimated has a large change compared with the elevation of the surrounding sampling stations.

10、对大气改正后的InSAR差分干涉图进行解缠10. Unwrap the InSAR differential interferogram after atmospheric correction

本发明设计了两种GPS辅助InSAR解缠算法：GPS辅助InSAR Goldstein枝切线解缠算法和引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法。这两种方法有各自的优缺点，GPS辅助InSAR Goldstein枝切线解缠算法采用的是精确的积分算法，解缠精度高且计算速度快，但只能在局部解缠，且解缠范围和精度受噪声水平影响较大；引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法是一种全局解缠算法，抑噪能力强，但由于采用的是模拟退火寻优搜索算法，解缠速度较慢，且不可避免地出现陷入局部最优的可能，从而导致解缠的误差。因此，用户可根据不同的需要选择不同的解缠算法。两种解缠算法的基本原理具体如下：The present invention designs two kinds of GPS-assisted InSAR unwrapping algorithms: GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm and GPS-assisted InSAR unwrapping algorithm based on Markov random field which introduces residual points. These two methods have their own advantages and disadvantages. The GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm uses an accurate integral algorithm, which has high unwrapping accuracy and fast calculation speed, but it can only be unwrapped locally, and the unwrapping range and accuracy It is greatly affected by the noise level; the GPS-assisted InSAR unwrapping algorithm based on Markov random field that introduces residual points is a global unwrapping algorithm with strong noise suppression ability. The speed is slow, and it is inevitable to fall into the local optimum, which will lead to unwrapping errors. Therefore, users can choose different unwrapping algorithms according to different needs. The basic principles of the two unwrapping algorithms are as follows:

(1)GPS辅助InSAR Goldstein枝切线解缠算法(1) GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm

在利用Goldstein枝切线法进行相位解缠时，会出现这样一种情况：某一小区域被枝切线包围而与干涉图的其他区域完全隔离，这样一小区域被形象地称为“解缠孤岛”。“解缠孤岛”的形成原因很多，如水中小岛、植被包围的空地等。对于“解缠孤岛”，传统的处理方法有：(a)当孤岛范围较小时，可以在解缠时放弃对其进行解缠而仅对其他区域进行解缠，当求得其他区域的最终产品(如DEM或形变量)后，再利用适当的插值算法对孤岛区域进行插值，从而得到孤岛区域的最终产品。这种方法的缺点是忽略了孤岛上的相位信息，且插值过程会引入较大误差。(b)当孤岛范围较大时，在孤岛区域内重新选择一解缠起算点，对孤岛区域进行独立的相位解缠。这种方法虽然用到了孤岛上的相位信息，但由于采用了新的解缠起算点，使得孤岛区域与干涉图其他区域解缠基准不统一，造成最终产品解译上的困难。针对传统处理方法中存在的问题，由于InSAR形变测量精度通常是厘米级，而GPS形变测量的精度可达到毫米级，本发明假设离散的GPS形变观测量为已知值，采用GPS辅助Goldstein枝切线法进行InSAR相位解缠。When using the Goldstein branch tangent method for phase unwrapping, there will be such a situation: a small area is surrounded by branch tangent lines and completely isolated from other areas of the interferogram, such a small area is vividly called "unwrapping island". ". There are many reasons for the formation of "untangle islands", such as small islands in the water, open spaces surrounded by vegetation, etc. For the "unwrapping island", the traditional processing methods are: (a) When the island is small, you can give up unwrapping it during unwrapping and only unwrap other areas. When the final product of other areas is obtained (such as DEM or deformation), and then use an appropriate interpolation algorithm to interpolate the island area, so as to obtain the final product of the island area. The disadvantage of this method is that the phase information on the island is ignored, and the interpolation process will introduce a large error. (b) When the island range is large, re-select a starting calculation point for unwrapping in the island area, and perform independent phase unwrapping on the island area. Although this method uses the phase information on the island, due to the adoption of a new starting point for unwrapping, the unwrapping references of the island area and other areas of the interferogram are inconsistent, resulting in difficulties in the interpretation of the final product. Aiming at the problems existing in the traditional processing method, since the InSAR deformation measurement accuracy is usually at the centimeter level, while the GPS deformation measurement accuracy can reach the millimeter level, the present invention assumes that the discrete GPS deformation observations are known values, and adopts the GPS-assisted Goldstein branch tangent line method for InSAR phase unwrapping.

GPS辅助InSAR Goldstein枝切线解缠算法的原理为：The principle of the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm is as follows:

引入离散GPS点，选一GPS点为起算点，求其它GPS点相对于起算点的相对高程或形变，并利用下式将其转化为解缠相位值：Introduce discrete GPS points, select a GPS point as the starting point, calculate the relative elevation or deformation of other GPS points relative to the starting point, and use the following formula to convert it into an unwrapped phase value:

$\{\begin{matrix} {φ φ}_{topo topo}^{u u} = = - - \frac{44 π π}{λ λ} \frac{{B B}_{&perp; &perp;}^{00}}{ρ ρ sin sin {θ θ}_{00}} Δh Δh \\ {φ φ}_{defo defo}^{u u} = = \frac{44 π π}{λ λ} Δr Δr \end{matrix} - - - - - - ((1313))$

其中，φ_topo ^u和φ_defo ^u分别为解缠地形干涉相位和解缠形变干涉相位，Δh和Δr分别为相对于总解缠起算点的相对高程和形变。Among them, φ _topo ^u and φ _defo ^u are the unwrapped terrain interferometric phase and unwrapped deformation interferometric phase, respectively, and Δh and Δr are the relative elevation and deformation relative to the starting point of total unwrapping, respectively.

然后以孤岛点为起算点重新解缠，直到所有孤岛点解算完成。Then use the island point as the starting point to re-unwrap until all the island points are solved.

此方法是一种局部解缠算法，具有解缠速度快、精度高、基准统一且能评价解缠精度的优点。This method is a local unwrapping algorithm, which has the advantages of fast unwrapping speed, high precision, uniform benchmark and the ability to evaluate unwrapping accuracy.

(2)引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法(2) GPS-assisted InSAR unwrapping algorithm based on Markov random field with the introduction of residual points

残差点的识别与连接是Goldstein枝切线解缠算法的最大特点，它可以有效地隔离误差点，避免误差的全局传播。但Goldstein枝切线解缠可能产生孤岛，且枝切线上的像元点也不能得到解缠，降低了解缠的有效范围。Gudmundsson提出的基于马尔科夫随机场的GPS辅助解缠算法，没有考虑残差点的影响，会造成误差的全局传播。由此，本发明将残差点识别技术引入基于马尔科夫随机场的GPS辅助解缠算法，提出了引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法，实现解缠精度与解缠范围的最佳综合。此方法将InSAR解缠过程转化为求能量函数U(K＝k|Y＝I_w)最小值的过程，即：The identification and connection of residual points is the biggest feature of Goldstein branch tangent unwrapping algorithm, which can effectively isolate error points and avoid the global propagation of errors. However, the unwrapping of Goldstein branch tangents may generate isolated islands, and the pixel points on the branch tangents cannot be unwrapped, which reduces the effective range of unwrapping. The GPS-assisted unwrapping algorithm based on Markov random field proposed by Gudmundsson does not consider the influence of residual points, which will cause global propagation of errors. Therefore, the present invention introduces the residual point identification technology into the GPS-assisted unwrapping algorithm based on the Markov random field, and proposes the GPS-assisted InSAR unwrapping algorithm based on the Markov random field that introduces residual points, so as to realize unwrapping accuracy and unwrapping accuracy. The best synthesis of the range. This method transforms the InSAR unwrapping process into a process of finding the minimum value of the energy function U(K=k|Y=I _w ), namely:

$P P ((K K = = k k | | Y Y = = {I I}_{w w})) &Proportional; &Proportional; exp exp ((- - \frac{U u ((K K = = k k | | Y Y = = {I I}_{w w}))}{T T})) - - - - - - ((1414))$

采用反复快速退火策略以较高速率和概率获得全局最优解，即给定一较大温度衰减率cool和终止温度T_s，当T＜T_s时，令T＝T₀，反复快速退火，直到达到全局最优解。迭代时所用到的初始整周数矩阵k，初始解缠干涉图I_u，解缠干涉图I′_u和比例系数r计算公式如下：Using repeated rapid annealing strategy to obtain the global optimal solution with higher rate and probability, that is, given a large temperature decay rate cool and termination temperature T _s , when T<T _s , let T=T ₀ , repeated rapid annealing, until the global optimal solution is reached. The calculation formulas of the initial integer cycle number matrix k, the initial unwrapped interferogram I _u , the unwrapped interferogram I′ _u and the scaling factor r used in the iteration are as follows:

$\{\begin{matrix} k k = = round round ((\frac{{I I}_{u u}^{GPS GPS} - - {I I}_{w w}}{\frac{λ λ}{22}})) \\ {I I}_{u u} = = {I I}_{w w} + + k k \frac{λ λ}{22} \\ {I I}_{u u}^{' '} = = {I I}_{w w} + + {k k}^{' '} \frac{λ λ}{22} \\ r r = = \frac{P P ((K K = = {k k}^{' '} | | Y Y = = {I I}_{u u}^{' '}))}{P P ((K K = = k k | | Y Y = = {I I}_{u u}))} = = exp exp ((- - \frac{U u ((K K = = {k k}^{' '} | | Y Y = = {I I}_{u u}^{' '})) - - U u ((K K = = k k | | Y Y = = {I I}_{u u}))}{T T})) \end{matrix} - - - - - - ((1515))$

其中，k为初始整周数矩阵，round(□)为取整算子，I_u ^GPS为由GPS插值结果反演的解缠干涉图，I_w为缠绕干涉图，λ为波长，I_u为初始解缠干涉图，I′_u为解缠干涉图，r为比例系数。Among them, k is the initial integer cycle number matrix, round(□) is the rounding operator, I _u ^GPS is the unwrapped interferogram inverted from the GPS interpolation result, I _w is the twisted interferogram, λ is the wavelength, and I _u is The initial unwrapped interferogram, I′ _u is the unwrapped interferogram, and r is the scaling factor.

与Goldstein枝切线解缠算法先进行残差点识别再进行枝切线连接的方法不同，本发明采用的残差点识别和标记策略为：计算2×2大小的封闭路径单元累积值S，若S不为零则将所有4个像元点都标记成残差点，从而完成残差点的识别和标记。这样标记残差点的优势在于：(1)避免了传统残差点识别方法中人为指定残差点位置(左上角点)可能导致的残差点位置标记错误，从而更加有效地抑制误差的传播；(2)不再进行复杂的枝切线连接，从而提高计算效率。Different from the method of Goldstein branch tangent unwrapping algorithm, which first identifies residual points and then connects branch tangents, the residual point identification and marking strategy adopted in the present invention is: calculate the cumulative value S of closed path units with a size of 2×2, if S is not Zero marks all 4 pixel points as residual points, thus completing the identification and marking of residual points. The advantages of marking the residual points in this way are: (1) avoiding the error of marking the position of the residual points that may be caused by artificially specifying the position of the residual points (upper left corner point) in the traditional residual point identification method, thereby more effectively suppressing the propagation of errors; (2) Complicated branch tangent connections are no longer required, thereby improving computational efficiency.

11、将InSAR解缠相位转换为形变信息11. Convert InSAR unwrapped phase into deformation information

其中，□R_d为视线向的形变量，φ_d为解缠相位。Among them, □R _d is the deformation in the line-of-sight direction, and φ _d is the unwrapping phase.

12、InSAR形变信息地理编码12. Geocoding of InSAR deformation information

这一步是将前面处理得到的雷达坐标系下的结果转化到地理坐标系下的过程，采用距离多普勒(R-D)定位模型进行。在R-D定位模型中，SAR影像上任意像元的位置由三个方程确定：距离方程、多普勒质心平面方程和地球模型方程。This step is the process of transforming the results in the radar coordinate system obtained from the previous processing into the geographic coordinate system, using the range-Doppler (R-D) positioning model. In the R-D positioning model, the position of any pixel on the SAR image is determined by three equations: the distance equation, the Doppler centroid plane equation and the earth model equation.

距离方程： $R = {[(\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}}) \cdot (\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}})]}^{\frac{1}{2}}$ Distance equation: $R = {[(\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}}) \cdot (\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}})]}^{\frac{1}{2}}$

多普勒方程： $f_{D} = \frac{2}{λR} (\overset{&RightArrow;}{V_{s}} - \overset{&RightArrow;}{V_{t}}) \cdot (\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}})$ Doppler equation: $f_{D.} = \frac{2}{λR} (\overset{&Right Arrow;}{V_{the s}} - \overset{&Right Arrow;}{V_{t}}) \cdot (\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}})$

地球模型方程： $\frac{x_{t}^{2}}{{(a + h)}^{2}} + \frac{y_{t}^{2}}{{(a + h)}^{2}} + \frac{z_{t}^{2}}{b^{2}} = 1$ Earth model equations: $\frac{x_{t}^{2}}{{(a + h)}^{2}} + \frac{{the y}_{t}^{2}}{{(a + h)}^{2}} + \frac{z_{t}^{2}}{b^{2}} = 1$

式中，R为SAR传感器到地面点间的距离，

分别为卫星和地面点的位置矢量，“·”为矢量点乘，f_D为多普勒频移，λ为雷达波长，

为卫星速度矢量，

为地面点速度矢量，当采用WGS-84坐标系统时，

的数值为零，x_t，y_t和z_t为地面点位置矢量

在x，y和z轴上的投影，a为地球长半轴，a＝6378.137km，b为地球短半轴，b＝6356.752km，h为平均正常高。通过对上述三个方程的求解即可完成测量成果的地理编码过程。In the formula, R is the distance between the SAR sensor and the ground point,

are the position vectors of the satellite and the ground point respectively, "·" is the vector dot product, f _D is the Doppler frequency shift, λ is the radar wavelength,

is the satellite velocity vector,

is the velocity vector of the ground point, when using the WGS-84 coordinate system,

The value of is zero, x _t , y _t and z _t are ground point position vectors

Projection on x, y and z axes, a is the semi-major axis of the earth, a=6378.137km, b is the semi-minor axis of the earth, b=6356.752km, h is the average normal height. The geocoding process of measurement results can be completed by solving the above three equations.

13、采用基于改进马尔科夫随机场的GPS和InSAR融合方法，融合InSAR和GPS的形变信息获取高时空分辨率的地表三维形变监测结果13. Adopt the GPS and InSAR fusion method based on the improved Markov random field, and integrate the deformation information of InSAR and GPS to obtain the three-dimensional deformation monitoring results of the surface with high temporal and spatial resolution

首先，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息，获得高空间分辨率的三维形变场。采用分解优化法获取最优结果，其中对GPS约束条件的权值确定采用交叉验证定权的方法。具体原理如下：First, the one-dimensional deformation results of the InSAR line of sight are decomposed into three-dimensional deformation information through GPS constraints, and the three-dimensional deformation field with high spatial resolution is obtained. Decomposition optimization method is used to obtain the optimal result, and the weight of GPS constraints is determined by cross-validation and weight determination. The specific principles are as follows:

构建InSAR视线向的一维形变结果与三维形变结果之间的关系：Construct the relationship between the one-dimensional deformation results and the three-dimensional deformation results of the InSAR line of sight:

d_los＝sinθcos(α_h-3π/2)d_n+sinθsin(α_h-3π/2)d_e-cosθd_u d _los ＝sinθcos(α _h -3π/2)d _n +sinθsin(α _h -3π/2)d _e -cosθd _u

(17)...

＝[u_n，u_e，u_u][d_n，d_e，d_u]^T =[u _n , u _e , u _u ][d _n , d _e , d _u ] ^T

其中，d_los为InSAR视线向形变测量，[u_n，u_e，u_u]为视线向单位矢量，[u_n，u_e，u_u]＝[0.34-0.095，0.935]，d_n、d_e和d_u为北东天方向优化值，θ为雷达视角，其值随像元点至卫星星下点的距离而改变；α_h为卫星飞行方向与北方向夹角，其值可以从SAR数据中读出。Among them, d _los is the InSAR line-of-sight deformation measurement, [u _n , u _e , u _u ] is the line-of-sight unit vector, [u _n , u _e , u _u ]=[0.34-0.095, 0.935], d _n , d _e and d _u are the optimized values of the north-east sky direction, θ is the radar viewing angle, and its value changes with the distance from the pixel point to the sub-satellite point of the satellite; α _h is the angle between the flying direction of the satellite and the north direction, and its value can be obtained from the SAR data read out.

确定利用马尔科夫随机场进行建模计算的总的能量函数：Determine the overall energy function for modeling computations using Markov random fields:

$U u = = {Σ Σ}_{i i = = 11}^{M m} {{{W W}_{ins ins}^{i i} {(({d d}_{los los}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i}))}^{22} - - - - - - ((1818))$

$+ + {W W}_{krig krig}^{ni ni} {(({d d}_{krig krig}^{ni ni} - - {d d}_{n no}^{i i}))}^{22} + + {W W}_{krig krig}^{ei ei} {(({d d}_{krig krig}^{ei ei} - - {d d}_{e e}^{i i}))}^{22} + + {W W}_{krig krig}^{ui ui} {(({d d}_{krig krig}^{ui ui} - - {d d}_{u u}^{i i}))}^{22}}}$

其中，in,

${W W}_{ins ins}^{i i} = = \frac{11}{22 {(({σ σ}_{ins ins}^{i i}))}^{22}}$

$\{\begin{matrix} {W W}_{krig krig}^{ni ni} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ni ni}))}^{22}} \\ {W W}_{krig krig}^{ei ei} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ei ei}))}^{22}} \\ {W W}_{krig krig}^{ui ui} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ui ui}))}^{22}} \end{matrix}$

式(18)中，U为总能量，d_los ⁱ为InSAR视线向形变测量，[u_n ⁱ，u_e ⁱ，u_u ⁱ]为视线向单位矢量，d_n ⁱ、d_e ⁱ和d_u ⁱ为北东天方向优化值，σ_ins ⁱ为InSAR观测误差，d_krig ⁿⁱ、d_krig ^ei和d_krig ^ui为北东天方向克里格插值，σ_krig ⁿⁱ、σ_krig ^ei和σ_krig ^ui为北东天方向克里格插值标准差。In formula (18), U is the total energy, d _los ⁱ is the InSAR line-of-sight deformation measurement, [ _un ⁱ , u _e ⁱ , u _u ⁱ ] is the line-of-sight unit vector, d _ni ^, d _e ⁱ and d _u ⁱ is the optimized value in the NE direction, σ _ins ⁱ is the InSAR observation error, d _krig ⁿⁱ , d _krig ^ei and d _krig ^ui are Kriging interpolation values in the NE direction, σ _krig ⁿⁱ , σ _krig ^ei and σ _krig ^ui are Standard deviation of kriging interpolation in the northeast direction.

式(18)中包含了4×M个非负项，因此，当所有非负项达到最小时，式(18)便达到了全局最小。于是对d_n ⁱ、d_e ⁱ和d_u ⁱ求偏导并令其为零得，可得解析优化法的计算公式，即：Formula (18) contains 4×M non-negative items, so when all non-negative items reach the minimum, formula (18) reaches the global minimum. Then take the partial derivatives ^of d _ni , d _e ⁱ and d _u ⁱ and make them zero, and the calculation formula of the analytical optimization method can be obtained, namely:

$[\begin{matrix} \frac{&PartialD; &PartialD; U u}{&PartialD; &PartialD; {d d}_{n no}^{i i}} \\ \frac{&PartialD; &PartialD; U u}{&PartialD; &PartialD; {d d}_{e e}^{i i}} \\ \frac{&PartialD; &PartialD; U u}{&PartialD; &PartialD; {d d}_{u u}^{i i}} \end{matrix}] = = [\begin{matrix} - - 22 {W W}_{ins ins}^{i i} {u u}_{n no}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ni ni} (({d d}_{krig krig}^{ni ni} - - {d d}_{n no}^{i i})) \\ - - 22 {W W}_{ins ins}^{i i} {u u}_{e e}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ei ei} (({d d}_{krig krig}^{ei ei} - - {d d}_{e e}^{i i})) \\ - - 22 {W W}_{ins ins}^{i i} {u u}_{u u}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ui ui} (({d d}_{krig krig}^{ui ui} - - {d d}_{u u}^{i i})) \end{matrix}] = = 00 - - - - - - ((1919))$

采用式(18)中的权值确定方式存在如何评价InSAR测量精度和克里格插值精度两个问题，为有效解决这两个问题，对GPS约束条件的权值，本发明采用交叉验证定权的方法确定权值，具体过程如下：There are two problems of how to evaluate InSAR measurement accuracy and Kriging interpolation accuracy in the weight determination method in formula (18). In order to effectively solve these two problems, the present invention uses cross-validation to determine the weight of the GPS constraints The method to determine the weight, the specific process is as follows:

●对GPS与InSAR观测值进行配准，并将GPS三维观测值换算到InSAR视线方向，按式(20)估计GPS观测站对应InSAR像元测量方差：●Register the GPS and InSAR observations, convert the GPS three-dimensional observations to the InSAR line of sight direction, and estimate the measurement variance of the InSAR pixel corresponding to the GPS observation station according to formula (20):

${σ σ}_{ins ins,, GPS GPS}^{22} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {(({d d}_{GPS GPS,, i i}^{los los} - - {d d}_{ins ins,, i i}^{los los}))}^{22} - - - - - - ((2020))$

式中，n为GPS观测站个数，d_GPS，i ^los为第i个GPS观测站三维观测值换算到InSAR视线向结果，d_ins，i ^los为InSAR第i个像元视线向观测值。In the formula, n is the number of GPS observation stations, d _{GPS, i} ^los is the conversion result of the three-dimensional observation value of the i-th GPS observation station to the InSAR line of sight, d _{ins, i} ^los is the line-of-sight observation value of the i-th pixel of InSAR.

●计算InSAR每一像元点的测量方差：● Calculate the measurement variance of each pixel point of InSAR:

${σ σ}_{ins ins . . i i}^{22} = = \frac{E E. (({γ γ}_{GPS GPS}))}{{γ γ}_{i i}} {σ σ}_{ins ins,, GPS GPS}^{22} - - - - - - ((21 twenty one))$

式中，σ_ins，i ²和γ_i分别为InSAR第i个像元点测量方差估值和相干系数，E(γ_GPS)为GPS观测站对应InSAR像元相干系数均值。In the formula, σ _{ins, i} ² and γ _i are the measurement variance estimate and coherence coefficient of the i-th pixel of InSAR, respectively, and E(γ _GPS ) is the mean value of the coherence coefficient of the InSAR pixel corresponding to the GPS observation station.

●将GPS观测值分成A、B两组，对A组进行普通克里格插值，求B组GPS三维观测值与其对应插值的方差：●Divide the GPS observations into two groups, A and B, perform ordinary kriging interpolation on group A, and calculate the variance of the three-dimensional GPS observations of group B and their corresponding interpolation values:

$\{\begin{matrix} {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, n no}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{n no} - - {d d}_{krig krig,, i i}^{n no}))}^{22} \\ {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, e e}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{e e} - - {d d}_{krig krig,, i i}^{e e}))}^{22} \\ {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, u u}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{u u} - - {d d}_{krig krig,, i i}^{u u}))}^{22} \end{matrix} - - - - - - ((22 twenty two))$

式中，n′为B组GPS观测站个数，d_GPS，i ⁿ、d_GPS，i ^e和d_GPS，i ^u为北东天方向B组第i个GPS观测站观测值，d_krig，i ⁿ、d_krig，i ^e和d_krig，i ^u为北东天方向B组第i个GPS观测站克里格插值。In the formula, n′ is the number of GPS observation stations in Group B, d _{GPS, i} ⁿ , d _{GPS, i} ^e and d _{GPS, i} ^u is the observation value of the i-th GPS observation station in Group B in the northeast direction, d _{krig, i} ⁿ , d _{krig, i} ^e and d _{krig, i} ^u are the kriging interpolation values of the i-th GPS observation station in group B in the northeast sky direction.

●计算克里格插值方差改正因子：●Calculate the kriging interpolation variance correction factor:

$\{\begin{matrix} {k k}^{n no} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, n no}^{22}}{E E. [[{(({σ σ}_{krig krig,, n no}^{22}))}_{B B}]]} \\ {k k}^{e e} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, e e}^{22}}{E E. [[{(({σ σ}_{krig krig,, e e}^{22}))}_{B B}]]} \\ {k k}^{u u} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, u u}^{22}}{E E. [[{(({σ σ}_{krig krig,, u u}^{22}))}_{B B}]]} \end{matrix} - - - - - - ((23 twenty three))$

式中，E[(σ_krig，n ²)_B]、E[(σ_krig，e ²)_B]和E[(σ_krig，u ²)_B]分别为B组GPS观测站北东天方向克里格插值方差均值。In the formula, E[(σ _{krig, n} ² ) _B ], E[(σ _{krig, e} ² ) _B ], and E[(σ _{krig, u} ² ) _B ] are respectively the gram Rigging interpolation variance mean.

●对克里格方差进行改正：●Correction for kriging variance:

$\{\begin{matrix} {(({σ σ}_{cor cor,, i i}^{n no}))}^{22} = = {k k}^{n no} {(({σ σ}_{krig krig,, i i}^{n no}))}^{22} \\ {(({σ σ}_{cor cor . . i i}^{e e}))}^{22} = = {k k}^{e e} {(({σ σ}_{krig krig,, i i}^{e e}))}^{22} \\ {(({σ σ}_{cor cor,, i i}^{u u}))}^{22} = = {k k}^{u u} {(({σ σ}_{krig krig,, i i}^{u u}))}^{22} \end{matrix} - - - - - - ((24 twenty four))$

式中，(σ_krig，i ⁿ)²、(σ_krig，i ^e)²和(σ_krig，i ^u)²为第i个像元点北东天方向克里格插值方差。In the formula, (σ _{krig, i} ⁿ ) ² , (σ _{krig, i} ^e ) ² and (σ _{krig, i} ^u ) ² are the kriging interpolation variance in the north-east direction of the i-th pixel point.

●确定权值：● Determine the weight:

${W W}_{ins ins,, i i} = = \frac{11}{22 {σ σ}_{ins ins,, i i}^{22}} - - - - - - ((2525))$

$\{\begin{matrix} {W W}_{krig krig}^{ni ni} = = \frac{11}{22 {(({σ σ}_{cor cor,, i i}^{n no}))}^{22}} \\ {W W}_{krig krig}^{ei ei} = = \frac{11}{22 {(({σ σ}_{cor cor,, i i}^{e e}))}^{22}} \\ {W W}_{krig krig}^{ui ui} = = \frac{11}{22 {(({σ σ}_{cor cor,, i i}^{u u}))}^{22}} \end{matrix} - - - - - - ((2626))$

然后，对高时间频率GPS数据(其所采集的地面沉降量)在时间域上对上述高空间分辨率的三维形变场进行克里格插值和加密，从而将准GPS结果在时间域内加密成准GPS时间序列，获取高时间分辨率的三维形变场。Then, Kriging interpolation and encryption are performed on the above-mentioned three-dimensional deformation field with high spatial resolution for the high-time-frequency GPS data (the collected ground subsidence) in the time domain, so that the quasi-GPS results are encrypted into quasi-GPS results in the time domain. GPS time series to obtain 3D deformation field with high time resolution.

最后利用卡尔曼滤波对格网中的所有点进行估计获得高时空分辨率的地表三维形变结果。Finally, the Kalman filter is used to estimate all the points in the grid to obtain the three-dimensional surface deformation results with high spatial and temporal resolution.

本发明基于InSAR与GPS数据融合的地表三维形变监测方法，由于综合了InSAR与GPS这两种技术的优点，突破单一技术应用的局限，为地表监测提供了一种高时空分辨率、高精度的新技术，适合在大范围地表形变监测中使用。The present invention is based on the three-dimensional surface deformation monitoring method based on the fusion of InSAR and GPS data. Due to the combination of the advantages of InSAR and GPS, it breaks through the limitations of a single technology application and provides a high temporal and spatial resolution and high precision for surface monitoring. New technology, suitable for use in large-scale surface deformation monitoring.

附图说明 Description of drawings

图1：本发明方法的数据处理的流程图；Fig. 1: the flow chart of the data processing of the inventive method;

图2：GPS改正InSAR数据轨道误差的流程图；Figure 2: Flow chart of GPS correcting InSAR data orbit error;

图3：GPS反演大气延迟改正InSAR大气延迟误差的原理图；Figure 3: Schematic diagram of GPS retrieval atmospheric delay correction for InSAR atmospheric delay error;

图4：“解缠孤岛”产生示意图；Figure 4: Schematic diagram of the generation of "unwrapping islands";

图5：考虑残差点的基于马尔科夫随机场GPS辅助InSAR相位解缠数据处理流程图；Figure 5: Flowchart of data processing for GPS-assisted InSAR phase unwrapping based on Markov random field considering residual points;

具体实施方式 Detailed ways

下面结合实施例，对本发明作进一步的详细说明。Below in conjunction with embodiment, the present invention is described in further detail.

实施例：Example:

1、布设GPS点1. Set up GPS points

本方法布设两类GPS点：第一类是监测区内稳定的特征点，如河流、道路的交叉处，这类点主要用于将GPS数据和InSAR数据统一到同一参考坐标系下、改正SAR卫星轨道误差以及反演大气延迟；第二类是监测区的重点监测部位。This method lays out two types of GPS points: the first type is stable feature points in the monitoring area, such as the intersection of rivers and roads. This type of point is mainly used to unify GPS data and InSAR data into the same reference coordinate system and correct SAR Satellite orbit error and inversion atmospheric delay; the second category is the key monitoring parts of the monitoring area.

采集时间同步的GPS信号与SAR图像，采集的GPS信息相对于SAR图像的获取时间，向前和向后顺延10分钟左右。The time-synchronized GPS signal and SAR image are collected, and the collected GPS information is delayed by about 10 minutes forward and backward relative to the acquisition time of the SAR image.

两类GPS分开解算，解算的第一类GPS点的三维坐标信息包括：平面位置信息和高程信息，第二类GPS点的三维坐标信息包括：平面位移和垂直沉降。The two types of GPS are calculated separately. The three-dimensional coordinate information of the first type of GPS point includes: plane position information and elevation information, and the three-dimensional coordinate information of the second type of GPS point includes: plane displacement and vertical settlement.

4、SAR图像的斑点噪声抑制4. Speckle noise suppression for SAR images

本方法采用改进的Lee算法抑制SAR图形啊的斑点噪声。This method uses the improved Lee algorithm to suppress the speckle noise of SAR graphics.

首先，搜索3个以上InSAR图像上特征点目标地物，并用GPS测定其地面坐标。然后，采用下式计算两套坐标系统之间的转换矩阵估计值：Firstly, search for feature points on more than three InSAR images, and use GPS to measure their ground coordinates. Then, an estimate of the transformation matrix between the two coordinate systems is calculated using:

$\overset{^^}{X x} = = {(({B B}^{T T} PB PB))}^{- - 11} {B B}^{T T} PL PL$

其中，L为特征点目标的地面坐标矩阵，B为其对应影像坐标矩阵，D_□是观测误差

的方差矩阵。Among them, L is the ground coordinate matrix of the feature point target, B is the corresponding image coordinate matrix, D _□ is the observation error

The variance matrix of .

最后，利用转换矩阵将GPS求得的大气延迟改正、轨道改正等数据平移、旋转转换为InSAR影像坐标系下的数据。Finally, the transformation matrix is used to translate and rotate the atmospheric delay correction, orbit correction and other data obtained by GPS into the data in the InSAR image coordinate system.

在轨道误差的改正中，首先，将噪声抑制后的SAR图像作为输入数据，搜索6个以上InSAR图像上特征点目标地物，并用GPS测定其地面坐标。In the correction of track error, firstly, the noise-suppressed SAR image is used as input data, and more than 6 InSAR images are searched for feature points and target objects, and their ground coordinates are measured by GPS.

其次，构建地形测图方程，即Second, construct the topographic mapping equation, namely

$[\begin{matrix} X x \\ Y Y \\ Z Z \end{matrix}] = = λ λ [\begin{matrix} {a a}_{11} & {a a}_{22} & {a a}_{33} \\ {b b}_{11} & {b b}_{22} & {b b}_{33} \\ {c c}_{11} & {c c}_{22} & {c c}_{33} \end{matrix}] [\begin{matrix} 00 \\ r r sin sin θ θ \\ r r cos cos θ θ \end{matrix}] + + [\begin{matrix} {X x}_{s the s} \\ {Y Y}_{s the s} \\ {Z Z}_{s the s} \end{matrix}]$

$θ θ = = \frac{π π}{22} - - arccos arccos ((\frac{λc λc}{44 πB πB} φ φ)) + + α α$

其中，X，Y，Z为特征点目标地物的地面坐标，a₁，a₂，…，c₃为SAR卫星姿态角

构成的旋转矩阵元素，λ是雷达波长，c为光速，r为特征点目标地物到SAR图像的斜距，θ为视角，X_s，Y_s，Z_s为SAR卫星的轨道位置，φ为解缠相位。Among them, X, Y, Z are the ground coordinates of the feature point target object, a ₁ , a ₂ ,..., c ₃ are the attitude angles of the SAR satellite

λ is the radar wavelength, c is the speed of light, r is the slant distance from the feature point target object to the SAR image, θ is the viewing angle, X _s , Y _s , Z _s are the orbital position of the SAR satellite, φ is Unwrap phase.

X_s＝X_s0+X_s(t-t₀)X _s ＝X _s 0+X _s (tt ₀ )

B＝B₀+B(t-t₀)B＝B ₀ +B(tt ₀ )

V＝AX-LV=AX-L

式中，In the formula,

V＝[V_X，V_Y，V_Z]^T V=[V _X , V _Y , V _Z ] ^T

L＝[L_X，L_Y，L_Z]L=[L _X , L _Y , L _Z ]

其中，ΔX_s0，ΔY_s0，ΔZ_s0和分别为t₀时刻的轨道位置和姿态改正数；

和

分别为相应的变化率改正数；Δφ₀为相位常数改正数；a₁₁，a₁₂，…a₃₁₇为各未知参数的偏导数；V_x，V_Y，V_Z分别为X，Y，Z的已知值与其近似值的差。Among them, ΔX _s0 , ΔY _s0 , ΔZ _s0 and are the orbital position and attitude correction number at time t ₀ , respectively;

and

are the correction numbers of the corresponding rate of change; Δφ ₀ is the correction number of the phase constant; a ₁₁ , a ₁₂ ,...a ₃₁₇ are the partial derivatives of each unknown parameter; V _x , V _Y , V _Z are X, Y, Z The difference between a known value and its approximate value.

采用最小二乘法估计观测误差方程中的17个未知参数，由此得到SAR卫星轨道数据误差的改正值。The least square method is used to estimate 17 unknown parameters in the observation error equation, and thus the correction value of the SAR satellite orbit data error is obtained.

7、形成InSAR差分干涉图7. Form InSAR differential interferogram

首先采用随机过程法获取GPS的对流层估计，然后采用站间时域双差法估计大气延迟改正，并将其统一到InSAR的坐标系统中。其中，站间时域双差法为：First, the stochastic process method is used to obtain the tropospheric estimate of GPS, and then the inter-station time-domain double-difference method is used to estimate the atmospheric delay correction, which is unified into the InSAR coordinate system. Among them, the inter-station time domain double difference method is:

${D D.}_{AB AB}^{jk jk} = = {D D.}_{AB AB}^{k k} - - {D D.}_{AB AB}^{j j}$

其中，A点是SAR图像上的一个GPS参考点，B是同一张SAR图像上的另外一GPS点，j和k分别为主副图像获取时刻，D_A ^j和A_B ^j分别为将GPS解算得到的A、B两点在历元j的对流层延迟。Among them, point A is a GPS reference point on the SAR image, B is another GPS point on the same SAR image, j and k are the acquisition times of the main and secondary images respectively, D _A ^j and A _B ^j are the GPS solution The calculated tropospheric delay of points A and B at epoch j.

首先采用融合GPS高程信息的改进反距离加权内插法对GPS反演的大气延迟进行加密，然后将加密后的大气延迟从InSAR干涉相位图中消除。其中，融合GPS高程信息的改进反距离加权内插法为：First, the improved inverse distance weighted interpolation method combined with GPS elevation information is used to encrypt the atmospheric delay retrieved by GPS, and then the encrypted atmospheric delay is eliminated from the InSAR interferometric phase map. Among them, the improved inverse distance weighted interpolation method fused with GPS elevation information is:

$\{\begin{matrix} Z Z (({x x}_{00},, {y the y}_{00},, {z z}_{00})) = = {Σ Σ}_{i i = = 11}^{n no} {λ λ}_{i i} [[Z Z (({x x}_{i i},, {y the y}_{i i},, {z z}_{i i})) - - {m m}_{i i}]] + + {m m}_{00} \\ {m m}_{i i} = = {b b}_{00} + + {b b}_{i i} {h h}_{i i} \\ {λ λ}_{i i} = = \frac{{d d}_{i i 00}^{- - p p}}{{Σ Σ}_{i i = = 11}^{N N} {d d}_{i i 00}^{- - p p}} \\ {d d}_{i i 00} = = \frac{\sqrt{{((α α {d d}_{i i 00}^{s the s}))}^{22} + + {((((11 - - α α)) {d d}_{i i 00}^{h h}))}^{22}}}{α α} \\ {d d}_{i i 00}^{s the s} = = \sqrt{{(({x x}_{i i} - - {x x}_{00}))}^{22} + + {(({y the y}_{i i} - - {y the y}_{00}))}^{22}} \\ {d d}_{i i 00}^{h h} = = {h h}_{i i} - - {h h}_{00} ((i i = = 11,, . . . . . .,, n no)) \\ α α = = \frac{abs abs ((k k 11))}{abs abs ((k k 11)) + + abs abs ((k k 22))} \end{matrix}$

其中，Z(x₀，y₀，z₀)是点(x₀，y₀，z₀)处的估计值，Z(x_i，y_i，z_i)是点(x_i，y_i，z_i)处的观测值，n是用于估计的观测值的个数，λ_i是点(x_i，y_i，z_i)处观测值对应的权，h_i为点(x_i，y_i，z_i)处的高程，b₀、b₁为待求拟合系数，d_i0 ^s为观测点i到估值点的水平距离，d_i0 ^h为观测点i与估值点的高程差，k₁、k₂分别为大气延迟的差值随水平距离和高程差变化关系的线性拟合曲线的斜率，abs(·)代表取绝对值。Among them, Z(x ₀ , y ₀ , z ₀ ) is the estimated value at point (x ₀ , y ₀ , z ₀ ), Z( _xi , y _i , z _i ) is the point (xi , y _i _, z _i ), n is the number of observations used for estimation, λ _i is the weight corresponding to the observation at point (xi _, y _i , z _i ), h _i is the point ( _xi , y _i , z _i ), b ₀ , b ₁ are fitting coefficients to be obtained, d _i0 ^s is the horizontal distance from observation point i to evaluation point, d _i0 ^h is the elevation difference between observation point i and evaluation point , k ₁ and k ₂ are the slopes of the linear fitting curves of the relationship between the difference of atmospheric delay and the variation of horizontal distance and elevation difference, and abs(·) represents the absolute value.

本方法设计了两种GPS辅助InSAR解缠算法：GPS辅助InSAR Goldstein枝切线解缠算法和引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法。In this method, two GPS-assisted InSAR unwrapping algorithms are designed: GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm and GPS-assisted InSAR unwrapping algorithm based on Markov random field with the introduction of residual points.

这两种方法有各自的优缺点，GPS辅助InSAR Goldstein枝切线解缠算法一种局部解缠算法，具有解缠速度快、精度高、基准统一且能评价解缠精度的优点，但解缠范围和精度受噪声水平影响较大；引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法是一种全局解缠算法，抑噪能力强，但解缠速度较慢。用户可根据不同的需要选择不同的解缠算法。These two methods have their own advantages and disadvantages. The GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm is a local unwrapping algorithm, which has the advantages of fast unwrapping speed, high precision, unified benchmark and the ability to evaluate unwrapping accuracy. However, the unwrapping range The sum accuracy is greatly affected by the noise level; the GPS-assisted InSAR unwrapping algorithm based on Markov random field with the introduction of residual points is a global unwrapping algorithm with strong noise suppression ability but slow unwrapping speed. Users can choose different unwrapping algorithms according to different needs.

1)GPS辅助InSAR Goldstein枝切线解缠算法1) GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm

在GPS辅助InSAR Goldstein枝切线解缠算法中，数据处理的步骤如下：In the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm, the steps of data processing are as follows:

(1)将GPS观测站配准到缠绕干涉图中。(1) Register the GPS observation stations to the winding interferogram.

(2)求缠绕相位梯度场主值、识别残差点并进行枝切线的连接。(2) Calculate the main value of the winding phase gradient field, identify the residual points and connect the branch tangent lines.

(3)在位于孤岛外的GPS观测站中选定一总解缠起算点，如点1，求其它GPS点相对于点1的相对高程或形变，并利用下式将其转化为解缠相位值：(3) Select a total unwrapping starting point in the GPS observation station located outside the island, such as point 1, find the relative elevation or deformation of other GPS points relative to point 1, and use the following formula to convert it into the unwrapping phase value:

${φ φ}_{topo topo}^{u u} = = - - \frac{44 π π}{λ λ} \frac{{B B}_{&perp; &perp;}^{00}}{ρ ρ sin sin {θ θ}_{00}} Δh Δh$

${φ φ}_{defo defo}^{u u} = = \frac{44 π π}{λ λ} Δr Δr$

上式中，φ_topo ^u和φ_defo ^u分别为解缠地形干涉相位和解缠形变干涉相位，Δh和Δr分别为相对于总解缠起算点的相对高程和形变。In the above formula, φ _topo ^u and φ _defo ^u are the unwrapped terrain interferometric phase and unwrapped deformation interferometric phase, respectively, and Δh and Δr are the relative elevation and deformation relative to the starting point of total unwrapping, respectively.

(4)以点1为解缠起算点用Goldstein枝切线法对干涉图进行解缠，直到所有可解缠点被解缠。以图4所示情况为例，在这个过程中，孤岛外的所有点将被成功解缠，而孤岛内的点不能被解缠。(4) Use the Goldstein branch tangent method to unwrap the interferogram with point 1 as the starting point for unwrapping until all unwrapping points are unwrapped. Taking the situation shown in Figure 4 as an example, in this process, all points outside the island will be successfully unwrapped, while points inside the island cannot be unwrapped.

(5)对干涉图进行扫描，找到尚未被解缠的GPS点，即点2，并以该点为新的解缠起算点进行新一轮相位解缠。(5) Scan the interferogram to find the GPS point that has not been unwrapped, that is, point 2, and use this point as the starting point for new unwrapping to perform a new round of phase unwrapping.

(6)重复步骤5，直到整个干涉图解缠完毕。(6) Repeat step 5 until the entire interferogram is unwrapped.

(7)找到已解缠但未做过起算点的GPS点，即点3和点4，对其由步骤3反演得到的解缠相位已知值和由步骤4-5得到的解缠相位估计进行统计分析，从而实现对解缠算法的精度评价。(7) Find the GPS points that have been unwrapped but have not been calculated as the starting point, that is, point 3 and point 4, for which the known value of the unwrapped phase obtained by step 3 inversion and the unwrapped phase obtained by steps 4-5 Statistical analysis is performed on the estimation, so as to realize the accuracy evaluation of the unwrapping algorithm.

(8)将GPS反演的解缠相位已知值付给解缠干涉图上相应的点，整个解缠过程结束。(8) Pay the known value of the unwrapping phase retrieved by GPS to the corresponding point on the unwrapping interferogram, and the whole unwrapping process ends.

2)引入残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法2) GPS-assisted InSAR unwrapping algorithm based on Markov random field with the introduction of residual points

此方法将InSAR解缠过程转化为求能量函数U(K＝k|Y＝I_w)最小值的过程，即：This method transforms the InSAR unwrapping process into a process of finding the minimum value of the energy function U(K=k|Y=I _w ), namely:

$P P ((K K = = k k | | Y Y = = {I I}_{w w})) &Proportional; &Proportional; exp exp ((- - \frac{U u ((K K = = k k | | Y Y = = {I I}_{w w}))}{T T})) - - - - - - ((2727))$

$\{\begin{matrix} k k = = round round ((\frac{{I I}_{u u}^{GPS GPS} - - {I I}_{w w}}{\frac{λ λ}{22}})) \\ {I I}_{u u} = = {I I}_{w w} + + k k \frac{λ λ}{22} \\ {I I}_{u u}^{' '} = = {I I}_{w w} + + {k k}^{' '} \frac{λ λ}{22} \\ r r = = \frac{P P ((K K = = {k k}^{' '} | | Y Y = = {I I}_{u u}^{' '}))}{P P ((K K = = k k | | Y Y = = {I I}_{u u}))} = = exp exp ((- - \frac{U u ((K K = = {k k}^{' '} | | Y Y = = {I I}_{u u}^{' '})) - - U u ((K K = = k k | | Y Y = = {I I}_{u u}))}{T T})) \end{matrix}$

其中，残差点识别和标记策略为：计算2×2大小的封闭路径单元累积值S，若S不为零则将所有4个像元点都标记成残差点，从而完成残差点的识别和标记。残差点标记完成后利用基于马尔科夫随机场GPS辅助InSAR解缠算法对非残差点进行解缠。该步解缠结束后，将所有已解缠点标记为固定域再进行一次只采用平滑约束的基于马尔科夫随机场GPS辅助InSAR解缠算法，从而完成对整个干涉图的解缠。可见考虑残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法首先对质量好的像元点进行解缠，再以质量好的像元点的解缠值为约束进一步对质量不好的点进行解缠，从而最大限度的抑制误差传播。考虑残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法实现了解缠精度与范围的较好综合。考虑残差点的基于马尔科夫随机场GPS辅助InSAR解缠算法的数据处理流程如图5所示。Among them, the residual point identification and marking strategy is: calculate the cumulative value S of the closed path unit with a size of 2×2, and if S is not zero, mark all four pixel points as residual points, thereby completing the identification and marking of residual points . After the residual points are marked, the non-residual points are unwrapped using the GPS-assisted InSAR unwrapping algorithm based on Markov random field. After this step of unwrapping, all the unwrapped points are marked as fixed domains, and then a GPS-assisted InSAR unwrapping algorithm based on Markov random field using only smoothness constraints is performed to complete the unwrapping of the entire interferogram. It can be seen that the Markov random field-based GPS-assisted InSAR unwrapping algorithm considering the residual points first unwraps the pixel points with good quality, and then uses the unwrapping value of the good-quality pixel points as constraints to further unwrap the poor-quality points. Unwrapping is performed to minimize error propagation. Considering the residual points, the GPS-assisted InSAR unwrapping algorithm based on Markov random field achieves a better synthesis of unwrapping accuracy and range. The data processing flow of the GPS-assisted InSAR unwrapping algorithm based on Markov random field considering residual points is shown in Fig. 5.

本方法采用距离多普勒(R-D)定位模型，将前面处理得到的雷达坐标系下的结果转化到地理坐标系下的过程。在R-D定位模型中，SAR影像上任意像元的位置由三个方程确定：距离方程、多普勒质心平面方程和地球模型方程。This method adopts the range-Doppler (R-D) positioning model, and transforms the results in the radar coordinate system obtained from the previous processing into the process in the geographic coordinate system. In the R-D positioning model, the position of any pixel on the SAR image is determined by three equations: the distance equation, the Doppler centroid plane equation and the earth model equation.

距离方程： $R = {[(\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}}) \cdot (\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}})]}^{\frac{1}{2}}$ Distance equation: $R = {[(\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}}) &Center Dot; (\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}})]}^{\frac{1}{2}}$

多普勒方程： $f_{D} = \frac{2}{λR} (\overset{&RightArrow;}{V_{s}} - \overset{&RightArrow;}{V_{t}}) \cdot (\overset{&RightArrow;}{R_{s}} - \overset{&RightArrow;}{R_{t}})$ Doppler equation: $f_{D.} = \frac{2}{λR} (\overset{&Right Arrow;}{V_{the s}} - \overset{&Right Arrow;}{V_{t}}) &Center Dot; (\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}})$

式中，R为SAR传感器到地面点间的距离，

为卫星速度矢量，

为地面点速度矢量，当采用WGS-84坐标系统时，

的数值为零，x_t，y_t和z_t为地面点位置矢量

在x，y和z轴上的投影，a为长半轴，b为短半轴，h为平均正常高，a＝6378.137km，b＝6356.752km。通过对上述三个方程的求解即可完成测量成果的地理编码过程。In the formula, R is the distance between the SAR sensor and the ground point,

is the satellite velocity vector,

The value of is zero, x _t , y _t and z _t are ground point position vectors

Projection on x, y and z axes, a is the semi-major axis, b is the semi-minor axis, h is the average normal height, a=6378.137km, b=6356.752km. The geocoding process of measurement results can be completed by solving the above three equations.

首先，将InSAR视线向的一维形变结果通过GPS约束分解成三维形变信息。具体算法如下：First, the one-dimensional deformation results of InSAR line-of-sight are decomposed into three-dimensional deformation information through GPS constraints. The specific algorithm is as follows:

(1)建立视线向的一维形变结果与三维形变结果之间的关系(1) Establish the relationship between the one-dimensional deformation results of the line of sight and the three-dimensional deformation results

(28)...

其中，d_los为InSAR视线向形变测量，[u_n，u_e，u_u]为视线向单位矢量，[u_n，u_e，u_u]＝[0.34，-0.095，0.935]，d_n、d_e和d_u为北东天方向优化值，θ为雷达视角，其值随像元点至卫星星下点的距离而改变；α_h为卫星飞行方向与北方向夹角，其值可以从SAR数据中读出。Among them, d _los is the InSAR line-of-sight deformation measurement, [u _n , u _e , u _u ] is the line-of-sight unit vector, [u _n , u _e , u _u ]=[0.34, -0.095, 0.935], d _n , d _e and d _u are the optimized values of the north-east sky direction, θ is the radar viewing angle, and its value changes with the distance from the pixel point to the sub-satellite point of the satellite; α _h is the angle between the flying direction of the satellite and the north direction, and its value can be changed from read out from the SAR data.

(2)确定总能量函数(2) Determine the total energy function

$U u = = {Σ Σ}_{i i = = 11}^{M m} {{{W W}_{ins ins}^{i i} {(({d d}_{los los}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i}))}^{22}$

其中，in,

${W W}_{ins ins}^{i i} = = \frac{11}{22 {(({σ σ}_{ins ins}^{i i}))}^{22}}$

其中，W_ins ⁱ为InSAR约束权值，W_krig ⁿⁱ，W_krig ^ei，W_krig ^ui为GPS约束权值，d_los ⁱ为InSAR视线向形变测量，[u_n ⁱ，u_e ⁱ，u_u ⁱ]为视线向单位矢量，d_n ⁱ、d_e ⁱ和d_u ⁱ为北东天方向优化值，σ_ins ⁱ为InSAR观测误差，d_krig ⁿⁱ、d_krig ^ei和d_krig ^ui为北东天方向克里格插值，σ_krig ⁿⁱ、σ_krig ^ei和σ_krig ^ui为北东天方向克里格插值标准差。Among them, _Wins ⁱ is InSAR constraint weight, W _krig ⁿⁱ , W _krig ^ei , W _krig ^ui are GPS constraint weight, d _los ⁱ is InSAR line-of-sight deformation measurement, [u _n ⁱ , u _e ⁱ , u _u ⁱ ] is the line-of-sight unit vector, d _n ⁱ , d _e ⁱ and d _u ⁱ are the optimized values of the NE sky direction, σ _ins ⁱ is the InSAR observation error, d _krig ⁿⁱ , d _krig ^ei and d _krig ^ui are the NE sky direction Kriging interpolation, σ _krig ⁿⁱ , σ _krig ^ei and σ _krig ^ui are standard deviations of Kriging interpolation in the north-east direction.

(3)采用交叉验证定权法确定约束权值，具体如下：(3) Use the cross-validation weighting method to determine the constraint weight, as follows:

●对GPS与InSAR观测值进行配准，并将GPS三维观测值换算到InSAR视线方向。●Register the GPS and InSAR observations, and convert the GPS three-dimensional observations to the InSAR line-of-sight direction.

●按下式估计GPS观测站对应InSAR像元测量方差：●Estimate the measurement variance of the InSAR pixel corresponding to the GPS observation station according to the following formula:

${σ σ}_{ins ins,, GPS GPS}^{22} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {(({d d}_{GPS GPS,, i i}^{los los} - - {d d}_{ins ins,, i i}^{los los}))}^{22}$

●InSAR每一像元点的测量方差可由下式计算：The measurement variance of each pixel point in InSAR can be calculated by the following formula:

${σ σ}_{ins ins . . i i}^{22} = = \frac{E E. (({γ γ}_{GPS GPS}))}{{γ γ}_{i i}} {σ σ}_{ins ins,, GPS GPS}^{22}$

●将GPS观测值分成A、B两组，对A组进行普通克里格插值。● Divide the GPS observations into two groups, A and B, and perform ordinary kriging interpolation on group A.

●求B组GPS三维观测值与其对应插值的方差：●Find the variance of group B GPS three-dimensional observations and their corresponding interpolation values:

$\{\begin{matrix} {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, n no}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{n no} - - {d d}_{krig krig,, i i}^{n no}))}^{22} \\ {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, e e}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{e e} - - {d d}_{krig krig,, i i}^{e e}))}^{22} \\ {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, u u}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{u u} - - {d d}_{krig krig,, i i}^{u u}))}^{22} \end{matrix}$

$\{\begin{matrix} {k k}^{n no} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, n no}^{22}}{E E. [[{(({σ σ}_{krig krig,, n no}^{22}))}_{B B}]]} \\ {k k}^{e e} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, e e}^{22}}{E E. [[{(({σ σ}_{krig krig,, e e}^{22}))}_{B B}]]} \\ {k k}^{u u} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, u u}^{22}}{E E. [[{(({σ σ}_{krig krig,, u u}^{22}))}_{B B}]]} \end{matrix} - - - - - - ((2929))$

●对克里格方差进行改正：●Correction for kriging variance:

$\{\begin{matrix} {(({σ σ}_{cor cor,, i i}^{n no}))}^{22} = = {k k}^{n no} {(({σ σ}_{krig krig,, i i}^{n no}))}^{22} \\ {(({σ σ}_{cor cor . . i i}^{e e}))}^{22} = = {k k}^{e e} {(({σ σ}_{krig krig,, i i}^{e e}))}^{22} \\ {(({σ σ}_{cor cor,, i i}^{u u}))}^{22} = = {k k}^{u u} {(({σ σ}_{krig krig,, i i}^{u u}))}^{22} \end{matrix} - - - - - - ((3030))$

式中，(σ_krig，i ⁿ)₂、(σ_krig，i ^e)²和(σ_krig，i ^u)₂为第i个像元点北东天方向克里格插值方差。In the formula, (σ _{krig, i} ⁿ ) ₂ , (σ _{krig, i} ^e ) ² and (σ _{krig, i} ^u ) ₂ are the kriging interpolation variance in the north-east direction of the i-th pixel point.

●定权：●Fixed power:

${W W}_{ins ins}^{i i} = = \frac{11}{22 {(({σ σ}_{ins ins}^{i i}))}^{22}}$

(4)采用直接解析法获取三维形变场的最大后验估计(4) Using the direct analytical method to obtain the maximum a posteriori estimation of the three-dimensional deformation field

$[\begin{matrix} \frac{&PartialD; &PartialD; U u}{&PartialD; &PartialD; {d d}_{n no}^{i i}} \\ \frac{&PartialD; &PartialD; U u}{&PartialD; &PartialD; {d d}_{e e}^{i i}} \\ \frac{&PartialD; &PartialD; U u}{&PartialD; &PartialD; {d d}_{u u}^{i i}} \end{matrix}] = = [\begin{matrix} - - 22 {W W}_{ins ins}^{i i} {u u}_{n no}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ni ni} (({d d}_{krig krig}^{ni ni} - - {d d}_{n no}^{i i})) \\ - - 22 {W W}_{ins ins}^{i i} {u u}_{e e}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ei ei} (({d d}_{krig krig}^{ei ei} - - {d d}_{e e}^{i i})) \\ - - 22 {W W}_{ins ins}^{i i} {u u}_{u u}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ui ui} (({d d}_{krig krig}^{ui ui} - - {d d}_{u u}^{i i})) \end{matrix}] = = 00$

Claims

1. A method for monitoring surface three-dimensional deformation based on InSAR and GPS data fusion, the steps of which are:

Step 1. Set up GPS points

There are two types of GPS points: the first type is the stable feature point in the monitoring area, which is used to unify the GPS data and InSAR data into the same reference coordinate system, correct the SAR satellite orbit error and invert the atmospheric delay; the second type is the monitoring point of the monitoring area;

Step 2. Collect and record GPS signals and collect SAR images

The collected GPS signal should be synchronized with the time of SAR image acquisition;

Step 3: Calculate the three-dimensional coordinate information of two types of GPS points according to the GPS signal

The three-dimensional coordinate information of the first type of GPS monitoring point includes: position information and elevation information, and the three-dimensional coordinate information of the second type of GPS monitoring point includes: position information and deformation information;

Step 4. Speckle noise suppression of SAR images

Step five, search for feature points on the SAR image, and establish the coordinate transformation relationship between the searched feature points and the first type of GPS points

Step 6. Use GPS data to correct the orbit error of SAR satellites

A method combining surveying equation method and least square estimation is adopted;

Step 7. Form InSAR differential interferogram

Step 8. Use the first type of GPS points to invert the atmospheric delay and unify it into the InSAR coordinate system

The stochastic process method is used to obtain the GPS tropospheric estimate, and the inter-station time-domain double-difference method is used to estimate the atmospheric delay correction;

Step 9. Correct the atmospheric delay error of the InSAR differential interferogram by using the atmospheric delay retrieved by GPS

The improved inverse distance weighted interpolation method fused with GPS elevation information is used to encrypt the atmospheric delay retrieved by GPS, and then the encrypted atmospheric delay is eliminated from the InSAR differential interferogram;

Step 10. Unwrap the atmospherically corrected InSAR differential interferogram

Using the GPS-assisted InSAR Goldstein branch tangent unwrapping algorithm or the Markov random field-based GPS-assisted InSAR unwrapping algorithm that introduces residual points, the InSAR differential interferogram after atmospheric correction is unwrapped;

Step 11. Convert InSAR unwrapped phase into deformation information

Among them, □R _d is the deformation in the line-of-sight direction, φ _d is the unwrapping phase, and λ is the radar wavelength;

Step 12. Geocoding of InSAR deformation information

Distance equation:

R = {[(\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}}) \cdot (\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}})]}^{\frac{1}{2}}

Doppler equation:

f_{D.} = \frac{2}{λR} (\overset{&Right Arrow;}{V_{the s}} - \overset{&Right Arrow;}{V_{t}}) &Center Dot; (\overset{&Right Arrow;}{R_{the s}} - \overset{&Right Arrow;}{R_{t}})

Earth model equations:

\frac{x_{t}^{2}}{{(a + h)}^{2}} + \frac{{the y}_{t}^{2}}{{(a + h)}^{2}} + \frac{z_{t}^{2}}{b^{2}} = 1

In the formula, R is the distance between the SAR sensor and the ground point,

are the position vectors of the satellite and the ground point respectively, "·" is the vector dot product, f _D is the Doppler frequency shift,

is the satellite velocity vector,

is the velocity vector of the ground point, x _t , y _t and z _t are the position vectors of the ground point

Projection on the x, y and z axes, a is the semi-major axis of the earth, b is the semi-minor axis of the earth, h is the mean normal height;

By solving the above three equations, the geocoding of InSAR deformation information can be obtained;

Step 13, using the GPS and InSAR fusion method based on the improved Markov random field, fusing the deformation information of InSAR and GPS to obtain high-precision three-dimensional deformation monitoring results of the surface;

Finally, the Kalman filter is used to estimate all the points in the grid to obtain the three-dimensional surface deformation results with high spatial and temporal resolution.

2. the surface three-dimensional deformation monitoring method based on InSAR and GPS data fusion according to claim 1, is characterized in that: the GPS and the InSAR fusion method based on improved Markov random field described in step 13:

First, the one-dimensional deformation results of the InSAR line of sight are decomposed into three-dimensional deformation information through GPS constraints; the specific algorithm is as follows:

(1) Establish the relationship between the one-dimensional deformation results of the line of sight and the three-dimensional deformation results

d _los ＝sinθcos(α _h -3π/2)d _n +sinθsin(α _h -3π/2)d _e -cosθd _u

=[u _n , u _e , u _u ][d _n , d _e , d _u ] ^T

Among them, d _los is the InSAR line-of-sight deformation measurement, [u _n , u _e , u _u ] is the line-of-sight unit vector, [u _n , u _e , u _u ]=[0.34, -0.095, 0.935], d _n , d _e and d _u are the optimized values of the north-east sky direction, θ is the radar viewing angle, and its value changes with the distance from the pixel point to the sub-satellite point of the satellite; α _h is the angle between the satellite flight direction and the north direction, and its value is changed from read from the data;

(2) Determine the total energy function

U u = = {Σ Σ}_{i i = = 11}^{M m} {{{W W}_{ins ins}^{i i} {(({d d}_{los los}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} - - {d d}_{u u}^{i i}))}^{22}

+ + {W W}_{krig krig}^{ni ni} {(({d d}_{krig krig}^{ni ni} - - {d d}_{n no}^{i i}))}^{22} + + {W W}_{krig krig}^{ei ei} {(({d d}_{krig krig}^{ei ei} - - {d d}_{e e}^{i i}))}^{22} + + {W W}_{krig krig}^{ui ui} {(({d d}_{krig krig}^{ui ui} - - {d d}_{u u}^{i i}))}^{22}}}

in,

{W W}_{ins ins}^{i i} = = \frac{11}{22 {(({σ σ}_{ins ins}^{i i}))}^{22}}

\{\begin{matrix} {W W}_{krig krig}^{ni ni} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ni ni}))}^{22}} \\ {W W}_{krig krig}^{ei ei} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ei ei}))}^{22}} \\ {W W}_{krig krig}^{ui ui} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ui ui}))}^{22}} \end{matrix}

in,

is the InSAR constraint weight,

is the GPS constraint weight,

For InSAR line-of-sight deformation measurement, is the line-of-sight unit vector,

and

is the optimized value for the north-east direction,

is the InSAR observation error,

and

Kriging interpolation for the NE direction,

and

Kriging interpolation standard deviation for the northeast direction;

(3) Use the cross-validation weighting method to determine the constraint weight, as follows:

Register the GPS and InSAR observations, and convert the GPS three-dimensional observations to the InSAR line of sight direction;

The measurement variance of the GPS observation station corresponding to the InSAR pixel is estimated by the following formula:

{σ σ}_{ins ins,, GPS GPS}^{22} = = \frac{11}{n no} {Σ Σ}_{i i = = 11}^{n no} {(({d d}_{GPS GPS,, i i}^{los los} - - {d d}_{ins ins,, i i}^{los los}))}^{22}

In the formula, n is the number of GPS observation stations,

is the result of converting the three-dimensional observation value of the i-th GPS station to the InSAR line of sight,

is the line-of-sight observation value of the i-th pixel of InSAR;

The measurement variance of each pixel point in InSAR is calculated by the following formula:

{σ σ}_{ins ins,, i i}^{22} = = \frac{E E. (({γ γ}_{GPS GPS}))}{{γ γ}_{i i}} {σ σ}_{ins ins,, GPS GPS}^{22}

In the formula,

and γ _i are the measurement variance estimate and coherence coefficient of the i-th pixel of InSAR respectively, and E(γ _GPS ) is the mean value of the coherence coefficient of the InSAR pixel corresponding to the GPS observation station;

Divide the GPS observations into two groups, A and B, and perform ordinary kriging interpolation on group A;

Calculate the variance of the GPS three-dimensional observations of group B and their corresponding interpolation values:

\{\begin{matrix} {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, n no}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{n no} - - {d d}_{krig krig,, i i}^{n no}))}^{22} \\ {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, e e}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{e e} - - {d d}_{krig krig,, i i}^{e e}))}^{22} \\ {\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, u u}^{22} = = \frac{11}{{n no}^{' '}} {Σ Σ}_{i i = = 11}^{{n no}^{' '}} {(({d d}_{GPS GPS,, i i}^{u u} - - {d d}_{krig krig,, i i}^{u u}))}^{22} \end{matrix}

In the formula, n' is the number of GPS observation stations in group B,

and

is the observation value of the i-th GPS observation station in group B in the northeast sky direction, and Kriging interpolation for the i-th GPS observation station in group B in the northeast direction;

Compute the kriging variance correction factor:

\{\begin{matrix} {k k}^{n no} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, n no}^{22}}{E E. [[{(({σ σ}_{krig krig,, n no}^{22}))}_{B B}]]} \\ {k k}^{e e} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, e e}^{22}}{E E. [[{(({σ σ}_{krig krig,, e e}^{22}))}_{B B}]]} \\ {k k}^{u u} = = \frac{{\overset{&OverBar; &OverBar;}{σ σ}}_{krig krig,, u u}^{22}}{E E. [[{(({σ σ}_{krig krig,, u u}^{22}))}_{B B}]]} \end{matrix}

In the formula,

and

Respectively, the average variance of Kriging interpolation in the north-east direction of GPS observation stations in group B;

Correct for kriging variance:

\{\begin{matrix} {(({σ σ}_{cor cor,, i i}^{n no}))}^{22} = = {k k}^{n no} {(({σ σ}_{krig krig,, i i}^{n no}))}^{22} \\ {(({σ σ}_{cor cor,, i i}^{e e}))}^{22} = = {k k}^{e e} {(({σ σ}_{krig krig,, i i}^{e e}))}^{22} \\ {(({σ σ}_{cor cor,, i i}^{u u}))}^{22} = = {k k}^{u u} {(({σ σ}_{krig krig,, i i}^{u u}))}^{22} \end{matrix}

In the formula,

and Variance of kriging interpolation for the i-th pixel point in the north-east sky direction; fixed weight:

{W W}_{ins ins}^{i i} = = \frac{11}{22 {(({σ σ}_{ins ins}^{i i}))}^{22}}

\{\begin{matrix} {W W}_{krig krig}^{ni ni} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ni ni}))}^{22}} \\ {W W}_{krig krig}^{ei ei} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ei ei}))}^{22}} \\ {W W}_{krig krig}^{ui ui} = = \frac{11}{22 {(({σ σ}_{krig krig}^{ui ui}))}^{22}} \end{matrix}

(4) Using the direct analytical method to obtain the maximum a posteriori estimation of the three-dimensional deformation field

[\begin{matrix} \frac{&PartialD; &PartialD; U u}{{&PartialD; &PartialD; d d}_{n no}^{i i}} \\ \frac{&PartialD; &PartialD; U u}{{&PartialD; &PartialD; d d}_{e e}^{i i}} \\ \frac{&PartialD; &PartialD; U u}{{&PartialD; &PartialD; d d}_{u u}^{i i}} \end{matrix}] = = [\begin{matrix} {- - 22 W W}_{ins ins}^{i i} {u u}_{n no}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ni ni} (({d d}_{krig krig}^{ni ni} - - {d d}_{n no}^{i i})) \\ {- - 22 W W}_{ins ins}^{i i} {u u}_{e e}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - {22 W W}_{krig krig}^{ei ei} (({d d}_{krig krig}^{ei ei} - - {d d}_{e e}^{i i})) \\ {- - 22 W W}_{ins ins}^{i i} {u u}_{u u}^{i i} (({d d}_{ins ins}^{i i} - - {u u}_{n no}^{i i} {d d}_{n no}^{i i} - - {u u}_{e e}^{i i} {d d}_{e e}^{i i} - - {u u}_{u u}^{i i} {d d}_{u u}^{i i})) - - 22 {W W}_{krig krig}^{ui ui} (({d d}_{krig krig}^{ui ui} - - {d d}_{u u}^{i i})) \end{matrix}] = = 00

When the coefficient matrix of the equation is not zero, a set of unique solutions is obtained, which is the maximum a posteriori estimation of the three-dimensional deformation field with high spatial resolution;

Then, Kriging interpolation and encryption are performed on the above-mentioned 3D deformation field on the high time frequency GPS data in the time domain, so that the quasi-GPS results are encrypted into a quasi-GPS time series in the time domain, and 3D deformation with high temporal resolution is obtained.