A method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology

Abstract

The present invention discloses a method for adaptively identifying the thickness of black soil layers based on a variational modal decomposition method and ground-penetrating radar technology. Ground-penetrating radar technology can be used to directly measure the thickness of black soil layers in large areas of black soil farmland, automatically filter out ground-penetrating radar signal clutter, and correct the ground position. The technical solution of the present invention is based on a center frequency method and a residual method, automatically obtaining variational modal decomposition parameters applicable to black soil ground-penetrating radar signals, then performing a counting analysis on the modal components of the variational modal decomposition, comparing the spectrum of the original ground-penetrating radar signal with the actual black soil layer thickness, establishing an empirical model based on peak area ratio, and automatically screening the filtered center frequency and frequency sweep width that can be used to identify the thickness of the black soil layer. After frequency filtering, the black soil layer position is adaptively identified, enabling rapid acquisition of the black soil layer thickness and improving the accuracy and efficiency of identification.

Description

Technical Field

[0001] The present invention relates to the field of black soil layer thickness detection, and in particular to a method for adaptively identifying the thickness of black soil layers based on a variational mode decomposition method and ground penetrating radar technology. Background Art

[0002] Black soil is a precious arable land resource. Its properties are loose and porous, with good structure, a deep, mostly granular humus layer, a heavy, sticky texture, and a relatively uniform mechanical composition. This type of soil has good properties, high fertility, a high organic matter content, and a deep soil layer, resulting in high productivity and ideal for growing crops. It produces a large variety of crops, including soybeans, corn, rice, millet, sugar beets, and flax. The Northeast Black Soil Region is one of my country's key agricultural production areas, with a total arable land area of ​​approximately 35.84 million hectares. 2 , grain production accounts for about 1 / 4 of the country's total grain production and 1 / 3 of commercial grain, and is a "stabilizer" and "ballast" to ensure national food security. Soil thickness, as a physical indicator that characterizes soil fertility, quality, and erosion level, is an important material basis for soil fertility and crop growth. Due to soil erosion and long-term heavy use and neglect of maintenance, various soil degradation problems such as thinning of the black soil layer in black soil arable land have occurred. Therefore, quickly and accurately obtaining the thickness of the black soil layer of black soil arable land is of great significance for evaluating the quality and fertility of black soil arable land and guiding agricultural production.

[0003] Soil thickness has a significant impact on crops. Traditional methods for obtaining soil thickness include soil profiling, drilling, drilling, and penetrometer methods. These methods have the advantages of being accurate, intuitive, and providing comprehensive data sets for detecting black soil thickness. However, they often require a large amount of manpower and material resources, causing damage to the soil structure. Sampling is also difficult, costly, and inefficient. These methods are suitable for small areas but difficult to apply to large-scale continuous testing. They are often used as a means of verification. Summary of the Invention

[0004] This application provides a method for adaptively identifying the thickness of black soil layers based on variational mode decomposition and ground penetrating radar technology, which solves the technical problems that the existing ground penetrating radar system needs to identify underground information through grayscale images, the recognition subjectivity has a great influence, and the traditional black soil layer thickness identification method cannot be measured quickly, continuously and over a large range.

[0005] The present application provides a method for adaptively identifying the thickness of a black soil layer based on a variational mode decomposition method and ground penetrating radar technology, which is characterized by comprising the following steps:

[0006] S1. Determine the thickness of the black soil layer of typical black soil farmland using traditional methods;

[0007] S2. Use ground-penetrating radar to detect and collect data from cultivated black soil land to obtain original signals. The original signals include the time window of signal propagation in the soil and the amplitudes of different time windows under two-way travel time. During the ground-penetrating radar detection process, a signal is collected at a certain interval (the interval between each channel in the present invention is 4.48-4.67 cm). The time window-amplitude data of the total number of channels collected along the ground-penetrating radar route are combined to form the original radar signal of the profile.

[0008] S3, pre-processing the obtained original signal, determining the ground position by filtering and finding the zero time, and eliminating signal noise and sampling error;

[0009] S4. performing variational mode decomposition on the preprocessed original signal, using a center frequency method and a residual method to determine the optimal number of decomposition layers and an update step size for each channel of the original signal, to obtain an optimal decomposition mode model, and to retain the original signal characteristics to the greatest extent possible;

[0010] S5. Analyze the spectrum of the original signal to obtain the frequency distribution of the original signal, obtain representative modal components and their corresponding center frequencies based on the optimal number of decomposition levels and update step size, perform statistics on the modal components using a center frequency counting method, and obtain the frequency characteristics of the original signal by combining the continuous change and discrete distribution of the signal frequency;

[0011] S6. Based on the frequency characteristics and the actual black soil layer thickness determined in step S1, the average peak area ratio is determined, and the center frequency and frequency span of each data channel are automatically screened;

[0012] S7. Filter the original signal based on the filtered center frequency and frequency span, automatically identify the corresponding position of the black soil layer, and calculate the average black soil layer thickness.

[0013] Preferably, step S1 obtains the actual thickness of the black soil layer by combining the traditional profiling method, the soil drilling method with the field judgment method and the sampling laboratory test, and the specific steps include:

[0014] S1.1. Excavate the soil profile of cultivated black soil land and identify the location and thickness of the black soil layer according to the field black soil layer identification standard. Observe the soil structure, divide the soil layers, and collect soil samples. Collect 3-5 kg ​​of soil samples from each layer, and collect three ring knife samples from each layer. Subsequently, measure the soil bulk density in the laboratory.

[0015] S1.2. Drive at least three soil augers sequentially around the soil profile. After pulling them out, preliminarily determine the thickness of the black soil layer according to field judgment standards.

[0016] S1.3. After the soil samples were air-dried, their organic matter content was determined using the potassium dichromate-oil bath method. Referring to the criteria for determining the dark fertile topsoil layer in the Chinese soil classification system, the layer with a soil organic matter content greater than 10.3 g / kg was determined to be the black soil layer. The average actual black soil layer thickness was calculated in combination with the field soil profile survey and soil drilling measurement.

[0017] Preferably, step S2 includes the following steps:

[0018] S2.1. First, level the ground to reduce the impact of ground unevenness on the calculation of black soil thickness and ground-coupled antenna detection.

[0019] S2.2. Set sampling parameters, including the detection time window, offset, and gain. The time window is set according to the detection depth. The offset parameter is set so that the first positive wave peak is completely collected to facilitate subsequent ground location search. A relatively large value is used when setting the gain.

[0020] S2.3. Sample 2 to 3 m along the top of the soil profile to include the entire profile. Repeat the sampling three times to obtain the original GPR data signal.

[0021] Preferably, the pre-processing method in step S3 comprises the following steps:

[0022] S3.1. Construct a signal filter to perform filtering preprocessing on the sampled raw GPR data signal. According to frequency analysis, the noise is in the mid-high frequency part, so a high-frequency filter is used to filter out the noise.

[0023] S3.2. Calculate the thickness of the black soil layer starting from the ground, and construct a method to find the ground zero time. The ground position is the first negative wave peak, and the first negative wave peak is taken as the zero time.

[0024] S3.3. Add zero values ​​of the same length to each eliminated GPR data to make the data length on the profile the same.

[0025] Preferably, step S4 includes the following steps:

[0026] S4.1. Construct a center frequency method for adaptively obtaining the number of decomposition layers. Iterate each channel of the original signal separately, set the upper limit of iteration to 30, and iterate from 1 to 30. Stop the iteration when the difference between two adjacent center frequencies is less than 5% of the antenna frequency to obtain the optimal number of decomposition layers.

[0027] S4.2. Based on the optimal number of decomposition levels, construct a residual method for adaptively obtaining an update step size. Perform variational mode decomposition on each trace of the original signal and iterate separately, with an update step size ranging from 0.01 to 1. Iterate from 0.01 to 1. Calculate the difference between the superposition value of the decomposed modal amplitudes and the original signal using the optimal number of decomposition levels. The update step size that minimizes the difference is the update step size for that trace.

[0028] S4.3. Based on the obtained optimal decomposition layer number and update step size, variational mode decomposition is performed on each signal to obtain different decomposition modes for each signal on the soil profile.

[0029] Preferably, step S5 includes the following steps:

[0030] S5.1. Use fast Fourier transform to calculate the frequency distribution of the original signal, superimpose the spectra of all traces of the profile, and preliminarily obtain the frequency distribution of the signal;

[0031] S5.2. Based on the frequency distribution range of the profile, construct a method for analyzing the decomposed modes of the profile. Perform statistical analysis on the decomposed modes of all traces within the profile. Then, divide the frequency range into different intervals at intervals of 5% of the antenna frequency and calculate the number of modal components in each interval.

[0032] S5.3. Based on the profile frequency distribution and the obtained modal component center frequency distribution, analyze the profile frequency characteristics, take the frequency peak obtained in step S5.1 as the frequency spectrum peak, find the interval to which the frequency peak in step S5.2 belongs, which is generally the maximum value of the center frequency count, and take this interval as the frequency spectrum peak interval.

[0033] Preferably, step S6 includes the following steps:

[0034] S6.1. Based on the modal center frequency count value and frequency spectrum peak interval obtained in step 5, partition the count value and take the first maximum value after the frequency spectrum peak interval as the frequency sub-peak value and its interval;

[0035] S6.2. Based on the frequency sub-peak value obtained in step S6.1, construct an empirical model for the spectrum area share. With this value as the center frequency and 0.5% of the antenna center frequency as the interval, iterate forward and backward from 0. Combined with the actual black soil layer thickness, the iteration stops when the predicted black soil layer thickness after filtering is close to the actual thickness, and the profile spectrum area share is obtained.

[0036] S6.3. Based on the profile spectral area ratios obtained in step S6.2, calculate the average spectral area ratios of all profiles to obtain an empirical model of spectral area ratios;

[0037] S6.4. Based on the empirical value of the spectrum area ratio obtained in step S6.3, use the frequency sub-peak obtained in step S6.2 to iterate the profile forward and backward from 0 with an interval of 0.5% of the antenna center frequency, and calculate the spectrum area ratio of each iteration. If the spectrum area ratio is greater than the empirical modulus of the spectrum area ratio, stop the iteration and filter to obtain the frequency span.

[0038] Preferably, step S7 includes the following steps:

[0039] S7.1. Filter the original signal based on the filtered center frequency and frequency span obtained in step S6 to obtain a signal for determining the position of the black soil layer;

[0040] S7.2. Based on the obtained black soil layer position signal, perform envelope processing on the signal, use a smooth envelope curve, take the second trough position of the upper envelope as the black soil layer position, and repeat the operation for all traces of the profile;

[0041] S7.3. Based on the positions of the black soil layers obtained for all the sections of the profile, calculate the average value to obtain the thickness of the black soil layer of the profile.

[0042] Preferably, step S4.1 includes the following steps:

[0043] S4.1.1. Assume that the multi-component signal is The modal components of limited bandwidth are composed of , the center frequency of each intrinsic mode function IMF is , the constraint condition is that the modal sum is equal to the input signal, which is obtained by Hilbert transform The analytical signal and its unilateral spectrum are calculated by using the operator Multiply, The center band is modulated to the corresponding baseband:

[0044]

[0045] in is the Dirac impulse function; ∗ denotes convolution; is an imaginary unit; is the angular frequency; is the natural index;

[0046] S4.1.2. Calculate the square norm of the demodulation gradient , and estimate the bandwidth of each mode component. The process is shown in the following formula:

[0047]

[0048] In the above formula, , Represents the decomposed IMF components, Represents the center frequency of each component; is the Dirac function; * is the convolution operator; is the original signal;

[0049] S4.1.3. Introduction of Lagrange multipliers and the second-order penalty factor , transforming the constrained variational problem into an unconstrained variational problem, the extended Lagrangian expression is as follows:

[0050]

[0051] S4.1.4. Use the alternating direction multiplier method to continuously update each component and its center frequency, and finally obtain the saddle point of the unconstrained model, which is the optimal solution to the original problem. All components can be obtained in the frequency domain space by the following formula:

[0052]

[0053] S4.1.5, yes After the Wiener filter, the algorithm re-estimates the centroid frequency based on the centroid of the power spectrum of each component. The specific process is as follows:

[0054] a) , , and ;

[0055] b) Cycle: ;

[0056] c) When updated according to S4.1.4 ;

[0057] d) Update according to the following formula ;

[0058]

[0059] e) Update according to the following formula ;

[0060]

[0061] Where: Tolerance to noise, meeting signal decomposition fidelity requirements; and Corresponding respectively and Fourier transform of

[0062] f) Repeat steps a to f until the iteration stop condition is met. The stop condition is:

[0063]

[0064] After stopping the iteration, we finally get IMF components and center frequencies.

[0065] As a preference, step S4.2 constructs a residual method for adaptively obtaining the update step size, and uses the residual method to determine the update step size. The steps for value include:

[0066] 4.2.1、Set the number of decomposition levels to ,set up Iterate from 0.01 to 1, and the IMF components obtained after each decomposition are , the total value of all component amplitudes is recorded as , as shown below:

[0067]

[0068] 4.2.2、 The absolute value of the difference from the original signal is recorded as the residual , the total residual values ​​obtained after iteration are , and record the minimum value as , as shown below:

[0069]

[0070]

[0071] in, The corresponding update step is the update step of the final decomposition, recorded as .

[0072] The technical solution provided by this application has at least the following technical effects or advantages:

[0073] 1. The technical solution of this invention uses ground-penetrating radar to directly measure soil over large areas, automatically filtering out signal clutter and correcting ground position. This enables non-destructive, continuous, and rapid detection of black soil thickness, efficiently obtaining the thickness of the black soil layer. This overcomes the shortcomings of existing ground-penetrating radar systems, which rely on grayscale images to identify underground information, resulting in significant subjectivity, and the inability of traditional black soil thickness identification methods to rapidly and continuously measure over large areas.

[0074] 2. The technical solution of the present invention is based on the center frequency method and the residual method to automatically obtain the variational modal decomposition parameters suitable for the black soil ground penetrating radar signal, then count and analyze the modal components of the variational modal decomposition, compare the frequency spectrum of the original ground penetrating radar signal and the actual black soil layer thickness, establish an empirical model based on the peak area ratio, and automatically screen the screened center frequency and frequency sweep width that can be used to identify the thickness of the black soil layer. After frequency filtering, the position of the black soil layer is adaptively identified, which can realize the accurate acquisition of the black soil layer thickness and improve the accuracy and efficiency of recognition. BRIEF DESCRIPTION OF THE DRAWINGS

[0075] Figure 1 It is a schematic diagram of model construction and actual measurement during the specific implementation of the present invention.

[0076] Figure 2 This is a rendering of the thickness of the black soil layer actually measured in the field during the specific implementation of the present invention;

[0077] Figure 3 This is the effect diagram after variational mode decomposition of a single channel of ground penetrating radar data;

[0078] Figure 4 This is the effect diagram after the iteration of the residual method;

[0079] Figure 5 This is the effect diagram of the center frequency counts after the variational mode decomposition of all traces of the sample points;

[0080] Figure 6 This is the effect diagram of frequency sweep width screening for a single channel of a sample point;

[0081] Figure 7 This is a comparison chart of the black soil layer thickness obtained by the traditional profile survey method and the black soil layer thickness predicted by the present invention. DETAILED DESCRIPTION

[0082] The specific embodiments of the present invention are described below to facilitate understanding of the present invention by those skilled in the art. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, as long as various changes are within the spirit and scope of the present invention as defined and determined by the appended claims, these changes are obvious, and all inventions and creations utilizing the concepts of the present invention are protected.

[0083] like Figure 1 As shown, the method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology provided in this embodiment includes the following steps:

[0084] S1. Determine the actual thickness of the black soil layer using traditional profile survey and soil drilling methods combined with field black soil layer determination criteria and laboratory soil testing after sampling. The specific operations in this embodiment include:

[0085] S1.1. Excavate a soil profile from a field black soil farmland, typically 2 m long, 1.5 m wide, and 1.2 m deep. Observe the profile and, according to field black soil layer identification criteria, identify the location and thickness of the black soil layer. Observe the soil structure and divide the soil into layers. Soil samples are collected based on the layers, weighing 3-5 kg ​​per layer. Three ring cutter samples are collected from each layer. The soil bulk density is then measured in the laboratory.

[0086] S1.2. Drive three or more soil augers sequentially around the profile and, after pulling them out, preliminarily determine the thickness of the black soil layer according to field judgment standards;

[0087] S1.3. Field criteria for identifying black soil layers are: 1) Lightness and chroma of the soil under moist conditions ≤ 3; 2) The soil structure is granular, with small angular blocks, and the soil layer is not hard, with less than half the volume of the rock;

[0088] S1.4. Bring the soil and samples collected in the field back to the laboratory and dry them. Then use the potassium dichromate-oil bath method to determine their organic matter content. The layer with a soil organic matter content greater than 10.5 g / kg is determined to be the black soil layer. Combine the profile measurement and soil drilling measurement to calculate the average actual black soil layer thickness, such as Figure 2 As shown, the yellow line is the lower boundary of the black soil layer.

[0089] S2. Use ground penetrating radar to detect and sample black soil farmland to obtain original signals. This embodiment specifically includes the following steps:

[0090] S2.1. Before conducting ground penetrating radar detection, level the ground to reduce the impact of ground unevenness on the calculation of black soil thickness and ground-coupled antenna detection;

[0091] S2.2 Before ground penetrating radar detection, it is necessary to set the channel settings, including detection time window, offset and gain. The time window is set according to the detection depth. Generally, the detection depth D is selected as 1.5 times the target depth, as shown in the following formula:

[0092] ;

[0093] D is the expected detection depth, in meters. V is the average radar wave velocity in the formation medium, in meters per nanosecond. W is the sampling window, in nanoseconds. A larger window increases the detection depth and reduces the resolution.

[0094] Depth is typically calculated by converting the time window into depth using the dielectric constant, estimating the depth of each boundary and the thickness of the generating layer. The dielectric constant is the ability of a substance to retain an electric charge. Its magnitude determines the medium's ability to absorb or reflect electromagnetic waves, and ranges from 1 (air) to 81 (water). Its calculation formula is shown below:

[0095]

[0096] Convert the above formula into: ;

[0097] When calculating the two-way electromagnetic wave, the depth (H, m) is calculated as follows:

[0098]

[0099] Where ε is the dielectric constant, c is the speed of light (3.00×108m / s), and v (m / s) is the speed of electromagnetic waves propagating in the medium.

[0100] Since the ground-coupled antenna is close to the ground, the offset is set to ensure that the air wave, that is, the first positive wave peak, is collected completely to facilitate subsequent search for the ground position. Since black soil severely reduces the amplitude, a larger value is used when setting the gain.

[0101] S2.3. Based on the sampling parameters set in S2.2, sample 2-3 m along the top of the profile, completely encompassing the entire profile. Repeat sampling three times to obtain raw GPR data, which includes the time window (time window) during which the signal propagates in the soil and the signal strength (amplitude) of each time window under two-way travel time. During the GPR detection process, the machine collects a signal at regular intervals (the interval between each channel in this invention is 4.48-4.67 cm). The time window-amplitude data of the total number of channels collected along the GPR route are combined to form the original radar signal of the profile.

[0102] S3. Preprocess the original signal data to eliminate signal noise and sampling errors by filtering and finding the ground position at time 0. This specifically includes the following steps:

[0103] S3.1. Construct a signal filter. Since the sampling signal of the black soil ground penetrating radar is severely attenuated, a larger gain parameter needs to be set to improve the readability of the grayscale image during sampling, which generates a large amount of machine noise. Therefore, the sampled signal needs to be pre-processed by filtering. According to frequency analysis, the noise is the mid-high frequency part of the frequency, so a high-frequency filter is used to filter out the noise part.

[0104] S3.2. After filtering out noise using the high-frequency filter constructed in S3.1, calculate the thickness of the black soil layer starting from the ground, and construct a method for finding the ground time 0. The ground position is the first negative wave peak, and the first negative wave peak is taken as the time 0.

[0105] S3.3. After filtering out noise and finding the ground zero moment based on S3.1 and S3.2, each channel of eliminated ground penetrating radar data needs to be added with zero values ​​of the same length to make the data length on the profile the same.

[0106] S4. Perform variational modal decomposition on the original signal, using the center frequency method and residual method to determine the optimal decomposition layer number and update step size for each channel of the original signal, obtain the optimal decomposition modal model, and retain the original signal characteristics to the greatest extent possible. The specific steps include:

[0107] S4.1. Construct a center frequency method that adaptively obtains the number of decomposition layers and iterates each channel of the original signal. The specific process of variational mode decomposition (VMD) can be understood as the optimal solution to the variational problem and can be converted into the construction and solution of the variational problem. The specific construction steps are as follows:

[0108] S4.1.1. Assume that the multi-component signal is The modal components of limited bandwidth are composed of , the center frequency of each intrinsic mode function (IMF) is , the constraint condition is that the modal sum is equal to the input signal, which is obtained by Hilbert transform The analytical signal and its unilateral spectrum are calculated by using the operator Multiply, The center band is modulated to the corresponding baseband:

[0109]

[0110] in is the Dirac impulse function; ∗ denotes convolution; is an imaginary unit; is the angular frequency; is the natural index;

[0111] S4.1.2. Calculate the square norm of the demodulation gradient , and estimate the bandwidth of each mode component. The process is shown in the following formula:

[0112]

[0113] In the above formula, , Represents the decomposed IMF components, Represents the center frequency of each component; is the Dirac function; * is the convolution operator; is the original signal;

[0114] S4.1.3. In order to find the optimal solution to the constrained variational problem, we first introduce the Lagrange multiplier and the second-order penalty factor , transforming the constrained variational problem into an unconstrained variational problem. The second-order penalty factor It can ensure the accuracy of signal reconstruction in Gaussian noise environment. Lagrange multiplier It can ensure that the strictness of the constraints is maintained. The extended Lagrangian expression is as follows:

[0115]

[0116] S4.1.4. Use the alternating direction method of multipliers (ADMM) to continuously update the components and their center frequencies, and finally obtain the saddle point of the unconstrained model, which is the optimal solution to the original problem. All components can be obtained in the frequency domain space by the following formula:

[0117]

[0118] S4.1.5, yes After the Wiener filter, the algorithm re-estimates the centroid frequency based on the centroid of the power spectrum of each component. The specific process is as follows:

[0119] a) , , and ;

[0120] b) Cycle: ;

[0121] c) When updated according to S4.1.4 ;

[0122] d) Update based on the data ;

[0123]

[0124] e) Update according to the following formula ;

[0125]

[0126] Where: Tolerance to noise, meeting signal decomposition fidelity requirements; and Corresponding respectively and The Fourier transform of .

[0127] f) Repeat steps a to f until the iteration stopping condition is met. The stopping condition is shown in the formula;

[0128]

[0129] After stopping the iteration, we finally get IMF components and center frequencies. Usually, when the number of decomposition layers is greater than 30, the center frequencies of adjacent decomposition modes are close and the amplitudes are small. The iteration limit is set to 30, and iterations are started from 1 to 30. The stopping condition is that the difference between two adjacent center frequencies is less than 5% of the antenna frequency. For example, Figure 3 shown.

[0130] In this embodiment, the central frequency method is used to determine The specific steps are as follows:

[0131] ① Determine the center frequency interval. This method sets the interval to 5% of the GPR antenna and denotes the center frequency of the antenna as , the center frequency interval is recorded as .

[0132] ② After setting the center frequency interval, decompose the signal according to the variational mode decomposition method steps and set Iterate from 3 to 20, and get components and center frequency values , when the difference between the last three digits of the center frequency is less than The iteration stops when . The judgment formula is as follows:

[0133]

[0134]

[0135] Where: ; The number of decomposition layers after iterative use of the center frequency method is , for IMF components, for The center frequency corresponding to the IMF component is shown in Table 1. The iterative steps of obtaining the decomposition level K using the center frequency method are shown in Table 1.

[0136] Table 1

[0137]

[0138] S4.2. Based on the number of decomposition layers obtained in S4.1, a residual method for adaptively obtaining the update step size is constructed. Each trace is subjected to variational mode decomposition and iterated separately. The update step size is usually 0.01~1, and it is iterated from 0.01 to 1. The difference between the superposition value of the decomposed modal amplitude and the original signal is calculated based on the number of decomposition layers obtained. The update step size with the smallest difference is the update step size of the trace. The update step size is determined using the residual method. The specific steps of the value are as follows:

[0139] 4.2.1、Set the number of decomposition levels to Iterate from 0.01 to 1, and the IMF components obtained after each decomposition are , the total value of all component amplitudes is recorded as , as shown below:

[0140]

[0141] 4.2.2、 The absolute value of the difference from the original signal is recorded as the residual , the total residual values ​​obtained after iteration are , and record the minimum value as , as shown below:

[0142]

[0143]

[0144] in, The corresponding update step is the update step of the final decomposition, recorded as ,like Figure 4 shown.

[0145] S4.3, based on the decomposition layer number and update step size obtained in S4.1 and S4.2, perform variational mode decomposition on each signal channel to obtain different decomposition modes for each signal channel on the profile. The signal data obtained by the ground penetrating radar is a three-dimensional data of channel number, time window and amplitude. Assume that the signal data is Channel, each signal (time window-amplitude two-dimensional data) Perform variational mode decomposition and use the center frequency method to determine the number of decomposition layers for each channel value.

[0146] S5. Analyze the spectrum of the original signal to obtain the original signal frequency distribution. After obtaining representative modal components and their corresponding center frequencies based on the number of decomposition levels and update step size obtained in step S4, use the center frequency counting method to count the modal components. Combined with the continuous change and discrete distribution of the signal frequency, the frequency characteristics of the original signal are obtained. Specifically, the following steps are included:

[0147] S5.1. Use fast Fourier transform to calculate the frequency distribution of the original signal, superimpose the spectrum of all channels of the profile, and preliminarily obtain the frequency distribution of the signal; determine the number of decomposition layers and update step size After that, the original signal Perform variational mode decomposition and calculate using fast Fourier transform The frequency distribution of , the peak value is , the frequency range is separated by the center frequency Divide, a total of Round off the interval. To record all The component of , its corresponding center frequency is:

[0148] ;

[0149] S5.2. Based on the profile frequency range obtained in S5.1, a method for analyzing the profile decomposition mode is constructed. The decomposition modes of all traces in the profile are statistically analyzed. The frequency range is then divided into different intervals at intervals of 5% of the antenna frequency as follows:

[0150] ;in is the floor function.

[0151] The center frequency intervals after division are , all center frequencies correspond Classify and calculate the number of IMF components in each interval: .

[0152] S5.3. Based on the cross-sectional frequency distribution obtained in S5.1 and the modal component center frequency distribution obtained in S5.2, analyze the cross-sectional frequency characteristics and convert the frequency distribution obtained by spectrum analysis into Corresponding to the frequency of the component count, the frequency peak in S5.1 is taken as the frequency spectrum peak, and the interval where the peak in S5.2 is located is found, which is generally also the maximum value of the center frequency count. This interval is taken as the frequency spectrum peak interval and recorded with The most recent peak count value The central frequency peak range of the variational mode decomposition is ,mark The peak value of the spectrum is , recorded as the secondary peak, the corresponding center frequency range is , recorded as the sub-peak interval, and the median of the range is recorded as the center frequency after screening ,like Figure 5 shown.

[0153] S6. Based on the frequency characteristics obtained in step S5 and the thickness of the black soil layer obtained on site, after determining the average peak area ratio, the center frequency and frequency span of each data channel are automatically screened. Specifically, the steps include:

[0154] S6.1. Based on the modal center frequency count values ​​and frequency spectrum peak intervals obtained in S5, partition the count values ​​and take the first maximum value after the frequency spectrum peak interval as the frequency sub-peak value and its interval;

[0155] S6.2. Based on the frequency sub-peak obtained in S6.1, construct an empirical model of spectrum area ratio. Take this value as the center frequency and 0.5% of the antenna center frequency as the interval, iterate forward and backward from 0, and combine it with the actual black soil layer thickness. The stopping condition is that the predicted black soil layer thickness after filtering is close to the actual thickness, and obtain the profile spectrum area ratio; obtain the center frequency after screening Afterwards, each signal Filtering is performed to filter the original signal into the upper and lower filter ranges. This range is the frequency span after filtering, which is recorded as ,filter The steps are as follows:

[0156] by As the center, the upper and lower limits of frequency are respectively Iterate towards both ends with an iteration step length of , the upper and lower limits of iteration are shown as follows: The spectrum of and Bandpass filtering is performed for the upper and lower frequencies, and the filtered signal is recorded as .

[0157]

[0158]

[0159] S6.3. Based on the profile spectrum area percentages obtained in S6.2, calculate the average spectrum area percentages of all profiles to obtain an empirical model for spectrum area percentages. Statistical analysis shows that the area of ​​the spectrum used to identify the thickness of the black soil layer is typically approximately 10.4% of the total area. Therefore, the iteration is terminated when the spectrum area exceeds 10.4% of the original spectrum area.

[0160] S6.4. Based on the empirical value of spectrum area ratio obtained in S6.3, use the frequency sub-peak obtained in S6.2 to iterate the profile forward and backward at intervals of 0.5% of the antenna center frequency, and calculate the spectrum area ratio of each iteration. The stopping condition is that the spectrum area ratio is greater than the empirical modulus of the spectrum area ratio, and the frequency span is obtained by screening. The calculation formula is as follows:

[0161] ; The effect of frequency sweep width of single-channel screening of sample points is as follows Figure 6 shown.

[0162] S7, based on the filtered center frequency and frequency span obtained in step S6, the original signal is filtered, the corresponding position of the black soil layer is automatically identified, and the average black soil layer thickness is calculated. The specific method of filtering the original signal based on the filtered center frequency and frequency span in this step includes the following steps:

[0163] S7.1, based on the center frequency and frequency span after screening obtained in S6, by Perform band-pass filtering to obtain a signal for determining the location of the black soil layer;

[0164] S7.2, based on the black soil layer position signal obtained in S7.1, the filtered signal Perform the envelope operation using a smooth envelope curve. The second trough of the upper envelope is taken as the location of the black soil layer. Repeat this operation for all traces in the profile. Generally, at least 20 traces are required for center frequency counting.

[0165] S7.3. Based on S7.2, obtain the position of the black soil layer in all sections of the profile, calculate the average value, and obtain the thickness of the black soil layer in the profile.

[0166] The above steps are performed on each radar data channel to obtain the position of the black soil layer on the profile, thereby achieving the purpose of predicting the thickness of the black soil layer along the radar measurement route. In one embodiment of the present invention, after performing ground-penetrating radar detection on black soil cultivated land, the predicted black soil layer thickness is highly accurate. The sampling method, ground-penetrating radar channel setting, and data processing must be consistent with this method. Preprocessing can address the problem of ground unevenness, find the ground position, and accurately predict the black soil layer thickness.

[0167] The comparison of the black soil layer thickness obtained by the traditional soil profile survey method and the black soil layer thickness predicted by the present invention is as follows: Figure 7 shown.

[0168] In summary, the present invention collects ground-penetrating radar data and other data, adopts a variety of data processing methods, and integrates radar data and field data features to conduct large-scale detection of the black soil layer thickness of black soil farmland. This can improve the efficiency and accuracy of obtaining the black soil layer thickness, and solves the problem that existing black soil layer thickness detection methods are difficult to achieve large-scale rapid and non-destructive detection.

[0169] The above is only a preferred embodiment of the present invention. It should be pointed out that for ordinary technicians in this technical field, several improvements can be made without departing from the principles of the present invention. These improvements should also be regarded as the scope of protection of the present invention.

Claims

1. A method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology, characterized in that: The steps include: S1. Determine the thickness of the black soil layer of typical black soil farmland using traditional methods, including traditional profiling method, soil drilling method combined with field judgment method and sampling laboratory testing to obtain the actual black soil layer thickness; S2. Use a ground-penetrating radar to detect and collect data from cultivated black soil to obtain an original signal. The original signal includes a time window for signal propagation in the soil and the amplitudes of different time windows under two-way travel time. During the ground-penetrating radar detection process, a signal is collected at each set interval. The time window-amplitude data of the total number of channels collected along the ground-penetrating radar route are combined to form the original radar signal of the profile. S3. Preprocess each acquired original signal, determine the ground position by filtering and finding the zero time, and eliminate signal noise and sampling error; S4. Perform variational modal decomposition on the preprocessed original signal, use the center frequency method and the residual method to determine the optimal decomposition layer number and update step size for each channel of the original signal, obtain the optimal decomposition modal model, retain the original signal characteristics to the maximum extent, and construct a residual method for adaptively obtaining the update step size based on the optimal decomposition layer number. Perform variational modal decomposition on each channel of the original signal, and iterate separately, with an update step size of 0.01 to 1, starting from 0.01 and iterating toward 1. Calculate the difference between the superposition value of the decomposed modal amplitude and the original signal with the optimal decomposition layer number, where the update step size with the smallest difference is the update step size of the channel; S5. Analyze the spectrum of the original signal to obtain the frequency distribution of the original signal, obtain representative modal components and their corresponding center frequencies based on the optimal number of decomposition levels and update step size, perform statistics on the modal components using a center frequency counting method, and obtain the frequency characteristics of the original signal by combining the continuous change and discrete distribution of the signal frequency; S6. Based on the frequency characteristics and the black soil layer thickness determined in step S1, the average peak area ratio is determined, and the center frequency and frequency span of each data channel are automatically screened; S7. Filter the original signal based on the filtered center frequency and frequency span, automatically identify the corresponding position of the black soil layer, and calculate the average black soil layer thickness.

2. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 1 is characterized in that: The step S1 comprises: S1.

1. Excavate the soil profile of cultivated black soil land and identify the location and thickness of the black soil layer according to the field black soil layer identification standard. Observe the soil structure, divide the soil layers, and collect soil samples. Collect 3-5 kg ​​of soil samples from each layer, and collect three ring knife samples from each layer. Subsequently, measure the soil bulk density in the laboratory. S1.

2. Drive at least three soil augers sequentially around the soil profile. After pulling them out, preliminarily determine the thickness of the black soil layer according to field judgment standards. S1.

3. After the soil samples were air-dried, their organic matter content was determined using the potassium dichromate-oil bath method. Referring to the criteria for determining the dark fertile topsoil layer in the Chinese soil classification system, the layer with a soil organic matter content greater than 10.3 g / kg was determined to be the black soil layer. The average actual black soil layer thickness was calculated in combination with the field soil profile survey and soil drilling measurement.

3. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 1 is characterized in that: Step S2 includes the following steps: S2.

1. First, level the ground to reduce the impact of ground unevenness on the calculation of black soil thickness and ground-coupled antenna detection. S2.

2. Set sampling parameters, including the detection time window, offset, and gain. The time window is set according to the detection depth. The offset parameter is set so that the first positive wave peak is completely collected to facilitate subsequent ground location search. A relatively large value is used when setting the gain. S2.

3. Sample 2 to 3 m along the top of the soil profile to include the entire profile. Repeat the sampling three times to obtain the original GPR data signal.

4. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 1 is characterized in that: The pre-processing method in step S3 includes the following steps: S3.

1. Construct a signal filter to perform filtering preprocessing on the sampled raw GPR data signal. According to frequency analysis, the noise is in the mid-high frequency part, so a high-frequency filter is used to filter out the noise. S3.

2. Calculate the thickness of the black soil layer starting from the ground, and construct a method to find the ground zero time. The ground position is the first negative wave peak, and the first negative wave peak is taken as the zero time. S3.

3. Add zero values ​​of the same length to each eliminated GPR data to make the data length on the profile the same.

5. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 4 is characterized in that: Step S4 includes the following steps: S4.

1. Construct a center frequency method for adaptively obtaining the number of decomposition layers. Iterate each channel of the original signal separately, set the upper limit of iteration to 30, and iterate from 1 to 30. Stop the iteration when the difference between two adjacent center frequencies is less than 5% of the antenna frequency to obtain the optimal number of decomposition layers. S4.

2. Calculate the difference between the superposition value of the decomposed modal amplitude and the original signal. The update step size with the smallest difference is the update step size of the trace. S4.

3. Based on the obtained optimal decomposition layer number and update step size, variational mode decomposition is performed on each signal to obtain different decomposition modes for each signal on the soil profile.

6. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 5 is characterized in that: Step S5 includes the following steps: S5.

1. Use fast Fourier transform to calculate the frequency distribution of the original signal, superimpose the spectra of all traces of the profile, and preliminarily obtain the frequency distribution of the signal; S5.

2. Based on the frequency distribution range of the profile, construct a method for analyzing the decomposed modes of the profile. Perform statistical analysis on the decomposed modes of all traces within the profile. Then, divide the frequency range into different intervals at intervals of 5% of the antenna frequency and calculate the number of modal components in each interval. S5.

3. Based on the profile frequency distribution and the obtained modal component center frequency distribution, analyze the profile frequency characteristics, take the frequency peak obtained in step S5.1 as the frequency spectrum peak, find the interval to which the frequency spectrum peak belongs in the center frequency count obtained in step S5.2, and take this interval as the frequency spectrum peak interval.

7. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 6, characterized in that: Step S6 includes the following steps: S6.

1. Based on the modal center frequency count value and frequency spectrum peak interval obtained in step 5, partition the count value and take the first maximum value after the frequency spectrum peak interval as the frequency sub-peak value and its interval; S6.

2. Based on the frequency sub-peak value obtained in step S6.1, construct an empirical model for the spectrum area share. With this value as the center frequency and 0.5% of the antenna center frequency as the interval, iterate forward and backward from 0. Combined with the actual black soil layer thickness, the iteration stops when the predicted black soil layer thickness after filtering is close to the actual thickness, and the profile spectrum area share is obtained. S6.

3. Based on the profile spectral area ratios obtained in step S6.2, calculate the average spectral area ratios of all profiles to obtain an empirical model of spectral area ratios; S6.

4. Based on the empirical value of the spectrum area ratio obtained in step S6.3, use the frequency sub-peak obtained in step S6.2 to iterate the profile forward and backward from 0 with an interval of 0.5% of the antenna center frequency, and calculate the spectrum area ratio of each iteration. If the spectrum area ratio is greater than the empirical modulus of the spectrum area ratio, stop the iteration and filter to obtain the frequency span.

8. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 7 is characterized in that: Step S7 includes the following steps: S7.

1. Filter the original signal based on the filtered center frequency and frequency span obtained in step S6 to obtain a signal for determining the position of the black soil layer; S7.

2. Based on the obtained black soil layer position signal, perform envelope processing on the signal, use a smooth envelope curve, take the second trough position of the upper envelope as the black soil layer position, and repeat the operation for all traces of the profile; S7.

3. Based on the positions of the black soil layers obtained for all the sections of the profile, calculate the average value to obtain the thickness of the black soil layer of the profile.

9. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 5, characterized in that: Step S4.1 includes the following steps: S4.1.

1. Assume that the multi-component signal is The modal components of limited bandwidth are composed of , the center frequency of each intrinsic mode function IMF is , the constraint condition is that the modal sum is equal to the input signal, which is obtained by Hilbert transform The analytical signal and its unilateral spectrum are calculated by using the operator Multiply, The center band is modulated to the corresponding baseband: in is the Dirac impulse function; ∗ denotes convolution; is an imaginary unit; is the angular frequency; is the natural index; S4.1.

2. Calculate the square norm of the demodulation gradient , and estimate the bandwidth of each mode component. The process is shown in the following formula: In the above formula, , Represents the decomposed IMF components, Represents the center frequency of each component; is the Dirac function; * is the convolution operator; is the original signal; S4.1.

3. Introduction of Lagrange multipliers and the second-order penalty factor , transforming the constrained variational problem into an unconstrained variational problem, the extended Lagrangian expression is as follows: S4.1.

4. Use the alternating direction multiplier method to continuously update each component and its center frequency, and finally obtain the saddle point of the unconstrained model, which is the optimal solution to the original problem. All components can be obtained in the frequency domain space by the following formula: S4.1.5, yes After the Wiener filter, the algorithm re-estimates the centroid frequency based on the centroid of the power spectrum of each component. The specific process is as follows: a) , , and ; b) Cycle: ; c) When updated according to S4.1.4 ; d) Update according to the following formula ; e) Update according to the following formula ; Where: Tolerance to noise, meeting signal decomposition fidelity requirements; and Corresponding respectively and Fourier transform of f) Repeat steps a to f until the iteration stop condition is met. The stop condition is: After stopping the iteration, we finally get IMF components and center frequencies.

10. The method for adaptively identifying the thickness of black soil layer based on variational mode decomposition method and ground penetrating radar technology according to claim 9, characterized in that: Step S4.2: Construct a residual method to adaptively obtain the update step size, and use the residual method to determine the update step size. The steps for value include: 4.2.1、Set the number of decomposition levels to ,set up Iterate from 0.01 to 1, and the IMF components obtained after each decomposition are , the total value of all component amplitudes is recorded as , as shown below: 4.2.2、 The absolute value of the difference from the original signal is recorded as the residual , the total residual values ​​obtained after iteration are , and record the minimum value as , as shown below: in, The corresponding update step is the update step of the final decomposition, recorded as .

Images