CN113552565A - 针对sar数据高噪声及大梯度变化区域的相位解缠方法 - Google Patents

针对sar数据高噪声及大梯度变化区域的相位解缠方法 Download PDF

Info

Publication number
CN113552565A
CN113552565A CN202110823529.3A CN202110823529A CN113552565A CN 113552565 A CN113552565 A CN 113552565A CN 202110823529 A CN202110823529 A CN 202110823529A CN 113552565 A CN113552565 A CN 113552565A
Authority
CN
China
Prior art keywords
pixel
phase
matrix
covariance
square root
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
CN202110823529.3A
Other languages
English (en)
Other versions
CN113552565B (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.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
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 China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN202110823529.3A priority Critical patent/CN113552565B/zh
Publication of CN113552565A publication Critical patent/CN113552565A/zh
Application granted granted Critical
Publication of CN113552565B publication Critical patent/CN113552565B/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

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Image Processing (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法，属于雷达卫星SAR数据领域，在高噪声区域及大梯度变化区域相位解缠问题中具有较好地应用。通过结合低通滤波技术及奇异值分解有效地抑制高频噪声，并且对频率估计结果的异常值进行检测和修正，提高相位梯度估计的精度；基于相位梯度估计结果采用自适应平方根无迹Kalman滤波技术结合最大堆解缠策略完成相位解缠。本方法步骤简单，使用效果好，对于高噪声区域和地形大梯度变化区域的相位解缠具有更好的解缠精度和模型鲁棒性，有利于提高最终InSAR数据处理结果的精度。

Description

针对SAR数据高噪声及大梯度变化区域的相位解缠方法
技术领域
本发明涉及一种相位解缠方法，尤其适用于处理雷达卫星SAR数据使用的一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法。
背景技术
合成孔径雷达干涉(Interferometric SyntheticAperture Radar,InSAR)测量技术近年来已广泛应用于地形三维重建和由地震活动的物理模型反演，尽管目前已经存在几十种相位解缠方法，但是在噪声量大的区域相位解缠方法的鲁棒性一般较差，在地形梯度大的区域相位解缠结果误差一般较大。这几十种相位解缠方法可以大致分为两类：一类是单基线相位解缠方法，另一类是多基线相位解缠方法。由于单基线相位解缠方法受相位连续性假设的限制，在地形梯度大的区域难以获得理想的结果，虽然多基线相位解缠方法可以摆脱相位连续性假设的限制，但是在解缠效率和基线配比等方面仍然存在一些困难。
由于修正矩阵束模型具有一定的频率估计精度和效率，所以本发明保留了修正矩阵束模型的框架，在其奇异值处理过程中引入了低通滤波技术，并且可以根据干涉图像素均方差值自适应选择局部频率估计窗口的大小，尽可能地实现方位向和距离向相位梯度的准确估计。这种低通滤波技术是由Butterworth滤波器来实现的，其特点是通频带内的频率响应曲线最大限度平坦，没有纹波，而在阻频带则逐渐下降为零，这意味着Butterworth低通滤波器可以减少高频噪声。此外，我们也引入了自适应平方根无迹Kalman技术到相位解缠问题中，平方根无迹Kalman算法只需要误差协方差的平方根，而不需要完整的协方差阵，其采用误差协方差阵的平方根形式替代协方差阵进行递推计算，保证了误差协方差阵的非负定性，从而提高了数值的稳定性，并且在计算中只需要存储和运算平方根因数，这样可以降低计算负担，以实现对无迹Kalman相位解缠技术性能的提高。
发明内容
针对现有技术的不足之处，提供一种步骤简单，使用效果好的一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法
为实现上述技术目的，本发明的针对高噪声及大梯度变化区域的相位解缠方法，所包括以下步骤：
步骤S1，收集需要提取高程量或微小形变量地区的SAR数据，对SAR数据进行预处理得到干涉相位图；通过改进的修正矩阵束模型对干涉相位图所有像元的局部频率进行估计，改进的修正矩阵束是在修正矩阵束模型的奇异值处理步骤中替用一阶Butterworth低通滤波器对分解后得到的每一个奇异值赋予不同的权值来实现信号重构，从而获取准确的局部频率，以此对干涉相位图中所有像元的方位向和距离向的相位梯度进行估计：局部频率估计窗口的大小根据干涉图像元相位均方差自适应选择，窗口包括小的窗口和大的窗口，像元的干涉相位均方差反映条纹的密集程度，条纹密集则均方差大，利用小的窗口进行局部频率估计以免破坏局部条纹平稳的假设；条纹稀疏则均方差小，利用大的窗口进行局部频率估计，既满足局部条纹平稳的假设，又增加了样本点数量，提高了频率估计结果的准确性；
步骤S2，由于局部频率估计结果不存在大的跳变，当局部频率估计结果出现明显的不连续时，则认定该局部像元的频率估计结果存在误差，需要进行相应地修正，利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测；
Figure BDA0003172748920000021
Figure BDA0003172748920000022
式中M为与判定连续性窗口大小相关的参数，例如M＝3，即表示判定连续性窗口大小为(2M+1)×(2M+1)＝7×7；u和v表示局部窗口内像元偏离k位置中心像元的程度，
Figure BDA0003172748920000023
和
Figure BDA0003172748920000024
分别表示干涉相位图中k位置中心像元的频率估计结果在距离向和方位向的可信值，fx和fy分别是距离向和方位向上的频率估计结果；
步骤S3，根据求得的可信值
Figure BDA0003172748920000025
和
Figure BDA0003172748920000026
设置阈值，分别设置为
Figure BDA0003172748920000027
和
Figure BDA0003172748920000028
其中
Figure BDA0003172748920000029
和
Figure BDA00031727489200000210
为
Figure BDA00031727489200000211
的最大值和
Figure BDA00031727489200000212
的最大值，利用下式将干涉相位图中k位置中心像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值：
Figure BDA00031727489200000213
Figure BDA00031727489200000214
其中
Figure BDA00031727489200000215
和
Figure BDA00031727489200000216
表示修正后干涉相位图中k位置中心像元在方位向和距离向的窗口范围内所有像元频率的平均值，从而解决了局部频率估计结果存在异常值的问题；
步骤S4，建立无迹Kalman滤波相位解缠观测方程的自适应平方根模型，并利用无迹Kalman滤波相位解缠观测方程的自适应平方根模型获得相位解缠起始像元；首先计算自适应平方根无迹Kalman相位解缠模型的Sigma采样点及其权值系数；将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传播后均值和误差协方差的平方根，利用误差协方差阵的平方根形式替代协方差阵进行递推计算，从而保证了误差协方差阵的非负定性，提高了数值的稳定性，并且在计算中只需要存储和运算平方根因数，从而降低计算负担；
步骤S5，计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差；其中采用测量噪声等价协方差平方根阵原理，引入自适应因子，实现对测量噪声协方差平方根的异常值分离和自适应修正，进而提高平方根无迹Kalman相位解缠模型的精度；
步骤S6，计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的滤波增益；利用求得的Kalman滤波增益矩阵对状态估计及其协方差平方根进行更新，获取鲁棒性更好的无迹Kalman相位解缠模型；
步骤S7，通过基于堆排序的有效质量指导策略对干涉相位图进行逐个相位展开，将干涉复数信号所有像元的相位质量，即质量图作为指导策略，以最大堆作为数据结构，通过不断更新最大堆的节点来实现对质量图的排序，保证了相位解缠过程一直从高质量像元向低质量像元计算；干涉相位图中每一个像元的质量是由相干系数表示，相干系数位于区间[0，1]，相干系数越大则代表像元的质量越高，相反若相干系数越小则代表像元的质量越低，所有像元的相干系数组合在一起即为质量图；重复步骤S4-步骤S6，直到干涉相位图中所有像元都被解缠。
进一步，使用改进的修正矩阵束模型对干涉图距离向和方位向梯度进行估计：
步骤a1，计算干涉图中所有像元干涉相位的均方差ξ，公式为
Figure BDA0003172748920000031
其中l×l为计算像元干涉相位均方差的一个窗口大小，即判断条纹密集程度的一个窗口大小，一般可设置l＝5，As表示局部窗口内第s个像元的干涉相位复数的振幅；
Figure BDA0003172748920000032
为A的平均值；
步骤a2，使用公式(1)求得的干涉相位均方差ξ，根据ξ值自适应选择待估像元的局部频率估计窗口大小，其中7×7或9×9为小的窗口，17×17或19×19为大的窗口，公式为：
Figure BDA0003172748920000033
其中，H×H表示干涉相位图中k位置中心像元进行局部频率估计计算的窗口大小；
步骤a3，假设Ik是以k像元为中心，窗口大小为H×H的干涉相位样本，利用公式：Ik＝UDVT对信号矩阵Ik进行奇异值分解，其中U和V表示H阶正交(酉)矩阵；D表示对应的特征值矩阵，D＝diag(σ1,σ2,…σH)，且σ1≥σ2≥…≥σH≥0，σi,i＝1,2,…,H表示信号矩阵Ik奇异值分解后对角线上的特征值；
步骤a4，由一阶Butterworth低通滤波器计算传递函数，即对信号矩阵Ik分解后对角线上的每个特征值赋予不同的权值以降低高频噪声对信号重构的影响，具体计算过程为：
Figure BDA0003172748920000041
式中Wt表示特征值矩阵D对角线上第t个奇异值即将被分配的权值系数，σa表示前t个奇异值中第a个奇异值，σt表示特征值矩阵D对角线上第t个奇异值；
步骤a5，利用公式(4)计算处理后的干涉相位样本矩阵
Figure BDA0003172748920000042
即尽量剔除信号中高频噪声成分，尽量完整地保留有用的相位信息：
Figure BDA0003172748920000043
步骤a6，利用公式(5)根据处理后的干涉相位样本矩阵
Figure BDA0003172748920000044
构建母Hankel矩阵以及两个子Hankel矩阵：
Figure BDA0003172748920000045
其中
Figure BDA0003172748920000046
及
Figure BDA0003172748920000047
分别为母Hankel矩阵以及两个子Hankel矩阵；
步骤a7，利用公式(6)对母Hankel矩阵
Figure BDA0003172748920000048
进行奇异值分解：
Figure BDA0003172748920000049
其中
Figure BDA00031727489200000410
和
Figure BDA00031727489200000411
仅具有与有用信号相关的主特征值对应的信息；
步骤a8，通过奇异值分解，并结合公式(5)和公式(6)得到降维处理后的相位数据信号，由于相位数据矩阵具有低秩结构，经过降维处理后以抑制噪声，推导获得如下公式：
Figure BDA00031727489200000412
其中
Figure BDA00031727489200000413
及
Figure BDA00031727489200000414
分别是通过母Hankel矩阵
Figure BDA00031727489200000415
以及两个子Hankel矩阵
Figure BDA00031727489200000416
和
Figure BDA00031727489200000417
重构后的信号矩阵；
步骤a9，利用公式(8)表示以k像元为中心窗口内样本干涉相位的距离向和方位向的局部频率：
Figure BDA00031727489200000418
其中
Figure BDA0003172748920000051
和
Figure BDA0003172748920000052
分别为
Figure BDA0003172748920000053
和
Figure BDA0003172748920000054
的伪逆矩阵，conj(·)为矩阵的共轭相乘因子；
Figure BDA0003172748920000055
和
Figure BDA0003172748920000056
分别是k位置像元在距离向和方位向上的频率估计结果；
步骤a10，由于局部频率估计结果不存在大的跳变，当频率估计结果出现明显的不连续时，则认定该局部像元的频率估计结果存在误差，需要进行相应地修正，利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测；
Figure BDA0003172748920000057
式中M为与判定连续性窗口大小相关的参数，例如M＝3，即表示判定连续性窗口大小为(2M+1)×(2M+1)＝7×7；u和v表示局部窗口内像元偏离k位置中心像元的程度，
Figure BDA0003172748920000058
和
Figure BDA0003172748920000059
分别表示干涉相位图中k位置中心像元的频率估计结果在方位向和距离向的可信值，fx和fy分别是像元在距离向和方位向上的频率估计结果；
步骤a11，根据求得的可信值
Figure BDA00031727489200000510
和
Figure BDA00031727489200000511
设置阈值，分别设置为
Figure BDA00031727489200000512
和
Figure BDA00031727489200000513
其中
Figure BDA00031727489200000514
和
Figure BDA00031727489200000515
为
Figure BDA00031727489200000516
的最大值和
Figure BDA00031727489200000517
的最大值，利用下式计算该像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值：
Figure BDA00031727489200000518
其中fx和fy分别是像元在距离向和方位向上的频率估计结果；
Figure BDA00031727489200000519
和
Figure BDA00031727489200000520
表示修正后像元k在距离向和方位向的频率估计结果，以解决局部频率估计结果存在异常值的问题；
步骤a12，根据步骤a11获得局部频率可以获得像素k的距离向和方位向梯度，公式如下：
Figure BDA00031727489200000521
其中
Figure BDA00031727489200000522
和
Figure BDA00031727489200000523
分别为像素k的距离向和方位向梯度；
重复步骤a2-步骤a12，直至估计出干涉相位图中所有像元距离向和方位向的梯度。
进一步，建立自适应平方根无迹Kalman滤波相位解缠模型，获得相位解缠起始像元的步骤为：
步骤b1，利用下式计算自适应平方根无迹Kalman相位解缠模型解缠起始像元的Sigma采样点及其权值系数，公式为：
Figure BDA0003172748920000061
其中ρi为起始像元的Sigma采样点，
Figure BDA0003172748920000062
是第k-1位置解缠相位的状态估计值，Sk-1是第k-1位置像元的误差协方差的平方根，N表示状态变量的维数(N＝1)；{·}i表示括号内矩阵的第i列，η表示尺度函数；
Figure BDA0003172748920000063
其中
Figure BDA0003172748920000064
和
Figure BDA0003172748920000065
为权值系数，η＝τ2(N+α)-N表示尺度函数，τ，α，ζ用来调节Sigma点，通常τ＝0.01,ζ＝2,α＝0，τ决定Sigma点在均值附近的散布程度，α为冗余量；
步骤b2，将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传递后均值和误差协方差的平方根，公式为：
Figure BDA0003172748920000066
其中f[·]表示状态方程转移函数；符号qr{·}和cholupdate{·}表示QR分解运算和Cholesky因子一阶更新运算，
Figure BDA0003172748920000067
表示第k个像元Sigma点预测值；
Figure BDA0003172748920000068
表示干涉相位图第k个像元的状态估计；Sψψ表示干涉相位图第k个像元状态预测误差协方差平方根；Q(k-1)是干涉相位图第k个像元的相位梯度估计误差方差。
进一步，计算自适应平方根无迹Kalman滤波相位解缠模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差；采用测量噪声等价协方差平方根阵原理对测量误差方差的平方根进行自适应修正的具体步骤为：
步骤c1，计算待解缠像元的量测预测值，公式为：
Figure BDA0003172748920000071
其中，
Figure BDA0003172748920000072
表示k-1像元的Sigma点经传递后得到的k像元的Sigma点预测值；
Figure BDA0003172748920000073
表示待解缠像元的量测预测值；h[·]为观测方程系数。
步骤c2，利用公式(16)和公式(17)计算待解缠像元的量测预测值协方差的平方根和状态与量测间的互协方差：
Figure BDA0003172748920000074
Figure BDA0003172748920000075
其中
Figure BDA0003172748920000076
表示待解缠像元k的量测预测值协方差的平方根；
Figure BDA0003172748920000077
表示待解缠像元k的状态与量测间的互协方差；
步骤c3，利用下式计算待解缠像元的预测残差信息：
Figure BDA0003172748920000078
式中
Figure BDA0003172748920000079
为k位置像元的预测残差，
Figure BDA00031727489200000710
表示k位置像元的真实缠绕相位，
Figure BDA00031727489200000711
为待解缠像元的量测预测值；
步骤c4，利用下式计算待解缠像元测量噪声协方差的等价协方差阵的对角线元素：
Figure BDA00031727489200000712
其中Ri表示测量噪声协方差阵R的第i个对角线元素；
Figure BDA00031727489200000713
表示测量噪声协方差的等价协方差阵第i个对角线元素；U0和U1为小数值的常数阈值，这里可设置U0＝0.45,U1＝3；vi为标准化预测残差。
Figure BDA00031727489200000714
表示预测观测值协方差阵的第i个对角线元素；median(·)表示取中位数算子。
步骤c5，将公式(19)求得的对角线元素
Figure BDA0003172748920000081
组成的测量噪声协方差的等价协方差阵
Figure BDA0003172748920000082
带入式(16)中替换测量噪声协方差阵R重新更新计算量测预测值的协方差平方根。
进一步，计算自适应平方根无迹Kalman滤波相位解缠模型的滤波增益，对状态估计及其协方差平方根进行更新，公式为
Figure BDA0003172748920000083
其中υk表示干涉相位图第k个像元的增益矩阵；
Figure BDA0003172748920000084
表示干涉相位图第k个像元的状态估计值；S(k)表示干涉相位图第k个像元状态估计的协方差平方根。
有益效果：
本发明通过在修正矩阵束模型中嵌入一个Butterworth低通滤波器，并可以根据干涉图的相位均方差自适应调整局部频率估计分析窗的大小，然后对频率估计结果异常值进行修正以获取高精度的距离向和方位向梯度；通过自适应平方根无迹Kalman相位解缠方法结合最大堆解缠策略得到解缠相位。相比于其它常规相位解缠方法，本方法对于高噪声区域和地形大梯度变化区域的相位解缠具有更好的解缠精度和模型鲁棒性，有利于提高最终InSAR数据处理结果的精度。
附图说明
图1是本发明的针对高噪声及大梯度变化区域的相位解缠方法流程图流程图；
图2为本发明研究区的地理位置区域图；
图3为本发明研究区的干涉相位图；
图4为本发明距离向的相位梯度估计结果示意图；
图5为本发明方位向的相位梯度估计结果示意图；
图6为本发明修正后的距离向相位梯度估计结果示意图；
图7为本发明修正后的方位向相位梯度估计结果示意图；
图8为本发明针对高噪声及大梯度变化区域的相位解缠方法的解缠结果示意图；
图9为MCF方法的相位解缠结果示意图；
图10为SNAPHU方法的相位解缠结果示意图；
图11为无迹Kalman相位解缠方法的相位解缠结果示意图；
图12为由高精度DEM生成的模拟解缠相位图。
具体实施方式
下面结合附图，对本发明作进一步的阐述：
如图1所示，本发明的一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法，输入信息为研究区的地理位置区域图的干涉相位图，其中实际地理位置区域图如图2所示，输入的干涉相位图如图3所示，包括一个干涉像对经过粗配准、精配准、去平/去地形、复共轭相乘等处理过程，所述方法包括以下步骤：
步骤S1，使用改进的修正矩阵束模型对干涉图所有像元的距离向和方位向梯度进行估计。
步骤S1.1，计算干涉图中所有像元干涉相位的均方差ξ，公式为：
Figure BDA0003172748920000091
其中l×l为计算像元干涉相位均方差的一个窗口大小，即判断条纹密集程度的一个窗口大小，一般可设置l＝5；As表示局部窗口内第s个像元的干涉相位复数的振幅；
Figure BDA0003172748920000092
为A的平均值；
步骤S1.2，使用公式(1)求得的干涉相位均方差ξ，根据ξ值自适应选择待估像元的局部频率估计窗口大小，其中7×7或9×9为小的窗口，17×17或19×19为大的窗口，公式为：
Figure BDA0003172748920000093
其中，H×H表示干涉相位图中k位置中心像元进行局部频率估计计算的窗口大小；
步骤S1.3，假设Ik是以k像元为中心，窗口大小为H×H的干涉相位样本，利用公式：Ik＝UDVT对信号矩阵Ik进行奇异值分解，其中U和V表示H阶正交(酉)矩阵；D表示对应的特征值矩阵，D＝diag(σ1,σ2,…σH)，且σ1≥σ2≥…≥σH≥0，σi,i＝1,2,…,H表示信号矩阵Ik奇异值分解后对角线上的特征值。
步骤S1.4，由一阶Butterworth低通滤波器计算传递函数，即对信号矩阵Ik分解后对角线上的每个特征值赋予不同的权值以降低高频噪声对信号重构的影响，具体计算过程为：
Figure BDA0003172748920000094
式中Wt表示特征值矩阵D对角线上第t个奇异值即将被分配的权值系数，σa表示前t个奇异值中第a个奇异值，σt表示特征值矩阵D对角线上第t个奇异值；
步骤S1.5，利用公式(4)计算处理后的干涉相位样本矩阵
Figure BDA0003172748920000101
即尽量剔除信号中高频噪声成分，尽量完整地保留有用的相位信息：
Figure BDA0003172748920000102
步骤S1.6，利用公式(5)根据处理后的干涉相位样本矩阵构建母Hankel矩阵以及两个子Hankel矩阵：
Figure BDA0003172748920000103
其中
Figure BDA0003172748920000104
及
Figure BDA0003172748920000105
分别为母Hankel矩阵以及两个子Hankel矩阵；
步骤S1.7，利用公式(6)对母Hankel矩阵
Figure BDA0003172748920000106
进行奇异值分解：
Figure BDA0003172748920000107
其中
Figure BDA0003172748920000108
和
Figure BDA0003172748920000109
仅具有与有用信号相关的主特征值对应得信息；
步骤S1.8，通过奇异值分解，结合公式(5)和公式(6)可以得到降维处理后的相位数据信号，由于相位数据矩阵具有低秩结构，经过降维处理后可以抑制噪声，推导获得如下公式：
Figure BDA00031727489200001010
其中
Figure BDA00031727489200001011
及
Figure BDA00031727489200001012
分别是通过母Hankel矩阵
Figure BDA00031727489200001013
以及两个子Hankel矩阵
Figure BDA00031727489200001014
和
Figure BDA00031727489200001015
重构后的信号矩阵；
步骤S1.9，利用公式(8)表示以k像元为中心窗口内样本干涉相位的距离向和方位向的局部频率：
Figure BDA00031727489200001016
其中
Figure BDA00031727489200001017
和
Figure BDA00031727489200001018
分别为
Figure BDA00031727489200001019
和
Figure BDA00031727489200001020
的伪逆矩阵，conj(·)为矩阵的共轭相乘因子；
Figure BDA00031727489200001021
和
Figure BDA00031727489200001022
分别是k位置像元在距离向和方位向上的频率估计结果；
步骤S1.10，根据步骤S1.9获得局部频率可以获得像素k的距离向和方位向梯度，公式如下：
Figure BDA0003172748920000111
其中
Figure BDA0003172748920000112
和
Figure BDA0003172748920000113
分别为像素k的距离向和方位向梯度；
重复步骤S1.2-S1.10，直至估计出所有像元距离向和方位向的梯度。方位向和距离向的相位梯度估计结果如图4和图5所示。
步骤S2，根据公式(9)，从距离向和方位向上分别对频率估计结果的连续性进行计算和检测。阈值分别设置为
Figure BDA0003172748920000114
和
Figure BDA0003172748920000115
其中
Figure BDA0003172748920000116
和
Figure BDA0003172748920000117
为
Figure BDA0003172748920000118
的最大值和
Figure BDA0003172748920000119
的最大值。当可信值超过阈值时，认为该像元的频率估计结果存在异常，需要进行修正。
Figure BDA00031727489200001110
式中M为与判定连续性窗口大小相关的参数，例如M＝3，即表示判定连续性窗口大小为(2M+1)×(2M+1)＝7×7；u和v表示局部窗口内像元偏离k位置中心像元的程度，
Figure BDA00031727489200001111
和
Figure BDA00031727489200001112
分别表示干涉相位图中k位置中心像元的频率估计结果在方位向和距离向的可信值，fx和fy分别是像元在距离向和方位向上的频率估计结果；
步骤S3，从距离向和方位向上分别对频率估计结果的异常值进行修正，异常像元的频率估计结果优化为以k位置像元为中心(2M+1)×(2M+1)范围内所有像元频率的平均，公式为
Figure BDA00031727489200001113
其中fx和fy分别是像元在距离向和方位向上的频率估计结果；
Figure BDA00031727489200001114
和
Figure BDA00031727489200001115
表示修正后像元在方位向和距离向的频率估计结果。由公式(10)的频率结果重新根据公式(28)计算修正后的方位向和距离向相位梯度。本发明实例干涉图在方位向和距离向上的相位梯度修正后的结果如图6和图7所示。
步骤S4，建立自适应平方根无迹Kalman相位解缠模型；获得相位解缠起始像元，假设解缠起始像元在干涉图中的位置为k-1，待解缠像元的位置为k。
步骤S4.1，计算解缠起始像元的Sigma采样点及其权值系数，公式为
Figure BDA0003172748920000121
其中ρi为起始像元的Sigma采样点，
Figure BDA0003172748920000122
是第k-1位置解缠相位的状态估计值，Sk-1是第k-1位置像元的误差协方差的平方根，N表示状态变量的维数(N＝1)；{·}i表示括号内矩阵的第i列，η表示尺度函数；
Figure BDA0003172748920000123
其中
Figure BDA0003172748920000124
和
Figure BDA0003172748920000125
为权值系数，，η＝τ2(N+α)-N表示尺度函数，τ，α，ζ用来调节Sigma点，通常τ＝0.01，α＝0，ζ＝2，τ决定Sigma点在均值附近的散布程度，α为冗余量。
步骤S4.2，将解缠起始像元的Sigma点通过自适应平方根无迹Kalman滤波相位解缠模型的状态方程计算传播后的均值和误差协方差的平方根，公式为
Figure BDA0003172748920000126
其中符号qr{·}和cholupdate{·}表示QR分解运算和Cholesky因子一阶更新运算，
Figure BDA0003172748920000127
表示第k个像元Sigma点预测值；
Figure BDA0003172748920000128
表示干涉相位图第k个像元的状态估计；Sψψ表示干涉相位图第k个像元状态预测误差协方差平方根；Q(k-1)是干涉相位图第k个像元的相位梯度估计误差方差。
步骤S5，计算待解缠像元的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差；采用测量噪声等价协方差平方根阵原理对测量误差方差的平方根进行自适应修正。
步骤S5.1，计算待解缠像元的量测预测值，公式为：
Figure BDA0003172748920000131
其中，
Figure BDA0003172748920000132
表示k-1像元的Sigma点经传递后得到的k像元的Sigma点预测值；
Figure BDA0003172748920000133
表示待解缠像元的量测预测值；h[·]为观测方程系数。
步骤S5.2，计算待解缠像元的量测预测值协方差的平方根和状态与量测间的互协方差，公式为
Figure BDA0003172748920000134
Figure BDA0003172748920000135
其中
Figure BDA0003172748920000136
表示待解缠像元k的量测预测值协方差的平方根；
Figure BDA0003172748920000137
表示待解缠像元k的状态与量测间的互协方差；
步骤S5.3，计算待解缠像元的预测残差信息，公式为
Figure BDA0003172748920000138
式中
Figure BDA0003172748920000139
为k位置像元的预测残差，
Figure BDA00031727489200001310
表示k位置像元的真实缠绕相位，
Figure BDA00031727489200001311
为待解缠像元的量测预测值；
步骤S5.4，计算待解缠像元的测量噪声协方差的等价协方差阵的对角线元素，公式为
Figure BDA00031727489200001312
其中Ri表示测量噪声协方差阵R的第i个对角线元素；
Figure BDA00031727489200001313
表示测量噪声协方差的等价协方差阵第i个对角线元素；U0和U1为小数值的常数阈值，这里可设置U0＝0.45,U1＝3；vi为标准化预测残差。
Figure BDA00031727489200001314
表示预测观测值协方差阵的第i个对角线元素；median(·)表示取中位数算子。
步骤S5.5，将公式(19)求得的对角线元素
Figure BDA00031727489200001315
组成的测量噪声协方差的等价协方差阵
Figure BDA00031727489200001316
带入式(16)中替换测量噪声协方差阵R重新更新计算量测预测值的协方差平方根。
步骤S6，计算自适应平方根无迹Kalman滤波相位解缠模型的滤波增益，对状态估计及其协方差平方根进行更新，公式为
Figure BDA0003172748920000141
其中υk表示干涉相位图第k个像元的增益矩阵；
Figure BDA0003172748920000142
表示干涉相位图第k个像元的状态估计值；S(k)表示干涉相位图第k个像元状态估计的协方差平方根。
步骤S7，通过基于堆排序有效质量指导策略对干涉图进行逐个相位展开；重复步骤S4-S6，直到所有像元都被解缠。最终得到干涉图的相位解缠结果。相位解缠的结果如图8所示。
为本发明针对SAR数据高噪声及大梯度变化区域的相位解缠方法的示例性实施例更加有说服力，我们对此同样的干涉图分别采用MCF、SNAPHU、无迹Kalman相位解缠等方法进行解缠处理。MCF方法的解缠结果如图9所示；SNAPHU方法的解缠结果如图10所示；无迹Kalman相位解缠方法的解缠结果如图11所示。从图中可以看出，本专利提出的方法可以获得较其它方法更好的解缠结果；为了定量描述相位解缠的质量，我们将不同方法的解缠结果与由高精度DEM生成的模拟干涉图(图12)进行差分处理，计算各种相位解缠方法的解缠结果的平均误差，其中，MCF方法的解缠结果平均误差为0.8411rad；SNAPHU方法的解缠结果的平均误差为0.7822rad，无迹Kalman相位解缠方法的解缠结果的平均误差为0.6981rad，本发明的一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法的解缠结果的平均误差为0.5948rad。结果表明，本发明的一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法具有更好的相位解缠精度。

Claims (5)

1.一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法，所述方法具体包括以下步骤：
步骤S1，收集需要提取高程量或微小形变量地区的SAR数据，对SAR数据进行预处理得到干涉相位图；通过改进的修正矩阵束模型对干涉相位图所有像元的局部频率进行估计，改进的修正矩阵束是在修正矩阵束模型的奇异值处理步骤中替用一阶Butterworth低通滤波器对分解后得到的每一个奇异值赋予不同的权值来实现信号重构，从而获取准确的局部频率，以此对干涉相位图中所有像元的方位向和距离向的相位梯度进行估计：局部频率估计窗口的大小根据干涉图像元相位均方差自适应选择，窗口包括小的窗口和大的窗口，像元的干涉相位均方差反映条纹的密集程度，条纹密集则均方差大，利用小的窗口进行局部频率估计以免破坏局部条纹平稳的假设；条纹稀疏则均方差小，利用大的窗口进行局部频率估计，既满足局部条纹平稳的假设，又增加了样本点数量，提高了频率估计结果的准确性；
步骤S2，由于局部频率估计结果不存在大的跳变，当局部频率估计结果出现明显的不连续时，则认定该局部像元的频率估计结果存在误差，需要进行相应地修正，利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测；
Figure FDA0003172748910000011
Figure FDA0003172748910000012
式中M为与判定连续性窗口大小相关的参数，例如M＝3，即表示判定连续性窗口大小为(2M+1)×(2M+1)＝7×7；u和v表示局部窗口内像元偏离k位置中心像元的程度，
Figure FDA0003172748910000013
和
Figure FDA0003172748910000014
分别表示干涉相位图中k位置中心像元的频率估计结果在距离向和方位向的可信值，fx和fy分别是距离向和方位向上的频率估计结果；
步骤S3，根据求得的可信值
Figure FDA0003172748910000015
和
Figure FDA0003172748910000016
设置阈值，分别设置为
Figure FDA0003172748910000017
和
Figure FDA0003172748910000018
其中
Figure FDA0003172748910000019
和
Figure FDA00031727489100000110
为
Figure FDA00031727489100000111
的最大值和
Figure FDA00031727489100000112
的最大值，利用下式将干涉相位图中k位置中心像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值：
Figure FDA00031727489100000113
Figure FDA00031727489100000114
其中
Figure FDA00031727489100000115
和
Figure FDA00031727489100000116
表示修正后干涉相位图中k位置中心像元在方位向和距离向的窗口范围内所有像元频率的平均值，从而解决了局部频率估计结果存在异常值的问题；
步骤S4，建立无迹Kalman滤波相位解缠观测方程的自适应平方根模型，并利用无迹Kalman滤波相位解缠观测方程的自适应平方根模型获得相位解缠起始像元；首先计算自适应平方根无迹Kalman相位解缠模型的Sigma采样点及其权值系数；将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传播后均值和误差协方差的平方根，利用误差协方差阵的平方根形式替代协方差阵进行递推计算，从而保证了误差协方差阵的非负定性，提高了数值的稳定性，并且在计算中只需要存储和运算平方根因数，从而降低计算负担；
步骤S5，计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差；其中采用测量噪声等价协方差平方根阵原理，引入自适应因子，实现对测量噪声协方差平方根的异常值分离和自适应修正，进而提高平方根无迹Kalman相位解缠模型的精度；
步骤S6，计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的滤波增益；利用求得的Kalman滤波增益矩阵对状态估计及其协方差平方根进行更新，获取鲁棒性更好的无迹Kalman相位解缠模型；
步骤S7，通过基于堆排序的有效质量指导策略对干涉相位图进行逐个相位展开，将干涉复数信号所有像元的相位质量，即质量图作为指导策略，以最大堆作为数据结构，通过不断更新最大堆的节点来实现对质量图的排序，保证了相位解缠过程一直从高质量像元向低质量像元计算；干涉相位图中每一个像元的质量是由相干系数表示，相干系数位于区间[0，1]，相干系数越大则代表像元的质量越高，相反若相干系数越小则代表像元的质量越低，所有像元的相干系数组合在一起即为质量图；重复步骤S4-步骤S6，直到干涉相位图中所有像元都被解缠。
2.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法，其特征在于使用改进的修正矩阵束模型对干涉图距离向和方位向梯度进行估计：
步骤a1，计算干涉图中所有像元干涉相位的均方差ξ，公式为
Figure FDA0003172748910000021
其中l×l为计算像元干涉相位均方差的一个窗口大小，即判断条纹密集程度的一个窗口大小，一般可设置l＝5，As表示局部窗口内第s个像元的干涉相位复数的振幅；
Figure FDA0003172748910000022
为A的平均值；
步骤a2，使用公式(1)求得的干涉相位均方差ξ，根据ξ的值自适应选择待估像元的局部频率估计窗口大小，其中7×7或9×9为小的窗口，17×17或19×19为大的窗口，公式为：
Figure FDA0003172748910000031
其中，H×H表示干涉相位图中k位置中心像元进行局部频率估计计算的窗口大小；
步骤a3，假设Ik是以k像元为中心，窗口大小为H×H的干涉相位样本，利用公式：Ik＝UDVT对信号矩阵Ik进行奇异值分解，其中U和V表示H阶正交(酉)矩阵；D表示对应的特征值矩阵，D＝diag(σ1,σ2,…σH)，且σ1≥σ2≥…≥σH≥0，σi,i＝1,2,...,H表示信号矩阵Ik奇异值分解后对角线上的特征值；
步骤a4，由一阶Butterworth低通滤波器计算传递函数，即对信号矩阵Ik分解后对角线上的每个特征值赋予不同的权值以降低高频噪声对信号重构的影响，具体计算过程为：
Figure FDA0003172748910000032
式中Wt表示特征值矩阵D对角线上第t个奇异值即将被分配的权值系数，σa表示前t个奇异值中第a个奇异值，σt表示特征值矩阵D对角线上第t个奇异值；
步骤a5，利用公式(4)计算处理后的干涉相位样本矩阵
Figure FDA0003172748910000033
即尽量剔除信号中高频噪声成分，尽量完整地保留有用的相位信息：
Figure FDA0003172748910000034
步骤a6，利用公式(5)根据处理后的干涉相位样本矩阵
Figure FDA0003172748910000035
构建母Hankel矩阵以及两个子Hankel矩阵：
Figure FDA0003172748910000036
其中
Figure FDA0003172748910000037
及
Figure FDA0003172748910000038
分别为母Hankel矩阵以及两个子Hankel矩阵；
步骤a7，利用公式(6)对母Hankel矩阵
Figure FDA0003172748910000039
进行奇异值分解：
Figure FDA00031727489100000310
其中
Figure FDA0003172748910000041
和
Figure FDA0003172748910000042
仅具有与有用信号相关的主特征值对应的信息；
步骤a8，通过奇异值分解，并结合公式(5)和公式(6)得到降维处理后的相位数据信号，由于相位数据矩阵具有低秩结构，经过降维处理后以抑制噪声，推导获得如下公式：
Figure FDA0003172748910000043
其中
Figure FDA0003172748910000044
及
Figure FDA0003172748910000045
分别是通过母Hankel矩阵
Figure FDA0003172748910000046
以及两个子Hankel矩阵
Figure FDA0003172748910000047
和
Figure FDA0003172748910000048
重构后的信号矩阵；
步骤a9，利用公式(8)表示以k像元为中心窗口内样本干涉相位的距离向和方位向的局部频率：
Figure FDA0003172748910000049
其中
Figure FDA00031727489100000410
和
Figure FDA00031727489100000411
分别为
Figure FDA00031727489100000412
和
Figure FDA00031727489100000413
的伪逆矩阵，conj(·)为矩阵的共轭相乘因子；
Figure FDA00031727489100000414
和
Figure FDA00031727489100000415
分别是k位置像元在距离向和方位向上的频率估计结果；
步骤a10，由于局部频率估计结果不存在大的跳变，当频率估计结果出现明显的不连续时，则认定该局部像元的频率估计结果存在误差，需要进行相应地修正，利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测；
Figure FDA00031727489100000416
式中M为与判定连续性窗口大小相关的参数，例如M＝3，即表示判定连续性窗口大小为(2M+1)×(2M+1)＝7×7；u和v表示局部窗口内像元偏离k位置中心像元的程度，
Figure FDA00031727489100000417
和
Figure FDA00031727489100000418
分别表示干涉相位图中k位置中心像元的频率估计结果在方位向和距离向的可信值，fx和fy分别是像元在距离向和方位向上的频率估计结果；
步骤a11，根据求得的可信值
Figure FDA00031727489100000419
和
Figure FDA00031727489100000420
设置阈值，分别设置为
Figure FDA00031727489100000421
和
Figure FDA00031727489100000422
其中
Figure FDA00031727489100000423
和
Figure FDA00031727489100000424
为
Figure FDA00031727489100000425
的最大值和
Figure FDA00031727489100000426
的最大值，利用下式计算该像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值：
Figure FDA0003172748910000051
其中fx和fy分别是像元在距离向和方位向上的频率估计结果；
Figure FDA0003172748910000052
和
Figure FDA0003172748910000053
表示修正后像元k在距离向和方位向的频率估计结果，以解决局部频率估计结果存在异常值的问题；
步骤a12，根据步骤a11获得局部频率可以获得像素k的距离向和方位向梯度，公式如下：
Figure FDA0003172748910000054
其中
Figure FDA0003172748910000055
和
Figure FDA0003172748910000056
分别为像素k的距离向和方位向梯度；
重复步骤a2-步骤a12，直至估计出干涉相位图中所有像元距离向和方位向的梯度。
3.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法，其特征在于建立自适应平方根无迹Kalman滤波相位解缠模型，获得相位解缠起始像元的步骤为：
步骤b1，利用下式计算自适应平方根无迹Kalman相位解缠模型解缠起始像元的Sigma采样点及其权值系数，公式为：
Figure FDA0003172748910000057
其中ρi为起始像元的Sigma采样点，
Figure FDA0003172748910000058
是第k-1位置解缠相位的状态估计值，Sk-1是第k-1位置像元的误差协方差的平方根，N表示状态变量的维数(N＝1)；{·}i表示括号内矩阵的第i列，η表示尺度函数；
Figure FDA0003172748910000059
其中
Figure FDA00031727489100000510
和
Figure FDA00031727489100000511
为权值系数，η＝τ2(N+α)-N表示尺度函数，τ，α，ζ用来调节Sigma点，通常τ＝0.01,ζ＝2,α＝0，τ决定Sigma点在均值附近的散布程度，α为冗余量；
步骤b2，将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传递后均值和误差协方差的平方根，公式为：
Figure FDA0003172748910000061
其中f[·]表示状态方程转移函数；符号qr{·}和cholupdate{·}表示QR分解运算和Cholesky因子一阶更新运算，
Figure FDA0003172748910000062
表示第k个像元Sigma点预测值；
Figure FDA0003172748910000063
表示干涉相位图第k个像元的状态估计；Sψψ表示干涉相位图第k个像元状态预测误差协方差平方根；Q(k-1)是干涉相位图第k个像元的相位梯度估计误差方差。
4.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法，其特征在于计算自适应平方根无迹Kalman滤波相位解缠模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差；采用测量噪声等价协方差平方根阵原理对测量误差方差的平方根进行自适应修正的具体步骤为：
步骤c1，计算待解缠像元的量测预测值，公式为：
Figure FDA0003172748910000064
其中，
Figure FDA0003172748910000065
表示k-1像元的Sigma点经传递后得到的k像元的Sigma点预测值；
Figure FDA0003172748910000066
表示待解缠像元的量测预测值；h[·]为观测方程系数。
步骤c2，利用公式(16)和公式(17)计算待解缠像元的量测预测值协方差的平方根和状态与量测间的互协方差：
Figure FDA0003172748910000067
Figure FDA0003172748910000068
其中
Figure FDA0003172748910000069
表示待解缠像元k的量测预测值协方差的平方根；
Figure FDA00031727489100000610
表示待解缠像元k的状态与量测间的互协方差；
步骤c3，利用下式计算待解缠像元的预测残差信息：
Figure FDA0003172748910000071
式中
Figure FDA0003172748910000072
为k位置像元的预测残差，
Figure FDA0003172748910000073
表示k位置像元的真实缠绕相位，
Figure FDA0003172748910000074
为待解缠像元的量测预测值；
步骤c4，利用下式计算待解缠像元测量噪声协方差的等价协方差阵的对角线元素：
Figure FDA0003172748910000075
其中Ri表示测量噪声协方差阵R的第i个对角线元素；
Figure FDA0003172748910000076
表示测量噪声协方差的等价协方差阵第i个对角线元素；U0和U1为小数值的常数阈值，这里可设置U0＝0.45,U1＝3；vi为标准化预测残差。
Figure FDA0003172748910000077
表示预测观测值协方差阵的第i个对角线元素；median(·)表示取中位数算子。
步骤e5，将公式(19)求得的对角线元素
Figure FDA0003172748910000078
组成的测量噪声协方差的等价协方差阵
Figure FDA0003172748910000079
带入式(16)中替换测量噪声协方差阵R重新更新计算量测预测值的协方差平方根。
5.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法，其特征在于：计算自适应平方根无迹Kalman滤波相位解缠模型的滤波增益，对状态估计及其协方差平方根进行更新，公式为：
Figure FDA00031727489100000710
其中υk表示干涉相位图第k个像元的增益矩阵；
Figure FDA00031727489100000711
表示干涉相位图第k个像元的状态估计值；S(k)表示干涉相位图第k个像元状态估计的协方差平方根。
CN202110823529.3A 2021-07-21 2021-07-21 针对sar数据高噪声及大梯度变化区域的相位解缠方法 Active CN113552565B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110823529.3A CN113552565B (zh) 2021-07-21 2021-07-21 针对sar数据高噪声及大梯度变化区域的相位解缠方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110823529.3A CN113552565B (zh) 2021-07-21 2021-07-21 针对sar数据高噪声及大梯度变化区域的相位解缠方法

Publications (2)

Publication Number Publication Date
CN113552565A true CN113552565A (zh) 2021-10-26
CN113552565B CN113552565B (zh) 2023-07-18

Family

ID=78103706

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110823529.3A Active CN113552565B (zh) 2021-07-21 2021-07-21 针对sar数据高噪声及大梯度变化区域的相位解缠方法

Country Status (1)

Country Link
CN (1) CN113552565B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114265062A (zh) * 2021-11-11 2022-04-01 电子科技大学 一种基于相位梯度估计网络的InSAR相位解缠方法
CN117269960A (zh) * 2023-09-12 2023-12-22 中国矿业大学 一种基于梯度优化的快速范数相位解缠方法
CN117368916A (zh) * 2023-10-12 2024-01-09 云南大学 一种InSAR相位解缠方法、系统、设备及介质

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020112855A1 (en) * 2001-02-21 2002-08-22 Arndt G. Dickey Medium frequency pseudo noise geological radar
US20080240525A1 (en) * 2007-03-29 2008-10-02 Martti Kalke Method and system for reconstructing a medical image of an object
CN104936531A (zh) * 2013-01-22 2015-09-23 株式会社东芝 超声波诊断装置、图像处理装置以及图像处理方法
CN109725333A (zh) * 2018-12-17 2019-05-07 中国空间技术研究院 一种场景自适应的卫星信号接收及处理方法
CN111983654A (zh) * 2020-08-24 2020-11-24 中国矿业大学 一种基于gnss的北极区域电离层相位闪烁因子构建方法
US20210109210A1 (en) * 2019-10-14 2021-04-15 The Boeing Company Stripmap synthetic aperture radar (sar) system utilizing direct matching and registration in range profile space

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20020112855A1 (en) * 2001-02-21 2002-08-22 Arndt G. Dickey Medium frequency pseudo noise geological radar
US20080240525A1 (en) * 2007-03-29 2008-10-02 Martti Kalke Method and system for reconstructing a medical image of an object
CN104936531A (zh) * 2013-01-22 2015-09-23 株式会社东芝 超声波诊断装置、图像处理装置以及图像处理方法
CN109725333A (zh) * 2018-12-17 2019-05-07 中国空间技术研究院 一种场景自适应的卫星信号接收及处理方法
US20210109210A1 (en) * 2019-10-14 2021-04-15 The Boeing Company Stripmap synthetic aperture radar (sar) system utilizing direct matching and registration in range profile space
CN111983654A (zh) * 2020-08-24 2020-11-24 中国矿业大学 一种基于gnss的北极区域电离层相位闪烁因子构建方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114265062A (zh) * 2021-11-11 2022-04-01 电子科技大学 一种基于相位梯度估计网络的InSAR相位解缠方法
CN114265062B (zh) * 2021-11-11 2023-11-10 电子科技大学 一种基于相位梯度估计网络的InSAR相位解缠方法
CN117269960A (zh) * 2023-09-12 2023-12-22 中国矿业大学 一种基于梯度优化的快速范数相位解缠方法
CN117269960B (zh) * 2023-09-12 2024-04-26 中国矿业大学 一种基于梯度优化的快速范数相位解缠方法
CN117368916A (zh) * 2023-10-12 2024-01-09 云南大学 一种InSAR相位解缠方法、系统、设备及介质

Also Published As

Publication number Publication date
CN113552565B (zh) 2023-07-18

Similar Documents

Publication Publication Date Title
CN113552565A (zh) 针对sar数据高噪声及大梯度变化区域的相位解缠方法
CN109782282B (zh) 一种集成对流层大气延迟改正的时间序列InSAR分析方法
CN111059998B (zh) 一种基于高分辨率的时序InSAR形变监测方法及系统
CN111273293B (zh) 一种顾及地形起伏的InSAR残余运动误差估计方法及装置
CN102270341B (zh) 一种自适应的高精度干涉sar相位估计方法
CN115856889B (zh) 一种误差自动校正的InSAR时序变形监测方法
CN112363161B (zh) 基于散射机制分解的植被垂直结构及林下地形反演方法及装置
CN109100720B (zh) 一种InSAR地表形变监测方法
CN113589286B (zh) 基于D-LinkNet的无迹卡尔曼滤波相位解缠方法
WO2020240720A1 (ja) 合成開口レーダの信号処理装置および信号処理方法
Gruszczynski et al. A filtering of incomplete GNSS position time series with probabilistic Principal Component Analysis
US5200754A (en) Fourth-order-product phase difference autofocus
Li et al. Least squares-based filter for remote sensingimage noise reduction
CN102289715A (zh) 基于前向线性预测的自适应小波神经网络去噪建模方法
Ichikawa et al. Singular unit restoration in InSAR using complex-valued neural networks in the spectral domain
CN116299247B (zh) 一种基于稀疏卷积神经网络的InSAR大气校正方法
CN113030964A (zh) 基于复Laplace先验的双基地ISAR稀孔径高分辨成像方法
CN111915570A (zh) 基于反向传播神经网络的大气延迟估计方法
CN112034457B (zh) 基于干涉条纹方向的多基线高程干涉相位估计方法
WO2000005598A9 (en) Averaging-area-constrained adaptive interferometric filter that optimizes combined coherent and noncoherent averaging
CN115221768A (zh) 基于多嵌套阵列的改进无迹粒子滤波直接跟踪方法
CN112986990B (zh) 一种大气相位改正方法及系统
CN114996653A (zh) 一种基于原子范数最小化的二维鲁棒自适应波束形成方法
CN110858309B (zh) 一种多基准时钟加权合成方法
CN114966687A (zh) 基于低秩和非局部自相似的稀疏isar成像方法及系统

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