Space-Based Passive Orbital Maneuver Detection Algorithm for High-Altitude Situational Awareness

Abstract

Orbital maneuver detection for non-cooperative targets in space is a key task in space situational awareness. This study develops a passive maneuver detection algorithm using line-of-sight angles measured by a space-based optical sensor, especially for targets in high-altitude orbit. Emphasis is placed on constructing a new characterization for maneuvers as well as the corresponding detection method. First, the concept of relative angular momentum is introduced to characterize the orbital maneuver of the target quantitatively, and the sensitivity of the proposed characterization is analyzed mathematically. Second, a maneuver detection algorithm based on the new characterization is designed in which sliding windows and correlations are utilized to determine the mutation of the maneuver characterization. Subsequently, a numerical simulation system composed of error models, reference missions and trajectories, and computation models for estimating errors is established. Then, the proposed algorithm is verified through numerical simulations for both long-range and close-range targets. The results indicate that the proposed algorithm is effective. Additionally, the sensitivity of the proposed algorithm to the width of the sliding window, accuracy of the optical sensor, magnitude and number of maneuvers, and different relative orbit types is analyzed, and the sensitivity of the new characterization is verified using simulations.

1. Introduction

With the rapid development of space technology, the frequency of launching satellites in various countries has increased rapidly. For example, SpaceX plans to launch 42,000 satellites for Starlink. The huge number of launches may have disastrous consequences for near-Earth space [1,2,3,4,5]. The falling off of objects produced in the launch process, and failed satellites at the end of their service life or after failure may affect the normal operation of active satellites or even cause collisions. Collisions and secondary collisions will produce countless space debris, which is very difficult to remove and poses a major threat to the safety of satellites [6,7,8,9]. In addition, the mastery of space requires huge national interests. Space provides an important support for national security. Fierce competition in space utilization technology has occurred, which makes the importance of space security increasingly prominent [9,10,11]. Therefore, to protect space resources, whether to avoid damage caused by space debris or to defend against space military threats, it is extremely critical to improve space situational awareness [12].

To date, many situational awareness methods have been proposed. Based on the environment in which the measuring equipment is located, the existing situational awareness methods can be divided into space-based awareness and ground-based awareness methods. The main method of ground-based sensing establishes an observation station on the ground, and the operational state of the target is observed using ground-based radar and photoelectric systems [13,14,15]. However, these methods are not only constrained by political factors but also have the problems of high construction costs and complex maintenance. Therefore, increasing attention has been paid to research focusing on space-based situational awareness. Space-based methods can overcome many of the shortcomings of ground-based methods, such as the direct observation of space targets in space which can overcome the influence of the atmosphere and the limitation of the field of regard. At present, researchers have proposed many space-based sensing methods. Several space-based sensing systems have been used in practical applications, like the MiTEx satellite of the United States [16], the GSSAP project [16,17,18], and the ESPAStar platform [19].

Sensors for space-based sensing generally include lidar, radio ranging equipment, and optical cameras [16,17,18,19,20,21,22,23,24]. Among space-based navigation methods, the optical angles-only navigation method has become the prominent development trend because of its relatively simple equipment and high reliability [21,22,23,24]. Considering the related measurement errors in practical applications, the Kalman filter method is the most commonly introduced to deal with the errors in measurement. However, the traditional Kalman filter method can only deal with known perturbations and common measurement errors. When the target satellite performs unknown orbit maneuvers due to collisions, perturbations, or orbit changes, the orbit determination error often diverges gradually over time. Therefore, navigation based on this method will lose its effectiveness in these cases. To solve this problem, it is necessary to find a way to detect the orbit maneuver and estimate the maneuver time.

Researchers have studied orbit maneuver detection from multiple angles. Some studies have started with improving the navigation accuracy or data processing methods [25,26,27,28,29,30]. For example, Jia et al. [25] proposed an improved Kalman filter algorithm using four cooperative satellites, which could achieve the goal of detecting the target orbit maneuver accurately, but the cost increases a lot, like for more individual cooperation or the high cost of high-precision sensors. Previous studies [26,27] proposed empirical mode decomposition (EMD) processing for the detected quantity, which allowed the abnormal points of parameters to be highlighted among the observed values to realize the identification of maneuvers or abnormal points. Scientists have proposed these methods because the state parameters of general space targets, such as velocity and position, are not sensitive to orbit maneuvers, and their changes caused by orbit maneuvers are easily submerged in measurement errors, particularly for long-range targets. Therefore, identifying new feature quantities that are sensitive to maneuvers is a new trend in this field [31,32,33,34,35,36]. For example, Yu et al. and Huang et al. [31,32] proposed the use of the semi-major axis and eccentricity as detection quantities for maneuvers, but these quantities have certain limitations in the use scenarios. Liu et al. [33] proposed a maneuvering weighted fusion of multiple hypothesis tests (WFMHTs) method for space-based angles-only measurement, in which the observation error in the extended Kalman filter calculation was set as the characterization, and target maneuvers were determined by judging whether the residual had the characteristics of zero-mean white noise. The method proposed by [13] used the orbit determination residual to evaluate the start and end times of maneuvers. Wang et al. [34] proposed using the calculated rate of change in the distance as the characterization and transformed the orbit maneuver detection problem into a hypothesis test problem. Although these methods are valuable, they all have certain application limitations. Some are not suitable for small orbit maneuver detection, whereas others have poor real-time performance or some limitations.

The objective of this work is to develop a maneuver detection algorithm for non-cooperative targets in high-altitude orbits, where angles-only measurements are used. The main contribution is focused on two aspects: a new characterization based on relative angular momentum is modeled to solve the non-sensitive problem of the existing characterization in high-altitude relative motion application; the concept of the correlation coefficient is introduced to design a mutation detection algorithm for the characterized time series. Compared to existing methods, the proposed algorithm is more sensitive to maneuvering where a high precision detection is achieved.

The rest of the manuscript is organized as follows: Following the problem statement in Section 2, the model of the maneuvering characterization and its sensitivity analysis are presented in Section 3. The maneuver detection algorithm is addressed in Section 4. The numerical simulation framework is described in Section 5. The results of the simulations and analysis are presented in Section 6. Finally, the conclusions are given in Section 7.

2. Problem Statement

Based on space-based observation, the target’s line-of-sight angle is measured to provide a continuous evaluation of whether the target is performing an orbital maneuver, and the orbital maneuvering time is estimated. The logic framework of the continuous detection algorithms is given in Section 2.1, and the kinematic model used in this study is introduced in Section 2.2; finally, the measurement model of the target line-of-sight angle is introduced in Section 2.3.

2.1. Framework of the Maneuver Detection Algorithm Based on Angles-Only Measurements

The basic logical relationship of the non-cooperative target orbit maneuver angles-only detection method proposed in this study is shown in Figure 1. First, the line-of-sight angle measurement (the azimuth az and elevation el) is filtered and estimated in combination with the dynamic model to obtain the relative position vector ${\mathit{r}}_{rel}$ and relative velocity vector ${\mathit{v}}_{rel}$ of the non-cooperative target. Then, the maneuver-detecting characterization (the relative angular momentum ${\mathit{h}}_{rel}$) is constructed based on the relative orbit parameters. Finally, the maneuver time is determined through the inspection and discrimination of the characterization, and the process of determining the initial orbit is restarted at the same time.

Because the focus of this study is the design of the maneuver detection algorithm, the angle measurement and orbit determination algorithm including the Unscented Kalman filter (UKF)-based filter in Figure 1 is not discussed in detail, and the method proposed by Gong et al. [37] is adopted directly. The mapping relationship between the history of the line-of-sight angle and the relative orbit parameters is given as follows:

$${\left[{\mathit{r}}^T,{\mathit{v}}^T\right]}^T=\mathrm{Mapping}(az,el)$$

(1)

where az and el are the history of the azimuth and elevation, respectively.

2.2. Relative Motion Dynamics Model

As shown in Figure 2, a rotating frame, i.e., a local vertical–local horizontal (LVLH) reference coordinate system, is established. The origin is located at the center-of-mass (COM) of the chaser satellite, the X-axis points to the COM of the satellite from the Earth’s center, the Z-axis is normal to the orbital plane, and the Y-axis completes the right-hand system.

When the orbital eccentricity, $e$, of the two satellites satisfies $0\le e$ <1, the nonlinear relative dynamics model proposed by Yamanaka et al. [38] can be used to model the relative state between the chaser satellite and the target satellite:

$$\left[\begin{array}{c}\ddot{x}\\ \ddot{y}\\ \ddot{z}\end{array}\right]=\left[\begin{array}{c}-\mathrm{k}{\omega }^{\frac{3}{2}}x+2\omega \dot{z}+\dot{\omega }z+{\omega }^{2}x\\ -\mathrm{k}{\omega }^{\frac{3}{2}}y\\ 2\mathrm{k}{\omega }^{\frac{3}{2}}z-2\omega \dot{x}-\dot{\omega }x+{\omega }^{2}z\end{array}\right]+{\mathit{a}}_f+{\mathit{a}}_{cd}-{\mathit{a}}_{td}$$

(2)

where ${[x,y,z]}^T$ is the position of the target’s COM relative to the COM of the chaser satellite in the LVLH frame, $\omega$ is the orbital angular velocity of the chaser satellite, ${\mathit{a}}_f$ is the acceleration caused by the propulsion force of the target satellite, ${\mathit{a}}_{cd}$ is the acceleration caused by forces other than the inverse square gravity term on the chaser satellite, and ${\mathit{a}}_{td}$ is the acceleration caused by forces other than the inverse square gravity term on the target satellite. The constant $\mathrm{k}$ is defined as follows:

$$\mathrm{k}\cong \frac{\mathsf{\mu }}{h^{3/2}}$$

(3)

where $\mathsf{\mu }$ is the gravitational constant, and $h$ is the orbital angular momentum.

In this paper, we assume the chaser satellite is flying in free motion and the external forces on the chaser satellite and the target satellite are identical.

$${\mathit{a}}_f=0, {\mathit{a}}_{cd}={\mathit{a}}_{td}$$

(4)

2.3. Measurement Model

In this study, it is assumed that the orbit and attitude of the chaser satellite are known, and only passive optical cameras on the satellite can be used to measure non-cooperative targets. According to the engineering constraints, the camera is installed away from the chaser’s center-of-mass (COM), and the measurement is shown in Figure 3. The offset vector of the camera is denoted by $\mathit{d}={[d_x,d_y,d_z]}^T$ in the LVLH frame. Then, the line-of-sight vector from the camera to the target in the LVLH frame can be modeled as follows:

$$i_{LOS}=r-d=\left[\begin{array}{c}x-d_x\\ y-d_y\\ z-d_z\end{array}\right]$$

(5)

where x, y, and z are the components of the target’s relative position $r$.

Next, the line-of-sight angles, i.e., azimuth az and elevation el, can be obtained.

