# Global Sensitivity Analysis of Earth-Moon Transfer Orbit Parameters Based on Sobol Method

Authors: Qi Wang; Haiyang Li; Lin Lu; Luyi Yang; Shenquan Wang
Journal: International Journal of Aerospace Engineering
 (2022)

---

## 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.


---

## Body

## 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 ΔvA 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 ΔvB is applied to enter the low lunar orbit. The whole flight process is shown in Figure 1.Figure 1
Schematic diagram of Earth-Moon transfer orbit.It can be seen from Figure1 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(1)r¨=−μErr3+AN+ANSE+ANSM+AR+AD,where μE is the Earth’s gravitational constant, r is the geocentric position vector, AN is gravitational perturbation acceleration of N body, ANSE is the nonspherical perturbation acceleration of the Earth, ANSM is the nonspherical perturbation acceleration of the Moon, AR is the perturbation acceleration of solar radiation pressure, and AD 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.


## 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 ΔvA 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 ΔvB is applied to enter the low lunar orbit. The whole flight process is shown in Figure 1.Figure 1
Schematic diagram of Earth-Moon transfer orbit.It can be seen from Figure1 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(1)r¨=−μErr3+AN+ANSE+ANSM+AR+AD,where μE is the Earth’s gravitational constant, r is the geocentric position vector, AN is gravitational perturbation acceleration of N body, ANSE is the nonspherical perturbation acceleration of the Earth, ANSM is the nonspherical perturbation acceleration of the Moon, AR is the perturbation acceleration of solar radiation pressure, and AD 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 n has a unique R-ary representation which is expressed as
(2)n=nmnm−1⋯n2n1n0R=n0+n1R+n2R2+⋯+nmRm,where m=lnn/lnR and ⋅ represents the rounding operation, n0,n1,n2,⋯,nm⊆N. In Sobol sequence, integers are expanded in base 2, n0,n1,n2,⋯,nm take 0 or 1, that is,
(3)n=n0+2n1+22n2+⋯+2mnm.The primitive polynomial of degreed is defined as
(4)pix=xd+h1xd−1+h2xd−2+⋯+hd−1x+1,where pix is the primitive polynomial of the i-th dimension, d is the degree of freedom of the polynomial, and h1,h2,⋯,hd−1 are 0 or 1. A set of integers m1,m2,⋯ can be obtained by recursion of the coefficient hj of the primitive polynomial; the recursion is as follows:
(5)mk=2h1mk−1⊕22h2mk−2⊕⋯⊕2d−1hd−1mk−d+1⊕2dmk−d⊕mk−d,where ⊕ represents a binary bitwise XOR operation and the initial values m1,m2,⋯,md can be chosen which freely provided that mkk=1,2,⋯ is an odd number less than 2k. Direction number is defined as
(6)vk=mk2kk=1,2,⋯.The Sobol sequence based on the direction number can be obtained by(7)xn=n0v1⊕n1v2⊕⋯,where xn is the n-th element of the generated Sobol sequence and nk is the k-th digit from the right when n is written in binary n=⋯n3n2n12; here, we use the notation ⋅2 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 functionY=fX; X is the state parameter of the TLI point of the Earth-Moon transfer orbit represented by vector
(8)X=xx,xy,xz,vx,vy,vzor
(9)X=rE,eE,iE,ΩE,ωE,fE,where xx,xy,xz represents the position parameters of the spacecraft at the TLI point, vx,vy,vz represents the velocity parameters of the spacecraft at the TLI point, and rE,eE,iE,ΩE,ωE,fE 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 X=X1,X2,⋯,Xd, d is the number of parameters, and Y represents perilune parameter of Earth-Moon transfer orbit. In general, any input parameter space can be converted to a unit cube Xi∈0,1,i=1,2,⋯,d. 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
(10)YX=f0+∑i=1dfiXi+∑i<jdfijXi,Xj+⋯+f1,2,⋯,dX1,X2,⋯,Xd,where f0 is constant, fi is the function of Xi, fij is the function of Xi and Xj, and the rest are analogized. One condition that satisfies this decomposition is
(11)∫01fi1i2⋯isXi1,Xi2,⋯,XisdXk=0,k=i1,i2,⋯,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(12)i1,i2,⋯,is≠j1,j2,⋯,jl,the following equation can be obtained
(13)∫Ωdfi1i2⋯is⋅fj1j2⋯jldX=0.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:
(14)f0=∫ΩdYXdX,fiXi=−f0+∫01⋯∫01YXdX~i,fijXi,Xj=−f0−fixi−fjxj+∫01⋯∫01YXdX~ij,where X∼i,X∼ij represents parameters other than Xi and Xij. Similarly, other higher-order terms can be solved.Further assuming thatYX is square integrable, the total variance is calculated as
(15)V=∫ΩdY2XdX−f02.The deviation ofYX can be obtained from each subterm of the above equation
(16)Vi1,i2,⋯,is=∫01⋯∫01Y2i1,i2,⋯,isXi1,i2,⋯,isdXi1⋯dXis.Then, the total variance can be decomposed into(17)V=∑i=1dVi+∑i<jdVij+⋯+V12⋯d,where
(18)Vi=VXiEX∼iYXi,Vij=VXijEX∼ijYXi,Xj−Vi−Vj.The expressions for the remaining terms are similar.X∼i represents the variables other than Xi; X∼ij is a variable of all parameters except Xi and Xj. The sensitivity index of each order can be obtained by dividing the items in the above formula by the total variance V(19)∑i=1dSi+∑1≤i≤j≤dSi,j+⋯+S1,2,⋯,d=1.The first-order sensitivity indexSi, the second-order sensitivity index Sij, and the global sensitivity index STi are defined as
(20)Si=ViV=VXiEX:iYXiV,Sij=VijV=VXijEX∼ijYXi,XjV−Si−Sj,STi=1−V∼iV=EX∼iVXiYX∼iV,where V⋅ represents variance and E⋅ represents expectation. The first-order sensitivity index Si represents the influence of a single model parameter on the model output, Sij represents the influence of the combination of two model parameters on the model output, and the global sensitivity index STi represents the influence of all the combination of parameters including the i-th model parameter on the model output. For example, consider a function with three parameters. The first-order sensitivity index S1,S2,S3 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. S12 represents the degree of influence on the results under the interaction between the parameters X1 and X2. Similarly, S13,S23 represents the degree of influence on the results under the interaction between the parameter X1,X3 and X2,X3, respectively. The global sensitivity index ST1 of X1 is
(21)ST1=S1+S12+S13+S123.This result reflects the degree of influence of the combined effect of parametersX1 and other parameters on the results. ST2,ST3 are the global sensitivity index of X2,X3, 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 performedN times within a given feasible range to generate a sample matrix N×2d in total. Then, two matrices A and B of size N×d are defined separately, each containing half of the samples(22)A=a11⋮aN1⋯⋱⋯a1d⋮aNd,B=b11⋮bN1⋯⋱⋯b1d⋮bNd(2)
Construct a new matrixABi from matrices A and B, where the i-th column is from matrix B and the remaining d−1 columns are from matrix A:(23)ABi=a11⋮aN1⋯⋯b1i⋱bNi⋯⋯a1d⋮aNd(3)
Substitute the aboved+2 sample matrices into the system model, and d+2 output vectors can be obtained:(24)YA=fA,YB=fB,YABi=fABi,i=1,2,⋯,d(4)
Substitute the above results into the simplified Monte Carlo calculation formula to calculate the variances of each order:(25)VXiEX∼iYXi=1N∑j=1NfBj∗fABij−fAj,(26)EX∼iVXiYX∼i=12N∑j=1NfAj−fABij2,(27)EX∼ijVXiXjYX∼ij=12N∑w=1NfABiw−fABjw2Equations (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 Figure2.Figure 2
Global sensitivity analysis calculation process for Earth-Moon transfer orbit.


## 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 n has a unique R-ary representation which is expressed as
(2)n=nmnm−1⋯n2n1n0R=n0+n1R+n2R2+⋯+nmRm,where m=lnn/lnR and ⋅ represents the rounding operation, n0,n1,n2,⋯,nm⊆N. In Sobol sequence, integers are expanded in base 2, n0,n1,n2,⋯,nm take 0 or 1, that is,
(3)n=n0+2n1+22n2+⋯+2mnm.The primitive polynomial of degreed is defined as
(4)pix=xd+h1xd−1+h2xd−2+⋯+hd−1x+1,where pix is the primitive polynomial of the i-th dimension, d is the degree of freedom of the polynomial, and h1,h2,⋯,hd−1 are 0 or 1. A set of integers m1,m2,⋯ can be obtained by recursion of the coefficient hj of the primitive polynomial; the recursion is as follows:
(5)mk=2h1mk−1⊕22h2mk−2⊕⋯⊕2d−1hd−1mk−d+1⊕2dmk−d⊕mk−d,where ⊕ represents a binary bitwise XOR operation and the initial values m1,m2,⋯,md can be chosen which freely provided that mkk=1,2,⋯ is an odd number less than 2k. Direction number is defined as
(6)vk=mk2kk=1,2,⋯.The Sobol sequence based on the direction number can be obtained by(7)xn=n0v1⊕n1v2⊕⋯,where xn is the n-th element of the generated Sobol sequence and nk is the k-th digit from the right when n is written in binary n=⋯n3n2n12; here, we use the notation ⋅2 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 functionY=fX; X is the state parameter of the TLI point of the Earth-Moon transfer orbit represented by vector
(8)X=xx,xy,xz,vx,vy,vzor
(9)X=rE,eE,iE,ΩE,ωE,fE,where xx,xy,xz represents the position parameters of the spacecraft at the TLI point, vx,vy,vz represents the velocity parameters of the spacecraft at the TLI point, and rE,eE,iE,ΩE,ωE,fE 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 X=X1,X2,⋯,Xd, d is the number of parameters, and Y represents perilune parameter of Earth-Moon transfer orbit. In general, any input parameter space can be converted to a unit cube Xi∈0,1,i=1,2,⋯,d. 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
(10)YX=f0+∑i=1dfiXi+∑i<jdfijXi,Xj+⋯+f1,2,⋯,dX1,X2,⋯,Xd,where f0 is constant, fi is the function of Xi, fij is the function of Xi and Xj, and the rest are analogized. One condition that satisfies this decomposition is
(11)∫01fi1i2⋯isXi1,Xi2,⋯,XisdXk=0,k=i1,i2,⋯,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(12)i1,i2,⋯,is≠j1,j2,⋯,jl,the following equation can be obtained
(13)∫Ωdfi1i2⋯is⋅fj1j2⋯jldX=0.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:
(14)f0=∫ΩdYXdX,fiXi=−f0+∫01⋯∫01YXdX~i,fijXi,Xj=−f0−fixi−fjxj+∫01⋯∫01YXdX~ij,where X∼i,X∼ij represents parameters other than Xi and Xij. Similarly, other higher-order terms can be solved.Further assuming thatYX is square integrable, the total variance is calculated as
(15)V=∫ΩdY2XdX−f02.The deviation ofYX can be obtained from each subterm of the above equation
(16)Vi1,i2,⋯,is=∫01⋯∫01Y2i1,i2,⋯,isXi1,i2,⋯,isdXi1⋯dXis.Then, the total variance can be decomposed into(17)V=∑i=1dVi+∑i<jdVij+⋯+V12⋯d,where
(18)Vi=VXiEX∼iYXi,Vij=VXijEX∼ijYXi,Xj−Vi−Vj.The expressions for the remaining terms are similar.X∼i represents the variables other than Xi; X∼ij is a variable of all parameters except Xi and Xj. The sensitivity index of each order can be obtained by dividing the items in the above formula by the total variance V(19)∑i=1dSi+∑1≤i≤j≤dSi,j+⋯+S1,2,⋯,d=1.The first-order sensitivity indexSi, the second-order sensitivity index Sij, and the global sensitivity index STi are defined as
(20)Si=ViV=VXiEX:iYXiV,Sij=VijV=VXijEX∼ijYXi,XjV−Si−Sj,STi=1−V∼iV=EX∼iVXiYX∼iV,where V⋅ represents variance and E⋅ represents expectation. The first-order sensitivity index Si represents the influence of a single model parameter on the model output, Sij represents the influence of the combination of two model parameters on the model output, and the global sensitivity index STi represents the influence of all the combination of parameters including the i-th model parameter on the model output. For example, consider a function with three parameters. The first-order sensitivity index S1,S2,S3 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. S12 represents the degree of influence on the results under the interaction between the parameters X1 and X2. Similarly, S13,S23 represents the degree of influence on the results under the interaction between the parameter X1,X3 and X2,X3, respectively. The global sensitivity index ST1 of X1 is
(21)ST1=S1+S12+S13+S123.This result reflects the degree of influence of the combined effect of parametersX1 and other parameters on the results. ST2,ST3 are the global sensitivity index of X2,X3, 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 performedN times within a given feasible range to generate a sample matrix N×2d in total. Then, two matrices A and B of size N×d are defined separately, each containing half of the samples(22)A=a11⋮aN1⋯⋱⋯a1d⋮aNd,B=b11⋮bN1⋯⋱⋯b1d⋮bNd(2)
Construct a new matrixABi from matrices A and B, where the i-th column is from matrix B and the remaining d−1 columns are from matrix A:(23)ABi=a11⋮aN1⋯⋯b1i⋱bNi⋯⋯a1d⋮aNd(3)
Substitute the aboved+2 sample matrices into the system model, and d+2 output vectors can be obtained:(24)YA=fA,YB=fB,YABi=fABi,i=1,2,⋯,d(4)
Substitute the above results into the simplified Monte Carlo calculation formula to calculate the variances of each order:(25)VXiEX∼iYXi=1N∑j=1NfBj∗fABij−fAj,(26)EX∼iVXiYX∼i=12N∑j=1NfAj−fABij2,(27)EX∼ijVXiXjYX∼ij=12N∑w=1NfABiw−fABjw2Equations (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 Figure2.Figure 2
Global sensitivity analysis calculation process for Earth-Moon transfer orbit.


## 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 N equal parts according to the number of samples to obtain N 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. Figure3 shows the distribution of sampling points in the sampling area generated by different sampling methods.Figure 3
Distribution of data points generated by different sampling methods: (a) 100 sampling points, (b) 1000 sampling points, and (c) 4000 sampling points.
(a)(b)(c)Then, the effect of different sampling methods on the results of parametric sensitivity analysis is further analyzed. Selecting the simple functiony=x1+x2, the values of x1 and x2 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 x1 and x2 on the output results, and the error of the calculated results relative to the reference value is counted at different numbers of sampling points.Figure4 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 104. 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.Figure 4
Effect of different sampling methods on the results of 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 Table1.Table 1
Orbital elements of Earth-Moon transfer orbit.
Orbit typeGeocentric distance/perilune distancer (m)EccentricityeInclinationi (deg)Right Ascension of ascending nodeΩ (deg)Argument of perigee/periluneω (deg)True anomalyf (deg)LEO6548142.30.9745228.0142.8199.890LLO1831332.71.4945148.42171.0223.93148.42In 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 Table2.Table 2
Spacecraft departure status in orbital coordinate system.
Tangential (x)Normal (y)Radial (z)Position (m)00-6548142.3Velocity (m/s)10963.269100At 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 Table3.Table 3
Range of values for TLI position.
xxxyxzrmin (m)-70-70-6548212.3rmax (m)7070-6548072.3In Table3, xx, xy, and xz 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.Figure 5
Effect of the number of sampling points on sensitivity results.
(a)
First-order sensitivity index(b)
Second-order sensitivity index(c)
Global sensitivity indexFigure5 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 Table4, the subscript M denotes the orbital elements at the perilune. In the histogram of global sensitivity index, subscripts xx, xy, and xz 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 xz 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 xy to the perilune altitude rM, eccentricity eM, and true anomaly fM is slightly greater than that of the tangential position xx, while the sensitivity of the tangential position xx to the orbital Inclination angle iM, right ascension of ascending node ΩM, and argument of perilune ωM are greater than that of the normal position xy, but both of these parameters have little effect on the results.Table 4
Sensitivity index of position parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexxx0.0005280.0022850.0087690.0078340.0021340.003946xy0.0008560.0042770.0027490.002870.0018660.009647xz0.9986040.9934230.9884640.9892890.9959910.986409Figure 6
Histogram of the global sensitivity index of the position parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyFigure7 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 xx and normal position xy has little effect on the results. Therefore, when considering the effect of position on the results, the radial position xz at the TLI point is the main parameter needed to be considered.Figure 7
Sensitivity index for combinations of position parameters.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyNext, 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 Table5.Table 5
The range of velocity of the TLI point.
vxvyvzvmin (m/s)10963.2091-0.06-0.06vmax (m/s)10963.32910.060.06In Table5, vx, vy, and vz 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 vx is the most dominant sensitivity parameter with a sensitivity coefficient of 99%, and among the other two velocity parameters, the radial velocity vz is more sensitive than the normal velocity vy, so the total sensitivity magnitude is ranked as tangential>radial>normal. Similarly, Figure 9 shows that the sensitivity results are larger for the parameter combination terms related to the tangential velocity vx. 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.Table 6
Sensitivity index of velocity parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexvx0.9994200.9976510.9910710.9920160.9978110.995974vy0.0000040.0000020.0000130.0000130.0000080.000045vz0.0005390.0023030.0088720.0079150.0021290.003988Figure 8
Histogram of the global sensitivity index of the velocity parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyFigure 9
Sensitivity index for combinations of velocity parameters.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyNext, the effect on the results when the position and velocity change simultaneously is considered, and the range of variation of position and velocity in Table7 is consistent with the above.Table 7
The range of position and velocity values.
xxxyxzvxvyvzrmin/vmin-70-70-6548212.310963.2091-0.06-0.06rmax/vmax7070-6548072.310963.32910.060.06The global sensitivity analysis method is used to analyze the six parameters. Table8 and Figure 10 show that when the position and velocity are changed simultaneously, the main influencing parameters are still the radial position xz and tangential velocity vx, 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.Table 8
Sensitivity index for position and velocity parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexxx0.0002610.0011300.0043810.0039070.0010600.001983xy0.0004230.0021140.0013720.0014310.0009300.004844xz0.4936580.4912200.4934570.4937530.4952230.492900vx0.5054090.5047450.4963680.4972700.5018730.498328vy0.0000020.0000100.0000070.0000060.0000040.000023vz0.0002720.0011650.0044440.0039640.0010700.002006Figure 10
Histogram of the global sensitivity index of the position and velocity parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyIn the sensitivity results of the mutual combination of parameters, Figure11 shows that the higher sensitivity index is also for the combined parameters related to radial position xz and tangential velocity vx, and the highest combined parameter is the combination of radial position xz and tangential velocity vx. Therefore, the radial position xz and tangential velocity vx are the main influencing parameters when velocity and position are considered together.Figure 11
Sensitivity index for combination of position and velocity parameters.
(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 Table9.Table 9
The range of orbital elements values.
rE (m)eEiE (deg)ΩE (deg)ωE (deg)fE (deg)Min6548072.30.97451469627.999505142.796823199.884422-0.0006056457Max6548212.30.97451570528.000511142.798823199.8862010.0006056915The influence of the orbital elements at the TLI point on the orbital elements at the perilune is further analyzed. The subscriptM in Table 10 denotes the orbital elements at the perilune, and the subscript E 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 eE for the perilune altitude rM and eE to rMwhich is 0.587674, the most sensitive parameter for the eccentricity eM is the right ascension of ascending node ΩE and ΩE to eMwhich is 0.545750, the most sensitive parameter for the inclination angle iM is the argument of perigee ωE and ωE to iMwhich is 0.531340, the most sensitive parameter for the right ascension of ascending node ΩM is also the argument of perigee ωE and ωE to ΩMwhich is 0.526958, the most sensitive parameter for the argument of perilune ωM is the eccentricity eE and eE to ωMwhich is 0.522090, and the most sensitive parameter for the true anomaly fM is the right ascension of ascending node ΩE and ΩE to fMwhich is 0.637671.Table 10
Sensitivity index for orbital elements.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexrE0.0126760.0039060.0126250.0132030.0231600.007696eE0.5876740.2533520.1901840.2116080.5220900.150813iE0.0026190.0054400.0025810.0031520.0050130.007428ΩE0.2788660.5457500.2536470.2346700.0827060.637671ωE0.1016760.1854990.5313400.5269580.3461870.189492fE0.0164700.0058550.0096000.0103770.0207820.007005Figure 12
Histogram of the global sensitivity index of orbital elements at the TLI point to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyThe specific sensitivity indexes are ranked in Table11. It shows that eE, ΩE, and ωE 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 eE, ΩE, and ωE may have a greater impact on the results.Table 11
Sorting of parametric sensitivity analysis results.In addition, from the sensitivity calculation results of the combination of parameters in Figure13, it can be seen that the most sensitive combination of parameters to the perilune altitude rM is eccentricity eE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the eccentricity eM is eccentricity eE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the Inclination iM is argument of perigee ωE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the right ascension of ascending node ΩM is argument of perigee ωE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the argument of perilune ωM is eccentricity eE and argument of perigee ωE. The most sensitive combination of parameters to the true anomaly fM is the argument of perigee ωE and right ascension of ascending node ΩE. 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.Figure 13
Sensitivity index for combination of orbital elements.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyIn 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;eE, ΩE, and ωE 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.


## 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 N equal parts according to the number of samples to obtain N 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. Figure3 shows the distribution of sampling points in the sampling area generated by different sampling methods.Figure 3
Distribution of data points generated by different sampling methods: (a) 100 sampling points, (b) 1000 sampling points, and (c) 4000 sampling points.
(a)(b)(c)Then, the effect of different sampling methods on the results of parametric sensitivity analysis is further analyzed. Selecting the simple functiony=x1+x2, the values of x1 and x2 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 x1 and x2 on the output results, and the error of the calculated results relative to the reference value is counted at different numbers of sampling points.Figure4 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 104. 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.Figure 4
Effect of different sampling methods on the results of 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 Table1.Table 1
Orbital elements of Earth-Moon transfer orbit.
Orbit typeGeocentric distance/perilune distancer (m)EccentricityeInclinationi (deg)Right Ascension of ascending nodeΩ (deg)Argument of perigee/periluneω (deg)True anomalyf (deg)LEO6548142.30.9745228.0142.8199.890LLO1831332.71.4945148.42171.0223.93148.42In 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 Table2.Table 2
Spacecraft departure status in orbital coordinate system.
Tangential (x)Normal (y)Radial (z)Position (m)00-6548142.3Velocity (m/s)10963.269100At 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 Table3.Table 3
Range of values for TLI position.
xxxyxzrmin (m)-70-70-6548212.3rmax (m)7070-6548072.3In Table3, xx, xy, and xz 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.Figure 5
Effect of the number of sampling points on sensitivity results.
(a)
First-order sensitivity index(b)
Second-order sensitivity index(c)
Global sensitivity indexFigure5 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 Table4, the subscript M denotes the orbital elements at the perilune. In the histogram of global sensitivity index, subscripts xx, xy, and xz 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 xz 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 xy to the perilune altitude rM, eccentricity eM, and true anomaly fM is slightly greater than that of the tangential position xx, while the sensitivity of the tangential position xx to the orbital Inclination angle iM, right ascension of ascending node ΩM, and argument of perilune ωM are greater than that of the normal position xy, but both of these parameters have little effect on the results.Table 4
Sensitivity index of position parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexxx0.0005280.0022850.0087690.0078340.0021340.003946xy0.0008560.0042770.0027490.002870.0018660.009647xz0.9986040.9934230.9884640.9892890.9959910.986409Figure 6
Histogram of the global sensitivity index of the position parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyFigure7 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 xx and normal position xy has little effect on the results. Therefore, when considering the effect of position on the results, the radial position xz at the TLI point is the main parameter needed to be considered.Figure 7
Sensitivity index for combinations of position parameters.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyNext, 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 Table5.Table 5
The range of velocity of the TLI point.
vxvyvzvmin (m/s)10963.2091-0.06-0.06vmax (m/s)10963.32910.060.06In Table5, vx, vy, and vz 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 vx is the most dominant sensitivity parameter with a sensitivity coefficient of 99%, and among the other two velocity parameters, the radial velocity vz is more sensitive than the normal velocity vy, so the total sensitivity magnitude is ranked as tangential>radial>normal. Similarly, Figure 9 shows that the sensitivity results are larger for the parameter combination terms related to the tangential velocity vx. 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.Table 6
Sensitivity index of velocity parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexvx0.9994200.9976510.9910710.9920160.9978110.995974vy0.0000040.0000020.0000130.0000130.0000080.000045vz0.0005390.0023030.0088720.0079150.0021290.003988Figure 8
Histogram of the global sensitivity index of the velocity parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyFigure 9
Sensitivity index for combinations of velocity parameters.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyNext, the effect on the results when the position and velocity change simultaneously is considered, and the range of variation of position and velocity in Table7 is consistent with the above.Table 7
The range of position and velocity values.
xxxyxzvxvyvzrmin/vmin-70-70-6548212.310963.2091-0.06-0.06rmax/vmax7070-6548072.310963.32910.060.06The global sensitivity analysis method is used to analyze the six parameters. Table8 and Figure 10 show that when the position and velocity are changed simultaneously, the main influencing parameters are still the radial position xz and tangential velocity vx, 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.Table 8
Sensitivity index for position and velocity parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexxx0.0002610.0011300.0043810.0039070.0010600.001983xy0.0004230.0021140.0013720.0014310.0009300.004844xz0.4936580.4912200.4934570.4937530.4952230.492900vx0.5054090.5047450.4963680.4972700.5018730.498328vy0.0000020.0000100.0000070.0000060.0000040.000023vz0.0002720.0011650.0044440.0039640.0010700.002006Figure 10
Histogram of the global sensitivity index of the position and velocity parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyIn the sensitivity results of the mutual combination of parameters, Figure11 shows that the higher sensitivity index is also for the combined parameters related to radial position xz and tangential velocity vx, and the highest combined parameter is the combination of radial position xz and tangential velocity vx. Therefore, the radial position xz and tangential velocity vx are the main influencing parameters when velocity and position are considered together.Figure 11
Sensitivity index for combination of position and velocity parameters.
(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 Table9.Table 9
The range of orbital elements values.
rE (m)eEiE (deg)ΩE (deg)ωE (deg)fE (deg)Min6548072.30.97451469627.999505142.796823199.884422-0.0006056457Max6548212.30.97451570528.000511142.798823199.8862010.0006056915The influence of the orbital elements at the TLI point on the orbital elements at the perilune is further analyzed. The subscriptM in Table 10 denotes the orbital elements at the perilune, and the subscript E 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 eE for the perilune altitude rM and eE to rMwhich is 0.587674, the most sensitive parameter for the eccentricity eM is the right ascension of ascending node ΩE and ΩE to eMwhich is 0.545750, the most sensitive parameter for the inclination angle iM is the argument of perigee ωE and ωE to iMwhich is 0.531340, the most sensitive parameter for the right ascension of ascending node ΩM is also the argument of perigee ωE and ωE to ΩMwhich is 0.526958, the most sensitive parameter for the argument of perilune ωM is the eccentricity eE and eE to ωMwhich is 0.522090, and the most sensitive parameter for the true anomaly fM is the right ascension of ascending node ΩE and ΩE to fMwhich is 0.637671.Table 10
Sensitivity index for orbital elements.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexrE0.0126760.0039060.0126250.0132030.0231600.007696eE0.5876740.2533520.1901840.2116080.5220900.150813iE0.0026190.0054400.0025810.0031520.0050130.007428ΩE0.2788660.5457500.2536470.2346700.0827060.637671ωE0.1016760.1854990.5313400.5269580.3461870.189492fE0.0164700.0058550.0096000.0103770.0207820.007005Figure 12
Histogram of the global sensitivity index of orbital elements at the TLI point to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyThe specific sensitivity indexes are ranked in Table11. It shows that eE, ΩE, and ωE 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 eE, ΩE, and ωE may have a greater impact on the results.Table 11
Sorting of parametric sensitivity analysis results.In addition, from the sensitivity calculation results of the combination of parameters in Figure13, it can be seen that the most sensitive combination of parameters to the perilune altitude rM is eccentricity eE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the eccentricity eM is eccentricity eE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the Inclination iM is argument of perigee ωE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the right ascension of ascending node ΩM is argument of perigee ωE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the argument of perilune ωM is eccentricity eE and argument of perigee ωE. The most sensitive combination of parameters to the true anomaly fM is the argument of perigee ωE and right ascension of ascending node ΩE. 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.Figure 13
Sensitivity index for combination of orbital elements.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyIn 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;eE, ΩE, and ωE 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.


## 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 Table1.Table 1
Orbital elements of Earth-Moon transfer orbit.
Orbit typeGeocentric distance/perilune distancer (m)EccentricityeInclinationi (deg)Right Ascension of ascending nodeΩ (deg)Argument of perigee/periluneω (deg)True anomalyf (deg)LEO6548142.30.9745228.0142.8199.890LLO1831332.71.4945148.42171.0223.93148.42In 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 Table2.Table 2
Spacecraft departure status in orbital coordinate system.
Tangential (x)Normal (y)Radial (z)Position (m)00-6548142.3Velocity (m/s)10963.269100At 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 Table3.Table 3
Range of values for TLI position.
xxxyxzrmin (m)-70-70-6548212.3rmax (m)7070-6548072.3In Table3, xx, xy, and xz 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.Figure 5
Effect of the number of sampling points on sensitivity results.
(a)
First-order sensitivity index(b)
Second-order sensitivity index(c)
Global sensitivity indexFigure5 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 Table4, the subscript M denotes the orbital elements at the perilune. In the histogram of global sensitivity index, subscripts xx, xy, and xz 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 xz 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 xy to the perilune altitude rM, eccentricity eM, and true anomaly fM is slightly greater than that of the tangential position xx, while the sensitivity of the tangential position xx to the orbital Inclination angle iM, right ascension of ascending node ΩM, and argument of perilune ωM are greater than that of the normal position xy, but both of these parameters have little effect on the results.Table 4
Sensitivity index of position parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexxx0.0005280.0022850.0087690.0078340.0021340.003946xy0.0008560.0042770.0027490.002870.0018660.009647xz0.9986040.9934230.9884640.9892890.9959910.986409Figure 6
Histogram of the global sensitivity index of the position parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyFigure7 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 xx and normal position xy has little effect on the results. Therefore, when considering the effect of position on the results, the radial position xz at the TLI point is the main parameter needed to be considered.Figure 7
Sensitivity index for combinations of position parameters.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyNext, 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 Table5.Table 5
The range of velocity of the TLI point.
vxvyvzvmin (m/s)10963.2091-0.06-0.06vmax (m/s)10963.32910.060.06In Table5, vx, vy, and vz 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 vx is the most dominant sensitivity parameter with a sensitivity coefficient of 99%, and among the other two velocity parameters, the radial velocity vz is more sensitive than the normal velocity vy, so the total sensitivity magnitude is ranked as tangential>radial>normal. Similarly, Figure 9 shows that the sensitivity results are larger for the parameter combination terms related to the tangential velocity vx. 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.Table 6
Sensitivity index of velocity parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexvx0.9994200.9976510.9910710.9920160.9978110.995974vy0.0000040.0000020.0000130.0000130.0000080.000045vz0.0005390.0023030.0088720.0079150.0021290.003988Figure 8
Histogram of the global sensitivity index of the velocity parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyFigure 9
Sensitivity index for combinations of velocity parameters.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyNext, the effect on the results when the position and velocity change simultaneously is considered, and the range of variation of position and velocity in Table7 is consistent with the above.Table 7
The range of position and velocity values.
xxxyxzvxvyvzrmin/vmin-70-70-6548212.310963.2091-0.06-0.06rmax/vmax7070-6548072.310963.32910.060.06The global sensitivity analysis method is used to analyze the six parameters. Table8 and Figure 10 show that when the position and velocity are changed simultaneously, the main influencing parameters are still the radial position xz and tangential velocity vx, 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.Table 8
Sensitivity index for position and velocity parameters.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexxx0.0002610.0011300.0043810.0039070.0010600.001983xy0.0004230.0021140.0013720.0014310.0009300.004844xz0.4936580.4912200.4934570.4937530.4952230.492900vx0.5054090.5047450.4963680.4972700.5018730.498328vy0.0000020.0000100.0000070.0000060.0000040.000023vz0.0002720.0011650.0044440.0039640.0010700.002006Figure 10
Histogram of the global sensitivity index of the position and velocity parameters to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyIn the sensitivity results of the mutual combination of parameters, Figure11 shows that the higher sensitivity index is also for the combined parameters related to radial position xz and tangential velocity vx, and the highest combined parameter is the combination of radial position xz and tangential velocity vx. Therefore, the radial position xz and tangential velocity vx are the main influencing parameters when velocity and position are considered together.Figure 11
Sensitivity index for combination of position and velocity parameters.
(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 Table9.Table 9
The range of orbital elements values.
rE (m)eEiE (deg)ΩE (deg)ωE (deg)fE (deg)Min6548072.30.97451469627.999505142.796823199.884422-0.0006056457Max6548212.30.97451570528.000511142.798823199.8862010.0006056915The influence of the orbital elements at the TLI point on the orbital elements at the perilune is further analyzed. The subscriptM in Table 10 denotes the orbital elements at the perilune, and the subscript E 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 eE for the perilune altitude rM and eE to rMwhich is 0.587674, the most sensitive parameter for the eccentricity eM is the right ascension of ascending node ΩE and ΩE to eMwhich is 0.545750, the most sensitive parameter for the inclination angle iM is the argument of perigee ωE and ωE to iMwhich is 0.531340, the most sensitive parameter for the right ascension of ascending node ΩM is also the argument of perigee ωE and ωE to ΩMwhich is 0.526958, the most sensitive parameter for the argument of perilune ωM is the eccentricity eE and eE to ωMwhich is 0.522090, and the most sensitive parameter for the true anomaly fM is the right ascension of ascending node ΩE and ΩE to fMwhich is 0.637671.Table 10
Sensitivity index for orbital elements.
Symbols of parametersrMeMiMΩMωMfMGlobal sensitivity indexrE0.0126760.0039060.0126250.0132030.0231600.007696eE0.5876740.2533520.1901840.2116080.5220900.150813iE0.0026190.0054400.0025810.0031520.0050130.007428ΩE0.2788660.5457500.2536470.2346700.0827060.637671ωE0.1016760.1854990.5313400.5269580.3461870.189492fE0.0164700.0058550.0096000.0103770.0207820.007005Figure 12
Histogram of the global sensitivity index of orbital elements at the TLI point to the orbital elements at the perilune.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyThe specific sensitivity indexes are ranked in Table11. It shows that eE, ΩE, and ωE 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 eE, ΩE, and ωE may have a greater impact on the results.Table 11
Sorting of parametric sensitivity analysis results.In addition, from the sensitivity calculation results of the combination of parameters in Figure13, it can be seen that the most sensitive combination of parameters to the perilune altitude rM is eccentricity eE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the eccentricity eM is eccentricity eE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the Inclination iM is argument of perigee ωE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the right ascension of ascending node ΩM is argument of perigee ωE and right ascension of ascending node ΩE. The most sensitive combination of parameters to the argument of perilune ωM is eccentricity eE and argument of perigee ωE. The most sensitive combination of parameters to the true anomaly fM is the argument of perigee ωE and right ascension of ascending node ΩE. 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.Figure 13
Sensitivity index for combination of orbital elements.
(a)
Perilune altitude(b)
Eccentricity(c)
Inclination(d)
Right ascension of ascending node(e)
Argument of perigee(f)
True anomalyIn 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;eE, ΩE, and ωE 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



---

*Source: 6587890-2022-05-29.xml*