A full waveform decomposition method with oscillatory signals in an airborne lidar

Abstract

The application discloses a full waveform decomposition method with an oscillation signal in an airborne laser radar, belongs to the technical field of laser radar data processing, and is used for waveform decomposition and comprises the following steps: reading and converting original data, performing data blocking, selecting a wavelet threshold for denoising of each block of data, constructing a mixed model, selecting initial parameters of the mixed model, performing particle swarm optimization, obtaining optimized parameters, constructing output data of each block of data, splicing all the output data, and completing full waveform decomposition. The application removes almost all high-frequency noises, has a good retention degree for signals, directly performs full waveform decomposition, fits an oscillation signal, reduces the influence caused by the oscillation signal, and extracts distance data of the airborne laser radar. The multi-point search and randomness of the particle swarm make it more possible to find a global optimal solution, the particle updating is naturally parallelized, is suitable for distributed computing, and the calculation speed is accelerated.

Description

TECHNICAL FIELD

[0001] The application discloses a full waveform decomposition method with oscillation signals in an airborne laser radar, and belongs to the technical field of laser radar data processing. BACKGROUND

[0002] In an airborne laser radar system, laser is first generated by a laser, and after unnecessary light is filtered through a telescope and a lens, the received laser signal is converted into an electric signal through a photomultiplier tube. When the laser reaches a target object, a backscattering signal is generated and is received by the airborne laser radar receiving system, and is represented as a wave peak in full waveform data. Due to impedance mismatching, the entire circuit will generate low-frequency oscillation due to the change of the input of the echo signal. When the echo electric signal is collected by the collection system, the oscillation signal will also be collected by the collection system, and the oscillation signal will be superimposed on the echo signal and form a longer tail. In the decomposition of the full waveform, an approximate function is mainly used to fit the echo signal hitting the target object, and the extreme point is extracted as the measured distance. The oscillation signal not only affects the shape of the entire waveform, but also causes deviation when the function optimizes the fitting of the echo signal, and the oscillation signal also causes the extreme point of the echo signal to deviate, thereby changing the distance and affecting the final data accuracy. SUMMARY

[0003] The application aims to provide a full waveform decomposition method with oscillation signals in an airborne laser radar, so as to solve the problem that the oscillation signal affects the waveform decomposition of the laser radar data in the prior art.

[0004] A full waveform decomposition method with oscillation signals in an airborne laser radar comprises the following steps: reading and converting original data, performing data blocking, selecting a wavelet threshold for denoising each block of data, constructing a hybrid model, selecting initial parameters of the hybrid model, performing particle swarm optimization, obtaining optimized parameters, constructing output data of each block of data, splicing all the output data, and completing full waveform decomposition.

[0005] The airborne laser radar comprises an optical system and a collection system, and the optical system and the collection system share the same square wave with a frequency of 5 kHz. The square wave controls the laser emission of the optical system, and the square wave controls the frequency collection of the collection system. The optical system and the collection system are triggered at the edges of the square wave, so that the echo signal corresponds to the collected data one by one. The optical system receives the echo signal, and the collected data is a quantized voltage value, which is a negative number.

[0006] Reading and converting the original data comprises reading binary numbers in groups of 16 bits into decimal numbers, and after reading is completed, performing bit inversion on the decimal numbers, converting the negative numbers into positive numbers for processing, and storing all the data in the form of tab character intervals.

[0007] All data blocks are divided, different blocks are numbered, and are distributed to the processes of the central processor for processing, after each block of data is processed, the processed data is sorted according to the number, and finally the data is spliced and stored, the distance data is extracted from the original data, and the distance data is used as the full waveform data of the airborne laser radar.

[0008] Wavelet threshold denoising includes selecting adaptive decomposition layers, discrete wavelet transform, determining wavelet threshold conforming to the system and wavelet inverse transform;

[0009] Determining the decomposition layer includes cutting and extracting the front data of the distance data, and determining the length of the collected data , the decomposition layer is:

[0010] ;

[0011] In the formula, is the laser pulse width, is the frequency of the airborne laser radar system;

[0012] Discrete wavelet transform includes performing discrete wavelet transform on the waveform data as a given signal:

[0013] ;

[0014] In the formula, is the detail coefficient obtained by discrete wavelet transform, is the scale of the stretching and shrinking of the wavelet, is the wavelet base function, is the translation amount for controlling the position of the wavelet base, is the given signal, is the time.

[0015] Determining the wavelet threshold includes calculating the energy density of the detail coefficients in each layer of the decomposition layer :

[0016] ;

[0017] ;

[0018] In the formula, is the detail coefficient of the first layer, represents the first detail coefficient in a layer, is the length of the first layer;

[0019] Let , find the density difference that meets the following judgment condition :

[0020] ;

[0021] In the formula, is an empirical value, represents the median, is the energy density of the first layer;

[0022] For the signal satisfying the judgment condition , the noise soft threshold is estimated from the sub-band of the corresponding first layer discrete wavelet transform :

[0023] ;

[0024] The wavelet threshold is:

[0025] ;

[0026] The wavelet inverse transform is performed on the detail coefficients higher than the wavelet threshold, and the given signal after filtering out the noise is reconstructed.

[0027] Constructing a hybrid coupling model conforming to an airborne laser radar system , comprising:

[0028] ;

[0029] In the formula, is a natural constant, is the amplitude of the Gaussian function, is the angular frequency, is the amplitude of the damped oscillation function, is the attenuation factor, is the Gaussian mean, is the signal standard deviation, is the time scale;

[0030] Selecting the initial parameters of the hybrid model includes selecting the estimated Gaussian function initial parameters and the damped oscillation function initial parameters, performing multi-level wavelet decomposition on the given signal after filtering out the noise, and extracting sub-signals of different frequency bands, and the multi-level wavelet decomposition is:

[0031] ;

[0032] In the formula, is a low-pass filter in the given signal after filtering out the noise, is a high-pass filter in the given signal after filtering out the noise, represents a scale coefficient, is the filtered approximation coefficient of the low-pass filter, for scale approximated signal, filtered detail coefficients of a low-pass filter; discrete time index, summation index, integer set;

[0033] low frequency convolution is performed on to obtain low frequency information, and high frequency convolution is performed on to obtain high frequency information, and the low frequency information is continuously subjected to low frequency convolution, and the high frequency information is fixed, and after times of decomposition, each sub-band detail coefficient from high frequency to low frequency is obtained, and the low frequency component and the intermediate frequency component are selected, and after inverse wavelet transform reconstruction, the component where the Gaussian peak is located and the component where the oscillation signal is located are obtained.

[0034] The region around the end point of the Gaussian peak is selected as the starting point of the oscillation signal, and the initial value of the Gaussian term is used as the parameter to extract the phase at the starting time of the oscillation;

[0035] Initial parameter correction is performed, and the corrected initial parameter is used as the initial particle of the particle swarm optimization algorithm.

[0036] The initial parameters of the Gaussian function include the Gaussian function :

[0037] ;

[0038] In the formula, is the independent variable of the Gaussian function, which is equivalent to the component where the Gaussian peak is located;

[0039] The intensity of the low frequency signal is used as , the center point of the low frequency signal is used as , and the half-height width of the low frequency signal is used as :

[0040] .

[0041] The initial value of the damped oscillation function includes , , and the initial phase , and the damped oscillation function is:

[0042] ;

[0043] The analytical signal is obtained:

[0044] ;

[0045] wherein, is the original signal, equivalent to the component where the oscillation signal is located;

[0046] the frequency domain after Euler transformation is:

[0047] ;

[0048] wherein, is a symbol function, is the original signal in the frequency domain;

[0049] Let , in the time domain:

[0050] ;

[0051] wherein, is an intermediate variable in the frequency domain, is an intermediate variable in the time domain, represents a convolution operation, is the original signal in the time domain;

[0052] find :

[0053] ;

[0054] wherein, is a limit symbol, represents Cauchy principal value processing of integration.

[0055] find , and then substitute it into the formula of the analytical signal to find ;

[0056] is:

[0057] ;

[0058] wherein, is a complex argument function;

[0059] Take the logarithmic transformation of the envelope signal of and perform linear fitting:

[0060] ;

[0061] obtain the slope and the intercept of the linear fitting result, , is:

[0062] ;

[0063] ;

[0064] for:

[0065] .

[0066] Particle swarm optimization involves randomly generating a swarm of particles, initializing particle positions and velocities, and calculating the mean squared error. The fitting effect, when When the mean square error of the fitting result is minimized, the current solution of the particle swarm optimization is the optimal solution. Record the individual optimal solution and the global optimal solution at this time.

[0067] calculate When fitting the results, the particles are updated as follows:

[0068] ;

[0069] ;

[0070] In the formula, For the first Subparticle update particle velocity Inertia factor and Let be the acceleration constant. Representing an interval Random numbers in the data, For the first The first particle The individual optimal solution updated by the secondary particle. Indicates the first The global optimal solution for sub-particle updates. For the first The particle position is updated by the secondary particle;

[0071] After particle swarm optimization is completed, the Gaussian function parameters and damped oscillation function parameters are converged through multiple iterations of learning, thus completing the full waveform decomposition.

[0072] Compared with existing technologies, the present invention has the following advantages: it removes a large amount of high-frequency noise and preserves the signal well; it directly performs full waveform decomposition, fits the oscillating signal, reduces its impact, and extracts the distance data of the airborne lidar; the multi-point search and randomness of the particle swarm make it more likely to find the global optimal solution, and the particle update is naturally parallelizable, suitable for distributed computing, and speeds up the computing speed. Attached Figure Description

[0073] Figure 1 for the flowchart of the present application;

[0074] Figure 2 for the fitted optimized waveform;

[0075] Figure 3 for the original waveform;

[0076] Figure 4 for the first waveform component after wavelet transform;

[0077] Figure 5 for the second waveform component after wavelet transform;

[0078] Figure 6 for the third waveform component after wavelet transform;

[0079] Figure 7 for the fourth waveform component after wavelet transform;

[0080] Figure 8 for the fifth waveform component after wavelet transform;

[0081] Figure 9 for the sixth waveform component after wavelet transform;

[0082] Figure 10 for the seventh waveform component after wavelet transform;

[0083] Figure 11 for the eighth waveform component after wavelet transform;

[0084] Figure 12 for the wavelet threshold denoising waveform. DETAILED DESCRIPTION

[0085] In order to make the objectives, technical solutions and advantages of the present application clearer, the technical solutions in the present application are described clearly and completely below. Obviously, the described embodiments are part of, but not all of the embodiments of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work fall within the scope of the present application.

[0086] A full waveform decomposition method with oscillation signals in an airborne laser radar, comprising reading and converting original data, data blocking, selecting wavelet threshold denoising for each block of data, constructing a mixed model, selecting initial parameters of the mixed model, particle swarm optimization, obtaining optimized parameters, constructing output data of each block of data, splicing all output data, and completing full waveform decomposition.

[0087] The airborne laser radar comprises an optical path system and a collection system, the optical path system and the collection system share a same square wave with a frequency of 5 kHz, the square wave controls laser emission of the optical path system and controls frequency collection of the collection system, the optical path system and the collection system are triggered at edges of the square wave, so that the echo signal corresponds to the collected data one by one, the optical path system receives the echo signal, and the collected data is quantized voltage value, which is negative.

[0088] Reading and converting the original data comprises reading and converting binary numbers in groups of 16 bits into decimal numbers, after reading is completed, inverting the decimal numbers by bit, converting negative numbers into positive numbers for processing, and storing all data in the form of tab character interval;

[0089] All data are divided into blocks, different blocks are numbered and distributed to processes of the central processor for processing, after processing of each block of data is completed, the processed data are sorted according to the number, and finally the data are spliced and stored, distance data are extracted from the original data, and the distance data are taken as full waveform data of the airborne laser radar.

[0090] Wavelet threshold denoising comprises selecting an adaptive decomposition layer number, discrete wavelet transform, determining a wavelet threshold conforming to the system and wavelet inverse transform;

[0091] Determining the decomposition layer number comprises cutting and extracting a part of the distance data in front of the data, and determining the length of the collected data , and the decomposition layer number is:

[0092] ;

[0093] In the formula, is a laser pulse width, is a collection frequency of the airborne laser radar system;

[0094] The discrete wavelet transform comprises performing the discrete wavelet transform on the waveform data as a given signal:

[0095] ;

[0096] In the formula, is a detail coefficient obtained by the discrete wavelet transform, is a scale for controlling stretching of the wavelet, is a wavelet base function, is a translation amount for controlling a position of the wavelet base, is the given signal, is time.

[0097] Determining the wavelet threshold comprises calculating energy density of the detail coefficient of each layer in the decomposition layer number :

[0098] ;

[0099] ;

[0100] wherein, is the detail coefficient of the layer, denotes the detail coefficient in the layer, is the length of the layer;

[0101] Let , find the density difference that meets the following judgment conditions:

[0102] ;

[0103] wherein, is an empirical value, denotes the median, is the energy density of the layer;

[0104] For the that meets the judgment conditions, estimate the noise soft threshold from the subband of the corresponding layer discrete wavelet transform:

[0105] ;

[0106] The wavelet threshold is:

[0107] ;

[0108] Perform wavelet inverse transform on the detail coefficients higher than the wavelet threshold, and reconstruct into the given signal filtered out of noise.

[0109] Construct a hybrid coupling model consistent with the airborne laser radar system , including:

[0110] ;

[0111] wherein, is a natural constant, is the amplitude of the Gaussian function, is the angular frequency, is the amplitude of the damped oscillation function, is the attenuation factor, is the Gaussian mean, is the signal standard deviation, is the time scale;

[0112] The initial parameters of the mixed model are selected, including selecting the initial parameters of the estimated Gaussian function and the initial parameters of the damped oscillation function, performing multi-level wavelet decomposition on the given signal after filtering noise, and extracting sub-signals in different frequency bands. The multi-level wavelet decomposition is as follows:

[0113] ;

[0114] In the formula, is a low-pass filter in the given signal after filtering noise, is a high-pass filter in the given signal after filtering noise, represents a scale coefficient, is an approximate coefficient filtered by a low-pass filter, is an approximate signal under a scale ; is a detail coefficient filtered by a low-pass filter; is a discrete time index, is a summation index, is an integer set;

[0115] Low-frequency convolution is performed on to obtain low-frequency information, and high-frequency convolution is performed on to obtain high-frequency information. The low-frequency information is continuously convolved, and the high-frequency information is fixed. After times of decomposition, the detail coefficients of each sub-band from high frequency to low frequency are obtained. The low-frequency component and the intermediate-frequency component are selected, and after inverse wavelet transform reconstruction, the component where the Gaussian peak is located and the component where the oscillation signal is located are obtained.

[0116] The region around the end point of the Gaussian peak is selected as the starting point of the oscillation signal, and the initial value of the Gaussian term is used as the parameter to extract the phase at the starting time of the oscillation.

[0117] The initial parameters are corrected, and the corrected initial parameters are used as the initial particles of the particle swarm optimization algorithm.

[0118] The initial parameters of the Gaussian function include the Gaussian function :

[0119] ;

[0120] In the formula, is the independent variable of the Gaussian function, which is equivalent to the component where the Gaussian peak is located;

[0121] The intensity of the low-frequency signal is used as , the center point of the low-frequency signal is used as , and the half-height width of the low-frequency signal is used as :

[0122] .

[0123] The initial value of the damped oscillation function includes , , and the initial phase , the damped oscillation function is:

[0124] ;

[0125] The analytical signal is obtained:

[0126] ;

[0127] In the formula, is the original signal, which is equivalent to the component in which the oscillation signal is located;

[0128] The frequency domain after Euler transformation is:

[0129] ;

[0130] In the formula, is a symbol function, is the original signal in the frequency domain;

[0131] Let , in the time domain:

[0132] ;

[0133] In the formula, is an intermediate variable in the frequency domain, is an intermediate variable in the time domain, represents the convolution operation, is the original signal in the time domain;

[0134] Find :

[0135] ;

[0136] In the formula, is a limit symbol, represents the Cauchy principal value processing of integration.

[0137] After finding , it is brought into the analytical signal formula to obtain ;

[0138] is:

[0139] ;

[0140] wherein, is a complex argument function;

[0141] The envelope signal of is logarithmically transformed and linearly fitted:

[0142] ;

[0143] The slope and the intercept of the linear fitting result are: , :

[0144] ;

[0145] ;

[0146] :

[0147] .

[0148] Particle swarm optimization includes randomly generating a group of particles, initializing particle positions and particle velocities, calculating the fitting effect of using mean square error, when the mean square error of the fitting result of is the smallest, the current particle swarm optimization solution is the optimal solution, and the particle individual optimal solution and global optimal solution at this time are recorded;

[0149] When calculating the fitting result of , the particle is updated as:

[0150] ;

[0151] ;

[0152] wherein, is the particle velocity of the th particle update, is an inertia factor, and are acceleration constants, represents a random number in the interval , is the individual optimal solution of the th particle in the th particle update, represents the global optimal solution of the th particle update, is the particle position of the th particle in the th particle update;

[0153] After the particle swarm optimization is completed, the Gaussian function parameters and the damped oscillation function parameters are converged through multiple iterations of learning to complete the full waveform decomposition.

[0154] In the power system, it is often affected by low-frequency oscillation signals, and in the airborne laser radar system used by the application, the photomultiplier tube (PMT) is also affected when converting optical-electrical signals, resulting in a decrease in data accuracy. The commonly used algorithms for processing low-frequency oscillation signals include prony algorithm, denoising Gaussian filter, mean value smoothing filter, and traditional Gaussian decomposition in the full waveform decomposition method, as well as the least square method and LM algorithm in the optimization method.

[0155] Gaussian filter is a linear smoothing filter based on Gaussian function, widely used in two-dimensional image processing field, when processing one-dimensional time series, one-dimensional Gaussian filter is used to smooth the waveform, taking the continuous odd points of one-dimensional time series data, weighting these points to get the center value, then sliding on the data according to this method, the core idea is to suppress noise through weighted average, the weight is determined by Gaussian function, the closer to the center pixel, the greater the weight, which can be used to filter out high-frequency noise. In the full waveform decomposition step, Gaussian decomposition is a process of fitting the original echo signal as a superposition of multiple Gaussian function curves. Laser height measurement system usually transmits Gaussian distributed pulse signal, and the full waveform echo signal reflected by complex terrain is a superposition of multiple echo signals, which can be considered as a superposition of multiple Gaussian signals. Full waveform decomposition is simply summarized as using Gaussian function model to fit the data, and finally obtaining the real data points.

[0156] LM algorithm is an optimization algorithm for solving nonlinear least square problems, especially suitable for curve fitting and parameter estimation. It combines the advantages of gradient descent method and Gauss-Newton method, and balances the characteristics of the two methods by dynamically adjusting the damping parameter, LM algorithm adjusts the parameter vector through iteration to minimize the objective function to update the parameters.

[0157] Fourier transform is a classical tool for signal analysis, which converts signal from time domain to frequency domain by linear combination of sinusoidal basis functions to reveal the global frequency components of the signal. However, Fourier transform has significant limitations, and wavelet transform is proposed to solve this problem. By introducing wavelet basis functions with finite support, wavelet transform can perform multi-resolution analysis on signals at different scales and positions, thus balancing the global trend and local details of the signal. This feature makes it particularly suitable for processing laser radar full waveform data, because full waveform data usually contains non-stationary noise and target reflection pulse mixed characteristics. The core idea of denoising is to decompose the signal into different scales and frequency subbands through wavelet transform, and use the difference between signal and noise in the wavelet domain to suppress noise coefficients and retain effective signals, especially for the removal of non-stationary noise such as Gaussian white noise and pulse noise. The main steps are wavelet transform and threshold selection.

[0158] Discrete wavelet transform is used to decompose the signal into wavelet signals. Wavelet is a function that is small in time and highly localized in oscillation. Fourier transform decomposes the signal into infinite long sine and cosine, thus losing all the time position information, while wavelet transform is a time-localized scaling and shifting version that well preserves the time position information. Unlike Fourier transform, which always uses complex exponential basis functions, wavelet transform uses time-localized oscillation functions as basis functions. This basis function can be scaled and translated, and scaling at different scales corresponds to low and high frequencies of the signal, while translation contains time information. Wavelet analysis is more complex than Fourier analysis, because it must specify the basic function from which it is constructed.

[0159] Wavelet transform can extract the local features of the signal at different frequencies and times. Strong reflection pulses in laser radar full waveform correspond to larger values, while noise is characterized by widely distributed low amplitude coefficients. The actual wavelet basis function used is the Symlets wavelet. Symlets wavelet has tight support and approximate symmetry, which is suitable for processing phase-sensitive signals, so this wavelet is used to process the full waveform data of airborne laser radar.

[0160] Oscillation tailing is caused by the fact that the photomultiplier tube receives the signal and produces a rapid change in voltage. It is an active oscillation, and the voltage generated by the echo signal of the full waveform rises rapidly as a stimulus, affecting the output signal. The final output can be regarded as a dynamic coupling process of laser pulse and system response. The damping sinusoidal term of the mixed model represents the impact response of the system or environment, and the Gaussian term represents the instantaneous energy distribution of the laser emission. The convolution operation describes the dynamic coupling process of laser pulse and system response. The convolution model can more truly reflect the delay effect and frequency modulation in signal propagation. The amplitude and phase of the oscillation change dynamically with time during the propagation of the Gaussian pulse. Due to the nonlinear characteristics of the convolution model, the initial values of the Gaussian term and the damping sinusoidal term need to be estimated step by step.

[0161] In the estimation of the initial value of the parameter, the key sub-band is first selected, and after the wavelet transform, the full waveform echo peak is usually located in the low-frequency sub-band, and the damped oscillation is usually located in the medium-frequency sub-band, and after the selection verification, the low-frequency sub-band and part of the medium-frequency sub-band are selected to perform the inverse wavelet transform to obtain the reconstructed signal component. After the low-frequency sub-band is reconstructed, the Gaussian peak signal is contained, and after the part of the medium-frequency sub-band is reconstructed, most of the oscillation signals are contained. In the reconstructed low-frequency signal component and the medium-frequency signal component, the initial parameter estimation of the mixed model is performed. In this way, the resources are reasonably utilized, and the initial parameter estimation method is simplified. The mixed model contains a Gaussian function and a damped oscillation function, and the two are estimated respectively.

[0162] The initial parameter amplitude of the damped sinusoidal signal decays exponentially with time, the oscillation frequency is related to the periodicity of the signal, and the instantaneous phase should have a linear relationship with time. In the estimation, the real signal is first converted into an analytic signal through Hilbert transform, from which the instantaneous amplitude and phase are extracted. Hilbert transform is a linear time-invariant transform that converts a real-valued signal into a signal orthogonal to the original signal. The current data waveform is not convenient for analysis, and it is difficult to obtain the initial parameters of the current data waveform, so the complex signal, that is, the analytic signal, needs to be obtained. The analytic signal of a signal is a complex signal composed of the original signal and its Hilbert transform.

[0163] The embodiment of the application uses airborne laser radar to perform data experiments, part of the binary original data collected by the airborne laser radar is hexadecimal, the data converted from the original file of the airborne laser radar is decimal, Figure 3 is one of the waveforms collected by the airborne laser radar fixed-point test, corresponding to the full waveform data collected by one laser emission.

[0164] Figure 4 、 Figure 5 、 Figure 6 、 Figure 7 、 Figure 8 、 Figure 9 、 Figure 10 、 Figure 11 A series of pictures of the wavelet transformed waveform components, Figure 4 is the approximation coefficient component of the waveforms with the layer number 0, Figure 5 、 Figure 6 、 Figure 7 、 Figure 8 、 Figure 9 、 Figure 10 、 Figure 11 are respectively the detail coefficient components of the waveforms with the layer number 1 to the layer number 7.

[0165] Figure 12For wavelet threshold denoising waveform, the number of decomposition layers and the threshold value are adaptively selected by algorithm, and most of the noise is filtered out. Figure 2 For fitting optimization waveform, the full waveform data is fitted, the function parameters are extracted, and the target distance representation is stable and more accurate. Part of the calibration distance includes three columns, the first column is the target position extracted in the waveform, the second column is the intensity of the target echo, and the third column is the measured distance. Figure 1 For the technical process of the present application, it includes reading binary data, data conversion, three data blocks respectively, the same processing is carried out on each data block, including wavelet threshold denoising, extracting the intermediate value of wavelet coefficient, then being divided into two branches, one branch inputs wavelet transform layer 1, obtains overall characteristics, obtains the initial value of Gaussian function, the other branch inputs wavelet transform layer 2, performs initial fitting, obtains the initial value of damping function, then constructs a mixed model for the two initial values, performs global optimization by particle swarm, obtains target parameters, outputs three data, sorts and splices the data, and outputs the data.

[0166] The above examples are only used to illustrate the technical solutions of the present application, and not to limit them. Although the present application has been described in detail with reference to the foregoing examples, those skilled in the art should understand that the technical solutions recorded in the foregoing examples can still be modified, or some or all of the technical features can be replaced by equivalents, without changing the essence of the corresponding technical solutions out of the scope of the technical solutions of the embodiments of the present application.

Claims

1. A method for full waveform decomposition with oscillating signals in airborne lidar, characterized in that, The method comprises reading and converting original data, dividing the data into blocks, selecting wavelet threshold denoising for each block of data, constructing a hybrid model, selecting initial parameters of the hybrid model, performing particle swarm optimization to obtain optimized parameters, constructing output data of each block of data, splicing all output data, and completing full waveform decomposition. Constructing a mixture model including: ； wherein is a natural constant, is a Gaussian function amplitude, is an angular frequency, is a damped oscillation function amplitude, is a decay factor, is a Gaussian mean value, is a signal standard deviation, is a time scale, is a time.

2. The method according to claim 1, wherein, The airborne laser radar comprises an optical path system and a collection system, the optical path system and the collection system share a same square wave with a frequency of 5 kHz, the square wave controls laser emission of the optical path system, the square wave controls frequency collection of the collection system, the optical path system and the collection system are triggered at edges of the square wave, so that echo signals correspond to collected data one by one, the optical path system receives the echo signals, and the collected data are quantized voltage values which are negative numbers.

3. The method of claim 2, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, Reading and converting the original data comprises reading binary numbers in groups of 16 bits into decimal numbers, after reading is completed, inverting the decimal numbers bit by bit, converting negative numbers into positive numbers for processing, and storing all data in the form of tab character intervals; Dividing all data into blocks, numbering different blocks, distributing the blocks to processes of a central processing unit for processing, sorting the processed blocks according to numbers after processing of each block is completed, finally splicing and storing the data, extracting distance data from the original data, and taking the distance data as full waveform data of the airborne laser radar.

4. The method of claim 3, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, Wavelet threshold denoising comprises selecting an adaptive decomposition layer number, discrete wavelet transform, determining a wavelet threshold conforming to the system, and inverse wavelet transform; The determining the decomposition layer number comprises cutting and extracting the part of the distance data in front of the data, and determining the length of the collected data , the decomposition layer number is: ； In the formula, is the laser pulse width, is the acquisition frequency of the airborne laser radar system; The discrete wavelet transform comprises performing discrete wavelet transform on waveform data as a given signal: ； wherein is a detail coefficient obtained by a discrete wavelet transform, is a scale controlling the stretching of the wavelet, is a wavelet basis function, is a shift controlling the position of the wavelet basis, is a given signal.

5. The method of claim 4, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, Determining the wavelet threshold includes calculating an energy density of the detail coefficients of each layer in the decomposition level : ； ； wherein is the first layer of the detail coefficient, denotes the first detail coefficient in the layer, is the length of the first layer. Set , find the density difference that meets the following conditions : ； wherein is an empirical value, denotes the median, is the energy density of the layer. For the judgment condition of meeting From the corresponding first Estimate noise soft threshold from the sub-band of the discrete wavelet transform of the layer : ； Wavelet thresholding is: ； Performing inverse wavelet transform on detail coefficients higher than the wavelet threshold to reconstruct the given signal from which noise is filtered.

6. The method of claim 5, wherein the full waveform decomposition with oscillatory signals in an airborne lidar is characterized by, Selecting initial parameters of the hybrid model comprises selecting initial parameters of a Gaussian function and initial parameters of a damped oscillation function, performing multi-level wavelet decomposition on the given signal from which noise is filtered to extract sub-signals of different frequency bands, and the multi-level wavelet decomposition is: ； wherein is a low-pass filter for the given signal after filtering out noise, is a high-pass filter for the given signal after filtering out noise, represents a scale coefficient, is an approximated coefficient after filtering by a low-pass filter, is an approximated signal at scale after filtering by a low-pass filter, is a detail coefficient after filtering by a low-pass filter; is a discrete time index, is a summation index, is a set of integers; right Perform low-frequency convolution to obtain low-frequency information. Perform high-frequency convolution to obtain high-frequency information, then perform low-frequency convolution on the low-frequency information, while keeping the high-frequency information unchanged. After the decomposition, the detail coefficients of each sub-band from high frequency to low frequency are obtained. The low frequency component and the mid frequency component are selected and reconstructed by inverse wavelet transform to obtain the component where the Gaussian peak is located and the component where the oscillation signal is located. Selecting a region around an end point of a Gaussian peak as a starting point of an oscillation signal, and selecting a Gaussian item initial value as a parameter to extract a phase at an oscillation starting time; Correcting the initial parameters, and taking the corrected initial parameters as initial particles of a particle swarm optimization algorithm.

7. The method of claim 6, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, The initial parameters of the Gaussian function include the Gaussian function are: ； wherein is the Gaussian function argument, equivalent to the component in which the Gaussian peak is located; using the intensity of the low-frequency signal as , the center point of the low-frequency signal as , the half-height width of the low-frequency signal as : 。 8. The method of claim 7, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, The initial values of the damped oscillation function include , , and the initial phase , the damped oscillation function is: ； obtaining a resolved signal : ； wherein is the original signal, equivalent to the component in which the oscillating signal is located; frequency domain after Euler transformation is: ； wherein is a symbol function, is the original signal in the frequency domain; Set , In the time domain: ； wherein is an intermediate variable in the frequency domain, is an intermediate variable in the time domain, denotes a convolution operation, is the original signal in the time domain; find out : ； wherein is the limit symbol, represents Cauchy principal value integration.

9. The method of claim 8, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, Solve After, into the analytical signal formula, get ; is: ； wherein is a complex argument function; right After taking a logarithmic transformation of the envelope signal, a linear fit is performed: ； The slope of the linear fit result and the intercept , , are: ； ； is: 。 10. The method of claim 9, wherein the full waveform decomposition method with oscillating signal in an airborne lidar is characterized by, Particle swarm optimization includes randomly generating a group of particles, initializing particle position and particle velocity, calculating the fitting effect using mean square error , when the mean square error of the fitting result is the smallest, the solution of the current particle swarm optimization is the optimal solution, and the particle individual optimal solution and the global optimal solution at this time are recorded; , when the mean square error of the fitting result is the smallest, the solution of the current particle swarm optimization is the optimal solution, and the particle individual optimal solution and the global optimal solution at this time are recorded; The computation of the fitting result is as follows: The particle is updated as follows: ； ； wherein is the particle velocity of the is an inertia factor, and is an acceleration constant, denotes a random number in the interval is the individual best solution of the particle of the particle update, is the global best solution of the particle update, is the particle position of the After the particle swarm optimization is completed, the Gaussian function parameters and the damped oscillation function parameters are converged through multiple iterations to complete full waveform decomposition.

Images