基于波数域相干因子的高速超分辨率驻点扫描成像方法与流程

本发明涉及合成孔径成像算法，尤其涉及基于波数域相干因子的高速超分辨率驻点扫描成像方法。

背景技术：

合成孔径(syntheticaperture，sa)成像算法是目前应用最广泛的一种雷达成像技术，该技术的特点是以脉冲压缩技术为根本，相对于传统的电磁成像方法，计算量更小，内存需求更低，是实时成像系统的首选方法。sa采用一种born近似的电场信号模型，为保证图像的精度，不考虑任何远场近似，则最经典的两种sa成像方法是反投影成像算法(backprojectionalgorithm,bpa)和距离偏移成像算法(rangemigrationalgorithm，rma)。两种方法各有特点。bpa的重要优势在于对天线的空间采样方式没有要求，允许非均匀采样。这一点对于非均匀排布的多输入多输出(multiinputmultioutput，mimo)阵列成像和一些复杂的双基成像雷达具有重要意义。另一方面，bpa的抗噪性能也更强。它的主要缺点是计算量非常大，难以应用在实时系统中。为克服这一问题，有人提出了快速bpa(fastbpa，fbpa)和快速分解bpa(fastfactorizedbpa，ffbpa)，这一类算法是以牺牲图像精度为代价换取提高速度，因此需要均衡考虑。相对来说，rma得益于快速傅里叶变换(fft)，其速度优势是与生俱来的，虽然精度比bpa稍低，但是整体性能仍然可接受。rma要求天线的空间采样必须是均匀的。

对于确定的系统参数，为了进一步提高rma的性能，提出了一类超分辨率成像算法。这里超分辨率的含义比较广泛，既包括降低主瓣宽度，也包括压低旁瓣和基底噪声电平。这类方法的主要原理是，将rma的最后一步逆快速傅里叶变换(ifft)用功率谱分析方法代替，因为某些功率谱分析方法得到的谱密度函数比信号的频谱具有更高的频谱分辨率，如capon谱分析和子空间类谱分析方法。但是谱分析方法的计算量和内存资源占用都很高，因此这一类方法多数局限于二维rma成像算法中。

随着压缩感知(compressivesensing,cs)技术的快速发展，雷达成像的一个重要的技术应用就是压缩感知成像。它的应用背景可以这样考虑：对于求解逆问题的优化问题，可以首先考虑将问题稀疏化，进一步在目标函数中加入稀疏项，就可以与传统的正则化项，如全变差项一起求解逆问题。从采用的信号模型看，cs成像问题大致分为两类，第一类采用与sa完全一致的模型，稀疏化以后求解；另一类引入广义反射率的概念，将反射率与探测器方向联系在一起，具有更高精度的同时，也增加了计算量。但是cs成像的重要特点就是具有超分辨率的特性。

由以上分析可见，目前的提出的超分辨率成像算法普遍特点是，计算量大，内存占用率高。对于实时成像系统，尚不存在可用的超分辨率成像算法。

技术实现要素：

针对以上所述超分辨成像方法的技术缺陷，本发明提出基于波数域相干因子的高速超分辨率驻点扫描成像方法及系统，该算法将空间域相干因子扩展到波数域，并与rma相结合，具有实时成像能力；另一方面，该方法在压低图像旁瓣和基底噪声的同时，进一步提高了分辨率。

本发明是通过如下技术方案实现的：

基于波数域相干因子的高速超分辨率驻点扫描成像方法，包括如下步骤：

步骤(1)：天线接收并记录回波信号；

步骤(2)：根据距离偏移成像算法rma处理回波信号，求得反射率函数；

步骤(3)：计算波数域相干因子；

步骤(4)：根据步骤(3)计算得到的相干因子，对步骤(2)得到的反射率函数图像进行校正，最后得到校正后的图像。

所述步骤(2)的步骤为：

步骤(2.1)：对回波信号做二维快速傅里叶变换fft，得到波数域回波信号；

