CN112924925A - 基于稀疏贝叶斯学习的机载三维异构阵doa估计方法 - Google Patents

Abstract

本发明公开了一种基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，包括步骤：建立机载三维异构阵的阵列模型；利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换；通过过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题；利用稀疏贝叶斯学习算法求解。本发明将欧拉旋转变换引入到机载三维异构阵的方向图求解中，得到更为通用的机载三维异构阵阵列流形模型，从而构建出适用于机载三维异构阵的DOA估计模型。

Description

基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法

技术领域

本发明涉及雷达信号处理技术领域，具体涉及基于稀疏贝叶斯学习的机载三维异构阵DOA(波达方向)估计方法。

背景技术

机载三维异构阵是由多个不同的曲面阵构成的三维复杂阵列，它具有与载机外形相一致的空气动力外形，可以有效减小载机雷达反射截面积、降低载机负荷及增大有效发射孔径，能够满足复杂战场环境下雷达对目标探测的需求。相比于传统的二维平面阵，三维异构阵是一种更加广义的构型，可以看作是多个共形阵的组合阵。因此，三维异构阵的DOA估计方法可以参考共形阵的研究方法，但是目前很多研究方法如MUSIC(多信号分类)和ESPRIT(旋转不变性子空间)方法虽然具有一定效果，但是都忽略了载体曲率带来的对天线极化方向图的影响，并且受到算法需要大量快拍数据支撑的限制。

近年来，随着稀疏恢复类算法的兴起，如OMP(正交匹配追踪)算法、IAA(迭代自适应)算法、SBL(稀疏贝叶斯学习)算法等被应用到阵列的DOA估计中，通过构建压缩感知模型，使得少量快拍的DOA估计成为可能。但是，由于机载三维异构阵的复杂性，线阵平面阵的稀疏DOA估计方法无法适用，仍需要对相应模型进行改进。

发明内容

针对现有技术存在的问题，本发明的目的在于提供一种基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，将欧拉旋转变换引入到机载三维异构阵的方向图求解中，得到更为通用的机载三维异构阵阵列流形模型，从而构建出适用于机载三维异构阵的DOA估计模型。

为了达到上述目的，本发明采用以下技术方案予以实现。

基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，包括以下步骤：

步骤1，建立机载三维异构阵的阵列模型；

步骤2，利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换，得到机载三维异构阵的全局天线方向图；

步骤3，获取机载三维异构阵的接收信号，并建立对应的过完备字典矩阵；通过所述过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题，并建立对应的压缩感知模型；

步骤4，利用稀疏贝叶斯学习算法求解所述压缩感知模型，得到机载三维异构阵的DOA估计结果。

与现有技术相比，本发明的有益效果为：

(1)本发明方法考虑了机载三维异构阵的载体曲率对方向图综合的影响，DOA估计模型更加符合实际情况。

(2)本发明方法可以实现少量快拍情况下机载三维异构阵的DOA估计；且本发明方法无需目标数量的先验信息。

(3)本发明方法首次将压缩感知模型以及稀疏贝叶斯方法引进到机载三维异构阵的DOA估计上面，对后续学者的相关研究具有一定的指引作用。

附图说明

下面结合附图和具体实施例对本发明做进一步详细说明。

图1是本发明提供一种基于稀疏贝叶斯学习的机载三维异构阵的DOA估计方法的流程示意图；

图2是本发明实施例提供的机载三维异构阵的几何特性示意图

图3是本发明实施例提供的机载三维异构阵的阵元布置示意图；

图4是本发明实施例提供的机载三维异构阵的阵元坐标图；

图5是本发明实施例提供的异构阵三轴欧拉旋转变换过程示意图；

图6本发明实施例提供的圆柱阵和圆锥阵的欧拉旋转变换过程示意图；其中，6(a)对应圆柱阵，6(b)对应圆锥阵；

