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

Abstract

本发明公开了一种基于快速稀疏贝叶斯学习的波达方向角估计方法，主要解决现有技术运算量大，定位估计误差大的问题，其实现步骤是：1)采用天线接收机形成均匀线阵；2)对空间信号进行采样并计算阵列协方差矩阵R；3)将R矢量化后得到稀疏模型向量y；4)将空域网格划分，根据稀疏模型向量y的结构构造超完备基Φ(θ)；5)根据稀疏模型向量和超完备基的稀疏表示关系，建立稀疏方程；6)定义超参数向量α，采用快速稀疏贝叶斯学习算法求解该稀疏方程；7)根据α的最优估计值绘制幅度谱图，获得波达方向角度值。本发明提高了目标侦察和无源定位在低信噪比和低快拍数条件下的估计精度，降低了运算复杂度，可用于目标侦察和无源定位。

Description

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

技术领域

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

背景技术

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

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

发明内容

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

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

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

2)由阵列天线接收机对空间信号进行采样，得到输出信号Y(t)，并根据该输出信号，计算阵列协方差矩阵R：

R＝E[Y(t)Y^H(t)]

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

3)根据阵列协方差矩阵R构造稀疏模型向量y：

y＝vec(R)，其中，vec(·)表示向量化运算；

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

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

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

B(θ)＝[b(θ₁),...,b(θ_q),...,b(θ_Q)]，

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

$b\left({\mathrm{\&theta;}}_q\right)={[e^{(M-1)\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},...,e^{\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},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)计算选择矩阵G：

$G=[\mathrm{vec}\left(J_{M-1}\right),...,\mathrm{vec}\left(J_1\right),\mathrm{vec}\left(J_0\right),\mathrm{vec}\left(J_1^T\right),...,\mathrm{vec}\left(J_{M-1}^T\right)],$

其中，J₀,J₁,…,J_M-1按下式计算：

$J_l=\left[\begin{array}{cc}{0}_{M-l,l}& I_{M-l}\\ 0_{l,l}& 0_{l,M-l}\end{array}\right],l=\mathrm{0,1},...,M-1;$

其中，I_M-l表示M-l阶的单位矩阵，0_m-l,l,0_l,l,0_l,M-l分别表示m-l×l,l×l,l×m-l维的零矩阵；

4d)根据选择矩阵G和导向矩阵B(θ)，得到超完备基Ф(θ)：

Ф(θ)＝G B(θ)，

其中，称为基向量；

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

y＝Φ(θ)w+σ²vec(I_M),

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

6)定义一个超参数向量α＝[α₁,...,α_q,...,α_Q]^T，α_q为控制w分布的未知先验方差，称为超参数，并采用快速稀疏贝叶斯学习算法求解该稀疏优化方程，得到超参数向量α的收敛解；

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

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

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

2)本发明采用矢量化建模方法，将稀疏模型从多测量矢量模型变为单测量矢量模型，大大降低了稀疏模型的维度，从而在运算量及重构速度上得到很好的改善。

3)本发明利用快速稀疏贝叶斯学习的统计优化算法求解DOA估计问题中的稀疏矩阵方程，该方案综合考虑了先验分布和观测数据，对噪声的初始化不敏感，因而避免了传统方法因噪声方差估计错误或误差较大时带来的噪声影响，从而提高了低信噪比条件下的波达方向估计效果，减小了目标侦察和无源定位的估计误差。

附图说明

图1是本分明的实现场景图；

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

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

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

图5是本发明与现有两种波达方向角估计方法在不同快拍数条件下的均方根误差对比图。

具体实施方式

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

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

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

步骤1：计算均匀线性天线阵列的协方差矩阵R。

用均匀线性阵列的M个天线接收机以固定的采样频率对空间信号进行并行采样，采样点数为N，得到天线接收机的输出信号Y(t)，并根据输出信号Y(t)，计算阵列协方差矩阵R：

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

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

步骤2：对阵列协方差矩阵R进行矢量化运算，得到稀疏模型向量y：

