A method for measuring the bulk density of topsoil based on the characteristics of sound wave penetration

Abstract

A method for measuring the bulk density of topsoil based on the sound wave penetration characteristics belongs to the field of intelligent agricultural sensing technology. In the sound wave penetration characteristic acquisition system of this invention, a digital power amplifier board is electrically connected to a subwoofer, a DC power supply, and a computer terminal; a data acquisition card is electrically connected to a microphone and the computer terminal; and a robotic arm is electrically connected to the DC power supply and the computer terminal. The robotic arm, through five servo motors in coordinated control, achieves automatic grasping, precise positioning, and vertical insertion of the ring cutter tube, ensuring consistent sample placement. The computer terminal is equipped with a soil bulk density measurement model based on an adaptive mutant particle swarm optimization variational mode decomposition model and an adaptive mutant particle swarm optimization BP neural network. This system can automatically acquire real-time soil penetration sound waves and realize soil bulk density inversion based on sound wave penetration characteristics. The soil bulk density measurement model of this invention has high timeliness and high accuracy.

Description

Technical Field

[0001] This invention belongs to the field of intelligent agricultural sensing technology, specifically relating to a method for measuring the bulk density of topsoil based on the penetration characteristics of sound waves. Background Technology

[0002] Exploring soil bulk density measurement is of significant importance for agricultural production, environmental protection, and soil resource management. Currently, soil bulk density determination is divided into direct and indirect methods. Direct methods often employ core sampling (ring sampler) to determine soil bulk density, calculating the bulk density based on the mass difference between dried and undried soil. However, this method suffers from poor timeliness and labor-intensive operation, making it unsuitable for real-time rapid soil bulk density detection. To address these shortcomings, some studies have used mid-infrared spectroscopy to invert soil bulk density (SHI LN, O'ROURKE S, BACHION F, et al. Prediction of soil bulk density in agricultural soils using mid-infrared spectroscopy[J].Geoderma, 2023, 434: 116487.). This method offers advantages such as high accuracy, but it requires soil drying, resulting in high cost and low timeliness. Therefore, designing a timely soil bulk density measurement device is a pressing technical challenge that needs to be addressed.

[0003] Soil bulk density is a three-dimensional characteristic, and general image information cannot characterize it. Some scholars have used acoustic wave penetration to determine soil information (SOMG KL, NIE J, LI Y, et al. Regional soil water content monitoring based on time-frequency spectrogram of low-frequency sweptacoustic signal[J]. Geoderma, 2024, 441: 116765.). This allows for the inversion of soil bulk density using the differences in the propagation properties of acoustic waves in soils with different structures. This process involves key technologies such as speech signal processing models and machine learning model construction. The parameters of the above models have a significant impact on model performance, but the performance of current parameter optimization algorithms needs improvement. Therefore, there is an urgent need to propose a high-performance acoustic wave penetration characteristic analysis method for soil bulk density inversion to improve the efficiency of soil bulk density calculation. Summary of the Invention

[0004] To address the shortcomings of existing technologies, this invention provides a method for measuring the bulk density of topsoil based on the characteristics of sound wave penetration, comprising the following steps:

[0005] S1 is equipped with a sound wave penetration characteristic acquisition system;

[0006] S2 collects n sets of acoustic wave data W;

[0007] S3 uses the AMPSSO-VMD model to reduce noise in n sets of acoustic data W acoustic signals;

[0008] S4 extracts the Mel-frequency inverse coefficients of the n groups of noise-reduced Wj acoustic signals;

[0009] S5. Establish the AMPSSO-BP neural network model;

[0010] S6 Real-time measurement of soil bulk density in the topsoil layer.

[0011] The acoustic wave penetration characteristic acquisition system described in S1 consists of an acquisition box A, a computer terminal 1, a data acquisition card 2, a digital power amplifier board 3, a DC power supply 4, and a robotic arm device C. The acquisition box A comprises a base assembly B and a top cover 5, with a pull ring 6 at the center of the top cover 5. The base assembly B consists of a base 7. The base 7 consists of sound insulation cotton I8, conduit I9, subwoofer 10, sound insulation cotton II11, ring cutter tube 12, sound insulation cotton III13, sound insulation cotton IV14, pickup 15, and conduit II16. The base 7 is a rectangular box structure. The ring cutter tube 12 is fixed to the center of the base plate; the cross-section of the ring cutter tube 12 is square. The pickup 15 is fixed to the rear of the rear plate of the ring cutter tube 12. The front end of the conduit II16 is connected to the rear end of the pickup 15, and the rear end of the conduit II16 is connected to the center of the rear plate of the base 7. The subwoofer 10 is fixed to the front of the front plate of the ring cutter tube 12. The rear end of the conduit I9 is ​​connected to the front end of the subwoofer 10, and the front end of the conduit I9 is ​​connected to the center of the front plate of the base 7. Sound insulation cotton I8 is placed in the space between the front plate of the ring cutter tube 12 and the front plate of the base 7. Sound insulation cotton IV14 is placed in the space between the rear plate of the ring cutter tube 12 and the rear plate of the base 7. The left plate of the ring cutter tube 12 is connected to the left plate of the base 7. Sound insulation cotton Ⅲ13 is placed in the space between the plates; sound insulation cotton Ⅱ11 is placed in the space between the right plate of the ring tube 12 and the right plate of the base 7; the digital power amplifier board 3 is electrically connected to the subwoofer 10, the DC power supply 4 and the computer terminal 1 respectively; the servo motor group 20 of the robotic arm device C is electrically connected to the computer terminal 1; the data acquisition card 2 is electrically connected to the microphone 15 and the computer terminal 1 respectively; the computer terminal 1 is equipped with a soil bulk density measurement model; the top cover 5 and the base 7 are made of acrylic material; the robotic arm device C consists of a base 17, a multi-degree-of-freedom robotic arm body 18, an end effector 19 and five servo motors of the servo motor group 20, wherein the base 17 is a cylindrical structure, the lower end of the multi-degree-of-freedom robotic arm body 18 is fixed to the base 17, and the upper end of the multi-degree-of-freedom robotic arm body 18 is connected to the end effector 19 through one servo motor of the servo motor group 20.

[0012] The acquisition of n sets of acoustic wave data W mentioned in S2 specifically refers to:

[0013] Using a ring cutter, n groups of soil samples with ring cutters are collected. The sound insulation cotton IV (14) and the top cover (5) are pulled out by the pull ring (6). Each soil sample is placed into the ring cutter tube (12) in sequence. Then, the sound insulation cotton IV (14) and the top cover (5) are put back into place by the pull ring (6). The computer terminal (1) emits a sinusoidal 200 Hz step sound wave signal. The digital power amplifier board (3) receives the sound wave signal and drives the subwoofer (10) to emit a step sound wave. The pickup (15) receives the sound wave data after passing through the soil sample. The sampling frequency of the data acquisition card (2) is set to 50000 Hz. The data acquisition card (2) collects the 0.1s sound wave data W of n groups of pickups (15) through the digital power amplifier board (3) and transmits it to the computer terminal (1).

[0014] S3 describes the use of the AMPSSO-VMD model to denoise n sets of acoustic data W acoustic signals, including the following steps:

[0015] S3.1 Set the mode function obtained after signal decomposition to I k I k Time Series I k The mathematical expression for (t) is:

[0016] ;

[0017] Among them: A k (t) is the amplitude, and A k (t)≥0;δ k (t) is a non-monotonic decreasing phase function;

[0018] S3.2 Initialize the VMD algorithm iteration count N, Lagrange multiplier λ, quadratic penalty factor α, and decomposition level K;

[0019] S3.3 The current number of VMD parameter iterations is m. Let m = m + 1, and the VMD algorithm will perform iterative calculations.

[0020] S3.4 Make the value of k continuously increase from 1 to K, I k ω k The update method is as follows:

[0021] ;

[0022] Where: Im+1 k is I k The updated value; IFTD<·> is the operator for Fourier transform and differentiation; ω k For I k The frequency center of (t); ωm+1 k is ω kThe updated value; ω is the variable obtained by performing a Fourier transform on t;

[0023] S3.5 Updates λ, and the update method for λ is as follows:

[0024] ;

[0025] Where: λ m+1 The updated value of λ; IHT<·> is the inverse Hilbert transform;

[0026] S3.6 Repeat steps S3.3 to S3.5 until the termination condition is met:

[0027] ;

[0028] Where: ε is the convergence criterion, and ε > 0;

[0029] S3.7 The parameters α and K of VMD are determined using the AMPSSO algorithm, including the following steps:

[0030] S3.7.1 Determine the particle dimension in the AMPSSO algorithm. The VMD parameters to be optimized are α and K, therefore the particle dimension is 2:

[0031] S3.7.2 Determine the fitness function of the AMPSSO algorithm and calculate the fitness of each individual. The mathematical expression for the individual fitness function, Fitness, is as follows:

[0032] ;

[0033] Where: PE<·> is the permutation entropy operator; m l This is the l-th signal component after VMD;

[0034] S3.7.3 Initialize the individual positions and velocities in the AMPSSO algorithm;

[0035] S3.7.4 Perform population individual velocity updates;

[0036] S3.7.5 Update the location of individuals in the population;

[0037] S3.7.6 Assign the result of the second position update of the individual to the parameters α and K of the VMD algorithm. When the fitness value generated by the current population iteration update is less than the fitness value generated by the previous generation population iteration update, update the individual extreme value and the population extreme value; otherwise, proceed to the termination condition judgment.

[0038] S3.7.7 If the population update iteration number meets the termination condition, stop updating and obtain the optimal parameters α and K; otherwise, repeat steps S3.1 to S3.7.6 to continue updating parameters α and K.

[0039] S3.8 Repeat steps S3.1 to S3.7.7 for the n groups of signals in step S2 to obtain the optimal parameter α. best K best and α best K best Assigning values ​​to the VMD model yields the AMPSSO-VMD model, α best K best The calculation method is as follows:

[0040] ;

[0041] Where: round<·> is the rounding operator; α W(i) The optimal quadratic penalty factor parameters for the AMPSSO-VMD model of the i-th sound wave; K W(i) The optimal decomposition layer parameter for the AMPSSO-VMD model of the i-th sound wave;

[0042] S3.9 The AMPSSO-VMD model obtained in step S3.8 is used to decompose the W acoustic signal to obtain K signal components. The correlation coefficient between each component and W is calculated, and the signals with correlation coefficients higher than 0.6 are summed to obtain W. j .

[0043] S4 describes the extraction of Mel-frequency inverse coefficients from the n groups of noise-reduced Wj acoustic signals, which includes the following steps:

[0044] S4.1 to W j The mathematical expression for pre-emphasis processing of acoustic signals is:

[0045] ;

[0046] Among them: W pe (p) is W j The signal value after pre-emphasis at position p; W j (p), W j (p-1) represent W j Signal values ​​at positions p and p-1; η is a constant close to 1;

[0047] S4.2 to W pe (p) The acoustic signal is framed, with each frame being F in size. size For 1250 data points, frame shift F move With 500 data points, 9 sets of frame signals are obtained. The mathematical expression is:

[0048] ;

[0049] Where: F(q) is the q-th frame; zeros(1,250) is a 1×250 all-zero matrix;

[0050] S4.3 Windowing is applied to F(q) to avoid signal boundary effects; a Hamming window is used.

[0051] ;

[0052] Where: H(F(q)) y The signal value at point y after windowing;

[0053] S4.4 performs a Fast Fourier Transform on each windowed signal and then converts the amplitude spectrum into a power spectrum to enhance the spectral characteristics.

[0054] S4.5 A set of Mel filters is used to filter the power spectrum to improve the resolution in the low-frequency region, resulting in n sets of 40×9 arrays. The Mel filter bank contains 40 Mel filters, and the Mel frequency and frequency conversion in the Mel filters are as follows:

[0055] ;

[0056] Where: Mel is the Mel frequency, and f is the frequency;

[0057] S4.6 Convert the values ​​of the array in step S4.5 into decibel values, then perform discrete cosine transform on the decibel values, extract the first 13 rows of data, and obtain n sets of feature vectors containing 13×9=117 feature values.

[0058] The establishment of the AMPSSO-BP neural network model described in S5 includes the following steps:

[0059] S 5.1 Determine the true bulk density of n groups of soil samples with ring cutters; dry the soil samples in an oven at 105℃ for 12 hours. The mathematical expression for the true bulk density is:

[0060] ;

[0061] Where: m w The mass of the soil sample with the ring cutter before drying, m d V represents the mass of the dried soil sample with the ring cutter. r For the volume of the ring cutter;

[0062] S5.2 Select the feature values ​​in the feature vector from step S4.6 that have a high correlation with BD, calculate the Pearson coefficient between each feature value in the feature vector from step S4.6 and BD, and select z feature values ​​with a Pearson coefficient greater than 0.5 to form n new feature vectors v. f ;

[0063] S5.3 Establish a 3-layer BP neural network topology, consisting of an input layer, a hidden layer, and an output layer. The input layer has z nodes, the hidden layer has H nodes, and the output layer has 1 node. Input v into the input layer. f Ultimately, there are corresponding expected and actual outputs BD; initialize the number of nodes, weights, and thresholds of each layer of the BP neural network;

[0064] S5.4 Optimize the BP neural network using the AMPSSO algorithm, including the following steps:

[0065] S5.4.1 Determine the particle dimension P in the AMPSSO algorithm v Its mathematical expression is:

[0066] ;

[0067] S5.4.2 Determine the particle fitness function and calculate the fitness of each particle. The mathematical expression for the particle fitness function is:

[0068] ;

[0069] Where: Y e (o) indicates the expected output of the o-th particle; Y a (o) indicates the actual output of the o-th particle;

[0070] S5.4.3 Repeat steps S3.7.3 to S3.7.5;

[0071] S5.4.4 Assign the results of the particle update to the weights and thresholds of the BP neural network; when the fitness value generated by the current particle swarm iteration update is less than the fitness value generated by the previous generation particle swarm iteration update, update the individual extreme value and the population extreme value; otherwise, proceed to the termination condition judgment.

[0072] S5.4.5 If the number of particle swarm update iterations meets the termination condition, the update stops and the BP neural network obtains the optimal weights and thresholds; otherwise, repeat steps S5.4.1 to S5.4.4 to continue updating the weights and thresholds of the BP neural network.

[0073] S5.4.6 Train the BP neural network with optimal weights and threshold values ​​to obtain the AMPSSO-BP neural network model.

[0074] The real-time measurement of topsoil bulk density mentioned in S6 is specifically as follows:

[0075] Soil samples with ring cutters were collected in real time using a ring cutter, and acoustic data W was collected according to the method described in step S2. a W aThe denoised signal W is obtained by inputting the AMPSSO-VMD model established in step S3. aj W aj Step S4 extracts the Mel-frequency inverse coefficients to obtain the feature vector. Finally, the feature vector is input into the AMPSSO-BP neural network model established in step S5 to obtain the actual bulk density value of the soil sample.

[0076] The population individual velocity update method described in step S3.7.4 is as follows:

[0077] If the population particle dimension is ≥5, then:

[0078] ;

[0079] Where: vc+1 p is the velocity of individual p in the updated population; vc p is the velocity of individual p in the current population; x best The optimal solution is dominated by individual P(te); xcp is the position of individual p in the current population; g best This represents the current globally optimal solution for the population; f<·> is the fitness operator; x P(te) Let ω1 be the position of individual P(te) in the previous population, where te takes the values ​​1, 2, 3, 4, and 5; ω1 is the weight coefficient of the velocity update variable; ω max1 Update the maximum inertia weight for velocity; ω min1 The minimum inertia weight is updated for velocity; c represents the c-th iteration; c max The maximum number of iterations is given; r1, r2, and r3 are three distinct random numbers between [0, 1].

[0080] The individual P(te) is determined as follows:

[0081] ;

[0082] Where: P(other) refers to all individuals in the population other than individual P(te);

[0083] If the population particle dimension is less than 5, then:

[0084] ;

[0085] Where: p βbest This represents the optimal solution for an individual in the current population.

[0086] The population individual location update method described in step S3.7.5 is as follows:

[0087] If the population particle dimension is ≥5, then:

[0088] ;

[0089] Where: vc+1 p is the position of individual p in the population after the first position update; r4 is a random number between [-1, 1]; r5 is a random number between [0, 1];

[0090] If the population particle dimension is less than 5, then:

[0091] ;

[0092] Where: xc othep is the position of any individual in the current population except individual p; r6 and r7 are random numbers between [0, 1]; r8 is a random number between [-1, 1].

[0093] The beneficial effects of this invention are as follows:

[0094] (1) In view of the problem that soil bulk density is a three-dimensional feature and cannot be measured by two-dimensional image information, this invention creates a method based on the penetration characteristics of step sound waves to determine soil bulk density. The method uses Mel reciprocal coefficient to extract the characteristics of sound waves penetrating the soil, which improves the timeliness of soil bulk density measurement.

[0095] (2) The VMD method was used to reduce the noise of the acoustic signal. The high-frequency noise of the data acquisition card sampling process and the environment was proposed. At the same time, a BP neural network was built to retrieve the soil bulk density value.

[0096] (3) To address the problem that swarm optimization algorithms are prone to getting trapped in local optima, this invention proposes an adaptive mutation particle swarm optimization starfish optimization (AMPSSO) algorithm. The adaptive mutation method is used to update the velocity of particles in the swarm, avoiding the 'convergence' effect of particles in the later stages of the algorithm and expanding the optimization range of the algorithm in the later stages. To address the problem of different optimization dimensions of the algorithm, optimization methods for different dimensions are proposed to improve the optimization range of the algorithm. The position and velocity of individuals in the swarm are updated using the starfish optimization algorithm. The starfish algorithm has better global optimization ability. Applying it to the particle swarm optimization algorithm can avoid the algorithm from getting trapped in local optima. The spiral hunting strategy in the whale algorithm is used to update the position of individuals, which improves the optimization accuracy of the algorithm.

[0097] (4) To address the problem of parameter optimization in different dimensions of VMD and BP neural networks, this invention uses the proposed AMPSSO algorithm to optimize the parameters and optimize VMD and BP neural networks. Attached Figure Description

[0098] Figure 1 This is a schematic diagram of a data acquisition system based on the penetration characteristics of acoustic waves.

[0099] Figure 2 This is a schematic diagram of the structure of data acquisition box A;

[0100] Figure 3 This is a structural schematic diagram of base component B;

[0101] Figure 4 This is a schematic diagram of the structure of the robotic arm device C.

[0102] Figure 5 This is a flowchart of the soil bulk density measurement model algorithm based on the acoustic wave penetration characteristic acquisition system;

[0103] Figure 6 The sound wave signal is without noise reduction processing;

[0104] Figure 7 The sound wave signal is processed by window sliding filter;

[0105] Figure 8 The sound wave signal is processed by AMPSSO-VMD;

[0106] Figure 9 To optimize the particle fitness curve of a BP neural network based on the PSO algorithm;

[0107] Figure 10 To optimize the particle fitness curve of a BP neural network based on the AMPSSO algorithm of this invention;

[0108] Among them: A. Acquisition box B. Base assembly C. Robotic arm device 1. Computer terminal 2. Data acquisition card 3. Digital power amplifier board 4. DC power supply 5. Top cover 6. Pull ring 7. Base 8. Sound insulation cotton I 9. Conduit I 10. Subwoofer 11. Sound insulation cotton II 12. Ring knife tube 13. Sound insulation cotton III 14. Sound insulation cotton IV 15. Pickup unit 16. Conduit II 17. Base 18. Multi-degree-of-freedom robotic arm body 19. End effector 20. Servo motor assembly. Detailed Implementation

[0109] The present invention will now be described in conjunction with the accompanying drawings.

[0110] like Figure 5 As shown, the present invention provides a method for measuring the bulk density of topsoil based on the sound wave penetration characteristics, comprising the following steps:

[0111] S1 is equipped with a sound wave penetration characteristic acquisition system;

[0112] S2 collects n sets of acoustic wave data W;

[0113] S3 uses the AMPSSO-VMD model to reduce noise in n sets of acoustic data W acoustic signals;

[0114] S4 extracts the Mel-frequency inverse coefficients of the n groups of noise-reduced Wj acoustic signals;

[0115] S5. Establish the AMPSSO-BP neural network model;

[0116] S6 Real-time measurement of soil bulk density in the topsoil layer.

[0117] like Figures 1 to 4 As shown, the acoustic wave penetration characteristic acquisition system described in S1 consists of an acquisition box A, a computer terminal 1, a data acquisition card 2, a digital power amplifier board 3, a DC power supply 4, and a robotic arm device C. The acquisition box A consists of a base assembly B and a top cover 5, with a pull ring 6 at the center of the top cover 5; the base assembly B consists of a base 7, ... The base 7 consists of sound insulation cotton I8, conduit I9, subwoofer 10, sound insulation cotton II11, ring cutter tube 12, sound insulation cotton III13, sound insulation cotton IV14, pickup 15, and conduit II16. The base 7 is a rectangular box structure. The ring cutter tube 12 is fixed to the center of the base plate; the cross-section of the ring cutter tube 12 is square. The pickup 15 is fixed to the rear of the rear plate of the ring cutter tube 12. The front end of the conduit II16 is connected to the rear end of the pickup 15, and the rear end of the conduit II16 is connected to the center of the rear plate of the base 7. The subwoofer 10 is fixed to the front of the front plate of the ring cutter tube 12. The rear end of the conduit I9 is ​​connected to the front end of the subwoofer 10, and the front end of the conduit I9 is ​​connected to the center of the front plate of the base 7. Sound insulation cotton I8 is placed in the space between the front plate of the ring cutter tube 12 and the front plate of the base 7. Sound insulation cotton IV14 is placed in the space between the rear plate of the ring cutter tube 12 and the rear plate of the base 7. The left plate of the ring cutter tube 12 is connected to the left plate of the base 7. Sound insulation cotton Ⅲ13 is placed in the space between the plates; sound insulation cotton Ⅱ11 is placed in the space between the right plate of the ring tube 12 and the right plate of the base 7; the digital power amplifier board 3 is electrically connected to the subwoofer 10, the DC power supply 4 and the computer terminal 1 respectively; the servo motor group 20 of the robotic arm device C is electrically connected to the computer terminal 1; the data acquisition card 2 is electrically connected to the microphone 15 and the computer terminal 1 respectively; the computer terminal 1 is equipped with a soil bulk density measurement model; the top cover 5 and the base 7 are made of acrylic material; the robotic arm device C consists of a base 17, a multi-degree-of-freedom robotic arm body 18, an end effector 19 and five servo motors of the servo motor group 20, wherein the base 17 is a cylindrical structure, the lower end of the multi-degree-of-freedom robotic arm body 18 is fixed to the base 17, and the upper end of the multi-degree-of-freedom robotic arm body 18 is connected to the end effector 19 through one servo motor of the servo motor group 20.

[0118] The acquisition of n sets of acoustic wave data W mentioned in S2 specifically refers to:

[0119] Using a ring cutter, n groups of soil samples with ring cutters are collected. The sound insulation cotton IV (14) and the top cover (5) are pulled out by the pull ring (6). Each soil sample is placed into the ring cutter tube (12) in sequence. Then, the sound insulation cotton IV (14) and the top cover (5) are put back into place by the pull ring (6). The computer terminal (1) emits a sinusoidal 200 Hz step sound wave signal. The digital power amplifier board (3) receives the sound wave signal and drives the subwoofer (10) to emit a step sound wave. The pickup (15) receives the sound wave data after passing through the soil sample. The sampling frequency of the data acquisition card (2) is set to 50000 Hz. The data acquisition card (2) collects the 0.1s sound wave data W of n groups of pickups (15) through the digital power amplifier board (3) and transmits it to the computer terminal (1).

[0120] S3 describes the use of the AMPSSO-VMD model to denoise n sets of acoustic data W acoustic signals, including the following steps:

[0121] S3.1 Set the mode function obtained after signal decomposition to I k I k Time Series I k The mathematical expression for (t) is:

[0122] ;

[0123] Among them: A k (t) is the amplitude, and A k (t)≥0;δ k (t) is a non-monotonic decreasing phase function;

[0124] S3.2 Initialize the VMD algorithm iteration count N, Lagrange multiplier λ, quadratic penalty factor α, and decomposition level K;

[0125] S3.3 The current number of VMD parameter iterations is m. Let m = m + 1, and the VMD algorithm will perform iterative calculations.

[0126] S3.4 Make the value of k continuously increase from 1 to K, I k ω k The update method is as follows:

[0127] ;

[0128] Where: Im+1 k is I k The updated value; IFTD<·> is the operator for Fourier transform and differentiation; ω k For I k The frequency center of (t); ωm+1 k is ω k The updated value; ω is the variable obtained by performing a Fourier transform on t;

[0129] S3.5 Updates λ, and the update method for λ is as follows:

[0130] ;

[0131] Where: λ m+1 The updated value of λ; IHT<·> is the inverse Hilbert transform;

[0132] S3.6 Repeat steps S3.3 to S3.5 until the termination condition is met:

[0133] ;

[0134] Where: ε is the convergence criterion, and ε > 0;

[0135] S3.7 The parameters α and K of VMD are determined using the AMPSSO algorithm, including the following steps:

[0136] S3.7.1 Determine the particle dimension in the AMPSSO algorithm. The VMD parameters to be optimized are α and K, therefore the particle dimension is 2:

[0137] S3.7.2 Determine the fitness function of the AMPSSO algorithm and calculate the fitness of each individual. The mathematical expression for the individual fitness function, Fitness, is as follows:

[0138] ;

[0139] Where: PE<·> is the permutation entropy operator; m l This is the l-th signal component after VMD;

[0140] S3.7.3 Initialize the individual positions and velocities in the AMPSSO algorithm;

[0141] S3.7.4 Perform population individual velocity updates;

[0142] S3.7.5 Update the location of individuals in the population;

[0143] S3.7.6 Assign the result of the second position update of the individual to the parameters α and K of the VMD algorithm. When the fitness value generated by the current population iteration update is less than the fitness value generated by the previous generation population iteration update, update the individual extreme value and the population extreme value; otherwise, proceed to the termination condition judgment.

[0144] S3.7.7 If the population update iteration number meets the termination condition, stop updating and obtain the optimal parameters α and K; otherwise, repeat steps S3.1 to S3.7.6 to continue updating parameters α and K.

[0145] S3.8 Repeat steps S3.1 to S3.7.7 for the n groups of signals in step S2 to obtain the optimal parameter α. best K best and α best K best Assigning values ​​to the VMD model yields the AMPSSO-VMD model, α best K best The calculation method is as follows:

[0146] ;

[0147] Where: round<·> is the rounding operator; α W(i) The optimal quadratic penalty factor parameters for the AMPSSO-VMD model of the i-th sound wave; K W(i) The optimal decomposition layer parameter for the AMPSSO-VMD model of the i-th sound wave;

[0148] S3.9 The AMPSSO-VMD model obtained in step S3.8 is used to decompose the W acoustic signal to obtain K signal components. The correlation coefficient between each component and W is calculated, and the signals with correlation coefficients higher than 0.6 are summed to obtain W. j .

[0149] S4 describes the extraction of Mel-frequency inverse coefficients from the n groups of noise-reduced Wj acoustic signals, which includes the following steps:

[0150] S4.1 to W j The mathematical expression for pre-emphasis processing of acoustic signals is:

[0151] ;

[0152] Among them: W pe (p) is W j The signal value after pre-emphasis at position p; W j (p), W j (p-1) represent W j Signal values ​​at positions p and p-1; η is a constant close to 1;

[0153] S4.2 to W pe (p) The acoustic signal is framed, with each frame being F in size. size For 1250 data points, frame shift F move With 500 data points, 9 sets of frame signals are obtained. The mathematical expression is:

[0154] ;

[0155] Where: F(q) is the q-th frame; zeros(1,250) is a 1×250 all-zero matrix;

[0156] S4.3 Windowing is applied to F(q) to avoid signal boundary effects; a Hamming window is used.

[0157] ;

[0158] Where: H(F(q)) y The signal value at point y after windowing;

[0159] S4.4 performs a Fast Fourier Transform on each windowed signal and then converts the amplitude spectrum into a power spectrum to enhance the spectral characteristics.

[0160] S4.5 A set of Mel filters is used to filter the power spectrum to improve the resolution in the low-frequency region, resulting in n sets of 40×9 arrays. The Mel filter bank contains 40 Mel filters, and the Mel frequency and frequency conversion in the Mel filters are as follows:

[0161] ;

[0162] Where: Mel is the Mel frequency, and f is the frequency;

[0163] S4.6 Convert the values ​​of the array in step S4.5 into decibel values, then perform discrete cosine transform on the decibel values, extract the first 13 rows of data, and obtain n sets of feature vectors containing 13×9=117 feature values.

[0164] The establishment of the AMPSSO-BP neural network model described in S5 includes the following steps:

[0165] S 5.1 Determine the true bulk density of n groups of soil samples with ring cutters; dry the soil samples in an oven at 105℃ for 12 hours. The mathematical expression for the true bulk density is:

[0166] ;

[0167] Where: m w The mass of the soil sample with the ring cutter before drying, m d V represents the mass of the dried soil sample with the ring cutter. r For the volume of the ring cutter;

[0168] S5.2 Select the feature values ​​in the feature vector from step S4.6 that have a high correlation with BD, calculate the Pearson coefficient between each feature value in the feature vector from step S4.6 and BD, and select z feature values ​​with a Pearson coefficient greater than 0.5 to form n new feature vectors v. f ;

[0169] S5.3 Establish a 3-layer BP neural network topology, consisting of an input layer, a hidden layer, and an output layer. The input layer has z nodes, the hidden layer has H nodes, and the output layer has 1 node. Input v into the input layer. f Ultimately, there are corresponding expected and actual outputs BD; initialize the number of nodes, weights, and thresholds of each layer of the BP neural network;

[0170] S5.4 Optimize the BP neural network using the AMPSSO algorithm, including the following steps:

[0171] S5.4.1 Determine the particle dimension P in the AMPSSO algorithm v Its mathematical expression is:

[0172] ;

[0173] S5.4.2 Determine the particle fitness function and calculate the fitness of each particle. The mathematical expression for the particle fitness function is:

[0174] ;

[0175] Where: Y e (o) indicates the expected output of the o-th particle; Y a (o) indicates the actual output of the o-th particle;

[0176] S5.4.3 Repeat steps S3.7.3 to S3.7.5;

[0177] S5.4.4 Assign the results of the particle update to the weights and thresholds of the BP neural network; when the fitness value generated by the current particle swarm iteration update is less than the fitness value generated by the previous generation particle swarm iteration update, update the individual extreme value and the population extreme value; otherwise, proceed to the termination condition judgment.

[0178] S5.4.5 If the number of particle swarm update iterations meets the termination condition, the update stops and the BP neural network obtains the optimal weights and thresholds; otherwise, repeat steps S5.4.1 to S5.4.4 to continue updating the weights and thresholds of the BP neural network.

[0179] S5.4.6 Train the BP neural network with optimal weights and threshold values ​​to obtain the AMPSSO-BP neural network model.

[0180] The real-time measurement of topsoil bulk density mentioned in S6 is specifically as follows:

[0181] Soil samples with ring cutters were collected in real time using a ring cutter, and acoustic data W was collected according to the method described in step S2. a W aThe denoised signal W is obtained by inputting the AMPSSO-VMD model established in step S3. aj W aj Step S4 extracts the Mel-frequency inverse coefficients to obtain the feature vector. Finally, the feature vector is input into the AMPSSO-BP neural network model established in step S5 to obtain the actual bulk density value of the soil sample.

[0182] The population individual velocity update method described in step S3.7.4 is as follows:

[0183] If the population particle dimension is ≥5, then:

[0184] ;

[0185] Where: vc+1 p is the velocity of individual p in the updated population; vc p is the velocity of individual p in the current population; x best The optimal solution is dominated by individual P(te); xcp is the position of individual p in the current population; g best This represents the current globally optimal solution for the population; f<·> is the fitness operator; x P(te) Let ω1 be the position of individual P(te) in the previous population, where te takes the values ​​1, 2, 3, 4, and 5; ω1 is the weight coefficient of the velocity update variable; ω max1 Update the maximum inertia weight for velocity; ω min1 The minimum inertia weight is updated for velocity; c represents the c-th iteration; c max The maximum number of iterations is given; r1, r2, and r3 are three distinct random numbers between [0, 1].

[0186] The individual P(te) is determined as follows:

[0187] ;

[0188] Where: P(other) refers to all individuals in the population other than individual P(te);

[0189] If the population particle dimension is less than 5, then:

[0190] ;

[0191] Where: p βbest This represents the optimal solution for an individual in the current population.

[0192] The population individual location update method described in step S3.7.5 is as follows:

[0193] If the population particle dimension is ≥5, then:

[0194] ;

[0195] Where: vc+1 p is the position of individual p in the population after the first position update; r4 is a random number between [-1, 1]; r5 is a random number between [0, 1];

[0196] If the population particle dimension is less than 5, then:

[0197] ;

[0198] Where: xc othep is the position of any individual in the current population except individual p; r6 and r7 are random numbers between [0, 1]; r8 is a random number between [-1, 1].

[0199] Traditional methods typically employ window smoothing filtering algorithms for signal denoising. Compared to these methods, the AMPSSO-VMD algorithm of this invention exhibits significantly superior denoising performance. Furthermore, while traditional methods use the PSO algorithm for parameter optimization, the AMPSSO algorithm of this invention demonstrates superior performance in terms of particle fitness. Specific experimental data and analysis are as follows: Regarding signal denoising, by… Figures 6 to 8 As shown, traditional window smoothing filtering algorithms can reduce noise in the original signal, but the AMPSSO-VMD denoised signal of this invention has fewer 'glitch' and a smoother waveform; in terms of algorithm performance, it is superior to traditional methods. Figure 9 and Figure 10 As shown, with the algorithm terminating at 300 iterations, the final fitness of the particles in the BP neural network optimized by the PSO algorithm is 0.1407, while the final fitness of the particles in the BP neural network optimized by the AMPSSO algorithm of this invention is 0.1303. Therefore, the soil bulk density measurement method of this invention is more effective.

Claims

1. A method for measuring the bulk density of topsoil based on the penetration characteristics of sound waves, characterized in that, Includes the following steps: S1 is equipped with a sound wave penetration characteristic acquisition system; S2 collects n sets of acoustic wave data W; S3 uses the AMPSSO-VMD model to reduce noise in n sets of acoustic data W acoustic signals; S4 extracts the Mel-frequency inverse coefficients of the n groups of noise-reduced Wj acoustic signals; S5. Establish the AMPSSO-BP neural network model; S6 Real-time measurement of soil bulk density in the topsoil layer.

2. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 1, characterized in that, The acoustic wave penetration characteristic acquisition system described in S1 consists of an acquisition box (A), a computer terminal (1), a data acquisition card (2), a digital power amplifier board (3), a DC power supply (4), and a robotic arm device (C). The acquisition box (A) consists of a base assembly (B) and a top cover (5). The top cover (5) has a pull ring (6) at its center. The base assembly (B) consists of a base (7), sound insulation cotton I (8), a wire conduit I (9), a subwoofer (10), sound insulation cotton II (11), a ring cutter tube (12), sound insulation cotton III (13), sound insulation cotton IV (14), a microphone (15), and a wire conduit II (16). The base (7) is a rectangular box structure. The ring cutter tube (12) is fixed to the center of the base plate. The cross-section of the ring cutter tube (12) is square. The microphone (15) is fixed to the back of the ring cutter tube (12). The front end of the wire conduit II (16) is connected to the rear end of the microphone (15). The rear end of conduit II (16) is connected to the center of the rear plate of the base (7); the front of the front plate of the ring tube (12) is fixed to the subwoofer (10), the rear end of conduit I (9) is connected to the front end of the subwoofer (10), and the front end of conduit I (9) is connected to the center of the front plate of the base (7); sound insulation cotton I (8) is placed in the space between the front plate of the ring tube (12) and the front plate of the base (7); sound insulation cotton IV (14) is placed in the space between the rear plate of the ring tube (12) and the rear plate of the base (7); the left plate of the ring tube (12) is connected to the base (7) Sound insulation cotton III (13) is placed in the space between the left and right plates; sound insulation cotton II (11) is placed in the space between the right plate of the ring knife tube (12) and the right plate of the base (7); the digital power amplifier board (3) is electrically connected to the subwoofer (10), the DC power supply (4) and the computer terminal (1) respectively; the servo motor group (20) of the robotic arm device (C) is electrically connected to the computer terminal (1); the data acquisition card (2) is electrically connected to the microphone (15) and the computer terminal (1) respectively; the computer terminal (1) is equipped with soil The bulk density measurement model; the top cover (5) and the base (7) are made of acrylic material; the robotic arm device (C) consists of a base (17), a multi-degree-of-freedom robotic arm body (18), an end effector (19) and five servo motors of the servo motor group (20). The base (17) is a cylindrical structure. The lower end of the multi-degree-of-freedom robotic arm body (18) is fixed to the base (17). The upper end of the multi-degree-of-freedom robotic arm body (18) is connected to the end effector (19) through one servo motor of the servo motor group (20).

3. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 1, characterized in that, The acquisition of n sets of acoustic wave data W mentioned in S2 specifically refers to: Using a ring cutter, n groups of soil samples with ring cutters are collected. The sound insulation cotton IV (14) and the top cover (5) are pulled out by the pull ring (6). Each soil sample is placed into the ring cutter tube (12) in sequence. Then, the sound insulation cotton IV (14) and the top cover (5) are put back into place by the pull ring (6). The computer terminal (1) emits a sinusoidal 200 Hz step sound wave signal. The digital power amplifier board (3) receives the sound wave signal and drives the subwoofer (10) to emit a step sound wave. The pickup (15) receives the sound wave data after passing through the soil sample. The sampling frequency of the data acquisition card (2) is set to 50000 Hz. The data acquisition card (2) collects the 0.1s sound wave data W of n groups of pickups (15) through the digital power amplifier board (3) and transmits it to the computer terminal (1).

4. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 1, characterized in that, S3 describes the use of the AMPSSO-VMD model to denoise n sets of acoustic data W acoustic signals, including the following steps: S3.1 Set the mode function obtained after signal decomposition to I k I k Time Series I k The mathematical expression for (t) is: ； Among them: A k (t) is the amplitude, and A k (t)≥0;δ k (t) is a non-monotonic decreasing phase function; S3.2 Initialize the VMD algorithm iteration count N, Lagrange multiplier λ, quadratic penalty factor α, and decomposition level K; S3.3 The current number of VMD parameter iterations is m. Let m = m + 1, and the VMD algorithm will perform iterative calculations. S3.4 Make the value of k continuously increase from 1 to K, I k ω k The update method is as follows: ； Where: Im+1 k is I k The updated value; IFTD<·> is the operator for Fourier transform and differentiation; ω k For I k The frequency center of (t); ωm+1 k is ω k The updated value; ω is the variable obtained by performing a Fourier transform on t; S3.5 Updates λ, and the update method for λ is as follows: ； Where: λ m+1 The updated value of λ; IHT<·> is the inverse Hilbert transform; S3.6 Repeat steps S3.3 to S3.5 until the termination condition is met: ； Where: ε is the convergence criterion, and ε > 0; S3.7 The parameters α and K of VMD are determined using the AMPSSO algorithm, including the following steps: S3.7.1 Determine the particle dimension in the AMPSSO algorithm. The VMD parameters to be optimized are α and K, therefore the particle dimension is 2: S3.7.2 Determine the fitness function of the AMPSSO algorithm and calculate the fitness of each individual. The mathematical expression for the individual fitness function, Fitness, is as follows: ； Where: PE<·> is the permutation entropy operator; m l This is the l-th signal component after VMD; S3.7.3 Initialize the individual positions and velocities in the AMPSSO algorithm; S3.7.4 Perform population individual velocity updates; S3.7.5 Update the location of individuals in the population; S3.7.6 Assign the result of the second position update of the individual to the parameters α and K of the VMD algorithm. When the fitness value generated by the current population iteration update is less than the fitness value generated by the previous generation population iteration update, update the individual extreme value and the population extreme value; otherwise, proceed to the termination condition judgment. S3.7.7 If the population update iteration number meets the termination condition, stop updating and obtain the optimal parameters α and K; otherwise, repeat steps S3.1 to S3.7.6 to continue updating parameters α and K. S3.8 Repeat steps S3.1 to S3.7.7 for the n groups of signals in step S2 to obtain the optimal parameter α. best K best and α best K best Assigning values ​​to the VMD model yields the AMPSSO-VMD model, α best K best The calculation method is as follows: ； Where: round<·> is the rounding operator; α W(i) The optimal quadratic penalty factor parameters for the AMPSSO-VMD model of the i-th sound wave; K W(i) The optimal decomposition layer parameter for the AMPSSO-VMD model of the i-th sound wave; S3.9 The AMPSSO-VMD model obtained in step S3.8 is used to decompose the W acoustic signal to obtain K signal components. The correlation coefficient between each component and W is calculated, and the signals with correlation coefficients higher than 0.6 are summed to obtain W. j .

5. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 1, characterized in that, The extraction of Mel-frequency inverse coefficients from the n groups of noise-reduced Wj acoustic signals described in S4 includes the following steps: S4.1 to W j The mathematical expression for pre-emphasis processing of acoustic signals is: ； Among them: W pe (p) is W j The signal value after pre-emphasis at position p; W j (p), W j (p-1) represent W j Signal values ​​at positions p and p-1; η is a constant close to 1; S4.2 to W pe (p) The acoustic signal is framed, with each frame being F in size. size For 1250 data points, frame shift F move With 500 data points, 9 sets of frame signals are obtained. The mathematical expression is: ； Where: F(q) is the q-th frame; zeros(1,250) is a 1×250 all-zero matrix; S4.3 Windowing is applied to F(q) to avoid signal boundary effects; a Hamming window is used. ； Where: H(F(q)) y The signal value at point y after windowing; S4.4 performs a Fast Fourier Transform on each windowed signal and then converts the amplitude spectrum into a power spectrum to enhance the spectral characteristics. S4.5 A set of Mel filters is used to filter the power spectrum to improve the resolution in the low-frequency region, resulting in n sets of 40×9 arrays. The Mel filter bank contains 40 Mel filters, and the Mel frequency and frequency conversion in the Mel filters are as follows: ； Where: Mel is the Mel frequency, and f is the frequency; S4.6 Convert the values ​​of the array in step S4.5 into decibel values, then perform discrete cosine transform on the decibel values, extract the first 13 rows of data, and obtain n sets of feature vectors containing 13×9=117 feature values.

6. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claims 1 and 4, characterized in that, The establishment of the AMPSSO-BP neural network model described in S5 includes the following steps: S 5.1 Determine the true bulk density of n groups of soil samples with ring cutters; dry the soil samples in an oven at 105℃ for 12 hours. The mathematical expression for the true bulk density is: ； Where: m w The mass of the soil sample with the ring cutter before drying, m d V represents the mass of the dried soil sample with the ring cutter. r For the volume of the ring cutter; S5.2 Select the feature values ​​in the feature vector from step S4.6 that have a high correlation with BD, calculate the Pearson coefficient between each feature value in the feature vector from step S4.6 and BD, and select z feature values ​​with a Pearson coefficient greater than 0.5 to form n new feature vectors v. f ; S5.3 Establish a 3-layer BP neural network topology, consisting of an input layer, a hidden layer, and an output layer. The input layer has z nodes, the hidden layer has H nodes, and the output layer has 1 node. Input v into the input layer. f Ultimately, there are corresponding expected and actual outputs BD; initialize the number of nodes, weights, and thresholds of each layer of the BP neural network; S5.4 Optimize the BP neural network using the AMPSSO algorithm, including the following steps: S5.4.1 Determine the particle dimension P in the AMPSSO algorithm v Its mathematical expression is: ； S5.4.2 Determine the particle fitness function and calculate the fitness of each particle. The mathematical expression for the particle fitness function is: ； Where: Y e (o) indicates the expected output of the o-th particle; Y a (o) indicates the actual output of the o-th particle; S5.4.3 Repeat steps S3.7.3 to S3.7.5; S5.4.4 Assign the results of the particle update to the weights and thresholds of the BP neural network; when the fitness value generated by the current particle swarm iteration update is less than the fitness value generated by the previous generation particle swarm iteration update, update the individual extreme value and the population extreme value; otherwise, proceed to the termination condition judgment. S5.4.5 If the number of particle swarm update iterations meets the termination condition, the update stops and the BP neural network obtains the optimal weights and thresholds; otherwise, repeat steps S5.4.1 to S5.4.4 to continue updating the weights and thresholds of the BP neural network. S5.4.6 Train the BP neural network with optimal weights and threshold values ​​to obtain the AMPSSO-BP neural network model.

7. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 1, characterized in that, The real-time measurement of topsoil bulk density mentioned in S6 is specifically as follows: Soil samples with ring cutters were collected in real time using a ring cutter, and acoustic data W was collected according to the method described in step S2. a W a The denoised signal W is obtained by inputting the AMPSSO-VMD model established in step S3. aj W aj Step S4 extracts the Mel-frequency inverse coefficients to obtain the feature vector. Finally, the feature vector is input into the AMPSSO-BP neural network model established in step S5 to obtain the actual bulk density value of the soil sample.

8. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 4, characterized in that, The population individual velocity update method described in step S3.7.4 is as follows: If the population particle dimension is ≥5, then: ； Where: vc+1 p is the velocity of individual p in the updated population; vc p is the velocity of individual p in the current population; x best The optimal solution is dominated by individual P(te); xcp is the position of individual p in the current population; g best This represents the current globally optimal solution for the population; f<·> is the fitness operator; x P(te) Let ω1 be the position of individual P(te) in the previous population, where te takes the values ​​1, 2, 3, 4, and 5; ω1 is the weight coefficient of the velocity update variable; ω max1 Update the maximum inertia weight for velocity; ω min1 The minimum inertia weight is updated for velocity; c represents the c-th iteration; c max The maximum number of iterations is given; r1, r2, and r3 are three distinct random numbers between [0, 1]. The individual P(te) is determined as follows: ； Where: P(other) refers to all individuals in the population other than individual P(te); If the population particle dimension is less than 5, then: ； Where: p βbest This represents the optimal solution for an individual in the current population.

9. The method for measuring the bulk density of topsoil based on the sound wave penetration characteristics according to claim 4, characterized in that, The population individual location update method described in step S3.7.5 is as follows: If the population particle dimension is ≥5, then: ； Where: vc+1 p is the position of individual p in the population after the first position update; r4 is a random number between [-1, 1]; r5 is a random number between [0, 1]; If the population particle dimension is less than 5, then: ； Where: xc othep is the position of any individual in the current population except individual p; r6 and r7 are random numbers between [0, 1]; r8 is a random number between [-1, 1].

Images