本发明涉及不确定性反问题技术领域，尤其涉及一种滑动轴承-转子系统不平衡量识别方法。
背景技术：
滑动轴承-转子系统材质的不均匀性、制造和安装过程引起的变形以及工作中的磨损等都会引起轴承转子系统的不平衡振动。利用不平衡量识别技术进行轴系动平衡来提高转子及其构成的产品质量，减小噪声和振动，提高轴承的试用寿命，以保证轴系运行的长期性和稳定性，是轴承转子系统经常使用的一种校准方式。但由于影响因素的多样性和复杂性，不平衡量识别结果存在一定的误差。对于滑动轴承-转子系统，转子几何特征、滑动轴承油膜特性系数和测量响应的随机性这些不确定性因素即使是在较小情况下，也很可能导致不平衡量识别结果产生较大的偏差。而工程实际中，工程师能根据经验和知识对不平衡量参数在尚未获取实验测量信息之前有一定的预先估计。如何利用这些先验信息减小不确定因素对待识别参数的影响成为不确定性反问题领域的研究热点。
贝叶斯理论对结构模型参数进行反求分析，同时考虑了已知参数概率密度样本信息和未知参数先验信息。但贝叶斯理论求解工程不确定性反问题常常涉及非常耗时的正问题计算，难以满足实际工程对计算效率的要求。同时，MCMC(Markov Chain Monte Carlo,马尔科夫链蒙特卡罗)法获取不平衡量边缘后验概率密度分布时要生成大量的抽样点，并要大量调用耗时的正问题计算，效率低的问题。
技术实现要素：
鉴于上述的分析，本发明旨在提供一种滑动轴承-转子系统不平衡量识别方法，用以解决贝叶斯理论和MCMC法在采样中要生成大量的抽样点，及求解正问题时耗时严重、效率低的问题。
本发明的目的主要是通过以下技术方案实现的：
提供了一种滑动轴承-转子系统不平衡量识别方法，包括以下步骤：
步骤S1、获取滑动轴承-转子系统不平衡量先验分布空间、初始样本及滑动轴承-转子系统的测试不平衡响应；
步骤S2，利用贝叶斯理论求解代价函数，并获得不平衡量的后验联合概率密度分布；
步骤S3，计算上述代价函数的最小值，当所述最小值大于收敛精度时，更新所述滑动轴承-转子系统不平衡量先验分布空间，并利用遗传智能采样技术获取新样本，进入步骤S2，否则，进入步骤S4；所述遗传智能采样技术是将旧样本遗传到新样本中；
步骤S4，采用MCMC方法识别得到滑动轴承-转子系统不平衡量的近似边缘后验概率密度分布，进而确定滑动轴承-转子系统不平衡量的均值、置信区间。
本发明有益效果如下：本发明基于遗传智能采样技术的滑动轴承-转子系统不平衡量识别方法，大大减小样本的计算次数，使采集样本和先验空间逐步向联合后验概率密度分布的高密度区域集中，提高采样效率，并遗传有效样本点，减少正问题调用次数。有效地对识别结果进行评价，减小不确定性因素对识别结果的影响。该方法在不牺牲计算精度的前提下能够获得较高的计算效率，能够准确而快速地识别不平衡量的均值和置信区间等信息。进而，为工程师进行转子动平衡提供辅助决策依据。
在上述方案的基础上，本发明还做了如下改进：
进一步，所述利用遗传智能采样技术获取新样本，为通过遗传拉丁超立方实验设计更新先验分布空间中的新样本，包括：
根据极大极小距离准则，对落入到下代先验分布空间的样本进行筛选；
采用模拟火优化求解器使新样本在各个设计变量上投影均匀且产生的新样本到遗传样本的距离最大。
采用上述进一步方案的有益效果是：通过遗传拉丁超立方实验设计(ILHD)不断更新先验空间和遗传信赖区中的样本，使采集样本和先验空间逐步向联合后验概率密度分布的高密度区域集中，提高采样效率。
进一步，所述求解代价函数，包括：
采用传递矩阵法将不平衡量作为输入变量通过滑动轴承-转子系统不平衡响应仿真分析，求取滑动轴承-转子系统不平衡响应；
将上述滑动轴承-转子系统测试不平衡响应与求取的不平衡响应之间的误差作为代价函数。
进一步，所述更新先验分布空间，为根据信赖域更新方法获取下代先验分布空间。
进一步，所述信赖域更新方法为：通过黄金分割法自适应调整不平衡量先验分布空间的上下边界XL和XR，进而获得下代采样点的先验分布空间。
进一步，所述获取不平衡量先验分布包括：
设定误差常数ε，给定滑动轴承-转子系统不平衡量先验均匀分布的初始上下边界和给定每一迭代步样本点数量Nsi。
进一步，所述不平衡量的后验联合概率密度分布p(X|q)为：
式中，X为滑动轴承-转子系统不平衡量；q为滑动轴承-转子系统的测试不平衡响应；a为常数；g(X)为代价函数；σ2为滑动轴承-转子系统随机噪声的方差。
进一步，代价函数计算公式为：
式中，X为滑动轴承-转子系统不平衡量，h为测点个数，q为滑动轴承-转子系统的测试不平衡响应，qj(X)为不平衡量作为输入变量通过滑动轴承-转子系统不平衡响应仿真分析获得的不平衡响应。
进一步，所述获得下代采样点的先验分布空间为：
式中，1-α为黄金分割点，α＝0.382；和为第i代采样先验分布空间的上、下边界；和为第i+1代采样先验分布空间的上、下边界。
进一步，在MCMC方法中，链长取105量级。
本发明中，上述各技术方案之间还可以相互组合，以实现更多的优选组合方案。本发明的其他特征和优点将在随后的说明书中阐述，并且，部分优点可从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过说明书、权利要求书以及附图中所特别指出的内容中来实现和获得。
附图说明
附图仅用于示出具体实施例的目的，而并不认为是对本发明的限制，在整个附图中，相同的参考符号表示相同的部件。
图1为本发明实施例1中滑动轴承-转子系统不平衡量识别方法流程图；
图2为本发明实施例2中转子实验台及测试方案；
图3为本发明实施例2中转子实验台参数模型；
图4为本发明实施例2中转子实验台传递矩阵法模型；
图5为本发明实施例2中转子实验台两测点x径向方向瞬态响应；
图6为本发明实施例2中转子实验台两测点y径向方向瞬态响应；
图7为本发明实施例2中样本点的遗传智能更新迭代过程(第一步)；
图8为本发明实施例2中样本点的遗传智能更新迭代过程(第二步)；
图9为本发明实施例2中样本点的遗传智能更新迭代过程(第三步)；
图10为本发明实施例2中抽样的二维随机游走图；
图11为本发明实施例2中抽样的一维随机游走图；
图12为本发明实施例3中两测点动平衡前后的X向位移时域响应比较图；
图13为本发明实施例3中两测点动平衡前后的Y向位移时域响应比较图。
具体实施方式
下面结合附图来具体描述本发明的优选实施例，其中，附图构成本申请一部分，并与本发明的实施例一起用于阐释本发明的原理，并非用于限定本发明的范围。
实施例1
为了在滑动轴承-转子动平衡应用时一方面利用由滑动轴承-转子系统服役环境变化、测试条件的局限和滑动轴承-转子系统信息中存在的不确定因素的样本信息，来减小不确定因素对不平衡量识别结果的影响。另一方面利用工程师根据经验和知识对不平衡量参数在尚未获取实验测量信息之前的预先估计信息来提高不平衡量的识别精度。本实施例采用能处理不确定因素的样本信息和不平衡量的先验信息的贝叶斯理论及MCMC法获得不平衡量的均值、置信区间。针对贝叶斯理论在轴承转子系统不平衡量识别中的计算效率问题，提出了一种基于遗传智能采样技术的滑动轴承-转子系统不平衡量识别方法，包括以下步骤：
步骤S1、获取滑动轴承-转子系统不平衡量先验分布空间、初始样本及滑动轴承-转子系统的测试不平衡响应；
步骤S2，利用贝叶斯理论求解代价函数，并获得不平衡量的后验联合概率密度分布；
步骤S3，计算上述代价函数的最小值，当所述最小值大于收敛精度时，更新所述滑动轴承-转子系统不平衡量先验分布空间，并利用遗传智能采样技术获取新样本，进入步骤S2，否则，进入步骤S4；所述遗传智能采样技术是将旧样本遗传到新样本中；
步骤S4，采用MCMC方法识别得到滑动轴承-转子系统不平衡量的近似边缘后验概率密度分布，进而确定滑动轴承-转子系统不平衡量的均值、置信区间。
与现有技术相比，本实施例提供的滑动轴承-转子系统不平衡量识别方法，大大减小了样本的计算次数，使采集样本和先验空间逐步向联合后验概率密度分布的高密度区域集中，提高计算效率。有效地对识别结果进行评价，减小不确定性因素对识别结果的影响。该方法在不牺牲计算精度的前提下能够获得较高的计算效率，能够准确而快速地识别不平衡量的均值和置信区间等信息。这些信息为工程师进行转子动平衡提供辅助决策依据。
转子不平衡量的识别是在转子动力学分析的基础上进行的，即根据转子动力学的基本分析方法建立转子的运动方程，进而由可测知的滑动轴承-转子振动响应(不平衡响应)，结合运动方程，建立包括未知系统参数的不平衡量识别方程。转子系统的角速度、阻尼和刚度矩阵已知，不平衡量(质量、偏心距、相位)为未知参数作为待识别参数，不平衡响应可以通过实验测量获得，通过不平衡响应来识别不平衡量。当未知参数或转子结构参数有不确定因素时，不平衡响应也是不确定的。这种情况下基于不平衡响应识别不平衡量为不确定性识别问题。
在贝叶斯理论中，构建未知参数联合概率密度获得代价函数时，要获得高精度的联合后验概率密度分布空间要进行大量的采样，一种有效的做法是能够使在先验分布空间产生的用于构建代价函数所需的有限样本点集中在呈现真实联合后验概率密度分布的高密度区域。实际工程中，通过给定的参数先验分布是相对真实区间较大的空间，因此，采样时要不断调整先验分布空间以便最大限度地集中有限的样本点反映真实后验空间高密度区域，进而保证联合后验概率分布的精度。本实施例基于遗传智能采样技术采样，通过信赖域模型管理方法探测非支配解区域，不断更新先验分布空间，来保证获得与真实解接近的区间。通过智能布点技术和样本遗传策略使各个信赖域上的样本均匀分布，通过遗传部分样本落入下代信赖域的点作为智能布点，减小正问题计算次数而提高信赖域上关键区域的精度，从而加快收敛速度。
具体来说，步骤S101，获取不平衡量先验分布；
设定误差常数(即收敛精度)ε，滑动轴承-转子系统不平衡量先验均匀分布的初始上下边界和给定每一迭代步样本点数量Nsi。
需要说明的是，给定初始样本点数量可以设定为相对较小，此时的样本点主要用于探测样本点后续加密的方向，同时也不易太小，较小的样本点将无法反应真实参数联合后验概率密度分布的主要特征，这会增加加密迭代的次数。
对于滑动轴承-转子系统不平衡量的先验分布用均匀分布来表示，其上下界由实践经验或专家知识来确定。
步骤S102，对样本进行采样；优选的，采用LHD方法获得滑动轴承-转子系统不平衡量的初始样本X1；在之后的每步迭代中，通过遗传拉丁超立方实验设计方法，获取滑动轴承-转子系统不平衡量的新样本。
步骤S103，测量得到滑动轴承-转子系统的测试不平衡响应(如：瞬时位移响应)；
在步骤S2中，求解代价函数；通过上述获取的滑动轴承-转子系统的测试不平衡响应，并计算求取不平衡量作为输入变量通过转子系统不平衡响应仿真分析获得的不平衡响应，得到代价函数，进而获得不平衡量的后验联合概率密度分布；
具体地，利用下述贝叶斯理论近似公式，求取滑动轴承-转子系统不平衡量的后验联合概率密度分布：
式中X为滑动轴承-转子系统不平衡量(质量m，偏心矩e，相位)；q为滑动轴承-转子系统的测试不平衡响应；p(X|q)为滑动轴承-转子系统不平衡量的后验联合概率密度分布；a为常数；g(X)为代价函数；h为测点个数；qj(X)为滑动轴承-转子系统不平衡量作为输入变量通过滑动轴承-转子系统不平衡响应仿真分析获得的不平衡响应(即不平衡响应计算值)。
测试不平衡响应的不确定性用独立的随机噪声表征，为均值为0，方差为σ2的正态分布。
需要强调的是，代价函数g(X)为滑动轴承-转子系统测试不平衡响应与其不平衡响应计算值之间的误差，即最小二乘意义下的误差，本实施例中用代价函数作为指示因子以调整不平衡量的先验空间以便其逐渐覆盖高密度后验空间。
其中，滑动轴承-转子系统的不平衡响应计算值是通过传递矩法，结合滑动轴承-转子运动方程求取的，用少量的有效样本构建高精度的代价函数，进而获得未知参数不平衡量的后验联合概率密度分布p(X|q)。
在步骤S3中，计算上述代价函数的最小值，并将此最小值与预先设定的收敛精度进行对比；当上述最小值大于收敛精度时，更新先验空间，并利用遗传智能采样技术获取新样本，进入步骤S2，利用新样本求取代价函数；否则，进入步骤S4；
由于先验分布不确定是否包含了其真实后验空间，有时只能反映出参数后验空间的局部信息，有时甚至相差较大，因此，需要调整先验分布空间以便最大限度地集中有限的样本点反应其真实后验空间高密度区域，进而保证联合后验概率密度分布的近似精度。具体地，
步骤S301，计算代价函数的最小值；并将此最小值与步骤S1中设定的收敛精度进行对比；
为第i步代价函数的最小值，获得此步中导致代价函数最小值的不平衡量向量由于此时误差最小，说明联合先验概率密度分布的高密度区域集中在样本点周围。将此最小值与收敛精度ε，如果ymin≤ε，则联合后验概率密度分布高密度区域获得，迭代结束，进入步骤S5。当大于收敛精度，则进入步骤S402。
步骤S302，根据信赖域更新方法确定下代信赖域区域(即下代的先验分布空间)；
当大于收敛精度，以黄金分割法自适应调整不平衡量先验均匀分布的上下边界和从而获得下代采样点的信赖域区域。
下一步与当前步的不平衡量先验空间关系为：
式中1-α为黄金分割点，α＝0.382；和为当前i步先验空间的上下边界；和为下一步(i+1步)先验空间的上下边界。
根据信赖域更新方法更新不平衡量先验均匀分布的上下边界和确保生成的有限样本点能够反映未知参数不平衡量联合后验概率密度分布的高密度区域。
步骤S303，在上述获取的新的先验空间中，基于遗传智能采样技术进行采样；
下代信赖域与当代的信赖域会有重合的区域。当代的样本点可能会落入下代信赖域里。本实施例将遗传的旧样本和遗传拉丁超立方实验设计(ILHD)产生的新样本组合起来作为下代代价函数构建立的样本点，这样可以大大减少需要做不平衡响应计算的总样本的个数，提高计算效率。
另外，如果将上一步得到的样本全部遗传给下一步，会产生一部分区域样本过于紧凑而不利于智能布点的样本均匀分布。为此，本实施例采用样本遗传策略在尽可能充分利用遗传的旧样本的原则下根据极大极小距离准则对拉丁超立方实验设计进行优化，筛选部分落入到下代信赖域的样本，同时，采用模拟退火优化求解器使新样本在各个设计变量上投影均匀且产生的新样本到遗传样本的距离最大，遗传样本和新产生新样本在下代信赖域区域保持空间均布性、投影均匀性。
其中，根据极大极小距离原则，计算样本集Si+1最小距离dmin(即相邻样本的最小距离)，如下式：
式中，Ns为下代信赖域样本总数；分别为第i+1代样本集中第j和第k个样本点。
落入下代信赖域的遗传样本到下代样本集Si+1的最大距离Dmax的求取公式如下所示：
式中，为第i代遗传样本中的第m个样本点，m＝1,2,…,Np,Np为遗传样本个数。为第i+1代样本集Si+1中的第n个样本点，n＝1,2,…,Ns。
将经过上述落入到新的信赖域中的遗传样本与利用优化后的拉丁超立方实验设计产生的新样本作为新一代的样本。
在步骤S4中，基于高精度的代价函数，利用上述求取的不平衡量的后验联合概率密度分布，采用MCMC方法识别获得不平衡量的近似边缘后验概率密度分布，进而确定不平衡量的均值、置信区间。
设置链长(抽样点数目)Nmax，优选的，链长Nmax取105量级。通过大量实验分析，链长选择此量级时MCMC法才可能得到后验分布较为准确的近似解。
具体地，根据未知参数不平衡量的后验联合概率密度分布，马尔科夫链从当前样本中获取候选样本；通过求解接受概率，判断候选样本点是否本接受；不断获取候选样本，直到满足马尔科夫链长条件。
实施例2
本实施例通过具体的实验，对实施例1中的方法进行验证，
首先，建立转子轴承系统模型
如图2所示为INV1612型转子实验台作为实施例1中方法检验模型，圆盘厚度为15mm，直径60mm处沿周向分布16个孔，可配置不同重量的螺钉来人为设置不平衡质量。转子的两端由滑动轴承支承，其油膜特性系数如表1。转速设计为2000r/min。
表1滑动轴承油膜刚度和阻尼系数
表2不确定变量的中点和区间
然后，依照实施例1中的方法，进行不平衡量识别；
为了避免传感器磁头发生磁场交叉耦合引起的误差，所以x，y传感器不要安装在同一平面内。布置两个测点，一个测点距离左轴承为73mm，用来测x方向径向位移，一个测点距离左轴承334mm，用来测y方向径向位移。通过涡流传感器测得两测点的位移瞬态响应如图5、图6所示。
假设不平衡量初始先验信息服从独立的均匀分布m～U(0，20)和误差常数取为1e-7，马尔科夫链长为50000，采用最优拉丁超立方实验设计(OLHD)获取11个初始步样本。遗传落在下一步信赖区域的样本点，采用ILHD实验设计遗传样本点并在信赖区域产生新的样本点。图7给出采用识别不平衡量质量和相位的具体过程。图7-图9显示用本方法通过3步样本点的遗传智能更新加密迭代过程，图7、8、9显示每一迭代步更新的先验空间中加密样本点的总体分布情况，用符号“*”、“□”和“☆”分别表示对应步中产生的样本，通过实施例1中遗传拉丁超立方实验设计更新先验分布空间中的新样本，第1步产生11个样本点，第2步产生12个新样本，遗传了11个样本，第3步产生16个新样本，遗传了23个样本。ILHD实验设计将遗传的旧样本和ILHD产生的新样本组合起来作为下代的样本点(本实施例中，经过刷选遗传样本，在前3步中，上一代的遗传样本全部遗传到了下一代中)，这样可以大大减少总样本的个数，达到提高效率的目的。第3步遗传了23个样本，减少了23次采样，这对样本获取耗时的工程实际问题来说意义重大。
图10给出不平衡量二维空间中，MCMC方法在未知参数不平衡量联合后验概率分布上50000次抽样的计算结果，图11给出了不平衡质量和相位角各自的马尔科夫链的收敛过程，表明其相应的识别结果有95％的可能性落在狭窄区间[1.2319,1.2322]和[46.0354,46.0358]中。其抽样结果的接受率为42.6％，这与Tarantola建议的接受概率范围为30％-50％的比较吻合。
实施例3
通过上述实施例识别出转子不平衡量后，可以直接对转子进行动平衡，以解决各种因素引起的转子不平衡。
为了验证不平衡量识别的准确性和精度，本实施例利用上面实施例2的转子实验装置进行不平衡度的识别，并对加了配重后的转子系统进行测试。配重的大小与识别出的不平衡量大小相同，相位角相反，即在原来的角度上加上180°或者减去180°，测试加了配重后的系统振动响应，并进行配重前后响应对比。从图12、13可以看出，测点在动平衡前后的振动幅度有了很大的下降，动平衡的效果比较好。
综上，本发明技术方案带来的有益效果：本发明利用了工程师不平衡量的先验信息，也考虑了转子系统的环境因素和测试噪声，符合实际工程条件，有工程实际应用价值。同时，本发明基于遗传智能采样技术减少贝叶斯理论中调用耗时的正问题计算次数，提高了计算效率，并使有限的样本集中反映真实后验空间高密度区域，进而保证了不平衡量识别的精度。
本领域技术人员可以理解，实现上述实施例方法的全部或部分流程，可以通过计算机程序指令相关的硬件来完成，所述的程序可存储于计算机可读存储介质中。其中，所述计算机可读存储介质为磁盘、光盘、只读存储记忆体或随机存储记忆体等。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。