CN104933673B - 基于解析搜索亚像素偏移量的干涉sar图像精配准方法 - Google Patents

基于解析搜索亚像素偏移量的干涉sar图像精配准方法 Download PDF

Info

Publication number
CN104933673B
CN104933673B CN201510364208.6A CN201510364208A CN104933673B CN 104933673 B CN104933673 B CN 104933673B CN 201510364208 A CN201510364208 A CN 201510364208A CN 104933673 B CN104933673 B CN 104933673B
Authority
CN
China
Prior art keywords
mrow
msub
sub
msup
picture
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
CN201510364208.6A
Other languages
English (en)
Other versions
CN104933673A (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 CN201510364208.6A priority Critical patent/CN104933673B/zh
Publication of CN104933673A publication Critical patent/CN104933673A/zh
Application granted granted Critical
Publication of CN104933673B publication Critical patent/CN104933673B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/14—Transformations for image registration, e.g. adjusting or mapping for alignment of images
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T3/00—Geometric image transformations in the plane of the image
    • G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
    • G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
    • G06T3/4069—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution by subpixel displacements
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/10—Image acquisition modality
    • G06T2207/10032—Satellite or aerial image; Remote sensing
    • G06T2207/10044—Radar image

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Theoretical Computer Science (AREA)
  • Medicines Containing Antibodies Or Antigens For Use As Internal Diagnostic Agents (AREA)
  • Image Analysis (AREA)

Abstract

本发明公开了一种基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，包括以下步骤：(1)采用干涉SAR获取两幅SAR复图像，分别是主图像I1和副图像I2；(2)从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口定义副图像窗口相对于主图像窗口的亚像素偏移量为(dx,dy)，通过对副图像窗口上的像素点(x′+dx,y′+dy)周围的像素点(m,n)进行插值，计算副图像窗口在像素点(x′+dx,y′+dy)的像素值并计算主图像窗口和副图像窗口的二维互相关函数构造连续代价函数为(3)优化连续代价函数，并求解连续代价函数的最大值，该最大值对应的亚像素偏移量为副图像窗口相对于主图像窗口的最佳的亚像素偏移量。本发明不仅提高了干涉SAR图像配准的运算速度，而且获得了更高的配准精度。

Description

基于解析搜索亚像素偏移量的干涉SAR图像精配准方法
技术领域
本发明属于干涉合成孔径雷达图像配准领域，进一步涉及一种基于解析搜索亚像素偏移量的干涉合成孔径雷达(SAR)图像精配准方法。
背景技术
干涉合成孔径雷达(InSAR)图像配准是将不同时间下获得的两幅或多幅含有相同场景或目标的SAR复图像进行匹配、叠加的过程，是InSAR数据处理过程中的一个关键步骤。InSAR利用SAR在不同轨道或不同时间对同一地区进行两次或多次观测，获得地面同一场景的两幅或多幅SAR复图像；InSAR对两幅或多幅SAR复图像完成图像间的配准后，再进行差相位处理产生干涉相位，然后从干涉条纹中获取地形的高度信息和变化信息。传统粗配准方法在像素级的基础上进行图像配准，不能达到更高的图像配准精度要求。为了保证在测量时干涉相位的精度，InSAR图像配准的精度必须达到亚像素级，也就是说所获得的SAR复图像间的偏移量的估计精度要达到更高，我们一般采用插值技术来实现。目前，基于图像数据的InSAR图像的亚像素级配准主要采用互相关函数法和最大谱估计法，这两种方法均利用SAR复图像的幅值信息进行图像配准，其中，互相关函数法是以两幅SAR复图像之间的离散像素偏移量为单位计算相关函数，然后对相关函数采用亚像素级插值，得到对应的最大偏移量；最大谱估计法基于这样一个前提：即当图像准确配准时，对应于复干涉图的频谱中主要条纹对应的模最大；最大谱估计法基于该前提得到一组频谱矩阵，其中每一个频谱值对应着一组偏移参数，再通过对频谱矩阵进行二维亚像素级插值，得出最大的频谱值，进而得到最大的偏移量。互相关函数法和最大谱估计法的实质都是在离散域搜索最大的偏移量，从而获得亚像素级的配准精度。但是，在离散域搜索时所获得的偏移量是亚最优的，一方面由于偏移量的计算精度取决于插值单元的尺寸，并不能完全消除配准误差，而且理论上偏移量的计算精度只能达到插值单元尺寸的一半；另一方面，在匹配之前进行的插值运算会导致很大的计算量，特别是在配准精度要求较高的时候，计算量会更大。
发明内容
基于上述现有技术的不足，本发明的目的在于提出一种基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，该方法构造以亚像素偏移量为变量的连续代价函数，能够获得更高的图像配准精度；同时该方法采用拟牛顿方法中的BFGS算法来优化连续代价函数，能够提高图像配准的运算速度。
为达到上述技术目的，本发明采用以下技术方案予以实现。
一种基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，其特征在于，包括以下步骤：
步骤1，利用干涉SAR获取两幅SAR复图像I1和I2，其中，I1为主图像，I2为副图像；设定副图像I2相对于主图像I1的位置偏移量为(u,v)，并主图像I1和副图像I2的二维互相关函数Ru,v，其中，u为副图像I2相对于主图像I1在水平方向上的偏移量，v为副图像I2相对于主图像I1在垂直方向上的偏移量；
步骤2，从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口定义副图像窗口相对于主图像窗口的亚像素偏移量为(dx,dy)，通过对副图像窗口上的像素点(x′+dx,y′+dy)周围的像素点(m,n)进行插值，计算得到副图像窗口在像素点(x′+dx,y′+dy)的像素值根据副图像窗口在像素点(x′+dx,y′+dy)的像素值计算得到主图像窗口和副图像窗口的二维互相关函数构造以亚像素偏移量(dx,dy)为变量的连续代价函数为其中，(x′,y′)为主图像窗口上给定像素点的坐标；
步骤3，采用拟牛顿方法中的BFGS算法优化连续代价函数设计迭代算法求解连续代价函数的最大值，该最大值对应的亚像素偏移量为副图像窗口相对于主图像窗口在方向X和方向Y上的最佳的亚像素偏移量。
本发明的特点和进一步改进在于：
(1)所述步骤1的具体子步骤为：
1.1利用干涉SAR获取两幅SAR复图像I1和I2，其中，I1为主图像，I2为副图像；设定像素点(x,y)在主图像I1的像素值为I1(x,y)，像素点(x,y)在副图像I2的像素值为I2(x,y)，I1(x,y)和I2(x,y)的表达式分别为：
其中，|·|表示取模；|I1(x,y)|和|I2(x,y)|分别为像素点(x,y)在主图像I1和副图像I2的幅度；和分别为像素点(x,y)在主图像I1和副图像I2的相位；x表示像素点(x,y)在水平方向上的分量，y表示像素点(x,y)在垂直方向上的分量；
1.2计算主图像I1和副图像I2的二维互相关函数R为：
其中，上标*表示共轭；
1.3设定P1和P2分别为主图像I1和副图像I2的能量函数，P1和P2的表达式分别为：
设定副图像I2相对于主图像I1的位置偏移量为(u,v)，则副图像I2相对于主图像I1的位置偏移量为(u,v)时，主图像I1和副图像I2的二维互相关函数Ru,v为：
其中，P2(u,v)表示副图像I2相对于主图像I1的位置偏移量为(u,v)时，副图像I2的能量函数。
(2)所述步骤2的具体子步骤为：
2.1从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口针对亚像素级的图像配准，将副图像窗口相对于主图像窗口的位置偏移的整数部分移除，则副图像窗口相对于主图像窗口仅存在小于一个像素的偏移量，称之为亚像素偏移量，用dx和dy分别表示方向X和方向Y上的亚像素偏移量，|dx|＜1，|dy|＜1，|·|表示取模；设定主图像窗口上给定像素点的坐标为(x′,y′)，则在副图像窗口上与其精确配准的像素点的坐标为(x′+dx,y′+dy)；通过对副图像窗口上的像素点(x′+dx,y′+dy)周围的像素点(m,n)进行插值，计算副图像窗口在像素点(x′+dx,y′+dy)的像素值为复数，其表达式为：
其中，x′,y′,m,n∈N，N表示正整数，dx,dy∈R，R表示实数；是副图像窗口在像素点(m,n)的像素值，为复数，m表示像素点(m,n)在水平方向上的分量，n表示像素点(x,y)在垂直方向上的分量；φ(x′,y′)是插值函数，l表示插值函数长度的一半；
将更新为副图像I2在像素点(x′,y′)的像素值，即将副图像窗口逐点精确配准到主图像窗口
2.2将副图像I2相对于主图像I1的位置偏移量为(u,v)时，副图像I2的能量函数P2(u,v)设定为常数，则副图像窗口相对于主图像窗口的亚像素偏移量为(dx,dy)时，主图像窗口和副图像窗口的二维互相关函数为：
其中，是主图像窗口在像素点(x′,y′)的像素值，为复数；
2.3将副图像窗口在像素点(x′+dx,y′+dy)的像素值的表达式，代入主图像窗口和副图像窗口的二维互相关函数的表达式为：
令x′-m＝p,y′-n＝q，将上式调整为：
其中，Cp,q为副图像窗口相对于主图像窗口存在整数倍像素偏移时的互相关表示，简称为互相关表示，互相关表示Cp,q为复数，其表达式为：
其中，*表示卷积，上标*表示共轭，和分别表示主图像窗口和副图像窗口的傅里叶变换，μ和κ分别表示p和q对应的傅里叶变换频率，表示逆傅里叶变换；
2.4根据步骤2.3，构造以亚像素偏移量为变量的连续代价函数为并将求解连续代价函数的最大值的问题表示为：
其中，|·|表示取模。
(3)所述步骤3的具体子步骤为：
3.1将互相关表示Cp,q表示为：
其中，和分别为互相关表示Cp,q的实部和虚部；
将步骤2.4的求解连续代价函数的最大值的问题简写为：
根据上式，令令
3.2分别求连续代价函数对方向X上的亚像素偏移量dx和方向Y上的亚像素偏移量dy的一阶偏导，表达式如下：
其中|·|表示取模；
3.3采用拟牛顿方法中的BFGS算法来优化步骤3.1的求解连续代价函数的最大值的问题；定义u为亚像素偏移量向量，u＝(dx,dy)T，将连续代价函数简写为|Ru|2；设计迭代公式为：
其中，k表示迭代次数，Sk表示第k次迭代的近似Hessian逆矩阵；δk＝uk+1-uk，uk为第k次迭代的亚像素偏移量向量；γk＝gk+1-gk，为连续代价函数|Ru|2在第k次迭代的亚像素偏移量向量uk处的梯度；
设计以下迭代算法求解连续代价函数|Ru|2的最大值，迭代的具体步骤为：
3.3.1设置初始点u0＝[0,0]T；设置近似Hessian逆矩阵的初值S0＝I2，I2为2×2维单位矩阵，允许误差ε＞0；
3.3.2设置迭代次数k的初始值为0；计算连续代价函数|Ru|2在第k次迭代的亚像素偏移量向量uk处的梯度gk；
3.3.3令dk＝-Skgk，从第k次迭代的亚像素偏移量向量uk出发，沿方向dk进行搜索，求解第k次迭代的步长因子αk，使第k次迭代的步长因子αk满足如下公式：
其中，α为步长因子；
根据求解得到的步长因子αk，令δk＝αkdk；计算第k+1次迭代的亚像素偏移量向量uk+1，uk+1＝uk+δk；
3.3.4计算连续代价函数|Ru|2在第k+1次迭代的亚像素偏移量向量uk+1处的梯度gk+1，若梯度gk+1的范数||gk+1||≤ε，则停止迭代，得到最佳的亚像素偏移量向量uopt＝uk+1，uopt对应的副图像窗口相对于主图像窗口在方向X和方向Y上的最佳的亚像素偏移量分别为dxopt和dyopt；否则，进行下一步操作；
3.3.5令迭代次数k增加1，计算第k+1次迭代的近似Hessian逆矩阵Sk+1，返回步骤3.3.3。
本发明提出了一种基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，该方法可以快速地找到最优的亚像素偏移量。首先，通过联合插值操作和互相关函数搜索过程构造一个以亚像素偏移量为变量的连续代价函数；然后，利用连续代价函数对于亚像素偏移量的梯度信息，采用拟牛顿方法中的BFGS算法来优化连续代价函数，得到最佳的亚像素偏移量；本发明不仅提高了干涉SAR图像配准的运算速度，而且获得了更高的配准精度。
附图表说明
下面结合附图说明和具体实施方式对本发明作进一步详细说明。
图1是本发明的流程图。
图2是TerraSAR-X卫星实测InSAR图像对，其中：
图2a是实验采用的主图像；
图2b是实验采用的副图像。
图3a是区域A的幅度图；
图3b是区域A配准后的干涉相位图；
图3c是本发明方法、粗配准法、互相关法以及最大谱估计法对区域A的相干直方图，横坐标为相干系数，纵坐标为像素数目；
图4a是区域B的幅度图；
图4b是区域B配准后的干涉相位图；
图4c是本发明方法、粗配准法、互相关法以及最大谱估计法对区域B的相干直方图，横坐标为相干系数，纵坐标为像素数目；
图5a是区域C的幅度图；
图5b是区域C配准后的干涉相位图；
图5c是本发明方法、粗配准法、互相关法以及最大谱估计法对区域C的相干直方图，横坐标为相干系数，纵坐标为像素数目。
具体实施方式
参照图1，本发明的基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，包括以下步骤：
步骤1：利用干涉SAR获取两幅SAR复图像I1和I2，其中，I1为主图像，I2为副图像；设定副图像I2相对于主图像I1的位置偏移量为(u,v)，并计算主图像I1和副图像I2的二维互相关函数Ru,v，其中，u为副图像I2相对于主图像I1在水平方向上的偏移量，v为副图像I2相对于主图像I1在垂直方向上的偏移量；
步骤1的具体子步骤为：
1.1利用干涉SAR获取两幅SAR复图像I1和I2，其中，I1为主图像，I2为副图像；设定像素点(x,y)在主图像I1的像素值为I1(x,y)，像素点(x,y)在副图像I2的像素值为I2(x,y)，I1(x,y)和I2(x,y)的表达式分别为：
其中，|·|表示取模；|I1(x,y)|和|I2(x,y)|分别为像素点(x,y)在主图像I1和副图像I2的幅度；和分别为像素点(x,y)在主图像I1和副图像I2的相位；x表示像素点(x,y)在水平方向上的分量，y表示像素点(x,y)在垂直方向上的分量；
1.2计算主图像I1和副图像I2的二维互相关函数R为：
其中，上标*表示共轭；
1.3设定P1和P2分别为主图像I1和副图像I2的能量函数，P1和P2的表达式分别为：
设定副图像I2相对于主图像I1的位置偏移量为(u,v)，则副图像I2相对于主图像I1的位置偏移为(u,v)时，主图像I1和副图像I2的二维互相关函数Ru,v为：
其中，P2(u,v)表示副图像I2相对于主图像I1的位置偏移量为(u,v)时，副图像I2的能量函数。
根据Cauchy-Schwarz不等式，如果副图像I2相对于主图像I1的位置偏移量为(u0,v0)时，副图像I2准确配准主图像I1，则主图像I1和副图像I2的二维互相关函数取得最大值为计算副图像I2相对于主图像I1的所有可能的位置偏移量时的主图像I1和副图像I2的二维互相关函数，找到其中的最大值，则该最大值所对应的偏移量，就是副图像I2相对于主图像I1的真实位置偏移量。
步骤2：从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口定义副图像窗口相对于主图像窗口的亚像素偏移量为(dx,dy)，通过对副图像窗口上的像素点(x′+dx,y′+dy)周围的像素点(m,n)进行插值，计算得到副图像窗口在像素点(x′+dx,y′+dy)的像素值根据副图像窗口在像素点(x′+dx,y′+dy)的像素值计算得到主图像窗口和副图像窗口的二维互相关函数构造以亚像素偏移量(dx,dy)为变量的连续代价函数为其中，(x′,y′)为主图像窗口上给定像素点的坐标。
步骤2的具体子步骤为：
2.1从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口针对亚像素级的图像配准，将副图像窗口相对于主图像窗口的位置偏移的整数部分移除，则副图像窗口相对于主图像窗口仅存在小于一个像素的偏移量，称之为亚像素偏移量，用dx和dy分别表示副图像窗口相对于主图像窗口在方向X和方向Y上的亚像素偏移量，|dx|＜1，|dy|＜1，|·|表示取模；设定主图像窗口上给定像素点的坐标为(x′,y′)，则在副图像窗口上与其精确配准的像素点的坐标为(x′+dx,y′+dy)；通过对副图像窗口上的像素点(x′+dx,y′+dy)周围的像素点(m,n)进行插值，计算副图像窗口在像素点(x′+dx,y′+dy)的像素值为复数，其表达式为：
其中，x′,y′,m,n∈N，N表示正整数，dx,dy∈R，R表示实数；是副图像窗口在像素点(m,n)的像素值，为复数，m表示像素点(m,n)在水平方向上的分量，n表示像素点(x,y)在垂直方向上的分量；φ(x′,y′)是插值函数，l表示插值函数长度的一半；
将更新为副图像I2在像素点(x′,y′)的像素值，即将副图像窗口逐点精确配准到主图像窗口
2.2将副图像I2相对于主图像I1的位置偏移量为(u,v)时，副图像I2的能量函数P2(u,v)设定为常数，则副图像窗口相对于主图像窗口的亚像素偏移量为(dx,dy)时，主图像窗口和副图像窗口的二维互相关函数为：
其中，是主图像窗口在像素点(x′,y′)的像素值，为复数；
2.3将副图像窗口在像素点(x′+dx,y′+dy)的像素值的表达式，代入主图像窗口和副图像窗口的二维互相关函数的表达式为：
令x′-m＝p,y′-n＝q，将上式调整为：
其中，Cp,q为副图像窗口相对于主图像窗口存在整数倍像素偏移时的互相关表示，简称为互相关表示，互相关表示Cp,q为复数，其表达式为：
其中，算符*表示卷积，上标*表示共轭，和分别表示主图像窗口和副图像窗口的傅里叶变换，μ和κ分别表示p和q对应的傅里叶变换频率，表示逆傅里叶变换；
2.4根据步骤2.3，构造以亚像素偏移量为变量的连续代价函数为并将求解连续代价函数的最大值的问题表示为：
其中，|·|表示取模。
步骤3，采用拟牛顿方法中的BFGS算法优化连续代价函数设计迭代算法求解连续代价函数的最大值，该最大值对应的亚像素偏移量为副图像窗口相对于主图像窗口在方向X和方向Y上的最佳的亚像素偏移量。
步骤3的具体子步骤为：
3.1将互相关表示Cp,q表示复数形式：
其中，和分别为互相关表示Cp,q的实部和虚部；
将步骤2.4的求解连续代价函数的最大值的问题简写为：
根据上式，令令
3.2分别求连续代价函数对方向X上的亚像素偏移量dx和方向Y上的亚像素偏移量dy的一阶偏导，表达式如下：
其中|·|表示取模；
3.3采用拟牛顿方法中的BFGS算法来优化步骤3.1的求解连续代价函数的最大值的问题；定义u为亚像素偏移量向量，u＝(dx,dy)T，将连续代价函数简写为|Ru|2；设计迭代公式为：
其中，k表示迭代次数，Sk表示第k次迭代的近似Hessian逆矩阵；δk＝uk+1-uk，uk为第k次迭代的亚像素偏移量向量；γk＝gk+1-gk，为连续代价函数|Ru|2在第k次迭代的亚像素偏移量向量uk处的梯度；
设计以下迭代算法求解连续代价函数|Ru|2的最大值，迭代的具体步骤为：
3.3.1设置初始点u0＝[0,0]T；设置近似Hessian逆矩阵的初值S0＝I2，I2为2×2维单位矩阵，允许误差ε＞0；
3.3.2设置迭代次数k的初始值为0；计算连续代价函数|Ru|2在第k次迭代的亚像素偏移量向量uk处的梯度gk；
3.3.3令dk＝-Skgk，从第k次迭代的亚像素偏移量向量uk出发，沿方向dk进行搜索，求解第k次迭代的步长因子αk，使第k次迭代的步长因子αk满足如下公式：
其中，α为步长因子；
根据求解得到的步长因子αk，令δk＝αkdk；计算第k+1次迭代的亚像素偏移量向量uk+1，uk+1＝uk+δk；
3.3.4计算连续代价函数|Ru|2在第k+1次迭代的亚像素偏移量向量uk+1处的梯度gk+1，若梯度gk+1的范数||gk+1||≤ε，则停止迭代，得到最佳的亚像素偏移量向量uopt＝uk+1，uopt对应的副图像窗口相对于主图像窗口在方向X和方向Y上的最佳的亚像素偏移量分别为dxopt和dyopt；否则，进行下一步操作；
3.3.5令迭代次数k增加1，计算第k+1次迭代的近似Hessian逆矩阵Sk+1，返回步骤3.3.3。
下面结合实验，对本发明的效果作进一步的介绍。
为进一步说明本发明方法相比于其它图像配准方法(如粗配准法、互相关法和最大谱估计法)的优越性，进行以下实验。
1)实验条件
实验采用的数据为TerraSAR-X卫星实测的InSAR图像对，成像区域是澳大利亚的艾尔斯岩石。从TerraSAR-X卫星实测的InSAR图像对中分别取出不同场景的子图像对来进行实验，分别为：岩石一角(区域A)、岩石一侧的陡坡(区域B)和一块平坦区域(区域C)，如图2a中的实现框所示。实验采用800×800个像素大小的子图像对，并将主副图像分解为大小为29×29像素的窗口对，然后对于互相关系数大于0.3的窗口对，采用本发明方法确定这些窗口对之间的亚像素偏移量。
2)实验内容及结果分析
针对TerraSAR-X卫星实测的InSAR图像对中的区域A、区域B以及区域C，分别绘制三个区域的幅度图、配准后的干涉相位图，并用分别用本发明方法、粗配准法、互相关法以及最大谱估计法绘制三个区域的相干直方图。图2a是实验采用的主图像；图2b是实验采用的副图像。
图3a是区域A的幅度图；图3b是区域A配准后的干涉相位图；图3c是本发明方法、粗配准法、互相关法以及最大谱估计法对区域A的相干直方图；从图3c中看出，本发明方法得到的区域A的相关系数高于其它三种方法，表明本发明方法具有更高的配准精度。
图4a是区域B的幅度图；图4b是区域B配准后的干涉相位图；图4c是本发明方法、粗配准法、互相关法以及最大谱估计法对区域B的相干直方图；从图4c中看出，本发明方法得到的区域B的相关系数高于其它三种方法，表明本发明方法具有更高的配准精度。
图5a是区域C的幅度图；图5b是区域C配准后的干涉相位图；图5c是本发明方法、粗配准法、互相关法以及最大谱估计法对区域C的相干直方图；从图5c中看出，本发明方法得到的区域C的相关系数高于其它三种方法，表明本发明方法具有更高的配准精度。
表2是本发明方法、粗配准法、互相关法以及最大谱估计法的对三个区域的相干系数均值和SNR的比较结果。以区域A为例，从表2中可以看出，采用本发明方法、粗配准法、互相关法以及最大谱估计法得到的区域A的相干系数均值分别为0.5755、0.5295、0.5691和0.5465，说明本发明方法获得了更高的相干系数均值，即本发明方法具有更高的配准精度；采用本发明方法、粗配准法、互相关法以及最大谱估计法得到的区域A的SNR(dB)分别为-48.7046、-49.3704、-48.8257和-48.9815，说明本发明方法获得了更高的SNR值，即本发明方法能够获得更清晰的干涉相位图。
表2

