National Academies Press: OpenBook

Spatial Statistics and Digital Image Analysis (1991)

Chapter: 3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis

« Previous: 2. Image Analysis and Computer Vision
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 37
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 38
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 39
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 40
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 41
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 42
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 43
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 44
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 45
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 46
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 47
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 48
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 49
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 50
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 51
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 52
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 53
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 54
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 55
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 56
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 57
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 58
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 59
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 60
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 61
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 62
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 63
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 64
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 65
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 66
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 67
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 68
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 69
Suggested Citation:"3. Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis." National Research Council. 1991. Spatial Statistics and Digital Image Analysis. Washington, DC: The National Academies Press. doi: 10.17226/1783.
×
Page 70

Below is the uncorrected machine-read text of this chapter, intended to provide our own search engines and external engines with highly rich, chapter-representative searchable text of each book. Because it is UNCORRECTED material, please consider the following text as a useful but insufficient proxy for the authoritative book pages.

3 Oceanographic and Atmospheric Applications of Spatial Statistics and Digital Image Analysis James ]. Simpson Scripps Institution of (oceanography 3.l Introduction Historically, remote sensing of the environment has helped process-oriented studies examine individual aspects of the physics, chemistry, and biology of the earth-ocean-atmosphere system. From such studies, several problems of global significance have emerged that are cross-disciplinary in nature. Ex- amples include acid rain, the increase in atmospheric carbon dioxide, antici- pated depletion of the ozone layer, ET Nino-related mollifications in weather and ocean circulation with their resultant effects in agriculture en c! fisheries, and tropical rain forest destruction by fires of human origin. It has been rec- ognized that the key to progress on these and other cross-clisciplinary issues in earth science during the decade of the 1990s probably will be addressing those questions that concern the integrated functioning of the earth as a system (EOS, 1984~. The hyclrologic cycle, the biogeochemical cycle, and climate processes are the fundamental processes that integrate the earth as a system, and thus each of these cycles must be examined in detail and on a global scale if meaningful progress is to be made on the problems cited above. I,arge-scaTe, synoptic observations of a winkle variety of phenomena (e.g., sea surface temperature, sea level, winct stress, ozone concentration, radi- 37

38 ation heat balance, sea ice, vegetation, and near-surface current velocity) are required to study these cycles properly, and remote sensing provides the only practical way to collect synoptically the necessary data on the global scale. The magnitude of the data set, however, may preclude understanding unless it can be distilled and synthesized into organized patterns that can be related meaningfully to the underlying governing physics, chemistry, or biology of these global cycles and processes. Spatial statistics and math- ematical methods of digital image analysis provide mechanisms for such a synthesis. This chapter discusses a few areas of spatial statistics and data represen- tation useful in digital image analysis that have immediate application in the area of remote sensing of the earth. Some special considerations needed for the correct analysis of remotely-sensed data are presented and a few critical research areas are identified. 3.2 Selected Analysis Areas In this section three mathematical methods useful in digital analysis of se- quences of remotely sensed images are presented. These methods were cho- sen because they have a broad range of applications in earth science. These methods are non-Bayesian in nature; a presentation on Bayesian methods used in image analysis is given in chapter 2 of this report. Several abbreviations commonly used in remote sensing and the earth sciences appear in this chapter. Each is defined at the place it first occurs in the text. For easy reference, there is an appendix containing abbreviations and their definitions at the end of this chapter. 3.2.] Principal Component Analysis General Concepts Principal component analysis (PCA) is a multivariate statistical technique that can be applied to all forms of multispectral or multi-temporal image data and is most commonly applied in the general context of arbitrary mul- tivariate data. Forms of PCA have been used to study pattern classification (Geladi e! al., 1989), sequential segmentation (Esbensen and Geladi, 1989), spatial patterns of variability in sea surface temperature and phytoplankton pigment observed in satellite data (LagerIoef, 1986), and in algorithms for cloud removal from satellite data (e.g., GalIaudet and Simpson, 1991a).

39 The principal component transformation (PCT) is the basis for all these analyses. The PCT is a linear transformation that isolates uncorrelated linear combinations of a given set of variables in such a way that each element in this combination represents a decreasing amount of variance in the original variables (Ingerbritsen and Lyon, 1985~. This linear transformation defines a new set of coordinate axes for the data such that (13 the transformed origin is at the mean of the data distribution (LiDesand and Kiefer, 1987), (2) the transformed coordinate axes are mutually orthogonal (Jensen, 1986), and (3) the transformed coordinate axes are in the directions of maximum variance (Jensen, 1986~. Below, two examples of PCA are developed, one using multispectral image data and the other using multi-temporal image data. An extensive treatment of the use of PCA techniques in atmospheric sciences and oceanography is given by Preisendorfer (1988~. Multispectral Data Application The ability to accurately and automatically segment clouds in remotely sensed imagery is critically important to a broad range of disciplines in earth science. Clouds significantly affect the net heating of the atmosphere and the underlying ocean-land surface by modifying solar and terrestrial radiation (Ohring and CIapp, 1980~. This net radiative heating governs the thermody- namics and dynamics of the atmosphere, which in turn influence the forma- tion and dissipation of clouds (e.g., Matveev, 1984~. The potential feedback effects associated with this cloud-radiation interaction are among the great- est sources of uncertainty in determining the relation between changes in external conditions such as solar radiation and atmospheric carbon dioxide concentration and changes in climate (e.g., Henderson-SeDers, 1982; Ra- manathan, 1987~. Clouds also affect our ability to remotely sense the prop- erties of the atmosphere, ocean, and land; such observations are needed, for example, in weather prediction (e.g., Pailleux, 1986), oceanography (e.g., Eckstein and Simpson, 1990a,b), and agriculture (e.g., rainfall, Browning, 1986~. The PCT can be used with multispectral image data to robustly segment clouds from natural images. An example of such an application is given below. The Principal Component Transformation Split-ancI-Merge Clus- tering (PCTSMC) Algorithm. Here, the development of the PCTSMC algorithm is presented in abbreviated form. Mathematical details of the al- gorithm are given in Gallaudet and Simpson (1991a). An AVHRR infrared

40 image, designated T. consisting of brightness temperatures calibrated to de- grees Celsius in bands 3, 4, and 5 (hereinafter referred to as T3, T4, Ts) is first differenced to construct a 2-banded differenced image D (i.e., T3—T4, T3—Thy. This differenced image now contains aD the information of the original infrared image needed for cloud detection, but only in two linearly independent bands. Next, this differenced image is transformed using the PCT. The result is a 2-banded transformed differenced image in which all interband correlation is destroyed (Masher, 1987), thus making it a Togi- cal preprocessor to the segmentation operation (step 3) of the PCTSMC algorithm. Step 3 of the PCTSMC algorithm performs image segmentation using a split-and-merge clustering procedure (e.g., PavIidis, 1977; Seddon and Hunt, 1985; Richards, 1986) on the PC transform of the di~erencecl image. This results in a segmented image in which the natural spectra] classes in the original image are separated into distinct groups (i.e., land versus ocean versus cloud). The method of clustering that is employed in the PCTSMC algorithm combines both the partitional and hierarchical approaches. It consists of a partitional clustering algorithm augmented by a splitting-and- merging step at each iteration. Combining a partitional with a hierarchical method has several advantages over the use of either method alone: (1) pure hierarchical methods are not appropriate for complex data (MuerIe and Allen, 1968; Fukada, 1980; Jain and Dubes, 1988~; (2) pure hierar- chical methods are more appropriate for data that is to be partitioned on the basis of both local and global information, rather than global informa- tion only plain and Dubes, 1988~; (3) pure hierarchical methods impose a taxonomic structure on the data (Anderberg, 1973), which is not character- istic of cloud-containing AVHRR imagery; (4) pure hierarchical methods are order dependent i.e., the resulting segmentation will vary ~lepending upon the order in which the regions are split and merged (Cheevasuvit et al., 1986~; this often results in less than optimal segmentations of the data; (5) pure partitional algorithms often converge to local minima of the clus- tering criterion function (Pairman and Kittler, 1986; Jain and Dubes, 19SS); (6) the combined approach is more efficient than pure merging or pure splitting methods of region detection (Paviidis, 1977; Richards, 1986~; and (7) the combined approach is less dependent on the initial segmentation, and therefore is more capable of recovering all of the natural clusters in the data (Seddon and Hunt, 1985; Jain and Dubes, 1988~; this is because the number of clusters in the initial segmentation need not be the same as those O ~ ~ is

41 that actuaby exist in the given data (Seddon and Hunt, 1985; Pairman en c! Kittler, 1986). In the final step of the PCTSMC cloud screening procedure, the data are retransformed back to the feature space, and each of the clusters is labeled as either a cloud or non-cloud. In the first three steps discussed above, the operations performed on the data were entirely unsupervised- i.e., no a priori knowledge was required. Hence, they apply to an arbitrary image segmentation. In this fourth step, expert knowledge is introduced to perform a boolean classification. Rules appropriate for land versus cloud versus ocean separation in AVHRR image data are given in Simpson and Humphrey (1990) and GaHaudet and Simpson (1991a). Plate 3.1a shows AVHRR Band 2 data; clouds appear as white or gray tones. The coastline is white in this panel. Plate 3.1b shows AVHRR Band 4 infrared temperature; coldest temperature is white and warmest tempera- ture is black. In this panel, the coastline is black. The segmented image produced by the PCTSMC algorithm is shown in Plate 3.1c and the final cloud-masked sea surface temperature is shown in Plate 3.16. In this final panel, the warmest ocean temperature is white, cooler ocean temperatures appear as shades of gray, and cloud contaminated pixels and land are shown as black. Land was masked from the Plate 3.1 images using a recursive polygon fill algorithm (Simpson, 1991) and is rendered either white or black depending on the gray scale mapping used in the individual panel. Multi-Temporal Data Application A major objective of earth science studies is to identify spatial patterns of variance in temporal sequences of images. Examples include the analysis of variability in sea surface temperature structure in oceanic current systems (e.g., Lagerioef, 1986) and in seasonal and interannual variation in phyto- plankton abundance (e.g., Strub et al., 1990~. The form of PCA used in such studies is generally referred to as empirical orthogonal function (EOF) analysis: "empirical" because the functions arise from the data themselves and "orthogonal" because they are uncorrelated. Note that closed-form mathematical functions generally cannot be used as the basis functions rep- resenting the complex images observed in nature. Empirical orthogonal functions are useful in work with large data sets. The method separates a data set D(x,t) into spatial components Fi(X) and

42 temporal components Alit) such that NT D(x, t) = ~ Fi(X) · Ai(t) i=1 (3.1) where ~ is the number of EOFs computed. Note that x and t repre- sent generalized spatial and temporal coordinates, i.e., x represents the set (x~,...,xArx) and ~ represents the set {ti,...,t~. In matrix notation, D(x,t) is an NX x IN matrix, with NX representing the number of posi- tions x and ~ the number of time steps t. Each column of the F matrix is an EOF. Each EOF has NX values, and there are JVT EOFs. Thus, F is an NX x JVT matrix. The time series matrix A has ~ rows and JVT columns. \ Forms of Normalization. The computations may be implemented with two different normalization schemes because the EOF representation decom- poses the space-time series, using separation of variables, into a sum of prod- ucts of temporal amplitudes, Ain't), modulating spatial patterns of variance, Fi(X). Note that, in physics, separation of variables occurs widely because of the form of the underlying differential equations. Here, the motivation is to provide a compact statistical representation of the data in which spatial patterns in variance can be distinguished from temporal patterns. Either the Apt) or the Fi(X) can retain the same physical units as the original data. Historically, a normalization scheme (method 2) has been used in which the temporal components of the decomposition retained the same units of the data. More recent studies in oceanography and atmospheric science (e.g., Barnett and Patzert, 1980) prefer to use EOFs with the same units as the original data for easier interpretation of the spatial patterns of variance re- sulting from the EOF analysis (method I). Both methods of normalization are equivalent and both are given here for completeness. (The summation notation is used for consistency with the vast majority of published studies.) In the first method Fi(X) · Fjord) NX ~ Fi~xk)Fj (Xk) = Aitij, (3.2, k=1 where Fi and Fj are columns of the F matrix and {ij is the Kronecker delta function. The coefficient Ai is an eigenvalue of the covariance matrix C = D D , where DT is the transpose of D. The rows of the time series

43 matrix are orthonormal, i.e., Ai(t) Aj(t) NT = JVT NT ~ Ai(tk)Aj(tk) = ttj ~ (3.3) where Ai and Aj are rows of A. The EOFs then give the variance in the same units as the data. The second normalization method forces the EOFs to be orthonormal, ye., Hence' Fi(x) Fj(x) = bij. TF F—~ _ , (3.4) (3.5) where ~ is the identity matrix. The rows of the time series matrix are orthogonal: Ai(t) Aj(t) = I. (3.6) With this normalization, the EOFs no longer have the same units as the data. Theoretical Basis. The objective of EOF analysis is to represent a given matrix of data D by the product F A. In the discussion that follows, the second method of normalization is used, and thus the matrix F satisfies equation (3.5~. The equations governing EOF theory can be derived from the eigen- equation of the covariance matrix, C · F = F · A , (3.7) where C is the NX x NX covariance matrix = D D , A is the diagonal matrix of eigenvalues, and F is the matrix of EOFs. Because C is a real symmetric positive definite matrix, the elements of A are positive. By (lefi- nition of the covariance matrix C, (3.7) can be rewritten as { ~ } (3.8) The eigenvectors in the matrix F are now used to define the matrix A = FT · D. This is the principal component transform. Thus, F · FT = ~ and necessarily D = F.A. (3.9)

44 At this point the number of EOFs is equal to JVT. The purpose of EOF analysis, however, is to produce a representation of the data that is more compact than the original data set. For this goal to be met, it is necessary for most of the variance (65 to 80~) in the original data set to be contained in the first few EOFs. (By convention, the eigenvalues are ordered with the largest being first.) If this is true, then F may be truncated to produce a new matrix F having just NT columns, where hT ~ Ad. Then D is the approximation to D computed from the EOF decomposition: D~D=F A. (3.10) The decomposition is the most efficient representation of D with regard to a mean square error criterion (Davis, 1976~. Finally, it should be noted that if D contains a noninteger number of cycles of a sinusoidally varying variance mode, then the variance represented by the associated eigenvalue may not agree with the actual variance. For example, the variance of a sine wave differs for one-half and one full cycle, but the covariance matrix C = D D is the same. Hence, the EOF will return the variance for an entire cycle when the data represent one-half cycle. A Computational Method. The empirical orthogonal functions can be computed from a singular value decomposition (SVD) of the data D. When NX > By, then (Press et al., 1986) D = U . W . VT (3.11) where U is an NX x NX orthogonal matrix with only NT linearly indepen- (lent columns, W is an NX x IVY matrix with an 1\7 x HI diagonal upper portion having positive or zero elements, and V is an NT x ~ orthogonal matrix. Note that the diagonal elements in the IVI x NT upper portion of W contain the singular values of D. Cad this portion of W the matrix W'. The U and V matrices are orthogonal in the sense that their columns are orthonormal, that is, TU U = I art Vl V = I. (3.12) Using this decomposition, the covariance eigen-equation can be written as C U= I NT } U = ~ MU · W · VT) {V · WT · UT) · U. (3.13)

45 which simplifies to C.U= ~ u.W2 NT . (3.14) This is equivalent to the eigen-equation (3.7) with F—the highest order columns of U and ~—W2/~, using the second method of normalization. Hence, the diagonal elements Wi of W' relate to the eigenvalues Ai of the covariance matrix C via the equation ~ = A. (3.15) Thus, the SVD returns the EOFs as the highest order columns of U. and A = W'~VT Generalizations to Two Dimensions. The EOF decomposition is not confined to space-time series where space is one-dimensional. If the posi- tions are actuary (x, y) coordinates, such as latitudes and longitudes or lines and samples from images, then if each time step covers the same locations, D(x,y,t) can be reduced to D(x',t). The sequence of positions x' is con- structed by setting x~ = (at, yip, X2 = (X2, yI), · - -, XNX = (XNX, YI ), · ~ ~ . XNX.NY = (XNX, ANY ). If one considers x and y to be rows and columns of an individual time step matrix, then x' is equivalent to concatenating the rows of D together. Whether one concatenates the rows or the columns is imma- terial. (Note, EOF analysis with two-dimensional images is actually three dimensional: x, y, and t. To render the problem tractable, the spatial di- mensions must be concatenated so that one can work with a two-climensional space-time matrix.) Then the individual EOFs Fi are also vectors in x' and, in order to map the patterns of variation associated with the EOFs, this vector must now be parsed back into its rows and columns. Example. The purpose of this example is to show that EOF analysis can be useful for determining the dominant patterns of spatial variance in a sequence of images. (The normalization scheme used in this example is that of method 1.) For this purpose, a sequence of eight images was constructed by superimposing an image of an inclined plane with that of a disk that is out of phase with the inclined plane. Note also that the range of data in the image of the inclined plane is about twice that of the image of the circle. One-half cycle of the image sequence so constructed is shown in Plate 3.2a. The EOF decomposition of this sequence (Plate 3.2b) identifies

46 two dominant patterns of spatial variance in the test sequence. EOF #1 is an inclined plane that accounts for 67.7% of the total variance in the sequence. EOF #2 is a circle that accounts for 32.2~o of the total variance. The corresponding time series for EOF #1 and EOF #2 (Plate 3.2b) correctly establish the phase relationship between the two patterns. Additional Considerations. In remote-sensing applications, multiple ob- servations are generally spread over different times. Moreover, data are gen- erally not evenly distributed either in time (e.g., due to variations in orbit) or in space (e.g., cloud cover may obscure some pixels in the scene). These circumstances can bias the results of an EOF analysis because the EOF analysis is predicated on evenly distributed data in both space and time. If the NT images in the data set D (equation (3.1~) are not equispaced in time, the resultant EOF decomposition may be biased toward specific time periods. A practical way to minimize such temporal bias is to construct a weighting scheme for the images in the data set. Generally, temporal weights are constructed in such a way as to preserve the total mean and to- tal variance of the data while simultaneously minimizing the temporal bias (e.g., Kelly, 1985~. Spatial data (i.e., one or more of the NX pixels in a given image) are often replaced by compositing data from other images that surround the given image in time. Care must be taken to ensure that the composite time scale is very much less than the time step between images in the sequence undergoing EOF analysis. Other interpolation schemes (e.g., kriging) can also be used to minimize spatial and temporal bias in EOF analyses resulting from imperfectly sampled data. There are three ways to compute EOFs: (1) large covariance matrix approach, (2) small covariance matrix approach, and (3) singular value de- composition. The first two methods involve a direct solution of the ei~en- -- r ~ o- ~ . · ~ . ~ · . · 4. . ~ ~ . . ~ rams value equation of the covarlance matrix of the data set L'. Mere are two ways to do this because the data in the image sequence can be stored in the two-dimensional data array in two different ways. If the data are stored such that there are NX rows and NT columns, the covariance matrix will be NX x NX. If, however, the data are stored such that there are NT rows and NX columns, then the covariance matrix will be NT x NT. For most remote sensing applications, the number of spatial points NX in a given image in the sequence is usually very much greater than the number of im- ages, NT, in the sequence. Thus the method of data storage giving rise to an NX x NX covariance matrix is often caned the large covariance ma-

47 trix approach and generally cannot be used in remote sensing applications. Likewise, the method of storage resulting in an NT x JVT covariance matrix is generally caped the small covariance matrix approach. It has been used often in oceanographic and atmospheric applications (e.g., Preisendorfer, 1988~. Solution of the eigenvalue equation of the covariance matrix by sin- gular value decomposition was discussed earlier in this chapter. A detailed discussion of the various methods of solution of the eigenvalue equation of the covariance matrix, in the context of EOF analysis, is given by GalIaudet and Simpson (199lb). 3.2.2 Velocity Estimates from Image Sequences General Comments AD methods of motion estimation based on image sequence analysis depend in some way on the detection of image brightness gradients. These gradients are defined as normal in direction to contours of constant brightness (or sea surface temperature [SST] for the case of the Advanced Very High Resolution Radiometer [AVHRR] flying on the NOAA series of operation satellites). The total velocity v at any point on a contour can be written as v = zJn+rt, (3.16) where ~ and T are the magnitudes of the normal and tangential velocity components, respectively, and n and t are the unit vectors in the directions normal and tangential to the contour, respectively. The total velocity vector can also be decomposed into ordinary Cartesian coordinates v= u'+v', (3.17) where u and v are the magnitudes of the x and y velocity components, respectively, and ~ and ~ are the unit vectors in the x and y directions, respectively. Below, objective methods for computing the total, normal, and tangential components of near-surface oceanic flow from sequences of AVHRR data are presented. The formulation follows closely that of Wahl and Simpson (199Oa,b). Note that atmospheric motions can be computed using these same methods from sequences of either AVHRR or Geostationary Operational Environmental Satellite (GOES) data.

48 Total Flow Fielc! From Pattern Matching Consider a pair of real, discrete, two-dimensional random functions, sum, n) and Pam, n). The expectation value of these random functions can be ap- prox~mated as the normalized sum over the tile occupied by the random function in the image. Thus, NL NS E[(s(m, n)] = Nip x NS ~ ~ s(m, n), (3.18) m=1 n=1 where NT and NS are the number of lines and the number of samples in the tile. Assuming stationarity of their first two cross moments, the auto- covariance and cross-covariance of the two functions are given by Css(mO'n0~) = 13~(s(m, n) - r~S)(s~m + mO' n + nO) - risen (3.19) Csp~mO, no) = Et~stm, n)—~S)(p~m + me, n + nO)—Opel, (3.20) respectively, where Et.] is the expected value, (mO,nO) is the spatial lag between the functions, and Is and lip are given by the function means 71s = Ets(m, n)] up = E[p~m, n)] . (3.21) (3.22) It can be seen that the function variances ~2 and up are the zero lag auto- covariances. The correlation coefficient is defined as rSp~mO' no) = csp~mo, nO) (3.23) Used such that rummy, note < 1 rSS(O, 0) = 1 . (3.24) (3.25) If the second signal is an exact spatially lagged version of the first signal, s(m, n) = p(m + mO, n + no), (3.26) then equations (3.20) and (3.25) require that the correlation coefficient achieve an absolute maximum value at this lag, or rSp(mO,nO) = 1. For any physical signals' the maximum correlation will be a value less than 1 since the second signal is not necessarily an exact lagged version of the first.

49 ~ ~. ~— ~ A detailed discussion of the two-dimensional cross-correlation function is given by Dudgeon and Mersereau (1984~. Now consider two consecutive satellite SST images mapped to the same spatial grid. These images can be thought of as two-dimensional discrete functions. Given an image subsection skm,n) from the second image, the problem is to determine if it contains a region similar to a smaller subsection finis from the first image. Let p(~,j) be caned the pattern tile, and let skm, n) be called the search tile. The pattern tile is a section of the first image that occupies the same spatial coordinates as the center region of the search tile in the second image. Then, the correlation matrix between the pattern and search tiles is given by k ~ — ski ~,j [ski + k' j + l) - tr75(k, [~] split j) - 71p~J rapt, ~— THE,' ~,jEsf i + k' j + i) _ 77S(k, b12 Hi Hi, j) _ spit / ~ ~ /O ~ (3.27) where r/5(k, I) is the average value of Dim, n) in the subregion coincident with pair j), and the summations are over the coordinates common to both s and p. The value of no is computed once outside the summations and is given — Or — _ ~ ,~ CAN ~ . ~ . . . · ~ ~ ~ ~ , ~ ~ , · ~ ~ ~ — ov ~ ~.-~-2 1. where the pattern tile vi m. n ~ replaces the search tile stm, n) in (3.18~. The ranges of k and ~ correspond to the regions of correlation in which phi, j) is completely contained in Him, n). The pattern matching method determines the spatial lag between the ~ ~ - ,, ~ , ~ , , ~ ma. . . . ~ . ~ 0 pattern tile from the first image and the search tile in the second image by finding the location of their maximum correlation. With this spatial lag and the time between the images, the average velocity of the features in the pattern tile can be computed. The most basic assumption of the method! is that the spatial displacement of thermal gradient features can be tracked as if the shape were time invariant. This assumption would be rigorously true if the correlation (3.27) were equal to 1 for each pattern-search tile pair. Unfortunately, this condition is never met. Hence, it becomes necessary to determine an acceptable minimum correlation or correlation threshold (see Wah! and Simpson, 1990a). Velocities obtained from the pattern matching technique and proper choice of correlation threshold show good correspon- dence with observations. The pattern matching method may also be implemented in the wave number domain. This implementation uses the discrete Fourier transform (DFT) property s(<m, n`) * prim, A) ~ ~ S'f~k=, ky)P(k=, ky`), (3.28)

50 where * represents the correlation operation, S'(k=,ky) is the conjugate of the DFT of s(m, n), and P(k=,ky) is the DFT of p(m, n). This is the case because the correlation in the space-time domain and the spectrum in the frequency-wavenumber domain are Fourier transform pairs (Dudgeon and Mersereau, 1984). Thus, the correlation between two discrete signals can be computed by taking the product of the DFT of a zero-padded pattern tile with the conjugate of the DFT of the search tile and taking the inverse DFT of this product. It may seem that computing the DFT of the discrete func- tions via the fast Fourier transform (FFT) is a more efficient approach. In most applications, however, the pattern tile usually occupies a much smaller area than the search tile. If the number of non-zero terms in the pattern tile is less than 132, it is more efficient to implement (3.27) than to use the FFT algorithm to compute the correlation function (Campbell, 1969~. (For odd size tiles the Winograd implementation of the FFT is recommencled.) The method used here is similar to other template matching schemes, such as area correlation and matched filtering (Jain, 1989~. Matched filtering involves the construction of a linear filter that maximizes the output signal- to-noise ratio. Using the matched filtering technique, the area surrounding the pattern is assumed to be colored noise. If the power spectral density of the noise is known, the signal-to-noise ratio can be maximized by passing the signals through a high-pass filter before performing the area correlation. Minimum distortion methods have been used to do interframe registra- tion in video camera systems (Jain and Jain, 1981~. Simpson and Bloom (1990) have applied this method to the computation of near-surface velocity from sequences of images and have shown that the distortion is simply re- lated to the correlation. In effect, maximizing the correlation is equivalent to minimizing the distortion assuming that the variance and standard devi- ations of the different search tiles remain the same. Both methods yield the same velocity fields (Simpson and Bloom, 1990) for a given image sequence. However, the minimum distortion method executes faster because it does not require a standard deviation computation. The basic assumption of these methods is shape invariance of the pat- tern under translation. Rotational motion of the pattern, however, also can occur. Pattern matching techniques that detect rotational motion and/or combined translational-rotational motion of the pattern have been devel- oped (e.g., Jain, 1989~. These methods are computationally very expensive but win perform well in the presence of curvature in the motion field of the pattern. Only translational techniques were used herein because they are sufficient to determine the general flow pattern in the image sequence under consideration.

51 The Normal Component of Flow Marr and UlIman (MU) Method. Early visual primitives can provide clues to establish the motion of elements in a visual field (Marr and UlIman, 1981~. The simplest such primitives are the image raw intensity values, but these provide no information about the shapes of objects. The next higher- order primitives are the zero crossing segments produced by the convolution of an image with the Laplacian of Gaussian (LOG) operator, V2G (Marr and Hildreth, 1980~. The LOG operator is defined by V2G = k [1 - ~ 202 ] exp [ 2~2 , (3.29) where x and y are the number of rows and columns from the function center, a is the Gaussian width parameter, and k is a normalization constant. The parameter ~ determines the spatial scales of intensity changes detectable by the V2G operator. The V2G operator is the optimal smoothing bandpass filter in the sense that it minimizes the product of bandwidth and spatial localization (e.g., Marr and Hildreth, 1980~. If ~ is the demeaned image function, then locations of zero crossings of the convolution of the V2G operator with ~ wiD correspond to locations of intensity changes (i.e., gradi- ents). Let this convolution be denotes! V2G * ~ = I', where I' is the output of the convolution. Note the units of ~ and I' are the same (°C for AVHRR data). Note also that the LOG operator (a commonly used edge detector) produces a set of zero crossings in the image. The locations of the zero crossings (i.e., the edges) are then determined by a zero crossing operator that detects the positions of the edges by finding the locations of changes in sign in the LOG-convolved image. The idea of directionally sensitive units, which can establish the direc- tion of movement of an edge detected by the V2G operator, was introduced by Marr and UlIman (1981, hereafter referred to as MU) to determine the motion of visual elements. Given an approximation to the time derivative of I' and the spatial rate of change of I', the direction of motion of a zero crossing in either the line or sample direction can be cletermined. The To- cations of zero crossings in the first convolved image of the sequence are determined by a zero crossing operator. Then, at each zero crossing pixel, the normal component of flow can be estimated using the equation __O_ a VI' I (3.30)

52 where the subscript represents differentiation with respect to time, and the unit normal vector n is in the direction of the gradient V]'. Equation (3.30) is the conservation of heat equation (= = 0) for the convolved image I'. If the magnitude of the gradient is negative and It' is positive, then motion of the edge is in the positive x direction. If at iS negative, motion is in the negative x direction. The opposite is true if the magnitude of the gradient is positive. Thus, the negative sign in (3.30) gives the correct direction of the normal component of flow. It can be seen that a directionally sensitive unit is represented by a transition in the sign of I~ combined with the sign of the gradient. Spatial-Scale Considerations. Methods used to compute the normal component of flow employ small neighborhoods, typically 4 x 4 pixels or smaller in size. The total velocity is computed over a much larger neighbor- hood. Typically, pattern tile sizes vary between 16 x 16 and 32 x 32 pixels. Thus, estimates of the total velocity via pattern recognition represent the mean motion of the centroid of the pattern as measured by the displacement of the two-dimensional cross-correlation maximum of the pattern. The nor- mal velocities, however, provide local estimates of the displacement of small spatial-scale gradient, typically computed over 4 x 4 tiles. This basic differ- ence in spatial scales further constrains the computation of the tangential flow. Other Representations Optical Flow (OF) Met hocI. Optical flow (hereafter referred to as OF) is an estimate of the motion of solid bodies based on a first-order variation of brightness patterns in an image (e.g., Horn and Schunck, 1981). This method of computing the velocity field from a sequence of images is based on the solution of two constraint equations. The first constraint equation relates the velocity in an image to the image brightness (or temperature) pattern and is called the "motion constraint equation": DT Dt =Tt+V-VT=o' ,... `. - (3.31) where Dt is the material derivative operator, T' is the partial derivative of brightness (or SST) with respect to time, and v is the total velocity vector. This equation can only estimate the normal component of velocity because the tangential component solution is an annihilator (i.e., aft ETA.

53 The form and physical interpretation of the second constraint used in OF methods is one that has been widely debated. In general, it seems that the second constraint imposed is not based on any physical characteristics of the flow field, but, rather, it is a mathematical constraint imposed to pro- duce a unique solution. Horn and Schunck (1981) used a constraint based on the smooth variation of the flow field to derive an iterative scheme for computing what they referred to as the "total flow." Various other schemes have been introduced to estimate the total flow (see Aggarwal and Nand- hakumar (1988) for a review of OF methods). It is noted (e.g., Horn and Schunck, 1981; Hildreth, 1983; Verri and Poggio, 1989) that the estimate of the total flow field may be very far from the actual velocity field, depending on various factors influencing the time series of images. Given the material derivative constraint (3.31), one assumes that the flow is continuous and varies smoothly over small spatial scales. One way to define the measure of smoothness is to exarn~ne the squares of the mag- nitudes of the spatial rate of change of the OF velocity. This can be written as a departure from smoothness error E2 = U2 + U2 + V2 + V2 (3.32) where (u,v) is the local total velocity vector, and the subscripts represent differentiation with respect to the spatial coordinates (x, y). There wiD also be errors in the estimation of the partial derivatives of brightness because noise is amplified by differentiation. Thus the equality of (3.31) will not be exact. Define this error term as Eb = Tt + uTx + vTy, (3.33) where subscripts indicate differentiation with respect to either a spatial (x, y) or temporal (t) coordinate. The objective function to be minimized can be written as the integral J = / i(Eb + ct2E2) a?x cly, (3.34) where Eb is the error in computing the material derivative, Ec is the measure of smoothness, and Ct2 iS a weighting parameter. The calculus of variations can be used to m~nimize this integral. Then the variation equations can be rewritten as spatial iterative equations Uk+} = Uk T=(Txuk + Tyvk + Tt) Ct2 + TX2 + Ty2

54 k+1 -k _ Ty(T-~u + Ty + ), (3.35) where the superscript k is an iteration index. The first step is to com- pute the brightness derivatives at Al points in the image using centered finite differences. Then, starting with an initial estimate of zero velocity, the method spatially iterates on the velocity values until velocity residual between estimates is small. The resulting velocity field is then used as the new initial velocity estimate for the next time step if more than one time step is available. It is interesting to note that updated velocity values in the iteration equations (3.35) do not rely solely on the previous values at a given point, but rather on the local averages of velocity. Note that the local averages of velocities typically are computed over small spatial domains and are computationally efficient. The parameter ar2 is seen to be important in regions of small gradient. If the gradients are small relative to O2 then (X2 wiD dominate any perturbations in the estimation of the derivatives at this point. Minimum Norm (MN) Solution for Normal Flow. Equation (3.31) gives one equation for the two unknowns (u,v) of the total velocity. This underconstrained system has fewer equations than unknowns and thus has an infinite number of solutions. One way to solve such underconstrained systems is to find the solution with the minimum vector length, or norm (Luenberger, 1969~. Equation (3.31) was solved for the normal component of velocity using the singular value decomposition to obtain the solution of minimum norm (hereafter referred to as MN) at every point in the given image subsection (Wah! and Simpson, 1990b). It can again be seen that the solution yields only the normal component of flow because the tangential component is an annihilator of (3.31~. The MN solution of normal velocity was done on a point-by-point basis on raw temperature data using a centered finite difference for the derivatives of temperature. Equivalence of Minimum Norm (MN) ant} Optical Flow (OF) for C`2 = 0. Both the MN and OF estimates of the normal component of flow are based on the same motion constraint equation (i.e., equation (3.31~. The first error equation used in the OF method (equation (3.32~) is a measure of the smoothness of the flow field. The second error equation (3.33) is the error in the estimation of the motion constraint equation. The objective function minimized in the OF method is the integral of the sum of these

55 errors- (equation (3.34~) where the parameter a!2 weights the smoothness error. If O!2 iS set to zero then the OF method essentially minimizes the error in the estimation of the motion constraint (3.31~. When Cl2 = 0 the solution for (3.34) using the OF method converges to the solution of (3.31) using the MN method. The Tangential Component of Velocity General Considerations. In the previous section, only the component of flow normal to isobrightness contours was discussed because the tangential component of flow cannot be calculated directly. The difficulty arises from what is known as the aperture problem and manifests itself in different ways for each motion estimation algorithm. In the MU case, the problem occurs when the motion of an oriented elide is detected by a direction-sensitiv unit that is small compared to the moving edge. Then the only information that can be extracted is the component of motion perpendicular to the lo- cal orientation of the edge. Hence the component of motion oriented along the edge is invisible. In both the OF and MN methods for estimating the normal component of velocity, the aperture problem manifests itself as the annihilation of the tangential component in the basic constraint equation (3.31~. An alternative interpretation of the aperture problem occurs if a point of brightness along an isobrightness contour moves along that contour from time to to time t2; this motion cannot be cletected. These consiclera- tions show that a direct method for computing the tangential component of flow is not possible. Vector decomposition of the known total velocity field, given a known normal component of flow, however, can yield an estimate of the tangential component of flow. These considerations are consistent with a proof given by Verri and Poggio (1989~. An Indirect Solution. Given the total flow field computed on a rect- angular grid from the pattern matching method mentioned previously, one can take the normal component of flow and perform a vector subtraction to obtain the tangential component. This decomposition was performed using the three normal component representations discussed above. The MU nor- mal component of velocity was chosen for this purpose because WahT and Simpson (199Ob) have shown that it produces better estimates of the normal component of velocities than either the OF or MN method. The MU method velocities are computed only at the points of zero cross- ings of edges. These vectors must be spatially aligned on the same grid as

56 that of the total velocity prior to the decomposition. To estimate the normal components of velocity at the locations of the tote] flow, the MU normal components were subsectioned into the same size tiles as the total flow. The normal velocities which fell within the region of the pattern tile area of the first image were then averaged to produce an overall average normal veloc- ity for a given tile. This procedure is consistent with the assumption that the total velocity vector produced by the pattern matching technique repre- sents the average velocity of the feature in a pattern tile. The mean normal component of flow was then decomposed into its Cartesian components ant! these components were used in the final vector decomposition to compute the tangential component of flow. Example A sequence of cloud-free AVHRR images for a region off the central CaTi- fornia coast was co-reg~stered to a latitude-Iongitude grid and calibrated to SST (Plate 3.3~. The image is stored as a matrix where, by tradition in image analysis, the row index is referred to as a line and the column index is referred to as a sample. Co-registration is the process of mapping images observed in the line and sample domain to the same latitude and longitude domain using an appropriate map projection (e.g., Brush, 1985; Jezclhing and Jain, 1989~. Calibration is the process of converting the raw measured brightness counts in one or more of the images to a geophysical variable (e.g., Kaufman, 1988~. This sequence is characterized by a cold-water filament ex- tending southward from the top of the sampled region. Therm structure edge maps for time step 2 of the image sequence were computed using the LOG operator with a value of it = 5 (Plate 3.4a). Motion inferred from these edge maps agrees wed with estimates of the total flow field computed using the pattern matching method (Plate 3.4b). These edge maps were then used to compute the normal component of velocity of the thermal structure over time using (3.30) at the zero-crossing points. A centered finite difference scheme was used to compute the spatial gradient of I', and a single time-centered difference was used to approxi- mate the temporal derivative of ·'. It is important to emphasize that the MU method (i.e., equation (3.304) only yields an estimate of the normal component of velocity near a well-defined edge (Plate 3.4a). The normal component of flow so obtained (Plate 3.4c) accurately represents motion in- ferred from the edge maps. Note especially the north-south oriented feature in the center right region (see region marked 3 in Plate 3.4c). The tangen-

57 tie] component of flow (Plate 3.41) again shows good correspondence with motion inferred from the edge maps. IdeaDy the estimated normal and tangential components of flow should be orthogonal. The need for spatial averaging of the normal components, however, may introduce directional errors in the approximation of the tan- gential component. Wah} and Simpson (199Ob) have shown that typically the angles between the normal and tangential components of flow are be- tween 75° and 80°. Thus, this method generally will not produce an ex- act tangential solution. It does, however, produce an approximate tangen- tial solution, which can be useful in many oceanographic applications (e.g., computation of the offshore transport of nutrients associated with coastal upwelling). It should be reemphasized that there is no direct method for computing the tangential component of motion from sequences of image data. 3.2.3 Ice Floe Iclentification ant} Principal Curves Overview of Banfield and Raftery Algorithm Knowledge of the spatial distribution, size, and shape of ice floes is critical to understanding physical processes in polar regions and the potential role of these processes in studies of global warming. Moreover, in high-latitude zones, shipping, naval operations, fishing, and the successful deployment of offshore structures are all strongly influenced by the distribution of the polar ice pack. Banfield and Raftery (1989) have developed an automated procedure for identifying ice floes in Landsat data. Automated procedures are needed for several operational reasons: (~) to handle the huge volume of data; (2) to eliminate intercalibration problems associated with subjective analyses; and (3) to improve on the poor performance records of human analysts working under the adverse weather conditions often associated with polar operations. The Banfield and Raftery method uses principal curves (Hastie and Stuetzle, 1989), an erosion propagation algorithm, and a method for clustering about principal curves to automatically identify the floes. Only the major elements of the method are reviewed here: the interested reader is referred to Banfield and Raftery (1989) for details of the procedure. Hastie and Stuetzle (1989) developed the concept of a principal curve. A principal curve can be thought of as a smooth one-dimensional curve that passes through the middle of an m-dimensional data set. It is nonparametric, and its shape is suggested by the data; it thus provides a nonlinear summary

58 of the data (Banfield and Raftery, 1989~. A one-dimensional curve in m- space is an m-vector consisting of m functions of a single variable A, called coordinate functions. The variable ~ parameterizes the curve and provides an ordering along it; ~ often is the arc length along the curve. Let X ~ Rm be a continuous random vector. Then f(~) is a principal curve if E[x~f-~(x) = A] = fall ~ where f-~ (x) = max ~ ~ : fix—f(~ = inf fix—fall ~~ ~ Given the distribution of X, Hastie and Stuetzle (1989) proposed the follow- ing algorithm for finding f: Pi+ = E[x~f' i(X) = >] where fi is the ith iterate. If the distribution of X is unknown, then an esti- mate of f can be obtained from the ciata set Exit by estimating E[x~f' ~(X) —A]. Hastie and Stuetzle (1989) obtain this estimate by scatterplot smooth- ~ng. Banfield and Raft ery (1989) noted that scatterplot smoothers generally produce curves that are biased toward the center of curvature. They modi- fied the Hastie and Stuetzle (1989) principal curve estimation algorithm to reduce bias and variance by using projections of the data rather than the data itself to mode} the principal curves when the distribution is unknown. Next, Banfield and Raftery (1989) used an erosion-propagation (EP) algorithm to select potential edge Pixels and to provide an initial grouping of the edge pixels into floe outlines. The EP algorithm operates on a binary image. Hence, the l.andsat data must be binarized by thresholding prior to EP analysis. Banfield and Raftery (1989) justified this procedure by noting that the marginal distribution of pixel intensities in the high-resolution polar Landsat data is highly bimodal. They further noted that the final result is relatively insensitive to the precise value of the threshold chosen. The erosion part of the EP algorithm identifies potential edge elements by using standard concepts from mathematical morphology (Serra, 1982), while the propagation part keeps track of the floe to which an etlge pixel belongs by locally propagating the information about the e(lge elements into the interior of the floe as it is eroded. The algorithm is applied iteratively to the binarized image. -

59 Banfield and Raftery (1989) noted that the EP algorithm tends to sub- divide floes. Therefore, they developed a method, based on clustering about the closed principal curves, for determining which of the floes identified by the EP algorithm should be merged. This final component of the over- all procedure to identify ice floes in polar Landsat data is hierarchical and agglomerative. Example . Shown in Plate 3.5a is a polar Landsat image. This image is 200 x 200 pixels, where each pixel is an SO-m square. The entire image represents a 15 x 15 km area. Ice floes appear as the gray features against the darker background. Ice floe outlines for this image, obtained using the Banfield and Raftery (1989) algorithm, are shown in Plate 3.5b. The algorithm accurately identifies the distribution, size, and shape of the floes on space scales of 200-300 m and larger. 3.3 Storage and Image Representation A typical AVHRR image consists of five channels of matrix data. The matrix size typically is in the range of 4 000 lines by 2qOOO samples. Generally the v ~ ~ ~ A ~ · . ~ . . ~ ~ ~ . . ~ . . . ~ . ~ ~ ~ . . ~ ~ lU-olt sample IS stored as a lu-olt integer Wltn the upper ~ bItS 01 each sample filled with zeros. All of these factors work out to about 80 to 100 megabytes (Mb) of archived data for a typical single scene. Use of satellite data sets in the analysis of problems related to gIobal-scaTe climate processes may require the analysis of literaBy thousands of such images. Thus, there is a need to have efficient storage and economical data structures for proper and efficient representation of the image data. 3.3.1 Storage Considerations The primary archive of satellite data is generally the raw digitized teleme- try stream directly received from the satellite. For multispectral images, the data are usually band interleaved rather than band sequential, usually con- tain embedded calibration information, and often include other data needed for proper Earth location of the scene. This data structure is inherently one- dimensional and has little resemblance to the two-dimensional image data structures normally associated with satellite images. (Note, however, that

60 the more familiar two-dimensional satellite images are subsequently con- structed from this telemetry stream using some set of mathematical trans- formations.) These data generally are called "level 1" data and are mandated for primary archives because higher-level data (e.g., the two-dimensional im- ages) often are produced by irreversible transformations. Moreover, level 1 data usually consist of 10-bit (e.g., AVHRR, CZCS) or less (e.g., DMSP) data strings packed into zero-filled 16-bit integer format for convenience. Data compression techniques and optical disc storage technology clearly are required if the data storage issue is to be adequately addressed. (Data compression is the process of reducing the number of bits required to store a given amount of information without loss of information.) For example, tests with a recent compression algorithm using Lempel-Ziv-Welch (LZW) coding (Welch, 1984) conservatively show that the average SO-Mb scene can be compressed to 32 Mb. This represents a 60~o reduction in size from the original data set. Some scenes may be compressed by as much as 75~o with EZW coding. The EZW compression algorithm has two main competitors currently in common use: Huffman coding and run length coding (R£C). The Hu~- man coding algorithm is not well suited to satellite data: preliminary tests indicate that only a 14-20% reduction in size is achievable (Jain, 1989~. Fur- thermore, it is much slower than the EZW algorithm. The REC algorithm was originally designed to vectorize bit maps. It is designed to work on strings of bits which are all Is or Os (Jain, 1989~. This makes it impractical for satellite data, which tends to vary too much (i.e., the strings of uniform Is or Os are too short). While the upper 6 bits of each 16-bit sample can be coded efficiently with the REC algorithm, the remaining 10 bits pose a se- rious problem for REC methods. Preliminary tests show that the estimated size reduction obtained from R£C algorithms wiB be in the 10-25~o range. The speed of the REC is comparable to that of EZW. These considerations show that the EZW compression algorithm is best suited to the proposed task. Preliminary tests also indicate that decompressing a typical satellite pass so that it can be used for analysis can take as much as 20 minutes per pass. This time factor depends on the speed of the disc and processing unit. No compression algorithm can get around this problem; increased access time is the trade-off for decreased space usage. Data compression techniques are also needed for the two-climensional higher level satellite images. Unfortunately, most two-dimensional algo- rithms either achieve speed by creating distortion in the data or achieve lack

61 of distortion by requiring excessively long execution times. Considerable re- search is needed to develop efficient two-dimensional compression/decompres- sion algorithms which do not distort the data. 3.3.2 Image Representation The data structure selected to represent the spatial data in an image wiD have a critical effect on the implementation and final performance of the analysis algorithm. The quadtree and octree are hierarchical data struc- tures often used to represent spatial data. The term quadiree is user! to describe a class of hierarchical data structures whose common property is that they are based on the principle of recursive decomposition of space (Samet, 1989~. They can be differentiated on the following bases: (1) the principle used to determine the decomposition process, (2) the resolution (variable or constant), and (3) the type of data they are used to represent. The prime motivation for the development of the quadtree is the need to reduce the amount of space necessary to store data through the use of ag- gregation of homogeneous blocks (Samet, 1989~. An important by-product of this aggregation is the reduction in execution time of an analysis process. Quadtrees have proved to be useful data structures for dithering algorithms, computing geometric properties of images, implementation of linear image transformations, development of hierarchical hidden-surface algorithms, and ray tracing. The quadtree is only one of several digital data structures useful in spatial statistics and digital image analysis. The constraints on and the need for the large-scale use of remotely sensed images in studies of global change is clear. Research is required in areas of both data storage and image representation if optimal use of remotely sensed data by the earth sciences community is to be achieved. 3.4 Special Considerations Remote sensing of the environment with earth-observing satellites poses some additional considerations beyond those normally encountered in labora- tory-based applications of digital image analysis. In the laboratory, both il- lumination and viewing geometry can be controlled. Moreover, the imaging detector is close to the object being cletectecI, and any interfering influence between imaging detector and object can be minimized. Finally, one gener- ally has a good notion of what constitutes the detected object. In contrast, earth observing satellites typically fly 800-900 km above the surface of the

62 earth. Viewing geometry and illumination are not controlled and can vary greatly from orbit to orbit. The SOO-900 km layer of atmosphere between target and detector acts as a filter that varies spatially and temporally, often partially corrupting image quality. there is the need to accurately project the image taken 900 km above the Earth onto a flat map projection suit- able to the particular application under study. Thus atmospheric correction algorithms (e.g., Curran and Dungan, 1989; Kaufman, 1988; Gratzki en cl GerstI, 1989; Simpson and Humphrey, 1990), sensor calibration algorithms (e.g., Gordon et al., 1983; Eckstein and Simpson, 199Oa), and earth location algorithms (e.g., Brush, 1985; Goshtasby et al., 1986; Jezching and Jain, 1989) are all pre-processing steps essential prior to meaningful mathemati- cal analysis. . . 3.5 Summary Remote sensing provides the only practical way to obtain the large-scale synoptic data sets necessary to address major problems of global signifi- cance in earth science that are fundamentally cross-disciplinary in nature. The magnitude of the data set, however, may preclude meaningful under- standing unless it can be distilled and synthesized into organized patterns of variance that can be meaningfully related to the underlying governing physics, chemistry, and biology of global-scale cycles and processes. Spa- tial statistics and mathematical methods of digital image analysis provicle mechanisms for such a synthesis. The examples cited herein included prin- cinal component analyses, which are useful for image segmentation and for determining spatial patterns of variance in large data sets; edge (letection; pattern matching; optical flow methods, which are useful for determining fields of motion from sequences of image data; and principal curves, which are useful for determining the spatial distribution, size, and shape of ice floes observed from spacecraft data. Atmospheric correction algorithms, sensor calibration algorithms, and earth location algorithms generally are required as pre-processes to digital image analysis of remotely-sensed images. Each of these pre-processing areas contains challenging mathematical problems which will have to be solved before earth sciences can benefit from the full potential of remote-sensing technology. O ~

63 Acknowledgments This work was sponsored in part by the Marine Life Research Group of the Scripps Institution of Oceanography and by a grant from the Office of Naval Research. AVHRR data (Plates 3.1, 3.3, and 3.4) were provided by the Scripps Satellite Oceanography Center. Landsat data (Plate 3.5) and permission to use material from the Banfield and Raftery ice floe study were generously provided by Professor Raftery of the University of Washington. Dr. Barbara Eckstein produced the sample EOF analysis (Plate 3.2) while working as an ONR-sponsored postdoctoral research assistant in the author's research group at Scripps. Several other present or former members of this group (I`. Al-Rawi, D. Atkinson, J. Bloom, T. Gallau(let, C. Humphrey, J. Toman, and D. WahI) contributed in one way or another to the completion of this chapter. S. McBride typed initial drafts of the manuscript, and F. Crowe and G. Tupper assisted with figure preparation. Bibliography [1] Aggarwal, J. K., and N. Nandhakuymar, On the computation of motion from sequences of images A review, Proceedings of ·EEE 76 (1988), 917-935. [2] Anderberg, M. R., Cluster Analysis for Applications, Academic Press, New York, 1973. t3] Banfield, J. D., and A. E. Raft ery, Ice Floe Identification in Satellite Images using Mathematical Morphology and Clustering about Principal Curves, Technical Report No. 172, Department of Statistics, University of Washington, Seattle, 1989.

64 [4] Barnett, T. P., and W. C. Patzert, Scales of temperature variability in the tropical Pacific, J. Phys. Oceanogre 10 (1980), 529-540. [5] Browning, K. A., Use of radar and satellite imagery for the mea- surement and short-term prediction of rainfall in the United King- dom, in Remote Sensing Applications in Meteorology and Climatology, R. A. Vaughn, ea., NATO AST Series, vol. 201, 1986. [6] Brush, R. J. H., A method for read time navigation of AVHRR imagery, IEEE Trans. Geosci. and Remote Sensing 23 (1985), 876-~87. [7] Campbell, J. D., Edge Structure and the Representation of Pictures, Ph.D. dissertation, Department of Electrical Engineering, University of Missouri, Columbia, 1969. [~] Cheevasuvit, F., H. Maitre, and D. Vidal-Madjar, A robust method for picture segmentation based on a split and merge procedure, Comput. Vision, Graph. Image Proc. 34 (1986), 268-281. [9] Curran, P. J., and J. A. Dungan, Estimation of signal-to-noise: A new procedure applied toAVIRIS data, IEEE Trans. on Geosci. and Remote Sensing 27 (1989), 620-627. [10] Davis, R. E., Predictability of sea surface temperature and sea level pressure anomalies over the North Pacific Ocean, '7. Phys. Oceanogr. 6 (1976), 249-266. [11] Davis, R. E., Techniques for statistical analysis and prediction of geo- physical fluid systems, Geophys. Astrophys. Fluid Dyn. 8 (1977), 245- 277. [12] Deschamps, P. Y., M. Herman, and D. Tanre, Modeling of the atmo- spheric effects and its application to remote sensing of ocean color, Appl. Opt. 22 (1983), 3751-3778. [13] Dudgeon, D. E., and R. M. Mersereau, Multidlimensional Digital Signal Processing, Prentice HaD, Englewood Cliffs, Hi., 1984. [14] Eckstein, B. A., and J. J. Simpson, Aerosol and RayTeigh radiance con- tributions to Coastal Zone Color Zone Scanner images, Int. ]. Remote Sensing, in press (199Oa). [15] Eckstein, B. A., and J. J. Simpson, Cloud screening Coastal Zone Color tone Scanner Images using channel 5, Int. ]. Remote Sensing, in press (199Ob).

65 [16] EOS, Science and Mission Requirements Working Group Report, Vol- ume ~ and Appendix, National Aeronautics and Space Administration' Goddard Space Flight Center, Greenbelt, Md., 1984. [17] Esbensen, K., and P. Geladi, Strategy of multivariate image analysis (MIA), Chemometrics and Intelligent Laboratory Systems 7 (1989), 67- 86. [18] Fukada, Y., Spatial clustering procedures for region analysis, Pattern Recognition 12 (1980), 395-403. [19] GalIaudet, T. C., and J. J. Simpson, Automated cloud screening of AVHRR imagery using split-and-merge clustering, Remote Sensing En- viron., submitted (1991a). [20] Galiaudet, T. C., and J. J. Simpson, Seasonal and non-seasonal vari- ability in sea surface temperature off Punta Eugenia, to be submitted to Remote Sensing Environ. ~ 199lb). [21] Geladi, P., H. Isaksson, L. Lin~quist, S. Vold, and K. Esbensen, Prin- cipal component analysis of multivariate images, Chemometrics and Intelligent Laboratory Systems 5 (1989), 209-220. [22] Gordon, H. R., J. W. Brown, O. B. Brown, R. H. Evans, and D. K. Clark, Nimbus 7 CZCS: Reduction of its radiometric sensitiv- ity with time, Appl. Opt. 22 (1983), 3929-3931. [23] Goshtasby, A., G. C. Stockman, and C. V. Page, A regional-based approach to digital image registration with subpixe! accuracy, IEEE Trans. Geosci. and Remote Sensing 24 (1986), 390-399. t24] Gratzki, A., and S. A. W. GerstI, Sensitivity of an atmospheric correc- tion algorithm for non-Lambertian vegetation surfaces to atmospheric parameters, IEEE Trans. Geosci. and Remote Sensing 27 (1989), 326- 331. [25] Hastie, T., and W. X. Stuetzle, Principal curves, A. Am. Stat. Assoc. 84 (1989), 502-516. t26] Henderson-Sellers, A., Defogging cloud determination algorithms, Na- t~re 298 (1982), 419-420. [27] Hildreth, E. C., The Measure of Visual Motion, MTT Press, Cambridge, Mass., 1983.

66 [28] Horn, B. K. P., and B. G. Schunck, Determining optical flow, Artif. Intell. 17 (1981), 185-203. [29] Ingerbritsen, S. E., and J. P. Lyon, Principal component analysis of multi-temporal image pairs, Int. J. Remote Sensing 6 (1985), 687-696. [30] Jain, A. K., Fundamentals of Digital Image Processing, Prentice Had, Englewood Cliffs, N.~., 1989. [31] Jain, A. K., and R. C. Dubes, Algorithms for Clustering Data, Prentice Hall, En~ewood Cliffs, N.~., 1988. [32] Jain, J. R., and A. K. Jain, Displacement measurement andits appli- cation to interframe image coding, IEEE Trans. Common. 29 (1981), 1789-1808. [33] Jensen, J. R., Introductory Digital Image Processing, Prentice Hall, Englewood Cliffs, N.J., 1986. [34] Jezching, T., and A. K. Jain, Registering I'andsat images by point matching, IEEE Trans. Geosci. and Remote Sensing 27 (1989), 642- 651. [35] Kaufman, Y. J., Atmospheric effects on spectral signature Measurements and corrections, IEEE Trans. Geosci. and Remote Sens- ing 26 (1988), 441-449. [36] Kelly, K. A., The influence of winds and topography on surface tem- perature patterns over the northern California slope, if. Geophys. Res. 90 (1985), 11,783-11,793. [37] Lagerioef, G. S. E., EOF analysis of AVHRR and CZCS data, Third Conference on Satellite Meteorology and Oceanography, 1986. [38] I.illesand, T. M., and R. W. Kiefer, Remote Sensing and Image Inter- pretation, 2nd ea., John Wiley and Sons, New York, 1987. [39] lLuenberger, D. G. Optimization by Vector Space Methods, John Wiley and Sons, New York, 1969. Marr, D., and E. C. Hildreth, Theory of edge detection, Proc. R. Soc. London, B 207 (1980), 187-217. [41] Marr, D., and S. Uliman, Directional sensitivity and its use in early visual primitives, Proc. R. Soc. London, B 211 (1981), 151-180.

67 [42] Mather, P. M., Computer Processing of Remotely-Sensed Images: An Introduction, John Wiley and Sons, New York, 1987. [43] Matveev, L. T., Computer Processing of Remotely Sensed Images: An Introduction, John Wiley and Sons, New York, 1984. [44] MuerIe, J. L., and D. C., Allen, Experimental evaluation of techniques for automatic segmentation of objects in a complex scene, pp. 3-13 in Pictorial Pattern Recognition, G. Cheng et al., eds., Thompson, Wash- ington, D.C., 1968. [45] Ohring, G., and P. F. CIapp, The effects of changes in cloud amount on the net radiation at the top of the atmosphere, A. Atmos. Sci. 37 (1980), 447-454. [46] Pailleux, F. S., The impact of satellite data on global numerical weather prediction, in Remote Sensing Applications in Meteorology and Clima- tology, R. A. Vaughn, ea., NATO AST Series, vol. 201, 1986. [47] Barman, D., and J. Kittler, Clustering algorithms for use with images of clouds, Int. ~7. Remote Sensing 7 (1986), 855-866. [48] Pavlidis, T., Structural Pattern Recognition, Springer-VerIag, Berlin, 1977. [49] Preisendorfer, R. W., Principal Component Analysis in Meteorology and Oceanography, C. D. Mobley, compiler ea., Developments in Atmo- spheric Science 17, Elsevier, Amsterdam, 1988. [50] Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes, Cambridge University Press, New York, 1986. [51] Ramanathan, V., The role of earth radiation budget studies in climate and global circulation research, '7. Geophys. Res. 92 (1987), 4075-4095. [52] Richards, J. A., Remote Sensing Digital Image Analysis: An Introduc- tion, Springer-VerIag, New York, 1986. [53] ~ -- gamut, tit' Applications of Spatial Data Structures, Computer Graphics, Image Processing and GIS, Addison-Wesley Publishing Corporation, New York, 1989. t54] Seddon, A. M., and G. E. Hunt, Segmentation of clouds using cluster analysis, Int. '7. Remote Sensing 6 (1985), 717-731.

68 [55] Serra, J., Image Analysis and Mathematical Morphology, Academic Press, New York, 1982. [56] Simpson, J. J., Image segmentation using recursive polygon fill opera- tions, Remote Sensing Environ., submitted (1991~. [57] Simpson, J. J., and J. Bloom, Objective estimates of velocity from sequences of satellite data using log search and minimum distortion methods, if. Geophys. Res., submitted (1990~. [58] Simpson, J. J., and C. Humphrey, An automated cloud screening al- gorithm for daytime AVHRR imagery, J. Geophys. Res. 95 (1990), 13,459-13,481. [59] Singh, A., and A. Harrison, Standardized Principal Components, Int. ]. Remote Sensing 6 (1985), 883-896. [60] Strub, P. T., C. James, A. C. Thomas, and M. A. Abbott, Seasonal and non-seasonal variability of sateDite-derived surface pigment concentra- tion in the California Current, '7. Geophys. Res. 95 (1990), 11,501- 11,530. [61] Verri, A., and T. Poggio, Motion field and optical flow: Qualitative properties, IEEE Trans. Pattern Anal. Machine Intell. 11 (1989), 490- 498. [62] WahI, D. D., and J. J. Simpson, Physical processes affecting the objec- tive determination of near-surface velocity from satellite data, '7. Geo- phys. Res. 95 (199Oa), 13,511-13,528. [63] WahI, D. D., and J. J. Simpson, SateDite-derived estimates of the nor- mal and tangential components of near-surface flow, Int. ]. Remote Sensing, in press (199Ob). [64] Welch, Terry A., A technique for high performance data compression, IEEE Computer 17 (June 1984), 8-19.

69 Appendix to Chapter 3 Listed below are definitions and explanations for abbreviations commonly used in remote sensing. AVHRR The Advanced Very High Resolution Radiometer on the Na- tional Oceanic and Atmospheric Administration's polar-orbiting weather satellite. It measures cloud cover and infrared sea surface temperature. DFT—Discrete Fourier transform. EOF Empirical orthogonal function. EOS Earth Observing System. A proposed National Aeronautics and Space Administration program for earth observing systems to be launched between 1997 and 2007. FFT—Fast Fourier transform. GOES Geostationary Operational Environmental Satellite. tional weather satellite used to measure cloud cover. solar radiation often are computed from GOES data. LOG Laplacian of the Gaussian operator. An opera- Estimates of EZW Lempel-Ziv-Welch coding used in data compression algorithms. MN Minimum norm solution. MU—The Marr-UlIman solution for the normal component of velocity. OF Optical flow method of computation. REC Run length coding used in data compression algorithms. SST Sea surface temperature.

Next: 4. Spatial Statistics in Environmental Science »
Spatial Statistics and Digital Image Analysis Get This Book
×
Buy Paperback | $65.00
MyNAP members save 10% online.
Login or Register to save!
Download Free PDF

Spatial statistics is one of the most rapidly growing areas of statistics, rife with fascinating research opportunities. Yet many statisticians are unaware of those opportunities, and most students in the United States are never exposed to any course work in spatial statistics. Written to be accessible to the nonspecialist, this volume surveys the applications of spatial statistics to a wide range of areas, including image analysis, geosciences, physical chemistry, and ecology.

The book describes the contributions of the mathematical sciences, summarizes the current state of knowledge, and identifies directions for research.

  1. ×

    Welcome to OpenBook!

    You're looking at OpenBook, NAP.edu's online reading room since 1999. Based on feedback from you, our users, we've made some improvements that make it easier than ever to read thousands of publications on our website.

    Do you want to take a quick tour of the OpenBook's features?

    No Thanks Take a Tour »
  2. ×

    Show this book's table of contents, where you can jump to any chapter by name.

    « Back Next »
  3. ×

    ...or use these buttons to go back to the previous chapter or skip to the next one.

    « Back Next »
  4. ×

    Jump up to the previous page or down to the next one. Also, you can type in a page number and press Enter to go directly to that page in the book.

    « Back Next »
  5. ×

    To search the entire text of this book, type in your search term here and press Enter.

    « Back Next »
  6. ×

    Share a link to this book page on your preferred social network or via email.

    « Back Next »
  7. ×

    View our suggested citation for this chapter.

    « Back Next »
  8. ×

    Ready to take your reading offline? Click here to buy this book in print or download it as a free PDF, if available.

    « Back Next »
Stay Connected!