METHOD FOR CHARACTERISING THE BEHAVIOUR OF A ROTATING SHAFT USING EQUIVALENT ELLIPSES

Abstract

A method for characterising the behaviour of a rotating shaft using equivalent ellipses. The method for characterising the behaviour of a rotating shaft comprises the following steps: acquiring at least two input signals over a predetermined period comprising measurements representing the movement of the rotating shaft (step 10); determining two components of the input signals, the components being orthogonal to each other (step 20); merging the two components into a bivariate signal (step 30); time-frequency decomposing the bivariate signal into spectral elements (step 40); determining an equivalent ellipse for each spectral element (step 50); and extracting at least one characterisation indicator of the rotating shaft based on parameters of the equivalent ellipse (step 80).

Description

CROSS-REFERENCE TO RELATED APPLICATION

This application claims priority to French Application No. 2310162, filed Sep. 25, 2023, the entirety of which is hereby incorporated by reference.

FIELD

The present disclosure relates to the characterisation of a shaft rotating about an axis of rotation.

More specifically, the present disclosure relates to the characterisation of the gyroscopic effects and deformations of a rotating shaft in order to facilitate the control, modelling, monitoring, diagnosis and maintenance thereof.

BACKGROUND

A perfectly balanced rotating system, such as a rotating shaft, rotates about an axis of rotation such that said axis of rotation is immobile relative to its theoretical position, and therefore relative to the environment close to the rotating shaft.

However, faults, deformations, gyroscopic effects or even effects due to the force of gravity can cause the axis of rotation to shift relative to its theoretical position during the rotation of said rotating shaft. This is the case, for example, for deformations of the rotating shaft due to flexible Eigen-modes, with these modes deforming the shaft as a function of the speed of rotation of said shaft.

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

The properties of a rotating shaft are conventionally understood through 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.

The changes in these quantities often need to be characterised depending on whether or not the shaft is rotating. However, many causes can influence the rotational behaviour of a rotating shaft and it is difficult to precisely characterise such behaviour.

SUMMARY

Therefore, the aim of the present disclosure is to overcome the aforementioned disadvantages.

The aim of the present disclosure is a method for characterising the behaviour of a rotating shaft, the method comprising the following steps:

• • a) acquiring at least two input signals over a predetermined period, the input signals comprising measurements representing the movement of the rotating shaft, the measurements, associated with a vector quantity, being non-collinear to each other and optionally coplanar with a plane intersecting the axis of rotation of the shaft;
  • b) determining two components of the input signals, the components being orthogonal to each other and optionally to the axis of rotation;
  • c) merging the two components into a bivariate signal;
  • d) time-frequency decomposing the bivariate signal into spectral elements;
  • e) determining an equivalent ellipse for each spectral element; and
  • f) extracting at least one characterisation indicator of the rotating shaft based on parameters of the equivalent ellipse.

The gyroscopic effects and the deformations of the rotor, notably the flexible modes, are thus simpler to consider subject to 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 characterised using a frequency decomposition that can highlight a rotating orbit, also called an equivalent ellipse, for each frequency, with each orbit rotating at a specific speed corresponding to the considered frequency.

The present disclosure therefore enhances the characterisation of a rotating shaft, allowing applications of the control, identification, modelling, monitoring or diagnostic type.

Advantageously, the input signals include measurements of the position, or speed, or acceleration, or force, or magnetic field, or magnetic induction of the rotating shaft.

In various embodiments, the time-frequency decomposition step is carried out using a Fourier, or wavelet, transform, or by filtering, or using a quaternionic Fourier transform.

In a particular embodiment, the step of determining an equivalent ellipse includes computing Euler parameters of said equivalent ellipse.

Advantageously, the method comprises a step of computing the spectral density of the time-frequency decomposition, and a step of computing Stokes parameters based on said spectral density.

In one embodiment, the at least one characterisation indicator includes a polarisation rate, and/or an indicator of circularity and/or amplitude and/or orientation and/or direction of rotation and/or sign and/or synchronisation of an equivalent ellipse.

Advantageously, the method comprises a step of introducing the at least one characterisation indicator into a neural network in order to classify said indicator as a function of the rotation of the shaft and to link it to a previously listed anomaly.

In a particular embodiment, the method comprises a step of displaying the at least one extracted characterisation indicator so that it can be read by an operator.

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 based on the images generated in the generation step, a step of providing a classification tool with said collection of images, and a step of classifying, using the classification tool, the images as a function of the rotation of the rotating shaft and/or of the anomalies present in the collection of images.

Advantageously, the classification tool is a neural network, preferably a convolutional neural network, even more preferably a convolutional neural network with dimensions that are greater than or equal to the dimensions of the images of the collection of images.

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

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 characterisation indicators of the planes and/or creating new characterisation indicators common to the two planes.

A further aim of the present disclosure is a computer-implemented method comprising the steps of the method as defined above. The method can be implemented on a real-time target, notably in a magnetic bearing control cabinet.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a schematic representation of the various steps of the method according to the present disclosure; and

FIG. 2 is a representation of an example of an equivalent ellipse and its Euler parameters as determined in the method according to the present disclosure.

DETAILED DESCRIPTION

FIG. 1 schematically shows the various steps of the method according to the present disclosure.

The method according to the present disclosure is used to characterise the behaviour of a rotating shaft. A rotating shaft is understood to be a shaft that is free to rotate about an axis of rotation. The rotating shaft is a shaft of a rotating machine, for example.

In order to implement the method, a step 10 of acquiring input signals is initially carried out.

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

The input signals include measurements representing the movement of the rotating shaft. For example, the input signals include measurements of the position, or speed, or acceleration, or force, or magnetic field, or even of the magnetic induction of the rotating shaft. The input signals are acquired, for example, using sensors configured to measure these quantities.

In general, a vector quantity is measured in terms of several spatially distributed components sensitive to the rotation of the shaft. The spatial distribution of these components is used to evaluate the movement of the rotor.

A particular case is that of measurements that are non-collinear to each other and are coplanar with a plane intersecting the axis of rotation of the shaft. Two non-collinear measurements can be used to characterise the behaviour of the rotating shaft in a plane intersecting the axis of rotation of the shaft. Additional measurements allow greater precision for said characterisation.

Measurements that are not coplanar with a plane intersecting the axis of rotation of the shaft nevertheless can be returned by projecting into a reference plane intersecting the axis of rotation of the shaft.

For M measurements with M≥2, M input signals u₁(t) to u_M(t) are obtained, with t denoting the time covering the acquisition of the measurements, the set of values of which is denoted by T, in other words dates, which t can assume.

