CN111091043B - 基于多尺度残余图正则化的sar全图变形场估计方法 - Google Patents

Abstract

本发明提供一种基于多尺度残余图正则化的SAR全图变形场估计方法，首先，利用SIFT算法提取的大量局部特征点作为候选控制点，设计基于相位相关和梯度求精的控制点生成策略，以生成稠密且分布相对均匀的控制点对；其次，在变形场估计阶段，设计基于多尺度残余图正则化一致性点漂移的畸变模型，以实现全图像的点对点亚像素配准。本发明在一致性点漂移算法的基础上引入残余图正则化项，解决了一致性点漂移算法中容易陷入局部极小值的缺点，并且能够快速收敛至全局最优解。

Description

基于多尺度残余图正则化的SAR全图变形场估计方法

技术领域

本发明涉及基于参考图的SAR图像高精度(亚像素)全图变形场估计(也称高精度稠密配准或者全图逐点配准)技术。该技术是遥感应用(图像变化检测、高精度数据融合)领域的关键技术。特别涉及一种稠密控制点对产生策略(提取，高精度匹配)及基于多尺度残余图正则化一致性漂移的变形场参数估计。

背景技术

随着遥感观测技术的发展，高分辨率、多平台(卫星、近地空间、机载)、多源(SAR、光学、高光谱)遥感图像越来越丰富，多源图像数据时空高精度配准、多源时空数据融合与变化检测等等相关技术变得越来越重要。要实现多传感器、多平台等数据的融合，关键技术是多源图像的高精度(尤其是亚像素)点点配准技术，其实际是全图像变形场的估计问题。

高分辨率合成孔径雷达(以下简称为“高分SAR”)能够提供全天候、全天时条件下、具有高分辨能力的、详细的地面测绘资料和图像，SAR成为一种极为重要的遥感观测传感器。目前SAR发展向多平台(星载、机载)、多极化(HH、VV、HV)、多视角、多波段(X波段、C波段、ku波段)、高分辨率(亚米级别，1m,0.5m或更高)发展。不同波段的SAR数据反映了地物不同的散射特性，具有互补性；因此融合多波段SAR图像，可以更好地对目标进行观测。高精度SAR逐点匹配(变形场的估计)是实现多时相SAR融合的关键与先决条件。

SAR图像变形场估计技术难点在于：①SAR图像受斑点噪声的影响；②SAR属于斜距成像，转化成地距图像时，局部畸变相当大，尤其在地形起伏大的区域，变形更加复杂，而且这种全图变形场是随着空间位置的变化而变化的。没有数学模型能够解析地描述SAR的局部形变。因此要实现SAR的点对点的配准，其实际是变形场的估计，目前很少有研究集中在SAR变形场的参数估计方面。

发明内容

为了克服现有全图变形场估计(点对点的配准)技术的不足，本发明提出了一种基于多尺度残余图能量最小正则化的全图SAR变形场估计算法。

本发明设计思路为：首先，利用SIFT算法提取的大量局部特征点作为候选控制点，设计基于相位相关和梯度求精的控制点生成策略，以生成稠密且分布相对均匀的控制点对；其次，在变形场估计阶段，设计基于多尺度残余图正则化一致性点漂移的畸变模型，以实现全图像的点对点亚像素配准。

为实现上述目的，本发明提供了基于多尺度残余图正则化的SAR全图变形场估计方法，用于SAR全图像变形场的估计，包括以下步骤：

S100，利用SIFT算法分别对参考图像和输入图像提取局部特征点作为候选控制点；

S200，对候选控制点采用基于相位相关和梯度求精的方法获得稠密且分布均匀的控制点对，本步骤进一步包括：

S210，对候选控制点采用基于相位相关的整像素匹配得到整像素坐标和粗亚像素坐标；

S220，以整像素和粗亚像素坐标为初始值，进行基于梯度求精的亚像素匹配；

S300，基于多尺度残余图正则化一致性点漂移对稠密控制点对进行变形场参数估计，本步骤进一步包括：

S310，基于稠密控制点对，采用残余图正则化一致性点漂移算法估计畸变模型参数；

S320，基于多尺度残余图正则化策略，完成由粗到精的变形场参数估计。

进一步的，步骤S210的具体实现方式如下，

对于每一对特征点

和

其中

表示输入图像中的第i个特征点，用

表示输入图像中以

为中心截取的图像块；

表示参考图像中的第j个特征点，用

