CN103744061A - 基于迭代最小二乘方法的mimo雷达doa估计方法 - Google Patents

Abstract

一种基于迭代最小二乘方法的MIMO雷达DOA估计方法，用迭代最小二乘方法求解经过降维处理的接收和发射阵列响应矩阵。首先对多个雷达发射脉冲的回波数据矩阵以及接收和发射阵列响应矩阵进行降维处理；然后在最小二乘条件下建立代价函数，并利用基于梯度下降的迭代方法求解代价函数；最后利用已知的收、发阵列流形估计目标方向。与传统的单基地MIMO雷达阵列DOA估计方法相比：本发明直接得到目标的DOA估计，不需要进行谱峰搜索；采用降维处理有效抑制了噪声，提高了低信噪比下的估计精度；避免了高维数据协方差矩阵的估计、求逆及特征值分解运算，克服了传统阵列DOA估计方法应用于单基地MIMO雷达时计算量大、需要样本数大的缺点。

Description

基于迭代最小二乘方法的MIMO雷达DOA估计方法

技术领域

本发明属于雷达信号处理技术领域，具体说是根据单基地MIMO（multiple-transmitmultiple-receive）雷达目标回波数据所包含的结构信息，运用迭代最小二乘（I-LS）方法对雷达目标进行波达方向（DOA,Direction-of-Arrival）估计的方法。

背景技术

从上世纪30年代第二次世界大战以来，现代雷达技术已经历了80多年突飞猛进的发展。21世纪以来，在MIMO无线通信理论所取得的巨大成功的推动下，一种新体制雷达——MIMO雷达逐渐成为雷达领域的研究热点。根据雷达天线在空间的分布，MIMO雷达可以分为集中式MIMO雷达和分布式MIMO雷达两类。其中，集中式MIMO雷达根据接收端和发射端之间的位置关系，又可以采用单基地和双基地两种工作模式。单基地MIMO雷达的接收端和发射端相距很近（相对于目标距离），位于同一个雷达站，而双基地MIMO雷达的接收端和发射端分别位于相距较远的两个雷达站，与目标所成的角度不同。

信号的波达方向（DOA）估计是阵列信号处理中的一个重要研究内容，它在雷达、无线通信、电磁场、声纳、医学成像和地震勘探等诸多领域的应用日益受到人们的重视。DOA估计的主要目的是确定同时处在空间某一区域内的多个感兴趣的目标源信号的空间位置。早期的DOA估计方法以常规的波束形成方法为代表。这种方法的实质是用阵元的空域采样代替时域数据，是传统时域傅里叶谱估计向空域的简单推广，因而其分辨力会受到阵列物理孔径(瑞利限制)的约束，使得位于同一个波束宽度内的空间目标无法分辨。通过增大阵列孔径可以提高分辨力，但是在实际情况中，阵列孔径的增大要受到各种因素的制约。因此，高分辨DOA估计技术应运而生，并在随后的几十年里得到了飞速的发展。Burg提出的最大熵法和Capon提出的最小方差法克服了常规波束形成算法分辨力较低的缺点，但是仍然难以突破阵列孔径的瑞利限制。Schimdt提出的多重信号分类（MUSIC）算法则是通过对阵列协方差矩阵进行特征值分解，将观测数据空间划分为相互正交的信号子空间和噪声子空间，然后利用阵列导向矢量与噪声子空间之间的正交性来获取信号的DOA估计，从而突破了以往谱估计算法中阵列孔径的瑞利限制。但是MUSIC算法需要进行谱峰穷尽搜索，运算量和存储量极大。Paulraj、Roy和Kailath基于阵列结构信号子空间的旋转不变特性提出了著名的ESPRIT算法。ESPRIT算法不需要进行谱峰搜索就可以直接给出超分辨的DOA估计，因而成为了最具应用价值的空间谱估计方法之一。这些经典的阵列高分辨DOA估计算法都在MIMO雷达中得到了应用，即MIMO-Capon算法，MIMO-MUSIC算法，接收ESPRIT算法和MIMO-ESPRIT算法等。但是在数据维数剧增的MIMO雷达系统中，这些方法不但计算量庞大，而且要得到好的性能需要极大的样本数，不利于工程实现。

发明内容

本发明的目的是：针对现有的阵列DOA估计方法应用于MIMO雷达时存在的诸多不足，如需要对数据协方差矩阵进行求逆或特征值分解运算而使得算法所涉及的计算量巨大，要获得好的性能需要极大的样本数等，本发明提出了一种迭代最小二乘（I-LS）方法对单基地MIMO雷达目标进行波达方向（DOA）估计。该方法大大地降低了所需的计算量，并提高了小样本条件下的估计性能。

本发明的技术方案概括为：首先对多个雷达发射脉冲的回波数据矩阵以及接收和发射阵列响应矩阵进行降维处理，从而使计算量减小，迭代速度加快；然后，在最小二乘条件下建立代价函数，并给出一种基于梯度下降的迭代方法求解代价函数，其中，本发明的迭代方法是交替地估计目标相对于接收、发射阵列的响应矩阵以及一组包含回波幅度信息的对角矩阵；最后，利用已知的收、发阵列流形估计目标方向。具体实现过程如下：