图7本发明实施例的单次快拍的DOA估计结果图；其中，7(a)对应目标1的DOA估计结果，7(b)对应目标2的DOA估计结果，7(c)对应目标3的DOA估计结果。

具体实施方式

下面将结合实施例对本发明的实施方案进行详细描述，但是本领域的技术人员将会理解，下列实施例仅用于说明本发明，而不应视为限制本发明的范围。

参考图1，本发明提供的一种基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，包括以下步骤：

步骤1，建立机载三维异构阵的阵列模型；

参考图2，机载三维异构阵由四个部分的共形阵组成，分别为机头的圆锥阵、机身圆柱阵、机身圆台阵和机翼阵。下面以其中的圆柱共形阵和圆锥共形阵组成的三维异构阵为例，进行本发明方法的介绍。载机沿y轴正方向飞行，速度矢量为v,飞行高度为H，异构阵的几何关系图如图3所示，阵元坐标图如图4所示。

参考图3，机载圆柱阵由M₁个圆环阵组成，每个圆环阵由关于圆心对称的N₁个天线阵元组成。M₁个圆环阵彼此平行，阵元间距和圆环阵间距都为d,且共同垂直于y轴。每个圆环阵半径为r，第1个圆环阵位于xoz平面。以圆环阵上顶部阵元为第1个阵元，阵元序号为逆时针排序，则圆柱阵第m₁个阵列圆环上第n₁个阵元的位置矢量

为：

圆柱阵对应的阵元位置坐标为：

圆柱阵整体的位置矢量：

与圆柱阵相比，圆锥阵由不同半径、不同阵元个数的M₂个圆环阵组成。将第m₂(m₂＝1,2,…,M₂)个圆环阵的半径记为

阵元个数记为

底层半径与阵元个数与圆柱阵一致，圆柱阵半径和阵元数的表达式为：

第1个圆环阵离原点距离为D₂，同样以圆环阵上顶部阵元为第1个阵元，阵元序号为逆时针排序，则圆锥阵第m₂个阵列圆环上第n₂个阵元的位置矢量

为：

其中，

n₂的最大取值

为随着圆环阵序号m₂变换。

圆锥阵对应的阵元位置坐标为：

则圆锥阵整体的位置矢量：

则整个异构阵的坐标为：

P＝[P₁，P₂]

其中，P∈C^3×N，N为几个共形阵的阵元数总和，对应的第n个阵元的位置矢量为：

其中，n∈1,2,…,N。

如果再加入其他共形阵，以列向量形式将坐标进行合并即可。

步骤2，利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换，得到机载三维异构阵的全局天线方向图；

由于三维异构阵面的非线性特性，即使每个阵元的辐射方向图具有相同极化，在阵列的全局坐标系下，由于载体曲面的弯曲，载体表面各阵元辐射方向图具有不同极化，亦即产生了交叉极化。根据交叉极化定义的极化基

和θ来计算极化阵列的辐射方向图。因为考虑极化因素，阵列天线在

方向的场强方向图可以表示为：

其中，

和

分别为

在

和θ方向的场分量，

和

分别表示表达式为

其中，Δφ_n表示第n个阵元的相位差，上标H表示共轭转置操作；

和w_nθ分别代表沿

方向极化分量加权和沿θ方向极化分量加权，它们的比值由第n个阵元姿态决定，不能独立调节。

接下来通过欧拉坐标旋转变换将局部坐标系下的天线方向图变换到全局直角坐标系下的天线方向图。

欧拉旋转变换是一种实现局部坐标和全局坐标之间变换的有效方向。不同坐标系、不同共形阵对应的欧拉旋转矩阵不同，本发明提到的异构阵要实现坐标变换需要进行三次欧拉旋转变换。下面说明一种变换情况，如图5所示，其旋转轴依次为z轴，y轴，x轴，欧拉旋转角分别为D_n、E_n和F_n，三次欧拉旋转矩阵表达式分别为：

沿z轴顺时针旋转：

沿y轴顺时针旋转：

沿x轴顺时针旋转：

