Patent Application
20230026000
The present invention relates to a particle filtering method, as well as a calculation unit and a navigation system using measurement correlation which implement such a method. More particularly, the particle filtering method is of the regularised box type.
The navigation function of an aircraft includes the estimation of its instantaneous position, its instantaneous speed and its instantaneous orientation, called attitude, in the navigation reference frame, also called local geographic trihedron. The set of instantaneous values of the coordinates of position, that is to say the latitude, longitude and altitude of the aircraft, of speed, comprising a speed coordinate in the northern direction, a speed coordinate in the eastern direction and descent speed coordinate, and attitude angles of the aircraft, comprising a roll angle, pitch angle and yaw angle, constitute the instantaneous state of the system formed by the aircraft. The aircraft can be an aeroplane, a drone or any self-propelled airborne vehicle, without limitation.
The acceleration and the angular speed of the aircraft are measured repeatedly and at high speed each along three axes, for example at a repetition frequency of approximately 1000 Hz (Hertz), using accelerometers and gyrometres of an inertial system which is on board the aircraft. The navigation system then delivers estimates of the position and speed coordinates, and the attitude angles of the aircraft, by integrating the results of the accelerometric and gyrometric measurements over time. However, each measurement of acceleration and angular speed is affected by an error, which is mainly composed of bias, a scale factor and a random process, and the accumulation of measurement errors results into an error which affects the estimation of the instantaneous state of the aircraft. This error on the instantaneous state which is estimated grows as a function of time, and is commonly called inertial drift. It is then necessary to associate the inertial system with at least one additional sensor, in order to correct or reduce the inertial drift.
Usual methods for correcting or reducing the inertial drift consist in using geolocation signals, such as satellite navigation signals, for example of the GNSS type, for “global navigation satellite system”, or in using signals which are produced by beacons located on the ground, for example radionavigation or GBAS for “ground-based augmentation system” type signals or else to receive a location of the aircraft by radio which was carried out using radar. But such methods are sensitive to jamming or decoys, to the availability of coverage of the area where the aircraft is located in location signals or in communication, etc. It is then desirable in certain circumstances to have on board the aircraft a method which is autonomous for correcting and/or reducing inertial drift. For this purpose, it is common to associate the inertial system with a telemetry probe, which measures a distance between the aircraft and the ground. This telemetry probe can be a radio altimeter, a laser telemetry, etc. without limitation. It measures the distance between the aircraft and the ground in a direction which may or may not be fixed in relation to the aircraft. When this distance measurement direction can be variable, its orientation relative to the aircraft is known. The inertial system is then also associated with a calculation unit which correlates the results of successive measurements which are carried out by the telemetry probe with the instantaneous state of the aircraft as estimated by the inertial system. More specifically, a characterisation of the area overflown by the aircraft is stored on board the aircraft, for example in the form of a relief map which associates a relief height value with each pair of values of latitude and of longitude. Such a relief map record may be in the form of a table, with the latitude and the longitude constituting the table entries, and the relief height values constituting the reading responses in the table. Alternatively, the characterisation of the overflown area can be stored in the form of an analytical function which allows to calculate the relief height values according to the values of latitude or longitude, or in any other suitable form. Then, with each new estimate of the instantaneous position of the aircraft which is produced by the navigation system, a value of the distance which should exist between the aircraft and the ground is obtained by interrogating the characterisation of the overflown area as stored on board the aircraft, in accordance with the estimated position of the aircraft. Optionally, the aircraft-ground distance value may result from a calculation which combines the estimated position of the aircraft with the result(s) of one or more interrogation(s) of the characterisation of the overflown area which is stored on board the aircraft, in particular when the measurement direction of the telemetry probe is oblique with respect to the altitude axis. Such a calculation is known to the person skilled in the art, so that it is not useful to repeat it here. The aircraft-ground distance value which is thus estimated is then compared with the measurement result which is delivered by the telemetry probe. Such a measurement correlation navigation method is commonly referred to as terrain correlation navigation. Several variants have been proposed, but some of them show great sensitivity to the existence of non-linearities in the terrain profiles. In other words, they have the disadvantage of a lack of robustness in their efficiency of convergence towards the true state of the aircraft, depending on the possible terrain profiles.
The methods using terrain correlation navigation which are based on a Box Regularised Particle Filter, or BRPF, allow to correct the inertial drift with greater robustness, compatible with the existence of non-linearities and field ambiguities. The article by Merlinge, N., Dahia, K., Piet-Lahanier, H., Brusey, J., & Horri, N. entitled “A Box regularised particle filter for state estimation with severely ambiguous and non-linear measurements”, Automatica (2019), Vol. 104, pp. 102-110, describes such a method. Each of these methods still consists in iteratively calculating an instantaneous state of the aircraft from a last state previously determined, and in correlating the distance measurement result which is obtained by the telemetry probe with a distance value which is moreover reconstituted from the characterisation of the terrestrial relief on board the aircraft and the position and attitude values. But, a box regularised particle filter proceeds by manipulating state intervals, of dimension nine when each state of the aircraft consists of three position coordinates, three speed coordinates and three attitude angles as described above. In addition, a weight value with probabilistic significance is associated with each state interval: the weight of each state interval corresponds to the probability that the true state of the aircraft lies in this state interval. However, these methods using terrain correlation navigation that are based on particle filters have not yet been implemented for real aircraft, because of the very large calculation resources that are required. Indeed, it is required for many aeronautical applications that the method using terrain correlation navigation which is used can be implemented by a calculation unit of the Field-Programmable Gate Array circuit type, or FGPA. However, circuits of this type have capacities which are still too limited.
In a known manner, a box regularised particle filtering method comprises 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, each state interval being possibly multidimensional;
/2/ a step of measuring a true state;
/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 this subsequent state interval as resulting from step /3/, with a size of the same subsequent state interval as resulting from step /1/ before step /3/, and a weight of the previous state interval from which this 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 ni sub-intervals originating from a division of the subsequent state interval i, i being an integer numbering index of the subsequent state intervals, and ni being a number of the sub-intervals which replace this 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/.
These steps comprise various calculations, some of which require significant calculation resources. This is the case in particular for those of the previous steps which require calculating a result of a function of the type Xα, where X is a variable whose value changes between different executions of the step, and α is a non-zero exponent whose value is constant for all executions of the same step. For this reason, the box regularised particle filtering method is not compatible with an implementation on board an autonomous vehicle. In particular, it cannot be implemented by a calculation unit of the FGPA type.
Document WO 2018/002932 describes a method for determining the position of a vehicle relative to a traffic lane, which is based on captured images of light sources, and which uses a particle filtering algorithm.
Document FR 2 832 796 describes a system using terrain correlation navigation for aircraft.
From this situation, an object of the present invention is to optimise the use of a box regularised particle filter, in particular so that this filter can be implemented on board an aircraft. Beyond that, the object of the invention is to contribute to the fact that a box regularised particle filtering method can be implemented by an FPGA-type circuit.
To achieve at least one of these objects or another, a first aspect of the invention proposes a box particle filtering method for predicting a state of a system, which comprises the steps /1/ to /5/ mentioned above, and wherein at least one of these steps /1/ to /5/ comprises at least one calculation of an estimate of the value of Xα, where X is a variable number for each execution of the step and α is a non-zero exponent value that is constant between different executions of the step, the number X being positive if the exponent α is not an integer and non-zero if the exponent α is negative. The concerned system can be a land, air, sea or space vehicle which is equipped with a navigation system using measurement correlation. The calculation of the estimate of the value of Xα comprises the following steps:
/a/ writing the number X in the form X=±(1+m)·2ex, 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·(m+ex+B), where L=2n 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α in the form:
|(Xα)=α·I(X)+L·(1−α)·(B−σ), where σ is a constant number whose value is recorded; and
/c/ obtaining the estimate of the value of Xα from the binary representation I(Xα).
Such a calculation, which replaces the estimation of the function of Xα by a calculation based on the binary representation of the number X, is particularly economical in calculation resources, and short in calculation time. Furthermore, it can be performed by an FPGA-type circuit.
In particular embodiments of the invention, the step of the box regularised particle filtering method which comprises the calculation of the estimate of the value of Xα may be the step /5/, in particular a sub-step of this step /5/ which consists in calculating a measure value of a part of a subsequent state interval of the system. This subsequent state interval part of the system may be a system position interval, a system speed interval, or a system attitude interval. Alternatively or in addition, the calculation of the estimate of the value of Xα can be part of 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/. The calculation that is the object of the invention can therefore be used on the one hand for the measure value of a part of the subsequent state interval of the system, and on the other hand to apply the probability kernel smoothing.
Generally for the invention, and in particular to obtain a more exact estimate of the value of Xα, the calculation of the estimate of the value of Xα can further comprise executing at least once the following additional step:
/d/ calculating a new estimate of the value of Xα from a previous estimate of the value of Xα, by applying a recursive algorithm for approximate equation solving to the equation Y1/α−X=0 of unknown Y, the estimate of the value of Xα which was obtained in step /c/ being used as a previous estimate for a first application of the algorithm, and the new estimate of the value of Xα which is produced by a qth application of the algorithm forming the previous estimate of the value of Xα for the (q+1)th application of the algorithm, if such a (q+1)th application of the algorithm is performed, q being an integer number greater than or equal to 1.
For such a refinement of the invention, the recursive algorithm for approximate equation resolution which is used in step /d/ may be Newton's method, or another equivalent algorithm or which performs the same function.
In preferred embodiments of the invention, at least one of the following additional features can be reproduced optionally, alone or in combination of several of them:
A second aspect of the invention proposes a calculation unit, this calculation unit comprising at least one first input which is adapted to receive results of repeated measurements of acceleration and angular speed of a system, and a second input which is adapted to receive results of repeated measurements of a true state of the system, additional to the acceleration and angular speed measurements. The calculation unit is then arranged to execute a box regularised particle filtering method which is in accordance with the first aspect of the invention. In this way, the calculation unit outputs 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 this state interval.
Such a calculation unit can be of one of the following types:
A third aspect of the invention relates to a navigation system using measurement correlation, which is adapted to be on board a vehicle and which comprises:
Finally, a fourth aspect of the invention proposes a vehicle which comprises a navigation system using measurement correlation in accordance with the third aspect of the invention. Such a vehicle may be an aircraft, in particular an aeroplane, a flying drone or any self-propelled airborne vehicle, or a vehicle capable of moving on the ground, in particular a mobile drone on the ground, or a ship, a submarine, a space machine, in particular a space probe, satellite, etc. without limitation. Depending on the case, the measurement system may be a telemetry probe which is intended to measure a distance between the vehicle and a point on the ground, a system for tracking landmarks or stars, a sonar for measuring a height of water under the ship or submarine, etc.
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:
For reasons of clarity, the dimensions of the elements which are represented symbolically in [
The mode of calculating an estimate of the value of Xα 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.
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)·2ex,
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.
Then, when X is positive, we have:
log2X=ex+m+σ,
where σ is a fixed real number allowing to minimise an error on the value of log2X, 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=log2X−σ.
Moreover, the number X can be represented in binary form by I(X) defined by:
I(X)=L·(m+ex+B),
where L=2n, 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 log2X, it becomes:
I(X)=L·[log2X−σ+B], that is to say: log2X=I(X)/L+σ−B.
However, in the same way as I(X) in the previous line, the binary representation of Xα is:
I(Xα)=L·[log2(Xα)−σ+B]
But log2(Xα)=α·log2(X), therefore: I(Xα)=L·[α·log2(X)−σ+B], and replacing log2(X) by its expression according to the binary representation I(X), it becomes:
I(Xα)=L·[α·(I(X)/L+σ−B)−σ+B], that is to say: I(Xα)=α·I(X)+L·(1−α)·(B−σ).
An approximate value of Xα, denoted Y0, can then be reconstructed from the binary representation of Xα 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α 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α 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 Y0 is satisfactorily suitable for replacing Xα, given the simplicity of the method for obtaining this approximate value Y0, as just described.
For the particular case where the exponent α is equal to 2:
I(X2)=2·I(X)−L·(B−σ).
For the particular case where the exponent α is equal to −2:
I(X−2)=−2·I(X)+3·L·(B−σ).
For the particular case where the exponent α is equal to ½:
I(X1/2)=I(X)/2+0.5·L·(B−σ).
For the particular case where the exponent α is equal to −½:
I(X−1/2)=−I(X)/2+1.5·L·(B−σ).
For applications where the approximate value Y0 does not constitute a sufficiently precise evaluation of Xα, 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)=Y1/α−X and to the equation f(Y)=0. Successive approximate values Yq, q being an integer index for numbering these values, can thus be obtained according to the formula: Yq+1=Yq−f(Yq)/f′(Yq), where f(Yq) is the value of the function derived from f, estimated for the value Yq. That is to say, by calculating the expression of f′(Y) from that of f(Y):
Y
q+1=(1−α)·Yq+α·X·Yq(α−1)α, for q=0,1,2, . . . .
For the particular case where the exponent α is equal to 2, it becomes:
Y
q+1
=−Y
q+2·X·Yq1/2.
In particular, the first-order approximate value of X2 is:
Y
1
=−Y
0+2·X·Y01/2.
The value of Yq1/2 can be estimated each time by using the formula I(Xα)=α·I(X)+L·(1−α)·(B−σ), and by replacing in this formula X by Yq and a by ½.
For the particular case where the exponent α is equal to −2, it becomes:
Y
q+1=3·Yq−2·X·Yq3/2.
In particular, the first-order approximate value of X−2 is:
Y
1=3·Y0−2·X·Y03/2.
The value of Yq3/2 can be estimated each time by again using the formula I(Xα)=α·I(X)+L·(1−α)·(B−σ), and replacing in this formula X by Yq and a by 3/2.
For the particular case where the exponent α is equal to ½, it becomes:
Y
q+1=0.5Yq+0.5X/Yq.
In particular, the first-order approximate value of X1/2 is:
Y
1=0.5·Y0+0.5·X/Y0.
For the particular case where the exponent α is equal to −½, it becomes:
Y
q+i=(3/2)Yq−0.5·X·Yq3.
In particular, the first-order approximate value of X−1/2 is:
Y
1=(3/2)·Y0−0.5·X·Y03.
The approximate value Y0, as obtained by using the binary representations of numbers, and the values Yq, 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.
In accordance with [
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.
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.
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.
Each sequence of steps comprises a prediction step, denoted /1/ in [
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.
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.
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.
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.
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 [
Step /5/ comprises first determining a division number ni 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 [
Each state interval box_i is then intended to be divided into ni sub-intervals, in sub-step /5-3/, for example in one of the ways that are described in the following articles:
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.
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 x1, x2 and x3, are state coordinates of a first type, the three speeds, denoted v1, v2 and v3, are state coordinates of a second type, and the three attitude angles, denoted θ1, θ2 and θ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:
Δxj_n=Δxj·(Δx12·Δx22+Δx32)−1/2,
Δvj_n=Δvj·(Δv12+Δv22+Δv32)−1/2, and
Δθj_n=Δθj·(Δθ12+Δθ22+Δθ32)−1/2,
for j=1, 2 and 3 in each case, where Δxj, Δvj and Δθ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 Δxj2, Δvj2 and Δθj2 can be calculated using the method that was presented above for Xα, with α equal to 2. Then each of the factors (Δx12+Δx22+Δx32)−1/2, (Δv12+Δv22+Δv32)−1/2 and (Δθ12+Δθ22+Δθ32)−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 Δxj_n, Δvj_n and Δθ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.
These normalised values Δxj_n, Δvj_n and Δθ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 Δxj_n, Δvj_n and Δθj_n, considering their absolute values. For example, the state interval box_i is the largest according to the position coordinate x1 if Δx1_n is the largest of the normalised values Δxj_n, Δvj_n and Δθj_n, or it is the largest according to the coordinate speed v2 if Δv2_n is the greatest of the nine normalised values, etc.
In sub-step /5-3/, the state interval box_i is divided into n contiguous sub-intervals, by dividing into ni 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 Δxj_n, Δvj_n and Δθ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 ni 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 ni. These new state intervals are then re-numbered by the index i, advantageously from 1 to N, for the continuation of the filtering method.
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.
In a known manner, the Epanechnikov kernel is defined by the probability density function f(x)=3·(1−x2)/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 ⅕.
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.
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·d. In the case of the aircraft 20, d=9 and d′=18.
First, N first random values are generated, denoted β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 xd/2−1·(1−x)·Γ(d/2+2)/[Γ(d/2)·Γ(2)], where x denotes the random variable, and Γ denotes the gamma function. A possible way to generate the values βi in accordance with the beta law uses Cheng's algorithm which will be recalled later, with reference to [
N·d′ second random values are then generated, denoted vk,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/π1/2)·exp[−x2/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.
The first N numbers, denote ξi, are then calculated from the random values vk,i as follows: ξi=[Σk=1 to d′ (vk,l)2]1/2. Advantageously, each of these numbers ξi can be calculated by applying the method for calculating an estimate of Xα which was described above, at X=vk,i with α=2, then at X=Σk=1 to d′ (vk,i)2 with α=½.
N second numbers, denoted αi, are then calculated according to the formula: αi=βi1/2/ξi. Advantageously, each of these numbers αi can be calculated by again applying the method for calculating an estimate of Xα which was described above, at X=βi with α=½.
Under these conditions, N·d′ third numbers which are calculated in the following way: εk,i=vk,i·αi, each satisfy the Epanechnikov statistical distribution law with zero expectation and variance equal to ⅕.
Moreover, a noise amplitude, denoted h, is calculated as follows:
h=μ·A·N
−1/(d′+4),
where A=[8·cd′−1·(d′+4)·(2·π1/2)d′]1/(d′+4). In this expression of A, cd′ 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, cd′=πd′/2/Γ(d′/2+1), where Γ still designates the gamma function. The value of cd′ 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.
The random modifications to be applied to the state interval box_i are then h·ε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 Ei, such that Ei=[h·εk,i]k=1, . . . , d′.
A possible method for applying the random modifications to the state intervals as resulting from sub-step /5-3/ will now be described.
For this purpose, each state interval box_i can be represented by a vector Ξ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 Ξi can be obtained by replacing this vector Ξi by Ξi+D×Ei, where D is a square matrix of dimension d′ such that the product of D by the transpose of D is equal to Σi=1 to N Ξi·wi·tΞi: D×tD=Σi=1 to N Ξi·wi·tΞi, where wi 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 Ξi is to be replaced by the vector Ξi+D×Ei, 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 Ξi+D×Ei determine all the new state intervals that result from the modification by smoothing.
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.
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/.
A method for generating each random value βi, so as to respect the statistical distribution law Beta(d, 2) will now be described with reference to [
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·a·b−α)]1/2, which is equal to [d/(3·d−2)]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 u1 and u2 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=β·log[u1/(1−u1)], where log denotes the logarithmic function with base e,
W=a·exp(V), where exp denotes the exponential function with base e,
Z=u
1
2
·u
2,
R=γ·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α which was described above, at X=u1 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+α·log[α/(b+W)] is less than log(Z), return to step /iii/;
step /viii/: if a is equal to d, then the value β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 β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 βi is predictable.
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:
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.
| Number | Date | Country | Kind |
|---|---|---|---|
| 1914339 | Dec 2019 | FR | national |
This application is the U.S. national phase of International Application No. PCT/EP2020/084817 filed Dec. 7, 2020 which designated the U.S. and claims priority to FR 1914339 filed Dec. 13, 2019, the entire contents of each of which are hereby incorporated by reference.
| Filing Document | Filing Date | Country | Kind |
|---|---|---|---|
| PCT/EP2020/084817 | 12/7/2020 | WO |
