Efficient Reachable Domain Search-Tracking for Cislunar Non-Cooperative Targets via Designed Quadrature

Abstract

To address the triple challenges of data sparsity, highly nonlinear dynamics, and maneuver uncertainty in tracking non-cooperative targets in cislunar space, we propose a collaborative framework combining Particle Filter (PF) and Unscented Kalman Filter (UKF). This framework optimizes search efficiency through a two-phase strategy: in the search phase, PF constructs the target reachable domain and leverages undetected information to dynamically shrink the search scope; upon target detection, the framework switches to UKF for high-precision and low-overhead tracking. To overcome the computational bottleneck in high-dimensional reachable domain integration, we integrate a non-product-type Designed Quadrature (DQ) method—one that generates minimal quadrature point sets to replace traditional Monte Carlo sampling by matching the moment conditions of mixed distributions via Gauss–Newton optimization. Distinct from existing single-filter or reachability modeling approaches, the key novelties of this work lie in a two-phase PF-UKF switching framework tailored to the unique cislunar environment resolving the trade-off between search capability and computational efficiency and integration of the non-product DQ method to break the dimensionality curse in high-dimensional reachable domain computation ensuring both moment-matching accuracy and real-time performance. This work holds potential to support space domain awareness and cislunar mission safety: reliable tracking of non-cooperative targets is a key prerequisite for avoiding collisions, safeguarding space assets, and enabling effective space defense, and the proposed framework provides a feasible technical path for this goal through simulation validation. Simulations demonstrate that on a three-dimensional Distant Retrograde Orbit (DRO) observation platform, successful recapture of cislunar transfer orbit targets can be achieved. Under fifth-order accuracy conditions, the system exhibits a position error of $3.745\times {10}^{-1}\mathrm{km}$ and a velocity tracking error of $9.703\times {10}^{-3}\mathrm{m}/\mathrm{s}$ for target search-and-tracking tasks, with a system response time of 1.8343 h. Compared with the traditional PF + numerical integration method, our proposed PF-UKF framework achieves an 86.7% reduction in time cost and a 24.1% reduction in position error.

1. Introduction

In recent years, cislunar space mission planning, exemplified by the Artemis Program and its lunar orbital platform “Lunar Orbital Platform-Gateway (LOP-G)”, has witnessed rapid development, marking humanity’s entry into a new era of deep-space exploration [1,2]. These mission architectures are highly complex and invariably traverse cislunar space, thereby significantly elevating the demand for Space Domain Awareness (SDA) capabilities, particularly in detecting, tracking, and discerning behavioral intentions of non-cooperative maneuvering targets [3,4]. However, the unique dynamical environment of cislunar space introduces new challenges for SDA technologies: due to limited sensor coverage and observational resource availability, prolonged monitoring of all targets is impractical, resulting in sparse and discontinuous measurement data [5]; targets are subject to the gravitational interaction of three bodies, leading to highly nonlinear characteristics in the evolution of system states [6]; unknown maneuvers executed by non-cooperative targets further amplify the complexity and uncertainty of state estimation problems [7].

Under the combined conditions of sparse measurements, strongly nonlinear dynamics, and pronounced maneuver uncertainty, conventional target tracking algorithms encounter substantial robustness challenges [8]. In the specific context of non-cooperative target tracking in cislunar space, the limited coverage and scarcity of sensor assets often produce observation gaps ranging from several hours to several days. During these data-starved intervals, unobserved target maneuvers can induce rapid deviations in the state trajectory, ultimately causing loss of sensor lock on the intended object. Consequently, there is an urgent imperative to develop novel theoretical frameworks and methodological approaches capable of delivering high-precision situational awareness for non-cooperative targets in the cislunar environment.

In the field of non-cooperative target tracking filtering algorithms, existing studies focus on diverse filtering frameworks, yet their adaptability fails to match the strong nonlinearity and high-dimensional uncertainty inherent to cislunar space. The Extended Kalman Filter (EKF) handles nonlinear systems via linearization and performs well in low-nonlinearity scenarios, but linearization errors easily trigger filter divergence in the strongly nonlinear Circular Restricted Three-Body Problem (CR3BP) of the Earth–Moon system [9]. The UKF avoids linearization defects using the Unscented Transformation (UT), improving nonlinear estimation accuracy, but it is inherently designed for target tracking with available observational data rather than initial large-scale searches. Lacking the capability to construct global reachable domains for blind search in data-sparse scenarios, it cannot accommodate the non-Gaussian uncertainties arising from unknown target maneuvers, thus limiting its adaptability in the search phase of cislunar non-cooperative target tracking [9,10]. The PF addresses non-Gaussian and strongly nonlinear problems, yet its “curse of dimensionality” is particularly pronounced in cislunar high-dimensional state spaces, with computational efficiency far below the real-time requirements of cislunar tasks [9,11]. Even improved algorithms, such as the robust information filter based on the maximum correntropy criterion [12], the Adaptive Consensus-Based Unscented Information Filter [10], the Cubature Kalman Filter (CKF) with adaptive noise estimation [13], and the Distributed Particle Filter (DPF) for multi-platform collaboration [14]—optimize robustness or distributed performance but fail to resolve the high-dimensional reachable domain modeling and dynamic search scope contraction issues specific to cislunar non-cooperative target tracking. For instance, Ye et al. [13] proposed an adaptive robust CKF to adjust process noise covariance for maneuvering targets, but it still relies on the Gaussian assumption and cannot handle cislunar non-Gaussian uncertainties. Han et al. [14] developed a DPF for multi-platform data fusion, yet its computational complexity accumulates with the number of platforms, making it unsuitable for real-time cislunar tasks.

Reachable domain computation, a core component of cislunar non-cooperative target search that complements filtering algorithms, faces a critical challenge in balancing accuracy and efficiency. The Reachability Set theory [15] introduces a novel paradigm for observation association in Space Situational Awareness (SSA) but lacks detailed characterization of internal reachable domain regions, limiting its applicability to refined cislunar target search tasks. Xue et al. [16] and Wen et al. [17] conducted analytical investigations on reachable domains for impulsive maneuvering spacecraft; however, they overlooked high-dimensional uncertainty propagation under the Earth–Moon’s complex gravitational field, failing to account for trajectory deviations induced by unknown maneuvers. Traditional Monte Carlo sampling requires extensive sample sizes to cover reachable domains, with computational complexity growing exponentially with dimensionality [11,18]. While the High-Order Sensitivity Matrix (HOSM) method [18,19] improves efficiency through non-product integration strategies, it suffers from insufficient moment matching accuracy for mixed distributions. Additionally, the traditional tensor product quadrature method [11] exhibits low computational efficiency, rendering it unsuitable for real-time cislunar missions. To address high-dimensional integration challenges, Hu et al. [20] proposed a Gaussian-measure sparse quadrature method, yet its Gaussian-specific design restricts applicability to cislunar mixed-distribution scenarios, leading to inadequate accuracy. For cislunar-specific reachable set approximation, Jones et al. [21] developed a method using pseudospectral techniques and homotopy, but it lacks generalization under the strong nonlinearity of CR3BP dynamics.

As a critical link between “reachable domain modeling” and “target detection,” existing sensor scheduling strategies are poorly adapted to the data sparsity and high-dimensional uncertainty of cislunar space. Li et al.’s Double Coverage Probability (DCP)-based sensor scheduling [22] accounts for orbital uncertainty and improves space debris tracking but ignores the high-dimensional reachable domain modeling and dynamic search needs specific to cislunar non-cooperative targets, preventing targeted allocation of observation resources. Adurthi et al.’s mutual information-based sensor task assignment [3] optimizes SSA observation efficiency but fails to utilize the negative information from undetected targets to dynamically shrink search scopes, leading to resource waste and target loss in cislunar scenarios with observation gaps of several hours [23]. To adapt to dynamic target states, Tang et al. [24] designed a deep reinforcement learning-based sensor scheduler, but it requires large-scale pre-training and cannot handle sudden maneuvers. For multi-satellite scheduling of moving targets, Qin et al. [25] proposed a method using an enhanced hybrid genetic simulated annealing algorithm and observation strip selection to balance observation conflicts and workload distribution, yet it does not integrate reachable domain information, resulting in suboptimal search coverage for cislunar targets.

To address these shortcomings in filtering adaptability, reachable domain computation, and sensor scheduling, this paper proposes a “filtering-modeling-scheduling” integrated solution based on DQ method for cislunar non-cooperative target search and tracking. First, a phased hybrid PF-UKF framework is developed: PF handles high-dimensional uncertainty and dynamically shrinks reachable domains during the search phase, while UKF enables high-precision, low-overhead estimation during the tracking phase—resolving the scenario adaptability limitations of single filters. Second, a non-product DQ method is integrated, which matches mixed distribution moment conditions via Gauss–Newton optimization and models high-dimensional reachable domains with minimal quadrature points, balancing computational efficiency and accuracy. Third, an adaptive grid search strategy based on reachable domain samples is designed, which updates search scopes using undetected information to adapt to cislunar data sparsity. Specifically, the hybrid PF-UKF framework’s two-phase design operates as follows: during the search phase, PF constructs reachable domains based on the distribution of unknown target maneuver parameters and performs global searches within these domains; when targets remain undetected, PF combines Bayesian updates and undetected target information to effectively shrink reachable domain scopes. Once targets are detected, the framework switches to UKF for tracking, avoiding the heavy computational burden of continuous large-scale reachable domain calculations during tracking. The detailed workflow is illustrated in Figure 1.

To contextualize the practical challenge of uncooperative target search and tracking in cislunar space, Figure 2 is included herein. As shown, an observation platform first estimates the target’s trajectory as the Predicted Orbit using prior data. However, the noncooperative target performs an unannounced maneuver at the Maneuver Point during an observation gap (common in cislunar missions due to limited sensor coverage), deviating from the Predicted Orbit to an Actual Orbit. This causes the platform to fail in detecting the target via the Predicted Orbit, making active search within the platform’s Observer Orbit essential to reestablish target lock and ensure subsequent tracking continuity.

2. Problem Formulation

This section defines the applicable orbital dynamical model, system state dynamic evolution model, and observation sensor model for the problem of searching and tracking non-cooperative spacecraft targets in cislunar space.

2.1. Dynamical Model

Under the framework of the CRTBP for the Earth–Moon system, considering the motion of a spacecraft with negligible mass in the gravitational fields of the Earth (mass $m_E$) and the Moon (mass $m_M$), the dimensionless potential function is defined as [26]

(1)$U\left(\mathit{r}\right)={\displaystyle \frac{1-\mu }{r_E}}+{\displaystyle \frac{\mu }{r_M}}+{\displaystyle \frac{1}{2}}(x^2+y^2)$

where $\mathit{r}={[x,y,z]}^T$ is the target’s position vector, $\mu =m_M/(m_E+m_M)$ represents the Moon’s mass fraction, $r_E=\sqrt{{(x+\mu )}^2+y^2+z^2}$ and $r_M=\sqrt{{(x+\mu -1)}^2+y^2+z^2}$ denote the Euclidean distances from the spacecraft to Earth and Moon, respectively. The spacecraft’s equations of motion are given by

(2)$\ddot{\mathit{r}}=\nabla U\left(\mathit{r}\right)+2\left[\begin{array}{c}\dot{y}\\ -\dot{x}\\ 0\end{array}\right]+\left[\begin{array}{c}x\\ y\\ 0\end{array}\right]+\left[\begin{array}{c}{u}_x\left(t\right)\hfill \\ u_y\left(t\right)\hfill \\ u_z\left(t\right)\hfill \end{array}\right]$

where $\mathit{u}\left(t\right)={[u_x\left(t\right),u_y\left(t\right),u_z\left(t\right)]}^T$ is the control input vector, and $\nabla U\left(\mathit{r}\right)$ represents the gradient of the potential function.

2.2. System Configuration

Section 2.1 establishes the basic dynamical model for target motion. Section 2.2 refines this model for the specific scenario of “spacecraft-non-cooperative target tracking”: it splits spacecraft motion into two parts—unperturbed orbital dynamics and control input effects. We consider a system made up of an observation platform and a non-cooperative target spacecraft. Their continuous-time dynamics are governed by the following ordinary differential equations.

(3)$\dot{\mathit{x}}=\mathit{f}(\mathit{x},t)+\mathit{g}(\mathit{u},t)$

where $\mathit{x}\in {\mathbb{R}}^n$ represents the state vector, $\mathit{f}(·)$ characterizes the unperturbed orbital dynamics, and $\mathit{g}(·)$ models the control input effects. For the non-cooperative target spacecraft, the control input $\mathit{u}$ is a priori unknown and is assumed to consist of a finite number M of impulsive maneuvers occurring at unknown epochs ${\left\{t_i\right\}}_{i=1}^M$ with corresponding unknown magnitudes and directions ${\left\{{\mathit{u}}_i\right\}}_{i=1}^M$. Consequently, the target’s control input set is expressed as $\mathit{u}=\{({\mathit{u}}_1,t_1),({\mathit{u}}_2,t_2),\dots ,({\mathit{u}}_M,t_M)\}$.

The target’s control input is mathematically represented as a linear combination of Dirac delta functions:

(4)$\mathit{g}(\mathit{u},t)=\sum _{i=1}^M\mathit{g}({\mathit{u}}_i,t)\delta (t-t_i)$

For numerical implementation and state estimation purposes, the continuous-time dynamical system is discretized over a sequence of time instants $\left\{t_k\right\}$. Assuming the existence of a state transition function from the initial epoch $t_0$ to the current time $t_k$, the discrete-time system evolution is described by:

(5)${\mathit{x}}_k=\mathsf{\chi }({\mathit{x}}_0,\mathit{u},t_k)$

where ${\mathit{r}}_0$ denotes the initial state at time $t_0$. The trajectory evolution is completely determined by the initial state conditions and the discrete input sequence $\mathit{u}$.

To facilitate parameter estimation and optimization procedures, all unknown parameters and initial conditions influencing the system trajectory are consolidated into an augmented state vector:

(6)${\mathit{z}}^T=\left[{\mathsf{r}}_0^T,{\mathit{u}}_1^T,t_1,{\mathit{u}}_2^T,t_2,\dots ,{\mathit{u}}_M^T,t_M\right]$

For enhanced numerical conditioning and computational stability, the augmented vector $\mathit{z}$ is typically normalized prior to estimation procedures.

2.3. Sensor Model

The geometric relationship between the sensor field-of-view (FOV), target position, and detection likelihood function is shown in Figure 3. When the target is within the FOV (${\mathbf{C}}_s<0$), the observation models for azimuth angle A and elevation angle E are expressed as

(7)$\left\{\begin{array}{cc}A\hfill & ={tan}^{-1}\left({\displaystyle \frac{{\rho }_y}{{\rho }_x}}\right)+v_1\hfill \\ E\hfill & ={tan}^{-1}\left({\displaystyle \frac{{\rho }_z}{\sqrt{{\rho }_x^2+{\rho }_y^2}}}\right)+v_2\hfill \end{array}\right.$

where $\mathsf{\rho }={\mathit{r}}_{tar}-{\mathit{r}}_{ob}={[{\rho }_x,{\rho }_y,{\rho }_z]}^T$ is the relative position vector, $\mathit{v}={[v_1,v_2]}^T$ is a measurement noise vector following a zero-mean Gaussian distribution.

The detection capability of sensors is constrained by their FOV. The FOV constraints are defined as

(8)${\mathit{C}}_s(\mathit{x},\mathsf{\theta })=\left\{\begin{array}{c}-\mathsf{\rho }·\mathit{a}\hfill \\ arctan\left({\displaystyle \frac{{\rho }_x}{\mathsf{\rho }·\mathit{a}}}\right)-T{H}_1\hfill \\ arctan\left({\displaystyle \frac{{\rho }_y}{\mathsf{\rho }·\mathit{a}}}\right)-T{H}_2\hfill \end{array}\right.$

where $\mathit{x}$ is the target state vector, $\mathsf{\theta }$ denotes the sensor pointing parameter vector, $T{H}_1$ and $T{H}_2$ represent the horizontal and vertical semi-FOV angles, and $\mathbf{a}$ is the sensor boresight direction. When ${\mathbf{C}}_s<0$, the target resides within the sensor FOV.

The piecewise detection likelihood function is given by

(9)${\pi }_d\left({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1}\right)=\left\{\begin{array}{cc}{\pi }_d^{\prime }\left({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1}\right)& \mathrm{if}{\mathit{C}}_s\left({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1}\right)<0\\ 0& \mathrm{otherwise}\end{array}\right.$

where ${\pi }_d^{\prime }\left({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1}\right)\in [0,1]$ quantifies the probability of successfully detecting a target with state ${\mathit{x}}_{k+1}$ under sensor configuration ${\mathsf{\theta }}_{k+1}$.

The detection optimization criterion aims to maximize the expected detection likelihood

(10)$\underset{{\mathsf{\theta }}_{k+1}}{max}{J}_d=E\left[{\pi }_d\left({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1}\right)\right]=\int {\pi }_d\left({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1}\right)\pi \left({\mathit{x}}_{k+1}\right)d\mathit{x}$

where $\pi \left({\mathit{x}}_{k+1}\right)$ denotes the probability density function of the system state, characterizing prior uncertainty in the target state.

3. Reachability Set Computation via Designed Quadrature

Computationally, direct reachable domain evaluation is prohibitively expensive, requiring explicit computation of numerous MC samples to accurately represent the search space. This section introduces the fundamental concepts and equations of the HOSM method, which efficiently approximates the reachable domain through polynomial modeling.

3.1. High-Order Sensitivity Matrix Method

Traditional MC methods demand individual dynamic model integration for each sample, incurring high computational costs. The HOSM method constructs polynomial models as surrogates for direct dynamics integration. While polynomial approximations introduce marginal errors compared to numerical integration, they significantly enhance computational efficiency.

The Taylor expansion of a function $\mathsf{\chi }\left(\mathsf{\zeta }\right)$ at $\mathsf{\zeta }=0$ (truncated to order d) is

(11)$\mathsf{\chi }\left(\mathsf{\zeta }\right)\approx \sum _{k=0}^d{\displaystyle \frac{1}{k!}}{(\mathsf{\zeta }·\nabla )}^k\mathsf{\chi }\left(0\right)$

where $\mathsf{\zeta }$ is a normalized random input vector with known probability density function (pdf), k denotes the expansion order, ${\displaystyle \frac{1}{k!}}$ is the factorial coefficient of the Taylor expansion, ${(\mathsf{\zeta }·\nabla )}^k$ represents the k-th order parameter-gradient operator, and $\mathsf{\zeta }·\nabla ={\mathsf{\zeta }}_{{\alpha }_1}{\displaystyle \frac{\partial }{\partial {\mathsf{\zeta }}_{{\alpha }_1}}}+{\mathsf{\zeta }}_{{\alpha }_2}{\displaystyle \frac{\partial }{\partial {\mathsf{\zeta }}_{{\alpha }_2}}}+\dots$ is a linear combination of parameters and partial derivatives.

The HOSM method numerically computes an approximation analogous to the Taylor series without explicitly evaluating partial derivatives of $\mathsf{\chi }\left(\mathsf{\zeta }\right)$. All quadrature methods adopt the same fundamental approach: approximating a function as a Taylor polynomial and integrating this polynomial. All quadrature methods approximate polynomial integration using weighted sums evaluated at specific nodes ${\mathsf{\zeta }}^{\left(i\right)}$ [27]

(12)$E\left[{\mathsf{\chi }}_j\left(\mathsf{\zeta }\right)\right]\approx \sum _{i=1}^N{w}^{\left(i\right)}{\mathsf{\chi }}_j\left({\mathsf{\zeta }}^{\left(i\right)}\right)$

Computational cost scales with the number of evaluation points N. Substituting the Taylor expansion yields the Moment Constraint Equations (MCEs)

(13)$\begin{array}{cc}\hfill E\left[1\right]& =\sum _{i=1}^N{w}^{\left(i\right)}\hfill \\ \hfill E\left[{\mathsf{\zeta }}_{{\alpha }_1}\right]& =\sum _{i=1}^N{w}^{\left(i\right)}{\mathsf{\zeta }}_{{\alpha }_1}^{\left(i\right)}\hfill \\ \hfill E\left[{\mathsf{\zeta }}_{{\alpha }_1}{\mathsf{\zeta }}_{{\alpha }_2}\right]& =\sum _{i=1}^N{w}^{\left(i\right)}{\mathsf{\zeta }}_{{\alpha }_1}^{\left(i\right)}{\mathsf{\zeta }}_{{\alpha }_2}^{\left(i\right)}\hfill \\ \hfill \dots \\ \hfill E\left[{\mathsf{\zeta }}_{{\alpha }_1}{\mathsf{\zeta }}_{{\alpha }_2}\dots {\mathsf{\zeta }}_{{\alpha }_d}\right]& =\sum _{i=1}^N{w}^{\left(i\right)}{\mathsf{\zeta }}_{{\alpha }_1}^{\left(i\right)}{\mathsf{\zeta }}_{{\alpha }_2}^{\left(i\right)}\dots {\mathsf{\zeta }}_{{\alpha }_d}^{\left(i\right)}\hfill \end{array}$

For non-cooperative targets, the normalized input vector $\mathsf{\zeta }$ follows a hybrid distribution combining two components. The maneuver time $t_m$ is uniformly distributed within the interval $\left[t_a,t_b\right]$, denoted as $t_m\sim \mathcal{U}(t_a,t_b)$; The maneuver velocity $\Delta \mathit{v}$ is spherically uniformly distributed with directional isotropy and a magnitude constraint $\left[0,\Delta v_{max}\right]$, expressed as $\Delta \mathit{v}\sim {\mathcal{U}}_s(0,\Delta v_{max})$. This model imposes bounded uncertainty constraints under the assumption of single-impulse maneuvers, where all directions and magnitudes (up to $\Delta v_{max}$) are equiprobable.

To overcome dimensionality inflation, we introduce the non-product-type DQ method, which optimizes nodes ${\left\{{\mathsf{\zeta }}_i\right\}}_{i=1}^N$ and weights ${\left\{w_i\right\}}_{i=1}^N$ to match moments with minimal points N.

(14)$\sum _{i=1}^N{w}_i\varphi \left({\mathsf{\zeta }}_i\right)=\int \varphi \left(\mathsf{\zeta }\right)\pi \left(\mathsf{\zeta }\right)d\mathsf{\zeta }$

where $\varphi (·)$ is a polynomial basis function, $w_i$ is a weight, and $\pi \left(\mathsf{\zeta }\right)$ is a mixed probability density function. This method optimally selects nodes ${\left\{{\mathsf{\zeta }}_i\right\}}_{i=1}^N$ and weights ${\left\{w_i\right\}}_{i=1}^N$, to accurately match the specified order moments with the minimum number of points N.

3.2. Designed Quadrature Algorithm

A concise overview of the DQ algorithm is presented herein; comprehensive theoretical discussions may be found in [28]. The algorithm calculates point sets that meet specified moment targets and constraints. As an iterative technique, DQ minimizes an augmented cost function through integration with Gauss–Newton optimization. It gradually reduces the cardinality of the point set until no smaller effective set can be formed. This method focuses on solving numerical integration problems for mixed distributions, with its core procedure outlined in Figure 4. The computational cost of the quadrature scales linearly with N. Although recursive DQ adds initial overhead for identifying minimal quadrature, the resulting set significantly improves efficiency in reachable domain computations.

The key advantage of the algorithm lies in its iterative handling of moment constraints through the Gauss–Newton framework. It directly optimizes node positions and weights without relying on symmetry assumptions, thereby fully preserving odd-order moment statistics. Using only known statistical moments and support constraints, DQ efficiently generates minimal point sets—this significantly speeds up numerical integration for high-dimensional mixed distributions.

4. Hybrid PF-UKF Tracking Framework

4.1. Sensor Parameter Optimization

After constructing the target reachable domain propagation model based on the HOSM method, the sensor task is transformed into an optimal search problem for this reachable domain. The core optimization objective is to maximize the expected value of the detection likelihood

(15)$\underset{\mathsf{\theta }}{max}{J}_d=E\left[{\pi }_d(\mathit{x},\mathsf{\theta })\right]={\int }_{\mathsf{\Omega }}{\pi }_d(\mathit{x},\mathsf{\theta })\pi \left(\mathit{x}\right)d\mathit{x}$

Due to the axisymmetric characteristics of the pyramid-shaped FOV and detection probability, the sensor pointing can be completely described by the azimuth angle ${\mathsf{\theta }}_1$ and elevation angle ${\mathsf{\theta }}_2$. An adaptive grid search strategy is proposed, which uses the sample set $\left\{{\mathit{x}}^{\left(i\right)}\right\}$ in the reachable domain to determine the azimuth angle range $\left[minA\left({\mathit{x}}^{\left(i\right)}\right),maxA\left({\mathit{x}}^{\left(i\right)}\right)\right]$ and the elevation angle range $\left[minE\left({\mathit{x}}^{\left(i\right)}\right),maxE\left({\mathit{x}}^{\left(i\right)}\right)\right]$, and constructs grid boundary constraints. The incremental search grid is generated as

(16)$\begin{array}{cc}\hfill {\mathsf{\Theta }}_1& =\left[A_{min},A_{min}+T{H}_1,\dots ,A_{max}\right]\hfill \\ \hfill {\mathsf{\Theta }}_2& =\left[E_{min},E_{min}+T{H}_2,\dots ,E_{max}\right]\hfill \\ \hfill \mathsf{\Phi }& ={\mathsf{\Theta }}_1\otimes {\mathsf{\Theta }}_2\hfill \end{array}$

This method avoids the computational burden of global grid search and improves efficiency while ensuring optimality. After determining the optimal pointing parameter ${\mathsf{\theta }}^{\ast }$, observations are made on the reachable domain.

4.2. Target Tracking

In target tracking tasks, the traditional filtering problem focuses on fusing measurement data with propagated state uncertainty to obtain the posterior distribution

(17)$\pi \left({\mathit{x}}^{+}\right)\propto \pi \left({\mathit{x}}^{-}\right)\pi (\tilde{\mathit{y}}\mid {\mathit{x}}^{-})$

where $\pi \left({\mathit{x}}^{-}\right)$ is the prior probability density function, and $\pi (\tilde{\mathit{y}}\mid {\mathit{x}}^{-})$ is the measurement likelihood function. This paper focuses on how to use state uncertainty to guide task sensors for target detection, with the switching between PF and UKF being a core design to balance search capability and tracking precision.

Specifically, the switching between PF and UKF is triggered by two mutually complementary conditions to ensure reliability: (1) the target detection probability ${\pi }_d(x_{k+1},{\theta }_{k+1})\ge 0.9$ (defined in Equation (9)), indicating the target is stably within the sensor’s FOV; (2) the Effective Sample Size (ESS) of the PF particle set $\ge N/2$ (where $N={10}^4$ is the initial number of particles), ensuring the PF state estimation is sufficiently reliable to initialize the UKF. When both conditions are satisfied, the framework switches to UKF for high-precision tracking; otherwise, it maintains PF for reachable domain construction and global search.

The switching execution process is implemented as follows: First, during state handover, the posterior state estimate of the PF and its corresponding covariance matrix are, respectively, given by

(18)$\begin{array}{cc}\hfill {\widehat{\mathit{x}}}_{\mathrm{PF}}& =\sum _{i=1}^N{w}^{\left(i\right)+}{\mathit{x}}^{\left(i\right)+}\hfill \\ \hfill {\mathit{P}}_{\mathrm{PF}}& =\sum _{i=1}^N{w}^{\left(i\right)+}\left({\mathit{x}}^{\left(i\right)+}-{\widehat{\mathit{x}}}_{\mathrm{PF}}\right){\left({\mathit{x}}^{\left(i\right)+}-{\widehat{\mathit{x}}}_{\mathrm{PF}}\right)}^T\hfill \end{array}$

where ${\widehat{\mathit{x}}}_{\mathrm{PF}}$ denotes the posterior state estimate and ${\mathit{P}}_{\mathrm{PF}}$ represents the associated covariance matrix. These two quantities are directly adopted as the initial prior state ${\widehat{\mathit{x}}}_{\mathrm{UKF},0}$ and initial covariance matrix ${\mathit{P}}_{\mathrm{UKF},0}$ of the UKF, respectively. To enhance the numerical stability of the UKF initialization, a small identity matrix $ϵ\mathit{I}$ (with $ϵ$ being a tiny positive scalar) is added to ${\mathit{P}}_{\mathrm{UKF},0}$, leading to

(19)${\mathit{P}}_{\mathrm{UKF},0}={\mathit{P}}_{\mathrm{PF}}+{10}^{-6}\mathit{I}$

where $\mathit{I}$ is the 6 × 6 identity matrix corresponding to 3D position and velocity states. Second, for UKF parameter initialization, the UT parameters are set as $\alpha =0.001$, $\beta =2$, and $\kappa =3-n$ (with $n=6$ as the state dimension), which are optimal for orbital dynamics estimation tasks [9]. Additionally, a back-switching condition is set: if the target detection probability ${\pi }_d({\mathit{x}}_{k+1},{\mathsf{\theta }}_{k+1})<0.1$ for three consecutive observation steps, the framework switches back to PF to reconstruct the reachable domain and resume global search, avoiding tracking divergence due to temporary target loss.

This paper adopts a particle filter framework for full-density estimation. As a subset of Sequential Monte Carlo (SMC) technology, this method approximates the probability density as a discrete sample set

(20)$\pi \left(\mathit{x}\right)\approx \sum _{i=1}^N{w}^{\left(i\right)}\delta ({\mathit{x}}^{\left(i\right)}-\mathit{x})$

Considering the prior MC system representation, the sensor parameter $\mathsf{\theta }$ updates the weights through the measurement value $\tilde{y}$

(21)${\widehat{w}}^{{\left(i\right)}^{+}}\propto w^{{\left(i\right)}^{-}}\pi (\tilde{\mathit{y}}\mid {\mathit{x}}^{{\left(i\right)}^{-}})$

The normalized posterior weights are given by

(22)$w^{{\left(i\right)}^{+}}={\displaystyle \frac{{\widehat{w}}^{{\left(i\right)}^{+}}}{{\displaystyle \sum _{i=1}^N}{\widehat{w}}^{{\left(i\right)}^{+}}}}$

If the target detection fails, the measurement likelihood function is defined as the complementary function of the detection probability

(23)$\pi (\tilde{\mathit{y}}\mid {\mathit{x}}^{-})=1-{\pi }_d({\mathit{x}}^{-},\mathsf{\theta })$

This mechanism ensures that samples outside the sensor’s field-of-view have their weights attenuated as detection probability rises. When the target is successfully detected, the UKF activates for tracking purposes. This algorithm linearizes nonlinear systems via the Unscented Transform (UT).

In cases where the target remains undetected, the particle filter tends to suffer from sample impoverishment during sequential iterations. When the particle measurement likelihood function value $\pi \left(\tilde{\mathit{y}}|{\mathit{x}}^{\left(i\right)-}\right)$ approaches zero, its posterior weight ${\omega }^{\left(i\right)+}$ decreases sharply, leading to the concentration of the reachable domain estimation on very few particles. To overcome this defect, the systematic resampling algorithm is used to update the particle set ${\left\{\left({\mathsf{\zeta }}_i,{\omega }_i\right)\right\}}_{i=1}^N({\displaystyle \sum _{i=1}^N}{\omega }_i=1)$.

The ESS criterion is used to trigger resampling

(24)$ESS={\left(\sum _{i=1}^N{\left({\omega }^{\left(i\right)}\right)}^2\right)}^{-1}$

The resampling process is activated when $ESS<N_t$. In this framework, the threshold $N_t=N/2$ is set, which significantly reduces the computational overhead of meaningless resampling.

The above framework combines the HOSM method and the maximum detection likelihood sensor task allocation, providing a feasible framework for propagating and searching the reachable domain of non-cooperative target satellites. The following sections will present numerical simulations and discussion results.

5. Numerical Simulations

Taking the single-impulse orbital maneuver of the cislunar transfer hybrid orbit as an example, the target is searched and tracked. Among them, the cislunar transfer hybrid orbit is also referred to as the circumlunar free-return orbit + mid-course orbital maneuver strategy. It means that after the target spacecraft departs along the circumlunar free-return orbit for a certain period of time, all systems of the spacecraft are tested and confirmed to work normally, and a new lunar transfer orbit is entered by applying a mid-course orbital maneuver.

In cislunar space, the DRO has high stability; however, the DRO and the cislunar transfer orbit are basically in the same orbital plane, resulting in poor target observability. Therefore, this paper deploys the observation platform in a three-dimensional DRO.

The observation platform adopts optical measurement. It is assumed that the measurement noise follows a Gaussian distribution, with an angle measurement error of 2 arcsec. The sensor’s FOV is pyramid-shaped, where both the horizontal and vertical half-FOV angles are 1.5°. The effective detection distance is 0.6 times the cislunar space range, which is used to calculate the detection probability, as shown in Table 1.

Assuming no observed target, after 12 h, the maneuver time is defined as uniformly distributed between 0 and 12 h, i.e., $t_i\in \mathcal{U}(0,12\mathrm{h})$. Starting from $t=12$ h, measurements are taken every 10 min ($\Delta t=600\mathrm{s}$)until the simulation time ends at $t_f=24$ h, as shown in Table 2. It is also assumed that the target has a maximum maneuvering capability of 100 m/s, so the unknown maneuver follows a spherical uniform distribution $u_{t,i}\in {\mathcal{U}}_s(0,100\mathrm{m}/\mathrm{s})$, as shown in Table 3. The initial states of the target and the platform at $t=0$ are given by the state vectors in the Earth–Moon centroid rotating coordinate system in Table 4.

The search and tracking scenario for targets in Earth–Moon transfer orbits is illustrated in Figure 5, where the coordinate axis unit ’L’ corresponds to the average Earth–Moon distance (L$=3.844\times {10}^5\mathrm{km}$). In this configuration, the non-cooperative target is assumed to have entered the latter segment of its Earth–Moon transfer trajectory.

5.1. Method Validation

Under the second-order accuracy condition, in the search and tracking scenario, after the target is found, the target reachable domain and the platform’s field of view pointing are shown in Figure 6.

Table 5 and Figure 7 present the search and tracking results under the second-order accuracy condition. The orange segment represents the phase where the platform is searching within the target’s reachable domain but has not yet detected the target, while the blue segment denotes the tracking error curve after the target is detected. Owing to the relatively low second-order accuracy, the target search duration is prolonged, and the tracking error remains substantial.

Under the fifth-order accuracy condition, in the search and tracking scenario, after the target is found, the target reachable domain and the platform’s field of view pointing are shown in Figure 8.

Table 6 and Figure 9 illustrate the search and tracking results under the fifth-order accuracy condition. Specifically, under the fifth-order accuracy condition, the position error is reduced by 84.04%, the velocity error is decreased by 66.55%, and the system response speed is improved by 56.00%.

5.2. Method Comparison

To further clarify the innovation and superiority of the PF-UKF framework proposed in this paper, comparative simulations with existing related methods are added. The comparison objects include the traditional PF + numerical integration method, the PF + DQ method that only combines reachable domain modeling, and the PF-UKF framework with different accuracies (second-order and fifth-order) proposed in this paper. All simulation scenarios are consistent with those in the previous sections to ensure the comparability of results.

For the comparative simulation, a unified hardware platform is adopted to eliminate the impact of hardware differences on time cost: specifically, a 12th Gen Intel(R) Core(TM) i5-12500H CPU, 16 GB memory, and an Intel(R) Iris(R) Xe Graphics GPU are used as the standard configuration. Meanwhile, the evaluation indicators are defined as follows: the time cost refers to the total calculation time from the start of the simulation to the end of target tracking, which is recorded via the tic-toc function in Matlab; each method is run independently 100 times, and the average value is taken as the final result to reduce random errors. The definitions of position error and velocity error remain consistent with those specified in the previous sections.

The search and tracking performance of the traditional PF + numerical integration method and the PF + DQ method in non-cooperative target search-tracking scenarios is visualized in Figure 10, respectively: Figure 10a illustrates the trajectory deviation and the system response speed of the PF + numerical integration method, while Figure 10b shows the large position fluctuation and slow error attenuation of the PF + DQ method. To quantify these performance differences, the comparative results of core indicators across all methods are summarized in Table 7. Notably, the fifth-order accuracy PF-UKF framework proposed in this paper achieves the optimal performance in key evaluation dimensions.

To summarize, the fifth-order PF-UKF framework reduces computation time by 86.7% relative to the traditional PF + numerical integration method, and outperforms the PF + DQ method in efficiency—validating the PF-UKF switching mechanism’s acceleration effect. In accuracy, its position error is smaller than both the PF + numerical integration method and the second-order PF-UKF, while its velocity error is the smallest across all methods. The PF + DQ method’s extreme errors stem from DQ’s inaccurate dynamic model fitting: this degrades particle quality in the PF process, leading to large position fluctuations and slow error attenuation. Overall, the fifth-order PF-UKF optimally balances efficiency and accuracy, fully demonstrating the practicality of the “reachable domain modeling + dual-filter switching” design.

To quantitatively verify the accuracy-efficiency trade-off of polynomial approximation models with different orders, additional simulations are conducted to compare the computational overhead and tracking performance of the DQ algorithm under second-order, third-order, fourth-order, and fifth-order accuracy. All simulations adopt the same hardware platform and simulation parameters as Section 5.2, with 100 independent runs for each order to ensure statistical reliability. The computational time costs are visually illustrated in Figure 11 and summarized in Table 8.

As shown in Figure 11 and Table 8, the computational cost of the DQ algorithm increases exponentially with the accuracy order: the execution time of the fifth-order model is 44.06 times that of the second-order model. Nevertheless, this cost is accompanied by a significant improvement in accuracy: compared with the second-order model, the fifth-order model reduces the position error by 84.04% and the velocity error by 66.55%. For real-time onboard processor applications, the feasibility of the fifth-order model depends on mission requirements. If high-precision tracking (e.g., maneuver reconstruction, collision avoidance) is prioritized, the fifth-order model is feasible provided that the mission observation interval is longer than the acceptable time cost of the fifth-order model. For scenarios where only rapid target localization is required, the second-order or third-order model represents a more balanced choice with lower overhead.

Notably, the order of the polynomial approximation model—a core design parameter for search-tracking algorithms and reachable domain modeling—directly influences this balance. Higher-order models introduce increased computational complexity but offer significant advantages in detection delay, tracking accuracy, and reachable domain calculation precision, thereby enhancing system convergence. In contrast, lower-order models sacrifice some accuracy to gain higher computational efficiency, making them more suitable for scenarios where only target localization is required without the need for subsequent maneuver reconstruction.

6. Discussion

The core academic contributions of this paper are outlined as follows: First, a hybrid PF-UKF framework is proposed to address the three key technical challenges in cislunar non-cooperative target tracking: data sparsity, strong nonlinear dynamics, and maneuver uncertainty. By synergizing the strengths of PF and UKF, the framework adopts a two-phase design that overcomes the inherent limitations of standalone filters in cislunar scenarios. Specifically, PF is exclusively utilized for reachable domain construction and global search during data-sparse periods, where it leverages undetected target information to dynamically shrink the search scope; UKF is activated only after target detection for high-precision, low-overhead tracking—given that UKF is inherently designed for tracking with observational data and lacks the capability to perform blind search in the absence of measurements, as its reliance on Gaussian distribution assumptions cannot accommodate non-Gaussian uncertainties from unknown target maneuvers. This phase-specific allocation resolves the long-standing trade-off between global search capability and computational efficiency that has plagued existing methods—overcoming PF’s “curse of dimensionality” in high-dimensional cislunar state spaces and UKF’s inherent lack of search capability without observational data. Simulation results validate this innovation: compared with the traditional PF + numerical integration method, the proposed fifth-order PF-UKF framework reduces time cost by 86.7% and position error by 24.1%, providing a feasible solution to the efficiency–accuracy dilemma in cislunar non-cooperative target tracking.

Second, the integration of the non-product-type DQ method breaks the computational bottleneck of high-dimensional reachable domain modeling, a critical complement to the two-phase framework. Instead of relying on sample-intensive traditional Monte Carlo sampling or low-efficiency tensor product quadrature, the DQ method constructs high-dimensional reachable domains with minimal quadrature points by matching the moment conditions of mixed distributions via Gauss–Newton optimization—ensuring both moment matching accuracy and computational efficiency. This design directly underpins the framework’s real-time performance, a prerequisite for time-sensitive cislunar missions. Practically, this capability aligns with the growing demand for SDA of non-cooperative targets such as space debris and unknown spacecraft: the framework is expected to provide fast and accurate tracking support for collision avoidance, sensor task scheduling, and mission safety, bearing potential value for safeguarding the safe operation of deep-space cislunar missions.

7. Conclusions

This work presents an effective framework for searching and tracking non-cooperative targets in cislunar space. Its main contributions are as follows: (1) An efficient reachable domain search and tracking framework based on DQ; (2) An intelligent switching mechanism between PF and UKF to optimize performance across different mission phases; (3) A flexible polynomial approximation approach for reachable domain modeling that balances accuracy and computational efficiency.

For future work, two research directions will be pursued: (1) Develop a Transformer-based approach for rapid reachable domain modeling, using attention mechanisms to capture nonlinear correlations between uncertainty parameters (e.g., maneuver time, velocity) and reachable domains, reducing computational complexity in high dimensions; (2) Introduce a multi-sensor federated filtering architecture fusing optical and radar data to enhance real-time detection and parameter estimation of target maneuvers (e.g., continuous pulse maneuvers), improving adaptability to complex cis-lunar scenarios.

Author Contributions

Conceptualization, K.L. and Y.W.; methodology, K.L.; validation, K.L., Y.W. and W.Z.; formal analysis, K.L.; investigation, K.L.; resources, Y.W.; data curation, K.L.; writing—original draft preparation, K.L.; writing—review and editing, K.L., Y.W. and W.Z.; visualization, K.L.; supervision, Y.W.; project administration, Y.W.; funding acquisition, Y.W. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data are unavailable due to privacy and ethical restrictions.

Acknowledgments

The authors would like to acknowledge the financial support provided by The National Natural Science Foundation of China. The financial support from the Natural Science Foundation for Distinguished Young Scholars of Hunan Province, China is also acknowledged.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1 Workflow of hybrid PF-UKF framework for target search and tracking.

Figure 2 Schematic of the uncooperative target search and tracking scenario in cislunar space.

Figure 3 Schematic diagram of the geometric relationship between platform sensor FOV and target. Note: The yellow arrow denotes the FOV center direction vector; the black arrow denotes the target direction vector; the dashed lines represent the boundary of the sensor FOV.

Figure 4 DQ algorithm flowchart.

Figure 5 Search and tracking scenario of Earth-Moon transfer orbit targets. (a) Observation schematic diagram. Note: the blue line represents the target trajectory, and the red line represents the platform trajectory. (b) Observation simulation scenario diagram.

Figure 6 Three-dimensional and planar (X-Y, X-Z, Y-Z) views of target reachable domain and platform FOV pointing at target discovery moment under second-order accuracy condition. Note: The blue cloud represents the target reachable domain, the yellow solid line denotes the platform’s pointing direction, the red dashed lines indicate the boundary of the FOV, the green rectangle represents the platform’s field of view, and the yellow triangle marks the platform position.

Figure 7 Search and tracking results of non-Cooperative targets in search and tracking scenarios under sceond-order accuracy condition. Note: In the figure, the blue curve represents the position/velocity error variation; the orange dashed line denotes the state corresponding to undetected target; the orange circles represents the “target undetected” state; and the blue circles represents the “target detected” state.

Figure 8 Three-dimensional and planar (X-Y, X-Z, Y-Z) views of target reachable domain and platform FOV pointing at target discovery moment under fifth-order accuracy condition.

Figure 9 Search and tracking results of non-Cooperative targets in search and tracking scenarios under fifth-order accuracy condition.

Figure 10 Search and tracking results of non-Cooperative targets in search and tracking scenarios. (a) The search and tracking performance of the traditional PF + numerical integration method. (b) The search and tracking performance of the PF + DQ method.

Figure 11 Time cost variation of DQ algorithm with different accuracy orders.

Table 1

Observation platform parameter settings.

Parameter	Value
Field of view angle	3°
Observation distance (detection probability attenuation scale)	0.6 L (L: average distance of Earth–Moon space)
Angle measurement accuracy	2 arcsec

Table 2

Parameter settings of simulation scenarios.

Parameter	Value
Simulation duration	0–24 h
Observation interval	10 min
Initial number of particles	$1\times {10}^4$

Table 3

Parameter settings of target maneuvering.

Parameter	Value
Maximum maneuvering capability	100 m/s
Maneuvering time range	[0, 24 h]

Table 4

Initial state parameter settings of target and platform in Earth-Moon centroid rotating coordinate System.

State Vector Under EMI	3D DRO Orbit Observation Platform	Earth–Moon Transfer Orbit Target
x [km]	303,649.49	263,948.22
y [km]	26.16	130,460.52
z [km]	−15,375.99	3673.20
vx [km/s]	0	0.91
vy [km/s]	0.54	−0.25
vz [km/s]	0.02	0.01

Table 5

Tracking error under second-order accuracy condition.

Parameter	Value
Position error [km]	$2.346\times {10}^0$
Velocity error [m/s]	$2.904\times {10}^{-2}$
System Response Time [T *]	$4.00\times {10}^{-2}$

* Here, $T={\displaystyle \frac{1}{\sqrt{\frac{G·(m_E+m_M)}{L^3}}}}$ denotes the characteristic time scale, where G is the gravitational constant, L represents the average Earth–Moon distance, and the corresponding value is $T=104.2194$ h.

Table 6

Tracking error under fifth-order accuracy condition.

Parameter	Value
Position error [km]	$3.745\times {10}^{-1}$
Velocity error [m/s]	$9.703\times {10}^{-3}$
System Response Time [T]	$1.76\times {10}^{-2}$

Table 7

Performance comparison of different methods.

Method	Time Cost [s]	Position Error [m]	Velocity Error [m/s]
PF + Numerical Integration	$9.013\times {10}^2$	$4.931\times {10}^2$	$1.640\times {10}^{-2}$
PF + DQ	$2.121\times {10}^2$	$3.306\times {10}^5$	$3.660\times {10}^0$
PF-UKF (Second-order)	$3.978\times {10}^1$	$2.346\times {10}^3$	$2.904\times {10}^{-2}$
PF-UKF (Fifth-order)	$1.201\times {10}^2$	$3.745\times {10}^2$	$9.703\times {10}^{-3}$

Table 8

Computational overhead comparison of DQ algorithm under different accuracy orders.

Accuracy Order	Time Cost [s]
Second-order accuracy	$1.798\times {10}^0$
Third-order accuracy	$5.008\times {10}^0$
Fourth-order accuracy	$2.328\times {10}^1$
Fifth-order accuracy	$7.920\times {10}^1$

References

1. Johnson, S.K.; Mortensen, D.J.; Chavez, M.A.; Woodland, C.L. Gateway—A communications platform for lunar exploration. IET Conf. Proc.; 2021; 2021, pp. 9-16. [DOI: https://dx.doi.org/10.1049/icp.2022.0544]

2. Smith, M.; Craig, D.; Herrmann, N.; Mahoney, E.; Krezel, J.; McIntyre, N.; Goodliff, K. The Artemis Program: An Overview of NASA’s Activities to Return Humans to the Moon. Proceedings of the 2020 IEEE Aerospace Conference; Big Sky, MT, USA, 7–14 March 2020; IEEE: New York, NY, USA, 2020; pp. 1-10. [DOI: https://dx.doi.org/10.1109/aero47225.2020.9172323]

3. Adurthi, N.; Singla, P.; Majji, M. Mutual Information Based Sensor Tasking with Applications to Space Situational Awareness. J. Guid. Control Dyn.; 2020; 43, pp. 767-789. [DOI: https://dx.doi.org/10.2514/1.G004399]

4. Evans, M.E.; Graham, L.D. A Flexible Lunar Architecture for Exploration (FLARE) supporting NASA’s Artemis Program. Acta Astronaut.; 2020; 177, pp. 351-372. [DOI: https://dx.doi.org/10.1016/j.actaastro.2020.07.032] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32834186]

5. Hou, X.; Li, B.; Liu, X.; Cheng, H.; Shen, M.; Wang, P.; Xin, X. Initial orbit determination of some cislunar orbits based on short-arc optical observations. Astrodynamics; 2024; 8, pp. 455-469. [DOI: https://dx.doi.org/10.1007/s42064-024-0210-z]

6. Black, A.; Frueh, C. Fragmentation characterization in the circular restricted three body problem for cislunar space domain awareness. Adv. Space Res.; 2025; 75, pp. 1177-1204. [DOI: https://dx.doi.org/10.1016/j.asr.2024.08.076]

7. Kaderali, S.; Misra, A. Detection and characterisation of unknown maneuvers in spacecraft. Acta Astronaut.; 2023; 205, pp. 310-318. [DOI: https://dx.doi.org/10.1016/j.actaastro.2023.01.018]

8. Zaborsky, S. Generating Solutions for Periodic Orbits in the Circular Restricted Three-Body Problem. J. Astronaut. Sci.; 2020; 67, pp. 1300-1319. [DOI: https://dx.doi.org/10.1007/s40295-020-00222-3]

9. György, K.; Kelemen, A.; Dávid, L. Unscented Kalman Filters and Particle Filter Methods for Nonlinear State Estimation. Procedia Technol.; 2014; 12, pp. 65-74. [DOI: https://dx.doi.org/10.1016/j.protcy.2013.12.457]

10. Li, Z.; Wang, Y.; Zheng, W. Adaptive Consensus-Based Unscented Information Filter for Distributed Space Target Tracking. Proceedings of the 2021 5th Chinese Conference on Swarm Intelligence and Cooperative Control; Lecture Notes in Electrical Engineering; Springer: Berlin/Heidelberg, Germany, 2023; pp. 964-974.

11. Hall, Z.; Singla, P.; Johnson, K. Reachability-Based Search for Tracking of Noncooperative Maneuvering Satellites in Data Sparse Environment. J. Astronaut. Sci.; 2023; 70, 9. [DOI: https://dx.doi.org/10.1007/s40295-023-00365-z]

12. Wang, Y.; Zheng, W.; Sun, S.; Li, L. Robust Information Filter Based on Maximum Correntropy Criterion. J. Guid. Control Dyn.; 2016; 39, 1124. [DOI: https://dx.doi.org/10.2514/1.G001576]

13. Ye, X.; Wang, J.; Wu, D.; Zhang, Y.; Li, B. A Novel Adaptive Robust Cubature Kalman Filter for Maneuvering Target Tracking with Model Uncertainty and Abnormal Measurement Noises. Sensors; 2023; 23, 6966. [DOI: https://dx.doi.org/10.3390/s23156966] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37571748]

14. Han, B.; Ge, Z.; Su, Z.; Hao, J. Research on a Particle Filtering Multi-Target Tracking Algorithm for Distributed Systems. Sensors; 2025; 25, 3495. [DOI: https://dx.doi.org/10.3390/s25113495] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/40969049]

15. Holzinger, M.; Scheeres, D. Applied Reachability for Space Situational Awareness and Safety inSpacecraft Proximity Operations. Proceedings of the AIAA Proceedings; Chicago, IL, USA, 10–13 August 2009.

16. Xue, D.; Li, J.; Baoyin, H.; Jiang, F. Reachable Domain for Spacecraft with a Single Impulse. J. Guid. Control Dyn.; 2010; 33, pp. 934-942. [DOI: https://dx.doi.org/10.2514/1.43963]

17. Wen, C.X.; Zhao, Y.S.; Shi, P. Precise determination of reachable domain for spacecraft with single impulse. J. Guid. Control Dyn.; 2014; 37, pp. 1767-1779. [DOI: https://dx.doi.org/10.2514/1.G000583]

18. Hall, Z.; Singla, P. Higher-order sensitivity matrix method for probabilistic solution to uncertain Lambert problem and reachability set problem. Celest. Mech. Dyn. Astron.; 2020; 132, 50. [DOI: https://dx.doi.org/10.1007/s10569-020-09988-y]

19. Wang, J.; Rui, X.; He, B.; Wang, X.; Xie, K. A Fast Calculation Method for Sensitivity Analysis Based on Multibody System Transfer Matrix Method. Int. J. Appl. Mech.; 2025; 17, 2550008. [DOI: https://dx.doi.org/10.1142/S1758825125500085]

20. Hu, L.; Ma, L.; Shen, J. Sparse Quadrature for High-Dimensional Integration with Gaussian Measure. ESAIM Math. Model. Numer. Anal.; 2018; 52, pp. 631-657.

21. Jones, R.; Curtis, D.; Zagaris, C. Reachable Set Approximation in Cislunar Space with Pseudospectral Method and Homotopy. Proceedings of the 2023 IEEE Aerospace Conference; Big Sky, MT, USA, 4–11 March 2023.

22. Li, Z.; Wang, Y.; Zheng, W.; Song, X. Sensor scheduling with orbital uncertainty for space target tracking. Acta Astronaut.; 2023; 205, pp. 23-32. [DOI: https://dx.doi.org/10.1016/j.actaastro.2023.01.030]

23. Wolf, T.; Zucchelli, E.; Jones, B.A. Multi-Fidelity Uncertainty Propagation for Objects in Cislunar Space. Proceedings of the AIAA SciTech Forum and Exposition; San Diego, CA, USA, 3–7 January 2022.

24. Tang, R.; Chai, R.; Li, P. Deep Reinforcement Learning-based Sensing and Communication Scheduling Algorithm for UAV-Assisted Target Detection Systems. Proceedings of the 2023 IEEE 98th Vehicular Technology Conference (VTC2023-Fall); Hong Kong, China, 10–13 October 2023.

25. Qin, J.; Bai, X.; Du, G.; Liu, J.; Peng, N.; Xu, M. Multi-Satellite Scheduling for Moving Targets Using Enhanced Hybrid Genetic Simulated Annealing Algorithm and Observation Strip Selection. IEEE Trans. Aerosp. Electron. Syst.; 2024; 60, pp. 5773-5800. [DOI: https://dx.doi.org/10.1109/TAES.2024.3397958]

26. Hall, Z.; Schwab, D.; Eapen, R.; Singla, P. Reachability-Based Approach for Search and Detection of Maneuvering Cislunar Objects. Proceedings of the AIAA SCITECH 2022 Forum; San Diego, CA, USA, 3–7 January 2022; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2022; [DOI: https://dx.doi.org/10.2514/6.2022-0853]

27. Gautschi, W. Orthogonal Polynomials: Computation and Approximation; Oxford University Press: New York, NY, USA, 2004; [DOI: https://dx.doi.org/10.1093/oso/9780198506720.001.0001]

28. Keshavarzzadeh, V.; Kirby, R.M.; Narayan, A. Numerical Integration in Multiple Dimensions with Designed Quadrature. SIAM J. Sci. Comput.; 2018; 40, pp. A2033-A2061. [DOI: https://dx.doi.org/10.1137/17M1137875]