则最终的变换矩阵为：R＝R_xR_yR_z

局部坐标(x,y,z)和全局坐标

之间的变换关系为：

不同的坐标系，不同共形阵的旋转情况都不一样。

如图6所示，以圆柱阵和圆锥阵为例，圆柱阵只需要进行一个沿y轴的旋转，如图6(a)所示，则第m个阵元对应的旋转角分别为：

D_n＝0

E_n＝π-θ₀+(m-1)β

F_n＝0；

如图6(b)所示，圆锥阵需要先进行一个y轴的旋转，然后沿x轴旋转，则第m个阵元对应的旋转角分别为：

D_n＝0

E_n＝π-θ₀+(m-1)β

F_n＝θ₀；

其中，θ₀为圆锥阵的锥角，m表示阵元序号，β表示相邻阵元之间的夹角；

首先需要知道阵元的局部方向图，然后进行旋转变换得到全局的方向图。阵元局部球坐标系下的方向图跟具体天线有关，一般选用圆形微带天线作为共形阵的阵元天线，其局部球坐标系下的第n个阵元的天线方向图表示为：

其中，

和

分别为局部球坐标下

和θ_l方向上的单位向量；

和

分别为第n个阵元局部方向图

在

和θ_l方向的场分量；J₂()和J₀()分别为2阶和0阶Bessel函数；

为入射信号方向在局部坐标系下的坐标，其与全局球坐标系下的坐标的转换关系为

实现步骤为：

(a1)将对应的全局球坐标

转换为全局直角坐标系下(x,y,z)的形式，转换关系即：

(b1)利用欧拉旋转变换将全局直角坐标系下的入射方向(x,y,z)变换为局部直角坐标系下的入射方向(x_l,y_l,z_l)，转换关系为：

(c1)由于阵元的辐射方向图一般由方位角和俯仰角所确定，因此需要将局部直角坐标(x_l,y_l,z_l)转换为局部球坐标

转换关系：

θ_l＝arccosz_l

将局部球坐标代入之前的微带天线的方向图函数就可得到局部坐标系下的阵元辐射方向图，然后再用上述欧拉变换的逆变换转换成全局球坐标系方向图，具体步骤如下：

(a2)利用球坐标系和直角坐标系的转换关系，得到天线方向图在局部直角坐标系下的分量，转换后的阵元方向图为：

转换关系：

其中，

分别为局部直角坐标下各正交方向上的单位向量。

(b2)将在局部直角坐标系下的方向图转换为全局直角坐标系下的方向图，表示成

使用欧拉旋转变换，即：

(c2)利用球坐标系和直角坐标系的转换关系将全局直角坐标系下天线方向图转换为全局球坐标系下的方向图

可表示为：

其中，

和

是全局球坐标系下方位向和俯仰向上的单位方向向量；

分别为局部直角坐标下各正交方向上的单位向量。

步骤3，获取机载三维异构阵的接收信号，并建立对应的过完备字典矩阵；通过所述过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题，并建立对应的压缩感知模型；

(3.1)接受信号获取：

假设有K个远场信号s_i以入射角

到达阵列，其中θ_i和

分别表示第i个入射信号的俯仰角和方位角。入射角为

的信号的单位方向矢量为

在得到阵元的位置矢量后，机载三维异构阵的阵元的空间相位差为：

对应的导向矢量为：

其中，

γ_n为第n个阵元对单位信号的响应；g_n为阵元方向图

g_nθ和

分别为

和

k_θ，

分别表示为信源来波极化状态，θ_ngk为g_n与

的夹角。

接收信号快拍数据模型可表示为：

各向量分别定义为：

X(t)＝[x₁(t),x₂(t),…,x_N(t)]^T

S(t)＝[s₁(t),s₂(t),…,s_K(t)]^T

n(t)＝[n₁(t),n₂(t),…,n_N(t)]^T

其中，N为机载三维异构阵的阵元总数，t为采样时间，长度即快拍数为N_t，

为机载三维异构阵的阵列流形矩阵，

为阵列入射信号，

为加性高斯白噪声。

(3.2)考虑到上述信号模型的DOA估计是一个二维非线性问题，本发明假设信号快拍数为N_t＝1，将阵元角度问题转变为稀疏表达问题，通过构建过完备字典矩阵建立压缩感知模型，通过稀疏恢复算法估计方位角和俯仰角。

将方位角度-俯仰角度组成的二维空间域沿方位角轴和俯仰角轴划分为

和N_θ个均匀格子，则对应的所有网格点的导向矢量组成的完备矩阵为：

均匀量化后的快拍数据模型可表示为：

其中，γ_i',j'为信号方向为

的接收信号幅度；

为高斯噪声；入射信号矩阵为

(3.3)在上述量化后的模型中，通过过完备字典矩阵Φ使得角度估计问题转化为γ的稀疏恢复问题。由于信号源个数远小于空间域的量化的网格点数，γ在空间域是稀疏的，可以构建成L1范数最小化的求解问题，其目标函数为：

步骤4，利用稀疏贝叶斯学习算法求解所述压缩感知模型，得到机载三维异构阵的DOA估计结果。

稀疏贝叶斯学习(SBL，Sparse Bayesian Learning)是在贝叶斯原理的基础上，加入最大后验概率估计和最大似然估计的一种稀疏恢复算法。稀疏贝叶斯在无需稀疏度先验信息和噪声干扰阈值信息的条件下，可以较为准确得到稀疏信号的估计值，本发明采用该算法来进行L1范数的求解。直接给出SBL步骤：

(4.1)为了便于迭代计算开始，需要初始化参数，该初始化过程并不会影响最后的计算结果。初始化噪声水平σ²＝0.001，初始化超参数α_i＝0.001，

设置收敛条件Δ＝0.001；

(4.2)根据下式计算稀疏信号后验概率的方差Σ和均值μ，并令μ^old＝μ；

∑＝(σ^-2Φ^TΦ+Λ)

μ＝σ^-2ΣΦ^TX

其中，

(4.3)根据最大似然估计方法计算超参

和噪声水平(σ²)^new的估计值：

其中，γ_i＝1-α_i∑_ii，∑_ii表示∑的对角线元素组成的向量，

(4.4)根据步骤(4.2)计算的

和(σ²)^new重新计算方差和均值，并记作Σ^new和μ^new；

(4.5)计算dμ＝abs(max(μ^new-μ^old))，并判断dμ是否小于等于Δ，若是，则停止迭代，μ^new即为重构信号γ；否则令μ^old＝μ^new，返回步骤(4.3)。

仿真实验

以下通过仿真实验验证本发明的有效性，以圆柱-圆锥阵的组成的异构阵为例，仿真参数如表所示：

表1系统仿真参数

目标1、2和3对应的估计结果如图7(a)、7(b)和7(c)所示，图7(a)估计的目标1的角度值(60°,100°)，(70°,110°)，幅度值为0.7dB和0.1dB；7(b)估计的角度值为(64.76°,105.2°)，(70°,110°)，幅度值为0.9999dB，1dB；7(c)估计的角度值为(64.76°,105.2°)，(74.76°,115.2°),幅度值为0.75dB，0.9dB，可以看出，三个目标幅度值的估计基本正确，目标2和目标3角度估计有一个小于网格精度0.476误差，这是由于网格划分的精度不够而产生的误差，在误差允许范围。综上，本发明方法能很好的进行机载三维异构阵的DOA估计。

虽然，本说明书中已经用一般性说明及具体实施方案对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (10)

1.基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，包括以下步骤：

步骤1，建立机载三维异构阵的阵列模型；

步骤2，利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换，得到机载三维异构阵的全局天线方向图；

步骤3，获取机载三维异构阵的接收信号，并建立对应的过完备字典矩阵；通过所述过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题，并建立对应的压缩感知模型；

