Ephemeris orbit calculation method suitable for orbit extrapolation on elliptical orbit satellite

Abstract

The invention relates to an ephemeris orbit calculation method suitable for orbit extrapolation on an elliptical orbit satellite. The ephemeris orbit calculation method comprises the following steps: calculating ephemeris orbit base point parameters on the ground; ephemeris orbit base point parameters which are injected from the ground are provided; and performing orbit extrapolation by using the ephemeris orbit base point parameters to calculate the position and the speed of the satellite in the inertial coordinate system. According to the method, orbit extrapolation is carried out on the basis of ephemeris orbit base point parameters injected on the ground, the calculation amount is low, the method is suitable for satellite-borne calculation, the calculation precision is high, the high-precision requirement of real-time orbit forecasting in the elliptic orbit satellite inter-satellite communication period can be met, the calculation capacity constraint of a satellite service system can be met, and the completion degree of an inter-satellite communication task is improved.

Description

Technical Field

[0001] The present invention relates to the field of orbit calculation technology, and in particular to an ephemeris orbit calculation method suitable for orbit extrapolation on an elliptical orbit star. Background Art

[0002] Currently, the real-time orbit data of target satellites can only be obtained by extrapolating the data onboard a satellite using a ground-based base point. Common methods for calculating satellite orbits onboard include analytical and numerical methods. Analytical methods offer simplicity and stable extrapolation, making them suitable for onboard calculations. However, their accuracy is limited and they cannot meet the high-precision requirements of real-time orbits for laser link establishment missions. While numerical integration (numerical methods) can meet these accuracy requirements, their frequent and complex calculations make them unsuitable for onboard calculations.

[0003] Currently, high-precision orbit calculations for target satellites often rely on ephemeris orbit extrapolation. Spaceborne ephemeris orbit extrapolation algorithms already support intersatellite laser link establishment missions for multiple circular orbit constellations. However, direct application of this algorithm to elliptical orbits results in significant errors and requires optimization. Therefore, a suitable ephemeris orbit calculation method for elliptical orbits is urgently needed. Summary of the Invention

[0004] To solve at least some of the above problems in the prior art, the present invention provides an ephemeris orbit calculation method applicable to orbit extrapolation on an elliptical orbit satellite, comprising:

[0005] Provide the ephemeris orbit base point parameters recorded by the ground; and

[0006] The ephemeris orbit base point parameters are used to perform orbit extrapolation to calculate the position and velocity of the satellite in the inertial coordinate system.

[0007] Calculate the fitting initial value σ of the ephemeris orbit base point parameters;

[0008] Calculate the position and velocity data, partial derivative matrix H and error vector L of each epoch in the fitting arc segment using the fitting initial values ​​of the orbital base point parameters of each ephemeris; and

[0009] Calculate the change Δσ of the ephemeris orbit base point parameter and add the change Δσ to the fitting initial value σ of each ephemeris orbit base point parameter to obtain the updated ephemeris orbit base point parameter σ 1 The update calculation is stopped after the loop is repeated until the iterative error vector is less than the target value or the absolute value of the change Δσ is less than the change threshold or the number of iterations reaches the iteration limit.

[0010] Furthermore, the ephemeris orbit base point parameters include ephemeris reference time, ephemeris fitting orbit time length, orbit semi-major axis at the reference time, orbit inclination at the reference time, mean anomaly angle at the reference time, ascending node longitude at the reference time, eccentricity, perigee argument, orbit inclination change rate, ascending node longitude change rate, difference between satellite average angular velocity and calculated value, amplitude of cosine harmonic correction term of latitude argument, amplitude of sine harmonic correction term of latitude argument, amplitude of cosine harmonic correction term of orbit inclination, amplitude of sine harmonic correction term of orbit inclination, amplitude of cosine harmonic correction term of orbit radius, amplitude of sine harmonic correction term of orbit radius, semi-major axis change rate and change rate of difference between satellite average angular velocity and calculated value.

[0011] Furthermore, calculating the position of the satellite using the ephemeris orbit base point parameters includes:

[0012] Calculate the coordinates of the satellite in the orbital plane;

[0013] Calculate the longitude of the ascending node; and

[0014] Calculate the coordinates of the satellite in the J2000 inertial coordinate system.

[0015] Furthermore, calculating the coordinates of the satellite in the orbital plane includes:

[0016] Calculate the time difference t between the onboard time and the reference time k :

[0017] t k =tt oe , where t oe is the ephemeris reference time (integrated seconds), t is the on-board time in integrated seconds at the calculation moment;

[0018] Calculate the semi-major axis A of the satellite orbit k :

[0019]

[0020] in, is the rate of change of the semi-major axis, A0 is the orbital semi-major axis at the reference time;

[0021] Calculate the satellite's average angular velocity n0 at the reference time:

[0022]

[0023] where μ = 3.986004418*10 14 ;

[0024] Calculate the deviation Δn of the satellite average angular rate based on the satellite average angular rate n0 at the reference time A :

[0025]

[0026] Where Δn0 is the difference between the satellite’s average angular velocity and the calculated value, is the rate of change of the difference between the satellite's average angular velocity and the calculated value;

[0027] Calculate the corrected satellite average angular velocity n based on the satellite average angular velocity n0 at the reference time A :

[0028] n A =n0+Δn A ;

[0029] Calculate the mean anomaly M based on the mean anomaly M0 at the reference time k :

[0030] M k =M0+n A ·t k ;

[0031] Based on the mean anomaly M k Iterative calculation of the eccentric anomaly angle E k :

[0032] E k =M k +coef_E_1·sinM k +coef_E_2·sin2M k +coef_E_3·sin3M k +coef_E_4·sin4M k

[0033] +coef_E_5·sin5M k +coef_E_6·sin6M k +coef_E_7·sin7M k +coef_E_8·sin8M k

[0034] +coef_E_9·sin9M k +coef_E_10·sin10M k

[0035] Among them, the calculation formulas of coef_E_1 to coef_E_10 are as follows:

[0036]

[0037] Where e is the eccentricity;

[0038] Based on the mean anomaly M k Calculate the true anomaly υk :

[0039] υ k =M k +coef_f_1·sinM k +coef_f_2·sin2M k +coef_f_3·sin3M k +coef_f_4·sin4M k

[0040] +coef_f_5·sin5M k +coef_f_6·sin6M k +coef_f_7·sin7M k +coef_f_8·sin8M k

[0041] +coef_f_9·sin9M k +coef_f_10·sin10M k

[0042] Among them, the calculation formulas of coef_f_1 to coef_f_10 are as follows:

[0043]

[0044] Based on the true anomaly υ k Calculate the latitude argument φ k :

[0045] φ k =υ k +ω,

[0046] where ω is the argument of perigee;

[0047] Calculate the latitude argument correction term δu k :

[0048] δu k =C us sin(2φ k )+C uc cos(2φ k )

[0049] Among them, C us is the amplitude of the sine harmonic correction term for the argument of latitude, C uc is the amplitude of the cosine harmonic correction for the argument of latitude.

[0050] Calculate the radial distance correction term δr k :

[0051] δr k=C rs sin(2φ k )+C rc cos(2φ k )

[0052] Among them C rs is the amplitude of the sinusoidal harmonic correction term for the orbital radius, C rc is the amplitude of the cosine harmonic correction term of the orbital radius;

[0053] Calculate the orbital inclination correction term δi k :

[0054] δi k =C is sin(2φ k )+C ic cos(2φ k ),

[0055] Among them C is is the amplitude of the sinusoidal harmonic correction term for orbital inclination, C ic is the amplitude of the cosine harmonic correction term for orbital inclination;

[0056] Using the latitude argument correction term δu k Calculate the corrected latitude argument u from the latitude argument φ k :

[0057] u k =φ k +δu k ;

[0058] Using the radial distance correction term δr k 、Semi-major axis of satellite orbit A k and the eccentric anomaly E k Calculate the corrected radial distance r k :

[0059] r k =A k ·(1-e·cosE k )+δr k ;

[0060] Based on the orbit inclination correction term δi k Calculate the corrected orbital inclination i k :

[0061]

[0062] where i0 is the orbital inclination at the reference time, is the rate of change of orbital inclination;

[0063] Based on the corrected radial distance rk and the corrected latitude argument u k Compute the coordinates of the satellite in the orbital plane:

[0064]

[0065] Furthermore, calculating the longitude of the ascending node includes:

[0066] Calculate the ascending node longitude Ω based on the ascending node longitude and the ascending node longitude change rate at the reference time k :

[0067]

[0068] Where Ω0 is the longitude of the ascending node at the reference time, is the rate of change of longitude of ascending node;

[0069] Calculating the satellite's coordinates in the J2000 inertial coordinate system includes:

[0070] Based on the longitude of the ascending node Ω k , the corrected orbital inclination i k The coordinates of the satellite in the J2000 inertial coordinate system are calculated using the coordinates of the satellite in the orbital plane:

[0071]

[0072] Furthermore, calculating the satellite's velocity using the ephemeris orbit base point parameters includes:

[0073] Based on the corrected satellite average angular velocity n A 、Eccentric anomaly angle E k The difference between the satellite's average angular velocity and the calculated value Δn0 is used to calculate the rate of change of the approach angle.

[0074]

[0075] Based on the eccentric anomaly angle E k , rate of change of approach angle and eccentricity e to calculate the rate of change of the angular distance of the ascending node

[0076]

[0077] Latitude argument φ k , the rate of change of the angular distance of the ascending node and the amplitude C of the sine harmonic correction term of the latitude argument us Calculate the rate of change of the latitude argument correction term

[0078]

[0079] Based on the latitude argument φ k , the rate of change of the angular distance of the ascending node and the amplitude C of the sinusoidal harmonic correction term of the orbital radius rs Calculate the radial correction distance change rate

[0080] Based on the latitude argument φ k , orbital inclination change rate and the amplitude C of the sinusoidal harmonic correction term of the orbital inclination is Calculate the satellite orbit inclination change rate

[0081]

[0082] Based on the semi-major axis change rate Eccentric anomaly E k , rate of change of approach angle Eccentricity e and satellite orbit semi-major axis A k Calculate the average angular velocity of the satellite

[0083]

[0084] Based on the rate of change of longitude of the ascending node Calculate the rate of change of the longitude of the ascending node while on board the star

[0085]

[0086] Based on the average angular velocity of the satellite Corrected radial distance r k and the corrected latitude argument u k Calculate the rate of change of the satellite within the orbital plane:

[0087]

[0088] The satellite's velocity in the J2000 inertial coordinate system is calculated based on the rate of change of the satellite in the orbital plane:

[0089]

[0090] Furthermore, the calculation of the fitting initial value σ of the ephemeris orbit base point parameter includes:

[0091] According to the time information of the ephemeris orbit base point parameters, the position and velocity of the satellite in the J2000 inertial coordinate system with the base point time as the starting point and the time x hours after the base point time as the end point is selected as the fitting arc segment;

[0092] Use the first epoch t of the fitted arc segment oe The position vector and velocity vector of the satellite in the J2000 system at this moment Calculate the corresponding t oe The six elements of Kepler's orbit at the moment;

[0093] t oe The six Kepler orbit elements at the time are assigned to the six parameters corresponding to the ephemeris orbit base point. The initial values ​​of the remaining parameters in the ephemeris orbit base point parameters are directly assigned to 0 to obtain the fitting initial values ​​of the ephemeris orbit base point parameters.

[0094] t oe The six Kepler orbital parameters at the time include the semi-major axis a0, eccentricity e0, inclination i0, longitude of ascending node Ω0, argument of perigee ω0, and mean anomaly M0. The specific calculations are as follows:

[0095]

[0096] Where R0 is the modulus of the position vector, h 0x h 0y h 0z are the components of angular momentum in the x, y, and z axes respectively, u0 is the argument of latitude, E0 is the eccentric anomaly, and f0 is the true anomaly;

[0097]

[0098] in, is the position vector, is the angular momentum, H0 is the modulus of the angular momentum, r 0x r 0y r 0z are the components of the position vector on the x, y, and z axes, respectively, and h 0z is the component of angular momentum on the z-axis, and the function atan2() is the four-quadrant inverse tangent function.

[0099] Furthermore, based on the position vector and velocity vector of the satellite corresponding to epoch j in the fitting arc segment data Calculate the partial derivative matrix H of the satellite position with respect to the ephemeris orbit base point parameters k ,include:

[0100] The partial differential derivative of the satellite position with respect to the orbital semi-major axis A0 at the reference time is:

[0101]

[0102] The semi-major axis The time difference between the onboard time and the reference time t k =tt oe ;

[0103] Partial derivative of satellite position with respect to eccentricity e:

[0104]

[0105] Among them E k is the eccentric anomaly angle, r k is the corrected radial distance, n A is the corrected average angular velocity of the satellite;

[0106] The partial derivative of the satellite position with respect to the orbital inclination i0 at the reference time is:

[0107]

[0108] where Ω k is the longitude of the ascending node, i k is the corrected orbital inclination, u k is the corrected latitude argument;

[0109] Partial derivative of satellite position with respect to parameter Ω0:

[0110]

[0111] Partial derivative of the satellite position with respect to the argument of perigee ω:

[0112]

[0113] The partial differential derivative of the satellite position with respect to the mean anomaly M0 at the reference time is:

[0114]

[0115] Satellite position change rate of semi-major axis Partial derivative of :

[0116]

[0117] where t k is the time difference between the onboard time and the reference time;

[0118] The partial derivative of the satellite position with respect to the satellite average angular rate Δn0 at the reference time is:

[0119]

[0120] The rate of change of the satellite position relative to the difference between the satellite's average angular velocity and the calculated value Partial derivative of :

[0121]

[0122] Rate of change of satellite position with respect to the longitude of ascending node Partial derivative of :

[0123]

[0124] Where Ω0 is the longitude of the ascending node at the reference time;

[0125] Satellite position versus orbital inclination rate of change Partial derivative of :

[0126]

[0127] The amplitude C of the sinusoidal harmonic correction of satellite position to latitude argument us Partial derivative of :

[0128]

[0129] The amplitude C of the cosine harmonic correction of the satellite position to the latitude argument uc Partial derivative of :

[0130]

[0131] The amplitude C of the sinusoidal harmonic correction of satellite position to orbit inclination is Partial derivative of :

[0132]

[0133] The amplitude C of the cosine harmonic correction of the satellite position to the orbital inclination ic Partial derivative of :

[0134]

[0135] The amplitude C of the sinusoidal harmonic correction of the satellite position to the parametric orbit radius rs Partial derivative of :

[0136]

[0137] The amplitude C of the cosine harmonic correction of the satellite position to the orbit radius rc Partial derivative of :

[0138]

[0139] Calculate the partial derivative matrix H of the satellite position with respect to the ephemeris orbit base point parameters k :

[0140]

[0141] Furthermore, the partial derivative matrix H of the satellite position to the ephemeris orbit base point parameters is k Iteratively calculate the ephemeris orbit base point parameters, including:

[0142] Calculate the satellite position calculation value of each epoch under the current fitting initial value conditions Where j is the epoch identifier and k is the iteration number identifier;

[0143] The iterative calculation equation is established as follows:

[0144] H k ·Δσ k =L k

[0145] In the above formula, H k Calculate the partial derivative matrix for k iterations, L k is the difference between the kth ephemeris calculation value and the fitting arc segment, Δσ k is the change of each parameter of the ephemeris orbit base point calculated in the kth iteration;

[0146] Next, we use Givens transform to calculate Δσ k Valuation, the formula is as follows:

[0147]

[0148] In the above formula, qr(H k ) represents the QR decomposition of the matrix in the brackets, and the subscripts in the above formula represent the elements at the corresponding positions in the matrix;

[0149] The ephemeris orbit base point parameters after k iterations are calculated as follows:

[0150] σ k+1 =σ k +Δσ k ;

[0151] Calculating the error vector involves:

[0152] Error matrix L k The error vector at each epoch in is calculated according to the following formula:

[0153]

[0154] In the above formula, is the position vector of the satellite at the jth epoch in the fitting arc; The satellite position vector at the jth epoch calculated based on the ephemeris orbit base point parameters calculated based on the k-th iteration results;

[0155] The conditions for the end of iterative calculation are as follows:

[0156] When the error matrix L k The mean square error is less than the target value, or the absolute value of the change is less than the change threshold |Δσk When |≤ε, or the number of iterations reaches the iteration limit, the update calculation is stopped.

[0157] The present invention has at least the following beneficial effects: the present invention provides an ephemeris orbit calculation method suitable for orbit extrapolation on an elliptical orbit satellite, which performs orbit extrapolation based on the ephemeris orbit base point parameters noted on the ground. The method has a low computational complexity and is suitable for onboard calculations, and has high computational accuracy. The method can not only meet the high-precision requirements of real-time orbit prediction during inter-satellite communication of elliptical orbit satellites, but also meet the computing power constraints of the satellite service system, thereby improving the completion of inter-satellite communication tasks. DETAILED DESCRIPTION

[0158] In the present invention, each embodiment is only intended to illustrate the aspects of the present invention and should not be construed as limiting.

[0159] In the present invention, unless otherwise specified, the quantifiers "a" and "an" do not exclude the presence of multiple elements.

[0160] It should also be pointed out that in the embodiments of the present invention, for the sake of clarity and simplicity, only a portion of the parts or components may be shown, but a person skilled in the art will understand that under the teachings of the present invention, the required parts or components may be added according to the needs of the specific scenario.

[0161] It should also be pointed out that within the scope of the present invention, the terms "same", "equal", "equal to" and the like do not mean that the two values ​​are absolutely equal, but allow a certain reasonable error, that is, the terms also cover "substantially the same", "substantially equal", and "substantially equal to".

[0162] In addition, the terms "first" and "second" are used for descriptive purposes only and should not be understood as indicating or implying relative importance.

[0163] In addition, the embodiments of the present invention describe the process steps in a specific order, but this is only for the convenience of distinguishing the steps, and does not limit the order of the steps. In different embodiments of the present invention, the order of the steps can be adjusted according to the adjustment of the process.

