CN104537249A - 基于稀疏贝叶斯学习的波达方向角估计方法 - Google Patents

Abstract

本发明公开了一种基于稀疏贝叶斯学习的波达方向角估计方法，主要解决现有技术运算量大，处理相干信号源性能差，造成无源定位估计误差大的问题，其实现步骤是：1)利用天线接收机形成均匀线阵，并对空间信号进行采样获得观测数据；2)将观测数据实值化并计算协方差矩阵；3)将空域网格划分，构造实值化的超完备基；4)根据协方差矩阵和超完备基的稀疏表示关系，构建稀疏矩阵方程；5)采用稀疏贝叶斯学习算法求解矩阵方程，获得未知矩阵方差的最稀疏解；6)根据稀疏解与空间角度一一对应的关系，绘制幅度谱图，获得波达方向角度值。本发明提高了无源测向的运算速度及低快拍数时对信号方向角的估计性能，可用于目标侦察和无源定位。

Description

基于稀疏贝叶斯学习的波达方向角估计方法

技术领域

本发明属于信号处理技术领域，特别涉及一种基于均匀线阵的波达方向角估计方法，可用于目标侦察与无源定位。

背景技术

信号的波达方向角DOA估计是阵列信号处理领域的一个重要分支，它是指利用天线阵列对空间信号进行感应接收，再运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。随着科技的不断进步，对信号波达方向估计的精确度和和分辨率也有越来越高的要求。

目前，超分辨DOA估计技术主要有子空间类方法和基于稀疏表示的方法。出现较早，应用较为广泛的是多重信号分类MUSIC等子空间类方法，然而，这些方法依赖于大量采样数据或较高的信噪比才能得到精确的DOA估计。近年来出现的基于稀疏表示的DOA估计方法基本是利用信号的空域稀疏性进行建模，以贪婪算法和凸优化方法为主要手段而展开的。其中贪婪算法在低信噪比情况下，估计性能大幅下降，已不能满足工程需求；而凸优化方法运算速度很慢，且在低信噪比情况下，估计精度不理想，对抗相干信号性能不强。在实际应用中，目标侦察与无源定位均需要在角度估计的基础上进行，以上算法中的缺陷将造成目标侦察和无源定位反应速度慢和估计误差较大的不足。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于稀疏贝叶斯学习算法的波达方向角度估计方法，以在降低运算量的情况下，提高目标侦察和无源定位在低信噪比、低快拍数条件下的估计精度和对相干信号的估计能力，避免因角度估计误差引起的目标侦察失误。

为实现上述目的，本发明的实现步骤包括如下：

1)采用M个天线接收机形成均匀线性阵列，并假设有K个信号入射到该均匀线性阵列，各天线接收机间距均为d，每个天线接收机称为一个阵元，其中，M≥2，K≥1，0<d≤λ/2，λ为入射窄带信号波长；

2)由阵列天线接收机对空间信号进行并行采样，得到输出信号Y(t)；

3)将阵列输出信号Y(t)转换为实值信号矩阵Y_r，并根据实值信号矩阵Y_r，计算阵列协方差矩阵R：

R＝E[Y_r(t)Y_r ^H(t)]，

其中，E[·]表示求数学期望，H表示共轭转置运算；

4)对观测空间进行网格划分，构造实值化的超完备基

4a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角度,定义为波达方向角范围θ＝[θ₁,θ₂,...,θ_q,...,θ_Q]，θ_q为目标信号的来波方向角，q＝1,2,...,Q，Q>>M；

4b)构造一个空域稀疏化后对应的M×Q维的导向矩阵A(θ)：

A(θ)＝[α(θ₁),...,α(θ_q),...,α(θ_Q)]，

其中，α(θ_q)表示方向角θ_q对应的导向矢量：

$\mathrm{\&alpha;}\left({\mathrm{\&theta;}}_q\right)={[1,e^{-\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},...,e^{-(M-1)\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q}]}^T,$

其中，表示相邻两个阵元间的相位差，T表示矩阵转置运算，j为虚数单位；

4c)计算实值化的超完备字典

$\tilde{A}\left(\mathrm{\&theta;}\right)=Q_M^{H}A\left(\text{\&theta;}\right)\mathrm{\&Lambda;}$

其中，Q_M是酉变换矩阵，Λ是一个Q阶对角矩阵,其第q行的对角元素为

5)根据步骤(4)和(5)得到的结果，将波达方向角估计问题转化为求解如下矩阵方程：

$R=\tilde{A}\left(\mathrm{\&theta;}\right)X+{\mathrm{\&sigma;}}^2{I}_M,$

其中X是一个Q×M维的未知矩阵，σ²是加性高斯噪声方差，I_M是M阶单位矩阵；

6)定义一个超参数向量γ＝[γ₁,…,γ_q,…,γ_Q]^T，γ_q为控制矩阵X第q行元素分布的未知先验方差，并采用稀疏贝叶斯学习算法求解该矩阵方程，得到超参数向量γ最稀疏的解γ^*；

7)以波达方向角范围θ＝[θ₁,θ₂,...,θ_q,...,θ_Q]的值为x轴坐标，以γ^*向量的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

本发明与现有技术相比具有以下优点：

1)本发明采用稀疏表示的思想将波达方向角估计问题转化为稀疏重构问题，是新理论技术与传统问题的结合，利用入射信号源的空域稀疏特性进行建模，避免了传统算法的角度搜索或角度匹配过程，同时可以用远低于奈奎斯特采样率所需的采样数据精确估计出波达方向角，极大地降低了信号处理系统的工作负担。

2)本发明利用稀疏贝叶斯学习的统计优化算法求解DOA估计问题中的稀疏矩阵方程，综合考虑了先验分布和观测数据，避免了传统方法只利用观测数据所带来的噪声影响，因而减小了目标侦察和无源定位的估计误差。本方法采用相对门限进行迭代判决，加快了参数收敛速度，在稀疏重构过程中大大提高了运算速度。

3)本发明无需事先估计入射信号的个数，同时可用于处理相干信号源和非相干信号源，在现实环境中具有更广泛的应用价值。

附图说明

图1是本发明的实现流程图；

图2是本发明与现有两种波达方向角估计方法在不同信噪比条件下均方根误差对比图；

图3是本发明与现有两种波达方向角估计方法在不同信噪比条件下检测率对比图；

图4是本发明与现有两种波达方向角估计方法的运算时间对比图；

图5是本发明与现有两种波达方向角估计方法对相干信号源的检测率对比图；

具体实施方式

以下参照附图，对本发明的技术方案和效果作进一步的详细说明。

参照图1，本发明的实现步骤如下：

步骤1：利用天线接收机形成均匀线阵。

每隔间距d放置1个天线接收机，共放置M个，形成一个均匀线阵，每个天线接收机称为一个阵元，假设有K个远场窄带信号入射到该均匀线阵上，且信号在传播过程中加入了均值为0的复高斯白噪声，其中，M≥2，K≥1，0<d≤λ/2，λ为入射窄带信号波长。

步骤2：对空间信号进行并行采样，得到输出信号Y(t)。

由均匀线阵的M个天线接收机以固定的采样频率对空间信号进行并行采样，采样点数为N，得到天线接收机的输出信号Y(t)，其中，N>0且为整数，Y(t)为复数信号。

步骤3：将输出信号Y(t)转换为实值信号矩阵Y_r，并根据实值信号矩阵Y_r，计算阵列协方差矩阵R。

3a)计算输出信号Y(t)的增广数据矩阵Y_aug：

Y_aug＝[Y(t),Π_MY^*(t)Π_N]，

其中，*表示共轭运算，Π_M和Π_N分别表示反对角线元素为1，其余元素均为0的M×M维置换矩阵和N×N维置换矩阵；

3b)根据矩阵理论，将增广数据矩阵Y_aug通过以下变换转换为实值信号矩阵Y_r：

$Y_r=Q_M^H{Y}_{\mathrm{aug}}{Q}_{2N},$

其中，H表示共轭转置运算，Q_2N是一个酉矩阵，Q_M称为酉变换矩阵，其按如下规则计算：

若M为偶数，则 $Q_M=\frac{1}{\sqrt{2}}\left|\begin{array}{cc}{I}_n& j{I}_n\\ {\mathrm{\&Pi;}}_n& -j{\mathrm{\&Pi;}}_n\end{array}\right|,$ 式中n＝M/2，

