CN103091674A - 基于hrrp序列的空间目标高分辨成像方法 - Google Patents

Abstract

本发明公开了一种基于HRRP序列的空间目标高分辨成像方法，其步骤包括：(1)雷达录取ISAR回波；(2)高速运动补偿；(3)获取高质量HRRP；(4)散射点航迹关联；(5)高分辨成像。本发明首先通过稀疏信号重构获得高质量HRRP序列，之后采用有效的航迹关联方法从HRRP序列中生成散射点航迹矩阵，进而通过带约束条件的矩阵分解实现高分辨成像。该方法克服了基于分段伪Keystone变换方法对PRF要求高、对目标运动模型要求严格、对高速旋转目标成像时精确平动补偿和初相校正难度大、运算量高等缺点，具有在低PRF甚至方位多普勒模糊条件下能够进行高分辨成像、对目标运动模型具有普适性、不需平动补偿、效率高、图像聚焦良好等优点。

Description

基于HRRP序列的空间目标高分辨成像方法

技术领域

本发明属于信号处理技术领域，更进一步涉及雷达成像领域中的基于高分辨一维距离像(high-resolution range profiles，HRRP)序列的空间目标高分辨成像方法。本发明可以在低重频条件下，有效地对空间目标进行高分辨成像，避免方位欠采样、平动量和随机初相误差对图像质量的影响，为后续目标识别提供有力保障。

背景技术

为获得对空间目标散射点分布特征更为精细和准确的描述，目前对雷达图像分辨率的要求在不断提高。其中距离向高分辨依靠发射大的时宽-带宽积信号获得；而方位向高分辨则依靠成像时间内目标相对于雷达视线运动所引起的多普勒变化获得。为满足距离向高分辨成像要求，可以通过信号处理方法实现多波段超宽带合成，或通过系统硬件升级等方法使雷达系统发射超宽带信号。但是在获得距离向高分辨率的同时，若仍采用传统的成像方法获得方位向高分辨率时会遇到新的问题：首先，距离分辨率提高的同时会使散射点更易发生越距离单元走动，从而限制了距离-瞬时多普勒、时频分析等成像方法的应用；其次，由于电磁波受大气扰动的影响，回波相位中会包含时变的随机初相误差。当散射点存在越距离单元走动时，对该误差的校正及平动补偿会变得更加困难；最后，由于地基雷达发射功率有限，当对相距较远的空间目标进行探测时，需要降低雷达脉冲重复频率以保证雷达接收机对信噪比的要求。在这种情况下，对于旋转频率较高的机动目标，其回波会出现方位欠采样，混迭的方位频谱将进一步增加对空间目标高分辨成像的难度。因此，有必要研究在低重频条件下的新的高分辨成像方法。

当采用宽带雷达对目标进行高分辨观测时，根据投影定理，目标的高分辨一维距离像(HRRP：high-resolution range profiles)将反映出散射点分布沿雷达视线方向的投影。当雷达视线连续变化时，不同观测时刻对应的HRRP将形成HRRP序列。当散射点不存在越距离单元走动时，HRRP序列仅描述散射点沿距离向的分布；当采用超宽带对机动目标进行观测时，由于发生越距离单元走动，实HRRP序列会反映出散射点在观测时间内的运动轨迹，从而包含散射点的二维或三维分布信息。因此，通过合适的算法设计能够实现基于HRRP序列的目标高分辨成像。同传统成像方法相比，这类方法不会受到方位欠采样和随机初相误差的影响，同时降低了成像对雷达PRF的要求。

Kai Huo等人在文献“A novel imaging method for fast rotating targets based on thesegmental pseudo Keystone transform”(IEEE Trans.Geosci.Remote Sens.，vol.49，no.4，pp.1464-1472，Apr.2011)中提出基于分段伪Keystone变换的快速旋转目标二维成像方法，该方法在平动可被准确补偿，并且目标的旋转频率较高且在观测时间内恒定前提下可以精确估计出目标的位置。但是，该方法仍存在的不足是，当散射点存在严重越距离单元走动而使传统相邻相关平动补偿方法失效时，没有考虑有效的平动补偿方法；当目标存在变加速旋转、或雷达进行重频抖动(即脉冲重复频率随时间变化)以实现ECCM(electronic counter-countermeasures)时，该方法假设的正弦模型已经不再满足，从而导致图像模糊；至少需要接收到目标一个旋转周期内的回波才能得到低旁瓣、聚焦良好的图像，当回波次数较少时，旁瓣很高；运算复杂度较高。

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种基于HRRP序列的空间目标高分辨成像方法。该方法弥补了基于分段伪Keystone变换方法具有对PRF要求高、对目标运动模型要求严格、对高速旋转目标成像时精确平动补偿和初相校正难度大、参数搜索运算量高的不足，不需要对回波进行精确的平动补偿和初相校正，在低重频条件下仅利用散射点在距离-慢时间平面上的航迹特征进行非参数化成像即可得到目标聚焦良好的高分辨图像。

实现本发明的基本思路是：首先，采用稀疏信号重构方法获得目标高质量的HRRP序列；其次，根据散射点的幅度和运动特征，采用航迹关联方法从HRRP序列中提取出目标强散射点对应的轨迹并生成散射点航迹矩阵；最后，从航迹矩阵中求解出散射点坐标并得到高分辨图像。

本发明的具体步骤如下：

(1)雷达录取ISAR回波；

(2)高速运动补偿

2a)将参考距离为雷达到场景中心距离，频率、调频率与雷达发射信号相同的线性调频信号作为参考信号；

2b)将参考信号取共轭后与雷达录取的ISAR回波相乘得到解线调频处理结果；

2c)按照下式构造相位补偿函数：

其中，

为相位补偿函数，exp(·)为指数函数运算符，j为虚数单位，c为光速，v为测定的目标速度，t为方位维的采样时间，R(t)为t时刻场景中心到雷达的参考距离，f₁为信号载频，γ为调频率，f₂为距离频率；

2d)将相位补偿函数与步骤2b)得到的解线调频处理结果相乘得到高速运动补偿后回波；

(3)获取高质量HRRP

3a)按照下式构造过冗余字典：

$\mathrm{\Φ}=[\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(1\right)\right),\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(2\right)\right),\·\·\·,\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(N\right)\right)]$

其中，Ф为包含N个基向量的过冗余字典，N取为距离频率点数的整数倍，exp(·)为指数函数运算符，j为虚数单位，c为光速，f₂为距离频率，

f₁为载频，γ为调频率，

为距离维采样时间，R为参考距离，R_Δ(1)，R_Δ(2)，…R_Δ(N)为构造的N个散射点斜距，

W为雷达观测场景长度，相邻两个距离采样点的距离步进量满足关系

M为距离频率点数，B为信号带宽；

3b)令m＝1，其中，m表示方位单元序号；

3c)将步骤2d)得到的高速运动补偿后回波的第m个方位单元的复回波作为残余信号，并生成一个空矩阵Ф₁；

3d)用残余信号向步骤3a)得到的过冗余字典中的基向量作投影，并且计算内积；

3e)将内积最大的基向量记录在步骤3c)生成的矩阵Ф₁的第m列中，按照投影方法计算内积最大的基向量对应的投影系数，记录得到的投影系数，同时将内积最大的基向量从过冗余字典中剔除；

3f)从残余信号中减去矩阵Ф₁与记录的投影系数的乘积，得到更新的残余信号；

3g)重复步骤3d)至步骤3f)的操作，直至更新的残余信号能量低于噪声门限，将步骤3e)中记录的投影系数取模值作为对应基向量所代表斜距处散射点的幅度，以基向量对应的散射点斜距为横坐标，幅度为纵坐标绘制HRRP；

3h)将m递增1个序号，重复步骤3c)至步骤3g)的操作，直至绘制完所有方位单元对应的HRRP；

(4)散射点航迹关联

4a)令k＝1，其中，k表示方位单元序号，生成一个大小为N₁×M的空的航迹矩阵，其中，N₁为第k个方位单元对应的HRRP中峰值点的个数，M表示散射点回波次数，将第k个方位单元对应的HRRP中所有峰值点对应的位置坐标和幅度记录在航迹矩阵第k列中，将k递增一个序号；

4b)在第k个方位单元对应的HRRP中寻找峰值点，并将所有峰值点对应的位置坐标、幅度记录在第k时刻的位置幅度集合中，根据小转角内散射点幅度近似不变且位置接近准则，将第k时刻位置幅度集合中的元素与航迹矩阵第k-1列中元素进行关联，将关联得到的第k时刻散射点位置坐标和幅度记录在航迹矩阵第k列；

4c)将k递增一个序号，执行步骤4b)；

4d)分别将各散射点在航迹矩阵中第k-1列记录的位置坐标与第k-2列记录的位置坐标相减，并除以脉冲重复周期得到速度，将航迹矩阵中第k列记录的位置坐标减去第k-1列记录的位置坐标的两倍，再加上第k-2列记录的位置坐标，将计算结果除以脉冲重复周期的平方得到加速度；

4e)用速度乘以脉冲重复周期并加上加速度与脉冲重复周期平方乘积的0.5倍，求得距离步进量；

4f)分别将各散射点在航迹矩阵中第k列记录的位置坐标、步骤4d)得到的速度与加速度作为第k时刻散射点状态向量；

4g)根据Kalman预测准则，对第k+1、k+2及k+3时刻的状态向量进行预测，并将航迹矩阵中记录的对应散射点幅度求均值作为预测幅度值；

4h)在第k+1个方位单元对应的HRRP中，分别以各散射点在航迹矩阵中第k列记录的位置坐标为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+1时刻的位置幅度集合中；