（1）对回波数据矩阵X(k)＝A(θ)Λ(k)B^H(θ)+W(k)进行降维处理，使其成为L×L维方阵Y(k)，用计算式表示为：

则N×L维的接收阵列响应矩阵A(θ)和M×L维的发射阵列响应矩阵B(θ)分别降维成L×L维方阵

其中，M为雷达发射通道数，N为雷达接收通道数，L为目标个数，k＝1,…,K，K表示在一个相干处理间隔内各个发射通道发射的脉冲数，U和V分别表示N×L维左降维矩阵和M×L维右降维矩阵，令β(k)＝[β₁(k),…,β_L(k)]^T，β(k)表示目标随机幅度向量，对角矩阵Λ(k)＝diag{β₁(k),…,β_L(k)}，W(k)为匹配滤波器的噪声输出，

上标^T和上标^H分别表示对矩阵或向量求转置和求共轭转置；

（2）寻找接收和发射阵列响应矩阵

以及一组对角矩阵

在最小二乘条件下建立出代价函数为 $f(\tilde{A};\mathrm{\Λ}\left(1\right),...,\mathrm{\Λ}\left(K\right);\tilde{B})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|Y\left(k\right)-\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\left|\right|}_F^2,$ 约束条件为

其中，

分别表示矩阵

的第l列，符号||·||_F和||·||分别表示矩阵或向量的F-范数和2-范数；

（3）将步骤（2）建立的代价函数

中的待定参量分为三组：基于梯度下降原理，运用迭代方法，通过交替估计

和

搜索代价函数

的最小值点，从而得到

的估计值

（4）利用步骤（3）求解出的

的估计值

恢复出目标的接收导向矢量矩阵和发射导向矢量矩阵

由于

从而得到L个目标的接收和发射导向矢量的估计值

，进而根据已知的阵列流形信息估计出各个目标的方位角θ₁,…,θ_L，即完成了目标的DOA估计。

本发明与现有技术相比具有以下特点：

1、传统的MIMO-Capon算法和MIMO-MUSIC算法都需要在空域进行全域搜索，计算量比较大且估计精度与搜索步长有关。而本发明方法可以直接得到目标的DOA估计，不需要进行谱峰搜索。

2、本发明方法避免了高维数据协方差矩阵的估计、求逆以及特征值分解运算，显著降低了计算量并提高了小样本条件下的估计性能，并且有较快的收敛速度。若单基地MIMO雷达的发射通道数和接收通道数分别为M和N，各个发射通道同时辐射由K个脉冲组成的脉冲串波形。传统的MIMO-MUSIC算法，接收ESPRIT算法和MIMO-ESPRIT算法的计算量主要来自于估计NM×NM维的观测数据协方差矩阵R并对其进行特征值分解，所涉及的计算复杂度约为O(KN²M²+N³M³)，其中KN²M²是估计R的计算量，N³M³是进行特征值分解的计算量，运算量比较大。而本发明使用迭代最小二乘（I-LS）方法，每个迭代周期的计算复杂度约为O(3KL³)。假设本发明方法经过I步迭代收敛，则整个的计算复杂度约为O(KM²N+KN²M+3IKL³)，有效地降低了计算量。图3为本发明I-LS方法的均方根误差（RMSE）随迭代次数变化的曲线图，可以看出经过几次迭代本方法即可收敛。图6表示接收ESPRIT算法，本发明I-LS方法和MIMO-ESPRIT算法的平均计算时间随样本数K变化的曲线图，从图中可以看出，I-LS方法的计算量明显小于MIMO-ESPRIT算法，且随着样本数的增加其优势更为明显。

3、本发明方法采用降维处理有效地抑制了噪声，提高了低信噪比下的估计精度。图4的（a）和（b）表示RMSE随信噪比变化的曲线图，图5的（a）和（b）表示RMSE随样本数变化的曲线图。从图4的（a）、（b）和图5的（a）、（b）中可以看出，接收ESPRIT算法仅利用了接收阵列孔径，因此其性能最差；本发明方法的估计精度优于MIMO-ESPRIT算法，这是因为：MIMO-ESPRIT算法利用两个子阵的平移不变结构对目标进行DOA估计，无法避免阵列孔径损失，而本发明方法的降维处理实质上是对接收信号进行了预滤波，有效地抑制了噪声，提高了信噪比；MIMO-ESPRIT算法需要估计MN×MN维的接收数据相关矩阵并对其进行特征值分解，这就需要大量独立同分布的样本来支撑，而I-LS方法只需要估计两个低维的导向矢量矩阵，因而在低信噪比和有限快拍条件下，其估计精度就要优于MIMO-ESPRIT算法。

附图说明

图1是单基地MIMO雷达系统结构图；

图2是本发明迭代最小二乘（I-LS）方法的流程图；

图3是试验一中本发明I-LS方法的均方根误差（RMSE）随迭代次数变化的曲线图；

图4是试验一中各方法的RMSE随信噪比（SNR）变化的曲线图：图4的（a）是样本数K＝200，方位角θ₁＝0°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随信噪比变化的曲线图；图4的（b）是样本数K＝200，方位角θ₁＝10°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随信噪比变化的曲线图；

图5是试验一中各方法的RMSE随样本数变化的曲线图：图5的（a）是SNR＝0dB，θ₁＝0°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随样本数K变化的曲线图；图5的（b）是SNR＝0dB，θ₁＝10°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随样本数K变化的曲线图；

