An Adaptive Variable-Parameter Calculation Method for Spacecraft Orbits

Abstract

The present invention discloses an adaptive variable parameter calculation method for spacecraft orbits, belonging to the field of aerospace technology. The upper bound of the number of collocation points within the total calculation interval is selected, and the upper bound of the number of collocation points is used as the number of collocation points for calculation, and the collocation point feedback iteration method is used for iterative calculation; the interval is divided according to the iterative calculation error situation and the number of collocation points in each sub-interval is determined; within the sub-interval where the orbit does not reach the target accuracy, the latest determined step size and the number of collocation points are used as the calculation parameters of the collocation point feedback iteration method for iterative calculation; according to the change of the iterative error in each sub-interval, it is determined whether each sub-interval needs to be further divided and whether to perform the iterative calculation of the next sub-interval, and the corresponding interval division operation or continuous calculation operation is executed, and the above process is repeated until the estimated orbits in all sub-intervals reach the target accuracy. The present invention uses the second derivative of the estimated orbit position state as the basis for determining the calculation parameters of the collocation point feedback iteration method, providing a high-efficiency calculation method for spacecraft orbit prediction that can adaptively select calculation parameters.

Description

Technical Field

[0001] The present invention relates to the technical field of aerospace, and particularly to an adaptive variable parameter calculation method for spacecraft orbits. Background Art

[0002] Nonlinear ordinary differential equations are one of the most common and important mathematical models, and many problems in the aerospace field are abstracted into this model. In actual aerospace engineering, the accuracy and efficiency of solving nonlinear ordinary differential equations have an important impact on the execution results of tasks. Among the numerical solution methods of nonlinear ordinary differential equations, a class of collocation iteration methods that combine the iteration method and the collocation method has received a great deal of attention due to its advantages of fast calculation speed and high calculation accuracy. The calculation efficiency of this method is usually severely affected by the step size and the number of collocation points. Good calculation parameters are the key to enabling such methods to always maintain a high calculation speed. However, the existing adaptive parameter selection strategies of collocation iteration algorithms often have various drawbacks and cannot meet the actual calculation needs. Summary of the Invention

[0003] Aiming at the problems existing in the prior art, the present invention provides an adaptive variable parameter calculation method for spacecraft orbits to improve the calculation efficiency of spacecraft orbits.

[0004] The present invention is realized through the following technical solutions:

[0005] An adaptive variable parameter calculation method for spacecraft orbits includes the following steps:

[0006] Step 1: Use the upper bound of the number of collocation points in the time interval of the orbit to be solved as the number of collocation points to be calculated, and determine the estimated orbit parameter values according to the collocation point time vector corresponding to the calculated number of collocation points;

[0007] Step 2: Construct a Chebyshev basis function square matrix according to the collocation point time vector corresponding to the calculated number of collocation points, and combine the collocation iteration method to iteratively correct the estimated orbit parameter values of the spacecraft at each collocation point time to obtain the orbit parameters at each collocation point time after correction;

[0008] Step 3: Calculate the two-norm of the difference vector of the orbit parameters of the spacecraft at each collocation point time after two adjacent iterative corrections. When the two-norm of the difference vector is less than the expected error, obtain the orbit parameters of the current time interval;

[0009] When the two-norm of the difference vector is greater than the expected error, and the count variable of the continuous increase times of the error is greater than the upper limit of the error increase times, the predicted orbit parameters of the current interval are obtained. According to the divergence state of the predicted orbit parameters, the number of collocation points in the current interval is adjusted. According to the adjusted number of collocation points, the corresponding collocation time vector, and the Chebyshev basis function square matrix, repeat step 3 to update the predicted orbit parameters until the two-norm of the difference vector is less than the set value, and the orbit parameters of the current time interval are obtained.

[0010] Step 4: Determine whether the current interval is the last time interval. If it is not the last time interval, repeat steps 2 and 3 to calculate the spacecraft orbit parameters of the next time interval until the orbit parameters of all time intervals are calculated, and the spacecraft orbit parameters of the entire time interval to be solved are obtained.

[0011] Preferably, in step 1, taking the spacecraft starting from the initial position and moving in a uniform straight line with the initial velocity as the initial condition, according to the initial condition and combined with the collocation time vector, the estimated values of the orbit states of the spacecraft at each collocation time are determined.

[0012] Preferably, the Chebyshev basis function square matrix in step 2 includes the Chebyshev basis function value square matrix, the Chebyshev basis function differential value square matrix, and the Chebyshev basis function integral value square matrix.

[0013] Preferably, the iteration formula is as follows:

[0014] x k+1 =x k -P(Dx k -g(x k ))

[0015] Where x is the orbit parameter, k is the number of iterations, the k + 1th is the current number of iterations, P is the Chebyshev basis function integral value square matrix, D is the Chebyshev basis function differential value square matrix, and g is the function value vector of the orbit dynamics system in the first-order form.

[0016] Preferably, the calculation method of the two-norm of the difference vector in step 3 is as follows:

[0017] e=||x k+1 -x k || 2

[0018] Preferably, in step 3, the time interval step size of the current iteration is less than or equal to the given minimum interval length, and the predicted orbit parameters continue to diverge from the first iteration. The number of collocation points in this time interval is increased to the number of collocation points M 1 , taking the spacecraft moving in a uniform straight line as the initial condition, and combined with the number of collocation points M 1And the corresponding collocation time vector and Chebyshev basis function square matrix are used to iteratively calculate the predicted orbit parameters again until the two-norm of the difference vector is less than the given error limit, and the orbit parameters of the current time interval are obtained.

[0019] Preferably, in step 3, the time interval step size of the current iteration is less than or equal to the given minimum interval length, and the predicted orbit parameters do not diverge continuously from the first iteration. The predicted orbit parameters are used to interpolate the spacecraft orbit parameter states at each newly added collocation time in the current time interval, and the interpolation result is used as the initial estimate of the spacecraft orbit parameters, combined with the number of collocations M 1 And the corresponding collocation time vector and Chebyshev basis function square matrix are used to iteratively calculate the orbit parameters of the spacecraft again until the two-norm e of the difference vector is less than the given error limit, and the orbit parameters of the current interval are obtained.

[0020] Preferably, in step 3, the time interval step size of the current time interval is greater than the given minimum interval length, and the predicted orbit parameters diverge continuously from the first iteration;

[0021] The current time interval is divided into two sub-intervals at the midpoint, and the number of collocations in the current interval is allocated to the two sub-intervals. Taking the spacecraft moving in a uniform straight line as the initial condition, starting from the first sub-interval after division, combined with the number of collocations and the corresponding collocation time vector and Chebyshev basis function square matrix of each sub-interval, the orbit parameters of the spacecraft are iteratively calculated again until the two-norm e of the difference vector is less than the given error limit, and the orbit parameters of the current interval are obtained.

[0022] Preferably, in step 3, the time interval step size of the current time interval is greater than the given minimum interval length, and the predicted orbit parameters do not diverge continuously from the first iteration;

[0023] According to the Chebyshev basis function square matrix, determine the number of maximum values of the second derivative vector of the orbit dynamics system state value. Divide the current time interval according to the number of maximum values of the second derivative vector, and ensure that the number of collocations in the divided intervals is not less than the minimum number of collocations. Use the predicted orbit parameters to interpolate the spacecraft orbit parameters at the collocation times in the divided intervals, and use the interpolation result as the initial value. Repeat step 2 until the two-norm e of the difference vector is less than the given error limit, and the orbit parameters of the current interval are obtained.

[0024] Preferably, the number of maximum value points of the second derivative vector of the current time interval is 0, and the current time interval is divided into two sub-intervals at the midpoint;

[0025] The number of maximum values of the second derivative vector of the current calculation interval is 1, and the current time interval is divided into two sub-intervals from the second derivative extreme point;

[0026] If the number of maxima of the second derivative vector in the current calculation interval is greater than 1, the current time interval is divided at the moments corresponding to the maxima and secondary maxima points of the second derivative to obtain three sub-intervals.

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

[0028] The present invention discloses an adaptive variable parameter calculation method for a spacecraft orbit. The upper bound of the number of collocation points within the total calculation interval is selected, and the upper bound of the number of collocation points is used as the number of collocation points for calculation, and the collocation point feedback iteration method is used for iterative calculation; according to the iterative calculation error situation in the previous step, the interval is divided and the number of collocation points in each sub-interval is determined; within the sub-interval where the target accuracy orbit has not been obtained, the latest determined step size and the number of collocation points are used as the calculation parameters of the collocation point feedback iteration method for iterative calculation; according to the change of the iterative error within each sub-interval, it is determined whether each sub-interval needs to be further divided and whether to perform the iterative calculation of the next sub-interval, and the corresponding interval division operation or the continued calculation operation is performed. Repeat the previous step and this step until the estimated orbits in all sub-intervals reach the target accuracy. The present invention uses the second derivative of the estimated orbit position state as the basis for determining the calculation parameters of the collocation point feedback iteration method, providing a high-efficiency calculation method that can adaptively select calculation parameters for spacecraft orbit prediction. BRIEF DESCRIPTION OF THE DRAWINGS

[0029] Figure 1 It is a schematic diagram of the principle of an adaptive variable parameter calculation method for a spacecraft orbit provided by the present invention;

[0030] Figure 2 It is a flow chart of an adaptive variable parameter calculation method for a spacecraft orbit provided by the present invention;

[0031] Figure 3 It is a schematic diagram of the calculation result of a spacecraft orbit provided by the present invention.

[0032] Figure 4 It is a schematic diagram of the position error of the calculation result of a spacecraft orbit provided by the present invention. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0033] The following further describes the present invention in detail with reference to the drawings, which is an explanation rather than a limitation of the present invention.

[0034] As Figure 1 and Figure 2 shown, the present invention provides an adaptive variable parameter calculation method for a spacecraft orbit, including the following steps:

[0035] Step 1, set the upper bound of the number of collocation points for the time interval of the orbit to be solved, set the number of collocation points for the calculation of the point feedback iteration method, and the number of collocation points for calculation is equal to the upper bound of the number of collocation points;

[0036] Step 2: Using the upper bound of the number of collocation points as the number of collocation points to be calculated and combining with the Chebyshev basis function matrix to calculate the orbital parameters, which specifically includes the following steps:

[0037] S2.1. Construct the collocation time vector according to the calculated number of collocation points The m-th element τ m in the collocation time vector τ has the following expression:

[0038]

[0039] S2.2. Construct the Chebyshev basis function matrix according to the collocation time vector. The Chebyshev basis function includes the Chebyshev basis function value matrix C, the Chebyshev basis function differential value matrix D, and the Chebyshev basis function integral value matrix P.

[0040] Construct the Chebyshev basis function value matrix C according to the collocation time vector. The construction rule is:

[0041] C = [cos(0arccos(τ T )), cos(1arccos(τ T )),..., cos((M - 1)arccos(τ T ))]

[0042] Construct the Chebyshev basis function differential value matrix D according to the collocation time vector and the Chebyshev basis function value matrix C. The construction rule is:

[0043]

[0044] Among them, each sub-matrix constituting the Chebyshev basis function differential value matrix D has the calculation formula When calculating the sub-matrix 1 , the Chebyshev basis function differential sub-matrix D 1 is constructed according to the rule D 1 1 = [D 1 2 , D 1 M ,...]. The formula for the differential Chebyshev basis function corresponding to the element in the first row of the k-th column of the Chebyshev basis function differential sub-matrix D 1 is (D 1 k ) 1,1 = (-1) k (k - 1) 2 . The formula for the differential Chebyshev basis function corresponding to the element in the M 1 -th row of the k-th column of the Chebyshev basis function differential sub-matrix D max is k = 1, L, M max ; The differential matrix D of Chebyshev basis functions 1 The k-th column D in 1 k Except for the first and last elements, the l-th element in the middle (D 1 k ) l,1 The calculation formula of the corresponding differential Chebyshev basis function is as follows:

[0045]

[0046] According to the collocation time vector τ and the Chebyshev basis function value matrix C, construct the Chebyshev basis function integral value matrix P, and the construction rule is:

[0047]

[0048] Among them, each sub-matrix that constitutes the Chebyshev basis function integral value matrix P The calculation formula is When calculating the sub-matrix 1 The construction rule of the Chebyshev basis function integral sub-matrix P Matrix P 1 The first column P in 1 1 The calculation formula of the corresponding integral Chebyshev basis function vector is P 1 1 = 1 + τ T P 1 The second column P in 1 2 The calculation formula of the corresponding integral Chebyshev basis function vector is P 1 For the r-th column P except the first two columns in 1 r The calculation formula of the corresponding integral Chebyshev basis function vector is

[0049] Calculate the first term 1 r in the formula of P The calculation formula is as follows:

[0050]

[0051] Calculate the second term 1 r in the formula of P The calculation formula is as follows:

[0052]

[0053] S2.3. Starting from the initial position of the spacecraft moving in a uniform straight line with an initial velocity as the initial condition, calculate the estimated values of the spacecraft's orbital states at each collocation time according to the initial condition and in combination with the collocation time vector.

[0054] Step 3. Set the iteration parameters of the collocation iteration method, and perform iterative correction on the estimated values of the spacecraft's orbital parameters at each collocation time according to the Chebyshev basis function matrix and in combination with the collocation iteration method to obtain the corrected orbital parameters at each collocation time. The iteration formula is as follows:

[0055] x k+1 = x k - P(Dx k - g(x k ))

[0056] where x k = [x 1 k (τ); x 2 k (τ); x 3 k (τ); x 4 k (τ); x 5 k (τ); x 6 k (τ); ], which is the estimated value of the six orbital parameters representing the position and velocity of the spacecraft at each collocation point after the k-th iterative correction. x k+1 is the estimated value of the spacecraft's orbital parameters at each collocation point after the (k + 1)-th iterative correction. P is the Chebyshev basis function integral matrix, D is the Chebyshev basis function differential matrix, g is the right-hand function value vector of the orbital dynamics system in the first-order form, and the elements of the g vector are the first-order derivative values of the corresponding elements in the x vector.

[0057] The iteration parameters include the iteration count variable k = 0, the error continuous increase count variable s = 0, the given error increase count upper limit s max , the minimum interval step size dt min , the minimum number of nodes M min , and the expected error limit e lim .

[0058] Step 4. Calculate the two-norm e of the difference vector of the spacecraft's orbital parameters at each collocation time after two adjacent iterative corrections. When the two-norm of the difference vector is less than the set value, obtain the orbital parameters of the current time interval.

[0059] When the two-norm of the difference vector is greater than the set value, and the count variable of the continuous increase times of the error is greater than the upper limit of the error increase times, the predicted orbit parameters of the current interval are obtained. According to the divergence state of the predicted orbit parameters, the number of collocation points in the current interval is increased. Based on the adjusted number of collocation points, the corresponding collocation point time vector, and the Chebyshev basis function square matrix, step 3 is repeated to update the predicted orbit parameters again until the two-norm of the difference vector is less than the set value, and the orbit parameters of the current time interval are obtained.

[0060] S4.1. Calculate the two-norm e of the orbit parameter difference vector of the spacecraft at each collocation point time after two adjacent iterative corrections;

[0061] The calculation formula of the two-norm of the difference vector is as follows:

[0062] e = ||x k+1 - x k || 2

[0063] where x is the orbit parameter, k is the number of iterations, and the (k + 1)th time is the current iteration number.

[0064] S4.2. Determine whether the value of the two-norm e of the orbit parameter difference vector obtained from this iterative calculation is greater than the two-norm e of the difference vector obtained from the previous iterative calculation:

[0065] If the two-norm e of the orbit parameter difference vector obtained from this iterative calculation is greater than the two-norm e of the orbit parameter difference vector obtained from the previous iterative calculation, let the value of the count variable s of the continuous increase times of the error increase by 1;

[0066] If the two-norm e of the orbit parameter difference vector obtained from this iterative calculation is less than the two-norm e of the orbit parameter difference vector obtained from the previous iterative calculation, let the value of the count variable s of the continuous increase times of the error be equal to 0.

[0067] S4.3. If the value of the count variable s + 1 of the continuous increase times of the error does not exceed the given upper limit s of the error increase times max , judge whether the two-norm e of the difference vector is less than the fixed error limit:

[0068] If the two-norm e of the difference vector is less than the fixed error limit, the orbit parameters of the current time interval are obtained.

[0069] If the two-norm e of the difference vector is greater than the given error limit, jump to step 3 to continue the iterative calculation until the two-norm e of the difference vector is less than the given error limit.

[0070] If the value of the count variable s of the continuous increase times of the error exceeds the given upper limit s of the error increase times max , and the two-norm of the difference vector is greater than the given error limit, record the predicted orbit parameters of the current interval, and execute step S4.4.

[0071] S4.4. If the step size of the current iterative calculation interval is less than or equal to the given minimum interval length dt min , determine whether the predicted orbital parameters diverge continuously starting from the first iteration according to the values of the iteration count variable K and the error continuous increase count variable S.

[0072] If the predicted orbital parameters diverge continuously starting from the first iteration, jump to step 3, increase the number of collocation points in this interval by 10% and record the number of collocation points in this interval as M 1 , take the spacecraft moving in a uniform straight line as the initial condition, and combine the number of collocation points M 1 and the corresponding collocation point time vector and Chebyshev basis function square matrix to re-iterate the calculation of the predicted orbital parameters until the two-norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

[0073] If the predicted orbital parameters do not diverge continuously starting from the first iteration, jump to step 3, interpolate the spacecraft orbital parameter state at each newly added collocation point time in the interval with the predicted orbital parameters, use the interpolation result as the initial estimate of the spacecraft state position and velocity, and combine the number of collocation points M 1 and the corresponding collocation point time vector and Chebyshev basis function square matrix to re-iterate the calculation of the spacecraft position until the two-norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

