CN111898508B - 一种基于听觉感知的电动冲击批零件缺损的检测方法 - Google Patents

Abstract

本发明涉及电动工具的质量检测技术领域，具体涉及一种基于听觉感知的电动冲击批零件缺损的检测方法。传统的电动冲击批零件缺损检测方法依靠听测工人完成，难以用稳定且一致的标准来评判电动冲击批是否存在零件缺损情况。本发明基于听觉感知模拟人工听音来分析电动冲击批打滑状态声音，对稳定的时变响度数据作离散傅里叶变换DFT处理得到响度变化周期，消除了人工听测的主观感受影响。本发明提出了时变特征响度频带能量信息熵下的权重系数，该系数对频带特征信息进行了增强，便于故障信息的特征提取。本发明中采用声音信号权重下的时变响度的峰均比统一衡量和评判，工程易实现。

Description

一种基于听觉感知的电动冲击批零件缺损的检测方法

技术领域

本发明涉及电动工具的质量检测技术领域，具体涉及一种基于听觉感知的电动冲击批零件缺损的检测方法。

背景技术

电动冲击批也称为电动冲击扳手，可将其内部冲击机构的冲击力转化为作用于螺丝的扭力。由于它以每分钟数百次脉冲对螺丝钉加以高强度的扭力，使用方便，可连续操作，被广泛应用于机械工业和建筑行业中。

电动冲击批在出厂前必须对其进行质量检测，以确保其内部零件没有缺损。目前生产线中常用的检测方法是依靠工人听测电动冲击批打滑时的声音，并对其打滑时发出的声音进行人耳听音来评判是否存在零件缺损。由于工人受技术成熟程度、身体与情绪状况和工作疲劳等因素影响，听测的主观感受差别较大，难以用一个稳定且一致的标准来评判电动冲击批是否存在零件缺损情况，导致企业生产效率和产品质量难以进一步提高。

振动信号分析也被用来对该类旋转机械进行故障诊断(俞昆，罗志涛，李鸿飞，马辉.广义参数化同步压缩变换及其在旋转机械振动信号中的应用[J].机械工程学报，2019，55(11)：149-159)；但是被测工具的壳体会屏蔽或者改变一部分机械振动信号特性(周瑞丽，姚春燕，杨涛，郑士先.手持冲击类工具壳体动力学特性分析与结构优化设计[J].机械工程师，2016，(7)：93-95)，同时在工程应用中需要将振动传感器安装在被测工具表面，不利于自动化检测系统的集成。

声信号分析也是旋转机械故障诊断的主要手段(武倩平.基于噪声的旋转机械故障诊断研究[D].徐州：中国矿业大学硕士学位论文，2017)，但相关特征多数是利用现代信号处理手段从信号的时域、频域、时频域等进行提取，未能从人的听觉感知进行判断，导致分析方法与从振动等手段无明显差异，未能体现出声信号分析的优势。

发明内容

为了解决上述问题，本发明提供了一种基于听觉感知的电动冲击批零件缺损的检测方法，基于人耳的非线性作用提出了时变特征响度频带能量信息熵权重系数，利用信息熵权重下的时变响度模拟了人的听觉感知对电动冲击批工作时发出的声音进行分析和评价，以此代替人工听测，一方面可降低工人的工作强度，保护听力健康，另一方面也可保证检测结果的客观性和一致性，提高生产企业的综合竞争力。本发明的具体技术方案如下：

一种基于听觉感知的电动冲击批零件缺损的检测方法，包括以下步骤：

S1：使用传声器采集电动冲击批打滑状态的声音信号，设传声器的采样频率为f_s，单位为Hz；声音信号时长为T，单位为s；按照德国DIN45631/A1标准对声音信号计算各个频带的时变特征响度v(n,z_i)，其中n是指将信号时长T按照2ms为分度划分的点数；