图6是试验一中接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的平均计算时间随样本数K变化的曲线图；

图7是试验二中各方法的RMSE随信噪比（SNR）变化的曲线图：图7的（a）是样本数K＝200，方位角θ₁＝0°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随信噪比变化的曲线图；图7的（b）是样本数K＝200，方位角θ₁＝3°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随信噪比变化的曲线图；

图8是试验二中各方法的RMSE随样本数变化的曲线图：图8的（a）是SNR＝0dB，θ₁＝0°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随样本数K变化的曲线图；图8的（b）是SNR＝0dB，θ₁＝3°时接收ESPRIT算法、I-LS方法和MIMO-ESPRIT算法的RMSE随样本数K变化的曲线图。

具体实施方式

下面参照附图说明本发明的方法实施过程。

为了更好地理解本发明，先介绍单基地MIMO雷达信号模型。图1为单基地MIMO雷达系统结构图，发射通道数和接收通道数分别为M和N。远场目标源相对于发射和接收阵元具有相同的角度θ。在一个相干处理间隔内，各个发射通道同时辐射由K个脉冲组成的脉冲串波形，且M个发射波形是相互正交的。假设第m(m＝1,2,…,M)个发射波形的基带信号为s_m(t)，发射信号矢量s(t)＝[s₁(t),s₂(t),…,s_M(t)]^T满足：

${\∫}_{T}s\left(t\right)s^H\left(t\right)=I_M---\left(1\right)$

其中，T表示发射脉冲的宽度，I_M表示M阶单位矩阵。对于方位角为θ的远场目标源，到达目标处的雷达信号可以表示为如下形式：

b^T(θ)s(t)   (2)

其中，b(θ)为发射导向矢量，即目标相对于发射阵列的导向矢量。若发射阵列是均匀线阵且阵元间距为d_T，则：

$b\left(\mathrm{\θ}\right)={[1,e^{-j2\mathrm{\π}(d_T/\mathrm{\λ})\sin \mathrm{\θ}},...,e^{-j2\mathrm{\π}(d_T/\mathrm{\λ})(M-1)\sin \mathrm{\θ}}]}^T---\left(3\right)$

因此，在第k个脉冲周期，经过单个目标反射之后到达接收阵列的回波信号可以表示为：

$x({\mathrm{kT}}_r+t+\frac{2r}{c})=\mathrm{\β}\left(k\right)a\left(\mathrm{\θ}\right)b^T\left(\mathrm{\θ}\right)s\left(t\right)+w_k\left(t\right),k=1,...,K---\left(4\right)$

其中，r表示雷达与目标之间的距离，T_r为脉冲重复周期，c表示光速；

，f_d代表由目标径向速度引起的多普勒频率，ρ表示目标回波复幅度；w_k(t)为第k个脉冲周期的接收机噪声，它服从均值为零，协方差矩阵为

的复高斯分布；a(θ)为接收导向矢量，即目标相对于接收阵列的导向矢量。若接收阵列是均匀线阵且阵元间距为d_R，则有：

$a\left(\mathrm{\θ}\right)={[1,e^{-j2\mathrm{\π}(d_R/\mathrm{\λ})\sin \mathrm{\θ}},...,e^{-j2\mathrm{\π}(d_R/\mathrm{\λ})(N-1)\sin \mathrm{\θ}}]}^T---\left(5\right)$

式(4)给出的是单个目标的回波信号模型。若观测距离单元存在L个目标，则式(4)可以扩展为如下形式：

$x({\mathrm{kT}}_r+t+\frac{2r}{c})=\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(k\right)a\left({\mathrm{\θ}}_l\right)b^T\left({\mathrm{\θ}}_l\right)s\left(t\right)+w_k\left(t\right)---\left(6\right)$

式(6)给出的是多个目标的回波信号模型。在各个接收通道，分别利用M个参考发射信号对K个脉冲的回波数据进行匹配滤波，那么所有的匹配滤波输出可表示为如下N×M维的矩阵：

$\begin{array}{c}X\left(k\right)={\∫}_{T}x({\mathrm{kT}}_r+t+\frac{2r}{c})s^H\left(t\right)\mathrm{dt}\\ =\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(k\right)a\left({\mathrm{\θ}}_l\right)b^T\left({\mathrm{\θ}}_l\right)+W\left(k\right)\end{array}---\left(7\right)$

式(7)中，

是匹配滤波器的噪声输出。如果将矩阵X(k)进行矩阵的向量化运算，即将X(k)的列向量依次堆砌，则可以得到MIMO雷达的虚拟阵列快拍数据为：

