CN104123462A - 实多项式求根实现均匀线阵的谱music方法 - Google Patents

Abstract

本发明属于雷达信号处理技术领域，特别涉及实多项式求根实现均匀线阵的谱MUSIC方法。该实多项式求根实现均匀线性阵列的谱MUSIC方法包括以下步骤：步骤1，求复多项式系数矢量；步骤2，计算窗矢量；步骤3，计算所有区间的多项式系数；步骤4，求解MUSIC谱的极值。本发明通过用傅里叶变换快速计算复多项式系数矢量，把求MUSIC谱的极点转化为求多组低阶多项式的根，对信号子空间相关矢量做加窗傅里叶变换得到所有组的多项式系数，克服了传统方法需要精细的角度搜索算法、计算复杂度高和数值稳定性差的不足，具有计算复杂度低、数值稳定性好的优点，在低复杂度估计波达方向方面有巨大的潜力。

Description

实多项式求根实现均匀线阵的谱MUSIC方法

技术领域

本发明属于雷达信号处理技术领域，特别涉及实多项式求根实现均匀线阵的谱MUSIC方法，更进一步涉及阵列信号处理中用谱MUSIC估计均匀线性阵列波达方向的方法。本发明运用多窗技术，可以有效地实现谱MUSIC。运用多窗傅里叶变换，把高阶复多项式分解为多组低阶实多项式，使得所有的根可以同时求解，并且具有低的计算复杂度和好的数值稳定性。本发明在低复杂度估计波达方向方面有巨大的潜力。

背景技术

波达方向估计是阵列信号处理的一项重要内容。在已提出的算法中，多重信号分类方法(MUSIC)，因为其优良的波达方向估计性能而被广泛应用，但是这种方法的不足是需要精细的角度搜索算法，具有较高的计算复杂度。

为了降低网格搜索的计算复杂度，B.D.Rao和K.V.S.Hari在文献“Performance analysis of root-music”(IEEE Trans.Acoust.,Speech,SignalProcess.37(12)(Dec.1989)1939–1949)中提出了求根MUSIC技术，用2N-2阶复多项式求根实现波达方向估计，其中N是阵元数。这种方法的缺点是高阶复多项式求根具有较高的计算复杂度，并且不能保证数值稳定性。

发明内容

本发明的目的在于克服现有技术的不足，提出实多项式求根实现均匀线性阵列的谱MUSIC方法。该方法弥补了现有方法需要精细的网格搜索算法、计算复杂度高和数值稳定性差的不足，充分利用在有限空间带宽里的任一导向矢量可以用低阶多项式矢量表示的性质，使MUSIC谱的极点由多组低阶多项式求根确定，因而有更低的计算复杂度和更好的数值稳定性，在低复杂度估计波达方向方面有巨大的潜力。

实现本发明的基本思路是：首先利用在有限空间带宽里的任一导向矢量可以用低阶多项式矢量表示的性质，把求MUSIC谱的极点转化为求多组低阶多项式的根。然后，通过对信号子空间相关矢量做加窗傅里叶变换得到所有组的多项式系数。最后用低阶多项式求根的方法得到MUSIC谱的极点，从而得到均匀线性阵列波达方向的估计。

为实现上述技术目的，本发明采用如下技术方案予以实现。

实多项式求根实现均匀线性阵列的谱MUSIC方法包括以下步骤：

步骤1，利用雷达接收回波数据，得出回波数据的协方差矩阵根据回波数据的协方差矩阵计算出MUSIC谱极值点对应的复多项式系数矢量b；

步骤2，根据MUSIC谱极值点对应的复多项式系数矢量b，构造MUSIC谱函数；把空间频率划分为M′段，M′为大于1的自然数，划分后的每个空间频率段的长度为ε；然后根据以下公式得出多项式适应矩阵B：

$\underset{B}{\min }\underset{g\∈[0,\mathrm{\ϵ}]}{\mathrm{\Σ}}\mathrm{\ζ}\left(g\right)$

其中，

$\mathrm{\ζ}\left(g\right)=\frac{{\left|\right|c\left(g\right)-\mathrm{Bg}\left|\right|}_2^2}{{\left|\right|c\left(g\right)\left|\right|}_2^2},g={\left[\begin{array}{cccc}{g}^L& g^{L-1}& ...& 1\end{array}\right]}^T$

c(g)＝[e^-j2π(N-1)g,e^-j2π(N-2)g,...,e^j2π(N-1)g]^T

其中，多项式适应矩阵B为(2N-1)×(L+1)维的矩阵，N为雷达的均匀线阵的阵元数，L为设定的多项式阶数，g∈[0,ε)，||·||₂表示l₂范数，上标T表示矩阵或向量的转置；

步骤3，根据以下公式计算划分后的每个空间频率段的每个多项式系数，划分后的第m个空间频率段的第l个多项式系数γ_ml的计算公式为：