i是指将24个临界频带Bark按照0.1Bark为分度划分为240个子频带，i为对应的频带号，i＝1,2,…,240；z_i为心理声学的24个临界频带Bark的频率尺度；

S2：计算时变特征响度v(n,z_i)各频带的能量占整个频带的总能量的比值η(z_i)以及η(z_i)中最大值η(z_i)_max，得到η(z_i)_max对应的频带号i：

η(z_i)_max＝max[η(z_i)]； (4)

其中N(z_i)是指第i个子频带的能量，N_total是指整个频带的总能量；

S3：对各频带的时变特征响度v(n,z_i)作离散傅里叶变换DFT处理后得r(k_z,z_i)，求取经离散傅里叶变换DFT处理后的时变特征响度r(k_z,z_i)各频带的能量信息熵E(z_i)：

其中，k_z是指频率信息；

为向下取整运算符；

S4：对时变特征响度r(k_z,z_i)各频带的能量信息熵E(z_i)进行多项式拟合，得到信息熵权重系数g(z_i)；

g(z_i)＝a₀+a₁*E(z_i)+a₂*E(z_i)^2+a₃*E(z_i)^3+a₄*E(z_i)^4+a₅*E(z_i)^5； (7)

其中a＝[a₀,a₁,a₂,a₃,a₄,a₅]为多项式拟合系数；

S5：利用人耳听觉的非线性作用，求取信息熵权重下的时变特征响度s(n,z_i)，即各频带时变特征响度乘以对应频带的信息熵权重系数：

s(n,z_i)＝v(n,z_i)·g(z_i)； (8)

S6：对s(n,z_i)进行运算，得到信息熵权重下的时变响度x(n)；

S7：对x(n)作离散傅里叶变换DFT得到y(f)；并计算|y(f)|中[c,d]频带内的峰均比e，具体方法为：

查找[c,d]频带内|y(f)|的幅值最大值A₁和对应频率f₁，分别计算下临界频点f₂和上临界频点f₃为：

|y(f)|是幅值谱，表征信号的幅值随频率f的分布情况；Δf是扩展常数，表征最大值所在波峰的统计半宽；

计算频带[c,f₂]和[f₃,d]内|y(f)|的平均值A₂，则峰均比e为：

S8：将各频带的能量占整个频带的总能量的比值中的最大值η(z_i)_max、对应频带号i、峰均比e和频率f₁构造输入特征向量m＝[η(z_i)_max,i,e,f₁]^T；

S9：重复步骤S1-S8，得到若干个特征向量，划分为训练样本集和测试样本集；

S10：利用逻辑回归分类模型，将训练样本集中的各个特征向量m作为逻辑回归分类模型的输入，对逻辑回归分类模型进行训练；

S11：将测试样本集中的各个特征向量m送到已经训练好的逻辑回归分类模型中进行分类，当得到的预测概率小于等于0.5时，即认为该样本为零件缺损样本；当得到的预测概率是大于0.5，即认为该样本为合格样本。

优选地，所述步骤S7中Δf的选取方法为：

(1)查找测试样本集中每个合格产品的声音信号的整个幅值谱|y(f)|的峰值对应的频率点f_0i，并找到峰值左右两边的第一个峰谷的对应频率，设为f_0il和f_0ir；i表示第i个合格产品；

(2)计算第i个合格产品的声音信号的波峰宽度δ_i＝f_0ir-f_0il；

其中，δ_i表示第i个合格产品的声音信号的整个波峰的宽度；

(3)计算测试样本集中采集的每个合格产品的声音信号的波峰宽度，并对计算得到的波峰宽度作统计平均，具体如下：

其中，samples表示测试样本集中的合格产品数；

(4)则

其中

表示向上取整。

优选地，所述频带[c,d]的确定方法为：

(1)对测试样本集中合格产品的声音信号的峰值频率f_0i做统计平均，得到测试样本集中合格产品的声音信号的峰值平均频率f_{0_}average，具体如下：

