Onboard satellite visibility prediction using metamodeling based framework

Aerospace Science and Technology 94 (2019) 105377

Contents lists available at ScienceDirect

Aerospace Science and Technology www.elsevier.com/locate/aescte

Onboard satellite visibility prediction using metamodeling based framework Xinwei Wang, Chao Han ∗ , Pengbin Yang, Xiucong Sun School of Astronautics, Beihang University, Beijing, 100191, China

a r t i c l e

i n f o

Article history: Received 12 November 2018 Received in revised form 1 July 2019 Accepted 2 September 2019 Available online 4 September 2019 Keywords: Onboard autonomous mission Satellite visibility Roots-finding Metamodeling Self-adaptive interpolation

a b s t r a c t Satellite autonomous systems are employed to address complex space applications through onboard data processing and mission planning. To take advantage of onboard autonomous systems, rapid onboard satellite visibility predictions are necessary for certain decision-making missions, including Earth observation resource allocation and satellite data transmission. We consider this visibility prediction process as a roots-finding problem for a multiple hump function, and design a metamodeling-based framework with a self-adaptive interpolation method. Metamodels are developed as surrogates for visibility prediction functions to reduce expensive computational costs. Our proposed framework has a broad range of applications for all orbital types and orbit propagators. We conduct the experiments using different metamodeling techniques, radial basis functions, Kriging, and support vector regression based upon real China’s satellites. Numerical simulations indicate that the proposed framework outperforms existing interpolation methods, efficiently reducing the onboard computational cost. © 2019 Elsevier Masson SAS. All rights reserved.

1. Introduction Onboard satellite autonomous systems have significant potential for space missions, since they are capable of higher performance and have better fault tolerance. As the complexity of future space missions increases, onboard autonomous systems will become more crucial. For the first mission of NASA’s New Millennium Program Deep Space 1, the Jet Propulsion Laboratory developed several onboard autonomous systems regarded as the new standard for future missions and plans [1]. The Air Force Research Laboratory initiated the TechSat-21 program to serve as an autonomous satellite constellation flight demonstration [2]. The spacecraft Earth Observing-1 applied autonomous software technologies for scientific analysis and mission planning [3]. Following NASA’s example, the European Space Agency launched satellite PROBA-1 [4], whose process of data transmission, resource management, and payload operations is controlled by an onboard autonomous system. Meanwhile, NASA created a planning and scheduling platform EUROPA [5], which has been used for various missions and demonstrations. A detailed technical report [6] described spacecraft autonomy challenges for next-generation space missions, which could be addressed by these systems. Through the

*

Corresponding author. E-mail address: [email protected] (C. Han).

https://doi.org/10.1016/j.ast.2019.105377 1270-9638/© 2019 Elsevier Masson SAS. All rights reserved.

above examples, the feasibility and superiority of onboard satellite autonomous systems have been demonstrated. To take advantage of onboard autonomous systems, rapid onboard satellite visibility predictions are necessary for certain decision-making missions, including Earth observation resource allocation and satellite data transmission. In regards to Earth observation, a growing number of orbiting satellites has raised observation demands [7]. With onboard autonomous systems, the satellite will be allocated to observe assigned ground sites, only under the conditions that the satellite-to-site visibilities are efficiently obtained in advance [8–11]. As during autonomous operations of a satellite constellation, data has to be transmitted within a specific visibility [12–14]. Therefore the satellites can only execute observation mission or data transmission within given visible intervals, indicating that rapid onboard satellite visibility prediction is of great significance. An intuitive brute force method has been used as a practical tool to generate satellite visibility [15]. With constant interpolation steps set by the user, the brute force method checks every minute whether the satellite has access to the associated satellite or site. In order to obtain visible intervals of satellite-to-site efficiently, several fast methods have been designed. On the basis of the Keplerian orbits, Escobal [16] derived a solution from a transcendental equation by an approximate technique. Mai and Palmer [17] further designed a fast access prediction method with coarse search and refinement phases, including two-body motion, secular perturbation, and atmospheric drag. Focusing on a circu-

2

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

