Method and device for detecting a substantially invariant rotation axis

[0001] This application claims priority to French Patent Application No. 0754693, filed Apr. 25, 2007, which is incorporated by reference herein.

[0002] The present invention concerns a method for detecting a substantially invariant rotation axis of the motion of a solid. It also concerns a method for estimating the motion of a solid and a method for calibrating an inertial or magnetic sensor with three sensitive axes. In a correlated way, the present invention concerns a device for detecting a substantially invariant rotation axis of the motion of a solid. Generally, the present invention concerns the field of motion capture. The fields of application are varied and cover in particular the biomedical field, for the observation of persons, and the sports field, for the analysis of the motion of a sportsperson or their equipment (racket, ball, etc.). It also applies to the field of automobiles or robotics, as well as to virtual reality, and generally to any application including a mobile body the motion whereof is to be determined or observed.

[0003] There are known in particular inertial systems consisting of one or more sensors selected in particular from accelerometers, magnetometers or rate gyros, having one or more sensitive axes.

[0004] These inertial systems have the advantage of being autonomous and of not requiring the environment in which the motion of a mobile object is observed to be fitted out with equipment beforehand.

[0005] There is known in particular a device for sensing rotation of a solid as described in the document FR 2 838 185.

[0006] The principle of observing a mobile body in motion is to seek magnitudes representing the displacement and orientation of the mobile body, or its speed and rotation speed, or its acceleration and angular acceleration.

[0007] In practice, the sensors are placed on the observed mobile body and the aim is to determine a motion that potentially has six degrees of freedom, i.e. three degrees of freedom corresponding to the orientation of the solid in space and three degrees of freedom corresponding to the position of the mobile body.

[0008] Solving systems of equations to determine magnitudes representing the motion of a solid is generally complex and necessitates solution systems in which certain constraints are known in advance.

[0009] An aim of the present invention is to overcome the drawbacks of earlier motion estimating systems and to propose a method that is particularly suitable for processing the motion of a mobile body having a substantially invariant rotation axis.

[0010] A first aspect of the present invention is directed to a method for detecting a substantially invariant rotation axis of a motion of a mobile body equipped with at least one inertial or magnetic sensor with three sensitive axes.

[0011] According to the invention, this detection method comprises the following steps:

[0012] acquisition of physical measurements with respect to the three sensitive axes of the sensor, the physical measurements comprising at least three samples at different times; estimation of a substantially invariant rotation axis in the physical measurements space; and identification of said estimated axis as the substantially invariant rotation axis of the motion.

[0013] The applicant has found that by determining the substantially invariant rotation axis among the measurements obtained at the output of the same sensor it is possible to determine the substantially constant rotation axis of the motion of the mobile body equipped with that sensor.

[0014] Estimating the rotation axis in the physical measurements space requires at least three different measurements at different times and with respect to the three sensitive axes of the sensor, to enable determination of the coordinates of the estimated rotation axis.

[0015] In practice, the detection method further comprises the following steps: calculation of an indicator of the variations of the rotation axis;

[0016] comparison of the value of the indicator with a predetermined threshold value; and

[0017] validation of the estimated rotation axis as the substantial invariant rotation axis of the motion if the value of the indicator is less than the predetermined threshold value which is advantageously equal to the value of the noise affecting the physical measurements from said sensor.

[0018] Accordingly, thanks to an indicator of variations of the estimated rotation axis, it is possible to obtain an indication of the variations of the rotation axis of the motion. If the value of this indicator is sufficiently low, i.e. of the same order as the noise affecting the physical measurements from the observed sensor, the motion can be considered to correspond faithfully to rotation about a substantially invariant axis. Thus it is possible to validate the estimated rotation axis in the physical measurements space as the substantially invariant rotation axis of the motion of the mobile body.

[0019] According to one practical embodiment of the invention, when the mobile body is equipped at least with a magnetometer having three sensitive axes and an accelerometer having three sensitive axes, the detection method comprises the following steps:

[0020] acquisition of the physical measurements with respect to the three sensitive axes of the magnetometer and the accelerometer, respectively, the physical measurements comprising at least three samples at different times;

[0021] estimation of the substantially invariant rotation axes in the said physical measurements spaces of said magnetometer and said accelerometer, respectively; and

[0022] validation as the substantially invariant rotation axis of said motion of a linear combination of the estimated rotation axes for said physical measurements from the magnetometer and from the accelerometer, the coefficients of the combination being not strictly between 0 and 1, their sum being equal to 1, and being a function of the values of the indicators of the variations of the rotation axes of the magnetometer and/or the accelerometer.

[0023] In particular, if one of the coefficients of the combination is equal to zero, it is possible to use either sensor (magnetometer or accelerometer) providing the most reliable estimate of the substantially invariant rotation axis.

[0024] In particular, if the rotation occurs about an axis of the physical field considered (the gravitational field in the case of accelerometers, the terrestrial magnetic field in the case of magnetometers), it is possible to determine the invariant rotation axis using the other sensor if the indicator of the variations of the rotation axis estimated by means of the physical measurements from the sensor associated with the invariant physical field will necessarily be inadequate, i.e. very much greater than the noise affecting the physical measurements obtained by that sensor.

[0025] In one practical embodiment of the invention, the detection method further includes a step of transformation of the coordinates of the estimated rotation axis in the measurement frame of reference of the sensor into coordinates in a geophysical frame of reference, such as an NED (North-East-Down) frame of reference.

[0026] A second aspect of the present invention concerns a method for estimating a motion of a mobile body equipped with at least one inertial or magnetic sensor with three sensitive axes, the motion comprising rotation about a substantially invariant axis.

[0027] According to the invention the estimation method comprises the following steps:

[0028] detection of a substantially invariant rotation axis by the detection method according to the invention;

[0029] determination of the rotation angle of the mobile body about the substantially invariant rotation axis.

[0030] If the sensor is a magnetometer or a goniometer and the mobile body is further equipped with an accelerometer, the estimation method further comprises a step of determining the acceleration of the mobile body.

[0031] Thanks to the method for the invention of detecting a substantially invariant rotation axis, the method for estimating the motion of a mobile body is facilitated if the number of unknowns to be determined is limited to the determination of the rotation angle, and where appropriate the acceleration of the mobile body.

[0032] Following detection in accordance with the invention of the invariant rotation axis, it is possible to calculate this motion in two steps, firstly seeking the rotation angle, and then calculating the acceleration of the mobile body.

[0033] A third aspect of the invention concerns a method for calibrating an inertial or magnetic sensor with three sensitive axes for determining a rotation matrix suitable for transforming coordinates expressed in a measurement frame of reference of the sensor into coordinates expressed in a predetermined frame of reference stationary with respect to a mobile body carrying the sensor.

[0034] This calculation method comprises the following steps:

[0035] acquisition of at least two series of physical measurements by the sensor respectively corresponding to the rotation of the mobile body with respect to at least two of the three axes of the predetermined frame of reference stationary with respect to the mobile body;

[0036] determination from each of the two series of coordinates of each of the rotation axes corresponding to one of said two of the three axes of the predetermined frame of reference stationary with respect to the mobile body in the frame of reference of the physical measurements from the sensor according to the detection method according to the invention; and

[0037] generation of a rotation matrix from the determined coordinates of said at least two of the three rotation axes.

[0038] This calibration method means that the measures obtained by the sensor can then be expressed in a predetermined frame of reference linked to the observed mobile body, such as a casing enclosing the sensor or an anatomical feature of a person on whom the sensor is placed.

[0039] Finally, a fourth aspect of the present invention concerns a device for detecting a substantially invariant rotation axis of a motion of a mobile body equipped with at least one inertial or magnetic sensor with three sensitive axes.

[0040] This detector device comprises means adapted to implement the detection method according to the invention.

[0041] This detection device has features and advantages analogous to those described hereinabove in relation to the detection method.

[0042] Other features and advantages of the invention will become apparent in the course of the following description.

[0051] A first embodiment of a method for detecting a substantially invariant rotation axis of a motion of a mobile body equipped with a sensor with three sensitive axes is described first with reference to FIG. 1.

[0052] This method can be applied to studying the motion of a mobile body in numerous applications where the motion of the mobile body corresponds to rotation about a fixed or substantially invariant axis.

[0053] There are many applications in which the motion of the observed mobile body corresponds to a rotation about an axis, generally also with no motion along this rotation axis, and so the motion is performed in a plane and can be considered a planar motion.

[0054] For example, such motion is observed in the case of a vehicle with zero roll and pitch. The rotation of the vehicle corresponds to a rotation about a fixed vertical axis, corresponding to a yaw motion. The vertical acceleration is generally zero and so the motion of the vehicle is a planar motion.

[0055] Similarly, motion of the members of a person (arms, legs, and the like) is observed in a plane if there is only one rotation axis of one joint and one angle that varies over time.

[0056] The motion of a wheel rolling without slipping on a plane and in a straight line also corresponds to motion about a rotation axis, perpendicular to the plane of the wheel.

[0057] The movement of a computer mouse in a plane about an axis orthogonal to that plane is also a planar motion that can be observed in the context of the present invention.

[0058] Generally speaking, the method for detecting the substantially invariant rotation axis of such motion is implemented by means of a detector device equipped with at least one sensor having three sensitive axes, for example an accelerometer, a magnetometer or a three-axis rate gyro.

[0059] In theory, the detection method consists in estimating a substantially invariant rotation axis in the space of the physical measurements obtained by the sensor with three sensitive axes.

[0060] In the embodiment illustrated in FIG. 1, it is considered that the substantially invariant rotation axis in the measurements space corresponds to the axis showing the least variation of the measurements for the magnetometer and the accelerometer and that showing the greatest variation for the rate gyro.

[0061] Each measurement at a time t obtained at the output of a sensor with three sensitive axes is considered as a point in the three-dimensional measurements space.

[0062] Thus the invariant axis is the axis with respect to which magnetometer or accelerometer measurements show the least variation on that axis, and in particular, for which the standard deviation of points projected onto that axis is the lowest.

[0063] In the case of a rate gyro, this will be the axis with the greatest variation, since the rate gyro measures the value of the rotation about each of its axes.

[0064] FIG. 2 is a simplified illustration in a two-dimensional space of measurements represented by a set of points.

[0065] Of course, this diagrammatic representation can be extended to a space with three dimensions.

[0066] The axis F showing the least variation and the axis G showing the greatest variation are shown in FIG. 2.

[0067] The arrow v in the direction of the axis F showing the least variation corresponds to the standard deviation of the measurements projected onto that axis and thus corresponds to an indicator of variations of the rotation axis. The lower this variance v, the more closely the observed motion corresponds to rotation about a fixed axis.

[0068] There will be described with reference to FIG. 1 a method for detecting a substantially invariant rotation axis enabling the axis showing the least variation to be determined from a principal component analysis.

[0069] In practice, a step E10 of acquiring physical measurements with respect to the three sensitive axes of the sensor is executed.

[0070] The physical measurements comprise N samples at different times with respect to the three sensitive axes of the sensor.

[0071] To enable the problem to be solved, N is greater than or equal to 3.

[0072] Thus N measurements (mx, my, mz)(t) are stored with t varying from 0 to N−1.

[0073] A step E11 of constructing a matrix M is then executed so that each row corresponds to a temporal sample of the measurements and each column corresponds to a component obtained with respect to the sensitive axes X, Y, Z of the sensor.

[0074] $\mathrm{Thus}\text{:}$ $\begin{array}{ccc}M=[\mathrm{mx}\left(0\right)& \mathrm{my}\left(0\right)& \mathrm{mz}\left(0\right)\\ \mathrm{mx}\left(1\right)& \mathrm{my}\left(1\right)& \mathrm{mz}\left(1\right)\\ \dots & & \\ \mathrm{mx}\left(N\right)& \mathrm{my}\left(N\right)& \mathrm{mz}\left(N\right)]\end{array}$

[0075] At the same time, a step E12 of calculating the mean value obtained for each component with respect to the sensitive axes X, Y, Z is executed.

[0076] Thus in practice a mean value of the measurements obtained is calculated for each column of the matrix M.

[0077] The following mean values are determined:

[0078] Mean(mx), Mean(my), Mean(mz)

[0079] A homogenization step E13 is then executed to subtract the mean value calculated in the step E12 from each element of the matrix M to obtain a matrix Mc:

[0080] $\begin{array}{ccc}Mc=[\mathrm{mx}\left(0\right)\text{-}\mathrm{Mean}\left(\mathrm{mx}\right)& \mathrm{my}\left(0\right)\text{-}\mathrm{Mean}\left(\mathrm{my}\right)& \mathrm{mz}\left(0\right)\text{-}\mathrm{Mean}\left(\mathrm{mz}\right)\\ \mathrm{mx}\left(1\right)\text{-}\mathrm{Mean}\left(\mathrm{mx}\right)& \mathrm{my}\left(1\right)\text{-}\mathrm{Mean}\left(\mathrm{my}\right)& \mathrm{mz}\left(1\right)\text{-}\mathrm{Mean}\left(\mathrm{mz}\right)\\ \dots & & \\ \mathrm{mx}\left(N\right)\text{-}\mathrm{Mean}\left(\mathrm{mx}\right)& \mathrm{my}\left(N\right)\text{-}\mathrm{Mean}\left(\mathrm{my}\right)& \mathrm{mz}\left(N\right)\text{-}\mathrm{Mean}\left(\mathrm{mz}\right)]\end{array}$

[0081] It is possible to effect an analysis into principal components of the matrix Mc obtained in an analysis step E14.

[0082] Determining the axis showing the least variation by a principal component analysis is known in the art and need not be described in detail here. See in particular the description of this technique in Jolliffe I. T., Principal component analysis, 2nd edition, Springer 2002 or Jackson J. E., A User's Guide to Principal Components, John Wiley and Sons, 1991, p. 4-25.

[0083] In theory, this principal components analysis amounts to seeking a 3×3 unitary matrix V in which each column vector is of unitary norm and orthogonal to the other column vectors, a 3×3 diagonal ordered matrix SIGMA with positive elements of the diagonal in decreasing order, and an N×N unitary matrix U such that:
Mc=U SIGMA V^T