步骤4，利用稀疏贝叶斯学习算法求解所述压缩感知模型，得到机载三维异构阵的DOA估计结果。

2.根据权利要求1所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述机载三维异构阵的阵列模型包括机头共形阵、机身圆柱共形阵、两侧机翼共形阵和机身圆台共形阵中的一种或多种。

3.根据权利要求2所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述建立机载三维异构阵的阵列模型具体为：

首先，对于圆柱-圆锥阵进行模型建立：以圆柱阵第1个圆环阵的中心点的正下方地面点为参考点即坐标原点，以载机飞行方向为y轴正方向，y轴正方向顺时针转90度为x轴正方向，竖直向上为z轴正方向，建立全局直角坐标系；

其次，载机速度矢量为v，飞行高度为H；圆柱阵部分：由M₁个等半径的圆环阵组成，每个圆环阵由关于圆心对称的N₁个天线阵元组成，阵元间距为d；M₁个圆环阵彼此平行，间距为d，且共同垂直于y轴；每个圆环阵半径为r；

圆锥阵部分：由M₂个半径按比例降低的圆环阵组成，圆锥顶部为单点阵，圆锥底层半径与圆柱底层半径相同，阵元数相等；以底层为第一层，依次排序，则圆锥阵的第m层圆环阵上的阵元个数为N_m；

则圆柱阵上第m₁层阵列环上第n₁个阵元的空间坐标

为：

其中，上标T为矩阵转置操作；

以圆环阵上z轴正向上的阵元为第1个阵元，阵元序号为逆时针排序，则圆柱阵对应的阵元位置的坐标为：

圆柱阵整体的位置矢量：

将第m₂个圆环阵的半径记为

阵元个数记为

底层半径阵元个数与圆柱阵一致，圆柱阵半径和阵元数的表达式为：

设第1个圆环阵离原点距离为D₂，以圆环阵上顶部阵元为第1个阵元，阵元序号为逆时针排序，则圆锥阵第m₂个阵列圆环上第n₂个阵元的位置矢量

为：

其中，

上标T表示矩阵的转置操作；

圆锥阵对应的阵元位置坐标为：

则整个圆锥阵的位置矢量为：

则整个三维异构阵的坐标为：

P＝[P₁，P₂]

其中，P∈C^3×N，N为整个三维异构阵的阵元数总和，对应的第n个阵元的位置矢量为：

其中，n∈1，2，...，N。

4.根据权利要求1所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述利用欧拉旋转变换，对机载三维异构阵的每个阵元进行局部直角坐标系到全局直角坐标系的转换，具体为：

(2.1)根据每个天线形式确定其局部球坐标系下的天线方向图函数；

(2.2)将入射信号方向的全局球坐标

转换为全局直角坐标系下的坐标(x，y，z)，转换关系为：

(2.3)利用欧拉旋转变换将全局直角坐标系下的入射方向(x，y，z)变换为局部直角坐标系下的入射方向(x_l，y_l，z_l)，转换关系为：

其中，R为欧拉旋转变换的变换矩阵；

(2.4)将局部直角坐标(x_l，y_l，z_l)转换为局部球坐标

转换关系为：

θ_l＝arccos z_l

(2.5)将步骤(2.4)的局部球坐标代入步骤(2.1)中的天线方向图函数，得到局部坐标系下的阵元辐射方向图；然后再用所述欧拉变换的逆变换将局部坐标系下的阵元辐射方向图转换成全局球坐标系下的天线方向图即可。

5.根据权利要求4所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，当阵元天线为圆形微带天线时，其局部球坐标系下的第n个阵元的天线方向图函数表示为：

其中，

和

分别为局部球坐标下方位角

和俯仰角θ_l方向上的单位向量；

和

分别为第n个阵元局部方向图

在

和θ_l方向的场分量；J₂()和J₀()分别为2阶和0阶Bessel函数；

