一种基于玉米粒碰撞声信号多域融合的PSO‑SVM优化方法与流程

本发明属于碰撞声检测领域，具体涉及一种基于玉米粒碰撞声信号多域融合的pso-svm优化方法。

背景技术：

粮食作为人类生存的必需品，人类的生活处处离不开粮食。但粮食产后损失是一个世界性难题，根据统计调查，粮食产后损失率在发达国家一般为3％，在发展中国家一般为5％，我国粮食产后损失率达到8％。每年产后粮食损失量为500亿斤，这相当于1亿多亩地生产的粮食被浪费。粮食产后损失中，储存环节损失所占比重最高，损失比例达到40.3％。这是由于粮食储存环节，极易受到虫害和霉菌的危害。这些危害不但会降低粮食的产量，而且也会降低粮食的质量。人类长期食用这类粮食，不但会出现营养不良，甚至会导致诱发各种疾病。在农产品中，由于玉米自身的特性更易受到虫类和霉菌的危害。玉米的胚占整个玉米颗粒体积的1/3左右，占重量的10％～20％。玉米的胚含有较多的蛋白质和大量可溶性糖，由于胚部营养成分较多，容易被微生物寄生和遭到虫类危害，导致玉米的带菌量远大于其他禾谷类粮食。因此，在玉米的储藏过程中预防虫害和霉菌感染的检测工作具有重要的现实意义。

传统的储粮检测方法常见有二氧化碳检测、x射线成像法、尿酸检测、机器视觉法等，这些方法有一定的优势，但也存在局限性。目前，利用碰撞声进行检测的方法因其效率高和方便等诸优势已经引起有关研究人员的注意。

技术实现要素：

针对现有技术中的不足，本发明提出一种基于玉米粒碰撞声信号多域融合的pso-svm优化方法，目的在于通过提取玉米粒碰撞声信号时域特征、频域特征和希尔伯特域特征，且将所述特征送入到粒子群优化的支持向量机进行识别和检测，以提高玉米完好粒、虫蛀粒和霉变粒的识别率。

一种基于玉米粒碰撞声信号多域融合的pso-svm优化方法，具体包括：

s1：对玉米粒碰撞声信号进行eemd分解；

s2：提取玉米粒碰撞声信号的时域特征；

所述步骤s2包括：

s21：计算平均振幅变化，所述平均振幅变化用于代表信号衰减，计算公式为：

s22：统计wilson振幅值，当碰撞声信号振幅的变化超过给定阈值时，所述wilson振幅的值就加1，否则不变，计算公式为：

s23：计算峰峰值，所述峰峰值为在信号的时间序列中最大值和最小值之间的差异，用于描述信号幅值总体的变化范围，计算公式为：

pk-pk＝max(ci，k)-min(ci，k)；

s24：计算样本熵，所述样本熵的计算公式为：

当n为有限值时，计算公式为：

以上各式中，n是信号采样点数，ci，k为第k次试验的第i阶imf分量；

s3：提取玉米粒碰撞声信号的频域特征；

所述步骤s3包括：

s31：计算频率均值，对原始碰撞声信号做fft变换，提取fft变换幅度均值作为所述频率均值，计算公式为：

其中，m代表采样点个数，x(k)表示经傅里叶变换之后的幅度；

s32：计算功率谱均方根，所述功率谱均方根用于反映功率谱能量的离散程度，计算公式为：

其中，m代表采样点个数，p代表功率谱的幅度；

s33：计算频带方差，所述频带方差用于表示信号在不同频带的能量方差，计算公式为：

其中，m代表采样点个数，x(k)表示经fft变换之后的幅度；

s4：提取玉米粒碰撞声信号的希尔伯特域特征；

所述步骤s4包括：

s41：计算低频带能量和高频带能量，所述高频带频率范围为10～15khz，计算公式为：

所述低频带频率范围为5～10khz，计算公式为：

ah，i为高频带瞬时强度，al，i为低频带瞬时强度，n为信号采样点数；

s42：计算包络平均值，具体步骤又包括：

