CN100435167C - 高分辨率sar影像配准处理方法及系统 - Google Patents

Abstract

本发明提供一种高分辨率合成孔径雷达(SAR)影像配准处理方法及系统，该方法包括步骤：(a)创建金字塔影像，进行影像纹理特征提取，得到边缘信息作为匹配输入数据；(b)基于所提取特征进行影像全局粗匹配；(c)提取粗配准后的主、辅图像的特征角点，并进行同名点交叉搜索，获得同名点对；(d)根据同名点对进行自适应区域分块处理。所提供的系统包括：影像数据提取单元，影像预处理单元，影像配准处理单元。采用该方法或系统能较好地达到全局子像素级配准精度，这就为SAR数据在地图制图、变化检测、三维地形重建等应用领域的进一步推广使用提供了基础技术支持。

Description

高分辨率SAR影像配准处理方法及系统

技术领域

本发明涉及影像处理技术领域，具体涉及遥感影像配准处理方法及系统，尤其是合成孔径雷达(SAR)影像的自动配准方法及系统。

背景技术

在遥感影像集成分析中，影像配准是遥感数据分析和处理的基础，是影像融合、动态变化检测、三维地形重建、遥感影像镶嵌等应用技术的关键环节。回顾近十年，遥感影像自动配准技术发展迅速，但有关合成孔径雷达(SAR)影像，尤其是高分辨率SAR影像自动配准的研究工作却正处于起步阶段。至今，已有一些有关SAR影像的自动配准算法被相继提出；2000年，Thepaut等从影像融合角度出发，提出了一种有关ERS和SPOT影像的自动配准算法，然而该算法基于地球形状、地球自转，卫星轨道参数及传感器几何成像方式等大量先验知识，在很多情况下条件难于得到满足；2001年，Dare等提出了一种有关SAR和SPOT影像自动配准模型，使用影像边缘进行特征匹配，然而其处理过程需要人为干预，对大图幅的影像支持不佳；2003年，Chen等提出了一种基于影像MI(Mutual Information)的匹配算法，然而该算法存在处理图幅小，配准精度不高等不足。

发明内容

鉴于此，为解决高分辨率机载、星载SAR影像自动/半自动配准问题。本发明提供一种高分辨率SAR影像配准处理方法，综合采用特征匹配和灰度匹配的方法，实施由粗到精的匹配策略；精匹配中，基于最小二乘影像区域自适应分块算法的提出，不仅使匹配达到了全局子像素级的精度，而且增加了算法的灵活性和稳健性。这就为SAR数据在地图制图、变化检测、三维地形重建等应用领域的进一步推广使用提供了基础技术支持。本发明所提供的高分辨率SAR影像配准处理方法，包括步骤：

(a)提取影像纹理特征，得到边缘信息并作为匹配输入数据；

(b)基于所提取的影像纹理特征进行影像全局粗匹配，得到配准后的主、辅影像；

(c)提取粗配准后主、辅影像的特征角点，并进行同名点交叉搜索，获得同名点对；

(d)根据同名点对进行自适应区域分块处理。

所述同名点为主、辅影像上对应同一地物的特征点；主、辅影像同一地物特征点对可构成同名点对，同名点对有时亦可称为控制点(Control Points or Tie-Points)。

其中步骤(a)包括：创建金字塔影像；采用Canny算法，对主、辅影像进行边缘提取；对边缘图像进行去除噪声处理；根据边缘长度对所述边缘图像进一步处理。

其中步骤(b)中采用快速Hausdorff距离算法进行粗匹配。

其中步骤(c)中采用均匀分布策略提取特征角点；

其中步骤(c)中采用主、辅影像交叉搜索方法，分别以主、辅影像特征角点为参考，根据特征粗匹配的变换关系，在影像搜索区中寻找同名点；然后把所得同名点对赋予相应的索引号；将所得同名点对进行合并处理。

其中步骤(d)中采用基于最小二乘区域自适应分块算法，所述自适应区域分块处理，包括步骤：

(d1)初始区域分块，根据公式(III)获得初始分块子图像同名点对误差方程式(IV)，