步骤(2.2)：对步骤(2.1)得到的波数域回波信号，进行自由空间衰减和信号脉冲波形补偿；

步骤(2.3)：对步骤(2.2)得到的结果，进行场景中心补偿；

步骤(2.4)：对采样点进行stolt插值，对插值后得到的数据采用三维ifft进行处理，得到反射率函数的图像。

所述步骤(3)的步骤为：

步骤(3.1)：对波数域的回波数据进行数据重排；

步骤(3.2)：对重排后的数据在垂直距离维度上做二维fft；

步骤(3.3)：场景中心补偿；

步骤(3.4)：进行stolt插值，对插值后得到的数据采用三维ifft进行处理，得到反射率函数的相干功率项；

步骤(3.5)：通过步骤(3.4)的到的相干功率项，计算相干因子。

驻点扫描成像是一种单站雷达成像方法，尤其适用于近距离目标成像。

所述步骤(1)的步骤为：单天线同时发射和接收电磁信号，在笛卡尔坐标系内，设天线在空间中形成二维扫描平面与x-y平面平行，用a(x′,y′)表示，且天线的扫描路线与x和y轴平行，并保证成像目标在扫描范围内。天线在x和y方向上均匀采样，采样间隔满足奈奎斯特采样定律。天线位于坐标(x′,y′,zl)接收到的回波信号s(x′,y′,zl；k)表示为

其中，d(x,y,z)是目标所在区域，由天线波束宽度和目标体共同决定；r表示位于(x,y,z)的散射点与天线位置(x′,y′,zl)之间的单程距离；p(k)是信号在波数域上的脉冲波形；σ(x,y,z)是待成像的目标反射率函数；k是信号波数，与频率f相对应。

所述步骤(2)利用距离偏移成像算法求取反射率函数。

所述步骤(2.1)的步骤为：

根据驻定相位原理，对回波信号s(x′,y′,zl；k)的x′和y′做fft，并取傅里叶变换对x′→kx′和y′→ky′，其中，→表示傅里叶变换中的变量关系，得到变换结果为

其中，kz是z轴方向上的波数；kx′、ky′是x′和y′对应的垂直距离向波数。

所述步骤(2.2)中，虽然合成孔径以相位聚焦为主，但是自由空间的衰减补偿有利于提高图像精度；信号调制脉冲波形体现了探测信号在各频率分量上的信噪比，大带宽信号有利于深度向分辨率的提高。

所述步骤(2.2)中自由空间衰减和信号脉冲补偿，实际上是完成距离向的脉冲压缩，补偿结果sc(kx′,ky′,k)为

由于扫描平面与坐标系中的x-y平面平行，因此zl是常数，在公式(5)中将s(kx′,ky′,zl；k)简写为s(kx′,ky′,k)。

所述步骤(2.3)由于傅里叶变换的周期性特点，波数kz对应的z方向的成像范围并不一定包含目标的真实距离位置，有可能周期性地反褶进入成像区，从而难以确定真实距离。

所述步骤(2.3)的步骤为：将需要探测的成像区域事先划定，设在z方向上场景中心与天线平面的距离为hc，则场景中心的补偿结果是

scs(kx′,ky′,k)＝sc(kx′,ky′,k)exp(ikzhc)(6)

所述步骤(2.4)的步骤为：由公式(4)发现，波数kz是k，kx′和ky′的非线性函数，其波数域分布是不均匀的，如果要采用fft加速计算，必须先完成采样点的插值，使得scs(kx′,ky′,k)在三维波数域(kx′,ky′,kz)中均匀排布。

所述插值采用一维的插值方法完成，所述一维的插值方法包括：最近点插值、拉格朗日插值或样条插值，根据需求精度自由选择插值方法，所述插值也称为stolt插值。

插值之后的数据直接利用三维ifft处理，求得反射率的图像为

所述步骤(3.1)的步骤为：

首先引入反射率函数的非相干功率表达式