An intersecting plane is understood to be a plane in which the measurements are taken that intersect the axis of rotation of the rotating shaft. The intersecting plane is, for example, orthogonal to the axis of rotation or is not orthogonal to the axis of rotation.

If the intersecting plane is not orthogonal to the axis of rotation, a step 20 is carried out of determining two components of the input signals, with 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 subsequent computations.

In the event that the intersecting plane is orthogonal to the axis of rotation, the previous step is intrinsically carried out.

In the event that the axis of rotation is not strictly and precisely defined, with the shaft not always being perfectly cylindrical and/or balanced, and/or non-deformable and/or known, said orthogonal intersecting plane is treated as a study plane that is assumed to be orthogonal to the axis of rotation, or at least sensitive to the rotational movement of the shaft.

In practice, the step 20 of determining the two orthogonal components u_x(t) and u_y(t) is carried out by averaging the projections of the input signals u₁(t) to u_M(t) over two orthogonal axes x and y, with these two axes ideally being orthogonal to the axis of rotation. Thus,

$\forall t\in T,\begin{array}{c}{u}_x\left(t\right)={\mathrm{moy}}_{i=\left[1,M\right]}\left({\mathrm{Proj}}_x\left(u_i\left(t\right)\right)\right)\\ u_y\left(t\right)={\mathrm{moy}}_{i=\left[1,M\right]}\left({\mathrm{Proj}}_y\left(u_i\left(t\right)\right)\right)\end{array}$

the components u_x(t) and u_y(t) are now in the study plane, ideally orthogonal to the axis of rotation.

A step 30 is then carried out of merging the two orthogonal components u_x(t) and u_y(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:

$\forall t\in T,$ $r\left(t\right)=u_x\left(t\right)+1i·u_y\left(t\right)$

It also can be expressed as follows:

$r\left(t\right)=\left\{u_x\left(t\right),u_y\left(t\right)\right\}$

Or even in vector form:

$r\left(t\right)=\left[\begin{array}{c}{u}_x\left(t\right)\\ u_y\left(t\right)\end{array}\right]$

A step 40 is then carried out of frequency decomposing the bivariate signal r(t) into spectral elements characterised by their time-frequency pair (t, f), also called time-frequency cells, with this step being carried out periodically in order to obtain both a temporal and a frequency decomposition, called time-frequency decomposition s(t, f), such that:

$s\left(t,f\right)={\mathrm{spectre}}_{t,f}\left(r\left(\tau \right)\right)$

This time-frequency decomposition step 40 can be carried out in various ways. These various embodiments are described hereafter, with time-frequency decompositions s, or ŝ, or {tilde over (s)}, or s̊ being provided in continuous form, in other words in analogue and not discrete form. A computer-implemented digital implementation then includes a time and frequency sampling step.

In a first embodiment, denoted 41 in FIG. 1, the time-frequency decomposition step 40 is carried out with a filter bank around each frequency f.

Thus, the time-frequency decomposition s(t, f) is expressed as follows:

$\forall f\in F_{+},$ $\forall t\in T,$ $s\left(t,f\right)=r\left(t\right)*h_f\left(t\right)$

with h_f being the impulse response of a bandpass filter centred around the frequency f, with F₊ being all the frequencies that are greater than or equal to 0.

In a second embodiment 42, the time-frequency decomposition step 40 is carried out using a Fourier transform, in particular a short-term single-sided Fourier transform.

As standard, the time-frequency decomposition is expressed as:

$\forall f\ge 0,\begin{array}{c}{\widehat{s}}_x\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^t{u}_x\left(\tau \right)·e^{-1i·2\pi ·f·\tau }·d\tau \\ {\widehat{s}}_y\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^t{u}_y\left(\tau \right)·e^{-1i·2\pi ·f·\tau }·d\tau \end{array}$

In practice, more efficient estimators are preferably used, based on a division into temporal intervals, with time weighting, averaging over several temporal intervals, and overlaps between successive intervals. The formulation is then:

$\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}{\widehat{s}}_{x_k}\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}{u}_x\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i·2\pi ·f·\tau }·d\tau \\ {\widehat{s}}_{y_k}\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}{u}_y\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i·2\pi ·f·\tau }·d\tau \end{array}$ $\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}{\hat{s}}_x\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\widehat{s}}_{x_k}\left(t,f\right)}{w_k^{\left(1\right)}}\right)\\ {\hat{s}}_y\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\widehat{s}}_{y_k}\left(t,f\right)}{w_k^{\left(1\right)}}\right)\end{array}$

With dt_k⁻ and dt_k⁺ for k=[1, p] defining a division into p intervals, with possible overlap, around a given date t. Conventionally, this division does not depend on t, and is identical for each cell. The temporal extent of interval k is dt_k=dt_k⁻+dt_k⁺.

The overall temporal extent of the corresponding cell (t,f) is:

${\mathrm{dt}}^{+}={\max }_{k=\left[1,p\right]}\left({\mathrm{dt}}_k^{+}\right)$ ${\mathrm{dt}}^{-}={\max }_{k=\left[1,p\right]}\left({\mathrm{dt}}_k^{-}\right)$ $\mathrm{dt}={\mathrm{dt}}^{-}+{\mathrm{dt}}^{+}$

With Et w_k(t) for τ∈[−dt_k⁻, dt_k⁺] for k=[1, p] defining a temporal weighting window for each of these intervals, for example, of the Hanning type. Conventionally, w_k does not depend on k, which is identical for each interval.

Depending on the estimated spectral quantity, the following compensations need to be introduced:

$\forall k\in \left[1,p\right],\begin{array}{c}{w}_k^{\left(1\right)}={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}{w}_k\left(\tau \right)·d\tau \\ w_k^{\left(2\right)}={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}{w}_k^2\left(\tau \right)·d\tau \end{array}\forall k\in \left[1,p\right],\begin{array}{c}{w}_k^{\left(1\right)}={\int }_{-{\mathrm{dt}}_k^{-}}^{{\mathrm{dt}}_k^{+}}{w}_k\left(\tau \right)·d\tau \\ w_k^{\left(2\right)}={\int }_{-{\mathrm{dt}}_k^{-}}^{{\mathrm{dt}}_k^{+}}{w}_k^2\left(\tau \right)·d\tau \end{array}$

In a particular case, the absence of weighting corresponds to w_k=1, in which case w_k⁽¹⁾=w_k⁽²⁾=dt_k.

In addition, a sub-set of dates is introduced: T_∇={t_i′, i=[1, m′]}

In a third embodiment 43, the time-frequency decomposition step 40 is carried out using a Fourier transform, in particular a short-term two-sided Fourier transform, also called “Full spectrum” transform, in other words, for a set F of frequencies taking values from the set of real numbers.