若M为奇数，则 $Q_M=\frac{1}{\sqrt{2}}\left|\begin{array}{ccc}{I}_n& 0& j{I}_n\\ 0& \sqrt{2}& 0\\ {\mathrm{\&Pi;}}_n& 0& -j{\mathrm{\&Pi;}}_n\end{array}\right|,$ 式中n＝(M-1)/2，

式中，j表示虚数单位，Π_n表示反对角线元素为1，其余元素均为0的n×n维置换矩阵，I_n和I_N均表示单位矩阵；

3c)根据实值信号矩阵Y_r，计算阵列协方差矩阵R：

R＝E[Y_r(t)Y_r ^H(t)]，

其中，E[·]表示求数学期望，H表示共轭转置运算。

步骤4：对观测空间进行网格划分，构造实值化的超完备基Φ(θ)。

根据稀疏信号重构理论，任意信号都可以由一个基矩阵线性表示，此处构造实值化的超完备基Φ(θ)矩阵的目的就是将阵列观测数据的协方差矩阵R，通过矩阵的形式表示出来，便于构建稀疏矩阵方程，其构造步骤如下：

4a)根据入射信号源具有的空域稀疏特性，对观测空域进行空间网格划分处理，即将观测空域[-90°,90°]等间隔划分成Q个角度,θ＝[θ₁,θ₂,...,θ_q,...,θ_Q]，θ表示波达方向角范围，θ_q为第q个角度区间，q＝1,2,...,Q，Q>>M，网格划分间隔的取值根据期望达到的角度估计精度进行设定，网格划分的间隔越小，最终得到的角度估计值精度越高；

4b)构造一个信号稀疏化后对应的M×Q维的导向矩阵A(θ)：

A(θ)＝[α(θ₁),...,α(θ_q),...,α(θ_Q)]，

其中，α(θ_q)表示方向角θ_q对应的导向矢量：

$\mathrm{\&alpha;}\left({\mathrm{\&theta;}}_q\right)={[1,e^{-\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},...,e^{-(M-1)\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q}]}^T,$

其中，表示相邻两个阵元间的相位差，T表示矩阵转置运算；

4c)将导向矩阵A(θ)作线性变换，得到实值化的超完备基Φ(θ)：

$\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)=Q_M^{H}A\left(\text{\&theta;}\right)\mathrm{\&Lambda;},$

其中称为基向量，Q_M是酉变换矩阵，Λ为一个Q阶对角矩阵,其第q个对角元素为，q＝1,2,…,Q。

步骤5：根据步骤(3)和(4)得到的结果，将波达方向角估计问题转化为求解如下稀疏矩阵方程：

R＝Φ(θ)X+σ²I_M，

其中X是一个Q×M维的未知矩阵，σ²是加性高斯噪声方差，I_M是M阶单位矩阵；

定义一个超参数向量γ＝[γ₁,...,γ_q,...,γ_Q]^T，γ_q为控制矩阵X第q行元素分布的未知先验方差，根据实际环境中噪声服从方差为σ²的高斯分布的特性，可知X的每一行均服从一个均值为0、方差为γ_q的高斯先验分布，通过贝叶斯准则，进一步将波达方向角估计问题的求解转化为对超参数向量γ的求解。

步骤6：采用稀疏贝叶斯学习算法求解上述稀疏矩阵方程，得到超参数向量γ最稀疏的解γ^*。

6a)将超参数向量γ初始化为所有元素均为1的向量，设定加性高斯噪声方差σ²的初始值为阵列协方差矩阵R的最小特征值；

6b)根据超完备基Φ(θ)和阵列协方差矩阵R,依据下式计算迭代过程中未知矩阵X的均值μ和方差V：

$\mathrm{\&mu;}=\mathrm{\&Gamma;}{\mathrm{\&Phi;}}^T\left(\mathrm{\&theta;}\right)V_R^{-1}R,$

$V=\mathrm{\&Gamma;}-\mathrm{\&Gamma;}{\mathrm{\&Phi;}}^T\left(\mathrm{\&theta;}\right)V_R^{-1}\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)\mathrm{\&Gamma;}$

其中，Γ＝diag(γ)，V_R＝σ²I_M+Φ(θ)ΓΦ^T(θ)，(·)^-1为矩阵求逆运算，diag表示构造对角矩阵；

