CN107544069B - 基于平面近似模型的多基线相位解缠绕方法 - Google Patents

基于平面近似模型的多基线相位解缠绕方法 Download PDF

Info

Publication number
CN107544069B
CN107544069B CN201710748937.0A CN201710748937A CN107544069B CN 107544069 B CN107544069 B CN 107544069B CN 201710748937 A CN201710748937 A CN 201710748937A CN 107544069 B CN107544069 B CN 107544069B
Authority
CN
China
Prior art keywords
phase
gradient
interference
pixel position
interference 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.)
Active
Application number
CN201710748937.0A
Other languages
English (en)
Other versions
CN107544069A (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.)
Xidian University
Original Assignee
Xidian 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 Xidian University filed Critical Xidian University
Priority to CN201710748937.0A priority Critical patent/CN107544069B/zh
Publication of CN107544069A publication Critical patent/CN107544069A/zh
Application granted granted Critical
Publication of CN107544069B publication Critical patent/CN107544069B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Analysis (AREA)
  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)

Abstract

本发明属于雷达信号处理领域，公开了一种基于平面近似模型的多基线相位解缠绕方法，包括：获取干涉相位图，计算干涉相位图的缠绕相位梯度，然后进行基于平面近似的模糊数梯度估计，再通过MCF模型求解模糊数，根据模糊数与缠绕相位计算绝对相位，得到相位解缠绕结果；能够解决在噪声干扰下多基线相位解缠绕的噪声鲁棒性差的问题。

Description