[0074] S4.5. If the current interval step size is greater than the given minimum interval length dt min , determine whether the calculation results diverge continuously starting from the first iteration:

[0075] If the predicted orbital parameters diverge continuously starting from the first iteration, divide the current interval into two sub-intervals at the midpoint, make the number of collocation points in the two sub-intervals equal to half of the number of nodes in the larger interval before division and round up. If the number of collocation points after division is less than M min , increase the number of collocation points to M min ; jump to step 3, take the spacecraft moving in a uniform straight line as the initial condition, start from the first sub-interval after division, and combine the number of collocation points in each sub-interval and the corresponding collocation point time vector and Chebyshev basis function square matrix to re-iterate the calculation of the spacecraft orbital parameters until the two-norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

[0076] If the predicted orbital parameters do not diverge continuously starting from the first iteration, multiply the function value vector g of the orbital dynamics system on the left by the Chebyshev basis function differential value square matrix to obtain the second derivative vector of the orbital parameters, and record the maximum value points of all the second derivative vectors of the orbital parameters;

[0077] Count the number of maximum points of the second derivative vector of the orbital parameters within the current time interval. If the number of maximum points of the second derivative vector in the current time interval is 0, then execute step S4.6;

[0078] S4.6: Divide the current interval into two sub-intervals at the midpoint, and set the number of collocation points in the two sub-intervals to be rounded up to half of the number of collocation points in the interval before division. If the number of collocation points in the two sub-intervals is less than the minimum number of nodes at this time, increase their number of collocation points to the minimum number of collocation points. Interpolate the spacecraft orbital parameters at each new collocation point time in the interval using the orbital parameters of the current interval, and use the interpolation result as the initial value. Jump to step 3 to recalculate until the two-norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

[0079] If the number of extreme points of the second derivative of the estimated solution of the spacecraft orbital parameters in the current calculation interval is 1, then execute step S4.7;

[0080] S4.7: Divide the current time interval into two sub-intervals at the extreme point of the second derivative, and set the number of collocation points in the two sub-intervals to be rounded up to half of the number of collocation points in the larger interval before division. If the number of collocation points in the two sub-intervals is less than the minimum number of collocation points at this time, increase the number of collocation points to the minimum number of collocation points. Interpolate the spacecraft orbital parameters at each new collocation point time in the interval using the orbital parameters of the current interval, and use the interpolation result as the initial value. Jump to step 3 to recalculate until the two-norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

[0081] If the number of extreme points of the second derivative of the estimated solution of the spacecraft orbital parameters in the current calculation interval is greater than 1, execute step S4.8;

[0082] S4.8: Divide the current interval into three small intervals at the corresponding times of the maximum point and the second maximum point of the second derivative, and proportionally divide and round up the number of collocation points in the current interval according to the length of each sub-interval. If the number of collocation points in a certain sub-interval is less than the minimum number of collocation points, increase the number of collocation points to the minimum number of collocation points. Interpolate the spacecraft orbital parameters at each collocation point time in the three new intervals using the orbital parameters of the current interval, and use the interpolation result as the initial value. Jump to step 3 to recalculate until the two-norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

[0083] Step 5: Determine whether the current interval is the last time interval. If it is not the last time interval, repeat steps 3 and 4 to calculate the spacecraft orbital parameters of the next time interval until the orbital parameters of all times are calculated to obtain the spacecraft orbital parameters of the entire time interval to be solved.

[0084] The embodiments of the present invention will be described below with reference to the accompanying drawings.

[0085] In view of the problem of calculating the orbit of a spacecraft, the present invention proposes a local collocation feedback iteration method, and its specific implementation includes the following steps:

[0086] The first step: For the orbit prediction problem of a low-earth satellite, select the maximum number of collocation points M according to engineering requirements max , iteration precision error limit, minimum interval length, minimum number of collocation points, maximum number of consecutive error increases.

[0087] The second step: Determine the collocation times according to the selected maximum number of collocation points, and construct a collocation time vector τ, a differential square matrix D, and an integral square matrix P. Set the orbit parameters to be iteratively corrected as a 6M max -dimensional column vector, where the first to the M max th elements of the vector correspond to the state estimation values of the spacecraft's x-axis coordinate at each collocation point, and the M max +1 to the 2M max th elements of the vector correspond to the state estimation values of the spacecraft's y-axis coordinate at each collocation point, and the 2M max +1 to the 3M max th elements of the vector correspond to the state estimation values of the spacecraft's z-axis coordinate at each collocation point, and the 3M max +1 to the 4M max th elements of the vector correspond to the state estimation values of the spacecraft's x-axis velocity component at each collocation point, and the 4M max +1 to the 5M max th elements of the vector correspond to the state estimation values of the spacecraft's y-axis velocity component at each collocation point, and the 5M max +1 to the 6M max th elements of the vector correspond to the state estimation values of the spacecraft's z-axis velocity component at each collocation point;

[0088] The third step: Substitute the iteratively corrected state variables into the equation x k+1 = x k - P(Dx k - g(x k )) for iterative calculation.

[0089] Step 4: After each iteration, calculate the 2-norm e of the difference vector between the orbital parameter estimation values of the previous and current iterations, record the number of iterative calculations, and determine whether the 2-norm e of the difference vector is less than or equal to the predetermined accuracy. If the 2-norm e of the difference vector is less than or equal to the predetermined accuracy, record the iterative result. If the 2-norm e of the difference vector continuously increases and the number of consecutive increases exceeds the maximum number of consecutive increases in error, terminate the iterative calculation, record the current calculation result for use as the initial value in subsequent calculations for warm start, and calculate the second derivative value of the current result and find the maximum point of the second derivative of the state vector except for the boundary points. If the step size of the above calculation interval is less than a given minimum interval length and the calculation results have been continuously divergent since the first iteration, increase the number of nodes in the interval, assume that the spacecraft moves in a uniform straight line, determine the initial estimated values of the state of each collocation point accordingly, and recalculate; if the step size in the above calculation is less than a given minimum interval length, but the calculation results have not been continuously divergent since the first iteration, increase the number of nodes in the interval, interpolate the new collocation points with the current interval calculation results, and use the interpolation results as the initial values and recalculate; for the case where the step size of the interval is greater than a given minimum interval length and the calculation results have been continuously divergent since the first iteration, divide the current interval into two sub-intervals at the midpoint, and set the number of nodes in the two sub-intervals to the collocation node numbers respectively. Starting from the first sub-interval after the division, assume that the spacecraft moves in a uniform straight line, determine the iterative initial values of the orbital parameters at each collocation time accordingly, and re-iterate and update the orbital parameters at each collocation time in this sub-interval; for the case where the step size of the interval is greater than a given minimum interval length, but the calculation results do not start to diverge from the first iteration, and there is no extreme point of the second derivative in the interval, divide the current interval into two sub-intervals at the midpoint, and set the number of nodes in the two sub-intervals to be rounded up to half of the number of nodes in the larger interval before the division. If the number of nodes in the two sub-intervals is less than the minimum number of nodes at this time, increase their number of nodes to the minimum number of nodes, interpolate each collocation point in the two new intervals with the previous calculation results, and use the interpolation results as the initial values and recalculate; for the case where the step size of the interval is greater than a given minimum interval length, but the calculation results do not start to diverge from the first iteration, and there is one extreme point of the second derivative in the interval, divide the current interval into two sub-intervals at the extreme point of the second derivative, and set the number of nodes in the two sub-intervals to be rounded up to half of the number of nodes in the larger interval before the division. If the number of nodes in a certain sub-interval is less than the minimum number of nodes at this time, increase the number of nodes in it to the minimum number of nodes, interpolate the orbital parameters corresponding to each collocation point in the two new intervals with the current calculation results, and use the interpolation results as the initial values and recalculate.For the case where the interval step size is greater than a given minimum interval length, but the calculation result does not diverge at the beginning of the first iteration and there is more than one extreme point of the second derivative within the interval, the current interval is divided into three sub - intervals through the maximum point and the second - maximum point of the second derivative. And let the nodes of the larger interval before division be proportionally divided according to the lengths of the sub - intervals and rounded up within the three sub - intervals. If the number of nodes in a certain sub - interval is less than the minimum number of nodes at this time, the number of nodes is increased to the minimum number of nodes. Interpolate the orbital parameters at each collocation time in the three new intervals with the current calculation results, and use the interpolation results as the initial values for calculation.

[0090] Step 5: Determine whether the number of sub - intervals for which the calculation results of the orbital parameters have reached the target calculation accuracy is equal to the total number of sub - intervals. If the number of sub - intervals for which the orbital calculation results have reached the target calculation accuracy is equal to the total number of sub - intervals, output the calculation results and terminate the calculation; if the number of sub - intervals for which the orbital calculation results have not reached the target calculation accuracy is equal to the total number of sub - intervals, use the latest determined step size and the number of collocation points as the calculation parameters for the collocation feedback iteration method in the first sub - interval where the target accuracy orbit has not been obtained for iterative calculation until the calculation results of the orbital parameters in all sub - intervals reach the target calculation accuracy.

[0091] Example 1: The dynamic equation for the problem of the orbit recurrence of a near - earth satellite under the 40 - order EGM2008 Earth gravity field model is:

[0092]

[0093] where \(r = [x,y,z]\) T is the position vector of the near - earth satellite, \(\mu=3986004.418\times10\) 8 \(m\) 3 / s 2 , representing the Earth's gravitational constant, \(r = ||r||\), representing the distance between the Earth's center and the center of mass of the near - earth satellite, \(t\) 0 is the initial time of orbit recurrence, \(a\) is the perturbation term. In this example, the perturbation of the 40 - order EGM2008 Earth gravity field model is considered.

[0094] This example adopts the initial calculation conditions as shown in Table 1, the calculation parameters as shown in Table 2, and the 24 - hour operating orbit of the near - earth satellite obtained by calculation is as Figure 3 and Figure 4 shown.

[0095] Table 1 Initial calculation conditions adopted in Example 1

[0096] Parameter Value <![CDATA[Initial position coordinate r 0 (m)]]> <![CDATA[[-0.3889×10 6 ,7.7388×10 6 ,0.6736×10 6 T > ​ <![CDATA[Initial velocity vector v 0 (m / s)]]> <![CDATA[[-3.5794×10 3 ,0,6.7997×10 3 T > ​

[0097] Table 2 Calculation parameters adopted in Example 1

[0098] Parameter Value Satellite operation time / s 7000 Upper bound of the number of collocation points / piece 200 Iterative precision error limit / m <![CDATA[10 -6 > Maximum number of consecutive increases in the maximum error 10 Minimum interval length / s 3 Minimum number of nodes / piece 4

[0099] The present invention discloses an adaptive variable parameter calculation method for spacecraft orbits. The calculation method includes the following steps: Select the upper bound of the number of collocation points within the total calculation interval, and use the upper bound of the number of collocation points as the number of collocation points for calculation, and perform iterative calculation using the collocation point feedback iteration method; Divide the interval according to the iterative calculation error situation and determine the number of collocation points in each sub-interval; In the sub-intervals where the orbit does not reach the target accuracy, use the latest determined step size and the number of collocation points as the calculation parameters of the collocation point feedback iteration method for iterative calculation; According to the change of the iterative error in each sub-interval, determine whether each sub-interval needs to be further divided and whether to perform iterative calculation in the next sub-interval, and perform the corresponding interval division operation or continue the calculation operation. Repeat the above process until the estimated orbits in all sub-intervals reach the target accuracy. The present invention uses the second derivative of the estimated orbit position state as the basis for determining the calculation parameters of the collocation point feedback iteration method, providing a high-efficiency calculation method for spacecraft orbit prediction that can adaptively select calculation parameters.

[0100] The above content is only to illustrate the technical idea of the present invention, and the protection scope of the present invention cannot be limited thereby. Any modification made on the basis of the technical solution according to the technical idea proposed by the present invention falls within the protection scope of the claims of the present invention.

Claims

1. An adaptive variable parameter calculation method for spacecraft orbits, characterized in that, it includes the following steps: Step 1: Use the upper bound of the number of collocation points in the time interval of the orbit to be solved as the number of collocation points to be calculated, and determine the orbit parameter estimation value according to the collocation point time vector corresponding to the calculated number of collocation points; Step 2: Construct a Chebyshev basis function square matrix according to the collocation point time vector corresponding to the calculated number of collocation points, and combine the collocation point iteration method to iteratively correct the orbit parameter estimation values of the spacecraft at each collocation point time to obtain the orbit parameters at each collocation point time after correction; Step 3: Calculate the two-norm of the difference vector of the orbit parameters of the spacecraft at each collocation point time after two adjacent iterative corrections. When the two-norm of the difference vector is less than the expected error, obtain the orbit parameters of the current time interval; When the two-norm of the difference vector is greater than the expected error, and the error continuous increase count variable is greater than the error increase count upper limit, obtain the estimated orbit parameters of the current interval. According to the divergence state of the estimated orbit parameters, adjust the number of collocation points of the current interval. According to the adjusted number of collocation points, the corresponding collocation point time vector, and the Chebyshev basis function square matrix, repeat Step 3 to update the estimated orbit parameters until the two-norm of the difference vector is less than the set value, and obtain the orbit parameters of the current time interval; Step 4: Determine whether the current interval is the last time interval. If it is not the last time interval, repeat Steps 2 and 3 to calculate the orbit parameters of the spacecraft in the next time interval until the orbit parameters of all time intervals are calculated to obtain the orbit parameters of the spacecraft in the entire time interval to be solved.

2. An adaptive variable parameter calculation method for spacecraft orbits according to claim 1, characterized in that, in Step 1, taking the spacecraft starting from the initial position and moving in a uniform straight line at the initial velocity as the initial condition, and determining the orbit state estimation values of the spacecraft at each collocation point time according to the initial condition and in combination with the collocation point time vector.

3. An adaptive variable parameter calculation method for spacecraft orbits according to claim 1, characterized in that, the Chebyshev basis function square matrix in Step 2 includes a Chebyshev basis function value square matrix, a Chebyshev basis function differential value square matrix, and a Chebyshev basis function integral value square matrix.

4. An adaptive variable parameter calculation method for spacecraft orbits according to claim 1, characterized in that, the iterative formula is as follows: x k+1 = x k -P(Dx k -g(x k )) where x is the orbit parameter, k is the number of iterations, the (k + 1)th time is the current number of iterations, P is the Chebyshev basis function integral value square matrix, D is the Chebyshev basis function differential value square matrix, and g is the function value vector of the orbit dynamics system in the first-order form.

5. An adaptive variable parameter calculation method for spacecraft orbits according to claim 1, characterized in that, the calculation method of the two-norm of the difference vector in Step 3 is as follows: e = ||x k+1 -x k || 2 。 6. An adaptive variable parameter calculation method for spacecraft orbits according to claim 1, characterized in that, In step 3, the time interval step size of the current iteration is less than or equal to the given minimum interval length, and the estimated orbital parameters continue to diverge since the first iteration. Increase the number of collocation points within this time interval to the number of collocation points M 1 , with the spacecraft moving in a uniform straight line as the initial condition, and combined with the number of collocation points M 1 and the corresponding collocation point time vector and Chebyshev basis function square matrix, re-iterate the calculation of the estimated orbital parameters until the two-norm of the difference vector is less than the given error limit, and obtain the orbital parameters of the current time interval.

7. An adaptive variable parameter calculation method for spacecraft orbits according to claim 6, characterized in that, In step 3, the time interval step size of the current iteration is less than or equal to the given minimum interval length, and the predicted orbital parameters do not continuously diverge starting from the first iteration. Interpolate the spacecraft orbital parameter states at each newly added collocation point moment in the current time interval using the predicted orbital parameters, and use the interpolation result as the initial estimate of the spacecraft orbital parameters. Combine the number of collocation points M 1 and the corresponding collocation point moment vector and Chebyshev basis function square matrix to re-iterate the calculation of the spacecraft's orbital parameters until the two-norm e of the difference vector is less than the given error limit, and obtain the orbital parameters of the current interval.

8. An adaptive variable parameter calculation method for spacecraft orbits according to claim 1, characterized in that, In step 3, the step size of the current time interval is greater than the given minimum interval length, and the estimated orbital parameters diverge continuously since the first iteration; Divide the current time interval into two sub - intervals at the mid - point, and distribute the number of collocation points of the current interval to the two sub - intervals. Taking the spacecraft moving in a uniform straight - line motion as the initial condition, starting from the first sub - interval after division, recompute the orbital parameters of the spacecraft iteratively by combining the number of collocation points of each sub - interval, the corresponding collocation time vector, and the Chebyshev basis function square matrix until the two - norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

9. An adaptive variable - parameter calculation method for a spacecraft orbit according to claim 8, characterized in that, In step 3, the step size of the current time interval is greater than the given minimum interval length, and the estimated orbital parameters do not diverge continuously since the first iteration; Determine the number of maximum values of the second - derivative vector of the orbital dynamics system state value according to the Chebyshev basis function square matrix, divide the current time interval according to the number of maximum values of the second - derivative vector, and ensure that the number of collocation points in the divided intervals is not less than the minimum number of collocation points. Interpolate the orbital parameters of the spacecraft at the collocation times in the divided intervals using the estimated orbital parameters, and use the interpolation result as the initial value. Repeat step 2 until the two - norm e of the difference vector is less than the given error limit to obtain the orbital parameters of the current interval.

10. An adaptive variable - parameter calculation method for a spacecraft orbit according to claim 9, characterized in that, The number of maximum - value points of the second - derivative vector of the current time interval is 0, and divide the current time interval into two sub - intervals at the mid - point; The number of maximum values of the second - derivative vector of the current calculation interval is 1, and divide the current time interval into two sub - intervals at the second - derivative extreme - value point; The number of maximum values of the second - derivative vector of the current calculation interval is greater than 1, and divide the current time interval at the times corresponding to the second - derivative maximum - value points and the secondary - maximum - value points to obtain three small intervals.

Images