x₂＝a₀+a₁x+a₂y    (III)

y₂＝b₀+b₁x+b₂y

$V_x^i=\frac{\∂x}{\∂a_0^i}d{a}_0^i+\frac{\∂x}{\∂a_1^i}d{a}_1^i+\frac{\∂x}{\∂a_2^i}d{a}_2^i-l_x^i$

$V_y^i=\frac{\∂y}{\∂b_0^i}d{b}_0^i+\frac{\∂y}{\∂b_1^i}d{b}_1^i+\frac{\∂y}{\∂b_2^i}d{b}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

其中l_x ⁱ表示观测值x₂ ⁱ同a₀ ⁱ+a₁ ⁱx+a₂ ⁱy计算值的拟合误差，l_y ⁱ表示观测值y₂ ⁱ同b₀ ⁱ+b₁ ⁱx+b₂ ⁱy计算值的拟合误差。

(d2)迭代计算变形参数改正值da₀ ⁱ，da₁ ⁱ，da₂ ⁱ，db₀ ⁱ，db₁ ⁱ，db₂ ⁱ及对应变形参数a₀ ⁱ，a₁ ⁱ，a₂ ⁱ，b₀ ⁱ，b₁ ⁱ，b₂ ⁱ；

(d3)子图像块分裂，回退，增加控制点处理；

(d4)修正同名点对的图像坐标值；

(d5)依次对辅影像中包含的子图像块使用对应变形参数进行几何变换，得到配准后的影像。

根据本发明的一个方案，还提供一种高分辨率SAR影像配准处理系统，包括：

影像数据提取单元，用于提取影像纹理特征数据，得到边缘信息，并将该信息送入到影像预处理单元作为粗匹配输入数据；

影像预处理单元，对已提取影像纹理特征的影像进行全局粗匹配；得到粗配准后的主、辅影像信息，并将相关信息提供给影像配准处理单元处理；

影像配准处理单元，用于提取粗匹配后得到的主、辅影像的特征角点，并进行同名点交叉搜索，获得同名点对，并根据同名点对进行自适应区域分块处理。

其中所述影像数据提取单元处理影像的步骤，包括：创建金字塔影像；采用Canny算法，对主、辅影像进行边缘提取；对边缘图像进行去除噪声处理；根据边缘长度对所述边缘图像进一步处理。

其中在影像预处理单元中采用快速Hausdorff距离算法进行粗匹配。

其中影像配准处理单元中采用均匀分布策略提取特征角点；

其中影像配准处理单元中采用主、辅影像交叉搜索方法，分别以主、辅影像特征角点为参考，根据特征粗匹配的变换关系，在影像搜索区中寻找同名点；然后把所得同名点对赋予相应的索引号；将所得同名点对进行合并处理。

其中影像配准处理单元中采用基于最小二乘区域自适应分块算法，所述自适应区域分块处理，包括步骤：

(d1)初始区域分块，根据公式(III)获得初始分块子图像同名点对误差方程式(IV)，

x₂＝a₀+a₁x+a₂y    (III)

y₂＝b₀+b₁x+b₂y

$V_x^i=\frac{\∂x}{\∂a_0^i}d{a}_0^i+\frac{\∂x}{\∂a_1^i}d{a}_1^i+\frac{\∂x}{\∂a_2^i}d{a}_2^i-l_x^i$

$V_y^i=\frac{\∂y}{\∂b_0^i}d{b}_0^i+\frac{\∂y}{\∂b_1^i}d{b}_1^i+\frac{\∂y}{\∂b_2^i}d{b}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

其中l_x ⁱ表示观测值x₂ ⁱ同a₀ ⁱ+a₁ ⁱx+a₂ ⁱy计算值的拟合误差，l_y ⁱ表示观测值y₂ ⁱ同b₀ ⁱ+b₁ ⁱx+b₂ ⁱy计算值的拟合误差。

