CN103983948A - 基于稀疏表示的目标角度距离联合估计方法 - Google Patents

Abstract

本发明公开一种基于稀疏表示的目标角度-距离联合估计方法，主要解决现有技术计算量大的问题。其技术方案是：1.对阵元接收信号进行离散化，构造接收信号矩阵；2.将探测区域的角度范围与距离范围分别进行等间隔划分，构造基于角度的稀疏基和基于距离的稀疏基；3.根据接收信号矩阵、基于角度的稀疏基和基于距离的稀疏基，建立稀疏优化模型，并用二维交替优化的方法对稀疏优化模型进行求解，得到最终优化解；4.通过对最终优化解的峰值搜索，得到目标的角度与距离。本发明具有计算量小和估计精确的优点，可用于雷达、声呐中的目标探测。

Description

基于稀疏表示的目标角度距离联合估计方法

技术领域

本发明属于信号处理技术领域，特别涉及一种目标估计方法，可以用于雷达、声呐中的目标探测。

背景技术

相控阵是利用电磁波的相干原理，通过计算机控制馈往各辐射阵元电流的相位，从而改变波束方向的阵列天线。相控阵采用电子扫描方式，可实现回波数据的实时更新，因此获得了广泛的关注。其中利用相控阵进行目标位置估计是相控阵应用的一个主要方面。

目前，目标位置的估计方法已有很多，例如多重信号分类MUSIC方法、旋转不变技术估计信号参数ESPRIT方法、基于稀疏表示的方法等，其中：

多重信号分类MUSIC方法，是将阵列接收信号的自相关矩阵进行特征值分解，然后利用信号子空间和噪声子空间之间的正交关系来估计目标的位置。例如，孙洪亮的硕士论文“近场源参数联合估计算法仿真研究”，就是用MUSIC方法进行目标位置的估计，这种方法的最大不足是分辨率不高。

旋转不变技术估计信号参数ESPRIT方法，是利用两个对称子阵自相关和互相关矩阵的广义特征值进行目标位置的估计。例如，张群飞，保铮，黄建国的论文“一种水下多目标方位、频率、距离联合估计新方法”(电子学报，Vol.32No.9Sep.2004)，就是用ESPRIT方法进行目标位置的估计，这种方法的最大不足是在噪声较大的情况下，该方法的精度不是很高。

基于稀疏表示的目标位置估计方法，是将探测区域进行网格化，利用目标在网格化探测区域中具有稀疏性的假设，将目标位置估计问题转化成目标具有稀疏性假设的优化求解问题。例如，王鹏的硕士论文“阵列天线雷达的目标方位距离联合估计和高分辨成像”，就是一种基于稀疏表示的方法，这种方法的最大不足是在稀疏优化模型中，目标的角度距离双参量是构造在一个稀疏域里，若网格划分过密，会导致该稀疏域维数巨大，使得优化求解的计算量很大，工程应用比较困难。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于稀疏表示的目标角度距离联合估计的方法，以减少求解计算量，便于实现工程应用。

本发明的思路是：将角度距离双参量稀疏基的联合构造转变为各参量对应的稀疏基级联。其实现步骤包括如下：

(1)将探测区域进行离散化，即将距离范围r_min～r_max和角度范围θ_min～θ_max分别进行等间隔划分，共得到G个距离单元r_g和Q个角度单元θ_q：

