本发明涉及雷达三维成像技术领域，特别是涉及一种窄带无源雷达三维成像方法及系统。
背景技术：
无源雷达成像技术是无源雷达研究领域的一个热点课题。目前无源雷达成像技术主要对二维成像进行了集中研究。早在本世纪初，美国伊利诺伊大学就利用民用外辐射源信号采用傅里叶逆变换技术成功对目标成像，在此基础上，二维窄带无源雷达成像的频域和时域算法也被提出。
目前三维成像是雷达成像领域的前沿和热点课题，现在的雷达三维成像技术主要包括合成孔径雷达、逆合成孔径雷达、太赫兹等技术，由以上技术组成的三维雷达成像系统多为主动雷达系统，为了得到高分辨成像结果，需要采用宽带或超宽带信号，但并没有对窄带进行研究，因此如何设计一种窄带无源雷达三维成像成为本领域亟需解决的技术问题。
技术实现要素：
本发明的目的是提供一种窄带无源雷达三维成像方法及系统，以实现利用窄带进行无源雷达三维成像。
为实现上述目的，本发明提供一种窄带无源雷达三维成像方法，其特征在于，所述方法包括：
根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系；
利用接收机获取成像目标的两种窄带调频回波信号；所述两种窄带调频回波信号分别为绕z轴旋转的窄带调频回波信号以及先绕z轴旋转再绕y轴旋转的窄带调频回波信号；其中，所述三维旋转成像模型直角坐标系中的x轴由接收机、转台转动中心和外辐射源的位置确定的，y轴为与所述x轴垂直的轴线，由接收机、转台转动中心和外辐射源构成的平面为xoy平面，z轴为与xoy平面垂直垂直的轴线；
对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵；
根据所述离散回波信号矩阵确定信号匹配矩阵；
根据所述离散回波信号矩阵和所述信号匹配矩阵得到聚焦成像。
可选的，所述对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵，具体包括：
对两种所述窄带调频回波信号分别进行混频和补偿处理，获得两种补偿后的窄带调频回波信号；
对两种所述补偿后的窄带调频回波信号进行合成处理，获得合成回波信号；
对所述合成回波信号进行离散采样，获得离散回波信号矩阵。
可选的，所述根据所述离散回波信号矩阵确定信号匹配矩阵，具体包括：
对成像目标所在的立体空间进行网格划分；
根据所述离散回波信号矩阵确定网格处散射点的离散回波信号矩阵；
根据所述网格处散射点的离散回波信号矩阵确定信号匹配矩阵。
可选的，利用接收机获取成像目标的两种窄带调频回波信号，具体公式为：
其中，a为窄带调频回波信号的幅度，f为窄带调频回波信号的频率，φ为窄带调频回波信号的初始相位，σ为散射点散射强度，c为光速，λ为波长，rt为外辐射源到坐标原点的距离，rr为接收机到坐标原点的距离，α为绕z轴转动总角度，β为1/2双基角，ωz为绕z轴旋转的角速度，tz为绕z轴旋转的时间，ωy为绕y轴旋转的角速度，ty为绕y轴旋转的时间，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标，ry(ty)为成像目标的散射点绕y轴旋转的斜距历程，rz(tz)为成像目标的散射点绕z轴旋转的斜距历程。
可选的，所述对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵，具体公式为：
其中，σ为散射点散射强度，λ为波长，β为1/2双基角，α为绕z轴转动总角度，ωz为绕z轴旋转的角速度，ωy为绕y轴旋转的角速度，δty、δtz分别为对合成回波信号离散化时ty、tz的离散时间采样步长，l为总采样点数，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标。
可选的，在获取接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系之前，所述方法还包括：
获取目标，并对目标进行运动补偿，获得能够在转台模型下研究的成像目标。
本发明还提供一种窄带无源雷达三维成像系统，所述系统包括：
构建模块，用于根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系；
获取模块，用于利用接收机获取成像目标的两种窄带调频回波信号；所述两种窄带调频回波信号分别为绕z轴旋转的窄带调频回波信号以及先绕z轴旋转再绕y轴旋转的窄带调频回波信号；其中，x轴由接收机、转台转动中心和外辐射源的位置确定的，y轴为与所述x轴垂直的轴线，由接收机、转台转动中心和外辐射源构成的平面为xoy平面，z轴为与xoy平面垂直垂直的轴线；
预处理模块，用于对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵；
匹配矩阵确定模块，用于根据所述离散回波信号矩阵确定信号匹配矩阵；
聚焦成像模块，用于根据所述离散回波信号矩阵和所述信号匹配矩阵确定得到聚焦成像。
可选的，所述预处理模块，具体包括：
混频和补偿处理单元，用于对两种所述窄带调频回波信号分别进行混频和补偿处理，获得两种补偿后的窄带调频回波信号；
合成处理单元，用于对两种所述补偿后的窄带调频回波信号进行合成处理，获得合成回波信号；
离散采样单元，用于对所述合成回波信号进行离散采样，获得离散回波信号矩阵。
可选的，所述匹配矩阵确定模块，具体包括：
网格划分单元，用于对成像目标所在的立体空间进行网格划分；
离散回波信号矩阵确定单元，用于根据所述离散回波信号矩阵确定网格处散射点的离散回波信号矩阵；
匹配矩阵单元，用于根据所述网格处散射点的离散回波信号矩阵确定信号匹配矩阵。
可选的，所述获取模块，用于利用接收机获取成像目标的两种窄带调频回波信号，具体公式为：
其中，a为窄带调频回波信号的幅度，f为窄带调频回波信号的频率，φ为窄带调频回波信号的初始相位，σ为散射点散射强度，c为光速，λ为波长，rt为外辐射源到坐标原点的距离，rr为接收机到坐标原点的距离，α为绕z轴转动总角度，β为1/2双基角，ωz为绕z轴旋转的角速度，tz为绕z轴旋转的时间，ωy为绕y轴旋转的角速度，ty为绕y轴旋转的时间，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标，ry(ty)为成像目标的散射点绕y轴旋转的斜距历程，rz(tz)为成像目标的散射点绕z轴旋转的斜距历程。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明先根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系；利用接收机获取成像目标的两种窄带调频回波信号，所述两种窄带调频回波信号分别为绕z轴旋转的窄带调频回波信号以及先绕z轴旋转再绕y轴旋转的窄带调频回波信号；其次对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵；然后根据所述离散回波信号矩阵确定信号匹配矩阵；最后根据所述离散回波信号矩阵和所述信号匹配矩阵确定得到聚焦成像。本发明不仅能够实现窄带无源雷达三维成像，同时利用在时间维上进行计算避开了现有方法需要建立四维回波信号矩阵的过程，减小了存储负担，降低了算法复杂度。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例窄带无源雷达三维成像方法流程图；
图2为本发明实施例窄带无源雷达三维成像系统结构框图；
图3为本发明实施例窄带无源雷达三维成像系统旋转结构图；
图4为本发明实施例第一聚焦成像仿真图；
图5为本发明实施例第二聚焦成像仿真图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种窄带无源雷达三维成像方法及系统，以实现利用窄带进行无源雷达三维成像。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明研究发现，如果转台绕着固定转动轴转动，系统不具有三维成像能力。既使转台同时绕着x轴、y轴、z轴做匀速三维转动，其效果相当于绕着一固定转动轴转动，所以也无法做三维成像。因此本发明提出了基于转台分时段绕不同转动轴转动然后将回波信号合成进行三维成像的模型。
图1为本发明实施例窄带无源雷达三维成像方法流程图；如图1所示，本发明提供一种窄带无源雷达三维成像方法，所述方法包括：
步骤100：根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系。
步骤200：利用接收机获取成像目标的两种窄带调频回波信号；所述两种窄带调频回波信号分别为绕z轴旋转的窄带调频回波信号以及先绕z轴旋转再绕y轴旋转的窄带调频回波信号；其中，x轴由接收机、转台转动中心和外辐射源的位置确定的，y轴为与所述x轴垂直的轴线，由接收机、转台转动中心和外辐射源构成的平面为xoy平面，z轴为与xoy平面垂直垂直的轴线。
步骤300：对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵；所述离散回波信号矩阵为行矩阵，即一维行向量。
步骤400：根据所述离散回波信号矩阵确定信号匹配矩阵。
步骤500：根据所述离散回波信号矩阵和所述信号匹配矩阵确定得到聚焦成像。
在步骤100之前还包括：获取目标，并对目标进行运动补偿，获得能够在转台模型下研究的成像目标。
下面对各个步骤进行具体详细论述：
步骤100：根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系。
如图3所示，所述三维旋转成像模型直角坐标系是以转台转动中心、外辐射源、接收机所在的平面为xoy平面，并以外辐射源、坐标原点、接收机的位置所构成角的角平分线为x轴，垂直于x轴的方向为y轴，经过o点且垂直于xoy平面的的轴线为z轴。外辐射源和接收机分别位于(rtcosβ,rtsinβ,0)、(rrcosβ,-rrsinβ,0)，其中rt、rr分别为外辐射源和接收机到坐标原点的距离，系统的双基角为2β。假设初始时刻成像目标上一散射点的三维坐标为(x0,y0,z0)，成像目标先绕z轴转动(偏航)、再绕y轴转动(纵摇)，转台绕z轴转动角速度为ωz，经过时间tz散射点坐标为(x0cosωztz-y0sinωztz,x0sinωztz+y0cosωztz,z0)，如果其绕z轴转动总角度为α，则其最终转动至(x0cosα-y0sinα,x0sinα+y0cosα,z0)。在此基础上，转台绕着y轴开始转动，转动角速度为ωy并重新开始计时，经过时间ty散射点转动至(xy(ty),yy(ty),zy(ty))，其具体形式为：
步骤200：利用接收机获取成像目标的两种窄带调频回波信号；具体公式为：
其中，a为窄带调频回波信号的幅度，f为窄带调频回波信号的频率，φ为窄带调频回波信号的初始相位，σ为散射点散射强度，c为光速，λ为波长，rt为外辐射源到坐标原点的距离，rr为接收机到坐标原点的距离，α为绕z轴转动总角度，β为1/2双基角，ωz为绕z轴旋转的角速度，tz为绕z轴旋转的时间，ωy为绕y轴旋转的角速度，ty为绕y轴旋转的时间，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标，ry(ty)为成像目标的散射点绕y轴旋转的斜距历程，rz(tz)为成像目标的散射点绕z轴旋转的斜距历程。
在上面两个转动过程中，由于外辐射源和接收机到目标距离远大于目标的转动尺寸，目标的散射点绕z轴旋转的斜距历程和绕y轴旋转的斜距历程可以近似表示为：
rz(tz)＝rr+rt-2cosβ(x0cosωztz-y0sinωztz)
ry(ty)＝rr+rt-2cosβ(x0cosαcosωyty-y0sinαcosωyty+z0sinωyty)
步骤300：所述对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵，具体包括：
步骤301：对两种所述窄带调频回波信号分别进行混频和补偿无用固定相位和幅度处理，获得两种补偿后的窄带调频回波信号；具体公式为：
其中，σ为散射点散射强度，λ为波长，α为绕z轴转动总角度，β为1/2双基角，ωz为绕z轴旋转的角速度，tz为绕z轴旋转的时间，ωy为绕y轴旋转的角速度，ty为绕y轴旋转的时间，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标。
步骤302：对两种所述补偿后的窄带调频回波信号进行合成处理，获得合成回波信号。
本发明对两种所述补偿后的窄带调频回波信号进行乘积，获得合成回波信号；具体公式为：
步骤303：对所述合成回波信号进行离散采样，获得离散回波信号矩阵；具体公式为：
其中，σ为散射点散射强度，λ为波长，β为1/2双基角，α为绕z轴转动总角度，ωz为绕z轴旋转的角速度，ωy为绕y轴旋转的角速度，δty、δtz分别为对合成回波信号离散化时ty、tz的离散时间采样步长，l为总采样点数，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标。
步骤400：所述根据所述离散回波信号矩阵确定信号匹配矩阵，具体包括：
步骤401：对成像目标所在的立体空间进行网格划分，获得多个网格。
本发明对成像目标所在的立体空间进行网格划分，分别在方位、距离、高度维上将成像目标划分为m列n行k页，其三轴方向的步长分别为δx、δy、δz。
步骤402：根据所述离散回波信号矩阵确定网格处散射点的离散回波信号矩阵；具体公式为：
m＝1,2,…,m；n＝1,2,…,n；k＝1,2,…,k；l＝1,2,…,l
其中，σm,n,k为网格处散射点的散射强度，λ为波长，β为1/2双基角，α为绕z轴转动总角度，δx、δy、δz分别为三轴方向的步长，ωz为绕z轴旋转的角速度，ωy为绕y轴旋转的角速度，δty、δtz分别为对合成回波信号离散化时ty、tz的离散时间采样步长，ωz＝ωzδtz、ωy＝ωyδty，l为总采样点数。
步骤403：根据所述网格处散射点的离散回波信号矩阵确定信号匹配矩阵；具体公式为：
m＝1,2,…,m；n＝1,2,…,n；k＝1,2,…,k；l＝1,2,…,l
其中，λ为波长，β为1/2双基角，α为绕z轴转动总角度，δx、δy、δz分别为三轴方向的步长，ωz为绕z轴旋转的角速度，ωy为绕y轴旋转的角速度，δty、δtz分别为对合成回波信号离散化时ty、tz的离散时间采样步长，ωz＝ωzδtz、ωy＝ωyδty，l为总采样点数。
步骤500：根据所述离散回波信号矩阵和所述信号匹配矩阵确定得到聚焦成像；具体公式为：
m＝1,2,…,m；n＝1,2,…,n；k＝1,2,…,k；l＝1,2,…,l
其中，σ为散射点散射强度，λ为波长，β为1/2双基角，α为绕z轴转动总角度，ωz为绕z轴旋转的角速度，ωy为绕y轴旋转的角速度，δty、δtz分别为对合成回波信号离散化时ty、tz的离散时间采样步长，l为总采样点数，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标，δx、δy、δz分别为三轴方向的步长，ωz＝ωzδtz、ωy＝ωyδty。
当mδx＝x0、nδy＝y0、kδz＝z0时，成像结果出现峰值lσ2，因此能够重建散射点位置。在x方向、y方向、z方向上，其成像结果分别为：
其中，j0(·)为零阶第一类贝塞尔函数。由上式可知，在x方向、y方向上的成像结果受到l、ωz、ωy、α、β等参数的影响。在z方向上为贝塞尔函数形式，其分辨率为0.18λ/cosβrad，峰值旁瓣比为-7.9db。
本发明对不同转动轴的两种窄带调频回波信号进行了合成，克服了现有方法只有z轴转动、未进行回波信号合成导致只能做二维成像的特殊情况。
图2为本发明实施例窄带无源雷达三维成像系统结构框图；如图2所示，本发明还提供一种窄带无源雷达三维成像系统，所述系统包括：
构建模块1，用于根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系。
获取模块2，用于利用接收机获取成像目标的两种窄带调频回波信号；所述两种窄带调频回波信号分别为绕z轴旋转的窄带调频回波信号以及先绕z轴旋转再绕y轴旋转的窄带调频回波信号；其中，x轴由接收机、转台转动中心和外辐射源的位置确定的，y轴为与所述x轴垂直的轴线，由接收机、转台转动中心和外辐射源构成的平面为xoy平面，z轴为与xoy平面垂直垂直的轴线。
预处理模块3，用于对两种所述窄带调频回波信号进行预处理，获得离散回波信号矩阵。
匹配矩阵确定模块4，用于根据所述离散回波信号矩阵确定信号匹配矩阵。
聚焦成像模块5，用于根据所述离散回波信号矩阵和所述信号匹配矩阵确定得到聚焦成像。
本发明系统还包括：补偿模块6，用于获取目标，并对目标进行运动补偿，获得能够在转台模型下研究的成像目标。
下面对各个模块进行具体分析：
所述构建模块1，用于根据获取的接收机、转台转动中心和外辐射源的位置，构建三维旋转成像模型直角坐标系。
如图3所示，所述三维旋转成像模型直角坐标系是以转台转动中心、外辐射源、接收机所在的平面为xoy平面，并以外辐射源、坐标原点、接收机的位置所构成角的角平分线为x轴，垂直于x轴的方向为y轴，经过o点且垂直于xoy平面的的轴线为z轴。外辐射源和接收机分别位于(rtcosβ,rtsinβ,0)、(rrcosβ,-rrsinβ,0)，其中rt、rr分别为外辐射源和接收机到坐标原点的距离，系统的双基角为2β。假设初始时刻目标上一散射点的三维坐标为(x0,y0,z0)，成像目标先绕z轴转动(偏航)、再绕y轴转动(纵摇)，转台绕z轴转动角速度为ωz，经过时间tz散射点坐标为(x0cosωztz-y0sinωztz,x0sinωztz+y0cosωztz,z0)，如果其绕z轴转动总角度为α，则其最终转动至(x0cosα-y0sinα,x0sinα+y0cosα,z0)。在此基础上，转台绕着y轴开始转动，转动角速度为ωy并重新开始计时，经过时间ty散射点转动至(xy(ty),yy(ty),zy(ty))，其具体形式为：
所述获取模块2，用于利用接收机获取成像目标的两种窄带调频回波信号，具体公式为：
其中，a为窄带调频回波信号的幅度，f为窄带调频回波信号的频率，φ为窄带调频回波信号的初始相位，σ为散射点散射强度，c为光速，λ为波长，rt为外辐射源到坐标原点的距离，rr为接收机到坐标原点的距离，α为绕z轴转动总角度，β为基角，ωz为绕z轴旋转的角速度，tz为绕z轴旋转的时间，ωy为绕y轴旋转的角速度，ty为绕y轴旋转的时间，(x0,y0,z0)为初始时刻目标上一散射点的三维坐标，ry(ty)为成像目标绕y轴旋转的斜距历程，rz(tz)为成像目标绕z轴旋转的斜距历程。
所述预处理模块3，具体包括：
混频和补偿处理单元，用于对两种所述窄带调频回波信号分别进行混频和补偿处理，获得两种补偿后的窄带调频回波信号；
合成处理单元，用于对两种所述补偿后的窄带调频回波信号进行合成处理，获得合成回波信号；
离散采样单元，用于对所述合成回波信号进行离散采样，获得离散回波信号矩阵。
所述匹配矩阵确定模块4，具体包括：
网格划分单元，用于对成像目标所在的立体空间进行网格划分；
离散回波信号矩阵确定单元，用于根据所述离散回波信号矩阵确定网格处散射点的离散回波信号矩阵；
匹配矩阵单元，用于根据所述网格处散射点的离散回波信号矩阵确定信号匹配矩阵。
利用本发明所提方法进行下面两组仿真实验。由前面内容可知，本方法对rt和rr不敏感，所以在下面仿真中没有给出其值。而前面推导过程中的参数α是转台绕z轴转动的角度。不失一般性，仿真过程中散射点的散射强度均取1。
图4为本发明实施例第一聚焦成像仿真图；如图4所示，外辐射源信号频率为1ghz，双基角为0.5πrad，绕z轴转动2.2πrad，绕y轴转动2πrad，对置于原点处的点目标进行成像；图4中的4(a)为x方向成像结果，4(b)为y方向成像结果，4(c)为z方向成像结果，4(d)为xoy面成像结果，4(e)为xoz面成像结果，4(f)yoz面成像结果。
由上可知，在z方向上为贝塞尔函数形式收敛，理论分辨率、峰值旁瓣比分别为0.0764m、-7.9db，其实验结果为其成像结果分别为0.0762m、-8.2db，与理论值吻合。实际上，其偏差是由两种转动之间的耦合作用造成的。从图4(d)到图4(f)可以看出，由于耦合作用的存在，在xoy平面、xoz平面、yox平面的成像结果也出现波动现象。
图5为本发明实施例第二聚焦成像仿真图；图5中的5(a)为实验1聚焦结果，5(b)为实验2聚焦结果，5(c)为实验3聚焦结果，5(d)为实验4聚焦结果，5(e)为实验5聚焦结果，5(f)为实验6聚焦结果；如图5所示，
在长宽高均为10m的立体空间内随机产生10个散射点。然后利用不同的仿真参数对这10个散射点进行聚焦成像。仿真参数如表1所示。仿真结果如图5所示。其中“o”表示散射点的真实位置，而“*”表示成像结果。
表1仿真实验参数
对比图5(a)、图5(b)和图5(c)可以看出，在其他参数相同的情况下，聚焦结果会随着外辐射源信号频率的提升而提升。对比图5(c)和图5(d)可以看出，在其他参数相同的情况下，双基角越小，其聚焦效果越好。对比图5(c)和图5(e)可以看出，当目标绕着z轴、y轴转动不到一周时，会影响其成像效果，转动的角度越小，成像效果越差。在实际情况下，成像目标一般为非合作对象，其转动角度一般较小，所以对小转动角成像的研究就尤为重要。从图5(f)可以看出，当外辐射源信号较高时，即使在双基角较大、转动角度较小时，也能取得较好的成像效果。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。对于实施例公开的系统而言，由于其与实施例公开的方法相对应，所以描述的比较简单，相关之处参见方法部分说明即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。