表示参考图像中以

为中心截取的图像块；其中

由

经过平移得到，

P(u,v)和Q(u,v)分别是

和

的频域形式，c和r分别表示图像块的宽度和高度，为得到平移参数，计算它们的互功率谱：

式(1)中，(x₀+Δx,y₀+Δy)表示待求的平移参数，(x₀,y₀)表示平移参数的整数部分，即整像素坐标，(Δx,Δy)表示平移参数的小数部分，即粗亚像素坐标，j表示虚数单位，u和v表示频域分量，(*)代表复共轭，||表示该项的幅值谱，频域与时域之间的转换均采用FFT进行计算，相位相关的方法就是求取式(1)的傅里叶反变换，找到峰值所在位置来确定平移参数，整像素坐标(x₀,y₀)和粗亚像素坐标(Δx,Δy))。

进一步的，步骤S220的具体实现方式如下，

针对相位相关法中得到的整像素坐标(x₀,y₀)和粗亚像素坐标(Δx,Δy)对亚像素坐标进行细化，首先对q(x+x₀+Δx,y+y₀+Δy)在(x+x₀,y+y₀)处计算一阶泰勒展开，去掉高阶小量后的一阶泰勒公式见式(2)：

q(x+x₀+Δx,y+y₀+Δy)＝q(x+x₀,y+y₀)+Δxq_x+Δyq_y (2)

式(2)中，q_x表示q(x+x₀,y+y₀)对x的偏导数，q_y表示q(x+x₀,y+y₀)对y的偏导数；

为了使控制点对之间的匹配误差最小，转化为最优化问题求解：

式(3)中，H(Δx,Δy)表示平方和误差函数，x'＝x+x₀+Δx，y'＝y+y₀+Δy；根据极值点特性求得H(Δx,Δy)取极小值时的精亚像素坐标(Δx,Δy)，结合相位相关法计算而来的整像素坐标(x₀,y₀)，更新参考图像中的特征点坐标，即得到输入图像与参考图像之间的亚像素级控制点对。

进一步的，步骤S310的具体实现方式如下，

所述残余图正则化一致性点漂移算法是在考虑一致性运动约束的情况下，将点集配准问题考虑为混合高斯模型的估计问题，采用EM算法优化迭代求解最大似然函数，具体为：

通过最大化目标函数f₂(W^s,σ²)求解畸变模型参数，

式(4)中，s∈{0,1,...,S-1}表示尺度参数，S表示多尺度的最优尺度参数，计算时需对原始图像及控制点对进行下采样获得s尺度下的图像和控制点对；f₁(W^s,σ²)表示负对数似然函数；W^s表示s尺度下的畸变模型参数，σ表示单高斯模型的标准差；

表示上一次迭代中，以s尺度下参考图像中的第n个控制点坐标为先验条件的第m个高斯分量的后验概率，且

表示s尺度下参考图像中第n个控制点坐标，

表示s尺度下输入图像中第m个控制点坐标；

表示在

变形后的坐标；M和N分别表示输入图像与参考图像的控制点数；

式(5)中，q_n表示参考图像中第n个控制点坐标；Γ(p_m,W^old)表示输入图像中的控制点p_m经过非刚性变换后的新位置；W^old表示上一次迭代求解的非刚性变换参数；σ^old表示上一次迭代求解的单高斯模型的标准差；ω为权重因子。

其中，

为正则化项，λ表示正则化参数，

表示s尺度下的参考图像，

表示变形后的输入图像；目标函数模型表示为f₂(W^s,σ²)，通过求解目标函数的极大值，得到s尺度下畸变模型参数W^s。

进一步的，步骤S320的具体实现方式如下，

多尺度变形场方程由式(6)表示：

其中：Γ^s表示s尺度下的变形场方程，Γ⁰表示原始图像的变形场方程，P^s表示待变形图像或控制点坐标，W^s表示s尺度下的变形场参数，G^s是一个核矩阵，其中

和

分别是s尺度下待变形图像中第i个和第j个控制点坐标，

表示核矩阵中第i行，第j列的值，

为变形后的图像或坐标；

通过变形场方程完成图像的局部畸变：

式(7)中，

即为s尺度下待变形的输入图像，

即为变形结果；这种由粗到精的方法降低了变形图像与参考图像之间的均方根误差，结合(4)式中用到的残余图正则化，进一步降低变形误差；

最后使用阈值η限制每一次迭代正则化项的变化程度，即将

