A Fractional-Order Total Variation Regularization-Based Method for Recovering Geiger-Mode Avalanche Photodiode Light Detection and Ranging Depth Images

: High-quality image restoration is typically challenging due to low signal–to–background ratios (SBRs) and limited statistics frames. To address these challenges, this paper devised a method based on fractional-order total variation ( FOTV ) regularization for recovering Geiger-mode avalanche photodiode (GM-APD) light detection and ranging (lidar) depth images. First, the spatial differential peak-picking method was used to extract the target depth image from low SBR and limited frames. FOTV regularization was introduced based on the total variation regularization recovery model, which incorporates the fractional-order differential operator, in order to realize FOTV -regularization-based depth image recovery. These frameworks were used to establish an algorithm for GM-APD depth image recovery based on FOTV . The simulation and experimental results demonstrate that the devised FOTV -recovery algorithm improved the target reduction degree, peak signal–to–noise ratio, and structural similarity index measurement by 76.6%, 3.5%, and 6.9% more than the TV , respectively, in the same SBR and statistic frame conditions. Thus, the devised approach is able to effectively recover GM-APD lidar depth images in low SBR and limited statistic frame conditions.


Introduction
Geiger-mode avalanche photodiode (GM-APD) light detection and ranging (lidar) is a laser active imaging radar that uses GM-APD photon-level detection sensitivity to achieve long-range target weak echo signal detection, which has the advantages of high detection sensitivity, long-acting distance, and high distance resolution [1,2] and has significant applications in target detection, remote sensing, military guidance, security monitoring, and other fields. However, GM-APD lidar is seriously disturbed by noise, such as daylight and atmospheric backscattering, and when the signal-to-background ratio (SBR, the ratio of target echo signal photons to background noise photons in the gate) of the laser echo signal is low, the target echo signal is easily drowned in the noise, resulting in the recovered target depth image containing a large amount of noise. In addition, GM-APD is the first photon trigger system, and a single detection only corresponds to a single echo signal photon; thus, it is impossible to distinguish it from a signal photon and noise photon, so it is necessary to make multiple statistical measurements and use the depth recovery algorithm to obtain the depth image of the target. In the field of target detection, security monitoring and other applications, the relative motion between the lidar itself and the target often occurs, and the integration time of the target detection is very short, which cannot guarantee a sufficient number of GM-APD lidar statistics, and the recovery performance of most depth recovery algorithms decreases at this point. Therefore, there is an urgent need to investigate the total-variance recovery provides higher quality recovering results by introducing fractionalorder regularization terms, making it a powerful and flexible tool in the field of image recovery, which can better balance the needs of noise suppression and detail preservation. Currently, image recovery based on fractional-order total-variance regularization has been widely used in 2D image recovery [21][22][23][24][25][26], but less research has been reported in the field of GM-APD lidar distance image recovery. In order to restore GM-APD lidar depth images with low SNRs and frame rates, this paper devises a FOTV-regularization-based method for depth image restoration. First, a spatial-domain differential peak-picking method is used to extract target depth images from low SNR and low frame rate GM-APD lidar data. Second, fractional-order differential operators are introduced into the TV-regularizationbased image recovery model to construct a GM-APD depth image recovery model based on FOTV regularization. This fractional-order recovery model is then used to process target depth images and restore the true values of noisy points in depth images. The effectiveness of the devised algorithm is verified through simulations and experiments.
The remainder of this paper is organized as follows. Section 2 describes the construction and solution of the FOTV regularization recovery model. Section 3 describes the FOTV-based method for GM-APD depth-image recovery. Section 4 describes the evaluation of the recovery performance of the devised algorithm, with the evaluation metrics being the K, PSNR, and SSIM. Section 5 presents the concluding remarks.

TV Regularization Recovery Model
The variational method is typically used to solve ill-posed inverse problems by converting the problem into a functional problem. Variational methods can be applied to image recovery because it is an inverse problem. The TV regularization method can be mathematically expressed as follows [27]: where |(∇u) i,j | = (∇ 1 u) 2 i,j + (∇ 2 u) 2 i,j and ∇ represents the gradient difference. For an image with a resolution of M × N, the horizontal gradient (∇ 1 u) i,j and vertical gradient (∇ 2 u) i,j can be represented as follows: The energy functional form of the TV regularization recovery model can be obtained by introducing a Lagrange multiplier λ to generate an unconstrained extremum model: In Equation (4), the first term is the regularization term with the differential operator, which aims to extract the inherent structural features of a noisy image during recovery and remove the noise points that do not belong to the image features. The second term is the data fidelity term, which is used to minimize the difference between the denoised and original images, thereby ensuring that the two images become infinitely close and that the fidelity of the image can be preserved without distortion. u and f represent the denoised and noisy images, respectively. The parameter λ > 0 is the Lagrange multiplier or regularization weight parameter, which is used to balance the regularization and data fidelity terms in the recovery process. The TV-recovery algorithm exhibits a satisfactory recovery performance and edge preservation ability due to its anisotropic diffusion characteristics. However, the denoised images typically suffer from the "staircase effect" or virtual edges. Moreover, as shown in Equations (2) and (3), the TV regularization term considers only the first-order differences at the pixel points, i.e., it establishes connections between neighboring pixels without considering distant pixels, and thus it cannot fully exploit the information of the neighboring pixels.

TV-Regularization Recovery Model
The introduction of fractional derivatives can enable more accurate descriptions of certain phenomena than is possible with the above-described method. Specifically, many physical phenomena with memory and hereditary characteristics can be effectively described by introducing fractional-order systems, and thus, such systems have attracted significant attention [28].
The commonly used definitions of fractional derivatives are the Grumwald-Letnikov, Riemann-Liouville, and Caputo definitions [29,30]. The Grumwald-Letnikov definition has a simple formula, which facilitates numerical computation; order flexibility and tunability, enabling adaptation to different signal processing needs; and insensitivity to noise, enabling noise suppression. Therefore, the Grumwald-Letnikov definition is suitable for imageprocessing applications.
The Grumwald-Letnikov fractional derivative of a real function f (x), with x ∈ [a, t] and a < t, a ∈ R, is defined as follows: where α k = Γ(α+1) Γ(k+1)Γ(α−k+1) and h represents the step size for differentiation. The equivalent expression for the v-order fractional derivative of a one-dimensional (1D) signal f (t) over the interval [a, t], with uniform partitioning using h = 1 and The 2D image signal f (x, y) is defined by assuming that the fractional derivatives in the xand y-axis directions are separable under certain conditions. Using the separability of the Fourier transform, the fractional calculus framework can be extended from one to two dimensions. The fractional derivatives for the xand y-axes can be obtained by uniformly partitioning the 2D image signal f (x, y) based on the pixel spacing, as follows: The coefficient w v m for the v-order fractional derivative can be obtained from Equations (3) and (4) (where N is the number of polynomial terms) as follows:

FOTV-Regularization Recovery Model
Because the TV regularization model uses only first-order differences and cannot fully exploit the information of neighboring pixels, it may result in "block artifacts" in images. Therefore, fractional-order differences must be used to process neighboring pixels. Fractional-order differences possess memory, which allows them to not only use information from adjacent pixels but also incorporate information from distant pixels. Therefore, theoretically, fractional-order differences can help capture more abundant pixel information and mitigate block artifacts.
Fractional calculus is an extension of integer-order calculus, and the fractional TV regularization term is obtained from the TV regularization term TV(u) = ||∇u|| 1 = ∑ i,j |(∇u) i,j | as follows: D v is a linear operator for the fractional-order derivatives; and D v 1 , D v 2 represent the fractional-order derivative operators in the horizontal and vertical directions, respectively. Equations (5)- (7) can be used to obtain the fractional-order finite forward difference scheme: Using the matrix approximation method, Equation (11) can be expressed as follows: Matrix D has the following form: Substituting Equations (10)- (13) into Equation (4) yields the following mathematical model for FOTV regularization: When v = 1, D is a sparse banded matrix consisting of two diagonal elements, and the gradient information is determined by only the two adjacent points. When v is not an integer, D is a lower triangular matrix, as shown in Equation (13). That is when calculating the fractional-order derivative at the kth point, all points before k are involved. Thus, the fractional-order derivative is a global operator with a long memory, which distinguishes it from integer-order derivatives.

Solution of the FOTV-Regularization Recovery Model
Traditional optimization algorithms cannot be used to solve the TV-regularization recovery model because it often exhibits non-smooth and non-convex characteristics. Thus, this model is typically solved using iterative algorithms. Among such algorithms, the split Bregman algorithm, which transforms the original problem into a series of subproblems, is commonly used, as it can efficiently and rapidly obtain solutions [31][32][33].
If H(·) and φ(·) are convex functions and φ(·) is non-differentiable, the constrained optimization problem can be formulated as follows: This problem can be transformed into an unconstrained problem: where λ is the penalty parameter, and the auxiliary variable b is introduced to facilitate computation.
The solution for the anisotropic FOTV-regularization recovery problem considered in this study is: The auxiliary variable d is introduced as follows: The penalty term is added as follows: Then, the solution is obtained using Bregman iterations as follows: During the Bregman iteration process, the solution approaches the optimal value b k . To solve this minimization problem, the iterative format for the subproblems of d and u can be obtained by solving them separately as follows: The subproblem of d can be solved using the shrink operation, namely: where shrink(x, y) = x |y| · max(|x| − y, 0). The subproblem of u can be solved as follows: Because Equation (24) is strictly diagonally dominant, it can be solved using the Gauss-Seidel method. Let u k+1 The flow of the algorithm is shown in Algorithm 1.
Algorithm 1: Split Bregman Algorithm for Solving Recovery Models for Anisotropic FOTV Regularization

GM-APD Depth Image FOTV Restoration Algorithm
3.1. Depth-Image Extraction from Low SBR and Few-Frame Data Using a Spatial-Domain Differential Peak-Picking Method Equation (26) defines the triggering probability of a single-photon detector in the jth bin within the gating window of GM-APD, assuming that the noise photons in the GM-APD from the lidar follow a uniform distribution and only noise photons trigger the single-photon detector in the absence of a target in the gating window. Figure 1 shows the probability density curve of triggering noise photons, where the probability of triggering noise photons is the probability that each time-interval photon in the GM-APD lidar range gate triggers the GM-APD to generate an avalanche event and output an avalanche current. In Figure 1, the gate width of GM-APD lidar is 70 ns, and the time resolution of GM-APD is 1 ns, then the gate is divided into 70 bins, and the width of each bin is 1 ns; thus, the vertical coordinates of each GM-APD bin represent the probability of being triggered.
The triggering probability function of GM-APD is a monotonically decreasing function, as shown in Equation (27): Figure 1. Probability curve of triggered noise photons.
In Figure 1, the gate width of GM-APD lidar is 70 ns, and the time resolution of GM-APD is 1 ns, then the gate is divided into 70 bins, and the width of each bin is 1 ns; thus, the vertical coordinates of each GM-APD bin represent the probability of being triggered. The triggering probability function of GM-APD is a monotonically decreasing function, as shown in Equation (27): Because the noise photons are uniformly distributed within the gating window, the triggering probability of noise photons can be defined as follows: Substituting Equation (28) into Equation (27) yields: According to Equation (29), when GM-APD is triggered only by noise, its triggering probability curve is a logarithmically decreasing curve (it decreases rapidly at first and gradually in later stages).
The triggering probability for each time interval can be calculated by setting a target at the 20th bin as follows: where S(i) and N(i) denote the number of signal and noise photons in the ith bin, respectively. Figure 2 shows the triggering probability density curve of the GM-APD under this condition.  The triggering probability density curve of the GM-APD transitions from a logarith mically decreasing function to a convex function near the bin in which the signal photons are located (determined by the laser pulse width). When the photon arrives at time t 1 ( ) ( 1) 0 P t P t Δ = − − > , the derivative of the GM-APD trigger probability density curve can be calculated, and its distribution curve is shown in Figure 3. The triggering probability density curve of the GM-APD transitions from a logarithmically decreasing function to a convex function near the bin in which the signal photons are located (determined by the laser pulse width). When the photon arrives at time t, ∆1 = P(t) − P(t − 1) > 0, the derivative of the GM-APD trigger probability density curve can be calculated, and its distribution curve is shown in Figure 3. The triggering probability density curve of the GM-APD transitions from a logarithmically decreasing function to a convex function near the bin in which the signal photons are located (determined by the laser pulse width). When the photon arrives at time t, 1 ( ) ( 1) 0 P t P t Δ = − − > , the derivative of the GM-APD trigger probability density curve can be calculated, and its distribution curve is shown in Figure 3. In Figure 3, the Y-coordinate is the value obtained by deriving the frequency-triggered histogram within the selective gate.
When the GM-APD is used for cumulative detection, the frequency-triggered histogram within the gating window is identical to its probability density curve distribution. At this point, the target position can be expressed as follows: where diff represents the first-order derivative. During the detection process, the target echo signal occupies multiple bins. However, the probability of triggering is highest at the time at which the first photon arrives. The value of the first-order derivative at this time is the largest. Therefore, this method can help decrease the ranging error caused by the laser pulse width.

FOTV-Regularization Recovery Algorithm
Because a target typically occupies multiple pixels in the focal plane of a detector, the recovery accuracy can be enhanced by introducing fractional-order operators to establish In Figure 3, the Y-coordinate is the value obtained by deriving the frequency-triggered histogram within the selective gate.
When the GM-APD is used for cumulative detection, the frequency-triggered histogram within the gating window is identical to its probability density curve distribution. At this point, the target position can be expressed as follows: where diff represents the first-order derivative. During the detection process, the target echo signal occupies multiple bins. However, the probability of triggering is highest at the time at which the first photon arrives. The value of the first-order derivative at this time is the largest. Therefore, this method can help decrease the ranging error caused by the laser pulse width.

FOTV-Regularization Recovery Algorithm
Because a target typically occupies multiple pixels in the focal plane of a detector, the recovery accuracy can be enhanced by introducing fractional-order operators to establish connections between multiple pixels. However, when the depth image of a target contains a large amount of noise, blindly establishing connections between pixels can increase the influence of noise on the current pixel. To address this problem, noise judgment is introduced to retain the depth values of target pixels and perform only fractional-order recovery on the noise points in the depth image. Furthermore, when multiple distance values exist in the target region, the establishment of too many connections between pixels may reduce the accuracy of the depth value of the current pixel. In this study, a 5 × 5 neighborhood pixel calibration method was used to strengthen the connections of the current pixel with the surrounding pixels while decreasing the influence of distant pixels with different depth values, thereby improving the recovery accuracy.
The recovery process involves the following steps: Step one. Identify the noise points. Generally, points with maximum or minimum values represent noise in an image, as an image is typically composed of pixels with similar and continuous values. However, certain extreme points may be edge points instead of noise points. To accurately detect noise and retain edge information in an image, a fractional-order gradient judgment is introduced. Figure 4 shows the neighborhood pixels of point f (x, y).
The recovery process involves the following steps: Step one. Identify the noise points. Generally, points with maximum or minimum values represent noise in an image, as an image is typically composed of pixels with similar and continuous values. However, certain extreme points may be edge points instead of noise points. To accurately detect noise and retain edge information in an image, a fractional-order gradient judgment is introduced. Figure 4 shows the neighborhood pixels of point ( , ) f x y .  To determine whether f (x, y) is a noise point, the fractional-order gradients in eight directions around the point must be calculated.
The direction gradient is calculated as follows: where If the derivative values in all eight directions are greater than a given threshold T, then the current pixel is a noise point; otherwise, it is a signal point.
Step two. Compute the local fractional-order gradient operator D v loc . Fractional-order local TV regularization and fractional-order global TV regularization are commonly used image-recovery methods, which differ in terms of the range of image structures they consider during recovery. The former method is typically performed within a small window, and thus it can better preserve the details in the image than the latter method. The latter method considers the structure of an entire image, and although it produces better overall smoothing effects than the local method, it may lose several details. Therefore, when recovering depth images, the focus is on image details, such as the depth value information. Thus, the local fractional-order gradient operator matrix D v loc is selected as follows: Step three. Solve the recovery model with anisotropic FOTV regularization using the split Bregman algorithm.

Simulation and Experimental Verification
Computer simulations and experiments were conducted to evaluate and analyze the performance of the devised approach. The evaluation metrics were K, PSNR, and SSIM.

Simulation and Experimental Verification
Computer simulations and experiments were conducted to evaluate and analyze the performance of the devised approach. The evaluation metrics were K, PSNR, and SSIM.

Evaluation Metrics
The target reduction degree is denoted as K. The higher its value, the better the depth image recovery.
where d s is the reconstructed distance value, d s is the standard distance value, d b is the allowed error distance value, n is the total number of target pixels, and m is the number of pixels with allowed error distance values. K measures the ratio of the number of correctly recovered target pixels to the total number of target pixels. K is calculated by preserving the target recovery rate up to the fourth decimal place and retaining its integer part when multiplied by the total number of target pixels.

PSNR
The PSNR evaluates the difference between two corresponding pixels in two images based on the mean square error. The distortion rate of the restored image is evaluated using a standard image as the reference. A higher PSNR corresponds to higher fidelity of the image. The PSNR can be expressed as follows: where x(i, j) represents the depth image data with noise and y(i, j) represents the denoised depth image data.

SSIM
The SSIM is a stability-optimized system evaluation metric that is used to quantify the structural similarity between two images. The SSIM can measure the distortion level of an image or the similarity between two images. The performance of an algorithm can be evaluated using a standard image as the reference. A high SSIM corresponds to a high similarity between the two images and shows that the algorithm can effectively denoise, restore, and maintain image fidelity. The SSIM can be expressed as follows: The SSIM is defined based on three comparison measurements between the x and y samples: luminance (l), contrast (c), and structure (s).

Depth Image Extraction
The Monte Carlo method [34] was used to simulate the GM-APD lidar echo data. The system parameters were set as follows: the imaging resolution of the GM-APD detector was 64 × 64 pixels; the transmittance values of the receiving and transmitting lenses were 0.8 and 0.9, respectively; the laser pulse energy was 10 nJ; the laser wavelength was 1064 nm; the target reflectivity was 0.1; and the receiving aperture was 50 mm. The time resolution of the detector was 1 ns, and the gate width was set to 70 m. A cup model was simulated, and the distance between the cup model and the detection system was set as 20 m, as shown in Figure 6.

Depth Image Extraction
The Monte Carlo method [34] was used to simulate the GM-APD lid system parameters were set as follows: the imaging resolution of the G was 64 × 64 pixels; the transmittance values of the receiving and transm 0.8 and 0.9, respectively; the laser pulse energy was 10 nJ; the laser wav nm; the target reflectivity was 0.1; and the receiving aperture was 50 m lution of the detector was 1 ns, and the gate width was set to 70 m. A simulated, and the distance between the cup model and the detection s 20 m, as shown in Figure 6. In the depth image shown in Figure 6, the horizontal and vertical c sent the pixel locations of the image, and the color values represent the the target from the GM-APD camera.
The peak-picking method and devised spatial-domain differen method were applied to process echo data with SBRs of 0.1, 0.11, and 0. the imaging results of the two methods for 20, 40, 60, 80, and 100 frames In the depth image shown in Figure 6, the horizontal and vertical coordinates represent the pixel locations of the image, and the color values represent the actual distance of the target from the GM-APD camera.
The peak-picking method and devised spatial-domain differential peak-picking method were applied to process echo data with SBRs of 0.1, 0.11, and 0.2. Figure 7 shows the imaging results of the two methods for 20, 40, 60, 80, and 100 frames.  When the SBR is 0.2, both the peak-picking method and devised algorithm can obtain the depth image of the cup, which becomes clearer with an increase in the number of statistical frames. For a given number of statistical frames, the depth image obtained by the devised algorithm is superior to that obtained by the peak-picking method. When the SBRs are 0.1 and 0.11, the peak-picking method cannot obtain a clear depth image of the cup, whereas the devised algorithm can obtain a relatively clear depth image. The depth image quality was evaluated by 1000 sets of Monte Carlo repetition experiments at different statistical frame numbers using the mean values of three metrics: target recovery, peak signal-to-noise ratio, and structural similarity, as shown in Figures 8-10.
The devised algorithm outperforms the peak-picking method in terms of K and SSIM, and this outperformance increases as the SBR decreases. Moreover, when there are more than 50 frames, the PSNR of the devised algorithm is significantly better than that of the peak-picking method. When the SBR is 0.2, both the peak-picking method and devised algorithm can obtain the depth image of the cup, which becomes clearer with an increase in the number of statistical frames. For a given number of statistical frames, the depth image obtained by the devised algorithm is superior to that obtained by the peak-picking method. When the SBRs are 0.1 and 0.11, the peak-picking method cannot obtain a clear depth image of the cup, whereas the devised algorithm can obtain a relatively clear depth image. The depth image quality was evaluated by 1000 sets of Monte Carlo repetition experiments at different statistical frame numbers using the mean values of three metrics: target recovery, peak signal-to-noise ratio, and structural similarity, as shown in Figures 8-10.    When the SBR is 0.2, both the peak-picking method and devised algorithm can obtain the depth image of the cup, which becomes clearer with an increase in the number of statistical frames. For a given number of statistical frames, the depth image obtained by the devised algorithm is superior to that obtained by the peak-picking method. When the SBRs are 0.1 and 0.11, the peak-picking method cannot obtain a clear depth image of the cup, whereas the devised algorithm can obtain a relatively clear depth image. The depth image quality was evaluated by 1000 sets of Monte Carlo repetition experiments at different statistical frame numbers using the mean values of three metrics: target recovery, peak signal-to-noise ratio, and structural similarity, as shown in Figures 8-10.    The devised algorithm outperforms the peak-picking method in terms of K and SSIM, and this outperformance increases as the SBR decreases. Moreover, when there are more than 50 frames, the PSNR of the devised algorithm is significantly better than that of the peak-picking method.

Selection of optimal fractional order for fractional calculus
In order to explore the influence of fractional order on the recovery performance of the same statistical frame number and different SBR simulated echo data, the statistical frame number is set to 50 frames, and the recovery algorithm in this paper is used to analyze the simulated echo data with SBRs of 0.1, 0.11, and 0.2, respectively. For recovery processing, 1000 groups of Monte Carlo repeated experiments were carried out, and the average value of the three metrics of target reduction degree, peak signal-to-noise ratio, and structural similarity was taken to evaluate the depth image quality. The results are shown in Figures 11-13.

Selection of optimal fractional order for fractional calculus
In order to explore the influence of fractional order on the recovery performance of the same statistical frame number and different SBR simulated echo data, the statistical frame number is set to 50 frames, and the recovery algorithm in this paper is used to analyze the simulated echo data with SBRs of 0.1, 0.11, and 0.2, respectively. For recovery processing, 1000 groups of Monte Carlo repeated experiments were carried out, and the average value of the three metrics of target reduction degree, peak signal-to-noise ratio, and structural similarity was taken to evaluate the depth image quality. The results are shown in Figures 11-13. The devised algorithm outperforms the peak-picking method in terms of K and SSIM, and this outperformance increases as the SBR decreases. Moreover, when there are more than 50 frames, the PSNR of the devised algorithm is significantly better than that of the peak-picking method.

Selection of optimal fractional order for fractional calculus
In order to explore the influence of fractional order on the recovery performance of the same statistical frame number and different SBR simulated echo data, the statistical frame number is set to 50 frames, and the recovery algorithm in this paper is used to analyze the simulated echo data with SBRs of 0.1, 0.11, and 0.2, respectively. For recovery processing, 1000 groups of Monte Carlo repeated experiments were carried out, and the average value of the three metrics of target reduction degree, peak signal-to-noise ratio, and structural similarity was taken to evaluate the depth image quality. The results are shown in Figures 11-13.  The results show that the fractional order of the devised algorithm considerably affects the recovery results of the simulated echo data with different SBRs but the same statistical frame number. Moreover, the effect of the fractional order on the evaluation indices depends on the SBRs. Table 1 shows the optimal orders corresponding to each evaluation index for SBRs of 0.1, 0.11, and 0.2. For the same statistical frame number, the optimal orders for different evaluation metrics of the echo signal with different SBRs are different. In contrast, the optimal orders are similar when the SBR is low.    The results show that the fractional order of the devised algorithm considerably affects the recovery results of the simulated echo data with different SBRs but the same statistical frame number. Moreover, the effect of the fractional order on the evaluation indices depends on the SBRs. Table 1 shows the optimal orders corresponding to each evaluation index for SBRs of 0.1, 0.11, and 0.2. For the same statistical frame number, the optimal orders for different evaluation metrics of the echo signal with different SBRs are different. In contrast, the optimal orders are similar when the SBR is low.   The results show that the fractional order of the devised algorithm considerably affects the recovery results of the simulated echo data with different SBRs but the same statistical frame number. Moreover, the effect of the fractional order on the evaluation indices depends on the SBRs. Table 1 shows the optimal orders corresponding to each evaluation index for SBRs of 0.1, 0.11, and 0.2. For the same statistical frame number, the optimal orders for different evaluation metrics of the echo signal with different SBRs are different. In contrast, the optimal orders are similar when the SBR is low.
In order to explore the influence of the fractional order on the recovery performance of the simulated echo data with the same SBR and different statistical frame numbers, first set the SBR = 0.11 of the echo signal and then set the fractional order to 0.1, 0.5, 0.9, 1.3, and 1.7, respectively. The depth image recovery is performed on the simulated echo data with a statistical frame number of 10-100 frames, and the average value of the three metrics of target reduction degree, peak signal-to-noise ratio, and structural similarity obtained by 1000 Monte Carlo repeated experiments is used to recover the depth image. This evaluation was carried out, and the results are shown in Figures 14-16. In order to explore the influence of the fractional order on the recovery performance of the simulated echo data with the same SBR and different statistical frame numbers, first set the SBR = 0.11 of the echo signal and then set the fractional order to 0.1, 0.5, 0.9, 1.3, and 1.7, respectively. The depth image recovery is performed on the simulated echo data with a statistical frame number of 10-100 frames, and the average value of the three metrics of target reduction degree, peak signal-to-noise ratio, and structural similarity obtained by 1000 Monte Carlo repeated experiments is used to recover the depth image. This evaluation was carried out, and the results are shown in Figures 14-16.   In order to explore the influence of the fractional order on the recovery performance of the simulated echo data with the same SBR and different statistical frame numbers, first set the SBR = 0.11 of the echo signal and then set the fractional order to 0.1, 0.5, 0.9, 1.3, and 1.7, respectively. The depth image recovery is performed on the simulated echo data with a statistical frame number of 10-100 frames, and the average value of the three metrics of target reduction degree, peak signal-to-noise ratio, and structural similarity obtained by 1000 Monte Carlo repeated experiments is used to recover the depth image. This evaluation was carried out, and the results are shown in Figures 14-16.   The fractional order of the devised algorithm considerably affects the recovery of the simulated echo data with the same SBR. Moreover, the effect of the fractional order on the evaluation indices depends on the statistical frame numbers. To clarify this aspect, we analyzed the effect of the fractional order on the evaluation indices for data containing 40 and 70 frames. The results are summarized in Table 2.   The fractional order of the devised algorithm considerably affects the recovery of the simulated echo data with the same SBR. Moreover, the effect of the fractional order on the evaluation indices depends on the statistical frame numbers. To clarify this aspect, we analyzed the effect of the fractional order on the evaluation indices for data containing 40 and 70 frames. The results are summarized in Table 2. The optimal orders corresponding to different evaluation metrics vary with the statistical frame numbers, and the optimal orders corresponding to K and SSIM are similar.
Overall, the optimal orders corresponding to different evaluation metrics for echo data with different statistical frame numbers and SBRs are different. The target reduction degree (K) is an evaluation index that better reflects the target detection and identification capability of a radar system, which is able to evaluate the performance of the radar system by comparing the difference between the target detected by the radar system and the real target. Because K is the most important evaluation metric for the considered application, the fractional order corresponding to the highest K is selected as the optimal fractional order for recovering different types of echo data.

2.
FOTV-recovery algorithm To evaluate the recovery performance of the devised algorithm for low SBR and fewframe target echo data, the devised algorithm and TV-recovery algorithm were applied to simulate echo data with a SBR of 0.1 and 30, 50, and 70 frames, respectively. The recovery results are shown in Figure 17. It can be seen that the target area is smoother after recovery via the devised algorithm. In order to evaluate the image quality using quantitative metrics, 1000 Monte Carlo repetition experiments were conducted at different statistical frame numbers, and the mean values of each obtained metric are shown in Table 3. It can be seen that the target area is smoother after recovery via the devised algorithm. In order to evaluate the image quality using quantitative metrics, 1000 Monte Carlo repetition experiments were conducted at different statistical frame numbers, and the mean values of each obtained metric are shown in Table 3. It can be seen that both the devised FOTV-recovery algorithm and the traditional TV-recovery algorithm can improve the three metrics of K, PSNR, and SSIM. When the number of statistical frames is 30, the K, PSNR, and SSIM values of the FOTV recovery are, respectively, 10.1%, 14.4%, and 2.99% better than those of the TV.
To verify the level of advancement of the algorithm, this paper was compared with the frame-less detection algorithm of [8] GM-APD lidar, as shown in Table 4. The lower the SBR, the fewer the statistical frames and the more difficult the target recovery is. It can be seen from Table 4 that when the SBR is 0.1, and the number of statistical frames is 100, the three indicators, namely K, PSNR, and SSIM, are better than [8] when the SBR is 0.12, and the number of statistical frames is 200. It can be seen that the performance of the algorithm devised in this paper improves by a certain degree.

Experimental Platform
A 64 × 64 array GM-APD lidar system was established to validate the performance of the devised algorithm, as shown in Figure 18. The lidar system consisted of a 1064 nm fiber laser, a 64 × 64 array GM-APD, and transmission and reception optical paths with a transmittance of 0.9. The fields of view for transmission and reception were 0.8 The maximum output energy of the laser was 100 µJ, the pulse width was 5 ns, and the repetition frequency was 10 kHz.
The laser emits a beam that is transmitted through the optical system and illuminates the target area. The feedback signal from the laser triggers the GM-APD to begin timing. The laser beam is diffusely reflected from the target surface and collected by the optical system onto the focal plane of the GM-APD, which stops the timing process. The readout circuitry completes the time-to-digital conversion of the laser photon flight time and transfers the data to the host computer. The host computer then extracts and denoises the depth image of the target and displays it.
A 64 × 64 array GM-APD lidar system was established to validate the performance of the devised algorithm, as shown in Figure 18. The lidar system consisted of a 1064 nm fiber laser, a 64 × 64 array GM-APD, and transmission and reception optical paths with a transmittance of 0.9. The fields of view for transmission and reception were 0.8° × 0.8°. The maximum output energy of the laser was 100 µJ, the pulse width was 5 ns, and the repetition frequency was 10 kHz. The laser emits a beam that is transmitted through the optical system and illuminates the target area. The feedback signal from the laser triggers the GM-APD to begin timing. The laser beam is diffusely reflected from the target surface and collected by the optical system onto the focal plane of the GM-APD, which stops the timing process. The readout circuitry completes the time-to-digital conversion of the laser photon flight time and transfers the data to the host computer. The host computer then extracts and denoises the depth image of the target and displays it.

Outdoor Experiment
To validate the recovery performance of the devised algorithm, imaging experiments were conducted on a residential building at a distance of 1.531 km under strong sunlight. The target scene is shown in Figure 19a. To obtain the ideal target depth image, the same target area was imaged and detected at night using the peak-picking method with 5000 frames to accumulate the data. The resulting image was used as the ideal target depth image of this scene, as shown in Figure 19b.

Outdoor Experiment
To validate the recovery performance of the devised algorithm, imaging experiments were conducted on a residential building at a distance of 1.531 km under strong sunlight. The target scene is shown in Figure 19a. To obtain the ideal target depth image, the same target area was imaged and detected at night using the peak-picking method with 5000 frames to accumulate the data. The resulting image was used as the ideal target depth image of this scene, as shown in Figure 19b. The SBR of the daytime imaging experiment was calculated to be 0.1 by dividing the average number of target photons detected by each pixel within the gating window by the total number of photons within the gating window. In a case involving 200 frames, the target depth image was reconstructed using the peak-picking method and the devised spatial-domain differential peak-picking method. The results are shown in Figure 20. The SBR of the daytime imaging experiment was calculated to be 0.1 by dividing the average number of target photons detected by each pixel within the gating window by the total number of photons within the gating window. In a case involving 200 frames, the target depth image was reconstructed using the peak-picking method and the devised spatial-domain differential peak-picking method. The results are shown in Figure 20. The SBR of the daytime imaging experiment was calculated to be 0.1 by dividing the average number of target photons detected by each pixel within the gating window by the total number of photons within the gating window. In a case involving 200 frames, the target depth image was reconstructed using the peak-picking method and the devised spatial-domain differential peak-picking method. The results are shown in Figure 20. The devised algorithm obtained a clear contour image of the target. The quality of the reconstructed depth image was quantitatively evaluated using various metrics, as summarized in Table 5. The devised method outperforms the classical peak-picking method, with improvements of 188%, 23.6%, and 87.9% in the K, PSNR, and SSIM values, respectively. The devised algorithm obtained a clear contour image of the target. The quality of the reconstructed depth image was quantitatively evaluated using various metrics, as summarized in Table 5. Table 5. Evaluation metrics for the peak-picking method and spatial-domain differential peakpicking method. The devised method outperforms the classical peak-picking method, with improvements of 188%, 23.6%, and 87.9% in the K, PSNR, and SSIM values, respectively. The depth image obtained by the spatial-domain differential peak-picking method was denoised through both TV recovery and FOTV recovery. Figure 21 shows the denoised results. The depth image obtained by the spatial-domain differential peak-picking method was denoised through both TV recovery and FOTV recovery. Figure 21 shows the denoised results. The devised FOTV-recovery algorithm restored more of the noise present in the target depth image and produced a smoother target area than the TV-recovery approach. The quality of the reconstructed depth image was quantitatively evaluated using various metrics, as summarized in Table 6.  The devised FOTV-recovery algorithm restored more of the noise present in the target depth image and produced a smoother target area than the TV-recovery approach. The quality of the reconstructed depth image was quantitatively evaluated using various metrics, as summarized in Table 6. Compared with the spatial-domain differential peak picking method, although the TV-recovery method improves the SSIM by 0.28%, K decreases by 23.7%. This decrease in K is attributable to the fact that in TV recovery, the noise points in the depth image are corrected, but the true points of the targets are also affected by the noise and thus deviate from the true depth values. In contrast, the devised FOTV-recovery method achieves improvements of 34.6%, 3.5%, and 7.2% in the K, PSNR, and SSIM values, respectively. The results highlight the potential of the devised FOTV-recovery method for effectively recovering GM-APD lidar depth images.

Discussion
This paper devised a FOTV-based method to restore GM-APD lidar depth images in low SBR and limited statistical frame conditions. First, the target depth image was extracted, and then, FOTV-based recovery was performed. The simulation results show that compared with the peak-picking method, the devised spatial-domain differential peak-picking method significantly improved the K, PSNR, and SSIM metrics when the SBR was 0.1 and the number of statistical frames was 30. Both the FOTVand TV-recovery methods enhanced these metrics, but compared with TV recovery, FOTV recovery achieved improvements of 10.1%, 14.4%, and 2.99% in the K, PSNR, and SSIM metrics, respectively. The experimental results demonstrate the effectiveness of the devised method: when the SBR was 0.1 and the number of statistical frames was 200, the K, PSNR, and SSIM values of the spatial-domain differential peak-picking method were 2.88, 1.236, and 1.87 times better than those of the peak-picking method, respectively. Combined with FOTV recovery, the devised method achieved improvements of 76.6%, 3.5%, and 6.9% in the K, PSNR, and SSIM metrics, respectively, compared with those when it was combined with TV recovery. These results indicate that the devised FOTV regularization method effectively restores, denoises, and preserves the fidelity of GM-APD depth images under low SBR and limited statistical frame conditions. Based on the research presented in this paper, there are several potential directions for future research. First, the depth image-recovery algorithm can be improved further. The current method improves the quality of the GM-APD lidar depth image to some extent, but there is still some noise and some artifacts present. Improving the quality of depth images acquired from GM-APD lidar under low signal-to-background ratio conditions will be of great value. Secondly, research on optimizing the extraction method of the target depth image under the condition of a low signal-to-background ratio and few frames, such as combining motion compensation technology or using the prior knowledge of the scene, can further improve the image restoration effect. In addition, since the application scenarios of GM-APD lidar require real-time performance, improving the efficient data processing performance of the algorithm is also an important direction for future research.