(d2)迭代计算变形参数改正值da₀ ⁱ，da₁ ⁱ，da₂ ⁱ，db₀ ⁱ，db₁ ⁱ，db₂ ⁱ及对应变形参数a₀ ⁱ，a₁ ⁱ，a₂ ⁱ，b₀ ⁱ，b₁ ⁱ，b₂ ⁱ；

(d3)子图像块分裂，回退，增加控制点处理；

(d4)修正同名点对的图像坐标值；

(d5)依次对辅影像中包含的子图像块使用对应变形参数进行几何变换，得到配准后的影像。

附图说明

图1为本发明的高分辨率SAR影像配准处理方法的流程图；

图2为根据本发明的具体实施例的处理过程示意图；

图3为本发明的方法中所采用的快速Hausdorff距离粗匹配的程序框图；

图4为本发明的方法中所采用的同名点交叉搜索的程序框图；

图5为本发明的方法中所采用的LS-MMSE算法流程的程序框图。

具体实施方式

现对本发明的高分辨率SAR影像配准处理方法的具体实现予以详细描述。

图2为根据本发明的具体实施例的处理过程示意图，参照图2。

纹理特征提取是为基于Hausdorff Distance全局特征粗匹配提供影像纹理信息，并作为匹配方法的输入数据。为了得到主要边缘信息，抑制噪声点，采用以下方法，其包括：

步骤1：创建金字塔影像的步骤；

步骤2：采用Canny算法，对主、辅影像进行边缘提取的步骤。值得注意的是，创建金字塔影像步骤已经达到了一定的影像平滑效应，因此，在Canny算法执行中，不需另包含Gaussian滤波处理过程；接着使用Hysteresis方法，对输入边缘图像进行下述处理去除噪声，具体步骤为：

①扫描所有边缘点，幅度值大于T₁得到保留，并标为正确点；

②边缘点幅度值小于T₂直接去除；

③扫描幅度值介于T₁，T₂的所有边缘点，如果该点的八邻域已经存在正确点，则亦标为正确点；否则，直接去除。重复对边缘图像进行该运算，直至稳定而中止。

步骤3：对经由Hysteresis处理的边缘图像进行边缘长度检测步骤。设定阈值T，边缘长度大于T的边缘得到保留，反之，直接去除。

图3为本发明的方法中所采用的快速Hausdorff距离粗匹配的程序框图，参照图3。

给定两维特征序列组A和B，Hausdorff距离h^f(B，A)定义如下：

$h^f(B,A)\stackrel{\mathrm{def}}{=}\mathrm{fth}\min \left|\right|b-a\left|\right|---\left(I\right)$

式中，

表示g(x)在序列X中的分位数，f∈(0，1)。譬如，1th代表最值，1/2th代表中值。||b-a||表示点a，b之间的欧氏距离。

本发明采用的基于Hausdorff距离算法简要描述如下：

使用快速算法，计算主影像E₁的距离变换。E₁距离变换的定义如下： $\mathrm{\Δ}(x,y)=\underset{i\∈E_1}{\min }\left|\right|(x,y)-i\left|\right|$

构造存储单元R的堆栈(Stack)。全局特征粗匹配的目的是寻找主、辅影像间的最佳平移变换t[T_x，T_y]，即满足关系式：

h^f(t(E₂)，E₁)≤τ+ε(II)

式中，E₁，E₂表示经由特征提取后的金字塔边缘图像；τ是一个预先设定的阈值，t(E₂)表示对E₂进行t变换， $\mathrm{\ϵ}=\sqrt{2}.$ R表示平移变换的二维空间。对于给定的R，t^l表示T_x、T_y都取极小值T_x ^l、T_y ^l时，即t^l(T_x ^l，T_y ^l)的情况；同理，t^h表示T_x、T_y都取极大值T_x ^h、T_y ^h时，即t^h(T_x ^h，T_y ^h)的情况。对于任意的边缘点P∈E₂，t∈R，变换t(P)必定落于E₁对应的一定矩形区域内。该区域被定义为P的不定区域。因此，对于P∈E₂，每一个点便对应一个不定区域。初始化堆栈(Stack)，使堆栈(Stack)包含一个初始单元R，令R初始等级为零级。值得注意的是，初始单元R对应不定区域应该选得足够大，使得它能包含目标转换函数t，来满足关系式(II)。

