一种利用特征点Voronoi图优选的星载合成孔径干涉雷达影像配准方法

摘要

本发明涉及一种特征点Voronoi图优选的星载合成孔径干涉雷达影像配准方法，其主要特点在于配准参考点的布设策略较为新颖，充分考虑了特征点的统计信息、专题信息、拓扑信息、度量信息等属性信息，包括以下步骤：（1）获取测区的星载InSAR复数影像对；（2）计算影像对中一幅复数影像的能量图并进行滤波处理；（3）在滤波后的能量图上使用Forstner算子提取特征点；（4）对提取的特征点构建Voronoi图，利用Voronoi图优选特征点；（5）使用优选后的特征点，采用最大谱相关的方法找寻其同名点；（6）根据找寻的同名点对，计算两幅影像间的模型参数；（7）利用计算出的模型参数，将两幅影像统一到同一影像坐标系下。所述方法提高了同名点找寻成功率与相干系数值，使得优选后的特征点在分布上更加合理，提高了模型系数的计算精度，进而提高了影像对的配准精度。

说明书

技术领域

本发明属于遥感影像的合成孔径雷达干涉测量领域，特别涉及星载合成孔 径干涉雷达遥感影像的地形变形监测和三维重建方法。

背景技术

星载合成孔径干涉雷达利用单颗雷达卫星对同一区域两次重访，或多颗雷 达卫星的编队飞行，得到某一区域的两幅复图像(包括强度信息和相位信息), 这两幅复图像之间的相位差包含了地面目标的形变信息和三维位置信息。星载 合成孔径干涉雷达影像配准是指通过一定的手段，使得影像对中同一分辨率单 元（即像元）对应地面上同一回波点单元的技术和方法。在实践中，考虑到应 用效率，并非对影像对重叠区域中的每一个像元逐一进行匹配，通常的做法是 首先在干涉影像对上选取一定数量的同名点，然后利用这些同名点使用有理多 项式模型拟合出影像对间对应像元间的RPC参数，最后，使用RPC参数将影像 对中的一幅影像内插重采样到另一幅影像的坐标系下。

一般来说，主流的干涉雷达复数影像配准的评判指标为相关系数、信噪比、 平均扰动函数三种，星载合成孔径干涉雷达影像配准方法主要有(1)在干涉影像 对中的主影像上布设规则网格，然后以每个网格中心点作为特征点，使用前文 提及的任一评判指标在从影像上搜索同名点，该算法的缺陷是算法选取的特征 点不具有代表性，有一部分甚至落在水域或雷达视线盲区中，从而导致后续同 名点无法找寻或找寻出错；(2)在干涉影像对中的主影像上提取特征点，并根据 特征点的权值对特征点进行排序，删除权值较小的特征点，然后使用筛选的特 征点，利用前文提及的任一评判指标在从影像上搜索同名点，该算法的不足是 仅根据特征点的权值筛选特征点，会导致优选后的特征点不够均匀，从而会降 低后续的多项式拟合精度；(3)基于特征点导引的影像金字塔的匹配策略，首先 在干涉影像对中的主影像上提取特征点，然后对主从影像构建一定层数的金字 塔影像，最后将特征点映射到不同层的金字塔影像上，再使用匹配评判指标在 从影像同级金字塔上找寻同名点，这种基于金字塔的多级匹配策略，极大的增 加了算法的复杂度，每个同名点的匹配涉及到大量的重复冗余计算。

Voronoi图是三角网的对偶图，该图是平分点集中任意具有拓扑关系的两个 点之间的区域而生成的一种几何图形，是点集的平面区域剖分的一种方法，利 用Voronoi图的拓扑关系可以追踪到点集里的每一个特征点，拓扑信息和几何 度量信息是该图固有的两大信息，而特征点具有专题信息和统计信息，将特征 点和其Voronoi图的固有信息综合起来，对特征点进行优选，更具合理性与严 密性。

发明内容

本发明的目的在于解决现有合成孔径干涉雷达影像配准方法的不足，提出 一种基于Forstner特征点提取，Voronoi图特征点优选的复数影像配准方法，克 服了以规则网格中心点为特征点进行同名特征点找寻的随机性与盲目性、基于 特征点权值取舍特征点算法造成的同名特征点分布不均匀性，以及基于金字塔 的多级匹配策算法的重复冗余计算的缺点。

