广播星历计算速度和加速度的扩展算法

Extension of GPS Broadcast Ephemeris to Determine Satellite Velocity and AccelerationWen Zhang, University of Leeds, Leeds, UK Mounir Ghogho, University of Leeds, Leeds, UKLuis Enrique Aguado, Advanced Digital Institute, Bingley, UKBIOGRAPHYWen Zhang received a B.Sc. degree in Physics from Central China Normal University (CCNU), Wuhan, China in July 2003, and an M.S. degree in Optical Engineering from National University of Defense Technology (NUDT), Changsha, China, in December 2005. Now she is a Joint-training doctorate candidate at both the University of Leeds, UK and the National University of Defense Technology, under the supervision of Professor Mounir Ghogho, Dr. Luis Enrique Aguado and Professor Bolong Gao. Her research currently focuses on the GPS/INS Integration.Prof. Mounir Ghogho received the M.S. degree in 1993 and the PhD degree in 1997, both in Signal and Image Processing from the National Polytechnic Institute (INP), Toulouse, France. He was an EPSRC Research Fellow with the University of Strathclyde, Glasgow, from September 1997 to November 2001. Since December 2001, he has been a faculty member with the school of Electronic and Electrical Engineering at the University of Leeds, where he is currently a Professor. His research interests are in signal processing for communications and networking. He is currently an Associate Editor of the IEEE TRANSACTIONS ON SIGNAL PROCESSING and a member of the IEEE SPCOM and SPTM Technical Committees.Dr. Luis Enrique Aguado has got over 15 years of research and project management experience in industrial and academic environments in the areas of satellite navigation and wireless communications. He is Sector Manager for the DSP Research Group at the Advanced Digital Institute in Bingley, UK. Before that he was Senior Research Fellow at the University of Leeds for 10 years, being Acting Director of the CAA Institute of Satellite Navigation from 2006 to 2007. Dr. Aguado has a PhD in Electrical Engineering from University of Manchester (UK) and a B.S., also in Electrical Engineering, from University of Mondragon (Spain).ABSTRACTThis paper presents an expansion of the algorithm provided in the GPS ICD to include the precise real-time computation in ECEF coordinate of the velocity and acceleration components of the satellite orbit besides the position components using the broadcast ephemeris. Precise satellite velocity andacceleration data are needed for applications using precise receiver Doppler measurements, such as GPS/INS integration. The expansion proposed in this paper improves the observed performance of an alternative derivation recently published. Ephemeris records are then expanded using this derivation and the resulting velocities and accelerations compared with the National Geodetic Survey (NGS) SP3 precise velocities and accelerations. This resulted in agreement in the order of one millimetre per second and one millimetre per second squared, respectively. It should be noted that using Kepler’s second law to calculate the first derivative of the true anomaly is not appropriate and it will decrease the agreement of the velocities to the order of one decimetre per second when compared with NGS SP3 precisevelocities. KEY WORDSGPS; Broadcast Ephemeris; Satellite Velocity; Satellite Acceleration.1. IntroductionGPS is used for applications where precise real-time satellite velocities and accelerations are needed aswell as satellite positions, such as GPS/INS integration. The satellite motion needs to be subtracted from the received signal dynamics to measure those of receiver antenna. Satellite’s ECEF positions can be derived from GPS Broadcast Ephemeris data by using the method described in the GPS interface control document ICD-GPS-200D [1]. However, it did not include algorithms for satellite ECEF velocity and acceleration determination. One method to derive the satellite velocity and acceleration is to use the rotation matrix method [2]. However, it should be noted that using Kepler’s second law to calculate the first derivative of the true anomaly is not appropriate as explained below and it will decrease the accuracy of velocity estimation. Another way to derive the satellite velocity and acceleration is to differentiate the satellite position with respect to time to calculate the velocity, and differentiate it a second time for the acceleration. This differential method can be implemented using the first order central differencing approximation of Taylor series [2], where one satellite velocity determination needs 3 successive satellite positions and one satellite acceleration determination needs 5 successive satellite’s positions.This paper shows how to calculate satellite velocities and accelerations using only the Keplerian elements and correction terms of the broadcast ephemeris, which is considered as an improvementof the rotation matrix method mentioned above. In the proposed method, just one set of Keplerian elements and correction terms of a specific epoch is required to compute both the satellite velocity and acceleration. The processing can be carried out using either post-processed/predicted ephemeris downloaded from the International GNSS Service (IGS), or real-time processing using GPS Broadcast Ephemeris parameters in the navigation data decoded by the receiver.In this paper, GPS ephemeris data definitions and computation of a satellite’s ECEF position vector are first introduced. Then, a set of equations to calculate the satellite velocity and acceleration is derived using Keplerian elements and correction terms included in the broadcast ephemeris. All the parameters and equations needed to compute a satellite position, velocity and acceleration are explicitly listed in tables for ease of use. An example is presented comparing the results obtained using the equations described with the results of precise NGS SP3 data. A comparison between the results obtained using Kepler’s second law to calculate the first derivative of true anomaly with NGS SP3 precise ephemeris is also carried out. Finally, a summary and conclusions are provided. 2.GPS Ephemeris Data Definitions andComputation of a Satellite’s ECEF PositionVectorLet us first list the symbols needed to compute the broadcast orbits. These are listed in Table 1 below using the notation adopted in [3].ucC Amplitude of cosine correction to argument of latitudeusC Amplitude of sine correction to argument of latitudercC Amplitude of cosine correction to orbital radiusrsC Amplitude of sine correction to orbital radiusicC Amplitude of cosine correction to inclination angleisC Amplitude of sine correction to inclination angleTable 1 GPS ephemeris data definitionsTable 2 provides the algorithm by which a GPS receiver computes the position vector of a satellite (sx,sy,sz) in the ECEF coordinate system from the orbital elements in Table 1 [3].(14) e e k e k t t 00)(Ω−Ω−Ω+Ω=Ω&&& Corrected longitude of node, eΩ&is the Earth rotation rate(15) k k p u r x cos = In-plane x position (16) k k p u r y sin =In-plane y position (17) k k p k p s i y x x Ω−Ω=sin cos cos ECEF x-coordinate (18) k k p k p s i y x y Ω+Ω=cos cos sin ECEF y-coordinate (19)k p s i y z sin =ECEF z-coordinateTable 2 Computation of a satellite’s ECEF Position VectorAs indicated in Equations (15) to (19) in Table 2, the in-plane (orbital plane) satellite position vectoris []Tp p p y x r 0=r and the satellite ECEF position vector is []Ts sss z y x r =r. s r r can beexpressed in terms of p r ras pp k x k z s r R r i R R r rr r =−Ω−=)()(, where the rotationmatrices are⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡ΩΩΩ−Ω=Ω−100cos sin 0sin cos )(kkk kk z R ,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−=−k kk kk x i i i i i R cos sin 0sin cos 0001)(, and ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡Ω−ΩΩΩΩ−Ω=−Ω−=k kk k kk k k k kk k k x k z i i i i i i i R R R cos sin 0sin cos cos cos sin sin sin cos sin cos )()( Thus, the ECEF coordinates Equations (17), (18)and (19) in Table 2 can be readily computed using the above matrix. Besides, the Keplerian orbital elements which define the orientation of the orbit are illustrated in Figure 1.Figure 1 The Keplerian Orbital Elements Defining the Orientation of the Orbit3. Computation of the Satellite ECEF VelocityThe first derivative of the rotation matrix R with respect to time, R&, is given by⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡−Ω−ΩΩΩ−ΩΩ−ΩΩΩ+ΩΩΩ+ΩΩ−ΩΩ−=k k k k k k k k k k k k kk k k k k k k kk k k k k kk k k k i i i ii i i i i i i i i i i i R sin cos 0cos cos sin sin sin cos cos sin cos cos sin sin cos sin sin cos cos sin &&&&&&&&&&&&&Differentiating Equations (4) and (5) with respect to time and substituting the result of the first derivationfor kM &in the second derivation, we obtain )cos 1(k k E e n E −=& (20) Differentiating with respect to time the expression for )sin(k ν in Equation (6) and then using the expression for )cos(k ν in Equation (6) yields)cos 1(12k k k E e e E −−=&&ν (21) In reference [2], the author presented a different equation for calculating k ν& using Kepler’s second law (the radius vector sweeps out the equal areas in equal time). In this case the derivative of the true anomaly k ν&is given by: 2221k k r n e a ⋅−=ν& (21’) Using Equation (21’) instead of (21) will result in asatellite ECEF velocity error when the time from ephemeris epoch k t is not equal to zero. The largerthe value of k t , the larger the satellite ECEF velocity error. The main reason for this is that the semi-major axis a could not be corrected and compensated for when used at ephemeris epochs k t beyond the reference epoch e t 0 in the algorithm. Meanwhile, the effect on the satellite ECEF acceleration is negligible. This is illustrated with anexample in Table 7 of Section 5. Assuming the argument of the perigee ω is time-invariant, differentiating Equation (7) with respectto time shows that k kνφ&&≈. Further, substituting the right hand side of Equation (10) for k i δ in equation(13) and taking the first-order derivative of the resulting equation, we obtain the first-order derivative of the corrected inclination as followskk ic k is kk ic k is k v C C i C C i i &&&&&))2sin()2cos((2))2sin()2cos((2φφφφφ−+≈−+= (22)where the correction terms can be readily calculatedfrom the expression for k ν& in Equation (21). Substituting the right hand side of Equation (8) for k δφ in Equation (11) and taking the first-order derivative of the resulting equation, we obtain thefirst-order derivative of the corrected latitude, wehave thatkk uc k us k kk uc k us k k v C C v C C u &&&&&))2sin()2cos((2))2sin()2cos((2φφφφφφ−+≈−+= (23)Substituting the right hand side of Equation (9) for k r δ in Equation (12) and taking the first-orderderivative of the resulting equation, we obtain the first-order derivative of the corrected radius, wehave thatk krc k rs k k kk rc k rs k k k v C C E E ae C C E E ae r &&&&&))2sin()2cos((2sin ))2sin()2cos((2sin φφφφφ−+≈−+= (24) Differentiating Equation (14) for the correctedlongitude of node with respect to time, we have that ek Ω−Ω=Ω&&& (25) The derivatives of the in-plane x and y positions in equations (15) and (16) with respect to time give[][]Tkk k k k k k k k k T pp p u u r u r u u r u r y x r 0cos sin sin cos 0&&&&&&&r +−==that isk k k k k p u u r u r x sin cos &&&−= (26) k k k k k p u u r u r y cos sin &&&+= (27) Differentiating the rotated vector p s r R r r r = with respect to time, we have []pp T s s sr R r R yx r &r r &&&&r +==0 that isk k p k p k k k kk k p kk p s i y xi i i y x x cos sin cos )sin sin cos cos (sin Ω−Ω+Ω−ΩΩ−ΩΩ−=&&&&&& (28) k k p k p k k k kk k p kk p s i y xi i i y x y cos cos sin )sin cos cos sin (cos Ω+Ω+Ω+ΩΩ−ΩΩ=&&&&&& (29) k p k k p s i y i i y zsin cos &&&+= (30)4. Computation of Satellite’s ECEFAccelerationTaking the double-differenced terms k i && and k Ω&&terms as zero, the second derivative of the rotation matrix R is⎥⎥⎥⎦⎤−ΩΩ+Ω+ΩΩΩ+Ω+Ω−⎢⎢⎢⎣⎡−ΩΩ+Ω+Ω−ΩΩ−ΩΩ+Ω+ΩΩΩ−=k k k k k k k k k k kk k k k k k k k k kk k k k k k k k k k k k k k k k k kk i i i i i i i i i i i ii i i i i i i i R cos cos sin 2sin cos )(cos cos 2sin sin )(sin 0sin sin 2cos cos )(sin sin cos 2cos sin )(cos 222222222222&&&&&&&&&&&&&&&&&&&&&&In Classical Mechanics, considering the two-body problem of Earth and the satellite, the satellite’s acceleration in the nominal orbit plane is[]33kp pp Tp p p r r r r yxr r rr &&&&&&r ⋅−=⋅−==μμwhere GM =μ, with G being the universal gravitational constant and M being the Earthmass, respectively. Therefore, the in frame acceleration components can be expressed as 3k p p r x x⋅−=μ&& (31) 3k p p r y y⋅−=μ&& (32) Differentiating the rotated velocityvector pp s r R r R r &r r &&r += with respect to time, we have[]pp p T s s s s r R r R r R z y x r &&r &r &r &&&&&&&&&&r ++==2 That isk k p k p k k k kk k p kkpkkkk k k k k p k k p s i yx i i i y x i i i i y x x cos sin cos )sin sin cos cos (2sin 2)sin cos 2cos sin )((cos 222Ω−Ω+Ω−ΩΩ−ΩΩ−ΩΩ+Ω+Ω+ΩΩ−=&&&&&&&&&&&&&&&& (33)k k p k p k k k kk k p k k p k k k k k k k k p k k p s i yx i i i y x i i i i y x y cos cos sin )sin cos cos sin (2cos 2)sin sin 2cos cos )((sin 222Ω+Ω+Ω+ΩΩ−ΩΩ+ΩΩ−Ω+Ω−ΩΩ−=&&&&&&&&&&&&&&&& (34) k p k k p k k p s i yi i y i i y z sin cos 2sin 2&&&&&&&++−= (35) Table 3 summarizes the algorithm to compute the velocity vector (s s s z y x&&&,,) and acceleration vector (s s s zy x &&&&&&,,) of a GPS satellite in the ECEF coordinate system from the orbital elements inTable 1.5.Performance ResultsThe above satellite position, velocity and acceleration determination algorithm, described from equations (1) to (35), has been implemented using Matlab code. Table 4 is a set of input broadcast ephemeris data corresponding to the simulated time as employed in the Matlab program.The data in table 5 has been downloaded from the NGS as given in [2]. It lists both velocities and accelerations for PRN 07 on 20th August 2002 from time 22:19:01 to 22:19:10. The velocities are from the SP3 precise ephemeris and the accelerations are obtained using the first-order central difference of a Taylor series approximation of the NGS SP3 velocities. We can use Table 5 as a precise reference of velocities and accelerations of these time specifications.Time)/(s m V x )/(s m V y)/(s m V z)/(2s m A x )/(2s m A y )/(2s m A z 22:19:01 329.9513 -888.5965 -2997.0555 0.2055 -0.3022 0.1413 22:19:02 330.1568 -888.8987 -2996.9141 0.2054 -0.3022 0.1414 22:19:03 330.3622 -889.2009 -2996.7727 0.2054 -0.3022 0.1414 22:19:04 330.5676 -889.5032 -2996.6312 0.2053 -0.3023 0.1415 22:19:05 330.7729 -889.8055 -2996.4897 0.2053 -0.3023 0.1415 22:19:06 330.9782 -890.1078 -2996.3481 0.2052 -0.3023 0.1416 22:19:07 331.1834 -890.4101 -2996.2064 0.2052 -0.3023 0.1417 22:19:08 331.3886 -890.7124 -2996.0646 0.2051 -0.3023 0.1418 22:19:09 331.5937 -890.0147 -2995.9228 0.2051 -0.3024 0.1418 22:19:10 331.7988 -891.3171 -2995.7809 0.2050-0.30240.1419Table 5 Satellite Velocity and acceleration from NGS SP3 precise position and velocityTable 6 are satellite velocities and accelerations calculated using the derived equations (28), (29), (30), (33), (34) and (35) in Table 3. The velocities and accelerations in Table 6 have been compared with the precise reference in Table 5 and the agreements are better than 1 mm/s and 1mm/s 2.Time)/(s m V x )/(s m V y )/(s m V z)/(2s m A x )/(2s m A y)/(2s m A z 22:19:01 329.9516 -888.5963 -2997.0556 0.2054 -0.3022 0.1413 22:19:02 330.1570 -888.8985 -2996.9143 0.2054 -0.3022 0.1413 22:19:03 330.3624 -889.2007 -2996.7728 0.2053 -0.3022 0.1414 22:19:04 330.5678 -889.5030 -2996.6314 0.2053 -0.3022 0.1415 22:19:05 330.7731 -889.8052 -2996.4898 0.2052 -0.3022 0.1415 22:19:06 330.9784 -890.1075 -2996.3482 0.2052 -0.3023 0.1416 22:19:07 331.1836 -890.4098 -2996.2065 0.2051 -0.3023 0.1417 22:19:08 331.3888 -890.7122 -2996.0648 0.2051 -0.3023 0.1417 22:19:09 331.5939 -891.0145 -2995.9229 0.2050 -0.3023 0.1418 22:19:10 331.7990 -891.3169 -2995.7811 0.2050-0.3023 0.1419Table 6 Satellite velocities and accelerations calculated using equations in Table 3Table 7 shows the satellite velocity andacceleration calculated using the formulas in Table 3, but using equation (21’) instead of (21) when calculating k ν&. The velocities and accelerations in Table 7 have been compared with the precise reference in Table 5 and the agreement of the velocities decreased to better than 1 dm/s, while the agreement of the accelerations are still better than 1mm/s 2. This decreased accuracy of the velocities is because the semi-major axis a change over the time epochs k t could not becorrected and compensated in the algorithm.Meanwhile, the impact of the change of the semi-major axis a on the acceleration calculation is negligible.Time)/(s m V x)/(s m V y )/(s m V z )/(2s m A x )/(2s m A y)/(2s m A z 22:19:01 329.9635 -888.5622 -2997.0121 0.2054 -0.3022 0.1413 22:19:02 330.1689 -888.8645 -2996.8708 0.2054 -0.3022 0.1413 22:19:03 330.3744 -889.1667 -2996.7294 0.2053 -0.3022 0.1414 22:19:04 330.5797 -889.4690 -2996.5879 0.2053 -0.3022 0.1415 22:19:05 330.7850 -889.7713 -2996.4464 0.2052 -0.3022 0.1415 22:19:06 330.9903 -890.0735 -2996.3048 0.2052 -0.3023 0.1416 22:19:07 331.1955 -890.3759 -2996.1631 0.2051 -0.3023 0.1417 22:19:08 331.4007 -890.6782 -2996.0214 0.2051 -0.3023 0.1417 22:19:09 331.6058 -890.9805 -2995.8796 0.2050 -0.3023 0.1418 22:19:10 331.8108 -891.2829 -2995.7377 0.2050-0.30230.1419Table 7 Satellite velocities and accelerations ignoring the change on the semi-major axis6.ConclusionsIn summary, this paper has presented an expansion of the ICD-GPS-200D algorithm to include the computation of the velocity components and acceleration components of the satellite additionally to the position components in the ECEF coordinate system using the broadcast ephemeris. The equations to calculate velocity and acceleration are summarized in Table 3. Computed acceleration and velocity for an example data set have been compared with National Geodetic Survey SP3 precise velocities and accelerations resulting in an agreement to one millimetre per second and one millimetre per second squared, respectively. With such a high accuracy satellite velocities and accelerations in real time, precise, Doppler based, GPS velocity and acceleration measurements can be obtained with simple single frequency stand-alone set-ups that may be used in different applications. We have proposed here a new algorithm for the derivation of the satellite velocity which was shown to provide millimetric performance, therefore outperforming existing methods.7.AcknowledgementsThe first author would like to thank Chinese Scholarship Council which sponsors her studies in the UK.8.References1.ARINC Engineering Services LLC., IGS-GPS-200, Revision D,Navstar GPS SpaceSegment/Navigation User Interfaces, 7thDecember 20042.Zhang, J. et al, GPS Satellite Velocity andAcceleration Determination using theBroadcast Ephemeris, Journal of Navigation(2006), 59:2:293-305.3.Kaplan, E. D. (Ed.), Understanding GPS:Principles and Application, 1996, ArtechHouse Publishers。