公式(8)重新整理成如下形式

其中，s(k)是自由空间衰减和信号脉冲波形补偿后的回波信号s(x′,y′,zl；k)的简写，*表示复数的共轭操作；波数k1和k2的取值范围与k相同，都是[kmin,kmax]。r是空间位置的函数，与k1和k2无关，因此通过数据重排后，公式(9)是能够通过傅里叶变换快速计算的。

数据重排的方法如下

ks＝k1-k2(10)

其中，ks是重排后的波数。在物理意义上，回波信号s(x′,y′,k)的相位因子具有exp[ikr(x′,y′；x,y)]的形式，其中r(x′,y′；x,y)指目标与发射和接收天线的双程距离。因此s(k1)s^*(k2)具有exp[i(k1-k2)r(x′,y′；x,y)]形式，即是k1-k2的函数。

由此在公式(10)中关于k1和k2的双重积分能够化简为ks的单重积分。

所述化简的步骤：

假设采样的回波数据包含n个波数域等间隔采样数据，则ki(i＝1,2)对应的n个离散的波数为ki0,ki1,…,ki(n-1)。以s(k1)的数据为行，s^*(k2)的数据为列，构成s(k1)s^*(k2)的二维矩阵。矩阵中的每个对角线元素即对应一个重排后的波数(k1-k2)，即ks。将各对角线的数据相加可得到重排后的新数据e(x′,y′,,zl；ks)。

所述步骤(3.2)的步骤为：

经过回波数据重排后，公式(9)能够写成如下形式：

考察公式(11)的形式完全符合公式(1)的反投影成像算法(backprojectionalgorithm,bpa)计算式，且ks也具有均匀分布格式，因此可以采用rma算法对(11)进行快速计算。根据匹配滤波原理，将公式(11)变换成回波形式，并在两端同时做二维的fft得到

e(kx′,ky′,zl；ks)＝∫∫∫d(x,y,z)ip(x,y,z)exp(-ikx′x-iky′y-ikz′|z-zl|)dxdydz(12)

其中，e(kx′,ky′,zl；ks)是重排后数据e(x′,y′,ks)在垂直距离向的傅里叶变换结果。

所述步骤(3.3)过程与步骤(2.3)相似，首先完成场景中心补偿，保证ip(x,y,z)的中心与场景中心对齐，补偿结果ec(kx′,ky′,zl；ks)为

ec(kx′,ky′,zl；ks)＝e(kx′,ky′,zl；ks)exp(ikz′hc)(14)

所述步骤(3.4)的步骤为：

根据公式(13)，在kz′上的数据并不满足均匀分布的采样要求，需要先进行插值，使得ec(kx′,ky′,zl；ks)在三维波数域(kx′,ky′,kz′)中均匀排布。

经过stolt插值后的回波数据，能够直接经过三维的ifft求得反射率的相干功率项ip(x,y,z)；

所述步骤(3.5)的相干因子cf(x,y,z)定义为

公式(16)中的参数α和β用于调节相干因子的性能，根据经验，一般取α∈[0.5,4.5],β∈[1/2,1/500]。不同系统中参数的取值不同，但是一旦确定，相干因子的性能会非常稳定，不需要进一步修改。

所述步骤(4)的步骤为：

由步骤(3)得到的相干因子，计算最终图像为

σm(x,y,z)＝σ(x,y,z)cf(x,y,z)(17)。

与现有技术相比，本发明的有益效果是：

(1)采用相干因子与rma相结合的成像算法，相对于传统的rma算法具有更高的分辨率，更低的最大旁瓣，并有利于压低基底噪声电平；公式(8)的快速计算是本方法的核心。

(2)采用波数域方法计算相干因子，结合rma具有快速成像能力，能够满足实时成像系统的要求；算法结构有利于并行运算，能够进一步提高成像速度；

(3)引入参数α和β用于调节相干因子的性能，能够平衡相干因子对弱散射点的压制作用，且一旦参数确定，算法性能非常稳定；

