CN106249236A - 一种星载InSAR长短基线图像联合配准方法 - Google Patents

Abstract

本发明公开了一种星载InSAR长短基线图像联合配准方法，步骤一：利用短基线下的两幅SAR图像反演出成像区域地面点的三维坐标信息；步骤二：选取短基线主图像的特征像素点，将每个特征点对应的地面三维坐标和长基线图像的卫星参数代入多普勒方程，寻找满足成像条件的最优解；步骤三：求解每个特征点在长基线SAR图像中的方位向像素位置和距离向像素位置；步骤四：以特征点在短基线主图像中的像素位置为基准，对特征点在长基线图像中的方位向及距离向偏移量进行多项式拟合，从而得到每个像素点的偏移量；步骤五：利用步骤四中得到的偏移量，对长基线SAR图像进行插值和重采样处理，最终生成与短基线主图像图幅一致的配准后图像。

Description

一种星载InSAR长短基线图像联合配准方法

技术领域

本发明涉及星载合成孔径雷达干涉处理领域，具体地说，是指一种针对干涉合成孔径雷达(简称InSAR)的长短基线图像联合配准方法。

背景技术

InSAR技术主要是利用两幅或两幅以上同一测绘区域的SAR图像，经干涉处理后得到不同信号之间的相位差，利用相位差和高程信息之间的定量关系来获取场景地形高程信息。因此，InSAR技术不仅保持了SAR的成像特征，而且根据其相干特性产生了新的特点。为了从干涉复图像中提取出有效干涉相位，进行干涉处理前需对两次获取的SAR图像进行配准。

图像配准的目的是使两幅图像中同一图像位置的像素对应于地面同一分辨单元，从而保证正确获取同一分辨单元的干涉相位。图像配准的关键在于距离向和方位向偏移量的确定。一般地，偏移量可直接利用两幅图像数据的相关性信息进行估计。然而，对于星载InSAR系统在不同基线下获取的SAR图像，由于时间去相干严重、成像角度差异大等因素，两幅SAR图像极易产生图像噪声、旋转、拉伸等现象，难以直接采用相关法进行图像配准。尤其在多基线干涉处理中，两种基线下的主图像一般不存在相干性。此时不妨联合短基线干涉所得的DEM信息，进行长短基线图像的联合配准。

发明内容

本发明的目的是为了实现星载InSAR长基线SAR图像与短基线SAR图像的联合配准。通过分析星载InSAR干涉处理技术的基本原理，利用基线相距较短的两幅SAR图像相关性较强的特点，首先通过星载InSAR系统短基线干涉处理流程反演出成像区域较精准的数字高程模型(简称DEM)。然后利用短基线主图像中各像素点对应的地面三维坐标，结合卫星位置、卫星速度等成像信息，反推出此地面点在长基线图像中的像素位置。最后对长基线图像进行插值和重采样处理，从而实现与短基线主图像的同名点配准。配准后的长基线图像和短基线主图像图幅一致，二者可联合进行干涉处理并产生质量更好的干涉条纹，并且利用长基线误差传递系数小的特点，获得高精度的DEM。一方面，本发明可用于常规多轨法干涉处理中，利用短基线干涉处理所得DEM来提高长基线图像对的相干性；另一方面，本发明也可用于多基线干涉处理中，通过进行无相干性的长短主图像的配准，为后续多基线解缠和高程反演提供先验信息和输入参数。

本发明采用的技术方案为：一种星载InSAR长短基线图像联合配准方法，包括以下几个步骤：

步骤一：利用短基线下的两幅SAR图像反演出成像区域地面点的三维坐标信息；

步骤二：选取短基线主图像的特征像素点，将每个特征点对应的地面三维坐标和长基线图像的卫星参数代入多普勒方程，寻找满足成像条件的最优解；

步骤三：求解每个特征点在长基线SAR图像中的方位向像素位置和距离向像素位置；