与η进行比较，其中

若满足

则继续进行迭代，否则停止迭代并输出变形结果。

进一步的，多尺度的最优尺度参数S为4，最优正则化参数λ是9.5。

与现有技术相比，本发明具有以下优点：

(1)本发明利用SIFT算法提取的大量局部特征点作为候选控制点，该控制点产生策略优势相对于常规结构特征点而言更加稠密，而且分布相对均匀，这些点具有局部特征丰富的特点；使用基于相位相关和梯度求精的块匹配算法既能快速完成匹配，也能提高控制点对匹配的精度；这样提高匹配精度的同时，也能获得大量的控制点对，为后续变形场的精确估计提供了参考依据。

(2)在一致性点漂移算法的基础上引入残余图正则化项，解决了一致性点漂移算法中容易陷入局部极小值的缺点，并且能够快速收敛至全局最优解。

(3)在残余图正则化一致性点漂移算法的基础上引入多尺度的策略，在大尺度(低分辨率)时，对形变进行较大尺度的纠正；在小尺度(高分辨率)较小时，对微小形变进行纠正，整个多尺度的变形场估计过程是一种由粗到精的过程，比单一尺度下的估计更为准确，其次，利用得到的残余图对区域变形场进行校正，提高了定位精度。

附图说明

图1为本发明的流程图；

图2为基于相位相关与梯度求精的块匹配算法流程图；

图3为基于相位相关与梯度求精的块匹配算法结果示意图；

图4为基于多尺度残余图正则化一致性点漂移算法流程图；

图5为基于多尺度残余图正则化的一致性点漂移结果示意图，(a)是参考图像，(b)是要估计变形的图像，(c)是常规算法的残差图，(d)本发明算法的残差图。I、II是高分辨率分析窗口的原图与残差图像；

图6为基于多尺度残余图正则化的一致性点漂移向量场与变形网格，(a)是参考图像，(b)是要估计变形的图像，(c)局部变形纠正后图像，(d)迭代收敛曲线，(e)变形场(f)局部窗口放大后的变形网格；

图7为全图变形场估计后控制点对误差方差，(a)x方向，(b)y方向；

图8为不同方法在笛卡尔坐标系中坐标误差比较，(a)TPS-RPM，(b)CPD，(c)本发明方法。

具体实施方式

下面结合附图对本发明的技术方案作进一步的描述。

如图1所示，本发明采用的技术方案包括以下几个关键部分与技术：

第一部分：稠密并且分布均匀的控制点对产生与高精度匹配策略。

本发明在稠密控制点对产生与匹配上设计了一个新方法。首先，我们利用大量的SIFT点作为候选控制点对，所以利用SIFT算子提取大量的候选点对。这些特征点比较稠密，而且相对分布均匀，这些点具有局部特征丰富的特点。

所述基于SIFT的特征点提取可概括为：

(1)构建尺度空间(2)极值点检测(3)关键点定位(4)主方向确定(5)特征描述符的构建(6)特征抽取与匹配。

在控制点对匹配上，块匹配新算法使用相位相关的方法计算特征匹配整像素和粗亚像素坐标，在计算出整像素和粗亚像素坐标后，使用基于梯度求精的方法可以计算出亚像素的匹配精度。该方法流程如图2所示。

本发明设计的基于相位相关与梯度求精的块匹配具体如下：

对于每一对特征点

和

其中

表示输入图像中的第i个特征点，用

表示输入图像中以

为中心截取的图像块；

表示参考图像中的第j个特征点，用

表示参考图像中以

为中心截取的图像块。其中

可由

经过平移得到，

P(u,v)和Q(u,v)分别是

和

的频域形式，c和r分别表示图像块的宽度和高度。为得到平移参数，计算它们的互功率谱：

式(1)中，(x₀+Δx,y₀+Δy)表示待求的平移参数，(x₀,y₀)表示平移参数的整数部分，即整像素坐标，(Δx,Δy)表示平移参数的小数部分，即粗亚像素坐标，j表示虚数单位，u和v表示频域分量，(*)代表复共轭，||表示该项的幅值谱。本专利中频域与时域之间的转换均采用FFT进行计算，相位相关的方法就是求取式(1)的傅里叶反变换，找到峰值所在位置来确定平移参数。只使用相位相关法得到的平移参数在位移上存在一定误差，为了提高亚像素的匹配精度，本发明引入了一种基于梯度求精的亚像素匹配方法：