基于平面近似模型的多基线相位解缠绕方法
技术领域
本发明属于雷达信号处理技术领域，尤其涉及一种基于平面近似模型的多基线相位解缠绕方法，适用于多基线干涉合成孔径雷达(InSAR)***在配准、滤波与去平地之后的多基线相位解缠绕问题。
背景技术
作为传统SAR的发展，InSAR***能够全天候、全天时工作，现已广泛应用于地形测绘、地表形变监测与自然灾害检测等诸多方面，因此，InSAR技术的发展一直受到了高度重视。
与单基线InSAR相比，多基线InSAR能够利用多幅干涉相位图，突破“相位连续性假设”，实现复杂地形(例如有峡谷、陡峭山崖等)的高精度测绘，这一优势使得多基线InSAR在军事及民用领域都具有很广泛的应用前景。
在多基线InSAR***进行地形测绘中，多基线相位解缠绕是一个关键步骤，在传统的多基线相位解缠绕方法中，因其理论基础——中国余数定理本身对噪声的影响十分敏感，这直接导致了多基线InSAR在实际中的应用受到很大限制，如果能提高多基线相位解缠绕方法的噪声鲁棒性，将极大的推动复杂地形的高精度测绘以及相关领域的发展。
发明内容
针对上述问题，本发明的目的在于提供一种基于平面近似模型的多基线相位解缠绕方法，能够解决在噪声干扰下多基线相位解缠绕的噪声鲁棒性差的问题。
为达到上述目的，本发明采用如下技术方案予以实现。
一种基于平面近似模型的多基线相位解缠绕方法，所述方法包括如下步骤：
步骤1，从多基线InSAR***中获取两个基线分别对应的干涉相位图，记为第一干涉相位图和第二干涉相位图，两幅干涉相位图的大小相同，且为m×n，其中，m、n分别为大于零的正整数；
步骤2，根据第一干涉相位图中每个像素的缠绕相位，确定第一干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第一相位梯度矩阵，其大小为m×(n-1)，并确定第一干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第一相位梯度矩阵，其大小为(m-1)×n；
根据第二干涉相位图中每个像素的缠绕相位，确定第二干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第二相位梯度矩阵，其大小为m×(n-1)，并确定第二干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第二相位梯度矩阵，其大小为(m-1)×n；
步骤3，设定预设窗口大小，根据预设窗口大小选定给定像素位置对应的估计窗口；
步骤4，给定像素位置(i,j)，1≤i≤m,1≤j≤n-1；
根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；以及估计得到第二干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；
步骤5，给定像素位置(i,j)，1≤i≤m-1,1≤j≤n；
根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；以及估计得到第二干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；
步骤6，将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素以及干涉相位图1在垂直方向的模糊数梯度矩阵中的元素带入设定的优化模型，得到第一干涉相位图的模糊数矩阵1；
将第二干涉相位图在水平方向的模糊数梯度矩阵的元素以及第二干涉相位图在垂直方向的模糊数梯度矩阵的元素带入所述设定的优化模型，得到第二干涉相位图的模糊数矩阵2；
步骤7，给定像素位置(i，j)，1≤i≤m，1≤j≤n；
根据第一干涉相位图在像素位置(i，j)处对应的缠绕相位和第一干涉相位图在像素位置(i，j)处的模糊数，计算得到第一干涉相位图在像素位置(i，j)处的绝对相位；
根据第二干涉相位图在像素位置(i，j)处对应的缠绕相位和第二干涉相位图在像素位置(i，j)处的模糊数，计算得到第二干涉相位图在像素位置(i，j)处的绝对相位；
步骤8，令像素位置(i，j)遍历集合{(i，j)|1≤i≤m，1≤j≤n}中的所有元素，重复执行步骤7，从而得到第一干涉相位图的绝对相位矩阵，以及第二干涉相位图的绝对相位矩阵。
本发明的有益效果为：(1)本发明所提出的方法能够有效降低噪声对于模糊数梯度估计的影响，从而提高了多基线相位解缠绕方法的噪声鲁棒性；(2)本发明所提出的方法通过有效估计模糊数梯度，从而能够得到相对稀疏的残点分布，为后续最小费用流(MCF)模型的求解提供了良好的铺垫，也可以移植到其它相位解缠绕方法中。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的一种基于平面近似模型的多基线相位解缠绕方法的流程示意图；
图2为本发明方法与传统TSPA方法的处理结果示意图，所有图的横坐标对应相应矩阵的列，所有图的纵坐标对应相应矩阵的行；
图2中，图2(a)为参考绝对相位图；图2(b)为仿真干涉相位图一；图2(c)为仿真干涉相位图二；图2(d)为由TSPA方法求解得到的相位解缠绕结果图；图2(e)为由本发明方法求解得到的相位解缠绕结果图；图2(f)为图2(d)与2(a)的差值图；图2(g)为图2(e)与2(a)的差值图；
图3为本发明方法与传统TSPA方法的处理结果示意图，所有图的横坐标对应相应矩阵的列，所有图的纵坐标对应相应矩阵的行；
图3中，图3(a)为Google Earth图像图3(b)为短基线干涉相位图；图3(c)为长基线干涉相位图；图3(d)为TSPA方法的相位解缠绕结果图；图3(e)为本发明方法的相位解缠绕结果图；图3(f)为TSPA方法所得到的残点分布局部图；图3(g)为本发明方法得到的残点分布局部图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
从多基线InSAR***可以获得多幅干涉相位图，以下以基线数为2的情况为例说明本发明，对于基线数大于2的情况，可以同理推得。
本发明实施例提供一种基于平面近似模型的多基线相位解缠绕方法，如图1所示，所述方法包括如下步骤：
步骤1，从多基线InSAR***中获取两个基线分别对应的干涉相位图，记为第一干涉相位图和第二干涉相位图，两幅干涉相位图的大小相同，且为m×n，其中，m、n分别为大于零的正整数。
设多基线InSAR***中，有两个基线，即为基线1和基线2，基线1的基线长度为B1，基线2的基线长度为B2，具有两个基线的多基线InSAR***可以获得两幅干涉相位图，基线1对应第一干涉相位图，基线2对应第二干涉相位图，两幅干涉相位图的尺寸均为m×n，干涉相位图中的元素称为缠绕相位，第一干涉相位图中位置(i,j)的缠绕相位记为
Figure GDA0002396247630000051
第二干涉相位图中位置(i,j)的缠绕相位记为
Figure GDA0002396247630000052
步骤2，根据第一干涉相位图中每个像素的缠绕相位，确定第一干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第一相位梯度矩阵，其大小为m×(n-1)，并确定第一干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第一相位梯度矩阵，其大小为(m-1)×n；
根据第二干涉相位图中每个像素的缠绕相位，确定第二干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第二相位梯度矩阵，其大小为m×(n-1)，并确定第二干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第二相位梯度矩阵，其大小为(m-1)×n。
步骤2具体包括如下子步骤：
(2a)获取第一干涉相位图中像素位置(i,j)处的缠绕相位
Figure GDA0002396247630000053
1≤i≤m,1≤j<n；获取第二干涉相位图中像素位置(i,j)处的缠绕相位
Figure GDA0002396247630000061
1≤i≤m,1≤j<n；像素位置(i,j)处的像素值为像素位置(i,j)处的缠绕相位；
(2b)记第一干涉相位图的水平缠绕第一相位梯度矩阵为
Figure GDA0002396247630000062
中第i行第j列元素
Figure GDA0002396247630000063
1≤i≤m，1≤j≤n-1；
记第二干涉相位图的水平缠绕第二相位梯度矩阵为
Figure GDA0002396247630000064
中第i行第j列元素
Figure GDA0002396247630000065
1≤i≤m，1≤j≤n-1；
(2c)记第一干涉相位图的垂直缠绕第一相位梯度矩阵为
Figure GDA0002396247630000066
中第i行第j列元素
Figure GDA0002396247630000067
1≤i≤m-1，1≤j≤n；
记第二干涉相位图的垂直缠绕第二相位梯度矩阵为
Figure GDA0002396247630000068
中第i行第j列元素
Figure GDA0002396247630000069
1≤i≤m-1，1≤j≤n。
步骤3，设定预设窗口大小，根据预设窗口大小选定给定像素位置对应的估计窗口。
步骤3中，设定预设窗口大小为(2p+1)×(2q+1)；
根据地形的先验信息确定正整数p和q作为预设窗口尺寸参数。预设窗口大小的选择(即参数p和q的大小)，通常需要根据对地形的先验信息进行选择，预设窗口选取的原则是：保证预设窗口内像素所对应的地形变化相对平缓，可以近似的认为是在同一个平面内。
(3a)对于水平缠绕第一相位梯度矩阵或者水平缠绕第二相位梯度矩阵中给定的像素位置(i,j)，1≤i≤m，1≤j≤n-1，其对应的估计窗口内所包含的水平缠绕相位梯度记为
Figure GDA00023962476300000610
或
Figure GDA00023962476300000611
其中，a0≤i'≤a1，b0≤j'≤b1，a0，a1，b0，b1的值由下式确定：
Figure GDA00023962476300000612
Figure GDA0002396247630000071
Figure GDA0002396247630000072
Figure GDA0002396247630000073
(3b)对于垂直缠绕第一相位梯度矩阵或者垂直缠绕第二相位梯度矩阵中给定的像素位置(i,j)1≤i≤m-1，1≤j≤n，其对应的估计窗口内所包含的垂直缠绕相位梯度记为
Figure GDA0002396247630000074
或
Figure GDA0002396247630000075
其中，c0≤i″≤c1，d0≤j″≤d1，c0，c1，d0，d1的值由下式确定：
Figure GDA0002396247630000076
Figure GDA0002396247630000077
Figure GDA0002396247630000078
Figure GDA0002396247630000079
步骤4，本步骤的目的是，利用窗内所有水平缠绕相位梯度，对给定位置的水平模糊数梯度进行估计。
给定像素位置(i,j)，1≤i≤m,1≤j≤n-1；
根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；以及估计得到第二干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)。
步骤4具体包括：
(4a)对于给定像素位置(i,j)，第一干涉相位图的水平缠绕相位梯度
Figure GDA00023962476300000710
和第二干涉相位图的水平缠绕相位梯度
Figure GDA00023962476300000711
以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗口内所有水平缠绕相位梯度位置(i',j')的集合
Figure GDA0002396247630000081
定义Rpq误差
Figure GDA0002396247630000082
Figure GDA0002396247630000083
其中，
Figure GDA0002396247630000084
Figure GDA0002396247630000085
Figure GDA0002396247630000086
Figure GDA0002396247630000087
Figure GDA0002396247630000088
其中，符号[·]表示取四舍五入，
Figure GDA0002396247630000089
为第一干涉相位图在像素位置(i,j)的水平方向模糊数梯度，
Figure GDA00023962476300000810
为第二干涉相位图在像素位置(i,j)的水平方向模糊数梯度，B1表示基线1的长度，B2表示基线2的长度，
Figure GDA00023962476300000813
和
Figure GDA00023962476300000811
为水平梯度修正因子，用于修正窗内像素位置(i,j)处的水平绝对相位梯度与窗内其它位置(i',j')处的水平绝对相位梯度因为相位缠绕所带来的梯度不一致；
(4b)求解如下优化模型：
Figure GDA00023962476300000812
其中，Z表示整数集合，第一干涉相位图在像素位置(i,j)的水平方向模糊数梯度估计值
Figure GDA0002396247630000091
和第二干涉相位图在像素位置(i,j)的水平方向模糊数梯度估计值
Figure GDA0002396247630000092
是该优化模型要求解的变量。
在实际应用中，可以根据关于地形的先验知识，确定
Figure GDA0002396247630000093
和
Figure GDA0002396247630000094
的取值范围，通过在可能的取值范围内遍历所有
Figure GDA0002396247630000095
和
Figure GDA0002396247630000096
的取值组合，计算相应的Rpq误差，经过比较即可获得满足优化模型的最优解。
步骤5，本步骤的目的是，利用窗内所有垂直缠绕相位梯度，对给定位置的垂直模糊数梯度进行估计。
给定像素位置(i,j)，1≤i≤m-1,1≤j≤n；
根据给定像素位置(i,j)，以及给定像素位置(i,j)对应的估计窗口，估计得到第一干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；以及估计得到第二干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n。
步骤5具体包括：
(5a)对于给定像素位置(i,j)，第一干涉相位图的垂直缠绕相位梯度
Figure GDA0002396247630000097
和第二干涉相位图的垂直缠绕相位梯度
Figure GDA0002396247630000098
以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗内所有垂直相位梯度位置(i″,j″)的集合
Figure GDA0002396247630000099
定义Rpq误差
Figure GDA00023962476300000910
Figure GDA00023962476300000911
其中，
Figure GDA00023962476300000912
Figure GDA0002396247630000101
Figure GDA0002396247630000102
Figure GDA0002396247630000103
Figure GDA0002396247630000104
其中，符号[·]表示取四舍五入，
Figure GDA0002396247630000105
为第一干涉相位图在像素位置(i,j)的垂直方向模糊数梯度，
Figure GDA0002396247630000106
为第二干涉相位图在像素位置(i,j)的垂直方向模糊数梯度，B1表示基线1的长度，B2表示基线2的长度，
Figure GDA0002396247630000107
和
Figure GDA0002396247630000108
为垂直梯度修正因子，用于修正窗内像素位置(i,j)处的垂直绝对相位梯度与窗内其它位置(i”,j”)处的垂直绝对相位梯度因为相位缠绕所带来的梯度不一致；
(5b)求解如下优化模型：
Figure GDA0002396247630000109
其中，Z表示整数集合，第一干涉相位图在像素位置(i,j)的垂直方向模糊数梯度估计值
Figure GDA00023962476300001010
和第二干涉相位图在像素位置(i,j)的垂直方向模糊数梯度估计值
Figure GDA00023962476300001011
是该优化模型要求解的变量。
在实际应用中，可以根据关于地形的先验知识，确定
Figure GDA00023962476300001012
和
Figure GDA00023962476300001013
的取值范围，通过在可能的取值范围内遍历所有
Figure GDA00023962476300001014
和
Figure GDA00023962476300001015
的取值组合，计算相应的Rpq误差，经过比较即可获得满足优化模型的最优解。
步骤6，将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素以及干涉相位图1在垂直方向的模糊数梯度矩阵中的元素带入设定的优化模型，得到第一干涉相位图的模糊数矩阵1；
将第二干涉相位图在水平方向的模糊数梯度矩阵的元素以及第二干涉相位图在垂直方向的模糊数梯度矩阵的元素带入所述设定的优化模型，得到第二干涉相位图的模糊数矩阵2。
步骤6具体包括：
(6a)将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素
Figure GDA0002396247630000111
以及第一干涉相位图在垂直方向的模糊数梯度矩阵中的元素
Figure GDA0002396247630000112
带入如下设定的优化模型：
Figure GDA0002396247630000113
其中，Z为整数集合，
Figure GDA0002396247630000114
和
Figure GDA0002396247630000115
表示设定的权值(它们可以利用InSAR图像的质量图生成，当不需要加权的时候，取
Figure GDA0002396247630000116
即可)，
Figure GDA0002396247630000117
和
Figure GDA0002396247630000118
为自由变量，k1(i,j)、k1(i+1,j)、k1(i,j+1)分别为第一干涉相位图在像素位置(i,j)、(i+1,j)、(i,j+1)处的模糊数，它们是优化模型的输出。
该优化模型可以利用线性规划优化工具进行求解，求解结果记为模糊数矩阵1，对应于第一干涉相位图，该矩阵的尺寸为m×n。
(6b)将第二干涉相位图在水平方向的模糊数梯度矩阵中的元素
Figure GDA0002396247630000119
以及第二干涉相位图在垂直方向的模糊数梯度矩阵中的元素
Figure GDA00023962476300001110
带入如下设定的优化模型：
Figure GDA0002396247630000121
其中，Z为整数集合，
Figure GDA0002396247630000122
和
Figure GDA0002396247630000123
表示设定的权值(它们可以利用InSAR图像的质量图生成，当不需要加权的时候，取
Figure GDA0002396247630000124
即可)，
Figure GDA0002396247630000125
和
Figure GDA0002396247630000126
为自由变量，k2(i,j)、k2(i+1,j)、k2(i,j+1)分别为第二干涉相位图在像素位置(i,j)、(i+1,j)、(i,j+1)处的模糊数它们是优化模型的输出。
该优化模型可以利用线性规划优化工具进行求解，求解结果为模糊数矩阵2，对应于第二干涉相位图，该矩阵的尺寸为m×n。
步骤7，给定像素位置(i,j)，1≤i≤m,1≤j≤n；
根据第一干涉相位图在像素位置(i,j)处对应的缠绕相位和第一干涉相位图在像素位置(i,j)处的模糊数，计算得到第一干涉相位图在像素位置(i,j)处的绝对相位；
根据第二干涉相位图在像素位置(i,j)处对应的缠绕相位和第二干涉相位图在像素位置(i,j)处的模糊数，计算得到第二干涉相位图在像素位置(i,j)处的绝对相位。
步骤7具体包括：
根据第一干涉相位图在像素位置(i,j)处对应的缠绕相位
Figure GDA0002396247630000127
和第一干涉相位图在像素位置(i,j)处的模糊数k1(i,j)，计算得到第一干涉相位图在像素位置(i,j)处的绝对相位
Figure GDA0002396247630000128
根据第二干涉相位图在像素位置(i,j)处对应的缠绕相位
Figure GDA0002396247630000129
和第二干涉相位图在像素位置(i,j)处的模糊数k2(i,j)，计算得到第二干涉相位图在像素位置(i,j)处的绝对相位
Figure GDA0002396247630000131
步骤8，令像素位置(i,j)遍历集合{(i,j)|1≤i≤m,1≤j≤n}中的所有元素，重复执行步骤7，从而得到第一干涉相位图的绝对相位矩阵(相位解缠绕结果)，以及第二干涉相位图的绝对相位矩阵(相位解缠绕结果)。
本发明的有效性可通过以下实验作进一步说明。
(一)InSAR仿真数据的实验
1、数据说明
仿真数据采用Colorado的Isolation Peak地区的数字高程模型DEM，分别采用基线长度112.1米和基线长度389.2米得到短基线和长基线的干涉相位图，并在干涉相位图中分别加入了相干系数为0.7和0.65的随机噪声。
2、仿真内容和结果分析
为了展示本发明方法的有效性，将本发明方法与传统TSPA方法进行了对比实验。图2(a)所示为由DEM得到的参考绝对相位。为了展示本发明方法克服噪声的能力，在进行干涉相位图仿真时加入了随机噪声，图2(b)为由图2(a)得到的干涉相位图，垂直基线(normalbaseline)长度为112.1米，加入了相干系数为0.7的随机噪声，图2(c)为由图2(a)得到的干涉相位图，垂直基线(normal baseline)长度为389.2米，加入了相干系数为0.65的随机噪声。图2(d)是采用TSPA方法求解的解缠绕结果，图2(e)是采用本发明方法求解的解缠绕结果，从这两个解缠绕相位图可以看出，TSPA方法的解缠绕结果中存在一些较为明显的色斑，而本发明方法几乎不存在类似的色斑，解缠绕结果看起来比TSPA方法相对“光滑”。图2(f)是TSPA方法的解缠绕结果与参考绝对相位的差值，差值范围为-9.3558～11.9075，图2(g)是本发明方法解缠绕结果与参考绝对相位的差值，差值范围为-6.8004～8.8865。图2(f)和图2(g)反映了求解结果与理想的参考绝对相位之间的误差大小，根据误差范围可以看出，本发明方法比传统TSPA方法的误差小。
为了能够定量的反映两种方法在精度上提高，本发明采用下式衡量解缠绕结果的精度ξ：
Figure GDA0002396247630000141
其中，h是由参考绝对相位得到的向量，
Figure GDA0002396247630000142
由解缠绕结果得到的绝对相位向量。符号||·||表示2范数。经过计算，图2(f)的ξ为0.1003，图2(g)的ξ值为0.0769，这说明了本发明方法比TSPA方法具有更高的相位解缠绕精度。
(二)实测数据实验
1、数据说明
实验数据采用我国甘肃省兰州地区的一组实测数据，基线长度分别为6.54米和40.29米。获取该干涉数据的SAR图像由TerraSAR-X***在2010年于2月25日(轨道14975)、3月8日(轨道15142)与3月19(轨道15309)分别获得，***轨道高度为512千米，入射角为28.7°，波长为0.03米。
2、仿真内容和结果分析
图3(a)为该地区的Google Eearth图像(该地区的纬度约为36.0405°，经度约为103.7802°)。该地区有许多陡峭的山崖与山谷，对于相位解缠绕而言，这样的复杂地形是有难度的。图3(b)是短基线干涉相位图，基线长度是6.51米，图3(c)为长基线干涉相位图，基线长度为40.29米。图3(d)为TSPA方法的解缠绕结果，图3(e)为本发明方法的解缠绕结果，为了方便对比，图3(d)与图3(e)采用了相同的颜色条(Colorbar)。通过对比可以看出，在TSPA方法的相位解缠绕图像中存在一些色斑(在图中用白色的矩形方框标识)，这说明在这些局部，解缠绕结果精度较低，而本发明方法在这些区域的图像质量却相对良好，并且通过将两种方法的解缠绕结果与图3(a)所示的谷歌地球(Google Earth)图像进行对比，本发明方法的求解结果更加接近于真实地形。为了进一步说明TSPA方法与本发明方法求解精度不同的原因，本发明根据两种方法得到的梯度估计值分别计算了残点，并在图3(f)和图3(g)中分别给出了白色方框内区域的残点分布图。对比图3(f)和图3(g)，TSPA方法估计的模糊数梯度产生了90078个残点，整个区域几乎布满了残点，本发明方法得到的残点数量为仅为18050，相对集中在干涉相位图中条纹分布变化不规律的位置。在TSPA方法和本发明方法中，在完成了模糊数梯度的估计之后，都采用了MCF模型对模糊数进行求解，然而，该优化模型在残点多为偶极子且分布较为稀疏的情况下求解的精度较高，对于残点密集的区域求解精度将会大大下降，因此，本发明方法通过众值滤波得到了更为合理的模糊数梯度估计值和相对稀疏的残点分布，因而有利于MCF模型提高求解精度。总之，通过本实验，可以说明对于实测的复杂地形，本发明方法相比于传统TSPA方法，相位解缠绕的精度得到了有效的提升。
综上，我们分别采用了仿真数据与实测数据验证了本发明的有效性。
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。