$$\left\{\begin{array}{l}a{z}_i={\tan }^{-1}\left(\frac{z-d_z}{y}\right)+n_1\\ e{l}_i={\tan }^{-1}\left(\frac{x-d_x}{\sqrt{{\left(y-d_y\right)}^2+{\left(z-d_z\right)}^2}}\right)+n_2\end{array}\right.$$

(6)

where ${[n_1,n_2]}^T$ is the measurement noise which is commonly modeled as zero-mean Gaussian noise.

The offset of the camera can increase the observability of relative orbit if the offset satisfies the observable criterion. Moreover, the increment of observability depends on the ratio of the effective offset to the relative range between two spacecraft. Detailed analysis for the angles-only observability from the offset can be found in a study by Gong et al. [37].

3. Modeling of the Maneuvering Characterization and Its Sensitivity Analysis

As analyzed in Section 1, the existing algorithms that use features such as semi-major axis and eccentricity to detect the orbit maneuver suffer from poor sensitivity in many applications, especially in the cases of high-altitude orbit, e.g., GEO. Therefore, a new maneuver detection feature based on the concept of relative angular momentum is proposed for high-altitude orbit cases in this study. The following sections first describe the definition of the characterization to be tested and then analyze and verify its sensitivity.

3.1. Concept of Relative Angular Momentum

The relative angular momentum is defined as follows:

$${\mathit{h}}_{rel}={\mathit{r}}_{rel}\times {\mathit{v}}_{rel}$$

(7)

where ${\mathit{r}}_{rel}={\left[x,y,z\right]}^T\in {\mathbb{R}}^3$, and ${\mathit{v}}_{rel}={\left[\dot{x},\dot{y},\dot{z}\right]}^T\in {\mathbb{R}}^3$.

A maneuver directly affects the position and velocity vectors of the satellite, and the calculation of the angular momentum is related to the position and velocity vector of the object [6]. The cross product of relative position and velocity vectors can not only reflect the changes in the values of relative position and relative velocity, but it can also directly reflect the changes in the direction of the relative position and relative velocity through the change in the value of the angular momentum, which makes the relative angular momentum more sensitive to satellite maneuvers than individual velocity or position. In a very short time after maneuver, the value of position is large compared to velocity, but its change is close to zero, and the value of velocity and of its change are both very small. Relative angular momentum amplifies the change in velocity by the cross product of velocity and position so that its change is easier to detect than that of velocity. Compared with the position almost unchanged in a very short time after maneuver, the relative angular momentum has more obvious changes because of the implicit velocity information, so it is easier to detect. In addition, using the relative angular momentum as a characterization to detect maneuvers has fewer restrictions, such as not restricting a specific orbit of the target satellite, such as some methods only applicable to targets in GEO [39], or a specific way of maneuver, such as maneuvering at perigee or apogee in the direction of velocity [32].

In addition, compared with the inertial orbital angular momentum in the geocentric coordinate system, the relative angular momentum defined in Equation (7) is more sensitive to maneuvers. The satellite carries limited energy, and thus the impulse of satellite maneuvers is small [13,40], resulting in a very small speed variation in a short time after the maneuver. However, the velocities of low-Earth-orbiting (LEO) satellites are on the order of 7.5 km/s. Therefore, the variation in the inertial position and velocity of the satellite in a short time after a maneuver is very small compared with its initial position and velocity before the maneuver. The angular momentum in a geocentric coordinate system consists of the position and velocity relative to the Earth’s center, which makes it difficult to identify the variation in the angular momentum within the noise. However, based on the assumption mentioned in the abstract that the two spacecraft are in high-altitude orbits and the relative range between them is far shorter than their geocentric distance, if the chaser satellite is taken as the reference point, the relative angular momentum in this coordinate system consists of the relative position and relative velocity, which greatly reduces the gap between the variation and the initial relative velocity. Therefore, the risk of the variation in the relative angular momentum being submerged by noise can be effectively reduced.

3.2. Sensitivity Analysis

To mathematically prove the quantitative analysis after the definition of the relative angular momentum in the previous section, the sensitivity analysis for this new characterization is conducted in this section. It should be noted that the following analysis is based on the assumption mentioned in the abstract that the two spacecraft are in high-altitude orbits and the relative range between them is far shorter than their geocentric distance.

According to Equation (7), the relative angular momentum at $t_1$ before the maneuver and $t_2$ after the maneuver can be obtained as follows:

$$\left\{\begin{array}{l}{\mathit{h}}_{rel1}={\mathit{r}}_{rel1}\times {\mathit{v}}_{rel1}\\ {\mathit{h}}_{rel2}={\mathit{r}}_{rel2}\times {\mathit{v}}_{rel2}\end{array}\right.$$

(8)

where ${\mathit{r}}_{rel1}\in {\mathbb{R}}^3$ and ${\mathit{r}}_{rel2}\in {\mathbb{R}}^3$ are the positions of the target satellite relative to the chaser satellite at $t_1$ and $t_2$, respectively, ${\mathit{v}}_{rel1}\in {\mathbb{R}}^3$ and ${\mathit{v}}_{rel2}$ $\in {\mathbb{R}}^3$ are the velocities of the target satellite relative to the chaser satellite at $t_1$ and $t_2$, respectively,$\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{h}}_{rel1}\in {\mathbb{R}}^3$ and ${\mathit{h}}_{rel2}\in {\mathbb{R}}^3$ are the relative angular momentums of the target satellite relative to the chaser satellite at $t_1$ and $t_2$, respectively.

For a single-impulse maneuver, the maneuver time is very short and the position and velocity at $t_2$ are as follows:

$$\left\{\begin{array}{l}{\mathit{r}}_{rel2}\approx {\mathit{r}}_{rel1}\\ {\mathit{v}}_{rel2}={\mathit{v}}_{rel1}+∆{\mathit{v}}_{rel}\end{array}\right.$$

(9)

where $∆{\mathit{v}}_{rel}$ is the change in velocity.

Then, the relative angular momentum at $t_2$ is

$${\mathit{h}}_{rel2}={\mathit{r}}_{rel2}\times {\mathit{v}}_{rel2}={\mathit{r}}_{rel1}\times \left({\mathit{v}}_{rel1}+\mathsf{\Delta }{\mathit{v}}_{rel}\right)={\mathit{r}}_{rel1}\times {\mathit{v}}_{rel1}+{\mathit{r}}_{rel1}\times \mathsf{\Delta }{\mathit{v}}_{rel}$$

(10)

Therefore, the change in the relative angular momentum is

$$\mathsf{\Delta }{\mathit{h}}_{rel}={\mathit{h}}_{rel2}-{\mathit{r}}_{rel1}\times {\mathit{v}}_{rel1}={\mathit{r}}_{rel1}\times \mathsf{\Delta }{\mathit{v}}_{rel}$$

(11)

Similar to the process of Equations (7)–(11), the change in angular momentum in the geocentric inertial frame can be derived from the following steps: First, define the angular momentum:

$$\left\{\begin{array}{l}{\mathit{h}}_{e1}={\mathit{r}}_{e1}\times {\mathit{v}}_{e1}\\ {\mathit{h}}_{e2}={\mathit{r}}_{e2}\times {\mathit{v}}_{e2}\end{array}\right.$$

(12)

where ${\mathit{r}}_{e1}\in {\mathbb{R}}^3$ and ${\mathit{r}}_{e2}\in {\mathbb{R}}^3$ are the positions of the target satellite in the geocentric inertial frame at $t_1$ and $t_2$, respectively, ${\mathit{v}}_{e1}\in {\mathbb{R}}^3$ and ${\mathit{v}}_{e2}$ $\in {\mathbb{R}}^3$ are the velocities of the target satellite relative to the geocentric inertial frame at $t_1$ and $t_2$, respectively,$\mathrm{a}\mathrm{n}\mathrm{d}{\mathit{h}}_{e1}\in {\mathbb{R}}^3$ and ${\mathit{h}}_{e2}\in {\mathbb{R}}^3$ are the angular momentums of the target satellite in the geocentric inertial frame at $t_1$ and $t_2$, respectively. Similar to Equation (10), the angular momentum at $t_2$ is

$${\mathit{h}}_{e2}={\mathit{r}}_{e2}\times {\mathit{v}}_{e2}={\mathit{r}}_{e1}\times \left({\mathit{v}}_{e1}+\mathsf{\Delta }{\mathit{v}}_e\right)={\mathit{r}}_{e1}\times {\mathit{v}}_{e1}+{\mathit{r}}_{e1}\times \mathsf{\Delta }{\mathit{v}}_e$$

(13)

Therefore, the change in the angular momentum in the geocentric inertial frame is

$$\mathsf{\Delta }{\mathit{h}}_e={\mathit{r}}_{e1}\times \mathsf{\Delta }{\mathit{v}}_e={\mathit{r}}_{e1}\times \left(C_{lvlh}^e·\mathsf{\Delta }{\mathit{v}}_{rel}\right)$$

(14)

where the subscript e is used to denote the ${\mathit{r}}_{e1}={[x_{e1},y_{e1},z_{e1}]}^T$, $\mathsf{\Delta }{\mathit{v}}_e={\mathit{v}}_{e2}-{\mathit{v}}_{e1}$, and ${\mathit{v}}_{rel}={[{\dot{x}}_{e1},{\dot{y}}_{e1},{\dot{z}}_{e1}]}^T$, and $C_{lvlh}^e$ is the direction cosine matrix from the LVLH frame to the geocentric inertial frame which can be calculated from the orbit position and velocity vector.

Because the change in velocity provided by a single-impulse is approximately 0.1~1 m/s [13,40], and the speed of the satellite in the geocentric coordinate system is at least close to 7.5 km/s, it is much larger than the single-impulse velocity increment:

$$‖\mathsf{\Delta }{\mathit{v}}_e‖=‖{\mathit{v}}_{e2}-{\mathit{v}}_{e1}‖\ll ‖{\mathit{v}}_{e1}‖$$

(15)

Then, substituting into Equations (12) and (14) yields the following:

$$h_{e1}\gg ∆h_e$$

(16)

where $h=norm\left(\mathit{h}\right)$.

Therefore, if the angular momentum in the geocentric inertial frame is selected as the maneuver characterization, it is difficult to distinguish the maneuver from the noisy estimation in a short time. Relative angular momentum, on the other hand, measures the relative change in the LVLH frame. For the common orbit types, the orbit shape of the two spacecraft with respect to the geocentric inertial coordinate system changes little before and after the maneuver. However, the relative orbits in the LVLH frame are very different, which also shows the superiority of relative angular momentum from the other side. Thus, the relative motion change is much larger than the absolute motion change. However, it can be found that the following inequality holds:

$$\frac{‖∆{\mathit{v}}_e‖}{‖{\mathit{v}}_{e1}‖}\ll \frac{‖∆{\mathit{v}}_{rel}‖}{‖{\mathit{v}}_{rel1}‖}$$

(17)

Then, the following inequality is established:

$$\frac{‖{\mathit{r}}_{e1}\times ∆{\mathit{v}}_e‖}{‖{\mathit{r}}_{e1}\times {\mathit{v}}_{e1}‖}=\frac{∆h_e}{h_{e1}}\ll \frac{∆h_{rel}}{h_{rel1}}=\frac{‖{\mathit{r}}_{rel1}\times ∆{\mathit{v}}_{rel}‖}{‖{\mathit{r}}_{rel1}\times {\mathit{v}}_{rel1}‖}$$

(18)

From Equation (18), it can be seen that the relative angular momentum works better as a characterization to detect maneuvers.

Next, when the relative orbit estimation error is considered, Equation (10) is switched to the following form:

$$\begin{array}{c}{\mathit{h}}_{rel2}=({\mathit{r}}_{rel2}+{\mathit{n}}_{\mathit{r}})\times ({\mathit{v}}_{rel2}+{\mathit{n}}_{\mathit{v}})=({\mathit{r}}_{rel1}+{\mathit{n}}_{\mathit{r}})\times \left({\mathit{v}}_{rel1}+\mathsf{\Delta }{\mathit{v}}_{rel}+{\mathit{n}}_{\mathit{v}}\right)\\ \hspace{1em} ={\mathit{r}}_{rel1}\times {\mathit{v}}_{rel1}+{\mathit{r}}_{rel1}\times \mathsf{\Delta }{\mathit{v}}_{rel}+{\mathit{r}}_{rel1}\times {\mathit{n}}_{\mathit{v}}+{\mathit{n}}_{\mathit{r}}\times {\mathit{v}}_{rel1}+{\mathit{n}}_{\mathit{r}}\times \mathsf{\Delta }{\mathit{v}}_{rel}+{\mathit{n}}_{\mathit{r}}\times {\mathit{n}}_{\mathit{v}}\end{array}$$

(19)

where ${\mathit{n}}_r$ and ${\mathit{n}}_v$ are the errors in the relative position and relative speed, respectively.

Due to the fact that the estimated error is much smaller than the state when the system is observable, the following inequality holds:

$$\left\{\begin{array}{l}‖{\mathit{n}}_r‖\ll ‖{\mathit{r}}_{rel1}‖\\ ‖{\mathit{n}}_v‖\ll ‖{\mathit{v}}_{rel}‖,‖\mathsf{\Delta }{\mathit{v}}_{rel}‖\end{array}\right.$$

(20)

Then, combining Equations (19) and (20) yields

$$‖{\mathit{r}}_{rel1}\times \mathsf{\Delta }{\mathit{v}}_{rel}+{\mathit{n}}_{\mathit{r}}\times \mathsf{\Delta }{\mathit{v}}_{rel}‖\gg ‖{\mathit{r}}_{rel1}\times {\mathit{n}}_{\mathit{v}}+{\mathit{n}}_{\mathit{r}}\times {\mathit{v}}_{rel1}+{\mathit{n}}_{\mathit{r}}\times {\mathit{n}}_{\mathit{v}}‖$$

(21)

Equation (21) shows that the change in the relative angular momentum caused by the orbit maneuver will be significantly different from the change caused by the relative orbit estimation error. Therefore, as long as the relative orbit estimation error satisfies the conditions in Equation (20), the relative angular momentum is quite sensitive to the target’s maneuver.

4. Maneuver Detection Algorithm

If the resultant force acting upon the target satellite does not point to the sensing satellite (the origin of the LVLH reference system), its relative angular momentum will change as follows [41]:

$$\frac{d{\mathit{h}}_{rel}}{dt}=\tau$$

(22)

where $\tau$ is the unit mass moment on the target satellite relative to the COM of the chaser satellite.

When satellites are flying in free motion and the target satellite does not execute a maneuver, $\tau$ is only provided by inertial forces, gravity, and related perturbation forces. As a result, the change in $\tau$ is very gentle, i.e., the change in ${\mathit{h}}_{rel}$ is continuous during the coasting flight. On the contrary, if there is a sudden change in ${\mathit{h}}_{rel}$, there must be a thrust or collision force put on the target. Therefore, when the orbit of the target changes, ${\mathit{h}}_{rel}$ will exhibit a sudden change, which can be naturally used to check the maneuver and determine the maneuver epoch.

However, it is difficult to use the individual time point information to determine the true single mutation point because the fluctuating local information value caused by the measurement is very likely to be identified as the true mutation point. Therefore, we aim to use more information to characterize the data and make decisions to reduce this possibility. According to Equation (15), the information segments before the mutation point maintain the characteristics of continuous change, and they have a high correlation. However, when a certain information segment contains the mutation point, its trend changes, and thus the correlations between the information segments containing the mutation point information and all the previous information segments are very low. Using this phenomenon, we can divide the information into two types of segments with different characteristics: the information segments before the mutation point and the information segments containing the mutation point. For example, Figure 4 shows a group of virtual data containing a mutation point. The data before and after the mutation point satisfy Gaussian distributions, and a, b, c, and d are the information segments of four sequential sliding windows with equal widths. In Figure 4, we can see that although some local trends (the trends of any number of sequential points, but the number is far less than the total) of a, b, and c are different owing to noise, the overall trends are almost the same, i.e., a, b, and c have a high correlation. However, the overall trend of d is very different from that of a, b, and c (i.e., the correlation is very low). Therefore, we can divide information segment d into a different category from that of information segments a, b, and c.

Therefore, we can transform the problem of finding a single mutation point into that of finding an information segment including the mutation point; therefore, the concept of the correlation coefficient is introduced in this study. The correlation coefficient is a quantity that reflects the degree of correlation between variables. For example, for the variables X and Y, the correlation coefficient is calculated as follows:

$${\rho }_{XY}=r(X,Y)=\frac{\mathrm{Cov}(X,Y)}{\sqrt{\mathrm{Var}[X]\times \mathrm{Var}[Y]}}$$

(23)

where $\mathrm{Cov}(X,Y)$ is the covariance of $X$ and $Y$, $\mathrm{Var}[X]$ is the variance of $X$, and $\mathrm{Var}[Y]$ is the variance of $Y$. The correlation coefficient, ${\rho }_{XY}\in \left[-1,1\right]$, reflects the correlation between $X$ and $Y$. The greater the value of $\left|{\rho }_{XY}\right|$, the higher the degree of correlation; the smaller the value of $\left|{\rho }_{XY}\right|$, the lower the degree of correlation.

Based on the characteristics of the correlation coefficient mentioned above, a simulated mutation array is established in this study. By calculating the correlation coefficients between each information segment of ${\mathit{h}}_{rel}$ and the simulated mutation array, the information segments before the mutation point can be distinguished from the information segments containing the information after the mutation point. Thus, the information segment containing the mutation information can be identified to determine the maneuver point. The specific calculation steps are described below.

Firstly, ${\mathit{h}}_{simulated}$, an N-width sliding window made of an array of numbers, is defined to simulate the sudden change in the data. It is noticed that the setting of ${\mathit{h}}_{simulated}$ not unique. The requirement for setting ${\mathit{h}}_{simulated}$ is that the trend of ${\mathit{h}}_{simulated}$ is obviously different from the gentle change trend of $\mathit{h}$ with no maneuvers (as mentioned after Equation (15)), and is as similar to the rapid change trend of $\mathit{h}$ after maneuvers. Therefore, in fact, any group of arrays showing rapid changes can be set as ${\mathit{h}}_{simulated}$, but it should be noted that when setting different ${\mathit{h}}_{simulated}$, ${\rho }_{threshold}$ will be different (the definition of ${\rho }_{threshold}$ will be given later).

Then, each series of the three axes components of the relative angular momentum in the sliding window can be used to calculate the correlation coefficient with the simulated mutation array $h_{simulated}$ as follows:

$$\left\{\begin{array}{l}\left|{\rho }_{h_X(i)}\right|=\left|\frac{\mathrm{Cov}(h_{simulated},h_{X,window})}{\sqrt{\mathrm{Var}[h_{simulated}]\times \mathrm{Var}[h_{X,window}]}}\right|\\ \left|{\rho }_{h_Y(i)}\right|=\left|\frac{\mathrm{Cov}(h_{simulated},h_{Y,window})}{\sqrt{\mathrm{Var}[h_{simulated}]\times \mathrm{Var}[h_{Y,window}]}}\right|\\ \left|{\rho }_{h_Z(i)}\right|=\left|\frac{\mathrm{Cov}(h_{simulated},h_{Z,window})}{\sqrt{\mathrm{Var}[h_{simulated}]\times \mathrm{Var}[h_{Z,window}]}}\right|\end{array}\right.$$

(24)

where $h_{X,window}$, $h_{Y,window}$, and $h_{Z,window}$ are the components of the estimated relative angular momentum in the sliding window.

Finally, the absolute values of the correlation coefficients, i.e., $\left|{\rho }_{h_{X/Y/Z}(i)}\right|$ (short for $\left|{\rho }_{h_X(i)}\right|$, $\left|{\rho }_{h_Y(i)}\right|$, and $\left|{\rho }_{h_Z(i)}\right|$) are used to compare with a designed threshold to determine the maneuver. In detail, if $\left|{\rho }_{h_{X/Y/Z}(i)}\right|>{\left|\rho \right|}_{threshold}$ is satisfied, the maneuver is deduced; on the contrary, it infers no maneuver. Since the measurements are continuously obtained, the above steps can be repeated. Moreover, the threshold ${\left|\rho \right|}_{threshold}$ is defined as the maximum of $\left|{\rho }_{h_{X/Y/Z}(i)}\right|$ when the target does not maneuver in a large number of simulation experiments, which changes corresponding to specific ${\mathit{h}}_{simulated}$.

5. Numerical Simulations

To verify the effectiveness of the proposed method, a standard Monte Carlo simulation system with two-body dynamics is created to model the dynamics of the orbital maneuver detection. As this study is focused on orbit maneuver detection, the attitude of the satellite is not discussed as it is directly assumed the attitude of the chaser satellite remains constant. Additionally, the effects of perturbations (J2 and higher-order gravity terms), vehicle mass and geometry (aerodynamic drag and solar radial pressure), and specific instrument error model structures are beyond the scope of this study as well.

Then, the truth models in the geocentric inertial frame for two spacecraft under two-body assumption can be given as follows:

$$\begin{array}{cc}\begin{array}{l}{\dot{\mathit{R}}}_c={\mathit{V}}_c\\ {\dot{\mathit{V}}}_c=-\frac{\mu }{{\left|{\mathit{R}}_c\right|}^3}{\mathit{R}}_c,\end{array}& \begin{array}{l}{\dot{\mathit{R}}}_t={\mathit{V}}_t\\ {\dot{\mathit{V}}}_t=-\frac{\mu }{{\left|{\mathit{R}}_t\right|}^3}{\mathit{R}}_t\end{array}\end{array}$$

(25)

where $\mu$ is the standard gravity constant, ${\mathit{R}}_c\in {\mathbb{R}}^3$ and ${\mathit{V}}_c\in {\mathbb{R}}^3$ are the position and velocity of the chaser satellite in the geocentric inertial frame, respectively, and ${\mathit{R}}_c\in {\mathbb{R}}^3$ and ${\mathit{V}}_c\in {\mathbb{R}}^3$ are the position and velocity of the target satellite in the geocentric inertial frame, respectively.

According to the measurement model in Equation (5), the state of the target relative to the chaser satellite is obtained.

$$\begin{array}{l}{\mathit{r}}_{rel}=C_e^{lvlh}\left({\mathit{R}}_t-{\mathit{R}}_c\right)\\ {\mathit{v}}_{rel}=C_e^{lvlh}\left({\mathit{V}}_t-{\mathit{V}}_c\right)-C_e^{lvlh}\left[\omega \times \right]\left({\mathit{R}}_t-{\mathit{R}}_c\right)\end{array}$$

(26)

where is $\omega$ the chaser satellite’s inertial angular velocity vector, $\left[\times \right]$ denotes the cross-product operator, $C_e^{lvlh}$ is the direction cosine matrix from the geocentric inertial frame to the LVLH frame which can be calculated from the orbit position and velocity vector. Finally, the true measurements are calculated from the true relative state and the camera offset as shown in Equation (6).

5.1. Error Models

It is assumed the attitude measurement error of the sensor satellite is made up of bias and Gaussian white noise. The bias for attitude is set to 0.001 rad/axis (X-axis, Y-axis, and Z-axis) while the standard deviation of the white noise is set to 0.001 rad/axis (X-axis, Y-axis, and Z-axis).

Meanwhile, as stated in Section 2.3, the line-of-sight angle measurement error is assumed to be Gaussian noise. Two values of the standard deviation of the LOS angle measurement noise are studied, i.e., σ_cam = 3 $\times$ 10⁻⁴ rad/axis or σ_cam = 3 $\times$ 10⁻⁵ rad/axis (azimuth-axis and elevation-axis). The specific selection of the sensor is described in a subsequent section.

5.2. Reference Mission and Trajectory

Because the proposed algorithm is developed for the application in high-altitude orbits and Geostationary Earth Orbit (GEO) is a kind of quite important high-altitude orbit, it is assumed that the chaser is orbiting near GEO as shown in Table 1, while two types of target’s orbit are set in terms of the relative range. The first type is the long-range case that is beyond one hundred kilometers, as shown in Table 2, where the target orbit is set by the differences with the parameters shown in Table 1 (the eccentricity and right ascension are the same). The second type is the kilometer-level close-range case as shown in Table 3, where four classical motions including the Vbar Stationary, Co-elliptic, Football, and Oscillating orbit for inspection are presented.

5.3. Computation Models for Estimation of the Error

Since the LOS angles and camera offset are measured values with noise, Monte Carlo simulation is employed to evaluate the accuracy of these statistics. The error models presented in Section 5.1 are then used to generate a set of measurements that are processed by the proposed algorithm to determine the maneuver for the i-th Monte Carlo sample.

To verify the validity and test the performance of the proposed algorithm, two parameters will be checked. The first one is the accuracy of the maneuver determination, defined by the true positive rate in all of the Monte Carlo runs.

$$P_m=\frac{m}{n}\times 100\%$$

(27)

where $m$ represents number of the correct determination cases while $n$ is the number of Monte Carlo runs. A total of 200 runs are selected for the Monte Carlo simulations which can roughly lead to 90% confidence (Buckland, 1984) [42].

The second one is the estimated error of the maneuver epoch for the determined cases, defined as follows:

$$e_{true}=\frac{1}{n}{\displaystyle \sum _{j=1}^n\left|t_e-t_r\right|}$$

(28)

where $t_e$ is the estimated maneuver epoch and $t_r$ is the reference one for the $j$-th Monte Carlo sample.

6. Performance Analysis

6.1. Settings of Key Parameters

As the calculation steps mentioned in Section 4, ${\mathit{h}}_{simulated}$ and ${\rho }_{threshold}$ need to be set at first.

As mentioned in the steps, ${\mathit{h}}_{simulated}$ is not unique, but after a large number of experimental attempts, it is found that the growth section of the quadratic parabola meets the requirements well, so this paper sets ${\mathit{h}}_{simulated}$ as follows:

$$h_{simulated}=[1^2,2^2,3^2\dots ,N^2]$$

(29)

And ${\rho }_{threshold}$ is set to 0.989 in this paper after 200 Monte Carlo runs corresponding to the setting of ${\mathit{h}}_{simulated}$ in Equation (22).

6.2. Results of Long-Range Cases

6.2.1. Advantage of the Relative Angular Momentum

To verify the advantage of using the relative angular momentum for maneuver detection, a single-impulse maneuver with the magnitude of 0.1 m/s and random direction is executed on the epoch of 3000th second by the three targets in Table 2. The ratio of the change in the angular momentum within 1 s after the maneuver of the target satellite relative to the angular momentum before the maneuver is presented in Table 4, where the results are the statistics of 200 Monte Carlo runs.

It can be seen that the relative change in the angular momentum is close to 0%, while the variation in the relative angular momentum is greater than 0.2%. The relative variation in the relative angular momentum is much larger than that of the angular momentum, which shows that the sensitivity of the relative angular momentum to orbital maneuvers is much higher than that of angular momentum.

To verify that the sensor accuracy requirements mentioned in Section 3.2 are applicable to remote targets, Target 1 is taken as an example. A single-impulse maneuver with a magnitude of 0.1 m/s and random direction is executed at t = 1000 s for the case of σ_cam = 3 $\times$ 10⁻⁵ rad/axis. Through the second-order polynomial fitting of the data before the maneuver, the expected values before and after the maneuver are obtained and compared to the measured values, as shown in Figure 5.

It can be seen that the measured value fluctuates around the expected value for the reason of sensor noise when the maneuver does not occur. After t = 1000 s, the measured value deviates markedly from the expected value, and the deviation magnitude is significantly greater than that caused by noises. Therefore, the change in the relative angular momentum caused by the maneuver can be easily distinguished. This verifies the correctness of the theoretical derivation results presented in Section 3.2.

Moreover, Target 1 and Target 2 are selected to show the estimated performance of the relative angular momentum for an orbit maneuver. Figure 6 shows the theoretical values of Target 1’s relative angular momentum after 0.1 m/s impulse maneuver at t = 3500 s and the corresponding estimation with noises is shown in Figure 7, while the results for 1 m/s impulse on Target 2 are shown in Figure 8 and Figure 9.

As shown in Figure 6 and Figure 8, the relative angular momentum exhibits an obvious mutation after maneuver which verifies that the relative angular momentum has very good observability for the target orbit maneuver. The changes in the relative angular momentum values in Figure 6 and Figure 9 are somewhat different after maneuvers because the filter for the estimation of relative position and velocity mentioned in Section 2.1 deals with the change caused by the orbital maneuver as noise, but they are generally similar and have the same characteristics. In Figure 6 and Figure 8, when there is no maneuver before t = 3500 s, the relative angular momentum exhibits a trend similar to a straight line after the filter converges smoothly. At t = 3500 s, the relative angular momentum in the Y-direction of the LVLH system has a larger mutation than the original foundation. After the mutation, the relative angular momentum values in the three directions tend to be a stable straight line. This shows that the relative angular momentum still has excellent observability for the orbital maneuvers of long-range targets when the noises are considered.

Compared with Figure 8, the curve in Figure 7 is smoother and closer to the corresponding theoretical curve, which is caused by the difference in the measurement accuracy of the sensor. This shows that the accuracy of the sensor will affect the observation performance, and the higher the accuracy, the better the observation performance achieved.

6.2.2. Maneuver Detection Error and the True Positive Rate

As discussed in the previous section, the estimation of the relative angular momentum is polluted by the sensor noise. Thus, an appropriate width sliding window is quite important for matching with the pre-defined mutation simulation array during the detection. Too small or too large of a sliding window will affect the ability of the window containing data to reflect the overall characteristics of the actual inspection quantity. To test the performance of different width windows, the following simulation is conducted: Target 1 is selected to test, σ_cam = 3 × 10⁻⁴ rad/axis, and a randomly directed 1 m/s impulse is executed at t = 3000 s.

The results of 20 runs are shown in Figure 10 to present the distribution of maneuver detection errors in terms of different sliding widths. It can be seen that when the width of the sliding window is too small, the detection errors are relatively large. In these cases, the false positive rate is high because the estimation data in the narrow window cannot reflect the overall actual characteristics of the trend. As the width increases, the false positive rate significantly reduces. According to the performance shown in Figure 10, a sliding window with 50 widths is a good choice for the detection, which will be used for the subsequent simulations.

Next, the sensitivities of the proposed algorithm to the accuracy of the optical sensor, magnitude of the maneuver, and different relative orbit types are checked in the following simulations. A total of 200 Monte Carlo runs are conducted for each case, where the maneuver is executed at t = 3000 s. The simulation results are shown in Figure 11, Figure 12 and Figure 13 while the average errors of the 200 runs are presented in Table 5. The six points from left to right on the curve represent the probability that the detection error is less than 60 s, 120 s, 180 s, 240 s, 300 s, and 360 s, respectively.

By comparing the detection results for the same target with different observation sensor accuracies in Figure 11, Figure 12 and Figure 13, it can be concluded that the detection performance is improved with higher detection accuracy. Taking Target 1 in Figure 12 as an example, the red solid line is the detection result under the high-precision sensor, and the red dotted line is the detection result under the low-precision sensor. It can be seen in Figure 12 that the distance between the two curves is large initially and then decreases. This shows that with the high-precision sensor, the detection error has a greater probability of being concentrated in a place with a small value, and the overall error is smaller.

Figure 11, Figure 12 and Figure 13 also indicate that the maneuver magnitude has a significant influence on the performance. For small maneuvers (e.g., 0.2 m/s), the detection with the low-precision sensor is less effective, but it is more effective with the high-precision sensor. With the high-precision sensor, the detections of coplanar Target 1 and Target 2 are better than that of the heterogeneous Target 3. The true positive rate of Target 1 with a detection error of less than 60 s is 88.5% and that of less than 120 s is 98%. The true positive rate of Target 2 with a detection error of less than 120 s is 96%. The detection error of Target 3 is greater than that of the coplanar targets overall; however, for a detection error of less than 360 s, the true positive rate is 88%. When the maneuver magnitude is 0.6 m/s, the probability that the detection error of Target 1 and Target 2 is less than 60 s reaches 98.5%, and the true positive rate for a detection error of Target 3 of less than 360 s increases to 94%. When the maneuver magnitude reaches 1.0 m/s, the probability that the detection errors of Target 1 and Target 2 are less than 60 s is close to 100%. This shows that the method proposed in this study requires a high-precision measurement sensor for small maneuvers; otherwise, the detection performance will be reduced or lost. For large maneuver magnitude cases, the detection error of the high-precision measurement sensor is smaller overall. However, if the magnitude of the detection error is not considered, but only whether the target orbit maneuver is detected, there is little difference in the detection true positive rate. In this case, a low-precision measurement sensor with a lower cost can be selected.

Additionally, the results show that the detection performance for the targets in coplanar orbits such as Target 1 and Target 2 is better than for the targets in non-coplanar orbits, such as Target 3. The detection errors for the targets in coplanar orbits have a greater probability than those of the targets in coplanar orbits to fall in range with small values, so the overall errors for the targets in coplanar orbits are smaller. Moreover, a comparison of Target 1 and Target 2 shows that among the coplanar targets, the performance for the closer targets is better, the reason for which is that the orbit determination of closer targets is more accurate and stable, thus the variation in characterization, i.e., ${\mathit{h}}_{simulated}$ caused by orbit maneuvers is more obvious.

The above analysis demonstrates that distance and coplanarity can impact the detection performance. Comparing the results for Target 2, which is more distant but coplanar, with those for Target 3, which is closer but not coplanar, indicates that the negative impact of not being coplanar is greater than that of the larger distance.

Table 5 shows the average detection errors for different targets with different settings of sensors. As expected, errors were smaller with high-precision measurement sensors. With the increase in $∆\mathrm{V}$, the errors reduce, except for Target 3. The reason is that when selecting ${\mathit{h}}_{simulated}$ and ${\rho }_{threshold}$, we decided to try to detect all the maneuvers possible sacrificing the size of errors, so ${\mathit{h}}_{simulated}$ may not be most suitable considering the size of errors when $∆\mathrm{V}$ is a certain value. And this also explains why some other expected phenomena, such as the obvious differences in the average detection errors between the targets in different distances or between the targets in coplanar orbits and in non-coplanar orbits do not appear in Table 5.

Furthermore, a comparison of the performance between the proposed approach in this paper and the algorithm in studies by Liu et al. [33] and Wang et al. [34] is conducted. Liu et al. [33] present a maneuver detection scheme based on space-based LOS angle measurements. To make a fair comparison, a magnitude of maneuver impulse, i.e., 1 m/s impulse as same as the one used in the study by Liu et al., is executed on Target 1. The detection true positive rate by using Liu’s scheme is approximately 90%, which is less than the result of the proposed algorithm in this study.

Wang et al. [34] developed a maneuver-detecting method based on space-based Lidar, which provides excellent ranging performance (the standard deviation of the noise is 0.05 m). When the standard deviation of the noise increases to 0.5 m, Wang’s method maintains excellent detection capability. However, when it comes up to 5 m, the detection performance greatly reduces, and if the standard deviation of the noise comes up to 8 m or even larger, the method will fail. On the contrary, the proposed method in this study works well in the case that the relative position uncertainty is more than 8 m when σ_cam = 3 × 10⁻⁵ rad/axis.

6.2.3. Double Single-Impulse Maneuver Detection and Tracking

To demonstrate the multiple-maneuver tracking ability of the proposed method, the following simulations based on Target 1 are conducted. The first maneuver, i.e., randomly directed impulse with a magnitude of 0.4 m/s is executed at t = 2500 s, and the second maneuver, i.e., randomly directed impulse with a magnitude of 1 m/s is executed at t = 5000 s. Moreover, σ_cam is set to 3 × 10⁻⁴ rad/axis. The filter for the estimation of relative position and velocity is set to be restarted after a maneuver is determined.

The results of the estimated error of the relative position and velocity are shown in Figure 14 and Figure 15. It can be seen that both the maneuvers are detected, and the detection error is quite small. In detail, the first maneuver detection result is t = 2579 s with an error of 79 s, and the second detection result is t = 5058 s with an error of 58 s. After the filter is restarted, the estimation error is going to converge and then keep convergent, as indicated by the red curves. In contrast, the blue curves presenting the errors without the proposed method are divergent after maneuvers.

6.3. Results of Close-Range Cases

The application in close-range cases is tested in this section. A randomly pointed impulse maneuver with a 0.1 m/s magnitude executed at t = 600 s and σ_cam = 3 × 10⁻⁴ rad/axis is selected for the simulation. The results based on 200 Monte Carlo runs are summarized in Table 6, where the initial relative orbit is randomly selected for each run in the four types shown in Table 3.

As can be seen from Table 6, the maneuver of 93.5% of the runs is successfully estimated within 30 s after the maneuver, and 98.5% is estimated within 60 s. However, the true positive rate does not change when the error is greater than 60 s. It indicates that the remaining 1.5% of the runs were not detected because of the sensor noise.

The continuous tracking performance for a close-range target initialized in a V-bar stationary orbit (which is considered as the most disadvantageous trajectory for angles-only navigation) is tested. There are two randomly directed impulses executed at t = 500 s and t = 700 s with the magnitudes of 0.1 m/s and 1 m/s, respectively. The results are shown in Figure 16 and Figure 17, where the red curves present the estimated errors with the proposed method and the blue curves present the errors without the proposed method. It can be seen that the maneuvers are detected and the estimate error is convergent after the restarting of the filter when each maneuver is detected. The errors of the two detection results are 23 s and 29 s, respectively.

7. Conclusions

In this study, a novel space-based passive orbital maneuver detection algorithm for a non-cooperative space target is developed, where the concept of relative angular momentum is introduced to characterize the maneuver. A standard Monte Carlo simulation system with two-body dynamics was created to validate the proposed algorithm and evaluate the performance in the context of near-GEO orbit space. The accuracy of the detection algorithm was determined and presented for a variety of parameters and trajectories. When σ_cam = 3 × 10⁻⁵ rad/axis, the true positive rates approach 100%, much better than when σ_cam = 3 × 10⁻⁴ rad/axis. The average detection error is greater when the orbit of the chaser satellite and the target satellite are not in the same plane, but there is no obvious difference in the true positive rate. As expected, the algorithm worked better for the close-range cases than for the cases where the target was kilometers away. Surprisingly, the detection algorithm worked well for the long-range cases that were hundreds of kilometers away while the magnitude of the maneuver impulse was as small as 0.2 m/s. In the future, the detection algorithm for continuous low-thrust maneuvers will be studied.