$\left\{\begin{array}{c}{r}_g=r_{\min }+\frac{g-1}{G-1}(r_{\max }-r_{\min })\\ {\mathrm{\θ}}_q={\mathrm{\θ}}_{\min }+\frac{q-1}{Q-1}({\mathrm{\θ}}_{\max }-{\mathrm{\θ}}_{\min })\end{array}\right.$

其中，r_min和r_max分别表示探测区域的最小距离和最大距离，θ_min和θ_max分别表示探测区域的最小角度和最大角度，g＝1,2,…，G，q＝1,2,…，Q；

(2)用M个阵元组成阵列，并将第一个阵元作为参考阵元，设参考阵元发射的窄带信号为t表示时间，p(t)是基带信号，f₀是信号载频，j是虚数单位，M≥2；

(3)设探测目标在离散化的探测区域中，则第m个阵元接收到目标反射信号为x_m(t)，m＝1,2,…，M；对目标反射信号x_m(t)进行离散采样，得到离散数据x_m(t_n)，n＝1,2,…，N，N是采样点数；并以该离散数据x_m(t_n)作为第m行，构成阵列信号接收矩阵X：

$X=\left[\begin{array}{c}{x}_1\left(t_n\right)\\ \·\\ \·\\ \·\\ x_m\left(t_n\right)\\ \·\\ \·\\ \·\\ x_M\left(t_n\right)\end{array}\right];$

(4)根据Q个角度单元θ_q和G个距离单元r_g，构造基于Q个角度单元的稀疏基Φ和基于G个距离单元的稀疏基Ψ；

(5)用距离单元的稀疏基Ψ和角度单元的稀疏基Φ作为约束条件，引入辅助矩阵Y，采用二维交替优化的方法，求解下述优化式，得到目标角度距离二维谱矩阵A：

$\underset{Y,F}{\min }\{\mathrm{\λ}{\left|\right|Y\left|\right|}_{\mathrm{2,1}}+{\left|\right|A\left|\right|}_{\mathrm{2,1}}\},s.t.\left\{\begin{array}{c}X=\mathrm{\ΦY}\\ Y^T=\mathrm{\ΨA}\end{array}\right.$

其中，λ是用户输入的正则化参量，Y^q表示辅助矩阵Y的第q行，q＝1,2,…，Q，A^g表示角度距离二维谱矩阵A的第g行，g＝1,2,…，G，||||₂表示求向量的2范数，(·)^T表示转置；

(6)根据步骤(5)中求的角度距离二维谱矩阵A，采用阈值比较法，得到目标的距离位置a与角度位置b，并通过下式确定目标的距离r_a与角度θ_b：

$\left\{\begin{array}{c}{r}_a=r_{\min }+\frac{a-1}{G-1}(r_{\max }-r_{\min })\\ {\mathrm{\θ}}_b={\mathrm{\θ}}_{\min }+\frac{b-1}{Q-1}({\mathrm{\θ}}_{\max }-{\mathrm{\θ}}_{\min })\end{array}\right..$

本发明与常规方法相比具有以下优点：

第一，由于本发明是将距离角度双参量稀疏基的联合构造转变为各参量对应的稀疏基级联，降低了优化求解的计算量，求解所用时间少，便于工程应用。

第二，由于本发明是将距离角度双参量稀疏基的联合构造转变为各参量对应的稀疏基级联，可以将探测区域划分的更密，最后得到目标的距离与角度更精确。

附图说明

图1是探测区域的离散化图；

图2是本发明的流程图；

图3是用现有技术获得的目标角度距离二维谱的仿真结果图；

图4是用本发明方法获得的目标角度距离二维谱的仿真结果图；

图5是用现有技术和本发明方法进行目标角度距离联合估计所用时间对比图。

具体实施方式

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

参照图2，本发明的具体实施步骤如下：

步骤一：构造接收信号矩阵。

1.1)探测区域进行离散化

参照图1，将探测区域的距离范围r_min～r_max和角度范围θ_min～θ_max分别进行等间隔划分，得到G个距离单元r_g和Q个角度单元θ_q：

$\left\{\begin{array}{c}{r}_g=r_{\min }+\frac{g-1}{G-1}(r_{\max }-r_{\min })\\ {\mathrm{\θ}}_q={\mathrm{\θ}}_{\min }+\frac{q-1}{Q-1}({\mathrm{\θ}}_{\max }-{\mathrm{\θ}}_{\min })\end{array}\right.,$

其中，r_min和r_max分别表示探测区域的最小距离和最大距离，θ_min和θ_max分别表示探测区域的最小角度和最大角度，g＝1,2,…,G，q＝1,2,…,Q；

1.2)用M个阵元组成阵列，并将第一个阵元作为参考阵元，设参考阵元发射的窄带信号为t表示时间，p(t)是基带信号，j是虚数单位，f₀是信号载频，M≥2；

1.3)设离散化的探测区域中共有L个探测目标，则第m个阵元接收到目标反射信号x_m(t)为:

$x_m\left(t\right)=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}s(t-{\mathrm{\τ}}_{\mathrm{ml}}),$

