基于多输出支持向量回归机的近场声源参数估计方法与流程

本发明属于阵列信号处理技术领域，尤其涉及一种多输出支持向量回归机近场声源参数估计方法。

背景技术

到达角(directionofarrival，doa)估计是阵列信号处理领域的一个重要研究方向。传统的近场源doa估计方法有两步music方法、广义esprit方法以及其它的改进方法等等。两步music方法利用信号子空间和噪声子空间的正交性来实现目标的参数估计，但传统两步music方法应用于近场源参数估计时存在一些不足之处，谱峰搜索计算量大、当信源为相干信源时会出现秩亏以致无法分辨、在低信噪比情形下对信源角度相隔较近时的估计精度迅速下降等。广义esprit方法对于相干信号也存在秩亏损以致无法分辨。

支持向量机(supportvectormachine，svm)是基于统计学习理论的一种机器学习方法，它能够有效的解决高维度、小样本、非线性等问题，该方法本身在训练过程中具有全局最优解、泛化能力强、收敛速度较快和小样本条件下预测性能较好等优点。并且支持向量机模型训练过程可以离线进行，而训练好的支持向量机模型具有优异的数据处理速度，可以对阵列接收数据实现快速处理，有利于实际工程中运用。近年来，多输入单输出的支持向量回归理论日渐成熟，但是由于到达角估计往往定位目标具有多个，是多输出问题，如果对于每一个定位目标都构建一个支持向量回归机(supportvectorregression，svr)模型，这就需要构建多个支持向量回归机模型，计算量大而且不能体现各个输出变量之间的关系。目前对多输出的支持向量回归模型研究还较少，将多输出的支持向量回归机模型运用到近场声源的参数估计中来还鲜见报道。

针对两步music方法、广义esprit方法存在的不足之处诸如非相干信号低信噪比下小间隔角度分辨率差以及无法处理相干信源等问题，本发明的多输出支持向量回归机的近场声源参数估计方法将接收数据得到的特征向量作为输入，定位目标的到达角作为输出，以支持向量回归机逼近得到输入与输出两者的非线性映射关系，然后利用得到的多输出支持向量回归机模型来估计目标的角度值，利用得到的角度值和music方法来估计近场声源的距离参数，这对于相干信号以及低信噪比下的非相干信号doa估计都有着良好的效果。

技术实现要素：

本发明的目的是提供一种基于多输出支持向量回归机的近场声源参数估计方法。

为了实现上述目的，本发明采取如下的技术解决方案：

基于多输出支持向量回归机的近场声源参数估计方法，采用普通的声压传感器阵列接收k个窄带、非高斯的平稳近场声源信号。阵列通过以下方式构成：在空间任意选取一点作为坐标轴的原点位置o，从左至右经过该原点的水平线为x轴，垂直于该水平线的直线为z轴，即假设声源从xoz平面入射，第k个入射声源于z轴的夹角为θk，θk的取值范围为[-π/2，π/2]，在坐标原点处两边分别以d＝λmin/4等间隔放置p个阵元，λmin为入射声源中最小波长，阵元总个数为m＝2p+1。对阵元从左到右依次标记为[-p，-p+1，...，m，...，p-1，p]，m取值范围为[-p，p]。

基于多输出支持向量回归机的近场声源参数估计方法步骤如下：

步骤一、利用阵元数为m＝2p+1的均匀对称阵列作为接收阵列，接收k个窄带、非高斯、平稳近场声源信号在训练区间为[-θ0，+θ0]中产生的w组接收数据y＝[y1，y2，...，yw]和声源信号到达角集合θ＝[θ1；θ2；…；θw；…；θw]；声源信号距离r＝[r1，r2…，rl，…，rl]；

