CN103901395A - 一种冲击噪声环境下相干信号波达方向动态跟踪方法 - Google Patents

Abstract

本发明涉及一种冲击噪声环境下相干信号波达方向动态跟踪方法。本发明包括：获取信号采样数据；对接收到的信号数据矢量进行类归一化预处理；将均匀线阵划分为若干个子阵；将经过前后向空间平滑的接收数据矢量带入PAST算法，得到信号子空间；对信号子空间的数据使用MUSIC算法进行处理，经谱峰搜索获得目标信号波达方向；令t=t+1，得到下一块拍数对应的信号波达方向角度，直到达到最大快拍数，实现冲击噪声环境下相干信号源的DOA动态跟踪。本发明的方法采用投影子空间逼近跟踪算法，能对入射信号的信号子空间进行实时地跟踪，然后采用MUSIC方法对信号波达方向进行精测，该方法稳健性高，跟踪效果好。

Description

一种冲击噪声环境下相干信号波达方向动态跟踪方法

技术领域

本发明涉及一种冲击噪声环境下相干信号波达方向（Direction of Arrival，简称DOA）动态跟踪方法。

背景技术

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

传统的DOA算法往往针对的是固定信源，然而，在实际应用中，信源的角度会随时间而变化，所以，针对入射角度随时间变化的信号源的动态DOA跟踪问题是空间谱估计理论应用中一个重要课题。另外，实际中的噪声往往具有很强的冲击性，如海杂波噪声，大气噪声，无线信道噪声等，因此，对于冲击噪声环境下的DOA跟踪是波达方向估计中一个重要的课题。

对于DOA估计，一般性方法是谱估计，即利用阵列形式对空间波达信号的角度功率谱，例如MUSIC（多重信号分类）谱和ESPRIT（旋转不变子空间）谱等进行估计，进而得到用户信号的波达角度。其中MUSIC算法和ESPRIT算法都属于子空间类算法，MUSIC算法属于噪声子空间类算法，ESPRIT算法属于信号子空间类算法。子空间类算法是通过对接收信息的处理得到信号子空间或噪声子空间，通过构造谱函数进而得到空间波达信号的角度功率谱。以MUSIC算法为代表的算法包括特征矢量法、MUSIC、求根MUSIC法及MNM等，以ESPRIT算法为代表的算法主要有TAM、LS-ESPRIT及TLS-ESPRIT等。

对于DOA跟踪，一般是通过对信号子空间的更新来实现对入射角度的跟踪。在众多的子空间更新算法中，Yang于1995年提出的PAST（投影逼近子空间跟踪）具有一定的代表性。该算法通过求解一个无约束最优化问题，实现了对信号子空间的跟踪。但是，PAST算法在冲击噪声的环境下难以实现对信号的有效跟踪，估计精度和成功概率随着噪声冲击性的增强而迅速下降。对于相干信号源，PAST算法随着信号源的相关性的增加，性能将迅速恶化，直至完全失效。在实际环境中，噪声往往具有很强的冲击性，同时，在无线信道的传播中，无线环境十分复杂，大量存在相关信号。因此，为了能对冲击噪声环境下的相干信号进行有效的跟踪，下述问题亟待解决：

1）冲击噪声环境下的测向。传统的DOA估计方法是对接收数据矢量的二阶距进行特征值分解，得到噪声子空间或信号子空间，进而采用子空间类算法进行谱峰搜索得到信号入射角度。但是由于噪声环境的冲击性，接收数据矢量的二阶距是不存在的，所以传统的DOA估计方法在冲击噪声环境下是失效的。

2）相干信号源的跟踪。当信号源完全相干时，阵列接收数据矢量的协方差矩阵的秩降为1，这就会导致信号子空间的维数小于信号源数，从而致使某些相干源的导向矢量与噪声子空间不完全正交，从而无法正确估计信号源方向。要实现对相干信号的有效跟踪，跟踪算法必须对接收信号矢量进行解相干。

发明内容

本发明的目的在于提出一种基于均匀线阵，应用于冲击噪声环境下的相干信号源的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）对接收到的信号数据矢量进行类归一化预处理：

z(t)为经过类归一化预处理的接收数据矢量，参数p的值由噪声的特征指数α确定；

（3）将均匀线阵划分为若干个子阵：

将阵元数为L的均匀线阵划分为K个子阵，每个子阵的阵元数为N，N=L-K+1，则第k个子阵的前向空间平滑的接收数据矢量为：

z_fk(t)＝[z_k(t),z_k+1(t),...,z_k+N-1(t)]^T

后向空间平滑的接收数据矢量为：

z_bk(t)＝[z_L-k+1(t),z_L-k(t),...,z_L-k-N+2(t)]^H

其中，k=1,2,…,K；

（4）将经过前后向空间平滑的接收数据矢量带入PAST算法，得到信号子空间：

参数y_fk(t)＝W^H(t-1)z_fk(t),y_bk(t)＝W^H(t-1)z_bk(t)，W(t)为第t个快拍时的权矩阵，更新参数：

$P_{\mathrm{fk}}\left(t\right)=P_{\mathrm{fk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)},$

$P_{\mathrm{bk}}\left(t\right)=P_{\mathrm{bk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)},$ β为遗忘系数，0＜β＜1，更新参数：

$g_{\mathrm{fk}}\left(t\right)=\frac{P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)},$ $g_{\mathrm{bk}}\left(t\right)=\frac{P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)},$ $h_{\mathrm{fk}}\left(t\right)=y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1),$

由先验信息得估计误差为e_fk(t)＝z_fk(t)-W(t-1)y_fk(t)，

e_bk(t)＝z_bk(t)-W(t-1)y_bk(t)，权矩阵更新为

$W\left(t\right)=W(t-1)+\frac{1}{2K}\{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{fk}}\left(t\right)g_{\mathrm{fk}}^H\left(t\right)+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{bk}}\left(t\right)g_{\mathrm{bk}}^H\left(t\right)\},$ 权矩阵W(t)即为快拍数为t时对应的信号特征矢量，由W(t)张成的空间为信号子空间，信号子空间U_S(t)＝orth[W(t)]；

（5）对信号子空间的数据使用MUSIC算法进行处理，经谱峰搜索获得目标信号波达方向：

信号子空间的MUSIC谱函数：

其中，I为单位阵，对MUSIC谱峰值进行搜索，得到入射信号的波达方向角度；

（6）令t=t+1，将下一块拍数的数据矢量带入步骤（1）中，得到下一块拍数对应的信号波达方向角度，直到达到最大快拍数，实现冲击噪声环境下相干信号源的DOA动态跟踪。

本发明的有益效果在于：

与现有的DOA跟踪方法相比，通过应用本发明的方法，采用类归一化预处理，能对冲击噪声进行有效地抑制；采用空间平滑的解相干方法，能对相干信号进行有效地分辨；采用投影子空间逼近跟踪方法，能对入射信号的信号子空间进行实时地跟踪，然后采用MUSIC方法对信号波达方向进行精测，该方法稳健性高，跟踪效果好。

附图说明

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

图2是本发明中子空间平滑拟合方法示意图。

图3是在信号为完全相干，阵元数为9，子阵数为6，广义信噪比为10dB，扫描间隔为0.1度，快拍数为1000，阵元间距为半波长，遗忘系数为0.99，特征指数为0.8时的DOA跟踪仿真图。

图4是在信号为完全相干，阵元数为9，子阵数为6，广义信噪比为10dB，扫描间隔为0.1度，快拍数为1000，阵元间距为半波长，遗忘系数为0.99，特征指数为1.5时的DOA跟踪仿真图。

具体实施方式

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

步骤（10）、对所述均匀线阵接收到的信号数据矢量，进行类归一化预处理，得到经过类归一化预处理的数据z(t)；

步骤（20）、将所述均匀线阵划分为若干个子阵，使用PAST算法对子阵的数据进行处理，得到快拍数为t时入射信号的信号子空间数据U_S(t)；

步骤（30）、对所述信号子空间的数据使用MUSIC算法进行处理，经谱峰搜索获得目标信号波达方向；

步骤（40）、输入下一快拍的信号数据矢量，转到步骤（10）。

其中，步骤（10）进一步包括：

步骤（110）、获取所述阵元接收的数据矢量，对所述接收数据矢量中的元素取绝对值，并找出所述数据绝对值中的最大值；

步骤（120）、根据噪声的冲击性，确定参数p；

步骤（130）、将所述阵元接收数据除以所述最大值的p次幂，得到的数据矢量Z(t)即为所述经过类归一化预处理的数据矢量。

其中，步骤（20）进一步包括：

步骤（210）、将所述均匀线阵划分为若干个子阵，用前一快拍权矩阵的共轭转置乘以各子阵的数据矢量，得到各子阵数据矢量的加权矢量y(t)；

步骤（220）、根据PAST算法原理得到当前快拍下的更新参数；

步骤（230）、用前一快拍的权矩阵乘以所述加权矢量，并用所述对应子阵数据矢量与所得矢量作差，得到估计误差；

步骤（240）、将所述各子阵的估计误差带入权矩阵的更新公式，得到当前快拍数的权矩阵W(t)；

步骤（250）、由所述当前快拍数的权矩阵得到入射信号的信号子空间数据。

其中，步骤（30）进一步包括：

步骤（310）、将所述信号子空间数据带入谱峰搜索公式，得到信号的空间谱；

步骤（320）、找出所述空间谱最大值所对应的角度，就得到信号的入射角度。

其中，步骤（220）中，更新参数由下式定义：

$P_{\mathrm{fk}}\left(t\right)=P_{\mathrm{fk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}$

$P_{\mathrm{bk}}\left(t\right)=P_{\mathrm{bk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}$

$g_{\mathrm{fk}}\left(t\right)=\frac{P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}$

$g_{\mathrm{bk}}\left(t\right)=\frac{P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}$

上式中，β为遗忘系数，0＜β＜1。y_fk(t)和y_bk(t)分别为快拍数等于t时的前向空间平滑加权矢量和后向空间平滑加权矢量。

其中，步骤（230）中，前后向空间平滑的估计误差的计算公式为

e_fk(t)＝z_fk(t)-W(t-1)y_fk(t)

e_bk(t)＝z_bk(t)-W(t-1)y_bk(t)

上式中，z_fk(t)和z_bk(t)分别为经过前向空间平滑后的数据矢量和经过后向空间平滑后的数据矢量。

其中，步骤（240）中，权矩阵的更新公式为

$W\left(t\right)=W(t-1)+\frac{1}{2K}\{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{fk}}\left(t\right)g_{\mathrm{fk}}^H\left(t\right)+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{bk}}\left(t\right)g_{\mathrm{bk}}^H\left(t\right)\}$

空间平滑算法是一种降低信号相关性的有效方法。其基本思想是通过将均匀线阵划分为几个互相重叠的子阵，各子阵的阵列流型都相同，各子阵的协方差矩阵可进行平均运算，从而实现去相关。前后向空间平滑是除了前向空间平滑外，对各子阵进行共轭反向重新构成一个子阵列，然后分别对经过前后向空间平滑处理的子阵的协方差矩阵进行平均处理，构成总的协方差矩阵，实现去相关。

传统的空间平滑算法需要对子阵的协方差矩阵进行平均运算，但是在PAST算法的迭代计算中，没有涉及到协方差矩阵的计算，所以，传统的空间平滑算法无法直接应用到PAST算法中。本发明将经过子阵划分后的数据矢量直接带入PAST算法的迭代计算过程，将空间平滑的思想和PAST算法相结合，在迭代计算过程中实现前后向的空间平滑。

如图1所示，本发明的方法首先对接收数据矢量进行类归一化预处理，之后将经过处理的数据矢量带入基于前后向空间平滑的PAST算法中，得到信号子空间，最后由MUSIC算法得到DOA估计。

本发明的阵列结构采用均匀线阵，阵元数为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维加性噪声矢量，噪声类型为复冲击噪声。此处的冲击噪声使用Alpha稳定分布来代表。

传统的DOA跟踪方法往往针对的是高斯噪声环境，由接收数据矢量的二阶矩或高阶矩来进行DOA估计。由于Alpha稳定分布没有二阶或高阶矩，因此在冲击噪声环境下传统算法失效。本发明采用对接收数据矢量类归一化预处理的方法，以降低冲击噪声对跟踪效果的影响。

令

$z\left(t\right)=\frac{X\left(t\right)}{{\left\{\max \right|x_i\left(t\right)\left|\right\}}^{(p-1)}}$

其中z(t)为经过类归一化预处理的接收数据矢量，参数p的值由噪声的特征指数α确定。由仿真可知，p的值随α的增大而减小。

对经过类归一化预处理的接收数据矢量采用前后向空间平滑。将阵元数为L的均匀线阵划分为K个子阵，每个子阵的阵元数为N,N=L-K+1，则第k个子阵的前向空间平滑的接收数据矢量可表示为

z_fk(t)＝[z_k(t),z_k+1(t),...,z_k+N-1(t)]^T

后向空间平滑的接收数据矢量可表示为

z_bk(t)＝[z_L-k+1(t),z_L-k(t),...,z_L-k-N+2(t)]^H

其中，k=1,2,…,K。

z_fk(t)为第k个子阵在第t个快拍时的前向平滑接收数据矢量，z_bk(t)为第k个子阵在第t个快拍时的后向平滑接收数据矢量。然后将经过空间平滑的数据带入PAST算法，此方法的示意图如图2所示。

定义:

y_fk(t)＝W^H(t-1)z_fk(t),y_bk(t)＝W^H(t-1)z_bk(t)

其中，W^H(t)为第t个快拍时的权矩阵。

定义更新参数P_fk(t)和P_bk(t):

$P_{\mathrm{fk}}\left(t\right)=P_{\mathrm{fk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}$

$P_{\mathrm{bk}}\left(t\right)=P_{\mathrm{bk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}$

其中，β为遗忘系数，0＜β＜1。

定义更新参数g_fk(t)，g_bk(t)，h_fk(t)和h_bk(t):

$g_{\mathrm{fk}}\left(t\right)=\frac{P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}$

$g_{\mathrm{bk}}\left(t\right)=\frac{P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}$

$h_{\mathrm{fk}}\left(t\right)=y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1)$

$h_{\mathrm{bk}}\left(t\right)=y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}^H(t-1)$

则，

$P_{\mathrm{fk}}\left(t\right)=P_{\mathrm{fk}}(t-1)-\frac{{\mathrm{\βH}}_{\mathrm{fk}}^H\left(t\right)h_{\mathrm{fk}}\left(t\right)}{1+\mathrm{\β}{y}_{\mathrm{fk}}^H\left(t\right)h_{\mathrm{fk}}^H\left(t\right)}$

$P_{\mathrm{bk}}\left(t\right)=P_{\mathrm{bk}}(t-1)-\frac{{\mathrm{\βH}}_{\mathrm{bk}}^H\left(t\right)h_{\mathrm{bk}}\left(t\right)}{1+\mathrm{\β}{y}_{\mathrm{bk}}^H\left(t\right)h_{\mathrm{bk}}^H\left(t\right)}$

由先验信息得到的估计误差为

e_fk(t)＝z_fk(t)-W(t-1)y_fk(t)

e_bk(t)＝z_bk(t)-W(t-1)y_bk(t)

因此，权矩阵更新公式为

$W\left(t\right)=W(t-1)+\frac{1}{2K}\{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{fk}}\left(t\right)g_{\mathrm{fk}}^H\left(t\right)+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{bk}}\left(t\right)g_{\mathrm{bk}}^H\left(t\right)\}$

W(t)即为快拍数为t时对应的信号特征矢量，由W(t)张成的空间即为信号子空间，W(t)W^H(t)即为信号子空间的投影矩阵。

由W(t)张成快拍数为t时刻的信号子空间U_S(t)，然后将信号子空间带入子空间类算法得到信号在快拍数为t时对应的入射角。

在本实施例中采用MUSIC谱估计方法，通过构造MUSIC谱函数来对信号的波达方向来进行估计。

求DOA可由最小化搜索来实现，即

$\widehat{\mathrm{\θ}}\left(t\right)=\arg \min a^H\left(\mathrm{\θ}\right)U_N\left(t\right)U_N^H\left(t\right)a\left(\mathrm{\θ}\right)$

$U_N\left(t\right)U_N^H=I-U_S\left(t\right)U_S^H\left(t\right)$

其中，I为单位阵，U_N(t)为噪声子空间。

构造信号子空间的MUSIC谱函数

$P_{\mathrm{MUSIC}}\left[\mathrm{\θ}\right]=\frac{1}{a^H\left(\mathrm{\θ}\right)[I-U_S\left(t\right)U_S^H\left(t\right)]a\left(\mathrm{\θ}\right)}$

然后对MUSIC谱峰值进行搜索，得到入射信号的波达方向角度。