寻找最佳平移变换量t。

a)计算待选单元R对应的平移变换差分和最大不定区域

平移变换差分：Δt＝t^h-t^l

不定区域： $\underset{(x_i,y_i)\∈E_2}{\max }(T_{x_i},T_{y_i})$

b)分裂单元R为四等分子单元R_i(i＝1，2，3，4)，R_i对应等级为其父单元等级数加一，R_i分别记为：{(T_x ^l，T_y ^l)，(T_x ^m，T_y ^m)}，{(T_x ^m，T_y ^l)，(T_x ^h，T_y ^m)}；{(T_x ^l，T_y ^m)，(T_x ^m，T_y ^h)}，{(T_x ^m，T_y ^m)，(T_x ^h，T_y ^h)}；

c)判断R是否包含目标平移转换函数。判断工作分两步进行：首先计算在待选单元R内，每一边缘点P∈E₂所对应的不定区域；若区域包含至少一个兴趣点(满足关系式(II)的边缘点)，则把该区域定为合格区域；其次，核查在单元R内，合格区域数占总边缘点的比例；若该比例值大于预先设定的fth值，则待选单元R被记为有效。

d)对每一层有效子单元R_i，根据fth值实施排序，选取最大和次大者压入堆栈，作为进一步分裂的待选单元。

e)迭代步骤a)-d)，直至不定区域的长度收敛于预先设定的τ值。步骤a)-d)的迭代处理，本质上等价于构建一棵搜索二叉树。

下面对特征角点的提取予以描述。

特征角点是影像典型点特征，它在影像相对定向中具有明显、几何位置精确等性质。可选用经典的Forstner算子对主、辅影像进行特征角点提取，该算子通过计算各像素的Robert’s梯度和像素(c，r)为中心的一个窗口(如5×5)的灰度协方差矩阵，在影像中寻找具有尽可能小而接近圆的误差椭圆的点作为特征点，其步骤包括：

(1)计算各像素的Robert’s梯度

$\left.\begin{array}{c}{g}_u=\frac{\∂g}{\∂u}=g_{i+1,j+1}-g_{i,j}\\ g_v=\frac{\∂g}{\∂v}=g_{i,j+1}-g_{i+1,j}\end{array}\right\}$

(2)计算l×l(如5×5或更大)窗口中灰度的协方差矩阵

$Q=N^{-1}={\left[\begin{array}{cc}\mathrm{\Σ}{g_u}^2& \mathrm{\Σ}{g}_u{g}_v\\ \mathrm{\Σ}{g}_v{g}_u& \mathrm{\Σ}{g_v}^2\end{array}\right]}^{-1}$

其中

$\mathrm{\Σ}{g_u}^2=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\Σ}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\Σ}}}{(g_{i+1,j+1}-g_{i,j})}^2$

$\mathrm{\Σ}{g_v}^2=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\Σ}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\Σ}}}{(g_{i,j+1}-g_{i+1,j})}^2$

$\mathrm{\Σ}{g}_u{g}_v=\mathrm{\Σ}{g}_v{g}_u=\underset{i=c-k}{\overset{c+k-1}{\mathrm{\Σ}}}\underset{j=r-k}{\overset{r+k-1}{\mathrm{\Σ}}}(g_{i+1,j+1}-g_{i,j})(g_{i,j+1}-g_{i+1,j})$

k＝INT(l/2)

(3)计算兴趣值q与w：

$q=\frac{4\mathrm{DetN}}{{\left(\mathrm{trN}\right)}^2}$

$w=\frac{1}{\mathrm{trQ}}=\frac{\mathrm{DetN}}{\mathrm{trN}}$

其中DetN代表矩阵N之行列式；trN代表矩阵N之迹。

(4)确定待选点。如果兴趣点大于给定的阈值，则该像元点为待选点。阈值为经验值，可参考如下：