其中，τ_ml是第l个目标与第m个阵元之间的距离引起的信号传播相对时间延迟，m＝1,2,…,M，l＝1,2,…,L，L远小于G×Q；

1.4)对目标反射信号x_m(t)进行离散采样，得到离散数据并以该离散数据x_m(t_n)作为第m行，构成接收信号矩阵X：

$X=\left[\begin{array}{c}{x}_1\left(t_n\right)\\ \·\\ \·\\ \·\\ x_m\left(t_n\right)\\ \·\\ \·\\ \·\\ x_M\left(t_n\right)\end{array}\right],$

其中，n＝1,2,…,N，N是采样点数；

步骤二：构造稀疏基。

根据Q个角度单元θ_q和G个距离单元r_g，按如下公式构造基于Q个角度单元的稀疏基Φ和基于G个距离单元的稀疏基Ψ：

其中，d为阵列中相邻阵元的间距，c是信号s(t)的传播速度，[·]^T表示转置，j是虚数单位。

步骤三：建立稀疏优化模型并求解。

3a)用距离单元的稀疏基Ψ和角度单元的稀疏基Φ作为约束条件，引入辅助矩阵Y，构建优化式：

$\underset{Y,F}{\min }\{\mathrm{\λ}{\left|\right|Y\left|\right|}_{\mathrm{2,1}}+{\left|\right|A\left|\right|}_{\mathrm{2,1}}\},s.t.\left\{\begin{array}{c}X=\mathrm{\ΦY}\\ Y^T=\mathrm{\ΨA}\end{array}\right.,$

其中，λ是用户输入的正则化参量，Y^q表示辅助矩阵Y的第q行，q＝1,2,…,Q，A^g表示角度距离二维谱矩阵A的第g行，g＝1,2,…,G，||||₂表示求向量的2范数，(·)^T表示转置；

3b)采用二维交替优化的方法，求解上述优化式，得到目标角度距离二维谱矩阵A：

3b1)设置角度距离二维谱矩阵A的初值A₀、辅助矩阵Y的初值Y₀、两个拉格朗日乘子Λ₁和Λ₂的初值为全1矩阵，输入正则化因子λ，初始化与第一个拉格朗日乘子Λ₁对应的惩罚因子β₁＞0，初始化与第二个拉格朗日乘子Λ₂对应的惩罚因子β₂＞0，初始化与第一个拉格朗日乘子Λ₁对应的步长0＜γ₁＜1.618，初始化与第二个拉格朗日乘子Λ₂对应的步长0＜γ₂＜1.618，循环迭代次数k＝1；

3b2)循环执行以下步骤：

3b2.1)按如下公式计算与角度距离二维谱矩阵A的第k-1次迭代值A_k-1对应的对角矩阵Θ₁：

${\mathrm{\Θ}}_1=\mathrm{diag}(1/{\left|\right|A_{k-1}^g\left|\right|}_2),$

其中，diag是构造对角矩阵的函数，是A_k-1的第g行，g＝1,2,…,G，G是距离单元的个数；

3b2.2)按如下公式计算角度距离二维谱矩阵A的第k次迭代值A_k：

A_k＝(Θ₁+β₂Ψ^HΨ)^-1(β₂Ψ^H(Y_k-1)^T+Ψ^HΛ₂)，

其中，(·)^H表示共轭转置，(·)^-1表示矩阵求逆，(·)^T表示转置；

3b2.3)按如下公式计算与辅助矩阵Y的第k-1次迭代值Y_k-1对应的对角矩阵Θ₂：

${\mathrm{\Θ}}_2=\mathrm{diag}(1/{\left|\right|Y_{k-1}^q\left|\right|}_2),$

其中，diag是构造对角矩阵的函数，是Y_k-1的第q行，q＝1,2,…,Q，Q是角度单元的个数；

3b2.4)按如下公式计算辅助矩阵Y的第k次迭代值Y_k：