y＝vec(R)，其中，vec(·)表示向量化运算。

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

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

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

4b)构造一个空域稀疏化后对应的(2M-1)×Q维的导向矩阵B(θ):

B(θ)＝[b(θ₁),...,b(θ_q),...,b(θ_Q)]，

其中，b(θ_q)表示B(θ)的第q列，b(θ_q)为角度θ_q对应的2M-1维导向向量：

$b\left({\mathrm{\&theta;}}_q\right)={[e^{(M-1)\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},...,e^{\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},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)计算选择矩阵G：

$G=[\mathrm{vec}\left(J_{M-1}\right),...,\mathrm{vec}\left(J_1\right),\mathrm{vec}\left(J_0\right),\mathrm{vec}\left(J_1^T\right),...,\mathrm{vec}\left(J_{M-1}^T\right)],$

其中，J₀,J₁,…,J_M-1按下式计算：

$J_l=\left[\begin{array}{cc}{0}_{M-l,l}& I_{M-l}\\ 0_{l,l}& 0_{l,M-l}\end{array}\right],l=\mathrm{0,1},...,M-1;$

其中，I_M-l表示M-l阶的单位矩阵，0_m-l,l,0_l,l,0_l,M-l分别表示m-l×l,l×l,l×m-l维的零矩阵；

4d)根据选择矩阵G和导向矩阵B(θ)，得到超完备基Ф(θ)：

Ф(θ)＝G B(θ)，

其中，称为基向量。

步骤5：根据波达方向的空域稀疏性，将波达方向角估计问题转化为超完备基Φ(θ)中各个基向量对应的系数的计算。

定义w是一个Q×1维的未知稀疏向量，则w中的大系数对应基向量的角度即为所求的波达方向角，将稀疏模型向量y用超完备基Φ(θ)表示，则波达方向角估计可转化为求解如下稀疏方程：

y＝Φ(θ)w+σ²vec(I_M),

其中σ²为加性高斯噪声方差，I_M是M阶单位矩阵。

步骤6：采用快速稀疏贝叶斯学习算法求解上述稀疏矩阵方程，得到超参数向量α的收敛解。

由于步骤5所述稀疏方程的求解通常需要较为准确的噪声估计，在低信噪比条件下性能急剧恶化，因此本发明采用快速稀疏贝叶斯学习方法，定义一个超参数向量α＝[α₁,...,α_q,...,α_Q]^T，α_q为控制w分布的未知先验方差，称为超参数，根据实际环境中噪声服从方差为σ²的高斯分布的特性，可知w是服从一个均值为0，方差为α的高斯先验分布，通过贝叶斯准则，将对未知稀疏向量w的求解转化为对超参数向量α的求解。一旦确定了α的收敛解中非零元素的位置，即可得到入射信号的DOA估计值，其求解步骤如下：

6a)设定噪声方差σ²的初始值为0.1var(y)，此处的初始值取值可以不固定，因为该方法对噪声方差初始值不敏感，其初始值对收敛的结果影响不大；初始化超参数向量α为第一个元素其余元素均为无穷大的向量，其中，var(·)表示求方差运算，为超完备基Φ中的第一个基向量，∥·∥表示求1范数；

6b)计算未知向量x的方差V和均值μ：

$V={(\frac{1}{{\mathrm{\&sigma;}}^2}{\mathrm{\&Phi;}}^T\left(\mathrm{\&theta;}\right)\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)+\mathrm{diag}\left(\mathrm{\&alpha;}\right))}^{-1}$

$\mathrm{\&mu;}=\frac{1}{{\mathrm{\&sigma;}}^2}V{\mathrm{\&Phi;}}^H\left(\mathrm{\&theta;}\right)y;$

6c)计算所有基向量对应的质量因子q_i和稀疏因子s_i：

其中，C＝σ^-2I_M-σ^-2I_MΦ(θ)VΦ^T(θ)σ^-2I_M，i＝1,2,...,Q，表示超完备基Φ(θ)矩阵的第i行向量，T表示转置运算；