Claims (3)

1.一种基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，其特征在于，包括以下步骤：
步骤1，利用干涉SAR获取两幅SAR复图像I1和I2，其中，I1为主图像，I2为副图像；设定副图像I2相对于主图像I1的位置偏移量为(u，v)，并计算主图像I1和副图像I2的二维互相关函数Ru，v，其中，u为副图像I2相对于主图像I1在水平方向上的偏移量，v为副图像I2相对于主图像I1在垂直方向上的偏移量；
步骤2，从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口定义副图像窗口相对于主图像窗口的亚像素偏移量为(dx，dy)，通过对副图像窗口上的像素点(x′+dx，y′+dy)周围的像素点(m，n)进行插值，计算得到副图像窗口在像素点(x′+dx，y′+dy)的像素值根据副图像窗口在像素点(x′+dx，y′+dy)的像素值计算得到主图像窗口和副图像窗口的二维互相关函数构造以亚像素偏移量(dx，dy)为变量的连续代价函数为其中，(x′，y′)为主图像窗口上给定像素点的坐标；
步骤3，优化连续代价函数并求解连续代价函数的最大值，该最大值对应的亚像素偏移量为副图像窗口相对于主图像窗口在方向X和方向Y上的最佳的亚像素偏移量；
其中，所述步骤3的具体子步骤为：
3.1将互相关表示Cp，q表示为：
<mrow> <msub> <mi>C</mi> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mo>=</mo> <msub> <mover> <mi>C</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mo>+</mo> <mi>j</mi> <msub> <mover> <mi>C</mi> <mo>~</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> </mrow>
其中，和分别为互相关表示Cp，q的实部和虚部；
将步骤2.4的求解连续代价函数的最大值的问题简写为：
<mrow> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </munder> <mo>|</mo> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>=</mo> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </munder> <mfrac> <mrow> <msup> <mrow> <mo>&amp;lsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mover> <mi>C</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mover> <mi>C</mi> <mo>~</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </mfrac> </mrow>
根据上式，令令
3.2分别求连续代价函数对方向X上的亚像素偏移量dx和方向Y上的亚像素偏移量dy的一阶偏导，表达式如下：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <mo>|</mo> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </mfrac> <mo>{</mo> <mn>2</mn> <msub> <mover> <mi>R</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>&amp;lsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mover> <mi>C</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> </mrow> </mfrac> <mo>&amp;rsqb;</mo> <mo>+</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mover> <mi>R</mi> <mo>~</mo> </mover> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>&amp;lsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mover> <mi>C</mi> <mo>~</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> </mrow> </mfrac> <mo>&amp;rsqb;</mo> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mfrac> <mrow> <mo>&amp;part;</mo> <mo>|</mo> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </mfrac> <mo>{</mo> <mn>2</mn> <msub> <mover> <mi>R</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>&amp;lsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mover> <mi>C</mi> <mo>&amp;OverBar;</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mo>&amp;rsqb;</mo> <mo>+</mo> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>2</mn> <msub> <mover> <mi>R</mi> <mo>~</mo> </mover> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>&amp;lsqb;</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mover> <mi>C</mi> <mo>~</mo> </mover> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mfrac> <mrow> <mo>&amp;part;</mo> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <mo>&amp;part;</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </mfrac> <mo>&amp;rsqb;</mo> <mo>}</mo> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中|·|表示取模；
3.3采用拟牛顿方法中的BFGS算法来优化步骤3.1的求解连续代价函数的最大值的问题；定义u为亚像素偏移量向量，u＝(dx，dy)T，将连续代价函数简写为|Ru|2；设计迭代公式为：
<mrow> <msub> <mi>S</mi> <mrow> <mi>k</mi> <mo>+</mo> <mn>1</mn> </mrow> </msub> <mo>=</mo> <msub> <mi>S</mi> <mi>k</mi> </msub> <mo>+</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <mfrac> <mrow> <msubsup> <mi>&amp;gamma;</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>S</mi> <mi>k</mi> </msub> <msub> <mi>&amp;gamma;</mi> <mi>k</mi> </msub> </mrow> <mrow> <msubsup> <mi>&amp;gamma;</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>&amp;delta;</mi> <mi>k</mi> </msub> </mrow> </mfrac> <mo>)</mo> </mrow> <mfrac> <mrow> <msub> <mi>&amp;delta;</mi> <mi>k</mi> </msub> <msubsup> <mi>&amp;delta;</mi> <mi>k</mi> <mi>T</mi> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;gamma;</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>&amp;delta;</mi> <mi>k</mi> </msub> </mrow> </mfrac> <mo>-</mo> <mfrac> <mrow> <msub> <mi>&amp;delta;</mi> <mi>k</mi> </msub> <msubsup> <mi>&amp;gamma;</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>S</mi> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>S</mi> <mi>k</mi> </msub> <msub> <mi>&amp;gamma;</mi> <mi>k</mi> </msub> <msubsup> <mi>&amp;delta;</mi> <mi>k</mi> <mi>T</mi> </msubsup> </mrow> <mrow> <msubsup> <mi>&amp;gamma;</mi> <mi>k</mi> <mi>T</mi> </msubsup> <msub> <mi>&amp;delta;</mi> <mi>k</mi> </msub> </mrow> </mfrac> </mrow>
其中，k表示迭代次数，Sk表示第k次迭代的近似Hessian逆矩阵；δk＝uk+1-uk，uk为第k次迭代的亚像素偏移量向量；γk＝gk+1-gk，为连续代价函数|Ru|2在第k次迭代的亚像素偏移量向量uk处的梯度；
设计以下迭代算法求解连续代价函数|Ru|2的最大值，迭代的具体步骤为：
3.3.1设置初始点u0＝[0，0]T；设置近似Hessian逆矩阵的初值S0＝I2，I2为2×2维单位矩阵，允许误差ε＞0；
3.3.2设置迭代次数k的初始值为0；计算连续代价函数|Ru|2在第k次迭代的亚像素偏移量向量uk处的梯度gk；
3.3.3令dk＝-Skgk，从第k次迭代的亚像素偏移量向量uk出发，沿方向dk进行搜索，求解第k次迭代的步长因子αk，使第k次迭代的步长因子αk满足如下公式：
<mrow> <mo>|</mo> <msub> <mi>R</mi> <mrow> <msub> <mi>u</mi> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>&amp;alpha;</mi> <mi>k</mi> </msub> <msub> <mi>d</mi> <mi>k</mi> </msub> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>=</mo> <munder> <mi>max</mi> <mrow> <mi>&amp;alpha;</mi> <mo>&gt;</mo> <mn>0</mn> </mrow> </munder> <mo>|</mo> <msub> <mi>R</mi> <mrow> <msub> <mi>u</mi> <mi>k</mi> </msub> <mo>+</mo> <msub> <mi>&amp;alpha;d</mi> <mi>k</mi> </msub> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow>
其中，α为步长因子；
根据求解得到的步长因子αk，令δk＝αkdk；计算第k+1次迭代的亚像素偏移量向量uk+1，uk+1＝uk+δk；
3.3.4计算连续代价函数|Ru|2在第k+1次迭代的亚像素偏移量向量uk+1处的梯度gk+1，若梯度gk+1的范数||gk+1||≤ε，则停止迭代，得到最佳的亚像素偏移量向量uopt＝uk+1，uopt对应的副图像窗口相对于主图像窗口在方向X和方向Y上的亚像素偏移量分别为dxopt和dyopt；否则，进行下一步操作；
3.3.5令迭代次数k增加1，计算第k+1次迭代的近似Hessian逆矩阵Sk+1，返回步骤3.3.3。
2.如权利要求1所述的基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，其特征在于，所述步骤1的具体子步骤为：
1.1利用干涉SAR获取两幅SAR复图像I1和I2，其中，I1为主图像，I2为副图像；设定像素点(x，y)在主图像I1的像素值为I1(x，y)，像素点(x，y)在副图像I2的像素值为I2(x，y)，I1(x，y)和I2(x，y)的表达式分别为：
其中，|·|表示取模；|I1(x，y)|和|I2(x，y)|分别为像素点(x，y)在主图像I1和副图像I2的幅度；和分别为像素点(x，y)在主图像I1和副图像I2的相位；x表示像素点(x，y)在水平方向上的分量，y表示像素点(x，y)在垂直方向上的分量；
1.2计算主图像I1和副图像I2的二维互相关函数R为：
<mrow> <mi>R</mi> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <munder> <mo>&amp;Sigma;</mo> <mi>x</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>y</mi> </munder> <msup> <msub> <mi>I</mi> <mn>1</mn> </msub> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <msub> <mi>I</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <mrow> <msqrt> <mrow> <munder> <mo>&amp;Sigma;</mo> <mi>x</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>y</mi> </munder> <mo>|</mo> <msub> <mi>I</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> </msqrt> <msqrt> <mrow> <munder> <mo>&amp;Sigma;</mo> <mi>x</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>y</mi> </munder> <mo>|</mo> <msub> <mi>I</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> </msqrt> </mrow> </mfrac> </mrow>
其中，上标*表示共轭；
1.3设定P1和P2分别为主图像I1和副图像I2的能量函数，P1和P2的表达式分别为：
<mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <mo>=</mo> <munder> <mo>&amp;Sigma;</mo> <mi>x</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>y</mi> </munder> <mo>|</mo> <msub> <mi>I</mi> <mn>1</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow>
<mrow> <msub> <mi>P</mi> <mn>2</mn> </msub> <mo>=</mo> <munder> <mo>&amp;Sigma;</mo> <mi>x</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>y</mi> </munder> <mo>|</mo> <msub> <mi>I</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow>
设定副图像I2相对于主图像I1的位置偏移量为(u，v)，则副图像I2相对于主图像I1的位置偏移量为(u，v)时，主图像I1和副图像I2的二维互相关函数Ru，v为：
<mrow> <msub> <mi>R</mi> <mrow> <mi>u</mi> <mo>,</mo> <mi>v</mi> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <munder> <mo>&amp;Sigma;</mo> <mi>x</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>y</mi> </munder> <msup> <msub> <mi>I</mi> <mn>1</mn> </msub> <mo>*</mo> </msup> <mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> </mrow> <msub> <mi>I</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>+</mo> <mi>u</mi> <mo>,</mo> <mi>y</mi> <mo>+</mo> <mi>v</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <msqrt> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>u</mi> <mo>,</mo> <mi>v</mi> <mo>)</mo> </mrow> </mrow> </msqrt> </mfrac> </mrow>
其中，P2(u，v)表示副图像I2相对于主图像I1的位置偏移量为(u，v)时，副图像I2的能量函数。
3.如权利要求1所述的基于解析搜索亚像素偏移量的干涉SAR图像精配准方法，其特征在于，所述步骤2的具体子步骤为：
2.1从主图像I1和副图像I2中分别提取主图像窗口和副图像窗口针对亚像素级的图像配准，将副图像窗口相对于主图像窗口的位置偏移的整数部分移除，则副图像窗口相对于主图像窗口仅存在小于一个像素的偏移量，称之为亚像素偏移量，用dx和dy分别表示方向X和方向Y上的亚像素偏移量，|dx|＜1，|dy|＜1，|·|表示取模；设定主图像窗口上给定像素点的坐标为(x′，y′)，则在副图像窗口上与其精确配准的像素点的坐标为(x′+dx，y′+dy)；通过对副图像窗口上的像素点(x′+dx，y′+dy)周围的像素点(m，n)进行插值，计算副图像窗口在像素点(x′+dx，y′+dy)的像素值为复数，其表达式为：
<mrow> <msub> <mover> <mi>I</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> </msub> <mrow> <mo>(</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mi>l</mi> </mrow> <mrow> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <mi>l</mi> </mrow> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mi>l</mi> </mrow> <mrow> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <mi>l</mi> </mrow> </munderover> <msub> <mover> <mi>I</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>-</mo> <mi>m</mi> <mo>,</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>-</mo> <mi>n</mi> <mo>)</mo> </mrow> </mrow>
其中，x′，y′，m，n∈N，N表示正整数，dx，dy∈R，R表示实数；是副图像窗口在像素点(m，n)的像素值，为复数，m表示像素点(m，n)在水平方向上的分量，n表示像素点(x，y)在垂直方向上的分量；φ(x′，y′)是插值函数，l表示插值函数长度的一半；
将更新为副图像I2在像素点(x′，y′)的像素值，即将副图像窗口逐点精确配准到主图像窗口
2.2将副图像I2相对于主图像I1的位置偏移量为(u，v)时，副图像I2的能量函数P2(u，v)设定为常数，则副图像窗口相对于主图像窗口的亚像素偏移量为(dx，dy)时，主图像窗口和副图像窗口的二维互相关函数为：
<mrow> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <munder> <mo>&amp;Sigma;</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> </munder> <munder> <mo>&amp;Sigma;</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> </munder> <msup> <msub> <mover> <mi>I</mi> <mo>&amp;OverBar;</mo> </mover> <mn>1</mn> </msub> <mo>*</mo> </msup> <mrow> <mo>(</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <msub> <mover> <mi>I</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> </msub> <mrow> <mo>(</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <msqrt> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </msqrt> </mfrac> </mrow>
其中，是主图像窗口在像素点(x′，y′)的像素值，为复数；
2.3将副图像窗口在像素点(x′+dx，y′+dy)的像素值的表达式，代入主图像窗口和副图像窗口的二维互相关函数的表达式为：
<mrow> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <munder> <mo>&amp;Sigma;</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> </munder> <munder> <mo>&amp;Sigma;</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> </munder> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>m</mi> <mo>=</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mi>l</mi> </mrow> <mrow> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <mi>l</mi> </mrow> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mo>=</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>-</mo> <mi>l</mi> </mrow> <mrow> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <mi>l</mi> </mrow> </munderover> <mo>&amp;lsqb;</mo> <msup> <msub> <mover> <mi>I</mi> <mo>&amp;OverBar;</mo> </mover> <mn>1</mn> </msub> <mo>*</mo> </msup> <mrow> <mo>(</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>,</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>)</mo> </mrow> <msub> <mover> <mi>I</mi> <mo>&amp;OverBar;</mo> </mover> <mn>2</mn> </msub> <mrow> <mo>(</mo> <mi>m</mi> <mo>,</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <msup> <mi>x</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>-</mo> <mi>m</mi> <mo>,</mo> <msup> <mi>y</mi> <mo>&amp;prime;</mo> </msup> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>-</mo> <mi>n</mi> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <msqrt> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </msqrt> </mfrac> </mrow>
令x′-m＝p，y′-n＝q，将上式调整为：
<mrow> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mi>C</mi> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <mo>|</mo> </mrow> <msqrt> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </msqrt> </mfrac> </mrow>
其中，Cp，q为副图像窗口相对于主图像窗口存在整数倍像素偏移时的互相关表示，简称为互相关表示，互相关表示Cp，q为复数，其表达式为：
其中，算符*表示卷积，上标*表示共轭，和分别表示主图像窗口和副图像窗口的傅里叶变换，μ和κ分别表示p和q对应的傅里叶变换频率，表示逆傅里叶变换；
2.4根据步骤2.3，构造以亚像素偏移量为变量的连续代价函数为并将求解连续代价函数的最大值的问题表示为：
<mrow> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </munder> <mo>|</mo> <msub> <mi>R</mi> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </msub> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>=</mo> <munder> <mrow> <mi>m</mi> <mi>a</mi> <mi>x</mi> </mrow> <mrow> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> </mrow> </munder> <mfrac> <mrow> <mo>|</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>p</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mo>-</mo> <mi>l</mi> </mrow> <mi>l</mi> </munderover> <msub> <mi>C</mi> <mrow> <mi>p</mi> <mo>,</mo> <mi>q</mi> </mrow> </msub> <mi>&amp;phi;</mi> <mrow> <mo>(</mo> <mi>p</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>x</mi> </msub> <mo>,</mo> <mi>q</mi> <mo>+</mo> <msub> <mi>d</mi> <mi>y</mi> </msub> <mo>)</mo> </mrow> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>P</mi> <mn>1</mn> </msub> <msub> <mi>P</mi> <mn>2</mn> </msub> </mrow> </mfrac> </mrow>
其中，|·|表示取模。
CN201510364208.6A 2015-06-26 2015-06-26 基于解析搜索亚像素偏移量的干涉sar图像精配准方法 Active CN104933673B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510364208.6A CN104933673B (zh) 2015-06-26 2015-06-26 基于解析搜索亚像素偏移量的干涉sar图像精配准方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510364208.6A CN104933673B (zh) 2015-06-26 2015-06-26 基于解析搜索亚像素偏移量的干涉sar图像精配准方法

Publications (2)

Publication Number Publication Date
CN104933673A CN104933673A (zh) 2015-09-23
CN104933673B true CN104933673B (zh) 2018-04-06

Family

ID=54120828

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510364208.6A Active CN104933673B (zh) 2015-06-26 2015-06-26 基于解析搜索亚像素偏移量的干涉sar图像精配准方法

Country Status (1)

Country Link
CN (1) CN104933673B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108022259A (zh) * 2016-11-04 2018-05-11 北京建筑大学 干涉sar复图像配准方法和系统
CN106526594B (zh) * 2016-12-27 2019-02-15 中国科学院电子学研究所 一种ati-sar海面单视复图像配准的方法
CN106990412B (zh) * 2017-05-11 2019-06-11 哈尔滨工业大学 一种干涉合成孔径激光雷达系统的图像配准方法
CN107909606A (zh) * 2017-12-25 2018-04-13 南京市测绘勘察研究院股份有限公司 一种sar图像配准联系点粗差剔除方法
CN109146941A (zh) * 2018-06-04 2019-01-04 成都通甲优博科技有限责任公司 一种基于网格区域划分的深度图像优化方法及系统
CN109100718B (zh) * 2018-07-10 2019-05-28 中国人民解放军国防科技大学 基于贝叶斯学习的稀疏孔径isar自聚焦与横向定标方法
CN109633639B (zh) * 2018-11-20 2020-07-31 上海无线电设备研究所 Topsar干涉数据的高精度快速配准方法
CN110264490B (zh) * 2019-08-15 2019-12-10 成都新西旺自动化科技有限公司 一种应用于机器视觉系统中的亚像素精度边缘提取方法
CN115423848B (zh) * 2022-11-07 2023-02-10 江苏省水利科学研究院 一种识别及去除像素偏移量追踪监测结果异常的方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101126809A (zh) * 2007-09-20 2008-02-20 西安电子科技大学 基于相关加权的干涉合成孔径雷达干涉相位估计方法
CN101706961A (zh) * 2009-11-10 2010-05-12 北京航空航天大学 一种图像的配准方法及装置
CN102955157A (zh) * 2011-08-26 2013-03-06 中国科学院空间科学与应用研究中心 一种用于干涉合成孔径雷达图像精配准的快速相干系数法
CN103325105A (zh) * 2013-02-20 2013-09-25 中国科学院电子学研究所 一种高精度合成孔径雷达图像自动配准方法及设备

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7420592B2 (en) * 2004-06-17 2008-09-02 The Boeing Company Image shifting apparatus for enhanced image resolution

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101126809A (zh) * 2007-09-20 2008-02-20 西安电子科技大学 基于相关加权的干涉合成孔径雷达干涉相位估计方法
CN101706961A (zh) * 2009-11-10 2010-05-12 北京航空航天大学 一种图像的配准方法及装置
CN102955157A (zh) * 2011-08-26 2013-03-06 中国科学院空间科学与应用研究中心 一种用于干涉合成孔径雷达图像精配准的快速相干系数法
CN103325105A (zh) * 2013-02-20 2013-09-25 中国科学院电子学研究所 一种高精度合成孔径雷达图像自动配准方法及设备

Also Published As

Publication number Publication date
CN104933673A (zh) 2015-09-23

Similar Documents

Publication Publication Date Title
CN104933673B (zh) 基于解析搜索亚像素偏移量的干涉sar图像精配准方法
Ye et al. Fast and robust matching for multimodal remote sensing image registration
US8355564B2 (en) Corresponding point searching method and three-dimensional position measuring method
CN105205858B (zh) 一种基于单个深度视觉传感器的室内场景三维重建方法
CN103455797B (zh) 航拍视频中运动小目标的检测与跟踪方法
CN106373087B (zh) 一种改进初始估计的图像超分辨率重建方法
CN105427298B (zh) 基于各向异性梯度尺度空间的遥感图像配准方法
CN106023303B (zh) 一种基于轮廓有效性提高三维重建点云稠密程度的方法
KR101475382B1 (ko) 광학적 3차원 측량의 자기 적응 윈도우 푸리에 위상추출방법
CN105631872B (zh) 基于多特征点的遥感图像配准方法
CN103034982B (zh) 一种基于变焦视频序列的图像超分辨率重建方法
CN107301664A (zh) 基于相似性测度函数的改进局部立体匹配方法
CN106875443B (zh) 基于灰度约束的三维数字散斑的整像素搜索方法及装置
CN110223330A (zh) 一种可见光和红外图像的配准方法及系统
CN107993258A (zh) 一种图像配准方法及装置
CN105654423B (zh) 基于区域的遥感图像配准方法
CN104933678A (zh) 一种基于像素强度的图像超分辨率重建方法
CN104517317A (zh) 一种车载红外图像三维重建方法
CN103400393B (zh) 一种图像匹配方法和系统
CN106709870A (zh) 一种近景影像直线段匹配方法
Ravanbakhsh et al. A comparative study of DEM registration approaches
CN103310482B (zh) 一种三维重建方法及系统
CN103914807B (zh) 一种缩放尺度补偿的非局部性图像超分辨率方法及系统
CN104252701B (zh) 校正视差图的方法和系统
CN104166977B (zh) 一种图像匹配相似性测度方法及其图像匹配方法

Legal Events

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