$\begin{array}{c}\tilde{x}\left(k\right)=\mathrm{vec}\left(X\left(k\right)\right)\\ =\underset{l=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_l\left(k\right)b\left({\mathrm{\θ}}_l\right)\⊗a\left({\mathrm{\θ}}_l\right)+\tilde{w}\left(k\right)\\ =C\left(\mathrm{\θ}\right)\mathrm{\β}\left(k\right)+\widetilde{w\left(k\right)}\end{array}---\left(8\right)$

上式中，C(θ)＝[b(θ₁)

a(θ₁),…,b(θ_L)

a(θ_L)]表示NM×L维的虚拟阵列响应矩阵，β(k)＝[β₁(k),…,β_L(k)]^T是目标随机幅度矢量，

是零均值的复高斯白噪声矢量，符号vec表示矩阵的向量化运算，符号表示克罗内克积。

由式(8)可知，虚拟阵列协方差矩阵的最大似然估计可以表示为：

$\hat{R}=\frac{1}{K}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\tilde{x}\left(k\right){\tilde{x}}^H\left(k\right)---\left(9\right)$

式(7)的信号模型可以等价地写成如下的矩阵形式：

X(k)＝A(θ)Λ(k)B^H(θ)+W(k)   (10)

上式中，A(θ)＝[a(θ₁),…,a(θ_L)]表示N×L维的接收阵列响应矩阵，B(θ)＝[b^*(θ₁),…,b^*(θ_L)]为M×L维的发射阵列响应矩阵，L×L维的对角矩阵Λ(k)＝diag{β₁(k),…,β_L(k)}。由式(10)可以看出，单基地MIMO雷达的接收数据具有天然的对角结构，在此基础上，本发明提出一种迭代最小二乘（I-LS）方法来交替地估计两个阵列响应矩阵A(θ)和B(θ)，从而得到L个目标的发射和接收导向矢量，进而根据已知的阵列流形直接得到目标的DOA估计。

根据图2，本发明的迭代最小二乘（I-LS）方法步骤如下：

步骤1：对回波数据矩阵X(k)、接收阵列响应矩阵A(θ)和发射阵列响应矩阵B(θ)进行降维处理。实际情况中，有L＜＜M,L＜＜N。为降低运算量和处理的方便，先将式(10)中的数据矩阵X(k)降维成L×L维方阵Y(k)，则N×L维矩阵A(θ)和M×L维矩阵B(θ)也分别降维成L×L维的方阵

降维过程用如下的式子表示为：

$Y\left(k\right)=U^{H}X\left(k\right)V=\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H+\tilde{W}\left(k\right)---\left(11\right)$

其中，U为N×L维左降维矩阵，V为M×L维右降维矩阵，

实际上，理想的降维矩阵U和V应满足如下条件：

$\left\{\begin{array}{c}\mathrm{span}\left(U\right)=\mathrm{span}\left(A\left(\mathrm{\θ}\right)\right)\\ \mathrm{span}\left(V\right)=\mathrm{span}\left(B\left(\mathrm{\θ}\right)\right)\end{array}\right.---\left(12\right)$

从而降维后的阵列响应矩阵

均为满秩方阵。式中span(U)和span(V)分别表示矩阵U和V的列空间。为处理上的方便，一般将矩阵U和V的列向量取为单位正交。无噪声情况下，可以将式(10)中某个数据矩阵X(k)的奇异值分解所得的信号子空间作为降维矩阵。但实际中，噪声总是存在的，这时可以利用多个数据矩阵

来获得高质量的降维矩阵。

令矩阵

$\left\{\begin{array}{c}\mathrm{\Φ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}X\left(k\right)X^H\left(k\right)\\ \mathrm{\Ψ}=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{X}^H\left(k\right)X\left(k\right)\end{array}\right.---\left(13\right)$

Φ和Ψ的特征值分解分别为：

$\mathrm{\Φ}=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{u}_i{u}_i^H,\mathrm{\Ψ}=\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{\μ}}_j{v}_j{v}_j^H,---\left(14\right)$

上式中，λ_i和μ_j分别是矩阵Φ和Ψ的特征值，u_i和v_j分别是相应的矩阵Φ和Ψ关于特征值λ_i和μ_j的特征向量，其中，i＝1,…,N，j＝1,…,M。分别取矩阵Φ和Ψ的L个主特征向量构成左、右降维矩阵，即：

U＝[u₁,u₂,…,u_L]，V＝[v₁,v₂,L,v_L]  (15)

由于未受到噪声污染的数据矩阵具有低秩结构，因此式(11)的降维处理不仅使得信号分量完全保留，还起到了抑制噪声的作用。

步骤2：寻找接收和发射阵列响应矩阵

及一组对角矩阵

在最小二乘条件下建立代价函数。利用式(11)中数据矩阵的对角结构，寻找接收和发射阵列响应矩阵

及一组对角矩阵

在最小二乘条件下最小化代价函数：

$f(\tilde{A};\mathrm{\Λ}\left(1\right),...,\mathrm{\Λ}\left(K\right);\tilde{B})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|Y\left(k\right)-\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\left|\right|}_F^2---\left(16\right)$

但是，式(16)的代价函数存在一个尺度模糊问题。该问题可以大致描述为：若存在L×L维矩阵分别满足

（F和G均为L×L维可逆的对角矩阵），则最小化代价函数(16)本质上等价于最小化代价函数：

$f(\underset{~}{A};\underset{~}{\mathrm{\Λ}}\left(1\right),...,\underset{~}{\mathrm{\Λ}}\left(K\right);\underset{~}{B})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|Y\left(k\right)-\underset{~}{A}\underset{~}{\mathrm{\Λ}}\left(k\right){\underset{~}{B}}^H\left|\right|}_F^2---\left(17\right)$

其中，为消除尺度模糊问题，不失一般性地，可约束矩阵

各列的模为1，进而将式(16)修正为如下形式：