lar low Earth orbit (LEO) satellite, Ali et al. [18] presented an algorithm to quickly determine the visibility-time function, where the ground trace of the satellite was simplified as a great-circle arc. Palmer and Mai [19] designed a fast method for LEO satellites, combining the perturbations and atmospheric drag. For the low-eccentricity orbit satellites around oblate-Earth, Lawton [15] developed a rapid algorithm with Fourier series. Notice that these methods are limited to certain orbital types. For all orbital eccentricities and perturbed satellite motions, Alfano et al. [20] proposed a fast interval prediction method via parabolic blending, while the waveform of the access function was described at a fixed time step of 250 seconds. To improve algorithmic efficiency, the Hermite interpolation method was proposed to design and analyze a constellation configuration optimization model [21]. As a continuation of this method [21], Han et al. [22] further improved the efficiency of the Hermite interpolation method, in which the approximation error was utilized to optimize the selection of interpolation points. However, it is observed that the improved method sacrifices some accuracy to guarantee algorithmic efficiency, especially for long-term visibility predictions. On the other hand, the satellite-to-satellite visibility calculation can be also regarded as a visibility prediction problem, which has attracted less attention. Considering low-eccentricity satellite orbits, Lawton [15] developed a fast Fourier-transform computation method for satellite-to-satellite in-view periods. Sun et al. [23] considered the satellite-to-satellite communication problem, where the line of sight between the satellites is obstructed by the Earth. A curve fitting method was introduced to approximate the waveforms of the visibility functions, though the prediction conditions with two inequalities are complicated. As shown in Fig. 1, we illustrate the visibility prediction functions for satellite-to-site and satellite-to-satellite, which will be explained in the next section. Clearly the prediction functions with multiple humps are erratic and without fixed periods. Therefore, the analytical expressions for satellite-to-site and satellite-tosatellite visibility prediction functions are difficult to find. In order to find the roots of a prediction function efficiently, metamodeling techniques are introduced in this work. Metamodeling techniques have been developed for decades in different fields of computer science, statistics, and engineering [24–27]. These metamodels are typically designed as surrogates for a time consuming process. Considering the visibility prediction function is complicated and computationally time consuming, three types metamodeling techniques – radial basis functions (RBFs), Kriging, and support vector regression (SVR) – are introduced to approximate intervals prediction function with satisfied precision. In this paper, we propose an metamodeling-based framework for rapid onboard visibility predictions that combine preprocessing, metamodels, and self-adaptive interpolation methods. The main contributions of this paper are threefold: (1) a general mathematical model of satellite-to-site and satellite-to-satellite visibility predictions is established; (2) for the first time, these metamodeling techniques are comprehensively introduced for onboard satellite visibility predictions; (3) we compare the proposed framework with traditional brute force and other interpolation algorithms. The remainder of this paper is structured as follows. In Section 2, the basic problem statement of visibility predictions and multiple hump function roots-finding is presented. The proposed framework, which includes preprocessing, metamodeling techniques, and a self-adaptive interpolation process, is detailed in Section 3. Section 4 contains a series of computational results which demonstrate the performance of the proposed framework. We then summarize our work in Section 5.

Fig. 1. Examples of visibility prediction as multiple hump function.

2. Problem statement 2.1. Satellite-to-site visibility prediction The definitions of satellite-to-site visibility function and corresponding threshold [17,21,22] vary in literature, although the ground elevation angle is always the key factor. As seen in Fig. 2, the Earth is considered as an ellipsoid, whose eccentricity and equatorial radius are denoted as e E and R e , and o is the Earth center. r sat and r site are the position vectors of the satellite and ground site respectively, r = r sat − r site , θ is the ground elevation angle, and α is the minimum elevation angle at which the satellite has access to the ground site. Following [22], we define the visibility function as

v (t ) =

r · r Lsite r r Lsite 

(1)

where r Lsite is the vector of local vertical direction. Operator · represents the dot product of two vectors and  ·  is the magnitude of a vector. Notice that the vector r Lsite can be determined by geographical method and given as input parameters. When the Earth is considered as round, r Lsite = r site . Notice that the threshold value λ is given as an input parameter and carefully selected according to different mission requirements. This prediction problem becomes even more difficult when perturbations and non-circular orbits are considered. 2.2. Satellite-to-satellite visibility prediction We assume that a satellite has access to other satellites once there is no Earth obstruction. A complicated prediction condition

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

3

where h(t ) is a multiple hump function without fixed periods, and h˙ (t ) represents the first derivative of h(t ). The interval rise and set times are easily distinguished in accordance with h˙ (t ) > 0 and h˙ (t ) < 0. Besides, one visibility interval would regress to a single point when h˙ (t ) = 0 at the root of h(t ), which is rarely seen in practice. Therefore, we intend to design a onboard metamodeling-based framework to efficiently solve these visibility prediction problems, which have been classified as a multiple hump function rootsfinding problem. 3. Method In this section, an onboard metamodeling-based framework for the roots-finding of a multiple hump function is explicitly described. Initially, Latin hypercube sampling (LHS) is adopted as preprocessing to generate the initial sampling points. For the satelliteto-site visibility prediction, we further employ the interval shrinking strategy during preprocessing. Then we introduce different metamodeling techniques to approximate the visibility prediction function. Finally, we conclude the proposed framework.

Fig. 2. Illustration of satellite-to-site view geometry.

3.1. Preprocessing

Fig. 3. Illustration of satellite-to-satellite view geometry.

with two inequalities is proposed in [23]. As illustrated in Fig. 3, r 1 and r 2 are satellite position vectors, and the satellites’ relative position vector is denoted as r 12 . The critical condition of visibility prediction is described in the figure, where r 12 is tangent to the Earth. Therefore, the satellite-to-satellite visibility prediction function is concisely expressed as 1

1

g (t ) = (r 1 2 − R c2 ) 2 + (r 2 2 − R c2 ) 2 − r 12 

(2)

where R c is the distance from the Earth center to r 12 and equal to R e

k2 +1−e 2E k2 +1

[28], in which k is the slope of r 12 . When the Earth

is considered as round, R c = R e . The satellites have access to each other if g (t ) > 0 and they do not otherwise. When the critical condition is satisfied, g (t ) = 0. 1

1

g (t ) = (r 1 2 − R e2 ) 2 + (r 2 2 − R e2 ) 2 − r 12 

(3)

where the satellites have access to each other if g (t ) > 0 and they do not otherwise. When the critical condition is satisfied, g (t ) = 0. 2.3. Roots-finding of multiple hump function As illustrated in Fig. 1, the problems of satellite-to-site and satellite-to-satellite visibility predictions can be generally expressed as

find t subject to h(t ) = 0 h˙ (t ) = 0

(4)

3.1.1. Latin hypercube sampling LHS, originally from statistical science, is a method used to produce input values for the estimation of function output variables [29,30]. In order to balance the efficiency and precision of the function approximation, one needs to sample limited points and consider the distribution of points that also affects the approximation accuracy. Given that LHS ensures uniform random sampling for each design variable and provides better sampling distribution than random sampling [31], we apply LHS to generate the initial samples for each searching interval. The number of initial sampling points is selected as three, according to the following equation [32].

min{

(n v + 1)(n v + 2) 2

, 5n v }

(5)

where n v is the dimension of the input variable, and subsequently n v = 1 with respect to f (t ) and g (t ). 3.1.2. Interval shrinking strategy Interval shrinking with an in-view cone [18,21,22] has been widely considered during preprocessing for the satellite-to-site visibility prediction. As seen in Fig. 1(a), the value of f (t ) is generally less than the threshold value λ. This is because the in-view cone generated by the ground site is fixed on the Earth’s surface and slowly rotates with the Earth in inertial space, while the satellite runs in the fixed orbital plane. Therefore, a satellite only has access to the ground site when it moves into the intersecting arc between the in-view cone and the orbit trajectory [22]. The non-intersecting periods are filtered through the interval shrinking strategy, reducing the search for intervals for visibility predictions. Typically, each interval is further divided into several smaller subintervals according to the orbital periods in practice. Interested readers are referred to [22] for more details. 3.2. Metamodeling techniques Given the high complexity of the visibility prediction function and the limited computational capacity of a satellite’s onboard system, it is more practical to approximate the visibility function using interpolation methods. Metamodeling is found to be a valuable tool to support efficient and effective function approximation [25]. Three types of metamodeling techniques are described as follows.

4

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

3.2.1. Radial basis functions RBFs are typically used to build up function approximations [33, 34]. Given a set of different measured points xi (i = 1, ..., n) and the corresponding measured values y i (i = 1, ..., n), RBF is expressed as

f r (x) =

n 

βir φ r (x − xi  ) = β T ϕ

(6)

i =1

with

T  β = β1r , β2r , · · · , βnr T  ϕ = φ r (x − x1 ) , φ r (x − x2 ) , · · · , φ r (x − xn )

(7) (8)

βir

r

where f (x) is the prediction value at x and is the corresponding coefficients for xi .  ·  denotes the Euclidean norm, and φ r (·) represents the radial basis function at xi . The coefficients βir (i = 1, ..., n) are determined as

β = A −1 y with

(9)

⎛

φ r (x1 − x1 ) · · · ⎜ .. A=⎝ . φ r (xn − x1 ) · · ·

⎞ φ r (x1 − xn ) ⎟ .. ⎠ . φ r (xn − xn )

y = ( y i , · · · , yn ) T

(10)

(11)

Following the above definitions, it has been verified that the RBF approximation passes through all measured points and reveals a high prediction accuracy around those points [35]. As the number of measured points increases, the prediction accuracy of the RBF approximation is further improved. More importantly, its implementation is ready enough for RBF to be utilized in different engineering applications [27,36]. There is a large class of radial basis functions including multiquadrics, inverse multiquadrics, Gaussian functions, thin-plate spline, and biharmonic spline. 3.2.2. Kriging Originally from geographic science, Kriging [37,38] is an advanced geostatistical procedure that generates an estimated surface function from a scattered set of points with certain values. Kriging weights the surrounding measured values to derive a prediction for an unmeasured location, where the weights are based not only on the distance between the measured points and the prediction location, but also on the overall spatial arrangement of the measured points. From this point of view, Kriging is a promising procedure for approximating the visibility function. Here Kriging is utilized to weight the surrounding measured values to derive a prediction for an unmeasured point. The general formulation is formed as

f k (x) =

n 

βik y i

(12)

i =1

where f k (x) is the prediction value at point x and βik is an unknown weight coefficient for xi . To appropriately arrange the coefficients, the point difference autocorrelation should be quantified. Then the experimental semivariogram, representing the difference relationships among the measured points, is expressed as

nd nd S V (d) =

i =1

j = i +1 ( y i C n2d

− y j )2

(13)

where d is the difference between pairs of points, nd is the number of pairs with the distance d, and S V (d) represents the value of the semivariance at time difference d. Following the experimental semivariogram, a semivariogram mathematical model is proposed. While the spherical, Gaussian, and exponential models are widely used in statistical geoscience, the selected model influences the prediction of the unknown points, particularly when the shape of the curve near the origin differs significantly. The weight coefficient for xi , whose difference from x is out of the range c e , is 0 according to the nonnegative property of the semivariogram. Then the error between the real and estimated n k semivariograms is solved with minimum variances and i =1 βi = 1. Eventually the values of βik are obtained and the visibility func-

tion prediction f k (x) is described in Equation (12).

3.2.3. Support vector regression The support vector method is a universal tool for solving function approximation problems. It was initially designed for pattern recognition [39], leading to a new technique where the function approximation is a linear expansion whose elements are nonlinear functions. When this method is employed for function approximation, the approaches are often referred to as SVR [40,41]. This novel combination opens new opportunities for solving various problems of function approximation. Given a set of points xi and corresponding measured function values y i (i = 1, ..., n), we map each point xi into a feature space where x → g (x). We then design the approximation function using SVR as

f s (x) = ω T g (x) + b

(14)

where ω is a coefficient vector and b is a threshold constant. In order to estimate the quality of f s (x), we further define a loss function as follows.

1 n

R (ω) =

n

| y i − f s (xi )|ε

(15)

i =1

with

| y i − f s (xi )|ε =

if | y i − f s (xi )| < ε ,

0

| y i − f (xi )| − ε otherwise, s

(16)

where ε is the approximation precision. By introducing slack variables and Lagrange multipliers [42], the problem of minimizing loss function R (ω ) can be solved to optimality when the approximation function f s (x) is determined as

