基于主元分析法的磨音信号特征提取.pdf

球磨机是物料粉碎的主要设备， 广泛应用于物 料加工和矿物选矿。 磨机工作时的磨音与其负荷具 有一定的关系[１]，因此常利用球磨机的磨音信号来 间接检测球磨机的负荷。 因球磨机磨音信号的频谱 信息存在非线性、高维性以及复杂性等特点，研究人 员常通过对频谱信息进行特征提取， 来提高检测精 度和稳定性。 姜德轩[２]使用小波阈值消噪和带通滤 波器对采集的磨音信号降噪， 以数字信号处理器 （ DSP） 芯片为基础设计了一套磨音检测判别装置， 能够实时对磨机的负荷状态进行检测。 方涛[３]建立 了球磨机料位的数学模型并进行验证， 然后采用小 波包分析对磨音信号进行处理， 并以此为基础应用 嵌入式技术设计球磨机料位检测系统。 王飞等[４]人 依据实时采集的磨音信号来提取出所需参数值，然 后基于 Stein 无偏风险估计 （ SURE）得到信号系数的 最优阈值，再进行去噪，得到的磨音信号能更加精确 地反映磨机负荷。 张景等[5]人对球磨机磨音和球磨 机负荷之间的关系进行分析， 将得到的磨音有效频 段的能量作为特征，采用径向基函数 （ RBF）神经网 络建立预测模型来对球磨机的负荷进行预测。 梁朝 霞[6]提出对采集的球磨机的振动和磨音信号的特征 频段分段，用分频段能量作为检测负荷的依据，基于 BP 神经网络建立负荷预测软测量模型对球磨机负 荷进行预测。 惠瑜[7]通过分析磨音强度与磨机负荷 之间的关系，采集球磨机磨音信号，并对磨音信号进 行频谱分析和压缩处理，再采用 RBF 神经网络建立 预测模型来预测球磨机的负荷。 本研究简述主元分 析法的实现过程，给出流程图，通过实验室试验采集 磨音信号， 采用 Welch 法对磨音信号进行频谱分 析，并利用主元分析法对功率谱进行特征提取，再通 过支持向量机预测模型对提取后的特征进行验证并 给出结论。 1功率谱估计和主元分析 1．1Wｅｌcｈ 功率谱估计 Welch 法实际上是对周期图法采用平均和平滑 的思想改进后的一种方法， 其谱估计是将原长度为 N 的音频信号数据分成 K 段，每段长 MN/K，并对 分段后数据进行加窗来克服矩形窗导致的分辨率差 的缺点[8]。 求解时，先对每段数据点加矩形窗，由此 可以得到平均周期图法，见式 （ １）。 （ １） 式中XiN（ n）为加窗后的数据段。 p -（ w） １ MK K i１ Σ M－１ n０ ΣXiNe－jwn ２ 第 35 卷第 3 期 ２０２０ 年 6 月 Ｖｏｌ．35，Ｎｏ．3 Jｕｎ．２０２０Ｃｈｉｎａ Ｔｕｎｇｓｔｅｎ Ｉｎｄｕｓｔｒｙ DOI１０.３969/j.issn.１００9－０6２２.２０２０.０３.０１１ 基于主元分析法的磨音信号特征提取 浦友尚，魏镜弢，吴张永，钱杰，冯先丁 （ 昆明理工大学 机电工程学院，云南 昆明 65０5０４） 摘要为有效提取球磨机磨音信号特征实现负荷识别，提出一种基于主元分析法 （ PCA）的特征提取方法。 以实验 室球磨机为对象，采集球磨机的工作磨音信号，经降噪处理后，通过 Welch 法对磨音信号进行功率谱估计，并通过主 元分析法对磨机磨音功率谱分段后得到的基本特征进行特征提取。 通过支持向量机预测模型对特征提取后的样本 集进行仿真验证。 结果表明，特征提取后的特征有更好的识别效果和更小的随机性，可作为负荷预测模型的输入，实 现球磨机负荷的预测。 关键词特征提取；磨音信号；主元分析法；球磨机；负荷预测 中图分类号TD9２１.４；TD４5３文献标识码A 收稿日期２０２０－０４－２5 资助项目国家自然科学基金 （ 5１１65０１２） 作者简介浦友尚 （ １996－），男，云南宣威人，硕士研究生，研究方向为复杂机电系统及集成。 通讯作者魏镜弢 （ １96４－），男，四川三台人，博士，教授，主要从事矿物加工理论及应用工作。 万方数据 第 ３ 期 对所求的结果归一化，得到 （ ２） 式中U 为归一化因子；d２（ n）为数据加窗函数。 1．２主元分析 （ PＣA） 主元分析法是特征提取和数据压缩中比较经典 的方法[9－１０]，它是基于多元统计的一种分析方法，其 基本思路是将高维样本数据信息投影到一个低维的 子空间中， 子空间中的低维信息能够表示样本数据 的主要信息。 它主要是将一个高维样本数据中的各 个随机分向量 x （ x１ ，x ２ ，x ３，xm） T通过正交变换变 换为不相关的随机分量 y （ y１ ，y ２ ，y ３，ym） T。 主元分 析法可以减小原始样本数据空间中分量间的相关 性，在这些分量中提取出独立的特征，减小原始样本 数据的冗余度和无效特征， 并且这些分量能够尽量 多地体现出原始样本数据中的信息[１１]。 假设原始数 据矩阵 Y 如公式 （ ３）所示，原始数据中有 m 个样本 数据，n 个观测量。 （ ３） 主元分析法的实现过程如下[１２] 对原始数据矩阵中的数据进行 Z－score 标准化 处理，得到标准化矩阵y *ij。 （ ４） 其中yij表示原始样本数据中的第 i 样本上的第 j 个观测量；y － j为矩阵 Y 中第 j 列均值，Sj为第 j 列的 标准差，计算公式分别如式 （ 5）、式 （ 6）所示。 （ 5） （ 6） 计算样本矩阵的相关系数矩阵。 （ 7） R 为矩阵 YS的相关系数矩阵，又称为协方差矩 阵。 在式 （ 7）中 rij为矩阵 YS中任意两个元素 （ yj和 yi） 之间的相关系数，计算公式如式 （ 8）所示。 （ 8） 求协方差矩阵 R 的特征根与特征向量。 （ 9） 计算主元矩阵 T。 TYmnLnn t１t２ t n []（ １０） 计算主元贡献率和累计贡献率。 δλiλi/ n i１ Σλi １００ （ １１） （ １２） 前 r 个主元的累计贡献率也就是前 r 个特征值 的累计方差贡献率。 主元分析法的流程图如图 １ 所 示。 （ １３） （ １４） Lnn α１ α２ αn 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 α１１α１２ α１n α２１α２２ α２n αn１αn２ αnn 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 rij m k１ Σy * ki－y －* i 移移y * kj－y －* j 移移 m k１ Σy * ki－y －* i 移移 ２ m k１ Σy * kj－y －* j 移移 ２ 姨 p（ w） １ MU K i１ Σ M－１ n０ ΣXiNd２（ n）e－jwn ２ r i１ Σ Sαi n i１ ΣSαi r i１ Σλi/S φr r i１ Σδλi r i１ Σλi/ n i１ Σλi １００ Ymn y１１y１２ y１n y２１y２２ y２n 移 移 移 移 移 移 移 移 移移 移 移 移 移 移 移 移 移 移 移移 移 Y１Y２ Yn [] Y（１） Y（２） Y（m） 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 移 y *ij yij－y － j 移移/Sj y － j １ m m i１ Σyij Sj １ m－１ m i１ Σyij－y － j 移移 ２ 姨 （ j１，２，，n） R r１１r１２ r１p r２１r２２ r２p rp１rp２ rpp 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 姨 λi/ n i１ Σλiλi/S 图 1主元分析算法的流程图 Fｉｇ．1Fｌｏw cｈａｒｔ ｏf ｔｈｅ pｒｉｎcｉpａｌ cｏmpｏｎｅｎｔ ａｎａｌｙｓｉｓ ａｌｇｏｒｉｔｈm 浦友尚，等基于主元分析法的磨音信号特征提取 原始数据矩阵 Y 标准化矩阵 YSYs 的协方差矩阵 R 求矩阵 R 的特征值 和特征向量 输出主元矩阵 选取主元的数目， 取对应 的特征向量组成主元矩阵 将特征向量按照特征 值的大小排序 69 万方数据 第 ３5 卷 ２磨音信号检测 ２．1磨音采集试验 采用型号为 QM－２SP１２ 的行星式球磨机作为 试验设备， 通过变频器调节球磨机筒体转速， 采用 TP6CT9４ kW 变频器调节转速， 传声器为 MPA４１6，数据采集仪是 NI 9２３４。采用云南某矿山的 锡矿石作为试验物料，粒度为 ３～7 mm。磨矿介质是 直径分别为 6 mm、１０ mm、２０ mm 的 ３０４ 不锈钢钢 球， 其密度为 7.8 t/m３。 试验中以 ３０ 的介质填充 率，配比为 6∶３∶１ 的钢球作为定量，给料量作为变量， 给定不同的物料重量来表示球磨机的负荷大小变 化。 相关的试验参数见表 １。 表 1每组试验的球磨机负荷参数 Ｔａb．1Bａｌｌ mｉｌｌ ｌｏａｄ pａｒａmｅｔｅｒｓ fｏｒ ｅａcｈ ｇｒｏｕp ｏf ｅxpｅｒｉmｅｎｔｓ ２．２磨音信号功率谱估计 磨矿介质和磨矿介质之间、 磨矿介质和物料之 间撞击产生的声音主要受到筒体内的负荷影响。 随 着球磨机内部负荷大小的变化， 采集到的磨音信号 的频率与幅值也各不相同。 试验中为避免磨机起磨、 停磨过程中磨音的干 扰，待磨机运行 ３０ s 后再进行数据采集。 每组试验 采集磨音信号数据时长为 １6０ s，以 ４ s 为单位时间， 将每个组的试验数据分为 ４０ 个样本，采样频率 Fs 5１ ２００ Hz，每个样本的数据长度为 NFsT5１ ２００ ４２０４ 8００。 为了确保较高的频率分辨率，将数据点 分为 ４ 段，每段数据长度为 MN/K5１ ２００，对数据 进行汉宁窗加窗处理， 每段数据功率谱估计累加求 均值以获得平均功率谱估计。 不同负荷下球磨机磨 音信号的功率谱如图 ２ 所示。 观察功率谱图后可知，功率谱幅值在 １5 ０００ Hz 后已经很小， 磨音信号数据有效频带确定为 ０～ １5 ０００ Hz，为了在特征提取时能提高计算效率和准 确性，需要对功率谱数据进行处理，将频段分为 １5０ 段，得到了每个样本的特征维度为 １5０ 维。为更好地 说明球磨机负荷的变化与磨机磨音信号之间关系， 选取球磨机在负荷较小、适中、较大 ３ 种情况时的功 率谱图，如图 ３ 所示。 观察图 ３ 可知，随着球磨机负荷增大，磨音信号 的功率谱值减小，在 b 点尤为明显。 从 a 点可知，负 荷增大， 有效频率段向低频靠近且功率谱值也随之 减小。 而在 d、e 两点，当磨机负荷从较小到适中时， 功率谱值并无明显变化； 当磨机负荷从适中到较大 时，功率谱值则明显减小。 由分析可知，球磨机磨音 信号功率谱变化存在随机性、 复杂性和非线性等特 点，且在有效频段中存在部分无效、冗余的特征。 试验 编号 介质填 充率/ 各级钢球重量/g 物料重 量/g 球磨机 负荷/ 6 mm１０ mm２０ mm １２5０１０ ２5００２０ ３75０３０ ４１ ０００４０ 5１ ２5０5０ 6１ 5００6０ 7１ 75０7０ 8２ ０００8０ 9２ ２5０9０ １０２ 5００１００ ３０ ３０ ３０ ３０ ３０ ３０ ３０ ３０ ３０ ３０ ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 ２ 9４8 １ ４7４ １ ４7４ １ ４7４ １ ４7４ １ ４7４ １ ４7４ １ ４7４ １ ４7４ １ ４7４ １ ４7４ ４9２ ４9２ ４9２ ４9２ ４9２ ４9２ ４9２ ４9２ ４9２ ４9２ 图 ２不同负荷时的功率谱图 Fｉｇ．２Pｏwｅｒ ｓpｅcｔｒａ ａｔ ｄｉffｅｒｅｎｔ ｌｏａｄｓ ０ 功率谱１０－３/dB ０ ２ ４ 6 ３ ０００ 6 ０００ 9 ０００ １２ ０００ 频率/Hz 负荷/ １０ ３０ 5０ 7０ 9０ １００ ２０ ４０ 6０ 8０ 7０ 万方数据 第 ３ 期 图 3负荷小、适中、大时磨音信号的功率谱图 Fｉｇ．3Pｏwｅｒ ｓpｅcｔｒａ ｏf ｓmａｌｌ，mｏｄｅｒａｔｅ ａｎｄ ｌａｒｇｅ ｔｉmｅ ｇｒｉｎｄｉｎｇ ｓｉｇｎａｌｓ 3磨音信号特征提取及验证 根据图 １ 给出的主元分析算法的流程图对分段 后的基础特征数据进行特征提取。将 １０ 组试验对应 的基础特征数据作为 PCA 算法的输入，组成输入矩 阵 Y４００１5０，求出矩阵对应的特征值和特征向量并进行 大小排序， 计算特征值对应的主元贡献率和累计贡献 率，如表 ２ 所示，主元累计贡献率的变化如图４ 所示。 表 ２各主元的特征值、贡献率和累计贡献率 Ｔａb．２Ｃｈａｒａcｔｅｒｉｓｔｉc vａｌｕｅ，cｏｎｔｒｉbｕｔｉｏｎ ｒａｔｅ ａｎｄ cｕmｕｌａｔｉvｅ cｏｎｔｒｉbｕｔｉｏｎ ｒａｔｅ ｏf ｅａcｈ pｒｉｎcｉpａｌ 图 4主元累计贡献率变化图 Fｉｇ．4Pｒｉｎcｉpａｌ cｕmｕｌａｔｉvｅ cｏｎｔｒｉbｕｔｉｏｎ ｒａｔｅ cｈａｎｇｅ 本研究采用累计贡献率方法来确定主元个数， 即要求累计贡献率达到 9０ 或者以上，为保证特征 提取后的数据能保留原始数据大部分的特征， 将累 计贡献率 96 作为确定主元个数的指标。 分析表 ２ 和图 ４ 后，将主元个数确定 ２１，将原始数据的基本 特征维数降低到 ２１ 维。其中，由图 ４ 可知，选取主元 中的第一主元、 第二主元和第三主元的累计贡献率 分别达到了 ４9.9３ 、6４.9２ 和 7２.77 。 分别把第 一主元和第二主元投影在二维平面中，将第一主元、 第二主元和第三主元投影到三维立体图中， 以便更 直观地观察经过特征提取后特征的整体识别效果， 如图 5 和图 6 所示。 分析图 5 和图 6 可得， 对球磨机不同负荷工作 时采集的磨音信号进行频率谱估计后， 再基于主元 分析法特征提取后的特征， 除了部分主元因工作负 荷接近导致特征区分度较小而出现混杂外，其余大 图 5第一、第二主元投影图 Fｉｇ．5Pｒｏjｅcｔｉｏｎ ｏf ｔｈｅ fｉｒｓｔ pｒｉｎcｉpａｌ cｏmpｏｎｅｎｔ ａｎｄ ｔｈｅ ｓｅcｏｎｄ pｒｉｎcｉpａｌ cｏmpｏｎｅｎｔ 图 6第一、第二、第三主元三维投影图 Fｉｇ．6Pｒｏjｅcｔｉｏｎｓ ｏf ｔｈｅ fｉｒｓｔ, ｓｅcｏｎｄ, ａｎｄ ｔｈｉｒｄ pｒｉｎcｉpａｌ cｏmpｏｎｅｎｔｓ －１０ －5 ０ 5 －１００１０２０ 第一主元 第二主元 １０ ２０ ３０ ４０ 5０ 6０ 7０ 8０ 9０ １００ －１０ ０ １０ １０ ０ －１０ －２０ －１０ ０ １０ ２０ 第一主元 第二主元 第三主元 １０ ２０ ３０ ４０ 5０ 6０ 7０ 8０ 9０ １００ 主元号特征值贡献率累计贡献率 １6３.9０9 9０.４99 ３０.４99 ３ ２１9.１9２ 9０.１４9 9０.6４9 ２ ３１０.０４6 ００.０78 5０.7２7 7 ４7.999 7０.０6２ 5０.79０ ２ 55.97４ 9０.０４6 7０.8３6 9 6３.686 5０.０２8 8０.865 7 ２１０.１65 6０.００１ ３０.96０ 8 ２２０.１6０ ００.００１ ３０.96２ １ １5００.０００ ００.０００ ０１ ０.４ ０.5 ０.6 ０.7 ０.8 ０.9 １ ０２０４０6０8０１００ 累计贡献率 主元个数/个 浦友尚，等基于主元分析法的磨音信号特征提取 ０ ０.０１ ０.０３ ０.０４ ０.０6 ０6００ 频率/Hz 负荷小/１０ 负荷适中/5０ 负荷大/１００ a 功率谱/dB b c d e ０.０２ ０.０5 5００４００３００２００１００ 7１ 万方数据 第 ３5 卷 部分主元都有很好的分类识别效果， 能很好地区分 不同负荷下球磨机的工作状态， 为预测模型的输入 数据做进一步的优化。 将通过 PCA 特征提取后的样本集数据进行归 一化处理， 使用支持向量机工具箱来建立 PCA－SVM 预测模型。用训练集数据训练预测模型，而将测试集 数据作为负荷预测模型的输入来进行预测。 然后建 立相同的预测模型， 通过未经过特征提取的原始的 训练集数据进行训练， 测试集数据作为负荷预测模 型的输入来进行预测。 引入均方根误差 （ RMSE）、 平均绝对误差 （ MAE）、平均百分比误差 （ MAPE）、最大误差 ４ 个评 价指标和运算时间， 以便更好地比较两个模型的优 劣性，２ 种预测模型的结果评价指标见表 ３。 表 3２ 种预测模型的结果评价指标 Ｔａb．3Ｔｈｅ ｒｅｓｕｌｔ ｅvａｌｕａｔｉｏｎ ｉｎｄｅxｅｓ ｏf ｔｈｅ ｔwｏ pｒｅｄｉcｔｉｏｎ mｏｄｅｌｓ 由表 ３ 可知，PCA－SVM 预测模型的 ４ 个评价 指标都高于 SVM 预测模型， 这说明了通过 PCA 特 征提取能降低特征之间的冗余度以及减少了干扰信 息和误差，提高了球磨机负荷预测模型的精确度，同 时缩短了负荷预测模型的计算时间。 4结语 研究探究了采用主元分析法 （ PCA）对球磨机磨 音信号进行特征提取， 将特征提取后的特征作为支 持向量机预测模型的输入， 对磨机工作负荷进行预 测。通过评估指标可知，经过特征提取后的特征相比 原始的特征分布更加确定、随机性更小，大大降低了 特征的维数和特征之间的冗余度， 球磨机负荷预测 模型的精度得到有效提高， 同时用特征提取后的特 征作为预测模型的输入能缩短预测模型的计算时 间。本研究是在实验室条件下进行磨音信号的采集， 该方法还需结合工业试验进行深入研究。 参考文献 [１]SPENCER S J，CAMPBELL J J，WELLER K R.Acoustic emissions monitoring of SAG mill perance [J].Intelligent Processing and Manufacturing of Materials，１999，２ （ ４）9３6－9４6. [２]姜德轩. 基于球磨机负荷优化控制的磨音检测[D].济南济南大 学，２０１１. JIANG Dexuan.The mill noise detection based on optimal control of mill load[D]. Jinan University of Jinan，２０１１. [３]方涛.基于音频信号的球磨机料位检测研究[D].桂林桂林电 子科技大学，２０１２. FANG Tao.The research on the detection of ball mill level based on audio signal[D].GuilinGuilin University of Electronic Technology， ２０１２. [４]王飞，李智勇，朱强.基于自适应阈值小波分析的磨音信号 去噪[J].矿山机械，２０１5 （ １２）76－79. WANG Fei，LI Zhiyong，ZHU Qiang.Processing of mill noise based on adaptive threshold wavelet analysis [J].Mining grinding signal; principal component analysis; ball mill; load prediction （ 编辑游航英） 浦友尚，等基于主元分析法的磨音信号特征提取7３ 万方数据