[0084] These matrices can be obtained by decomposing into singular values the matrix Mc or decomposing into eigenvalues and eigenvectors the matrix S=(N−1) Mc^TMc, so that S V=V LAMBDA, where LAMBDA is a diagonal matrix that contains the eigenvalues of S. Moreover, LAMBDA (i, i)=(N−1)⁻¹ SIGMA (i, i)².

[0085] The columns of the unitary matrix V are called factorial axes and the rotation axis looked for is expressed in the frame of reference associated with the captors, corresponding to the projection axis on which the measurements have the lowest variance, corresponds to the third column of the unitary matrix V.

[0086] The value LAMBDA (3, 3) in the third column of the diagonal matrix LAMBDA corresponds to the variance associated with the axis determined in this way.

[0087] Thus it is possible to estimate the substantially invariant rotation axis in an estimation step E15 from the factorial axis of the third column of the unitary matrix V, the rotation axis therefore being expressed in the frame of reference associated with the sensors. Hereinafter, this rotation axis is denoted by its coordinates in the frame of reference of the sensor, as follows:
(nx,ny,nz)

[0088] To validate the rotation axis estimated in this way, the detection method further includes a step E16 of calculating a rotation axis variation indicator.

[0089] Thus the value d is calculated such that:

[0090] $d=\frac{\mathrm{LAMBDA}\left(3,3\right)}{\mathrm{LAMBDA}\left(1,1\right)+\mathrm{LAMBDA}\left(2,2\right)+\mathrm{LAMBDA}\left(3,3\right)}$

[0091] This value d of the indicator is an indication as to the validation of the estimated rotation axis as the invariant rotation axis of the motion.

[0092] The smaller the value d, the more valid the hypothesis whereby the motion of the mobile occurs around a constant or substantially invariant rotation axis.

[0093] In practice, the value of LAMBDA (3, 3) of the indicator is compared with the value of the noise affecting the physical measurements obtained by the sensor. If the value of LAMBDA (3, 3) of the indicator is substantially equal to or of the same order of magnitude as the noise value, the estimated rotation axis (nx, ny, nz) is validated as the substantially invariant rotation axis of the motion of the mobile body.

[0094] Also, the variations with respect to the other axes must be high: LAMBDA (3, 3)+LAMBDA (2, 2)+LAMBDA (1, 1) must be large, which guarantees that the motion is of sufficiently large amplitude.

[0095] In practice, if d is close to ⅓, the three axes have equivalent variations, and nothing can be concluded. The “constant rotation axis” hypothesis is retained only if d<threshold, which threshold is taken as equal to 0.1, for example.

[0096] Another indicator can be d*=LAMBDA (3, 3)/LAMBDA (2, 2), which should be small compared to 1. The “constant rotation axis” hypothesis is retained only if d*<threshold, which threshold is taken as equal to 0.1, for example.

[0097] Other methods can be used for estimating the lowest variation axis in the measurements obtained by an inertial or magnetic sensor.

[0098] The least squares method can also be used, for example.

[0099] In theory, the least squares method amounts to estimating the axis n on which the projection of the measurements M is constant. That axis n can be calculated from the following formula:
n=pinv(M)O

[0100] where pinv(M) is the generalized inverse of the matrix M as constructed in the step E11 described with reference to FIG. 1,

[0101] n is a column vector with three components n=(nx, ny, nz)^T, and

[0102] O is a vector of N rows and 1 column, all elements of which are equal to 1.

[0103] The vector n is normalized to obtain the director vector of the normal to the plane.

[0104] Note that the rotation axis detection method described above can generally be executed for different types of inertial or magnetic sensors, for example an accelerometer with three sensitive axes, a magnetometer with three sensitive axes or a rate gyro with three sensitive axes.

[0105] However, in the case of a rate gyro type inertial sensor, the estimated rotation axis correspond not to the axis showing the least variation of the measurements but on the contrary to the axis showing the greatest variation. The rotation axis variation indicator used then corresponds to the variance or the standard deviation calculated for the two axes orthogonal to the estimated rotation axis.

[0106] Of course, a single inertial or magnetic sensor can detect the constant rotation axis in the physical measurements at the output of this sensor if the value of the rotation axis variation indicator is substantially equal to the value of the noise affecting the associated physical measurements.

[0107] In particular, when the sensor is an accelerometer, if the standard deviation of the measurements projected onto the estimated rotation axis is substantially equal to the noise affecting the sensor, the estimated rotation axis is validated as the substantially invariant rotation axis of the motion and additionally the motion is considered planar.

[0108] In practice, if the mobile body is equipped with a plurality of sensors and, for example, at least a magnetometer with three sensitive axes and an accelerometer with three sensitive axes, it is possible to estimate the rotation axis by the method described above in the physical measurements space associated with each of the sensors, i.e. with the magnetometer and with the accelerometer.

[0109] There is then validated as the substantially invariant rotation axis of the motion a linear combination of the estimated rotation axes for the physical measurements of the magnetometer and the accelerometer.