$\left.\begin{array}{c}{T}_q=0.5-0.75\\ T_w=\begin{array}{}\end{array}\left\{\begin{array}{cc}f\stackrel{\‾}{w}& (f=0.5-1.5)\\ c{w}_c& (c=5)\end{array}\right.\end{array}\right\}$

其中w为权平均值；w_c为权的中值。当q＞T_q同时w＞T_w时，该像元为待选点。

(5)以权值w为依据，选择极值点，即在一定适当窗口中选择w最大的待选点。

图4为本发明的方法中所采用的同名点交叉搜索的程序框图，参照图4。

影像特征是指影像灰度曲面的不连续点。点特征和线特征是影像常见特征信息。点特征主要指明显点，如角点、圆点等。线特征是指影像的“边缘”。所谓特征角点即是影像两“边缘”交叉点；在遥感影像上，特征角点一般对应道路交叉点、建筑角点、河流交叉点等等。

选用经典的Forstner算子对主、辅影像进行特征角点提取，该算子通过计算各像素的Robert’s梯度和以像素(c，r)为中心的一个窗口(如5×5)的灰度协方差矩阵，在影像中寻找具有尽可能小而接近圆的误差椭圆的点作为特征点，详细算法参见相关文献。为了提高效率，只对特征粗匹配所对应的主、辅重叠区域进行特征角点提取；提取过程中，采用了均匀分布策略，并同时给每一特征角点赋予相应的覆盖区域索引号，均匀分布能使特征角点达到全局意义上的均衡化分布。均匀分布策略包括：首先，对主、辅影像重叠区域均匀布置规则格网；其次，使用特征点提取算子，在每个规则子网中提取一个特征角点。覆盖区域索引号在数值上等价于规则子格网编号。

所谓同名点，即是主、辅影像上对应同一地物的特征点；主、辅影像同一地物特征点对可构成同名点对，在基于影像匹配文献中，同名点对有时亦可称为控制点对(Control Points or Tie-Points)。同名点搜索采用区域特征匹配的序贯相似性测度法。由于在之前已进行全局特征粗匹配，同名点搜索范围已大大缩小，进而带来计算复杂度的大幅度下降和匹配正确性的提高。为了丰富同名点信息，采用主、辅影像交叉搜索方法，包括以下步骤：

以主影像特征角点为参考，根据特征粗匹配的变换关系，在辅影像搜索区中寻找同名点；然后把所得同名点对赋予相应的索引号，即同名点图像坐标在主、辅影像中所处的覆盖区域索引号。

以辅影像特征角点为参考，同上一步骤交叉处理。

对所得同名点对合并处理。合并过程中，两组同名点对索引号相同的，则只取其一；反之，两者皆被入选，并纳入最终同名点对列表中。

交叉搜索方法的采用，全面等权地考虑了主、辅影像特征信息，进一步丰富了同名点对，为基于最小二乘的自适应区域分块算法的实施做好了准备。

对于高分辨率SAR数据，阴影是影像的一大重要特征。经典特征角点提取算法不能很好的去除阴影特征点(阴影特征点非硬性地物目标，不满足作为同名点的要求)，这显然会在一定程度上影响总体同名点对的精度性。因此，去除高分辨率SAR影像同名点对中阴影效应将显得格外重要。此外，该方法还采用概率统计的思想，对同名点对进行了粗差剔除处理，从而进一步保证了最终获取同名点对的准确性。

图5为本发明的方法中所采用的LS-MMSE算法流程的程序框图。

基于最小二乘区域自适应分块算法是确保子像素级精度指标的关键技术。高分辨率SAR影像图幅覆盖范围大，局部不同入射角引入的SAR影像近距收缩和传感器系统误差等都是引起影像发生局部变形的直接误差源。因此，对辅影像整体进行单一的几何变换将变得不现实。

辅影像几何变形纠正方法采用仿射变换法，见式(III)：

x₂＝a₀+a₁x+a₂y    (III)

y₂＝b₀+b₁x+b₂y

之所以选择仿射变换，这是有其特性决定的。首先，该变换简单直观；其次，仿射变换不仅能用来处理机载影像，而且也适用于星载影像(高阶的多项式纠正一般只适用于星载数据)；最后，分块子图像局部仿射变换的综合叠加，在影像全局尺度上，等价于拟合一个高阶几何变换函数式。

基于最小二乘区域自适应分块算法的具体步骤如下：

1)初始区域分块，使用公式(III)，获得初始分块子图像同名点对误差方程式(IV)：

$V_x^i=\frac{\∂x}{\∂a_0^i}d{a}_0^i+\frac{\∂x}{\∂a_1^i}d{a}_1^i+\frac{\∂x}{\∂a_2^i}d{a}_2^i-l_x^i$

$V_y^i=\frac{\∂y}{\∂b_0^i}d{b}_0^i+\frac{\∂y}{\∂b_1^i}d{b}_1^i+\frac{\∂y}{\∂b_2^i}d{b}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

其中l_x ⁱ表示观测值x₂ ⁱ同a₀ ⁱ+a₁ ⁱx+a₂ ⁱy计算值的拟合误差，l_y ⁱ表示观测值y₂ ⁱ同b₀ ⁱ+b₁ ⁱx+b₂ ⁱy计算值的拟合误差。

式(IV)表示成矩阵形式，得到式(V)：

V＝AX-L    (V)

式中 $V=\left[\begin{array}{c}{V}_x^i\\ V_y^i\end{array}\right],$ $A=\left[\begin{array}{cccccc}1& x& y& 0& 0& 0\\ 0& 0& 0& 1& x& y\end{array}\right],$ $X=\left[\begin{array}{c}d{a}_0^i\\ d{a}_1^i\\ {\mathrm{da}}_2^i\\ d{b}_0^i\\ d{b}_1^i\\ d{b}_2^i\end{array}\right],$ $L=\left[\begin{array}{c}{l}_x^i\\ l_y^i\end{array}\right].$

2)法化误差方程式(V)：

(A^TA)X＝A^TL

X＝(A^TA)^-1A^TL

迭代计算变形参数改正值da₀ ⁱ，da₁ ⁱ，da₂ ⁱ，db₀ ⁱ，db₁ ⁱ，db₂ ⁱ及对应变形参数a₀ ⁱ，a₁ ⁱ，a₂ ⁱ，b₀ ⁱ，b₁ ⁱ，b₂ ⁱ。经推导，当前变形参数a₀ ⁱ，a₁ ⁱ，a₂ ⁱ，b₀ ⁱ，b₁ ⁱ，b₂ ⁱ可由前一次变形参数a₀ ^i-1，a₁ ^i-1，a₂ ^i-1，b₀ ^i-1，b₁ ^i-1，b₂ ^i-1经由关系式(VI)计算获得：

$a_0^i=a_0^{i-1}+{\mathrm{da}}_0^i+a_0^{i-1}d{a}_1^i{+b}_0^{i-1}{\mathrm{da}}_2^i$

$a_1^i=a_1^{i-1}+a_1^{i-1}{\mathrm{da}}_1^i+b_1^{i-1}d{a}_2^i$

$a_2^i=a_2^{i-1}+a_2^{i-1}{\mathrm{da}}_1^i+b_2^{i-1}{\mathrm{da}}_2^i$

$b_0^i=b_0^{i-1}+d{b}_0^{i-1}+a_0^{i-1}{\mathrm{db}}_1^i+b_0^{i-1}{\mathrm{db}}_2^i---\left(\mathrm{VI}\right)$

$b_1^i=b_1^{i-1}+a_1^{i-1}{\mathrm{db}}_1^i+b_1^{i-1}{\mathrm{db}}_2^i$

$b_2^i=b_2^{i-1}+a_2^{i-1}{\mathrm{db}}_1^i+b_2^{i-1}{\mathrm{db}}_2^i$

3)子图像块分裂，回退，增加控制点处理。

分裂：当步骤2)中变形参数解趋于收敛时，计算该子图像块对应的l_x ⁱ，l_y ⁱ，若计算所得数值大于指定阈值时(子像素级精度，此阈值为1)，则以该图像块为父图像块，一分为四，分裂生成覆盖区域更小的四个子图像块。