为入射信号方向在局部坐标系下的坐标；λ为波长，d为阵元间距，j为虚数单位。

6.根据权利要求4所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述欧拉旋转变换的变换矩阵的获取过程：

首先，根据三维异构阵的阵列模型，确定实现坐标变换需要进行的欧拉旋转变换；

然后，一般性地，设旋转轴依次为z轴、y轴、x轴，对应的欧拉旋转角分别为D_n、E_n、F_n，则每次旋转对应的欧拉旋转矩阵表达式分别为：

沿z轴顺时针旋转：

沿y轴顺时针旋转：

沿x轴顺时针旋转：

则该欧拉旋转变换的变换矩阵为：R＝R_xR_vR_z。

7.根据权利要求6所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，当异构阵为圆柱阵时，对应的欧拉变换为沿y轴旋转，则圆柱阵某一圆环阵上的第m₁个阵元对应的旋转角分别为：

D_n＝0

E_n＝π-θ₀+(m₁-1)β

F_n＝0；

当异构阵为圆锥阵时，对应的欧拉变换为先进行一个y轴的旋转，然后沿x轴旋转，则圆锥阵某一圆环阵上的第m₂个阵元对应的旋转角分别为：

D_n＝0

E_n＝π-θ₀+(m₂-1)β

F_n＝θ₀；

其中，θ₀为圆锥阵的锥角，β为相邻的阵元夹角。

8.根据权利要求1所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述获取机载三维异构阵的接收信号，具体为：

设有K个远场信号以不同入射角到达阵列，其中入射角包含方位角和俯仰角；入射角为

的入射信号的单位方向矢量为：

则机载三维异构阵的阵元的空间相位差为：

对应的导向矢量为：

其中，

为第n个阵元的位置矢量，

γ_n为第n个阵元对单位信号的响应；

接收信号的快拍数据模型表示为：

其中，t为采样时间，

N_t为快拍数，N为机载三维异构阵的阵元总数，

为机载三维异构阵的阵列流形矩阵，

为阵列入射信号，

为加性高斯白噪声。

9.根据权利要求7所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述通过所述过完备字典矩阵将目标的角度估计问题转化为入射信号矩阵的稀疏恢复问题，具体过程为：

(3.1)将方位角度-俯仰角度组成的二维空间域沿方位角轴和俯仰角轴划分为

个均匀网格，则对应所有网格点的导向矢量组成的过完备字典矩阵为：

则均匀量化后的快拍数据模型为：

其中，γ_i′，j′为信号方向为

的接收信号幅度；

为高斯噪声；入射信号矩阵为

(3.2)由于信号源个数远小于空间域的量化网格点数，Υ在空间域是稀疏的，构建成L1范数最小化的求解问题，其目标函数为：

其中，|| ||₁为求L1范数，|| ||₂为求L2范数，ε为常数。

10.根据权利要求9所述的基于稀疏贝叶斯学习的机载三维异构阵DOA估计方法，其特征在于，所述利用稀疏贝叶斯学习算法求解所述压缩感知模型，具体过程为：

(4.1)初始化噪声水平σ²＝0.001，初始化超参数

设置收敛条件Δ＝0.001；

(4.2)根据下式计算稀疏信号后验概率的方差∑和均值μ，并令μ^old＝μ；

∑＝(σ^-2Φ^TΦ+Λ)

μ＝σ^-2∑Φ^TX

其中，

(4.3)根据最大似然估计方法计算超参

和噪声水平(σ²)^new的估计值：

其中，γ_i＝1-α_i∑_ii，∑_ii表示∑的对角线元素组成的向量，

(4.4)根据步骤(4.2)计算的

和(σ²)^new重新计算方差和均值，并记作∑^new和μ^new；

(4.5)计算dμ＝abs(max(μ^new-μ^old))，并判断dμ是否小于等于Δ，若是，则停止迭代，μ^new即为重构信号Υ；否则令μ^old＝μ^new，返回步骤(4.3)。

Images