Explosive source signature determination: Two unequal shots in the same hole

Buried explosive charges generate seismic waves with unrivalled bandwidth. The source time function is not always repeatable and is difficult to measure, but is required for subsequent data processing including reverse time migration and full waveform inversion. I present a new method for estimating the explosive seismic source time function for every shot, using two unequal shots in the same shot hole. The measured data for each charge are deconvolved for the modelled monopole seismic source time function to recover estimated earth impulse responses or Green's functions. The source model parameters are adjusted to maximize the likeness of the estimated Green's functions from the two shots, applying physical constraints to obtain the prior probability density functions for each parameter. Ten parameters specify the two source models, but only five of these are independent and a grid search is used to find them. The method is tested on an experimental data set, obtained for a different purpose, consisting of a line of single geophones and three in‐line shot holes, each with two unequal charges, the larger one at twice the depth of the smaller. The method works well where the charge size ratio is at least 2: the estimated source time functions are recovered and the estimated Green's functions are more similar than the seismograms before deconvolution. To apply the method in practice, the smaller charge should be deeper than the larger. If the small charge is small enough, the vertical separation between the charges can be reduced to less than 1 m. The Green's functions for the two charges would then be almost identical for frequencies up to 250 Hz, and the differences in the measured data from the two charges would be caused principally by the differences in the two source time functions.


INTRODUCTION
The source time function of a buried explosive source in a homogeneous elastic medium is a damped sinusoid (Jeffreys, 1931;Sharpe, 1942;Blake, 1952), containing energy down to DC, or 0 Hz. For charge sizes of about 1 kg at depths of the Such a high resonant frequency may be a surprise for geophysicists acquainted with seismic data acquisition using explosives. Lansley (2021b), for example, shows explosive source data in which the resonant frequency of a 1 kg charge at 24 m depth is about 22 Hz. In defence of my statement, I present the data below, in which the source time functions and resonant frequencies of six explosives of different sizes and different depths are derived from seismic data by inversion, using Blake's (1952) explosive source model. I also derive the earth impulse responses, or Green's functions, which all show a decay above about 50 Hz. The reason we do not see the high-frequency resonance in the recorded data is that the data are the result of the convolution of the source time function with the low-pass earth filter.
The bandwidth of the source can exceed the bandwidth of the recording system used in seismic data acquisition. If low frequencies are required, as they are for full waveform inversion (FWI), for example, explosives outperform other sources on land. It is possible to generate frequencies down to 1 Hz with vibrators, at a cost (Lansley, 2021a), but not down to DC. The source time function of explosives has energy down to DC, as can be inferred from Figure 1. FWI and reverse time migration involve forward modelling of the wavefield, starting with the source time function.
It is well known that the explosive source signature is not always repeatable and is difficult to measure (O'Brien, 1957(O'Brien, , 1969Sixta, 1982). The purpose of this paper is to propose a method of seismic data acquisition with explosive sources that enables the source time function of every shot to be determined.
Recently, I proposed a method to invert explosive source data obtained with two explosive charges of unequal size fired separately in essentially the same place, using Blake's (1952) theory to model the source time functions (Ziolkowski, 2021). The earth impulse response from the source to a receiver is the same if the sources are in the same place. In the absence of noise, the spectral ratio of the seismogram of the larger source to the smaller equals the spectral ratio of the larger source time function to the smaller, because the spectrum of the earth impulse response cancels. The inversion problem was to find the parameters of the two sources that minimized the difference between the spectral ratio of the modelled source time functions and the spectral ratio obtained by averaging the spectral ratios of the measured seismogram pairs. The results were reasonable but somewhat disappointing.
There are at least two reasons for this. First, in the available data set, the explosive charges in each pair were not close enough to be regarded as 'in the same place.' This conclusion was obvious and was stated in the paper. Second, no regularization was applied in the inversion. Wang (2017) argues that seismic inversion is inherently ill-posed and requires proper regularization to produce a well-behaved model estimate. Thus, the lack of regularization is a flaw in the paper.
In this paper, I use a Bayesian approach, instead of regularization, applying physical constraints as prior information. Source parameters are estimated using Blake's model, but the method is different. Each shot record is deconvolved for the source time function estimated for that shot, resulting in estimated Green's functions. The differences in the Green's functions for the two shots in the same shot hole are minimized by adjusting the parameters for the two source time functions.
Using explosive seismic exploration data from the 1970s, I first establish that the radius of the nonlinear zone surrounding the explosive source (the elastic radius) is small compared with a wavelength and the source can, therefore, be treated as a monopole. I show that the seismograms of two shots in the same hole are almost identical if the shots are the same size and close enough together. That is, the source time functions are the same and the Green's functions are the same. If the shots are close together, but not the same size, the seismograms are different: the Green's functions are the same, but the source time functions are different. I introduce Blake's (1952) source model and explain the choice of parameters. The data set is the same as that used in Ziolkowski (2021). The method outlined above is then described in detail. The results obtained are better than before, but not as good as we would like, largely because the shots are too far apart. A clear conclusion is reached: The two-shot method should have the smaller charge below the larger, separated by a small distance, of the order of 1 m.

