PARTICLE FILTERING AND NAVIGATION SYSTEM USING MEASUREMENT CORRELATION

Abstract

Disclosed is a box-regularized particle filtering process which includes an Epanechnikov kernel smoothing step. For this purpose, the process uses a special method for generating random numbers that follow an Epanechnikov probability density function. The process can be performed autonomously in a navigation system using correlation measurement, in particular on board an aircraft such as an aircraft, a flying drone or any self-propelled aerial carrier.

Claims

1. A box-regularized particle filtering method, for predicting a state of a system by a set of state intervals with weights associated with said state intervals, so as to form a probability distribution that characterizes the state of the system, said system being a land, air, sea or space vehicle that is provided with a navigation system using measurement correlation, the method comprising repetitively applying a sequence of steps to the set of state intervals with the associated weights for updating the state intervals and said associated weights, the sequence of steps comprising a so-called smoothing step that consists in modifying at least one of the state intervals by applying random modifications to a set of interval bounds, or center values and interval lengths, that determine the state interval according to state coordinates of the system, wherein the random modifications relating to each state interval to be modified, which is identified by an integer index i, are determined by executing the following steps: generating a first random value, denoted β.sub.i and comprised between 0 and 1, the values 0 and 1 being permitted, according to a statistical-distribution beta law with a first parameter equal to d and a second parameter equal to 2, where d is a number of state coordinates of the system; generating 2.Math.d second random values, denoted v.sub.k,i, each according to a normal statistical distribution law with zero mean value and standard deviation equal to unity, where k is another integer index that varies from 1 to 2.Math.d and identifies the interval bounds, or center values and interval lengths, for each state interval; calculating a first number, denoted ξ.sub.i, according to the first formula: ξ.sub.i = [Σ.sub.k=1 .sub.to .sub.2.Math..sub.d (V.sub.k,i).sup.2].sup.½., calculating a second number, denoted α.sub.i, according to the second formula: α.sub.i = β.sub.i.sup.½/ξ.sub.i; and calculating 2.Math.d third numbers, denoted ε.sub.k,.sub.i, according to the third formula: ε.sub.k,.sub.i = V.sub.k,i•α.sub.i, and wherein the random modifications that are applied to the state interval i are proportional one-to-one to the third numbers ε.sub.k,i, with a proportionality coefficient that is non-zero and common to said random modifications.

2. The method according to claim 1, wherein all the interval bounds, or center values and interval lengths, which determine the state interval i, are modified by implementing the following steps: combining the random modifications that relate to said state interval i using a square matrix of dimension 2.Math.d, so as to obtain 2.Math.d combinations of random modifications; then adding said combinations of random modifications that are thus obtained one-to-one to the interval bounds, or center values and interval lengths, of the state interval i.

3. The method according to claim 2, wherein the matrix that is used for combining the random modifications is such that the product of said matrix multiplied by a transpose of said matrix is equal to a mean product matrix, said mean product matrix being square of dimension 2.Math.d, and having as coefficients mean values that are calculated over all the state intervals, of products of the interval bounds, or center values or interval lengths, taken in pairs separately for each state interval.

4. The method according to claim 1, wherein each first random value β.sub.i is generated using an algorithm that combines: a generation of two random numbers each according to a uniform statistical distribution law; and at least one acceptance criterion that is based on the two random numbers, such that, if said at least one acceptance criterion is satisfied, a first of the two random numbers is used for calculating the first value β.sub.i, otherwise the generation of the two random numbers is recommenced.

5. The method according to claim 4, wherein each of the two random numbers is generated using a method of the shift register type with linear feedback.

6. The method according to claim 1, wherein each second random value V.sub.k,i is calculated as a sum of several initial random values, each of said initial random values being generated according to a uniform statistical distribution law.

7. The method according to claim 1, wherein respective estimations of each first number ξ.sub.i and of each second number α.sub.i are obtained using the following steps at least once, where X is a positive or zero variable number and α is an exponent value equal to 2 or ½: /a/ writing the number X in the form X = (1+m).2.sup.ex, where ex is a negative, positive or zero integer, and m is a mantissa comprised between 0 and 1, the value 0 and being permitted, so that a binary representation of the number X is: I(X) = L.Math.(m + ex + B), where L=2.sup.n with n being a number of bits of a binary writing of the mantissa m, and B being a positive or zero constant number, referred to as bias; /b/ calculating a binary representation of X.sup.α in the form: I(X.sup.α) = α.Math.I(X) + L.Math.(1 - α).Math.(B - σ), where σ is a constant number the value of which is recorded ; and /c/ obtaining the estimation of the value of X.sup.α from the binary representation I(X.sup.α), steps /a/-/c/ being applied to X = E.sub.k=1 à .sub.2.Math.d (V.sub.k,.sub.i).sup.2 with α=½, to obtain an estimation of the first number ξ.sub.i , and steps /a/-/c/ being applied to X = β.sub.i with α=½, to obtain an estimation of the second number (α.sub.i as a result of a division of the estimation of (β.sub.i.sup.½ by the estimation of the first number ξ.sub.i.

8. The method according to claim 7, wherein the obtaining of the estimation of the value of X.sup.α is completed by at least once executing the following additional step, after the step /c/: /d/ calculating a new estimation of the value of X.sup.α from a prior estimation of the value of X.sup.α, by applying an approximate equation resolution recursive algorithm to the equation Y.sup.⅟α - X =0 of unknown Y, the estimation of the value of X.sup.α that was obtained at the step /c/ being used as prior estimation for a first application of said algorithm, and the new estimation of the value of X.sup.α that is produced by a q.sup.th application of the algorithm forming the prior estimation of the value of X.sup.α for the (q+1).sup.th application of said algorithm, if such (q+1).sup.th application of the algorithm is carried out, q being an integer greater than or equal to 1.

9. The method according to claim 1, wherein the sequence of steps that is applied repetitively for updating the state intervals with the weights that are associated with said state intervals comprises the following steps /1/ to /5/: /1/ a prediction step, comprising predicting subsequent state intervals, each subsequent state interval being obtained by applying at least one propagation rule to one of a plurality of prior state intervals; /2/ a step of measuring a true state of the system; /3/ a step of contracting at least one of the subsequent state intervals, according to at least one result of measuring the true state that was obtained at step /2/; /4/ a weight updating step, comprising allocating a weight to each subsequent state interval according to a size of said subsequent state interval as resulting from step /3/, a size of said subsequent state interval as resulting from step /1/ before step /3/, and a weight of the prior state interval from which the subsequent state interval resulted during step /1/; and /5/ a step of redistributing the state intervals, comprising replacing at least one of the subsequent state intervals by a plurality of sub-intervals that result from a division of the subsequent state interval, each sub-interval forming a new state interval, said redistribution step comprising applying the smoothing step to at least each new state interval, the state intervals as resulting from an implementation of the sequence of steps /1/-/5/, comprising the new state intervals and subsequent state intervals that were kept without being replaced by a plurality of new state intervals, being intended to constitute the prior state intervals for a following implementation of the sequence of steps /1/ to /5/.

10. A computing unit, comprising at least one first input adapted for receiving results of repeated measurements of acceleration and angular velocity of a system, and a second input adapted for receiving results of repeated measurements of a true state of the system, supplementary with respect to the acceleration and angular velocity measurements, and the computing unit being arranged for implementing a box-regularized particle filtering method that is in accordance with claim 1, so as to output a series of state intervals with respective weights, the weight that is associated with each of the state intervals corresponding to a probability of the true state of the system being in said state interval.

11. The computing unit according to claim 10, of a type consisting of a circuit with an array of programmable ports, a circuit with an array of fixed ports, or a central processor unit.

12. A navigation system using measurement correlation, adapted for being installed on board a vehicle, comprising: an inertial unit, adapted for iteratively measuring accelerations and angular velocities of the vehicle, and for deducing, by using results of measurements of accelerations and angular velocities, subsequent state intervals respectively from a plurality of prior state intervals, each state of the vehicle comprising position, speed and attitude coordinates of said vehicle; a measurement system, adapted for iteratively measuring at least one feature of a true state of the vehicle vehicle and a computing unit in accordance with claim 10, and adapted for reducing at least one drift of the inertial unit, from a position drift, a speed drift and an attitude drift, using the results of measurement of the feature of the true state of the vehicle that are delivered by the measurement system.

13. A vehicle, comprising a navigation system using measurement correlation that is in accordance with claim 12.

14. The method of claim 7, wherein steps /a/-/c/ are applied previously to an absolute value of each second random value, according to X = |V.sub.k,i|, with α=2.

15. The method according to claim 2, wherein each first random value βi is generated using an algorithm that combines: a generation of two random numbers each according to a uniform statistical distribution law; and at least one acceptance criterion that is based on the two random numbers, such that, if said at least one acceptance criterion is satisfied, a first of the two random numbers is used for calculating the first value βi, otherwise the generation of the two random numbers is recommenced.

16. The method according to claim 3, wherein each first random value βi is generated using an algorithm that combines: a generation of two random numbers each according to a uniform statistical distribution law; and at least one acceptance criterion that is based on the two random numbers, such that, if said at least one acceptance criterion is satisfied, a first of the two random numbers is used for calculating the first value βi, otherwise the generation of the two random numbers is recommenced.

17. The method according to claim 2, wherein each second random value ν.sub.k,.sub.i is calculated as a sum of several initial random values, each of said initial random values being generated according to a uniform statistical distribution law.

18. The method according to claim 3, wherein each second random value ν.sub.k,.sub.i is calculated as a sum of several initial random values, each of said initial random values being generated according to a uniform statistical distribution law.

19. The method according to claim 4, wherein each second random value ν.sub.k,.sub.i is calculated as a sum of several initial random values, each of said initial random values being generated according to a uniform statistical distribution law.

20. The method according to claim 5, wherein each second random value ν.sub.k,.sub.i is calculated as a sum of several initial random values, each of said initial random values being generated according to a uniform statistical distribution law.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0056] The features and advantages of the present invention will emerge more clearly in the following detailed description of non-limitative example embodiments, with reference to the accompanying figures, among which:

[0057] FIG. 1 shows an aircraft that is equipped with an inertial unit using terrain correlation according to the invention;

[0058] FIG. 2 is a diagram that shows the concatenation of steps of a method according to the invention;

[0059] FIG. 3 is a diagram that details the implementation of a smoothing step according to the invention; and

[0060] FIG. 4is a flow diagram of an algorithm that can be used for implementing the invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0061] For clarity sake, the dimensions of the elements that are shown symbolically in FIG. 1 correspond neither to real dimensions nor to ratios of real dimensions. Furthermore, the invention is described by way of non-limitative example for a case of application to an aircraft, but it is understood that it can be applied to any vehicle that is provided with a navigation system using measurement correlation, whether this vehicle be terrestrial, airborne, maritime, space, etc.

[0062] A method for calculating an estimation of the value of X.sup.α that can be used in the invention is described first. X is a number of variable value, positive or zero, and α designates an exponent that can be equal to 2 or ½.

[0063] In a known manner, the number X can be written uniquely in the following form, in accordance with IEEE 754:

$\text{X = (1+m)}\cdot {\text{2}}^{\text{ex}},$

where ex is a positive or zero integer, and m is a mantissa comprised between 0 and 1, the value zero also being permitted. The number ex and the mantissa m thus depend on the value of the number X.

[0064] This thus gives

${\text{log}}_2\text{X = ex + m +}\sigma \text{,}$

where σ is a fixed real number making it possible to minimize an error of the value of log.sub.2X, in particular when a numeric interval thereof is known a priori for the number X. For example, the value of the number σ may be taken to be equal to 0.043036. And therefore:

${\text{ex + m = log}}_2\text{X -}\sigma .$

[0065] Moreover, the number X can be represented in binary form by I(X) defined by:

$\text{I(X) = L}\cdot \text{(m + ex + B),}$

where L=2.sup.n, with n being a number of bits set for writing the mantissa m in binary form, and B is a positive or zero constant number, which is referred to as bias. In the representation I(X) of the number X, L, m, ex and B are expressed in binary form. For example, n may be equal to 23, and B may be equal to 127. By transferring into the binary representation I(X) the expression of ex + m as coming from log.sub.2X, this gives:

$\text{I(X) = L}\cdot \left[{\text{log}}_2\text{X -}\sigma +\text{B}\right],\text{i}\text{.e}{\text{.: log}}_2\text{X = I(X)/L +}\sigma \text{- B}\text{.}$

[0066] However, in the same way as I(X) in the previous line, the binary representation of X.sup.α is:

${\text{I(X}}^{\alpha })=\text{L}\cdot \left[{\text{log}}_2\left({\text{X}}^{\alpha }\right)-\sigma +\text{B}\right]$

However, log.sub.2(X.sup.α) = α.Math.log.sub.2(X), therefore: I(X.sup.α) = L.Math.[α.Math.log.sub.2(X) - σ + B], and, by replacing log.sub.2(X) with its expression as a function of the binary representation I(X), this gives:

$\begin{array}{l}\text{I}\left({\text{X}}^{\text{a}}\right)=\text{L}\cdot \left[\text{a}\cdot \left(\text{I}\left(\text{X}\right)/\text{L +}\sigma \text{- B}\right)\text{-}\sigma \text{+ B}\right],\text{i}\text{.e}\text{.: I}\left({\text{X}}^{\text{a}}\right)=\\ \text{a}\cdot \text{I}\left(\text{X}\right)\text{+ L}\cdot \left(1\text{-}\text{a}\right)\cdot \left(\text{B -}\sigma \right)\end{array}$

An approximate value of X.sup.α, denoted Y.sub.0, can then be reconstructed from the binary representation of X.sup.α that was thus obtained, using a method that is the inverse of the one that supplies the binary representation of a number from this number. The difference between this approximate value Y.sub.0 and the true value of X.sup.α depends on the value that was adopted for the number σ. For many applications, the approximate value Y.sub.0 is suitable in a satisfactory manner for replacing X.sup.α, considering the simplicity of the method for obtaining this approximate value Y.sub.0, as has just been described.

[0067] For the particular case where the exponent α is equal to 2:

$\text{I}\left({\text{X}}^2\right)=2\cdot \text{I}\left(\text{X}\right)\text{- L}\cdot \left(\text{B -}\sigma \right)$

[0068] For the particular case where the exponent α is equal to ½:

$\text{I}\left({\text{X}}^{1/2}\right)=\text{I}\left(\text{X}\right)/2\text{+ 0,5}\cdot \text{L}\cdot \left(\text{B -}\sigma \right)$

[0069] For applications where the approximate value Y.sub.0 does not constitute a sufficiently precise evaluation of X.sup.α, it is possible to improve this evaluation by using one of the methods for refining approximate values that are known to a person skilled in the art. Newton‘s algorithm, also referred to as Newton’s method, may be used in particular, by applying it to the function f(Y) = Y.sup.⅟α - X and to the equation f(Y) = 0. Successive approximate values Y.sub.q, q being an integer index for numbering these values, can thus be obtained according to the formula: Y.sub.q+1 = Y.sub.q - f(Y.sub.q)/f’(Y.sub.q), where f‘(Y.sub.q) is the value of the function derived from f, estimated for the value Y.sub.q. That is, when calculating the expression of f’(Y) from the expression of f(Y):

${\text{Y}}_{\text{q+1}}=\left(\text{1 - a}\right)\cdot {\text{Y}}_{\text{q}}+\text{a}\cdot \text{X}\cdot {\text{Y}}_{\text{q}}^{\left(\text{a-1}\right)/\text{a}},\text{for q = 0, 1, 2,}$

[0070] For the particular case where the exponent α is equal to 2, this gives:

${\text{Y}}_{\text{q+1}}{\text{= -Y}}_{\text{q}}\text{+ 2}\cdot \text{X}\cdot {\text{Y}}_{\text{q}}^{1/2}$

In particular, the first-order approximate value of X.sup.2 is:

${\text{Y}}_1=-{\text{Y}}_0+2\cdot \text{X}\cdot {\text{Y}}_0{}^{1/2}.$

The value of Y.sub.q.sup.½ can be estimated each time using the formula I(X.sup.α) = α.Math.I(X) + L.Math.(1 - α).Math.(B - σ), and by replacing in this formula X by Y.sub.q and α by ½.

[0071] For the particular case where the exponent α is equal to ½, this gives:

${\text{Y}}_{\text{q+1}}=0,5\cdot {\text{Y}}_{\text{q}}+0,5\cdot {\text{X/Y}}_{\text{q}}.$

In particular, the first-order approximate value of X.sup.½ is:

${\text{Y}}_1=0,5\cdot {\text{Y}}_0+0,5\cdot {\text{X/Y}}_0.$

[0072] The approximate value Y.sub.0, as obtained using the binary representations of numbers, and the values Y.sub.q, .sub.q≥1, as obtained using one of the methods for refining approximate values such as Newton’s method, does not require extensive computing resources. They can therefore be calculated easily by a computing unit of the FPGA, DSP, CPU or RISC type.

[0073] In accordance with FIG. 1, an aircraft 20 is equipped with a navigation unit using terrain correlation, designated by the reference 10. The navigation unit 10 comprises an inertial unit 1, a telemetry sensor 2 and a computing unit 3. In a known manner, the inertial unit 1 repetitively makes measurements of three acceleration coordinates and of three angular velocity coordinates of the aircraft 20, by means of accelerometers and gyrometers, not shown. Moreover, the telemetry sensor 2 repetitively makes measurements of the distance H that exists between the ground 100 and the aircraft 20. This distance H is measured in a direction that may be fixed with respect to the aircraft 20, for example perpendicular to a reference plane of the aircraft. Optionally, the telemetry measurement direction may be variable with respect to the aircraft 20, but in such case this measurement direction is controlled and taken into account in a suitable manner, and known moreover. The distance H therefore varies depending on the terrestrial relief over which the aircraft 20 is flying, and also depending on its altitude and attitude. Hereinafter, the present description is limited solely to the case where the measurements that are supplementary with respect to those of the inertial unit 1 consist of the measurements of the distance H. This case corresponds to a navigation with terrain correlation. However, it is understood that other supplementary measurements may be used alternatively to those of the distance H, or in addition to them. Finally, the invention is compatible with inertial unit and telemetry sensor models as commercially available. In particular, the inertial unit may be of an MEMS type, standing for “Micro-Electro-Mechanical System”, of the quadrason type, of the gyrolaser type, etc., and the telemetry sensor may be a radio altimeter, a laser telemeter, etc.

[0074] Each possible state for the aircraft 20 may be composed of three spatial coordinate values that identify a position for the aircraft, for example in the terrestrial reference frame, three velocity values, each according to one of the spatial coordinates, and three angular values for identifying an orientation of the aircraft, for example three Euler angle values, for a total of nine state coordinates. Under these conditions, a state interval for the aircraft 20 is formed by a combination of nine unidimensional intervals that each relate separately to one of the state coordinates. Such a state interval is referred to as a box in the jargon of a person skilled in the art.

[0075] The box-regularized particle filtering method is initialized by the supply of several initial state intervals, each associated with a weight that indicates a probability of the true state of the aircraft 20 being located initially in this initial state interval. Thus, N initial state intervals are supplied, N being an integer, for example equal to 16 or 32, preferably less than or equal to 128. Each initial state interval is individually associated with a weight value, which may be equal to 1/N.

[0076] The method next consists in successive iterations of a sequence of steps, each new implementation of the sequence of steps producing an updating of state intervals, with updated values of weights that are associated one-to-one with the updated state intervals. Furthermore, each new implementation of the sequence of steps is executed from the state intervals and their associated weight values as supplied by the just prior implementation of the sequence of steps.

[0077] Each sequence of steps comprises a prediction step denoted /1/ in FIG. 2, a measurement step denoted /2/, a contraction step denoted /3/, a step /4/ of updating the weight of each state interval, and a step /5/ of redistribution of the state intervals. Steps /1/ and /3/ to /5/ are implemented for each state interval. Since step /5/ produces a redistribution of the state intervals as resulting from step /3/, it is preferably implemented so as to keep a constant number of state intervals. Then the computing unit 3 can be designed and sized to process N state intervals at each implementation of the sequence of steps /1/ to /5/. i is an integer index from 1 to N, which counts the state intervals that are processed at each iteration of this sequence of steps /1/ to /5/.

[0078] For each state interval i, denoted box_i, the prediction step /1/ consists in collecting the results of last measurements of acceleration and angular velocity, as delivered by the inertial unit 1. Optionally, the results of several measurements that were made by the inertial unit 1 since the previous implementation of the sequence of steps /1/ to /5/ can be collected. For each of the box_i state intervals, this step /1/ also comprises calculating a change in this state interval during the period of time that has elapsed between the two implementations of the step /1/, relating to the box_i state interval, from the results of the measurements of acceleration and angular velocity, is assumed to be known to a person skilled in the art. In this regard reference can be made to the article by Merlinge, N., Dahia, K., Piet-Lahanier, H., Brusey, J., & Horri, N., which is entitled “A Box Regularized Particle Filter for state estimation with severely ambiguous and non-linear measurements”, Automatica (2019), Vol. 104, pp. 102-110. This step /1/ results in a translation, usually accompanied by a change in length, of each unidimensional interval that relates to one of the state coordinates of the aircraft 20. In the general part of the present description, each box_i state interval as existing at the moment of starting the execution of step /1/ has been called prior state interval, and this box_i state interval as modified by the step /1/ has also been called subsequent state interval.

[0079] The measurement step /2/ consists in collecting the result of the last distance measurement H as delivered by the telemetry sensor 2. Optionally, the results of a plurality of last measurements that were made by the telemetry sensor 2 since the previous execution of step /2/, at the time of the previous iteration of the sequence of steps /1/ to /5/, may be collected.

[0080] The contraction step /3/ consists in reducing the size of each box_i state interval according to any incompatibilities that might exist between parts of this state interval and the last result or results of distance measurement H collected at step /2/. Typically, the six position and attitude coordinates of the aircraft 20, such as possibly varying within each box_i state interval resulting from step /1/, are combined with relief height values read in a terrain map that is stored, in order to obtain an estimation of the distance H associated with each state. Such calculation may be carried out in the known manner that has been reminded at the start of this description. The distance H thus estimated is compared with the result of the measurement of step /2/. Generally, an implicit equation can be used for converting each state of the aircraft 20 into a distance value H, using the terrain map stored. However, such equation, referred to as an observation equation by a person skilled in the art, may be difficult to invert locally, so that modeling its inverse function by an analytical function part and a tabulated function part may be advantageous. A person skilled in the art will also be able to refer to the article by Merlinge et al., Automatica 104(2019), cited above, with regard to step /3/, to which the improvement to the method procured by the invention does not directly relate.

[0081] The updating step /4/ consists in updating the value of the weight of each state interval as resulting from step /3/. For example, a new value of the weight of the box_i state interval may be equal to the value of the weight of this box_i state interval as existing before implementing this update, multiplied by the size of the box_i state interval as resulting from the contraction step /3/, and divided by the size of the box_i state interval as resulting from the prediction step /1/ before applying the contraction step /3/. Other formulae for updating the weight values may be adopted alternatively. Optionally, the weight values as resulting from one of these formulae may also be corrected, for example by multiplying them by a non-zero common factor, to ensure that the sum thereof is equal to unity.

[0082] The purpose of step /5/ is to redistribute the state intervals as resulting from step /3/, to obtain a better statistical representativeness of the states that are possible for the aircraft 20. The implementation of this step /5/ may be subject to the result of an optional test designated by CR in FIG. 2. This test consists in determining whether a representativity criterion is satisfied for the state intervals with their respective weight values. It relates to all the values of the weights W.sub.i as updated at step /4/, where W.sub.i is the weight that is associated with the box_i state interval. A first representativeness criterion that may be used is the one known by the term N-effective criterion. This criterion is satisfied if (Σ.sub.i=1....sub.N w.sub.i.sup.2).sup.-1 < θ.sub.eff.Math.N, where θ.sub.eff is an adjustment parameter of the N-effective criterion, comprised between 0 and 1. Another representativity criterion that is also possible is the one that is known by the term entropic criterion: it is satisfied if log(N) + Σ.sub.i=1...N w.sub.i.Math.log(w.sub.i) > θ.sub.ent.Math.N, where θ.sub.ent is an adjustment parameter of the entropic criterion, comprised between 0 and 1. However, other representativity criteria that are also known to a person skilled in the art may also be used alternatively. Step /5/ of redistribution of the state intervals is then applied if the representativeness criterion is not satisfied.

[0083] Step /5/ first comprises determining a division number n.sub.i that is attributed to the box_i state interval, with the index i also numbering the state intervals from 1 to N. This sub-step is denoted /5-1/ in FIG. 2. In a known manner, it can be executed using a multinomial drawing method. Such a method consists in drawing points randomly and repetitively within a unidimensional segment, to determine the number n.sub.i of sub-intervals that are intended to replace the box_i state interval, according to the values of the weights w.sub.i of all the state intervals as updated at step /4/. The drawing segment is formed by a juxtaposition of base segments, the individual lengths of which correspond one-to-one to the values of the weights of the state intervals. The number n.sub.i of sub-intervals that will replace the box_i state interval is then proportional to the number of randomly drawn points that belong to the base segment whose the length is equal to the weight value of the box_i state interval. If the value of n.sub.i that is determined for one of the state intervals is zero, this state interval is suppressed, otherwise the box_i state interval is divided into n.sub.i sub-intervals, the value obtained for n.sub.i being statistically all the greater, the higher the weight of the box_i state interval as resulting from step /4/. Optionally, the values of n.sub.i that are thus obtained may be multiplied by a constant factor and rounded, in order to limit the total number of sub-intervals used, and/or to ensure that the sum of all the values of n.sub.i is equal to N. However, this way of implementing sub-step /5-1/ is provided only by way of example, and different methods may be used alternatively in variant embodiments of the invention.

[0084] Each box_i state interval is then intended to be divided into n.sub.i sub-intervals, at sub-step /5-3/, for example in one of the ways that are described in the following articles: “An introduction to box particle filtering [lecture notes]” by Gning, A., Ristic, B., Mihaylova, L., & Abdallah, F., IEEE Signal Processing Magazine (2013), Vol. 30(4), pp. 166-171; and “A Box Regularized Particle Filter for state estimation with severely ambiguous and non-linear measurements” by Merlinge, N., Dahia, K., Piet-Lahanier, H., Brusey, J., & Horri, N., Automatica (2019), Vol. 104, pp. 102-110, already cited above.

[0085] These division methods consist in first determining which one of the state coordinates of the system is that where the box_i state interval is the most extended. This determination is the object of sub-step /5-2/, which is now described.

[0086] Each of the state intervals box_i relates to a plurality of state coordinates of at least two different types. For example, in case of the navigation unit 10 that was described above for the aircraft 20, the three position coordinates, denoted x.sub.1, x.sub.2 and x.sub.3, are state coordinates of a first type, the three speeds, denoted v.sub.1, v.sub.2 and v.sub.3, are state coordinates of second type, and the three attitude angles, denoted θ.sub.1, θ.sub.2 and θ.sub.3, of the aircraft 20 are state coordinates of a third type. Thus normalized dimensionless values of the lengths of the unidimensional intervals of the box_i state interval, according to each of the state coordinates, are:

$\Delta {\text{x}}_{\text{j\_n}}=\Delta {\text{X}}_{\text{j}}\cdot {\left(\Delta {\text{X}}_1{}^2+\Delta {\text{X}}_2{}^2+\Delta {\text{X}}_3{}^2\right)}^{-1/2},$

$\Delta {\text{v}}_{\text{j\_n}}=\Delta {\text{V}}_{\text{j}}\cdot {\left(\Delta {\text{V}}_1{}^2+\Delta {\text{V}}_2{}^2+\Delta {\text{V}}_3{}^2\right)}^{-1/2},$

and

$\Delta {\theta }_{\text{j\_n}}=\Delta {\theta }_{\text{j}}\cdot {\left(\Delta {\theta }_1{}^2+\Delta {\theta }_2{}^2+{\theta }_3{}^2\right)}^{-1/2},$

for j=1, 2 and 3 in each case, where ∆x.sub.j, ∆v.sub.j and ∆θ.sub.j are lengths of the respective unidimensional intervals of the nine state coordinates for the box_i state interval. Each of the values ∆x.sub.j.sup.2, ∆v.sub.j.sup.2 and ∆θ.sub.j.sup.2 can be calculated using the method that was presented above for X.sup.α, with α equal to 2. Then each of the factors (∆x.sub.1.sup.2+∆x.sub.2.sup.2+∆x.sub.3.sup.2).sup.-½, (∆v.sub.1.sup.2+∆v.sub.2.sup.2+∆v.sub.3.sup.2).sup.-½ and ∆θ.sub.1.sup.2+∆θ.sub.2.sup.2+∆.sub.3.sup.2).sup.-½ can also be calculated using the same method, but with α then being equal to -½. The normalized values ∆x.sub.j_n, ∆v.sub.j_n and Vθ.sub.j_n of the lengths of the new unidimensional intervals of the box_i state interval are then obtained in accordance with previous formulae, by calculation of products.

[0087] These normalized values ∆x.sub.j_n, ∆v.sub.j_n and ∆θ.sub.j_n can be compared with each other, and the state coordinate where the box_i state interval is the most extended is the one that corresponds to the highest of the normalized values ∆x.sub.j_n, ∆v.sub.j_n and ∆θ.sub.j_n, considering their absolute values. For example, the box_i state interval is the most extended according to the position coordinate x.sub.1 if ∆x.sub.1_n is the highest of the normalized values ∆x.sub.j_n, ∆v.sub.j_n and ∆θ.sub.j_n, or it is the most extended according to the speed coordinate v.sub.2 if ∆v.sub.2­_n is the highest of the nine normalized values, etc.

[0088] At sub-step /5-3/, the box_i state interval is divided into n.sub.i contiguous sub-intervals, by dividing, into n.sub.i segments with the same lengths, that one of the unidimensional intervals of box_i that is the most extended, within the meaning of the normalized values ∆x.sub.j_n, ∆v.sub.j_n and ∆θ.sub.j_n. Each sub-interval therefore has one of these segments as a unidimensional interval along the state coordinate that corresponds to the maximum extension of the box_i state interval, and the same unidimensional intervals as this box_i state interval along the other state coordinates. The n.sub.i sub-intervals that are thus constructed therefore constitute a partitioning of the box_i state interval. They form new state intervals from which the box-regularized particle filtering is continued. A weight value is attributed to each of them, which is equal to that of the box_i state interval divided by the division number n.sub.i. These new state intervals, with those of the state intervals that have not been divided, are then re-numbered by the index i, advantageously from 1 to N, to continue the filtering method.

[0089] Finally, sub-step /5-4/, which was referred to as smoothing step in the general part of the present description, consists in correcting the new state intervals as resulting from sub-step /5-3/, so that they produce, with their associated weight values, an even better statistical representation of the state of the system, i.e. of the state of the aircraft 20 for the example in question. Such a modification of the state intervals is also commonly referred to as regularization of the statistical representation of the state of the system, in the jargon of a person skilled in the art. It can be applied not only to the new state intervals that each resulted from a division of one of the subsequent state intervals, but also to all the state intervals, including those that have not been divided. Corrections that are applied for this to the state intervals may consist of random shifts of the bounds of the unidimensional intervals that constitute the edges of each box_i state interval. Preferably, the respective weight values that are associated with the state intervals are not modified in this sub-step /5-4/. For the invention, a smoothing method by Epanechnikov kernel is applied, in particular as described in the article by Merlinge et al., Automatica 104(2019), already cited.

[0090] In a known manner, the Epanechnikov kernel is defined by the probability density function f(x) = 3.square-solid.(1 - x.sup.2)/4, where x is the random variable comprised between -1 and 1, the values -1 and 1 being permitted. Its average value is zero, and its variance is equal to ⅕.

[0091] Sub-step /5-4/ therefore initially comprises generating random corrections to be applied to each unidimensional state-coordinate interval that each determines box_i state intervals, and then applying these corrections. A detailed implementation of sub-step /5-4/ is shown in FIG. 3.

[0092] The number of state coordinates of the system in question is denoted d, and the number of bounds that determine each state interval of this system, i.e. each box used in the regularized particle filtering method, is denoted d‘=2.Math.d. In the case of the aircraft 20, d=9 and d’=18.

[0093] N first random values are first generated, denoted β.sub.i and comprised between 0 and 1, the values 0 and 1 being permitted, i being the integer index of numbering of the state intervals as used previously, each according to a statistical distribution beta law, with first parameter equal to d and second parameter equal to 2. In a known manner, the beta law the parameters of which are d and 2 is defined by the probability density function x.sup.d/2-1.square-solid.(1-x).square-solid.Γ(d/2+2)/[Γ(d/2).square-solid.Γ(2)], where x designates the random variable, and Γ designates the gamma function. One possible way of generating the values β.sub.i in accordance with the beta law uses Cheng’s algorithm, which will be reminded below, with reference to FIG. 4.

