PARTICLE FILTERING METHOD AND NAVIGATION SYSTEM USING MEASUREMENT CORRELATION

Abstract

A box regularized particle filtering method implements a binary representation of numbers. This implementation can be used to determine a box division coordinate and/or to modify state intervals according to a fixed probability kernel, for example according to an Epanechnikov kernel. The method can be executed autonomously within a navigation system using measurement correlation, in particular on board an aircraft such as an airplane, a flying drone or any self-propelled airborne vehicle.

Claims

1. A box regularised particle filtering method for predicting a state of a system, said system being a land, air, sea or space vehicle which is equipped with a navigation system using measurement correlation, the method comprising at least one execution of a sequence which includes the following steps: /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 previous 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 measurement result of the true state which was obtained in step /2/; /4/ a step of updating weights, comprising assigning a weight to each subsequent state interval according to a size of said subsequent state interval as resulting from step /3/, with a size of said subsequent state interval as resulting from step /1/ before step /3/, and a weight of the previous state interval from which said subsequent state interval resulted in step /1/; and /5/ a step of redistributing the subsequent state intervals, comprising replacing at least one of the subsequent state intervals i by n.sub.i sub-intervals originating from a division of the subsequent state interval i, i being an integer numbering index of the subsequent state intervals, and n.sub.i being a number of the sub-intervals which replace said subsequent state interval i, each sub-interval being intended to constitute a previous state interval for a subsequent iteration of the sequence of steps /1/ to /5/, wherein at least one of the steps /1/ to /5/ comprises at least one calculation of an estimate of the value of X.sup.α, where X is a variable number for each execution of said step and α is a non-zero exponent value that is constant between different executions of said step, the number X being positive if the exponent α is not an integer and non-zero if the exponent α is negative, the calculation comprising the following steps: /a/ writing the number X in the form X=±(1+m).Math.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 being allowed, so that a binary representation of the number X is: I(X)=L.Math.(m+ex+B), where L=2.sup.n with n which is a number of bits of a binary writing of the mantissa m, and B is a positive or zero constant number, called 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 whose value is recorded; and /c/ obtaining the estimate of the value of X.sup.α from the binary representation I(X.sup.α).

2. The method according to claim 1, wherein the step of the box regularised particle filtering method which comprises the calculation of the estimate of the value of X.sup.α is step /5/.

3. The method according to claim 1 or 2, wherein the calculation of the estimate of the value of X.sup.α further comprises executing at least once the following additional step: /d/ calculating a new estimate of the value of X.sup.α from a previous estimate of the value of X.sup.α, by applying a recursive algorithm for approximate equation solving to the equation Y.sup.1/α−X=0 of unknown Y, the estimate of the value of X.sup.α which was obtained in step /c/ being used as a previous estimate for a first application of said algorithm, and the new estimate of the value of X.sup.α which is produced by a q.sup.th application of the algorithm forming the previous estimate of the value of X.sup.α for the (q+1).sup.th application of said algorithm, if such a (q+1).sup.th application of the algorithm is performed, q being an integer number greater than or equal to 1.

4. The method according to claim 3, wherein the recursive algorithm for approximate equation solving that is used in step /d/ is Newton's method.

5. The method according to claim 1, wherein the number n of the bits of the binary writing of the mantissa m is equal to 23, and the bias B is equal to 127, or the number n is equal to 52 bits and the bias B is equal to 1023.

6. The method according to claim 1, wherein the constant number σ is comprised between 0 and 1.

7. The method according to claim 1, wherein the exponent value α is equal to +2, −2, +½ or −½.

8. A calculation unit, comprising at least one first input adapted to receive results of repeated measurements of acceleration and angular speed of a system, and a second input adapted to receive results of repeated measurements of a true state of the system, additional to the acceleration and angular speed measurements, and the calculation unit being arranged to execute a box regularised particle filtering method which 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 value that the true state of the system is within said state interval.

9. The calculation unit according to claim 8, of field-programmable gate array circuit, fixed gate array circuit, or central unit processor type.

10. A navigation system using measurement correlation, adapted to be on board a vehicle, comprising: an inertial system, adapted to iteratively measure accelerations and angular speeds of the vehicle, and to deduce, using the results of measurements of the accelerations and angular speeds, subsequent state intervals respectively from several previous state intervals, each state of the vehicle comprising position, speed and attitude coordinates of said vehicle; a measurement system, adapted to iteratively measure at least one feature of a true state of the vehicle; and a calculation unit according to claim 8, and adapted to reduce at least one drift of the inertial system, among a position drift, a speed drift and an attitude drift, by using measurement results of the feature of the true state of the vehicle which are delivered by the measurement system.