4i)在第k+2个方位单元对应的HRRP中，分别以各散射点在第k+1时刻的位置幅度集合中所有位置坐标的均值为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+2时刻的位置幅度集合中；

4j)在第k+3个方位单元对应的HRRP中，分别以各散射点在第k+2时刻的位置幅度集合中所有位置坐标的均值为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+3时刻的位置幅度集合中；

4k)根据归一化的最近邻准则函数对散射点航迹进行关联，求得散射点在第k+1时刻的真实位置坐标和幅度，并记录在航迹矩阵第k+1列中；

41)将步骤4g)预测得到的k+1时刻状态向量中的速度和加速度按照步骤4e)进行计算，得到距离步进量，将k递增一个序号，重复步骤4g)至步骤4k)，直至完成散射点在所有方位单元对应的HRRP中的位置关联，从航迹矩阵中抽取散射点的位置坐标并且剔除对应幅度值，得到更新的航迹矩阵；

(5)高分辨成像

5a)将航迹矩阵中每个元素都减去该元素所在行所有元素的均值，得到平动校正后的航迹矩阵；

5b)对平动校正后的航迹矩阵进行带约束的矩阵分解得到散射点分布矩阵和雷达视线矩阵；

5c)用散射点分布矩阵重构目标图像。

与现有技术相比，本发明具有以下优点。

第一，本发明通过采用目标回波对应的实HRRP序列进行高分辨成像，克服了现有技术中对雷达PRF要求高，在方位欠采样条件下无法得到聚焦良好的高分辨图像等不足，使得本发明具有在低PRF甚至方位多普勒模糊条件下能够进行高分辨成像。

第二，本发明通过采用航迹关联方法从HRRP序列中提取出目标强散射点对应的轨迹并生成散射点航迹矩阵，通过对航迹矩阵的简单处理实现对目标的平动补偿，克服了现有技术中对高速旋转目标成像时精确平动补偿和初相校正难度大的不足，使得本发明在回波存在剩余平动量和随机初相误差的条件下能够得到目标聚焦良好的高分辨图像。

第三，本发明通过带约束条件的矩阵分解，采用非参数方法得到目标散射点坐标和等效雷达视线的分布矩阵，克服了现有技术中对目标运动模型要求严格、高质量成像时要求目标回波数较多，参数搜索时运算量较大的不足，使得本发明具有对目标运动模型具有稳健性、成像要求观测回波次数少、运算效率高、图像聚焦良好的优点。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图1，对本发明具体实施方式作进一步的详细描述：

步骤1，雷达录取目标的逆合成孔径ISAR回波，雷达以脉冲重复频率发射并接收脉冲，得到以距离为行向量以方位为列向量的ISAR回波。

步骤2，高速运动补偿

2a)将参考距离为雷达到场景中心距离，频率、调频率与雷达发射信号相同的线性调频信号作为参考信号，其中参考信号如下：

$S(\hat{t},t_1)=\mathrm{rect}\left(\frac{\hat{t}-2R^{\′}/c}{T}\right)e^{j2\mathrm{\π}(f(t-2R^{\′}/c)+\frac{1}{2}\mathrm{\γ}{(\hat{t}-2R^{\′}/c)}^2)}$

其中，为参考信号，

为距离维采样时间，t₁为方位维采样时间，rect(·)为矩形窗函数， $\mathrm{rect}\left(u\right)=\left\{\begin{array}{cc}1& \left|u\right|\≤\frac{1}{2}\\ 0& \left|u\right|>\frac{1}{2}\end{array}\right.,$ R′为参考距离，T为参考信号的脉宽，j为虚数，f为中心频率，t为距离维采样时间与方位维采样时间之和，c为光速，γ为调频率。

2b)将参考信号取共轭后与雷达录取的ISAR回波相乘得到解线调频处理结果；

2c)按照下式构造相位补偿函数：

其中，

为相位补偿函数，exp(·)为指数函数运算符，j为虚数单位，c为光速，v为测定的目标速度，t为方位维的采样时间，R(t)为t时刻场景中心到雷达的参考距离，f₁为信号载频，γ为调频率，f₂为距离频率。

2d)将相位补偿函数与步骤2b)得到的解线调频处理结果相乘得到高速运动补偿后回波。

步骤3，获取高质量HRRP

3a)按照下式构造过冗余字典：

$\mathrm{\Φ}=[\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(1\right)\right),\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(2\right)\right),\·\·\·,\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(N\right)\right)]$

其中，Ф为包含N个基向量的过冗余字典，N取为距离频率点数的整数倍，exp(·)为指数函数运算符，j为虚数单位，c为光速，f₂为距离频率，

f₁为载频，γ为调频率，为距离维采样时间，R为参考距离，R_Δ(1)，R_Δ(2)，…R_Δ(N)为构造的N个散射点斜距，

W为雷达观测场景长度，相邻两个距离采样点的距离步进量满足关系

M为距离频率点数，B为信号带宽。

