CN105425234A - 基于多任务贝叶斯压缩感知的距离-多普勒成像方法 - Google Patents

Abstract

本发明公开一种基于多任务贝叶斯压缩感知的距离-多普勒成像方法，主要解决现有技术复杂度高分辨率低的问题。其技术方案是：1)雷达向目标场景发射脉冲调制序列，并同时得到回波信号；2)将目标场景距离范围和多普勒范围进行等间隔划分，对回波信号进行离散化采样，构造回波矩阵；3)引用辅助变量，构造优化求解函数；4)对优化函数进行求解，得到场景的距离-多普勒成像结果。本发明能直接对场景距离-多普勒维进行重构，在降低算法复杂度的同时提高了分辨率，而且不需要复杂的参数设置，可用于战场监视、导弹制导、空间探测等领域。

Description

基于多任务贝叶斯压缩感知的距离-多普勒成像方法

技术领域

本发明属于雷达信号处理技术领域，具体地说是一种距离-多普勒成像方法。可用于雷达目标位置和速度的估计。

背景技术

雷达成像的出现扩展了雷达的原始概念，使雷达具有了对区域目标和运动目标进行成像和识别能力，同时它在微波遥感应用方面表现出巨大的潜力，为人们提供越来越多的有用信息。雷达通常具有距离、多普勒以及角度三维分辨信息。随着雷达技术的不断发展和人们对高分辨雷达的迫切需求，对多维信息量的获取越发显得重要。

雷达成像的根本原理是采用各种方法提高雷达各维的分辨率。为了提高距离分辨率，发射脉冲序列，其中每个脉冲由调制的子脉冲组成。在此基础上，有三种方法用于距离-多普勒成像。

第一种失配滤波器法，这种方法由于用脉冲序列代替单脉冲，虽然能得到更高的距离分辨率，但却不能得到高的多普勒分辨率。

第二种迭代自适应方法，是一种数据自适应方法，这种方法虽然克服了失配滤波器在多普勒分辨率上的缺陷，能够得到很高的多普勒分辨率，但是这种方法在有噪声存在时成像不稳定，而且计算量大。

第三种方法是利用距离多普勒像的稀疏先验性，将距离-多普勒成像转换为稀疏信号恢复问题进行求解，求解速度快，分辨率高，但却由于其复杂的参数设置而影响了该方法的实际使用。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于多任务贝叶斯压缩感知的距离-多普勒成像方法，以在提高距离-多普勒频谱分辨率的同时，降低成像的复杂度。

为实现上述目的，本发明的技术方案包括如下步骤：

(1)雷达向目标场景H发射Q个脉冲调制序列并在发射脉冲的同时接收该场景H的回波脉冲，得到第q个脉冲的回波信号为y(q)；

(2)对回波信号y(q)进行离散化采样，构造回波矩阵：Y＝DXΘ+E，其中，D为距离维基矩阵，X为散射系数矩阵，Θ为多普勒域基矩阵，E为噪声矩阵；

(3)引入辅助变量矩阵Z，令Z＝XΘ，将回波矩阵变形为下式：

$\left\{\begin{array}{c}Y=DZ+E\\ Z^T={\mathrm{\Θ}}^T{X}^T\end{array}\right.;$ 其中(·)^T表示矩阵转置；

(4)将步骤(3)中第一个方程记作下式：

y_i＝Dz_i+e_i,i＝1,…,Q

其中，y_i为Y的第i列，z_i为Z的第i列，e_i为E的第i列，Q为回波矩阵Y的列数；

(5)对步骤(4)中的式子进行求解，得到辅助变量矩阵Z，设定检测门限Δ＝0.1，筛选出辅助变量矩阵Z中的非零行组成新矩阵：Z_select＝X_selectΘ，对该式进行转置，得到：

$Z_{select}^T={\mathrm{\Θ}}^T{X}_{select}^T$

其中，X_select为散射系数矩阵X剔除全零行后的矩阵；

(6)将步骤(5)中方程记作：z_i＝Θ^Tx_i,i＝1,…,R，其中，z_i为的第i列，x_i为的第i列；

(7)对步骤(6)中的式子进行求解，得到散射系数矩阵X剔除全零行后的矩阵再对补全零行并取转置，得到散射系数矩阵X，由散射系数矩阵X即可得到距离-多普勒像。

本发明与现有技术相比具有如下优点：

第一，本发明采用多任务贝叶斯压缩感知方法，无需对信号的噪声方差进行估计，且省略了调参的步骤。

第二，本发明采用二维场景图像进行重构成像，避免了传统稀疏成像建模中，场景需向量化处理带来的巨大运算量，且具有更高的分辨率。

附图说明

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

图2是本发明使用的距离-多普勒二维谱结果图；

图3是用现有的稀疏成像法得到的距离-多普勒二维谱的仿真结果图；

图4是用本发明方法得到的距离-多普勒二维谱的仿真结果图。

具体实施方法

下面结合附图对本发明做进一步的详细描述。

参照图1，本发明的实施步骤如下：

步骤1：获得回波信号y(q)。

雷达发射Q个脉冲，每个脉冲包含N个子脉冲，这N个子脉冲为Barker序列得到第q个脉冲的回波信号y(q)：

其中，x_r,l是距离为r多普勒频率为l的目标的散射强度，τ_r为延迟时间，为多普勒相移，q＝1,2,…,Q。

步骤2：对回波信号进行离散化采样，构造回波矩阵。

(2a)在距离-多普勒二维探测区域将探测距离范围r_min～r_max进行等间隔划分，得到R个距离单元，将探测多普勒范围-f_r/2～f_r/2进行等间隔划分，得到L个多普勒单元，其中，r_min表示探测区域的最小探测距离，r_max表示探测区域的最大探测距离，f_r是脉冲重复频率；

(2b)构造维数为(N+R-1)×R的距离维基矩阵D：

(2b1)以最小探测距离r_min作为参考距离，不考虑多普勒项，得到位于第r个距离单元处的目标的回波为d_r＝[0_r-1s₁s₂…s_N0_R-r]^T，其中，0_r-1表示长度为r-1的零向量，0_R-r表示长度为R-r的零向量，r＝1,2,...,R，当r＝1时，0_r-1＝0₀，当r＝R时，0_R-r＝0₀，此时表示d_r中无此项；

(2b1)以d_r为列得出维数为(N+R-1)×R的距离维基矩阵D；

(2c)构造维数为L×Q的多普勒域基矩阵Θ：

(2c1)对步骤(1)得到的回波信号y(q)中的多普勒相移进行离散化处理，得到离散化后的多普勒相移进而得到第q个回波的第l个多普勒相移项为其中f_l＝-f_r/2+f_r(l-1/2)/L，l＝1,2,...,L，q＝1,2,...,Q，B是子脉冲的带宽，f_r是脉冲重复频率；

(2c2)以θ(l,q)为元素得出维数为L×Q的多普勒域基矩阵Θ；

(2d)由步骤(2b)中的距离维基矩阵D，步骤(2c)中的多普勒域基矩阵Θ得到维度为(N+R-1)×Q回波矩阵：Y＝DXΘ+E，其中，X为散射系数矩阵，E为噪声矩阵。

步骤3：引入辅助矩阵，构造优化求解模型。

(3a)引入辅助矩阵Z，令Z＝XΘ，将由步骤(2d)中得到的回波矩阵Y变形为下式：

$\left\{\begin{array}{c}Y=DZ+E\\ Z^T={\mathrm{\Θ}}^T{X}^T\end{array}\right.;$ 其中(·)^T表示矩阵转置；

(3b)将步骤(3a)中得到的第一个方程Y＝DZ+E记作下式：

y_i＝Dz_i+e_i,i＝1,…,Q

其中，y_i为Y的第i列，z_i为Z的第i列，e_i为E的第i列，Q为回波矩阵Y的列数；

步骤4：对优化函数进行求解，得到场景的距离-多普勒像。

传统的求解方法有自适应迭代算法，基于稀疏信号重构的平滑零范数法等，本发明采用的是多任务贝叶斯压缩感知算法进行求解：

(4a)对步骤(3b)中得到的式子进行求解，得到辅助变量矩阵Z，本步骤采用现有的多任务贝叶斯压缩感知算法进行求解：

(4a1)设置第一辅助参数a＝1，第二辅助参数b＝1，超参数矩阵A的初始值A⁰为全1矩阵，初始循环迭代次数k＝1；

(4a2)循环执行以下步骤，得到均值矩阵μ：

(4a2.1)按如下公式计算方差矩阵Σ的第k次迭代矩阵Σ^k：

Σ^k＝(D^TD+A^k-1)^-1

其中，(·)^-1表示矩阵取逆；

(4a2.2)按如下公式计算均值矩阵μ的第i列μ_i的第k次迭代向量

${\mathrm{\μ}}_i^k={\mathrm{\Σ}}_i^k{D}_i^T{y}_i,$

其中，Σ_i表示方差矩阵Σ的第i列，D_i表示距离维基矩阵D的第i列，y_i表示回波矩阵Y的第i列，i＝1,2,…,Q；

(4a2.3)按如下公式计算超参数矩阵A的元素α_j的第k次迭代值

${\mathrm{\α}}_j^k=\frac{Q}{\underset{i=1}{\overset{Q}{\Σ}}{\left({\mathrm{\μ}}_{i,j}^k\right)}^2(N+R-1+2a)/(y_i^T{\left(B^{k-1}\right)}^{-1}{y}_i+2b)+{\mathrm{\Σ}}_{j,j}^k}$

其中，μ_i,j表示均值矩阵μ的第i行第j列的元素值，Σ_j,j表示方差矩阵Σ的第j行第j列的元素值，B表示辅助矩阵，B＝I+D(A^k-1)^-1D^T，A＝diag(α₁,α₂,...,α_j,…,α_Q)，diag是构造对角矩阵的函数，N表示子脉冲个数，R表示距离单元个数，j＝1,2,...,Q；

(4a2.4)由步骤(4a2.3)中得到的超参数矩阵A的元素α_j的第k次迭代值组成新的超参数矩阵： $A^k=diag({\mathrm{\α}}_1^k,{\mathrm{\α}}_2^k,...,{\mathrm{\α}}_j^k,...,{\mathrm{\α}}_Q^k);$

(4a2.5)判断是否满足若满足，则停止迭代，并将均值矩阵μ的第i列μ_i的第k次迭代向量作为最终迭代结果；否则，更新迭代次数k＝k+1，重复执行步骤(4a2.1)到(4a2.5)，其中，eps为迭代终止值，eps＝10^-5，表示矩阵Frobenius范数的平方，i＝1,2,...,Q；

(4a3)用均值矩阵μ作为辅助变量矩阵Z的估计：

$Z_i={\mathrm{\μ}}_i^k={\mathrm{\Σ}}_i^k{D}_i^T{y}_i$

其中，Z_i表示辅助变量矩阵Z的第i列；

(4b)设定检测门限Δ＝0.1，筛选出辅助变量矩阵Z中的非零行组成新矩阵：Z_select＝X_selectΘ，对该式进行转置，得到：

$Z_{select}^T={\mathrm{\Θ}}^T{X}_{select}^T$

其中，Z_select为辅助变量矩阵Z剔除全零行后的矩阵，X_select为散射系数矩阵X剔除全零行后的矩阵；

(4c)将步骤(4b)中方程记作：z_i＝Θ^Tx_i,i＝1,…,R，其中，z_i，为的第i列，x_i为的第i列；

(4d)对步骤(4c)中的式子z_i＝Θ^Tx_i进行求解，得到散射系数矩阵X剔除全零行并进行转置后的矩阵本步骤采用现有的多任务贝叶斯压缩感知算法进行求解：

(4d1)设置第一辅助参数a＝1，第二辅助参数b＝1，超参数矩阵G的初始值G⁰为全1矩阵，初始循环迭代次数k＝1；

(4d2)循环执行以下步骤，得到均值矩阵η：

(4d2.1)按如下公式计算方差矩阵Ω的第k次迭代矩阵Ω^k：

Ω^k＝(ΦΦ^T+G^k-1)^-1

其中，(·)^-1表示矩阵取逆；

(4d2.2)按如下公式计算均值矩阵η的第i列η_i的第k次迭代向量

${\mathrm{\η}}_i^k={\mathrm{\Ω}}_i^k{\mathrm{\Φ}}_i{z}_i,$

其中，Ω_i表示方差矩阵Ω的第i列，Φ_i表示多普勒域基矩阵Φ的第i列，i＝1,2,…,R；

(4d2.3)按如下公式计算超参数矩阵G的元素β_j的第k次迭代值

${\mathrm{\β}}_j^k=\frac{R}{\underset{i=1}{\overset{R}{\Σ}}{\left({\mathrm{\η}}_{i,j}^k\right)}^2(Q+2a)/(z_i^T{\left(W^{k-1}\right)}^{-1}{z}_i+2b)+{\mathrm{\Ω}}_{j,j}^k}$

其中，η_i,j表示均值矩阵η的第i行第j列的元素值，Ω_j,j表示方差矩阵Ω的第j行第j列的元素值，W表示辅助矩阵，W＝I+Φ^T(G^k-1)^-1Φ，G＝diag(β₁,β₂,...,β_j,…,β_L)，diag是构造对角矩阵的函数，j＝1,2,...,L；

(4d2.4)由步骤(4d2.3)中得到的超参数矩阵G的元素β_j的第k次迭代值组成新的超参数矩阵： $G^k=diag({\mathrm{\β}}_1^k,{\mathrm{\β}}_2^k,...,{\mathrm{\β}}_j^k,...,{\mathrm{\β}}_L^k);$

(4d2.5)判断是否满足若满足，则停止迭代，并将均值矩阵η的第i列η_i的第k次迭代向量作为最终迭代结果；否则，更新迭代次数k＝k+1，重复执行步骤(4d2.1)到(4d2.5)，其中，eps为迭代终止值，eps＝10^-5，表示矩阵Frobenius范数的平方，i＝1,2,...,L；

(4d3)用均值矩阵η作为的估计：

$x_i={\mathrm{\η}}_i^k={\mathrm{\Ω}}_i^k{\mathrm{\Φ}}_i{z}_i$

其中，x_i表示的第i列；

(4d4)对补全零行并取转置，得到散射系数矩阵X，由散射系数矩阵X即可得到距离-多普勒像。

本发明的效果可以通过下述仿真实验加以说明：

1.仿真条件

(1a)运行平台配置：

CPU：Intel(R)Core(TM)i7-3770CPU3.40GHz；内存：8GB；

操作系统：Windows7旗舰版64位操作系统；

仿真软件：MATLABR2012a。

(1b)仿真参数设置

发射信号参数以及实验仿真参数设置如表1所示。场景的距离向网格个数R＝40，多普勒向网格个数L＝80。

表1发射信号参数以及实验仿真参数

参数	参数值
		Barker序列	[1 1 1 1 1 -1 -1 1 1 -1 1 -1 1]T
脉冲个数	16
		脉冲重复频率	2KHz
子脉冲带宽	10MHz
		载波频率	1GHz
信噪比	10db
		目标个数	5
目标距离	15、15、20、20、22
		目标速度	-37.5、22.5、-7.5、0、0

2.仿真内容与结果

仿真1，根据表1的仿真参数，用现有的稀疏成像法对目标进行距离-多普勒成像，得到目标的距离-多普勒二维谱如图3所示，得到目标的距离与速度如表2所示：

表2

距离	15	15	20	22
					速度	-37.5	22.5	-3.75	0

仿真2，根据表1的仿真参数，用本发明方法对目标进行距离-多普勒成像，得到目标的距离-多普勒二维谱如图4所示，得到目标的距离与速度如表3所示：

表3

距离	15	15	20	20	22
						速度	-37.5	22.5	-7.5	0	0

3.仿真结果分析

由图3、图4、表2和表3的比较可以看出，用本发明方法可以有效地分辨出相距较近，速度也相近的目标，即具有更高的分辨率。

Claims (4)

1.基于多任务贝叶斯压缩感知的距离-多普勒成像方法，包括如下步骤：

(1)雷达向目标场景H发射Q个脉冲调制序列并在发射脉冲的同时接收该场景H的回波脉冲，得到第q个脉冲的回波信号为y(q)：

其中，x_r,l是距离为r多普勒频率为l的目标的散射强度，τ_r为延迟时间，为多普勒相移，q＝1,2,…,Q，R是距离单元个数，L是多普勒单元个数；

(2)对回波信号y(q)进行离散化采样，构造回波矩阵：Y＝DXΘ+E，其中，D为维数为(N+R-1)×R的距离维基矩阵，X为维数为R×L的散射系数矩阵，Θ为维数为L×Q的多普勒域基矩阵，E为噪声矩阵；

(3)引入辅助变量矩阵Z，令Z＝XΘ，将回波矩阵变形为下式：

其中(·)^T表示矩阵转置；

(4)将步骤(3)中第一个方程记作下式：

y_i＝Dz_i+e_i,i＝1,…,Q

其中，y_i表示回波矩阵Y的第i列，z_i表示辅助变量矩阵Z的第i列，e_i为噪声矩阵E的第i列，Q为回波矩阵Y的列数；

(5)对步骤(4)中的式子进行求解，得到辅助变量矩阵Z，设定检测门限Δ＝0.1，筛选出辅助变量矩阵Z中的非零行组成新矩阵：Z_select＝X_selectΘ，对该式进行转置，得到：

其中，X_select为散射系数矩阵X剔除全零行后的矩阵；

(6)将步骤(5)中方程记作：z_i＝Θ^Tx_i,i＝1,…,R，其中，z_i为，的第i列，x_i为的第i列；

(7)对步骤(6)中的式子进行求解，得到散射系数矩阵X剔除全零行并进行转置后的矩阵再对补全零行并取转置，得到散射系数矩阵X，由散射系数矩阵X即可得到距离-多普勒像。

2.根据权利要求1所述的基于多任务贝叶斯压缩感知的距离-多普勒成像方法，其中步骤(2)中对回波信号y(q)进行离散化采样，构造回波矩阵Y，按如下步骤进行：

(2a)在距离-多普勒二维探测区域将探测距离范围r_min～r_max进行等间隔划分，得到R个距离单元，将探测多普勒范围-f_r/2～f_r/2进行等间隔划分，得到L个多普勒单元，其中，r_min表示探测区域的最小探测距离，r_max表示探测区域的最大探测距离，f_r是脉冲重复频率；

(2b)构造维数为(N+R-1)×R的距离维基矩阵D：

以最小探测距离r_min作为参考距离，不考虑多普勒项，则位于第r个距离单元处的目标的回波为d_r＝[0_r-1s₁s₂…s_N0_R-r]^T，其中，0_r-1表示长度为r-1的零向量，0_R-r表示长度为R-r的零向量，r＝1,2,...,R，当r＝1时，0_r-1＝0₀，当r＝R时，0_R-r＝0₀，此时表示d_r中无此项。

以d_r为列得出距离维基矩阵D；

(2c)构造维数为L×Q的多普勒域基矩阵Θ：

对步骤(1)得到的回波信号y(q)中的多普勒相移进行离散化处理，得到离散化后的多普勒相移进而得到第q个回波的第l个多普勒相移项为其中f_l＝-f_r/2+f_r(l-1/2)/L，l＝1,2,...,L，q＝1,2,...,Q，B是子脉冲的带宽，f_r是脉冲重复频率。

以为元素得出多普勒域基矩阵Θ；

(2d)由步骤(2b)中的距离维基矩阵D，步骤(2c)中的多普勒域基矩阵Θ得到维度为(N+R-1)×Q回波矩阵：Y＝DXΘ+E。

3.根据权利要求1所述的基于多任务贝叶斯压缩感知的距离-多普勒成像方法，其中步骤(5)中对步骤(4)中的式子进行求解，得到辅助变量矩阵Z，依照如下步骤进行：

(5a)设置第一辅助参数a＝1，第二辅助参数b＝1，超参数矩阵A的初始值A⁰为全1矩阵，初始循环迭代次数k＝1；

(5b)循环执行以下步骤，得到均值矩阵μ：

(5b1)按如下公式计算方差矩阵Σ的第k次迭代向量Σ^k：

Σ^k＝(D^TD+A^k-1)^-1

其中，(·)^-1表示矩阵取逆；

(5b2)按如下公式计算均值矩阵μ的第i列μ_i的第k次迭代向量

其中，Σ_i表示方差矩阵Σ的第i列，D_i表示距离维基矩阵D的第i列，y_i表示回波矩阵Y的第i列，i＝1,2,…,Q；

(5b3)按如下公式计算超参数矩阵A的元素α_j的第k次迭代值

其中，μ_i,j表示均值矩阵μ的第i行第j列的元素值，Σ_j,j表示方差矩阵Σ的第j行第j列的元素值，B表示辅助矩阵，B＝I+D(A^k-1)^-1D^T，A＝diag(α₁,α₂,...,α_j,…,α_Q)，diag是构造对角矩阵的函数，j＝1,2,...,Q；

(5b4)由步骤(5b3)得到的超参数矩阵A的元素α_j的第k次迭代值组成新的超参数矩阵：

(5b5)判断是否满足若满足，则停止迭代，并将均值矩阵μ的第i列μ_i的第k次迭代向量作为最终迭代结果；否则，更新迭代次数k＝k+1，重复执行步骤(5b1)到(5b5)，其中，eps为迭代终止值，eps＝10^-5，表示矩阵Frobenius范数的平方，i＝1,2,...,Q；

(5c)用均值矩阵μ作为辅助变量矩阵Z的估计：

其中，Z_i表示辅助变量矩阵Z的第i列。

4.根据权利要求1所述的基于多任务贝叶斯压缩感知的距离-多普勒成像方法，其中对步骤(7)中步骤(6)中的式子进行求解，得到散射系数矩阵X剔除全零行并进行转置后的矩阵按如下步骤进行：

(7a)设置第一辅助参数a＝1，第二辅助参数b＝1，超参数矩阵G的初始值G⁰为全1矩阵，初始循环迭代次数k＝1；

(7b)循环执行以下步骤，得到均值矩阵η：

(7b1)按如下公式计算方差矩阵Ω的第k次迭代矩阵Ω^k：

Ω^k＝(ΦΦ^T+G^k-1)^-1

其中，(·)^-1表示矩阵取逆；

(7b2)按如下公式计算均值矩阵η的第i列η_i的第k次迭代向量

其中，Ω_i表示方差矩阵Ω的第i列，Φ_i表示多普勒域基矩阵Φ的第i列，i＝1,2,…,R；

(7b3)按如下公式计算超参数矩阵G的元素β_j的第k次迭代值

其中，η_i,j表示均值矩阵η的第i行第j列的元素值，Ω_j,j表示方差矩阵Ω的第j行第j列的元素值，W表示辅助矩阵，W＝I+Φ^T(G^k-1)^-1Φ，G＝diag(β₁,β₂,...,β_j,…,β_L)，diag是构造对角矩阵的函数，j＝1,2,...,L；

(7b4)由步骤(7b3)中得到的超参数矩阵G的元素β_j的第k次迭代值组成新的超参数矩阵：

(7b5)判断是否满足若满足，则停止迭代，并将均值矩阵η的第i列η_i的第k次迭代向量作为最终迭代结果；否则，更新迭代次数k＝k+1，重复执行步骤(7b1)到(7b5)，其中，eps为迭代终止值，eps＝10^-5，表示矩阵Frobenius范数的平方，i＝1,2,...,L；

(7c)用均值矩阵η作为的估计：

其中，x_i表示的第i列。