- Email: [email protected]
A simple analytic method for the evaluation of the performance of a radiolocation system
Acta AstronauticaVol. 15, No. 1, pp. 35M4, 1987 Printed in Great Britain. All rights reserved
0094-5767/87 $3.00+ 0.00 © 1987 Pergamon Journals Ltd
A SIMPLE ANALYTIC METHOD FOR THE EVALUATION OF THE PERFORMANCE OF A RADIOLOCATION SYSTEMt O. ZARROUATI Centre National d'Etudes Spatiales (CNES), 18 Av. E. Belin, 31400 Toulouse, France
(Received 11 February 1986)
Abstract--The problem is, knowing the distribution over the world of a set of radiolocation beacons, and the individual performance of any beacon, to evaluate the performance of the global system so obtained for the localization of a satellite on a particular orbit. Under some very realistic assumptions on the orbit, the author gives an analytical formulation of the covariance of the computation of the satellite position. The formulation is illustrated in one particular case: localization of a GEO with a ground network (distance and angular measurements), but may be used in much different situations. The results take into account the "white" measurement noise and the "colored" modelization noise. The method can be very helpful in mission analysis for the optimal dimensioning of a radiolocation system.
i. INTRODUCTION 1.1 Hypotheses The method presented here allows the size and form of the geostationary spacecraft ranging error volume to be evaluated. The ranging hypotheses and conditions are as follows: (i) Range and angular measurements (elevation, bearing, Cardan a n g l e s , . . . ) only. (ii) Unlimited number of Earth tracking stations. (i~i) Orbit calculation is performed (batch processing) using a set of measurements obtained over a given period of time (a few days: this period is called an "adjustment arc"). (iv) The spacecraft position and velocity (or its orbital parameters, which is equivalent) are adjusted at a given date (here: the middle of the measuring period) so that the resulting ephemeris "fit best" through the measurements of the period. The method used is that of the least squares fit. (v) The orbit is near geostationary. That is, the spacecraft and the tracking stations relative geometry does not vary during the adjustment arc. The spacecraft motion is Keplerian.
ex = e cos(co + ~): ey = e sin(o) + f~): components of the "eccentricity vector" e, which belongs to the equatorial plane and is directed towards the perigee. Its modulus is the eccentricity e. hx = i cos ~: hy = i sin ~: components of the "inclination vector" h, which belongs to the equatorial plane and is directed about the rotation axis of the orbital plane with respect to the equatorial plane. Its modulus is the orbit inclination i. 2 = M + co + f~: spacecraft mean longitude counted from the reference axis 07 (Vernal axis). During orbit calculation we know the nominal orbital parameters (a0, e0, ho, 2o) at date to but we do not know the increments (6a, fie, fib, 62) by which these parameters must be modified to "fit best through the measurements" in the least squares fit sense. Conversely, we know how to evaluate, at any date • according to Keplerian motion hypothesis (v), • using the definition of the adapted parameters, the impact of these orbital parameters increments upon the spacecraft positioning. It is a positioning deviation 6r given by:
fir = afh ^ U M - a b e - (VMV r ) are
1.2 Orbital-dynamic parameters of the problem
+ VM[ar~. -- 3 (2 -- 20)6a] + UM 6a
Parameters adapted to the geostationary orbit are used. These parameters are expressed, using conventional Keplerian parameters (a, e, i, co, ~, M), by:
(1)
with: U M = u n i t vector directed about the " E a r t h ' s center-mean spacecraft" axis, k = unit vector directed about the pole axis towards the north, V M = k ^ UM.
a: semi-major axis tPaper presented at the 36th Congress of the International Astronautical Federation, Stockholm, Sweden, 7-12 October 1985. 35
36
O. ZARROUATI
In the expression of fir, the following term appears: --3 VM(2 -- 20)~a; this term does not directly derive from the orbital parameters definition. This is a term induced by Kepler's 3rd law. The expression ( 2 - 20) that appears is the longitude excursion (in radians) followed by the geostationary spacecraft as from the adjustment date. A date when its longitude is 20. Equation (1) immediately gives the projection relations: (6r UM) = -U~t(afe) + 6a,
(2.1)
(fir VM) = -2V~t(afe) + a62 _ 3 (2 - 20)6a, (2.2) (fir k) = - V r ( a f h ) .
(2.3)
1.3 Measurement modelling Any range measurement between the spacecraft and a station and any angular measurement can be modelled by a form: rues = (6r.I)
(3)
where I is a unit vector depending on the type of measurement. • The range measurement is a projection of the deviation 6r on the "station-spacecraft" axis. • An angular measurement is equivalent to a projection of the deviation fir on an axis normal to the "station-spacecraft" axis. Formulation (3) thus allows any type of measurement (range or angular) to be assigned to a standard deviation a in distance units (in meters). This formulation is a rather general approximation which only relies on the assumption of the spacecraft/tracking station distance, a long distance with respect to a. But the geostationary orbit benefits from an additional property: For each tracking station and for each type of measurement, I is fixed within the Earth's reference frame. We may also write: F ) = l~i) U u + 1~i) ~2 v V M + l~i) k,
1]~, l~~, l~~)=constants,
l]m + l~i~2 + l~m = 1;
(4)
(i) is the generic index of a type of measurement on a given tracking station. Finally, each measuring instrument of each tracking station is defined by 4 parameters: I(i) 1(i) I(i) O-0)). ~1 , ~ 2 , ~ 3 ,
Appendix 1 shows how these parameters are computed for elevation, bearing and range measurements as a function of the station coordinates (2s, ~b3, of the spacecraft positioning longitude 2p (longitudes 2~ and 2p are calculated within the Earth's reference frame,
unlike the spacecraft longitude 2 on its orbit, calculated from the vernal axis). The partial derivatives of the measurements obtained in relation to the adjusted parameters (also called state vector parameters) are directly inferred from formulae (2) to (4): - 8
mesj c36a
8mesj a86e
Bj =
8mesj at36h
l, - ~ 6 (~ - ;.o)~
- Ii U us - 2/2 VMj (5) - 6vMj
8mesj a062
6
• ~. is the spacecraft longitude on its orbit at the date when the measurement: mesj was carried out. The spacecraft is then on the axis directed by Uuj. Its velocity is along V~j. • The measuring instrument which performed mesj is now characterized by values (ll, 12,/3). Some important orders of magnitude can be inferred from the calculations given in Appendix 1. • For the range measurement, we have 1 ~ 1 , ~12, l 3.
• For the angular measurement: 12 ~13 ~ l I n 0 ;
These figures result from the high value of the geosynchronous semi-major axis with respect to the Earth's equatorial radius. 1.4 Some results o f the least squares theory
Orbit modelling and measurement simulation depend on a certain number of unknown or uncertain parameters. Unknown parameters will be identified by the least squares fit method--this is the state vector (normalized by a): X = (6a, a6e, a f h , a62).
Errors will be made in the knowledge of uncertain parameters. As these parameters are not adjusted (one cannot adjust all of them or the observability problem will be impaired), errors made in their evaluation will inevitably impair identification of X. The set of uncertain parameters is noted 7: 7 = (calibration bias, station coordinates, e t c . . . ) One writes respectively B and K as the partial derivatives matrices for all the measurements carried out in relation to the X state vector parameters and to the uncertain parameters.
Evaluation of the performance of a radiolocation system
37
That is, by using, in addition, formula (5) notations: 8mes~
B=
8mes~
BI
8X .
87
= B y , K =
c3mesn
8X
mesj = B j X - Kj7 - Ej,
(7)
Ej= creation of a Gaussian white noise with a standard deviation a for measurement j. By incorporating all the measurements (mesj)j = 1, n in a single-column vector M: m = [mes, . . . . mesj . . . . mes,] r. The X state vector identification problem is reduced to the resolution of a redundant linear system: (8)
B Y = M.
If the weights (reliability) attributed to each measurement are grouped together within a n × n diagonal matrix W, the "best solution" of the system (8) in the least squares fit sense is known: (9)
~B T W M .
If we knew Ej and ? (which is not the case), we could directly evaluate, by means of (7), the error made by the least squares (9) estimator: - X=(BrWB)
I BTW(BX_
= ( B r W B ) -l B r W ( - E
K?_E)_
K.
(6)
87
One obtains a set of n measurements (mesj) j = 1, n which are reduced to position deviation 5r projections on axes I(i). The non-zero value of these measurements (existence of a deviation in position) may be due to the state vector (the real orbit is not the nominal orbit) or to a lack of knowledge of the parameters 7- Linearization enables each measurement to be written in the form:
= (BrWB)
.
c3mes,
o.
KI
X
- K?).
(10)
In fact, we only know:
This covariance is calculated under the hypothesis: E(E~ q = 0. It is the sum of two terms: • A so-called "Measurement noise" term which describes the effect of measurement error 57 on the adjustment of the state vector. • A so-called "Model noise" which describes the effect of lack of knowledge of y on the adjustment of the state vector. If the optimum weighting matrix is selected, the value of the measurement noise covariance is: W = 5 7 -1. Formulae (l l) and (12) are used subsequently for evaluation of the areas of orbital parameter adjustment error on a geostationary orbit: cov,,(X) = ( B r W B ) '.
(12)
2. GENERAL RESULTS
2.1 Basis o f the analytical method--measurement noise Formulae (11) and (12) show that all covariance analyses (measurement noise, model noise, different types of ? sets) are reduced to a calculation of the matrix: (BrWB),
(BrWK).
The procedure to form these matrices will be described in the specific case of B r W B . For B r W K , the situation will not greatly differ. We have [with the notations of formula (5)]:
E(EE r ) = Z: measurement errors covariance, E(77 r ) = cov(7): covariance of lack of knowledge on parameters of ?. We therefore only have access to the covariance of the state vector evaluation error: cov(X) = E [ ( 2 - X ) ( 2 - X) r] = (B T W B ) - 1 (B rW57 W B ) (B r W B ) - ' + (BrWB)
~( B T W K ) coy(7)
x [ ( B r W B ) -' ( B r W K ) ] r.
(1 l)
(B WB) =
Bf
Bj;
j=l Wj is then the elementary weighting for each measurement. Wc can notice that such weighting is dimensioned m-2 (convention of Section 1.3). We separate the sum o n j ( j varying from 1 to n in all of the available measurements--all types--all stations) in sub-sums per type of measurement and per tracking station. By designating by (i) the generic index of a measuring instrument in a tracking station, by W(0,
38
O. ZARROUATI
the weighting--assumed to be constant--of all the measurements performed by this instrument and by N(,) the number of measurements performed by this instrument, one obtains:
Hypothesis 2: The duration of the arc is an integral number of days. Some mean values can be deduced directly from these hypotheses:
(BrWB) -__ ~ B(iT ) W(i ) B(,,), (i)
(UMUM)=I1;2
[0
N(i)
B(,r) W,) B(1)= W(,)~ BfBj.
(13)
j=l
(UMV~)=
1/20J =(VMV~')
_1/2
,
(UM)=(VM)=0.
It is clear in formula (5) that the expression of To calculate the other mean values, more data are necessary:
N(i) j=l
will result in the appearance of partial sums in the form: N(i)
N(i)
Z (4-2o)j, 2 (4-2o)~, j=l
j=l
Z u~jv~,~,
E uM~v;,~,
i
)
EvMjv~,~, Eu~,, EvM~, J
/
J
sums which can always be expressed through a method of writing, using mean values:
• The duration of the adjustment arc is 2M--that is to say that the spacecraft (at the rate of 2x radians daily) performs, on its orbit, a longitude excursion of 2M during the measurement period. • The orbital parameters are adjusted at a given date (longitude 2 of the spacecraft is variable in time). This date is situated in the middle of the adjustment arc in such a way that the spacecraft, during a period of measurement, describes the longitude section: [40 - 2n,:, 20 + 2M,':].
N(i)
(Xj)
=
We may then calculate:
N(,)X.
j=l
To calculate mean values, we must, of necessity, know the date of each measurement. In the majority of cases, measurement plans (station programming) require an almost uniform distribution of these dates on the adjustment arc.
(2
(2 - 4 0 ) = 0 ,
-
20)2 = 2 2 / 1 2
Hence, the expression of elementary symmetrical 6 × 6 matrices (remember that for elements e, h, 2, the state vector is normalized by a and deviations ale, a6h, a62 are considered):
+ l 2 2M 2
0
l~/2+21~
0
0
l~/2+2l~
N(i)
(14)
2 Bj B f = N(i)
-1~6/2 ll/2
66
j=l
0
Ii 13/2
~ l3
o
ll/2
1112
0
0
0
0
We may then dispense with detailed knowledge of the dates of each measurement by means of two hypotheses of distribution:
Hypothesis 1: For each station and for each type of measurement, the dates on which the measurements are performed have an uniform distribution on the adjustment arc. This hypothesis is all the more reaslistic as the orbit is geostationary, which means that visibility of the spacecraft from each station is continuous.
Here, l~, /2, 13, are those of the instrument (i). The global matrix (BTWB) is a weighted sum of similar elementary matrices, therefore: • The zeros present in each elementary matrix will normally be found in the global matrix. • By reordering the status vector, we make B r W B diagonal by blocks one block concerns the parameters
(6a, a62), the other block concerns the parameters
(a6 e, a5 h).
Evaluation of the performance of a radiolocation system • In inverting the B r W B matrix, the two systems (6a, a62) and (a6e, a6h) may be treated independently--the covariances of these two systems are completely decorrelated. The covariance of the system (Sa, a62) is obtained immediately by inversion of the symmetrical 2 x 2 matrix. The covariance of the system (a6e, afh) is obtained by inversion of a 4 x 4 matrix. Inversion is not immediate a priori but the particular form of the matrix to be inversed makes the operation analytically possible. The block concerning (a6e, a6h) of ( B r W B ) still has the following structure:
39
C
mes'~ the partial derivatives of the measure--~--7 } = ments performed by the instrument (i), l 1,12, l3 = those of the instrument (i),
We can see that the calculation of B r W K only causes intervention of the following mean quanties: (d mes'],
VL(2 - 20) dines ];
~3rues \
//V d mes'~
)' t
)'
(17.1)
(17.2)
This situation requires some discussion: (B r WB)/(a6 e, a6 h) =
We then check:
i°z l X
-W
-W
Y
Z
0
[o0 z o]
( B r W B ) - l/(a6e, a6h)
-W
1 A
Y
W -Z
W X 0
A = X Y - Z 2 -- W 2.
(15)
2.2 Model noise Model noises evaluation requires the formation of matrices of the B r W K type. We use the technique defined for the calculation of B r WB. Matrices B r W K are expressed [with formula (5) notations] as follows: BTWK=
y=l
B7
N(i)
B(r) W(i) K(i) = W(i) E B f Kj, j~l
L
/=I
(0mes'~
3/2 [ - Y
• constant (projection on 1--referred to as constant bias), • secular [projection on (2 -20)--referred to as drift], • periodical at the day's period (projection on UM, VM--referred to as medium-term bias),
Furthermore, the diagonal structure by blocks of ( B r W B ) -I noticed in (2.1) and formula (16) allow
separate analysis of the effects of model errors on the systems (6a, a52) and (a6e, a6h). Indeed, we notice in formula (11) that ( B r W B ) -~ ( B r W K ) is the basic matrix allowing model noises to be evaluated. One may split this matrix into two sub-matrices concerning the two blocks (6a, a62) and (a6e, a6h):
T W(i) K(i), B r W K __ ~, BU) (~)
l,
Consequently, only the components which are:
of the function (t~mes/0v) may impair orbit restitution.
where the sum on index j may be split up into sub-sums.
N(i)
• The quantities (17) to be evaluated are, in fact, scalar products (in functional space) of the function (0mes/aT) by the functions 1, (2 - 20), UM and VM. • In functional space as well as in a finite dimension space, the scalar product is maximum if the functions are proportional and nil if the functions are orthogonal. • Therefore, it is the projection of the function dmes/d~, (measurement error due to unknown parameters) on the functions 1, (2, 20)UM, VM which results in errors in the restitution of orbital parameters.
dmes(2 - 20)
( B r W B ) -1 B r W K / ( f a , a62) = ~ N(oW~o(BrWB)-I/(6a, a62) (~)
( U MOmes'~ Omes Oy j - 212 ( v u - - ~ - - y ) = N(i)
-1 {v
Omes'~
l ( mes
l~0( ~mes'] \-UC/
(16)
40
O. ZARROUATI
(BrWB) 1BVWK/(afe, a6h)
In addition, data (19.1) and (19.2) provide, with Appendix (1), a calculation of the components l~,/2,/3 of vectors I(i) relative to each type of measurement. For the distance (note, for further use, the definition of the generic indices), we obtain:
= ~ N(,)W(i)(BTWB)-I/(a6e, afh) (o _l]i) /
8mes\
(i) /
r:0l
Omes~ i~~ =
× - l ~ i~(VM ~mes'~8?, ,]
0.0 , ? 0 . 1 lJ
Fo.ol
and for the elevation and bearing:
(18) 1~s'= 10.34[
Formula (11) enables us to observe immediately that: • the quantities (17.1) only affect the "semi-major axis--longitude" system, • the quantities (17.2) only affect the "eccentricity inclination" system. From which, we can draw the following qualitative conclusions: • the constant or secular measurement biases disturb restitution of the semi-major axis and longitude of the spacecraft ONLY, • medium-term biases (periodic at the day's period), only disturb restitution of the eccentricity and inclination vectors. Finally, orbit restitution is a filter activated by functions (Smes/07) and whose outputs are the errors in the computed orbital parameters. It appears that only constant, linear and periodical at the day's period signals make the filter in resonance. In what follows, we shall therefore consider only these signals.
Ie ' =
[0.93J
0.94 •
(19.4)
L-0.35J
Finally, evaluation (19.3) concerning the elevation and bearing measurement accuracy can be expressed in metres [formulation (3)]. The elevation and bearing angular measurements are equivalent to range measurements about axes l(s) and l(g) with the following respective accuracies: c~ ~ % ~ 7500 m!
(19.5)
It can be seen that if the orbit restitution algorithm uses a weighting "close to" optimum weighting (which is the case within 1 or 2 orders of magnitude):
angular measurements are very under-weighted with respect to the range measurements:
Wd~210-Im 2, W,~ Wg~210 8m-2
3. CASE OF THE SINGLE STATION It is considered here that the geostationary spacecraft, at the stationing longitude 2p, is tracked by a single station with the geographical coordinates
(L, 4,3. The station performs range, elevation and azimuth measurements on the spacecraft.
3.1 Some orders of magnitude The numerical data given hereafter relate to the TDF1 satellite: 2p = -- 19°,
(19.1)
localized by the Aussaguel station
S P
2~= 1.5 °, qS~= 43.43.
For range, elevation and bearing measurements, the short-term noise standard-deviations are evaluated, respectively, at: ad = 2.5 m, a, = ~rg = 2 10 4 rad.
~s'
(19.2)
Xp
Stationing point longitude Zenithal unit vector at point S
ar
Earth's radius
a
Geosynchronous semi-major axis
N
(19.3)
~
Tracking station Stationing point Station coordinates
Fig. I. Calculation of I(i).
41
Evaluation of the performance of a radiolocation system There are almost 7 orders of magnitude of deviation between the weightings. 3.2 Measurement noises The measurement noise covariance matrix is evaluated using (13) and (14) by summing, in BrWB, the contributions of the 3 types of measurement. The observations made in (3.1) concerning the weightings of the various types of measurements would allow the determining character of the range measurement contribution to be declared beforehand. Angular measurements only seem to intervene as corrective terms. 3.3 Semi-major axis--longitude system: measurement
noises The BrWB block relative to this system must be calculated. The approximation is:
BrWB /(fa, a~52)
• It is Kepler's 3rd law which provides observability of the (6a, a62) system using range measurements performed by a single station. Restitution of this system relies on geometrical and dynamic considerations. • If only a geometrical resolution is desired (2 g ---0; the longitude is not given the time to drift), the angular measurements become necessary for BrWB to remain reversible. From (20) and (1 !), we deduce the matrix (BrWB) -I ( B r W E W B ) (BrWB) -b for the system (6a, a~52), hence the measurement noise standard deviations (here, weighting is not supposed to be optimum):
1
a~"-
4aa
(i)
=NdWd F(/la)2+ 3(12a)2 (d~d2l ) 2 l l
(approximation valid for an arc of a few days)
(12) l~l~ d2
In chapter (2.2) we have seen that only the constant or linear measurement biases can distort restitution of the system (6a, a62). We shall deal here only with constant range or angular measurement biases. We therefore define a set of parameters which are not well known:
Hence the system's measurement noise covariance matrix:
(BrWB) - L/(6a, a62) 1
= (bd, bs, bg)
a ~s
,
(22.2)
(t,~) ~ + 3(lg) ~
Nd: a number of ranging measurements carried out. We can see that this matrix only depends (3.2) on the performance of the ranging system. We notice that: • The system (6a, a62) is best restituted when 2M is higher (i.e. when the arc is longer); if we take 2M=0, the matrix BrWB here carried out would not be reversible. The angular measurement would then be determining. • Physically, the filter operations can be schematized as follows: - - a mean drift of the spacecraft/station range measurements can be observed; --this drift is geometrically interpreted as a drift of the spacecraft's longitude; --the semi-major axis is inferred from the longitude drift using Kepler's 3rd law; --the longitude is geometrically inferred from the semi-major axis. D
(22.1)
where the biases are represented on the 3 types of measurements. The poor knowledge of these biases is characterized by the following matrix:
cov(v) =
A.A. 151i
4lldltrd
3x/3-ffd2M(12d)2 (21)
3.4 Semi-major axis--longitude system: model noise
= ~ B(ri) W¢~ B¢i)/(~a, a62) ~ BTd WdBd/(fa, a62)
-t~/g
1
3x~ ,~Mll~I atr~a~
with the following realistic values of the standard deviations (corresponding to the Aussaguel station): trhd~ 20 m, ab~ ~ ab~ ~ 7000 m. We have seen in section 3 that angular measurements do not participate in the determination of the system (6a, a62) but the biases on these measurements are so considerable that they may distort restitution. Definition (22) allows the partial derivatives to be calculated: (1, 0, 0) if mes is a distance,
~mes _ 67
f
(0, 1,0) if mes is an elevation, (0, 0, 1) if mes is a bearing.
Manipulating (20), (18) and (11) then gives the covariance increment of the orbital parameters due to the measurement biases.
42
O. ZARROUATI Acov(fia, a62) =
~
Formula (15) shows that in order to reverse the 2 matrix, the following quantity must be calculated:
2 r ~ b(i) X(i) X ( i ) ,
(i) = d,s,g
1
A = X Y - Z 2 - W 2.
N(o Wo~
We can immediately see that when ignoring the angular measurements, A is zero. The distance measurements only do not provide inversion of the eccentricity-inclination system. It is (and only) when reversing this system that the angular measurements are useful. Inversion of the eccentricity-inclination system can
3(l:a)4(_~y Nd Wd F I d [1 d l(i) __ I d l(i)'~
7
i_-1,-1-2 - - 2 , 1 , - 3(19: l~')
]
Hence the standard deviations [to be summed quadratically between themselves and with those given in (21)]
ff 6a
aft6;.
0
a bd/I l~1
Elevation bias
abs N~ W~ II~l 3(2M/4)2 Nd Wa Ildl3
~rh~ N~ W~ Ila~l~l 3(ZM/4)2Nd Wd Ildl4
Bearing bias
abg Ng W e I1~1 3(2M/4)2 Na Wd I/zdl3
abe Ng W e II~ IS3[ 3(2M/4)2Nd Wd 1/2dl4
Range bias
We notice that: • the range measurement biases have zero or negligible impact on the restitution of (6a, a62); • the angular measurement biases have non negligible effects on the restitution of the system. As the elementary effects are directly proportional to the weightings of the angular measurements W~ and Wg, it will often be advisable to underweight these measurements with respect to the optimum weighting.
(23)
thus be approximated through:
( B r W B ) 1~(ale, afh)
1
A
io oZ] Y
W
W -Z
X 0
3.5 Eccentricity-inclination system: measurement noise The B r W B block relative to this system must be ~,calculated. Like in the case of the "semi-major axis-longitude" system, the current orders of magnitude seem to make it possible to have an approximation:
l 0Zi 0
X
-W Z
-W
A: calculated with the contribution of all measurements (distance and angular), i.e.:
z N(O Wo~~"'//z L(i) ~ /L(i) [~ -
N(0
(0
=
2
]
=2
Y 0
F(lla)z + 2(12d)2] Na Wa L 2 j
l(i)l(i)] 2
The calculations developed in Appendix (2) allow the expression of A to be simplified to:
A = NaWaNgWgN, W~
[
Z = [l~l~] N a Wa W = [la~l~/2] Na Wd
N(i) W~,)
l(i) 1(072 '3 /
W~ '1
~
X
1(027
1(i)27 I-
A=
BrWB/(a6e, a3h) ~ B~" WaBd/(a6e, afh)
=
X, Y, Z, W: correspond to definitions (24) (contribution of the sole distance measurements).
(24)
l: ~v
IU]
x Zk L-4N,~.,W, - - W(k)d . (k) (*) F Nck)
Evaluation of the performance of a radiolocation system Hence, taking the current orders of magnitude into account:
A ~ NdWaNgWsN~Ws
standard deviation, but on the weighting they are granted. The physical performance of the angular measurement instruments in short term noise does not affect orbit calculation accuracy.
x [(l~)2/4Ng Wg + (l~)2/4N~ W~],
3.6 Eccentricity-inclination system: model noise
this formula gives, under a realistic hypothesis of elevation and bearing measurements achieved in equal
number and identically weighted: A
43
N~ W~Ns, W~, 4
N~W.= N, We= N,,Ws,.
We have seen in Section 2.2 that only medium-term biases (periodical at day's period) may distort restitution of the system (ale, afh). We shall therefore only be concerned with these measurement error signatures. Hence three ill-known sets of parameters are defined:
It is then possible to calculate the block relative to the eccentricity-inclination system of the matrix
( B T W B ) - I ( B T W E W B ) (BTWB) 1. The value of the system adjustment covariance matrix (measurement noise only) is therefore (here again, the optimum weighting hypothesis is not used):
y(i) = (6~°, 6~°) = 8(0,
(26.1)
so that the error on the measurement of type (i) can be written: 3 mes(i) = 8(i). U r ,
(26.2)
(lg):/2 0 4(tr2a Wa)
cov(afe, afih)
N~gG
(13a)2/2
d d -1213
lfl~/2 (1~)2/2 + 2(12a)2
-l¢l~/2 -l~l~
0 (t¢):/2+2(l~) ~
Hence, finally the standard deviations: =
aaex
la U 2
aaey= 3 4 N s g W s g ~ W d
aahx = aa.y ~ l~
~
~
(25)
We can clearly see that it is the angular measurements which, by allowing system inversion, provide eccentricity and inclination. The parameters can be observed all the better as weighting granted to the angular measurements is higher. We notice that, in case of non optimum weighting: W~g# 1/a~,
#: 1/o'g2.
The eccentricity and inclination calculation performance does not depend on the angular measurements
Thus periodical measurement errors at day's period have been well defined. The partial derivatives of each type of measurements with respect to parameters of the corresponding y set are: 0mes(i) = U~, t3T(i)
(26.3)
so that calculation of mean quantities already encountered (UMU~) and (UMVrM) and application of formula (18) gives the orbital parameters covariance increment due to the periodical error: Acov (a6 e, a6 h)
=
~
X(i)[a~ f
9~,f]X(i) r
(i) = d,s,g
1d[Id l(i)
l~(lal I(3° -- I~0 la3)
• 3k~3 ~2
l~.~ 1~o_ l~O1~) 3 k ~ l ~3
21~(I~ I(3° -- I~ I(2°) X(i) = N(,~ W(O
I d I(i)3 -- ~2"3 )
d d (i) 412(12 13 - 13d 12(i)) "~+ ld(l d Ito l g IIO)~
tdttd "3~.'2 ~(o "1 -- lal l(2i))
-
f
l• d1[\ i) 3d I(i) I d l(i)$ ~1 - - ~ I L 3 1
)~.
+ 4ld(ld31~~)-- I d lt,))j
-
d d (i) 213(lz ll -- lal l(2i))
(27)
44
O. ZARROUATI
with ~r~,?and cr~ the respective standard deviations of parameters 61!) and 61.'L The results given by the analytical method have been indicated, taking into account only the prevailing terms (above formulae) as well as those given by the analytical method, taking all the terms into account (the most important corrective terms are due to angular measurement biases). Simplified analytical method
Complete analytical method
Numerical method
28 0.59 3.1 1.30
31 0.59 3.1 1.375
31 0.60 3.0 1.40
cr6~(m) ~%.,.= ~6,,~( × 105)
tY6hx=t76hy(lO 3d'~) O',s;(10 -3 d °)
which can be approximated by: l~
= R r S I,
(A1.7.1)
av
1~ = S. cos 4~,.sin (2s - 2e),
(AI.7.2)
1~ = S. sin 4)~.
(AI.7.3)
It is then important to notice that for a same station and for the measurements considered above, we always have: I(i) A I ( j ) = _+l(k).
(AI.8)
APPENDIX 2
Calculation of A APPENDIX I
A2.1 Preamble
Calculation of l(i) (Fig. 1)
Let us consider quantities of the type:
AI.I Range measurement Range measurement positions the spacecraft in a plane orthogonal to axis SP. We thus have: I d = SP/ISPI.
(ml.1)
This formula provides a rigorous calculation of coordinates lid, 12a, 1~. A 2nd order approximation in (Rr/a~) gives: l~ = 1,
(A1.2.2)
Rv l~ = - - - cos ~bs sin (2 s - £p),
(A1.2.2)
as
Rr
1~ = - - -
Q =~c, cia~b~-~c, cja, b~a,bi i, I
= ~
I,j
cz cj(a~2 b i2 aiajbibj )
= ~ ci cj(a, bj - aj bi)2.
(A2.1)
I>j
sin q~s,
(A1.2.3)
as
A2.2 Application
Al.2 Bearing measurement A constant-elevation bearing deviation displaces the spacecraft along an axis parallel to Z ^ SP.
The quantity to be calculated is of the " Q " type defined above. We thus have: A=
We thus have: F=(Z
by developing we have:
A S P ) / I Z A SPI
(A1.3)
x [1~') l~j) - l~J ) l~i)]:
Here again, this formula provides a rigorous calcug g lation of l~, 12, l 3. However, we propose the following approximation: l~ = 0,
(A1.4.1)
l~ = S. sin q~s,
(A1.4.2)
l~ = - S . cos ~ sin (2s - 2p),
(A1.4.3)
with S = [1 --cos 2 ~bscosZ(£s - 2e)] ,.'2 (A1.5)
+ y N~,)~,~N,j)~j) l >/
x r/if)/(J) ,.2 . 3
_
1~ I l~'q 2.
This is the case of a single tracking station performing measurements of the range, elevation and bearing types. Formula (AI.8) thus makes it possible to write:
A= ~
N(i)W~,)N(j)W(j){[l~k)]2/4+ [l~k)]2}
i>j k # (i,j)
AI.3 Elevation measurement The motions of the "constant elevation" or "constant bearing" antenna are orthogonal. We thus take: IS= ia A is,
1
~ i~>iN(1)W(1)N(j) W(j)
(A1.6)
= UdWdUsW~U~W.
×~
c [l~q = [li~q~ - - 4 - - . k ~4N(k) W(k) N(k)W(k)t
(A2.2)
0094-5767/87 $3.00+ 0.00 © 1987 Pergamon Journals Ltd
A SIMPLE ANALYTIC METHOD FOR THE EVALUATION OF THE PERFORMANCE OF A RADIOLOCATION SYSTEMt O. ZARROUATI Centre National d'Etudes Spatiales (CNES), 18 Av. E. Belin, 31400 Toulouse, France
(Received 11 February 1986)
Abstract--The problem is, knowing the distribution over the world of a set of radiolocation beacons, and the individual performance of any beacon, to evaluate the performance of the global system so obtained for the localization of a satellite on a particular orbit. Under some very realistic assumptions on the orbit, the author gives an analytical formulation of the covariance of the computation of the satellite position. The formulation is illustrated in one particular case: localization of a GEO with a ground network (distance and angular measurements), but may be used in much different situations. The results take into account the "white" measurement noise and the "colored" modelization noise. The method can be very helpful in mission analysis for the optimal dimensioning of a radiolocation system.
i. INTRODUCTION 1.1 Hypotheses The method presented here allows the size and form of the geostationary spacecraft ranging error volume to be evaluated. The ranging hypotheses and conditions are as follows: (i) Range and angular measurements (elevation, bearing, Cardan a n g l e s , . . . ) only. (ii) Unlimited number of Earth tracking stations. (i~i) Orbit calculation is performed (batch processing) using a set of measurements obtained over a given period of time (a few days: this period is called an "adjustment arc"). (iv) The spacecraft position and velocity (or its orbital parameters, which is equivalent) are adjusted at a given date (here: the middle of the measuring period) so that the resulting ephemeris "fit best" through the measurements of the period. The method used is that of the least squares fit. (v) The orbit is near geostationary. That is, the spacecraft and the tracking stations relative geometry does not vary during the adjustment arc. The spacecraft motion is Keplerian.
ex = e cos(co + ~): ey = e sin(o) + f~): components of the "eccentricity vector" e, which belongs to the equatorial plane and is directed towards the perigee. Its modulus is the eccentricity e. hx = i cos ~: hy = i sin ~: components of the "inclination vector" h, which belongs to the equatorial plane and is directed about the rotation axis of the orbital plane with respect to the equatorial plane. Its modulus is the orbit inclination i. 2 = M + co + f~: spacecraft mean longitude counted from the reference axis 07 (Vernal axis). During orbit calculation we know the nominal orbital parameters (a0, e0, ho, 2o) at date to but we do not know the increments (6a, fie, fib, 62) by which these parameters must be modified to "fit best through the measurements" in the least squares fit sense. Conversely, we know how to evaluate, at any date • according to Keplerian motion hypothesis (v), • using the definition of the adapted parameters, the impact of these orbital parameters increments upon the spacecraft positioning. It is a positioning deviation 6r given by:
fir = afh ^ U M - a b e - (VMV r ) are
1.2 Orbital-dynamic parameters of the problem
+ VM[ar~. -- 3 (2 -- 20)6a] + UM 6a
Parameters adapted to the geostationary orbit are used. These parameters are expressed, using conventional Keplerian parameters (a, e, i, co, ~, M), by:
(1)
with: U M = u n i t vector directed about the " E a r t h ' s center-mean spacecraft" axis, k = unit vector directed about the pole axis towards the north, V M = k ^ UM.
a: semi-major axis tPaper presented at the 36th Congress of the International Astronautical Federation, Stockholm, Sweden, 7-12 October 1985. 35
36
O. ZARROUATI
In the expression of fir, the following term appears: --3 VM(2 -- 20)~a; this term does not directly derive from the orbital parameters definition. This is a term induced by Kepler's 3rd law. The expression ( 2 - 20) that appears is the longitude excursion (in radians) followed by the geostationary spacecraft as from the adjustment date. A date when its longitude is 20. Equation (1) immediately gives the projection relations: (6r UM) = -U~t(afe) + 6a,
(2.1)
(fir VM) = -2V~t(afe) + a62 _ 3 (2 - 20)6a, (2.2) (fir k) = - V r ( a f h ) .
(2.3)
1.3 Measurement modelling Any range measurement between the spacecraft and a station and any angular measurement can be modelled by a form: rues = (6r.I)
(3)
where I is a unit vector depending on the type of measurement. • The range measurement is a projection of the deviation 6r on the "station-spacecraft" axis. • An angular measurement is equivalent to a projection of the deviation fir on an axis normal to the "station-spacecraft" axis. Formulation (3) thus allows any type of measurement (range or angular) to be assigned to a standard deviation a in distance units (in meters). This formulation is a rather general approximation which only relies on the assumption of the spacecraft/tracking station distance, a long distance with respect to a. But the geostationary orbit benefits from an additional property: For each tracking station and for each type of measurement, I is fixed within the Earth's reference frame. We may also write: F ) = l~i) U u + 1~i) ~2 v V M + l~i) k,
1]~, l~~, l~~)=constants,
l]m + l~i~2 + l~m = 1;
(4)
(i) is the generic index of a type of measurement on a given tracking station. Finally, each measuring instrument of each tracking station is defined by 4 parameters: I(i) 1(i) I(i) O-0)). ~1 , ~ 2 , ~ 3 ,
Appendix 1 shows how these parameters are computed for elevation, bearing and range measurements as a function of the station coordinates (2s, ~b3, of the spacecraft positioning longitude 2p (longitudes 2~ and 2p are calculated within the Earth's reference frame,
unlike the spacecraft longitude 2 on its orbit, calculated from the vernal axis). The partial derivatives of the measurements obtained in relation to the adjusted parameters (also called state vector parameters) are directly inferred from formulae (2) to (4): - 8
mesj c36a
8mesj a86e
Bj =
8mesj at36h
l, - ~ 6 (~ - ;.o)~
- Ii U us - 2/2 VMj (5) - 6vMj
8mesj a062
6
• ~. is the spacecraft longitude on its orbit at the date when the measurement: mesj was carried out. The spacecraft is then on the axis directed by Uuj. Its velocity is along V~j. • The measuring instrument which performed mesj is now characterized by values (ll, 12,/3). Some important orders of magnitude can be inferred from the calculations given in Appendix 1. • For the range measurement, we have 1 ~ 1 , ~12, l 3.
• For the angular measurement: 12 ~13 ~ l I n 0 ;
These figures result from the high value of the geosynchronous semi-major axis with respect to the Earth's equatorial radius. 1.4 Some results o f the least squares theory
Orbit modelling and measurement simulation depend on a certain number of unknown or uncertain parameters. Unknown parameters will be identified by the least squares fit method--this is the state vector (normalized by a): X = (6a, a6e, a f h , a62).
Errors will be made in the knowledge of uncertain parameters. As these parameters are not adjusted (one cannot adjust all of them or the observability problem will be impaired), errors made in their evaluation will inevitably impair identification of X. The set of uncertain parameters is noted 7: 7 = (calibration bias, station coordinates, e t c . . . ) One writes respectively B and K as the partial derivatives matrices for all the measurements carried out in relation to the X state vector parameters and to the uncertain parameters.
Evaluation of the performance of a radiolocation system
37
That is, by using, in addition, formula (5) notations: 8mes~
B=
8mes~
BI
8X .
87
= B y , K =
c3mesn
8X
mesj = B j X - Kj7 - Ej,
(7)
Ej= creation of a Gaussian white noise with a standard deviation a for measurement j. By incorporating all the measurements (mesj)j = 1, n in a single-column vector M: m = [mes, . . . . mesj . . . . mes,] r. The X state vector identification problem is reduced to the resolution of a redundant linear system: (8)
B Y = M.
If the weights (reliability) attributed to each measurement are grouped together within a n × n diagonal matrix W, the "best solution" of the system (8) in the least squares fit sense is known: (9)
~B T W M .
If we knew Ej and ? (which is not the case), we could directly evaluate, by means of (7), the error made by the least squares (9) estimator: - X=(BrWB)
I BTW(BX_
= ( B r W B ) -l B r W ( - E
K?_E)_
K.
(6)
87
One obtains a set of n measurements (mesj) j = 1, n which are reduced to position deviation 5r projections on axes I(i). The non-zero value of these measurements (existence of a deviation in position) may be due to the state vector (the real orbit is not the nominal orbit) or to a lack of knowledge of the parameters 7- Linearization enables each measurement to be written in the form:
= (BrWB)
.
c3mes,
o.
KI
X
- K?).
(10)
In fact, we only know:
This covariance is calculated under the hypothesis: E(E~ q = 0. It is the sum of two terms: • A so-called "Measurement noise" term which describes the effect of measurement error 57 on the adjustment of the state vector. • A so-called "Model noise" which describes the effect of lack of knowledge of y on the adjustment of the state vector. If the optimum weighting matrix is selected, the value of the measurement noise covariance is: W = 5 7 -1. Formulae (l l) and (12) are used subsequently for evaluation of the areas of orbital parameter adjustment error on a geostationary orbit: cov,,(X) = ( B r W B ) '.
(12)
2. GENERAL RESULTS
2.1 Basis o f the analytical method--measurement noise Formulae (11) and (12) show that all covariance analyses (measurement noise, model noise, different types of ? sets) are reduced to a calculation of the matrix: (BrWB),
(BrWK).
The procedure to form these matrices will be described in the specific case of B r W B . For B r W K , the situation will not greatly differ. We have [with the notations of formula (5)]:
E(EE r ) = Z: measurement errors covariance, E(77 r ) = cov(7): covariance of lack of knowledge on parameters of ?. We therefore only have access to the covariance of the state vector evaluation error: cov(X) = E [ ( 2 - X ) ( 2 - X) r] = (B T W B ) - 1 (B rW57 W B ) (B r W B ) - ' + (BrWB)
~( B T W K ) coy(7)
x [ ( B r W B ) -' ( B r W K ) ] r.
(1 l)
(B WB) =
Bf
Bj;
j=l Wj is then the elementary weighting for each measurement. Wc can notice that such weighting is dimensioned m-2 (convention of Section 1.3). We separate the sum o n j ( j varying from 1 to n in all of the available measurements--all types--all stations) in sub-sums per type of measurement and per tracking station. By designating by (i) the generic index of a measuring instrument in a tracking station, by W(0,
38
O. ZARROUATI
the weighting--assumed to be constant--of all the measurements performed by this instrument and by N(,) the number of measurements performed by this instrument, one obtains:
Hypothesis 2: The duration of the arc is an integral number of days. Some mean values can be deduced directly from these hypotheses:
(BrWB) -__ ~ B(iT ) W(i ) B(,,), (i)
(UMUM)=I1;2
[0
N(i)
B(,r) W,) B(1)= W(,)~ BfBj.
(13)
j=l
(UMV~)=
1/20J =(VMV~')
_1/2
,
(UM)=(VM)=0.
It is clear in formula (5) that the expression of To calculate the other mean values, more data are necessary:
N(i) j=l
will result in the appearance of partial sums in the form: N(i)
N(i)
Z (4-2o)j, 2 (4-2o)~, j=l
j=l
Z u~jv~,~,
E uM~v;,~,
i
)
EvMjv~,~, Eu~,, EvM~, J
/
J
sums which can always be expressed through a method of writing, using mean values:
• The duration of the adjustment arc is 2M--that is to say that the spacecraft (at the rate of 2x radians daily) performs, on its orbit, a longitude excursion of 2M during the measurement period. • The orbital parameters are adjusted at a given date (longitude 2 of the spacecraft is variable in time). This date is situated in the middle of the adjustment arc in such a way that the spacecraft, during a period of measurement, describes the longitude section: [40 - 2n,:, 20 + 2M,':].
N(i)
(Xj)
=
We may then calculate:
N(,)X.
j=l
To calculate mean values, we must, of necessity, know the date of each measurement. In the majority of cases, measurement plans (station programming) require an almost uniform distribution of these dates on the adjustment arc.
(2
(2 - 4 0 ) = 0 ,
-
20)2 = 2 2 / 1 2
Hence, the expression of elementary symmetrical 6 × 6 matrices (remember that for elements e, h, 2, the state vector is normalized by a and deviations ale, a6h, a62 are considered):
+ l 2 2M 2
0
l~/2+21~
0
0
l~/2+2l~
N(i)
(14)
2 Bj B f = N(i)
-1~6/2 ll/2
66
j=l
0
Ii 13/2
~ l3
o
ll/2
1112
0
0
0
0
We may then dispense with detailed knowledge of the dates of each measurement by means of two hypotheses of distribution:
Hypothesis 1: For each station and for each type of measurement, the dates on which the measurements are performed have an uniform distribution on the adjustment arc. This hypothesis is all the more reaslistic as the orbit is geostationary, which means that visibility of the spacecraft from each station is continuous.
Here, l~, /2, 13, are those of the instrument (i). The global matrix (BTWB) is a weighted sum of similar elementary matrices, therefore: • The zeros present in each elementary matrix will normally be found in the global matrix. • By reordering the status vector, we make B r W B diagonal by blocks one block concerns the parameters
(6a, a62), the other block concerns the parameters
(a6 e, a5 h).
Evaluation of the performance of a radiolocation system • In inverting the B r W B matrix, the two systems (6a, a62) and (a6e, a6h) may be treated independently--the covariances of these two systems are completely decorrelated. The covariance of the system (Sa, a62) is obtained immediately by inversion of the symmetrical 2 x 2 matrix. The covariance of the system (a6e, afh) is obtained by inversion of a 4 x 4 matrix. Inversion is not immediate a priori but the particular form of the matrix to be inversed makes the operation analytically possible. The block concerning (a6e, a6h) of ( B r W B ) still has the following structure:
39
C
mes'~ the partial derivatives of the measure--~--7 } = ments performed by the instrument (i), l 1,12, l3 = those of the instrument (i),
We can see that the calculation of B r W K only causes intervention of the following mean quanties: (d mes'],
VL(2 - 20) dines ];
~3rues \
//V d mes'~
)' t
)'
(17.1)
(17.2)
This situation requires some discussion: (B r WB)/(a6 e, a6 h) =
We then check:
i°z l X
-W
-W
Y
Z
0
[o0 z o]
( B r W B ) - l/(a6e, a6h)
-W
1 A
Y
W -Z
W X 0
A = X Y - Z 2 -- W 2.
(15)
2.2 Model noise Model noises evaluation requires the formation of matrices of the B r W K type. We use the technique defined for the calculation of B r WB. Matrices B r W K are expressed [with formula (5) notations] as follows: BTWK=
y=l
B7
N(i)
B(r) W(i) K(i) = W(i) E B f Kj, j~l
L
/=I
(0mes'~
3/2 [ - Y
• constant (projection on 1--referred to as constant bias), • secular [projection on (2 -20)--referred to as drift], • periodical at the day's period (projection on UM, VM--referred to as medium-term bias),
Furthermore, the diagonal structure by blocks of ( B r W B ) -I noticed in (2.1) and formula (16) allow
separate analysis of the effects of model errors on the systems (6a, a52) and (a6e, a6h). Indeed, we notice in formula (11) that ( B r W B ) -~ ( B r W K ) is the basic matrix allowing model noises to be evaluated. One may split this matrix into two sub-matrices concerning the two blocks (6a, a62) and (a6e, a6h):
T W(i) K(i), B r W K __ ~, BU) (~)
l,
Consequently, only the components which are:
of the function (t~mes/0v) may impair orbit restitution.
where the sum on index j may be split up into sub-sums.
N(i)
• The quantities (17) to be evaluated are, in fact, scalar products (in functional space) of the function (0mes/aT) by the functions 1, (2 - 20), UM and VM. • In functional space as well as in a finite dimension space, the scalar product is maximum if the functions are proportional and nil if the functions are orthogonal. • Therefore, it is the projection of the function dmes/d~, (measurement error due to unknown parameters) on the functions 1, (2, 20)UM, VM which results in errors in the restitution of orbital parameters.
dmes(2 - 20)
( B r W B ) -1 B r W K / ( f a , a62) = ~ N(oW~o(BrWB)-I/(6a, a62) (~)
( U MOmes'~ Omes Oy j - 212 ( v u - - ~ - - y ) = N(i)
-1 {v
Omes'~
l ( mes
l~0( ~mes'] \-UC/
(16)
40
O. ZARROUATI
(BrWB) 1BVWK/(afe, a6h)
In addition, data (19.1) and (19.2) provide, with Appendix (1), a calculation of the components l~,/2,/3 of vectors I(i) relative to each type of measurement. For the distance (note, for further use, the definition of the generic indices), we obtain:
= ~ N(,)W(i)(BTWB)-I/(a6e, afh) (o _l]i) /
8mes\
(i) /
r:0l
Omes~ i~~ =
× - l ~ i~(VM ~mes'~8?, ,]
0.0 , ? 0 . 1 lJ
Fo.ol
and for the elevation and bearing:
(18) 1~s'= 10.34[
Formula (11) enables us to observe immediately that: • the quantities (17.1) only affect the "semi-major axis--longitude" system, • the quantities (17.2) only affect the "eccentricity inclination" system. From which, we can draw the following qualitative conclusions: • the constant or secular measurement biases disturb restitution of the semi-major axis and longitude of the spacecraft ONLY, • medium-term biases (periodic at the day's period), only disturb restitution of the eccentricity and inclination vectors. Finally, orbit restitution is a filter activated by functions (Smes/07) and whose outputs are the errors in the computed orbital parameters. It appears that only constant, linear and periodical at the day's period signals make the filter in resonance. In what follows, we shall therefore consider only these signals.
Ie ' =
[0.93J
0.94 •
(19.4)
L-0.35J
Finally, evaluation (19.3) concerning the elevation and bearing measurement accuracy can be expressed in metres [formulation (3)]. The elevation and bearing angular measurements are equivalent to range measurements about axes l(s) and l(g) with the following respective accuracies: c~ ~ % ~ 7500 m!
(19.5)
It can be seen that if the orbit restitution algorithm uses a weighting "close to" optimum weighting (which is the case within 1 or 2 orders of magnitude):
angular measurements are very under-weighted with respect to the range measurements:
Wd~210-Im 2, W,~ Wg~210 8m-2
3. CASE OF THE SINGLE STATION It is considered here that the geostationary spacecraft, at the stationing longitude 2p, is tracked by a single station with the geographical coordinates
(L, 4,3. The station performs range, elevation and azimuth measurements on the spacecraft.
3.1 Some orders of magnitude The numerical data given hereafter relate to the TDF1 satellite: 2p = -- 19°,
(19.1)
localized by the Aussaguel station
S P
2~= 1.5 °, qS~= 43.43.
For range, elevation and bearing measurements, the short-term noise standard-deviations are evaluated, respectively, at: ad = 2.5 m, a, = ~rg = 2 10 4 rad.
~s'
(19.2)
Xp
Stationing point longitude Zenithal unit vector at point S
ar
Earth's radius
a
Geosynchronous semi-major axis
N
(19.3)
~
Tracking station Stationing point Station coordinates
Fig. I. Calculation of I(i).
41
Evaluation of the performance of a radiolocation system There are almost 7 orders of magnitude of deviation between the weightings. 3.2 Measurement noises The measurement noise covariance matrix is evaluated using (13) and (14) by summing, in BrWB, the contributions of the 3 types of measurement. The observations made in (3.1) concerning the weightings of the various types of measurements would allow the determining character of the range measurement contribution to be declared beforehand. Angular measurements only seem to intervene as corrective terms. 3.3 Semi-major axis--longitude system: measurement
noises The BrWB block relative to this system must be calculated. The approximation is:
BrWB /(fa, a~52)
• It is Kepler's 3rd law which provides observability of the (6a, a62) system using range measurements performed by a single station. Restitution of this system relies on geometrical and dynamic considerations. • If only a geometrical resolution is desired (2 g ---0; the longitude is not given the time to drift), the angular measurements become necessary for BrWB to remain reversible. From (20) and (1 !), we deduce the matrix (BrWB) -I ( B r W E W B ) (BrWB) -b for the system (6a, a~52), hence the measurement noise standard deviations (here, weighting is not supposed to be optimum):
1
a~"-
4aa
(i)
=NdWd F(/la)2+ 3(12a)2 (d~d2l ) 2 l l
(approximation valid for an arc of a few days)
(12) l~l~ d2
In chapter (2.2) we have seen that only the constant or linear measurement biases can distort restitution of the system (6a, a62). We shall deal here only with constant range or angular measurement biases. We therefore define a set of parameters which are not well known:
Hence the system's measurement noise covariance matrix:
(BrWB) - L/(6a, a62) 1
= (bd, bs, bg)
a ~s
,
(22.2)
(t,~) ~ + 3(lg) ~
Nd: a number of ranging measurements carried out. We can see that this matrix only depends (3.2) on the performance of the ranging system. We notice that: • The system (6a, a62) is best restituted when 2M is higher (i.e. when the arc is longer); if we take 2M=0, the matrix BrWB here carried out would not be reversible. The angular measurement would then be determining. • Physically, the filter operations can be schematized as follows: - - a mean drift of the spacecraft/station range measurements can be observed; --this drift is geometrically interpreted as a drift of the spacecraft's longitude; --the semi-major axis is inferred from the longitude drift using Kepler's 3rd law; --the longitude is geometrically inferred from the semi-major axis. D
(22.1)
where the biases are represented on the 3 types of measurements. The poor knowledge of these biases is characterized by the following matrix:
cov(v) =
A.A. 151i
4lldltrd
3x/3-ffd2M(12d)2 (21)
3.4 Semi-major axis--longitude system: model noise
= ~ B(ri) W¢~ B¢i)/(~a, a62) ~ BTd WdBd/(fa, a62)
-t~/g
1
3x~ ,~Mll~I atr~a~
with the following realistic values of the standard deviations (corresponding to the Aussaguel station): trhd~ 20 m, ab~ ~ ab~ ~ 7000 m. We have seen in section 3 that angular measurements do not participate in the determination of the system (6a, a62) but the biases on these measurements are so considerable that they may distort restitution. Definition (22) allows the partial derivatives to be calculated: (1, 0, 0) if mes is a distance,
~mes _ 67
f
(0, 1,0) if mes is an elevation, (0, 0, 1) if mes is a bearing.
Manipulating (20), (18) and (11) then gives the covariance increment of the orbital parameters due to the measurement biases.
42
O. ZARROUATI Acov(fia, a62) =
~
Formula (15) shows that in order to reverse the 2 matrix, the following quantity must be calculated:
2 r ~ b(i) X(i) X ( i ) ,
(i) = d,s,g
1
A = X Y - Z 2 - W 2.
N(o Wo~
We can immediately see that when ignoring the angular measurements, A is zero. The distance measurements only do not provide inversion of the eccentricity-inclination system. It is (and only) when reversing this system that the angular measurements are useful. Inversion of the eccentricity-inclination system can
3(l:a)4(_~y Nd Wd F I d [1 d l(i) __ I d l(i)'~
7
i_-1,-1-2 - - 2 , 1 , - 3(19: l~')
]
Hence the standard deviations [to be summed quadratically between themselves and with those given in (21)]
ff 6a
aft6;.
0
a bd/I l~1
Elevation bias
abs N~ W~ II~l 3(2M/4)2 Nd Wa Ildl3
~rh~ N~ W~ Ila~l~l 3(ZM/4)2Nd Wd Ildl4
Bearing bias
abg Ng W e I1~1 3(2M/4)2 Na Wd I/zdl3
abe Ng W e II~ IS3[ 3(2M/4)2Nd Wd 1/2dl4
Range bias
We notice that: • the range measurement biases have zero or negligible impact on the restitution of (6a, a62); • the angular measurement biases have non negligible effects on the restitution of the system. As the elementary effects are directly proportional to the weightings of the angular measurements W~ and Wg, it will often be advisable to underweight these measurements with respect to the optimum weighting.
(23)
thus be approximated through:
( B r W B ) 1~(ale, afh)
1
A
io oZ] Y
W
W -Z
X 0
3.5 Eccentricity-inclination system: measurement noise The B r W B block relative to this system must be ~,calculated. Like in the case of the "semi-major axis-longitude" system, the current orders of magnitude seem to make it possible to have an approximation:
l 0Zi 0
X
-W Z
-W
A: calculated with the contribution of all measurements (distance and angular), i.e.:
z N(O Wo~~"'//z L(i) ~ /L(i) [~ -
N(0
(0
=
2
]
=2
Y 0
F(lla)z + 2(12d)2] Na Wa L 2 j
l(i)l(i)] 2
The calculations developed in Appendix (2) allow the expression of A to be simplified to:
A = NaWaNgWgN, W~
[
Z = [l~l~] N a Wa W = [la~l~/2] Na Wd
N(i) W~,)
l(i) 1(072 '3 /
W~ '1
~
X
1(027
1(i)27 I-
A=
BrWB/(a6e, a3h) ~ B~" WaBd/(a6e, afh)
=
X, Y, Z, W: correspond to definitions (24) (contribution of the sole distance measurements).
(24)
l: ~v
IU]
x Zk L-4N,~.,W, - - W(k)d . (k) (*) F Nck)
Evaluation of the performance of a radiolocation system Hence, taking the current orders of magnitude into account:
A ~ NdWaNgWsN~Ws
standard deviation, but on the weighting they are granted. The physical performance of the angular measurement instruments in short term noise does not affect orbit calculation accuracy.
x [(l~)2/4Ng Wg + (l~)2/4N~ W~],
3.6 Eccentricity-inclination system: model noise
this formula gives, under a realistic hypothesis of elevation and bearing measurements achieved in equal
number and identically weighted: A
43
N~ W~Ns, W~, 4
N~W.= N, We= N,,Ws,.
We have seen in Section 2.2 that only medium-term biases (periodical at day's period) may distort restitution of the system (ale, afh). We shall therefore only be concerned with these measurement error signatures. Hence three ill-known sets of parameters are defined:
It is then possible to calculate the block relative to the eccentricity-inclination system of the matrix
( B T W B ) - I ( B T W E W B ) (BTWB) 1. The value of the system adjustment covariance matrix (measurement noise only) is therefore (here again, the optimum weighting hypothesis is not used):
y(i) = (6~°, 6~°) = 8(0,
(26.1)
so that the error on the measurement of type (i) can be written: 3 mes(i) = 8(i). U r ,
(26.2)
(lg):/2 0 4(tr2a Wa)
cov(afe, afih)
N~gG
(13a)2/2
d d -1213
lfl~/2 (1~)2/2 + 2(12a)2
-l¢l~/2 -l~l~
0 (t¢):/2+2(l~) ~
Hence, finally the standard deviations: =
aaex
la U 2
aaey= 3 4 N s g W s g ~ W d
aahx = aa.y ~ l~
~
~
(25)
We can clearly see that it is the angular measurements which, by allowing system inversion, provide eccentricity and inclination. The parameters can be observed all the better as weighting granted to the angular measurements is higher. We notice that, in case of non optimum weighting: W~g# 1/a~,
#: 1/o'g2.
The eccentricity and inclination calculation performance does not depend on the angular measurements
Thus periodical measurement errors at day's period have been well defined. The partial derivatives of each type of measurements with respect to parameters of the corresponding y set are: 0mes(i) = U~, t3T(i)
(26.3)
so that calculation of mean quantities already encountered (UMU~) and (UMVrM) and application of formula (18) gives the orbital parameters covariance increment due to the periodical error: Acov (a6 e, a6 h)
=
~
X(i)[a~ f
9~,f]X(i) r
(i) = d,s,g
1d[Id l(i)
l~(lal I(3° -- I~0 la3)
• 3k~3 ~2
l~.~ 1~o_ l~O1~) 3 k ~ l ~3
21~(I~ I(3° -- I~ I(2°) X(i) = N(,~ W(O
I d I(i)3 -- ~2"3 )
d d (i) 412(12 13 - 13d 12(i)) "~+ ld(l d Ito l g IIO)~
tdttd "3~.'2 ~(o "1 -- lal l(2i))
-
f
l• d1[\ i) 3d I(i) I d l(i)$ ~1 - - ~ I L 3 1
)~.
+ 4ld(ld31~~)-- I d lt,))j
-
d d (i) 213(lz ll -- lal l(2i))
(27)
44
O. ZARROUATI
with ~r~,?and cr~ the respective standard deviations of parameters 61!) and 61.'L The results given by the analytical method have been indicated, taking into account only the prevailing terms (above formulae) as well as those given by the analytical method, taking all the terms into account (the most important corrective terms are due to angular measurement biases). Simplified analytical method
Complete analytical method
Numerical method
28 0.59 3.1 1.30
31 0.59 3.1 1.375
31 0.60 3.0 1.40
cr6~(m) ~%.,.= ~6,,~( × 105)
tY6hx=t76hy(lO 3d'~) O',s;(10 -3 d °)
which can be approximated by: l~
= R r S I,
(A1.7.1)
av
1~ = S. cos 4~,.sin (2s - 2e),
(AI.7.2)
1~ = S. sin 4)~.
(AI.7.3)
It is then important to notice that for a same station and for the measurements considered above, we always have: I(i) A I ( j ) = _+l(k).
(AI.8)
APPENDIX 2
Calculation of A APPENDIX I
A2.1 Preamble
Calculation of l(i) (Fig. 1)
Let us consider quantities of the type:
AI.I Range measurement Range measurement positions the spacecraft in a plane orthogonal to axis SP. We thus have: I d = SP/ISPI.
(ml.1)
This formula provides a rigorous calculation of coordinates lid, 12a, 1~. A 2nd order approximation in (Rr/a~) gives: l~ = 1,
(A1.2.2)
Rv l~ = - - - cos ~bs sin (2 s - £p),
(A1.2.2)
as
Rr
1~ = - - -
Q =~c, cia~b~-~c, cja, b~a,bi i, I
= ~
I,j
cz cj(a~2 b i2 aiajbibj )
= ~ ci cj(a, bj - aj bi)2.
(A2.1)
I>j
sin q~s,
(A1.2.3)
as
A2.2 Application
Al.2 Bearing measurement A constant-elevation bearing deviation displaces the spacecraft along an axis parallel to Z ^ SP.
The quantity to be calculated is of the " Q " type defined above. We thus have: A=
We thus have: F=(Z
by developing we have:
A S P ) / I Z A SPI
(A1.3)
x [1~') l~j) - l~J ) l~i)]:
Here again, this formula provides a rigorous calcug g lation of l~, 12, l 3. However, we propose the following approximation: l~ = 0,
(A1.4.1)
l~ = S. sin q~s,
(A1.4.2)
l~ = - S . cos ~ sin (2s - 2p),
(A1.4.3)
with S = [1 --cos 2 ~bscosZ(£s - 2e)] ,.'2 (A1.5)
+ y N~,)~,~N,j)~j) l >/
x r/if)/(J) ,.2 . 3
_
1~ I l~'q 2.
This is the case of a single tracking station performing measurements of the range, elevation and bearing types. Formula (AI.8) thus makes it possible to write:
A= ~
N(i)W~,)N(j)W(j){[l~k)]2/4+ [l~k)]2}
i>j k # (i,j)
AI.3 Elevation measurement The motions of the "constant elevation" or "constant bearing" antenna are orthogonal. We thus take: IS= ia A is,
1
~ i~>iN(1)W(1)N(j) W(j)
(A1.6)
= UdWdUsW~U~W.
×~
c [l~q = [li~q~ - - 4 - - . k ~4N(k) W(k) N(k)W(k)t
(A2.2)