$\left\{\begin{array}{c}f(\tilde{A};\mathrm{\Λ}\left(1\right),...,\mathrm{\Λ}\left(K\right);\tilde{B})=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|Y\left(k\right)-\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\left|\right|}_F^2\\ s.t.{\left|\right|{\tilde{a}}_l\left|\right|}^2={\left|\right|{\tilde{b}}_l\left|\right|}^2=1,l=1,...,L\end{array}\right.---\left(18\right)$

其中，

分别表示矩阵

的第l列。

由于目标的DOA信息包含在矩阵

各列元素的相位中，归一化只是对幅度进行操作，因此对矩阵

进行归一化并不影响对目标的DOA估计。

步骤3：基于梯度下降原理，运用迭代方法，通过交替估计

搜索步骤2建立的代价函数的最小值点。下面介绍本发明的迭代方法：

将代价函数

中的待定参量分为三组：若分别固定其中的任意两组参量，那么代价函数就简化为关于第三组参量的二次函数，相应的求解代价函数关于这组参数的最小值点就成为一个最小二乘问题。根据这一思想，本发明基于梯度下降原理提出一种迭代方法，通过交替估计

来寻找代价函数的最小值点。

（1）估计对角矩阵组

由于式(16)等号右边的K个分量均大于等于零，如果固定矩阵

那么最小化

就等价于分别最小化如下的K个子函数：

$f_k\left(\mathrm{\Λ}\left(k\right)\right)={\left|\right|Y\left(k\right)-\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\left|\right|}_F^2---\left(19\right)$

将

和Λ(k)＝diag{β₁(k),…,β_L(k)}分别代入式(19)，有：

$\begin{array}{c}{f}_k\left(\mathrm{\Λ}\left(k\right)\right)={\left|\right|Y\left(k\right)-\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_i\left(k\right){\tilde{a}}_i{\tilde{b}}_i^H\left|\right|}_F^2\\ =\mathrm{tr}\{\[Y\left(k\right)-\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_i\left(k\right){\tilde{a}}_i{\tilde{b}}_i^H\]{\[Y\left(k\right)-\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_j\left(k\right){\tilde{a}}_j{\tilde{b}}_j^H\]}^H\}\\ =\mathrm{tr}\left\{Y\left(k\right)Y^H\left(k\right)\right\}-\mathrm{tr}\left\{Y\left(k\right)\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_j^{*}\left(k\right){\tilde{b}}_j{\tilde{a}}_j^H\right\}-\\ \mathrm{tr}\left\{\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_i\left(k\right){\tilde{a}}_i{\tilde{b}}_i^H{Y}^H\left(k\right)\right\}+\mathrm{tr}\{\[\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_i\left(k\right){\tilde{a}}_i{\tilde{b}}_i^H\]\[\underset{j=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_j^{*}\left(k\right){\tilde{b}}_j{\tilde{a}}_j^H\]\}\end{array}---\left(20\right)$

令式(20)关于β_l(k)的共轭导数等于零，可以得到：

${\tilde{a}}_l^{H}Y\left(k\right){\tilde{b}}_l={\tilde{a}}_l^H\[\underset{i=1}{\overset{L}{\mathrm{\Σ}}}{\mathrm{\β}}_i\left(k\right){\tilde{a}}_i{\tilde{b}}_i^H\]{\tilde{b}}_l---\left(21\right)$

将l分别取为1,…,L，可以得到如下矩阵组：

$\left\{\begin{array}{c}{\tilde{a}}_1^{H}Y\left(k\right){\tilde{b}}_1={\mathrm{\β}}_1\left(k\right){\tilde{a}}_1^H{\tilde{a}}_1{\tilde{b}}_1^H{\tilde{b}}_1+{\mathrm{\β}}_2\left(k\right){\tilde{a}}_1^H{\tilde{a}}_2{\tilde{b}}_2^H{\tilde{b}}_1+...+{\mathrm{\β}}_L\left(k\right){\tilde{a}}_1^H{\tilde{a}}_L{\tilde{b}}_L^H{\tilde{b}}_1\\ {\tilde{a}}_2^{H}Y\left(k\right){\tilde{b}}_2={\mathrm{\β}}_1\left(k\right){\tilde{a}}_2^H{\tilde{a}}_1{\tilde{b}}_1^H{\tilde{b}}_2+{\mathrm{\β}}_2\left(k\right){\tilde{a}}_2^H{\tilde{a}}_2{\tilde{b}}_2^H{\tilde{b}}_2+...+{\mathrm{\β}}_L\left(k\right){\tilde{a}}_2^H{\tilde{a}}_L{\tilde{b}}_L^H{\tilde{b}}_2\\ \·\\ \·\\ \·\\ {\tilde{a}}_L^{H}Y\left(k\right){\tilde{b}}_L={\mathrm{\β}}_1\left(k\right){\tilde{a}}_L^H{\tilde{a}}_1{\tilde{b}}_1^H{\tilde{b}}_L+{\mathrm{\β}}_2\left(k\right){\tilde{a}}_L^H{\tilde{a}}_2{\tilde{b}}_2^H{\tilde{b}}_L+...+{\mathrm{\β}}_L\left(k\right){\tilde{a}}_L^H{\tilde{a}}_L{\tilde{b}}_L^H{\tilde{b}}_L\end{array}\right.---\left(22\right)$

将式(22)写成如下矩阵形式：

d(k)＝Eβ(k)   (23)

其中，向量 $d\left(k\right)={\[{\tilde{a}}_1^{H}Y\left(k\right){\tilde{b}}_1,{\tilde{a}}_2^{H}Y\left(k\right){\tilde{b}}_2,...,{\tilde{a}}_L^{H}Y\left(k\right){\tilde{b}}_L\]}^T,$ β(k)＝[β₁(k),…,β_L(k)]^T，矩阵

其中，符号⊙表示Hadamard积运算，即将矩阵的对应元素相乘。固定矩阵

由式(23)可以求得对角矩阵Λ(k)为：

Λ(k)＝diag{β(k)}＝diag{E^-1d(k)}  (25)

（2）估计接收阵列响应矩阵

固定对角矩阵组

和矩阵

代价函数

对于矩阵

求最小值。那么关于的二次代价函数可以表示为如下形式：

$\begin{array}{c}f\left(\tilde{A}\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|Y\left(k\right)-\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\left|\right|}_F^2\\ =\mathrm{tr}\left\{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}Y\left(k\right)Y^H\left(k\right)\right\}-\mathrm{tr}\{\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}Y\left(k\right)\widetilde{B{\mathrm{\Λ}}^H\left(k\right)\]{\stackrel{~}{A}}^H\}-}\\ \mathrm{tr}\{\tilde{A}\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\mathrm{\Λ}\left(k\right){\tilde{B}}^H{Y}^H\left(k\right)\]\}+\mathrm{tr}\{\tilde{A}\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\tilde{B}{\mathrm{\Λ}}^H\left(k\right)\]{\tilde{A}}^H\}\end{array}---\left(26\right)$

