Open Access
Add to BibSonomy
Abstract
It is extremely challenging to rapidly and accurately extract target echo photon signals from massive photon point clouds with strong background noise without any prior geographic information. Herein, we propose a fast surface detection method realized by combining the improved density-dimension algorithm (DDA) and Kalman filtering (KF), termed the DDA-KF algorithm, for photon signals with a high background noise rate (BNR) to improve the extraction of surface photon signals from spacecraft platforms. The results showed that the algorithm exhibited good adaptability to strong background noise and terrain slope variations, and had real-time processing capabilities for massive photon point clouds in large-scale detection range without prior altitude information of target. Our research provides a practical technical solution for single-photon lidar applications in deep space navigation and can help improve the performance in environments characterized by strong background noise.
© 2024 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
In the field of deep space exploration, lidar can provide range, velocity, and topographic information that is critical for spacecraft navigation, which is increasingly being considered as the primary relative navigation sensor for astronomical observation, deep space navigation, spacecraft rendezvous and landing [1–3]. Photon-counting lidar based on single-photon detection technology can be characterized by its long range, high resolution, and high accuracy. It is extensively applied in fields such as long-range remote sensing and mapping [4,5]. Using a high repetition rate, low-energy pulse lasers and photon counting technology can accurately describe the surface information of targets. This technology can effectively reduce the energy consumption of lasers and the feasibility of spaceborne applications to obtain elevation information has been demonstrated [6–8]. Single-photon lidar will become the next-generation of miniaturized, low-power, high-precision space navigation device.
Single-photon detectors have a high sensitivity and can detect extremely weak photon signals at the single-photon level. However, they are influenced by factors such as solar background radiation, atmospheric scattering, and their dark counts, including after-pulse effects. In addition, the presence of large amount of high-energy particle radiation in the space environment makes the increase in dark counts of single-photon avalanche diode (SPAD) more significant, resulting in a more severe background noise rate (BNR) and a lower signal-to-noise ratio (SNR) [9–11]. These factors result in significant noise in the photon detectors, making it crucial for space applications to effectively distinguish between noise and signal photons in environments with strong background noise. It is currently challenging to perform a rapid and accurate extraction of target echo photon signals from massive photon point clouds with strong background noise.
For the applications of spaceborne or airborne earth surface elevation mapping, due to the extremely large detection range, the usual method to limit the surface search region is using the spacecraft onboard knowledge of position and attitude to obtain the geolocation information of the target, afterwards the prior altitude information corresponding to the target geolocation is obtained with the help of onboard database such as Digital Elevation Model (DEM), then the search region can be set to reduce the amount of point cloud data [12]. Another simpler method to obtain the prior altitude information is accumulating the photon point by histogram to obtain the distance range of the signal, which is a common method in photon counting technology. However, for deep space navigation applications such as spacecraft rendezvous and landing, neither position knowledge nor DEM database can be used, and due to the high BNR, the noise almost drowns out the signal, so the prior altitude information of the target cannot be obtained by the above methods. Therefore, it is necessary to extract signals from massive point cloud data in the entire detection range. In addition, as the measurement device for deep space navigation, single-photon lidar needs to obtain the altitude information of the target in real time, so real-time performance is also one of the most important features in photon signal extraction method.
Many scholars, both domestically and internationally, have conducted in-depth research on the development of filtering algorithms for photon point clouds. They proposed excellent algorithms for ICESat-2 (Ice, Cloud, and Land Elevation Satellite-2) [6–8] and MABEL (Multiple Altimeter Beam Experimental Lidar) [13,14] photon point cloud datasets from the micropulse photon-counting lidar. These algorithms can be mainly categorized into three types: grid-based point cloud filtering methods [15–17], local statistical information-based point cloud filtering methods [18–21], and clustering-based point cloud filtering methods [22–26]. The latter two methods are based on the original point cloud. Filtering methods based on the original point cloud have a high recognition rate for signal photons and yield good extraction results; however, they are associated with higher time complexity, resulting in lower operating efficiency and longer computation time when dealing with large datasets. Clustering-based point cloud filtering methods have a relatively simple computational process and can run faster than methods based on local statistical information. However, their filtering performance is significantly affected by nonuniform point cloud distributions.
Photon point cloud filtering algorithms have seen significant developments and can achieve high processing accuracy. However, many algorithms face limitations to above challenges. The large volume of original photon point cloud data makes it challenging for these algorithms to rapidly extract surface photon signals over a large-scale detection area without any prior altitude information. Additionally, some have limited adaptability to photon signals with strong background noise, which is not conducive to subsequent applications of photon-counting lidars.
This study explored a fast surface signal extraction method for high-BNR photon point clouds in the context of extracting signal photon point clouds from large detection range. We propose a novel photon point cloud filtering method that combines the improved density-dimension algorithm (DDA) with the Kalman filtering (KF), referred to as the DDA–KF algorithm. The improved DDA algorithm takes advantage of the speed of grid-based point cloud filtering methods and the accuracy of local statistical information-based filtering methods. In the processing of photon point cloud stream data, Kalman filtering is introduced to track the surface altitude, thereby dynamically adapting the signal retrieval range. Multilevel background noise photon point clouds were generated through simulation experiments, and our method was compared with DBSCAN-related algorithms to assess its adaptability to strong background noise environments. Subsequently, a system-on-chip (SoC)-based real-time implementation solution was derived, and its real-time performance was analyzed. Finally, a photon-counting lidar prototype was designed for airborne flight experiments to validate the practical processing performance and adaptability of the proposed algorithm in complex environments. The experimental results were analyzed and quantitatively evaluated.
2. Algorithm principles
The surface-detection algorithm comprises four parts: preprocessing, signal-noise separation, post-processing, and Kalman filtering, and the algorithm flow chart is shown in Fig. 1.
- (1) Preprocessing rasterizes the raw point cloud data into 2D pseudo-image data. Pseudo-image data were transformed into density-dimensional data using an anisotropic radial basis density kernel function.
- (2) The signal-noise threshold is adaptively determined along the orbit to remove noise.
- (3) The image data after the signal-noise separation step still contains residual small-cluster noise, which requires post-processing to perform precise filtering of the data and finally restore it to the original point cloud.
- (4) The Kalman filtering algorithm is combined with the DDA to track and retrieve surface altitude.
2.1 Photon point cloud to pseudo-image
The preprocessing method of the photon point cloud filtering algorithm determines the processing efficiency of the algorithm. Preprocessing can be done using three main methods: histogram, pseudo-image, and raw-point cloud-based methods. Among these, the raw point cloud-based filtering algorithm has the highest accuracy but the lowest processing efficiency, making it unsuitable for the real-time extraction of photon point cloud data. Histogram-based filtering algorithms offer the highest efficiency but are less sensitive to terrain details and slopes, resulting in a lower extraction accuracy.
Point cloud rasterization, which transforms the original photon point cloud into a pseudo-image, where the processing units of the algorithm change from points to pixels, is an efficient method to process point clouds. Typically, rasterization can lead to a loss of point cloud accuracy; however, the degree of accuracy loss often depends on the rasterization method, specifically the grid size, which corresponds to the pixel resolution. Smaller grid sizes imply lower accuracy losses. In addition, the processed pseudo-image data can be restored to the original point cloud, which significantly reduces the accuracy loss caused by rasterization.
First, the photon point cloud streaming data were segmented along the orbit direction to determine the point cloud processing units. Grid sizes were then determined, and finally, the point cloud processing unit was rasterized to generate a pseudo-image of the photon point cloud.
2.2 Improved density-dimension algorithm
The DDA is a density-layer detection algorithm developed by Herzfeld et al. [19,20] to analyze the ICESat-2 dataset. It uses an anisotropic radial basis function (RBF) as the kernel function to quantify the photon density. These density data are treated as dimensions, enhancing the contrast between the signal photons and background noise through clustering to extract the signal photons. The proposed method is an improvement based on the DDA for photon signals extraction of strong background noise and special terrain environment. The specifics of the method details of this method can be found in Ref. [19,20], and only the improvement are described here.
For 2D data, the transformation matrix of anisotropic norm in original method is defined as:
where a represents the horizontal stretching scale.The anisotropic direction of the anisotropic RBF (hereinafter referred to as the density kernel function) is closely related to the radial motion trend of the spacecraft and terrain trend. To match the terrain slope and radial motion, it is necessary to perform a rotation operation on the density kernel function concerning the central point. The transformation matrix for this purpose is as follows:
The density kernel function with direction of θ can be expressed as follows:
where $\phi ()$ is the RBF, $\textrm{||} \cdot \textrm{ ||}$ is the Euclidean distance, c is the central point, and x is the neighboring point of c.2.2.1 Density-dimensional data generation
To identify points within the high-density clusters in a photon point cloud, the following specific steps are performed to transform them into density-dimensional data:
- (1) Rasterization of the photon point cloud: The point cloud is discretized into a grid with uniform intervals in the 2D plane. The number of photon points within each small grid is counted and used as the corresponding pixel value for that grid. To minimize the loss of geometric information between point clouds, the vertical dimension of each grid should capture the smallest terrain details, whereas the horizontal dimension should match the laser spot size.
- (2) Generation of the density kernel matrix: The size of the density kernel matrix is first determined, and its proportions are determined using the anisotropic norm. The choice of the matrix size involves choosing the neighborhood size in the density calculation for the pseudo-image. The density kernel functions are sampled and rasterized to generate a density kernel matrix.
- (3) Transformation of pseudo-image into density-dimensional data: The density kernel matrix is used as a convolution kernel, and convolution operations are performed with a sliding window on the pseudo-image of the point cloud. The results are finally normalized to generate the density-dimensional data.
The density value of point c in the direction of θ is calculated as:
where Rc represents the rectangular window region in which the density kernel matrix centered at c is located, i.e., the neighborhood range of the central point c; x is the pixel point in the pseudo-image, and z(x) is the corresponding pixel value. Finally, it is necessary to normalize the density kernel weights to ensure uniformity of the total weights.2.2.2 Multidirectional density value fusion
Multiple density kernel directions were set to match the terrain slope and radial motion of spacecraft. In deep space navigation applications, it is necessary to determine as many density kernel directions as possible to accommodate multiple signal slope trends, but each density kernel direction requires a lot of computational resources. Since here we only designed scheme for airborne validation experiments, only three angle density kernel matrices at 0° and ±30° are generated to save computational resources. First, each density kernel matrix was applied separately to the pseudo-image data, generating three corresponding density-dimensional datasets. Subsequently, the maximum density value among the corresponding pixels of the three sets of density-dimensional data was used to obtain the final density-dimensional data:
2.3 Signal-noise separation
The signal-noise separation is a crucial step in the proposed algorithm, as it determines the completeness of signal photon extraction and the amount of residual noise. The magnitude of background noise can vary due to the differences in the magnitude of background photons between day and night, variations in the reflectivity of different land surface types, differences in the densities of atmospheric aerosols and molecules, and even variations in the terrain slope. Land-surface detection algorithms must automatically determine the density threshold for noise filtering along a track.
In the original DDA, the density threshold is only determined by the maximum density value or the density value corresponding to the maximum value of the density histogram, and the density threshold is obtained by adding offset or taking the quantile of the value, which adapts to the change of noise level and has good adaptability to the data with high SNR. However, in areas with strong background noise, signal is almost drowned out by noise, resulting in an unstable SNR, it would be difficult to determine the density threshold accurately only by the maximum density value, which will cause signal cluster loss, or more noise clusters be retained as signal clusters, making the post-processing part more difficult. The histogram of photon density with background noise rate of 1 MHz and 6 MHz are shown in Fig. 3.
Fig. 3. Histogram of photon density with background noise rate of 1 MHz (a) and 6 MHz (b). there is a clear boundary between noise density and signal density in histogram (a), while noise density and signal density are almost fused together in histogram (b)
Therefore, we chose to match the density threshold with both the signal and noise intensities (where the default maximum-density region represents the signal region, and the noise intensity is equal to the density mean). The density threshold can be calculated as follows:
Here, thi is the density threshold for the i-th column along the track direction, and q is the threshold factor, typically set between 0.6 and 0.8, {di} is the set of density values in the i-th column before initializing the Kalman filter. This set spans the entire detection range, and when the Kalman filter is initialized, to achieve adaptive signal retrieval range, this set is represented by the following expression:
where j is the row index, $\hat{h}_{k + 1}^ -$ is the density-dimension-height prior estimate obtained from the Kalman filter step (described in Section 1.4.2), and r is the size of the signal retrieval window.2.4 Post-processing
2.4.1 Noise cluster removal
Due to the uneven distribution of background noise, even after signal-noise separation, small noise clusters may remain in some noise-dense areas, and they must be removed. Typically, the areas of each cluster’s connected regions are detected to eliminate small clusters. However, selecting a high area threshold may lead to the removal of some small signal clusters as noise, whereas selecting a small threshold may not completely remove the small noise clusters. Therefore, an initial smaller threshold t is chosen to remove clusters ci with an area ${\varpi _i}$ smaller than this threshold, forming a cluster set C1:
For each cluster in set C1, the area of each cluster ${\varpi _i}$ is used as a weight to calculate the weighted average $\bar{h}_w^{{C_1}}$ and weighted standard deviation $\sigma _w^{{C_1}}$ of the height of the cluster’s center $h_i^c$. When the land surface signal is extremely flat, the weighted standard deviation approaches zero. Therefore, the minimum value of ${\sigma _{\min }}$ is set to ensure the completeness of the signal clusters. The calculation formulae are as follows:
Clusters with centers within the weighted standard deviation distance from the weighted average are selected as the final cluster set ${C_2}$, as follows:
Finally, the weighted average and weighted variance of ${C_2}$ are calculated for updating in the subsequent Kalman filtering steps:
2.4.2 Identification of surface signal and atmospheric backscatter signal
This step is designed for airborne flight experiment and can be omitted due to the fact that the atmosphere backscatter is virtually non-existent in deep space applications. To extract the surface altitude signal from a spacecraft platform, it is necessary to apply the DDA to retrieve surface signals throughout the detection range because the initial surface altitude is unknown. The Kalman filter should be initialized if the surface signal is retrieved. To ensure that the retrieved surface signal is not a signal backscattered from clouds or aerosols and that the Kalman filter can be effectively initialized, the signal type must be identified.
Owing to the wider vertical distribution and weaker local density of the atmospheric backscatter signal relative to the surface signal, the average thickness T and average signal strength S of the final cluster ${C_2}$ can be calculated as follows:
The SID reflects the degree of signal divergence and can be used to distinguish surface signals from atmospheric backscatter signals.
2.5 Kalman filter
During the processing of photon point cloud streaming data, when the SNR of the point cloud is low, many noise clusters may be misidentified as signal clusters, resulting in the loss of signal photons. If the signal detection range is limited according to the signal distance extraction results obtained in the past, the loss of signal photons can be reduced and the interference of clouds in signal extraction in airborne flight experiments can be avoided.
The Kalman filter [27] is an optimal estimation method for linear systems, and its superiority has been demonstrated in fields such as control and navigation. The Kalman filter can also be used to process photon point cloud data; however, the processing flow must be designed specifically for density-dimensional data and DDA [28,29].
Therefore, the Kalman filter method was introduced to track surface targets, adaptively retrieve photon signals, and reduce the data processing volume. In navigation tasks such as spacecraft rendezvous and landing, the Kalman filter can also estimate radial altitude as well as radial velocity in real time. Therefore, the introduction of Kalman filter into the photon point cloud extraction algorithm will play a crucial role in the navigation task.
2.5.1 Combination of Kalman filter and DDA
We adopted a simple spacecraft dynamics model in which the system state vector xk includes two state variables: ranging height H and radial velocity V. The Kalman filter is combined with the DDA by taking the DDH and its rate of change as the state variables H and V, respectively. When the radial velocity of the spacecraft is zero, V is proportional to the rate of change in the surface height, that is, the slope of the terrain. If the radial velocity is not zero, V is the superposition of the radial velocity of the spacecraft and the rate of change in the surface height. The linear system model can be expressed as follows:
where $A = \left( {\begin{array}{cc} 1&{\Delta T}\\ 0&1 \end{array}} \right)$ is the system state matrix, $\Delta T$ is the time interval of the point cloud streaming segmentation data, $H\textrm{ = }\left( {\begin{array}{*{20}{c}} 1\\ 0 \end{array}} \right)$ is the measurement matrix, ${w_k}$ and ${v_k}$ are the process noise vector and measurement noise vector, respectively, and the corresponding covariance matrices are ${Q_k}$ and ${R_k}$, respectively.The following equations are involved in the prediction and updation processes of the Kalman filter:
2.5.2 Initial state, state updation, and state prediction of Kalman filter
We set the initial speed of the spacecraft to zero; therefore, the initial state of the Kalman filter is ${\hat{{\boldsymbol x}}_0} = \left( {\begin{array}{*{20}{c}} {\bar{h}_w^{{C_2}}}\\ 0 \end{array}} \right)$, where $\bar{h}_w^{{C_2}}$ is the DDH obtained from Eq. (14). The measurement value${{\mathbf Z}_k}$ was also chosen as $\bar{h}_w^{{C_2}}$ for updation:
The measurement noise covariance matrix ${{\mathbf R}_k}$ is a scalar corresponding to the weighted variance of ${C_2}$ obtained from Eq. (15):
The updation of the process noise covariance matrix ${{\boldsymbol Q}_k}$ requires calculating the variance of the state estimation values obtained from the previous n update steps as follows:
After obtaining the prior-estimated value of DDH in the next stage using Eq. (22), it is transmitted to the signal-noise separation step in the next stage to achieve an adaptive signal retrieval range.
3. Algorithm performance analysis
We analyzed the performance of the algorithm in terms of its adaptability to the signal extraction of high-BNR surface photon point cloud, its adaptability to terrain slopes, and its processing efficiency for large-scale spatial photon point clouds. We compared our algorithm with the DBSCAN-related point cloud clustering algorithms to demonstrate its superiority. DBSCAN-related algorithms contain classical DBSCAN with circular filter kernel [30], modified DBSCAN with elliptic filter kernel [22–24] and adaptive DBSCAN with slope adaptive elliptic filter kernel [25].
3.1 Photon point cloud simulation and analysis methods
To evaluate the performance of the algorithm, multiple sets of point cloud test datasets were generated through simulations using the MABEL photon-counting LiDAR raw photon point cloud dataset. MABEL is an airborne photon-counting lidar system designed by NASA and is used for the simulation and verification of the ICESat-2 spaceborne laser altimetry. The MABEL operates at an altitude of 20 km and is equipped with 532/1064 nm dual-band lasers with adjustable laser repetition rates in the range of 5–25 kHz. We only selected data in the 1064 nm band with a laser repetition rate of 5 kHz for the simulation. Table 1 presents the instrument parameters for the MABEL in the 1064 nm band [13].
Table 1. Instrument parameters for MABEL in the 1064 nm band
The MABEL dataset separates signal photons from noise photons and labels each photon type. We first used the signal marker photon as the target photon and then added different levels of noise for the surface signal extraction experiments.
According to the photon-detection model of the photon-counting lidar, the number of noisy photon events per unit time interval follows a Poisson distribution.
where k is the number of noise photon events per unit time interval, and Nn is the average number of noise photon events per unit time interval. The probability of noise photon events occurring per unit time interval can be expressed as follows:We first obtain the probability of noise photon events based on the average number of noise photons and then generate a uniformly distributed random number with a value range of [0,1] to determine the occurrence of noise photon events at each time interval:
MABEL uses a Geiger-mode avalanche photodiode (GM-APD) detector with a center wavelength of 1064 nm, and signal detection is affected by the dead-time effect. This implies that after detecting a photon signal, there is a waiting period before it can continue to detect the next photon signal. Moreover, with the increase in the background noise, the impact of the dead-time effect is pronounced. Therefore, dead-time corrections are necessary. To simulate the effects of the dead time, photons within the dead time were removed from the data based on the dead time data. Owing to the dead-time effect, both the signal and noise photons decrease in number, making it necessary to recalculate the noise level. The dead time of the silicon GM-APD is typically within 100 ns. Thus, the dead time for the detector was set to 50 ns.
We used recall (R), precision (P), and F-measure (F) as the quantitative evaluation metrics. The recall represents the proportion of correctly identified signal photons relative to the actual number of signal photons. The precision represents the proportion of correctly identified signal photons relative to those identified as signal photons. The F-measure (F) is the harmonic mean of the precision and recall. The equations used to calculate these metrics are as follows:
where TP represents the number of signal photons correctly classified as a signal, FP represents the noise photons detected incorrectly by the algorithm, and FN represents the number of signal photons incorrectly classified as noise.3.2 Adaptability of algorithm to the extraction of high-BNR point cloud
To explore the extraction adaptability of the algorithm in complex environments, we selected two types of point cloud data, one from flat terrain and the other from rough terrain, to simulate the algorithm processing effects under different background noise levels. The original data used for the simulation test were two sets of MABEL flight data from September 27, 2013.
As shown in Fig. 4, 5 and Table 2, the processing efficiencies of the DDA–KF and DBSCAN-related algorithms on relatively flat surface areas are compared under different noise count rates.
Fig. 4. Simulated point cloud, density-dimensional point cloud, and extracted signal point cloud of flat terrain under (a) 1 MHz (b) 2 MHz (c) 3 MHz (d) 4 MHz (e) 5 MHz (f) 6 MHz noise count rates. The colorbar represents the density-dimensional data corresponding to each photon point.
Fig. 5. Comparison of F-Measure between DDA–KF and DBSCAN-related algorithms on flat surface areas under different noise count rates.
Table 2. Processing efficiencies of the DDA–KF and DBSCAN-related algorithms on flat surface areas under different noise count rates
As shown in Fig. 6, 7 and Table 3, the processing effects of the DDA–KF and DBSCAN-related algorithms for rough terrain areas are compared under different noise count rates.
Fig. 6. Simulated point cloud, density-dimensional point cloud, and extracted signal point cloud of rough terrain under (a) 1 MHz (b) 2 MHz (c) 3 MHz (d) 4 MHz (e) 5 MHz (f) 6 MHz noise count rates. The colorbar represents the density-dimensional data corresponding to each photon point.
Fig. 7. Comparison of F-measure between DDA–KF and DBSCAN-related algorithms on rough terrain areas under different noise count rates.
Table 3. Processing efficiencies of the DDA–KF and DBSCAN-related algorithms for data from rough terrain areas under different noise count rates
For the flat terrain areas, the DBSCAN-related algorithms have a slight advantage under 1–4 MHz background noise, whereas under high background noise count rate, the performance of DBSCAN-related algorithms decrease, the DDA–KF algorithm has a higher F-measure than DBSCAN under 5–6 MHz background noise rate. For the rough terrain areas in mountainous regions, the F-measure of the DDA–KF algorithm was higher than that of the DBSCAN-related algorithm under different noise levels. With increasing noise levels, the advantages of the DDA–KF algorithm became more evident. Even in environments with strong background noise at 6 MHz, the F-measure was approximately 16% higher than that of the best performing algorithm in DBSCAN-related algorithms, indicating that the DDA–KF algorithm has strong adaptability in environments with strong background noise.
3.3 Real-time performance analysis
The rasterization preprocessing method is a key step in improving the real-time performance of the DDA–KF algorithm. The convolution operation in generating density-dimensional data is a computationally intensive step in the algorithm, and it is difficult to ensure real-time processing by performing these calculations on a CPU. We can take advantage of the massively parallel processing capabilities of convolution operations and apply the field-programmable gate array (FPGA) hardware acceleration to generate density-dimensional data from raw point clouds. Computationally complex post-processing and Kalman filtering were deployed on the CPU. The algorithm was implemented using the Xilinx ZYNQ UltraScale + MPSoC architecture, which integrates a quad-core ARM Cortex-A53 processor and an FPGA programmable logic section. These two components were interconnected via the AXI bus, with an ARM Cortex-A53 having a clock frequency of 1.2 GHz and the FPGA operating at 150 MHz.
To efficiently process point cloud streaming data on a chip, we designed a logical structure for the convolution pipeline in the FPGA section, allowing density-dimensional data to be generated in sync with the raw point cloud streaming data. Subsequently, the raw point cloud and density-dimensional data were transferred to the ARM DDR via the AXI bus for further processing. The main bottleneck is the processing time of the algorithm.
We compared the runtime of the DDA–KF algorithm deployed on the ZYNQ UltraScale + MPSoC against that of the classical DBSCAN based on the AMD Ryzen 7 5800 H CPU platform, because classical DBSCAN runs most efficiently of all DBSCAN-related algorithms. To simulate the large-scale surface signal detection mode on a spacecraft platform, we generated simulated point cloud data with different vertical window sizes based on the MABEL dataset, which included varying levels of background noise. This helped assess the efficiencies of both methods in handling different noise levels and window sizes. The horizontal segmentation interval for the point cloud data was fixed at 500 shots to simulate the real-time segmented processing of the photon point cloud streaming data at a processing frequency of 10 Hz and a repetition rate of 5 kHz. The runtime results of the two algorithms are listed in Tables 4 and 5.
Table 4. Processing time (ms) of the DDA–KF algorithm for point cloud data with different noise levels and vertical window sizes before initializing the Kalman filter; the processing interval is 500 shots for horizontal segmentation
Table 5. Processing time (ms) of the classical DBSCAN for point cloud data with different noise levels and vertical window sizes; the processing interval is 500 shots for horizontal segmentation
Figure 8 illustrates the runtime of the DDA–KF algorithm was 1–4 orders of magnitude lower than that of DBSCAN, and it did not decrease with increasing background noise density; it only slightly increased with increasing vertical window size. When the Kalman filter was initialized for surface tracking and signal retrieval, the runtime of the DDA–KF algorithm remained at approximately 6 ms and did not change with increasing background noise rate or vertical window size. Therefore, the proposed algorithm exhibited better performance and stability in terms of processing efficiency than the DBSCAN algorithm.
Fig. 8. Comparison between the DDA–KF and classical DBSCAN algorithms in terms of their processing speeds with logarithmic coordinates on the vertical axis.
4. Airborne flight experiment
To further validate the real-time adaptability and suitability of the DDA–KF surface signal detection method in complex environments on a spacecraft platform, a single-photon LiDAR prototype was constructed using the system parameters listed in Table 6, and the block diagram of LiDAR system is shown in Fig. 9. The start signal generated by the laser module and the echo signal generated by the SPAD module are received by the FPGA high-precision time measurement and main control module to calculate the photon flight time, the photon point cloud data stream is received by the ZYNQ UltraScale + MPSoC data processing module, and the data is processed in real time by the onboard DDA-KF algorithm, and finally the point cloud data is visualized by the upper computer. Airborne flights were conducted in the Hainan Province, China (19.8°N, 110.5°E) at an altitude of approximately 4 km, the vertical detection range was 10 km.
Table 6. System parameters of the prototype
Figure 10 shows a segment of the surface photon point cloud data collected during the flight and the surface photon signals extracted by the algorithm, with noise levels ranging from 50 kHz to 400 kHz within the point cloud data. Due to the difference between the earth's atmospheric environment and the space environment, it is difficult to obtain a higher background noise count rate in the experiment. The maximum slope of the terrain was approximately 31°. Figure 11 shows the signal photon extraction results of the photon point cloud data surface altitude trajectory obtained using the Kalman filter. Clearly, this method can adapt well to changes in noise levels along the track direction and variations in the terrain slope. These results were obtained through onboard real-time processing.
When the surface is covered by thin clouds and high-density aerosol layers, the photon echo intensity is attenuated, leading to a weaker SNR, and the signal-noise separation process is susceptible to interference from the cloud layers, causing signal loss. Figure 12 shows point cloud data in which both the surface and cloud-scattering signals are detected when a thin cloud layer is covered above the surface. The cloud and surface signals have similar morphologies, and cloud-scattering signals are easily mistaken as surface signals. The tracking of the surface altitude through the Kalman filter ensured that the surface photon signals could be extracted without interference from clouds.
Fig. 12. Raw photon point cloud data obtained during airborne flight, with the surface covered by clouds.
Fig. 13. Diagram of DBSCAN searching neighborhood. (a) Circular searching neighborhood. (b) Horizontal elliptical searching neighborhood. (c) Rotating elliptical searching neighborhood.
5. Discussion
In terms of computational efficiency, filtering algorithms based on raw point clouds experience exponential growth in computational complexity with an increase in the number and density of photons, rendering them almost impractical for larger target detection scales. However, from the above experiments, it can be seen that the grid-based DDA algorithm deployed on the SoC platform can achieve real-time signal photon extraction at an extremely high speed.
In terms of adaptability in rough terrain signal extraction with a high background noise rate, we take the DBSCAN-related algorithm as an example. Many researchers have improved this algorithm to make it more adaptable in complex environments, such as changing the circular searching neighborhood [30] to a horizontal elliptical searching neighborhood (modified DBSCAN) [22–24] and subsequently improving it to a rotating elliptical searching neighborhood (adaptive DBSCAN) [25] to match the trends of complex terrains. Figure 13 shows the diagram of DBSCAN searching neighborhood. This improvement method shares similarities with the anisotropic RBF and multidirectional anisotropic RBF in our algorithm presented in this paper.
DBSCAN-related algorithms only search for signals in the spatial dimension of the raw point cloud, however DDA increases the weight of the center point and reduces the weight of the neighboring points through RBF to enhance the contrast between the center point and the neighboring points, thereby converting the raw point cloud into density dimension, compared with DBSCAN-related algorithms, DDA is more conducive to extracting signal photons from strong background noise, and has a better effect on slope direction adaptability. Figure 14 shows the comparison of signal retrieval methods between DDA and DBSCAN.
The combination of the Kalman filter with the DDA algorithm is a new approach we have attempted. The Kalman filter serves as an auxiliary measure for processing photon point cloud streaming data from the spacecraft platform and is a critical step in tracking and retrieving ground signals. Owing to the complex terrain conditions and frequent cloud layer cover, tracking ground signals becomes challenging. To avoid tracking failures, we have considered nearly all possible scenarios, and the application of the Kalman filter has proven successful in onboard flight experiments.
6. Conclusions
In this study, we developed a real-time surface signal extraction method for a spacecraft platform based on DDA and Kalman filtering for photon point cloud streaming data with strong background noise without prior altitude information. The results of simulation experiments confirmed that the proposed method outperformed DBSCAN-related algorithms in extracting surface photon signals in environments with strong background noise and complex rough terrain. Under a strong 6 MHz background noise, the DDA–KF method achieved a signal photon extraction rate of approximately 80% in flat terrain environments and approximately 50% in rough terrain environments, demonstrating its strong adaptability to environments with strong background noise.
By leveraging the advantages of the rasterization pre-processing method and convolution with massive parallel processing, this method maintained a runtime of less than 20 ms for 500-shot along-track segmented data with a large-scale vertical detection range of 10 km. It exhibited a stable computational efficiency, which did not vary with changes in the noise level. Through airborne flight experiments conducted with a photon-ranging platform on a spacecraft operating at a laser repetition rate of 5 kHz and a vertical detection range of 10 km, the method could meet real-time processing requirements, adapt to noise level variations along the track, and facilitate surface altitude tracking and adaptive signal retrieval range after initializing the Kalman filter.
In summary, our method can successfully extract surface echo photon signals from strong background noise in large-scale photon point cloud detection ranges without requiring any prior geographic information. It can be potentially applied in rapid terrain acquisition, spacecraft landing, deep space navigation, and other related fields.
Funding
Innovation Program for Quantum Science and Technology (2021ZD0300304).
Acknowledgments
We thank the Goddard Space Flight Center (GSFC) for distributing the MABEL data.
Disclosures
The authors declare no conflicts of interest.
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. H. Zhou, Y. Chen, J. Hyyppä, et al., “An overview of the laser ranging method of space laser altimeter,” Infrared Phys. Technol. 86, 147–158 (2017). [CrossRef]
2. J. O. Woods and J. A. Christian, “Lidar-based relative navigation with respect to non-cooperative objects,” Acta Astronaut. 126, 298–311 (2016). [CrossRef]
3. M. Vetrisano and M. Vasile, “Autonomous navigation of a spacecraft formation in the proximity of an asteroid,” Adv. Sp. Res. 57(8), 1783–1804 (2016). [CrossRef]
4. R. Lussana, F. Villa, A. D. Mora, et al., “Enhanced single-photon time-of-flight 3D ranging,” Opt. Express 23(19), 24962–24973 (2015). [CrossRef]
5. Z.-P. Li, J.-T. Ye, X. Huang, et al., “Single-photon imaging over 200 km,” Optica 8(3), 344–349 (2021). [CrossRef]
6. W. Abdalati, H. J. Zwally, R. Bindschadler, et al., “The ICESat-2 laser altimetry mission,” Proc. IEEE 98(5), 735–751 (2010). [CrossRef]
7. D. D. Mclennan, “Ice, Clouds and Land Elevation (ICESat-2) Mission,” Proc. SPIE 7826, 782610 (2010). [CrossRef]
8. T. Markus, T. Neumann, A. Martino, et al., “The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2): Science requirements, concept, and implementation,” Remote Sens. Environ. 190, 260–273 (2017). [CrossRef]
9. J. Kodet, I. Prochazka, J. Blazej, et al., “Single photon avalanche diode radiation tests,” Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers Detect. Assoc. Equip. 695, 309–312 (2012). [CrossRef]
10. F. Moscatelli, M. Marisaldi, P. Maccagnani, et al., “Radiation tests of single photon avalanche diode for space applications,” Nucl. Instruments Methods Phys. Res. Sect. A Accel. Spectrometers Detect. Assoc. Equip. 711, 65–72 (2013). [CrossRef]
11. E. Anisimova, B. L. Higgins, J.-P. Bourgoin, et al., “Mitigating radiation damage of single photon detectors for space applications,” EPJ Quantum Technol. 4(1), 10–14 (2017). [CrossRef]
12. J. F. McGarry, C. C. Carabajal, J. L. Saba, et al., “ICESat-2/ATLAS Onboard Flight Science Receiver Algorithms: Purpose, Process, and Performance,” Earth Sp. Sci. 8(4), 2 (2021). [CrossRef]
13. M. Mcgill, T. Markus, V. S. Scott, et al., “The Multiple Altimeter Beam Experimental Lidar (MABEL): An Airborne Simulator for theICESat-2Mission,” J. Atmos. Ocean. Technol. 30(2), 345–352 (2013). [CrossRef]
14. K. M. Brunt, T. A. Neumann, J. M. Amundson, et al., “MABEL photon-counting laser altimetry data in Alaska for ICESat-2 simulations and development,” in Eos Transactions American Geophysical Union (2016), pp. 1707–1719.
15. L. A. Magruder, M. E. Wharton III, K. D. Stout, et al., “Noise filtering techniques for photon-counting ladar data,” in Laser Radar Technology and Applications XVII (International Society for Optics and Photonics, 2012), 8379, p. 83790Q. [CrossRef]
16. M. Awadallah, S. Ghannam, L. Abbott, et al., “Active contour models for extracting ground and forest canopy curves from discrete laser altimeter data,” in Proceedings: 13th International Conference on LiDAR Applications for Assessing Forest Ecosystems (2013), pp. 129–136.
17. S. C. Popescu, T. Zhou, R. Nelson, et al., “Photon counting LiDAR: An adaptive ground and canopy height retrieval algorithm for ICESat-2 data,” Remote Sens. Environ. 208, 154–170 (2018). [CrossRef]
18. K. H. Horan and J. P. Kerekes, “An automated statistical analysis approach to noise reduction for photon-counting lidar systems,” in IGARSS (2013), pp. 4336–4339.
19. U. C. Herzfeld, B. W. McDonald, B. F. Wallin, et al., “Algorithm for Detection of Ground and Canopy Cover in Micropulse Photon-Counting Lidar Altimeter Data in Preparation for the ICESat-2 Mission,” IEEE Trans. Geosci. Remote Sens. 52(4), 2109–2125 (2014). [CrossRef]
20. U. C. Herzfeld, T. M. Trantow, D. Harding, et al., “Surface-Height Determination of Crevassed Glaciers—Mathematical Principles of an Autoadaptive Density-Dimension Algorithm and Validation Using ICESat-2 Simulator (SIMPL) Data,” IEEE Trans. Geosci. Remote Sens. 55(4), 1874–1896 (2017). [CrossRef]
21. X. Zhu, S. Nie, C. Wang, et al., “A Ground Elevation and Vegetation Height Retrieval Algorithm Using Micro-Pulse Photon-Counting Lidar Data,” Remote Sens. 10(12), 1962 (2018). [CrossRef]
22. S. Nie, C. Wang, X. Xi, et al., “Estimating the vegetation canopy height using micro-pulse photon-counting LiDAR data,” Opt. Express 26(10), A520 (2018). [CrossRef]
23. J. Zhang and J. P. Kerekes, “An Adaptive Density-Based Model for Extracting Surface Returns From Photon-Counting Laser Altimeter Data,” IEEE Geosci. Remote Sens. Lett. 12(4), 726–730 (2015). [CrossRef]
24. J. Zhang, J. Kerekes, B. Csatho, et al., “A clustering approach for detection of ground in micropulse photon-counting LiDAR altimeter data,” in 2014 IEEE Geoscience and Remote Sensing Symposium (IEEE, 2014), pp. 177–180.
25. Z. Zhang, X. Liu, Y. Ma, et al., “Signal photon extraction method for weak beam data of ICESat-2 using information provided by strong beam data in mountainous areas,” Remote Sens. 13(5), 863 (2021). [CrossRef]
26. J. Huang, Y. Xing, H. You, et al., “Particle Swarm Optimization-Based Noise Filtering Algorithm for Photon Cloud Data in Forest Area,” Remote Sens. 11(8), 980 (2019). [CrossRef]
27. R. E. Kalman, “A New Approach To Linear Filtering and Prediction Problems,” J. Basic Eng. 82(1), 35–45 (1960). [CrossRef]
28. C. Wu, W. Xing, Z. Feng, et al., “Moving target tracking in marine aerosol environment with single photon lidar system,” Opt. Lasers Eng. 127(Apr.), 105967 (2020). [CrossRef]
29. M. Vacek, M. Peca, V. Michalek, et al., “Single photon laser altimeter data processing, analysis and experimental validation,” Adv. Sp. Res. Off. J. Comm. Sp. Res. 56(7), 1307–1318 (2015). [CrossRef]
30. M. Ester, H.-P. Kriegel, J. Sander, et al., “A density-based algorithm for discovering clusters in large spatial databases with noise,” Kdd 96(34), 226–231 (1996).