6c)采用均值最大EM准则分别更新超参数向量γ的第i个元素γ_i和噪声方差σ²，得到更新后的元素γ_i′和噪声方差(σ²)′：

${{\mathrm{\&gamma;}}_i}^{\&prime;}=\frac{(1/M){\left|\right|{\mathrm{\&mu;}}_i\left|\right|}_2^2}{1-{\mathrm{\&gamma;}}_i^{-1}{V}_{\mathrm{ii}}},{\left({\mathrm{\&sigma;}}^2\right)}^{\&prime;}=\frac{(1/M){\left|\right|R-\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)X\left|\right|}_F^2}{M-Q+{\mathrm{\&Sigma;}}_{i=1}^Q(V_{\mathrm{ii}}/{\mathrm{\&gamma;}}_i)},$

其中，V_ii为方差V的第i行和第i列对应的元素，μ_i为均值μ的第i行元素组成的向量，||·||₂,||·||_F分别表示求2范数和F范数；

6d)计算更新后的元素γ_i′与超参数向量γ中最大元素值的相对残差ξ：

ξ＝10lg(γ_i′/max(γ))，

将相对残差ξ与一个确定门限值进行比较，该门限值的大小根据应用环境进行设置,通常其值设置为-30dB即可较好地平衡算法的收敛速度和精度，若ξ<-30dB，则将元素γ_i′及其在超完备基Φ(θ)中对应的基向量分别置零，若ξ>-30dB，则保留元素γ_i′及其在超完备基Φ(θ)中对应的基向量进入下一次迭代计算；

6e)迭代计算步骤6b)到步骤6d)，直至满足max(max|μ′-μ|)<ε时停止，此时超参数向量γ收敛到最稀疏的解γ^*，其中μ′为上次迭代过程中的均值，ε为迭代停止门限，其值设定为10^-8。

步骤7：根据稀疏解γ^*向量绘制幅度谱图，得到波达方向角的估计值。

步骤6得到的最稀疏的解γ^*向量是一个K稀疏向量，即其中只有K个值为非零值，其余值均为零，这K个非零值对应的空间方向角就是入射信号的波达方向角，因此，以波达方向角范围θ＝[θ₁,θ₂,...,θ_q,...,θ_Q]的值为x轴坐标，以γ^*向量的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

本发明的效果可通过以下仿真说明：

1.仿真条件：

采用10个天线接收机形成均匀线阵，各天线接收机的间距d为入射信号波长的一半，采样点数为100，观测空域角度范围为[-90°,90°]，空间网格划分间隔为1°。

均方根误差的计算式为：

$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{KJ}}\underset{j=1}{\overset{J}{\mathrm{\&Sigma;}}}\underset{k=1}{\overset{K}{\mathrm{\&Sigma;}}}{({\widehat{\mathrm{\&theta;}}}_{\mathrm{kj}}-{\mathrm{\&theta;}}_k)}^2}$

其中，J表示实验次数，J＝100，表示第j次实验的DOA估计值，θ_k表示信号的DOA真实值。此处的检测率表示100次独立实验中检测成功的次数，其中一次实验检测成功定义为该次所有DOA估计值均在真实值±1°范围内。

2.仿真内容与结果：

仿真1：假设有2个非相干窄带信号分别以角度6°和18°入射到均匀线阵，信噪比由-8dB增加到4dB，利用本发明和现有的L1-SVD、L1-SRACV算法分别进行100次独立的波达方向角估计实验，分别计算不同信噪比下三种方法的均方根误差RMSE和检测率，实验结果如图2和图3所示，其中：

图2中横坐标表示信噪比值，纵坐标表示均方根误差；

图3中横坐标表示信噪比值，纵坐标表示检测率。

从图2可以看出，在低信噪比的情况下，本发明的均方根误差均小于其它两个算法；

从图3可以看出，本发明的检测率明显高于其它两个算法。

仿真2：假设有2个非相干窄带信号分别以角度16°和26°入射到均匀线阵，信噪比设定为0dB，阵列的阵元数由4增加到12，利用本发明和现有的L1-SVD、L1-SRACV算法分别进行100次独立的波达方向角度估计实验，分别统计每次实验条件下的运算时间，结果如图4所示。图4中横坐标表示阵元数，纵坐标表示运算时间。