(2)以测试样本集中合格产品的声音信号的峰值平均频率f_{0_}average为中心，分别向左和向右各扩展(2～5)Δf即可得到频带[c,d]。

优选地，所述步骤S9中的逻辑回归分类模型为：

其中θ＝(θ₀,θ₁,θ₂,…,θ_t)为逻辑回归分类模型中的参数向量，其中y＝{0，1}代表样本类别，其中1代表合格类产品，0代表零件缺损类产品，x为特征向量，在本发明中x＝m；t为特征向量的个数，T表示向量转置符号，P(y|x；θ)代表样本的预测概率值，P(y＝1|x；θ)代表合格类产品的概率值，P(y＝0|x；θ)代表零件缺损类产品的概率值。

优选地，采用对数似然函数作为逻辑回归分类模型的损失函数J(θ)，具体如下：

单个样本的损失函数：

cost(h_θ(x),y)＝-y_jlog(h_θ(x))-(1-y_j)log(1-h_θ(x))； (16)

其中h_θ(x)是指P(y|x；θ)，j是指样本序号，y_j是指第j个样本；如果单个样本预测正确，单个样本的损失函数值就趋近于0，相反，单个样本预测不正确，单个样本的损失函数值就趋近于无穷大；

全体样本的损失函数：

其中j是指样本序号，y_j是指第j个样本，x_j是第j个样本的特征向量；如果全体样本预测正确，全体样本的损失函数值就趋近于0，相反，全体样本预测不正确，全体样本的损失函数值就趋近于无穷大；

当全体样本的损失函数的值无限趋近0的时候，则说明模型最优，即所有样本都预测准确，则逻辑回归分类模型训练完毕。

本发明采用对数似然函数J(θ)作为损失函数来计算代价值，并采用MATLAB的子函数fminun函数对上述全体样本的损失函数代价值进行优化，使得全体样本的损失函数值达到最小。

本发明的有益效果为：

传统的电动冲击批零件缺损检测方法依靠听测工人完成，难以用稳定且一致的标准来评判电动冲击批是否存在零件缺损情况。为了代替人工检测，本发明基于听觉感知模拟人工听音来分析电动冲击批打滑状态声音，对稳定的时变响度数据作离散傅里叶变换DFT处理得到响度变化周期，消除了人工听测的主观感受影响。与现有方法相比，本发明具有如下突出优点：

(1)传统检测方法依靠人工听测进行判断，而本发明基于听觉感知模拟人耳听音过程，消除了主观感受对检测结果的影响；

(2)在该发明中我们提出了时变特征响度频带能量信息熵下的权重系数，该系数对频带特征信息进行了增强，便于故障信息的特征提取。

(3)传统检测方法难以用稳定且一致的标准来评判电动冲击批是否存在零件缺损情况，本发明中采用声音信号权重下的时变响度的峰均比和各个频带的能量占整个频带总能量的比值来统一衡量和评判，工程易实现。

附图说明

图1为合格产品打滑状态下声音信号时域图。

图2为零件缺损产品打滑状态下声音信号时域图。

图3为合格产品打滑状态下声音信号频谱图。

图4为零件缺损产品打滑状态下声音信号频谱图。

图5为合格产品打滑状态下声音信号时频图。

图6为零件缺损产品打滑状态下声音信号时频图。

图7为合格产品打滑状态下声音信号时变特征响度图。

图8为零件缺损产品打滑状态下声音信号时变特征响度图。

图9为合格产品打滑状态下声音信号时变特征响度各频带能量占整个频带的总能量的比值变化图。

图10为零件缺损产品打滑状态下声音信号时变特征响度各频带能量占整个频带的总能量的比值变化图。

图11合格产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后的子带能量图。

图12零件缺损产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后的子带能量图。

图13合格产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后各子带能量图的信息熵。

图14零件缺损产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后的子带能量图的信息熵。

图15合格产品打滑状态下声音信号时变特征响度频带能量信息熵权重图。

