CN104155629B - 一种冲击噪声背景下小快拍数信号波达方向估计方法 - Google Patents

Abstract

本发明涉及一种冲击噪声环境下小快拍数信号波达方向估计方法。本发明包括：采集信号采样数据；对接收到的信号数据矢量进行去冲击预处理；对去冲击后的矩阵的协方差矩阵；将经过重采样构建后的协方差矩阵分解子空间得到其噪声子空间；对B次重采样的协方差矩阵分别进行求解，得到由B个入射角度估计值构成的矢量矩阵，利用置信空间对矩阵进行校正后求均值即可得到最终波达方向角度的估计值本发明经过对接收数据的预处理，能对冲击噪声进行有效地抑制；采用block‑bootstrap方法对数据矩阵进行重采样构建，在小快拍数下能对信号进行有效地扩展；该方法稳健性高，适用于冲击噪声背景下的小快拍数DOA估计，具有良好的测向效果。

Description

一种冲击噪声背景下小快拍数信号波达方向估计方法

技术领域

本发明涉及一种冲击噪声环境下小快拍数(snapshot)信号波达方向(Direction ofArrival，简称DOA)估计方法。

背景技术

DOA估计是阵列信号处理中的一个重要研究领域，在电子战、雷达、声呐、射电天文学、地震学、卫星以及移动通信系统等方面有着广泛的应用。

传统的DOA算法针对的是固定信源，设定为已知分布为高斯白噪声条件下的大快拍数估计，然而，在实际应用中噪声往往不服从高斯分布或具有冲击性，如海杂波噪声，大气放电噪声，无线信道瞬时干扰噪声等，所以，对于未知参数特性的非高斯白噪声环境下的DOA算法是空间谱估计理论中一个重要的课题。另外，实际中由于种种原因，并不一定能够采样到足够的信号形成大的快拍数，因此，在小快拍数下的DOA估计问题也是波达方向估计中一个重要的课题。

对于DOA估计，一般性方法是谱估计，即利用阵列形式对空间波达信号的角度功率谱，例如1979年,Schmit等人在“Multiple emitter location and signal parameterestimation”中提出了著名的MUSIC(多重信号分类)方法，后来Roy等人在“Esprit-asubspace ration approach to estimation of parameters of cissoids in noise”中提出的ESPRIT(旋转不变子空间)方法，都是利用谱进行估计，进而得到用户信号的波达角度。MUSIC算法和ESPRIT算法都属于子空间类算法，子空间类算法是通过对接收信息的处理得到信号子空间或噪声子空间，通过构造谱函数进而得到空间波达信号的角度功率谱。其中MUSIC算法利用噪声子空间而ESPRIT利用的是信号子空间。以MUSIC算法为代表的算法包括特征矢量法、MUSIC、root-MUSIC法及MNM等，以ESPRIT算法为代表的算法主要有TAM、LS-ESPRIT及TLS-ESPRIT等。

对于复杂电磁环境下的测向研究，近年来先后出现了利用共变矩阵(ROC)、分数低阶矩(FOLM)等来估计冲击噪声环境下的测向算法，以及利用特殊阵列模型(Toeplitz重构等)、参数化的噪声模型(将噪声视为AR、MA等模型)来估计色噪声下的测向算法。以上算法较为复杂，运算量大，需要已知或者预先假设参数。对于小快拍数下的DOA估计，一般采用智能优化算法(粒子群算法(PSO))、阵列插值、阵列变换、空间重采样等统计学方法来计算。其算法的复杂程度、迭代次数和分辨力高低都各有缺点。因此为了更好在复杂电磁环境下更有效的DOA估计，下述问题亟待解决：

1)复杂噪声环境下的测向。传统算法和子空间类算法基本上都是假设噪声是白高斯，或者是已知统计特性的色噪声。然而,未知统计特性的色噪声环境在实际应用中却经常出现。此时传统的MUSIC等子空间算法性能急剧下降。

2)小快拍数下的测向。现实生活中电磁环境越来越复杂，尤其是在电子信息战中。由于各种干扰的存在，单位时间内获得的快拍数有限，已有的大多数算法在此情况下性能严重下降，甚至不能测向。所以小块拍数下的测向研究必须借助于其他参数估计等手段来完成。

