1.

INTRODUCTION

Following the Nobel prizing winning work of Steven Chu et al,¹ the use of cold atom matter-wave interferometry to measure inertial forces has developed from lab-based measurements of fundamental physical constants to deployable sensors for high precision measurement of acceleration and rotation. The current capabilities of atom interferometers for measuring inertial forces are comparable to their state-of-the-art classical counterparts.² However, in principle their precision can be improved by several orders of magnitude.³ With the potential for greatly improved sensitivity to tiny changes in acceleration, applications of gravimeters and gravity gradiometers include: improved identification of oil reservoirs; the detection of sewage services and sink-holes, and civil engineering more generally; detection of underground buildings; the study of geophysical phenomena; and geophysical mapping as a means to navigate in GPS denied environments.⁴

In order to meet the requirements of these applications, signal processing of the raw gravity data is necessary; this is therefore a geophysical inversion problem. Many excellent texts describe the problems and challenges of determining the composition and morphology of the sub-surface (see for example^{5, 6}). For the problems under consideration it is clear that a method that accounts for all sources of uncertainty and is flexible enough to accommodate prior knowledge is critical. A Bayesian approach to this problem is therefore a natural choice. As an example, even for the simple case of a gravitational anomaly caused by a single void, it is difficult to differentiate between the void’s volume, its depth beneath the surface and the density contrast with respect to the soil. This is complicated greatly by the inhomogeneity of the soil; the location of interfaces between different rock and soil types; and numerous other environmental sources: tidal variation, water-table height variation, atmospheric pressure, seismic waves, terrain effects and other buildings, noise generated by the sensor itself, and even vehicles and living creatures.⁷

For the Defence and Security led challenge of imaging underground buildings and tunnels at depths not amenable to current survey technology such as ground penetrating radar, the probability of excavation (POE) map was introduced.⁸ This paper builds on this work by attempting to account for variations in the composition of the soil. We show that variations in soil density will generate spatially correlated gravitational fields at the surface of the Earth. By inferring a spatial length scale that accounts for this, we compare the posterior distributions, and resulting POES maps, with and without this additional parameter and speculate on its utility in future measurement campaigns.

2.

MEASURING SMALL CHANGES IN GRAVITY

To assess the capability of atom interferometers to locate underground voids, it is instructive to take the demonstrated precision of cold-atom gravimeters at 8ng/Hz^1/2 (1ng is 9.81 × 10⁻⁹m/s²), and compare this to typical signal strengths.⁹ For example, the gravitational acceleration exerted by a 150kg mass at 1m is approximately 1ng (so it is possible to detect the sensor operator if the sensor makes an integrated measurement for approximately a minute). For the detection of voids beneath the Earth a simple comparison with the analytic expression for a spherical void can be made (see⁷ page 49). Figure 1 shows the expected magnitude of acceleration anomaly caused by a 100m³ spherical void at 10m and 5m depth, compared to the noise floor of the atom interferometer of approximately 1ng. The local density of earth is assumed to be 1800kg/m³ (see⁷ for representative values). From this simple assessment it can be concluded that atom interferometry based technology is capable of detecting the presence of such voids. However, this result is without due consideration of sources of environmental and geophysical noise encountered in the real world. The spatial resolution of more complex cuboid voids, and the consideration of background noise sources, will be described in terms of a Bayesian model in the following section.

Figure 1.

Magnitude of the anomalous gravitational acceleration, g_z, of a spherical void as a function of horizontal distance from the cavity’s centroid, x, and depth beneath the Earth, z. The cavity has a volume, V, of 100m³, the density contrast, Δρ, is 1800kg/m³, and G is the gravitational constant. The noise floor of the atom interferometer gravimeter, 1.03ng, is calculated assuming a minute of integration at 8ng/Hz^1/2 precision.

3.

A BAYESIAN APPROACH TO GRAVITY ANOMALIES

To determine the spatial extent of a void using geophysical measurements such as gravity, there exists extensive literature on the use of methods in linear algebra.⁶ Whilst these are powerful techniques and can quickly provide a solution, they often suffer from great sensitivity to small changes in the measurement data. The geophysical inversion problem is generally known to be ill-posed so methods are applied to regularise the equations to provide stability (see for example Tikhonov regularisation⁵). In contrast, a Bayesian formulation of the problem provides a framework for accounting for all sources of error and prior knowledge; it is not limited by linear dependence on model parameters.