图16零件缺损产品打滑状态下声音信号时变特征响度频带能量信息熵权重图。

图17信息熵权重下的合格产品打滑状态下声音信号时变特征响度图。

图18信息熵权重下的零件缺损产品打滑状态下声音信号时变特征响度图。

图19合格产品打滑状态下声音信号时变权重响度图。

图20零件缺损产品打滑状态下声音信号时变权重响度图。

图21为截取得到的合格产品时变权重响度数据离散傅里叶变换DFT频谱图。

图22为截取得到的零件缺损产品时变权重响度数据离散傅里叶变换DFT频谱图。

图23为本实施例中迭代100次的全体样本的损失函数的曲线图；

图24为本发明的检测实施流程图。

具体实施方式

为了更好的理解本发明，下面结合附图和具体实施例对本发明作进一步说明：

如图24所示，一种基于听觉感知的电动冲击批零件缺损的检测方法，包括以下步骤：

S1：使用传声器采集电动冲击批打滑状态的声音信号，设传声器的采样频率为f_s，单位为Hz；声音信号时长为T，单位为s；按照德国DIN45631/A1标准对声音信号计算各个频带的时变特征响度v(n,z_i)，其中n是指将信号时长T按照2ms为分度划分的点数；

i是指将24个临界频带Bark按照0.1Bark为分度划分为240个子频带，i为对应的频带号，i＝1,2,…,240；z_i为心理声学的24个临界频带Bark的频率尺度；本发明实施例中采用的电动冲击批是良明牌18V锂电池带冲击式手电钻。

在本发明实施例中，所采用的传声器是自由场测量传声器套装CRY333，其频率响应为3.15Hz～20KHz，灵敏度是50mV/Pa，采样率f_s为51200Hz，信号总时长为1.6s，打滑状态的信号时长T为1.5s。在本次发明中采用NI公司USB 6343型号采集卡和LabVIEW软件作为采集系统，将传声器信号输出端连接在采集卡的输入端，并将传声器放置于距离电动冲击批0.2m处的麦克风支架的位置，其中麦克风支架起固定作用，任意形式的三角形支架都可以。传声器采集电动冲击批打滑状态下的声音信号，采集得到信号时域图如图1和图2，其中图1为合格产品打滑声音信号，图2为零件缺损产品打滑声音信号。

对图1中的时域信号作离散傅里叶变换DFT和STFT时频分析，得到结果分别如图3和图5；对图2信号作相同分析的结果分别如图4和图6。图3为合格产品打滑状态下声音信号频谱图。图4为零件缺损产品打滑状态下声音信号频谱图。图5为合格产品打滑状态下声音信号时频图。图6为零件缺损产品打滑状态下声音信号时频图。对比图3和图4，可以发现从合格产品和零件缺损产品的频谱图中分辨不出产品有无缺损。对比图5和图6时频分析结果，可以发现频率随时间变化的特征二者区别不大。所以单从信号的时域、频域和时频域分析进行零件缺损产品检测的难度较大。

对图1和图2中的时域信号按照德国DIN 45631/A1标准计算时变特征响度，得到时变特征响度v(n,z_i)随时间变化的结果如图7和图8所示。图7为合格产品打滑状态下声音信号时变特征响度图。图8为零件缺损产品打滑状态下声音信号时变特征响度图。

S2：计算时变特征响度v(n,z_i)各频带的能量占整个频带的总能量的比值η(z_i)以及η(z_i)中最大值η(z_i)_max，得到η(z_i)_max对应的频带号i：

η(z_i)_max＝max[η(z_i)]； (4)

其中N(z_i)是指第i个子频带的能量，N_total是指整个频带的总能量；

对图7和图8中的时变特征响度v(n,z_i)计算各频带的能量占整个频带的总能量的比值，得到各频带能量比值变化结果分别如图9和图10所示。图9为合格产品打滑状态下声音信号时变特征响度各频带能量占整个频带的总能量的比值变化图。图10为零件缺损产品打滑状态下声音信号时变特征响度各频带能量占整个频带的总能量的比值变化图。

