SPICE算法与基于四阶累积量的SPICE算法

YHCANDOU

已于 2022-11-26 10:39:55 修改

阅读量3.4k

点赞数 13

分类专栏：阵列信号处理文章标签：算法

于 2020-10-31 16:52:58 首次发布

本文链接：https://blog.csdn.net/yhcwjh/article/details/109275170

版权

阵列信号处理专栏收录该内容

16 篇文章 157 订阅

订阅专栏

SPICE（协方差稀疏迭代估）算法介绍

假设有 $K$ 个远场窄带信号，入射到由 $M\left( M>K \right)$ 个阵元组成的均匀线上，同时信号和噪声互不相关，令 $\Omega$ 表示信源可能的方向， $\left\{ \theta _{k=1}^{Q} \right\}$ 是包含 $\Omega$ 的集合，并且 $Q\gg K$ ，即真实的信源方向包含在 $Q$ 中，那么阵列的输出信号可以表示为
$\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$
其中 $L$ 表示快拍数大小， $\mathbf{Y}\left( t \right)\in {{C}^{M\times 1}}$ 表示第 $t$ 个快拍时的观察数据， $\mathbf{a}\left( {{\theta }_{k}} \right)\in {{C}^{M\times 1}}$ 为入射方向 ${{\theta }_{k}}$ 对应的导向矢量， ${{s}_{k}}\left( t \right)\in C$ 表示信源可能方向 ${{\theta }_{k}}$ 的入射信号， $\mathbf{e}\left( t \right)\in {{C}^{M\times 1}}$ 是零均值的加性高斯噪声。
$\left\{ \mathbf{a}\left( {{\theta }_{k}} \right) \right\}_{k=1}^{Q}$ 表示阵列的过完备基，可以表示为
$\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]$
其中 $\lambda$ 表示入射信号的波长， $d$ 表示阵元之间的间距大小。
那么，阵列的协方差矩阵可以表示为
${{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)$
其中 ${{R}_{s}}\text{=}E\left[ \mathbf{S}\left( t \right){{\mathbf{S}}^{H}}\left( t \right) \right]$ 表示信号的协方差矩阵，
$\mathbf{S}\left( t \right)={{\left[ {{s}_{1}}\left( t \right),{{s}_{2}}\left( t \right),\cdots ,{{s}_{Q}}\left( t \right) \right]}^{T}}$ ， $\mathbf{B}\left( \theta \right)$ 是 $M\times Q$ 的阵列流形矩阵。
对噪声进行求期望可以得到
$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}}}$
其中，
${{\delta }_{t,\bar{t}}}\text{=}\left\{ \begin{aligned} & 1\ \ \ \ \ \ \ t=\bar{t} \\ & 0\ \ \ \ \ \ t\ne \bar{t} \\ \end{aligned} \right.$
协方差矩阵可以写成
$\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}$
为了估计 $\left\{ {{p}_{k}} \right\}_{k=1}^{Q+M}$ ，可以通过以下得到
$f=\ {{\left\| {{R}_{a}}^{-1/2}\left( {{{\hat{R}}}_{a}}-{{R}_{a}} \right){{{\hat{R}}}_{a}}^{-1/2} \right\|}^{2}}$
同时为了求解上式，通过简单的运算，可以将其转化为
$\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}$
其中 ${{R}_{a}}^{-1/2}$ 表示正定矩阵 ${{R}_{a}}^{-1}$ 的平方根， ${{\hat{R}}_{a}}\text{=}\mathbf{Y}{{\mathbf{Y}}^{H}}/M$ ,等价于最小化函数
$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}}$
为了求解上式可以将其转化为带有限制的最小化问题
$\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.$
其中 ${{\omega }_{k}}={{\mathbf{a}}_{k}}^{*}R_{a}^{-1}{{\mathbf{a}}_{k}}/M$ 。
求解约束条件下的方程就可以得到最终的解，可以利用迭代算法来进行求解上述方程，这就是SPICE算法。
因此SPICE算法可以总结为
1）通过时间平均来初始化功率谱估计
$p_{k}^{0}={{\mathbf{a}}_{k}}^{*}{{\hat{R}}_{a}}{{\mathbf{a}}_{k}}/{{\left\| {{\mathbf{a}}_{k}} \right\|}^{4}}\ \ \ \ \ \ k=1,2,\cdots ,Q+M$
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.$
其中 ${{R}_{a}}\left( i \right)$ 上述计算公式重新得到的协方差矩阵。
基于协方差稀疏迭代估计法的波达方向估计，是一种计算效率很高的稀疏迭代算法，该方法以最小化协方差矩阵拟合误差准则得到的迭代算法，在多快拍的情况下具有很好的效果，在少快拍的情况下也具有很好的效果。与其他稀疏算法相比，概算法具有以下特点：
 简单且完美的统计基础
 自然地考虑了数据中的噪声
 不需要设置一些超参数
 具有很好的收敛性

基于四阶累积量的SPICE算法

根据四阶累积量的定义，可以得到阵列输出信号的四阶累积量矩阵，如下
${{C}_{y}}=D\left( \theta \right){{C}_{s}}D{{\left( \theta \right)}^{H}}+{{C}_{e}}$

其中 ${{C}_{y}}$ 为输出信号的四阶累积量矩阵， ${{C}_{s}}$ 为信源的四阶累积量矩阵， ${{C}_{e}}$ 为噪声的四阶累积量矩阵， $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]$ 为阵列流形矩阵。
通过将SPICE算法中的 ${{R}_{a}}$ 和 ${{R}_{a}}$ 替换成 ${{C}_{y}}$ 和 $\left\{ {{\mathbf{d}}_{k}} \right\}$ 就可以得到基于四阶累积量的SPICE算法。但是这种方法的计算量过大，因为四阶累积量矩阵中包含了大量的冗余数据，为了简化计算量，令 $X\left( t \right)={{\left[ {{\mathbf{y}}_{1}}\left( t \right),{{\mathbf{y}}_{M}}\left( t \right) \right]}^{T}}$ ，可以得到
$\mathbf{X}\left( t \right)=\mathbf{\bar{B}}\left( \theta \right)\mathbf{S}\left( t \right)+\mathbf{e}\left( t \right)$
式中， $\mathbf{\bar{B}}\left( \theta \right)\text{=}\left[ \mathbf{\bar{a}}\left( {{\theta }_{1}} \right),\cdots ,\mathbf{\bar{a}}\left( {{\theta }_{Q}} \right) \right]$ ， $\mathbf{\bar{a}}\left( {{\theta }_{k}} \right)={{\left[ 1,{{e}^{-j2\pi \left( M-1 \right)d\sin {{\theta }_{k}}/\lambda }} \right]}^{T}}$ 。
那么四阶累积量矩阵可以变化为
${{\bar{C}}_{y}}=\bar{D}\left( \theta \right){{\bar{C}}_{s}}\bar{D}{{\left( \theta \right)}^{H}}+{{\bar{C}}_{e}}$
其中， $\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]$ 表示导向矢量， ${{\bar{C}}_{y}}$ 是 $2M\times 2M$ 的四阶累积量矩阵。
将SPICE算法中的 ${{R}_{a}}$ 和 $\left\{ {{\mathbf{a}}_{k}} \right\}$ 替换成 ${{\bar{C}}_{y}}$ 和 $\left\{ {{{\mathbf{\bar{d}}}}_{k}} \right\}$ 就可以得到基于四阶累积量的SPICE算法，其步骤如下：
1）通过时间平均来初始化功率谱估计
$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}}$
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.$
基于四阶累积量的协方差稀疏迭代估计算法，不仅提高了原本的超分辨性能，而且在低信噪比的条件下也具有很好的性能，同时通过对均匀线性阵列非冗余数据的提取，降低了计算复杂度。

实验参数设置

参数名称	参数值
阵元数	10
阵元间距	$\lambda/2$
信源数目	3
信源角度	$\theta_1$ =-10, $\theta_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

参考文献

P. Stoica, P. Babu and J. Li, “SPICE: A Sparse Covariance-Based Estimation Method for Array Processing,” in IEEE Transactions on Signal Processing, vol. 59, no. 2, pp. 629-638, Feb. 2011, doi: 10.1109/TSP.2010.2090525.

YHCANDOU

关注

13
点赞
踩
48

收藏

觉得还不错? 一键收藏
10
评论
SPICE算法与基于四阶累积量的SPICE算法

SPICE（协方差稀疏迭代估）算法介绍假设有KKK个远场窄带信号，入射到由 M(M>K)M\left( M>K \right)M(M>K)个阵元组成的均匀线上，同时信号和噪声互不相关，令 Ω\OmegaΩ表示信源可能的方向， {θk=1Q}\left\{ \theta _{k=1}^{Q} \right\}{θk=1Q}是包含Ω\OmegaΩ 的集合，并且Q≫KQ\gg KQ≫K，即真实的信源方向包含在QQQ 中，那么阵列的输出信号可以表示为Y(t)=∑k=1Qa(θk)sk(t)+
复制链接

扫一扫