Thus, as standard the time-frequency decomposition š (t, f) is expressed as:

$\forall f\in R,\tilde{s}\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^{t}r\left(\tau \right)·e^{-1i·2\pi ·f·\tau }·d\tau$

Or with a division into temporal blocks:

$\forall t\in T_{\downarrow },\forall f\in F,\widehat{s}\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1i·2\pi ·f·\tau }·d\tau$ $\forall t\in T_{\downarrow },\forall f\in F,\widehat{s}\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\tilde{s}}_k\left(t,f\right)}{w_k^{\left(1\right)}}\right)$

In a fourth embodiment 44, the time-frequency decomposition step 40 is carried out using a Fourier transform, in particular a short-term quaternionic Fourier transform, using the unitary imaginary numbers verifying 1i²=1j²=1k²=1i.1j.1k=−1.

Thus, as standard the time-frequency decomposition s(t, f) is expressed as:

$\forall f\ge 0,\stackrel{o}{s}\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-\mathrm{dt}}^{t}r\left(\tau \right)·e^{-1j·2\pi ·f·\tau }·d\tau$

Or with a division into temporal blocks:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\stackrel{o}{s}}_k\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·w_k\left(\tau -t\right)·e^{-1j·2\pi ·f·\tau }·d\tau$ $\forall t\in T_{\downarrow },\forall f\in F_{+},\stackrel{o}{s}\left(t,f\right)={\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\stackrel{o}{s}}_k\left(t,f\right)}{w_k^{\left(1\right)}}\right)$

The expressions of the Fourier transforms provided in embodiments 41 to 44 above form applicable formulae, producing efficient estimators, notably in terms of their formulation with temporal weighting, averaging over several temporal intervals, overlaps between successive intervals, a more general expression of which can be provided with:

$\overbrace{s_k}\left(t,f\right)=\frac{1}{{\mathrm{dt}}_k}·{\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·w_k\left(\tau \right)·e^{-1\mu ·2\pi ·f·\tau }·d\tau$

with w_k being the temporal weighting, and u being defined such that 1μ²=−1. Finally,

$\overbrace{s}\left(t,f\right)={\mathrm{moy}}_k\left(\overbrace{s_k}\left(t,f\right)\right)$

The index k applied to w_k(above) or p_k(below) means that the weighting can differ from one interval to another.

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

A more general expression is:

$\overbrace{s_k}\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}}^{t+{\mathrm{dt}}_k^{+}}r\left(\tau \right)·p_k\left(\tau ,f\right)·d\tau$

Even more generally, the time division and the weighting can be varied as a function of the time-frequency cell (t, f):

$\left(t,f\right)={\int }_{t-{\mathrm{dt}}_k^{-}\left(t,f\right)}^{t+{\mathrm{dt}}_k^{+}\left(t,f\right)}r\left(\tau \right)·p_k\left(t,f\right)\left(\tau \right)·d\tau$

The trajectory (orbit) is obtained by inverse transformation, of the following type:

$s\left(t,f\right)\left(\tau \right)={\int }_{f-{\mathrm{df}}^{-}\left(t,f\right)}^{f+{\mathrm{df}}^{+}\left(t,f\right)}\left(\overbrace{s}\left(t,v\right)·q\left(t,f\right)\left(\tau ,v\right)+\overbrace{s}\left(t,-v\right)·q\left(t,f\right)\left(\tau ,-v\right)\right)·\mathrm{dv}$

With the frequency extent of the cell (t, f) then being df=df⁺+df⁻, the temporal extent of the cell (t, f) being dt=dt⁺+dt⁻.

From each spectral element s(t, f), at a fixed date t and frequency f, it is possible to deduce an equivalent ellipse from the instantaneous orbit. This can be referred to as local in the temporal and frequency sense, in other words, at the date t and frequency f.

Therefore, a step 50 is carried out of determining the equivalent ellipse for each spectral element. Irrespective of 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 s_xo, s_yo, and Ø_x, Ø_y being the instantaneous amplitudes and phases:

$s_x\left(t,f\right)=s_{x_o}\left(t,f\right)·\cos \left({\varnothing }_x\left(t,f\right)\right)$ $s_y\left(t,f\right)=s_{y_o}\left(t,f\right)·\cos \left({\varnothing }_y\left(t,f\right)\right)$

For a fixed frequency f, covering the time t allows the trajectory to be obtained at this frequency.

This instantaneous orbit assumes the form of an equivalent ellipse, provided that s_x(t, f), s_y(t, f)), {dot over (Ø)}_x (t, f), {dot over (Ø)}_y(t, f) vary slowly as a function of time, for example, during a complete revolution at the frequency f, and in any case for the duration dt of the temporal interval defining the spectral element.

In this case, the date t of the considered spectral element can be distinguished from the time t that evolves in this spectral element and, by introducing the frequency f into the expression of the phases that are assumed to be quasi-linear, the equivalent ellipse is obtained in the form of a temporal parametric equation for each component x and y:

${\stackrel{\prime }{s}}_x\left(t,f\right)\left(\tau \right)=a_x\left(t,f\right)·\cos \left(2·\pi ·f·\tau +{\phi }_x\left(t,f\right)\right)$ ${\stackrel{\prime }{s}}_y\left(t,f\right)\left(\tau \right)=a_y\left(t,f\right)·\cos \left(2·\pi ·f·\tau +{\phi }_y\left(t,f\right)\right)$

These parametric equations also can be written as follows, for example:

${\stackrel{\prime }{s}}_x\left(t,f\right)\left(\tau \right)=a_x\left(t,f\right)·\left(\cos \left(2·\pi ·f·\tau \right)·\cos \left({\phi }_x\left(t,f\right)\right)-\sin \left(2·\pi ·f·\tau \right)·\sin \left({\phi }_x\left(t,f\right)\right)\right)$ ${\stackrel{\prime }{s}}_y\left(t,f\right)\left(\tau \right)=a_y\left(t,f\right)·\left(\cos \left(2·\pi ·f·\tau \right)·\cos \left({\phi }_y\left(t,f\right)\right)-\sin \left(2·\pi ·f·\tau \right)·\sin \left({\phi }_y\left(t,f\right)\right)\right)$

In order to obtain these parametric equations, notably in the case of a frequency decomposition using a filter bank, the parameters a_x(t, f), φ_x(t, f), a_y(t, f), φ_y(t, f) of the equivalent ellipse can be obtained by estimates, such as direct or quadratic averaging or by RMS value, from the instantaneous quantities (τ, f) pour τ∈[t−dt⁻, t+dt⁺].