[0110] The coefficients of the linear combination are not strictly between 0 and 1 and their sum is equal to 1.

[0111] These coefficients are a function of the values of the rotation axis variation indicators of the magnetometer and/or the accelerometer.

[0112] If one of the coefficients of the combination is equal to 0, this amounts to using the sensor (either the magnetometer or the accelerometer) producing the most reliable estimate of the substantially invariant rotation axis.

[0113] It is therefore possible to circumvent the situation in which the rotation is about the axis of the physical field concerned and associated with the sensor, i.e. the gravitational field in the case of accelerometers and the terrestrial magnetic field in the case of magnetometers.

[0114] In this case, the estimate of the rotation axis will be highly inadequate for the sensor concerned and the detection method will automatically validate the estimated rotation axis for the other sensor as the substantially invariant rotation axis of the motion.

[0115] For example, in the case of motion about a vertical axis (or in a horizontal plane), like that of a mouse on a desk, using an accelerometer as the inertial sensor in the detection method would not be satisfactory. On the other hand, the vertical rotation axis of the motion could be estimated using a magnetometer.

[0116] A second embodiment of the invention is described with reference to FIG. 3 in which the rotation axis detection method begins with a step of vectorial calculation applied to the physical measurements obtained for the three sensitive axes of the sensor.

[0117] In theory, each measurement is considered as a vector in a three-dimensional space.

[0118] The constant rotation axis in the measurements space is determined by obtaining the intersection of a number of planes each corresponding to all possible rotation axes for moving from one vector to another. In practice, the vectors are grouped and analyzed two by two in order to determine each plane.

[0119] The distance between the planes at the intersection is an indicator of the variations of the rotation axis.

[0120] As in the previous embodiment, by calculating this indicator and comparing the value of this indicator with the noise value of the physical measurements from the sensor, it is possible to validate the estimated rotation axis, and in particular, the axis corresponding to the intersection of the planes, as the substantially invariant rotation axis of the motion if the value of the indicator is substantially equal to the noise value of the physical measurements from the sensor.

[0121] In practice, a storage step E31 acquires physical measurements for the three sensitive axes of the sensor.

[0122] Thus N measurements are obtained for the three sensitive axes:
mes(t)=(mx,my,mz)(t), t=0, . . . , N−1

[0123] A step E32 of calculating all the solution planes P(t) is then executed in the following manner.

[0124] The following three vectors are defined for each time t:

[0125] $\overrightarrow{v}\left(t\right)=\frac{\left(\mathrm{mes}\left(0\right)+\mathrm{mes}\left(t\right)\right)}{\mathrm{mes}\left(0\right)+\mathrm{mes}\left(t\right)};$ $\overrightarrow{w}\left(t\right)=\frac{\left(\mathrm{mes}\left(0\right)\bigwedge \mathrm{mes}\left(t\right)\right)}{\mathrm{mes}\left(0\right)\bigwedge \mathrm{mes}\left(t\right)};$ $\begin{array}{c}\overrightarrow{p}\left(t\right)=\overrightarrow{v}\left(t\right)\bigwedge \overrightarrow{w}\left(t\right)\\ =|\begin{array}{c}\mathrm{px}\left(t\right)\\ \mathrm{py}\left(t\right)\\ \mathrm{pz}\left(t\right)\end{array}\end{array}$

[0126] In practice, at each time t, the axis of the rotation that the mobile body has effected between the time t=0 and the time t necessarily belongs to the plane defined by the vectors {right arrow over (v)}(t) and {right arrow over (w)}(t), which again corresponds to the equation of a plane {right arrow over (p)}(t). {right arrow over (n)}=0.

[0127] Thus if there is a constant rotation axis in the motion of the observed mobile body, it corresponds to the intersection of all the planes defined by the vectors {right arrow over (v)}(t) and {right arrow over (w)}(t), with t varying from 0 to N−1.

[0128] On the other hand, if the intersection exists, the hypothesis that the motion is performed about a constant rotation axis can be considered to be validated and that rotation axis can be considered to correspond to the intersection of the planes.

[0129] In practice, a step E33 of solving the system is executed starting with the construction of a matrix A the rows whereof correspond to the N samples and the columns whereof correspond to the three components of the vector {right arrow over (p)}:

[0130] $A=\left(\begin{array}{ccc}\mathrm{px}\left(0\right)& \mathrm{py}\left(0\right)& \mathrm{pz}\left(0\right)\\ \dots & & \dots \\ \mathrm{px}\left(N-1\right)& \mathrm{py}\left(N-1\right)& \mathrm{pz}\left(N-1\right)\end{array}\right)$

[0131] Seeking the rotation axis {right arrow over (n)}(nx, ny, nz) amounts to solving the equation A{right arrow over (n)}=0, that is to say:

[0132] $A\overrightarrow{n}=\left(\begin{array}{ccc}\mathrm{px}\left(0\right)& \mathrm{py}\left(0\right)& \mathrm{pz}\left(0\right)\\ \mathrm{px}\left(1\right)& \mathrm{py}\left(1\right)& \mathrm{pz}\left(1\right)\\ \dots & \dots & \dots \\ \mathrm{pz}\left(N-1\right)& \mathrm{py}\left(N-1\right)& \mathrm{pz}\left(N-1\right)\end{array}\right)\left(\begin{array}{c}\mathrm{nx}\\ \mathrm{ny}\\ \mathrm{nz}\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ \dots \\ 0\end{array}\right)$

[0133] The coordinates of the vector {right arrow over (n)} can be obtained by decomposing the matrix A into singular values as described hereinabove so that A=USV^T.

[0134] In the same way as for the principal component analysis, the diagonal terms are classified in decreasing order.

[0135] The rotation axis {right arrow over (n)} then corresponds to the third column of the unitary matrix V.

[0136] The values of the diagonal of the matrix S are used as previously to calculate the value d of an indicator of the variations of the rotation axis, which corresponds to the mean distance with the hypothesis A{right arrow over (n)}=0.

[0137] It will be noted that solving the system can be facilitated if the number N of samples is high, by first effecting orthogonal triangularization of the matrix A and thereafter applying decomposition into singular values directly to the triangularized matrix.

[0138] Moreover, the estimate of the vector n can be improved, in particular by eliminating from the matrix the measurements for which mes(t) and mes(0) are close.

[0139] If the distance between the vectors is small, the orientation of the vector w is strongly disturbed by the noise affecting the measurements, and so the general estimate of the vector {right arrow over (n)} by the matrix A is also disturbed.

[0140] To limit this effect, it is possible to eliminate from the matrix A the rows corresponding to observations such that the distance between the vectors mes(t) and mes(0) is less than a threshold value s that is predetermined as a function of the noise affecting the measurements obtained at the output of the associated sensor.

[0141] Note that the determination method using vectorial analysis can be applied to measurements obtained from magnetometers in the absence of magnetic disturbances or measurements obtained at the output of accelerometers in the absence of acceleration.

[0142] On the other hand, this type of vectorial calculation on the physical measurements is not applicable to a rate gyro type inertial vector.

[0143] In the embodiments described hereinabove, it is necessary to have a number N of samples at least equal to 3 in order to detect the rotation axis.

[0144] Moreover, the distance between the measurements must be greater than the noise affecting the measurements.

[0145] Thus for a magnetometer for which the noise is substantially equal to 0.1, it is necessary to have a distance between the measurements of at least 0.3.

[0146] In practice, this means that the detection method for determining the rotation axis is particularly suitable if the N samples cover sufficiently large rotation displacements about that rotation axis.

[0147] The detection method described hereinafter can be complemented by a step of transformation of the coordinates of the estimated rotation axis n=(nx, ny, nz) in the measurement frame of reference of the sensor into coordinates in a geophysical frame of reference.

[0148] The geophysical frame of reference can typically be the North-East-Down (NED) frame of reference.

[0149] For example, the coordinates of the rotation axis in the NED geophysical frame of reference can be calculated from measurements projected onto the rotation axis expressed in the frame of reference of the sensor in the following manner.

[0150] In the geophysical frame of reference NED, the gravitational field is written g=(0, 0, 1) and the magnetic field b=(bx, 0, bz).

[0151] The magnetic field varies as a function of the measurement location, the values bx and bz being the reference of the terrestrial magnetic field at the location where the measurement is effected.

[0152] For example, in France, bx≅1/2 and

[0153] $\mathrm{bz}\simeq \frac{\sqrt{3}}{2}.$

[0154] It is therefore the coordinates nz_NED, nx_NED, ny_NED of the vector in the geophysical frame of reference NED that are required.

[0155] If the accelerometer measurements are expressed in the new frame of reference by effecting the matrix product newmesacc=mes_aV, the third component should be seen to be quasi-constant. It reflects the inclination of the plane in the inertial frame of reference and corresponds to the coordinate nz_NED.
nz_NED=mean(newmesaccz)

[0156] Similarly, if the accelerometer measurements are expressed in the new frame of reference by effecting the matrix product newmesmag=Mes_bV, the third component should be seen to be quasi-constant. It is the conjoint reflection of the inclination of the plane and its orientation relative to local magnetic north in the inertial frame of reference, and corresponds to the coordinate z.

[0157] We have tmp=mean (newmesmagz)
tmp=0.5nx_NED+0.86nz_NED
or nx_NED=2*tmp−sqrt(3)*nz_NED

[0158] Finally, the coordinate ny_NED can be determined by seeking the normed vector with the coordinates determined previously:
ny_NED=±√{square root over (1−nx_NED²−nz_NED²)}

[0159] Starting from detecting a substantially invariant rotation axis by the method described above, it is possible to estimate the complete motion of a mobile body equipped with at least one inertial or magnetic sensor with three sensitive axes.