本发明的技术方案具体如下面所描述：一种星载合成孔径干涉雷达影像配 准方法，依靠Forstner算子提取特征点，Voronoi图优选特征点，包括：（1）获 取测区的星载合成孔径干涉雷达单视复数影像对；（2）计算影像对中一幅复数 影像的能量图，并对能量图进行滤波处理；（3）在滤波后的能量图上使用 Forstner算子提取特征点；（4）对提取的海量特征点构建Voronoi图，利用 Voronoi图对特征点进行优选；（5）使用优选后的特征点，在频率域中找寻其 同名点；（6）根据找寻的同名点对，计算两幅影像坐标系间的RPC参数；（7） 利用计算出的RPC参数，将一幅影像重采样到另一幅影像的影像坐标系下，以 实现两幅影像的相对配准。

进一步地，优选的方法是，所述步骤（1）中，使用的复数影像对是执行过 概略配准的两幅具有重叠区域的复数影像。

进一步地，优选的方法是，所述步骤（2）中，当选取一幅复数影像为主影 像时，通常另一幅复数影像称为辅影像或从影像，一般计算主影像的能量图， 能量图的计算按表达式其中re_i,j，im_i,j，A_i,j分别表 示像元(i,j)的实部、虚部和强度值，使用高斯核函数执行滤波操作，二维高斯 核函数的表达式为i,j表示二维平面坐标，δ表示高斯核宽 度，A是核函数的系数，e是自然指数的根，π是圆周率，这里取 A＝1,π＝3.1415926535。

进一步地，优选的方法是，所述步骤（3）中，特征点提取的详细操作方法 包括：

（a）将整幅主影像的幅度影像划分为n个大小为9×9的均匀网格，在每个网 格中，计算各像素点的Robert梯度

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

g表示像素点的灰度值，下标表示像素点的位置，g_u,g_v分别为像素(i,j)处在u,v 方向的Robert梯度，u指向右下方，v垂直于u指向左下方；

（b)计算网格内像素灰度的协方差矩阵

$Q=N^{-1}={\left(\begin{array}{cc}\Sigma g_u^2& \Sigma g_u{g}_v\\ \Sigma g_v{g}_u& \Sigma g_v^2\end{array}\right)}^{-1}$

其中

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

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

$\Sigma g_u{g}_v=\underset{j=r-k}{\overset{r+k-1}{\Sigma }}\underset{i=t-k}{\overset{t+k-1}{\Sigma }}(g_{i+1,j+1}-g_{i,j})(g_{i,j+1}-g_{i+1,j})$

k＝INT(m/2)

m表示网格宽度，k是网格大小的一半，(r,t)是网格中心位置，i,j是求和的循 环变量；

(c)计算像素点的兴趣值q与权w

$q=\frac{4\mathrm{Det}\left(N\right)}{{\mathrm{tr}}^2\left(N\right)}$

$w=\frac{1}{\mathrm{tr}\left(Q\right)}=\frac{\mathrm{Det}\left(N\right)}{\mathrm{tr}\left(N\right)}$

其中符号Det(*)表示行列式；tr(*)表示矩阵的迹。

(d)确定备选点,兴趣值阈值T_q和权重阀值T_w均为经验值，T_q的大致范围为 [0.5～0.75]，取T_q＝0.65，T_w采用下面公式计算：其中为窗口内权值 的均值，f的取值范围为[0.5～1.5]。当q＞T_q且w＞T_w时，该像元被列入备选点列 表。

(e)选定极值点,以权值w为依据，选择极值点，即在一个适当窗口中选择权 值w最大的备选点作为特征点。

进一步地，优选的方法是，所述步骤（4）中，若步骤（3）中提取的特征点 集合用点集P(n)表示，则特征点的优选方法包括：

(a)对特征点集P(n)构建不规则三角网（TIN)，根据三角网的拓扑信息追踪包含 所有特征点的边界多边形，边界多边形的顶点就是边界点。

(b)根据边界多边形按如下原则生成虚拟边界：

I虚拟边界上的顶点位于点集外包中心点到边界多边形的相应顶点的射线 上;

II虚拟边界点和初始边界点的距离等于初始边界上边长的平均值。

(c)将虚拟边界点作为伪特征点加入特征点集合中组成新点集，重新构建TIN， 并在此基础上生成Voronoi图。

(d)特征点优选，遍历整个Voronoi图所有网格，选取一个属性为“可删除”且 优选评判指标值最小的网格对应的特征点，删除该特征点，并将与该特征点具 有拓扑邻接关系的点的属性设置为“不可删除”，其中特征点的优选指标由下式 给出：