第一组声源信号到达角为θ1＝[θ1¹，θ1²，...，θ1^l，...，θ1^l]入射到接收阵列上，经过t次采样后得到维度为m×t的第一组接收数据y1；相邻声源间隔为[δθ1¹²，δθ1²³，...，δθ1^(l-1)l，...，δθ1^(l-1)l]，其中δθ1^(l-1)l表示第一组声源信号中第l和第l-1个声源信号的间隔，θ1^l是第一组声源中第l个信号的到达角，rl是第l个信号的到坐标原点的距离，在声源的距离和相邻声源信号间隔不变的前提下，使整体旋转角度δφ，此时的声源信号为θ2＝[θ2¹，θ2²，...，θ2^l，...，θ2^l]，此时阵列接收信号经过t次采样后得到维度为m×t的第二组接收数据y2，按照这样的方式得到接收数据集合为y＝[y1，y2，...，yw，...，yw]；

步骤二、对接收数据集合y中的每一组接收数据，利用阵列对称位置上接收数据作协方差之后得到接收数据集合y相对应的协方差矩阵γ＝[γ1，γ2，...，γw，...，γw]；

对集合y中的第w组接收数据yw，构造一个矩阵γw，yw，p表示yw中第p个阵元的接收数据，yw，-p表示yw中第-p个阵元的接收数据，(·)^*表示取共轭，其中-p≤p≤p，δ(·)为狄利克雷函数，代表第l个声源信号的方差，为噪声方差，其中表示第w组接收数据中的l个信号的到达角，λl表示第l个测试信号的波长，γw矩阵：

对于接收数据集合y中的其它组接收数据也做同样的处理得到相对应的协方差矩阵，最后得到w组接收数据协方差矩阵γ＝[γ1，γ2，...，γw，...，γw]。

步骤三、对数据协方差矩阵γ进行矩阵运算处理和归一化处理可得训练数据集的特征矩阵h＝[h1，h2，…，hw，…，hw]；

选取协方差矩阵中的元素γw(p，-p)，得到特征向量pw＝[γw(0，0)，γw(-1，1)，...，γw(p，-p)，…，γw(-p，p)]，接收数据协方差矩阵γ＝[γ1，γ2，...，γw，...，γw]中w协方差矩阵都按照这样的方式构造得到特征向量集合p＝[p1，p2，...，pw，...，pw]，对pw归一化处理可得hw＝pw/||pw||，γ1，γ2，...，γw，...，γw都做相同的处理可得训练数据的特征向量矩阵h＝[h1，h2，...，hw，...，hw]；

步骤四、根据得到的训练数据的特征向量矩阵h＝[h1，h2，...，hw，...，hw]和信号源集合θ＝[θ1；θ2；…；θw；…；θw]训练得到多输出支持向量最优回归参数βomp；

1)将训练数据的特征向量矩阵h＝[h1，h2，…，hi，…，hw]代入高斯核函数得到基于训练样本的核矩阵k，其中kij＝exp(-||hi-hj||²/(2σ^２))，σ为高斯核函数的系数，kij表示核矩阵k的第i行第j列的元素，hi和hj分别表示数据特征向量集合h的第i和第j个特征向量；

2)初始化回归参数矩阵β⁰为一个w×t的全零矩阵、常数惩罚因子c和管壁误差ε；初始误差矩阵e⁰＝θ，根据初始误差矩阵e⁰得到初始拉格朗日参数因子矩阵和初始拉格朗日参数系数矩阵表示e⁰的第w行，其中寻找u⁰矩阵中满足的位置并且按照其在u⁰中的先后顺序存放于支持向量矩阵和支持向量位置矩阵根据矩阵和位置矩阵计算出损失函数矩阵

3)第k步循环，步长因子η^k＝1，计算出误差矩阵e^k＝θ-kβk^-1，根据误差矩阵表示e^k的第w行，得到拉格朗日参数因子矩阵和拉格朗日参数系数矩阵其中寻找u^k矩阵中满足的位置从而得到支持向量位置矩阵和对应位置的拉格朗日参数因子矩阵根据矩阵、位置矩阵和管壁误差ε计算出损失函数从而得到损失函数矩阵

特别地，u^k矩阵中满足的矩阵元素是则支持向量矩阵支持向量位置矩阵