$Y_k={({\mathrm{\λ\Θ}}_2+{\mathrm{\β}}_1{\mathrm{\Φ}}^H\mathrm{\Φ}+{\mathrm{\β}}_{1}I)}^{-1}({\mathrm{\Φ}}^H{\mathrm{\Λ}}_1+{\mathrm{\Λ}}_2^T+{\mathrm{\β}}_1{\mathrm{\Φ}}^{H}X+{\mathrm{\β}}_2{\left(A_{k-1}\right)}^T{\mathrm{\Ψ}}^T),$

其中，(·)^H表示共轭转置，(·)^-1表示矩阵求逆，(·)^T表示转置，I是Q×Q的单位矩阵；

3b2.5)按如下公式更新两个拉格朗日乘子Λ₁和Λ₂：

$\left\{\begin{array}{c}{\mathrm{\Λ}}_1={\mathrm{\Λ}}_1-{\mathrm{\γ}}_1{\mathrm{\β}}_1({\mathrm{\ΦY}}_k-X)\\ {\mathrm{\Λ}}_2={\mathrm{\Λ}}_2-{\mathrm{\γ}}_2{\mathrm{\β}}_2({\mathrm{\ΨA}}_k-{\left(Y_k\right)}^T)\end{array}\right.,$

3b2.6)判断是否满足迭代终止值δ＝10^-5，若满足，则停止迭代；否则，继续迭代，更新迭代次数k＝k+1，直到满足其中，表示矩阵Frobenius范数的平方；

迭代最终得到的角度距离二维谱矩阵A_k，就是本发明要求的角度距离二维谱矩阵A。

步骤四：计算目标的角度与距离。

4a)令向量表示转置，A^g是角度距离二维谱矩阵A的第g行，g＝1,2,…,G，G是距离单元的个数；

4b)对向量进行归一化处理，得到归一化向量A；

4c)设置检测阈值ε＝0.1，按下式得到归一化向量A峰值元素索引值z：

z＝{i|A_i＞ε,i＝1,2,…,G×Q}，

其中，A_i是归一化向量A的第i个元素，G是距离单元的个数，Q是角度单元的个数；

4d)根据峰值元素索引值z和距离单元的个数G，计算目标的角度位置：b＝z％G，％表示取余；

4e)根据目标的角度位置b、峰值元素索引值z和角度单元个数Q，计算目标的距离位置： $a=\frac{z-b}{Q}+1;$

4f)根据目标的距离位置a、距离单元个数G、探测区域的最小距离_rmin和最大距离r_max，计算目标的距离： $r_a=r_{\min }+\frac{a-1}{G-1}(r_{\max }-r_{\min });$

4g)根据目标的角度位置b、角度单元个数Q、探测区域的最小角度θ_min和最大角度θ_max，计算目标的角度： ${\mathrm{\θ}}_b={\mathrm{\θ}}_{\min }+\frac{b-1}{Q-1}({\mathrm{\θ}}_{\max }-{\mathrm{\θ}}_{\min }).$

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

1、仿真条件

本实验的硬件测试平台是：Intel Core i5CPU，主频3.10GHz，内存8.00GB；软件仿真平台是：64位windows8操作系统和Matlab R(2013b)，仿真参数设置如下表1所示：

表1仿真参数

参数	参数值
		系统载频	300MHz
阵元个数	16
		阵元间距	0.5m
时间采样频率	600MHz
		时间采样点数	276
信噪比	10dB
		目标个数	3
目标角度	25、25、30
		目标距离	1020、1030、1020

2.仿真内容与结果

仿真1，用现有技术进行目标角度距离联合估计，得到目标的角度距离二维谱如图2所示，搜索峰值得到目标的角度与距离如下表2所示：

表2

角度	24	24	30
				距离	1030	1020	1020

仿真2，用本发明方法进行目标角度距离联合估计，得到目标的角度距离二维谱如图3所示，搜索峰值得到目标的角度与距离如下表3所示：

表3

角度	25	25	30
				距离	1030	1020	1020

从图2、图3、表2和表3的比较可以看出，用本发明方法得到的目标距离与角度更精确。

仿真3，分别用现有技术和本发明方法进行目标角度距离联合估计，得到两种方法随采样点数的时间对比图，仿真结果图如图4所示。

从图4可以看出，随着采样点数的增大，本发明方法所用时间远远小于现有技术，说明本发明方法降低了求解计算量。