(4)波数域数据重排，将反射率函数非相干功率项变换成适应于rma处理的固定形式，并且不需要插值操作，有利于提高成像效率；

附图说明

图1是单天线驻点扫描成像原理图；

图2是基于波数域相干因子的高速超分辨率驻点扫描成像方法原理框图；

图3是x轴方向的点扩展函数；

图4是y轴方向的点扩展函数；

图5是用于散射场计算的八叶金属扇形目标的示意图；

图6是传统rma电磁仿真成像结果；

图7是相干因子rma电磁仿真成像结果；

图8是用于成像实验的八叶金属扇形目标的图片；

图9是传统rma的实验成像结果；

图10是相干因子rma的实验成像结果；

图11是传统bpa的实验成像结果；

图12是波数域数据重排的数据阵列。

图中：1、用于发射和接收电磁波的天线，2、目标域。

具体实施方式

下面结合附图与实施例对本发明作进一步说明。

单天线驻点扫描成像系统的工作方式，如图1所示，用于发射和接收电磁波的天线(1)用于同时发射和接收电磁信号，用于发射和接收电磁波的天线(1)在与x-y平行的平面上进行二维驻点扫描，形成的二维孔径表示为a(x′,y′)。网格的节点表示天线停驻扫描的位置，相邻节点应满足nyquist采样定律。天线扫描应覆盖目标域(2)，任何待成像的目标都应位于目标域d(x,y,z)的范围内。

实施例一：点扩展函数

本实施例分别计算了传统rma和相干因子rma两种成像算法的点扩展函数性能。用于计算的主要参数在表1中列出。

表1点扩展函数计算采用的主要参数

点扩展函数即是采用理想的点目标，由sa的信号模型产生回波信号，并利用成像算法得到的点目标的图像。由图2所示，相干因子rma的详细步骤如下：

步骤1：天线接收并记录回波信号；

步骤2：反射率计算。

这一步实际上就是采用rma直接重构目标的三维图像，主要包括如下几步：

步骤2.1：对回波信号做二维fft；

步骤2.2：自由空间衰减和信号脉冲波形补偿；

步骤2.3：场景中心补偿；

步骤2.4：stolt插值和三维ifft。

stolt插值是rma中制约其分辨率的主要因素，同时影响成像精度。根据成像速度和精度要求，一般采用线性插值，如拉格朗日插值或sinc插值等。

步骤3：波数域相干因子计算。主要包括如下几步：

步骤3.1：波数域数据重排；

这一步是快速计算相干因子的关键，方法是首先将回波数据s(k1)s^*(k2)写成如图12所示的矩阵形式。假设采样的回波数据包含n个波数域等间隔采样数据，则ki(i＝1,2)对应的n个离散的波数为ki0,ki1,…,ki(n-1)。以s(k1)的数据为行，s^*(k2)的数据为列，构成s(k1)s^*(k2)的二维矩阵。该矩阵中的每个对角线元素即对应一个重排后的波数(k1-k2)，即ks。将各对角线的数据相加可得到重排后的新数据e(ks)。

步骤3.2：在垂直距离向两个维度上做fft；

由公式(9)可见，反射率的相干功率项ip(x,y,z)是在垂直距离向的空间域上定义的，因此用于重排的数据也是进行自由空间衰减和信号脉冲波形补偿后的空间域数据，需要变换到波数域才能进一步处理。

步骤3.3：场景中心补偿；

步骤3.4：stolt插值和三维ifft；

步骤3.5：计算相干因子；

步骤4：计算相干因子校正后的图像。

分别采用传统的rma和相干因子rma两种方法计算得到的点扩展函数结果如图3和图4所示。其中图3是x方向上的点扩展函数，图4是y方向上的点扩展函数，各图中十字曲线是传统rma方法的处理结果，点状曲线是相干因子rma方法的处理结果。由于图3和图4得到的曲线结果完全相同，只采用图3结果进行进一步分析，结果如表2所示。