s421：对前几个imf做hilbert变换定义为：

其中，h[]为hilbert变换；

s422：分别以c(i)和hilbert变换h(ci(t))为实部和虚部，由此可以将原始信号解析表示如下：

z(t)＝c(i)+jh[ci(t)]；

s423：信号的包络b(t)就是上述解析信号的幅值函数：

s424：包络谱则可以由对包络函数b(t)进行傅里叶变换求出：

a(f)＝f(b(t))；

s5：利用pso算法优化svm参数，构建pso-svm分类器，用于对所提取的玉米粒碰撞声信号的时域特征、频域特征和希尔伯特域特征进行分类识别，以获得更高的玉米颗粒完好粒、霉变粒和虫蛀粒的识别率；

所述步骤s5包括：

s51：初始化粒子群，随机产生粒子的位置和速度；

s52：利用碰撞声信号样本，训练svm分类器并计算适应度值，采用交叉验

证的方法计算适应度值；

s53：根据适应度值更新个体极值和全局最优值；

s54：更新当前粒子的速度和位置；

s55：当达到迭代次数或满足适应度条件时，则终止迭代；否则返回步骤s52。

发明效果：

根据本发明具有如下效果：多域融合可以提高碰撞声信号的特征分类效果，有助于识别玉米完好粒，霉变粒和虫蛀粒。

附图说明

图1(a)-图1(c)是本发明实施例提供的玉米完好粒、虫蛀粒和霉变粒三类碰撞声信号时域图；

其中，图1(a)为完好粒，图1(b)为虫蛀粒，图1(c)为霉变粒；

图2(a)-图2(d)是本发明实施例提供的三类玉米粒碰撞声信号的平均振幅变化散点图；

其中，图2(a)为imf1的平均振幅变化散点图，图2(b)为imf2的平均振幅变化散点图，图2(c)为imf3的平均振幅变化散点图，图2(d)为imf4的平均振幅变化散点图；

图3是本发明实施例提供的三类碰撞声信号的wilson振幅均值图；

图4是本发明实施例提供的三类碰撞声信号的峰峰值的均值图；

图5是本发明实施例提供的三类碰撞声信号的样本熵总数柱状图；

图6(a)-图6(c)是本发明实施例提供的玉米完好粒、虫蛀粒和霉变粒三类碰撞声信号的频域图；

其中，图6(a)为完好粒，图6(b)为虫蛀粒，图6(c)为霉变粒；

图7是本发明实施例提供的三类碰撞声信号的频率均值图；

图8是本发明实施例提供的三类碰撞声信号的功率谱均方根幅值图；

图9是本发明实施例提供的三类碰撞声信号的频带方差图；

图10是本发明实施例提供的三类碰撞声信号的包络平均值对比图；

图11是本发明实施例提供的pso优化svm参数的流程图；

图12是本发明实施例提供的基于多域融合的玉米受损颗粒检测系统流程图；

图13是本发明实施例提供的pso适应度曲线图。

具体实施方式

下面结合附图和实施例对本发明作进一步的详细描述。

本发明的基本思想在于提取玉米粒碰撞声信号的时域特征、频域特征和希尔伯特域特征，利用pso优化svm分类器，对上述多域特征进行分类识别。

本发明实施例提供一种基于玉米粒碰撞声信号多域融合的pso-svm优化方法，具体技术方案如下：

1、对玉米粒碰撞声信号进行eemd分解

eemd又称总体经验模态分解，该方法是将白噪声添加到原始信号中作为一个新的信号序列，对此信号序列进行emd得到imf，通过多次添加白噪声获取imf，然后对所得imf进行平均处理，经过多次平均运算之后，加入原始信号的白噪声可以相互抵消。

2.提取玉米粒碰撞声信号时域特征

如图1(a)-图1(c)所示的时域图中，我们看出玉米完好粒、虫蛀粒和霉变粒三种类型的碰撞声信号的幅值的波动范围不同，各类信号的衰减程度有区别。提取四个时域特征，即平均振幅变化(aac)、wilson振幅(wamp)、峰峰值(pk-pk)、样本熵。

