一种基于脉冲复用的子孔径SAR图像配准方法

一种基于脉冲复用的子孔径sar图像配准方法
技术领域
1.本发明属于图像配准技术领域，特别是涉及一种基于脉冲复用的子孔径sar图像配准方法。


背景技术：

2.合成孔径雷达(synthetic aperture radar,sar)是一种二维高分辨率成像雷达，能够全天时、全天候工作，具有较高军事和民用价值。sar图像的应用领域广泛，例如目标检测识别、变化检测和数据融合等，作为前期准备技术的sar图像配准是很大程度上影响了后续应用的精度和效果。sar图像配准是指通过寻找参考图像和待配准图像之间的几何变换关系，将它们校正到统一空间坐标系的过程。
3.目前，按照选用的配准因素，图像配准方法主要分为两类：基于灰度的配准方法和基于特征的配准方法。基于灰度的配准方法通过图像像素点的灰度信息来度量参考图像与待配准图像间的相似性，多次迭代找出相似性最高时对应的几何变换参数，较易实现，但该方法计算量大、易受强噪声干扰。基于特征的配准方法构建特征描述子并分别从参考图像和待配准图像中提取出特征点，通过特征点对匹配数量从而估计出最优的几何变换参数。该方法计算量小、鲁棒性强，但sar图像中的相干斑噪声会给配准效果带来较大影响。传统方法sar图像配准方法都是将光学图像配准的方法直接应用于sar图像上，sar图像与光学图像有很大不同，在飞机飞行过程中，雷达以不同的视角对目标区域进行连续成像，同一目标在不同成像结果中的灰度值可能会有很大差异，同时sar图像中随机分布的相干斑噪声导致图像间的相关性很差，光学图像的配准方法因此会受到影响。如何能够削弱sar图像中相干斑噪声带来的弊端，并且与光学图像配准方法结合起来，是一个值得研究的问题。


技术实现要素：

