基于振动的齿形磨床简易预知维修系统设计

基于振动的齿形磨床简易预知维修系统设计,基于,振动,齿形,磨床,简易,预知,维修,系统,设计 摘要齿形磨床加工齿形时，其机械性能状态会随着工作时间的延长而下降。齿形磨床振动所引起的零件精度问题是零件安装、装配的关键，它对整个机器的使用寿命和产品质量存在至关重要的影响。零件精度误差会导致生产的零件使用寿命下降甚至报废。预知维修系统广泛应用于各行各业中，其也成为当前国际上研究的重点课题。本文设计的预知维修系统由测试系统、计算机软件两部分组成。首先应用测试系统进行振动信号提取特征参数。然后通过 VB 编写的程序，进行频域分析比较和时域细化分析。然后运用最小二乘法原理建立预知维修预测数学模型。最后，依照预知维修模型拟合的曲线上的警戒线进行预知维修状态监控。振动的齿形磨床简易预知维修系统可以对机械设备在振动故障潜伏期时就对机械进行预知维修，在故障潜伏期内，将其消除，保证机械设备的正常工作状态。达到齿形磨床的预知维修要求。通过该毕业设计可以实现振动的齿形磨床简易预知维修的功能。关键词：最小二乘法；预知维修；VBIABSTRACTWith the extension of working time,the mechanical properties of the gear grinding machine will decrease . The precision of part precision caused by the vibration of tooth grinding machine is the key to the parts installation and assembly, which has a crucial influence on the service life and product quality of the whole machine. The precision error of the parts can lead to the decrease of the service life of the parts and even the scrap. The predictive maintenance system is widely used in all walks of life, and it has become the focus of international research.The predictive maintenance system designed in this paper consists of two parts: test system and computer software.Firstly, the test system is applied to extract characteristic parameters of vibration signal. Thenthe VB program accomplishes the frequency domain analysis comparison and time domain refinement analysis. Then the mathematical model of predictive maintenance is established by using the least square method. Finally, the detection and maintenance condition monitoring is conducted according to the warning line on the curve fitting according to the predicted maintenance model. Vibration of the tooth profile grinding machine easy to predict maintenance system to the vibration fault of mechanical equipment in the incubation period is to predict the mechanical maintenance, the failure latent period, its elimination, guarantee the normal work of the mechanical equipment state. The anticipation maintenance requirement of gear grinding machine is achieved.Through this graduation design can realize the vibration of the toothshape grinder simple predictive maintenance function.Keywords : leastsquaremethod ;predictivemaintenance ;VBII目录摘要IABSTRACTII1. 绪论11.1 研究背景11.2 选题的理论意义和应用价值11.3 预知维修的国内外发展趋势22. 齿形磨床简易预知维修概述32.1 齿形磨床的振动32.1.1 齿形磨床的振动机理32.1.2 齿形磨床的振动原因32.1.3 齿形磨床的振动测试位置42.2 预知维修原理研究52.3 预知维修系统组成53. 振动信号特征参数提取方法研究63.1 计算机测试系统的组成63.2 特征参数的选择63.3 振动信号传感器的选择73.4 信号调理电路73.5 离散傅里叶变换83.6 提取数字化的振动信号特征参数124. 预知维修预测数学模型建立和计算机实现134.1 最小二乘法数学模型研究134.1.1 最小二乘法原理：134.1.2 一元线性回归模型134.1.2 多元线性回归模型154.2 预知维修预测数学模型计算机实现185. 系统软件编程225.1 VB 编制预知维修软件22结论与展望41参考文献42附录 1：外文翻译43附录 2：外文原文52致谢59I基于振动的齿形磨床简易预知维修系统1. 绪论1.1 研究背景随着国家经济的不断快速发展，在信息数字化的智能化时代，各行各业的技术水平不断提高，传统的机械维修方式太过于落后，已经不能满足机械行业的生产效率和生产效益。机械维修经历了“事后维修”、“定期维修”两个阶段。“事后维修”是一种消极和不得已的维修策略，在工业生产发展高度现代化和复杂化的今天，已经不能适应市场激烈的竞争。“定期维修”虽然对企业的设备维修保养和专业化大生产曾起过不可忽视的积极作用，但随着现代工业的发展、生产规模的逐年扩大以及对降低生产成本和提高企业经济效益要求的提高，特别是对于企业全面推行“精益生产”方式和实现“敏捷制造”系统的要求，定期维修体制似乎已走到了尽头。为了应对新时代下的大国竞争， 急需一种快速、提前预知机械设备故障的维修方式，在设备未发生故障之前，就做好准备。预知维修能使设备或部件的使用寿命得到最充分的利用，而故障可减少到最低限度。机器设备在实际工作中不可避免地会出现由于零部件损坏而影响生产的情况，提前预知故障并采取维修措施可以有效避免这一问题，为机器的应用降低了一部分隐患。预知维修系统应用越来越广泛，尤其是“工业2.0”和“中国制造2025”概念提出之后，更成为智能化、现代化工厂企业提高经济生产效率，获得最大额度生产效益的利剑。增强国家竞争力的国之重器，因此此研究具有了强烈的现实指导意义。1.2 选题的理论意义和应用价值随着社会的快速发展，科技的飞速进步，机电一体化的不断完善。在计算机智能控制下，预知维修有了新的发展前景，成为了人们关注的焦点。机械振动引起的零件精度问题是零件安装、装配的关键，乃至于影响整个机器的寿命和产品质量问题。振动是影响齿形磨床生产力和表面质量的一个最常见的局限操作。预知维修是一种动态维修制度，也称为视情维修。它可以在不拆卸、不破坏设备的条件下，对设备的工作性能、状态变化及故障趋势进行定量的描述。预知维修能使设备或部件的使用寿命得到最充分的利用，而故障可减少到最低限度。利用设备振动分析技术，在设备运行中使用振动分析仪器进行信号测量，通过对设备监测点进行间断或连续的信号数据采集，并对振动信号的诊断分析，有效判定设备状态以及故障的部位和原因，并预测设备状态发展趋势，实现对设备预知维修，使故障隐患消除在萌芽状态，达到设备状态监测与诊断的目的。且简易预知维修具有投资少、见效快、技术水平要求较低和易于推广等特点。经过实践证9明，预知维修了很好的效果。预知维修的最大意义在于：利用计算机技术实现简易预知维修，可大大提高设备有效利用率，降低维修成本。减少振动故障造成的损失。 （1）避免不必要的维修；（2） 节省返工 （修）的费用； （3）减少故障造成的生产损失；（4）减少备件库存； （5） 提高加工过程的效率； （6）提高产品质量； （7）延长设备系统的使用寿命； （8） 提高生产能力； （9）节约维修费用； （10）获得好的经济效益。目前，预知维修系统作为一项新的方式，正在逐渐被认识并在被广大的工业企业中应用。研究表明，在很多的工厂企业已经取得了很高的效益，很好的应用价值。预知维修满足了机器设备故障的检测与维修，提高了社会生产能力，促进了社会的进步与发展。1.3 预知维修的国内外发展趋势预知维修一直是国内外相关领域的重点关注对象。被广泛用于机械、兵工武器、煤矿、石油产业等。大约几十年前，国外已经开始发展了预知维修，并被运用了实际工程。到目前，它已经涵盖了军事、天文、航母、飞机等重要领域。应用于武器、导弹、天文望远镜、飞机发动机等各个方面。与国外相比，我国预知维修比国外先进国家制造行业发展晚了几十年，但在不断地追赶努力下，已取得了重要成果。我国自主研发的维修系统还不能完全满足生产预知维修需要的精度与准确度。但随着传感器技术、信号处理技术、计算机技术的发展，更为设备技术诊断学提供了良好的平台。通过采集到的信号的幅值、频率、冲击个数等，通过多种判断方式综合给出故障的严酷程度。已能很好的做到故障诊断和预知维修。近年来，机器人应用越来越广泛，尤其是“工业2.0”和“中国制造2025”概念提出之后，更是掀起了机器人应用的大潮，用机器人替代人工完成繁重的体力劳动甚至是脑力劳动是当前工业发展的必然趋势。机器人在实际工作中不可避免地会出现由于零部件损坏而影响生产的情况，提前预知机器人故障并采取维修措施可以有效避免这一问 题，为机器人的应用降低了一部分隐患。工业无线技术为机器人故障诊断提供了灵活、方便的传输方式。未来基于工厂自动化WIA-FA无线网络的机器人故障诊断和预知系统， 为工业企业和机器人生产厂家提前预知机器人故障提供了一种可执行的思路。预知维修系统将在智能化、便利化、效益化等方面大显身手。因此，预知维修系统有一个广阔的发展前景。2. 齿形磨床简易预知维修概述2.1 齿形磨床的振动2.1.1 齿形磨床的振动机理机械设备中任何一个运动部件与之相关的零件出现故障，必然破坏机械运动的平稳性，在传递力的参与下，这种力和运动的非平稳性现象表现为振动。磨床加工中的振动主要有强迫振动和自激振动。强迫振动是在外力作用下的振动； 自激振动是系统自身所产生的激励而引起的振动，也称颤振。任何振动工程问题可用下图1.1概括说明。输入是指作用在系统上的激励或干扰， 输出也称为响应。对于机械振动而言，激励大多为力，而常用的响应物理量一般为位移、速度和加速度。输入输出图2.1 振动问题简化系统模型2.1.2 齿形磨床的振动原因磨床震颤是磨削加工中最重要的误差源之一，对工件最终几何精度有直接的影响。自激振动是齿形磨床的主要振动来源，也称颤振。对于磨削加工，主要体现为再生型颤振。再生型震颤是指砂轮或工件表面上一转留下的振纹对瞬时磨削力的影响，进而影响磨削震颤。(1) 砂轮颤振砂轮由砂轮轴承支撑，经砂轮主轴、皮带轮与砂轮电机连接，其中任一部分的振动都有可能导致砂轮的颤振。齿形磨床磨削时，引起砂轮强迫振动的主要原因是:砂轮不平衡由于砂轮不平衡产生的离心力使砂轮产生振动， 振动的频率等于砂轮回转频率， 振动幅值由砂轮不平衡量的大小决定。该振动直接作用在工件上，产生等于砂轮转频的振纹。砂轮驱动电机的转子不平衡该组件的不平衡产生等于砂轮电机转频的振动，该振动通过传动链传给砂轮而作用在工件上，产生频率等于电机转频的振纹。砂轮电机传动皮带长短不一致对于三角皮带传动，由于长短不一致而引起皮带的拍打，形成固定频率的强迫振动源，对于磨削用量很小的精磨和光磨作业，此类振动的影响不可忽略。（2）工件震颤齿轮被精确夹夹具固定，砂轮进行磨齿时，由于砂轮与齿轮之间磨削力的作用， 工件会产生微小变形，引起工件震颤齿形磨床在磨削时，砂轮与齿轮间的接触存在跳动。由于齿轮被精确夹夹具固定、刚性好，齿轮出现颤振的可能性非常小，因此主要是砂轮颤振。2.1.3 齿形磨床的振动测试位置在磨齿过程中,砂轮轴的轴向振动和砂轮的修整质量对被磨齿轮的齿廓偏差影响较大,由于砂轮轴高速旋转,因此不能直接用接触式传感器进行测量,故通过测量砂轮轴前端盖的振动来近似地反映砂轮轴的轴向振动。此外 ,砂轮修整器的振动也会影响到砂轮的修整质量,从而引起被磨齿轮的齿廓偏差,所以将砂轮修整器的导轨和砂轮轴的前端盖作为振动的测试位置如图2。图2.2 测试系统连接通过加速度传感器即时采样得到砂轮修整器导轨(通道 CH1)和砂轮轴前端盖(通道CH2)处的加速度信号,提取振动信号的特征参数。2.2 预知维修原理研究图2.3为机械设备性能状态与运行时间之间的关系曲线,从图中可以看到机械设备在使用过程中,其性能状态随着使用时间的推移而逐步下降。通过研究发现,大多数故障发生前都会有一些明显征兆, 即所谓的潜在故障, 如图中的A 点,如所谓的潜在故障, 如图中的 A 点,图2.3 机械设备工作性能变化曲线这时故障依然可以识别, 但此时设备依然能够“正常运转 ”。经过一定的潜伏期后形成功能性故障 , 如图中的 B点 ,同时设备也失去既定功能。为了避免功 能故障停产造成的损失 ,很显然有必要寻求给出设备到达 A、 B点状态时有关参数的标准 ,力求在故障潜伏期内通过机械设备性能的变化规律，较为正确地预测其运行状态10。因此，我们需要对提取的振动信号的特征参数设置维修标准。先设定一基准特征参数值，依据基准特征参数值设置机械设备正常工作的上下界限，即是警戒线。通过把提取的振动信号的特征参数与标准参数进行比较，来判断齿形磨床能否正常工作。当提取的数据反映在预知维修模型拟合的曲线上时，如果其介于上限和下限之间，则齿形磨床仍能够正常运行，但已在警戒线内，即在故障潜伏期内。我们此时应当适当维修。如果其高于上限，则齿形磨床不能进行正常工作。若低于下限，则表示齿形磨床能正常工作。2.3 预知维修系统组成此系统由测试系统、计算机软件两部分组成。其各部分的功能介绍如下：测试部分包括：传感器、中间变换装置、显示记录装置等。传感器是能感受被测对象并按照一定规律将其转换成可用以输出信号的器件。中间变换装置进行信号的调理与转换并对信号进行处理等，将来自传感器的信号转换成更适合于传输和处理的形式，且对信号进行各种运算和分析。显示记录装置是测试部分的最终环节，信号显示是将来自信号处理的信号记录，以观察者易于观察和分析的形式来显示或储存。计算机软件主要是由VB编写，面向对象程序设计作为应用程序软件开发的主流技术 , 不仅具有封装性、 继承性、 多态性、 动态链接性等多个优点 , 更重要的是它能够简单地实现友好的人机交互界面。在实际软件的编制过程中 , 本统开发平台为VB6.0, 主要考虑其界面实现和管理方便灵活 , 且数据库维护简易。对于较为复杂的信号分析处理 , 先用 C +实现 , 然后以动态链接库的方式链接到主程序中。3. 振动信号特征参数提取方法研究3.1 计算机测试系统的组成图3.1 测试系统的组成图框如图3.1所示，被测振动信号由传感器测得，将振动信号转换成电信号，对传感器输出的模拟信号进行放大、滤波、调制与解调。 然后 将模拟信号进行A/D转换为数字信号，并实现数据采集功能的计算机将其数据传到计算机，完成振动信号提取。3.2 特征参数的选择描述信号的时域特征参数通常有峰值、峰峰值、均值、绝对均值、方差、均方差和有效值等。信号的时域分析：直接观测或记录的信号一般随时间变化的物理量，这种以时间为独立变量，用信号的幅值随时间变化的函数或图形来描述信号的方法称为时域描述。信号的时域波形时是时域描述的一种重要形式。峰值 xp 用于描述信号 x(t )在时域中出现的最大瞬时幅值，是指时域波形上与零线的最大偏离值，即xp =x(t )max（3-1）峰值在实际应用中有它的价值。对信号的峰值应该有足够的估计，以便确定测试系统的动态测量范围，不会产生削波的现象，从而能真实地反映被测信号的最大值。均值mx 是指信号在时间间隔T内幅值对时间的平均值，即m = 1 T x(t )dtxT 0（3-2）均值表示了信号变化的中心趋势，也称固定分量或直流分量。经过综合分析考虑，本次试验选取峰值与均值作为信号的研究特征参数3.3 振动信号传感器的选择压电式加速度传感器图3.2 压电式加速度传感器的结构图1- 螺栓2- 压电元件3- 预压弹簧4- 外壳5- 质量块6- 基座如图3.2所示，它主要由压电元件、质量块、预压弹簧、螺栓、基座及外壳等组成。整个部件装在外壳内并由螺栓予以固定。当加速度传感器和被测物一起受到冲击振动 时，压电元件受质量块惯性力的作用，根据牛顿第二定律，此惯性力是加速度的函数，即F = ma式中： F 为质量块产生的惯性力； m 质量块的质量； a 为加速度。由于压电元件的压电效应，压电元件所产生的电荷与力 F 成正比，亦即与加速度传感器所感受到的加速度成正比。因此，测得加速度传感器输出的电荷，便可知加速度的大小。本实验通过压电式加速度传感器测量振动信号的情况。3.4 信号调理电路信号的放大：在机械信号测量中，传感器的输出信号电压是很微弱的，一般不能直接用于显示、记录或数模转换，需要进行信号放大。因此需采用电荷放大器对压电式加速度传感器获得的信号进行放大。信号的滤波：在进行信号分析与处理时，经常会遇到有用信号叠加上无用信号的问题，这些噪声有些是与信号同时产生的，有些是在信号传输时混入的，噪声有时会大于有用信号从而淹没有用信号。所以，从原始信号中消除或减弱干扰噪声是信号处理的一个重要问题。故而需要用滤波器对信号进行滤波提取有用信号。信号的调制与解调：由于被测参数，经过传感器交换以后，多为低频缓变的微弱信号，因此在传输中被测信号有可能存在被淹没的问题，为了解决微弱缓变信号的放大问题，信号处理技术中需要采用一种对信号进行调制，即先将微弱的缓变信号加载到高频交流信号中去，然后利用交流放大器进行放大，最后再从放大器的输出信号中取出放大了的缓变信号。这种信号传输中的变换过程称为调制与解调。信号的数字化处理：在用计算机分析处理连续测量的机械信号时，都需要对连续测量的动态信号经过数字化处理，转换成离散的数字序列。数字化的过程主要包括对模拟信号离散采样和幅值量化和编码。首先，采用保持器把预处理后的模拟信号按选定的采样间隔采样为离散序列，此时的信号变为时间离散而幅值连续的采样信号。然后，量化编码装置将每一个采样信号的幅值转换为数字码，最终把采样信号变为数字序列。3.5 离散傅里叶变换对连续时间信号进行离散傅里叶变换。图3.3 时域采样采样过程如图3.3，在数字信号处理前，把一个模拟信号转化为计算机可以接受的信号，这一转变通过对模拟信号的采样来完成，信号的采样是由模-数转换电路（A/D）来实施的。把连续时间信号转换离散数字信号的过程称为模-数（A/D）转换过程。模-数（A/D）转换过程包括采样、量化、编码，工作原理如下图3.4。图3.4 模-数（A/D）转换过程采样是在模/数转换过程中以一定时间间隔对连续时间信号进行取值的过程。从理论上讲，采样的实质是用间隔为Ts 的周期单位脉冲序列函数 g(t) 与时域信号 x(t) 相乘。g(t) 可写为20g(t) =+d(t - nTs ),n=-n = 0,1,2,（3-9）由d函数的筛选性质知+- x(t)d(t - nTs )dt = x(nTs ),n = 0,1,2,（3-10）式表明，经过时域采样后，各采样点的信号幅值为 x(nTs ) 函数 g(t) 称为采用函数。根据傅里叶变换性质，相乘结果在频域上表现为两者卷积，表示为F x(t)g(t) = X ( f ) * G( f ) = 1+nX ( f -)Ts n=-Ts（3-11）卷积结果是把 X ( f ) 平移到各脉冲所在频率位置上，即G( f ) 各个脉冲发生位置上重新构图， X ( f ) * G( f ) 在频域上则表现为周期连续函数。在采用间隔合适时，一个周期上的频谱与原信号 x(t) 的频谱 X ( f ) 将保持一致。时域截断图3.5 时域截断过程在信号测试中，所采集的信号 x(t) 通常经历时间很长，以至于长度超出了计算机所能处理的范围，所以采样后数据量就很大，而计算机只能处理其中一部分数据，这在理论上就相当于用一个宽度为T的窗函数w(t) 去乘 x(t) g(t) ，其时域内函数x(t) g(t) w(t) = 1T+x(nTs )d(t - nTs )在频域内则三者的卷积，表示为s n=-（3-12） = 1 +- n *pF x(t)g(t)w(t)Ts n=-X ( f) T sin c(TsTf )（3-13）矩形窗函数的频谱是抽象函数，当函数的宽度T越宽抽样函数中间的尖峰区域就越窄，形状越尖，能量愈集中。若宽度T 趋近无穷大，此时， x(t) g(t) w(t) 的乘积在频域中的结果将与图3.3所示一样。如果宽度T 足够宽，则时域截断前后的频谱图形基本一致。由于时域截断，所以图3.5中频谱图出现了皱纹波。频域采样图3.6 频率采样过程图3.6说明频域采样过程。经过时域采样和截断处理后，信号 x(t) 变成了“有限长”的离散信号。但从频域来看，这一有限离散信号的频谱仍是一连续信号。实际上，计算机处理后，所输出的频率一定是离散的，这表明计算机对连续频谱进行了频域采样。时域函数*= 1 +d - * + d-x(t)g (t)w(t)d (t)Tx(nTs ) (tnTs )(nnT ) s n=- n=-（3-14）频域采样理论上就是对3.5中的连续频谱乘上频率采样G( f ) 的过程，这样，其连续频谱在频域内也离散化了。频率离散化的时域函数表示为 = 1 +- n *p 1 + d- nF x(t) g(t)w(t)TX ( f) T sin c(TTf ) T( f)T s n=-ss n=-（3-15）这个过程是离散傅里叶变换的演变过程。3.6 提取数字化的振动信号特征参数在计算机上，应用VB编写的预知维修程序软件，进行提取数字化的振动信号特征参数峰值 xp 与均值mx ，将其数据传送至计算机进行最小二乘法模型建立函数关系。4. 预知维修预测数学模型建立和计算机实现4.1 最小二乘法数学模型研究4.1.1 最小二乘法原理：最小二乘法通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。4.1.2 一元线性回归模型给定数据点 pi (xi , yi ) ，其中i=1,2,3，m。数据点的两个变量 y 与 x 之间存在某种关系，可假设如下结构：y = b0 + b1 x +e（4-1）其中，b0 、b1 是两个未知参数，e为其他随机因素对 y 的影响， x 是非随机可精确 观察的，e是均值为零的随机变量，是不可观察的。一般的，由（4-1）确定的模型为一元线性回归模型，记为 y = b0 + b1 x +eEe= 0, De= s2固定的未知参数b0 、b1 称为回归系数，自变量 x 称为回归变量。（4-1）两边同时取期望得： Y = b0 + b1 x ，称 y 对 x 的回归直线方程。（4-2）运用最小二乘法估计模型中的未知参数b0 、b1 。假设有n 组独立观测值：(x1 , y1 ), (x2 , y2 ),.(xn , yn ) ， 则由（4-2）yi = b0 + b1 xi +ei , i = 1,2, nEe= 0, De = s2且e，e，e 相互独立ii记12nnn2(4-3)Q = (b ,b) = e2 = ( y- b - bx )01ii=1ii=10 1 i1称Q(b0 ,b1 )为偏离真实直线的偏差平方和。最小二乘法就是选择b0 和b 的估计b0 ,b1 ，使得Q(,)= min Q(b ,b)b0 b1b0 ,b101为此，将上式分别对b0 、b1 求偏导数，得 Q= -2( y - b - bx ) bi01 ini=1 0 Q = -2nx ( y - b - bx )bii01 i1i=1(4-4)令上两式为0，并用b0 、b1 取代b0 ， b1 ，得 n ( yi - b0 - b1 xi ) = 0i=1 i=1n于是有 xi ( yi - b0 - b1 xi ) = 0(4-5) nnn b0 + b1 xi = yi ni=1 ni=1nb x + b x2 ) = x y 0ii=11 ii=1i ii=1(4-6)此方程组称为正规方程。由正规方程解得 b0 = y - b1 x xy - x yb1 = 2x2 - x1 n 1 n21 n21 n(4-7)其中， x = n xi ， y = n yi ， x = n x i ， xy = n xi yii=1i=1i=1i=1用此最小二乘法可求出用b0 、b1 即得方程x = y (x - x)y = b0 + b1+ b1(4-8)显然，b1 为拟合直线的斜率，b0 是拟合直线在 x = x0 的截距。n个点 pi (xi , yi )(i = 1,2, n)的几何重心(x, y)落在拟合直线上。本题研究两个变量（x,y）之间的相互关系时，得到一系列成对的数据（x1,y1.x2,y2. xm,ym）；将这些数据描绘在x -y直角坐标系中，发现这些点在某条直线附近，可以令这条线方程如（式4-9）。其中： a0 、b1 是任意实数。y = a0 x + b1（4-9）试验提峰值和均值数据，运用最小二乘法原理拟合曲线，由计算机提取的多组数据：峰值 xp 与均值mx按偏差平方和最小的原则选取拟合曲线，并且采取一次函数方程为拟合曲线，完成维修预测模型的建立。4.1.2多元线性回归模型一般地，影响试验指标的因素往往不止一个，即有多个因素，可假设如下结构：y = b0 + b1 x1 + + bk xk +e(4-10)其中， y 可精确观察随机变量，即因变量。 x1 , x2 , xk为非随机精确观察变量，称自变量。b0 ,b1 ,b2 ,bk 为k +1 个未知参数，e是随机变量，一般假设 Ee= 0, De= s2 0 。 为了估计未知参数b0 ,b1 ,b2 ,bk 及s2 ，我们对 y 与 x1 , x2 , xk 同时做n 次观察得n 组观 察值， t = 1, n(n k +1)，它们满足yt = b0 + b1 xt1 + + bk xtk +et ,t = 1, n(4-11)其中，e1 ,en 互不相关且均是与e同分布随机变量。用矩阵表示上式，令y1 1x11x12K x1k b1 e1 M 1xxK x b e Y = X = 21222k b= 2 e= 2 M MMMK M M M y n ，1xn1xn 2K xnk ，bn ，en In于是（4-10）式变为Y = Xb+e(4-12)其中 X 为已知的n (k +1) 矩阵，称为回归设计矩阵或资料矩阵。b为k +1 维未知列向量，e是满足E(e) = 02COV (e,e) = s In(4-13)tn的n 维随机列向量，其中列向量s2 是未知参数， D(e) = s2 ， t = 1, n ， I 为n 阶单位矩阵，即对随机误差e1 ,en 作无偏、等方差与互不相关的假定。Y 是n 维观察列向量。一般称Y = Xb+enE(e) = 0, COV (e,e) = s2 I为高斯-马尔科夫线性模型（k元线性回归模型）。对（4-10）式取期望(4-14)称为回归平面方程。y = b0 + b1 x1 + + bk xk(4-15)运用最小二乘法求b0 ,b1 ,b2 ,bk 的估计量。作离差平方和n2Q = ( yi - b0 - b1 xi1 -L- bk xik )i=1(4-16)选择b0 ,b1 ,b2 ,bk ，使Q 达到最小，即Q = min 。根据微积分学中求最小值得方法，只 需求方程组= Qb0 0 Q = 0 b1M Qkb = 0b=的解。解方程组得到的不是b0 ,b1 ,b2 ,bk 的真值，而是估计值b0 ,L,bk ，故可将方程 组化简为Ab= B 或其中(X T X ) X TY 1x11x12K x1k y1 b0 1xxK x M X = 21222k Y = b= b1 MMMK M M M 1xn1xn 2K xnk ， yn ，bk 解得估计值b = A-1B = (X T X )-1 (X TY )bbx = xi, i = 1,2,L, k.这样计算得到的i 代入回归平面方程。i 称为经验回归系数。令 i多项式回归模型，变为多元线性回归模型。设 x, y 的回归模型为y = b + bx + bx2 + + b xk +e011k(4-17)其中b0 ,b1 ,b2 ,bk 是未知参数，e服从正态分布 N (0,m2 )。因而Ey = b + bx + bx2 + + b xk011k(4-18)本题研究两个变量（x,y）之间的相互关系时，得到一系列成对的数据（x1,y1.x2,y2. xm,ym）；将这些数据描绘在x -y直角坐标系中，若发现这些点在某条直线或曲线附近，可以令这条线方程如（式4-18）。y = a x2 + b x + c其中： a0 、b1 、c2i012是任意实数。（4-19）为建立这条线方程就要确定a0 、b1 、c2 ，应用最小二乘法原理，将实测值Yi与利用y = a x2 + b x + c计算值Yj（ i012 ）（式1-1）的偏差（Yi-Yj）的平方和最小为“优化判据”。试验提峰值和均值数据，运用最小二乘法原理拟合曲线，由计算机提取的多组数据：峰值 xp 与均值mx按偏差平方和最小的原则选取拟合曲线，并且采取二次函数方程为拟合曲线，完成维修预测模型的建立。4.2 预知维修预测数学模型计算机实现本实验通过VB编写的程序，提取振动信号的特征参数峰值和有效值，进行频域分析比较和时域细化分析。频谱参数提取频谱比较时域细化运用最小二乘法原理对提取的提取振动信号的特征参数进行曲线一次拟合和二次拟合。得到预知维修预测数学模型。一次曲线拟合二次曲线拟合29预测出齿形磨床的故障潜伏期。依据画出的振动信号正常工作的上限和下限来判断齿形磨床能否正常工作。当提取的数据反映在预知维修模型拟合的曲线上时，如果其介于上限和下限之间，则齿形磨床能够正常运行，但状态在故障潜伏期。如果其高于上限， 则齿形磨床不能进行正常工作。若低于下限则表示齿形磨床能正常工作。5. 系统软件编程5.1 VB 编制预知维修软件各程序代码：/定义全局变量/Public Type MeaStruct不加密 MEA 文件结构EquipPath9As String * 9mea文件路径EquipCode11 As String * 11设备编号EquipType11 As String * 11设备型号EquipName5As String * 5设备文件名标识EquipPoint3 As String * 3设备测点编号MeasureNum3 As String * 3设备测点序号RunSpeed6As String * 6测量时设备转速Chanel2As String * 2测量通道数Anyfreq6As String * 6分析频率DataLength3 As String * 3采样长度MeasureDate9 As String * 9测量日期WinType6As String * 6加窗窗口类型FftLength6As String * 6FFT长度Overratio5As String * 5重叠系数AlarmOne6As String * 6 警 戒 线 1 AlarmTwo6As String * 6 警 戒 线 2 DirectOne2As String * 2A通道方向DirectTwo2As String * 2B通道方向UnitOne6As String * 6 A通道单位UnitTwo6 As String * 6 B通道单位SensorOne2 As String * 2 A通道传感器类型SensorTwo2 As String * 2 B通道传感器类型SensiOne6 As String * 6 A通道灵敏度SensiTwo6 As String * 6 B通道灵敏度GainOne6As String*6A通道增益GainTwo6As String*6B通道增益End TypePublic Type MeaEncrypt加密后 MEA 文件结构EquipPath(8)As Doublemea文件路径EquipCode(10) As Double设备编号EquipType(10) As Double设备型号EquipName(4)As Double设备文件名标识EquipPoint(2) As Double设备测点编号MeasureNum(2) As Double设备测点序号RunSpeed(5)As Double测量时设备转速Chanel(1)As Double测量通道数Anyfreq(5)As Double分析频率DataLength(2) As Double采样长度MeasureDate(8) As Double测量日期WinType(5)As Double加窗窗口类型FftLength(5)As DoubleFFT长度Overratio(4)As Double重叠系数AlarmOne(5)As Double 警 戒 线 1 AlarmTwo(5)As Double 警 戒 线 2 DirectOne(1)As DoubleA通道方向DirectTwo(1)As DoubleB通道方向UnitOne(5)As Double A 通 道 单 位UnitTwo(5) As Double B 通 道 单 位SensorOne(1) As Double A通道传感器类型SensorTwo(1) As Double B通道传感器类型SensiOne(5) As Double A通道灵敏度SensiTwo(5) As Double B通道灵敏度GainOne(5) As Double A 通 道 增 益GainTwo(5) As Double B 通 道 增 益End Type / 画坐标Public Sub Drawzbiao(zbiao As PictureBox) Dim x1, y1, x2, y2 As SingleDim i As Integer DrawStyle = vbDash zbiao.ScaleMode = 6x1 = 0x2 = zbiao.ScaleWidth - 10 y1 = zbiao.ScaleHeight - 10 y2 = 0zbiao.Scale (x1, y1)-(x2, y2) For i = 1 To 15 Step 1zbiao.DrawStyle = 2zbiao.ScaleMode = 0zbiao.Line (10 * i, 0)-(10 * i, -zbiao.ScaleHeight), RGB(192, 192, 192)zbiao.Line (0, 7 * i)-(zbiao.ScaleWidth, 7 * i), RGB(192, 192, 192) NextEnd Sub画时域波形Public Sub Drawtime(bxing As PictureBox, m_meastru As MeaStruct, mydata() As Double)bxing.Cls Drawzbiao bxing Dim n As LongDim rameter As Double Dim cof As DoubleDim i As Integer Dim a As DoubleDim maxvalue As Double bxing.DrawStyle = vbslod bxing.ScaleMode = 6n = 1024 * CInt(m_meastru.DataLength3)rameter = 5000# / 32768# / Val(m_meastru.GainOne6) /Val(m_meastru.SensiOne6)For i = 1 To n - 1mydata(i) = mydata(i) * rameter Next icof = max2(mydata(), m_meastru)a = 0.8 * bxing.ScaleHeight / 2 / cof For i = 1 To n - 1bxing.Line (bxing.ScaleLeft + (i - 1) * bxing.ScaleWidth / n, a * mydata(i - 1) + bxing.ScaleHeight / 2) _-(bxing.ScaleLeft + i * bxing.ScaleWidth / n, a * mydata(i)+ bxing.ScaleHeight / 2), RGB(0, 255, 0)Next画基准线bxing.Line (bxing.ScaleLeft, bxing.ScaleTop + bxing.ScaleHeight /2)-(bxing.ScaleLeft + bxing.ScaleWidth, bxing.ScaleTop + bxing.ScaleHeight / 2), RGB(192, 0, 192)End Sub画频域波形Public Sub Drawspec(Drawtimewave As PictureBox, m_meastru As MeaStruct, datasp() As Double)Dim a As Double Dim cof As Double Dim n As LongDim m As Integer Drawtimewave.Cls Drawtimewave.DrawStyle = vbslod Drawtimewave.ScaleMode = 6cof = max1(datasp(), m_meastru)a = 0.8 * Drawtimewave.ScaleHeight / cof n = 512For i = 1 To n - 1Drawtimewave.Line (Drawtimewave.ScaleLeft + (i - 1) * Drawtimewave.ScaleWidth / n, -a * datasp(i - 1) + Drawtimewave.ScaleHeight) _-(Drawtimewave.ScaleLeft + i * Drawtimewave.ScaleWidth / n, -a * datasp(i) + Drawtimewave.ScaleHeight), RGB(0, 255, 0)NextLabel19 = smax(datasp() End SubPublic Sub Drawqushi(Drawtimewave As PictureBox, m_meastru As MeaStruct, jspecdata() As Double, cspecdata() As Double)Dim a As Double Dim cof As Double Dim n As LongDim b As Double Dim m As Integer Drawtimewave.ClsDrawtimewave.DrawStyle = vbslod Drawtimewave.ScaleMode = 6cof1 = max1(jspecdata(), m_meastru) Val (Form2.Text5) *cof2 = max1(cspecdata(), m_meastru)If cof1 cof2 Then b = cof1Elseb = cof2 End Ifa = 0.98 * Drawtimewave.ScaleHeight / b n = 512