非均匀旋转目标运动参数的快速估计方法与流程

本发明属于雷达信号处理技术领域，尤其涉及一种非均匀旋转目标运动参数的快速估计方法，用于消除传播误差的影响并同时降低计算复杂度。

背景技术：

逆合成孔径雷达(isar)可以通过地基雷达获取空域非合作运动目标的高分辨雷达图像，实现对目标的观测成像与识别监视。在许多应用中，目标总是具有高机动性和不均匀旋转。对于典型的非均匀旋转目标，由于目标轮廓上的每一个散射点均具有时变性，isar图像质量会出现明显退化。

为了解决这个问题，现有的距离瞬时多普勒(rid)算法，作为一种参数化方法，rid通常是对数据单个距离单元内的回波信号进行建模，随后对模型中的参数进行估计，重建目标isar图像。但该方法会引入较为明显的传递误差，造成参数估计精度与成像质量下降。此外，该方法通常需要数次积累与搜索来估计参数的值，这会导致算法运算量巨大难以在实际工程中应用。

由于描述机动目标的优越性，二次频率调制(qfm)信号通常用于对接收信号回波建模并在isar成像中进行运动参数估计。基于该模型，开发了许多isar参数估计算法，如乘积立方相函数(pcpf)，积分立方相函数(icpf)和相干积分广义立方相函数(cigcpf)等。这些算法在估计各次项参数时，往往会不可避免地引入传递误差，这会对最终的参数估计精度造成一定影响；而且，这些方法总是需要多次集成或搜索，使实现过程的计算量巨大，限制其实际应用。而时频分析类(tfd)方法利用信号在时频分布平面的几何关系信息，能够独立地估计isar信号的各项运动参数，有效降低了传递误差对估计精度的影响，因此在isar成像处理中得到了广泛的应用。然而，当目标进行非均匀旋转时，由于高阶项的存在，传统基于tfd的方法的参数估计精度与成像质量都会出现明显下降。

技术实现要素：

针对上述问题，本发明的目的是提出一种非均匀旋转目标运动参数的快速估计方法。本发明能够显著降低目标运动参数估计过程的运算量，同时，有效避免了传递误差对各项运动参数估计的影响，保证了最终目标参数的估计精度与isar成像质量。

本发明的技术原理：将单个距离单元内的回波信号建模为qfm信号。构造相关函数，将qfm信号变换至时间-调频率分布(tcd)平面内，并对二次调频率与时间-调频率轨迹斜率之间的关系进行分析。考虑到搜索会引入较大的计算量，提出了基于tcd的最小二乘radon变换(rt变换)来估计线性调频率的导数，即二次调频率，利用数次rt变换的结果即可实现对二次调频率的估计。考虑到交叉项和低信噪比(snr)可能降低估计精度，提出了基于最小熵准则的加权最小二乘radon变换(weightedleastsquareradontransform，wlsrt-tcd)估计方法，能够有效改善估计精度与低信噪比环境下算法的稳健性。类似地，提出基于tfd的加权最小二乘radon变换(wlsrt-tfd)来估计线性调频率，并且通过滤波来消除交叉项干扰。随后根据距离-多普勒分布平面中的几何信息，通过wlsrt方法直接估计出中心频率。

为了达到上述目的，本发明采用以下技术方案予以解决。

非均匀旋转目标运动参数的快速估计方法，包括以下步骤：

步骤1，雷达发射线性调频(lfm)信号，即为发射信号，并接收回波信号，对回波信号进行基带调制和距离压缩后，得到距离压缩后的回波信号，将距离压缩后的回波信号转换为二次调频信号形式，得到每个距离单元内的方位回波信号，即为二次调频率形式的回波信号s(ta)。

步骤2，构建二次调频率形式的回波信号s(ta)的自相关函数r(ta，τ)，对自相关函数r(ta，τ)进行快速傅里叶变换，得到时间-调频率分布平面内的回波信号ws(ta，γ)，即将回波信号s(ta)转换到时间-调频率分布平面。

步骤3，将二次调频率形式的回波信号s(ta)中的二次调频率去除，再进行维纳分布(wigner-villedistribution，wvd)，得到维纳分布形式的回波信号即将回波信号s(ta)转换到时频分布平面。