当导向矢量属于信号子空间时，对于构造的空间谱函数P_MUSIC[θ]，分母是一个趋于0的值，P_MUSIC[θ]在此处有一尖峰。使θ的值变化，通过搜索谱峰值来估计快拍数为t时的信号波达方向角度。

令t=t+1,将下一块拍数的数据矢量带入所述算法中，得到下一块拍数对应的信号波达方向角度，以此来实现冲击噪声环境下相干信号源的DOA动态跟踪，跟踪效果如图3和图4所示。

与现有的DOA跟踪方法相比，采用本发明基于类归一化预处理的投影逼近子空间跟踪方法，能够有效地减少冲击噪声对DOA跟踪的影响，并且在噪声冲击性较强时有较好的跟踪性能。采用本发明将前后向空间平滑思想与投影逼近子空间跟踪算法相结合的跟踪方法，能够有效地降低接收信号之间的相关性，实现对相干信号的有效跟踪。由于本发明采用了投影逼近子空间跟踪方法，无需进行特征值分解计算，降低了计算量。综上，本发明能够对冲击噪声环境下的相干信号实现有效的跟踪，该方法的稳定性较好。

除了上述的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）对接收到的信号数据矢量进行类归一化预处理：

z(t)为经过类归一化预处理的接收数据矢量，参数p的值由噪声的特征指数α确定；

（3）将均匀线阵划分为若干个子阵：

将阵元数为L的均匀线阵划分为K个子阵，每个子阵的阵元数为N，N=L-K+1，则第k个子阵的前向空间平滑的接收数据矢量为：

z_fk(t)＝[z_k(t),z_k+1(t),...,z_k+N-1(t)]^T

后向空间平滑的接收数据矢量为：

z_bk(t)＝[z_L-k+1(t),z_L-k(t),...,z_L-k-N+2(t)]^H

其中，k=1,2,…,K；

（4）将经过前后向空间平滑的接收数据矢量带入PAST算法，得到信号子空间：

参数y_fk(t)＝W^H(t-1)z_fk(t),y_bk(t)＝W^H(t-1)z_bk(t)，W(t)为第t个快拍时的权矩阵，更新参数：

$P_{\mathrm{fk}}\left(t\right)=P_{\mathrm{fk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)},$

$P_{\mathrm{bk}}\left(t\right)=P_{\mathrm{bk}}(t-1)-\frac{\mathrm{\β}{P}_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}^H(t-1)}{1+{\mathrm{\βy}}_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)},$ β为遗忘系数，0＜β＜1，更新参数：

$g_{\mathrm{fk}}\left(t\right)=\frac{P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}(t-1)y_{\mathrm{fk}}\left(t\right)},$ $g_{\mathrm{bk}}\left(t\right)=\frac{P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)}{1-\mathrm{\β}+Y_{\mathrm{bk}}^H\left(t\right)P_{\mathrm{bk}}(t-1)y_{\mathrm{bk}}\left(t\right)},$ $h_{\mathrm{fk}}\left(t\right)=y_{\mathrm{fk}}^H\left(t\right)P_{\mathrm{fk}}^H(t-1),$

由先验信息得估计误差为e_fk(t)＝z_fk(t)-W(t-1)y_fk(t)，

e_bk(t)＝z_bk(t)-W(t-1)y_bk(t)，权矩阵更新为 $W\left(t\right)=W(t-1)+\frac{1}{2K}\{\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{fk}}\left(t\right)g_{\mathrm{fk}}^H\left(t\right)+\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{e}_{\mathrm{bk}}\left(t\right)g_{\mathrm{bk}}^H\left(t\right)\},$ 权矩阵W(t)即为快拍数为t时对应的信号特征矢量，由W(t)张成的空间为信号子空间，信号子空间U_S(t)＝orth[W(t)]；

（5）对信号子空间的数据使用MUSIC算法进行处理，经谱峰搜索获得目标信号波达方向：

信号子空间的MUSIC谱函数：

其中，I为单位阵，对MUSIC谱峰值进行搜索，得到入射信号的波达方向角度；

（6）令t=t+1，将下一块拍数的数据矢量带入步骤（1）中，得到下一块拍数对应的信号波达方向角度，直到达到最大快拍数，实现冲击噪声环境下相干信号源的DOA动态跟踪。

Images