针对相位相关法中得到的整像素坐标(x₀,y₀)和粗亚像素坐标(Δx,Δy)对亚像素坐标进行细化。首先对q(x+x₀+Δx,y+y₀+Δy)在(x+x₀,y+y₀)处计算一阶泰勒展开，去掉高阶小量后的一阶泰勒公式见式(2)：

q(x+x₀+Δx,y+y₀+Δy)＝q(x+x₀,y+y₀)+Δxq_x+Δyq_y (2)

式(2)中，q_x表示q(x+x₀,y+y₀)对x的偏导数，q_y表示q(x+x₀,y+y₀)对y的偏导数。

为了使控制点对之间的匹配误差最小，转化为最优化问题求解：

式(3)中，H(Δx,Δy)表示平方和误差函数，x'＝x+x₀+Δx，y'＝y+y₀+Δy。根据极值点特性求得H(Δx,Δy)取极小值时的精亚像素坐标(Δx,Δy)。结合相位相关法计算而来的整像素坐标(x₀,y₀)，更新参考图像中的特征点坐标，即可得到输入图像与参考图像之间的亚像素级控制点对。

基于梯度求精的亚像素匹配算法具有较低的算法复杂度，易于实现，同时匹配结果精确，控制点对的精度可达亚像素级。本发明将两种方法结合起来。首先利用相位相关法得到相对粗糙的像素匹配位置即整像素坐标，然后利用整像素坐标和粗亚像素坐标作为梯度求精法的初值，利用一阶泰勒级数展开减小亚像素级匹配的误差，大大提高了亚像素匹配定位精度。图3表示了基于相位相关与梯度的块匹配配准效果，从图3中可以看出，经过相位相关与梯度求精的配准方法之后，得到经过匹配的特征点均匀分布在图像中，这些均匀分布的特征点也为后续的变形场估计提供了初始化条件。

第二部分：基于多尺度残余图正则化一致性点漂移的变形场参数估计阶段。

这个阶段主要设计了基于多尺度残余图正则化一致性的点漂移算法来完成变形场的估计。通过多尺度方法的引入，在大尺度(低分辨率)时，求解大尺度下的图像及相应的控制点对，对形变进行较大尺度的纠正；在小尺度(高分辨率)时，求解小尺度下的图像及控制点对，对微小形变进行纠正，整个多尺度的变形场估计过程是一种由粗到精的过程，这样一个由粗尺度到精尺度的图像变形估计过程，降低了陷入局部极小值的可能性，提高了变形场参数精确度，同时，为了指导算法进行变形场估计时的变形趋势，使用了残余图作为正则化项，提高了算法对于局部形变的鲁棒性。该算法流程图如图4所示。

本发明所设计的基于多尺度残余图正则化一致性点漂移算法具体如下：

一致性点漂移算法在考虑一致性运动约束的情况下，将点集配准问题考虑为混合高斯模型的估计问题，采用EM算法优化迭代求解最大似然函数。EM算法的思想是先初始化参数的值(“旧的”参数值)，然后利用贝叶斯定理计算混合分量的后验概率分布P^old(m|q_n)，即为算法的E步骤；通过最大化负对数似然函数，找到“新”参数值，即为算法的M步骤。EM算法在E步骤和M步骤之间交替迭代进行，直到所求模型收敛。

通过最大化目标函数f₂(W^s,σ²)求解畸变模型参数：

式(4)中，s∈{0,1,...,S-1}表示尺度参数，计算时需对原始图像及控制点对进行下采样获得s尺度下的图像和控制点对；f₁(W^s,σ²)表示负对数似然函数；W^s表示s尺度下的畸变模型参数，σ表示单高斯模型的标准差；

表示上一次迭代中，以s尺度下参考图像中的第n个控制点坐标为先验条件的第m个高斯分量的后验概率，且

表示s尺度下参考图像中第n个控制点坐标，

表示s尺度下输入图像中第m个控制点坐标；

表示在

变形后的坐标；M和N分别表示输入图像与参考图像的控制点数；

式(5)中，q_n表示参考图像中第n个控制点坐标；Γ(p_m,W^old)表示输入图像中的控制点p_m经过非刚性变换后的新位置；W^old表示上一次迭代求解的非刚性变换参数；σ^old表示上一次迭代求解的单高斯模型的标准差；ω为权重因子；