[0160] The method for estimating the motion of a mobile is illustrated by FIG. 4 in particular.

[0161] After a step E40 of storing N measurements, corresponding to the storage steps E10 or E31 in FIGS. 1 and 3, a detection step E41 is executed as described hereinabove in order to detect a rotation axis represented by its coordinates n=(nx, ny, nz).

[0162] It is then possible to calculate successively the rotation angle in a calculation step E42 and then the acceleration of the mobile body in a calculation step E43, for each measurement sample at times t between 0 and N−1 and for the subsequent samples if the rotation axis no longer varies.

[0163] In practice, the step E42 of calculating the rotation angle can be implemented in the following manner.

[0164] Once the rotation axis has been estimated, the angle θ of the rotation effected can be calculated at any time by an analytical calculation:

[0165] $\theta =+2\arcsin \frac{\mathrm{mes}\left(0\right)-\mathrm{mes}\left(t\right)}{\sqrt{{\mathrm{mes}\left(0\right)-\mathrm{mes}\left(t\right)}^2+{\mathrm{mes}\left(0\right)+\mathrm{mes}\left(t\right)}^2{\left(\overrightarrow{w},\overrightarrow{n}\right)}^2}}$

[0166] The symbol θ corresponds to that for reconstructing the measurement mes(t) as well as possible.

[0167] A second method estimating θ by optimization could be used:
min_θ(∥mes(t)−M({right arrow over (n)},θ)mes(0)∥²)

[0168] where M({right arrow over (n)},θ) is the rotation matrix of axis n and angle θ.

[0169] The matrix V obtained by principal component analysis as described with reference to step E14 in FIG. 1 can also be used to change the data frame of reference.

[0170] Using the same notation as before, measurements in the new frame of reference are given by M′=VM. The third component will be constant in this frame of reference.

[0171] For each time (row of matrix M′ of N rows by 3 columns), the angle of the rotation is the angle formed by the two 2D vectors (their coordinates correspond to the first two columns).

[0172] In practice, what is required is the angle between two vectors corresponding to two rows of the matrix that correspond to the rotation effected between those two times.

[0173] To this end, these vectors can be expressed in polar coordinates (length, orientation), for example, and then the difference of their orientations obtained. Their scalar product can also be used.

[0174] In the case of an accelerated movement, once the orientation of the mobile object has been estimated at each time t using a magnetometer or a goniometer, the acceleration of the mobile body can be estimated using an accelerometer.

[0175] Knowing that the accelerometer measurements are written as follows:
mes(t)=M(n,θ)(g−1)=M(n,θ)(mes(0)−a)

[0176] where M(n,θ) is the rotation matrix for going from the fixed frame of reference to the mobile frame of reference,

[0177] θ is the rotation angle,

[0178] a is the required acceleration, and

[0179] g is the gravitational field equal to mes(0) in the reference frame of reference at time 0.