The Bayesian solution to the inversion problem is given by the posterior probability density function (pdf): the space of all possible solutions describable by the model, m, (see for example^{5, 10}). The posterior pdf, p(θ | y, m), over model parameters, θ, given the data y, is related to the prior, p(θ | m), through:

where p(y | θ, m) is the likelihood function, and the denominator can be regarded as a normalising constant for the purpose of computing marginal distributions and, for example, probability of excavation maps.

3.1

Bayesian model

The model used throughout this paper is given in Figure 2. It describes the dependencies between the measurement data, measured stochastic variables, and the parameters, unobserved stochastic variables. These latent variables describe a single cuboid void within a heterogeneous medium of soil specified by a soil density fluctuation (see section 5 for details of the soil noise model). Whilst not accounting for access, this would be typical of a simple underground building and is a good starting point before moving to more complex geometries and correspondingly complex models. A simpler version of this model was first reported in.⁸

Figure 2.

Gravimeter model with spatially correlated soil noise. Circles define stochastic variables: white are inferred, grey are observed. The diamond denotes the deterministic function g_z (see Equation 4). The observed data, y, is modelled as a multi-variate normal distribution, , of mean vector μ and covariance matrix Σ. The mean vector is determined by the background offset, η, the parameters of the void, and implicitly the measurement locations associated with y. Similarly the covariance matrix has an implicit dependence on the measurement locations and is defined by the soil noise parameter, ξ, and the total noise parameter σ². The void parameters are: cuboid centroid x, cuboid dimensions ℓ, cuboid orientation in the horizontal plane ϕ, and void density contrast Δρ.

Whilst this model does not directly consider terrain effects and rock soil interfaces that would give rise to very strong spatial variations in the background gravity field, it does account for random variations in soil density through the delta-correlated soil noise model described in section 5. Essentially random uncorrelated spatial variations in the soil give rise to spatial correlations in the gravitational field at the surface. These spatial correlations can be described in terms of the soil density fluctuation parameter.

3.2

The prior distribution

To proceed with Bayesian inference we define distributions for each stochastic variable of θ. Following work described in,^{8, 11} one simple choice is to assign uniform priors for each parameter. However, despite their apparent attractiveness as non-informative (of course they are not really non-informative since a range must be defined), uniform priors can become a problem if the true posterior distribution is non-vanishing outside the prior’s support. In the event of this occurring the convergence of the Markov chain Monte Carlo (MCMC) sampler is ruined: it is effectively not possible for the sampler to explore the entire parameter space and thereby the simulation is not ergodic. Instead we choose to use informative priors and compare with the resulting posterior distribution.

For the scenario to be analysed, a set of 121 measurements over a 20m×20m flat Earth grid serves as the gravity data to be processed, see subsection 3.4 for more detail. The target gravity anomaly is assumed to be a man-made structure and corresponds with a cuboid void. This simplistic model assumes no other localised gravity anomalies other than that caused by natural variations in the soil density. The cuboid void is described in terms of its centroid and its height, length and width. However this parametrisation suffers from strong correlations between its variables, which is known to slow convergence in MCMC Metropolis-Hastings sampling (see for example¹⁰). To alleviate this problem, and also incorporate more knowledge about the target anomaly, we re-parametrise the model as:

where is the distance from the co-ordinate system origin, defined as the centre of mass of the measurement locations, to the centroid of the cuboid anomaly. Thus, z₀, the z component of the centroid is defined to be a fraction, α, of the distance d. Similarly the x and y components are calculated using the angle θ. The cuboid volume, V, is defined to be quadratically dependent on the distance d, following the inverse square law, and the scaling parameter ν. The height of the cuboid, l_z, is some fraction, β, of the maximum height 2z₀ (defined such that the cuboid doesn’t extend above the surface of the Earth). The cuboid length, l_x is defined to be larger than the width, l_y, and this is ensured by the positive definite parameter γ. We specify the length parameter to be larger than the width parameter in order to prevent multiple posterior maxima with respect to cuboid orientation in the horizontal plane ϕ.

