Satellite Orbit computing

 

1. GPS observation volume

While observing the phase and pseudo-distance data, the receiver also records the broadcast and forecast calendars. GPS signals can also be received to obtain time reference signals with nanosecond precision.

Due to the wide range of receiver models, the vendor designed different data formats. In order to uniformly use the data of different receivers internationally, a RINEX (REceiver
InDependentExChange format). Version 2 is used currently. The following describes the broadcast calendar files, observation data files, and ground meteorological data files in RINEX 2 format.

The naming rules for GPS data files in RINEX 2 format are as follows:

Where :~ The station name in four bytes;

~ The accumulation day of the first group of data observation time in the file (for example, January 1, and January 2 );

~ The sequence number of a certain type of files received by the station on that day. 0 indicates that only one file exists;

~ Year in double digits (for example, 96 represents 1996 );

~ File Type:

~ Observation data file;

~ Broadcast the calendar file;

~ Meteorological data files on the ground.

For ease of communication, the GPS data files in RINEX 2 format are marked with ① no logo; ② ASCII code; ③ each record is 80 characters in length and the block size is 8000; recorded on tape, the first file on the tape is the directory of all files. However, at present, international organizations such as IGS use communication methods (Internet) to quickly extract data from global GPS observation stations for years and store the data in large computers, attackers can extract data over the Internet.

It should be noted that in the GPS data in RINEX 2 format, the time is measured in gpst, that is, an integer jump seconds to UTC.

 

Broadcast a calendar File

After the receiver locks the satellite and solves the c/a code, it can obtain the broadcast calendar, that is, the satellite coordinate calculation parameter. in real-time GPS applications, it is essential, broadcast calendars are also used for post-processing of most engineering network observation data. RINEX

 

The broadcast calendar file in the format of 2 is shown in Table 2.1.1 below. As an example, the table shows the broadcast calendar data of the prn9 and prn17 satellites. PRN indicates the pseudo-random number of the GPS satellite, in some cases, GPS satellites use NASA (NRational
AEronautics andSPaceADministration.

Table 2.1.1 broadcast calendar files in RINEX 2 format

2 navigation data RINEX version/Type Ephtornx version 1.09 29-nov-95 PGM/run by/date . 1676d-07. 2235d-07-. 1192d-06-. 1192d-06 ion alpha . 1208d + 06. 1310d-07-. 1310d + 06-. 1966d + 06 ion Beta . 133179128170d-06. 107469588780d-12 552960 39 DELTA-UTC: A0, A1, T, W 10 leap seconds End of Header 9 94 10 21 8 0 0.0-0.103851780295d-04-0.909494701773d-12 0.20.0000d + 00 0.720000000000d + 02 0.106062500000d + 03 0.476841277575d-08 0.132076112444d + 01 0.548548996449d-05 0.312971079256d-02 0.747293233871d-05 0.515371790504d + 04 0.460800000000d + 06 0.558793544769d-07-0.229012694900d + 01-0.130109160446d-07 0.950477774712d + 00 0.229593750000d + 03-0.491558992251d + 00-0.819034084998d-08 0.233938313166d-09 0.20.0000000d + 01 0.77000000000d + 03 0.20.0000d + 00 0.700000000000d + 01 0.000000000000d + 00 0.1396980000192d-08 0.328000000000d + 03 0.20.00000000d + 00 17 94 10 21 8 0 0.0-0.635907053947d-04-0.909494701773d-12 0.20.0000d + 00 0.228000000000d + 03 0.167187500000d + 02 0.0000946255961d-08 0.104717256943d + 01 0.566244125366d-06 0.790101150051d-02 0.883266329765d-05 0.515109790649d + 04 0.460800000000d + 06 0.279425772385d-07 0.925235566518d + 00-0.203028321266d-06 0.97044258460d + 00 0.214781250000d + 03 0.199825790573d + 01-0.789747200969d-08 0.40000097220000d-09 0.00000000000d + 01 0.77000000000d + 03 0.000000000000d + 00 0.700000000000d + 01 0.000000000000d + 00 0.1396980000192d-08 0.484000000000d + 03 0.20.00000000d + 00

The first few behavior headers in the table above, 60th ~ The 80 characters indicate the corresponding line. For example, "2" in the first line is the RINEX version number, and "navigation data" indicates that the file type is broadcast calendar; the second line is the unit, executor, and date on which the file is generated. The header ends with "End of Header", and the header ends with a blank line. The two lines marked as "Ion Alpha" and "Ion beta" in the header indicate that the corresponding parameter is the correct parameter of the electron delta (this parameter is useless when the differential model is used for relative positioning ), the line annotated as "Leap seconds" provides the number of seconds between gpst and UTC. An infinite number of annotation lines can be inserted before the end of the header. The description of the annotation line is "comment ".

After the header ends, each eight rows indicates the broadcast star calendar of a satellite. The first row contains the PRN number (pseudo-random number) of the satellite ). For example, rows 8 to 15th in Table 2.1.1 are the broadcast calendars of the satellite prn9. The meanings of the corresponding quantities are shown in Table 2.1.2.

 

