SPICE(协方差稀疏迭代估)算法介绍
假设有
K
K
K个远场窄带信号,入射到由
M
(
M
>
K
)
M\left( M>K \right)
M(M>K)个阵元组成的均匀线上,同时信号和噪声互不相关,令
Ω
\Omega
Ω表示信源可能的方向,
{
θ
k
=
1
Q
}
\left\{ \theta _{k=1}^{Q} \right\}
{θk=1Q}是包含
Ω
\Omega
Ω 的集合,并且
Q
≫
K
Q\gg K
Q≫K,即真实的信源方向包含在
Q
Q
Q 中,那么阵列的输出信号可以表示为
Y
(
t
)
=
∑
k
=
1
Q
a
(
θ
k
)
s
k
(
t
)
+
e
(
t
)
t
=
1
,
2
,
⋯
,
L
\mathbf{Y}\left( t \right)=\sum\limits_{k=1}^{Q}{\mathbf{a}\left( {{\theta }_{k}} \right){{s}_{k}}\left( t \right)+\mathbf{e}\left( t \right)}\ \ \ \ \ \ \ t=1,2,\cdots ,L
Y(t)=k=1∑Qa(θk)sk(t)+e(t) t=1,2,⋯,L
其中
L
L
L表示快拍数大小,
Y
(
t
)
∈
C
M
×
1
\mathbf{Y}\left( t \right)\in {{C}^{M\times 1}}
Y(t)∈CM×1表示第
t
t
t个快拍时的观察数据,
a
(
θ
k
)
∈
C
M
×
1
\mathbf{a}\left( {{\theta }_{k}} \right)\in {{C}^{M\times 1}}
a(θk)∈CM×1 为入射方向
θ
k
{{\theta }_{k}}
θk 对应的导向矢量,
s
k
(
t
)
∈
C
{{s}_{k}}\left( t \right)\in C
sk(t)∈C 表示信源可能方向
θ
k
{{\theta }_{k}}
θk 的入射信号,
e
(
t
)
∈
C
M
×
1
\mathbf{e}\left( t \right)\in {{C}^{M\times 1}}
e(t)∈CM×1是零均值的加性高斯噪声。
{
a
(
θ
k
)
}
k
=
1
Q
\left\{ \mathbf{a}\left( {{\theta }_{k}} \right) \right\}_{k=1}^{Q}
{a(θk)}k=1Q表示阵列的过完备基,可以表示为
a
(
θ
k
)
=
[
1
,
e
−
j
2
π
d
sin
θ
k
/
λ
,
⋯
,
e
−
j
2
π
(
M
−
1
)
d
sin
θ
k
/
λ
]
\mathbf{a}\left( {{\theta }_{k}} \right)=\left[ 1,{{e}^{-j2\pi d\sin {{\theta }_{k}}/\lambda }},\cdots ,{{e}^{-j2\pi \left( M-1 \right)d\sin {{\theta }_{k}}/\lambda }} \right]
a(θk)=[1,e−j2πdsinθk/λ,⋯,e−j2π(M−1)dsinθk/λ]
其中
λ
\lambda
λ表示入射信号的波长,
d
d
d 表示阵元之间的间距大小。
那么,阵列的协方差矩阵可以表示为
R
a
=
E
[
Y
(
t
)
Y
H
(
t
)
]
=
B
(
θ
)
R
s
B
H
(
θ
)
+
e
(
t
)
{{R}_{a}}=E\left[ \mathbf{Y}\left( t \right){{\mathbf{Y}}^{H}}\left( t \right) \right]=\mathbf{B}\left( \theta \right){{R}_{s}}{{\mathbf{B}}^{H}}\left( \theta \right)+\mathbf{e}\left( t \right)
Ra=E[Y(t)YH(t)]=B(θ)RsBH(θ)+e(t)
其中
R
s
=
E
[
S
(
t
)
S
H
(
t
)
]
{{R}_{s}}\text{=}E\left[ \mathbf{S}\left( t \right){{\mathbf{S}}^{H}}\left( t \right) \right]
Rs=E[S(t)SH(t)]表示信号的协方差矩阵,
S
(
t
)
=
[
s
1
(
t
)
,
s
2
(
t
)
,
⋯
,
s
Q
(
t
)
]
T
\mathbf{S}\left( t \right)={{\left[ {{s}_{1}}\left( t \right),{{s}_{2}}\left( t \right),\cdots ,{{s}_{Q}}\left( t \right) \right]}^{T}}
S(t)=[s1(t),s2(t),⋯,sQ(t)]T ,
B
(
θ
)
\mathbf{B}\left( \theta \right)
B(θ) 是
M
×
Q
M\times Q
M×Q的阵列流形矩阵。
对噪声进行求期望可以得到
E
[
e
(
t
)
e
∗
(
t
ˉ
)
]
=
[
σ
1
0
⋯
0
0
σ
2
⋯
0
⋮
⋮
⋱
⋮
0
⋯
⋯
σ
M
]
δ
t
,
t
ˉ
E\left[ \mathbf{e}\left( t \right){{\mathbf{e}}^{*}}\left( {\bar{t}} \right) \right]\text{=}\left[ \begin{matrix} {{\sigma }_{1}} & 0 & \cdots & 0 \\ 0 & {{\sigma }_{2}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & {{\sigma }_{M}} \\ \end{matrix} \right]{{\delta }_{t,\bar{t}}}
E[e(t)e∗(tˉ)]=⎣⎢⎢⎢⎡σ10⋮00σ2⋮⋯⋯⋯⋱⋯00⋮σM⎦⎥⎥⎥⎤δt,tˉ
其中,
δ
t
,
t
ˉ
=
{
1
t
=
t
ˉ
0
t
≠
t
ˉ
{{\delta }_{t,\bar{t}}}\text{=}\left\{ \begin{aligned} & 1\ \ \ \ \ \ \ t=\bar{t} \\ & 0\ \ \ \ \ \ t\ne \bar{t} \\ \end{aligned} \right.
δt,tˉ={1 t=tˉ0 t=tˉ
协方差矩阵可以写成
R
a
=
E
[
Y
(
t
)
Y
(
t
)
H
]
=
∑
k
=
1
Q
p
k
a
k
a
k
∗
+
[
σ
1
0
⋯
0
0
σ
2
⋯
0
⋮
⋮
⋱
⋮
0
⋯
⋯
σ
M
]
=
[
a
1
,
⋯
,
a
Q
,
I
]
[
p
1
0
⋯
⋯
⋯
⋯
0
0
p
2
0
⋯
⋯
⋯
0
⋮
0
⋱
⋮
⋮
⋮
⋮
0
⋯
⋯
p
Q
⋯
⋯
0
0
⋯
⋯
⋯
σ
1
⋯
0
⋮
⋮
⋮
⋮
⋮
⋱
⋮
0
⋯
⋯
⋯
⋯
⋯
σ
M
]
[
a
1
∗
⋮
a
Q
∗
I
]
≜
A
P
A
H
\begin{aligned} & {{R}_{a}}=E\left[ \mathbf{Y}(t)\mathbf{Y}{{(t)}^{\text{H}}} \right]=\sum\limits_{k=1}^{Q}{{{p}_{k}}}{{\mathbf{a}}_{k}}\mathbf{a}_{k}^{*}+\left[ \begin{matrix} {{\sigma }_{1}} & 0 & \cdots & 0 \\ 0 & {{\sigma }_{2}} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & {{\sigma }_{M}} \\ \end{matrix} \right] \\ & =\left[ {{\mathbf{a}}_{1}},\cdots ,{{\mathbf{a}}_{Q}},\mathbf{I} \right]\left[ \begin{matrix} {{p}_{1}} & 0 & \cdots & \cdots & \cdots & \cdots & 0 \\ 0 & {{p}_{2}} & 0 & \cdots & \cdots & \cdots & 0 \\ \vdots & 0 & \ddots & \vdots & \vdots & \vdots & \vdots \\ 0 & \cdots & \cdots & {{p}_{Q}} & \cdots & \cdots & 0 \\ 0 & \cdots & \cdots & \cdots & {{\sigma }_{1}} & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & \cdots & \cdots & \cdots & \cdots & \cdots & {{\sigma }_{M}} \\ \end{matrix} \right]\left[ \begin{matrix} {{\mathbf{a}}_{1}}^{*} \\ \vdots \\ {{\mathbf{a}}_{Q}}^{*} \\ \mathbf{I} \\ \end{matrix} \right] \\ & \triangleq AP{{A}^{\text{H}}} \end{aligned}
Ra=E[Y(t)Y(t)H]=k=1∑Qpkakak∗+⎣⎢⎢⎢⎡σ10⋮00σ2⋮⋯⋯⋯⋱⋯00⋮σM⎦⎥⎥⎥⎤=[a1,⋯,aQ,I]⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡p10⋮00⋮00p20⋯⋯⋮⋯⋯0⋱⋯⋯⋮⋯⋯⋯⋮pQ⋯⋮⋯⋯⋯⋮⋯σ1⋮⋯⋯⋯⋮⋯⋯⋱⋯00⋮00⋮σM⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎡a1∗⋮aQ∗I⎦⎥⎥⎥⎤≜APAH
为了估计
{
p
k
}
k
=
1
Q
+
M
\left\{ {{p}_{k}} \right\}_{k=1}^{Q+M}
{pk}k=1Q+M ,可以通过以下得到
f
=
∥
R
a
−
1
/
2
(
R
^
a
−
R
a
)
R
^
a
−
1
/
2
∥
2
f=\ {{\left\| {{R}_{a}}^{-1/2}\left( {{{\hat{R}}}_{a}}-{{R}_{a}} \right){{{\hat{R}}}_{a}}^{-1/2} \right\|}^{2}}
f= ∥∥∥Ra−1/2(R^a−Ra)R^a−1/2∥∥∥2
同时为了求解上式,通过简单的运算,可以将其转化为
f
=
tr
[
R
a
−
1
(
R
^
a
−
R
a
)
R
^
a
−
1
(
R
^
a
−
R
a
)
]
=
tr
[
(
R
a
−
1
R
^
a
−
I
)
(
I
−
R
^
a
−
1
R
a
)
]
=
tr
(
R
a
−
1
R
^
a
)
+
tr
(
R
^
a
−
1
R
a
)
−
2
M
\begin{aligned} & f=\operatorname{tr}\left[ {{R}_{a}}^{-1}({{{\hat{R}}}_{a}}-{{R}_{a}}){{{\hat{R}}}_{a}}^{-1}({{{\hat{R}}}_{a}}-{{R}_{a}}) \right] \\ & =\operatorname{tr}\left[ \left( {{R}_{a}}^{-1}{{{\hat{R}}}_{a}}-\mathbf{I} \right)\left( \mathbf{I}-{{{\hat{R}}}_{a}}^{-1}{{R}_{a}} \right) \right] \\ & =\operatorname{tr}\left( {{R}_{a}}^{-1}{{{\hat{R}}}_{a}} \right)+\operatorname{tr}\left( {{{\hat{R}}}_{a}}^{-1}{{R}_{a}} \right)-2M \end{aligned}
f=tr[Ra−1(R^a−Ra)R^a−1(R^a−Ra)]=tr[(Ra−1R^a−I)(I−R^a−1Ra)]=tr(Ra−1R^a)+tr(R^a−1Ra)−2M
其中
R
a
−
1
/
2
{{R}_{a}}^{-1/2}
Ra−1/2 表示正定矩阵
R
a
−
1
{{R}_{a}}^{-1}
Ra−1的平方根,
R
^
a
=
Y
Y
H
/
M
{{\hat{R}}_{a}}\text{=}\mathbf{Y}{{\mathbf{Y}}^{H}}/M
R^a=YYH/M ,等价于最小化函数
g
=
tr
(
R
^
a
1
/
2
R
−
1
R
^
a
1
/
2
)
+
∑
k
=
1
Q
+
M
(
a
k
∗
R
^
a
−
1
a
k
)
p
k
g=\operatorname{tr}\left( {{{\hat{R}}}_{a}}^{1/2}{{R}^{-1}}{{{\hat{R}}}_{a}}^{1/2} \right)+\sum\limits_{k=1}^{Q+M}{\left( a_{k}^{*}{{{\hat{R}}}_{a}}^{-1}{{a}_{k}} \right)}{{p}_{k}}
g=tr(R^a1/2R−1R^a1/2)+k=1∑Q+M(ak∗R^a−1ak)pk
为了求解上式可以将其转化为带有限制的最小化问题
{
min
p
k
≥
0
t
r
(
R
^
a
1
/
2
R
−
1
R
^
a
1
/
2
)
s
.
t
∑
k
=
1
Q
+
M
ω
k
p
k
=
1
\left\{ \begin{aligned} & \underset{{{p}_{k}}\ge 0}{\mathop{\min }}\,\ \ \ \ \ tr\left( {{{\hat{R}}}_{a}}^{1/2}{{R}^{-1}}{{{\hat{R}}}_{a}}^{1/2} \right) \\ & s.t\ \ \ \ \ \ \ \sum\limits_{k=1}^{Q+M}{{{\omega }_{k}}{{p}_{k}}=1} \\ \end{aligned} \right.
⎩⎪⎪⎪⎨⎪⎪⎪⎧pk≥0min tr(R^a1/2R−1R^a1/2)s.t k=1∑Q+Mωkpk=1
其中
ω
k
=
a
k
∗
R
a
−
1
a
k
/
M
{{\omega }_{k}}={{\mathbf{a}}_{k}}^{*}R_{a}^{-1}{{\mathbf{a}}_{k}}/M
ωk=ak∗Ra−1ak/M 。
求解约束条件下的方程就可以得到最终的解,可以利用迭代算法来进行求解上述方程,这就是SPICE算法。
因此SPICE算法可以总结为
1) 通过时间平均来初始化功率谱估计
p
k
0
=
a
k
∗
R
^
a
a
k
/
∥
a
k
∥
4
k
=
1
,
2
,
⋯
,
Q
+
M
p_{k}^{0}={{\mathbf{a}}_{k}}^{*}{{\hat{R}}_{a}}{{\mathbf{a}}_{k}}/{{\left\| {{\mathbf{a}}_{k}} \right\|}^{4}}\ \ \ \ \ \ k=1,2,\cdots ,Q+M
pk0=ak∗R^aak/∥ak∥4 k=1,2,⋯,Q+M
2) 重复下式迭代过程直至收敛
{
p
k
i
+
1
=
p
k
i
∥
a
k
R
a
−
1
(
i
)
R
^
a
1
/
2
∥
ω
k
1
/
2
ρ
(
i
)
k
=
1
,
2
,
⋯
,
Q
+
M
ρ
(
i
)
=
∑
m
=
1
Q
+
M
ω
m
1
/
2
p
m
i
∥
a
m
R
a
−
1
(
i
)
R
^
a
1
/
2
∥
\left\{ \begin{aligned} & p_{k}^{i+1}=p_{k}^{i}\frac{\left\| {{\mathbf{a}}_{k}}{{R}_{a}}^{-1}\left( i \right){{\hat{R}}_{a}}^{1/2} \right\|}{\omega _{k}^{1/2}\rho \left( i \right)}\ \ \ \ \ \ \ \ k=1,2,\cdots ,Q+M \\ & \rho \left( i \right)=\sum\limits_{m=1}^{Q+M}{\omega _{m}^{1/2}p_{m}^{i}}\left\| {{\mathbf{a}}_{m}}{{R}_{a}}^{-1}\left( i \right){{\hat{R}}_{a}}^{1/2} \right\| \\ \end{aligned} \right.
⎩⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎧pki+1=pkiωk1/2ρ(i)∥∥∥akRa−1(i)R^a1/2∥∥∥ k=1,2,⋯,Q+Mρ(i)=m=1∑Q+Mωm1/2pmi∥∥∥amRa−1(i)R^a1/2∥∥∥
其中
R
a
(
i
)
{{R}_{a}}\left( i \right)
Ra(i) 上述计算公式重新得到的协方差矩阵。
基于协方差稀疏迭代估计法的波达方向估计,是一种计算效率很高的稀疏迭代算法,该方法以最小化协方差矩阵拟合误差准则得到的迭代算法,在多快拍的情况下具有很好的效果,在少快拍的情况下也具有很好的效果。与其他稀疏算法相比,概算法具有以下特点:
简单且完美的统计基础
自然地考虑了数据中的噪声
不需要设置一些超参数
具有很好的收敛性
基于四阶累积量的SPICE算法
根据四阶累积量的定义,可以得到阵列输出信号的四阶累积量矩阵,如下
C
y
=
D
(
θ
)
C
s
D
(
θ
)
H
+
C
e
{{C}_{y}}=D\left( \theta \right){{C}_{s}}D{{\left( \theta \right)}^{H}}+{{C}_{e}}
Cy=D(θ)CsD(θ)H+Ce
其中
C
y
{{C}_{y}}
Cy 为输出信号的四阶累积量矩阵,
C
s
{{C}_{s}}
Cs 为信源的四阶累积量矩阵,
C
e
{{C}_{e}}
Ce 为噪声的四阶累积量矩阵,
D
(
θ
)
=
[
a
(
θ
1
)
⊗
a
∗
(
θ
1
)
,
⋯
,
a
(
θ
Q
)
⊗
a
∗
(
θ
Q
)
]
=
[
d
1
,
⋯
,
d
Q
]
D\left( \theta \right)=\left[ \mathbf{a}\left( {{\theta }_{1}} \right)\otimes {{\mathbf{a}}^{*}}\left( {{\theta }_{1}} \right),\cdots ,\mathbf{a}\left( {{\theta }_{Q}} \right)\otimes {{\mathbf{a}}^{*}}\left( {{\theta }_{Q}} \right) \right]\text{=}\left[ {{\mathbf{d}}_{1}},\cdots ,{{\mathbf{d}}_{Q}} \right]
D(θ)=[a(θ1)⊗a∗(θ1),⋯,a(θQ)⊗a∗(θQ)]=[d1,⋯,dQ] 为阵列流形矩阵。
通过将SPICE算法中的
R
a
{{R}_{a}}
Ra 和
R
a
{{R}_{a}}
Ra替换成
C
y
{{C}_{y}}
Cy和
{
d
k
}
\left\{ {{\mathbf{d}}_{k}} \right\}
{dk} 就可以得到基于四阶累积量的SPICE算法。但是这种方法的计算量过大,因为四阶累积量矩阵中包含了大量的冗余数据,为了简化计算量,令
X
(
t
)
=
[
y
1
(
t
)
,
y
M
(
t
)
]
T
X\left( t \right)={{\left[ {{\mathbf{y}}_{1}}\left( t \right),{{\mathbf{y}}_{M}}\left( t \right) \right]}^{T}}
X(t)=[y1(t),yM(t)]T ,可以得到
X
(
t
)
=
B
ˉ
(
θ
)
S
(
t
)
+
e
(
t
)
\mathbf{X}\left( t \right)=\mathbf{\bar{B}}\left( \theta \right)\mathbf{S}\left( t \right)+\mathbf{e}\left( t \right)
X(t)=Bˉ(θ)S(t)+e(t)
式中,
B
ˉ
(
θ
)
=
[
a
ˉ
(
θ
1
)
,
⋯
,
a
ˉ
(
θ
Q
)
]
\mathbf{\bar{B}}\left( \theta \right)\text{=}\left[ \mathbf{\bar{a}}\left( {{\theta }_{1}} \right),\cdots ,\mathbf{\bar{a}}\left( {{\theta }_{Q}} \right) \right]
Bˉ(θ)=[aˉ(θ1),⋯,aˉ(θQ)] ,
a
ˉ
(
θ
k
)
=
[
1
,
e
−
j
2
π
(
M
−
1
)
d
sin
θ
k
/
λ
]
T
\mathbf{\bar{a}}\left( {{\theta }_{k}} \right)={{\left[ 1,{{e}^{-j2\pi \left( M-1 \right)d\sin {{\theta }_{k}}/\lambda }} \right]}^{T}}
aˉ(θk)=[1,e−j2π(M−1)dsinθk/λ]T 。
那么四阶累积量矩阵可以变化为
C
ˉ
y
=
D
ˉ
(
θ
)
C
ˉ
s
D
ˉ
(
θ
)
H
+
C
ˉ
e
{{\bar{C}}_{y}}=\bar{D}\left( \theta \right){{\bar{C}}_{s}}\bar{D}{{\left( \theta \right)}^{H}}+{{\bar{C}}_{e}}
Cˉy=Dˉ(θ)CˉsDˉ(θ)H+Cˉe
其中,
D
ˉ
(
θ
)
=
[
a
ˉ
(
θ
1
)
⊗
a
∗
(
θ
1
)
,
⋯
,
a
ˉ
(
θ
Q
)
⊗
a
∗
(
θ
Q
)
]
=
[
d
ˉ
1
,
⋯
,
d
ˉ
Q
]
\bar{D}\left( \theta \right)=\left[ \mathbf{\bar{a}}\left( {{\theta }_{1}} \right)\otimes {{\mathbf{a}}^{*}}\left( {{\theta }_{1}} \right),\cdots ,\mathbf{\bar{a}}\left( {{\theta }_{Q}} \right)\otimes {{\mathbf{a}}^{*}}\left( {{\theta }_{Q}} \right) \right]\text{=}\left[ {{{\mathbf{\bar{d}}}}_{1}},\cdots ,{{{\mathbf{\bar{d}}}}_{Q}} \right]
Dˉ(θ)=[aˉ(θ1)⊗a∗(θ1),⋯,aˉ(θQ)⊗a∗(θQ)]=[dˉ1,⋯,dˉQ] 表示导向矢量,
C
ˉ
y
{{\bar{C}}_{y}}
Cˉy 是
2
M
×
2
M
2M\times 2M
2M×2M的四阶累积量矩阵。
将SPICE算法中的
R
a
{{R}_{a}}
Ra和
{
a
k
}
\left\{ {{\mathbf{a}}_{k}} \right\}
{ak} 替换成
C
ˉ
y
{{\bar{C}}_{y}}
Cˉy 和
{
d
ˉ
k
}
\left\{ {{{\mathbf{\bar{d}}}}_{k}} \right\}
{dˉk} 就可以得到基于四阶累积量的SPICE算法,其步骤如下:
1)通过时间平均来初始化功率谱估计
p
k
0
=
d
ˉ
k
∗
C
ˉ
^
y
d
ˉ
k
/
∥
d
ˉ
k
∥
4
k
=
1
,
2
,
⋯
,
Q
+
M
2
p_{k}^{0}={{\mathbf{\bar{d}}}_{k}}^{*}{{\hat{\bar{C}}}_{y}}{{\mathbf{\bar{d}}}_{k}}/{{\left\| {{{\mathbf{\bar{d}}}}_{k}} \right\|}^{4}}\ \ \ \ \ \ k=1,2,\cdots ,Q+{{M}^{2}}
pk0=dˉk∗Cˉ^ydˉk/∥∥dˉk∥∥4 k=1,2,⋯,Q+M2
2)重复下式迭代过程直至收敛
{
p
k
i
+
1
=
p
k
i
∥
d
ˉ
k
k
C
ˉ
y
−
1
(
i
)
C
ˉ
^
y
1
/
2
∥
ω
k
1
/
2
ρ
(
i
)
k
=
1
,
2
,
⋯
,
Q
+
M
2
ρ
(
i
)
=
∑
m
=
1
Q
+
M
2
ω
m
1
/
2
p
m
i
∥
d
ˉ
m
C
ˉ
y
−
1
(
i
)
C
ˉ
^
y
1
/
2
∥
\left\{ \begin{aligned} & p_{k}^{i+1}=p_{k}^{i}\frac{\left\| {{{\mathbf{\bar{d}}}}_{k}}_{k}{{{\bar{C}}}_{y}}^{-1}\left( i \right){{{\hat{\bar{C}}}}_{y}}^{1/2} \right\|}{\omega _{k}^{1/2}\rho \left( i \right)}\ \ \ \ \ \ \ \ k=1,2,\cdots ,Q+{{M}^{2}} \\ & \rho \left( i \right)=\sum\limits_{m=1}^{Q+{{M}^{2}}}{\omega _{m}^{1/2}p_{m}^{i}}\left\| {{{\mathbf{\bar{d}}}}_{m}}{{{\bar{C}}}_{y}}^{-1}\left( i \right){{\hat{\bar{C}}}_{y}}^{1/2} \right\| \\ \end{aligned} \right.
⎩⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎧pki+1=pkiωk1/2ρ(i)∥∥∥∥dˉkkCˉy−1(i)Cˉ^y1/2∥∥∥∥ k=1,2,⋯,Q+M2ρ(i)=m=1∑Q+M2ωm1/2pmi∥∥∥∥dˉmCˉy−1(i)Cˉ^y1/2∥∥∥∥
基于四阶累积量的协方差稀疏迭代估计算法,不仅提高了原本的超分辨性能,而且在低信噪比的条件下也具有很好的性能,同时通过对均匀线性阵列非冗余数据的提取,降低了计算复杂度。
实验参数设置
参数名称 | 参数值 |
---|---|
阵元数 | 10 |
阵元间距 | λ / 2 \lambda/2 λ/2 |
信源数目 | 3 |
信源角度 | θ 1 \theta_1 θ1=-10, θ 2 \theta_2 θ2=10 |
信噪比 | 10dB |
快拍数 | 512 |
基于四阶累积量的SPICE算法相比SPICE计算量较大,这里只给出最终的结果
上图给出的SPICE算法与基于四阶累积量的SPICE(QSPICE)归一化空间谱的对比,从上图可以看出,基于四阶累积量的SPICE算法较SPICE算法有了改善,将上图以对数空间谱画出可以得到如下结果
可以看出四阶累积量的SPICE算法较SPICE算法的改善还是比较大的,虽然说在非目标有较大的起伏,但是是负几百dB,也就是说可以忽略不计,坐标调整后可以得到如下的结果
这样就更能显示出基于四届累计的的SPICE算法的优点了。
这里只给出常规的SPICE算法的代码,关于基于四阶累积量的SPICE算法的代码只需要按照上面的方法对其进行替换即可,这里不再叙述。
关于SPICE算法的代码如下:
clear;
close all;
clc;
derad = pi/180; % 角度->弧度
M = 10; % 阵元个数
theta = [-10 10]; % 待估计角度
K = size(theta,2); % 信源数目
snr = 10; % 信噪比
L = 512; % 快拍数
angleinterv=1; % 角度间隔
angle=-90:angleinterv:90;
dd = 0.5; % 阵元间距,半波长
d=0:dd:(M-1)*dd;
A=exp(-1i*2*pi*d.'*sin(theta*derad)); %导向矢量
iter=100;
%% 构建信号模型
S=randn(K,L); % 信源信号,入射信号
X=A*S; % 构造接收信号
X1=awgn(X,snr,'measured'); % 将白色高斯噪声添加到信号中
%% 计算协方差矩阵
Rxx=X1*X1'/L;
x1=my_SPICE(Rxx,dd,angleinterv,iter);
x11=20*log10(x1/max(x1));
figure(1);
h=plot(angle,x11);
xlim([-40 40]);ylim([-70 0])
set(gca,'fontsize',10);set(h,'Linewidth',2);
xlabel('Angle(degree)','linewidth',2,'fontsize',20,'fontname','Times New Roman');
ylabel('Nomalized power spectrum (dB)','linewidth',2,'fontsize',20,'fontname','Times New Roman')
grid on;hold on;
SPICE算法函数
function x = my_SPICE(R,d,angleinterv,iter)
%R: 协方差矩阵
M0 = size(R,1); %相关矩阵的行数,阵元数目
invR = pinv(R); %相关矩阵的广义逆
theta = -90:angleinterv:90; %角度迭代范围
N = length(theta); %角度范围长度
A = zeros(M0,N+M0);
%产生导向矢量矩阵
A= exp(-1i*2*pi*(0:M0-1)'*d*sin((pi/180)*theta));
A(:,N+1:N+M0) = eye(M0); %后面M0列为单位阵
%初始化: p_ini
p_ini = zeros(1,N+M0);
for k = 1:N+M0
a = zeros(M0,1);
a(:) = A(:,k);
p_ini(k) = (a'*R*a)/(norm(a,2))^4;
end
%初始化迭代
p = zeros(1,N+M0);
W = zeros(1,N+M0);
for k = 1:N+M0
a = zeros(M0,1);
a(:) = A(:,k);
W(k) = (a'*invR*a)/M0;
end
R_ini = A*diag(p_ini.')*A'; %求出协方差矩阵
invR_ini = pinv(R_ini); %求广义逆
rho = 0;
for k = 1:N+M0
a = zeros(M0,1);
a(:) = A(:,k);
rho = rho+sqrt(W(k))*p_ini(k)*norm(a'*invR_ini*R^(1/2),2);
end
for k = 1:N+M0
a = zeros(M0,1);
a(:) = A(:,k);
p(k) = p_ini(k)*(norm(a'*invR_ini*R^(1/2),2))/(sqrt(W(k))*rho);
end
%迭代至收敛
i = 1;
while i <= iter
p_ini(:) = p(:);
R_ini = A*diag(p_ini.')*A';
invR_ini = pinv(R_ini);
rho = 0;
for k = 1:N+M0
a = zeros(M0,1);
a(:) = A(:,k);
rho = rho+sqrt(W(k))*p_ini(k)*norm(a'*invR_ini*R^(1/2),2);
end
for k = 1:N+M0
a = zeros(M0,1);
a(:) = A(:,k);
p(k) = p_ini(k)*(norm(a'*invR_ini*R^(1/2),2))/(sqrt(W(k))*rho);
end
i = i+1;
end
x = zeros(1,N);
x(1,:) = p(1:N);
x = abs(x);
end