2.1计算平均振幅变化(aac)

由时域图可以看出，玉米完好粒、虫蛀粒和霉变粒三类碰撞声信号的衰减程度有差别。用振幅变化代表信号的衰减，计算公式如下：

其中，n取值2048。

图2(a)-图2(d)展示了三类玉米粒碰撞声信号的平均振幅变化散点图，星号代表的是完好粒，五角星代表的是霉变粒，圆圈代表的是虫蛀粒，图2(a)、图2(b)、图2(c)、图2(d)分别代表了imf1、imf2、imf3、imf4的平均振幅变化散点图，其中，imf1、imf2、imf3、imf4分别为eemd在第1阶、第2阶、第3阶、第4阶的固有模态函数。从图中可以看出，imf1和imf2中，除了有少部分相互重叠，三类碰撞声信号具有良好的可分性，但在imf3、imf4中，使用此特征的重叠部分相对比较多，特别是imf4的平均振幅改变的散点图相互混淆，无法区分三类碰撞声信号。

2.2统计wilson振幅值(wamp)

wilson振幅是指碰撞声信号中相邻两个时间段之间的振幅差异超过了预定义阈值的数目。当碰撞声信号振幅的变化超过给定阈值时，wilson振幅的值就加1，否则不变，计算公式如下：

其中，n取值2048。

对于每一个imf分量，图3分别展示三类碰撞声信号的wilson振幅均值图的变化情况。从均值图中可以看出，三类碰撞声的wilson振幅的变化趋势是相似的，对于第1阶imf分量，完好粒的wilson振幅平均值大于虫蛀粒和霉变粒的wilson振幅平均值，但是由于三种碰撞声信号wilson振幅平均值曲线的趋势相似，在第3阶imf分量之后，完好粒、虫蛀粒和霉变粒碰撞声信号的wilson振幅均值就相互重叠，所以本实施例提取前三阶imf分量的wilson振幅作为判别特征向量。

2.3计算峰峰值(pk-pk)

峰峰值是指在信号的时间序列中最大值和最小值之间的差异，即最大值和最小值之间的范围。峰峰值描述了信号幅值总体的变化范围，计算公式如下：

pk-pk＝max(ci，k)-min(ci，k)(3)

对于每一阶imf，图4展示了三类碰撞声信号的峰峰值的均值图变化趋势。在前4阶imf分量之前，完好粒碰撞声信号的峰峰值的均值大于其他两种碰撞声信号的峰峰值的均值，三类信号的峰峰值曲线相互不交叉，说明峰峰值可以用来区分三类信号。

2.4计算样本熵

对于每一个具有n个采样点的时间序列{x(n)＝x(1)，x(2)，…，x(n)}，样本熵可以定义为如下：

①按照序号形成一组大小为m的向量序号，xm(1)，…xm(n-m+1)，其中xm(i)＝{x(i)，x(i+1)，…，x(i+m-1)}，1≤i≤n-m+1，这些向量表示从i点开始的选取m个连续的x值。

②d[xm(i)，xm(j)]是向量xm(i)与xm(j)之间最大差值的绝对值，即两者之间的最大距离：

d[xm(i)，xm(j)]＝maxk＝0，1，…，m-1(|x(i+k)-x(j+k)|)(4)

③给定一个xm(i)，统计xm(i)与xm(j)之间距离小于或等于r的j(1≤j≤n-m，j≠i)的数目，并记为bi。对于1≤i≤n-m，定义：

④定义b^(m)(r)为：

⑤增加数组的维数到m+1，计算xm+1(i)与xm+1(j)(1≤j≤n-m，j≠i)距离小于等于r的个数，记为ai。定义为：

⑥定义a^m(r)为：

这样，b^m(r)是两个序列在相似容限r下匹配m个点的概率，而a^m(r)是两个序列匹配m+1个点的概率。样本熵定义为：

当n为有限值时，可以用下式估计：