In the case of a frequency decomposition using a Fourier transform, the parametric equation of the equivalent ellipse is obtained using an inverse Fourier transform of the frequency element:

${\stackrel{\prime }{s}}_x\left(t,f\right)\left(\tau \right)=2·\mathrm{Re}\left({\hat{s}}_x\left(t,f\right)·e^{1i\text{.2}\pi .f.\tau }\right)$ ${\stackrel{\prime }{s}}_y\left(t,f\right)\left(\tau \right)=2·\mathrm{Re}\left({\hat{s}}_y\left(t,f\right)·e^{1i\text{.2}\pi .f.\tau }\right)$

By forming the associated complex quantity ś(t, f)={acute over (x)}_x(t, f)+1i. ś_y(t, f), in other words, the bivariate representation, it is possible to highlight the decomposition of a parametric equation of the equivalent ellipse in two circles in opposite directions:

$\stackrel{\prime }{s}\left(t,f\right)\left(\tau \right)=a_{-}\left(t,f\right)·e^{1i.{\theta }_{+}\left(t,f\right)}·e^{1i\text{.2}.\pi ,f.\tau }+a_{+}\left(t,f\right)·e^{1i.{\theta }_{+}\left(t,f\right)}·e^{1i\text{.2}.\pi ,f.\tau }$

With the following parameters obtained following a computation step 59:

$a_{+}=\frac{a}{2}·\sqrt{a_x^2+a_y^2+2·a_x·a_y·\sin \left({\phi }_x-{\phi }_y\right)}$ $a_{-}=\frac{a}{2}·\sqrt{a_x^2+a_y^2+2·a_x·a_y·\sin \left({\phi }_y-{\phi }_x\right)}$ ${\theta }_{+}=\arctan 2\left(a_x·\sin \left({\phi }_x\right)+a_y·\cos \left({\phi }_y\right),a_x·\cos \left({\phi }_x\right)-a_y·\sin \left({\phi }_y\right)\right)$ ${\theta }_{-}=\arctan 2\left(a_x·\sin \left({\phi }_x\right)-a_y·\cos \left({\phi }_y\right),a_x·\cos \left({\phi }_x\right)+a_y·\sin \left({\phi }_y\right)\right)$

The characterisation of the equivalent ellipse is straightforward using the two-sided (“full spectrum”) Fourier transform, which can be directly written in polar form:

$\forall f>0,\text{}\begin{array}{cc}\tilde{s}\left(t,f\right)& =a_{+}\left(t,f\right)·e^{1i.{\theta }_{+}\left(t,f\right)}\\ \tilde{s}\left(t,-f\right)& =a_{-}\left(t,f\right)·e^{-1i.{\theta }_{-}\left(t,f\right)}\end{array}$

In general, an equivalent ellipse can be determined and plotted by obtaining its Euler parameters.

FIG. 2 shows an example of an ellipse and its Euler parameters a(t, f), θ(t, f), χ(t, f), φ(t, f), which fully define an equivalent ellipse, with:

$a\ge 0$ $\theta \in \left[-\pi /2,\pi /2\right]$ $\chi \in \left[-\pi /4,\pi /4\right]$ $\phi \in \left[-\pi ,\pi \right]$

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 provides the direction of rotation of the equivalent ellipse.

The parametric equation of the trajectory of the equivalent ellipse is more practical in complex form, since:

• • an equation of a circle is notably expressed as cos and sin of the same phase;
  • the axes are dilated in order to obtain an ellipse;
  • a rotation is applied thereto in order to obtain the desired direction, thus:

$\begin{array}{cc}\stackrel{\prime }{s}\left(t,f\right)\left(\tau \right)=& a(t,f·e^{1i.\theta \left(t,f\right)}·(\cos \left(\chi \left(t,f\right)\right)·\cos \left(2·\pi ·f·\tau +\phi \left(t,f\right)\right)\end{array}+1i·\sin \left(\chi \left(t,f\right)\right)·\sin \left(2·\pi ·f·\tau +\phi \left(t,f\right)\right))$

Direction is understood to mean the angle θ of the equivalent ellipse between the major axis and the x-axis. Orientation is understood to mean the direction of rotation describing the trajectory of said ellipse.

This characterisation of the equivalent ellipse is natural with the time-frequency decomposition using a quaternionic Fourier transform, which can be directly placed into Euler form:

$\stackrel{o}{s}\left(t,f\right)=❘\stackrel{o}{s}\left(t,f\right)❘·e^{1i.\theta \left(t,f\right)}·e^{-1k.\chi \left(t,f\right)}·e^{1j.\phi \left(t,f\right)}$

This parametric equation has the particular advantage of rejecting the temporal parametrisation towards a single exponential term e^1j.2.π.f.τ when the inverse Fourier transform is applied to the frequency element in order to obtain the temporal parametric equation of the trajectory of the equivalent ellipse:

$\stackrel{\prime }{s}\left(t,f\right)\left(\tau \right)\approx ❘\stackrel{o}{s}\left(t,f\right)❘={\mathrm{Proj}}_{{\subset }_{1i}}\left(e^{1i.\theta \left(t,f\right)}·e^{-1k.\chi \left(t,f\right)}·e^{1j.\phi \left(t,f\right)}·e^{1j\text{.2}.\pi .f.\tau }\right)$

with Proj_C1i being a projection into the set of complex numbers C_1i.

More specifically, the step 50 of determining the equivalent ellipse includes a step of computing the Euler parameters of said equivalent ellipse.

Starting from the time-frequency decomposition using a filter bank, a step 51 of computing a_x and a_y is obtained such that:

$\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},& \begin{array}{cc}{a}_x\left(t,f\right)& =\sqrt{\frac{2}{\mathrm{dt}}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_x^2\left(\tau ,f\right)·d\tau }\\ a_y\left(t,f\right)& =\sqrt{\frac{2}{\mathrm{dt}}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_y^2\left(\tau ,f\right)·d\tau }\end{array}\end{array}$

• • with T_∇ being a sub-set of T.

Starting from the time-frequency decomposition using a short-term single-sided Fourier transform, a step 52 of computing a_x, a_y, φ_x and φ_y is obtained such that:

$$\begin{gathered} \begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},& \end{array}\begin{array}{c}\begin{array}{cc}{a}_x\left(t,f\right)=2·❘{\tilde{s}}_x\left(t,f\right)❘& {\phi }_x\left(t,f\right)=\mathrm{Arg}\left({\tilde{s}}_x\left(t,f\right)\right)\end{array} \\ \\ \begin{array}{cc}{a}_y\left(t,f\right)=2·❘{\tilde{s}}_y\left(t,f\right)❘& {\phi }_y\left(t,f\right)=\mathrm{Arg}\left({\tilde{s}}_y\left(t,f\right)\right)\end{array}\end{array} \end{gathered}$$

Starting from the time-frequency decomposition using a short-term two-sided Fourier transform, a step 53 of computing a₊, a₋, θ₊ and θ₋ is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{a}_{+}\left(t,f\right)=❘\tilde{s}\left(t,f\right)❘{\theta }_{+}\left(t,f\right)=\mathrm{Arg}\left(\tilde{s}\left(t,f\right)\right)\\ a_{-}\left(t,f\right)=❘\tilde{s}\left(t,f\right)❘{\theta }_{-}\left(t,f\right)=-\mathrm{Arg}\left(\tilde{s}\left(t,-f\right)\right)\end{array}$

An additional computation step 54 or 55 yields the following Euler parameters:

$a=\sqrt{2}·\sqrt{a_{+}^2+a_{-}^2}a=\sqrt{a_x^2+a_y^2}$ $\phi =\frac{{\theta }_{+}+{\theta }_{-}}{2}\mathrm{ou}\phi =\frac{{\phi }_x+{\phi }_y}{2}$ $\theta =\frac{{\theta }_{+}+{\theta }_{-}}{2}\theta =\frac{1}{2}·\mathrm{arc}\tan 2\left(2·a_x·a_y·\cos \left({\phi }_x-{\phi }_y\right),a_x^2-a_y^2\right)$ $\chi =\arctan 2\left(a_{+}-a_{-},a_{+}+a_{-}\right)\chi =\frac{1}{2}·\mathrm{arc}\sin \left(\mathrm{ratio}*\left(2·a_x·a_y,a_x^2+a_y^2\right)·\sin \left({\phi }_x-{\phi }_y\right)\right)$

An inverse computation step 56 is also possible, such that:

$a_{+}=\frac{a}{2}·\left(\cos \left(\chi \right)+\sin \left(\chi \right)\right)$ $a_{-}=\frac{a}{2}·\left(\cos \left(\chi \right)-\sin \left(\chi \right)\right)$ ${\theta }_{+}=\phi +\theta$ ${\theta }_{-}=\phi -\theta$

The method according to the present disclosure optionally comprises a step 60 of computing the spectral density of the time-frequency decomposition. The spectral densities of each component and the cross-spectral densities between components are used to estimate the energy distribution between the x and y components.

Starting from the time-frequency decomposition using a filter bank, a step 61 of computing the spectral density S is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}{S}_x\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_x^2\left(\tau ,f\right)·d\tau \\ S_y\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_y^2\left(\tau ,f\right)·d\tau \\ S_{\mathrm{xy}}\left(t,f\right)=\frac{1}{\mathrm{dt}}·{\int }_{t-{\mathrm{dt}}^{-}}^{t+{\mathrm{dt}}^{+}}{s}_x\left(\tau ,f\right)·s_y\left(\tau ,f\right)·d\tau \end{array}$

Starting from the time-frequency decomposition using a short-term single-sided Fourier transform, a step 62 of computing the spectral density Ŝ is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{\widehat{S}}_x\left(t,f\right)=2·\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘{\hat{s}}_{x_k}\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)\\ {\widehat{S}}_y\left(t,f\right)=2·\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘{\hat{s}}_{y_k}\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)\\ {\widehat{S}}_{\mathrm{xy}}\left(t,f\right)=2·\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\hat{s}}_{x_k}\left(t,f\right)·\stackrel{\_}{{\hat{s}}_{y_k}\left(t,f\right)}}{w_k^{\left(2\right)}}\right)\end{array}$

Starting from the time-frequency decomposition using a short-term two-sided Fourier transform, a step 63 of computing the spectral density S is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F^{*},\tilde{S}\left(t,f\right)=\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘{\hat{s}}_k\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)$

The following also will be used:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\left(t,f\right)=\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{\tilde{s}}_k\left(t,f\right)·{\tilde{s}}_k\left(t,-f\right)}{w_k^{\left(2\right)}}\right)$

Starting from the time-frequency decomposition using a short-term quaternionic Fourier transform, a step 64 of computing the spectral density S̊ is obtained such that: ∀t∈T_∇, ∀f∈F₊,

$\stackrel{o}{S}\left(t,f\right)=\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{{❘\stackrel{o}{s_k}\left(t,f\right)❘}^2}{w_k^{\left(2\right)}}\right)+\mathrm{df}·{\mathrm{moy}}_{k=\left[1,p\right]}\left(\frac{\stackrel{o}{s_k}\left(t,f\right)·1j·\stackrel{\_}{\stackrel{o}{s_k}\left(t,f\right)}}{w_k^{\left(2\right)}}\right)$

Based on these spectral densities S, Ŝ, Ś, and S̊, a step 70 is carried out of computing the Stokes parameters S₀, S₁, S₂, and S₃ of at least one equivalent ellipse.

The Stokes parameters are four real numbers and are conventionally used to describe the polarisation state of an electromagnetic wave.

S₀ is commonly the total intensity of an electromagnetic wave and in this case by extension represents the total intensity of the bivariate signal, S₀>0.

S₁ and S₂ provide the intensity of rectilinear polarisation √{square root over (S₁²+S₂²)}, as well as the direction of the polarisation of the ellipse

$\theta =\frac{1}{2}.$

atan2(S₂,S₁).

S₃ provides the circular polarisation intensity |S₃|, with its sign providing the direction of rotation.

The reduced Stokes parameters are also defined:

$\mathrm{si}{S}_0\ne 0,\forall i\in \left\{1,2,3\right\}{s}_i\stackrel{\mathrm{def}}{=}\frac{s_i}{s_0}$

Starting from the spectral density obtained for the time-frequency decomposition using a filter bank, a step 71 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{S}_0\left(t,f\right)=2·\left(S_x\left(t,f\right)+S_y\left(t,f\right)\right)\\ S_1\left(t,f\right)=2·\left(S_x\left(t,f\right)-S_y\left(t,f\right)\right)\end{array}$

Starting from the spectral density obtained for the time-frequency decomposition using a short-term single-sided Fourier transform, a step 72 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{S}_0\left(t,f\right)=2·\left({\hat{S}}_x\left(t,f\right)+{\hat{S}}_y\left(t,f\right)\right)\\ S_1\left(t,f\right)=2·\left({\hat{S}}_x\left(t,f\right)+{\hat{S}}_y\left(t,f\right)\right)\\ S_2\left(t,f\right)=4·\mathrm{Re}\left({\hat{S}}_{\mathrm{xy}}\left(t,f\right)\right)\\ S_3\left(t,f\right)=4·\mathrm{Im}\left({\hat{S}}_{\mathrm{xy}}\left(t,f\right)\right)\end{array}$