回退：当父图像块分裂后，新得到子图像块包含的同名点对(控制点)个数小于必要观测数时，则直接给该新子图像块赋予父图像块对应的变形参数值。为了保证重采样影像子图像块边缘过渡平滑，算法误差精度达到MMSE(Minimum Mean Square Error)，往往需要一定数量的多余观测值。“回退”过程对最终匹配精度是有害的。其实质暗示着，预定匹配精度越高，则需求的控制点数量也越多。因此，充足的控制点数是避免“回退”效应的最直接、有效手段。

增加控制点：为了尽可能避免回退效应，当子图像块内同名点对(控制点)个数小于必要观测数时，重新提取特征角点来加密同名点对(控制点)，进而可继续子图像块分裂、回退和增加控制点整体迭代处理。值得注意的是，鉴于子图像块包含信息量有限，增加控制点处理往往只需迭代1～2次；这也意味着，对于影像信息相对贫乏区域，回退效应不可避免。

4)修正同名点对的图像坐标值。当辅影像对应子图像块分裂、回退和增加控制点过程中止时，重新迭代运算步骤1)，2)；不过，此时需要对x₂ ⁱ，y₂ ⁱ做如下修正处理，直至所有同名点对处理完毕：

$x_2^i=x_2^{i-1}+V_x^{i-1}$

$y_2^i=y_2^{i-1}+V_y^{i-1}$

5)采用间接法，依次对辅影像中包含的子图像块使用对应变形参数进行几何变换，最终得到配准后的影像(重采样选择双线性插值)。

本发明的方法综合考虑了特征匹配的鲁棒性和区域匹配的精确性，总体采用由粗到精的匹配策略实施影像间自动配准：基于二叉树搜索的Hausdorff距离测度全局特征粗匹配算法，很好地解决了Hausdorff距离测度的计算效率瓶颈问题，算法稳健，达到了子像素级的匹配精度，从而缩小了同名点搜索的不定区域，提高了处理效率。同名点匹配采用了序贯相似性检测；交叉搜索方法的引入，等权考虑了主、辅影像的特征信息，丰富了同名点对；同名点匹配获得了影像间像素级的匹配精度。在获取像素级精度的同名点对前提下，基于最小二乘区域自适应分块算法的提出，是保证达到影像全局意义上子像素级(采用2003年7月江苏兴化地区9个航带3米机载SAR影像作为实验数据，配准精度为(0.810076，0.833336)像素)匹配精度的关键环节；自适应区域分块策略的实施，进一步增加了算法的灵活性和鲁棒性。

本发明的高分辨率SAR影像配准处理方法具体应用

应用实例一：

数据构成：2003年7月江苏兴化地区9个航带机载SAR影像(中国科学院电子学研究所机载SAR数据，分辨率3m)。

处理平台：中国科学院中国遥感卫星地面站微波组自行研CAESAR-RegSAR系统软件。

作用：实现了江苏兴化地区9个航带机载高分辨率SAR影像自动/半自动配准，生成大图幅镶嵌影像。

效果：达到了全局子像素级配准精度，配准精度为(0.810076，0.833336)像素。

应用实例二：

数据构成：2005年1月海南三亚地区Ku、X、L、P四波段5个条带机载SAR影像(中国电子科技集团公司第38研究所机载SAR数据，分辨率1m)。

处理平台：中国科学院中国遥感卫星地面站微波组自行研发的CAESAR-RegSAR系统软件。

作用：分别实现了海南三亚地区Ku、X、L、P四波段5个条带机载SAR影像自动/半自动配准，生成大图幅镶嵌影像。

效果：达到了全局子像素级配准精度，配准精度为(0.730816，0.582958)像素。

Claims (6)

1、一种高分辨率SAR影像配准处理方法，其特征是，包括步骤：

(a)创建金字塔影像，提取影像纹理特征，得到边缘信息作为匹配输入数据；

(b)基于所提取的影像纹理特征进行影像全局粗匹配，得到粗配准后的主、辅影像；

(c)提取粗配准后主、辅影像的特征角点，并进行同名点交叉搜索，获得同名点对；

