Method for characterizing the behavior of a rotating shaft using equivalent ellipses

Abstract

[Method for characterizing the behavior of a rotating shaft by using equivalent ellipses] This method for characterizing the behavior of a rotating shaft comprises the following steps: - a) acquisition over a predetermined period of at least two input signals comprising measurements representative of the movement of the rotating shaft (step 10); - b) determination of two components of the input signals, the components being orthogonal to each other (step 20); - c) merging the two components into a bivariate signal (step 30); - d) time-frequency decomposition of the bivariate signal into spectral elements (step 40); - e) determination of an equivalent ellipse for each spectral element (step 50); and - f) extraction of at least one characterization indicator of the rotating shaft from parameters of the equivalent ellipse (step 80). Figure for abstract: Fig. 1

Description

Title of the invention: Method for characterizing the behavior of a rotating shaft using equivalent ellipses Technical field of the invention

[0001] The present invention relates to the characterization of a shaft rotating around an axis of rotation.

[0002] The present invention relates more particularly to the characterization of the gyroscopic effects and the deformations of a rotating shaft in order to facilitate its control, modeling, monitoring, diagnosis and maintenance. State of the prior art

[0003] A rotating system, such as a rotating shaft, perfectly balanced rotates around an axis of rotation so that said axis of rotation is stationary relative to its theoretical positioning, and therefore relative to the environment close to the rotating shaft.

[0004] However, defects, deformations, gyroscopic effects, or even effects due to the force of gravity can induce an offset of the axis of rotation relative to its theoretical position during the rotation of said rotating shaft. This is the case, for example, of deformations of the rotating shaft due to flexible natural modes, these modes deforming the shaft as a function of the rotation speed of said shaft.

[0005] When such defects are present, the axis of rotation is not stationary and vibrates. It moves away from its theoretical position along a particular trajectory. More precisely, each point of the axis of rotation can move away from its theoretical position, two points of the axis of rotation at two different ends of the shaft can have a different trajectory.

[0006] The properties of a rotating shaft are conventionally understood by a lateral analysis of parameters. The parameters are values ​​of radial quantities such as the position or the radial speed of the shaft, having a non-zero component in a plane orthogonal to the axis of rotation.

[0007] It is often necessary to characterize the changes in these quantities depending on whether the shaft is rotating or not. However, multiple causes can influence the rotational behavior of a rotating shaft and it is complex to characterize such behavior precisely. Summary of the invention

[0008] The present invention therefore aims to overcome the aforementioned drawbacks.

[0009] The present invention relates to a method for characterizing the behavior of a rotating shaft, the method comprising the following steps:

[0010] - a) acquisition over a predetermined period of at least two input signals, the input signals comprising measurements representative of the movement of the rotating shaft, the measurements, linked to a vector quantity, being non-collinear with each other and optionally coplanar with a plane intersecting the axis of rotation of the shaft;

[0011] - b) determination of two components of the input signals, the components being orthogonal to each other and optionally to the axis of rotation;

[0012] - c) merging the two components into a bivariate signal;

[0013] - d) time-frequency decomposition of the bivariate signal into spectral elements;

[0014] - e) determination of an equivalent ellipse for each spectral element; and

[0015] - f) extraction of at least one characterization indicator of the rotating shaft at from parameters of the equivalent ellipse.

[0016] The gyroscopic effects and the rotor deformations, in particular the flexible modes, are thus simpler to consider under a frequency decomposition of the input signals. Indeed, the rotating effects of these deformations are asynchronous, in other words their frequencies are not directly linked to the rotation frequency of the shaft, and are therefore more easily characterized using a frequency decomposition which can highlight a rotating orbit, also called an equivalent ellipse, for each frequency, each orbit rotating at a specific speed corresponding to the frequency considered.

[0017] The present invention therefore makes it possible to enhance the characterization of a rotating shaft, allowing applications of the control, identification, modeling, monitoring or diagnostic type.

[0018] Advantageously, the input signals comprise measurements of position, or speed, or acceleration, or force, or magnetic field, or magnetic induction of the rotating shaft.

[0019] In various embodiments, the time-frequency decomposition step is performed by Fourier transform, or by wavelet, or by filtering, or by quaternion Fourier transform.

[0020] In a particular embodiment, the step of determining an equivalent ellipse comprises the calculation of Euler parameters of said equivalent ellipse.

[0021] Advantageously, the method comprises a step of calculating the spectral density of the time-frequency decomposition, and a step of calculating Stokes parameters from said spectral density.

[0022] In one embodiment, the at least one characterization indicator comprises a polarization rate, and / or an indicator of circularity and / or amplitude and / or orientation and / or direction of rotation and / or sign and / or synchronization of an equivalent ellipse.

[0023] Advantageously, the method comprises a step of inserting at least one characterization indicator into a neural network to classify said indicator as a function of the rotation of the shaft and link it to a previously listed anomaly.

