CN113552565B - 针对sar数据高噪声及大梯度变化区域的相位解缠方法 - Google Patents
针对sar数据高噪声及大梯度变化区域的相位解缠方法 Download PDFInfo
- Publication number
- CN113552565B CN113552565B CN202110823529.3A CN202110823529A CN113552565B CN 113552565 B CN113552565 B CN 113552565B CN 202110823529 A CN202110823529 A CN 202110823529A CN 113552565 B CN113552565 B CN 113552565B
- Authority
- CN
- China
- Prior art keywords
- phase
- pixel
- covariance
- matrix
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 82
- 230000008859 change Effects 0.000 title claims abstract description 14
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 21
- 238000001914 filtration Methods 0.000 claims abstract description 19
- 238000012545 processing Methods 0.000 claims abstract description 14
- 230000002159 abnormal effect Effects 0.000 claims abstract description 9
- 239000011159 matrix material Substances 0.000 claims description 93
- 238000005259 measurement Methods 0.000 claims description 51
- 230000003044 adaptive effect Effects 0.000 claims description 18
- 238000012937 correction Methods 0.000 claims description 17
- 238000010587 phase diagram Methods 0.000 claims description 16
- 238000004364 calculation method Methods 0.000 claims description 14
- 238000010586 diagram Methods 0.000 claims description 11
- 238000005070 sampling Methods 0.000 claims description 8
- 230000008569 process Effects 0.000 claims description 7
- 230000009467 reduction Effects 0.000 claims description 6
- 238000012546 transfer Methods 0.000 claims description 5
- 238000012163 sequencing technique Methods 0.000 claims description 4
- 230000007480 spreading Effects 0.000 claims description 3
- 238000004804 winding Methods 0.000 claims description 3
- 230000005540 biological transmission Effects 0.000 claims description 2
- 238000007781 pre-processing Methods 0.000 claims description 2
- 238000000926 separation method Methods 0.000 claims description 2
- 238000005516 engineering process Methods 0.000 abstract description 4
- 230000009286 beneficial effect Effects 0.000 abstract description 3
- 230000000694 effects Effects 0.000 abstract description 3
- 230000004044 response Effects 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 239000000470 constituent Substances 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000005305 interferometry Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000012876 topography Methods 0.000 description 1
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数据高噪声及大梯度变化区域的相位解缠方法。
背景技术
合成孔径雷达干涉(Interferometric SyntheticAperture Radar,InSAR)测量技术近年来已广泛应用于地形三维重建和由地震活动的物理模型反演,尽管目前已经存在几十种相位解缠方法,但是在噪声量大的区域相位解缠方法的鲁棒性一般较差,在地形梯度大的区域相位解缠结果误差一般较大。这几十种相位解缠方法可以大致分为两类:一类是单基线相位解缠方法,另一类是多基线相位解缠方法。由于单基线相位解缠方法受相位连续性假设的限制,在地形梯度大的区域难以获得理想的结果,虽然多基线相位解缠方法可以摆脱相位连续性假设的限制,但是在解缠效率和基线配比等方面仍然存在一些困难。
由于修正矩阵束模型具有一定的频率估计精度和效率,所以本发明保留了修正矩阵束模型的框架,在其奇异值处理过程中引入了低通滤波技术,并且可以根据干涉图像素均方差值自适应选择局部频率估计窗口的大小,尽可能地实现方位向和距离向相位梯度的准确估计。这种低通滤波技术是由Butterworth滤波器来实现的,其特点是通频带内的频率响应曲线最大限度平坦,没有纹波,而在阻频带则逐渐下降为零,这意味着Butterworth低通滤波器可以减少高频噪声。此外,我们也引入了自适应平方根无迹Kalman技术到相位解缠问题中,平方根无迹Kalman算法只需要误差协方差的平方根,而不需要完整的协方差阵,其采用误差协方差阵的平方根形式替代协方差阵进行递推计算,保证了误差协方差阵的非负定性,从而提高了数值的稳定性,并且在计算中只需要存储和运算平方根因数,这样可以降低计算负担,以实现对无迹Kalman相位解缠技术性能的提高。
发明内容
针对现有技术的不足之处,提供一种步骤简单,使用效果好的一种针对SAR数据高噪声及大梯度变化区域的相位解缠方法
为实现上述技术目的,本发明的针对高噪声及大梯度变化区域的相位解缠方法,所包括以下步骤:
步骤S1,收集需要提取高程量或微小形变量地区的SAR数据,对SAR数据进行预处理得到干涉相位图;通过改进的修正矩阵束模型对干涉相位图所有像元的局部频率进行估计,改进的修正矩阵束是在修正矩阵束模型的奇异值处理步骤中替用一阶Butterworth低通滤波器对分解后得到的每一个奇异值赋予不同的权值来实现信号重构,从而获取准确的局部频率,以此对干涉相位图中所有像元的方位向和距离向的相位梯度进行估计:局部频率估计窗口的大小根据干涉图像元相位均方差自适应选择,窗口包括小的窗口和大的窗口,像元的干涉相位均方差反映条纹的密集程度,条纹密集则均方差大,利用小的窗口进行局部频率估计以免破坏局部条纹平稳的假设;条纹稀疏则均方差小,利用大的窗口进行局部频率估计,既满足局部条纹平稳的假设,又增加了样本点数量,提高了频率估计结果的准确性;
步骤S2,由于局部频率估计结果不存在大的跳变,当局部频率估计结果出现明显的不连续时,则认定该局部像元的频率估计结果存在误差,需要进行相应地修正,利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测;
式中M为与判定连续性窗口大小相关的参数,例如M=3,即表示判定连续性窗口大小为(2M+1)×(2M+1)=7×7;u和v表示局部窗口内像元偏离k位置中心像元的程度,和分别表示干涉相位图中k位置中心像元的频率估计结果在距离向和方位向的可信值,fx和fy分别是距离向和方位向上的频率估计结果;
步骤S3,根据求得的可信值和设置阈值,分别设置为和其中和为的最大值和的最大值,利用下式将干涉相位图中k位置中心像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值:
其中和表示修正后干涉相位图中k位置中心像元在方位向和距离向的窗口范围内所有像元频率的平均值,从而解决了局部频率估计结果存在异常值的问题;
步骤S4,建立无迹Kalman滤波相位解缠观测方程的自适应平方根模型,并利用无迹Kalman滤波相位解缠观测方程的自适应平方根模型获得相位解缠起始像元;首先计算自适应平方根无迹Kalman相位解缠模型的Sigma采样点及其权值系数;将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传播后均值和误差协方差的平方根,利用误差协方差阵的平方根形式替代协方差阵进行递推计算,从而保证了误差协方差阵的非负定性,提高了数值的稳定性,并且在计算中只需要存储和运算平方根因数,从而降低计算负担;
步骤S5,计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差;其中采用测量噪声等价协方差平方根阵原理,引入自适应因子,实现对测量噪声协方差平方根的异常值分离和自适应修正,进而提高平方根无迹Kalman相位解缠模型的精度;
步骤S6,计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的滤波增益;利用求得的Kalman滤波增益矩阵对状态估计及其协方差平方根进行更新,获取鲁棒性更好的无迹Kalman相位解缠模型;
步骤S7,通过基于堆排序的有效质量指导策略对干涉相位图进行逐个相位展开,将干涉复数信号所有像元的相位质量,即质量图作为指导策略,以最大堆作为数据结构,通过不断更新最大堆的节点来实现对质量图的排序,保证了相位解缠过程一直从高质量像元向低质量像元计算;干涉相位图中每一个像元的质量是由相干系数表示,相干系数位于区间[0,1],相干系数越大则代表像元的质量越高,相反若相干系数越小则代表像元的质量越低,所有像元的相干系数组合在一起即为质量图;重复步骤S4-步骤S6,直到干涉相位图中所有像元都被解缠。
进一步,使用改进的修正矩阵束模型对干涉图距离向和方位向梯度进行估计:
步骤a1,计算干涉图中所有像元干涉相位的均方差ξ,公式为
其中l×l为计算像元干涉相位均方差的一个窗口大小,即判断条纹密集程度的一个窗口大小,一般可设置l=5,As表示局部窗口内第s个像元的干涉相位复数的振幅;为A的平均值;
步骤a2,使用公式(1)求得的干涉相位均方差ξ,根据ξ值自适应选择待估像元的局部频率估计窗口大小,其中7×7或9×9为小的窗口,17×17或19×19为大的窗口,公式为:
其中,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分解后对角线上的每个特征值赋予不同的权值以降低高频噪声对信号重构的影响,具体计算过程为:
式中Wt表示特征值矩阵D对角线上第t个奇异值即将被分配的权值系数,σa表示前t个奇异值中第a个奇异值,σt表示特征值矩阵D对角线上第t个奇异值;
步骤a5,利用公式(4)计算处理后的干涉相位样本矩阵即尽量剔除信号中高频噪声成分,尽量完整地保留有用的相位信息:
步骤a6,利用公式(5)根据处理后的干涉相位样本矩阵构建母Hankel矩阵以及两个子Hankel矩阵:
其中及分别为母Hankel矩阵以及两个子Hankel矩阵;
步骤a7,利用公式(6)对母Hankel矩阵进行奇异值分解:
其中和仅具有与有用信号相关的主特征值对应的信息;
步骤a8,通过奇异值分解,并结合公式(5)和公式(6)得到降维处理后的相位数据信号,由于相位数据矩阵具有低秩结构,经过降维处理后以抑制噪声,推导获得如下公式:
其中及分别是通过母Hankel矩阵以及两个子Hankel矩阵和重构后的信号矩阵;
步骤a9,利用公式(8)表示以k像元为中心窗口内样本干涉相位的距离向和方位向的局部频率:
其中和分别为和的伪逆矩阵,conj(·)为矩阵的共轭相乘因子;和分别是k位置像元在距离向和方位向上的频率估计结果;
步骤a10,由于局部频率估计结果不存在大的跳变,当频率估计结果出现明显的不连续时,则认定该局部像元的频率估计结果存在误差,需要进行相应地修正,利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测;
式中M为与判定连续性窗口大小相关的参数,例如M=3,即表示判定连续性窗口大小为(2M+1)×(2M+1)=7×7;u和v表示局部窗口内像元偏离k位置中心像元的程度,和分别表示干涉相位图中k位置中心像元的频率估计结果在方位向和距离向的可信值,fx和fy分别是像元在距离向和方位向上的频率估计结果;
步骤a11,根据求得的可信值和设置阈值,分别设置为和其中和为的最大值和的最大值,利用下式计算该像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值:
其中fx和fy分别是像元在距离向和方位向上的频率估计结果;和表示修正后像元k在距离向和方位向的频率估计结果,以解决局部频率估计结果存在异常值的问题;
步骤a12,根据步骤a11获得局部频率可以获得像素k的距离向和方位向梯度,公式如下:
其中和分别为像素k的距离向和方位向梯度;
重复步骤a2-步骤a12,直至估计出干涉相位图中所有像元距离向和方位向的梯度。
进一步,建立自适应平方根无迹Kalman滤波相位解缠模型,获得相位解缠起始像元的步骤为:
步骤b1,利用下式计算自适应平方根无迹Kalman相位解缠模型解缠起始像元的Sigma采样点及其权值系数,公式为:
其中ρi为起始像元的Sigma采样点,是第k-1位置解缠相位的状态估计值,Sk-1是第k-1位置像元的误差协方差的平方根,N表示状态变量的维数(N=1);{·}i表示括号内矩阵的第i列,η表示尺度函数;
其中和为权值系数,η=τ2(N+α)-N表示尺度函数,τ,α,ζ用来调节Sigma点,通常τ=0.01,ζ=2,α=0,τ决定Sigma点在均值附近的散布程度,α为冗余量;
步骤b2,将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传递后均值和误差协方差的平方根,公式为:
其中f[·]表示状态方程转移函数;符号qr{·}和cholupdate{·}表示QR分解运算和Cholesky因子一阶更新运算,表示第k个像元Sigma点预测值;表示干涉相位图第k个像元的状态估计;Sψψ表示干涉相位图第k个像元状态预测误差协方差平方根;Q(k-1)是干涉相位图第k个像元的相位梯度估计误差方差。
进一步,计算自适应平方根无迹Kalman滤波相位解缠模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差;采用测量噪声等价协方差平方根阵原理对测量误差方差的平方根进行自适应修正的具体步骤为:
步骤c1,计算待解缠像元的量测预测值,公式为:
其中,表示k-1像元的Sigma点经传递后得到的k像元的Sigma点预测值;表示待解缠像元的量测预测值;h[·]为观测方程系数。
步骤c2,利用公式(16)和公式(17)计算待解缠像元的量测预测值协方差的平方根和状态与量测间的互协方差:
其中表示待解缠像元k的量测预测值协方差的平方根;表示待解缠像元k的状态与量测间的互协方差;
步骤c3,利用下式计算待解缠像元的预测残差信息:
式中为k位置像元的预测残差,表示k位置像元的真实缠绕相位,为待解缠像元的量测预测值;
步骤c4,利用下式计算待解缠像元测量噪声协方差的等价协方差阵的对角线元素:
其中Ri表示测量噪声协方差阵R的第i个对角线元素;表示测量噪声协方差的等价协方差阵第i个对角线元素;U0和U1为小数值的常数阈值,这里可设置U0=0.45,U1=3;vi为标准化预测残差。表示预测观测值协方差阵的第i个对角线元素;median(·)表示取中位数算子。
步骤c5,将公式(19)求得的对角线元素组成的测量噪声协方差的等价协方差阵带入式(16)中替换测量噪声协方差阵R重新更新计算量测预测值的协方差平方根。
进一步,计算自适应平方根无迹Kalman滤波相位解缠模型的滤波增益,对状态估计及其协方差平方根进行更新,公式为
其中υk表示干涉相位图第k个像元的增益矩阵;表示干涉相位图第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,计算干涉图中所有像元干涉相位的均方差ξ,公式为:
其中l×l为计算像元干涉相位均方差的一个窗口大小,即判断条纹密集程度的一个窗口大小,一般可设置l=5;As表示局部窗口内第s个像元的干涉相位复数的振幅;为A的平均值;
步骤S1.2,使用公式(1)求得的干涉相位均方差ξ,根据ξ值自适应选择待估像元的局部频率估计窗口大小,其中7×7或9×9为小的窗口,17×17或19×19为大的窗口,公式为:
其中,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分解后对角线上的每个特征值赋予不同的权值以降低高频噪声对信号重构的影响,具体计算过程为:
式中Wt表示特征值矩阵D对角线上第t个奇异值即将被分配的权值系数,σa表示前t个奇异值中第a个奇异值,σt表示特征值矩阵D对角线上第t个奇异值;
步骤S1.5,利用公式(4)计算处理后的干涉相位样本矩阵即尽量剔除信号中高频噪声成分,尽量完整地保留有用的相位信息:
步骤S1.6,利用公式(5)根据处理后的干涉相位样本矩阵构建母Hankel矩阵以及两个子Hankel矩阵:
其中及分别为母Hankel矩阵以及两个子Hankel矩阵;
步骤S1.7,利用公式(6)对母Hankel矩阵进行奇异值分解:
其中和仅具有与有用信号相关的主特征值对应得信息;
步骤S1.8,通过奇异值分解,结合公式(5)和公式(6)可以得到降维处理后的相位数据信号,由于相位数据矩阵具有低秩结构,经过降维处理后可以抑制噪声,推导获得如下公式:
其中及分别是通过母Hankel矩阵以及两个子Hankel矩阵和重构后的信号矩阵;
步骤S1.9,利用公式(8)表示以k像元为中心窗口内样本干涉相位的距离向和方位向的局部频率:
其中和分别为和的伪逆矩阵,conj(·)为矩阵的共轭相乘因子;和分别是k位置像元在距离向和方位向上的频率估计结果;
步骤S1.10,根据步骤S1.9获得局部频率可以获得像素k的距离向和方位向梯度,公式如下:
其中和分别为像素k的距离向和方位向梯度;
重复步骤S1.2-S1.10,直至估计出所有像元距离向和方位向的梯度。方位向和距离向的相位梯度估计结果如图4和图5所示。
步骤S2,根据公式(9),从距离向和方位向上分别对频率估计结果的连续性进行计算和检测。阈值分别设置为和其中和为的最大值和的最大值。当可信值超过阈值时,认为该像元的频率估计结果存在异常,需要进行修正。
式中M为与判定连续性窗口大小相关的参数,例如M=3,即表示判定连续性窗口大小为(2M+1)×(2M+1)=7×7;u和v表示局部窗口内像元偏离k位置中心像元的程度,和分别表示干涉相位图中k位置中心像元的频率估计结果在方位向和距离向的可信值,fx和fy分别是像元在距离向和方位向上的频率估计结果;
步骤S3,从距离向和方位向上分别对频率估计结果的异常值进行修正,异常像元的频率估计结果优化为以k位置像元为中心(2M+1)×(2M+1)范围内所有像元频率的平均,公式为
其中fx和fy分别是像元在距离向和方位向上的频率估计结果;和表示修正后像元在方位向和距离向的频率估计结果。由公式(10)的频率结果重新根据公式(28)计算修正后的方位向和距离向相位梯度。本发明实例干涉图在方位向和距离向上的相位梯度修正后的结果如图6和图7所示。
步骤S4,建立自适应平方根无迹Kalman相位解缠模型;获得相位解缠起始像元,假设解缠起始像元在干涉图中的位置为k-1,待解缠像元的位置为k。
步骤S4.1,计算解缠起始像元的Sigma采样点及其权值系数,公式为
其中ρi为起始像元的Sigma采样点,是第k-1位置解缠相位的状态估计值,Sk-1是第k-1位置像元的误差协方差的平方根,N表示状态变量的维数(N=1);{·}i表示括号内矩阵的第i列,η表示尺度函数;
其中和为权值系数,,η=τ2(N+α)-N表示尺度函数,τ,α,ζ用来调节Sigma点,通常τ=0.01,α=0,ζ=2,τ决定Sigma点在均值附近的散布程度,α为冗余量。
步骤S4.2,将解缠起始像元的Sigma点通过自适应平方根无迹Kalman滤波相位解缠模型的状态方程计算传播后的均值和误差协方差的平方根,公式为
其中符号qr{·}和cholupdate{·}表示QR分解运算和Cholesky因子一阶更新运算,表示第k个像元Sigma点预测值;表示干涉相位图第k个像元的状态估计;Sψψ表示干涉相位图第k个像元状态预测误差协方差平方根;Q(k-1)是干涉相位图第k个像元的相位梯度估计误差方差。
步骤S5,计算待解缠像元的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差;采用测量噪声等价协方差平方根阵原理对测量误差方差的平方根进行自适应修正。
步骤S5.1,计算待解缠像元的量测预测值,公式为:
其中,表示k-1像元的Sigma点经传递后得到的k像元的Sigma点预测值;表示待解缠像元的量测预测值;h[·]为观测方程系数。
步骤S5.2,计算待解缠像元的量测预测值协方差的平方根和状态与量测间的互协方差,公式为
其中表示待解缠像元k的量测预测值协方差的平方根;表示待解缠像元k的状态与量测间的互协方差;
步骤S5.3,计算待解缠像元的预测残差信息,公式为
式中为k位置像元的预测残差,表示k位置像元的真实缠绕相位,为待解缠像元的量测预测值;
步骤S5.4,计算待解缠像元的测量噪声协方差的等价协方差阵的对角线元素,公式为
其中Ri表示测量噪声协方差阵R的第i个对角线元素;表示测量噪声协方差的等价协方差阵第i个对角线元素;U0和U1为小数值的常数阈值,这里可设置U0=0.45,U1=3;vi为标准化预测残差。表示预测观测值协方差阵的第i个对角线元素;median(·)表示取中位数算子。
步骤S5.5,将公式(19)求得的对角线元素组成的测量噪声协方差的等价协方差阵带入式(16)中替换测量噪声协方差阵R重新更新计算量测预测值的协方差平方根。
步骤S6,计算自适应平方根无迹Kalman滤波相位解缠模型的滤波增益,对状态估计及其协方差平方根进行更新,公式为
其中υk表示干涉相位图第k个像元的增益矩阵;表示干涉相位图第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,由于局部频率估计结果不存在大的跳变,当局部频率估计结果出现明显的不连续时,则认定该像元的频率估计结果存在误差,需要进行相应地修正,利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测;
式中M为与判定连续性窗口大小相关的参数,M=3,即表示判定连续性窗口大小为(2M+1)×(2M+1)=7×7;u和v表示局部窗口内像元偏离k位置中心像元的程度,和分别表示干涉相位图中k位置中心像元的频率估计结果在距离向和方位向的可信值,fx和fy分别是距离向和方位向上的频率估计结果;
步骤S3,根据求得的可信值和设置阈值,分别设置为和其中和为的最大值和的最大值,利用下式将干涉相位图中k位置中心像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值:
其中和表示修正后干涉相位图中k位置中心像元在方位向和距离向的窗口范围内所有像元频率的平均值,从而解决了局部频率估计结果存在异常值的问题;
步骤S4,建立无迹Kalman滤波相位解缠观测方程的自适应平方根模型,并利用无迹Kalman滤波相位解缠观测方程的自适应平方根模型获得相位解缠起始像元;首先计算自适应平方根无迹Kalman相位解缠模型的Sigma采样点及其权值系数;将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传播后均值和误差协方差的平方根,利用误差协方差阵的平方根形式替代协方差阵进行递推计算,从而保证了误差协方差阵的非负定性,提高了数值的稳定性,并且在计算中只需要存储和运算平方根因数,从而降低计算负担;
步骤S5,计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差;其中采用测量噪声等价协方差平方根阵原理,引入自适应因子,实现对测量噪声协方差平方根的异常值分离和自适应修正,进而提高平方根无迹Kalman相位解缠模型的精度;
步骤S6,计算无迹Kalman滤波相位解缠观测方程的自适应平方根模型的滤波增益;利用求得的Kalman滤波增益矩阵对状态估计及其协方差平方根进行更新,获取鲁棒性更好的无迹Kalman相位解缠模型;
步骤S7,通过基于堆排序的有效质量指导策略对干涉相位图进行逐个相位展开,将干涉复数信号所有像元的相位质量,即质量图作为指导策略,以最大堆作为数据结构,通过不断更新最大堆的节点来实现对质量图的排序,保证了相位解缠过程一直从高质量像元向低质量像元计算;干涉相位图中每一个像元的质量是由相干系数表示,相干系数位于区间[0,1],相干系数越大则代表像元的质量越高,相反若相干系数越小则代表像元的质量越低,所有像元的相干系数组合在一起即为质量图;重复步骤S4-步骤S6,直到干涉相位图中所有像元都被解缠。
2.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法,其特征在于使用改进的修正矩阵束模型对干涉图距离向和方位向梯度进行估计:
步骤a1,计算干涉图中所有像元干涉相位的均方差ξ,公式为
其中l×l为计算像元干涉相位均方差的一个窗口大小,即判断条纹密集程度的一个窗口大小,设置l=5,As表示局部窗口内第s个像元的干涉相位复数的振幅;为A的平均值;
步骤a2,使用公式(1)求得的干涉相位均方差ξ,根据ξ的值自适应选择待估像元的局部频率估计窗口大小,其中7×7或9×9为小的窗口,17×17或19×19为大的窗口,公式为:
其中,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分解后对角线上的每个特征值赋予不同的权值以降低高频噪声对信号重构的影响,具体计算过程为:
式中Wt表示特征值矩阵D对角线上第t个奇异值即将被分配的权值系数,σa表示前t个奇异值中第a个奇异值,σt表示特征值矩阵D对角线上第t个奇异值;
步骤a5,利用公式(4)计算处理后的干涉相位样本矩阵即尽量剔除信号中高频噪声成分,尽量完整地保留有用的相位信息:
步骤a6,利用公式(5)根据处理后的干涉相位样本矩阵构建母Hankel矩阵以及两个子Hankel矩阵:
其中及分别为母Hankel矩阵以及两个子Hankel矩阵;
步骤a7,利用公式(6)对母Hankel矩阵进行奇异值分解:
其中和仅具有与有用信号相关的主特征值对应的信息;
步骤a8,通过奇异值分解,并结合公式(5)和公式(6)得到降维处理后的相位数据信号,由于相位数据矩阵具有低秩结构,经过降维处理后以抑制噪声,推导获得如下公式:
其中及分别是通过母Hankel矩阵以及两个子Hankel矩阵和重构后的信号矩阵;
步骤a9,利用公式(8)表示以k像元为中心窗口内样本干涉相位的距离向和方位向的局部频率:
其中和分别为和的伪逆矩阵,conj(·)为矩阵的共轭相乘因子;和分别是k位置像元在距离向和方位向上的频率估计结果;
步骤a10,由于局部频率估计结果不存在大的跳变,当频率估计结果出现明显的不连续时,则认定该像元的频率估计结果存在误差,需要进行相应地修正,利用下式从距离向和方位向上分别对频率估计结果的连续性进行计算和检测;
式中M为与判定连续性窗口大小相关的参数,M=3,即表示判定连续性窗口大小为(2M+1)×(2M+1)=7×7;u和v表示局部窗口内像元偏离k位置中心像元的程度,和分别表示干涉相位图中k位置中心像元的频率估计结果在方位向和距离向的可信值,fx和fy分别是像元在距离向和方位向上的频率估计结果;
步骤a11,根据求得的可信值和设置阈值,分别设置为和其中和为的最大值和的最大值,利用下式计算该像元的频率估计结果优化为以k位置像元为中心的窗口(2M+1)×(2M+1)范围内所有像元频率的平均值:
其中fx和fy分别是像元在距离向和方位向上的频率估计结果;和表示修正后像元k在距离向和方位向的频率估计结果,以解决局部频率估计结果存在异常值的问题;
步骤a12,根据步骤a11获得局部频率可以获得像素k的距离向和方位向梯度,公式如下:
其中和分别为像素k的距离向和方位向梯度;
重复步骤a2-步骤a12,直至估计出干涉相位图中所有像元距离向和方位向的梯度。
3.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法,其特征在于建立自适应平方根无迹Kalman滤波相位解缠模型,获得相位解缠起始像元的步骤为:
步骤b1,利用下式计算自适应平方根无迹Kalman相位解缠模型解缠起始像元的Sigma采样点及其权值系数,公式为:
其中ρi为起始像元的Sigma采样点,是第k-1位置解缠相位的状态估计值,Sk-1是第k-1位置像元的误差协方差的平方根,N表示状态变量的维数(N=1);{·}i表示括号内矩阵的第i列,η表示尺度函数;
其中和为权值系数,η=τ2(N+α)-N表示尺度函数,τ,α,ζ用来调节Sigma点,通常τ=0.01,ζ=2,α=0,τ决定Sigma点在均值附近的散布程度,α为冗余量;
步骤b2,将Sigma点通过无迹Kalman滤波相位解缠状态方程计算传递后均值和误差协方差的平方根,公式为:
其中f[·]表示状态方程转移函数;符号qr{·}和cholupdate{·}表示QR分解运算和Cholesky因子一阶更新运算,表示第k个像元Sigma点预测值;表示干涉相位图第k个像元的状态估计;Sψψ表示干涉相位图第k个像元状态预测误差协方差平方根;Q(k-1)是干涉相位图第k个像元的相位梯度估计误差方差。
4.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法,其特征在于计算自适应平方根无迹Kalman滤波相位解缠模型的量测预测值、量测预测值协方差平方根和状态与量测间的互协方差;采用测量噪声等价协方差平方根阵原理对测量误差方差的平方根进行自适应修正的具体步骤为:
步骤c1,计算待解缠像元的量测预测值,公式为:
其中,表示k-1像元的Sigma点经传递后得到的k像元的Sigma点预测值;表示待解缠像元的量测预测值;h[·]为观测方程系数;
步骤c2,利用公式(16)和公式(17)计算待解缠像元的量测预测值协方差的平方根和状态与量测间的互协方差:
其中表示待解缠像元k的量测预测值协方差的平方根;表示待解缠像元k的状态与量测间的互协方差;
步骤c3,利用下式计算待解缠像元的预测残差信息:
式中为k位置像元的预测残差,表示k位置像元的真实缠绕相位,为待解缠像元的量测预测值;
步骤c4,利用下式计算待解缠像元测量噪声协方差的等价协方差阵的对角线元素:
其中Ri表示测量噪声协方差阵R的第i个对角线元素;表示测量噪声协方差的等价协方差阵第i个对角线元素;U0和U1为小数值的常数阈值,这里可设置U0=0.45,U1=3;vi为标准化预测残差;表示预测观测值协方差阵的第i个对角线元素;median(·)表示取中位数算子;
步骤e5,将公式(19)求得的对角线元素组成的测量噪声协方差的等价协方差阵带入式(16)中替换测量噪声协方差阵R重新更新计算量测预测值的协方差平方根。
5.根据权利要求1所述针对SAR数据高噪声及大梯度变化区域的相位解缠方法,其特征在于:计算自适应平方根无迹Kalman滤波相位解缠模型的滤波增益,对状态估计及其协方差平方根进行更新,公式为:
其中υk表示干涉相位图第k个像元的增益矩阵;表示干涉相位图第k个像元的状态估计值;S(k)表示干涉相位图第k个像元状态估计的协方差平方根。
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 CN113552565A (zh) | 2021-10-26 |
CN113552565B true 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) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114265062B (zh) * | 2021-11-11 | 2023-11-10 | 电子科技大学 | 一种基于相位梯度估计网络的InSAR相位解缠方法 |
CN117269960B (zh) * | 2023-09-12 | 2024-04-26 | 中国矿业大学 | 一种基于梯度优化的快速范数相位解缠方法 |
CN117368916B (zh) * | 2023-10-12 | 2024-06-14 | 云南大学 | 一种InSAR相位解缠方法、系统、设备及介质 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
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的北极区域电离层相位闪烁因子构建方法 |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US6651739B2 (en) * | 2001-02-21 | 2003-11-25 | The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration | Medium frequency pseudo noise geological radar |
US8335358B2 (en) * | 2007-03-29 | 2012-12-18 | Palodex Group Oy | Method and system for reconstructing a medical image of an object |
US11333753B2 (en) * | 2019-10-14 | 2022-05-17 | The Boeing Company | Stripmap synthetic aperture radar (SAR) system utilizing direct matching and registration in range profile space |
-
2021
- 2021-07-21 CN CN202110823529.3A patent/CN113552565B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
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的北极区域电离层相位闪烁因子构建方法 |
Also Published As
Publication number | Publication date |
---|---|
CN113552565A (zh) | 2021-10-26 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113552565B (zh) | 针对sar数据高噪声及大梯度变化区域的相位解缠方法 | |
CN109782282B (zh) | 一种集成对流层大气延迟改正的时间序列InSAR分析方法 | |
CN104459696B (zh) | 一种基于平地相位的sar干涉基线精确估计方法 | |
CN105974414B (zh) | 基于二维自聚焦的高分辨聚束sar自聚焦成像方法 | |
CN115856889B (zh) | 一种误差自动校正的InSAR时序变形监测方法 | |
CN111059998B (zh) | 一种基于高分辨率的时序InSAR形变监测方法及系统 | |
CN104123464B (zh) | 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法 | |
CN111273293B (zh) | 一种顾及地形起伏的InSAR残余运动误差估计方法及装置 | |
CN108007476B (zh) | 一种天基干涉成像雷达高度计的干涉定标方法及系统 | |
CN110018476B (zh) | 时间差分基线集时序干涉sar处理方法 | |
CN109100720B (zh) | 一种InSAR地表形变监测方法 | |
CN111025294B (zh) | 基于均方容积卡尔曼滤波器的InSAR相位解缠绕方法 | |
CN109061641A (zh) | 一种基于序贯平差的InSAR时序地表形变监测方法 | |
CN103018729B (zh) | 金属圆柱定标体雷达散射截面的计算方法 | |
CN112269192B (zh) | 一种快速自适应的动态北斗监测实时解算去噪方法 | |
CN105678716B (zh) | 一种地基sar大气干扰相位校正方法及装置 | |
CN112685819A (zh) | 一种用于大坝及滑坡变形gb-sar监测的数据后处理方法及系统 | |
CN107977939B (zh) | 一种基于可靠度的加权最小二乘相位展开计算方法 | |
CN106932773B (zh) | 基于修正嵌入式容积卡尔曼滤波的相位展开算法 | |
CN106707283B (zh) | 基于无味信息滤波的相位展开算法 | |
CN110658521B (zh) | 一种基于缠绕相位的GBInSAR大气校正方法及系统 | |
CN109143234B (zh) | 基于FFT高精度估计干涉相位梯度的InSAR滤波方法及系统 | |
CN111915570A (zh) | 基于反向传播神经网络的大气延迟估计方法 | |
CN112034457B (zh) | 基于干涉条纹方向的多基线高程干涉相位估计方法 | |
CN104808205A (zh) | 基于PhaseLift自聚焦算法的稀疏微波成像方法 |
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 |