步骤四：以特征点在短基线主图像中的像素位置为基准，对特征点在长基线图像中的方位向及距离向偏移量进行多项式拟合，从而得到每个像素点的偏移量；

步骤五：利用步骤四中得到的偏移量，对长基线SAR图像进行插值和重采样处理，最终生成与短基线主图像图幅一致的配准后图像。

其中，步骤二具体步骤如下：

(1)选取短基线主图像的特征像素点，从第一个图像点开始，在方位向和距离向每隔4个像素点选取一个特征点，使得特征点均匀分布在整幅图像；

(2)寻找特征点对应的长基线成像条件的最优解，具体的将长基线图像对应的卫星位置和卫星速度代入由多普勒方程而衍生的函数，并将每个方位向参数所得函数值排序，取最小的那组函数值对应的位置和速度，即为此特征点在长基线成像条件下的最优卫星速度、位置解向量；

(3)求出所有特征点的最优卫星速度、位置解向量。

本发明的优点在于：

(1)本发明提出了一种全新的长短基线图像联合配准方法，该方法不需要先验DEM。在此之前，关于InSAR的图像配准方法只集中在几何配准和像素配准，没有人做过短基线主辅图像干涉所得DEM与长基线主图像联合配准的研究工作。因此，本发明为星载InSAR领域长短基线图像的联合处理提供了一种新的途径。

(2)本发明可以实现星载InSAR长短基线图像的联合配准，可用于常规多轨法干涉处理中，利用短基线干涉处理所得DEM来提高长基线图像对的相干性。

(3)本发明还可用于多基线干涉处理中，将短基线主图像与长基线主图像配准后，所得结果可为后续干涉相位滤波、长短基线相位解缠等研究工作奠定基础，例如，在精度允许范围内，短基线干涉处理所得DEM信息、长短基线联合处理所得干涉相位信息均可作为多基线干涉处理的先验信息和输入参数；

(4)本发明中利用寻找特征点多普勒方程最优解的方式实现对特征点的精确配准，利用多项式拟合偏移量的方式实现对每个像素点的快速配准，最终可得出配准后图像。

附图说明

图1是本发明的方法流程图；

图2是本发明中短基线SAR主辅图像的干涉处理流程；

图3是本发明中采用的星载干涉SAR测高模型；

图4是本发明中实施示例生成的短基线主图像；

图5是本发明中实施示例生成的短基线辅图像；

图6是本发明中实施示例生成的长基线图像；

图7是本发明中实施示例生成的短基线主辅图像干涉相位；

图8是本发明中实施示例生成的短基线滤波后相位；

图9是本发明中实施示例生成的短基线解缠相位；

图10是本发明中短基线干涉处理生成的数字高程信息；

图11是采用本发明方法经长短基线图像联合配准后的长基线图像。

具体实施方式

下面将结合附图和具体实施方式对本发明作进一步的详细说明。

本发明是一种星载InSAR长短基线图像联合配准方法，总流程如图1所示，包括以下几个步骤：

步骤一、利用短基线下的两幅SAR图像反演出成像区域地面点的三维坐标信息。

通过对扫描区域进行SAR成像仿真(或直接获取真实星载SAR图像)，可得星载SAR在两种不同基线下的单视复图像数据。假设基线较短的两幅SAR图像分别记为s₀和s₁，基线较长的SAR图像记为s₂。首先将s₀作为主图像，s₁作为辅图像，进行短基线干涉处理。短基线干涉处理流程如图2所示。

配准时，采用常规的相关函数法对主辅图像进行精配准和粗配准。两幅SAR图像的相关函数分为实相关函数和复相关函数两类。实相关函数在实际数据处理时，尤其在相干性较低的区域，表现出更为稳健的特性。因此，粗配准时采用实相关函数求取偏移量。又因为在散射特性较为一致且相干性高的区域复相关函数的配准精度高于实相关函数。因此，精配准时首先需对SAR图像进行分块处理以满足每一小块的散射特性较为一致，同时采用复相关函数求取主辅图像中对应像素点的偏移量。其中，粗配准采用的实相关函数计算公式为：

