CN112090994B - 一种基于最小二乘法的轴类工件最大弯曲点的检测方法 - Google Patents

一种基于最小二乘法的轴类工件最大弯曲点的检测方法 Download PDF

Info

Publication number
CN112090994B
CN112090994B CN202010871515.4A CN202010871515A CN112090994B CN 112090994 B CN112090994 B CN 112090994B CN 202010871515 A CN202010871515 A CN 202010871515A CN 112090994 B CN112090994 B CN 112090994B
Authority
CN
China
Prior art keywords
fitting
circle
detection point
particle
point
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.)
Active
Application number
CN202010871515.4A
Other languages
English (en)
Other versions
CN112090994A (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.)
Shaoxing Keqiao Zhejiang University Of Technology Innovation Research Institute Development Co ltd
Original Assignee
Shaoxing Keqiao Zhejiang University Of Technology Innovation Research Institute Development Co ltd
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 Shaoxing Keqiao Zhejiang University Of Technology Innovation Research Institute Development Co ltd filed Critical Shaoxing Keqiao Zhejiang University Of Technology Innovation Research Institute Development Co ltd
Priority to CN202010871515.4A priority Critical patent/CN112090994B/zh
Publication of CN112090994A publication Critical patent/CN112090994A/zh
Application granted granted Critical
Publication of CN112090994B publication Critical patent/CN112090994B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B21—MECHANICAL METAL-WORKING WITHOUT ESSENTIALLY REMOVING MATERIAL; PUNCHING METAL
    • B21D—WORKING OR PROCESSING OF SHEET METAL OR METAL TUBES, RODS OR PROFILES WITHOUT ESSENTIALLY REMOVING MATERIAL; PUNCHING METAL
    • B21D3/00—Straightening or restoring form of metal rods, metal tubes, metal profiles, or specific articles made therefrom, whether or not in combination with sheet metal parts
    • B—PERFORMING OPERATIONS; TRANSPORTING
    • B21—MECHANICAL METAL-WORKING WITHOUT ESSENTIALLY REMOVING MATERIAL; PUNCHING METAL
    • B21C—MANUFACTURE OF METAL SHEETS, WIRE, RODS, TUBES OR PROFILES, OTHERWISE THAN BY ROLLING; AUXILIARY OPERATIONS USED IN CONNECTION WITH METAL-WORKING WITHOUT ESSENTIALLY REMOVING MATERIAL
    • B21C51/00—Measuring, gauging, indicating, counting, or marking devices specially adapted for use in the production or manipulation of material in accordance with subclasses B21B - B21F
    • 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/004—Artificial life, i.e. computing arrangements simulating life
    • G06N3/006—Artificial life, i.e. computing arrangements simulating life based on simulated virtual individual or collective life forms, e.g. social simulations or particle swarm optimisation [PSO]
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Mechanical Engineering (AREA)
  • General Engineering & Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Software Systems (AREA)
  • Computational Mathematics (AREA)
  • Probability & Statistics with Applications (AREA)
  • Biomedical Technology (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Databases & Information Systems (AREA)
  • Operations Research (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Health & Medical Sciences (AREA)
  • Artificial Intelligence (AREA)
  • Algebra (AREA)
  • Biophysics (AREA)
  • Computational Linguistics (AREA)
  • Evolutionary Computation (AREA)
  • General Health & Medical Sciences (AREA)
  • Molecular Biology (AREA)
  • Computing Systems (AREA)
  • Length Measuring Devices With Unspecified Measuring Means (AREA)

Abstract

本发明公开了一种基于最小二乘法的轴类工件最大弯曲点的检测方法，包括：沿轴类工件的轴向预设多个检测点，获取轴类工件周向旋转一周在各检测点处的采样数据；利用最小二乘法，根据获取的所述采样数据计算每个检测点对应的拟合圆的圆心坐标；利用最小二乘法，根据各检测点对应的拟合圆的圆心坐标计算拟合基准轴的直线方程；计算各拟合圆的圆心坐标到拟合基准轴的径向距离，取径向距离最远的圆心坐标对应的检测点作为轴类工件的最大弯曲点。本发明的基于最小二乘法的轴类工件最大弯曲点的检测方法，可精准、快速确定轴类工件的最大弯曲点。

Description

一种基于最小二乘法的轴类工件最大弯曲点的检测方法
技术领域
本申请属于工件校直技术领域，具体涉及一种基于最小二乘法的轴类工件最大弯曲点的检测方法。
背景技术
齿轮、轴杆类零件是制造业中最常见、最常使用的零部件之一。据统计，我国轴杆类零件的年产量高达十多亿件，但是这些零件中大约有70％在生产加工过程中以及热处理后会出现弯曲变形。如果这些零件得不到有效的校直处理，会严重影响产品的质量，导致产品故障率高，甚至会导致废品数量增加，造成材料以及加工成本的浪费。因此，对发生弯曲变形的轴杆类零件进行校直处理在制造行业是必须的过程。
目前，轴杆类零件校直主要包括人工手动校直和自动校直设备校直。人工手动校直过程中使用的设备简单，但对操作者的要求较高，而且校直效率低，校直精度及稳定性难以保证。自动校直设备在校直精度和校直效率上都远远高于人工手动校直，而且能够大大提高生产效率，但是在自动校直设备上，单纯通过检测数据分析比较获得的最大弯曲位置并不准确，行业内亟需寻找一种更有效的获取精准最大弯曲位置的方法。
发明内容
本申请的目的在于提供一种基于最小二乘法的轴类工件最大弯曲点的检测方法，可精准、快速确定轴类工件的最大弯曲点。
为实现上述目的，本申请所采取的技术方案为：
一种基于最小二乘法的轴类工件最大弯曲点的检测方法，所述基于最小二乘法的轴类工件最大弯曲点的检测方法，包括：
步骤S1、沿轴类工件的轴向预设多个检测点，获取轴类工件周向旋转一周在各检测点处的采样数据；
步骤S2、利用最小二乘法，根据获取的所述采样数据计算每个检测点对应的拟合圆的圆心坐标；
步骤S2.1、取一个检测点对应的采样数据，初始化圆拟合的次数；
步骤S2.2、根据该检测点对应的最新的采样数据利用最小二乘法计算拟合圆的圆心坐标，并累加一次圆拟合的次数；
步骤S2.3、根据本次拟合圆的圆心坐标，剔除采样数据中的干扰数据点；
步骤S2.4、判断剔除干扰数据点后剩余的采样数据中的数据点个数，若数据点个数大于个数阈值，则以剔除干扰数据点后的采样数据作为最新的采样数据继续执行步骤S2.2；否则保存本次拟合圆的圆心坐标作为该检测点对应的拟合圆的圆心坐标，并执行步骤S2.5；
步骤S2.5、判断是否完成所有检测点对应拟合圆的圆心坐标计算，若完成则执行步骤S3；否则继续执行步骤S2.1；
步骤S3、利用最小二乘法，根据各检测点对应的拟合圆的圆心坐标计算拟合基准轴的直线方程；
步骤S4、计算各拟合圆的圆心坐标到拟合基准轴的径向距离，取径向距离最远的圆心坐标对应的检测点作为轴类工件的最大弯曲点。
以下还提供了若干可选方式，但并不作为对上述总体方案的额外限定，仅仅是进一步的增补或优选，在没有技术或逻辑矛盾的前提下，各可选方式可单独针对上述总体方案进行组合，还可以是多个可选方式之间进行组合。
作为优选，所述步骤S1，沿轴类工件的轴向预设多个检测点，获取轴类工件周向旋转一周在各检测点处的采样数据，包括：
沿轴类工件的轴向预设N个检测点，则轴类工件的轴向截面数n为N，以Zi表示第i个检测点，即检测点i，其中i＝1,2,3,…,n；
令轴类工件以固定角度旋转，获取每次旋转后各检测点处的采样数据，直至旋转M次后轴类工件旋转一周，得到M次旋转获取的采样数据，所述采样数据包括AD值和角度值；
以每次旋转后每个检测点处获取的AD值和角度值作为一对数据，将所有采样数据保存至数组中，数组中的每个数据表示为(Zi,Rij,θij)，i＝1,2,3,…,n，j＝1,2,3,…,M，Rij表示在第j次旋转后检测点i处获取的AD值，θij表示在第j次旋转后检测点i处获取的角度值。
作为优选，所述步骤S2.2中，根据该检测点对应的最新的采样数据利用最小二乘法计算拟合圆的圆心坐标，包括：
一个检测点对应的采样数据有M对，取检测点i对应的采样数据为(Rij,θij)，i＝1，2，3，…，n，j＝1，2，3，…，M，则可计算得到检测点i对应的坐标Pi(Rijsinθij，Rijcosθij)，记为Pi(xij，yij)；
令检测点i对应的拟合圆的圆心坐标为Oi(Xi，Yi，Zi)，拟合圆的半径为R，则拟合圆的圆方程可表示为：
(xij-Xi)2+(yij-Yi)2＝R2
令残差为ei，并且ei可表示为：
ei＝(xij-Xi)2+(yij-Yi)2-R2
以函数∑ei 2取最小值为目标，采用最小二乘法求解拟合圆的Xi，Yi，如下：
设F(Xi，Yi，R)＝min∑ei 2＝min∑[(xij-Xi)2+(yij-Yi)2-R2]2；
则得到F(Xi，Yi，R)＝min∑[xij 2+yij 2-2Xixij-2Yiyij+Xi 2+Yi 2-R2]2；
令Q(a，b，c)＝∑(xij 2+yij 2+axij+byij+c)2，其中，a＝-2Xi，b＝-2Yi，c＝Xi 2+Yi 2-R2；
式中Q(a，b，c)的偏导数满足最小化的条件为：
Figure BDA0002651271490000031
其中，
Figure BDA0002651271490000035
表示偏导数；
根据Q(a，b，c)的偏导数满足最小化的条件可得：
∑2(xij 2+yij 2+axij+byij+c)Xi＝0
∑2(xij 2+yij 2+axij+byij+c)Yi＝0
∑2(xij 2+yij 2+axij+byij+c)＝0
令A＝N∑xij 2-∑xij∑xij；
B＝N∑xijyij-∑xij∑yij；
C＝N∑xij 3+N∑xijyij 2-∑xij∑(xij 2+yij 2)；
D＝N∑yij 2-∑yij∑yij；
E＝N∑yij 3+N∑xij 2yij-∑yij∑(xij 2+yij 2)；
可得到
Figure BDA0002651271490000032
由此可以解得
Figure BDA0002651271490000033
解得检测点i对应的拟合圆的圆心坐标为
Figure BDA0002651271490000034
作为优选，所述步骤S2.3，根据本次拟合圆的圆心坐标，剔除采样数据中的干扰数据点，包括：
令当前圆拟合的次数为m次，并且该检测点对应的最新的采样数据中每个数据点到本次拟合计算得到的拟合圆的圆心坐标的距离为dij如下：
Figure BDA0002651271490000041
取dij中最远的距离dmmax，设置本次剔除数据的阈值为dm，并且dm＝kdmmax，其中k为预设的比例系数；
剔除该检测点对应的最新的采样数据中dij＞dm的数据点，保留dij≤dm的数据点，并统计剔除干扰数据点后剩余的采样数据中的数据点个数Nm。
作为优选，所述步骤S2.4，判断剔除干扰数据点后剩余的采样数据中的数据点个数，包括：
令个数阈值为Nf，并且Nf＝pM，其中p为预设的比例系数，M为轴类工件以固定角度旋转一周的旋转次数；
若Nm＞Nf，则剔除干扰数据点后剩余的采样数据中的数据点个数大于个数阈值；否则剔除干扰数据点后剩余的采样数据中的数据点个数小于或等于个数阈值。
作为优选，所述步骤S2.3和步骤S2.4，比例系数k和比例系数p的预设，包括：
(1)初始化粒子群，设置种群规模、粒子的初始速度和位置、最大迭代次数以及目标搜索空间维度，设置初始惯性权重、学习因子；
(2)通过每个粒子的位置计算每个粒子对应的适应度值，对于每个粒子，用粒子的适应度值和个体极值比较，若粒子的适应度值优于个体极值，则用粒子的适应度值替换个体极值；对于每个粒子，用粒子的适应度值和全局极值比较，若粒子的适应度值优于全局极值，则用粒子的适应度值替换全局极值；
(3)更新粒子的速度和位置：在H维空间中，由I个粒子构成的搜索集合中，第q个粒子在H维空间中的位置表示为xq＝(xq1，xq2，…，xqH)，第q个粒子的飞行速度表示为vq＝(vq1，vq2，…，vqH)，第q个粒子经过的最好位置记为pq＝(pq1，pq2，…，pqH)，所有粒子经过的最好位置记为pg＝(pg1，pg2，…，pgH)，q＝1，2，…，I，每一代粒子根据如下公式更新自己的速度和位置：
vqd＝wvqd+c1r1(pqd-xqd)+c2r2(pgd-xqd)
xqd＝xqd+vqd
其中，xqd是第q个粒子的第d维位置，vqd是第q个粒子的第d维速度，pqd是第q个粒子在第d维搜索到最优位置，pgd是整个粒子群在第d维搜索到的最优位置，w是惯性权重，c1和c2为学习因子，r1和r2为[0，1]范围内的均匀随机数，d＝1，2，3，…，H；
(4)若达到最大迭代次数，则结束迭代并输出全局极值pg，将最优的全局极值pg映射至比例系数k、p，获得最优的k、p值进行预设。
作为优选，所述步骤S3，利用最小二乘法，根据各检测点对应的拟合圆的圆心坐标计算拟合基准轴的直线方程，包括：
令拟合基准轴过点P0(x0，y0，z0)，其方向向量为(e，s，t)，则拟合基准轴的直线方程为：
Figure BDA0002651271490000051
对直线方程进行等价转换，得到：
Figure BDA0002651271490000052
Figure BDA0002651271490000053
其中，
Figure BDA0002651271490000054
令残差为
Figure BDA0002651271490000055
则
Figure BDA0002651271490000056
i＝1，2，3，…，n，可表示为：
Figure BDA0002651271490000057
Figure BDA0002651271490000058
其中，Xi，Yi，Zi为检测点i对应的拟合圆的圆心坐标，以
Figure BDA0002651271490000059
取最小值为目标，采用最小二乘法求解拟合基准轴的k1，k2，b1，b2如下：
设Q1(k1，b1)＝∑(Xi-k1Zi-b1)2；
Q2(k2，b2)＝∑(Yi-k2Zi-b2)2；
式中，Q1(k1，b1)和Q2(k2，b2)的偏导数满足最小化的条件为：
Figure BDA00026512714900000510
Figure BDA00026512714900000511
根据Q1(k1，b1)和Q2(k2，b2)的偏导数满足最小化的条件可得：
∑2(Xi-k1Zi-b1)*(-Zi)＝0
∑2(Xi-k1Zi-b1)*(-1)＝0
∑2(Yi-k2Zi-b2)*(-zi)＝0
∑2(Yi-k2Zi-b2)*(-1)＝0
因此求解得到：
Figure BDA0002651271490000061
Figure BDA0002651271490000062
则解得拟合基准轴的直线方程为：
Figure BDA0002651271490000063
Figure BDA0002651271490000064
完成计算拟合基准轴的直线方程。
作为优选，所述步骤S4，计算各拟合圆的圆心坐标到拟合基准轴的径向距离，包括：
计算各拟合圆的圆心坐标到拟合基准轴的径向距离为：
Figure BDA0002651271490000065
其中，Hi表示检测点i对应的拟合圆的圆心坐标(Xi，Yi，Zi)到拟合基准轴的径向距离，并且i＝1，2，3，…，n。
本申请提供的基于最小二乘法的轴类工件最大弯曲点的检测方法，基于最小二乘法求解轴类工件的最大弯曲点，提高求解效率，并且在求解拟合圆的圆心坐标时，对采样数据进行过滤，剔除干扰数据点后，求解得到更加准确的拟合圆的圆心坐标，以有效提高最大弯曲点的确定精度。
附图说明
图1为本申请基于最小二乘法的轴类工件最大弯曲点的检测方法的流程图；
图2为本申请基于最小二乘法拟合圆的模型；
图3为本申请基于最小二乘法拟合基准轴的模型。
具体实施方式
下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。
除非另有定义，本文所使用的所有的技术和科学术语与属于本申请的技术领域的技术人员通常理解的含义相同。本文中在本申请的说明书中所使用的术语只是为了描述具体的实施例的目的，不是在于限制本申请。
其中一个实施例中，提供一种基于最小二乘法的轴类工件最大弯曲点的检测方法，用于对待校直轴杆类工件进行拟合，能够更准确地找到最大弯曲点位置，方便校直。
如图1所示，本实施例中基于最小二乘法的轴类工件最大弯曲点的检测方法，包括：
步骤S1、沿轴类工件的轴向预设多个检测点，获取轴类工件周向旋转一周在各检测点处的采样数据。
在自动校直领域中，采集的数据主要是通过工作台多个检测点位置处的位移传感器检测的当前AD值，根据AD值进行最大弯曲点的判断。本实施例为了提高最大弯曲点的判断准确性，在工作台上布置N路位移传感器采集轴类工件的数据，即沿轴类工件的轴向预设N个检测点，则轴类工件的轴向截面数n为N，以Zi表示第i个检测点，即检测点i，其中i＝1,2,3,…,n。N的取值根据轴类工件的长度或者弯曲程度进行设置，例如可以是5、10等。
令轴类工件以固定角度旋转，获取每次旋转后各检测点处的采样数据，直至旋转M次后轴类工件旋转一周，得到M次旋转获取的采样数据，所述采样数据包括AD值和角度值。一般设置每次旋转的固定角度为1.2度，则轴类工件旋转一周大约是旋转M＝300次。
以每次旋转后每个检测点处获取的AD值和角度值作为一对数据，将所有采样数据保存至数组中，数组中的每个数据表示为(Zi,Rij,θij)，i＝1,2,3,…,n，j＝1,2,3,…,M，Rij表示在第j次旋转后检测点i处获取的AD值，θij表示在第j次旋转后检测点i处获取的角度值。
步骤S2、利用最小二乘法，根据获取的所述采样数据计算每个检测点对应的拟合圆的圆心坐标。
步骤S2.1、取一个检测点对应的采样数据，初始化圆拟合的次数。本实施例中每次初始化圆拟合的次数为0次。
步骤S2.2、根据该检测点对应的最新的采样数据利用最小二乘法计算拟合圆的圆心坐标，并累加一次圆拟合的次数。
对于一个检测点而言，轴类工件旋转一周采集的数据为沿轴类工件的周向依次采集的数据，因此所拟合的圆相当于轴类工件的一个横截面。
现有技术中存在利用一个检测点的数据进行最小二乘法拟合得到拟合圆的圆心坐标，但拟合的效果并不是很理想，因此如图2所示，本实施例提供一种优选的拟合方式如下：
一个检测点对应的采样数据有M对，取检测点i对应的采样数据为(Rij，θij)，i＝1，2，3，…，n，j＝1，2，3，…，M，则可计算得到检测点i对应的坐标Pi(Rijsinθij，Rijcosθij)，记为Pi(xij，yij)。
令检测点i对应的拟合圆的圆心坐标为Oi(Xi，Yi，Zi)，拟合圆的半径为R，则拟合圆的圆方程可表示为：
(xij-Xi)2+(yij-Yi)2＝R2
令残差为ei，并且ei可表示为：
ei＝(xij-Xi)2+(yij-Yi)2-R2
以函数∑ei 2取最小值为目标，采用最小二乘法求解拟合圆的Xi，Yi，如下：
设F(Xi，Yi，R)＝min∑ei 2＝min∑[(xij-Xi)2+(yij-Yi)2-R2]2。
则得到F(Xi，Yi，R)＝min∑[xij 2+yij 2-2Xixij-2Yiyij+Xi 2+Yi 2-R2]2。
令Q(a，b，c)＝∑(xij 2+yij 2+axij+byij+c)2，其中，a＝-2Xi，b＝-2Yi，c＝Xi 2+Yi 2-R2。
式中Q(a，b，c)的偏导数满足最小化的条件为：
Figure BDA0002651271490000081
其中，
Figure BDA0002651271490000082
表示偏导数。
根据Q(a，b，c)的偏导数满足最小化的条件可得：
∑2(xij 2+yij 2+axij+byij+c)Xi＝0
∑2(xij 2+yij 2+axij+byij+c)Yi＝0
Σ2(xij 2+yij 2+axij+byij+c)＝0
令A＝N∑xij 2-∑xij∑xij。
B＝N∑xijyij-∑xij∑yij。
C＝N∑xij 3+N∑xijyij 2-∑xij∑(xij 2+yij 2)。
D＝N∑yij 2-∑yij∑yij。
E＝N∑yij 3+N∑xij 2yij-∑yij∑(xij 2+yij 2)。
式中N表示检测点的个数，∑表示求和公式，例如∑yij表示以i和j进行循环对yij求和。
可得到
Figure BDA0002651271490000091
由此可以解得
Figure BDA0002651271490000092
解得检测点i对应的拟合圆的圆心坐标为
Figure BDA0002651271490000093
需要说明的是，以上为本实施例提供的一种优选的拟合方式，圆的拟合度高，准确率较高。在其他实施例或其他实施场景中，根据需要可以采用现有的最小二乘法进行拟合。
步骤S2.3、根据本次拟合圆的圆心坐标，剔除采样数据中的干扰数据点。
由于传感器的采集数据容易受环境影响，因此要进一步提高拟合圆的准确度，需要剔除受环境影响的采样数据。
在剔除数据时可采用现有的方法实现，例如狄克逊准则、罗曼诺夫斯基准则等。但本实施例为了提高干扰数据点的剔除率，同时避免有效数据点被剔除，提供一种优选的剔除方法如下：
对于一个检测点的采样数据而言，一次性剔除所有干扰数据点的可能性较小，因此本实施例根据圆拟合的次数进行动态剔除，以提高剔除的准确性。
令当前圆拟合的次数为m次，并且该检测点对应的最新的采样数据中每个数据点到本次拟合计算得到的拟合圆的圆心坐标的距离为dij如下：
Figure BDA0002651271490000094
取dij中最远的距离dmmax，设置本次剔除数据的阈值为dm，并且dm＝kdmmax，其中k为预设的比例系数。根据每次剔除时每个数据点与圆心坐标的距离，设置适用于本次剔除的阈值，从而可有效提高干扰数据点的剔除率，同时避免有效数据点被剔除。
剔除该检测点对应的最新的采样数据中dij＞dm的数据点，保留dij≤dm的数据点，并统计剔除干扰数据点后剩余的采样数据中的数据点个数Nm。
步骤S2.4、判断剔除干扰数据点后剩余的采样数据中的数据点个数，若数据点个数大于个数阈值，则以剔除干扰数据点后的采样数据作为最新的采样数据继续执行步骤S2.2；否则保存本次拟合圆的圆心坐标作为该检测点对应的拟合圆的圆心坐标，并执行步骤S2.5。
在判断剔除干扰数据点后剩余的采样数据中的数据点个数时，本实施例采用的方式如下：
令个数阈值为Nf，并且Nf＝pM，其中p为预设的比例系数，M为轴类工件以固定角度旋转一周的旋转次数。
若Nm＞Nf，则剔除干扰数据点后剩余的采样数据中的数据点个数大于个数阈值；否则剔除干扰数据点后剩余的采样数据中的数据点个数小于或等于个数阈值。
由于本实施例中在剔除干扰数据点时，比例系数k和比例系数p的影响较大，因此需要预设合适的k、p值，以保证数据剔除效果。
预设合适的k、p值时，可进行实验确定最佳的k、p值，但该方法局限性较大，无法全面反映是否为最佳的取值。因此本实施例采用粒子群算法对k、p值进行优化，获得最优的k、p值。具体的本实施例中k、p值预设过程如下：
(1)初始化粒子群，设置种群规模、粒子的初始速度和位置、最大迭代次数以及目标搜索空间维度，设置初始惯性权重、学习因子。
(2)通过每个粒子的位置计算每个粒子对应的适应度值，对于每个粒子，用粒子的适应度值和个体极值比较，若粒子的适应度值优于个体极值，则用粒子的适应度值替换个体极值；对于每个粒子，用粒子的适应度值和全局极值比较，若粒子的适应度值优于全局极值，则用粒子的适应度值替换全局极值。
(3)更新粒子的速度和位置：在H维空间中，由I个粒子构成的搜索集合中，第q个粒子在H维空间中的位置表示为xq＝(xq1，xq2，…，xqH)，第q个粒子的飞行速度表示为vq＝(vq1，vq2，…，vqH)，第q个粒子经过的最好位置记为pq＝(pq1，pq2，…，pqH)，所有粒子经过的最好位置记为pg＝(pg1，pg2，…，pgH)，q＝1，2，…，I，每一代粒子根据如下公式更新自己的速度和位置：
vqd＝wvqd+c1r1(pqd-xqd)+c2r2(pgd-xqd)
xqd＝xqd+vqd
其中，xqd是第q个粒子的第d维位置，vqd是第q个粒子的第d维速度，pqd是第q个粒子在第d维搜索到最优位置，pgd是整个粒子群在第d维搜索到的最优位置，w是惯性权重，c1和c2为学习因子，r1和r2为[0，1]范围内的均匀随机数，d＝1，2，3，…，H。
(4)若达到最大迭代次数，则结束迭代并输出全局极值pg，将最优的全局极值pg映射至比例系数k、p，获得最优的k、p值进行预设。本实施例中得到最优的比例系数k、p为k＝0.93、p＝0.89。需要说明的是，本实施例中关于粒子群算法未具体说明的部分，为现有技术中常规粒子群算法的执行逻辑，这里不再进行赘述。
步骤S2.5、判断是否完成所有检测点对应拟合圆的圆心坐标计算，若完成则执行步骤S3；否则继续执行步骤S2.1。
步骤S3、利用最小二乘法，根据各检测点对应的拟合圆的圆心坐标计算拟合基准轴的直线方程。
现有技术中存在利用多个坐标进行最小二乘法拟合得到直线方程的方法，但拟合的效果并不是很理想，因此如图3所示，本实施例提供一种优选的拟合方式如下：
令拟合基准轴过点P0(x0，y0，z0)，其方向向量为(e，s，t)，则拟合基准轴的直线方程为：
Figure BDA0002651271490000111
对直线方程进行等价转换，得到：
Figure BDA0002651271490000112
Figure BDA0002651271490000113
其中，
Figure BDA0002651271490000114
令残差为
Figure BDA0002651271490000115
则
Figure BDA0002651271490000116
i＝1，2，3，…，n，可表示为：
Figure BDA0002651271490000117
Figure BDA0002651271490000118
其中，Xi，Yi，Zi为检测点i对应的拟合圆的圆心坐标，以
Figure BDA0002651271490000119
取最小值为目标，采用最小二乘法求解拟合基准轴的k1，k2，b1，b2如下：
设Q1(k1，b1)＝∑(Xi-k1Zi-b1)2。
Q2(k2，b2)＝∑(Yi-k2Zi-b2)2。
式中，Q1(k1，b1)和Q2(k2，b2)的偏导数满足最小化的条件为：
Figure BDA0002651271490000121
Figure BDA0002651271490000122
根据Q1(k1，b1)和Q2(k2，b2)的偏导数满足最小化的条件可得：
∑2(Xi-k1Zi-b1)*(-Zi)＝0
∑2(Xi-k1Zi-b1)*(-1)＝0
∑2(Yi-k2Zi-b2)*(-zi)＝0
∑2(Yi-k2Zi-b2)*(-1)＝0
因此求解得到：
Figure BDA0002651271490000123
Figure BDA0002651271490000124
则解得拟合基准轴的直线方程为：
Figure BDA0002651271490000125
Figure BDA0002651271490000126
步骤S4、计算各拟合圆的圆心坐标到拟合基准轴的径向距离，取径向距离最远的圆心坐标对应的检测点作为轴类工件的最大弯曲点。
根据点到直线的距离计算原理可得，计算各拟合圆的圆心坐标到拟合基准轴的径向距离为：
Figure BDA0002651271490000127
其中，Hi表示检测点i对应的拟合圆的圆心坐标(Xi，Yi，Zi)到拟合基准轴的径向距离，并且i＝1，2，3，…，n。
求得个检测点对应的径向距离后，取径向距离最远的圆心坐标对应的检测点作为轴类工件的最大弯曲点。并且该检测点对应的圆心坐标到拟合基准轴的径向距离即为最大弯曲点的最大弯曲形变量。
由于每个检测点的位移传感器采集的AD值和角度值一一对应，因此找到最大弯曲点的最大弯曲形变量后，可以相应地确定最大弯曲点的方向角，取径向距离最远的圆心坐标对应的检测点中所有采样数据中的最大角度值作为最大弯曲点的方向角。
本实施例基于最小二乘法求解轴类工件的最大弯曲点，提高求解效率，并且在求解拟合圆的圆心坐标时，对采样数据进行过滤，剔除干扰数据点后，求解得到更加准确的拟合圆的圆心坐标，以有效提高最大弯曲点的确定精度。
应该理解的是，虽然流程图中的各个步骤按照箭头的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，图中的至少一部分步骤可以包括多个子步骤或者多个阶段，这些子步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些子步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤的子步骤或者阶段的至少一部分轮流或者交替地执行。
以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
以上所述实施例仅表达了本申请的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本申请构思的前提下，还可以做出若干变形和改进，这些都属于本申请的保护范围。因此，本申请专利的保护范围应以所附权利要求为准。

Claims (7)

1.一种基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述基于最小二乘法的轴类工件最大弯曲点的检测方法，包括：
步骤S1、沿轴类工件的轴向预设多个检测点，获取轴类工件周向旋转一周在各检测点处的采样数据，包括：
沿轴类工件的轴向预设N个检测点，则轴类工件的轴向截面数n为N，以Zi表示第i个检测点，即检测点i，其中i＝1，2，3，…，n；
令轴类工件以固定角度旋转，获取每次旋转后各检测点处的采样数据，直至旋转M次后轴类工件旋转一周，得到M次旋转获取的采样数据，所述采样数据包括AD值和角度值；
以每次旋转后每个检测点处获取的AD值和角度值作为一对数据，将所有采样数据保存至数组中，数组中的每个数据表示为(Zi，Rij，θij)，i＝1，2，3，…，n，j＝1，2，3，…，M，Rij表示在第j次旋转后检测点i处获取的AD值，θij表示在第j次旋转后检测点i处获取的角度值；
步骤S2、利用最小二乘法，根据获取的所述采样数据计算每个检测点对应的拟合圆的圆心坐标；
步骤S2.1、取一个检测点对应的采样数据，初始化圆拟合的次数；
步骤S2.2、根据该检测点对应的最新的采样数据利用最小二乘法计算拟合圆的圆心坐标，并累加一次圆拟合的次数；
步骤S2.3、根据本次拟合圆的圆心坐标，剔除采样数据中的干扰数据点；
步骤S2.4、判断剔除干扰数据点后剩余的采样数据中的数据点个数，若数据点个数大于个数阈值，则以剔除干扰数据点后的采样数据作为最新的采样数据继续执行步骤S2.2；否则保存本次拟合圆的圆心坐标作为该检测点对应的拟合圆的圆心坐标，并执行步骤S2.5；
步骤S2.5、判断是否完成所有检测点对应拟合圆的圆心坐标计算，若完成则执行步骤S3；否则继续执行步骤S2.1；
步骤S3、利用最小二乘法，根据各检测点对应的拟合圆的圆心坐标计算拟合基准轴的直线方程；
步骤S4、计算各拟合圆的圆心坐标到拟合基准轴的径向距离，取径向距离最远的圆心坐标对应的检测点作为轴类工件的最大弯曲点。
2.如权利要求1所述的基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述步骤S2.2中，根据该检测点对应的最新的采样数据利用最小二乘法计算拟合圆的圆心坐标，包括：
一个检测点对应的采样数据有M对，取检测点i对应的采样数据为(Rij，θij)，i＝1，2，3，…，n，j＝1，2，3，…，M，则可计算得到检测点i对应的坐标Pi(Rijsinθij，Rijcosθij)，记为Pi(xij，yij)；
令检测点i对应的拟合圆的圆心坐标为Oi(Xi，Yi，Zi)，拟合圆的半径为R，则拟合圆的圆方程可表示为：
(xij-Xi)2+(yij-Yi)2＝R2
令残差为ei，并且ei可表示为：
ei＝(xij-Xi)2+(yij-Yi)2-R2
以函数∑ei 2取最小值为目标，采用最小二乘法求解拟合圆的Xi，Yi，如下：
设F(Xi，Yi，R)＝min∑ei 2＝min∑[(xij-Xi)2+(yij-Yi)2-R2]2；
则得到F(Xi，Yi，R)＝min∑[xij 2+yij 2-2Xixij-2Yiyij+Xi 2+Yi 2-R2]2；
令Q(a，b，c)＝∑(xij 2+yij 2+axij+byij+c)2，其中，a＝-2Xi，b＝-2Yi，c＝Xi 2+Yi 2-R2；
式中Q(a，b，c)的偏导数满足最小化的条件为：
Figure FDA0003605115270000021
其中，
Figure FDA0003605115270000022
表示偏导数；
根据Q(a，b，c)的偏导数满足最小化的条件可得：
∑2(xij 2+yij 2+axij+byij+c)Xi＝0
∑2(xij 2+yij2+axij+byij+c)Yi＝0
∑2(xij 2+yij 2+axij+byij+c)＝0
令A＝N∑xij 2-∑xij∑xij；
B＝N∑xijyij-∑xij∑yij；
C＝N∑xij 3+N∑xijyij 2-∑xij∑(xij 2+yij 2)；
D＝N∑yij 2-∑yij∑yij；
E＝N∑yij 3+N∑xij 2yij-∑yij∑(xij 2+yij 2)；
可得到
Figure FDA0003605115270000023
由此可以解得
Figure FDA0003605115270000031
解得检测点i对应的拟合圆的圆心坐标为
Figure FDA0003605115270000032
3.如权利要求2所述的基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述步骤S2.3，根据本次拟合圆的圆心坐标，剔除采样数据中的干扰数据点，包括：
令当前圆拟合的次数为m次，并且该检测点对应的最新的采样数据中每个数据点到本次拟合计算得到的拟合圆的圆心坐标的距离为dij如下：
Figure FDA0003605115270000033
取dij中最远的距离dmmax，设置本次剔除数据的阈值为dm，并且dm＝kdmmax，其中k为预设的比例系数；
剔除该检测点对应的最新的采样数据中dij＞dm的数据点，保留dij≤dm的数据点，并统计剔除干扰数据点后剩余的采样数据中的数据点个数Nm。
4.如权利要求3所述的基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述步骤S2.4，判断剔除干扰数据点后剩余的采样数据中的数据点个数，包括：
令个数阈值为Nf，并且Nf＝pM，其中p为预设的比例系数，M为轴类工件以固定角度旋转一周的旋转次数；
若Nm＞Nf，则剔除干扰数据点后剩余的采样数据中的数据点个数大于个数阈值；否则剔除干扰数据点后剩余的采样数据中的数据点个数小于或等于个数阈值。
5.如权利要求4所述的基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述步骤S2.3和步骤S2.4，比例系数k和比例系数p的预设，包括：
(1)初始化粒子群，设置种群规模、粒子的初始速度和位置、最大迭代次数以及目标搜索空间维度，设置初始惯性权重、学习因子；
(2)通过每个粒子的位置计算每个粒子对应的适应度值，对于每个粒子，用粒子的适应度值和个体极值比较，若粒子的适应度值优于个体极值，则用粒子的适应度值替换个体极值；对于每个粒子，用粒子的适应度值和全局极值比较，若粒子的适应度值优于全局极值，则用粒子的适应度值替换全局极值；
(3)更新粒子的速度和位置：在H维空间中，由I个粒子构成的搜索集合中，第q个粒子在H维空间中的位置表示为xq＝(xq1，xq2，...，xqH)，第q个粒子的飞行速度表示为vq＝(vq1，vq2，...，vqH)，第q个粒子经过的最好位置记为pq＝(pq1，pq2，...，pqH)，所有粒子经过的最好位置记为pg＝(pg1，pg2，...，pgH)，q＝1，2，...，I，每一代粒子根据如下公式更新自己的速度和位置：
vqd＝wvqd+c1r1(pqd-xqd)+c2r2(pgd-xqd)
xqd＝xqd+vqd
其中，xqd是第q个粒子的第d维位置，vqd是第q个粒子的第d维速度，pqd是第q个粒子在第d维搜索到最优位置，pgd是整个粒子群在第d维搜索到的最优位置，w是惯性权重，c1和c2为学习因子，r1和r2为[0，1]范围内的均匀随机数，d＝1，2，3，...，H；
(4)若达到最大迭代次数，则结束迭代并输出全局极值pg，将最优的全局极值pg映射至比例系数k、p，获得最优的k、p值进行预设。
6.如权利要求2所述的基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述步骤S3，利用最小二乘法，根据各检测点对应的拟合圆的圆心坐标计算拟合基准轴的直线方程，包括：
令拟合基准轴过点P0(x0，y0，z0)，其方向向量为(e，s，t)，则拟合基准轴的直线方程为：
Figure FDA0003605115270000041
对直线方程进行等价转换，得到：
Figure FDA0003605115270000042
Figure FDA0003605115270000043
其中，
Figure FDA0003605115270000044
令残差为
Figure FDA0003605115270000045
则
Figure FDA0003605115270000046
可表示为：
Figure FDA0003605115270000047
Figure FDA0003605115270000048
其中，Xi，Yi，Zi为检测点i对应的拟合圆的圆心坐标，以
Figure FDA0003605115270000049
取最小值为目标，采用最小二乘法求解拟合基准轴的k1，k2，b1，b2如下：
设Q1(k1，b1)＝∑(Xi-k1Zi-b1)2；
Q2(k2，b2)＝∑(Yi-k2Zi-b2)2；
式中，Q1(k1，b1)和Q2(k2，b2)的偏导数满足最小化的条件为：
Figure FDA0003605115270000051
Figure FDA0003605115270000052
根据Q1(k1，b1)和Q2(k2，b2)的偏导数满足最小化的条件可得：
∑2(Xi-k1Zi-b1)*(-Zi)＝0
∑2(Xi-k1Zi-b1)*(-1)＝0
∑2(Yi-k2Zi-b2)*(-zi)＝0
∑2(Yi-k2Zi-b2)*(-1)＝0
因此求解得到：
Figure FDA0003605115270000053
Figure FDA0003605115270000054
则解得拟合基准轴的直线方程为：
Figure FDA0003605115270000055
Figure FDA0003605115270000056
完成计算拟合基准轴的直线方程。
7.如权利要求6所述的基于最小二乘法的轴类工件最大弯曲点的检测方法，其特征在于，所述步骤S4，计算各拟合圆的圆心坐标到拟合基准轴的径向距离，包括：
计算各拟合圆的圆心坐标到拟合基准轴的径向距离为：
Figure FDA0003605115270000057
其中，Hi表示检测点i对应的拟合圆的圆心坐标(Xi，Yi，Zi)到拟合基准轴的径向距离，并且i＝1，2，3，…，n。
CN202010871515.4A 2020-08-26 2020-08-26 一种基于最小二乘法的轴类工件最大弯曲点的检测方法 Active CN112090994B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010871515.4A CN112090994B (zh) 2020-08-26 2020-08-26 一种基于最小二乘法的轴类工件最大弯曲点的检测方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010871515.4A CN112090994B (zh) 2020-08-26 2020-08-26 一种基于最小二乘法的轴类工件最大弯曲点的检测方法

Publications (2)

Publication Number Publication Date
CN112090994A CN112090994A (zh) 2020-12-18
CN112090994B true CN112090994B (zh) 2022-08-09

Family

ID=73756721

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010871515.4A Active CN112090994B (zh) 2020-08-26 2020-08-26 一种基于最小二乘法的轴类工件最大弯曲点的检测方法

Country Status (1)

Country Link
CN (1) CN112090994B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113843308B (zh) * 2021-08-25 2022-07-12 北京科技大学 一种金属型材压力矫直策略方法
CN115752344B (zh) * 2022-11-15 2023-09-05 上海羿弓精密科技有限公司 一种rv减速器曲柄轴相位夹角的检测方法

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1237326C (zh) * 2003-03-08 2006-01-18 东风汽车公司 曲轴弯曲变形的测量方法
CN101196394A (zh) * 2007-09-27 2008-06-11 上海大学 小段圆弧圆度的优化最小二乘评价方法
US7813599B2 (en) * 2009-02-23 2010-10-12 The United States Of America As Represented By The Administrator Of The National Aeronautics And Space Administration Method and apparatus for shape and end position determination using an optical fiber
CN110064680B (zh) * 2019-04-12 2021-07-16 太原科技大学 一种棒材大弯曲变形快速测量方法
CN110849291B (zh) * 2019-11-26 2021-05-11 二重（德阳）重型装备有限公司 大型弯管弯曲半径检测方法

Also Published As

Publication number Publication date
CN112090994A (zh) 2020-12-18

Similar Documents

Publication Publication Date Title
CN112090994B (zh) 一种基于最小二乘法的轴类工件最大弯曲点的检测方法
CN109975793B (zh) 一种激光二维测距传感器的标定方法
CN110929347A (zh) 一种基于梯度提升树模型的热连轧带钢凸度预测方法
CN111667156B (zh) 一种多点生产的卷烟物理质量一致性评价的方法
CN103488561B (zh) 一种在线升级主样本模型的kNN故障检测方法
CN111582123B (zh) 一种基于信标识别与视觉slam的agv定位方法
CN110068274B (zh) 涂胶传感器胶条检测示教方法
CN113688534B (zh) 一种基于多特征融合模型寻找最优铣削参数的研究方法
CN117148784A (zh) 一种多轴多通道数控系统运行故障分析方法
CN112518425B (zh) 基于多源样本迁移强化学习的智能加工刀具磨损预测方法
CN110908361B (zh) 一种基于在线孪生支持向量回归机的发酵过程软测量方法
CN111626359A (zh) 数据融合方法、装置、控制终端和船舶
CN114061957A (zh) 一种柴油机主轴承的健康评估方法
CN108446452B (zh) 一种混流泵叶轮鲁棒优化设计方法
CN111546337B (zh) 基于自由曲面的工业机器人全覆盖路径生成方法及系统
CN108638060B (zh) 多自由度机器人参数标定中冗余参数分析剔除方法
CN106078509B (zh) 用于外圆切入磨削工件圆度误差实时计算方法
CN113468794A (zh) 一种小型密闭空间温湿度预测及反向优化方法
CN111599016A (zh) 点云误差计算方法
CN114487284A (zh) 一种测量空气中重金属浓度的方法及系统
CN113399344B (zh) 高压喷射清洗机工艺参数优化方法和计算装置
CN109614758A (zh) 具有空间相关性的圆形形状误差的监控方法
CN111967489B (zh) 一种基于质量数据流形特征的制造过程异常监测方法
CN108133261B (zh) 基于bp神经网络的汽车胶管芯轴质量评价方法
CN113537454A (zh) 基于ADw-CLPSO径向基情感神经网络的出水氨氮软测量方法

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