将式(26)关于

求偏导并令偏导数等于零，可以求得：

$\tilde{A}=\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}Y\left(k\right)\tilde{B}{\mathrm{\Λ}}^H\left(k\right)\]{\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\tilde{B}{\mathrm{\Λ}}^H\left(k\right)\]}^{-1}---\left(27\right)$

（3）估计发射阵列响应矩阵

固定对角矩阵组和矩阵

代价函数对于矩阵

求最小值。那么关于矩阵

的二次代价函数可以表示为如下形式：

$\begin{array}{c}f\left(\tilde{B}\right)=\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\left|\right|Y\left(k\right)-\tilde{A}\mathrm{\Λ}\left(k\right){\tilde{B}}^H\left|\right|}_F^2\\ =\mathrm{tr}\left\{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{Y}^H\left(k\right)Y\left(k\right)\right\}-\mathrm{tr}\{\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{Y}^H\left(k\right)\tilde{A}\mathrm{\Λ}\left(k\right)\]{\tilde{B}}^H\}-\\ \mathrm{tr}\{\tilde{B}\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\Λ}}^H\left(k\right){\tilde{A}}^{H}Y\left(k\right)\]\}+\mathrm{tr}\{\tilde{B}\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\Λ}}^H\left(k\right){\tilde{A}}^H\tilde{A}\mathrm{\Λ}\left(k\right)\]{\tilde{B}}^H\}\end{array}---\left(28\right)$

将式(28)关于

求偏导并令偏导数等于零，可以求得：

$\tilde{B}=\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{Y}^H\left(k\right)\tilde{A}\mathrm{\Λ}\left(k\right)\]{\[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\Λ}}^H\left(k\right){\tilde{A}}^H\tilde{A}\mathrm{\Λ}\left(k\right)\]}^{-1}---\left(29\right)$

根据上述分析，给出本发明I-LS方法的具体流程：

首先，给定初值矩阵

当迭代步数i＝1,2,…时，迭代过程如下：

1）将

分别代入式(25)中，更新对角矩阵组

${\mathrm{\Λ}}_{\left(i\right)}\left(k\right)=\mathrm{diag}\left\{E_{(i-1)}^{-1}{d}_{(i-1)}\left(k\right)\right\},$

其中，向量 $d_{(i-1)}\left(k\right)={[{\tilde{a}}_{1,(i-1)}^{H}Y\left(k\right){\tilde{b}}_{1,(i-1)},{\tilde{a}}_{2,(i-1)}^{H}Y\left(k\right){\tilde{b}}_{2,(i-1)},...,{\tilde{a}}_{L,(i-1)}^{H}Y\left(k\right){\tilde{b}}_{L,(i-1)}]}^T$ ，矩阵

2）将

分别代入式(27)中，更新矩阵

${\tilde{A}}_{\left(i\right)}=\left[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}Y\left(k\right){\tilde{B}}_{(i-1)}{\mathrm{\Λ}}_{\left(i\right)}^H\left(k\right)\right]{\left[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\Λ}}_{\left(i\right)}\left(k\right){\tilde{B}}_{(i-1)}^H{\tilde{B}}_{(i-1)}{\mathrm{\Λ}}_{\left(i\right)}^H\left(k\right)\right]}^{-1},$

同时将中的各列归一化；

3）将

分别代入式(29)中，更新

