CN109613532B - 一种机载雷达实时多普勒波束锐化超分辨成像方法 - Google Patents

Abstract

本发明提出了一种机载雷达实时多普勒波束锐化超分辨成像方法，属于雷达探测与成像领域。本发明通过边扫描‑边获取‑边处理的方式，对当前相邻相干处理间隔(CPI)的回波信号，采用快速迭代加权最小二乘估计实现多普勒波束锐化(DBS)超分辨成像，再利用相邻CPI方位回波的空间相关性，对下一CPI加权最小二乘估计结果进行递归更新，从而提升收敛速度，极大降低计算复杂度。与传统DBS成像技术相比，本发明不仅可以有效提高成像分辨率，还可以减少迭代次数，提高收敛速度，保证成像分辨率的同时，满足了成像的实时性。

Description

一种机载雷达实时多普勒波束锐化超分辨成像方法

技术领域

本发明属于雷达探测与成像领域，特别涉及一种机载雷达实时多普勒波束锐化超分辨成像方法。

背景技术

多普勒波束锐化(DBS)成像具有成像视角宽，实时性好的特点，广泛应用于地面侦察、地形测绘、飞机导航等。传统DBS成像主要基于快速傅里叶变换(FFT)谱分析方法，分辨率低下，在慢速情况下，如直升机自主着陆等应用中无法满足成像分辨率要求，因此，基于超分辨谱估计算法的DBS成像技术被广泛研究。

为了提高DBS成像分辨率，现有技术中有一种DBS超分辨成像方法，该方法能够提高DBS成像分辨率，该方法需要矩阵乘法和矩阵求逆等运算，无法实现DBS实时超分辨成像；还有一种DBS超分辨成像方法，通过快速傅里叶变换和低位移秩表示，实现了矩阵快速相乘，避免了矩阵求逆，该方法虽然提高了DBS超分辨成像的运算效率，但其过多的迭代次数为实际工程应用带来运算负担。

发明内容

本发明的目的在于克服现有技术中存在计算复杂度高及冗余迭代等导致无法实现DBS实时超分辨成像的缺陷，提出了一种机载雷达实时多普勒波束锐化超分辨成像方法，通过边扫描-边获取-边处理的方式，对当前相邻相干处理间隔(CPI)的回波信号，采用快速迭代加权最小二乘估计实现DBS超分辨成像，再利用相邻CPI方位回波的空间相关性，对下一CPI加权最小二乘估计结果进行递归更新，从而提升收敛速度，极大降低计算复杂度。与传统DBS成像技术相比，本发明不仅可以有效提高成像分辨率，还可以实现边扫描-边成像实时信号处理。

一种机载雷达实时多普勒波束锐化超分辨成像方法，包括以下步骤：

S1、获取当前相干处理间隔(Coherent process interval，CPI)回波信号，对回波信号进行脉冲压缩和距离走动校正处理；

S2、判断当前回波信号是否为第一个CPI序号，若当前回波信号为第一个CPI回波信号，流程进入步骤S3；否则流程进入步骤S4；

S3、计算当前第一CPI回波信号的自相关矩阵，流程进入步骤S5；

S4、更新前一CPI回波信号迭代收敛后的自相关矩阵得到当前CPI回波信号的自相关矩阵，流程进入步骤S5；

S5、对当前CPI回波信号的自相关矩阵进行求逆运算，得到自相关矩阵的逆矩阵；

S6、投影当前CPI回波信号，得到该CPI回波信号的多普勒波束锐化(Doppler beamsharpening，DBS)成像结果。

进一步地，所述步骤S1包括以下流程：

载机平台的飞行高度为H，运动方向沿X轴方向，速度为V，初始位置的波束俯仰角和水平方位角分别为

和θ₀，空间方位角为α₀，有

雷达波束以速度w对成像场景进行扫描，载机平台初始时刻与目标之间的距离为R₀，则空间中一目标P(x,y)，在初始时刻，载机平台与目标的瞬时距离为