f s (x) =

n  (αˆ i − αi )κ (x, xi ) + b

(17)

i =1

ˆ i and αi are Lagrange multipliers, κ (x, xi ) is the kernel where α function equivalent to the dot product in the feature space, and b is robustly expressed as b=

1 2

(min{ y i −

+ max{ y i −

n  (αˆ i − αi )κ (x, xi )} i =1 n  (αˆ i − αi )κ (x, xi )})

(18)

i =1

Similar to RBF, a series of kernel functions has been developed for SVR, e.g. the polynomial kernel, the hyperbolic tangent kernel, and the anova spline kernel [43].

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

5

Algorithm 1 Metamodeling based framework for rapid visibility prediction.

Algorithm 2 Self-adaptive interpolation method.

Input: Satellite and site information, mission horizon [t s , t e ], and error tolerance; Output: Visible timetables R; 1: Obtain subintervals [t si , t ei ](i = 1, 2, ..., n1 ) using the interval shrinking strategy (only for the satellite-to-site visibility prediction); 2: for each i = 1 to n1 do ij ij 3: Further divide [t si , t ei ] into subintervals [t s , t e ]( j = 1, 2, ..., n2 ) according to orbital period, and initiate sampling points for each subinterval; 4: for each j = 1 to n2 do 5: Generate new sampling points with the self-adaptive interpolation method and add these points into sampling set P ; 6: Obtain approximate visibility function with metamodeling techniques, predict rise and set times and add these points into P ; 7: if Prediction function value at time t p satisfies error tolerance then 8: Add t p into R and continue; 9: else 10: Go to line 6. 11: end if 12: end for 13: end for

Satellite and site information, subinterval [t s , t e ], initial interpolation step D k ; Output: Generated sampling points;

Input: ij

ij

ij

1: Set k = 1, t 0 = t s ;

ij

2: Generate sampling point tk = min(tk−1 + D k , t e ) and obtain value of h(tk ); 3: if h(tk ) · h(tk−1 ) < 0 then t +t 4: Generate sampling point k−12 k and obtain corresponding value of visibility function; 5: end if ij 6: if tk < t e then 7: Update interpolation step, set k = k + 1 and go to line 2; 8: else 9: Method terminates; 10: end if 11: Output the entire generated sampling points.

Table 1 Orbital parameters of satellites. Name

a (km)

e

(◦ )

i (◦ )

ω(◦ )

M (◦ )

Gaojing-1 Beidou-3M

6,903.67 27,906.60

0.00165 0.00078

50.508 166.091

97.583 55.011

97.8446 286.749

2.0288 36.965

3.3. Framework In this section, the proposed algorithm framework is described in detail. As seen in Algorithm 1, the initial mission horizon [t s , t e ] is separated into several subintervals during the preprocessing, where the non-intersecting intervals are removed for the satelliteto-site visibility prediction. The intervals are then further divided into smaller intervals in the orbital period. We select the smaller orbital period of two satellites to generate subintervals for the satellite-to-satellite visibility prediction. For each divided subinterval, five sampling points are initially selected through LHS. The self-adaptive interpolation method, which will be expatiated later, is then applied to generate more sampling points. The sampling set P now consists of points from LHS and the self-adaptive interpolation method, respectively. By utilizing sampling set P which is dynamically updated, the metamodeling techniques are iteratively invoked to approximate visibility function and predict rise and set times within the divided subinterval. For each iteration, current prediction rise and set times with real value of visibility function would be added to P as new sampling points, since larger number of sampling points ensure higher accuracy for possible further visibility prediction. The metamodeling techniques continue to work until the accuracy of the entire prediction rise and set times is satisfied. Notice that the connection between the proposed framework and metamodeling techniques is established by the procedure described in line 6 of Algorithm 1. In the next section of experimental study, different metamodeling techniques would be tested through the proposed framework. The accuracy of each prediction moment t p is checked as follows. If the absolute value of the prediction function is within the acceptable error tolerance , then t p is determined as one of the rise and set times; otherwise, we add this point into the set of sampling points, restarting the function approximation and visibility prediction process guided by the metamodeling techniques. The proposed framework continues generating rise and set times until all subintervals are considered. The self-adaptive interpolation method, shown in Algorithm 2, is designed to efficiently generate sampling points for the whole framework. Here h(t ) represents the visibility prediction function. The interpolation step D k is updated iteratively as

⎧ ⎨ max(d1 D k−1 , d0 ) D k = min(d2 D k−1 , t ei j − t si j ) ⎩ D k −1

r < r1 r > r2 r1 ≤ r ≤ r2

(19)

where d0 , d1 , d2 , r1 , and r2 are controllable parameters, and r = hˆ (tk ) | h(tkh)− |, where hˆ (t ) is the prediction approximation function. (t ) k

According to the above settings, d0 represents the minimal interpolation step, d1 stands for the step increase ratio (typically > 1) and d2 is the reduction ratio (typically < 1). A smaller r1 corresponds with a lower possibility of increasing the interpolation step, while a higher r2 correlates with a smaller possibility of shrinking the step. 4. Simulations 4.1. Data generation and experiment setup In order to comprehensively test the performance of the proposed framework, two satellites, Gaojing-1 [44] and Beidou3M [45], representing different orbit altitudes, are adopted for data generation. Gaojing-1 [44], also known as SuperView-1, is one of China’s high-resolution LEO satellites; Beidou-3M, in medium Earth orbit (MEO), belongs to the global and regional navigation system proposed by China. Details of the satellites are provided in Table 1. As seen in Fig. 4, we randomly generate 100 ground sites in the range of 60◦ N–60◦ S. The scenario’s initial time is fixed at 1 Jan 2018 00:00:00.000 UTCG with a 24-hour (86,400 seconds) horizon. The minimum elevation angle α and error tolerance for obtaining rise and set times are respectively set as 15◦ and 10−4 . The simplified general perturbation model (SGP4) [46] is consulted for generating satellite to ground site velocities and positions. The default parameters for the self-adaptive interpolation methods are set as d0 = 10, d1 = 2, d2 = 0.5, r1 = 0.25 and r2 = 0.75. Various metamodeling techniques with different surrogate models are selected to conduct the experiments. The multiquadric, inverse multiquadric, and thin plate spline models are implemented in the RBF. Similarly, cubic, spherical, and spline surrogate models are adopted for the Kriging procedure. As for the SVR technique, we utilize anova spline as its kernel function. We note the above metamodels as R-MQ, R-IMQ, R-TPS, K-CUB, K-SPH, K-SPL, and S-SPL respectively. For each combination of parameter settings, 10 runs are conducted. The computational experiments are carried out on a Windows 10 64-bit OS with Intel Core i5-7200 CPU at 2.5 GHz and 8 GB of RAM. The programs are implemented in Matlab. Since Han et al. [22] reported that the Hermite interpolation method (HIM) outperforms previous algorithms, we only

6

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

Fig. 4. Illustration of sites distribution. Table 2 Metamodeling-based framework and HIM results comparison for Gaojing-1. Method

Cost

Error(s)

PNE (%)

HIM R-MQ R-IMQ R-TPS K-CUB K-SPH K-SPL S-SPL

184.7 95.0 327.5 84.2 82.1 96.0 82.0 84.1

8.48 1.46 1.31 1.45 1.44 1.47 1.47 1.43

3.33 0.53 0.46 0.53 0.52 0.51 0.55 0.49

where t p and t a denotes the prediction and actual rise and set time, and dt a represents the corresponding actual in-view duration. 4.2. Results for satellite-to-site visibility prediction 4.2.1. Results for interval shrinking strategy Fig. 5 illustrates an example of the interval shrinking strategy employed for the different satellites. Using the interval shrinking strategy as preprocessing, the original horizon of 24 hours is cut into subintervals, which are within green lines in the figure. The non-intersecting time intervals between the site viewing cone and orbit trajectory are successfully removed. This strategy appears to be more successful in reducing wasted period for the LEO satellite; the entire subinterval to each site for Gaojing-1, after shrinking, is 4.26 hours on average, while only 3.32 hours could be removed for Beidou-M3, an MEO satellite.

Fig. 5. Illustration of interval shrinking strategy for different satellites. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

compare the results of these various metamodels against the performance of HIM. In the tables below, the column Method lists seven metamodeling methods as well as HIM. The column Cost shows the average number of samples for each run. Columns Error and PNE represent the mean prediction absolute time error and the percentage of normalized error [18] respectively. For each rise and set time, PNE is calculated as

PNE =

|t p − t a | dt a

(20)

4.2.2. Results for Gaojing-1 The results of the rapid visibility prediction for Gaojing-1 from various metamodeling techniques and HIM are reported in Table 2. Overall, metamodeling-based techniques have superior performance for Gaojing-1 as compared with HIM. The average number of sampling points for HIM is 184.7, while the corresponding numbers for metamodeling-based techniques never exceed 100 (except R-IMQ), indicating that the computational cost is decreased by 50%. As for prediction accuracy, the average absolute error of HIM for a single rise and set time is 8.48 seconds, and its mean PNE reaches 3.33%. The metamodeling-based interpolation methods (with the exception of R-IMQ) have similar performances with around 1.5 seconds of time error and about 0.50% PNE. Notably, R-IMQ represents the best accuracy for visibility prediction, although it has the highest computational cost. Also, using five second trajectory checking, the metamodeling-based framework achieves a decrease in computational cost of about 99.5% compared with the brute force method. K-SPL is recommended to tackle satellite-to-site visibility prediction problem considering LEO satellites.

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

Fig. 6. Comparisons between HIM and K-SPL for Gaojing-1.

7

Fig. 7. Comparisons between HIM and K-SPL for Beidou-M3.

Table 3 Metamodeling-based framework and HIM results comparison for Beidou-M3. Method

Cost

Error(s)

PNE (%)

HIM R-MQ R-IMQ R-TPS K-CUB K-SPH K-SPL S-SPL

49.8 89.1 — 80.3 78.8 89.2 79.0 82.4

38.22 17.55 — 17.93 17.08 17.87 17.77 17.60

0.37 0.20 — 0.19 0.19 0.20 0.19 0.19

In Fig. 6, we provide an example of visibility prediction results for Gaojing-1 in an subinterval obtained by the interval shrinking strategy. The true curve of f (t ) and the interpolation results of K-SPL and HIM are plotted in blue, red, and green lines respectively. The sampling points for K-SPL and HIM are marked in squares and circles. In general, the approximation function of HIM is closer to the true curve of f (t ), while the interpolation results of K-SPL cannot effectively track the true curve when f (t ) arrives at peak range. However, it is worth noting that we attempt to find the roots of f (t ) rather than track its whole curve. When we zoom into the first visibility rise time, the K-SPL samples point very closely to the root point of f (t ) with less than 1 second of error, while HIM predicts the root point with about 8 seconds of error. It provides an intuitive explanation of why K-SPL shows better prediction accuracy. Additionally, fewer points are sampled for K-SPL to track the peak range of f (t ), reducing the computational cost compared with HIM. 4.2.3. Results for Beidou-M3 We describe the visibility prediction comparisons between metamodeling-based techniques and HIM for Beidou-M3 in Table 3. Notice that R-IMQ does not work for the visibility prediction of Beidou-M3 as it is unable to terminate iterations within an allowable error tolerance. Although a prediction error of 38.22 seconds is observed, HIM requires the least computational cost with no more than 50 sampling points. Following HIM, K-CUB ranks second in algorithm efficiency and provides the best solution accuracy with 17.08 seconds of error and 0.19% PNE. According to the table, other metamodeling-based methods (except R-IMQ) have similar performances for Beidou-M3. Both the average prediction error and PNE decreases by 50% compared with HIM. In addition, the computational cost decreases about 99.5% for the metamodeling-based frameworks over the brute force method using five second trajectory checking. Therefore it is better select K-CUB to deal with satellite-to-site visibility prediction with MEO satellites. Similar to Fig. 6, the comparisons between HIM and K-SPL for Beidou-M3 are illustrated in Fig. 7. The sampling points for K-SPL and HIM are marked in red squares and green circles, re-

Table 4 Metamodeling-based framework and HIM results comparison for satellite-tosatellite visibility. Method

Cost

Error(s)

PNE (%)

HIM R-MQ R-IMQ R-TPS K-CUB K-SPH K-SPL S-SPL

905.0 449.3 — 391.7 477.4 474.4 377.4 410.0

0.05 0.03 — 0.05 0.02 0.03 0.02 0.03

0.008 0.005 — 0.008 0.003 0.005 0.003 0.005

spectively. The approximation function of K-SPL is very close to f (t ), while errors exist in HIM’s results, especially in the bottom range of f (t ). This is because K-SPL generates more sampling points to precisely approximate f (t ). More information is obtained when zooming into the area around the rise time; as seen in the subfigure, K-SPL samples the point very close to the actual visibility rise time, ensuring high-precision results. HIM does not sample points near the rise time, resulting in a larger prediction error. 4.3. Results for satellite-to-satellite visibility prediction The results of the rapid visibility prediction for satellites Gaojing-1 and Beidou-3M using different metamodeling techniques are reported in Table 4. Notice that the function value of g (t ) has been normalized to improve prediction accuracy. R-IMQ is unable to perform the calculation, as it has difficulty terminating iterations within an acceptable error tolerance. The remaining metamodeling-based methods have superior performance in the computational cost compared with HIM. K-SPL ranks first for algorithm efficiency, followed by R-TPS, S-SPL, R-MQ, KSPH and K-CUB. HIM applies more than 900 samples to find the roots of g (t ), more than twice as many samples as most of the metamodeling-based methods. Moreover, all of the methods provide high precision results according to columns Error and PNE. We further compare the proposed framework with the brute force method using five second trajectory checking, and decreases of 97.7% are observed in computational cost. Overall, K-SPL has been demonstrated to be the best metamodeling based method for the satellite-to-satellite visibility prediction. In Fig. 8, we illustrate examples of satellite-to-satellite visibility prediction results for K-SPL and HIM over a period of 24 hours. The true curve of g (t ) and the approximation results from K-SPL and HIM are plotted in blue, red, and green lines respectively. The sampling points for K-SPL and HIM are marked in squares and circles. Generally, both of the approximation functions are close to the true curve of g (t ), though HIM utilizes a larger number of sampling points. We zoom in on the third peak range of g (t ) in the figure,

8

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

Table 5 Results for adjusting parameters of self-adaptive interpolation method.

Gaojing-1 Beidou-3M

Cost

Cvar

Error(s)

Evar

PNE (%)

Pvar

83.0 79.1

1.33 0.25

1.35 17.53

0.06 0.23

0.49 0.19

< 0.01 < 0.01

Table 6 Results comparisons using different tolerance errors.



Gaojing-1

Beidou-3M

Cost

Error(s)

PNE (%)

Cost

Error(s)

PNE (%)

10−5 5 × 10−5 10−4 5 × 10−4 10−3 5 × 10−3

84.2 83.3 82.0 80.6 80.2 78.2

1.35 1.35 1.47 1.49 1.99 5.87

0.49 0.50 0.55 0.55 0.72 1.92

79.4 79.3 79.0 78.8 78.8 78.3

17.42 17.65 17.77 18.68 18.19 18.98

0.19 0.18 0.19 0.19 0.18 0.18

Overall

81.4

2.25

0.79

79.0

18.12

0.19

parameters are in default and K-SPL is enabled. As seen in Table 6, the computational cost gradually decreases, and the values of prediction error and PNE rise, correlating with a tolerance error increase from 10−5 to 5 × 10−3 . For satellite Gaojing-1, the prediction accuracy is acceptable until the tolerance error reaches 5 × 10−3 . For Beidou-3M, however, the results change minimally and retain high-quality precision. This is because the visibility prediction function f (t ) of Beidou-M3 has a gentle slope compared with that of Gaojing-1 (see Figs. 6 and 7), making it easier to maintain prediction accuracy for Beidou-M3. In conclusion, the tolerance setting is controllable and does not heavily influence the prediction results. Fig. 8. Illustrations of K-SPL and HIM results for satellite-to-satellite visibility.