Claims (7)

1.一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，所述方法包括如下步骤：
步骤1，从多基线InSAR***中获取两个基线分别对应的干涉相位图，记为第一干涉相位图和第二干涉相位图，两幅干涉相位图的大小相同，且为m×n，其中，m、n分别为大于零的正整数；
步骤2，根据第一干涉相位图中每个像素的缠绕相位，确定第一干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第一相位梯度矩阵，其大小为m×(n-1)，并确定第一干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第一相位梯度矩阵，其大小为(m-1)×n；
根据第二干涉相位图中每个像素的缠绕相位，确定第二干涉相位图中水平方向相邻像素间的水平缠绕相位梯度，记为水平缠绕第二相位梯度矩阵，其大小为m×(n-1)，并确定第二干涉相位图中垂直方向相邻像素间的垂直缠绕相位梯度，记为垂直缠绕第二相位梯度矩阵，其大小为(m-1)×n；
步骤3，设定预设窗口大小，根据预设窗口大小选定给定像素位置对应的估计窗口；
步骤4，给定像素位置(i，j)，1≤i≤m，1≤j≤n-1；
根据给定像素位置(i，j)，以及给定像素位置(i，j)对应的估计窗口，估计得到第一干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；以及估计得到第二干涉相位图的水平方向模糊数梯度矩阵，其大小为m×(n-1)；
步骤5，给定像素位置(i，j)，1≤i≤m-1，1≤j≤n；
根据给定像素位置(i，j)，以及给定像素位置(i，j)对应的估计窗口，估计得到第一干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；以及估计得到第二干涉相位图的垂直方向模糊数梯度矩阵，其大小为(m-1)×n；
步骤6，将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素以及干涉相位图1在垂直方向的模糊数梯度矩阵中的元素带入设定的优化模型，得到第一干涉相位图的模糊数矩阵1；
将第二干涉相位图在水平方向的模糊数梯度矩阵的元素以及第二干涉相位图在垂直方向的模糊数梯度矩阵的元素带入所述设定的优化模型，得到第二干涉相位图的模糊数矩阵2；
步骤7，给定像素位置(i，j)，1≤i≤m，1≤j≤n；
根据第一干涉相位图在像素位置(i，j)处对应的缠绕相位和第一干涉相位图在像素位置(i，j)处的模糊数，计算得到第一干涉相位图在像素位置(i，j)处的绝对相位；
根据第二干涉相位图在像素位置(i，j)处对应的缠绕相位和第二干涉相位图在像素位置(i，j)处的模糊数，计算得到第二干涉相位图在像素位置(i，j)处的绝对相位；
步骤8，令像素位置(i，j)遍历集合{(i，j)|1≤i≤m，1≤j≤n}中的所有元素，重复执行步骤7，从而得到第一干涉相位图的绝对相位矩阵，以及第二干涉相位图的绝对相位矩阵。
2.根据权利要求1所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤2具体包括如下子步骤：
(2a)获取第一干涉相位图中像素位置(i，j)处的缠绕相位
Figure FDA0002396247620000021
Figure FDA0002396247620000022
获取第二干涉相位图中像素位置(i，j)处的缠绕相位
Figure FDA0002396247620000023
像素位置(i，j)处的像素值为像素位置(i，j)处的缠绕相位；
(2b)记第一干涉相位图的水平缠绕第一相位梯度矩阵为
Figure FDA0002396247620000031
中第i行第j列元素
Figure FDA0002396247620000032
记第二干涉相位图的水平缠绕第二相位梯度矩阵为
Figure FDA0002396247620000033
中第i行第j列元素
Figure FDA0002396247620000034
(2c)记第一干涉相位图的垂直缠绕第一相位梯度矩阵为
Figure FDA0002396247620000035
中第i行第j列元素
Figure FDA0002396247620000036
记第二干涉相位图的垂直缠绕第二相位梯度矩阵为
Figure FDA0002396247620000037
中第i行第j列元素
Figure FDA0002396247620000038
3.根据权利要求2所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤3中，设定预设窗口大小为(2p+1)×(2q+1)；
(3a)对于水平缠绕第一相位梯度矩阵或者水平缠绕第二相位梯度矩阵中给定的像素位置(i，j)，1≤i≤m，1≤j≤n-1，其对应的估计窗口内所包含的水平缠绕相位梯度记为
Figure FDA0002396247620000039
或者
Figure FDA00023962476200000310
其中，a0≤i＇≤a1，b0≤j＇≤b1，a0，a1，b0，b1的值由下式确定：
Figure FDA00023962476200000311
Figure FDA00023962476200000312
Figure FDA00023962476200000313
Figure FDA00023962476200000314
(3b)对于垂直缠绕第一相位梯度矩阵或者垂直缠绕第二相位梯度矩阵中给定的像素位置(i，j)，1≤i≤m-1，1≤j≤n，其对应的估计窗口内所包含的垂直缠绕相位梯度记为
Figure FDA00023962476200000315
或者
Figure FDA00023962476200000316
其中，c0≤i″≤c1，d0≤j″≤d1，c0，c1，d0，d1的值由下式确定：
Figure FDA0002396247620000041
Figure FDA0002396247620000042
Figure FDA0002396247620000043
Figure FDA0002396247620000044
4.根据权利要求3所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤4具体包括：
(4a)对于给定像素位置(i，j)，第一干涉相位图的水平缠绕相位梯度
Figure FDA0002396247620000045
和第二干涉相位图的水平缠绕相位梯度
Figure FDA0002396247620000046
以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗口内所有水平缠绕相位梯度位置(i＇，j＇)的集合
Figure FDA0002396247620000047
定义Rpq误差
Figure FDA0002396247620000048
Figure FDA0002396247620000049
其中，
Figure FDA00023962476200000410
Figure FDA0002396247620000051
Figure FDA0002396247620000052
Figure FDA0002396247620000053
Figure FDA0002396247620000054
其中，符号[·]表示取四舍五入，
Figure FDA0002396247620000055
为第一干涉相位图在像素位置(i，j)的水平方向模糊数梯度，
Figure FDA0002396247620000056
为第二干涉相位图在像素位置(i，j)的水平方向模糊数梯度，B1表示基线1的长度，B2表示基线2的长度，
Figure FDA0002396247620000057
和
Figure FDA0002396247620000058
为水平梯度修正因子，用于修正窗内像素位置(i，j)处的水平绝对相位梯度与窗内其它位置(i＇，j＇)处的水平绝对相位梯度因为相位缠绕所带来的梯度不一致；
(4b)求解如下优化模型：
Figure FDA0002396247620000059
其中，Z表示整数集合，第一干涉相位图在像素位置(i，j)的水平方向模糊数梯度估计值
Figure FDA00023962476200000510
和第二干涉相位图在像素位置(i，j)的水平方向模糊数梯度估计值
Figure FDA00023962476200000511
是该优化模型要求解的变量。
5.根据权利要求3所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤5具体包括：
(5a)对于给定像素位置(i，j)，第一干涉相位图的垂直缠绕相位梯度
Figure FDA00023962476200000512
和第二干涉相位图的垂直缠绕相位梯度
Figure FDA00023962476200000513
以及尺寸为(2p+1)×(2q+1)的估计窗口，记估计窗内所有垂直相位梯度位置(i″，j″)的集合
Figure FDA00023962476200000514
定义Rpq误差
Figure FDA0002396247620000061
Figure FDA0002396247620000062
其中，
Figure FDA0002396247620000063
Figure FDA0002396247620000064
Figure FDA0002396247620000065
Figure FDA0002396247620000066
Figure FDA0002396247620000067
其中，符号[·]表示取四舍五入，
Figure FDA0002396247620000068
为第一干涉相位图在像素位置(i，j)的垂直方向模糊数梯度，
Figure FDA0002396247620000069
为第二干涉相位图在像素位置(i，j)的垂直方向模糊数梯度，B1表示基线1的长度，B2表示基线2的长度，
Figure FDA00023962476200000610
和
Figure FDA00023962476200000611
为垂直梯度修正因子，用于修正窗内像素位置(i，j)处的垂直绝对相位梯度与窗内其它位置(i″，j″)处的垂直绝对相位梯度因为相位缠绕所带来的梯度不一致；
(5b)求解如下优化模型：
Figure FDA00023962476200000612
其中，Z表示整数集合，第一干涉相位图在像素位置(i，j)的垂直方向模糊数梯度估计值
Figure FDA00023962476200000613
和第二干涉相位图在像素位置(i，j)的垂直方向模糊数梯度估计值
Figure FDA00023962476200000614
是该优化模型要求解的变量。
6.根据权利要求4或5所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤6具体包括：
(6a)将第一干涉相位图在水平方向的模糊数梯度矩阵中的元素
Figure FDA0002396247620000071
以及第一干涉相位图在垂直方向的模糊数梯度矩阵中的元素
Figure FDA0002396247620000072
带入如下设定的优化模型：
Figure FDA0002396247620000073
其中，Z为整数集合，
Figure FDA0002396247620000074
和
Figure FDA0002396247620000075
表示设定的权值，
Figure FDA0002396247620000076
和
Figure FDA0002396247620000077
为自由变量，k1(i，j)、k1(i+1，j)、k1(i，j+1)分别为第一干涉相位图在像素位置(i，j)、(i+1，j)、(i，j+1)处的模糊数；
(6b)将第二干涉相位图在水平方向的模糊数梯度矩阵中的元素
Figure FDA0002396247620000078
以及第二干涉相位图在垂直方向的模糊数梯度矩阵中的元素
Figure FDA0002396247620000079
带入如下设定的优化模型：
Figure FDA00023962476200000710
其中，Z为整数集合，
Figure FDA00023962476200000711
和
Figure FDA00023962476200000712
表示设定的权值，
Figure FDA00023962476200000713
和
Figure FDA00023962476200000714
为自由变量，k2(i，j)、k2(i+1，j)、k2(i，j+1)分别为第二干涉相位图在像素位置(i，j)、(i+1，j)、(i，j+1)处的模糊数。
7.根据权利要求6所述的一种基于平面近似模型的多基线相位解缠绕方法，其特征在于，步骤7具体包括：
根据第一干涉相位图在像素位置(i，j)处对应的缠绕相位
Figure FDA0002396247620000081
和第一干涉相位图在像素位置(i，j)处的模糊数k1(i，j)，计算得到第一干涉相位图在像素位置(i，j)处的绝对相位
Figure FDA0002396247620000084
根据第二干涉相位图在像素位置(i，j)处对应的缠绕相位
Figure FDA0002396247620000082
和第二干涉相位图在像素位置(i，j)处的模糊数k2(i，j)，计算得到第二干涉相位图在像素位置(i，j)处的绝对相位
Figure FDA0002396247620000083
CN201710748937.0A 2017-08-28 2017-08-28 基于平面近似模型的多基线相位解缠绕方法 Active CN107544069B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710748937.0A CN107544069B (zh) 2017-08-28 2017-08-28 基于平面近似模型的多基线相位解缠绕方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710748937.0A CN107544069B (zh) 2017-08-28 2017-08-28 基于平面近似模型的多基线相位解缠绕方法

