CN111310569A - 基于加工表面二维变分模态分解的空间频率辨识方法 - Google Patents

基于加工表面二维变分模态分解的空间频率辨识方法 Download PDF

Info

Publication number
CN111310569A
CN111310569A CN202010047284.5A CN202010047284A CN111310569A CN 111310569 A CN111310569 A CN 111310569A CN 202010047284 A CN202010047284 A CN 202010047284A CN 111310569 A CN111310569 A CN 111310569A
Authority
CN
China
Prior art keywords
dimensional
decomposition
data
spatial frequency
hanning window
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202010047284.5A
Other languages
English (en)
Other versions
CN111310569B (zh
Inventor
李星占
高炜祥
冯艳冰
李加胜
魏巍
陈刚利
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Mechanical Manufacturing Technology of CAEP
Original Assignee
Institute of Mechanical Manufacturing Technology of CAEP
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Mechanical Manufacturing Technology of CAEP filed Critical Institute of Mechanical Manufacturing Technology of CAEP
Priority to CN202010047284.5A priority Critical patent/CN111310569B/zh
Publication of CN111310569A publication Critical patent/CN111310569A/zh
Application granted granted Critical
Publication of CN111310569B publication Critical patent/CN111310569B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2218/00—Aspects of pattern recognition specially adapted for signal processing
    • G06F2218/08—Feature extraction
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques
    • G06F18/241—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches
    • G06F18/2415—Classification techniques relating to the classification model, e.g. parametric or non-parametric approaches based on parametric or probabilistic models, e.g. based on likelihood ratio or false acceptance rate versus a false rejection rate
    • G06F18/24155—Bayesian classification

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Computational Mathematics (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Artificial Intelligence (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Evolutionary Computation (AREA)
  • Signal Processing (AREA)
  • Computing Systems (AREA)
  • Probability & Statistics with Applications (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Evolutionary Biology (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Discrete Mathematics (AREA)
  • Length Measuring Devices By Optical Means (AREA)

Abstract

本发明公开了基于加工表面二维变分模态分解的空间频率辨识方法；提出利用镜像延拓和自卷积Hanning窗对采集到的三维形貌数据进行预处理。然后，采用二维变分模态分解方法对表面形貌空间频率误差进行分解。并且，针对BVMD中对分解结果影响最大的两个参数惩罚因子α和分解层数K，采用粒子群退火优化算法进行了寻优；通过对实际加工表面的分析，对比离散小波分解和BEMD算法，以模态混叠作为指标，验证了本文算法的优势及其适用性；本发明提出了一种基于超精密加工表面形貌二维变分模态分解的空间频率辨识方法；解决了现有分解方法中存在的模态混叠、重构误差大的问题，利于有效解决超精密加工表面空间频率分解不彻底，分解产生信号失真的问题。

Description

基于加工表面二维变分模态分解的空间频率辨识方法
技术领域
本发明涉及超精密加工领域，具体地，涉及基于加工表面二维变分模态分解的空间频率 辨识方法。
背景技术
超精密加工是光学元件加工的主要手段，以单点金刚石飞切加工为例，机床主轴摆动、 刀具振动、液压和气浮系统压力波动、环境扰动等因素会导致加工工件表面产生不同空间频 段的形貌误差。在光学研究领域，光学元件表面空间频率误差严重影响其光学性能。高频误 差影响光学元件散射损耗和薄膜损伤性能；中频误差导致光线的小角度散射；低频误差主要 影响激光束聚焦性能。针对不同光学性能的研究，需要对特定频段误差进行分析。因此，准 确分离并辨识表面形貌中的空间频率成分是深入研究和提高工件光学性能的重要前提。
当前研究中，离散小波分析是实现表面形貌空间频率分解与辨识的主要手段之一。但是， 该方法需要对小波基以及小波分解层数做出选择，不同小波基和分解层数对分解效果会有较 大差异。此外，经验模态分析也在加工表面三维形貌的分解与辨识中得到了应用。该方法理 论上适用于任何类型的信号，并能够自适应地将复杂信号分解为多个本征模态分解，避开了 小波分解中对小波基和分解层数选择的需求。但是，该方法存在的严重模态混叠现象，不能 将不同空间频率误差进行有效分离。变分模态分解具有可靠的理论基础，并且能够针对传统 分解方法的模态混沌现象做出较大改善。但是，该方法由于FFT和IFFT的存在，会引入较 大重构误差，导致三维形貌失真。因此，为了更准确的对特定频段的空间频率进行分析，需 要采用有效分解手段，对表面空间频率做出分解，尽可能精确地辨识各频段的频率成分与幅 值。一个有效分解与辨识方法的提出，对加工表面空间频率的分析至关重要。
发明内容
本发明所要解决的技术问题是：现有超精密加工表面空间频率分解与辨识方法存在模态 混叠及重构误差较大的缺陷，导致加工表面空间频率不能得到较好分离，并可能导致分解后 的数据在原始数据上产生失真。本发明提供了解决上述问题的一种基于超精密加工表面形貌 二维变分模态分解的空间频率辨识方法，有效解决了现有方法模态混叠严重导致空间频率分 离效果不理想，重构误差大导致原始信号失真等问题。
本发明通过下述技术方案实现：
一种基于超精密加工表面形貌二维变分模态分解的空间频率辨识方法，包括以下步骤：
根据采集加工后工件三维形貌数据的大小，应用逐步延拓方法对三维形貌数据进行延拓， 实现截断误差初处理；采用二维自卷积Hanning窗函数对延拓后数据进行加窗处理，在减小 截断误差的同时改善延拓数据与原始数据变化趋势不同的缺陷；采用粒子群退火优化算法对 二维变分模态分解算法中的惩罚参数α以及分解层数K进行寻优；在算法寻优的过程中，利 用最小风险贝叶斯决策理论，将分解形貌之间的频谱KL散度与分解形貌的重构误差权衡， 得到优化算法中的适应度函数；采用优化的二维变分模态分解算法对采集三维形貌数据进行 分析，实现超精密加工表面空间频率的低混叠分离与辨识。
进一步地，所述逐步延拓方法，包括以下步骤：
步骤a，将超精密加工后的工件面形数据记为A0(x,y)，其中，x、y分别为超精密加工后的 工件面形数据的行和列的采样点，(x＝1,2,3,…,M；y＝1,2,3,…,N)；
步骤b，对A0(x,y)以其下边界为轴线，进行镜像翻转，将翻转后数据与A0(x,y)结合得到第 一次延拓结果A1(x,y)，其中，(x＝1,2,3,…,2M；y＝1,2,3,…,N)；
步骤c，对A1(,x)y以其右边界为轴线，进行镜像翻转，将翻转后数据与A1(x,y)结合得到第二 次延拓结果A2(x,y)，其中，(x＝1,2,3,…,2M；y＝1,2,3,…,2N)；
步骤d，选择A2(x,y)上边界为轴线，以(x＝1,2,3,…,M；y＝1,2,3,…,2N)范围内的数据为对 象，进行镜像翻转，将翻转后数据与A2(x,y)结合得到第三次延拓结果A3(x,y)，其中，(x＝-M,…,1,…,2M；y＝1,2,3,…,2N)；
步骤e，选择A3(x,y)左边界为轴线，以(x＝-M,…,1,…,2M；y＝1,2,3,…,N)范围内的数据 为对象，进行镜像翻转，将翻转后数据与A3(x,y)结合得到第四次延拓结果A4(x,y)，其中， (x＝-M,…,1,…,2M；y＝-N,…,1,…,2N)。
进一步地，所述加窗处理方法，包括以下步骤：
步骤a，根据延拓数据大小，确定初始Hanning窗函数，其表达式为
Figure BDA0002369887910000023
步骤b，对Hanning窗函数进行卷积计算，以2阶自卷积Hanning窗函数为例，其计算表 达式为
Figure BDA0002369887910000024
步骤c，轴向旋转法实现一维Hanning窗函数到二维Hanning窗函数转换；
步骤d，归一化自卷积Hanning窗函数，并对延拓数据加窗处理。
进一步地，窗函数的频谱形状对频谱泄漏等具有直接影响，主要考察指标为窗函数的旁 瓣特性。采用自卷积方法计算得到的窗函数具有比原窗函数更好的旁瓣特性。随着自卷积阶 数增大，虽然窗函数的旁瓣特性增强，但过高自卷积阶数的窗函数会造成原始数据失真。因 此，针对不同采集三维形貌数据大小，需要根据去旁瓣特性和原始数据的真实性进行择优选 取自卷积阶数。
进一步地，所述粒子群退火优化二维变分模态分解方法，包括以下步骤：
步骤a，输入三维形貌数据X，初始化参数(局部搜索能力c1，全局搜索能力c2，最大种群数size，速度系数k，速率弹性系数wV，种群弹性系数wP，最大最小分解层数Kmax、Kmin，最大最小惩罚系数αmax、αmin)；
步骤b，产生初始粒子pop和速度V，计算初始适应度值fitness，并初始化温度T；
步骤c，更新速度与种群，进行个体最优更新和群体最优更新，然后随机更新速度和种 群，再次对个体最优更新和群体最优更新，最后更新温度；
步骤d，循环进行步骤c，直至温度更新至T＜Tmin，输出最优分解层数bestK和最优惩罚 系数bestα。
进一步地，所述适应度函数建立方法，包括以下步骤：
步骤a，利用FFT计算超精密加工表面三维形貌经二维变分模态分解算法分解得到的各 分解成分的频谱；
步骤b，根据各分解成分的频谱，计算相邻形貌频率之间的KL散度，KL散度的计算表 达式为
Figure BDA0002369887910000031
由于三维形貌数据及其频谱数据皆为矩阵形式，因此将上述表达式转化为矩阵形式为
Figure BDA0002369887910000032
步骤c，计算二维变分模态分解后的重构误差，并取均值；
步骤d，引入最小风险贝叶斯决策理论对KL散度和重构误差进行权衡。最小风险贝叶斯 决策中，某一决策规则的期望损失表达式为
R(α)＝∫R(α(x)|x)p(x)dx；
其中，R(α(x)|x)表示每个决策的期望损失；p(x)表示总体概率密度。表达式分别如下 所示
Figure BDA0002369887910000041
Figure BDA0002369887910000042
其中，c表示类别总数；λ(αi,ωi)表示决策αi对类别ωj的决策损失；P(ωj|x)为ωj在条 件x下，所属各类别的概率，即贝叶斯公式求解结果，如下文公式所示；p(x|ωj)为类条件密 度，表示ωj在x条件方面的概率密度；P(ωi)表示所属各类别的先验概率。
Figure BDA0002369887910000043
步骤e，根据步骤d中的权衡结果，求得适应度函数。
进一步地，所述步骤d中，ω1表示重构误差，ω2表示模态混叠，λ(α,ω1)为重构误差值 均值，λ(α,ω2)为KL散度值均值的倒数。由于无条件干涉，先验概率P(ω1)、P(ω2)都设定 为0.5。采取优化措施时需要在保证高精度重构误差的基础上进行。另外，由于二维变分模态 分解在模态混叠上的改善，以及该算法会在一定程度上影响重构误差精度，模态混叠在评价 指标中的占比应远小于重构误差。通过对多组预实验分析，总结得出，KL散度倒数值比重构 误差值在数量级上大10-9倍。由此，拟定重构误差和模态混叠的影响力概率密度p(x|ω1)为 1-10-9、p(x|ω2)为10-9。
本发明具有如下的优点和有益效果：
本发明提供了一种基于超精密加工表面形貌二维变分模态分解的空间频率辨识方法，能 够有效分离与辨识超精密加工表面的空间频率，解决了现有分解方法中存在的模态混叠、重 构误差大的问题，利于有效解决超精密加工表面空间频率分解不彻底，分解产生信号失真的 问题。
本发明针对超精密加工表面空间频率分解中存在的模态混叠和重构误差大的问题，提出 了有效的抑制模态混叠和降低重构误差的方法，有效地消除了现有分解方法中模态混叠和失 真问题，从而提高了分解后对面形评价的准确性。采用二维的变分模态分解对工件表面三维 形貌进行分解。但是，该算法存在两个对分解效果影响较大的参数，惩罚因子α和分解层数K。为了得到较好分解效果，采用粒子群退火算法对以上两个参数进行优化处理。优化算法 的性能不仅与算法理论有关，还与所优化的适应度函数有较大关系。根据空间频率误差分解 中对模态混叠以及重构误差的要求，将KL散度作为模态混叠的数值指标，并引入最小风险 贝叶斯理论权衡模态混叠与重构误差在适应度函数中的占比。三维形貌数据采样是对离散点 进行采集，存在截断误差。二维变分模态分解算法中多次采用了FFT和IFFT算法，截断误 差对这两种算法的精度影响很大。因此，本文采用边界延拓以及Hanning自卷积窗对采集数 据进行预处理，降低分析误差。
附图说明
此处所说明的附图用来提供对本发明实施例的进一步理解，构成本申请的一部分，并不 构成对本发明实施例的限定；
图1为本发明超精密加工表面空间频率误差分解方法流程图；
图2为本发明粒子群退火优化算法流程图；
图3为本发明三维形貌数据延拓示意图；
图4为本发明超精密加工实验三维形貌数据1面形图；
图5为本发明对三维形貌数据1延拓后的面形图；
图6为本发明超精密加工实验三维形貌数据1频谱图；
图7为本发明构造的二维Hanning窗示意图；
图8a-图8d为本发明对三维形貌数据1三维形貌延拓数据加窗后，原始数据范围对应面 形图，图8a为1阶Hanning窗处理，图8b为2阶Hanning窗处理，图8c为3阶Hanning 窗处理，图8d为4阶Hanning窗处理；
图9a-图9d为本发明对三维形貌数据1采用二维变分模态分解得到的分解面形图；
图10a-图10d为本发明对三维形貌数据1采用二维变分模态分解得到的分解频谱图；
图11为本发明超精密加工实验三维形貌数据2面形图；
图12为本发明超精密加工实验三维形貌数据2频谱图；
图13a-图13e为本发明对三维形貌数据2采用二维变分模态分解得到的分解面形图；
图14a-图14e为本发明对三维形貌数据2采用二维变分模态分解得到的分解频谱图；
图15a-图15e为本发明对三维形貌数据1采用二维离散小波分解得到的分解面形图；
图16a-图16e为本发明对三维形貌数据1采用二维离散小波分解得到的分解频谱图；
图17a-图17f为本发明对三维形貌数据1采用二维经验模态分解得到的分解面形图；
图18a-图18f为本发明对三维形貌数据1采用二维经验模态分解得到的分解频谱图。
具体实施方式
为了能够更清楚地理解本发明的上述目的、特征和优点，下面结合附图和具体实施方式 对本发明进行进一步的详细描述。需要说明的是，在相互不冲突的情况下，本申请的实施例 及实施例中的特征可以相互组合。
在下面的描述中阐述了很多具体细节以便于充分理解本发明，但是，本发明还可以采用 其他不同于在此描述范围内的其他方式来实施，因此，本发明的保护范围并不受下面公开的 具体实施例的限制。
本领域技术人员应理解的是，在本发明的揭露中，术语“纵向”、“横向”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、 “外”等指示的方位或位置关系是基于附图所示的方位或位置关系，其仅是为了便于描述本 发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位 构造和操作，因此上述术语不能理解为对本发明的限制。
可以理解的是，术语“一”应理解为“至少一”或“一个或多个”，即在一个实施例中，一个元件的数量可以为一个，而在另外的实施例中，该元件的数量可以为多个，术语“一”不能理解为对数量的限制。
实施例1
本实施例提供了一种基于超精密加工表面形貌二维变分模态分解的空间频率辨识方法， 具体步骤如下所示：
步骤1，将采集三维形貌数据转换为矩阵形式，并确定矩阵大小，其初始三维形貌如图4 所示。
步骤2，对步骤1中的矩阵数据进行逐步延拓，其延拓示意图如图3所示，延拓后三维 形貌数据如图5所示：
步骤21，将超精密加工后的工件面形数据记为A0(x,y)，其中，x、y分别为超精密加工后 的工件面形数据的行和列的采样点，(x＝1,2,3,…,M；y＝1,2,3,…,N)；
步骤22，对A0(x,y)以其下边界为轴线，进行镜像翻转，将翻转后数据与A0(x,y)结合得到第 一次延拓结果A1(x,y)，其中，(x＝1,2,3,…,2M；y＝1,2,3,…,N)；
步骤23，对A1(x,y)以其右边界为轴线，进行镜像翻转，将翻转后数据与A1(x,y)结合得到第 二次延拓结果A2(x,y)，其中，(x＝1,2,3,…,2M；y＝1,2,3,…,2N)；
步骤24，选择A2(x,y)上边界为轴线，以(x＝1,2,3,…,M；y＝1,2,3,…,2N)范围内的数据为 对象，进行镜像翻转，将翻转后数据与A2(x,y)结合得到第三次延拓结果A3(x,y)，其中，(x＝-M,…,1,…,2M；y＝1,2,3,…,2N)；
步骤25，选择A3(x,y)左边界为轴线，以(x＝-M,…,1,…,2M；y＝1,2,3,…,N)范围内的数据 为对象，进行镜像翻转，将翻转后数据与A3(x,y)结合得到第四次延拓结果A4(x,y)，其中， (x＝-M,…,1,…,2M；y＝-N,…,1,…,2N)。
步骤3，根据步骤2中延拓数据大小，确定一维窗函数的长度，然后设计二维自卷积Hanning窗函数对延拓后数据进行加窗处理。其中，1～4阶自卷积Hanning窗函数的处理后的 三维形貌(原始数据范围)如图8a-图8d所示。
步骤31，建立与延拓数据长度相同的一维Hanning窗函数。
步骤32，采用轴向旋转法，获得中心对称的二维曲面形式的二维Hanning窗，如图7所 示。
步骤33，对步骤32中的二维Hanning窗进行归一化处理，然后将待处理的延拓数据点 乘Hanning窗函数矩阵的转置，即获得加窗后延拓数据。
步骤34，采用步骤33方法，获得1～4阶自卷积Hanning窗处理后的延拓数据，并分析 原始数据范围对应的三维形貌数据变化。由图8a-图8d可以看出，随着自卷积阶数增大，原 始数据范围对应三维形貌的边缘数据逐渐弱化。因此，为了在保证原始数据的真实性的前提 下获得较好的截断误差改善效果，本实施例中选择3阶自卷积Hanning窗对延拓数据进行处 理。
步骤4，采用粒子群退火优化的二维变分模态分解方法对数据1和数据2进行分解。其 中，粒子群退火优化二维变分模态分解如图2所示。
步骤41，输入经延拓及加窗处理后的三维形貌数据X，初始化参数(局部搜索能力c1， 全局搜索能力c2，最大种群数size，速度系数k，速率弹性系数wV，种群弹性系数wP，最大最小分解层数Kmax、Kmin，最大最小惩罚系数αmax、αmin)；
步骤42，产生初始粒子pop和速度V，计算初始适应度值fitness，并初始化温度T。其 中，适应度值由各分解成分之间的频谱KL散度和重构误差综合表示，具体如下：
(1)利用FFT计算超精密加工表面三维形貌经二维变分模态分解算法分解得到的各分解 成分的频谱；
(2)根据各分解成分的频谱，计算相邻形貌频率之间的KL散度，KL散度的计算表达式 为
Figure BDA0002369887910000081
由于三维形貌数据及其频谱数据皆为矩阵形式，因此将上述表达式转化为矩阵形式为
Figure BDA0002369887910000082
(3)计算二维变分模态分解后的重构误差，并取均值；
(4)引入最小风险贝叶斯决策理论对KL散度和重构误差进行权衡。最小风险贝叶斯决策 中，某一决策规则的期望损失表达式为
R(α)＝∫R(α(x)|x)p(x)dx；
其中，R(α(x)|x)表示每个决策的期望损失；p(x)表示总体概率密度。表达式分别如下 所示
Figure BDA0002369887910000083
Figure BDA0002369887910000084
其中，c表示类别总数；λ(αi,ωi)表示决策αi对类别ωj的决策损失；P(ωj|x)为ωj在条 件x下，所属各类别的概率，即贝叶斯公式求解结果，如下文公式所示；p(x|ωj)为类条件密 度，表示ωj在x条件方面的概率密度；P(ωi)表示所属各类别的先验概率。
Figure BDA0002369887910000085
(5)根据步骤d中的权衡结果，求得适应度函数。
进一步地，所述步骤d中，ω1表示重构误差，ω2表示模态混叠，λ(α,ω1)为重构误差值 均值，λ(α,ω2)为KL散度值均值的倒数。由于无条件干涉，先验概率P(ω1)、P(ω2)都设定 为0.5。采取优化措施时需要在保证高精度重构误差的基础上进行。另外，由于二维变分模态 分解在模态混叠上的改善，以及该算法会在一定程度上影响重构误差精度，模态混叠在评价 指标中的占比应远小于重构误差。通过对多组预实验分析，总结得出，KL散度倒数值比重构 误差值在数量级上大10-9倍。由此，拟定重构误差和模态混叠的影响力概率密度p(x|ω1)为 1-10-9、p(x|ω2)为10-9。
步骤43，更新速度与种群，进行个体最优更新和群体最优更新，然后随机更新速度和种 群，再次对个体最优更新和群体最优更新，最后更新温度；
步骤44，循环进行步骤c，直至温度更新至T＜Tmin，输出最优分解层数bestK和最优惩 罚系数bestα。
步骤5，通过步骤4方法寻找最优bestK和最优惩罚系数bestα，实现最终数据1(图4) 和数据2(图11)，其分解结果展示如图9a-图9d和图13a-图13e，分解频谱图分别展示于图 10a-图10d和图14a-图14e。可以看出，二维变分模态分解能够较好地根据空间频率误差将三 维形貌进行分解。通过图9a-图9d分析得出，原始三维形貌中所存在的主要空间频率误差被 清晰地划分到了四个分量之中，并且几乎没有模态混叠存在。因此，本发明方法对数据1中 的主要频率误差实现了较为有效地分离。另外，从数据2的分解结果分析得出，本发明方法 能够有效地将各主要频率误差进行分离，并基本不存在模态混叠。因此，本发明方法在保证 分解效果的前提下，对实测超精密加工表面空间频率误差分解与辨识具有普遍适用性。
步骤6，采用传统二维离散小波分解和二维经验模态分解方法对数据1进行分解，其结 果展示于图15a-图15e和图17a-图17f，频谱图分别展示于图16a-图16e和图18a-图18f。其 中，方括号对应坐标数据表示存在混叠现象的主要空间频率。图16a-图16e为BDWT分解所 得各形貌的一维频谱图，可以看出，各分解成分中存在较大混叠现象，主要频率误差不能有 效分离。图18a-图18f为BEMD分解所得各形貌的一维频谱图，分析得到，该方法分解结果 也存在较大模态混叠。并且，最后三个分量虽然主要频率一样，但再椭圆圈出的部分存在差 异，因此也属于模态混叠。更重要的是，该方法在中高频误差成分的分解效果很差，大量中 高频误差存在一个分解成分中，没有达到理想分解效果。
通过上述分析，本实施例实现了超精密加工表面形貌二维变分模态分解的空间频率辨识。
尽管已描述了本发明的优选实施例，但本领域内的技术人员一旦得知了基本创造性概念， 则可对这些实施例作出另外的变更和修改。所以，所附权利要求意欲解释为包括优选实施例 以及落入本发明范围的所有变更和修改。
显然，本领域的技术人员可以对本发明进行各种改动和变型而不脱离本发明的精神和范 围。这样，倘若本发明的这些修改和变型属于本发明权利要求及其等同技术的范围之内，则 本发明也意图包含这些改动和变型在内。

Claims (9)

1.基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，所述方法包括：
获得工件加工后表面三维形貌数据；对工件加工后表面三维形貌数据进行延拓；采用二维自卷积Hanning窗函数对延拓后数据进行加窗处理；采用粒子群退火优化算法对二维变分模态分解算法中的惩罚参数α以及分解层数K进行寻优；在二维变分模态分解算法寻优的过程中，引入最小风险贝叶斯决策理论，对分解形貌之间的频谱KL散度与分解形貌的重构误差进行权衡，得到优化二维变分模态分解算法中的适应度函数；采用优化后二维变分模态分解算法对延拓和加窗后的三维形貌数据进行分解，实现工件加工表面空间频率的低混叠分离及辨识。
2.根据权利要求1所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，采用逐步延拓方法对工件加工后表面三维形貌数据进行延拓。
3.根据权利要求2所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，采用逐步延拓方法对工件加工后表面三维形貌数据进行延拓，具体包括：
将加工后的工件面形数据记为A0(x,y)，其中，x、y分别为加工工件面形数据的行和列的采样点，(x＝1,2,3,…,M；y＝1,2,3,…,N)；
对A0(x,y)以其下边界为轴线，进行镜像翻转，将翻转后数据与A0(x,y)结合得到第一次延拓结果A1(x,y)，其中，(x＝1,2,3,…,2M；y＝1,2,3,…,N)；
对A1(x,y)以其右边界为轴线，进行镜像翻转，将翻转后数据与A1(x,y)结合得到第二次延拓结果A2(x,y)，其中，(x＝1,2,3,…,2M；y＝1,2,3,…,2N)；
选择A2(x,y)上边界为轴线，以(x＝1,2,3,…,M；y＝1,2,3,…,2N)范围内的数据为对象，进行镜像翻转，将翻转后数据与A2(x,y)结合得到第三次延拓结果A3(x,y)，其中，(x＝-M,…,1,…,2M；y＝1,2,3,…,2N)；
选择A3(x,y)左边界为轴线，以(x＝-M,…,1,…,2M；y＝1,2,3,…,N)范围内的数据为对象，进行镜像翻转，将翻转后数据与A3(x,y)结合得到第四次延拓结果A4(x,y)，其中，(x＝-M,…,1,…,2M；y＝-N,…,1,…,2N)。
4.根据权利要求2所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，所述加窗处理具体包括：
根据延拓数据大小，确定初始Hanning窗函数；
对Hanning窗函数进行卷积计算；
轴向旋转法实现一维Hanning窗函数到二维Hanning窗函数转换；
归一化自卷积Hanning窗函数，并对延拓数据加窗处理。
5.根据权利要求4所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，
根据延拓数据大小，确定初始Hanning窗函数，其表达式为：
Figure FDA0002369887900000021
对Hanning窗函数进行卷积计算，以2阶自卷积Hanning窗函数为例，其计算表达式为：
Figure FDA0002369887900000022
轴向旋转法实现一维Hanning窗函数到二维Hanning窗函数转换；
归一化自卷积Hanning窗函数，并对延拓数据加窗处理。
6.根据权利要求5所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，针对三维形貌数据的大小，根据旁瓣特性和原始数据的真实性择优选取自卷积阶数。
7.根据权利要求1所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，粒子群退火优化二维变分模态分解方法，包括以下步骤：
输入三维形貌数据X，初始化参数，包括：局部搜索能力c1，全局搜索能力c2，最大种群数size，速度系数k，速率弹性系数wV，种群弹性系数wP，最大最小分解层数Kmax、Kmin，最大最小惩罚系数αmax、αmin；
产生初始粒子pop和速度V，计算初始适应度值fitness，并初始化温度T；
更新速度与种群，进行个体最优更新和群体最优更新，然后随机更新速度和种群，再次对个体最优更新和群体最优更新，最后更新温度；
循环进行更新步骤，直至温度更新至T＜Tmin，输出最优分解层数bestK和最优惩罚系数bestα。
8.根据权利要求1所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，适应度函数的建立方法包括以下步骤：
利用二维变分模态分解算法对超精密加工表面三维形貌进行分解，并基于FFT计算各分解成分的频谱；
根据各分解成分的频谱，计算相邻形貌之间的频率KL散度并求均值，KL散度的计算表达式为：
Figure FDA0002369887900000031
将上述表达式转化为矩阵形式为：
Figure FDA0002369887900000032
计算二维变分模态分解后的重构误差，并取均值；
引入最小风险贝叶斯决策理论对KL散度和重构误差进行权衡，最小风险贝叶斯决策中，某一决策规则的期望损失表达式为：
R(α)＝∫R(α(x)|x)p(x)dx；
其中，R(α(x)|x)表示每个决策的期望损失；p(x)表示总体概率密度，表达式分别为：
Figure FDA0002369887900000033
Figure FDA0002369887900000034
其中，c表示类别总数；λ(αi,ωi)表示决策αi对类别ωj的决策损失；P(ωj|x)为ωj在条件x下，所属各类别的概率，即贝叶斯公式求解结果，如下文公式所示；p(x|ωj)为类条件密度，表示ωj在x条件方面的概率密度；P(ωi)表示所属各类别的先验概率；
Figure FDA0002369887900000035
根据权衡结果，求得适应度函数。
9.根据权利要求8所述的基于加工表面二维变分模态分解的空间频率辨识方法，其特征在于，ω1表示重构误差，ω2表示模态混叠，λ(α,ω1)为重构误差值均值，λ(α,ω2)为KL散度值均值的倒数；先验概率P(ω1)、P(ω2)设定为0.5；重构误差和模态混叠的影响力概率密度p(x|ω1)为1-10-9、p(x|ω2)为10-9。
CN202010047284.5A 2020-01-16 2020-01-16 基于加工表面二维变分模态分解的空间频率辨识方法 Active CN111310569B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010047284.5A CN111310569B (zh) 2020-01-16 2020-01-16 基于加工表面二维变分模态分解的空间频率辨识方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010047284.5A CN111310569B (zh) 2020-01-16 2020-01-16 基于加工表面二维变分模态分解的空间频率辨识方法

Publications (2)

Publication Number Publication Date
CN111310569A true CN111310569A (zh) 2020-06-19
CN111310569B CN111310569B (zh) 2022-04-26

Family

ID=71148866

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010047284.5A Active CN111310569B (zh) 2020-01-16 2020-01-16 基于加工表面二维变分模态分解的空间频率辨识方法

Country Status (1)

Country Link
CN (1) CN111310569B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117726133A (zh) * 2023-12-29 2024-03-19 国网江苏省电力有限公司信息通信分公司 一种基于强化学习的分布式能源实时调度方法及系统

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070240548A1 (en) * 2006-04-12 2007-10-18 Xpel Technologies Corporation Cutting pattern manipulation and methods
CN103367085A (zh) * 2012-04-05 2013-10-23 Fei公司 提供深度分辨图像的带电粒子显微镜
CN104887266A (zh) * 2015-05-29 2015-09-09 西安交通大学 基于面阵的小区域三维被动空化成像及三维复合成像方法
CN106845792A (zh) * 2016-12-27 2017-06-13 中国电子科技集团公司第五十四研究所 基于星地资源模型的任务规划通用平台
CN108344642A (zh) * 2018-02-02 2018-07-31 中国工程物理研究院机械制造工艺研究所 一种基于重力加载三点弯的蠕变实验装置及测试方法
CN109583621A (zh) * 2018-10-12 2019-04-05 中国农业大学 一种基于改进变分模态分解的pso-lssvm短期负荷预测方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20070240548A1 (en) * 2006-04-12 2007-10-18 Xpel Technologies Corporation Cutting pattern manipulation and methods
CN103367085A (zh) * 2012-04-05 2013-10-23 Fei公司 提供深度分辨图像的带电粒子显微镜
CN104887266A (zh) * 2015-05-29 2015-09-09 西安交通大学 基于面阵的小区域三维被动空化成像及三维复合成像方法
CN106845792A (zh) * 2016-12-27 2017-06-13 中国电子科技集团公司第五十四研究所 基于星地资源模型的任务规划通用平台
CN108344642A (zh) * 2018-02-02 2018-07-31 中国工程物理研究院机械制造工艺研究所 一种基于重力加载三点弯的蠕变实验装置及测试方法
CN109583621A (zh) * 2018-10-12 2019-04-05 中国农业大学 一种基于改进变分模态分解的pso-lssvm短期负荷预测方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
F. PINTO: "Applications of short space-time fourier analysis in digital acoustics", 《2011 IEEE INTERNATIONAL CONFERENCE ON ACOUSTICS, SPEECH AND SIGNAL PROCESSING (ICASSP)》 *
李艳娜: "基于小波变换的三维表面形貌分离技术研究", 《中国优秀博硕士学位论文全文数据库(硕士)信息科技辑》 *
魏巍: "电液伺服振动台的流量非线性补偿控制", 《华南理工大学学报(自然科学版)》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117726133A (zh) * 2023-12-29 2024-03-19 国网江苏省电力有限公司信息通信分公司 一种基于强化学习的分布式能源实时调度方法及系统

Also Published As

Publication number Publication date
CN111310569B (zh) 2022-04-26

Similar Documents

Publication Publication Date Title
US7653646B2 (en) Method and apparatus for quantum clustering
JPWO2006085537A1 (ja) 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体
JP5956624B1 (ja) 異常音の検出方法及びその検出値を用いた構造物の異常判定方法、並びに、振動波の類似度検出方法及びその検出値を用いた音声認識方法
CN106772227A (zh) 一种基于声纹多谐波识别的无人机方向估计方法
CN104849757B (zh) 消除地震信号中随机噪声系统及方法
CN102708373A (zh) 融合空间信息及光谱信息的遥感图像的分类方法及装置
CN109871855B (zh) 一种自适应的深度多核学习方法
CN102096819A (zh) 利用稀疏表示与字典学习进行图像分割的方法
CN111062425B (zh) 基于c-k-smote算法的不平衡数据集处理方法
CN111436926A (zh) 基于统计特征和卷积循环神经网络的房颤信号检测方法
CN114782745B (zh) 一种基于机器学习的海洋声速剖面分类方法及装置
CN111310569B (zh) 基于加工表面二维变分模态分解的空间频率辨识方法
CN103854660A (zh) 一种基于独立成分分析的四麦克语音增强方法
CN112906158A (zh) 一种基于多传感器多元数据融合的机械故障诊断方法
CN115345909B (zh) 基于深度空间光谱卷积融合特征的高光谱目标跟踪方法
CN116619136A (zh) 一种多工况多源数据的刀具磨损预测方法
CN114548272A (zh) 一种离心泵空化状态识别方法
CN107961005B (zh) 少通道脑机接口eeg信号的特征提取方法
CN112001256B (zh) 混合信号去工频干扰方法及系统
CN109940462B (zh) 铣刀切削振动变化特性的检测与高斯过程模型构建方法
JP4630203B2 (ja) 信号分離装置、信号分離方法、信号分離プログラム及び記録媒体、並びに、信号到来方向推定装置、信号到来方向推定方法、信号到来方向推定プログラム及び記録媒体
CN114062306B (zh) 一种近红外光谱数据分段预处理方法
JP7014682B2 (ja) 音源分離の評価装置および音源分離装置
JP2010250670A (ja) 演算装置及びプログラム
Lundberg Classifying dialects using cluster analysis

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant