Robust Noise Suppression Technique for a LADAR System via Eigenvalue-Based Adaptive Filtering

The laser detection and ranging system (LADAR) is widely used in various fields that require 3D measurement, detection, and modeling. In order to improve the system stability and ranging accuracy, it is necessary to obtain the complete waveform of pulses that contain target information. Due to the inevitable noise, there are distinct deviations between the actual and expected waveforms, so noise suppression is essential. To achieve the best effect, the filters’ parameters that are usually set as empirical values should be adaptively adjusted according to the different noise levels. Therefore, we propose a novel noise suppression method for the LADAR system via eigenvalue-based adaptive filtering. Firstly, an efficient noise level estimation method is developed. The distributions of the eigenvalues of the sample covariance matrix are analyzed statistically after one-dimensional echo data are transformed into matrix format. Based on the boundedness and asymptotic properties of the noise eigenvalue spectrum, an estimation method for noise variances in high dimensional settings is proposed. Secondly, based on the estimated noise level, an adaptive guided filtering algorithm is designed within the gradient domain. The optimized parameters of the guided filtering are set according to an estimated noise level. Through simulation analysis and testing experiments on echo waves, it is proven that our algorithm can suppress the noise reliably and has advantages over the existing relevant methods.


Introduction
LADAR is widely used in unmanned aerial vehicles (UAV), terrain mapping, robots, and automobile auxiliary driving. Recently, through acquiring the distance information between the receiver and the targets and through capturing the intensity signals of the echo, the LADAR system can accomplish four dimensional imaging for the targets or scenes [1]. The system emits a pulsed signal and receives the echo, which is a Gaussian pulse waveform, then a Gaussian voltage is obtained after circuit transformation [2]. Before transmitting this analog Gaussian signal to the processing system, it is necessary to use a high-speed analog-to-digital converter (ADC) to digitize the signal into discrete point data, which are related to the intensity values [1,3]. In order to obtain environment and target information, echoes need to be further analyzed to get the amplitude, pulse width, and integral intensity information [4]. All the analyzing methods are based on noise suppression, and a good noise suppression effect guarantees effective analysis [5,6].
The noise of LADAR signal is mainly Gaussian white noise, which is generally divided into internal and external noise, introduced by the avalanche photodiode (APD) in the process of photoelectric conversion [7,8]. Internal noise, which mainly consists of shot noise, thermal noise [9], generation-recombination noise, and flicker noise, comes from the electronic noise itself. The generation mechanism is the random fluctuation of voltage or current caused by the random movement of carriers [10][11][12]. External noise that comes from the external interference of the detector includes circuit noise and background light noise. In the actual photoelectric detection process, the noise of the amplification circuit is small and usually ignored [4]. Especially in the daytime, background light noise is usually much larger than the system's natural noise. Background light noise mainly comes from natural and artificial light sources, such as the radiation from the Sun, ground, stars, atmosphere, clouds, and other sources of incidence or reflection to the detector. It also includes the interfering signal caused by the atmosphere in the transmission process [13].
For Gaussian white noise in an image, the spatial domain and transform domain filtering methods are usually used to reduce the noise. After one echo has been digitized into a one-dimensional signal, the image filtering algorithm can be used to suppress the noise of the LADAR signal. The most widely-used noise suppression methods in image processing are Gaussian algorithm [14], the bilateral filtering algorithm [15,16], and the guided filtering algorithm [17,18].
Compared with image noise, there are two special characteristics of LADAR. Firstly, multi-echo superposition will result in multiple positive and negative gradient values. The waveform is usually integrated in the time domain to analyze the material, texture, roughness, and other information of the target. Therefore, peaks and valleys need to be protected. Secondly, the light conditions, temperature, and objects around the LADAR are constantly changing [19]. As a result, it needs dynamic filter parameters that are associated with noise to adapt to the environment.
To solve the aforementioned problems, we propose a noise suppression algorithm that can accurately estimate the noise variances in two dimensions. Based on the requirement of preserving waveform boundaries, the guided filtering algorithm is used to suppress the noise of echoes, and the actual echoes are taken as the guide images. As a novel contribution, we extract the characteristic values of noise signals after separating the subspaces between the signals and the noise based on random matrix theory (RMT) [20], then we set the regularization parameters and window radius in gradient domain guided filtering. The practical usefulness of our method is illustrated by experiments.

Related Works
The echoes of LADAR are converted into digital signals after sampling and digitization by ADC. Since each echo signal is processed by an ADC of the same frequency, the number of data points of every echo is the same. Taking all data points of one echo as a one-dimensional orderly array, it can be processed by image filtering.
Image filtering is a basic module of image preprocessing, which can remove noise from images on the premise of preserving image details and textures as much as possible. Many spatial image filtering algorithms have been used for one-dimensional signal filtering [17,21,22]. The key problem of image filtering is to get a good balance between effectively removing noise and preserving detail and texture.
Image filtering generally includes the following categories: spatial filtering and transform domain filtering. Spatial image filtering calculates the neighborhood with each pixel in the image as the center and then replaces the original pixel value with the results obtained. Linear spatial filtering mainly includes box filtering, mean filtering, and Gaussian filtering, and nonlinear spatial filtering mainly includes median filtering, bilateral image filtering (BIL), and guided image filtering (GIF).
The box filter is a quick summation of the pixel values in each window within a given sliding window [23]. In the mean filter, the normalized box filter, each pixel of the output images is the average of the corresponding pixels of the input image in the window [24]. The mean filtering algorithm is simple and fast, which has a good suppression effect on the periodic interference noise. However, it destroys the detailed part of the image, while suppressing noise, thus making the image blurred.
Gaussian filtering is very effective at suppressing noise that obeys a normal distribution and is widely used in image processing. It is a process of the weighted average of the whole image. The value of each pixel is obtained by the weighted average of itself and other pixel values in the neighborhood. The median filter replaces each pixel of the image with the median value of the neighborhood pixel [25]. It can protect the edge of the signal while filtering the noise, especially the pulse noise, but it will lose the texture in the uniform medium area. BIF, which has the characteristics of simple, non-iterative, and local processing, combines the spatial proximity and pixel value similarity of the image. It considers the spatial information and gray similarity at the same time and achieves the goal of preserving noise suppression [15]. The BIF is a Gaussian filter function based on the spatial distribution, so the pixels far away from the edge will not affect the pixel values on the edge much. However, because it saves high-frequency information, it has a poor effect on high-frequency noise in the color image, and pulse noise is retained as the edge. Many scholars have proposed a variety of improvement methods, such as multi-resolution bilateral filtering [26][27][28], kernel regression [29], cross bilateral filtering [30], etc. GIF is used to filter the input image through a guide image, so that the final output image is roughly similar to the initial image, but the texture part is similar to the guide image [17]. GIF can improve the computing speed while maintain the edge, so it is widely used in image enhancement, high dynamic range (HDR) compression, image matting, and haze removal [21,22].
Transform domain filtering transforms images from the spatial domain to other domains. Utilizing the properties of the transform domain, the transform coefficient is further processed to reduce the noise, and then, the image is reconstructed by inverse transformation, such as discrete cosine transform (DCT) [31], wavelet transform [32][33][34], and so on. However, due to its huge computational cost, the processed image cannot not be acquired in real time. To improve its computing speed, Stockwell proposed the discrete cosine S transform [35]. Wavelet transform was used for image noise suppression, which is a more widely-used image noise suppression method in the transform domain [36].
At present, wavelet filtering methods can be divided into the Bayesian method and non-Bayesian method. Among them, non-Bayesian methods can be divided into three types: modular maximum reconstruction filtering, wavelet threshold filtering, and the wavelet spatial correlation filtering. The mode-maximum noise suppression method mentioned by Mallat has a good noise suppression effect on both Gaussian white noise and impulse noise, but it has a large amount of computation, and the iteration is prone to instability [37]. The spatial correlation noise suppression method proposed by Xu selects the coefficients according to the strong and weak performance of the correlation between the signal and noise in each scale, and the noise filtering effect is very significant [36]. The spatial correlation noise suppression method is easy to implement, but the error is large. The wavelet threshold noise suppression method was first mentioned by Donoho [33]. The wavelet coefficient of the signal is larger than the noise, so selecting the appropriate critical threshold can remove the noise. Because threshold setting is the key to this method, many scholars devoted themselves to improving the algorithm by setting various thresholds. To overcome the visual illusion of interference caused by the critical sampling wavelet coefficient, the over-complete wavelet was used to improve this defect [38].
In combination with the characteristics of the spatial domain and the transform domain, scholars proposed the gradient domain guided image filter (GGIF) [21,39,40]. It can obtain better effects than GIF and weighted guided image filtering (WGIF) by using the gradient domain information in guiding the image as the weight.
However, noise would change the edge part of the image, and the impact is unpredictable, so the above algorithms suffer from the fixed noise parameter. To address this, we propose an adaptive gradient-guided filtering (AGGF) algorithm, which was inspired by [21,22], based on noise estimation. Its framework is shown as Figure 1.

Adaptive Gradient-Guided Noise Suppression Technique
In [21], there were a guidance image G and an image to be filtered X, and they could be identical. Let Ω δ 1 (p) be a square window centered at a pixel p of a radius δ 1 . It was assumed that the output imageẐ is a linear transform of the guidance image G in the window Ω δ 1 (p ): where a p and b p are two constants in the window Ω δ 1 (p ). Cost function E a p , b p with edge-aware weightingΓ G (p ) is as follows: ψ is a regularization parameter that restrains large a p . µ χ,∞ is the mean value of all χ (p). η is an adjustment coefficient that is defined as 4/ (µ χ,∞ − min (χ (p))); it makes a p less sensitive.
Γ G (p ) is an edge-aware weighting. The radius of the window is three and δ 1 for p and p, is a small positive constant. χ (p ) is defined as σ G,1 (p ) σ G,δ 1 (p ). It is usually set to 16 in detailed manipulation applications. The weightingΓ G (p ) measures the importance of pixel p with respect to the whole guidance image [21].
where is the element-wise product of two matrices. µ G X,δ 1 (p ), µ G,δ 1 (p ), and µ X,δ 1 (p ) are the mean values of G X, G, and X in the window Ω δ 1 (p ), respectively [22]. According to the influence of noise on the parameters, a new filtering algorithm is proposed.

Gradient-Guided Filtering
The same as the GGIF, the key assumption of the AGGF is a local linear model between the guidance image Gand the filtering outputZ (p). It is shown that Z (p) = a p G, and the smoothness ofZ (p) in Ω δ depends on a p . The optimal values of a p and b p are computed by minimizing a new cost function, which is defined as: where ψ n is the sensitivity coefficient of guided filtering and γ p is the same as Equation (3). The GGIF [21] is less sensitive to the selection of ψ, and edges could be preserved better than both the GIF and the WGIF. However, GGIF does not define where the edge is; when the overall gradient of the image does not change much, it will have a negative impact on the filtering effect, shown as Figure 2a.
To avoid this, we add the parameters α to modify γ p . When the maximum gradient is less than a threshold TH, α n is zero, and when the maximum gradient is greater than a threshold value, α n is one.
In other words, it uses γ p to modify ψ n when the overall gradient changes a little. TH is defined as: where δ is the window radius, which is defined as Equation (10). H is the peak value of the echo signal. A δ , B δ , and C δ are coefficients, which are 2.367, −0.82, and 0.286, respectively.σ 2 n is the noise level parameter that will be introduced in detail later. K δ is the normalized coefficients that can adjust δ as the pixel gradient changes. Grad max , Grad min , and Grad max are the maximum gradient value of all pixels, the minimum gradient value, and the gradient value of pixel p in the guided image, respectively. In addition, due to the box filter in [17], the complexity ofΓ G (p ) is O(N) for an image with N pixels. ψ n defined as Equation (12) can be adjusted according to the noise level to obtain a better filtering effect.
where A ψ n , B ψ n , and C ψ n are fitting coefficients, which are 5.5, 11, and 66, respectively.
The final value ofZ (p) is given as follows: whereā p andb p are the mean values of a p and b p in the window. |Ω δ (p)| is the cardinality of |Ω δ (p )|. When the pixel p is at an edge, the value ofΓ G (p ) is usually much larger than one, which means edges can be preserved. In addition, the complexity of the AGGF is O(N) for an image with N pixels, which is the same as that of the GIF [17].
In order to illustrate the influence of threshold TH on the algorithm, we designed a comparison scheme. The maximum gradient in Figure 2a is less than the threshold; γ p is zero. At this time, the filtering effect is similar to WGIF, but it is not exactly the same. In GGIF, γ p is one, which causes excessive correction of a p , and the filtering effect is very different.
In Figure 2b, the maximum gradient is larger than the threshold. the filtering effect is similar to GGIF, when the window is not adjusted adaptively. It should be noted that in the LADAR system, echo is expected by Gaussian noise, which rarely exceeds the threshold, so γ p needs to be corrected.

Noise Level Estimation via Eigenvalue Analysis
In this section, we analyze the asymptotical distribution of the ratio of extreme eigenvalues of a sample covariance matrix based on the limiting RTMlaw, then the noise level is estimated. To derive the relationship between the eigenvalues and noise level, we first construct the sample covariance matrix Σ.
Each element d i in the matrix Σ defined as Equation (17) is a data point obtained by sampling the echo. N is the number of columns in the matrix; S is the number of elements in each column. All data points of one echo is a column. N echos make the matrix Σ, and dimension of the matrix Σ is S. This matrix can be equivalent to an image to be processed, each data point being a pixel of the image. N is greater than S. The pulse width of the commonly-used pulsed laser is 5 ns, and the effective wavelength of an echo can be approximately 10 ns. The sampling frequency of the commonly-used ADC is 2 GHz, and an effective echo can obtain 20 sampling points. Therefore, in this work, N is at least 20.
According to the symmetric property, this matrix is decomposed into the product of three matrices: an orthogonal matrix U, matrix Σ, and a transpose matrix U T , which can be selected by satisfying U T U = I. Here, this transform process is written as: Given λ 1 ≥ λ 2 ≥ · · · ≥ λ S , we exploit the sequence of eigenvalues to enable the separation of signals from noise. To be more specific, we divide the eigenvalues into two sets S 0 = S 1 ∪ S 2 by finding the appropriate bound in a spiked population model [41]. The echo signal lies in low-dimensional subspace, and thus, the leading eigenvalues in set is dominated by the noise. Because the signals contribute very little to this latter portion, we take all the eigenvalues of S 2 into consideration to estimate the noise variance while eliminating the influence of the echo signal. Moreover, the random variables {λ i } S i=m+1 can be considered as the eigenvalues of the pure noise covariance matrix. Suppose the sample matrix Σ has the form (N − m) Σ = HH T , where the sample entries of H are independently generated from the distribution N 0,σ 2 n . Then, the real matrix M = HH T follows a standard Wishart distribution [42]. In the high dimensional situation: N/S → γ ∈ [0, ∞) as S, N → ∞, the Tracy-Widom law gives the limiting distribution of the largest eigenvalue of the large random matrix M [43]. Then, we have the following asymptotic expression: where F TW1 (Z) indicates the cumulative distribution function with respect to the Tracy-Widom random variable. In order to improve both the approximation accuracy and convergence rate, even only with a few signal samples, we need to choose the suitable centering and scaling parameters µ, ξ [44]. By the comparison between different values, such parameters are defined as: The empirical distribution of the eigenvalues of the large sample matrix converges almost surely to the Marchenko-Pastur distribution on a finite support. Based on the generalized result in [45], when N → ∞ and γ ∈ [0, ∞), with probability one, we derive the limiting value of the smallest eigenvalue as: According to the asymptotic distributions described in Equations (19) and (21), we further quantify the distribution of the ratio of the maximum eigenvalue to the minimum eigenvalue in order to detect the noise eigenvalues. Let T n be a detection threshold, then we find T n by the following expression: Note that there is no closed-form expression for the function F TW1 . Fortunately, the values of F TW1 and the inverse F −1 TW1 can be numerically computed at certain percentile points [41]. For a required detection probability α 1 , this leads to: Plugging the definitions of µ and ξ into Equation (23), we finally obtain the threshold.
When the detection threshold T n is known in the given probability, it means that an asymptotic upper bound can also be obtained for determining the noise eigenvalues of the matrix Σ. In general, the noise eigenvalues in the set S 2 surround the true noise variance as it follows the Gaussian distribution. The estimated largest eigenvalue λ m+1 should be no less thanσ 2 n . The known smallest eigenvalue λ S is no more thanσ 2 n by the theoretical analysis [46]. The location and value of λ m+1 in S are obtained by checking the bound λ m+1 ≤ T n · λ S with high probability α 1 .
Without requiring the knowledge of the signal, the threshold T n can provide good detection performance for finite S, N, even when the ratio N/S is not too large. Based on this result, more accurate estimation can be obtained by averaging all elements in S 2 . Hence, the maximum likelihood estimator ofσ 2 n is: Based on the above content, we propose a noise suppression algorithm as Algorithm 1.

Algorithm 1 Noise suppression algorithm
1: procedure (Input: X, N, S, α 1 , α 2 , f p ) X is input waveform; N is the number of echoes; S is the number of data points in each echo; α 1 , α 2 are probability levels; f p is frequency of the ADC.

Experimental Results
We studied the influence of the noise on filter parameters through a large number of data simulations. Experiments on land and simulation under water were carried out to verify the effectiveness of the algorithm. As shown in Figure 3a, all measured data in this paper were obtained by the LADAR system, which used an 8 × 8 array APDas the detector and a 1064-nm pulsed laser as the detection source to realize target distance detection based on the time-of-flight method. The system uses a multi-channel preamplifier circuit based on the variable gain amplifier (VGA) and implements full waveform sampling analysis based on a 2-GHz sampling rate ADC.

Evaluation
As shown in Figure 4a, the maximum deviation of estimation was 1.2%, while the SNR was in the range of 0-35. When the SNR is low, the eigenvalue of noise in signal matrix is large, which can be extracted more accurately, and the estimation is better. In Figure 4b, it can be seen that a good noise suppression effect can be achieved, and the fittings were more than 99% at different SNRs. Table 1 shows the comparison of the filtering effects of these algorithms with different sensitivities at different SNR. The optimal sensitivity of WGIF, GGIF, AGGF, and GIF were approximately equal at the same SNR, and this can be obtained by comparing Equations (2) and (7). In other words, all of them optimized the edge filtering effect by weighting the ψ n . Table 1. Influence of SNR on ψ n for different filtering algorithms, f p = 5G. SNR = 35, ψ n = 69, δ = 4 in WGIF, GGIF and GIF, δ = 3-5 in AGGF; SNR = 30, ψ n = 134, δ = 6 in WGIF, GGIF and GIF, δ = 3-7 in AGGF; SNR = 25, ψ n = 295, δ = 7 in WGIF, GGIF and GIF, δ = 4-8 in AGGF; SNR = 20, ψ n = 1400, δ = 8 in WGIF, GGIF and GIF, δ = 6-10 in AGGF; SNR = 15, ψ n = 5800, δ = 9 in WGIF, GGIF and GIF, δ = 8-12 in AGGF; SNR = 10, ψ n = 1,270,000, δ = 20 in WGIF, GGIF and GIF, δ = 12-16 in AGGF. The best results have been bolded.  Figure 5a shows the relationship between the optimal ψ n and SNR of the AGGF algorithm. Due to ψ n changing with SNR, the proposed AGGF algorithm can achieve a good noise suppression effect at different noise levels. It can be obtained that the larger the SNR, the smaller the optimal sensitivity. Figure 5b shows the relationship between ψ n , window width (2δ+1), f p , and SNR.   For the LADAR waveform, the sampling rate of the ADC determines the number of data points of one echo, which is directly related to the window radius. The higher the sampling rate, the larger the window width should be.

Analysis of the Relationship Between the Noise and Parameters
At the same time, in order to achieve better filtering, the noise level should also be considered when setting the window width.
The proposed algorithm can adjust filter parameters according to different noise levels. To make the comparison fairer, all window radii were optimized by Equation (10). It can be observed that the effect of indoor experiment was not as good as that of the underwater experiment, and the detailed reasons will be analyzed later.
Although the filtering effect can be reflected by the two parameters of mean variance and fitting degree, they cannot fully reflect the filtering effect. For example, at the sampling rate of 2.5 GHz, WGIF and GIF had the best parameters when δ was three and five, respectively.
The unsmoothed waveform needs to be processed again to obtain the centroid position, which was used for the accurate calculation of the target distance. In contrast, the proposed adaptive window width achieved a good filtering effect at different sampling rates and SNR. Specifically, GGIF can enhance the boundary effect with parameter γ p , but in the case of Gaussian echo, where the boundary is not obvious or the overall gradient change is not dramatic, excessive boundary judgment will occur.
In order to obtain a better filtering effect, we make the window radius larger in the part with a high slope and smaller in the part with a low slope, as Figure 3b shows. Figure 6a,b is the waveforms measured by the LADAR system on land, in which (a) is about one object and (b) is the waveform measured by two objects 30 cm apart. Since there was no obvious gradient change, the noise suppression effect of GGIF was poor, and the effect of the other algorithms was similar. Figures 7 and 8 show the underwater waveforms simulated at 15 m and 25 m, respectively. The two echoes before and after were the surface signal and underwater signal, separately. Because the light beam reflected continuously in the water, the first echo had trailing.

Analysis of Measurement and Simulation Data
It can be seen in Section 3 that WGIF and GIF had a good filtering effect on waveforms with a small gradient, and they can maintain edge features well at the same time, but the edge protection effect was poor in waveforms with a large gradient.
GGIF enhanced the edge detection based on WGIF and GIF, and it was more sensitive to the edge. When the gradient was small, there would be excessive protection of the edge, which can be seen in Figure 2. We hoped that the algorithm can better distinguish the two cases, so we increased the threshold to judge the gradient change. In the case of a large gradient, the proposed algorithm approximated the GGIF, as shown in the Figures 7 and 8. In the case of small gradient, the proposed algorithm approximated the WGIF and GIF, as shown in the Figure 6.     In the seabed detection, due to the change of the seabed height covered by the divergence of the beam and the scattering effect of sea water, the echo gradient will change to different degrees. When the seawater is shallower, the laser energy absorbed by seawater is less, and the echo power is larger, which leads to a big gradient of echoes. The gradient of the echo signal varies with the fluctuation of the seabed. Therefore, the proposed algorithm can achieve better results in different conditions for seabed echoes. In the laboratory experiment, considering that the echo power is too large to produce a saturated signal, it is impossible to carry out effective analysis. Therefore, in the indoor test, it is necessary to add an attenuation lens to the laser transmitting lens to reduce the laser power and at the same time reduce the gain of the signal processing circuit, which will reduce the echo gradient. Therefore, in the comparison of the indoor signal filtering effect, the proposed algorithm has no obvious advantages.
As shown in Figure 8, when the underwater echo is very strong, the advantage of the proposed algorithm is more obvious, as a result of adjusting parameters reasonably according to the noise level and gradient. As shown in Figure 7, when the underwater echo is very small, the AGGF noise suppression effect is still better. Therefore, this proposed method can be better applied to noise suppression of the LADAR waveform.
Detailed simulation results are shown in Figures 9-12. It can be seen that the algorithm proposed is more advantageous in the application of a large gradient, such as underwater.

Conclusions
To infer the noise level from a signal matrix, the eigen-spaces of the signal and noise are transformed and separated well by determining the eigen-spectrum interval. In addition, the developed estimator can effectively eliminate the estimation bias of a noise variance in a high-dimensional context. The experimental results have demonstrated that the proposed noise suppression technique for the LADAR system via eigenvalue-based adaptive filtering can adjust parameters according to the noise level, and it outperforms the relevant existing methods over a wide range of noise level conditions for LADAR system. Author Contributions: Conceptualization, X.X. and R.C.; Data curation, P.W.; Formal analysis, X.X.; Investigation, Y.Z.; Methodology, X.X.; Project administration, X.X. and Y.Z.; Resources, X.X.; Software, X.X. and P.W.; Writing-original draft, X.X.; Writing-review & editing, X.X., R.C. and Y.Z.