High accuracy ranging for space debris with spaceborne single photon Lidar

Yuan Tian; Yuan Tian; Yuan Tian; Yuan Tian; Xiaodong Hu; Songmao Chen; Songmao Chen; Songmao Chen; Yixin Zhao; Yixin Zhao; Yixin Zhao; Yixin Zhao; Xuan Zhang; Xuan Zhang; Xuan Zhang; Xuan Zhang; Dingjie Wang; Dingjie Wang; Dingjie Wang; Dingjie Wang; Weihao Xu; Weihao Xu; Weihao Xu; Weihao Xu; Meilin Xie; Meilin Xie; Meilin Xie; Wei Hao; Wei Hao; Wei Hao; Xiuqin Su; Xiuqin Su; Xiuqin Su

doi:10.1364/OE.519002

1. Introduction

The escalating challenge of space debris necessitates robust mitigation strategies, with current tracking technologies reliant on ground-based sensors limited by environmental and locational constraints [1,2]. In contrast, spaceborne systems, offering greater maneuverability and coverage, emerge as a promising solution, enhancing the efficiency and safety of space activities through advanced techniques like laser ranging [3–5]. High-sensitivity detectors such as single photon avalanche diodes (SPADs) are gaining significance in detection for non-cooperative dynamic targets, including space debris [6]. Time-correlated single photon counting (TCSPC)-based Lidar cumulates time-based signals to construct discrete reflected waveforms, effectively compensating for intensity response loss attributable to the single-photon sensitivity of SPADs. This technique efficiently gathers target distance and reflectivity data through repeated pulse sampling, extracting significant information from limited reflected photons without necessitating high signal-to-noise ratios (SNR) or laser power [7]. It is evident that dynamic targets, especially highly dynamic targets, pose significant challenges to the cumulation process. The cumulative distribution of reflected photons deviates from the Gaussian or quasi-Gaussian distribution typically observed under static conditions, evolving into a more complex pattern. This leads to suboptimal performance of traditional histogram-based peak extraction methods [8].

Current research in single photon Lidar (SPL) ranging for dynamic targets generally falls into two main categories. The first involves generating point clouds via single-shot ranging and deriving distance trajectories using algorithms like mean-shift [9], RANSAC [10], and the Hough transform [11]. Bao et al. employed coincidence photon-counting, which improved the SNR [12]. Du et al. developed a photon-counting laser ranging system for object tracking, combining mean-shift and RANSAC algorithms to accurately extract target distance trajectory [13]. Similarly, Xue et al. and Liu et al. proposed a feature extraction method based on Hough Transform to ascertain distance [14,15]. Peng et al. introduced dynamic coincidence which partially addressed the problem by matching detection windows [16]. Overall, these methodologies encounter obstacles in achieving absolute ranging accuracy, which are attributed to the point cloud generation process under this scheme not distinguish between effective signals and noise. The second devides the photon acquisition into multiple sub-intervals and gathers photon cumulation histograms for each through multiple measurements. Subsequent statistical analysis and processing correct waveform distortion and broadening caused by target motion and facilitate the estimation of distance. Jonsson et al. divided measurements into time intervals and cross-correlated histograms to estimate target movement parameters [17,18]. Yu et al. enhanced the capability to acquire distance via time-shift cumulation with a multi-resolution time grid [8]. Hou et al. proposed a structural acquisition method which combined with waveform correction and Gaussian fitting to solve the distoration problem [19]. In essence, these methodologies, generally suitable for linear motion, lack adequate robustness for targets characterized by nonlinear and high-dynamic motion, leading to restricted ranging accuracy and potentially sacrificing measurement resolution. Our preceding research introduced an dynamic impulse response function (IRF) specifically designed for dynamic targets, which significantly improved the ranging accuracy. Nonetheless, the study did not address targets characterized by nonlinear motion [20].

In this study, anchored in non-parametric Bayesian inference [21], a distance extraction method based on DPMM [22,23] has been proposed to address the complexities arising from the nonlinear and high-dynamic motion characteristics of LEO space debris in SSPL ranging. By establishing state transition and observation models for space debris and employing Bayesian filtering exemplified by EKF, the approach refines range gating guidance. This integration enhances space debris detection by improving distance predictions, thereby boosting detection efficiency. The probability density function (PDF) of the temporal distribution of reflected photons from space debris was characterized as a mixture distribution controlled by several prior distributions, controlling the mixture weights through the Dirichlet process and employing MCMC [24] for posteriori estimation. Experimental results demonstrate that this method, particularly effective for typical LEO space debris, not only provides more accurate distance than traditional methods based on point cloud feature extraction and matched filtering, but also exhibits a correlation between ranging error and target speed that aligns with intuitive expectations. The remainder of the paper is organized as follows. In Section 2, the composition and on-orbit detection mode of the SSPL are introduced. In Section 3, the theoretical model is introduced, and a novel distance extraction method based on DPMM is proposed. In Section 4, simulation scenarios are established based on actual conditions, and simulation validation for selected debris was carried out. Concluding remarks are drawn in Section 5.

2. Structure of proposed system

The schematic diagram of the on-orbit operation of the SSPL for space debris, as involved in this study, is shown in Fig. 1(a) and the composition of the system is illustrated in Fig. 2. The system predominantly comprises the telescope system, tracking system, laser, signal generator, TCSPC, and SPAD. In traditional ground-based detection systems, the orbital state of the target debris is estimated using observational data from ground observation stations and forecast information [25] about the debris. Initial orbit data from the target debris’s orbital elements is processed with suitable dynamical models to calculate the tracking ephemeris, which is then converted into position data within the observation satellite’s coordinate system. In this study, a state transition model for the target debris is constructed based on orbital dynamics, with spatial geometric model forming the basis of the observational model. The state prior of the target derived through EKF is employed to set the range gating for space debris. During its on-orbit operation, once the tracking system successfully maintains stable tracking of the target debris, the SSPL initiates ranging, and comprehensive data analysis and meticulous processing procedures are employed to improve the ranging accuracy and further refine the forecast information based on the state prior.

Fig. 1. The on-orbit operation and spatial geometric model schematic diagram for the proposed SSPL system. (a) LEO on-orbit operation. (b) Spatial geometric model.

Download Full Size | PDF

Fig. 2. System architecture of the proposed SSPL system.

Download Full Size | PDF

3. Model and proposed method

3.1 Model

In this study, a spatial geometric model has been established to characterize the relative positioning of space debris and the SSPL, as demonstrated in Fig. 1(b). Comprehensive analysis of noise sources within the detection scenario has also been conducted. Due to forecasting inaccuracies and the debris’s high dynamics, errors in tracking and pointing arise, affecting the reflected photon count from the debris and thus the detection probabilities. For a detailed analysis, refer to Appendix A.

The range gating employed in this study is guided by the state prior of the target debris obtained with the EKF, with details on the gate width and EKF parameters settings provided in Section 4. The state transition model is constructed based on the orbital dynamics of space debris, with the detailed process outlined in Appendix B.

3.2 Proposed method

The comparative analysis of Fig. 3(a) and (b) reveals that the radial motion of the target relative to the SSPL restricts the cumulation process, resulting in distortion of the histogram and a reduced SNR. Traditional methods, which aggregate reflected photons into cumulative histograms for peak detection, experience significant challenges in addressing this issue. The nonlinear and dynamic orbital motion of space debris results in suboptimal performance when applying a uniform system instrument response function (IRF) for photon distribution characterization [20]. Consequently, conventional matched filtering techniques fail to provide accurate distance measurements for targets with complex motion profiles, as reflected photons exhibit multi-peaked and broadened distributions. By leveraging the DPMM [26,27], this study presents a novel approach to accurately model the dynamic accumulation of reflected photons and address the multimode characteristics of their temporal distribution. The DPMM’s flexibility in modeling the complex, nonlinear motion of space debris significantly enhances the accuracy of distance measurements, overcoming the limitations of traditional methods. The integration of the DPMM with EKF and iterative refinement through MCMC allow for effective handling of measurement uncertainty and noise, while dynamically updating the model to accurately reflect the complex distribution of reflected photons, leading to more accurate ranging of space debris.

