SPATIAL AND SPECTRAL ANALYSES OF REMOTELY SENSED IMAGES USING SCALE-SPACE TECHNIQUES Brian Zuerndorfer 1991 RL-939 = RL-939

To my wife Celia and my daughters Carol and Sarah ii

ACKNOWLEDGEMENT This dissertation is the culmination of my formal education, and involved the help of many people. I am grateful to my committee for their suggestions and encouragement. I am indebted to my chairmen, Professors Wakefield and England, for their guidance and understanding. I consider myself fortunate to have these gentlemen chair my committee, and I am proud to call each my friend and colleague. I thank Professor Fawwaz Ulaby and Craig Dobson for their aid and guidance throughout the course of my work. I thank Leland Pierce for the use of his software and expertise, which have improved the quality of my work. I thank fellow and former students Richard Austin, Joe Landry, John Galantowizc, Andy DeJaco, Huili Wang, Brent Edwards, John Feng, Chih-Ming Chiang, Gail Feinman, and K.R. "Rags" Raghavan. Their friendship and support have eased an arduous process. I also thank Patty Humphrey, who cheerfully helped produce manuscripts, and their numerous revisions, over the past four years. Finally, none of this would have been possible without the constant love and patience of my family. To my parents, who are an unceasing source of encouragement. And most of all, to my wife Celia and daughters Carol and Sarah, who have sacrificed so much and are my greatest source of happiness. iii

TABLE OF CONTENTS DEDICATION................................................................................................... ii ACKNOWLEDGEMENTS....................................... iii LIST OF FIGURES....................................................... vi LIST OF TABLES............................................. ix LIST OF APPENDICES....................................................... x CHAPTER 1. REMOTE SENSING SYSTEMS.................................................... 1 1.1 Passive Microwave Measurement Systems................................ 2 1.2 Data Fusion and Multiresolution Systems...................................... 10 1.3 Thesis Overview.................................................... 14 2. SCALE-SPACE: REVIEW AND EXTENSIONS............................. 16 2.1 Scale-Space Filtering....................................................... 16 2.2 The Function-Crossing Problem..................................................... 22 2.3 Non-Gaussian Filtering...................................................... 38 3. MULTIRESOLUTION PROCESSING............................................. 54 3.1 Background........................................... 54 3.2 Boundary Migration...................................................... 56 3.3 Boundary Extrapolation....................................................... 74 4. FREEZE/THAW CLASSIFICATION............................................... 90 4.1 Passive Microwave Radiation...................................................... 90 4.2 Emissions From Frozen Soil.................................................... 100 4.3 Threshold Determination.................................. 112 5. FREEZE/THAW BOUNDARY PROCESSING.................... 119 5.1 Boundary Localization Procedure.................................................... 119 5.2 SMMR images.................................................... 123 5.3 Automated Boundary Localization................................................... 131 6. SUMMARY AND FUTURE WORK................................................... 147 iv

6.1 Summary.............................................................. 147 6.2 Future Work....................................................... 149 APPENDICES................................................................................................... 151 BIBLIOGRAPHY............................................................................................. 176 v

LIST OF FIGURES Figure 1.1 Radiometer-centered coordinate system........................4 1.2 Radiometer beampattern............................. 6 1.3 Multiresolution system....................................13 2.1 Scale-space fingerprints having embedded structure for the 1-D case............................................ 20 2.2 Input signal, and associated fingerprints, for frozen terain, warm terrain at a constant temperature, and constant threshold........30 2.3 Input signal, and associated fingerprints, for frozen terain, warm terrain at a variable temperature, and constant threshold........31 2.4 Input signal, and associated fingerprints, for frozen terain, warm terrain at a variable temperature, and piecewise linear threshold................................................................................................. 32 2.5 Family of fixed-support filters...............................43 2.6 Family of relative-support filters.............................46 3.1 Scale-space system diagram................................ 55 3.2 Migration coefficient value as a function of normalized edge.........65 3.3 Ideal signal and threshold..................................66 3.4 Even and odd components for right edge of ideal signal............. 67 3.5 Fingerprint and migration bound for right edge of ideal signal. 67 3.6 Even and odd components for left edge of ideal signal..............68 3.7 Fingerprint and migration bound for right and left edges of ideal signal............................................ 68 3.8 Ideal signal for non-symmertic thresholding......................................... 72 vi

3.9 Ideal single-senor and indicator signals................................................. 75 3.10 Fingerprints and scale limits for ideal single-sensor and indicator signals........................................ 81 4.1 Specific intensity.................................................................................... 91 4.2 Radiative transfer m odel........................................................................ 92 4.3 A pparent tem perature............................................................................. 99 4.4 Test area with air and soil temperature stations..................................... 101 4.5 TB(37) versus measured air temperature at meteorological sites in North Dakota and the surrounding region.......................................... 107 4.6 3TB/3f versus measured air temperature at meteorological sites in North Dakota and the surrounding region.......................................... 108 4.7 aTa/3f versus TB(37) at meteorological sites in North Dakota and the surrounding region.............................................. 109 4.8 Scatter diagram of aTB/3f versus TB(37) throughout North Dakota and the surrounding region............................................... 113 4.9 Single class ellipses of aTB/af versus TB(37) throughout North Dakota and the surrounding region.............................................. 116 5.1 Normalized Seasat SMMR and Gaussian beampatterns................ 122 5.2 A comparison of rerported air and soil temperatures with images of North Dakota and the surrounding region; unrefined criteria, O ctober 24, m idnight.............................................................................. 127 5.3 A comparison of rerported air and soil temperatures with images of North Dakota and the surrounding region; unrefined criteria, D ecem ber 11, noon................................................................................ 128 5.4 A comparison of rerported air and soil temperatures with images of North Dakota and the surrounding region; refined criteria, October 24, midnight..................................... 129 5.5 A comparison of rerported air and soil temperatures with images of North Dakota and the surrounding region; refined criteria, D ecem ber 11, noon............................................................................. 130 5.6 Automated boundary migration and classification for images of North Dakota and the surrounding region; October 24, midnight......... 134 vii

5.7 Automated images of North Dakota and the surrounding region; October 24, midnight........................................................ 138 5.8 Automated images of North Dakota and the surrounding region; October 30, midnight....................................................... 139 5.9 Automated images of North Dakota and the surrounding region; November 1, noon........................................................ 140 5.10 Automated images of North Dakota and the surrounding region; November 5, midnight....................................................... 141 5.11 Automated images of North Dakota and the surrounding region; November 27, midnight.................................................. 142 5.12 Automated images of North Dakota and the surrounding region; November 29, noon................................................. 143 5.13 Automated images of North Dakota and the surrounding region; December 3, midnight................................................. 144 5.14 Automated images of North Dakota and the surrounding region; December 9, midnight................................................. 145 5.15 Automated images of North Dakota and the surrounding region; December 1 1, noon................................................. 146 viii

LIST OF TABLES Table 4.1 SMMR performance characteristics..................................................... 100 4.2 Contributions of thermal gradient to radiobrightness and spectral gradients for noon and midnight data..................................................... 104 4.3 Cluster centrod in decision space........................................................... 115 4.4 Freeze/thaw criteria in decision space.................................................... 115 5.1 Filter bandwidths for resolution compensation................................... 121 5.2 Time summary for images of Figure 5.7 through Figure 5.15............ 135 5.3 Air temperatures for October through December measurement dates........................................................................................................... 136 5.4 Soil temperatures for October through December measurement dates........................................................................................................ 137 ix

LIST OF APPENDICES Appendix A Conditions for Embedded Structure in Function-Crossing Fingerprints............................................................................................ 152 B Proof of Theorem 3-1 (1-D)...................................... 161 C Proof of Theorem 3-2 (2-D).................................................. 165 D Proof of Lemma 3-3......................................... 172 x

CHAPTER 1 REMOTE SENSING SYSTEMS This thesis presents an image processing technique applied to remote sensing data. Remote sensing commonly refers to the remote monitoring of the atmospheres and surfaces of planetary objects, with special emphasis on earth observations [74]. A fundamental use of remote sensing is in the classification and extraction of information of land surfaces. In classification, images from remotely sensed land surfaces are segmented into regions, and the surface type of each region is identified. For example, surfaces can be classified as cultivated, forested, pasture, urban, or water-covered using synthetic aperture data [1 1], and as desert, jungle, wet, snow-covered, and ice-covered using passive microwave data [3,38]. That is, classification is a detection process where algorithms detect changes in surface types, and locate boundaries between surfaces. In the extraction of information, the characteristics of surface regions are inferred from remotely sensed data. For example, remotely sensed data can be used to estimate snow depth [16,46], the extent of surface vegetation [19,33], and soil moisture content [42, 57,64]. That is, the extraction of surface information is an estimation process where surface information is inferred from previously classified surfaces. This thesis is concerned with surface classification on the basis of the state of soil moisture (i.e., frozen ya thawed terrain). Soil moisture contributes to the energy exchange between the air and ground through latent heats of fusion and vaporization. 1

2 Whether as boundary conditions for mesoscale climate modelling, or as inputs to an agricultural productivity model, the amount and state of soil moisture are regional parameters that one would like to determine through satellite remote sensing. There is a large body of literature that addresses the estimation of soil moisture from passive remotely sensed measurements [9,12,14,38,42,57,62-64,79]. In this thesis, we infer moisture state from passive satellite measurements, and classify frozen terrain in the northern Great Plains. Data used are obtained from the Nimbus-7 scanning multichannel microwave radiometer satellite. 1.1 Passive Microwave Measurement Systems The remotely sensed data considered in this thesis are thermally emitted, microwave energy of a surface. The systems used to measure these emissions are called radiometers. The amount of energy emitted by a surface at a given frequency can be represented as a brightness temperature. The microwave brightness temperature, or radiobrightness, of a surface has the units of Kelvin and is defined as the physical temperature of an equally bright blackbody radiator. An ideal surface (i.e., a perfect radiator) would have a radiobrightness equal to its physical temperature. The ratio of radiobrightness to physical temperature for a real surface ranges between zero and one, and is determined by the electrical and geometrical properties of the surface. A more detailed treatment of radiobrightness is given in chapter 4. Radiometers cannot measure radiobrightness directly, but measure antenna temperature instead. Antenna temperature represents the total received energy at the output of the radiometer antenna, and includes the effects of atmospheric emissions, atmospheric attenuation, surface reflections, radiometer antenna beamshape and ohmic losses, as well as surface emissions. The antenna temperature of a radiometer pointed in the direction (00, O,) can be expressed as [74],

3 TA(9,,0, X) =- dO f do sin OTe,(0, 4;X)F,(0, 4;e,,( X) (1.1) where X is the receive wavelength of the radiometer, A, is the (constant) receiving area of the radiometer antenna, Tp(.,.) is apparent temperature, and Fn(.,.) is the normalized radiation pattern, or antenna pattern. The spatial variables 0 and ( are radiometer centered, spherical coordinates (Figure 1.1). The apparent temperature, TAp,( 0,), represents the total received energy at the radiometer antenna aperture, at wavelength X, along the direction (0, ~). The apparent temperature includes terms that account for atmospheric effects, surface reflections, and surface brightness. The radiometer antenna pattern is normalized with a peak value F,(O, 0) = 1. A more detailed discussion of apparent temperature is also given in chapter 4. At the relatively low frequencies considered in this thesis, atmospheric effects and surface reflections are ignored1. As a result, apparent temperature is approximately the brightness temperature, TAP(Q, O;X) TB,(A;X)9 and (1.1) can be written as, A 2 TAM, X4)= =2 dO do sin OTB(,;X)F^(;,,) (1.2) 1 See [74] for a more detailed treatment of atmospheric effects.

4 Figure 1.1. Radiometer-centered, spherical coordinate system. For narrow-beam antennas, (1.2) can be approximated by a convolution integral. Consider a radiometer antenna with a beampattern F^(a, f3;X), where angles a and 3 are taken with respect to beam-center (Figure 1.2). By approximating narrow-beam antennas as having limited spatial support, (1.2) is written as, TA (A, o, X) = A +D, *o+Dp =X2 do d sin OTB(O,(;X)F,.(O -,( -,0) sin;X). (1.3) o -D, o -DP where the limits of integration, Dt and Dp, are determined from the beampattem F^(a, P;X), and 00 >> D > 0. We further assume 0o +D, < n/2 so that TB(O, 4;.) is determined solely by surface emissions for all (0, 4). We define offset angles t and p by,

5 0e=00o+ and, 0 >>D, o00 >,' l e6 [-D,,+DJ. Thus, sin(Oo +t ) = sin 00 + tcos00 = sin 00, and (1.3) can be re-written as, TA(00, 00, X) = +D, +D sin OoA I f D = -J d f dpTB(0o +,, 0 + p;oX)F(no(t, p;X), -D, -Do (1.4) for, F.(t, p;X) - F,,(t, p sin 0,0;).

6 x Y I z Figure 1.2. Radiometer beampattern. In (1.4), TA is approximated as a convolution in spherical coordinates. To express TA as a convolution in carteasian coordinates for a satellite at altitude H above a planar surface, a variable r can be given by (Figure 1.1), r = H tan0, so that for small T, H tan(,o +t ) - ro + a where, r, = H tan 0, ca H tan t.

7 As a result, (1.4) can be written as, TA(ro 50, )) = sinOA,+Dr _ do +Dr X2 if H(l +(ro/H)2) dpTB(ro + a, 4o + p;X)Fo(a, p;X) +,, +D +D - 2 Jdo dPT (ro + <oto + p;))Po(o, p;X), (1.5) -Dr -Dp where, -1 I ro + Oh TB(rO +o(, o + p;X)= TB(an-f -L H ) +.P,(a, p;X) Fftan - pa D, H tanD,. The radiometer data considered in this thesis is obtained from a satellite that can be modelled as forward-looking with a fixed incidence angle (i.e., 00 and r, are constant), narrow field-of-view, and platform motion in the forward (i.e., x) direction. As a result, we have, x r Y roO, so that (1.5) can be approximated as,

8 +D +D I C(OS2'0d sin 0TArB( d8, YO +, -A-(-O )r odys -ir TBr ( (1.6) H i o r o x Xo = ro D =D,, and, Yo =rolo Dy = rODp. In the convolution integral of (1.6), variations in y0 are produced by the scanning of the radiometer antenna in the ( direction, and variations in x0 are produced by satellite motion in the x direction. Two-dimensional antenna patterns are usually modelled as circularly symmetric or separable into products of orthogonal, one-dimensional patterns. As a result, antenna patterns are either given as F.((Kc sin 0)/x), for circular symmetric patterns, or by F^((Kra sin a)/X), (Krb sin P)/X)), for separable patterns, where Kc, K., and Kr are constants and (Figure 1.2), tan a = cos ) tan 0 tan 3 = sin ( tan 0. In the case of narrow beam antennas,

9 sin a = a sin P = P sin 0= 0 = (a2 + P2)l/2 and the antenna patterns are given by, F(K, sin ) = (K)(+ ) or, F (K,, sin a K^sin P)f ' (Ka Kb By substituting either of these antenna patterns into the narrow-beam antenna temperature expression (1.3), the filter kernel in (1.6) becomes, F^o h;X) y(h5 ) X2 X2 XX for some function T(.,.). As a result, (1.6) can be written as, +Dx +D -pf2 TA(xo, y, X) y d6T(x + y, yo + 6;5) (1.7) -D, -D7 where T(.,.;X) is a function derived from the radiobrightness Tf(.,.;), and Cis a constant. In the case where the radiobrightness is independent of X (i.e., T(.,.;X) = T(.,.)), the filter

10 output is said to be a linearly-scaled with respect to wavelength. That is, X is linearly related to the equivalent low-pass spatial bandwidth of the antenna. Equivalently, there is an inverse relationship between the antenna resolution and wavelength. The importance of such a linear scaling property is seen in the scale-space filtering developments of chapter 2. 1.2 Data Fusion and Multiresolution Systems The spectral characteristics of surfaces are fundamental parameters in the classification of surfaces using radiobrightness. For the fixed altitude, constant incidence angle radiometer considered in this thesis, the radiobrightness of a surface is a function of frequency and surface composition. That is, a surface of a particular composition will have a unique spectral response, and sensors that measure surfaces at multiple frequencies are used to estimate the spectral characteristics of surfaces [35,40]. Region classification can be performed unambiguously using such multispectral processing if the spectral response of each region is unique for the combination of frequencies used. Multispectral systems are a subset of a class of systems referred to as multisensor or multisource systems. These systems are the subject of much current research in computer vision and remote sensing, under the heading of multisensor or multisource data fusion (e.g., see [43,61,66,80]). Multisensor systems, as used in remote sensing systems, geologic exploration, nuclear medicine, and robotics, process signals from a variety of sensors to classify regions within the sensed environment. As with multispectral systems, classification is unambiguous for multisensor systems if the response of each region is unique for the combination of sensors used. The performance of multisensor systems is limited by the response properties of the individual sensors, as well as by our ability to integrate information across different classes of sensors. For example, antennas, cameras, or acoustic arrays have point spread

11 functions (PSF's) that limit the resolution of the sensor systems that use them. The resolution of a sensor affects the accuracy with which regional boundary or "edge" information can be extracted. The relationship between sensor resolution and system performance is further complicated in multiresolution sensing systems, where systems integrate information across sensors with different PSF's. In multispectral systems, data from M sensors are integrated, where each sensor operates at a different (center) frequency [58]. In these systems, each sensor is scanned over a surface, and the output signal of each sensor is a function of spatially varying radiation intensity of the surface and the PSF of the sensor. Because surface radiation and sensor PSF's are functions of frequency (section 1.1), multispectral (classification) systems are also multiresolution systems. In present data fusion research, effects of PSF (resolution) differences between sensors on (multiresolution) system performance are addressed in two ways. One approach is to assume that all sensor PSF's are adequate for the job-at-hand, and are not the limiting factor in system performance [8] (i.e., PSF differences are neglected). A more common approach is to incorporate the effect of PSF differences indirectly through "uncertainty" functions, which differ from one sensor to another. These functions, assigned one per sensor, reflect the measurement uncertainty at each sensor, and have the forms of probability density functions [23,48], belief functions [69,71], or fuzzy membership functions [1]. In using these uncertainty functions, the measurements from all sensors are processed directly, and components common to all sensor measurements are estimated. These common components are then related to characteristics of the environment being sensed. The effect of PSF variation between sensors is to reduce the amount of information common to all sensor measurements. This thesis considers an alternative approach that first maximizes the information common to all sensor

12 measurements, by compensating (equalizing) the PSF's of all sensors prior to multisensor system processing, and then recovers information that may have been lost in compensating the PSF's. A model of a general multiresolution system is shown in Figure 1.3. In the figure, the input signal vector I(K(x)) from a sensed region is measured by M different sensors. Position within the sensed environment is denoted by the independent variable x, and K(x) is the vector signature (response) of the environment at x, kl(x) K(x)) gkM(x) The jh sensor input, ij(kj(x)), is a measurement of the j" component of I(K(x)), and is dependent on the jr component of K. The ideal (i.e., noiseless) output, rj(kj(x),sj), of the j sensor is dependent on kj(x) as well as the sensor resolution, Sj(x). Uncorrelated noise in the jth sensor output is shown as nj. The goal of the system is to estimate the signature, K(x), on the basis of measurements from the M sensors. In the case of multispectral processing, ij(kj(.)) is the surface radiobrightness at wavelength Xj, and rj(kj(.),sj) is the corresponding antenna temperature. In Figure 1.3, direct processing of the outputs rj(kj(x),sj) produce ambiguous (uncertain) signature estimates. That is, differences in kj(x) between sensors can generate an identical K(x) as that generated by differences in sj. To eliminate the ambiguity, and increase the common components in the sensor measurements, sensor outputs are resolution compensated (pre-processed) so that PSF differences between sensors are

13 equalized prior to signature processing. The inputs for signature processing are then, effectively, the outputs from sensors having identical PSF's, at some pre-determined resolution. Sensor Input Sensor Output n, i,(kl(x)) rl(kl(x), s) Processor Output A K 0 S S iM(km(X)) rM (kM(x), SM) Figure 1.3. Multiresolution system. One common approach practice in multispectral processing is to compensate sensor PSF's to the coarsest-resolution sensor, and perform spectral signature processing at coarse-resolution [32]. Subsequent to spectral signature processing, regional boundary estimates are obtained at coarse-resolution. Such boundary estimates are non-parametric (no modelling is needed), but the resulting resolution degradation can produce significant localization errors in boundary estimation. In a second approach, sensord ro resolution can be compensated to the PSF of the finest-resolution sensor, and regional boundaries estimated at fine-resolution. Such resolution compensation is performed through

14 model-based estimation techniques, such as linear programming [53] or ARMA estimation [47]. However, the estimation accuracy of these techniques is strongly dependent on modelling accuracy and on the global applicability of a single model. In remote sensing, globally applicable surface models do not exist. Alternatively, we note that when measurement differences between sensors are modeled exactly by differences in PSF (i.e., there are no signature differences between sensors), the multiresolution data-integration problem reduces to a problem for which scale-space filtering [76,77,84,85] can be used. This suggests that scale-space filtering can be adapted for signature processing in multiresolution systems. That is, by applying scale-space filtering to the resolution-degradation approach above, fine-resolution boundary information may be recovered. 1.3 Thesis Overview The thesis describes a scale-space approach for boundary estimation with multiresolution systems. This approach is applied to multispectral processing in remote sensing, and the localization of freeze/thaw boundaries in the northern Great Plains using SMMR data. In chapter 2, the standard results of scale-space filtering are first reviewed. Extensions are derived that allow the standard scale-space results to be applied to more general multiresolution systems. Some, but not all, of these results are used in the multispectral application that follows. In chapter 3, scale-space results are applied to boundary estimation in multiresolution systems. A technique for boundary estimation technique is proposed, and the migration of boundary estimates are shown as a function of resolution. It is shown that boundary estimates made at coarse resolution can have significant localization errors. A modification of the estimation scheme is presented which is derived from scale-space

15 filtering and can mitigate the localization errors. Limitations of this approach are also developed. Both the boundary estimation and migration problems are posed with respect to 1-D contours along a 2-D surface. That is, we develop l-D theory assuming that 2-D surfaces can be projected onto 1-D contours with minimal loss of information. More general 2-D estimation results are also developed. In the chapter 4, we develop the multispectral processing used for identifying frozen and thawed terrain in the northern Great Plains. We review aspects of passive microwave remote sensing, and develop a multispectral algorithm and decision criteria for surface classification. In chapter 5, boundary localization is applied to the multispectral classification data of chapter 4. The results of boundary localization are shown and discussed relative to local ground data measurements. Chapter 6 presents summaries as well as discussions of future work derived from results of the dissertation.

CHAPTER 2 SCALE-SPACE: REVIEW AND EXTENSION Scale-space filtering was introduced in the early 1980's as a technique for signal analysis over multiple scales [4,76,77,83,84]. The origins of scale-space filtering lie in the edge detection concerns of computer vision. Since its introduction, numerous papers on the theoretical application of scale-space filtering to problems in computer vision have been published (e.g., [7,20,41,45,49,60,75,78,81,85]). The potential success of a scale-space approach to these problems has been tempered by questions of implementation. For example, Lu and Jain [51] consider the use of Gaussian filters with finite spatial support. Schunck [65] and Lindeberg [49] consider effects of filtering scale-space processing of discrete, rather than continuous, images. In this chapter, we consider the application of scale-space processing to multiresolution systems that use function-crossings and non-Gaussian filtering. We develop our results for the case of 1-D signals, and extend them to the 2-D case. 2.1 Scale-Space Filtering 2.1.1 The Original Development As developed by Witkin [76,77] and Yuille and Poggio [84,85], scale-space theory describes the effects of filter scale, s, on functions e(x,s) of the form, e(x,s) L{r(x,s)}, 16

17 where L{. } is a linear, differential operator and r(x,s) is the output of a linear shift invariant (LSI) filter, h (x, s), with input signal i(x), r(x,s) - h(x,s)*i(x). In general, x is any continuous, M-dimensional independent variable, but we restrict our consideration to 1-D and 2-D cases. No restrictions are placed on the filter input, i(x), which is an arbitrary function of x. The LSI filter, h(x,s), is modeled as a family of spatial filters parametrized by the filter scale, s, where s is inversely proportional to filter bandwidth. The dependence of filter and operator outputs on filter scale is explicitly represented by r(x,s) and e(x,s), respectively. In standard scale-space theory, those x which satisfy e(x,s) = 0 (2.1) indicate positions where the Laplacian of the filtered s al is zero (i.e., L is the Laplacian operator), the x-s (hyper)plane is the scale-space, and the function e(x,s) is the scale-space image. The resulting (x,s) zero-crossing contours are referred to as fingerprints, due to the similarity of the shapes of these contours to the loops and whorls of human fingerprints. Scale-space theory establishes conditions under which all zero-crossings in e(x,s) are directly related to the zero-crossings of i(x), while minimizing the number of zero-crossings of i(x) that are eliminated by the smoothing effects of h(x,s)1. Yuille and Poggio [84,85] formalized these conditions over variable scale by requiring that the number of 1 These two conditions were examined from the standpoint of regularization by Torre and Poggio for a fixed scale [72] where it was shown that the Gaussian kernel satisfied both conditions for the Laplacian operator.

18 zero-crossings of e(x,s) increase monotonically as s varies from oo (zero filter-bandwidth) to 0 (infinite filter-bandwidth). They showed that three criteria must be met for this monotonic increase in zero-crossings with decreasing scale. {SS1} Thefilter impulse response h(x,s) must satisfy the following constraints: {SSla} h(x,s) is linear and shift invariant so that the definition of r(x,s) is valid, {SSlb} h(x,s) may be written in terms of another linear function p(.) as, h(x,s) = 1 p,) where x is M-Dimensional, so that the Fourier transform of h(x,s) is a linear function of scale with a d.c. response that is constant over all scale, {SSlc} limh(x,s) =(x), s to so that i(x) is exactly recovered as s approaches 0, {SSld} the center location of h(x,s) with respect to x is independent of s so that the value of s does not effect the spatial position of h(x,s), {SSle} lim h(x,s) = 0 and I h(x,s) j< oo for all s so that h(x,s) has Ixl-+o bounded volume. {SS2} The first and second derivatives of e(x,s) exist at all extrema. This condition holds for a continuous filter impulse response by virtue of the implicit function theorem.

19 {SS3} No extrema of fingerprints are minima. To enforce this condition for the 1-D case of x=x, it is sufficient to constrain the scale-space image, e(x,s), such that, e(x,s) =0 and e,(x,s) = 0 implies -ec(x, s) <0, e.(x,s) where, xs) e(x,s) e. (x, s) = —b-43X?e (x,s) e (x,) ae (xs) as In the 2-D case of x=(x,y), the scale-space image is constrained such that, e(u,v,s) = O and eM(u,v,s)=e(u,v,s) = 0, implies, -e.(u,v,s) -e.(u,v,s) -e(u,v,s) <0 and e(u,v <0 es(u,v,s) eS(u,v,s) where u - (u, v), and the u and v axes are taken along the directions of principal curvature in the (x,y) plane; u is a vector rotation x [37].

20 Given these three criteria, the unique class of filters h(x,s) that guarantees a monotonic increase in the number zero-crossings with decreasing scale is the Gaussian [4,84,85], h(x,s)= 2k x e where M is the dimension of x (M=1,2), and k is a constant. As a result, the fingerprints for a Gauss-Laplacian operator are smooth, continuous curves, and fingerprint extrema are never minima. Such fingerprints exhibit, what we call, an embedded structure with respect to scale. Figure 2. 1, depicting fingerprints for 1-D x, illustrates that these fingerprints have, at most, a single extremum per contour, and that these extrema are maxima. (0 Cn U) CD X Position Figure 2.1. Scale-space fingerprints having embedded structure for the 1-D case. Given the general multiresolution case in which the kernel violates the Gaussian assumption, we wish to know if is it possible to relax the condition on L{.} and/or the zero-crossing condition without severely restricting the utility of the scale-space analysis.

21 Relaxation of the conditions on L. } has been addressed by Yuille and Poggio [84,85] and Babaud et. al. [4] under the condition of Gaussian kernels. They show that standard scale-space results generalize to any operator L{.) that is a linear combination of integer-order differentiators. Further, Yuille and Poggio note that the zero-crossing criterion can be relaxed to include arbitrary level-crossings, and Geiger and Poggio [34] apply level-crossing to stereopsis. In this chapter, we consider scale-space analysis in greater detail for arbitrary function-crossings, and then address the issue of non-Gaussian kernels. 2.1.2 Remarks Similar to Yuille and Poggio [83,85], we are interested in 2-D Gaussian scale-space filters, h(x,s), that are rotationally symmetric and, therefore, separable in the (x,y) plane. Because of the advantages afforded in filter implementation and signal analysis, separability of 2-D filters and images are often assumed in image analysis [e.g., 2.24, 2.25]. Moreover, many remote sensing systems can be approximated as separable (chapter 1), and we exploit the advantages of filter separability in the development of this chapter. We note, however, that rotationally asymmetric Gaussian filters can be mapped into symmetric Gaussians, so that little generality is lost by assuming rotationally symmetric filters. A second remark concerns the effect of noise on scale-space filtering, where noise is the uncertainty between a filter input signal, i(x), and attributes of the object which generates it. Such noise can be treated as sensor noise, shown as ni for the i" sensor in Figure 1.3, and, for the purposes of scale-space filtering, is a component of i(x)2. For example, a "noisy" step-edge in computer vision is an ideal step-edge of image intensity with additive Gaussian noise [13]. Such edge noise could be introduced either by measurement uncertainties of the sensor or by imperfections in the underlying step-edge. The position of the underlying step-edge might be estimated by taking either first or second derivatives of i(x), and locating 2 The scale-space filter, h(x,s), is assumed to be noiseless.

22 extrema of first derivatives or zero-crossings of second derivatives. The effect of the noise, however, is to produce extrema, or zero-crossings, that are not associated with the location of the edge. Thus, i(x) must be processed (regularized [72]) in concert with differentiation to estimate the ideal step-edge location. The requirements for Gaussian scale-space filtering were originally developed deterministically. That is, the effects of noise were not explicitly considered. However, the edge detection work of Canny [13] and Torre and Poggio [72] showed that Gaussian filtering is nearly optimum for edge detection in the presence of noise. In this chapter, we also develop scale-space filtering results deterministically, and appeal to [13] and [72] to assure minimal degradation of system performance for signals in the presents of noise. 2.2 The Function-Crossing Problem Consider the case of an input signal i(x) = a (x) + n (x) (2.2) where a(x) and n(x) are the "known" and "unknown" components of i(x), respectively. For such a signal and L {. } = If{.} (I{. I is the identity operator), the fingerprint of the unknown component are derived from the expression, h(x,s)*n(x) =0, (2.3) provided n(x) can be extracted from the input signal i(x). Alternatively, the fingerprints of (2.3) can be approximated by fingerprints derived from, h(x,s)*i(x) = i(x,s)*a(x) (2.4)

23 where fi(x,s) is an estimate of the impulse response of the scale-space filter h(x,s). The use of a "good" estimate of h(x,s) will result in a close match between the fingerprints of (2.4) and those of (2.3). However, this may present problems in those multiresolution systems where h(x,s) is not accurately estimated for all sensor classes. Moreover, a(x) may also be estimated, and (2.4) have to be re-calculated for numerous a(x) estimates. Thus, we consider alternative approaches for approximating the scale-space fingerprints of (2.3). One such approach approximates the fingerprints of (2.3) by fingerprints derived from, h(x,s)*i(x) = d(x) (2.5) where d(x) is an estimate of a(x). That is, the true fingerprints could be approximated by fingerprints based on crossings of the function d (x). The extent to which this approximation holds is reflected, in part, by whether the fingerprints of (2.5) satisfy the standard scale-space properties, as reviewed in ~2.1. Appendix A presents a somewhat different, more general derivation for the conditions under which the embedding property of fingerprints is met exactly. Included are the cases of a general function-crossing and an arbitrary linear operator O {f.), rather than L {.). These conditions are summarized by (2.6a), (2.6b), and (2.7), a b d 2 y(x)+-y'(x)+cy"(x)-cm2 d 0 (2.6a) S S S for es(x,s)=l, a b 2 (2.6b) 'Y(x)+ y'(x)+cy"(x) + cm2 +- (2.6b) S S S for es(x,s)=-l, and

24 a b d a-h(x,s) +-h (x,s)+ch,(x,s) -h,(x,s)= 0. s s s (2.7) That is, embedded structure exists if (2.6a), (2.6b), and (2.7) are satisfied, where y(x) is a desired thresholdfunction, e(x,s) is given by e(x,s) = O{h(x,s)*i(x)}, m is an arbitrary constant, and constants a, b, c, and d are to be determined. 2.2.1 Analysis Equation (2.7) can be solved as a generalized heat equation (d ~ 0) [84,85], resulting in constraints on the coefficients, a=b=0 cd >0, and a Gaussian filter, h(x,s) - 4( (XS)2 h(x,s)= I- I= e _ J VS24) (2.8) From (2.6a) and (2.6b), the coefficient constraints require the threshold function, y(x), to satisfy, cy "(x)-cm2 —~0 for e,(x,s)=1 S d cy"(x)+cm2 +-0 for e,(x,s)=-l. S (2.9a) (2.9b)

25 For T(x) to satisfy (2.9a) and (2.9b), it is necessary that, y"(x)<O for e,(x,s)= 1 (2.10a) y"(x)>O for e,(x,s)=-1. (2.10b) However, obtaining a general function y(x) that satisfies (l Oa) and (1 Ob) for all input signals, i(x), may be impossible. To obtain a function y(x) that satisfies (2.10 a) and (2.10 b), and is independent of i(x), the sign of es(x,s) must be independent of s. Expressing es(x,s) as, e,(x,s) =s ({J h(y,s)i(x -y)dy} =j[ (h(y,s))]O{i(x - y)}dy, a Gaussian h(x,s) yields, e.(x, s) = (s2-1 s)e - O{i(x-y)}dy. (2.11) Thus, (x) is independent of s only if the sign of e,(x,s), given by (2.11), is independent of s. The difficulty in generating such an es(x,s) is seen by considering the relatively simple example of an intensity image (i.e., i(x) > 0) and the identity operator (O{i(x) }=i(x)). In this example, all sign changes in es(x,s) are introduced by h(.,.), but the sign of es(x,s) will still have an s-dependence unless,

26 y< S Thus, the maximum extent of O{i(x-y)}, in (2.11), must be less than the spatial width of the highest resolution filter, which eliminates virtually all reasonable inputs.3 As a result of the coupling of y(x) to the sign of e,(x,s) and to s, general threshold functions must be input-specific. This coupling can be avoided directly by using a constrained y(x), and forcing y"(x) = 0. In this case, either a constant threshold, y(x)= ro (2.12a) or a threshold that is a linear function of x, y(x) = rllx + rio (2.12b) satisfies the conditions of (2.9a) and (2.9b). Thus, threshold functions, y(x), given by either (2.12a) or (2.12b) are of practical use with an arbitrary input signal i(x). The alternative method for de-coupling y(x) from s, is to consider (2.7) not as a generalized heat equation, as before, but as a differential equation with scale-dependent coefficients (i.e., forcing d=O in the original constraints on the coefficients). By substitution of d=O into (2.7), a second-order differential equation results, a b h(x,s) +-h(x,s)+ ch,(x,s) = 0, (2.13) S S 3 The use of a more general operator introduces additional sign changes on O{i(x-y)) and does not eliminate the s-dependence of the sign of e,(x,s).

27 which has the solution, 1 bx h(x,s) e 2c 2+ ~ ac e +e In (2.13), the values of variables x and s are independent, so that x may be an arbitrarily large positive or negative real number relative to some positive real-valued s. To maintain a bounded h(x,s) for all x, it is necessary that, bx e2c <oo Since the constants b and c are independent of x and s, the above inequality is satisfied for all x and s only for b=O and c ~ 0. Therefore, the solution to (2.13) reduces to, h(x,s)=-(e s e j S \ + e The boundedness of h(x,s) also requires that (x/s)(a/c)'/ be real, so that a/c>0 and the solution to (2.13) becomes, 2 (xr h(x,s) =-cosl- a s s ) with the coefficient constraints, b=d=O ac >0. From (2.6a) and (2.6b), the effect of these coefficient constraints on the threshold function again results in a coupling between y(x) and the sign of es(x,s), for all y(x). Thus, no threshold function can satisfy (2.6a) and (2.6b) in this case.

28 In summary, only Gaussian solutions of (2.7) yield a threshold function which can generate fingerprints with embedded structure for general input signals. Moreover, the only forms of threshold function, y(x), that are of practical use with general inputs are the constant or linear thresholds given by either (2.12a) and (2.12b), respectively. Although it may be possible to construct highly constrained input signals which allow the use of nonlinear thresholds for Gaussian solutions to (2.7), such input signals may be so strongly constrained that scale-space analysis becomes merely an exercise in confirming that which is already known. As a result, if function-crossings are to be used in analyzing the fingerprints of signals with known and unknown components, the threshold functions should consist of linear and constant segments and the scale-space filter should be Gaussian. 2.2.2 Remarks The results of the previous section can be demonstrated using a simplified remote sensing example, where the objective is to detect frozen ground within a warmer land surface using antenna temperature (radiobrightness) at a single frequency. As discussed in chapter 4, physically colder surfaces emit less energy so that where effects of temperature change on the dielectric constant of the surface is relatively small, frozen surfaces can be located where antenna temperatures fall below a specified threshold. The antenna temperatures for a 1-D surface of relatively warm land containing a frozen region is simply modelled in Figure 2.2a4. Figure 2.2a shows the surface antenna temperature as a function of surface position, x. The warm land surface extends from -400 Km to +400 Km and, except for the frozen region, has a constant antenna temperature of 245 K. The frozen region extends from -80 Km to +80 Km, and has a minimum antenna temperature of 235 K at x=0. Also shown in 4 The parameter values used in the model are derived from frozen terrain results described in chapters 4 and 5.

29 Figure 2.2a is the constant threshold function used to detect the presence of frozen ground; the threshold value is 244 K. The resultant fingerprints from spatial Gaussian filtering are shown in Figure 2.2b. Now consider the case of frozen terrain contained in a relatively warm surface with non-constant antenna temperature. In the case of Figure 2.3, the antenna temperature of the warm surface has a known, or easily estimated, decrease in the vicinity of the frozen surface. Such an antenna temperature decrease could be caused either by a decrease in the physical temperature of the terrain or by an increase in soil moisture. The antenna temperature for such a surface is shown Figure 2.3a. Since the objective in this example is to detect changes in the state of the surface moisture, which are identical in Figures 2.2a and 2.3a, we seek a threshold function that eliminates the effects of the surface variations outside the frozen region. That is, the effects of freezing on radiobrightness form the n(x) signal of (2.2), and the radiobrightness variations outside the frozen region for the a(x) signal. A threshold function is sought where fingerprints of the antenna temperature of Figure 2.3a approximate fingerprints shown in Figure 2.2b. Using the constant threshold function of Figure 2.2a (repeated in Figure 2.3a), Gaussian filtering of the antenna temperature of Figure 2.3a produces the fingerprints shown in Figure 2.3b. Comparing Figure 2.2b with Figure 2.3b, it is seen that the fingerprints match at fine scales, but not at coarse. To better match the fingerprints of Figure 2.2b, we use the non-constant threshold function shown in Figure 2.4a.

30 260.0.... 257.0 254.0 251.0 248.0 S 245.0 ^ -------------------------- - — y - —. --- —----------- Threshold 242.0 239.0 236.0 233.0 2 3 0.0................................................ -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 240.0 320.0 400.0 X (km) Figure 2.2a. Input signal and threshold for warm terrain, at constant antenna temperature, and frozen region. 160.0...... 150.0 140.0 130.0 -. 120.0 110.0 100.0 90.0 _q 80.0 70.0 60.0 50.0 40.0 30.0 20.0 10.0 0.0.......... 400.0-320.0-240.0-160.0 -80.0 0.0 x (km) 80.0 160.0 240.0 320.0 400.0 Figure 2.2b. Resulting fingerprints.

31 260.0...... 257.0 2154. O - 251.0 E 248.0 | 245.0 S 2450: --- —-----—. —.- ^ -- --.- -.- -- 242.0 Threshold X 239.0 236.0 233.0 230.0 -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 240.0 320.0 400.0 X (km) Figure 2.3a. Input signal and constant threshold for warm terrain, with varying antenna temperature, and frozen region. 160.0........ i...............;.... 150.0 / 140.0 130.0 120.0 110.0 -s 100.0 -^ 90.0 Q) 80.0 c) 70.0 60.0 50.0 40.0 30.0 20.0 10.0 -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 240.0 320.0 400.0 x (km) Figure 2.3b. Resulting fingerprints.

32 260.0...................... 257.0 254.0, 251.0 248.0 245.0 Threshold 242.0 239.0 236.0 -233.0 230.0.......tI....... l l.............. 1" -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 240.0 320.0 400.0 X (km) Figure 2.4a. Input signal and piecewise linear threshold for warm terrain, with varying antenna temperature, and frozen region. 160.0 150.0 140.0 130.0 120.0.... 110.0 / 100.0. 90.0 80.0 70.0 60.0 50.0 40.0 30.0 20.0 10.0 0.0 0.0..............!.......,.......................... -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 240.0 320.0 400.0 x (km) Figure 2.4b. Resulting fingerprints.

33 In Figure 2.4a, the threshold function is composed of linear segments that approximate the antenna temperature variation of the warm land surface. For this example, the linear segments are connected so that the threshold function, y(x), is C~ continuous. The antenna temperature curve of Figure 2.4a is a copy of the antenna temperature curve of Figure 2.3a. The resulting fingerprints for the antenna temperature and threshold function of Figure 2.4a are shown in Figure 2.4b. We see that at coarse scales, the fingerprints of Figure 2.4b better match the fingerprints of Figure 2.2b, than those of Figure 2.3b. However, we also see fingerprint approximation errors in Figure 2.4b as fingerprints cross over the connections between linear segments in the threshold function of Figure 2.4a. As a result, the fingerprint structure of the unknown component, n(x), of a signal, i(x), may be derived from the fingerprints of, h(x,s)*i(x) =y (x), (2.14) where, i(x) = a(x) + n(x), y(x) = (x), da(x) is an estimate of a(x), and da(x) is a piecewise linear function of x. However, the degree of correspondence between the fingerprints of (2.3) and (2.14) depends on the number of linear segments that comprise y(x), as well as on the accuracy with which a(x) can be approximated by combinations of linear and constant functions. Previously, it was shown that Gaussian scale-space filtering produces fingerprints with embedded structure for crossings of linear threshold functions. However, embedded structure does not necessarily exist at the knots between linear segments. A breakdown in embedded structure at knots results in approximation errors between the fingerprints of (2.14) and those of (2.3). Such a breakdown is evident in the fingerprints of Figure 2.4b.

34 Despite breakdowns in embedded structure at knots between linear segments of y(x), the resulting fingerprint approximation errors may not impose severe limitations. First, the location of the knots in y(x) are known apriori. Thus, the locations of potential approximation errors in the fingerprints of (2.14) are known and may be avoided. Second, the fingerprints of (2.3) have known behavior (i.e., they have embedded structure). Thus, it may be possible to interpolate the fingerprint of (2.14) over the regions of connecting points to reduce approximation errors. Third, the form and, subsequently, the level of approximation errors depend on how linear segments of y(x) are connected. For example, if the linear segments are joined at points so that y(x) is C~, then the implicit function theorem breaks down at connecting points and the fingerprints of (2.14) may have discontinuities. If, however, the linear segments of y(x) are smoothly connected so that y(x) is Cn for n > 1, then the fingerprints of (2.14) will also be smooth. Thus, it may be possible to use connecting functions between linear segments of y(x) that minimize approximation errors. 2.2.3 Two-dimensional Extensions Previously, the fingerprints of, O{H(x,s)*i(x)} =y(x) (2.15) for one-dimensional signals were shown to have embedded structure if h(x,s) is Gaussian, Of.} is a linear operator, and the threshold function is given by, y(x)=rrx +ro, where r, and ro are constants. This result can be extended to two dimensions, so that the fingerprints of

35 O{h(x,y,s)**i(x,y)} e(x,y,s) =(x,y) (2.16) have embedded structure if O{.} is a 2-D linear operator, h(x,y,s) is a 2-D Gaussian filter, 1 x + y2 h(x,y, s)= 2k e. 2 2itks2 and the threshold function is given by, ^x, y) = ro + rx + r2y + rYy. The coefficients ro, rl, r2, and r3 are constants, and the constant k is non-negative; "**" represents 2-D convolution. The linearity requirement for the 2-D operator, Of{.}, follows directly from the same arguments as used for 1-D signals, which are not repeated here. The Gaussian filter requirement follows from the arguments of Yuille and Poggio. However, to see that the above threshold function, y(x,y), generates fingerprints with embedded structure, it is necessary to repeat previous derivations for 2-D signals. Following the approach of Yuille and Poggio, the fingerprints of (2.16) have embedded structure if, e(u,v,s) = y(u,v) (2.17a) e,(u, v,s) = yW(u, v) (2.17b) e,(u, v,s) = ^,(u,v) (2.17c) implies,

36 y~ (u,v)- e,(u,v) -e(uv,s) <O e,(u,v,s) Yw(u, v)- e,(u,v,s) <0, e3(u,v,s) (2.17d) (2.17e) are satisfied, where the u and v axes are taken along the lines of principal curvature. Equivalently, coefficients a, b1, b2, c1, c2, C3, and d must be found such that, a bl b2 d 2e +- e +- e + cle, +c2e +c3eu + es = 0, S S S S (2.18) and, a bl b2 2 d a7+ -y+-y, + c2y + c y + c2Ym - c2m2 - 0 for e, = 1 a b1 b2 2 2 d "Y +,-v + Cluu +C2Yw + C3uv cm + C2M2 +- ~0 for e =-1, S S S S (2.19a) (2.19b) where m, and m2 are arbitrary constants. The Gaussian filter of (2.16) is generated by, a = bl = b2 = 3 = 0 and, cc2 > 0 and cdd > 0. Using these coefficients, (2.19a) becomes,

37 c1.- + czyw - cumC - c2m2 - - 0 (2.20a) and (2.19b) becomes, Z d c + c2Y M + c2mc + c - 0. (2.20b) s As a result, a threshold function given by y(u,v) = q + qu +q2v +q3uv, (2.21) where qj, q2, and q3 are constants, solves (2.20a) and (2.20b) and generates fingerprints with embedded structure. Because the u and v axes are taken along directions of principal curvature, they are orthogonal [37] and there exists an invertable matrix A such that the u and v axes can be mapped onto the x and y axes, [u =Apx]. Thus, a threshold function given by, (x,y) = ro + rx + ry + rxy will give rise to threshold functions given by (2.21) for all u and v, and therefore, generate fingerprints with embedded structure. Two remarks are warranted. First, the use of a symmetric Gaussian filter h(x,y,s) in generating fingerprints of (2.16) with embedded structure produces no loss in generality,

38 because effects of an asymmetric Gaussian filter can be produced by the use of an asymmetric operator, 0 (. }, on a symmetric Gaussian filter. Second, while the general form of threshold function, y(x, y), is given above, perhaps a more useful form of threshold function is that of a plane, y(x, y) = ro + rlx + r2y. Therefore, the use of planar threshold functions allows a desired threshold function in 2-D to be approximated by combinations of planar surfaces. 2.3 Non-Gaussian Filtering Multiresolution systems commonly use sensors with apertures to collect received energy. Since apertures are always of finite size, they are mathematically equivalent to bandlimited filters [10] and, therefore, cannot be Gaussian. This raises questions of the extent to which the properties of scale-space filtering obtained in the Gaussian case generalize to the non-Gaussian case. In the sub-sections that follow, we consider the scale-space fingerprints generated by non-Gaussian filters. In the standard development of scale-space theory, as reviewed in ~2.1, it was shown that non-Gaussian filters cannot generate fingerprints with embedded structure. In [81] it was shown that non-Gaussian filtering can be used provided filter inputs are constrained to be polynomials. Alternatively, we show that if the standard constraints on the filter impulse response are relaxed, then non-Gaussian filters can generate fingerprints with embedded structure. We also show that while non-Gaussian fingerprints do not match Gaussian fingerprints at all scales, the non-Gaussian fingerprints are either related to Gaussian fingerprints, or match them exactly over a range of scales.

39 Two approaches are considered for generating fingerprints with embedded structure for non-Gaussian scale-space filters. In ~2.3.1, a spatial-domain approach of nonlinear scaling is considered for generating non-Gaussian fingerprints. By relaxing the impulse response condition { SS lb }, we show that the resulting fingerprints have embedded structure. In ~2.3.2, fingerprints from scale-space filters with limited frequency support are examined. By the relaxation of condition { SS Ic), we show that these fingerprints also have embedded structure. 2.3.1 Nonlinear Scaling Consider the case of non-Gaussian filters in which the kernels are functionally related to Gaussian kernels by continuous mappings in the spatial and scale domains. One such filter is the exponential, h, (x, s) - a(s)e-xl/s where a(s) is some function of scale, s. The exponential filter he(x,s) can be transformed to the Gaussian, go(u,t), 1 -2(u/t)2 go(, t) -- e(2 by the mappings, s =2t2 u2 for u <OJ -(u2) for u <J where,

40 a(s) = 1 g That is, the exponential filter is a Gaussian filter that has been "warped" in the spatial and scale domains, and is a nonlinearly scaled scale-space filter5. More generally, if there exists a non-Gaussian filter h(x,s) and continuous, monotonic mappings X(.) and S(.), (i.e., X-1(.) and S-'(.) exist6), and if h(x,s)=g(X-(x),S-'(s))=g(u,t), (2.22) where g(u,t) is Gaussian, 1 _ —(u/t)2 g(u,t)= — e then the non-Gaussian fingerprints given by, then the non-Gaussian fingerprints given by, (x)= h(x,s)*i(x), (2.23) are warped versions of the Gaussian fingerprints given by, f(u) = g(u, t)*j(u). (2.24) In (2.22), k is a positive constant, and 5 The case of a(s) = l/s is considered linear scaling due to the frequency response of he(x,s). 6 The existence of X-'(.) and S-'(.) is easily shown using iterations on the chain rule, and is not developed here.

41 ( = (x-'()) y(x) j(u) =j(X-(x))= i(x). From the results of [84,85], the fingerprints of (2.24) have embedded structure with respect to variables u and t. Since X(.) and S(.) are monotonic and invertable, the fingerprints of (2.23) have embedded structure with respect to variables x and s, even though h(x,s) is non-Gaussian. In addition, only the threshold function 8(u) must meet the linearity constraints of ~2.2 for embedded structure to occur. Thus, y(x) may be nonlinear, but will depend on the mapping X(.). The exponential filter clearly meets the requirements on X(.) and S(.), and would produce fingerprints with embedded structure with respect to u and t. However, standard scale-space theory [84,85] is predicated on linear scaling (condition {SSlb}), so while scale-space theory admits embedded structure for the linear fingerprints of (2.24) (with respect to u and t), it precludes embedded structure for the nonlinear fingerprints of (2.23) (with respect to x and s). For scale-space theory to admit embedded structure for the fingerprints of (2.23), the linearity condition {SSlb) must be relaxed to accept nonlinear scaling. 2.3.2 Limited-support Filtering As opposed to warping Gaussian scale-space filters in the frequency domain, truncating Gaussian filters in the frequency domain will generate non-Gaussian filters having fingerprints with embedded structure. Scale-space filters with limited frequency support, or bandlimited filters, are of particular interest in multiresolution systems since sensor apertures exhibit these types of characteristics. Two types of bandlimited filters are considered: fixed-support filters and relative-support filters.

42 2.3.2.1 Fixed-support Filters In fixed-support filters, the bandlimits of the scale-space filter are fixed and independent of scale. In Figure 2.5, the frequency response of a fixed-support scale-space filter is depicted, where the bandlimit extends over frequency from -fo to +fo. Fingerprints generated by fixed-support filters differ from those generated by ideal Gaussian filters, since fixed-support filters do not have the infinite frequency-domain support. However, fixed-support filters generate fingerprints with embedded structure. The solution to (2.7) that generates fingerprints with embedded structure has a frequency response, H(f,s), given by7, H (f,s)=forfE H(f,s)= for f e Sc where, G(f,s)-e k is a constant, and S is any set of real numbers. As a result, the fixed-support filter solution to (2.7) can be given by, H(f,s) = G (f,s)nf, (2.26) 2foj where II(.) is the rect function, 7 This result is derived from the work Yuille and Poggio [84,85]. However, they consider only the Gaussian solution to (2.7).

43 n(f ) for wise ^2/0) J 0 otherwise J In the spatial domain, h (x, s) = g (x, s)* [2fosinc (xfo)] (2.27) where, g (x, S) - 2(xs)2 so~n Figure 2.5. Family of fixed-support filters; filter response flattens with decreasing s. Fixed-support scale-space filters are an implicit, but generally ignored, component of the original scale-space filter results in computer vision. In the perspective of computer

44 vision, the optical systems that provide images for vision processing also provide the rectangular filtering of (2.26). It is usually assumed that the bandwidth of such rectangular filtering is adequate to ensure that features (edges) in the images being processed correspond to features in the physical scene. This is not necessarily the case in multiresolution processing, so that (2.26) must be used in analyzing features that exist in signals prior to sensor processing. The embedded structure results of [84,85] do not include the fixed-support scale-space filter due to the limiting value of h(x,s) in (2.27), where, limh(x,s) = 2fosinc(xfo) ~ 8(x). (2.28) s — o The impulse response condition {SS Ic} is violated by fixed-support filtering; {SS ic} must be relaxed to admit the use of fixed-support scale-space filters. In the perspective of standard scale-space filtering, if the "true" fingerprints of a signal i(x) are given by, g(x,s)*i(x) =y(x) for Gaussian g(x,s), the fingerprints given by, h (x,s)*i(x) =y(x) are identical to those of g(x,s)*i(x) =y(x)

45 where i(x) is a "corrupted" version of i(x), i(x) = i(x)*[2fosinc(xf,)]. If fo is large relative to the bandwidth of the input signal i(x), the corrupted fingerprints generated by the fixed-support filter of (2.27) may approximate those generated by an ideal Gaussian filter. 2.3.2.2 Relative-support Filters In relative-support filtering, the bandlimits of the scale-space filters depend on scale, and the shape of the filter is unchanged as a function of scale. Such a scale-space response would be characteristic of a multiple-frequency, single-aperture system. In Figure 2.6, the frequency response of a relative-support filter is depicted. The frequency response of a relative-support filter is given by, H(f,s)=F(f, s) (2f- (2.29) where F(f,s) is some filter function of interest that approximates the Gaussian over some limited region of space and scale. In multiresolution systems, F(f,s) is an aperture illumination function, which often allows the impulse response condition {SSlb) to be satisfied (see [74] for example), lim h(x, s) = (x). s -,o

46 In order to have a non-Gaussian filter, as described by (2.29), yield fingerprints with embedded structure, a double threshold-crossing (i.e., a hysteresis process) can be used to calculate threshold-crossings. Canny [13] successfully demonstrated the use of threshold hysteresis to reduce the effects of edge noise. Lu and Jain [51] apply threshold hysteresis to scale-space filtering, but do not discuss it in detail. In our application of threshold hysteresis, truncation errors, e, are determined from the difference between the relative-support (truncated) filter and ideal Gaussian responses. These truncation errors are then used to generate two threshold levels for defining fingerprints: {y(x) + ~} and {y(x) - E}. Figure 2.6. Family of relative-support filters; filter bandwidth broadens with decreasing s.

47 To determine an expression for the truncation error, a difference function between h(x,s) and an ideal Gaussian impulse response, g(x,t(s)) is calculated, b (x, s) h (x, s)-g(x, t(s)). (2.30) The variable t(s) is the variable of scale for g(.,.), and is a function of the scale variable s; the scales of h(.,.) and g(.,.) are functionally related. At each value of s, a function g(x,t(s)) is sought that minimizes the square error, E -J b(x,s) 2dx. Once g(x,t(s)) is found, the minimum E is used to bound the truncation error e, and to define the thresholder. The thresholder is determined from the function-crossings of a filter output r(x,s) with r(x,s) = i(x)*h(x,s), (2.31) for a filter input function i(x), a filter impulse response h(x,s), and some threshold function y(x)8. The thresholder defines a negative-goingfunction-crossing of r(x,s) to occur at x0 if, r(xo,s)=y(xo)+e and 3x, >x0 3 r(x, s)=y(x )-E and y(x)-e<r(x,s) <(x)+e Vx e (Xo,xl). (2.32a) 8 Without loss of generality, the operator O{.} is the identity (appendix A)

48 A positive-going function-crossing of r(x,s) is defined to occur at Xo if, r(x, s)=)y(Xo)-E and 3xl>xo 3 r(xl,s)-=(xl)+~ and y(x) - < r(x,s) < y(x)+ e Vx e (Xo, X). (2.32b) To guarantee that no "false" threshold-crossings occur, the truncation error ~ used in the thresholder can be given by i( |b(x,s)l2dx i(x-y)b(y,s) dy, (2.33) where, l/p 11 i(xl - i(x)lPdx and b(x,s) is the function that minimizes E. By excluding all false threshold-crossings, function-crossings given by (2.32a) and (2.32b) are also function-crossings between y(x) and a Gaussian filtered i(x), i(x)*g (x, t(s)) = y(x). (2.34) Thus, as c -> 0, the number and location of function-crossings given by (2.32a) and (2.32b) and those given by (2.34) become identical.

49 In light of the definition of the relative-support filter (e.g., (2.29)), the Fourier transform of (2.30) gives, B(f,s) - H(f,s)- G(f, t(s)) (2.35) where B(.,.), H(.,.), and G(.,.) are the spatial Fourier transforms of b(.,.), h(.,.), and g(.,.), respectively. By the Rayleigh Energy Theorem, the function g(x,t(s)) that minimizes the square error is the Fourier transform of the function G(f,t(s)) that minimizes, E= flB(f,s)I 2df Thus, the minimization of e can be performed in the frequency domain, where truncation limits of relative-support filters are defined. For example, consider a parabolic relative-support filter9, H(f)= -(fs)2 for I fS 1~1 0 otherwise J where E is given by, 1 56 2 1 2 1/2 e v 15 2t/s tipas (t/s)3 r ilerf/s) (t/s)2 ex 9 Some sensor (antenna) models used in microwave remote sensing are mathematically equivalent to parabolic relative-support filters; see [74] for example.

50 and, -< II I( fI I(f)l PdfJ I(.) is the spatial Fourier transform of the input signal, i(.). Thus, the truncation error for minimum E is dependent on the power of the input signal and the scale of the filter, 1.9 ~ III(fA 2 2.3.3 Two-dimensional Extensions Previously, one-dimensional, non-Gaussian scale-space filters produced by nonlinear scaling and bandlimiting required the relaxation of several impulse response condition (~2. 1) to generate fingerprints with embedded structure. The 2-D filter formulation of nonlinear scaling follows directly from the 1-D arguments. If there exists a separable 2-D map X(.,.), X(x,y) = (X,(x),X2(y)) (u,v) (2.36a) and a map S(.), S(t)=s (2.36b) such that X(.), X2(.), and S(.) are monotonic and continuous, then a filter h(x,y,s) can be expressed as,

51 h(x,y,s)= g(X,(x),X2(y),S-l(s))= g(u,v,t) (2.37) where [84,85], (21 kt2 J g(u, v,t) e t2 k>O and the fingerprints of, h(x,y,s)**i(x,y) = y(x,y) (2.38) will have embedded structure. As in the 1-D case, the scale-space linearity requirements of {SSlb) must be relaxed to admit the fingerprints given by (2.38). The separability of X(.,.) is a direct extension of the separability of 2-D Gaussian scale-space filters, as given by (2.37). The 2-D filter formulation of limited support filters follows directly from 1-D1) arguments for fixed-support and relative-support filters. For fixed-support filters, a filter with 2-D frequency response H(f,,fy,s) will produce fingerprints with embedded structure if, H(AVf.,s) = {GO,fy,s) for (f,,f)E s} f,ls) for (ff) e S } where, G(f fy, s)-e

52 k is a constant, and S c 9R2. As a result, a filter with frequency response, H(ffy,s) = G(fxfy,s)R(ffy), (2.39) where, {1 for Ifxl- <f1fyl <f2} 0 otherwise for specified f, and f2, will generate fingerprints with embedded structure, provided condition {SSlc} is relaxed. The fixed-support filter H(f,,fy,s) could have rectangular or elliptical support, amongst others, in the 2-D frequency domain. For relative-support filtering, the frequency response H(ft,fy,s) can be expressed as, H(ffs) = F(ffs)R(sfsf), (2.40) where F(ffy,s) is a given filter function. F(ffy,s) is usually determined by the physical devices being used and, as in remote sensing applications, is often separable. As in the 1-I) case, hysteresis in the thresholder must account for the non-Gaussian form of H(fx,fy,s). The level of threshold hysteresis is determined by the truncation error, e, derived from a minimum mean-square error, E - lB (f(,f,s)l 2dfxdfy, and, B(f,,fys) = H(f,,,,s) - G(ffyt(s)). (2.41)

53 The calculation of E is made using a 2-D, frequency domain form of (2.33). In fixed-support and relative-support filtering, the use of separable filters, G(.,.,s) and F(.,.,s) respectively, requires that R(.,.) be separable to make the resulting filter H(.,.,s) separable; see equations (2.39) and (2.40).

CHAPTER 3 MULTIRESOLUTION PROCESSING In this chapter, we present a multiresolution system that is derived from multispectral processing, and show that boundary estimates made at coarse resolution can have significant localization errors. We also present an estimation scheme that is derived from scale-space filtering and can mitigate the localization errors. Limitations of this approach are also developed. Both boundary migration and estimation problems are posed with respect to 1 -D contours along a 2-D surface. In addition, some extensions for more general 2-D signals are discussed. 3.1 Background As mentioned in chapter 1, scale-space filtering can provide a mathematical framework for data integration in multiresolution systems where sensors vary along a dimension called scale. That is, if signal differences between sensors can be characterized exactly by PSF (width) differences between sensors, then data from each sensor represents scale samples of the scale-space representation of the imaged object. In this chapter, we extend this scale-space formulation to multiresolution systems where PSF is not the only characteristic that differs between sensors. 54

55 We will be concerned with the scale-dependent locations of threshold crossings, x(s), calculated from e(x(s),s) = (x(s)), (3.1) where, e(x,s)- O{r(x,s)}, r(x,s) h(x,s)*i(x), O{.) is a linear operator, and y(.) is a specified threshold function (this process is modelled in Figure 3.1). Without loss of generality, 01. is the identity operator so that e (x, s) = r (x, s). Standard practice is to denote the threshold crossing contours (fingerprints) in scale-space as s(x), since scale-space is typically displayed with the scale axis as the ordinate (Figure 2.1). We also display scale-space with the scale axis as the ordinate, but our theoretical developments are concerned with threshold crossing locations as functions of scale and the more appropriate notation, x(s), is used. Vertical fingerprints in scale-space displays correspond to cases where x(s) - 0, and horizontal fingerprints in scale-space displays correspond to x(s) - oo. Input Filter Output Operator Output Linear Filter so- Operator.. Threshold i(x) h(x,s) r(xS.) } W) e(x,s) Figure 3.1. Scale-space system diagram.

56 In chapter 2 we showed that the original requirements of the theory can be relaxed and still preserve major points of scale-space theory. However for analytical convenience, we assume Gaussian filtering. Given the use of Gaussian kernels, h(x,s), Gaussian filtering can be used to degrade the resolution (broaden the PSF) of fine-resolution sensors in multiresolution systems to match that of a coarse-resolution sensor prior to processing. Alternatively, given the features of a signal measured by a coarse-resolution sensor, it is useful to register their locations with signal features measured by fine-resolution sensors. These two problems represent fingerprint interpolation and extrapolation, respectively, with respect to scale. In the sections that follow, we present a formal development of fingerprint extrapolation in scale-space. 3.2 Boundary Migration In multiresolution processing, regional boundaries can be estimated from threshold crossings of what we call an indicator function il(.). Indicator functions are derived from physical models so that regional boundaries occur between regions having different ii(.) values. That is, regional boundaries occur at x values where ii(.) crosses a threshold, i,(x) = a(x), (3.2) and ac(.) is a specified threshold function1. Indicator functions are commonly used in remote sensing (see [67], for example) and are often composed of linear combinations of single sensor data, 1 a(.) is a linear function of x.

57 M 6(.) - ajij(.), j=l where ij(.) is the ideal output signal from the jd sensor, and the aj's are coefficients. The function ij(.) would be the output from the j' sensor if it had infinitesimal (ideal) spatial resolution. In practice, ii(.) cannot be thresholded directly due to the finite (non-infinitesimal) PSF' s of the sensors. Instead, boundary locations are estimated by threshold crossings given by, M X airi(x,sj) = a(x), (3.3) j=1 where rj(x, sj) is the actual output of the jh sensor, ri(x, sj) h(x,s s)*ij(x). When resolution compensation is used, the effective scale (resolution) of all sensors are equalized and si = s, Vj, where s, is constant. Thus, (3.3) becomes, r,(x(s),s,) = a(x(s,)), (3.4) where, M r,(x(s),s,)- ~ ajr,(x(s,),s,), j=1 and the threshold crossing locations given by (3.4), x(s, are the regional boundary estimates. The kernel h(.,s) approximates a delta function as s -e 0 and, in the absence of noise, the

58 threshold crossings of ri(.,s) approach those of ix(.). To optimize estimates of boundary locations, we seek x(se) for as small a scale as possible. If sj is the (finest) scale at which data from the j channel is available, and s, <... < sM, then x(s,) can only be determined for s, > SM without parametric or model-based estimation. As a result, x(se) will be in error from the actual boundary locations, x(O), and first-order bounds on these estimation errors can be derived from boundary migration analysis. 3.2.1 Analysis Consider the case of threshold crossing contours, x(s), given by2, r(x(s),s) =0, (3.5) where, r(x(s),s) h(x(s),x)*i(x(s)). The true boundary location of i(.), x0, is given by, r(xo, O)= i(xo)= O, so that the boundary location estimate, x(s), migrates from xO as a function of s. This migration, D(s), can be approximated using a Taylor expansion as, D(s) = -r ) (3.6) 2 To simplify the analysis, the threshold function in (3.1) is set to 0.

59 where, D(s) -x(s)-Xo. A bound on D(s) can be determined from, {h(xo-y,s)i ID(s)l =I (y)dy (3.7) h(xo - y,s)i(y)dy and where Gaussian sensor PSF's are given by, 1 — (z/s)2 h(z,s)=g(z,s)- e 2 Even (symmetric) and odd (anti-symmetric) components of i(y) dictate the behavior of the boundary as a function of resolution [50]. We use a signal model given by, i(y) = w(x- y)io(x- y)+ [1l-w(o - y)] i,(x0-y), (3.8) where, io(.) - odd component of i(.) i,(.) - even component of i(.) w(.) - even weighting factor, and, 0~<w(.)~< 1.

60 As a result, (3.7) can be written as, g(-zs) [1-w(z)] ie(z)dz 1D(s)1 =- -0 f g(-z, s)w(z)io(z)dz fe(z [1 - w()] ie()dz =00 (3. 9) z -(Z/s)2 - e 2 w(z)io(z)dz so that a first-order migration bound can be given by, * 1<~ \ e-2(z's) 1 - w(z)] dz i (z j ID(s)I <, (3.10) r e 2(z/)w (z)io(z)dz 0 where, llf(xl - sup If(x)I. We are concerned with anti-symmetric boundaries, so that io(.) is guaranteed to exist at a boundary. Ifie(.) also exists, (3.10) shows that boundary migration occurs. The weighting function, w(.), controls the mixture of io(.) and ie(.) at the boundary. We consider migration bounds for two types of weighting functions. In both weighting functions, ij(.) is excluded close to the boundary, and io(.) is excluded away from the boundary. However, one weighting function makes a sharp distinction between these odd-only and even-only conditions, while the other is more gradual.

61 Equation (3.10) also shows that the shape of io(.) determines boundary migration. We consider migration bounds for two shapes of io(.). In both shapes, io(.) is constant at each side of the boundary. However, in one case the transition from one region to the other is piecewise linear, while in the other case the transition from one side of the boundary to the other is smooth. Case 1: The weighting function is given by W (){ for Iz 1<L 0 otherwise, J so that io(.) exists only in a well-defined neighborhood about the regional boundary, and (3.10) is given by, [s( 2t) erfc (L/Is 2)1 II i,(z)I \D(s)\ < -, (3.11) 2 yz s22i (z )dz where, erf (z)- $2 e-Y dy. The constant L determines where the change from io(.) to ie(.) occurs. Case la: i,(.) is piecewise linear, A for z>E ] i(z)= for I z <E. -A for z < E

62 Case Ib: io(.) is smooth, i0(z) =A erf (z/E) =A 2 fe-/E)2dy. In cases la and Ib, the constant E determines the width of the transition in io(.) from one region to the other, and the constant A determines the amplitude of the transition. The migration bound for case la is given by, [s( 2) erfc (L/s 2)]\ I i,(zl.1 \D(s)l < - L E (3.12) 2AE 2 e )-e2 The migration bound for case lb is given by, [s( 2 ) erfc (L/sx[2)] 11 i e(Z ID(s)l < 1 (3.13) JD) er ([(L/E)2 +(1/2) (L/s)2f ()2 -2A -- - erf Ye-) [1 + (1/2) (E/s)2] Case 2: The weighting function is given by -1 (zL)2 w(z)=e 2 so there is a smooth transition between io(.) and ie(.) outside the regional boundary, and (3.10) is given by,

63 ID (s)l < \ i \ / z - (z/s)2 - (z/L)2 2 J e 2 e 2 io(z)dz Case 2a: io(.) is the identical piecewise linear function used in case la. Case 2b: i(.) is the identical smooth function used in case lb. The migration bound for case 2a is given by, (3.14) ( [- + (/L)] 1 i(z - ID(s)l < I - - (3.15) 12 erf ( I {(EIs) + (E/L)1^}] [1 +(s/L)2 r[(E/s)2+ (EIL)2] The migration bound for case 2b is given by, s(\ 2 ) [1 - + (s/L)2] )1 I (Z ID(s)l < 2A [1 + (slL)2][l +1 (Es2 +(EIL)2 (3.16) In all cases, the expressions for bound migration can be written in the form, ID(s)l <K --- A A (3.17)

64 where the coefficient Kb is case-dependent, but is always a function of E/s and L/E. Values of Kb for the four cases are plotted as functions of E/s for L/E=1 in Figure 3.2a and forL/E=2 in Figure 3.2b. As expected, Figures 3.2a and 3.2b demonstrate that migration bounds decrease as scale decreases with respect to the transition width, E. In addition, migration bounds decrease as ie(.) is eliminated from the neighborhood of the boundary (i.e., as L/E increases). However, the shape of the weighting function, w(.), is clearly more influential on the migration bound than the shape of io(.). 3.2.2 Applications To compare boundary migration with the first-order theoretical bounds of (3.12), (3.13), (3.15) and (3.16), we simulate the migration for a simple example. Consider the case of the signal i(x) shown in Figure 3.3. In the figure, i(x) is seen to have a boundary (i.e., cross the zero threshold level) at two points, x=O and x=-1.6. The even and odd components (iO(.) and ie(.), respectively) for the boundary at x=O are calculated and shown in Figure 3.4. From Figure 3.4, it is seen that the signal boundary at x=O is an example of case la (~3.1), (3.12) can be used to calculate bounds for boundary migration where, A =i, (.)j =2 E=L=1.6. The resulting migration bounds are compared with simulated boundary migration in Figure 3.5.

65 Migration Bounds LE = 1 10 9 8 7 6 g 5 4 3 2 1 0 — Case la -— Caselb Case2a......... Case 2b 0 0.20.4 0.60.8 1 12 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 E/s (a) L/E=1. Migration Bounds UE=2 10 9 8 7 6 g 5 4 3 2 1 0 - Case la -- Case lb -Case2a......... Case2b - i I t - I I II i ti. - -- - 1 - I I 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 E/s (b) L/E=2. Figure 3.2. Migration coefficient value as a function of normalized edge width for two odd-signal widths. L is the length of the odd-signal, E is the edge transidtion width, and s is the scale value.

66 In Figure 3.5, boundary locations are shown as a function of filter scale; the abscissa is the boundary position and the ordinate is the scale value. The grey region on the plot represents possible boundary positions within the migration bounds calculated from (3.12). The solid curves are the boundary migrations determined from simulation. Figure 3.5 shows that for i(x) of Figure 3.3, the boundary at x=O migrates at its theoretical limit (i.e., maximum migration). Figure 3.5 also shows that despite the fact that the expressions for migration bounds in ~3.2.1 are first order estimates, they can hold over large changes in scale level. Signal and Threshold 4 Signal 3 ---------........Threshold 2 CD 0. c. - 1 CD) C -2 J. -.IL/ -3 4. | I t Is f i | l, | I -5 -4 -3 -2 -1 0 1 2 3 4 5 Position X Figure 3.3. Ideal signal and threshold.

67 Even and Odd Components Right Edge 4 Odd 3 -Even -2 -1- -------- -2 -3 -4 -5 -4 -3 -2 -1 0 1 2 3 4 5 Position X Figure 3.4. Even and odd components for right edge of ideal signal. Figure 3.5. Fingerprint and migration bound for right edge of ideal signal.

68 Figure 3.6. Even and odd components for left edge of ideal signal. 2.04 1.84 1.64 1.44 1.24 1.04 0.84 0.64 0.44 0.24 0.04:(::: -I T- - f,-I. -5.0 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 X Position 3.0 4.0 5.0 Figure 3.7. Fingerprint and migration bound for right and left edges of ideal signal.

69 Now consider the signal boundary at x=- 1.6. At this boundary location, the even and odd signal components are calculated and shown in Figure 3.6. From Figure 3.6, it is seen that the signal boundary at x=- 1.6 is not modelled precisely by any case considered in ~3.2.1. However, the boundary can be approximated by case 2a, so that (3.15) can be used to calculate bounds for boundary migration where, A =IIi(.l. =2 E=0 L=1.3. The resulting migration bounds are compared with simulated boundary migration in Figure 3.7. As in Figure 3.5, Figure 3.7 shows boundary locations as a function of filter scale, where the grey regions represent possible boundary positions within the migration bounds calculated from (3.12) and (3.15). The right-hand grey region in Figure 3.7 is determined from (3.12) for the boundary at x=O, and is identical to the grey region in Figure 3.5. The left-hand grey region in Figure 3.7 is determined from (3.15) for the boundary at x=-1.6. The solid curves are the boundary migration. As expected, Figure 3.7 shows that for i(x) of Figure 3.3, both boundaries migrate at their theoretical limits. Figure 3.7 again shows that the first-order, migration analysis of ~3.1 can generate accurate bounds over significant changes in scale level, even when using an approximate (non-exact) weighting function, w(x).

70 In addition to the example of Figure 3.3, the boundary migration analysis of ~3.1 is applicable to edge detection applications in computer vision. For example, an ideal "smooth" edge may be given by [75], i erf (ax)forx~ 0 - ) -erf (ax) for x< o' where the parameter "a" controls the sharpness of the edge. Using the Marr-Hildreth edge detector, edges in i(x) correspond to zero-crossings of i"(x) [52]. Using i(x) above, i' (x)= = (ax)e-<ax? so that zero-crossings of i"(x) are examples of case 2b boundaries and, for the ideal edge above, 11 ij(.j p = 0. Thus, the ideal edge produces no boundary migration. Another application of the migration analysis of ~3.1 is non-symmetric thresholding of boundary regions. Consider the signal and threshold shown in the top panel of Figure 3.8. Although the signal is strictly odd (about zero), the threshold-crossing of i(x) will migrate because the threshold is offset from zero. This is seen by determining the odd and even signal components with respect to the signal boundary (threshold-crossing) at x=l. These signal components are shown in the bottom panel of Figure 3.8, where the even component produces boundary migration. A migration bound can be calculated using the analysis of ~3.1. However, as opposed to the cases considered in ~3.1, the odd signal component for the signal of Figure 3.8 is not diminished outside a neighborhood of the boundary. Thus, the migration bounds given by the analysis of ~3.1 would be conservative

71 for the migration of the signal of Figure 3.8. However, it is a straightforward exercise to reformulate the expressions of ~3.1 to address cases where odd signal components are not limited to neighborhoods of the boundary. 3.2.3 Remarks In some edge detection analyses of computer vision [5,6,17,21,5 1], edge migration as a function of scale is calculated exactly by making the (strong) assumptions that all edges can be composed of ideal components, such as step and staircase edge, and that edge operators, 0. }, are first or second order derivatives. Alternatively, we attribute boundary migration to the non-ideal nature of edges. Moreover, migrations can be bounded by making weak assumptions, or by using easily obtained a priori information, about the sensors and symmetric properties of the underlying signal. For example, migration bounds are dependent on the peak level of the even signal component, 11 i,(.j)l, which is generally straightforward to estimate. In addition, migration bounds are calculated for signals that are odd-only in neighborhoods of boundaries. This is satisfied by signals that are approximately linear in a neighborhood of boundaries. Moreover, Figure 3.2 shows that bounds are relatively insensitive to modelling variations of the linearity. The weighting functions used in the migration analysis were considered to be "worst case", in that odd signal components, i(.), were excluded from regions far from boundaries. As seen in ~3.2 for non-symmetric thresholds, io(.) may extend to infinity, giving a boundary migration that is much smaller than the bounds predicted from the analysis. However even for non-symmetric thresholding, the weighting functions of ~3.1 will bound the migration. More accurate modelling of signals, in the vicinity of boundaries, is necessary only if tighter bounds are required. A more significant problem may be in the shape of the weighting function in the vicinity of regional boundaries.

72 -- Non-Symmetric Threshold 4 3 -J 2 ) 1 a) 0 c) -1 Cn -2 -3 3 2 1 0 CZ.O -1 cO -2 -3 -4 Signal Threshold Odd Even -5 -4 -3 -2 -1 0 1 2 3 4 5 Position X Figure 3.8. Ideal signal, and even and odd components, for non-symmetric thresholding.

73 It may be possible that a weighting function exists which gives migration bounds that are wildly different from those obtained from (3.17) using the curves of Figure 3.2. However, we assume such weighting functions are unusual for natural signals, and that there is a monotonic transition from odd-only signals, in the vicinity of a boundary, to a mixture of odd and even signals, away from a boundary. A "worst case" condition would have odd-only signals near a boundary, and even-only signals away from a boundary. Such assumptions are not unusually restrictive, since similar assumptions are the basis for regularization techniques in edge detection [72]. Under a monotonic transition assumption, the resulting weighting function should produce a migration bound that is similar to those determined by (3.17) using the curves of Figure 3.2. Finally, we note that the migration bounds of this section are strictly valid for 1-D signals. In the case of 2-D signals, this requires that signal characteristics in the neighborhood of boundaries can be described using 1-D transects, with minimal loss of information. This is often the case, as illustrated in chapter 5. In considering more general signals, a 2-D form of equation (3.6) is given by [36], -r(x0,y0,s) = rX(xOy0,s)(x -x) + ry(xO,yO,s)(y - yO), (3.18) where the point (xO, yo) is defined by the zero-crossing of r(xO, yo, 0), and the 2-D boundary migration is defined as, D(s) (x -X)2 +(y - y)2.

74 We see from (3.18) that D(s) can become unbounded if, for example, rX(XOYOs) -ry(xOyOs). Thus, the general 2-D migration must be calculated exactly on a case-by-case basis, which is the approach used in [5,6,17,21,51]. 3.3 Boundary Extrapolation We now consider a multiresolution system, derived from multispectral processing, where the objective is to accurately determine the boundary of a frozen surface3. In this application, frozen and wet surfaces can look identical to a sensor at a single frequency, and multiple (frequncy) sensors are used to distinguish frozen from wet surfaces. The top panel of Figure 3.9 shows the ideal (infinitesimal PSF) response, il(x), from a single sensor to a dry surface containing frozen and wet region. The wet region is centered about x=-138 Km, while the frozen region is centered about x=+138 Km. Both regions a approximatly 138 Km wide. From i,(x), it is impossible to determine which are frozen and wet regions. Alternatively, an indicator function, ii(x), is constructed from the difference between two sensors, i,(x) = i(x) - i2(x), (3.19) where() and.) are the i(.eal responses of the ideal responses of the response of ii(.) to these surfaces is shown in the bottom panel of Figure 3.9, where ii(x)<O indicates the presence of a frozen surface. Thus, iI(x) clearly shows the existence of the frozen surface. 3 This model is derived, and simplified, from results of chapter 4

75 Single-Sensor Output 250.0 "...... I...........%..e... I I I...........',,.,'...................... '. 247.0 [ V — -. -4 cn 244.0.............. I.......I....................... I....... ii Threshold I I.......I I I 1....... I....... 241.0 238.0 ^,% ^ I 235............................ -400.0-320.0-240.0-160.0 -80.0.L L.. L L LI...lll..... I.l......l. 0.0 80.0 160.0 240.0 320.0 400.0 x (km) Indicator Output,.,;;.~....;.,....,......................I......,.................,........ -...... 6.0 5.0 4.0 3.0 2.0 1.0 0.0 -1.0 -2.0 -3.0 i, Threshold.I.............. l..... l....... 4.0 c n% -3.U................. "". -400.0-320.0-240.0-160.0 -80.0....................................I......] ' 0.0 80.0 160.0 240.0 320.0 400.0 x (km) Figure 3.9. Ideal single-sensor and indicator signals, and thresholds.

76 The problem with the indicator is that the functions i1(.) and i2(.) are not observed directly. Instead, their corresponding spatial filtered signals rl(x,sl) and r2(x,s2) are used, where the sensor PSF's, hl(x,s,) and h2(x,s2), provide the spatial filtering. In this example, the sensor PSF's are Gaussian, 1 ~(x/s.)2 h,(x, i) = e 2 for i=1,2, with s1=30.0 Km and s2=100 Km. Since sensor 2 has coarser resolution than sensor 1, the resolution at which ii(x) can locate the boundaries of a frozen surface is limited to the (coarse) resolution of sensor 2. The error in estimating the boundary location at the resolution of sensor 2 can be derived from the migration analysis of ~3.2. In this case, (3.12) is used to determine the migration bound (estimation error) where, from Figure 3.9, A = i(.)l =3K E=0 L=138 As a result, the migration bound at s2=100 Km is, 100(2 1 -n) erfc (138/1002) Km D(s)\ < — -34.2 Km, n e so that the estimation error is almost a third of the width of the region.

77 The error in estimating the boundary location can be reduced by incorporating fine-resolution information of sensor 1. Such a technique for extrapolation-in-scale is described in the following. 3.3.1 Theoretical Result Let threshold crossing contours, x(s), be defined by, h(x(s),s)*i,(x(s)) = a(x(s)), (3.20) where ax(.) is a linear threshold function (chapter 2) and ii(.) is an indicator function composed of a linear combination of single sensor data, M i,(.)- ajij(.). j=1 The contours x(s) in (3.20) are estimates of boundary locations in ii(.), where actual boundaries occur at x values where ii(.) crosses the threshold function, i1(x) = a(x). According to migration analysis (~3.2), we must determine x(s) for as small a scale as possible, to minimize estimation errors in boundary localization. If Sj is the finest scale at which data from e channel is a theen x(s) can only be determined for s > sM without introducing resolution ambiguities. However, the boundary estimate can be improved if there exists a threshold function P(.) such that,

78 h(u(s),s)*i,(u(s)) = 5(u(s)), (3.21) where u(s) are the fingerprints derived from (3.21), and u(s) =x(s) for sm 2 s M> s. (3.22) That is, x(s) can be extrapolated-in-scale from S=SM to s=sl through the function u(s). Note that u(s), for a particular threshold function P(x), need only approximate x(s) over part of a single contour and that different threshold functions are used to approximate x(s) for different contours. The requirements for exact extrapolation-in-scale are given in the following theorem: 3-1 Theorem (1-D Extrapolation Theorem) Given contours x(s) and u(s) defined by (3.20) and (3.21) respectively, then u(s) is identical to x(s) in a neighborhood of the point x0 if, h (xo, so)* i1(Xo) dx A-h(xo,So)* i,(xo) h (xo, So)* - i (Xo) d n=2,3,4,... B-h(xo, So)*d^il(Xo) (A proof of the theorem is given in Appendix B). In multispectral processing, data from M sensors are integrated, where each sensor operates at a different (center) frequency and, as each sensor is scanned, each sensor output

79 signal is the convolution of the spatially varying radiation intensity of the surface and the PSF (antenna pattern) of the sensor. For the jh sensor, the radiation intensity of the surface is equivalent to the scale-space input, ij(.), and the PSF of the sensor is equivalent to a scale-space filter impulse response of the jh sensor, h(.,sj). The scale of the impulse response, Sj, is the width of the PSF (i.e., the beamwidth), and is proportional to the wavelength of the sensor. A standard product model of source imaging used in image processing and remote sensing is given by, ij(x) =I(x)Rj(x). In active systems, I(.) is surface illumination and Rj(.) is surface reflection. In passive systems, I(.) is surface temperature and Rj(.) is surface emissivity. In both systems, Rj(.) is a function of the surface type and is dependent on frequency. Using this model, the indicator is given by, i,(x)= XI(x)Rj(x), (3.23) A class of functions that satisfies the theorem are quadratic functions. Using the product model of (3.23), such signals occur where Ij(.) and Rj(.) vary linearly in the vicinity of a boundary, and we have, I(x) = acx +, Rj(x)= a=RX + iRj so that,

80 i(x)= a)x2+ (bj)+(S1) (3.24a) and i1(x) = ax2 + bx + c1 (3.24b) where, a, = a,a, bj = a]OR; + OlaRj cj = aIRj this reduces to a linear case if a, and aRj are much smaller than p1 and AR3, respectively. 3.3.2 Application We now return to the example considered at the beginning of ~3.3. Comparing the signals i,(x) and ij(x) in Figure 3.9, the conditions of theorem 3-1 hold in the vicinity of the frozen surface because i,(x) and i,(x) are translated and scaled versions of each other in the neighborhood of x=+138 Km. Thus, il(x) may be used to extrapolate i1(x) in scale. The extrapolation process can be seen from Figure 3.10. The bottom panel of Figure 3.10 shows the threshold crossing contours, x(s), of the indicator, ii(x), and associated threshold (the bottom panel of Figure 3.9); x(s) is shown as an arched curve. It is seen that at scales below approximately 60 km the position x(s) is at the actual boundary locations of the frozen surface. However, the minimum scale at which x(s) can be calculated is 100 km (i.e., the resolution of sensor 2), as indicated by the horizontal line in the bottom panel of Figure 3. 10. Thus, the estimated boundary of the frozen surface using strictly x(s) are the positions where x(s) crosses the horizontal line in the bottom of Figure 3.10. As expected, there are considerable differences between the estimated and actual boundary positions.

81 Single-Sensor Fingerprints E c) ut 120.0...................... 110.0 100.0 90.0 80.0 70.0 60.0 50.0 40.0 30.0 Scale Limit 20.0 10.0 0.0........................... -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 X (km) Indicator Fingerprints ) 240.0 320.0 400.0 Co a) 0 "jQ 120.0.... |.. 110.0 100.0..... Scale Limit 90.0 80.0 - 70.0 - 60.0 50.0 40.0 30.0 20.0 10.0 0.0............0.... -400.0-320.0-240.0-160.0 -80.0 0.0 80.0 160.0 240.0 320.0 400.0 x (km) Figure 3.10. Fingerprints and scale limits for ideal single-sensor and indicator signals.

82 The top panel of Figure 3.10 shows the threshold crossing contours, u(s), for the single sensor signal, il(x), and associated threshold (the top panel of Figure 3.9); u(s) is shown as a double-arched curve. While u(s) generates an extra (erroneous) contour as compared with x(s), the u(s) contour centered at x=+138 km registers with the x(s) contour. Additionally, the u(s) contours can be calculated to a minimum scale level of 30 km (i.e., the resolution of sensor 1), as indicated by the horizontal line in the top panel of Figure 3.10. Thus, the frozen surface boundary estimates of x(s) can be extrapolated from a scale level of 100 km to a scale level of 30 km using the contour u(s) centered at x=+138 km. Moreover, the extrapolated scale level is small enough so that the extrapolated estimate of boundary position is equal to the true surface boundary position. 3.3.3 Two-dimensional Extensions A 2-D version of the extrapolation theorem is formulated in the following theorem4: 3-2 Theorem (2-D Extrapolation Theorem) Let surfaces (h,s) and (x,t) in 3-space be defined by the function-crossings, r(x,s) = c() (3.25a) p (x, t) = P(_), (3.25b) where x - [xl xj2 T, and r(x,s) = h(x,s)**i,(x) (3.26a) p (x,t) = h ( t)**i() (3.26b) Cx() =Ao+Alx +A2x2+A 3xx2 (3.27a) 3(,) = Bo + Blx, + B22 + Baxx2. (3.27b) 4 In this section, the (x,y) plane is referred to as the (x1,x2) plane.

83 The coefficients Ao, Al, A2, A3, Bo, B1, B2, and B3 are constant. We define initial conditions by, r (o s~) = a(x~) P (o, to)= (xo The surface (s,s) is identical to the surface (f,t) in a neighborhood of the point x~if, a [r, SO)- a0)] oa [p ro to)- Po)] r,,,(x, s0~)+ r(X~,s~) Pxxl 0t))+px~, ) and, a~, r..%., atc axil... axiM) r1,X,,(X, s) ) -+ rp11~, st ) +p,(o,) o) where, M=2,3,... i,,i2,...,iM= 1,2 * (A proof of the theorem is given in Appendix C). To place the 2-D extrapolation theorem in a form similar to that of the 1-D theorem (~3.3.1), we use the following lemma: 3-3 Lemma

84 Conditions, [r (x ajD( SO) - j)][p A to) - PX! r11x1(Xj, s0) + rXO s 0) Px'x1 (X0, t0) + Px2r2 (-0 9 to) and 1 'M~~ax~ lp-o I M(Xl rx~1, (Xj0, S O) + rzr (X~s 0) Pxx(X0 to) + Pxx(X0 to) where, are satisfied if and only if, a [r (x, so) - ax.!)] ~ p (Xj03 to) - RXO)] ___ __ __a AI,,O2oO a 0 0) a x a W [r(~x,s) - (XO)] I. p (X- oO) - p(and, C ' I M I Mi _ [r9O s) - a(X(O)] I p (X-j0, S o) (- f(!)] (A proof of the lemma is given in Appendix D).

85 From (3.27a) and (3.27b), -a(xo)=A +A. (3.28a) axl a a(x )=A2+A x (3.28b) ax2 -a )=B +B3x0 (3.28c) ax, aO(x)=B,+BS3, (3.28d) aX2 so that, ax, x 32(3.29a) [r~,s~)-a(~)] = h(X~)**f i@l -At-A3x) (3.29a) a[r(~0s~)- ~)] = h@~,s)** xi,(x)] -A2-A3x) (3.29b) ] ax, (3.29c) a [P~~'t~)- [(~)] = h(-~,t~)**( xr il, ~ -Bn -BA~) x[p, (~ t~) - P~)] = h (~, t)**[ ax, -2-B3x (3.29d) By substitution of (3.29a)-(3.29d) into lemma 3-3, and using theorem 3-2, we obtain the following corollary: 3-4 Corollary

86 The surface (x,s) is identical to the surface (x,t) if five conditions hold: h(jso)**J Vixo] -A1 -A~x~j h~xo,so)** Ji(&)1 -A2- A34.x - hx -B1 -Bx) * h(x~o,tO)** i(X-O)] -B1 -B3xo) h(xO, sO)** Ji~jxO)1 - A) (3 ~,so)** J (i@1o) -Al-A,4= * h@xo,to)** ~i(O B h(xj0,tO)** j(xV) -B2 -B$xo (4)

87 h~,s~*'~] -A2-A3x~) h(O~, to)** j ) h(x~, t~)** (ilx ~)] -B2-Bx~) where, M =3,4,... JiJ2 *", jm= 1,2. Thus, if there exists signals i1(x) and il(x) which satisfy corollary 3-4, then the fingerprints of i1(O) match those of i1(x) in a neighborhood of the point x~. A 2-D formulation of the quadratic functions of ~3.3.1 illustrates an application of corollary 3-4. Using the standard product model of (3.23), we consider planar variations of I(.,.) and Rj(.,.) in the vicinity of a boundary. (Recall that surface reflection or emissivity, R(.,.), is indexed according to sensor.) For this case, I(x,y) = a, + Px +Y y Rj(x,y) = aRj + PRIX + RjY, and the ideal (infinitesimal-resolution) signal from the jm sensor is given by, ij(x,y) =I(x,y)Rj(x,y) = a + bx +cj y +djxy +ejx2+fj y2, (3.30) where,

88 a, alocRj bj atlfRj + Of3JORj Cj aYRj +7JlaRj dI3JYRj + 'YIf3Rj ej IPRj From (3.30), the indicator function is given by, i1(x, y)= XI(x, y)Rj(x, y) = a, +b~x + c.y + dxy +x+fy, (3.31) where, a,: Xaj ii ii ii de, =Xe. iJ e.=I and the signal used for extrapolation to fine-resolution is given by, il(x, y) = a, +b~x +cly +dxy +e +y2. (.2 (3.32)

89 Applying the conditions of corollary 3-4 to (3.31) and (3.32), the fingerprints of i1(x,y) will track those of ii(x,y) if, en el fe fi kf-" (3.33) Thus, when I(x,y) and Rj(x,y) vary along the same (x,y) direction, equation (3.33) and corollary 3-4 are satisfied.

CHAPTER 4 FREEZE/THAW CLASSIFICATION This chapter presents characteristics of the remotely sensed signals used for this thesis. We review the general expressions describing thermally emitted microwave radiation, and those for the specific case of frozen terrain. We present an indicator signal for classifying frozen terrain, and discuss decision criteria. Performance characteristics of the satellite radiometer used are also given. 4.1 Passive Microwave Radiation Thermally emitted (passive) microwave radiation from land surfaces can be described by radiative transfer theory [15]. The fundamental quantity in radiative transfer theory is the specific intensity. The power dP of non-coherent radiation flowing within a solid angle dQ2 through an area da in a frequency interval (f,f + df) is given by [73], dP =I(7r,)cosOdad2df (4.1) where, I(r,S) specific intensity, 7 is the point where intensity is determined, f is the direction of energy propagation at r, and 8 is the angle between & and the normal n of da (Figure 4.1). Power has the units of 90

91 watts and specific intensity has the units of watts m-2 Sr2 Hz'. Specific intensity is sometimes referred to as spectral brightness at microwave frequencies [74], and as spectral radiance at optical frequencies1. A radiation field is isotropic if the intensity is independent of direction s, and a field is homogeneous if the intensity is the same for all points r. Throughout this section, we will refer to specific intensity as intensity. Let a land surface be modelled as an isotropic dielectric, with complex dielectric constant ~ = E' - je", occupying a negative halfspace bounded at z=O (Figure 4.2) [25]. The surface is specular over the frequencies of concern and the temperature within the dielectric, T(z), is function of depth, -z. The dielectric contains scattering centers of cross section a, and density N, which neither absorb nor emit energy. The positive halfspace is occupied by free space. A n A S Fig. 4.1. Specific intensity. 1 Spectral radiance usually denotes watts m-2 Sr2 per unit wavelength, rather than unit frequency [82].

92 z/ z=O Free Space Dielectric ds -z Fig. 4.2. Radiative transfer model. By assuming rotational symmetry with respect to the z-axis, equations of transfer can be written for horizontally and vertically polarized intensities at wavelength AX as [25], dI(z, = I(z, ') + E,(z) + ds + I X fp ', ")Iq (z, ")dgj,z<0 (4.2) \= Hi - where,

93 a - extinction coefficient E - emissive spectral power co, scattering albedo PpqI scattering phase function, p is polarization (either V or H), 0' is the direction of energy propagation and g' is the direction cosine (g' = cos 0'), and ds is an interval along the direction of propagation. The extinction coefficient is the sum of the absorption and scattering coefficients, a 2p + Nca, and the absorption coefficient is given by, = 2 (4.3) The albedo is the ratio of scattering coefficient to extinction coefficient, Na, C,)~ 0 a Ep(z) represents the spectral power emitted within ds, while Pp q"(', p") relates the radiation over all directions p" to that scattered in the one direction J'. That is, equation (4.2) shows that the change in intensity along ds in decreased by the absorption and scattering of the material in ds, but increased by the self-emission of the material in ds and by the energy scattered into ds from all other directions.

94 The emitted intensity, at microwave frequencies, rom a blackbody at thermal temperature T is approximated by the Rayleigh-Jeans Law as, Ib(T) = CT (4.4) where C is a constant. To obtain intensity per unit wavelength, C is written as C, where [74] 2ck CX-= x- 4 ' (4.5a) and, k - Boltzmann's constant c - speed of light. To obtain intensity per unit frequency, C is written as Cf where [74] 2k Cf = For either form of C, the emissive spectral power is given by, For either form of C, the emissive spectral power is given by, (4.5b) E = (2p)( ~')CT. (4.6)

95 The radiation emitted into the positive (free space) halfspace at a viewing angle 0 is given by [25], Ip(+z,) =[ 1 -R (0')] (4.7) where 0 and ' are related by Snell's Law (Figure 4.2), andRp(g') are the Fresnel coefficients, ga' + 'i 1- ~'+ E~(Ct')22 i1'+ l- 1 '+c'(i')2 4 f1 - c' + c'(g')2 Ip(O, p1') is the solution to (4.2) at z=0, and is the intensity just below the dielectric/free space interface. By replacing the thermal temperature of the dielectric by it's first-order approximation, r(z) = To+ zl, (4.9) az=0 where, To surface temperature, the scatter-free (o0 = 0) solution to (4.2) as z t 0 is given by [25], (0,.t),o = o C[ + ( = J. (4.10)

96 Substituting (4.10) into (4.7) and noting, from (4.4), that I/(+z, 0)C is a temperature, the scatter-free brightness temperature is given by, TBP(O)o0 = [1 - Rp([)] To +2 z (4.11) That is, the intensity observed in the positive halfspace at viewing angle 0 is produced by substituting TBP(0)o, = into (4.4). In the case of isotropic scattering, Ppq4"( ', ") = 6pq" (4.12) where 8pq, is the Kronecker delta, Gaussian quadrature is used to approximate the integro-differential equation (4.2) with a system of first-order differential equations [ 15,25], dI (,r i' -=,p (tI )[C o (a)j (j t) l d~c ~ ( a, a -I = ~j 2 / = ~ (4.13) where,i = +1, 2,...,~n. x is refered to as optical depth and is given by X — az, a, are samples of the direction cosines p' such that 4,' = - g-', and the weighting coefficients aj are Christoffel numbers such that,

97 n ai =a aand I aj = 1. j=l The solution to (4.13) is [25], Ip ('CTg (:) (a-'jC To+ + { L e) L-cc ft ),,= O j = I l+rjlij (4.14) where rj are the n positive roots of the characteristic equation, 1 n ai j=C 1 — 2r(pi')2 (4.15) The coefficients Lpj are determined by applying (4.14) to the boundary condition, I (0,-[') = R? (C') (O, C ), (4.16) and solving (inverting) the resulting matrix equations. The resulting angle-sampled brightness temperature for isotropic scattering is given by, Trp(0i) ={ (') {T + ( zzLJ +j C(l+rigi)}' r,,(L) — -^^+lj ~zl~j i =I _^ + j^ (4.17) so that TBP(O) is interpolated from the samples TBP(9i). Comparing (4.11) and (4.17), the temperature difference due to scattering is given by, ATB, (O) = TB, (i) - TB, (RO), = O. (4.18)

98 This temperature difference, ATBp (O), is always negative and is referred to as scatter-induced emission darkening [25]. The intensity (per unit frequency) observed at a receiver in the positive (free space) halfspace due to the emissions from the negative (dielectric) halfspace can be expressed as, 2k IP(+z,e)=- TBP(Oi) (4.19) If the positive halfspace is an atmosphere and not free space, the received intensity is expressed as, Ip(+z,0)=- TAPp(9O), (4.20) where T^pp(O) is the apparent temperature at polarization p. Apparent temperature includes effects of atmospheric emission, atmospheric loss, and reflected atmospheric emission from the atmosphere/dielectric interface, as well as emission from the dielectric halfspace. To determine TApP(), equations of transfer are solved for the atmosphere halfspace, with the dielectric emission of (4.19) as a boundary conditions. Expressions for apparent temperature, assuming a scatter-free atmosphere, have been developed in [74] and are sumarized by, TAPp(0) = L( ) TBP() + Rp(0)TDNP(O) + TUPp(O;H)} (4.21) where H is the height of the receiver in the atmosphere, La(8) is the atmospheric loss factor (La(o) > 1), TDN(O) is the downward-emitted apparent temperture of the atmosphere, and Tup(0;H) upward-emitted apparent temperture of the atmosphere from z = 0 to z = H. To paraphrase (4.21), apparent temperature is increased by surface emissions (TB), upward

99 emissions of the atmosphere between the receiver and the surface (Tup) and reflected downward emissions from the entire atmosphere (R TDN). Apparent temperature is decreased by atmospheric absorption. From Snell's Law, the upward emissions of the atmosphere and the reflected downward emissions of the atmosphere are functions of the same (viewing) angle 0 (Figure 4.3). In the case of a lossless atmosphere, TDNP(0) = Tupp(;H)} = 0 La(0;H) = 1, so that, TAp (0) = TB(0). TDN Fig. 4.3. Apparent temperature.

100 4.2 Emissions from Frozen Soil Soil moisture contributes to the energy exchange between the air and the ground through latent heats of fusion and vaporization. The processes of thawing frozen ground or of evaporating soil moisture cause soil thermal inertias to appear anomalously high. There is a large body of literature about deriving soil moisture from radiobrightness [e.g., 4.6,4.7,4.8,4.9,4.10,4.11,4.12,4.13]. In addition, moisture state can also be inferred from radiobrightness. Using data from the Nimbus-7 scanning multichannel microwave radiometer (SMMR), a combination of low 37 GHz radiobrightness, TB(37), and a low spectral gradient of radiobrightness, aTB(f/af where f is frequency, becomes an effective freeze indicator, or discriminant, for classifying frozen terrain. The Nimbus-7 satellite, launched on October 24,1978, flew a 955 Km, sun-synchonous polar orbit and had local noon (ascending) and midnight (descending) equator crossings with 26.1~ of longitude separation. Its orbital period was 104.16 minutes [55]. SMMR, one of the components of the Nimbus-7, was a 10 channel, dual polarization (H and V) radiometer with a forward-looking, conical antenna scan of~ 25~, relative to the direction of flight. The incidence angle was constant on the surface of the earth at approximately 50~. Performance characteristsics for the five microwave wavelengths are shown in Table 4.1. Table 4.1. SMMR performance characteristics. I l Channel Parameter 1 2 3 4 5 Wavelength (cm) 4.54 2.8 1.66 1.36 0.81 Frequency (GHz) 6.6 10.69 18.0 21.0 37.0 RF Bandwidth (MHz) 250 250 250 250 250 Dynamic Range (K) 10-330 10-330 10-330 10-330 10-330 Absolute Accuracy (K rms) <2.0 <2.0 <2.0 <2.0 <2.0 RMS Temp. Resolution (K) 0.9 0.9 1.2 1.5 1.5 Beamwidth (Deg) ~ 2~ 4.2 2.6 1.6 1.4 0.8 Beam Efficiency (%) 87.0 87.0 87.0 87.0 87.0 + 25~ Scan Cycle (sec) 4.096 4.096 4.096 4.096 4.096

101 The SMMR radiobrightness data were derived from CELL tapes, obtained from the National Space Science Data Center. The SMMR data on these tapes were averaged and re-sampled [56], so that cell data dimensions (i.e., data sample intervals) were 30 km x 30 km for 37 GHz data, 60 km x 60 km for 21 GHz and 18 GHz data, 97.5 km x 97.5 km for 10.7 GHz data, and 156 km x 156 km for 6.6 GHz data. A single image from CELL tape data covers a 780 Km x 780 Km area on the ground. The radiobrightness data used at each frequency were averages of H and V channel data. The 21 GHz radiobrightness data were not used due to equipment problems on SMMR, and the 6.6 GHz data were not used due to the large cell dimensions. I 0-Willistoi D Miles City Newell 6 Rapid City Test Area Cavalier n D Thief River Falls Bismarck Fargo O-Eureka Waubay Aberdeen-O b 'Morris a Redfield Milesville Highmore H d db 0Huron -- a Brookings Cottonwood Fig. 4.4. Test area with air and soil temperature data stations. Air and soil stations are indicated by circles and squares, respectively.

102 Using these data, a surface is classified as frozen only if both TB(37) and aTB(/)af are sufficiently low, where aTB(f/)f is a linear, least-square fit to the 10.7, 18, and 37 GIHz radiobrightnesses. Frozen surfaces appear cold at 37 GHz, and exhibit a negative spectral gradient that is largely caused by volume scatter darkening at the shorter wavelengths. This two parameter freeze indicator (FI) has been applied to SMMR data from a test area that includes North Dakota, about half of the surrounding states, and part of Canada. This greater North Dakota test area is shown in Figure 4.4 with air and soil temperature data stations. Data was collected from August 1, 1984 through December 31, 1984. 4.2.1 Theory By modeling soil as a scatter-free, homogeneous halfspace, the radiobrightness temperature at a constant viewing angle, 0, can be approximated as a function of frequency f from (4.11) as, TB(f) = e (f)To + 8TB(f) (4.22) where, e () 1 - R (') at frequency f 6TB(f) e (f)[ \ 'z(f) zc(f) — at frequency f. e(f) is referred to as emissivity and z,(f) is the effective emitting depth of the soil, where (1-e1) of the total emission originates above z,(f) [30]. When compared to thawed surfaces,

103 frozen surfaces exhibit characterisitcs of (1) lower thermal temperatures, To, (2) higher emissivity, e (), (3) larger emitting depths, z,(f), and (4) decreasing radiobrightness, TB(f), with frequency caused principally by volume scattering. Characteristics (1) and (2) are well understood, but are generally ambiguous indicators of frozen surfaces. Ambiguities arise because changes in radiobrightness that result from freezing the surface may be either positive or negative, depending upon the surface moisture content. For example, a dry soil emissivity of 0.9 will yield a 9~ decrease in radiobrightness for a decrease in soil temperatures from +5~ C to -5~ C. Because the soil is dry, there is relatively little change in soil emissivity with freezing. In moist soils, freezing causes an increase in soil emissivity because water molecules in frozen plants and soils are not free to align themselves with microwave electric fields. This constraint upon the rotational freedom of water causes a decrease in the real part of the dielectric constant, E'(f). A moist soil emissivity would increase from 0.8 to 0.9 with freezing, so that a decrease in soil temperatures from +5~ C to -5~ C would produce a 19~ increase in TB(f). Because TB(f) can either increase or decrease with freezing, misclassification will results if TB(f) at a single frequency were solely used to discriminate between frozen soils and soils that are warmer or drier. These variations in emissivity with freezing are most pronounced at lower microwave frequencies. Signature (3) arises because freezing reduces the imaginary part of the dielectric constant, e"(f), proportionally more than it does the real part, E'(J). Reduced e"(f) means reduced absorption, so that thermally emitted photons originate deeper within emitting media. Thus, the effective depth of emission, z,(/), becomes a larger fraction of the free-space wavelength, Xo [26-29]. The effective emission depth of moist soils is typically 10% of the free-space wavelength. Frozen soils have effective emission depths that may be 30% or more of free-space wavelengths. As a result, the subsurface temperature gradient of frozen soil contributes more to radiobrightness than does an equivalent gradient in thawed soil.

104 The contribution can be several degrees at the lower microwave frequencies. Table 4.2 shows the contribution of thermal gradient on radiobrightness, aT8(f), and on spectral gradient, 6(aTB/af), as a function of soil moisture content and time of day (noon vs midnight). The data of Table 4.2 are derived from a model of typical soil near Bismark, North Dakota for September 22 [30]. The model used an incidence angle of 53.1~ and radiometer frequencies of 10.7 GHz, 19.35 GHz, and 37.0 GHz (horizontal and vertical polarizations have been averaged). 5TB(f) is calculated at these frequencies, and 6(aTBlaf) is the least-square regression slope to the three aTB(f) values. In the model, the soil contains 7% bound water, so that a soil moisture content of 10% has a mixing ratio of 0.03 free water. We see that thermal gradients exert the strongest influence on emissions from frozen or dry soil. The thermal gradient produced noon-to-midnight shift in spectral gradient, 6(aTBaf)NM, is 0.12 K/GHz for frozen soil and, as will be shown, is consistent with SMMR observations. Table 4.2. Contributions of thermal gradient to radiobrightness and spectral gradients for noon and midnight data. Radiobrightnesses, 56TB(), are in K and spectral gradients, 8(aTB/af), are in K/GHz; 6(aTB/af)NM is 5(aTB/af) at noon minus 6(aTB/laf) at midnight. 5TB(f), NOON 5TB(f), MIDNIGHT moisture 10.7 19.35 37.0 10.7 19.35 37.0.aor content, % GHz GHz GHz S) GHz GHz GHz e L Aar,) 7 (Frozen) -3.5 -1.8 -1.0 0.091 1.3 0.7 0.4 -0.028 0.12 10 -2.3 -1.0 -0.5 0.061 0.6 0.4 0.3 -0.013 0.07 15 -1.4 -0.6 -0.3 0.039 0.3 0.2 0.1 -0.007 0.05 20 -1.0 -0.4 -0.2 0.029 0.2 0.1 0.1 -0.005 0.03 25 -0.8 -0.3 -0.1 0.023 0.2 0.1 0.1 -0.003 0.03 As a consequence of signatures (1), (2), and (3), the 37 GHz SMMR radiobrightness is more strongly correlated with air temperature than are the 10.7 GHz and 18 GHz SMMR radiobrightnesses. That is, the 37 GHz radiobrightness should serve effectively as one

105 discriminant among frozen and thawed soils. Through diurnally heating, the spectral gradients ofradiobrightness produced by (4.22) are always positive for thawed soils, and are slightly negative (typically -0.1 K/GHz for frozen soils) [30]. However, observed spectral gradients in frozen soils -- signature (4) -- may be more negative than - 1.0 K/GHz. The likely cause of such strongly negative gradients is emission darkening at shorter wavelengths, caused by volume scattering within the frozen soil. That is, (4.22) does not adequately model frozen soil. In relatively transparent emitting media such as frozen soil, ze is large and the greater average thermal photon path length provides a greater opportunity for volume scattering of photons. Scattering is more severe at shorter wavelengths because soils and plants appear increasingly heterogeneous at the scales of these wavelengths [24,27]. In moist soils, Ze is smaller and the volume scattering is less pronounced. Thus, a negative spectral gradient should correlate with frozen soil, and the radiobrightness spectral gradient should serve as a second discriminant among frozen and thawed soils. 4.2.2 Freeze Indicator The data in Figures 4.5-4.7 are derived from SMMR radiobrightness measurements made from August 1984 to December 1984 over seven air temperature (meteorological) data sites -- Miles City, MT; Bismarck, Fargo, and Williston, ND; and Aberdeen, Huron, and Rapid City, SD -- within the greater North Dakota test site. The radiobrightness data are averages of horizontally and vertically polarized SMMR radiobrightness measurements, and are spatially interpolated to the latitude and longitude of each meteorological site using a bi-cubic approximation to a sinc function2 [54]. 2 Data have been compensated to the (coarse) resolution of the 10.7 GHz channel.

106 Scatter diagrams of TB(37) as a function of meteorological air temperature (Figures 4.5a and 4.5b) show a nominal tracking of air temperature by TB(37). However, there is an approximate 4~ decrease between the noon and the midnight regression lines caused by air temperatures that lag surface temperatures. A simple regression model for TB(37) is, TB(37) = e (37) [TaR + TBAs(t)] (4.23) where, e (37) - 37 GHz Emissivity TAIR Air Temperature (K) TBAs(.) Temperature Bias (K) t - Time, so that, (12) T () = 4.35, (424) TBIAS(12) - Tins(0) = e (7 4.35K, (4.24) B'S e(37) where t=O for midnight and t=12 for noon.

107 37 GHz Radiobrightness Noon Data 19 U'd 0 1 0 A.0'A L-a'.0A 500 290 280,, 270.- -, --- 260.250- -— 5aLi - 240 - i I-D --- - - 230 --- 220, i n I i. I - dc IV i -- - () I I I I I --- I I 2 n I.... I I I. -1 I I I I I I I I I I I. I I I I I. I I. I I I I I I t.i I, 240 250.260. 260 270 280 290 Air Temperature (K) 300 310 320 o Data - Regression (a) Noon data.-r% Y. I.(A a4 37 GHz Radiobrightness Midnight Data 300 290 280 -270 240 Ca 230:220 - 210: 200,.,.T...,....... -.....,'. -... 240 250 260 270 280 290 Air Temperature (K) 300 310 320 Data -- Regression J (b) Midnight data Fig. 4.5. TB(37) versus measured air temperature at meteorological sites in North Dakota and the surrounding region. Data were collected from 8/1/84 to 12/31/84.

108 a 1. 1.4 1.2 1 O.E 0 -I ' 0.6 (g 0.4 Y 0.2 -oc ''.'~ -0.4 a. o —.E -1 -1. -1.I Spectral Gradient Noon Data 3 I. -- ===== or --- - -- - -- - -- - [ ~ --- - --- - -- - 1 -- ---------- I 240 250 260 270 280 290 300 310 320 Air Temperature (K) (a) Noon data Spectral Gradient Midnight Data 1. 1. 1. 0. a O. NO I '- 0. -0. ' -0. o -1. -1.,.1 c - 4 2 1 8 6 4 2: 0 1 M - -I- -- 1-1-1-1-1 5~'-8:i --- -- -- -- -- -- -- --,2 o -- _m __ —__4 o C 6 -1 2 -A7 -1,.6 -1.6 2 240 I r......... >..... 250 260 270 280 290 300 310 320 Air Temperature (K) (b) Midnight data Fig. 4.6. aTBlaf versus measured air temperature at meteorological sites in North Dakota and the surrounding region. Data were collected from 8/1/84 to 12/31/84.

109 Decision Space Noon Data 1.6 1.4 -1.2 '0.6 0 6 0.4 SIR_ >0.2 __ 0:1 ~E-0.2 -o-0.6 c5-0.8 __ _ -12 -1.2 -1.4 200 210 220 230 240 250 260 270 280 290 300 37 GHz Radiobrightness (K) C3 Thawed X Mixed A Frozen (a) Noon data Decision Space Midnight Data 1.4 1.: "I 0.i 0 c -0.: 0 -0.4 o -0.1 e5- -0.1 -1.1 -1. -1 - 6 4 -2 B B 4 2 4 * I 6 B-2 _ 2 --- —-...-.- I.... 2006 210 220 230 240 250 260 270 280 290 300 37 GHz Radiobrightness (K) 0 Thawed x Mixed A Frozen (b) Midnight dataI Fig. 4.7. oiT8/af versus TB(37) at meteorological sites in North Dakota and the surrounding region. Data were collected from 8/1/84 to 12/31/84.

110 Equivalent scatter diagrams of aTB/af as a function of air temperature are shown in Figures 4.6a and 4.6b for noon and midnight SMMR data, respectively. The values of aTB/If are the slopes of linear-least-square regressions, as functions of frequency, to SMMR 10.7, 18, and 37 GHz radiobrightnesses at each meteorological site. Specifically, at each meteorological site we calculate, aTB COv[T (f),fl V= r (4.25) Of Var [] where Cov[TB(f),f is an estimate of the covariance of TB(f) andf, and Var f] is an estimate of the variance off. We use the estimates [70], Cov[T8(f),f] = {(TB(37)- T) (37 -f)+ (TB(18)- TB) (18-f)+ + (TB(10.7) -TB) (10.7 -)} (4.26a) I _h) +2 VarU] =5{(37-f)2+(18- )2+(10.7-J)2} (4.26b) where f and TB are estimated mean values, f - {37 + 18 + 10.7} = 21.9 (GHz) -- 1 TB- {TB(37)+TB(18)+TB(10.7)} (K), and TB(37), TB(18), and TB(10.7) are radiobrightness values at 37 GHz, 18 GHz, and 10.7 GHz, respectively. As a result, we obtain,

111 aTB j{15.1TB(37)-3.9TB(18)- 11.2TB(10.7)} f 1{368.66} 0.041 TB(37) - 0.0106 TB(18)- 0.0304 TB(10.7) (K/GHz). (4.27) The data of Figures 4.6a and 4.6b show the predicted decrease in spectral gradient with decreasing air temperature (i.e., as surfaces freeze). There are also the anticipated 0.1 K/GHz increase in noon gradients relative to midnight gradients caused by diurnal heating and cooling. In addition, a negative tendency in aTB/af is observed at noon for high air temperatures. This may be caused by volume scattering by dried surface vegetation. Heat and plant senescence in late-summer decrease the moisture in the vegetation canopy. This dry vegetation will act as a scattering layer -- particularly at higher microwave frequencies. These data (Figures 4.5 and 4.6) yield scatter diagrams of aTBI/f as a function of TB(37) (Figures 4.7a and 4.7b). In Figures 4.7a and 4.7b, data labelled as "frozen" have air temperatures less than 270 K, and "thawed" data have air temperatures greater than 274 K. Data labelled as "mixed" have air temperatures between 270 K and 274 K, inclusive. While air temperature is an imperfect indicator of frozen terrain, we see that low TB(37) and aTB/af correspond to frozen surfaces. As a result, a two-parameter (two-dimensional) binary decision criteria is used to indicate the presence of frozen terrain. This freeze indicator classifies a surface as frozen if TB(37) is below a specified threshold P37, and if aTB/lf is below a threshold PSG. We note that this freeze indicator is not strictly a linear combination of signals, as the indicator functions discussed in chapter 3. Nevertheless, we see in chapter 5 that the freeze indicator can be extrapolated over scale, despite the non-linear relation between the indicator and the signal used for extrapolation (i.e, the 37 GHz radiobrightness).

112 4.3 Threshold Determination Clustering determines decision criteria for coarse-resolution classification, and the location of freeze/thaw boundaries. The location of these freeze/thaw boundaries can be refined by imposing spatial constraints. 4.3.1 Data Clustering The decision criteria for TB(37) and aT/laf were based on clustering and unsupervised classification. Unsupervised classification, rather than supervised classification, is used because of the dearth of accurate ground measurements in our test area. The eighteen air and soil temperature recording sites (Figure 4.4) provided the ground data for the entire test area. Soil temperatures were measured at 5 cm depth and were made at dawn and dusk, whereas SMMR overflights were at noon and midnight. To increase the number of data available for clustering beyond that used for Figures 4.5-4.7, all data from SMMR satellite passes that covered more than 67% of the test area were incorporated in the scatter diagrams of Figures 4.8-4.9. Sixteen noon SMMR passes and thirteen midnight passes met this criterion during our August to December test period. As before, the 18 and 37 GHz averages of vertical and horizontal polarized radiobrightness data were resolution compensated to the (coarse) resolution of the 10.7 GHz channel, but unlike the data of Figures 4.5-4.7, these compensated data were re-sampled on a 97.5 Km grid (i.e., at the resolution of the 10.7 GHz channel). Spectral gradients were computed as the slopes of least-square linear regressions to the 10.7, 18, and 37 GHz radiobrightnesses. Scatter diagrams for the noon and midnight data are shown by month in Figures 4.8a and 4.8b, respectively.

113 SMMR Cluster Data, Noon 1.600 1.200 I I I.* I.. I I I I I. I I.. N I C Q. 0 u, 0.800 0.400 0.000 -0.400 -0.800 A i A A * a ~ ~ - - Y J4A V o August o September A October v November o December -1.200 -1.600. l............ 1 200.0 210.0 220.0 230.0 240.0 250.0 260.0 270.0 280.0 290.0 300.0 37 GHz Brightness (K) (a) Noon data SMMR Cluster Data, Midnight N I 0 12 C c. 4) (I2 O 1.600 1.200 0.800 0.400 0.000 -0.400 -0.800 -1.200 A A A o Ate * October v November o December............... I...I.......I...... I........ -1.bUU% 1 -1. ' ' ' ' ' '................................................. 200.0 210.0 220.0 230.0 240.0 250.0 260.0 270.0 280.0 290.0 300.0 37 GHz Brightness (K) (b) Midnight data Fig. 4.8. Scatter diagram of aTBlaf versus TB(37) throughout North Dakota and the surrounding region. Data were collected from 8/1/84 to 12/31/84.

114 Migrating means clustering determined cluster centroids for the data of Figures 4.8a and4.8b [58]3. On the basis of ground measurements and the data of Figures 4.5-4.7, surfaces were classified into three distinct types -- frozen, hot (and dry), and wet (and cool) -- and a fourth type that we call mixed. A frozen surface is characterized by relatively low spectral gradient and a low 37 GHz radiobrightness. Due to the influence of liquid water in the surface, a wet surface is characterized by a high spectral gradient [30] and low 37 GlHz radiobrightness [39]. A hot surface has relatively less surface moisture, producing a "dry" surface dielectric constant similar to a frozen surface. Moreover, the relaxation frequency of free-water increases with temperature, further reducing the spectral gradients of hot surfaces. As a result, a hot surface has a relatively low spectral gradient and a high 37 GIHz radiobrightness. A mixed surface has a combination of frozen, wet, and hot characteristics. Prior to freezing, a surface region is a combination of wet (and cool) and hot (and dry). As freezing begins, the region includes locally frozen surfaces, and would be classified as mixed. A freeze/thaw criteria lies within the mixed surface cluster in decision space whose components, P37 and PSG (section 4.2.2), represents maximum TB(37) and aTBlIf values along the freeze/thaw boundary. That is, any surface point on the freeze/thaw boundary has either TB(37) equal to P37 or 7TBl/f equal to PSG. Equivalently, the freeze indicator (FI) algorithm requires any surface point classified as frozen to have TB P37 and aTB/af< PSG Cluster centroids determined for the data of Figures 4.8a and 4.8b are given in Table 4.3. Due to SMMR recording problems, limited midnight data were available during December of our test period. Within this limitation, the frozen surface cluster centroid has a lower spectral gradient and 37 GHz radiobrightness at noon than at midnight. Furthermore, because there were few wet surfaces at midnight during the test period, the wet and mixed surface types were not separable for the midnight data. 3 Clustering and classification was performed on a Sun-4 workstation using EASI software, version 4.1, from PCI, Inc. of Richmond Hill, Ontario (Canada).

115 Table 4.3. Cluster centroid in decision space NOON DATA MIDNIGHT DATA 37 GHz Gradient 37 GHz Gradient Surface Type (K) (K/GHz) (K) (K/GHz) Frozen 227 -.43 234 -.35 Hot 277 0.11 258 -.01 Wet 238 0.37 Mixed 250 0.14 243 0.015 Bivariate normal distributions were fit to the cluster data. All data within three standard deviations of a cluster centroid were classified using a Mahanalobis minimum distance classifier (maximum likelihood classification). No preferential weightings of surface types were used. Constant-deviation, single-class ellipses were drawn in decision space for frozen, hot, and wet surfaces (at noon) and for frozen and hot surfaces (at midnight) using the classified data. The freeze/thaw criteria were determined by allowing the deviation of all ellipses to expand equally until all ellipses intersected. Thus, the freeze/thaw criteria represent a point in decision space that is equally likely to be classified as frozen, wet, or hot. The resulting classification ellipses for noon and midnight SMMR data are shown in Figures 4.9a and 4.9b, respectively. The corresponding freeze/thaw criteria in decision space are shown in Table 4.4. Table 4.4. Freeze/thaw criteria in decision space; a are standard deviations of data within the ellipses. Refined Deviation at P_ 37 (K) P37 (K) PSG (K/GHz) Intersection Noon 252 249 0.0625 3.1 l Midnight 247 244 -0.044 2.55 c

116 SMMR Cluster Data, Noon 1.600............,.......... 1.200 ^ 0.800 -, * ^ --, 0.800 - c -0.400 -1.200... -1. ~. '.. ' /.., ~. -1.600,...,..,,.., 200.0 210.0 220.0 230.0 240.0 250.0 260.0 270.0 280.0 290.0 300.0 37 GHz Brightness (K) 1.600 —....... -; (s -0.800 -1.200 -1.600 I I I 200.0 210.0 220.0 230.0 240.0 250.0 260.0 270.0 280.0 290.0 300.0 37 GHz Brightness (K) (a) Midnight data Fig. 4.9. Single class ellipses of lust ver Data, Midnight 1the surrounding region. Data were collected from 8 /84 to 12/31/84. 1.200 0.800 v2 0.400 -0.400 - 0 -0.800 -1.200 -1.600 I. I. i 200.0 210.0 220.0 230.0 240.0 250.0 260.0 270.0 280.0 290.0 300.0 37 GHz Brightness (K) (b) Midnight data the surrounding region. Data were collected from 8/1/84 to 12/31/84.

117 4.3.2 Refining the Decision Criteria If we view the freeze/thaw criteria derived from clustering as initial estimates for determining the freeze/thaw boundary, the boundary can be refined by requiring a minimum scatter of TB(37) along that boundary. This constraint ensures that boundaries in 37 GHz radiobrightness images correspond closely to freeze/thaw boundaries. The TB(37) boundary pixels are those whose 37 GHz radiobrightness equals P37 of the freeze/thaw criteria. The scatter of TB(37) along the freeze/thaw boundary is minimized by changing the TB(37) component of the freeze/thaw criteria to give a minimum sum of squared error. Using the freeze/thaw criteria derived from clustering, a square error of 37 GHz radiobrightness is calculated at each pixel along the freeze/thaw boundary, SE5 = [TBs(37)-P37], where, SE - Square Error at the id Boundary Pixel TBi(37) 37 GHz Radiobrightness at the ith Boundary Pixel. The sum of squared error, SSE, is calculated by, N N SSE = X SSi = X [Tsi(37)- P3] P37 i T(37)(4.25) P 1=i1

118 The process is first-order since we do not reiterate SSE minimization with the refined criteria. Applying this process to midnight data from October 24 and to noon data from December 11 yields the refined P37 values shown in Table 4.4.

CHAPTER 5 FREEZE/THAW BOUNDARY PROCESSING This chapter presents a boundary localization procedure derived from the extrapolation-in-scale results of chapter 3. The indicator -- the freeze indicator -- and signal used for extrapolation -- 37 GHz radiobrightness -- are characterized in chapter 4. 5.1 Boundary Localization Procedure Spectral gradients are regression slopes to SMMR 37 GHz, 18 GHz, and 10.7 GHz radiobrightness measurements, as given by (4.27). The nominal resolutions of these data are 30 km, 60 km, and 97.5 km, respectively (Chapter 4). Without compensating for the resolution differences between the channels, the spectral gradient estimates can be in error. Forexample, a non-zero gradient estimate can result from radiobrightnesses that are spatially variant but are locally constant over frequency. To avoid such errors, the image data were compensated to one common resolution -- the (coarse) resolution of the lowest frequency channel used in gradient estimation -- prior to clustering. Freeze/thaw boundaries combine 37 GHz threshold crossings and spectral gradient threshold crossings. Corresponding 37 GHz threshold crossings occur in fine-resolution 37 GHz images, but not all 37 GHz threshold crossings represent freeze/thaw boundaries. Some, as previously noted, are boundaries between moist and dry terrain. Boundary localization is a three-step process that identifies pixels in fine-resolution, 37 GHz images that correspond to freeze/thaw boundaries at coarse-resolution. 119

120 Step 1: Uncompensated 10.7 GHz, 18 GHz, and 37 GHz SMMR data are compensated to the resolution of the 10.7 GHz channel using Gaussian filtering. Gaussian spatial filtering is appropriate for resolution compensation if the Fourier transform of the SMMR beampattern is (approximately) Gaussian [10]. Performance data for the Nimbus-7 SMMR antenna are limited. Some investigators have assumed an antenna pattern based upon a uniformly illuminated SMMR aperture [18]. The half-power beamwidth for a uniformly illuminated circular aperture is approximately X/D, where D is the diameter of the aperture and X is the signal wavelength [74]. The aperture diameter for SMMR is 80 cm [55], so that the half-power beamwidth at 37 GHz (X = 0.81 cm) for a uniformly illuminated aperture would be 10.125 mradians, or 0.58~. Thus, uniform illumination produces beamwidths that are much narrower than those of the SMMR. Alternatively, we assume that the Seasat SMMR beampattern [56] approximates the Nimbus-7 SMMR beampattern. Figure 5.1 shows samples of the Seasat SMMR beampattem and a Gaussian beampattern at a transmit frequency of 6.6 GHz. Figure 5.1a shows normalized SMMR and Gaussian beampatterns as functions of (D/l) sin(0), where e is the offset angle from boresite. Figure 5.lb shows normalized beampatterns as functions of normalized ground position for a nadir-boresited antenna at altitude H (the point X = 0 is directly below the antenna). The Gaussian model uses a circularly symmetric illumination function. Figure 5.la shows that the Gaussian model matches Seasat antenna data down to approximately -25 dB of the peak level. Equivalently, the linear response of the beampattern (Figure 5.b) suggests that the (spatial) filtering of radiobrightness by the beampattern can be adequately modelled as Gaussian. Errors resulting from modelling a non-Gaussian filter by a Gaussian cause the channels used for generating the indictor signal and associated thresholds, P37 and PSG. to have non-indentical spatial filtering. The resulting errors in determining threshold criteria are mitigated by refining the thresholding process (~4.3.2). The alternative process is to use

121 exact, non-Gaussian filtering in resolution compensation, and threshold-crossing hysteresis in tracking boundary locations from coarse-to-fine resolution (step 3). The former process was used due to the antenna data limitations of the Nimbus-7 SMMR, as well as ease of implementation. The Gaussian filters used to synthesize compensated data at resolution s2 from uncompensated data at resolution s, are, H(f,S) = e (5.1) where the filter width, S, is, S=(S2- s1 ) s, > s, andfis spatial frequency. Values ofSfor different configurations of resolution compensation are shown in Table 5.1. Table 5.1. Filter bandwidths for resolution compensation. Nominal Synthesized Filter Resolution, s, Resolution, s2 Bandwidth, S 30 km (Fine) 60 km (Medium) 51.96 km 30 km (Fine) 97.5 km (Coarse) 93.77 km 60 km (Medium) 97.5 km (Coarse) 78.885 km

122 0-= -5:: m -20= 'a- 3 I * -35 - j -405-1 m -4 N -50 -55:= -60: -ic Beampatter ns Power (dB) If% P% ).uu — 6.u00 -2-00 2.00 b.00 1 U.UU (DIX) sin 0 -4f- MZf~~Gassian (a) Beampatterns as functions of angle, with dB scale. Beampatterns Power (Linear) I I.d a E p 0.3 0.7 0.1 03 SMMR Gausmin -01 0.00i I I I I -020 -0.15 -0.06 0.05 NormiaizW Poeiian (X44) 0.10 0.20 0.15 (b) Beampattern as a function of normalized position on the ground, with linear scale. Figure 5.1. Normalized Seasat SMMR and Gaussian beampatterns.

123 Step 2: Freeze/thaw boundaries are identified in coarse-resolution, 37GHz radiobrightness images. Using resolution compensated data, TB(37) and aTB/af are calculated for each image pixel at coarse-resolution. Boundaries in coarse-resolution, 37 GHz images are identified where 37 GHz data equal the TB(37) freeze/thaw criteria P37 (Table 4.4). Pixels along these 37 GHz image boundaries with 7TB/af at or below that of the freeze/thaw criteria PsG are identified as freeze/thaw boundary pixels. Step 3: Freeze/thaw boundaries are identified infine-resolution, 37 GHz radiobrightness images. Fine-resolution, freeze/thaw boundaries are determined by identifying those pixels in fine-resolution, 37 GHz data that equal P37 and correspond to coarse-resolution freeze/thaw pixels of step 2. This process involves tracking boundary locations in 37 GHz images as the amount of resolution compensation is reduced. The resulting boundary locations in the fine-resolution 37 GHz images are best estimates of freeze/thaw boundaries, in the sense that they are directly traceable to the coarse-resolution boundaries generated by clustering and maximum likelihood classification. The key to this process is that Gaussian image degradation of step 1 uniquely permits recovery of some fine-resolution information, as discussed in chapter 2. 5.2 SMMR Images Unrefined freze/thaw criteria (Table 4.4) were applied to SMMR data for midnight October24 (Figure 5.2) andfornoon December 11 (Figure 5.3). Refined freeze/thaw criteria were also applied to the October and December data (Figures 5.4 and 5.5, respectively). The dark pixels in thefreeze maps (Figures 5.2a, 5.3a, 5.4a, and 5.5a) correspond to surfaces

124 with low Fl value -- surfaces which are most likely frozen -- and freeze/thaw boundaries appear as a fuzzy white lines around these frozen regions. The dark pixels in the 37 GHz images (panels b, d, and f in Figures 5.2, 5.3, 5.4, and 5.5) correspond to surfaces of low 37 GHz radiobrightness. The fuzzy white lines around these dark regions are the boundary pixels that equal P37. Some or all of these boundaries correspond with the coarse-resolution freeze/thaw boundaries of the freeze maps. Similarly, the dark pixels in the spectral gradient images (Figures 5.2c, 5.3c, 5.4c, and 5.5c) correspond to surfaces with low spectral gradient, and the fuzzy white lines are boundary pixels that equal PsG. In all images, regions of no data are shown as white. Freeze maps are composite TB(37) and aTBlaf coarse-resolution images, and all boundaries on freeze maps are freeze/thaw boundaries. The grey-scale values at each pixel on a freeze map, IFM, are given by, IFM = max(137,ISG), (5.2) where 137 and IsC are the corresponding 37 GHz radiobrightness and spectral gradient grey-scale values. 137 and IsG are given by, TB(37)- P37 137= - (5.3a) 10 037 Or /af)ISG aTB )-PSG (5.3b) where P37 and P are the 37 GHz radiobrightness and spectral gradient decision criteria (Table 4.4), and 037 and GsG are RMS values of 37 GHz radiobrightness and spectral gradient measurement uncertainty (i.e., RMS Temperature Resolution in Table 4.1). I37,7 IS, and IFM

125 are simultaneously determined at coarse-resolution for each pixel, and are clipped at ~1. As a result, whenever I37, ISG, or IFM <- 1, the corresponding pixel is black, and whenever I37, ISG, or IFM >+1, the corresponding pixel is white. As mentioned in chapter 4, the sampling intervals for the 37 GHz, 18 GHz, and 10.7 GHz SMMR channels are 30 km, 60 km, and 97.5 km respectively, and are approximately the SMMR beamwidth (i.e., resolution) on the surface. As a result, resolution compensation combines approximately (97.5/30)2 = 3.252 fine-resolution data samples, or (97.5/60)2 = 1.6252 medium-resolution data samples, to form each coarse-resolution data sample. Assuming independent measurement uncertainties for each sample, we estimate RMS values for measurement uncertainty in coarse-resolution 37 GHz, 18 GHz and 10.7 GHz SMMR channels by, TRS(37) 037 = 0.5 K TR (1 8) 018 = 16 -075 K 1.625 o10.7 = TR(10.7) = 0.9 K. TRSs(37), TR7S(18), and TRs(10.7) are the temperature resolutions of the 37 GHz, 18 GHz, and 10.7 GHz SMMR channels (Table 4.1). Assuming independent measurement uncertainties in each channel, equation (4.27) yields, (15.12a37 + 3.92018 + 11.22o0y7)l/2 SG =368.66 = 0.035 K/GHz. Thus, the fuzzy white boundaries in the spectral gradient images (Figures 5.2c, 5.3c, 5.4c,

126 and 5.5c) correspond to aTBlaf within ~0.035 K/GHz of PSG, and the boundaries in the 37 GHz images (panels b, d, and f in Figures 5.2,5.3,5.4, and 5.5) correspond to TB(37) within ~0.5 K of P37. The white boundaries in the freeze map are composites of these TB/laf and TB(37) boundaries. Comparing the "unrefined" images of (Figures 5.2 and 5.3) with the "refined" images (Figure 5.4 and 5.5) shows that refined criteria generate coarse-resolution, 37 GHz boundaries that are located more closely to freeze map and spectral gradient boundaries. Moreover, refined fine-resolution, 37 GHz boundaries (Figures 5.4f and 5.5f) are more consistent with ground data' than are unrefined 37 GHz boundaries (Figures 5.2f and 5.3f). Thus, freeze/thaw boundaries derived from refined criteria should be more accurate than those derived from unrefined criteria. In the refined images of Figure 5.4, most sections of the coarse-resolution, 37 GHz boundary in the northwest corner of the test area (Figure 5.4b) correspond with boundaries of the freeze map (Figure 5.4a). These sections of 37 GHz boundary would be designated as freeze/thaw boundaries. None of the two other boundaries in Figure 5.4b correspond to any freeze map boundary, and are probably wet/dry boundaries. The freeze/thaw boundary in the coarse-resolution, 37 GHz radiobrightness image also corresponds to boundaries in medium-resolution and fine-resolution, 37 GHz images. That is, medium-resolution and fine-resolution freeze/thaw boundaries are the convoluted boundaries in the northwest corner of Figures 5.4d and 5.4f, respectively. Some boundaries are formed at fine-resolution that do not correspond to any boundary observed at coarse-resolution. These boundaries appear around dark radiobrightness "islands" in Figure 5.4f, and cannot be identified on the basis of the available information. Such boundaries are not part of the freeze/thaw boundary estimates. 1 See Tables 5.3 and 5.4.

127 -- - -9. N (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution - - -- F --- (c) Spectral gradient at (d) 37 GHz radiobrightness at coarse-resolution medium-resolution AIR TEP Thaw 0 Mixed Freeze ~ 0 QSOL TP Thaw 0 Mixed D Freeze ~ hi~d i 13 P r.... 0 II m * U (e) Air and soil temperatures (f) 37 GHz radiobrightness at fine-resolution Figure 5.2 A comparison of reported air and soil temperatures with images of North Dakota and the surrounding region. Boundaries were determined using unrefined freeze/thaw criteria. Data were collected at midnight, October 24, 1984.

128 'I im (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution II. ___ ___ ___ ___ _is__ J ___1. *4 (c) Spectral gradient at coarse-resolution (d) 37 GHz radiobrightness at medium-resolution (e) Air and soil temperatures ( 37 GHz radiobrightness at I fine-resolution Figure 5.3. A comparison of reported air and soil temperatures with images of North Dakota and the surrounding region. Boundaries were determined using unrefined freeze/thaw criteria. Data were collected at noon, December 11, 1984.

129 4' (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution - iT- - - _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ I (c) Spectral gradient at coarse-resolution (d) 37 GHz radiobrightness at medium-resolution * (e) Air and soil temperatures (f) 37 GHz radiobrightness at fine-resolution l Figure 5.4. A comparison of reported air and soil temperatures with images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, October 24, 1984.

130 -~ (a) Freeze map at coarse-resolution I-:.i;.:..:.:...:.:.:..:. I (b) 37 GHz radiobrightness at coarse-resolution..l ------- ". I, (c) Spectral gradient at coarse-resolution (d) 37 GHz radiobrightness at medium-resolution II 9 \U Oo a0 0 ~ L'q ---- -- 0____________,,,_ AIR TEMP Thaw 0 Freeze Frous 0 i (e) Air and soil temperatures (f) 37 GHz radiobrightness at fine-resolution Figure 5.5 A comparison of reported air and soil temperatures with images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at noon, December 11, 1984.

131 In Figure 5.5, all coarse-resolution 37 GHz radiobrightness boundaries (Figure 5.5b) correspond the freeze map boundaries (Figure 5.5a), and all fine-resolution boundaries (Figure 5.5f) become freeze/thaw boundary estimates. Midnight and noon fine-resolution freeze/thaw boundaries (Figures 5.4f and 5.5f) are consistent with midnight and noon ground data (Figures 5.4e and 5.5e). 5.3 Automated Boundary Processing Figure 5.6 represents automated boundary localization for SMMR data of midnight, October 24. Figure 5.6a repeats the freeze map of Figure 5.4a, and Figure 5.6b shows the associated coarse-resolution, 37 GHz radiobrightness image. As before, 37 GHz boundaries are composed of pixels whose 37 GHz radiobrightness equals P37. However, the 37 GiHz boundaries in Figure 5.6b consists of (fuzzy) white and black sections. Pixels along white boundaries have spectral gradients that are less than or equal toPSG. That is, white boundaries are most likely to be freeze/thaw boundaries. Pixels along black boundaries have larger spectral gradients and are less likely to be freeze/thaw boundaries. In the medium-resolution, 37 GHz radiobrightness images (Figure 5.6c), white boundaries are 37 GHz boundaries that are migrations of white boundaries at coarse-resolution (Figures 5.6b). White boundaries at medium-resolution are (spatially) near white boundaries at coarse-resolution. To determine boundary migration, we solve a 2-D Taylor expansion in an analogous process to edge focussing [5]. We first determine x and y coordinates of each pixel along all boundaries in 37 GHz, medium-resolution images. These coordinates are denoted as (xi,yi) for the i" boundary pixel. We also determine coordinates of each pixel along allfreeze/thaw boundaries in 37 GHz, coarse-resolution images. These coordinates are denoted as (xj,yj) for the jh

132 freeze/thaw boundary pixel. Using the development of Appendix C, a medium-resolution boundary pixel, (xi,yi), is a migration of coarse-resolution, freeze/thaw pixel, (x;, yj), if it satisfies the Taylor expansion, ds(xi, y) ds(xi, y) S2 =l+ dx (i-x)+ dy (Yj-Yy), (5.4) where s1 and s2 are the medium and coarse resolution levels, respectively (Table 5.1). The function s(x, y) is defined implicitly by (Appendix C), r(x, y,s(x,y)) = P37, (5.5) where r(.,.,S) is the 37 GHz radiobrightness at resolution S, so that by the implicit function theorem [68], ds(xi,yi) -rx,(x s (5.6a) =- (5.6a) dx r,(xi, Yi, S1) ds(xi,yi) -ry(xi,yi,sl) - (5.6b) dy r,(xi,Yi,sl) where variable subscripts denote partial differentiation. Rotationally symmetric Gaussian filtering yields, r(Xi, yj, S1) = s[r,(xj, yj, S1) + ryy(xj, yj,)] (5.7) so that the medium-resolution pixel (x,, yi) is a migration of the coarse-resolution pixel (xi, yj)

133 if = Sx(xi, i, )j xi r(xirYi,is) SI [r(xiyisS)+ryy(xi(yis) j S -slr,(xi, y,,s) + r(xi, ys)] (i - ) (5.8) is satisfied. That is, white pixels in the medium-resolution image satisfy (5.8) for white pixels the coarse-resolution image. An analogous process is performed for medium-to-fine resolution boundary migration. In this case, (xi, y,) values are determined for boundaries at fine-resolution and (xj, yj) values are determined for white boundaries at medium-resolution. Also, sI is the resolution level at fine-resolution, and s2 is the medium resolution level at medium-resolution. The accuracy of this boundary migration process depends on how accurately s(x,y) can be approximated as planar in the vicinity of boundaries. This process becomes more accurate with increased sampling with respect to scale. That is, the use of more images at smaller scale changes improves the accuracy of boundary migration, but increases the computational load. All boundary pixels of Figure 5.6 are determined using a hysteresis process to account for measurement uncertainty. That is, black and white boundary pixels are within ~a of P37 where, from section 5.2, a = 0.9 K for the fine-resolution 37 GHz image (Figure 5.6d), a = 0.75 K for the medium-resolution 37 GHz image (Figure 5.6c), and a = 0.5 K for the coarse-resolution 37 GHz image (Figure 5.6b). Moreover, 37 GHz radiobrightness values must be greater than P37 + a on one side of a boundary, and less than P37 - a on the other side.

134 ~XCP:ZI"I: ~:~.: : i jj ~:~:~::.b.. iii::::::;: i::: ~:~:~:~:~:~:::7r:ii~:::: i f:::::::: i.:i............ 1 (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution -— 1 -0 tn: i. I (c) 37 GHz radiobrightness at medium-resolution,i ii (e) Classified frozen ground at fine-resolution (d) 37 GHz radiobrightness at fine-resolution Figure 5.6. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, October 24, 1984.

135 Frozen terrain is identified iteratively using fine-resolution, 37 GHz data. First, pixels along white, fine-resolution boundaries (Figure 5.6d) are identified as "frozen" pixels. Second, pixels whose 37 GHz radiobrightnesses are less than or equal to P37, and are contiguous to frozen pixels, are also identified as frozen. Third, the previous step is repeated until no additional pixels are identified as frozen. Fourth, the resulting collections of frozen pixels constitute regions of frozen terrain. Using this procedure, terrain identified as frozen are indicated by whitened regions in the northwest corner of the Figure 5.6e. Because freeze/thaw boundaries must be closed contours, the final freeze/thaw boundary (i.e., the edge of the identified frozen region) contains boundary sections that did not, previously, show strong freeze/thaw boundary indications. Nonetheless, the final freeze/thaw boundaries of Figure 5.6e are the best fine-resolution estimates of the actual freeze/thaw boundaries using available data. Table 5.2. Time summary for images of Figure 5.7 through Figure 5.15; measurement interval is the time interval between present and previous measurements. Measurement Measurement Measurement Figure Date Time-of-Day Interval (Days) 5.7 October 24 Midnight 5.8 October 30 Midnight 6 5.9 November 1 Noon 2.5 5.10 November 5 Midnight 3.5 5.11 November 27 Midnight 22 5.12 November 29 Noon 2.5 5.13 December 3 Midnight 3.5 5.14 December 9 Midnight 6 5.15 December 11 Noon 2.5 Boundary migration and automated classification was applied to the images of Figure 5.7 through Figure 5.15. These data were obtained at irregular intervals from October 24, 1984 through December 11, 1984. Time summaries of the data are given in Table 5.2, and

136 weather summaries are given in Tables 5.3 and 5.4. These images show the growth and contraction of ground-freeze from October 24 to November 5, and again from November 27 to December 9. After December 9, the area remains frozen through the end of December. Table 5.3. Air temperatures for October through December measurement dates. (N) and (D) indicate night and day temperature measurements. Measurement times are shown with data site names. (S) indicates snowfall at measurement time. I,,,,,,,, A TD) TlIU l4DD A rtr TDTEiC 1O0^\ A I Il1lr.IA-r A I U L"i- % b) Data 10/24 10/30 11/1 11/5 11/27 11/29 12/3 12/9 12/11 Site (N) (N) (D) (N) (N) (D) (N) (N) (D) Aberdeen -0.5 -3.3 -6.1 -3.3 -0.5 0.0 -15.5 -1.1 +2.8 (N): 12 AM (S) (D): 12 PM Bismark +0.5 -7.2 -8.9 -3.9 -2.8 -2.2 -9.4 -0.5 -5.5 (N): 12 AM (S) (S) (S) (D): 12 PM Fargo +0.5 -4.4 -9.4 -4.4 +0.5 -1.7 -17.7 -1.7 -0.5 (N): 12 AM (S) (D): 12 PM Huron +2.2 -0.5 -3.9 -3.9 -0.5 +1.7 -15.0 +1.7 +4.4 (N): 12 AM (S) (D): 12 PM Miles City +1.7 -10.5 -9.4 +2.2 -7.8 +0.5 -19.4 -3.3 -8.9 (N): 11 (S) (D): 11 AM Rapid Cid. +2.2 -4.4 -2.2 +0.5 -2.2 +8.3 -11.7 +6.1 +2.8 (N): 1 (D): 11 AM Williston -1.1 -12.8 -13.9 -5.0 -7.8 -6.6 -12.8 -2.2 -13.3 (N): 12 AM (S) (S) (S) (D): 12 PM

137 Table 5.4. Soil temperatures for October through December measurement dates. (N) and (D) indicate night and day temperature measurements. Measurement times are shown with data site names. Measurements at Morris were maximum and minimum daily temperatures. Measurements at Thief River were at 2 inch depth; all other sites were at 5 cm depth. | __-SOIL TEMPERATURES (~C) Data 10/24 10/30 11/1 11/5 11/27 11/29 12/3 12/9 12/11 Site (N) (N) (D) (N) (N) (D) (N) (N) (D) Thief River +3.9 0.0 -- -2.2 -3.3 -- -5.0 -- (N): 7 AM (D): — ___________ Morris 0.0 -1.1 -1.1 -2.8 -1.1 -1.1 -2.2 -4.4 -2.8 (N): Min (D): Max_______ __ Cavalier +4.4 +3.3 -- +0.5 -1.1 -- -1.1 -1.1 -- (N): 8 AM (D): — _____ _ Brookine +2.8 +2. 2 +0.5 +1.7 +1.7 +0.5 -1.7 0.0 +0.5 (N): 8 AM (D): 5 PM Cottonwood +1.7 +1.1 +7.8 0.0 -1.1 -2.2 -- -- -- (N): 7 PM (D): 5 AM____ Eureka 0.0 -1.1 -2.2 -2.2 -3.3 -3.3 -6.6 -6.6 -5.5 (N): 8 PM (D): 6 AM____ __________ Highmore +1.7 -1.1 -1.1 0.0 0.0 -1.1 -2.8 -2.8 -1.7 (N): 8 AM (D): 5 PM____ Milesville +3.3 -2.2 +5.0 +2.8 -0.5 -1.7 -1.1 -1.7 -1.7 (N): 8 AM (D): 6 PM_____________ Newell +1.1 -0.5 +1.7 0.0 -2.2 -2.2 -2.2 -2.8 -2.8 (N): 8 AM (D): 6 PM_ _____ __ _ Redfield -2.2 -3.3 -1.1 -3.9 -3.9 -3.9 -5.0 -5.0 -5.0 (N): 8 AM (D): 5 PM_____ Waubay 0.0 -1.7 -2.8 -3.9 -3.3 -3.9 -3.9 -5.0 -3.9 (N): 8 AM (D): 5 PM

138 I............. (a) Freeze map at coarse-resolution.0 (b) 37 GHz radiobrightness at coarse-resolution I i (c) 37 GHz radiobrightness at fine-resolution (e) Air and soil temperatures (d) Classified frozen ground at fine-resolution Figure 5.7. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, October 24, 1984.

139 t _ -----— I-I-I --- —I-I-I - I ~ 1 I a1 (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution - -- -- -- (c) 37 GHz radiobrightness at fine-resolution EN - -- i (e) Air and soil temperatures (d) Classified frozen ground at fine-resolution Figure 5.8. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, October 30, 1984.

140 -1 (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at Ii coarse-resolution (c) 37 GHz radiobrightness at fine-resolution A IR TEM P:::.::.:..:.:::::..:.:..:....:.::::::::::.::..............T P Thaw 0.....................Q............. Fruss ~ -----— 0 I ~ I............iii..i h A x * * *le'l - **M1d C......... (e) Air and soil temperatures (d) Classified frozen ground at fine-reesolution.,,...........................U.......................... -...............0......*.....................B............. (e) Ar an soiltempraturs (dClasifie frozen...ground at... fine....... resolution................. i Figure 5.9. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at noon, November 1, 1984.

141 Is. - I - (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution so (c) 37 GHz radiobrightness at fine-resolution;. * I * 0 B O - - -0 O B 0 AI TEMP Thaw O Mxed 9 Freeze ~ SOL TEP Thaw D MFrzd U Freeze ~ (e) Air and soil temperatures I (d) Classified frozen ground at fine-resolution Figure 5.10. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, November 5, 1984.

142 (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at _____ coarse-resolution_ Ng I I I (c) 37 GHz radiobrightness at fine-resolution 1. AIR TEMP Thaw 0 Mixed e Freeze 0 SOL TEMP Thaw o Mixed f Fromm ~ II1 (e) Air and soil temperatures (d) Classified frozen ground at fine-resolution Figure 5.11. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, November 27, 1984.

143 -- (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution — l (c) 37 GHz radiobrightness at fine-resolution AR TEMP Thaw O lixed l fn-Freeze.es ~..........T...........P (e) Air and soil temperatures (d) Classified frozen ground at fine-resolution Figure 5.12. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at noon, November 29, 1984.

144 (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution -4 _a (c) 37 GHz radiobrightness at fine-resolution (e) Air and soil temperatures (d) Classified frozen ground at fine-resolution Figure 5.13. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, December 3, 1984.

145 (a) Freeze map at coarse-resolution I I (b) 37 GHz radiobrightness at coarse-resolution (c) 37 GHz radiobrightness at fine-resolution (e) Air and soil temperatures (d) Classified frozen ground at fine-resolution Figure 5.14. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at midnight, December 9, 1984.

146 11I I (a) Freeze map at coarse-resolution (b) 37 GHz radiobrightness at coarse-resolution (c) 37 GHz radiobrightness at fine-resolution.....AIR TEMP Thaw 0 Fr.ze... Thma 0 -... F..wze... fine-resolution ()A a si C frozen.ground.at fine resol uion i.................. i i:ei Figure 5.15. Automated images of North Dakota and the surrounding region. Boundaries were determined using refined freeze/thaw criteria. Data were collected at noon, December 11, 1984.

CHAPTER 6 SUMMARY AND FUTURE WORK 6.1 Summary This isi a scale-space based approach for analyzing spatial features in multiresolution satellite images. Our objective was to estimate locations of freeze/thaw boundaries in the northern Great Plains using multispectral data from the scanning multichannel microwave radiometer (SMMR). Scale-space techniques were used to integrate (fuse) the multispectral SMMR data ata coarse-resolution, and estimate freeze/thaw boundaries at fine-resolution. Image data collected by devices such as SMMR were modelled as multiple resolution data, where the receive apertures approximate linearly-scaled, spatial filters. Variations in filter scale occur as a result of variations in receive frequency, so that the different SMMR channels have different resolution. Scale-space filtering is a multiresolution analysis technique, and was used as a framework for integrating SMMR channel data. Because standard scale-space theory requires Gaussian spatial filtering, scale-space theory was extended to signals generated by physically realizable (i.e., non-Gaussian) filters. Without such extensions, scale-space filtering could not form the basis of data integration. Extrapolation-in-scale was developed by applying scale-space filtering to boundary estimation in multiresolution systems. Boundary estimates made at coarse resolution have significant location errors. Extrapolation-in-scale was used to extrapolate boundary estimates from coarse-to-fine resolution, and mitigate boundary location errors. The conditions for exact extrapolation were given in the form of theorems. 147

148 The freeze indicator was developed to classify frozen and thawed terrain, and to locate freeze/thaw boundaries. The freeze indicator is a combination of spectral gradient and 37 GHz radiobrightness from SMMR data. Decision criteria for surface classification using were derived from data clustering. Using the freeze indicator and 37 GHz radiobrightness from SMMR data, freeze thaw boundary locations were estimated at fine-resolution. In this process, coarse-resolution, freeze/thaw boundary locations were determined using the freeze indicator, and then extrapolated to fine-resolution using 37 GHz radiobrightness data. While there were inadequate ground data to verify the accuracy of this process, the process was consistent with trends in the ground data. There are three specific contributions of this thesis. First, in the area of remote sensing, the freeze indicator is a new parameter for indicating the presence of frozen terrain. While the SMMR is no longer operational, the techniques for determining freeze/thaw boundaries presented in this dissertation can be adapted to the newly operational Special Sensor Microwave Imager (SSM/I) [88], which flies on a series of Defense Meteorological Satellites. Using the SSM/I with these modified techniques, it should be possible to monitor the state of surface moisture over the large surface areas required by climatological and hydrological models. Second, in the area of signal processing, scale-space techniques were extended to signals generated by non-Gaussian filters. By requiring Gaussian filtering in its original development, scale-space filtering was constrained to sensing problems where filtering was computer-implemented (e.g., computer vision). By extending scale-space filtering to include non-Gaussian filtering, scale-space techniques can be applied to sensing problems where spatial filtering is implemented by a variety of physical devices, such as antennas and lenses.

149 Third, also in signal processing, extrapolation-in-scale was developed. In general data fusion problems, data is integrated from sensors of different resolution and fine-resolution information is often lost. Using extrapolation-in-scale, it is possible to analyze and retain fine-resolution information in these general data fusion problems. 6.2 Future Work In this dissertation, decision criteria for locating freeze/thaw boundaries using the freeze indicator were determined from clustering and unsupervised classification (chapter 4). An alternative and, potentially, more accurate approach would use supervised classification. Accurate ground data necessary to train a supervised classifier were not available for the dissertation. Thus, a more detailed examination of emissions from freezing ground would provide such ground data and, therefore, improve decision criteria for locating frozen ground [31]. Some freeze/thaw boundary contours produced by the freeze indicator did not register well with 37 GHz boundary contours at coarse-resolution (chapter 5). That is, threshold-crossings of 37 GHz radiobrightness did not adequately extrapolate all boundary locations identified by the freeze indicator. However, it may be possible to formulate other signals for extrapolating those boundary contours where the 37 GHz threshold-crossings are inadequate. For instance, zero-crossings of 37 GHz gradients (i.e., the spatial gradients of 37 GHz radiobrightness) might be used to extrapolate boundaries between frozen and very moist (thawed) surfaces. Thus, a series of such "complementary" signals could be studied for extrapolating freeze/thaw boundaries where conditions of the underlying surface (e.g., moisture content) vary. However, such complementary signals should be studied in concert with a detailed examination of surface emissions.

150 In addition, estimates of freeze/thaw boundaries did not always form closed-contours at fine-resolution (chapter 5). While closed boundary contours were constructed during the classification of frozen and thawed ground, all such contours may not have been optimal. Thus, it would be useful to study techniques for determining optimal closed-contour, fine-resolution estimates of freeze/thaw boundaries. Such techniques have been developed for locating the auroral oval in satellite images [59], and are derived from edge contour models of computer vision [44]. It should be possible to apply similar techniques to the determination of freeze/thaw boundaries. Finally, future work should expand the application of boundary estimation and extrapolation techniques developed in this dissertation. For example, these techniques can be applied to the data fusion of synthetic aperture radar (SAR) data gathered over multiple frequencies, as well as the fusing of radar and radiometer data [87].

APPENDICES 151

152 APPENDIX A: Conditions for embedding structure in function-crossing fingerprints In considering the function-crossing problem, we follow the procedure, of Yuille and Poggio [84,85] but consider equations of the form, e(x,s) =(x) (Al) where x) is a real, continuous function of x -- which we call the thresholdfunction -- and e(x,s) is given by, e(x,s) =O{r(x,s)}, where O{.} is a general linear operator and r(x,s) is a filter output signal. We extend scale-space results for the function-crossing problem by first considering an arbitrary linear operator O{. }, rather than the differential operator L {.). The linearity of the operator, rather than any specific form, is seen as an important restriction by considering O{r(x,s)} =0{h(x,s)*i(x)} (A2) = h(x,s)*O{i(x)} = h(x,s)*j(x) for some function j(x) (the existence of j(x) is assumed). From (A2) O{r(x,s)}= = = h(x,s)*j(x)=0 O (A3)

153 so that applying standard scale-space fingerprint results to (A3), it follows that fingerprints with embedded structure are generated for any linear operator in the same manner as they are for any linear differential operator. Therefore without loss of generality, the use of a general, linear scale-space operator is assumed. The necessary conditions for the existence of fingerprints with embedded structure, {SS1)-{SS3} (chpater 2), are assumed to hold for y(x)= O. Of conditions {SS 1 -{SS3}, only (SS2) and (SS3) are affected by the shape of the threshold function, and must be re-derived for y(x) ~ 0. In the case of (SS2 }, (Al) can be equivalently written as, f(x,s) - e(x,s)- y(x) = O, (A4) so that by the implicit function theorem [36], if f(x,s) f,(xs) 0, at the fingerprint extrema of (A4), then the fingerprints of (A4) are differentiable and { SS2 } holds for y(x)0 O. Since y(x) is a real, continuous function (i.e., |y'(x) < oo ) and, by assumption, the implicit function theorem holds for y(x) = 0, we have, f,(x,s)= e,(x,s)O (A) at an extremum of (A4), so the fingerprints are differentiable and (SS2) holds for y(x) ~ 0.

154 In the case of {SS3}, a set of constraints must be established to guarantee that all extrema of the fingerprints of (Al) are maxima (i.e., (SS3 is guaranteed for 'yx) 0). Following the development of Yuille and Poggio, (Al) is differentiated with respect to x, yielding the expression, y (X) = e(x,s) + e(x,s) d(( )) (A6) where, 3e(x,s) e ) e (x,s)ax ) ae(x,s) as,(x) d=(x) dx for a fingerprint s(x). Extrema of s(x) exist at the points where ds(x)/dx=O, so that evaluating (A6) at an extremum xo yields, y'(xo) = ex(xo, so) (A7) where, So = S (Xo).

155 The existence of an extremum at Xo is assured provided that, I '(x) 1oo (A8a) I e,(x,s) < oo (A8b) e,(x,s) 0 (A8c) in a neighborhood of the extremum Xo. These conditions, in turn, impose restrictions on the threshold function, filter, and input signal. The inequality of (A8a) is satisfied for all x by using a real, continuous threshold function. The inequality of (A8b) is satisfied by using a filter impulse response h(x,s) or an input signal i(x) having non-infinite bandwidth, so that 3e(x,s)lax j< oo. Condition (A8c) is given by (A5). Differentiating (A6) with respect to x gives, ds d2s y"(x) =e,(x,s)+ e(x (x,s) e,(,) ds ds + [e,,(xs)+ e,(xs) dJ (A9) where d2^x) r "(x) - -(x2 dx2 ( a (ae(.), for arbitrary variables u and v. Evaluating (A9) at the extremum point x0 gives,

156 d2s (Xo) y"(xo) = e,.(xo, o) + e,(xo,sO) ). (A10) dx2 All extrema of s(x) are maxima if d2s(x)/dx2<0 for all x where ds(x)/dx=0. Applying this condition to (Al), (A7), and (A10) provides the constraints necessary for all extrema of s(x) to be maxima, e(x,s) =(x) (Alla) ex(x,s) = '(x) (All b) (x)- ex,s) X-,- s) < 0. (A 1lc) e,(x,s) Thus, the conditions of (Alla)-(Allc) must hold for embedded structure to exist. Furthermore, the conditions can be normalized with respect to les(x,s)l without loss of generality since, by (A8c), e,(x,s) O0. The result is a set of constraints for the fingerprints of (Al) such that for eS(x,s)=l, e(x,s) =(x) (A12a) e,(x,s) = y'(x) (A12b) e,(x,s) > y"(x) (A12c) and for eS(x,s)=-l, e(x,s) =y(x) (A13a) ex(x,s) = y'(x) (A13b)

157 e,(x,s) <y "(x) (A13c) Therefore, whenever a threshold function y(x), operator Of{.), and filter h(x,s) exist which meet the conditions of expressions (A12a)-(A12c) or expressions (A13a)-(A13c) -- depending on the sign of e,(x,s) -- the fingerprints of (Al) have embedded structure with respect to scale. The existence of embedded structure for a given h(x,s) is demonstrated by negating a set of conditions that are complementary to (A12a)-(A12c) or (A13a)-(A13c). A set of complementary conditions to (A12a)-(A12c) is, e(x,s) =T(x) (A14a) ex(xs) = yT(x) (A14b) e.(x,s)= "(x)-m2 (A14c) e(x,s) = 1. (A14d) Similarly, a set of complementary conditions to (A13a)-(A13c) is, e(x,s) =y(x) (A15a) e,(x,s) = y'(x) (A15b) e, (xs) = Y"(x) + m2 (A15c) e,(x,s) = -1 (A15d) In both (A14c) and (A15c), m is an arbitrary real number.

158 Negating the complementary set (A14a)-(A14d) is accomplished by introducing constants k,, k2, k3, and k4 and combining (A14a)-(A14d) to give, kle(xs)+ (xs+) +k3ex, (xs)+k = klx)+k2y'(x) +k3y"(x)- k3m+ k4, such that, k,y(x) + k2y'(x) + k3y"(x) - k3m2 + k4 ~ 0 (A16a) for a function e(x,s) that satisfies, kle(x,s) + k2ex(x,s) + k3e,(x,s) + k4e,(x,s) = 0. (A16b) Thus, finding a solution to (A16a) and (A16b) implies that embedded structure exists for the threshold function y(x) when es(x,s)=l. A similar manipulation of the expressions (A15a)-(A15d) results in kle(x,s)+ k2ex(x,s) + k3e(x,s) + k4es(x,s) = ky(x) + k2y'(x)+ k3^"(x) + k3m2- k4 Thus, embedded structure exists for e,(x,s)=- if constants k, k2, k3, and kand a function e,(x,s) can be found such that, kyx)+ k2(x+k() + kx) + k(x) + m2- k4 O (A171a)

159 and kle(x,s) + k2e(x,s)+ k3e,(x,s) + k4e,(x,s) = 0. (A17b) Further restrictions on the constants k,, k2, k3, and k4 are imposed by condition { SS 1b) on the filter impulse response (chapter 2). In order that (A 6b) and (A17b) are dimensionally correct, k,, k2, k3, and k4 must be written as a/s2, b/s, c, and -d/s, respectively, where a, b, c, and d are arbitrary constants. This substitution explicitly incorporates the scaling parameter s to yield one equation and two constraints for (A16a), (A16b), (A17a), and (A17b). Thus, the function e(x,s) must be found which solves, a b d -e(s)+e(x,s)+-ex,s)ce,(x,s) —e(x,s) = (A18) S S S while meeting the constraint, a b d Y(X)+ y'(X) + Cy"(X) cm2-__ 0 (A19a) S S S for e,(x,s)=l and, a b d -y(X) +-'(x)+cy"(x)+cm2 + 0 (A19b) sS S S for e,(x,s)=- 1.

160 Given a general linear operator O{. } and input signal i(x), the identity e(x,s) O {h(x,s)*i(x)} = h(y,s)O{i(x-y)}dy. reduces (A 18) to, a b d a-h(x,s) +- h(x,s)+ch(x,s) — h,(x,s) = 0 s2 s xxs (A20) Thus, the solution of (A20) for h(x,s) under the conditions of either (A19a) or (A19b), depending on the sign of es(x,s), will ensure that tt the fingerprints of (Al) will have embedded structure with respect to scale.

161 APPENDIX B: Proof of Theorem 3-1 (1-D) Define the functions r(x(s),s) and p(u(s),s) r(x(s),s) h(x(s),s)*i,(x(s)) p(u(s),s) - h(u(s),s)*i,(u(s)), where the contours x(s) and u(s) are implicitly defined by, r(x(s),s) = a(x(s)) p(u(s),s) = j(u(s)). A Taylor expansion of x(s) about so yields, (s -So) x(s) = x(So) + (s - So)x'(So) + 2 x ) +.... By the implicit function theorem [36] and repeated differentiation, r,(XO, S0) ax(xo) - rX(xo, So) x "(so) (r=(xo, So) - ac,(x)) (x'(so))2 + 2rX(xo, so)x (so) + rSs(Xo, so) ao(xo) - rx(xo, so) (B 1 a) (Blb) (B2) (B3a) (B3b)

162 and similarly for x()(so), where, Xo x(So). (Subscripted variables indicate partial differentiation.) The threshold function a(.) is linear, so that o(x) =Aix + A (B4) for constants A1 and Ao, and ac,(.) = 0. The filter kernel h(x,s) is Gaussian, so by solution to the heat equation, r (xO, so) = sor(Xo, So), and, x(so) - sor,=(xo, so) x A - r(xo, So)' x "(o) = (B5a) (B5b) r,(Xo, So) (x'(so))2 + 2sorm(xo, so)x'(So) + so2r (Xo, So) Al - rx(xo, so) and similarly for x(n)(so).

163 Repeating the steps above for u(s) yields, (s -So)2 u(s) =u(so)+( (s -So)u'(so) + 2 u "(So) +... where, sop.(Uo, so) () B1 - p(u, So) u "(s) = p.(uo, So) (u '(so))2 + 2sopp(uo, So)u '(so) + so2p(uo, So) B1 -Px(uo,so) and similarly for u(n)(so), where, Uo u (So). The threshold function,3(.) is linear, so that P(x) =Bx +Bo (B6) (B7a) (B7b) (B8) for constants B1 and Bo, and = 0. The constant Bo is chosen so that uo=xo. Thus, for u(s) to approximate x(s), it is sufficient that, x'(sO) = u (sO), x"(s) = u (s)...

164 so that (B5a), (B5b),.. and (B7a), (B7b),.. yields, A - h(xo~ so)* di,(xo) h (X0, So)* dx- i4(x0) B~~o~o*ixo)- n=2,3,4. QED

165 APPENDIX C: Proof of Theorem 3-2 (2-D) We are concerned with contours (x,s) and (_,t) in 3-space that satisfy the function-crossing relation, r (,s) = a() p(x,t) = x), (Cla) (Clb) where, XFunctions r(, s) and p (, t) are defined as', r(jx,s)- h(x,s)**i,(x p(x, t) = h(x,t)**i(x), and threshold functions a(.) and P(.) are given by, ax() =Ao +Alx, + Ax2 +A3xx2 Px) = BoB + BlXl + B2 +3xx2, (C2a) (C2b) 1 "**" denotes 2-D convolution.

166 where Ao, A1, A2, A3, Bo, B1, B2, and B3 are constants. Equivalently, the implicit function theorem [68] can be used to define contours s(s) and tx), in the vicinity of the point x~, so that (Cla) and (Clb) can be written as, r(,s(x)) = a(x) (C3a) P (x, t (x) = ). (C3b) The constants Ao and Bo are chosen so that, s(~) = t). (C4) Taylor expansions of s() and t(x) about x~ yield [36], s Ux = s C~) + I ) (Xi - Xi~) + i = 1 Xi 2 t) +2! Sa )(X i- zo)(x.-x + (C5a) tX=t(~)+ )(xi-x~) +2 a (Xi - Xio) i=1 Xi; ~'+ - -Xi ' (C5b) 2!ii;Taxtla(~ where s (~) = t(x~). Thus, if V K, aS(X ) t~) * * axilax2... axi,, i ax,...aX li-12 aX

167 then s Ux =tx. We define the function IukO.,.) by F-OkX0 s 0) =r (, s 0) - (~ ) which satisfies the initial condition, s0= =s(X0). By the implicit function theorem [68], O=17- (4 so) + S, (j) uk! s0)3, o O=F Cx, s 0) + sx2x0) Ji.0@10 0 s 0), (C6a) (C6b) so that', s1 0 - (I-Io, sX0) I M0(!0 s 0) F1(X02s0) where subscripted variables indicate partial differentiation. From the definition of JIG'(,.), and,

168 r~ X~s)= r,-O For a rotationally symmetric, Gaussian kernel h (xj, s (Ux), and (C6a) and (C6b) can be written as, o = [r~iY s0) + s0 [r(o s)+ s)]o o = [~r2@O, s0 - +s~ [r(~, 0) + r.2X2(-j0, SO)]o where, 0 I0 Repeating this process for tCx),, and equating, yields, (C7a) (C7b) rx, (x 0, S 0) - ax W PXI(CYJ SO) J3(xi ) r~l, (Xj03, S0) + r,@i0, 0 0) = X px1I1@_0, s0) + pzrX 01S0) as (x0)- at(xj0) a, ax, I (C8a) and, r 2~0, 0O - N(X@j0) p(0, 0O - I32(X0) rx1,@j, CIO s) + r"21G- (X-, S ) =px1@il WIO ) + PW, @Y 0O) aX2 DX2 (C8b)

169 We now define functions rF(.,.) and r12(.,.), r1,,oso) = i,(oso) +s(x~) r(X-, s), r2l 0, s0) =,( s0) + (oSo)+ o) rXo,0 s). Applying the implicit function theorem to (C6a) and (C6b) yields, 0=__,. o(o)+sx0),(~o,, o = Fe(x0,so)+sx~to) Jj(~,s0), o =. _ s, 0) + s( ~.O) rl20, s), and expanding terms yields, 0 = ^sxl~~) @i) + sx~,) Fj, s) + x~. _~,)+x,) l + + slo)[^l(XOSO)+Sx(O) r.12(, o)] 0=s ) ij= = 2, s~) )+ O) 0) + sx(oo)+ ) l.o)+ +SX (o)[r, c@o, So) + SX l(o) r. cXo, So)] o= ~~) So,~) +X20 t _o) r 0,S~ Io, o)+ rs ~,oX so)+ (C9a) (C9b) (C9c) (C9d) (ClOa) (ClOb) + S,(x0) [rsxo ~, s ) + sx(o2) rs, ~, sO)] (C10c)

170 O = S"(XO) r~(OrC- s 0) + S,(XO) J,,O,(X-O. s ) + F0@10Xo s 0) + + S, Gio)[F,, @Y, s 0) + S, @io) F(-Y, S)] (Cl0d) As a result, we obtain, &Sx- F) G(x-0, s0) ~ 1,212 where the function G(.,.) is composed of first-order partial derivatives of s(xO) and up-to fourth-order partials of Jj(,*,with respect to x1 and X2. Repeating this process for t(x), and equating, yields, rx1,(., W, O) + r'~(W0 s 0) Px1~@YI O + Px(! t) and 1 hxIs-ax' e hI(X" ail aii and aWS1a) a~t1a< where,

171 M =2,3, By continuing this process for higher order derivatives, we see that if, [r (x SO0) - C(x4O)] [p @Y 0) - 0) rx1~(x, Cx0O ) +r X2 (X02 ) PX i(-3 t0) + PX,1"2(X-0' t0) and I Mx,')a~' 1 (- M RZ rx1~(x, CX "s0) +r ~x2XO.) Pxx XOt0) + pI(X-0,t) then, and aKS (XO) tXO ai a~. DX& DXX.' DXK where, M =2, 3,. K = 2,3, 11,12 '2) 1. K = 1,2 Q% ETD

172 APPENDIX D: Proof of Lemma 3-3 We are given, 4i- [j0,O)- ~ x [p (x- t0) - p(-~ i1 1 r1(.0,O) + O 0) + P. t0) and (Dla) Ae [r @0 0) x ~ +liI. hi ax.... axbri (x03 t0) + p 0 t0) 1 (Dlb) where, M =2,3,. ip '29..gm= 1,2. From (Dila) we know, aDa. ~ rxx(X, so)r(so - [p, @x0 t0)+p- XO)] Da a r 1 (x 0, s 0) + a ~x) a 1 '0 t)+ t) ~ ~O)- f3~O)](D2b) rx1x1(xj0, s0) +,(2s0) px1j1(x0 t0) +,(!tO), Dividing the fight-hand side of (D2a) by the right-hand side of (D2b), and the left-hand side of (D2a) by the left-hand side of (D2b), we obtain,

173 a[r (x 0, SO) - ax0 a [P (X- ol to) - PCX-o)l c-k2 a lp (X-o 3, s 0) - 0 (X-O)l, axi (D3a) Dividing the right-hand side of (D2a) by the right-hand side of (Dl1 b), and the left-hand side of (D2a) by the left-hand side of (Dib), we obtain,, a ai &I (D3b) and, 'll 1..3 m= 1, 2. A similar division of (Dilb) by (D2b) yields', &a I... &axI. cb (D3c,) We define functions RM and PM, M - elp[r@-0, so) - cX-O)I P = ax.. a

174 for M = 2,3,....and 1~~o o)-(~o R1 so by assumption, a RM PM [r (xi, sO) - ~xx0 [p 1 ( so ) - P~O )] [r (xOs0)~ -O) ~ pxOo -fx~) where, Differentiating (D4a) with respect to x, yields, (D4a) (DOb) _ _ _ _ _ _ _ _R a- [r (x0, s 0) - a(xx0)] rxlx, (Xj0, so) pM PM a lp (xO s0) -P(3e)] Px1x1(x.0o SO)' (D5a) and differentiating (D4b) with respect toX2 yields, Rm RMPM PM a~[r(xO~so>..(x.x)] r,(@j0,S0) = 3 ~[p(-O~sO)...](XO)] p,(X!,s0) (D5b)

175 By assumption', xi xi RMX2 pM a [r(x0, s')~ -a(xi)] a II 09Oso) - so that (D5a) and (D5b) yield, Rm =-PM Da Rm =-PM Db Combining (D6a) and (D6b) yields, Rm PM (D7) r11(.0., (XO s) + cx!,O) Px1x1(x02 O) + P.1X1(-0' O)' QED

BIBLIOGRAPHY 176

177 [1] Abidi, M.A., 1989, Sensor fusion: a new approach and its applications, Sensor Fusion II: Human and Machine Strategies, P.S. Schenker, Ed., SPIE, Philadelphia, PA. [2] Andrews, H.C., and B.R. Hunt, 1977, Digital Image Restoration, Prentice-Hall, Inc., Englewood Cliffs, New Jersey. [3] Aschbacher, J., 1990, Microwave signatures from land surface radiometry, Intern. Geosci. Rem. Sens. Symp., College Park, MD, pp. 1149-1152. [4] Babaud, J., A. Witkin, M. Baudin, R. Duda, 1986, Uniqueness of the Gaussian kernel for scale-space filtering, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-8, No. 1, pp. 26-33. [5] Bergholm, F., 1987, Edge focussing, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-9, No. 6, pp. 726-741. [6] Berzins, V., 1984, Accuracy of laplacian edge detectors, Comp. Vision, Graphics, Image Process., Vol. 27, pp. 195-210. [7] Bischof, W., and T. Caelli, 1988, Parsing scale-space and spatial stability analysis, Comp. Vis. Graph. and Image Proc., Vol. 42, pp. 192-205. [8] Bixler, J.P. and D.P. Miller, 1987, A sensory input system for autonomoous mobile robots, Spatial Reasoning and Multi-Sensor Fusion, A. Kak and S.-S. Chan, Ed., Morgan Kaufman, Inc., Los Altos, CA. [9] Blanchard, B.J., and A.T.C. Chang, 1983, Estimation of soil moisture from Seasat SAR data, Water Res. Bull. 19, pp. 803-810. [10] Bracewell, R., 1986, The Fourier Transform and Its Application, McGraw-Hill, Inc., New York. [11] Brisco, B., F.T. Ulaby, and M.C. Dobson, 1983, Spaceborne SAR data for land-cover classification and change detection, 1983 IEEE Intern. Geosci. Rem. Sens. Symp., San Fransisco, CA. [12] Burke, W.J., T. Schmugge, and J.F. Paris, 1979, Comparison of 2.8- and 21-cm microwave radiometer observations over soils with emission model calculations, JGR 84, pp. 287-294. [13] Canny, J., 1986, A computational approach to edge detection, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-8, No. 6, pp. 679-698. [14] Camillo, P.J., and T.J. Schmugge, 1984, Correlating rainfall with remotely sensed microwave radiation using physically based models, IEEE Trans. on Geosc. and Rem. Sens. GE-22, pp. 415-423. [15] Chandrasekhar, S., 1960, Radiative Transfer, Dover, Inc., New York. [16] Chang, A.T.C., J.L Foster, and D.K. Hall, 1990, Satellite sensor estimates of northern hemisphere snow volume, Int. J. Remote Sensing, Vol. 11, No. 1, pp. 167-171. [17] Chen, J.S. and G. Medioni, 1989, Detection, localization, and estimation of edges, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-11, No. 2, pp. 191-198.

178 [18] Chin, R.T., Yeh, C., and Olson, W.S., 1985, Restoration of multichannel microwave radiometric images, IEEE Trans. Patt. Anal. Mach. Intell., PAMI-7, pp. 475-484. [19] Choudhury, B.J., C.J. Tucker, R.E. Golus, and W.W. Newcomb, 1987, Monitoring vegetation using Nimbus-7 scanning multichannel microwave radiometer data, Int. J. Remote Sensing, Vol. 8, No. 3, pp. 533-538. [20] Clark, J.J., 1989, Authenticating edges produced by zero-crossing algorithms, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-11, No. 1, pp. 43-57. [21] De Micheli, E., B Caprile, P Ottonello, and V. Torre, 1989, Localization and noise in edge detection, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-11, No. 10, pp. 1106-1117. [22] Dudgeon, D.E., and R.M. Mersereau, 1984, Multidimensional Digital Signal Processing, Prentice-Hall, Inc., Englewood Cliffs, New Jersey. [23] Durrant-Whyte, H.F., 1987, Sensor models and multi-sensor integration, Spatial Reasoning and Multi-Sensor Fusion, A. Kak and S.-S. Chan, Ed., Morgan Kaufman, Inc., Los Altos, CA. [24] Edgerton, A.T., A. Stogryn, and G. Poe, 1971, Microwave Radiometric Investigations of Snowpacks, Final Rept. 1285R-4 of Contract 14-08-001-11828 between Aerojet-General Corp., El Monte, CA, and the U.S. Geological Survey. [25] England, A.E., 1974, Thermal microwave emission from a halfspace containing scatterers, Radio Science, Vol. 9, No. 4, pp. 447-454. [26] England, A.W., 1974, The effect upon microwave emissivity of volume scattering in snow, in ice, and in frozen soil, Proc. URSI Spec Mtg on Microwave Scattering and Emission from the Earth, Berne, Switzerland, 23-26 Sept., 1974. [27] England, A.W., 1975, Thermal microwave emission from a scattering layer, JGR 80, pp. 4484-4496. [28] England, A.W., 1976, Relative influence upon microwave emissivity of fine-scale stratigraphy, internal scattering, and dielectric properties, Pageoph 114, pp. 287-299. [29] England, A.W., 1977, Microwave brightness spectra of layered media, Geophysics 42, pp. 514-521. [30] England, A.W., 1990, Radiobrightness of diurnally heated, freezing soil, IEEE Trans. Geosc. and Rem. Sens., Vol. GE-28, No. 4, pp. 464-476. [31] England, A.W., M.C. Dobson, and F.T.Ulaby, 1990, "Mapping regional freeze/thaw patterns with satellite microwave radiometry," Continuation proposal for NASA Grant NAGW-1983, NASA Hydrology Program. [32] Ferraro, R.R., N.C. Grody, and J.A. Kogut, 1986, Classification of geophysical parameters using passive microwave satellite measurements, IEEE Trans. Geosci. Rem. Sens., GE-24, No. 6, pp. 1008-1013. [33] Gallo, K.P. and J.F. Brown, 1990, Satellite-derived indices for monitoring global phytoclimatology, Intern. Geosci. Rem. Sens. Symp., College Park, MD, pp. 261-264.

179 [34] Geiger, D., and T. Poggio, 1987, Level crossings and the panum area, Proc. Workshop on Computer Vision, Miami, FL, pp. 211-214. [35] Gloerson, P. and F.T. Barath, 1977, A scanning multichannel microwave radiometer for Nimbus-G and SeaSat-A, IEEE J. Oceanic Engin., OE-2, No. 2, pp. 172-178. [36] Goffman, C., 1965, Calculus of Several Variables, Harper & Row, New York. [37] Graustein, W.C., 1966, Differential Geometry, Dover Publishing, Inc., New York. [38] Grody, N.C., 1988, Surface identification using satellite microwave radiometers, IEEE Trans. Geosci. Rem. Sens., GE-26, No. 6, pp. 850-859. [39] Hoehstra, P. and Delaney, A., 1974, Dielectric properties mof soils at UHF and microwave frequencies, J. Geophys. Res. 79, pp. 1699-1708. [40] Hollinger, J., R. Lo, G. Poe, R. Savage, and J. Pierce, 1987, Special Sensor Microwave/Imager User's Guide, Naval Research Laboratory, Washington, DC. [41] Hummel R. and R. Moniot, 1989, Reconstructions from zero crossings in scale space, IEEE Trans. Acous. Speech Sig. Process., Vol. ASSP-37, No. 12, pp. 2111-2130. [42] Jackson, T.J. and T.J. Schmugge, 1989, Passive microwave remote sensing system for soil moisture: some supporting research, IEEE Trans. Geosci. Rem. Sens., GE-27, No. 2, pp. 225-235. [43] Kak, A. and S.-S. Chan, 1987, Spatial Reasoning and Multi-Sensor Fusion, Morgan Kaufman, Inc., Los Altos, CA. [44] Kass, M., A. Witkin, and D. Terzopoulos, 1988, "Snakes: active contour models," Intern. J. Comp. Vision, pp. 321-331. [45] Koenderink, J., 1989, A hitherto unnoticed singularity of scale-space, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-11, No. 11, pp. 1222-1224. [46] Kunzi, F.K., S. Patil, and H. Rott, 1982, Snow-cover parameters retreived from Nimbus-7 scanning multichannel microwave radiometer (SMMR) data, IEEE Trans. Geosci. and Remote Sensing, GE-20, No. 4, pp. 452-467. [47] Lagendijk, R., J. Biemond, and R. Boekee, "Hierarchical Blur Identification," Proc. Inter. Conf. Acoustics, Speech, Signal Processing, Albuquerque, NM, 1990, pp. 1889-1892. [48] Lee, T., J.A. Richards, and P.H. Phillips, 1987, Probabalistic and evidential approaches for multiscource data analysis, IEEE Trans. Geosci. Rem. Sens., GE-25, No. 3, pp. 283-293. [49] Lindeberg, T., 1990, Scale-space for discrete signals, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-12, No. 3, pp. 234-254. [50] Liu, L., B. Schunck, and C. Meyer, 1990, On robust edge detection," EECS Technical Report, University of Michigan. [51] Lu, Y., and R. Jain, 1989, Behavior of edges in scale space, IEEE Trans. Pantt. Anal. Mach. Intell., Vol. PAMI- 11, No. 4, pp. 337-356.

180 [52] Marr, D., and E. Hildreth, 1980, Theory of edge detection, Proc. R. Soc. London B, Vol. 207, pp. 187-217. [53] Mammone, R., "Image restoration using linear programming," Image Recovery: Theory and Application, H. Stark, Ed., Academic Press, Inc., Orlando, FL, 1987, pp. 127-156. [54] Moik, J., 1980, Digital Processing of Remotely Sensed Images, NASA, NASA SP-431. [55] NASA, 1978, The Nimbus-7 Users Guide, The Landsat/Nimbus Project, Goddard Space Flight Center, NASA. [56] Njoku, E.G., E.J. Christensen, and R.E. Cofield, 1980, The Seasat scanning multichannel microwave radiometer (SMMR): Antenna pattern corrections -- development and implementation, IEEE Trans. Ocean Engin. Vol. OE-5, No. 2, pp. 125-137. [57] Paloscia, S., P Pampaloni, L. Chiarantini, P. Coppo, S. Gagliani, and G. Luzi, 1990, Multifrequency microwave radiometric measurements of soil moisture, Intern. Geosci. Rem. Sens. Symp., College Park, MD, pp. 1837-1840. [58] Richards, J., 1986, Remote Sensing Digital Image Analysis, Springer-Verlag, Berlin. [59] Samadani, R., D. Mihovilovic, C.R. Clauer, G. Wiederhold, J.D. Craven, and L.A. Frank, 1990, "Evaluation of an elastic curve technique for finding the auroral oval from satellite images automatically," IEEE Trans. Geosci. Remote Sensing, Vol. GE-28, No. 4, pp. 590-597. [60] Saund, E., 1990, Symbolic construction of a 2-D scale-space image, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-12, No. 8, pp. 817-830. [61] Schenker, P.S., 1989, Sensor Fusion II: Human and Machine Strategies, SPIE, Philadelphia, PA. [62] Schmugge, T.J., 1983, Remote sensing of soil moisture: Recent advances, IEEE Trans. on Geosc. and Rem. Sens. GE-21, pp. 336-344. [63] Schmugge, T.J., 1987, Remote sensing applications in hydrology, Rev. Geophys. 25, pp. 148-152. [64] Schmugge, T.J., P.E. O'Neill, and J.R. Wang, 1989, Passive microwave soil moisture research, IEEE Trans. Geosci. Rem. Sens., GE-24, No. 1, pp. 12-22. [65] Schunck, B., 1986, Gaussian Filters and Edge Detection, Research Publication, General Motors Res. Lab., Warren, MI, Tech. Rep. GMR-5586. [66] Session TP-10, 1990, Multisensor data fusion/GIS, Intern. Geosci. Rem. Sens. Symp., College Park, MD. [67] Soliday, B., 1989, Identification of data redundancies in vegetation index models, Proc. Int. Geosci. Remote Sensing Symp., Vancouver, B.C., pp. 1339-1342. [68] Spivak, M., 1965, Calculus on Manifolds, Addison-Wesley Publ. Co., Reading, MA.

181 [69] Srinivasan, A. and J.A.Richards, 1990, Knowledge-based techniques for multi-source classification, Int. J. Remote Sensing, Vol. 11, No. 3, pp. 505-525. [70] Stark, H., and J.W. Woods, 1986, Probability, Random Processes, and Estimation Theory for Engineers, Prentice-Hall, Englewood Cliffs, NJ. [71] Stephanou, H.E. and A.M. Erkmen, 1989, Multiresolution sensor fusion by conductivity analysis, Sensor Fusion II: Human and Machine Strategies, P.S. Schenker, Ed., SPIE, Philadelphia, PA. [72] Torre, V., and T. Poggio, 1986, On edge detection, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-8, No. 2, pp. 147-163. [73] Tsang, L., J.A. Kong, and R.T. Shin, 1985, Theory of Microwave Remote Sensing, John Wiley & Sons, Inc., New York. [74] Ulaby, F.T., R.K. Moore, and A.K. Fung, 1981, Microwave Remote Sensing: Active and Passive, Addison-Wesley, Reading, MA. [75] Van Warmerdam, W., and V. Algazi, 1989, Describing 1-D intensity transitions with Gaussian derivatives at the resolution matching the transition widths, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-11, No. 9, pp. 973-977. [76] Witkin, A., "Scale-space filtering," Proc. Int. Joint Conf. Artif. Intell., Karsruhe, West Germany, 1983, pp. 1019-1021. [77] Witkin, A., "Scale space filtering: A new approach to multiscale description", Image Understanding 1984, S. Ullman and W. Richards, Ed., Ablex Publ. Co., Norwood, NJ, 1984. [78] Witkin, A., D. Terzopoulos, and M. Kass, 1988, Signal matching through scale space, Int. J. Computer Vision, Vol. 1, pp. 134-144. [79] Wang, J.R., T.J. Schmugge, W.I. Gould, W.S. Glazar, and J.E. Fuchs, 1982, A multi-frequency radiometric measurement of soil moisture content over bare and vegetated fields, Geophys. Res. Let. 9, p. 416-419. [80] Weaver, C.B., 1989, Sensor Fusion II, SPIE, Orlando, FL. [81] Wu, L. and Z. Xie, 1990, Scaling theorems for zero-crossings, IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-12, No. 1, pp. 46-54. [82] Young, M., 1986, Optics and Lasers Including Fibers and Integrated Optics, Springer-Verlag, Berlin. [83] Yuille, A., and T. Poggio, 1985, Fingerprint theorems for zero crossings, J. Opt. Soc. Amer. A, Vol. 2, No. 5, pp. 683-692. [84] Yuille, A. and T. Poggio, "Scaling theorems for zero crossings," IEEE Trans. Patt. Anal. Mach. Intell., Vol. PAMI-8, No. 1, Jan. 1986, pp. 15-25. [85] Yuille, A. and T. Poggio, "Scaling and fingerprint theorems for zero crossings", Advances in Computer Vision, Vol. 2, C. Brown, Ed., Lawrence Erlbaum Associates, Hillsdale, NJ, pp. 47-79, 1988.

182 [86] Jain, A.K., 1986, Fundamentals of Digital Image Processing, Prentice-Hall, Inc., Englewood Cliffs, NJ. [87] England, A.W. and B. Zuerndorfer, 1990, "Mapping temporal freeze/thaw patterns in Alaska with ERS-1," Proposal for NASA Grant NRA 90-OSSA-17, Polar Oceans, Land, and Ice Cover Using Data From the Alaska SAR Facility. [88] Hollinger, J., R. Lo, G. Poe, R. Savage, and J. Pierce, 1987, Special Sensor Microwave/Imager User's Guide, Naval Research Laboratory, Washington, D.C..

ABSTRACT SPATIAL AND SPECTRAL ANALYSES OF REMOTELY SENSED IMAGES USING SCALE-SPACE TECHNIQUES by Brian Zuerndorfer Co-Chairs: Gregory H. Wakefield, Anthony W. England The thesis concerns multiple sensor processing and the development of a unified representation for multiple-sensor data. Considered are the images of land surfaces generated at different (center) frequencies from a single satellite. Through comparisons of these multispectral images, surfaces are segmented into regions, and regions are classified. All frequency channels share a common aperture, so that the images generated at different frequencies are of different resolution. To avoid mis-classification, images at each frequency are processed to compensate for resolution differences prior to classification. With no a priori scene information available, resolution compensation is performed by synthesizing images at all frequencies at the resolution of the lowest frequency (coarsest) image used in classification. Thus, fine-resolution information is lost. The thesis presents a technique to recover fine-resolution information in surface classification. Specifically, if the spatial transfer function of the sensor is approximately Gaussian, then scale-space techniques of computer vision can be used to estimate boundaries between different surfaces at fine-resolution. First, surface scenes are segmented and classified using multispectral processing at coarse-resolution. Second, boundaries between different surfaces, as derived from multispectral processing, are registered to boundaries on high-frequency images which have been resolution 1

2 compensated (i.e., coarse-resolution images derived from high-frequency, fine-resolution images). Third, boundaries on compensated, high-frequency images are tracked from coarse-to-fine resolution as resolution compensation is reduced. Estimates of surface boundaries are the resulting boundary locations in uncompensated, fine-resolution images.