A Method for Depth Estimation of Long-Range LiDAR Based on Intensity Guidance

Abstract

The present invention discloses a method for depth estimation of a long-distance lidar based on intensity guidance, including: in data preprocessing, fully utilizing the spatio-temporal correlation and intensity information of echo signals, and alternately establishing multi-scale superpixels and performing fast time-domain windowing on the lidar three-dimensional point cloud data; adopting an optimization framework, based on the Poisson distribution model of the preprocessed lidar three-dimensional point cloud data, establishing a depth estimation cost function that combines a data fidelity term and a regularization term introduced by intensity information; using the alternating direction multiplier algorithm to estimate the depth image from the cost function. The present invention addresses the problems of sparse and large-volume three-dimensional point clouds obtained by long-distance lidars, estimates the depth image with smaller memory requirements and lower computational complexity, and effectively improves the quality of the depth image and enhances the robustness to low photon levels.

Description

Technical Field

[0001] The present invention belongs to the technical field of lidar depth estimation, and specifically relates to a method for lidar depth estimation at a long distance based on intensity guidance. Background Art

[0002] Long-distance three-dimensional imaging has been widely applied in fields such as topographic surveying, underwater exploration, and three-dimensional remote sensing. The lidar system uses an avalanche photodiode (GM-APD) operating in Geiger mode, which has single-photon sensitivity and picosecond-level time resolution. The photon technology statistical histogram of each pixel is obtained by time-correlated single-photon counting (TCSPC) technology, where the position and quantity of the counts respectively correspond to depth and intensity information, thereby reconstructing the 3D scene.

[0003] The continuous development of long-distance photon-counting lidar enables it to work at a farther distance, with a larger field of view, and higher resolution requirements. This means that the size of the lidar three-dimensional data obtained is larger, and the corresponding memory requirements and computational complexity are greater. At the same time, in the outdoor scenes where long-distance lidar is applied, due to the existence of huge attenuation, the signals returned by the targets are extremely weak and are often overwhelmed by strong background noise. Therefore, accurate lidar depth estimation is challenging.

[0004] Currently, in addition to optimizing lidar hardware systems, there already exist some back-end processing algorithms to improve the depth estimation quality and robustness under low photon levels. Halimi et al. utilized non-local spatial correlations to improve the robustness of the Alternating Direction Method of Multipliers (ADMM) for estimating depth images. Although this method adopts a non-uniform sampling strategy to reduce computational costs, it only estimates from 2D images while ignoring the influence of the 3D point cloud's time dimension (1. Chen S, Halimi A, Ren X, et al. Learning Non-Local Spatial Correlations To Restore Sparse 3D Single-Photon Data[J]. IEEE Transactions on Image Processing, 2019, PP(99).). Ferstl et al. fused the high-resolution RGB map from the camera with the low-resolution depth map from the lidar to generate a dense depth image, but the registration between the camera and the lidar is complex, and for long-distance imaging, the camera is expensive and large in size (2. Ferstl D, Reinbacher C, Ranftl R, et al. Image Guided Depth Upsampling Using Anisotropic Total Generalized Variation[C] / / Proceedings of the 2013 IEEE International Conference on Computer Vision. IEEE, 2013.). Tachella et al. combined the reversible jump Markov chain Monte Carlo method with a spatial point process to accurately reconstruct the 3D point cloud in the case of multiple surfaces existing at each pixel, but the computational time is too long (3. Tachella J, Altmann Y, X Ren, et al. Bayesian 3D Reconstruction of Complex Scenes from Single-Photon Lidar Data[J]. Siam Journal on Imaging Sciences, 2018, 12(1):521-550.). Summary of the Invention

[0005] The object of the present invention is to provide a method for estimating the depth of a long-distance lidar based on intensity guidance, so as to solve the problem that the depth map estimated from the three-dimensional point cloud of a long-distance lidar with a large amount of data and extremely sparse data has a large error, thereby estimating a more accurate depth image with smaller memory requirements and computational complexity.

[0006] The technical solution for achieving the object of the present invention is: a method for estimating the depth of a long-distance lidar based on intensity guidance, and the specific steps are as follows:

[0007] Step 1: Utilize the spatio-temporal correlation of intensity information and echo signals to alternately perform the establishment of multi-scale superpixels and fast time-domain windowing on the three-dimensional lidar point cloud data to complete data preprocessing;

[0008] Step 2: Adopt an optimization framework, and based on the Poisson distribution model of the preprocessed three-dimensional lidar point cloud data, establish a depth estimation cost function that combines a data term and a regularization term introduced by intensity information;

[0009] Step 3: Use the alternating direction multiplier algorithm to estimate the depth image from the cost function.

[0010] Preferably, the specific method for alternately performing the establishment of multi-scale superpixels and fast time-domain windowing on the three-dimensional lidar point cloud data is as follows:

[0011] Step 1.1: Let the current stage be l, 1 ≤ l ≤ l max , where l max is the preset maximum number of stages;

[0012] Step 1.2: Use the Canny operator for edge detection on the intensity map r l-1 to obtain a binary matrix D l that marks the edges, where D l (i,j) is an element in the matrix D l ; and use the matrix F l to mark the empty pixels and noise pixels in the depth map t l-1 , where F l (i,j) is an element in the matrix F l ;

[0013] Step 1.3: Perform the establishment of superpixels. For each pixel (i,j), the superpixel N i,j The specific formula is:

[0014] N i,j ={(x,y)∈{i-l...,i+l}×{j-l...,j+l}:

[0015] F l-1(x, y) = 0, D l-1 (x, y) = D l-1 (i, j),

[0016]

[0017] Wherein, is the intensity difference between pixel (x, y) and pixel (i, j), and R is the threshold;

[0018] For pixel (x, y) ∈ N i,j , which is on the same edge as pixel (i, j) and is neither an empty pixel nor a noise pixel. When the intensity difference is less than the threshold R, they are merged into a new superpixel n. At this time, the photon counting histogram is expressed as Regarding these pixels as the new superpixel n, T is the total number of time units;

[0019] For pixel (x, y) ∈ {i - l..., i + l} × {j - l..., j + l} and F l-1 (x, y) = 0, but D l-1 (x, y) ≠ D l -1 (i, j) or Retain the superpixel at the original scale;

[0020] For pixel (x, y) ∈ {i - l..., i + l} × {j - l…, j + l} and F l-1 (x, y) = 1, they are empty pixels or noise pixels, and the photon counting histogram is expressed as The quotient of the number of pixels in the superpixel n;

[0021] Step 1.4: Perform fast time-domain windowing on the histogram with the number of time units T for each pixel (i, j) After windowing, the data is obtained Wherein, is the photon count of the t-th time unit on pixel (i, j) after being processed in stage l, and T m l (i, j) is the first time unit of the windowed interval, and T w is the preset number of time units;

[0022] Step 1.5: Use maximum likelihood estimation to obtain the intensity map r l and the depth map t l , and the specific formula is:

[0023]

[0024]

[0025] In the formula, is the photon count at the t-th time unit after processing in the l-th stage for the pixel (i, j), and N r is the number of row pixels, and N c is the number of column pixels;

[0026] Step 1.6: If the decision condition is satisfied: If there are no empty pixels and noise pixels in the current stage l, that is, F l (i, j) = 0, or the preset maximum number of stages l max is reached, the processed lidar point cloud data is obtained. Among them, is the photon count at the t-th time unit for the superpixel n, and N superpixel is the total number of superpixels, T m (n) is the first time unit of the windowed interval for the superpixel n. Otherwise, l = l + 1, and steps 1.1 to 1.6 are repeated.

[0027] Preferably, the specific method for constructing the depth estimation cost function is as follows:

[0028] Based on the Poisson distribution model of the point cloud data obtained by preprocessing, the data term L of the cost function is constructed, and the specific formula is:

[0029]

[0030] In the formula, g(t) is the lidar system impulse response function, is the photon count in the t-th time interval for the superpixel n, t ∈ {T m (n),..., T m (n) + T w}}, n ∈ {1,..., N superpixel}}, where T m (n) is the first time unit of the windowed interval for the superpixel n, T w is the number of time units, N superpixel is the number of superpixels, and t n is the depth value to be solved for the superpixel n;

[0031] The edge amplitude matrix M of the intensity map is obtained by using the Canny operator, and the edge amplitude matrix M is introduced as an indicator function into the total variation regularization model TV to establish the depth estimation cost function C combined with the data term, specifically:

[0032]

[0033] where λ is the regularization parameter, and the total variation model TV(t) = ||Dt|| 1 = ||D r t|| 1 + ||D c t|| 1 , where D r and D c represent the first-order differences in the horizontal and vertical directions respectively, M is the edge indicator function, α is a constant greater than 0, and t is the depth map to be solved.

[0034] Compared with the prior art, the present invention has the following remarkable advantages: (1) Both the data preprocessing and the depth image estimation fully utilize the intensity information. The intensity information and the depth information are both from the lidar three-dimensional point cloud data, having spatial consistency. Moreover, compared with the edges in the depth map being only determined by the target contour, the edges in the intensity map are affected by various factors such as the target contour, surface texture, and different illuminations, so it has more significant edge and texture features. (2) In data preprocessing, the establishment of multi-scale superpixels and pixel-by-pixel time-domain windowing are alternately performed. On the one hand, the method for establishing intensity-guided multi-scale superpixels has adaptability to natural scenes, improving the robustness of the present invention to low photon levels and minimizing the image size as much as possible while retaining target details. On the other hand, the pixel-by-pixel time-domain windowing method separates signals from noise and, inspired by the dichotomy method, uses a fast search method to accelerate its operation speed, which can greatly reduce the length of the time dimension of the point cloud data. In summary, the present invention significantly reduces the data size of the lidar three-dimensional point cloud, correspondingly reducing the memory requirements and computational complexity of the subsequent depth estimation algorithm, and filtering out noise pixels or empty pixels. (3) In depth image estimation, the present invention introduces an edge indicator function on the basis of TV regularization, and uses ADMM to estimate the depth map from the improved cost function model, effectively preserving the details of the edge part and filtering out the noise in the smooth area.

[0035] The present invention will be further described in detail below with reference to the accompanying drawings. Description of the Drawings

[0036] Figure 1 is an experimental target scene diagram for verifying the present invention.

[0037] Figure 2 is the target reference depth map and three-dimensional point cloud map obtained when the sampling integration time is 50 ms.

[0038] Figure 3 is the target depth map and three-dimensional point cloud map obtained by the traditional imaging method based on the maximum likelihood estimation algorithm when the sampling integration time is 10 ms.

[0039] Figure 4The target depth map and 3D point cloud map are obtained by the traditional imaging method based on the maximum likelihood estimation algorithm when the sampling integration time is 0.5 ms.

[0040] Figure 5 It is the flow chart of the depth estimation method for long-distance lidar based on intensity guidance of the present invention.

[0041] Figure 6 It is the example diagram of the fast time-domain windowing of the present invention.

[0042] Figure 7 The target depth map and 3D point cloud map are obtained by the depth estimation method for long-distance lidar based on intensity guidance of the present invention when the sampling integration time is 10 ms.

[0043] Figure 8 The target depth map and 3D point cloud map are obtained by the depth estimation method for long-distance lidar based on intensity guidance of the present invention when the sampling integration time is 0.5 ms. Detailed implementation manners

[0044] As Figure 1 shown, it is the schematic diagram of the target scene for verifying a depth estimation method for long-distance lidar based on intensity guidance of the present invention. The target scene is about 1.4 km away from the lidar system. Centering on the peak of the histogram, the data with a total length of 180 m before and after is retained, so as to obtain the original 3D point cloud with a size of 200×200×3000.

[0045] As Figure 2 shown, it is the target reference depth map and 3D point cloud map obtained when the sampling integration time is 50 ms. Figure 2 (a) Depth map, Figure 2 (b) is the 3D point cloud map.

[0046] As Figure 3 shown, it is the target depth map and 3D point cloud map obtained by the traditional imaging method based on the maximum likelihood estimation algorithm when the sampling integration time is 10 ms. Figure 3 (a) Depth map, Figure 3 (b) is the 3D point cloud map.

[0047] As Figure 4 shown, it is the target depth map and 3D point cloud map obtained by the traditional imaging method based on the maximum likelihood estimation algorithm when the sampling integration time is 0.5 ms. Figure 4 (a) Depth map, Figure 4 (b) is the 3D point cloud map.

[0048] Combined with Figure 5, a method for long-range LiDAR depth estimation based on intensity guidance, firstly, preprocess the LiDAR 3D point cloud data, on the one hand, perform pixel-by-pixel time domain windowing, and use the difference in the probability distribution characteristics of noise response and signal response: the signal response is mainly concentrated in the narrow half-maximum full width of the laser pulse, while the noise response is evenly distributed in the entire time dimension and is independent of the laser pulse, so it can be assumed to be a constant. According to the time correlation of the above signal response, a fixed time interval T is obtained on the photon count histogram of each pixel. w The set of responses with the largest sum of photon numbers is considered to be the set of signal responses. Inspired by the dichotomy method, a fast search method is proposed to optimize pixel-by-pixel time domain windowing. The main principle is that two-thirds of the values ​​can be discarded each time to narrow the range: for each pixel, in each search, the time dimension is divided into three intervals, and the sum of the photon count values ​​of the signal response in the three intervals is compared, and the interval with the largest sum is retained as the next search range. Then perform the same search in the new search range. Repeat this process until the length of the new search range is less than or equal to T. w The response set obtained by each pixel is judged. If the sum of the photon counts is less than the preset threshold, it is judged as a noise response and all responses of the pixel are ignored. Otherwise, the judgment length is T w The response within the interval is the signal response. On the other hand, assuming that the edge of the depth map corresponds to the edge of the intensity map, an intensity-guided multi-scale superpixel is established: pixels with small intensity differences are in the target smooth area, and they are merged into larger-scale (lower resolution) superpixels to reduce the image size; while pixels with large intensity differences are in the target edge part, and the original smaller-scale superpixels are retained, thereby maintaining details and clear edges. In particular, for empty pixels and noisy pixels located at the edge of the target, only the values ​​of larger-scale superpixels are used to repair these pixels. In summary, the establishment of multi-scale superpixels and fast time-domain windowing are performed alternately until the superpixel reaches the set maximum scale, or there are no empty pixels or noise pixels at the current scale.

[0049] Secondly, an optimization framework was adopted to establish a cost function, which transformed the depth estimation problem into an optimal solution problem: the data terms of the cost function were constructed according to the Poisson distribution model of the preprocessed lidar 3D point cloud data, and then the edge information of the intensity map was introduced to form a weighted TV regularization model.

[0050] Finally, the ADMM algorithm is used to estimate the depth image from the cost function: the optimization problem constructed in the second step is a convex function with separated functions and variables. Therefore, the fast-converging ADMM algorithm can decompose the problem into several independent parts, thereby iteratively solving the optimal depth image.

[0051] Using the above three steps, an accurate depth image can be estimated.

[0052] The specific implementation method is as follows:

[0053] Step 1: In data preprocessing, fully utilize the intensity information and spatio-temporal correlation of the echo signal, and alternately perform the establishment of multi-scale superpixels and fast time-domain windowing on the original three-dimensional point cloud s with a size of N r ×N c ×T, where N = N r ×N r is the total number of pixels, N r is the number of row pixels, N c is the number of column pixels, and T is the total number of time units.

[0054] The basic principle of establishing multi-scale superpixels is as follows: Both intensity information and depth information come from the three-dimensional point cloud data of the lidar, and they have spatial consistency. Therefore, the edges of the depth map and the intensity map correspond to each other. Moreover, compared with the edges in the depth map that are only determined by the target contour, the edges of the intensity map are affected by various factors such as the target contour, surface texture, and different illuminations. Therefore, they have more prominent edge and texture features.

[0055] The basic principle of time-domain windowing lies in the different probability distribution characteristics of the noise response and the signal response: The signal response is mainly concentrated within the relatively narrow full width at half maximum of the laser pulse, while the noise response is uniformly distributed over the entire time dimension and is independent of the laser pulse, so it can be assumed to be a constant. According to the above time correlation of the signal response, obtain the set of responses with the largest sum of the number of photons within the fixed time interval T w in the photon count histogram of each pixel, and consider it as the set of signal responses. Inspired by the dichotomy method, a fast search method is proposed to optimize the per-pixel time-domain windowing. The main principle is that two-thirds of the values can be discarded in each search to narrow the range. For the three-dimensional point cloud data with a size of N r ×N c ×T, compared with windowing by sliding an interval with a length of T w over T time units, the computational complexity of the algorithm is reduced from O(N r ×N c ×T) to O(N r ×N c ×log(T)).

[0056] Based on the above principles, alternately perform the establishment of multi-scale superpixels and fast time-domain windowing. The specific steps are as follows:

[0057] (1) Let the current stage be l (1 ≤ l ≤ l max ), where l max is the preset maximum number of stages. The point cloud data intensity map rl-1 and the depth map t l-1 as input data.

[0058] Specifically, when l = 1, is the raw data without processing, where s i,j,t is the photon count of the t-th time unit at the pixel (i, j), N = N r × N r is the total number of pixels, and T is the total number of time units. At this time, the intensity map r l-1 and the depth map t l-1 are expressed as:

[0059]

[0060]

[0061] In the formula, g(t) is the pulse response function of the lidar system.

[0062] (2) Apply the Canny operator for edge detection to the intensity map r l-1 to obtain the binary matrix D that marks the edges l , where and use the matrix F l to mark the empty pixels and noise pixels in the depth map t l-1 (if the intensity value of a pixel is greater than the average value of its 8-neighborhood pixels, then this pixel is considered noise), where

[0063] (3) Establish superpixels. At this stage l, the size of the superpixels is (2l + 1) × (2l + 1). For each pixel (i, j), the superpixel N i,j The specific formula is:

[0064] N i,j = {(x, y) ∈ {i - l..., i + l} × {j - l..., j + l}:

[0065] F l-1 (x, y) = 0, D l-1 (x, y) = D l-1 (i, j),

[0066]

[0067] In the formula, R is the preset intensity threshold.

[0068] For the pixels (x, y) ∈ N i,j , they are on the same edge as the pixel (i, j) and are neither empty pixels nor noise pixels. When the intensity difference If less than the threshold R, these pixels are considered to be in the smooth region of the target, and they are merged into a new superpixel n of a larger scale, thereby reducing the image size. At this time, the photon counting histogram is expressed as Regard these pixels as the new superpixel n,

[0069] For pixels (x, y) ∈ {i - l…, i + l} × {j - l…, j + l} and F l-1 (x, y) = 0, but D l-1 (x, y) ≠ D l-1 (i, j) or These pixels are considered to be in the edge part of the target, and the original smaller-scale superpixels are retained to maintain details and clear edges.

[0070] For pixels (x, y) ∈ {i - l..., i + l} × {j - l..., j + l} and F l-1 (x, y) = 1, they are empty pixels or noise pixels, and the photon counting histogram is expressed as The quotient of the number of pixels in the superpixel n, preventing the destruction of the characteristics of the target.

[0071] (4) Perform the fast time-domain windowing method on the histogram with the number of time units T for each pixel (i, j) After windowing, the data is obtained Among them, Is the photon count of the t-th time unit on the pixel (i, j) after being processed in stage l, T m l (i, j) is the first time unit of the windowed interval, T w Is the preset number of time units. The specific method of fast time-domain windowing is:

[0072] Define five variables for one search, which are respectively And:

[0073]

[0074]

[0075]

[0076] At this time, the histogram is divided into three subsets of equal length, and the sum of their photon counts is calculated: And As Figure 6 Shown, the histogram is divided into three subsets, so that even if the center of the signal response Or Happens to be at the edge of a subset, there is always another subset that completely covers this response.

[0077] Compare Retain the subset with the largest sum of photon counts And update the search value

[0078]

[0079]

[0080] Thus Used to index the first time unit of the set after windowing

[0081] For the length Perform the next search on the subset within, and calculate using the above formula And Generate three new subsets. Retain the subset with the most photon counts in each search. Stop until the condition is met, then the search stops. Obtain the set with the most photon counts on pixel (i, j)

[0082] Set the threshold K according to the noise level of the histogram for separating signal and noise: If the sum of photon counts of the set is less than the threshold K, it is determined as a noise response and all counts of this pixel are ignored; otherwise, it is determined as a signal response and retain the photon counts with an interval length of T w

[0083] Finally, after the above steps, the response set obtained on pixel (i, j) is

[0084]

[0085] (5) Obtain the intensity map r l and the depth map t l using maximum likelihood estimation, and the specific formulas are

[0086]

[0087]

[0088] In the formula is the photon count of the t-th time unit after stage processing l on pixel (i, j).

[0089] (6) If the decision condition is met: there are no empty pixels and noise pixels in the current stage l, that is, F l (i, j) = 0 or reach the preset maximum number of stages l max ​, the processed lidar point cloud data is obtained Among them, is the photon count of the t-th time unit on superpixel n, and N superpixel is the total number of superpixels, and T m (n) is the first time unit of the windowed interval on superpixel n. Otherwise, in order to further reduce the image size, enter the next stage (lower resolution): l = l + 1, and repeat the above steps (1)-(5).

[0090] Finally, the preprocessed lidar three-dimensional point cloud data has a data size of N superpixel ×T w , where T w is the total number of time intervals, and N superpixel is the number of superpixels. Compared with the original point cloud data s, the data size is reduced from N r ×N c ×T to N superpixel ×T w . The size of N superpixel is affected by the intensity threshold R: the larger R is, the smaller N superpixel is, and the more target details are lost, and vice versa. In practical applications, users need to set R according to the characteristics of the target and the size of the image.

[0091] Step 2: Using an optimization framework, based on the Poisson distribution model of the preprocessed lidar three-dimensional point cloud data, a depth estimation cost function combining the data term and the regularization term introduced by the intensity information is established:

[0092] Based on the Poisson distribution model of the point cloud data obtained by preprocessing , construct the data term L of the cost function, and the specific formula is:

[0093]

[0094] In the formula, is the photon count of the t-th time interval on superpixel n, t ∈ {T m (n),..., T m (n)+T w}, n ∈ {1,..., N superpixel}, and t n is the depth value to be solved on superpixel n.

[0095] Use the Canny operator to obtain the edge amplitude matrix M of the intensity map. M is introduced as an indicator function into the TV (Total Variation, TV) regularization model, and finally establish the depth estimation cost function C combined with the data term, and the specific formula is:

[0096]

[0097] where λ is the regularization parameter, t is the depth map to be solved, and TV(t)=||Dt|| 1 =||D r t|| 1 +||D c t|| 1 , where D r and D c represent the first-order differences in the horizontal and vertical directions respectively, M is the edge indicator function, and 0 < α ≤ 1: when the superpixel n is located in the edge part, M(n) is larger, and at this time λ / (1 + αM(n)) is smaller, thus reducing the diffusion in the edge part and retaining clear detail features. On the contrary, a smaller M(n) represents that the superpixel n is located in the smooth area, and at this time a larger λ / (1 + αM(n)) enhances the diffusion in the smooth area and effectively filters out noise.

[0098] Step 3: Use the Alternating Direction Method of Multipliers (ADMM) to estimate the depth image from the cost function. The memory requirement of the ADMM algorithm is proportional to 6 times the data size. Since the data size is reduced from N r ×N c ×T to N superpixel ×T w during data preprocessing, the corresponding memory requirement and calculation time can be significantly reduced. The specific method is as follows:

[0099] ADMM usually solves the equality-constrained optimization problem of two variables, as shown below:

[0100] minimize f(u)+g(v)

[0101] subject to Au + Bv = c

[0102] where u ∈ R m×n , v ∈ R p×q , A ∈ R d×m , B ∈ R d×p , c ∈ R d×1 , and g and f are appropriate closed convex functions. Referring to the structure of the optimization problem and introducing the Lagrange operator, the cost function C is extended to the augmented Lagrangian form and simplified as follows:

[0103]

[0104] where ρ > 0 is the penalty parameter, t = {t n : 1 ≤ n ≤ Nsuperpixel} is the depth map to be solved, N superpixel is the number of superpixels, v = t is the constructed equation constraint condition, and d is the Lagrange operator.

[0105] The augmented Lagrangian function is solved using the ADMM algorithm, and the solution process includes the following three independent sub-problems:

[0106]

[0107]

[0108] d k+1 = d k + t k+1 - v k+1

[0109] For superpixel n, the solution formula for t k+1 can be expressed as:

[0110]

[0111] In the formula, Therefore, it can be solved by traversing its range. Since time-domain windowing reduces the length of the time dimension, the computational cost of this scheme is acceptable. v k+1 The solution formula for is the weighted TV regularization model. d k+1 The solution formula for can be directly updated. Iteratively solve the three sub-problems. When max(||t k+1 - t k || 2 , ||v k+1 - v k || 2 , ||d k+1 - d k || 2 ) < tol, where tol is the tolerance value, stop the iteration. Finally, the solution t that meets the accuracy is obtained, representing the estimated depth image. As Figure 3 and 4 shown, respectively, are the target depth map and 3D point cloud map obtained by the intensity-guided long-range lidar depth estimation method of the present invention when the sampling integration time is 10 ms and 0.5 ms. Figure 7 is the target depth map and 3D point cloud map obtained by the intensity-guided long-range lidar depth estimation method of the present invention when the sampling integration time is 10 ms. Figure 8It is the target depth map and 3D point cloud map obtained by the intensity-guided long-distance lidar depth estimation method of the present invention at a sampling integration time of 0.5 ms. At integration times of 10 ms and 0.5 ms, the sizes of the point cloud data obtained by data preprocessing are 22360×12 and 24728×12 respectively, which are 0.22% and 0.25% of the size of the original data.

Claims

1. A method for depth estimation of a long-distance lidar based on intensity guidance, characterized in that, the specific steps are as follows: Step 1: Utilize the spatio-temporal correlation of intensity information and echo signals to alternately perform the establishment of multi-scale superpixels and fast time-domain windowing on the lidar three-dimensional point cloud data to complete data preprocessing. The specific method is as follows: Step 1.1: Set the current stage as l, where 1 ≤ l ≤ l max , and l max is the preset maximum number of stages; Step 1.2: For the intensity map r l-1 Use the Canny edge detection operator to obtain a binary matrix D that marks the edges, where l , D l (i,j) is an element in matrix D l ; and use matrix F l to mark the empty pixels and noise pixels in the depth map t, where l-1 , F l (i,j) is an element in matrix F l ; Step 1.3: Establish superpixels. For each pixel (i, j), the superpixel N i,j has the following specific formula: N i,j = (x, y) ∈ {i - l..., i + l} × {j - l..., j + l}: F l-1 (x,y) = 0, D l-1 (x,y) = D l-1 (i,j), wherein, is the intensity difference between pixel (x, y) and pixel (i, j), and R is the threshold value; For pixels (x, y) ∈ N i,j , being on the same edge as pixel (i, j) and neither being an empty pixel nor a noise pixel, when the intensity difference is less than the threshold R, merge them into a new superpixel n. At this time, the photon counting histogram is expressed as Regard these pixels as the new superpixel n, where T is the total number of time units; For pixels (x,y) ∈ {i-l...,i+l} × {j-l...,j+l} and F l-1 (x,y) = 0, but D l-1 (x,y) ≠ D l-1 (i,j) or Retain the superpixels at the original scale; For pixels (x, y) ∈ {i - l..., i + l} × {j - l..., j + l} and F l-1 (x, y) = 1, which are empty pixels or noise pixels, the photon counting histogram is expressed as The quotient of the number of pixels in the superpixel n; Step 1.4: Perform fast time-domain windowing on the histogram with T time units for each pixel (i, j). After windowing, the data is obtained where is the photon count of the t-th time unit on pixel (i, j) after being processed in stage l, and T m l (i, j) is the first time unit of the windowed interval, and T w is the preset number of time units; The specific method of fast time-domain windowing is as follows: For a histogram s with the number of time units being T for each pixel (i, j) i,j = {s i,j,t : 1 ≤ t ≤ T}, the following steps are performed: Step 1.4.1: Set the initial search value: Step 1.4.2: Divide the histogram into three subsets of equal length and calculate the sum of their photon counts: In the formula, Step 1.4.3: Compare Retain the subset with the largest sum of photon counts and update the search value Thus, is used to index the first time unit of the set after windowing; Step 1.4.4: If the condition is satisfied Then obtain the set with the most photon counts at pixel (i, j) T w Is the number of time units, and go to Step 1.4.5, otherwise go to Step 1.4.2; Step 1.4.5: Set the threshold K according to the noise level of the histogram. If the sum of the photon counts of the set is less than K, then ignore all the counts of this pixel; otherwise, retain the photon counts within the interval of length T . Obtain the response set on the pixel (i, j): w ​ Step 1.5: Obtain the intensity map r l and the depth map t l , and the specific formula is as follows: In the formula, is the photon count of the t-th time unit after the l-th stage of processing on the pixel (i, j), N r is the number of row pixels, N c is the number of column pixels; When l = 1, is the raw data without being processed, where s i,j,t is the photon count of the t-th time unit on the pixel (i, j), N r is the number of row pixels, N c is the number of column pixels, T is the total number of time units. At this time, the intensity map r l-1 and the depth map t l-1 are respectively: In the formula, g(t) is the pulse response function of the lidar system; Step 1.6: If the judgment condition is satisfied: If there are no empty pixels and noise pixels in the current stage l, that is or the preset maximum number of stages l is reached max , the processed lidar point cloud data is obtained where is the photon count of the t-th time unit on the superpixel n, and N superpixel is the total number of superpixels, and T m (n) is the first time unit of the windowed interval on the superpixel n. Otherwise, l = l + 1, and steps 1.1 to 1.6 are repeated; Step 2: Adopt an optimization framework to establish a depth estimation cost function that combines a data term with a regularization term introduced by intensity information based on the Poisson distribution model of the preprocessed lidar three-dimensional point cloud data; Step 3: Use the alternating direction multiplier algorithm to estimate the depth image from the cost function.

2. The method for depth estimation of a long-distance lidar based on intensity guidance according to claim 1, characterized in that, the specific method for constructing the depth estimation cost function is: Based on the point cloud data obtained by preprocessing , a Poisson distribution model is used to construct the data term L of the cost function. The specific formula is as follows: where \(g(t)\) is the pulse response function of the lidar system, is the photon count in the \(t\)-th time interval on superpixel \(n\), \(t\in\{T m (n),\ldots,T m (n)+T w}\), \(n\in\{1,\ldots,N superpixel}\), where \(T m (n)\) is the first time unit of the windowed interval on superpixel \(n\), \(T w is the number of time units, \(N superpixel is the number of superpixels, \(t n is the depth value to be solved on superpixel \(n\); Use the Canny operator to obtain the edge amplitude matrix M of the intensity map, introduce the edge amplitude matrix M as an indicator function into the total variation regularization model TV, and establish a depth estimation cost function C combined with the data term, specifically: where λ is the regularization parameter, and the total variation model TV(t) = ||Dt|| 1 = ||D r t|| 1 + ||D c t|| 1 , where D r and D c represent the first-order differences in the horizontal and vertical directions respectively, M is the edge indicator function, α is a constant greater than 0, and t is the depth map to be solved.

3. The method for depth estimation of a long-distance lidar based on intensity guidance according to claim 1, characterized in that, the specific method for using the alternating direction multiplier algorithm to estimate the depth image from the cost function is: Expand the cost function C into an augmented Lagrangian form and simplify it, specifically: where \(t = \{t n : 1\leq n\leq N superpixel \}\) is the depth map to be solved, \(N superpixel \) is the number of superpixels, \(v = t\) is the constructed equality constraint, \(d\) is the Lagrangian operator, \(\rho>0\) is the penalty parameter, \(\lambda\) is the regularization parameter, \(\alpha\) is a constant greater than 0, \(g(t)\) is the lidar system impulse response function, \(TV\) is the total variation model, \(M\) is the edge indicator function, \) is the photon count in the \(t\)-th time interval on superpixel \(n\), \(t\in\{T m (n),\ldots,T m (n)+T w \}\), \(n\in\{1,\ldots,N superpixel \}\), where \(T m (n)\) is the first time unit of the windowed interval on superpixel \(n\), \(T w \) is the number of time units, \(t n \) is the depth value to be solved on superpixel \(n\); Use the ADMM algorithm to solve the augmented Lagrangian function, and the solution process includes the following three independent sub-problems: d k+1 = d k + t k+1 - v k+1 Iteratively solve three sub-problems. When max(||t k+1 - t k || 2 , ||v k+1 - v k || 2 , ||d k+1 - d k || 2 ) < tol, where tol is the tolerance value, stop the iteration, and obtain the solution t that meets the accuracy. t is the estimated depth image.

Images