其中，t为方位向慢时间变量；

雷达天线发射线性调频(linear frequency modulation，LFM)信号，在回波信号的频域乘以脉冲压缩匹配函数

和距离走动校正因子

再变换到时域，其中，f为距离向频率变量，K_r为发射信号调频斜率，α为各方位目标空间方位角度，c为电磁波传播速度，经过脉冲压缩和距离走动校正后的回波信号的时域表达式为

其中，τ为距离向快时间变量，σ为目标散射系数，A(θ)为天线方向图图函数在方位向的调制，B为发射线性调频信号的带宽，f₀为载波频率。

进一步地，所述步骤S3包括以下流程：

当前CPI回波信号为第一个CPI回波信号时，将经过处理的回波信号表示为离散形式，s(m,n)表示第m个距离单元和第n个方位角单元的回波，m＝1,2,...,M，n＝1,2,...,N，第m个距离单元的回波向量S＝[s(1,m),s(2,m),...,s(N,m)]^T为

其中，a_k＝[e^j2πk/K,e^j2π2k/K,...,e^j2πNk/K]^T，k＝0,1...,K-1，[a₁,a₂,...,a_k,...,a_K]为多普勒相位矩阵，K为多普勒频域采样点数，T为转置运算，[σ(1,m),σ(2,m),...,σ(K,m)]^T为第m个距离单元内每一个多普勒频率点的目标散射系数，e＝[e(1,m),e(2,m),...,e(N,m)]^T为系统加性高斯白噪声；

根据加权最小二乘准则，自相关矩阵为

其中，

表示初始迭代时，第一个CPI的自相关矩阵；

n＝0,1,...,N-1，*表示共轭运算，且

n的其值范围为-N+1,...,-1,0,1,...,N-1，通过k点快速傅里叶变换(FFT)得到p的初始值，R的首行向量通过IFFT得到，根据自相关矩阵R具有Hermitian特征和Toeplitz矩阵结构，可通过一行元素的计算，获得自相关矩阵的计算结果；

根据协方差矩阵和自相关矩阵的关系Q_k＝R-p_ka_k(a_k)^H，及矩阵求逆原理，得到σ_k的最优估计

其中，

且P＝diag(p)，p_k＝||σ_k||²表示在频率网格点2πk/K处的功率估计值，每一次迭代更新p_k。

进一步地，所述步骤S4包括以下流程：

根据所述步骤S3，当迭代加权最小二乘估计到收敛时，相邻迭代的自相关矩阵有

其中，

表示第i次迭代后，第n-1个CPI的自相关矩阵，i＝1,...,ξ，n＝2,...,N；

根据相邻CPI回波信号的空间相关性，从第二个CPI起，对回波信号处理序列中散射系数的初始化可迭代更新为

其中，

表示第i次迭代后，第n-1个CPI的最终散射系数，n＝2,...,N；

第n个CPI的自相关矩阵的初值可迭代更新为

其中，

为第n个CPI的自相关矩阵的迭代初始值。

进一步地，所述步骤S5包括以下流程：

通过对角加载对自相关矩阵进行处理从而保持矩阵满秩，更新后自相关矩阵表示为

其中，I为单位矩阵，λ为正则化参数；

采用Gohberg-Semencul分解算法，更新后的自相关矩阵的逆矩阵表示为

其中，L(u,Z)＝(u,Du,D²u,...,D^N-1u)，

D为转移矩阵，且有

根据自相关矩阵特性，利用Levinson-Durbin算法求得向量u和

根据式

直接计算自相关矩阵的逆矩阵R'^-1；根据

分子为对R'^-1S做FFT处理，分母为

通过Toeplitz矩阵向量运算可获取ρ的值，对ρ做IFFT变换，可获得分母的值；最终计算得到第一次迭代

的值；

相邻迭代归一化相对误差为

其中，

为L₂范数，归一化相对误差门限设置为Δσ，当相邻迭代结果的归一化误差小于误差门限值，即迭代收敛，获得该次CPI的

结果。

进一步地，所述步骤S6包括以下流程：

根据设定的归一化相对误差门限Δσ判断了迭代是否收敛，得到

在回波信号距离-多普勒域中的运算结果，沿多普勒频率方向，根据目标成像区域的回波多普勒分布范围，投影该CPI回波成像结果至回波距离-方位域，获得该CPI回波信号的DBS成像结果。

本发明的有益效果：本发明提供了一种机载雷达实时多普勒波束锐化超分辨成像方法，根据机载扫描雷达回波录取规律，基于DBS成像中相邻CPI回波之间的空间相关性，用前一个CPI回波的自相关矩阵更新当前CPI回波自相关矩阵的迭代初始值，此时可避免自相关矩阵的运算，接着采用Levinson-Durbin算法、Toeplitz矩阵与向量快速相乘完成自相关矩阵快速求逆，实现目标散射系数的快速估计。本发明可以减少迭代次数，提高收敛速度，在保证成像分辨率的同时，满足了成像的实时性。

附图说明

图1为本发明实施例提供的流程图。

图2为本发明实施例提供的机载扫描雷达几何模型示意图。

图3为本发明实施例提供的原始成像场景图。

图4为机载雷达实波束扫描成像结果图。

图5为现有技术中超分辨成像结果图。

图6为本发明实施例提供的方法的成像结果图。

图7为不同成像方法沿序列CPI回波收敛次数对比图。

具体实施方式

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

请参阅图1，本发明提供了一种机载雷达实时多普勒波束锐化超分辨成像方法，通过以下步骤实现：

S1、获取当前相干处理间隔(Coherent process interval，CPI)回波信号，对回波信号进行脉冲压缩和距离走动校正处理。

本实施例中，机载扫描雷达几何模型示意图如图2所示，雷达平台参数如下表(一)所示，目标原始场景图如图3所示，设置的菱形阵列点目标，证明本发明方法的有效性。

参数	符号	数值
			载频	f<sub>0</sub>	30.75GHz
发射信号时宽	T<sub>p</sub>	2μs
			发射信号带宽	B	40MHz
平台高度	H	1km
			脉冲重复频率	PRF	4000Hz
天线扫描速度	w	30°/s
			天线波束宽度	θ<sub>β</sub>	4°
扫描范围	Φ	6°～14°
			机载平台速度	V	100m/s
距离范围	R<sub>0</sub>	19900m～20100m
			归一化相对误差门限	Δσ	10<sup>-4</sup>

表(一)雷达平台参数表

本实施例中，载机平台的飞行高度为H，运动方向沿X轴方向，速度为V，初始位置的波束俯仰角和水平方位角分别为

和θ₀，空间方位角为α₀，由几何关系可知：

雷达波束以速度w对成像场景进行扫描，载机平台初始时刻与目标之间的距离为R₀，则空间中一目标P(x,y)，在初始时刻，载机平台与目标的瞬时距离为

其中，t为方位向慢时间变量；

雷达天线发射大时宽带宽积的线性调频(linear frequency modulation，LFM)信号，在回波信号的频域乘以脉冲压缩匹配函数

和距离走动校正因子

再变换到时域，其中，f为距离向频率变量，K_r为发射信号调频斜率，α为各方位目标空间方位角度，c为电磁波传播速度。为了提高处理效率，脉冲压缩和距离走动校正可同时进行。

根据雷达参数仿真获取回波数据，当天线扫过斜前视区域Ω时，根据回波数据的序列特性，对于回波数据采取边扫描边处理的成像工作模式。当获取到第一个CPI回波数据，经过脉冲压缩和距离走动校正后的回波信号的时域表达式为

其中，τ为距离向快时间变量，σ为目标散射系数，A(θ)为天线方向图图函数在方位向的调制，B为发射线性调频信号的带宽，f₀为载波频率。