Starting from the spectral density obtained for the time-frequency decomposition using a short-term two-sided Fourier transform, a step 73 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\begin{array}{c}{S}_0\left(t,f\right)=2·\left(\tilde{S}\left(t,f\right)+\tilde{S}\left(t,-f\right)\right)\\ S_1\left(t,f\right)=4·\mathrm{Re}\left(\left(t,f\right)\right)\\ S_2\left(t,f\right)=4·\mathrm{Im}\left(\left(t,f\right)\right)\\ S_3\left(t,f\right)=2·(\tilde{S}\left(t,f\right)-\tilde{S}\left(t,-f\right)\end{array}$

Starting from the spectral density obtained for the time-frequency decomposition using a short-term quaternionic Fourier transform, a step 74 of computing the Stokes parameters of at least one equivalent ellipse is obtained such that:

$t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}{S}_0\left(t,f\right)=\mathrm{Re}\left(\stackrel{o}{S}\left(t,f\right)\right)\\ S_1\left(t,f\right)={\mathrm{Im}}_{1j}\left(\stackrel{o}{S}\left(t,f\right)\right)\\ S_2\left(t,f\right)={\mathrm{Im}}_{1k}\left(\stackrel{o}{S}\left(t,f\right)\right)\\ S_3\left(t,f\right)={\mathrm{Im}}_{1i}\left(\stackrel{o}{S}\left(t,f\right)\right)\end{array}$

As an alternative embodiment, the Stokes parameters S₁, S₂, and S₃ are directly computed using the time-frequency decomposition using a quaternionic Fourier transform:

$\forall t\in T_{\downarrow },\forall f\in F_{+},\begin{array}{c}\stackrel{o}{\sigma }\left(t,f\right)=\stackrel{o}{s}\left(t,f\right)·1j·\stackrel{\_}{\stackrel{o}{s}\left(t,f\right)}\\ S_1\left(t,f\right)={\mathrm{Im}}_{1j}\left(\stackrel{o}{\sigma }\left(t,f\right)\right)\\ S_2\left(t,f\right)={\mathrm{Im}}_{1k}\left(\stackrel{o}{\sigma }\left(t,f\right)\right)\\ S_3\left(t,f\right)={\mathrm{Im}}_{1i}\left(\stackrel{o}{\sigma }\left(t,f\right)\right)\end{array}$

Based on these Stokes parameters, a step 75 is carried out of computing the polarisation rate Ø, as well as the aforementioned reduced Stokes parameters, as well as the polarised part of the intensity S_p.

$\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+},& \forall i\in \left\{1,2,3\right\}{s}_i\left(t,f\right)=\mathrm{ratio}*\left(S_i\left(t,f\right),S_0\left(t,f\right)\right)\\ & \varnothing \left(t,f\right)=\sqrt{s_1^2\left(t,f\right)+s_2^2\left(t,f\right)+s_3^2\left(t,f\right)}\\ & S_p\left(t,f\right)=\varnothing \left(t,f\right)·S_0\left(t,f\right)\end{array}$

The polarisation expressed by the polarisation rate Ø highlights the cyclo-stationarity of an equivalent ellipse, namely, the regularity of the trajectory of said ellipse relative to a period of revolution. The polarisation reflects the structuring of said trajectory relative to defined frequencies. The polarisation rate thus can be computed 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.

In addition, a computation step 57 can be implemented in order to determine the Euler parameters a, θ, and χ based on the polarisation spectral density S_p and the reduced Stokes parameters:

$\begin{array}{cc}\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},& a\left(t,f\right)=\sqrt{S_p\left(t,f\right)}\\ & \theta \left(t,f\right)=\frac{1}{2}·\arctan 2\left(s_2\left(t,f\right),s_1\left(t,f\right)\right)\\ & \chi \left(t,f\right)=\frac{1}{2}·\arcsin \left(\mathrm{ratio}*\left(s_3\left(t,f\right),S_p\left(t,f\right)\right)\right)\end{array}$

A fourth Euler parameter φ providing the phase that is the source of the trajectory also can be deduced from the following computation step 58:

$\forall t\in T_{\downarrow },\forall f\in F_{+}^{*},\varphi \left(t,f\right)=\mathrm{Arg}\left({\mathrm{Proj}}_{C_{1j}}\left(e^{1k·\chi \left(t,f\right)}·e^{-1i·\theta \left(t,f\right)}·\stackrel{o}{s}\left(t,f\right)\right)\right)$

Following the step 50 of determining an equivalent ellipse, notably by computing the Euler parameters of said equivalent ellipse, a step 80 of extracting at least one characterisation indicator of the rotating shaft from parameters of the equivalent ellipse, for example, Euler or Stokes parameters or parameters derived therefrom, is carried out.

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 initially carried out. In this embodiment, it is thus possible to carry out a step of comparing the characterisation indicators of the planes.

For example, the at least one characterisation indicator includes the previously described polarisation rate Ø, and/or a circularity indicator (linked to |χ|), and/or an intensity indicator and/or an amplitude indicator, and/or a direction of rotation indicator (linked to the sign of χ), and/or a direction indicator (linked to θ), and/or synchronisation indicator (linked to φ) of the equivalent ellipse.

Each elementary indicator helps to provide a practical and useful characteristic of the rotating system, for example, a high rate of polarisation and a high degree of circularity are probably due to the rotating shaft. Distinguishing this part of the signal, with the complementary part being treated as disturbances, can be useful, for example, for identifying said shaft (for modelling purposes) or for better controlling it in the case of a shaft equipped with magnetic bearings.

Optionally, the at least one characterisation indicator is normalised, in other words, its value ranges between −1 and 1.

The characterisation indicator is provided for one plane, it is then denoted A, or by a comparison between two planes, it is then denoted I.

The characterisation indicator is referred to as a criterion, denoted C, when it is constructed from any characteristic quantity re-evaluated at each instant t. The characterisation indicator is referred to as stationarity, denoted S, when it is based on the stationarity of a quantity over time (maximum when the variation of said variable over time is zero).

Thus, the polarisation rate indicator is denoted:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{polarisation}}\left(t,f\right)=\Phi \left(t,f\right)$

For one plane, the intensity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{intensité}}\left(t,f\right)={\mathrm{normalise}}_{\tau \in T_{\downarrow }v\in F_{+}}\left(S_0\left(t,f\right)\right)$

With the normalisation function:

$\forall t\in \mathrm{TT},\forall f\in \mathrm{FF},{\mathrm{normalise}}_{\tau \in \mathrm{TT},v\in \mathrm{FF}}\left(x\left(t,f\right)\right)=\frac{x\left(t,f\right)}{{\max }_{\tau \in \mathrm{TT},v\in \mathrm{FF}}^{*}\left(❘x\left(\tau ,v\right)❘\right)}$

TT and FF denote sets of generic dates and frequencies, taking their values from the set of real values R.

For one plane, the amplitude indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{amplitude}}\left(t,f\right)={\mathrm{normalise}}_{\tau \in T_{\downarrow }v\in F_{+}}\left(a\left(t,f\right)\right)$

For one plane, the sign indicator is −1 or 1 depending on the sign of x:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{signe}}\left(t,f\right)=\mathrm{signe}\left(\chi \left(t,f\right)\right)$

For one plane, the circularity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AC}}_{\mathrm{circularité}}\left(t,f\right)=\frac{❘\chi \left(t,f\right)❘}{\pi /4}$

For one plane, the sign stationarity indicator, in other words, the direction of rotation, in other words, the orientation, is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{signe}}\left(t,f\right)=\mathrm{autostat}\left({\mathrm{AC}}_{\mathrm{signe}}\left(t,f\right),\mathrm{tol}=0\right)$

With the stationarity function:

$\forall t\in \mathrm{TT},\begin{array}{cc}\mathrm{pour}\mathrm{tol}>0:& \mathrm{autostat}\left(x\left(t\right),\mathrm{tol}\right)=e^{-\frac{❘\mathrm{dx}\left(f\right)❘}{\mathrm{tol}}}\\ \mathrm{pour}\mathrm{tol}=0:& \mathrm{autostat}\left(x\left(t\right),0\right)=\left(\mathrm{dx}\left(t\right)==0\right)\end{array}$

For one plane, the circularity stationarity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{circularité}}\left(t,f\right)=\mathrm{autostat}\left({\mathrm{AC}}_{\mathrm{circularité}}\left(t,f\right),\mathrm{tol}=0.1\right)$

For one plane, the direction stationarity indicator, in other words, the geometric inclination of the ellipse, is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{direction}}\left(t,f\right)=\mathrm{autostat}\left(\theta \left(t,f\right),\mathrm{tol}=10·\pi /180\right)$

For one plane, the synchronisation stationarity indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{AS}}_{\mathrm{synchro}}\left(t,f\right)=\mathrm{autostat}\left(\phi \left(t,f\right),\mathrm{tol}=10·\pi /180\right)$

For two planes, each parameter or indicator is redefined for an additional plane, by denoting these new parameters with an apostrophe, for example, a′.

The amplitude comparative characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{amplitude}}\left(t,f\right)={\mathrm{normalise}}_{\tau \in T_{\downarrow }v\in F_{+}}\left(\sqrt{a\left(t,f\right)·a^{\prime }\left(t,f\right)}\right)$

The sign comparative characterisation indicator is a Boolean between 0 and 1 depending on the result of the following test:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{signe}\left(t,f\right)=\left(A{C}_{signe}\left(t,f\right)==A{C}_{signe}^{\prime }\left(t,f\right)\right)$

The circularity comparative characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{circularit}é}\left(t,f\right)=1-\mathrm{reldif}\left(A{C}_{\mathrm{circularit}é}\left(t,f\right),A{C}_{\mathrm{circularit}é}^{\prime }\left(t,f\right)\right)$

With the relative error function such that:

$\mathrm{reldif}\left(x,y\right)={\mathrm{ratio}}^{*}\left(2.❘x-y❘,❘x❘+❘y❘\right)$

And with the division function with exception such that:

${\mathrm{ratio}}^{*}\left(x,y\right)=\begin{array}{c}x/y\mathrm{si}y\ne 0\\ 0\mathrm{sinon}\end{array}$

The direction comparative characterisation indicator, in other words, the geometric inclination of the ellipse, is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{\mathrm{direction}}\left(t,f\right)=1-\mathrm{reldif}\left(\theta \left(t,f\right),{\theta }^{\prime }\left(t,f\right)\right)$

The synchronisation comparative characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IC}}_{synchro}\left(t,f\right)=1-reldif\left(\varphi \left(t,f\right),{\phi }^{\prime }\left(t,f\right)\right)$

The stationarity and comparative circularity characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IS}}_{\mathrm{circularit}é}\left(t,f\right)=autostat\left(A{C}_{\mathrm{circularit}é}\left(t,f\right)-A{C}_{\mathrm{circularit}é}^{\prime }\left(t,f\right),\mathrm{tol}=0.1\right)$

The stationarity and comparative direction characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IS}}_{di\mathrm{rection}}\left(t,f\right)=autostat\left(\theta \left(t,f\right)-{\theta }^{\prime }\left(t,f\right),\mathrm{tol}=10.\pi /180\right)$

The stationarity and comparative synchronisation characterisation indicator is:

$\forall t\in T_{\downarrow },\forall f\in F_{+},{\mathrm{IS}}_{\mathrm{synchro}}\left(t,f\right)=autos\mathrm{tat}\left(\phi \left(t,f\right)-{\phi }^{\prime }\left(t,f\right),\mathrm{tol}=10.\pi /180\right)$

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

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

Optionally, a step 100 is carried out of displaying the at least one extracted characterisation indicator so that it can be read by an operator and so that they can understand potential anomalies in the rotating shaft.

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

The images are, for example, orbits, or time-frequency representations, or representations of phases according to the frequency, or coloured images, with the colour representing the intensity, the direction of rotation or another indicator.

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

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

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, with each cell being arranged according to a Cartesian coordinate representation. Composite or hybrid images also can 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, with the various quantities optionally being distinguished or supported by shapes or graphic symbols, for example, arrows for representing a vector or a direction.

A step 102 is then carried out of creating a collection of images from the images generated in the generation step 101, with the images being able to be ordered in a selected and/or predetermined order.

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

In a particularly advantageous embodiment, the classification tool is a neural network, preferably a convolutional neural network, which is particularly effective for classifying images. Even more preferably, the convolutional neural network has dimensions that are greater than or equal to the dimensions of the images of the collection of images.

For example, if the collection of images contains only one two-dimensional image, then this is equivalent to classifying only one two-dimensional image and a two-dimensional convolutional neural network can be used.

In another example, if the collection of images is used to represent a three-dimensional image, a three-dimensional convolutional neural network can be used.