and it shows that differences occur between the results of K-SPL and g (t ) using a few samples. Meanwhile, HIM applies a large number of points to interpolate g (t ) and achieves precise track in the peak range. This cloud explain that why HIM costs more computational time compared with K-SPL. However, when we further zoom in around one root of g (t ) in Fig. 8(a), we find that K-SPL closely samples points around actual rise and set times. Therefore, K-SPL realizes the same prediction accuracy as HIM, with fewer sampling points. 4.4. Sensitivity analysis The parameters of self-adaptive interpolation methods, especially for r1 and r2 , may impact the prediction results calculated by the metamodeling-based framework. According to Equation 19, a smaller r1 corresponds with a lower possibility of increasing the interpolation step, while a larger r2 correlates with a smaller possibility of shrinking the step. To test the influence of these parameters, we adjust r1 from 0.1 to 0.3 and r2 from 0.7 to 0.9 at increments of 0.05. Then 25 combinations of r1 and r2 are applied to the satellite-to-site visibility prediction for each site and satellite in which K-SPL is selected for function approximation. The results are reported in Table 5. Columns Cvar, Evar and Pvar respectively represent the variance for the number of sampling points, absolute prediction time error, and PNE with different parameter settings. For both satellites Gaojing-1 and Beidou-3M, a change in the parameters does not have a significant effect on the prediction results, which implies that the self-adaptive interpolation method is parameter-insensitive. We conduct additional experiments for satellite-to-site visibility prediction, using different tolerance errors in which the remaining

5. Conclusion Onboard satellite autonomous systems have significant potential for space missions. In order to take advantage of these systems, rapid onboard satellite visibility predictions are necessary for certain decision-making missions. We have developed mathematical models for satellite-to-site and satellite-to-satellite visibility predictions. Then a metamodeling-based framework is proposed for solving these models, which can be regarded as roots-finding process of multiple hump functions. We conduct a series of computational experiments based upon actual data from two China’s satellites and, compared with traditional brute force and existing interpolation algorithms, the proposed framework has an overall superior performance. Moreover, the metamodeling techniques are easily used with more complicated satellite onboard visibility prediction problems, such as considering limitations with the solar zenith angle and field-of-camera view. Declaration of competing interest None declared. Acknowledgement This research was supported by the Academic Excellence Foundation of BUAA for PhD students. References [1] M.D. Rayman, P. Varghese, D.H. Lehman, L.L. Livesay, Results from the deep space 1 technology validation mission, Acta Astronaut. 47 (2–9) (2000) 475–487. [2] S. Chien, R. Sherwood, G. Rabideau, R. Castano, A. Davies, M. Burl, R. Knight, T. Stough, J. Roden, P. Zetocha, et al., The Techsat-21 autonomous space science agent, in: Proceedings of the First International Joint Conference on Autonomous Agents and Multiagent Systems, 2002, pp. 570–577.

X. Wang et al. / Aerospace Science and Technology 94 (2019) 105377

[3] D. Tran, S. Chien, R. Sherwood, R. Castano, B. Cichy, A. Davies, G. Rabideau, The autonomous sciencecraft experiment onboard the eo-1 spacecraft, in: Proceedings of the Fourth International Joint Conference on Autonomous Agents and Multiagent Systems, 2005, pp. 163–164. [4] M.J. Barnsley, J.J. Settle, M.A. Cutter, D.R. Lobb, F. Teston, The PROBA/CHRIS mission: a low-cost smallsat for hyperspectral multiangle observations of the Earth surface and atmosphere, IEEE Trans. Geosci. Remote Sens. 42 (7) (2004) 1512–1520. [5] J. Barreiro, M. Boyce, M. Do, J. Frank, M. Iatauro, T. Kichkaylo, P. Morris, J. Ong, E. Remolina, T. Smith, et al., EUROPA: a platform for AI planning, scheduling, constraint programming, and optimization, in: 4th International Competition on Knowledge Engineering for Planning and Scheduling, ICKEPS, 2012, pp. 1–14. ¸ I.A. Nesnas, M. Pavone, Spacecraft Autonomy Chal[6] J.A. Starek, B. Açıkmese, lenges for Next-Generation Space Missions, Springer, 2016. [7] S. Nag, J. LeMoigne, O. de Weck, Cost and risk analysis of small satellite constellations for Earth observation, in: 2014 IEEE Aerospace Conference, 2014, pp. 1–16. [8] H. Kim, Y.K. Chang, Mission scheduling optimization of sar satellite constellation for minimizing system response time, Aerosp. Sci. Technol. 40 (2015) 17–32. [9] X. Wang, Z. Chen, C. Han, Scheduling for single agile satellite, redundant targets problem using complex networks theory, Chaos Solitons Fractals 83 (2016) 125–132. [10] X. Wang, R. Leus, C. Han, Fixed interval scheduling of multiple Earth observation satellites with multiple observations, in: 2018 9th International Conference on Mechanical and Aerospace Engineering, ICMAE, 2018, pp. 28–33. [11] Y. She, S. Li, Y. Zhao, Onboard mission planning for agile satellite using modified mixed-integer linear programming, Aerosp. Sci. Technol. 72 (2018) 204–216. [12] S. Jin, et al., Recent progresses on Beidou/COMPASS and other global navigation satellite systems (GNSS) – I, Adv. Space Res. 51 (6) (2013) 941. [13] D. Yang, J. Yang, P. Xu, Timeslot scheduling of inter-satellite links based on a system of a narrow beam with time division, GPS Solut. 21 (3) (2017) 999–1011. [14] Z. Zheng, J. Guo, E. Gill, Onboard mission allocation for multi-satellite system in limited communication environment, Aerosp. Sci. Technol. 79 (2018) 174–186. [15] J. Lawton, Numerical method for rapidly determining satellite-satellite and satellite-ground station in-view periods, J. Guid. Control Dyn. 10 (1) (1987) 32–36. [16] P.R. Escobal, Rise and set time of a satellite about an oblate planet, AIAA J. 1 (10) (1963) 2306–2310. [17] Y. Mai, P. Palmer, Fast algorithm for prediction of satellite imaging and communication opportunities, J. Guid. Control Dyn. 24 (6) (2001) 1118–1124. [18] I. Ali, N. Al-Dhahir, J.E. Hershey, Predicting the visibility of LEO satellites, IEEE Trans. Aerosp. Electron. Syst. 35 (4) (1999) 1183–1190. [19] P. Palmer, Y. Mai, A fast prediction algorithm of satellite passes, in: 14th Annual AIAA/USU Conference on Small Satellite, 2000. [20] S. Alfano, D. Negron Jr, J.L. Moore, Rapid determination of satellite visibility periods, J. Astronaut. Sci. 40 (1992) 281–296. [21] H. Cui, C. Han, Satellite constellation configuration design with rapid performance calculation and ordinal optimization, Chin. J. Aeronaut. 24 (5) (2011) 631–639. [22] C. Han, X. Gao, X. Sun, Rapid satellite-to-site visibility determination based on self-adaptive interpolation technique, Sci. China, Technol. Sci. 60 (2) (2017) 264–270.