Claims (4)

1.一种基于稀疏表示的目标角度距离联合估计方法，包括如下步骤：

(1)将探测区域进行离散化，即将距离范围r_min～r_max和角度范围θ_min～θ_max分别进行等间隔划分，共得到G个距离单元r_g和Q个角度单元θ_q：

$\left\{\begin{array}{c}{r}_g=r_{\min }+\frac{g-1}{G-1}(r_{\max }-r_{\min })\\ {\mathrm{\θ}}_q={\mathrm{\θ}}_{\min }+\frac{q-1}{Q-1}({\mathrm{\θ}}_{\max }-{\mathrm{\θ}}_{\min })\end{array}\right.,$

其中，r_min和r_max分别表示探测区域的最小距离和最大距离，θ_min和θ_max分别表示探测区域的最小角度和最大角度，g＝1,2,…，G，q＝1,2,…，Q；

(2)用M个阵元组成阵列，并将第一个阵元作为参考阵元，设参考阵元发射的窄带信号为t表示时间，p(t)是基带信号，j是虚数单位，f₀是信号载频，M≥2；

(3)设探测目标在离散化的探测区域中，则第m个阵元接收到目标反射信号为x_m(t)，m＝1,2,…，M；对目标反射信号x_m(t)进行离散采样，得到离散数据x_m(t_n)，n＝1,2,…,N，N是采样点数；并以该离散数据x_m(t_n)作为第m行，构成阵列信号接收矩阵X：

$X=\left[\begin{array}{c}{x}_1\left(t_n\right)\\ \·\\ \·\\ \·\\ x_m\left(t_n\right)\\ \·\\ \·\\ \·\\ x_M\left(t_n\right)\end{array}\right];$

(4)根据Q个角度单元θ_q和G个距离单元r_g，构造基于Q个角度单元的稀疏基Φ和基于G个距离单元的稀疏基Ψ；

(5)用距离单元的稀疏基Ψ和角度单元的稀疏基Φ作为约束条件，引入辅助矩阵Y，采用二维交替优化的方法，求解下述优化式，得到目标角度距离二维谱矩阵A：

$\underset{Y,F}{\min }\{\mathrm{\λ}{\left|\right|Y\left|\right|}_{\mathrm{2,1}}+{\left|\right|A\left|\right|}_{\mathrm{2,1}}\},s.t.\left\{\begin{array}{c}X=\mathrm{\ΦY}\\ Y^T=\mathrm{\ΨA}\end{array}\right.$

其中，λ是用户输入的正则化参量，Y^q表示辅助矩阵Y的第q行，q＝1,2，…，Q，A^g表示角度距离二维谱矩阵A的第g行，g＝1,2,…G，||||₂表示求向量的2范数，(·)^T表示转置；

(6)根据步骤(5)中求得的角度距离二维谱矩阵A，采用阈值比较法，得到目标的距离位置a与角度位置b，并通过下式确定目标的距离r_a与角度θ_b：

$\left\{\begin{array}{c}{r}_a=r_{\min }+\frac{a-1}{G-1}(r_{\max }-r_{\min })\\ {\mathrm{\θ}}_b={\mathrm{\θ}}_{\min }+\frac{b-1}{Q-1}({\mathrm{\θ}}_{\max }-{\mathrm{\θ}}_{\min })\end{array}\right..$

2.根据权利要求1所述的基于稀疏表示的目标角度距离联合估计方法，其中所述步骤(4)中根据Q个角度单元θ_q和G个距离单元r_g，构造基于Q个角度单元的稀疏基Φ和基于G个距离单元的稀疏基Ψ，按如下公式进行：

其中，d为阵列中相邻阵元的间距，c是信号s(t)的传播速度，[·]^T表示转置，j是虚数单位。

3.根据权利要求1所述的基于稀疏表示的目标角度距离联合估计方法,其中所述步骤(5)中采用二维交替优化的方法，求解优化式，得到目标角度距离二维谱矩阵A，按如下步骤进行：

(5a)设置角度距离二维谱矩阵A的初值A₀、辅助矩阵Y的初值Y₀、两个拉格朗日乘子Λ₁和Λ₂的初值为全1矩阵，输入正则化因子λ，初始化与第一个拉格朗日乘子Λ₁对应的惩罚因子β₁＞0，初始化与第二个拉格朗日乘子Λ₂对应的惩罚因子β₂＞0，初始化与第一个拉格朗日乘子Λ₁对应的步长0＜γ₁＜1.618，初始化与第二个拉格朗日乘子Λ₂对应的步长0＜γ₂＜1.618，循环迭代次数k＝1；

(5b)循环执行以下步骤：

(5b1)按如下公式计算与角度距离二维谱矩阵A的第k-1次迭代值A_k-1对应的对角矩阵Θ₁：

${\mathrm{\Θ}}_1=\mathrm{diag}(1/{\left|\right|A_{k-1}^g\left|\right|}_2),$

其中，diag是构造对角矩阵的函数，是A_k-1的第g行，g＝1,2,…，G，G是距离单元的个数；

(5b2)按如下公式计算角度距离二维谱矩阵A的第k次迭代值A_k：

A_k＝(Θ₁+β₂Ψ^HΨ)^-1(β₂Ψ^H(Y_k-1)^T+Ψ^HΛ₂)，

其中，(·)^H表示共轭转置，(·)^-1表示矩阵求逆，(·)^T表示转置；

(5b3)按如下公式计算与辅助矩阵Y的第k-1次迭代值Y_k-1对应的对角矩阵Θ₂：

${\mathrm{\Θ}}_2=\mathrm{diag}(1/{\left|\right|Y_{k-1}^q\left|\right|}_2),$

其中，diag是构造对角矩阵的函数，是Y_k-1的第q行，q＝1，2，…，Q，Q是角度单元的个数；

(5b4)按如下公式计算辅助矩阵Y的第k次迭代值Y_k：

$Y_k={({\mathrm{\λ\Θ}}_2+{\mathrm{\β}}_1{\mathrm{\Φ}}^H\mathrm{\Φ}+{\mathrm{\β}}_{1}I)}^{-1}({\mathrm{\Φ}}^H{\mathrm{\Λ}}_1+{\mathrm{\Λ}}_2^T+{\mathrm{\β}}_1{\mathrm{\Φ}}^{H}X+{\mathrm{\β}}_2{\left(A_{k-1}\right)}^T{\mathrm{\Ψ}}^T),$

其中，(·)^H表示共轭转置，(·)^-1表示矩阵求逆，(·)^T表示转置，I是Q×Q的单位矩阵；

(5b5)按如下公式更新两个拉格朗日乘子Λ₁和Λ₂：

$\left\{\begin{array}{c}{\mathrm{\Λ}}_1={\mathrm{\Λ}}_1-{\mathrm{\γ}}_1{\mathrm{\β}}_1({\mathrm{\ΦY}}_k-X)\\ {\mathrm{\Λ}}_2={\mathrm{\Λ}}_2-{\mathrm{\γ}}_2{\mathrm{\β}}_2({\mathrm{\ΨA}}_k-{\left(Y_k\right)}^T)\end{array}\right.;$

(5b6)判断是否满足迭代终止值δ＝10^-5，若满足，则停止迭代；否则，继续迭代，更新迭代次数k＝k+1，直到满足其中，表示矩阵Frobenius范数的平方。

4.根据权利要求1所述的基于稀疏表示的目标角度距离联合估计方法，其中所述步骤(6)中采用阈值比较法，得到目标的距离位置a与角度位置b，按如下步骤进行：

(6a)令向量表示转置，A^g是角度距离二维谱矩阵A的第g行，g＝1,2,…，G，G是距离单元的个数；

(6b)对向量进行归一化处理，得到归一化向量A；

(6c)设置检测阈值ε＝0.1，按下式得到归一化向量A峰值元素索引值z：

z＝{i|A_i＞ε,i＝1,2,…，G×Q}，

其中，A_i是归一化向量A的第i个元素，G是距离单元的个数，Q是角度单元的个数；

(6d)按如下公式求得目标的角度位置b与距离位置a：

$\left\{\begin{array}{c}b=z\%G\\ a=\frac{z-b}{Q}+1\end{array}\right.,$

其中，％表示取余。