图5是对每一类碰撞声信号的样本熵进行求和所得总值的柱状图。从图中可以看到，在第1阶imf时，三种碰撞声信号的样本熵有较大的不同。随着imf阶数的增加，玉米虫蛀粒碰撞声信号的样本熵总数下降比较快，而玉米霉变粒碰撞声信号的样本熵总数下降比较慢，差别逐渐减小。这是因为信号经过eemd分解之后，前几阶imf包含了信号的主要信息。在本发明实施例中，取m＝2，r＝0.25std(std代表原始数据的标准差)，n＝2048。

3.提取玉米粒碰撞声信号频域特征

图6(a)-图6(c)为玉米完好粒、虫蛀粒和霉变粒三类碰撞声信号的频域图。从图中可以看出，就整体趋势而言，霉变粒有别于完好粒和虫蛀粒，频域最大幅值也有较大的差异，特别是霉变粒碰撞声信号与其余两类玉米碰撞声信号相比，其最大幅值偏低，较大幅值集中区域在较低频率处，这与玉米粒被霉菌感染有关。

3.1计算频率均值

对原始碰撞声信号做fft变换，提取fft变换幅度均值作为特征值，计算公式如下：

其中，m代表采样点个数，取值2048。x(k)表示经傅里叶变换之后的幅度。

三类碰撞声信号的频率均值对比，如图7所示。将三类碰撞声信号按照以上公式求取碰撞声信号的频率均值，可以看出：霉变粒碰撞声信号的频率均值低于完好粒和虫蛀粒，完好粒的频率均值最高，此特征可以用于识别受损颗粒。

3.2计算功率谱均方根(rmsps)

提取功率谱的均方根作为特征向量，以反映功率谱能量的离散程度，计算公式如下：

其中，m代表采样点个数，取值2048，p代表功率谱的幅度。

如图8所示，三类碰撞声信号的功率谱均方根幅值波动区域有差别。三条曲线互不混淆，说明此特征向量可以检测和识别受损的玉米颗粒(虫蛀粒，霉变粒)。

3.3计算频带方差(fbv)

频带方差表示信号在不同频带的能量方差，物理意义是各频带间的波动特性和其短时能量特性。

其中，m代表采样点个数，取值2048，x(k)表示经fft变换之后的幅度。

三类碰撞声信号的频带方差对比，如图9所示。由图中可以看出，霉变粒的频带方差低于完好粒和虫蛀粒的频带方差，完好粒碰撞声信号的频带方差最高。这是因为玉米内部被虫啃食之后，质量变轻，导致虫蛀粒碰撞声信号的频带方差低于完好粒碰撞声信号的频带方差。霉变粒虽然内部并没有被啃食，但因其表面受到霉菌感染，导致其内部结构发生变化，因而和完好粒碰撞声信号有区别。

4.提取玉米粒碰撞声信号希尔伯特域特征

4.1计算低频带能量和高频带能量

将原始碰撞声信号经过hilbert变换之后得到hh谱，在5～15khz频率范围内，三类信号有差别，将此频段分为两个部分，即定义5～10khz为低频带，定义10～15khz为高频带。根据帕斯维尔定理，信号x(t)的总能量为：

其中，xk表示信号第k个采样点的幅值。

假设10～15khz为高频带瞬时强度为ah，i，5～10khz为低频带瞬时强度为al，i，高频带和低频带总能量分别为enh和enl：

其中，n取值2048。

三类信号的高频带能量和低频带能量的平均值、标准差如表1所示：

表1三类信号的高频带能量和低频带能量的平均值、标准差

由上表分析可知，完好粒的高、低频带能量均大于虫蛀粒和霉变粒，这是由于虫蛀粒和霉变粒受到损害，内部结构发生变化。

4.2计算包络平均值(ave)

信号的包络通常含有丰富信息，hilbert变换是求取信号包络的主要方法之一。先对信号进行eemd分解，然后求取信号的包络，具体步骤如下：

(1)对前几个imf做hilbert变换定义为：

其中，h[]为hilbert变换；