表2中的几种参量的计算方法是：

半峰全宽(fwhm)：靠近0db处的两个负峰值被认为是主瓣的两端，主瓣幅值下降到3db处对应的宽度是fwhm；

背景噪声电平：主瓣采样点之外，所有采样点的幅值平均值；

最大旁瓣幅值：主瓣外区域中采样点的最大幅值。

表2点扩展函数的分析结果

对表2的结果进行总结。采用相干因子后，相对于传统的rma：

1)点扩展函数主瓣的fmhw减少了0.48mm，相当于分辨率提高了大约20％；

2)背景噪声的平均电平下降了22.190db。该参量对于稀疏阵列成像格外有意义；

3)最大旁瓣幅值下降了11.300db。

实施例二：电磁仿真成像

该实施例主要用于验证相干因子rma对连续目标体的成像性能。仿真中采用的主要参数如表1所示，唯一的不同之处是α和β参数取值不同，即α＝1.5，β＝0.00375。

仿真包括两步：

步骤一：利用矩量法产生目标的散射场作为回波数据；

成像的目标体如图5所示，是八叶扇形理想电导体的二维模型。电磁仿真采用赫兹电偶极子源，源和点探测器位于同一位置，且只考虑同方向极化电场。

步骤二：分别采用传统rma和相干因子rma对回波数据进行处理，完成图像重构。

成像结果如图6和图7所示。图6是传统rma的成像结果，图7是相干因子rma的成像结果。图像采用灰度图，动态范围是30db。图像的背景噪声在-30db以下，因此在图中被滤去。图6中可以明显看到在扇叶附近存在较多的旁瓣，使得边缘处变得模糊。对比于图7，采用相干因子后，扇叶附近的旁瓣被极大地压制，边缘处非常清晰，另外可以发现图7中的扇叶形状更为尖锐，目标中心处轮廓更清晰，这些都是相干因子rma的分辨率被提高的原因。

实施例三：成像实验结果

本实施例通过单天线驻点扫描成像实验验证相干因子rma的性能。实验平台由两部分构成，分别是扫描平台和矢量网络分析仪。扫描平台用于按照规定方式平动被成像样品，允许在平面的两个维度上定义移动范围和移动步进；矢量网络分析仪(agilent，n5247a)用于产生、发射并接收回波信号，由于仪器结构是相干的，因此不需要进行载波的波形校正，采样数据s11或s22可直接处理。

成像的目标体是八叶金属扇形体，如图8所示。目标体被固定在扫描平台上，实验的参数在表3中列出。

表3实验参数

实验的成像结果如图9,10和11所示。其中图9是传统rma的成像结果，图10是相干因子rma成像结果，图11是传统bpa成像结果。各图的动态范围是20db。在图9中，基底噪声和旁瓣都比较明显，分布于图像各处。相对于仿真结果，图像质量下降的原因主要是仪器噪声因素。相比之下，图11的质量要好一些，这是由于bpa的抗噪声性能更强。图10的质量是最好的，图像的背景非常干净，基底噪声和旁瓣的压制效果明显，而且扇叶形状更为尖锐，证明图像的分辨率有所提高。

比较三种成像算法的效率。实验中采用计算机服务器运行成像程序，计算机配置为：双核intel(r)xeon(r)cpue5-2650，64gbram。程序运行时间如表4所示。相干因子rma的运行时间大约是传统rma运行时间的2倍，可以推断波数域数据重排操作几乎不影响算法效率，主要运算量在于两次rma运算，因此如果采用并行运算能很好优化该程序。传统bpa的运算时间非常长，几乎不可能用于实时成像系统。

表4各种成像方法所用的时间

上述虽然结合附图对本发明的具体实施方式进行了描述，但并非对本发明保护范围的限制，所属领域技术人员应该明白，在本发明的技术方案的基础上，本领域技术人员不需要付出创造性劳动即可做出的各种修改或变形仍在本发明的保护范围以内。