3b)令m＝1，其中，m表示方位单元的序号；

3c)将步骤2d)得到的高速运动补偿后回波的第m个方位单元的复回波作为残余信号，并生成一个空矩阵Ф₁；

3d)用残余信号向步骤3a)得到的过冗余字典中的基向量作投影，并且计算内积；

3e)将内积最大的基向量记录在步骤3c)生成的矩阵Ф₁的第m列中，按照投影方法计算内积最大的基向量对应的投影系数，记录得到的投影系数，同时将内积最大的基向量从过冗余字典中剔除。其中，投影方法可描述为如下表达式：

$\mathrm{\α}={\left({\mathrm{\Φ}}_1^H{\mathrm{\Φ}}_1\right)}^{-1}{\mathrm{\Φ}}_1^{H}g$

其中，α为投影系数，Ф₁为记录最大投影系数对应基向量的矩阵，

为Ф₁的共轭转置，(·)^-1为矩阵求逆运算，g为残余信号。

3f)从残余信号中减去矩阵Ф₁与记录的投影系数的乘积，得到更新的残余信号；

3g)重复步骤3d)至步骤3f)的操作，直至更新的残余信号能量低于噪声门限，将步骤3e)中记录的投影系数取模值作为对应基向量所代表斜距处散射点的幅度，以基向量对应的散射点斜距为横坐标，幅度为纵坐标绘制HRRP。其中噪声门限可以设定为对步骤3c)得到的残余信号进行逆傅里叶变换并取模值，将取模值后的图像中远离峰值点的其余点对应的平均功率作为噪声门限。

3h)将m递增1个序号，重复步骤3c)至步骤3g)的操作，直至绘制完所有方位单元对应的HRRP。

步骤4，散射点航迹关联

4a)令k＝1，其中，k表示方位单元的序号，生成一个大小为N₁×M的空的航迹矩阵，其中，N₁为第k个方位单元的HRRP中峰值点的个数，M表示散射点回波次数，将第k个方位单元的HRRP中所有峰值点对应的位置坐标、幅度记录在航迹矩阵的第k列中，将k递增一个序号；

4b)在第k个方位单元的HRRP中寻找峰值点，并将所有峰值点对应的位置坐标、幅度记录在第k时刻的位置幅度集合中，根据小转角内散射点幅度近似不变且位置接近准则，将第k时刻位置幅度集合中的元素与航迹矩阵的第k-1列中元素进行关联，将关联得到的第k时刻散射点的位置坐标及幅度记录在航迹矩阵的第k列。其中小转角内散射点幅度近似不变位置接近准则为当前时刻散射点幅度应该满足如下等式：

A_k＝arg min||A′_k-A_k-1||₂

其中，A_k为当前散射点在第k时刻对应的最优幅度值，argmin(·)为求使目标函数取最小值时对应变量值的运算，||·||₂求为2范数，A′_k为第k时刻的位置幅度集合中当前散射点对应的所有幅度值，A_k-1为航迹矩阵中第k-1列中当前散射点对应的幅度值。

4c)将k递增一个序号，执行步骤4b)；

4d)分别将各散射点在航迹矩阵中第k-1列记录的位置坐标与第k-2列记录的位置坐标相减，并除以脉冲重复周期得到速度，将航迹矩阵中第k列记录的位置坐标减去第k-1列记录的位置坐标的两倍，再加上第k-2列记录的位置坐标，将计算结果除以脉冲重复周期的平方得到加速度；

4e)用速度乘以脉冲重复周期并加上加速度与脉冲重复周期平方乘积的0.5倍，求得距离步进量；

4f)分别将各散射点在航迹矩阵中第k列记录的位置坐标、步骤4d)得到的速度及加速度作为第k时刻散射点状态向量；

4g)根据Kalman预测准则，对第k+1、k+2及k+3时刻的状态向量进行预测，并将航迹矩阵中记录的对应散射点幅度求均值作为预测幅度值。其中，Kalman预测准则如下：

${\hat{x}}_{k+1/k}=\mathrm{\Φ}{\hat{x}}_{k/k-1}+{\mathrm{\ΦG}}_k(y_k-H_k\mathrm{\Φ}{\hat{x}}_{k-1})$

其中，

为由k时刻预测的k+1时刻的状态，Ф为状态转移矩阵， $\mathrm{\Φ}=\left(\begin{array}{ccc}1& T& {0.5T}^2\\ 0& 1& T\\ 0& 0& 1\end{array}\right),$ T＝1/PRF为脉冲重复间隔，PRF为雷达脉冲重复频率，