Fig. 3. Comparison of the cumulation distributions of reflected photons from different targets under the same background noise conditions. (a) The cumulation distribution of reflected photons from a static target. (b) The cumulation distribution of reflected photons from a dynamic target.

Download Full Size | PDF

3.2.1 Preprocessing

The temporal attributes of reflected photons of non-cooperative targets in space usually have nonlinear and non-stationary characteristics coupled with a generally low SNR. Therefore, to effectively address these complexities, a mode decomposition approach that combines empirical mode decomposition (EMD) [28] and variational mode decomposition (VMD) [29] is employed here for the preliminary denoising of the signal. These techniques are crucial for isolating and mitigating noise components, thereby enhancing the clarity and robustness of the signal analysis. The return signal $s\left ( t \right )$ is initially decomposed into multiple IMFs utilizing EMD, which is articulated as

(1)$$s\left( t \right) = \sum_{i = 1}^N {{{\cal C}_i}\left( t \right) + r\left( t \right)}$$

where ${{\cal C}_i}\left ( t \right )$ denotes the $i$th mode component, $i = 1,2,\ldots,N$, and $r\left ( t \right )$ is the residual. For each ${{\cal C}_i}\left ( t \right )$, it is targeted to have a specific narrow bandwidth around a center frequency ${\omega _{ij}}$, and can be further decomposed into multiple IMFs via VMD, which is articulated as

(2)$${{\cal C}_i}\left( t \right) = \sum_{j = 1}^M {{\hat \nu_{ij}}\left( t \right)}$$

where ${\hat \nu _{ij}}\left ( t \right )$ denotes the $j$th mode component, $j = 1,2,\ldots,M$. The variational problem associated with this procedure can be formulated as

(3)$$\mathop {\arg \min }_{{\hat \nu_{ij}},{\hat{\omega}_{ij}}} \sum_{j = 1}^M {\int {{{\left( {{\partial _t}\left( {\frac{{{\mathop{\rm Re}\nolimits} \left( {{\hat \nu_{ij}}\left( t \right){e^{ - j{\hat{\omega}_{ij}}t}}} \right)}}{{1 + j{\hat{\omega}_{ij}}t}}} \right)} \right)}^2}} } dt,{\rm{ s}}{\rm{.t}}{\rm{. }}{\mathcal{C}_{i}}\left( t \right) = \sum_{j = 1}^M {{\mathop{\rm Re}\nolimits} \left( {{e^{j{\hat{\omega}_{ij}}t}}{\hat \nu_{ij}}} \right)}$$

where $M$ denotes the number of IMFs decomposed using VMD. The output of Eq. (3) constitute the solutions to the optimization problem, specifically ${\hat \nu _{ij}}\left ( t \right )$ and ${\hat \omega _{ij}}$. This optimization problem is solved using the alternating direction method of multipliers (ADMM) [30]. Following the extraction of mode components, an assessment is conducted to identify noise-impacted modes. This assessment involves an analysis of each mode’s energy, frequency spectrum, and statistical characteristics. A threshold criterion is established to differentiate between noise and signal-bearing modes. Modes falling below this threshold are deemed to be noise-dominated and are thus excluded from signal reconstruction. Reconstruction of the signal is undertaken using only those mode components identified as carriers of authentic signal data. Designating $I$ as the collection of indices for the selected modes, the formula for the reconstructed signal is thus expressed as

(4)$$\hat x\left( t \right) = \sum_{i \in I} {{\nu_{ij}}\left( t \right)}$$

For target debris, a cumulative unit within the detection window was selected to compare the denoised signals with the original signals, as depicted in Fig. 4. Fluctuations in the average signal-to-background ratio (SBR) of reflected signals within the detection window, as well as variations in the average SBR of denoised signals, are illustrated in Fig. 5. The average SBR of the original signals within the detection window were recorded as 4.82, 4.82, 4.61 and 4.68, respectively. Following the denoising process, the average SBRs of the denoised signals improved to 29.69, 28.84, 27.56 and 27.30, respectively. Analysis of Fig. 4 and Fig. 5 reveals that the proposed method is capable of effectively reducing noise in reflected signal representations, thereby enhancing the SBR. Combining EMD and VMD leverages the advantages of both methods to achieve a more accurate time-frequency representation and more effective denoising results. EMD offers an adaptive approach to capture the primary characteristics of the signal, while VMD further separates and refines these modes, providing a more precise foundation for denoising.

Fig. 4. The original and denoised signals with proposed method. (a)(c)(e)(g) Original signal of target debris 49607, 50078, 50578, 28742. (b)(d)(f)(h) Denoised signal of target debris 49607, 50078, 50578, 28742.

Download Full Size | PDF

Fig. 5. The average SBR within the detection window of target debris. (a) Target debris 49607. (b) Target debris 50078. (c) Target debris 50578. (d) Target debris 28742.

Download Full Size | PDF

3.2.2 DPMM

Given a sample space $\mathcal {L}$ and a base measure $G_0$, let Dir(.) denotes the Dirichlet PDF and ${L_{1:K}}$ be $K$ disjoint subsets of $\mathcal {L}$ satisfying $L_1 \cup L_2 \cup \ldots \cup L_k = \mathcal {L}$. The DP with a concentration parameter $\alpha _c$ is defined as

(5)$$(G(L_1), G(L_2), \ldots, G(L_k)) \sim \text{Dir}(\alpha_c G_0(L_1), \alpha_c G_0(L_2), \ldots, \alpha_c G_0(L_k))$$

which implies that the random probability distribution $G$ is sampled from the DP. The joint distribution of $N$ latent variables ${\epsilon }_{1:N}$ sampled from $\mathcal {L}$ with the distribution $G$ can be expressed as

(6)$$p({\epsilon}_{1:N}|\alpha_c, G_0) = \int p({\epsilon}_{1:N}|G) p(G|\alpha_c, G_0) dG$$

which demonstrates a clustering effect [27]. Gibbs sampling [26,31], which is one type of MCMC sampling [32,33], is carried out for Eq. (6) to obtain the discrete joint distribution of ${\epsilon }_{1:N}$. Maximizing the posterior estimation on this joint distribution produces the most probable values of ${\epsilon }_{1:N}$, which are used to partition ${\epsilon }_{1:N}$ [26]. The conditional probability distribution used to express the Gibbs sampling can be expressed as [34]

(7)$$p(\epsilon_i | \epsilon_{1:i-1}) \propto \alpha_c G_0(\epsilon_i) + \sum_{j=1}^{i-1} \delta(\epsilon_i - \epsilon_j)$$

where $\delta (\varphi )$ is the Dirac delta function. Eq. (7) suggests that given the preceding $i-1$ values of $\epsilon _{i-1}$, the value of $\epsilon _i$ is contingent upon either an extant value from $\epsilon _{i-1}$ or a novel value derived from the base measure $G_0$. The introduction of a novel value, signifying the inception of a new cluster, facilitates the autonomous determination of the requisite number of clusters. Eq. (7) can be further expressed as

(8)$$p(\epsilon_i | \epsilon_{{-}i}) \propto \alpha_c G_0(\eta_i) + \sum_{j \in \mathcal-i} \delta(\epsilon_i - \epsilon_j)$$

if the random variables ${\epsilon }_{1:N}$ are exchangeable, where $-i$ denotes the indices in $1:N$ except for $i$. Here the DP is employed as a non-parametric piror and the DPMM can be expressed as

(9)$$\begin{aligned} & \tau_i \sim p({\cdot}|\epsilon_i) \\ & \epsilon_i \sim G \\ & G \sim \text{DP}(\alpha_c, G_0) \end{aligned}$$

where $\tau _i$ is the $i$th observation. Each observation $\tau _i$ has its own latent variable $\epsilon _i$, correspond to the cluster label for $\tau _i$. In other words, the Gibbs sampling employed here is to determine the values of the latent variables $\{\epsilon _i\}$ given the observation$\{\tau _i\}$.

3.2.3 DPMM based distance extraction

For space debris and SSPL, given observed data $\mathcal {T} = \left \{ {{\tau _1},{\tau _2},\ldots,{\tau _N}} \right \}$, the posterior probability $p(\epsilon _{1:N} | \tau _{1:N})$ is approximated by sampling from $p(\epsilon _i | \epsilon _{1:i-1},\tau _{1:N})$ iteratively with the Gibbs sampler. The probability $p(\epsilon _i | \epsilon _{-i}, \tau _{1:N})$ is represented as

(10)$$p(\epsilon_i | \epsilon_{{-}i}, \tau_{1:N}) \propto f(\tau_i | \epsilon_i) p(\epsilon_i | \epsilon_{{-}i})$$

according to Bayes rule. Where $\tau _i$ denotes the time of flight (ToF) of the $i$th reflected photon, $f(\tau _i | \epsilon _i)$ denotes the likelihood, which can be conceptualized with a mixture distribution comprising multiple Gaussian or quasi-Gaussian and an uniform distribution, formulated as

(11)$$f(\tau_i | \epsilon_i) = \sum_{k = 1}^{K,K \le {K^*}} {{\pi _k}{\cal N}} \left( {{\epsilon_i}\left| {{\mu _k},{\sigma ^2}} \right.} \right) + {\pi _{K + 1}}{\cal U}\left( {{a,b}} \right)$$

and the selection of the base measure ${G_0}$, as expressed by the mixture of Gaussian and uniform distributions

(12)$$G_0 = p\mathcal{N}(X | \mu_0, \sigma^2) + (1 - p)\mathcal{U}(X | a, b)$$

is congruent with the data’s inherent characteristics and integrate prior knowledge concerning the parameters. Where ${\cal N}$ represnts the Gaussian distribution, ${\cal U}$ represnts the uniform distribution and [a, b] is the range gating set based on the EKF illustrated in Appendix B. $p$ is the prior probability of the data coming from a Gaussian distribution, and $\mu _0$ is the prior mean of the Gaussian distribution. The means $\mu _k$ of the Gaussian components in the mixture model are to be estimated, yet it is assumed here that they all exhibit a common variance $\sigma ^2$, which is equated to the pulse width of the return signal under ideal conditions. $K$ denotes the quantity of clusters and $K^*$ is introduced here as the supremum of $K$, which is determined by the relative radial velocity of the target debris and the resolution of TCSPC, and can be expressed as

(13)$$K^* = \sup \{K \mid K \in \mathbb{K}\} = \left\lceil \frac{(b - a) \cdot \bar{\nu}}{(c \cdot \Delta \tau / 2)} \right\rceil$$

where $\mathbb {K}$ denotes the set of cluster counts, $\bar {\nu }$ denotes the maximum mean velocity within different intervals, $\Delta \tau$ denotes the resolution of TCSPC and $c$ is the speed of light. The DP is employed as a piror for the mixture model, where each $\epsilon _i$ correndsponds to a potential cluster center $\mu _k$. A practical representation of the DP is the stick-breaking construction [35], for the mixing proportions $\pi _k$, we have

(14)$${\pi _k} = {\beta _k}\prod_{j = 1}^{k - 1} {\left( {1 - {\beta _j}} \right)} ,{\rm{ }}{\beta _k} \sim {Beta} \left( {1,\alpha_c } \right)$$

where $\beta _k$ are independent Beta-distributed random variables. In Gibbs sampling, the conditional distributions

(15)$$p(\epsilon_i = k | \epsilon_{{-}i}, \tau_i, \mu_k, \sigma^2) \propto (n_{{-}i,k} + \frac{\alpha_c}{K}) \mathcal{N}(\tau_i | \mu_k, \sigma^2)$$

are sampled for existing clusters, and

(16)$$p(\epsilon_i = K+1 | \epsilon_{{-}i}, \tau_i, \mu, \sigma^2) \propto \alpha_c \mathcal{U}(\tau_i | a, b)$$

for new clusters to infer the posterior distribution of $\epsilon$. $n_{-i,k}$ is the number of photons in cluster $k$ excluding $\tau _i$. The updates for the cluster centers $\mu _k$ and the mixing proportions $\pi _k$ can be performed given the cluster assignments as

(17)$$\mu_k | \tau, \epsilon, \sigma^2 \sim \mathcal{N}\left(\frac{\sum_{i : \epsilon_i = k} \tau_i}{n_k}, \frac{\sigma^2}{n_k}\right)$$

(18)$$\pi_k | \epsilon, \alpha_c \sim {Beta}(1 + n_k, \alpha_c + \sum_{j=k+1}^{K} n_j)$$

where $n_k$ is the number of photons in cluster $k$. The estimated distance for each cumulative unit is equal to the the means of these cluster centers, which can be expressed as

(19)$${\hat d} = \frac{1}{2K}\sum_{k = 1}^K {{\mu_k}\cdot c}$$

Figure 6 illustrates the flowchart of the proposed distance extraction based on DPMM. Wherein, ${N_{iter}}$ denotes the number of iterations for the MCMC. By iterating the above steps with Gibbs sampling, the DPMM can adaptively identify the number of Gaussian distributions present in the reflected photons, as well as their respective parameters. Ultimately, this yields the distribution of clusters, alongside estimations of each cluster’s mean and mixing proportions. The DPMM approach effectively acquires structural information from reflected photons with undetermined clusters, leading to accurate distance extraction of target debris.

Fig. 6. The flowchart of the proposed DPMM based distance extraction.

Download Full Size | PDF

4. Experiment results

4.1 System parameters

The parameter configuration for proposed system is delineated in accordance with the details presented in Table 1. The ranging system employed a laser operating at a wavelength of 1064 nm, with 10 ns duration pulses (full width at half maximum, FWHM) at a repetition of 1 MHz and the mean power of the laser is 5 W. The transmitting divergence angle is designed to be 35 $\mu$rad. Accounting for the filter’s transmittance and the lens assembly’s transmission in the receiving optical system, the overall transmission efficiency of the transmitting optical system is 0.8. The transmission efficiency of the receiving optical system is also 0.8. Additionally, the receiving aperture diameter is established at 1 m, while the receiving field of view (FoV) is designed to be 20 urad. The time jitter of the system is 2 ns, with a dark count rate of 300 cps (counts per second), and a detector dead time of 20 ns. The filter has a bandwidth of 1 nm@1064 nm. The TCSPC offers a time resolution of 1 ns, corresponding to a distance resolution of 0.15 m. The target debris assumes the form of a regular rectangular piece measuring 1 m $\times$ 1 m. The BRDF parameters are adopted from the values specified in [36].

Table 1. System parameters

View Table | View all tables in this article

4.2 Scenario

Following the completion of system model development and parameter configuration, the next step involves setting up the detection scenario. The SSPL was envisioned as being mounted aboard the International Space Station (ISS, NORAD CAT ID: 25544). The simulation spanned from 04:00:00 on May 11, 2023, to 04:00:00 on May 12, 2023. Within a 24-hour period, four space debris on different orbits near their respective time of closest approach (TCA) [5] were tracked with a 10s interval to exploit their nonlinear motion, particularly in distance, showcasing the proposed method’s accuracy. It was hypothesized that during the detection interval, the visibility of the debris would not be hindered by Earth’s shadow, and conditions of ambient light from the ground, daylight, and moonlight were conducive to observation. As shown in Table 2, the detection windows for the target debris is presented. These debris were represented for convenience by their respective catalog numbers: COSMOS_1408_DEB (49607, 50078, 50578) and CZ-2D_DEB (28742).

Table 2. The detection windows for target debris during the observation period.

View Table | View all tables in this article

Fig. 7(a)-(d) and Fig. 8(a)-(d) respectively depict the position (indicated with azimuth-elevation-range, AER) and rate of position change (indicated with azimuth-elevation-range rate, AER rate) of the target debris within their respective detection windows. It is evident that the motion of the target debris within these windows demonstrates nonlinear and high-dynamic motion characteristics, especially in the distance dimension. To harmonize computational efficiency with data accuracy, the study uses data sampled at 1ms intervals, subjected to high-order polynomial interpolation to meet the strict temporal resolution requirements of high-repetition-rate SSPL in space debris detection, ensuring accurate data processing for tracking and ranging.

Fig. 7. The AER of target debris. (a) The AER of 49607. (b) The AER of 50078. (c) The AER of 50578. (d) The AER of 28742.

Download Full Size | PDF

Fig. 8. The AER rate of target debris. (a) The AER rate of 49607. (b) The AER rate of 50078. (c) The AER rate of 50578. (d) The AER rate of 28742.

