Calibration of Acoustic-Soil Discrete Element Model and Analysis of Influencing Factors on Accuracy

: To obtain accurate soil parameters and improve the accuracy of the acoustic-soil discrete element simulation model, we studied the model’s parameter calibration. The simulation test was carried out using the measured acoustic velocity and dominant frequency as the response value (index). Firstly, the Plackett — Burman test scheme was used to obtain the sensitivity ranking of soil parameters to the dominant frequency and velocity of the acoustic wave. The parameters that significantly affect the acoustic wave were obtained: Shear modulus, Poisson’s ratio, and coefficient of restitution. Then the Box — Behnken test scheme was used to establish the regression relationship between the dominant frequency and the velocity of the sound wave and Shear modulus, Poisson’s ratio, and Coefficient of restitution. The results shows that the indexes that researchers focus on are different in different scenarios, and the sensitivity of soil parameters to different indicators is different, which results in different soil parameter values after calibration. This study analyzed the main factors affecting the accuracy of the acoustic-discrete element model in constructing the model, provided a method for improving the construction accuracy of the acoustic wave — soil discrete element model and provided a reference for the construction of discrete elements models in other fields.


Introduction
There are more and more studies on the detection of farmland soil properties by acoustic waves, such as soil water content detection [1][2][3], soil compaction detection [4,5], and seed spatial position detection in soil [6], etc.To study the process of acoustic wave propagation in soil and clarify the propagation mechanism of the acoustic wave in discrete bodies to better guide practical applications, some scholars have introduced the discrete element method [7,8].This method takes a single soil particle as an independent unit for calculation, can analyze information such as the force and motion of a single particle, and can visualize the acoustic wave propagation process.Therefore, it can promote the study of acoustic wave propagation in soil.
Parameter setting is a prerequisite for discrete element simulation, and the accuracy of the input parameters affects the error of simulation results, which is the key to whether the simulation can genuinely reflect the actual situation [9][10][11].Parameter calibration can make soil particle simulation parameters conform to actual physical parameters, thus reflecting soil properties.Therefore, parameter calibration is critical in establishing an accurate acoustic-soil discrete element model.Due to the complexity of the soil structure, some constitutive parameters and contact parameters such as the shear modulus of soil particles, the rolling friction parameters between particles, and the restitution coefficient of collision between particles cannot be directly obtained in practical experiments [12,13].At the same time, there is no detailed theory to establish the quantitative relationship between soil particle parameters and macroscopic mechanical response.Therefore, these parameters are currently obtained through parameter calibration methods [14][15][16][17][18][19].That is, by continuously changing the soil particle parameters to obtain a set of simulation parameter combinations that can make the simulation experiment and the actual experiment macroscopic phenomenon consistent.
At present, there are two main methods for calibrating simulation parameters.One is the traditional calibration method with repose angle or sliding friction angle as a macro index.Nakashima analyzed the effect of gravity on different parameters (such as friction coefficient and rolling friction coefficient) by comparing the angle of repose of simulations and experiments to study the terrain-mechanical interaction under the condition of gravity on the lunar surface [20].Marigo used the angle of repose as an index to calibrate cylindrical pellets [16].The obtained parameters were used for the simulation of a rotating drum, and the simulation and experimental results were compared.While qualitative comparisons were good in most cases, there were some errors in quantitative comparison results.To get the soil particle parameters of the soil DEM (discrete element method) simulation model, Song Zhanhua calibrated the soil particle contact parameters with the soil repose angle and sliding friction angle as indexes [21].The above calibrations are all calibrated by the traditional repose or sliding friction angle, and the same calibration indexes (angle of repose, interactive friction angle) are used for different application scenarios.Due to the significant differences in using scenarios, researchers pay attention to different indexes when establishing models, which may lead to that the models can not satisfy different scenarios and purposes.Therefore, it is necessary to search for the characteristic values that researchers pay special attention to in specific scenarios as calibration indexes according to the research purpose.The second is the method of calibration based on the actual scenario, calibrated by macro-indexes that can represent the main features of the material in the scenario.Coetzee determined the particles' internal friction angle and stiffness values with a combination of shear and compression test results [17].The calibration process is validated by modelling silo discharge and bucket filling.To develop a 3D numerical simulation model of burden distribution for the pre-reduction shaft furnace of COREX 3000, Li, Qiang determined the parameters of DEM (Discrete element method) by quantitatively comparing the angle of repose and discharge time of the hopper of simulation and experiment [19].Furthermore, they verified the parameters' feasibility.Shen Haohan calibrated the rock's elastic modulus, Poisson's ratio, and tensile strength, obtained various parameters of the particles, constructed a discrete element model of the rock, and verified the model by experiments [22].In the above calibration experiments, according to different application scenarios, in addition to selecting the angle of repose index, other macro indexes that can represent the characteristics of the material and the critical indexes that researchers pay attention to during the experiment are also selected for calibration.This method can help researchers find the main particle parameters that affect macroscopic indexes, narrow the range of input parameters, and improve calibration accuracy.Therefore, the selected macro indexes for the calibration should be different for different application scenarios and research purposes.
In addition, when a single macro index is used for calibration, different particle parameter sets can be obtained; that is, one macro index can obtain many different particle parameter combinations.Each parameter combination satisfies the macroscopic phenomenon, so it is impossible to distinguish which parameter combinations are reliable.When the number of macroscopic indexes is increased to calibrate the particle parameters, the union of each parameter combination can be obtained to narrow the range of parameters to obtain a more accurate calibration effect.Therefore, the indexes that can characterize the material characteristics should be found according to the actual application scenario when calibrating.At the same time, the number of indexes should be increased as much as possible to improve the accuracy of parameter calibration.
According to the acoustic theory, the main characteristics of the wave are the velocity transmission and the dominant frequency, which are also two critical indexes in the acoustic analysis.Moreover, these two macroscopic indexes can be accurately measured in simulation and experiments.They are the macroscopic manifestation of the interaction between the soil medium and the acoustic wave, which can reflect the change of the acoustic wave in the medium to a certain extent.Therefore, in this research, the parameters of soil particles are fused and calibrated using the velocity and dominant frequency of the acoustic as indexes.In this research, the sandy soil is used as the calibration object, the acoustic velocity, and the dominant frequency of the received signal are used as the response, and the soil particle parameters that significantly affect the simulation results are screened through the Plackett-Burman test.Then, the significant parameters were calibrated by the Box-Behnken design test, and the optimal parameter combination was obtained, which was compared with the measured value to provide reliable soil particle parameters for building an accurate model of acoustic-soil interaction.
In addition, in the previous study, we have established a physical model of soundsoil propagation, so the introduction of the process of model establishment is relatively simple.

Actual Scene Experiment
To obtain the calibration indexes, a self-made acoustic wave test bench (Figure 1) [8] was used to get the value of acoustic velocity and dominant frequency.It includes the RGM-4005 universal material testing machine, ZBL-U5200 nonmetal ultrasonic tester, and self-made soil acoustic wave testing device.RGM-4005 Universal Material Tester is manufactured by Rigel Instruments Co., Ltd. and consists of a computer and a test machine host.The maximum load is 5 kN, and the relative error of the test force is within <1%.ZBL-U5200 non-metallic ultrasonic tester is manufactured by Beijing ZBL Science & Technology Co., Ltd.(Beijing, China).Ultrasonic time accuracy of this detector is ±0.025 μs, the sampling interval is 0.05 μs.The self-made soil acoustic wave testing device mainly included the box, excitation and receiving transducer, etc.The height, length, and thickness of the box were 130, 100 and 40 mm, respectively, which were the same as in the simulation.The soil was dried by drying oven and the mass moisture content is 0%.After screening, soil particles with the particle size range of 1~2 mm were taken for the experiment.In this experiment the soil compression ratio was 12%.When starting the experiment, we filled the box with soil, and the size and shape of the box were the same as those of the discrete element model box.Then the universal material machine compressed the soil through the clamping block above the box until the compression ratio reached 12%, and the clamping block remained stationary.The acoustic experiment was carried out in a pressurized state.The transducer frequency was about 40 kHz.The excitation transducer excited the acoustic wave.The receiving transducer received the acoustic wave and collected the received signal waveform.Then, we calculated the wave velocity c and the dominant frequency F of the received signal.The measurement result was that the acoustic velocity was 507 m/s, and the dominant frequency was 18 kHz.The results were used as the target value (index) for discrete element calibration.

DEM Calculation Method
An ideal model is used to reduce the influence of complex soil environments, such as different soil particle sizes and water content, on acoustic wave attenuation.The Hertz-Mindlin (no slip) model [23,24] was chosen for this study, which is mainly used for discrete-body simulations.Normal damping, tangential damping, normal stress, tangential stress, and friction forces and moment among soil particles are calculated as follows.
Normal force Fn: (1 ) ( 1) where where Sn represents the normal stiffness; m* represents the equivalent mass; vn rel represents the normal component of the relative velocity; β represents the damping coefficient.
The tangential force Ft depends on the tangential stiffness St and the tangential overlap δt.
where vt rel represents the relative velocity; G* represents the equivalent tangential modulus.
where Ri represents the distance of the contact point from the center of mass, ωi represents the unit angular velocity vector of the object at the contact point and μr represents the coefficient of rolling friction.

Discrete Element Model Building
Since this model has been established in previous research [8,25], the critical parameters and settings are briefly introduced.The DEM model is established according to the experimental bench, which only retains the soil, box, pressure plate, excitation transducer and receiving transducer related to acoustic wave propagation so as to obtain the propagation law of acoustic waves in soil.As shown in Figure 2a,b, the box is the same as the experiment.The soil is compressed by the pressure plate and the h1 is the downward distance of the pressure plate.The excitation transducer transmits acoustic signals, and the receiving transducer receives acoustic signals.

Simulation Parameter Settings
The particle size is 2 mm.The fixed time step is 1.0 × 10 −8 s, 0.812% of the contact time.The grid cell sizes is 2 mm and the gravity is 9.8 m/s 2 .The number of soils is 600,000.The pressing plate compressed the soil until the compression ratio reached 12% (Equation ( 12)).The excitation transducer sends a single period sine wave to the soil in the y-axis direction (Figure 2b) [26].The excitation amplitude A1 is set to 10 μm [27], and the excitation frequency f is 40 kHz, the same as the experiment.The relationship between the excitation transducer amplitude, frequency, phase angle, and period in the DEM model is shown in Equations ( 13)- (15).
where h1 represents the soil compression height, mm; α represents the soil compression ratio; h represents the height of soil, mm.
1/ Tf = (15) where A1 represents the maximum amplitude (displacement) of the transducer motion, um; y represents the instantaneous displacement, um; t represents the time, s; f represents the vibration frequency, Hz; T represents the motion period, s; θ represents the angle of phase, rad.
The aim of this research is to calibrate the parameters in soil particles that significantly affect the acoustic wave characteristics.The function of the transducer is to apply an external excitation to the soil, thereby sending out acoustic wave signals.The parameters of the transducer and the contact parameters between the transducer and the soil particles will not affect the propagation law of the acoustic wave in the soil.Furthermore, the box plays the role of containing the soil and does not affect the acoustic wave propagation process.For the above reasons, the parameters of the box and transducer and the contact parameters between them and soil particles are not focused.Determining parameter values mainly refers to previous literature [27][28][29].The soil particle parameters were determined by calibration as factors of the experiment and will be described in detail in Section 3. The material of the pressure plate and the box is steel, and the transducer's material is piezoelectric ceramics (pzt).Physical parameters and contact parameters are listed in Tables 1 and 2, respectively.

Index Measurement
(1) Acoustic velocity The acoustic velocity is one of the main characteristics of the acoustic wave timedomain signal, which can reflect the physical properties of the transmission medium so that the acoustic velocity can be used to calibrate parameters.Figure 3a is the time-domain curve of the acoustic pressure P of the receiving transducer in the discrete element model.
Where d represents the transmission distance, in m; Δt represents the transmission time, in s. (

2) Dominant frequency
The dominant frequency of an acoustic wave can reflect a large number of physical and mechanical properties of the medium.Therefore, the dominant frequency is taken as the calibration index.Figure 3 is the Fourier transform of the first wave of the received signal, and the calculation method is based on formula (17).The dominant frequency value F is the frequency corresponding to the maximum amplitude in the spectrum curve, as shown in Figure 3b.Similarly, the experimental calculation method is the same.

Test Scheme
The constitutive parameters of soil particles reflect soil's inherent properties, including particle shear modulus, Poisson's ratio, and density.Since acoustic waves propagate in elastic media, the constitutive parameters of soil particles are essential parameters that affect the propagation characteristics.The contact parameters include the coefficient of restitution, the coefficient of static friction and the coefficient of kinetic friction of the particles.The contact parameters are related to the surface smoothness of the particles and the particles shape.In order to simplify the model, This paper replaces the soil particle shape with a spherical shape, so the contact parameters must be calibrated.Therefore, the experimental factors selected in this paper are particle density, Poisson's ratio, shear modulus, coefficient of restitution, coefficient of static friction and coefficient of kinetic friction.Not all factors significantly affect acoustic wave propagation, so this paper uses the Plackett-Burman design screening test to analyze the sensitivity of the constitutive particle parameters and contact parameters and obtains the factors that significantly affect the soil acoustic wave propagation test.Subsequently, in order to obtain the best combination of soil parameters, the Box-Behnken test was carried out.

Plackett-Burman Test Scheme
The simulation test was carried out with the Poisson's ratio, shear modulus, density, restitution coefficient, static friction coefficient and rolling friction coefficient of soil particles as the factors, and the acoustic velocity v and the dominant frequency F as the indexes.There are a total of 6 factors, each with a high and a low level.According to the literature [30][31][32], the range of the parameters value of soil particles is roughly obtained.The factors and levels are listed in Table 3.

Box-Behnken Test Scheme
According to the results of the Plackett-Burman test, it can be seen that X3, X5 and X6 are insignificant factors, so we take the median value of soil parameters, which are 2650 kg/m 3 , 0.5, and 0.25, respectively.The Box-Behnken design test was applied to the significant factors to conduct response surface analysis and find the optimal solution.X1 (Poisson's ratio), X2 (Shear modulus) and X4 (Coefficient of restitution) were used as test factors.The values of specific factors are shown in Table 4.The velocity and dominant frequency of acoustic waves are the test indexes.

Plackett-Burman Test Results and Sensitivity Analysis
Table 5 shows the simulation test results of acoustic velocity and dominant frequency under each factor level and analyzed by Design-Expert 8.0.6 software according to the results.
Table 6 shows the analysis of the variance of the Plackett-Burman experiment.It can be seen from the p-value that within the selected parameter value range, Poisson's ratio (X1), Shear modulus (X2), and Coefficient of restitution (X4) have a significant impact on the acoustic velocity and dominant frequency.Table 7 is the order of the contribution of each factor to the two indexes of the acoustic velocity and the dominant frequency of the received signal.It can be seen from the table that the order of the contribution of each factor to the acoustic velocity from large to small is X2 > X1 > X4 > X5 > X6 > X3; the order of influence on the dominant frequency of the acoustic wave is X2 > X1 > X4 > X6 > X3 > X5.In order to make the model optimization simple and feasible, Poisson's ratio (X1), Shear modulus (X2), and Coefficient of restitution (X4), which have a significant impact on both indexes, are used for parameter calibration.(1) Significance analysis of acoustic wave dominant frequency Y1 Table 9 shows that the regression model of the dominant frequency Y1 is p < 0.0001 (highly significant), and the lack of fit term is p = 0.7221 > 0.05 (not significant), indicating that the regression model is reliable and can be used for the optimization of parameter calibration.In addition, the coefficient of determination of the model R 2 = 0.9987, indicating that the fitting degree is good; the correction coefficient of determination Adj.R 2 = 0.9970, indicating that the predicted value has a high correlation with the actual value.Therefore, the regression model can be used to analyze and predict the dominant frequency Y1 of the acoustic wave.In the regression model, the primary items x1, x2, x4, the quadratic items x2 2 and the interaction items x1x2, x2x4 have a very significant effect on the dominant frequency of the acoustic wave (p < 0.01), while the quadratic items x1 2 , x4 2 and the interaction items x1x3 are all no significant effect (p > 0.05).Therefore, remove x1 2 , x4 2 and the interaction term x1x4, and obtain the regression equation between various factors and indexes (Equation ( 18)).(2) Significance analysis of acoustic velocity Y2 Table 9 shows that the regression model of the acoustic velocity Y2 is p < 0.0001 (highly significant), and the lack of fit term p = 0.0957 > 0.05 (not significant), indicating that the regression model is reliable and can be used for the optimization of parameter calibration.In addition, the coefficient of determination of the model R 2 = 0.9980, indicating a good degree of fitting; the correction coefficient of determination Adj.R 2 = 0.9954, indicating that the predicted value has a high correlation with the actual value.Therefore, the regression model can be used to analyze and predict the acoustic velocity Y2.In the regression model, the primary items x1, x2, x4, and the quadratic items x1 2 and x2 2 have extremely significant effects on the acoustic velocity (p < 0.01); the quadratic items x4 2 and interaction items x2x4 are significant levels (0.05 < p< 0.01); the interaction items x1x2 and x1x4 had no significant effect (p > 0.05).Therefore, the interaction terms x1x2 and x1x4 are eliminated, and the regression equation between various factors and indexes is obtained (Equation ( 19)).The measured parameter value of the actual material will fluctuate within a certain range, but the fluctuation amount is small.When performing parameter calibration, the smaller the variable range of material parameters, the more successful the calibration.In order to compare the influence of the number of indexes on the calibration effect, we compared the desirability of the single-index and dual-index calibration results, respectively.The higher the desirability of parameter calibration prediction, the better the prediction effect and the more consistent with the actual material parameter value.We will try to select the parameter combination with higher desirability for calibration.The smaller the range of parameters that can be selected, the better it can reflect the essential properties of the actual material.Figures 4-6 are the contour plots of the desirability of the predicted values of the regression equation after the orthogonal simulation test.
Figure 4 is the contour map of the predicted value desirability calibrated with the dominant frequency of the acoustic wave (18 kHz) as the target value.Since the relationship between the factors is similar, here, as an example, we fix the value of Poisson's ratio to show the influence of shear modulus and coefficient of restitution on the confidence.It can be seen from the figure that the parameters in the red part have the highest fitting desirability, which is greater than 0.9.In addition, the figure shows that the red area is wide; that is, the shear modulus and coefficient of restitution can be selected in a wide range, which is 1700~2500 MPa and 0.35~0.7,respectively.However, the actual material parameters are uniquely fixed, so the reliability of the parameter combination is poor.Similarly, when the acoustic velocity (507 m/s) is used as an index for calibration (Figure 5), the parameter combination of soil particles can be selected in a wide range.Therefore, the selected simulation parameter values may differ from the measured material properties.The coefficient of restitution in the range of 0.3~0.7 can make the simulated acoustic velocity equal to the actual acoustic velocity (507 m/s).This will seriously affect the accuracy of parameter calibration and cause difficulties for subsequent simulation analysis.Figure 6 is the contour map of the desirability that is jointly calibrated with the dominant frequency and the velocity as the target.The figure must simultaneously satisfy the dominant frequency = 18 kHz and the velocity = 507 m/s.It can be seen from the figure that the area where the dominant frequency of the acoustic wave and the predicted value of the velocity desirability is above 0.9 is small (the red area is small).The range of parameters that can be selected is small, but there is no parameter combination with a desirability level of 1.Moreover, the optimal target value of dual-index calibration is acoustic velocity = 517.0m/s, dominant frequency = 17.5 kHz.The predicted optimal value of the acoustic velocity is greater than the actual value, and the predicted optimal value of the dominant frequency is smaller than the actual value.The above shows that the value range of parameter combinations is limited when the two indexes work together.Most parameter combinations cannot satisfy the two indexes simultaneously, thus narrowing the value range of parameter combinations and then excluding parameter combinations that do not meet the two indexes.It can be seen from the above analysis that when the single index is used for calibration, the desirability of the predicted value can reach 1; when the dual-index is used for calibration, the desirability of the predicted value can reach more than 0.9, but no parameter combination reaches 1.The predicted values of the dual-index calibration are acoustic velocity = 517.0m/s, dominant frequency = 17.5 kHz; the corresponding parameter combination is Poisson's ratio 0.14, shear modulus 1663 MPa, coefficient of restitution 0.7, and the desirability level is 0.964.We conducted a univariate analysis of the model to analyze the reasons for the low desirability of dual-index calibration.And then to improve the desirability of dual-index co-calibration.The specific method is: to fix the horizontal coding value of two of the factors to 0 and find the regression equation of the dominant frequency Y1 and velocity Y2 and the coding value of a single factor, as shown in Equations ( 20)-( 25): =15.52 1.81 Figures 7 and 8 are the influence curves of various factors on the dominant frequency Y1 and the velocity Y2 of the acoustic wave.Comparing the two figures, we can see that with the increase of Poisson's ratio (x1), the dominant frequency Y1 and the wave velocity Y2 both increase; with the increase of shear modulus (x2), the dominant frequency of acoustic wave Y1 and velocity Y2 all increase; with the increase of the coefficient of restitution (x4), the dominant frequency decreases, and the acoustic velocity Y2 increases.Univariate analysis shows that Poisson's ratio (x1), shear modulus (x2) and coefficient of restitution (x4) had significant effects on the dominant frequency and wave velocity.
The above results show that the desirability of the predicted value cannot reach 1 when the two indexes are jointly calibrated.The reason is that after calibration, the predicted acoustic velocity of the model is greater than the actual acoustic velocity, and the predicted dominant frequency is smaller than the actual dominant frequency.Then it is necessary to reduce the acoustic velocity and increase the dominant frequency at the same time.It can be seen from the single factor analysis results (Figures 7 and 8) that increasing the coefficient of restoration can reduce the acoustic velocity and increase the dominant frequency, so the range of the coefficient of restoration can be increased.We increased the coefficient of restitution from 0.3-0.7 to 0.3-0.8,resulting in a contour plot of desirability (Figure 9). Figure 9 that the desirability is increased to 1.At this time, the predicted value of dual-index calibration is acoustic velocity = 507 s, dominant frequency = 18 Hz, which is equal to the target value (index).Correspondingly, the parameter combination is Poisson's ratio 0.3, shear modulus 1470 MPa, and collision recovery coefficient 0.79.It shows that the optimized parameters can make the calibration result reach the target value.(2) Validation of the optimal combination of parameters This section verifies the optimal prediction results of two single and one dual index, respectively, and compares three.The optimal parameter combination was obtained by the optimization module, with the actual measured acoustic wave dominant frequency of 18 kHz and acoustic velocity of 507 m/s as the target values (index).We inputted the above optimal parameter combination in the EDEM 2018 software (A software for discrete element simulation.) and conducted a simulation experiment to verify the optimization results.The optimal parameter combination, model prediction results, simulation verification results, and errors between the simulation and experimental results are shown in Table 10.The calculation method of the error between the simulation verification and the experimental result is shown in Equation (26).It can be seen from the table that when the dominant frequency is used as the index.The errors between the simulation verification and the experimental results are as follows: the dominant frequency is 2.8%, and the velocity is 25.2%.The dominant frequency error of the acoustic wave is small, but the acoustic velocity error is significant, so the calibration result is poor.When the acoustic velocity is used as the index, the errors of the simulation verification and the measured results are, as follows: the dominant frequency of the acoustic wave is 35.6%, and the acoustic velocity is 0.14%.The dominant frequency error is small, but the velocity error is significant, so the calibration result is poor.In the case of dual-index calibration, the errors of simulation verification and experimental results are, respectively, as follows: the dominant frequency is 4.4%, and the acoustic velocity is 2.6%.The dominant frequency and velocity errors are minor, below 5%, within the acceptable range.
By observing the results of the two calibration methods, it can be seen that the singleindex calibration can satisfy a single target well.However, the error between the other index and the experimental value may be significant.The dual-index calibration can limit the two indexes to a small range, significantly reducing the range of parameter values and improving the accuracy of the calibration.

The Necessity of Constitutive Parameter Calibration
For granular bodies such as soil, only contact parameters are generally calibrated, while constitutive parameters such as Poisson's ratio, shear modulus, and density are generally obtained by measurement.In this paper, the constitutive and contact parameters are given a range value by referring to the measured experiments, and then the parameters are calibrated.That is, we still use the parameter calibration method to determine the parameters that can be directly measured.It is determined by the complexity of the soil medium and the inconsistency between the discrete element simulated particles and the actual soil particles.Soil is a complex mixture of various substances, including other substances such as organic matter, and we cannot include all substances in the model.The properties of soil particles and aggregates of various substances can only be represented by a single particle in the simulation.The properties exhibited by this aggregate are undoubtedly different from the actual individual soil particles.Therefore, even if the parameters of actual soil particles can be obtained, the phenomena displayed in the simulation may be different from those in the experiment.At the same time, to reduce the calculation, the particle diameter is generally enlarged, and uniform particle size is used.However, the diameter distribution of actual soil particles is random.Mapping each particle with actual soil particles one-to-one in the simulation is impossible so the simulation results may differ from the experiment.Therefore, it is necessary to adjust the parameter values in a certain range through calibration and finally obtain the macroscopic properties of the soil material and, at the same time, reflect the particle characteristics.

Difference of Calibration Effect
Here we discuss the influence of different scenarios, different indexes, and the different number of indexes on the calibration effect.Rui et al. [31] calibrated the parameters of sandy soil particles with the angle of repose as the index by combining experiments and simulations.The results are Poisson's ratio 0.3, shear modulus 1.15 × 10 7 , and soil-soil coefficient of restitution 0.15.Fei et al. [33] calibrated the soil parameters with the angle of repose as the index.The results are Poisson's ratio 0.3, shear modulus 5 × 10 7 , and soilsoil coefficient of restitution 0.21.The above are all calibrated with the angle of repose as the index.However, the optimal parameter combination is different, especially the shear modulus, which is about four times the former.For the same medium (soil), the final calibration parameters should be the same or have a slight difference.The above are all calibrated with the angle of repose as the index.However, the optimal parameter combination is different, especially the shear modulus, which is about four times the former.For the same medium (soil), the final calibration parameter combination should be the same or have a slight difference.This is because the parameters that significantly impact the angle of repose are contacting parameters such as the coefficient of restitution, coefficient of static friction, and coefficient of rolling friction.In contrast, constitutive parameters such as shear modulus have no significant effect on the angle of repose.Therefore, significant changes in shear do not significantly affect the angle of repose.The calibration results of Rui et al. [31] and Fei t al. [33] are also different from the results of this paper (Poisson's ratio 0.3, shear modulus 1.407 × 10 9 Pa, coefficient of restitution 0.79).It is because the indexes and the number of indexes are different in different scenarios, so the parameter combinations after calibration are different.Poisson's ratio, shear modulus, and collision recovery coefficient significantly impact acoustic wave propagation, so these parameters must meet the two indexes of acoustic dominant frequency and velocity.Furthermore, other parameters that have no significant effect on acoustic wave propagation are not the focus of attention, and their values can be selected within a reasonable range.It is different from the parameters (coefficient of restitution, coefficient of static friction, and coefficient of rolling friction) that are focused on when calibrating the angle of repose.Therefore, the indexes and the number of indexes to be satisfied in different scenarios are different, and the calibrated parameter values may differ.
It can be seen from the above that the angle of repose calibration method cannot meet the requirements of all application scenarios and cannot be used as a widely applicable calibration method for all studies.Those calibrations that do not consider application scenarios have little reference to practical applications and appropriate calibration indexes should be selected according to the application scenarios.Of course, the method in this paper cannot satisfy all scenarios but can only satisfy the application of acoustic wave simulation.We consider that this phenomenon of inconsistent parameter values after calibration is acceptable.Those parameters that have no significant effect on the research purpose are allowed to have a certain difference from the parameters of the actual material, provided that these parameters are determined to have no or no significant effect on the research purpose.

Conclusions
In this research, the Hertz-Mindlin (no-slip) model was used to calibrate the parameters of the acoustic-soil propagation mode.Sensitivity analysis was carried out with the characteristic value of acoustic wave (acoustic velocity and dominant frequency) in the experiment as the calibration index to determine the parameters that significantly impact the acoustic velocity and dominant frequency.The sensitive parameters were optimized to obtain the optimal parameter combination.The following conclusions can be drawn: The optimal solutions of the two models are obtained as follows: shear modulus is 1407 MPa, Poisson's ratio is 0.3, and coefficient of restitution is 0.79.The relative error values of the simulated dominant frequency and velocity and the actual values were 4.4% and 2.2%, respectively.It shows that the model constructed by the climbing test and the response surface test can well optimize the parameter values, and the optimal parameter combination can be used to study the acoustic wave propagation in the soil.(3) Comparing the calibration effect of the single and dual-indexes, it is found that a single index can only meet the calibrated index, while other indexes have significant errors.Moreover, the optional range of parameter values is wide and cannot be limited to a reasonable range.Dual-indexes can narrow the range of parameter values and meet the requirements of the two indexes.Therefore, the number of indexes should be increased as reasonably as possible to make the calibration effect more curate and the parameter values more aligned with the actual materials.(4) The calibration of different scenarios is compared and discussed in this paper.It can be seen that the calibrated parameter combinations are different in different scenarios, so we should select the significant parameters according to the scenarios, then calibrate the significant parameters, and finally determine the final parameter combination.
This study provides a calibration method for building a model of acoustic propagation in soil.It can improve the accuracy of model construction and provide help for analyzing the interaction mechanism between soil and acoustic waves.This method can also select specific indexes according to the scene and apply it to the calibration of other materials.

Figure 3 .
Figure 3. Fourier transform of the received signal.(a) Time-domain signal; (b) Frequency-domain signal.

Figure 4 .
Figure 4. Desirability contour map dominant frequency as an indicator.

Figure 6 .
Figure 6.Desirability contour map dominant frequency and velocity as indicators.3.2.3.Parameter Optimization and Verification of Desirability Based on Dual-Index (1) Parameter optimization for desirability

Figure 7 .
Figure 7. Curves of various factors on the dominant frequency.

Figure 8 .
Figure 8. Curves of various factors on the velocity.

Figure 9 .
Figure 9. Contour plot of desirability with the coefficient of restitution 0.3-0.8.

( 1 )
Based on EDEM software, the sensitivity analysis was carried out by the Plackett-Burman test.It was concluded that the sensitivity of each parameter to the dominant frequency of acoustic waves was ranked as Shear modulus, Poisson's ratio, Coefficient of restitution, Coefficient of rolling friction, Density, and Coefficient of static friction.The order of sensitivity to acoustic velocity is Shear modulus, Poisson's ratio, Coefficient of restitution, Coefficient of static friction, Coefficient of rolling friction, and Density.(2)Through the Box-Behnken test, the quadratic regression model of the three sensitive parameters on the dominant frequency and velocity is constructed and optimized.
R* represents the equivalent radius; δn represents the normal overlap; E* represents the Young's modulus, Rj and Ri represent the radius of the soil particles in contact; vj and vi represent the Poisson's ratios of the soil particles in contact; Ej and Ei represent the

Table 3 .
Factors and Levels of simulation test.

Table 4 .
Factors and levels of the Box-Behnken experimental.

Table 5 .
Design and results of Plackett-Burman experiments.

Table 6 .
Analysis of variance of Plackett-Burman experiment.

Table 7 .
The results of the sensitivity analysis.Establishmentof Regression Model and Significance AnalysisIn order to obtain the optimal combination of soil parameters, the Box-Behnken test was carried out, and the program and results are shown in Table8.Quadratic regression fitting was performed on the experimental data in Table8by Design-Expert 8.0.6 software.The regression equation with the dominant frequency Y1 and the acoustic velocity Y2 as the response values and Poisson's ratio (X1), Shear modulus (X2), and Coefficient of restitution (X4) as the variables are obtained.

Table 8 .
The design scheme and results of Box-Behnken.

Table 9 .
ANOVA of Box-Behnken design quadratic model.

Table 10 .
Model prediction optimal parameter combination, prediction results, simulation verification results and errors between simulation verification and experimental results.