In a particular embodiment, each image from the collection is made up of one or more pixels grouped together as a function of the selected resolution, representing the position of the rotor in a fixed orthogonal plane relative 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 rotor orbit 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 collection of images contains all the images described above but for all the planes orthogonal to the longitudinal axis of the rotor ordered along this axis and for all the instants ordered according to the increasing time, then the collection of images would describe a four-dimensional image, which could be classified by a four-dimensional convolutional neural network.

The collections of images supplying the neural network allow a learning and training base to be formed for the neural network, the classification function of which improves as it receives collections of images.

The classification step 104 optionally allows alarms to be activated for controlling or protecting the rotating shaft and/or its environment, notably when significant anomalies are detected during said classification step 104.

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

In a particular embodiment, all the steps of the method are computer-implemented, with the possibility of real time execution, notably in a magnetic bearing control cabinet.

Some terms and functions used in the computations described above are redefined below:

The set of strictly positive frequencies is: F₊*={f_j, j=[1, n]}

The set of strictly negative frequencies is: F₋*=−F₊

The following is also defined: F₊={0}∪F;

The following is also defined: F₋=F₋*∪ {0}

The complete set of frequencies is: F=F₋∪F₊=F₋*∪{0} u F;

The generic stationarity function is:

$\forall t\in \mathrm{TT},\begin{array}{cc}\mathrm{pour}\mathrm{tol}>0:& \mathrm{autostat}\left(x\left(t\right),\mathrm{tol}\right)=e^{-\frac{❘\mathrm{dx}\left(t\right)❘}{\mathrm{tol}}}\\ \mathrm{pour}\mathrm{tol}=0:& \mathrm{autostat}\left(x\left(t\right),0\right)=\left(\mathrm{dx}\left(t\right)==0\right)\end{array}$

The maximum function with exception is:

${\max }^{*}\left(x\right)=\begin{array}{c}\max \left(x\right)\mathrm{si}\max \left(x\right)\ne 0\\ 1\mathrm{sinon}\end{array}$

The normalisation function is:

$\forall t\in TT,\forall f\in \mathrm{FF},{\mathrm{normalise}}_{\tau \in \mathrm{TT},\nu \in \mathrm{FF}}\left(x\left(t,f\right)\right)=\frac{x\left(t,f\right)}{{\max }_{\tau \in \mathrm{TT},\nu \in \mathrm{FF}}^{*}\left(❘x\left(\tau ,\nu \right)❘\right)}$

The division function with exception is:

${\mathrm{ratio}}^{*}\left(x,y\right)=\begin{array}{c}x/y\mathrm{si}y\ne 0\\ 0\mathrm{sinon}\end{array}$

The relative error function is:

$\mathrm{reldif}\left(x,y\right)={\mathrm{ratio}}^{*}\left(2.❘x-y❘,❘x❘+❘y❘\right)$

The definition of the arctan 2 function is also provided by way of a reminder, with a formulation equivalent to the argument of the complex number x+1μ. y:

$Arg\left(x+1\mu .y\right)=\mathrm{arc}\tan 2\left(y,x\right)$

Claims

1. A method for characterising the behaviour of a rotating shaft, the method comprising the following steps: acquiring at least two input signals over a predetermined period, the input signals comprising measurements representing the movement of the rotating shaft, the measurements being non-collinear to each other;determining two components of the input signals, the components being orthogonal to each other;merging the two components into a bivariate signal;time-frequency decomposing the bivariate signal into spectral elements;determining an equivalent ellipse for each spectral element; andextracting at least one characterisation indicator of the rotating shaft based on parameters of the equivalent ellipse.

2. The method according to claim 1, wherein the input signals include measurements of the position, or speed, or acceleration, or force, or magnetic field, or magnetic induction of the rotating shaft.

3. The method according to claim 1, wherein the time-frequency decomposition step is carried out using a Fourier, or wavelet, transform, or by filtering, or using a quaternionic Fourier transform.

4. The method according to claim 1, wherein the step of determining an equivalent ellipse includes computing Euler parameters of said equivalent ellipse.

5. The method according to claim 1, comprising a step of computing the spectral density of the time-frequency decomposition, and a step of computing Stokes parameters based on said spectral density.

6. The method according to claim 1, wherein the at least one characterisation indicator includes a polarisation rate, and/or an indicator of circularity and/or amplitude and/or direction of rotation and/or direction and/or synchronisation of an equivalent ellipse.

7. The method according to claim 1, further comprising a step of introducing the at least one characterisation indicator into a neural network in order to classify said indicator as a function of the rotation of the shaft and to link it to a previously listed anomaly.

8. The method according to claim 1, further comprising a step of displaying the at least one extracted characterisation indicator so that it can be read by an operator.

9. The method according to claim 1, comprising 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 based on the images generated in the generation step, a step of providing a classification tool with said collection of images, and a step of classifying, using the classification tool, the images as a function of the rotation of the rotating shaft and/or of the anomalies present in the collection of images.

10. The method according to claim 9, wherein the classification tool is a neural network, preferably a convolutional neural network, even more preferably a convolutional neural network with dimensions that are greater than or equal to the dimensions of the images of the collection of images.

11. The method according to claim 1, comprising a step of displaying an image of at least one equivalent ellipse.

12. The method according to claim 1, wherein the steps of claim 1 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 characterisation indicators of the planes and/or creating new characterisation indicators common to the two planes.

13. The method according to claim 2, wherein the time-frequency decomposition step is carried out using a Fourier, or wavelet, transform, or by filtering, or using a quaternionic Fourier transform.

14. The method according to claim 13, wherein the step of determining an equivalent ellipse includes computing Euler parameters of said equivalent ellipse.

15. The method according to claim 14, comprising a step of computing the spectral density of the time-frequency decomposition, and a step of computing Stokes parameters based on said spectral density.

16. The method according to claim 15, wherein the at least one characterisation indicator includes a polarisation rate, and/or an indicator of circularity and/or amplitude and/or direction of rotation and/or direction and/or synchronisation of an equivalent ellipse.

17. The method according to claim 16, further comprising a step of introducing the at least one characterisation indicator into a neural network in order to classify said indicator as a function of the rotation of the shaft and to link it to a previously listed anomaly.

18. The method according to claim 17, further comprising a step of displaying the at least one extracted characterisation indicator so that it can be read by an operator.

19. The method according to claim 18, comprising 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 based on the images generated in the generation step, a step of providing a classification tool with said collection of images, and a step of classifying, using the classification tool, the images as a function of the rotation of the rotating shaft and/or of the anomalies present in the collection of images.

20. The method according to claim 19, wherein the classification tool is a neural network, preferably a convolutional neural network, even more preferably a convolutional neural network with dimensions that are greater than or equal to the dimensions of the images of the collection of images.