6d)计算偏移角度如果β_i>0并且α_i≤∞，则更新超参数如果β_i≤0并且α_i<∞，则更新超参数α_i＝∞；

6e)更新噪声方差σ²，得到更新后的噪声方差(σ²)′：

${\left({\mathrm{\&sigma;}}^2\right)}^{\&prime;}=\frac{{\left|\right|y-\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)\mathrm{\&mu;}\left|\right|}^2}{M-Q+V_i{\mathrm{\&alpha;}}_i{V}_{\mathrm{ii}}},$

其中，M为阵元数目，Q为空间网格划分数目，V_ii为方差V的第i行和第i列对应的元素，V_i为的第i行元素组成的向量，i＝1,2,...,Q；

6f)从超完备基Ф中任意选取一个基向量作为候选基向量，并根据更新后的超参数向量α再次计算均值和方差，得到更新后的均值μ′和方差V′；

6g)判断是否满足max(|μ-μ′|)<ε，满足则算法结束，得到超参数向量α的收敛解，否则转到步骤6c)继续迭代计算，其中ε为迭代停止门限，其值可取为10^-8。

步骤7：根据α的收敛解绘制幅度谱图，得到波达方向角的估计值。

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

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

1.仿真条件与方法：

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

2.仿真内容与结果：

仿真1：假设有2个非相干窄带信号分别以角度6°和18°入射到均匀线阵，信噪比为-4dB,利用本发明和现有的L1_SRACV，L1_SVD算法分别进行100次独立的波达方向角估计试验，分别计算不同阵元数条件下三种方法的运算时间，结果如图3所示，图3中的横坐标为阵元数目，单位为个，纵坐标为运算时间，单位为秒。

从图3可以看出，在相同条件下，本发明需要的运算时间要少于现有的两种方法，且随着阵元数目的增加，本发明并不需要额外的运算增加。

仿真2：假设有2个非相干窄带信号分别以角度6°和18°入射到均匀线阵，信噪比SNR由-8dB增加到4dB，利用本发明和现有的PL1_SRACV，MTCS算法分别进行100次独立的波达方向角估计试验，分别计算不同信噪比下三种方法的均方根误差RMSE和检测率，此处的均方根误差按下式计算:

$\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真实值。试验结果如图4所示，图4中横坐标均表示信噪比值，单位dB，图4中纵坐标表示均方根误差，单位为度。

从图4可以看出，在低信噪比的情况下，本发明的均方根误差明显低于其它两个算法，本发明在低信噪比条件下表现出良好的性能。

仿真3：假设有2个非相干窄带信号分别以角度16°和26°入射到均匀线阵，信噪比设定为0dB，采样快拍数由10增加至90，利用本发明和现有的bcs fast rvm，MTCS方法分别进行100次独立的波达方向角度估计实验，分别统计每次实验条件下的运算时间，结果如图5所示。图5中横坐标表示快拍数，单位为个，纵坐标表示均方根误差，单位为度。

从图5可以看出，在低快拍的情况下，本发明的均方根误差要小于其他的两种方法。

综上，本发明在降低了DOA估计运算量的同时，降低了低信噪比，低快拍情况下角度估计的估计误差，保证了目标侦察和无源定位的快速反应和准确有效。

Claims (2)

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

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

2)由阵列天线接收机对空间信号进行采样，得到输出信号Y(t)，并根据该输出信号，计算阵列协方差矩阵R：

R＝E[Y(t)Y^H(t)]

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

3)根据阵列协方差矩阵R构造稀疏模型向量y：

y＝vec(R)，其中，vec(·)表示向量化运算；

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

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

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

B(θ)＝[b(θ₁),...,b(θ_q),...,b(θ_Q)]，

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