9

[23] X. Sun, H. Cui, C. Han, G. Tang, APCHI technique for rapidly and accurately predicting multi-restriction satellite visibility, in: Proceedings of the 22nd AAS/AIAA Space Flight Mechanics Meeting, Charleston, South Carolina, 2012. [24] M.F. Hussain, R.R. Barton, S.B. Joshi, Metamodeling: radial basis functions, versus polynomials, Eur. J. Oper. Res. 138 (1) (2002) 142–154. [25] G.G. Wang, S. Shan, Review of metamodeling techniques in support of engineering design optimization, J. Mech. Des. 129 (4) (2007) 370–380. [26] I. Rosenbaum, J. Staum, Multilevel Monte Carlo metamodeling, Oper. Res. 65 (4) (2017) 1062–1077. [27] J. Niu, J. Lei, J. He, Radial basis function mesh deformation based on dynamic control points, Aerosp. Sci. Technol. 64 (2017) 122–132. [28] E.H. Lockwood, A Book of Curves, Cambridge University Press, 1967. [29] M.D. McKay, R.J. Beckman, W.J. Conover, Comparison of three methods for selecting values of input variables in the analysis of output from a computer code, Technometrics 21 (2) (1979) 239–245. [30] M.D. Shields, J. Zhang, The generalization of Latin hypercube sampling, Reliab. Eng. Syst. Saf. 148 (2016) 96–108. [31] D. Wang, Collaboration pursuing method using Latin hypercube sampling, in: 47th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, 2006, p. 1713. [32] T. Long, L. Wang, D. Wu, X. Guo, L. Liu, A simultaneous computing framework for metamodel-based design optimization, in: ASME 2014 International Design Engineering Technical Conferences, 2014, p. V02BT03A022. [33] J. Park, I.W. Sandberg, Approximation and radial-basis-function networks, Neural Comput. 5 (2) (1993) 305–316. [34] G.-B. Huang, P. Saratchandran, N. Sundararajan, A generalized growing and pruning RBF (GGAP-RBF) neural network for function approximation, IEEE Trans. Neural Netw. 16 (1) (2005) 57–67. [35] D. Wang, G.G. Wang, G.F. Naterer, Collaboration pursuing method for multidisciplinary design optimization problems, AIAA J. 45 (5) (2007) 1091–1103. [36] M.E. Biancolini, Fast Radial Basis Functions for Engineering Applications, Springer, 2017. [37] M.L. Stein, Interpolation of Spatial Data: Some Theory for Kriging, Springer Science & Business Media, 2012. [38] I. Couckuyt, T. Dhaene, P. Demeester, ooDACE toolbox: a flexible object-oriented Kriging implementation, J. Mach. Learn. Res. 15 (1) (2014) 3183–3186. [39] G. Baudat, F. Anouar, Kernel-based methods and function approximation, in: International Joint Conference on Neural Networks Proceedings, Vol. 2, 2001, pp. 1244–1249. [40] A.J. Smola, B. Schölkopf, A tutorial on support vector regression, Stat. Comput. 14 (3) (2004) 199–222. [41] R. Khelif, B. Chebel-Morello, S. Malinowski, E. Laajili, F. Fnaiech, N. Zerhouni, Direct remaining useful life estimation based on support vector regression, IEEE Trans. Ind. Electron. 64 (3) (2017) 2276–2285. [42] S.M. Clarke, J.H. Griebsch, T.W. Simpson, Analysis of support vector regression for approximation of complex engineering analyses, J. Mech. Des. 127 (6) (2005) 1077–1087. [43] C.-C. Chuang, S.-F. Su, J.-T. Jeng, C.-C. Hsiao, Robust support vector regression networks for function approximation with outliers, IEEE Trans. Neural Netw. 13 (6) (2002) 1322–1330. [44] SpaceView, Superview-1, available at http://www.spaceview.com/SuperView1English/index.html, 2018. [45] Beidou, Beidou navigation satellite system, available at http://en.beidou.gov.cn, 2018. [46] N.Z. Miura, Comparison and design of simplified general perturbation models, available at http://hpiers.obspm.fr/eop-pc, 2018.