S3：对各频带的时变特征响度v(n,z_i)作离散傅里叶变换DFT处理后得r(k_z,z_i)，求取经离散傅里叶变换DFT处理后的时变特征响度r(k_z,z_i)各频带的能量信息熵E(z_i)：

其中，k_z是指频率信息；

为向下取整运算符；

S4：对时变特征响度r(k_z,z_i)各频带的能量信息熵E(z_i)进行多项式拟合，得到信息熵权重系数g(z_i)；

g(z_i)＝a₀+a₁*E(z_i)+a₂*E(z_i)^2+a₃*E(z_i)^3+a₄*E(z_i)^4+a₅*E(z_i)^5； (7)

其中a＝[a₀,a₁,a₂,a₃,a₄,a₅]为多项式拟合系数；多项式拟合系数利用MATLAB中的函数polyfit来进行确定即可。

图11和图12是各频带时变特征响度作离散傅里叶变换DFT后的结果，图11合格产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后的子带能量图；图12零件缺损产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后的子带能量图。对图11和图12计算能量信息熵，结果如图13和图14所示；图13合格产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后各子带能量图的信息熵；图14零件缺损产品打滑状态下声音信号时变特征响度作离散傅里叶变换DFT后的子带能量图的信息熵。在对图13和图14作多项式拟合得到信息熵权重系数曲线，结果如图15和图16；图15合格产品打滑状态下声音信号时变特征响度频带能量信息熵权重图；图16零件缺损产品打滑状态下声音信号时变特征响度频带能量信息熵权重图。

S5：利用人耳听觉的非线性作用，求取信息熵权重下的时变特征响度s(n,z_i)，即各频带时变特征响度乘以对应频带的信息熵权重系数：

s(n,z_i)＝v(n,z_i)·g(z_i)； (8)

S6：对s(n,z_i)进行运算，得到信息熵权重下的时变响度x(n)；

由于人耳的具有的非线性频谱分辨能力，对时变特征响度中的频带能量值进行信息熵拟合，得到频带能量的权重系数，如图15和图16所示。为了更好的体现其非线性特性，将时变特征响度与权重系数相乘得到权重下的时变特征响度，如图17和图18所示，图17信息熵权重下的合格产品打滑状态下声音信号时变特征响度图；图18信息熵权重下的零件缺损产品打滑状态下声音信号时变特征响度图。

将得到的权重下的时变特征响度进行积分运算得到权重下时变响度，如图19和图20所示，图19合格产品打滑状态下声音信号时变权重响度图；图20零件缺损产品打滑状态下声音信号时变权重响度图。从图19和图20中可以看出，合格产品和零件缺损产品的最大响度值相近，但响度随时间变化的趋势不同，合格产品响度随时间变化幅度较小，零件缺损产品的变化幅度较大。

S7：对x(n)作离散傅里叶变换DFT处理，得到y(f)；并计算|y(f)|中[c,d]频带内的峰均比e，具体方法为：

查找[c,d]频带内|y(f)|的幅值最大值A₁和对应频率f₁，分别计算下临界频点f₂和上临界频点f₃为：

|y(f)|是幅值谱，表征信号的幅值随频率f的分布情况，y(f)加绝对值是因为FFT之后得到的结果是复数；Δf是扩展常数，表征最大值所在波峰的统计半宽；

计算频带[c,f₂]和[f₃,d]内|y(f)|的平均值A₂，则峰均比e为：

其中，Δf的选取方法为：

(1)查找测试样本集中每个合格产品的声音信号的整个幅值谱|y(f)|的峰值对应的频率点f_0i，并找到峰值左右两边的第一个峰谷的对应频率，设为f_0il和f_0ir；i表示第i个合格产品；