其中，m取1至M′，l取1至L+1，表示Hadamard积，t_l为多项式适应矩阵B的第l列，w_m为：

w_m＝[e^{-j2π(N-1)[-M′/2+(m-1)]ε},e^{-j2π(N-2)[-M′/2+(m-1)]ε},...,e^{j2π(N-1)[-M′/2+(m-1)]ε}]^T；

步骤4，根据划分后的每个空间频率段的所有L+1个多项式系数，得出划分后的每个空间频率段的MUSIC谱函数的极点，根据划分后的每个空间频率段的MUSIC谱函数的极点，得出划分后的每个空间频率段的MUSIC谱函数的极值。

本发明的特点和进一步改进在于：

所述步骤1的具体子步骤为：

(1.1)将入射至雷达均匀线阵的信号的个数表示为P，P为大于1的自然数；则雷达第q次快拍获取的回波数据矢量x(q)为：

$x\left(q\right)=\underset{p-1}{\overset{P}{\mathrm{\Σ}}}{s}_p\left(q\right)a\left(\mathrm{\θ}\right)+n\left(q\right)$

其中，q取1至Q，Q为雷达的快拍次数；θ为入射至雷达均匀线阵的每个信号的入射角，s_p(q)是入射至雷达均匀线阵的第p个信号的复幅度，p取1至P；a(θ)表示入射至雷达均匀线阵的每个信号的导向矢量，n(q)表示高斯白噪声矢量；

(1.2)利用雷达每次快拍获取的回波数据矢量，估计出回波数据的协方差矩阵为：

$\hat{R}=\frac{1}{Q}\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}x\left(q\right)x^H\left(q\right)$

其中，上标H表示取共轭转置；

对回波数据的协方差矩阵进行特征分解，在回波数据的协方差矩阵的特征值中，选取P个最大的特征值，将选取的P个最大的特征值对应的特征矢量分别表示为u_s1至u_sP；利用选取的P个最大的特征值组成信号子空间U_s，U_s＝[u_s1,...,u_sP]，在回波数据的协方差矩阵的其余特征值中，选取P个最大的特征值对应的特征矢量组成噪声子空间U_n；

(1.3)在特征矢量u_s1至u_sP中，针对每个特征矢量作傅里叶运算，特征矢量u_sp的傅里叶运算过程表示为：

${\mathrm{\δ}}_p=F_M^{*}\left[\begin{array}{c}{u}_{\mathrm{sp}}\\ 0_{(M-N)\×1}\end{array}\right]$

其中，F_M是M×M维的傅里叶变换矩阵，*表示共轭，ceil(·)表示向上取整，u_sp表示回波数据的协方差矩阵的P个最大的特征值中第p个特征值对应的特征矢量；0_(M-N)×1表示M-N维的全零列向量；

根据以下公式得出MUSIC谱极值点对应的复多项式系数矢量b，

η＝[-jπ(N-1)，-jπ(N-2) … 0 … jπ(N-2)，jπ(N-1)]^Τ

${\stackrel{\‾}{b}}_0={\left[\begin{array}{ccccccc}0& ...& 0& N& 0& ...& 0\end{array}\right]}^T$

其中，表示Hadamard积，上标T表示矩阵或向量的转置，N为雷达的均匀线阵的阵元数，为2N-1维的列向量，的第N个元素为N，其余元素均为0；为矢量中非零元素构成的矢量，即：

${\stackrel{\‾}{b}}_1=\left[\begin{array}{c}{\stackrel{\‾}{b}}_2\\ 0_{(M-(2N-1)\×1)}\end{array}\right]$

矢量为：

${\mathrm{\Δ}}_M=\mathrm{diag}\begin{array}{cccc}\left(\right[\exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(-\frac{M}{2})\right)& \exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(-\frac{M}{2}+1)\right)& ...& \exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(\frac{M}{2}-1)\right)\left]\right)\end{array}$

其中，*表示共轭，表示M×M维的傅里叶逆变换矩阵，diag(·)表示求对角矩阵。

在步骤2中，构造的MUSIC谱函数表示为h(f)，h(f)＝b^Τc(f)，c(f)为：

c(f)＝[e^-j2π(N-1)f e^-j2π(N-2)f … 1 … e^j2π(N-1)f]^Τ

其中，f表示频率变量；

在步骤2中，划分后的第m个空间频率段中频率点f_m(g)的表达式为：

f_m(g)＝[-M′/2+(m-1)]ε+g

其中，m取1至M′，g∈[0,ε)；

然后将划分后的第m个空间频率段导向矢量定义为：

其中，上标T表示矩阵或向量的转置，表示Hadamard积，并且有：

c(g)＝[e^-j2π(N-1)g,e^-j2π(N-2)g,...,e^j2π(N-1)g]^T

则划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

则划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

其中，Φ＝diag(b)，Φ是矢量b对角线元素组成的对角矩阵。

在步骤3中，将划分后的第m个空间频率段的第l个多项式系数γ_ml的计算公式替换为：

