# Spectral and wavelet analysis of gilgai patterns from air photography.

IntroductionGilgais are essentially Australian phenomena; nowhere else in the world, it seems, are they as common and as well developed. They fascinated Australian pedologists for decades. Studies on their occurrence, physical chemistry, and mechanical behaviour were reported in the international literature. Their topography and their various repeating patterns also attracted attention and were described in qualitative terms (see, for example, Hallsworth et al. 1955). By the 1970s new computing power and developments in spatial statistics made it possible to analyse the spatial distribution, at least in 1 dimension. Russell and Moore (1972) estimated parameters of gilgai microrelief by Fourier analysis of transect data. Webster (1977) took the analysis a stage further with soil data from a 1-km transect across gilgai-patterned ground on the Bland Plain of New South Wales. He had recorded the height of the soil surface above a local datum and the morphology, pH, and electrical conductivity of the soil down to 1 m on cores at 4-m intervals. All of the properties showed strong spatial dependence at that scale. More interesting, however, was that the correlograms of several fluctuated in an apparently periodic manner and that their Fourier transforms, the power spectra, had peaks at a frequency corresponding approximately to the distances between the centres of neighbouring gilgais, ~34m. The results suggested that there was a degree of regularity in the spatial pattern. Webster and Oliver (2001) came to the same view after modelling the variogram of the soil's electrical conductivity, EC. Lark (2005) analysed the categorical variables on morphology that Webster (1977) had collected on his transect. There were 3 classes: 'puff', 'depression', and 'plain'. Geostatistical analysis showed that the classes 'plain' and 'depression' varied periodically with wavelengths of 36 and 33 m, respectively.

The patterns are 2-dimensional, of course, and the results of the 1-dimensional analysis led to the question: is the 2-dimensional arrangement of the gilgais regular? Some of the photographs in the paper by Hallsworth et al. (1955) of gilgai stripes on sloping land suggested periodic repetition across the stripes but not in the perpendicular direction. On aerial photographs of fiat plains, the gilgais appear as isometric patches, in some instances roughly circular, of similar size and with similar distances between neighbours--see figure 1 in Webster (1977). One could imagine that the gilgais represent distortions of a hexagonal close packing arrangement.

[FIGURE 1 OMITTED]

Measuring properties of the soil, such as its EC, at enough sites for a 2-dimensional spectral analysis was prohibitively expensive. An alternative was to analyse the photographic images. At the time, the 1970s, however, the available digitising equipment proved too unreliable to provide sufficient surrogate data. An alternative was to process the images optically, a technique then being explored by Preston and Davis (1976) for analysing photomicrographs of sedimentary rocks.

Transparencies were made of several aerial photographs, including those from which Fig. 1 is derived; they were placed on a light bench and the images converted to their Fourier transforms in this way. The results were in all instances the same. The transforms were dominated by bright centres away from which the light gradually diminished towards the peripheries; there were no evident subsidiary peaks to indicate periodicity.

Modem digitising equipment hugely increased computing power, and new mathematical developments such as wavelet analysis have enabled us to revisit the question. Here we describe what we have done and our results.

The Bland Plain and its photography

The Bland Plain lies in central New South Wales. It is traversed by the Bland Creek, which flows north into Lake Cowal and from there to the Lachlan River. The Plain is virtually level, and in years of exceptionally heavy rain it is deeply flooded. Its soil, formed in alluvium, is dominantly clay, alkaline, and more or less saline. At depth it is calcareous and contains gypsum. The topsoil is somewhat more sandy, but the texture contrast is typically gradual rather than abrupt. Locally the Plain is dotted with gilgais, mainly of the 'crab-hole' kind (McKenzie et al. 2004, page 48), in distinct patterns including ones like that illustrated on page 371 of McKenzie et al. (2004). The soil in the gilgais themselves is cracking clay, Vertosol in the Australian system of soil classification.

The Plain was photographed at a scale of ~1 : 40 000 by the New South Wales Department of Lands and Surveys in the late 1960s. Webster (1977) selected individual photographs showing patterned ground, had them enlarged to ~1:20 500, and then visited the ground to confirm that the patterns were indeed of gilgai. These photographs provide the images for our new analysis. They include the photograph of Caragabal station, which Webster (1977) surveyed and sampled for his original study.

For present purposes we have chosen 2 photographs displaying gilgai patterns, one of Caragabal and the other immediately to the west of Back Creek lying some 23 km east of West Wyalong (between Caragabal and West Wyalong). On each we selected a rectangular patch of side 30-40 mm, corresponding to ~25 ha on the ground, and scanned it electronically at a resolution of 400 pixels per inch, equivalent to 15.7 pixels per mm or ~1.3 m as the diameter of a pixel. The optical density of each pixel was recorded as a grey level in the range 0-255. Figure 1 shows the 2 rectangles as digital images, and Table 1 summarises their statistics.

Correlograms and spectral analysis

We computed variograms and autocorrelograms from the pixel data both for individual rows and columns in the images and then for the full 2 dimensions by the usual method of moments. Thus, for the rows and columns of length N, we computed semivariances for lags h = 1,2,...:

[??](h) = 1/2N[N-h.summation over (i=1][{[Z.sub.i] - [Z.sub.i+h]}.sup.2], (1)

where [z.sub.i] and [z.sub.i+h] are the optical densities of the/th and i+hth pixels in the row or column. The corresponding autocorrelations are

[??](h) = 1 - [??](h)/[s.sup.2], (2)

where [s.sup.2] is the variance of the data. We computed the 2-dimensional variograms on grids of size [N.sub.x] x [N.sub.y] by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (3)

where p and q are, respectively, the lags along the rows and down the columns of the grid. These equations enable half the variogram to be computed for lags from -q to q to and from 0 to p. The variogram is symmetrical about its centre, and the full set of semivariances is obtained as:

[??](-p, q) = [??](p, -q) and [??](-p, -q) = [??](p, q).

Following convention we have computed the autocorrelation using the biased estimator for 1 dimension. For 2 dimensions we retain the unbiased formula. Again, the autocorrelations are obtained as the complements of [??](.)/[s.sup.2] as in Eqn (2).

The 1-dimensional spectra for individual rows and columns of the images were computed as the Fourier transforms of the autocorrelograms (see Chatfield 1984):

[??](f) = 1 + 1 [K.summation over (k=1)][??](k)w(k)cos(2[pi]fk), (4)

for frequency f in the range 0 to 1/2 cycle and lags k from 1 to a maximum of K. This maximum effectively defines the width of a window in the lag domain through which the correlogram is viewed; the narrower it is, the more the spectrum is smoothed. The quantity w(k) defines the shape of the window. We chose the Parzen window because the 'leakage' is small. It is given by

w(k) = l-6 [(k/K).sup.2] + 6 [([absolute value of k]/K).sup.2] for 0 [less than or equal to] [absolute value of k] [less than or equal to] K/2

= 2[(1-[absolute value of k]/K).sup.2] for K/2 [less than or equal to] [absolute value of k] [less than or equal to] K

= 0 for [absolute value of k] > K. (5)

We computed the 90% confidence interval of the estimate for the 1-dimensional spectra using the formula given in Priestley (1981).

For 2 dimensions we computed

[??](u, v) = [K.summation over (q = -K)] [K.summation over (p = -K)] [??](p, q)w([square root of [p.sup.2] + [q.sup.2]])cos{2[pi](up + vq)}, (6)

where u and v are, respectively, the frequencies along the rows and down the columns of the correlogram. The argument of the window, [square root of [p.sup.2] + [q.sup.2]], is simply the distance between the origin of the correlogram and the point (p, q).

Results

One-dimensional spectra

We computed variograms, correlograms, and spectra for numerous individual rows and columns from the images using the formulae above. Some of the variograms were distinctly wavy and some were less so. The spectra of the former had pronounced peaks, especially when computed with wide lag windows (110 pixels). Figures 2 and 3 show examples; Fig. 2a shows the wavy correlogram for row 400 of the Caragabal image, and in Fig. 2b there is a marked peak in the corresponding spectrum. The peak is at [approximately equal to] 0.04 cycles, equivalent to a wavelength of 25 pixels. The digitising resolution (see above) was 157 pixels per cm and the scale 1 :20 500, so the wavelength on the ground is approximately

[lambda] [approximately equal to] 25 x 20 500/157 x 100 m [approximately equal to] 32 m.

In Table 2 we have tabulated some values which occur in the paper in terms of frequency (cycles/pixel), wavelength (pixels), and wavelength (m) to help the reader. Row 400 from Back Creek is another example with a strong peak in its spectrum (Fig. 3b). This peak is at [approximately equal to] 0.025 cycles, corresponding to a wavelength of [approximately equal to] 52 m.

Column 300 from the Caragabal image and column 400 from the image at Back Creek have more erratic correlograms, and less pronounced peaks--see Figs 2d and 3d. Of the 16 one-dimensional spectra, 8 from each photograph, approximately half, had pronounced peaks in the range 0.025-0.05 cycles and more or less wavy correlograms.

To assess visually the significance of the spectral density peaks observed in Figs 2 and 3 we calculated the 90% confidence intervals for each peak and plotted them on the graphs. In the 2 cases with pronounced peaks (Figs 2b and 3b), the confidence intervals indicate that it is unlikely that the peaks have arisen by chance. The evidence is weaker for the 2 less pronounced peaks (Figs 2d and 3d).

Two-dimensional analyses

Figures 4 and 5 show the correlograms computed with a lag window of 110, corresponding to a bandwidth of [approximately equal to] 0.017 cycles per pixel, for the 2 images. Both correlograms decay smoothly and approximately exponentially in all directions but with perceptible waves on them. The corresponding spectra (Figs 6 and 7) have sharp peaks at the origins and subsidiary rings of large spectral density at 0.025-0.04 cycles per pixel, suggesting periodic fluctuation at those frequencies.

Wavelet analysis

One-dimensional wavelet transforms

Spectral analysis as described above treats data such as the grey levels in the gilgai images as the outcomes of stationary random processes, that is, stationary in both the mean and the variance. It helps us to identify periodicity and to estimate the frequency or frequencies of the periodicity. A shortcoming is that it loses all information on the positions of features in the images. Wavelet analysis can make good this deficiency by decomposing the data into separate components of both frequency (or scale) and location. Further, it requires no assumption of stationarity; instead, it can reveal where the characters of an image or series of data change.

[FIGURE 2 OMITTED]

Wavelets are relatively new. They were devised in the late 1980s for processing signals from various sources and were introduced to soil science by Lark and Webster (1999). Since then they have been applied in several pedological case studies by Lark et al. (2004a, 2004b), Milne et al. (2005), Lark (2006), and Biswas et al. (2008); Percival and Walden (2000) have written an excellent text book on the subject.

[FIGURE 3 OMITTED]

Equation 1 shows that the autocorrelogram is computed from comparisons over the whole length of a row or column of the data, apart from near the ends. The spectrum, its Fourier transform, is similarly a function of all the data therefore. In contrast, wavelet transform functions are non-zero for only a local window (they have compact support). So the coefficient obtained for a particular wavelet function corresponds to the variation of the data within that restricted window. The set of wavelet functions

[[PSI].sub.[lambda],x] = 1/[square root of [lambda]] [PSI](u - x /[lambda]), [lambda] > 0, (7)

is obtained from what is known as the mother wavelet, [psi](x). In this equation, [lambda], is a scale parameter that controls the width over which the wavelet takes non-zero values, x is the position at which the wavelet is centred, and u represents the displacement from that central position. By integrating the variable of interest with the wavelet function we obtain wavelet coefficients:

W([lambda],x] = [[integral].sup.[infinity].sub.-[infinity]] z(u) 1/[square root of [lambda]] [PSI] (U - X/[lambda]) du. (8)

[FIGURE 4 OMITTED]

Varying x moves the wavelet over the space occupied by the data. More interesting in the present context is what is achieved by changing [lambda]; the smaller it is, the finer is the scale at which it describes the variation about x; increasing it dilates the wavelet and describes variation about x at an increasingly coarse scale. It remains then to choose values for [lambda].

Lark and Webster (1999) described the discrete wavelet transform (DWT) whereby [lambda] is incremented in a series of powers of 2: [lambda]=[2.supj], j=1, 2, ..., J, to give a series of orthogonal sequences with locations moved in steps of [2.sup.j]. There is a restriction that the data must be an integer power of 2 in length, and that [2.sup.J] [less than or equal to] N. In practice, the wavelet coefficients for a particular scale are calculated by convolution of the data with a scale-specific filter [a.sub.j] (u) of length [L.sub.j], where u= l, 2, ..., [L.sub.j] The output is sampled at intervals of [2.sup.j] to give the N/[2.sup.j] wavelet coefficients. Scaling coefficients, which capture the trend in the data, result from filtering them with the filter [b.sub.J](u). If N= [2.sup.J] then we obtain a single scaling coefficient, which is the mean of the data.

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

[FIGURE 8 OMITTED]

Wavelet variance

Lark and Webster (1999) gave the formula for the sample variance, [[??].sup.2.sub.j], associated with each scale interval [2.sup.j]l to [2.sup.j+1]l, where l is the interval between observations. Percival and Guttorp (1994) had already explored an alternative approach to the estimation of wavelet variances, which was adopted by Lark and Webster (2001). In this technique, all N coefficients from the wavelet convolution are retained. In other words, the wavelet coefficients at scale j and locations x = 1, 2, ..., N are given by

[d.sub.j](x) = [[L.sub.j].summation over (i=1)][a.suib.j](i)z(x - i) [mod N]. (9)

[FIGURE 9 OMITTED]

The notation [mod N] implies circular filtering, i.e. that z(N+K)=z(k). Percival and Guttorp called it the maximal overlap discrete wavelet transform, MODWT. The resulting sequences are no longer orthogonal to one another, but there are 2 advantages of MODWT: the number of data, N, need no longer be an integer power of 2; and we achieve more efficient estimation of the wavelet variance, which is given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (10)

where [n.sub.j], which may be less than N, is the number of data involved in the computation. We note that in examples such as the one discussed here, the circular filtering does not make any physical sense. Therefore, in cases where the filter overlaps the ends of the data, we pad the data by reflection.

[FIGURE 10 OMITTED]

The wavelet coefficients capture local components of variance at specific scales, and so not only can we see how variance is partitioned from scale to scale, but we can see also how it changes with location. Lark and Webster (2001) describe a method of detecting changes in wavelet variance at a given scale based on a proposal by Whitcher et al. (2000), and testing to see if the changes are significant based on a null hypothesis that the data are a realisation of a random variable that is stationary in the variance.

Multi-resolution analysis

The data contain components of variation at several frequencies, and these components can be difficult to identify by visual inspection. The multi-resolution analysis (MRA) provides a way of dividing the data so that behaviours at specific frequencies (or scales) can be seen. The original data can be reconstructed from the wavelet and scaling coefficients. If, however, we set all coefficients to zero, with the exception of those for a given scale, and then do the reconstruction we obtain the component of the original signal corresponding to that scale. The original signal is given by the sum of all the components.

Wavelet packet analysis

We have already pointed out some of the advantages of wavelets, but using wavelets to describe data involves compromise. It is a compromise between (a) the original data in which each sampling point has its unique value but no information on frequency and (b) the spectrum which gives the best resolution in the frequency domain but no information on the locations of features. The DWT and MODWT representations lie between these extremes. As frequency increases, the DWT and MODWT have decreasing resolution in frequency, and it follows increasing resolution in space. This is illustrated in Fig. 8a where we show how these wavelet transforms divide the frequency interval (we note that scale is simply the inverse of frequency).

A further development of the wavelet techniques is the maximum overlap discrete wavelet packet transform (MODWPT), described admirably by Percival and Walden (2000). This transform gives more flexibility in the way space and frequency are resolved. At each level of decomposition k the frequency interval between 0 and 1/2l can be divided into [2.sup.k] equal intervals, known as wavelet packets. Here k = 1, 2 , ..., m, where m is an integer such that [2.sup.m] < N. We denote the ith packet at level k as [P.sub.k,i]. For example, [P.sub.4,1] is the first packet at level 4 (see Fig. 8b). A set of non-overlapping wavelet packets that cover the whole frequency interval forms a complete basis. Figure 8b illustrates such a complete basis by the shaded packets {[P.sub.3,1], [P.sub.3,2], [P.sub.2,2], [P.sub.4,9] [P.sub.4,10], [P.sub.3,6] [P.sub.2,4]

Each wavelet packet has an associated filter, which when convolved with the data produces N wavelet packet coefficients, rather like the MODWT filter for a given value of j.

One can select the MODWPT basis to suit the data by trading resolution in space for resolution in the frequency domain. Usually, it is desirable to select a basis that concentrates variance in as few wavelet coefficients as possible (Percival and Walden 2000; Lark 2006). This results in good frequency resolution over stationary stretches and good spatial resolution at frequencies with short-range episodic noise. In practice, we identify the 'best' basis, i.e. the basis that minimises a suitable cost function. Here we use a cost function that is, in effect, the DWPT Shannon entropy averaged over all possible positions, x, in space (Constantine and Reinhall 2001).

Like the DWT and MODWT described above, the MODWPT coefficients can be used in a multi-resolution analysis to estimate the partition of variance between frequency intervals (the wavelet packet variance) and to test for significant changes in variance across the space for a given packet. Convention is that results are presented in terms of frequency, and where appropriate 'wavelength', the reciprocal of frequency. This is because, unlike wavelet functions, wavelet packet functions are not simply dilations and translations of one another.

Testing a hypothesis that data are stationary in the variance

The MODWPT basis offering the best frequency resolution is given by [P.sub.m,i], where i= 1, 2, ..., [2.sup.m]. If we regard the data as a realisation of a stationary random process then we expect the best basis to comprise mostly packets from level m, i.e. the one with best frequency resolution. Lark (2007) proposed that we use the ratio of the entropy of the best basis coefficients to the entropy of the coefficients from the basis [P.sub.m,i], where i- 1, 2, ..., [2.sup.m], as a test statistic for the null hypothesis that the data are a realisation of a stationary process. A sampling distribution of this ratio under the null hypothesis is generated by Monte Carlo simulation of datasets specifying the variogram fitted to the empirical variogram of real data under study.

[FIGURE 11 OMITTED]

Two-dimensional wavelet transforms

Lark and Webster (2004) showed how the MODWT analysis can be extended to 2 dimensions. They then applied the technique to study the scale-specific behaviour of soil thickness, potential solar radiation, and slope gradient of a landscape in south-west England. They revealed the structure at several scales by inspecting the 2-dimensional multi-resolution analysis, and they showed how to estimate the partition of variation (and correlation) across the scales by calculating wavelet variances.

[FIGURE 12 OMITTED]

The 2-dimensional wavelet transform of [N.sub.x] x [N.sub.y] gridded data results in 4 sets of coefficients at each level of decompositionj. The first set comprises the scaling function coefficients, [s.sub.j](x,y), where x = 1, 2, ..., [N.sub.x], and y = 1, 2, ..., [N.sub.y] where [N.sub.x] and [N.sub.y] are the numbers grid nodes in the x and y dimensions, respectively. Each coefficient is a weighted average of data near the position of the coefficient. The second set, [d.sub.Hj](x, y), describes the variation at scale j along the rows of the grid, and the third set, [d.sub.vj](x, y), describes that down the columns. The fourth set, [d.sub.Gj](x, y), comprises the diagonal wavelet coefficients.

Conceptually the word 'diagonal' is somewhat misleading because it implies that these coefficients capture the diagonal variation in the data. However, if we consider the first level in the 2-dimensional Haar wavelet transform (as discussed in Lark and Webster 2004) we see this is not the case.

The formulae for the wavelet variances (i) along the rows, (ii) down the columns, and (iii) for the diagonal components are given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], (11)

for scale interval number j, where m is H for row, V for column, and G for diagonal components. To estimate isotropic wavelet variances we add together the horizontal, vertical, and diagonal estimates. Where variation is isotropic we expect the horizontal and vertical components to be similar. In the appendix we use the Haar wavelet to show that if the data can be thought of as a realisation of an isotropic, stationary process with decreasing covariance function, then we expect the variance of the diagonal component to be less than the horizontal and vertical components.

[FIGURE 13 OMITTED]

Results

We used the MODWT in 1 and 2 dimensions and the MODWPT to investigate the data from the 2 aerial photographs. In particular, we used the 2-dimensional multi-resolution analysis to partition the data into scale-dependent components, which gave us a qualitative indication of the characteristics of the data across the scales. To get a more quantitative measure we also calculated the 2-dimensional wavelet variance and identified the dominant scales of variation. We studied 5 sample transects in more detail using the 1-dimensional wavelet analysis. We chose one of the transects (Back Creek, row 100) to cross the north-east corner of the rectangle where the gilgais are much darker, and we presume because they were wetter, than elsewhere. In particular, we investigated whether the underlying processes can be thought of as stationary in the variance. First we tested this hypothesis using the wavelet packet method proposed by Lark (2007).

Second we looked for significant changes in the variation across the space. We also calculated wavelet variances and wavelet packet variances to identify the dominant scales (or frequencies) of variation. In all of the analysis we used Daubechies's extremal phase wavelet with 2 vanishing moments (Daubechies 1992). We selected this wavelet because it has a very compact support (i.e. it takes non-zero values over a narrow interval and so is particularly suitable for identifying local features in the data).

[FIGURE 14 OMITTED]

Two-dimensional multi-resolution analysis

Figures 9 and 10 display the results of the multi-resolution analysis of Caragabal and Back Creek, respectively. Considering first the smooth components we see that Back Creek (Fig. 10a) has a clear coarse pattern (scales >64). It arises because of the 2 fairly large regions, e.g. the north-east comer of the scene, where the gilgais are darker than elsewhere. The gilgais were probably wetter there at the time and the vegetation lusher. The isotropic detail components (Figs 9b and 10b) show the regular uniform pattern of the gilgais at the scale interval 16-32 pixels. For both images, variation is most pronounced at this scale interval.

Two-dimensional wavelet variance

Figure 11 shows the horizontal, vertical, and diagonal wavelet variances for Caragabal and Back Creek. The results for the 2 sites are broadly similar. There is a pronounced peak in the wavelet variance at the 16-32 pixel scale. This supports the observation made on the detail components of the multi-resolution analysis, namely that the 16-32 pixel scale has the most pronounced variation. The horizontal and vertical components of variation are similar in magnitude across the scales, and the diagonal variation is consistently less than the other two at scales <32 pixels. As we note above, this accords with an isotropic spatial pattern.

One-dimensional wavelet variance

Figure 12 shows the partition of wavelet variance for the 5 transects. There is a similar pattern of decomposition across them all. The largest contribution to variation in all but column 300 at Caragabal is at the 16-32 pixel scale.

The wavelet packet (WP) variance was also estimated for each transect decomposed on its own best basis. As the frequency interval changes from packet to packet, a simple graph of WP variance against frequency can mislead. Therefore, we have plotted the WP variance divided by the associated frequency interval against the midpoint of the frequency interval. This quantity is a standardised WP variance for the frequency interval. Figure 13 shows that the largest contributions to variation occur at frequency intervals 0.023-0.027 for row 400 and column 400 at Back Creek and column 300 at Caragabal. For row 400 at Caragabal the largest contribution occurs in the frequency interval 0.031-0.039, and for row 100 at Back Creek in the interval 0.008-0.016. This last result supports the above observation on the MODWT wavelet variance, where we saw that the relative contribution of the 64-128 pixel scale variation is inflated compared to the other transects. This large peak is caused by the low frequency variation between the darker gilgais and the remaining area. If we consider the series from location 95 onward only, we get results more consistent with the other transects (see Figs 12 and 13). In particular, the peak in wavelet packet variance is now at the frequency interval (0.027-0.031).

One-dimensional best basis and detection of change in variance

The best bases for each of the sample transects are shown in Fig. 14. For all transects except row 100 at Back Creek, at frequencies <0.047 (periods >21.33 pixels) the best bases achieve good resolution in the frequency domain (i.e. the packets are from levels 6 and 7) at the expense of resolution in space. For row 100 at Back Creek, the best basis has poorer resolution at these low frequencies. This is not surprising, because, as we observed in the 2-dimensional multiresolution analysis, the changes associated with the darker gilgais affect the coarser scales (i.e. low frequencies). Good resolution of frequency would reduce the spatial resolution of this feature, and so the compression of its contribution to the variation of the series into the smallest possible number of wavelet packet coefficients. Therefore, the best basis appropriately selects packets with greater spatial resolution.

For each transect, the probability for the null hypothesis that the underlying process is stationary is small (<0.03), so we reject the hypothesis. We found significant changes in wavelet packet variation on all transects. The packets where changes were detected are identified on Fig. 14 by the darker shading. For all but row 100 at Back Creek, we detected no significant changes in variation in packets corresponding to frequencies <0.094 (periods >10.66 pixels). Row 100 at Back Creek does have significant changes in variation at low frequencies. These are at positions 101 ([P.sub.5,2]) and 94 ([P.sub.4,2]), i.e. bounding the region with the darkest gilgais.

Discussion and conclusion

Our Introduction looked forward to the insight into gilgai distribution that we might gain with the aid of modem digitising equipment, hugely increased computing power over the last 35 years, and new developments in wavelet analysis. We have not been disappointed. The first impression from the visual images is that the gilgais repeat in apparently regular patterns. This is confirmed by the peaks in the 2-dimensional spectra and reinforced in the wavelet variances. This is despite the evidence of non-stationarity of both short and, at Back Creek, of long range (trend).

The 2-dimensional MODWT analysis shows the largest contribution to variation in both horizontal and vertical dimensions is at the 16-32 pixel scale. At both Caragabal and Back Creek, horizontal and vertical variation is similar; the variation is evidently isotropic. This enables us to draw conclusions about the 2 dimensions from 1-dimensional analyses, for which we have at our disposal a broader range of techniques.

The 5 sample transects we investigated with the discrete wavelet transform showed large contributions to variation at the 16-32 pixel scale, as in the 2-dimensional analysis. The wavelet packet analysis provides more detail. It suggests that we should not treat the data as realisations of a stationary process, particularly at frequencies >0.094. These frequencies are greater than the frequency we expect to represent the features of central interest, the gilgais, which are ~34 m apart (corresponding to frequency 0.038 per pixel). In fact, the interval of periods with the largest wavelet packet variation approximates the distance between the centres of neighbouring gilgais. These large wavelet packet variances correspond to the peaks in the spectrum.

We conclude that the gilgais do indeed occur in regular uniform patterns across the landscapes studied, even though the data on which our analysis is based exhibit more complex and irregular variation at both higher and lower frequencies than the gilgai pattern.

Appendix 1

Let us consider the first level decomposition of an [N.sub.x] x [N.sub.y] grid of data using the Haar wavelet. We identify the datum from the xth column and yth row by z(x, y), but for brevity we denote it by [p.sub.1] =z(x, y). In like manner [p.sub.2] =z(x+ 1, y), [p.sub.3] =Z(x, y+ 1) and [p.sub.4]=Z(x + 1, y+ 1). The 1-dimensional Haar wavelet is effectively the difference between 2 neighbouring points. Lark and Webster (2004) showed that the 2-dimensional wavelet coefficients [d.sub.M,1](x, y) are given by

[d.sub.M,1] (x,y) = [4.summation over (i=1)][[lambda].sub.M,1],1[p.sub.i], (12)

where M is H for row, V for column or G for the diagonal coefficient and:

-[[lambda].sub.H,1] =-[[lambda].sub.H,2] = [[lambda].sub.H,3] = [[lambda].sub.H,4] = [lambda]

[[lambda].sub.V,1] = [[lambda].sub.V,3] = [[lambda].sub.V,2] = [[lambda].sub.V,4] = [lambda]

and

-[[lambda].sub.G,1] = -[[lambda].sub.G,4] = [[lambda].sub.G,2] = [[lambda].sub.G,3] = [lambda]

The wavelet variance is given by

[[sigma].sup.2.sub.M,1] = [4.summation over (i=1)] [4.summation over (j=1)] [[lambda].sub.M,i] [[lambda].sub.M,j] C(i, j), (13)

where C(i, j) is the covariance between points [p.sub.i] and [p.sub.j]. If the underlying process is isotropic and stationary then we define the covariances at lags 0, 1 and [square root of 2] as [C.sub.0], [C.sub.1], and [C.sub.[square root of 2]], respectively. This means that:

[[sigma].sup.2.sub.M,1] = [[lambda].sup.2]([C.sub.0] - [C.sub.[square root of 2]])

for M = H or V and:

[[sigma].sup.2.sub.G,1] = [[lambda].sup.2]([C.sub.0] - 2[C.sub.1] + [C.sub.[square root of 2]]).

If the data are random noise then [C.sub.0] = [C.sub.[square root of 2]] = 0 and the wavelet variances are equal. If, however, the data have a decreasing covariance function such that

[C.sub.0] > [C.sub.1] > [C.sub.[square root of 2]] (14)

then [[sigma].sup.2.sub.H,1] = [[sigma].sup.2.sub.G,1]. The wavelet coefficients [d.sub.M,j] (x, y) at j > 1 are calculated in a similar way to Eqn (12), except that in these cases [p.sub.i] is the scaled average of a [2.sup.j] x [2.sup.j] grid of data centred near x, y. Therefore similar arguments hold.

doi: 10.1071/SR09189

Acknowledgments

This work was undertaken within Rothamsted Research's programme in Computational and Mathematical Biology, which is funded by the Biotechnology and Biological Sciences Research Council (BBSRC).

Manuscript received 26 October 2009, accepted 3 February 2010

References

Biswas A, Si BC, Walley F (2008) Spatial relationship between 815N and elevation in agricultural landscapes. Nonlinear Processes in Geophysics 15, 397-407. doi:10.5194/npg-15-397-2008

Chatfield C (1984) 'The analysis of time series: an introduction.' (Chapman and Hall: London)

Constantine WLB, Reinhall PG (2001) Wavelet-based in-band denoising technique for chaotic sequences. Journal of Bifurcation and Chaos 11, 483-495. doi:10.1142/S0218127401002201

Daubechies 1(1992) 'Ten lectures on wavelets.' (Society for Industrial and Applied Mathematics (SIAM): Philadelphia)

Hallsworth EG, Robertson GK, Gibbons FR (1955) Studies in pedogenesis. VII. The 'gilgai' soils. Journal of Soil Science 6, 1-31. doi:10.1111/ j.1365-2389.1955.tb00826.x

Lark RM (2005) Spatial analysis of categorical soil variables with the wavelet transform. European Journal of Soil Science 56, 779-792.

Lark RM (2006) The representation of complex soil variation on wavelet packet bases. European Journal of Soil Science 57, 868-882. doi: 10.1111/j.1365-2389.2005.00779.x

Lark RM (2007) Inference about soil variability from the structure of the best wavelet packet basis. European Journal of Soil Science 58, 822-831. doi: 10.1111/j.1365-2389.2006.00872.x

Lark RM, Milne AE, Addiscott TM, Goulding KWT, Webster CP, O'Flaherty S (2004a) Analysing spatially intermittent variation of nitrous oxide emissions from soil with wavelets and the implications for sampling. European Journal of Soil Science 55, 611-627. doi: 10.1111/j.1365-2389.2004.00620.x

Lark RM, Milne AE, Addiscott TM, Goulding KWT, Webster CP, O'Flaherty S (2004b) Scale- and location-dependent correlation of nitrous oxide emissions with soil properties: an analysis using wavelets. European Journal of Soil Science 55, 601-610. doi:10.1111/ j.1365-2389.2004.00629.x

Lark RM, Webster R (1999) Analysis and elucidation of soil variation using wavelets. European Journal of Soil Science 50, 185-206. doi: 10.1046/ j.1365-2389.1999.t01-1-00234.x

Lark RM, Webster R (2001) Changes in variance and correlation of soil properties with scale and location: analysis using an adapted maximal overlap discrete wavelet transform. European Journal of Soil Science 52, 547-562. doi:10.1046/j.1365-2389.2001.00420.x

Lark RM, Webster R (2004) Analysing soil variation in two dimensions with the discrete wavelet transform. European Journal of Soil Science 55, 777-797. doi:10.1111/j.1365-2389.2004.00630.x

McKenzie N, Jacquier D, Isbell R, Brown K (2004) 'Australian soils and landscapes: an illustrated compendium.' (CSIRO Publishing: Collingwood, Vic.)

Milne AE, Lark RM, Addiscott TM, Goulding KWT, Webster CP, O'Flaherty S (2005) Wavelet analysis of the scale- and location-dependent correlation of modelled and measured nitrous oxide emissions from soil. European Journal of Soil Science 56, 3-17. doi:10.1111/j.1365-2389.2004.00650.x

Percival DB, Guttorp P (1994) Long-term memory processes, the Allen variance and wavelets. In 'Wavelets in geophysics'. (Eds E Foufoula-Georgiou, P Kumar) pp. 325-344. (Academic Press: New York)

Percival DB, Walden AT (2000) 'Wavelet methods for time series analysis.' (Cambridge University Press: Cambridge, UK)

Preston FW, Davis JC (1976) Sedimentary porous materials as a realization of a stochastic process. In 'Random processes in geology'. (Ed. DF Merriam) pp. 63-86. (Spfinger-Vedag: New York)

Priestley MB (1981) 'Spectral analysis and time series.' (Academic Press: London)

Russell JS, Moore AW (1972) Some parameters of gilgai microrelief. Soil Science 114, 82-87. doi:10.1097/00010694-197208000-00002

Webster R (1977) Spectral analysis of gilgai soil. Australian Journal of Soil Research 15, 191-204. doi:10.1071/SR9770191

Webster R, Oliver MA (2001) 'Geostatistics for environmental scientists.' (John Wiley and Sons: Chichester, UK)

Whitcher BJ, Guttorp P, Percival DB (2000) Multiscale detection and location of multiple variance changes in the presence of long memory. Journal of Statistical Computation and Simulation 68, 65-87. doi:10.1080/00949650008812056

A. E. Milne (A), R. Webster (A,B), and R. M. Lark (A)

(A) Rothamsted Research, Harpenden, Hertfordshire AL5 2JQ, Great Britain.

(B) Corresponding author. Email: richard.webster@bbsrc.ac.uk

Table 1. Summary statistics of grey level of gilgai images Caragabal Back Creek Rows 563 498 Columns 570 458 Mean 129.4 147.1 Median 127.0 150.0 Variance 808.0 1212.0 s.d. 28.4 34.8 Skewness 0.10 -1.53 Table 2. Conversion between frequency (pixel) and scale (in units of pixels and metres) for values mentioned in the text Note that the units of scale and wavelength are equivalent, but it is convention to discuss the DWT and MODWT in terms of scale, and spectral analysis and wavelet packets in terms of frequency and wavelength Frequency Scale Scale (pixel 1) (pixel) (m) 0.500 2.0 2.6 0.250 4.0 5.2 0.125 8.0 10.4 0.094 10.6 13.9 0.063 16.0 20.9 0.047 21.3 27.8 0.040 25.0 32.6 0.031 32.0 41.8 0.025 40.0 52.2 0.016 64.0 83.6 0.008 128.0 167.1 0.004 256.0 334.3

Printer friendly Cite/link Email Feedback | |

Author: | Milne, A.E.; Webster, R.; Lark, R.M. |
---|---|

Publication: | Australian Journal of Soil Research |

Article Type: | Report |

Geographic Code: | 8AUST |

Date: | Jul 1, 2010 |

Words: | 6544 |

Previous Article: | Soils fieldwork, analysis, and interpretation to support hydraulic and hydrodynamic modelling in the Murray floodplains. |

Next Article: | Temporal and spatial patterns of salinity in a catchment of the central wheatbelt of Western Australia. |

Topics: |