Download Full Size | PDF

The analytical assessment of the detection system’s tracking and pointing proficiency posits that the PDFs for both tracking and pointing inaccuracies along the azimuth and elevation axes are characterized by Gaussian distribution, with a mean value of zero and a standard deviation $\sigma _t$ set to $0.2''$, denoted as: ${\delta _{\alpha _p}} \sim \mathcal {N}(0,\sigma _t^2)$, ${\delta _{\beta _p}} \sim \mathcal {N}(0,\sigma _t^2)$, ${\theta _{\alpha _p}} \sim \mathcal {N}(0,\sigma _t^2)$, and ${\theta _{\beta _p}} \sim \mathcal {N}(0,\sigma _t^2)$. This presumption enables accurate statistical modeling of the system’s performance, particularly in regards to the fidelity of angular measurements.

Within the range gating guided by EKF, the fixed prior value $d_r$ is set to 50 m, ensuring that at least 90% of the reflected photons within the detection window are captured. The initial probabilistic model for the distance factor of the observation vector adheres to a uniform distribution across the interval of (20, 50) m, reflecting the initial uncertainty in ranging for the state estimation process. This statistical parameterization is integral for the subsequent update phases within the filter algorithm. The stochastic properties of process noise are characterized by a Gaussian PDF as delineated by ${\boldsymbol{w}_k} \sim \mathcal {N}\left (0, \boldsymbol{Q}_k\right )$, with the matrix $\boldsymbol{Q}_k = I(n)$ denotes the process noise covariance matrix (PNCM), where $I$ is the identity matrix. In parallel, the observational noise adheres to a Gaussian distribution as conveyed by ${\boldsymbol{v}_k} \sim \mathcal {N}\left (0, \boldsymbol{R}_k\right )$, with the matrix $\boldsymbol{R}_k = [\sigma _t^2, \sigma _t^2, \Delta d]^T$ denotes the measurement noise covariance matrix (MNCM), where $\Delta d$ quantifying the distance error attributable to detector’s timing jitter. For a detailed derivation and analysis of range gating guided by EKF, refer to Appendix B.

4.3 Experiment results

In this study, four LEO space debris were selected as target debris to evaluate the performance of the SSPL system proposed. It employed a comprehensive set of simulated data for ranging, employing the Monte Carlo simulation [37], considering the system configuration and detection scenarios outlined in Sections 4.1 and 4.2. The study compared the results of the proposed method with established point cloud feature extraction based statistical methods like density filtering [38], the Hough Transform [14] and local distance statistics [39] and matched filtering. Computational implementation was carried out using MATLAB on an Intel-i7 9700 processor at 3.0GHz with 16GB of RAM.

As illustrated in Fig. 9(a)-(f) and Fig. 10(a)-(f), the point cloud data of target debris 49607 are presented for two distinct time intervals within the detection window, along with the target distance obtained using the aforementioned varied methodologies. Specifically, the absolute mean AER rate of the target during the time interval represented in Fig. 9 are 0.12 rad/s, 0.04 rad/s and $5.95 \times {10^3}$ m/s, while for the interval depicted in Fig. 10, they are 0.18 rad/s, $1.74 \times {10^{-3}}$ rad/s and $1.34 \times {10^2}$ m/s.