${\tilde{B}}_{\left(i\right)}=\left[\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{Y}^H\left(k\right){\tilde{A}}_{\left(i\right)}{\mathrm{\Λ}}_{\left(i\right)}\left(k\right)\right][\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\Λ}}_{\left(i\right)}^H\left(k\right){\tilde{A}}_{\left(i\right)}^H{\tilde{A}}_{\left(i\right)}{\mathrm{\Λ}}_{\left(i\right)}\left(k\right){]}^{-1},$

并将

中的各列归一化；

4）将 $e=|f({\tilde{A}}_{\left(i\right)};{\left\{{\mathrm{\Λ}}_{\left(i\right)}\left(k\right)\right\}}_{k=1}^K;{\tilde{B}}_{\left(i\right)})-f({\tilde{A}}_{(i-1)};{\left\{{\mathrm{\Λ}}_{(i-1)}\left(k\right)\right\}}_{k=1}^K;{\tilde{B}}_{(i-1)})|$ 与门限ε(0＜ε＜＜1)比较大小。若e＜ε，则迭代停止，并令

的估计值分别为

否则，将

作为初值返回步骤1）继续迭代。其中，e表示相邻两次迭代的代价函数值之间的差值。

步骤4：运用

的估计值

恢复出目标导向矢量矩阵

从而得到L个目标的接收和发射导向矢量的估计值

进而根据已知的阵列流形信息估计目标的波达方向θ1,…,θL。当得到

的估计值

时，可根据如下式子恢复目标导向矢量矩阵A(θ)和B(θ)：