Table 2.1.2 meaning of broadcast calendar data in RINEX 2 format

Satellite PRN, year, month, day, hour, minute, second ,,,

,,,

,,,

,,,

,,,

,,,

,,,

 

Table:

~ Satellite clock difference constant, time second;

~ Satellite Clock Deviation drift, in seconds/in seconds;

~ Satellite clock difference drift rate, in seconds/in seconds;

~ Long-Term Variation of the flat angle (near-location parameter), radians/time seconds;

~ Reference point and radian of the moment;

~ Flat heart rate;

~ Square root of the long half axis, meters;

~ Refer to the time when the accesskey is raised to the chijing, radian;

~ Refer to the time-point orbital inclination, radians;

~ Near-location corner, radians;

~ The Long-Term Variation of the elevation point chijing in the Equator plane (mainly caused by the 2nd Order Harmonic Coefficient () of the Earth's gravitational field), radians/hour seconds;

~ The amplitude and radians of the cosine are corrected periodically in the direction of orbital trace at the reference time of the star calendar;

~ Corrected the amplitude and radians of the sine in the direction of orbital trace at the reference time of the star calendar;

~ The amplitude of the cosine, meter, is corrected on the radial direction of the orbit at the reference time of the star calendar;

~ Corrected the amplitude of the sine in the radial direction of the orbit at the reference time of the star calendar;

~ The amplitude and radians of the cosine are corrected in the cycle of the orbital tilt (similar to the normal direction) at the reference time of the star calendar;

~ Corrected the amplitude and radians of the sine in the cycle of the orbital tilt (similar to the normal direction) at the reference time of the star calendar;

~ The reference time (the number of seconds in the week) and the second of the calendar;

~ The age of the calendar data;

~ Orbital dip change rate, radians/hour seconds;

~ Indicates the code that exists on;

~ Number of GPS days per week;

~ Upper p-code pseudo-distance indication;

~ The accuracy indicators of the current broadcast stars, meters;

~ Satellite health indicators;

~ The time in seconds for modifying the group delay;

~ Age of satellite clock data;

~ Information Transmission Time (related to the satellite signal decoding received by the receiver), in seconds.

 

 

✓ Observation File

The phase and pseudo-range observed values measured by the receiver are recorded in the observation file. The following table 2.1.3 is the example file. There are two historical observations.

 

Table 2.1.3 observation data files in RINEX 2 format

2 Observation Data gps rinex version/Type Tb2rnx xxxxxxxxxxxx 95-07-20 22:20:20 PGM/run by/date Turbo ⅱ RINEX formatter version: 95.5.19 comment Mode: static comment Comment 1 marker name Wan observer/Agency 153575902 Turbo ⅱ Production Unit rec #/type/vers 153575902 Turbo ⅱ ant #/Type -2852389.2261 4650364.5453 3293350.3434 approx position XYZ 0.1120 0.0000 0.0000 antenna: Delta H/e/n original slant height (m): 0.1120 comment 1 1 wavelength fact L1/2 5 C1 L1 L2 D1 P2 #/types of observ End of Header 95 7 19 0 21 0.0000000 0 5 27 28 15 31 19 21805891.33516-0.28716 0.07713 0.00000 21805895.08513 23733969.59714-0.19314 0.66912 0.00000 23733973.16712 20746137.25719 0.01719 0.54216 0.00000 20746140.99416 20903455.21318-0.28018 0.99414 0.00000 20903458.41214 20681238.04917-0.32817-0.25610 0.00000 20681238.04910 95 7 19 0 21 30.0000000 0 5 27 28 15 31 19 21794932.63606-57587.90306-44873.38602 0.00000 21794935.76502 23750223.73904 85421.79604 66563.21801 0.00000 23750225.60401 20736273.69909-51834.40309-40389.90405 0.00000 20736277.71705 20913805.23508 54389.67808 42382.77504 0.00000 20913808.93304 20682314.76407 4906.19307 3823.00700 0.00000 20682314.76400

 

Like broadcast stars, the observation file also has a header that ends with an "End of Header" or a blank line ~ The 80 characters indicate the content of the row, as shown in table 2.1.3, the header of the observation file lists the RINEX version number, the units, personnel, point number, name, observer and unit, receiver number type version, antenna number type, station approximate coordinates, and antenna L1. relationship between the phase center and the point (number of deviations from the above, to the east, and to the north), sampling interval, wavelength factor, number of observed data types and observed data types, the first time of observation, the last time of observation, and the number of the measured satellite, the annotation line described as "comment.

The wavelength factor is 1, which indicates that the carrier is restored using related technologies. The full-week ambiguity and the number of weeks of loss can only be integers. The wavelength factor is 2, which indicates that the carrier is generated using the square technology, the integer and week loss may be multiples of the week 0.5. Of the observed types:

L1 ~ Carrier Phase on L1;

L2 ~ Carrier Phase on L2;

C1 ~ The pseudo-distance of the C/A code on L1;

P1 ~ P-code pseudo-distance on L1;

P2 ~ P-code pseudo-distance on L2;

D1 ~ Frequency Variation on L1;

D2 ~ L2 frequency change.

 

After the header ends, it is the observation data of each calendar element. Each calendar element data is composed of one calendar time row and multiple observation data rows. The content of each calendar element time row is:

Year, month, day, hour, minute, second, quality mark, Satellite Number,..., Clock Difference

Where:

① The number of satellites refers to the total number of satellites observed by this element, followed by the number sequence of the observed satellites.

When the quality flag is:

~ The Historical Observation Data is normal;

~ There is a power outage between the previous calendar and the current calendar;

2 ~ Start to move the antenna;

~ Start observation at the new station (then a new point number will appear );

~ New header information will appear below;

~ Other external events;

~ Weekly failure indication.

② Satellite Number:

: Satellite System:

~ GPS;

~ GLONASS system (GPS-like positioning system developed by the former Soviet Union );

~ Doppler satellite positioning system.

: If it is a GPS system, the PRN number;

If it is a GLONASS system, the channel number;

If other satellite systems have two digits.

③ The Clock Difference (option) in the calendar yuan line is the receiver clock difference, which is between 68 and ~ 80 bits. If this option exists, the following correction should be made to the historical time and the measured pseudo-range phase:

Calendar hour = give the calendar hour-Clock Difference

Pseudo-distance = measured pseudo-distance-Clock Difference x light speed

Phase = measured phase-Clock Difference x frequency

 

The rows after the element row (the number of rows equals to the total number of satellites) are the observed values of each satellite. The number of data in each row is equal to the observed types in the header. The data types are arranged according to the observed types in the header, each observed value is followed by a one-digit week loss indicator and a one-digit signal strength indicator. The unit of phase observation is Week, and the unit of pseudo distance is meter. When a data is not tested, the value is calculated as 0.0 or a space.

Weekly loss indicator (value range: 0 ~ 7) data is normal when it is null or empty.

When the signal strength indicator is:

Or space ~ Normal;

~ Weakest Signal;

~ Ideal Signal-to-Noise Ratio;

~ The strongest signal.

 

Meteorological documents

In high-precision observation, meteorological data is sometimes required to be recorded. The RINEX 2 meteorological data file is shown in Table 2.1.4 below:

 

Table 2.1.4 weather data files in RINEX 2 format

2 meteorological data RINEX version/Type Tb2rnx xxxxxxxxxxxx 95-07-20 22:20:20 PGM/run by/date 1 marker name 3 pr td hr #/types of observ End of Header 95 7 19 0 10 00 987.1 10.6 89.5 95 7 19 1 10 00 987.2 10.9 90.0 95 7 19 2 10 00 987.1 11.5 89.0

 

The table header and terminator are similar to those in the calendar and observation data files. They record the meteorological observation time (year, month, day, hour, minute, second), pressure (maba) around the station, and dry temperature (degree Celsius), relative humidity (percentage ).

 

2. Computing satellite positions by broadcast calendars

If you want to calculate the spatial coordinates of a satellite at a time point, read the relevant parameters of the Satellite Broadcast calendar according to the previous broadcast calendar format, and perform the following steps:

1. Find the long Half Axis

2. Calculate the angular speed

It is the gravity constant of the Earth.

3. Calculate the time difference from the time needed to the reference time

4. Correct the normal angle speed

5. Calculate the approximate point angle

6. iterative computation based on the following formula: near point

7. Calculate the real-near-point angle from the following two formula:

8. Calculate latitude Parameters

9. Period correction items

10. Calculate the latitude parameter after correction.

11. Calculate the correct forward longitude

12. Calculate the corrected inclination

13. Coordinate Calculation of satellites in the orbital plane

14. Correct the longitude of the elevation point

The formula is the earth's rotation rate.

15. Finally, the coordinates of the satellite in the earth-solid system are calculated.

When calculating the satellite position based on the above process, it should be noted that the satellite coordinate is obtained to represent the time when the satellite signal is transmitted in the earth-solid system. The time when the calculation is based on gpst, it is obtained after the signal receiving time minus the propagation time delay (iteration is required, as shown in the following figure. When the value is greater than 302400 seconds, because the number of seconds corresponding to a GPS week is reduced to 604800 seconds in the middle, 604800 seconds should be added when the value is less than 0 seconds.

Generally, the coordinates obtained from the preceding process are represented in WGS84. However, due to the large error, the coordinates may exceed 100 in some cases, therefore, it can be considered that the coordinates are represented in the cts of any earth-solid reference system.

 

16. Calculate the time when a GPS signal is sent from a satellite (minus the delay of the receiver transmission time)

GPS observations are the functions of satellite and receiver oscillator frequencies and propagation time latencies. propagation time latencies include geometric propagation time latencies and atmospheric propagation time latencies, namely:

2.2.15

The parameters required for GPS measurement are included. The following describes the determination.

Strictly speaking, the determination should be calculated in the coordinate system with the center of the solar system as the origin in the framework of general relativity. However, in practical use, the computation in the central gravity inertial coordinate system can meet the precision requirements. The computation can also be equivalent in the coordinate system.

In positioning applications, satellite orbit is known, station coordinates are unknown, but approximate values are known (even if there is no other known information, the approximate coordinates of stations can also be obtained from the observation file). The iterative process of data processing is to calculate the station coordinates from the approximate value to the accurate value. When positioning and determining the orbit, the approximate value of the satellite orbit is also known and can be obtained from the broadcast calendar file without other known information.

If the geometric propagation time delay of the satellite signal received by the station is calculated in the earth-solid coordinate system, the electromagnetic wave propagation law is as follows:

In the formula, the coordinates of the satellite at the time of signal transmission are represented in the coordinate system of the Earth, the coordinates of the station in the earth-solid coordinate system, and the propagation speed of the electromagnetic wave in the vacuum (light speed, 299792458 meters/second ). Since the time when the signal is included in formula 2.2.16 is unknown, it must be calculated before it can be obtained. Therefore, an iteration process is required. The specific steps are as follows:

① Based on the geometric propagation time delay value obtained in the first iteration (in the first iteration, the second can be assumed because the height of the satellite from the ground is approximately 20000 kilometers), and the satellite signal transmission time is obtained:

The atmospheric propagation time delay in the formula can be calculated based on the atmospheric correction model (see later ).

② Obtain the satellite coordinate at the moment based on the satellite calendar table interpolation.

③ Calculation based on the above formula:

④ When the difference is greater than the limit difference, return to ①. When the difference is less than the limit difference, it is deemed that it has been converged and replaced by the value, plus the instant, the limit deviation should be a value less than or equal to, because the influence of the phase is the product of the phase and the frequency, and the magnitude of the frequency is, the measurement accuracy of GPS reception phase can reach 0.01 weeks.

The above propagation time delay is calculated in the earth-solid system, and can also be calculated in the inertial system. It can be easily selected. If the star calendar is expressed in the earth-Solid System (for example: broadcast stars and precision calendars in NGS format). If no orbit determination is performed, it is easier to choose the Geosystem for calculation. If the orbit or calendar is required to be expressed in the inertial system, the calculation in the inertial system is more convenient. During the calculation, the coordinates of the station must be rotated to the inertial system in each calendar element, this avoids the computation of rotating the star calendar from the inertial series to the earth-solid series with a large amount of computing.

 

17. Central atmospheric correction

Atoms and molecules are in neutral state from the earth's surface to the atmosphere 80 km higher than the ground. Therefore, they are called neutral atmosphere and can also be called the trotters, some documents have divided it into two layers, however, since GPS observation only uses observations with an elevation angle greater than 15 degrees, there is no difference in effect between one layer and two layers. The neutral atmosphere increases the propagation time of electromagnetic waves, which is called the neutral atmospheric delay. At the top of the sky, it can reach about 2.5 meters and increase with the increase of the height angle. The neutral atmospheric delay is divided into two parts, which are called dry term due to the polarization displacement of all atmospheric molecules in the atmosphere, and wet term caused by the dipolar distance of water molecules. The items are relatively stable and can be well corrected using a suitable model. The additional delay caused by wet items is much smaller, only dozens of centimeters, but the changes are irregular, without a high-precision correction model, an expensive moisture radiation machine can be used for determination. In most cases, parameters to be evaluated are introduced for processing.

Models such as Saastamoinen 1973, Marini, and Chao can be used to correct the atmospheric delay at the top of the sky. In general, correction must also be multiplied by a ing function whose distance is the variable, ing functions include Lanyi, CFA, and other models. In GPS data processing, the observed data volume is very large. Generally, we only use the observed data with an elevation angle greater than 15. Therefore, most of the existing Skytop correction and ing function models can be used. The commonly used Saastamoinen 1973 model is as follows:

The maximum latency of the task is as follows:

Maximum humidity delay:

Because the wet term model in formula 4.1 is inaccurate, A undetermined parameter can be introduced to change formula 4.1:

4.1.2

In the above formula:

The above is the ground pressure (maba), the ground temperature (Celsius), the water pressure (maba), the relative humidity, the central latitude of the station, it is the geolevel height (km) of the station ).

The water vapor spectrometer can be used to accurately determine the influence of accumulated water vapor and clouds in the atmosphere on the Growth of electromagnetic wave propagation paths. The model is as follows:

To indicate the light temperature in the sky in any direction measured by the vapor analytics, the electromagnetic wave propagation path growth (that is, the product of the front and the speed of light) caused by the atmospheric humidity can be calculated by the following formula:

In the preceding formula, the sum is the two-band frequency, and corresponds to the bright temperature. It can be approximately equivalent to the radiation temperature of the cosmic background. It is the weight function corresponding to the propagation path, and for the absorption coefficient of water vapor and oxygen, and can be determined by a standard atmospheric model and the temperature measured on the ground, which is similar to (ground temperature ,). The accuracy measured by the water vapor spectrometer can reach about 1.5mm.

 

18. Correction of the lower extremity

Earth Surface 60 ~ 1000 km of the atmosphere, due to the solar radiation, the atoms in the atmosphere are ionization into a large number of positive ions and electrons, forming the electron. The time delay for electromagnetic waves to pass through the electron sphere is as follows:

4.1.5

The formula is the observed frequency, the speed of light, the total Electronic Content in the transmission path, and the classic electronic radius.

Because the 4.1.5 formula is related to the solar radiation pressure, the day and night may differ by an order of magnitude, and it is difficult to use a model for exact representation. However, we can see that the propagation time delay is inversely proportional to the square of the frequency. Therefore, we can use dual-band observation to eliminate its impact.

For example, if, and are the pseudo-distance and phase of a GPS receiver measured at the same time and in the two bands, the effect of the electrons on them is ,,,, here, the number of observations of a receiver on the same satellite in the same calendar element is constant, and is the frequency of two bands. Therefore, in the relative positioning of long-distance GPS, the observed amount after linear combination of phase observations in two bands is usually used as the observed value.

Formula, that is, the ratio of L1 to L2 frequency.

 

3. Pseudo-Range Differential Positioning Model

The pseudo distance from the reference station R to the GPS satellite J is

The real distance between the reference station R and the J-th satellite is the distance deviation caused by the GPS satellite calendar error, which is the deviation between the receiver clock and the GPS time system; it is the deviation between the clock of the J satellite and the GPS time system, the distance deviation caused by the time delay of the elliptic plane, the distance deviation caused by the frequency delay of the trotters, and the propagation speed of the electromagnetic waves.

The real distance can be accurately calculated based on the known coordinates of the Reference Station and the GPS satellite calendar. While the pseudo distance is measured by the reference station receiver, the pseudo distance is changed to a positive value.

At the same time as the benchmark receiver's pseudo-distance measurement, the receiver K of the mobile station also measured the pseudo-distance of the J-satellite.

We bring the pseudo-distance correction value measured by the benchmark station to the upper formula (that is, the upper two formula are added:

When the distance between the mobile station and the benchmark station is within a certain range, we can consider it as follows:

 

The above formula is changed:

(6) There are four unknowns in formula. They are the three-dimensional coordinates of the mobile station K and the correction items caused by the clock difference of the GPS receiver. If the reference station and the mobile station view more than four satellites, the error equation can be established based on the above formula:

Formula () is the three-dimensional approximate coordinate of the station K; () is the three-dimensional coordinate of the moment when the satellite J emits a signal, which can be obtained according to the calculation of the satellite calendar; j = ,.... N =.

Then, based on the least square method, the three-dimensional coordinates of the calendar element and the errors in the Three-dimensional coordinates are solved.