γ_l＝[γ_1l,γ_2l,...,γ_ml,...,γ_M′l]^T

其中，F_M′是M′×M′维的傅里叶变换矩阵，b为MUSIC谱极值点对应的复多项式系数矢量，表示Hadamard积，0_{(M′-2N+1)×1}表示M′-2N+1维的全零列向量，M′表示划分的空间频率段的个数，l取1至L+1，L为设定的多项式阶数，划分后的第m个空间频率段的第l个多项式系数γ_ml为列向量γ_l的第m个元素。

所述步骤4的具体子步骤为：

(4.1)根据第m个空间频率段的所有L+1个多项式系数，用Newton-Raphson方法求解以下方程的实根：

$\underset{l=1}{\overset{L+1}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ml}}{g}^{l-1}=0$

其中，l取1至L+1，L为设定的多项式阶数；

求解出的以上方程的实根为第m个空间频率段的MUSIC谱函数的极点；将m从1至M′进行遍历，得出M′个空间频率段的MUSIC谱函数的极点；

(4.2)将子步骤(4.1)得出的第m个空间频率段的MUSIC谱函数的每个极点代入划分后的第m个空间频率段中频率点f_m(g)的表达式中，得到第m个空间频率段的对应频率点；所述划分后的第m个空间频率段中频率点f_m(g)的表达式为：

f_m(g)＝[-M′/2+(m-1)]ε+g

然后，在MUSIC谱函数h(f)中，令f分别取第m个空间频率段的每个频率极点，得出第m个空间频率段的MUSIC谱极值；

将m从1至M′进行遍历，得出M′个空间频率段的MUSIC谱函数的MUSIC谱极值。

本发明的有益效果为：第一，本发明用傅里叶变换快速实现复多项式系数矢量的计算，与传统方法相比，具有高效的优点。第二，本发明利用在有限空间带宽里的任一导向矢量可以用低阶多项式矢量表示的性质，把求MUSIC谱的极点转化为求多组低阶多项式的根。由于低阶多项式求根具有计算高效、数值稳定的特点，所以本发明与现有方法相比，具有低的计算复杂度和好的数值稳定性的优点。第三，本发明通过对信号子空间相关矢量做加窗傅里叶变换得到所有组的多项式系数，因此，大大降低估计多项式系数的计算复杂度，具有高效的优点。第四，本发明不仅适用于阵元域求根MUSIC，而且适用于波束空间求根MUSIC或子阵求根MUSIC。并且，对于波束形成矩阵或子阵合成矩阵没有特殊要求，因此，提出的新方法在低复杂度估计波达方向方面有巨大的潜力。

附图说明

图1为本发明的实多项式求根实现均匀线性阵列的谱MUSIC方法的流程图；

图2a为仿真实验1第一种情况下拟合误差随着归一化空间频率变化的曲线图；

图2b为仿真实验1第二种情况下(划分的空间频率段的个数为23)拟合误差随着归一化空间频率变化的曲线图；

图3为仿真实验2利用本发明得出的MUSIC谱函数的极点和MUSIC谱图；

图4为仿真实验2采用两种不同方法时得出的信噪比和MUSIC谱函数的极点的均方根误差的关系示意图；

图5为仿真实验2采用本发明得出的设定的多项式阶数和MUSIC谱函数的极点的均方根误差的关系图。

具体实施方式

下面结合附图对本发明作进一步说明：

参照图1，为本发明的实多项式求根实现均匀线性阵列的谱MUSIC方法的流程图。本发明实施例中，雷达的天线为均匀线阵。对于雷达的均匀线阵，其阵元数为N，阵元间隔为d。本发明的实多项式求根实现均匀线性阵列的谱MUSIC方法包括以下步骤：

步骤1，利用雷达接收回波数据，得出回波数据的协方差矩阵根据回波数据的协方差矩阵计算出MUSIC谱极值点对应的复多项式系数矢量b。

其具体子步骤为：

(1.1)建立雷达接收的回波数据模型。

具体地，入射至雷达均匀线阵的信号的个数表示为P，P为大于1的自然数。入射至雷达均匀线阵的信号为远场信号，则雷达第q次快拍获取的回波数据矢量x(q)为：

$x\left(q\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{s}_p\left(q\right)a\left(\mathrm{\θ}\right)+n\left(q\right),q=\mathrm{1,2},...,Q$

其中，q取1至Q，Q为雷达的快拍次数，Q>P；θ为入射至雷达均匀线阵的每个信号的入射角(指对应信号与雷达均匀线阵法线方向的夹角)，每个信号的入射角相同，s_p(q)是入射至雷达均匀线阵的第p个信号的复幅度，p取1至P；a(θ)表示入射至雷达均匀线阵的每个信号的导向矢量，n(q)表示高斯白噪声矢量，n(q)∈C^N×1。

并且有，

$a\left(\mathrm{\θ}\right)={\left[\begin{array}{cccc}1& e^{j2\mathrm{\π}\frac{d}{\mathrm{\λ}}\sin \mathrm{\θ}}& ...& e^{j2\mathrm{\π}(N-1)\frac{d}{\mathrm{\λ}}\sin \mathrm{\θ}}\end{array}\right]}^T$

其中，上标T表示矩阵或向量的转置，λ是雷达发射信号的波长。

(1.2)利用雷达每次快拍获取的回波数据矢量，估计出回波数据的协方差矩阵为：

$\hat{R}=\frac{1}{Q}\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}x\left(q\right)x^H\left(q\right)$

