CN115453534B - 一种顾及解缠误差的序贯InSAR时序形变解算方法 - Google Patents

一种顾及解缠误差的序贯InSAR时序形变解算方法 Download PDF

Info

Publication number
CN115453534B
CN115453534B CN202211137798.5A CN202211137798A CN115453534B CN 115453534 B CN115453534 B CN 115453534B CN 202211137798 A CN202211137798 A CN 202211137798A CN 115453534 B CN115453534 B CN 115453534B
Authority
CN
China
Prior art keywords
historical
deformation
sequential
phase
matrix
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.)
Active
Application number
CN202211137798.5A
Other languages
English (en)
Other versions
CN115453534A (zh
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.)
Sun Yat Sen University
Original Assignee
Sun Yat Sen 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 Sun Yat Sen University filed Critical Sun Yat Sen University
Priority to CN202211137798.5A priority Critical patent/CN115453534B/zh
Publication of CN115453534A publication Critical patent/CN115453534A/zh
Application granted granted Critical
Publication of CN115453534B publication Critical patent/CN115453534B/zh
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
    • 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/904—SAR modes
    • 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
    • G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • G01S7/40—Means for monitoring or calibrating
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Software Systems (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Electromagnetism (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computing Systems (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Radar Systems Or Details Thereof (AREA)
  • Image Processing (AREA)

Abstract

本发明提供一种顾及解缠误差的序贯InSAR时序形变解算方法，如下：对历史N+1幅SAR影像序列，预处理得到解缠差分干涉相位图，建立性方程组；估计出历史形变相位，并引入M估计量，采用选权迭代法进行历史形变相位的求解，迭代得到协因数阵；获取新增N1景SAR影像，并与历史影像通过预处理，得到解算对象；对新增M1个观测值建立观测方程并写成误差方程的形式；根据历史形变相位建立虚拟观测方程，根据贝叶斯平差原理，构造目标函数；对目标函数求解得到序贯平差的参数估值；引入M估计量，采用选权迭代法，更新得到序贯平差的参数估值的协因数阵；通过对获得更新的序贯平差拟合获取形变速率参数。本发明方法避免参数估值受异常误差影响而歪曲，保证了时序InSAR序贯形变估计结果的可靠性。

Description

一种顾及解缠误差的序贯InSAR时序形变解算方法
技术领域
本发明涉及SAR数据处理技术领域，更具体地，涉及一种顾及解缠误差的序贯InSAR时序形变解算方法。
背景技术
随着当代SAR卫星重返频率逐步提高，海量观测数据不仅为更加客观地了解形变的发展动态和规律提供了契机，同时也给快速、高效的形变解算带来了新的挑战。对于标准时序InSAR数据处理，每新增一幅观测数据需要对历史数据重新解算，计算效率低，滞后性严重。这使得在实现大范围InSAR准实时监测要求和处理效率上变得难以兼顾。目前，时序InSAR技术已进入序贯数据处理阶段，其核心思想是以历史数据的形变解算结果为基础，结合新增影像观测数据，利用贝叶斯序贯平差对历史形变进行增量解算，进而达到整体解算的作用，实现高效运算。
相位解缠是InSAR技术中有待攻克的难题，在时序InSAR相位解缠过程中，不可避免的会存在解缠误差，即以一种粗差的形式存在于观测值当中。而现有时序InSAR序贯平差方法没有考虑解缠误差对形变解算的影响，不具备抵抗粗差干扰的能力。当历史观测数据存在解缠误差时，历史形变解算结果不可靠，同时新增观测数据中也可能存在解缠误差，导致序贯平差解算出的参数估值严重歪曲。因此，在进行序贯估计中，必须同时考虑历史和新增观测数据受异常干扰时的质量控制及其数据处理问题。
发明内容
本发明为克服上述现有时序InSAR序贯平差方法没有考虑解缠误差会严重歪曲序贯解算结果的问题，提出了一种顾及解缠误差的序贯InSAR时序形变解算方法，能避免参数估值受异常误差影响而歪曲，保证时序InSAR序贯形变估计结果的可靠性。
为解决上述技术问题，本发明的技术方案如下：
一种顾及解缠误差的序贯InSAR时序形变解算方法，所述的方法包括步骤如下：
S1：对历史N+1幅SAR影像序列，通过预处理得到其解缠差分干涉相位图上单个高相干点，根据单个高相干点的M个历史解缠相位观测值建立关于设计矩阵、历史形变相位、历史形变相位的协因数阵的线性方程组；
S2：根据最小二乘原理估计出历史形变相位，并引入M估计量来抵制历史数据集中解缠误差的影响，采用选权迭代法进行历史形变相位的求解，迭代得到历史形变相位、历史形变相位的协因数阵；
S3：获取研究区域新增N1景SAR影像，将新增影像与邻近时刻的历史影像通过预处理，得到以历史解缠差分干涉相位图中选取的高相干点目标作为解算对象；
S4：对新增M1个观测值建立包含历史形变相位及新增形变相位的观测方程；将观测方程写成误差方程的形式；
S5：根据历史形变相位中的先验信息及历史形变相位的协因数阵，建立虚拟观测方程，结合步骤S5所述的误差方程，根据贝叶斯平差原理，构造目标函数；对目标函数求解得到序贯平差的参数估值；
S6：引入M估计量来抵制新增观测数据中解缠误差对参数估值的影响，采用选权迭代法进行新增形变相位的求解，更新得到序贯平差的参数估值的协因数阵；
S7：根据步骤S1～S6获得更新的序贯平差，然后通过对更新后的序贯平差拟合获取形变速率参数。
优选地，所述的预处理是按时空基线集方法处理，得到M幅解缠差分干涉相位图，经过大气误差、轨道误差及DEM误差校正之后选取高相干点作为解算对象。
进一步地，所述的线性方程组；
Figure BDA0003852904880000021
其中，G为M×N的设计矩阵，与历史解缠差分干涉相位图组成的时空基线网络有关，其第l行代表的第l幅解缠差分干涉相位图由第s，t幅SAR影像干涉生成，0≤s＜t≤N，当s＞0时，对应矩阵G的第l行的第s个元素为-1；当s＝0时，矩阵G的第t个元素为1，其余元素均为0；
Figure BDA0003852904880000031
表示第n个时刻的形变相位，
Figure BDA0003852904880000032
表示历史解缠相位序列，n＝1,…,N，m＝1,…,M。
再进一步地，若M＞N，并且设计矩阵G是列满秩的，根据最小二乘原理估计出历史形变相位
Figure BDA0003852904880000033
P0为历史观测值L0的权矩阵；残差向量V0＝[v1,v2,…,vM]T通过
Figure BDA0003852904880000034
计算得到，其中下标0表示变量为初始化计算过程中的变量。
再进一步地，步骤S2，所述的采用选权迭代法进行历史形变相位的求解，其迭代计算公式如下：
Figure BDA0003852904880000035
其中，
Figure BDA0003852904880000036
分别为第k次迭代计算的历史形变相位，历史残差向量，历史形变相位的协因数阵，历史观测值L0的等价权矩阵；w(k-1)为第k-1次迭代计算的权因子；
当前后两次迭代计算得到的历史形变相位估值之差
Figure BDA0003852904880000037
满足迭代收敛条件时，即获得历史形变相位的抗差解。
再进一步地，所述的观测方程：
Figure BDA0003852904880000038
其中，G0为与历史形变相位序列有关的系数阵，G1为与新增形变相位序列有关的系数阵，二者的构造与新增解缠差分干涉相位图的网形有关；
假设
Figure BDA0003852904880000039
由tN+1时刻获取的新增影像和tN+1时刻获取的历史影像差分干涉得到，则G0(1,N-1)＝-1，G1(1,1)＝1，其它元素都为0；
Figure BDA0003852904880000041
为新增的观测向量，N1表示新增影像的数量；M1表示新增影像与邻近历史影像生成的新增解缠差分干涉相位图的数量。
再进一步地，所述的误差方程表达式如下：
Figure BDA0003852904880000042
其中，V1为新增观测的残差向量，P1为新增观测的权矩阵；
Figure BDA0003852904880000043
表示序贯平差的参数估值。
再进一步地，所述的虚拟观测方程：
Figure BDA0003852904880000044
Figure BDA0003852904880000045
其中，
Figure BDA00038529048800000413
为历史形变相位
Figure BDA00038529048800000414
的第二次估计，
Figure BDA00038529048800000415
为其对应的残差，
Figure BDA00038529048800000416
为其对应的权矩阵；结合所述的误差方程，根据贝叶斯平差原理，构造如下目标函数：
Figure BDA0003852904880000046
求解可得序贯平差的参数估值为：
Figure BDA0003852904880000047
再进一步地，步骤S6，通过迭代求解以确定新增观测L1的等价权矩阵
Figure BDA0003852904880000048
最后，根据矩阵求逆引理，得到序贯解的递推形式：
Figure BDA0003852904880000049
序贯平差的参数估值
Figure BDA00038529048800000417
的协因数阵为：
Figure BDA00038529048800000410
其中，
Figure BDA00038529048800000411
为增益矩阵，
Figure BDA00038529048800000412
为新增观测对应的更新的协因数阵。
再进一步地，建立第i期新增影像获取时刻的序贯平差更新公式：
Figure BDA0003852904880000051
其中，
Figure BDA0003852904880000052
分别为第i期新增影像获取之前的形变相位序列的第i次更新结果以及第i期影像获取时刻的形变相位估计结果；Jx(i)表示第i期序贯更新公式中的增益矩阵；Li表示第i期的新增观测值；G0(i)表示第i期时空基线网络中与历史序贯平差有关的设计矩阵；G1(i)表示第i期时空基线网络中与新增序贯平差有关的设计矩阵；
Figure BDA0003852904880000053
表示第i期新增观测值对应的更新的协因数阵的逆。
与现有技术相比，本发明技术方案的有益效果是：
本发明在求解历史形变相位序列时，引入M估计量采用选权迭代法进行形变解算，以抵制历史观测值中可能存在的解缠误差的影响，保证在序贯平差序贯更新时先验参数的可靠性。
本发明在利用序贯平差更新序贯平差时，进一步根据M估计量，通过迭代计算获取新增观测等价权，以抵制新增观测值中解缠误差的影响，保证序贯平差更新结果的可靠性；此外，通过对更新的序贯平差拟合获取形变速率，可减少直接利用序贯平差更新形变速率参数时形变模型偏差造成的估计偏差。
本发明同时考虑了在历史和新增解缠差分干涉相位图存在解缠误差时的序贯解算方法。当历史观测和新增观测值不含解缠误差时，该方法接近于最优无偏估计；当历史观测和新增观测值中含有解缠误差时，该方法能同时抵制参数先验异常和观测异常的影响，避免参数估值受异常误差影响而歪曲，极大程度的保证了时序InSAR序贯形变估计结果的可靠性。
附图说明
图1为顾及解缠误差的InSAR时序形变解算方法的更新流程图。
图2为解缠差分干涉相位图时空基线图；三角形代表历史SAR影像，实线代表历史解缠差分干涉相位图；圆圈代表新增SAR影像，虚线代表新增解缠差分干涉相位图。
图3为标准序贯估计与本发明所述的方法的序贯平差结果对比；(a)、(b)线性形变模型；(c)、(d)指数形变模型；(e)、(f)周期形变模型；左栏(a)(c)(e)为不含解缠误差的估计结果；右栏(b)(d)(f)为添加了15％的历史解缠误差以及10％的新增解缠误差的估计结果
图4为蒙特卡洛模拟结果；圆圈和三角形分别对应标准序贯估计方法和本发明序贯估计方法1000次蒙特卡洛模拟的RMSE。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，仅用于示例性说明，不能理解为对本专利的限制。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
下面结合附图和实施例对本发明的技术方案做进一步的说明。
实施例1
如图1所示，一种顾及解缠误差的序贯InSAR时序形变解算方法，所述的方法包括步骤如下：
S1：对历史N+1幅SAR影像序列，通过预处理得到其解缠差分干涉相位图上单个高相干点，根据单个高相干点的M个历史解缠相位观测值建立关于设计矩阵、历史形变相位、历史形变相位的协因数阵的线性方程组。
在一个具体的实施例中，所述的预处理是按时空基线集方法处理，得到M幅历史解缠差分干涉相位图，经过大气误差、轨道误差及DEM误差校正之后选取高相干点作为解算对象。
S2：根据最小二乘原理估计出历史形变相位，并引入M估计量来抵制历史数据集中解缠误差的影响，采用选权迭代法进行历史形变相位的求解，迭代得到历史形变相位、历史形变相位的协因数阵；图中的抗差估计指的是步骤S2中求解历史形变序列时所采用的估计方法。
S3：获取研究区域新增N1景SAR影像，将新增影像与邻近时刻的历史影像通过预处理，得到以历史解缠差分干涉相位图中选取的高相干点目标作为解算对象。
在一个具体的实施例中，所述的预处理是按时空基线集方法进行差分干涉得到M1幅新增解缠差分干涉相位图，经过大气误差、轨道误差及DEM误差校正之后选取高相干点作为解算对象。
S4：对新增M1个观测值建立包含历史形变相位及新增形变相位的观测方程；将观测方程写成误差方程的形式。
S5：根据历史形变相位中的先验信息及历史形变相位的协因数阵，建立虚拟观测方程，结合步骤S5所述的误差方程，根据贝叶斯平差原理，构造目标函数；对目标函数求解得到序贯平差的参数估值；
S6：引入M估计量来抵制新增观测数据中解缠误差对参数估值的影响，采用选权迭代法进行新增形变相位的求解，更新得到序贯平差的参数估值的协因数阵；图1中M-M抗差序贯平差指的是步骤S6中历史形变序列参数估值与新增观测值联合平差时所采用的估计方法。
S7：根据步骤S1～S6获得更新的序贯平差，然后通过对更新后的序贯平差拟合获取形变速率参数。
在一个具体的实施例中，所述的线性方程组；
Figure BDA0003852904880000071
其中，G为M×N的设计矩阵，与历史解缠差分干涉相位图组成的时空基线网络有关，其第l行代表的第l幅解缠差分干涉相位图由第s，t幅SAR影像干涉生成，0≤s＜t≤N，当s＞0时，对应矩阵G的第l行的第s个元素为-1；当s＝0时，矩阵G的第t个元素为1，其余元素均为0；
Figure BDA0003852904880000075
表示第n个时刻的形变相位，
Figure BDA0003852904880000076
表示历史解缠相位序列，n＝1,…,N，m＝1,…,M。
本实施例中，s大于0指的是G中某一行元素值的一般构成情况，当s＝0时，表示和第一景影像进行干涉，而第一景影像对应时刻的形变相位通常假设
Figure BDA0003852904880000072
不考虑在
Figure BDA0003852904880000077
内，即只需确定第t个元素的位置赋值为1即可(如矩阵G第一行)。
在一个具体的实施例中，若M＞N，并且设计矩阵G是列满秩的(即所有历史解缠差分干涉相位图都在一个集合里)，根据最小二乘原理估计出历史形变相位
Figure BDA0003852904880000073
P0为历史观测值L0的权矩阵；残差向量V0＝[v1,v2,…,vM]T通过
Figure BDA0003852904880000074
计算得到，其中下标0表示变量为初始化计算过程中的变量。
一般地，对于N+1幅SAR影像自由组合生成M个解缠差分干涉相位图，满足如下关系：
Figure BDA0003852904880000081
对于历史形变解算，当M>N时，观测方程存在多余观测，多余观测数越多，可提高参数估计的可靠性。
在一个具体的实施例中，步骤S2，考虑到最小二乘估计对粗差很敏感，引入M估计量来抵制历史数据集中解缠误差的影响，采用选权迭代法(IRLS)进行历史形变相位的求解，其迭代计算公式如下：
Figure BDA0003852904880000082
其中，
Figure BDA0003852904880000083
分别为第k次迭代计算的历史形变相位，历史残差向量，历史形变相位的协因数阵，历史观测值L0的等价权矩阵；w(k-1)为第k-1次迭代计算的权因子，可通过Huber权函数构造得到：
Figure BDA0003852904880000084
其中，wii(k)为权因子矩阵w中(i,i)位置在第k次迭代计算得到的权因子；
Figure BDA0003852904880000085
为标准化残差，c为常数，一般取为2.0～2.5。
本实施例也可以采用其他权函数构造w(k-1)，如IGGⅢ法。
当前后两次迭代计算得到的历史形变相位估值之差
Figure BDA0003852904880000086
满足迭代收敛条件(如
Figure BDA0003852904880000087
或历史形变相位估值之差小于0.0002等其他值也行)时，即获得历史形变相位的抗差解。
在一个具体的实施例中，所述的观测方程：
Figure BDA0003852904880000091
其中，G0为与历史形变相位序列有关的系数阵，G1为与新增形变相位序列有关的系数阵，二者的构造与新增解缠差分干涉相位图的网形有关；
例如，假设
Figure BDA0003852904880000092
由tN+1时刻获取的新增影像和tN+1时刻获取的历史影像差分干涉得到，则G0(1,N-1)＝-1，G1(1,1)＝1，其它元素都为0；
Figure BDA0003852904880000093
为新增的观测向量，N1表示新增影像的数量；M1表示新增影像与邻近历史影像生成的新增解缠差分干涉相位图的数量。
在一个具体的实施例中，将观测方程写成误差方程的形式，所述的误差方程表达式如下：
Figure BDA0003852904880000094
其中，V1为新增观测的残差向量，P1为新增观测的权矩阵；
Figure BDA0003852904880000099
表示序贯平差的参数估值。
在一个具体的实施例中，由于历史形变相位
Figure BDA00038529048800000910
含有先验信息
Figure BDA00038529048800000911
及其协因数阵
Figure BDA00038529048800000912
，建立虚拟观测方程：
Figure BDA0003852904880000095
Figure BDA0003852904880000096
其中，
Figure BDA00038529048800000913
为历史形变相位
Figure BDA00038529048800000914
的第二次估计，
Figure BDA00038529048800000915
为其对应的残差，
Figure BDA00038529048800000916
为其对应的权矩阵；结合所述的误差方程，根据贝叶斯平差原理，构造如下目标函数：
Figure BDA0003852904880000097
求解可得序贯平差的参数估值为：
Figure BDA0003852904880000098
在一个具体的实施例中，步骤S6，当获取到初始的序贯解时，为了抵制新增观测数据中解缠误差对参数估计的影响，同样引入M估计量来抵制新增观测数据中解缠误差对参数估值的影响，采用选权迭代法进行新增形变相位的求解，通过迭代求解以确定新增观测L1的等价权矩阵
Figure BDA0003852904880000101
最后，根据矩阵求逆引理，得到序贯解的递推形式：
Figure BDA0003852904880000102
序贯平差的参数估值
Figure BDA0003852904880000108
的协因数阵为：
Figure BDA0003852904880000103
其中，
Figure BDA0003852904880000104
为增益矩阵，
Figure BDA0003852904880000105
为新增观测对应的更新的协因数阵。
注意到新增观测的等价权矩阵
Figure BDA0003852904880000109
取代原始的权矩阵P1应用于参数的抗差序贯估计当中，因此可有效地抵制新增观测值中粗差的影响。
对于形变速率参数的动态更新，常规的方法是对形变相位用速率建模，进而按照序贯平差的方式更新形变速率，然而当地表运动先验知识缺乏的情况下，难以对地表形变有效建模，因此估计出的形变参数会有较大的偏差。本发明对于速率更新的策略是先按照上述步骤S1～S6获得更新的序贯平差，然后通过对序贯平差拟合获取形变速率参数。该方法无需对地表形变建模，因此在地表运动先验知识缺乏的情况下，也可有效的估计出形变速率。
在一个具体的实施例中，建立第i期新增影像获取时刻的序贯平差更新公式：
Figure BDA0003852904880000106
其中，
Figure BDA0003852904880000107
分别为第i期新增影像获取之前的形变相位序列的第i次更新结果以及第i期影像获取时刻的形变相位估计结果；Jx(i)表示第i期序贯更新公式中的增益矩阵；Li表示第i期的新增观测值；G0(i)表示第i期时空基线网络中与历史序贯平差有关的设计矩阵；G1(i)表示第i期时空基线网络中与新增序贯平差有关的设计矩阵；
Figure BDA00038529048800001010
表示第i期新增观测值对应的更新的协因数阵的逆。
因此，每次有新SAR影像集加入时，将新的观测数据结合前期平差的结果重复上述步骤S4～S6，可实现动态高效的形变参数解算，并且在序贯过程中抵制解缠误差的影响，保证形变参数更新结果的可靠性，为地质灾害实时监测和解译提供有力保障。
实施例2
基于实施例1所述的顾及解缠误差的序贯InSAR时序形变解算方法，其效果可以通过以下模拟实验进一步说明：
模拟实验描述：①根据2017/3/14～2020/2/3覆盖日本Echigo平原地区的45景Sentinel-1A数据的时空基线信息模拟得到169个解缠差分干涉相位图。如图2所示，以前35景影像构成的129个解缠差分干涉相位图作为历史观测数据进行序贯平差的初始化，后10景影像构成的40个解缠差分干涉相位图作为新增观测数据用于序贯平差的序贯更新；②分别用三种不同的形变模式来模拟形变信号，在模拟的形变信号中加入标准差4mm的随机噪声作为形变解算的观测值；③为了评估解缠误差对InSAR序贯平差序贯估计的影响，以及本发明所述的序贯InSAR时序形变解算方法的抗差性能，解缠误差通过对干涉相位添加±2π来实现。
图3分别显示了标准序贯估计方法以及本发明所述的序贯InSAR时序形变解算方法估计的序贯平差结果。从图3(a、c、e)左栏可以看出，当历史和新增观测数据中都不存在解缠误差时，两种方法获取的序贯平差结果几乎一致，且与真实的形变信号非常接近；从图3(b、d、f)可以看出，当在历史观测数据和新增观测数据分别加入15％和10％的解缠误差时，对于标准序贯估计方法，其历史序贯平差估计结果受解缠误差影响严重偏离真实形变情况，以扭曲的历史序贯平差结果为基础进行标准序贯更新的序贯平差估计结果仍严重偏离真值；本发明方法保证解算的历史序贯平差参数估值可靠，同时采用等价权抑制新增观测值中解缠误差的影响，以此为基础序贯更新的序贯平差结果接近真实的形变情况。
为了得到更多量化的结果，1000次蒙特卡洛模拟被用来评估本发明所述的方法在不同数量解缠误差下的抗差效果。本次模拟实验所设置的分布概率为0％，3％，6％，9％，12％，15％，18％，21％的历史解缠误差以及0％，5％，10％，15％，20％，25％的新增解缠误差。1000次蒙特卡洛模拟的均方根误差(RMSE)被用来衡量最终结果。如图4所示，本发明所述的方法估计结果的精度整体上优于标准序贯估计方法，在解缠误差比例较大的情况下，本发明所述的方法仍有较好的表现，说明本发明所述的方法有较好的抗差性，可以有效抵制一定数量解缠误差的影响。此外，我们注意到当解缠误差比例接近于0％的时候，标准序贯估计方法优于本发明所述的方法，这表明本发明所述的方法以牺牲很小的最优性为代价，换取较强的抗差性能，从而获得可靠的参数估值。
实施例3
一种计算机设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，所述的处理器执行所述的计算机程序时，实现如实施例1所述的顾及解缠误差的序贯InSAR时序形变解算方法的步骤。
其中，存储器和处理器采用总线方式连接，总线可以包括任意数量的互联的总线和桥，总线将一个或多个处理器和存储器的各种电路连接在一起。总线还可以将诸如外围设备、稳压器和功率管理电路等之类的各种其他电路连接在一起，这些都是本领域所公知的，因此，本文不再对其进行进一步描述。总线接口在总线和收发机之间提供接口。收发机可以是一个元件，也可以是多个元件，比如多个接收器和发送器，提供用于在传输介质上与各种其他装置通信的单元。经处理器处理的数据通过天线在无线介质上进行传输，进一步，天线还接收数据并将数据传送给处理器。
本实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，所述的计算机程序被处理器执行时，实现如实施例1所述的顾及解缠误差的序贯InSAR时序形变解算方法的步骤。
即，本领域技术人员可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，该程序存储在一个存储介质中，包括若干指令用以使得一个设备(可以是单片机，芯片等)或处理器(processor)执行本申请各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (10)

1.一种顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：所述的方法包括步骤如下：
S1：对历史N+1幅SAR影像序列，通过预处理得到其解缠差分干涉相位图上单个高相干点，根据单个高相干点的M个历史解缠相位观测值建立关于设计矩阵、历史形变相位、历史形变相位的协因数阵的线性方程组；
S2：根据最小二乘原理估计出历史形变相位，并引入M估计量来抵制历史数据集中解缠误差的影响，采用选权迭代法进行历史形变相位的求解，迭代得到历史形变相位、历史形变相位的协因数阵；
S3：获取研究区域新增N1景SAR影像，将新增影像与邻近时刻的历史影像通过预处理，以历史解缠差分干涉相位图中选取的高相干点目标作为解算对象；
S4：对新增M1个观测值建立包含历史形变相位及新增形变相位的观测方程；将观测方程写成误差方程的形式；
S5：根据历史形变相位中的先验信息及历史形变相位的协因数阵，建立虚拟观测方程，结合步骤S4所述的误差方程，根据贝叶斯平差原理，构造目标函数；对目标函数求解得到序贯平差的参数估值；
S6：引入M估计量来抵制新增观测数据中解缠误差对参数估值的影响，采用选权迭代法进行新增形变相位的求解，更新得到序贯平差的参数估值的协因数阵；
S7：根据步骤S1～S6获得更新的序贯平差，然后通过对更新后的序贯平差拟合获取形变速率参数。
2.根据权利要求1所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：所述的预处理是按时空基线集方法处理，得到M幅解缠差分干涉相位图，经过大气误差、轨道误差及DEM误差校正之后选取高相干点作为解算对象。
3.根据权利要求1所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：所述的线性方程组；
Figure FDA0004131869880000021
其中，G为M×N的设计矩阵，与历史解缠差分干涉相位图组成的时空基线网络有关，其第l行代表的第l幅解缠差分干涉相位图由第s，t幅SAR影像干涉生成，0≤s＜t≤N，当s＞0时，对应矩阵G的第l行的第s个元素为-1；当s＝0时，矩阵G的第t个元素为1，其余元素均为0；
Figure FDA0004131869880000022
表示第n个时刻的形变相位，
Figure FDA0004131869880000023
表示历史解缠相位序列，n＝1，...，N，m＝1，...，M。
4.根据权利要求3所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：若M＞N，并且设计矩阵G是列满秩的，根据最小二乘原理估计出历史形变相位
Figure FDA0004131869880000024
P0为历史观测值L0的权矩阵；残差向量V0＝[v1，v2，...，vM]T通过
Figure FDA0004131869880000025
计算得到，其中下标0表示变量为初始化计算过程中的变量。
5.根据权利要求3所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：步骤S2，所述的采用选权迭代法进行历史形变相位的求解，其迭代计算公式如下：
Figure FDA0004131869880000026
其中，
Figure FDA0004131869880000027
分别为第k次迭代计算的历史形变相位，历史残差向量，历史形变相位的协因数阵，历史观测值L0的等价权矩阵；w(k-1)为第k-1次迭代计算的权因子；
当前后两次迭代计算得到的历史形变相位估值之差
Figure FDA0004131869880000028
满足迭代收敛条件时，即获得历史形变相位的抗差解。
6.根据权利要求5所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：所述的观测方程：
Figure FDA0004131869880000031
其中，G0为与历史形变相位序列有关的系数阵，G1为与新增形变相位序列有关的系数阵，二者的构造与新增解缠差分干涉相位图的网形有关；
假设
Figure FDA00041318698800000312
由tN+1时刻获取的新增影像和tN+1时刻获取的历史影像差分干涉得到，则G0(1，N-1)＝-1，G1(1，1)＝1，其它元素都为0；
Figure FDA0004131869880000032
为新增的观测向量，N1表示新增影像的数量；M1表示新增影像与邻近历史影像生成的新增解缠差分干涉相位图的数量。
7.根据权利要求6所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：所述的误差方程表达式如下：
Figure FDA0004131869880000033
其中，V1为新增观测的残差向量，P1为新增观测的权矩阵；
Figure FDA0004131869880000034
表示序贯平差的参数估值。
8.根据权利要求6所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：所述的虚拟观测方程：
Figure FDA0004131869880000035
Figure FDA0004131869880000036
其中，
Figure FDA0004131869880000037
为历史形变相位
Figure FDA0004131869880000038
的第二次估计，
Figure FDA0004131869880000039
为其对应的残差，
Figure FDA00041318698800000310
为其对应的权矩阵；结合所述的误差方程，根据贝叶斯平差原理，构造如下目标函数：
Figure FDA00041318698800000311
式中，V1为新增观测的残差向量，P1为新增观测的权矩阵；
求解可得序贯平差的参数估值为：
Figure FDA0004131869880000041
9.根据权利要求8所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：步骤S6，通过迭代求解以确定新增观测向量L1的等价权矩阵
Figure FDA0004131869880000042
最后，根据矩阵求逆引理，得到序贯解的递推形式：
Figure FDA0004131869880000043
序贯平差的参数估值
Figure FDA0004131869880000044
的协因数阵为：
Figure FDA0004131869880000045
其中，
Figure FDA0004131869880000046
为增益矩阵，
Figure FDA0004131869880000047
为新增观测对应的更新的协因数阵。
10.根据权利要求2～9任一项所述的顾及解缠误差的序贯InSAR时序形变解算方法，其特征在于：建立第i期新增影像获取时刻的序贯平差更新公式：
Figure FDA0004131869880000048
其中，
Figure FDA0004131869880000049
分别为第i期新增影像获取之前的形变相位序列的第i次更新结果以及第i期影像获取时刻的形变相位估计结果；Jx(i)表示第i期序贯更新公式中的增益矩阵；Li表示第i期的新增观测值；G0(i)表示第i期时空基线网络中与历史序贯平差有关的设计矩阵；G1(i)表示第i期时空基线网络中与新增序贯平差有关的设计矩阵；
Figure FDA00041318698800000410
表示第i期新增观测值对应的更新的协因数阵的逆。
CN202211137798.5A 2022-09-19 2022-09-19 一种顾及解缠误差的序贯InSAR时序形变解算方法 Active CN115453534B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211137798.5A CN115453534B (zh) 2022-09-19 2022-09-19 一种顾及解缠误差的序贯InSAR时序形变解算方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211137798.5A CN115453534B (zh) 2022-09-19 2022-09-19 一种顾及解缠误差的序贯InSAR时序形变解算方法

Publications (2)

Publication Number Publication Date
CN115453534A CN115453534A (zh) 2022-12-09
CN115453534B true CN115453534B (zh) 2023-05-16

Family

ID=84304584

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211137798.5A Active CN115453534B (zh) 2022-09-19 2022-09-19 一种顾及解缠误差的序贯InSAR时序形变解算方法

Country Status (1)

Country Link
CN (1) CN115453534B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116299245B (zh) * 2023-05-11 2023-07-28 中山大学 一种时序InSAR形变速率结果自适应镶嵌校正方法

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107909606A (zh) * 2017-12-25 2018-04-13 南京市测绘勘察研究院股份有限公司 一种sar图像配准联系点粗差剔除方法
CN109061641B (zh) * 2018-07-06 2020-01-17 中南大学 一种基于序贯平差的InSAR时序地表形变监测方法
EP3866105A1 (en) * 2020-02-17 2021-08-18 Paris Sciences et Lettres - Quartier Latin Method for processing insar images to extract ground deformation signals
CN112986993B (zh) * 2021-02-07 2022-10-25 同济大学 一种基于空间约束的InSAR形变监测方法

Also Published As

Publication number Publication date
CN115453534A (zh) 2022-12-09

Similar Documents

Publication Publication Date Title
CN107193003B (zh) 一种稀疏奇异值分解扫描雷达前视成像方法
Kim Development of track to track fusion algorithms
CN102597701B (zh) 用于补偿错误测量的系统和方法
CN109116293B (zh) 一种基于离格稀疏贝叶斯的波达方向估计方法
CN113792411B (zh) 一种基于中心误差熵准则无迹卡尔曼滤波的航天器姿态确定方法
CN111046591B (zh) 传感器幅相误差与目标到达角度的联合估计方法
CN115453534B (zh) 一种顾及解缠误差的序贯InSAR时序形变解算方法
CN113176532B (zh) 基于波束空间的阵列误差和阵元失效的dnn鲁棒性doa估计方法、装置及存储介质
CN113655444B (zh) 一种阵元失效下基于重加权先验的mimo雷达doa估计方法
CN112766304A (zh) 一种基于稀疏贝叶斯学习的机动阵列方位估计方法
CN114626307B (zh) 一种基于变分贝叶斯的分布式一致性目标状态估计方法
Kulikova et al. NIRK-based accurate continuous–discrete extended Kalman filters for estimating continuous-time stochastic target tracking models
CN115453528A (zh) 基于快速sbl算法实现分段观测isar高分辨成像方法及装置
CN112083457A (zh) 一种神经网络优化的imm卫星定位导航方法
CN114139109A (zh) 一种目标跟踪方法、系统、设备、介质及数据处理终端
CN110146084B (zh) 面向卫星故障的多星编队系统分布式协同导航滤波方法
CN114567288B (zh) 基于变分贝叶斯的分布协同非线性系统状态估计方法
CN114070262B (zh) 附加扰动的集成混合集合Kalman滤波天气预报同化方法及装置
CN113740802B (zh) 以自适应噪声估计进行矩阵补全的信号源定位方法和系统
CN110716219A (zh) 一种提高定位解算精度的方法
CN109035301B (zh) 一种基于斥力模型修正随机矩阵算法的群组目标跟踪方法
CN114741659A (zh) 一种自适应模型在线重构建鲁棒滤波方法、设备及系统
CN109684771B (zh) 基于交互式多模型的机动目标状态预测优化方法
CN118519174B (zh) 面向高级接收机自主完好性监测的智能选星方法
CN110988790B (zh) 一种无源目标的定位方法及装置

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