为由k-1时刻预测的k时刻的状态，G_k为k时刻增益矩阵，ФG_k＝ФP_k/k-1H^T[HP_k/k-1H^T+R_k]^-1，P_k/k-1为k时刻预测误差的协方差矩阵，P_k/k-1＝ΦP_k-1Ф^T+Q_k-1，P_k-1为k-1时刻平滑估计的协方差矩阵，Ф^T为Ф的转置，Q_k-1为k-1时刻动态噪声的协方差矩阵， $Q_{k-1}=\left[\begin{array}{c}0.5T\\ T\\ 1\end{array}\right]n(k-1),$ n(k-1)是均值为零的高斯白噪声序列在k-1时刻的取值，H为观测矩阵，H＝[1 0 0]，H^T为H的转置，R_k为k时刻观测噪声的协方差矩阵，(·)^-1为矩阵求逆运算，y_k为k时刻的观测值，为k-1时刻状态的估计值。

4h)在第k+1个方位单元的HRRP中，分别以各散射点在航迹矩阵中第k列记录的位置坐标为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+1时刻的位置幅度集合中；

4i)在第k+2个方位单元的HRRP中，分别以各散射点在第k+1时刻的位置幅度集合中所有位置坐标的均值为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+2时刻的位置幅度集合中；

4j)在第k+3个方位单元的HRRP中，分别以各散射点在第k+2时刻的位置幅度集合中所有位置坐标的均值为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+3时刻的位置幅度集合中；

4k)根据归一化的最近邻准则函数对散射点航迹进行关联，求散射点在第k+1时刻的真实位置坐标和幅度，并记录在航迹矩阵的第k+1列中。其中，归一化的最近邻准则函数如下：

$r_{k+1}=\arg \min \sqrt{{\left(\frac{A_{k+1}^{\′}-\stackrel{\‾}{A_{k+1}}}{\stackrel{\‾}{A_{k+1}}}\right)}^2+{\left(\frac{r_{k+1}^{\′}-{\hat{x}}_{k+1/k}\left(1\right)}{{\hat{x}}_{k+1/k}\left(1\right)}\right)}^{\text{2}}+{\left(\frac{r_{k+2}^{\′}-{\hat{x}}_{k+2/k+1}\left(1\right)}{{\hat{x}}_{k+2/k+1}\left(1\right)}\right)}^2+{\left(\frac{r_{k+3}^{\′}-{\hat{x}}_{k+3/k+2}\left(1\right)}{{\hat{x}}_{k+3/k+2}\left(1\right)}\right)}^2}$

其中，r_k为当前散射点在第k时刻的最优位置坐标，argmin(·)为求使目标函数取最小值时对应变量值的运算，A′_k+1为第k+1时刻的位置幅度集合中当前散射点对应的所有幅度值，为当前散射点在航迹矩阵中前k列的幅度均值，r′_k+1为第k+1时刻的位置幅度集合中当前散射点对应的所有位置坐标，

为由k时刻预测k+1时刻的状态向量中的位置坐标，r′_k+2为第k+2时刻的位置幅度集合中当前散射点对应的所有位置坐标，

为由k+1时刻预测k+2时刻的状态向量中的位置坐标，r′_k+3为第k+3时刻的位置集合中当前散射点对应的所有位置坐标，

为由k+2时刻预测k+3时刻的状态向量中的位置坐标。

41)将步骤4g)预测得到的k+1时刻状态向量中的速度和加速度按照步骤4f)进行计算，得到距离步进量，将k递增一个序号，重复步骤4g)至步骤4k)，直至完成散射点在所有方位单元的HRRP中的位置关联，从航迹矩阵中抽取散射点的位置坐标并且剔除对应幅度值，得到更新的航迹矩阵。

步骤5，高分辨成像

5a)将航迹矩阵中每个元素都减去该元素所在行所有元素的均值，得到平动校正后的航迹矩阵；

5b)对平动校正后的航迹矩阵进行带约束的矩阵分解得到散射点分布矩阵和雷达视线矩阵，其中带约束的矩阵分解步骤如下：

第一步，对平动校正后的航迹矩阵进行矩阵分解：

$W=\hat{R}\hat{S}$

其中，W为平动校正后的航迹矩阵，

为对W矩阵分解后得到的两个酉矩阵；

第二步，对下列方程求最小二乘解，得到尺度变换矩阵：

$l_n{\mathrm{AA}}^T{l}_n^T=I$

其中，l_n为

的第n个行向量，为第一步奇异值分解得到的酉矩阵，n∈[1，N_a]，N_a为方位向采样点数，A为尺度变换矩阵，A^T为A的转置，为l_n的转置，I为单位矩阵；

第三步，对下列非线性最小二乘方程组求解，得到旋转矩阵：

$\left\{\begin{array}{c}\hat{l}{A}_1=l\\ A_1{A}_1^T=I\end{array}\right.$

其中，l为初始时刻雷达视线矢量，

为从

中得到的对初始时刻雷达视线矢量的估计，

为第一步奇异值分解得到的酉矩阵，A₁为旋转矩阵，

为A₁的转置，I为单位矩阵；

第四步，获取真实的等效雷达视线旋转矩阵与散射点分布矩阵：

$R=\hat{R}{\mathrm{AA}}_1,$ $S=A_1^T{A}^{-1}\hat{S}$

其中，R为等效雷达视线旋转矩阵，S为散射点分布矩阵，

为第一步奇异值分解得到的两个酉矩阵，A为第二步得到的尺度变换矩阵，A₁为第三步得到的旋转矩阵。

5c)用散射点分布矩阵重构目标图像。

