背景介绍
在信号处理领域内,高斯分布假设一直占据着主导地位,因为其概率密度函数可以仅由均值和方差两个变量描述,这给信号的理论分析提供了极大的方便;而且基于高斯假设条件下的信号处理方法通常是线性的并且一般可以得到闭式最优解。然而,实际环境中存在许多非高斯特性的冲击噪声,如通信线路上的瞬间尖峰、由大气放电而引起的大气噪声以及雷达产生的部分杂波等。相比于高斯噪声,这类冲击噪声在时域上出现大量显著的尖峰脉冲特性,在概率密度上具有更加厚重的拖尾现象。Shao 和 Nikias 发现这种带有冲击性质的非高斯噪声可以用
α
\alpha
α稳定分布来描述。相对于高斯噪声,冲击噪声在实际中的用处更为广泛,因此关于冲击噪声背景下的研究也开始受到人们的关注,而DOA估计成为了研究的热点。
传统测向的方法较为突出的代表是MUSIC,ESPRIT算法等则在高斯噪声背景下表现出较好的性能,但是在冲击噪声背景下表现得往往不尽如人意,所以需要给出新的方法来实现DOA估计。
基于共变矩阵的ROC-MUSIC算法
采用阵元数目为
M
M
M的均匀线阵对
N
N
N个信号进行测向,则阵元接收的信号可以表示为
[
x
1
(
t
)
x
2
(
t
)
⋮
x
M
(
t
)
]
=
[
e
−
j
ω
0
τ
11
e
−
j
ω
0
τ
12
⋯
e
−
j
ω
0
τ
1
N
e
−
j
ω
0
τ
21
e
−
j
ω
0
τ
22
⋯
e
−
j
ω
0
τ
2
N
⋮
⋮
⋮
e
−
j
ω
0
τ
M
1
e
−
j
ω
0
τ
M
2
⋯
e
−
j
ω
0
τ
MN
]
[
s
1
(
t
)
s
2
(
t
)
⋮
s
N
(
t
)
]
+
[
n
1
(
t
)
n
2
(
t
)
⋮
n
M
(
t
)
]
\left[ \begin{matrix} {{x}_{1}}(t) \\ {{x}_{2}}(t) \\ \vdots \\ {{x}_{M}}(t) \\ \end{matrix} \right]=\left[ \begin{matrix} {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{11}}}} & {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{12}}}} & \cdots & {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{1N}}}} \\ {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{21}}}} & {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{22}}}} & \cdots & {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{2N}}}} \\ \vdots & \vdots & {} & \vdots \\ {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{M1}}}} & {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{M2}}}} & \cdots & {{\text{e}}^{-\text{j}{{\omega }_{0}}{{\tau }_{\text{MN}}}}} \\ \end{matrix} \right]\left[ \begin{matrix} {{s}_{1}}(t) \\ {{s}_{2}}(t) \\ \vdots \\ {{s}_{N}}(t) \\ \end{matrix} \right]+\left[ \begin{matrix} {{n}_{1}}(t) \\ {{n}_{2}}(t) \\ \vdots \\ {{n}_{M}}(t) \\ \end{matrix} \right]
⎣⎢⎢⎢⎡x1(t)x2(t)⋮xM(t)⎦⎥⎥⎥⎤=⎣⎢⎢⎢⎡e−jω0τ11e−jω0τ21⋮e−jω0τM1e−jω0τ12e−jω0τ22⋮e−jω0τM2⋯⋯⋯e−jω0τ1Ne−jω0τ2N⋮e−jω0τMN⎦⎥⎥⎥⎤⎣⎢⎢⎢⎡s1(t)s2(t)⋮sN(t)⎦⎥⎥⎥⎤+⎣⎢⎢⎢⎡n1(t)n2(t)⋮nM(t)⎦⎥⎥⎥⎤
将上式用矢量的形式表示则变为
X
(
t
)
=
A
(
θ
)
S
(
t
)
+
N
(
t
)
\mathbf{X}\left( t \right)=\mathbf{A}\left( \theta \right)\mathbf{S}\left( t \right)+\mathbf{N}\left( t \right)
X(t)=A(θ)S(t)+N(t)
其中,
X
(
t
)
\mathbf{X}\left( t \right)
X(t)表示阵列接收到的
M
×
1
M\times 1
M×1维的快拍数据,
A
(
θ
)
\mathbf{A}\left( \theta \right)
A(θ)表示
M
×
N
M\times N
M×N维的导向矢量,
S
(
t
)
\mathbf{S}\left( t \right)
S(t)表示
N
×
1
N\times 1
N×1维的信号矢量,
N
(
t
)
\mathbf{N}\left( t \right)
N(t)表示
M
×
1
M\times 1
M×1维的噪声矢量。
假设接收到的信号和噪声都是服从
S
α
S
S\alpha S
SαS分布的随机变量,彼此之间独立,非相干。信号和噪声的共变矩阵分别为
{
Γ
s
=
d
i
a
g
(
γ
s
1
,
γ
s
2
,
⋯
,
γ
s
N
)
Γ
N
=
d
i
a
g
(
γ
n
,
γ
n
,
⋯
,
γ
n
)
\left\{ \begin{aligned} & {{\mathbf{\Gamma }}_{s}}=diag\left( {{\gamma }_{{{s}_{1}}}},{{\gamma }_{{{s}_{2}}}},\cdots ,{{\gamma }_{{{s}_{N}}}} \right) \\ & {{\mathbf{\Gamma }}_{N}}=diag\left( {{\gamma }_{n}},{{\gamma }_{n}},\cdots ,{{\gamma }_{n}} \right) \\ \end{aligned} \right.
{Γs=diag(γs1,γs2,⋯,γsN)ΓN=diag(γn,γn,⋯,γn)
阵列的共变矩阵为
[
Γ
i
k
]
=
[
x
i
(
t
)
,
x
k
(
t
)
]
α
=
∑
j
=
1
N
a
i
(
θ
j
)
a
k
<
α
−
1
>
(
θ
j
)
γ
s
j
+
γ
n
δ
i
k
i
,
k
=
1
,
2
⋯
,
M
\left[ {{\mathbf{\Gamma }}_{ik}} \right]={{\left[ {{x}_{i}}(t),{{x}_{k}}(t) \right]}_{\alpha }}=\sum\limits_{j=1}^{N}{{{a}_{i}}}\left( {{\theta }_{j}} \right)a_{k}^{<\alpha -1>}\left( {{\theta }_{j}} \right){{\gamma }_{{{s}_{j}}}}+{{\gamma }_{n}}{{\delta }_{ik}}\quad i,k=1,2\cdots ,M
[Γik]=[xi(t),xk(t)]α=j=1∑Nai(θj)ak<α−1>(θj)γsj+γnδiki,k=1,2⋯,M
矢量表示形式为
Γ
X
=
A
Γ
s
A
H
+
γ
n
I
{{\mathbf{\Gamma }}_{X}}=\mathbf{A}{{\mathbf{\Gamma }}_{s}}{{\mathbf{A}}^{H}}+{{\gamma }_{n}}\mathbf{I}
ΓX=AΓsAH+γnI
根据共变矩阵的定义可以求出
Γ
i
k
=
E
{
x
i
(
t
)
∣
x
k
(
t
)
∣
p
−
2
x
k
∗
(
t
)
}
E
{
∣
x
k
(
t
)
∣
p
}
1
<
p
≤
2
{{\mathbf{\Gamma }}_{ik}}\text{=}\frac{E\left\{ {{x}_{i}}(t){{\left| {{x}_{k}}(t) \right|}^{p-2}}x_{k}^{*}(t) \right\}}{E\left\{ {{\left| {{x}_{k}}(t) \right|}^{p}} \right\}}\quad 1<p\le 2
Γik=E{∣xk(t)∣p}E{xi(t)∣xk(t)∣p−2xk∗(t)}1<p≤2
通过上式可以看出,矩阵
Γ
\mathbf{\Gamma }
Γ的
M
−
N
M-N
M−N个较小的特征值对应
γ
n
{{\gamma }_{n}}
γn,其特征矢量张成的噪声子空间与矩阵
A
\mathbf{A}
A的列矢量正交,信号子空间则是由
N
N
N个较大特征值对应的特征矢量张成的空间组成。因此ROC-MUSIC算法的流程可以简述为以下过程
1)计算
M
×
M
M\times M
M×M维的矩阵
Γ
\mathbf{\Gamma }
Γ,其元素为
Γ
i
k
=
∑
t
=
1
K
x
i
(
t
)
∣
x
k
(
t
)
∣
p
−
2
x
k
∗
(
t
)
∑
t
=
1
K
∣
x
k
(
t
)
∣
p
1
<
p
<
α
≤
2
{{\mathbf{\Gamma }}_{ik}}\text{=}\frac{\sum\limits_{t=1}^{K}{{{x}_{i}}(t){{\left| {{x}_{k}}(t) \right|}^{p-2}}x_{k}^{*}(t)}}{\sum\limits_{t=1}^{K}{{{\left| {{x}_{k}}(t) \right|}^{p}}}}\quad 1<p<\alpha \le 2
Γik=t=1∑K∣xk(t)∣pt=1∑Kxi(t)∣xk(t)∣p−2xk∗(t)1<p<α≤2
2)对矩阵
Γ
\mathbf{\Gamma }
Γ进行奇异值分解,构成
M
×
(
M
−
N
)
M\times \left( M-N \right)
M×(M−N)矩阵
E
n
=
(
e
^
N
+
1
,
e
^
N
+
2
,
⋯
,
e
^
M
)
{{E}_{n}}=\left( {{{\hat{e}}}_{N+1}},{{{\hat{e}}}_{N+2}},\cdots ,{{{\hat{e}}}_{M}} \right)
En=(e^N+1,e^N+2,⋯,e^M) ,其中
e
^
N
+
1
,
e
^
N
+
2
,
⋯
,
e
^
M
{{\hat{e}}_{N+1}},{{\hat{e}}_{N+2}},\cdots ,{{\hat{e}}_{M}}
e^N+1,e^N+2,⋯,e^M是
Γ
\mathbf{\Gamma }
Γ中最小的
M
−
N
M-N
M−N个奇异值对应的奇异向量
3)计算谱
S
(
θ
)
=
1
a
H
(
θ
)
E
n
E
n
H
a
(
θ
)
S(\theta )=\frac{1}{{{\mathbf{a}}^{\text{H}}}(\theta ){{\mathbf{E}}_{n}}\mathbf{E}_{n}^{\text{H}}\mathbf{a}(\theta )}\quad
S(θ)=aH(θ)EnEnHa(θ)1
其中
a
(
θ
)
=
[
1
,
e
−
j
2
π
d
2
λ
sin
θ
,
e
−
j
2
π
d
3
λ
sin
θ
,
⋯
,
e
−
j
2
π
d
M
λ
sin
θ
]
T
\mathbf{a}(\theta )={{\left[ 1,{{\text{e}}^{-\text{j}2\pi \frac{{{d}_{2}}}{\lambda }\sin \theta }},{{\text{e}}^{-\text{j}2\pi \frac{{{d}_{3}}}{\lambda }\sin \theta }},\cdots ,{{\text{e}}^{-\text{j}2\pi \frac{{{d}_{M}}}{\lambda }\sin \theta }} \right]}^{\text{T}}}
a(θ)=[1,e−j2πλd2sinθ,e−j2πλd3sinθ,⋯,e−j2πλdMsinθ]T
4)
S
(
θ
)
S\left( \theta \right)
S(θ)中
N
N
N个局部极大值即为波达方向角
{
θ
1
,
θ
2
,
⋯
,
θ
N
}
\left\{ {{\theta }_{1}},{{\theta }_{2}},\cdots ,{{\theta }_{N}} \right\}
{θ1,θ2,⋯,θN}的估计。
仿真实验
仿真实验所用到的信号和噪声分别是远场窄带信号和冲击噪声,采用阵元数为10的均匀线阵,阵元间距为半波长,三个独立的入射信号,入射角度分别为-5°,15°和35°,广义信噪比为3dB,具体参数如下:
仿真参数 | 参数值 |
---|---|
阵元数目 | 10 |
信源数目 | 3 |
阵元间距 | 半波长 |
快拍数 | 1000 |
α \alpha α | 1.4 |
γ \gamma γ | 1 |
对于服从对称
α
\alpha
α稳定分布的冲击噪声来说,它的分布特性由特征指数
α
\alpha
α和尺度参数
γ
\gamma
γ 决定,但是当
α
∈
(
0
,
2
]
\alpha \in (0,2]
α∈(0,2]时,其方差是无界的,因而可以用信号的平均功率和冲击噪声的尺度参数的比值即广义信噪比(GSNR)代替一般的信噪比,其表达式如下
GSNR
=
10
lg
(
1
γ
K
∑
t
=
1
K
∣
s
(
t
)
2
∣
)
\text{GSNR}=10\lg \left( \frac{1}{\gamma K}\sum\limits_{t=1}^{K}{\left| s{{\left( t \right)}^{2}} \right|} \right)
GSNR=10lg(γK1t=1∑K∣∣∣s(t)2∣∣∣)
其中,
K
K
K为快拍数。当
α
=
2
\alpha \text{=}2
α=2时,上式退化为一般的信噪比。
冲击噪声的具体参考:冲击噪声介绍
实验结果如下:
冲击噪声背景下,传统MUSIC算法不能对入射角度进行准确的估计,这是因为冲击噪声会对入射信号进行较大的干扰,从而对结果造成影响甚至是误判,而经过对噪声进行重新建模分析的ROC-MUSIC算法很好地区分出了入射信号的角度。
令
α
=
2
\alpha=2
α=2,此时的冲击噪声退化为高斯噪声,相同仿真条件下得到如下结果
当冲击噪声退化成高斯噪声时,无论是传统的MUSIC算法还是ROC-MUSIC算法能较好地估计出入射信号的角度,但是此时传统MUSIC算法表现出较好性能,这也和前面的理论分析一致。
参考文献:
[1]P. Tsakalides and C. L. Nikias, “The robust covariation-based MUSIC (ROC-MUSIC) algorithm for bearing estimation in impulsive noise environments,” in IEEE Transactions on Signal Processing, vol. 44, no. 7, pp. 1623-1633, July 1996, doi: 10.1109/78.510611.
[2]赵大勇,刁鸣,杨丽丽,陈超.冲击噪声背景下的动态DOA跟踪[J].山东大学学报(工学版),2010,40(01):133-138.
[3]李洪升,张瑞丰,杜宇,王永孝.冲击噪声环境中波束域DOA估计方法研究[J].现代防御技术,2017,45(03):93-97.