Publications (2)

Publication Number Publication Date
CN107544069A CN107544069A (zh) 2018-01-05
CN107544069B true CN107544069B (zh) 2020-07-03

Family

ID=60957503

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710748937.0A Active CN107544069B (zh) 2017-08-28 2017-08-28 基于平面近似模型的多基线相位解缠绕方法

Country Status (1)

Country Link
CN (1) CN107544069B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108663678B (zh) * 2018-01-29 2022-01-18 西北农林科技大学 基于混合整数优化模型的多基线InSAR相位解缠算法
CN113791413B (zh) * 2021-09-14 2023-08-01 华北水利水电大学 多基线InSAR分枝定界纯整数规划相位解缠算法
CN117269960B (zh) * 2023-09-12 2024-04-26 中国矿业大学 一种基于梯度优化的快速范数相位解缠方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102621549B (zh) * 2011-10-14 2014-04-16 中国人民解放军国防科学技术大学 多基线/多频段干涉相位解缠频域快速算法
CN103605107B (zh) * 2013-12-03 2016-05-25 西安电子科技大学 基于多基线分布式阵列的波达方向估计方法
CN103885059B (zh) * 2014-01-26 2017-04-05 中国测绘科学研究院 一种多基线干涉合成孔径雷达三维重建方法

Also Published As