[0180] The acceleration a of the mobile object can be estimated in the fixed frame of reference by:
{tilde over ({right arrow over (a)}=−M({right arrow over (n)},θ)^Tmes(t)+mes(0)
where a and mes(0) are expressed in the same fixed frame of reference.

[0181] In the particular case of planar motion, the acceleration a must be contained in the plane orthogonal to n. The problem to be solved then becomes an optimization with constraints:
min_ã(∥mes(t)−M({right arrow over (n)},θ)(mes(0)−{tilde over ({right arrow over (a)}∥²)
where {right arrow over (n)}·{tilde To find H, the triangular orthogonal decomposition of n can be effected. Let {tilde over (H)}=H(:2:3) be the matrix consisting of the last two columns of H.

[0182] The acceleration of the mobile object is then obtained from:
{tilde over ({right arrow over (a)}_t={tilde over (H)}{tilde over (H)}^T(−M({right arrow over (n)},θ)mes(t)+mes(0))

[0183] By way of nonlimiting example, the observed motion can be walking. The sensor is placed on the tibial platform, for example. Precise alignment of the sensor is no longer necessary if only the angle described in the sagittal plane of walking is of interest.

[0184] The motion estimation method described hereinabove can be applied to each time segment corresponding to a stride or a step.

[0185] Moreover, thanks to the method for detecting an invariant axis described hereinabove, it is possible, by using this method, to calibrate a system of sensors, i.e. to determine a rotation matrix adapted to transform a measurement frame of reference of the sensor into a predetermined frame of reference stationary with respect to a mobile object carrying the sensor.

[0186] As shown in FIG. 5, if the sensor 50 is integrated into a casing 51, the axes xb, yb, zb of the predetermined frame of reference Rb are defined by edges of the casing 51.

[0187] The calibration method to be described hereinafter determines the rotation matrix for moving from the frame of reference Rc associated with the sensor 50 to the predetermined frame of reference Rb stationary with respect to the casing.

[0188] Similarly, as shown in FIG. 6, the sensor 50 can be placed on an anatomical feature, for example the forearm of a person.

[0189] The predetermined frame of reference then corresponds to an anatomical frame of reference defined by three anatomical rotation axes of that anatomical feature.

[0190] As shown clearly in FIG. 6, the axes of the anatomical frame of reference can be defined by the respective rotation axes of a bending movement of the forearm, a pronosupination movement of the forearm, and an abduction/adduction movement of the forearm.

[0191] In theory, the calibration method effects successive rotations about the three axes of the predetermined frame of reference Rb.

[0192] By then seeking for each of the three motions the coordinates of the corresponding rotation axis in the frame of reference of the sensor Rc it is possible to generate a rotation matrix in which the coordinates of each of the rotation vectors found correspond to a column of the rotation matrix for moving from the frame of reference of the sensor Rc to the predetermined frame of reference Rb stationary with respect to the mobile body.

[0193] As shown in FIG. 7, the calibration method includes in particular steps E71, E72, E73 for acquisition of physical measurements by the sensor, respectively corresponding to each rotation about the axes xb, yb, zb of the frame of reference Rb.

[0194] For each rotation motion, N measurement samples are stored as described hereinabove with reference to the storage step E10 shown in FIG. 1, for example.

[0195] From each stored set of N measurements a determination step E74, E75, E76 is executed in order to determine the coordinates of the three rotation axes n1, n2, n3 corresponding to the axes of the predetermined frame of reference Rb in the frame of reference Rc of the physical measurements from the sensor. These steps of determining the rotation axis E74, E75, E76 are executed as described hereinabove by the method for detecting a substantially invariant rotation axis.

[0196] A step E77 of generating a matrix Rbc is then executed using the coordinates of the rotation axes determined in the physical measurement frame of reference Rc:
Rbc=[nx1ny1nz1nx2ny2nz2nx3ny3nz3]

[0197] During the step E77 of generating the matrix Rbc, it is also possible to orthogonalize the matrix formed in this way by the three vectors n1, n2, n3 to render it unitary. The matrix Rbc then corresponds to a rotation matrix enabling measurements in the sensor frame of reference Rc to be expressed in the predetermined frame of reference Rb stationary with respect to the mobile body.

[0198] Accordingly, the measurements in the new frame of reference Rb can be expressed in the following manner:
mes_b=Rbcmes_c

[0199] in which mes_b corresponds to measurements in the frame of reference associated with the casing and mes_c corresponds to measurements in the frame of reference associated with the sensor.

[0200] It will be noted that in the method described hereinabove only two of the three rotations about axes of the predetermined frame of reference Rb can be effected, the coordinates of the third rotation axis then being obtained by orthogonalizing the matrix formed by the coordinates of the first two rotation axes determined by the detection method for the invention.

[0201] Also, the initial frame of reference may not be orthogonal (the axes are not mutually orthogonal). If the final frame of reference is orthogonal, then the matrix Rbc that effects the change of frame of reference then changes from data expressed in a non-orthogonal frame of reference to data expressed in an orthogonal frame of reference. For this it is essential to effect the three rotations and not to orthogonalize the matrix Rbc obtained.

[0202] The choice can also be made to pass through an intermediate orthogonal frame of reference and then to effect a change between two orthogonal frames of reference (in which case only two of the three rotations are necessary).

[0203] Accordingly, thanks to the method for the invention of detecting a substantially constant rotation axis, it is possible to reduce the number of degrees of freedom in the motion of an observed mobile object and thus to facilitate solving the systems for determining the parameters of that motion.

[0204] This method uses an inertia center comprising fewer inertial or magnetic sensors than prior art methods and in particular no rate gyro.

[0205] The method for estimating a motion using the detection of a constant rotation axis is particularly suitable for the situation in which the rotation axis varies slowly provided that the N measurement samples are stored over a small time window. It is thus possible, over each time window, to determine a substantially invariant rotation axis that is substantially invariant over that time window provided that there exist measurements that are sufficiently far apart over that time window, i.e. having a sufficient distance with respect to the noise affecting the measurements from the sensor.

[0206] Thus for each small time window the greatest distance between the measurement vectors taken two by two is calculated.

[0207] If that distance is less than a predetermined threshold that depends on the noise affecting the sensors, then the time window is enlarged and a sample of N measurements is stored over this new window.

[0208] If not, if the maximum distance between the measurement vectors is greater than the noise affecting the sensors, the rotation axis is calculated as indicated hereinabove with reference to the method for detecting a substantially invariant rotation axis and the indicator of variations of that rotation axis in order to validate the existence of a constant rotation axis in the observed motion.

[0209] Where appropriate, the angle and the acceleration of the motion over this time window can also be estimated if the motion of a mobile body in time is to be observed.

[0210] The set of operations is then repeated over a time window shifted in time.