(2)计算第i个合格产品的声音信号的波峰宽度δ_i＝f_0ir-f_0il；

其中，δ_i表示第i个合格产品的声音信号的整个波峰的宽度；

(3)计算测试样本集中采集的每个合格产品的声音信号的波峰宽度，并对计算得到的波峰宽度作统计平均，具体如下：

其中，samples表示测试样本集中的合格产品数；

(4)则

其中

表示向上取整。

频带[c,d]的确定方法为：

(1)对测试样本集中合格产品的声音信号的峰值频率f_0i做统计平均，得到测试样本集中合格产品的声音信号的峰值平均频率f_{0_}average，具体如下：

(2)以测试样本集中合格产品的声音信号的峰值平均频率f_{0_}average为中心，分别向左和向右各扩展(2～5)Δf即可得到频带[c,d]。

表1测试样本集中合格样品的δ_i

本实施例的Δf由表1计算得到，计算结果具体如下：

δ_average＝4.6129Hz；

在本实施例中f_{0_}average和[c,d]具体为：

f_{0_}average＝172.9692≈173Hz，分别向左和向右各扩展4Δf即可得到频带[c,d]为[133Hz，213Hz]。

对图19和图20数据作离散傅里叶变换DFT结果如图21和图22所示，图21为截取得到的合格产品时变权重响度数据离散傅里叶变换DFT频谱图；图22为截取得到的零件缺损产品时变权重响度数据离散傅里叶变换DFT频谱图。在变化频率小于50Hz的区间内，合格产品和零件缺损产品变化相似，幅值存在较为明显的差异，但根据幅值进行判断，方法鲁棒性较差。而在较高频段内，合格产品明显存在一个峰值，且该频段内零件缺损产品不存在峰值。

本实例对图21和图22中频段[133Hz，213Hz]内时变权重响度频谱数据进行分析。可发现与零件缺损产品相比，合格产品在188Hz频点明显具有峰值。取Δf为10Hz。合格产品频段[133Hz，213Hz]内最大值A₁为110.75，对应f₁为188Hz，f₂为178Hz，f₃为198Hz，平均值A₂为8.403，峰均比e为13.18。零件缺损产品频段[133Hz，213Hz]内最大值A₁为31.25，对应f₁为155Hz，f₂为145Hz，f₃为165Hz，平均值A₂为9.615，峰均比e为3.52。

S8：将各频带的能量占整个频带的总能量的比值中的最大值η(z_i)_max、对应频带号i、峰均比e和频率f₁构造输入特征向量m＝[η(z_i)_max,i,e,f₁]^T，本实施例中图21中的合格产品的特征向量为[0.89,153,13.18,188]^T，本实施例中图22中零件缺损产品的特征向量为[0.81,182,3.52,155]^T；

S9：重复步骤S1-S8，得到若干个特征向量，划分为训练样本集和测试样本集；

S10：利用逻辑回归分类模型，将训练样本集中的各个特征向量m作为逻辑回归分类模型的输入，对逻辑回归分类模型进行训练；逻辑回归分类模型为：

其中θ＝(θ₀,θ₁,θ₂,…,θ_t)为逻辑回归分类模型中的参数向量，y＝{0，1}代表样本类别，其中1代表合格类产品，0代表零件缺损类产品，x为特征向量，在该发明中x＝m；t为特征向量的个数，T表示向量转置符号，P(y|x；θ)代表样本的预测概率值，P(y＝1|x；θ)代表合格类产品的概率值，P(y＝0|x；θ)代表零件缺损类产品的概率值。

采用对数似然函数作为逻辑回归分类模型的损失函数J(θ)，具体如下：

单个样本的损失函数：

cost(h_θ(x),y)＝-y_jlog(h_θ(x))-(1-y_j)log(1-h_θ(x)) ；(16)

其中h_θ(x)是指P(y|x；θ)，j是指样本序号，y_j是指第j个样本，取值0或1。

全体样本的损失函数：