步骤4，分别对时间-调频率分布平面内的回波信号ws(ta，γ)和维纳分布形式的回波信号进行radon变换，并采用加权最小二乘法估计出对应的二次调频率轨迹斜率角θ和线性调频率轨迹的斜率角α。

步骤5，根据目标在时间-调频率分布平面内轨迹的斜率角、目标在时频分布平面内轨迹的斜率角、目标在距离-多普勒分布平面内轨迹的斜率角，分别估计出对应的目标运动参数，即目标的二次调频率en，3、目标的线性调频率en，2和目标的中心频率en，1。

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

(1)本发明通过基于时间-调频率分布的最小二乘radon变换估计二次调频率，降低了计算负担，提高了非均匀旋转目标成像的精度和鲁棒性。

(2)本发明通过基于时频分布的加权最小二乘radon变换估计线性调频率和基于距离-多普勒分布平面的加权最小二乘radon变换估计中心频率，消除了非均匀旋转目标成像的交叉项干扰，同时消除了传播误差的影响。

附图说明

下面结合附图和具体实施例对本发明做进一步详细说明。

图1是本发明的使用场景图；

图2是本发明的实现流程图；

图3是radon变换对应的几何关系图；

图4是本发明在模拟环境下对非均匀旋转目标的成像效果图；

图5是本发明在真实环境下对非均匀旋转目标的成像效果图。

具体实施方式

下面结合附图对本发明的实施例及效果作进一步详细描述。

参照图1，将目标的旋转中心作为坐标原点o，构造笛卡尔坐标系o-xyz。设定雷达视线方向的单位向量目标的三维角速度向量选定目标观测成像的投影平面。其中，目标的三维角速度向量可分解为：平行于雷达视线方向的分量和垂直于雷达视线方向的分量且有由于水平分量不会产生任何旋转运动，因此其对信号回波的多普勒没有贡献，这里我们仅需分析讨论垂直分量的作用影响。

参照图2，本发明的实现步骤如下：

步骤1，雷达发射线性调频(lfm)信号，即为发射信号，并对回波信号进行基带调制和距离压缩后，得到距离压缩后的回波信号，将距离压缩后的回波信号转换为二次调频信号形式，得到每个距离单元内的方位回波信号，即为二次调频率形式的回波信号s(ta)；

子步骤1.1，设定目标上的散射点p的坐标为则其多普勒频率fd为：

其中，表示叉乘运算，⊙表示内积运算。表示从原点到散射点p的方向矢量，vr为散射点p的平动速度，λ是发射信号的波长。

由于逆合成孔径雷达成像中合成孔径时间短，通常假设成像投影平面在成像周期内保持不变。因此散射点p的多普勒频率简化式为：

子步骤1.2，由于散射点p的平动速度vr、散射点p垂直于雷达视线方向的角速度分量ωpx、ωpy、ωpz都是与非均匀旋转目标的时间变量，将其通过泰勒展开扩展到二阶，得到对应的泰勒展开式为：

其中，ta是方位慢时间，v0表示初始径向速度，a0表示初始径向加速度，γ0表示初始径向加加速度；目标的三维角速度向量在垂直于雷达视线方向的分量投影至笛卡尔坐标系中的三条坐标轴上的分量的常系数分别为μx，μy，μz，且目标的三维角速度向量在垂直于雷达视线方向的分量投影至笛卡尔坐标系中的三条坐标轴上的分量的一阶系数分别为αx，αy，αz，且目标的三维角速度向量在垂直于雷达视线方向的分量投影至笛卡尔坐标系中的三条坐标轴上的分量的二阶系数分别为βx，βy，βz，且

子步骤1.3，将子步骤1.2中的泰勒展开式代入散射点p的多普勒频率简化式，得到散射点p的二阶多普勒频率表达式为：

子步骤1.4，计算散射点p到雷达的距离，即为目标斜距项rp(ta)：

其中，r0表示在初始时间t0时散射点到雷达的斜距，γ0表示初始加加速度，表示散射点p的等效位置向量，且

子步骤1.5，距离压缩后的回波信号的表达式为：