其中，上标H表示取共轭转置，q取1至Q，Q为雷达的快拍次数；

对回波数据的协方差矩阵进行特征分解，在回波数据的协方差矩阵的特征值中，选取P个最大的特征值，将选取的P个最大的特征值对应的特征矢量分别表示为u_s1至u_sP；利用选取的P个最大的特征值组成信号子空间U_s，U_s＝[u_s1,...,u_sP]，在回波数据的协方差矩阵的其余特征值(回波数据的协方差矩阵的特征值中排除P个最大的特征值)中，选取P个最大的特征值对应的特征矢量组成噪声子空间U_n。即为：

$\hat{R}=U_s{\mathrm{\Λ}}_s{U}_s^H+U_n{\mathrm{\Λ}}_n{U}_n^H$

其中，Λ_s是由信号子空间U_s对应的特征值组成的对角矩阵(信号子空间U_s对应的特征值按顺序排列在Λ_s的主对角线上)，Λ_n是由噪声子空间U_n对应的特征值组成的对角矩阵(噪声子空间U_n对应的特征值按顺序排列在Λ_n的主对角线上)。

(1.3)计算MUSIC谱极值点对应的复多项式系数矢量b。

具体地说，在特征矢量u_s1至u_sP中，针对每个特征矢量作傅里叶运算，得出对应特征矢量的空域响应。特征矢量u_sp的傅里叶运算过程表示为：

${\mathrm{\δ}}_p=F_M^{*}\left[\begin{array}{c}{u}_{\mathrm{sp}}\\ 0_{(M-N)\×1}\end{array}\right]$

其中，F_M是M×M维的傅里叶变换矩阵，*表示共轭，ceil(·)表示向上取整，u_sp表示回波数据的协方差矩阵的P个最大的特征值中第p个特征值对应的特征矢量，u_sp是N×1维矢量；0_(M-N)×1表示M-N维的全零列向量。

然后，根据以下公式得出MUSIC谱极值点对应的复多项式系数矢量b，

η＝[-jπ(N-1)，-jπ(N-2) … 0 … jπ(N-2)，jπ(N-1)]^Τ

${\stackrel{\‾}{b}}_0={\left[\begin{array}{ccccccc}0& ...& 0& N& 0& ...& 0\end{array}\right]}^T$

表示Hadamard积，上标T表示矩阵或向量的转置，N为雷达的均匀线阵的阵元数，∈R^(2N-1)×1(为2N-1维的列向量)，的第N个元素为N，其余元素均为0；为矢量中非零元素构成的矢量，即：

${\stackrel{\‾}{b}}_1=\left[\begin{array}{c}{\stackrel{\‾}{b}}_2\\ 0_{(M-(2N-1)\×1)}\end{array}\right]$

矢量 ${\stackrel{\‾}{b}}_1\∈C^{M\×1},$ 矢量为：

${\mathrm{\Δ}}_M=\begin{array}{cccc}\left(\right[\exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(-\frac{M}{2})\right)& \exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(-\frac{M}{2}+1)\right)& ...& \exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(\frac{M}{2}-1)\right)\left]\right)\end{array}$

其中，*表示共轭，表示傅里叶逆变换矩阵(M×M维)，diag(·)表示求对角矩阵。

将MUSIC谱极值点对应的复多项式系数矢量b表示为：

b＝[b₁,...,b_i...,b_2N-1]

其中，b_i为多项式矢量b的第i个元素，i取1至2N-1，则有*表示共轭，则MUSIC谱极值点对应的复多项式系数矢量b表示为：

$b={\left[\begin{array}{ccccccc}{b}_1& ...& b_{N-1}& b_N& b_{N-1}^{*}& ...& b_1^{*}\end{array}\right]}^T.$

步骤2，根据MUSIC谱极值点对应的复多项式系数矢量b，构造MUSIC谱函数；把空间频率划分为M′段，M′为大于1的自然数，划分后的每个空间频率段的长度为ε；然后根据以下公式得出多项式适应矩阵B：

$\underset{B}{\min }\underset{g\∈[0,\mathrm{\ϵ}]}{\mathrm{\Σ}}\mathrm{\ζ}\left(g\right)$

其中，

$\mathrm{\ζ}\left(g\right)=\frac{{\left|\right|c\left(g\right)-\mathrm{Bg}\left|\right|}_2^2}{{\left|\right|c\left(g\right)\left|\right|}_2^2},g={\left[\begin{array}{cccc}{g}^L& g^{L-1}& ...& 1\end{array}\right]}^T$

c(g)＝[e^-j2π(N-1)g,e^-j2π(N-2)g,...,e^j2π(N-1)g]^T

其中，多项式适应矩阵B为(2N-1)×(L+1)维的矩阵，N为雷达的均匀线阵的阵元数，L为设定的多项式阶数，g∈[0,ε)，||·||₂表示l₂范数，上标T表示矩阵或向量的转置。

其具体子步骤为：

(2.1)构造MUSIC谱函数。

具体地说，构造MUSIC谱函数h(f)，h(f)＝b^Τc(f)，其中，b为MUSIC谱极值点对应的复多项式系数，c(f)为：

c(f)＝[e^-j2π(N-1)f … 1 … e^j2π(N-1)f]^Τ

其中，f表示频率变量，N为雷达的均匀线阵的阵元数。

(2.2)把空间频率划分为M′段，M′为大于1的自然数，划分后的每个空间频率段的长度相等，划分后的每个空间频率段的长度表示为ε。则划分后的第m个空间频率段中频率点f_m(g)的表达式为：

f_m(g)＝[-M′/2+(m-1)]ε+g

其中，m取1至M′，g∈[0,ε)；

然后将划分后的第m个空间频率段导向矢量定义为：

其中，上标T表示矩阵或向量的转置，表示Hadamard积，并且有：

w_m＝[e^{-j2π(N-1)[-M′/2+(m-1)]ε},e^{-j2π(N-2)[-M′/2+(m-1)]ε},...,e^{j2π(N-1)[-M′/2+(m-1)]ε}]^T

c(g)＝[e^-j2π(N-1)g,e^-j2π(N-2)g,...,e^j2π(N-1)g]^T

则划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

如果α₁、α₂和α₃是三个同维列向量，则有因此，划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

其中，Φ＝diag(b)，Φ是矢量b对角线元素组成的对角矩阵(矢量b对角线元素按顺序排列在Φ的主对角线上。

(2.3)因为c(g)只与雷达的均匀线阵的阵元数N有关，与雷达接收的回波数据无关，所以根据c(g)，并采用最小二乘方法得出多项式适应矩阵B，多项式适应矩阵B为：

$\underset{B}{\min }\underset{g\∈[0,\mathrm{\ϵ}]}{\mathrm{\Σ}}\mathrm{\ζ}\left(g\right)$

其中，g∈[0,ε)，并且有：

$\mathrm{\ζ}\left(g\right)=\frac{{\left|\right|c\left(g\right)-\mathrm{Bg}\left|\right|}_2^2}{{\left|\right|c\left(g\right)\left|\right|}_2^2}$

g＝[g^L g^L-1 … 1]^Τ

$c\left(g\right)={[e^{-j2\mathrm{\π}(N-1)g},e^{-j2\mathrm{\π}(N-2)g},...,e^{j2\mathrm{\π}(N-1)g}]}^T$

B＝[t₁，t₂，...，t_L+1]

其中，||·||₂表示l₂范数，上标T表示矩阵或向量的转置，L为设定的多项式阶数，L为大于1的自然数；ζ(g)表示g的估计误差，t_l表示待求解的第l个窗矢量(维数为2N-1)，l取1至L+1。

在求解出多项式适应矩阵B之后，根据多项式适应矩阵B，得出第1个窗矢量t₁至第L+1个窗矢量t_L+1。

在以上求解多项式适应矩阵B的过程中，由于则多项式适应矩阵B为：B＝c(g)g^T(gg^T)^-1，其中，上标T表示矩阵或向量的转置，上标-1表示矩阵的逆。则有：c(g)＝Bg，在得出多项式适应矩阵B之后，即可得出每个窗矢量。

步骤3，得出划分后的每个空间频率段的每个多项式系数，得出划分后的第m个空间频率段的第l个多项式系数γ_ml的过程包括为：

由于c(g)＝Bg，则划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

$\begin{array}{c}{h}_m\left(g\right)={W_m}^T\mathrm{\ΦBg}\\ =\underset{l=1}{\overset{L+1}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ml}}{g}^{l-1}\end{array}$

其中，γ_ml表示划分后的第m个空间频率段的第l个多项式系数，l取1至L+1，L为设定的多项式阶数，Φ＝diag(b)，B为多项式适应矩阵，t_l为多项式适应矩阵B的第l列。

直接计算γ_ml复杂度高，由于w_m是一个离散傅里叶变换向量，这样，所有M′个空间频率段的第l个系数可以通过FFT运算(傅里叶变换)一次得到。因此，可以根据以下公式计算划分后的第m个空间频率段的第l个多项式系数γ_ml：

γ_l＝[γ_1l,γ_2l,...,γ_ml,...,γ_M′_l]^T