其中j是指样本序号，y_j是指第j个样本，xj是第j个样本的特征向量；

当全体样本的损失函数的值无限趋近0的时候，则说明模型最优，即所有样本都预测准确，则逻辑回归分类模型训练完毕。如图23所示，本实施例中迭代100次即可满足。

S11：将测试样本集中的各个特征向量m送到已经训练好的逻辑回归分类模型中进行分类，当得到的预测概率小于等于0.5时，即认为该样本为零件缺损样本；当得到的预测概率是大于0.5，即认为该样本为合格样本。最终完成基于听觉感知的电动冲击批零件缺损的目标分类。

本发明所示的这种检测方法，完整的计算过程如图24所示。首先对采集到的声音按照德国DIN 45631/A1标准计算时变特征响度，得到特征响度随时间变化的数据，计算时变特征响度中各频带能量占总能量的比值η_i，并找出能量比的最大值η(z_i)_max来反映人耳对于时变信号的总体感知情况；接下来求取时变特征响度频带能量信息熵，对其进行多项式拟合得到信息熵权重系数，并将权重系数和时变特征响度进行运算得到权重下的时变特征响度；接着对权重下的时变特征响度进行积分运算得到权重下的时变响度数据，对得到的权重下的时变响度数据作傅里叶变换，计算频谱图中[c,d]频带内的峰均比e；最后将最大值η(z_i)_max、对应频带号i、峰均比e和频点f₁形成特征向量m，并将m作为逻辑回归的输入，建立逻辑回归模型进行训练，实现电动冲击批零件缺损的分类。

本发明不局限于以上所述的具体实施方式，以上所述仅为本发明的较佳实施案例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (5)

1.一种基于听觉感知的电动冲击批零件缺损的检测方法，其特征在于：包括以下步骤：

S1：使用传声器采集电动冲击批打滑状态的声音信号，设传声器的采样频率为f_s，单位为Hz；声音信号时长为T，单位为s；按照德国DIN45631/A1标准对声音信号计算各个频带的时变特征响度v(n,z_i)，其中n是指将信号时长T按照2ms为分度划分的点数；

i是指将24个临界频带Bark按照0.1Bark为分度划分为240个子频带，i为对应的频带号，i＝1,2,…,240；z_i为心理声学的24个临界频带Bark的频率尺度；

S2：计算时变特征响度v(n,z_i)各频带的能量占整个频带的总能量的比值η(z_i)以及η(z_i)中最大值η(z_i)_max，得到η(z_i)_max对应的频带号i：

η(z_i)_max＝max[η(z_i)]；(4)

其中N(z_i)是指第i个子频带的能量，N_total是指整个频带的总能量；

S3：对各频带的时变特征响度v(n,z_i)作离散傅里叶变换DFT处理后得r(k_z,z_i)，求取经离散傅里叶变换DFT处理后的时变特征响度r(k_z,z_i)各频带的能量信息熵E(z_i)：

其中，k_z是指频率信息；

为向下取整运算符；

S4：对时变特征响度r(k_z,z_i)各频带的能量信息熵E(z_i)进行多项式拟合，得到信息熵权重系数g(z_i)；

g(z_i)＝a₀+a₁*E(z_i)+a₂*E(z_i)^2+a₃*E(z_i)^3+a₄*E(z_i)^4+a₅*E(z_i)^5； (7)

其中a＝[a₀,a₁,a₂,a₃,a₄,a₅]为多项式拟合系数；

S5：利用人耳听觉的非线性作用，求取信息熵权重下的时变特征响度s(n,z_i)，即各频带时变特征响度乘以对应频带的信息熵权重系数：

s(n,z_i)＝v(n,z_i)·g(z_i)； (8)

S6：对s(n,z_i)进行运算，得到信息熵权重下的时变响度x(n)；

S7：对x(n)作离散傅里叶变换DFT得到y(f)；并计算|y(f)|中[c,d]频带内的峰均比e，具体方法为：

