A method for determining the continuous spatial variation of the gley layer based on variational mode decomposition and ground-penetrating radar technology

Abstract

This invention relates to the field of soil barrier layer detection. It discloses a method for determining the spatially continuous variation of gleyed layers based on variational mode decomposition (VMD) and ground-penetrating radar (GPR) technology. The technical solution of this invention is based on VMD, obtaining low-frequency mode components through the center frequency, constructing a Pearson correlation matrix for the original signal and each mode component, extracting the optimal combination of mode components, and combining the grayscale image of the optimal combination of mode components with the location of the gleyed layer using a soil drill to extract the characteristic wave information of the gleyed layer. The location of the gleyed layer on the radar survey line is then extracted. Using the location of the gleyed layer on the radar survey line, the location of the gleyed layer in the test area is obtained through inverse range weighted interpolation. Finally, the R-value of the linear fitting is used... 2 This invention demonstrates the accuracy of obtaining the gill layer thickness along the survey line. It enables rapid acquisition of the location and thickness of the gill layer on a large scale, improving identification efficiency while ensuring accuracy.

Description

Technical Field

[0001] This invention relates to the field of soil barrier layer detection, specifically to a method for determining the spatial continuous changes of the gley layer based on variational mode decomposition and ground penetrating radar technology. Background Technology

[0002] Since the 1930s, Chinese soil researchers have conducted extensive studies on the formation and physicochemical properties of paddy soils. Gleyed paddy soils, in the Chinese soil classification system, belong to the Anthropoid Soils, Subclass Anthropoid Hydrophytes, and subclass Paddy Soils. They are reductive paddy soils developed under long-term flooding conditions. Due to prolonged flooding, soil aggregate structure is destroyed, soil permeability deteriorates, and the root environment experiences multi-layered anoxic conditions, leading to the accumulation of organic acids, Fe2+, Mn2+, and H2S. Gleyed paddy soils are among the low-yielding paddy soils in my country, but also among those with significant yield-increasing potential. Existing research has revealed the causes, hazards, and improvement pathways of gleyed paddy soils. However, quickly identifying the location and distribution of the gleyed layer in paddy soils is a prerequisite and foundation for their improvement.

[0003] Soil barrier layer thickness has a significant impact on crops. Traditional methods for obtaining soil barrier layer thickness include soil profiling, drilling, probing, and permeability measurement. These methods have advantages such as accuracy, intuitiveness, and comprehensive data sets for detecting gley layer thickness. However, they often require a lot of manpower and resources, damage the soil structure, and are difficult, costly, and inefficient to sample. They are suitable for small areas but difficult to apply to large-scale continuous detection and are often used as verification methods. Summary of the Invention

[0004] This application provides a method for determining the continuous spatial variation of gleyed layers based on variational mode decomposition and ground penetrating radar technology. This method solves the technical problems of existing ground penetrating radar systems, such as the inability to identify underground information by relying solely on grayscale images, the significant influence of subjective identification, and the inability of soil obstacle identification methods to perform rapid, continuous, and large-scale measurements.

[0005] This application provides a method for determining the continuous spatial variation of the gley layer based on variational mode decomposition and ground-penetrating radar technology, comprising the following steps:

[0006] S1. Select the test area and lay out a grid-shaped radar survey line with two horizontal and two vertical lines. Extend the survey line from the center to both sides and drill a set of soil drills at preset intervals as verification soil drills. Drill a soil drill in the center of the test area as a reference soil drill. Add o-phenanthroline reagent to determine the location of the gleyed layer in typical gleyed paddy soil.

[0007] S2. Ground penetrating radar is used to detect and collect data on the soil gley layer. At least two longitudinal and two transverse survey lines are laid out in the identified area to obtain the ground penetrating radar data signals of the soil borer profile and the original signals of the ground penetrating radar data of the survey lines on the four survey lines. The original signals include the time window of the signal propagation in the soil and the amplitude of different time windows under two-way travel time. During the ground penetrating radar detection process, a signal is collected at a set distance every time interval. The time window-amplitude data of the total number of signals collected on the ground penetrating radar route are combined to form the original radar signals of the soil borer profile and the survey line profile.

[0008] S3. Preprocess each acquired raw radar signal to determine the ground wave position and perform filtering to eliminate sampling errors and signal noise;

[0009] S4. Perform variational mode decomposition on the preprocessed original signal, set a fixed number of decomposition layers, process each signal one by one to obtain the optimal decomposition model, retain the original signal features to the maximum extent, analyze the center frequency of the decomposed signal, analyze the correlation between the modal components after variational mode decomposition and the original signal through the correlation matrix, select the modal components with high correlation with the original signal as the optimal modal component combination, and obtain the radar signal for interpreting soil structure based on the optimal modal component combination and the component center frequency.

[0010] S5. Based on the radar signal obtained to interpret the soil structure, a grayscale map is drawn. Combining the location of the soil gleyed layer and the grayscale map determined in step S1, the characteristic wave corresponding to the location of the gleyed layer is obtained. The peak position of the radar signal after each variational mode decomposition is automatically filtered. Based on the filtered characteristic wave, the corresponding positions of the upper and lower boundaries of the gleyed layer are identified on each measurement line, and the thickness of the gleyed layer is calculated.

[0011] S6. Based on the upper and lower boundaries of the gecko layer obtained from each survey line, assign spatial attributes to the survey line data, set three-dimensional coordinate axes, (X, Y) is the location of the survey line, Z1 and Z2 represent the upper and lower boundaries of the gecko layer, and use inverse range weight interpolation to obtain the upper and lower boundaries of the gecko layer in the area measured by the radar.

[0012] S7. Based on the actual value of the gley layer thickness obtained from the soil drill in S1, the gley layer thickness at the corresponding location of the radar survey line is used as the predicted value. The correlation coefficient R of the linear regression model is then used to determine the gley layer thickness. 2 To evaluate the accuracy of the predicted gill layer thickness.