其中，F_M′是M′×M′维的傅里叶变换矩阵，b为MUSIC谱极值点对应的复多项式系数矢量，表示Hadamard积，t_l表示第l个窗矢量，0_{(M′-2N+1)×1}表示M′-2N+1维的全零列向量，M′表示划分的空间频率段的个数，l取1至L+1，L为设定的多项式阶数，划分后的第m个空间频率段的第l个多项式系数γ_ml为列向量γ_l的第m个元素。

步骤4，根据划分后的每个空间频率段的所有L+1个多项式系数，得出划分后的每个空间频率段的MUSIC谱函数的极点，根据划分后的每个空间频率段的MUSIC谱函数的极点，得出划分后的每个空间频率段的MUSIC谱函数的极值。

(4.1)根据第m个空间频率段的所有L+1个多项式系数，用Newton-Raphson方法求解以下方程的实根：

$\underset{l=1}{\overset{L+1}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ml}}{g}^{l-1}=0$

其中，l取1至L+1，L为设定的多项式阶数。

求解出的以上方程的实根为第m个空间频率段的MUSIC谱函数的极点。将m从1至M′进行遍历，得出M′个空间频率段的MUSIC谱函数的极点。

(4.2)将子步骤(4.1)得出的第m个空间频率段的MUSIC谱函数的每个极点代入划分后的第m个空间频率段中频率点f_m(g)的表达式中，得到第m个空间频率段的对应频率点。所述划分后的第m个空间频率段中频率点f_m(g)的表达式为：

f_m(g)＝[-M′/2+(m-1)]ε+g

然后，在MUSIC谱函数h(f)中，令f分别取第m个空间频率段的每个频率极点，得出第m个空间频率段的MUSIC谱极值(将第m个空间频率段的每个频率极点代入MUSIC谱函数中，得出的MUSIC谱函数的值为第m个空间频率段的MUSIC谱极值)。

将m从1至M′进行遍历，得出M′个空间频率段的MUSIC谱函数的MUSIC谱极值。

本发明实施例中，得出的M′个空间频率段的MUSIC谱函数的MUSIC谱极值用于雷达接收的回波的波达方向。其过程为：设定极值门限，根据大于极值门限的MUSIC谱极值，确定雷达接收的回波的波达方向。

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

仿真实验1，本实验研究拟合误差与多项式阶数和空间频率分段数的关系。

实验条件

仿真在MATLAB7.0软件下进行，仿真参数如下：第一种情况，设定的多项式阶数L分别为2、5和8，雷达的均匀线阵的阵元数N为12，划分的空间频率段的个数M′为92。第二种情况，设定的多项式阶数L分别为2、5和8，雷达的均匀线阵的阵元数N为12，划分的空间频率段的个数M′为23。

实验结果

参照图2a，为仿真实验1第一种情况下(划分的空间频率段的个数为92)拟合误差随着归一化空间频率变化的曲线图。参照图2b，为仿真实验1第二种情况下(划分的空间频率段的个数为23)拟合误差随着归一化空间频率变化的曲线图。在图2a和图2b中，横轴表示归一化空间频率，纵轴表示拟合误差，单位为dB。从图2a和图2b可以看出，拟合误差与设定的多项式阶数L成反比，同时与划分的空间频率段的个数成反比。也就是说，设定的多项式阶数或划分的空间频率段的个数越大，则拟合误差越小

仿真实验2，本实验研究本发明提出的实多项式求根实现均匀线性阵列的谱MUSIC的性能。

实验条件

雷达天线为一个均匀线阵，阵元数为16，阵元间距为半波长。有四个互不相关的窄带信源分别从-20°、-2°、2°、35°入射。所有信源阵元级的信噪比是0dB，其中信噪比定义为 和分别是阵元级信号功率和噪声功率。

实验结果

利用本发明和传统的求根MUSIC估计器分别得出MUSIC谱函数的极值。参照图3，为仿真实验2利用本发明得出的MUSIC谱函数的极点和MUSIC谱图。图3中，横轴表示角度，单位为度，纵轴表示MUSIC谱函数的值，单位为dB。图3中，MUSIC谱用实线表示，MUSIC谱函数的极点用星号表示。从图3可以看出，MUSIC谱函数的各个极点均被解出，证明了本发明的有效性。

参照图4，为仿真实验2采用两种不同方法时得出的信噪比和MUSIC谱函数的极点的均方根误差的关系示意图。图4中，横轴表示信噪比，单位为dB。纵轴表示MUSIC谱函数的极点的均方根误差，单位为dB。图4中，本发明(对应图4中的多窗MUSIC)用带圆圈的线条表示，传统的求根MUSIC估计器(对应图4中的传统MUSIC)用带星号的线条表示。从图4中可以看出，对于所有的信噪比，本发明和传统的求根MUSIC估计器得出的MUSIC谱函数的极点的均方根误差几乎没有差别，因此，本发明可以有效地实现MUSIC谱函数的极点的估计，并且有较低的计算复杂度和更好的数值稳定性。