查找[c,d]频带内|y(f)|的幅值最大值A₁和对应频率f₁，分别计算下临界频点f₂和上临界频点f₃为：

|y(f)|是幅值谱，表征信号的幅值随频率f的分布情况，Δf是扩展常数，表征最大值所在波峰的统计半宽；

计算频带[c,f₂]和[f₃,d]内|y(f)|的平均值A₂，则峰均比e为：

S8：将各频带的能量占整个频带的总能量的比值中的最大值η(z_i)_max、对应频带号i、峰均比e和频率f₁构造输入特征向量m＝[η(z_i)_max,i,e,f₁]^T；

S9：重复步骤S1-S8，得到若干个特征向量，划分为训练样本集和测试样本集；

S10：利用逻辑回归分类模型，将训练样本集中的各个特征向量m作为逻辑回归分类模型的输入，对逻辑回归分类模型进行训练；

S11：将测试样本集中的各个特征向量m送到已经训练好的逻辑回归分类模型中进行分类，当得到的预测概率小于等于0.5时，即认为该样本为零件缺损样本；当得到的预测概率是大于0.5，即认为该样本为合格样本。

2.根据权利要求1所述的一种基于听觉感知的电动冲击批零件缺损的检测方法，其特征在于：所述步骤S7中Δf的选取方法为：

(1)查找测试样本集中每个合格产品的声音信号的整个幅值谱|y(f)|的峰值对应的频率点f_0i，并找到峰值左右两边的第一个峰谷的对应频率，设为f_0il和f_0ir；i表示第i个合格产品；

(2)计算第i个合格产品的声音信号的波峰宽度δ_i＝f_0ir-f_0il；

其中，δ_i表示第i个合格产品的声音信号的整个波峰的宽度；

(3)计算测试样本集中采集的每个合格产品的声音信号的波峰宽度，并对计算得到的波峰宽度作统计平均，具体如下：

其中，samples表示测试样本集中的合格产品数；

(4)则

其中

表示向上取整。

3.根据权利要求2所述的一种基于听觉感知的电动冲击批零件缺损的检测方法，其特征在于：所述频带[c,d]的确定方法为：

(1)对测试样本集中合格产品的声音信号的峰值频率f_0i做统计平均，得到测试样本集中合格产品的声音信号的峰值平均频率f₀_average，具体如下：

(2)以测试样本集中合格产品的声音信号的峰值平均频率f₀_average为中心，分别向左和向右各扩展2～5Δf即可得到频带[c,d]。

4.根据权利要求1所述的一种基于听觉感知的电动冲击批零件缺损的检测方法，其特征在于：所述步骤S9中的逻辑回归分类模型为：

其中θ＝(θ₀,θ₁,θ₂,…,θ_t)为逻辑回归分类模型中的参数向量，其中y＝{0，1}代表样本类别，其中1代表合格类产品，0代表零件缺损类产品，x为特征向量，x＝m；t为特征向量的个数，T表示向量转置符号，P(y|x；θ)代表样本的预测概率值，P(y＝1|x；θ)代表合格类产品的预测概率值，P(y＝0|x；θ)代表零件缺损类产品的预测概率值。

5.根据权利要求4所述的一种基于听觉感知的电动冲击批零件缺损的检测方法，其特征在于：采用对数似然函数作为逻辑回归分类模型的损失函数J(θ)，具体如下：

单个样本的损失函数：

cost(h_θ(x),y)＝-y_jlog(h_θ(x))-(1-y_j)log(1-h_θ(x))； (18)

其中h_θ(x)是指P(y|x；θ)，j是指样本序号，y_j是指第j个样本；

全体样本的损失函数：

其中j是指样本序号，y_j是指第j个样本，x_j是第j个样本的特征向量；

当全体样本的损失函数的值无限趋近0的时候，则说明模型最优，即所有样本都预测准确，则逻辑回归分类模型训练完毕。

Images