Cokriging with ground-based radiometry

REMOTE SENS. ENVIRON. 41:45-60 (1992)

Cokriging with Ground-Based Radiometry P. M. Atkinson Department of Geography, University of Bristol, Bristol, United Kingdom

R. Webster INRA Science du Sol, MontpeUier, France.

P. J. Curran Department of Geography, University College of Swansea, Swansea, United Kingdom

T h e soil and crop cover of agricultural land vary spatially in a way that is both random and autocorrelated. This enables them to be estimated and mapped from sparse sample data by kriging. These properties (primary variables) are usually related to the radiation they reflect: They are coregionalized with it. In many circumstances the primary variables can be estimated more precisely by measuring, in addition, the radiation sparsely, using a ground-based radiometer and combining the two by cokriging. The coregionalization must be formalized in a coherent set of variograms, one for the primary variable, one for each variable derived from the radiometry, and the cross variograms between all pairs of variables involved in the estimation. Given this set, it is possible to determine estimation variances for any configuration of sampling and to design an optimal scheme that will achieve a desired precision for least effort. The formulae for cokriging are presented, as are the conditions for a coherent model of the coregionalization, and the article shows how these can be used to design sampling schemes that combine survey of the primary variables and radiation to best Address correspondence to Dr. Peter Atkinson, Department of Geography, University of Bristol, University Road, Bristol, BS81SS, United Kingdom. Received 13 May 1991; revised 8 February 1992. 0034-4257 / 92 / $5. O0 ©Elsevier Science Publishing Co. Inc., 1992 655 Avenue of the Americas, New York, NY 10010

advantage. Three examples from intensive agriculture in Britain illustrate the technique. In one example where the aim was to estimate and map the cover of clover in pasture, cokriging using measured radiation was nine times as efficient as kriging the cover alone. INTRODUCTION In almost all applications of remote sensing for estimating properties of the Earth's surface, such as the amount of vegetation, the moisture content of the soil, and the suspended sediment load in a body of water, the radiation is not itself of prime interest. The radiation is usually closely related, however, to properties of the Earth's surface that we do wish to know about, and we measure radiation so that we can estimate the latter from it. The relations are usually expressed in the form of a regression equation, which is then used to predict the unknown property at the ground from the measured radiation. Data from sensors in aircraft and satellites can provide complete cover, and so the procedure gives estimates for every position on the ground. In detailed studies of vegetation, the radiation can be measured in a much more controlled and

45

46

Atkinson et al.

accurate way from just above the canopy. Groundbased radiometers are potentially more accurate both geometrically and radiometrically than airborne or spaceborne sensors. In addition: i. They enable the measurement of radiation at time intervals to be set by the investigator rather than the aircraft or satellite. ii. The field-of-view is small and can be chosen to include only the vegetation canopy of interest at the time. iii. They record a signal that is largely free of atmospheric interference or attentuation (Milton, 1988). On the other hand, complete cover is very costly and may not be feasible in the time available: investigators must usually be satisfied with a sample. So now, if we wish to estimate some property of vegetation everywhere, we need to interpolate. The usual situation will be as follows. We shall have measured at least one property of the vegetation, such as the green leaf area index (GLAI, one-sided area of green leaves per unit of ground) at some sites on the ground. We shall have measured the radiation from the vegetation at many more sites, including those where we measured the vegetation. Our task then is to estimate the properties of the vegetation from the radiation both at the site where we have measured it and at intermediate sites. We could do it by first interpolating and then using a regression. If we kriged the radiation, as Webster et al. (1989) did, we should use our knowledge of the spatial correlation in the radiation optimally, but we should not take into account any spatial dependence in the primary variables. If we are to make the best possible use of our limited ground data, then we must consider the spatial dependence of the two sets of variables and their interdependence simultaneously. The technique for doing this is cokriging, and it is based on a model of coregionalization (Myers, 1982; 1984; Journel, 1989). The cokriging itself is readily automated, but the modeling requires understanding and judgement. Once a satisfactory model has been found it may be used to study the consequences of particular sampling strategies and to plan sampling to achieve given precision optimally, that is, with the least effort to achieve a given precision. In this article, we describe acceptable models

of coregionalization and how to select them, the cokriging formulation, and the way in which they may be applied to design sampling strategies for surveys using a ground-based radiometer. We illustrate these procedures with three examples. COI:tEGIONAIJZATION AND COKRIGING

Coregionalization is embodied in the more general theory of regionalized variables (Matheron, 1965; 1971). Essentially, it concerns the relations between two or more spatially correlated random functions (regionalized variables). Each of the variables is assumed to be continuous and to take a value z(x) in a space of one, two, or three dimensions and to be a realization of a random function Z(x) at that point. In the present context the space is two-dimensional. The spatial correlation in a single regionalized variable may be expressed in the variogram, the function that relates the variance of the variables at two points to their separation: 3'(h) = 1AE[[Z(x) - Z(x + h)]2],

(1)

where the vector h is the separation in distance and direction, the lag. If the variable is second-order stationary, i.e. has both constant mean and variance, then the spatial covariance function exists: C(h) = E[IZ(x) - #] [Z(x + h) - #/],

(2)

where # is the mean of the variable. It is related to the variogram by C(h) = C(0) - v(h),

(3)

where C(0) is the covariance at zero lag, the a priori variance a 2 of the random variable. These functions are finding use in remote sensing (e.g., Atkinson and Danson, 1988; Atkinson et al., 1990; Atkinson, 1991), and we have already described them in this journal (Curran, 1988; Webster et al., 1989) as have others (e.g., Woodcock et al., 1988a,b; Smith et al., 1989; Cohen et al., 1990; Bhatti et al., 1991). Rather than dwell on them here we want to extend the ideas to two or more variables. Suppose we have two variables, u and v. By extension of the above we can define a cross variogram: vn,(h) = I/2E[Z,(x) - Zu(x + h)/[Z,)(x) - Z~(x + h)/l,

(4)

Cokrigingwith Ground-BasedRadiometry 47 where Z,(x), Zv(x), Z.(x + h) and Zv(x + h) are the values o f u and v at places x and x + h, respectively. Similarly, we can consider the cross covariance function: C~o(h) = E[[Z.(x) - t~uJ[Zo(x + h) - #v]].

(5)

T h e cross variogram is related to it by

n

r = ZX,z(x,),

(9)

i=l

w h e r e the z(xi) are the values of Z at x, and the hi are any arbitrary weights attributed to them. The sum Y is itself a random variable with variance n

23,,~(h) = 2C.~(0) - C,v(h) - Co,(h)

~kAjC(xi,xj),

var[Y] = ~

(10)

i=lj=l

= 2C.o(0) -

C.o( -

h) - C~.( - h).

(6)

Again, this exists provided both variables are second-order stationary. However, unlike the cross variogram, it is not usually symmetric: C,~(h) only incidentally equals Co,(h). Any substantial difference b e t w e e n the two would m e a n that one variable lags b e h i n d the other, in some direction, and in these circumstances the coregionalization can be described fully only by the cross covariance function. Small observed differences b e t w e e n C,~(h) and Co,(h) are likely to be sampling effects. They can be ignored, and the cross variogram may then be used to express the coregionalization. T h e usual formula for estimating the semivariance at any particular lag, h, is 1 m(h) ~(h) = 2m(h-----)i=IZ [Z(Xi) -- Z(Xi "}" h)} 2, (7) where the z(xi) are the observed values at sampling points x. i = 1, 2 . . . . . n, and m(h) is the n u m b e r of paired comparisons at lag h. By varying h we can obtain an o r d e r e d set of values for "~(h), the sample or experimental variogram. T h e cross variogram can be estimated similarly by 1 ,,(h) ~uv(h) = 2m(h) i=~Z{z,(xi) - z,(xi + h)] x/z~(x,) -

zv(xi+ h)].

(8)

Note that both u and v must have been m e a s u r e d at x~ and xi + h to use this formula. The Linear Coregionalization

Model

Provided that we have m e a s u r e m e n t s of all the variables at enough positions, we can c o m p u t e the sample variograms using Eqs. (7) and (8). To p r o c e e d further, we must choose and fit models that will both describe t h e m reasonably faithfully and are acceptable as variance functions. For a single variable, Z(x), we can create from it another variable Y as the linear sum o f n values:

where C(xi,xj) denotes the covariance of Z between the points x, and xj. Clearly, var[Y] cannot be negative, and the covariance function in Eq. (10) must be such as to guarantee this: It must be positive semidefinite. However, since a variance of zero is most unusual, acceptable functions are usually said to be positive defnite. Given the equivalence of the covariance and semivariance [Eq. (3)], we ~an express the variance of Y in terms of semivariances: tl

var[Y] =

n

rl

C(O)~X,EXj- ~ EkikD,(x,xj). (11) i=l

j=l

i=lj=l

If the variation is u n b o u n d e d , we must express the variance of Y in this way since the covariances are undefined; and we must also choose weights that sum to zero so that the first term on the right h a n d side of the equation disappears. This leaves us with vartY] = - ~, ~hihjT(x,,x#).

(12)

i=lj=l

The variance will be positive provided that the variogram is negative definite. If we allow for zero variance, then the variogram must be conditional negative semidefinite (CNSD), the condition being that the weights sum to zero. There are only a few simple functions that are CNSD (Webster and Oliver, 1990). These are now well known, and McBratney and Webster (1986) list those most frequently encountered when describing environmental data. The same restricted set of models applies to cross variograms. In addition, however, w h e n we choose models for two or more coregionalized variables, we have also to take into account linear combinations of the variables. Any such combination is itself a regionalized variable, and its variances must also be positive or zero. This is assured as follows. Each variable Z(x) with which we are con-

48

Atkinson et al.

cerned is assumed to be a linear sum of K orthogonal random variables, Y(x), whose variograms we may denote by qck(h), k = 1, 2 . . . . . K. In the linear model these combine to p r o d u c e the variograms for any pair of variables u and v as

where nt is the n u m b e r of sampling points for variable l. To avoid bias, the weights must sum as follows: nl

~)xa=l

K

~,=(h) = ~]~v3,~(h), k=l

and

where the superscript k is simply an index, not a power. In this expression b~v =/~w for all values of k, and for each k the matrix of coefficients

bkb ] must be positive definite. Since the matrix is symmetric, it is sufficient that ~ . i> 0 and b~ >t 0 and that its determinant is positive or zero:

Ib l = Ib21 < ~/bk,~,

forl--u,

i=1

(13)

(14)

which is Schwarz's inequality. If there are V variables in coregionalization, then the full matrix of coefficients, [~], will be of the order V, and its determinants and all its principal minors must be positive or zero. Schwarz's inequality means that every basic c o m p o n e n t or variogram, "yk(h), appearing in a cross variogram, 3'u,(h), must also be present in the two auto variograms 3',,(h) and 3~v,(h). However, the reverse is not necessary. A c o m p o n e n t may appear in the two auto variograms without being present in the cross variogram. Thus, we may have b~v = 0, where b~, > 0. In practice, it means that the actual variograms chosen to describe the coregionalization must be compatible in the above sense and the matrices [b~] must be checked for positive definiteness. In the research described below we made this check.

nl

~ h a -- 0

for all 1 ¢: u.

i=1

The first condition means that there must be at least one observation of u for cokriging to be possible. Otherwise, unlike Eq. (8), where we c o m p u t e d the cross variogram, there are no constraints on w h e t h e r the different variables are recorded. Subject to these conditions, we minimize the estimation variance: a~(B) = E[[z,(B) - ~,(B)} 2].

This leads to equations of the form V

Z

nl

Z "~/v(X//,Xjv) + Cv = ~uv(B,xjv),

(16)

l=l i=l

for all v = 1 to V and all j = 1 to nv. Here ~/tv(xa,Xjv) is the cross semivariance b e t w e e n variables l and v at sites i and j, and 3~,v(B,xjv)is the average cross semivariance b e t w e e n the two variables b e t w e e n the block B and site j. The quantity ¢v is a Lagrange multiplier for the vth variable. There is one such equation for each combination of variables and sampling points, and in addition there are V equations for the nonbias constraints. These equations are solved to find the weights h~l, which are then inserted into Eq. (15) to obtain the estimate of zu(B). In addition, we obtain the cokriging variance as v

nt~

a~(B) = Z ~,hj~uv(B,xjv) - ¢, - ~u~(B,B), (17) v=lj=l

Local Estimation The cokriging equations for estimating one variable from a coregionalized set of variables are simply extensions of those for autokriging, which we described earlier in this journal (Webster et al., 1989). The notation, however, is unavoidably more complex. Let there be V variables, l = 1, 2, . . . . V, and suppose that we wish to estimate one of them, u, in a block B. Our estimate will be V

nl

= Z ZX z(x ), l=1i=1

(15)

where ~uu(B,B) is the within-block variance of variable u. The above equations are for blocks of size B. If B is r e d u c e d to the support of the sample at a point x0, we perform punctual kriging. The average ~,~(B,x~) becomes simply 2(,0 (x0,x~), and ~ ( B , B ) , the within-block variance, disappears. Myers (1982) and McBratney and W e b s t e r (1983) give the matrix formulation of these equations, and examples in the environmental sciences are provided by Vauclin et al. (1983), Aboufirassi and Marifio (1984), and Leenaers et al. (1989).

Cokriging with Ground-Based Radiometry 49

In practice, we are likely to have many measurements of radiation and fewer measurements of the variable that really interests us, u. We choose and fit a suitable model and check that it satisfies our rule. We can then use it and our data to estimate local values of u at as many places as we wish within the region sampled. The kriging is automatic once we have made our choices of model and block size, and there are now many examples reported in the literature (Webster and Oliver, 1990). Suppose, however, that we have a model for coregionalization, either from a preliminary reconnaissance or from an earlier investigation of a similar kind and that we have yet to perform our survey. Can we use the model to design a sampling scheme that will give tolerable estimates with least effort (and will be in this sense optimal), in the same way as we did for a single variable (Webster et al., 1989; Atkinson, 1991)? Briefly the answer is "yes," and we devote the remainder of this article to it.

REMOTE SENSING EXAMPLES

Summarizing the procedure in Webster et al. (1989) for a single variable, the maximum kriging variance, O~ma~,is USually minimized by sampling on an equilateral triangular grid (provided that variation is isotropic), although a square grid gives variances that are almost as small and in general is preferred. For punctual kriging tr~, always occurs at the centers of grid cells whereas for block kriging it can occur either where the blocks are centered over the centers of the grid cells or over the grid nodes. The quantity OZu~, is computed for several sampling intervals and plotted against the sampling interval. The designer decides what maximum value of ~m~x is tolerable and then reads from the graph the corresponding sampling interval. A sampling grid with this interval is optimal. The same principle applies to two or more coregionalized variables. Two factors complicate the choice of strategy, however. The first is that not only do the sampling intensities of the variables vary but their intensities differ from one another. The other complicating factor is that for variable u, the one we wish to estimate, the maximum kriging variance, au~m~x,does not neces-

sarily occur at the centre of the sampling configuration as it does where we have just the single variable. The position of ~m~ in relation to the primary and subsidiary grids depends on: i. the form of the variograms and the cross variogram(s); ii. the strenth of the cross correlation; iii. the sampling interval of the primary grid; iv. the relative intensities of the primary and subsidiary grids. For square grids with sampling ratios of 2 and 4, the positions of o ~ x lie on the diagonals of the grids. For larger ratios there are other possibilities, and Figure 1, taken from McBratney and Webster (1983), shows these. For block kriging, there are even more possible solutions. For the two smallest grid ratios tr~um~xcan be found by computing the the kriging variance along the relevant diagonals. For all other sampling grid ratios O'2max c a n be found only by calculating the kriging variance at the nodes of a very fine grid within a representative quadrant of a primary grid cell. The maximum variance of £,(xo) or £,(B) thereby obtained can be plotted against sampling interval, and the optimal sampling strategy can be determined by reading from the graph as before. However, in the multivariable case, we have still not obtained the final solution because there may be more than one combination of variables and sampling configuration: Which we choose will depend on the costs, and the final choice will be the one that achieves the desired precision for least cost. To illustrate the above procedures, we describe three experiments using measured radiation to predict the GLAI of barley and dry biomass and percentage cover of clover in pasture. In each experiment, five grid ratios (1:1, 1:4, 1:9, 1:16, and 1:25) were compared. The number of observations used for each estimate was determined by selecting all observations within a search radius of two lags of the primary variable with a maximum number of secondary observations for any variable set to 80. For all grid ratios 12 pr/mary observations were used and for grid ratios of 1:1, 1:4, 1:9, 1:16, and 1:25, 12, 48, 80, 80, and 80 subsidiary observations were used. Our first experiment was at Broom's Barn Experimental Station, lying 5 km to the west of Bury St. Edmunds in Suffolk, England. The crop was spring barley growing in a field approximately

50 Atkinson et al.

•

•

I I I I

&%

• / \ \

•

x

~

•

x

x

I

•

•

•

•

•\

•

•

x

• AN

x

x

•

•

•

•

• /

x

\ e\

•

X

• X

•

/

"to' /

• \

" x

•

•

x

•

/

"

•\

•

•

"

%x

"

x

Main and subsidiary variables observed

t

Subsidiary variables observed

.

.

• X

/•

/ .

•

•

/N 0/

g

/

o\/o

•

•

•

%.

• X •

•

/•

•

.

• /

•

•

.

•

•

oN

•

\

•

•

•

" \ %

\

x

X "

•

\

X

•

•

X

/0%

• X

\ •

/

•

\ • /

•

/

~

/o /

•

/•

/

/•"

\

/

•

\ •N

/

• x

x

X

•

"•

/ N

•

I I :

a ( e r o p ) ~--

\

•

•

I

x

\

/

x

II • ~--~•

• . . . .

\

/

'

•

/ / % / / \

/ I I

measurements w e r e made over a Kodak grey card (reference) with a second sensor for conversion to reflectance. The reflectance (R) of the crop was then c o m p u t e d from these data by the equation

/•

\

.

"•

Figure 1. Places of maximum estimation variance for sampling schemes on square grids with different intensifies for the subsidiary variable. Maximum variance can occur on the pecked lines and near the crosses (McBratney and Webster, 1983).

450 m x 200 m on almost level ground. Its soil is a rich sandy loam overlying chalky drift known to vary appreciably over distances of 5 - 2 0 m. On 2 May 1987, two sets of measurements were made on the crop. O n e h u n d r e d contiguous plots, each 1 m x 1 m, were demarcated with dark brown strings and arranged in line to form a transect. The transect was orientated north-south, parallel to the drills of the crop. The relative radiance (Lr) was measured using a Milton multiband radiometer (Milton, 1988) over the center of each 1 m 2 plot, giving an instantaneous-field-of-view ( I F • V ) of 0.78 m. Three measurements were made in rapid succession and averaged. Simultaneous

[Zr(crop)a(refe . . . . . )l /ILr( refe. . . . . . ) / ,

(18)

w h e r e the reflectance of the grey card was determined using a G E R infrared intelligent spectroradiometer (IRIS) calibrated to the British national standard. The measurements of relative radiance were made b e t w e e n 11:37 h and 12:55 h Greenwich Mean Time (GMT). The whole of the barley above ground was cut from each plot and weighed. The area of the green leaves was then measured using a Licor AT leaf area meter and used to calculate the GLAI. Our second experiment was performed on pasture at Morrell Wood, 3 km east of Belper in Derbyshire, England. It lies some 100 m above sea level on a gentle slope. The site had been reclaimed after coal mining, and the sward was a recently sown mixture of rye grass (Lolium perenne, Lolium multiflorum) and clover (Trifolium repens). Cattle and sheep grazed the field intermittently throughout the year. W e sampled this field on two occasions. On 6 May 1988 we laid out a 200 m transect of 100 1 m x 1 m plots. The grass was about 10 cm high, but it had not b e e n grazed since the previous autumn. On 6 August 1988 we chose a shorter transect of 100 m, but with 200 20 cm x 20 cm plots located regularly. This time the field had b e e n grazed, but otherwise the conditions were similar to those on 6 May. Using the same kind of radiometer as above, we measured the relative radiance once over the center of each 1 m 2 plot from a height of 2 m; there was bright sunshine and a clear blue sky on both dates. Simultaneous measurements were made over a Halon panel so that the relative radiance could be converted to reflectance using Eq. (18). Similarly the reflectance of the Halon panel had been measured by an IRIS to provide reflectance calibrated against the UK national standard. These measurements were made between 11:12 h and 12:02 h GMT on 6 May 1988 and 10:55 h and 11:37 h GMT on 6 August 1988. On 6 May 1988 we harvested the crop, weighed it, and determined the dry biomass by drying and reweighing a sample. This value was used to estimate the dry biomass for each plot.

Cokriging with Ground-Based Radiometry 51

On 6 August 1988 color positive photographs were taken of each 1 m ~ plot. These were projected onto a square grid of 156 dots within a circle that represented the IFOV of the sensor. The occurrences of the dots in clover and grass were counted and converted to the percentage cover of clover. To study the eoregionalization, we converted the reflectance measured in all three experiments into the normalized difference vegetation index (NDVI) (near i n f r a r e d - r e d ) / ( n e a r infrared + red), and used this as the subsidiary variable for cokriging.

18

•

16

Z

14

• •

•

•

•

•

dl

q

•

•

12 10-

'~ .~

8-

>

64-

GLAI

2o

i 2

i 4

i 6

i 8

i 10

i 12

i 14

i 15

i 15

i 20

22

26

Lag (m) 4.0

3.6

E s t i m a t i n g t h e G L A I at B r o o m ' s B a r n o n 2 M a y 1987

The experimental variograms and the cross variogram representing the coregionalization of the GLAI and NDVI for 2 May 1987 at Broom's Barn were computed using Eqs. (7) and (8) (Fig. 2). We fitted several nonlinear models by weighted least squares approximation with the weights proportional to the numbers of paired comparisons using Ross's (1987) program, MLP. The variogram of the GLAI and the cross variogram (GLAI × NDVI) appeared to have similar forms and could be fitted well by a spherical model with nugget: 3'(h) = Co+ ~2a

_

)

2\a/

9

N~ ~ o ~v°2.5 v

~ 2.0,~ ,.~~ ,.oi

0.5-

0.o

GLAI x NDVI i

i

I

I

I

i

i

i

i

2

4

6

8

10

12

14

16

18

i 20

i 22

i 24

i 28

Lag (m)

30~ zsx

for 0 < b ~< a,

~zoCD O

"y(b) = Co+ c

.§ ,.5- / /

for b > a,

>

-r(o) = o.

(19)

These models are shown as continuous lines in Figures 2a-c. The ranges of the spherical model from the best least-squares solution were 18.4 m (GLAI) and 18.23 m (GLAI x NDVI). To convert these into a single model of coregionalization, a common range of 18 m was chosen and the models refitted (Table 1). The variogram of the NDVI was also fitted with this model using a range fixed at 18 m (Fig. 2c). The fit appeared acceptable, although the experimental semivariance at a lag of 1 was overestimated. Therefore, a second spherical component with a short range, 2.76 m, and a sill variance of 0.802 was introduced for the variogram of the NDVI (Fig. 2c) to increase the accuracy of the fit. The coefficient matrices of these functions were positive definite.

.',.oE

co

0.5-

NDVl

00

i

i

i

i

i

i

i

i

i

=

i

2

4

6

8

10

12

14

16

18

20

22

i 24

i 26

Lag (m)

Figure 2. Variograms of GLAI and NDVI for Broom's Barn data of 2 May 1987 for GLAI, GLAI × NDVI, and NDVI with two models.

The coefficients of this model were inserted into Eqs. (16) with the right-hand sides replaced by "yuv (xou, x~) to obtain the punctual cokriging variances. Figures 3a shows the resulting graphs of maximum kriging variance against sampling interval of the principal variable, GLAI, at five ratios with the subsidiary variable NDVI. These ratios are 1:1 (effectively univariate kriging), 1:4,

52

Atkinson et al.

Table 1. Models of Variograms and Cross Variograms for Estimating the GLAI at Broom's Barn on 2 May 1987 Variables

Model

CLAI x 104 GLAI x NDVI x 10 ~ NDVI x 104

Spherical Spherical Spherical

Double spherical

co

cl

al

10.8 0.534 1.67 1.0

7.53 2.43 0.932 0.786

18.0 18.0 18.0 18.0

1:9, 1:16, and 1:25. The graphs show how the precision of estimation increases not only with increased sampling of the GLAI but also with increasing intensity of the subsidiary NDVI. As above, we can use the graphs to choose a sampling scheme to achieve some specified tolerance. For example, if a tolerable cokriging variance was set at (18 %)2, then the sample spacing of the primary grid could be increased from 14 m to 22 m with a grid ratio of 1:4, to 29 m for a grid ratio of 1:9,

e2

a~

Maximum Lag

2.76

25 25 25

25

0.802

to 36 m for a grid ratio of 1:16, and to 45 m for a grid ratio of 1:25. Bear in mind that these intervals are in one dimension for a two-dimensional grid, and so they represent quite substantial decreases in sampling effort for the primary variable. Of course, the investigator might want more precision. If a maximum tolerable estimation variance were (15%) 2, then the investigator would have to sample very intensely and so could not

(a)

20-

x

N o

15-

O O ¢-t'~ I,,t~ >

Grid ratios 10-

1:1 1:4 1:9

..........

1:16

~.--

1:25

......

C -e-

S-

6

0

0 0

I

I

I

I

I

10

20

30

40

50

Spacing

6-

x

(m)

(b)

5-

~o 4O ¢t~ •z.--

3-

c-

2-

f

.~,..'"

~¢'Y-

(5 I

I

I

I

10

20

30

40

Spacing ( m )

Figure 3. Maximum cokriging variance against sample spacing for mapping the G L A I at Broom's Barn on 2 May 1987 with ground-based radiometry: a) punctual cokriging and b) cokriging over 10 m × 10 m blocks.

Cokriging with Ground-Based Radiometry 53

expect to diminish the intensity of the primary grid much by cokriging. The nuggelf variance of the primary variable sets a lower limit to the estimation variance, and the investigator must accept a realistic limit to what can be achieved. In this instance, the nugget variance is 10.8 (GLAI × 104)3, and the punctual estimation variance cannot be less than this. As implied in the theoretical section, an investigator may wish to estimate the average value of a property over areas larger than the support of the sample. To do this, the area is chosen and Eqs. (16) are used as they stand. The size of the area depends on the nature of the investigation. To illustrate what happens, we chose 10 m × 10 m squares, and we repeated the calculations for such blocks. Figure 3b shows the results. The variances of block estimates are much less than those for points. This is because the large nugget variance in the variogram of the GLAI is always encompassed within a block. The shape of the graph is much the same, however. Figure 3b has two features that we mention here and that recur in Figures 5b and 7b. First, at the shortest spacings of the primary variable the maximum kriging variance appears to increase as the sampling of the subsidiary variable intensifies, and the lines on the graph cross as the spacing increases. This is not what you would expect. In principle, the kriging variance must decrease as the amount of information on which estimation is based increases. However, because of computing limitations we have to restrict the n u m b e r of subsidiary sampling points, and as the spacing narrows so the points are concentrated increasingly around the center of the block being estimated. This may cause the maximum estimation variance to increase in two ways. First, the observations of the subsidiary variable may become so concentrated around the centre that they do not cover the block evenly. For example, at a grid ratio of 1:25 and a primary grid spacing of 5 m, the 10 m × 10 m block is not evenly covered by the subsidiary variable. Second, whether the estimate is made at a point or over a block, the concentration of observations near the location of the estimate causes some duplication of information, which may result in a small increase in the maximum estimation variance. Second, the asymptotes of the lines on the graph are different. This is caused by the differ-

ence in the n u m b e r of subsidiary observations that are used at each grid ratio. For example, at a very large grid spacing, say 1000 m, the grid ratio of 1:1 (12 primary observations and 12 subsidiary observations) will produce a maximum cokriging variance greater than that produced by a grid ratio of 1:4 (12 primary observations and 48 subsidiary observations), even though all the observations contain uncorrelated information. The estimation variances are not the only factors that determine a sampling strategy. Different kinds of observation require different amounts of effort, and these differing costs must b e t a k e n into account so that the investigator can choose the combination that achieves the desired precision for least cost. Suppose it takes approximately 30 min to take 100 measurements of relative radiance (two personnel) and 4-5 h to harvest the crop (four personnel). The measurements of relative radiance must be converted to reflectance and the NDVI calculated, but this is quick once the necessary calibration equation is available (no more than a few hours). By contrast the GLAI is time-consuming to determine, and it would take one person 2-3 days to obtain 100 observations. In all, the. GLAI is around two orders of magnitude more expensive to estimate than the NDVI. However, as the n u m b e r of measurements increases, the relative costs of setting up equipment, traveling to each observation point in the field and programming decrease because of economies of scale. A more subjective estimate of the difference in cost may be more realistic, therefore. Judging from experience (e.g., Curran and Williamson, 1985), the ratio in cost between radiometric and GLAI measurements should be closer to 1:20. There are more pragmatic considerations. The maximum n u m b e r of observations of the GLAI that can be made before leaves start to decay is about 100 if only one leaf area meter is used. A sample of 100 is equivalent to a sample spacing of 30 m. From Figure 3 it is evident that in these circumstances cokriging is the only means by which to reduce the error. However, at a grid ratio of 1:25, 100 observations of the primary variable correspond to 2500 measurements of reflectance. It is not possible to make 2500 measurements of reflectance (12.5 h according to the earlier estimate) in 1 day using a single radiometer.

54

Atkinson et al.

It is not easy to estimate the cost of measurements. Nevertheless, a position on the graph (Figs. 3a and b) must be chosen, and this position depends on the relative costs of measuring the primary and subsidiary variables. In practice, it will be a matter of judgement for the investigator. For example, if we wished to estimate the values of the GLAI over 10 m × 10 m blocks with a O'2maxof no more than (4 % )2, then the primary grid spacings for the grid ratios in parentheses would be 14 m (1:1), 19 m (1:4), 25 m (1:9), 29 m (1:16), and 36 m (1:25). If the relative cost of measurement is taken to be about 1:20, then an optimal sampling strategy may be found by computing the total cost at each grid ratio (Table 2). The table shows how the use of subsidiary variables, in this case NDVI, decreases to around a third the total sampling cost. Estimating the Dry Biomass at Belper on 6 May 1988 Choosing models to represent the variograms of dry biomass on 6 May 1988 at Belper was straightforward. The spherical model was the most suitable with ranges of 16.9 m, 9.3 m, and 12.9 m when fitted separately. We compromised by substituting a common range of 13 m and refitted the model (Figs. 4a-c) (Table 3). The coefficient matrices of these functions are positive definite. The maximum cokriging variances were computed using the coefficients of these variograms for both points and 10 m x 10 m blocks, and the results are displayed in Figure 5. As for Broom's Barn, we can explore the consequences of these graphs for planning sampling. Suppose that the investigator wishes to produce estimates and can

afford a primary sample spacing of 15 m. If punctual univariate kriging is used, O~umaxis 1726 (g/ m2)2. By cokriging with a grid ratio of 1:4, O~umaxis reduced to 1492 (g / m~)2. Increasing the grid ratio to 1:9 reduces o ~ x further to 1415 (g/m2) 2, and at a grid ratio of 1:16 O~umaxis 1365 (g/m2) 2. The cost of measuring the secondary variable has not been accounted for, however. To do this, we need to know the relative cost of measuring the two variables. As for the GLAI in the previous example, the primary variable, the dry biomass, may be taken to cost 20 times as much to measure as the secondary variable, NDVI. The cost of univariate kriging at a 15 m grid spacing may be evaluated as the relative cost multiplied by the size of the sample (20 x 400 = 8000), which is the maximum affordable. Clearly, if cokriging is to be useful, we must be more precise for the same cost. At a grid ratio of 1:4 the same cost is achieved by sampling with a primary sample size of 334, equivalent to a primary grid spacing of ~/270 = 16.4 m. At this primary grid spacing and with a grid ratio of 1:4, the cokriging variance is about 1480 (g / m~)~. The cokriging variances may be evaluated similarly for grid ratios of 1:9, 1:16, and 1:25 (Table 4). Evidently, when the cost of measuring the secondary variable is taken into account, cokriging at grid ratios greater than 1:4 provides only a small increase in precision for the same cost. Estimating the Percentage Cover of Clover at Belper on 6 August 1988 Both the variogram of the NDVI and the cross variogram (percentage cover of clover × NDVI) were regular and similar, and the linear model

Table 2. R e l a t i v e C o s t o f S a m p l i n g for E s t i m a t i n g t h e G L A I at P o i n t s U s i n g F i v e G r i d Sizes at B r o o m ' s B a r n o n 2 M a y 1 9 8 7 a Primary Secondary Grid Grid Primary Secondary Grid Spacing Spacing Sample Sample Ratio (m) (m) Size Size 1:1 1:4 1:9 1:16 1:25

14 19 25 29 36

-9.5 8.33 7.25 7.2

460 250 144 108 70

-998 1296 1713 1737

Total Cost Primary Sample (Secondary + 20 × Primary Sample Size Size) × 20 9200 5000 2880 2160 1400

9200 5998 4176 3873 3137

a Primary samples are GLAI, secondary samples are NDVI, and m a x i m u m cokriging variance is set to 4%.

Cokriging with Ground-Based Radiometry

55

18.

% X

II

16-

1412-

~,m i

°~

4-

GO

2-

Dry Biomass i

I

,

4

6

8

I

I

I

I

i

i

i

10

12

14

16

18

20

22

Lag (m) 16O 14>,¢ A

12-

.ei

6-

4-

Dry Biomass x NDVI 0

Lag (m) 30-

•

•

•

25-

t-t¢1 >

15-

10-

i

NDVl . 2

.

.

. 4

~ 6

8

1

~ 1

, 14

¢ 6

,

i

Figure 4. Variograms of dry biomass and

i

NDVI for Belper data of 6 May 1988 for dry biomass, dry biomass × NDVI, and NDVI.

18

Lag(m)

Table 3. Models of Variograms and Cross Variograms for Estimating the Dry Biomass at Belper on 6 May 1988 Model

Co

cl

(m)

Maximum Lag (m)

Spherical Spherical Spherical

6.82 4.23 15.8

9.42 9.21 13.3

13.0 13.0 13.0

18 20 15

a

Variables Dry biomass x 10 -2 Dry biomass x NDVI x 10 NDVI × 104

56

Atkinson et al.

(a)

20O

15E

Grid ratios C ~ O~

e-

10-

1

1

1:4

5-

1:9

..........

1:16

-----

1:25

......

6

0

0

I 10

0

I 20

I 30

I 40

I 50

Spacing (m)

6O

(b)

x / I

5,/ f f

E .i.

0

tt~

. ,.=..-- • . .----*"

=........==.. - ~ "

. ......

.--'"

//

3-

/ > iT) c-

..-"I

................

...-" . ...'"'........-"""

..:",~-'"" ..",fir" ""

2-

-t'0

0

I

I

I

[

I

10

20

30

40

50

Spacing (m)

with a sill component provided the best fit (Table 5). By contrast the variogram for the percentage cover of clover was best described by a pentaspherical model. The circular model provided the most accurate compromise between these two models with ranges of 12.8 m (NDVI) and 13.3 m (cross variogram). Therefore, a range of 13 m was chosen and the circular model refitted to

Table 4. Relationship between Grid Size and Punctual Cokriging Variance for Estimating the Dry Biomass at Belper on 6 May 1988

Grid Ratio

Primary Grid Spacing (m)

Cokriging Variance (g / m2)2

1:1 1:4 1:9 1:16

15 16.4 18 20.1

1620 1480

1:25

22.5

1470 1450 1430

Figure 5. Maximum cokriging variance against sample spacing for mapping the dry biomass at Belper on 6 May 1988 with ground-based radiometry: a) punctual cokriging and b) cokriging over 10 m x 10 m blocks.

all three variograms (Figs. 6a-c) (Table 5). Each variogram was represented satisfactorily with the same basic model as in intrinsic coregionalization (Journal and Huijbregts, 1978). The coefficient matrices of these functions are positive definite. The variograms were then used to compute the maximum cokriging variances for both points (Fig. 7a) and 10 m x 10 m blocks (Fig. 7b). In this example cokriging provides substantial savings in sampling effort. Figures 7b shows a large difference in cokriging variance between the curves representing the grid ratios of 1:1 and 1:25. Assuming a maximum tolerable cokriging variance of (13%) 2, grid ratios of 1:1, 1:4, 1:9, 1:16, and 1:25 would require grid spacings of 11 m, 21 m, 31 m, 41 m, and 50 m, respectively, which represent sample sizes of 744, 204, 94, 56, and 36. Thus, the size of sample for the primary variable (percentage cover of clover), which is

Cokriging with Ground-Based Radiometry

57

Table 5. Models of Variograms and Cross Variogramsfor Estimating the Percentage Cover of Clover at Belper on 6 August 1988a co

Cl

(m)

Maximum Lag (m)

37.3 -1.53 37.9 - 2.44 0.802

39.3 26.1 23.1 37.3 29.0 23.1

18.7 11.2 11.1 13.0 13.0 13.0

18 25 16 18 16 16

a

Variables

Model

% cover (dover) % c o v e r x N D V I x 102 N D V I x 104 % c o v e r (clover) % cover x NDVI NDVI

Penta-sph. Lin. + sill Lin. + sill Circular Circular Circular

N o t e t h a t t h e % c o v e r x N D V I cross v a r i o g r a m was c o n s t r a i n e d to pass t h r o u g h the origin.

Figure 6. Variograms of percent clover cover and NDVI

tedious and time-consuming to determine, can be reduced from 744 to 36 by adding (25 x 36--) 900 measurements of NDVI. The cost of these sampling strategies may be evaluated as in Table 6 if the relative cost of measuring the percentage cover of clover is taken to be the same as that for the GLAI in the preceding example, that is, 20 times that of the NDVI. Clearly cokriging at a grid ratio of 1:25 is the most economical solution given that the relative cost of measurement is 1:20. For this example, cokriging is a very attractive alternative to univariate kriging.

for Belper data of 6 August 1988 for percentage cover of clover, percentage cover of clover x NDVI, and NDVI.

70-

20 10

Percentage Cover of Clover

0

i 1

i 2

, 3

i 4

i 5

i 6

J 7

i 8

i ~

i 10

, 11

i 12

i 13

J 14

i 15

~ 16

J 17

i 18

Lag (m)

DISCUSSION AND CONCLUSIONS

x

i

15" 10-

ver of Clove x NDVI 0 f

1

2

3

4

$

6

7

8

9

10

11

12

13

14

15

16

, 11

i 12

i t3

i 14

i 15

J 16

17

18

Lag (rn) 28.

24.

1816'

•

12. 10.

E

8.

4"

NDVI

2. 0

i 1

i 2

i 3

i 4

i 5

i 6

i 7

i 8

t 9

1 10

Lag (m)

i 17

i 18

The first problem that must be faced when attempting to economize on sampling using cokriging is to model each of the variograms within the constraints of the linear model of coregionalization. It might not be possible to do this, and a solution is not guaranteed. Second, we may be able to choose and fit a suitable set of models, but these might not enable us to improve on the precision that we could obtain by univariate kriging. It depends where the envelope bounded by univariate kriging and cokriging at large grid ratios lies in the space of the graph of variance against spacing. Provided that the above problems can be solved, however, cokriging can often provide an economical method for estimating a variable of interest from a sparse sample of that variable and a denser sample of reflectance. In our examples cokriging resulted in a precision of around three times that achievable with univariate kriging for the same effort (see example 1; Table 2). This gain in precision comes from the extra information

58

Atkinson et al.

100 "

(a)

80-

e-

~

6o

°~

~ 4o

1:4

~

1:16

----

1:25

......

1:9

2o

0

Grid

i

I 20

10

ratios

..........

i

t

q

30

40

50

Spacing (m)

30

(b)

v /

2o j

>

I

O) r-.

..'"°

t

IO

.~1

°,-°"

.i

°°,°"

.o.'"

°o..-

. .,.~. •

..-'°

. I ..~-°

.

~.---°°

Figure 7. Maximum cokrigingvari-

..,n~,~.~t..T~... . . . . . . . . . . . . . . .

.~ v,

6

0

o

'

10

'

2'0

3'0

I

,0

Spacing (m)

provided by the observations of reflectance. When the cost of measuring the secondary variable is taken into account, the gain for a given effort is reduced, but not by much because the difference in the cost of measuring the two variables is so large. Similarly the saving in effort provided by cokriging to a given precision was as

s'0

ance against sample spacing for mapping the p e r c e n t a g e cover of clover at B e l p e r on 6 August 1988: a) punctual cokriging and b) cokriging over 10 m x 10 m blocks.

much as nine times that achievable with univariate kriging (see example 3; Table 6). That is, if cokriging with ground-based radiometry were used in place of kriging with only the primary variable, then the investigator would require as little as one-ninth of the total effort. The data for our analysis were, of course,

Table 6. Relative Cost of Sampling for Estimating the Percentage Cover o f Clover over 10 m x 10 m Blocks Using Five Grid Sizes at Belper on 6 August 1988"

Grid Ratio

Primary Grid Spacing (m)

Secondary Grid Spacing (m)

Primary Sample Size

Secondary Sample Size

Primary Sample x 20

Total Cost (Secondary + 20 x Primary Sample Size)

1:1 1:4 1:9 1:16 1:25

11 21 31 41 50

-10.50 10.33 10.25 10.00

744 204 94 56 36

-816 843 857 900

14880 4080 1880 1120 720

14880 4896 2723 1977 1620

° Primary samples are percentage cover of clover; secondary samples are NDVI, and maximum cokriging variance is set to 13%.

Cokriging with Ground-Based Radiometry

one-dimensional, and we used them to design sampling in two dimensions. Thus, we proceeded as though variation were isotropic, which we have every reason to believe was the case at our working scale both at Broom's Barn and Belper. It would be equally justified in many agricultural scenes elsewhere. However, it is also very easy to extend the principles to situations where there is geometric anisotropy. The analysis is performed for the direction in which the variogram has the shortest range or steepest gradient to obtain the spacings in that direction. The spacings in the perpendicular direction are obtained from them simply by multiplying by the anisotropy ratio. Burgess et al. (1981) describe the procedure fully. Interpolation is necessary to map a variable of interest at the ground from a sample of that variable. Kriging does this optimally in the sense that it estimates intermediate values with minimum variance. However, even with the best sampling scheme, this minimum estimation variance might be too large to be acceptable. It might be diminished by incorporating observations of a subsidiary variable reflectance, using the multivariable extension ofunivariate kriging, cokriging. In all of the examples presented here the precision of estimating the variable of interest was reduced greatly by adding observations of reflectance. Even when the cost of measuring this reflectance is taken into account, cokriging can still offer substantial increases in precision over univariate kriging for a given amount of effort. In one of our examples above this amounted to a ninefold saving in effort compared with univariate kriging. Cokriging coupled with ground-based radiometry provides an economical and operational technique for using reflectance when estimating properties of the Earth's surface. We advise an investigator who wishes to map a variable that is likely to be related to reflectance to compute the variograms for the variable of prime interest and the reflectance and their cross variogram. Then, by cokriging, the investigator will be able to design an optimal sampling strategy to achieve a desired precision. Most of this research was undertaken while P. M. Atkinson was at the University of Sheffleld and Rothamsted Experimental Station, R. Webster was at Rothamsted Experimental Station, and P. J. Curran was at NASA Ames Research Center (California) and the University of Sheffield. We thank the Natural

59

Environment Research Council for its CASE award to P. M. Atkinson; the Leverhulme Trust for its research award to P. M. Atkinson and our colleagues who supported and assisted us at Rothamsted and the University of Sheffteld, notably, Mrs. P. M. Ashcroft, Professor J. A. Catt, Dr. E M. Danson, Mrs. J. Munden, Mr. R. Patel, and Dr. S. E. Plummer.

REFERENCES Aboufirassi, M., and Marifio, M. (1984), Cokriging of aquifer transmissivity and specific capacity, Math. Geol. 16:1936.

Atkinson, P. M. (1991), Optimal ground-based sampling for remote sensing investigations: estimating the regional mean, Int. 1. Remote Sens. 12:559-567. Atkinson, P. M., and Danson, E M. (1988), Spatial resolution for remote sensing of forest plantations, in Proceedings, IGARSS "88 Symposium, ESA Publications Divisions, ESTEC, Noordwijk, The Netherlands, pp. 221-223. Atkinson, E M., Curran, P. J., and Webster, R. (1990), Sampling remotely sensed imagery for storage, retrieval and reconstruction, Prof. Geographer 42:345-353. Bhatti, A. U., Mulla, D. J., and Frazier, B. E. (1991), Estimation of soil properties and wheat yields on complex eroded hills using geostatistics and Thematic Mapper images, Remote Sens. Environ. 37:181-191. Burgess, T. M., Webster, R., and McBratney, A. B. (1981), Optimal interpolation and isarithmic mapping of soil properties IV. Sampling strategy, J. Soil Sci. 32:643-659. Cohen, W. B., Spies, T. A., and Bradshaw, G. A. (1990), Semi-variograms of digital imagery for analysis of conifer canopy structure, Remote Sens. Environ. 34:167-178. Curran, P. J. (1988), The semi-variogram in remote sensing: an introduction, Remote Sens. Environ. 24:483-507. Curran, P. J., and Williamson, H. D. (1985), The accuracy of ground data used in remote sensing investigations, Int. J. Remote Sens. 6:1637-1651. Journel, A. G. (1989), Fundamentals of Geostatistics in Five Lessons, American Geophysical Union, Washington, DC. Journel, A. G., and Huijbregts, C. J. (1978), Mining Geostatistics, Academic, London. Leenaers, H., Burrough, P. A., and Okx, J. P. (1989), Efficient mapping of heavy metal pollution on floodplains by cokriging from elevation data, in Three Dimensional Applications in Geographical Information Systems (J. Raper, ed.) Taylor and Francis, London, pp. 37-50. Matheron, G. (1965), Les Variables R~gionalis~es et leur Estimation, Masson, Paris. Matheron, G. (1971), The Theory of Regionalized Variables and Its Applications, Ecole des Mines de Paris, Fontainebleau. McBratney, A. B., and Webster, R. (1983), Optimal interpolation and isarithmic mapping of soil properties V. Coregionalization and multiple sampling strategy, J. Soil Sci. 34:137-162.

60 Atkinson et al.

McBratney, A. B., and Webster, R. (1986), Choosing functions for semivariograms and fitting them to sampling estimates, J. Soil Sci. 37:617-639.

Vauclin, M., Vieira, S. R., Vachaud, G., and Nielsen, D. R. (1983), The use of cokriging with limited field soil observations, Soil Sci. Soc. Am. J. 47:175-184.

Milton, E. J. (1988), Principles of field spectroscopy, Int. J. Remote Sens. 8:1807-1827.

Webster, R., and Oliver, M. A. (1990), Statistical Methods for Soil and Land Resource Survey, Oxford University Press, Oxford.

Myers, D. E. (1982), Matrix formulation of cokriging, Math. Geol. 14:249-257. Myers, D. E. (1984), Cokriging-new developments, in Geostatistics for Natural Resources Characterizations (Part 1),

(Verly, G., David, M., Journel, A. G., and Marechal, A., eds.), Reidel, Dordrecht, pp. 295-305. Ross, G. J. S. (1987), Maximum Likelihood Program, Numerical Algorithms Group, Oxford. Smith, R. C. G., Pratehrapar, S. A., Barrs, H. D., and Slavich, P. (1989), Use of a thermal scanner image of a water stressed crop to study soil spatial variability, Remote Sens. Environ. 29:111-120.

Webster, R., Curran, P. J., and Munden, J. w. (1989), Spatial correlation in reflected radiation from the ground and its implications for sampling and mapping by ground-based radiometry, Remote Sens. Environ. 29:67-78. Woodcock, C. E., Strahler, A. H., and Jupp, D. L. B. (1988a), The use of variograms in remote sensing I: Scene models and simulated images, Remote Sens. Environ. 25:323-348. Woodcock, C. E., Strahler, A. H., and Jupp, D. L. B. (1988b), The use of variograms in remote sensing Ih Real digital images, Remote Sens. Environ. 25:349-379.