| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
a Soil-Plant Dynamics Unit, Scottish Crop Research Institute, Invergowrie, Dundee, DD2 5DA, UK
b SIMBIOS, University of Abertay Dundee, Kydd Building, Bell Street, Dundee DD1 1HG, Scotland
c National Soil Resources Institute, Cranfield University, Silsoe, Bedfordshire MK45 4DT, UK
* Corresponding author (j.crawford{at}simbios.ac.uk).
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONMATERIALS AND METHODSCONCLUSIONS AND DISCUSSIONREFERENCES |
|---|
Abbreviations: MCMC, Markov chain Monte Carlo • MRF, Markov random fields
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSCONCLUSIONS AND DISCUSSIONREFERENCES |
|---|
The most useful of these models has been used to interpret the impact of structure on physical properties and processes; but comparatively little work has examined the impact on biology. Some attempt has been made to link biological processes with soil structure; these have generally been limited to N transformations (Young and Ritz, 2000). These studies have clearly demonstrated the importance of understanding the relative spatial distribution of the physical and biotic elements of soil structure in determining the larger-scale properties of the resultant biological process (Arah et al., 1997; Rappoldt and Crawford, 1999). Therefore there is a need to further develop models of soil structure that are capable of integrating the physical and biological heterogeneities that occur in most soils. However there are a number of challenges that place constraints on appropriate methodologies, and one of the most significant of these is the limitation of existing imaging technology.
The only reliable methods for visualizing soil in three dimensions are 
and X-ray tomography (Rogasik et al., 1999). While the technology is rapidly improving, it is not a trivial matter to differentiate between pore and solid matrix in these visualizations. Although resolutions of 5 µm or higher are now possible, it is still not possible to directly image soil microbes in situ and in three dimensions over comparable scales. The only method for simultaneously imaging soil structure and the distribution of microbes is by using biological thin sections (e.g., Nunan et al., 2001, 2003). Therefore any model of structure must be capable of being parameterized from two-dimensional data and extrapolated to three dimensions.
The requirements of a useful model capable of describing the heterogeneity of physical and biological elements in soil are three-fold. First, the models must be able to describe the spatial structure of multiphase media (matrix, pore, microbe etc.) at the scale of individual pores and microbes. Second, the model must also accommodate any spatial anisotropy inherent in soil. Finally, the model should be capable of using parameters determined from two-dimensional sections. Currently, no method exists that has been demonstrated to simultaneously satisfy these requirements. Yeong and Torquato (1998) state that their method can be developed to satisfy these constraints, although to date no such modification has appeared. The development of their method based on correlation functions, to multiphase anisotropic media is far from a trivial extension of the existing methodology.
In this paper, we introduce an efficient method for modeling the architecture of soil, based on two-dimensional images obtained from soil thin sections. The algorithm is based on a Markov chain process that permits the macroscopic (centimeter scale) properties to be determined from local (pore scale) conditions. A Monte Carlo implementation is used to reproduce the stochastic properties of the architecture and to estimate features of the posterior or predictive distribution of interest by using samples drawn from images derived from real soil.
Markov chain Monte Carlo methods usually employ an iterative scheme to obtain the final spatial description (Besag and Green, 1993), and where correlated structures of the order of the image size exist (such as macropores in soil thin sections) the iterative schemes are slow to converge or may even fail. The development of MCMC presented in this paper avoids these problems by a novel choice of neighborhood for the Markov chain and the use of a scanning scheme, based on Elfeki and Dekking (2001) but without the requirement for preconditioning, in which the model image is constructed from a single-pass raster. The resulting structure models are quantitatively compared with the original images to test the validity of the approach. A brief discussion of the applicability of the method to extrapolate to three-dimensions is presented, although a full treatment is deferred to a forthcoming paper. Finally, the trivial extension to simultaneously modeling the biological and physical components of soil architecture in three-dimensions is outlined.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSCONCLUSIONS AND DISCUSSIONREFERENCES |
|---|
In using the MRF method to model visualizations such as those in the current application, a central assumption is that the state of the structure at some point conditionally depends on only a relatively small number of points in a predefined neighborhood. Implementation of the algorithm commences with the derivation of the conditional probabilities from direct measurements of the probabilities of different neighborhood structures in an image of soil structure. The generation of model structures starts with an initial estimate of the spatial distribution of states (e.g., a spatially random distribution). The structure at each point is then updated in accordance with the conditional probabilities derived from the image. This update is then repeated by successive applications of the conditional probabilities at each point, until the statistical properties of the resulting spatial distribution converge (i.e., do not change significantly between successive iterations). Larger-scale correlations emerge as a consequence of the local dependency built into the conditional probabilities, but only after a number of iterations of the algorithm. Indeed a major limitation of the application of the conventional MRF method is the tendency to a convergence slowly and an associated high computational demand. This limits the application to systems where the number of states of each point is low, and where there are only relatively short-range spatial correlations in the structure. Furthermore, the standard implementation of the method cannot reproduce nonisotropic structures, nor can it recognize concavity or convexity of shapes (Tjelmeland and Besag, 1998). Thus, this renders the method inappropriate for modeling most soil structures.
Elfeki and Dekking (2001) employed a raster scheme within a MCMC approach to simulate geological strata. Their model was parameterized from data collected from a set of wells, and they assumed that the state of a particular point in the strata was dependent on the points immediately to one side and to the top. Therefore if the state of the points are known along a transect at the surface and in a vertical direction at a single well location, it is possible to use the model to predict the state of the remaining points in the transect-well plane. The model calculates the state of a particular point on the basis of the state of the point to the left and above, and the conditional probabilities associated with the Markov process. This is repeated in a step-wise way as the chain moves from left to right across the domain in a raster fashion. Crucially, however, to get an accurate representation the chain must be conditioned on several wells across the transect. Therefore the probabilities are adjusted to guarantee agreement with the well data (Elfeki and Dekking, 2001). This approach is appropriate for data on geological strata since the number of wells that can be drilled always limits one, and the aim is to reproduce the actual structure as closely as possible. The purposes are quite different to those of the algorithm presented here. In our case, we have complete information about the state of the points in a two-dimensional domain (i.e., soil thin section), and we aim to reproduce the functionally important statistical properties of the structure rather than a literal copy of the structure itself. The method should also be extendable to three dimensions, and for the reasons outlined above, we are restricted to methods that can extrapolate from two-dimensional data. The method of Elfeki and Dekking (2001) is then inappropriate.
Multidimensional Markov Chain Model
We consider an image made up of pixels arranged in a rectangular array. The standard implementation of the MRF method assumes that the state of a particular point in an image depends on the state of an isotropic neighborhood centered about the point. As this is unsuitable for modeling anisotropic soil structures, we proceed by removing this constraint.
The Potential Function
The standard MRF method calculates the required conditional probability from a potential function defined in an isotropic neighborhood. The fundamental framework contains the following components:
S} of (pixel) sites; N} attached to the site;
 | [1] |
|
 | [2] |
x, if all states x in the neighborhood nx are equal, and –x otherwise. For higher-order interactions, the potential functions are assumed to be a linear function of parameters, which can be derived from an image only under such simplified assumptions. The advantage of determining the probabilities from the potential function is that the neighborhood interactions of an entire image can be represented by relatively few parameters. However, the apparent compactness is a consequence of the assumption that the potential functions can be expressed as a linear function of the parameters. Because of this, the parameters can be estimated using maximum-likelihood methods (Qian and Titterington, 1991). This appears to work for many kinds of images, but its applicability to images of soil, where relatively long-range correlations exist, has not been verified. We attempted to model the structure of our soil samples by calculating the potential function. However, the resulting model structures were inadequate. For moderate-sized images, the modeled structure showed substantial departures from the original (parent) structure, and as the size of the modeled domain increased, the agreement deteriorated further (Fig. 2) . Thus, the simplifying assumptions underlying the linear formulation of the potential functions are incompatible with the image structures associated with soil. We therefore developed an alternative method based on direct measurement of the 64 conditional probabilities associated with all possible configurations of the five-pixel neighborhood (Fig. 1).
|
The standard methods for implementing MCMC are iterative, for example, the Gibbs sampler version of the Metropolis Hastings algorithm (Geman and Geman, 1984). However, these suffer from long convergence times as discussed above, and are not suitable for modeling large correlated structures such as are found in soil. Here, we use a more efficient method based on the scanning scheme algorithm proposed by Qian and Titterington (1991), modified to cope with long-range correlations in the structural heterogeneity found in soil images. The modification replaces the potential function with a more explicit determination of the transition probabilities as detailed below.
The Scanning Scheme Algorithm
After considering several different forms for the neighborhood, we determined that the neighborhood given in Fig. 1 was the smallest capable of reproducing the observed soil properties. The modeling proceeds in two steps. First, the state of the pixel located at the point (i,j) is determined from knowledge of those at (i,j – 1), (i – 1, j – 1), (i – 1,j), and (i – 1,j + 1) using the associated four-neighborhood conditional probability. Second, the state of the pixel at (i,j + 1) is obtained from knowledge of the new state at (i,j) together with the state of the those at (i,j – 1), (i – 1,j – 1), (i – 1, j), and (i – 1, j + 1) using the associated five-neighborhood conditional probability. These probabilities are obtained from the original image by sampling the four- and five-neighborhoods, and enumerating the different realizations of the state of the point (i,j) and (i,j + 1) respectively for each configuration of the neighborhood.
To initiate the model we need to assign states to all the cells in the first row, and the first cell of the next row. The cells in the first row are obtained using a two-neighborhood Markov chain where the state of a cell is conditionally dependent on the state of the cell to its left. The parameters for this Markov chain are obtained as above, but using the smaller two-cell neighborhood. The state of the first cell of the next row is determined using the same two-neighborhood conditional probability. Using these boundary values, the chain runs from the left-hand corner of the image and progresses in raster fashion across the image to the right-hand side. First, the state of cell (i,j) in the neighborhood is evaluated, followed by the state of cell (i,j + 1). The neighborhood is then advanced two cells to the left and the process is repeated. This continues until the last cell on the right-hand side is reached and its state is evaluated. At this point, the cells in the first two rows have been evaluated. Next, the neighborhood remains at the right-hand side but moves one row down. The chain now reverses direction, and instead of deriving the states of the (i,j) and (i,j + 1) cells in terms of the state of the others, it is the states of cells (i,j) and (i, j – 1) that are determined. However, before the chain can proceed leftwards, the state of the first cell on the right-hand side of the third row must be evaluated. This is done in the same way as when the neighborhood was at the left-hand side of the domain, using the two-neighborhood conditional probability. The chain can now advance leftwards until the left-hand side of the domain is reached. The neighborhood then moves down one row, and the chain reverses as before. Thus, the whole domain is scanned in a raster-like fashion on this basis, until the states of all the required cells are obtained.
The scanning scheme algorithm converges rapidly. In the examples reported here, we observed the transition kernel (i.e., the matrix of conditional probabilities for all possible neighborhood configurations), calculated from the reconstructed image as the chain progressed. Almost all the probabilities had become stable after the chain had completed 200 rows, which is equivalent to a depth of 4.0 mm in the original soil section and takes about 10 s of computing time on a 1.7-GHz Pentium IV computer. In other words, the minimum size of a simulated representative soil image should be 0.6 by 0.4 cm2, and it takes only a few minutes to generate an image covering several squared.
Validation
The method was validated using images obtained from soil thin sections and selected to represent a broad contrast in structural properties.
Sampling
Soil cores were collected from an arable field and thin sections were produced using the method described in Nunan et al. (2001). Soil pore maps were obtained by subtracting images obtained with cross-polarized light from images captured using transmitted bright-field light. The resultant images were then segmented into solid and void. The images employed in this study were binary pore maps of dimension 760 by 570 pixels, representing an area of 1.6 by 1.2 cm. Images were selected to represent a range of characteristic soil structural properties, as shown below.
Comparison of Real and Simulated Images
To compare the simulated and real images we selected a range of quantitative metrics that characterize the heterogeneity and connectivity of the structures under investigation. The most obvious of these is the porosity and this is readily determined from both the real and simulated structures. The corresponding values are listed in Table 1 and range from 7 to 24% in the real structures. There was no significant difference between the porosities in the simulated and real structures (P > 0.05, paired t test).
|
The third and fourth metrics chosen characterize the pore space, where visually obvious differences between the samples were present. The third is the variance of the porosity as measured in a 0.4 by 1.6 mm sampling window placed at 50 random locations in each image. For a given porosity, this is a measure of the connectivity of the pore space (Mandelbrot, 1985). We used the Chi-Square goodness-of-fit test, and tested the null hypothesis that the variances in each pair of simulated and original images were different. This could be rejected at the 95% confidence level indicating that this aspect of the structure of the pore space was not significantly different in the real and simulated images. The fourth metric measured the spatial correlation of the pore space by determining the semivariogram. Figure 3 shows the variograms for the different soil samples used in this study, and there are clear differences in these between the different soil sections. The figure shows the comparison between the variograms for the real and simulated structures. There is no formal statistical way of comparing the properties of semivariograms, however the high degree of correspondence between the curves for the measured and simulated structures is good, adding further support for the modeling methodology.
|
|
CONCLUSIONS AND DISCUSSION |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSCONCLUSIONS AND DISCUSSIONREFERENCES |
|---|
The method is easily extendable to simultaneously model the spatial distribution of a variety of components of soil architecture. In the current paper, we verified the approach by modeling the physical elements hence each pixel could be in one of two states—pore or solid. In a forthcoming publication, we have extended the approach to model the relative spatial distribution of microbes in soil, parameterized directly from soil thin sections that have been prepared in a manner that preserves the microbes in situ (Nunan et al., 2003).
While it is of interest to be able to produce two-dimensional models of the different components of soil architecture that have the same statistical properties as real soil, the most important advantage of the current approach is that is has the ability to extrapolate to three dimensions. The efficiency of the algorithm together with a Markov process that is based on local neighborhood, without the need to condition the chain on existing data, makes it possible to consider modeling large three-dimensional structures. This means it is possible to produce models for the three-dimensional distribution of microbes in structured soil for the first time. This can be combined with efficient algorithms for modeling unsaturated flow through porous media (e.g., Zhang et al., 2002) and the distribution of O2 (Rappoldt and Crawford, 1999) to begin to understand how the physical and biological processes in soil interact. Again, the ability to relate pore- to core-scale means we can approach the issue of possible scaling laws from first principals. By linking microbes and their microhabitats directly in this way, the potential for theoretical and experimental soil ecology is significantly advanced.
|
ACKNOWLEDGMENTS |
|---|
Received for publication June 4, 2002.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONMATERIALS AND METHODSCONCLUSIONS AND DISCUSSIONREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  J. M. Blair, R. E. Falconer, A. C. Milne, I.M. Young, and J. W. Crawford ModelingThree-Dimensional Microstructure in Heterogeneous Media Soil Sci. Soc. Am. J., October 29, 2007; 71(6): 1807 - 1812. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
W. Li Transiograms for Characterizing Spatial Variability of Soil Classes Soil Sci. Soc. Am. J., May 16, 2007; 71(3): 881 - 893. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
E. Park, A. M. M. Elfeki, Y. Song, and K. Kim Generalized Coupled Markov Chain Model for Characterizing Categorical Variables in Soil Mapping Soil Sci. Soc. Am. J., May 16, 2007; 71(3): 909 - 917. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
W. Li and C. Zhang A Random-Path Markov Chain Algorithm for Simulating Categorical Soil Variables from Random Point Samples Soil Sci. Soc. Am. J., April 5, 2007; 71(3): 656 - 668. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
W. Li, C. Zhang, J. E. Burt, and A-X. Zhu A Markov Chain-Based Probability Vector Approach for Modeling Spatial Uncertainties of Soil Classes Soil Sci. Soc. Am. J., October 27, 2005; 69(6): 1931 - 1941. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| The SCI Journals | Agronomy Journal | Crop Science | |||
| Journal of Natural Resources and Life Sciences Education |
Vadose Zone Journal | ||||
| Journal of Plant Registrations | Journal of Environmental Quality |
The Plant Genome | |||