BAYESIAN APPROACH TO INVERSION
'From a statistical point of view, the minimisation problem with a model constraint is equivalent to maximising a conditional probability, given the data measurements. In other words, any probability function that can appropriately describe the behaviour of the expected model can be adopted as a regularisation in the inversion' (Wang, 2017, p. 68). This is the approach I adopt here, using a Bayesian formulation (Tarantola, 2005;Sen & Stoffa, 2013).
Bayes' rule is where, in this problem, the data are the two shot records, the elements of are the model parameters of the two sources, ( | ) is the conditional probability density function (pdf) of the model, given the data, ( | ) is the conditional pdf of the data, for a given model , ( ) is the pdf of the model, and ( ) is the pdf of the data. The denominator ( ) is independent of the model and is usually regarded as impossible to compute for complex problems (e.g. Lambert, 2018, p. 115). Sen and Stoffa (2013, p. 63) say that, in geophysical problems, ( ) is a constant and equation (1) then becomes The constant of proportionality is unknown. The ( | ) is still a probability density function, however, and is known as the posterior probability density function (PPD; Sen & Stoffa, 2013). When the measured data obs are substituted into the expression for the conditional pdf ( | ), it is called the likelihood function ( obs | ) (Sen & Stoffa, 2013, p. 63), and Equation (2) becomes The goal is then to find the model parameters that maximize the PPD ( | ).
In this formulation, I use Blake's (1952) theory to model the two source time functions, as before (Ziolkowski, 2021). The rest of the scheme is different. The measured seismograms for each shot are deconvolved for the corresponding modelled source signature to recover the estimated earth impulse responses or Green's functions. The normalized root-mean square difference (NRMS) (Kragh & Christie, 2002) between the Green's functions of the two shot records is a measure of the error. The reciprocal of NRMS is taken as the likelihood function. The prior probability of the model p(m) is the product of the probability density functions of the independent model parameters m, derived from physical constraints, and is independent of the data. The PPD of the model, given the data, is the product of the likelihood function and the prior probability of the model. The goal is to adjust the model parameters to maximize the PPD. Figure 2, first used in Ziolkowski and Lerwill (1979), shows the result of firing three separate shots in a single hole, starting with the deepest one and progressing upwards. The first two shots were each a single detonator; the third shot was 120 g of explosive plus a detonator. Detonators vary in size, but normally contain not more than 10 g of explosive (SLP-12, 2018). The second shot record is almost identical to the first. It follows that the influence of the first shot on the second is negligible and the nonlinear zone surrounding the first shot must be small, say less than 30 cm radius and small compared with a wavelength. It can, therefore, be modelled as a monopole. The third record is different because the charge size is different. The increase in charge size shifts the bandwidth to lower frequencies. The path from each source to a given receiver is almost the same, so the earth impulse response from each source to the receiver is almost the same. The appearance of the seismograms depends on charge size, which affects only the source time function. Figure 3 shows the model for Blake's (1952) spherical monopole source, of radius , embedded in a homogeneous isotropic full space. The pressure on the inside of the sphere is ( ), where is time. Blake (1952) solves the problem for ( ) = 0 ( ), with 0 a constant and ( ) the Heaviside unit step function. Any arbitrary pressure function can then be incorporated in the solution using Duhamel's integral.

THE MONOPOLE EXPLOSIVE SOURCE
The origin is at the centre of the sphere and the motion is radial. Following the argument in Ziolkowski (2021), the displacement as a function of radial distance and time is where the prime denotes differentiation with respect to the argument, and the compressional wave speed F I G U R E 2 Left: the setup (not to scale), showing the three charges in the hole with 60 cm separation between charges; right: the three shot records, with 1 ms sampling, the order of firing, bottom-to-top, corresponding to records left-to-right. Blake's (1952) model, in which an internal pressure ( ) is applied uniformly over the interior surface of a cavity of radius , buried in an infinite homogeneous isotropic elastic medium with density and Lamé elastic constants and . Equation (4) may be written as

F I G U R E 3
where ( ) is the Dirac impulse function, the asterisk denotes convolution and the function ( ) has dimensions of volume and is known as the reduced displacement potential (Werth et al., 1962) and is reduced, because it is a function of time only. The factor in square brackets is the Green's function for a whole space and consists of two terms. The 1∕ 2 term dominates in the near-field at distances small compared with a wavelength, and the 1∕ term dominates in the far-field at distances large compared with a wavelength. In a heterogeneous medium, the Green's function is more complicated. We assume that the scattering from heterogeneities does not affect the behaviour of the source, which depends, therefore, on the medium in the immediate neighbourhood. The particle velocity is the time derivative of the displacement and may be written as where the source time function ′ ( ) has dimensions of volume divided by time and is known as the reduced velocity potential. For explosive source data measured with particlevelocity sensors, or geophones, this is the source time function we need. Blake (1952) derived the reduced displacement potential ( ), which can be differentiated to yield the reduced velocity potential: where = (1∕2)(1 − )∕(1 − 2 ), Poisson's ratio = (1∕2) ∕( + ), the damping factor 0 = ∕2 and the resonant angular frequency 0 = 0 (4 − 1) 1 2 . The Fourier transform of Equation (9) is (Ziolkowski, 2021)

The new idea
We consider two unequal explosive sources in the same shot hole, with reduced velocity potential source time function ( ) for the small source and ( ) for the big source. There are receivers. The velocity seismograms at the kth receiver are for the small source, where ( , ) is the Green's function from the small source to the kth receiver and ( , ) is the noise, and for the big source, where ( , ) is the Green's function from the big source to the kth receiver and ( , ) is the noise. If the distance between the sources is small compared with a wavelength, the Green's functions should be nearly identical: In the frequency domain, Equations (11) and (12) may be written aŝ( If we have estimates of the source spectra,̂( ) for the small source and̂( ) for the big source, we can estimate the Fourier transforms of the Green's functions aŝ . (17) The goal is to choose the parameters for the two source time functions so as to maximize the likeness of the recovered Green's functions ( , )and ( , ) for all receivers. The two Green's functions for the kth receiver include the receiver impulse response and ground coupling.
To quantify the likeness of two traces and , Kragh and Christie (2002) use the normalized root-mean square difference (NRMS) as a metric, expressed as where the rms operator is defined as and is the number of samples in the interval 2 − 1 . Kragh and Christie (2002) state that the values of nrms are non-intuitive and give examples. 200% is the theoretical maximum. For the kth receiver, the NRMS obtained using the recovered Green's functions ( , )and ( , ) is NRMS( ). NRMS increases with the difference and is zero if the traces are identical. As mentioned in the Introduction, we choose the reciprocal of NRMS as the likelihood function. For receivers this is where the data d are the pairs of recovered Green's functions for the two shot records, and m are the model parameters, defined below.
The data set Figure 4 shows part of an experiment to measure seismograms from explosive charges of different sizes and depths. The geophone spread consisted of single geophones with geophones 1 to 48 at 2.54 m (8.33 ft) intervals and then geophones 49 to 96 at 3.81 m (12.5 ft) intervals. The shot holes were at 7.62 m (25 ft) intervals. Figure 4 shows the first three shot holes at shot stations 501, 502 and 503, coincident with geophone station numbers 1, 4 and 7. Each of the shot holes contained two charges, with parameters given in Figure 4. In each case, the deeper shot was fired first. Figure 5 shows the three pairs of shot records, with both time-dependent and offset-dependent gain applied. The deeper shot is to the right for each pair. There are significant differences in data quality between the shallower and deeper shots. The shallower shots generate relatively more surfacewave energy and less reflection energy. It is clear that the earth responses for two shots in the same hole are not identical.

Preprocessing
Three preprocessing steps were applied to the data before the application of the method: geophone response removal, source-to-receiver distance correction and first-break picking correction. The geophone response removal and first-break picking are described in Ziolkowski (2021).
The source-to-receiver correction is a first-order correction for the difference in distance to a receiver between the deeper and shallower shots in each pair. Figure 6 shows the geometry.
The Fourier transform of the full-space Green's function of Equations (7) and (8) F I G U R E 5 Vertical pairs of shot records with the shallower shot on the left and the deeper on the right. Vertical axis is seconds.

F I G U R E 6
Showing source-receiver geometry in the ( , ) plane: small source at ( , ), big source at ( , ) and receiver at ( , 0). This is in the radial direction. Approximate correction of the amplitude and phase of the spectrum of the big source, to account for the extra travel distance of the vertical component of the direct wave to the receiver compared with the vertical component of the direct wave of the small source, may be performed by a multiplication of where = ( 2 + ( − ) 2 ) 1 2 , = ( 2 + ( − ) 2 ) 1 2 and and are the elastic radii of the small source and big source, respectively. We do not know the exact values of and . They can be estimated, as described below. This correction is important at small offsets, where the direct wave arrives first. At larger offsets, where refracted arrivals arrive first, I make no correction. For these data, refracted arrivals arrive first at offsets greater than about 20 m.
Because the correction is not exact, it is still necessary to line up the first breaks of the big source to be equal to those of the small source.

Estimation of source model parameters
The model parameters for each source are density , internal pressure 0 , P-wave speed , elastic radius and Poisson's ratio . It is seen from Equations (9) and (10) that density and internal pressure 0 are scalar quantities that affect only the amplitude and not the variation with time or frequency. The density is known approximately and is expected to be about the same, 2000 kg m −3 , for both shots in the same hole.
The internal pressure is, in principle, not known. If the rock has negligible tensile strength, it fails in tension as the gaseous products of the explosive expand, forming radial tensile fractures. This has important consequences for the estimates of internal pressure and elastic radius. The internal pressure must then be equal to the overburden pressure. That is, 0 = + atm for the small source, where is gravity, 9.81 m s −2 , and atm is atmospheric pressure, 10 5 Pa. Then 0 = 0 for the big source, with the parameter to be found. If the rock has negligible tensile strength, we expect the value of to be For each shot hole in this data set, the big source is at a greater depth than the small source, so > 1. The absolute scale of the source time functions is not known. Our source amplitudes are underestimates, because the internal pressures are the minimum possible, if the density is correct.
The remaining parameters are the P-wave speed , elastic radius , and Poisson's ratio , for each source. The P-wave speed can be estimated from the data and was estimated to be 2018 m s −1 from uphole times for this data set (Ziolkowski, 2021). Of course, the P-wave speed varies slightly from place to place, both laterally and vertically, and we do not know what this variation is. In Equations (9) and (10), the P-wave speed occurs with the elastic radius as the ratio ∕ in both the resonant angular frequency 0 and the damping factor 0 . Thus, any error in P-wave speed is transferred to a corresponding error in elastic radius . The two elastic radii, for the small source and for the big source, and the two Poisson's ratios, and , are independent parameters. We do not know the values of Poisson's ratio and , but do know that 0 < < 0.5. At the clay site, O'Brien found Poisson's ratio to be 0.47 (p. 515).
To obtain an initial estimate of the elastic radius, we can start with measurements made by O'Brien (1969). Using a hydrophone in a borehole at a distance of about 30 m from an 8.5 kg explosive shot at a depth of 60 m in wet clay, O'Brien (1969) found the elastic radius to be about 3.7 m, with the tangential stress 'commensurate with the overburden pressure' (p. 535). He wrote, 'We may imagine that the surface of the equivalent radiator divides the rock into an inner part which acts as a fluid and an outer part which acts as a solid' (p. 536). That is, the pressure inside the elastic sphere is uniform and equal to hydrostatic pressure, or + atm . O'Brien found a density of 2000 kg m −3 at the wet clay site. Therefore, the internal pressure at 60 m depth was 2000 × 9.81 × 60 + 10 5 Pa, or about 12.8 × 10 5 Pa. The gas produced by the explosion expands very quickly, and therefore adiabatically, to reach hydrostatic pressure. Assuming the gas to be ideal, the pressure-volume relation for a given mass of gas is given by the well-known equation where is the ratio of the specific heats. The detonation products of mining explosives, including dynamite, are predominantly carbon dioxide with small amounts of carbon monoxide and nitrogen oxides (Zawadzka-Małota, 2015). Some steam is also produced, with a specific heat ratio of 1.3. The ratio of specific heats of carbon dioxide is 1.3022 at 290 K (Partington, 1921). We assume = 1.3. The mass of gas produced is equal to the charge mass. For a given rock and type of explosive, the elastic radius scales as the cube root of the charge mass (O'Brien, 1957;Werth & Herbst, 1963;O'Brien, 1969;Ziolkowski, 1993). For a rock with negligible tensile strength, the elastic radius of an explosive source of mass at depth may, therefore, be estimated as where 0 is the elastic radius of a charge of mass 0 at a depth 0 . We do not know the rock type for the Texas data set. Using O'Briens's values of 0 = 8.5 kg, 0 = 60 m and 0 = 3.7 m in wet clay, we may apply Equation (25) to the data of Figure 3, to obtain the estimates of elastic radius presented in Table 1. This is an approximation, since the Texas rock may not be wet clay, but gives an indication of the relative sizes of the elastic radii of the sources and helps to define the starting model.
In summary, we search for five independent parameters: , , , and , for the pair of shots in each borehole. These are the elements of the model vector = [ , , , , ].

Prior probabilities
For each model parameter, , say, there is a physically possible range, min to max , which can be divided into a number of equal increments Δ , such that the kth value in the range is = min + ( − 1)Δ . If we know only a reasonable estimate of the range and nothing about the probability of any value within the range, we assume a uniform distribution. The probability of any value in the range is then ( ) = Δ ∕( max − min ).
There is considerable uncertainty about the values of Poisson's ratio, but we need constraints to obtain physically reasonable answers. Poisson's ratio for soils and rocks must lie in the range 0 < < 0.5. Gercek (2006, Figure 4) gives typical ranges of values for Poisson's ratio for some rock types, from 0.05 to 0.4. A Poisson's ratio of 0.25 is possible for all rocks in Gercek's Figure 4. If = 0.5, the shear modulus is zero and the medium is fluid, not solid. We know, however, that Poisson's ratios approaching 0.5 are possible (e.g. O'Brien, 1969). If = 0, it means a cylinder compressed axially does not expand radially. Rocks do not behave like this.
To exclude values of = 0 and = 0.5, we adopt the beta probability distribution, scaled to have bounds 0 and 0.5, as shown in Figure 7, with, in this case, parameters 1.1 and 1.1. The parameters are arbitrary, the bounds are not. I have found that physically unrealistic answers are possible if Poisson's ratio is not constrained slightly. The beta distribution shown in Figure 7 is zero at = 0 and = 0.5 and almost flat in the middle of the range, with a maximum at = 0.25. beta(1.1,1.1). The area under the curve is 1.

Grid search
The five model parameters = [ , , , , ] are varied within nested loops, each with its own range and increment. For each combination of parameters, the source time functions ( ) and ( ) are calculated and the shot records of the corresponding shots are deconvolved in the frequency domain, according to Equations (16) and (17) Green's functions ( , ) and ( , ) and the NRMS of the estimated Green's functions is calculated for that combination. We choose the reciprocal of NRMS as the likelihood, as defined in Equation (20). Then the posterior probability distribution of this particular model, given the data, is We search for the combination of parameters that maximizes ( | ).

RESULTS
The results of inversion of shots F630 and F631 are shown in Figure 8. Figure 8a shows the recovered reduced velocity potential source time functions at 2 ms sampling; Figure 8b shows the result of deconvolving a random trace (60) of shot records F630 and F631 to recover the estimated Green's function; Figure 8c compares the mean amplitude spectra of the shot records; Figure 8d compares the mean amplitude spectra of the deconvolved records. To obtain the mean spectra, each trace was scaled to have a standard deviation equal to the mean standard deviation of the record; the mean amplitude spectrum of all traces in the record was then calculated.
The results of inversion of shots F628 and F629 are shown in Figure 9 and the results of inversion of shots F626 and F627 are shown in Figure 10. Table 2 presents the source parameters of all six shots found by inversion.

DISCUSSION
The idea under test in this paper is that the source time functions of two unequal explosive charges, fired separately in the same shot hole, can be obtained from the seismic data measured with the same receivers. The charges must be close enough together for the Green's functions of the two shots to be essentially the same. They must also be unequal, such that their natural frequencies are different. From Figure 4 it is clear that the shots in each hole in this data set are not really close enough together to satisfy the requirements of the method: the shallower shots produce proportionately more surface-wave energy and less reflection energy. Nevertheless, reasonable results are obtained for two pairs of shots: pair 630 and 631, and pair 626 and 627.
Pair 630 and 631. From Figure 8b and 8c, it is clear that the larger, deeper 7.5 kg shot generates higher amplitudes and has proportionately more low-frequency energy than the shallower 2.5 kg shot. It is no surprise, therefore, that the recovered source time function of the larger source has a higher amplitude and a longer period, as shown in Figure 8a. Deconvolution of the shot records for the respective source time functions recovers estimated Green's functions that are more similar than the input data, as shown in Figure 8b. Figure 8d shows that the average spectra of the two sets of estimated Green's functions are very close, except at frequencies well above 100 Hz.
Pair 626 and 627. This is a similar story, as shown in Figure 10, with the estimated Green's functions being very similar above about 10 Hz, but less similar at lower frequencies. Again, the larger charge generates higher amplitudes and has proportionately more low-frequency energy, reflected in the recovered source time functions, shown in Figure 10a.
Pair 628 and 629. Here the method clearly does not work. The ratio of the charge sizes is 2.5:1.5. If they were at the same depth, we would expect cube-root scaling of the source time functions (e.g. Werth & Herbst, 1963;O'Brien, 1969;Ziolkowski, 1993). Because the larger charge is deeper, however, the overburden pressure is higher and the recovered elastic radius is actually smaller than for the smaller charge. The recovered source time functions are almost equal, apart from an amplitude scaling factor and the recovered natural frequency of the larger charge is actually higher (240 Hz) than for the smaller charge (216 Hz), as given in Table 2. Putting the larger charge deeper tends to equalize the natural frequencies of the sources.
One reviewer has suggested that the proximity of the free surface may have an effect on the recovered source time function, especially for the shallowest shot, 629, which is 9 m deep. The recovered elastic radius is 1.4 m (see Table 2). The reflected wave, therefore, travels a distance of 15.2 m, or roughly 11 elastic radii and appears to come from a virtual source 9 m above the surface, delayed by about 7.5 ms. For the lowest frequencies (< 2 Hz) this is a near-field or 1∕ 2 effect and is less than 1 % of the emitted amplitude. At the highest frequencies (> 200 Hz) this is a 1∕ effect and is about 10% of the emitted amplitude. For most of the bandwidth, the effect is between 1% and 10%. The interaction with the sur-face is not a major issue for these data, but it would be for larger charges at shallower depths and must certainly be the case for underground nuclear tests, in which the devices are normally placed just deep enough to contain the products of the explosion, although it is not completely obvious what this depth should be (Carrigan et al., 2016).
Comparing the recovered elastic radii of Table 2 with the estimated elastic radii of Table 1, we see an approximate correspondence within about a factor of 2. We should bear in mind two considerations. First, the estimates of Table 1 were based on results from O'Brien (1969) in wet clay in the UK and were applied to the experimental data in unknown soil or rock in Texas. Second, because the sources in the same hole were not close enough together, the algorithm compensates for differences in the Green's functions as well as differences in the source time functions. The application of the scaled beta function probability density function for the prior probability of Poisson's ratio is an important physical constraint on the resulting solution. Where the seismic data generated by the two shots in the same hole are sufficiently different (pair 630 and 631 and pair 626 and 627), the method produces reasonable results and the deconvolved data are much more similar than the raw data.
To put the idea into practice, the small charge should be deeper than the large charge and fired first, of course. The data of Figure 1 show that the vertical separation between the charges can be small compared with a wavelength if the small charge is small enough. Then the Green's functions for the two charges are almost identical and the differences in the data from the two charges would be caused only by differences in the source time functions, and not also by differences in the Green's functions. A vertical separation of 1 m between the charges would create a negligible difference between the Green's functions for frequencies up to 250 Hz.

CONCLUSIONS
I present a new method for estimating the source time functions of buried explosive source seismic data. Two unequal charges are fired separately in the same shot hole, the resulting wavefields being recorded by an array of seismic receivers. The measured data for each charge are deconvolved for the modelled source signature to recover estimated earth impulse responses or Green's functions. The parameters of the modelled monopole source time functions are adjusted to maximize the likeness of the estimated Green's functions, applying physical constraints to obtain prior probability density functions for each parameter. Ten parameters specify the models of the two source monopoles, but only five parameters are independent and a grid search may be used to find them.
The method is tested on an experimental data set obtained for a different purpose. Since soils and rocks are weaker in tension than in compression, the expanding explosive gases must cause radial tensile fractures to form in the nonlinear zone around the explosive cavity. The pressure at the elastic radius is therefore assumed to be equal to the overburden pressure as a guide to start the model. Since the gases expand adiabatically, the resulting gas volume is related directly to the pressure and proportional to the mass of explosive. The volume of the monopole sources is assumed to be proportional to the gas volume to give a starting estimate of the elastic radius. The method works well where the charge size ratio is at least 2: The estimated source time functions are recovered and the estimated Green's functions are more similar than the seismograms before deconvolution.
To apply the method in practice, the smaller charge should be deeper than the larger and fired first. The vertical separation between the charges can be small compared with a wavelength, perhaps about 1 m, if the small charge is small enough. The Green's functions for the two charges would be almost identical up to 250 Hz and the differences in the data from the two charges is then caused principally by differences in the source time functions.

A C K N O W L E D G E M E N T S
I am very grateful to Bruce Hobbs, Guido Baeten, Paul Stoffa and Andrew Curtis for reading the paper and for constructive suggestions. I thank the three reviewers, Malcolm Lansley, Peter Pecholcs and Gabrielle Busanello for very helpful comments and suggestions and editors Claudio Bagaini and Tijmen Jan Moser for a very speedy turnaround of my paper. I am also very grateful to Shell for providing the test data set.

D A T A AVA I L A B I L I T Y S T A T E M E N T
The data used in this paper do not belong to the author and are not shared.