从图4可以看出，本发明与现有L1-SVD、L1-SRACV算法相比，大幅度降低了DOA估计的运算时间，即使阵元数增加，运算时间仍基本保持不变。

仿真3：假设2个相干窄带信号分别以角度16°和26°入射到均匀线阵，信噪比固定为0dB。采样快拍数由20增加至200，利用本发明和现有的L1-SVD、L1-SRACV算法分别进行100次独立的波达方向角度估计实验，统计不同快拍数下的检测率，结果如图5所示，图5中横坐标表示快拍数，纵坐标表示检测率。

从图5可以看出，在快拍数低于100时，本发明的检测率明显高于其它两种方法，说明本发明处理相干信号时有明显的优势。

综上，本发明在降低了DOA估计运算量的同时，提高了低信噪比情况下角度估计的检测率，降低了估计误差，并且兼具非相干和相干信号源的估计能力，保证了目标侦察和无源定位的快速反应和准确有效，避免了相干信号源背景下角度估计失误的问题。

Claims (3)

1.一种基于稀疏贝叶斯学习的波达方向角估计方法，包括以下步骤：

1)采用M个天线接收机形成均匀线性阵列，并假设有K个信号入射到该均匀线性阵列，各天线接收机间距均为d，每个天线接收机称为一个阵元，其中，M≥2，K≥1，0<d≤λ/2，λ为入射窄带信号波长；

2)由阵列天线接收机对空间信号进行并行采样，得到输出信号Y(t)；

3)将阵列输出信号Y(t)转换为实值信号矩阵Y_r，并根据实值信号矩阵Y_r，计算阵列协方差矩阵R：

R＝E[Y_r(t)Y_r ^H(t)]，

其中，E[·]表示求数学期望，H表示共轭转置运算；

4)对观测空间进行网格划分，构造实值化的超完备基Φ(θ)：

4a)根据信号源的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个角度,定义为波达方向角范围θ＝[θ₁,θ₂,...,θ_q,...,θ_Q]，θ_q为目标信号的来波方向角，q＝1,2,...,Q，Q>>M；

4b)构造一个空域稀疏化后对应的M×Q维的导向矩阵A(θ)：

A(θ)＝[α(θ₁),...,α(θ_q),...,α(θ_Q)]，

其中，α(θ_q)表示方向角θ_q对应的导向矢量：

$\mathrm{\&alpha;}\left({\mathrm{\&theta;}}_q\right)={[1,e^{-\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},...,e^{-(M-1)\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q}]}^T,$

其中，表示相邻两个阵元间的相位差，T表示矩阵转置运算，j为虚数单位；

4c)计算实值化的超完备基Φ(θ)：

$\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)=Q_M^{H}A\left(\mathrm{\&theta;}\right)\mathrm{\&Lambda;}$

其中，q＝1,2,…,Q，称为基向量，Q_M是酉变换矩阵，Λ是一个Q阶对角矩阵,其第q行的对角元素为

5)根据步骤(4)和(5)得到的结果，将波达方向角估计问题转化为求解如下矩阵方程：

R＝Φ(θ)X+σ²I_M,

其中X是一个Q×M维的未知矩阵，σ²是加性高斯噪声方差，I_M是M阶单位矩阵；

6)定义一个超参数向量γ＝[γ₁,...,γ_q,...,γ_Q]^T，γ_q为控制矩阵X第q行元素分布的未知先验方差，并采用稀疏贝叶斯学习算法求解该矩阵方程，得到超参数向量γ最稀疏的解γ^*；

7)以波达方向角范围θ＝[θ₁,θ₂,...,θ_q,...,θ_Q]的值为x轴坐标，以γ^*向量的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

2.根据权利要求1所述的一种基于稀疏贝叶斯学习的波达方向角估计方法，其中步骤3)所述的将阵列输出信号Y(t)转换为实值信号矩阵Y_r，按如下步骤进行：

3a)计算输出信号Y(t)的增广数据矩阵Y_aug：

Y_aug＝[Y(t),Π_MY^*(t)Π_N]，

其中，*表示共轭运算，N表示采样点数，Π_M和Π_N分别表示反对角线元素为1，其余元素均为0的M×M维置换矩阵和N×N维置换矩阵；