$P_i=\frac{W_i\times A_i}{\underset{k=1}{\overset{n}{\Sigma }}{W}_k\times A_k}$

式中，P_i，W_iA_i分别表示第i个特征点的优选指标值，权重，和对应的Voronoi 多边形面积，n代表特征点总数。

(e)重复第(4)步操作，直到剩余特征点数n_remain小于等于预期值n_evaluate，这 里n_evaluate的计算公式如下：

其中rows,cols分别表示主影像的行数（方位向采样数）和列数（斜距向采 样数），gridSize表示影像划分网格时的网格尺寸。

若第(e)步执行到所有剩余特征点的属性值都为“不可删除”，且仍然满足 n_remain＞n_evaluate条件，则清除剩余特征点的所有标记，进入下一次迭代计算，对剩 余特征点集合重复(a)～(e)步骤。

进一步地，优选的方法是，所述步骤（5）中，使用的匹配指标为影像谱相 关法，相关函数计算表达式为：

$R(u,v)=\frac{\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}(f_{j+u,k+v}-\bar{f})(g_{j,k}-\bar{g})}{\sqrt{\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{(f_{j+u,k+v}-\bar{f})}^2}\sqrt{\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{(g_{j,k}-\bar{g})}^2}}$

式中

$\bar{f}=\frac{1}{\mathrm{MN}}\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{f}_{j,k}$

$\bar{g}=\frac{1}{\mathrm{MN}}\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{g}_{j,k}$

M,N表示影像块的长宽，f,g分别是主影像和辅影像能量图的灰度值，是窗口内像素灰度均值。

进一步地，优选的方法是，所述步骤（6）中，RPC模型采用二元二次有理 多项式，多项式模型的计算表达式为：

$\left(\begin{array}{c}{x}_1=a_0+a_{1}x+a_{2}y+a_3\mathrm{xy}+a_4{x}^2+a_5{y}^2\\ y_1=b_0+b_{1}x+b_{2}y+b_3\mathrm{xy}+b_4{x}^2+b_5{y}_2\end{array}\right)$

a₀～a₅，b₀～b₅分别表示多项式系数，(x,y)是辅影像坐标，(x1,y1)是主影像 坐标。

进一步地，优选的方法是，所述步骤（7）中，重采样操作使用sinc三次卷 积法作为影像的内插方法，插值核计算表达式为：

$\sin c\left(\omega \right)=\frac{\sin \omega }{\omega }$

ω表示在插值区域内距离插值位置的距离。

采取了本发明所述的星载合成孔径干涉雷达影像配准方法以后，该种方法 通过特征点的Forstner算子提取，Voronoi图优选，能够充分利用测区地物地貌 特征点的显著性，考虑了特征点分布对影像对几何关系拟合的影响，可以获取 到更高的特征点匹配成功率，提高了干涉相干系数值，进而提高了配准的精度， 提供的方案可在地形形变监测和地表信息三维重建中起到重要的作用，具有较 好的技术效果。

附图说明

下面结合附图对本发明进行进一步详细的描述，以使得本发明的上述优点更 加明确。

图1是本发明所述的星载合成孔径干涉雷达影像配准方法的流程图；

图2是本发明所述星载合成孔径干涉雷达特征点优选示意图；

图3是特征点提取与简化的结果；

图4是采用本发明方法处理后形成的干涉图和相关系数图。

具体实施方式

下面结合附图对本发明作进一步详细的描述。

该种星载合成孔径干涉雷达影像配准方法，通过特征点的Forstner算子提 取，Voronoi图优选，能够充分利用测区地物地貌特征点的显著性与高相干性， 考虑了特征点分布对影像对几何关系拟合的影响，可以获取到更高的特征点匹 配成功率，提高了干涉相干系数值，进而提高了影像配准的精度。

首先，从同一星载干涉雷达传感器对同一区域的两次航过获取的一对干涉 影像对，或两颗不同星载干涉雷达传感器对同一地区的单次访问获得的两幅复 数影像中选取一幅影像作为主影像，另一幅作为从影像。

根据表达式计算主影像的强度图，其中re_i,j，im_i,j， A_i,j分别表示像元(i,j)的实部、虚部和强度值，对强度图进行高斯平滑滤波处理， 高斯平滑滤波器的核大小取为5×5，滤波器核表达式为i,j 表示二维平面坐标，δ表示高斯核宽度，A是核函数的系数，e是自然指数的根， π是圆周率，这里取A＝1,π＝3.1415926535。