根据梯度下降法计算出回归系数下降方向矩阵p^k，p^k＝inv(g^k)θ，inv表示对矩阵求逆，矩阵矩阵p^k的每一列表示回归系数的下降方向，表示只保留核矩阵k支持向量位置上的元素其他元素全部为零的矩阵，表示只保留拉格朗日参数系数矩阵α^k支持向量位置上的元素其他元素全部为零的矩阵，β^k＝β^k-1+p^k再根据回归参数矩阵β^k、核矩阵k和损失函数矩阵构造出迭代前的优化目标函数π＝diag((β^k)^hkβ^k)，πt，t表示矩阵π的第t行第t列元素，t＝1，…，t，diag(·)是取矩阵的对角线元素运算，其中

4)比较l^k和l^k-1的大小，如果l^k＞l^k-1，η^k＝γη^k，γ＜1，η^k＝γη^k表示将η^k缩小为原来的γ倍，根据矩阵p^k、γ＜1和β^k-1得到本次迭代β^k＝γp^k+(1-γ)β^k-1，重复步骤四中的3)；如果l^k＜l^k-1跳到步骤四中的5)；

5)验证是否满足则k＝k+1继续步骤3)其中k＝k+1表示将k增加1；如果满足则跳到步骤四中的6)；

6)如果满足则训练过程结束；此时的β^k＝βomp为多输出支持向量最优回归参数，εmin为误差门限；

步骤五、将训练数据特征矩阵h＝[h1，h２，…，hw，…，hw]和测试数据特征矢量he带入高斯核函数得到测试核矩阵ke，利用训练好的近场声源多输出支持向量最优回归参数βomp和ke估计声源的角度

阵列接收的测试数据ye，利用阵列对称位置上阵元的接收数据作协方差之后得到与测试数据ye相对应的协方差矩阵re，再按照步骤三中的方法由re得到归一化特征矢量he，将训练数据特征矩阵h＝[h1，h２，…，hw，…，hw]和测试数据特征矢量he代入高斯核函数得到测试核矩阵ke，核矩阵ke与步骤四得到的最优回归参数矩阵βomp相乘得到声源的角度的估计值其中ke(w)＝exp(-||hw-he||²/(2σ²))，

其中-p≤p≤p，δ(·)为狄利克雷函数，代表第l个测试声源信号的方差，σ²为测试噪声协方差，其中θel表示第l个测试的到达角，λel表示第l个测试的波长，利用γe中数据元素可以构造出类似远场的接收数据协方差矩阵；取数据协方差矩阵re中的元素构成特征向量pe＝[γe(0，0)，γe(-1，1)，…，γe(p，-p)，…，γw(-p，p)]，对pe归一化处理可得he＝pe/||pe||；

步骤六、求阵列接收的测试数据ye的数据相关矩阵通过数据相关矩阵的特征分解得到噪声子空间eyn，对每一个到达角通过music谱函数进行距离项的谱峰搜索估计其对应的距离得到距离的估计

根据步骤一中的阵列结构和步骤五估计得到的到达角信息构造music谱函数通过谱峰搜索估计测试信号的距离对每一个到达角通过music谱函数进行距离项的谱峰搜索估计其对应的距离本方法无需额外参数配对运算；其中为music谱函数，为将到达角信息代入后的阵列导向矢量，它是距离r的函数；

前述步骤中的k＝1，2，…，kx表示循环次数，l＝1，......，l表示信号个数，w＝1，2，......，w表示样本数，p＝-p，...，p表示阵元的标号，i＝1，2，......w表示元素在矩阵中的位置，j＝1，2，......，w表示元素在矩阵中的位置，f＝1，2，.....，f表示支持向量的个数。

本发明利用支持向量回归机方法进行阵列参数估计，通过对称阵元数据对消处理运算，将距离项对消掉，先由支持向量回归机方法回归出到达角信息，再利用到达角信息进行距离项参数估计；本发明方法在角度估计过程是通过数据的学习得到到达角的信息，对于近场相干和非相干信号都可以估计得到到达角不需要解相干的繁琐过程。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中需要使用的附图做简单介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明实施简化传感器阵列的示意图；

图2为本发明的流程图；

图3为两步music方法到达角估计均方根误差随信噪比的变化曲线图；

图4为本发明方法到达角估计均方根误差随信噪比变化曲线图；

图5为本发明方法到达角估计成功概率随信噪比的变化曲线图；

图6为本发明方法到达角估计散布图；

