Abstract
In the process of Earth-Moon transfer orbit, many parameters are involved and the degree of influence varies among the parameters. How to accurately distinguish the influence relationship between different parameters is of great significance to engineering missions. Based on Sobol sequence sampling method and Sobol global sensitivity analysis method, a calculation process of global sensitivity analysis is proposed in this paper for high-fidelity Earth-Moon transfer orbit. A numerical simulation method is used to verify that the Sobol sequence sampling method has better convergence and higher precision than other sampling methods and has better adaptability in global sensitivity analysis. The effects of different state parameters and the combination of different parameters on the perilune parameters of Earth-Moon transfer orbit are given by simulation examples, which verifies the effectiveness and feasibility of the calculation process proposed in this paper. Simulation results show that the radial position and tangential velocity at the trans-lunar injection point are the main sensitivity parameters, and the other parameters have little effect on the results. The sensitivity of the orbital elements at the trans-lunar injection point to the perilune parameters is different and needs to be determined according to the specific parameters. The results of this study can provide important reference for future Earth-Moon transfer orbit design and related engineering missions.
1. Introduction
As the nearest celestial body to the Earth, the Moon is the first choice for human exploration and a suitable transfer station for deep space exploration [1, 2]. In recent years, along with the advance of space detection technique, application value, and development potential of the Moon appears gradually, many countries have put forward their own lunar exploration plans [3, 4], including “Artemis” program in the United States [5] and China’s “International Lunar Research Station” project plan [6], which starts a new round of worldwide lunar exploration boom.
In lunar exploration, especially manned lunar exploration related missions, Earth-Moon transfer orbit is a channel enabling the spacecraft to fly from the Earth to the Moon, and the rationality of its design is related to the safety of the whole mission directly [7, 8]. In practical missions, the transfer orbit in the Earth-Moon space is affected by a combination of the Earth’s gravitation, the Moon’s gravitation, the Sun’s gravitation, and other disturbing forces and also tends to involve a lot of engineering constraints and orbital constraints. In order to ensure the accuracy and rationality of the orbit design results, a high-fidelity dynamics model is often used for calculation. However, because the high-fidelity dynamics model itself has a strong nonlinearity, some minor changes in parameters may lead to great differences in results. The uncertainty of influence degree among parameters not only makes it difficult to further understand the model but also makes it not beneficial for engineers to adjust and optimize the important influence parameters in the course of orbit design. Therefore, the parametric sensitivity analysis of the Earth-Moon transfer orbit under the high-fidelity model is carried out to discuss the influence degree of the initial state parameters of the Earth-Moon transfer orbit on the output state parameters, which is helpful for engineers to identify the important influence parameters, and it is of great significance for engineering practice and optimization of Earth-Moon transfer orbit.
A parametric sensitivity analysis method is a method to study mathematical models based on mathematical statistics, which is used to analyze the impact of each input uncertainty of the system on the output uncertainty of the model [9–11]. It is widely used in hydrological process simulation [12], material structure analysis [13], risk assessment [14], and other fields. Parametric sensitivity analysis methods are mainly divided into local sensitivity analysis method and global sensitivity analysis method [15]. The local sensitivity analysis method is simple and has a small amount of calculation, but it can only analyze the influence of a single parameter on the output results, ignoring the influence of the combination between parameters on the results, and cannot measure the influence of input variables on the output variables within the whole feasible range, so the application range is limited [16, 17]. The global sensitivity analysis method is used to analyze the influence degree of input parameters on output parameters within the feasible range and allows the parameters to change at the same time. Therefore, it can not only analyze the influence of individual parameters but also analyze the influence of combination between parameters on output results [18, 19]. Compared with local sensitivity analysis, global sensitivity analysis can reflect the characteristics of the system more accurately. Common global sensitivity analysis methods include Sobol method [20, 21], Morris method [22, 23], Fourier Amplitude Sensitivity Test (FAST) method [24], and Regression Analysis method [25]. The Sobol method is mainly based on variance decomposition to calculate the sensitivity of each parameter, whose form is relatively simple, and accurate sensitivity analysis results can be given by counting the input and output results of the system. Therefore, it has become a widely used global sensitivity analysis method nowadays [26–28].
In the field of aerospace, Dell’Elce and Kerschen studied the problem of quantifying probabilistic uncertainty in the estimation of the orbital lifetime of LEO satellites using parametric sensitivity analysis and obtained the effect of uncertainty on the orbital lifetime [29]. Somma et al. studied the effect of the low Earth orbit environment on the changes in the number and distribution of newly launched satellites by means of parametric sensitivity analysis [30]. Yuan et al. used a covariance-based method to analyze the sensitivity of dynamic parameters on landing accuracy for landing dynamics model of a small celestial body and finally obtained qualitative results [31]. Taking the Inertial Measurement Unit (IMU) and radio ranging integrated navigation scheme as the background, Lou et al. established the sensitivity matrix and perturbation matrix of the detector state to the percentage parameter and analyzed the sensitivity of the state estimation error of the detector on the uncertainty of the lift-drag ratio and the law of its distribution [32]. Bai analyzed the characteristics of the Earth-Moon transfer orbit under the two-body model, the relationship curves between each design parameter and the input parameters are given by changing the values of individual design parameters, and the degree of influence of each design parameter on the output results is evaluated qualitatively [33]. In the study of Mars entry problem, Huang used the Sobol index calculation method to analyze the sensitivity of opening point parameters to the uncertainty of key input parameters comprehensively, but the influence of the combination between parameters on the results is not further analyzed [34]. Zhao et al. used the state transition matrix between constraint variables and design variables to analyze the sensitivity of various parameters in the study of multiflight trajectory of launch vehicle, but the state transition matrix only reflected the analysis results at the selected local state points [35].
To sum up, some scholars have carried out studies on parametric sensitivity analysis in the aerospace field, but there are still some shortcomings: First, most studies still use methods such as state transition matrix for analysis. The analysis process can only reflect the influence of single parameter on the results, so the influence of combination between parameters on the output results cannot be further obtained. Secondly, some studies only provide qualitative conclusions for sensitivity analysis results, without further quantitative analysis results, which has limited guiding significance for engineering. Thirdly, related studies on the design of Earth-Moon transfer orbit have only analyzed the orbit characteristics, without systematic and in-depth analysis of the sensitivity among parameters, and there are few sensitivity analysis studies on high-fidelity Earth-Moon transfer orbit. In view of the above shortcomings, this paper carries out global sensitivity analysis of high-fidelity Earth-Moon transfer orbit based on the Sobol method, and all of the above three shortcomings are improved by the method proposed. The main contributions are as follows: First, a calculation process of global sensitivity analysis for high-fidelity Earth-Moon transfer orbit is proposed by using the Sobol sequence sampling method and Sobol global sensitivity analysis method, and the solution process is intuitive. Second, the sensitivity of each parameter is ranked qualitatively by the global sensitivity analysis, and the results of the quantitative sensitivity analysis of each parameter are also given, which reflects how much each parameter can influence the results in the given interval range for given parameters. Third, the influence of different sampling methods on the calculation results of global sensitivity analysis is compared, and the influence of different trans-lunar injection point parameters and the combination of each parameter on the perilune parameters is given.
This paper is divided into four sections. The first section introduces the scheme of Earth-Moon transfer orbit and the high-fidelity dynamic model used in calculating the transfer orbit. The second section gives the mathematical principles of Sobol sequence sampling method and Sobol global sensitivity analysis method. Based on the two methods, a global sensitivity analysis calculation process for high-fidelity Earth-Moon transfer orbit is given. In the third section, according to the proposed calculation process, the global parametric sensitivity simulation analysis is carried out for a specific Earth-Moon transfer orbit case. The last section summarizes the whole study and gives the corresponding conclusions.
2. Problem Description
2.1. Description of Earth-Moon Transfer Orbit Scheme
The Earth-Moon transfer orbit is the orbital segment in which a spacecraft starts from the Earth’s parking orbit and then escapes the Earth’s gravitation to the Moon, until it is finally captured by the Moon’s gravitation and brakes near the Moon. At present, there are mainly four orbital transfer modes proposed or adopted: Hohmann Transfer Orbit, Transfer Orbit based on Collinear Translational Point L1 Corridor, Weak Stability Boundary (WSB) Transfer Orbit, and Low Thrust Transfer Orbit [8]. The Hohmann transfer method is the main transfer method for Earth-Moon transfer orbit at present because of its simple process and low requirements for engineering implementation. The specific implementation of this scheme is as follows: Spacecraft is firstly placed in a low Earth orbit, and then an impulse is applied at the appropriate trans-lunar injection (hereinafter referred to as TLI) point A to accelerate the spacecraft into an Earth-Moon transfer orbit. After reaching the perilune B, a coplanar impulse is applied to enter the low lunar orbit. The whole flight process is shown in Figure 1.
It can be seen from Figure 1 that in practical engineering problems, the TLI point, as a main control point in the project, has great influence on the whole transfer process of Earth-Moon transfer orbit and terminal arrival parameters. Therefore, the study of the influence of various state parameters at the TLI point on the Earth-Moon transfer orbit results is helpful for us to have a deeper understanding of the characteristics of the Earth-Moon transfer orbit. Based on the above considerations, the state parameters of spacecraft at the TLI point are chosen as the research object of sensitivity analysis in this paper.
2.2. Dynamic Model
For the Earth-Moon transfer orbit problem studied in this paper, the high-fidelity dynamic model is mainly chosen as the dynamic model, which is also the recognized accurate orbit calculation model. The high-fidelity dynamic model of spacecraft in geocentric J2000 coordinate system is expressed as where is the Earth’s gravitational constant, is the geocentric position vector, is gravitational perturbation acceleration of body, is the nonspherical perturbation acceleration of the Earth, is the nonspherical perturbation acceleration of the Moon, is the perturbation acceleration of solar radiation pressure, and is the perturbation acceleration of atmospheric drag. The calculation method of various perturbation terms is described in detail in literature [36], and the specific formula used in this paper is shown in Appendix A. In this paper, the RKF78 method is used for integral calculation. The RKF78 method gives two sets of RK formulas of 7 and 8 orders at the same time. The local truncation error is estimated by the difference of the integral calculated by the two sets of formulas, so the step size of the next step can be given and the step size can be automatically controlled.
3. Global Sensitivity Analysis Method
The Sobol global sensitivity analysis method is a method to study the mathematical model based on mathematical statistics; within a given feasible range, the sensitivity analysis of parameters can be performed directly by sampling the sample points of input and output parameters of the model. Therefore, the Sobol global sensitivity analysis mainly consists of two stages. The first stage is to select the sample points of initial input parameters. The second stage is to substitute the selected sample points into the model to obtain the corresponding output results and then conduct sensitivity analysis according to the calculation results. This section firstly introduces the Sobol sequence sampling method and then proposes a calculation process of high-fidelity global sensitivity analysis for Earth-Moon transfer orbit based on the Sobol global sensitivity analysis method.
3.1. Sobol Sequence Sampling Method
A global sensitivity analysis method needs to use a large number of sample points to ensure the convergence and the accuracy of the calculation results, and the random selection of sampling method may lead to nonconvergence of computational results or slow convergence process and waste a lot of computing resources. Therefore, it is necessary to select an appropriate sampling method before conducting global sensitivity analysis. Sobol sequence has become a common sampling method for global sensitivity analysis due to its homogeneity of distribution in multidimensional sampling space and better simulation accuracy [37, 38].
Sobol sequence is a quasirandom low-discrepancy sequence proposed by mathematician I. M. Sobol in 1967 [39]. Different from ordinary random sequences, quasirandom low-discrepancy sequences are generated based on specific algorithms. It has good uniformity of distribution in sampling space, which can be extended to any dimension, and has better computational efficiency and accuracy than the traditional Monte Carlo method in large-scale sampling simulation. A quasirandom low-discrepancy sequence represents a series of integers as bits based on a certain number; then the integer has a unique -ary representation which is expressed as where and represents the rounding operation, . In Sobol sequence, integers are expanded in base 2, take 0 or 1, that is,
The primitive polynomial of degree is defined as where is the primitive polynomial of the -th dimension, is the degree of freedom of the polynomial, and are 0 or 1. A set of integers can be obtained by recursion of the coefficient of the primitive polynomial; the recursion is as follows: where represents a binary bitwise XOR operation and the initial values can be chosen which freely provided that is an odd number less than . Direction number is defined as
The Sobol sequence based on the direction number can be obtained by where is the -th element of the generated Sobol sequence and is the -th digit from the right when is written in binary ; here, we use the notation to denote the binary representation of numbers [40]. In this paper, we use six and twelve dimensional Sobol sequences for sensitivity calculations, respectively, and the polynomials used are different for each dimensional parameter. The primitive polynomials used in this paper are given in Appendix B.
3.2. Sobol Global Sensitivity Analysis Method
The Sobol global sensitivity analysis method was first proposed by a Russian scholar Sobol in 1993 [21]. The essence of this method is to use the variance decomposition method to decompose the implicit function of the system into functions of single parameter and the combination of parameters, and then the total variance of the system can be expressed as the variance of each single parameter and the combination of parameters. The importance of each parameter and the combination between parameters on the output of the system model can be analyzed by calculating these variances. The basic principles of Sobol global sensitivity analysis method for high-fidelity Earth-Moon transfer orbit model are as follows:
For the complex high-fidelity Earth-Moon transfer orbit model, the system model is treated as a black box, and the model is expressed as a function ; is the state parameter of the TLI point of the Earth-Moon transfer orbit represented by vector or where represents the position parameters of the spacecraft at the TLI point, represents the velocity parameters of the spacecraft at the TLI point, and represents the geocentric distance, eccentricity, inclination, right ascension of ascending node, argument of perigee, and true anomaly of the TLI point. The state parameter vectors can be expressed as , is the number of parameters, and represents perilune parameter of Earth-Moon transfer orbit. In general, any input parameter space can be converted to a unit cube . For the Earth-Moon transfer orbit model, the state parameters at TLI point can be regarded as independent of each other. Therefore, the function can be decomposed according to its parameters and the combination of parameters and express it as the sum of its subterms where is constant, is the function of , is the function of and , and the rest are analogized. One condition that satisfies this decomposition is
The formula states that the integral of each subterm over any contained parameter is zero. According to the above two formulas, each subterm is orthogonal to each other; if the following equation can be obtained
Sobol has proved in his paper [11] that the decomposition of Equation (10) is unique, and each term can be evaluated by the following multiple integrals: where represents parameters other than and . Similarly, other higher-order terms can be solved.
Further assuming that is square integrable, the total variance is calculated as
The deviation of can be obtained from each subterm of the above equation
Then, the total variance can be decomposed into where
The expressions for the remaining terms are similar. represents the variables other than ; is a variable of all parameters except and . The sensitivity index of each order can be obtained by dividing the items in the above formula by the total variance
The first-order sensitivity index , the second-order sensitivity index , and the global sensitivity index are defined as where represents variance and represents expectation. The first-order sensitivity index represents the influence of a single model parameter on the model output, represents the influence of the combination of two model parameters on the model output, and the global sensitivity index represents the influence of all the combination of parameters including the -th model parameter on the model output. For example, consider a function with three parameters. The first-order sensitivity index represents the degree of influence on the results when the parameters are acting individually, which reflects how much each parameter can influence the results in the given range of the three parameters. represents the degree of influence on the results under the interaction between the parameters and . Similarly, represents the degree of influence on the results under the interaction between the parameter and , respectively. The global sensitivity index of is
This result reflects the degree of influence of the combined effect of parameters and other parameters on the results. are the global sensitivity index of , and the meaning is the same as above. According to the above results, the impact of each parameter and the combination of parameters on the model output results can be analyzed quantitatively.
In addition, it can be seen from the above expression that it needs to use multiple integrals to solve for each partial variance, which will make the actual calculations very complicated. Therefore, Saltelli et al. proposed a simplified Monte Carlo method [41], which required only single-layer sampling of sample points, greatly reducing the amount of calculation. The calculation steps of this method are as follows: (1)Firstly, the sampling is performed times within a given feasible range to generate a sample matrix in total. Then, two matrices and of size are defined separately, each containing half of the samples(2)Construct a new matrix from matrices and , where the -th column is from matrix and the remaining columns are from matrix :(3)Substitute the above sample matrices into the system model, and output vectors can be obtained:(4)Substitute the above results into the simplified Monte Carlo calculation formula to calculate the variances of each order:
Equations (25), (26), and (27) are the simplified formulas needed to calculate the first-order, global, and second-order sensitivity indices, respectively. The above calculation method avoids the multiple integration of the partial variance, simplifies the calculation process of Sobol global sensitivity analysis method, and also has high accuracy.
From the principle of the above sensitivity analysis method, it can be seen that the sensitivity analysis first needs to sample the sample points within a given range. Due to the many parameters of the Earth-Moon transfer orbit and the complex constraints, it is not only impossible to meet the engineering requirements by arbitrarily choosing the range of parameters but also very easy to cause the calculation results to be so different from the real situation that the sensitivity analysis cannot be performed. Therefore, based on the above-mentioned Sobol global sensitivity analysis method, a global sensitivity analysis process under the high-fidelity Earth-Moon transfer orbit model is established. Firstly, the required Earth-Moon transfer orbit can be obtained by preliminary calculation, and the spacecraft state parameters at the TLI point are selected as the input parameters. Then, the deviation of the initial parameters at the TLI point is calculated under the premise of the engineering constraints. Based on this, the range of values of each input parameter is determined and sample points are sampled to ensure the reasonableness of the calculation results. Then, the sample points are substituted into the high-fidelity orbit model to calculate the output parameters of the spacecraft at the perilune, and finally, the parametric sensitivity analysis is performed based on the calculated input and output of the sample points. The basic flow of the whole analysis process is shown in Figure 2.
4. Simulation Analysis
In this section, the applicability of the Sobol sampling method in the global sensitivity analysis is first verified through simulation analysis, and then, the global sensitivity characteristics under different state parameters at the TLI point are analyzed for the Earth-Moon transfer orbit.
4.1. Sampling Method Validation
Different sampling methods have different effects on the sensitivity analysis results. In order to compare the effects of each sampling method on the sensitivity analysis results more clearly, four sampling methods are compared: simple random sampling, Latin hypercube sampling (LHS) [42], Halton sequence [43], and Sobol sequence. The Latin hypercube sampling method is a method based on the idea of stratified sampling, which divides the sampling space into equal parts according to the number of samples to obtain sublayers. Sampling points are randomly selected in each layer, and finally, the obtained sampling points are disrupted and reordered. Halton sequences are similar to Sobol sequences in that they are essentially quasirandom low-difference sequences, with the difference that the bases of Halton sequences are a series of prime numbers.
The distribution of each sampling method in the sampling area is compared firstly. The sampling method is chosen as two-dimensional sampling and the value spaces is [0,1]. The number of sampling points are selected as 100, 1000, and 4000, where the base number of Halton sequence is chosen as 2 and 3, respectively. Figure 3 shows the distribution of sampling points in the sampling area generated by different sampling methods.
(a)
(b)
(c)
Then, the effect of different sampling methods on the results of parametric sensitivity analysis is further analyzed. Selecting the simple function , the values of and are both in the range [0,1], and the standard result of global sensitivity index for both parameters is easy to know which is 0.5. The Sobol global sensitivity analysis method is used to analyze the parametric sensitivity of parameters and on the output results, and the error of the calculated results relative to the reference value is counted at different numbers of sampling points.
Figure 4 shows the results obtained from sensitivity calculations using the sampling points generated by the Latin hypercube sampling method and the simple random sampling method are poorly converged, and an error value of 0.02 still exists when the number of sampling points reaches . The results obtained by using Sobol sequence and Halton sequence converge better than the above two methods. The Sobol sequence sampling method is significantly better in that it can converge to results with less error when using a lower number of sampling points. The above analysis shows that the Sobol sequence sampling method has higher accuracy and better convergence and can be better applied to global sensitivity analysis.
4.2. Global Sensitivity Analysis
4.2.1. Analysis of Position and Velocity Parameters
In the simulation analysis, the 2029 Earth-Moon transfer orbit window is selected, and the spacecraft enters a circular lunar orbit of about 100 km after departing from Earth parking orbit of about 170 km altitude. After preliminary calculations, the TLI moment of one of the Earth-Moon transfer orbits is selected as 2029-07-11 22 : 17 : 1.86, and the moment of the perilune is 2029-07-14 19 : 10 : 00.00. The orbital elements of the TLI point in the geocentric J2000 coordinate system and the perilune in the selenocentric J2000 coordinate system are shown in Table 1.
In practical engineering, the direction of the spacecraft impulse application is mainly applied in the tangential, radial, and normal directions of the spacecraft orbital plane. Therefore, the state of the spacecraft under the orbital coordinate system at the TLI point is more concerned. Converting the spacecraft orbital elements at the TLI point to the orbital coordinate system, the position and velocity of the spacecraft in the orbital coordinate system can be obtained in Table 2.
At the TLI point, both position and velocity changes affect the results, so the sensitivity of the position and velocity parameters need to be analyzed separately. In addition, when the orbital elements are adjusted in the actual project, the position and velocity will change at the same time. So further comprehensive analysis of the sensitivity of position and velocity is required.
First, the sensitivity of the position parameters to the orbital elements of perilune point is considered. To ensure that the calculated target circumlunar orbit still meets the engineering constraints, the circumlunar orbit height deviation is considered to be set within -50 km and 50 km. Then, the initial deviation analysis is conducted in the high-fidelity model, and the minimum deviation value of the spacecraft position in tangential, radial and normal directions is 0.07 km, so the position variation range is set in Table 3.
In Table 3, , , and denotes the tangential, normal, and radial positions of the spacecraft in the orbital coordinate system, respectively. The velocity parameter at the departure moment is kept constant, and the output parameter is chosen as the orbital elements at the perilune. Use the global sensitivity analysis calculation process of Earth-Moon transfer orbit proposed in this paper to analyze the sensitivity of position parameters. First, the convergence of sensitivity results under different sampling points is analyzed. The initial sampling point is set as 200, the maximum number of sampling point is set as 16,000, and the sampling interval is set as 200. The results are illustrated in Figure 5.
(a) First-order sensitivity index
(b) Second-order sensitivity index
(c) Global sensitivity index
Figure 5 shows that in the case of a small number of sampling points, the calculated results of sensitivity index of each order show unstable fluctuation. However, with the increase of the number of sampling points, the calculation results of the sensitivity index of each order gradually tend to be stable. Among them, the calculation results of the global sensitivity coefficient and the second order sensitivity coefficient converge rapidly, and the calculation accuracy has reached 0.002 when the sampling points are about 5000. In contrast, the first order sensitivity index converges more slowly, with a convergence accuracy of about 0.02 when the sampling points reach about 15000, and at the same time, due to the influence of the algorithm, the first-order sensitivity index of some parameters will appear slightly larger than the global sensitivity index or less than 0. To continue to improve the computational accuracy of the first-order sensitivity index, a larger sampling is required, which will lead to a sharp increase in the computational volume. Since the comprehensive influence of the parameters of the Earth-Moon transfer orbit on the results within the whole range of engineering constraints is more concerned and the influence is mainly reflected in the global sensitivity index of the parameters, it is not meaningful to increase the calculation cost to calculate the first order sensitivity index. In order to ensure the accuracy of the calculation results, the number of sampling points is chosen as 15000.
In Table 4, the subscript denotes the orbital elements at the perilune. In the histogram of global sensitivity index, subscripts , , and represent the tangential, normal, and radial positions of position parameters, respectively. Table 4 and Figure 6 show the global sensitivity index of position parameters that the radial position is the most sensitive parameter for the orbital elements at the perilune, and the sensitivity index for all six orbital element reaches more than 98%. Among the other two position parameters, the sensitivity of the normal position to the perilune altitude , eccentricity , and true anomaly is slightly greater than that of the tangential position , while the sensitivity of the tangential position to the orbital Inclination angle , right ascension of ascending node , and argument of perilune are greater than that of the normal position , but both of these parameters have little effect on the results.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
Figure 7 shows the sensitivity results of the combinations between the parameters that the sensitivity results are larger for the combinations of parameters related to the radial direction, while the combination of parameters of tangential and normal position has little effect on the results. Therefore, when considering the effect of position on the results, the radial position at the TLI point is the main parameter needed to be considered.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
Next, the sensitivity of the velocity parameters to the orbital elements of perilune point is considered. Keeping the position state unchanged, the altitude deviation of the circumlunar orbit is also set in the range of -50 km and 50 km. The initial deviation of the velocity parameters at the TLI point is calculated, and the minimum deviation in the three velocity directions is obtained as 0.06 m/s. Therefore, the range of the velocity is set as shown in Table 5.
In Table 5, , , and denotes the tangential, normal, and radial velocity of the spacecraft in the orbital coordinate system, respectively. Calculate the parametric sensitivity of the velocity parameters in the three directions to the orbital elements at the perilune. Table 6 and Figure 8 show that among the three directions of velocity in the orbital coordinate system, the tangential velocity is the most dominant sensitivity parameter with a sensitivity coefficient of 99%, and among the other two velocity parameters, the radial velocity is more sensitive than the normal velocity , so the total sensitivity magnitude is ranked as . Similarly, Figure 9 shows that the sensitivity results are larger for the parameter combination terms related to the tangential velocity . Considering the above results, the tangential velocity at the TLI point is the main influencing parameter to be considered when only the effect of velocity is considered.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
Next, the effect on the results when the position and velocity change simultaneously is considered, and the range of variation of position and velocity in Table 7 is consistent with the above.
The global sensitivity analysis method is used to analyze the six parameters. Table 8 and Figure 10 show that when the position and velocity are changed simultaneously, the main influencing parameters are still the radial position and tangential velocity , and the sensitivity index of the two parameters is similar. The remaining parameters have a small and almost negligible effect on the results, which is also consistent with the results above for the separate analysis of position and velocity.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
In the sensitivity results of the mutual combination of parameters, Figure 11 shows that the higher sensitivity index is also for the combined parameters related to radial position and tangential velocity , and the highest combined parameter is the combination of radial position and tangential velocity . Therefore, the radial position and tangential velocity are the main influencing parameters when velocity and position are considered together.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
4.2.2. Analysis of Orbital Element Parameters
Finally, the above range of position and velocity are converted into the range of orbital elements, and the range of orbital elements is shown in Table 9.
The influence of the orbital elements at the TLI point on the orbital elements at the perilune is further analyzed. The subscript in Table 10 denotes the orbital elements at the perilune, and the subscript denotes the orbital elements at the TLI point. Table 10 and Figure 12 show the calculation results that the sensitivity of the orbital elements at the perilune to the orbital elements at the TLI point varies compared to the parametric sensitivity of the position and velocity. Comparing the global sensitivity index of each parameter, the most sensitive parameter is the eccentricity for the perilune altitude and to which is 0.587674, the most sensitive parameter for the eccentricity is the right ascension of ascending node and to which is 0.545750, the most sensitive parameter for the inclination angle is the argument of perigee and to which is 0.531340, the most sensitive parameter for the right ascension of ascending node is also the argument of perigee and to which is 0.526958, the most sensitive parameter for the argument of perilune is the eccentricity and to which is 0.522090, and the most sensitive parameter for the true anomaly is the right ascension of ascending node and to which is 0.637671.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
The specific sensitivity indexes are ranked in Table 11. It shows that , , and are the main sensitivity elements and the other elements are less sensitivity. The result also shows that although any small change in orbital elements at TLI point will lead to a large error at perilune, but , , and may have a greater impact on the results.
In addition, from the sensitivity calculation results of the combination of parameters in Figure 13, it can be seen that the most sensitive combination of parameters to the perilune altitude is eccentricity and right ascension of ascending node . The most sensitive combination of parameters to the eccentricity is eccentricity and right ascension of ascending node . The most sensitive combination of parameters to the Inclination is argument of perigee and right ascension of ascending node . The most sensitive combination of parameters to the right ascension of ascending node is argument of perigee and right ascension of ascending node . The most sensitive combination of parameters to the argument of perilune is eccentricity and argument of perigee . The most sensitive combination of parameters to the true anomaly is the argument of perigee and right ascension of ascending node . Therefore, when it is necessary to also consider the effect of the combination between orbital elements of the TLI point on the results, it can be determined based on the results of the above sensitivity analysis.
(a) Perilune altitude
(b) Eccentricity
(c) Inclination
(d) Right ascension of ascending node
(e) Argument of perigee
(f) True anomaly
In the simulation analysis, the global sensitivity analysis method is used to calculate the degree of influence of the position, velocity, and the orbital elements of the TLI point on the orbital elements of the perilune, respectively. The sensitivity results obtained using the Sobol method integrally reflect the degree of influence of the parameter on the results in the interval range. Compared with the previous analysis methods, it is more comprehensive and accurate, and we can know the specific proportional size of the influence of each parameter on the results, which verify the validity and feasibility of the calculation process proposed in this paper. From the analysis results, the radial position and tangential velocity are the main influencing parameters of TLI point; , , and are the main sensitivity orbital elements. On the one hand, sensitive parameters need to be controlled with increased precision to ensure the safe implementation of the task. On the other hand, when optimizing the design of the orbit, the adjustments are made mainly for sensitive parameters. For insensitive parameters, the adjustment range can be limited or fixed to constant values directly, thus reducing the number and range of adjusted parameters and optimizing the design process. The related results can provide references for manned Earth-Moon transfer orbit projects.
5. Conclusion
This paper proposes a calculation process of global sensitivity analysis for high-fidelity Earth-Moon transfer orbits based on the Sobol method and discusses the degree of influence of each state parameter at the TLI point on the orbital elements at the perilune. The influence of different sampling methods on the global sensitivity analysis results is also compared, and the conclusions obtained in this paper are as follows: (1)The data points generated by the Sobol sequence exhibit good uniformity of distribution in the sampling area, and this uniformity allows for faster convergence and higher accuracy in the global sensitivity analysis(2)Among the position and velocity parameters at the TLI point, the radial position and tangential velocity are sensitivity parameters and have the greatest influence on the output results, while the rest of the parameters have almost negligible influence on the results. In the actual engineering control and orbit design, it is necessary to focus on these two sensitivity parameters(3)The sensitivity of the orbital elements of the departure moment to the individual orbital elements of the perilune varies in magnitude. When the influence of the number of orbital elements on the results needs to be considered specifically, it can be calculated by the global sensitivity calculation process proposed in this paper and then adjust it for the main influencing parameters based on the results of the sensitivity analysis
Appendix
A. Perturbation Acceleration Terms
A.1. Gravitational Perturbation Acceleration of Body
where is the number of perturbation celestial bodies, is the position vector of the spacecraft relative to the -th celestial body, is the position vector of the -th celestial body relative to the earth, and is the gravity constants of perturbed celestial bodies. In the Earth-Moon transfer orbit problem of this paper, only the solar and lunar perturbations are considered and the relative positions between the stars are obtained from the JPL DE405 ephemeris.
A.2. Nonspherical Perturbation Acceleration of the Earth
The potential function of the Earth’s nonspherical perturbation can be expanded into the following spherical harmonic function in the Earth-fixed coordinate system: where is the Earth equatorial radius; is the Earth’s gravitational constant; is the geocentric distance, geocentric longitude, and geocentric latitude of the spacecraft in the Earth-fixed coordinate system; is the normalized geopotential coefficient; is the normalized Legendre polynomial, and is the order of the Earth’s gravity field model; here, 6 degrees by 6 orders is adopted using the WGS84 model.
Then, the Earth’s nonspherical perturbation acceleration can be expressed as where is the transformation matrix from Earth-fixed coordinate system to geocentric J2000 coordinate system and is the gradient of to Earth solid coordinates .
A.3. Nonspherical Perturbation Acceleration of the Moon
Similar to the method for solving the nonspherical acceleration of the Earth, the spherical harmonic function in the Moon-fixed coordinate system is where is the Moon equatorial radius; is the Moon’s gravitational constant; is the distance, longitude, and latitude of the spacecraft in the Moon-fixed coordinate system; is the gravity coefficient of the moon; is the normalized Legendre polynomial; is the order of the Moon’s gravity field model; here, 6 degrees by 6 orders is adopted using the LP165P model.
Then, the Moon’s nonspherical perturbation acceleration can be expressed as where is the transformation matrix from Moon-fixed coordinate system to the selenocentric J2000 coordinate system and is the gradient of to Moon solid coordinates .
A.4. Perturbation Acceleration of Solar Radiation Pressure
The acceleration of solar radiation pressure perturbation on the spacecraft orbit is where is the reflective property of the spacecraft surface; this paper takes the value of 1.33, is the radiation pressure acting on a blackbody from an astronomical unit of the Sun, the value is , is the area to mass ratio, here the value is 0.004, is the position vector of spacecraft relative to the Sun, and is the shined factor and calculated by using cylindrical shadow model. The algorithm of ground shadow is
The spacecraft is considered to be in the ground shadow when the following conditions are met, and ; otherwise, .
The lunar shadow calculation algorithm is the same as the Earth shadow calculation algorithm and here will not be repeated.
A.5. Perturbation Acceleration of Atmospheric Drag
The atmospheric drag perturbation acceleration is where is the drag coefficient, the nominal value can be taken as 2.2, is cross-section to mass ratio in the pneumatic direction, and here, the value is 0.004, is the atmospheric density, is the velocity of the spacecraft relative to the atmosphere, and it is generally assumed that the atmosphere is synchronized with the rotation of the Earth and can be written as
Atmospheric density models include exponential model, Jacchia model series (J64, J71, J73, and J77), density-temperature model, and international reference atmospheric model. The relative errors of various models are about 15%. The effect of atmospheric drag is generally below 1000 km and has little effect on high orbiting spacecraft. Therefore, this paper only considers the influence of the atmosphere when the altitude of the spacecraft is less than 1000 km; the exponential model is used in the model.
B. Primitive Polynomial of Sobol Sequence
We can use an integer to identify the coefficients of a primitive polynomial
So each primitive polynomial is uniquely specified by its degree together with the number . For example, from and , we can get polynomial . Therefore, the primitive polynomials we used in this paper are as follows: (1)Six dimensions:(2)Twelve dimensions:
In the above equations,represents initial values.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work is supported by the National Natural Science Foundation of China (No. 12072365).