[0164] The existing inter-satellite communication has a high-precision demand for real-time orbit prediction of the target satellite, which is usually achieved by performing on-orbit real-time orbit prediction based on ephemeris orbit parameters, but the current ephemeris orbit prediction algorithm is only applicable to circular orbits. In response to the high-precision real-time orbit prediction accuracy demand for the target satellite in elliptical orbit inter-satellite communication, the present invention proposes an ephemeris orbit calculation method suitable for orbit prediction on elliptical orbit satellites, so that it can not only meet the high-precision requirements of real-time orbit prediction during inter-satellite communication of elliptical orbit satellites, but also meet the computing power constraints of the satellite service system, thereby improving the completion of inter-satellite communication tasks. The method of the present invention fills the technical gap in the aerospace field for high-precision orbit extrapolation algorithms on elliptical orbit satellites, has a corresponding calculation method, increases the universality of the ephemeris orbit extrapolation algorithm to satellite orbit types, and improves the completion of inter-satellite communication tasks.

[0165] A method for calculating an ephemeris orbit for orbit extrapolation on an elliptical orbit satellite includes:

[0166] Provide the ephemeris orbit base point parameters recorded by the ground; and

[0167] The orbital extrapolation is performed using the ephemeris orbital base point parameters to calculate the position and velocity of the satellite in the inertial coordinate system.

[0168] Specifically, the ephemeris orbit base point parameters are calculated on the ground and then uploaded to the satellite. The satellite's onboard computer calculates the position and velocity of the satellite in the inertial coordinate system (J2000.0 coordinate system) at any time based on the ephemeris orbit base point parameters.

[0169] Table 1 Ephemeris orbit base point parameters

[0170]

[0171]

[0172] The angles in the ephemeris orbit base point parameters are converted to radians before use. The units in Table 1 are π, π / s, and π / s. 2 The parameters are all angle-related quantities.

[0173] The specific process of using the ephemeris orbit base point parameters recorded by the ground station to perform orbit extrapolation to calculate the position and velocity of the satellite in the inertial coordinate system is as follows:

[0174] Step 1, using the ephemeris orbit base point parameters to calculate the satellite position includes:

[0175] Step 1.1, calculate the time difference t between the satellite time and the reference time k :

[0176] t k =tt oe ,

[0177] where t oe is the ephemeris reference time (integrated seconds), and t is the on-board time in integrated seconds at the calculation moment.

[0178] Step 1.2, calculate the semi-major axis A of the satellite orbit k :

[0179]

[0180] in, is the rate of change of the semi-major axis, and A0 is the orbital semi-major axis at the reference time.

[0181] Step 1.3, calculate the average satellite angular rate n0 at the reference time:

[0182]

[0183] where μ = 3.986004418*10 14 .

[0184] Step 1.4: Calculate the deviation Δn of the satellite average angular rate based on the satellite average angular rate n0 at the reference time. A :

[0185]

[0186] Where Δn0 is the difference between the satellite’s average angular velocity and the calculated value, is the rate of change of the difference between the satellite's average angular velocity and the calculated value.

[0187] Step 1.5: Calculate the corrected satellite average angular velocity n based on the satellite average angular velocity n0 at the reference time. A :

[0188] n A =n0+Δn A .

[0189] Step 1.6, calculate the mean anomaly M based on the mean anomaly M0 at the reference time k :

[0190] M k =M0+n A ·t k .

[0191] Step 1.7, based on the mean anomaly M k Iterative calculation of the eccentric anomaly angle E k :E k =M k +coef_E_1·sinM k +coef_E_2·sin2M k+coef_E_3·sin3M k +coef_E_4·sin4M k

[0192] +coef_E_5·sin5M k +coef_E_6·sin6M k +coef_E_7·sin7M k +coef_E_8·sin8M k

[0193] +coef_E_9·sin9M k +coef_E_10·sin10M k

[0194] Among them, the calculation formulas of coef_E_1 to coef_E_10 are as follows:

[0195]

[0196] Where e is the eccentricity.

[0197] Step 1.8, based on the mean anomaly M k Calculate the true anomaly υ k :

[0198] υ k =M k +coef_f_1·sinM k +coef_f_2·sin2M k +coef_f_3·sin3M k +coef_f_4·sin4M k

[0199] +coef_f_5·sin5M k +coef_f_6·sin6M k +coef_f_7·sin7M k +coef_f_8·sin8M k

[0200] +coef_f_9·sin9M k +coef_f_10·sin10M k

[0201] Among them, the calculation formulas of coef_f_1 to coef_f_10 are as follows:

[0202]

[0203] Step 1.9, based on the true anomaly υk Calculate the latitude argument φ k :

[0204] φ k =υ k +ω,

[0205] where ω is the argument of perigee.

[0206] Step 1.10, calculate the latitude argument correction term δu k :

[0207] δu k =C us sin(2φ k )+C uc cos(2φ k ),

[0208] Among them, C us is the amplitude of the sine harmonic correction term for the argument of latitude, C uc is the amplitude of the cosine harmonic correction for the argument of latitude.

[0209] Step 1.11, calculate the radial distance correction term δr k :

[0210] δr k =C rs sin(2φ k )+C rc cos(2φ k ),

[0211] Among them C rs is the amplitude of the sinusoidal harmonic correction term for the orbital radius, C rc is the amplitude of the harmonic correction term for the cosine of the orbital radius.

[0212] Step 1.12, calculate the orbital inclination correction term δi k :

[0213] δi k =C is sin(2φ k )+C ic cos(2φ k ),

[0214] Among them C is is the amplitude of the sinusoidal harmonic correction term for orbital inclination, C ic is the amplitude of the cosine harmonic correction for orbital inclination.

[0215] Step 1.13, use the latitude argument correction term δu k and latitude argument φ k Calculate the corrected latitude argument uk :

[0216] u k =φ k +δu k .

[0217] Step 1.14, use the radial distance correction term δr k 、Semi-major axis of satellite orbit A k and the eccentric anomaly E k Calculate the corrected radial distance r k :

[0218] r k =A k ·(1-e·cosE k )+δr k .