在高斯平滑滤波处理后的主影像强度图上提取特征角点或圆点，通过以下 步骤可以实现这一操作：

（a）将整幅主影像的幅度影像划分为n个大小为9×9的均匀网格，在每个网 格中，计算各像素点的Robert梯度

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

g表示影像灰度值，g_u,g_v分别为像素(i,j)处在u,v方向的Robert梯度，u指向右 下方，v垂直于u指向左下方；

（b)计算网格内像素灰度的协方差矩阵

$Q=N^{-1}={\left(\begin{array}{cc}\Sigma g_u^2& \Sigma g_u{g}_v\\ \Sigma g_v{g}_u& \Sigma g_v^2\end{array}\right)}^{-1}$

其中

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

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

$\Sigma g_u{g}_v=\underset{j=r-k}{\overset{r+k-1}{\Sigma }}\underset{i=t-k}{\overset{t+k-1}{\Sigma }}(g_{i+1,j+1}-g_{i,j})(g_{i,j+1}-g_{i+1,j})$

k＝INT(m/2)

m表示影像块大小，k是网格大小的一半，(r,t)是网格中心位置，i,j是求和的 循环变量；

(c)计算像素点的兴趣值q与权w

$q=\frac{4\mathrm{Det}\left(N\right)}{{\mathrm{tr}}^2\left(N\right)}$

$w=\frac{1}{\mathrm{tr}\left(Q\right)}=\frac{\mathrm{Det}\left(N\right)}{\mathrm{tr}\left(N\right)}$

其中符号Det(*)表示行列式；tr(*)表示矩阵的迹；

(d)确定备选点,兴趣值阈值T_q和权重阀值T_w均为经验值，T_q的大致范围为 [0.5～0.75]，取T_q＝0.65，T_w采用下面公式计算：其中为窗口内权值 的均值，f的取值范围为[0.5～1.5]。当q＞T_q且w＞T_w时，该像元被列入备选点列 表；

(e)选定极值点,以权值w为依据，选择极值点，即在一个适当窗口中选择 权值w最大的备选点作为特征点。

在提取的大量特征点中，根据特征点的统计信息、专题信息、拓扑信息和 几何信息，优选特征点。实际操作过程中可以按以下步骤进行：

(1)对特征点集P(n)构建不规则三角网（TIN)，根据三角网的拓扑信息追踪包 含所有特征点的边界多边形，边界多边形的顶点就是边界点，见图2（b）。

(2)根据边界多边形按如下原则生成虚拟边界：

a虚拟边界上的顶点位于点集外包中心点到边界多边形的相应顶点的射线 上;

b虚拟边界点和初始边界点的距离等于初始边界上边长的平均值。

虚拟边界点计算如图2（c）所示。

(3)将虚拟边界点作为伪特征点加入特征点集合中组成新点集，重新构建TIN， 并在此基础上生成Voronoi图。

(4)特征点优选，遍历整个Voronoi图所有网格，选取一个属性为“可删除”且 优选评判指标值最小的网格对应的特征点，删除该特征点，并将与该特征点具 有拓扑邻接关系的点的属性设置为“不可删除”，其中特征点的优选指标由下式 给出：

$P_i=\frac{W_i\times A_i}{\underset{k=1}{\overset{n}{\Sigma }}{W}_k\times A_k}$

式中，P_i，W_iA_i分别表示第i个特征点的优选指标值，权重，和对应的Voronoi 多边形面积，n代表特征点总数，示意图如图2（d）所示。

(5)重复第(4)步操作，直到剩余特征点数n_remain小于等于预期值n_evaluate，这里 n_evaluate的计算公式如下：

其中rows,cols分别表示SAR影像块的行数（方位向采样数）和列数（斜距向采样 数），gridSize表示影像划分网格时的网格尺寸。

若第五步执行到所有剩余特征点的属性值都为“不可删除”，且仍然满足 n_remain＞n_evaluate条件，则清除剩余特征点的所有标记，进入下一次迭代计算，对剩 余特征点集合重复(1)～(5)步骤。

使用的匹配指标为影像谱相关法，相关函数计算表达式为：

$R(u,v)=\frac{\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}(f_{j+u,k+v}-\bar{f})(g_{j,k}-\bar{g})}{\sqrt{\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{(f_{j+u,k+v}-\bar{f})}^2}\sqrt{\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{(g_{j,k}-\bar{g})}^2}}$