本发明为防止EM算法求解的混合高斯模型陷入局部极小值，对负对数似然函数添加正则化项

因此目标函数模型表示为f₂(W^s,σ²)；λ表示正则化参数；

表示s尺度下的参考图像；

表示变形后的输入图像；通过求解目标函数的极大值，得到s尺度下畸变模型参数W^s；由于该正则化项计算前未知，需要对其进行初始化(置0)。

多尺度变形场方程可由式(6)表示：

式(6)中，Γ^s表示s尺度下的变形场方程，Γ⁰表示原始图像的变形场方程，且不需要进行下采样。P^s表示待变形图像或控制点坐标；W^s表示s尺度下的畸变模型参数，G^s是一个核矩阵，其中

和

分别是s尺度下待变形图像中第i个和第j个控制点坐标，

表示核矩阵中第i行，第j列的值，

为变形后的图像或坐标。

通过变形场方程完成图像的局部畸变：

式(7)中，

即为s尺度下待变形的输入图像，

即为变形结果。这种由粗到精的方法降低了变形图像与参考图像之间的均方根误差，结合(4)式中用到的残余图正则化，进一步降低变形误差。在本发明中，我们使用阈值η限制每一次迭代正则化项的变化程度，即将

与η进行比较，其中

若满足

则继续进行迭代，否则停止迭代并输出变形结果。在本发明方法的实施过程中，我们发现每次迭代后，残差的范数总缩小一个常数因子，从而收敛至零。换句话说，它是参数空间中一种由粗到精的策略。

图5和图6分别表示了基于多尺度残余图正则化的一致性点漂移的变形场估计结果。其中图5是基于多尺度残余图正则化一致性点漂移结果示意图，图(a)是参考图像，图(b)是要估计变形的输入图像，图(c)是常规算法的残余图，图(d)是本发明算法的残余图，矩形框放大显示出变形前后的局部形变，I、II是高分辨率分析窗口的原图与残差图像。图6是多尺度一致性点漂移向量场与变形网格，着重给出了变形场的矢量图与变形网格，通过矢量图与变形网格可以更加直接的看出局部变形的趋势，图(a)是参考图像，图(b)是要估计变形的图像，图(c)是局部变形纠正后图像，图(d)是迭代收敛曲线，图(e)是变形场，图(f)是局部窗口放大后的变形网格。

通过对算法进行大量的实验验证，得到最优的正则化参数与尺度层次，多尺度的最优尺度参数S为4；而最优正则化参数λ是9.5；这两个参数值也是专利申请保护点。

除了定性的算法进行分析外，下面则是定量的对算法进行了验证，主要是计算了经过变形前后的方差与坐标误差的极坐标如图7和图8所示。从下图可以看出，本发明的优势远远超过常规的算法(TPS-RPM，CPD)

从图8中可以看出，基于多尺度一致性点漂移算法的坐标误差主要集中在0.5像素以内，相比于TPS-RPM与CPD算法而言，提出的算法在变形结果的精度上有了很大的提升。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (4)

1.基于多尺度残余图正则化的SAR全图变形场估计方法，其特征在于，包括如下步骤：

S100，利用SIFT算法分别对参考图像和输入图像提取局部特征点作为候选控制点；

S200，对候选控制点采用基于相位相关和梯度求精的方法获得稠密且分布均匀的控制点对，本步骤进一步包括：

S210，对候选控制点采用基于相位相关的整像素匹配得到整像素坐标和粗亚像素坐标；

S220，以整像素和粗亚像素坐标为初始值，进行基于梯度求精的亚像素匹配；

S300，基于多尺度残余图正则化一致性点漂移对稠密控制点对进行变形场参数估计，本步骤进一步包括：

S310，基于稠密控制点对，采用残余图正则化一致性点漂移算法估计畸变模型参数；

步骤S310的具体实现方式如下，

所述残余图正则化一致性点漂移算法是在考虑一致性运动约束的情况下，将点集配准问题考虑为混合高斯模型的估计问题，采用EM算法优化迭代求解最大似然函数，具体为：

通过最大化目标函数f₂(W^s,σ²)求解畸变模型参数，

式(4)中，s∈{0,1,...,S-1}表示尺度参数，S表示多尺度的最优尺度参数，计算时需对原始图像及控制点对进行下采样获得s尺度下的图像和控制点对；f₁(W^s,σ²)表示负对数似然函数；W^s表示s尺度下的畸变模型参数，σ表示单高斯模型的标准差；