其中，an为第n个散射点的幅度常数项，tr为距离快时间，c为光速，br为发射信号带宽，ta表示散射点积累时间，n(tr，ta)表示方差为的加性复高斯白噪声项，sinc[·]表示辛格函数，且有sinc(t)＝sin(πt)/πt，rect[·]表示为矩形脉冲，

子步骤1.6，目标斜距信息项rp(ta)会引起包括目标平动、跨距离单元走动以及多普勒频率偏移在内的距离徙动问题，本方法主要讨论对多普勒频率偏移信息的估计与补偿。由于距离徙动的存在会对后续目标运动参数估计的精度产生明显的影响，因此利用包络对齐与初相校正首先对信号进行平动补偿。设定经过平动补偿后，目标的所有散射点均位于正确的距离单元之内，且每个距离单元内共有n个散射点，将每个散射点到雷达的距离代入距离压缩后的回波信号的表达式，并合并同类项，得到单个距离单元内的方位回波信号，即为二次调频率形式的回波信号s(ta)，其表达式为：

其中，en，0表示单个距离单元内的方位回波信号的初始相位，且en，0＝r0；en，1表示单个距离单元内的方位回波信号的中心频率，且en，2表示单个距离单元内的方位回波信号的线性调频率，且对应目标的一阶运动参数；en，3表示单个距离单元内的方位回波信号的二次调频率，且对应目标的二阶运动参数；n(ta)表示合并同类项后的n(tr，ta)。

步骤2，构建二次调频率形式的回波信号s(ta)的自相关函数r(ta，τ)，对自相关函数r(ta，τ)进行快速傅里叶变换，得到时间-调频率分布平面内的回波信号ws(ta，γ)，即将回波信号转换到时间-调频率分布平面。

子步骤2.1，构建二次调频率形式的回波信号s(ta)的自相关函数r(ta，τ)：

其中，*是取共轭操作，rc，n(ta，τ)表示交叉项，其由多组待估计的分量和加性噪声引起；τ表示滞后时间，rect[·]为矩形脉冲，

子步骤2.2，沿滞后时间τ对自相关函数r(ta，τ)进行快速傅里叶变换，得到时间-调频率分布平面内的回波信号ws(ta，γ)：

其中，tcdc，n(ta，γ)是快速傅里叶变换后的交叉项；γ表示发射信号的调频率，从上式可以看出，由于时间-调频率分布平面内的回波信号的表达式为关于方位慢时间的一个线性函数，且该线性函数的斜率仅由二次调频率所决定，而不受线性项与二次项的影响，因此在一定程度上，能够有效避免传递误差对各阶运动参数估计精度的影响。

步骤3，将二次调频率形式的回波信号s(ta)中的二次调频率去除，得到去除二次调频率的回波信号

对去除二次调频率的回波信号进行维纳分布(wigner-villedistribution，wvd)，得到维纳分布形式的回波信号为：

其中，fτ是滞后时间τ的频率变量，即滞后时间τ对应的频率，wvdc，n(ta，fd)表示维纳分布的交叉项。从上式可以看出线性函数的斜率包含线性调频率，中心频率项引起的传播误差不会影响该斜率。

步骤4，分别对时间-调频率分布平面内的回波信号ws(ta，γ)和维纳分布形式的回波信号进行radon变换，并采用加权最小二乘法估计出对应的二次调频率轨迹斜率角θ和线性调频率轨迹的斜率角α。

子步骤4.1，在时间-调频率分布平面，对步骤2得到时间-调频率分布平面内的回波信号ws(ta，γ)进行距离为ρ，角度为的radon变换，得到radon变换后的时间-调频率分布平面回波信号

其中，δ(·)表示冲击函数。

在时频分布平面，对维纳分布形式的回波信号的进行距离为ξ，角度为β的radon变换，得到radon变换后的时频分布平面回波信号wr(ξ，β)：

子步骤4.2，采用加权最小二乘法，分别在时间-调频率分布平面和时频分布平面上估计目标在时间-调频率分布平面内轨迹的斜率角θ和目标在时频分布平面内轨迹的斜率角α。

子步骤4.2.1，采用加权最小二乘法，在时间-调频率分布平面上估计目标在时间-调频率分布平面内轨迹的斜率角θ。