式中

$\bar{f}=\frac{1}{\mathrm{MN}}\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{f}_{j,k}$

$\bar{g}=\frac{1}{\mathrm{MN}}\underset{j=1}{\overset{M}{\Sigma }}\underset{k=1}{\overset{N}{\Sigma }}{g}_{j,k}$

M,N表示影像块的长宽，f,g分别是主影像和辅影像能量图的灰度值，是 窗口内像素灰度均值。

采用二元二次有理多项式作为RPC模型，多项式模型的计算表达式为：

$\left(\begin{array}{c}{x}_1=a_0+a_{1}x+a_{2}y+a_3\mathrm{xy}+a_4{x}^2+a_5{y}^2\\ y_1=b_0+b_{1}x+b_{2}y+b_3\mathrm{xy}+b_4{x}^2+b_5{y}_2\end{array}\right)$

a₀～a₅，b₀～b₅分别表示多项式系数，(x,y)是辅影像坐标，(x1,y1)是主影像坐 标。

重采样操作使用sinc三次卷积法作为影像的内插方法，插值核计算表达式 为：

$\sin c\left(\omega \right)=\frac{\sin \omega }{\omega }$

ω表示在插值区域中，距离插值位置的距离。

一般来说,利用多项式S(ω)来逼近理论上的最佳插值函数sinc(ω),多项式 S(ω)的表达式为:

$S\left(\omega \right)=\left(\begin{array}{c}1-2{\left|\omega \right|}^2+{\left|\omega \right|}^3,\\ 4-8\left|\omega \right|+5{\left|\omega \right|}^2-{\left|\omega \right|}^3,\\ 0,\end{array}\right)\left(\begin{array}{c}\left|\omega \right|<1\\ \left|\omega \right|\le 2\\ \left|\omega \right|\ge 2\end{array}\right)$

利用插值位置周围16个邻点的值按下式进行内插：

Q(u,v)＝A·B·C

其中Q(u,v)是(u,v)处的插值结果。

$A={\left(\begin{array}{c}S(1+t)\\ S\left(t\right)\\ S(1-t)\\ S(2-t)\end{array}\right)}^T,{\left(\begin{array}{c}S(1+s)\\ S\left(s\right)\\ S(1-s)\\ S(2-s)\end{array}\right)}^T$

$B=\left(\begin{array}{cccc}p(i-1,j-1)& p(i-1,j)& p(i-1,j+1)& p(i-1,j+2)\\ p(i,j-1)& p(i,j)& p(i,j+1)& p(i,j+2)\\ p(i+1,j-1)& p(i+1,j)& p(i+1,j+1)& p(i+1,j+2)\\ p(i+2,j-1)& p(i+2,j)& p(i+2,j+1)& p(i+2,j+2)\end{array}\right)$

s,t分别表示插值位置处在水平和竖直方向不足一个像素的小数部分，i,j是 插值位置处下水平和竖直方向的整像素部分。

本发明所述的利用特征点Voronoi图优选的星载合成孔径干涉雷达配准方 法，通过使用Forstner算子提取特征点，Voronoi图特征点优选的处理策略，能够 充分利用明显地物地貌特征点在雷达影像上的高相干性，极大地提高了同名特 征点的找寻成功率，增大了匹配指标的数值，从而提高现有合成孔径干涉雷达 配准处理的精度，提供的方案可在地形形变监测和地表三维重建中起到重要的 作用，具有较好的技术效果。

图3示出了特征点提取与简化结果，其中（a）是使用Forstner算子提取的 2025个特征点，（b）使用Voronoi图化简后的保留的100个点，从图（b）可以看 出优选后的待配准点都是明显地物地貌特征点，且几何分布上较均匀，有利于 后期联系点的找寻与主辅影像对间几何变换模型的解算。图4示出了本算法的执 行结果，其中(a)是干涉条纹图，(b)是相关系数图，两幅图均未作多视处理。 实验结果显示，本算法获得了条纹较清晰的干涉图和相关系数值较高的相关系 数图。

上述具体实施例仅仅是示例性的，在本发明的上述教导下，本领域技术人 员可以在上述实施例的基础上进行各种改进和变形，而这些改进或者变形落在 本发明的保护范围内。本领域技术人员应该明白，上面的具体描述只是为了解 释本发明的目的，并非用于限制本发明。本发明的保护范围由权利要求及其等 同物限定。