${\mathrm{\ρ}}_r=\frac{\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}\left|s_0(m,n)\right|\left|s_1(m+u,n+v)\right|}{\sqrt{\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}\left(\right|s_0\left(m,n\right){|}^2)}\sqrt{\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}\left(\right|s_1\left(m+u,n+v\right){|}^2)}}---\left(1\right)$

其中：ρ_r为实相关函数，s₀和s₁分别表示主图像和辅图像，(m,n)表示像素中心点的坐标，(u,v)表示滑窗中心点的坐标，M、N表示进行相关计算的窗口的长度和宽度；|·|表示取模操作。精配准采用的复相关函数计算公式为：

${\mathrm{\ρ}}_c=\frac{\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}\left|{s_0}^{\′}(m,n){s_1}^{\′}{(m+u,n+v)}^{*}\right|}{\sqrt{\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}\left({s_0}^{\′}\right(m,n\left){s_0}^{\′}{\left(m,n\right)}^{*}\right)}\sqrt{\underset{m=0}{\overset{M-1}{\Σ}}\underset{n=0}{\overset{N-1}{\Σ}}\left({s_1}^{\′}\right(m+u,n+v\left){s_1}^{\′}{\left(m+u,n+v\right)}^{*}\right)}}---\left(2\right)$

其中：ρ_c为复相关函数，s₀′和s₁′分别表示粗配准后的主图像和粗配准后的辅图像，(m,n)表示图像小块的像素中心点坐标，(u,v)表示滑窗中心点坐标，M、N表示进行相关计算的窗口的长度和宽度，|·|表示取模操作，*表示对复数数据进行复共轭操作。

将精配准后的主辅图像进行复共轭相乘，便可得到干涉条纹图。对干涉条纹图采用基于局部频率估计的斜坡自适应算法进行相位滤波处理，并对滤波后干涉相位采用最小费用流法进行解缠绕处理。这样便得出了解缠相位。可利用控制点的先验信息计算出解缠相位需补偿的相位常数。解缠相位补偿该常数后即得出了绝对相位。将短基线干涉处理所得的绝对相位记为φ_abs。高程反演时，采用的星载干涉SAR测高模型的空间几何关系如图3。图中：A₁为主星，A₂为辅星，P为地面目标，B为InSAR基线长度，B_v为垂直基线长度，B_h为水平基线长度，ξ为基线倾角(基线与水平方向夹角)，θ为雷达视角，H为主天线距地面的高度，r为距地面上目标点的斜距，r+Δr为辅天线距地面上同一目标点的斜距，h为地面目标高度，R_e为地球半径。高程反演时，需利用牛顿迭代法求解以下方程组：

