METHOD FOR PREDICTING THE ORBIT OF A SATELLITE AND CORRESPONDING SATELLITE SIGNAL RECEIVER

Abstract

A signal receiver method to achieve satellite position fix by improving satellite orbit prediction includes: acquiring satellite signals and navigation data and calculating a position solution, which includes predicting the state or orbit of one or more satellites. The prediction includes using a model of the solar radiation pressure operating on a selected satellite. The method includes: expressing the model of the solar radiation pressure operating on the satellite as a Fourier series having a frequency function of the satellite-Earth-Sun angle and having respective Fourier coefficients, calculating position approximation errors comparing true satellite positions at given time points against predicted satellite positions at corresponding time points, estimating said Fourier coefficients as a function of the position approximation errors, and using the estimated Fourier coefficients in the model of the solar radiation pressure operating on the satellites as a Fourier series used in the model to predict the state or orbit of the satellite.

Claims

1. A satellite signal receiver method to predict an orbit of a satellite to achieve a satellite signal receiver position fix, comprising: acquiring satellite signals with a satellite signal receiver; acquiring navigation data with the satellite signal receiver; and calculating a position solution, wherein calculating the position solution includes: predicting at least one of a state or an orbit (r(t)) of a selected satellite that will be used by the satellite signal receiver, said prediction including using a model ({umlaut over (r)}.sub.SRP(t)) of solar radiation pressure operating on said selected satellite, the predicting including: expressing the model of solar radiation pressure ({umlaut over (r)}.sub.SRP(t)) operating on the selected satellite as a first Fourier series (ë′(t)=[ë′.sub.x.sub.s ë′.sub.y.sub.s ë′.sub.z.sub.s]) having a frequency function of an angle defined by a solar position with respect to a body fixed frame of the selected satellite and having respective first Fourier coefficients of said first Fourier series; calculating position approximation errors (a.sub.pae(t.sub.i,j)) by comparing true satellite positions r.sub.eph(t.sub.i,j) at given time points against predicted satellite positions (r.sub.pred(t.sub.i,j)) at corresponding time points; estimating said first Fourier coefficients of said first Fourier series as a function of said position approximation errors (a.sub.pae(t.sub.i,j)); and using said estimated first Fourier coefficients of said first Fourier series in the model of solar radiation pressure operating on the selected satellite expressed as a second Fourier series and used in the model of solar radiation pressure to predict the at least one of the state or the orbit of the selected satellite.

2. The method according to claim 1, wherein the satellite signal receiver is a Global Navigation Satellite System (GNSS) receiver.

3. The method according to claim 1, wherein the angle defined by the solar position with respect to the body fixed frame of the selected satellite is a satellite-Earth-Sun angle (α).

4. The method according to claim 1, wherein calculating position approximation errors (a.sub.pae(t.sub.i,j)) by comparing true satellite positions r.sub.eph(t.sub.i,j) at given time points against predicted satellite positions (r.sub.pred(t.sub.i,j)) at corresponding time points includes: obtaining a verification set of ephemerides (eph.sub.v) at given time instants; obtaining time end points of time intervals around said given time instants; obtaining a prediction set of ephemerides (eph.sub.p); using said prediction set of ephemerides (eph.sub.p) to obtain initial conditions; generating a predicted orbit as a function of said initial conditions (r.sub.0, v.sub.0); calculating predicted positions (r.sub.pred(t.sub.i,j) of said predicted orbit of the selected satellite at said time end points; obtaining true positions (r.sub.eph(t.sub.i,j)) of the selected satellite at said time end points from said verification set of ephemerides (eph.sub.v); and comparing said predicted positions (r.sub.pred(t.sub.i,j)) of said predicted orbit of the selected satellite at said time end points and true positions (r.sub.eph(t.sub.i,j)) of the selected satellite at said time end points to obtain position approximation errors (a.sub.pae(t.sub.i,j)) at said time end points.

5. The method according to claim 1, comprising: truncating the first Fourier series (ë′(t)=[ë′.sub.x.sub.s ë′.sub.y.sub.s ë′.sub.z.sub.s]) of the model of the solar radiation pressure ({umlaut over (r)}.sub.SRP(t)) at a given element (I.sub.max).

6. The method according to claim 1, wherein estimating said first Fourier coefficients of said first Fourier series as a function of said position approximation errors (a.sub.pae(t.sub.i,j)) includes choosing values minimizing said position approximation errors (a.sub.pae(t.sub.i,j)).

7. The method according to claim 1, wherein said estimating said first Fourier coefficients of said first Fourier series as a function of said position approximation errors (a.sub.pae(t.sub.i,j)) includes: selecting a coefficient among said first Fourier coefficients of the first Fourier series and setting the selected coefficient at a first given value; for a determined number of iterations, setting all other unselected first Fourier coefficients to fixed values and calculating an estimator of the approximation errors; setting the selected coefficient at a second given value; calculating a second estimator of the approximation errors to obtain a respective sequence of set values and a sequence of estimator values; fitting said sequence of set values and said sequence of estimator values for the selected coefficients to identify an estimation function (y.sub.est); and selecting a value associated with the minimum value of the estimation function as an estimated value for the selected coefficient.

8. The method according to claim 7, comprising: repeating the estimating of said first Fourier coefficients of said first Fourier series as a function of said position approximation errors a determined number of times.

9. The method according to claim 1, wherein generating the predicted orbit as a function of said initial conditions (r.sub.0, v.sub.0) includes: using an approximated second derivative of a position error (ë′(t)) between the orbit of the selected satellite that does take the model of solar radiation pressure ({umlaut over (r)}.sub.SRP(t)) into account and the orbit of the selected satellite that does not take the model of solar radiation pressure ({umlaut over (r)}.sub.SRP(t)) into account.

10. The method according to claim 1, wherein calculating position approximation errors (a.sub.pae(t.sub.i,j)) by comparing true satellite positions r.sub.eph(t.sub.i,j) at given time points against predicted satellite positions (r.sub.pred(t.sub.i,j)) at corresponding time points includes: calculating velocity approximation errors; and calculating acceleration approximation errors.

11. A satellite signal receiver, comprising: a front-end to acquire satellite signals and to acquire navigation data; and a processor arranged to calculate a position solution that includes prediction of a state or an orbit of a selected navigation satellite, wherein the prediction solution applies a model of solar radiation pressure operating on said selected navigation satellite, the position solution calculation including: expressing the model of solar radiation pressure operating on the selected navigation satellite as a first Fourier series having a frequency function of an angle defined by a solar position with respect to a body fixed frame of the selected navigation satellite and having respective first Fourier coefficients of said first Fourier series; calculating position approximation errors by comparing true satellite positions at given time points against predicted satellite positions at corresponding time points; estimating said first Fourier coefficients of said first Fourier series as a function of said position approximation errors; and using said estimated first Fourier coefficients of said first Fourier series in the model of solar radiation pressure operating on the selected navigation satellite expressed as a second Fourier series used in the model of solar radiation pressure to predict the state or the orbit of the selected navigation satellite.

12. A satellite signal receiver according to claim 11, wherein the satellite signal receiver is a Global Navigation Satellite System (GNSS) receiver.

13. A satellite signal receiver according to claim 11, wherein the position solution calculation further includes: obtaining a verification set of ephemerides at given time instants; obtaining time end points of time intervals around said given time instants; obtaining a prediction set of ephemerides; using said prediction set of ephemerides to obtain initial conditions; generating a predicted orbit as a function of said initial conditions; calculating predicted positions of said predicted orbit of the selected satellite at said time end points; obtaining true positions of the selected satellite at said time end points from said verification set of ephemerides; and comparing said predicted positions of said predicted orbit of the selected satellite at said time end points and true positions of the selected satellite at said time end points to obtain position approximation errors at said time end points.

14. A satellite signal receiver according to claim 11, wherein the position solution calculation further includes: selecting a coefficient among said first Fourier coefficients of the first Fourier series and setting the selected coefficient at a first given value; for a determined number of iterations, setting all other unselected first Fourier coefficients to fixed values and calculating an estimator of the approximation errors; setting the selected coefficient at a second given value; calculating a second estimator of the approximation errors to obtain a respective sequence of set values and a sequence of estimator values; fitting said sequence of set values and said sequence of estimator values for the selected coefficients to identify an estimation function (y.sub.est); and selecting a value associated with the minimum value of the estimation function as an estimated value for the selected coefficient.

15. A satellite signal receiver according to claim 11, wherein the position solution calculation further includes: calculating velocity approximation errors; and calculating acceleration approximation errors.

16. A satellite navigation system, comprising: a satellite signal receiver coupleable to a plurality of navigation satellites, the satellite signal receiver arranged to calculate a position solution that includes at least one prediction of a state or an orbit of a selected satellite of the plurality of navigation satellites, wherein the prediction solution applies a model of solar radiation pressure operating on said selected satellite, the position solution calculation including: expressing the model of solar radiation pressure operating on the selected satellite as a first Fourier series having a frequency function of an angle defined by a solar position with respect to a body fixed frame of the selected satellite and having respective first Fourier coefficients of said first Fourier series; calculating position approximation errors by comparing true satellite positions at given time points against predicted satellite positions at corresponding time points; estimating said first Fourier coefficients of said first Fourier series as a function of said position approximation errors; and using said estimated first Fourier coefficients of said first Fourier series in the model of solar radiation pressure operating on the selected satellite expressed as a second Fourier series used in the model of solar radiation pressure to predict the state or the orbit of the selected satellite.

17. A satellite navigation system according to claim 16, wherein the satellite signal receiver is a Global Navigation Satellite System (GNSS) receiver.

18. A satellite navigation system according to claim 16, wherein the position solution calculation further includes: obtaining a verification set of ephemerides at given time instants; obtaining time end points of time intervals around said given time instants; obtaining a prediction set of ephemerides; using said prediction set of ephemerides to obtain initial conditions; generating a predicted orbit as a function of said initial conditions; calculating predicted positions of said predicted orbit of the selected satellite at said time end points; obtaining true positions of the selected satellite at said time end points from said verification set of ephemerides; and comparing said predicted positions of said predicted orbit of the selected satellite at said time end points and true positions of the selected satellite at said time end points to obtain position approximation errors at said time end points.

19. A satellite navigation system according to claim 16, wherein the position solution calculation further includes: selecting a coefficient among said first Fourier coefficients of the first Fourier series and setting the selected coefficient at a first given value; for a determined number of iterations, setting all other unselected first Fourier coefficients to fixed values and calculating an estimator of the approximation errors; setting the selected coefficient at a second given value; calculating a second estimator of the approximation errors to obtain a respective sequence of set values and a sequence of estimator values; fitting said sequence of set values and said sequence of estimator values for the selected coefficients to identify an estimation function (y.sub.est); and selecting a value associated with the minimum value of the estimation function as an estimated value for the selected coefficient.

20. A satellite navigation system according to claim 16, wherein the position solution calculation further includes: calculating velocity approximation errors; and calculating acceleration approximation errors.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Non-limiting and non-exhaustive embodiments are described with reference to the following drawings, wherein like labels refer to like parts throughout the various views unless otherwise specified. One or more embodiments are described hereinafter purely by way of a non-limiting example with reference to the annexed drawings, in which:

[0032] FIG. 1 is a flow diagram showing operations of the method here described;

[0033] FIG. 2 represents quantities used by the method here described.

DETAILED DESCRIPTION

[0034] The ensuing description illustrates various specific details aimed at an in-depth understanding of the embodiments. The embodiments may be implemented without one or more of the specific details, or with other methods, components, materials, etc. In other cases, known structures, materials, or operations are not illustrated or described in detail so that various aspects of the embodiments will not be obscured.

[0035] Reference to “an embodiment” or “one embodiment” in the framework of the present description is meant to indicate that a particular configuration, structure, or characteristic described in relation to the embodiment is comprised in at least one embodiment. Likewise, phrases such as “in an embodiment” or “in one embodiment”, that may be present in various points of the present description, do not necessarily refer to the one and the same embodiment. Furthermore, particular conformations, structures, or characteristics can be combined appropriately in one or more embodiments.

[0036] The references used herein are intended merely for convenience and hence do not define the sphere of protection or the scope of the embodiments.

[0037] Now the method here described will be detailed.

[0038] As already mentioned, the four main forces acting on artificial satellites (like GNSS satellites) are Earth gravity, Sun gravity, Moon gravity and solar radiation pressure.

[0039] Therefore, the overall acceleration {umlaut over (r)}(t) acting on the satellite can be approximated to the sum {umlaut over (r)}.sub.a(t) of the Earth gravity {umlaut over (r)}.sub.Earth(t) the Sun gravity {umlaut over (r)}.sub.Sun(t) and Moon gravity {umlaut over (r)}.sub.Moon(t) and solar radiation pressure acceleration {umlaut over (r)}.sub.SRP(t), because the other forces acting on the satellite are negligible for the purpose of the method here described. Overall acceleration {umlaut over (r)}(t) is therefore.

{umlaut over (r)}(t)={umlaut over (r)}.sub.a(t)={umlaut over (r)}.sub.Earth(t)+{umlaut over (r)}.sub.Sun(t)+{umlaut over (r)}.sub.Moon(t)+{umlaut over (r)}.sub.SRP(t)  (1)

[0040] Satellite position as a function of time t, r(t), i.e., the satellite orbit, can be found by integrating the equation (1) expressing the overall acceleration {umlaut over (r)}(t) of the satellite twice:

r(t)≅r.sub.a(t)=r.sub.0+v.sub.0t+∫∫.sub.t.sub.0.sup.t({umlaut over (r)}.sub.Earth(t)+{umlaut over (r)}.sub.Sun(t)+{umlaut over (r)}.sub.Moon(t)+{umlaut over (r)}.sub.SRP(t)dt  (2)

r.sub.0 and v.sub.0 being initial conditions on position and velocity. If it is not taken into account the solar radiation pressure {umlaut over (r)}.sub.SRP(t) operating on the satellite, it is obtained the following approximation r′(t) for the satellite position:

r′(t)=r.sub.0+v.sub.0t+∫∫.sub.t.sub.0.sup.t({umlaut over (r)}.sub.Earth(t)+{umlaut over (r)}.sub.Sun(t)+{umlaut over (r)}.sub.Moon(t))dt  (3)

[0041] Hence, from the equation (3) a position error e(t) can be calculated, as difference of the satellite position as a function of time t, r(t), and of approximation of the satellite position r′(t) not taking into account the solar radiation pressure {umlaut over (r)}.sub.SRP(t) operating on the satellite, that can be used to approximate the missing solar radiation pressure term:

e(t)=r(t)−r′(t)≅r.sub.a(t)−r′(t)=r.sub.SRP(t)=r.sub.SRP,0+v.sub.SRP,0t+∫∫.sub.t.sub.0.sup.t({umlaut over (r)}.sub.SRP(t))dt  (4)

[0042] Therefore, the second derivative of the position error e(t) can be used as good approximation of the solar radiation pressure term of the position r(t):

{umlaut over (e)}(t)≅{umlaut over (r)}.sub.SRP(t)  (5)

[0043] Equations (1) (5) detail the relations between quantities used in the method here described

[0044] Now, the method for predicting the orbit of a satellite will be described in more detail with reference to the diagram flow of FIG. 1, which depicts possible operations of such method.

[0045] Since it is known that solar radiation pressure acceleration {umlaut over (r)}.sub.SRP(t) (and therefore the second derivative of the position error ë(t)) can be approximated to a periodic function of the satellite-Earth-Sun angle α if a short period (e.g., few days) is considered, the method provides a first step 110 of expressing the model of the solar radiation pressure, i.e., the acceleration {umlaut over (r)}.sub.SRP(t)) operating on the satellite due to the solar radiation as a Fourier series (ë′(t)=[ë′.sub.x.sub.s ë′.sub.y.sub.s ë′.sub.z.sub.s]) having a frequency function of the satellite-Earth-Sun angle α and having respective Fourier coefficients CF.

[0046] It has to be specified that the angle α can be any angle defined by the solar position with respect to the satellite body fixed frame. In particular, it can be the angle between the direction to which the satellite body (excluding the solar arrays) is pointed (i.e., the direction of an arbitrary axis fixed in the satellite body) and the direction to which the satellite solar arrays are pointed (i.e., toward the Sun). In case of antennas which are integral with the body of the satellite, it can be chosen as the angle between the direction to which the antennas point and the direction to which the satellite solar arrays are pointed. Since the face containing the antennas of a GNSS satellite body is always pointed toward the Earth then the satellite-Earth-Sun angle can been chosen in this case. Of course, a different choice may be needed in the case of non-GNSS satellites.

[0047] In other words, it can be used a Fourier series ë′(t)=[ë′.sub.x.sub.s ë′.sub.y.sub.s ë′.sub.z.sub.s] to approximate the second derivative of the position error ë(t)=[ë.sub.x.sub.s ë.sub.y.sub.s ë.sub.z.sub.s], as follows:

ë.sub.x.sub.s≅ë′.sub.x.sub.s=A.sub.0+Σ.sub.l=0.sup.∞(A.sub.cl cos α+A.sub.sl sin α)

ë.sub.y.sub.s≅ë′.sub.y.sub.s=B.sub.0+Σ.sub.l=0.sup.∞(B.sub.cl cos α+B.sub.sl sin α)

ë.sub.z.sub.s≅ë′.sub.z.sub.s=C.sub.0+Σ.sub.l=0.sup.∞(C.sub.cl cos α+C.sub.sl sin α)  (6)

[0048] Equation (6) expresses the components ë.sub.x.sub.s, ë.sub.y.sub.s, ë.sub.z.sub.s of the second derivative of the position error ë(t) along the threes axis x.sub.s, axis y.sub.s, axis z.sub.s. I is the index of the Fourier coefficients, varying between zero and infinite, the set of Fourier coefficients CF includes cosine coefficients A.sub.cl, B.sub.cl, C.sub.cl, and sine coefficients A.sub.sl, B.sub.sl, C.sub.sl as well as constants A.sub.0, B.sub.0, C.sub.0.

[0049] Axis x.sub.s by way of example is oriented along the boresight direction of the solar panels (i.e., is pointed toward the Sun) of the satellite. Axis y.sub.s is parallel to the rotation axis of the satellite's solar panels. Axis z.sub.s is orthogonal to the plane defined by x.sub.s, y.sub.s completing the right handed system.

[0050] It is also possible to choose a different reference system for axes x.sub.s, y.sub.s and z.sub.s, where axes x.sub.s and z.sub.s lies on the satellite-Earth-Sun plane and y.sub.s is orthogonal to that plane, for instance.

[0051] As indicated in FIG. 1, it is preferable to perform then a step of truncation 120 of the Fourier series expressed by the equations (6), choosing a maximum value I.sub.max for the index I of the series.

[0052] Then the method here described, after obtaining a set of Fourier coefficient CF of the Fourier series ë.sub.x.sub.s, ë.sub.y.sub.s ë.sub.z.sub.s, such set CF including Fourier coefficients A.sub.0, A.sub.cl, A.sub.sl, B.sub.0, B.sub.cl, B.sub.sl, C.sub.0, C.sub.cl, C.sub.sl, with I=0 . . . I.sub.max, which value is not yet assigned, includes a procedure 200 to evaluate the values of such Fourier coefficients in the set CF, or part of them, as better detailed in the following, which better fit the pressure radiation model, in particular the second derivative of the error ë′(t).

[0053] To this end the procedure 200 includes first a procedure 300 of calculating position approximation errors α.sub.pae comparing true satellite positions at given time points against predicted satellite positions at the same given time points.

[0054] Such procedure 300 includes a step 310 of obtaining two sets of the satellite's ephemerides, containing at least one ephemeris each: [0055] a first set, indicated as prediction set P of ephemerides eph.sub.p, is obtained. Such prediction set P is used in following steps 312 and 314 to obtain an orbit prediction r.sub.e(t); [0056] a second set, indicated as verification set V of ephemerides eph.sub.v, is used to obtain true positions and in a step 320 to measure the accuracy of a such predicted orbit.

[0057] m indicates the number of ephemerides contained in the prediction set and n the number of ephemerides contained in the verification set. As said m and n are greater than or equal to one.

[0058] Ephemerides from both sets, prediction set P and verification set V, have to be enough close in time, i.e., the time of each ephemeris eph.sub.v in the verification set V has to be in the range where the prediction generated from one or more ephemerides eph.sub.p in the prediction set can be considered accurate.

[0059] For what regards the prediction set of ephemerides eph.sub.p, in a step 312 is used to calculate the initial conditions r.sub.0 and v.sub.0.

[0060] Then in a step 314 it is generated an orbit prediction r.sub.e(t) using the equation (2), in which are used such initial conditions r.sub.0 and v.sub.0 and the approximated second derivative error ë′(t) instead of solar radiation pressure term {umlaut over (r)}.sub.SRP(t):

r.sub.e(t)=r.sub.0+v.sub.0t+∫∫.sub.t.sub.0.sup.t({umlaut over (r)}.sub.Earth(t)+{umlaut over (r)}.sub.Sun(t)+{umlaut over (r)}.sub.Moon(t)+{umlaut over (e)}′(t))dt  (7)

[0061] As indicated in the following the orbit prediction r.sub.e(t) at step 314 is calculated at specific time points t.sub.i,j obtained through the verification set of ephemerides eph.sub.p.

[0062] Such time points t.sub.i,j are calculated from the verification set in a step 316 as follows. Defining a time when the prediction begins (i.e., the initial conditions time) as t.sub.0, and t.sub.1, t.sub.2 . . . t.sub.n times of the different verification ephemerides eph.sub.v in the verification set V i.e., verification times t.sub.i (as mentioned n can be 1, thus there is only one verification time t.sub.1), for each ephemeris eph.sub.v in the verification set, it can be defined a set of a number p of points in time, t.sub.i,1 . . . t.sub.i,p, around each verification time t.sub.i, which are included in that ephemeris time window ETW. Thus i is the index of the verification times while t.sub.i, is the index of the end points around a given verification time t.sub.i.

[0063] For ephemeris time window ETW, it is meant the period in which the ephemeris is accurate. For instance, for the GPS broadcast ephemeris lasts 4 hours and it is centered on the time t.sub.i of the ephemeris itself.

[0064] Thus the method includes, after obtaining the verification set at step 310, of obtaining said p points in step 316, whose quantities are also in part indicated in FIG. 2. In said step 316, for example, each time window ETW can be divided in equal steps SP, and a time point t.sub.i,j can be set for each end point of every time interval defined by the step SP (of course for overlapping ends it is set only one point). As shown in the following in the example with reference to FIG. 2, a window ETW of 4 hours can be divided in steps SP of 15 minutes). As another example, for instance, it is possible to divide a 240-minute GPS broadcast ephemeris window ETW in steps SP of 30 minutes each, so that it is obtained a set of 9 time points t.sub.i,j:

[0065] t.sub.i,1=t.sub.i−120 minutes,

[0066] t.sub.i,2=t.sub.i−90 minutes,

[0067] t.sub.i,3=t.sub.i−60 minutes,

[0068] t.sub.i,4=t.sub.i−30 minutes,

[0069] t.sub.i,5=t.sub.i,

[0070] t.sub.i,6=t.sub.i+30 minutes,

[0071] t.sub.i,7=t.sub.i+60 minutes,

[0072] t.sub.i,8=t.sub.i+90 minutes,

[0073] t.sub.i,9=t.sub.i+120 minutes

[0074] where t.sub.i is the verification time of that ephemeris eph.sub.v.

[0075] After step 316, which produces time points t.sub.i,j, it is possible to calculate in the step 314 the position predictions for each of those time points t.sub.i,j, indicated as r.sub.pred(t.sub.i,j).

[0076] Then it is used a position from the ephemerides eph.sub.v in the verification set to obtain in a step 318 the true position of the satellite's spacecraft at the same time points t.sub.i,j—indicated as r.sub.eph(t.sub.i,j). In this case positions are not calculated, but directly taken from the verification ephemerides at the given time coordinate, i.e., time point t.sub.i,j.

[0077] Then in a step 320 position approximation errors a.sub.pae(t.sub.i,j) at said time end points t.sub.i,j are obtained comparing said predicted positions (r.sub.pred(t.sub.i,j)) and true positions (r.sub.eph(t.sub.i,j)).

[0078] The comparing function of step 320 which produces the position approximation errors a.sub.pae can be defined in different ways, for example in one embodiment here described is a difference function:

a.sub.pae(t.sub.i,j)=|r.sub.pred(t.sub.i,j)−r.sub.eph(t.sub.i,j)|

[0079] Therefore, summarizing, after last step 320 of the procedure 300, it is calculated the true satellite positions for all the time points defined as described above, and then they are compared against the positions predicted using the equation (7) for the same points.

[0080] The results from such comparison is used in a procedure 400 of estimation of the Fourier coefficients CF of the equation (6): A.sub.0, A.sub.cl, A.sub.sl, B.sub.0, B.sub.cl, B.sub.sl, C.sub.0, C.sub.cl, C.sub.sl.

[0081] In other words, coefficients CF of the equation (6) are chosen which minimize the position approximation errors a.sub.pae.

[0082] Alternatively, it is also possible calculate the predicted satellite velocity and acceleration along with the position, use them for the definition of the approximation errors, and then compare them to the corresponding true satellite data from the ephemerides. In this latter case, coefficients CF of the equation (6) are chosen which minimize either the position approximation errors, or the velocity approximation errors, or the acceleration approximation errors, or a combination of them.

[0083] To perform the estimation 400, in coefficient selection step 411 a first coefficient among Fourier coefficients CF is chosen, for instance coefficient A.sub.0. In the list of coefficients CF, index w going from 0 to a maximum value w.sub.max is used for sake of clarity to indicate the Fourier coefficients in order, i.e., CF=(CF.sub.0, . . . CF.sub.1 . . . CF.sub.wmax). In equation (6) with 9 coefficients if I.sub.max=0 is chosen, then w.sub.max=8 and w=0 . . . 8, and CF.sub.0 corresponds to A.sub.0, CF.sub.1 to A.sub.cl and so on till CF.sub.8 which corresponds to C.sub.cl. The maximum value of the index I of the series, I.sub.max, can be chosen as any value greater or equal than zero.

[0084] Then in a set 412 a value x.sub.1 is set for the selected coefficient CF.sub.w and all other Fourier coefficients CF are set to fixed values. Then in a step 413 it is calculated an estimator y.sub.est,1 of the approximation errors. For instance, it is possible to use the average modulus of the position approximation errors a.sub.pae for all the points (t.sub.1,1 . . . t.sub.n,p) of a verification set as estimator y.sub.est,1 of the approximation errors:

[00001]$y_{\mathrm{est},1}=\frac{\underset{i=1}{\overset{n}{.Math.}}.Math..Math.\underset{j=1}{\overset{p}{.Math.}}.Math..Math.a_{\mathrm{pae}}\left(t_{i,j}\right)}{p.Math.n}$

[0085] Step 412 and 413 are then repeated for k−1 times, selecting at step 411 another different value (x.sub.2 . . . x.sub.k) for the selected coefficient CF.sub.w. k is greater than or equal to three, k≧3. Then it is calculated the correspondent estimator (y.sub.est,2 . . . y.sub.est,k) for each of those values. k is the index of the calculated estimator values for different k value of the selected Fourier coefficient.

[0086] The procedure including step 411 and the k times repeated steps 412, 413 is indicated as a whole as procedure 410.

[0087] k≧3 different values for the selected coefficient, each one with its own estimator value, are obtained. In other words as output of procedure 410 it is obtained a sequence X of set values x.sub.1, x.sub.2 . . . x.sub.k for the selected coefficient CF.sub.w under estimation, X=(x.sub.1, x.sub.2 . . . x.sub.k)—and a sequence Y of estimator values y.sub.est,1, y.sub.est,2 . . . y.sub.est,k associated to them, Y=(y.sub.est,1, y.sub.est,2 . . . y.sub.est,k).

[0088] Then, such sequence X of set values and sequence Y of estimator values are fitted to find (i.e., identify) an estimation function v.sub.est=f(x) which express the estimator in function of the coefficient values. For instance, a least squares approach can be used to find a second-order polynomial curve fitting those data.

[0089] In step 420 it is found a value x.sub.min associated to the minimum value of the estimation function v.sub.est=f(x). This value x.sub.min corresponds to an estimated value CF*.sub.w for the selected coefficient CF.sub.w.

[0090] Block 430 indicates the repetition of the procedure 400 for each coefficient CF.sub.w it is desired to estimate, producing a corresponding estimated coefficient CF*.sub.w for each selected coefficient CF.sub.w.

[0091] When repeating the procedure 410 the fixed values to which set in step 412 the other Fourier coefficients of the set CF with index w different from the index of the coefficient selected have to be chosen. For example, for each coefficient, the value could be zero or it could be the estimated value if that coefficient has already undergone the estimation process corresponding to procedure 400.

[0092] After all the coefficients CF* have been estimated, in a step 500 such estimated Fourier coefficients CF* are introduced in the model of the solar radiation pressure operating on the satellites as Fourier series used to predict the state or orbit of the satellite, according to equation (6).

[0093] Subsequently, the state or orbit r(t) of the satellite is predicted using a model {umlaut over (r)}.sub.SRP(t) of the solar radiation pressure operating on said one or more satellite which corresponds to the model obtained at step 500.

[0094] Finally, more sets of prediction and verification ephemerides can be used together to improve the accuracy of the coefficients estimation. In this latter case, all the time points are joined and all the approximation errors from all the verification sets are joined before calculating the estimators.

[0095] An example of implementation will be given in the following, also with reference to the quantities shown in FIG. 2.

[0096] For example, two GPS ephemerides are considered: the first one is put into the prediction set P (i.e., it is used to predict the positions of the spacecraft) and the second one is put into the verification set V (i.e., it is used to verify the prediction in its own four hour time window).

[0097] Then, it is needed to choose which coefficients CF of the equation (6) which have to been used.

[0098] A non-limiting example of choice is the following one: [0099] to choose a small maximum value I.sub.max for the Fourier index I, for instance 0, to reduce the computational effort, [0100] to set the Fourier components f.sub.y along axis y.sub.s, parallel to the rotation axis of the satellite's solar panels, f.sub.y.sub.s=0 (i.e., setting all B coefficients to 0 in equation (6)) since the component of the solar pressure orthogonal to the satellite-Earth-Sun plane is negligible due to the satellite's geometry and attitude control.

[0101] In that example case, the Fourier series expressing the pressure radiation model in equations (6) can be rewritten as:

ë.sub.x.sub.s=A.sub.0+A.sub.c0 cos α+A.sub.s0 sin α

ë.sub.y.sub.s=0

ë.sub.z.sub.s=C.sub.0+C.sub.c0 cos α+C.sub.s0 sin α  (8)

[0102] Then the coefficients CF of equations (8) are estimated, which in this example are six. It must be noted that since acceleration is equal to the ratio of force to mass, then it is not needed to estimate the satellite mass separately, but the estimated coefficients already take it into account.

[0103] A first coefficient CF.sub.w to be estimated is selected in step 411 and the other ones are set to zero. A good choice, as already mentioned, can be the A.sub.0 coefficient, since it can be expected that solar radiation pressure will be higher in the axis x.sub.s direction because it is parallel to the Sun-satellite line and perpendicular to the solar panels.

[0104] A sequence X.sub.A0 of values in a predefined interval with a predefined step is selected for the A.sub.0 coefficient, and the first value of the sequence is assigned in step 412 to A.sub.0.

[0105] As shown in FIG. 2, which schematically shows a satellite S, sending satellite signals and navigation data D to a GNSS receiver R, and its orbit predicted r.sub.e(t) from prediction ephemeris (FIG. 2 is not to scale), the predicted positions r.sub.pred(t.sub.i,j)r.sub.pred(t.sub.i,j) are compared against the true positions and true positions r.sub.eph(t.sub.i,j)r.sub.eph(t.sub.i,j) from the verification ephemeris over the predefined points of the time window ETW, which, in this case, has been divided into 15 minute steps SP, and the position errors a.sub.pae(t.sub.i,j) are calculated. Of course, if more verification ephemerides are used, then the comparison procedure described above will be repeated for each verification ephemeris eph.sub.v. That will allow the improvement of the accuracy of the Fourier coefficients estimation. FIG. 2 shows a magnification of the time window ETW to better appreciate the steps SP and the position errors a.sub.pae(t.sub.i,j).

[0106] Then the estimator Y of the errors over the verification set window ETW (or windows, if they are more than one) is calculated in step 413 and stored along with the value of the coefficient used to obtain it. In this example, it is chosen the average modulus of the calculated position errors a.sub.pae(t.sub.i,j) of the satellite as estimator.

[0107] The above two steps 412, 413 are repeated assigning to the Fourier coefficient A.sub.0 each other solar pressure coefficient value of X.sub.A0 sequence (one at time). As mentioned, sequence X.sub.A0 has to contain at least three terms (of course, the accuracy of the estimation will increase as the number of terms grows).

[0108] At the end of procedure 410, for the A.sub.0 coefficient of equation (8) it is obtained a finite sequence of coefficient values X.sub.A0, and an associated sequence of error estimators Y.sub.A0.

[0109] Now it can be used in step 410 the least squares method to find a second-order polynomial curve fitting to find an estimate sequence y.sub.est=f(x) function, which fits the error sequence Y.sub.A0. The value x.sub.min which correspond to the minimum value of that curve (i.e., min(y.sub.est)=f(x.sub.min)) is the estimated value for the coefficient.

[0110] Once the A.sub.0 coefficient has been estimated, it is possible to set such estimated value in equation (8) and estimate all other coefficients of the equation (8) in the same way, one at a time, repeating the procedure 400 under the control of repetition block 430.

[0111] For instance, it is possible to continue the coefficients estimation in the following way: [0112] set coefficient A.sub.0 to the value that has been previously estimated and estimate coefficient C.sub.0 with same procedure, while A.sub.c0, A.sub.s0, C.sub.c0, C.sub.s0 are set to zero; [0113] set A.sub.0 and C.sub.0 to the values that have been previously estimated and estimate A.sub.c0 with the same procedure, while A.sub.s0, C.sub.c0, C.sub.s0 are set to zero. [0114] set A.sub.0 and C.sub.0 to the values that have been previously estimated and estimate A.sub.s0 with same procedure, while A.sub.s0, C.sub.c0, C.sub.s0 are set to zero; [0115] set A.sub.0, C.sub.0, A.sub.c0 and A.sub.s0 to the values that have been previously estimated and estimate C.sub.c0 with same procedure, while C.sub.s0 is set to zero; [0116] set A.sub.0, C.sub.0, A.sub.c0 and A.sub.s0 to the values that have been previously estimated and estimate C.sub.s0 with same procedure, while C.sub.c0 is set to zero.

[0117] At the end of these operations, i.e., at the end of procedure 400, all the coefficients CF=(A.sub.0, A.sub.c0, A.sub.s0, C.sub.0, C.sub.c0, C.sub.s0) of the equations (8) are estimated, i.e., is produced a corresponding estimate CF*.

[0118] The following tests have been executed using the solar pressure model detection described in the previous example. Moreover the coefficients of the model have been calculated as described in the very same example.

[0119] A first test compares the prediction accuracy of the GPS block IIR satellite solar radiation pressure model proposed by Y. Bar-Sever and D. Kuang in “New Empirically Derived Solar Radiation Pressure Model for Global Positioning System Satellites” JPL Interplanetary Network Progress Vol. 24, 159, November 2004 to the accuracy of the method proposed here using a GPS block IIR satellite.

[0120] A second test compares the prediction accuracy of the same Bar-Sever and Kuang GPS block IIR satellite solar radiation pressure model to the accuracy of the method proposed here using a GPS block IIF satellite (simulating, in that way, the substitution of an IIR satellite with an IIF one).

[0121] In both tests, six different broadcast ephemerides from the same day have been used for verification. Moreover, results from seven different days have been used together to esteem the coefficients of the model.

[0122] Then the maximum and average position errors modulus from 12 orbit predictions from both models have been compared.

[0123] The prediction method used for all the tests is the same, except for the solar radiation pressure model used, of course.

[0124] To highlight the differences due to the adoption of a solar radiation pressure model expressed and calculated as Fourier series, the results are reported as ratio between the maximum and average errors of the proposed model and the maximum and average errors of the Bar-Sever and Kuang model. Therefore, values lower than 100% indicates that the method here disclosed provides lower errors.

[0125] The following Table I shows the results for such ratio (second column) corresponding to average and maximum position errors (first column).

[0126] The results from the first test are:

TABLE-US-00001 TABLE I Position error Ratio Average 93.06% Maximum 84.84%

[0127] As it can be seen, the proposed model provides a lower value of the average and maximum position errors with respect to the Sever-Kuang model in the case of correct identification of the satellite block. The results from the second test are shown herebelow in Table II:

TABLE-US-00002 TABLE II Position error Ratio Average 94.17% Maximum 77.53%

[0128] As it can be seen, the model detection provide a lower average and maximum errors than Sever Kuang model in the case of wrong identification of the satellite block, e.g., when an IIR block satellite has been substituted with a newer IIF block satellite and the prediction algorithm is not updated for the new satellite block.

[0129] The method according to the various embodiments here described presents one or more advantages with respect to previous solutions that are easily applicable in practice, i.e.: [0130] it works with a receiver that is either connected or not connected to the Internet;

[0131] it works with a receiver that does not need to operate all the time, i.e., can be switched off for some time; and [0132] it needs just a small amount of data to be stored into the receiver.

[0133] In particular, the method according to the various embodiments here described needs only to collect periodically two or more ephemerides for each satellite for a certain period. For instance (but not limited to) it can be fed with four ephemerides per day for 7 days.

[0134] The method here described is a stand-alone method to detect dynamically the spacecraft type and therefore to avoid the progressive loss of performance in orbit prediction without any external update of the receiver, even in case new types of satellites, with unknown solar pressure models, are put into use.

[0135] Of course, without prejudice to the principle of the embodiments, the details of construction and the embodiments may vary widely with respect to what has been described and illustrated herein purely by way of example, without thereby departing from the scope of the present embodiments, as defined the ensuing claims.

[0136] The various embodiments described above can be combined to provide further embodiments. These and other changes can be made to the embodiments in light of the above-detailed description. In general, in the following claims, the terms used should not be construed to limit the claims to the specific embodiments disclosed in the specification and the claims, but should be construed to include all possible embodiments along with the full scope of equivalents to which such claims are entitled. Accordingly, the claims are not limited by the disclosure.