3b)根据矩阵理论，将增广数据矩阵Y_aug通过以下变换转换为实值信号矩阵Y_r：

$Y_r=Q_M^H{Y}_{\mathrm{aug}}{Q}_{2N},$

其中，H表示共轭转置运算，Q_2N是一个酉矩阵， $Q_{2N}=\frac{1}{\sqrt{2}}\left|\begin{array}{cc}{I}_N& {\mathrm{jI}}_N\\ {\mathrm{\&Pi;}}_N& -{\mathrm{j\&Pi;}}_N\end{array}\right|,$

Q_M称为酉变换矩阵，其按如下规则计算：

若M为偶数，则 $Q_M=\frac{1}{\sqrt{2}}\left|\begin{array}{cc}{I}_n& {\mathrm{jI}}_n\\ {\mathrm{\&Pi;}}_n& -{\mathrm{j\&Pi;}}_n\end{array}\right|,$ 式中n＝M/2，

若M为奇数，则 $Q_M=\frac{1}{\sqrt{2}}\left|\begin{array}{ccc}{I}_n& 0& {\mathrm{jI}}_n\\ 0& \sqrt{2}& 0\\ {\mathrm{\&Pi;}}_n& 0& -{\mathrm{j\&Pi;}}_n\end{array}\right|,$ 式中n＝(M-1)/2，

式中，Π_n表示反对角线元素为1，其余元素均为0的n×n维置换矩阵，I_n和I_N均表示单位矩阵。

3.根据权利要求1所述的一种基于稀疏贝叶斯学习的波达方向角估计方法，其中步骤6)所述的采用稀疏贝叶斯学习算法求解该矩阵方程，得到超参数向量γ最稀疏的解γ^*，按如下步骤进行：

6a)将超参数向量γ初始化为所有元素均为1的向量，设定加性高斯噪声方差σ²的初始值为阵列协方差矩阵R的最小特征值；

6b)根据超完备基Φ(θ)和阵列协方差矩阵R,计算迭代过程中未知矩阵X的均值μ和方差V：

$\mathrm{\&mu;}={\mathrm{\&Gamma;\&Phi;}}^T\left(\mathrm{\&theta;}\right)V_R^{-1}R$

$V=\mathrm{\&Gamma;}-\mathrm{\&Gamma;}{\mathrm{\&Phi;}}^T\left(\mathrm{\&theta;}\right)V_R^{-1}\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)\mathrm{\&Gamma;}$

其中：Γ＝diag(γ)，V_R＝σ²I_M+Φ(θ)ΓΦ^T(θ)，(·)^-1为矩阵求逆运算，diag表示构造对角矩阵；

6c)采用均值最大EM准则分别更新超参数向量γ的第i个元素γ_i和噪声方差σ²，得到更新后的元素γ_i′和噪声方差(σ²)′：

${{\mathrm{\&gamma;}}_i}^{\&prime;}=\frac{(1/M){\left|\right|{\mathrm{\&mu;}}_i\left|\right|}_2^2}{1-{{\mathrm{\&gamma;}}_i}^{-1}{V}_{\mathrm{ii}}},$ ${\left({\mathrm{\&sigma;}}^2\right)}^{\&prime;}=\frac{(1/M){\left|\right|R-\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)X\left|\right|}_F^2}{M-Q+{\mathrm{\&Sigma;}}_{i=1}^Q(V_{\mathrm{ii}}/{\mathrm{\&gamma;}}_i)},$

其中，V_ii为方差V的第i行和第i列对应的元素，μ_i为均值μ的第i行元素组成的向量，i＝1,…,Q,||·||₂,||·||_F分别表示求2范数和F范数；

6d)计算更新后的元素γ_i′与超参数向量γ中最大元素值的相对残差ξ：

ξ＝10lg(γ_i′/max(γ))，

若ξ<-30dB，则将元素γ_i′及其在超完备基Φ(θ)中对应的基向量分别置零，若ξ>-30dB，则保留元素γ_i′及其在超完备基Φ(θ)中对应的基向量进入下一次迭代计算；

6e)迭代计算步骤6b)到步骤6d)，直至满足max(max|μ′-μ|)<ε时停止，得到最稀疏的解γ^*，其中μ′为上次迭代过程中的均值，ε为迭代停止门限，其值为10^-8。