(d)根据同名点对进行自适应区域分块处理；

其中步骤(c)中采用主、辅影像交叉搜索方法，分别以主、辅影像特征角点为参考，根据特征粗匹配的变换关系，在影像搜索区中寻找同名点；然后把所得同名点对赋予相应的索引号；将所得同名点对进行合并处理。

2、如权利要求1所述的高分辨率SAR影像配准处理方法，其特征在于，其中步骤(a)包括：创建金字塔影像；采用Canny算法，对主、辅影像进行边缘提取；对边缘图像进行去除噪声处理；根据边缘长度对所述边缘图像进一步处理。

3、如权利要求1所述的高分辨率SAR影像配准处理方法，其特征在于，其中步骤(b)中采用快速Hausdorff距离算法进行粗匹配。

4、如权利要求1所述的高分辨率SAR影像配准处理方法，其特征在于，其中步骤(c)中采用均匀分布策略提取特征角点，包括：对主、辅影像重叠区域均匀布置规则格网；使用特征点提取算子，在每个规则子网中提取一个特征角点。

5、如权利要求1所述的高分辨率SAR影像配准处理方法，其特征在于，其中在步骤(d)中包括步骤：

(d1)初始区域分块，根据公式(III)获得初始分块子图像同名点对误差方程式(IV)，

x₂＝a₀+a₁x+a₂y                            (III)

y₂＝b₀+b₁x+b₂y

$V_x^i=\frac{\∂x}{{\∂a}_0^i}{\mathrm{da}}_0^i+\frac{\∂x}{{\∂a}_1^i}{\mathrm{da}}_1^i+\frac{\∂x}{{\∂a}_2^i}{\mathrm{da}}_2^i-l_x^i$

$V_y^i=\frac{\∂y}{{\∂b}_0^i}{\mathrm{db}}_0^i+\frac{\∂y}{{\∂b}_1^i}{\mathrm{db}}_1^i+\frac{\∂y}{{\∂b}_2^i}{\mathrm{db}}_2^i-l_y^i---\left(\mathrm{IV}\right)$

$l_x^i=x_2^i-(a_0^i+a_1^{i}x+a_2^{i}y)$

$l_y^i=y_2^i-(b_0^i+b_1^{i}x+b_2^{i}y)$

其中l_x ⁱ表示观测值x₂ ⁱ同a₀ ⁱ+a₁ ⁱx+a₂ ⁱy计算值的拟合误差，l_y ⁱ表示观测值y₂ ⁱ同b₀ ⁱ+b₁ ⁱx+b₂ ⁱy计算值的拟合误差；

(d2)迭代计算变形参数改正值da₀ ⁱ，da₁ ⁱ，da₂ ⁱ，db₀ ⁱ，dv₁ ⁱ，db₂ ⁱ及对应变形参数a₀ ⁱ，a₁ ⁱ，a₂ ⁱ，b₀ ⁱ，b₁ ⁱ，b₂ ⁱ；

(d3)子图像块分裂，回退，增加控制点处理；

(d4)修正同名点对的图像坐标值；

(d5)依次对辅影像中包含的子图像块使用对应变形参数进行几何变换，得到配准后的影像。

6、一种高分辨率SAR影像配准处理系统，其特征在于，包括：

影像数据提取单元，用于提取影像纹理特征数据，得到边缘信息并将该信息送入到影像预处理单元作为粗匹配输入数据；

影像预处理单元，对已提取影像纹理特征的影像进行全局粗匹配，得到配准后的主、辅影像信息，并将相关信息提供给影像配准处理单元处理；

影像配准处理单元，用于提取粗匹配后得到的主、辅影像的特征角点，并进行同名点交叉搜索，获得同名点对；根据同名点对进行自适应区域分块处理；

其中影像配准处理单元具体用于采用主、辅影像交叉搜索方法，分别以主、辅影像特征角点为参考，根据特征粗匹配的变换关系，在影像搜索区中寻找同名点；然后把所得同名点对赋予相应的索引号；将所得同名点对进行合并处理。

Images