Satellite signal receiver with speed computing integrity control

A satellite positioning receiver with velocity calculation integrity check. In order to increase the reliability of the velocity measurements by the receiver, it is proposed to measure the divergence between the measurements taken on n satellite director axes (n>4) and to given an alarm if this divergence exceeds a predetermined threshold. The divergence is measured by the mean square error between the calculated velocity vector and the velocities as measured along the director axes and projected on to the calculated velocity vector.

BACKGROUND OF THE INVENTION 
 1. Field of the Invention 
 The invention relates to satellite positioning receivers such as GPS 
 (Global Positioning System) receivers. 
 2. Discussion of the Background 
 The GPS system uses a constellation of satellites which move around the 
 earth on very precisely determined orbits, that is to say it is possible 
 to know the position of an arbitrary satellite at any time. The satellites
 transmit radiofrequency signals, containing navigation data and codes 
 which make it possible to identify each satellite. These codes phase 
 modulate a carrier frequency. A GPS receiver, on the ground or on a land, 
 air or sea vehicle, can receive the signals from several satellites 
 simultaneously, precisely calculate its distance from each of the 
 satellites, and deduce therefrom its precise position in latitude, 
 longitude and altitude in a terrestrial frame. It can also deduce 
 therefrom the precise date and time of the reception in the time frame of 
 the GPS system. It can lastly deduce therefrom, by Doppler measurements, 
 its own velocity vector in the terrestrial frame (the case of a receiver 
 mounted on a moving vehicle). 
 In the GPS system, each satellite is identified by a pseudo-random code 
 which is individual to it and repetitively (for example every millisecond)
 modulates a carrier frequency which the satellite transmits. There are 
 systems similar to GPS, in particular the GLONASS system, in which this 
 pseudo-random code also exists even though it is not used to identify 
 individual satellites. The invention which will be described is directly 
 applicable to the GLONASS system, but for the sake of simplicity reference
 will be made below only to the GPS system, and more precisely the "civil" 
 part of the GPS system which also has a military part. 
 The pseudo-random code is a long code (1023 bits at 1.023 MHz, i.e. 1 
 millisecond), and one of its main functions is to make it possible to 
 extract the satellite's signal from a noise level much higher (for example
 30 dB) than the level of the signal. This technique is now widely known as
 spread spectrum transmission. The signal is extracted from the noise using
 an operation, in the receiver, of correlation between the received signal 
 and a periodic pseudo-random code which is identical to the one expected 
 to be found in the signal. If the codes do not coincide temporally, there 
 is no correlation between the received signals and the local code 
 generated by a local code generator; if they almost coincide, there is 
 some degree of correlation, the correlation energy becoming stronger as 
 the coincidence becomes more exact. It is therefore possible to establish 
 a correlation signal making it possible to slave a local code generator 
 until exact coincidence is obtained between the local code and the code 
 modulating the signal which the satellite transmits. A code slaving loop 
 then makes it possible to maintain this coincidence. 
 The pseudo-random code is transmitted by the satellite at extremely precise
 times which are known at the receiver. Use is made of the correlation 
 operation to determine the arrival time of this code in the receiver: the 
 characteristic time or epoch of transmission of the local code is 
 determined, and since this local code coincides with the received code 
 when the maximum correlation is established, this time represents the 
 arrival time of the received code. The difference between a time at which 
 the code is transmitted via the satellite and a time at which the code is 
 received by the receiver makes it possible to determine a propagation time
 of the signals between the satellite and the receiver. Knowing that the 
 propagation velocity of the signals is the velocity of light, the distance
 between the receiver and a given satellite can be calculated. The same 
 operation performed on two other satellites makes it possible, by 
 triangulation, to determine the exact position of the receiver relative to
 the three satellites. 
 By using a fourth satellite, the clock errors of the receiver are 
 eliminated, the clock of the receiver not being as precise as that of the 
 satellites. Further to the position of the receiver, it is then possible 
 to calculate a precise time for the position measurement, in the time 
 frame of the GPS satellites. 
 The position of each of the satellites is known at any time: it is 
 calculated from tables which are stored in the receiver and are updated by
 the navigation message broadcast by the satellites. The velocity of the 
 satellites at any time can also be calculated on the basis of these 
 tables. 
 It is possible to determine, on the basis of the signals sent by four 
 satellites, the time and the position of the receiver relative to the four
 satellites. Furthermore, by changing co-ordinates, the position of the 
 receiver in a fixed terrestrial frame is obtained. 
 Similarly, the velocity of the receiver is calculated on the basis of a 
 Doppler-effect measurement on the carrier frequency of the radiofrequency 
 signal sent by the satellites. It is therefore possible to calculate the 
 relative velocity of the receiver with respect to each of the satellites, 
 along the director axis which joins this satellite to the receiver. Four 
 satellites are needed to eliminate the time ambiguity. Four different 
 relative velocity vectors are obtained, along the director axes joining 
 the receiver to the four satellites. Simple calculations make it possible 
 to determine the temporal drift of the clock of the receiver relative to 
 the GPS time, and the velocity of the receiver along three axes in the 
 terrestrial frame on the basis of these four velocity vectors and the 
 following information: 
 the directions of the receiver-satellite director axes with respect to a 
 fixed terrestrial frame (longitude, latitude, altitude); these directions 
 are themselves obtained by knowledge of the position of the receiver at a 
 given time and the position of each satellite at the same time; 
 the individual velocities of the satellites in the terrestrial frame at 
 this time. 
 However, if more than four satellites are used, redundant information is 
 obtained. This is the case for the professional-quality receivers used, in
 particular, in aeronautics. 
 If using all the redundant information led exactly to the same velocity 
 calculation results, it would be possible to make do with taking any four 
 satellites from the satellites in the constellation observed at a given 
 time. 
 However, the measurements are affected by various imprecisions, so that the
 redundancy is not perfect. Furthermore, a satellite may be operating 
 defectively at a given time and therefore give aberrant information 
 interfering with the velocity determination. 
 In certain applications of GPS receivers, it may be important to determine 
 the velocity with precision and certainty. This is the case, for example, 
 for receivers used to assist in the landing of aircraft. 
 SUMMARY OF THE INVENTION 
 This is why the present invention proposes to calculate the relative 
 velocity between the receiver and the satellites using more than four 
 satellites, to deduce information about the integrity of the calculated 
 velocity from the various measurements, and to give an indication of the 
 fact that the measurement is not valid if the integrity of the velocity 
 measurement is not satisfactory. 
 The term "integrity" of the measurement is intended to mean a quantity 
 which represents the variable compatibility of the measurements taken on a
 set of more than four satellites, that is to say which represents the 
 extent to which the velocity measurements obtained by taking one group of 
 four satellites among n are identical with the velocity measurements 
 obtained by choosing other groups of four satellites among the n. 
 More precisely, the invention proposes a satellite positioning receiver, 
 comprising position calculation means which simultaneously use a number of
 n satellites at least equal to 4, receiver velocity calculation means 
 which can determine the velocity of the receiver along the director axes 
 joining the receiver to each of the n satellites, and means for 
 calculating a velocity vector of the receiver in a fixed terrestrial frame
 on the basis of the velocities along the director axes and on the basis of
 a matrix of vectors representing the directions of these director axes in 
 the fixed frame, this receiver being characterized in that it furthermore 
 comprises means for calculating a value representative of the divergence 
 between the velocity measurements and the n director axes, this value 
 representing information about the validity of the velocity vector of the 
 receiver in the fixed frame. 
 If the integrity is not sufficient, that is to say if the measurement of 
 the divergence exceeds a determined threshold, the receiver may give an 
 indication that the measurement is invalid, but may also initiate a search
 for the satellite which has caused an abnormal divergence, and temporary 
 elimination of this satellite. 
 Preferably, the quantity which represents the divergence, or lack of 
 integrity, is a residue .DELTA.V equal to the norm of a vector which is 
 the product S.multidot.VD, where S is an n.times.n matrix defined below, n
 being the number of satellites used and VD being a measured velocity 
 vector whose components are the n velocity measurements along the director
 axes between the receiver and each of the n satellites; the matrix S is 
 calculated on the basis of the following formula, where H is the director 
 cosine matrix, that is to say an n.times.4 matrix representing the 
 directions of the director axes in a fixed terrestrial spatio-temporal 
 frame; H.sup.T is the transpose of this matrix, and I is the n.times.n 
 identity matrix: 
EQU S=I-H(H.sup.T.multidot.H).sup.-1.multidot.H.sup.T. 
 The norm of the matrix product S.multidot.VD is 
 .DELTA.V=(.vertline.S.multidot.VD.vertline.).sup.2. 
 The integrity threshold is preferably selected on the basis of a false 
 alarm probability value acceptable in the application in question. This 
 threshold depends on the number of satellites used. 
 The quantity which represents the divergence and which is compared with a 
 threshold is preferably a quantity normalized with respect to an estimated
 measurement noise of energy .sigma..sup.2. In this case, the calculated 
 residue is the normalized residue .DELTA.V.sub.nr =.DELTA.V/.sigma..sup.2.
 The estimated noise may be a fixed value selected a priori as a function of
 the application, or a value measured on the basis of all the velocity 
 measurements taken previously. In one advantageous solution, a moving 
 average of the ratio .DELTA.V/n-4 between the residue and the number n-4 
 of degrees of freedom in the velocity calculation is used as an estimate 
 of the noise .sigma..sup.2. The number of degrees of freedom represents 
 the number of satellites used less the number of satellites strictly 
 necessary for this calculation. 
 Means for calculating a velocity error limit, representing the quality of 
 the integrity check made, are moreover provided. This limit is a maximum 
 velocity error which the system has a probability PND of not detecting. In
 other words, it is a velocity error linked with a probability threshold 
 selected such that, when a certain probability referred to as probability 
 non-detection is picked, there is this probability of the system not 
 seeing this velocity error. 
 This limit may be displayed for the user, or else if it exceeds a 
 determined threshold it may be used to give an indication of lack of 
 reliability in the velocity measurement.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 
 The microprocessor 60 has two functions: 
 on the one hand, it works out digital data used by the digital signal 
 processing circuit 50, on the basis of digital data delivered by this 
 circuit; in particular, it performs numerical calculations necessary for 
 the digital slaving loops present in the digital processing circuit 50; 
 and on the other hand it gives final results of position, time and velocity
 calculation to the user, that is to say either on a digital display screen
 or on a digital bus to other equipment which need the results. 
 It could clearly be possible to have two separate processes for performing 
 these two functions. In the embodiment with a single microprocessor 60, a 
 bus 70 has been represented for exchanges between the microprocessor 60, 
 the processing circuit 50, an input/output peripheral 80, the working 
 memories 90, and the program memories 100 which contain the programs 
 needed for the microprocessor to function. 
 Very briefly, the digital signal processing circuit has either a single 
 processing channel, with the information from the various satellites being
 multiplex-processed, or preferably several channels which each work in 
 parallel on a determined satellite. 
 Each channel has a double slaving loop: carrier phase slaving and code 
 phase slaving. 
 The carrier phase loop essentially uses a local oscillator with digital 
 phase control, delivering a periodic (sawtooth) digital phase at a 
 frequency corresponding to the transposed carrier frequency, taking into 
 account the Doppler effect to which the carrier frequency broadcasted by a
 satellite is subjected. The Doppler effect is taken into account by the 
 very fact of the existence of the slaving loops. The microprocessor 60 
 calculates a carrier phase error signal; this signal is used to control 
 the local oscillator in order to slave a local carrier frequency to the 
 carrier frequency received from the satellite. 
 The code phase slaving loop has a local code generator driven by an 
 oscillator with digital phase control. It makes it possible to slave the 
 local codes to the code received from the satellite and then to be able to
 determine the exact temporal position of the local codes thus slaved. The 
 local code is correlated with the code received from the satellite; the 
 correlation signal is calculated by the microprocessor and is used to 
 slave the loop in order to bring the local code into synchrony with the 
 code received from the satellite. 
 The two slaving loops, for code and carrier, take into account the Doppler 
 frequency shift on the carrier frequency and on the code, which result 
 from the relative motion of the aircraft and the detected satellite. This 
 Doppler shift can be measured in the loops. 
 The GPS time and position calculations are performed on the basis of the 
 status of the slaving loops at a determined measurement time. At this 
 time, the exact status of the phase of the two oscillators with digital 
 phase control are read. 
 The slaving loops provided in the receiver act to lock a local frequency 
 onto the carrier frequency received from the satellites. The shift between
 this local frequency and the stable and known frequency transmitted by the
 satellites gives an indication of Doppler shift and therefore the 
 difference between the velocity of the satellite and the velocity of the 
 receiver along the axis joining the satellite to the receiver. As will be 
 seen below, this Doppler indication needs to be corrected for the local 
 clock frequency error of the receiver, which error can be measured by the 
 temporal drift of the local clock relative to the GPS time. 
 The receiver therefore calculates, for the n different satellites observed 
 at a given time (n greater than or equal to 4), its relative position with
 respect to these satellites. It deduces therefrom the directions of the 
 observation axes of each of the satellites. 
 It moreover measures the pseudo-velocity, or relative velocity PV.sub.i 
 between the receiver and the satellite of rank i, along the axis joining 
 the receiver to the satellite of rank i (i varying from 1 to n). The 
 absolute velocity of the satellite of rank i in the terrestrial frame is 
 known. The receiver then calculates the satellite velocity projected onto 
 the director axis of rank i: this is the velocity Vrsat.sub.i. A clock 
 drift correction c.Dersat.sub.i for the satellite can be calculated on the
 basis of the navigation data sent by the satellite; c is the velocity of 
 light; Dersat.sub.i is the clock drift for the satellite of rank i, that 
 is to say the variation, in the course of time, of the time bias of the 
 clock of the satellite: the time bias of the clock of the satellite is the
 discrepancy between the time which this clock gives and the time defined 
 on the ground for the entire GPS system. This bias can vary in the course 
 of time because of frequency discrepancies of the clocks of the satellites
 relative to a theoretical value. 
 The receiver's absolute velocity, in the fixed terrestrial frame but 
 projected onto the director axis of rank i, is deduced therefrom. Let 
 Vdop.sub.i be this absolute velocity. 
EQU Vdop.sub.i =PV.sub.i +c.Dersat.sub.i -Vrsat.sub.i 
 Taken together, all the velocities Vdop.sub.i for i=1 to n form a measured 
 absolute velocity vector VD which is an n.times.1 matrix of n individual 
 velocities. 
 In a fixed three-dimensional terrestrial frame, the absolute velocity of 
 the receiver could be expressed by an absolute velocity vector with three 
 components V1, V2, V3. However, given the local clock drift of the 
 receiver, which has a direct influence on the velocity measurement since 
 the velocities are determined by the Doppler effect on the frequencies, it
 is preferable to express the absolute velocity of the receiver in a 
 four-dimensional frame: three space dimensions in a fixed terrestrial 
 frame and one time dimension in the form of a calculable velocity drift 
 V4=c.multidot.Dhr, where Dhr is a variation in the clock bias of the 
 receiver in the course of time. The clock bias of the receiver, which the 
 GPS receiver itself can calculate so long as there are at least four 
 observed satellites, is the discrepancy between the time given by the 
 clock of the receiver and the time defined by the GPS system in its 
 entirety. Dhr is the time variation of this bias. The velocity correction 
 is the product of Dhr times the velocity of light c. 
 The absolute velocity of the receiver in the terrestrial frame can 
 therefore be expressed in the form of a vector VABS with four components 
 V1, V2, V3, V4, respectively representing the longitude, latitude and 
 altitude velocities, and the correction due to the clock drift of the 
 receiver. 
 The absolute velocities V1, V2, V3, V4 can be deduced from the velocity 
 vector VD with n components. 
 If H denotes the (n.times.4) matrix of the director cosines, that is to say
 a matrix of n vectors (with four dimensions in the fixed terrestrial frame
 with one unitary time component) representing the directions of the 
 director axes of rank i=1 to n, then the following equation should be 
 satisfied: 
EQU VD=H.multidot.VABS (1) 
 which means that the velocity measurements for the receiver along the 
 director axes (which are considered to form a fixed frame) can be 
 projected into the four-dimensional terrestrial frame to obtain a 
 four-dimensional velocity vector uniquely. 
 The director cosine matrix is a matrix of n rows of four coefficients 
 C.sub.i,x, C.sub.i,y, C.sub.i,z, 1, where C.sub.i,x, C.sub.i,y, C.sub.i,z 
 represent the cosines of the angles between the i.sup.th axis and the axes
 Ox, Oy, Oz (longitude, latitude, altitude) of the terrestrial frame, O 
 being the position of the receiver. 
 In reality, because of the measurement noise, there is a discrepancy 
 between VD and H.multidot.VABS, and the integrity of the velocity 
 measurement using more than four satellites can be represented by a 
 measurement of this discrepancy between VD and H.multidot.VABS. 
 If an optimization criterion is picked to find the vector VABS which is 
 most representative of the velocity in the terrestrial frame on the basis 
 of the vector VD, then VABS can be calculated. Mathematically, if a "least
 squares" criterion is adopted, the above matrix equation 
 VD=H.multidot.VABS admits the following solution: 
EQU VABS=(H.sup.T.multidot.H).sup.-1.multidot.H.sup.T.multidot.VD (2) 
 H.sup.T is the transpose of H; it is a 4.times.n matrix. 
 The integrity of the calculation of the velocity can be evaluated 
 quantitatively on the basis of the following value .DELTA.V, which will be
 referred to as the "calculation error residue": 
EQU .DELTA.V=(.vertline.VD-H.multidot.VABS.vertline.).sup.2 (3) 
 which is to say the residue .DELTA.V is the norm of the difference vector 
 between the measured velocity vector VD and the vector product 
 H.multidot.VABS. It is moreover this quantity which is minimized when 
 using the least squares criterion. 
 This residue at .DELTA.V increases as the measurements taken by the n 
 satellites coincide less. It is a measure of the divergence between the n 
 measurements along the director axis. 
 If S denotes the n.times.n matrix equal to 
 I-H(H.sup.T.multidot.H).sup.-1.multidot.H.sup.T, where I is the n.times.n 
 matrix, then the following may be written on the basis of (2) and (3): 
EQU .DELTA.V=(.vertline.S.multidot.VD.vertline.).sup.2 (4) 
 With this assumption that the value representative of the velocity 
 calculation integrity is the residue .DELTA.V defined above, the invention
 essentially relies on calculating the norm of the vector S.multidot.VD on 
 the basis of the director cosine matrix H and the individual velocity 
 measurements on the n director axes, and comparing the result with a 
 threshold to determine whether or not the measurement is acceptable. There
 is actually a risk of the threshold being exceeded whenever one satellite 
 gives a measurement which is abnormal, that is to say deviates from the 
 normal (Gaussian) statistical distribution of the measurement noise. 
 A threshold is therefore defined which is not to be exceeded for the norm 
 of the vector product S.multidot.VD. An indication that it has been 
 exceeded is given to the user in order to indicate to him that the 
 velocity measurement is not reliable. Alternatively, the user is provided 
 with a .DELTA.V value indication, and therefore an indication of the 
 quality of the measurement. 
 The integrity thresholds may, of course, vary from one application to 
 another. 
 In practice, it will more often be a "normalized" residue which is 
 calculated, that is to say one expressed in proportion to an estimated 
 measurement noise .sigma..sup.2. The normalized residue will be 
 .DELTA.V.sub.nr =.DELTA.V/.sigma..sup.2. The way in which a noise value is
 determined in order to perform this calculation will be indicated below. 
 In this case, it can be shown that an acceptable threshold T can be picked 
 for the normalized residue, this threshold depending essentially on a 
 maximum probability of false alarm which will be selected and depends on 
 the number n of satellites used. 
 An indication will be given below by way of example of a calculation method
 allowing a realistic alarm threshold T to be defined. 
 The following assumption is made in this example: the velocity measurement 
 errors are statistically distributed according to Gaussian laws. It can 
 then be demonstrated that the residue .DELTA.V.sub.nr, which can be 
 identified with a sum of squares of Gaussian random variables, follows a 
 distribution law referred to as a .chi..sup.2 law with n-4 degrees of 
 freedom. The number of degrees of freedom is the number of satellites 
 observed, n, less the number of satellites indispensable for the velocity 
 measurement, here 4 if it is assumed that four satellites are necessary 
 for giving a four-dimensional velocity vector. 
 In the case of such a Gaussian probability law, it is known to calculate a 
 probability of false alarm PFA, that is to say the probability that the 
 residue .DELTA.V.sub.nr will be greater than the threshold T even though 
 there is no velocity calculation error. The probability of false alarm is:
EQU PFA=Q(T.vertline.n-4) (5) 
 where Q is a probability law with n-4 degrees of freedom which is deduced 
 from the .chi..sup.2 probability law P with n-4 degrees of freedom by the 
 formula Q=1-P. The probability law P will be expressed by the following 
 symbolic representation: P(.chi..sup.2.vertline.n-4) 
 which gives, with the same notation 
EQU Q(.chi..sup.2.vertline.n-4)=1-P(.chi..sup.2.vertline.n-4) (6) 
 If the curves representing the .chi..sup.2 probability laws as a function 
 of T for each value n are plotted, and if a given probability PFA (of 
 false alarm) is picked, a respective threshold value T will be found for 
 each value of n which is the abscissa of the curve corresponding to the 
 value of n, for an ordinate equal to PFA. These thresholds are collated in
 a table and will be the thresholds used for the integrity check. 
 Mathematically, it may be assumed that, for a given probability of false 
 alarm PFA, the following threshold T is calculated: 
EQU T=Q.sup.-1 (PFA) (7) 
 Q.sup.-1 is the inverse function of Q(.chi..sup.2.vertline.n-4), that is to
 say if Q(x)=a, then x=Q.sup.-1 (a). 
 By way of example, having plotted the .chi..sup.2 probability curves, the 
 thresholds were calculated for an arbitrarily selected (but realistic) 
 probability of false alarm of 10.sup.-3. The normalized values of the 
 threshold T for a number n of satellites of between 5 and 10 are then the 
 following (approximately): 
 
 n = 5 6 7 8 9 10 
 T = 10 14 16 18 20 22 
 The thresholds T which are thus calculated will therefore be used according
 to the invention: when, in a velocity measurement, the matrix H, then the 
 matrix S, then the norm .DELTA.V.sub.nr of the vector product 
 S.multidot.VD are calculated, this norm will be compared with the 
 threshold T (selected as a function of the number of satellites used, from
 a table containing one threshold value for each value of n). 
 The calculations which make it possible to determine the divergence 
 .DELTA.V.sub.nr require an estimate of the velocity measurement noise. In 
 general, these parameters are predetermined on the basis of statistical 
 measurements. A realistic value for .sigma. is 0.2 m/s. In the GPS system,
 this noise actually results principally from selective availability, which
 is a variation in clock frequency of the satellites introduced to 
 intentionally degrade the position precision that can be obtained. 
 However, an estimate of the measurement noise can be obtained by a 
 statistical calculation on the basis of the velocity error residues 
 .DELTA.V, because it can be demonstrated that the mathematical expectancy 
 of the velocity error residue has a value (n-4).sigma..sup.2, n-4 being 
 the number of degrees of freedom of the velocity resolution, and .sigma. 
 being the velocity measurement noise. Following a certain number of 
 measurements, it is therefore possible to calculate measurement noise 
 statistics, and these statistics are used to calculate the normalized 
 residue .DELTA.V.sub.nr. 
 The simplest case is to calculate .DELTA.V every 200 ms, for example (the 
 delivery rate of the velocity measurements) then to calculate a moving 
 average of the terms .DELTA.V/n-4 every few minutes, and this average 
 represents the estimate of the square of the measurement noise. It 
 naturally assumes integrity of the system. Initially, the value of .sigma.
 can be arbitrarily fixed, for example at 0.2 m/s, before a moving average 
 value is available. 
 Besides the indication of the lack of integrity in the velocity measurement
 (.DELTA.V.sub.nr exceeds the threshold T), an indication of the quality of
 the integrity monitoring performed may also be provided. 
 It is actually not sufficient to confirm the measurement's integrity with a
 determined probability of false alarm. It would be better to know as well 
 what the probability of non-detection of a velocity error is. 
 This is why it is proposed, according to the invention, to fix a maximum 
 probability of non-detection which can be tolerated in the receiver, and 
 to deduce therefrom what is the velocity error limit which results 
 therefrom and beyond which the measurement will be considered as 
 unreliable. 
 For example, a maximum probability of non-detection of 10.sup.-2 is picked,
 that is to say at most one chance in 100 that the system will not see a 
 velocity error is tolerated, and this velocity error which gives rise to 
 one chance in 100 of not being detected is calculated. 
 When such a PND probability is picked, a geometrical parameter .THETA. may 
 be calculated which will itself be used to calculate a check indicator of 
 the quality of the integrity surveillance performed. 
 The following assumption may be made: it is assumed that the velocity error
 is due to a velocity bias on one of the satellites (for example resulting 
 from an error in its clock), and that the distribution of this velocity 
 bias is Gaussian; it can be demonstrated that the residue .DELTA.Vhd nr 
 follows an uncentred .chi..sup.2 distribution law with n-4 degrees of 
 freedom, of parameter .THETA. proportional to the bias, with n-4 degrees 
 of freedom. The parameter .THETA. of the uncentred probability law 
 P'(.chi..sup.2.vertline.n-4, .THETA.) is then: 
EQU .THETA.=(B.sup.2 /.sigma..sup.2)S.sub.i,i (8) 
 B is the velocity bias; 
 .sigma. is the estimate of the measurement noise standard deviation for the
 velocity; 
 S.sub.i,i is the sensitivity to velocity errors along the axis i: 
 measurement of the error generated on the axis i by a bias on the axis i 
 itself. This is the element of rank i,i of the matrix S mentioned above 
 with reference to equation (4). 
 Under these conditions, there is a probability of non-detection, that is to
 say a probability that the velocity residue .DELTA.V.sub.nr will be less 
 than the threshold T (fixed beforehand as a function of PFA and of n) even
 though there is actually a velocity calculation error due to a bias on one
 of the satellites. 
 The curves representing this probability of non-detection as a function of 
 .THETA. for various possible values of n can be plotted. If a probability 
 of non-detection is then fixed, a table of values of .THETA. for various 
 values of n is obtained. For example, if PND=10.sup.2 is picked, the 
 following table is obtained (approximate values) for n lying between 5 and
 10 satellites: 
 
 n 5 6 7 8 9 10 
 .THETA. 31.5 35 38 40 42 45 
 Now, it can be shown that .THETA. is also equal to 
 (RPE/.sigma..DELTA.H.sub.i).sup.2, where RPE is the velocity radius of 
 protection, that is to say the velocity error for which the probability of
 non-detection is equal to PND, .DELTA.H.sub.i is a difference between two 
 values of "dilution of precision" XDOP which are linked with the satellite
 constellation used; more precisely, it is the difference between the 
 dilution of precision XDOP of the full constellation of n satellites, and 
 the dilution of precision XDOP of the constellation of n-1 satellites 
 excluding the axis i. The dilution of precision XDOP is derived from the 
 matrix H.sup.T.multidot.H.sup.-1 (4.times.4 matrix) . It is the square 
 root of the sum of the terms on the diagonal of the matrix. 
 If it is primarily the horizontal velocity which is of interest, it will 
 instead be the horizontal dilution of precision HDOP that is used, which 
 is the square root of the sum of the first two terms on the diagonal. If 
 only the vertical velocity is of interest, only the third term on the 
 diagonal will be taken (VDOP). 
 The term XDOP or HDOP or VDOP represents the quality of the constellation 
 of satellites at the time of the measurement, and does not have the same 
 value according to whether all the satellites or all the satellites except
 the satellite of rank i are taken. .DELTA.H.sub.i represents the 
 difference in XDOP or HDOP or VDOP calculated in the two cases. 
 The result of all this is that, if a probability of non-detection PND is 
 picked, the values of geometrical criteria .THETA. which link the radius 
 of protection and the .DELTA.H.sub.i can be deduced therefrom. Then, as a 
 function of these values of .THETA. and the status of the satellite 
 constellation observed at a given time, it is possible to calculate the 
 radius of protection for the velocity calculation integrity surveillance: 
 this radius RPE is the maximum value found for the product of the 
 following three values: 
 the square root of .THETA., 
 .DELTA.H.sub.i, 
 .sigma.. 
 This product is calculated for all the satellites i from 1 to n, and the 
 maximum value found represents the radius of protection RPE. 
EQU RPE=max(.THETA..sup.1/2.multidot..DELTA.H.sub.i.multidot..sigma.) for i=1 
 to n 
 This radius of protection RPE, or velocity error limit for a given 
 probability of non-detection, is expressed in velocity units (for example 
 m/s) and can be displayed for the user or give rise to an alarm if it 
 exceeds a predetermined value. The alarm may be an indication of 
 "unreliable velocity measurement". 
 The invention therefore makes it possible to give both an alarm in case of 
 insufficient integrity in the velocity measurement on more than four 
 satellites and a value of velocity radius of protection. 
 When insufficient integrity is observed, it is desirable to at least 
 temporarily suppress the data from the satellite which is causing the 
 integrity loss. To do this, it is necessary to determine which this 
 satellite is. 
 To do this, the velocity error residue is calculated for all the groups of 
 n-1 satellites, and the satellite rank i, for which the residue resulting 
 from the n-1 other satellites does not cause the integrity threshold to be
 exceeded whereas the residue for n satellites causes it to be exceeded, is
 determined. 
 Once the satellite which has caused the error has been identified, the data
 from this satellite are suppressed.