通过脉冲压缩和距离走动校正预处理，实现距离向的高分辨率。

S2、判断当前回波信号是否为第一个CPI序号，若当前回波信号为第一个CPI回波信号，流程进入步骤S3；否则流程进入步骤S4。

本实施例中，通过步骤S1，可依次获得不同的CPI回波信号。判断当前处理的CPI回波信号是否为第一个CPI回波信号，若是，则流程进入步骤S3；否则流程进入步骤S4.

S3、计算当前第一CPI回波信号的自相关矩阵，流程进入步骤S5。

本实施例中，当前CPI回波信号为第一个CPI回波信号时，将经过处理的距离向高分辨二维回波信号表示为离散形式，s(m,n)表示第m个距离单元和第n个方位角单元的回波，m＝1,2,...,M，n＝1,2,...,N，第m个距离单元的回波向量S＝[s(1,m),s(2,m),...,s(N,m)]^T为

其中，a_k＝[e^j2πk/K,e^j2π2k/K,...,e^j2πNk/K]^T，k＝0,1...,K-1，[a₁,a₂,...,a_k,...,a_K]为多普勒相位矩阵，K为多普勒频域采样点数，T为转置运算，[σ(1,m),σ(2,m),...,σ(K,m)]^T为第m个距离单元内每一个多普勒频率点的目标散射系数，e＝[e(1,m),e(2,m),...,e(N,m)]^T为系统加性高斯白噪声；

根据加权最小二乘准则，自相关矩阵为

其中，

表示初始迭代时，第一个CPI的自相关矩阵；

n＝0,1,...,N-1，*表示共轭运算，且

n的其值范围为-N+1,...,-1,0,1,...,N-1，通过k点快速傅里叶变换(FFT)得到p的初始值，R的首行向量通过IFFT得到，根据自相关矩阵R具有Hermitian特征和Toeplitz矩阵结构，可通过一行元素的计算，获得自相关矩阵的计算结果；

根据协方差矩阵和自相关矩阵的关系Q_k＝R-p_ka_k(a_k)^H，及矩阵求逆原理，得到σ_k的最优估计

其中，

且P＝diag(p)，p_k＝||σ_k||²表示在频率网格点2πk/K处的功率估计值，每一次迭代更新p_k。

S4、更新前一CPI回波信号迭代收敛后的自相关矩阵得到当前CPI回波信号的自相关矩阵，流程进入步骤S5。

本实施例中，根据步骤S3，当迭代加权最小二乘估计到收敛时，相邻迭代的自相关矩阵有如下关系

其中，

表示第i次迭代后，第n-1个CPI的自相关矩阵，i＝1,...,ξ，n＝2,...,N；

根据相邻CPI回波信号的空间相关性，从第二个CPI起，对回波信号处理序列中散射系数的初始化可迭代更新为

其中，

表示第i次迭代后，第n-1个CPI的最终散射系数，

n＝2,...,N；

第n个CPI的自相关矩阵的初值可迭代更新为

其中，

为第n个CPI的自相关矩阵的迭代初始值。

根据该步骤，利用公式(8)，当获得

可直接进行下一CPI回波处理。

S5、对当前CPI回波信号的自相关矩阵进行求逆运算，得到自相关矩阵的逆矩阵。

本实施例中，通过上述步骤，可以得到当前CPI回波信号的起始迭代自相关矩阵。同时，可通过对角加载对自相关矩阵进行处理从而保持矩阵满秩，更新后自相关矩阵表示为

其中，I为单位矩阵，λ为正则化参数；

采用Gohberg-Semencul分解算法，更新后的自相关矩阵的逆矩阵表示为

其中，L(u,Z)＝(u,Du,D²u,...,D^N-1u)，

D为转移矩阵，且有

根据自相关矩阵特性，利用Levinson-Durbin算法求得向量u和

根据式

直接计算自相关矩阵的逆矩阵R'^-1；根据

分子为对R'^-1S做FFT处理，分母为