(2)分别以c(i)和hilbert变换h(ci(t))为实部和虚部，由此可以将原始信号解析表示如下：

z(t)＝c(i)+jh[ci(t)](18)

(3)信号的包络b(t)就是上述解析信号的幅值函数：

(4)包络谱则可以由对包络函数b(t)进行傅里叶变换求出：

a(f)＝f(b(t))

如图10所示，是三类碰撞声信号的包络平均值对比图。可以看出，完好粒碰撞声信号包络平均值波动比较大，而霉变粒碰撞声信号包络平均值的波动比较小，三类碰撞声信号的包络平均值波动的区域不同，相互之间并不混淆。

5.pso-svm分类器构建

5.1svm参数选择

svm(supportvectormachine)又称支持向量机，是一类新型机器学习方法，具有很好的学习性能。线性svm的原理是找到一个分类超平面作为决策曲面，最大化正例和反例的隔离边缘。先考虑分两类的情况，我们假设测试集数据是n维向量形成d维特征空间x。对每一个向量xi，与之相关联的目标函数yi∈{-1，+1}。在非线性情况下，通常将两类用核方法投影到高维空间，使得原空间中的非线性问题转化为新空间中的线性问题，它的决策规则是基于功能函数sign[f(x)]，其中f(x)代表在超平面中的判定函数，被定义如下：

f(x)＝ω*φ(x)+b^*(20)

定义最优超平面中的权向量和偏差是一个最小化代价函数，它表示为如下两个标准：利润最大化和误差最小化。

代价功能最小化有如下约束：

yi(ωφ(xi)+b)≥1-ξi，i＝1，2，…，n，ξi≥0，i＝1，2，…，n

其中，ξi表示的是不可分数据的松弛变量，常数c代表误差项的惩罚因子。求解svm中的超平面问题可以通过建立拉格朗日函数将最初求解问题转换为求解如下最优问题：

有如下约束：

其中，α＝[α1，α2，…，αn]是拉格朗日函数的向量，k(xi，xj)是核函数，计算公式如下：

k(xi，xj)＝exp(-g||xi-xj||²)，g>0(24)

核函数类型确定之后，求解最优问题还需要确定两个未知数，即式(23)中的惩罚参数c及式(24)中核函数参数g，两者的选取对于分类器的性能很重要。本文采用pso来选择惩罚参数c以及核函数参数g，以提高支持向量分类器的分类效果。

5.2利用pso算法优化svm参数

所谓pso(particleswarmoptimization)，又称粒子群优化算法或鸟群觅食算法，是一种随机优化技术。设pso是在一个n维空间进行搜索目标，群落中共有m个粒子构成，定义x＝(x1，x2，…，xm)，那么如果将i个粒子用n维向量表示，则可以表示为xi＝(xi1，xi2，…xin)^t，第i个粒子飞行的速度表示为vi＝(vi1，vi2，…vin)^t，pi＝(pi1，pi2…，pin)t代表第i个粒子迄今为止搜索到的最佳位置，也称个体极值。整个粒子群搜索到的最佳位置为pb＝(pb1，pb2，…，pbn)^t，也称全局最优。在迭代过程中，粒子在搜索空间内单位迭代次数是由vi决定，粒子按照式(25)和式(26)更新自身的速度和位置，计算公式如下：

其中，rand1和rand2为属于[0，1]的随机分布函数，c1、c2为学习因子，相对于局部最优和全局最优调整相对速度。进一步地说，这些参数被认为是相对于局部最佳位置和全局最佳位置的比例因子，它们决定粒子的最佳位置。ω为惯性权值，ω∈[0.8，1.2]，实验中ω的调整公式为：

ω＝ωmax-iter*(ωmax-ωmin)/itermax(27)

其中，iter为当前迭代次数，itermax为最大迭代次数。

本实施例利用pso算法对svm参数进行优化，用两个粒子分别代表式(23)中的惩罚参数c及式(24)中核函数参数g来优化，如图11所示，具体步骤如下：

(1)初始化粒子群，随机产生粒子的位置和速度；

(2)利用碰撞声信号样本，训练svm分类器并计算适应度值，采用交叉验证的方法计算适应度值；

(3)根据适应度值更新个体极值和全局最优值；

(4)更新当前粒子的速度和位置；

(5)当达到迭代次数或满足适应度条件时，则终止迭代；否则返回步骤(2)。

本实施例利用pso算法优化svm参数，构建pso-svm分类器，目的在于对所提取的玉米粒碰撞声信号的时域特征、频域特征和希尔伯特域特征进行分类识别，以期获得更高的玉米完好粒、霉变粒和虫蛀粒识别率。

6.结果分析

6.1实验检测流程

本实施例展示了基于多域融合的玉米受损颗粒检测系统流程，如图12所示。对信号进行eemd分解，本文中加入白噪声的幅值系数为0.2，总体平均次数为100，计算得到每一阶imf。分别从时域、频域和希尔伯特域提取可区分信号的特征，随后将这些特征输入到pso-svm进行检测识别。本实施例共采用玉米碰撞声信号840个，其中包括280个完好粒、280霉变粒和280个虫蛀粒。从每类碰撞声信号中随机抽取140个，共420作为训练集，剩下的玉米碰撞声信号作为测试集，数据集的详细说明见表2。

表2训练集和测试集

6.2基于时域结果分析

通过前面的分析，时域特征提取前3个平均振幅变化、前3个威尔逊振幅、前4个峰峰值、前2个样本熵作为判别特征。为了全面分析所提特征的有效性，表3展示了特征子集多种组合的识别率，目的是为了找到最优的组合。从表中可以看出，采用单个特征时，霉变粒的识别率相对比较低，最高识别率达到85.7％。随着组合的特征数目增加，可看出霉变粒识别率在增加，即当使用全部的12个特征作为判别特征时，完好粒、虫蛀粒和霉变粒的识别率分别达到97.1％、95.7％、93.6％，由此可知，提取的时域特征能够检测出受损玉米颗粒。

表3时域不同特征子集下检测识别率

6.3基于频域结果分析

表4给出了频域特征子集下的检测识别率，可以看出，随着特征数目的增加，识别率不断提高，当使用频域全部特征时，完好粒、虫蛀粒和霉变粒的识别率分别可以达到98.6％、98.6％、97.8％。

表4频域不同特征子集下检测识别率

6.4基于希尔伯特域结果分析

表5描述了当使用希尔伯特特征时，各个特征子集的检测识别率。单独使用一类特征时，霉变粒的识别率比较低，即对于霉变粒的可分性低于完好粒和虫蛀粒。而随着特征数目的不断增加，三类碰撞声信号的可分性增强。

表5希尔伯特域不同特征子集下检测识别率

6.5基于多域结果分析

本实施例的分类方法使用pso-svm，采用pso算法来选择svm中的径向基核函数中的参数g和惩罚系数c，通过200次迭代，5折交叉验证，最优粒子的适应度值变化如图13所示。最佳适应度曲线表示的是最优粒子的适应度变化，平均适应度表示所有粒子适应度值的平均，根据最初设定的参数，最优粒子的适应度从初始化参数时呈上升趋势，直到迭代次数到达一定数时，适应度值才呈平稳不变的状态。获得的参数c＝9.15、g＝0.01，将训练后的pso-svm用于测试样本的判别，分类结果如表6。

表6多域不同特征子集下检测识别率

表6展示了不同域特征子集下的检测识别率，通过将三个域的特征作为pso-svm的输入，得到完好粒、虫蛀粒和霉变粒的识别率分别是98.6％、99.2％、99.3％。

从表中还可以看出，通过将时域、频域和希尔波特域的特征进行融合后送入到svm和pso-svm分类器中，通过比较两者之间的识别率，得出pso优化svm可以提高各类识别率。pso优化svm的分类模型比单独运用svm算法能更快地找到最佳参数，也避免了随机设定svm参数的盲目性，获得更高的分类识别率。