子步骤4.2.1.1，在时间-调频率分布平面，设定表示任一角度下的radon变换结果，即radon变换后的时间-调频率分布平面回波信号轨迹在角度上的投影，即归一化投影长度(lnp)；设定目标在时间-调频率分布平面内轨迹的斜率角为θ，根据图3，利用几何关系，得到几何关系式：

其中，lθ表示目标在时间-调频率分布平面的归一化长度，表示当radon变换角度为时，回波信号轨迹在角度上的投影长度。

子步骤4.2.1.2，设定当radon变换角度为时，引入的测量误差为ε，则得到时间-调频率分布平面内的几何模型为：

其中，表示当radon变换角度为时，实际测量的回波信号轨迹在角度上的投影长度。

子步骤4.2.1.3，将时间-调频率分布平面内的几何模型中的一次radon变换扩展为n次radon变换，则时间-调频率分布平面内的几何模型的扩展式为：

其中，k＝1，2…，k，k为正整数；表示第k个radon变换的角度，εk表示第k个radon变换的测量误差，lk表示第k个radon变换的归一化投影长度；将几何模型的扩展式写成矩阵方程为：

ax+ε＝l；

其中，

采用最小二乘法求解所述矩阵方程，得到估计的目标运动参数矩阵x为：

x＝[lθsinθlθcosθ]^t＝[a^ta]^-1a^tl；

子步骤4.2.1.4，采用radon变换的熵作为权重，采用加权最小二乘法求解矩阵方程，得到精确估计的目标运动参数矩阵，来提高参数估计精度。

首先，设定wr(m，b)表示旋转角为b的第m个距离采样点的radon变换结果，计算radon变换的熵i：

其中，m表示距离采样点数。

其次，将radon变换的熵作为权重值，得到权重矩阵w为：

w＝diag[i]；

其中，diag[·]是对角化操作。

最后，采用radon变换的熵i表示几何模型中测量误差的模型，即加权最小二乘法求解矩阵方程，得到精确估计的目标运动参数矩阵x′：

得到目标在时间-调频率分布平面内轨迹的斜率角的正切为：

tanθ＝lθsinθ/lθcosθ＝sinθ/cosθ；

进而得到二次调频率轨迹的斜率角θ，即为二次调频率轨迹与x轴的夹角。

子步骤4.2.2，在时频分布平面，设定为任一角度为β0下的radon变换结果，即radon变换后的时频分布平面回波信号wr(ξ，β)在角度β0上的投影；重复子步骤4.2.1.1-子步骤4.2.1.4，采用加权最小二乘法，估计目标在时频分布平面内轨迹的斜率角的正切tanα，进而得到线性调频率轨迹的斜率角α，即线性调频率轨迹与x轴的夹角。

步骤5，根据目标在时间-调频率分布平面内轨迹的斜率角、目标在时频分布平面内轨迹的斜率角、目标在距离-多普勒分布平面内轨迹的斜率角，分别估计出对应的目标运动参数en，3、en，2、en，1。

子步骤5.1，在时间-调频率分布平面，估计目标的二次调频率en，3。

从的表达式可以看出，将二次频率调制信号变换至时间-调频率分布平面后，其变换轨迹是一关于方位慢时间ta的线性函数，由于该线性函数的斜率仅由qfm信号的二次调频率所决定，而不受线性项与二次项的影响，因此在一定程度上，能够有效避免传递误差对各阶次运动参数估计精度的影响。

其具体步骤为：在积分时间ta内，时间-调频率分布平面中，机动目标在x轴上的投影长度δx3和机动目标在y轴上的投影长度δy3分别表示为：

其中，prf是脉冲重复频率，fsγ是时间-调频率域的采样频率；

则目标在时间-调频率分布平面内轨迹的斜率角的正切表示为：

结合子步骤4.2得到的目标在时间-调频率分布平面内轨迹的斜率角的正切tanθ，从而得到二次调频率en，3：

en，3＝prf·fsγ/tanθ；

子步骤5.2，在时频分布平面，估计目标的线性调频率en，2。

从的表达式可以看出其线性函数的斜率包含调频率en，2，中心频率项引起的传播误差不会影响该斜率。