通过Toeplitz矩阵向量运算可获取ρ的值，对ρ做IFFT变换，可获得分母的值；最终计算得到第一次迭代

的值；

定义相邻迭代归一化相对误差为

其中，

为L₂范数，归一化相对误差门限设置为Δσ，当相邻迭代结果的归一化误差小于误差门限值，即迭代收敛，获得该次CPI的

结果。

S6、投影当前CPI回波信号，得到该CPI回波信号的多普勒波束锐化(Doppler beamsharpening，DBS)成像结果。

本实施例中，利用上述步骤，根据设定的归一化相对误差门限Δσ判断是否收敛，获得了

在回波信号距离-多普勒域中的运算结果，沿多普勒频率方向，根据目标成像区域的回波多普勒分布范围，投影该CPI回波成像结果至回波距离-方位域，获得该CPI回波信号的DBS成像结果。雷达系统持续扫描成像场景，自第二个CPI回波数据起，将由前一个CPI回波数据得到自相关矩阵的迭代收敛值作为当前CPI回波数据的迭代初值，即：

得到自相关矩阵的迭代更新值，前n个CPI回波数据按照其实际距离-方位域坐标进行成像结果投影，直至雷达系统停止工作。

从图4的实波束成像结果、图5的现有技术中超分辨成像结果和图6的本发明方法成像结果的对比可知，本发明方法实现了方位的高分辨，在提升成像效率的同时，保证成像性能与传统FIAA方法一致，均实现了场景中目标的有效分辨。

请参阅图7，不同的成像方法沿序列CPI回波收敛迭代次数对比图，FIAA方法的平均迭代次数为20余次，本发明方法的平均迭代次数约为6次。验证了本发明方法利用相邻波束间的空间相关性，能够进一步减少了序列CPI回波数据成像处理的迭代次数，有效的提高迭代收敛速度，与FIAA实现DBS成像相比，本发明方法可将迭代收敛次数降低约4倍，进一步地提高了成像的计算效率。

本领域的普通技术人员将会意识到，这里的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (3)

1.一种机载雷达实时多普勒波束锐化超分辨成像方法，其特征在于，包括以下步骤：

S1、获取当前相干处理间隔(Coherent process interval，CPI)回波信号，对回波信号进行脉冲压缩和距离走动校正处理；包括以下流程：

载机平台的飞行高度为H，运动方向沿X轴方向，速度为V，初始位置的波束俯仰角和水平方位角分别为

和θ₀，空间方位角为α₀，有

雷达波束以速度w对成像场景进行扫描，载机平台初始时刻与目标之间的距离为R₀，则空间中一目标P(x,y)，在初始时刻，载机平台与目标的瞬时距离为

其中，t为方位向慢时间变量；

雷达天线发射线性调频(linear frequency modulation，LFM)信号，在回波信号的频域乘以脉冲压缩匹配函数

和距离走动校正因子

再变换到时域，其中，f为距离向频率变量，K_r为发射信号调频斜率，α为各方位目标空间方位角度，c为电磁波传播速度，经过脉冲压缩和距离走动校正后的回波信号的时域表达式为

其中，τ为距离向快时间变量，σ为目标散射系数，A(θ)为天线方向图图函数在方位向的调制，B为发射线性调频信号的带宽，f₀为载波频率；

S2、判断当前回波信号是否为第一个CPI序号，若当前回波信号为第一个CPI回波信号，流程进入步骤S3；否则流程进入步骤S4；

S3、计算当前第一CPI回波信号的自相关矩阵，流程进入步骤S5；包括以下流程：

当前CPI回波信号为第一个CPI回波信号时，将经过处理的回波信号表示为离散形式，s(m,n)表示第m个距离单元和第n个方位角单元的回波，m＝1,2,...,M，n＝1,2,...,N，第m个距离单元的回波向量S＝[s(1,m),s(2,m),...,s(N,m)]^T为