图7为本发明方法到达角估计成功概率图；

图8为本发明方法距离参数估计的归一化功率谱；

具体实施方式

为了让本发明的上述和其它目的、特征及优点能更明显，下文特举本发明实施例，并配合所附图示，做详细说明如下。

本发明的目的是提供一种基于多输出支持向量回归机的近场声源参数估计方法。

为了实现上述目的，本发明采取如下的技术解决方案：

基于多输出支持向量回归机的近场声源参数估计方法，采用普通的标量传感器阵列接收k个窄带、非高斯的平稳近场声源，所述数据接收阵列模型通过以下方式构成：在空间任意选取一点作为坐标轴的原点位置o，从左至右经过该原点的水平线为x轴，垂直于该水平线的直线为z轴，即假设声源从xoz平面入射，第l个入射声源于z轴的夹角为θl，θl的取值范围为[-π/2，π/2]，在坐标原点处两边分别以d＝λmin/4等间隔放置p个阵元，λmin为入射声源中最小波长，阵元总个数为m＝2p+1，对阵元从左到右依次标记为[-p，-p+1，...，...，m，...，p-1，p]，m取值范围为[-p，p]；

基于多输出支持向量回归机的近场声源参数估计方法步骤如下：

步骤一、利用阵元数为m＝2p+1的均匀对称阵列作为接收阵列，接收k个窄带、非高斯、平稳近场声源信号在训练区间为[-θ0，+θ0]中产生的w组接收数据y＝[y1，y2，...，yw]；声源信号到达角集合θ＝[θ1；θ2；…；θw；…；θw]；声源信号距离r＝[r1，r2…，rl，…，rl]；

第一组声源信号到达角为θ1＝[θ1¹，θ1²，...，θ1^l，...，θ1^l]入射到接收阵列上，经过t次采样后得到维度为m×t的第一组接收数据y1；相邻声源间隔为[δθ1¹²，δθ1²³，...，δθ1^(l-1)l，...，δθ1^(l-1)l]，其中δθ1^(l-1)l表示第一组声源信号中第l和第l-1个声源信号的间隔，θ1^l是第一组声源中第l个信号的到达角，rl是第l个信号的到坐标原点的距离，在声源的距离和相邻声源信号间隔不变的前提下，使整体旋转角度δφ，此时的声源信号为θ2＝[θ2¹，θ2²，...，θ2^l，...，θ2^l]，此时阵列接收信号经过t次采样后得到维度为m×t的第二组接收数据y2，按照这样的方式得到接收数据集合为y＝[y1，y2，...，yw，...，yw]；

步骤二、对接收数据集合y中的每一组接收数据，利用阵列对称位置上接收数据作协方差之后得到接收数据集合y相对应的协方差矩阵γ＝[γ1，γ2，...，γw，...，γw]；

对集合y中的第w组接收数据yw，构造一个矩阵γw，yw，p表示yw中第p个阵元的接收数据，yw，-p表示yw中第-p个阵元的接收数据，(·)^*表示取共轭，其中-p≤p≤p，δ(·)为狄利克雷函数，代表第l个声源信号的方差，为噪声方差，其中表示第w组接收数据中的l个信号的到达角，λl表示第l个测试信号的波长，γw矩阵：

对于接收数据集合y中的其它组接收数据也做同样的处理得到相对应的协方差矩阵，最后得到w组接收数据协方差矩阵γ＝[γ1，γ2，...，γw，...，γw]。

步骤三、对数据协方差矩阵γ进行矩阵运算处理和归一化处理可得训练数据集的特征矩阵h＝[h1，h2，…，hw，…，hw]；

选取协方差矩阵中的元素γw(p，-p)，得到特征向量pw＝[γw(0，0)，γw(-1，1)，...，γw(p，-p)，…，γw(-p，p)]，接收数据协方差矩阵γ＝[γ1，γ2，...，γw，...，γw]中w协方差矩阵都按照这样的方式构造得到特征向量集合p＝[p1，p2，...，pw，...，pw]，对pw归一化处理可得hw＝pw/||pw||，γ1，γ2，...，γw，...，γw都做相同的处理可得训练数据的特征向量矩阵h＝[h1，h2，...，hw，...，hw]；