$\left\{\begin{array}{c}r=|{\stackrel{\→}{R}}_{s_1}-{\stackrel{\→}{R}}_p|\\ r+\frac{\mathrm{\λ}}{4\mathrm{\π}}{\mathrm{\φ}}_{abs}=|{\stackrel{\→}{R}}_{s_1}+\stackrel{\→}{B}-{\stackrel{\→}{R}}_p|\\ f_{d_1}=\frac{2}{\mathrm{\λ}r}({\stackrel{\→}{V}}_{s_1}-{\stackrel{\→}{V}}_P)\·({\stackrel{\→}{R}}_{s_1}-{\stackrel{\→}{R}}_p)\end{array}\right.---\left(3\right)$

其中：r为A₁距地面上目标点的斜距，为主星A₁的位置矢量，为地面目标的位置矢量，λ为雷达信号波长，φ_abs为通过干涉处理所得的绝对相位，为A₁和A₂之间的干涉基线矢量，为短基线主图像成像时的多普勒中心频率，为主星A₁的位置矢量，为地面目标的位置矢量。

根据以上方程组，即可求得短基线主图像上每个像素点对应的地面点三维坐标

步骤二：选取短基线主图像的特征像素点，将每个特征点对应的三维坐标和长基线图像的卫星参数代入多普勒方程，寻找满足成像条件的最优解。

为减少计算量，加快程序执行效率，先对短基线主图中的特征点进行配准，步骤四中再拟合特征点偏移量得出短基线主图所有像素点在长基线图像中的偏移量。首先，选取短基线主图像的一些像素点作为特征点，特征点的选取不可过于稀疏且需遍布整幅图像。例如，短基线主图中，由第一个像元开始，按方位向和距离向每五个像素点均匀选取一个作为特征点。所有特征像素点对应的地面点三维坐标均可由步骤一的结果得出。以第一个特征点为例，假设其地面位置矢量为将长基线图像s₂所有方位向对应的卫星位置和卫星速度代入以下函数：

$F({\stackrel{\→}{V}}_{s_2},{\stackrel{\→}{R}}_{s_2})=f_{d_2}-\frac{2}{\mathrm{\λ}r}({\stackrel{\→}{V}}_{s_2}-{\stackrel{\→}{V}}_{P1})\·({\stackrel{\→}{R}}_{s_2}-{\stackrel{\→}{R}}_{p1})---\left(4\right)$

其中，代表与和有关的函数，为长基线图像成像时的多普勒中心频率，为长基线图像s₂所有方位向对应的卫星位置，为s₂所有方位向对应的卫星速度，为该特征点的地面速度，为该特征点的位置矢量。

即为长基线图像成像时满足的多普勒方程。所以将长基线图像s₂所有方位向对应的卫星位置和卫星速度代入后，使得函数最接近0的那组即为此特征点在长基线成像条件下的最优卫星速度、位置解向量。同理可得其他特征点在长基线成像条件下的最优卫星速度、位置解向量。

步骤三：求解每个特征点在长基线SAR图像中的方位向像素位置和距离向像素位置。

特征点在长基线成像条件下的最优卫星速度、位置解向量，对应着该点在长基线图像中唯一的方位向像素位置。依然考虑步骤二中地面位置为的特征点，假设其在长基线成像条件下的最优卫星速度、位置为长基线图像s₂对应的成像卫星到该地面点的斜距为：

$R=|{\stackrel{\→}{R}}_{s_{2}1}-{\stackrel{\→}{R}}_{p1}|---\left(5\right)$

找出解向量后，此特征点在长基线图像中的方位向像素位置可由对应关系直接得到，记方位向像素位置为i。此特征点在长基线图像中的距离向像素位置则可由下式求得：

$R=R_0+\frac{c}{2f_s}j---\left(6\right)$

其中：R₀为成像时的最小斜距，c为光速，f_s为信号采样率，j为所求的特征点距离向像素位置。

求得第一个特征点在长基线图像中的位置(i,j)之后，同理可求得其余特征点在长基线图像中的方位向像素位置和距离向像素位置。

步骤四：以特征点在短基线主图像中的像素位置为基准，对特征点在长基线图像中的方位向及距离向偏移量进行多项式拟合，从而得到每个像素点的偏移量。

首先，将所有特征点在长基线主图像中的像素位置减去特征点在短基线主图像中的像素位置，得到特征点的位置偏移量。短基线主图像每个像素点在长基线图像中的偏移量可通过拟合特征点的偏移量与特征点像素位置的关系来求解。拟合参数可由以下二阶多项式确定：

$\left\{\begin{array}{c}u=a_0+a_{1}m+a_{2}n+a_3{m}^2+a_4{n}^2+a_{5}mn\\ v=b_0+b_{1}m+b_{2}n+b_3{m}^2+b_4{n}^2+b_{5}mn\end{array}\right.---\left(7\right)$

其中，a₀～a₅、b₀～b₅为拟合参数，(m,n)为短基线主图像中的像素坐标，(u,v)为对应点在长基线图像中的偏移量。

确定拟合参数后，可将主图像中的所有像素位置依次代入式(7)，则可得到短基线主图像所有像素点在长基线图像中的偏移量。

步骤五：对长基线SAR图像进行插值和重采样处理，最终生成与短基线主图像图幅一致且可进行联合干涉处理的配准后图像。

首先，对长基线SAR图像进行三次样条插值，采用的插值核函数为：

$i\left(x\right)=\left\{\begin{array}{cc}(\mathrm{\&alpha;}+2)|x{|}^3-(\mathrm{\&alpha;}+3\left)\right|x{|}^3+1& 0\&le;\left|x\right|<1\\ \mathrm{\&alpha;}|x{|}^3-5\mathrm{\&alpha;}|x{|}^2+8\mathrm{\&alpha;}\left|x\right|-4\mathrm{\&alpha;}& 1\&le;\left|x\right|<2\\ 0& \left|x\right|\&GreaterEqual;2\end{array}\right.---\left(8\right)$

其中，α＝-1，x为偏移量，i(x)为四点插值时的偏移量权重函数。

然后，根据短基线主图像所有像素点在长基线图像中的偏移量，对插值后的长基线图像进行重采样，即可得到与短基线主图像图幅一致的配准后长基线图像。

实施例：

本发明为一种星载InSAR长短基线图像联合配准方法，具体实施例为：

步骤一：进行长短基线下的SAR图像仿真，并对短基线下的两幅SAR图像进行干涉处理生成较高精度的DEM。该处理过程主要包括场景设置与SAR成像仿真、短基线SAR主辅图像干涉处理。具体为：

1、SAR回波与成像仿真

(1)根据图3的空间几何关系，设置SAR成像的仿真场景。短基线InSAR系统和长基线InSAR系统的仿真及成像参数如表1所示。

表1InSAR系统仿真参数

参数	数值
		长基线长度	1000m
短基线长度	500m
		波长	0.03125m
传感器速度	7604.175655m/s
		带宽	130MHz
采样率	145MHz

(2)对地面山区场景进行回波仿真并利用Chirp Scaling(CS)算法进行成像处理，得到1700×1000的SAR图像。成像后的短基线主图像如图4所示，短基线辅图像如图5所示，长基线SAR图像如图6所示。

2、短基线SAR主辅图像干涉处理

(1)首先对短基线SAR主辅图像进行基于实相关函数的粗配准，然后将粗配准后的辅图像分成128×128的小块，对每一小块进行基于复相关函数的精配准。

(2)对配准后的主辅图像进行复共轭处理，然后提取其相位信息得到产生的干涉相位如图7所示。

(3)对干涉相位图采用基于局部频率的斜坡自适应滤波算法进行相位滤波，滤波后的干涉相位如图8所示。

(4)对滤波后相位采用最小费用流法进行相位解缠，解缠后相位如图9所示。

(5)根据控制点信息得出绝对相位，根据绝对相位和轨道数据进行高程反演，所得精度较高的数字高程信息如图10所示。

步骤二：选取短基线主图像的特征像素点，将每个特征点对应的三维坐标和长基线图像对应的卫星参数代入多普勒方程，寻找满足成像条件的最优解。

1、选取短基线主图像的特征像素点

在短基线主图像中，从第一个图像点开始，在方位向和距离向每隔4个像素点选取一个特征点，使得特征点均匀分布在整幅图像。

2、寻找特征点对应的长基线成像条件的最优解

(1)对于第一个特征点，将长基线图像s₂对应的卫星位置和卫星速度代入由多普勒方程而衍生的函数，并求取

(2)将每个方位向参数所得到的排序。

(3)使得最小的那组即为此特征点在长基线成像条件下的最优卫星速度、位置解向量。

(4)求出所有特征点的最优卫星速度、位置解向量。

步骤三：求解特征点在长基线SAR图像中的方位向像素位置和距离向像素位置。

(1)根据上一步骤所得的最优卫星速度、位置解向量，得出特征点在长基线图像中的方位向像素位置。

(2)根据特征点的地面三维坐标和对应的卫星位置，得出特征点的成像斜距。

(3)根据特征像素点的成像斜距和仿真时的最小斜距得出特征点在长基线图像中的距离向像素位置。

步骤四：利用多项式拟合，得出短基线主图像中每个像素点在长基线图像中的对应位置。

(1)以特征点在短基线主图像中的像素位置为基准，得出特征点在长基线图像中的方位向及距离向偏移量。

(2)对特征点的像素位置和偏移量进行多项式拟合，得出拟合系数。

(3)将主图像中所有像素位置代入拟合多项式，得出短基线主图像中所有像素点在长基线图像中的偏移量。

步骤五：对长基线SAR图像进行四点三次样条插值和重采样处理，得到配准后图像。经验证知，配准后的长基线图像可以与短基线主图像进行干涉，且生成了较好的干涉条纹。

(1)插值和重采样后的长基线图像如图11所示。

(2)采用本发明所述星载InSAR长短基线图像联合配准方法，对配准后的长基线图像与短基线主图像进行评估,可得出长基线图像与短基线主图像的配准精度如表2所示。

表2短基线主图像与长基线主图像的配准精度

指标	数值
		行配准精度(像素个数)	0.3804
列配准精度(像素个数)	0.5403

本发明主要针对星载InSAR长基线SAR图像与短基线SAR图像因相干性低而无法配准的问题。利用星载InSAR基线越短越易进行干涉处理的特点，首先通过InSAR短基线干涉处理流程反演出成像区域精准的地面三维坐标。然后利用短基线主图像中各像素点对应的地面三维坐标，结合卫星位置、卫星速度等成像信息，反推出此地面点在长基线图像中的像素位置。最后对长基线图像进行插值处理并与短基线主图像进行同名点配准。由于匹配特定像素点方位向位置时的运算量较大，为提升配准效率，本发明中先求解一些特征点的偏移量，再根据特征点偏移量拟合出所有像素点的偏移量。通过实例分析，进一步详述了本发明方法的完整实施过程，验证了本发明方法的合理性与正确性。

Claims (2)

1.一种星载InSAR长短基线图像联合配准方法，其特征在于：包括以下几个步骤：

步骤一：利用短基线下的两幅SAR图像反演出成像区域地面点的三维坐标信息；

步骤二：选取短基线主图像的特征像素点，将每个特征点对应的地面三维坐标和长基线图像的卫星参数代入多普勒方程，寻找满足成像条件的最优解；

步骤三：求解每个特征点在长基线SAR图像中的方位向像素位置和距离向像素位置；

步骤四：以特征点在短基线主图像中的像素位置为基准，对特征点在长基线图像中的方位向及距离向偏移量进行多项式拟合，从而得到每个像素点的偏移量；

步骤五：利用步骤四中得到的偏移量，对长基线SAR图像进行插值和重采样处理，最终生成与短基线主图像图幅一致的配准后图像。

2.根据权利要求1所述的一种星载InSAR长短基线图像联合配准方法，其特征在于：步骤二具体步骤如下：

(1)选取短基线主图像的特征像素点，从第一个图像点开始，在方位向和距离向每隔4个像素点选取一个特征点，使得特征点均匀分布在整幅图像；

(2)寻找特征点对应的长基线成像条件的最优解，具体的将长基线图像对应的卫星位置和卫星速度代入由多普勒方程而衍生的函数，并将每个方位向参数所得函数值排序，取最小的那组函数值对应的位置和速度，即为此特征点在长基线成像条件下的最优卫星速度、位置解向量；

(3)求出所有特征点的最优卫星速度、位置解向量。