[0094] Next N.Math.d‘ second random values are generated, denoted v.sub.k,i, i also being the same index previously and k being another integer index that varies from 1 to d’, each according to a normal statistical distribution law with zero mean value and standard deviation equal to unity, commonly referred to as reduced normal law. The index k counts the degrees of freedom in the definition of each state interval. It identifies two unidimensional interval bounds for each state coordinate or, in an equivalent manner, a center value and an interval length for each unidimensional state-coordinate interval. In a known manner, the normal law with zero mean value and standard deviation equal to unity is defined by the probability density function (⅟Π.sup.½).square-solid.exp[-x.sup.2/2], where x also designates the random variable, but in this case positive, zero or negative. It can be simulated by a sum of initial random values that are each generated according to a uniform statistical distribution law, and such a uniform statistical distribution law can be produced by a method of the LFSR type, standing for “linear feedback shift register”, for example.

[0095] N first numbers, denoted ξ.sub.i, are then calculated from the random values v.sub.k,i in the following manner: ξ.sub.i = [Σ.sub.k=1 .sub.à .sub.d′ (v.sub.k,i).sup.2].sup.½. Advantageously, each of these numbers ξ.sub.i can be calculated by applying the method for calculating an estimation of x.sup.α that was described above, to X = v.sub.k,i with α=2, then to X = Σ.sub.k=1 .sub.à .sub.d′ (v.sub.k,i).sup.2 with α=½.