步骤四、根据得到的训练数据的特征向量矩阵h＝[h1，h2，...，hw，...，hw]和信号源集合θ＝[θ1；θ2；…；θw；…；θw]训练得到多输出支持向量最优回归参数βomp；

1)将训练数据的特征向量矩阵h＝[h1，h2，…，hi，…，hw]代入高斯核函数得到基于训练样本的核矩阵k，其中kij＝exp(-||hi-hj||²/(2σ²))，σ为高斯核函数的系数，kij表示核矩阵k的第i行第j列的元素，hi和hj分别表示数据特征向量集合h的第i和第j个特征向量；

2)初始化回归参数矩阵β⁰为一个w×t的全零矩阵、常数惩罚因子c和管壁误差ε；初始误差矩阵e⁰＝θ，根据初始误差矩阵e⁰得到初始拉格朗日参数因子矩阵和初始拉格朗日参数系数矩阵表示e⁰的第w行，其中寻找u⁰矩阵中满足的位置并且按照其在u⁰中的先后顺序存放于支持向量矩阵和支持向量位置矩阵根据矩阵和位置矩阵计算出损失函数矩阵

3)第k步循环，步长因子η^k＝1，计算出误差矩阵e^k＝θ-kβ^k-1，根据误差矩阵表示e^k的第w行，得到拉格朗日参数因子矩阵和拉格朗日参数系数矩阵其中寻找u^k矩阵中满足的位置从而得到支持向量位置矩阵和对应位置的拉格朗日参数因子矩阵根据矩阵、位置矩阵和管壁误差ε计算出损失函数从而得到损失函数矩阵

特别地，u^k矩阵中满足的矩阵元素是则支持向量矩阵支持向量位置矩阵

根据梯度下降法计算出回归系数下降方向矩阵p^k，p^k＝inv(g^k)θ，inv表示对矩阵求逆，矩阵矩阵p^k的每一列表示回归系数的下降方向，表示只保留核矩阵k支持向量位置上的元素其他元素全部为零的矩阵，表示只保留拉格朗日参数系数矩阵α^k支持向量位置上的元素其他元素全部为零的矩阵，β^k＝β^k-1+p^k再根据回归参数矩阵β^k、核矩阵k和损失函数矩阵构造出迭代前的优化目标函数π＝diag((β^k)^hkβ^k)，πt，t表示矩阵π的第t行第t列元素，t＝1，…，t，diag(·)是取矩阵的对角线元素运算，其中

4)比较l^k和l^k-1的大小，如果l^k＞l^k-1，η^k＝γη^k，γ＜1，η^k＝γη^k表示将η^k缩小为原来的γ倍，根据矩阵p^k、γ＜1和β^k-1得到本次迭代β^k＝γp^k+(1-γ)β^k-1，重复步骤四中的3)；如果l^k＜l^k-1跳到步骤四中的5)；

5)验证是否满足则k＝k+1继续步骤3)其中k＝k+1表示将k增加1；如果满足则跳到步骤四中的6)；

6)如果满足则训练过程结束；此时的β^k＝βomp为多输出支持向量最优回归参数，εmin为误差门限；

步骤五、将训练数据特征矩阵h＝[h1，h２，…，hw，…，hw]和测试数据特征矢量he带入高斯核函数得到测试核矩阵ke，利用训练好的近场声源多输出支持向量最优回归参数βomp和ke估计声源的角度

阵列接收的测试数据ye，利用阵列对称位置上阵元的接收数据作协方差之后得到与测试数据ye相对应的协方差矩阵re，再按照步骤三中的方法由re得到归一化特征矢量he，将训练数据特征矩阵h＝[h1，h２，…，hw，…，hw]和测试数据特征矢量he代入高斯核函数得到测试核矩阵ke，核矩阵ke与步骤四得到的最优回归参数矩阵βomp相乘得到声源的角度的估计值其中ke(w)＝exp(-||hw-he||²/(2σ²))，