下面结合附图2对本发明的效果做进一步说明。

附图2所示的仿真图在MATLAB7.0软件下进行的，仿真数据的参数如下：雷达带宽为8GHz，载频为96GHz，PRF为150Hz，观测时间为0.1s，方位向回波为15次，接收回波的信噪比为10dB，仿真数据包含5个强散射点，目标旋转频率在[3.5π，4.5π]rad/s间随机变化，速度为4000m/s。按照仿真参数，方位无模糊采样所需的PRF应为9651Hz，因此仿真的PRF远远小于传统成像方法所需的PRF。

图2(a)是采用基于HRRP序列的空间目标高分辨成像方法对散射点航迹进行关联的结果，其中水平坐标表示回波序号，垂直坐标表示散射点航迹，分别用加号、乘号、星号、圆圈及实点代表了不同的航迹，从图中可看出，共五条散射点航迹。

图2(b)是目标旋转频率随回波的变化示意图，其中水平坐标表示回波序号，垂直坐标表示旋转角频率，单位为rad/s，由图中可看出旋转角频率随时间呈随机变化。

图2(c)是采用基于HRRP序列的空间目标高分辨成像方法对图2(a)进行高分辨成像的结果图，其中，水平坐标表示散射点的水平位置，垂直坐标表示散射点的垂直位置，单位均为米，实际的散射点位置用星号表示，估计出的散射点位置用圆圈表示，由图中可看出估计的所有散射点位置与实际散射点位置均重合到一起，说明虽然目标存在平动，旋转角频率随机变化并且脉冲重复频率较低，但本方法仍能够准确估计散射点位置，从而可以获得高分辨的图像，证明了该方法的有效性。

Claims (7)

1.基于HRRP序列的空间目标高分辨成像方法，包括如下步骤：

(1)雷达录取ISAR回波；

(2)高速运动补偿

2a)将参考距离为雷达到场景中心距离，频率、调频率与雷达发射信号相同的线性调频信号作为参考信号；

2b)将参考信号取共轭后与雷达录取的ISAR回波相乘得到解线调频处理结果；

2c)按照下式构造相位补偿函数：

其中，

为相位补偿函数，exp(·)为指数函数运算符，j为虚数单位，c为光速，v为测定的目标速度，t为方位维的采样时间，R(t)为t时刻场景中心到雷达的参考距离，f₁为信号载频，γ为调频率，f₂为距离频率；

2d)将相位补偿函数与步骤2b)得到的解线调频处理结果相乘得到高速运动补偿后回波；

(3)获取高质量HRRP

3a)按照下式构造过冗余字典：

$\mathrm{\Φ}=[\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(1\right)\right),\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(2\right)\right),\·\·\·,\exp \left(j\frac{4\mathrm{\π}}{c}{f}_2{R}_{\mathrm{\Δ}}\left(N\right)\right)]$

其中，Ф为包含N个基向量的过冗余字典，N取为距离频率点数的整数倍，exp(·)为指数函数运算符，j为虚数单位，c为光速，f₂为距离频率，

f₁为载频，γ为调频率，为距离维采样时间，R为参考距离，R_Δ(1)，R_Δ(2)，…R_Δ(N)为构造的N个散射点斜距，

W为雷达观测场景长度，相邻两个距离采样点的距离步进量满足关系

M为距离频率点数，B为信号带宽；

3b)令m＝1，其中，m表示方位单元序号；

3c)将步骤2d)得到的高速运动补偿后回波的第m个方位单元的复回波作为残余信号，并生成一个空矩阵Ф₁；

3d)用残余信号向步骤3a)得到的过冗余字典中的基向量作投影，并且计算内积；

3e)将内积最大的基向量记录在步骤3c)生成的矩阵Ф₁的第m列中，按照投影方法计算内积最大的基向量对应的投影系数，记录得到的投影系数，同时将内积最大的基向量从过冗余字典中剔除；

3f)从残余信号中减去矩阵Ф₁与记录的投影系数的乘积，得到更新的残余信号；

3g)重复步骤3d)至步骤3f)的操作，直至更新的残余信号能量低于噪声门限，将步骤3e)中记录的投影系数取模值作为对应基向量所代表斜距处散射点的幅度，以基向量对应的散射点斜距为横坐标，幅度为纵坐标绘制HRRP；

3h)将m递增1个序号，重复步骤3c)至步骤3g)的操作，直至绘制完所有方位单元对应的HRRP；

(4)散射点航迹关联