11. A vehicle, comprising a navigation system using measurement correlation which is in accordance with claim 10.

12. The method according to claim 1, wherein the step of the box regularised particle filtering method which comprises the calculation of the estimate of the value of X.sup.α is a sub-step of step /5/ which consists in calculating a measure value of a part of a subsequent state interval of the system, corresponding to a position interval of the system, a speed interval of said system or an attitude interval of said system, and/or another sub-step of step /5/ which consists in applying a probability kernel smoothing to the sub-intervals to obtain the previous state intervals of the subsequent iteration of the sequence of steps /1/ to /5/.

13. The method according to claim 1, wherein the constant number σ is comprised between 0 and 0.5.

14. The method according to claim 2, wherein the calculation of the estimate of the value of X.sup.α further comprises executing at least once the following additional step: /d/ calculating a new estimate of the value of X.sup.α from a previous estimate of the value of X.sup.α, by applying a recursive algorithm for approximate equation solving to the equation Y.sup.1/α−X=0 of unknown Y, the estimate of the value of X.sup.α which was obtained in step /c/ being used as a previous estimate for a first application of said algorithm, and the new estimate of the value of X.sup.α which is produced by a q.sup.th application of the algorithm forming the previous estimate of the value of X.sup.α for the (q+1).sup.th application of said algorithm, if such a (q+1).sup.th application of the algorithm is performed, q being an integer number greater than or equal to 1.

15. The method according to claim 2, wherein the number n of the bits of the binary writing of the mantissa m is equal to 23, and the bias B is equal to 127, or the number n is equal to 52 bits and the bias B is equal to 1023.

16. The method according to claim 3, wherein the number n of the bits of the binary writing of the mantissa m is equal to 23, and the bias B is equal to 127, or the number n is equal to 52 bits and the bias B is equal to 1023.

17. The method according to claim 4, wherein the number n of the bits of the binary writing of the mantissa m is equal to 23, and the bias B is equal to 127, or the number n is equal to 52 bits and the bias B is equal to 1023.

18. The method according to claim 2, wherein the constant number σ is comprised between 0 and 1.

19. The method according to claim 3, wherein the constant number σ is comprised between 0 and 1.

20. The method according to claim 4, wherein the constant number σ is comprised between 0 and 1.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0039] The features and advantages of the present invention will appear more clearly in the detailed description below of non-limiting examples of implementation, with reference to the appended figures, among which:

[0040] FIG. 1 shows an aircraft which is equipped with an inertial system using terrain correlation in accordance with the invention;

[0041] FIG. 2 is a diagram which shows the sequence of steps of a method in accordance with the invention; and

[0042] FIG. 3 is a flowchart of an algorithm that can be used to implement the invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0043] For reasons of clarity, the dimensions of the elements which are represented symbolically in [FIG. 1] correspond neither to actual dimensions nor to actual dimensional ratios. Furthermore, the invention is described by way of non-limiting example for a case of application to an aircraft, but it is understood that it can be applied to any vehicle which is provided with a navigation system using measurement correlation, whether this vehicle is a land, air, sea, space vehicle, etc.

[0044] The mode of calculating an estimate of the value of X.sup.α which is used in the invention is described first. α denotes an exponent which is non-zero, whose value is constant. X is a number of variable value, which is positive when the value of the exponent α is not an integer, and non-zero when the value of the exponent α is negative.

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

X=±(1+m).Math.2.sup.ex,

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

[0046] Then, when X is positive, we have:

log.sub.2X=ex+m+σ,

where σ is a fixed real number allowing to minimise an error on the value of log.sub.2X, in particular when a belonging numerical interval is known a priori for the number X. For example, the value of the number σ can be taken equal to 0.043036. And therefore:

ex+m=log.sub.2X−σ.

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

I(X)=L.Math.(m+ex+B),

where L=2.sup.n, with n being a fixed number of bits to write the mantissa m in binary form, and B being a positive or zero constant number, which is called bias. The role of the bias B is to allow the use of a representation of the type I(X) also for negative values of the number X. 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. To have a calculation precision that is higher, n may alternatively be equal to 52, and in this case B may be equal to 1023. By transferring in the binary representation I(X) the expression of ex+m as derived from log.sub.2X, it becomes:

I(X)=L.Math.[log.sub.2X−σ+B], that is to say: log.sub.2X=I(X)/L+σ−B.

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

I(X.sup.α)=L.Math.[log.sub.2(X.sup.α)−σ+B]

But log.sub.2(X.sup.α)=α.Math.log.sub.2(X), therefore: I(X.sup.α)=L.Math.[α.Math.log.sub.2(X)−σ+B], and replacing log.sub.2(X) by its expression according to the binary representation I(X), it becomes:

I(X.sup.α)=L.Math.[α.Math.(I(X)/L+σ−B)−σ+B], that is to say: I(X.sup.α)=α.Math.I(X)+L.Math.(1−α).Math.(B−σ).

An approximate value of X.sup.α, denoted Y.sub.0, can then be reconstructed from the binary representation of X.sup.α which was thus obtained, using a reverse method to that which provides the binary representation of a number from this number. The difference between this approximate value Y0 and the true value of X.sup.α depends on the value that has been adopted for the number σ. When X can be negative and the value of the exponent α is integer, the sign of X.sup.α can be determined according to that of X, retained since the start of the calculation method, and the parity of the exponent α. For many applications, the approximate value Y.sub.0 is satisfactorily suitable for replacing X.sup.α, given the simplicity of the method for obtaining this approximate value Y.sub.0, as just described.

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

I(X.sup.2)=2.Math.I(X)−L.Math.(B−σ).

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

I(X.sup.−2)=−2.Math.I(X)+3.Math.L.Math.(B−σ).

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

I(X.sup.1/2)=I(X)/2+0.5.Math.L.Math.(B−σ).

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

I(X.sup.−1/2)=−I(X)/2+1.5.Math.L.Math.(B−σ).

[0053] 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 which are known to the person skilled in the art. Newton's algorithm, also called Newton's method, may be used in particular, by applying it to the function f(Y)=Y.sup.1/α−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 to say, by calculating the expression of f′(Y) from that of f(Y):

Y.sub.q+1=(1−α).Math.Y.sub.q+α.Math.X.Math.Y.sub.q.sup.(α−1)α, for q=0,1,2, . . . .

[0054] For the particular case where the exponent α is equal to 2, it becomes:

Y.sub.q+1=−Y.sub.q+2.Math.X.Math.Y.sub.q.sup.1/2.

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

Y.sub.1=−Y.sub.0+2.Math.X.Math.Y.sub.0.sup.1/2.

The value of Y.sub.q.sup.1/2 can be estimated each time by 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 a by ½.

[0055] For the particular case where the exponent α is equal to −2, it becomes:

Y.sub.q+1=3.Math.Y.sub.q−2.Math.X.Math.Y.sub.q.sup.3/2.

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

Y.sub.1=3.Math.Y.sub.0−2.Math.X.Math.Y.sub.0.sup.3/2.

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

[0056] For the particular case where the exponent α is equal to ½, it becomes:

Y.sub.q+1=0.5Y.sub.q+0.5X/Y.sub.q.

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

Y.sub.1=0.5.Math.Y.sub.0+0.5.Math.X/Y.sub.0.

[0057] For the particular case where the exponent α is equal to −½, it becomes:

Y.sub.q+i=(3/2)Y.sub.q−0.5.Math.X.Math.Y.sub.q.sup.3.

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

Y.sub.1=(3/2).Math.Y.sub.0−0.5.Math.X.Math.Y.sub.0.sup.3.

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

[0059] In accordance with [FIG. 1], an aircraft 20 is equipped with a navigation system using terrain correlation, designated by the reference 10. The navigation system 10 comprises an inertial system 1, a telemetry probe 2 and a calculation unit 3. In a known manner, the inertial system 1 repeatedly performs measurements of three acceleration coordinates and three angular speed coordinates of the aircraft 20, using accelerometers and gyrometres, not shown. Moreover, the telemetry probe 2 repeatedly performs measurements of the distance H which exists between the ground 100 and the aircraft 20. This distance H is measured in a direction which may be fixed with respect to the aircraft 20, for example perpendicularly to a reference plane of the aircraft. Optionally, the telemetric measurement direction can be variable with respect to the aircraft 20, but in such a case this measurement direction is controlled and taken into account in an adequate manner, which is also known. The distance H therefore varies according to the terrestrial relief which is overflown by the aircraft 20, as well as according to its altitude and its attitude. In the following, the present description is limited to the only case where the measurements additional to those of the inertial system 1 consist of the measurements of the distance H. This case corresponds to a terrain correlation navigation. However, it is understood that other additional measurements can be used alternatively to those of the distance H, or in addition to them. Finally, the invention is compatible with inertial system and telemetry probe models as commercially available. In particular, the inertial system can be of a MEMS, for “Micro-Electro-Mechanical System”, quadrason, gyrolaser type, etc., and the telemetry probe can be a radio altimeter, a laser telemeter, etc.

