Airborne laser radar full-waveform data real-time acquisition and processing method

Abstract

The invention provides an airborne laser radar full-waveform data real-time acquisition and processing method, and relates to the technical field of laser radar data processing, and the method specifically comprises the steps: collecting original full-waveform data through an acquisition card, and changing the acquisition and storage mode of the acquisition card into the mode that 5000 lines of data are collected and stored as an original data file, setting a thread count according to a computer configuration condition, and setting a thread task; each thread executes a Gaussian decomposition algorithm to decompose each echo component, and waveform fitting of full waveform data is carried out; using an LM algorithm to optimize the waveform parameters; and carrying out data merging on the full-waveform data according to the file time information of the processed full-waveform data to enable the distance information to be complete, and finally outputting the processed data. According to the technical scheme, the problems that in the prior art, the data size of airborne laser radar full-waveform data is large, a large amount of data needs to be copied to a computer, and the airborne laser radar full-waveform data processing efficiency is low are solved.

Description

Technical Field

[0001] The present invention relates to the technical field of laser radar data processing, and in particular to a method for real-time acquisition and processing of full waveform data of an airborne laser radar. Background Art

[0002] The waveform of a full-waveform laser radar refers to the function of the energy of a laser emission pulse or a received pulse changing over time. At present, the amount of full-waveform data of airborne laser radar is very large, collecting 9.6G of data per minute. Each flight generally takes 20-30 minutes, and more than ten minutes of data are collected. Therefore, the amount of collected data is very large, and it is impossible to achieve fast and real-time processing, that is, it is impossible to achieve simultaneous collection and processing. The data processing methods in the prior art include deconvolution, waveform decomposition, LM waveform decomposition, and some improved algorithms based on them. However, these algorithms all require the data to be copied from the radar system hard disk to the computer before waveform extraction and processing, which greatly reduces work efficiency. In addition, the algorithm processing time is relatively long, and it takes ten to twenty minutes to process one minute of data.

[0003] Therefore, there is a need for a real-time acquisition and processing method for airborne laser radar full waveform data that does not require copying the data to a computer and has high processing efficiency. Summary of the invention

[0004] The main purpose of the present invention is to provide a real-time acquisition and processing method for airborne laser radar full waveform data, so as to solve the problems in the prior art that the data volume of airborne laser radar full waveform data is large, a large amount of data needs to be copied to a computer, and the processing efficiency of airborne laser radar full waveform data is low.

[0005] To achieve the above object, the present invention provides a method for real-time acquisition and processing of full waveform data of an airborne laser radar, which specifically comprises the following steps: S1, using the acquisition card to collect the original full waveform data, and changing the acquisition card acquisition storage mode to store every 5000 lines of data collected as a raw data file, that is, storing one raw data file per second in chronological order.

[0006] S2, set the number of threads and thread tasks according to the computer configuration and CUDA version.

[0007] S3, each thread executes the Gaussian decomposition algorithm to decompose each echo component and perform waveform fitting of the full waveform data.

[0008] S4, use the LM algorithm to optimize the waveform parameters until the residual between the original waveform and the fitted waveform is less than the standard deviation of the noise Or the waveform decomposition parameter reaches the upper limit , stop entering initial waveform parameters.

[0009] S5, merging the full waveform data according to the file time information of the processed full waveform data to make the distance information complete, and finally outputting the processed data.

[0010] Furthermore, in step S2, the number of threads is set according to the computer configuration, specifically: the number of threads cannot exceed the computer configuration and meets the processing speed.

[0011] Further, step S3 includes the following steps: S3.1, the beginning and end waveforms are used as the basis for evaluating the background noise. If the noise mean value at the front end of the original data waveform is <Noise standard deviation of the front end of the raw data waveform , then the background noise ,otherwise = , background noise standard deviation = , is the noise mean at the tail of the original data waveform, the noise threshold Calculation formula: +3 .

[0012] S3.2, use Gaussian smoothing filter to filter out background noise.

[0013] S3.3, using peak detection and inflection point detection, estimate the initial values ​​of the waveform parameters in the ideal Gaussian mixture model. The waveform parameters include: the amplitude, center position and half-height width of the echo component. The ideal Gaussian mixture model is: ; , , Respectively represent The amplitude, center position and half-height width of each echo component; is the background noise of the original full waveform signal.

[0014] S3.4, sort the echo components by time or intensity using a quick sorting method.

[0015] Furthermore, step S3.3 specifically includes the following steps: S3.3.1, calculate the first-order derivative of each sampling point, and use the convolution kernel [-1, 0, 1] as a differential operator to approximate the first-order derivative of the signal. Slide the convolution kernel [-1, 0, 1] along the sampling points of the signal, and perform a weighted summation operation at each position. The zero-crossing value point of the first-order derivative from positive to negative is regarded as the peak position.

[0016] S3.3.2, calculate the second-order derivative of each sampling point, and use the convolution kernel [-1, -2, 1] as the differential operator to approximate the second-order derivative of the signal. The zero-crossing value of the second-order derivative is regarded as the inflection point, and the value less than the noise threshold is eliminated. The peak point and inflection point.

[0017] S3.3.3, for each peak point, if there is no inflection point in the monotonic interval around the peak point, then the peak point does not belong to any echo component, otherwise calculate and , is the mean of the inflection point positions in the monotonic interval to the left of the peak point, is the mean of the inflection point positions in the right monotonic interval of the peak point, and then calculate and , is the distance between the inflection point and the peak point in the monotonic interval to the left of the peak point, is the distance between the inflection point and the peak point in the right monotonic interval of the peak point. and Both are smaller than the half-width of the transmitted wave When the peak point is half of the value of the echo component, the peak point does not belong to any echo component; otherwise, the peak point belongs to a certain echo component.

[0018] Furthermore, step S3.3 also includes the following steps: S3.3.4, estimate the initial value of the echo component where the peak point is located, and use the amplitude of the peak point as the initial value of the pulse amplitude The peak point position is taken as the center position, and the initial value of the center position is , initial value of half-height width The calculation formula is as follows: = ; In the formula, the half-width of the emission wave is ,when , Not less than When, take , The smaller value is twice as The initial value of , One of which is not less than , the other is smaller than When, take , The larger value is twice as The initial value of .

[0019] S3.3.5, for the echo component in the initial judgment, detect the inflection points in the left and right intervals of the echo component, if there is a point that does not belong to the interval [ - ] continuous inflection points, where is the horizontal coordinate position of the current echo component, and the maximum distance of the continuous inflection points is greater than , it means that there may be other echo components nearby, and the larger amplitude value among the continuous inflection points is used as the initial value of the pulse amplitude of the new echo component. , the position of the inflection point with a larger amplitude is taken as the initial value of the center position of the new echo component , take twice the maximum distance between the continuous inflection points as the initial value of the half-height width of the new echo component .

[0020] Further, step S4 is specifically as follows: input the initial waveform parameters into the ideal Gaussian mixture model in the order of the echo components obtained in step S3.4, and remove the waveform parameters whose amplitude is less than the noise threshold or whose half-height width is less than and merge the echo components with a relative distance less than Until the residual between the original waveform and the fitted waveform is less than the standard deviation of the noise , or the waveform parameters of the ideal Gaussian mixture model fit reach the upper limit , stop inputting the initial waveform parameters and output the fitted waveform parameters.

[0021] Furthermore, in step S5, the full waveform data with high waveform fitting accuracy is converted into distance information as follows: ; in, is the distance information, is the round trip time of the laser, is the speed of light; ; in, is the mean value of the echo component.

[0022] The present invention has the following beneficial effects: The present invention changes the acquisition and storage method, and changes the full waveform data of the original .out file with a large amount of data stored every minute to a .out file with a small amount of data stored every second, and stores the file in the form of time naming, so that the acquired data is not easy to overflow, and provides effective support for the merging of data after subsequent Gaussian decomposition processing.

[0023] It adopts multi-threaded parallel processing mode. Since the tasks have been decomposed during the collection, each thread performs the same task, which increases the processing efficiency and improves the real-time performance.

[0024] Finally, based on the processed files, the time-marked file names are merged front and back into the overall data information. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] In order to more clearly illustrate the specific implementation of the present invention or the technical solutions in the prior art, the following will briefly introduce the drawings required for the specific implementation or the prior art description. Obviously, the drawings described below are some implementations of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work. In the drawings: Figure 1 A flow chart of a method for real-time acquisition and processing of full waveform data of an airborne laser radar according to the present invention is shown. DETAILED DESCRIPTION

[0026] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0027] like Figure 1 The method for real-time acquisition and processing of full waveform data of an airborne laser radar shown in the figure specifically comprises the following steps: S1, using the acquisition card to collect the original full waveform data, and changing the acquisition and storage mode of the acquisition card to store every 5000 lines of data collected as a raw data file, that is, storing a raw data file in chronological order every second.

[0028] Use ADQ7DC acquisition card to collect full waveform data. The sampling frequency of the acquisition card is 10GHz. The sampling rate of the digital-to-analog converter ADC is set to 16834hz in the system. Set the acquisition and saving mode of the ADQ7DC acquisition card to 1 storage file per second. In this way, the .out file stores one data file per second, with a data volume of approximately 0.16G. Then, through multi-threading, 9.6G of data can be split into 60 .out files with a data volume of 0.16G. Each .out file stores 5000 rows and 16384 columns of data.

[0029] S2, according to the computer configuration and CUDA version, set the number of threads and thread tasks. According to the order of the acquisition card to transfer and save files, in C++ <thread>Or use Python's threading standard library module to set thread tasks. These tasks are independent and do not depend on the order of processing. Using more than 16 threads to process simultaneously can achieve fast real-time processing effects.

[0030] S3, each thread executes the Gaussian decomposition algorithm to decompose each echo component and perform waveform fitting of the full waveform data.

[0031] Each thread executes the Gaussian decomposition algorithm to decompose each echo component, calculates the distance information corresponding to each component to achieve the purpose of Gaussian decomposition, and then merges the decomposed data according to the file naming time (file information named by hours, minutes and seconds) to obtain continuous distance information. Since the .out file is 0.16G in size, the speed of processing a file is increased by 60 times. A .out full waveform file can be processed in 5 seconds. There are 16 threads processing at the same time, which is faster. The amount of data when the data is finally merged is the processed data, so the efficiency is improved.

[0032] S4, use the LM algorithm to optimize the waveform parameters until the residual between the original waveform and the fitted waveform is less than the standard deviation of the noise Or the waveform decomposition parameter reaches the upper limit , stop entering initial waveform parameters.

[0033] S5, merging the full waveform data according to the file time information of the processed full waveform data to make the distance information complete, and finally outputting the processed data.

[0034] Specifically, in step S2, the number of threads is set according to the computer configuration, specifically: the number of threads cannot exceed the computer configuration and meets the processing speed.

[0035] Specifically, step S3 includes the following steps: S3.1, the range of the detected target is uncertain, which may cause the presence of echo components at the beginning or end, so the noise level may be overestimated. However, the possibility of echo components at the beginning and end is very low. Therefore, the waveform at the beginning and end is used as the basis for evaluating the background noise. If the noise mean value at the front end of the original data waveform is <Noise standard deviation of the front end of the raw data waveform , then the background noise ,otherwise = , background noise standard deviation = , is the noise mean at the tail of the original data waveform, the noise threshold Calculation formula: +3 .

[0036] S3.2, use Gaussian smoothing filter to filter out background noise.

[0037] S3.3, using peak detection and inflection point detection, estimate the initial values ​​of the waveform parameters in the ideal Gaussian mixture model. The waveform parameters include: the amplitude, center position and half-height width of the echo component. The ideal Gaussian mixture model is: ; , , Respectively represent The amplitude, center position and half-height width of each echo component; is the background noise of the original full waveform signal. The optimal solution of 3*k echo components and one background noise is calculated by Gaussian mixture model. The peak detection and inflection point detection can be used to estimate the initial value of the full waveform more accurately, which is conducive to the fitting waveform being closer to the real waveform.

[0038] S3.4, sort the echo components by time or intensity using the quick sort method. The echo components can be sorted by time or intensity using the quick sort method, which is convenient for understanding the echo characteristics and subsequent processing (such as target identification) is more efficient, and some outliers can be easily eliminated during sorting.

[0039] Specifically, step S3.3 includes the following steps: S3.3.1, calculate the first-order derivative of each sampling point, and use the convolution kernel [-1, 0, 1] as a differential operator to approximate the first-order derivative of the signal. Slide the convolution kernel [-1, 0, 1] along the sampling points of the signal, and perform a weighted summation operation at each position. The zero crossing value (from positive to negative) of the first-order derivative is regarded as the peak position. Calculate the first-order derivative of each sampling point, and use the convolution kernel [-1, 0, 1] as a differential operator to approximate the first-order derivative of the signal. For each sampling point in the signal , where i is the index of the sampling point. If the sampling point If there are two adjacent sampling points, calculate the sampling points The weighted sum of two adjacent sampling points, where the weight is given by the middle value of the convolution kernel.

[0040] S3.3.2, calculate the second-order derivative of each sampling point, and use the convolution kernel [-1, -2, 1] as the differential operator to approximate the second-order derivative of the signal. The second-order derivative zero-crossing value (from positive to negative and from negative to positive) is regarded as the inflection point, and the value less than the noise threshold is eliminated. The peak point and inflection point.

[0041] S3.3.3, for each peak point, if there is no inflection point in the monotonic interval around the peak point, then the peak point does not belong to any echo component, otherwise calculate and , is the mean of the inflection point positions in the monotonic interval to the left of the peak point, is the mean of the inflection point positions in the right monotonic interval of the peak point, and then calculate and , is the distance between the inflection point and the peak point in the monotonic interval to the left of the peak point, is the distance between the inflection point and the peak point in the right monotonic interval of the peak point. and Both are smaller than the half-width of the transmitted wave When the peak point is half of the value of the echo component, the peak point does not belong to any echo component; otherwise, the peak point belongs to a certain echo component.

[0042] Specifically, step S3.3 also includes the following steps: S3.3.4, estimate the initial value of the echo component where the peak point is located, and use the amplitude of the peak point as the initial value of the pulse amplitude The peak point position is taken as the center position, and the initial value of the center position is , initial value of half-height width The calculation formula is as follows: = ; In the formula, the half-width of the emission wave is ,when , Not less than When, take , The smaller value is twice as The initial value of , One of which is not less than , the other is smaller than When, take , The larger value is twice as The initial value of .

[0043] S3.3.5, for the echo of the overlapping part that may be missed in the initial judgment, make a second judgment to ensure the accuracy of the initial value estimation. For the echo component in the initial judgment, detect the inflection point in the left and right intervals of the echo component. If there is an inflection point that does not belong to the interval [ - ] continuous inflection points, where is the horizontal coordinate position of the current echo component, and the maximum distance of the continuous inflection points is greater than , it means that there may be other echo components nearby, and the larger amplitude value among the continuous inflection points is used as the initial value of the pulse amplitude of the new echo component. , the position of the inflection point with a larger amplitude is taken as the initial value of the center position of the new echo component , take twice the maximum distance between the continuous inflection points as the initial value of the half-height width of the new echo component .

[0044] Specifically, step S4 is as follows: use the LM algorithm to optimize the estimation of the waveform parameters, input the initial waveform parameters into the ideal Gaussian mixture model in the order of the echo components obtained in step S3.4, iteratively calculate the optimal solution, and update the waveform parameters fitted by the ideal Gaussian mixture model every time a set of waveform parameters is added, and eliminate the waveform parameters whose amplitude is less than the noise threshold or whose half-height width is less than and merge the echo components with a relative distance less than Until the residual between the original waveform and the fitted waveform is less than the standard deviation of the noise , or the waveform parameters of the ideal Gaussian mixture model fit reach the upper limit , stop inputting the initial waveform parameters and output the fitted waveform parameters.

[0045] Specifically, in step S5, the full waveform data with high waveform fitting accuracy is converted into distance information as follows: ; in, is the distance information, is the round trip time of the laser, is the speed of light; ; in, is the mean value of the echo component.

[0046] Of course, the above description is not a limitation of the present invention, and the present invention is not limited to the above examples. Changes, modifications, additions or substitutions made by technicians in this technical field within the essential scope of the present invention should also fall within the protection scope of the present invention.< / thread>

Claims

1. A method for real-time acquisition and processing of full waveform data of an airborne laser radar, characterized in that: The specific steps include: S1, using an acquisition card to collect the original full waveform data, and changing the acquisition card acquisition storage mode to store every 5000 lines of data collected as a raw data file, that is, storing one raw data file per second in chronological order; S2, set the number of threads according to the computer configuration and set the thread tasks; S3, each thread executes the Gaussian decomposition algorithm to decompose each echo component and perform waveform fitting of the full waveform data; S4, use the LM algorithm to optimize the waveform parameters until the residual between the original waveform and the fitted waveform is less than the standard deviation of the noise Or the waveform decomposition parameter reaches the upper limit , stop inputting initial waveform parameters; S5, merging the full waveform data according to the file time information of the processed full waveform data to make the distance information complete, and finally outputting the processed data.

2. The method for real-time acquisition and processing of full waveform data of an airborne laser radar according to claim 1, characterized in that: In step S2, the number of threads is set according to the computer configuration, specifically: the number of threads cannot exceed the computer configuration and meets the processing speed.

3. The method for real-time acquisition and processing of full waveform data of an airborne laser radar according to claim 1, characterized in that: Step S3 includes the following steps: S3.1, the beginning and end waveforms are used as the basis for evaluating the background noise. If the noise mean value at the front end of the original data waveform is <Noise standard deviation of the front end of the raw data waveform , then the background noise ,otherwise = , background noise standard deviation = , is the noise mean at the tail of the original data waveform, the noise threshold Calculation formula: +3 ； S3.2, use Gaussian smoothing filter to filter out background noise; S3.3, using peak detection and inflection point detection, estimate the initial values ​​of the waveform parameters in the ideal Gaussian mixture model. The waveform parameters include: the amplitude, center position and half-height width of the echo component. The ideal Gaussian mixture model is: ； , , Respectively represent The amplitude, center position and half-height width of each echo component; is the background noise of the original full waveform signal; S3.4, sort the echo components by time or intensity using a quick sorting method.

4. The method for real-time acquisition and processing of full waveform data of an airborne laser radar according to claim 3 is characterized in that: Step S3.3 specifically includes the following steps: S3.3.1, calculate the first-order derivative of each sampling point, use the convolution kernel [-1, 0, 1] as the difference operator to approximate the first-order derivative of the signal, slide the convolution kernel [-1, 0, 1] along the sampling points of the signal, and perform a weighted summation operation at each position, and regard the zero-crossing value point of the first-order derivative from positive to negative as the peak position; S3.3.2, calculate the second-order derivative of each sampling point, and use the convolution kernel [-1, -2, 1] as the differential operator to approximate the second-order derivative of the signal. The zero-crossing value of the second-order derivative is regarded as the inflection point, and the value less than the noise threshold is eliminated. The peak point and inflection point; S3.3.3, for each peak point, if there is no inflection point in the monotonic interval around the peak point, then the peak point does not belong to any echo component, otherwise calculate and , is the mean of the inflection point positions in the monotonic interval to the left of the peak point, is the mean of the inflection point positions in the right monotonic interval of the peak point, and then calculate and , is the distance between the inflection point and the peak point in the monotonic interval to the left of the peak point, is the distance between the inflection point and the peak point in the right monotonic interval of the peak point. and Both are smaller than the half-width of the transmitted wave When the peak point is half of the value of the echo component, the peak point does not belong to any echo component; otherwise, the peak point belongs to a certain echo component.

5. The method for real-time acquisition and processing of full waveform data of an airborne laser radar according to claim 4, characterized in that: Step S3.3 also includes the following steps: S3.3.4, estimate the initial value of the echo component where the peak point is located, and use the amplitude of the peak point as the initial value of the pulse amplitude , the peak point position is taken as the center position, and the initial value of the center position is , initial value of half-height width The calculation formula is as follows: = ； In the formula, the half-width of the emission wave is ,when , Not less than When, take , The smaller value is twice as The initial value of , One of which is not less than , the other is smaller than When , The larger value is twice as The initial value of S3.3.5, for the echo component in the initial judgment, detect the inflection points in the left and right intervals of the echo component, if there is a point that does not belong to the interval [ - ] continuous inflection points, where is the horizontal coordinate position of the current echo component, and the maximum distance of the continuous inflection points is greater than , it means that there may be other echo components nearby, and the larger amplitude value among the continuous inflection points is used as the initial value of the pulse amplitude of the new echo component. , the position of the inflection point with a larger amplitude is taken as the initial value of the center position of the new echo component , take twice the maximum distance between the continuous inflection points as the initial value of the half-height width of the new echo component .

6. The method for real-time acquisition and processing of full waveform data of an airborne laser radar according to claim 1, characterized in that: Step S4 is specifically as follows: input the initial waveform parameters into the ideal Gaussian mixture model in the order of the echo components obtained in step S3.4, and remove the waveform parameters whose amplitude is less than the noise threshold or whose half-height width is less than and merge the echo components with a relative distance less than Until the residual between the original waveform and the fitted waveform is less than the standard deviation of the noise , or the waveform parameters of the ideal Gaussian mixture model fit reach the upper limit , stop inputting the initial waveform parameters and output the fitted waveform parameters.

7. The method for real-time acquisition and processing of full waveform data of an airborne laser radar according to claim 1, characterized in that: In step S5, the full waveform data with high waveform fitting accuracy is converted into distance information as follows: ； in, is the distance information, is the round trip time of the laser, is the speed of light; ； in, is the mean value of the echo component.

Images