A GNSS navigation satellite PCO correction method based on sun position and dual-frequency combination

Abstract

The invention relates to the field of satellite positioning technology, in particular to a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination, which comprises the following steps: reading GNSS satellite PCO information, calculating the sun's position in a geocentric inertial system by an analytical method, calculating Z, S, Y, and X vectors, calculating a rotation matrix S2C, reading PCO information of a corresponding satellite, rotating the PCO to an ECI coordinate system by S2C, calculating the combined PCO of the GNSS satellite by dual-frequency combination, and then correcting the combined PCO to the GNSS satellite position. The invention improves the accuracy of calculating GNSS satellites, and has important application value in fields with high accuracy requirements for GNSS satellite positions in remote sensing detection such as GNSS occultation detection, GNSS reflection detection, and other low-orbit satellite payloads.

Description

Technical Field

[0001] The invention relates to the technical field of satellite positioning, in particular to a GNSS navigation satellite PCO correction method based on solar position and dual-frequency combination. Background Art

[0002] GNSS is the Global Navigation Satellite System, which includes the United States' GPS, China's BDS, Russia's GLONASS and Europe's GALILEO.

[0003] GNSS satellites can broadcast GNSS navigation signals and messages to the outside world. GNSS receivers on low-orbit satellites can receive signals transmitted by GNSS satellites and perform capture and tracking to obtain observation information such as pseudorange, carrier phase, and signal-to-noise ratio.

[0004] In the field of GNSS remote sensing, for example, in technologies such as GNSS radio occultation detection and GNSS reflection detection, it is necessary to determine the geometric relationship or altitude angle relationship between low-orbit satellites and navigation stars, which requires the precise position of GNSS navigation satellites.

[0005] The official website of the International GNSS Service IGS or the Wuhan University IGS website publishes the precise orbit and clock error information of GNSS satellites. Usually, the orbit is the position of the satellite's center of mass. Generally, GNSS satellites are large, and the point where they actually receive signals is the phase center of the receiving antenna. The satellite's center of mass is about 3 meters away from the phase center of the receiving antenna. Therefore, when accurately calculating the GNSS satellite position, it is necessary to correct the deviation of the GNSS satellite's antenna phase center relative to the center of mass. At the same time, it is also necessary to consider the impact of the GNSS satellite's attitude on it.

[0006] In the prior art, the GNSS satellite orbits published on the IGS official website or the Wuhan University IGS Center are all accurate position information at the centroid of each satellite, which does not represent the position at that point when it receives the signal. Using the position at the satellite centroid instead of the position at the GNSS signal receiving antenna will produce a large error. In the GNSS occultation detection technology, the GNSS occultation position accuracy requirement is 20cm. The deviation of the GNSS antenna phase center relative to the satellite centroid will cause the calculated position deviation of the received signal to exceed 20cm, thereby causing the inversion error of the neutral atmosphere occultation event to be too large. This patent provides a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination.

[0007] Therefore, in order to solve the above problems, the present application proposes a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination, which solves the problem of inaccurate GNSS satellite calculation and improves the accuracy of calculating GNSS satellites. It has important application value in remote sensing detection such as GNSS occultation detection, GNSS reflection detection and other low-orbit satellite payloads, which have high accuracy requirements for GNSS satellite positions. Summary of the invention

[0008] The purpose of the present invention is to overcome the shortcomings of the prior art and provide a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination, which solves the problem of inaccurate GNSS satellite calculation and improves the accuracy of calculating GNSS satellites. It has important application value in the fields of GNSS occultation detection, GNSS reflection detection and other remote sensing detection such as low-orbit satellite payloads that have high accuracy requirements for GNSS satellite positions.

[0009] In order to achieve the above object, the present invention provides a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination, comprising the following steps:

[0010] S1, read GNSS satellite PCO information:

[0011] S1-1, download the GNSS satellite antenna phase center deviation file;

[0012] S1-2, select the corresponding phase center deviation file according to time;

[0013] S1-3, the phase center deviation file stores the antenna phase center deviation information at different frequencies of each GNSS satellite, and the storage method is based on the navigation system, frequency, and three-axis direction mode;

[0014] S2, calculate the position of the sun in the geocentric inertial system by analytical method:

[0015] S2-1, calculate TT time:

[0016] Formula 1: TT-TAI=32.184;

[0017] Formula 2: GPST-TAI=-19.0;

[0018] TAI is the world atomic time, in seconds; GPST is the GPS time, in seconds; TT time corresponding to the known GPST time is calculated by equation 1 and equation 2;

[0019] S2-2, calculate the position of the sun:

[0020] Convert TT time into the form of simplified Julian day, calculate the rotation angle θ at TT time, rotate the 3×3 matrix of all zeros about the X axis according to the rotation angle θ to obtain the matrix R, calculate the equatorial plane position vector W, the ecliptic coordinates, and rotate W to the equatorial coordinates through the matrix R, and then obtain the solar position S0, as shown in the following formula 3:

[0021] Formula 3: S0 = R·W;

[0022] S3, calculate the Z vector:

[0023] S3-1, using the Lagrange interpolation algorithm to obtain the current position P (x, y, z) vector of the GNSS navigation satellite; the P vector is a coordinate vector under the inertial coordinate system ECI, with the center of the earth as the coordinate origin;

[0024] S3-2, calculate the conversion matrix from the star-fixed system to the ECI, and the conversion matrix is ​​calculated in the celestial system; S3-3, calculate the Z matrix through Z=0~P;

[0025] S3-4, convert the Z matrix into a unit vector, the formula is as follows:

[0026] Formula 4:

[0027] S4, calculate the S vector:

[0028] Based on the known positions of S2-2 and GNSS satellites, the position vector S from the sun to the GNSS satellite is calculated, and the S vector is obtained by subtracting the two vectors. The formula is as follows:

[0029] Formula 5: S = S0-P;

[0030] According to formula 4, the S matrix is ​​converted into a unit vector;

[0031] S5, calculate the Y vector:

[0032] The Y vector represents the vertical vector of the plane formed by the Z vector and the S vector. The Y vector can be obtained by using the vector cross product calculation. The formula is as follows:

[0033] Formula 6: Y = Z × S;

[0034] According to formula 4, convert the Y matrix into a unit vector;

[0035] S6, calculate the X vector:

[0036] The X vector represents the vertical vector of the plane formed by the Y vector and the Z vector. The X vector can be obtained by using the vector cross product calculation. The formula is as follows:

[0037] Formula 7: X = Y × Z;

[0038] Convert the X matrix into a unit vector according to formula 4;

[0039] S7, calculate the rotation matrix S2C:

[0040] S2C is used to represent the conversion matrix from the star-fixed system to the ECI system. The formula is as follows:

[0041] Formula 8:

[0042] S8, read the PCO information of the corresponding satellite:

[0043] According to the current GNSS navigation satellite PRN number, find the PCO at the two frequency points of the corresponding satellite, which are represented by pco1(x,y,z) and pco2(x,y,z) respectively;

[0044] S9, rotate PCO to ECI coordinate system through S2C:

[0045] Use the S2C matrix to rotate pco1(x,y,z) and pco2(x,y,z) to the ECI coordinate system. The formula is as follows:

[0046] C1[0]=S2C[0]*pco1[0]+S2C[1]*pco1[1]+S2C[2]*pco1[2]

[0047] C1[1]=S2C[3]*pco1[0]+S2C[4]*pco1[1]+S2C[5]*pco1[2]

[0048] Formula 9: C1[2]=S2C[6]*pco1[0]+S2C[7]*pco1[1]+S2C[8]*pco1[2];

[0049] C2[0]=S2C[0]*pco2[0]+S2C[1]*pco2[1]+S2C[2]*pco2[2]

[0050] C2[1]=S2C[3]*pco2[0]+S2C[4]*pco2[1]+S2C[5]*pco2[2]

[0051] Formula 10: C2[2]=S2C[6]*pco2[0]+S2C[7]*pco2[1]+S2C[8]*pco2[2];

[0052] C1 and C2 are PCOs at two frequency points after rotation;

[0053] S10, dual-frequency combination calculation of the combined PCO of GNSS satellites:

[0054] The dual-frequency combination method is used to calculate an antenna phase center deviation C(x, y, z) to represent the point where the GNSS navigation star receives the signal. The formula is as follows:

[0055] Formula 11:

[0056] C[0]=(Freq_1*Freq_1*C1[0]-Freq_2*Freq_2*C2[0]) / (Freq_1*Freq_1-Freq_2*Freq_2)

[0057] C[1]=(Freq_1*Freq_1*C1[1]-Freq_2*Freq_2*C2[1]) / (Freq_1*Freq_1-Freq_2*Freq_2)

[0058] C[2]=(Freq_1*Freq_1*C1[2]-Freq_2*Freq_2*C2[2]) / (Freq_1*Freq_1-Freq_2*Freq_2);

[0059] Freq_1 and Freq_2 represent the frequencies of the two frequency points of GNSS navigation satellites, L1 and L2 represent the frequencies of GPS satellites, B1 and B3 represent the frequencies of BDS satellites, E1 and E5 represent the frequencies of GALILEO satellites, and G1 and G2 represent the frequencies of GLONASS satellites.

[0060] S11, correct the combined PCO to the GNSS satellite position:

[0061] By adding vectors, C is corrected to the GNSS navigation satellite position P to obtain the position P2 (x, y, z) after correcting the antenna phase center deviation. The formula is as follows:

[0062] P2[0]=P[0]+C[0]

[0063] P2[1]=P[1]+C[1]

[0064] P2[2]=P[2]+C[2].

[0065] Compared with the prior art, the present invention has the following beneficial effects:

[0066] The present invention improves the accuracy of calculating GNSS satellites, and has important application value in fields with high accuracy requirements for GNSS satellite positions in remote sensing detection such as GNSS occultation detection, GNSS reflection detection and other low-orbit satellite payloads. BRIEF DESCRIPTION OF THE DRAWINGS

[0067] Figure 1It is a schematic diagram of the technical route of the method of the present invention. DETAILED DESCRIPTION

[0068] The present invention will now be further described with reference to the accompanying drawings.

[0069] See also Figure 1 , a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination:

[0070] 1. Read GNSS satellite PCO information:

[0071] Download the GNSS satellite antenna phase center deviation file igs14.atx from the Internet, and select the corresponding file according to the time. Select igs14.atx before 344th of the year 2022, and select igs20.atx after that. This file stores the antenna phase center deviation information of each GNSS satellite at different frequencies, and stores this information according to the navigation system, frequency, and three-axis direction mode.

[0072] 2. Calculate the position of the sun:

[0073] The position of the sun is used in the PCO correction process of the GNSS navigation satellite, and the position of the sun needs to be calculated. The present invention uses an analytical method to calculate the position of the sun in the geocentric inertial system.

[0074] The sun's position is calculated in the geocentric time TT system using the following formula:

[0075] TT-TAI=32.184

[0076] GPST-TAI=-19.0

[0077] Calculate the TT time, where TAI is the world atomic time, the unit is seconds, and GPST is the GPS time, the unit is seconds. The TT time corresponding to the known GPST time can be calculated through two formulas.

[0078] Convert TT time into the form of simplified Julian day. Calculate the rotation angle θ at TT time, rotate the all-zero 3×3 matrix about the X axis according to the rotation angle θ to obtain the matrix R, calculate the equatorial position vector W, the ecliptic coordinates, and rotate W to the equatorial coordinates through the matrix R, and then obtain the sun position S0, as shown in the following formula:

[0079] S0=R·W

[0080] 3. Calculate the Z vector:

[0081] The Lagrange interpolation algorithm is used to obtain the current position P (x, y, z) vector of the GNSS navigation satellite (hereinafter referred to as P). The P vector is a coordinate vector in the inertial coordinate system ECI, with the center of the earth as the coordinate origin. The PCO information of the GNSS satellite is defined in the satellite-fixed system. The position of the GNSS navigation satellite is in the ECI coordinate system, so it is necessary to calculate the conversion matrix from the satellite-fixed system to the ECI. The conversion matrix is ​​calculated in the celestial system, so the Z matrix is ​​calculated by Z = 0-P. The Z matrix is ​​converted to a unit vector by the following formula:

[0082]

[0083] 4. Calculate the S vector:

[0084] The position of the sun to the GNSS satellite is represented by a vector S. Based on the above-mentioned known positions of the sun and the GNSS satellite, the position vector S from the sun to the GNSS satellite is calculated, and the S vector is obtained directly by subtracting the two vectors: S = S0-P. The S matrix is ​​converted into a unit vector by the formula in 3.

[0085] 5. Calculate the Y vector:

[0086] The Y vector represents the vertical vector of the plane formed by the Z vector and the S vector. The Y vector can be obtained by using the vector cross multiplication: Y = Z × S. And the Y matrix is ​​converted into a unit vector by the formula in 3.

[0087] 6. Calculate the X vector:

[0088] The X vector represents the vertical vector of the plane formed by the Y vector and the Z vector. The X vector can be obtained by using the vector cross multiplication: X = Y × Z. And the X matrix is ​​converted into a unit vector by the formula in 3.

[0089] 7. Calculate the rotation matrix S2C:

[0090] The antenna PCO of the GNSS navigation satellite needs to be rotated from the satellite-fixed system to the ECI inertial system in order to be corrected to the position of the GNSS navigation satellite. S2C is used to represent the conversion matrix from the satellite-fixed system to the ECI system. The calculation is shown in the following formula.

[0091]

[0092] 8. Read the PCO information of the corresponding satellite:

[0093] The GNSS satellite antenna PCO is stored according to the navigation system, frequency, and three-axis direction mode. The PCOs on the two frequencies of the corresponding satellites are represented by pco1(x, y, z) and pco2(x, y, z) respectively according to the current GNSS navigation satellite PRN number.

[0094] 9. Rotate PCO to ECI coordinate system through S2C:

[0095] The pco1(x, y, z) and pco2(x, y, z) are rotated to the ECI coordinate system through the S2C matrix, as shown in the following formula, where C1 and C2 are the PCOs at the two frequency points after rotation;

[0096] C1[0]=S2C[0]*pco1[0]+S2C[1]*pco1[1]+S2C[2]*pco1[2]

[0097] C1[1]=S2C[3]*pco1[0]+S2C[4]*pco1[1]+S2C[5]*pco1[2]

[0098] C1[2]=S2C[6]*pco1[0]+S2C[7]*pco1[1]+S2C[8]*pco1[2]

[0099] C2[0]=S2C[0]*pco2[0]+S2C[1]*pco2[1]+S2C[2]*pco2[2]

[0100] C2[1]=S2C[3]*pco2[0]+S2C[4]*pco2[1]+S2C[5]*pco2[2]

[0101] C2[2]=S2C[6]*pco2[0]+S2C[7]*pco2[1]+S2C[8]*pco2[2]

[0102] 10. Dual-frequency combination calculation of the combined PCO of GNSS satellites:

[0103] A GNSS position is usually used in GNSS remote sensing technology processing, but the rotated C1 and C2 on the two frequency points are calculated in i), and an antenna phase center deviation C(x, y, z) is calculated by dual-frequency combination to represent the point where the GNSS navigation star receives the signal. The specific calculation is shown in the following formula, where Freq_1 and Freq_2 represent the frequencies on the two frequency points of the GNSS navigation star, respectively. For GPS satellites, these two frequencies are L1 and L2, respectively, for BDS satellites, these two frequencies are B1 and B3, respectively, for GALILEO satellites, these two frequencies are E1 and E5, respectively, and for GLONASS satellites, these two frequencies are G1 and G2, respectively.

[0104] C[0]=(Freq_1*Freq_1*C1[0]-Freq_2*Freq_2*C2[0]) / (Freq_1*Freq_1-Freq_2*Freq_2)

[0105] C[1]=(Freq_1*Freq_1*C1[1]-Freq_2*Freq_2*C2[1]) / (Freq_1*Freq_1-Freq_2*Freq_2)

[0106] C[2]=(Freq_1*Freq_1*C1[2]-Freq_2*Freq_2*C2[2]) / (Freq_1*Freq_1-Freq_2*Freq_2);

[0107] 11. Correct the combined PCO to the GNSS satellite position:

[0108] By adding vectors, C is corrected to the GNSS navigation satellite position P to obtain the position P2 (x, y, z) after the antenna phase center deviation is corrected, as shown in the formula:

[0109] P2[0]=P[0]+C[0]

[0110] P2[1]=P[1]+C[1]

[0111] P2[2]=P[2]+C[2].

[0112] The above are only preferred embodiments of the present invention, which are only used to help understand the method and core ideas of the present application. The protection scope of the present invention is not limited to the above embodiments. All technical solutions under the concept of the present invention belong to the protection scope of the present invention. It should be pointed out that for ordinary technicians in this technical field, some improvements and modifications without departing from the principle of the present invention should also be regarded as the protection scope of the present invention.

[0113] The present invention solves the problem in the prior art that the orbit disclosed by the global GNSS service is the position at the satellite's center of mass, which is caused by the calculation deviation from the position of the receiving antenna, which further affects the excessive error in the inversion of the occultation event. By introducing a GNSS navigation satellite PCO correction method based on the sun position and dual-frequency combination, the accuracy of calculating GNSS satellites is improved, which has important application value in the fields of GNSS occultation detection, GNSS reflection detection and other remote sensing detections such as low-orbit satellite payloads that have high accuracy requirements for GNSS satellite positions.

Claims