Publication number Publication date
CN107544069A (zh) 2018-01-05

Similar Documents

Publication Publication Date Title
CN106772342B (zh) 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法
CN106093939B (zh) 一种基于相位差统计模型的InSAR图像相位解缠方法
CN109633648B (zh) 一种基于似然估计的多基线相位估计装置及方法
CN102621549B (zh) 多基线/多频段干涉相位解缠频域快速算法
US8879871B2 (en) Method for processing images using automatic georeferencing of images derived from a pair of images captured in the same focal plane
CN106249236B (zh) 一种星载InSAR长短基线图像联合配准方法
CN107544069B (zh) 基于平面近似模型的多基线相位解缠绕方法
CN104730519B (zh) 一种采用误差迭代补偿的高精度相位解缠方法
CN109100720B (zh) 一种InSAR地表形变监测方法
CN112051572A (zh) 一种融合多源sar数据三维地表形变监测方法
CN107014313B (zh) 基于s变换脊值的加权最小二乘相位展开的方法及***
CN113589286B (zh) 基于D-LinkNet的无迹卡尔曼滤波相位解缠方法
CN108562900B (zh) 一种基于高程校正的sar图像几何配准方法
CN116338607B (zh) 时间域和空间域两步式InSAR对流层延迟矫正方法
CN114814839B (zh) 基于InSAR相位梯度堆叠的广域地表形变探测方法及***
CN106157258B (zh) 一种星载sar图像几何校正方法
CN113282695B (zh) 一种基于遥感影像的矢量地理信息采集方法和装置
CN110988876B (zh) 闭合鲁棒双基线InSAR相位解缠方法、***及可读存储介质
CN112665529A (zh) 基于条纹密度区域分割和校正的物体三维形貌测量方法
CN113866765A (zh) 基于多成分时间相干模型的PS-InSAR测量方法
CN112946647A (zh) 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置
CN116609779A (zh) 两阶段InSAR多观测处理方法，***及相关设备
CN115546264A (zh) 一种星载InSAR图像精细配准与立体测量方法
CN113124834B (zh) 一种结合多源数据的区域网平差方法、***及存储介质
CN114926524A (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