其具体步骤为：在积分时间ta内，时频分布(tfd)平面中机动目标在x轴上的投影长度δx2和在y轴上的投影长度δx2分别表示为：

其中，表示多普勒频率的采样率。

则目标在时频分布平面内轨迹的斜率角的正切表示为：

结合子步骤4.3得到的线性调频率轨迹的斜率角α，从而得到机动目标的线性调频率：

子步骤5.3，将二次调频率形式的回波信号s(ta)中的二次调频率和线性调频率去除后，回波信号为关于方位慢时间的线性函数，其斜率包含中心频率en，1。

设定中心频率轨迹与x轴的夹角为ψ，在积分时间ta内，距离-多普勒平面中，机动目标在x轴上的投影长度δx1和y轴上的投影长度δy1分别表示为：

其中，δr表示距离分辨率，δr＝c/2·fsr，fsr是距离采样频率。

根据距离-多普勒分布平面中的几何关系，可以得到：

得到机动目标的中心频率en，1为：

仿真实验

本发明的效果可通过以下仿真实验进一步说明。

仿真实验1，验证在模拟环境下本发明方法对非均匀旋转目标的成像效果。

(1)仿真参数：

模拟机动目标模型是一架有39个散射点的飞机。目标与雷达的初始距离为3.5km，目标的平移速度为30m/s。

参数设置如表1：

表1系统仿真参数

(2)仿真内容：

仿真1，在上述仿真参数下，利用本发明方法，在信噪比snr＝2db，利用12个rt结果通过不同方法估计参数，如图4所示。

由图4可以看出，图4(a)为基于pcpf方法的isar成像结果图，图4(b)为基于cigcpf方法的isar成像结果图，图4(c)为基于二维搜索方法的isar成像结果图，图4(d)为本发明方法的的isar成像结果图，从图4可以看出，pcpf方法、cigcpf方法、二维搜索方法和本发明方法都可以准确估计目标的运动参数并获得聚焦良好的isar图像。与本发明方法相比，二维搜索方法、pcpf方法、cigcpf方法获得图像性能稍有劣化。在低信噪比的情况下，pcpf方法的估计精度比其他方法稍差，随着信噪比的降低，精度会严重恶化，对于cigcpf，尽管它具有很好的噪声鲁棒性，但存在相当大的传播误差影响估计精度。

仿真实验2，验证真实环境下本发明方法对非均匀旋转目标的成像效果。

(1)仿真参数：

表2真实系统参数

(2)仿真内容：

在上述仿真参数下，采用本发明方法，实验中观察到非合作飞机，snr约为-4db，利用实际数据来验证本发明方法的可行性和有效性，结果如图5所示。

图5(a)微微真实环境下pcpf方法获得的图像，图5(b)为真实环境下cigcpf方法获得图像，图5(c)为真实环境下二维搜索方法获得的图像，图5(d)为真实环境下本发明方法获得的图像。由图5(a)可以看出，pcpf方法获得的图像具有严重的散焦，并且由于其对噪声的鲁棒性差而在低信噪比情况下遭受更多的散射损失。由图5(b)可以看出，cigcpf方法通过相干积分提高了对噪声的鲁棒性，但其本身引起相当大的传播误差，仍然会引起一些散焦。由图5(c)可以看出，二维搜索方法可以获得准确的参数估计和良好的聚焦图像，但其计算负担繁重限制了其实际应用。由图5(d)可以看出，本发明方法获得了准确的参数估计和良好的聚焦图像。

仿真实验3.计算仿真实验1中4种方法的isar成像的熵和所花费的成像时间，结果如表3所示。

表3isar成像的参数性能

由表3可以看出，二维搜索方法具有最低的熵，但具有巨大的计算负担。本发明方法的计算时间远远小于其他方法，且具有较低的熵，说明本发明方法对目标的运动参数的估计精度较高，同时大大减少了成像所需时间，易于工程实现。

本领域普通技术人员可以理解：实现上述方法实施例的全部或部分步骤可以通过程序指令相关的硬件来完成，前述的程序可以存储于一计算机可读取存储介质中，该程序在执行时，执行包括上述方法实施例的步骤；而前述的存储介质包括：rom、ram、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。