参照图5，为仿真实验2采用本发明得出的设定的多项式阶数和MUSIC谱函数的极点的均方根误差的关系图。图5中，横轴表示设定的多项式阶数，纵轴表示MUSIC谱函数的极点的均方根误差，单位为dB。图5中，实线表示采用本发明得出的MUSIC谱函数的极点的均方根误差，圆圈表示采用传统的求根MUSIC估计器得出的MUSIC谱函数的极点的均方根误差。从图5中可以看出，MUSIC谱函数的极点的均方根误差随着多项式阶数增大而减小。并且，当多项式阶数大于4时，采用本发明和统的求根MUSIC估计器得出的MUSIC谱函数的极点的均方根误差几乎没有差别。

显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则本发明也意图包含这些改动和变型在内。

Claims (5)

1.实多项式求根实现均匀线性阵列的谱MUSIC方法，其特征在于，包括以下步骤：

步骤1，利用雷达接收回波数据，得出回波数据的协方差矩阵根据回波数据的协方差矩阵计算出MUSIC谱极值点对应的复多项式系数矢量b；

步骤2，根据MUSIC谱极值点对应的复多项式系数矢量b，构造MUSIC谱函数；把空间频率划分为M′段，M′为大于1的自然数，划分后的每个空间频率段的长度为ε；然后根据以下公式得出多项式适应矩阵B：

$\underset{B}{\min }\underset{g\∈[0,\mathrm{\ϵ}]}{\mathrm{\Σ}}\mathrm{\ζ}\left(g\right)$

其中，

$\mathrm{\ζ}\left(g\right)=\frac{{\left|\right|c\left(g\right)-\mathrm{Bg}\left|\right|}_2^2}{{\left|\right|c\left(g\right)\left|\right|}_2^2},g={\left[\begin{array}{cccc}{g}^L& g^{L-1}& ...& 1\end{array}\right]}^T$

c(g)＝[e^-j2π(N-1)g,e^-j2π(N-2)g，...，e^j2π(N-1)g]^T

其中，多项式适应矩阵B为(2N-1)×(L+1)维的矩阵，N为雷达的均匀线阵的阵元数，L为设定的多项式阶数，g∈[0,ε)，||·||₂表示l₂范数，上标T表示矩阵或向量的转置；

步骤3，根据以下公式计算划分后的每个空间频率段的每个多项式系数，划分后的第m个空间频率段的第l个多项式系数γ_ml的计算公式为：

其中，m取1至M′，l取1至L+1，表示Hadamard积，t_l为多项式适应矩阵B的第l列，w_m为：

w_m＝[e^{-j2π(N-1)[-M′/2+(m-1)]ε},e^{-j2π(N-2)[-M′/2+(m-1)]ε},...,e^{j2π(N-1)[-M′/2+(m-1)]ε}]^T；

步骤4，根据划分后的每个空间频率段的所有L+1个多项式系数，得出划分后的每个空间频率段的MUSIC谱函数的极点，根据划分后的每个空间频率段的MUSIC谱函数的极点，得出划分后的每个空间频率段的MUSIC谱函数的极值。

2.如权利要求1所述的实多项式求根实现均匀线性阵列的谱MUSIC方法，其特征在于，所述步骤1的具体子步骤为：

(1.1)将入射至雷达均匀线阵的信号的个数表示为P，P为大于1的自然数；则雷达第q次快拍获取的回波数据矢量x(q)为：

$x\left(q\right)=\underset{p=1}{\overset{P}{\mathrm{\Σ}}}{s}_p\left(q\right)a\left(\mathrm{\θ}\right)+n\left(q\right)$

其中，q取1至Q，Q为雷达的快拍次数；θ为入射至雷达均匀线阵的每个信号的入射角，s_p(q)是入射至雷达均匀线阵的第p个信号的复幅度，p取1至P；a(θ)表示入射至雷达均匀线阵的每个信号的导向矢量，n(q)表示高斯白噪声矢量；

(1.2)利用雷达每次快拍获取的回波数据矢量，估计出回波数据的协方差矩阵为：

$\hat{R}=\frac{1}{Q}\underset{q=1}{\overset{Q}{\mathrm{\Σ}}}x\left(q\right)x^H\left(q\right)$

其中，上标H表示取共轭转置；

对回波数据的协方差矩阵进行特征分解，在回波数据的协方差矩阵的特征值中，选取P个最大的特征值，将选取的P个最大的特征值对应的特征矢量分别表示为u_s1至u_sP；利用选取的P个最大的特征值组成信号子空间U_s，U_s＝[u_s1,...,u_sP]，在回波数据的协方差矩阵的其余特征值中，选取P个最大的特征值对应的特征矢量组成噪声子空间U_n；

(1.3)在特征矢量u_s1至u_sP中，针对每个特征矢量作傅里叶运算，特征矢量u_sp的傅里叶运算过程表示为：