$\left\{\begin{array}{c}\hat{A}\left(\mathrm{\θ}\right)=U\hat{\tilde{A}}\\ \hat{B}\left(\mathrm{\θ}\right)=V\hat{\tilde{B}}\end{array}\right.---\left(30\right)$

由此，可以得到L个目标的接收导向矢量的估计值

和发射导向矢量的估计值

然后，根据已知的阵列流形信息估计目标的波达方向θ₁,…,θ_L，即完成了对目标的DOA估计。

仿真试验对比：

为了进一步说明本发明的迭代最小二乘（I-LS）方法较传统单基地MIMO雷达DOA估计方法（如MIMO-ESPRIT算法和接收ESPRIT算法）的优越性，做如下两个仿真试验。

系统模型：单基地MIMO雷达的发射和接收阵列均为均匀线阵且阵元数为M＝N＝10。载波波长λ＝0.3m，脉冲重复周期T_r＝5×10^-4s。为防止接收端出现栅瓣，设定接收阵元间距为半波长，即d_R＝λ/2＝0.15m，发射阵元间距d_T＝Nd_R＝1.5m。

试验一：同一距离分辨单元内存在两个雷达目标，方位角分别为θ₁＝0°和θ₂＝10°。

图3是样本数K＝200，目标信噪比SNR=0dB时I-LS方法的均方根误差(RMSE)随迭代次数变化的曲线图。从图3可以看出，I-LS方法的RMSE经过5～7次迭代就已经接近最小值，说明I-LS方法有较快的收敛速度。

图4的（a）和（b）是样本数K＝200时三种方法的RMSE随SNR变化的曲线图，图5的（a）和（b）为SNR＝0dB时三种方法的RMSE随样本数K变化的曲线图。上述四个图中还同时给出了估计值RMSE的克拉-美罗界（CRB）。从图4的（a）、图4的（b）和图5的（a）、图5的（b）中可以看出，接收ESPRIT算法性能最差（这是由于该算法仅利用了接收阵列孔径），而本发明的I-LS方法的估计精度优于MIMO-ESPRIT算法。图6给出了三种方法的平均计算时间随样本数K变化的曲线图。从图6可以看出，本发明的I-LS方法的计算量明显小于MIMO-ESPRIT算法，并且随着样本数的增加其优势更为明显。

试验二：两目标角度间隔较小，目标2的角度变为θ₂＝3°，其余实验参数同试验一。

图7的（a）和图7的（b）给出了三种方法的RMSE随SNR的变化曲线，图8的（a）和图8的（b）给出了三种方法的RMSE随样本数K的变化曲线。在目标角度间隔较小的情况下，与试验一的结果相比，三种方法的估计精度都有一定的下降，尤其是接收ESPRIT算法，其RMSE增大了近5个dB。但是本发明的I-LS方法的估计精度仍明显优于其他两种算法。

Claims (5)

1.基于迭代最小二乘方法的MIMO雷达DOA估计方法，其特征是：用迭代最小二乘方法求解经过降维处理的接收和发射阵列响应矩阵；首先对多个雷达发射脉冲的回波数据矩阵以及接收和发射阵列响应矩阵进行降维处理，然后在最小二乘条件下建立代价函数，并利用基于梯度下降的迭代方法求解代价函数，最后利用已知的收、发阵列流形估计目标方向；具体实现过程如下： 

1）对回波数据矩阵X(k)＝A(θ)Λ(k)B^H(θ)+W(k)进行降维处理，使其成为L×L维方阵Y(k)，用计算式表示为：

则N×L维的接收阵列响应矩阵A(θ)和M×L维的发射阵列响应矩阵B(θ)分别降维成L×L维方阵 

和

其中，M为雷达发射通道数，N为雷达接收通道数，L为目标个数，K表示在一个相干处理间隔内各个发射通道发射的脉冲数，k＝1,…,K，U和V分别表示N×L维左降维矩阵和M×L维右降维矩阵，令β(k)＝[β₁(k),…,β_L(k)]^T，β(k)表示目标随机幅度向量，对角矩阵Λ(k)＝diag{β₁(k),…,β_L(k)}，W(k)为匹配滤波器的噪声输出，

上标^T和上标^H分别表示对矩阵或向量求转置和求共轭转置； 

2）寻找接收和发射阵列响应矩阵

以及一组对角矩阵

在最小二乘条件下建立出代价函数为

约束条件为 

其中，

分别表示矩阵

的第l列，符号||·||_F和||·||分别表示矩阵或向量的F-范数和2-范数； 

3）将步骤2）建立的代价函数

中的待定参量分为三组： 

和基于梯度下降原理，运用迭代方法，通过交替估计

和 

搜索代价函数

的最小值点，从而得到的估计值

4）运用步骤3）所得的

的估计值

估计目标的方位角θ₁,…,θ_L，即完成了对目标的DOA估计。 

2.根据权利要求1所述的基于迭代最小二乘方法的MIMO雷达DOA估计方法，其特征是：对回波数据矩阵X(k)、接收阵列响应矩阵A(θ)和发射阵列响应矩阵B(θ)进行降维处理，具体过程如下： 

{1}降维处理方法用如下的式子表示为： 

理想的降维矩阵U和V满足如下条件： 

降维后的阵列响应矩阵

均为满秩方阵，式中span(U)和span(V)分别表示矩阵U和V的列空间，为处理上的方便，将矩阵U和V的列向量取为单位正交； 

{2}无噪声情况下，将K个数据矩阵中某个数据矩阵X(k)的奇异值分解所得的信号子空间作为降维矩阵，但实际中，噪声总是存在的，则利用多个数据矩阵

来获得高质量的降维矩阵，令矩阵 

Φ和Ψ的特征值分解分别为： 

式（4）中，λ_i和μ_j分别是矩阵Φ和Ψ的特征值，u_i和v_j分别是相应的矩阵Φ和Ψ关于特征值λ_i和μ_j的特征向量，其中，i＝1,…,N，j＝1,…,M； 

{3}分别取矩阵Φ和Ψ的L个主特征向量构成左、右降维矩阵，即： 

U＝[u₁,u₂,…,u_L]，V＝[v₁,v₂,…,v_L]   (5) 。

3.根据权利要求2所述的基于迭代最小二乘方法的MIMO雷达DOA估计方法，其特征是：用最小二乘方法求解经过降维处理的接收和发射阵列响应矩阵及其迭代求解的过程如下： 

{1}最小二乘条件下的代价函数及约束条件为： 

将代价函数中的待定参量分为三组：

本发明的迭代方法是基于梯度下降原理，交替估计这三组待定参量

搜索代价函数

的最小值点； 

{2}估计对角矩阵组

由于式(6)等号右边的K个分量均大于等于零，固定矩阵最小化 

就等价于分别最小化如下的K个子函数： 

将

，

和Λ(k)＝diag{β₁(k),…,β_L(k)}分别代入式(7)，得： 

令式(8)关于β_l(k)的共轭导数等于零，得： 

将l分别取为1,…,L，得到如下矩阵组： 

将式(10)写成如下的矩阵形式： 

d(k)＝Eβ(k)， 0(11) 

其中，向量

β(k)＝[β₁(k),…,β_L(k)]^T，矩阵 

其中，符号⊙表示Hadamard积运算，即将矩阵的对应元素相乘，固定矩阵

由式(11)求得对角矩阵Λ(k)为： 

Λ(k)＝diag{β(k)}＝diag{E^-1d(k)}，  (13) 

{3}估计接收阵列响应矩阵

固定对角矩阵组

和矩阵

代价函数对于矩阵

求最小值，将关于的二次代价函数表示为如下形式： 

令式(14)关于

的偏导数等于零，得到： 

求解式(15)，得到： 

{4}估计发射阵列响应矩阵

固定对角矩阵组和矩阵

代价函数

对于矩阵

求最小值，将关于矩阵

的二次代价函数表示为如下形式： 

令式(17)关于的偏导数等于零，得到： 

求解式(18)，得到： 

4.根据权利要求3所述的基于迭代最小二乘方法的MIMO雷达DOA估计方法，其特征是：用迭代方法求解接收和发射阵列响应矩阵的流程如下： 

[1]给定初值矩阵

[2]将

分别代入式(13)中，更新对角矩阵组

其中，向量

矩阵 

[3]将

分别代入式(16)中，更新矩阵

同时将中的各列归一化； 

[4]将

分别代入式(19)中，更新

并将

中的各列归一化； 

[5]将

与门限ε进行比较：若e＜ε，则迭代停止，令

的估计值分别为

否则，将

作为初值返回步骤[2]继续迭代，其中，e表示相邻两次迭代的代价函数值之间的差值，0＜ε＜＜1。 

5.根据权利要求4所述的基于迭代最小二乘方法的MIMO雷达DOA估计方法，其特征是：利用求解出的的估计值恢复出目标的接收和发射导向矢量矩阵由于

从而得到L个目标的接收导向矢量的估计值

和发射导向矢量的估计值 

进而根据已知的阵列流形信息估计出目标的方位角θ₁,…,θ_L，即完成了对目标的DOA估计。 

Images