发明内容

本发明的目的在于提供一种稳定性更高的冲击噪声背景下小快拍数信号波达方向估计方法。

本发明的目的是这样实现的：

(1)采集信号采样数据：

X(t)＝AS(t)+N(t)，

其中，X(t)＝[x₁(t),x₂(t),...,x_L(t)]^T为快拍数为t时的L×1维接收数据矢量，A为阵列流型矩阵，S(t)为Μ×1维窄带信号矢量，N(t)为L×1维加性噪声矢量，噪声类型为复冲击噪声，A＝[a(θ₁),a(θ₂),...,a(θ_M)]，为导向矢量，θ_i为第i个信源的入射角度，v_i＝exp(-jπsin(θ_i)),i＝1,2,...,M；

(2)对接收到的信号数据矢量进行去冲击预处理：

判断接收阵列数据矢量矩阵绝对值最大值X_imax与门限值之间的大小，其中门限值为为接收阵列数据矢量矩阵绝对值的均值，如果X_imax大于等于门限值，则用替换矩阵中最大值，如小于门限值则执行下一步骤；重复上述比较过程，直至矩阵中的值均小于门限值；

$\left\{\begin{array}{c}{X}_{i\max }\&GreaterEqual;a{\stackrel{\&OverBar;}{X}}_i,X_{i\max }=b{\stackrel{\&OverBar;}{X}}_i\\ X_{i\max }<a{\stackrel{\&OverBar;}{X}}_i,X_{i\max }=X_{i\max }\end{array}\right.,$

其中，b＝0.6(2+α)，α为冲击噪声特征系数，α∈(0,2]，GSNR是广义信噪比；

(3)对去冲击后的矩阵的协方差矩阵R_X：

R_X＝E[XX^H]，

以协方差矩阵中每一列为单位进行可放回的随机均匀重采样B次后重构协方差矩阵R_X ^*，B＝200～1000，B∝1/N，N为快拍数；

(4)将经过重采样构建后的协方差矩阵R_X ^*分解子空间得到其噪声子空间U_N，

$f\left(z\right)=e_i^{H}p\left(z\right),i=M+1,...,L$

式中e_i ^H是数据协方差矩阵R_X ^*中小特征值对应的L-M个特征矢量即代表着噪声子空间，其中

p(z)＝[1 z…z^L-1]^Τ

求解多项式的根：

z_i＝exp(jw_i)(1≤i≤M)

得到入射角度：

${\mathrm{\&theta;}}_i=\arcsin \left(\frac{\mathrm{\&lambda;}}{2\mathrm{\&pi;d}}\arg \right\{{\hat{z}}_i\left\}\right),i=1,...,M;$

(5)对B次重采样的协方差矩阵分别进行求解，得到由B个入射角度估计值构成的矢量矩阵，利用置信空间对矩阵进行校正后求均值即可得到最终波达方向角度的估计值

本发明的有益效果在于：

与现有的DOA方法相比，通过应用本发明的方法，经过对接收数据的预处理，能对冲击噪声进行有效地抑制；采用block-bootstrap方法对数据矩阵进行重采样构建，在小快拍数下能对信号进行有效地扩展；然后采用root-MUSIC算法对信号波达方向进行精测，重复采样及root-MUSIC算法B次，利用大数据理论逼近真实值。该算法稳健性高，适用于冲击噪声背景下的小快拍数DOA估计，具有良好的测向效果。

附图说明

图1是本发明的方法示意图。

图2是本发明中去冲击预处理方法示意图。

图3是在信号独立，阵元数为4，广义信噪比为5dB，阵元间距为半波长，快拍数为2，特征指数为0.9下的DOA测向仿真图。

图4是在信号独立，阵元数为10，阵元间距为半波长，快拍数为20，特征指数为1.2，不同广义信噪比下的DOA测向成功率仿真图。

具体实施方式

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

本发明要解决的技术问题是，克服现有技术的不足，提出一种能在冲击噪声环境下对小快拍数的DOA测向方法。

为解决上述问题，本发明基于均匀线阵(ULA)，提供一种新的DOA测向方法，应用于冲击噪声环境下的小快拍数的DOA估计，

(1)获取信号采样数据：

X(t)＝AS(t)+N(t)，

其中，X(t)＝[x₁(t),x₂(t),...,x_L(t)]^T为快拍数为t时的L×1维接收数据矢量，A为阵列流型矩阵，S(t)为Μ×1维窄带信号矢量，N(t)为L×1维加性噪声矢量，噪声类型为复冲击噪声，A＝[a(θ₁),a(θ₂),...,a(θ_M)]，为导向矢量，θ_i为第i个信源的入射角度，v_i＝exp(-jπsin(θ_i)),i＝1,2,...,M；

(2)对接收到的信号数据矢量进行去冲击预处理：判断接收阵列数据矢量矩阵绝对值最大值X_imax与门限值之间的大不小于小，其中门限值为为接收阵列数据矢量矩阵绝对值的均值，如果X_imax大于等于门限值，则用替换矩阵中最大值，如小于门限值则进入下一步骤；重复上述比较过程，直至矩阵中的值均小于门限值。

$\left\{\begin{array}{c}{X}_{i\max }\&GreaterEqual;a{\stackrel{\&OverBar;}{X}}_i,X_{i\max }=b{\stackrel{\&OverBar;}{X}}_i\\ X_{i\max }<a{\stackrel{\&OverBar;}{X}}_i,X_{i\max }=X_{i\max }\end{array}\right.,$

其中，b＝0.6(2+α)，α为冲击噪声特征系数，α∈(0,2]，GSNR是广义信噪比。

(3)对去冲击后的矩阵的协方差矩阵R_X：

R_X＝E[XX^H]，

利用Block-bootstrap方法以协方差矩阵中每一列为单位进行可放回的随机均匀重采样B次后重构协方差矩阵R_X ^*，B的数值取决于快拍数多少，B＝200～1000，B∝1/N，N为快拍数。

(4)将经过重采样构建后的协方差矩阵R_X ^*分解子空间得到其噪声子空间U_N，

定义：

$f\left(z\right)=e_i^{H}p\left(z\right),i=M+1,...,L$

式中e_i ^H是数据协方差矩阵R_X ^*中小特征值对应的L-M个特征矢量即代表着噪声子空间，其中

p(z)＝[1 z…z^L-1]^Τ

求解多项式的根：

z_i＝exp(jw_i)(1≤i≤M)

最终得到入射角度：

${\mathrm{\&theta;}}_i=\arcsin \left(\frac{\mathrm{\&lambda;}}{2\mathrm{\&pi;d}}\arg \right\{{\hat{z}}_i\left\}\right),i=1,...,M;$

(5)对B次重采样的协方差矩阵分别进行上述root-MUSIC方法求解，得到由B个入射角度估计值构成的矢量矩阵，利用置信空间对该矩阵进行校正后求其均值即可得到最终波达方向角度的估计值

此方法包括：

步骤(10)、对所述均匀线阵接收到的信号数据进行预处理；

步骤(20)、将处理后的数据阵列求解其协方差R_X；

步骤(30)、利用block-bootstrap方法对方差进行重采样得到R_X ^*；

步骤(40)、对重采样后的数据R_X ^*进行root-MUSIC运算，得到波达方向角度θ；

步骤(50)、重复步骤(30)到(40)B次，利用置信区间对得到的数据进行校正。

其中步骤(10)进一步包括：

步骤(110)、获取所述阵元接收的数据矢量，求其绝对值，并找出绝对值中的最大值X_imax和均值

步骤(120)、判断步骤(110)中最大值X_imax是否大于如大于则用替代，如不大于则保持原数据不变；

步骤(130)、重复步骤(110)至(120)，直到所有数据均小于

其中步骤(30)进一步包括：

步骤(310)、确立块采样(block-bootstrap)中每一个采样块(block)的大小；

步骤(320)、利用block-bootstrap规则对R_X进行随机重采样，构建新的R_X ^*。

其中，步骤(40)进一步包括：

步骤(410)、根据R_X ^*矩阵求出特征子空间后带入root-MUSIC公式求出多项式的根；

步骤(420)、找出单位圆上的根，求出其对应的信号入射角度。

其中，步骤(50)进一步包括：

步骤(510)、根据原始接收数据的快拍数大小，确立采样次数B；

步骤(520)、根据B次计算结果，利用置信空间校正最终结果。

其中，步骤(120)中，参数a和参数b由下式定义：

$a=\frac{20(1+\mathrm{\&alpha;})}{0.5+\sqrt{\mathrm{GSNR}}}$

b＝0.6(2+α)

上式中，α为冲击噪声特征系数，α∈(0,2]，GSNR是广义信噪比。

其中，步骤(310)中，根据接收阵元特征，确立以协方差矩阵的一列大小为采样数据的一个采样块。

其中，步骤(510)，根据接收数据的快拍数N的大小，确立采样次数B＝200～1000；B∝1/N。

Bootstrap(自助法)算法是一种在小样本数下进行参数估计的有效的统计方法。其基本思想是通过多次(200次以上)有放回的随机重采样进行重构数据，利用统计学规律逼近真实值。当样本来自总体，其分布能用正态分布来描述时，它的抽样分布(Sampling Distribution)也为正态分布(The Normal Distribution)；但当样本来自的总体无法以正态分布来描述时，则以渐进分析法、自助法等来分析。采用随机可置换抽样(random sampling with replacement)的方法。对于小数据集，自助法效果很好。

传统冲击噪声环境下的测向方法无论是运算量较大的最大似然估计法还是计算量较小的ROC-MUSIC或者是基于分数低阶矩的FLOM-MUSIC算法在小块拍数下均不能有效的测向。本发明通过冲击预处理后再采用Bootstrap方法重构数据矩阵，不光可以降低计算复杂度，还可以在小快拍数下有效的测向。

下面结合附图对本发明进行详细描述。

如图1所示，本发明的方法首先对接收数据矢量进行去冲击预处理，之后将经过处理的数据矢量求解数据协方差，然后利用块采样规则对协方差数据进行重采样重构，最后由root-MUSIC算法得到DOA估计。

本发明的阵列结构采用ULA(均匀线阵)，阵元个数为L，阵元间隔为入射波长一半。假设在L源各向同性阵列远场处，有M个窄带信号源以平面波入射，入射角分别为θ₁(t),θ₂(t),...,θ_M(t)。

构造阵列流型矩阵：

A＝[a(θ₁),a(θ₂),...,a(θ_M)]

其中，为导向矢量矩阵，v_i＝exp(-jπsin(θ_i)),i＝1,2,...,M。

M×1维窄带信号矢量矩阵为

S(t)＝[s₁(t),s₂(t),...,s_M(t)]^T

则接收数据矢量矩阵为

X(t)＝AS(t)+N(t)

其中，X(t)＝[x₁(t),x₂(t),...,x_L(t)]^T为L×1维接收数据矢量矩阵，N(t)＝[n₁(t),n₂(t),...,n_L(t)]^T为L×1维加性噪声矢量矩阵，噪声类型为复冲击噪声。此处的冲击噪声使用α稳定分布来建模。

传统的DOA测向方法往往针对的是高斯噪声背景，利用接收数据矢量的二阶矩或高阶矩来进行DOA估计。由于α稳定分布没有二阶或高阶矩，因此在冲击噪声背景下传统方法失效。本发明采用对接收数据矢量去冲击预处理的方法(如图2所示)，以降低冲击噪声对测向效果的影响。

首先对接收数据矩阵求绝对值后求出其均值然后求出阵列最大值X_imax。判断X_imax是否大于如大于则用替换最大值，而后重新搜索矩阵最大值，继续比较替换，直到矩阵中没有大于门限值的数值存在。

其中参数a和参数b由下式定义：

$a=\frac{20(1+\mathrm{\&alpha;})}{0.5+\sqrt{\mathrm{GSNR}}}$

b＝0.6(2+α)

上式中，α为冲击噪声特征系数，α∈(0,2]，GSNR是广义信噪比。

对经过去冲击预处理的接收数据矢量计算其协方差矩阵R_X:

R_X＝E[XX^H]

以R_X中每一列为单位进行重采样重构矩阵R_X ^*，对矩阵R_X ^*进行子空间分解得到其噪声子空间U_N，

定义：

$f\left(z\right)=e_i^{H}p\left(z\right),i=M+1,...,L$

式中e_i ^H数据协方差矩阵R_X ^*中小特征值对应的L-M个特征矢量，其中

p(z)＝[1 z…z^L-1]^Τ

求解多项式的根：

z_i＝exp(jw_i)(1≤i≤M)

最终得到入射角度：

${\mathrm{\&theta;}}_i=\arcsin \left(\frac{\mathrm{\&lambda;}}{2\mathrm{\&pi;d}}\arg \right\{{\hat{z}}_i\left\}\right),i=1,...,M$

根据重采样规则重复上述采样并求解B次后得到θ_i矢量矩阵，根据其置信空间校正后求其均值即为波达方向角度。测向效果如图3和图4所示。

与现有的非高斯环境下的DOA测向方法相比，采用本发明基于去冲击预处理的方法，能够有效地减少冲击噪声对DOA算法的影响，并且在噪声冲击性较强时有较好的测向性能。采用本发明将去冲击与块采样算法相结合的测向方法，能够在小快拍数下有效地提高测向性能，实现对小块拍数的有效测向。综上，本发明能够有效的实现在冲击噪声环境下的小块拍数DOA测向，该方法的稳定性较好。

除了上述的root-MUSIC谱估计外，本发明的去冲击后的块采样重构矩阵测向方法得出的噪声子空间也适用于任何子空间类的DOA方法。

Claims (1)

1.一种冲击噪声背景下小快拍数信号波达方向估计方法，其特征在于：

(1)采集信号采样数据：

X(t)＝AS(t)+N(t)，

其中，X(t)＝[x₁(t),x₂(t),...,x_L(t)]^T为快拍数为t时的L×1维接收数据矢量，A为阵列流型矩阵，S(t)为Μ×1维窄带信号矢量，N(t)为L×1维加性噪声矢量，噪声类型为复冲击噪声，A＝[a(θ₁)，a(θ₂)，...，a(θ_M)]，为导向矢量，θ_i为第i个信源的入射角度，v_i＝exp(-jπsin(θ_i)),i＝1,2,...,M；

(2)对接收到的信号数据矢量进行去冲击预处理：

判断接收阵列数据矢量矩阵绝对值最大值X_imax与门限值之间的大小，其中门限值为 为接收阵列数据矢量矩阵绝对值的均值，如果X_imax大于等于门限值，则用替换矩阵中最大值，如小于门限值则执行下一步骤；重复上述比较过程，直至矩阵中的绝对值均小于门限值；

$\left\{\begin{array}{cc}{X}_{i\max }\&GreaterEqual;a{\stackrel{\&OverBar;}{X}}_i,& X_{i\max }=b{\stackrel{\&OverBar;}{X}}_i\\ X_{i\max }<a{\stackrel{\&OverBar;}{X}}_i,& X_{i\max }=X_{i\max }\end{array}\right.,$

其中，b＝0.6(2+α)，α为冲击噪声特征系数，α∈(0,2]，GSNR是广义信噪比；

(3)对去冲击后的矩阵的协方差矩阵R_X：

R_X＝E[XX^H]，

以协方差矩阵中每一列为单位进行可放回的随机均匀重采样B次后重构协方差矩阵R_X ^*，B＝200～1000，B∝1/N，N为快拍数；

(4)将经过重采样构建后的协方差矩阵R_X ^*分解子空间得到其噪声子空间U_N，

$f\left(z\right)=e_i^{H}p\left(z\right),i=M+1,...,L$

式中e_i ^H是数据协方差矩阵R_X ^*中小特征值对应的L-M个特征矢量即代表着噪声子空间，其中

p(z)＝[1 z…z^L-1]^Τ

求解多项式的根：

z_i＝exp(jw_i) 1≤i≤M

得到入射角度：

${\mathrm{\&theta;}}_i=arcsin\left(\frac{\mathrm{\&lambda;}}{2\mathrm{\&pi;}d}\arg \right\{{\hat{z}}_i\left\}\right),i=1,...,M;$

(5)对B次重采样的协方差矩阵分别进行求解，得到由B个入射角度估计值构成的矢量矩阵，利用置信空间对矩阵进行校正后求均值即可得到最终波达方向角度的估计值