The model prior takes the following form:

where denotes a normal distribution of mean, μ, and variance, σ². The uniform distribution, , is bounded between L and U. Many of the stochastic variables are described in terms of Equation 2; the remaining parameters are cuboid orientation in the horizontal plane ϕ, void density contrast Δρ, and background offset η.

The choice of the form of the priors is based on the concept of reposing the problem in terms of bounded multiplicative factors where possible. Distributions of this form are generally easier to sample and therefore improve convergence. The remaining unbounded parameters have been assigned to the simplest possible log-normal distribution. The rationale for this is twofold: these parameters should not take on the value of 0 as this is not part of the model space; and the expected value for these parameters is of order 1. The normal distributions have been used for the parameters where relatively strong prior knowledge is assumed to be available.

The distribution for Δρ is based on an assumed mean density of 1800kg.m⁻³ with a relatively small uncertainty of ±50kg.m⁻³. This is a deliberately informative prior as, in the event of no knowledge about the mean density of the soil, there is a strong degeneracy with anomaly volume and depth beneath the Earth. This would make the resulting inversion more ambiguous and widen the solution space. Future work will consider this factor in more detail. However, in practice, we are assuming that we have core samples of the local environment and therefore have some knowledge of the local mean density.

The total noise parameter, σ², corresponds with the diagonal elements of the covariance matrix shown in Figure 2. This includes all sources of noise, specifically sensor noise and soil noise due to density fluctuations. The soil noise parameter, ξ, is defined to be some fraction of the total noise parameter. It accounts for the magnitudes of the off-diagonal elements of the covariance matrix (see also Equation 5).

The background offset, η, accommodates systematic shifts on all measured data as caused by, for example, water table and tidal variations. Such effects typically occur on time-scales much greater than the sensor integration time and their spatial length scales are typically much greater than the domain over which inference is being made. This may not be true for large scale geophysical inversion problems, but is a reasonable approximation for inferring the location of underground voids as described in this paper. This background offset is not expected to vary by more than a few ng, so the assigned normal distribution’s variance is set to 100ng².

3.3

The likelihood function

The likelihood function is given by a multivariate normal distribution centred on the forward model.¹² The forward model gives the expected z-component of acceleration due to a cuboid anomaly in the absence of other sources, at the ith measurement location x_i = (x_i, y_i, z_i), as

The gravitational constant is given by G, η is the background offset, x = (x₀, y₀, z₀) is the centroid of the cuboid and the length, width and height of the cuboid are specified by the components of ℓ. Dependence on the cuboid’s orientation in the horizontal plane, ϕ, is implemented by rotating the co-ordinates of the measurement locations with respect to the co-ordinate system of the cuboid.

The measurement covariance matrix is given by

where W_ij is the off-diagonal matrix element between the ith and jth measurement locations x_i and x_j. The functional form for W_ij is derived from the delta-correlated soil noise model described in section 5. Note that stricty all z_is must be equal for this equation to be valid i.e. all sensor measurements should be made in a plane. The likelihood function is then given as

where Y is the measurement location matrix, R is the rotation operator in the horizontal plane, and g_z is the vectorised form of Equation 4 acting on the measurement location matrix.

3.4

The data

A small computer experiment was implemented to address the question: how important is the effect of spatially correlated gravitational noise for resolving underground structures? A set of computer generated surveys were created with the gravitational field composed of a single cuboid anomaly embedded in a soil with density variations described in section 5.

The computer-generated gravity signal was produced from an underground region with rectangular shape and dimensions 30m by 30m in the horizontal plane and 20m deep. The Earths field is not included, so the assumption is that this has been corrected for in the raw data. The underground region is divided into cubes with sides of 20cm; thus the model has 2.25 × 10⁶ volume elements. Each element is assigned a density value and its gravity signal is given by Equation 4. The data is taken in a horizontal plane 1m above the ground and each data point is a sum of the gravity signals from all of the volume elements. In order to reduce the time taken for the calculation, a convolution method was used whereby an array of density values is convolved with the individual cuboid signal of Equation 4; this approach only requires the formula Equation 4 to be calculated once and benefits from the speed advantages of fast Fourier transforms that are used to compute the convolution in the spatial frequency domain. It is possible to do a full 3D convolution, but in order to reduce memory requirements, and taking advantage of the fact that the data is only required in a single horizontal plane, the calculation was split into a series of 2D convolutions: one for each of the 100 horizontal planes within the underground region.