其中-p≤p≤p，δ(.)为狄利克雷函数，代表第l个测试声源信号的方差，σ²为测试噪声协方差，其中θel表示第l个测试的到达角，λel表示第l个测试的波长，利用γe中数据元素可以构造出类似远场的接收数据协方差矩阵；取数据协方差矩阵re中的元素构成特征向量pe＝[γe(0，0)，γe(-1，1)，...，γe(p，-p)，…，γw(-p，p)]，对pe归一化处理可得he＝pe/||pe||；

步骤六、求阵列接收的测试数据ye，的数据相关矩阵通过数据相关矩阵的特征分解得到噪声子空间eyn，对每一个到达角通过music谱函数进行距离项的谱峰搜索估计其对应的距离得到距离的估计

根据步骤一中的阵列结构和步骤五估计得到的到达角信息构造music谱函数通过谱峰搜索估计测试信号的距离对每一个到达角通过music谱函数进行距离项的谱峰搜索估计其对应的距离本方法无需额外参数配对运算；其中为music谱函数，为将到达角信息代入后的阵列导向矢量，它是距离r的函数；

前述步骤中的k＝1，2，…，kx表示循环次数，l＝1，......，l表示信号个数，w＝1，2，......，w表示样本数，p＝-p，...，p表示阵元的标号，i＝1，2，......w表示元素在矩阵中的位置，j＝1，2，......，w表示元素在矩阵中的位置，f＝1，2，.....，f表示支持向量的个数。

阵的行数，j＝1，2，......，w表示矩阵的列数，f＝1，2，.....，f表示支持向量的个数。

本发明利用支持向量回归机方法进行阵列参数估计，通过对称阵元数据对销处理运算，将距离项对销掉，先由支持向量回归机方法回归出到达角信息，再利用到达角信息进行距离项参数估计；本发明方法在角度估计过程中通过模型参数回归进行到达角估计，可以直接估计相干信号的到达角不需要解相干处理，过程简单。

本发明的效果可以通过以下的仿真结果进一步说明：

两个近场、窄带、非高斯平稳声源信号入射到图1所示的对称线性传感器阵列上，该接收阵列由9个阵元组成，阵元间隔为d＝λmin/4，相干时信号频率设为[fs/8，fs/8]，非相干时信号频率设为[fs/8，fs/12]，信号距离为[3λmin，6min]，快拍数为100，噪声为高斯白噪声。训练过程中对角度间隔为δθ1¹²＝5°、δθ2¹²＝10°、δθ3¹²＝15°、δθ4¹²＝20°四种间隔进行训练，训练角度的区间设在[-90°，+90°]；图3为两步music的近场到达角和距离估计方法，从图中可以看出music方法对相干源失效，由于music方法是利用信号子空间和噪声子空间的正交性来估计参数信息的，当信号源相干时信号子空间会出现秩亏损，不能准确的估计出信号的参数，如果要估计出信号源的参数信息需要额外的解相干处理，图4和图5分别为本发明方法(多输出支持向量回归机方法，简称msvr方法)的到达角估计均方根误差和成功概率，说明本发明方法能够有效的估计相干信号的到达角；

图6和图7验证两个测试角度为10°和18.5°时，两角度间隔为8.5°不在训练间隔5°、10°和15°之内时的估计性能，由图6可知对于角度为10°的信号估计值与真实值偏离了0.5°；角度为18.5°的信号估计值与真实值基本重合，这说明本发明方法(msvr方法)训练得到的近场声源msvr模型具有较好的泛化能力，由图7可知在信噪比大于5db时可以很高的成功概率实现信号估计；图8给出了本发明方法(msvr方法)的距离估计结果，从图上可以看出距离的估计基本是准确的，说明了本发明方法的有效性；

以上所述，仅是本发明的较佳实施例而已，并非对本发明做任何形式上的限制，虽然本发明已以较佳实施例揭露如上，然而并非用以限定本发明，任何熟悉本专业的技术人员，在不脱离本发明技术方案范围内，当可利用上述揭示的技术内容做出些许更动或修饰为等同变化的等效实施例，但凡是未脱离本发明技术方案的内容，依据本发明的技术实质对以上实施例所作的任何简单修改、等同变化与修饰，均仍属于本发明技术方案的范围内。