表示上一次迭代中，以s尺度下参考图像中的第n个控制点坐标为先验条件的第m个高斯分量的后验概率，且

表示s尺度下参考图像中第n个控制点坐标，

表示s尺度下输入图像中第m个控制点坐标；

表示在

变形后的坐标；M和N分别表示输入图像与参考图像的控制点数；

式(5)中，q_n表示参考图像中第n个控制点坐标；Γ(p_m,W^old)表示输入图像中的控制点p_m经过非刚性变换后的新位置；W^old表示上一次迭代求解的非刚性变换参数；σ^old表示上一次迭代求解的单高斯模型的标准差；ω为权重因子；

其中，

为正则化项，λ表示正则化参数，

表示s尺度下的参考图像，

表示变形后的输入图像；目标函数模型表示为f₂(W^s,σ²)，通过求解目标函数的极大值，得到s尺度下畸变模型参数W^s；S320，基于多尺度残余图正则化策略，完成由粗到精的变形场参数估计；

步骤S320的具体实现方式如下，

多尺度变形场方程由式(6)表示：

其中：Γ^s表示s尺度下的变形场方程，Γ⁰表示原始图像的变形场方程，P^s表示待变形图像或控制点坐标，W^s表示s尺度下的变形场参数，G^s是一个核矩阵，其中

和

分别是s尺度下待变形图像中第a个和第b个控制点坐标，

表示核矩阵中第a行，第b列的值，

为变形后的图像或坐标；

通过变形场方程完成图像的局部畸变：

式(7)中，

即为s尺度下待变形的输入图像，

即为变形结果；最后使用阈值η限制每一次迭代正则化项的变化程度，即将

与η进行比较，其中

若满足

则继续进行迭代，否则停止迭代并输出变形结果。

2.如权利要求1所述的基于多尺度残余图正则化的SAR全图变形场估计方法，其特征在于：步骤S210的具体实现方式如下，

对于每一对特征点

和

其中

表示输入图像中的第i个特征点，用

表示输入图像中以

为中心截取的图像块；

表示参考图像中的第j个特征点，用

表示参考图像中以

为中心截取的图像块；其中

由

经过平移得到，

P(u,v)和Q(u,v)分别是

和

的频域形式，c和r分别表示图像块的宽度和高度，为得到平移参数，计算它们的互功率谱：

式(1)中，(x₀+Δx,y₀+Δy)表示待求的平移参数，(x₀,y₀)表示平移参数的整数部分，即整像素坐标，(Δx,Δy)表示平移参数的小数部分，即粗亚像素坐标，j”表示虚数单位，u和v表示频域分量，*代表复共轭，||表示该项的幅值谱，频域与时域之间的转换均采用FFT进行计算，相位相关的方法就是求取式(1)的傅里叶反变换，找到峰值所在位置来确定平移参数，整像素坐标(x₀,y₀)和粗亚像素坐标(Δx,Δy)。

3.如权利要求2所述的基于多尺度残余图正则化的SAR全图变形场估计方法，其特征在于：步骤S220的具体实现方式如下，

针对相位相关法中得到的整像素坐标(x₀,y₀)和粗亚像素坐标(Δx,Δy)对亚像素坐标进行细化，首先对q(x+x₀+Δx,y+y₀+Δy)在(x+x₀,y+y₀)处计算一阶泰勒展开，去掉高阶小量后的一阶泰勒公式见式(2)：

q(x+x₀+Δx,y+y₀+Δy)＝q(x+x₀,y+y₀)+Δxq_x+Δyq_y (2)

式(2)中，q_x表示q(x+x₀,y+y₀)对x的偏导数，q_y表示q(x+x₀,y+y₀)对y的偏导数；

为了使控制点对之间的匹配误差最小，转化为最优化问题求解：

式(3)中，H(Δx,Δy)表示平方和误差函数，x'＝x+x₀+Δx，y'＝y+y₀+Δy；根据极值点特性求得H(Δx,Δy)取极小值时的精亚像素坐标(Δx,Δy)，结合相位相关法计算而来的整像素坐标(x₀,y₀)，更新参考图像中的特征点坐标，即得到输入图像与参考图像之间的亚像素级控制点对。

4.如权利要求1所述的基于多尺度残余图正则化的SAR全图变形场估计方法，其特征在于：多尺度的最优尺度参数S为4，最优正则化参数λ是9.5。

Images