Ground Penetrating Radar (GPR) Identification Method for Agricultural Soil Stratification in a Typical Mollisols Area of Northeast China

In order to achieve a rapid and accurate identification of soil stratification information and accelerate the development of smart agriculture, this paper conducted soil stratification experiments on agricultural soils in the Mollisols area of Northeast China using Ground Penetrating Radar (GPR) and obtained different types of soil with frequencies of 500 MHz, 250 MHz, and 100 MHz antennas. The soil profile data were obtained for 500 MHz, 250 MHz, and 100 MHz antennas, and the dielectric properties of each type of soil were analyzed. In the image processing procedure, wavelet analysis was first used to decompose the pre-processed radar signal and reconstruct the high-frequency information to obtain the reconstructed signal containing the stratification information. Secondly, the reconstructed signal is taken as an envelope to enhance the stratification information. The Hilbert transform is applied to the envelope signal to find the time-domain variation of the instantaneous frequency and determine the time-domain location of the stratification. Finally, the dielectric constant of each soil horizon is used to obtain the propagation velocity of the electromagnetic wave at the corresponding position to obtain the stratification position of each soil horizon. The research results show that the 500 MHz radar antenna can accurately delineate Ap/Ah, horizon and the absolute accuracy of the stratification is within 5 cm. The effect on the soil stratification below the tillage horizon is not apparent, and the absolute accuracy of the 250 MHz and 100 MHz radar antennas on the stratification is within 9 cm. The overwhelming majority of the overall calculation errors are kept to within 15%. Based on the three central frequency antennas, the soil horizon detection rate reaches 93.3%, which can achieve accurate stratification of soil profiles within 1 m. The experimental and image processing methods used are practical and feasible; however, the GPR will show a missed detection for soil horizons with only slight differences in dielectric properties. Overall, this study can quickly and accurately determine the information of each soil stratification, ultimately providing technical support for acquiring soil configuration information and developing smart agriculture.


Introduction
Soil stratification results from migration, transformation, and accumulation of materials within the soil body under physical, chemical, and biological effects during the formation process. Soils have apparent differences in their soil horizons under different natural conditions. Different soil horizons have other regulatory effects on fertility factors such as water, fertilizer, gas, heat, and water and salt transport that affect crop growth and development, which in turn affects crop yield, so soil horizons are an essential indicator for arable land quality evaluation (Song et al., 2018;Santos-Francés et al., 2019). In addition, soil water movement is further complicated by differences in soil texture and factors associated with soil water movement between different soil horizons, which affect the transport properties of soil water and solutes (Li et al., 2014). Accurate acquisition of soil layers can improve the accuracy of soil moisture simulation and enable the scientific and rational use of water resources in farmland. However, it is a challenge to obtain farmland soil stratification information quickly and efficiently when using the model mentioned above for quality evaluation and scientific management of cropland.
Visual soil evaluation (VSE) is a standard means to obtain soil stratification information and is a very effective tool in soil research and agricultural management (Emmet-Booth et al., 2019). VSE relies on hand-dug profiles or the use of hand-held augers, core drills, soil probes or direct observation of natural soil profiles to determine the stratified properties of soils (Morse et al., 2012;Han et al., 2016). Although these methods are more intuitive and do not require complex data interpretation, and can accurately obtain accurate stratification information for individual points, they are timeconsuming and inefficient in investigating soil stratification over large areas, and stratification requires personal experience to make judgments. In addition, soils have spatial variability, and a single sampling point can not represent the characteristics of a particular area. Therefore, efficient, non-destructive, continuous, and accurate soil stratification detection technology is critical. Geophysical methods, especially Ground Penetrating Radar (GPR), have great potential for application in soil and peatland detection and monitoring vegetation bodies. GPR is a non-destructive technique that uses highfrequency radio waves to detect the distribution patterns of subsurface media. It has been applied in environmental, earth exploration, and municipal engineering because of its continuous detection, high resolution, and easy operation (Jol, 2009). One of the main advantages of GPR over VSE is the ability to collect a large amount of profile information quickly and without any destructive effects on the soil.
Since the last century, numerous scholars have researched soil stratification using GPR. In the USA, GPR was applied in 1979 to a soil quality survey in Florida, where sandy soils were predominant, and the stratification structure of the soil horizons was more evident (Johnson et al., 1980). The physical, chemical, and soil moisture differences between the different soil horizons can cause strong reflections in the data recorded by GPR, making it a reliable and effective tool for estimating argillic, spodic, placic, and peat horizons (Collins and Doolittle, 1987;Moorman et al., 2003;Lebron et al., 2004). The thin and spatially heterogeneous soils of the Yili New Reclamation area in Xinjiang, China allow the GPR to clearly and accurately identify the location of the soil and gravel boundary (Yu et al., 2011). Spodic and C horizons have strong reflective properties in GPR images, and the four representative sandy soils of Campinas are distinctly different in GPR images and can be accurately identified (de Mendonça et al., 2014). Soil moisture varies seasonally, and the ground-penetrating radar signal was analyzed in different seasons to identify two different soil stratification boundaries and the best time for ground-penetrating radar detection (Zhang et al., 2014). Ground-penetrating radar has the advantage of being able to detect the thickness of the active horizon of permafrost well. However, retrievals are more challenging in soils with large amounts of gravel (Jafarov et al., 2017). The processing of GPR echo signals based on the envelope detector method and short-time Fourier transform spectral analysis enables accurate stratification of the four soil configurations in tidal soils (Li et al., 2020). Mollisols have a loose surface horizon and accumulate a large amount of humus, making them suitable for crop cultivation. The mollic epipedon thickness of the loess parent material is measured accurately by GPR with an error of between 3.72% and 10.64% . The basic parameters for modeling organic carbon stocks in forest soils require organic (O) and organomineral (A) horizons. GPR can accurately identify the O + A horizon/subsoil boundary under different soil moisture conditions, which is promising for applications (Zajícová and Chuman, 2022). In the existing GPR soil stratification studies, much research has focused on the significant differences between soil horizons and the apparent stratification boundaries between different soil horizons seen in the GPR profiles. However, there are fewer studies on the stratification of agricultural soils. In most studies, only a single center frequency antenna is used for detection, which has a limited detection depth and accuracy. Furthermore, the differences in soil stratification are less evident for agricultural soils in Northeast China. It is inconclusive whether a quantitative representation of the typical soil stratification in the region can be accurately made.
The main objectives of this study are 1) to decode the time-domain features containing soil stratification information in GPR profiles and develop a decoding method based on different central frequency antennas to acquire data; 2) to determine the location of soil stratification in specific agricultural fields in Northeast China based on different frequency antennas combined with different soil horizon dielectric properties; 3) to assess the accuracy of GPR for soil stratification in agricultural fields. Our results provide an effective method for rapidly acquiring soil stratification information and technical support for conserving and assessing farmland soils in the Mollisols region of Northeast China.

Overview of the study area
The study area is located in the Youyi Farm of Shuangyashan City and the Shuguang Farm of Huainan County, Jiamusi City, Heilongjiang Province, in the core of the Sanjiang Plain. The landscape of the Sanjiang Plain is vast and low, with elevations ranging from 10 to 1632 m, with most areas below 200 m. The weather is consistent with a mid-temperate monsoon climate, with an annual rainfall of about 350 to 800 mm (Chen et al., 2018;Zhang et al., 2020), and the soil types are mainly Phaeozems, Luvisols, Cambisols, and Gleysols. The geology is characterized as a basin formed by Tertiary argillaceous subsidence on a basement composed of pre-Paleogene metamorphic, Paleozoic, and Mesozoic sedimentary rocks . Youyi Farm is located between 131°28′E-132°15′E and 46°31′N-46°59′N, with a total area of about 1696 km 2 and an annual average temperature and average annual precipitation of 3.1°C and 514 mm, respectively (Guo et al., 2023). The main crops are corn, soybeans, and rice. Shuguang Farm is located between 129°56′E-130°38′E and 46°14′N-46°23′N, with a total area of about 177 km 2 and a coldtemperate continental monsoon climate, with an annual average temperature of 3.6°C and an average annual precipitation of 537 mm (Liu et al., 2019), and the main crop cultivation types are maize and soybean. The experimental sampling locations are shown in Fig. 1.

Data acquisition
In Northeast China, the temperature is often below 0°C in winter, and the moist soil is frozen, so the difference in dielectric properties between soil horizons is slight.  Fig. 1 Location of Northeast China and sampling sites. Panel a shows the location of the study area, while b and c show the location of the sampling sites at Youyi Farm and Shuguang Farm, respectively GPR can not operate during the crop growth period, so it usually conducts experiments after the autumn harvest or before the spring planting. In addition, rainfall significantly impacts soil surface moisture, which generally returns to the level before the downpour after seven days of rain, so we choose to conduct GPR experiments after a week of continuous clear weather. Based on the above considerations, we conducted a field experiment in the study area at the end of October 2021 to collect ground-penetrating radar data, soil profile stratification, and soil temperature and moisture data for several soil types commonly found in the study area. The penetration depth of ground-penetrating radar in the soil is related to the antenna center frequency; the higher the frequency, the deeper the penetration depth, but the lower the resolution. The soil moisture in the northeastern region is high, and the penetration depth of GPR is limited, so the antennas with three center frequencies can meet the detection requirements of different depths. Ground-penetrating radar includes three detection methods: profile method, common center point method, and transmitted wave method (Daniels, 1996). This study used the most commonly used profiling method for detection. The ground-penetrating radar used in this experiment was a ProEx model system unit manufactured by MALA GeoScience, Sweden, with antenna frequencies of 500 MHz, 250 MHz, and 100 MHz. The sampling frequency of the 500 MHz shielded antenna was set to 8000 MHz, and the sampling time window was set to 40 ns; the sampling frequency of the 250 MHz shielded antenna was set to 3000 MHz, and the sampling window was set to 100 ns; the sampling frequency of the 100 MHz shielded antenna was set to 1200 MHz, and the sampling window was set to 200 ns. Round-trip measurements were made on typical soil types in the study area using a ranging wheel trigger, with a line length of 20 m for each size. GPR reduced multiple reflections of the signal when leveling the ground for data collection, and the land was relatively level after autumn tilling, in which the tilling depth was 35-40 cm, and the harrowing depth was about 16 cm.
To verify the accuracy of the ground-penetrating radar (GPR) soil stratification and to collect soil moisture data at the corresponding stratification locations, excavation of the soil profiles was carried out in conjunction with the GPR data collection, with soil science experts on-site to stratify the soil profiles and identify the stratification locations. Four soil subclasses were collected at Youyi Farm, including Albic Phaeozems, Phaeozems, Cambisols, and Calcic Cambisols; additionally, Luvisols were collected at Shuguang Farm, yielding five soil subclasses in total. Soil profiles were excavated to a depth of 1 m, and five soil profiles were excavated (Fig. 2). Soil moisture was collected at 20 cm intervals using a Pro Check handheld multifunctional readout meter from DECAGON, USA.

Test methodology 2.3.1 Main processes for the identification of soil stratification information
Two primary parameters are required to obtain the location of soil stratification using GPR: the time-domain location for each soil stratification and the propagation velocity of the GPR signal in each stratification. To determine the time-domain location of soil stratification, it is necessary to pre-process the GPR signals obtained from the three frequency antennas, then perform a multilayer decomposition of the pre-processed signals using one-dimensional discrete wavelet analysis to reconstruct the high-frequency information of each layer, and finally determine the time-domain location of the stratification based on the envelope signal for the reconstructed signals. For the acquisition of the electromagnetic wave propagation velocity in each horizon of the soil, the GPR signal propagation velocity was obtained from the soil moisture data collected in the field in each hori-

Fig. 2
Soil profiles at the sampling site in Northeast China zon combined with the classical Topp formula to obtain the dielectric constant of the soil in each horizon. The specific experimental procedure is shown in Fig. 3. The data collected by GPR in the field contains noise, which needs to be filtered and enhanced to improve the signal-to-noise ratio and get a clearer image reflecting the subsurface stratification information. The specific steps of data pre-processing are as follows: 1) de-DC drift, eliminating the DC component or DC offset in the radar signal; 2) removing the start time, as the antenna is a certain distance from the surface, the zero moments of the GPR profile needs to be adjusted to the soil surface so that the multi-channel data can obtain a uniform time zero point; 3) energy attenuation gain, mainly to overcome the energy attenuation of electromagnetic waves in the propagation process. The GPR signal is a broadband signal and wavelet analysis has multi-scale characteristics. Using wavelet analysis to reconstruct high frequency information at different scales, then we obtain accurate stratification time-domain position based on envelope signal at different scales, before finally calculating the soil stratification thickness information by combining each stratification wave speed.

One-dimensional discrete wavelet analysis of GPR signals
GPR captures the variation of the subsurface medium with spatial location by receiving the reflected signal when electromagnetic waves encounter different media where the amplitude intensity, phase, and frequency of the reflected signal change accordingly. Differences in the physical, chemical, and mineral characteristics between different horizons of the soil cause changes in the instantaneous parameters of the reflected signal. The GPR signal is a broadband signal, and the signals in each frequency band affect each other, so direct timefrequency analysis of broadband signals is susceptible to interference from high-frequency signals, which is not conducive to interpreting the results. One-dimensional discrete wavelet analysis can decompose the signal at multiple scales, decomposing the signal within each frequency band, allowing both overall signal information and detailed information to be seen, as shown in Eq. (1) (Jevrejeva et al., 2003;Furon et al., 2008).
where DWT (m, n) is the wavelet coefficient of the GPR signal, a 0 and b 0 are constants that depend on the choice of wavelet basis, m is the scale factor and the translation, n is the wavelet basis, ψ(m, n) is the wavelet basis, the superscript * denotes the complex conjugate of a complex number, the sum in the equation is taken over all integer values of k, representing the different translations of the wavelet function,and f is the GPR data.
Wavelet analysis divides the original signal into highfrequency and low-frequency information at different scales. Since electromagnetic waves passing through different soil horizons will have instantaneous frequency jumps, i.e., manifesting as high-frequency information, the high-frequency coefficients at different scales are reconstructed to obtain high-frequency information containing layered information. The wavelet bases need to be selected before wavelet analysis can be carried out. After repeated trials, the db4 wavelet base was determined suitable for the GPR data obtained using 500 MHz and 100 MHz radar antennas. The coif4 wavelet base was suitable for GPR data obtained from the 250 MHz radar antennas.

Soil hierarchy based on envelope signals
The GPR signal reconstructed by wavelet analysis is the high-frequency information in the data. It has a narrow frequency band, which can better reflect the details of the original signal at different scales, but it also contains noise. The envelope signal can retain the waveform characteristics of the signal well, can demodulate  the low-frequency signal containing soil stratification information in the high-frequency signal, and requires a lower signal-to-noise ratio. The Hilbert transform is an essential tool for analyzing and processing signals, and the signal's instantaneous parameters can be obtained after the Hilbert transform. For the acquired envelope signal, the Hilbert transform is performed: where is the convolution of f(t) with 1/πt, τ is the intermediate integral variable introduced by the convolution operation and t is the travel time of the GPR signal.
In turn, the resolved signals of the and f(t) construction envelope signals are obtained as follows: Then, the instantaneous amplitude of the envelope signal a(t) can be expressed as: The instantaneous phase of the envelope signal, θ(t), can be expressed as: The instantaneous frequency ω(t) of the envelope signal is the time rate of change of the phase, and the instantaneous frequency can be expressed as: The instantaneous frequency is the temporal rate of phase change. It will change significantly when the electromagnetic waves pass through different media partitions. In this way, it can indicate changes between the different soil horizons and helps identify the soil partition boundaries.

Soil dielectric constant and estimation of the electromagnetic wave velocity
Soil dielectric properties differ due to the difference in physio-chemical properties between soil layers. Since the dielectric loss of soils is low and the soil dielectric constant is mainly related to the volumetric water content, the relative dielectric constants of the corresponding soil horizons were obtained by applying the Topp relation (Topp et al., 1980): where ε r is the relative permittivity and θ v is the volumetric water content percentage of the medium. For low-loss soils, the speed of electromagnetic wave propagation in the soil can be expressed by the following equation (Knight, 2001): where v is the speed of electromagnetic wave propagation in soil, c is the speed of electromagnetic wave propagation in a vacuum, c = 0.3 m/ns, and ε r is the relative permittivity.

GPR detection accuracy and determination of soil stratification location
The GPR resolution consists of two parts: the lateral resolution is the angular resolution, and the vertical resolution is the distance resolution. The main focus of this experiment is on soil stratification information, so only the vertical resolution is considered. In GPR applications, where one-fourth of the wavelength is generally used as the lower limit of the longitudinal resolution, the resolution of GPR in the soil can be expressed by the following equation: where ∆r is the resolvable longitudinal length; λ is the wavelength at which the electromagnetic wave propagates through the soil. Based on the time-domain position of the soil stratification obtained by the above method and the speed of electromagnetic wave propagation in each horizon, the stratification position of the soil can be obtained as follows: where h is the soil stratification position; v is the propagation speed of the electromagnetic wave in the soil; t is the two-way time-domain position of the electromagnetic wave at the soil stratification position.

Comparative analysis before and after wavelet analysis
For the five soil subcategories sampled, this paper uses wavelet analysis to obtain each layer's reconstructed high-frequency signal profiles (Fig. 4). Since the processing effect of each sample point is similar, one example is used to illustrate the effect before and after processing. Fig. 4a shows the instantaneous frequency profile based on a 500 MHz antenna after data pre-processing, and Figs. 4b-4h show the instantaneous frequency profile of the high-frequency signals reconstructed by wavelet analysis for layers 1 to 7. Fig. 4i shows the instantaneous frequency profile of the fifth layer reconstructed by wavelet analysis based on the 250 MHz antenna. Fig. 4j displays the instantaneous frequency profile of the fourth layer reconstructed by wavelet analysis based on a 100 MHz antenna. As shown in Fig. 4a, the high-frequency signal increases gradually with increased time depth. With the increase of time depth, the electromagnetic wave scattering process becomes more complex, resulting in more signal noise received. The higher the center frequency of the antenna, the shallower the penetration depth. In addition, the GPR signal is a broadband signal, and the noise signal is mixed with helpful information. The transient frequency mutation signal of the shallow surface layer will be annihilated, and the transient frequency change can not accurately determine the soil stratification boundary. In this paper, the wavelet analysis of frequency separation characteristics is used to reconstruct the high-frequency signal of each layer to obtain highfrequency signal profiles at different scales. From Figs. 4c-4h, it can be seen that as the scale of signal decomposition reconstruction increases, the high-frequency noise signal is continuously removed, and clear stratification information is obtained for the shallow surface layer, as shown in Fig. 4g. Fig. 4h shows the highfrequency signal reconstructed in the seventh layer, from which no abrupt information can be seen. This is because the high-frequency signal containing the stratification information is removed entirely. Therefore, the high-frequency signal reconstructed in the sixth layer is used as the final soil stratification boundary. Based on the same method, the data obtained based on the 250 MHz antennas and the 100 MHz antennas were processed, and the instantaneous frequencies of the fifth reconstructed layer of the 250 MHz antennas (Fig. 4i) and the fourth reconstructed layer of the 100 MHz antennas (Fig. 4j) were finally determined to give more explicit soil stratification boundaries. The instantaneous frequencies obtained directly from the reconstructed signals are susceptible to interference from neighboring signals, resulting in significant period domain positions of the stratification boundaries and weak local signals (Fig. 4i). Therefore, further processing of the discontinuous and inconspicuous stratification boundaries (Figs. 4g, 4j) is required to the reconstruct signals from the sixth layer to enhance the display and weaken the interference signals to obtain more precise and more accurate stratification boundaries.

Analysis of soil hierarchy results based on envelope signals
Based on wavelet analysis and reconstruction results, the corresponding frequency GPR signal is enveloped, and the low-frequency signal containing hierarchical information is demodulated from the high-frequency signal. To more clearly show the effect of taking the envelope line and analyzing the causes of discontinuity and unclear layering, the discontinuous single channel data at the layering position obtained by the three frequency antennas and the unclear layering data are exported. The time-frequency diagram before and after the singlechannel data processing of Luvisols is obtained (Fig. 5).   For the data acquired by the 500 MHz and 100 MHz antennas, a db4 wavelet basis is used. For the data acquired by the 250 MHz antenna, a coif4 wavelet basis is used, which has better symmetry than db4, resulting in a different denoising effect. 250 MHz GPR signals are smoother than 500 MHz and 100 MHz, and the envelopes behave more smoothly and compactly. From Figs. 5b and 5h, the instantaneous frequencies obtained from the direct wavelet analysis reconstructed signals show a random nature, with no obvious variation pattern and no information on the signal jumps at the stratification locations, thus causing discontinuities in the stratification boundaries of the instantaneous frequency profiles obtained after wavelet analysis. The instantaneous frequencies after taking the envelope show clear jump information in Figs. 5b and 5h. Fig. 5e shows that although the signal instantaneous frequency jump is shown, the signal instantaneous frequency jump has a more extended period and changes more slowly, followed by a wider stratification boundary and unclear stratification boundary in the instantaneous frequency profile. After solving the instantaneous frequency for its envelope signal, it can be seen in Fig. 5f that the time of the jump information narrows and becomes more prominent. Overall, the envelope signal's instantaneous frequency can more accurately identify the signal's jump time-domain location.

Soil dielectric properties and detection accuracy analysis
Soil dielectric properties remained almost constant from 0 to 20°C (Du et al., 2022). The five soil subclasses showed little change in temperature with increasing profile depth (Fig. 6a), so the effect of soil temperature on the dielectric constant was small. As the profile's depth changed, the moisture of each soil subclass differed significantly (Fig. 6b). The dielectric constant was observed to change accordingly along with the changing soil moisture (Fig. 6c), and the difference in dielectric constant between soil layers is helpful to distinguish the changes between different layers. In addition, differences between dielectric constants can result in different speeds of electromagnetic wave propagation between soil horizons. Changes in the dielectric properties of each soil subclass in different horizons also led to changes in longitudinal resolution. Taken together, the longitudinal resolution of each soil subclass decreases with decreasing central frequency. The longitudinal resolution is about 0.03 m for an antenna transmitting a center frequency of 500 MHz (Fig. 6d), about 0.06 m for an antenna transmitting a center frequency of 250 MHz (Fig. 6e), and about 0.16 m for an antenna having a transmitting center frequency of 100 MHz (Fig. 6f). For shallow surface soils with thin layers, it is more difficult to distinguish the differences between different layers for low-frequency signals.

Analysis of soil stratification results
Based on the results of wavelet analysis reconstruction, the high-frequency information reconstruction signals of five soil subclasses at three frequencies were taken as envelopes to obtain their instantaneous frequencies. Finally, the corresponding wave velocities were obtained with the Topp formula to obtain the GPR soil stratification and field soil profiles, as shown in Fig. 7.
The stratification positions are generally more accurate for Luvisols, Phaeozems, and Albic Phaeozems. However, for meadow soils, the accuracy of the stratification positions was relatively poor due to soil sesquioxides and salinity. For Luvisols (Fig. 7a4), the actual stratification positions were at 25 cm and 52 cm. The accuracy of the stratification position data obtained for the 100 MHz antenna (Fig. 7a3) could not be judged because the profile was excavated to a depth of 1 m. However, the stratification positions were more accurate for the data obtained based on the 500 MHz and 250 MHz antennas, which were stratified at approximately 29 cm (Fig. 7a1) and 56 cm (Fig. 7a2). Based on the triple frequency antenna data, the stratification positions for the Phaeozems (Fig. 7b4) are more accurate at 35 cm, 52 cm, and 97 cm. For Albic Phaeozems (Fig. 7c4), the stratification results obtained from the 500 MHz ( Fig. 7c(1)) and 250 MHz (Fig. 7c2) antennas were more accurate, with some breaks in the 100 MHz stratification boundaries (Fig. 7c3), due to the lower vertical resolution at 100 MHz, which resulted in some of the boundaries not being identified. For meadow soils (Fig. 7d4), the soil profile was stratified at 20 cm, 40 cm, and 90 cm. The GPR profile did not identify the stratification at 40 cm, noting that 20-40 cm is a transitional layer containing a small number of rust spots, with a slight variation in physical and chemical proper-ties, and the difference in dielectric constants between its layers is slight, hence the missed detection phenomenon. For the Calcic Cambisols (Fig. 7e4), the soil profile is stratified at 20 cm, 33 cm, and 71 cm. The actual stratification at 33 cm differs significantly from the results of the GPR stratification, as can be seen from the 250 MHz antenna stratification results (Fig. 7e2), which is approximately 45 cm, mainly because of the medium amount of lime spots contained in this layer. The stratification limits of the 100 MHz antenna show some fluctuations (Fig. 7e3), mainly due to the clay loam texture of the soil, which has the property of adsorbing ions, and the small number of lime spots, medium rust spots, and calcareous nodules in this layer, ultimately resulting in an uneven stratification.
As a whole, as the antenna center frequency de-creases, the depth of GPR detection becomes more profound, and as the depth increases, the signal noise increases, and the accuracy of the soil stratification location decreases. To quantitatively represent the accuracy of the test method, the measured data of the five soil subclasses were compared with the test result data (Table 1) to obtain the absolute and relative errors for each soil type, with the absolute error equivalent to the difference between the calculated and measured values and the relative error equivalent to the percentage of absolute error to the measured value. From the final results, the maximum relative error was 25.67%, and the absolute error remained within 9 cm. The stratification accuracy based on the 500 MHz radar antenna was the highest and gradually decreased as the center frequency of the radar antenna decreased.  Fig. 6 Plots of the data variation with soil depth for the five subclasses. Panels a-c show the variation of temperature, moisture, and dielectric constant with soil profile depth for each of the five soil subclasses, respectively, and panels d-f show the longitudinal resolution of the three center frequency antennas with soil profile depth for each of the five soil subclasses

Discussion
The electromagnetic wave propagation velocity is determined by the dielectric constant, which is influenced by soil moisture and soil physicochemical properties, and the accurate acquisition of electromagnetic wave propagation velocity can directly affect the accuracy of stratification information. Due to the complex soil en-  (4) shows the stratified profiles for each soil subclass vironment in agricultural fields, the sources affecting the dielectric constant are different for different soil subclasses. This experiment was conducted at the end of October when the cultivated land in the northeast was plowed up to a depth of 35-40 cm and harrowed to a depth of about 16 cm when autumn tilling was performed. This will reduce the surface soil macropore space and soil moisture evaporation, and allow the lower layer to rise, increasing the water content of the surface soil and reducing the speed of electromagnetic wave propagation. In addition, to protect the ecological environment and realize the perpetual use of black soil resources, the straw return policy has been implemented in Northeast China. However, due to the low temperature in Northeast China, the decomposition rate of straw is slower, thus causing the soil to be mixed with plant residues (Fig. 7b), which in turn affects the propagation speed of electromagnetic waves in the soil (Pan et al., 2012); when conducting the profile soil moisture data collection, we collected every 20 cm on the profile, but because soil moisture is a dynamic variation on the profile, soil moisture within 20 cm was con-sidered as homogeneous in the actual calculation of wave velocity. The thickness estimation by Zajicova et al. using the average dielectric constant of soil is less accurate, which is due to the spatial variability of soil and will bring some errors, which is similar to the results of this experiment; as the calcium carbonate content of the soil increases, the soil dielectric constant also increases. In addition, soil iron oxides have a larger specific surface area and water retention capacity than quartz sand, which can reduce radar wave velocity and cause reflective properties of electromagnetic waves in soil (Van Dam and Schlager, 2002;Lebron et al., 2004). Meadow soils and calcareous meadow soils contain these substances, which can lead to deviations between the calculated and actual values of electromagnetic wave velocity and relatively low stratification accuracy (Table 1). The discontinuities of the instantaneous frequencies at the three frequency points after wavelet analysis are caused by the low-frequency resolution of wavelet analysis. Soil background is complex, and the primary data processing process is challenging to obtain the location of soil stratification. Wavelet analysis has the character- Notes: In a typical soil profile, A represents the topsoil horizon, also known as the 'surface horizon'. Ap represents the plow horizon, while Ah represents the humus surface horizon. B represents a horizon of material accumulation, and C represents the parent material horizon. E represents the eluviation horizon. The AB and BC horizons are transitional zones between adjacent horizons, characterized by changes in soil properties such as color, texture, or mineral composition istic of frequency separation, and the analysis of data in different frequency bands can effectively identify the information of structural properties hidden in the data (Fig. 4). However, wavelet analysis has low time-frequency resolution, and some critical soil stratification information sometimes can not be presented, which is consistent with the view of He et al. (2019). The GPR signal after wavelet analysis was taken as an envelope to obtain more obvious stratification information in the low-frequency signal. Li et al. analyzed the GPR signals of four tidal soil types. Finally, they obtained the stratification information of the four soil types by taking the envelope of the preprocessed GPR signals and obtaining the location of the transient phase jump included in the envelope signal, and their stratification results were similar to the results in Table 1. GPR relies on the signals reflected from different dielectric layers to identify the target, and when the difference in dielectric properties between different soil horizons is slight, it will not be detected. The AB layer of Cambisols showed a missed detection phenomenon in this experiment. Winkelbauer et al. (2011) found that a clear boundary line between the Ah and B layers is usually not shown but rather a fuzzy transition zone, similar to the results in this paper (Fig. 7d). Soil thickness in the Mollisols region of Northeast China determines soil productivity and is essential to national food security (Liu et al., 2023). However, the available information on the spatial variation of black soil thickness is insufficient. Therefore, obtaining accurate soil thickness for national food security is vital. The traditional method to get soil thickness information is by excavating profiles, which could be more efficient and can only obtain point information. Zhang et al. (2021) performed interpolation based on limited soil sample profiles, and this method has a significant bias due to the large variety of soil thicknesses at different scales. GPR data acquired based on three center frequency antennas can accurately obtain soil thickness information after a series of data processing. For farmland plowing, the depth of deep plowing of the land does not exceed 40 cm at most and based on the results of this experiment, the accuracy can be satisfied. In addition, for the low-yielding soil white pulp soil, the main obstacle level affecting crop growth is the white pulp layer below the soft soil layer. Therefore, accurately acquiring the white pulp layer's location is essential for treating this soil. From the experimental results, the GPR based on 500 MHz and 250 MHz antennas can accurately acquire the barrier level of the white pulp soil (Fig. 7a).
These factors affecting the electromagnetic wave velocity will be fully considered in the next step of the work to obtain a more accurate stratification accuracy.

Conclusions
In this paper, for the experimental study of the stratification of the main soil types in the typical Mollisols area oof Northeast China, GPR data were pre-processed using traditional data processing methods. The time-domain location of the soil stratification was determined using wavelet analysis combined with the envelope method. Finally, the electromagnetic wave propagation velocity at the corresponding stratification location was calculated in combination with the Topp formula to determine the location of the soil stratification, and the following conclusions were obtained.
(1) The ground-penetrating radar technique can be used to determine the location of soil stratification in the primary soils of the typical Mollisols area. The accuracy of identification of different soil stratification varies. The maximum relative error is 26%.
(2) Using wavelet analysis to reconstruct the high-frequency information by combining the envelope signal in the high-frequency information to highlight its stratification information, which can accurately show the soil stratification location. Combined with the field stratification results, it is basically consistent with the actual soil stratification. GPR can miss detection for soils with slight differences in nodal properties between horizons.
(3) The high moisture content of the soil in the northeast region results in a greater degree of attenuation of electromagnetic waves and limited penetration depth. This experiment combines three center frequency antennas to enable detection at different depths and accurate detection of soil horizons at a depth of 1 m.
The research results provide technical support for the rapid determination of the thickness of the Mollisols horizon, the soil configuration of cultivated soils in the Mollisols area, and supports intelligent agriculture needs such as soil nutrient estimation and guidance on variable fertilization.