[0096] N second numbers, denoted α.sub.i, are next calculated according to the formula: α.sub.i = (β.sub.i.sup.½/ξ.sub.i. Advantageously, each of these numbers α.sub.i can be calculated by once again applying the method for calculating an estimation of x.sup.α that was described above, to X = β.sub.i with α=½.

[0097] Under these conditions, N.square-solid.d′ third numbers that are calculated in the following manner: ε.sub.k,i = v.sub.k,i.square-solid.α.sub.i, each satisfy the Epanechnikov statistical distribution law of zero average value and of variance equal to ⅕.

[0098] Moreover, a noise amplitude, denoted h, is calculated in the following manner:

$\text{h}\text{=}\mu \cdot \text{A}\cdot {\text{N}}^{-1/(\text{d′+4)}},$

where A = [8.square-solid.c.sub.d.sup.,-1.square-solid.(d‘+4).square-solid.(2.square-solid.Π.sup.½).sup.d′].sup.⅟(d′+4). In this expression of A, c.sub.d′ designates the volume of the hypersphere of dimension d’ and of unitary radius, and .Math. is an adjustment parameter that is comprised between 0 et 1. In a known manner, c.sub.d′ = Π.sup.d′/2/Γ(d′/2 + 1), where Γ again designates the gamma function. The value of c.sub.d′ can either be precalculated and stored so as to be available for the computing unit 3, or be calculated by the latter, for example using a prerecorded table of values of the gamma function. The adjustment parameter .Math. makes it possible to control a trade-off between the efficiency of the random smoothing and efficacy of the particle filter. This is because the efficiency of the particle filter results from a continuity of the possible trajectories for the system, described by the trajectories of the state intervals that are maintained during a plurality of successive implementations of the sequence of steps /1/ to /5/. On the other hand, random smoothing produces jamming of these trajectories. The value of the parameter .Math. can be fixed initially for the method to be executed by the computing unit 3. It depends in particular on the numbers N and d. For example, the adjustment parameter .Math. may be taken equal to 0.1. The use of the adjustment parameter .Math. is in particular described in the thesis of Merlinge, N., entitled “State estimation and trajectory planning using box particle kernels”, Paris-Saclay University, 2018.

[0099] The random modifications to be applied to the box_i state interval are then h.square-solid.ε.sub.k,i, where k identifies the state coordinate bounds, or the center values or lengths of state coordinate intervals. These modifications can be arranged in the form of a vector E.sub.i, such that E.sub.i = [h.square-solid.ε.sub.k,i].sub.k=1,...,d′.

[0100] A possible method for applying the random modifications to the state intervals that result from sub-step /5-3/ is now described. For this purpose, it is possible to represent each box_i state interval by a vector of height d‘, the coordinates of which contain the bounds of all the unidimensional state coordinate intervals for this state interval, or the center values and the lengths of its unidimensional intervals. Then the randomly modified values of the vector can be obtained by replacing this vector by + D x E.sub.i, where D is a square matrix of dimension d’ such that the product of D multiplied by the transpose of D is equal to .Math..sub.i=1 .sub.à .sub.N

${\Xi }_{\text{i}}\cdot {\text{w}}_{\text{i}}{\cdot }^{\text{t}}{\Xi }_{\text{i}}$

: D x .sup.tD = .Math..sub.i=1 .sub.à .sub.N

${\Xi }_{\text{i}}\cdot {\text{w}}_{\text{i}}{\cdot }^{\text{t}}{\Xi }_{\text{i}},$

where w.sub.i is the weight value of the box_i state interval. For example, the matrix D can be determined by Cholesky’s method, which is well-known to a person skilled in the art so that it is not necessary to describe it again here. Under these circumstances, the vector

${\Xi }_{\text{i}}$

is to be replaced by the vector

${\Xi }_{\text{i}}$

+ D x E.sub.i, for applying the random variations according to the Epanechnikov kernel to the unidimensional intervals of state coordinates of the box_i state interval. All the vectors + D x E.sub.i determine all the new state intervals that result from the modification by smoothing. In the general part of the present description, the square matrix .Math..sub.i=1 .sub.to .sub.N

${\Xi }_{\text{i}}\cdot {\text{w}}_{\text{i}}{\cdot }^{\text{t}}{\Xi }_{\text{i}},$

of dimension d′ has been called the mean product matrix. Its coefficients are the mean values calculated over all the state intervals in question, of products of interval bounds, or center values or interval lengths, taken in pairs separately for each state interval.

[0101] The method that has just been described for applying, at sub-step /5-4/, a smoothing by Epanechnikov kernel to the sub-intervals as resulting from sub-step /5-3/, is implemented by the computer unit 3 of the navigation system using measurement correlation 10.

[0102] The updated state intervals, as resulting from the sub-step /5-3/ or from sub-steps /5-4/, with their associated weight values, constitute the result of the box-regularized particle filtering method for an implementation of the sequence of steps /1/ to /5/. A plurality of iterations are chained recurrently, each new iteration from the results of the previous iteration. The iteration that is implemented last supplies an updated probability distribution that characterizes the true state of the aircraft 20. Furthermore, the updated state intervals that result therefrom, associated with their respective weight values, are to be used as prior state intervals for a new implementation of the sequence of steps /1/ to /5/.

[0103] A description is now given, with reference to FIG. 4, of a method for generating each random value β.sub.i, so as to satisfy the Beta(d, 2) statistical distribution law. This method corresponds to Cheng’s algorithm, published in the article entitled “Generating beta variates with non-integral shape parameters”, Communications of the ACM, 21(4), pp. 317-322, 1978.

[0104] For this purpose, the following numbers a and b are first determined at step /i/: [0105] a, which is equal to the minimum value between d and 2: a = min(d, 2), and [0106] b which is equal to the maximum value between d and 2: b = max(d, 2).

[0107] For most of the applications of the invention, d is greater than 2, so that a=2 and b=d. Next, at step /ii/, the three numbers α, β and γ are calculated such that:

$\alpha =\text{a + b, which is equal to d+2,}$

$\beta ={\left[\left(\alpha -2\right)/\left(2\cdot \text{a}\cdot \text{b -}\alpha \right)\right]}^{1/2},\text{which is equal to}{\left[\text{d/}\left(3\cdot \text{d - 2}\right)\right]}^{1/2};$

and

$\gamma =\text{a + 1/}\beta$

[0108] The numbers α, β and γ may have been precalculated and stored in order to be directly accessible to the computing unit 3.

[0109] At the step /iii/, two numbers u.sub.1 and u.sub.2 are generated randomly each according to the uniform statistical distribution law, between 0 and 1, for example using a method of the LFSR type. Then the numbers V, W, Z, R and S are calculated at step /iv/, according to the following formulae: [0110] V = β.square-solid.log[u.sub.1/(1-u.sub.1)], where log designates the base-e logarithmic function, [0111] W = a.square-solid.exp(V), where exp designates the base-e exponential function, [0112] Z = u.sub.1.sup.2.square-solid.u.sub.2, [0113] R = γ.square-solid.V - log(4), and [0114] S = a + R - W.

[0115] The value of the number Z can advantageously be calculated each time by applying the calculation method of an estimation of x.sup.α that was described above, to X = u.sub.1 with α=2. Furthermore, the values of the logarithm and exponential functions can be obtained from tables of values prerecorded for these functions.

[0116] Steps /v/ to /viii/ are then implemented successively, forming a sequence of tests that are applied in series: [0117] step /v/: if S + 1 + log(5) is greater than or equal to 5.square-solid.Z, go directly to step /viii/; [0118] step /vi/: if S is greater than or equal to log(Z), go directly to step /viii/; [0119] step /vii/: if R + α.square-solid.[log[α/(b+W)] is less than log(Z), go back to step /iii/; [0120] step /viii/: if a is equal to d, then the value β.sub.i that is generated randomly according to the beta law of parameters d and 2, is equal to W/(b+W), and if a is different from d, it is equal to b/(b+W). For most of the applications of the invention, where d is greater than 2, the random value β.sub.i is equal to d/(d+W).

[0121] One of the advantages of this algorithm is that the number of implementations of the sequence of steps /iii/ to /viii/ that is necessary for obtaining the N random values β.sub.i is predictable.

[0122] It is clear that the invention can be reproduced by modifying secondary aspects of the implementations that have been described in detail above, while keeping at least some of the advantages mentioned. Among such possible modifications, the following are cited non-limitedly: [0123] alternative algorithms may be used for some steps or sub-steps; [0124] the inertial unit may be used as a source for measuring the complete state of the vehicle. The state of the system as considered for the invention then comprises three acceleration values and three angular velocity values, according to the three spatial coordinates, in addition to the three values of spatial coordinates that identify a position for the vehicle, three speed values and three angular attitude values. In this case, each state of the system comprises values for fifteen state coordinates; [0125] in step /⅟, a dynamic state-evolution model may be used for the system, this model being able to take account, when it is a vehicle, of commands that are applied to one or more engines and to an attitude control system of the vehicle; and [0126] the computing unit may be a processor with an ARM structure, a processor with several computing cores, one or more graphics processors, etc., in place of a chip of the FPGA, DSP, CPU or RISC type.

[0127] Finally, the invention may be applied to fields other than aeronautics. For example, a box-regularized particle filtering method that is in accordance with the invention may also be implemented for a vehicle on the ground, a surface seagoing ship, a submarine, a satellite or a space probe, using in each case a reference frame and true state measurements that are suitable.