其中，a_k＝[e^j2πk/K,e^j2π2k/K,...,e^j2πNk/K]^T，k＝0,1...,K-1，[a₁,a₂,...,a_k,...,a_K]为多普勒相位矩阵，K为多普勒频域采样点数，T为转置运算，[σ(1,m),σ(2,m),...,σ(K,m)]^T为第m个距离单元内每一个多普勒频率点的目标散射系数，e＝[e(1,m),e(2,m),...,e(N,m)]^T为系统加性高斯白噪声；

根据加权最小二乘准则，自相关矩阵为

其中，

表示初始迭代时，第一个CPI的自相关矩阵；

n＝0,1,...,N-1，^*表示共轭运算，且

n的其值范围为-N+1,...,-1,0,1,...,N-1，通过k点快速傅里叶变换(FFT)得到p的初始值，R的首行向量通过IFFT得到，根据自相关矩阵R具有Hermitian特征和Toeplitz矩阵结构，可通过一行元素的计算，获得自相关矩阵的计算结果；

根据协方差矩阵和自相关矩阵的关系Q_k＝R-p_ka_k(a_k)^H，及矩阵求逆原理，得到σ_k的最优估计

其中，

且P＝diag(p)，p_k＝||σ_k||²表示在频率网格点2πk/K处的功率估计值，每一次迭代更新p_k；

S4、更新前一CPI回波信号迭代收敛后的自相关矩阵得到当前CPI回波信号的自相关矩阵，流程进入步骤S5；包括以下流程：

根据所述步骤S3，当迭代加权最小二乘估计到收敛时，相邻迭代的自相关矩阵有

其中，

表示第i次迭代后，第n-1个CPI的自相关矩阵，i＝1,...,ξ，n＝2,...,N；

根据相邻CPI回波信号的空间相关性，从第二个CPI起，对回波信号处理序列中散射系数的初始化可迭代更新为

其中，

表示第i次迭代后，第n-1个CPI的最终散射系数，

n＝2,...,N；

第n个CPI的自相关矩阵的初值迭代更新为

其中，

为第n个CPI的自相关矩阵的迭代初始值；

S5、对当前CPI回波信号的自相关矩阵进行求逆运算，得到自相关矩阵的逆矩阵；

S6、投影当前CPI回波信号，得到该CPI回波信号的多普勒波束锐化(Doppler beamsharpening，DBS)成像结果。

2.如权利要求1所述的机载雷达实时多普勒波束锐化超分辨成像方法，其特征在于，所述步骤S5包括以下流程：

通过对角加载对自相关矩阵进行处理从而保持矩阵满秩，更新后自相关矩阵表示为

其中，I为单位矩阵，λ为正则化参数；

采用Gohberg-Semencul分解算法，更新后的自相关矩阵的逆矩阵表示为

其中，L(u,Z)＝(u,Du,D²u,...,D^N-1u)，

D为转移矩阵，且有

根据自相关矩阵特性，利用Levinson-Durbin算法求得向量u和

根据式

直接计算自相关矩阵的逆矩阵R'^-1；根据

分子为对R'^-1S做FFT处理，分母为

通过Toeplitz矩阵向量运算可获取ρ的值，对ρ做IFFT变换，可获得分母的值；最终计算得到第一次迭代

的值；

相邻迭代归一化相对误差为

其中，

为L₂范数，归一化相对误差门限设置为Δσ，当相邻迭代结果的归一化误差小于误差门限值，即迭代收敛，获得该次CPI的

结果。

3.如权利要求2所述的机载雷达实时多普勒波束锐化超分辨成像方法，其特征在于，所述步骤S6包括以下流程：

根据设定的归一化相对误差门限Δσ判断了迭代是否收敛，得到

在回波信号距离-多普勒域中的运算结果，沿多普勒频率方向，根据目标成像区域的回波多普勒分布范围，投影该CPI回波成像结果至回波距离-方位域，获得该CPI回波信号的DBS成像结果。

Images