${\mathrm{\δ}}_p=F_M^{*}\left[\begin{array}{c}{u}_{\mathrm{sp}}\\ 0_{(M-N)\×1}\end{array}\right]$

其中，F_M是M×M维的傅里叶变换矩阵，*表示共轭，ceil(·)表示向上取整，u_sp表示回波数据的协方差矩阵的P个最大的特征值中第p个特征值对应的特征矢量；0_(M-N)×1表示M-N维的全零列向量；

根据以下公式得出MUSIC谱极值点对应的复多项式系数矢量b，

η＝[-jπ(N-1)，-jπ(N-2) … 0 … jπ(N-2)，jπ(N-1)]^T

${\stackrel{\‾}{b}}_0={\left[\begin{array}{ccccccc}0& ...& 0& N& 0& ...& 0\end{array}\right]}^T$

其中，表示Hadamard积，上标T表示矩阵或向量的转置，N为雷达的均匀线阵的阵元数，为2N-1维的列向量，的第N个元素为N，其余元素均为0；为矢量中非零元素构成的矢量，即：

${\stackrel{\‾}{b}}_1=\left[\begin{array}{c}{\stackrel{\‾}{b}}_2\\ 0_{(M-(2N-1)\×1)}\end{array}\right]$

矢量为：

$\begin{array}{cccc}{\mathrm{\Δ}}_M=\mathrm{diag}\left(\right[\exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(-\frac{M}{2})\right)& \exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(-\frac{M}{2}+1)\right)& ...& \exp \left(j\frac{2\mathrm{\π}(N-1)}{M}(\frac{M}{2}-1)\right)\left]\right)\end{array}$ 其中，*表示共轭，表示M×M维的傅里叶逆变换矩阵，diag(·)表示求对角矩阵。

3.如权利要求1所述的实多项式求根实现均匀线性阵列的谱MUSIC方法，其特征在于，在步骤2中，构造的MUSIC谱函数表示为h(f)，h(f)＝b^Τc(f)，c(f)为：

c(f)＝[e^-j2π(N-1)f e^-j2π(N-2)f … 1 … e^j2π(N-1)f]^Τ

其中，f表示频率变量；

在步骤2中，划分后的第m个空间频率段中频率点f_m(g)的表达式为：

f_m(g)＝[-M′/2+(m-1)]ε+g

其中，m取1至M′，g∈[0,ε)；

然后将划分后的第m个空间频率段导向矢量定义为：

其中，上标T表示矩阵或向量的转置，表示Hadamard积，并且有：

c(g)＝[e^-j2π(N-1)g,e^-j2π(N-2)g,...,e^j2π(N-1)g]^T

则划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

则划分后的第m个空间频率段的MUSIC谱函数h_m(g)为：

其中，Φ＝diag(b)，Φ是矢量b对角线元素组成的对角矩阵。

4.如权利要求1所述的实多项式求根实现均匀线性阵列的谱MUSIC方法，其特征在于，在步骤3中，将划分后的第m个空间频率段的第l个多项式系数γ_ml的计算公式替换为：

γ_l＝[γ_1l,γ_2l,...,γ_ml,...,γ_M′l]^T

其中，F_M′是M′×M′维的傅里叶变换矩阵，b为MUSIC谱极值点对应的复多项式系数矢量，表示Hadamard积，0_{(M′-2N+1)×1}表示M′-2N+1维的全零列向量，M′表示划分的空间频率段的个数，l取1至L+1，L为设定的多项式阶数，划分后的第m个空间频率段的第l个多项式系数γ_ml为列向量γ_l的第m个元素。

5.如权利要求1所述的实多项式求根实现均匀线性阵列的谱MUSIC方法，其特征在于，在步骤2中，构造的MUSIC谱函数表示为h(f)，f表示频率变量；

所述步骤4的具体子步骤为：

(4.1)根据第m个空间频率段的所有L+1个多项式系数，用Newton-Raphson方法求解以下方程的实根：

$\underset{l=1}{\overset{L+1}{\mathrm{\Σ}}}{\mathrm{\γ}}_{\mathrm{ml}}{g}^{l-1}=0$

其中，l取1至L+1，L为设定的多项式阶数；

求解出的以上方程的实根为第m个空间频率段的MUSIC谱函数的极点；将m从1至M′进行遍历，得出M′个空间频率段的MUSIC谱函数的极点；

(4.2)将子步骤(4.1)得出的第m个空间频率段的MUSIC谱函数的每个极点代入划分后的第m个空间频率段中频率点f_m(g)的表达式中，得到第m个空间频率段的对应频率点；所述划分后的第m个空间频率段中频率点f_m(g)的表达式为：

f_m(g)＝[-M′/2+(m-1)]ε+g

然后，在MUSIC谱函数h(f)中，令f分别取第m个空间频率段的每个频率极点，得出第m个空间频率段的MUSIC谱极值；

将m从1至M′进行遍历，得出M′个空间频率段的MUSIC谱函数的MUSIC谱极值。