4a)令k＝1，其中，k表示方位单元序号，生成一个大小为N₁×M的空的航迹矩阵，其中，N₁为第k个方位单元对应的HRRP中峰值点的个数，M表示散射点回波次数，将第k个方位单元对应的HRRP中所有峰值点对应的位置坐标和幅度记录在航迹矩阵第k列中，将k递增一个序号；

4b)在第k个方位单元对应的HRRP中寻找峰值点，并将所有峰值点对应的位置坐标、幅度记录在第k时刻的位置幅度集合中，根据小转角内散射点幅度近似不变且位置接近准则，将第k时刻位置幅度集合中的元素与航迹矩阵第k-1列中元素进行关联，将关联得到的第k时刻散射点位置坐标和幅度记录在航迹矩阵第k列；

4c)将k递增一个序号，执行步骤4b)；

4d)分别将各散射点在航迹矩阵中第k-1列记录的位置坐标与第k-2列记录的位置坐标相减，并除以脉冲重复周期得到速度，将航迹矩阵中第k列记录的位置坐标减去第k-1列记录的位置坐标的两倍，再加上第k-2列记录的位置坐标，将计算结果除以脉冲重复周期的平方得到加速度；

4e)用速度乘以脉冲重复周期并加上加速度与脉冲重复周期平方乘积的0.5倍，求得距离步进量；

4f)分别将各散射点在航迹矩阵中第k列记录的位置坐标、步骤4d)得到的速度与加速度作为第k时刻散射点状态向量；

4g)根据Kalman预测准则，对第k+1、k+2及k+3时刻的状态向量进行预测，并将航迹矩阵中记录的对应散射点幅度求均值作为预测幅度值；

4h)在第k+1个方位单元对应的HRRP中，分别以各散射点在航迹矩阵中第k列记录的位置坐标为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+1时刻的位置幅度集合中；

4i)在第k+2个方位单元对应的HRRP中，分别以各散射点在第k+1时刻的位置幅度集合中所有位置坐标的均值为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+2时刻的位置幅度集合中；

4j)在第k+3个方位单元对应的HRRP中，分别以各散射点在第k+2时刻的位置幅度集合中所有位置坐标的均值为中心，以距离步进量的三倍为搜索宽度寻找峰值点，将所有峰值点对应的位置坐标和幅度记录在第k+3时刻的位置幅度集合中；

4k)根据归一化的最近邻准则函数对散射点航迹进行关联，求得散射点在第k+1时刻的真实位置坐标和幅度，并记录在航迹矩阵第k+1列中；

41)将步骤4g)预测得到的k+1时刻状态向量中的速度和加速度按照步骤4e)进行计算，得到距离步进量，将k递增一个序号，重复步骤4g)至步骤4k)，直至完成散射点在所有方位单元对应的HRRP中的位置关联，从航迹矩阵中抽取散射点的位置坐标并且剔除对应幅度值，得到更新的航迹矩阵；

(5)高分辨成像

5a)将航迹矩阵中每个元素都减去该元素所在行所有元素的均值，得到平动校正后的航迹矩阵；

5b)对平动校正后的航迹矩阵进行带约束的矩阵分解得到散射点分布矩阵和雷达视线矩阵；

5c)用散射点分布矩阵重构目标图像。

2.根据权利要求1所述的基于HRRP序列的空间目标高分辨成像方法，其特征在于，步骤3e)中所述的投影方法如下：

$\mathrm{\α}={\left({\mathrm{\Φ}}_1^H{\mathrm{\Φ}}_1\right)}^{-1}{\mathrm{\Φ}}_1^{H}g$

其中，α为投影系数，Ф₁为记录最大投影系数对应基向量的矩阵，为Ф₁的共轭转置，(·)^-1为矩阵求逆运算，g为残余信号。

3.根据权利要求1所述的基于HRRP序列的空间目标高分辨成像方法，其特征在于，步骤3g)中所述的噪声门限为，对步骤3c)得到的残余信号进行逆傅里叶变换并取模值，将取模值后的图像中远离峰值点的其余点对应的平均功率作为噪声门限。

4.根据权利要求1所述的基于HRRP序列的空间目标高分辨成像方法，其特征在于，步骤4b)中所述的小转角内散射点幅度近似不变且位置接近准则为，当前时刻散射点幅度应该满足如下等式：

A_k＝argmin||A′_k-A_k-1||₂

其中，A_k为当前散射点在第k时刻对应的最优幅度值，argmin(·)为求使目标函数取最小值时对应变量值的运算，||·||₂求为2范数，A′_k为第k时刻的位置幅度集合中当前散射点对应的所有幅度值，A_k-1为航迹矩阵中第k-1列中当前散射点对应的幅度值。

5.根据权利要求1所述的基于HRRP序列的空间目标高分辨成像方法，其特征在于，步骤4g)中所述的Kalman预测准则如下：

${\hat{x}}_{k+1/k}=\mathrm{\Φ}{\hat{x}}_{k/k-1}+{\mathrm{\ΦG}}_k(y_k-H_k\mathrm{\Φ}{\hat{x}}_{k-1})$