[0219] Step 1.15, based on the orbital inclination correction term δi k Calculate the corrected orbital inclination i k :

[0220]

[0221] where i0 is the orbital inclination at the reference time, is the rate of change of orbital inclination.

[0222] Step 1.16, based on the corrected radial distance r k and the corrected latitude argument u k Compute the coordinates of the satellite in the orbital plane:

[0223]

[0224] Step 1.17, calculate the longitude of the ascending node Ω k :

[0225]

[0226] Where Ω0 is the longitude of the ascending node at the reference time, is the rate of change of longitude of the ascending node.

[0227] Step 1.18, based on the longitude of the ascending node Ω k , the corrected orbital inclination i k The coordinates of the satellite in the J2000 inertial coordinate system are calculated using the coordinates of the satellite in the orbital plane:

[0228]

[0229] Step 2, based on the above satellite position calculation method, using the ephemeris orbit base point parameters to calculate the satellite velocity includes:

[0230] Step 2.1, based on the corrected satellite average angular rate n A 、Eccentric angle E k The difference between the satellite's average angular velocity and the calculated value Δn0 is used to calculate the rate of change of the approach angle.

[0231] Step 2.2, based on the eccentric anomaly angle E k , rate of change of approach angle and eccentricity e to calculate the rate of change of the angular distance of the ascending node

[0232]

[0233] Step 2.3, latitude argument φ k , the rate of change of the angular distance of the ascending node and the amplitude C of the sine harmonic correction term of the latitude argument us Calculate the rate of change of the latitude argument correction term

[0234] Step 2.4, based on the latitude argument φ k , the rate of change of the angular distance of the ascending node and the amplitude C of the sinusoidal harmonic correction term of the orbital radius rs Calculate the radial correction distance change rate

[0235]

[0236] Step 2.5, based on the latitude argument φ k , orbital inclination change rate and the amplitude C of the sinusoidal harmonic correction term of the orbital inclination is Calculate the satellite orbit inclination change rate

[0237] Step 2.6, based on the semi-major axis change rate Eccentric anomaly E k , rate of change of approach angle Eccentricity e and satellite orbit semi-major axis A k Calculate the average angular velocity of the satellite

[0238] Step 2.7, based on the rate of change of longitude of the ascending node Calculate the rate of change of the longitude of the ascending node while on board the star

[0239]

[0240] Step 2.8, based on the average angular rate of the satellite Corrected radial distance r k and the corrected latitude argument u k Calculate the rate of change of the satellite within the orbital plane:

[0241]

[0242] Step 2.9, calculate the satellite's velocity in the J2000 inertial coordinate system based on the satellite's rate of change in the orbital plane:

[0243]

[0244] The ephemeris orbit base point parameters calculated at the ground station include:

[0245] The ephemeris orbit base point parameters are calculated through multiple fitting iterations. The initial fitting value σ of each ephemeris orbit base point parameter is first calculated. Then, the position and velocity data, partial derivative matrix H and error vector L of each epoch in the fitting arc are calculated using the initial fitting value of each ephemeris orbit base point parameter. The change Δσ of the ephemeris orbit base point parameter is further calculated. The change Δσ is added to the initial fitting value σ of each ephemeris orbit base point parameter to obtain the updated ephemeris orbit base point parameter σ. 1 , and the update calculation is stopped after the cycle repeats until the iteration accuracy meets the requirements or the number of iterations reaches the iteration limit.

[0246] The steps for calculating the ephemeris orbit base point parameters include: (1) calculating the initial fitting value; (2) calculating the partial derivative matrix; (3) iteratively calculating the ephemeris orbit base point parameters; (4) calculating the error vector; and (5) ending the iterative calculation.

[0247] The specific relevant parameters are calculated as follows:

[0248] Step 3.1, calculate the initial values ​​of the fitting of the ephemeris orbit base point parameters:

[0249] According to the time information of the ephemeris orbit base point parameters as needed, the position and velocity of the satellite in the J2000 inertial coordinate system are selected as the fitting arc segment, with the base point time as the starting point and the time x hours (the ephemeris fitting orbit time length) after the base point time as the end point.

[0250] Use the first epoch t of the fitted arc segment oe The position vector and velocity vector of the satellite in the J2000 system at this moment Calculate the corresponding t oeThe six Kepler orbital elements at the time (that is, the initial six Kepler orbital elements). The angles of the six Kepler orbital elements need to be converted to the range [0, 2π). The first epoch t oe The time is also the ephemeris reference time. oe The six Kepler orbital parameters at the moment include: orbital semi-major axis a0, eccentricity e0, orbital inclination i0, ascending node longitude Ω0, argument of perigee ω0, and mean anomaly M0. The specific calculation method is as follows:

[0251]

[0252] Where R0 is the modulus of the position vector, h 0x h 0y h 0z are the components of angular momentum in the x, y, and z axes, u0 is the argument of latitude, E0 is the eccentric anomaly, and f0 is the true anomaly.

[0253]

[0254] in, is the position vector, is the angular momentum, H0 is the modulus of the angular momentum, r 0x r 0y r 0z are the components of the position vector on the x, y, and z axes, respectively, and h 0z is the component of angular momentum on the z-axis, and the function atan2() is the four-quadrant inverse tangent function.

[0255] The t calculated above oe The six Kepler orbit numbers at the time are assigned to the six parameters corresponding to the ephemeris orbit base point. The initial values ​​of other parameters in the ephemeris orbit base point parameters are directly assigned to 0. The ephemeris orbit base point parameters at this time are the fitting initial values ​​of the ephemeris orbit base point parameters.

[0256] Step 3.2, calculate the partial derivative matrix. They are respectively the position vector and velocity vector of the satellite corresponding to epoch j in the fitting arc segment data, and the partial derivative matrix H of the satellite position with respect to the ephemeris orbit base point parameters k The specific calculation method is as follows. The definitions and calculation methods of the relevant parameters in the formula are consistent with the algorithm formula involved in the orbit extrapolation process using the ephemeris orbit base point parameters annotated by the ground station.

[0257] Step 3.2.1, calculate the partial differential of the satellite position with respect to the semi-major axis A0 of the orbit at the reference time:

[0258]

[0259] The semi-major axis The time difference between the onboard time and the reference time t k =tt oe .

[0260] Step 3.2.2, calculate the partial differential of the satellite position with respect to the eccentricity e:

[0261]

[0262] Among them E k is the eccentric anomaly angle, r k is the corrected radial distance, n A is the corrected average angular velocity of the satellite.

[0263] Step 3.2.3, calculate the partial derivative of the satellite position with respect to the orbital inclination i0 at the reference time:

[0264]

[0265] where Ω k is the longitude of the ascending node, i k is the corrected orbital inclination, u k is the corrected latitude argument.

[0266] Step 3.2.4, calculate the partial differential of the satellite position with respect to the parameter Ω0:

[0267]

[0268] Step 3.2.5, calculate the partial derivative of the satellite position with respect to the argument of perigee angle ω:

[0269]

[0270] Step 3.2.6, calculate the partial differential of the satellite position with respect to the mean anomaly M0 at the reference time:

[0271]

[0272] Step 3.2.7, satellite position change rate of semi-major axis Partial derivative of :

[0273]

[0274] where t k It is the time difference between the onboard time and the reference time.

[0275] Step 3.2.8, calculate the partial derivative of the satellite position with respect to the satellite average angular rate Δn0 at the reference time:

[0276]

[0277] Step 3.2.9, the rate of change of the satellite position to the difference between the satellite's average angular velocity and the calculated value Partial derivative of :

[0278]

[0279] Step 3.2.10, the rate of change of satellite position to the longitude of the ascending node Partial derivative of :

[0280]

[0281] Where Ω0 is the longitude of the ascending node at the reference time.

[0282] Step 3.2.11, satellite position versus orbital inclination rate of change Partial derivative of :

[0283]

[0284] Step 3.2.12, the amplitude C of the sinusoidal harmonic correction of the satellite position to the latitude argument us Partial derivative of :

[0285]

[0286] Step 3.2.13, the amplitude C of the cosine harmonic correction of the satellite position to the latitude argument uc Partial derivative of :

[0287]

[0288] Step 3.2.14, the amplitude C of the sinusoidal harmonic correction of satellite position to orbit inclination is Partial derivative of :

[0289]

[0290] Step 3.2.15, the amplitude C of the cosine harmonic correction of the satellite position to the orbital inclination ic Partial derivative of :

[0291]

[0292] Step 3.2.16, the amplitude C of the sinusoidal harmonic correction of the satellite position to the parametric orbit radius rs Partial derivative of :

[0293]

[0294] Step 3.2.17, the amplitude C of the cosine harmonic correction of the satellite position to the orbit radius rc Partial derivative of :

[0295]

[0296] Steps 3.2.1 to 3.2.17 are the calculation process of the partial differential of the satellite position with respect to the parameters of each ephemeris orbit base point.

[0297] Step 3.2.18, calculate the partial derivative matrix H of the satellite position with respect to the ephemeris orbit base point parameters k :

[0298]

[0299] Step 3.3, the partial derivative matrix H of the satellite position to the ephemeris orbit base point parameters k Perform iterative calculation of ephemeris orbit base point parameters.

[0300] Based on the calculation formula involved in the orbit extrapolation based on the ephemeris orbit base point parameters, the satellite position calculation value of each epoch under the current initial value conditions can be calculated (j is the epoch identifier, k is the iteration number identifier).

[0301] The iterative calculation equation is established as follows:

[0302] H k ·Δσ k =L k

[0303] In the above formula, H k Calculate the partial derivative matrix for k iterations; L k is the difference between the kth ephemeris calculation value and the fitting arc segment; Δσ k The change in each parameter of the ephemeris orbit base point calculated in the kth iteration.

[0304] Next, we use Givens transform to calculate Δσ k Valuation, the formula is as follows:

[0305]

[0306] In the above formula, qr(H k ) represents the QR decomposition of the matrix in the brackets. The subscripts in the above formula represent the elements at the corresponding positions in the matrix. For example: R ij Represents the (i, j)th element in the R matrix.

[0307] The ephemeris orbit base point parameters after k iterations can be further calculated as:

[0308] σ k+1 =σ k +Δσ k .

[0309] Step 3.4, calculate the error vector

[0310] Error matrix L k The error vector at each epoch in is calculated as follows:

[0311]

[0312] In the above formula, is the position vector of the satellite at the jth epoch in the fitting arc; The satellite position vector at the jth epoch is calculated based on the ephemeris orbit base point parameters calculated based on the k iteration results.

[0313] The conditions for the end of iterative calculation are as follows:

[0314] When the error matrix L k The mean square error is less than the target value (a small amount, determined according to actual needs), or |Δσ k The iterative calculation ends when |≤ε (the change threshold ε is a small amount and is determined according to actual needs). If the above requirements are too strict and the accuracy cannot be converged to the required range, the iterative calculation can be terminated after multiple iterations (for example, 30 times).

[0315] Although certain embodiments of the present invention have been described in this application, those skilled in the art will appreciate that these embodiments are provided by way of example only. Numerous variations, alternatives, and modifications will be contemplated by those skilled in the art in light of the teachings of this disclosure without departing from the scope of the present invention. The appended claims are intended to define the scope of the present invention and are intended to encompass methods and structures within the scope of these claims and their equivalents.

Claims

1. A method for calculating ephemeris orbits for orbit extrapolation on elliptical orbit satellites, characterized in that: include: Provide the ephemeris orbit base point parameters recorded on the ground; as well as The ephemeris orbit base point parameters are used to perform orbit extrapolation to calculate the position and velocity of the satellite in the inertial coordinate system.

2. The method according to claim 1, characterized in that It also includes calculating the ephemeris orbit base point parameters on the ground, which include: Calculate the fitting initial value σ of the ephemeris orbit base point parameters; Calculate the position and velocity data, partial derivative matrix H and error vector L of each epoch in the fitting arc segment using the fitting initial values ​​of the orbital base point parameters of each ephemeris; and Calculate the change Δσ of the ephemeris orbit base point parameter and add the change Δσ to the fitting initial value σ of each ephemeris orbit base point parameter to obtain the updated ephemeris orbit base point parameter σ 1 The update calculation is stopped after the loop is repeated until the iterative error vector is less than the target value or the absolute value of the change Δσ is less than the change threshold or the number of iterations reaches the iteration limit.