4.本发明提供了一种基于脉冲复用的子孔径sar图像配准方法，从sar原始回波信号出发，区别于纯粹利用光学图像配准技术，合理利用了自聚焦技术和脉冲复用的成像方式，相同脉冲成像所得的子孔径图像之间具有强相关性，利用它们进行配准获取的仿射变换关系更为准确，配准效果得到了提高。
5.为实现上述目的，本发明的具体步骤如下：一种基于脉冲复用的子孔径sar图像配准方法，包括以下步骤：
6.1)获取雷达原始回波数据；
7.2)对步骤1)中获取的原始回波数据进行处理得到子孔径回波信号，通过pfa成像算法对子孔径回波信号进行成像处理得到脉冲复用图像序列并保留距离向加窗、方位向不加窗的中间数据；
8.3)通过相位梯度自聚焦算法对步骤2)获取的图像序列进行处理自动估计并补偿残余的相位误差，得到聚焦效果好的图像并保留迭代累积的总相位误差φ
ε
；
9.4)将步骤3)保留的总相位误差补偿至步骤2)中距离向加窗、方位向不加窗的中间
数据；
10.5)将经过总相位误差补偿后的中间数据在方位向上等分为两个子孔径数据，并各自加汉明窗；
11.6)对距离向和方位向各自加窗后的子孔径数据进行二维傅里叶变换，得到两帧子孔径图像，即实际的参考图像和待配准图像；
12.7)利用sift尺度不变特征变换算法对步骤6)获取的参考图像和待配准图像分别提取特征点；
13.8)根据最近邻距离比原则对参考图像和待配准图像中提取出的特征点进行匹配得到初始匹配点对集合；
14.9)利用ransac随机抽样一致算法剔除初始匹配点对集合中的误匹配点对得到最优匹配点对集合；
15.10)根据最优匹配点对集合得到仿射变换关系；
16.11)将得到的仿射变换关系传递给实际的待配准图像得到配准结果，并将配准结果与参考图像进行非相干叠加得到融合图像。
17.优选的，步骤2)中第k个子孔径回波信号表示为：
18.其中，c为一个复常数，τ和t分别为距离向时间和方位向时间，r(t)为天线相位中心到场景中心点的瞬时斜距，c为电磁波传播速度，t
sub
为子孔径脉冲宽度，k为线性调频斜率，fc为发射信号载频，rect()是矩形窗函数，exp()是指数函数，j是虚数单位，(xm,ym)为场景中的某一点目标，φ和θ分别为天线相位中心的瞬时俯仰角和方位角；利用pfa成像算法对子孔径回波信号进行成像处理，遍历k(k＝1,2,
……
，k)个子孔径回波数据，得到一组连续的脉冲复用图像序列，相邻孔径之间存在重叠的相位历程，孔径重叠率为1/2；在pfa成像算法的二维傅里叶变换步骤前保留距离向加窗、方位向不加窗的中间数据s
wr
(t,τ)，表示为：
19.s
wr
(t,τ)＝wr(τ)
·
s(t,τ)
[0020][0021]
其中，wr(τ)为距离向所加的汉明窗，τ为距离向时间，tr为距离向汉明窗的长度。
[0022]
优选的，步骤3)中总相位误差φ
ε
表示为：其中gn(u)和gn(u-1)分别为相邻两个方位时刻的数据，∫为积分求和运算，∠为内积运算，*为共轭运算，步骤4)中将总相位误差φ
ε
相位补偿至步骤2)中的中间数据，补偿后的数据处于二维时域：
[0023]srpc
(t,τ)＝s
wr
(t,τ)
·
exp(-jφ
ε
)。
[0024]
优选的，步骤5)中两个子孔径数据均为：
[0025]ssub
(t,τ)＝wa(t)
·srpc
(t,τ)，
[0026][0027]
其中，wa(t)为方位向所加的汉明窗，t为方位向时间，ta为方位向汉明窗的长度。
[0028]
优选的，步骤10)中仿射变换关系的模型为：
[0029][0030]
其中(x1,y1)是参考图像上的点，(x2,y2)是对应的待配准图像上的点，是旋转及尺度变换因子，是移位因子。
[0031]
本发明与现有技术相比，具有以下优点：
[0032]
1.本发明考虑了子孔径图像散焦的情况，对全孔径回波数据进行相位梯度自聚焦处理估计出相位误差并补偿至原始数据，避免了加权窗对子孔径划分造成的影响。
[0033]
2.本发明从sar回波信号出发，在图像配准前利用脉冲复用的成像方式，削弱了随机分布的相干斑噪声给sar图像配准带来的弊端，增强了相邻图像之间的相关性，能够提取到更多的匹配点对。
[0034]
3.本发明结合sar原始回波信号和自聚焦处理得到聚焦良好的子孔径图像的过程中充分利用了sar复图像的相位信息，弥补了传统方法忽略相位信息的缺陷。
附图说明
[0035]
图1是本发明按照脉冲复用的读取方式进行子孔径sar图像配准的流程；
[0036]
图2(a)是本发明按照全孔径脉冲复用的方式划分孔径的示意图；
[0037]
图2(b)是本发明按照不重叠子孔径脉冲复用的方式划分孔径的示意图；
[0038]
图2(c)是本发明按照重叠子孔径脉冲复用的方式划分孔径的示意图；
[0039]
图3(a)是根据本发明sar图像配准的步骤得到的第一帧原始成像结果；
[0040]
图3(b)是根据本发明sar图像配准的步骤得到的第二帧原始成像结果；
[0041]
图4(a)是第一帧原始成像结果经过pga处理后的效果图；
[0042]
图4(b)是第二帧原始成像结果经过pga处理后的效果图；
[0043]
图5(a)是由第一帧pga处理结果经子孔径划分得到的子孔径图像1；
[0044]
图5(b)是由第一帧pga处理结果经子孔径划分得到的子孔径图像2；
[0045]
图5(c)是由第二帧pga处理结果经子孔径划分得到的子孔径图像3；
[0046]
图5(d)是由第二帧pga处理结果经子孔径划分得到的子孔径图像4；
[0047]
图6是本发明sar图像参考图像(子孔径图像2)和待配准图像(子孔径图像3)的匹配点对连线图；
[0048]
图7(a)是本发明得到的子孔径图像4的配准结果；
[0049]
图7(b)是本发明将子孔径图像1、子孔径图像2和子孔径图像4的配准结果进行非相干叠加得到的融合图像。
具体实施方式
[0050]
下面将结合附图对本发明的具体步骤进行详细、完整地描述，附图中的实施例仅是本发明在一个场景的应用，并不局限于单一场景，仍可应用于其它场景。
[0051]
结合图1，本发明阐述了一种基于脉冲复用的子孔径sar图像配准方法,实施步骤如下：
[0052]
步骤1、获取原始回波数据；
[0053]
步骤2、按一半脉冲数复用的方式进行连续成像，并保留二维傅里叶变换前方位向不加窗数据；
[0054]
步骤3、相位梯度自聚焦处理估计残余相位误差；
[0055]
步骤4、将相位误差补偿至方位向不加窗数据；
[0056]
步骤5、将数据等分为前后两个子孔径数据并各自在方位向加窗；
[0057]
步骤6、二维傅里叶变换得到两帧子孔径图像；
[0058]
步骤7、利用sift算法提取特征点；
[0059]
步骤8、利用最近邻距离比法进行特征点匹配；
[0060]
步骤9、利用ransac算法剔除误匹配点对；
[0061]
步骤10、获取仿射变换关系；
[0062]
步骤11、得到配准结果并进行非相干叠加得到融合图像。
[0063]
步骤2中按复用一半脉冲数的脉冲进行连续成像，并保留二维傅里叶变换前方位向不加窗数据，实现步骤为：
[0064]
(2a)以图像脉冲数m为规定的脉冲长度，从步骤1获取的雷达原始回波数据的第一个脉冲开始读取m个脉冲得到第一个子孔径回波数据，以m/2个脉冲数为步进读取m个脉冲得到第二个子孔径回波数据，以此类推直至第k个子孔径回波数据，第k个子孔径回波信号可表示为
[0065][0066]
式中，c为一个复常数，τ和t分别为距离向时间和方位向时间，r(t)为天线相位中心(apc)到场景中心点的瞬时斜距，c为电磁波传播速度，t
sub
为子孔径脉冲宽度，k为线性调频斜率，fc为发射信号载频，r
ref
为天线波束中心穿越过场景中心时天线相位中心到场景中心点的斜距，即参考距离。上式的最后一项为残留视频相位项(rvp)，各个脉冲对应的rvp项不同，若不去除会影响场景有效聚焦范围，并会引入额外的几何失真效应。因此需对该项进行补偿，实现如下式,
[0067][0068]
式中，fr和分别为距离向傅里叶变换和傅里叶逆变换，fr为距离向频率。该信号为pfa处理前的输入信号。
[0069]
在平面波前假设条件下，差分距离r(t)-r
ref
表示为
[0070]
r(t)-r
ref
≈-x
m cosφcosθ-y
m sinφcosθ
[0071]
将其带入上式，得到第k个子孔径回波信号的近似表示为
[0072][0073]
利用pfa成像算法对子孔径回波信号进行成像处理，遍历k个子孔径回波数据，最终得到一组连续的脉冲复用图像序列，相邻孔径之间存在重叠的相位历程，孔径重叠率为1/2；
[0074]
(2b)pfa成像算法的最后一步二维傅里叶变换前会在处理数据的方位向和距离向各自加汉明窗，此时保留距离向加窗、方位向不加窗的中间数据供步骤4使用；
[0075]
虽然pfa成像算法包含有运动补偿，但是成像算法中涉及到的近似部分等原因仍会产生相位误差，图像从而散焦。为了弥补此不足，步骤3中利用相位梯度自聚焦算法对(2a)的前两帧图像进行处理自动估计并补偿残余的相位误差，得到聚焦效果好的图像，相位梯度自聚焦处理的实现步骤为：
[0076]
(3a)循环移位：基于相位误差的冗余性的前提，即同一方位向上每个距离单元的相位误差近似相同，相位误差只与方位向有关，将每个距离单元上的最强散射点循环移位到图像中心，以去除散射点的多普勒频率偏移；
[0077]
(3b)加窗：对循环移位后所有距离单元中的最强散射点加窗，保留了相位误差造成的模糊区域，去除了对相位误差估计无用的区域，在保证相位误差估计精度的同时，通过尽可能减小窗的宽度来减轻无用区域信息对相位估计的干扰；
[0078]
(3c)对加窗后的点目标g(x)沿方位向进行傅立叶逆变换，此时数据处于距离压缩相位历史域，记作
[0079]gn
(u)＝|gn(u)|
·
exp(j[φ
ε
(u)+θn(u)])
[0080]
其中n为第n个距离单元，u为方位向坐标，φ
ε
为相位误差，θn(u)为点目标真实相位。
[0081]
(3d)利用方位向相邻数据对相位误差梯度进行估计，如下：
[0082][0083]
其中δφ
ε
(u)＝φ
ε
(u)-φ
ε
(u-1)为相位误差梯度，δθn(u)＝θn(u)-θn(u-1)为真
实相位梯度。经过循环移位后的θn(u)为一固定相位，所以δθn(u)＝0，由此得到了相位误差梯度δφ
ε
(u)为：
[0084][0085]
由于同一方位向上每个距离单元的相位误差近似相同，所以取所有距离单元的相位误差梯度δφ
ε
(u)的平均值作为估计结果仍记为δφ
ε
(u)。获取所有方位向相邻数据的相位梯度误差并求和得到相位误差φ
ε
(u)。
[0086]
(3e)判断估计的相位误差是否小于门限值，若小于门限值，则执行(3f)，若大于门限值，则将(3d)的相位误差补偿至(3c)的距离压缩相位历史域数据，并沿方位向进行傅里叶变换，此时数据由距离压缩相位历史域数据变换至复图像域，将该数据作为新的输入，重复步骤(3a)～(3e)，直至估计的相位误差小于门限值；
[0087]
(3f)将通过迭代最终得到的小于门限值的相位误差补偿至(3c)的距离压缩相位历史域数据，并沿方位向进行傅里叶变换，得到自聚焦结果并保留迭代累积的总相位误差；
[0088]
步骤4中将(3f)保留的最终相位误差φ
ε
补偿至(2b)中距离向加窗、方位向不加窗数据，此时数据处于二维时域：
[0089]srpc
(t,τ)＝wr(τ)
·
s(t,τ)
·
exp(-jφ
ε
)
[0090]
其中
[0091][0092]
步骤5中将相位误差补偿后的距离向加窗、方位向不加窗数据在方位向上等分为两个子孔径数据，并各自加汉明窗：
[0093]ssub
(t,τ)＝wa(t)
·srpc
(t,τ)
[0094]
其中
[0095][0096]
步骤6中对距离向和方位向加窗后的子孔径数据进行二维傅里叶变换，得到两帧子孔径图像，即本发明中实际的参考图像和待配准图像；
[0097]
步骤7中利用sift算法提取特征点，实现步骤为：
[0098]
(7a)利用高斯差分算子构建高斯差分金字塔并提取出特征点；
[0099]
(7b)利用梯度统计直方图确定特征点的主方向；
[0100]
(7c)生成具有尺度不变与旋转不变特性的128维特征描述符；
[0101]
步骤8中根据最近邻距离比原则对参考图像和待配准图像中提取出的特征点进行匹配得到初始匹配点对集合，具体实现步骤为：计算参考图像某一特征点与待配准图像中所有特征点的欧式距离，若最邻近距离与次邻近距离之比小于设定的阈值，则该点与待配准图像中的最邻近距离点构成一对匹配点对，遍历参考图像所有特征点，从而得到初始匹配点对集合；
[0102]
步骤9中利用ransac算法剔除初始匹配点对集合中的误匹配点对，实现步骤为：
[0103]
(9a)从初始匹配点对集合中随机选取p(p≥4)个样本数据，根据p个匹配点对计算它们之间的仿射变换关系，得到仿射变换矩阵m；
[0104]
(9b)计算初始匹配点对集合中剩余匹配点对的待配准图像点经过仿射变换矩阵m变换后与参考图像点之间的欧式距离，若距离小于设定的阈值，则保留该匹配点对，否则，视为误匹配点对进行剔除；
[0105]
(9c)遍历所有剩余匹配点对，筛选得到较为准确的匹配点对集合并记录匹配点对对数；
[0106]
(9d)重复步骤(9a)～(9c)若干次，找到匹配点对对数最多的集合作为最优匹配点对集合；
[0107]
步骤10根据最终匹配点对集合得到仿射变换关系，实现步骤为：将最终匹配点对集合代入参考图像与待配准图像的仿射变换模型，如下式所示，采用最小二乘法得到最优仿射变换模型参数；
[0108][0109]
其中(x1,y1)是参考图像上的点,(x2,y2)是对应的待配准图像上的点。
[0110]
步骤11将得到的仿射变换关系传递给实际的待配准图像得到配准结果，并将配准结果与参考图像进行非相干叠加得到融合图像。
[0111]
下面利用实测数据进行实验，来对本发明的效果做更深刻的说明。
[0112]
1、仿真条件：
[0113]
按照1024个脉冲连续读取雷达回波数据生成大小为1024*512的sar图像序列，硬件平台为：intel(r)core(tm)i5-7200u、16gb ram,软件平台：matlab r2013b。
[0114]
2、实验内容及结果分析：
[0115]
首先按照图2所示的孔径划分方式读取16k脉冲数据的第1个至第1024个脉冲作为第一段数据，接着读取第512个至1536个脉冲作为第二段脉冲数据，然后利用pfa成像算法成像得到连续两帧图像，如图3所示，可以看出图像质量较差，存在散焦现象。为了解决散焦问题，必须对图像进行自聚焦处理，两帧图像经过pga处理后的效果如图4所示，可以看出图像质量得到了明显的改善。图4中的两帧图像并非是本发明中的参考图像和待配准图像，第2个子孔径图像与第3个子孔径图像是实际的参考图像和待配准图像，从图中可以看出它们都是由第一段脉冲和第二段脉冲的重叠部分经成像处理得到的结果，如图5所示，由于相同脉冲所成的两张图像之间的相关性较强，有利于后面的特征点提取。特别需要注意的是，直接由图4中的两帧图像沿方位向进行傅里叶逆变换，随后去除方位向所加的汉明窗，最后等分并通过傅里叶变换得到的两个子孔径图像会存在散焦，这是因为经过自聚焦之后方位向所加的窗形状发生了改变，不能直接去除造成了散焦。为了解决这个问题，本发明将pga处理得到的相位误差补偿给pfa成像前的距离向加窗、方位向不加窗数据，随后在方位向上等分成两个子孔径数据并各自加窗，最后通过傅里叶变换得到了两个聚焦良好的子孔径图像。利用sift算法和ransac算法找出最优匹配点对集合，所有匹配点对间的连线如图6所示。图7给出了本发明得到的配准结果，由于第3个子孔径图像和第4个子孔径图像在同一成
像孔径内，选取的参考点一致，因此图像的几何位置信息完全一致，第3个子孔径图像与第2个子孔径图像之间的仿射变换关系可以直接传递给第4个子孔径图像，第4个子孔径图像经过仿射变换后得到最终的配准结果。为了更直接地看出配准效果，对第一个子孔径图像、第二个子孔径图像及第四个子孔径图像的配准结果进行非相干叠加得到了融合图像，图7中的融合图像清晰地表明了配准结果中的点基本都与参考图像中对应的点重合，配准准确性高。
[0116]
以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。