Fig. 9. The point cloud data for target debris 49607 and the target distance derived from various methodologies (absolute mean velocity: $5.95 \times {10^3}$ m/s, absolute mean azimuth / elevation angular velocity: 0.12 / 0.04 rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Fig. 10. The point cloud data for target debris 49607 and the target distance derived from various methodologies (absolute mean velocity: $1.34 \times {10^2}$ m/s, absolute mean azimuth / elevation angular velocity: 0.18 / $1.74 \times {10^{-3}}$ rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Similarly, as shown in Fig. 11(a)-(f) and Fig. 12(a)-(f), the point cloud data and target distance for debris 50078 are illustrated for two different periods within the detection window. The absolute mean AER rate in Fig. 11 are 0.12 rad/s, 0.02 rad/s and $5.99 \times {10^3}$ m/s, and in Fig. 12, they are 0.15 rad/s, $5.17 \times {10^{-4}}$ rad/s and $1.50 \times {10^2}$ m/s. Likewise, for debris 50578, Fig. 13(a)-(f) and Fig. 14(a)-(f) depict point cloud data and distance with mean AER rate of 0.13 rad/s, 0.02 rad/s, $4.76 \times {10^3}$ m/s and 0.15 rad/s, $5.60 \times {10^{-4}}$ rad/s, $1.53 \times {10^2}$ m/s respectively. Finally, for debris 28742, as illustrated in Fig. 15(a)-(f) and Fig. 16(a)-(f), the corresponding data for two distinct intervals show mean AER rate of 0.06 rad/s, 0.11 rad/s, $7.01 \times {10^3}$ m/s and 1.15 rad/s, 0.01 rad/s, $1.64 \times {10^2}$ m/s.

Fig. 11. The point cloud data for target debris 50078 and the target distance derived from various methodologies (absolute mean velocity: $5.99 \times {10^3}$ m/s, absolute mean azimuth / elevation angular velocity: 0.12 / 0.02 rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Fig. 12. The point cloud data for target debris 50078 and the target distance derived from various methodologies (absolute mean velocity: $1.50 \times {10^2}$ m/s, absolute mean azimuth / elevation angular velocity: 0.15 / $5.17 \times {10^{-4}}$ rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Fig. 13. The point cloud data for target debris 50578 and the target distance derived from various methodologies (absolute mean velocity: $4.76 \times {10^3}$ m/s, absolute mean azimuth / elevation angular velocity: 0.13 / 0.02 rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Fig. 14. The point cloud data for target debris 50578 and the target distance derived from various methodologies (absolute mean velocity: $1.53 \times {10^2}$ m/s, absolute mean azimuth / elevation angular velocity: 0.15 / $5.60 \times {10^{-4}}$ rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Fig. 15. The point cloud data for target debris 28742 and the target distance derived from various methodologies (absolute mean velocity: $7.01 \times {10^3}$ m/s, absolute mean azimuth / elevation angular velocity: 0.06 / 0.11 rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

Fig. 16. The point cloud data for target debris 28742 and the target distance derived from various methodologies (absolute mean velocity: $1.64 \times {10^2}$ m/s, absolute mean azimuth / elevation angular velocity: 1.15 / 0.01 rad/s). (a) True. (b) Local distance statistics. (c) Density-based filtering. (d) Hough-transform. (e) Matched filtering. (f) Proposed method.

Download Full Size | PDF

From Fig. 9 to Fig. 16, a qualitative analysis indicates that under varying target velocities, the error in the target debris distance obtained through the method proposed in this paper is notably superior to that derived from several other methodologies. The quantitative evaluation of measurement results, as shown in Table 3 and Table 4, utilized root mean square error (RMSE) and mean absolute error (MRE) within the detection window to represent the mean measurement accuracy for typical LEO space debris, thereby assessing the ranging accuracy of various methods applied to such debris. Specifically, the RMSE can be calculated as

(20)$$RMSE = \sqrt{\frac{1}{N_c} \sum_{i=1}^{N_c} (\hat d_i - d_i)^2}$$

and the MRE can be calculated as

(21)$$MRE = \frac{1}{N_c} \sum_{i=1}^{N_c} \left| \frac{\hat d_i - d_i}{d_i} \right|$$

where $N_c$ represents the total number of cumulative units in the detection window, $d_i$ denotes the mean true distance in the $i$th cumulative unit.

Table 3. RMSE (m) of different methods for target debris

View Table | View all tables in this article

Table 4. MRE ($10^{-10}$) of different methods for target debris

View Table | View all tables in this article

The results indicate that the proposed method demonstrates higher mean measurement accuracy compared to traditional methods such as local distance statistics, density-based filtering, Hough Transform-based approaches, and matched filtering. The analysis of Table 3 reveals that target debris 49607, 50078, 50578, and 28742, with respective absolute mean radial velocities of $4.11 \times {10^{3}}$ m/s, $4.49 \times {10^{3}}$ m/s, $4.60 \times {10^{3}}$ m/s and $4.91 \times {10^{3}}$ m/s in their detection windows, demonstrate a positive correlation between mean velocity and RMSE, corroborating intuitive expectations.

5. Conclusion

The SSPL detection emerges as a potent technological solution to the escalating challenge of space debris. This paper presents the development of a monostatic SSPL ranging system and emphasizes the unveiling of a novel SPL distance extraction method based on DPMM. This study formulated a state-space transition model based on the orbital dynamics of LEO space debris, utilizing state priors obtained from the EKF for efficient range gating. It also implemented a denoising strategy for the initial reflected signals, employing a collaborative blend of EMD and VMD based on modal decomposition. Overall, the study improves the accuracy of distance extraction for nonlinear, high-dynamic space debris with DPMM. Experimental results reveal that the method outperforms conventional techniques such as the density-based filtering method, local distance statistics method, Hough transform method, and matched filtering method in terms of mean ranging error, particularly for various target debris at different velocities. Additionally, this study preliminarily reveals the relationship between the mean radial velocity of space debris and the accuracy of range measurement. In summary, this method demonstrates enhanced ranging accuracy for space debris with nonlinear and high-dynamic motion characteristics, establishing a robust foundation for accurate space debris localization in future research and applications.

Appendix A: Spatial geometric model and noise and tracking and point error

A.1 Spatial geometric model

In the spatial geometric model, as illustrated in Fig. 1(b), $S$ denotes the observation satellite hosting the detection system and $T$ represents the target debris. The $O-XYZ$ coordinates correspond to the Earth-centered inertial (ECI) coordinate system, whereas the ${O_0} - {x_0}{y_0}{z_0}$ coordinates define the satellite coordinate system. The point of origin ${O_0}$ is situated at the center of mass of the observation satellite. The ${y_0}$ axis tangentially aligns with the orbital trajectory, directed along the velocity vector, while the ${x_0}$ axis extends from the Earth’s center to the satellite, and the ${z_0}$ axis is oriented orthogonally to the orbital plane, forming a Cartesian coordinate system along with the ${x_0}$ and ${y_0}$ axes. ${\boldsymbol{r}}_s$ and ${\boldsymbol{r}}_t$ respectively represent the position vectors of the observation satellite and the target debris concerning the Earth’s center within the ECI coordinate system, ${\boldsymbol{r}} = {\boldsymbol{r}}_t-{\boldsymbol{r}}_s$ represents the position vector of the target debris concerning the observation satellite.

Hypothetically, within a designated detection window during their respective orbital periods, the inclination, right ascension of the ascending node (RAAN), eccentricity, argument of perigee, and semi-major axis of both the observation satellite and the target debris remain constant, while the true anomaly varies over time. The lidar equation provides a form to estimate return signal energy levels within a specified range as [40]

(22)$${P_r} = \frac{{4E\left( t \right){\eta _T}}}{{\pi {R^2}{\varphi ^2}}}\frac{{{D^2}{\eta _R}}}{{4{R^2}}}4\pi {f_r}\cos {\theta _i}{\omega _s}^2\frac{\pi }{2}erf\left( {\frac{{\sqrt 2 l}}{{2{\omega _s}}}} \right)erf\left( {\frac{{\sqrt 2 w}}{{2{\omega _s}}}} \right)$$

where ${P_r}$ denotes the return optical power, $E$ is the transmit output optical power, ${\eta _T}$ is the transmission efficiency of transmitting optical system, $R$ is the target distance, $\varphi$ is the beam divergence angle, $D$ is the diameter of the receiving aperture, ${\eta _R}$ is the receiving optical system transmission efficiency, ${\omega _s}$ is the beam spot radius. It is assumed that the target debris have a regular rectangular shape with a length of $l$ and a width of $w$, $erf\left ( x \right ) = 1/\sqrt \pi \int _{ - x}^x {{e^{ - {t^2}}}dt}$ represents the error function.

The target’s reflectance is represented with a five-parameter bidirectional reflectance distribution function (BRDF) model [36,41] as

(23)$${f_r} = \frac{{{k_b}}}{{\pi \ln 2}} \cdot \frac{{{k_z}^2\cos \alpha_n }}{{1 + \left( {{k_z}^2 - 1} \right)\cos \alpha_n }} \cdot \exp \left[ { - \left| {{b_r}} \right| \cdot {{\left( {1 - \cos \gamma } \right)}^{{a_r}}}} \right] \cdot \frac{{G\left( {{\theta _i},{\theta _r},\varphi } \right)}}{{\cos {\theta _i}\cos {\theta _r}}} + \frac{{{k_d}}}{{\cos {\theta _i}}}$$

where the first captures specular reflection and the second accounts for diffuse reflection, with ${k_b}$ and ${k_d}$ indicating their proportions. Other variables such as normal inclination angle $\alpha _n$, surface slope distribution ${k_z}$, angle between incident vectors and the surface normal $\gamma$, local Fresnel reflection coefficient $\exp \left [ { - \left | b \right | \cdot {{\left ( {1 - \cos \gamma } \right )}^a}} \right ]$, parameters related to the refractive index of the medium ${a_r}$ and ${b_r}$ and shadowing function $G\left ( {{\theta _i},{\theta _r},\varphi } \right )$ further define the model’s applicability and robustness.

A.2 Noise

Sunlight, atmospheric scattering, lunar reflection, Earth’s thermal and lunar thermal radiation are prevalent background lights in near-Earth space. The radiation from various celestial bodies and the deep space background is ideally modeled as blackbody radiation with the Balckbody radiation law, and their peak wavelengths, being far removed from the tracking system’s sensitive band, allow for the exclusion of their radiative influence [4]. Similarly, due to the disparate thermal and spectral distributions, the blackbody radiation emitted by Earth and the Moon is deemed negligible for the detection system. Thus, solar radiation and its reflection off the Earth’s atmosphere are identified as the primary noise sources in the detection mechanism. In line with the spatial target irradiance model, the levels of solar irradiance on the entire target debris at the detector’s entry pupil and the radiant flux received by the detector are defined as

(24)$$I_{{Dsun}} = \frac{I_{{sun}} \left( \lambda_{ls}, T_b \right)' f_r}{\pi R^2} \int_A \cos \varphi_S \cos \varphi_D \, dA$$

(25)$$\phi _{sun} = {I_{Dsun}} \cdot {\eta _R} \cdot \frac{\pi }{4}{D^2}$$

where ${I_{{sun}}} \left ( \lambda _{{ls}}, T_b \right )' = \left ( \int _{\lambda _1}^{\lambda _2} I\left ( \lambda _{{ls}}, T_b \right ) d\lambda _{{ls}} \right ) \Delta F_{\Delta \lambda }$ denotes the solar irradiance transmitted through a narrowband filter, ${I_{{sun}}}\left ( \lambda _{{ls}}, T_b \right ) = R_s^2 M\left ( \lambda _{{ls}}, T_b \right )/d_{{SE}}^2$ is the solar irradiance reaching the Earth’s surface, $M\left ( \lambda _{{ls}}, T_b \right )$ is the spectral radiance of the blackbody, $\lambda _{{ls}}$ is the wavelength of the radiation source, $T_b$ is the blackbody’s absolute temperature, $R_s$ is the mean radius of the sun, $dA$ is laser scattering facet, $d_{{SE}}$ is the mean Earth-Sun distance, $\left [ \lambda _1, \lambda _2 \right ]$ is the primary visible spectrum range of the Sun, $\varphi _S$ is the angle between the normal of the space target element $dA$ and the incident light ray from the source, $\varphi _D$ is the angle between the normal of the space debris element $dA$ and the line of sight of the detector.

A.3 Tracking and pointing error

The influence of tracking and pointing errors on the reflected photons count is articulated as [42,43]

(26)$$\begin{aligned}\langle n_s \rangle &= n_s\int_{ - \infty }^\infty d\left( \delta_{\alpha_p} \right)P\left( \delta_{\alpha_p} \right)\int_{ - \infty }^\infty d\left( \delta_{\beta_p} \right)P\left( \delta_{\beta_p} \right) \exp \left[ - 2\left( \frac{\alpha_p + \delta_{\alpha_p} + \beta_p + \delta_{\beta_p}}{\theta_t} \right)^2 \right]\\ &= n_s\frac{\theta_{t\alpha_p}\theta_{t\beta_p}}{\theta_{\alpha_p}\theta_{\beta_p}}\exp \left[ - 2\left( \frac{\alpha_p^2 + \beta_p^2}{\theta_t^2} \right) \right]\exp \left[ \frac{2}{\theta_t^4}\left( \alpha_p^2\theta_{t\alpha_p}^2 + \beta_p^2\theta_{t\beta_p}^2 \right) \right] \end{aligned}$$

where $\left \langle {{n_s}} \right \rangle$ represents the mean detected photon count accounting for tracking and pointing errors, while ${n_s}$ denotes the mean count without such errors, ${\theta _t} = \varphi /2$ is laser divergence half angle, ${\alpha _p}$ is the azimuth pointing deviation and ${\beta _p}$ is the elevation pointing deviation of the telescope. The PDF of random pointing errors in the azimuth and elevation axes can be expressed as

(27)$$P\left( {{\delta _{{\alpha _p}}}} \right) = \sqrt {\frac{2}{\pi }} \frac{1}{{{\theta _{{\alpha _p}}}}}\exp \left[ { - 2{{\left( {\frac{{{\delta _{{\alpha _p}}}}}{{{\theta _{{\alpha _p}}}}}} \right)}^2}} \right]$$

(28)$$P\left( {{\delta _{{\beta _p}}}} \right) = \sqrt {\frac{2}{\pi }} \frac{1}{{{\theta _{{\beta _p}}}}}\exp \left[ { - 2{{\left( {\frac{{{\delta _{{\beta _p}}}}}{{{\theta _{{\beta _p}}}}}} \right)}^2}} \right]$$

where ${\delta _{{\alpha _p}}}$ and ${\delta _{{\beta _p}}}$ are the instantaneous random error in azimuth and elevation axes, ${\theta _{{\alpha _p}}}$ and ${\theta _{{\beta _p}}}$ respectively correspond to the tracking errors in both axes, under the following assumptions

(29)$$\frac{1}{{{\theta _{t{\alpha _p}}}^2}} = \frac{1}{{{\theta _t}^2}} + \frac{1}{{{\theta _{{\alpha _p}}}^2}}$$

(30)$$\frac{1}{{{\theta _{t{\beta _p}}}^2}} = \frac{1}{{{\theta _t}^2}} + \frac{1}{{{\theta _{{\beta _p}}}^2}}$$

Appendix B: Range gating guided by EKF

The orbital dynamics model for space debris, focusing on the mathematical representation and simulation of their orbital behavior, with the orbital dynamics equations for both the observation satellite and the debris expressed as

(31)$${\ddot{\boldsymbol r}}_{{{SatECI}}} ={-} \frac{\mu_c }{{{{\left\| {{\ddot{\boldsymbol r}}_{{{SatECI}}}} \right\|}^3}}}{\boldsymbol{{r}}_{{{SatECI}}}} + \boldsymbol{f}_{{Sat}}$$

(32)$${\ddot{\boldsymbol r}}_{{{TarECI}}} ={-} \frac{\mu_c }{{{{\left\| {{\ddot{\boldsymbol r}}_{{{TarECI}}}} \right\|}^3}}}{\boldsymbol{{r}}_{{{TarECI}}}} + \boldsymbol{f}_{{Tar}}$$

where $\mu _c$ denotes the gravitational constant of Earth, ${\boldsymbol{{r}}}_{{{SatECI}}}$ and ${\boldsymbol{{r}}}_{{{TarECI}}}$ respectively represent the position vectors of the observation satellite and the target debris in the ECI coordinate system. During the observation period, perturbing forces such as gravitational forces, magnetic Lorentz forces, solar radiation, the Poynting-Robertson effect, and plasma drag acting on the observation satellite $\boldsymbol{f}_{Sat}$ and the target debris $\boldsymbol{f}_{Tar}$ are often neglected for simplicity over short durations. The motion of the target debris relative to the observation satellite in the satellite coordinate system can be represented as [44]

(33)$${\ddot{\boldsymbol r}}_{{TSSat}} = 2 \boldsymbol{\omega}_{{TSSat}} \times {\dot{\boldsymbol r}}_{{TSSat}} + \boldsymbol{\omega}_{{TSSat}} \times (\boldsymbol{\omega}_{{TSSat}} \times {\dot{\boldsymbol r}}_{{TSSat}}) + \boldsymbol{\dot{\omega}}_{{TSSat}} \times \boldsymbol{r}_{{TSSat}}$$

where $\boldsymbol{\omega }_{{TSSat}}$ and $\boldsymbol{\dot {\omega }}_{{TSSat}}$ respectively represent the angular velocity and angular acceleration of the observation satellite. Further simplification yields [45]

(34)$$\left\{ {\begin{array}{c} {\ddot x = 2\dot \omega \dot y + {{\dot \omega }^2}x + \ddot \omega y + \frac{{2\mu_c }}{{{{\left\| {\boldsymbol{r}_{{TSSat}}} \right\|}^3}}}x}\\ {\ddot y ={-} 2\dot \omega \dot x + {{\dot \omega }^2}y - \ddot \omega y - \frac{\mu_c }{{{{\left\| {\boldsymbol{r}_{{TSSat}}} \right\|}^3}}}y}\\ {\ddot z ={-} \frac{\mu_c }{{{{\left\| {\boldsymbol{r}_{{TSSat}}} \right\|}^3}}}z} \end{array}} \right.$$

The orbit is considered circular, let $q = \dot \omega = \sqrt {\mu _c /{{\left \| {\boldsymbol{r}_{{s}}} \right \|}^3}}$ and $\ddot \omega = 0$, resulting in Eq. (34) being expressed as

(35)$$\left\{ {\begin{array}{c} {\ddot x = 3{q^2}x + 2q\dot y}\\ {\ddot y ={-} 2q\dot x}\\ {\ddot z ={-} {q^2}z} \end{array}} \right.$$

where the state vector $\boldsymbol{X}_{{k}} = {\left [ {{x},{y},{z},{{\dot x}},{{\dot y}},{{\dot z}}} \right ]^T} \in \mathbb {R}^n$ denotes the tri-dimensional positional and velocity of the target debris relative to the observation satellite, $n$ is the dimension of state vector and $k$ is the discrete time index. Integrating Eq. (35) results in the Clohessy-Wiltshire equation [46], expressed as

(36)$$\left\{ {\begin{array}{c} {x(t) = \frac{{{{\dot x}}}}{q}\sin qt - \left( {\frac{{2{{\dot y}}}}{q} + 3{x}} \right)\cos qt + 2\left( {\frac{{{{\dot y}}}}{q} + 2{x}} \right)}\\ {y(t) = 2\left( {\frac{{2{{\dot y}}}}{q} + 3{x}} \right)\sin qt + \frac{{2{{\dot x}}}}{q}\cos qt - 3\left( {{{\dot y}} + 2q{x}} \right)t + {y} - \frac{{2{{\dot x}}}}{q}}\\ {z(t) = \frac{{{{\dot z}}}}{q}\sin qt + {z}\cos qt}\\ {\dot x(t) = \left( {2{{\dot y}} + 3q{x}} \right)\sin qt + {{\dot x}}\cos qt}\\ {\dot y(t) ={-} 2{{\dot x}}\sin qt + 2\left( {2{{\dot y}} + 3h{x}} \right)\cos qt - 3\left( {{{\dot y}} + 2q{x}} \right)}\\ {\dot z(t) ={-} q{z}\sin qt + {{\dot z}}\cos qt} \end{array}} \right.$$

and the corresponding state space discretized expression can be expressed as

(37)$${\boldsymbol\Phi _{k\left| {k - 1} \right.}} = \left[ {\begin{array}{@{}cccccc@{}} {4 - 3\cos qT_d} & 0 & 0 & {\frac{1}{q}\sin qT_d} & {\frac{2}{q}(1 - \cos qT_d)} & 0\\ {6(\sin qT_d - qT_d)} & 1 & 0 & {\frac{2}{q}(\cos qT_d - 1)} & { - 3T_d + \frac{4}{q}\sin qT_d} & 0\\ 0 & 0 & {\cos qT_d} & 0 & 0 & {\frac{1}{q}\sin qT_d}\\ {3h\sin qT_d} & 0 & 0 & {\cos qT_d} & {2\sin qT_d} & 0\\ {6h(\cos qT_d - 1)} & 0 & 0 & { - 2\sin qT_d} & {4\cos qT_d - 3} & 0\\ 0 & 0 & { - q\sin qT_d} & 0 & 0 & {\cos qT_d} \end{array}} \right]$$

where $T_d$ denotes the discretized time interval. The relationship between the observation vector ${\boldsymbol{z}_k} = [\alpha, \beta, d]^T \in \mathbb {R}^m$ and the spatial position of the target debris $[x,y,z]$ can be expressed as

(38)$$\left\{ {\begin{array}{c} {\alpha = {{\tan }^{ - 1}}\left( {\frac{{{x}}}{{{y}}}} \right)}\\ {\beta = {{\tan }^{ - 1}}\left( {\frac{{{z}}}{{\sqrt {{x}^2 + {y}^2} }}} \right)}\\ {{d} = \sqrt {{x}^2 + {y}^2 + {z}^2} } \end{array}} \right.$$

where $\alpha$ denotes the azimuth angle, $\beta$ is the elevation angle and $d$ is the true distance of target. $\boldsymbol {v}_k \in \mathbb {R}^n$ is measurement noise vector and $m$ is the dimension of observation vector. The observation equation can be expressed as

(39)$${\boldsymbol{z}_k} = {h_k}\left( {\boldsymbol{X}_k^{k - 1}} \right) + {\boldsymbol{v}_k}$$

The nonlinear system observation equation is linearized at the posterior estimate $\hat {\boldsymbol{x}}_k$ via a first-order Taylor series expansion, with the corresponding Jacobian matrix of first-order partial derivatives, also referred to as observation matrix ${\boldsymbol{H}_k} \in \mathbb {R}^{m \times n}$, expressed as

(40)$${\boldsymbol{H}_k} = {\left. {\left[ {\frac{{\partial h}}{{\partial {x_n}}}} \right]} \right|_{{x_k} = {{\hat x}_k}}} = \left[ {\begin{array}{cccccc} {\frac{{{y}}}{{{x}^2 + {y}^2}}} & {\frac{{ - {x}}}{{{x}^2 + {y}^2}}} & 0 & 0 & 0 & 0\\ {\frac{{ - {x}{z}}}{{\sqrt {{y}^2 + {z}^2} \left( {{x}^2 + {y}^2 + {z}^2} \right)}}} & {\frac{{ - {y}{z}}}{{\sqrt {{y}^2 + {z}^2} \left( {{x}^2 + {y}^2 + {z}^2} \right)}}} & {\frac{{\sqrt {{x}^2 + {y}^2} }}{{{x}^2 + {y}^2 + {z}^2}}} & 0 & 0 & 0\\ {\frac{{{x}}}{{\sqrt {{x}^2 + {y}^2 + {z}^2} }}} & {\frac{{{y}}}{{\sqrt {{x}^2 + {y}^2 + {z}^2} }}} & {\frac{{{z}}}{{\sqrt {{x}^2 + {y}^2 + {z}^2} }}} & 0 & 0 & 0 \end{array}} \right]$$

In this study, the state prior at each step of the EKF, $\boldsymbol{X}_k^{-}$, serves as the center for range gating, while the range of gating is set to a fixed prior value $d_r$. Thus, range gating can be expressed as

(41)$$\left[ {a,b} \right] = \left[ {{{\left\| {\boldsymbol{X}_{k,1:3}^ - } \right\|}_2} - {d_r},{{\left\| {\boldsymbol{X}_{k,1:3}^ - } \right\|}_2} + {d_r}} \right]$$

where ${{{\left \| {\boldsymbol{X}_{k,1:3}^ - } \right \|}_2}}$ denotes Euclidean norm of the first three dimensions of the state prior at step $k$, $a$ and $b$ are respectively the lower and upper bounds of range gating.

Funding

Strategic High Tech Innovation Project of the Chinese Academy of Sciences (GQRC-19-19); National Natural Science Foundation of China (62305375); China Postdoctoral Science Foundation (2020M683600).

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

1. Z.-P. Zhang, F.-M. Yang, H.-F. Zhang, et al., “The use of laser ranging to measure space debris,” Res. Astron. Astrophys. 12(2), 212–218 (2012). [CrossRef]

2. R. Tang, Z. Li, Y. Li, et al., “Light curve measurements with a superconducting nanowire single-photon detector,” Opt. Lett. 43(21), 5488–5491 (2018). [CrossRef]

3. C. R. Phipps, “A laser-optical system to re-enter or lower low earth orbit space debris,” Acta Astronaut. 93, 418–429 (2014). [CrossRef]

4. G. Campiti, M. Tagliente, G. Brunetti, et al., “Debris detection and tracking through on-board lidar,” in International Conference on Applications in Electronics Pervading Industry, Environment and Society, (Springer, 2022), pp. 235–241.

5. M. Tagliente, G. Campiti, G. Brunetti, et al., “Spaceborne lidar for debris detection and tracking,” in Annual Meeting of the Italian Electronics Society, (Springer, 2022), pp. 172–177.

6. A. McCarthy, X. Ren, A. Della Frera, et al., “Kilometer-range depth imaging at 1550 nm wavelength using an ingaas/inp single-photon avalanche diode detector,” Opt. Express 21(19), 22098–22113 (2013). [CrossRef]

7. W. Becker, Advanced time-correlated single photon counting applications, vol. 111 (Springer, 2015).

8. Y. Yu, B. Liu, Z. Chen, et al., “A macro-pulse photon counting lidar for long-range high-speed moving target detection,” Sensors 20(8), 2204 (2020). [CrossRef]

9. D. Comaniciu and P. Meer, “Mean shift: A robust approach toward feature space analysis,” IEEE Trans. Pattern Anal. Machine Intell. 24(5), 603–619 (2002). [CrossRef]

10. R. Schnabel, R. Wahl, and R. Klein, “Efficient ransac for point-cloud shape detection,” in Computer graphics forum, vol. 26 (Wiley Online Library, 2007), pp. 214–226.

11. J. Illingworth and J. Kittler, “A survey of the Hough transform,” Computer Vision, Graphics, and Image Processing 44(1), 87–116 (1988). [CrossRef]

12. Z. Bao, Z. Li, Y. Shi, et al., “Coincidence photon-counting laser ranging for moving targets with high signal-to-noise ratio,” IEEE Photonics Technol. Lett. 26(15), 1495–1498 (2014). [CrossRef]

13. X. Du, J. Xing, and H. Huang, “Method to simulate the object tracking with photon-counting laser ranging system,” Opt. Eng. 54(11), 114103 (2015). [CrossRef]

14. W. Xue, L. Liu, X. Dai, et al., “Moving target ranging method for a photon-counting system,” Opt. Express 26(26), 34161–34178 (2018). [CrossRef]

15. D. Liu, J. Sun, W. Lu, et al., “3d reconstruction of the dynamic scene with high-speed targets for GM-APD lidar,” Opt. Laser Technol. 161, 109114 (2023). [CrossRef]

16. H. Peng, Y.-r. Wang, W.-d. Meng, et al., “Dynamic time-correlated single-photon counting laser ranging,” Optoelectron. Lett. 14(2), 129–132 (2018). [CrossRef]

17. P. Jonsson, J. Hedborg, M. Henriksson, et al., “Reconstruction of time-correlated single-photon counting range profiles of moving objects,” in Electro-Optical Remote Sensing, Photonic Technologies, and Applications IX, vol. 9649 (SPIE, 2015), pp. 11–18.

18. J. Hedborg, P. Jonsson, M. Henriksson, et al., “Time-correlated single-photon counting range profiling of moving objects,” in EPJ Web of Conferences, vol. 119 (EDP Sciences, 2016), p. 06010.

19. A.-H. Hou, Y.-H. Hu, N.-X. Zhao, et al., “Acquisition of moving target structural characteristics via photon detection,” in Ninth Symposium on Novel Photoelectronic Detection Technology and Applications, vol. 12617 (SPIE, 2023), pp. 770–775.

20. Y. Zhao, W. Hao, S. Chen, et al., “Ranging analysis of a moving target based on the dynamic instrument response function,” Opt. Lett. 48(21), 5487–5490 (2023). [CrossRef]

21. F. Wood, S. Goldwater, and M. J. Black, “A non-parametric bayesian approach to spike sorting,” in 2006 international conference of the IEEE engineering in medicine and biology society, (IEEE, 2006), pp. 1165–1168.

22. W. Hu, X. Li, G. Tian, et al., “An incremental dpmm-based method for trajectory clustering, modeling, and retrieval,” IEEE Trans. Pattern Anal. Mach. Intell. 35(5), 1051–1065 (2013). [CrossRef]

23. Y. Li, E. Schofield, and M. Gönen, “A tutorial on dirichlet process mixture modeling,” J. Mathematical Psychology 91, 128–144 (2019). [CrossRef]

24. C. Andrieu and J. Thoms, “A tutorial on adaptive mcmc,” Stat. Comput. 18(4), 343–373 (2008). [CrossRef]

25. J. Sang, B. Li, J. Chen, et al., “Analytical representations of precise orbit predictions for earth orbiting space objects,” Adv. Space Res. 59(2), 698–714 (2017). [CrossRef]

26. Y. W. Teh, “Dirichlet process,” Encyclopedia of machine learning 1063, 280–287 (2010).

27. T. S. Ferguson, “A bayesian analysis of some nonparametric problems,” Ann. Statist. 1(2), 209–230 (1973). [CrossRef]

28. Y. Kopsinis and S. McLaughlin, “Development of emd-based denoising methods inspired by wavelet thresholding,” IEEE Trans. Signal Process. 57(4), 1351–1362 (2009). [CrossRef]

29. N. Ur Rehman and H. Aftab, “Multivariate variational mode decomposition,” IEEE Trans. Signal Process. 67(23), 6039–6052 (2019). [CrossRef]

30. Z. Lin, H. Li, and C. Fang, Alternating Direction Method of Multipliers for Machine Learning (Springer, 2022).

31. A. E. Gelfand, “Gibbs sampling,” J. Am. Stat. Assoc. 95(452), 1300–1304 (2000). [CrossRef]

32. R. M. Neal, “Markov chain sampling methods for dirichlet process mixture models,” Journal of Computational and Graphical Statistics 9, 249–265 (2000). [CrossRef]

33. S. Chib, “Markov chain Monte Carlo technology,” Handb. Comput. Stat. Concepts Methods pp. 73–104 (2012).

34. D. Blackwell and J. B. MacQueen, “Ferguson distributions via pólya urn schemes,” The Annals of Statistics 1(2), 353–355 (1973). [CrossRef]

35. D. B. Dunson and J.-H. Park, “Kernel stick-breaking processes,” Biometrika 95(2), 307–323 (2008). [CrossRef]

36. W. Zhensen, X. Donghui, and X. Pinhua, “Modeling reflectance function from rough surface and algorithms,” Acta Opt. Sin. 22, 897–901 (2002).

37. S. Raychaudhuri, “Introduction to Monte Carlo simulation,” in 2008 Winter simulation conference, (IEEE, 2008), pp. 91–100.

38. S. Xia, C. Wang, X. Xi, et al., “Point cloud filtering and tree height estimation using airborne experiment data of icesat-2,” J. Remote Sens 18(6), 1199–1207 (2014). [CrossRef]

39. W. Lian, S. Li, G. Zhang, et al., “Denoising algorithm based on local distance weighted statistics for photon counting lidar point data,” in IGARSS 2019-2019 IEEE International Geoscience and Remote Sensing Symposium, (IEEE, 2019), pp. 8988–8991.

40. C. Grönwall, O. Steinvall, F. Gustafsson, et al., “Influence of laser radar sensor parameters on range-measurement and shape-fitting uncertainties,” Opt. Eng. 46(10), 106201 (2007). [CrossRef]

41. L. Bai, Z. Wu, X. Zou, et al., “Seven-parameter statistical model for brdf in the uv band,” Opt. Express 20(11), 12085–12094 (2012). [CrossRef]

42. J. J. Degnan, “Millimeter accuracy satellite laser ranging: a review,” Contributions of space geodesy to geodynamics: technology 25, 133–162 (1993). [CrossRef]

43. L. Tan, Y. Yang, J. Ma, et al., “Pointing and tracking errors due to localized deformation in inter-satellite laser communication links,” Opt. Express 16(17), 13372–13380 (2008). [CrossRef]

44. Y. Jiang and J. Schmidt, “Motion of dust ejected from the surface of asteroid (101955) bennu,” Heliyon 6(10), e05275 (2020). [CrossRef]

45. S.-G. Kim, J. L. Crassidis, Y. Cheng, et al., “Kalman filtering for relative spacecraft attitude and position estimation,” J. Guid. Control. Dyn. 30(1), 133–143 (2007). [CrossRef]

46. T. Carter and M. Humi, “Clohessy-wiltshire equations modified to include quadratic drag,” J. Guid. Control. Dyn. 25(6), 1058–1063 (2002). [CrossRef]

Wavelength	1064nm	Receiving FoV	20 $μ$ rad
Pulse Width (FWHM)	10ns	Time Resolution	1 ns
Optical Power	5W	Dark Count Rate	300 cps
Repetition Frequency	1 MHz	Timing Jitter	2 ns
Transmitting Divergence Angle	35 $μ$ rad	Dead Time	20 ns
Transmission Efficiency	0.8	Receiving Aperture Diameter	1 m
Mean Count Rate of the Noise	$1.12 \times 10^{7}$ cps	Target’s Reflectance	0.22

ID of Target Debris	Detection Window	Duration /s
COSMOS_1408_DEB_49607	11 May 2023 21:36:11.000 $\sim$ 21:36:21.000	10.000
COSMOS_1408_DEB_50078	12 May 2023 03:28:03.000 $\sim$ 03:28:13.000	10.000
COSMOS_1408_DEB_50578	11 May 2023 06:37:03.000 $\sim$ 06:37:13.000	10.000
CZ-2D_DEB_28742	11 May 2023 14:19:45.000 $\sim$ 14:19:55.000	10.000

ID of Target Debris	Density-based filtering	Local distance statistics	Hough transform	Matched filtering	Proposed method
49607	37.69	35.32	44.06	1.47	1.06
50078	36.56	33.49	47.20	1.25	1.08
50578	38.70	36.46	41.71	1.67	1.18
28742	37.68	36.61	51.80	1.48	1.24

ID of Target Debris	Density-based filtering	Local distance statistics	Hough transform	Matched filtering	Proposed method
49607	$27.5$	$33.9$	$81.1$	$5.38$	$3.29$
50078	$13.3$	$131$	$95.4$	$7.14$	$2.24$
50578	$40.5$	$557$	$71.1$	$7.11$	$2.32$
28742	$2.91$	$148$	$113$	$7.23$	$1.72$

Wavelength	1064nm	Receiving FoV	20 $μ$ rad
Pulse Width (FWHM)	10ns	Time Resolution	1 ns
Optical Power	5W	Dark Count Rate	300 cps
Repetition Frequency	1 MHz	Timing Jitter	2 ns
Transmitting Divergence Angle	35 $μ$ rad	Dead Time	20 ns
Transmission Efficiency	0.8	Receiving Aperture Diameter	1 m
Mean Count Rate of the Noise	$1.12 \times 10^{7}$ cps	Target’s Reflectance	0.22

High accuracy ranging for space debris with spaceborne single photon Lidar

Abstract

1. Introduction

2. Structure of proposed system

3. Model and proposed method

3.1 Model

3.2 Proposed method

3.2.1 Preprocessing

3.2.2 DPMM

3.2.3 DPMM based distance extraction

4. Experiment results

4.1 System parameters

4.2 Scenario

4.3 Experiment results

5. Conclusion

Appendix A: Spatial geometric model and noise and tracking and point error

A.1 Spatial geometric model

A.2 Noise

A.3 Tracking and pointing error

Appendix B: Range gating guided by EKF

Funding

Disclosures

Data availability

References

Data availability

Cited By

Figures (16)

Tables (4)

Equations (41)

Optics Express