[0013] Preferably, step S1 includes:

[0014] S1.1 Drill a soil drill in the center of the gleyed paddy soil test area as a reference soil drill. The soil drill should be 1m deep. Scrape off the surface soil of the soil column and add o-phenanthroline reagent evenly and continuously from top to bottom. The area that turns red is the gleyed layer. Take pictures of the soil drill and record the location of the gleyed layer. Drill two more verification soil drills on both sides 20cm away from the reference soil drill to determine the location of the gleyed layer. If the location is not much different from the reference soil drill, the soil drill is reliable. Otherwise, find other locations and drill soil drills again.

[0015] S1.2. Four survey lines, two horizontal and two vertical, are laid out in the test area in a grid pattern. The surface of the survey lines is flattened using a roller. A soil drill is driven into the center of each survey line. Using this soil drill as the origin, another soil drill is driven outward every 5 meters. O-phenanthroline reagent is dripped onto each soil drill. The location where the soil turns red is the location of the gleyed layer. The soil drills are numbered, photographed, and their gleyed layer locations are recorded as verification points.

[0016] Preferably, step S2 includes:

[0017] S2.1. Lay out two horizontal and two vertical radar survey lines in the test area, which overlap with the survey lines in step S1.2, in a well-shaped pattern, running through the entire test area, and then level the ground along the survey lines.

[0018] S2.2 Setting sampling parameters, including 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 fully acquired to facilitate subsequent ground location, and a relatively large value is used when setting the gain;

[0019] S2.3. Sample 1 m to 2 m along the top of the center position of the soil drill, completely including the reference soil drill and the verification soil drill, repeat the sampling 3 times to obtain the original soil drill ground-penetrating radar data signal;

[0020] S2.4. Sample along the established radar survey line, repeating the sampling twice. During the repeated sampling, ensure that the same starting point and the same ending point are used, and that the radar vehicle moves at a constant speed to obtain the original ground-penetrating radar data signal of the survey line.

[0021] Preferably, the preprocessing method in step S3 includes the following steps:

[0022] S3.1 Construct a signal filter. Perform filtering preprocessing on the sampled raw ground-penetrating radar data signal. According to frequency analysis, the noise is in the high-frequency part of the frequency. Use a high-frequency filter to filter out the noise part.

[0023] S3.2 Calculate the location of the gleyed layer starting from the ground, construct a method to find the ground at time 0, the ground location is the lowest value of the radar signal, set an initial threshold of -8000, use the threshold increment method to find the lowest value, and use the lowest value of the radar signal as the ground location;

[0024] S3.3. Add zero values ​​of the same length to each filtered ground-penetrating radar data to make the length of each radar signal data the same.

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

[0026] S4.1 Construct an adaptive method for obtaining the center frequency of the decomposition layer. Iterate each channel of the original signal separately, fix the number of iterations at 6, that is, fix the number of decomposition layers at 6, and obtain the center frequency of 6 iterations. Based on the obtained fixed number of decomposition layers and center frequency, perform variational mode decomposition on each signal to obtain 6 modal components of each signal on the soil profile.

[0027] S4.2. The frequency distribution of the original signal is calculated using Fast Fourier Transform, and the spectra of all channels in the profile are superimposed to obtain the initial frequency distribution of the signal. Based on the data of each channel, the correlation matrix between the original signal and the 6 modal components is constructed using the Pearson correlation coefficient. Combining the center frequency and the correlation matrix, the modal components with high correlation to the original signal and low frequency are selected and superimposed to obtain the new radar signal. The above steps are repeated for all data on the soil drill profile and the survey line to obtain the optimal modal combination dataset of the soil drill profile and the optimal modal combination dataset of the survey line.

[0028] S4.3. Based on the modal components obtained in step S4.2, calculate the Pearson correlation coefficients between each component and between the components and the original signal.

[0029] Preferably, step S4.1 includes:

[0030] S4.1.1, Assuming the multi-component signal is composed of... Composed of modal components with finite bandwidth With a fixed K value of 6, the center frequency of each intrinsic mode function (IMF) is... The constraint is that the modal sum equals the input signal, which is obtained through Hilbert transform. The analysis signal is obtained, and its one-sided spectrum is calculated by comparing it with the operator. Multiply, The center band is modulated to the corresponding baseband:

[0031]

[0032] in It is the Dirac impulse function; Represents convolution; The imaginary unit; Angular frequency; For natural index;

[0033] S4.1.2 Calculate the square norm of the demodulation gradient. And estimate the bandwidth of each modulus component, as shown in the following formula:

[0034]

[0035] In the above formula, , This represents the decomposed IMF components. Represents the center frequency of each component; * represents the Dirac function; * represents the convolution operator. The original signal;

[0036] S4.1.3, Introducing Lagrange multipliers and second-order penalty factor The constrained variational problem is transformed into an unconstrained variational problem, and the extended Lagrange expression is as follows:

[0037]

[0038] S4.1.4. By continuously updating each component and its center frequency using the alternating direction multiplier method, the saddle point of the unconstrained model is finally obtained, which is the optimal solution to the original problem. All components can be obtained from the frequency domain space using the following formula:

[0039] ;

[0040] S4.1.5 yes After Wiener filtering, the remaining quantity is used by the algorithm to re-estimate the barycenter frequency based on the barycenter of the power spectrum of each component. The specific process is as follows:

[0041] a) , , and ;

[0042] b) Period: ;

[0043] c) Updated according to S4.1.4. ;

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

[0045]

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

[0047]

[0048] In the formula: To meet the noise tolerance and signal fidelity requirements for decomposition; and Corresponding to and Fourier transform;

[0049] f) Repeat steps a~f until the iteration stopping condition is met. The stopping condition is:

[0050]

[0051] Finally obtained after stopping iteration IMF components and center frequency.

[0052] Preferably, step 4.3 includes the following steps:

[0053] S4.3.1 Based on the six modal components obtained in step S4.2, a correlation matrix {GPR IMF1 IMF2 IMF3 IMF4 IMF5 IMF6} is formed with the original GPR signal. By calculating the ratio of the covariance to the standard deviation between the variables, the influence of dimensions is eliminated, and the standardized linear correlation index is obtained. The formula principle is as follows:

[0054]

[0055] in, The closer the correlation is to ±1, the stronger the linear correlation; the closer it is to 0, the weaker the linear correlation. and is a column vector in the matrix.

[0056] S4.3.2. Based on the correlation between the original signal and the modal components in the correlation matrix, a correlation of less than 0.4 is generally considered to be weak or non-correlated. The strongly correlated modal components are retained as the optimal modal component combination to obtain the decomposed spectrum on the profile. The radar spectrum on the survey line is processed based on steps S4.1, S4.2 and S4.3.1.

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

[0058] S5.1. Based on the radar signal on the soil drill profile obtained in step S4.2, the RMS value is obtained using a sliding window. The gain is normalized by division, and a grayscale image is plotted on the soil drill profile. The amplitude change of the radar signal will cause the color change of the grayscale image. Based on the location of the gleyed layer obtained by the soil drill in step S1.1, the representative part of the gleyed layer in the grayscale image is found by the soil drill. The signal characteristic peak positions representing the upper and lower boundaries of the gleyed layer are screened by combining the grayscale image and the soil drill data to obtain the location signal of the gleyed layer.

[0059] S5.2 Based on the original ground-penetrating radar data of the survey line obtained in step S2.4, obtain the horizontal resolution of the antenna, that is, the number of radar waves when the measurement distance is 1 meter, set the initial radar wave position, randomly select a radar wave every 10 centimeters of the survey line, and repeat the processing for all survey lines.

[0060] S5.3. Based on the radar waves on the survey line in step S5.2, apply a smooth envelope to the single-channel radar wave. Based on the instantaneous phase change position of the radar signal envelope and the signal characteristic peak positions of the upper and lower boundaries of the gryllium obtained in step S5.1, obtain the upper and lower positions of the gryllium, calculate the difference, and obtain the gryllium thickness position. Repeat the process for all survey lines.

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

[0062] S6.1. Taking one corner of the test area as the origin, the north-south direction as the Y-axis and the east-west direction as the X-axis, assign the plane position coordinates of the survey line, that is, the coordinates (X, Y) of the survey line position relative to the origin. Combined with the upper and lower boundaries Z1 and Z2 of the gleyed layer obtained on the survey line in step S5.3, obtain the spatial position of the upper and lower boundaries of the gleyed layer on the survey line in the three-dimensional coordinate system.

[0063] S6.2. Based on the spatial positions of the upper and lower boundaries of the gleyed layer on the survey line obtained in step S6.1, perform inverse distance weight interpolation on the upper and lower boundaries of the gleyed layer on the survey line. The distance weight is 3 meters. Start iterating from 0 and stop iterating when the error is less than the preset threshold to obtain a three-dimensional spatial dataset of the gleyed layer position on the plane.

[0064] S6.3. Based on the three-dimensional spatial dataset of the gryllium location obtained in step S6.2, draw a three-dimensional spatial coordinate system. The location of the gryllium is between the upper and lower boundaries of the gryllium.

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

[0066] S7.1 Based on the radar survey line laid out in step S1.2, starting from the center of the survey line and extending to both sides, drill a soil drill every 5 meters and number it, record its position on the survey line and the thickness of the gley layer as the observation value, and based on the thickness of the gley layer on the radar survey line obtained in step S5.3, pick the thickness of the gley layer at the corresponding soil drill position as the prediction value.

[0067] S7.2. Based on the linear regression model, perform linear fitting between the observed and predicted values, using R0... 2 To evaluate model accuracy, R 2 The higher the precision, the higher the accuracy.

[0068] Beneficial technical effects of the present invention:

[0069] 1. The technical solution of this invention uses ground-penetrating radar to directly measure large-scale soil, and automatically filters out ground-penetrating radar signal clutter and corrects the ground position. It can detect the location and thickness of the gley layer non-destructively, continuously, and rapidly, and efficiently obtain the location and thickness of the gley layer. It solves the shortcomings of existing ground-penetrating radar systems that only identify underground information through grayscale images, which are greatly affected by the subjectivity of identification, and traditional gley layer identification methods that cannot perform rapid, continuous, and large-scale measurements.

[0070] 2. This invention utilizes ground-penetrating radar (GPR) technology to directly measure the location of the gleyed layer in large-scale gleyed paddy soil, automatically filtering out GPR signal clutter and correcting the ground position. The technical solution of this invention is based on variational mode decomposition (VMD), setting K to 6. Low-frequency mode components are obtained through the center frequency. A Pearson correlation matrix is ​​constructed for the original signal and each mode component. The optimal combination of mode components is extracted. The grayscale image of the optimal combination of mode components and the location of the gleyed layer using a soil drill are combined to extract the characteristic wave information of the gleyed layer. The location of the gleyed layer on the radar survey line is then extracted. Using the location of the gleyed layer on the radar survey line, the location of the gleyed layer in the test area is obtained through inverse range weighted interpolation. Furthermore, the linear fitting R... 2 This technique demonstrates the accuracy of obtaining gill thickness along the survey line. It enables rapid acquisition of gill thickness on a large scale, improving identification efficiency while ensuring accuracy. Attached Figure Description

[0071] Figure 1 This is a schematic diagram illustrating the model construction and actual measurement process during the specific implementation of this invention;

[0072] Figure 2 This is a schematic diagram of the radar survey line layout and soil drill position in an embodiment of the present invention;

[0073] Figure 3 This is a schematic diagram illustrating the correlation between the modal components after variational mode decomposition and the original signal in an embodiment of the present invention.

[0074] Figure 4 This is a schematic diagram of ground zero-point correction in an embodiment of the present invention and an effect diagram of each modal component after variational mode decomposition of a single channel of ground penetrating radar data.

[0075] Figure 5 This is a schematic diagram illustrating the search for characteristic waves by comparing grayscale images with soil drills in an embodiment of the present invention.

[0076] Figure 6 This is a rendering showing the position of the gleyed layer in two gleyed paddy fields in a three-dimensional coordinate system according to an embodiment of the present invention.

[0077] Figure 7 This is a graph showing the linear fitting effect of the gleyed layer thickness in the experimental field according to an embodiment of the present invention. Detailed Implementation

[0078] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0079] like Figure 1 As shown in the figure, this embodiment provides a method for determining the continuous spatial variation of the gley layer based on variational mode decomposition and ground penetrating radar technology, which includes the following steps:

[0080] S1. Select the test area, lay out radar survey lines, drill mapping drills at the center of the test area and along the survey lines, and add o-phenanthroline reagent to determine the location of the gleyed layer in typical gleyed paddy soils, such as... Figure 2 As shown. Specific operations include:

[0081] S1.1. Drill a soil drill at the center of the gleyed paddy soil study area as a reference drill. The drill should be 1m deep. Scrape away the surface soil of the column and continuously and evenly drip o-phenanthroline reagent from top to bottom. The area that turns red is the gleyed layer. Take pictures of the drill and record the location of the gleyed layer. Drill two more verification drills 20cm away from the original drill on both sides to determine the location of the gleyed layer. If the location is not significantly different from the original drill, the original drill is reliable; otherwise, find other locations and drill new drills.

[0082] S1.2. Four survey lines, two horizontal and two vertical, are laid out in the test area. The surface of the survey lines is flattened with a roller. A soil drill is driven into the center of each survey line. Using this soil drill as the origin, a soil drill is driven outward every 5 meters. O-phenanthroline reagent is dripped onto each soil drill. The location where the soil turns red is the location of the gleyed layer. The soil drills are numbered, photographed, and their gleyed layer locations are recorded as verification points.

[0083] S2. Use ground-penetrating radar to detect and sample the gleyed paddy soil to obtain raw signals; specifically including the following steps:

[0084] S2.1 Before conducting ground-penetrating radar detection, the ground should be leveled to reduce the impact of uneven ground on the calculation of black soil layer thickness and on ground-coupled antenna detection.

[0085] S2.1 First, the ground is leveled to reduce the impact of uneven ground on the calculation of the gryllium location and the detection of the ground-coupled antenna;

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

[0087] ;

[0088] D represents the expected detection depth in meters (m). V represents the average radar wave velocity in the formation medium in m / ns. W represents the sampling time window in ns. The larger the time window, the greater the detection depth, but the lower the resolution.

[0089] Depth is typically calculated by converting time windows into depth using the dielectric constant, estimating the depth of each boundary line and the thickness of the generation layer. The dielectric constant is the ability of a material to retain electrical charge; its magnitude determines the medium's ability to absorb or reflect electromagnetic waves, ranging from 1 (air) to 81 (water). Its calculation formula is shown below:

[0090]

[0091] Transform the above equation into: ;

[0092] The formula for calculating the depth (H, m) for calculating the two-way travel time of electromagnetic waves is as follows:

[0093]

[0094] In the formula, ε is the dielectric constant, and c is the speed of light (3.00 × 10⁻⁶). 8 m / s); v (m / s) is the speed at which electromagnetic waves propagate in the medium.

[0095] Since the ground-coupled antenna is close to the ground, the offset setting is to ensure that the air wave, i.e., the first positive wave peak, is completely collected to facilitate subsequent location of the ground. Because the paddy soil in the gley layer significantly reduces the amplitude, a larger value is used when setting the gain.

[0096] S2.3. Based on the sampling parameters set in S2.2, sample 1-2 m along the top of the soil drill profile to completely include the entire profile. Repeat the sampling 3 times to obtain the original ground-penetrating radar data signal of the soil drill profile.

[0097] S2.4. Sample along the established radar survey line, repeating the sampling twice. During repeated sampling, ensure the same starting point and ending point, and maintain a constant speed for the radar vehicle to obtain the raw ground-penetrating radar data signal for the survey line. The raw ground-penetrating radar data signal includes the time window (time window) of signal propagation in the soil and the signal intensity (amplitude) of different time windows during the two-way travel. During the ground-penetrating radar detection process, the machine collects a signal at regular intervals (in this invention, the interval between each signal is approximately 2 cm). The time window-amplitude data of the total number of signals collected along the ground-penetrating radar route are combined to form the raw radar signal of the profile.

[0098] S3. Preprocess the raw radar signal data by filtering and determining the ground position at time 0 to eliminate signal noise and sampling errors; specifically including the following steps:

[0099] S3.1 Constructing a signal filter: Since the sampling signal of the ground-penetrating radar in the black soil is severely attenuated, in order to improve the readability of the grayscale image during sampling, a large gain parameter needs to be set, which generates a large amount of machine noise. Therefore, the sampled signal needs to be preprocessed by filtering. According to frequency analysis, the noise is in the high-frequency part of the frequency. A high-frequency filter is used to filter out the noise part.

[0100] S3.2 Calculate the location of the gleyed layer starting from the ground, construct a method to find the ground at time 0, the ground location is the lowest value of the radar signal, set an initial threshold of -8000, use the threshold increment method to find the lowest value, and use the lowest value of the radar signal as the ground location;

[0101] S3.3. Each eliminated ground-penetrating radar data is then replaced with a zero value of the same length to make the length of each radar signal data the same.

[0102] S4. Perform variational mode decomposition on the preprocessed original radar signal, fixing the number of decomposition layers at 6, to obtain mode components that retain the characteristics of the original signal. Use the Pearson correlation matrix to retain components with high correlation to the original signal (correlation greater than 0.4 is retained), thus obtaining the optimal decomposed mode model and preserving the characteristics of the original signal to the maximum extent. Specific steps include:

[0103] S4.1. Construct a center frequency method with a fixed number of decomposition levels, iterating each channel of the original signal separately. The specific process of variational mode decomposition (VMD) can be understood as finding the optimal solution to a variational problem, and can be correspondingly transformed into constructing and solving a variational problem. The specific construction steps are as follows:

[0104] S4.1.1, Assuming the multi-component signal is composed of... Composed of modal components with finite bandwidth The center frequency of each intrinsic mode function (IMF) is The constraint is that the modal sum equals the input signal, which is obtained through Hilbert transform. The analysis signal is obtained, and its one-sided spectrum is calculated by comparing it with the operator. Multiply, The center band is modulated to the corresponding baseband:

[0105]

[0106] in It is the Dirac impulse function; Represents convolution; The imaginary unit; Angular frequency; For natural index;

[0107] S4.1.2 Calculate the square norm of the demodulation gradient. And estimate the bandwidth of each modulus component, as shown in the following equation:

[0108]

[0109] In the above formula, , This represents the decomposed IMF components. Represents the center frequency of each component; * represents the Dirac function; * represents the convolution operator. The original signal;

[0110] S4.1.3 To find the optimal solution to the constrained variational problem, we first introduce Lagrange multipliers. and second-order penalty factor This transforms the constrained variational problem into an unconstrained variational problem. The second-order penalty factor is involved. This ensures the accuracy of signal reconstruction in Gaussian noise environments. Lagrange multipliers. This ensures that the constraints remain strict. The extended Lagrange expression is as follows:

[0111]

[0112] S4.1.4. Using the Alternating Direction Multiplier Method (ADMM), each component and its center frequency are continuously updated to finally obtain the saddle point of the unconstrained model, which is the optimal solution to the original problem. All components can be obtained from the frequency domain space using the following formula:

[0113]

[0114] S4.1.5 yes After Wiener filtering, the remaining quantity is used by the algorithm to re-estimate the barycenter frequency based on the barycenter of the power spectrum of each component. The specific process is as follows:

[0115] a) , , and ;

[0116] b) Period: ;

[0117] c) Updated according to S4.1.4. ;

[0118] d) Update according to the formula ;

[0119]

[0120] e) Update according to formula ;

[0121]

[0122] In the formula: To meet the noise tolerance and signal fidelity requirements for decomposition; and Corresponding to and Fourier transform.

[0123] f) Repeat steps a~f until the iteration stopping condition is met, as shown in the equation;

[0124]

[0125] Finally obtained after stopping iteration The IMF components and center frequency are determined. Typically, when the number of decomposition layers is greater than 6, the correlation with the original signal is extremely low; therefore, the first six low-frequency mode components (K=6) are selected as the research object.

[0126] S4.2. Based on the decomposition layer number obtained in S4.1, variational mode decomposition is performed on each channel of the signal to obtain different decomposed modes for each channel on the profile. The signal data acquired by ground penetrating radar is three-dimensional data of channel number-time window-amplitude. Assume the signal data has a total of... Each channel represents a two-dimensional data set of time window and amplitude. Variational mode decomposition is performed, and the center frequency method is used to extract the six low-frequency components of each channel, such as... Figure 4 As shown.

[0127] S4.3. Based on the modal components obtained in S4.2, calculate the Pearson correlation coefficients between each component and between each component and the original signal. For example... Figure 3 As shown, this illustrates the correlation between modal components and the original spectrum. The closer the absolute value of the correlation coefficient is to 1, the stronger the correlation. A correlation less than 0.4 is generally considered weak or nonexistent. Strongly correlated modal components are retained as the optimal modal component combination. The optimal modal combinations for each trace are then superimposed to obtain new signal data in three dimensions: trace number, time window, and amplitude. The above steps are repeated for all data from the soil drill profile and survey line to obtain the optimal modal combination datasets for the soil drill profile and survey line. The main function of the Pearson correlation coefficient is to quantify the degree and direction of linear correlation between two continuous variables. Its value range is [-1, 1], and it is used to determine the positive / negative correlation and its strength between variables. The specific steps are as follows:

[0128] S4.3.1 Based on the six modal components obtained in step S4.2, a correlation matrix {GPR IMF1 IMF2 IMF3 IMF4 IMF5 IMF6} is formed with the original GPR signal. By calculating the ratio of the covariance to the standard deviation between the variables, the influence of dimensions is eliminated, and the standardized linear correlation index is obtained. The formula principle is as follows:

[0129]

[0130] in, The closer the correlation is to ±1, the stronger the linear correlation; closer to 0, there is no linear correlation. and is the column vector in the matrix. As shown in Table 1, the correlation matrix between the original signal and mode components of a 1GHz antenna single channel is used. {IMF1IMF2 IMF3 IMF4} in the matrix is ​​selected as the optimal mode combination.

[0131]

[0132] S4.3.2. Based on the correlation between the original signal and the modal components in the correlation matrix, a correlation of less than 0.4 is generally considered to be weak or non-correlated. The strongly correlated modal components are retained as the optimal modal component combination to obtain the decomposed spectrum on the profile. The radar spectrum on the survey line is processed based on steps S4.1, S4.2 and S4.3.1.

[0133] S5. Combining the grayscale image, soil drill profile, and optimal mode combination, extract the characteristic waves and characteristic positions of the gleyed layer location from the soil drill profile. Then, using the radar signal from the survey line and the extracted characteristic waves and characteristic positions, obtain the gleyed layer location on the survey line. Specifically, this includes the following steps:

[0134] S5.1. Based on the radar signal on the soil drill profile obtained in step S4.2, the RMS value is obtained using a sliding window. Gain is normalized by division, and a grayscale image of the soil drill profile is plotted. Changes in the radar signal amplitude will cause color changes in the grayscale image. Combining the gleyed layer location obtained from the reference soil drill in S1.1, the representative gleyed layer portion in the grayscale image is found using the soil drill. The grayscale image and soil drill data are used to filter the signal characteristic peak positions representing the upper and lower boundaries of the gleyed layer, i.e., the characteristic peak or trough position in the signal spectrum (e.g., the fifth peak). Simultaneously, a smooth envelope is used to reflect the instantaneous changes in the spectrum, such as... Figure 5 As shown in the figure, the spectral curve is composed of the optimal mode combination, and the grayscale image is drawn from the optimal mode combination dataset. The changes in the grayscale image reflect the changes in soil structure. The characteristic peaks on the spectral curve are found by comparing the grayscale image with the soil drill.

[0135] S5.2 Based on the radar data of the survey line obtained in S4.3, obtain the horizontal resolution of the antenna, that is, the number of radar waves when the measurement distance is 1 meter, set the initial radar wave position, and randomly select a radar wave every 10 centimeters of the survey line (in this embodiment, the radar is about 2 centimeters apart, that is, one wave is taken every 5 waves). Repeat the processing on all survey lines to obtain the optimal mode combination dataset of the survey line with a horizontal resolution of 10 centimeters.

[0136] S5.3, based on the optimal mode combination dataset of the survey line in S5.2, applies a smooth envelope to the single-channel radar wave, filters the original signal based on the characteristic peak positions of the upper and lower boundaries of the gleyed layer extracted from the soil drill profile, extracts the peak positions of the survey line radar data (e.g., the fifth peak) as the gleyed layer boundary positions, automatically identifies the corresponding positions of the upper and lower boundaries of the gleyed layer, calculates the difference, obtains the gleyed layer thickness position, and repeats the processing for all survey lines.

[0137] S6. By extracting the location of the gleyed layer along the radar survey line, and using inverse range weighted interpolation, the location of the gleyed layer in the experimental field is obtained. This specifically includes the following steps:

[0138] S6.1. Taking one corner of the test area as the origin, the north-south direction as the Y-axis and the east-west direction as the X-axis, assign the plane position coordinates of the survey line, that is, the coordinates (X, Y) of the survey line position relative to the origin. Combined with the upper and lower boundaries Z1 and Z2 of the gleyed layer obtained on the survey line in S5.3, obtain the spatial position of the upper and lower boundaries of the gleyed layer on the survey line in the three-dimensional coordinate system.

[0139] S6.2. Based on the spatial location of the upper and lower boundaries of the gley layer on the survey line obtained in step S6.1, Matlab software is used to perform inverse distance weight interpolation on the upper and lower boundaries of the gley layer on the survey line. The distance weight is set to 3 meters, and the iteration starts from 1. The iteration stops when the error is less than 1 cm, and a three-dimensional spatial grid dataset of the location of the gley layer on the field of the experimental area is obtained with a spatial resolution of 0.1 meters.

[0140] S6.3. Based on the three-dimensional spatial dataset of the gryllus location obtained in step S6.2, draw a three-dimensional spatial coordinate system. The location of the gryllus is between the upper and lower boundaries of the gryllus. Figure 6 As shown.

[0141] S7. Based on the linear fit between the gleyed layer thickness obtained on the survey line in step S6.1 and the gleyed layer thickness obtained by the soil drill in step S1.2, and using R... 2 To describe the model accuracy of the method for determining the spatial continuity of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology, the following steps are specifically included:

[0142] S7.1 Based on the radar survey line laid out in step S1.2, starting from the center of the survey line and extending to both sides, drill a soil drill every 5 meters and number it, recording its position on the survey line and the thickness of the gley layer as observed values. Based on the gley layer thickness on the radar survey line obtained in step S5.3, pick the gley layer thickness at the corresponding soil drill position as the predicted value;

[0143] S7.2. Based on the linear regression model, perform linear fitting between the observed and predicted values, using R0... 2 To evaluate model accuracy, R 2 Higher precision means higher accuracy. For example, in this embodiment, a 1 GHz antenna can detect the R value of the gryllium thickness. 2 Up to 0.9, such as Figure 7 As shown, the data comes from three experimental areas of gleyed paddy soil. Figure 7 (a) The figure shows the linear fit of the observed and predicted values ​​of a radar antenna at a frequency of 1 GHz using a linear regression model, with a fitting accuracy of 0.83. Figure 7 (b) The figure shows the linear fitting of the observed and predicted values ​​by a linear regression model for a radar antenna with a frequency of 700 MHz, with a fitting accuracy of 0.9.

[0144] Therefore, the method constructed based on variational mode decomposition and ground-penetrating radar technology to determine the spatial continuous change of the gleyed layer can efficiently and accurately identify the location of the gleyed layer, providing data support for reducing the gleyed paddy soil barrier layer and increasing grain yield.

[0145] In summary, this invention, by collecting ground-penetrating radar data and other data, employing various data processing methods, and integrating radar data and field data characteristics, enables large-scale detection of the gleyed layer spatial location in gleyed paddy soil. This improves the efficiency and accuracy of obtaining gleyed layer thickness and solves the problem that existing gleyed layer thickness detection methods struggle to achieve large-scale, rapid, and non-destructive detection.

[0146] The above description is only a preferred embodiment of the present invention. It should be noted that for those skilled in the art, several improvements can be made without departing from the principle of the present invention, and these improvements should also be considered within the scope of protection of the present invention.

Claims

1. A method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology, characterized in that, Includes the following steps: S1. Select the test area and lay out a grid-shaped radar survey line with two horizontal and two vertical lines. Extend the survey line from the center to both sides and drill a set of soil drills at preset intervals as verification soil drills. Drill a soil drill in the center of the test area as a reference soil drill. Add o-phenanthroline reagent to determine the location of the gleyed layer in typical gleyed paddy soil. S2. Ground penetrating radar is used to detect and collect data on the soil gley layer. At least two longitudinal and two transverse survey lines are laid out in the identified area to obtain the ground penetrating radar data signals of the soil borer profile and the original signals of the ground penetrating radar data of the survey lines on the four survey lines. The original signals include the time window of the signal propagation in the soil and the amplitude of different time windows under two-way travel time. During the ground penetrating radar detection process, a signal is collected at a set distance every time interval. The time window-amplitude data of the total number of signals collected on the ground penetrating radar route are combined to form the original radar signals of the soil borer profile and the survey line profile. S3. Preprocess each acquired raw radar signal to determine the ground wave position and perform filtering to eliminate sampling errors and signal noise; S4. Perform variational mode decomposition on the preprocessed original signal, set a fixed number of decomposition layers, process each signal one by one to obtain the optimal decomposition model, retain the original signal features to the maximum extent, analyze the center frequency of the decomposed signal, analyze the correlation between the modal components after variational mode decomposition and the original signal through the correlation matrix, select the modal components with high correlation with the original signal as the optimal modal component combination, and obtain the radar signal for interpreting soil structure based on the optimal modal component combination and the component center frequency. S5. Based on the radar signal obtained to interpret the soil structure, a grayscale map is drawn. Combining the location of the soil gleyed layer and the grayscale map determined in step S1, the characteristic wave corresponding to the location of the gleyed layer is obtained. The peak position of the radar signal after each variational mode decomposition is automatically filtered. Based on the filtered characteristic wave, the corresponding positions of the upper and lower boundaries of the gleyed layer are identified on each measurement line, and the thickness of the gleyed layer is calculated. S6. Based on the upper and lower boundaries of the gecko layer obtained from each survey line, assign spatial attributes to the survey line data, set three-dimensional coordinate axes, (X, Y) is the location of the survey line, Z1 and Z2 represent the upper and lower boundaries of the gecko layer, and use inverse range weight interpolation to obtain the upper and lower boundaries of the gecko layer in the area measured by the radar. S7. Based on the actual value of the gley layer thickness obtained from the soil drill in S1, the gley layer thickness at the corresponding location of the radar survey line is used as the predicted value. The correlation coefficient R of the linear regression model is then used to determine the gley layer thickness. 2 Evaluate the accuracy of gill layer thickness prediction; Step S4 includes the following steps: S4.1 Construct an adaptive method for obtaining the center frequency of the decomposition layer. Iterate each channel of the original signal separately, fix the number of iterations at 6, that is, fix the number of decomposition layers at 6, and obtain the center frequency of 6 iterations. Based on the obtained fixed number of decomposition layers and center frequency, perform variational mode decomposition on each signal to obtain 6 modal components of each signal on the soil profile. S4.

2. The frequency distribution of the original signal is calculated using Fast Fourier Transform, and the spectra of all channels in the profile are superimposed to obtain the initial frequency distribution of the signal. Based on the data of each channel, the correlation matrix between the original signal and the 6 modal components is constructed using the Pearson correlation coefficient. Combining the center frequency and the correlation matrix, the modal components with high correlation to the original signal and low frequency are selected and superimposed to obtain the new radar signal. The above steps are repeated for all data on the soil drill profile and the survey line to obtain the optimal modal combination dataset of the soil drill profile and the optimal modal combination dataset of the survey line. S4.

3. Based on the modal components obtained in step S4.2, calculate the Pearson correlation coefficients between each component and between the components and the original signal.

2. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 1, characterized in that, Step S1 includes: S1.1 Drill a soil drill in the center of the gleyed paddy soil test area as a reference soil drill. The soil drill should be 1m deep. Scrape off the surface soil of the soil column and add o-phenanthroline reagent evenly and continuously from top to bottom. The area that turns red is the gleyed layer. Take pictures of the soil drill and record the location of the gleyed layer. Drill two more verification soil drills on both sides 20cm away from the reference soil drill to determine the location of the gleyed layer. If the location is not much different from the reference soil drill, the soil drill is reliable. Otherwise, find other locations and drill soil drills again. S1.

2. Four survey lines, two horizontal and two vertical, are laid out in the test area in a grid pattern. The surface of the survey lines is flattened using a roller. A soil drill is driven into the center of each survey line. Using this soil drill as the origin, another soil drill is driven outward every 5 meters. O-phenanthroline reagent is dripped onto each soil drill. The location where the soil turns red is the location of the gleyed layer. The soil drills are numbered, photographed, and their gleyed layer locations are recorded as verification points.

3. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 2, characterized in that, Step S2 includes: S2.

1. Lay out two horizontal and two vertical radar survey lines in a well-shaped pattern in the test area, running through the entire test area, and then level the ground along the survey lines. S2.2 Setting sampling parameters, including 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 fully acquired to facilitate subsequent ground location, and a relatively large value is used when setting the gain; S2.

3. Sample 1 m to 2 m along the top of the center position of the soil drill, completely including the reference soil drill and the verification soil drill, repeat the sampling 3 times to obtain the original soil drill ground-penetrating radar data signal; S2.

4. Sample along the established radar survey line, repeating the sampling twice. During the repeated sampling, ensure that the same starting point and the same ending point are used, and that the radar vehicle moves at a constant speed to obtain the original ground-penetrating radar data signal of the survey line.

4. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 1, characterized in that, The preprocessing method in step S3 includes the following steps: S3.1 Construct a signal filter. Perform filtering preprocessing on the sampled raw ground-penetrating radar data signal. According to frequency analysis, the noise is in the high-frequency part of the frequency. Use a high-frequency filter to filter out the noise part. S3.2 Calculate the location of the gleyed layer starting from the ground, construct a method to find the ground at time 0, the ground location is the lowest value of the radar signal, set an initial threshold of -8000, use the threshold increment method to find the lowest value, and use the lowest value of the radar signal as the ground location; S3.

3. Add zero values ​​of the same length to each filtered ground-penetrating radar data to make the length of each radar signal data the same.

5. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 4, characterized in that, Step S4.1 includes: S4.1.1, Assuming the multi-component signal is composed of... Composed of modal components with finite bandwidth With a fixed K value of 6, the center frequency of each intrinsic mode function (IMF) is... The constraint is that the modal sum equals the input signal, which is obtained through Hilbert transform. The analysis signal is obtained, and its one-sided spectrum is calculated by comparing it with the operator. Multiply, The center band is modulated to the corresponding baseband: ； in It is the Dirac impulse function; * denotes convolution; The imaginary unit; Angular frequency; For natural index; S4.1.2 Calculate the square norm of the demodulation gradient. And estimate the bandwidth of each modulus component, as shown in the following formula: ； In the above formula, , This represents the decomposed IMF components. Represents the center frequency of each component; * represents the Dirac function; * represents the convolution operator. The original signal; S4.1.3, Introducing Lagrange multipliers and second-order penalty factor The constrained variational problem is transformed into an unconstrained variational problem, and the extended Lagrange expression is as follows: ； S4.1.

4. By continuously updating each component and its center frequency using the alternating direction multiplier method, the saddle point of the unconstrained model is finally obtained, which is the optimal solution to the original problem. All components can be obtained from the frequency domain space using the following formula: ； S4.1.5 yes After Wiener filtering, the remaining quantity is used by the algorithm to re-estimate the barycenter frequency based on the barycenter of the power spectrum of each component. The specific process is as follows: a) , , and ; b) Period: ; c) Updated according to S4.1.

4. ; d) Update according to the following formula ; ； e) Update according to the following formula ; ； In the formula: To meet the noise tolerance and signal fidelity requirements for decomposition; and Corresponding to and Fourier transform; f) Repeat steps a~f until the iteration stopping condition is met. The stopping condition is: ； Finally obtained after stopping iteration IMF components and center frequency.

6. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 4, characterized in that, Step 4.3 includes the following: S4.3.1 Based on the six modal components obtained in step S4.2, a correlation matrix {GPR IMF1 IMF2 IMF3 IMF4 IMF5 IMF6} is formed with the original GPR signal. By calculating the ratio of the covariance to the standard deviation between the variables, the influence of dimensions is eliminated, and the standardized linear correlation index is obtained. The formula principle is as follows: ； in, The closer the correlation is to ±1, the stronger the linear correlation; the closer it is to 0, the weaker the linear correlation. and For column vectors in a matrix; S4.3.

2. Based on the correlation between the original signal and the modal components in the correlation matrix, a correlation of less than 0.4 is generally considered to be weak or non-correlated. The strongly correlated modal components are retained as the optimal modal component combination to obtain the decomposed spectrum on the profile. The radar spectrum on the survey line is processed based on steps S4.1, S4.2 and S4.3.

1.

7. The method for determining the spatial continuity of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 4, characterized in that, Step S5 includes the following steps: S5.

1. Based on the radar signal on the soil drill profile obtained in step S4.2, the RMS value is obtained using a sliding window. The gain is normalized by division, and a grayscale image is plotted on the soil drill profile. The amplitude change of the radar signal will cause the color change of the grayscale image. Based on the location of the gleyed layer obtained by the soil drill in step S1.1, the representative part of the gleyed layer in the grayscale image is found by the soil drill. The signal characteristic peak positions representing the upper and lower boundaries of the gleyed layer are screened by combining the grayscale image and the soil drill data to obtain the location signal of the gleyed layer. S5.2 Based on the original ground-penetrating radar data of the survey line obtained in step S2.4, obtain the horizontal resolution of the antenna, that is, the number of radar waves when the measurement distance is 1 meter, set the initial radar wave position, randomly select a radar wave every 10 centimeters of the survey line, and repeat the processing for all survey lines. S5.

3. Based on the radar waves on the survey line in step S5.2, apply a smooth envelope to the single-channel radar wave. Based on the instantaneous phase change position of the radar signal envelope and the signal characteristic peak positions of the upper and lower boundaries of the gryllium obtained in step S5.1, obtain the upper and lower positions of the gryllium, calculate the difference, and obtain the gryllium thickness position. Repeat the process for all survey lines.

8. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 7, characterized in that, Step S6 includes the following steps: S6.

1. Taking one corner of the test area as the origin, the north-south direction as the Y-axis and the east-west direction as the X-axis, assign the plane position coordinates of the survey line, that is, the coordinates (X, Y) of the survey line position relative to the origin. Combined with the upper and lower boundaries Z1 and Z2 of the gleyed layer obtained on the survey line in step S5.3, obtain the spatial position of the upper and lower boundaries of the gleyed layer on the survey line in the three-dimensional coordinate system. S6.

2. Based on the spatial positions of the upper and lower boundaries of the gleyed layer on the survey line obtained in step S6.1, perform inverse distance weight interpolation on the upper and lower boundaries of the gleyed layer on the survey line. The distance weight is 3 meters. Start iterating from 0 and stop iterating when the error is less than the preset threshold to obtain a three-dimensional spatial dataset of the gleyed layer position on the plane. S6.

3. Based on the three-dimensional spatial dataset of the gryllium location obtained in step S6.2, draw a three-dimensional spatial coordinate system. The location of the gryllium is between the upper and lower boundaries of the gryllium.

9. The method for determining the continuous spatial variation of the gleyed layer based on variational mode decomposition and ground-penetrating radar technology according to claim 8, characterized in that, Step S7 includes the following steps: S7.1 Based on the radar survey line laid out in step S1.2, starting from the center of the survey line and extending to both sides, drill a soil drill every 5 meters and number it, record its position on the survey line and the thickness of the gley layer as the observation value, and based on the thickness of the gley layer on the radar survey line obtained in step S5.3, pick the thickness of the gley layer at the corresponding soil drill position as the prediction value. S7.

2. Based on the linear regression model, perform linear fitting between the observed and predicted values, using R0... 2 To evaluate model accuracy, R 2 The higher the precision, the higher the accuracy.

Images