[0024] In a particular embodiment, the method comprises a step of displaying the at least one extracted characterization indicator to be readable by an operator.

[0025] Advantageously, the method comprises a step of generating a representation of the at least one indicator in the form of images, then a step of creating a collection of images from the images generated in the generation step, a step of feeding a classification tool with said collection of images, and a step of classifying the images according to the rotation of the rotating shaft and / or the anomalies present in the collection of images by the classification tool.

[0026] Advantageously, the classification tool is a neural network, preferably a convolutional neural network, even more preferably a convolutional neural network of dimensions greater than or equal to the dimensions of the images in the image collection.

[0027] Advantageously, the method comprises a step of displaying an image of at least one equivalent ellipse.

[0028] In one embodiment, steps a) to f) are carried out at least a second time for an additional plane intersecting the axis of rotation of the shaft and distinct from said plane, the method further comprising a step of comparing the characterization indicators of the planes and / or the creation of new characterization indicators common to the two planes.

[0029] The present invention also relates to a computer-implemented method comprising the steps of the method as defined above. The implementation can be done on a real-time target, in particular in a magnetic bearing control cabinet. Brief description of the figures

[0030] [Fig.l]

[0031] is a schematic representation of the different steps of the method according to the invention; and

[0032] [Fig.2]

[0033] is a representation of an example of an equivalent ellipse and its Euler parameters as determined in the method according to the invention. Detailed description of the invention

[0034] The different stages of the method according to the invention are shown schematically in [Fig.l].

[0035] The method according to the invention makes it possible to characterize the behavior of a rotating shaft. A rotating shaft is understood to mean a shaft that is mobile and rotates around an axis of rotation. The rotating shaft is, for example, a shaft of a rotating machine.

[0036] To implement the method, a step 10 of acquiring input signals is first carried out.

[0037] At least two input signals are acquired over a predetermined period. The predetermined period is, for example, several seconds.

[0038] The input signals comprise measurements representative of the movement of the rotating shaft. For example, the input signals comprise measurements of position, or speed, or acceleration, or force, or magnetic field, or even magnetic induction of the rotating shaft. The acquisition of the input signals is carried out for example using sensors configured to measure these quantities.

[0039] Generally, a vector quantity is measured according to several components distributed in space sensitive to the rotation of the shaft. The distribution in space of these components makes it possible to evaluate the movement of the rotor.

[0040] A special case is that of measurements that are non-collinear with each other and coplanar with a plane intersecting the axis of rotation of the shaft. From two non-collinear measurements, a characterization of the behavior of the rotating shaft is possible in a plane intersecting the axis of rotation of the shaft. Additional measurements allow greater precision of said characterization.

[0041] Non-coplanar measurements with a plane intersecting the axis of rotation of the shaft can nevertheless be brought back by projection into a reference plane intersecting the axis of rotation of the shaft.

[0042] For M measurements with M>2, we obtain M input signals U\(t) to uM(t), t designating the time during the acquisition of the measurements, the set of values ​​of which is designated by T, in other words dates, which t can take.

[0043] A secant plane is understood to mean a plane in which the measurements are made which intersects the axis of rotation of the rotating shaft. The secant plane is, for example, orthogonal to the axis of rotation or non-orthogonal.

[0044] In the case where the secant plane is not orthogonal to the axis of rotation, a step 20 is carried out for determining two components of the input signals, the components being orthogonal to each other and to the axis of rotation. In other words, this step 20 is a change of basis facilitating the rest of the calculations.

[0045] In the case where the secant plane is orthogonal to the axis of rotation, the previous step is intrinsically carried out.

[0046] In the case where the axis of rotation is not rigorously and precisely defined, the shaft not always being perfectly cylindrical and / or balanced, and / or non-deformable and / or known, said orthogonal secant plane is assimilated to a study plane assumed to be orthogonal to the axis of rotation, or at least sensitive to the rotational movement of the shaft.

[0047] In practice, step 20 of determining the two orthogonal components ux(t) and uy(t) is carried out by averaging projections of the input signals at uM{t) onto two orthogonal axes x and y, these two axes being ideally orthogonal to the axis of rotation. Thus,

[0048] ux ( t) = avg{ ( Projx yt^T ' uy(t) = moy^(Projy

[0049] The components ux ( t ) and uy ( t ) are therefore now in the study plane, ideally orthogonal to the axis of rotation.

[0050] Then, a step 30 is carried out to merge the two orthogonal components llx ( / ) and uy(t) into a bivariate signal r(t) in order to combine these two components in an equivalent manner in the form of a complex signal, such that:

[0051] VteT, r(t) = ux(t)+ ]Luy(t)

[0052] It can also be expressed as follows:

[0053] r( / ) = {ux(t),uy(t)}

[0054] Or in vector form:

[0055]

[0056]

[0057]

[0058]

[0059] We then carry out a step 40 of frequency decomposition of the bivariate signal r(t) into spectral elements characterized by their time-frequency pair (t, f), also called time-frequency cells, and this periodically in order to obtain a decomposition that is both temporal and frequency, called time-frequency decomposition s(t, f) such that: s(t,f) = spectrum There are different ways to perform this time-frequency decomposition step 40. These different implementation modes are presented below, the time-frequency decompositions s, or orA, or "being given in continuous form, in other words analog and not discrete. A digital implementation implemented by computer then includes a time and frequency sampling step. In a first mode of implementation, noted 41 in [Fig.l], the time-frequency decomposition step 40 is carried out with a filter bank around each frequency f.

[0060]

[0061]

[0062]

[0063]

[0064]

[0065]

[0066]

[0067]

[0068]

[0069]

[0070]

[0071] Thus, the time-frequency decomposition is expressed: with hf the impulse response of a bandpass filter centered around the frequency f, F+ being the set of frequencies greater than or equal to 0. In a second embodiment 42, the time-frequency decomposition step 40 is carried out by Fourier transform, in particular short-term single-sided Fourier transform. Basically, the time-frequency decomposition is expressed: ( af ) = i -Ca (T ) V / >0, t In practice, we prefer to use more efficient estimators based on a division into time intervals, with time weightings, averaging over several time intervals, and overlaps between successive intervals. The formulation is then: Sxk (t, f) = J UX(.T) Wk (Tt) ¥ t Y f eF^ J +' çt+dt,. Syk Uy .e^^T sM f) = vrer.V / eF w With dtk and dtk for k = [ 1, / ?] defining a division into p intervals, with possible overlap, around a given date t. Classically, this division does not depend on t, which is identical for each cell. The temporal extent of the interval k is dtk = dtk + dtk. The overall time extent of the corresponding cell ( t,f ) is dt+ = max^ipy ( dtk ) dt = maxk=M( dtk ) dt = dt + dt+ And wk(r) for rE [ -dtk, dtk] for k = [1, p] defining a time weighting window for each of these intervals, for example of Hanning type. Classically, does not depend on k, identical for each interval.

[0072]

[0073]

[0074]

[0075]

[0076]

[0077]

[0078]

[0079]

[0080]

[0081]

[0082]

[0083]

[0084]

[0085]

[0086]

[0087] Depending on the estimated spectral quantity, it is necessary to introduce the following compensations: yke [1,p], [1, p], In a special case, the absence of weighting corresponds to w*=l, in which case 0 == . kkk In addition, we introduce a subset of dates: = {4 i = [ 1, m']} In a third embodiment 43, the time-frequency decomposition step 40 is carried out by Fourier transform, in particular short-term double-sided Fourier transform, also called “Fullspectrum” in English terms, in other words for a set F of frequencies taking values ​​in the set of real numbers. So, basically the time-frequency decomposition is expressed: Vf^R, (T)£-^ / ^7 Or with a division into time blocks: d+dtk YteTi, Y f eF sk(t,f) =^7(7)™^ VtETt,yfEF, s(t, f) = avg,.. In a fourth embodiment 44, the time-frequency decomposition step 40 is carried out by Fourier transform, in particular short-term quaternion Fourier transform, involving the unitary imaginary numbers verifying li2 = lj2 = 1k2 = li.lj.lk = -1. Thus, the basic time-frequency decomposition y) is expressed as: V / > 0, S = i ( T ) £^-^7 Or with a division into time blocks: fi d+dtk Vtg V f F+, Sk(t,f)=]t_dtkr(.7}.wk(7-t}£r^fr.d7 Yf^F+, s(t, f) The expressions for the Fourier transforms given in implementations 41 to 44 above constitute applicable formulas, producing estimators efficient, particularly in their formulation with time weightings, averaging over several time intervals, overlaps between successive intervals, for which a more general expression can be given:

[0088] . j (t+dtfc

[0089] being the time weighting, and being defined such that l / / 2 = - 1.

[0090] Finally,

[0091] "FU / ) ) ) n

[0092] The index k applied to wk (above) or p k (below) means that the weighting can differ from one interval to another.

[0093] In addition, other frequency transforms can be used, such as wavelet decompositions.

[0094] A more general expression is:

[0095] , rt+dt^ % (^f ) = kdfkr

[0096] More generally still, it is possible to vary the time division and the weighting as a function of the time-frequency cell (t, f):

[0097] 4+dt^tf\ hf) = h-dfk^r^-Pk^ f)

[0098] The trajectory (orbit) is obtained by inverse transformation, of the type:

[0099] rf+df(t,p__ __ s(t,f)(T) ( s (^)4((, / )(7^) + s (t, -v)4(tf)(T, -v))dv

[0100] The frequency range of the cell (t, f) then being df — df+ + df\ the time range of the cell (t, f) being dt = dt+ + dt

[0101] From each spectral element s(t, f), at a fixed date t and frequency f, we can deduce from the instantaneous orbit an equivalent ellipse. We can describe it as local in the temporal and frequency sense, in other words at date t and frequency f.

[0102] A step 50 is therefore carried out to determine the equivalent ellipse for each spectral element. Whatever the time-frequency decomposition method, the instantaneous orbit can be expressed in the form of a parametric equation such that, for each component x and y, and with Sx\ and 0X, 0y the instantaneous amplitudes and phases:

[0103] sx(t,f) — sXa(t, f)£os(0x(t, f))

[0104] (t,f)=Sya(t,f) £OS (0y (t,f))

[0105] For a fixed frequency, going through time t allows obtaining the trajectory at this frequency.

[0106] This instantaneous orbit takes the form of an equivalent ellipse, provided that sx(t,f),sy(tf)), 0x(t,f ), 0y(t, f ) vary slowly as a function of time, by

[0107]

[0108]

[0109]

[0110] [YES]

[0112]

[0113]

[0114]

[0115]

[0116]

[0117]

[0118]

[0119]

[0120] example during a complete revolution at frequency / , and in all cases during the duration dt of the time interval defining the spectral element. In this case, we can distinguish the date t of the spectral element considered from the time T which evolves in this spectral element and, by introducing the frequency / in the expression of the phases assumed to be quasi-linear, we obtain the equivalent ellipse in the form of a temporal parametric equation for each component x and y: (T) = ax(t, f )£Os(2jr.fT + (pjt, f ) ) Sy(t,f)(T) -ay(t,f}£OS{2jT.fT + (py{t,f) ) Another writing of these parametric equations is for example: sxdJ)(r) = axcos(2rr.fr)£os(rpx(tf)) -sin(2rr.fr)sin(rpx(t, f ) ) ) sy (t,f)(r)=ay( t, f ). ( cos ( 2jr.fr ) £os (>v ( t, f))- sin ( 2jr.fr ) sin ( <py ( t, f ) ) ) To obtain these parametric equations, particularly in the case of a frequency decomposition by filter bank, the parameters ax(t, f), Çx(t, f), ay(^ f),<Pv(t> f ) of the equivalent ellipse can be obtained by estimations, such as direct or quadratic averaging or by effective value, from the instantaneous quantities ( t, f ) for T &[t-dt,t + dt+}- In the case of a frequency decomposition by Fourier transform, we obtain the parametric equation of the equivalent ellipse by inverse Fourier transform of the frequency element: sx (t. f) (r) - 2Re(sx(t,f)£^-fr) Sy (t, f) (T) = ZRe ) By forming the associated complex quantity s ( t, f ) = sx{t, f) + Usy {t. f), in other words the bivariate representation, we highlight the decomposition of a parametric equation of the equivalent ellipse into two circles of opposite directions: sY f) (Tï +a+(j, f) With the following parameters obtained following a calculation step 59: a+ = f .\jâŸÏ-d^^ has. = + a$ + Z^MySin ( (py -(px^ 3+ = ardanZ^ axsin ( (px ) + ay£os (cpy), a^cos ( <px ) - aysin [(py)') 6- = arctanZ(üxSin( <px) -ay£os((py), ax£os(<px) + aysin(<PyY) The characterization of the equivalent ellipse is natural with the double-sided Fourier transform (“fullspectrum”) which can be directly written in polar form: s(t,f) =«+( W) > ' s(t, - f ) = u(t, f

[0121] Generally, an equivalent ellipse can be determined and plotted by obtaining its Euler parameters.

[0122] Figure 2 shows an example of an ellipse and its Euler parameters a(t, f), δ(t, f), f ), which completely define an equivalent ellipse, with:

[0123] ">0 Oe [- y2, y2 ] Z ^ [ - V4, V4 ] (pe [-TT, tf]

[0124] Furthermore, b is the semi-major axis, in other words the major radius, and c is the semi-minor axis, in other words the minor radius, the sign of which gives the direction of rotation of the equivalent ellipse.

[0125] The parametric equation of the trajectory of the equivalent ellipse is more convenient in complex form, because:

[0126] - we express a circle equation, notably in cos and sin of the same phase,

[0127] - the axes are dilated to obtain an ellipse,

[0128] - a rotation is applied to it to obtain the desired direction, thus (r> =a(t, ff ) )eos(ln.f:r+<p(t,f) ) + )sin(2;r.f.T + f[> (t, f) ) )

[0129] Direction means the angle 0 of the equivalent ellipse between the major axis and the x axis. Orientation means the direction of rotation describing the trajectory of said ellipse.

[0130] This characterization of the equivalent ellipse is natural with the time-frequency decomposition by quaternion Fourier transform which can be directly put in Euler form:

[0131] f) = | s(f,f) | .ele

[0132] This parametric equation has the particular interest of rejecting the time parameterization towards a single exponential term when applying the inverse Fourier transform to the frequency element in order to obtain the temporal parametric equation of the trajectory of the equivalent ellipse:

[0133] (r) ~ f ) LProf

[0134] with Projç a projection into the set of complex numbers C lr

[0135] More precisely, step 50 of determining the equivalent ellipse includes a step of calculating Euler parameters of said equivalent ellipse.

[0136] Starting from the time-frequency decomposition by filter bank, we obtain a calculation step 51 of ax and ay such that:

[0137] / _ rt+dt+ v ax(t, / ) = ^Jtdf si(r, f)dr -- - ! rt+dt Uy(F f) = Jt_dt s}{r,f)dr

[0138] Ti being a subset of T.

[0139]

[0140]

[0141]

[0142]

[0143] Starting from the time-frequency decomposition by short-term single-sided Fourier transform, we obtain a step 52 of calculation of ax, ay, Vx and (Py such that: <Px(Ff) ) V / eM / ef* / )=2^(( / )| v^f)=Arg(.sy(lf)-) Starting from the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 53 of calculation of a+, a\ G+ and such that: * a+(t,f) = \s(t,f)\ e+(t,f)=Arg(s(t,f)) WinJet* e_^ f) = . Arg(S(t,-f)) By an additional calculation step 54 or 55, we obtain the following Euler parameters:

[0144] a = ^2 + a? (^+6. <p = ~ x="4" - arctanz (a+-a^a+ + a.) ou ‘iy+v, (p="~T~" q — 4 .arctanz za^ycos ((px-tpyyyày-ajy arcsin ( ratio* znxny,   + aj ) sin      )

[0145]

[0146]

[0147]

[0148]

[0149] l’étape 56 de calcul inverse est également possible, tel que : a+~i.{cos(xy +sin(x) a.="l.(cos(xy" -ùn(x) 3+="(p" +8 8.="(p-6" le procédé selon l’invention comprend optionnellement une étape 60 la densité spectrale décomposition temps-fréquence. les densités spectrales chaque composante et croisées entre composantes permettent d’estimer répartition énergétique les x y. en partant temps-fréquence par banc filtre, on obtient 61 5’ telle । a+dp sx u j t, f dr „                                   , (t+df sy(af) çt+dt*' sxy(t, f)="df" sx(t, f)sy(t>f}dr

[0150]

[0151]

[0152]

[0153]

[0154]

[0155]

[0156]

[0157]

[0158]

[0159]

[0160]

[0161]

[0162]

[0163]

[0164]

[0165] Starting from the time-frequency decomposition by short-term single-sided Fourier transform, we obtain a step 62 of calculation of the spectral density s such that: ^ / )-24 / ^^^(^-) 5 / ( / )=24 / ^^(^) has / \ f) ) Starting from the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 63 of calculation of the spectral density S such that: VZeTt, V / ef S(tf)=dfmoyk^-^'} We will also use: V fe V / e S* (t, f) = dfjnoyk^p] Starting from the time-frequency decomposition by short-term quaternion Fourier transform, we obtain a step 64 of calculation of the spectral density § such that: yf^F+, S(t,f)=djjn0yHi^----) From these spectral densities^', s, S, and $, we carry out a step 70 of calculating the Stokes parameters So, S2, and $3 of at least one equivalent ellipse. The Stokes parameters are four real numbers and are classically used to describe the polarization state of an electromagnetic wave. S 0 is commonly the total intensity of an electromagnetic wave and represents here by extension the total intensity of the bivariate signal, S 0 >0. S i and S 2 give the linear polarization intensity »2 as well as the management of the polarization of the ellipse # = ^tanl^ • S 3 gives the circular polarization intensity |53|, its sign giving the direction of rotation. We also define the reduced Stokes parameters: Life {1,2,3} | Starting from the spectral density obtained for the time-frequency decomposition by filter bank, we obtain a step 71 of calculating the Stokes parameters of at least one equivalent ellipse such that:

[0166]

[0167]

[0168]

[0169]

[0170]

[0171]

[0172]

[0173]

[0174]

[0175] . 50((, / )=2(5,((, / )+5,((, / )) Starting from the spectral density obtained for the time-frequency decomposition by short-term single-sided Fourier transform, we obtain a step 72 of calculating the Stokes parameters of at least one equivalent ellipse such that: V t £ TV f F. Starting from the spectral density obtained for the time-frequency decomposition by short-term double-sided Fourier transform, we obtain a step 73 of calculating the Stokes parameters of at least one equivalent ellipse such that: S o ((, / ) =2(5((, / )+S(A - / )) „ 5,((, / ) =4JZe(5*((, / ) ) V t EF, V f Ê F+. , ~ , * S2(t, / )=4Jm(S.(t, / )) S3(t, / )=Z(S(t, / )-S(t, - / )) Starting from the spectral density obtained for the time-frequency decomposition by short-term quaternion Fourier transform, we obtain a step 74 of calculating the Stokes parameters of at least one equivalent ellipse such that: / o \ 5,((. / ) = / m, / s((, / )) (e^V / eF^ , S2((, / )= / <n,*(S(^ / )) S3(t, / ) =Imÿ(s(t, / )) Alternatively, the Stockes parameters S2, and are calculated directly using the time-frequency decomposition by quaternion Fourier transform: 5((, / ) =s(Z. / ).l / ((,7) 5,((, / ) = / 1(1, / 5((, / )) S2(t,f) = Im[k(a(t, f) ) From these Stokes parameters, we carry out a step 75 of calculating the polarization rate 0 as well as the reduced Stokes parameters previously mentioned, as well as the polarized part of the intensity Sp. [°1761 Life {1,2,3} s,(t, / ) =^1(0^(5 / 1, / ), S0(t, / )) VleT^ / eF* 0( / , / ) = ^(t, / ) +s^t, f) +s / t, f) S„(,t, / ) = 0(t, / )So(t,f)

[0177] The polarization expressed by the polarization rate 0 makes it possible to highlight the cyclo-stationarity of an equivalent ellipse, namely the regularity of the trajectory of said ellipse with respect to a period of revolution. The polarization reflects the structuring of said trajectory relative to defined frequencies. The polarization rate can thus be calculated for each time-frequency cell resulting from the time-frequency decomposition in order to highlight the cyclo-stationarity of the rotating shaft for a given frequency.

[0178] Furthermore, a calculation step 57 can be implemented to determine Euler parameters a, and from the polarization spectral density Sp and the reduced Stokes parameters: 101791 a( / , / )=^Sp(t, / ) V tG T\, Vf eF*, 3(t, f) = ^.arctan2(s2(t, f), sff f) ) / (A / ) - 5arcsin(ratw (s3(t, f ), Sp(t, f ) ) )

[0180] A fourth Euler parameter giving the phase at the origin of the trajectory can also be deduced by the following calculation step 58:

[0181] Vt^7\, VfE F*, (p(t, f) - Arg Pro( / , / ]) j

[0182] Following step 50 of determining an equivalent ellipse, in particular by calculating Euler parameters of said equivalent ellipse, a step 80 is carried out of extracting at least one characterization indicator of the rotating shaft from parameters of the equivalent ellipse, for example Euler parameters, Stokes parameters or parameters derived from the latter.

[0183] In addition, the preceding steps are optionally carried out at least a second time for an additional plane intersecting the axis of rotation of the shaft and distinct from the plane for which these steps were carried out in the first place. In this mode of implementation, it is thus possible to carry out a step of comparing the characterization indicators of the planes.

[0184] For example, the at least one characterization indicator comprises the previously described polarization rate 0, and / or a circularity indicator (linked to | / |), and / or an intensity indicator and / or an amplitude indicator, and / or a rotation direction indicator (linked to the sign of ^), and / or a direction indicator (linked to 0), and / or synchronization (linked to the equivalent ellipse.

[0185] Each elementary indicator contributes to a concrete and useful characteristic of the rotating system, for example a high polarization rate and a high circularity are probably due to the rotating shaft. Distinguishing this part of the signal, the complementary being assimilated to disturbances, can be useful for example to identify said shaft (for modeling purposes) or better control it in the framework of a shaft equipped with magnetic bearings.

[0186] Optionally, the at least one characterization indicator is normalized, in other words, its value is between -1 and 1.

[0187] The characterization indicator is given for a plan, it is then noted A, or by comparison between two plans, it is then noted I.

[0188] The characterization indicator is qualified as a criterion, noted C, when it is constructed from any characteristic quantity re-evaluated at each instant t. The characterization indicator is qualified as stationarity, noted S, when it is based on the stationarity over time of a quantity (maximum when the variation over time of said quantity is zero).

[0189] Thus, we note the indicator of the polarization rate:

[0190] Vre^V / eF*

[0191] For a plan, the intensity indicator is:

[0192] VtET^ V / eF^ ACintavi}té{t, f) ■= normalizes f Y)

[0193] With the normalization function:

[0194] vte TT fe FF normalized _„Ax(t,f)) =

[0195] TT and FF denote sets of generic dates and frequencies, taking their value in the set of real numbers R.

[0196]

[0197]

[0198]

[0199]

[0200]

[0201]

[0202]

[0203]

[0204]

[0205]

[0206] For a plane, the amplitude indicator is: V te TV fe F+, ACampiitude (t,f) = normalize^T^ (a(t,f)) For a plane, the sign indicator is -1 or 1 depending on the sign of X: V te T^, V / eF^ ACsign (t, f)^ sign (fffY For a plane, the circularity indicator is: VteT, For a plane, the indicator of stationarity of sign, in other words of direction of rotation, in other words of orientation, is: V te TV f eF+, ASsign (t,f) = autostat (AC sign (t, ff you = 0) With the stationarity function: V tE TT autostat (x(JY you) - e-^r pourtol-Q: autostat(x(t), G) = (dx(t) = =0) For a plane, the circularity stationarity indicator is:

[0207] ASdrMpj, f) ^aidostatt^

[0208] For a plane, the indicator of direction stationarity, in other words of geometric inclination of the ellipse, is:

[0209] V t EV F+, ASdirection (t, f) = autostatf ), you = W / 180)

[0210] For a plan, the synchronization stationarity indicator is:

[0211] VZeT|, V f &F+, ASsynclm, (t, f) = autostat (tp(t, f), you =10 jï / 180)

[0212] For two planes, we redefine each parameter or indicator for an additional plane, by noting these new parameters using an apostrophe, for example a

[0213] The comparative characterization indicator of the amplitude is:

[0214] yt^yfeF., ICampütude(t,f) = normalizes^T„i^F+^a(t, f)a'(t, f))

[0215] The comparative characterization indicator of the sign is a boolean between 0 and 1 according to the result of the following test:

[0216] vtevfe F„ IC^ (1, / ) = (AC^ ( / / ) = = AC'^ ((, / ))

[0217] The comparative circularity characterization indicator is:

[0218] yteT» v / eF, ICAMdt. f)

[0219] With the relative error function such as -.reldif(x,y) = ratio(2. |x - y|, I xi + |y|)

[0220] And with the division function with exception such as: ax IV if y & 0 ratio (x, y ) = ' J 0 otherwise

[0221] The indicator of comparative characterization of direction, in other words of geometric inclination of the ellipse, is:

[0222] v / eT„ v / eF., ) =

[0223] The comparative synchronization characterization indicator is:

[0224] V tg 7}, V / eF* ICsynchro(t, f) = 1 -reldif (tp(t, f),tp (t, f))

[0225] The indicator for characterizing stationarity and comparative circularity is:

[0226] VfEF¥ ISareularityO'f) = OUtOStat(AC^^^t, f) - ACcimilar^

[0227] The indicator of stationarity characterization and comparative direction is:

[0228] vte V / eF+, ISdMn(U) = autostat f)-0(t, f), tol = IŒt / 180)

[0229] The indicator for characterizing stationarity and comparative synchronization is:

[0230] V tg Ty y ISsynchm (t,f) = autostat f) - <p'(t, f), toi- 10jf / 180)

[0231] Other indicators can be constructed from the quantities constructed in steps 70 and 50.

[0232] Optionally, a step 90 is carried out of inserting the at least one characterization indicator into a neural network in order to classify said indicator into function of the shaft rotation and link it to a previously listed anomaly. Anomalies similar to the learning base are then detected.

[0233] Optionally, a step 100 is carried out for displaying the at least one extracted characterization indicator to be readable by an operator and to enable his understanding of potential anomalies of the rotating shaft.

[0234] Optionally, a step 101 is carried out for generating a representation of the at least one indicator in the form of images.

[0235] The images are for example orbits, or time-frequency representations, or phase representations according to frequency, or colored images, the color representing the intensity, the direction of rotation or another indicator.

[0236] The images may also be composite or hybrid images, for example in the form of a grid, for example with several levels of decomposition, or non-uniform. The image thus has two main axes, namely abscissa and ordinate, corresponding to two variation parameters, for example time and frequency, while the two axes of the cells of the grid correspond to two other variation parameters, for example Cartesian coordinates of a plane orthogonal to the axis of rotation.

[0237] Thus, in one embodiment, the grid comprises a set of cells each associated with an image representing one or more indicators, each cell being arranged according to a time-frequency mapping.

[0238] In another embodiment, the grid comprises a set of cells each associated with an image representing one or more indicators, such as orbits superimposed on force vectors, each cell being arranged according to a representation in Cartesian coordinates. The composite or hybrid images can also be obtained by superimposing several quantities, for example orbits and / or force fields, in the same system of axes, for example the Cartesian coordinates of a plane orthogonal to the axis of rotation, the different quantities possibly being distinguished or supported by shapes or graphic symbols, for example arrows to materialize a vector or a direction.

[0239] Then a step 102 is carried out for creating a collection of images from the images generated in the generation step 101, the images being able to be ordered according to a chosen and / or predetermined order.

[0240] A step 103 of feeding a classification tool with said collection of images is then carried out, as well as a step 104 of classifying the images by the classification tool as a function of the rotation of the rotating shaft and / or the anomalies present in the collection of images.

[0241] In a particularly advantageous embodiment, the classification tool is a neural network, preferably a neural network convolutional, particularly effective for image classification. Even more preferably, the convolutional neural network has dimensions greater than or equal to the dimensions of the images in the image collection.

[0242] For example, if the image collection contains only one two-dimensional image, then this amounts to classifying a single two-dimensional image and a two-dimensional convolutional neural network can be used.

[0243] In another example, if the image collection is used to represent a three-dimensional image, a three-dimensional convolutional neural network may be used.

[0244] In a particular embodiment, each image of the collection consists of one or more pixels grouped according to the chosen resolution, representing the position of the rotor in a plane orthogonal fixed to the longitudinal axis of the rotor at a given instant, and the collection of images contains all the images ordered along the axis of increasing time. In this case, if for example the orbit of the rotor described in the fixed orthogonal plane is circular, the collection of images will describe a three-dimensional image representing a helical shape and a three-dimensional convolutional neural network can be used.If the image collection contains all the images described above but for all planes orthogonal to the rotor's longitudinal axis ordered in order along that axis and for all instants ordered in increasing time, then the image collection would describe a four-dimensional image, which can be classified by a four-dimensional convolutional neural network.

[0245] The collections of images feeding the neural network make it possible to form a learning and training base for the neural network, the classification function of which improves as it receives collections of images.

[0246] The classification step 104 optionally makes it possible to activate alarms making it possible to control or protect the rotating shaft and / or its environment, in particular when significant anomalies are detected during said classification step 104.

[0247] A step 110 of displaying an image of the trajectory of at least one equivalent ellipse can also be carried out.

[0248] In a particular mode of implementation, all the steps of the method are implemented by computer, with the possibility of execution in real time, in particular in a magnetic bearing control cabinet.

[0249] Certain terms and functions used in the previously detailed calculations are redefined below:

[0250] The set of strictly positive frequencies is: / 7* = | yj — [ n] j

[0251] The set of strictly negative frequencies is: p* - _p

[0252] We also define = {0} u F*

[0253] We also define: F- = FU {0}

[0254] The complete set of frequencies is: F — F. U F+ = F. U {0} UF^

[0255] The generic stationarity function is y For you >0: autostat(x(t),tol) = pourtol = Q: autostat (x (t), 0) = (dx(t) = = 0)

[0256] The maximum function with exception is: . max(x) if max(x)*0 max(x) - , . 1 otherwise

[0257] The normalization function vte TT, Y f G FF, normalizes (x(t,f)) =

[0258] The division function with exception is: * xly gj y^Q ratio (x, y) = ' J 0 otherwise

[0259] The relative error function is:

[0260] reldif(x, y) -ratw (2.\xy\, IxI + |y | )

[0261] Let us also recall the definition of the arctan2 function, a formulation equivalent to the argument of the complex number X+ \ny:

[0262] Arg (x + \(iy) - arctanl (x)

Claims

Claims

1. Method for characterizing the behavior of a rotating shaft, characterized in that it comprises the following steps: - a) acquisition over a predetermined period of at least two input signals, the input signals comprising measurements representative of the movement of the rotating shaft, the measurements being non-collinear to each other (step 10); - b) determination of two components of the input signals, the components being orthogonal to each other (step 20); - c) merging of the two components into a bivariate signal (step 30); - d) time-frequency decomposition of the bivariate signal into spectral elements (step 40); - e) determination of an equivalent ellipse for each spectral element (step 50); and - f) extracting at least one characterization indicator of the rotating shaft from parameters of the equivalent ellipse (step 80), wherein the at least one characterization indicator comprises a polarization rate.

2. A method according to claim 1, wherein the input signals comprise measurements of position, or velocity, or acceleration, or force, or magnetic field, or magnetic induction of the rotating shaft.

3. Method according to one of claims 1 and 2, in which the time-frequency decomposition step (40) is carried out by Fourier transform, or by wavelet, or by filtering, or by quaternionic Fourier transform.

4. A method according to any one of claims 1 to 3, wherein the step (50) of determining an equivalent ellipse comprises calculating Euler parameters of said equivalent ellipse.

5. Method according to any one of claims 1 to 4, comprising a step (60) of calculating the spectral density of the time-frequency decomposition, and a step (70) of calculating Stokes parameters from said spectral density.

6. Method according to any one of claims 1 to 5, wherein the at least one characterization indicator further comprises an indicator of circularity and / or amplitude and / or direction of rotation and / or sign and / or synchronization of an equivalent ellipse.

7. Method according to any one of claims 1 to 6, further comprising a step (90) of inserting the at least one characterization indicator into a neural network to classify said indicator according to the rotation of the shaft and link it to a previously listed anomaly.

8. Method according to any one of claims 1 to 7, further comprising a step (100) of displaying the at least one extracted characterization indicator to be readable by an operator.

9. Method according to any one of claims 1 to 8, comprising a step (101) of generating a representation of the at least one indicator in the form of images, then a step (102) of creating a collection of images from the images generated in the generation step (101), a step (103) of feeding a classification tool with said collection of images, and a step (104) of classifying the images according to the rotation of the rotating shaft and / or the anomalies present in the collection of images by the classification tool.

10. Method according to claim 9, in which the classification tool is a neural network, preferably a convolutional neural network, even more preferably a convolutional neural network of dimensions greater than or equal to the dimensions of the images in the image collection.

11. Method according to any one of claims 1 to 10, comprising a step (110) of displaying an image of at least one equivalent ellipse.

12. Method according to any one of claims 1 to 11, in which steps a) to f) are carried out at least a second time for an additional plane intersecting the axis of rotation of the shaft and distinct from said plane, the method further comprising a step of comparing the characterization indicators of the planes and / or the creation of new characterization indicators common to the two planes.

Images