其中，

为由k时刻预测的k+1时刻的状态，Ф为状态转移矩阵， $\mathrm{\Φ}=\left(\begin{array}{ccc}1& T& {0.5T}^2\\ 0& 1& T\\ 0& 0& 1\end{array}\right),$ T＝1/PRF为脉冲重复间隔，PRF为雷达脉冲重复频率，

为由k-1时刻预测的k时刻的状态，G_k为k时刻增益矩阵，ФG_k＝ФP_k/k-1H^T[HP_k/k-1H^T+R_k]^-1，P_k/k-1为k时刻预测误差的协方差矩阵，P_k/k-1＝ΦP_k-1Ф^T+Q_k-1，P_k-1为k-1时刻平滑估计的协方差矩阵，Ф^T为Ф的转置，Q_k-1为k-1时刻动态噪声的协方差矩阵， $Q_{k-1}=\left[\begin{array}{c}0.5T\\ T\\ 1\end{array}\right]n(k-1),$ n(k-1)是均值为零的高斯白噪声序列在k-1时刻的取值，H为观测矩阵，H＝[1 0 0]，H^T为H的转置，R_k为k时刻观测噪声的协方差矩阵，(·)^-1为矩阵求逆运算，y_k为k时刻的观测值，

为k-1时刻状态的估计值。

6.根据权利要求1所述的基于HRRP序列的空间目标高分辨成像方法，其特征在于，步骤4k)中所述的归一化的最近邻准则函数如下：

$r_{k+1}=\arg \min \sqrt{{\left(\frac{A_{k+1}^{\′}-\stackrel{\‾}{A_{k+1}}}{\stackrel{\‾}{A_{k+1}}}\right)}^2+{\left(\frac{r_{k+1}^{\′}-{\hat{x}}_{k+1/k}\left(1\right)}{{\hat{x}}_{k+1/k}\left(1\right)}\right)}^{\text{2}}+{\left(\frac{r_{k+2}^{\′}-{\hat{x}}_{k+2/k+1}\left(1\right)}{{\hat{x}}_{k+2/k+1}\left(1\right)}\right)}^2+{\left(\frac{r_{k+3}^{\′}-{\hat{x}}_{k+3/k+2}\left(1\right)}{{\hat{x}}_{k+3/k+2}\left(1\right)}\right)}^2}$

其中，r_k为当前散射点在第k时刻的最优位置坐标，argmin(·)为求使目标函数取最小值时对应变量值的运算，A′_k+1为第k+1时刻的位置幅度集合中当前散射点对应的所有幅度值，

为当前散射点在航迹矩阵中前k列的幅度均值，r′_k+1为第k+1时刻的位置幅度集合中当前散射点对应的所有位置坐标，

为由k时刻预测k+1时刻的状态向量中的位置坐标，r′_k+2为第k+2时刻的位置幅度集合中当前散射点对应的所有位置坐标，为由k+1时刻预测k+2时刻的状态向量中的位置坐标，r′_k+3为第k+3时刻的位置集合中当前散射点对应的所有位置坐标，

为由k+2时刻预测k+3时刻的状态向量中的位置坐标。

7.根据权利要求1所述的基于HRRP序列的空间目标高分辨成像方法，其特征在于，步骤5b)中所述的带约束的矩阵分解步骤如下：

第一步，对平动校正后的航迹矩阵进行分解：

$W=\hat{R}\hat{S}$

其中，W为平动校正后的航迹矩阵，

为对W矩阵分解后得到的两个酉矩阵；

第二步，对下列方程求最小二乘解，得到尺度变换矩阵：

$l_n{\mathrm{AA}}^T{l}_n^T=I$

其中，l_n为

的第n个行向量，

为第一步奇异值分解得到的酉矩阵，n∈[1，N_a]，N_a为方位向采样点数，A为尺度变换矩阵，A^T为A的转置，

为l_n的转置，I为单位矩阵；

第三步，对下列非线性最小二乘方程组求解，得到旋转矩阵：

$\left\{\begin{array}{c}\hat{l}{A}_1=l\\ A_1{A}_1^T=I\end{array}\right.$

其中，l为初始时刻雷达视线矢量，

为从

中得到的对初始时刻雷达视线矢量的估计，

为第一步奇异值分解得到的酉矩阵，A₁为旋转矩阵，

为A₁的转置，I为单位矩阵；

第四步，获取真实的等效雷达视线旋转矩阵与散射点分布矩阵：

$R=\hat{R}{\mathrm{AA}}_1,$ $S=A_1^T{A}^{-1}\hat{S}$

其中，R为等效雷达视线旋转矩阵，S为散射点分布矩阵，

为第一步奇异值分解得到的两个酉矩阵，A为第二步得到的尺度变换矩阵，A₁为第三步得到的旋转矩阵。

Images