1. A GNSS navigation satellite PCO correction method based on sun position and dual-frequency combination, characterized in that: The following steps are involved: S1, read GNSS satellite PCO information: S1-1, download the GNSS satellite antenna phase center deviation file; S1-2, select the corresponding phase center deviation file according to time; S1-3, the phase center deviation file stores antenna phase center deviation information at different frequencies of each GNSS satellite, and the storage method is based on the navigation system, frequency, and three-axis direction mode; S2, calculate the position of the sun in the geocentric inertial system by analytical method: S2-1, calculate TT time: Formula 1: TT-TAI=32.184; Formula 2: GPST-TAI=-19.0; The TAI is the world atomic time, in seconds; the GPST is the GPS time, in seconds; Calculate the TT time corresponding to the known GPST time by using the formula 1 and the formula 2; S2-2, calculate the position of the sun: The TT time is converted into the form of a simplified Julian day, the rotation angle θ at the TT time is calculated, the 3×3 matrix of all zeros is rotated about the X axis according to the rotation angle θ to obtain the matrix R, the equatorial plane position vector W and the ecliptic coordinates are calculated, and W is rotated to the equatorial coordinates through the matrix R, thereby obtaining the solar position S0, as shown in the following formula 3: Formula 3: S0 = R·W; S3, calculate the Z vector: S3-1, using the Lagrange interpolation algorithm to obtain the position P (x, y, z) vector of the GNSS navigation satellite at the current moment; the P vector is a coordinate vector under the inertial coordinate system ECI, with the center of the earth as the coordinate origin; S3-2, calculating the conversion matrix from the star-fixed system to the ECI, wherein the conversion matrix is ​​calculated in the celestial system; S3-3, calculate the Z matrix by Z=0-P; S3-4, convert the Z matrix into a unit vector, the formula is as follows: Formula 4: S4, calculate the S vector: On the basis of the S2-2 and the known position of the GNSS satellite, the position vector S from the sun to the GNSS satellite is calculated, and the S vector is obtained by subtracting the two vectors. The formula is as follows: Formula 5: S = S0-P; According to the formula 4, the S matrix is ​​converted into a unit vector; S5, calculate the Y vector: The Y vector represents the vertical vector of the plane formed by the Z vector and the S vector. The Y vector can be obtained by using the vector cross product calculation. The formula is as follows: Formula 6: Y = Z × S; According to the formula 4, the Y matrix is ​​converted into a unit vector; S6, calculate the X vector: The X vector represents the vertical vector of the plane formed by the Y vector and the Z vector. The X vector can be obtained by using the vector cross product calculation. The formula is as follows: Formula 7: X = Y × Z; According to the formula 4, the X matrix is ​​converted into a unit vector; S7, calculate the rotation matrix S2C: S2C is used to represent the conversion matrix from the star-fixed system to the ECI system. The formula is as follows: Formula 8: S8, read the PCO information of the corresponding satellite: According to the current GNSS navigation satellite PRN number, find the PCO at the two frequency points of the corresponding satellite, which are represented by pco1(x,y,z) and pco2(x,y,z) respectively; S9, rotate PCO to ECI coordinate system through S2C: Use the S2C matrix to rotate pco1(x,y,z) and pco2(x,y,z) to the ECI coordinate system. The formula is as follows: C1[0]=S2C[0]*pco1[0]+S2C[1]*pco1[1]+S2C[2]*pco1[2] C1[1]=S2C[3]*pco1[0]+S2C[4]*pco1[1]+S2C[5]*pco1[2] Formula 9: C1[2]=S2C[6]*pco1[0]+S2C[7]*pco1[1]+S2C[8]*pco1[2]; C2[0]=S2C[0]*pco2[0]+S2C[1]*pco2[1]+S2C[2]*pco2[2] C2[1]=S2C[3]*pco2[0]+S2C[4]*pco2[1]+S2C[5]*pco2[2] Formula 10: C2[2]=S2C[6]*pco2[0]+S2C[7]*pco2[1]+S2C[8]*pco2[2]; The C1 and C2 are PCOs at two frequency points after rotation respectively; S10, dual-frequency combination calculation of the combined PCO of GNSS satellites: The dual-frequency combination method is used to calculate an antenna phase center deviation C(x, y, z) to represent the point where the GNSS satellite receives the signal. The formula is as follows: Formula 11: C[0]=(Freq_1*Freq_1*C1[0]-Freq_2*Freq_2*C2[0]) / (Freq_1*Freq_1-Freq_2*Freq_2) C[1]=(Freq_1*Freq_1*C1[1]-Freq_2*Freq_2*C2[1]) / (Freq_1*Freq_1-Freq_2*Freq_2) C[2]=(Freq_1*Freq_1*C1[2]-Freq_2*Freq_2*C2[2]) / (Freq_1*Freq_1-Freq_2*Freq_2); Freq_1 and Freq_2 represent the frequencies of the two frequency points of GNSS satellites. The two frequencies of GPS satellites are L1 and L2, the two frequencies of BDS satellites are B1 and B3, the two frequencies of GALILEO satellites are E1 and E5, and the two frequencies of GLONASS satellites are G1 and G2. S11, correct the combined PCO to the GNSS satellite position: By adding vectors, C is corrected to the GNSS navigation satellite position P to obtain the position P2 (x, y, z) after correcting the antenna phase center deviation. The formula is as follows: P2[0]=P[0]+C[0] P2[1]=P[1]+C[1] P2[2]=P[2]+C[2].

Images