[go: up one dir, main page]

CN115629384A - A correction method and related equipment for timing InSAR errors - Google Patents

A correction method and related equipment for timing InSAR errors Download PDF

Info

Publication number
CN115629384A
CN115629384A CN202211568937.XA CN202211568937A CN115629384A CN 115629384 A CN115629384 A CN 115629384A CN 202211568937 A CN202211568937 A CN 202211568937A CN 115629384 A CN115629384 A CN 115629384A
Authority
CN
China
Prior art keywords
insar
error
equation
deformation
phase
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202211568937.XA
Other languages
Chinese (zh)
Other versions
CN115629384B (en
Inventor
刘计洪
胡俊
李志伟
朱建军
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN202211568937.XA priority Critical patent/CN115629384B/en
Publication of CN115629384A publication Critical patent/CN115629384A/en
Application granted granted Critical
Publication of CN115629384B publication Critical patent/CN115629384B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
    • G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
    • G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明提供了一种时序InSAR误差的改正方法及相关设备，改正方法，包括：步骤1，获取研究区域的InSAR时序相位；步骤2，基于多源信号的时空特征，在InSAR时序相位中对空间趋势误差、地形相关误差和形变信号建立联合模型，并对联合模型中的模型参数进行求解，得到空间趋势误差和地形相关误差；步骤3，将空间趋势误差和地形相关误差在InSAR时序相位中减去，得到改正后的InSAR时序形变；有效避免了各类误差单独处理时易受形变或其他误差信号影响的问题，实现了大气延迟和轨道误差的准确估计，显著提高了时序InSAR地表形变测量的精度及可靠性。

Figure 202211568937

The present invention provides a correction method and related equipment for time series InSAR errors. The correction method includes: step 1, obtaining the InSAR time series phase of the research area; The trend error, terrain-related error and deformation signal are used to establish a joint model, and the model parameters in the joint model are solved to obtain the spatial trend error and terrain-related error; step 3, the spatial trend error and terrain-related error are subtracted from the InSAR time series phase To get the corrected InSAR time-series deformation; effectively avoid the problem that various errors are easily affected by deformation or other error signals when processed separately, realize accurate estimation of atmospheric delay and orbit error, and significantly improve the accuracy of time-series InSAR surface deformation measurement. precision and reliability.

Figure 202211568937

Description

一种时序InSAR误差的改正方法及相关设备A correction method and related equipment for timing InSAR errors

技术领域technical field

本发明涉及地表形变测量技术领域，特别涉及一种时序InSAR误差的改正方法及相关及设备。The invention relates to the technical field of surface deformation measurement, in particular to a method for correcting time series InSAR errors and related equipment.

背景技术Background technique

时序合成孔径雷达干涉测量（Interferometric Synthetic Aperture Radar，SAR，InSAR）技术通过分析时序SAR影像可获取研究区域长时序、高空间分辨率的地表形变结果，为相关地质灾害解译、分析和预防提供了重要的数据支撑。然而，受大气延迟（包括湍流大气和分层大气）和轨道误差等因素影响，时序InSAR技术的形变测量精度难以得到有效的保障，进而影响后续的应用分析。在准确的外部气象数据（如温度、气压等）或者其他大地测量资料（如（Global Navigation Satellite System，GNSS）全球卫星导航系统和水准数据）可用的情况下，可以有效地削弱InSAR观测值中的大气误差和轨道误差。但是，一般情况下，这些外部数据难以获取，仅能基于形变、大气延迟、轨道误差等信号在InSAR观测值中的不同时空特性来进行误差改正。例如，轨道误差在空间上可以利用二阶多项式进行建模和改正，分层大气可以通过与地形进行相关分析进行建模和改正，湍流大气可以通过时空滤波进行抑制。尽管现有研究针对各类误差的改正和抑制已开展了相关工作，但是大多数研究是在假定其他误差影响较小的前提下，针对某一类误差进行建模和改正。但是，在实际情况中，轨道误差、湍流大气和分层大气往往是同时存在的，若仅对某一类误差进行建模和改正，其改正精度勃然会受到其他误差信号的影响，最终降低时序InSAR形变测量的可靠性。Time-series synthetic aperture radar interferometry (Interferometric Synthetic Aperture Radar, SAR, InSAR) technology can obtain long-term, high-spatial-resolution surface deformation results in the study area by analyzing time-series SAR images, which provides a basis for the interpretation, analysis and prevention of related geological disasters. Important data support. However, due to factors such as atmospheric delay (including turbulent atmosphere and stratified atmosphere) and orbit errors, it is difficult to effectively guarantee the deformation measurement accuracy of time-series InSAR technology, which will affect subsequent application analysis. When accurate external meteorological data (such as temperature, air pressure, etc.) or other geodetic data (such as (Global Navigation Satellite System, GNSS) global satellite navigation system and leveling data) are available, it can effectively weaken the Atmospheric and orbital errors. However, in general, these external data are difficult to obtain, and error correction can only be performed based on the different temporal and spatial characteristics of signals such as deformation, atmospheric delay, and orbital error in InSAR observations. For example, orbital errors can be modeled and corrected spatially using second-order polynomials, stratified atmospheres can be modeled and corrected through correlation analysis with terrain, and turbulent atmospheres can be suppressed through spatiotemporal filtering. Although the existing research has carried out relevant work on the correction and suppression of various types of errors, most of the research is based on the modeling and correction of a certain type of error under the assumption that other errors have less influence. However, in actual situations, orbital errors, turbulent atmosphere, and stratified atmosphere often exist at the same time. If only one type of error is modeled and corrected, the correction accuracy will be affected by other error signals, which will eventually reduce the timing. Reliability of InSAR deformation measurements.

发明内容Contents of the invention

本发明提供了一种时序InSAR误差的改正方法及相关设备，其目的是为了避免各类误差单独处理时易受形变或其他误差信号影响的问题，提升时序InSAR形变测量的精度和可靠性。The present invention provides a method for correcting time-series InSAR errors and related equipment, the purpose of which is to avoid the problem that various errors are easily affected by deformation or other error signals when processed separately, and improve the accuracy and reliability of time-series InSAR deformation measurement.

为了达到上述目的，本发明提供了一种时序InSAR误差的改正方法，包括：In order to achieve the above object, the present invention provides a method for correcting timing InSAR errors, including:

步骤1，获取InSAR干涉图的InSAR时序相位；Step 1, obtain the InSAR timing phase of the InSAR interferogram;

步骤2，基于多源信号的时空特征，在InSAR时序相位中对空间趋势误差、地形相关误差和形变信号建立联合模型，并对联合模型中的模型参数进行求解，得到空间趋势误差和地形相关误差；Step 2. Based on the spatio-temporal characteristics of multi-source signals, a joint model is established for the spatial trend error, terrain-related error and deformation signal in the InSAR time series phase, and the model parameters in the joint model are solved to obtain the spatial trend error and terrain-related error ;

步骤3，将空间趋势误差和地形相关误差在InSAR时序相位中减去，得到改正后的InSAR时序形变。Step 3: Subtract the spatial trend error and terrain-related error from the InSAR time series phase to obtain the corrected InSAR time series deformation.

进一步来说，步骤2包括：Further, step 2 includes:

利用一阶多项式对空间趋势误差进行建模，得到多项式拟合系数方程；Using the first-order polynomial to model the spatial trend error, the polynomial fitting coefficient equation is obtained;

利用线性模型对地形相关误差进行表征，得到线性方程系数方程；The linear model is used to characterize the terrain-related error, and the coefficient equation of the linear equation is obtained;

利用三阶多项式函数对形变信号构建虚拟观测方程；Construct a virtual observation equation for the deformation signal by using a third-order polynomial function;

联立多项式拟合系数方程、线性方程系数方程和虚拟观测方程，得到联合模型。Simultaneous polynomial fitting coefficient equations, linear equation coefficient equations and virtual observation equations are used to obtain a joint model.

进一步来说，联合模型中的模型参数包括所有像元的时序形变、所有像元处空间趋势误差的多项式拟合系数以及所有像元处地形相关误差的线性方程系数。Furthermore, the model parameters in the joint model include the time-series deformation of all pixels, the polynomial fitting coefficients of spatial trend errors at all pixels, and the linear equation coefficients of terrain-related errors at all pixels.

进一步来说，所述步骤2还包括：Further, the step 2 also includes:

从所有像元中任意选择一个像元作为目标像元；Randomly select a cell from all the cells as the target cell;

在时刻

Figure 20214DEST_PATH_IMAGE001
，以目标像元
Figure 468513DEST_PATH_IMAGE002
处为中心取大小为
Figure 513830DEST_PATH_IMAGE003
的窗口，建立多项式拟合系数方程为：at the moment
Figure 20214DEST_PATH_IMAGE001
, with the target pixel
Figure 468513DEST_PATH_IMAGE002
Take the center as the center and take the size as
Figure 513830DEST_PATH_IMAGE003
window, the polynomial fitting coefficient equation is established as:

Figure 174618DEST_PATH_IMAGE005
Figure 174618DEST_PATH_IMAGE005

其中，

Figure 520149DEST_PATH_IMAGE006
为窗口内每个像元处的InSAR相位观测值，
Figure 608190DEST_PATH_IMAGE007
为每一个点处的形变相位，
Figure 406382DEST_PATH_IMAGE008
为
Figure 870862DEST_PATH_IMAGE010
区间内的整数，
Figure 805319DEST_PATH_IMAGE012
为地形相关误差相位，
Figure 329842DEST_PATH_IMAGE014
、
Figure 615330DEST_PATH_IMAGE015
、
Figure 617921DEST_PATH_IMAGE016
为多项式拟合系数；in,
Figure 520149DEST_PATH_IMAGE006
is the InSAR phase observation value at each pixel in the window,
Figure 608190DEST_PATH_IMAGE007
is the deformation phase at each point,
Figure 406382DEST_PATH_IMAGE008
for
Figure 870862DEST_PATH_IMAGE010
integers in the interval,
Figure 805319DEST_PATH_IMAGE012
is the terrain-related error phase,
Figure 329842DEST_PATH_IMAGE014
,
Figure 615330DEST_PATH_IMAGE015
,
Figure 617921DEST_PATH_IMAGE016
is the polynomial fitting coefficient;

将上式中的观测值向量

Figure 672464DEST_PATH_IMAGE018
纳入总的观测值向量
Figure 367888DEST_PATH_IMAGE019
中，则相应的系数矩阵
Figure 875093DEST_PATH_IMAGE020
会增加
Figure 681375DEST_PATH_IMAGE022
行，每一行与上式中观测值向量
Figure 590425DEST_PATH_IMAGE017
的元素相对应，每一行中大部分元素为0，只有与模型参数向量
Figure 191170DEST_PATH_IMAGE023
中元素相对应的列不为0，其数值为上式中与系数矩阵
Figure 454180DEST_PATH_IMAGE024
相应行对应的元素。The observation vector in the above formula
Figure 672464DEST_PATH_IMAGE018
Include total observations vector
Figure 367888DEST_PATH_IMAGE019
, then the corresponding coefficient matrix
Figure 875093DEST_PATH_IMAGE020
will increase
Figure 681375DEST_PATH_IMAGE022
row, each row is related to the observation vector in the above formula
Figure 590425DEST_PATH_IMAGE017
Corresponding to the elements, most of the elements in each row are 0, only the model parameter vector
Figure 191170DEST_PATH_IMAGE023
The column corresponding to the element in is not 0, and its value is the coefficient matrix in the above formula
Figure 454180DEST_PATH_IMAGE024
The element corresponding to the corresponding row.

进一步来说，所述步骤2还包括：Further, the step 2 also includes:

从所有像元中任意选择一个像元作为目标像元；Randomly select a cell from all the cells as the target cell;

在时刻

Figure 798574DEST_PATH_IMAGE001
，以目标像元
Figure 562130DEST_PATH_IMAGE002
为中心取大小为
Figure 599357DEST_PATH_IMAGE025
的窗口，建立线性方程系数方程：at the moment
Figure 798574DEST_PATH_IMAGE001
, with the target pixel
Figure 562130DEST_PATH_IMAGE002
Take the size of the center as
Figure 599357DEST_PATH_IMAGE025
window, set up the linear equation coefficient equation:

Figure 81153DEST_PATH_IMAGE027
Figure 81153DEST_PATH_IMAGE027

其中，

Figure 229238DEST_PATH_IMAGE006
为窗口内每个像元处的InSAR相位观测值，
Figure 847301DEST_PATH_IMAGE007
为形变相位，
Figure 321008DEST_PATH_IMAGE028
为趋势相位，
Figure 24522DEST_PATH_IMAGE029
为窗口区间内的整数，
Figure 710718DEST_PATH_IMAGE030
为像元
Figure 714446DEST_PATH_IMAGE002
与像元
Figure 93475DEST_PATH_IMAGE031
之间的高程之差，
Figure 284285DEST_PATH_IMAGE032
、
Figure 39751DEST_PATH_IMAGE033
为线性方程系数；in,
Figure 229238DEST_PATH_IMAGE006
is the InSAR phase observation value at each pixel in the window,
Figure 847301DEST_PATH_IMAGE007
is the deformation phase,
Figure 321008DEST_PATH_IMAGE028
is the trend phase,
Figure 24522DEST_PATH_IMAGE029
is an integer in the window interval,
Figure 710718DEST_PATH_IMAGE030
for the pixel
Figure 714446DEST_PATH_IMAGE002
with the pixel
Figure 93475DEST_PATH_IMAGE031
The difference in elevation between,
Figure 284285DEST_PATH_IMAGE032
,
Figure 39751DEST_PATH_IMAGE033
is the linear equation coefficient;

将上式中的观测值向量

Figure 366827DEST_PATH_IMAGE017
纳入总的观测值向量
Figure 916757DEST_PATH_IMAGE019
中，则相应的系数矩阵
Figure 594863DEST_PATH_IMAGE020
会增加
Figure 154021DEST_PATH_IMAGE021
行，每一行与上式中观测值向量
Figure 601182DEST_PATH_IMAGE017
的元素相对应，每一行中大部分元素为0，只有与模型参数向量
Figure 56435DEST_PATH_IMAGE023
中元素相对应的列不为0，其数值为上式中系数矩阵
Figure 221837DEST_PATH_IMAGE024
相应行对应的元素。The observation vector in the above formula
Figure 366827DEST_PATH_IMAGE017
Include total observations vector
Figure 916757DEST_PATH_IMAGE019
, then the corresponding coefficient matrix
Figure 594863DEST_PATH_IMAGE020
will increase
Figure 154021DEST_PATH_IMAGE021
row, each row is related to the observation vector in the above formula
Figure 601182DEST_PATH_IMAGE017
Corresponding to the elements, most of the elements in each row are 0, only the model parameter vector
Figure 56435DEST_PATH_IMAGE023
The column corresponding to the element in is not 0, and its value is the coefficient matrix in the above formula
Figure 221837DEST_PATH_IMAGE024
The element corresponding to the corresponding row.

进一步来说，所述步骤2还包括：Further, the step 2 also includes:

从所有像元中任意选择一个像元作为目标像元；Randomly select a cell from all the cells as the target cell;

在时刻

Figure 319106DEST_PATH_IMAGE001
，以目标像元
Figure 886353DEST_PATH_IMAGE002
为中心取大小为的时间窗口，构建虚拟观测方程为：at the moment
Figure 319106DEST_PATH_IMAGE001
, with the target pixel
Figure 886353DEST_PATH_IMAGE002
Taking a time window of size as the center, the virtual observation equation is constructed as:

Figure 512507DEST_PATH_IMAGE034
Figure 512507DEST_PATH_IMAGE034

其中，

Figure 165205DEST_PATH_IMAGE035
为中间模型参数变量，
Figure 63235DEST_PATH_IMAGE036
与
Figure 219410DEST_PATH_IMAGE037
之差等于0，
Figure 282044DEST_PATH_IMAGE038
为时序形变，
Figure 749934DEST_PATH_IMAGE039
为时序形变
Figure 454585DEST_PATH_IMAGE040
的三次多项式拟合值，
Figure 996425DEST_PATH_IMAGE041
、
Figure 229960DEST_PATH_IMAGE042
、
Figure 591671DEST_PATH_IMAGE043
分别为第
Figure 834434DEST_PATH_IMAGE044
景SAR影像与第
Figure 965201DEST_PATH_IMAGE045
景SAR影像时刻之间的时间差、第
Figure 369637DEST_PATH_IMAGE046
景SAR影像与第
Figure 218645DEST_PATH_IMAGE045
景SAR影像时刻之间的时间差、第
Figure 999519DEST_PATH_IMAGE047
景SAR影像与第
Figure 250372DEST_PATH_IMAGE045
景SAR影像时刻之间的时间差，
Figure 825709DEST_PATH_IMAGE048
。in,
Figure 165205DEST_PATH_IMAGE035
is the intermediate model parameter variable,
Figure 63235DEST_PATH_IMAGE036
and
Figure 219410DEST_PATH_IMAGE037
The difference is equal to 0,
Figure 282044DEST_PATH_IMAGE038
is the time series deformation,
Figure 749934DEST_PATH_IMAGE039
time series deformation
Figure 454585DEST_PATH_IMAGE040
The cubic polynomial fit value of ,
Figure 996425DEST_PATH_IMAGE041
,
Figure 229960DEST_PATH_IMAGE042
,
Figure 591671DEST_PATH_IMAGE043
respectively
Figure 834434DEST_PATH_IMAGE044
SAR image and the first
Figure 965201DEST_PATH_IMAGE045
The time difference between scene SAR image moments, the first
Figure 369637DEST_PATH_IMAGE046
SAR image and the first
Figure 218645DEST_PATH_IMAGE045
The time difference between scene SAR image moments, the first
Figure 999519DEST_PATH_IMAGE047
SAR image and the first
Figure 250372DEST_PATH_IMAGE045
The time difference between the scene SAR image moments,
Figure 825709DEST_PATH_IMAGE048
.

进一步来说，所述步骤2还包括：Further, the step 2 also includes:

利用先验信息建立额外的虚拟观测方程为：Using the prior information to establish an additional virtual observation equation is:

Figure 162013DEST_PATH_IMAGE049
Figure 162013DEST_PATH_IMAGE049

将该虚拟观测方程加入总的观测值向量

Figure 480999DEST_PATH_IMAGE019
中，并在对应的系数矩阵
Figure 851937DEST_PATH_IMAGE020
中加入一行。Add this dummy observation equation to the total observation vector
Figure 480999DEST_PATH_IMAGE019
, and in the corresponding coefficient matrix
Figure 851937DEST_PATH_IMAGE020
Add a line to the .

进一步来说，联立多项式拟合系数方程、线性方程系数方程和虚拟观测方程，得到联合模型，包括：Furthermore, the simultaneous polynomial fitting coefficient equation, linear equation coefficient equation and virtual observation equation obtain a joint model, including:

基于多项式拟合系数方程、线性方程系数方程、虚拟观测方程以及额外的虚拟观测方程，建立InSAR时序相位中空间趋势误差、地形相关误差和形变信号的联合模型为：Based on polynomial fitting coefficient equations, linear equation coefficient equations, virtual observation equations and additional virtual observation equations, the joint model of spatial trend error, terrain-related error and deformation signal in InSAR time series phase is established as follows:

Figure 598176DEST_PATH_IMAGE050
Figure 598176DEST_PATH_IMAGE050

其中，

Figure 156197DEST_PATH_IMAGE019
为总的观测值向量，
Figure 278873DEST_PATH_IMAGE020
为系数矩阵，
Figure 507248DEST_PATH_IMAGE051
为模型参数。in,
Figure 156197DEST_PATH_IMAGE019
is the total observation vector,
Figure 278873DEST_PATH_IMAGE020
is the coefficient matrix,
Figure 507248DEST_PATH_IMAGE051
is the model parameter.

本发明还提供了一种计算机可读存储介质，用于存储计算机程序，通过执行计算机程序，用于实现上述的时序InSAR误差的改正方法。The present invention also provides a computer-readable storage medium, which is used to store a computer program, and is used to implement the above-mentioned method for correcting timing InSAR errors by executing the computer program.

本发明还提供了一种时序InSAR误差的改正设备，用于实现上述的时序InSAR误差的改正方法，包括：The present invention also provides a correction device for timing InSAR errors, which is used to implement the above correction method for timing InSAR errors, including:

存储器和处理器；memory and processor;

存储器用于储存计算机程序；Memory is used to store computer programs;

处理器用于执行存储器存储的计算机程序。The processor is used to execute the computer program stored in the memory.

本发明的上述方案有如下的有益效果：Said scheme of the present invention has following beneficial effect:

本发明基于形变、大气延迟和轨道误差等多源信号的时空特征，构建了InSAR时序相位中空间趋势误差、地形相关误差和形变信号的联合模型，有效避免了各类误差单独处理时易受形变或其他误差信号影响的问题，实现了大气延迟和轨道误差的准确估计，显著提高了时序InSAR地表形变测量的精度及可靠性。Based on the spatio-temporal characteristics of multi-source signals such as deformation, atmospheric delay and orbit error, the present invention constructs a joint model of spatial trend error, terrain-related error and deformation signal in InSAR time series phase, effectively avoiding various types of errors that are susceptible to deformation when processed separately Or other problems affected by error signals, the accurate estimation of atmospheric delay and orbit error has been realized, and the accuracy and reliability of time-series InSAR surface deformation measurement have been significantly improved.

本发明的其它有益效果将在随后的具体实施方式部分予以详细说明。Other beneficial effects of the present invention will be described in detail in the following specific embodiments.

附图说明Description of drawings

图1为本发明实施例的流程原理图；Fig. 1 is the schematic flow chart of the embodiment of the present invention;

图2为本发明实施例采用不同方法得到的漏斗中心处InSAR时序形变结果对比图。Fig. 2 is a comparison diagram of InSAR time-series deformation results at the center of the funnel obtained by using different methods according to the embodiment of the present invention.

具体实施方式Detailed ways

为使本发明要解决的技术问题、技术方案和优点更加清楚，下面将结合附图及具体实施例进行详细描述。显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。In order to make the technical problems, technical solutions and advantages to be solved by the present invention clearer, the following will describe in detail with reference to the drawings and specific embodiments. Apparently, the described embodiments are some, but not all, embodiments of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by persons of ordinary skill in the art without making creative efforts belong to the protection scope of the present invention.

在本发明的描述中，需要说明的是，术语“中心”、“上”、“下”、“左”、“右”、“竖直”、“水平”、“内”、“外”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。此外，术语“第一”、“第二”、“第三”仅用于描述目的，而不能理解为指示或暗示相对重要性。In the description of the present invention, it should be noted that the terms "center", "upper", "lower", "left", "right", "vertical", "horizontal", "inner", "outer" etc. The indicated orientation or positional relationship is based on the orientation or positional relationship shown in the drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the referred device or element must have a specific orientation, or in a specific orientation. construction and operation, therefore, should not be construed as limiting the invention. In addition, the terms "first", "second", and "third" are used for descriptive purposes only, and should not be construed as indicating or implying relative importance.

在本发明的描述中，需要说明的是，除非另有明确的规定和限定，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是锁定连接，也可以是可拆卸连接，或一体地连接；可以是机械连接，也可以是电连接；可以是直接相连，也可以通过中间媒介间接相连，可以是两个元件内部的连通。对于本领域的普通技术人员而言，可以具体情况理解上述术语在本发明中的具体含义。In the description of the present invention, it should be noted that unless otherwise specified and limited, the terms "installation", "connection" and "connection" should be understood in a broad sense, for example, it can be a locking connection or a detachable connection. Connected, or integrally connected; it can be mechanically connected or electrically connected; it can be directly connected or indirectly connected through an intermediary, and it can be the internal communication of two components. Those of ordinary skill in the art can understand the specific meanings of the above terms in the present invention in specific situations.

此外，下面所描述的本发明不同实施方式中所涉及的技术特征只要彼此之间未构成冲突就可以相互结合。In addition, the technical features involved in the different embodiments of the present invention described below may be combined with each other as long as there is no conflict with each other.

本发明针对现有的问题，提供了一种时序InSAR误差的改正方法及相关设备。Aiming at the existing problems, the present invention provides a time series InSAR error correction method and related equipment.

具体来说，本发明实施例是对时序InSAR大气延迟与轨道误差进行联合改正，轨道误差在空间上均为认为是低频信号，同时湍流大气误差在小空间范围内，如1×1km2也可认为是低频信号，因此在本发明实施例中将轨道误差和湍流大气统称为趋势误差。对于某一像素点处的空间趋势误差而言，可利用空间一定范围内，如5×5像素的观测值进行多项式拟合建模。分层大气主要与局部地形相关，因此被称为地形相关误差，对于某一像素点处的地形相关误差而言，可对空间一定范围内的InSAR观测值和地形进行线性拟合建模，此外，排除地震等突发性事件，一般情况下，地表形变在时间上是一个缓慢变化的过程，在一定时间范围内，可认为时序地表形变满足与时间相关的三次方程，将此假设作为先验约束条件，在InSAR时序相位中对大气延迟误差、轨道误差和形变信号进行时空统一建模，即可有效避免各类误差单独处理时易受到形变或其他误差信号影响的问题，进而提升时序InSAR形变测量的精度和可靠性。Specifically, the embodiment of the present invention is to jointly correct the time-series InSAR atmospheric delay and orbital error. The orbital error is regarded as a low-frequency signal in space, and the turbulent atmospheric error is within a small spatial range, such as 1 ×1km2. It is considered to be a low-frequency signal, so in the embodiment of the present invention, the orbit error and the turbulent atmosphere are collectively referred to as the trend error. For the spatial trend error at a certain pixel point, polynomial fitting modeling can be carried out by using observations within a certain range of space, such as 5×5 pixels. The layered atmosphere is mainly related to local terrain, so it is called terrain-related error. For the terrain-related error at a certain pixel point, linear fitting modeling can be performed on InSAR observations and terrain within a certain range of space. In addition , excluding sudden events such as earthquakes, in general, surface deformation is a slowly changing process in time, within a certain time range, it can be considered that the time-series surface deformation satisfies the time-related cubic equation, and this assumption is taken as a priori Constraint conditions, the unified modeling of atmospheric delay error, orbit error and deformation signal in InSAR time series phase can effectively avoid the problem that various errors are easily affected by deformation or other error signals when processed separately, and then improve the time series InSAR deformation Measurement accuracy and reliability.

如图1所示，本发明的实施例提供了一种时序InSAR误差的改正方法，包括：As shown in Figure 1, the embodiment of the present invention provides a method for correcting timing InSAR errors, including:

步骤1，获取InSAR干涉图的InSAR时序相位；Step 1, obtain the InSAR timing phase of the InSAR interferogram;

步骤2，基于多源信号的时空特征，在InSAR时序相位中对空间趋势误差、地形相关误差和形变信号建立联合模型，并对联合模型中的模型参数进行求解，得到空间趋势误差和地形相关误差；Step 2. Based on the spatio-temporal characteristics of multi-source signals, a joint model is established for the spatial trend error, terrain-related error and deformation signal in the InSAR time series phase, and the model parameters in the joint model are solved to obtain the spatial trend error and terrain-related error ;

步骤3，将空间趋势误差和地形相关误差在InSAR时序相位中减去，得到改正后的InSAR时序形变。Step 3: Subtract the spatial trend error and terrain-related error from the InSAR time series phase to obtain the corrected InSAR time series deformation.

具体来说，本发明实施例基于SBAS（Satellite-Based Augmentation System）星基增强系统或SqueeSAR分布永久散射体卫星雷达监测技术方法，无需改正InSAR干涉图中的轨道误差和大气延迟误差，直到得到InSAR时序相位，后续的空间趋势误差和地形相关误差改正则是针对InSAR时序相位进行的。Specifically, the embodiment of the present invention is based on the SBAS (Satellite-Based Augmentation System) satellite-based augmentation system or the SqueeSAR distributed permanent scatterer satellite radar monitoring technology method, without correcting the orbit error and atmospheric delay error in the InSAR interferogram until the InSAR The timing phase, subsequent spatial trend error and terrain-related error corrections are performed for the InSAR timing phase.

InSAR干涉图中存在轨道误差和大气延迟误差，其根本是因为组成干涉图的SAR影像中包含相应的轨道误差和大气延迟误差，因此，基于不进行轨道误差和大气延迟误差改正的InSAR干涉图所得到的InSAR时序相位中，不仅包含时序形变还包含每一个时刻所所对应的轨道误差和大气延迟误差因此，直接针对InSAR时序相位进行相应的误差改正是更为合理的。同时，相较于InSAR干涉图而言，InSAR时序相位的数据量可大大降低，进而有利于提高数据处理效率。There are orbit errors and atmospheric delay errors in the InSAR interferogram, which is basically because the SAR images that make up the interferogram contain the corresponding orbit errors and atmospheric delay errors. Therefore, based on the InSAR interferogram without orbit error and atmospheric delay error correction The obtained InSAR timing phase includes not only the timing deformation but also the orbit error and atmospheric delay error corresponding to each moment. Therefore, it is more reasonable to directly correct the corresponding error for the InSAR timing phase. At the same time, compared with the InSAR interferogram, the data volume of the InSAR timing phase can be greatly reduced, which in turn helps to improve the data processing efficiency.

具体来说，本发明实施例的步骤2，主要考虑三个方面的函数模型建立：Specifically, step 2 of the embodiment of the present invention mainly considers three aspects of function model establishment:

1对于空间趋势误差而言，在1km×1km的小窗口范围内，可利用与位置相关的一阶多项式进行建模，得到多项式拟合系数方程；1 For the spatial trend error, within the small window range of 1km×1km, the first-order polynomial related to the position can be used for modeling, and the polynomial fitting coefficient equation can be obtained;

2对于地形相关误差而言，在一定的小窗口范围内，可利用与地形相关的线性模型进行表征，得到线性方程系数方程；2 For terrain-related errors, within a certain small window range, the linear model related to terrain can be used to characterize and obtain the linear equation coefficient equation;

3对于形变信号而言，其在一定的空间范围内，往往与空间趋势误差表现为相似的空间特征，不同的是，形变信号在一定的时间窗口内也是相关的，因此可以通过与时间相关的三阶多项式函数，建立与时序形变相关的虚拟观测方程，进而同1和2中的函数模型联立，实现空间趋势误差、地形相关误差和时序形变信号的时空同一建模。3 For the deformation signal, within a certain spatial range, it often exhibits similar spatial characteristics to the spatial trend error, but the difference is that the deformation signal is also correlated within a certain time window, so it can The third-order polynomial function establishes a virtual observation equation related to time-series deformation, and then combines it with the function model in 1 and 2 to realize the same time-space modeling of spatial trend error, terrain-related error and time-series deformation signal.

在本发明实施例中，假设现有N个时刻的InSAR时序相位，每个时刻的相位为I×J大小的矩阵，则本发明实施例的联合模型中需要求解的模型参数包括：所有像元的时序形变

Figure 424389DEST_PATH_IMAGE052
，所有像元处空间趋势误差的多项式拟合系数
Figure 735284DEST_PATH_IMAGE053
，所有像元处地形相关误差的线性方程系数
Figure 396073DEST_PATH_IMAGE054
，其中
Figure 476024DEST_PATH_IMAGE055
表示第
Figure 829645DEST_PATH_IMAGE056
个时刻。In the embodiment of the present invention, assuming that there are InSAR time-series phases at N moments, and the phase at each moment is a matrix of size I×J , the model parameters that need to be solved in the joint model of the embodiment of the present invention include: all pixels time series deformation
Figure 424389DEST_PATH_IMAGE052
, the polynomial fit coefficient for the spatial trend error at all cells
Figure 735284DEST_PATH_IMAGE053
, the coefficients of the linear equation for the terrain-dependent error at all cells
Figure 396073DEST_PATH_IMAGE054
,in
Figure 476024DEST_PATH_IMAGE055
Indicates the first
Figure 829645DEST_PATH_IMAGE056
moment.

上述模型参数是通过建立一个函数模型进行求解的，因此所建立的函数模型中模型参数

Figure 362258DEST_PATH_IMAGE057
的大小为
Figure 826737DEST_PATH_IMAGE058
（一般选第一景影像为参考影像），相应系数矩阵
Figure 495616DEST_PATH_IMAGE020
的列数与模型参数向量
Figure 285717DEST_PATH_IMAGE059
的长度一致。此时，若要求解模型参数
Figure 305626DEST_PATH_IMAGE057
的数值，则关键在于构建系数矩阵
Figure 308217DEST_PATH_IMAGE020
和观测值向量
Figure 97181DEST_PATH_IMAGE019
。The above model parameters are solved by establishing a function model, so the model parameters in the established function model
Figure 362258DEST_PATH_IMAGE057
is of size
Figure 826737DEST_PATH_IMAGE058
(Generally, the first scene image is selected as the reference image), and the corresponding coefficient matrix
Figure 495616DEST_PATH_IMAGE020
The number of columns and model parameter vector
Figure 285717DEST_PATH_IMAGE059
of the same length. At this point, if the model parameters are to be solved
Figure 305626DEST_PATH_IMAGE057
value, the key is to construct the coefficient matrix
Figure 308217DEST_PATH_IMAGE020
and the observation vector
Figure 97181DEST_PATH_IMAGE019
.

具体来说，本发明实施例的步骤2还包括：Specifically, step 2 of the embodiment of the present invention also includes:

从所有像元中任意选择一个像元作为目标像元为例，介绍建立与该时刻该像元模型参数相关的系数矩阵和观测值向量的具体过程。Randomly select a pixel from all the pixels as the target pixel, and introduce the specific process of establishing the coefficient matrix and observation value vector related to the model parameters of the pixel at that moment.

首先，建立空间趋势误差的多项式拟合系数方程，在

Figure 58184DEST_PATH_IMAGE060
时刻、以目标像元
Figure 299810DEST_PATH_IMAGE061
为中心取大小为
Figure 106091DEST_PATH_IMAGE062
的窗口，则窗口内的InSAR相位可认为是形变相位、空间趋势误差相位和地形相关误差相位之和，其中，窗口内的模型参数包括每一个像元的形变相位
Figure 15142DEST_PATH_IMAGE063
、地形相关误差相位
Figure 881467DEST_PATH_IMAGE064
以及窗口内的空间趋势误差相位的多项式拟合系数
Figure 875967DEST_PATH_IMAGE065
，其中
Figure 220361DEST_PATH_IMAGE066
为
Figure 983918DEST_PATH_IMAGE067
区间内的整数，因此，对于窗口内每一个像元处的InSAR相位观测值
Figure 755565DEST_PATH_IMAGE068
可写为：First, establish the polynomial fitting coefficient equation of the spatial trend error, in
Figure 58184DEST_PATH_IMAGE060
time, target pixel
Figure 299810DEST_PATH_IMAGE061
Take the size of the center as
Figure 106091DEST_PATH_IMAGE062
window, the InSAR phase in the window can be considered as the sum of the deformation phase, spatial trend error phase and terrain-related error phase, where the model parameters in the window include the deformation phase of each pixel
Figure 15142DEST_PATH_IMAGE063
, terrain-related error phase
Figure 881467DEST_PATH_IMAGE064
and the polynomial fit coefficients for the phase of the spatial trend error within the window
Figure 875967DEST_PATH_IMAGE065
,in
Figure 220361DEST_PATH_IMAGE066
for
Figure 983918DEST_PATH_IMAGE067
An integer in the interval, therefore, for the InSAR phase observations at each pixel in the window
Figure 755565DEST_PATH_IMAGE068
can be written as:

Figure 502941DEST_PATH_IMAGE069
（1）
Figure 502941DEST_PATH_IMAGE069
(1)

考虑到

Figure 651025DEST_PATH_IMAGE070
时刻窗口内共有
Figure 269088DEST_PATH_IMAGE071
个InSAR相位观测值，则有：considering
Figure 651025DEST_PATH_IMAGE070
Total time window
Figure 269088DEST_PATH_IMAGE071
InSAR phase observations, then there are:

Figure 477216DEST_PATH_IMAGE072
（2）
Figure 477216DEST_PATH_IMAGE072
(2)

其中，

Figure 443379DEST_PATH_IMAGE073
为窗口内每个像元处的InSAR相位观测值，
Figure 129575DEST_PATH_IMAGE074
为每一个点处的形变相位，
Figure 867724DEST_PATH_IMAGE075
为
Figure 981174DEST_PATH_IMAGE076
区间内的整数，
Figure 171984DEST_PATH_IMAGE077
为地形相关误差相位，
Figure 927450DEST_PATH_IMAGE078
、
Figure 785685DEST_PATH_IMAGE079
、
Figure 335615DEST_PATH_IMAGE080
为多项式拟合系数，T表示矩阵转置。in,
Figure 443379DEST_PATH_IMAGE073
is the InSAR phase observation value at each pixel in the window,
Figure 129575DEST_PATH_IMAGE074
is the deformation phase at each point,
Figure 867724DEST_PATH_IMAGE075
for
Figure 981174DEST_PATH_IMAGE076
integers in the interval,
Figure 171984DEST_PATH_IMAGE077
is the terrain-related error phase,
Figure 927450DEST_PATH_IMAGE078
,
Figure 785685DEST_PATH_IMAGE079
,
Figure 335615DEST_PATH_IMAGE080
is the polynomial fitting coefficient, and T represents the matrix transpose.

将上式（2）中的观测值向量

Figure 13721DEST_PATH_IMAGE081
纳入总的观测值向量
Figure 307299DEST_PATH_IMAGE082
中，则相应的系数矩阵
Figure 754461DEST_PATH_IMAGE083
会增加
Figure 475292DEST_PATH_IMAGE084
行，每一行与上式（2）中观测值向量
Figure 640694DEST_PATH_IMAGE017
的元素相对应，每一行中大部分元素为0，只有与模型参数向量
Figure 472384DEST_PATH_IMAGE023
中元素相对应的列不为0，其数值为上式（2）中与系数矩阵
Figure 39632DEST_PATH_IMAGE024
相应行对应的元素。The observation vector in the above formula (2)
Figure 13721DEST_PATH_IMAGE081
Include total observations vector
Figure 307299DEST_PATH_IMAGE082
, then the corresponding coefficient matrix
Figure 754461DEST_PATH_IMAGE083
will increase
Figure 475292DEST_PATH_IMAGE084
row, each row is the same as the observed value vector in the above formula (2)
Figure 640694DEST_PATH_IMAGE017
Corresponding to the elements, most of the elements in each row are 0, only the model parameter vector
Figure 472384DEST_PATH_IMAGE023
The column corresponding to the element in is not 0, and its value is the coefficient matrix in the above formula (2)
Figure 39632DEST_PATH_IMAGE024
The element corresponding to the corresponding row.

具体来说，本发明实施例的步骤2还包括：Specifically, step 2 of the embodiment of the present invention also includes:

建立地形相关误差的线性拟合系数方程，在

Figure 931364DEST_PATH_IMAGE060
时刻、以目标像元
Figure 584062DEST_PATH_IMAGE061
为中心取大小为
Figure 953864DEST_PATH_IMAGE085
的窗口，则窗口内的InSAR相位同样可认为是形变相位、趋势相位和地形相关相位之和，其中不同的是，窗口内的模型参数包括每一个点处的形变相位
Figure 641197DEST_PATH_IMAGE086
、趋势相位
Figure 703831DEST_PATH_IMAGE087
以及窗口内的地形相关误差的线性拟合系数
Figure 578246DEST_PATH_IMAGE088
，其中，
Figure 17318DEST_PATH_IMAGE089
为
Figure 559158DEST_PATH_IMAGE090
区间内的整数，因此，对于窗口内每一个像元处的InSAR相位观测值
Figure 792693DEST_PATH_IMAGE091
可写为：Establish the linear fitting coefficient equation of the terrain-related error, in
Figure 931364DEST_PATH_IMAGE060
time, target pixel
Figure 584062DEST_PATH_IMAGE061
Take the size of the center as
Figure 953864DEST_PATH_IMAGE085
, the InSAR phase in the window can also be considered as the sum of the deformation phase, trend phase and terrain-related phase. The difference is that the model parameters in the window include the deformation phase at each point
Figure 641197DEST_PATH_IMAGE086
, trend phase
Figure 703831DEST_PATH_IMAGE087
and a linear fit coefficient for the terrain-dependent error within the window
Figure 578246DEST_PATH_IMAGE088
,in,
Figure 17318DEST_PATH_IMAGE089
for
Figure 559158DEST_PATH_IMAGE090
An integer in the interval, therefore, for the InSAR phase observations at each pixel in the window
Figure 792693DEST_PATH_IMAGE091
can be written as:

Figure 154404DEST_PATH_IMAGE092
（3）
Figure 154404DEST_PATH_IMAGE092
(3)

考虑到

Figure 134517DEST_PATH_IMAGE070
时刻窗口内共有
Figure 530863DEST_PATH_IMAGE093
个InSAR相位观测值，则有：considering
Figure 134517DEST_PATH_IMAGE070
Total time window
Figure 530863DEST_PATH_IMAGE093
InSAR phase observations, then there are:

Figure 935300DEST_PATH_IMAGE094
（4）
Figure 935300DEST_PATH_IMAGE094
(4)

其中，

Figure 784307DEST_PATH_IMAGE095
为窗口内每个像元处的InSAR相位观测值，
Figure 299602DEST_PATH_IMAGE096
为形变相位，
Figure 816034DEST_PATH_IMAGE097
为趋势相位，
Figure 391372DEST_PATH_IMAGE098
为窗口区间内的整数，
Figure 727675DEST_PATH_IMAGE099
为像元
Figure 312240DEST_PATH_IMAGE100
与像元
Figure 417600DEST_PATH_IMAGE101
之间的高程之差，
Figure 163839DEST_PATH_IMAGE102
、
Figure 987438DEST_PATH_IMAGE103
为线性方程系数；in,
Figure 784307DEST_PATH_IMAGE095
is the InSAR phase observation value at each pixel in the window,
Figure 299602DEST_PATH_IMAGE096
is the deformation phase,
Figure 816034DEST_PATH_IMAGE097
is the trend phase,
Figure 391372DEST_PATH_IMAGE098
is an integer in the window interval,
Figure 727675DEST_PATH_IMAGE099
for the pixel
Figure 312240DEST_PATH_IMAGE100
with the pixel
Figure 417600DEST_PATH_IMAGE101
The difference in elevation between,
Figure 163839DEST_PATH_IMAGE102
,
Figure 987438DEST_PATH_IMAGE103
is the linear equation coefficient;

将上式（4）中的观测值向量

Figure 110115DEST_PATH_IMAGE104
纳入总的观测值向量
Figure 69981DEST_PATH_IMAGE105
中，则相应的系数矩阵
Figure 987121DEST_PATH_IMAGE106
会增加
Figure 563596DEST_PATH_IMAGE107
行，每一行与上式（4）中观测值向量
Figure 224385DEST_PATH_IMAGE104
的元素相对应，每一行中大部分元素为0，只有与模型参数向量
Figure 38757DEST_PATH_IMAGE023
中元素相对应的列不为0，其数值为上式（4）中系数矩阵
Figure 657957DEST_PATH_IMAGE024
相应行对应的元素。The observation vector in the above formula (4)
Figure 110115DEST_PATH_IMAGE104
Include total observations vector
Figure 69981DEST_PATH_IMAGE105
, then the corresponding coefficient matrix
Figure 987121DEST_PATH_IMAGE106
will increase
Figure 563596DEST_PATH_IMAGE107
row, each row is the same as the observation vector in the above formula (4)
Figure 224385DEST_PATH_IMAGE104
Corresponding to the elements, most of the elements in each row are 0, only the model parameter vector
Figure 38757DEST_PATH_IMAGE023
The column corresponding to the element in is not 0, and its value is the coefficient matrix in the above formula (4)
Figure 657957DEST_PATH_IMAGE024
The element corresponding to the corresponding row.

需要说明的是，由于空间趋势误差和地形相关误差在空间上的尺度一般是不相同的，比如用于地形相关误差建模的空间尺度往往与实际地形相关，而相比之下空间趋势误差的空间尺度往往范围更大，因此上述公式（2）和（4）建立过程中窗口尺寸模型参数

Figure 190569DEST_PATH_IMAGE108
一般是不同的。此外，从公式（2）和（4）中可以看出，两类方程中的观测值向量
Figure 389470DEST_PATH_IMAGE109
和
Figure 323928DEST_PATH_IMAGE110
中会包含相同的观测值，也就是说对于InSAR时序相位中某一个观测值
Figure 114029DEST_PATH_IMAGE111
可能会在最终的观测值向量
Figure 868358DEST_PATH_IMAGE112
中出现多次。并且，在公式（1）中对窗口内的趋势相位建模时，其空间位置的参考点即为目标像元
Figure 136529DEST_PATH_IMAGE100
，这种情况下多项式系数中的常数项即为该像元处的趋势相位数值，进而目标像元
Figure 188143DEST_PATH_IMAGE100
处的趋势相位可能也会作为模型参数出现在其他像元处建立的方程中，即公式（2）；同理，在公式（3）中对窗口内的地形相关误差建模时，其高程的参考点即为目标像元
Figure 617987DEST_PATH_IMAGE100
处的高程，这种情况下线性方程系数中的常数项即为该像元处的地形相关误差数值，进而目标像元
Figure 125192DEST_PATH_IMAGE113
处的地形相关误差可能也会作为模型参数出现在其他像元处建立的方程中，即公式（4）。从上述分析可以看出，针对不同像元建立的公式（2）和（4），会包含其他像元处的模型参数，进而无法通过逐像元进行每个像元处的模型参数求解，只能建立一个关于所有时刻、所有像元处的所有模型参数的大型函数模型，然后进行模型参数的整体求解，求解出模型参数即为空间趋势误差和地形相关误差。It should be noted that the spatial scales of spatial trend errors and terrain-related errors are generally different. For example, the spatial scale used for terrain-related error modeling is often related to the actual terrain, while the spatial trend error The spatial scale often has a larger range, so the window size model parameters in the above formulas (2) and (4) during the establishment process
Figure 190569DEST_PATH_IMAGE108
Generally different. Furthermore, from equations (2) and (4), it can be seen that the vector of observations in the two types of equations
Figure 389470DEST_PATH_IMAGE109
and
Figure 323928DEST_PATH_IMAGE110
will contain the same observation value, that is to say, for a certain observation value in the InSAR time series phase
Figure 114029DEST_PATH_IMAGE111
may be in the final observation vector
Figure 868358DEST_PATH_IMAGE112
appears multiple times in . And, when modeling the trend phase in the window in formula (1), the reference point of its spatial position is the target pixel
Figure 136529DEST_PATH_IMAGE100
, in this case the constant term in the polynomial coefficient is the trend phase value at the pixel, and then the target pixel
Figure 188143DEST_PATH_IMAGE100
The trend phase at may also appear as a model parameter in the equations established at other pixels, that is, formula (2); similarly, when modeling the terrain-related error in the window in formula (3), the height of its The reference point is the target pixel
Figure 617987DEST_PATH_IMAGE100
In this case, the constant term in the coefficient of the linear equation is the terrain-related error value at the pixel, and then the target pixel
Figure 125192DEST_PATH_IMAGE113
The terrain-related error at may also appear as a model parameter in the equations established at other pixels, that is, formula (4). From the above analysis, it can be seen that the formulas (2) and (4) established for different pixels will include the model parameters at other pixels, so it is impossible to solve the model parameters at each pixel pixel by pixel, only It is possible to establish a large-scale function model for all model parameters at all time points and all pixels, and then perform an overall solution to the model parameters, and the model parameters that are solved are spatial trend errors and terrain-related errors.

具体来说，本发明实施例中的步骤2还包括：Specifically, step 2 in the embodiment of the present invention also includes:

假定时序形变在一定的时间窗口范围内通过三次多项式进行拟合，即为时序形变提供了外部约束条件，用于建立与时序形变相关的虚拟观测方程。一般情况下，虚拟观测方程的形式为“0=系数矩阵×模型参数向量”。对于时序形变而言，模型参数向量即为时序形变，因此构建虚拟观测方程的关键在于如何根据外部约束条件确定相应的系数矩阵。It is assumed that the time-series deformation is fitted by a cubic polynomial within a certain time window, which provides external constraints for the time-series deformation and is used to establish a virtual observation equation related to the time-series deformation. In general, the form of the virtual observation equation is "0 = coefficient matrix × model parameter vector". For the time-series deformation, the model parameter vector is the time-series deformation, so the key to constructing the virtual observation equation is how to determine the corresponding coefficient matrix according to the external constraints.

在时刻

Figure 931474DEST_PATH_IMAGE070
，以目标像元
Figure 840524DEST_PATH_IMAGE113
处为例，来说明构建虚拟观测方程的过程如下：at the moment
Figure 931474DEST_PATH_IMAGE070
, with the target pixel
Figure 840524DEST_PATH_IMAGE113
Take this place as an example to illustrate the process of constructing the virtual observation equation as follows:

首先，以时刻

Figure 706849DEST_PATH_IMAGE070
为中心，取大小为
Figure 435770DEST_PATH_IMAGE114
的时间窗口
Figure 45743DEST_PATH_IMAGE115
，则相应的模型参数（时序形变）向量为
Figure 74879DEST_PATH_IMAGE116
。当
Figure 908843DEST_PATH_IMAGE117
已知时，则有以下关系：First, take the moment
Figure 706849DEST_PATH_IMAGE070
as the center, take the size as
Figure 435770DEST_PATH_IMAGE114
time window
Figure 45743DEST_PATH_IMAGE115
, then the corresponding model parameter (time series deformation) vector is
Figure 74879DEST_PATH_IMAGE116
. when
Figure 908843DEST_PATH_IMAGE117
When known, the following relationship exists:

Figure 452957DEST_PATH_IMAGE118
（5）
Figure 452957DEST_PATH_IMAGE118
(5)

系数矩阵

Figure 601041DEST_PATH_IMAGE119
中，
Figure 219104DEST_PATH_IMAGE120
、
Figure 958390DEST_PATH_IMAGE121
、
Figure 661904DEST_PATH_IMAGE122
分别为第
Figure 879259DEST_PATH_IMAGE123
景SAR影像与第
Figure 351829DEST_PATH_IMAGE124
景SAR影像时刻之间的时间差、第
Figure 733787DEST_PATH_IMAGE125
景SAR影像与第
Figure 190176DEST_PATH_IMAGE124
景SAR影像时刻之间的时间差、第
Figure 680063DEST_PATH_IMAGE126
景SAR影像与第
Figure 272719DEST_PATH_IMAGE124
景SAR影像时刻之间的时间差，三次多项式系数中的常数项即为
Figure 822649DEST_PATH_IMAGE070
时刻所对应的形变。由于系数矩阵
Figure 500755DEST_PATH_IMAGE119
是已知的，则基于最小二乘方法可得coefficient matrix
Figure 601041DEST_PATH_IMAGE119
middle,
Figure 219104DEST_PATH_IMAGE120
,
Figure 958390DEST_PATH_IMAGE121
,
Figure 661904DEST_PATH_IMAGE122
respectively
Figure 879259DEST_PATH_IMAGE123
SAR image and the first
Figure 351829DEST_PATH_IMAGE124
The time difference between scene SAR image moments, the first
Figure 733787DEST_PATH_IMAGE125
SAR image and the first
Figure 190176DEST_PATH_IMAGE124
The time difference between scene SAR image moments, the first
Figure 680063DEST_PATH_IMAGE126
SAR image and the first
Figure 272719DEST_PATH_IMAGE124
The time difference between the scene SAR image moments, the constant term in the cubic polynomial coefficient is
Figure 822649DEST_PATH_IMAGE070
The deformation corresponding to the moment. Since the coefficient matrix
Figure 500755DEST_PATH_IMAGE119
is known, then based on the least squares method, it can be obtained

Figure 794333DEST_PATH_IMAGE127
（6）
Figure 794333DEST_PATH_IMAGE127
(6)

进而可得时序形变

Figure 241495DEST_PATH_IMAGE128
的三次多项式拟合值
Figure 962326DEST_PATH_IMAGE129
为Then the time series deformation can be obtained
Figure 241495DEST_PATH_IMAGE128
The cubic polynomial fit value of
Figure 962326DEST_PATH_IMAGE129
for

Figure 393307DEST_PATH_IMAGE130
（7）
Figure 393307DEST_PATH_IMAGE130
(7)

式中，

Figure 224997DEST_PATH_IMAGE131
和
Figure 526665DEST_PATH_IMAGE132
均为已知量。In the formula,
Figure 224997DEST_PATH_IMAGE131
and
Figure 526665DEST_PATH_IMAGE132
are known quantities.

理论上，

Figure 418398DEST_PATH_IMAGE133
与
Figure 71096DEST_PATH_IMAGE134
之差等于0，则有In theory,
Figure 418398DEST_PATH_IMAGE133
and
Figure 71096DEST_PATH_IMAGE134
The difference is equal to 0, then there is

Figure 706477DEST_PATH_IMAGE135
（8）
Figure 706477DEST_PATH_IMAGE135
(8)

其中，

Figure 128231DEST_PATH_IMAGE136
为中间模型参数变量，
Figure 190865DEST_PATH_IMAGE137
与
Figure 330859DEST_PATH_IMAGE138
之差等于0，
Figure 504352DEST_PATH_IMAGE139
为时序形变，
Figure 780612DEST_PATH_IMAGE140
为时序形变
Figure 14147DEST_PATH_IMAGE141
的三次多项式拟合值，
Figure 641438DEST_PATH_IMAGE120
、
Figure 618621DEST_PATH_IMAGE121
、
Figure 14967DEST_PATH_IMAGE122
分别为第
Figure 419404DEST_PATH_IMAGE123
景SAR影像与第
Figure 2832DEST_PATH_IMAGE124
景SAR影像时刻之间的时间差、第
Figure 69793DEST_PATH_IMAGE125
景SAR影像与第
Figure 320646DEST_PATH_IMAGE124
景SAR影像时刻之间的时间差、第
Figure 895984DEST_PATH_IMAGE126
景SAR影像与第
Figure 966708DEST_PATH_IMAGE124
景SAR影像时刻之间的时间差，
Figure 551273DEST_PATH_IMAGE142
。in,
Figure 128231DEST_PATH_IMAGE136
is the intermediate model parameter variable,
Figure 190865DEST_PATH_IMAGE137
and
Figure 330859DEST_PATH_IMAGE138
The difference is equal to 0,
Figure 504352DEST_PATH_IMAGE139
is the time series deformation,
Figure 780612DEST_PATH_IMAGE140
time series deformation
Figure 14147DEST_PATH_IMAGE141
The cubic polynomial fit value of ,
Figure 641438DEST_PATH_IMAGE120
,
Figure 618621DEST_PATH_IMAGE121
,
Figure 14967DEST_PATH_IMAGE122
respectively
Figure 419404DEST_PATH_IMAGE123
SAR image and the first
Figure 2832DEST_PATH_IMAGE124
The time difference between scene SAR image moments, the first
Figure 69793DEST_PATH_IMAGE125
SAR image and the first
Figure 320646DEST_PATH_IMAGE124
The time difference between scene SAR image moments, the first
Figure 895984DEST_PATH_IMAGE126
SAR image and the first
Figure 966708DEST_PATH_IMAGE124
The time difference between the scene SAR image moments,
Figure 551273DEST_PATH_IMAGE142
.

公式（8）即为本发明方法中针对时刻

Figure 656632DEST_PATH_IMAGE060
目标像元
Figure 402871DEST_PATH_IMAGE113
处构建的虚拟观测方程。对于每个像元的每一个时刻，均可建立公式（8）的虚拟观测方程，此时在观测值向量
Figure 226471DEST_PATH_IMAGE019
中加入一个0值元素，并且在对应的系数矩阵
Figure 349148DEST_PATH_IMAGE020
中加入一行，这一行中对应公式（8）中模型参数的元素即为公式（8）中系数矩阵相应的元素。可以看出，类似空间上的建模过程，时间上的虚拟观测方程中，对于某一时刻而言，也会包含其他时刻对应像元的模型参数，因此也需要对所有时刻的模型参数进行同时建模和求解。Formula (8) is exactly the time point in the method of the present invention
Figure 656632DEST_PATH_IMAGE060
target cell
Figure 402871DEST_PATH_IMAGE113
The virtual observation equation constructed at . For each moment of each pixel, the virtual observation equation of formula (8) can be established. At this time, the observation value vector
Figure 226471DEST_PATH_IMAGE019
Add a 0-value element in , and in the corresponding coefficient matrix
Figure 349148DEST_PATH_IMAGE020
Add a row in , and the elements corresponding to the model parameters in formula (8) in this row are the corresponding elements of the coefficient matrix in formula (8). It can be seen that, similar to the modeling process in space, in the virtual observation equation in time, for a certain moment, the model parameters corresponding to the pixel at other moments will also be included, so it is also necessary to simultaneously carry out the model parameters at all moments Modeling and Solving.

除了上述的三种方程建立方法之外，一般也可以利用其他先验信息建立额外的虚拟观测方程。比如说，在一些研究区域，发生形变的范围往往小于整个SAR影像的覆盖范围，因此就会存在远场区域，认为远场区域是没有发生地表形变的，也就是说：In addition to the above three equation establishment methods, generally other prior information can also be used to establish additional virtual observation equations. For example, in some research areas, the range of deformation is often smaller than the coverage of the entire SAR image, so there will be a far-field area, which is considered to have no surface deformation, that is to say:

Figure 574593DEST_PATH_IMAGE143
（9）
Figure 574593DEST_PATH_IMAGE143
(9)

类似公式（8），上述虚拟观测方程也可以加入最终的模型参数解算函数模型中（即

Figure 757312DEST_PATH_IMAGE144
）。这种先验假设条件往往也是成立的，因为在InSAR数据处理过程中，不可避免的需要在空间上以某个点为参考点进行空间相位解缠，即本身InSAR数据的测量结果是一个相对结果。对于目标研究区域而言，基于先验信息人工选取一些区域，假定其形变为0，这种情况下可进一步提高本发明实施例中联合模型的可靠性。Similar to formula (8), the above virtual observation equation can also be added to the final model parameter solution function model (ie
Figure 757312DEST_PATH_IMAGE144
). This prior assumption is often true, because in the process of InSAR data processing, it is inevitable to use a certain point as a reference point in space to perform spatial phase unwrapping, that is, the measurement result of InSAR data itself is a relative result . For the target research area, some areas are artificially selected based on prior information, and the deformation is assumed to be 0. In this case, the reliability of the joint model in the embodiment of the present invention can be further improved.

综上，基于公式（2）、（4）、（8）、（9）即可得到用于求解模型参数

Figure 802629DEST_PATH_IMAGE051
的观测值向量
Figure 463417DEST_PATH_IMAGE019
和系数矩阵
Figure 543369DEST_PATH_IMAGE020
，其中
Figure 896990DEST_PATH_IMAGE145
满足以下公式In summary, based on the formulas (2), (4), (8), and (9), the parameters used to solve the model can be obtained
Figure 802629DEST_PATH_IMAGE051
The observation vector of
Figure 463417DEST_PATH_IMAGE019
and coefficient matrix
Figure 543369DEST_PATH_IMAGE020
,in
Figure 896990DEST_PATH_IMAGE145
satisfy the following formula

Figure 429602DEST_PATH_IMAGE146
（10）
Figure 429602DEST_PATH_IMAGE146
(10)

式（10）即为建立的InSAR时序相位中空间趋势误差、地形相关误差和形变信号的联合模型，其中，

Figure 894082DEST_PATH_IMAGE019
为总的观测值向量，
Figure 828540DEST_PATH_IMAGE020
为系数矩阵，
Figure 353062DEST_PATH_IMAGE051
为模型参数。Equation (10) is the joint model of spatial trend error, terrain-related error and deformation signal in the established InSAR time series phase, where,
Figure 894082DEST_PATH_IMAGE019
is the total observation vector,
Figure 828540DEST_PATH_IMAGE020
is the coefficient matrix,
Figure 353062DEST_PATH_IMAGE051
is the model parameter.

具体来说，在本发明实施例中，由于公式（10）中的系数矩阵

Figure 372970DEST_PATH_IMAGE020
为大型系数矩阵，因此在步骤3中采用现有的稀疏最小二乘迭代算法对模型参数
Figure 641141DEST_PATH_IMAGE051
进行求解，通过matlab中的lsmr函数实现公式（10）中模型参数的快速高精度解算，得到在每个时刻每一个像元处相应的空间趋势误差和地形相关误差。Specifically, in the embodiment of the present invention, since the coefficient matrix in formula (10)
Figure 372970DEST_PATH_IMAGE020
is a large coefficient matrix, so in step 3, the existing sparse least squares iterative algorithm is used to modify the model parameters
Figure 641141DEST_PATH_IMAGE051
To solve it, the fast and high-precision calculation of the model parameters in formula (10) is realized through the lsmr function in matlab, and the corresponding spatial trend error and terrain-related error at each pixel at each time are obtained.

具体来说，步骤4在原始InSAR时序相位中将空间趋势误差和地形相关误差减去，得到改正后的高精度InSAR时序形变。Specifically, in step 4, the spatial trend error and terrain-related error are subtracted from the original InSAR time series phase to obtain the corrected high-precision InSAR time series deformation.

综上所述，本发明实施例在利用传统方法得到InSAR时序相位的基础上，认为在一定的空间范围内（如1km×1km），湍流大气和轨道误差可以利用与位置相关的一阶多项式进行建模；可利用分层大气与局部地形相关进行线性建模；在一定的时间范围内（如1个月），绝大多数地表形变可认为是与时间相关的平滑过程，可利用与时间相关的三次函数进行拟合；在InSAR影像中，远场区域的形变等于0，进而可将这一先验信息作为约束条件辅助建模。本发明基于这4种先验信息，实现了InSAR时序相位中形变、大气延迟和轨道误差等不同信号的联合建模，在此基础上，利用稀疏最小二乘迭代算法进行求解，即可实现InSAR时序相位中大气延迟和轨道误差的高精度估计与改正，进而显著提高InSAR时序形变的测量精度；突破了传统方法对各类误差单独处理的思路，充分考虑了多源信号的时空特征，极大地丰富了时序InSAR地表形变测量的理论与方法体系。In summary, on the basis of obtaining the InSAR time series phase by using the traditional method, the embodiment of the present invention considers that within a certain space range (such as 1km×1km), the turbulent atmosphere and the orbit error can be calculated by using the position-related first-order polynomial Modeling; layered atmosphere can be used to perform linear modeling related to local terrain; within a certain time range (such as 1 month), most surface deformation can be considered as a time-related smooth process, and time-related In the InSAR image, the deformation of the far field area is equal to 0, and this prior information can be used as a constraint to assist in modeling. Based on these four kinds of prior information, the present invention realizes the joint modeling of different signals such as deformation, atmospheric delay and orbit error in the InSAR time series phase. On this basis, the sparse least squares iterative algorithm is used to solve the InSAR The high-precision estimation and correction of atmospheric delay and orbital error in the time-series phase can significantly improve the measurement accuracy of InSAR time-series deformation; it breaks through the traditional method of separately processing various errors, and fully considers the spatio-temporal characteristics of multi-source signals, greatly improving It enriches the theory and method system of time-series InSAR surface deformation measurement.

下面通过以下模拟实验对本发明实施例作进一步说明：The embodiment of the present invention will be further described below through the following simulation experiments:

通过SBAS方法模拟得到22个时刻的InSAR时序形变，形变场的空间尺寸为600×600像素，每个像素的空间分辨率为100m×100m。其中，形变的空间特征为漏斗形，每个点在时间上的变化趋势为对数形，即

Figure 695684DEST_PATH_IMAGE147
，形变
Figure 125529DEST_PATH_IMAGE148
的单位为米，时间
Figure 632733DEST_PATH_IMAGE149
表示相对于参考时刻的时间差，单位为天；利用模拟了每一个时刻整个形变场范围内的轨道误差，其中轨道误差的最大量级为10弧度（对于c波段而言，约为4.4厘米），轨道误差在空间上展示三维趋势方向在每个时刻是随机的；每个时刻的湍流大气误差可利用分型维度为2.2的分形函数进行模拟，最大量级为10弧度；对于地形相关的分层乏汽而言，考虑到分层大气与地形之间的比例系数在空间不同区域可能不一致，模拟实验中同样利用分形维度为2.2的分形函数模拟空间上不同区域处分层大气相位与地形之间的比例系数，然后通过模拟得到的比例系数与真实DEM数据（Digital Elevation Model，DEM）全球数字高程数据相乘得到模拟数据中的分层大气分量；将模拟的InSAR时序形变、轨道误差、湍流大气和分层大气相加，即可得到模拟试验所用到的InSAR时序相位。The InSAR time-series deformation at 22 moments is simulated by the SBAS method. The spatial size of the deformation field is 600×600 pixels, and the spatial resolution of each pixel is 100m×100m. Among them, the spatial characteristics of the deformation are funnel-shaped, and the change trend of each point in time is logarithmic, namely
Figure 695684DEST_PATH_IMAGE147
,deformation
Figure 125529DEST_PATH_IMAGE148
The unit is meter, time
Figure 632733DEST_PATH_IMAGE149
Indicates the time difference relative to the reference moment, and the unit is day; by simulating the orbit error within the entire deformation field range at each moment, the maximum magnitude of the orbit error is 10 radians (for the c-band, it is about 4.4 cm), The orbital error shows that the three-dimensional trend direction in space is random at each moment; the turbulent atmospheric error at each moment can be simulated by using a fractal function with a fractal dimension of 2.2, and the maximum magnitude is 10 radians; for terrain-related stratification For exhaust steam, considering that the proportional coefficient between the stratified atmosphere and the terrain may be inconsistent in different regions in space, the fractal function with a fractal dimension of 2.2 is also used in the simulation experiment to simulate the relationship between the stratified atmosphere phase and the terrain in different regions in space. , and then multiplied by the simulated proportional coefficient and the real DEM data (Digital Elevation Model, DEM) global digital elevation data to obtain the stratified atmospheric component in the simulated data; the simulated InSAR time series deformation, orbit error, turbulent atmosphere The InSAR timing phase used in the simulation experiment can be obtained by adding it to the layered atmosphere.

为了对比本发明实施例的优势，将模拟的InSAR时序相位作为观测值，分别利用本发明实施例所提到的方法与传统方法进行InSAR时序相位中的大气延迟误差和轨道误差改正。其中，传统方法是指在整个影像中，将形变区域掩膜，利用非形变区与进行趋势相位和地形相关相位联合建模和求解，然后利用求解得到的模型参数进行形变区域的趋势相位和地形相关相位正演，即可实现整个影像中相关误差的改正，如图2所示展示了漏斗中心处不同方法获取的InSAR时序形变结果可看出本发明实施例所提到的方法比传统方法得到的InSAR时序形变结果更为精确。In order to compare the advantages of the embodiments of the present invention, the simulated InSAR time series phase is used as the observed value, and the method mentioned in the embodiment of the present invention and the traditional method are used to correct the atmospheric delay error and orbit error in the InSAR time series phase. Among them, the traditional method refers to masking the deformed area in the entire image, using the non-deformed area to jointly model and solve the trend phase and terrain-related phase, and then use the model parameters obtained from the solution to calculate the trend phase and terrain of the deformed area. Correlation phase forward modeling can realize the correction of correlation errors in the entire image. As shown in Figure 2, the InSAR time series deformation results obtained by different methods at the center of the funnel are shown. It can be seen that the method mentioned in the embodiment of the present invention is better than the traditional method. The deformation results of InSAR time series are more accurate.

本发明实施例还提供了一种计算机可读存储介质，用于存储计算机程序，通过执行计算机程序，用于实现上述的时序InSAR误差的改正方法。An embodiment of the present invention also provides a computer-readable storage medium for storing a computer program, and implementing the above-mentioned method for correcting timing InSAR errors by executing the computer program.

本发明实施例还提供了一种时序InSAR误差的改正设备，用于实现上述的时序InSAR误差的改正方法，包括：The embodiment of the present invention also provides a timing InSAR error correction device, which is used to implement the above timing InSAR error correction method, including:

存储器和处理器；memory and processor;

存储器用于储存计算机程序；Memory is used to store computer programs;

处理器用于执行存储器存储的计算机程序。The processor is used to execute the computer program stored in the memory.

以上所述是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明所述原理的前提下，还可以作出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。The above description is a preferred embodiment of the present invention, it should be pointed out that for those of ordinary skill in the art, without departing from the principle of the present invention, some improvements and modifications can also be made, and these improvements and modifications can also be made. It should be regarded as the protection scope of the present invention.

Claims (10)

1.一种时序InSAR误差的改正方法，其特征在于，包括：1. A method for correcting a time series InSAR error, characterized in that it comprises: 步骤1，获取InSAR干涉图的InSAR时序相位；Step 1, obtain the InSAR timing phase of the InSAR interferogram; 步骤2，基于多源信号的时空特征，在所述InSAR时序相位中对空间趋势误差、地形相关误差和形变信号建立联合模型，并对所述联合模型中的模型参数进行求解，得到空间趋势误差和地形相关误差；Step 2, based on the spatio-temporal characteristics of multi-source signals, establish a joint model for the spatial trend error, terrain-related error and deformation signal in the InSAR time series phase, and solve the model parameters in the joint model to obtain the spatial trend error and terrain-related errors; 步骤3，去除所述InSAR时序相位中的所述空间趋势误差和所述地形相关误差，得到改正后的InSAR时序形变。Step 3, removing the spatial trend error and the terrain-related error in the InSAR time series phase to obtain the corrected InSAR time series deformation. 2.根据权利要求1所述的时序InSAR误差的改正方法，其特征在于，步骤2还包括：2. the correction method of time series InSAR error according to claim 1, is characterized in that, step 2 also comprises: 利用一阶多项式对所述空间趋势误差进行建模，得到多项式拟合系数方程；Using a first-order polynomial to model the spatial trend error, a polynomial fitting coefficient equation is obtained; 利用线性模型对所述地形相关误差进行表征，得到线性方程系数方程；Using a linear model to characterize the terrain-related error to obtain a linear equation coefficient equation; 利用三阶多项式函数对形变信号构建虚拟观测方程；Construct a virtual observation equation for the deformation signal by using a third-order polynomial function; 联立所述多项式拟合系数方程、线性方程系数方程和所述虚拟观测方程，得到联合模型。A joint model is obtained by combining the polynomial fitting coefficient equation, the linear equation coefficient equation and the virtual observation equation. 3.根据权利要求2所述的时序InSAR误差的改正方法，其特征在于，所述联合模型中的模型参数包括所有像元的时序形变、所有像元处空间趋势误差的多项式拟合系数以及所有像元处地形相关误差的线性方程系数。3. The method for correcting time series InSAR errors according to claim 2, wherein the model parameters in the joint model include the time series deformation of all pixels, the polynomial fitting coefficients of spatial trend errors at all pixel places, and all Coefficients of the linear equation for the terrain-dependent error at the cell. 4.根据权利要求3所述的时序InSAR误差的改正方法，其特征在于，所述步骤2还包括：4. the correction method of time series InSAR error according to claim 3, is characterized in that, described step 2 also comprises: 从所有像元中任意选择一个像元作为目标像元；Randomly select a cell from all the cells as the target cell; 在时刻
Figure 796490DEST_PATH_IMAGE001
，以所述目标像元
Figure 120155DEST_PATH_IMAGE002
为中心取大小为
Figure 368734DEST_PATH_IMAGE003
的窗口，建立所述多项式拟合系数方程为：at the moment
Figure 796490DEST_PATH_IMAGE001
, with the target cell
Figure 120155DEST_PATH_IMAGE002
Take the size of the center as
Figure 368734DEST_PATH_IMAGE003
window, establish the polynomial fitting coefficient equation as:
Figure 704556DEST_PATH_IMAGE004
Figure 704556DEST_PATH_IMAGE004
 其中，
Figure 987770DEST_PATH_IMAGE005
为窗口内每个像元处的InSAR相位观测值，
Figure 747915DEST_PATH_IMAGE006
为每一个点处的形变相位，
Figure 483790DEST_PATH_IMAGE007
为
Figure 151532DEST_PATH_IMAGE008
区间内的整数，
Figure 554831DEST_PATH_IMAGE009
为地形相关误差相位，
Figure 17037DEST_PATH_IMAGE010
、
Figure 505787DEST_PATH_IMAGE011
、
Figure 711640DEST_PATH_IMAGE012
为多项式拟合系数；in,
Figure 987770DEST_PATH_IMAGE005
is the InSAR phase observation value at each pixel in the window,
Figure 747915DEST_PATH_IMAGE006
is the deformation phase at each point,
Figure 483790DEST_PATH_IMAGE007
for
Figure 151532DEST_PATH_IMAGE008
integers in the interval,
Figure 554831DEST_PATH_IMAGE009
is the terrain-related error phase,
Figure 17037DEST_PATH_IMAGE010
,
Figure 505787DEST_PATH_IMAGE011
,
Figure 711640DEST_PATH_IMAGE012
is the polynomial fitting coefficient; 将上式中的观测值向量
Figure 703867DEST_PATH_IMAGE013
纳入总的观测值向量
Figure 868132DEST_PATH_IMAGE014
中，则相应的系数矩阵
Figure 578599DEST_PATH_IMAGE015
会增加
Figure 588143DEST_PATH_IMAGE016
行，每一行与上式中观测值向量
Figure 700456DEST_PATH_IMAGE013
的元素相对应，每一行中大部分元素为0，只有与模型参数向量
Figure 770043DEST_PATH_IMAGE017
中元素相对应的列不为0，其数值为上式中与系数矩阵
Figure 967806DEST_PATH_IMAGE018
相应行对应的元素。The observation vector in the above formula
Figure 703867DEST_PATH_IMAGE013
Include total observations vector
Figure 868132DEST_PATH_IMAGE014
, then the corresponding coefficient matrix
Figure 578599DEST_PATH_IMAGE015
will increase
Figure 588143DEST_PATH_IMAGE016
row, each row is related to the observation vector in the above formula
Figure 700456DEST_PATH_IMAGE013
Corresponding to the elements, most of the elements in each row are 0, only the model parameter vector
Figure 770043DEST_PATH_IMAGE017
The column corresponding to the element in is not 0, and its value is the coefficient matrix in the above formula
Figure 967806DEST_PATH_IMAGE018
The element corresponding to the corresponding row. 5.根据权利要求4所述的时序InSAR误差的改正方法，其特征在于，所述步骤2还包括：5. the correcting method of timing InSAR error according to claim 4, is characterized in that, described step 2 also comprises: 从所有像元中任意选择一个像元作为目标像元；Randomly select a cell from all the cells as the target cell; 在时刻
Figure 781041DEST_PATH_IMAGE001
，以所述目标像元
Figure 744931DEST_PATH_IMAGE002
为中心取大小为
Figure 250998DEST_PATH_IMAGE019
的窗口，建立所述线性方程系数方程：at the moment
Figure 781041DEST_PATH_IMAGE001
, with the target cell
Figure 744931DEST_PATH_IMAGE002
Take the size of the center as
Figure 250998DEST_PATH_IMAGE019
window, build the coefficient equation of the linear equation:
Figure 670478DEST_PATH_IMAGE020
Figure 670478DEST_PATH_IMAGE020
 其中，
Figure 21825DEST_PATH_IMAGE005
为窗口内每个像元处的InSAR相位观测值，
Figure 108730DEST_PATH_IMAGE006
为形变相位，
Figure 520120DEST_PATH_IMAGE021
为趋势相位，
Figure 692475DEST_PATH_IMAGE022
为窗口区间内的整数，
Figure 581934DEST_PATH_IMAGE023
为像元
Figure 523345DEST_PATH_IMAGE002
与像元
Figure 777740DEST_PATH_IMAGE024
之间的高程之差，
Figure 171812DEST_PATH_IMAGE025
、
Figure 864962DEST_PATH_IMAGE026
为线性方程系数；in,
Figure 21825DEST_PATH_IMAGE005
is the InSAR phase observation value at each pixel in the window,
Figure 108730DEST_PATH_IMAGE006
is the deformation phase,
Figure 520120DEST_PATH_IMAGE021
is the trend phase,
Figure 692475DEST_PATH_IMAGE022
is an integer in the window interval,
Figure 581934DEST_PATH_IMAGE023
for the pixel
Figure 523345DEST_PATH_IMAGE002
with the pixel
Figure 777740DEST_PATH_IMAGE024
The difference in elevation between,
Figure 171812DEST_PATH_IMAGE025
,
Figure 864962DEST_PATH_IMAGE026
is the linear equation coefficient; 将上式中的观测值向量
Figure 129721DEST_PATH_IMAGE013
纳入总的观测值向量
Figure 882913DEST_PATH_IMAGE014
中，则相应的系数矩阵
Figure 778930DEST_PATH_IMAGE015
会增加
Figure 275770DEST_PATH_IMAGE016
行，每一行与上式中观测值向量
Figure 395036DEST_PATH_IMAGE013
的元素相对应，每一行中大部分元素为0，只有与模型参数向量
Figure 53551DEST_PATH_IMAGE017
中元素相对应的列不为0，其数值为上式中系数矩阵
Figure 422215DEST_PATH_IMAGE018
相应行对应的元素。The observation vector in the above formula
Figure 129721DEST_PATH_IMAGE013
Include total observations vector
Figure 882913DEST_PATH_IMAGE014
, then the corresponding coefficient matrix
Figure 778930DEST_PATH_IMAGE015
will increase
Figure 275770DEST_PATH_IMAGE016
row, each row is related to the observation vector in the above formula
Figure 395036DEST_PATH_IMAGE013
Corresponding to the elements, most of the elements in each row are 0, only the model parameter vector
Figure 53551DEST_PATH_IMAGE017
The column corresponding to the element in is not 0, and its value is the coefficient matrix in the above formula
Figure 422215DEST_PATH_IMAGE018
The element corresponding to the corresponding row. 6.根据权利要求5所述的时序InSAR误差的改正方法，其特征在于，所述步骤2还包括：6. the correcting method of timing InSAR error according to claim 5, is characterized in that, described step 2 also comprises: 从所有像元中任意选择一个像元作为目标像元；Randomly select a cell from all the cells as the target cell; 在时刻
Figure 926009DEST_PATH_IMAGE001
，以所述目标像元
Figure 962098DEST_PATH_IMAGE002
为中心取大小为的时间窗口，构建虚拟观测方程为：at the moment
Figure 926009DEST_PATH_IMAGE001
, with the target cell
Figure 962098DEST_PATH_IMAGE002
Taking a time window of size as the center, the virtual observation equation is constructed as:
Figure 791514DEST_PATH_IMAGE027
Figure 791514DEST_PATH_IMAGE027
 其中，
Figure 913053DEST_PATH_IMAGE028
，
Figure 751696DEST_PATH_IMAGE029
为中间模型参数变量，
Figure 111134DEST_PATH_IMAGE030
与
Figure 642609DEST_PATH_IMAGE031
之差等于0，
Figure 454707DEST_PATH_IMAGE032
为时序形变，
Figure 362620DEST_PATH_IMAGE033
为时序形变
Figure 310985DEST_PATH_IMAGE034
的三次多项式拟合值，
Figure 744853DEST_PATH_IMAGE035
、
Figure 309826DEST_PATH_IMAGE036
、
Figure 755851DEST_PATH_IMAGE037
分别为第
Figure 89880DEST_PATH_IMAGE038
景SAR影像与第
Figure 166421DEST_PATH_IMAGE040
景SAR影像时刻之间的时间差、第
Figure 218690DEST_PATH_IMAGE041
景SAR影像与第
Figure 937248DEST_PATH_IMAGE040
景SAR影像时刻之间的时间差、第
Figure 656942DEST_PATH_IMAGE042
景SAR影像与第
Figure 435542DEST_PATH_IMAGE043
景SAR影像时刻之间的时间差，
Figure 240687DEST_PATH_IMAGE044
。in,
Figure 913053DEST_PATH_IMAGE028
,
Figure 751696DEST_PATH_IMAGE029
is the intermediate model parameter variable,
Figure 111134DEST_PATH_IMAGE030
and
Figure 642609DEST_PATH_IMAGE031
The difference is equal to 0,
Figure 454707DEST_PATH_IMAGE032
is the time series deformation,
Figure 362620DEST_PATH_IMAGE033
time series deformation
Figure 310985DEST_PATH_IMAGE034
The cubic polynomial fit value of ,
Figure 744853DEST_PATH_IMAGE035
,
Figure 309826DEST_PATH_IMAGE036
,
Figure 755851DEST_PATH_IMAGE037
respectively
Figure 89880DEST_PATH_IMAGE038
SAR image and the first
Figure 166421DEST_PATH_IMAGE040
The time difference between scene SAR image moments, the first
Figure 218690DEST_PATH_IMAGE041
SAR image and the first
Figure 937248DEST_PATH_IMAGE040
The time difference between scene SAR image moments, the first
Figure 656942DEST_PATH_IMAGE042
SAR image and the first
Figure 435542DEST_PATH_IMAGE043
The time difference between the scene SAR image moments,
Figure 240687DEST_PATH_IMAGE044
. 7.根据权利要求6所述的时序InSAR误差的改正方法，其特征在于，所述步骤2还包括：7. the correction method of time series InSAR error according to claim 6, is characterized in that, described step 2 also comprises: 利用先验信息建立额外的虚拟观测方程为：Using the prior information to establish an additional virtual observation equation is:
Figure 762935DEST_PATH_IMAGE045
Figure 762935DEST_PATH_IMAGE045
 将该虚拟观测方程加入所述总的观测值向量
Figure 71557DEST_PATH_IMAGE014
中，并在对应的系数矩阵
Figure 286638DEST_PATH_IMAGE015
中加入一行。Add this dummy observation equation to the total observation vector
Figure 71557DEST_PATH_IMAGE014
, and in the corresponding coefficient matrix
Figure 286638DEST_PATH_IMAGE015
Add a line to the . 8.根据权利要求7所述的时序InSAR误差的改正方法，其特征在于，所述联立所述多项式拟合系数方程、线性方程系数方程和所述虚拟观测方程，得到联合模型，包括：8. The correcting method of time series InSAR error according to claim 7, is characterized in that, described simultaneous polynomial fitting coefficient equation, linear equation coefficient equation and described virtual observation equation, obtain joint model, comprise: 基于所述多项式拟合系数方程、线性方程系数方程、所述虚拟观测方程以及额外的虚拟观测方程，建立所述InSAR时序相位中所述空间趋势误差、所述地形相关误差和所述形变信号的联合模型为：Based on the polynomial fitting coefficient equation, the linear equation coefficient equation, the virtual observation equation and the additional virtual observation equation, establish the spatial trend error, the terrain-related error and the deformation signal in the InSAR time series phase The joint model is:
Figure 313500DEST_PATH_IMAGE046
Figure 313500DEST_PATH_IMAGE046
 其中，
Figure 639439DEST_PATH_IMAGE014
为总的观测值向量，
Figure 802567DEST_PATH_IMAGE015
为系数矩阵，
Figure 188549DEST_PATH_IMAGE047
为模型参数。in,
Figure 639439DEST_PATH_IMAGE014
is the total observation vector,
Figure 802567DEST_PATH_IMAGE015
is the coefficient matrix,
Figure 188549DEST_PATH_IMAGE047
is the model parameter. 9.一种计算机可读存储介质，用于存储计算机程序，其特征在于，通过执行所述计算机程序，用于实现上述权利要求1-8中任意一项所述的时序InSAR误差的改正方法。9. A computer-readable storage medium for storing a computer program, characterized in that, by executing the computer program, it is used to implement the method for correcting the timing InSAR error described in any one of claims 1-8. 10.一种时序InSAR误差的改正设备，用于实现上述权利要求1-8任意一项所述的时序InSAR误差的改正方法，其特征在于，包括：10. A timing InSAR error correction device, used to implement the timing InSAR error correction method described in any one of claims 1-8, characterized in that it comprises: 存储器和处理器；memory and processor; 所述存储器用于储存计算机程序；The memory is used to store computer programs; 所述处理器用于执行所述存储器存储的计算机程序。The processor is configured to execute a computer program stored in the memory.
CN202211568937.XA 2022-12-08 2022-12-08 Correction method of time sequence InSAR error and related equipment Active CN115629384B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211568937.XA CN115629384B (en) 2022-12-08 2022-12-08 Correction method of time sequence InSAR error and related equipment

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211568937.XA CN115629384B (en) 2022-12-08 2022-12-08 Correction method of time sequence InSAR error and related equipment

Publications (2)

Publication Number Publication Date
CN115629384A true CN115629384A (en) 2023-01-20
CN115629384B CN115629384B (en) 2023-04-11

Family

ID=84910886

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211568937.XA Active CN115629384B (en) 2022-12-08 2022-12-08 Correction method of time sequence InSAR error and related equipment

Country Status (1)

Country Link
CN (1) CN115629384B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068511A (en) * 2023-03-09 2023-05-05 成都理工大学 Deep learning-based InSAR large-scale system error correction method
CN116299245A (en) * 2023-05-11 2023-06-23 中山大学 An adaptive mosaic correction method for time-series InSAR deformation rate results
CN116609781A (en) * 2023-05-25 2023-08-18 北京理工大学 Beidou InSAR DEM error compensation method combining multiple star data

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103675790A (en) * 2013-12-23 2014-03-26 中国国土资源航空物探遥感中心 Method for improving earth surface shape change monitoring precision of InSAR (Interferometric Synthetic Aperture Radar) technology based on high-precision DEM (Digital Elevation Model)
CN112051572A (en) * 2020-09-14 2020-12-08 广东省核工业地质局测绘院 Method for monitoring three-dimensional surface deformation by fusing multi-source SAR data
US20210011149A1 (en) * 2019-05-21 2021-01-14 Central South University InSAR and GNSS weighting method for three-dimensional surface deformation estimation
CN113064170A (en) * 2021-03-29 2021-07-02 长安大学 Expansive soil area surface deformation monitoring method based on time sequence InSAR technology
CN114415178A (en) * 2022-01-13 2022-04-29 姚鑫 A Fast InSAR Processing Method for Deformation Geological Hazard Identification—GHR-InSAR

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103675790A (en) * 2013-12-23 2014-03-26 中国国土资源航空物探遥感中心 Method for improving earth surface shape change monitoring precision of InSAR (Interferometric Synthetic Aperture Radar) technology based on high-precision DEM (Digital Elevation Model)
US20210011149A1 (en) * 2019-05-21 2021-01-14 Central South University InSAR and GNSS weighting method for three-dimensional surface deformation estimation
CN112051572A (en) * 2020-09-14 2020-12-08 广东省核工业地质局测绘院 Method for monitoring three-dimensional surface deformation by fusing multi-source SAR data
CN113064170A (en) * 2021-03-29 2021-07-02 长安大学 Expansive soil area surface deformation monitoring method based on time sequence InSAR technology
CN114415178A (en) * 2022-01-13 2022-04-29 姚鑫 A Fast InSAR Processing Method for Deformation Geological Hazard Identification—GHR-InSAR

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
LUO, XIAN GANG 等: "Dynamic analysis of urban ground subsidence in Beijing based on the permanent scattering InSAR technology" *
韦建超 等: "附加系统参数的多时相InSAR时空建模和形变估计" *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116068511A (en) * 2023-03-09 2023-05-05 成都理工大学 Deep learning-based InSAR large-scale system error correction method
CN116299245A (en) * 2023-05-11 2023-06-23 中山大学 An adaptive mosaic correction method for time-series InSAR deformation rate results
CN116299245B (en) * 2023-05-11 2023-07-28 中山大学 An adaptive mosaic correction method for time-series InSAR deformation rate results
CN116609781A (en) * 2023-05-25 2023-08-18 北京理工大学 Beidou InSAR DEM error compensation method combining multiple star data
CN116609781B (en) * 2023-05-25 2025-11-04 北京理工大学 A method for error compensation of BeiDou InSAR DEM based on multi-satellite data

Also Published As

Publication number Publication date
CN115629384B (en) 2023-04-11

Similar Documents

Publication Publication Date Title
CN115629384B (en) Correction method of time sequence InSAR error and related equipment
US11226431B2 (en) Method and device for filling invalid regions of terrain elevation model data
CN111650579B (en) InSAR mining area three-dimensional deformation estimation method and device for rock migration parameter adaptive acquisition and medium
CN117724098B (en) A time-series InSAR tropospheric delay correction method for alleviating atmospheric seasonal oscillations
CN111724465B (en) Satellite image adjustment method and device based on plane constraint optimization virtual control point
CN106780321B (en) CBERS-02 satellite HR sensor image overall tight orientation and correction splicing method
CN103310487B (en) A kind of universal imaging geometric model based on time variable generates method
CN112991535B (en) Three-dimensional space situation representation method and device of height information enhanced ink cartoo map
WO2021012592A1 (en) Coseismic and postseismic temporal-spatial slip distribution joint inversion method based on multisource sar data with additional logarithmic constraints
CN109636912A (en) Tetrahedron subdivision finite element interpolation method applied to three-dimensional sonar image reconstruction
CN106885585A (en) A kind of satellite borne photography measuring system integration calibration method based on bundle adjustment
CN118191839A (en) Surface three-dimensional deformation inversion method and system
CN109886910B (en) DEM (digital elevation model) correction method and device for external digital elevation model
CN105023281B (en) Asterism based on point spread function wavefront modification is as centroid computing method
CN115469308B (en) Multi-track InSAR interseismic deformation rate field splicing method, device, equipment and medium
Pi et al. Large-scale planar block adjustment of GaoFen1 WFV images covering most of mainland China
CN110310370B (en) Method for point-plane fusion of GPS (Global positioning System) and SRTM (short Range TM)
CN112835043A (en) A Monitoring Method for Two-dimensional Deformation in Arbitrary Direction
CN114969227B (en) Method and device for constructing regional scale atmospheric knowledge base based on latitude and longitude grid
CN112464433B (en) A Recursive Least Squares Solution for FPGA Hardware-Oriented Parameter Optimization Algorithm for RFM Models
CN115951353A (en) Multi-track InSAR extraction high-precision three-dimensional deformation method
CN114964169B (en) Remote sensing image adjustment method for image space object space cooperative correction
CN111879440A (en) Method and device for calculating surface temperature
CN120337779B (en) A method for oilfield modeling integrating InSAR 3D deformation and geophysical model
CN118097147B (en) Image shielding detection method and device

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant