Non-Contact in Situ Estimation of Soil Structure Using Acoustic Reflections

1. Introduction

Soil is the connection between the atmosphere and ground water, root systems, and sub-soil biological and chemical processes. Central to this connection is the degree to which air and water can penetrate the soil. This is primarily measured by the soil porosity, 𝜙, which is the ratio of the volume of spaces between soil particles to the total volume of a soil sample. Soil porosity significantly affects many important soil functions such as water, and air transmission, as well as water storage and availability to crops and pasture [1]. These soil functions in turn determine movement of nutrients and contaminants to water bodies, crops and pasture growth, and greenhouse gas emission [2]. Soil compaction, caused by activities like machinery use or animal traffic, reduces soil porosity by reducing the size and number of pore spaces [3]. Compacted soils have lower infiltration rates, poor drainage, and limited root penetration, leading to decreased plant productivity [4].

The spaces, or pores, between soil particles are not simple cylinders, leading to further soil physical parameters tortuosity, α_∞ (the length of a pore following its path, divided by the shortest distance between the two ends of the pore), and the pore radius r_pore. In practice there will be a wide range of pore geometries but, to first order, a single tortuosity and characteristic pore radius can be assigned to a soil sample.

Soil porosity is known to vary strongly and unpredictably in space and time as a result of management, such as the timing of tillage, crop management, and grazing practices [3,5]. Traditional methods for measuring soil porosity include: measuring the extra mass of water which is required to saturate a soil sample; use of a pycnometer to measure the air volume in the pore space for a soil sample in a gas-tight chamber; and compression of a soil sample to estimate the solid volume of soil [6]. Such methods require significant investments in time and labour [7]. Moreover, these existing methods provide fragmentary information in time and space and offer a limited capacity to assess temporal and spatial variation in soil health. In contrast, proximal sensing technologies have the potential to generate vastly more data at lower cost, providing that the complexity of such technologies do not present a barrier to adoption by land managers such as farmers.

Currently, soil pore size and pore connectivity measurement methods also rely on invasive methods [8]. For example, details of pore connectivity can be measured via X-ray imagery [9] but this requires specialised equipment in a laboratory environment. Nevertheless, pore size is a significant parameter in estimation of hydraulic flow [10].

Bradley et al. [11,12] developed a non-contact method for estimation of soil porosity and tortuosity based on the strength of reflections of ultrasound pulses. Their instrument could in principle be mounted on a farm vehicle, typically 1 m above the ground. At lower transmitted acoustic frequencies viscous interactions within the pores are dominant, whereas at high frequencies the dependence on air viscosity and on flow resistivity vanishes. These interactions of sound waves with the porous soil lead to a range of in-situ and laboratory methods which can be used to measure the morphological characteristics of porous materials [8]. Of these, only the reflection methods, which are applicable at high frequencies, lend themselves to proximal sensing. In this regime, the plane wave reflection coefficient depends only on angle of incidence, soil porosity, and soil tortuosity (a measure of the soil pore geometry), and not on the frequency of sound transmitted. By performing measurements at two or more angles of incidence, it is possible to estimate both porosity and tortuosity [13,14,15].

At low acoustic frequencies another parameter, the flow resistivity, can be estimated via reflection of sound [16]. The flow resistivity is related to the characteristic or average pore radius [8]. This suggests that a combination of ultrasound and low acoustic frequencies could yield estimates of porosity, tortuosity, and pore radius. However, the acoustic frequencies used for flow resistivity estimation are in the range 35– 75 Hz so that measurements are in the low-frequency viscous regime [16] and the propagation path was 50 m. This is problematic operationally for a proximal device operating at only 1 m from the target soil surface, given the wavelengths are in the range 4 – 10 m.

The purpose of the current work is to explore the potential of proximal sensing of soil pore radius in addition to porosity and tortuosity. This requires an extension of the previous work by Bradley et al. [11,12] to include frequency-dependence and, based on the above discussion, working with frequencies which are operationally practical. The goal is to develop a technology which can be applied at paddock or farm scale, based on a sensor mounted on a small farm vehicle.

2. Theory and Simulations

The aim is to estimate soil physical properties by using sound reflected from the soil surface, without disturbing the soil. This potentially will allow rapid assessment of soil physical health from a farm vehicle or drone. The proposed operational setup is shown schematically in Figure 1. A sound pulse of frequency f is incident on the soil surface at an angle of incidence θ from a directional speaker, and the reflected sound measured in the specular direction by a directional microphone. The ratio of the reflected pulse and the transmitted pulse is the plane wave reflection coefficient R. Some sound is transmitted into the soil so the amplitude |R| ≤ 1. In general, there is also a phase change upon reflection, including absorption.

2.1. Approximations for the Plane Wave Reflection Coefficient R

Bradley et al. [11,12] express R in the form

$${\displaystyle \frac{1-R}{1+R}}={\displaystyle \frac{{\rho }_0}{{\rho }_e}}\sqrt{{\displaystyle \frac{{\rho }_e}{{\rho }_0}}{\displaystyle \frac{K_a}{K_e}}\left(1+X\right)-X},$$

(1)

where

$$X={\tan }^2\theta$$

(2)

𝜙 is the soil porosity, ρ_e is an equivalent soil density and K_e the effective bulk density (allowing for the soil porosity), K_a is the bulk modulus of air, and ρ₀ the air density.

Brennan and To [17] show that

$${\displaystyle \frac{{\rho }_e}{{\alpha }_{\infty }{\rho }_0}}\approx 1+{\displaystyle \frac{1}{i}}$$

(3)

and

$${\displaystyle \frac{K_e}{K_a}}\approx 1$$

(4)

where the normalized frequency is

$$\Omega ={\displaystyle \frac{f}{f_c}}.$$

(5)

The critical frequency f_c is given by

$$f_c={\displaystyle \frac{4}{\pi {\rho }_0{r}_{pore}^2}}$$

(6)

where η is the dynamic viscosity of air and r_pore is the characteristic (or mean) pore radius [6,8]. Note that, when the focus of research is on acoustics, the flow resistivity is generally used instead of r_pore. The result of these approximations is

$${\displaystyle \frac{1-R}{1+R}}=\sqrt{c\gamma +d{\gamma }^2}.$$

(7)

where

$$c=A\left(1+X\right)$$

(8)

$$d=-A{\displaystyle \frac{X}{{\alpha }_{\infty }}}$$

(9)

$$\gamma ={\displaystyle \frac{i}{1+i}}$$

(10)

and

$$A={\displaystyle \frac{{𝜙}^2}{{\alpha }_{\infty }}}.$$

(11)

For Ω >> 1

$${\left({\displaystyle \frac{1-{\left|R\right|}_{\to \infty }}{1+{\left|R\right|}_{\to \infty }}}\right)}^2=c+d=a=A+BX$$

(12)

where

$$B=A{\displaystyle \frac{{\alpha }_{\infty }-1}{{\alpha }_{\infty }}}.$$

(13)

This approximation was successfully used by Bradley and Ghimire [11] and Bradley et al. [12] to estimate soil porosity 𝜙 and soil tortuosity α_∞ using the slope and intercept of a linear regression from measurements |R| at several angles of incidence θ and at an ultrasonic frequency of 25 kHz (see also [13,14,15]).

For Ω << 1

$${\left({\displaystyle \frac{1-{\left|R\right|}_{\to 0}}{1+{\left|R\right|}_{\to 0}}}\right)}^2=c$$

(14)

This was the asymptotic solution used by Sebaa et al. [16] for which Ω was in the range 0.035 to 0.075. Equation (14) suggests that, once 𝜙 and α_∞ have been estimated from transmissions such that Ω >> 1 using Equation (12), a linear regression involving measurements at several transmissions such that Ω << 1 will allow an estimation of Ω and hence r_pore.

2.2. Typical Values of Model Parameters

The strategy outlined above requires use of well-separated lower acoustic frequencies, f_low, and higher acoustic frequencies, f_high, satisfying

$$f_{low}\ll f_c\ll f_{high}$$

(15)

Table 1 gives typical values of parameters appearing in the above theory. Being simply a fraction of total volume, 𝜙 can range between 0 and 1, typically falling between 0.35 and 0.7 for soils [6]. Typical values for grasslands quoted by Salomons [18] are α_∞ = 1.35, 𝜙 = 0.3, and 5 ≤ Ω ≤ 20 when the sound has a frequency of 25 kHz, equating to r_pore = 63 – 125 μm.

The limitation on the angle of incidence arises from keeping the instrumentation reasonably compact. The limitations on transmitted frequency arise from aiming for directionality whilst not being too sensitive to scattering of sound by vegetation and surface roughness at higher frequencies. For a propagation path from speaker to ground to microphone of length L = 2 m (which would be a typical range from a small agricultural vehicle) and for an angle of incidence θ = 30°, the time for sound to travel horizontally between speaker and microphone is L/(2c₀) = 3 ms if the sound speed c₀ = 340 m s^-1 (see Fig. 1). The time of first arrival of reflected sound is L/c₀ = 6 ms. This means the duration of the acoustic pulse should be less than 3 ms. There are 12 cycles of 4 kHz sound in a 3 ms pulse, which would give a reasonable pulse shape for detection and peak analysis. For these parameters, one measurement would take L/c₀ + L/(2c₀) = 9 ms, so that around 100 measurements per second would be feasible if transmissions at all frequencies and angles of incidence simultaneous. However, an acoustic frequency of 4 kHz is considerably higher than the lowest likely critical frequency of 0.5 kHz. Similarly, the high acoustic frequency of 25 kHz is not much higher than the highest likely critical frequency of 22 kHz.

For example, assuming estimates of tortuosity and porosity have been made at higher frequencies, Figure 2 shows significant errors in estimated pore radius if Equation (14) is used with transmitted acoustic frequencies in the range of 1 – 5 kHz. There would also be concerns about external acoustic noise from, for example, engines, at these lower acoustic frequencies.

2.3. Approximation at Higher Frequencies

An alternative is to use higher acoustic frequencies and allow for the frequency-dependence of Y, where

$$Y={\left({\displaystyle \frac{1-\left|R\right|}{1+\left|R\right|}}\right)}^2$$

(16)

Figure 3 shows the variation of Y with f for frequencies in the range 4 to 16 kHz. Also shown are the Y values at 1/3-octave frequencies in this band. The asymptotic high-frequency limit of these curves should give A and B and hence porosity 𝜙 and tortuosity α_∞, and the curvature should give an estimate of Ω and hence f_c and r_pore.

The curvature can be found from a Taylor series expansion of the right-hand side of (7) in terms of Ω ^-1, Ω ^-2, ..., followed by expanding R from (7) as a series in in terms of Ω ^-1, Ω ^-2, ..., then expanding |R|, and finally Y from (16). This expansion is non-trivial. The result is

$$Y_{\to \infty }=a+{\displaystyle \frac{k}{2^2}}+\dots$$

(17)

where

$$k=2a{\displaystyle \frac{8a^2-c^2\left(1-3a\right)-4ca\left(1+a\right)}{4\left(a-1\right)a^2}}$$

(18)

2.4. Procedure for Estimation of Soil Physical Parameters

We measure $\left|R_{mn}\right|$ at an angle of incidence θ_m and at a transmitted frequency f_n, where m = 1, 2, …, M and n = 1, 2,…, N.

For each θ_m perform a regression of Y vs f ^-2 as in (17), and using (5), to find estimates ${\widehat{a}}_m$ and ${\widehat{b}}_m$ of the coefficients, m = 1, 2, …, M, where

$$Y_{mn}={\left({\displaystyle \frac{1-\left|R_{mn}\right|}{1+\left|R_{mn}\right|}}\right)}^2=a_m+b_m{f}_n^{-2}+\dots .$$

(19)

There are further expansion terms in (17) and (19) and, in order that κ or b_m are representative of the first expansion term, we perform a polynomial fit of order 4. Given there are N = 7 frequencies, there is still sufficient redundancy in the regression. Using the M estimated values of ${\widehat{a}}_m$ and the M known values X_m do a single linear regression of a vs X

$$a_m=A+B{X}_m$$

(20)

and from the two coefficients $\widehat{A}$ and $\widehat{B}$ estimate

$$\widehat{𝜙}={\displaystyle \frac{\widehat{A}}{\sqrt{\widehat{A}-\widehat{B}}}}$$

(21)

$${\widehat{\alpha }}_{\infty }={\displaystyle \frac{\widehat{A}}{\widehat{A}-\widehat{B}}}$$

(22)

It is now possible to estimate $\widehat{k}$ from (18), and

$${\widehat{f}}_{c,m}=\sqrt{2{\displaystyle \frac{{\widehat{b}}_m}{{\widehat{k}}_m}}}.$$

(23)

Hence estimate the characteristic pore radius

$${\widehat{r}}_{pore,m}=\sqrt{{\displaystyle \frac{4}{\pi {{\rho }_0\widehat{f}}_{c,m}}}}.$$

(24)

Note that M estimates are made of the pore radius, which are simply averaged.

3. Results

As an example, θ_m = 0°, 10°, 20°, and 30° for m = 1, 2, 3, 4 (M = 4), and f_n = 4.00, 5.04, 6.35, 8.00, 10.08, 12.7, and 16.00 kHz for n = 1, 2, 3, 4, 5, 6, 7 (N = 7), giving MN = 28 measurements $\left|R_{mn}\right|$. The acoustic frequencies are in 1/3 octave steps. Values of $\left|R_{mn}\right|$ with this set of θ and f values are given in Table 2, for porosity = 0.6, tortuosity = 1.6, and pore radius = 80 μm.

Simulations are performed for 5 soil porosity values of 0.4, 0.5, 0.6, 0.7, and 0.8, 5 tortuosity values of 1.4, 1.5, 1.6, 1.7, and 1.8, and 130 pore radius values from 50 μm to 179 μm in 1 μm steps. Values of $\left|R_{mn}\right|$ are calculated for these 3250 simulated soil profiles, and for the M angles of incidence and N frequencies. Random gaussian noise ε_g with g = 1, 2, ….G, and having zero mean and standard deviation σ_R is added to each simulated $\left|R_{mn}\right|$

$${\left|R_{mn}\right|}_g=\left|R_{mn}\right|\left(1+{\epsilon }_g\right)$$

(25)

where ε is a random gaussian number of zero mean and standard deviation σ_R. With G = 1000 there are 3.2 x 10⁶ simulated reflectivity measurements. A value σ_R = 0.01 (1% variability in measured |R|) was used, consistent with actual measurements [11].

For each soil “profile”, or group of simulation 𝜙, α_∞, and r_pore values, estimates $\widehat{𝜙}$, ${\widehat{\alpha }}_{\infty }$, and ${\widehat{r}}_{pore}$ are made from the G noisy “measurements”. Mean values and standard deviations of estimates are found for each of the 3250 simulated soil profiles.

There are two main sources of error: model error and random measurement noise error (systematic measurement errors can be made negligible with good hardware design and calibration). Model error manifests as bias in estimates, whereas noise error manifests as both bias and random fluctuations in estimates. Bias β is estimated from

$$\beta ={\displaystyle \frac{\widehat{𝜙}}{𝜙}}-1$$

(26)

and random uncertainty s from

$$s={\displaystyle \frac{{\sigma }_{\widehat{𝜙}}}{𝜙}}$$

(27)

and similarly for tortuosity and pore radius. Values of these measures averaged across all soil profiles are shown in Table 3.

The parameter estimation errors shown in Table 3 give no indication of whether more accurate estimates might be obtained from certain combinations of soil parameters. In order to answer this, β and s are averaged over 2 of the 3 soil parameters giving, for example, β_{𝜙, ∞}(r_pore) for β averaged over all 𝜙 and α_∞ as a function of r_pore . The results are plotted in Figure 4 and Figure 5. It is clear that the bias and random errors are low (a few %) and largely independent of soil properties except for a strong dependence on lower r_pore values. These 3 (solid line) bias curves reach 2% for r_pore smaller than about 60 μm. Estimates of all three soil parameters have significant errors for low pore size. Bias errors are also simulated with zero random noise, giving a very similar plot to Figure 4. This indicates that the bias errors are almost entirely due to model errors. Since these errors are systematic, they could be corrected by fitting curves to the error curves, although we have not attempted that here since the errors are very small.

Random fluctuation errors only appear if σ_R > 0. Tortuosity estimates are sensitive to low pore radius, reaching around 1.5% fluctuation errors for r_pore = 50 μm. Random errors of this magnitude are acceptable.

4. Discussion

The approach taken here to estimate soil physical properties is very different from traditional laboratory-based methods. Our method allows fast non-contact in situ simultaneous estimates of soil porosity, tortuosity and soil pore radius. This method also differs dramatically from the asymptotic acoustic methods used by others in either a very high-frequency domain or a very low-frequency domain. The key new element is the expansion of the full expression for the plan wave reflection coefficient R in a Taylor series with terms in 1/Ω ², 1/Ω ⁴,… Only the coefficient of the first expansion term (in 1/Ω ²) has been developed, yielding an expression for the pore radius. A further linear regression is required to obtain estimates of porosity and tortuosity from the same set of measurements.

A Monte Carlo simulation shows that the primary source of uncertainty in soil parameter estimation is model error. These errors are very small, being typically less than 1 %, and arise from truncation of the Taylor series expansion to 4 terms (proportional to 1, 1/Ω ², 1/Ω ⁴, and 1/Ω ⁶). The expansion is truncated so that the polynomial regression is over-determined with 7 acoustic frequencies and the 4 expansion coefficients.

In practice this analysis compares one model (the Taylor series approach) with another which assumes (1) – (6) are valid. The robustness of this approach therefore depends on

• Accuracy of (1) to (5). This has been validated by numerous investigations (see, for example, [17,19,20,21,22,23]).
• Accuracy of (6). Horoshenkov [8] used this relationship to estimate the radius of glass beads based on known packing, obtaining agreement to within 9%. This is equivalent to a 9% standard deviation in estimated pore radius.
• The assumed value σ_R = 0.01 for the Monte Carlo simulations. Bradley et al. [12] show an example of a received pulse in which the rms amplitude noise is around 3 mV on a signal of amplitude 80 mV. A group of 16 such pulses are used to fit a known pulse shape to the data. If the pulse peak is used, without any pulse shape fitting, the relative standard deviation is σ_R = (3/80)/4 or 0.9%. In practice fitting a known shape to the pulse data will give a smaller uncertainty in |R|. Uncertainties in estimated parameters simply scale with σ_R at these low levels.
• The assumption that errors due to signal loss in sound passing through grass and roughness are small. The relevant scattering parameter is 2πfh/c₀ where f is the acoustic frequency and h is the dimension of a scattering object. For pasture swards h is of order 5 mm and 2πfh/c₀ is around 1. At this value of the scattering parameter some loss of reflected energy will occur [24]. For soil particles forming a rough soil surface, the scattering parameter will generally be smaller [12]. For pugging due to animal hoof prints, major scattering can be expected. In all these cases the best operational strategy may be to detect scattering at each angle of incidence using the multiple microphones and, if scattered energy is more than a few percent of the incident energy, discard that data point [12].

These caveats do not appear to limit the scope and credibility of the new approach taken here. Of course, design of hardware and software generally also involves compromises, but previous work [11,24] gives confidence that a suitable design should be practical. This is particularly true given the availability of a much wider range of transmitting and receiving devices in the audio spectrum compared with the ultrasonic region. The new methodology also lends itself much more to being controlled entirely from a mobile phone because the audio interfaces on modern phones have the requisite bandwidth and sampling rates. The ultimate goal is therefore a compact rigid array of speakers and microphones which can be mounted on a small farm vehicle such as a quad-bike, with data recorded on a mobile phone and able to be readily uploaded to a farm management system.