The cuboid void is taken to have sides parallel to those of the underground rectangular region. In order to allow an arbitrary angle in the horizontal plane between the cuboid and the data array, the results are interpolated and then resampled on an appropriately rotated grid of dimensions 20m×20m. For convenience, all volume elements lying within the void were given a negative density of -1800kg/m³ and elements outside of the void had an average density of zero. Compared to the alternative of zero density for the void and positive density for the surrounds, the only effect is a uniform offset on the data which has no influence on the final result. A very simple model was used for random underground density variations (i.e. soil noise): uncorrelated, Gaussian variations of density in each volume element. This approach has the advantages of being simple to implement and allowing for straightforward calculation of the correlations produced in the gravity signal by the soil noise (see section 5) which can then be used in the Bayesian inference. However, it is unlikely to be a good description of real underground density variations. These would be expected to have some degree of spatial correlation and probably vary widely in different locations with different geology, soil types, etc. There is some evidence in the literature for power-law correlations in underground density variations,¹³ but much of this data comes from long vertical bore-holes and it is not clear how relevant this is to the near-surface scenario considered here.

As discussed in section 5 below, the soil noise can be characterised by the parameter d₀, which has units of kg.m^−3/2. Figure 3 shows the void used in these calculations and the resulting gravity anomaly both with and without the soil noise. In this case the magnitude of the soil noise is 300kg.m^−3/2. This level of soil noise was chosen to give a relatively large perturbation to the signal from the void, which has a peak value of 2.6 × 10⁻⁸m.s⁻². In contrast, the standard deviation of the additive sensor noise was much less than the peak from the void, at 5 × 10⁻⁹m.s⁻².

Figure 3.

Underground void (top left). Gravity signal from void alone (top right). Gravity signal with 300kg.m^−3/2 soil noise (bottom left). Signal down-sampled to 2m spacing with 0.5ng sensor noise (bottom right).

Ten realisations of this scenario were generated for two different values of the soil density fluctuation parameter: 50kg.m^−3/2 and 300kg.m^−3/2. For each realisation, a regular 20m×20m array of 121 measurements of the z-component of the gravitational field at 1m from the surface of the Earth was simulated.

3.5

Estimating the posterior distribution

The computation of an approximation to the posterior, p(θ | y, m), was performed using standard Metropolis-Hastings Markov chain Monte Carlo (MCMC) as described in,^{10, 14} implemented in PyMC.¹⁵ Each posterior pdf was calculated using 6 independent MCMC chains of 60,000 MCMC iterations. For each chain, a starting location was randomly selected from the prior distribution. A burn-in period of 30,000 simulation steps with adaptive sampling was performed, followed by Metropolis-Hastings without adaptation. A thinning of 1 in 10 left 3,000 samples per chain. The 3,000 samples from each independent chain were summed together providing a total of 18,000 samples to compute the posterior distribution and to test convergence. Figure 4 compares the traces resulting from running the inversion algorithm with two different models: ξ ≠ 0 and ξ = 0, on the same data generated with a soil density fluctuation of 300kg.m^−3/2.

Figure 4.

MCMC traces of volume (left) and noise standard error (right). The top traces corresponds with the model ξ ≠ 0, the bottom traces ξ = 0. The processed measurement data was generated with a soil density fluctuation of 300kg.m^−3/2. For each parameter the mean value is shown by the black line and the dotted lines are the 5% and 95% population boundaries. The true value is shown in red.

In order to draw valid conclusions from these calculations the MCMC simulations must approach convergence. Therefore it is necessary to test whether we are approximately sampling from the posterior distribution and whether we have enough independent samples. We compute the potential scale reduction factor (PSRF) statistic as an indicator of convergence. ¹⁰ From Table 1 and Figure 4, there is strong evidence that the inference where ξ ≠ 0 has approached convergence. In contrast, the model with ξ = 0 has not fully converged. Whilst some parameters, such as the location of the cuboid in the xy plane are stable, the more difficult parameters associated with the volume and depth of the anomaly show some variability. This can be seen clearly in the volume parameter, where PSRF > 1.1 and the trace of ξ = 0 shows autocorrelation which is much greater than in the ξ ≠ 0 model.

Table 1.

Potential scale reduction factors (PSRF) calculated for two MCMC simulations consisting of 6 independent chains of 3,000 samples each. The processed measurement data was generated with a soil density fluctuation of 300kg.m−3//2 and corresponds with the traces in Figure 4. It has been found that a good indicator of convergence is a PSRF < 1.1.10, 14

Model	ξ	V	σ	ℓx	ℓy	ℓz	x0	y0	z0	ϕ	η
ξ ≠ 0	1.0003	1.0081	1.0002	1.0010	1.0124	1.0033	1.0003	1.0000	1.0114	1.0008	1.0029
ξ = 0	-	1.1261	1.0017	1.0179	1.1621	1.1118	1.0009	1.0001	1.1623	1.0015	1.0781

As described in subsection 3.4, the Bayesian inversion algorithm was tested with ten different realisations of spatially correlated soil noise for d₀ = 50kg.m^−3/2 and d₀ = 300kg.m^−3/2. Table 2 shows the mean PSRF as calculated from the different realisations of soil noise. A number of features reported in Table 1 can also be seen in this table. It is clear that the convergence properties of the ξ ≠ 0 model are generally better compared to ξ = 0. Surprisingly this result is also observed for simulations where no soil noise is present in the simulated data (this data is excluded from this paper for brevity). It appears that including the possibility of spatial correlation between measurements broadens the posterior distribution, making it easier for the MCMC sampler to explore the space irrespective of the fact that the ξ parameter may not required for good agreement with the data.

Table 2.

Mean potential scale reduction factors (PSRF) averaged from 10 realisations of soil density fluctuation d0 = 50kg.m−3/2 and d0 = 300kg.m−3/2. For each set of 10 realisations, inversion models, ξ ≠ 0 and ξ = 0 were tested against the data. It has been found that a good indicator of convergence is a PSRF < 1.1.10, 14

Model	d0 / kg.m−3/2	V	σ	ℓx	ℓy	ℓz	x0	y0	z0	ϕ	η
ξ ≠ 0	300	1.019	1.001	1.005	1.025	1.019	1.001	1.000	1.025	1.001	1.012
ξ ≠ 0	50	1.124	1.002	1.022	1.138	1.100	1.000	1.000	1.136	1.000	1.116
ξ = 0	300	1.277	1.004	1.094	1.318	1.517	1.001	1.002	1.339	1.014	1.330
ξ = 0	50	1.654	1.002	1.197	1.836	1.804	1.000	1.000	1.798	1.001	1.672

Figure 5 and Figure 6 show a selection of histograms, approximate marginal posterior distributions, from the simulation analysed in Table 1 and Figure 4. From the histograms, it is clear that the ability of the inversion algorithm to model spatially correlated noise is important. In contrast to the ξ ≠ 0 model, the ξ = 0 model is over confident and the true values of the parameters are not contained within the 5% and 95% population boundaries. As might be expected, the ξ = 0 model is too simplistic and cannot fit the data. The prior distributions, shown by the green curves in Figure 5 and and Figure 6, do not appear to dominate the marginal posterior distributions: the data is informative. However future work should assess the affect of varying the prior.

Figure 5.

Selected MCMC histograms (approximate marginal distributions) for the inversion algoritm with ξ ≠ 0. The processed measurement data was generated with a soil density fluctuation of 300kg.m^−3/2. For each parameter the mean value is shown by the black line and the dotted lines are the 5% and 95% population boundaries. The true value is shown in red. The green curve is the corresponding prior distribution.

Figure 6.

Selected MCMC histograms (approximate marginal distributions) for the inversion algoritm with ξ = 0. The processed measurement data was generated with a soil density fluctuation of 300kg.m^−3/2. For each parameter the mean value is shown by the black line and the dotted lines are the 5% and 95% population boundaries. The true value is shown in red. The green curve is the corresponding prior distribution.

4.

PROBABILITY OF EXCAVATION HEAT MAPS

Whilst trace and histogram information is useful to understand the performance and convergence of the MCMC sampler, it is not a particularly useful output for an operator of a gravitational imaging system. This is exacerbated when the model dimension becomes large and, as a consequence, review of a few marginal distributions could easily lead to erroneous conclusions. A considerable weakness of inspecting marginal distributions is the inability to easily see correlations between the inferred parameters. A more useful visualisation, probability of excavation (POE) heat maps, was introduced in.^{8, 16}

Using MCMC samples that approximate the posterior probability density, each POE heat map pixel, a_ij, is given by

where m_k is a sample from the simulated Markov chain of length N (this could of course be generalised to voxels in 3 dimensions). The indicator function I_xy returns 1 if the pixel is partially or completely contained by the cuboid as projected into the xy plane (or otherwise) of the 2D heat map. Thus, each sample from the chain describes a cuboid; if every cuboid in the chain overlays a given particle pixel, when projected into the specified plane, then that pixel has a probability of 1 of being part of the void.

Figure 7 shows x-y plane and y-z plane POE heat maps for a single realisation of d₀ = 300kg.m^−3/2 and a single realisation of d₀ = 50kg.m^−3/2. For the anomaly with d₀ = 50kg.m^−3/2, there is no appreciable difference between the quality of the map for the ξ ≠ 0 and ξ = 0 model. This is also despite the fact that PSRF< 1.1 convergence criteria had not been met for the ξ = 0 model. This result was found to be true for all 10 realisations of the case where d₀ = 50kg.m^−3/2. For the case where d₀ = 300kg.m^−3/2, the soil noise does introduce a marked visual difference between the resulting POEs of the two models. It can be seen that the ξ = 0 model under-predicts the overall model uncertainty and generates maps which are essentially overconfident. A statistically significant result is difficult to obtain due to the low number of samples; instead we choose to display the resulting probability of excavation heat maps for the d₀ = 300kg.m^−3/2 realisations in Figures 9 and 10 at the end of this paper. Despite the lack of a comprehensive analysis, empirically it can be seen that the ξ = 0 model choice generates more diffuse probability clouds that are better estimates of the true excavation.

Figure 7.

Probability of excavation in the x-y plane (left) and y-z plane (right). The processed measurement data was generated with a soil density fluctuation of 300kg.m^−3//2 (rows 1 and 3) and 50kg.m^−3/2 (rows 2 and 4). The top 4 images are derived from the posterior distribution that includes inference on spatial correlation due to soil noise, ξ ≠ 0. The bottom 4 images are calculated from the posterior of the model where spatial correlation is neglected, ξ = 0. The projection of the true cuboid void in the plane is identified by the blue rectangle. Sensor measurment locations are on a regular grid identified by points.

Figure 8.

Solid line, autocorrelation function for infinite depth; dotted line, autocorrelation function for 20m depth.

Figure 9.

Probability of excavation in the x-y plane for 8 different soil noise realisations generated with a soil density fluctuation of 300kg.m^−3//2. The images are derived from posterior distributions that includes inference on spatial correlation due to soil noise, ξ ≠ 0. The projection of the true cuboid void in the plane is identified by the blue rectangle. Sensor measurment locations are on a regular grid identified by points.

Figure 10.

Probability of excavation in the x-y plane for 8 different soil noise realisations generated with a soil density fluctuation of 300kg.m^−3/2. The images are derived from posterior distributions that does not include inference on spatial correlation due to soil noise, ξ = 0. The projection of the true cuboid void in the plane is identified by the blue rectangle. Sensor measurment locations are on a regular grid identified by points.

5.

DELTA-CORRELATED SOIL NOISE MODEL

The model used here for the soil density variations is a delta-correlated, Gaussian random variable. This corresponds to white, spectrally flat noise in the spatial-frequency domain. This model assumes uncorrelated density fluctuations at all spatial scale sizes. This is, essentially, the model used for the numerical calculations described earlier (subsection 3.4), save that in the numerical model there is a lower limit to the size of the fluctuations equal to the size of the volume element in the model (a 20cm cube in this case). In fact, as will be shown below, gravity itself acts as a spatial filter, and, consequently, density variations at such small scales have negligible effect on the gravity signal.

Consider a region of ground of volume V divided into N smaller, identical regions with volumes of δV = V/N. Assume the density in these regions deviates randomly from the mean ground density. Let the density deviation be a zero-mean random variable and let the deviations in the different regions be mutually uncorrelated. The total excess mass in the region V is

where d_i is the density deviation in the ith volume. Clearly, the average value of the excess mass is zero (the excess mass may, of course, be positive or negative). An excess density can be defined as D_V = M/V. The variance of the excess density is given by

where the angled brackets denote ensemble averages. The average on the right hand side of Equation 9 is simply the correlation function of d, that is ρ_ij. Using the fact that this correlation function is zero unless i = j, one obtains

The delta-correlated model takes the limit in which the volume δV becomes infinitesimally small. In this case, the summation in Equation 8 can be replaced by an integral and the analogous equation to Equation 9 is

Note that the correlation function of the density variations, ρ_d, has been taken to depend on the difference coordinate only. The correlation function is given by

Here, , is the Dirac delta function, which has the following property

for any well-behaved function f. Thus, with f = 1, Equation 11 reduces to

From this it can be seen that the single parameter that characterises the magnitude of the random density variations, d₀, has units of density times square root of volume, kg.m^−3/2.

Consider a gravity measurement carried out at a point (x,y,0). The z component of the gravitational field produced by random density variations d within a volume V is given by

Clearly, when the density variations, d, are delta-correlated the mean of this randomly varying part of the gravitational field is zero. The autocorrelation function is given by a double volume integral, which can be reduced back to a single volume integral by use of Equation 13. Let the ground surface be in the x-y plane at distance z=z₀ from the measurement point. The integral for the autocorrelation function is then

Note that this is a two dimensional autocorrelation, confined to the z = 0 plane. Call this ρF_z(x, y). The volume V has been taken to be an infinite horizontal slab of thickness z₂ — z₀ and, for convenience, the z₁ axis is taken as pointing downward.

The integrals in Equation 16 can be evaluated by changing the order of integration to do the x₁, y₁ integrals first and noting that these integrals form a 2D convolution which can be calculated by performing a multiplication in the Fourier domain. After carrying out the required forward and inverse Fourier transformations, the final integral with respect to z₁ is an elementary one. The forward and inverse Fourier transforms can be written as Hankel transforms and the required results can be found in standard tables of Hankel transforms.¹⁷ The final result is

In the limit of an infinitely thick slab, this becomes

This is the formula used to describe the spatial correlations in the Bayesian inference. Comparing with Equation 5 it can be seen that .

Also of interest is the power spectral density, which is the Fourier transform of Equation 18 with respect to x and y and is a function of the modulus of the spatial frequency vector only (i.e. it is circularly symmetric).

It is worth examining the ways in which the numerical model differs from the theory. The differences are twofold: first, the volume is of finite extent, whereas Equation 18 is for an infinite half-plane, and, second, there is a minimum size of volume element. The effect of the minimum volume element can be investigated by considering the inverse Fourier transform of Equation 19 when there is an upper cut-off in spatial frequency. Using the fact that the variance is the correlation function with x = y = 0, we can write

where ω_c is the spatial cut-off. This gives

It can be seen that the contribution from spatial frequencies greater than ω_c falls exponentially with ω_c. Consider a spatial scale x_c such that there are two elements of length x_c within a single cycle of the cut-off spatial frequency ω_c, i.e. ω_cx_c = π. If x_c=z₀ the exponent in Equation 21 is equal to −2π and contributions from length scales smaller than x_c contribute less than 1/500 of the variance. In the numerical model, we have x_c = z₀/5 for the volume element of width 0.2m and therefore the effect of the neglected spatial frequencies is completely negligible. It is worth noting that for the modelled volume elements of 20cm×20cm×20cm and d₀ = 300kg.m^−3/2 Equation 14 gives a standard deviation of 3354kg².m^–3 for the density variations, which is actually greater than the mean density. However, because of the exponential cut-off, the resulting non-physical, negative values of density have no effect on the gravitational field variations.

The effect of the finite depth of the modelled region can be investigated using Equation 17. Figure 8 plots the autocorrelation function (normalised by the variance for infinite depth) for both the case of infinite depth and for the 20 metre depth used for the numerical calculations. It can be seen that there is a difference, but it is primarily just a constant offset between the two curves. The offset comes from the contribution of deep density variations which produce gravity fluctuations with characteristic length-scales wider than the 20m region plotted in Figure 8. As such, their effect simply gets included in the background offset parameter, η, so equation Equation 18 can be used for the covariance matrix with negligible error.

The effect of the limited horizontal extent of the volume in the numerical calculation was also investigated and it was found that this only produced errors in areas a metre or so from the edges of the region; however, these area are removed when the original 30m by 30m region is cut down to 20m by 20m.

6.

DISCUSSION

A Bayesian inversion model for the detection of highly localised underground structures at depths not amenable to conventional surveying technology has been presented. As quantum interferometry technology improves the ability to perform rapid gravity surveys, the limiting factor in detection performance will become the ability to intelligently process the gravitational signals. This work therefore represents another step towards developing suitable inversion algorithms.

In particular we have presented a model for spatially correlated soil noise and identified a simple scenario where it is beneficial to account for this spatial variation in the inference model. This can be seen clearly in the probability of excavation maps where the inversion algorithm that neglects spatial correlation calculates a map that is over confident and predicts larger false positive and false negative regions. For this specific scenario it could easily be argued that this error is tolerable, or good enough for practical purposes. However, it is expected that more realistic scenarios will present far more difficult challenges and therefore real situations may occur where the error is much larger. Indeed an important question to be addressed is how large do we expect soil density fluctuations to be in realistic environments? Or perhaps better, to what extent can this simple model of spatial correlation accommodate natural inhomogeneities in the local surface geology? These questions we will seek to address in future work.

The use of a spatially correlated background model can serve as a filter for gravitational anomalies that are below some specified level of interest. For example, there is a relationship between the soil density fluctuation parameter (SDFP) and the expected maximum spherical anomaly volume at some specified depth. Thus a user might calibrate the SDFP to be large enough to accommodate, for example, 1m³ voids at 5m depth. This method may prove advantageous when attempting to infer multiple voids as it may allow the inversion algorithm to avoid generating large numbers of small gravitational features. Alternatively, rather than calibrate the SDFP, we might see larger values of the SDFP preferred to the information theoretic cost of adding new “void particles”.

The utility of inferring the background off-set and the sensor noise have been discussed in previous work.⁸ Inference of the background off-set accommodates unknown systematic shifts due to, for example changes in the water table, variations due to the tides, or indeed soil inhomogeneities on spatial scales larger than the domain of inference. The variable may also be able to accommodate errors such as sensor temperature dependence and linear creep on sensor springs; assuming these parameters do not change significantly during the measurement campaign. Inference on the total noise naturally accommodates many sources of normally distributed noise, without necessarily knowing the origin or magnitude of the source. Thus, provided spatially correlated noise has been accounted for, then an imaging system should be able to resolve structure if the signal strength is comparable to the standard deviation of the inferred noise. It may therefore be possible to calibrate a system by first making measurements of the local environment in an area where it is known that no local anomaly exists. Then, assuming approximate homogeneity of nearby regions, an exploration for anomalies could be usefully undertaken. Comparing the inferred total noise with the previously characterised sensor noise may prove a useful diagnostic for detecting, and potentially correcting, model error.

Perhaps the most surprising result discussed in this paper is the finding that inferring spatial correlation between sensor measurements improves convergence even in the absence of spatially correlated data. A potential interpretation of this is that the MCMC sampler is able to navigate more freely through the posterior space because it is less strongly peaked. For example, changes in depth of the anomaly are less correlated with volume, and therefore the acceptance probabilities of the Metropolic-Hastings accept-reject step are higher. Thus is could be argued that the spatial inference parameter is important purely for computational efficiency. This finding will be explore in greater detail in future work.

With the ultimate aim of developing a gravity imaging system, the key processing challenge to be addressed is accounting for all sources of noise and better understanding the trade-off between number and quality of measurements. Future work will focus on semi-realistic models of the environment including geophysical effects and all other sources of gravitation. Optimised sampling strategies will be investigated: for example using the current probability of excavation to direct sampling to more informative regions of space. The probability of excavation heat maps have been shown as an appropriate way of presenting multi-dimensional information encapsulated within the posterior distribution. It is believed that this form of visualisation could be usefully employed in a gravitational imaging system.