$b\left({\mathrm{\&theta;}}_q\right)=[{e^{(M-1)\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},...,e^{\frac{j2\mathrm{\&pi;d}}{\mathrm{\&lambda;}}\sin {\mathrm{\&theta;}}_q},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)计算选择矩阵G：

$G=[\mathrm{vec}\left(J_{M-1}\right),...,\mathrm{vec}\left(J_1\right),\mathrm{vec}\left(J_0\right),\mathrm{vec}\left(J_1^T\right),...,\mathrm{vec}\left(J_{M-1}^T\right)$

其中，J₀,J₁,…,J_M-1按下式计算：

$J_l=\left[\begin{array}{cc}{0}_{M-l,l}& I_{M-1}\\ 0_{l,l}& 0_{l,M-l}\end{array}\right],l=\mathrm{0,1},...,M-1;$

其中，I_M-l表示M-l阶的单位矩阵，0_m-l,l,0_l,l,0_l,M-l分别表示m-l×l,l×l,l×m-l维的零矩阵；

4d)根据选择矩阵G和导向矩阵B(θ)，得到超完备基Ф(θ)：

Ф(θ)＝G B(θ)，

其中，q＝1,2,…,Q，称为基向量；

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

y＝Φ(θ)w+σ²vec(I_M)

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

6)定义一个超参数向量α＝[α₁,...,α_q,...,α_Q]^T，α_q为控制w分布的未知先验方差，称为超参数，并采用快速稀疏贝叶斯学习算法求解该稀疏优化方程，得到超参数向量α的收敛解；

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

2.根据权利要求1所述的一种基于快速稀疏贝叶斯学习的波达方向角估计方法，其中所述步骤6)中采用快速稀疏贝叶斯学习算法求解稀疏优化方程，按如下步骤进行：

6a)设定噪声方差σ²的初始值为0.1var(y)，定义超参数向量α，其第i个值为超参数α_i，i＝1,2,…,Q，初始化其余超参数均为无穷大，其中，var(··)表示求方差运算，为超完备基Φ中的第一个基向量，||·||表示求矩阵2范数；

6b)计算未知向量x的方差V和均值μ：

$V={(\frac{1}{{\mathrm{\&sigma;}}^2}{\mathrm{\&Phi;}}^T\left(\mathrm{\&theta;}\right)\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)+\mathrm{diag}\left(\mathrm{\&alpha;}\right))}^{-1},\mathrm{\&mu;}=\frac{1}{{\mathrm{\&sigma;}}^2}V{\mathrm{\&Phi;}}^H\left(\mathrm{\&theta;}\right)y,$

其中diag(i)表示对角化操作；

6c)计算所有基向量对应的质量因子q_i和稀疏因子s_i：

其中，C＝σ^-2I_M-σ^-2I_MΦ(θ)VΦ(θ)^Tσ^-2I_M，i＝1,2,...,Q，表示超完备基Φ(θ)矩阵的第i列向量，I_M为M阶单位矩阵,σ²为噪声方差；

6d)计算偏移角度如果β_i>0并且α_i≤∞，则更新超参数如果β_i≤0并且α_i<∞，则更新超参数α_i＝∞；

6e)更新噪声方差σ²，得到更新后的噪声方差(σ²)′：

${\left({\mathrm{\&sigma;}}^2\right)}^{\&prime;}=\frac{{\left|\right|y-\mathrm{\&Phi;}\left(\mathrm{\&theta;}\right)\mathrm{\&mu;}\left|\right|}^2}{M-Q+V_i{\mathrm{\&alpha;}}_i{V}_{\mathrm{ii}}},$

其中，Q为空间网格划分数目，V_ii为方差V的第i行和第i列对应的元素，V_i为V的第i行元素组成的向量，i＝1,2,...,Q；

6f)从超完备基Φ(θ)中任意选取一个基向量作为候选基向量，并根据更新后的超参数向量α返回到6b)再次计算均值和方差，得到更新后的均值μ′和向量x的方差V′；

6g)判断是否满足max(|μ-μ′|)<ε，若满足，则结束迭代，得到超参数向量α的收敛解，否则，返回步骤6c)继续迭代计算，其中ε为迭代停止门限，其取值为10^-8。