3. The method according to claim 2, characterized in that The ephemeris orbit base point parameters include the ephemeris reference time, the ephemeris fitting orbit time length, the orbit semi-major axis at the reference time, the orbit inclination at the reference time, the mean anomaly angle at the reference time, the ascending node longitude at the reference time, the eccentricity, the perigee argument, the orbit inclination change rate, the ascending node longitude change rate, the difference between the satellite average angular velocity and the calculated value, the amplitude of the cosine harmonic correction term of the latitude argument, the amplitude of the sine harmonic correction term of the latitude argument, the amplitude of the cosine harmonic correction term of the orbit inclination, the amplitude of the sine harmonic correction term of the orbit inclination, the amplitude of the cosine harmonic correction term of the orbit radius, the amplitude of the sine harmonic correction term of the orbit radius, the semi-major axis change rate, and the change rate of the difference between the satellite average angular velocity and the calculated value.

4. The method according to claim 3, characterized in that Calculating the satellite position using ephemeris orbit base point parameters includes: Calculate the coordinates of the satellite in the orbital plane; Calculate the longitude of the ascending node; and Calculate the coordinates of the satellite in the J2000 inertial coordinate system.

5. The method according to claim 4, characterized in that Calculating the coordinates of a satellite in the orbital plane involves: Calculate the time difference t between the onboard time and the reference time k : t k =tt oe , where t oe is the ephemeris reference time (integrated seconds), t is the on-board time in integrated seconds at the calculation moment; Calculate the semi-major axis A of the satellite orbit k : in, is the rate of change of the semi-major axis, A0 is the orbital semi-major axis at the reference time; Calculate the satellite's average angular velocity n0 at the reference time: where μ = 3.986004418*10 14 ; Calculate the deviation Δn of the satellite average angular rate based on the satellite average angular rate n0 at the reference time A : Where Δn0 is the difference between the satellite’s average angular velocity and the calculated value, is the rate of change of the difference between the satellite's average angular velocity and the calculated value; Calculate the corrected satellite average angular velocity n based on the satellite average angular velocity n0 at the reference time A : n A =n0+Δn A ; Calculate the mean anomaly M based on the mean anomaly M0 at the reference time k : M k =M0+n A ·t k ; Based on the mean anomaly M k Iterative calculation of the eccentric anomaly angle E k : AND k =M k +coef_E_1·sinM k +coef_E_2·sin2M k +coef_E_3·sin3M k +coef_E_4·sin4M k +coef_E_5·sin5M k +coef_E_6·sin6M k +coef_E_7·sin7M k +coef_E_8·sin8M k +coef_E_9·sin9M k +coef_E_10·sin10M k Among them, the calculation formulas of coef_E_1 to coef_E_10 are as follows: Where e is the eccentricity; Based on the mean anomaly M k Calculate the true anomaly υ k : υ k =M k +coef_f_1·sinM k +coef_f_2·sin2M k +coef_f_3·sin3M k +coef_f_4·sin4M k +coef_f_5·sin5M k +coef_f_6·sin6M k +coef_f_7·sin7M k +coef_f_8·sin8M k +coef_f_9·sin9M k +coef_f_10·sin10M k Among them, the calculation formulas of coef_f_1 to coef_f_10 are as follows: Based on the true anomaly υ k Calculate the latitude argument φ k : f k =υ k +oh, where ω is the argument of perigee; Calculate the latitude argument correction term δu k : you k =C us sin(2φ k )+C uc cos(2φ k ) Among them, C us is the amplitude of the sine harmonic correction term for the argument of latitude, C uc is the amplitude of the cosine harmonic correction for the argument of latitude. Calculate the radial distance correction term δr k : δr k =C rs sin(2φ k )+C rc cos(2φ k ) Among them C rs is the amplitude of the sinusoidal harmonic correction term for the orbital radius, C rc is the amplitude of the cosine harmonic correction term of the orbital radius; Calculate the orbital inclination correction term δi k : δi k =C is sin(2φ k )+C ic cos(2φ k ), Among them C is is the amplitude of the sinusoidal harmonic correction term for orbital inclination, C ic is the amplitude of the cosine harmonic correction term for orbital inclination; Using the latitude argument correction term δu k and latitude argument φ k Calculate the corrected latitude argument u k : you k =φ k +du k ; Using the radial distance correction term δr k 、Semi-major axis of satellite orbit A k and the eccentric anomaly E k Calculate the corrected radial distance r k : r k =A k ·(1-e·cosE k )+δr k ; Based on the orbit inclination correction term δi k Calculate the corrected orbital inclination i k : where i0 is the orbital inclination at the reference time, is the rate of change of orbital inclination; Based on the corrected radial distance r k and the corrected latitude argument u k Compute the coordinates of the satellite in the orbital plane:

6. The method according to claim 5, characterized in that Calculating the longitude of the ascending node involves: Calculate the ascending node longitude Ω based on the ascending node longitude and the ascending node longitude change rate at the reference time k : Where Ω0 is the longitude of the ascending node at the reference time, is the rate of change of longitude of ascending node; Calculating the satellite's coordinates in the J2000 inertial coordinate system includes: Based on the longitude of the ascending node Ω k , the corrected orbital inclination i k The coordinates of the satellite in the J2000 inertial coordinate system are calculated using the coordinates of the satellite in the orbital plane:

7. The method according to claim 6, characterized in that Calculating the satellite's velocity using ephemeris orbit base point parameters includes: Based on the corrected satellite average angular velocity n A 、Eccentric anomaly angle E k The difference between the satellite's average angular velocity and the calculated value Δn0 is used to calculate the rate of change of the approach angle. Based on the eccentric anomaly angle E k , rate of change of approach angle and eccentricity e to calculate the rate of change of the angular distance of the ascending node Latitude argument φ k , the rate of change of the angular distance of the ascending node and the amplitude C of the sine harmonic correction term of the latitude argument us Calculate the rate of change of the latitude argument correction term Based on the latitude argument φ k , the rate of change of the angular distance of the ascending node and the amplitude C of the sinusoidal harmonic correction term of the orbital radius rs Calculate the radial correction distance change rate Based on the latitude argument φ k , orbital inclination change rate and the amplitude C of the sinusoidal harmonic correction term of the orbital inclination is Calculate the satellite orbit inclination change rate Based on the semi-major axis change rate Eccentric anomaly E k , rate of change of approach angle Eccentricity e and satellite orbit semi-major axis A k Calculate the average angular velocity of the satellite Based on the rate of change of longitude of the ascending node Calculate the rate of change of the longitude of the ascending node while on board the star Based on the average angular velocity of the satellite Corrected radial distance r k and the corrected latitude argument u k Calculate the rate of change of the satellite within the orbital plane: The satellite's velocity in the J2000 inertial coordinate system is calculated based on the rate of change of the satellite in the orbital plane:

8. The method according to claim 7, characterized in that The fitting initial value σ for calculating the ephemeris orbit base point parameters includes: According to the time information of the ephemeris orbit base point parameters, the position and velocity of the satellite in the J2000 inertial coordinate system with the base point time as the starting point and the time x hours after the base point time as the end point is selected as the fitting arc segment; Use the first epoch t of the fitted arc segment oe The position vector and velocity vector of the satellite in the J2000 system at this moment Calculate the corresponding t oe The six elements of Kepler's orbit at the moment; t oe The six Kepler orbit elements at the time are assigned to the six parameters corresponding to the ephemeris orbit base point. The initial values ​​of the remaining parameters in the ephemeris orbit base point parameters are directly assigned to 0 to obtain the fitting initial values ​​of the ephemeris orbit base point parameters. t oe The six Kepler orbital parameters at the time include the semi-major axis a0, eccentricity e0, inclination i0, longitude of ascending node Ω0, argument of perigee ω0, and mean anomaly M0. The specific calculations are as follows: Where R0 is the modulus of the position vector, h 0x h 0y h 0z are the components of angular momentum in the x, y, and z axes respectively, u0 is the argument of latitude, E0 is the eccentric anomaly, and f0 is the true anomaly; in, is the position vector, is the angular momentum, H0 is the modulus of the angular momentum, r 0x r 0y r 0z are the components of the position vector on the x, y, and z axes, respectively, and h 0z is the component of angular momentum on the z-axis, and the function atan2() is the four-quadrant inverse tangent function.

9. The method according to claim 8, characterized in that The position vector and velocity vector of the satellite corresponding to epoch j in the fitted arc segment data Calculate the partial derivative matrix H of the satellite position with respect to the ephemeris orbit base point parameters k ,include: The partial differential derivative of the satellite position with respect to the orbital semi-major axis A0 at the reference time is: The semi-major axis The time difference between the onboard time and the reference time t k =tt oe ; Partial derivative of satellite position with respect to eccentricity e: Among them E k is the eccentric anomaly angle, r k is the corrected radial distance, n A is the corrected average angular velocity of the satellite; Orbital inclination of satellite position relative to reference time Partial derivative of : where Ω k is the longitude of the ascending node, i k is the corrected orbital inclination, u k is the corrected latitude argument; Partial derivative of satellite position with respect to parameter Ω0: Partial derivative of the satellite position with respect to the argument of perigee ω: The partial differential derivative of the satellite position with respect to the mean anomaly M0 at the reference time is: Satellite position change rate of semi-major axis Partial derivative of : where t k is the time difference between the onboard time and the reference time; The partial derivative of the satellite position with respect to the satellite average angular rate Δn0 at the reference time is: The rate of change of the satellite position relative to the difference between the satellite's average angular velocity and the calculated value Partial derivative of : Rate of change of satellite position with respect to the longitude of ascending node Partial derivative of : Where Ω0 is the longitude of the ascending node at the reference time; Satellite position versus orbital inclination rate of change Partial derivative of : The amplitude C of the sinusoidal harmonic correction of satellite position to latitude argument us Partial derivative of : The amplitude C of the cosine harmonic correction of the satellite position to the latitude argument uc Partial derivative of : The amplitude C of the sinusoidal harmonic correction of satellite position to orbit inclination is Partial derivative of : The amplitude C of the cosine harmonic correction of the satellite position to the orbital inclination ic Partial derivative of : The amplitude C of the sinusoidal harmonic correction of the satellite position to the parametric orbit radius rs Partial derivative of : The amplitude C of the cosine harmonic correction of the satellite position to the orbit radius rc Partial derivative of : Calculate the partial derivative matrix H of the satellite position with respect to the ephemeris orbit base point parameters k :

10. The method according to claim 9, characterized in that The partial derivative matrix H of the satellite position to the ephemeris orbit base point parameters k Iteratively calculate the ephemeris orbit base point parameters, including: Calculate the satellite position calculation value of each epoch under the current fitting initial value conditions Where j is the epoch identifier and k is the iteration number identifier; The iterative calculation equation is established as follows: H k ·Board k =L k In the above formula, H k Calculate the partial derivative matrix for k iterations, L k is the difference between the kth ephemeris calculation value and the fitting arc segment, Δσ k is the change of each parameter of the ephemeris orbit base point calculated in the kth iteration; Next, we use Givens transformation to calculate Δσ k Valuation, the formula is as follows: In the above formula, qr(H k ) represents the QR decomposition of the matrix in the brackets, and the subscripts in the above formula represent the elements at the corresponding positions in the matrix; The ephemeris orbit base point parameters after k iterations are calculated as follows: s k+1 =s k +Ds k ; Calculating the error vector involves: Error matrix L k The error vector at each epoch in is calculated according to the following formula: In the above formula, is the position vector of the satellite at the jth epoch in the fitting arc; The satellite position vector at the jth epoch calculated based on the ephemeris orbit base point parameters calculated based on the k-th iteration results; The conditions for the end of iterative calculation are as follows: When the error matrix L k The mean square error is less than the target value, or the absolute value of the change is less than the change threshold |Δσ k When |≤ε, or the number of iterations reaches the iteration limit, the update calculation is stopped.

Images