[0060] Each possible state for the aircraft 20 can be composed of three spatial coordinate values which identify a position for the aircraft, for example in the terrestrial reference frame, three speed values each according to one of the spatial coordinates, and three angular values to identify 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 one-dimensional intervals which each relate separately to one of the state coordinates. Such a state interval is called a box in the language of the person skilled in the art.

[0061] The box regularised particle filtering method is initialised by supplying several initial state intervals, each associated with a weight which indicates a probability value for the true state of the aircraft 20 to initially be in this initial state interval. Thus, N initial state intervals are provided, N being an integer number, 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 can be equal to 1/N.

[0062] The method then consists of successive iterations of a sequence of steps, each new execution of the sequence of steps producing an update of the state intervals, with updated weight values which are associated one-by-one with the updated state intervals. Furthermore, each new execution of the sequence of steps is performed from the state intervals and their associated weight values as provided by the immediately preceding execution of the sequence of steps.

[0063] 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 redistributing the state intervals. Steps /1/ and /3/ to /5/ are executed for each state interval. Since step /5/ produces a redistribution of the state intervals as resulting from step /3/, it is preferably executed so as to keep a constant number of state intervals. Then the calculation unit 3 may be designed and dimensioned to process N state intervals at each execution of the sequence of steps /1/ to /5/. i is an integer index, from 1 to N, which numbers the state intervals that are processed at each iteration of this sequence of steps /1/ to /5/.

[0064] For each state interval i, denoted box_i, the prediction step /1/ consists in collecting the results of the last acceleration and angular speed measurements, as delivered by the inertial system 1. Optionally, the results of several measurements which have been carried out by the inertial system 1 since the previous execution of the sequence of steps /1/ to /5/can be collected. For each of the state intervals box_i, this step /1/ also comprises calculating a change of this state interval during the period of time which has elapsed between the two executions of step /1/, relating to the previous iteration and the current new iteration of the sequence of steps /1/ to /5/. The calculation principle of such change in each state interval box_i, from the results of the acceleration and angular speed measurements, is assumed to be known to the person skilled in the art. Reference may be made on this subject 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, most often accompanied by a change in length, of each one-dimensional interval which relates to one of the state coordinates of the aircraft 20. In the general part of this description, each state interval box_i as it exists at the time of starting the execution of step /1/ has been called a previous state interval, and this state interval box_i as modified by step /1/ has been called the subsequent state interval.

[0065] The measurement step /2/ consists in collecting the result of the last measurement of distance H as delivered by the telemetry probe 2. Possibly, the results of several last measurements which have been carried out by the telemetry probe 2 since the previous execution of step /2/, during the previous iteration of the sequence of steps /1/ to /5/, can be collected.

[0066] The contraction step /3/ consists in reducing the size of each state interval box_i according to incompatibilities which would exist between parts of this state interval and the last distance H measurement result(s) collected in step /2/. Typically, the six position and attitude coordinates of the aircraft 20, as may vary within each state interval box_i resulting from step /1/, are combined with relief height values read in a terrain map which is stored, in order to obtain an estimate of the distance H associated with each state. Such a calculation can be performed in the known manner which was recalled at the beginning 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 implemented to convert each state of the aircraft 20 into a distance H value, using the stored terrain map. But such an equation, called an observation equation by the person skilled in the art, can be difficult to invert locally, so that modelling its inverse function by an analytical function part and a tabulated function part can be advantageous. The person skilled in the art can also refer to the article by Merlinge and al., Automatica 104 (2019), mentioned above, about step /3/ which is not directly concerned with improving the method provided by the invention.

[0067] The update 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 state interval box_i may be equal to the value of the weight of this state interval box_i as it existed before performing this update, multiplied by the size of the state interval box_i as resulting from the contraction step /3/, and divided by the size of the state interval box_i as resulting from the prediction step /1/ before applying the contraction step /3/. Other formulae for updating weight values may be adopted alternatively. Optionally, the weight values as resulting from one of these formulas can also be corrected, for example by multiplying them by a non-zero common factor, to ensure that their sum is equal to the unit.

[0068] 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 which are possible for the aircraft 20. The execution of this step /5/ may be subjected to the result of an optional test designated by CR in [FIG. 2]. This test consists in determining if a representativeness criterion is satisfied by the state intervals with their respective weight values. It relates to all the values of the weights w.sub.i as updated in step /4/, where w.sub.i is the weight that is associated with the state interval box_i. A first representativeness criterion that can be used is the one known as the N-effective criterion. This criterion is satisfied if (Σ.sub.i=1 . . . 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. Also another possible representativeness criterion is the one known as the 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. But other representativeness criteria which are also known to the person skilled in the art can still be used alternatively. Step /5/ of redistributing state intervals is then applied if the representativeness criterion is not satisfied.

[0069] Step /5/ comprises first determining a division number n.sub.i which is assigned to the state interval box_i, with the index i which again numbers 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 repeatedly inside a one-dimensional segment, to determine the number n.sub.i of the sub-intervals which are intended to replace the state interval box_i, according to the values of the weights w.sub.i of all state intervals as updated in step /4/. The drawing segment is made up of a juxtaposition of elementary segments, the individual lengths of which correspond one-by-one to the values of the weights of the state intervals. The number n.sub.i of sub-intervals which will replace the state interval box_i is then proportional to the number of randomly drawn points which belong to the elementary segment whose length is equal to the weight value of the state interval box_i. If the value of n.sub.i which is thus determined for one of the state intervals is zero, this state interval is deleted, otherwise the state interval box_i is divided into n.sub.i sub-intervals, the value obtained for n.sub.i being statistically greater as the weight of the state interval box_i as resulting from step /4/ is high. Optionally, the values of n.sub.i which are thus obtained can be multiplied by a constant factor and rounded, in order to limit the total number of sub-intervals implemented, and/or to ensure that the sum of all the values of n.sub.i is equal to N. However, this way of executing sub-step /5-1/ is only given by way of example, and different methods can be used alternatively in variant implementations of the invention.

[0070] Each state interval box_i is then intended to be divided into n.sub.i sub-intervals, in sub-step /5-3/, for example in one of the ways that are described in the following articles: [0071] “An introduction to box particle filtering [reading notes]” by Gning, A., Ristic, B., Mihaylova, L., & Abdallah, F., IEEE Signal Processing Magazine (2013), Vol. 30(4), pp. 166-171; and [0072] “A Box regularised 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 mentioned above.

[0073] These division methods consist in first determining that of the system state coordinates, according to which the state interval box_i is the largest. This determination is the subject of sub-step /5-2/ which is now described.

[0074] Each of the state intervals box_i relates to several state coordinates of at least two different types. For example, in the case of the navigation system 10 which 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 a 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. Then, normalised and dimensionless values of the lengths of the one-dimensional intervals of the state interval box_i, according to each of the state coordinates, are:

Δx.sub.j_n=Δx.sub.j.Math.(Δx.sub.1.sup.2.Math.Δx.sub.2.sup.2+Δx.sub.3.sup.2).sup.−1/2,

Δv.sub.j_n=Δv.sub.j.Math.(Δv.sub.1.sup.2+Δv.sub.2.sup.2+Δv.sub.3.sup.2).sup.−1/2, and

Δθ.sub.j_n=Δθ.sub.j.Math.(Δθ.sub.1.sup.2+Δθ.sub.2.sup.2+Δθ.sub.3.sup.2).sup.−1/2,

for j=1, 2 and 3 in each case, where Δx.sub.j, Δv.sub.j and Δθ.sub.j are the lengths of the respective one-dimensional intervals of the nine state coordinates for the state interval box_i. 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.−1/2, (Δv.sub.1.sup.2+Δv.sub.2.sup.2+Δv.sub.3.sup.2).sup.−1/2 and (Δθ.sub.1.sup.2+Δθ.sub.2.sup.2+Δθ.sub.3.sup.2).sup.−1/2 can also be calculated using the same method, but with a which is then equal to −1/2. These factors have been referred to as measure values of state interval parts in the general part of this description. The normalised values Δx.sub.j_n, Δv.sub.j_n and Δθ.sub.j_n of the lengths of the nine one-dimensional intervals of the state interval box_i are then obtained according to the previous formulas, by calculating products.

[0075] These normalised values Δx.sub.j_n, Δv.sub.j_n and Δθ.sub.j_n can be compared with each other, and the state coordinate according to which the state interval box_i is the largest is that which corresponds to the largest of the normalised values Δx.sub.j_n, Δv.sub.j_n and Δθ.sub.j_n, considering their absolute values. For example, the state interval box_i is the largest according to the position coordinate x.sub.1 if Δx.sub.1_n is the largest of the normalised values Δx.sub.j_n, Δv.sub.j_n and Δθ.sub.j_n, or it is the largest according to the coordinate speed v.sub.2 if Δv.sub.2_n is the greatest of the nine normalised values, etc.

[0076] In sub-step /5-3/, the state interval box_i is divided into n contiguous sub-intervals, by dividing into n.sub.i segments of the same lengths that of the one-dimensional intervals of box_i which is the largest, within the meaning of the normalised values Δx.sub.j_n, Δv.sub.j_n and Δθ.sub.j_n. Each sub-interval therefore has one of these segments as one-dimensional interval according to the state coordinate which corresponds to the maximum extension of the state interval box_i, and the same one-dimensional intervals as this state interval box_i according to other state coordinates. The n.sub.i sub-intervals which are thus constructed therefore constitute a partition of the state interval box_i. They form new state intervals from which the box regularised particle filtering method is continued. A weight value is assigned to each of them, which is equal to that of the state interval box_i divided by the division number n.sub.i. These new state intervals are then re-numbered by the index i, advantageously from 1 to N, for the continuation of the filtering method.

[0077] Finally, the optional sub-step /5-4/ consists in correcting the new state intervals as resulting from the sub-step /5-3/, so that they produce, with their values of associated weights, an even better statistical representation of the state of the system, that is to say of the state of the aircraft 20 for the considered example. Such a modification of the state intervals is commonly called smoothing or regularising the statistical representation of the state of the system, in the language of the person skilled in the art. Corrections which are applied for this to the state intervals can consist of random displacements of the limits of the one-dimensional intervals which constitute the edges of each state interval box_i. Preferably, the respective weight values which are associated with the state intervals are not modified in this sub-step /5-4/. Kernel smoothing methods, also called kernel regularisation, can be applied, in particular as described in the article by Merlinge and al., Automatica 104 (2019), already mentioned. According to these methods, the random displacement which is applied to each limit of one-dimensional interval is generated according to a determined statistical distribution, which is called kernel. In particular, smoothing methods by uniform kernel, or smoothing methods by Gaussian kernel, can be applied. However, an Epanechnikov kernel smoothing method, as described now, may be preferred.

[0078] In a known manner, the Epanechnikov kernel is defined by the probability density function f(x)=3.Math.(1−x.sup.2)/4, where x is the random variable comprised between −1 and 1, the values −1 and 1 being allowed. Its expectation is zero, and its variance is equal to ⅕.

[0079] The sub-step /5-4/ therefore firstly comprises generating random corrections to be applied to each one-dimensional state coordinate interval which determines each of the state intervals box_i, then applying these corrections.

[0080] The number of state coordinates of the considered system is noted d, and the number of limits which determine each state interval of this system, that is to say each box used in the regularised particle filtering method is noted d′=2.Math.d. In the case of the aircraft 20, d=9 and d′=18.

[0081] First, N first random values are generated, denoted β.sub.i and comprised between 0 and 1, the values 0 and 1 being permitted, i being the integer numbering index of the state intervals as used previously, each according to a beta law of statistical distribution, with a first parameter equal to d and a second parameter equal to 2. In a known manner, the beta law whose parameters are d and 2 is defined by the probability density function x.sup.d/2−1.Math.(1−x).Math.Γ(d/2+2)/[Γ(d/2).Math.Γ(2)], where x denotes the random variable, and Γ denotes the gamma function. A possible way to generate the values β.sub.i in accordance with the beta law uses Cheng's algorithm which will be recalled later, with reference to [FIG. 3].

[0082] N.Math.d′ second random values are then generated, denoted v.sub.k,i, i still being the same index as previously and k being another integer index which varies from 1 to d′, each according to a normal law of statistic distribution with zero mean value and standard deviation equal to the unit, commonly called reduced normal law. The index k counts the degrees of freedom in the definition of each state interval. It identifies two one-dimensional interval limits for each state coordinate or, equivalently, a centre value and interval length for each one-dimensional state coordinate interval. In a known manner, the normal law with zero mean value and standard deviation equal to the unit is defined by the probability density function (1/π.sup.1/2).Math.exp[−x.sup.2/2], where x again designates the random variable, but in this case positive, zero or negative. It can be simulated by a sum of initial random values which are each generated according to a uniform law of statistical distribution, and such a uniform law of statistical distribution can be produced by a method of the LFSR for “linear feedback shift register” type, for example.

[0083] The first N numbers, denote ξ.sub.i, are then calculated from the random values v.sub.k,i as follows: ξ.sub.i=[Σ.sub.k=1 to d′ (v.sub.k,l).sup.2].sup.1/2. Advantageously, each of these numbers ξ.sub.i can be calculated by applying the method for calculating an estimate of X.sup.α which was described above, at X=v.sub.k,i with α=2, then at X=Σ.sub.k=1 to d′ (v.sub.k,i).sup.2 with α=½.

[0084] N second numbers, denoted α.sub.i, are then calculated according to the formula: α.sub.i=β.sub.i.sup.1/2/ξ.sub.i. Advantageously, each of these numbers α.sub.i can be calculated by again applying the method for calculating an estimate of X.sup.α which was described above, at X=β.sub.i with α=½.

[0085] Under these conditions, N.Math.d′ third numbers which are calculated in the following way: ε.sub.k,i=v.sub.k,i.Math.α.sub.i, each satisfy the Epanechnikov statistical distribution law with zero expectation and variance equal to ⅕.

[0086] Moreover, a noise amplitude, denoted h, is calculated as follows:

h=μ.Math.A.Math.N.sup.−1/(d′+4),

where A=[8.Math.c.sub.d′.sup.−1.Math.(d′+4).Math.(2.Math.π.sup.1/2).sup.d′].sup.1/(d′+4). In this expression of A, c.sub.d′ designates the volume of the hyper-sphere of dimension d′ and unit radius, and μ is an adjustment parameter which is comprised between 0 and 1. In a known manner, c.sub.d′=π.sup.d′/2/Γ(d′/2+1), where Γ still designates the gamma function. The value of c.sub.d′ can either be pre-calculated and stored so as to be available for the calculation unit 3, or be calculated by the latter, for example by using a pre-recorded table of values of the gamma function. The adjustment parameter μ controls a compromise between the efficiency of the random smoothing and that of the particle filter. Indeed, the efficiency of the particle filter results from a continuity of possible trajectories for the system, described by those of the state intervals which are maintained during several successive executions of the sequence of steps /1/ to /5/. Conversely, random smoothing produces an interference of these trajectories. The value of the parameter μ can be fixed initially for the method to be executed by the calculation unit 3. It depends in particular on the numbers N and d, as well as on the type of kernel used for the sub-step /5-4/. For example, the adjustment parameter p can be taken equal to 0.1. The use of the adjustment parameter μ is in particular described in the thesis of Merlinge, N., entitled “State estimation and trajectory planning using box particle kernels”, Paris-Saclay University, 2018.

[0087] The random modifications to be applied to the state interval box_i are then h.Math.ε.sub.k,i, where k locates the state coordinate limits, or else the central values or lengths of state coordinate intervals. These modifications can be disposed in the form of a vector E.sub.i, such that E.sub.i=[h.Math.ε.sub.k,i].sub.k=1, . . . , d′.

[0088] A possible method for applying the random modifications to the state intervals as resulting from sub-step /5-3/ will now be described.

[0089] For this purpose, each state interval box_i can be represented by a vector Ξ.sub.i of height d′, the coordinates of which group together the limits of all the one-dimensional intervals of state coordinates for this state interval, or else the central values and the lengths of its one-dimensional intervals. Then, the randomly modified values of the vector Ξ.sub.i can be obtained by replacing this vector Ξ.sub.i by Ξ.sub.i+D×E.sub.i, where D is a square matrix of dimension d′ such that the product of D by the transpose of D is equal to Σ.sub.i=1 to N Ξ.sub.i.Math.w.sub.i.Math..sup.tΞ.sub.i: D×.sup.tD=Σ.sub.i=1 to N Ξ.sub.i.Math.w.sub.i.Math..sup.tΞ.sub.i, where w.sub.i is the weight value of the state interval box_i. For example, the matrix D can be determined by the Cholesky method, which is well known to the person skilled in the art so that it is not necessary to describe it again here. Under these conditions, the vector Ξ.sub.i is to be replaced by the vector Ξ.sub.i+D×E.sub.i, to apply the random variations according to the selected kernel, of Epanechnikov in the present case, to the one-dimensional intervals of state coordinates of the state interval box_i. The set of vectors Ξ.sub.i+D×E.sub.i determine all the new state intervals that result from the modification by smoothing.

[0090] The method which has just been described for applying, in sub-step /5-4/, a smoothing by Epanechnikov kernel to the sub-intervals as they result from sub-step /5-3/, is implemented by the calculation unit 3 of the navigation system using measurement correlation 10.

[0091] The updated state intervals, as resulting from the sub-step /5-3/ or possibly from the sub-step /5-4/, with their associated weight values, constitute the result of the box regularised particle filtering method for an execution of the sequence of steps /1/ to /5/. Several iterations are chained in a recurrent way, each new iteration from the results of the previous iteration. The iteration that is performed last provides a refreshed probability distribution that characterises the true state of the aircraft 20. Furthermore, the resulting updated state intervals, associated with their respective weight values, are to be used as previous state intervals for a new execution of the sequence of steps /1/ to /5/.

[0092] A method for generating each random value β.sub.i, so as to respect the statistical distribution law Beta(d, 2) will now be described with reference to [FIG. 3]. 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.

[0093] For this purpose, the following numbers a and b are first determined in step /i/:

a which is equal to the minimum value between d and 2: a=min(d, 2), and
b which is equal to the maximum value between d and 2: b=max(d, 2).
For most applications of the invention, d is greater than 2, so that a=2 and b=d.
In step /ii/, the three numbers α, β and γ are then calculated such that:

α=a+b, which is equal to d+2,

β=[(α−2)/(2.Math.a.Math.b−α)].sup.1/2, which is equal to [d/(3.Math.d−2)].sup.1/2; and

γ=a+1/β.

The numbers α, β and γ may have been pre-calculated and stored to be directly accessible by the calculation unit 3.
In step /iii/, two numbers u.sub.1 and u.sub.2 are each randomly generated according to the uniform law of statistical distribution, between 0 and 1, for example by using an LFSR type method. Then the numbers V, W, Z, R and S are calculated in step /iv/, according to the following formulas:

V=β.Math.log[u.sub.1/(1−u.sub.1)], where log denotes the logarithmic function with base e,

W=a.Math.exp(V), where exp denotes the exponential function with base e,

Z=u.sub.1.sup.2.Math.u.sub.2,

R=γ.Math.V−log(4), and

S=a+R−W.

The value of the number Z can advantageously be calculated each time by applying the method for calculating an estimate of X.sup.α which was described above, at X=u.sub.1 with α=2. Furthermore, the values of the logarithm and exponential functions can be obtained from tables of pre-recorded values for these functions.
Steps /v/ to /viii/ are then carried out successively, forming a sequence of tests which are applied in series:
step /v/: if S+1+log(5) is greater than or equal to 5.Z, go directly to step /viii/;
step /vi/: if S is greater than or equal to log(Z), go directly to step /viii/;
step /vii/: if R+α.Math.log[α/(b+W)] is less than log(Z), return to step /iii/;
step /viii/: if a is equal to d, then the value β.sub.i which is randomly generated according to the beta law with 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 applications of the invention, where d is greater than 2, the random value β.sub.i is equal to d/(d+W).
One of the advantages of this algorithm is that the number of executions of the sequence of steps /iii/ to /viii/ which is necessary to obtain the N random values β.sub.i is predictable.

[0094] It is understood that the invention may be reproduced by modifying secondary aspects of the implementations which have been described in detail above, while retaining at least some of the advantages mentioned. Among such possible modifications, the following are mentioned without limitation: [0095] alternative algorithms may be used for some steps or sub-steps; [0096] the inertial system 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 speed values, according to the three spatial coordinates, in addition to the three spatial coordinate values which identify a position for the vehicle, of the three speed values and the three attitude angular values. In this case, each state of the system comprises values for fifteen state coordinates; [0097] in step /1/, a dynamic state change model may be used for the system, this model being able to take account, when it is a vehicle, of commands which are applied to one or multiple engine(s) and to a vehicle attitude control system; and [0098] the calculation unit may be a processor with an ARM structure, a processor with several calculation cores, one or more graphics processor(s), etc., instead of an FPGA, DSP, CPU or RISC type chip.

[0099] Finally, the invention may be applied to fields other than aeronautics. For example, a box regularised particle filtering method which is in accordance with the invention may also be implemented for a ground vehicle, a maritime surface ship, a submarine, a satellite or a space probe, using each time a baseline and true state measures that are appropriate.