RL923 ELECTROMAGNETIC WAVE SCATTERING BY POWER-LAW SURFACES RICHARD TILLMAN AUSTIN 1994 RL-923 = RL-923

ELECTROMAGNETIC WAVE SCATTERING BY POWER-LAW SURFACES by Richard Tillman Austin A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1994 Doctoral Committee: Professor Anthony W. England, Co-Chairman Professor Fawwaz T. Ulaby, Co-Chairman Assistant Professor Kamal Sarabandi Professor John F. Vesecky

T ~ Richard Tillman Austin 1994 All Rights Reserved

To my parents, Bob and Betty, and my sister, Michelle ii

ACKNOWLEDGEMENTS I thank all the members of my committee for their support and helpful suggestions. Special thanks go to my advisor, Prof. Tony England, who has encouraged me in my work and given me a taste of professorial life by allowing me to lecture, TA, and collaborate in writing. Special thanks also go to Prof. Fawwaz Ulaby, who has always shown great confidence in me and from whom I have learned much about remote sensing, both in the classroom and in the course of our research projects. Many others have generously contributed their time and expertise in helping with my work. I thank Prof. Greg Wakefield for his help with the spectral estimation and filtering problems that are a central part of this work. I thank Mr. Rajendra Aggarwala of the Department of Civil and Environmental Engineering for the use of his surveying equipment at Mount St. Helens, and Mr. Bruce Babb and the staff of the Mount St. Helens National Volcanic Monument for their cooperation. Special thanks go to Mr. Ron Hartikka for his assistance and advice in the laboratory. I also thank Mr. Steve Thiel, Mr. David Carter, and Mr. D. Alexander van Oosten Slingeland for their help in my research as part of their directed study projects. In my years at the Radiation Laboratory, I have come to know many students who have helped me learn difficult concepts, helped me with experiments, helped with programming, and helped me to relax when my attitude faltered. I thank Dr. Kyle McDonald, Dr. Leland Pierce, Dr. Emilie van Deventer, Dr. Jeff Collins, Dr. Mike Whitt, Dr. Brian Zuerndorfer, Mr. Jim Stiles, Mr. John Galantowicz, Mr. Joe Landry, iii

Mr. John Kendra, Mr. Adib Nashashibi, Mr. Ed Kim, and Mr. Roger DeRoo for their technical, theoretical, and experimental assistance. I also thank them and my other friends in the Radiation Laboratory for their encouragement through a long and difficult process. I thank the National Science Foundation, the National Aeronautics and Space Administration, the Michigan Space Grant program, the Rackham School of Graduate Studies, and the Department of Electrical Engineering and Computer Science for financial support. Finally, I thank my parents, Bob and Betty Austin, and my sister, Michelle Wilford, for their unlimited support, love, patience, and encouragement during the last eight years. iv

TABLE OF CONTENTS DEDICATION................................. ACKNOWLEDGEMENTS.......................... LIST OF FIGURES............................... ii iii 111 Vlll. LIST OF TABLES...... LIST OF APPENDICES.. CHAPTER I. INTRODUCTION....................... X ll xiv 1.1 1.2 1.3 Interest in Power-Law Surfaces......... Questions to be Answered in This Dissertation Format of This Dissertation............ 1 4 5 II. STATISTICAL DESCRIPTION OF ROUGH SURFACES 2.1 Gaussian Surfaces....................... 2.2 Power-Law Surfaces....................... III. MEASUREMENTS OF SURFACE ROUGHNESS...... 3.1 Mount St. Helens........................ 3.2 2D Laser Profilometer System.................. 3.3 Profilometer Data Processing.................. 3.3.1 Overheating Errors.................. 3.3.2 Level Shift Errors................... 3.4 Surveying Measurements.................... 3.5 Measured Profiles........................ IV. SPECTRAL ESTIMATION.................... 4.1 Classical Spectral Estimators.................. 7 9 10 16 16 20 25 26 35 39 41 42 45 v

4.2 Estimation of Power-Law Spectra................ 48 4.2.1 Fourier-based estimators............... 48 4.2.2 Pre-whitening..................... 54 4.2.3 Capon's estimator................... 56 4.2.4 Variance comparison................. 59 4.2.5 Formulas for real frequency.............. 61 4.2.6 Aliasing......................... 63 4.2.7 Applicability to real data............... 64 4.3 Estimation of Two-Dimensional Spectra.................64 4.3.1 Two-dimensional estimators................. 65 4.3.2 2D spectral estimation of isotropic surfaces..... 65 4.4 Roughness Spectra of Mount St. Helens Debris Flows........ 70 4.4.1 Verification of median interpolation..............70 4.4.2 Primary debris flow surfaces..................71 4.4.3 EDM noise floor........................ 76 4.4.4 Mudflow- or water-deposited surfaces............ 76 4.5 Summary................................. 79 V. DESIGN AND MANUFACTURE OF ARTIFICIAL POWERLAW SURFACES................................ 80 5.1 Analogue design history.......................... 81 5.2 Dielectric materials........................... 82 5.3 Surface Design Algorithm.........................85 5.4 Milling System............................. 88 5.4.1 Original Milling System................ 88 5.4.2 Failure of the Techno System............. 91 5.4.3 Upgraded Milling System.................. 93 5.5 Milling procedure............................. 95 5.5.1 Generation of Toolpath Files................... 95 5.5.2 Milling machine operation...................97 5.6 Surface Completion and Temperature Sensitivity..........100 5.7 Surface Verification....................... 102 VI. RADAR MEASUREMENTS.................... 105 6.1 35 GHz Radar.......................... 105 6.2 Measurement Configuration.....................108 6.3 Calibration............................ 111 6.4 Measurement Procedure..................... 114 6.5 Post-Processing and Data Reduction.............. 115 6.5.1 Scattering Coefficients................. 116 6.5.2 Independent Samples................. 118 vi

VII. ROUGH SURFACE SCATTERING MODELS AND COMPARISONS TO EXPERIMENT................... 121 7.1 7.2 7.3 7.4 7.5 Experimental Results.......................... 121 Comparison to Physical Optics Model............. 122 Comparison to Geometric Optics Model............... 125 Comparison to Small Perturbation Model............ 128 Comparison to PWH model....................... 130 VIII. CONCLUSIONS, CONTRIBUTIONS, AND RECOMMENDATIONS............................... 135 8.1 Conclusions of This Dissertation................. 135 8.2 Comments on Methods and Theory............. 137 8.3 Contributions of This Dissertation............... 138 8.4 Recommendations for Future Work................ 139 APPENDICES................................... 141 BIBLIOGRAPHY...............................214 vii

LIST OF FIGURES Figure 1.1 Variation of normalized power spectral density G(1/A)/k with wavelength A. (Figure reproduced from Sayles and Thomas [52].).... 3 2.1 A power-law roughness spectrum in the form of (2.8)........... 11 2.2 Examples of Gaussian and power-law surfaces............ 13 3.1 Map of Mount St. Helens showing debris flows and measurement sites 18 3.2 Two dominant surface types observed on the debris flows...... 19 3.3 The 2D Laser Profilometer System.................... 21 3.4 Circuit diagram of external serial port switch............ 23 3.5 One of the rougher surfaces scanned by the EDM.......... 25 3.6 Sample surface elevation grid with outlier............... 26 3.7 Bimodal distribution of median residuals............... 27 3.8 Histogram of actual median residuals................. 28 3.9 Elevation distributions for 5 x 5 subregions on two surfaces, one of which has an outlier........................... 29 3.10 Successive stages of the median residual and lower quartile difference filtering procedure for scan 3a...................... 32 3.10 (continued).................................... 33 3.11 Median interpolation for scan 3a................. 36 viii * * Vlll

3.12 A profilometer scan showing the level shifts. (The EDM did not move during this scan.)............................. 37 3.13 Top view of survey configuration............................ 40 4.1 Spectral domain window functions for the periodogram and for the modified periodogram with Hanning window............. 46 4.2 Expected value of the periodogram spectral estimator for two powerlaw surfaces............................... 50 4.3 Expected value of PPER for five 256-point power-law profiles evaluated at frequencies corresponding to the case of no zero-padding... 51 4.4 Exact values of the power-law spectra used in calculations of the expected values of the spectral estimates............... 51 4.5 Expected values of PHAN(fA) for five different 256-point power-law profiles.................................. 52 4.6 Expected values of Ppw(fa) for five different 257-point power-law profiles.................................. 56 4.7 Expected values of PCAP(fA) for synthetic profiles with various /d using an autocorrelation matrix of dimension p = 70........ 58 4.8 Capon's, modified periodogram, and pre-whitened periodogram spectral estimates for ten different power-law profiles............... 60 4.9 Aliasing in the power-law spectrum case.................. 63 4.10 41 x 41-point sample of a 2048 x 2048-point synthetic power-law surface 69 4.11 Averaged Capon's estimates from 2D synthetic surface height data at three scales with power-law fit.................... 69 4.12 Spectra of primary debris flow surfaces................ 74 4.13 EDM noise floor............................. 77 4.14 Spectra of smooth surface scans..................... 78 4.15 Survey spectra.............................. 79 ix

5.1 Configuration for measuring UHMW dielectric constant....... 84 5.2 The original milling system supplied by Techno............ 89 5.3 A load-bearing spring was installed on the z-axis to help support its w eight................................... 92 5.4 The upgraded mill configuration....................... 93 5.5 Edges of the artificial surface are rounded by subtracting the surface elevation from a rounded reference surface with a 1/4-inch radius of curvature................................. 96 5.6 The spacing between adjacent finishing rows (0.0137 inch) was slightly less than the radius of the 1/32-inch ball-end endmill (0.0156 inches). 99 5.7 Grayscale file of estimated error in UHMW1............. 103 5.8 A histogram of differences between the specified surface elevations and the scanned surface height.........................103 6.1 Block diagram of the 35 GHz radar................... 107 6.2 Components of the radar measurement system...............109 6.3 Schematic of measurement setup.................... 110 6.4 Test measurement of calibration sphere (1.75 inches in diameter).. 113 6.5 Test measurement of 1.188-inch diameter sphere............113 6.6 Arrangement of footprints on slab................... 119 7.1 Measured co-polarized backscattering coefficients, ao and chh, of the three artificial surfaces....................... 122 7.2 Measured cross-polarized backscattering coefficients, a~Vh and Ao, of the three artificial surfaces....................... 123 7.3 Physical optics results for (a) HH and (b) VV backscatter compared to the measured values for the three artificial surfaces...........126 7.4 Geometrical optics results for backscattering compared to the measured values for the three artificial surfaces............... 128 x

7.5 Small perturbation model results for (a) HH and (b) VV backscatter compared to the measured values....................... 131 7.6 Phased Wiener-Hermite model results for (a) HH and (b) VV backscatter for the roughest and smoothest surfaces compared to the measured values............................... 134 A.1 Grayscale plots of profiles at scan site 2................... 147 A.2 Grayscale plots of profiles at scan site 3................... 148 A.3 Grayscale plots of profiles at scan site 4................... 149 A.4 Grayscale plots of profiles at scan site 5................ 150 A.5 Grayscale plots of profiles at scan site 6................ 151 A.6 Grayscale plots of the surveys..................... 152 B.1 Coordinate system of transmitter polarizers.............. 154 C.1 Geometry of PWH model........................ 168 xi

LIST OF TABLES Table 2.1 Values of kl and ka for the artificial power-law surfaces described in Table 5.1, calculated using three different methods......... 14 3.1 Rows excluded from spectral estimates due to visible level shifts..38 4.1 Spectral statistics derived from fits to modified periodograms, periodograms based on pre-whitened data, and Capon's estimators of ten synthetic 64-point profiles...................... 61 4.2 Power-law parameters derived from fits to averaged Capon's estimates at three scales........................... 70 4.3 Verification of median interpolation.................... 71 4.4 Spectral estimation parameters for the profilometer and survey data sets............................ 73 5.1 Spectral parameters (a, y, a, and a in (5.8)) used and fractal dimension Df of the surface analogues.................. 88 5.2 Types and sizes (in MB) of Mastercam files............. 97 6.1 Estimated number of independent samples............... 120 A.1 Filtering and interpolation steps for scan 2b (Johnston Ridge site).. 142 A.2 Filtering and interpolation steps for scan 2c (Johnston Ridge site).. 143 A.3 Filtering and interpolation steps for scan 3a (Johnston Ridge site).. 143 A.4 Filtering and interpolation steps for scan 3b (Johnston Ridge site).. 143 A.5 Filtering and interpolation steps for scan 3c (Johnston Ridge site).. 144 xii

A.6 Filtering and interpolation steps for scan 4a (Johnston Ridge site).. 144 A.7 Filtering and interpolation steps for scan 4c (Johnston Ridge site).. 145 B.1 Values of t and Et for various transmit polarizations.............156 C. 1 Comparison of moments of WH and perturbation expansions of eik((x) with moments of exact expression................... 166 xiii

LIST OF APPENDICES Appendix A. MOUNT ST. HELENS DATA........................142 B. COHERENT-ON-RECEIVE RADAR CALIBRATION........ 153 C. THE PHASED WIENER-HERMITE MODEL...............162 C.1 General Outline of PWH Formulation................... 167 C.2 Derivation of the Backscattering Coefficient a~............. 181 C.3 Evaluation of (C.111) for the Power-Law Case...........186 D. FLUCTUATION STATISTICS OF MILLIMETER-WAVE SCATTERING FROM DISTRIBUTED TARGETS..................189 E. MILLIMETER-WAVE RADAR SCATTERING FROM SNOW.... 204 xiv

CHAPTER I INTRODUCTION Rough surface scattering has long been a topic of interest in remote sensing. Improved models of electromagnetic wave scattering at rough surface interfaces benefit both investigators who study natural surfaces (e.g., geologists and volcanologists) and those for whom surface scattering effects are an undesirable corruption of signals from other targets above or below the surface. Measurements of backscatter from forest canopies, for instance, must make assumptions about the surface scattering properties before the canopy scatter can be determined. If the surface scattering is modeled incorrectly, interpretation of canopy scatter can be adversely affected. 1.1 Interest in Power-Law Surfaces Models of scattering by rough surfaces typically model surface elevation as a random process of position having some distribution of surface height and some correlation function representing the dependence between surface elevations at two surface locations. Over the last few decades, a large majority of studies of scattering by land surfaces have assumed that surfaces possess a Gaussian correlation function (and corresponding Gaussian roughness spectrum). Use of non-Gaussian surface models has been increasing, however. In particular, numerous studies have shown 1

2 that many natural surfaces are better described by power-law roughness spectra of the form Sz(f) = clfl-, (1.1) where f is the spatial frequency, and c and /3 are constants, with 1 <, < 3. Powerlaw surfaces are scaling; they may be described using a fractal dimension Df, which varies between 2.0 and 3.0 for a surface, or between 1.0 and 2.0 for a linear surface profile. (The relation between Df for a linear profile and /3 is given by (2.10).) Examples of natural surfaces that are well-described by power-law spectra include land surfaces [21, 30, 31, 52], rock surfaces [10], and seafloor topography [6, 7, 22, 26]. Ocean surfaces may also be described by a power-law spectrum having a form related to (1.1) [25]. Although synthetic topographies appear most realistic with Df 2.2 [34, 60], yielding linear profiles of Df e 1.2 and,/ e 2.6, several studies [6, 7, 30, 31, 52] report measured values of #/ near 2.0. For example, Sayles and Thomas [52] show measured roughness spectra for 23 surfaces ranging from lunar terrain to the surface of a hip joint (reproduced here as Figure 1.1). Sayles and Thomas conclude that natural topography has a characteristic power-law exponent of /3 = 2.0; this claim is also made by Huang and Turcotte [30]. If true, the fixed value of f3 would allow rough surface characterization using a single parameter, c in (1.1). Two current research projects in the Radiation Laboratory are studies of scattering by volcanic terrains and the diffuse radar backscatter from Mars. Both terrestrial volcanic terrains and Martian terrains (formed largely through volcanism and aeolian processes) seem likely candidates for modeling by a power-law roughness spectrum. The interest in these specific surfaces led to the study of power-law surfaces described in this dissertation.

3 105 1- 0 -' --- --- --- --- -- ' - --- --- -- E*.. _ _ _ _, 10' E -. 1. -2 4+. o 0 7- 10 -7 10 — Key Surface k(m) ref. A Ptriphety of honed berng roctway tOx10-12 '108 I - Periphery of ground disc 1tx-O'12 Part periphery of ground disc 30xT1016 4 Periphery of ground bearig rmceay K)Oxf13 Z3 O Pftiprwy o ground disc '7xlO1 I 10 - - I Fiphery of ring-lapped stil ball 5. 70- - Part pe phery oa ring-lapped stal blal 1 0 xl13 - SuIoca ground lacross -icy ) 30x11" -to Plunge groud taring Facwaybawss-lay) 20xlO13 10 ~~ -. - /-r 10I2 L I I - I 0X)' 10'5 104 1 0'3 10'2 10 0 10I 102 103 Wavelength (.Am) Figure 1.1: Variation of normalized power spectral density G(1/A)/k with wavelength A. (Figure reproduced from Sayles and Thomas [52].)

4 1.2 Questions to be Answered in This Dissertation The studies described in this dissertation were designed to answer five principal questions: 1. Does the power-law roughness model apply to targets of interest in remote sensing? 2. If the power-law model applies to volcanic terrain, what are typical values of the power-law roughness parameters? 3. Given a power-law roughness model, how does the backscattering coefficient vary with incidence angle as a function of surface parameters? 4. How well do various rough surface scattering models predict scattering by a power-law surface? 5. In particular, is Eftimiu's Phased Wiener-Hermite model [16] useful for powerlaw surfaces? The first question has been partially answered by the cited sources; the interest here is in verifying whether volcanic terrains belong to the class of surfaces having power-law roughness spectra. The parameter values sought in the second question drove the design of synthetic power-law surfaces used in laboratory scattering measurements; these measurements provided answers to the third question. The fourth question was approached by applying the more common rough surface scattering models to the cases represented by the manufactured surface analogues. A recent scattering model by Eftimiu, the Phased Wiener-Hermite (PWH) model, was examined in some depth (question 5). The PWH model uses a novel approach

5 that seemed promising for modeling scattering by extremely rough dielectric surfaces; however, there were no published comparisons of experimental results to the PWH model at the outset of this study. 1.3 Format of This Dissertation Chapter II discusses the random process representation of a rough surface and compares Gaussian and power-law surface properties. Potential ambiguities associated with the description of power-law surfaces by use of Gaussian parameters are discussed. The field work at Mount St. Helens is described in Chapter III. The laser profilometer system used for collecting small-scale surface profiles is discussed, along with descriptions of the filtering and interpolation processes made necessary by problems with the profilometer. Large-scale roughness was surveyed; the surveying techniques are described here. A summary of the measured surface profiles is given in Appendix A. Chapter IV discusses spectral estimation and the special problems encountered in obtaining estimates of the surface roughness spectrum of a power-law surface. The chapter includes a description of a spectral estimator designed for use with limited data sets and shows the application of this estimator on the surface profiles from Mount St. Helens. The design and manufacture of artificial power-law surfaces is the subject of Chapter V. The design algorithm, selection of materials, and milling procedure are discussed, along with a summary of difficulties encountered during milling operations. Chapter VI describes the scattering measurements, in which backscattering from the surface analogues was measured using a 35 GHz scatterometer. The radar com

6 ponents, measurement configuration, calibration procedure, and measurement procedures are detailed here. The chapter concludes with a description of the data reduction process, in which the measured fields are converted into estimates of the backscattering coefficient. A brief derivation of the calibration technique is included in Appendix B. Measured values of oa~ are compared with those predicted by rough surface scattering models in Chapter VII. Results using classical models are presented, and results derived from the PWH model are also given. A brief derivation and description of the PWH model is given in Appendix C. Chapter VIII presents the conclusions of the dissertation, applications of the results, and recommendations for future work. The study of scattering by power-law surfaces (the subject of this dissertation) is the third project in the remote sensing of natural terrains undertaken by the author. The first study, Fluctuation Statistics of Millimeter-Wave Scattering from Distributed Targets, was completed in 1987. Results of the study were reported by Ulaby, Haddock, and Austin [57] in 1988; this paper is included here as Appendix D. The author participated in a second study, Millimeter-Wave Radar Scattering from Snow, in 1988-89. Results of this study were reported in 1991 by Ulaby, Haddock, Austin, and Kuga [58]; this paper is reproduced as Appendix E.

CHAPTER II STATISTICAL DESCRIPTION OF ROUGH SURFACES Most rough surface scattering models describe natural rough surfaces as realizations of random processes. Natural rough surfaces often have a random appearance, making the random process description plausible even if the underlying processes that produce the terrain are not random. While a random process description does not readily allow for surfaces with overhangs (because the surface elevation becomes a multiple-valued function), the majority of rough surface scattering studies do not consider such complex surfaces. The rough surface elevation process Z(x,y) is typically described in terms of its surface height probability function (pdf), fz(z), its mean function or ensemble average, Jz(x, y), its correlation function, Rz(xl, yi; x2, y2), and its surface roughness spectrum, Sz(f, fy), also called the power spectral density or PSD. Rough surface random processes are usually assumed to be stationary; i.e., the surface statistics are independent of position. (The correlation function, for instance, would depend on the distance between two points and not on their absolute positions.) An additional and usually necessary assumption is that the surface process is ergodic, meaning that ensemble averages may be replaced by spatial averages. 7

8 The surface height pdf, fz(z), is nearly always assumed to be Gaussian. Estimation of the pdf is extremely difficult for rough surface random processes (the interesting ones, at least) due to the very high quantity of data required. The mean value of the random process,,z(x, y), which is more commonly written E[Z(x,y)] or (Z(x,y)), E[Z(x,y)]= J zfz(z)dz, (2.1) -00 is usually set to zero by choosing a proper coordinate system. For a zero-mean surface, the covariance function, Kz, Kz(xl, yi; x2, Y2) = E[(Z(x, yi) - Iz(Xl, y1))(Z(x2, Y2) - xZ(x2, Y2))], (2.2) reduces to the identical form of the correlation function, Rz: Rz(X, yl; X2, Y2) = E[Z(x1, Y1)Z(x2, y2)], (2.3) so we will often use the terms covariance and correlation interchangeably. The covariance of a surface point with itself is called the variance and is usually denoted by o2: '2 = Rz(Xi, y; x, i) = E[Z2(x1, y)]. (2.4) As mentioned above, for stationary surfaces, the correlation function is a function of the offset between two points and not their absolute locations. If rx = x2- xl and Ty = Y2 - Y1, then Rz(rx, ry) = E[Z(x, y)Z(x + r-, y + rT)]. (2.5) The correlation function, R(,r, ry), is a measure of the correlation between surface elevations separated by a given distance. The Fourier transform of the correlation function is the surface roughness spectrum, Sz(f, fy), which is a measure of the

9 spectral content of the surface roughness as a function of spatial frequency. (Unless noted otherwise, roughness spectra will be shown as a function of spatial frequency f in m-1 rather than spatial wavenumber k = 27r/A, which has units of rad/m.) We shall see in later chapters that the roughness spectrum is often of interest because it allows comparisons between the spatial wavelength of surface roughness components and the electromagnetic wavelength of a remote sensing radar. 2.1 Gaussian Surfaces In this dissertation, the term Gaussian surface denotes a surface random process having a Gaussian correlation function. For surfaces with isotropic surface statistics (the surfaces considered in this dissertation all share this property), the correlation function has the form of a Gaussian:. Rz(r) = a2exp(-r2/12), (2.6) where r is the radial lag (r2 = r2 +Tr2), a2 is the surface variance, and I is a parameter called the correlation length. For lag r = 1, the correlation function drops to a value equal to e-1 times its maximum. Thus the correlation length becomes a figure of merit; surfaces with small I decorrelate over short distances. A two-dimensional Fourier transform gives the roughness spectrum for a Gaussian surface: Sz(fr) =2127r exp(-7Ir212f2) (2.7) where f- is the radial spatial f the radial spatial frequency ( = f + ). The familiar property that the Fourier transform of a Gaussian is also a Gaussian function and the continuity of the Gaussian roughness spectrum make Gaussian surfaces amenable to analytical solutions in some cases.

10 As shown in (2.6), the Gaussian correlation function is completely specified by two parameters, the standard deviation a (equal to the square root of the variance), and the correlation length, 1. In rough surface scattering studies, o and 1 often appear in products with the electromagnetic wavenumber k, so they are usually specified as the normalized quantities ka and kl. While ka and kl are convenient parameters for the specification of, for example, regions of validity for a rough surface scattering model, care must be exercised when determining ka and kl (or a and 1) from measurements of a surface. This problem will be illustrated in the next section. 2.2 Power-Law Surfaces The studies cited in Chapter I found that linear profiles of many natural surfaces have (one-dimensional) roughness spectra that are well-described (over some measurable range of spatial frequencies) by a power-law function of the form Sz(f) = clfl-, (2.8) where f is the spatial frequency, and c and #/ are constants, with 1 < 3 < 3. P3 is sometimes called the spectral slope; a power-law spectrum is linear with slope -/ when plotted on a log-log scale. The roughness amplitude c is a multiplicative factor scaling the roughness at all spatial frequencies. Power-law surfaces with isotropic statistics have a two-dimensional roughness spectrum that has a similar form: Sz(fr) = afr-, (2.9) where fr is the radial spatial frequency and 2 < y < 4. Power-law surfaces have significant multi-scale roughness, that is, the variations in surface elevation have components with long, intermediate, and short spatial wave

11 10 ''''' 'E 103- '''s,> E 10- 'E 10-4-''-. 10 ) 10. 10- '. 0.1 1 10 spatial frequency, m'1 Figure 2.1: A power-law roughness spectrum in the form of (2.8). lengths, as is evident in the plot of the power-law roughness spectrum shown in Figure 2.1. It is this multi-scale roughness that makes power-law surfaces resemble natural rough surfaces. Just as a hammer or other object is included in photographs of geological formations to show the scale, power-law surfaces also seem independent of scale to the eye. Rough surfaces which have structure over a wide range of spatial scales may also be described using the concepts of fractal geometry introduced by Mandelbrot [34]. Fractal objects are continuous but not differentiable; the fractal dimension Df is a real-valued measure of how a line (surface) fill a plane (space). For a one-dimensional (ID) surface profile, Df takes on values between 1.0 (smooth and differentiable) and 2.0 (plane-filling). Two-dimensional (2D) surfaces have Df between 2.0 and 3.0. Adler [1] shows that the surface profile created by the intersection of a plane and a 2D fractal surface is itself fractal with a fractal dimension equal to that of the 2D surface decreased by one [28]. Random rough fractal surfaces have power-law spectra; Mandelbrot and Van Ness [35] and Voss [60] derive the relation between Df and the spectral slope / of a ID profile of a rough fractal surface: 5- (2.10) Df 2 ' (2.10) 2

12 Thus the spectral exponent satisfies 3 > / > 1 for 1 < Df < 2. Synthetic topographies appear most realistic with Df w 2.2 [34, 60], yielding linear profiles of Df - 1.2 and /3 - 2.6. Estimation of /3 from measured surface profiles is not straightforward; problems associated with the estimation of power-law spectra are described in Chapter IV. The correlation function for a power-law surface is obtained in the usual manner, by inverse transforming the roughness spectrum, but difficulties arise in the powerlaw case. Spectra of natural topography clearly cannot conform to a power-law model at all spatial frequencies-roughness features on the scale of millions of meters are unphysical, as are those at subatomic scales. We therefore must specify the form of the roughness spectrum for both very high and very low spatial frequencies. In general, the piecewise-defined roughness spectrum will not yield an analytical expression for the surface correlation. In addition, there will most likely be some uncertainty as to the accuracy of the assumed form of the roughness spectrum at very high and very low frequencies because these portions of the spectrum have not been measured. Synthetic surface profiles generated from Gaussian and power-law spectra are shown in Figure 2.2. The multi-scale roughness of the power-law surface is evident. A question that often arises when comparing the power-law and Gaussian surfaces is "What are the kao and kl for the power-law surface?" Defining a and 1 for a multi-scale surface is difficult because there are several possible methods that give very different answers. Some examples of the problem are summarized in Table 2.1, in which ka and kl are computed in three different ways for the three artificial power-law surfaces manufactured as part of this study. (The surfaces are described in Chapter V.) The first method used the rigorous definition of the surface corre

13 co L distance (a) 0) distance (b) Figure 2.2: Sample surface profiles generated from (a) a Gaussian roughness spectrum, (2.7), and (b) a power-law roughness spectrum, (2.8).

14 Theoretical Bandlimited Experimental Rz(0,0) O.lf, < f < 10f, (from sample rows) Surface ki ka o kl ka ki k| A (UHMW1) 5768 61.6 9.02 0.986 47.4 3.24 B (UHMW3) 6317 108.8 9.79 0.774 84.1 4.07 C (UHMW2) 6782 194.7 10.43 0.616 108.6 5.56 Table 2.1: Values of kl and ka for the artificial power-law surfaces described in Table 5.1, calculated using three different methods. lation to obtain the variance (equal to the value of the correlation function at zero argument). The correlation was evaluated numerically at zero by integrating (Fourier transforming at zero frequency) the modified power-law spectrum given in (5.8). Using a frequency of 9.75 GHz to obtain k, values of ka and kl for the theoretical method are two to three orders of magnitude higher than those commonly seen in scattering studies. Because some rough surface scattering models specify that ka and kl model inputs should be calculated based on the part of the spectrum that influences the scattering, values of ka and kl were calculated using the same formulation but using a bandlimited roughness spectrum that had been set to zero outside the (arbitrary) limits of 0.1 to 10 times the spatial frequency corresponding to the radar wavelength. These values, listed in the table under "Bandlimited," are much more reasonable when compared to values commonly reported, but it must be noted that the frequency band was arbitrarily chosen. A different set of bandpass limits would produce different values of ka and kl. For an additional comparison, values of kr and kl were "experimentally" deter

15 mined by calculating estimates of the surface correlation function based on sample rows of the synthetic profiles described in Chapter V. This method was used to simulate the results obtained from field data in which neither the roughness spectrum nor the surface correlation function are known beforehand. For each surface, estimates of the correlation function were calculated for three rows, and the resultant values of ka and kl were averaged. As shown in the table, these estimates of kcr and kl are markedly different from estimates produced by either of the other two methods. The experimental estimates varied widely from row to row. Furthermore, the nature of a multi-scale surface suggests that the estimated values of ka and kl will vary with the length of the sample profiles, although this dependence was not verified in this study. The parameters ka and kl are not particularly useful as descriptors of power-law surfaces due to the variety of values they can assume for a single surface. For this reason, power-law surfaces in this study are characterized using the constants c and fi in (2.8) or the corresponding 2D parameters a and 7y in (2.9). The difficulties in estimating these statistical parameters from (often) limited samples of the surface elevation are described further in Chapter IV. Estimation of the roughness spectrum is particularly difficult in the power-law surface case.

CHAPTER III MEASUREMENTS OF SURFACE ROUGHNESS In the last chapter, we saw how rough surfaces may be modeled as realizations of random processes. This chapter and the following chapter describe the methods used in finding the proper random process model for a specific group of volcanic terrains. The first task was the collection of elevation profiles of debris flow surfaces near Mount St. Helens. 3.1 Mount St. Helens Mount St. Helens is the youngest major volcano in the Cascade Range of the northwestern United States. It is located in the southwest corner of the state of Washington, about 70 kilometers northeast of Portland, Oregon. Its rim is 2549 meters above sea level. Mount St. Helens was chosen for this study because it has extremely young terrains and because access to the surrounding debris flows via logging roads is relatively easy. (Mount St. Helens is located within the Gifford Pinchot National Forest.) Mount St. Helens was named by Captain George Vancouver, a British navigator and explorer, in 1792 after the British ambassador in Spain, Baron St. Helens [23]. 16

17 The volcano was dormant for 123 years before seismic activity began on March 20, 1980. On March 27, a system of fractures formed high on the north flank, marking the upper edge of a bulge. By April 12, the bulge measured as much as 100 meters over an area nearly 2 kilometers in diameter. Measurements starting on April 25 showed that the bulge was expanding up to 1.5 meters per day [11]. On May 18, 1980, an earthquake of about magnitude 5 was followed by the breaking away of the bulge and the north flank of the volcano in an enormous avalanche. The debris avalanche consisted of about 2 cubic kilometers of rock debris and glacial ice, fluidized by exploding steam and entrapped air. Traveling at speeds up to 250 kilometers per hour, the debris avalanche extended northeast to Spirit Lake, north to Johnston Ridge, and northwest and west for about 21 kilometers along the North Fork Toutle River valley (Figure 3.1). The debris was deposited in large mounds, with depths up to 150 meters in the center of the valley near Johnston Ridge [11]. The loose debris mixture is quite susceptible to erosion by water, as evidenced by the deep channel cut by the Toutle River since 1980. A geological description of the various debris flows is given by Glicken [27]. The eruption also produced a huge mudflow composed of a mixture of volcanic debris, water from Spirit Lake, the Toutle River, and trapped snow and ice. The mudflow continued down the North Fork Toutle River valley, past the limit of the debris avalanche, and on to the Cowlitz and Columbia Rivers [11]. Debris flow surface topographies were measured in two areas on the debris flows extending along the North Fork Toutle River valley. The principal site (profiled on September 19-22, 1990) was near the west end of Johnston Ridge, about 10 km northwest and within view of the crater (JRS in Figure 3.1). Two distinct surface types were apparent here; elevation profiles were collected on each. The primary

18 British I NJ Columbia \ Wash. aMSH 0 1 2 3 4 5 km / North regon Mount St. Helens ( Figure 3.1: The Johnston Ridge site (JRS) is located about 10 km northwest of the crater. The Elk Rock site (ERS) is located farther west, about 19 km west-northwest of the crater. debris flow surfaces consisted of a loose mixture of rocks, cinders, ash, and dust. These surfaces were extremely rough, making equipment transport difficult on foot. They were formed as part of the original landslide-avalanche that accompanied the main eruption on May 18, 1980. Since then, they have been modified only by erosion. The other surface type observed at Johnston Ridge was much smoother and occurred in low-lying areas between debris mounds. These surfaces were formed either by the mudflows during the eruption or by the deposition of water-borne ash and sediments that were eroded from the primary debris flow surfaces in the years since the eruption. These two terrain types are compared in Figure 3.2. The surfaces at Johnston Ridge were rocky and exhibited less than 15% vegetation cover (less than 5% in many areas). The rougher surfaces supported only sparsely scattered weeds; the low-lying, smoother surfaces had scattered short grasses. Additional surface profiles were collected at a site farther down the North Fork

19 ^ --- —— ^ --- — ^-, ---rr primary secondary primary Figure3.2: The smoother, low-lying secondary surfaces occurred between the rougher (primary) debris mounds. Toutle River valley on September 25 and 26. The Elk Rock site (ERS in Figure 3.1) was located on the valley floor due west of Elk Rock and about 19 kilometers west-northwest of the crater. The terrain here was also deposited by the landslide-avalanche associated with the May 18, 1980 eruption; however, the surface here differs from the Johnston Ridge site in that the smoother mudflow- or waterdeposited surfaces dominate. The floor of the valley, over two kilometers wide, is a plain covered with the smooth deposits accompanied by hummocky mounds of rougher debris. These mounds generally measured less than 25 meters across. Vegetation at the Elk Rock site was less sparse than that nearer the crater, approaching 40% coverage. Grasses on the smooth surfaces were a bit thicker (though certainly not dense), and weeds were more noticeable on the rocky debris mounds. One goal of this study is to investigate scattering effects due to multi-scale surface roughness. The surface roughness of the debris flows was characterized at spatial scales smaller than, on the order of, and larger than the electromagnetic wavelengths of common remote sensing radars (66 cm for a P-band radar, 3 cm for an X-band system) using two measurement techniques. A computer-driven, 2D laser profilometer recorded surface height profiles of square grids with sides ranging from 8 cm to

20 1 m in length. Larger areas measuring 32 meters square were surveyed at one-meter intervals using a surveying level and stadia rod. 3.2 2D Laser Profilometer System A two-dimensional laser profilometer system was designed and assembled for field use at Mount St. Helens. A laser-based system was chosen because it allows measurements of surface elevation without requiring physical contact that may disturb the surface. Easy adaptability to computer control was an additional advantage over contact-probe profiling methods. Requirements for the profiling system were as follows: 1. Transportable by two or three people up to one mile away from vehicle 2. Operable from battery power 3. Unattended operation, freeing operators for other tasks 4. 1 square meter scan area 5. Reasonable cost The 2D Laser Profilometer System is pictured in Figure 3.3. A list of components and accessories follows: * Pulsar 50 Surveying Electronic Distancemeter (EDM) (manufactured by GEO Fennel) * XY table and stepper motor controllers (manufactured by Jasta, Inc.) * Zenith PC-compatible laptop computer * Two 12 V, 105 ampere-hour deep-cycle marine batteries

21 Figure 3.3: The 2D Laser Profilometer System.

22 * Four heavy-duty photographic tripods * Protective frame for mounting and transport of XY table * Case for external switches, stepper motor controllers, and power distribution * Army surplus stretcher for carrying computer, switch box, batteries, tripods, and EDM * Counterweight for EDM mount The Zenith computer and EDM were owned by the Radiation Laboratory. The EDM had been used previously in a linear (ID) profilometer system built by Roger De Roo. The EDM is the principal component of the profilometer system. It uses an infrared laser to measure the distance from a reference plane to a target surface. The EDM and counterweight are mounted on the XY table with the laser directed downward through the table frame. The XY table consists of three ball-screw assemblies that are driven by two stepper motors. The profilometer software, written in Microsoft QuickBasic, controls the EDM and XY table through the serial port of the Zenith laptop computer and records surface elevation data from the EDM in disk files. The computer has only one serial port (and no means of adding expansion cards), so it was necessary to build an external switching circuit, controlled through the parallel port, to allow communication with either the EDM or the XY table as needed. The switching circuit is shown in Figure 3.4. Automatic operation of the EDM was complicated by the random occurrence of "HARDWARE ERROR" messages. I was unable to determine the cause of the errors myself or through inquiries to the manufacturer. The error could be cleared by

23 Figure 3.4: Circuit diagram of external serial port switch.

24 cycling the power, but there was no means of doing so remotely through the serial interface. A workaround was devised whereby the EDM battery leads were routed through an external relay controlled through the computer's parallel port. A lack of response within seven seconds of a measurement command prompts an EDM power interrupt by the computer. The scan location of the hardware error is noted in the data file of measured surface heights, and the measurement is attempted again. On a grid of 2500 points, it was not uncommon to have 15 to 30 hardware errors. No correlation was observed between the error frequency and duration of operation. The hardware errors have not been found to affect measurement accuracy. Two 12 V marine batteries provided DC power for the computer, stepper motors, and EDM, while a third battery was held as a spare. A full day of operation did not seem to weaken the batteries significantly; nevertheless, the batteries were charged nightly and rotated with the spare. The XY table was mounted in a protective frame and supported by four heavyduty tripods. The table was carefully leveled at a height of approximately 1.5 meters-near the minimum working range of the EDM. The EDM laser has a spot diameter of about 1.5 mm. It was used in its most precise mode, in which the standard deviation of the measured surface height is 3 mm. Measurements using this mode require two to three seconds per surface point, depending on the sampling interval, A. A typical scan measuring 10 cm x 10 cm with A = 2 mm (2601 points) takes 1.8 hours. To reduce on-site times to acceptable levels, each surface was scanned along grids of 40 x 40 or 50 x 50 points using at least two values of A (e.g., 2 mm and 1 cm), rather than scanning a single large grid using a very small A. One of the rougher scanned surfaces is shown in Figure 3.5. The surfaces had poor reflectivity at the laser wavelength, causing EDM measurements to fail. Painting

25 Figure 3.5: One of the rougher surfaces scanned by the EDM (scan site 3, Johnston Ridge). The surface has been spray-painted to increase its reflectivity. the surfaces (including the surface in Figure 3.5) with white or yellow spray paint increased the reflectivity. The spray can was held sufficiently far from the surface to avoid disturbing the small rock fragments. 3.3 Profilometer Data Processing Errors in the profilometer data prevented the direct application of standard spectral estimation techniques. Replacing the data by repeating the on-site work was not feasible; instead, error effects were minimized through a combination of error detection and correction algorithms and a modified spectral estimation technique. The errors were of two types: incorrect surface elevation values due to overheating of the EDM, and intermittent level shifts due to an instability in the EDM.

26 Figure 3.6: Grayscale plot of a surface elevation grid (scan 2b, Johnston Ridge site) showing a single outlier due to overheating. High points are light; low areas are darker. The black pixel is the outlier. 3.3.1 Overheating Errors Two factors combined to produce an overheating of the EDM that had not been observed indoors: unusually hot weather (sunny and clear with highs of 27-31 'C) and placement of the EDM such that the underside and battery connection module (both colored black) were exposed to direct sunlight. A sun shield installed for the last scan at Johnston Ridge reduced the frequency of these errors. Overcast skies during the Elk Rock debris flow scans prevented overheating errors at that site. The overheating introduced negative outliers in the measured surface elevation grid, i.e., measured elevations that were markedly lower than those at adjacent points. The outliers resembled holes in the surface and were easily distinguishable in plots of the surface elevation (see Figure 3.6). A variety of filtering schemes for the removal of these bad pixels were examined. The simplest was the imposition of a lower limit beyond which all pixels were classified as outliers. This technique failed for surfaces not having low relief.

27 25 20 o 15 -CD E 10 --0.10 -0.08 -0.06 -0.04 -0.02 0.00 0.02 residual (m) Figure 3.7: A histogram of a bimodal distribution of median residuals. The left group, corresponding to outliers, is clearly separable from the main distribution. The majority of the outliers were removed using median residual filtering, in which a statistic called the median residual, rm, is calculated for each pixel by subtracting the median, m50, of the N x N surrounding subregion from the surface elevation, Z(i,j), at the pixel of interest: rm(i,j) = Z(i,j) - m50. (3.1) Pixels with overheating errors have large (negative) rm because their elevations differ greatly from the surrounding pixels. Residuals outside some lower bound are deemed unphysical, and pixels associated with these residuals are designated as outliers. A' might range from 5 to 9-I used the smallest N for which the median was stable. If enough outliers were present, the median of the entire scan was used for the first iteration. A data set ideally suited for median residual filtering would have a bimodal distribution of residuals in which the modes were widely separated. For example, there would be little uncertainty in labeling the left group of residuals in Figure 3.7 as being due to outliers. Of course, such a clear separation of residuals is rarely the case. More often one faces a distribution like that shown in Figure 3.8, which

28 30 -25 -20 -15 0' 10 --0.6 -0.4 -0.2 0.0 residual (m) Figure 3.8: A histogram of median residuals (scan 3a, Johnston Ridge) for which values due to outliers are not easily distinguishable. is from a measured data set (scan 3a, Johnston Ridge site). Rather than choosing some arbitrary cutoff (and risking the elimination of good pixels or the acceptance of outliers), one can use an iterative approach, setting the cutoff at a conservative, low value such that confidence is high that pixels with residuals below the cutoff are indeed outliers. A new set of residuals based on the remaining pixels is then calculated, and the process is repeated. While only a small portion of the erroneous pixels are removed in each iteration, the majority of the outliers may be eliminated in this manner. After a number of iterations, the median residual filter becomes impractical because residuals due to overheating errors become indistinguishable from those due to edges of rocks or pixels having few neighbors because they are located near the edge of a scan. One can decrease the area over which the median is calculated, but cutoff selection eventually becomes subjective, and it becomes difficult to determine when to terminate the filtering process. The problem at this stage is one of distinguishing between the two distributions of N x N subregion surface elevations (not residuals) shown in Figure 3.9. If mi50

29 10 8 I 6 1r m50 4 2 x 0 I -0.4 -0.2 0.0 0.2 elevation (m) (a) 10 m0O 2 25 -0.4 -0.2 0.0 0.2 elevation (m) (b) Figure 3.9: Elevation distributions for 5 x 5 subregions on two surfaces. Point X in (a) is probably not an outlier, but point X in (b) probably is, yet both have the same median residual.

30 designates the median height for the N x N region, then we see that the median residual calculated at point X is the same for both surface A and surface B (-0.3 m). However, we can see from the from tsurfae surface height distribution that pixel X on surface A is less likely to be an outlier than pixel X on surface B. One way to quantify this difference is to compile a lower quartile difference statistic for each pixel. The lower quartile point, m25, is the pixel for which 25% of the other pixels have lower elevations and 75% have higher elevations. The lower quartile difference, dq, is calculated for each 5 x 5 pixel subregion similarly to the median residual: dq(i j) = Z(ij )-m25 (3.2) A high value of Idq I indicates that the pixel elevation differs significantly from most of the other points in the subregion and should therefore be discarded. The cutoff value of dq in this study was determined in a reproducible manner by plotting a histogram of all dq values in a data set and discarding those pixels having a dq value outside the main lobe of the distribution. In summary, the median and quartile difference filtering procedure consists of the following steps: 1. Compute median residuals, rm, using the smallest subgrid (9 x 9, 7 x 7, 5 x 5) that does not corrupt the median. (For some data sets, it will be necessary to use the median of the whole grid in the first iteration.) 2. Filter out pixels whose residuals exceed some conservative limit-whose residuals are sufficiently large that there is a high probability that the pixels are outliers. Increase the lower residual limit only moderately in successive iterations.

31 3. Repeat steps 1 and 2 until the obvious outliers are fairly sparse; i.e., there remain no more than two in any 5 x 5 subregion. Steps 1 and 2 may be unnecessary in data sets where the outliers are already sparse (e.g., scans 2b and 2c). 4. Compute the lower quartile difference dq for each remaining pixel based on valid pixels (those not filtered out in previous steps) in a 5 x 5 subregion. (The scan site 4 data sets had so many points missing that 7 x 7 subregions were necessary.) 5. Determine the lower cutoff value for the quartile differences using a histogram of dq values. Set the limit equal to the first minimum on the left of the peak or the point on the left at which the histogram has a definite corner. 6. Filter out pixels whose dq values are below the limit set in step 5. The median residual and lower quartile difference filtering procedures are illustrated in Figure 3.10. Figure 3.10(a) shows an unfiltered profilometer scan measuring 100 cm x 100 cm in which the pixel spacing is 2 cm. In Figure 3.10(b), each pixel has been replaced by the median of a 9 x 9 subregion; we see that the 9 x 9 subregion medians are corrupted badly in the last few rows. (This is not unexpected because surface scans are completed from top to bottom; overheating errors become more frequent as the EDM becomes warmer.) Five rows and five columns were discarded because these rows were considered uncorrectable. Figure 3.10(c) shows the result. Figures 3.10(d), 3.10(f), and 3.10(h), show histograms of median residuals based on 9 x 9 pixel subregions. Lower cutoff values of -30, -20, and -14 cm led to the filtered data sets shown in Figures 3.10(e), 3.10(g), and 3.10(i), respectively. (Filtered pixels are displayed as white.) The isolated dark pixels in Figure 3.10(i) are the re

32 II U S*F Ux (a) (b) p 1 U. I; i'.r UI*U' ~ I U (c) 30 -25 -20 -1 5 -10 -5-.....,_.I..li -.. I i. I- I —............... -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 residual (m) (d) 30 -25 -> 20 -1 5 -10 -5 - rrfrrtf~rvflrlbnn taMMMN Oi n IIIlllli,,,,,,, -0.3 -0.2 -0.1 residual (m) 0.0 0.1 0.0 0.1 (e) (f) Figure 3.10: Successive stages of the median residual and lower quartile difference filtering procedure for scan 3a, Johnston Ridge site. The steps leading to the final filtered set (k) are explained in the text.

33 (g) 30 2 5 2 20 1 15 1 0 5 0 --0.20 -0. 15 -0.10 -0.05 0.00 0.05 0.10 residual (in) (h) 30 -2 5 2 0 1 15 10 5 -0 -0.10 -0.05 0.00 0.05 residual (in) (j ) (i (k) Fig ure 3.10: (continued).

34 maining outliers, for which the surface height appears artificially low. Lower quartile differences were computed for each remaining pixel in Figure 3.10(i). A histogram of the 1871 dq values appears in Figure 3.10(j). (The peak of the histogram is 561, well off the top of the figure.) Using the histogram, a lower cutoff was set at the first minimum and all pixels with dq < -15 cm were discarded. The resulting surface elevation data set is shown in Figure 3.10(k). Note that the isolated dark pixels have been removed while dark pixels in low areas are intact. The median residual and lower quartile difference filtering procedures seem to eliminate outliers accurately with minimal loss of valid pixels. Because the modified spectral estimation technique described in Chapter IV is most direct when using equally spaced data points, it was necessary to replace the pixels removed by the filtering process. The next task, then, was to find a method of assigning a value to each missing pixel such that the spectral estimate based on the "patched" profile is as accurate as possible. Median interpolation was selected as the interpolation method. Tests in which holes were introduced and then replaced (via median interpolation) in synthetic data sets showed that median interpolation has a minimal effect on the resultant roughness statistics when spectra from different rows are averaged. (These tests are described more fully in Chapter IV.) The median interpolation algorithm consists of two steps: 1. Find the median of the 3 x 3-pixel subregion around each hole and count the number of valid pixels. 2. Replace the hole by the median of the 3 x 3 pixel subregion if at least 5, 3, or 2 non-holes are present in the inner, edge, and corner-pixel cases, respectively.

35 For data sets with numerous adjacent holes, multiple passes of the interpolation algorithm are necessary. Steps 1 and 2 are repeated until all holes are filled. The procedure is illustrated in Figure 3.11. The filtering and interpolation steps performed on each data set are listed in Tables A.1 through A.7 in Appendix A. 3.3.2 Level Shift Errors The second source of error in the profilometer scans was an occasional shift in measured elevations due to an instability in the EDM. After 20 to 30 minutes of normal operation, the EDM reference plane would shift by 2-7 mm and then become stable again. The shifts occurred at irregular intervals and did not seem to be correlated with overheating errors or with the occasional hardware errors. The level shifts are manifested as horizontal bands across the profilorneter scans. Because the magnitudes of the level shifts were rather small (only slightly larger than the standard deviation of the EDM), they were noticed only in gra;yscale plots of very smooth surfaces at Mount St. Helens. The presence of the level shifts was verified by running a 51 x 51-point profilometer scan with the XY table disabled to prevent movement. The sequence of 2601 points was mapped onto a square grid as if the EDM had been moving. (The total scan time was one hour and 46 minutes.) The grayscale plot of this grid is shown in Figure 3.12. The pixel at the upper left corner represents the first measurement. Subsequent pixels follow along this row to the right, then continue along row 2 from right to left, and then continue along successive rows in alternating directions. The only variation in Figure 3.12 should be due to the limited precision of the instrument, but there are obvious level shifts between (for example) rows 15 and 16 and between

36 (a) (b) (c) (d) Figure 3.11: Median interpolation for scan 3a. (a) The data set immediately after medial residual and lower quartile difference filtering. (b) The data set after one iteration of median interpolation. (c) After three iterations. (d) After six iterations.

37 Figure 3.12: A profilometer scan showing the level shifts. (The EDM did not move during this scan.) rows 21 and 22; these shifts measure 4.42 mm and 7.62 mm, respectively. The level shifts were small compared to changes in elevation due to topography; nevertheless, I felt that they might corrupt the spectral estimates enough to make their mitigation or elimination worthwhile. To avoid the effects of the level shifts, I assumed isotropic statistics and employed a procedure of spectral estimation by linear sampling; i.e., I used linear profiles (rows of the profilometer scan) as samples of the surface. By using individual scan rows and removing the mean surface height from each, the level shifts are avoided because the reference level is stable within most rows. The shifts usually occurred over a sequence of 10 to 15 points, so rows with visible level shifts were discarded. (The discarded rows for each data set are listed in Table 3.1.) Averaging the spectral estimates obtained from different rows within a scan reduces the effect of any remaining (undetected) errors as well as the variance inherent to the surface random process. The procedure for spectral estimation by linear sampling is described in detail in the next chapter. The assumption of isotropic surface roughness statistics seems reasonable due to

38 Scan Data file Rows with visible level shifts 2b lp2bdata.fql 1, 2, 5, 8, 13, 14, 24, 25 2c lp2cdata.fql 4, 5, 6, 7, 8, 13, 22 3a lp3axdata.fq4 13, 14, 24, 28 3b lp3bdata.fq4 3, 4, 24, 25, 35, 36 3c lp3cxdata.fq6 None discernible 4a lp4axdata.fq5 1, 2, 4, 5, 10, 13, 18, 29, 30, 31, 40 4c lp4cdata.fq5 1, 2, 8, 33, 34 5a lp5adata.all 1, 2, 16, 21, 38, 41, 59, 61, 71 5b lp5bdata.all 1, 19, 32, 33, 42, 43 6a lp6adata.all 1, 2, 3, 4, 23, 24, 25, 26, 28, 29, 30, 38 6b lp6bdata.all 1, 2, 6,11, 24, 25, 38 Table 3.1: Rows excluded from spectral estimates due to visible level shifts

39 the structure and origin of the surfaces. The flow units examined were formed in a debris avalanche when the north flank of the volcano collapsed. The avalanche was composed of a loose mixture of rock and ash that had no obvious directional structure, unlike a viscous lava flow which has a definite flow direction. While the assumption of isotropic statistics can be neither proven nor disproven from the data, collected as part of this study, the isotropic statistical model provides a useful firstorder approximation for initial scattering studies such as the present work. 3.4 Surveying Measurements Surveying equipment was used to collect elevation data over larger areas. The surveying tools included a self-leveling surveying level, a stadia rod (both on loan from the Department of Civil and Environmental Engineering), three 100-meter lengths of airplane cable (multi-strand stainless steel cable about 1/16 inch in diameter) with markers (crimp-on wire splice connectors) installed at one-meter intervals, and eight pitons that were used to anchor the cables in the ground. The surveying level is essentially a telescopic sight mounted on a tripod. The sight may be rotated about the vertical axis; its plane of rotation becomes a reference plane. Surveying measurements consisted of the following steps: 1. Using the angular scale on the surveying level to obtain right angles, lay cables on three sides of the area to be surveyed, as shown in Figure 3.13. The area surveyed was 32 x 32 meters for each of the four surveys, although the first survey was initially configured for a 50 x 50-meter area. 2. Anchor the two cables on opposite sides using pitons (cables A and B in Figure 3.13).

40 A B surveying level ' Figure 3.13: Top view of survey configuration. 3. Anchor the third cable (C) on top of A and B. This cable will be moved along cables A and B as each row is completed. 4. Mount the surveying level securely in place. 5. The stadia rod bearer places the stadia rod on the ground at the marker and holds the rod vertically. 6. The surveyor reads the elevation on the stadia rod (the distance between the reference plane and the base of the rod) by noting where the cross-hair in the level's telescopic sight crosses the rod. He announces the elevation to the recorder, who records it and signals the rod bearer. 7. The rod bearer then moves along C to the next position (1-4 meters away), and the last two steps are repeated. This continues until the end of the row. Every fifth marker has been painted-the rod bearer signals the recorder when moving to these markers as a position check.

41 8. After the row is complete, cable C is moved 1-4 meters along cables A and B to mark the next row to be surveyed. The surveying continues along the new row. The surface height was measured at intervals of 1, 2, or 4 meters, depending on the surface roughness. Surveying was very slow. In spite of the time advantage gained through use of the self-leveling feature of the surveying level, each data point required 20 to 40 seconds to measure and record. In surveys 1 and 4, the level was used in multiple locations either because the surface height variation exceeded the length of the stadia rod (survey 1) or because the entire grid was not visible from a single location (survey 4). The principal source of error in the surveying measurements was difficulty in reading the stadia rod when the wind caused it to sway. When wind was a problem, the surveyor took readings more slowly. Slight departures of the rod from the vertical should not introduce large errors. 3.5 Measured Profiles Grayscale plots of all profilometer and survey data are shown in Figures A.1 through A.6 in Appendix A.

CHAPTER IV SPECTRAL ESTIMATION The surface elevation profiles described in the previous chapter provide quantitative measures of the surface roughness at specific locations on selected debris flows near Mount St. Helens. This data has practical value only through its ability to represent statistical properties of large debris flow areas. In Chapter II, rough surfaces were modeled as random processes of surface elevation; we therefore seek a description of the statistical properties of the random process. One key descriptor is the roughness spectrum or power spectral density of surface elevation. The roughness spectrum is useful because it allows surface structure to be interpreted as a sum of sinusoidal components with differing wavelengths. Some rough surface scattering models utilize the roughness spectrum directly, while others use its inverse Fourier transform, the surface autocorrelation function. This chapter describes methods of finding estimates of the roughness spectrum based upon finite samples of single realizations of rough surface random processes. We shall see that this task, called spectral estimation, is not trivial. Kay [32] states the spectral estimation problem (for a iD data set) formally: "Based on the N contiguous observations {x[0],x[1],...,x[N - 1]} of a single realization of a WSS [wide sense stationary] random process, it is desired to estimate the PSD [power 42

43 spectral density] for -1/2 < f& < 1/2." (I have departed from Kay's notation and written [normalized] frequency as fi for reasons that will be explained later.) I should probably italicize the word estimate in Kay's definition to emphasize that we cannot obtain an exact solution for the PSD in general because an infinite amount of information would be required. A good estimate is the best one can hope to achieve; the degree of "goodness" will be dependent on the quantity of sample data, the form of the true PSD, and the suitability of the estimation algorithm. Signal processing and spectral estimation texts (e.g., [32], [36], and [47]) list a wide variety of spectral estimation algorithms ranging from classical Four-ier methods to autoregressive spectral estimators to minimum variance estimators. The novice who asks "Which is the best spectral estimator?" will quickly learn that there is no simple answer. Various spectral estimators have different strengths and weaknesses. Selection of an estimator is very application-dependent, involving trade-offs among desirable characteristics such as low bias, low variance, high resolution, low leakage, and minimal spurious features. Many spectral estimators are based upon assumptions regarding the form of the underlying spectrum. Such an estimator may provide enhanced accuracy if the assumed model is correct (or may fail miserably if it is not). It is therefore advantageous to have some knowledge of the form of the spectrum when selecting an estimation algorithm. Other factors influencing the selection of an estimator are the quantity and character of sample data (e.g., irregularly sampled data sets require special measures) and the specific questions one seeks to answer. For example, we may want to determine whether a signal is present among background noise and may not need to know its strength precisely. Conversely, as in the present work, we might assume that no narrow signals are likely to be present and may only want to estimate the magnitude

44 of the spectrum across a band of frequencies. We would therefore be less concerned with high resolution and would seek an estimator with low bias. Another issue, perhaps less appreciated than it should be, is the importance of the sample size N. Most spectral estimators assume large values of N-very large values, from a remote sensing viewpoint. This dependence on N may be overlooked when consulting a text on spectral estimation theory. The importance of sample size is perhaps best stated in Numerical Recipes [48, p. 456]: in the description of variance given in the chapter on the statistical description of data, the authors explain the presence of a factor of 1/(N - 1) rather than 1/N and remark: We might also comment that if the difference between N and N - 1 ever matters to you, then you are probably up to no good anyway-e.g., trying to substantiate a questionable hypothesis with marginal data. For purposes of spectral estimation, then, we would like to have surface profiles with very large N. This is not always possible: in remote-sensing studies, surface profiles are often only one of many "ground-truth" parameters to be collected. As a result, investigators are often faced with the problem of estimating the surface spectrum using a limited data set. In the case of ID profiles, a typical record might consist of at most 500 points; 2D profiles are even shorter, rarely exceeding 50 x 50 points. Numbers of this magnitude are considered small from a statistical inference viewpoint, where samples of thousands or hundreds of thousands of points are more common. The profiles in this study are limited; we shall consider this when selecting an estimator. In this chapter, I review some spectral estimators, examine some difficulties in spectral estimation of power-law spectra, determine an estimation algorithm that avoids these difficulties, and apply the algorithm to the surface profiles described in

45 the previous chapter. (Much of the material here has been summarized in recent papers by the author [2, 3]). For simplicity, one-dimensional spectral estimation will be considered first. 4.1 Classical Spectral Estimators Fourier-based spectral estimators (e.g., the periodogram) render an estimate of the power spectral density whose expectation or mean value is a convolution of the spectrum of the sampled profile and the Fourier transform of a window function introduced by the finite extent of the sample profile. Several window functions are used; common windows include rectangular, triangular, and Hanning windows. The selection of a particular window is based on a trade-off between spectral resolution and spectral leakage. Spectral leakage refers to the inaccuracy of a spectral estimate at a given frequency due to convolution with the transformed window function; i.e., spectral power "leaks" from nearby frequencies according to the shape of the transformed window. Spectral domain window functions corresponding to a periodogram and a modified periodogram with Hanning window are shown in Figure 4.1. At each point or frequency at which the spectrum is to be estimated, the true spectrum is weighted by the window function centered on that frequency and then averaged over frequency. Contributions from frequencies other than the frequency of interest constitute spectral leakage. Spectral leakage can be reduced by choosing an estimator whose window has side lobes which decay quickly (such as WH in Figure 4.1), but such a window will have a broader main lobe, decreasing the spectral resolution. (See [48] for a discussion.) Spectral estimates derived from sampled data suffer from aliasing if the sampled

46 NE 1U - -.- - WT - 1 — WH ' - - power-law, — 2.8 C 0 10- \ i i i;., 10.& I[-.,F..._! X,, l, 'la' 'gtl, 0.001 0.01 0.1 10 I" I I I I m I, result in a spectral bias an inaccuracy in the mean value of the estimate). 0.001 0.01 0(4.1 where A is the sampling interval. In the sections that follow, it is assumed that Waliasing is not present. Modificated perions due to aliasth Hanning will be presented in a later valuedwith N = 256. The convolufunction of(i.e., there are no overhangs). Suppose withat the surface height has.spectrum (e.g., the power-lawiodogram spectrumal estimator, PpE(f)hown asfor the samplshed linsurface heightmay reis given bya spectral bias (an inaccuracKay [32]:in the mean value of the estimate) The hat in PPER denotes an estimate. The frequency is written NA as a reminder that (4.2) is givthe sampling interms of a normal. In the sections that follow, it is relatssumed to theat aliasing is not present. Modifica tions due t o aliasing will be presented in a later is given by Kay [32]:

47 spatial frequency in m-1 by fA = fA, (4.3) and fa satisfies -1/2 < fA < 1/2. Most spectral estimation texts use normalized frequency; I shall follow this convention until a later section when considering how to combine spectra obtained from sampled profiles with differing A. While the periodogram may be evaluated for any IfAl less than 1/2 (the normalized Nyquist frequency), it is often calculated using a fast Fourier transform (FFT) and consequently is evaluated only at discrete frequencies fAi = i/N, where i = -N/2,..., -1,0, 1,..., N/2. If additional frequencies are desired, the original data series should be zero-padded to length M by adding M - N zeros to the end. The periodogram will then be evaluated at frequencies fAi = i/M. Most FFT algo-. rithms require that N (or M in the case of an extended series) be equal to an integer power of two (e.g., 256, 1024). Series for which N is not equal to an integer power of two must be zero-padded up to the next higher power of two. The periodogram is an unreliable estimator of the power spectral density because it has a variance which is equal to the square of its expected value and independent of N [32]. The usual practice in the field of spectral estimation is to average multiple periodograms to reduce the variance, but this may not be possible if the quantity of data is limited.1 A modified periodogram PHAN(fa) may be obtained by multiplying the data series by a Hanning window before calculating the spectral estimate: 1 N-1 2 PHAN(fA) = NU Z Z[x]wH[n] exp(-j27rfan), (4.4) n=O lFor example, suppose the data record is 1000 points long. To reduce the variance by a factor of 8, we can divide the data record into 8 segments of length 125, calculate 8 spectral estimates, and average the results. The bias of the estimate will be increased (due to the broader main and side lobes of the window), and the reduction in variance may be less than eightfold if the segments are not uncorrelated. If the entire data record has only 125 points, this procedure is not useful.

48 where the Hanning window, H[n]= [1-cos( 2 ] (4.5) replaces the rectangular window (wR[n] = 1) implicit in (4.2), and an extra normalization factor is introduced: 1 N-1 U -N E w[n]. (4.6) n --- —O U=NZw [n]. n=0 (Equation (4.4) follows the formulation of Welch [61] using a single data segment.) The transform of the Hanning window wH[n] has side lobes which are lower and a main lobe which is roughly twice as wide as that of the transform of a rectangular window WR[n] of the same width. Convolution with a transformed Hanning window results in less leakage from low spatial frequencies (because the side lobes are smaller) at the cost of decreased spectral resolution (smoothing due to the wider main lobe). One must consider the smoothness of the spectrum and its roll-off rate and then decide whether the leakage reduction is worth the loss in resolution. 4.2 Estimation of Power-Law Spectra As stated in the beginning of this chapter, the selection of a spectral estimator should be guided by the expected form of the unknown spectrum. The previous studies mentioned in Chapter II suggested that the debris flow roughness spectra might conform to a power-law spectral model. The consequent task was to determine the best estimator for a power-law spectrum using a limited data set. 4.2.1 Fourier-based estimators We begin by examining the expected values or ensemble averages of the classical, Fourier-based estimators described in the previous section, when applied to a onedimensional surface profile having a power-law roughness spectrum. The expectation

49 of the periodogram (4.2) is given by [32]: E [PPER(f1)] = l/2 WT (fA - 0)SZA (0) d~, (4.7) where WT(fA) is the transform of a triangular window: WT(fA) = (sin 7rf) ' (4.8) and SzA (fA) is the normalized form of the true power spectral density of Z(x) chosen such that f SzA(fA) dfA = fSz(f) df: Sza(f A) = 1Sz(f) = cf= A-lc(fA)-) Sza(fa) = cf2, (4.9) where we have defined ca = AO-lc. The expected value of the modified periodogram estimator PHAN (4.4) has a form similar to (4.7) [47, p. 553]: Af 1/2 E [PHAN(fA)] = / WH (fA -() Sza () d, (4.10) where WH(fA) = 4U (+ bl + b + b3) (4.11) sin(rfAN) (4.12) = sin{lr[fA - N1}iN} sinlir[fa - N-i]J b2 = sin{[/a- N]} ' (4.13) sin{lr[fA + N1 1]N} b3 = sin{r[fA+ ]}' (4.14) and U is given by (4.6).

50. 102- -- 1=2.0 - = 2.4 E 1 * 2.0, nozero-paddin ~ 100 10 1-3 10'2 M 0.001 0.01 0.1 normalized spatial frequency, dimensionless Figure 4.2: Expected value of the periodogram spectral estimator for power-law surfaces with / = 2.0 and 2.4 and ca = 1.0 x 10-6. The estimates are based on 256-point profiles which were zero-padded to 16 384 points so that the periodogram would be evaluated at many frequencies. Solid dots indicate the expected values of the periodogram when evaluated without zero-padding. The performance of these Fourier-based estimators may be determined by calculating their expected values using (4.7) and (4.10). The singularity at the low end of a pure power-law spectrum is eliminated by modeling the exact surface spectrum as a Rayleigh function at low spatial frequencies: (a/la2)lflexp(- f/[2a2]), IfAl 0.0001 SzA(fA) = CAlfAl-, 0.0001 < IfAl < 1/2, (4.15) 0, If > 1/2 where the parameters a and a are chosen by enforcing the continuity of Sza and its first derivative at IfA| = 0.0001. Aliasing has been avoided by setting the spectrum to zero for IfAl > 1/2. Figure 4.2 shows the expected value of the periodogram spectral estimator calculated from theoretical spectra with spectral exponents of 2.0 and 2.4. The oscillatory behavior is due to the side lobes of WT(fA). These oscillations are not visible if the periodogram is evaluated only at frequencies ft = i/N (i.e., without zero-padding) as indicated in the figure. Since it is difficult to fit power-law functions to such

51 E....... 1=1.6 = 2.8. 10~ ''.. =.=2.0 0.01 0.1 Figure 4.3: Expected value of PPER for five 256-point power-law profiles evaluated at frequencies corresponding to the case of no zero-padding. Nc 10~ -- 1.2 -- P=2.4 - - - --- 1.6 ---- P2.8 10.1. '. ' 10 -1 0. 0.01 0.1 normalized spatial frequency, dimensionless Figure 4.4: Exact values of the power-law spectra used in calculations of the expected values of the spectral estimates. oscillatory estimates, one may sample the expectations of periodograms of surfaces with five different spectral exponents at frequencies corresponding to the case of no zero-padding. zero-padding (Figure 4.3).2 These curves may be compared to the exact spectra in Figure 4.4.3 While the exact spectra in Figure 4.4 vary in slope, the periodogram — 1. -13=-2.8 2It is emphasized that Figure 4.3 shows expected values, E[PpER(f)], (mean values) of the -2 periodogram estimator, calculated using equation (4.7) at frequencies corresponding to the case of no zero-padding for random profiles Z(z) having power-law spectra given by (4.15) with five different values of,8. Plots of the estimator itself, PPER, would show wide variance and departures from linearity (similar to the estimates shown in Figure 4.8) and would vary for different realizations r-' 0.01 0.1 of the random processa i.e., for differequencyt samples of a given surface. The expected values of the estimator appear linear at the low-frequency end because the assumed valform of the true surfac the spectrum (4.1) is power-law (with linear slope) down to a frequency of oscillato = 0.0001a frequency too low to be accurately measured by the sampled profile segment. This he expected values of the periodogram estimator in Figure.3 have a higher total energya in level than the corresponding exact spectra in Figure 4.4 due to the spectral leakage that is also of no zero-padding for random profiles Z(x) having power-law spectra given by (4.15) with five different values of fi. Plots of the estimator itself, PPER, would show wide variance and departures of the random process, i.e., for different samples of a given surface. assumption was used because surface profiles are often too short to show very long wavelength roughness components. (This was certainly true for the field work at Mount St. Helens.) 3The expected values of the periodogram estimator in Figure 4.3 have a higher total energy level than the corresponding exact spectra in Figure 4.4 due to the spectral leakage that is also

52 104. E 102- -1.2 -- 132.4 * -_ \ — 1.6 - - 3=2.8 10~ ---,-2.0 104 1. 10 0.01 0.1 normalized spatial frequency, dimensionless Figure 4.5: Expected values of PHAN(fA) for five different 256-point power-law profiles. The profiles were zero-padded to 16 384 points. estimates are insensitive to the slope for several values of /3. Fitting power-law functions to the linear portions of E[PPER], fA < 0.2, we see (Table 4.2.1) that the estimated spectral exponent /3 clusters near 2.0 for 2 < /3 < 3.4 These values of /3 correspond to fractal dimensions that are typical for natural terrains (1.0-1.5 for profiles, 2.0-2.5 for surfaces). This slope insensitivity due to spectral leakage may well explain the number of studies [6, 30, 31, 52] reporting /3 = 2.0 for a variety of natural topographies. Slope insensitivity is not a problem with PHAN(fA), whose expected value is shown for various /3 in Figure 4.5. The expected values of these estimates closely approximate the exact values for frequencies for which the main lobe and first two side lobes do not overlap the low-frequency peak (fA > 0.015625). (Power-law fits responsible for the slope insensitivity. Remember that E[PPER] is the convolution of the true spectrum Sz, with the transform of a window function introduced by the finite surface sample (WT in (4.7)). If the estimator were perfect, WT would be a delta function, and the estimate at frequency f would reflect only power in the true spectrum at that frequency. In the periodogram case, however, WT is given by (4.8) (a sinc2 function), so power from many frequencies above and below f affect PPER(f). Leakage from lower frequencies is particularly significant in the power-law case. 4For a given single periodogram, /3 may indeed take values between 2.0 and 3.0, due to the variance of the periodogram. Long data records do not reduce this variance; the variance is independent of N (see [20], page 422). One may reduce the variance by segmenting the long data record as described previously, but this process will lead to the expected values shown in Figure 4.3. Sampling PPER at frequencies other than fAi = i/N could lead to a variety of incorrect slopes, as can be seen in Figure 4.2.

Theoretical E[PPER] E[PHAN] E[Ppw] E[PCAp] ca = 1.0 x 10-6 Ca- = 9. 10-7 99 x 10-7 A = 8.42x 10-7 ia = 1.01 X 10-6 / =1.2 / =1.231 = 1.201 / = 1.255 p= 1.198 ca = 1.0 x 10-6 =1.04x 10-6 CA = 9.99 x 10-7 C = 8.63 x 10-7 = 1.01 x 10-6 / = 1.6 = 1.656 / = 1.601 = 1.643 = 1.598 c = 1.0 x 10-6 A = 2.10 x 10-6 c = 9.96 x 10-7 C = 8.66 x 10-7 Ca = 1.02x 10-6 = 2.0 = 1.982 = 2.001 = 2.042 /= 1.998 = 1.0 10-6 a = 1.87 x 10-5 Ca = 9.98 x 10-7 = 8.65 x 10-7 C = 1.02 x 10-6 p=2.4 / = 2.032 = 2.402 / = 2.443 / = 2.397 ca= 1.x 10-6 CA =3.83 x 10-4 a = 9.96 x 10-7 = 8.43 x 10-7 C = 1.02 x 10-6 / =2.8 = 1.986 = 2.804 3 = 2.859 1 = 2.796 Table 4.1: Power-law parameters derived from fits to the expectations of periodogram (1/N < fA < 0.2), modified periodogram (4/N < fA < 0.5), pre-whitened (1/N < fa < 0.2), and Capon's (1/[2p] < fA < 0.5) estimators. cn 00

54 to the upper regions of these estimates are listed in Table 4.2.1.) While PHAN has a mean value which closely approximates the exact spectrum, this estimator suffers from a deficiency which is common to Fourier-based estimators: undesirably high variance. The variance of PHAN(fA) is examined in a later section.5 4.2.2 Pre-whitening The periodogram (Figure 4.2) and, to a lesser degree, the modified periodogram (Figure 4.5) suffer from spectral leakage in the power-law case, especially for spectra with 2 < / < 3. Fox and Hayes [22] and Gilbert and Malinverno [26] avoid the leakage problem by using pre-whitening procedures in which they modify the surface height data to have a relatively flat spectrum, perform spectral estimation, and correct the estimate according to the original modification. We can derive the simplest pre-whitening approach similar to that used by Gilbert and Malinverno as follows: If the power spectral density of the surface height is defined as [64, p. 270] 1 T/2 2 Sz(f) = lin -E [ /2 Z(x) exp(-j27rfx) dx, (4.16) oo T T/2 `(4.16) then the derivative of the surface height dZ(x)/dx will have a related spectrum: 1 TI />/2 dZ(x} 2 Sz'(f) lim 1E [T2 d ) exp(-j27rfx) dx T —oo T J-T/2 dx 5The Blackman-Tukey Spectral Estimator [9], PBT, was also examined. For processes other than white noise, the BT estimator decreases variance but increases bias as compared to the periodogram (Kay [20], p. 80). PBT was not included in this comparison because it is even more susceptible to spectral leakage (and therefore, slope insensitivity) than PPER. The expected value of PBT is equal to the expected value of PPER convolved with another window, as derived in Kay [20] (p. 98):.1/2 E[PBT(fa)] = / w(fA - )E[PpER(~)]d< J-1/2 Since the BT estimator has an expected value which is the result of the convolution of the true spectrum with two window functions, its value at a given frequency can be greatly corrupted by leakage from lower frequencies.

55 T= lim E j27rf T/ Z(x)exp(-j27rfx) dx T-oo T J-T/2 = 47r2f2Sz(f) = 47r2cIf -+2 Sz'(f) = c'lfl-J, (4.17) and we see that the spectral exponent #/' of the derivative process takes on values between -1 and 1 for 1 < #/ < 3. The spectrum of the derivative process is more nearly flat (it is white noise for /' = 0). Since there is no pronounced peak at low frequencies, spectral leakage is less of a problem. Gilbert and Malinverno approximate the derivative by taking first differences of the sampled data: Z'[xi] f Z[xi+l] - Z[x]. They apply a Hanning window and use a periodogram to obtain a spectral estimate of Sz'A(f^). After averaging multiple periodograms, an estimate of the spectral exponent /3 is obtained from an estimate of /': d =- '+2. (4.18) This procedure is now examined for the limited data case in which multiple periodograms are not available for averaging. The Hanning window is omitted for simplicity. Define a pre-whitened spectral estimate Ppw(fA) in terms of the first differences of the sampled data: IN-1 2 PPw(f^) = {Z[Xn+l] - Z[xn]}exp(-j27rfAn). (4.19) n=O The expected value of the pre-whitened estimator may be derived by a procedure similar to the derivation of (4.7): A 1/2 E Pw(fa = 2 1/2 WT(fA - 0)SzA(E) [1 - cos(27r)] d]. (4.20) [ J-1/2

56 10'2 P.8 -- =0.4 E.-....... 1=-0.4 ---- '=0.8., 10-3 I- -. - P '=-0.0 _ E! ' ---- 0.0 —. —8=0. 104. -.1.0........................ 1~ 10 10 ' 0.01 0.1 normalized spatial frequency, dimensionless Figure 4.6: Expected values of Ppw(f^) for five different 257-point power-law profiles. The profiles were zero-padded to 16 384 points. The expected values of the pre-whitened estimator are calculated using (4.20) and the modified power-law spectrum defined in (4.15) and are shown in Figure 4.6. These estimates are based on 257-point profiles. Power-law functions are fit to the linear portions of these estimate expectations (fA < 0.2), and the parameters c' and f' are converted to c and,3 using (4.18) and c = c'/(47r2). The results, given in Table 4.2.1, show that the expected value of the spectral slope is reasonably accurate because the spectral leakage has been largely eliminated. A Hanning window may still be useful in practice (since one does not know a priori that a spectrum is power-law), but we see that leakage is not a problem for the pre-whitened estimator in the power-law case. The variance of Ppw(fA) will be examined in a later section. 4.2.3 Capon's estimator An alternative procedure which is more direct than the pre-whitened periodogram is the use of Capon's estimator [13] (the so-called "minimum variance spectral estimator") described in Chapter 11 of Kay [32]. This estimator was developed originally for geological signal processing problems in which the number of sensors (and there

57 fore the number of spatial samples) were extremely limited and posed restrictions on what could be inferred from standard spectral analysis. Capon's estimator essentially customizes a filter at each frequency to minimize the total power output, subject to the constraint that the gain at the frequency of interest is unity. Therefore, the filter may be asymmetric in side lobe level according to the shape of the signal spectrum; in the power-law case, the side lobes are adjusted to reduce the leakage from low-frequency components. Capon's estimator, PcAP(fA), is obtained by first calculating an estimate of the autocorrelation matrix Rzz, whose elements are defined as [Rzz]ij = E [Z*[n]Z[n + i - j]]. (4.21) We can estimate the elements of the autocorrelation matrix using the modified covariance method, described in [32]: [N-1 N-1-p [Rzz]ij = 2(N1- p)L I Z[n - i]Z[n - j] + Z[n + i]Z[n + j]. (4.22) ( P) Ln=p n=O Capon's estimator is then given by PCAP(fA) = -t (4.23) where p is the dimension of the autocorrelation matrix and e = [1 ej2rjfA ej4rfA... ej27r(p-i)fJ (4.24) Since (4.22) is an unbiased6 and consistent estimator of the covariance matrix, we can evaluate the expected value of Capon's estimator by substituting the exact covariance matrix (obtained from the exact theoretical spectrum) into (4.23). This 6(4.22) is approximately but not strictly unbiased because all lags are not weighted equally and overlapping lags are used. A variety of covariance estimators have been examined, and (4.22) was found to be the preferred estimator for short data records. See Kay [32] or other texts for more information.

58 E 10~ 1- 1.2 - -=2.4 -------- P=1.6 P=2.8 -4 10 0.01 0.1 normalized spatial frequency, dimensionless Figure 4.7: Expected values of PcAP(f^) for synthetic profiles with various 1l using an autocorrelation matrix of dimension p = 70. Estimates are evaluated at the same 16 384 frequencies used in Figures 4.2, 4.5, and 4.6. procedure was used to calculate the expected values of PCAP(feA) for five values of f/ using a covariance matrix of dimension p = 70. The expected values of these estimates (Figure 4.7) have good agreement with the exact spectra over a wide range of spectral slopes. Parameters of power-law functions fit to the expected values of PCAP(fA) for fA > 1/(2p) are compared with the exact spectral parameters in Table 4.2.1. The bias of Capon's estimator is independent of N because (4.22) is an unbiased estimator of Rzz. However, the bias does depend on p (pA is the longest lag for which the covariance is estimated). For a given data set, an increase in p will result in reduced bias at the cost of increased variance. Thus N indirectly influences the bias because the range of p is constrained by N. I generally found p e 0.3N to be a good compromise between variance and bias. In the power-law case, the bias of Capon's estimator was found to increase noticeably for spatial frequencies corresponding to wavelengths much longer than pA. Calculated values of PcAp(f&) for spatial frequencies below l/(2p) were therefore discarded.

59 4.2.4 Variance comparison Capon's estimator, the modified periodogram with Hanning window, and the pre — whitened periodogram all produce parameter estimates whose expected values closely approximate the exact values for a power-law spectrum. It is instructive to compare the performance of these estimators in terms of variance. The exact expression for the variance of the modified periodogram involves higher-order moments and is usually' evaluated only for a Gaussian white-noise case [61]. The statistical properties of Capon's estimator for time series data are not known [32], although some results have been derived for array data. To avoid these problems and test these estimators under actual-use conditions, short synthetic topographic profiles of known statistics were generated using a spectral synthesis algorithm. Parameters obtained from Capon's, the modified periodogram, and the pre-whitened periodogram estimates computed from the synthetic profiles are then compared with those of the spectra used to create the profiles. The spectral synthesis algorithm closely follows [51] and consists of the following steps: (i) Generate a set of discrete Fourier amplitudes which, when squared and multiplied by 1/N, satisfy the desired power law. While [51] used a pure power-law spectrum, this study uses the modified spectrum (4.15) to assure a zero-mean surface height. (ii) Multiply the amplitudes by a Gaussian random variable such that the mean value of the amplitude satisfies the power law. (iii) Generate a set of complex discrete Fourier coefficients (the FFT of a real surface) using the randomized amplitudes and a uniformly distributed random phase, enforcing symmetry conditions such that the inverse FFT will be real-valued. (iv) Calculate the inverse FFT of the coefficients, resulting in a synthetic surface profile. Using this algorithm, 10 synthetic 64-point surface profiles were generated with

60 10 A Z100 ~i0-^^^'^^^ -' 0I Hanning "* ^ - ' Pre-whitened1;$;*, r ^ ' 0.1 normalized spatial frequency, dimensionless Figure 4.8: Capon's, modified periodogram, and pre-whitened periodogram spectral estimates for ten different power-law profiles. The profiles were 64 points long. The modified periodogram and Capon's estimates have been offset for clarity. The solid lines represent the exact spectrum. /3 = 2.4 and cA, a, and oa chosen as before. I then calculated (i) Capon's estimates using a covariance matrix of dimension p = 20, (ii) modified periodogram estimates using a Hanning window, and (iii) periodogram estimates based on the first differences of the surface profiles. We see in Figure 4.8 that the modified periodogram and pre-whitened periodogram estimates have noticeably greater variance. Estimation of power-law parameters c and /3 from the spectral estimates is an estimation problem in itself, independent of the spectral estimation problem described thus far. In this analysis, estimates of c and /3 are obtained by performing a minimum absolute deviation fit of a power-law function to the estimated values of the power density spectrum. The minimum absolute deviation fit, which is performed in log-log space, is more robust than a least-squares fit. To make the power-law fit independent of the frequency sampling of the spectral estimates, the estimates were evaluated at frequencies spaced very closely together so that the estimator would be a smooth curve between the sampled points. Values of the spectral estimates outside their regions of high accuracy (4/N < fA < 0.5 for PHAN, 1/N < fA < 0.2 for Ppw, and l/(2p) < fA < 0.5 for PCAP) were

61 A Estimator mean cA j std. dev. ca mean / std. dev. p PHAN 0.724 x 10-6 0.819 x 10-6 2.767 0.631 Ppw 1.310 x 10-6 1.690 x 10-6 2.469 0.567 PCAP 0.891 x 10-6 0.433 x 10-6 2.371 0.313 Table4.1: Spectral statistics derived from fits to modified periodograms, periodograms based on pre-whitened data, and Capon's estimators of ten synthetic 64-point profiles with ca = 1.0 x 10-6 and / = 2.4. discarded as before. The mean and standard deviation of the estimated roughness amplitude ca and spectral slope /3 are given in Table 4.1. We see that the spectral parameters predicted using Capon's estimator have roughly half the standard deviation or one-quarter the variance of estimates derived from either the modified periodogram with the Hanning window or the periodogram based on pre-whitened data. Capon's estimator is therefore selected as the preferred estimation algorithm for use with short data records in the power-law spectrum case. 4.2.5 Formulas for real frequency Expressions for the spectral estimators listed in the preceding sections were given as functions of the normalized frequency fA that is dependent on the sampling interval of a measured surface profile. In remote sensing studies, estimates of the surface spectrum in terms of real spatial frequency f (in m-1) are necessary so that surface features may be compared to the electromagnetic wavelength and so that spectral estimates derived from profiles with different sampling intervals A may be compared or combined into a composite spectrum. The following expressions correspond to (4.2), (4.7), (4.8), (4.4), (4.10), (4.11),

62 (4.12), (4.13), (4.14), (4.19), (4.20), (4.23), and (4.24), respectively, for the case of real (non-normalized) frequency f, where f satisfies -1/(2A) < f < 1/(2A): PPER(f) = Z[x,z] exp(-j27r f An), (4.25) E [PPER(f)] =A /2A WT(f - )Sz() d, (4.26) 1(sin WT(f) = si n) (4.27) AN-1 2 PHAN(f) = NU Z Z[xn]wH[n] exp(-j2rf An), (4.28) n=O J 1/2A E [PHAN(f)] = A 2_ WH (f- )Sz() d, (4.29) WH(f)= 4NU (b+ - b2 + b3 (4.30) = sin(rfAN) (4.31) b1 =. ( fA^ 1 (4-31) sin(7rfA) sin{r[f - (N7i)]AN} b2= sin[f (N-l)A]A} (4.32) sin{ir[f+ (Ni)A]AN} (4.) b3 = { (N1) (4.33) sin{ir[f~ (N-i)A]A} Ppw() = - { [Xn+l] - Z[] } exp(-j2,rf An) (4.34) 7Vwn=0 f - o E [P (f)] = 2A /2 WT(f - )Sz() [1 - cos(2'7r)] d, (4.35) PCAP(f) = - ' (4-36) Te = e=[1 ej2r'fA' j41rfA... ej2i'(p-I)-A 1(437)

63 -E -- spectral estimate - - real signal E 10'9 -... --- aliased high freq. component.. x discarded points _.=L |I Nyquist frequency| 10.101~..... --- —-— \ 2 3 4 5 6 7 8 9 2 0.1 1 spatial frequency, m' Figure 4.9: Aliasing in the power-law spectrum case. Points marked (x) are discarded to reduce errors due to aliasing. 4.2.6 Aliasing The estimates presented in the previous sections assumed that aliasing was not present, and the synthetic surface profiles were constructed to be bandlimited; i.e., they had no spatial frequency components above the Nyquist frequency f, = 1/(2A). Real surfaces will generally have frequency components at spatial frequencies above f,. As a result of a necessarily insufficient sampling rate, power at frequencies higher than fc will be aliased into the range -fc < f < f,. We can derive an equation showing this result in terms of real frequency by extending the derivation of Oppenheim and Schafer [47, pp. 26-28], obtaining 00(4.38) SZaliased(f) = E Sz(f + ). (4.38) n=-oo In the power-law case, the power spectral density falls off quickly with increasing frequency. The n = 0 (exact) and n = 1 (first alias) terms and their sum are plotted for a simulated power-law case in Figure 4.9. We see that the effects of aliasing are most significant at the high-frequency end of the spectral estimate. As a first-order correction to avoid the aliasing problem, values of the spectral estimate computed for frequencies greater than fc/2 may be discarded, as indicated in the figure.

64 4.2.7 Applicability to real data Real geological data possess a variety of irregularities and surface statistics. Spectra of natural topography clearly cannot conform to a power-law model at all spatial frequencies: roughness features on the scale of millions of meters are unphysical, as are those at subatomic scales. However, in the studies cited in Chapter II as well as in the present investigation, measured spectra of topography are well modeled by a power-law spectrum in the form of (2.8) over some range of spatial frequencies. The comparison of spectral estimators here has been based on an idealized spectrum with constant slope (except for the very low frequency roll-off introduced in (4.15) to avoid an unphysical singularity) in order to have an easy reference for comparison. Because both the Fourier estimators with windowing and Capon's estimator are local-neighborhood estimators, this comparison also illustrates the relative performance of these estimators for spectra which have variations in slope with frequency or local perturbations in level that violate monotonicity. 4.3 Estimation of Two-Dimensional Spectra While the previous techniques allow the estimation of the spectrum of a linear profile, we would in general prefer a two-dimensional spectrum of surface topography. Such a spectrum might show hidden anisotropy or may indicate that the surface is isotropic in the statistical sense. In either case, the full 2D spectrum provides valuable additional information. Two-dimensional spectral analyses are precluded when 2D data are not available. Spectral studies of seafloor morphology, for example, usually utilize profiles collected by bathymetric instruments which are towed by a ship, producing a seafloor profile along the ship's path. Some directional insight can be obtained by obtaining profiles

65 in orthogonal directions. On land, 2D data are somewhat more accessible. Huang and Turcotte [30, 31] and England [21] use digital elevation models to estimate surface spectra, but their methods average over azimuth angle in the spectral domain. Two-dimensional data at scales of centimeters to tens of meters can be tedious and costly to obtain, necessitating smaller dimensions of a sampling grid, i.e., fewer than 100 x 100 measured points per 2D profile. 4.3.1 Two-dimensional estimators Two-dimensional Fourier-based spectral estimators such as the 2D periodogram suffer from the same problems as the corresponding ID estimators: leakage in the power-law case and high variance. A 2D Capon's estimator is one solution to these problems. The ID Capon's estimator is readily extended to the two-dimensional case [32, ~15.8], at the cost of increased complexity. Inversion of the covariance matrix may become problematic because the matrix expands from dimension p x p to p2 x p2. For a similar level of resolution, the 2D Capon's estimator will require a much greater volume of data (as will the other 2D estimators). If such high quantities of data are not available, as is often the case for 2D data sets, the variance of the Capon's estimate will be correspondingly higher. 4.3.2 2D spectral estimation of isotropic surfaces If a surface area is known to have or can be assumed to have isotropic statistics (from knowledge of the origin of the surface, or perhaps from measured 1D spectra in several directions), then an estimate of the 2D roughness spectrum can be obtained from estimates of the ID spectrum. These ID spectral estimates may be calculated

66 from individual rows of the 2D surface elevation grid; we can calculate a spectral estimate for several or all rows and then average the estimates together. The resulting estimate will have a variance which is much lower than that of a single row estimate even though the grid rows are not completely independent. Suppose that an isotropic surface has a 2D power spectral density Sz(f, fy) which has the form of a power law down to some low spatial frequency flom and some unspecified finite form below that: af7-', f,> flo~ Sz(f, fy) = Sz(f ) = f f(4.39) finite, fr < flow where the radial spatial frequency fr = (f2 + fy)1/2. We now need to relate the parameters a and -y to the parameters of the ID spectrum of surface heights measured along a straight line on the same surface. This ID spectrum will also have the form of a power law for f > flow and can therefore be described by the parameters c and To obtain this relation, begin with the general expression for the Wiener-Khintchine relation in terms of the Fourier transform and write Rz(Tr, Ty) for the surface correlation E [Z(x, y)Z(x + r, y + Ty)]: Rz(rT, y) = j / Sz(f., f) exp(j2ir[frZ + fyTy]) dfi dfy. (4.40) -00 J-00 Changing to radial coordinates (rx = rT cos Tr, Ty = rr sin o) and noting that Sz(fr, fe) is independent of fo, we obtain Rz(r, re) = Sz(fr)fr f exp(j2rf7r, cos () d( df,, (4.41) where ( = fe - To, leading to ) = S(f)Jo(2f)f df (4.42) Rz(Tr) = 27 Sz(fr)Jo(2W7fr)fr dfr, (4.42)

67 where Jo is the zero-order Bessel function. This expression differs from a similar (but incorrect) expression given by Voss [60, equation (1.52)]. The ID and 2D roughness spectra of an isotropic random field are related by an Abel transform [28]. This relation is more easily visualized using the surface autocorrelation functions (the inverse Fourier transforms of the ID and 2D surface spectra). The ID autocorrelation function is a slice of the 2D autocorrelation function: SzlD(fx/) = / / RZ(Tr, rTy) 6(ry) e-j2r(fxrT+f f d r) dry. (4.43) 00oo J-oo Using Fourier transform properties:,00 SZ1D(fM) = Sz(f fy) * (f,) = Sz(f, fy) df4. (4.44) Let Sz(fx, fy) be given by (4.39), and let fx be restricted to the power-law region, i.e. If| l > flow: SZlD(fr) = o [(f)' + (fy)2]/ dfy. (4.45) Changing variables, we have SZ1D(fx) = 2a 0 f — dfr, (4.46) which can be integrated analytically [4, p. 201]: a,%WFr F "I — 1 Sz1D(f) = 2) If (4.47) where r(x) is the gamma function. Therefore, if SzlD(fx) = clfxl- for If,[ > flow, the ID and 2D spectral parameters are related by? = + 1 (4.48) a = (-y(2) c. (4.49)

68 Equation (4.48) agrees with the result derived by Voss; he does not give an expression comparable to (4.49). As discussed in Chapter II, ID profiles of a fractal surface with 1 < Df < 2 have power-law spectra with exponent / given by /3 = 5 - 2Df, so that 3 > / > 1. From (4.48), the 2D power-law exponent 7y satisfies 4 > 7 > 2, corresponding to 2D surfaces with 2 < Df < 3, where D/ = (8- -)/2. The performance of the linear, averaged Capon's estimator in the determination of the 2D power-law parameters a and y is examined next. A variation of the 2D spectral synthesis algorithm described by [51] was used to generate 2048 x 2048 -point synthetic topographic surfaces with spectral exponents of 3.0, 3.4, and 3.6 with sampling intervals A of 1 cm. The exact surface spectrum was modeled as isotropic, with a form similar to (4.15): (a/a2)frexp(-fr/[2{2]), fr < fiou Sz(fr) aft 7i fow fr < fc, (4-50) 0, fr > fc where fr is the radial spatial frequency, a = 2.384 x 10-11, flow = 0.01 m-1, and a and oa are chosen as before. The cutoff frequency, fc = 1/2A = 50 m-1, is equal to the Nyquist frequency along the fx and fy axes. As noted previously, arrays of data containing over a million points are rarely available. In practice, in situ surface height measurements yield data sets which are much smaller, say, 40 x 40 points. Since such a data set is useful over a relatively narrow band of frequencies, several sparse grids may be collected at different scales, a spectral estimate calculated for each, and a composite spectrum formed from the individual spectral estimates. Scale factors are particularly important in the calculation of these multi-scale spectra.

69 0 5 10 15 20 25 3.0 35 E 2.5 E 2.0 - c 1.5 -o 1.0 -co 0.5 0.0 1 0 5 10 15 20 25 30 35 distance, m Figure 4.10: 41 x 41-point sample of a 2048 x 2048-point synthetic power-law surface E 10scales wh power-law fit 1 l 1 1- A l cm c 10 -A =4cm 0-12 a10 ___ 0.1 1 10 spatial frequency, aim Figure 4.11: Averaged Capon's estimates from 2D synthetic surface height data at three scales with power-law fit. Each 2048 x 2048-point synthetic surface is sampled at sampling intervals of' 50 cm, 10 cm, and 4 cm, yielding three 41 x 41-point data grids. (A sample data grid is shown in Figure 4.10.) At each scale, a ID Capon's estimate (with p = 12) is calculated for each grid row at frequencies l/(2pL) < f < 1/(4A), and the estimates are averaged over the rows. A power-law function is then fit to the averaged spectral estimates at the three scales. (A sample fit is shown in Figure 4.11.) The upper frequency limit of 1/(4A) was chosen to reduce the effect of aliasing. Estimates of ID power-law parameters c and 13 are then converted into estimates of a and -y using (4.48) and (4.49). This process was performed on ten synthetic surfaces at each of three spectral

70 Surface (exact) Estimated parameters a 7 mean a std. dev. a mean ] std. dev. y 2.384 x 10-11 3.0 2.34 x 10-11 2.83 x 10-12 2.971 0.041 2.384 x 10-11 3.4 2.28 x 10O-11 1.55 x 1012 3.386 0.041 2.384 x 10-11 3.6 2.17 x o10-11 1.14 x 10-12 3.555 0.061 Table 4.2: Power-law parameters derived from fits to averaged Capon's estimates at three scales. exponents. The mean and standard deviation of a and 5 are listed in Table 4.2. We see that the mean a and ' are quite close to the exact values. 4.4 Roughness Spectra of Mount St. Helens Debris Flows The linear, row-averaged Capon's estimator has been shown to be the preferred method for obtaining roughness spectra from data sets like those collected at Mount St. Helens. Before applying the estimation algorithm to the surface elevation profiles, it is appropriate to use this estimator to re-examine the median interpolation procedure used in the last chapter to replace holes in the surface profiles due to filtered outliers. 4.4.1 Verification of median interpolation The effects of median interpolation on spectral estimates were tested using three synthetic 2D surface profiles generated using the spectral synthesis technique described in the last section, with a = 2.384 x 10-11 and y = 3.0. The synthetic profiles were sampled at three sampling rates. Holes were introduced into the grid samples in positions corresponding to the holes in the measured debris flow surface

71 Exact spectra 4.768 x 10-11 2.0 Original sampled profiles: mean 1.974 x 10-11 2.013 std. dev. 0.232 x 10-1 0.040 Interpolated profiles: mean 1.693 x 10-11 2.006 std. dev. 0.129 x 10-11 0.032 Table 4.3: Verification of median interpolation. profiles. These holes were then filled in using multiple passes of median interpolation, and the linear, row-averaged Capon's estimator was used to produce spectral estimates over three ranges of spatial frequency. The mean elevation was subtracted from each row before calculating the spectral estimate. A power-law function was then fit to the three spectral estimates for each surface. A second set of spectral estimates were obtained directly from the original grid samples for comparison to those arising from the interpolated profiles. The mean and standard deviation of c and,3 are compared in Table 4.3. Values derived from interpolated profiles are quite similar to those obtained from the original data. 4.4.2 Primary debris flow surfaces A Capon's spectral estimate was calculated for each row of each profilometer and survey profile using a covariance matrix of dimension p; 0.3N, where N is the length of the profile row. (The mean row height was subtracted before calculating the estimate.) The spectral estimates were then averaged over the rows of each grid. The high-frequency half of each averaged estimate was discarded to reduce aliasing

72 effects. Parameters used in preparing spectral estimates from each profile are listed in Table 4.4. Averaged estimates at a single scale or at multiple scales were then fit with a power-law function (2.8), where appropriate, using a minimum absolute deviation fit. Profile (ID) spectra fitting the power-law model were then converted to surface (2D) spectra using the formulation given in section 4.3, in which the 1D and 2D roughness spectra are related through their corresponding autocorrelation functions and the 2D autocorrelation function is assumed isotropic. The first terrains examined were at the Johnston Ridge site, a debris flow area adjacent to and due west of Johnston Ridge (see Figure 3.1). This debris flow had areas exhibiting large variations in elevation (e.g., 8 m in a 32 m x 32 m survey grid). The roughest surfaces that we measured in this area were the primary debris flow surfaces-surfaces formed during the landslide-avalanche and modified through the removal of material by wind and water. Spectra from the roughest of these are shown in Figure 4.12(a). The survey spectrum is the curve extending from 0.05 to 0.25 m-l; the higher frequency spectra were computed from profilometer scans. The linear character of the A = 1 m, A = 2 cm, and A = 5 mm spectral estimates suggests that a power-law spectrum is appropriate for this surface over some range of spatial frequencies. The equation of the power-law fit to these three estimates is Sz(f) = (3.21 x 10-4)fl-2 34, (4.51) where the hat denotes an estimate. This fit is shown as a dashed line in Figure 4.12(a). While we have no evidence that the surface spectrum follows this power-law at the intermediate spatial frequencies, it does seem significant that both the meter- and centimeter-scale spectra are well fit by the same power law. The spectral estimate in Figure 4.12(a) derived from the A = 2 mm profile seems

73 data set N I p 1/(2pA) 1/(4A) Now, Nhigh scan 2b 51 0.002 16 15.625 125 16 128 scan 2c 41 0.01 13 3.84615 25 20 128 scan 3a 46 0.02 14 1.78571 12.5 19 128 scan 3b 41 0.005 13 7.69231 50 20 128 scan 3c 41 0.002 13 19.23077 125 20 128 scan 4a 55 0.005 17 5.88235 50 16 128 scan 4c 51 0.002 16 15.625 125 16 128 scan 5a 71 0.005 22 4.54545 50 12 128 scan 5b 51 0.002 16 15.625 125 16 128 scan 6a 51 0.002 16 15.625 125 16 128 scan 6b 41 0.01 13 3.84615 25 20 128 survey 1 33 1 10 0.05 0.25 26 128 survey 2 9 4 3 0.04167 0.0625 86 128 survey 3 9 4 3 0.04167 0.0625 86 128 survey 4 17 2 6 0.04167 0.125 43 128 Table 4.4: Spectral estimation parameters: N is the dimension of the sample grid, A is the sampling interval in meters, p is the dimension of the autocorrelation matrix, l/(2pA) and 1/(4A) are the low- and high-frequency limits for the spectral estimator in m-1, and Nl,, and Nhigh are the indices of the low and high estimator data points used.

74 E E -.? U) 10 CD t a 10'1 A m 10'2 10-3 104 105 A =2cm 10'6- EDM noise floor A =5mm 10-8 108 A =2 mm 0.1 1 10 100 spatial frequency (m1) (a), A ----------------------------------------------—. -- 10 E. — cmJ a) -4 *o 10. 10'6 0 a 10's A = 1 m.. EDM noise floor..A = c im 0.1...... I........ I 1 10 spatial frequency (m'1) (b) 100 E r'CO C -a ). 4 -CO Q. 100 10-2 10-4 106 10-8 N- A=2m EDM noise floor A=5 mm............... -.. ~~............................... 0.1 0.1 1 10 spatial frequency (m'1) 100 I00 (c) Figure 4.12: Spectra of primary debris flow surfaces: (a) JRS, survey 1, scans 3a, 3b, and 3c; (b) JRS, survey 1, scan 2c; (c) ERS, survey 4, scan 5a.

75 inconsistent with the A = 5 mm estimate for the range of spatial frequencies where both apply (19-50 m-1). This is true for all spectral estimates based on scans with A = 2 mm. The time interval between EDM measurements was shortest for these scans because the EDM moved only a short distance. It is suspected that the high measurement rate adversely affected the elevation measurements and resultant spectra. The blurred appearance of the A = 2 mm grayscale images in Appendix A supports this interpretation. Another possibility is that the laser spot size was measured inaccurately. Because of these uncertainties, I have omitted spectral estimates based on A = 2 mm scans from the remaining figures and from the spectral analysis. No claims will be made regarding the spectrum at frequencies estimable by the A = 2 mm scans alone. We can convert (4.51) to an estimate of the surface (2D) spectrum using (4.48) and (4.49), obtaining the following power law: SZ2D(fr) = (1.77 x 10-4)f73'34. (4.52) In Figure 4.12(b), the same survey spectrum is shown together with a profilometer spectrum derived from profilometer scan 2c, which was collected on a similar flow unit adjacent to the one surveyed. A power-law fit to these estimates gives the following spectra: Sz(f) = (3.57 10-4)1f-231 (4.53) SZ2D(fr) = (1.95 x 10-4)f331. (4.54) The other group of terrain measurements were at the Elk Rock site, in the North Fork Toutle River valley due west of Elk Rock. The valley was a mixture of hummocky debris mounds scattered across a flat, sedimented plain. A profilometer scan from one of the debris mounds yields the spectral estimate shown at the

76 high-frequency end of Figure 4.12(c) (A = 5 mm, scan 5a). The figure also shows an estimate from a survey grid (survey 4) containing both terrain types. (The debris mounds were typically smaller than the 32 m x 32 m survey area.) A power-law fit may be appropriate for the A = 2 m and A = 5 mm spectra here as well; the fit is shown in the figure. Equations for this fit are Sz(f) = (4.98 x 10-0) If j25 (4.55) SZ2D(fr) = (1.95 x 10V4)fr31. (4.56) 4.4.3 EDM noise floor While investigating the problems with the A = 2 mm scans, it became necessary to determine the minimum spectral amplitude detectable by the EDM using the selected estimation algorithm. The minimum detectable spectral magnitude was determined by generating spatial domain traces having single-frequency roughness (i.e., a sinusoidal surface profile) with amplitude equal to the minimum vertical resolution of the EDM (3 mm). Traces were generated for single-frequency roughness components of 1, 4, 10, and 40 m-1. Spectral estimates of these four surfaces are shown in Figure 4.13. Each spectral estimate has a single-frequency peak; the peaks trace out a linear equivalent noise floor for the EDM. Portions of a surface roughness spectrum that lie below this level will not be observable above the spurious spectrum of noise due to the limited precision of the EDM. 4.4.4 Mudflow- or water-deposited surfaces The smoother surfaces formed by the original mudflows or by water-deposited sediments had spectra that were markedly different from those of the primary surfaces. For example, in Figure 4.14(a), we see that the spectral estimate from a scan

77,-:,.......... I/ o)" " " ".,-.A A r,0-8 j — = -simulated noise, 40 m'1 | V\ - simulated noise 4 Q)............................ ~,.- ',,~ ' 10e '.......... simulated noise, 40 m'l,'.. v. D- - simulated noise, 1 m1 ~. -1: - EDM noise floor I I I gill I-I I I I I' ' I" 1 I I I I 1gil I I I I 'I''11 I 0.1 1 10 100 spatial frequency (m') Figure 4.13: The noise floor of the EDM was determined by calculating spectral estimates of surfaces with single-frequency roughness. (scan 6a) of such a surface (at the Elk Rock site) lies on the EDM noise floor discussed in the previous section. The true spectrum of this surface must lie below this level and therefore cannot form a power-law trend with the survey spectrum of the local mixed terrain (A = 2 m from survey 4 in the figure). This difference in spectral behavior is consistent with the visible character of the sedimented surface, which appeared quite smooth and flat. Similar remarks apply to the profilometer spectrum of a mudflow- or water-deposited surface at the Johnston Ridge site, which is shown in Figure 4.14(b) along with the survey spectrum of the surrounding primary debris flow. Spectral estimates derived from two low-resolution (A = 4 m) surveys on sedimented (survey 3) and mixed (survey 2) surfaces at the Johnston Ridge site are shown in Figure 4.15 together with the primary debris flow survey from Johnston Ridge and the mixed, lumpy terrain survey at the Elk Rock site. These spectral estimates are very short due to the small number of sample points (9 x 9 in each). The mixed terrain spectrum is similar to the primary debris flow spectrum at low spatial frequencies, while the sedimented surface has a much lower roughness.

78 -..0 - 10 -Or E cE 10 c - 4 0 10' -a 0) Ca 10 0 10 NA = 2 m EDM noise floor A= 1A-......................11...I. 0.1 0.1 1 10 spatial frequency (m'1) 100 (a) E E 0. n. sn *~ s 10'1 10-2 10'3 10-4 10'5 106 10-7 10-8 A=1m EDM noise floor i 1 I Wii l i i Ii i i I f1i i i wil! ^^ - - l A=5m. a, 1.. 0.1 1 10 spatial frequency (m"1) 100 (b) Figure 4.14: Spectra from scans of smooth surfaces shown together with survey spectra of surrounding areas: (a) ERS, survey 4, scan 6b; (b) JRS, survey 1, scan 4a.

79 E 10 Lumpy terrain v. Sedimented 0) 10'2 Primary debris flow. 5 6 7 8 9 2 0.1 spatial frequency (mr1) Figure 4.15: Survey spectra 4.5 Summary The preceding analysis shows how leakage causes some Fourier-based estimators to yield an estimated spectral exponent of 2.0 for surface profiles having power-law spectra with spectral exponents between 2 and 3. This may explain the number of studies reporting power-law exponents of 2.0 for natural terrains. It has also been shown that Capon's estimator may be used to avoid the leakage problem and measure the surface spectrum more accurately. In addition, Capon's estimator has reduced variance, which is useful when surface data records are short, as is often the case in remote-sensing studies. For surfaces which are known or can be assumed to be isotropic, a linear, averaged Capon's estimator can be used to estimate the 2D power-law parameters. This procedure was used on the debris flow profiles collected at Mount St. Helens. The rough surfaces associated with the 10-year-old debris flows possess spectra that are well modeled by a power-law function over a range of spatial frequencies. These measured spectra will drive the design of artificial power-law surfaces, which is the subject of the next chapter.

CHAPTER V DESIGN AND MANUFACTURE OF ARTIFICIAL POWER-LAW SURFACES The spectral estimates in the previous chapter show that certain volcanic debris flow surfaces at Mount St. Helens belong to the class of power-law surfaces discussed in Chapter II. The cited references and the present study show that power-law surfaces are common and often of interest in remote sensing. A better understanding of electromagnetic scattering by these surfaces may be obtained through scattering measurements and application of rough surface scattering models. In situ scattering measurements of natural power-law surfaces are difficult to perform and even more difficult to interpret for purposes of model verification. Aside from common logistical problems (e.g., operation of radars in dusty environments, cost and time constraints for on-site work), there are several measurement uncertainties which interfere with model verification. Natural rough surfaces rarely conform to a single statistical description over extended areas due to the variety of processes that contribute to their formation. For land surfaces, volume scattering by subsurface inhomogeneities is usually present. Both of these factors yield effects that are difficult to decouple from those due to surface roughness. More important, interpretation of scattering measurements is hindered by the difficulty of obtaining an 80

81 accurate statistical description of surface roughness. As discussed in Chapters III and IV, estimation of the surface roughness spectrum is particularly problematic for power-law surfaces, requiring a large quantity of surface height data and specialized spectral estimation algorithms. The uncertainties just discussed may be largely avoided through the use of man — ufactured surface analogues. The analogues are homogeneous, conforming to a single specified statistical model. They are constructed from a homogeneous material to eliminate subsurface scattering, and their surface statistics are well known. Because the surface roughness is specified by the investigator, particular cases of roughness may be studied, or a single roughness parameter may be varied in order to study its effect on surface scattering. Finally, the scattering measurements may be performed under more controlled conditions, typically indoors, allowing greater accuracy and repeated measurements if necessary. This chapter describes the design and manufacture of artificial power-law surfaces. The design history, material selection, and surface generation algorithm are discussed, and the milling system is described. The milling procedure is explained, including time and disk space requirements, and difficulties in milling complex surfaces are discussed. The chapter ends with a description of surface verification efforts. 5.1 Analogue design history The original project plan called for the manufacture of power-law surfaces measuring one meter square. The surfaces were to be composed of a dielectric material with permittivity near that of soil. Scattering by the surfaces was to be measured using an X-band (9.75 GHz) radar and radars of other frequencies if available. Surface construction was planned to be similar to that used by Nance [44], who

82 presented scattering results from artificial Gaussian conducting surfaces [45]. As in his study, the surface fabrication was to be performed by a contractor; however, requirements for the present study were more sophisticated due to the more complex surface roughness (power-law instead of Gaussian) and the use of a homogeneous dielectric rather than structural foam that was to be coated with conducting paint. After an extensive search of machining companies in southeast Michigan, I was unable to find a contractor who was able to fabricate the surfaces for a manageable cost. One obstacle was the desired surface size: one meter square was too large for several companies' equipment. A more significant impediment was insufficient data capacity. The design specification for the artificial surfaces called for the replication of surface features on the scale of one-tenth or one-twentieth of the radar wavelength. A high number of surface points were therefore required to represent the surface for fabrication (the surfaces eventually required 1318 x 1318 points). This quantity of data was quite beyond the capabilities of numerous machining firms, several of whom specified their data capacities in feet of paper tape. One company was found in Brighton which would probably have been able to fabricate the surfaces. However, the eventual cost per surface and inconvenience of working with an out-of-town company, as well as the likelihood of needing several attempts to get a usable surface, led to the decision to purchase milling equipment and manufacture the surfaces at the Radiation Laboratory. 5.2 Dielectric materials Although a dielectric surface introduces greater complexity in construction and in backscatter modeling, the greater similarity to natural terrains led to the selection of a dielectric surface rather than a conducting surface. There were multiple conflicting

83 requirements for the dielectric material: * homogeneous (to avoid volume scattering) * high permittivity (to assure a measurable return signal) * lossy (to eliminate reflections from the bottom surface of the dielectric) * machinable using standard milling equipment * inexpensive * low density (to eliminate need for elaborate target mount or special handling) A variety of candidate materials were examined [55]. Machinable wax and mixtures of machinable wax and graphite were rejected as too expensive and insufficiently strong to support their weight (a one meter square slab might crack when being moved). A mixture of graphite, titanium dioxide, and silicone resin suggested by engineers at Texas Instruments was attempted but was never made to solidify. A mixture of plaster and conductive paint was rejected due to concerns about homogeneity and durability (cracking under tension). The material ultimately selected, ultra-high molecular weight polyethylene (UHMW), was chosen for its good machinability, homogeneity, and reasonable cost. The search for dielectric materials led to a change in the dimensions of the artificial surfaces. This change was principally due to mass and cost constraints. A 1 m x 1 m x 15 cm slab of these materials could measure from 130 to 260 kg, requiring a very strong (and probably highly scattering) surface mount and special support to prevent the material from cracking under its own weight. By changing the surface scale by a factor of 3.538, the X-band measurement could be simulated using a 35 GHz radar with a smaller and shallower target surface. Surface strength was less

84 A A 35 GHz radar UHMW slabs metal plate Figure 5.1: Configuration for measuring UHMW dielectric constant. critical at the new scale, and obtaining a measurable return was less a problem at 35 GHz than at X-band. The decision to reduce the surface scale was fortuitous in view of problems with the milling table that will be described later. A surface size of 18 inches (45.72 cm) square was selected based on availability of UHMW slabs. This size corresponds to an X-band surface measuring 63.7 in. (161.8 cm) square and is sufficiently large to allow two radar looks with minimal footprint overlap. Because UHMW polyethylene is not lossy, elimination of reflections from the bottom surface of the slabs was accomplished using the time-gating feature of the network analyzer that processes the radar signals. Several 6-inch slabs of UHMW were placed beneath the milled slab to move the bottom interface away from the milled surface. This configuration is described and illustrated in Chapter VI. The dielectric constant of the UHMW slabs was measured using the 35 GHz radar as shown in Figure 5.1. The radar was pointed vertically downward at a stack of UHMW slabs (two 6-inch and one 2-inch) and the radar bandwidth set to 3 GHz for maximum resolution. A metal plate was inserted at various levels in the stack to

85 identify the reflections from the top, bottom, and internal interfaces. Removing the plate, the time interval between the returns from the top and bottom was found to be 3.397 ns. The round-trip time for travel through free space at the same distance (0.355 m) is 2.37 ns; the ratio of these times gives the index of refraction, n = 1.433. The dielectric constant is then obtained directly: e, = n2 =2.05. (5.1) Polyethylene has extremely low loss in tabulated values below and above 35 GHz; it was believed that the loss tangent would be too small to measure at 35 GHz using techniques available in the Radiation Laboratory. The dielectric constant was therefore assumed to be completely real. 5.3 Surface Design Algorithm The surface height was specified over a 1318 x 1318-point grid with a point spacing of 1.23 mm (0.04843 in.) at the X-band scale or 0.3477 mm (0.01369 in.) at the 35 GHz (milling) scale. The surface profile was generated using a variation of the two-dimensional spectral synthesis technique described by Saupe [51] and explained as follows: The desired surface spectrum is the two-dimensional power-law spectrum: Sz(fr) = a f,-, (5.2) where fr = (f + fy2)1/2 is the radial spatial frequency. First, a set of 2048 x 2048 discrete Fourier amplitudes Zl(fx, fy) is generated such that 1 IA2Z(f,,)12 = Sz( ff) (5.3) N2A2 (5 3) for all points within the circle fr < f, where fc, the maximum spatial frequency, was (2A)-' = 406.5 m-1. Point spacing in the frequency domain was Af = (NA)-'1 =

86 0.39698 m-1. In other words, Z (fx fy) = V(fr)-12 (5.4) A dimension of 2048 is chosen because the FFT algorithm requires N to be an integer power of two, and 1024 is not large enough. The amplitudes Z1 are then multiplied by a randomness factor and a phase factor: Z(f, fA) = Z 1(f, Lf)Ge'*, (5.5) where G is a Gaussian random variable with zero mean and unity variance, and 1 is a phase random variable uniformly distributed over [0, 27r]. In this way, the mean value of the normalized square of the spectral amplitudes still satisfies (5.3). To ensure that the surfac heights are real-valued, the above procedure is used to generate half the surface amplitudes, and the other half are obtained by enforcing Hermitian symmetry: Z(ff ) = Z*(-f=,-f), (5.6) The randomized coefficien s Z(fxm, fy4) constitute the fast Fourier transform (FFT) of one realization of a rough surface random process having a roughness spectrum given by (5.2). To obtain the two-dimensional surface profile, an inverse FFT is performed on the coefficient array: I N-1 N-1 Z(xk, yI) = N2 E E Z(fm fyn) exp(-i27r[km + ln]/N). (5.7) m=O n=O The resultant synt synthetic surface array has 2048 x 2048 points; profiles were truncated to 1318 x 1318 points for surfaces in the present study. The surface roughness parameters, a and y in (5.2), were selected to be similar to those measured at Mount St. Helens and to be physically realizable on UHMW slabs using the milling machine. At the time the first artificial surface was designed,

87 the estimated two-dimensional roughness spectrum based upon survey site 1 and profilometer scans 3a and 3b had a = 1.52 x 10-4 and a = 3.29. (The surface spectral estimate was subsequently changed to a = 1.7672 x 10-4 and 7 = 3.345.) The UHMW slabs were 2 inches thick; reserving a minimum thickness of one-half inch for stability left a possible total relief of 1.5 inches. In order to obtain a surface profile having this amount of relief after scale reduction by a factor of 3.538, it was necessary to reduce the roughness amplitude, a, by a factor of 7. The final roughness spectrum used for the first surface analogue, UHMW1, had the mixed form of (4.50): (a/a )fr exp(-f/[2]), f < 0.01 m-1 Sz(fr) af= y 0.01 m1 fr < fc ' (5.8) O, fr > fc where a = 2.17143 x 10-5 7 = 3.29 a = 1.6438 a = 0.00482805 fc = 406.5 m-1 The parameters a and a were chosen as in Chapter IV by enforcing continuity of the spectrum and its first derivative. In effect, however, these parameters had no effect (except at frequency zero) due to the limited resolution of the spectral synthesis technique; the lowest non-zero frequency (equal to 1//NA) was 0.39698 m-1, which fell outside the Rayleigh region of the mixed spectrum. Spectral parameters for the second and third surface analogues were selected to span the range of spectral slopes (or alternatively, fractal dimensions) that natural

88 Surface UHMW1 UHMW2 UHMW3 a 2.17143 x 10-5 2.17143 x 10-5 2.17143 x 10-5 y 3.29 3.8 3.545 a 1.6438 19.8523 5.70356 a 0.00482805 0.00456435 0.00469065 Df 2.355 2.1 2.2275 Table 5.1: Spectral parameters (a, y, a, and a in (5.8)) used and fractal dimension Df of the surface analogues. terrains reasonably assume. The roughness amplitude, a, was held constant; the spectral slope, 7y, was set to 3.8 for UHMW2 and to an intermediate value of 3.545 for UHMW3. The values of the four spectral parameters and the associated fractal dimensions are summarized in Table 5.1. 5.4 Milling System The milling system was ordered in November 1991 and arrived in early 1992. The system has existed in two forms: the original configuration, as supplied by Techno, and the upgraded configuration made necessary due to reliability problems with the Techno components. 5.4.1 Origi al Milling System n The system as originally purchased (see Figure 5.2) consisted of a Techno threeaxis gantry positioning table, a Kress variable-speed electronic grinder, and a Techno MAC 100 programmable position controller, all purchased from Techno, New Hyde Park, New York; CNC Software's Mastercam CAD/CAM software, purchased from

89 PDS 486/33 Techno MAC 100 3-axis milling table computer x drive!-..,- x-axis motor serial port y drive y-axis motor z drive ---- z-axis motor I 1 Figure 5.2: The original milling system supplied by Techno. CIM Solutions of Canton, Michigan; and a Peer Data Systems 486/33 MHz DOScompatible computer, purchased from Peer Data Systems of Ann Arbor, Michigan. The components of the system are described in this section. The three-axis gantry positioning table holds the construction material in place and allows the milling tool to be moved to designated xyz coordinates. The table and x- and z-axis supports were constructed from extruded aluminum section. Precise motion in the x, y, and z directions was obtained via the rotation of ball screw assemblies by stepper motors. The ball screw mechanisms were covered by expanding vinyl bellows to protect them from dust and shavings. Limit switches at one end of each axis's range of motion provided a mechanism for zeroing the table. Table positioning was theoretically accurate to 0.01 mm (0.0004 inch), which was the linear motion obtained by rotating a stepper motor one step. This precision far exceeded the requirements for this project. The table was ordered in a special size, 48 x 49 inches, and with extra clearance (8 inches) for the z axis to allow milling a fairly thick slab measuring one meter square in a single piece. One problem that became apparent soon after starting work with the table was that the table sagged somewhat under its own weight due to support beams that were

90 insufficient for the table size. A test load of 227 kg (500 pounds) distributed over a 1 m2 area in the center caused a sag of 4.3 mm. To reduce the sag to acceptable levels, Techno provided an iron L-beam that was attached to table surface. The L-beam corrected the sag problem but reduced the usable table area. Fortunately, the desired surface size had been reduced by this time. The Kress electronic grinder rotates the milling tool, removing surface material as it travels along a programmed path. The grinder was attached to the z-axis and held an endmill (milling tool) in a collet of 3/8, 1/4, or 1/8-inch diameter. The rotation speed was adjustable from 8,000-20,000 rpm. The MAC 100 programmable position controller accepts toolpath coordinates from a post-processing program running on the 486 computer and translates these into properly timed currents in the windings of the stepper motors on the three axes. The MAC 100 controller was built by Techno, using some Techno-designed boards and other components built by Compumotor. The controller operated as an openloop system; i.e., there was no position feedback from the axes to the controller, so the controller had no way to verify that the stepper motors and axes had moved to specified positions. Open-loop control is not uncommon in stepper motor systems, but it is less reliable than closed-loop control, as will be discussed later. The CNC Software Mastercam Mill package is an extensive program for designing complex objects, producing machine toolpaths for the fabrication of these objects, and controlling the machining equipment while executing the machine toolpath instructions. (However, the Techno controller must be operated by a Techno postprocessor program and not from within Mastercam.) Because the surfaces in this study were extremely large and complex, and because they were produced by an external program on another computing platform, many of the features of the Mas

91 tercam package were not used. The principal use of the Mastercam package was to read a file of surface height coordinates and write files of toolpath coordinates for different milling tools. (Several sizes of tools were used in different stages of the milling process.) Techno sells a version of Mastercam designed for 286 computers; the faster and more sophisticated 386 version was needed in this study, so it was purchased from CIM Solutions. Both Mastercam and the post-processor program were run on the Peer Data Systems 486/33 MHz DOS-compatible computer. The computer had 16 megabytes of random access memory and a 210 megabyte hard disk drive, and was augmented with an Ethernet board to allow file transfer with the Unix workstations on campus. The communications link was necessary because Mastercam could process only a section of the artificial surface at a time. In the process, very large disk files were generated; it was not possible to retain all the surface files on the hard disk. 5.4.2 Failure of the Techno System The Techno system began having problems in October 1992. The z-axis motor either slipped or rotated in the wrong direction and seemed to be overheating during some test runs on wood surfaces. Upon inspection, the solder joints on the connector card inside the motor housing had heated up enough to melt the solder. Engineers at Techno were unable to determine a cause for the failure. In January 1993, it was proposed that the z-axis motor might be overheating due to excess weight on the z axis. The z axis was reverse-mounted on the Techno table to allow the full 8 inches of clearance; i.e., it was installed so that the movable carrier stayed fixed while the rest of the z-axis assembly moved up and down. To relieve part of the load on the z-axis motor, an 11.5-inch tension spring and support rod

92 11.5'0 z-axis assembly endmill Figure5.3: A load-bearing spring was installed on the z-axis to help support its weight. were installed as illustrated in Figure 5.3. The spring constant and position were selected such that the spring force balanced the weight of the z-axis assembly in its mid-position. The spring seemed to solve the slippage problems. Milling on the first UHMW slab commenced in February 1993. By early March, the machine seemed sufficiently reliable to leave it operating overnight, which was a great advantage, considering the very lengthy milling times. However, the machine began to malfunction again in late March: this time, the y axis became confused during a finishing run, resulting in a deep gouge in the slab. Manual commands to the mill were sporadically executed incorrectly. At this point, AC line noise was suspected, so an isolation transformer was installed. Machine operation was then normal again. In mid-April, the problem recurred. Both the y and z axes malfunctioned during an overnight roughing operation. The 1/4-inch endmill had milled down through the slab and through part of the table itself. The simultaneous failures of two axes led

93 Figure 5.4: The upgraded mill configuration. to the conclusion that the Techno MAC 100 controller was not reliable. This was a reasonable assumption in view of the other engineering problems with the Techno equipment. 5.4.3 Upgraded Milling System The upgraded milling system is shown in Figure 5.4. The Techno MAC 100 was replaced by a Compumotor AT6400 Indexer and a Compumotor SC30 Motherboard Rack, purchased from Amerinetics of Novi, Michigan. Optical encoders (part number MOD5641-25-100-L) were purchased from BEI Motion Systems of Plymouth, Michigan. Motion Architect software, included with the AT6400 from Compumotor, was used with Microsoft Visual Basic to write new control software. The indexer,

94 drive rack, and encoders arrived in late July 1993. The AT6400 indexer is the central controller of the upgraded system. The indexer receives ASCII-format move instructions from the control software, for example, "move x axis 100 steps right." The indexer then converts these instructions to step pulses and direction signals that are sent to the stepper motor drive card for each axis. These drive cards were removed from the MAC 100 but were manufactured by Compumotor; they were installed in the Compumotor SC30 Motherboard Rack for use in the upgraded system. The three stepper motor drivers use the step pulses and direction signals to drive the windings in the stepper motors in sequence, producing rotation in the specified direction. The optical encoders were installed to provide feedback to the indexer. They were mounted on the stepper motor cases using epoxy such that the motor shafts projected through the encoders. The encoders consist of two LED's and a rotating glass disc with 100 equally spaced lines. As the encoder turns, light from an LED is alternately blocked and then transmitted through the glass disc. The two detectors are in quadrature, resulting in 400 encoder pulses per revolution. The encoders are relative sensors-they do not provide absolute coordinates, so the indexer must keep a count of encoder pulses. The control software was written so that after each mill move, the encoder step count was compared to the stepper motor step count. If the counts disagreed by more than two steps, the program stopped, and the mill was deactivated. In addition to the encoder inputs, the indexer also had inputs from a home limit switch on each axis and a limit switch on the z axis that stopped the mill if the z axis went too low. Additional inputs from jog switches allowed manual control of the mill position. The indexer also had programmable outputs, one of which was used

95 to operate a relay that activated and deactivated the grinder motor. One other component of the MAC 100 was retained: the large transformer that converted 115 VAC down to 18 and 26 VAC to supply the motherboard rack and stepper motor drives. A transformer malfunction seemed unlikely. 5.5 Milling procedure 5.5.1 Generation of Toolpath Files Surface fabrication begins with the spectral synthesis procedure, described earlier in this chapter, which results in an 22.6 MB ASCII file of surface elevations arranged in a 1318 x 1318-point grid. Because the Mastercam software cannot accept input of more than 100 lines of ASCII surface elevations, the surface must be segmented into thirteen sections of 99 lines plus a fourteenth section of 31 lines. The FORTRAN program ds2mcinv.f performs this segmentation and applies the scale factor converting the surface from the X-band scale down to the 35 GHz scale. The program writes the surface sections as long lists of xyz coordinates. The resultant files measure about 5 MB each. The edges of the surface analogues were rounded by subtracting the random scaled height from a 1/4-inch radius along each edge, as shown in Figure 5.5. The rounding was performed so that no non-random sharp edges would be illuminated by the radar. Most rays striking the rounded edges should be scattered away from the radar. The FORTRAN program edgemax.f reads the surface file and writes a 3.8 MB file of edge coordinates. A couple of sections at a time are downloaded to the 486 computer. Mastercam is then used to convert the surface elevation files into toolpath files at two resolutions. Roughing files contain toolpaths for a 1/4-inch ball-end endmill. This tool is used to

96 reference surface 1/4 inch radius Figure5.5: Edges of the artificial surface are rounded by subtracting the surface elevation from a rounded reference surface with a 1/4-inch radius of curvature. remove the bulk of material at a rate much faster than the smaller tools. Finishing toolpaths are much more detailed (and much larger) and are used for the passes using the 1/16-inch and 1/32-inch ball-end endmills. Cross toolpaths are for the orthogonal finishing pass with the 1/32-inch endmill. The toolpath (. NCI) files are produced in several steps [54]. An ASCII file is input, and the resultant surface description is saved in Mastercam format as a geometry (.GE3) file. The next step is to chain the surface, i.e., to define the contours to be milled. This creates a chaining (.NCS) file. After chaining, which need be done only once per section, the tool parameters and desired resolution are entered, and Mastercam writes surface (.CDB), offset surface (P.CDB), and toolpath (.NCI) files, in addition to some miscellaneous small intermediate files. When all three toolpaths have been written, the roughing toolpath is filtered by removing points within a specified distance from the splines used to define the surface. This significantly

97 File type extension roughing finishing cross] geometry *. GE3 12 chaining *. NCS 12.5 12.5 12.5 surface *. CDB 0.4 3.2 3.2 offset *P.CDB 0.4 3.2 3.2 toolpath *.NCI 0.02 4.6 4.6 Table 5.2: Types and sizes (in MB) of Mastercam files. decreases the size and execution time of the roughing toolpath. Entering, processing, and transferring the generated files for a single section takes about 3 hours. The sizes of the various Mastercam files are given in Table 5.2. Some of these may be deleted (the offset surface and all but one of the chaining files), but the rest are uploaded to a Unix workstation and saved on magnetic tape, to avoid lengthy re-processing if a toolpath file is lost. 5.5.2 Milling machine operation The surface slab, which has been milled flat using a flat-end endmill, is mounted on the milling table after aligning its edges with the x and y table directions. It is held in place using aluminum angles tightened against each side of the slab. The mounting brackets are tightened carefully because re-aligning the slab would be extremely difficult if it came loose. In several instances, it was necessary to mill the aluminum brackets down a few millimeters in places where the surface was low near its edge. A Visual Basic program called mill controlled the entire milling process. Visual Basic was quite useful because it allowed easy construction of virtual control

98 panels and readouts describing the milling machine status. During operation, the program displays the current position of the mill (in steps and in inches), the encoderdetermined position (in steps and in inches), and an estimate of time remaining until completion. To begin roughing, a 1/4-inch ball-end endmill is placed in the grinder collet and tightened firmly. After homing the mill to zero using the limit switches on the three axes, the endmill is positioned over the starting corner of the surface. It is lowered carefully until a sheet of paper can no longer slide between the endmill and the surface. This position is recorded as the reference level. Next, the endmill is moved over a small block of metal attached to the table near the slab, and the process is repeated. The z coordinate is recorded and is used later to verify the endmill position. For example, the indexer and program have no way to detect when an endmill slips within the collet, but the metal block test will show the slippage. The block also provides a way to determine the starting position after changing tools. The first roughing pass is offset above the reference plane so that the maximum depth milled is 0.2 inches. Subsequent passes remove 0.2-inch layers; roughing is continued until only 0.1 inch of material to be removed remains. The electronic grinder is set to speed 2 for roughing. Roughing takes about 45 minutes per section per pass. An intermediate pass is next made using the finishing toolpath and a 1/16-inch ball-end endmill. The intermediate pass is offset 0.025 inches above the final surface. The grinder is set to speed 1 for the intermediate and finishing passes. The slower speed is necessary due to the reduced rate at which shavings are cleared from the flutes in the smaller endmills. (Shaving buildup causes the endmill to become hot.) The intermediate passes last about 5 hours per section.

99 radius = 0.0156 inch - L< 1/32 inch endmill spacing = 0.0137 inch Figure 5.6: The spacing between adjacent finishing rows (0.0137 inch) was slightly less than the radius of the 1/32-inch ball-end endmill (0.0156 inches). The edges are milled after the intermediate pass using the 1/16-inch endmill and the toolpath written for the edges. Two edge passes are completed: one with offset 0.1 inch, and another with an offset of 0.025 inches. The finishing passes are run using a 1/32-inch ball-end endmill. The miniature endmills were purchased from R & F Micro Tool Co. of Pembroke, Massachusetts. Finishing passes last 5 hours as well. In milling test wood surfaces and the first UHMW surface, it was determined that often the surface material was not completely removed in the finishing passes. This was partially due to the limited overlap between adjacent rows when using the 1/32-inch ball-end tools (as shown in Figure 5.6), and partially due to the endmills becoming less sharp with use. The residue remaining could sometimes be pulled off in strips, but often it formed individual fibers that could be described as plastic "peach fuzz." Sanding, singeing, and brushing were tried in attempts to remove the fibers, but the procedure found to produce the fewest fibers was to make the initial 1/32-inch finishing pass in the cross direction (i.e., orthogonal

100 to the section's long dimension) using zero offset, followed by a lengthwise finishing pass (i.e., along the section) with a positive offset of 0.002 inches. This procedure reduced the fibers to acceptable levels. Keeping the tools sharp was a constant task. The endmills tended to become dull quickly due either to the material or to the high rotation speed. Tool dullness was particularly troublesome with the 1/32-inch endmills, which have very little cutting area, and which need to be the most sharp. Some of these endmills seemed dull when new. After consulting with R & F Micro Tool, endmills coated with titanium nitride (TiN) were found to keep a sharp edge longer. The finishing process was quite lengthy, with 28 section passes at 5 hours each. The mill ran unattended, but was checked by an operator every few hours. Three sections per day was the effective maximum rate of completion. 5.6 Surface Completion and Temperature Sensitivity Programming and testing of the upgraded system continued through September 1993. Improvements during this period included the addition of a blower fan to keep the z-axis motor cool and optimization of mill travel speeds for roughing, intermediate, and finishing passes. In late September, work was started on the first artificial surface, UHMW1. A new problem appeared during the intermediate-scale milling on UHMW1. After noticing that the slab was slightly loose within the holding clamps early in the morning, a heat gun was used to heat the surface slightly. The slab soon became immobile, having expanded due to the increase in temperature. The thermal expansion and contraction became a significant problem in the milling process. When the temperature dropped, the slab became loose; when the temperature rose too high,

101 the slab bowed upward because the mounting brackets prevented lateral expansion. Either of these changes left visible ridges in the milled surface. The thermal contraction and expansion tended to lag behind the room temperature by several hours. A heat lamp and space heater were tried as means of keeping the surface warm, but neither could keep the slab at a constant temperature. In October, the steam supply to the room heater was turned on, but the radiator was unable to keep the room above 20 0C at night even at its maximum setting. (The milling machine was located in room 422-4 in the Aerospace Engineering Building; this room is quite large with a 20-foot ceiling and has windows on three sides which are poorly sealed.) The ambient temperature had been 24-27 0C when the surface was first mounted; maintaining this temperature became more difficult as outside temperatures fell through October. After repeated maintenance work on the room heater and improved sealing on the windows failed to help, a cubical enclosure ("the Greenhouse") was constructed from light wood and plastic sheeting around the milling machine., A compact electric heater with thermostat and a blower fan were placed inside in an attempt to stabilize the the temperature. The electric heater could not maintain an 80 'F air temperature around the clock, but it did provide sufficient heat during the day to allow UHMW1 to be completed on 19 November. A verification scan was performed on the surface before removing it from the mounting clamps. The scan is described in the next section. Milling on UHMW2 began on 11 March 1994 and was completed on 14 April. Temperature stability was less problematic because the milling was started at a lower temperature that was maintainable by the electric heater. The milling machine was moved to room 422-15 before starting UHMW3 on 14 June. Room 422-15 had an air conditioner, which was necessary to avoid slab

102 expansion due to overheating. UHMW3 was completed on 9 July. 5.7 Surface Verification After each surface was completed, the electronic grinder was removed and a small ranging laser mounted on the milling machine. The new laser, which was purchased to replace the surveying laser used in Chapter III, produced elevation measurements over a limited range at a very rapid rate. The laser generated an output voltage proportional to the measured distance; this voltage was read by a sampling multimeter and relayed to the 486 computer through an HP-IB link. The surface scanning program, millscan, was written in Visual Basic and was quite similar to the milling program. The laser beam was carefully adjusted to be vertical and was aligned as closely as possible to the origin of the milled surface. The surface scan consisted of 264 x 264 points with a point spacing equal to five times the point spacing used in the original surface specification. The elevation at each point was based on the average of 100 samples measured at two-millisecond intervals. An error estimate was calculated by subtracting the scanned elevations from the elevations specified in the original surface file. The error estimate from UHMW1 is shown in grayscale in Figure 5.7. No systematic features are visible in the file. The surface in the figure is actually a subset of the error estimate; several rows around the edge were removed because the laser range estimates were visibly wrong near the mounting brackets. The corrupted readings were most likely due to the angled orientation of the detector which senses the reflected laser light. This configuration seemed to decrease the laser accuracy on sloped surfaces, suggesting that the width of the error distribution (shown in Figure 5.8) was inaccurately large. The repeatability

103 Figure 5.7: Grayscale file of estimated error in UHMW1, calculated by subtracting the scanned height from the specified height in the original surface file. () 0 () L_ L_ L__.0 E 1000 -800 -600 -400 -200 - o -1 — - -0.3 -0.2 -0.1 0.0 surface error (inches) 0.1 0.2 Figure 5.8: A histogram of differences between the specified the scanned surface height. surface elevations and

104 of the endmill position on the metal block position to within 3-6 steps (0.03-0.06 mm) after numerous milling sessions supports the hypothesis that errors in the milled surface elevations are below the resolution of the scanning laser. Fabrication of the surface analogues was by far the most time-consuming part of this study. Measurements of radar backscatter from these surfaces is described in the next chapter.

CHAPTER VI RADAR MEASUREMENTS Backscatter from the artificial power-law surfaces was measured at 35 GHz using the Radiation Laboratory Millimeter-wave Polarimeter. Measured values of the scattering coefficients of the scaled surfaces were then interpreted as X-band scattering coefficients of equivalent larger-scale surfaces. (The scale factor turns out to be unity.) In this chapter, the radar components and measurement configuration are described. The calibration and measurement algorithms are discussed, and final post — processing and data reduction steps are detailed. 6.1 35 GHz Radar The Radiation Laboratory Millimeter-wave Polarimeter (MMP) was first assembled in 1986 [62] and was most recently rebuilt in 1994. The MMP is a dual-antenna, network-analyzer-based scatterometer designed for operation at 35 GHz. A phaselocked on-board local oscillator signal is mixed with the 2-4 GHz IF sweep received through cables from the network analyzer to produce a transmitted wave sweep from 34 to 36 GHz. The received signal is downconverted by mixing with the signal from a second local oscillator, resulting in a 2-4 GHz IF signal that is transmitted to the 105

106 network analyzer through coaxial cables. The transmit and receive local oscillators are phase-locked by split signals produced by an on-board X-band oscillator that are multiplied up to the LO frequency. Motorized rotating wave plates in the transmit channel allow the transmission of arbitrary polarizations, and separate horizontal (H) and vertical (V) receive channels allow the simultaneous measurement of H and V return signals when using a network analyzer with dual-channel reception capability. The transmit and receive assemblies may be separated for use as a bistatic radar. In addition, the removable transmitter module may be installed in the receiver assembly allowing monostatic operation using a single antenna. A block diagram of the radar is shown in Figure 6.1. For this study, the radar was configured as a dual-antenna, quasi-monostatic instrument operating at the near edge of the far-field region. The transmit and receive antennas were installed adjacent to each other and offset by 4.8~ so that the transmit and receive antenna boresights converge at the target range. The 3 dB product beamwidth was 30 and the 6 dB product beamwidth measured 4.20. While the original radar specifications called for components with a bandwidth of 2 GHz, a slightly wider bandwidth of 2.4 GHz (33.6-36 GHz) was used for the measurements in this study. An HP 8720 network analyzer was used to process the returned signals and translate them for output to a Gateway 2000 DOS-compatible computer. The HP 8720 was preferred over an HP 8753 because its frequency range allows input and output of signals up to 4 GHz, making IF sweeps of 2-4 GHz (or 1.6-4 GHz) possible. This advantage was offset by the test set of the HP 8720, which has only two ports, precluding reception of both V and H channels simultaneously. An external microwave switch was used to connect either the V or H receive channel to the network analyzer.

107 r- -* -* -.- -.- - -* -* - -* - - - - - -* - -_ - -* - -_ - -* - - -^ - -.- -_ -.- -* -* -.-.-.- -* - - - - - -* - -.- - -< -* -*, -f RECEIVER MODULE IF OUT V --------------- ( ) 2-4 GHz IF OUT H I 9 I NOUT PIN iIF UP 2-4 GHz w Polarization Switches I FOU --- —------ ----------- - ---------- -------- Figure 6.1: Block diagram of the 352-4 GHz radar. -E^O-^I! --- —Polarzatio Swtce TRNMTERMDL Figue 6.: Bock iagam o the35 ~z rdar

108 This switch was operated by the computer using an HP 3488A switch/control unit. The time-gating capability of the network analyzer allows the operator to isolate the response from a selected range interval in the time domain. This feature was used in the present study to mask out reflected signals from the bottom of the UHMW stack, the turntable, the wall behind the target, and other objects in the room, making accurate cross section measurements possible without an anechoic chamber. The network analyzer, transmit polarization switches, receive channel switch, and target turntable were all controlled by a Gateway 2000 4DX2-66V DOS-compatible computer using control software written in Visual Basic. Raw and time-gated trace data were stored on the computer's hard disk drive. The radar and measurement control system is shown in Figure 6.2. 6.2 Measurement Configuration The scattering measurements were performed in a side area of the Radiation Laboratory Bistatic Scattering Facility, Room 422-4, Aerospace Engineering Building. Figure 6.3 shows a measurement schematic. The milled surface was placed atop an 18-inch stack of UHMW slabs (two 6 -inch slabs on top of three 2-inch slabs) in order to increase the delay between the reflections from the top and bottom surfaces, allowing use of a time gate to isolate the upper surface reflection. The slabs were flat except for a small air gap around the edges of the two 6-inch slabs and a slightly larger air gap (maximum thickness: 1 mm) between the milled surface slab and the top of the UHMW stack. This larger air gap was due to a slight warping of the milled slab caused by the release of internal stresses during the milling process. A dielectric filler (Vaseline petroleum jelly, r, = 2.16 - j0.0022 at 10 GHz [29]) was used to fill the air gaps. Plots of

109 Polarization switch control Figure 6.2: Components of the radar measurement system.

110 Figure 6.3: Schematic of measurement setup. the time-domain response before and after applying the filler showed a reduction in reflections from the internal interfaces. The UHMW stack was placed on a square piece of Eccosorb AN-PXP-74 foam absorber. Surfaces UHMW1 and UHMW2 were measured with the slab stack and absorber placed on the large turntable used for the sandbox in the bistatic system. A much smaller turntable was used for the UHMW3 measurement; a square of 1/4-inch plexiglass was placed between the turntable and absorber to increase stack stability. The target turntable was placed about a meter in front of an absorber wall to eliminate reflections from the corner of the room. Other objects were removed so that the artificial surface was the sole object at target range near the main beam of the antenna. A large section of absorber was placed in possible sidelobe directions, but no significant reflections due to sidelobes were detectable. The radar transmitter and receiver assemblies rotated about a front axis that was mounted on the front of a large steel frame. The incidence angle was adjusted by pivoting the radar around the front axis, using a small electric motor to reel in a cable attached to the rear plate of the receiver. The steel frame was attached to a manually operated, mechanically telescoping lift (a Genie Superlift, manufactured by Genie Industries of Redmond, Washington). The lift was used to raise the radar to

111 heights up to 3.5 m, making measurements at near-normal incidence angles possible. The system had a high signal-to-noise ratio. The signal peak corresponding to a 1.75-inch sphere was about 58 dB above the noise level before calibration. Isolation between the co- and cross-polarized channels was about 25 dB before calibration and about 42 dB after calibration over most of the bandwidth. 6.3 Calibration The radar was calibrated using a recent technique for polarimetric coherent-onreceive radar systems [46]. The technique uses two calibration targets (a metallic sphere and a depolarizing target) for the initial system calibration and a sphere alone for subsequent gain corrections. A 1.75-inch (4.445 cm) diameter steel sphere and a pair of 0.028 cm diameter x 5.4 cm lengths of wire mounted in styrofoam were used as the calibration targets for these measurements. The calibration technique is described in Appendix B. The radar system stability needed to be established before measurements of scattering by the artificial surfaces could begin. Repeated calibrations using a sphere and depolarizing target were used to test system stability as reflected in the calibration coefficients. Several problems were fixed in order to get repeatable values for the time-invariant calibration constants cl, c2, C3, rI, and r2: * The use of trace subtraction to remove the response from the calibration target mount was discontinued because the mount position changed when placing or removing the calibration target. Fortunately, the scattering cross sections of the calibration targets were much higher than that of the styrofoam mount. * A sturdier styrofoam column was used. The steel sphere used as a calibration target was heavy enough to compress and give a slight tilt to the original

112 styrofoam mount. * A phase-lock problem was fixed by slightly changing the frequency of the Xband oscillator used to phase lock the two local oscillators. * The radar temperature was stabilized, first by installing sheets of clear plastic to make the radar assembly near-airtight, and later by installing a heater inside the enclosure. * The calibration target was allowed to settle for at least 15 minutes after placing it on the styrofoam target mount. Repeated tests indicated a non-visible sway introduced by placing the target (especially the heavy sphere) on even the sturdier styrofoam column. Target stability is particularly important while measuring C1, C2, c3, TI, and r2; the settling period increased the repeatability of these factors. * A minimum warm-up period of one hour was selected for both the radar and network analyzer. These measures, together with more mundane practices like avoiding contact with the radar cables or mount, brought the system to a sufficient level of stability, as reflected in the time-invariant constants and the channel gains R1 and R2. Test measurements of the calibration sphere (Figure 6.4) and a smaller sphere measuring 1.188 inches in diameter(Figure 6.5) showed that deviations from theoretical scattering cross sections were less than 0.5 dB (usually less than 0.3 dB) and that cross-polarized isolation was good, usually 40-45 dB or more. Although the system stability seemed good over fairly short periods, it was suspected that longer-term variations would be greater. A conservative approach of calibrating before and after

113.i IItI -27.6 - 1C",~ -,- Theoretical, E...... VV measured 7. 8 HH measured X28 - 2.0 -cn a) -,.;..'/ o, 3.. // -28.2 - -:...' ~,,...*,. co 33.6 34.2 34.8 35.4 36.0 frequency (GHz) Figure 6.4: Test measurement of calibration sphere (1.75 inches in diameter). -30.8 -' c -31.0 -E -31.2 -31.4 - 'retical 0..~ VV measured.... -32.0 Wmeasured ' --- HH measured '-32.2.. / 'I '"..," 33.6 34.2 34.8 35.4 36.0 frequency (GHz) Figure 6.5: Test measurement of 1.188-inch diameter sphere.

114 each rough surface data set was adopted. During the processing of the measured data from UHMW1 and UHMW2, it was noticed that the magnitudes of R1 and R2 computed after a measurement run were slightly different from those measured before the run. The relative phases also varied as a function of frequency. Linear interpolation of the magnitudes and phases of R1 and R2 with time was used as a first-order correction to minimize the effect of the changing R parameters. Values of R1 and R2 for a given trace are based on linear interpolations computed with respect to the actual times of the previous and subsequent calibrations and the start and stop times of the measurement set. Test measurements in which the phases of R1 and R2 were 180~ off showed that errors in o were only about 0.6 dB; the linear-interpolation-based values should be considerably more accurate in phase. 6.4 Measurement Procedure After characterizing the system by performing a full calibration to determine the time-invariant calibration coefficients, the UHMW stack was placed upon the turntable, and the surface to be measured was placed on the stack, as discussed in Section 6.2. Dielectric gel was used to fill in air gaps under the edges of the surface to be measured. Before beginning measurement, the calibration sphere (1.75 inches in diameter) was placed on its mount. The radar was moved to the operating range for this study, 2.7 m. (The range was measured using a 2.7-meter length of wire attached to the radar faceplate.) After carefully peaking the sphere signal by adjusting the radar azimuth and elevation angles, a red targeting laser (attached below the receiver antenna) was adjusted to point at the center of the antenna beam at the operating range. The laser made later pointing much simpler.

115 The artificial surfaces were measured at incidence angles ranging from 15~ to 60' at 5~ increments. The measurement sequence at each incidence angle consisted of the following steps: * with radar lowered, set incidence angle, 0 * set h, radar height, and d, distance from lift base to target surface * raise radar in elevation to point at sphere and calibrate * reset incidence angle, 0 * raise back to height h and position using target laser and range wires * measure 60 samples of surface backscatter using VV, HV, VH, and HH polarizations, rotating the target 5~ between measurements * raise radar in elevation to point at sphere and recalibrate A few additional sphere measurements were inserted between calibration sets. 6.5 Post-Processing and Data Reduction During the measurement runs, raw trace data is saved to disk files in binary form for reasons of speed (it is much faster than saving gated data) and flexibility (the data may be gated later in various ways). Gating refers to the process of transforming a frequency-domain data trace to the time domain; applying a time gate (i.e., a bandpass filter of specified width that is used to eliminate signals from objects at ranges outside the gate), and transforming the result back to the frequency domain. The network analyzer was set to its standard gate shape; the gate width was set to 5 ns in order to include reflections from both

116 near and far parts of the illuminated surface while excluding reflections from the bottom of the stack. The gate width was reduced to 4 ns for incidence angles of 25~, 20~, and 15~ due to the smaller range difference between the closest and most distant surface points. The gate center was set at the mean target range, usually 18.01 ns. After the surface and calibration target files have been gated, they are uploaded to a Unix workstation for further processing. Measurements of the calibration target are first converted into values of R1 and R2 as a function of frequency (i.e., for each of the 401 discrete frequencies that make up a data trace). Next, files of gated surface data are converted into average radar cross sections, a, for each polarization as if the targets were point targets (rather than distributed targets) using ot(fi) = (|2V/rSRt(f)I2)1 (6.1) where r and t indicate the receive and transmit polarizations and may be either V or H, and Srt is an element of the scattering matrix given by (B.43). The angle brackets (.) in (6.1) denote averaging over all traces in a set (i.e., over all surface measurements at various azimuth angles). Averaging is done separately for each of the 401 trace frequencies. 6.5.1 Scattering Coefficients Values of the scattering coefficient cra~(fi) are obtained by applying an illumination integral correction to a(f1). The correction may be derived from the forms of the radar equation used for point and distributed targets. The point target form for a target at boresight is [62] P (O) = PtGGA 1 (), (6.2) (4ir)3 R

117 and the distributed target form is P() PtGotGOrA2 (6.3) Pr()= (47)3 I(h, )a(0), (6.3) where I(h, 0) is called the illumination integral: I(h, ) = m t (t Ot)gr (Or, ) R4 dA. (6.4) um. area R4 (f>) The illumination integral accounts for the variation in antenna gain and range across the distributed target. The illumination integral correction is obtained by setting the ratio Pr(0)/Pt equal in (6.2) and (6.3) and solving for o~ in terms of a: (f) = R4(h a(i) (65) The program measillum.f was used to convert values of oa(fi) to values of o~(f,) for four polarizations at 401 frequencies per data set based upon the range and incidence angle for each measurement. Normalized antenna patterns were based upon data supplied by the manufacturer. Final values of a~ were obtained by averaging over frequency. Errors in the measured cross section of test spheres were greatest at the upper and lower ends of the frequency band; frequencies outside the range 33.8-35.8 GHz were therefore omitted from the frequency averaging. The values of the scattering coefficient were obtained using a surface model scaled for use with a 35 GHz radar instead of an X-band radar. Using Table 11-6 of Ruck [50], we see that the scaled quantities (denoted by primes) of length, conductivity, frequency, wavelength, permittivity, permeability, and radar cross section are related to the unscaled quantities as follows: ' = l/p, (6.6)

118 a' = pa, (conductivity) (6.7) f' = pf, (6.8) ' = A/p, (6.9) 6' =, (6.10) / = -, (6.11) a = oa/p2, (radar cross section) (6.12) where p is the scale factor (3.538 in this study). The conductivity scaling is not a problem because UHMW has negligible conductivity. The scale relation for the scattering coefficient of a distributed target, c~, is obtained by noting that the surface area, A, scales as length squared: A'=(') p2 (6.13) Writing the scattering coefficient as the average radar cross section normalized by area, we obtain )1 Al2 A C(6.14) (ao),- (CT') - (ti/p2) = (T) = 0 (6.14) A' Alp2 A Because a~ is a normalized quantity, it is independent of the model scale. 6.5.2 Independent Samples Independent samples of the surface scattering were obtained through a combination of target rotation, spatial sampling, and frequency averaging. The radar 6-dB footprint was positioned as close to the surface edge as possible to minimize overlap between footprints (Figure 6.6). Overlap of the 3-dB footprints was smaller still. The surface was rotated 5~ between measurements. A few measurement runs were executed using smaller angular increments, AO, in order to study the decorrelation as a function of target rotation. Based on limited angular decorrelation sets,

119 6 dB footprint Figure 6.6: Arrangement of footprints on slab decorrelation angles for UHMW3 are approximately 1.1~ for 60~ incidence and 4.5~ for 15~ incidence (based on 60-sample traces with AOq = 0.5~). For UHMW1, the decorrelation angle is approximately 1.5~ for 60~ incidence (based on a 90-sample trace with AX = 2~). A data set with fine angular increments was not taken on UHMW1 at 15~ incidence, but estimates based on the regular data sets indicate a higher decorrelation angle, though probably still less than 10~. No decorrelation sets were collected on UHMW2. The number of independent samples obtained through surface sampling Nss therefore varied from 300~/7.5~ = 40 (estimate for UHMW1 at 15~) to 60 (the nimber of samples collected) for UHMW1 at 60~ and UHMW3 at all incidence angles. Values of ao were averaged over a 2 GHz bandwidth, increasing the effective number of independent samples by a factor dependent on the footprint's extent in range. Ulaby et al. [57] give an expression for the effective number of independent samples obtained through continuous integration over a swept-frequency bandwidth

120 surface 0 Nss Nfavg Ntotal UHMW1 15~ 40 (est.) 1.3 52 60~ 60 5.1 306 UHMW3 15~ 60 1.3 78 60~ 60 5.1 306 Table 6.1: Estimated number of independent samples B: B fa [(1 sin d,a (6.15) where a = 27rD/c and D is the difference in range between the nearest and most distant points in the radar footprint. Using the dimensions and bandwidth from the present study, Nfavg can be shown to vary from 1.3 at 15~ incidence to 5.1 at 60~ incidence. The total number of independent samples, Ntotal, is calculated by multiplying the numbers gained by spatial sampling and frequency averaging: Ntotal = NsNf avg (6.16) The resulting numbers are listed in Table 6.1. Values of Ntotal for UHMW2 were probably slightly lower than those of UHMW3 due to decreased roughness. Values of the backscattering coefficient measured for the three surfaces are given in Chapter VII.

CHAPTER VII ROUGH SURFACE SCATTERING MODELS AND COMPARISONS TO EXPERIMENT The surface scattering measurements described in the previous chapter were performed in May and July of 1994. In this chapter, results of those measurements are shown and compared to predicted values of a~ resulting from several rough surface scattering models. 7.1 Experimental Results Measured values of h7 and ahh for the three artificial surfaces are shown in Figure 7.1. The open symbols represent VV polarization and the filled symbols represent HH. Backscatter dependence on the surface roughness behaves as expected, with the roughest surface (UHMW1) having the highest backscatter, followed by the intermediate surface (UHMW3) and the smoothest surface (UHMW2). The backscattering coefficient decreases with increasing incidence angle for all three surfaces, as expected, although there is a small upturn in o~ values for both polarizations at 5.5~ for the two smoother surfaces. There is little difference between the ao and Chh for a given surface —avo is usually slightly higher, but never more than a couple of decibels. The difference becomes small at small incidence angles, as it should. 121

122 -5- - 0 * Surface A (UHMW1) A A Surface B (UHMW3) i S O U Surface C (UHMW2) -1 0- W HH t ) -15- - a) i ) i * ^ * ~ o 03 0) 0 - A 0 A o oA 0 0 25 A A * -30-,_,___________________________ 20 30 40 50 60 incidence angle (degrees) Figure 7.1: Measured co-polarized backscattering coefficients, ao and oh, of the three artificial surfaces: Surface A (UHMW1), Df = 2.355; Surface B (UHMW3), Df = 2.2275; and Surface C (UHMW2), Df = 2.1. Measured values of the cross-polarized backscatter, ah and o, are shown in Figure 7.2. The measured ch and a0h also follow expected trends, decreasing with increasing incidence angle, and decreasing for smoother surfaces. The measured values are generally 20-23 dB lower than the co-polarized backscattering coefficients. Values of 0vh and ahv are roughly equal, with the exception of a couple of anomalous points for surface UHMW1 (the roughest surface) at low incidence angles. The backscattering coefficients seem somewhat noisier in the cross-polarized case. Because there was no test target that allowed easy and direct verification of measurement accuracy for cross-polarized returns, it is difficult to state probable measurement errors for o0h and ahv. 7.2 Comparison to Physical Optics Model The Kirchhoff or physical optics model is one of the most widely used surface scattering models [59]. The Kirchhoff model represents the scattered field in terms

123 0 -25 - o 0 t. E -30 CD. o -35 -C0 O * Surface A (UHMW1) A A Surface B (UHMW3) o * Surface C (UHMW2) HV VH C 0 -40 - 20 30 40 50 60 incidence angle (degrees) Figure 7.2: Measured cross-polarized backscattering coefficients, crh and chr, of the three artificial surfaces: Surface A (UHMW1), Df = 2.355; Surface B (UHMW3), Df = 2.2275; and Surface C (UHMW2), Df = 2.1. of the tangential field on the rough surface. It makes use of the tangent plane approximation to the surface field, in which the surface field at a point is approximated by the field which would be present if the rough surface were replaced by a planar surface tangential to the surface at that point [56]. The tangential plane approximation is valid if the radius of curvature at every point on the surface is large with respect to the radar wavelength. Use of the Kirchhoff model is usually limited to small incidence angles [14]. An additional simplification usually performed to make the physical optics solution more tractable is obtained by expanding the integrand of the scattered field integral into a series in terms of the surface slope and keeping only the lowest-order terms. This simplification is called the scalar approximation because it reduces to the scalar formulation of Beckmann and Spizzichino [5]; it requires that the radius of curvature be larger than the radar wavelength and that the rms surface slope be

124 small. Commonly cited constraints are [56]: kl > 6 (7.1) Rc > A (7.2) m < 0.25 (7.3) where k is the wavenumber of the incident field, 1 is the surface correlation length, Rc is the average radius of curvature of the surface, and m is the rms slope. The scattering coefficient under the scalar approximation is given by [59].0 = or0 0 +0 pq = qc + pqn + pqs (7.4) where the subscripts p and q indicate the scattered and incident polarizations and the subscripts c, n, and s indicate the coherent, non-coherent, and slope terms. For backscatter at non-normal incidence, we omit the coherent term. The non-coherent term is given for the backscatter case by ' I 11 E (4k2 cos2 a2)n I()(-2k sin 0, 0) (75) n=l n! 0n = 0 (7.6) where Rfo and R110 are the Fresnel reflection coefficients for horizontal and vertical polarization, respectively, and 4L '(q, q) = J J (R u v) (e dud (7.7) For the power-law case, Rz(u, v) was calculated numerically from the modified roughness spectrum given by (5.8) using an inverse fast Fourier transform (FFT). The result was raised to the power n and then re-transformed to obtain I(n).

125 The slope term c~qs was omitted for the power-law case because the non-analytical correlation function prevented the direct evaluation of the slope terms. The OC'qs is much smaller than c~ in the Gaussian case, so it is assumed that the term is of little consequence. The physical optics results are compared with the measured a~ in Figure 7.3. The model values underestimate the measured backscatter by about 5 dB in the HH case, although they do follow the trend of the measured values out to about 45~ or so. The VV predictions start about 5 dB low at 15~, but rapidly grow worse due to the Brewster angle effect which the model has but which is not found in the measured data. The first-order physical optics formulation does not predict cross-polarized backscatter. The roughness spectra used in the physical optics power-law solution had a lowfrequency cutoff of fl, = (10A)-1, where A is the radar wavelength. The low-frequency cutoff was imposed due to programming constraints in the calculation of the I(n) terms. Letting fic vary from (5A)-~ to (25A)-1 resulted in a change in r~ of less than 1 dB at 15~ and still less at higher incidence angles. 7.3 Comparison to Geometric Optics Model An alternative simplification to the Kirchhoff formulation involves the use of the stationary phase approximation, which means that scattering can occur only along directions for which there are specular points on the surface [59]. This highfrequency solution is valid when the average radius of curvature and the surface standard deviation are large compared to the radar wavelength. Commonly cited

126 0 IA 10 -10 -o 0 -15 --20 -0) C.-) c- 25 -a) u,) -35 - 0 A O Surface A (UHMW1) A Surface B (UHMW3) O Surface C (UHMW2) --- Phys. Optics, UHMW1 --- Phys. Optics, UHMW2 0 0 A 0 0 0 0 0 0 A A 0 0l 0 ~L HH I I I I I 20 30 40 50 incidence angle (degrees) 60 6 0 (a) -5 --10 -0 t -15-.a) O. -20 -OD 0 0 c -25 -a) 1-30 -Cn -35 --40 - 0 8 0 A O Surface A (UHMW1) A Surface B (UHMW3) 0 Surface C (UHMW2) -- Phys. Optics, UHMW1 -- Phys. Optics, UHMW2 0 A ~ 0 ( A 0 0A 0 A \ o I 0 0 A A El 0 A VV 20 I I I 30 40 50 incidence angle (degrees) 60 (b) Figure 7.3: Physical optics results for (a) HH and (b) VV backscatter compared to the measured values for the three artificial surfaces.

127 criteria for this solution are [56] kl > 6 (7.8) R, > A (7.9) ko > (7.10) cos Oi - cos 0s where a is the standard deviation of surface height, Oi and Os are the incident and scattered angles, and k, 1, RC, and A are defined as before. The geometrical optics solution is given by Ulaby et al. [59]: 0 IR(0)12 PP (7.11) pp 2 cos4 02lp1,(0) |) 0p = (7.12) where R is the Fresnel reflection coefficient and p"(O) is the second derivative of the normalized correlation function, p(r), P(R) r) (7.13) evaluated at zero. It is the p"(O) term that provides the main difficulty in applying (7.11) to the power-law case. Because the power-law correlation function has no analytical form, it must be evaluated numerically at discrete points. A discrete approximation was used to evaluate p"(O): _ 2 — P-O Po -ro 1 P2 - pl Pi - Po = (P2 - 2pi + po) (7.14) where the first three correlation values po, pi, and P2 are numerically calculated using the same bandlimited spectrum used for the physical optics case. For the three

128 -5-o O Surface A (UHMW1) 1 ~ - - A Surface B (UHMW3) ".", o Surface C (UHMW2) o '- 1 5.......- Geom. Optcs )- 2 - \ -200,30 20 30 40 50 60 artificial surfaces, / = 1.23 mm, and estimates of the rms slope, m = \ jp~(O)j, were 0.31, 0.13, and 0.20. The geometrical optics results are compared with the measured backscatter in Figure 7.4. The model performs quite poorly for the power-law surfaces. 7.4 Comparison to Small Perturbation Model The small perturbation model was first applied to electromagnetic scattering by Rice [49], who used the Rayleigh hypothesis to represent the scattered and transmitted fields near a rough surface in terms of waves traveling away from the surface only. These fields are expanded in a power series about some small parameter such as surface height or slope. The fields are then calculated to some low order, usually first order, by neglecting series terms of higher order [56]. Validity limits on the small perturbation model are small perturbation model are ka < 0.3, (7.15)

129 m < 0.3, (7.16) where oa, the standard deviation of surface elevation, is computed "from frequency components of the surface responsible for scattering at a given electromagnetic wavelength" [59]. The small perturbation model is usually used for larger incidence angles (0; > 20~) [14]. In a study of non-Gaussian surfaces by Chen et al. [15], the authors found that the small perturbation model has a wider range of validity when applied to a surface with a non-Gaussian (in their case, exponential) correlation function for large angles of incidence. The first-order small perturbation model solution for the backscatter case is given by 0a = 8k442 cos4 lapq12W(2k sin, 0), (7.17) where ahh = RI(0), sin2 - cr(1 + sin2 ) ~V - (= r~ ) [er cos 9 + ([O - sin2 0)1/2]2' aCvh = Chv 0, and W(k~, ky) is the normalized roughness spectrum, W(k, ky) = J- j J p(u, v) exp[-j(ku + kyv)]dudv (7.18) which is written in terms of the normalized surface correlation, p(u, v). The Fourier transform in (7.18) uses a different form than that used elsewhere in this dissertation, where the non-normalized correlation function, Rz(T', ry), is defined as r00 o00 Rz(rT, r,) = J Sz(f, fy)expP[j2r(fTrT + fyy)]dfdfy. (7.19) -00 -OO The normalized spectrum W(k,, ky) with k, and ky in rad/m is therefore W(k, k) = 2^ 2 (7.20) 2i) r _U2 2r' 27

130 and we see that W(2k sin 0, 0) = 2-Sz ( s, ~) X (7.21) and therefore pq = 8k4 cos4 ~ 2pq12 SZ (s ),0 (7.22) where the spectrum, Sz, is given by (5.8) or an equivalent expression. Note that the a2 term drops out of the final SPM equation. The SPM results are shown in Figure 7.5 along with the measured backscatter. The SPM values track the measured data trends fairly well but have a negative offset of about 5 dB. In spite of this offset, the SPM results are the most useful of any of the models examined thus far. Because the SPM values are calculated based on a single frequency of the roughness spectrum, the offset cannot be attributed to an incorrect bandwidth choice. Because an error in the dielectric measurement was possible, different values of ~r were tried. Changing Er from 2.05 to 2.3 raised 7~ by about 1 dB. Errors in the er measurement of more than 0.3 are considered extremely unlikely. Adding an imaginary component of c" = 2 raises a2 by about 5 dB, but it is difficult to believe that the UHMW or the dielectric filler has a lossy part that is so radically different from tabulated values. In the absence of other explanations, the offset in the SPM predictions is viewed as a limitation of the model in the power-law case. 7.5 Comparison to PWH model The last model examined is the Phased Wiener-Hermite (PWH) model for dielectric interfaces by Eftimiu [16], which is an extension of earlier work by Eftimiu [17, 18, 19] and Eftimiu and Pan [20]. The model, which is relatively recent and which had no experimental verification, is based on an expansion of the surface

131 -1 0 t -1 5 20 0 Ca) -35 -40 c-25 cn - 3 50 -40' O Surface A (UHMW1) A Surface B (UHMW3) 3 Surface C (UHMW2) --.-............ SPM 0 A 0 20 30 40 50 incidence angle (degrees) 60 (a) -10 o O 0. —1 -5.) a -20 - -25 co 0.,) (n -30 -35 O Surface A (UHMW1) A Surface B (UHMW3) o Surface C (UHMW2) --.. —.......... SPM 20 30 40 50 incidence angle (degrees) 60 (b) Figure 7.5: Small perturbation model results for (a) HH and (b) VV backscatter compared to the measured values.

132 current in a series of orthogonal Wiener-Hermite functionals. Earlier uses of WienerHermite functional expansions in the study of turbulence led to their application to the problem of rough surface scattering by Nakayama et al. [41, 42, 43] and Meecham and Lin [38], among others. The expansion is attractive because it represents the surface current in terms of functionals based on the rough surface random process, ((x, y), and because lower-order terms in the expansion contain reflections of all orders [38]. An overview of the PWH model is the subject of Appendix C. The PWH formulation is algebraically intensive in comparison with the other models in this chapter, perhaps because it relies on few simplifying assumptions. Aside from use of the Wiener-Hermite functional expansion (and the assumption that a first-order expansion is adequate), Eftimiu's formulation uses the extinction theorem and the tangentiality of surface currents. The formulation is shown to reduce to the physical optics and small perturbation solutions in the appropriate limits. The PWH expression for a~ is rather lengthy and involves a number of intermediate terms, and is therefore left in the Appendix. The PWH solutions for a~ of the three surface analogues are shown for HH and VV polarizations in Figure 7.6. The modeled cr, which were calculated using the bandlimited roughness spectrum truncated at fie = (10A)-1, display a moreor-less linear dependence on incidence angle and predict higher backscattering for the rougher surfaces. However, the theoretical curves do not match the measured values very well, either in magnitude, difference between the rougher and smoother surfaces, or in angular trends. Recalculations in which a wider roughness spectrum (down to (25A)-1) was used resulted in rough and smooth oa~ values that were closer together in magnitude, but the values overestimated the measured values worse than in the narrower spectrum case. Values of ao were closer, differing from the measured

133 values by 6 dB or less. The trends of ao, vs. 0 are closer than those of Cah, but still not as close as those predicted by the small perturbation method. While there are no obvious flaws in the PWH approach, one can imagine several possible reasons why the model performs poorly in the power-law case. The algebraic complexity of the model is quite high, increasing the probability of mathematical or programming errors. (The cited articles by Eftimiu contained numerous typographical errors; the author was supplied with a copy of Eftimiu's derivation notes to aid in reproducing the derivations.) The model requires the evaluation of a number of double integrals. These may be analytically reduced to one dimension in the Gaussian case but not in the power-law case, where the dynamic range (several orders of magnitude) may introduce problems in numerical precision. Finally, the power-law surface may simply be too complex to represent by a first-order WienerHermite expansion, due to the presence of surface structure over a wide range of spatial scales.

134 0 Surface A (UHMW1) A Surface B (UHMW3) o Surface C (UHMW2) - PWH, UHMW1 -m- PWH, UHMW2 o -10 O.I:, so c Co = -20 CD 0 0 'i — 3 0, com en CO -3 20 30 40 50 incidence angle (degrees) 60 (a) -10 o 0). -20 0 -25 0 A E 20 30 40 50 incidence angle (degrees) 60 (b) Figure 7.6: Phased Wiener-Hermite model results for (a) HH and (b) VV backscatter for the roughest and smoothest surfaces compared to the measured values.

CHAPTER VIII CONCLUSIONS, CONTRIBUTIONS, AND RECOMMENDATIONS This chapter summarizes the conclusions of the present work and describes some of the lessons learned while carrying out the investigations described in the previous chapters. The principal contributions of the dissertation are reviewed, and recommendations for future research are given. 8.1 Conclusions of This Dissertation The conclusions reached through the efforts described in the previous chapters are perhaps best related as answers to the questions originally posed in Chapter I: 1. Q: Does the power-law roughness model apply to targets of interest in remote sensing? A: Yes. In addition to the cited reports, the surface profile measurements at Mount St. Helens (Chapter III) and the resulting spectral estimates (Chapter IV) show that some (but not all) volcanic debris flow terrains there have a power-law roughness spectrum over some range of spatial frequencies. Lowand high-frequency spectral estimates from the primary debris flow surfaces seemed to conform to a single power law from about 0.05 to 25 or 50 m-1. It 135

136 is expected that the intermediate-frequency region (roughly 0.25 to 2 m'1) follows the same power law, but this cannot be verified without additional surface profiles with A\ in the 12.5 cm range. 2. Q: If the power-law model applies to volcanic terrain, what are typical values of the power-law roughness parameters? A: Equations for the one- and two-dimensional power-law fits to the spectral estimates derived from the volcanic debris flow profiles are given by (4.51) through (4.56). Values of d ranged from 2.31 to 2.51, corresponding to surface fractal dimensions of 2.245 to 2.345, which are close to the value of 2.2 cited by Mandelbrot [34] and Voss [60] as producing realistic synthetic profiles. Estimates of the spectral slope having values greater than 2.0 would not be measurable using a spectral estimator subject to spectral leakage. The leakage problem described in Chapter IV may be responsible for the fixed values of /3 = 2.0 reported by Bell [6, 7], Huang and Turcotte [30, 31], and Sayles and Thomas [52]. 3. Q: Given a power-law roughness model, how does the backscattering coefficient vary with incidence angle as a function of surface parameters? A: Measured values of ao and ahh are shown in Figure 7.1 as a function of incidence angle for three artificial dielectric surfaces having power-law roughness spectra with the spectral amplitude, a, held constant and the spectral slope, 7y, allowed to vary. Measured cross-polarized backscatter is shown in Figure 7.2. In both cases, Oa~ decreases with increasing incidence angle. The co-polarized difference, a - 0hh, is largest between 45~ and 50~. Measured a'r and ahh seem most sensitive to the spectral slope between 40~ and 50~; ao and a0h are

137 most sensitive between 45~ and 50~. 4. Q: How well do various rough surface scattering models predict scattering by a power-law surface? A: None of the models examined were outstanding. The small perturbation model was most accurate in following the trends of a~ vs. incidence angle and in reflecting the difference in ao between the rougher and smoother surfaces, but the model underestimates the backscatter by roughly 4 dB for all angles. 5. Q: In particular, is Eftimiu's Phased Wiener-Hermite model [16] useful for power-law surfaces? A: The PWH model did not seem to perform well for the power-law surfaces, although cr, was much closer than arh. While the model was not grossly wrong, it did not model the trends vs. angle nor the rough/smooth backscatter difference as well as the small perturbation model or the physical optics model (oh7 case only). Possible reasons for the poor performance of the PWH model are discussed in Section 7.5. 8.2 Comments on Methods and Theory One of the principa lessons of the research performed for this project is a much deeper appreciation of the data requirements of statistical and spectral estimation. In the multi-faceted aspects of a remote-sensing study, it is easy to underestimate the amount of data needed for the estimation of some target parameter. The author has certainly learned that if spectral estimation, for example, is difficult, spectral estimation based on limited data sets is extremely difficult. Attempts to infer statistical information based on limited knowledge or data may do more than just spread

138 out the error bars a bit-the answers may be quite wrong. Many remote sensing studies would benefit from the advice of signal processing engineers. The use of manufactured surface analogues in scattering studies has great benefits but high costs as well. The time and data requirements for surface construction are quite high, especially for complex surfaces like the power-law surfaces. The development of new fabrication techniques, such as a recent technique in which a desired object is fabricated by solidifying a layer at a time using a laser scan, may reduce the costs in money and time to more manageable levels. Increased stability and automation would allow the radar measurements to have greater productivity. The use of surface analogues contributes to these goals, which are much more difficult to achieve in a field setting. While the accuracy of the Phased Wiener-Hermite model was disappointing for the power-law surfaces measured here, the concept of a Wiener-Hermite expansion of currents in terms of the surface random process still has appeal. Perhaps an application to a simpler surface might reveal how to adapt the model to the powerlaw case with more success. 8.3 Contributions of This Dissertation The research described in this thesis comprises five principal contributions. First, estimates of the roughness spectra of volcanic debris flows were measured and found to have, in some cases, a power-law roughness spectrum. The parameters kl and ka were shown to be poor descriptors of power-law terrains. Second, special problems in the spectral estimation of power-law surfaces were explored and documented. Spectral leakage was shown to be particularly problematic, giving false estimates of the spectral slope. Third, a computer-driven, mechanical system was assembled for

139 the manufacture of artificial surface analogues, and the system was used to produce three artificial surfaces. Fourth, radar measurements of two-dimensional, dielectric, power-law surfaces were performed. Fifth, a comparison of measured backscatter with that predicted by several rough surface scattering models was made, including an initial comparison with the Phased Wiener-Hermite model. 8.4 Recommendations for Future Work Because of the statistical nature of most of the measurements in this study, an obvious recommendation would be the repetition and expansion of the various measurements, i.e., collection of more volcanic (and other) surface profiles, calculation of improved spectral estimates of roughness spectra, and construction of more surfaces for indoor scattering measurements. Each of these tasks can become a major enterprise. A more targeted list of recommended research might include a few different items. For example, construction of larger (in area) and thicker surface analogues, if it could be made practical, would allow radar measurements with more independent samples, fewer edge concerns, and fewer problems with matching separate slabs to reduce reflections. If power-law surface construction could be made simpler or faster, additional surfaces could be used to determine the effect of the roughness amplitude (c in (2.8) or a in (2.9)) on ros, and whether the performance of scattering models improves for certain c and 6. Knowledge of the c dependence would aid in the selection of angles and polarizations for radar measurements designed to map roughness features, such as the SIR-C study of volcanic terrain associated with the work described here. Knowledge of the c dependence, and information on whether the local (high-frequency) c and d parameters are physically dependent on slope (due to

140 low-frequency roughness), would also contribute to another current research project, the modeling and interpretation of the diffuse radar backscatter from Mars. Another research goal should be to try additional scattering models on the powerlaw surfaces, such as the Integral Equation Model by Li and Fung [24, 33], or one of the two-scale models used in other scattering studies. The question of measurement scale and footprint size should be examined in detail. The different ranges of surface scales observed by orbiting, airborne, and truck-mounted radars may produce different estimates of co. If roughness of many scales, both smaller and larger than the radar footprint, is present, then the scattering mechanisms may be completely incoherent (since the wavelength is very small compared to the surface relief), and differences in observed a~ may be entirely due to differences in the sample size. If this were the case, averaging over more samples should make the truck-based oa~ approach the estimate derived from space-borne instruments. Finally, the Phased Wiener-Hermite model deserves further analysis, but experimental comparison might be better started with a simpler Gaussian surface, perhaps one with roughness in the intermediate range between the regions of validity of the classical scattering models.

APPENDICES 141

142 APPENDIX A MOUNT ST. HELENS DATA The filtering and interpolation steps performed on each data set are listed in Tables A.1 through A.7. The tables should be read from left to right. For example, the original file for scan 2b (described in Table A.1) was lp2bdata.all. A lower quartile difference (LQD) was then compiled for each 5 x 5 subregion and stored in file lp2bdata.all.qdf 25. (A median residual statistic is indicated by MR.) Pixels with lower quartile differences lower than -0.01 m were then discarded, resulting in file lp2bdata.fql. One iteration of median interpolation was then performed, resulting in the final profile, file lp2bdata.fqlel. file subregion size statistic residual file limit (m) lp2bdata.all 5x 5 LQD lp2bdata. all. qdf25 -0.01 lp2bdata.fql 1 iteration of median interpolation lp2bdata. fqlel final surface profile Table A.1: Filtering and interpolation steps for scan 2b (Johnston Ridge site).

143 file subregion size statistic residual file limit (m) lp2cdata.all 5 x 5 LQD lp2cdata.all.qdf25 -0.02 lp2cdata.fql 1 iteration of median interpolation lp2cdata.fqlel final surface profile Table A.2: Filtering and interpolation steps for scan 2c (Johnston Ridge site). file subregion size statistic residual file limit (m) lp3adta.all Remove last 5 rows and columns lp3axdata.all 9 x 9 MR lp3axdata.res81 -0.3 lp3axdata.fl 9 x9 MR lp3axdata.fl.res8l -0.2 lp3axdata.f2 9 x9 MR lp3axdata. f2.res81 -0.14 lp3axdata.f3 5 x 5 LQD lp3axdata. f3. qdf25 -0.015 lp3axdata.fq4 6 iterations of median interpolation lp3axdata.fq4e6 final surface profile Table A.3: Filtering and interpolation steps for scan 3a (Johnston Ridge site). file subregion size statistic residual file limit (m) lp3bdata.all 41 x 41 MR lp3bdata.restt -0.4 lp3bdata.fl 9 x 9 MR lp3bdata.fq.res81 — 0.25 lp3bdata.f2 7 x 7 MR lp3bdata.f2.res49 -0.15 lp3bdata.f3 5 x5 LQD lp3bdata. f3. qdf25 -0.01 lp3bdata.fq4 11 iterations of median interpolation lp3bdata. fq4el I final surface profile Table A.4: Filtering and interpolation steps for scan 3b (Johnston Ridge site).

144 file subregion size statistic residual file limit (m) lp3cdata.all Remove last 10 rows and columns. lp3cxdata.all 41 x 41 MR lp3cxdata.restt -0.4 lp3cxdata.f 9 x 9 MR lp3cxdata.fl.res8l -0.25 lp3cxdata.f2 7 x 7 MR lp3cxdata.f2.res49 -0.15 lp3cxdata.f3 5 x 5 MR lp3cxdata.f3.res25 -0.1 lp3cxdata.f4 5 x5 MR lp3cxdata.f4.res25 -0.06 lp3cxdata.f5 5 x 5 LQD lp3cxdata.f5.qdf 25 -0.01 lp3cxdata.fq6 14 iterations of median interpolation lp3cxdata.fq6el4 final surface profile Table A.5: Filtering and interpolation steps for scan 3c (Johnston Ridge site). file subregion size statistic residual file limit (m) lp4adata.all Remove last 10 rows and columns. lp4axdata.all 55 x 55 MR lp4axdata.restt -0.4 lp4axdata.fl 9 x 9 MR lp4axdata.f1.res8! -0.25 lp4axdata.f2 7 x 7 MR lp4axdata.f2.res49 -0.15 lp4axdata.f3 5 x 5 MR lp4axdata.f3.res25 -0.1 lp4axdata. f4 7 x 7 LQD lp4axdata. f4. qdf49 -0.01 lp4axdata.fq5 5 iterations of median interpolation lp4axdata.fq5e5 final surface profile Table A.6: Filtering and interpolation steps for scan 4a (Johnston Ridge site).

145 file subregion size statistic residual file limit (m) lp4cdata.all 51 x 51 MR lp4cdata.restt — 0.4 lp4cdata.f 9 x 9 MR lp4cdata.fl.res81 — 0.25 lp4cdata.f2 7 x 7 MR lp4cdata.f2.res49 — 0.15 lp4cdata.f3 5 x 5 MR lp4cdata.f3.res25 — 0.1 lp4cdata.f4 7 x 7 LQD lp4cdata.f4.qdf49 -0.01 lp4cdata.fq5 7 iterations of median interpolation lp4cdata.fq5e7 final surface profile Table A.7: Filtering and interpolation steps for scan 4c (Johnston Ridge site).

146 Grayscale plots of all profilometer and survey data are shown in Figures A.1 through A.6. Filtering and median interpolation have been performed on the profilometer data sets where necessary.

147 (a) (b) Figure A.1: (a) Scan 2b, Johnston Ridge site, 51 x 51 points, A = 2 mm. (b) Scan 2c, Johnston Ridge site, 41 x 41 points, A = 1 cm.

148 (a) (b) (c) Figure A.2: (a) Scan 3a, Johnston Ridge site, 46 x 46 points, A = 2 cm. (b) Scan 3b, Johnston Ridge site, 41 x 41 points, A = 5 mm. (c) Scan 3c, Johnston Ridge site, 41 x 41 points, A = 2 mm.

149 (a) (b) Figure A.3: (a) Scan 4a, Johnston Ridge site, 55 x 55 points, A = 5 mm. (b) Scan 4c, Johnston Ridge site, 51 x 51 points, A = 2 mm.

150 (a) (b) Figure A.4: (a) Scan 5a, Elk Rock site, 71 x 71 points, A = 5 mm. (b) Scan 5b, Elk Rock site, 51 x 51 points, A = 2 mm.

151 (a) (b) Figure A.5: (a) Scan 6a, Elk Rock site, 51 x 51 points, A = 2 mm. (b) Scan 6b, Elk Rock site, 41 x 41 points, A = 1 cm.

152 (a) (b) (c) (d) Figure A.6: (a) Survey 1, Johnston Ridge site, 33 x 33 points, A = 1 m. (b) Survey 2, Johnston Ridge site, 9 x 9 points, A = 4 m. (c) Survey 3, Johnston Ridge site, 9 x 9 points, A = 4 m. (d) Survey 4, Elk Rock site, 17 x 17 points, A = 2 m.

153 APPENDIX B COHERENT-ON-RECEIVE RADAR CALIBRATION The radar was calibrated using a new calibration technique designed for polarimetric coherent-on-receive scatterometers. The new procedure is convenient because it uses a sphere and an unspecified depolarizing target to determine timre-invariant distortion parameters and afterward uses only a sphere for gain calibrations. Knowledge of the scattering matrix of the depolarizing target is not required. The calibration technique is fully described in Nashashibi et al. [46]; a short derivation is given here because the notation and polarizations used differ from the cited article. A basic schematic of the radar system was given in Figure 6.1. The received signal due to scattering of the transmitted field from a target may be written: Er e-i2kor R1 0 1 cl 1 C3 =,,; E [ = T9 2 S TTT2E (B. 1) Eh 0 P, C2 1 C3 1 where Tgi = transmitter gain/loss, /l, R -= receiver losses in V and H channels, c1, c2 = receiver cross-talk (distortion) factors,

154 optical axis h Figure B.I: Coordinate system of transmitter polarizers S = target scattering matrix, C3 = transmitter distortion factor, T, T2 = polarizer transmission matrices, 1 Ei = initial field input to polarizers. 0 Let the transmission coefficients along and orthogonal to the optical axes of the transmitter polarizers be denoted by -y, and 7-, respectively. Using the coordinate system shown in Figure B.1, we can write the following relationships among the unit vectors: o = vcosa +hsina, (B.2) n = -vsina + hcosa, (B.3) v = 6 cos ca-n sin a, (B.4) h =6 sin a + ni cos a. (B.5) We can therefore write the polarizer transmission matrix as a function of the

155 rotation angle, a: = cos a -sina 7 0 cosa sin a T - sin a cosa 0 -sina cos a cos2a +r sin2 a sin a cos a(1-r)) = 'o sin acos a(1-r) sin2 a+rcos2a = -yT, (B.6) where r = 7n/7o, and the above equation holds separately for each polarizer, i.e. T1 = 7olTl(al, r), T2 = Yo2T2(c2,T2). Now define tl = = t Ti(Ci, r1)T2(a2, r2)Ei, (B.7) t2 and 1 C3 Et= t. (B.8) C3 1 We can prepare a table of t and Et for the principal six polarizations (Table B.1). Let R1 = R\Tglyol and R2 = R'Tgyo2. We now have -i2kor R 0 1 CI C3 0 R2 [C21 C3 [ t2 (B.9) There are seven unknowns in (B.9): R1, R2, C1, C2, C3, and r1 and r2, which are contained in t. Of these, only R1 and R2 vary with time; the others are fixed for a given system. Notice also that only t varies with transmit polarization. Let R1I 0 1 cl R1 Rlc R- -, (B.10) 0 R2 C2 1 R2c2 R2

156 Table B.1: Values of t and Et for various transmit polarizations.

157 and 1 C3 T=. (B.11) C31 We can write the received field vector due to transmitted polarization p as e-i2kor.. tl Er = 2RST. (B.12) t2 Let the first calibration target be a sphere, making 1 0 Ssphere = -SO (B.13) 0 1 and e-i2kor = tl E =,RsoT (B.14) r2 t2 Examining the cases of transmitted V, LHC, RHC, 450 linear, and 1350 linear with received V, we can write e-i2koro Evv -- oRi(1 + c1C3) (B.15) and similar expressions apply for the other polarizations. Now define A1 = (1 + T+ ) (1-)] (B.16) Bv E _[(1 +i (C+C3) (1-T1)]' (B.17) By inspection, 71 = A + BV-1. (B.18) Similarly, define ( + C3 ) (1 (B.19) T 1+ -- A 722 B)1. (B.19) Evv 2C L1 + C3C1/ J F = E =5 2 [(L T2 ) - ( 3 )(1 -2) (B.20) E 2vv 2 + 1 + ClC3J

158 and obtain rT2= D. + Fv-1. (B.21) The remaining five unknowns will be obtained using vertical and horizontal polarizations only. Write (B.14) as a matrix composed of the transmit vertical and horizontal cases: EVV Evh -i2koro = = tly tlh: = 2o.Ro (B.22) Ehv Ehh R t2v t2h j Post-multiply both sides by so1 times the inverse of the t matrix and define the left side as Ao: — 1 1 Ev Evh tpv tlh Ao Ehv Ehh t2v t2h ~-i2koro - = = RT e-i2koro R (1 + c1c3) Rl (ci + c3) 2 ' (B.23) 0 L R2(c2 + 3) R2(1 + C2c3) and noting that Ao is composed of known quantities, we can write four equations using the components of Ao: ~-i2koro all = 2 R-(1 l+ClC3), (B.24) e-i2koro a12 = 2 R-l(cl + 3), (B.25) e-i2koro a21 = -— 2 R2(2 + C3), (B.26) e-i2koro a22 = — 2 R2(1 +C2C3). (B.27) We still need one more independent equation. We obtain it by writing similar expressions for the depolarizing target and enforcing reciprocity. Note that we need not know the scattering matrix of the depolarizing target.

159 Writing an expression similar to (B.23): E= E Eh [ tlV tlH e-i2korl B== 2 RSDTT (B.28, Ehv Ehh t2V t2H 1 where the scattering matrix of the depolarizing target SDT has replaced. that of the sphere. Solving (B.28) for SDT: S = r i2korl -1 SDT = r2e 2 lR BT (B.29) where -1 1 R2 -R1C1 R R1R2 (1 - (B30)) -R2c2 R1 and = -1 1 -C3 T = (B.31) -C3 1 The fifth equation is obtained by setting expressions for Shv and Svh equal: - R2c2(bil - c3b12) + Rl(b21 - C3b22) = R2(bl2 - c3bll) - R1Cl(b22 - c3b2l). (B.32) This equation is independent of rl, the range to the depolarizing target. Solving (B.25) and (B.26) for R1 and R2, substituting into (B.24) and (B.27) and solving for c1 and c2, then substituting all four into (B.32), we obtain a quartic equation for C3: - 2-c3 + 2- 3- 1 =0, (B.33) where o = a22b12 - a12b22 + a21bll - b21all i1 = a22b11 + a21b12 - al2b21 - allb22.

160 Two of the roots are ~1; these are discarded as unphysical. The remaining two roots are given by c2 - 21 + 1 = 0. (B.34) tPo We select the root with Ic31 < 1. The remaining unknowns follow directly: 3all - a2 (B.35) c = --, (B.35) 3a12 - all C3a22 - a21 C2 = (B.36) c3a21 - a22 R = alz2 ri2koro (B.37) C1 e (B37) R2= a2 ri2koro. (B.38) C2 + C3 Values of cl, C2, C3, T1, and T2 for each frequency are stored in a disk file; these quantities are constant unless the system components are altered or replaced. Subsequent calibrations need only involve a metal sphere, measured with vertical and horizontal polarizations only. From (B.15), the equation for the VV sphere response is e-i2kors Evv= -- soRl1(1 + ciC3), (B.39) and the HH response is -i2kor, Ehh = 2- o R2[c2( + 71T2) + C2C3(1 - T12) + c3(1 + T1i2) + (1 - 'T2)] r2 e-i2kor SO = -R2D (B.40) rs2 2 Solving for R1 and R2: R _Evv r 2 2kor. (B.41) I + ClC3 a Ehh 2 2 i2kor, R2 -r -0o8 (B.42) After calibration, the target scattering matrix may be obtained directly from

161 measurements of the scattered field: -1 i2kort = - Evv Evh tIv tih = -1 Starget = rt i2k~rt v l t - (B.43) Ehv Ehh t2v t2h = -where 1 and are given by (B.30) and (B.31). where R and T are given by (B.30) and (B.31).

162 APPENDIX C THE PHASED WIENER-HERMITE MODEL The Wiener-Hermite (WH) functional expansion [12, 63] has been applied to a number of problems in turbulence and wave propagation in random media. Applications of the WH functional expansion in studies of rough surface scattering include works by Nakayama et al. [41, 42, 43] and Meecham and Lin [38, 39]. The paucity of rough surface scattering models that were applicable to dielectric interfaces and that had potential for use on very rough surfaces led to an examination of the Phased Wiener-Hermite (PWH) model described by Eftimiu [16], which is an extension of earlier work by Eftimiu [17, 18, 19] and Eftimiu and Pan [20]. The use of functional or series expansions is widespread in science and engineering. Expansion of a complicated function in terms of a series of simpler functions with convenient properties (e.g., odd or even symmetry, smoothness) permits operations to be carried out term-by-term, and with luck (or skill), the number of terms needed will be small. Consider the Fourier series expansion of a periodic function f(t) [53]: f( ) = E Cne (C.1) n=-oo The coefficients c, are constants and are determined by multiplying each side of (C.1) by the corresponding e-inwt term and integrating over a period. Because the einwt

163 terms are orthogonal, all the integrals on the right side equal zero save one, and the coefficient is given by c,=l f(t)e-inwdt, (C.2) where r is the fundamental period of f(t). Observe that if the function f(t) has a form that is similar to or mathematically related to the form of the expansion functions (einwt in the above example), the series expansion will have fewer terms. The function sinwt is perhaps a trivial example; its Fourier series has only two non-zero coefficients, cl and c-1. The Wiener-Hermite functional expansion is another form of series expansionone that is designed to work with functions of a random process. The WH functional expansion is a series representation of a function of a random process written in terms of Hermite polynomial combinations of some random function. For example, the first four Wiener-Hermite functionals based on a one-dimensional stationary random process ((x) are [38]: H()() = 1, (C.3) H()()= (x), (C.4) H(2)() (x)((x2)- R(xi- X2), (C.5) H(3)(C) = C(Xl)((2)((X3)- R(x - x2)((X3) -R(x2 - 3)(x1) - R(x3 - Xl)((x2), (C.6) where R is the correlation function of the random process. The WH functionals are constructed to be statistically orthogonal; i.e., they are orthogonal under an ensemble average or expectation operation: (H(m)(()H(")(()) = 0 for m $ n. (C.7)

164 The WH functional expansion of f(x) = eikC(x) is a simple example. The expansion will use WH functionals that are a special case of those given by (C.3) to (C.6): H()() = 1, (C.8) H()() = x), (C.9) H(2)() = (2()- 2, (C.10) H(3)(C) = 3(Xl)- 32((x)), (C.i1) where ao2 is the variance of the random process. The WH expansion has the form f(x) = eikC(x) = aoH(o)(C) + alH()(() + a2H(2)(C)+** = ao + al((x)+ a2[C2(x) - o2] +... (C.12) The coefficients an are non-random constants and are found using a procedure like that used to find the coefficients Cn of the Fourier series. The coefficient an is found by multiplying both sides of (C.12) by the WH functional H(n)(() and ensemble averaging both sides. For example, the n = 1 term is found as follows: (1. eikC()) = (1. ao) + (1 * al((x)) + (1. a2[((x)- a2]) +.., (C.13) and all terms except one are zero on the right side due to orthogonality, giving a = (eikC(x)) = e kCfz((C)d J -00 ao = e-202/2, (C.14) where we have used the normal Gaussian probability distribution function to calculate the ensemble average. The al and a2 terms may be calculated similarly, giving

165 the WH expansion fwH(x) = e-2/2{1 + ik(x) + (ik)2 [2() - 2] +.}. (C.15) Note that the WH expansion is still a function of random quantities; however, the random process ( is known (i.e., its moments are known) and the ( are now in the form of a polynomial rather than in the argument of an exponential. While the distribution function of the C function is usually specified to be Gaussian, the correlation of the ( random process is left unspecified. The chief attraction of the WH expansion is its ability to closely approximate a function f(x) of a random process using a small number of terms. Meecham and Lin [38] compared the mean and variance of eikC(() to the mean and variance calculated using the second-order WH expansion (first three terms of (C.15)) and the first three terms of a perturbation (power series) expansion fPT of eikC(x): fpT(z) = 1 + ik((x) + [k((x)]2 (C.16) 2 Values of the moments were compared for several values of ka (Table C.1). Note that the expected value of f(x) is exactly equal to the first term of the WH expansion, therefore (f) - (fwH), regardless of ka. The perturbation expansion mean is badly off for ka = 2. The variance is also significantly closer for fwn than for fPT. The most general form of the WH functional expansion is [40] f(x) = H + K(1)(x - x)H(1)(xl)dxl + J K(2)(x -x, x -3)H(2(x2, 3)dx2dx3 + JJ K(3)(x -X4, x - x5, X6)H(3)(x4 X, x6)dX4dx:5dx6 +., I (C.17)

166 Exact WH PT kc I (f) a2 fwH) 2 (fPT) a2 0.5 0.882 -0.172 0.882 -0.170 0.87 -0.250 0.75 0.755 -0.245 0.755 -0.230 0.719 -0.563 1.0 0.606 -0.233 0.606 -0.183 0.500 -1.00 2.0 0.135 -0.018 0.135 0.073 -1.00 -4.00 Table C.1: Comparison of moments of WH and perturbation expansions of eikC(x) with moments of exact expression (data from Meecham and Lin [38]). where the "coefficients" K(n) are non-random integral kernels. (The kernels leading to (C.12) were products of a, and delta functions.) This general form, which can be shown to be complete [63], shows the dependence of f(x) on surface quantities at other points (X1, X2, X3, etc.). The original form of the WH expansion utilized Hermite polynomial combinations of the "ideal random function"-the Gaussian white noise function [63]. Later works introduced basis functions that vary with time [37] and the use of the surface random process itself as the expansion basis [38], as in the example given earlier. In his formulation of surface scattering by a rough dielectric interface [16], Eftimiu expands the surface electric current, J, the surface magnetic current, M, and an exponential function of the surface height into WH expansions in terms of the surface random process, ((x, y): J(x,y) = exp(ik[ysin0 - C(x,y)cosO ])[Fo + 1 F(x - x', y - y')((x', y')dx'dy'] (C.18) M(x, y) = exp(ik[y sin - ((x, y) cos 6])[Go

167 + G(x - x', y -y')C(', y')dx'dy'] (C.19) exp[-iu((',y')] =e- 2/2{1 - iU((x',y')- Z2 (-2(x') )- ]}. (C.20) We see that the forms of (C.18) and (C.19) are based on the general form (C.17) written to first order, modified by a multiplicative phase term. The phase factor is the origin of the name "Phased Wiener-Hermite model." This exponential term represents the value of the phase of the incident field on the surface. The exponential in (C.20) is written using terms up to the second order, but second- and higher-order products of the exponential and current terms are later reduced to zero or first-order expressions. C.1 General Outline of PWH Formulation The random surface z(x, y) is considered to be a single realization of the twodimensional random process C(x, y). The random process is assumed to be at least wide-sense stationary, having a zero-mean, Gaussian distribution: (((x,y)) = 0 (C.21) ((2(x,y)) = a2, (C.22) and having a normalized correlation function c(ra, rT) defined by C('r, ry) = 2C((, yY)C(X + r, y + ry)), (C.23) whose form is not specified. The unspecified correlation allows the model to be applied to both Gaussian and power-law surfaces, as will be shown later. The surface geometry and principal vectors are shown in Figure C.1. The upper medium is assumed to be free space, while the lower medium has a relative permittivity ~r and free-space permeability gHo. The incident field is defined in the yz-plane

168 Os z= C(x,y) EOgO y x Figure C.1: Geometry of PWH model

169 as E'(r) = e' exp[ik(y sin - z cos 0)], (C.24) where the incidence unit vector is for horizontal incidence, e = y sin 0 + z cos 0 for vertical incidence, and a time dependence of e-iwt is suppressed. The scattered field in the upper medium is given by the electric field integral equation: E'(r) = V x [ioV x J/ 4(r, r')j(r') dS' - ff (r, r')m(r') dS']. (C.25) Expressing the free space Green function in an eigenfunction expansion:, exp(ikr - r) 1 ff exp[i ~ (r- r')]d3 (C26) 47rlr - r'l (2r)3 K111 k K - ' ) where Kr = Kcx + Ksyy + Kz', and normalizing the current by the differential surface area gives another form of E3: E 1 // d3K ieidr E(r (2r)3y I K2- k dx' dy' x e-iKr' [ x J(', y') + M(x', y)]. (C.27) Similarly, the transmitted field below the interface is 1 [fff a r ' Er) (27r)3 ()2 (k)2 eiK'r J dx' dy'' x etK' [7Ki x J(x', y') + M(x', y')] (C.28) where c' = KxX + KcyY + /cz'z. We now have two integral equations in which the surface currents are unknowns to be determined. Once these currents are known, the scattered field above or below

170 the surface may be obtained directly from (C.27) and (C.28). Because the surface, and therefore the surface current, is a random function, the scattered fields will also be random functions. We therefore seek a WH functional representation of the surface currents so that moments of the scattered fields may be written as moments of the surface currents, which in turn (using the WH expansion) may be written in terms of moments of the surface random process, which are known. The surface electric current J is expanded into a WH expansion of first order: J(x', y') = exp[iky'sin 0 - ik((x', y') cos 0] {Fo + F(' - x", y' - y")((x", y") dx" dy"}, (C.29) where Fo and F are nonrandom vector constants or functions which must be determined. The inclusion of the ((x', y') term in the exponential differentiates the Phased Wiener-Hermite formulation from Eftimiu's previous models [17, 18, 20]. The magnetic current M is handled similarly: M(x', y') = exp[iky' sin 0 - ik((x', y') cos 0] {Go + G(x - x", y' y")((x", y") dx" dy"} (C.30) The exponential exp will also be e, xpanded into a WH series. Products of this series and the current expansion will have terms of order H(2) that must be removed in order to retain a model consistent to the first order. This process will be shown for the electric current, which we now rewrite using a transformed integral: J(x', y') = exp[iky'sin0 - ik((x', y') cos 0] {Fo + (21 J| F(a (, )ei(x'+y') da d/} = exp[-ik cos O(x', y')]{J(O) + J(1)} (C.31)

171 where J(o) = exp[iky'sinO]Fo (C.32) J(1) = exp[iky'sin 0] (2)2 F(a, /)((a, )ei(ax'+/y') dc d3 (C.33) and the tildes represent the Fourier transforms of the indicated variables. We now write the product J(x', y') exp[-iK,((x', y')] = exp[-iu((x', y')]{J() + J(1)}, (C.34) where u = rK + k cos 0, and note that the WH expansion of the exponential term on the right side was given as (C.20). Substituting in the first three terms of (C.20), the product becomes J(x', y') exp[-iKz(x', y')] = e-u202/2{J(~) + J(1) - iu((x', y')J(O) -iu((x', y')J(1) U[C2(1, y-,2]J(O) - [2(x', y') - 2]J(1)}. (C.35) The last three terms can be written as linear combinations of H(n); the H(2) and H(3) parts are omitted, resulting in the first-order expression J(x', y') exp[-iKz((x', y')] = e-22/2{J(O) + J(1) - iu(x', y')J(o) -iu(((x',y')J()} -u2((S(x', y')J(l))((x, y')). (C.36) Substituting in the definitions of J(O) and J(1) gives J(x', y') exp[-iK((x', y')] = e -a/2exp[iky'sin 0]{Fo

172 + (2 f F(C, /)((a, /)ei(ax'+6Y')dad/ -iu((x', y')Fo -iu (2)2 J F(,.)C(a, /3)dad/ -_U2 (, F(,, /)Z(a, )dad/((x', y')}, (C.37) where u = rc, + k cos 0. The magnetic current M has an identical expression with Fo and F replaced by Go and G. To repeat, Fo, Go, F, and G are the non-random coefficients of the WH expansions of the electric and magnetic surface currents. We can find these coefficient functions by writing the equations for conditions on the surface currents and applying orthogonality as before. The first two conditions are applications of the extinction theorem below and above the surface: Ei(r) + E8(r) = 0 for z < (C(x,y), (C.38) Et(r) = 0 for z > ((x, y). (C.39) The incident field Et is expressed in an eigenfunction expansion: E = e- exp[ik(y sin 0 - z cos 0)] 1 _ _ ei~'r = (27r)3 2;k2 (2r)2ei (Kr)6(ry - k sin 0)(i2k cos O)d3K (C.40) and then substituted into (C.38), resulting in (27r)e6(r)b(y - k sin 0)2k cos0 = f K x e-r' [k x J(x', y') + M(x', y')] dx'dy' (C.41) The transmitted field yields a similar equation, from (C.39): 0 = K x eir [K' x J(x', y') + M(x', y')] dx'dy'. (C.42) xJ yl y

173 The WH expansion for the electric current, (C.37) and the similar expression for the magnetic current are substituted into (C.41); both sides are multiplied by H() = 1; and both sides are ensemble averaged. Equation (C.41) leads to - e2 cos 0 = rt(eh. Fo) + r7e,(et - Fo) - e G)(e + ( Go), (C.43) where the delta functions in Et have forced Kx = 0, y = k sin0, Kz = - zk2 - K2 - -kcos0, u = -kcos0+kcos0=, K = yksin0 -kcos0. The current-exponential product similar to (C.37) for the dielectric medium case is J(x', y') exp[-izSC(x', y') = e-V2 2/2 exp[iky'sin 0]{Fo + (2r)2 F J F(, $)~(a, /)et(x'+8Y')dad,/ -iv((x', y')Fo -iv 2 jF(a, fl)(a, /)dad/ 2 -v2 (2 )2 F(a, /3)(a, /)dad ((x', y')} (C.44) where v = C' + k cos a, and a similar expression may be obtained for the magnetic current-exponential product. Substituting these into (C.42), multiplying both sides by H(O) = 1 and ensemble averaging leads to 0 = r'[y sin w(k'sin wFoy + k'cos wF,)

174 +z cos w(k' sin wFO + k' cos wFOz) -k'FO - ^-yk 'F - zk'FOz] +x(k'sin wGoz -k' cos wFo) + Syk'cos wGo - zk'sin G', (C.45) where F' = Fo- ivz2A, G' = G- ivorB, and A = (2r.)2 JJ F(a, )c(a, )dad/, (C.46) B =(2 r)2 JJ G(a, )c(a, )dad/, (C.47) and the delta functions have forced ' = 0, I = ck'sinw, = (k')2 -(2 - ( 'cosw, v = k'cosw + kcos, K' = yk'sin w + zk'cos w. Returning to (C.41), multiplying both sides by an H(1) functional ( (c, f) and ensemble averaging yields (after some manipulation) 0 = x K x [F(Kx, Ky - k sin 0) - iu(Fo - iuaA)] +K x [G(C(K, y - k sin 0) - iu(Go - iu2B)] (C.48) where =z = - - k2-+ 2 + K U = Kz+ kcos0.

175 A similar process using (C.42) yields o = k-K' x n' [F(K' - k sin 0)- iu'(Fo - iu'cA)] +K' x [G(, ' - k sin 0) - iu'(Go- iu'a2B)] (C.49) where = - /(k')2 + Kr2 + K 2 u' = cK + kcos0. The extinction theorem equations are not sufficient to completely specify the surface currents. We supplement (C.38) and (C.39) by explicitly requiring that the surface currents be tangential, resulting in two additional equations: f.nJ(x,y) = 0 (C.50) n M(x,y) = 0 (C.51) where 9 2 C U '\ 2 f1/2 n= z- - y 1+ + L a (y ) Y\Qy)J Expanding the electric current and neglecting H(2) and higher terms in the current dot product as before, we obtain 01 j ( Ir B)ff'( ( ))ei(a+yd)dad / 0 - Foz + (2r) F -o F o - (2wr)2 F (a, )() e (a) +Y)da+ d/\) -FOy (c.1r)2 F:((a, i d- + dd) (C.52) We now apply orthogonality as before: multiply (C.52) by H(O) = 1 and ensemble average, leading to Foz = (2I)2 [aFx(, f) + Fy(a, f)] C(a, P)dadd. (C.53) F"' = (C.53)J

176 A parallel procedure gives the magnetic analogue: Goz = (2r [aGx /)c +,/3Gy(a, )] (a,/)dad/. (C.54) Multiplying by H(1) functional ( (a', /') and ensemble averaging gives two more equations: Fz(a,,/3) = i(aFox + Foy), (C.55) Gz(a,/3) = Z(aG0o + Go0). (C.56) Equations (C.43), (C.45), (C.48), (C.49), (C.53), (C.54), (C.55), and (C.56) contain sufficient information to solve for F0, Go, F, and G. Once these are available, the surface currents and functions dependent on them may be obtained directly, usually in the form of integrals of the surface currents. The algebra is somewhat daunting; use of a symbolic math program such as Mathematica can ease the task significantly (although Mathematica also needed assistance from time to time). We now review the steps leading to Fo, Go, F, and G. First, notice that Az = Bz= 0 by inspection of (C.46) and (C.47) together with (C.55) and (C.56) and knowledge that c(a, /3) is even in both a and /3. Taking the x and y components of (C.48) and (C.49) gives four equations which are solved for F, Fy, Gj, and Gy in terms of ten other variables: FX(Kz, y -k sin O) =- >AAx + q->AyAy + OB B. + OBY By + toFoxFOx + qFoyFoy + rFo Foz +qGo:Goox + OG0oy Goy + OGo, Goz (C.57) Fy (Kz, %y -k sin 0) = A AXA +,AY Ay + OB BB + /B,, By +^o F+ Fo + Fo'yFoy + oFoz + 'Go. Gox + OG0o, Goy + IGo, Goz (C.58)

177 G.(K,, Ky - k sin90) = -/A., A+ -yAAy+ BBx + -IBBy +-YFox Fox + YFOy Foy + 'YF0,FOz +)Gox Gox + 2GOY Goy+ 'YGox Goz (C.59) Gy (K, r. - k sinG0) = Ur Ax + ~ AyAy + ~BBB + ~ByBy + F0, F01 + cFoy FOy + &Fo, Foz + Go, Gox + 5Goy Goy + ~GoZGOz (C.60) where the O(Kx, Ky - k sin90), 0(tcx, Ky - k sin90), 7y(ix, tc - k sin0), and c(tx, Krk sin 9) (whose arguments were suppressed) are complicated rational expressions of' k, k', IC1, /C, IC, K and. These four equations are multiplied by ~(icx, rvy - k sin 9), the transform of the normalized correlation function, and integrated with respect to rvx and' tj. Using (C.46) and (C.47), we obtain Ax = Q0ATAx + (DAYAy + lBhBx + 4OyBy + ~F0&FOx +SF'hFOf + 'Fo,ZFOz + "GoX Gox + 'Goy GOy + 4GGozGOZ (C.61) A = T ATAx + 'PA A + TBTBx + TByB y + 'I F0IFox +TFoyFOy + TIFo0 Foz + TIIGoTGOx + TIGoyGOy + TGoG,3Z (C.62) B3 = FA.Ax + FA Ay + FB, B + FBIBBy + rFF0 Fox + FFoy Fo + J'FOzFOz + FGo0TGox + rGoy GOy + rGoz Goz (C.63) By= 6A, Ax +""AAy + ZB.Bx + ZByBB + EFoX F0x +ZFo,FOy + EFo.FOz + =Go0 CGo + =Goy GOy + ZGOZ Goz (C.64) where (e.g.) A (2) Jf A(KAc - k sin9) )(tcx, r. - k sin90)dKrdc,. (C.65) X1 (27r)l

178 Exactly half of these terms are zero due to the odd symmetry of the integrands in K_,: IAy = IBs = B Foy = Fo- = IGo- = 0, (C.66) TA= = By = Fo. = Goy = -Goz = 0, (C.67) FA- = FBy = rFo, = rGoy = FGoz = 0, (C.68) Ay= ZB = Foy = ZFoz = -Go = 0. (C.69) Therefore, the four equations can be split into two uncoupled sets: Ay = ayAyA + BxBx + FoyFOy + ^FozFOz + GoGOx, (C.70) B, = rAy Ay + rBB, + rFOyFoy + rFoFz + rGGo,Go (C.71) and Ax = 4A.Ax + By+ $FoFOxx + 4GGoyGOy + 4GGozGo, (C.72) By= ZAxA + ZBBy + ByB Fo, Fo + ZGoyGoy + Goz GOz. (C.73) We can solve for the A and B components in terms of the others: Ay = fl(Foy,Foz, Go), (C.74) Bx = f2(Foy Foz, Gox), (C.75) A, = f3(Fox,Goy,Goz), (C.76) By = f4(Fo, Goy, Goz). (C.77) We now have two sets of three unknowns remaining. The horizontal and vertical incidence cases are considered separately. Let ei = e" = x. The z components of (C.43) and (C.45) give two (polarizationspecific) equations in Foy, Foz, and Gox: (Go: - iBa2v) +,[Fo cosw + (Foy - iAa2v) sinw] = 0, (C.78) kI

179 Go: +r [F sFoos + Foysin w]= 0. (C.79) To get a third equation, we examine (C.53), which is in terms of Fx, etc., which have already been solved for. Some terms of the integrand drop out due to odd symmetry, leaving the following equation, which is not specific to vertical or horizontal polarization: Fo = (A + + I)A+()Go(B + B)+(GO) + (C.80) where (e.g.) y = - (27r)2 af A,(a, /)(a,)dad/3 (C.81) and (Ay = PPAy / A (a, /)C(a, /3)dad/. (C.82) Equations (C.79), (C.78), and (C.80) make a homogeneous system of equations that has a non-vanishing determinant; the equations are therefore independent, and the only solution is therefore Foz = Gox = Foy = 0, implying Ay = Bx = 0 for the horizontal polarization case. We now need to solve for Fox, Goy, Go,, Ax, and By, which will lead to Fx, Gy, and G2, and then to Fy, F, and Gx. We have four equations: (C.76), (C.77), and the x components of (C.43) and (C.45). (C.54) will provide the fifth, similar to (C.53). This equation is also not polarization-specific: Go = (PA. +A.) + (rB, +B) + (r^Fo +F ) + (rGOY + G) + (rGo. + Goz), (C.83) where (e.g.) FAX = (27r)2 a 7AI (a, /)(a, /3)dad3 (C.84)

180 and A = f2 )21 (,A) (a, ).(a, ()dadC.) We now have five equations in five unknowns, which may be solved for Fo0, Goy, GoZ, AX, and By. Fx, Fy, G, and Gy follow from (C.57), (C.58), (C.59), (C.60), and the surface currents follow directly. The vertical polarization case is complementary. It should be noted that, while the normalized surface correlation, c(7r, ry), and its transform, c(a, /3), are still unspecified, it is better in practice to use an analytical expression for c(a, f) in the hope that some of the two-dimensional integrals (DA. etc.) might be solved analytically (not likely) or at least reduced to one dimension. Such a reduction is extremely advantageous from a numerical analysis viewpoint. For example, in the Gaussian correlation case, the two-dimensional integrals were converted into radial coordinates and reduced to one dimension analytically, with integrands in terms of Bessel functions. The power-law case has no similar reduction. The use of a symbolic algebra program like Mathematica is highly recommended due to the number of algebraic manipulations and the high probability of algebraic errors. Use of an algebra package alone does not guarantee usable solutions, however. The author found that Mathematica often needed very specific hints before it could reduce complicated formulas to simplest form. Also, the direct solution of the three-equation linear system in which the coefficients were rational functions led to incredibly complex solutions. The problem was avoided by introducing dummy coefficients and implementing Cramer's Rule via explicit Mathematica commands.

181 C.2 Derivation of the Backscattering Coefficient a~ The bistatic scattering cross section in direction (08, /3) for the scattered field with polarization es is given by k2 (9, S) = 4A (s.) ~-(e.)2) (C.86) where A is the illuminated area (which drops out of the final expressions) and E is the amplitude of the scattered far field: ikr E(r) = ik (8O,s), (C.87) 47r where ~(0e, s ) = -r x J[M(x', y') + rAI x J(x', y')] exp(-ikr r')dx'dy'. (C.88) Defining unit vector r in the direction of the scattered field, r = x sin 08 cos Os + y sin 80 sin Os + z cos 98, (C.89) and vector r' to a general position on the surface, r'= x+x' + y+ (x', y'), (C.90) we can show that, for scattering in the plane of incidence (&S = 7r/2), eh = x (C.91) e = -ycos O + z sin (C.92) r = ysin 0 + cos 0 (C.93) r r' = ysin 8 + ((x', y')cos 0. (C.94) The expression inside (C.86) may be written e. ( 7 2 )-(e *(0(9, )) e e (O E( 2 e \i

182 -= Jf e J(x', y') exp[-ik((zx' y) cos Os]e-iY'sindx'dy + ffe ' rx M(', y') exp[-ikC(zx', y') cos s]e-ikY' sin' dx'dy' +7 JJ e'S (J(x', y') exp[-ik((x', y') cos 0S])e-iky sin dx'dy' - fJe x (M(x', y) exp[-ik((x', y') cos O])e-ik' sindx'dy'. (C.95) If the surface electric current J(x,y) is written in a first-order Wiener-Hermite expansion as in (C.31), then the current-exponential product in (C.95), written to first order, becomes J(x', y)e-ikcos'C(I',Y') = e-ikp()isin {Fo + (2r)2- (aI, 3)((a, )e(Yx'+3Y')dadf}, (C.96) where p = cos 0 + cos 8O, and the ensemble average of the product becomes (J(x', yl)e-ikcosO C(x',y)) = e-k22/2iky' sine {FO ikp )2 F(a, F)c(a, a/)dad/. (C.97) Substituting these expressions and their magnetic current counterparts back into (C.95) and simplifying results in es.~(s, )-(es e(o8 )) - - f e - 'Y{[-r,* Fo + & ir x Go][To(', y') - (To(x', y))] + (2r )21 J/[- es F(a, 3) + es rx G(a, /,)] [[Tl (cr,; x', y') - (T (ax, 3; x', y'))]dad/}dx'dy', (C.98) where q = sin 0- sin (C.99) To(x', y') = exp[-ikp(x, y')] (C.100)

183 (To (x', y')) = exp (-sp'l/2) (C.1O1) T1(a, /3; x', y') = To(x', y')~(a, 3) exp[i(ax + 3y)] (C. 102) (Ti (a,/3; x',y')) = -ikpa(a,/3)oc2e- 8p2/2 (C. 103) s = k2 (C. 104) The next step is to calculate the ensemble average of the squared absolute value of (C.98). The following ensemble averages are calculated separately: ([0-(o]T*-(To*)]) = _ep2 [e8p2C(XX',Y~y) -1,(C.-105') ([T - (TO)][ T, * - (T 1)]= -ikpor2a(c, /3)C_,p2 {esp2c(X-X'YiY') eia(X-X')+i,(Y-Y') - i] + I} (C.106) (IT, - (,]T*- (T1*)]) = (27r )2o2 (a, /3)6(a - a'60- /31),:(cx-a1X'+iOY —3'Y') 7-SP2e8P2c(~X'"Y-Y') _S22 fl')e-p2{ +ea(-'+,'YY)- ei(a+a')( X-')+i(,8+3O')(Y-Y')j} (C. 107) Let Qo=-Iie.Fo~&$.^xGo (C.108) and Q(a, /3) = -?es- F(a, /3) + e&- r' x G(a, /3.(C.109) The ensemble average of the squared absolute value of (C.98) is the sum of four integrals: (J^,s _ E(,97,!) _-,, e(O.9, w))12) =

184 jJ e-ikq(Y-y')QoQoQ*([To(x, y) - (To(x, y))] [To*(x', y') - (To*(x', y'))])dxdydx'dy' + (2)2 j e ikY QoQQ* (a/')([To(x, y) - (To(x, y))] [Ti*(Ca, /'; x', y') - (T*(a', /'; x', y'))])da'd3'dxdydx'dy' + 1(2 )2 e-ikq(y-Y')QO*Q(a, p)([To*(x, y') - (To*(x, y'))] + (2r)2 JJJJJJ - [Ti(a, /; x, y) - (T1 (a, /; x, y))])dad/3dxdydx'dy' 1 eikq(y-Y') Q(a, /)Q* (a', 0/')([T1(a, /3; x, y) - (T1 (a, /3; x, y))] [Ti*(a', /'; x', y') - (T*(a', /3'; x', y'))])dad/dca'd/3'dxdydx'dy'. (C.110) After much algebra, some changes of variables, and some intermediate integrations, the bistatic cross section in the plane of incidence is given by the following general expression: k2 2 a~(0o, 2) = -Sep { IQol2Mo 01 ": It/' +2kpa2r 2 12 Im{QoQ*(a, l/)}C(a, /3)MI (a, 3)dad/3 o 2 + (27r)2 || IQ(aI v) 12i(a )M2(a, I )dad/ s0'2p2 f//f -(2i) ]] Q(a, /)Q*(a, /')c(a, /)c(a', /3') M3(a, /; a', 3')dad/3da'd'}, (( D.111) where Mo Mi(a, /3) M2(a, /) M (a, /; a', /3') F(a,/ ) f e-ikqv[ePc(uv) - l]dudv = (27r)26(a)6(/3 - kq) + J~(a, /3) - M m M(, /3) + Mo = M(a,/3)+ (a', /') - M (a +, / + /3') = e- i [esp2c(u'v) -_ ]eiat+i',vdudv. (C.112) (C.113) (C.114) (C.115) (C.116)

185 Eftimiu shows that (C.111) reduces to the physical optics and small perturbation expressions in the appropriate limits. For example, consider a long correlation surface, where (a, /) = (27r)26(a)6(). (C.117) All terms in (C.111) go to zero except the first, leaving k 2 aO~(O, ) = - e-P2IQoIMo. (C.118) For horizontal incidence, ei = e, and Qo =-7Fo - cos Goy + sin O8Goz. (C.119) The integrals 4Ax, etc., appearing in the expressions for For, Goy, and Go, become much simpler in this limit (through the appearance of (C.117) in the integrands), and we can show that, for the backscatter case (80 = -0), -2k' cosw cos 8 rl Fox (C.120) k' cos w + k cos 0 G 2k cos 0 kc' cos w + k cos (C121) Go, = 0; (C.122) therefore ( k' cos w - k cos 0\ Qo=2cos k'cosw-+kcos ) (C.123) and for the backscatter case, a~(- ) = ke-2e 2cos cos2 0IR 1Mo. (C.124) For c(T,, ry)= 1 and 08 = -0, Mo= eik2sinOv [e4k22cos20 - 1] dudv, (C.125) and (C.124) reduces exactly to (7.5), the physical optics solution.

186 C.3 Evaluation of (C.111) for the Power-Law Case The expression for the scattering coefficient, (C.111), is evaluated by substituting the explicit expressions for Qo and Q and performing the indicated integrations. Eftimiu remarks: "This calculation is, by far, easier said than done" [16]. Use of a power-law roughness spectrum (and the associated correlation function for which there is no analytical expression) only serves to increase the difficulty. The evaluation of (C.111), shown here for HH polarization, generally follows the procedure used by Eftimiu. We begin by calculating all the integrals similar to DA, using the normalized roughness spectrum SZ(Q,/3) C(a,) Z 2', (C.126) where Sz(a, 3) is given by (2.9), modified for use with arguments in rad/m instead of m-1. The next step is to expand the exponential containing the normalized correlation function in (C.116) into a power series. We therefore write oo 2 F(ka, kf) = z (sp 2) cn(u, )ei[kau+(k-kq)]dudv (C.127) n=l n! Each term of this function is calculated by filling a 2D array with values of c, using an inverse FFT to calculate c, raising c to power n, and transforming back to the frequency domain, resulting (after convergence) in a table of values for F. Values not in the table may be found by interpolation. The expression for a~ is divided into two parts: a~ = a~7+ ~, (C.128) where a(0, ) = -e Q{lo2Mo

187 +2kpa2 f( f ) Im{ QoQ*(a, /)}c(a, /3)Ml(a, P)dad/ + ()2 IQ(a, l )12(a, 3)M2(a, )dad, (C.129) and the 7~ part contains the term in (C.111) with the quadruple integral involving M3. Evaluation of this integral is impractical in the power-law case. Eftimiu remarks that this term has little or no effect in the (Gaussian) cases he has considered, so its omission here seems reasonable. Rewriting a~ as a~ after some manipulations: o(0e, ~ p) { e-sp IQo (0,)+ (, ) (0,q) [pQo + Q(0, kq) 2 +2kpo2 (2 )2 Im{QoQ*(ka, k/3)}c(ka, k/) *[F(ka, k/3) - F(O, O)]k2dad/ + (2)2 JJ IQ(kQ, k I) l2(ka, k/3)F(ka, k)k2dadS }, (C.130) where '- (sp2)n 'ei[ka+(k-kqv] -= Cs! |I (u, v)et +( v ) ]dudv. (C.131) n=2 n! For HH polarization, Qo = — Fo -cos 9Goy + sin 08Goz (C.132) Q(ka, ko3) = -7IFP(ka, kfl) - cos 08Gy(kac, k3) + sin OsGz(ka, k$). (C.133) The scattering coefficient in (C.130) was computed by the program plhang4.f. The integrals are rather messy because the integrands are products of radially symmetric functions with offset centers. However, the most taxing calculation by far was the summation for F(ka, k/$). This function was calculated by the program convbl5.f; values were stored in files for use in the other program. The sum for F

188 did not converge when the entire spectrum from (5.8) was used. Use of a bandlimited spectrum formed by enforcing a low-end cutoff of (10A)-1 allowed F to converge, but not always quickly (as many as 18 terms were required). The slow convergence is believed to be accurate, however, because F is essentially an n-fold self-convolution of c, which has a very high peak. Results of the PWH calculations are shown in Section 7.5.

189 APPENDIX D FLUCTUATION STATISTICS OF MILLIMETER-WAVE SCATTERING FROM DISTRIBUTED TARGETS

90 268 IEF.E TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 'S. NO. 3. MAY 1988 Fluctuation Statistics of Millimeter-Wave Scattering From Distributed Targets FAWWAZ T. ULABY, FELLOW, IEEE, THOMAS F. HADDOCK, MEMBER. IEEE, AND RICHARD T. AUSTIN. STUDENT MIE.BER, IEE Abstract-The applicability of the Ravleigh fading model for characterizint radar scattering trom terrain is examined at 35 (ilt for both backscattering and bistatic scattering. The model is found to be in excellent agreement with experimental observations for sinule-frequency observations of uniform targets such as asphalt and sno,-co\sered ground. The use of frequency aleraling to reduce sirnali fadintl sariations was examined experimentally b) se-ceping the radar siciual from 34-36 (;11z in 401 steps. [lhe results shoa that the formull;lti based on the RaK leiih lmodel relatingl the reduction in sitinal Iluctuationl to the bandwidth used provides a reasonable estimate of' the iml proement provided by frequency averaging. I. INTRODUCTION T O fide. as defined in Webster's dictionary [11] is 'to change gradually in loudness, strength, or visibility, when used (in connection) with a motion picture image or an electronics signal." In radio communications [2], signal fading refers to fluctuations in the received signal caused by multipath interference, and in radar sensing of terrain the terms fading, scintillation, and fluctuation have all been used interchangeably to describe random-like intensity variations corresponding to signals backscattered from cells at different locations on a distributed target [3. pp. 463-495, 1803-18041. If the radar is of the imaging type, the random variations produce a "speckle" pattern or appearance on the image, which complicates the image interpretation problem and reduces the effectiveness of information extraction algorithms. Consider, for example, the two image segments shown in Fig. 1. These two segments, one of which corresponds to a corn field and the other to a forest parcel and which were part of the same strip of X-band radar imagery. have different average tones, exhibit significantly different tertures, and both exhibit large pixel-to-pixel intensity variations. The average tone of an image is the average value of the image intensity for all pixels contained in that image. (Each image segment contains approximately 104 pixels.) This average tone is proportional to the average received power, which, in turn, is directly proportional to the backscattering coefficient oa of the imaged target. In Manuscript received October 6. 1987. The authors arc with the Radiation Laboratory, EECS Department. The University of Michigan. Ann Arbor. M1l 48109-2122. IEEE Log Number 8820738. I -300- 3 m —0 (a) (b) Fig. 1..-band SAR images ot (a) a corn field and (b) a forested area. Note the textural differences between the two inmages. other words, a0 of the imaged target is. by definition, the mean value of the random process characterizing the intensity variations in the image. Texture refers to the low spatial-frequency variations of intensity across the image [41; the corn field, being more spatially uniform than the forest parcel, exhibits the same type of random variations in all regions of the image, whereas the image of the forest parcel contains "clumps" of dark and bright regions, on which the random variation is superimposed. If we adopt the strict definition that the concept of "a backscattering coefficient for a distributed tarEet' is meaningful only for targets with uniform electromagnetic properties, then texture becomes the spatial variation of a') from one region of an image to another. In the case of the forest parcel, these variations are related to the spatial nonuniformity of tree density. Unlike textural variations, which may or may not have specific directional properties and which are governed by the spatial variation of the target scattering properties relative to the dimensions of the radar resolution cell., the random variations that give the image its speckled appearance are due to phase-interterence effects and are a characteristic feature of the scattering pattern for any distributed target (provided the target satisfies certain conditions, as we shall discuss later). Image speckle is simply a visual manifestation of fading statistics, which is the central topic of this paper. Thus. there are three types of intensity variations that one may observe in a radar image: 1) variations in average tone from one distributed target (such as a bare-soil field) 0196-2892/88/0500-268S01.00 ~ 1988 IEEE

191 ULABY. HADDOCK et al.: FLUCTUATION STATISTICS OF mm-WAVE SCATTERING 269 to another (such as a forest parcel), 2) textural variations from one region of a distributed target to another, and 3) random fading variations at the pixel-to-pixel scale. These variations are governed by different processes and are characterized by different probability density functions (pdf's). In some radar applications, these three types of variations are lumped together, treated ais a single variation, and characterized as terrain clutter. To determine the statistics of the clutter random variable for a given terrain type or geographic area, the area is imaged and then a pdf of the received voltage or power is generated. Next, the data is tested against theoretical pdfs to determine which fits best. Such an empirical approach may produce a statistical description appropriate to the imaged area, but it has some severe limitations. The empirically generated pdf is, in essence, a convolution of the three pdf's characterizing the three types of variations referred to above. Hence, it is both target-specific and sensor-specific. It is target-specific in that it pertains to the specific mix of terrain categories and the specific conditions of those categories at the time the radar observations were made. Most terrain surfaces exhibit dynamic variations with time of day, season, and weather history. The pdf is sensor-specific because one of the underlying variations, namely that due to signal fading, is governed by the detection scheme used in the receiver (linear or square-law) and the type of filtering or smoothing technique employed in the signal processor. Filtering techniques are used to reduce fading variations; they may include spatial averaging and/or frequency averaging schemes and may be performed coherently or incoherently [4]-[9]. To characterize the fading statistics associated with a terrain surface of uniform electromagnetic properties, the usual approach is to model the surface as an ensemble of independent, randomly located scatterers, all of comparable scattering strengths. Such a model leads to the result that the amplitude of the backscattered signal is Rayleighdistributed [3, pp. 476-481]. If the return is dominated by backscatter from one or a few strong scatterers, the fading process is characterized by the Nakagami-Rice distribution [10]. Some experimental observations support the Rayleigh behavior [4], [11], [12] while others, particularly those measured for complex terrain categories, are in closer agreement with the lognormal or the Weibull pdf's [13]-[17], or other more complicated distributions [18]. The purpose of this paper is to: 1) examine the applicability of Rayleigh fading at 35 GHz for both backscattering and bistatic scattering from uniform terrain media, 2) examine the statistics associated with the use of frequency averaging to reduce fading variations, and 3) determine if the statistical character of the backscatter is affected by the size of the ground cell (antenna footprint) illuminated by the radar. To this end, both experimental measurements and theoretical analyses were performed. A Fig. 2. The illuminated area A contains N, randomly distributed scatterers. II. RAYLEIGH FADING STATISTICS A. Underlying Assumptions The Rayleigh fading model used for describing radar scattering from an area-extended (distributed) target is essentially the same as the model used for random noise and is based on the same mathematical assumptions. A review of these assumptions will prove useful in later sections. The sketch shown in Fig. 2 depicts a radar beam illuminating an area A of an area-extended target. The illuminated area contains N, point scatterers designated by the index i = 1, 2, * * *, Ns. For simplicity, we shall confine our present discussion to the backscatter case. The field intensity at the input of the receiving antenna due to backscatter by the ith scatterer may be expressed as Ei = KiEio exp [j(wt - 2kri + oi)] (1) where E,o is the scattering amplitude and i0 is the scattering phase of the ith scatterer; ri is the range from the antenna to the scatterer; k = 27r/X is the wavenumber; and Ki is a system constant that accounts for propagation losses to and from the scatterer, antenna gain, and other radar system factors. The expression given by (1) may be abbreviated as Ei = KiEio ej' where (2) (3),i = wt - 2kri + O, is the instantaneous phase of Ei. Assumption 1: The scatterers are statistically independent. This assumption allows us to express the total instantaneous field due to the N, scatterers contained in the area A as a simple sum NI E = Z KiEoei' i, 'I (4) and it implies that interaction effects between adjacent scatterers may be ignored. Assumption 2: The maximum range extent of the target Ar = | ri - rj Ima is much smaller than the mean range to the target area A, and the antenna gain is uniform across A. This allows us to set Ki = K for all i. For convenience,

270 Im E 192 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 26. NO. 3. MAY 1988, E sino E 0 ---.. --- —----- i0 2 E cos. ' - * Re E Ee sin 0 e Fig. 3. The vector E is the phasor sum of N, fields. we shall set K = 1. Hence N, E= Z Eoe'. (5) The total field E is a vector sum of N, phasors. If we express these phasors graphically (Fig. 3) with the first one starting at the origin and the successive one starting each at the tip of the preceding one, the resultant is a vector from the origin to the tip of the last phasor. The length of this vector and its phase angle are denoted E, and 0, respectively. That is noiselike statistics do not apply and the statistics developed by Rice [20] for one or more large signals contained in a background of noise should be used instead. Use of Assumptions 3-6 can be shown to lead to the following properties [3, p. 479]: p(E,) = exp (-E/2s2), E, > 0 p(0) = 1/(2r) E.= 2 S E' = 2s2 (7) (8) (9) E = Eeej. (6) (10) Assumption 3: N, is a large number. This assumption allows us to use the central-limit theorem, which, in turn, allows us to assume that the x- and y-components of E, Er, and E,., are normally distributed. However, it can be shown through computer simulation that this condition can be satisfied (approximately) for N, as small as 10. (The same conclusion was reached by Kerr [19] in the 1940's.) Assumption 4: The scattering amplitude Ej0 and the instantaneous phase Oi are independent random variables. This condition is easily satisfied if E,o is independent of the range ri, which would be the case if the scatterers are randomly distributed in range. Assumption 5: The phase,i is uniformly distributed over the range [0, 2r ]. To satisfy this condition it is not only necessary that the scatterers be randomly distributed in range, but the maximum range extent of the target Ar must be several wavelengths across also. Assumption 6: No one individual scatterer produces a field intensity of magnitude commensurate with the resultant field from all scatterers. In other words, the field E is not dominated by one (or few) very strong scatterer(s). If this condition is not satisfied, the Rayleigh where p(E,) and p(0) denote the pdf's of E, and 0, respectively, E, is the ensemble average (mean value) of Ee, and s is the standard deviation of E. and Ev. Equation (7) is known as the Rayleigh distribution. B. Output Voltage 1) Linear Detection: If the receiver uses a linear detector, its output voltage VL is directly proportional to E, VL- KIE, - E = KIK2(a~)/2f = VLf (11) where K, is a system constant, K2 relates the mean field Ee to the backscattering coefficient of the target a~ (actually, a~ is directly proportional to E2, but E2 = 4/'ar ~2 ), and fis the normalized fading random variable given

193 ULABY. HADDOCK cr a(t. FLUCTUATION STATISTICS OF mm-WAVE SCAT'I1RING 271 by o, —. f= E,/E,. (12) Using the relation p( E,) dE, = p( f ) df, we obtain -rf p(f) = fexp ( -f /4), f>_ 0 7= 1 Sf = s,,/ VL = 0.523. IL 2 -a o a 5 (13) (14) ( 15) 08 -6 - 04 -0.2 -00 - exponential ',-DP p(F) Rayleigh PDF p(1) 0 1 2 3 4 Because the output voltage VL is a product of its mean value VL = K, K2(o~) )2 and the random variable f, the process is sometimes referred to as a nultiplicative noise model. 2) Square-Law Detection: The voltage output of a square-law detector is directly proportional to the power of the input signal rather than to its field intensity E,. Thus Vs = K, P E2 = K K4Ee f or F (a),, -I LL. -0 v %L 0 8 -6 0 4 0.2 P(F<Fo) / P(f<fo) = K K4 KS "F where F = E,/E~ is the normalized fading random variable for power pdf characterizing F is the exponential distribution I 4801 p(F) = e-F, F 0 with F= I (16) (17).The [3, p. n r, -20 -10 0 f. (dB) or F, (dB) 10 I. I I. I I I. I.. - I I *..01. 1 10 F. I' ' ' I3.0 I.. 1.3 1.0 3.0 1. (18) (b) Fig. 4. Plots of (a) probability density functions and (b) cumulative distributions for fand F. (19) and F = Sp/P = 1. (20) C. Interpretation What do these statistics tell us? To answer this question we start by examining Fig. 4(a), which shows plots of p(f) and p(F) for the Rayleigh and exponential distributions, respectively, and Fig. 4(b), which shows the corresponding cumulative distributions. We observe that the range of fading associated with these distributions is very large. That is, if one takes a single sample of the signal from a Rayleigh-distributed or exponentially distributed ensemble, one has very little chance of selecting a value close to the mean. To illustrate this with a specific example, according to the Rayleigh distribution in Fig. 4(b) the value off that exceeded 5 percent of the time is 1.95 (relative to the mean) and that exceeded 95 percent of the time is 0.25. In decibels, these levels correspond to +5.8 and - 11.9 dB, respectively. If we select a sample at random, the probability is 90 percent (95 to 5 range) that its value will be within the range extending from 11.9 dB below the mean to 5.8 dB above the mean. We may think of this as the 90-percent confidence interval associated with our measurement. The important point to note here is the fact that this interval ( 17.7 dB) is very large indeed. The situation is not much different when square-law detection is used; the 5- and 95-percent levels of the cumulative distribution for the exponential pdf are +4.8 and -12.9 dB, also totalling to 17.7 dB. Now let us illustrate the fading behavior with measured data. Fig. 5(b) presents a trace of radar backscatter measurements made by a 35-GHz truck-mounted scatterometer as the truck was driven across an asphalt surface with the radar beam pointing downward along the aft direction at an incidence angle of 40~ relative to normal incidence (Fig. 5(a)). The antennas were mounted atop a telescopic boom at a height of 10.3 m above the asphalt surface. The sampling rate was such that the footprints (on the asphalt surface) corresponding to adjacent samples were totally independent (no overlap). More detailed information on the system and measurement procedure is given in Section IV. The vertical axis in Fig. 5(b) represents F, the ratio of the received power to the average value computed for all 1000 measurements, expressed in decibels. (It is assumed

194 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 26. NO. 3. MAY 1988 (a) D. Independent Samples To improve the uncertainty of a radar measurement of the backscatter from a terrain surface, it is necessary to average many independent samples together. An easy way to increase the number of independent samples N contained in an estimate of the radar backscatter is through spatial averaging, which amounts to trading spatial resolution for improved radiometric resolution. Other ways to increase N are discussed in Section III. 1) Linear Detection: If N randomly selected samples of a Rayleigh-distributed voltage VL are averaged together, the average value VL, has the following properties: I, VLV= - L K,(i = KK(o)(:0) I j - ' r Variation of F(dB) with spatial position 10 0. *' 'o — 30 40 -50 0 200 400 600 800 1000 Posltlon Index (b) Fig. 5. The sketch in (a) shows how the measurement of the backscatt from asphalt (shown in (b)) were acquired. The incidence angle was the platform height 10.3 m. and the polarization VV. The mea' backscattering coelhcient (corresponding to P ) was -5.25 dB. Comparison of Measured and Exponential pdf's 0.2 r ering 40~. sured = K, K2((a) f where we defined NI V (21 ) (22) as the fading random variable corresponding to the average of N independent samples. Its properties are N = 1 LL -1 'a 0 1 e-F. F n 1 0.523 S/N = (23) O P(F) (meas) - P(F) (epon) and its pdf may be obtained by N-successive convolutions of the Rayleigh distribution (13). Plots ofp ( fN) are shown in Fig. 7(a) for several values of N. As expected. as N increases the distribution becomes more peaked and narrow (the standard deviation decreases as N-'/2 ). 2) Square-Law Detection: If the receiver uses squarelaw detection n n oll!.. - *S * *#0, F Fig. 6. Comparison of the measured pdf's with the exponential pdf. The quantity F = 0.2 is the bin size. that the mean of 1000 independent samples is a good estimate of the true mean.) We observe that I) The measured values of F extend over a range of 50.2 dB, and that 90.8 percent of these data points are within the -12.9- to +4.8-dB range (which corresponds to the 90-percent interval for the exponential distribution). 2) The standard deviation SF = 0.97, which is in close agreement with the value of 1 predicted by (20). 3) The measured pdf of F closely resembles the exponential distribution (Fig. 6); an acceptance hypothesis test using the chi-square goodness of fit test shows agreement with a probability of 86 percent. VSN = K3 K4K5a F with (24) (25) I N FN = - Z F,. Ni=I The mean value of FN is 1, its standard deviation is SP,v I "F = p h (26) and its pdf is a x2 distribution with 2N degrees of freedom [3, p. 1914] FN NN e-NFS p(FN) = Ir _ _ O. e swn in Fi. 7(b)! Plots of p(FN) are shown in Fig. 7(b). (27)

195 ULABY, HADDOCK et at.: FLUCTUATION STATISTICS OF mm-WAVE SCATTERING 273 2.5 -20 -N = 10 1.5 1.0 N= 4 * N = 4 \ 0.5 -, 0.0 0 1 f 2 N - (a) 3 1 2 1 0 -- 06 0.4 0.2 -0.0 FN- 2 (b) Fig. 7. Probability density functions for N = 1, 4, and 10 for (a)f, (linear detection) and (b) F, (square-law detection). E. Applicability of the Rayleigh Model Does the Rayleigh fading model provide an appropriate approach for characterizing the statistics of radar backscatter from terrain? The answer is a qualified yes. If the assumptions underlying the Rayleigh fading model are reasonably satisfied, the available experimental evidence suggests that the Rayleigh model is quite applicable [4], [11], [12]. Terrain targets satisfying the Rayleigh assumptions include bare ground surfaces, agricultural fields, dense forest canopies, and snow-covered ground. In all cases the target has to have stationary statistics, which requires that its "local-average" electromagnetic properties be uniform across the extent of the target. Rayleigh fading is inapplicable for a sparse forest observed by a high-resolution radar because the high spatial variations in tree density at the scale of the radar resolution violate the stationarity assumption. Thus, a very important parameter governing applicability of Rayleigh statistics to backscatter from terrain is the size of the radar resolution cell relative to the spatial frequency spectrum characterizing the scattering from the terrain target under consideration. An urban scene is another target class/condition for which Rayleigh statistics may not apply. If the resolution cell size is such that the backscatter is likely to be dominated by the return from one or a few strong scatterer(s), such as a building or a comer reflector formed bv two intersecting flat surfaces, the Rayleigh pdf is no longer applicable. In a recent study on image texture [4], Seasat SAR data was examined for five land use categories in a test site in Northeastern Oklahoma. Comparison of pdf's based on the data from the digital SAR image with the Rayleigh pdf revealed a good fit between data and theory for backscatter from a lake surface, a fair fit for grasslands and cultivated terrain, and poor agreement for forests and urban areas, particularly for the latter. III. WAYS TO INCREASE THE NUMBER OF INDEPENDENT SAMPLES According to the preceding section, if a radar is used to measure the backscattering from a uniform. randomly distributed target with backscattering coefficient a0, the voltage observed at the receiver output will be proportional to ( ao)", with n = 1 /2 for linear detection and n = I for square-law detection. However, associated with the measurement process there will be a multiplicative error represented by the random variable fN (for linear detection) or FN (for square-law detection). These random variables both have means of I and standard deviations proportional to N-1/2. Hence, the key to improving the precision of the measurement process is to make N as large as possible. Fundamentally, increasing N is equivalent to trading off spatial resolution for improved radiometric resolution. This statement is true when discrete measurements (corresponding to discrete resolution cells) are averaged together after detection, as well as when the averaging process is an integral part of the detection process (as we shall discuss later). A. Spatial Averaging 1) Discrete Samples: If N measurements corresponding to statistically independent nonoverlapping footprints are averaged together, then the number of independent samples characterizing the average value is simply N. Statistical independence requires that the spacing between adjacent footprints be greater than the spatial correlation length of the random surface L,. Thus, reflections from two nonoverlapping footprints on a very smooth surface are not considered independent because the correlation length of a smooth surface is very long (it is infinite for a specular surface). Conversely, for a random surface the returns from two footprints may be considered independent even if the footprints do overlap, provided that the spacing between the centers of the two footprints is greater than a certain distance which we shall call the fading decorrelation distance L,. Expressions for Ld are given in succeeding sections for specific antenna pointing configurations. In all cases the condition Ld > L, has to be satisfied in order for the samples to be statistically independent. 2) Continuous Averaging in Azimuth: Consider the an

196 IFEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 26. NO. 3. MAY 1988 274 z decorrelation distance is simply L, =- l,./2. (30) The result given in (29) is equally applicable to a pencil-beam scatterometer and to a fan-beam side-looking real-aperture-radar (RAR). In the case of a side-looking fully focused synthetic aperture radar (SAR), the Doppler bandwidth is used to improve the azimuth resolution from r,. = 3,.R to the resolution r, = 1,./2 corresponding to a synthetic aperture of length I, = 3vR. Thus, for the fully focused SAR Iv = rs/(1//2) = I. Looked at another way, N,. represents the degradation in spatial resolution from the best achievable (1,./2) down to r,,. 3) Continuous Averaging in Range: For a narrow pencil-beam scatterometer traveling in the x-direction, consideration of the time-bandwidth product leads to H J/ N, - rr/Ld and (31) (32) Ld - (l,/2) sec2 0 Fig. 8. Antenna with effective beamwidth I, illuminating a target at range R and incidence angle 6. tenna beam shown in Fig. 8; the boresight direction is in the x-z plane, pointing at an angle 0. and the effective beamwidth is 3,. in the v-direction. The antenna is moving along the v-direction (azimuth) at a velocity u, the nominal range to the antenna footprint is R. and the nominal azimuth resolution (width of the footprint in the v direction is) r R =,R. (28) If the radar output voltage is recorded as a function of time as the beam traverses the ground surface at the velocity u., the beam performs a form of continuous averaging equivalent to low-pass filtering. From considerations of the time it takes to travel over a distance r, and the Doppler bandwidth of the signal backscattered from the illuminated cell. it can be shown 13. pp. 585-586] that the output voltage represents an average of N,-equivalent discrete independent samples, and that N, is given by the approximate expression, - r,/(l/2) (29) where,. is the length of the antenna along the v-direction. The above result, which is independent of u,, may be interpreted as saying that the fading signal decorrelates whenever the antenna moves a distance 1,./2 in the v-direction. and therefore a resolution cell of width r, contains r,/(,./2) independent samples. Thus. the fading where r, is the ground resolution in the x-direction, and S. is the height of the antenna in the elevation plane. B. Frequency Averaging The criteria used to decide whether or not a pair of signals V, and V2 backscattered from two ground footprints may be treated as statistically independent observations is based on the magnitude of the correlation coefficient between them, p( V,, V,). If, on the average, p is smaller than some specified value, such as 0.2, the two observations may be regarded as statistically independent. Decorrelation is a consequence of differences in the instantaneous phases of the scatterers present in the observed cells. The phase of a given scatterer, as given by (3) 0i = wt - 2kri + O, 47r = ct - - vr + 9 c may be changed by altering the range r, between the scatterer and the antenna, or by changing the wave frequency v. Birkemeier and Wallace [21] derived an expression for the correlation function for two signals (one at frequency v, and the other at frequency v2) scattered from the same randomly distributed target as a function of the illumination geometry and the frequency separation Av = v, - v,. If Vs is the output voltage after square-law detection (i.e., Vs = KP, where P is the input power), the autocorrelation function for Vs ( v ) and Vs (v, ) is given by R(v,, Pv,) = Vs(v) Vs(v,) = P(,v) P(v,) (33)

197 tlI.ABY. HADDOCK <(r hi FLUCTUATION STATISTICS OF mm-WAVE SCATTERING 275 where K has been set equal to unity for convenience. For the randomly distributed target depicted in Fig. 9, Birkemeier and Wallace 121] argued that the process is stationary; i.e., R(v\, v,) = R(Av), and showed that the autocovanance function, defined as is given by R,.(Av) = R(Av) - P R,.(,v) = P2 sin av L i cAv (34) (35) D = r sin 0 x. L.. — where P = P( v) P(P v) is the mean value of the input power (assumed constant over the frequency separation Aiv), and 2 7rD 2r = = - r, sin. (36) c c The correlation coefficient is the normalized autocovariance function R?,.(Av) (sin aA \ p(Av) -= R.(O) sin c v) (37) R,,(O) \ cuA The two signals P( v ) and P( v ) may be regarded as statistically uncorrelated, and therefore independent, if the separation Av corresponds to the first zero of p(Av), which occurs at crAv = 7r. This was called the critical frequency change by Birkemeier and Wallace 121], but we shall refer to it as the decorrelation bandwidth Atd, and it is given by Fig. 9. Backscattering geometry for an illuminated cell with ground-range dimension r,. This result will be evaluated in Section V through comparison with measured data. If cB >> 1, the term / B is negligible over the region where the autocovariance function is of significant size in the integrand, which allows us to integrate the function analytically and obtain the approximate solution 2D N -B C = B/Av. (43) Here D is the slant-range resolution of the radar system. We may show the equivalence of the above result to the chirped pulse-radar case (as in a RAR or SAR) by noting that B is the chirp bandwidth and 2D/c is the de-chirped pulse length r. Hence N Br c 150 Av, = 2D - D MHz 2D D (38) = B/B, (44) with D in meters. For continuous integration over a swept-frequency bandwidth B extending from v, to v, the variance of P(B) = P(v) dv (39) is given by [201 s2(B) ( R,.() dt (40) where t is Av. Use of (35) in (40) leads to sP(B) = B- (- BI ( ) d (41) The effective number of independent samples realized as a result of frequency averaging may be obtained by relating the variance of P to its mean value as in (26) / p 2\2 = BB ~S ( - (\ sin ar) a' 2 d- (42) 2~ ~ o B ot - where Br = I /T is the receiver bandwidth. If the transmitted pulse is de-chirped in the receiver to obtain the narrowest possible pulse length, the receiver bandwidth Br has to be equal to the modulation bandwidth B. Hence, N = 1. However, if it is desired to have N be larger than 1, the pulse may be de-chirped only partially, thereby using the excess bandwidth to provide frequency averaging. This is referred to as coherent frequency averaging [8], in contrast with incoherent frequency averaging wherein the averaging operation is performed after the detection and sampling operations. That is, full de-chirping is performed to retrieve the best possible range resolution possible, and then after the image is produced, several range pixels are averaged together to increase N. IV. EXPERIMENT DESCRIPTION Two types of experiments were conducted in support of this study, one involving backscatter measurements using a truck-mounted platform and another experiment involving bistatic scattering measurements conducted in the laboratory. To maintain continuity in this presentation, only the backscatter measurement system will be described in this section, and description of the bistatic configuration will be deferred to Section VI.

198 ti-EE TRANSACTIONS ON GEOSCIE-NCE AND REMOTE SENSING, VOL.. 26. NO. 3. MAY 1988 The backscattered data analyzed in succeeding sections was measured by a.35-GHz scatterometer that was mounted on a truck-mounted telescopic boom as depicted in Fig. 5(a). The system, which is part of The University of Michigan's millimeter wave polarimeter (22], uses an HP8510A vector network analyzer to sweep frequency from 34 to 36 GHz in 401 discrete steps. Subsequent to calibration against a metal sphere of known radar cross section, the output is presented in the form of a frequency spectrum of the measured backscattering cross section per unit illuminated area (i.e., a) or in the form of a plot of the received power versus round-trip delay. A more detailed description of the system's operation and signal processing capabilities is given in Ulaby et al. [22]. The scatterometer uses a pair of 15-cm diameter lenscorrected horn antennas mounted onto a common positioner one above the other in the elevation plane. The antenna far-field distance is approximately 5.2 m and the effective beamwidth of the product gain pattern is 3~. The antenna positioner may be set at an angle of incidence 0 from 0~ (nadir) to 90~. and the platform height may be extended up to a maximum of 20 m above the ground surface. Two types of terrain targets were selected: 1) an asphalt surface, as a representative of targets from which the backscatter is due primarily to surface scattering, and 2) a layer of dry snow over a soil surface, as a representative of media from which the backscatter is due primarily to volume scattering. Several experiments were conducted for each of these targets to evaluate the statistical variability of the backscattered power for various combinations of incidence angle and platform height. The measurements for asphalt were acquired with the antennas pointing in the aft direction as shown in Fig. 5(a). To insure that measurements from adjacent footprints were statistically independent, the truck was moved a distance greater than the extent of the antenna footprint between successive measurements. The arrangement for snow was similar to that employed for asphalt except that the truck remained stationary and the boom was made to move in azimuth in order to avoid disturbing the snow surface. The rms height of the asphalt surface was measured to be 0.4 mm (from a surface mold), and the snow was 15 cm deep and had an average temperature of - 1 ~C. To limit the scope of the data-collection segment of this investigation, all observations were made with the VV polarization only. In addition to the scatterometer system, the truckmounted platform earned three microwave radiometers that were mounted on the same platform and their beams pointed along the same direction as that of the scatterometer. Their center frequencies were 35, 94, and 140 GHz. and all three had temperature resolutions better than 1 K. At the time of this investigation, however, only the two upper-frequency radiometers were in operating condition. These instruments proved extremely useful in verifying that I) the targets were uniform, and 2) the snow was dry (i.e., it contained no water in liquid form). Table I provides a summary of the statistics of the radiometric observations. At 94 GHz, the mean value of the brightness temperature TR based on measurements from 10 spatially independent footprints was 252.8 K and the standard deviation was only 1.2 K, which is an excellent indicator that the asphalt surface was electromagnetically uniform. For snow, the radiometric observations were made at both 94 and 140 GHz, and from heights of 11 and 19 m. The two 94-GHz sets of observations (each consisting of 50 measurements from spatially independent footprints) had mean value that were within 1 K of one another and standard deviations of only a few kelvins each. In spite of the slighly greater difference between the mean values of the 140-GHz observations (which is attnbuted to the greater sensitivity of the 140-GHz radiometer (relative to the 94-GHz radiometer) to variations in cloud conditions between the times corresponding to the 11- and 19-m experiments), the results again indicate that the snow medium was fairly uniform from one location to another. The magnitude of TBg 166 K at 94 GHz is characteristic of dry snow [23], and considering that a change in liquid water content by only 2 percent would cause TR to increase by about 100 K [23], the measured standard deviation of only a few kelvins is a clear inciicator that the snow layer was indeed dry everywhere. By way of comparisons, we show in Fig. 10 radiometric observations that were made later in the season for wet snow. We observe that TB of wet snow is about 266 K at 94 GHz (compared to 166 K for dry snow) and 270 K at 140 GHz (compared to about 208 K for dry snow), and again the standard deviations are only on the order of 1-2 K. A. Single-Frequency Observations As was mentioned in the previous section, the scatterometer measures the backscattered power at 401 equally spaced frequencies (channels) extending from 34 to 36 GHz. In this section we shall consider only the statistics associated with single-frequency measurements, namely the 35-GHz channel. It should be noted, however, that the results and conclusions realized at 35 GHz are statistically indistinguishable from those found at lower and higher frequencies in the 34-36 GHz range. Our first example showing the variability of the backscattered power as a function of spatial position was presented earlier in Fig. 5(b) for an asphalt surface, and the associated probability density function was compared to the exponential distribution in Fig. 6. Similar results were obtained for snow and a summary of the observed statustics is given in Table II. The asphalt results given in Table II are divided into two groups: (a) the near-nadir group (0~ and 4~), and (b) the higher-incidence-angle group (20~ and 40~). This division is necessary because the mechanics of signal fading are different in these two angular regions. At incidence angles near normal incidence, the backscattered power consists of a coherent component PC and an incoherent component Pi [3, p. 1812], and only the latter is subject

199 ULABY. HADDOCK et il. FLUCTUATION STATISTICS OF mm-VWAVE SCATFRING 277 Brightness Temperature of Wet Snow 272 2 2 T C,0 'DE 00, 0 c z 8 to signal tading fluctuations. Thus P = PC + P, Sp = Sp, and (45) (46) (47) Sp Sp, P PC + Pt For the incoherent component, Rayleigh fading suggests that Sp, = P, (20). Hence Spatial Position Index Fig. 10. Measured variation of the brightness temperaturc ot wet snow at 94 and 140 GHz. Each data point represents an independent footprint. TABLE I SUMMARY OF VERTICALLY POLARIZED RADIOMNETRIC OBSERVATIONS MADE CONTEMPORANEOUSLY WITH THE RADAR OBSERVATIONS (The incidence angle was 40~. and the statistics are based on observations of N spatially independent footprints. Sr, is the measured standard deviation of To.) TARGET FREQUENCY HEIGHT N TB STB Asphalt 94 GHz 5m I0 252.8 2 Snow 94GHz 11n 50 166.8 35 Snow 94GHz 19 m 50 1658 9 Snow 140 GHz 1 m 50 206.8: 3 Snow 140GHz 19m 50 210.8 2 TABLE 11 SUMMARY OF THE STATISTICS ASSOCIATED WITH THE BACKSCATTERING MEASURESMENTS FROM ASPHALT AND SNOW (INDFPENDENT FOOTPRINTS). - = P,/(P, + i,) P (48) which is always significantly smaller than 1 if P, is significant in magnitude relative to P,. The coherent component Pc is largest at 0 = 0. decreases exponentially with increasing 0, and becomes negligible in comparison with P, (for most natural surfaces) at angles greater than a few degrees [241. Consequently. the value of sp/P computed on the basis of the experimental data was found to be 0.35 at 0 = 0~, 0.59 at 60 = 4~, and close to I at 20~ and 40~. The major conclusions reached on the basis of the single-frequency observations are 1) The Rayleigh model is a reasonable descriptor of signal fading for uniform targets. This is supported by the good agreement shown in Fig. 6 between the measured pdf and the exponential distribution and by the result that sp/P 1I for both asphalt and snow (the deviation from an exact value of 1 is attributed to the fact that the sample size is only 50. and therefore the values of Sp and P given in Table II are merely measured estimates of the true values). 2) No discernable difference between the statistics for the backscatter from snow and those for asphalt is observed. 3) No discernable dependence on footprint size is observed over the range of values examined in this study, which varied in footprint area from 0.07 to 3.24 m2. The corresponding dimensions of the major and minor axes of the elliptically shaped footprint were 0.29 m x 0.29 m for the smallest footprint and 1.75 m x 2.31 m for the largest. B. Frequency Averaging Fig. 11 displays a typical example of the frequency spectrum of the measured power for a given footprint. We observe that P varies relatively slowly as a function of frequency, implying high correlation between adjacent frequency points, but the overall variation across the 34 -36 GHz band is on the order of 23 dB. The improvement (reduction) in spatial variability of the return provided by frequency averaging is demonstrated in Fig. 12, which shows both single-frequency measurements and the 2-GHz averaged measurements of the return from snow as a function of spatial position. The a -4 - T Sp/P ASPHALT Incidence Angle Helghlt(m) a x b (cmxcm) Areadm ) o. 4. 20' 20' 40* 40* 40 29x29 98 71 x 72 0.07 0 35 0.41 0.59 4.1 9.9 4.2 10.3 31 x33 76 x 82 39 x 52 98 x 129 0.08 0.50 0.16 1 01 Area(m ) 0.16 1 08 3 24 3 99 1 31 0 82 1 07 sp'P 1.10 1 21 1 00 Incidence Angle Helght(m) a x b (cmxcm) 40. 40 -40. 4.2 10.7 18.5 40 x 52 101 x 134 175 x 231

200 278 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 26. NO 3. MAY 1988 Measured vs. Theoretical Autocorrelations - Snow -20 -30 0. -40 - -. 34.0 34.5 35.0 35.5 36.0 Frequency (GHz) Fig. 11. Typical trace of the frequency variation from 34 to 36 GHz of the received power for a given footprint of snow. 10 SNOW incidence angie. 40 deqrees H =10.3 meters, 5' i -20 ---- CW(35GHz). S P = 1.0 -30 1 0 10 2 0 0 40 50 Spatia Position Index Fig. 12. Reduction of signal variability through frequency averaging. associated normalized standard deviation is 1.0 for the single-frequency data, compared to 0.27 for the frequency-averaged data. I) Correlation Function: Now we shall examine the role of frequency averaging relative to the theoretical expectations presented in Section III-B. We have 50 traces corresponding to 50 independent footprints, each consisting of measurements at 401 frequencies. Let us denote P, (vj) as the measured power corresponding to spatial position i (with i = 1, 2, *. 50) and frequency v,, where vl = 34 GHz, Vj = [34 + 0.005(j - 1)] GHz. andj = 1, 2, * *, 401. For each position i, we compute the autocovariance 1,vf R,, (A) =.- [P,(v) - P][P,(v+k) - P] (49) Nf = I and the correlation coefficient Pi (A) = R. (AP) (50) R(., (0) where k is the displacement index. Nf = 401 - k, Av = 0.005k (GHz), and P is the mean value of P, (j ) averaged over both i and j. The correlation coefficient is computed for integer values of k from 0 to 200, corresponding to a range of Av from 0 to 1 GHz. Once this process has been completed for each position i. the correlation function pi (Av) is averaged over all i to obtain a better estimate of its frequency spectrum. Thus, the measured cor so c 5 0 ~ < 1.0" 0.8 0 6 - \\ Eqn. 37 0.4 02 / Measured 0.0 - - 0 50 100 150 200 25 Shift(MHz) 0 Fig. 13. Companson of theoretical autocorrelation function given by (37) with that computed on the basis of the spectral measurements ot the radar backscatter. relation coefficient is given by 1 50.( ) Pm,(~v) = - Z,(Al ). 50 i= i (51) A plot of p, (Av) is shown in Fig. 13 for snow. The figure also includes a plot of the expression given by (37). We observe that the measured correlation coefficient decreases with increasing frequency shift Av in an exponential-like manner and at a rate somewhat faster than the theoretical function. Similar results were observed for asphalt. 2) Normalized Standard Deviation: From (37) and (40), the normalized standard deviation associated with the received power P, when averaged over a bandwidth B, is given by (,B) 2 8( \() P i =[ L\ '~ B-)Pfdj/ (52) where f = Av. Fig. 14 shows plots of the normalized standard deviation as given by (52) with p(Av) = p,(Av), the measured correlation function, and with p(Av) as given by (37). The figure also includes a plot of the normalized standard deviation as computed directly from measured data. For a given bandwidth B. Sp is based on the values of P measured at all frequencies between B/2 below and B/2 above 35 GHz. We observe that the 'measured" normalized standard deviation is close to the curve calculated using the theoretical expression for p ( Av) given by (37) and that using the experimental function p, (Av ). To provide a simple formula for estimating sp(B)/P, we propose to use Sp = 1- 1V P NV B ' ll, for B > Avd for B < Avd (53) with AvP selected to provide a good fit to the data. This process led to AvP = 138/D, (in megahertz) (54)

201 UI.AUY. H-DDOCK l, (, i F.LIC 'UAI\i()N STATISTICS (F imnl-WAVE SC.AI lFR.KN(; "'9 1 2 la. 0. 0L 1.( 0.1 0.1 0. O. Snow 0 6 -Theoretical p(Av) 4 Data Measured pm(Tv),2-.,-T l. 500 1000 Bandwidth (MHz) (a) 1500 2000 2 r.0 - 08 -0.6 - 0.4 - Asphalt Theoretical p(av) Data Measured pm(Av) 0.2 '.,, 500 1000 Bandwidth (MHz) (b) 1500 2000 Fig. 14. Normalized standard deviation versus bandwidth B. "Thcorctical' refers to (52) with p( Av) given by (37); "Measured' rclers to 152) with p(Av) = p,,,( v). (a) Snow and (b) asphalt. 1.2 -Snow 1.0. ^ --- Eq. (53) 0.8 i0 -Data,L 0.6 0.4 0.2 0 500 1000 1500 2000 Bandwidth (MHz) Fig. 15. Comparison of measured nornialized standard deviation ith single model given by (53) and (54). and is shown graphically in Fig. 15. In the above expressions, B is in megahertz and D (the slant-range resolution defined in Fig. 9) is in meters. These formulas, which provide an excellent fit to the measured normalized standard deviation for N - 2, indicate that the effective decorrelation bandwidth is approximately equal to the theoretical value predicted by (38). VI. RESULTS OF THE BISTATIC SCATTERING OBSERVATIONS The scatterometer system that was used to acquire the backscattering data reported in the preceding section had been designed to operate in a bistatic mode as well [22]. Bistatic scattering measurements were made for several sand and gravel surfaces using the arrangement shown in Fig. 16. Details of the results and their significance are given elsewhere [25]; our present interest pertains to the variability of the bistatically scattered signal only. Moreover, we shall limit the discussion to a typical example. In one of the bistatic scattering experiments, the received power was measured at many azimuth angles a ranging between 10~ and 180~ for fixed and equal values of the incidence angle Oi and scattered elevation angle 0,, namely 0i = 0, = 66~. The configuration with 4 = 180~ corresponds to the specular case. At each angle 0, the

202 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 26. NO. 3. MAY 1988 Z 20O RIECEIVE X y Fig. 16. Bistatic arrangement. Smooth Sand (HH polarization) I I.1 0, Average of points between azimuth angles of 10~-150 -0 198.. 66' I I I I.. I.,', 0 30 0O 90 120 150 180 Azimuth (deg) Fig. 17. Standard-deviation to mean ratio versus the azimuth angle tor the bistatic scattering measurements for smooth sand. Each data point is based on measurements of 10 spatial positions (tootprints). received power was measured for 10 independent (nonoverlapping) footprints on the target surface. This was accomplished by partially rotating the target platform between successive measurements. In each case the recorded power was the received power averaged over the 34-36 GHz band. Fig. 17 presents a plot of sp/P as a function of 0. At each value of 0, Sp, and P were calculated using the 10 observations described above. Except for the near-specular directions ( -= 180~), the normalized standard deviation exhibited an approximately constant value of 0.2. The equivalent total number of independent samples is N = (1/0.2)2 = 25. Thus, frequency averaging provides about 2.5 independent samples per spatial sample. The power scattered in the specular direction was dominated by coherent scattering that is not subject to fading. Consequently, the measured normalized standard deviation was found to be only 0.015. This result is analogous with the backscattering result for the normal incidence case (see (48) and Table II). VII. CONCLUSIONS This paper has shown that the Rayleigh fading model is indeed appropriate for characterizing the fluctuation statistics associated with radar scattering from terrain, provided the target's properties satisfy the model's underlying assumptions. One of these properties is spatial uniformity. If the terrain target is an asphalt surface. snow-covered ground, or a grass surface, the Rayleigh model gives results in good agreement with experimental observations at 35 GHz for both backscattering and bistatic scattering. When frequency averaging is used to reduce the variability of the radar return, however, the formulations based on the Rayleigh model provide a good estimate of the improvement provided by frequency averaging. REFERENCES Il1 Webster's New Collegiate Dictionary. Springfield. MA: G. and C. Merriam Pub. Co., 1975, p. 411. [2] R. L. Freeman, Telecommunication Transmission Handbook. New York: Wiley, 1975. pp. 197-198. 131 F. T. Ulaby, R. K. Moore, and A. K. Fung, Microwave Remote Sensing: Active and Passive. Vols. I-Ill. Reading, MA: Addison-Wesley/Artech House, 1981-1987. 141 F. T. Ulaby. F. Kouyate, B. Brisco, and T. H. L. Williams, "Textural information in SAR images," IEEE Trans. Geosci. Remote Sensing, vol. GE-24. no 2. pp. 235-245, 1986. 51] J. S. Zelenka. "Comparison of continuous and discrete mixed-intr

203 ULABY. HADDOCK et ul.: FLUCTUATION STATISTICS OF mm-WAVE SCATTERING 281 grator processors.' J. Opt. Soc. Amer.. vol. 66. pp. 1295-1304. 1976. 161 L. J. Porcello. N. G. Massey. R. B. Innes. and J. M Marks, "Speckle reduction in synthetic aperture radar." J. Opt. Soc. Amer.. vol. 66, pp. 1305-131 1. 1976. 171 L. J. Cutrona. E. N. Lcith. C. J. Palermo. and L. J. Porcell. "Optical data processing and filtering systems," IEEE Trans. Inform. Theory. vol. 2. no. 6. pp. 386-400. 181 F. K. Li. C. Croft. and D. N. Held, "Comparison of several techniques to obtain multiple-look SAR imagery, " IEEE Trans. Geosci. Remote Sensing,. vol. GE-21, pp. 370-375. 191 C. Y. Chi. D. G. Long, and F. K. Li. 'Radar backscatter measurement accuracies using digital Doppler processors in spaceborne scatterometers," IEEE Trans. Geosci. Remote Sensing, vol. GE-24, pp. 426-437, 1986. 1101 P. Beckmann. Probability in Communication Engineering. New York: Harcourt, Brace, and World. 1967, p. 122. 1111 W. W. Weinstock, "Radar cross-section target models." and "Illustrative problems in radar detection analysis." in Modern Radar Analysis. Evaluation, and System Design. R. S. Berkowitz. Ed. New York: Wiley, 1965. 1121 T. F. Bush and F. T. Ulaby, "Fading characteristics of panchromatic radar backscatter from selected agricultural targets," IEEE Traits. Geosci. Electron.. vol. GE-13. no. 4. pp. 149-157. Oct. 1975. [131 G. R. Valenzuela and M. B. Laing, "Point-scatterer formulatiorn of terrain clutter statistics," Naval Res. Lab. Rep. 7459. Washington. DC, Sept. 27. 1972. 1141 M. P. Warden, '"An experimental study of some clutter characteristics," in AGARD Cotf. Proc. No. 66 Advanced Radar Systems, Nov. 1970. 1151 R. R. Boothe. "The Weibull distribution applied to the ground clutter hackscatter coefficient." in Automatic Detection and Radar Data Processing. D. C. Schleher. Ed. Dedham. MA: Anech House, 1980. 1161 D. C. Schleher. "Radar detection in log normal clutter." presented at the IEEE Int. Radar Conf.. Washington, DC. 1975. 1171 H. Kashihara. K. Nakada, M. Murata. M. Hiroguchi. and H. Abia. 'A study of amplitude distribution of space-borne synthetic aperture radar(SAR) data." in Proc. ISNCR Symp. (Tokyo), Oct. 22-24. 1984. 1181 J. K. Jao, "Amplitude distribution of composite terrain radar clutter and the K-distribution," IEEE Trans. Antennas Propagat., vol. AP32, pp. 1049-1062, 1984. 1191 D. E. Kerr, Ed., Propagation of Short Radio Waves. New York: McGraw-Hill. 1951, p. 554. 1201 S. 0. Rice, "Mathematical analysis of random noise." Bell. Svst. Tech. J., vol. 23, p. 282, 1944. 1211 W. P. Birkemeier and N. D. Wallace, "Radar tracking accuracy improvement by means of pulse-to-pulse frequency modulation." IEEE Trans. Commun. Electron., pp. 571-575. Jan. 1963. 1221 F. T. Ulaby. T. F. Haddock. J. East. and M. Whitt. "A millimeterwave network analyzer based scatterometer," IEEE Trans. Geosct. Remote Sensing, vol. GE-26. Jan. 1988. 123] R. Hofer and W. Good, "Snow parameter determination by multichannel microwave radiometry," Remote Sensing Environ., vol. 8, pp. 211-224, 1979. (24] F. T. Ulaby. C. T. Allen. and A. K. Fung. "Method for retrieving the true backscattering coefficient from measurements with real antenna." IEEE Trans. Geosci. Remote Sensing. vol. GE-21, no. 3, pp. 308-313. July 1983. 1251 F. T. Ulaby, T. E. van Deventer. T. F. Haddock. M. Coluzzi. and J. R. East. "Millimeter-wave bistatic scattering from ground and vegetation targets." IEEE Trans. Geosci. Remote Sensing, vol. GE-26, May 1988. * Fawwaz T. Ulaby (M'68-SM'74-F'80) was born in Damascus. Syria. on February 4. 1943. He reccived the B.S. degree in physics from the American University of Beirut. Lebanon. in 1964 and the M.S.E.E. and Ph.D. degrees in electrical engineering from the University of Texas. Austin, in 1966 and 1968, respectively.; From 1968 to 1984, he was with the Electrical Engineering Department at the University of Kans' as, Lawrence, where he was the J. L. Constant Distinguished Professor, and the University of Kansas Center for Research, where he was Director of the Remote Scnsing Laboratory. He is currently with the Radiation Laboratory and the Department of Electrical and Computer Engineering. University of Michigan, Ann Arbor. His current research interests involve microwave propagation and active and passive microwave remote sensing. Along with R. K. Moore and A. K. Fung. he is a coauthor of the three-volume series Microatrve Remote Sensing: Actire anid Passive (Reading, MA: Addison-Wesley). In addition, he is coeditor of the Manual of Remore Sensing. 2nd ed.. vol. 1. American Society of Photogrammetry. Dr. Ulaby is a member of Eta Kappa Nu. Tau Beta Pi, and Sigma Xi. He has been named the Executive Editor tor IEEE TRANSACTIONS ON GFoSCIENCE AND RF.M)TF. SENSING. 1984-1985. and was the Geoscience and Remote Sensing Society's Distinguished Lecturer for 1987. He was named an IEEE Fellow in 1980 "for contributions to the application of radar to remote sensing for agriculture and hydrology, received the GRS Society's Outstanding Service Award in 1982. and its Distinguished Service Award in 1983. In 1984, he also received a Presidential Citation for Meritorious service from the American Service of Photograinmetry. He received the University of Kansas Chancellor's Award for Excellence in Teaching in 1980, the University of Kansas Gould Award for "distinguished service to higher education" in 1973. and the Eta Kappa Nu MacDonald Award as an "outstanding electrical engincering professor in the United States ol.Anerica" in 1975. Thomas F. Haddock (M'86) was born in Washington. DC. on November 2. 1949. He received the B.A. degree in mathematics and the M.S. and Ph.D. degrees in physics from the University of 1Michigan, Ann Arbor. in 1972. 1977, and 1984. respectively. Fron 1984 to 1985 he was Manager of Development Projects at Applied Intelligent Systems, a \machine vision firm involved in real-time optical.,./. infrared. and X-ray vision systems. He is currently with the Radiation Laboratory and the Department of Electrical Engineering and Computer Science, University of Michigan. He has conducted research in the fast flux density variations of quasi-stellar objects at a wavelength of 12.5 nmm. Other research has included development of real-time alphanumeric character recognition algorithms and ultrasonic weld inspection algorithms. Prior to receiving the Ph.D. degree, he worked as Applications Engineer for Sams/3M. a manufacturer of hean-lung machines and cardiac assist devices, where he developed electrodes for manufacturing applications. His current research interests are millimeter-wave scattering and emission from natural targets. Dr. Haddock is a member of the American Astronomical Society. - I Richard T. Austin (S'84) was born in Maryville. TN, on September 2. 1964. He received the B.S. degree in electrical engineering front the University of Kentucky. Lexington. in 1986 and the M.S. degree in electrical engineering from the University of Michigan, Ann Arbor, in 1987. He is currently working toward the Ph.D. degree in the Radiation Laboratory at the University of Michigan. He holds a National Science Foundation Graduate Fellowship. His research interests include millimeter-wave radiomctry, millimeter-wave radar, and scattering and emission from natural targets.

204 APPENDIX E MILLIMETER-WAVE RADAR SCATTERING FROM SNOW

205 Radio Science, Volume 26, Number 2, Pages 343-351, March-April 1991 Millimeter-wave radar scattering from snow: 2. Comparison of theory with experimental observations Fawwaz T. Ulaby, Thomas F. Haddock, Richard T. Austin, and Yasuo Kuga Radiation Laboratory, Department of Electrical Engineering and Computer Science. University of Michigan, Ann Arbor (Received May 17, 1990; revised October 30, 1990; accepted November 6, 1990.) Using a truck-mounted platform, backscatter measurements were made at 35, 95, and 140 GHz for a variety of snow conditions to evaluate the radar response to incidence angle, surface roughness, and liquid water content. Good agreement was obtained between the experimental observations and theoretical calculations based on the numerical solution of the radiative transfer equation presented in the preceding paper. A notable exception is when the snowpack is in the refreezing phase of the diurnal cycle, during which the snowpack is characterized by a dry surface boundary with wet layers underneath. To accommodate this type of condition, a hybrid first-order numerical solution is proposed. The hybrid approach provides excellent agreement between theory and experiment. 1. INTRODUCTION In Part 1 [Kuga et al., this issue] of this two-part sequence we proposed a radiative transfer model for characterizing radar scattering from snow at millimeter wavelengths. The purpose of the present paper is to (1) provide a summary of the observed behavior of the radar backscattering coefficient ar~ at 35, 95, and 140 GHz, including its dependence on incidence angle, surface roughness, and liquid water content, (2) compare the measured data to theoretical calculations, where possible, and (3) propose a hybrid first-order numerical model for explaining the diurnal variation of r~. The measurements reported in this paper were acquired by the University of Michigan's truckmounted millimeter-wave scatterometer [Ulaby et al., 1988a; Haddock and Ulaby, 1990]. A summary of the system's specifications is given below. frequencies 35, 94, 140 GHz; IF bandwidth 0 to 2.0 GHz; transmit power 35 GHz, +3 dBm; 94 GHz, 0 dBm; 140 GHz, -4 dBm; sweep rate I m-s/frequency, 51, 101, 201, 401; noise eq crossp pha near-fit uivalent o~ 35 GHz, -22 dB; 94 GHz, -28 dB; 140 GHz, -21 dB; )ol isolation 35 GHz, 23 dB; 94 GHz, 20 dB; 140 GHz, 10 dB; ase stability 35 GHz, -1 deg/hour; 94 GHz, -1 deg/min; 140 GHz, -10 to 50 deg/s; eld distance 35 GHz, 2.7 m; 94 GHz, 7.3 m; 140 GHz, 2.7 m; beamwidth 35 GHz, R: 4.2' T: 4.2~; 94 GHz, R: 1.4' T: 2.8~; 140 GHz, R: 2.2~ T: 11.8~; na diameter 35 GHz, R: 6 in T: 6 in; 94 GHz, R: 6 in T: 3 in; 140 GHz, R: 3 in T: 0.36 in; I processing HP 8510A/8511A based; anten signal output products received power versus range; received power versus frequency (at fixed R); phase and amplitude for each frequency. polarization incidence angles platform height frequencies/sweep; HH, HV, VV, VH; 0~ to 70~; 3 m minimum, to 18 m maximum; Ground-truth observations were made for the following parameters: (I) air temperature, (2) nearsurface snow temperature, as well as the temperature at deeper locations in the snow layer, (3) snow density p, (g/cm3), (4) height profile of the snow surface using a graded metal plate inserted edgewise into the snow, from which the rms height s is calculated, (5) depth of the snow layer, (6) volumetric liquid water content of the surface 5-cm layer, m,(%), Copyright 1991 by the American Geophysical Union. Paper number 90RS02559. 0048-6604/91/90RS-02559S08.00 343

206 344 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW. 2 measured using the freezing calorimeter technique, and (7) microscope photographs of thin snow samples, from which the ice crystal size is estimated. The standard measurement procedure consisted of measuring the backscattered power for hh, hv, vh, and vv polarizations at each of 30 or more spatial locations. The backscattering coefficient was obtained by averaging the measurements for the different spatial locations. In addition to spatial averaging, frequency averaging was used to further improve measurement precision [Ulaby et al., 1988b]. The estimated uncertainty associated with the values reported in this paper is +0.5 dB. Analysis of the data shows that the hv and vh measurements are essentially identical (within a fraction of I dB), which is expected from the reciprocity relation. In almost all cases the copolarized responses. crhh and,,, were within 1-2 dB of each other. 2. ANGULAR RESPONSE The data shown in Figure I were measured for a 12-cm-thick layer of dry, freshly fallen, unmetamorphosed snow composed of ice crystals with diameters on the order of 1 mm. The measured rms height was 1.4 mm. Only hv- and vv-polarized data are shown because the difference between ahh and or is I dB or less across the entire angular range at all three frequencies. The curves shown in the figure were calculated according to the theoretical model described in the preceding paper. For vv polarization, theory and experiment are in good agreement at 35 and 95 GHz. At 140 GHz, however, the level predicted by theory is lower than the experimental observations for vv polarization by about 4 dB. We attribute the difference to the backscattering enhancement effect, which the model does not take into account. A typical example of the angular dependence for wet snow is shown in Figure 2. The theoretical curves were computed assuming a mean ice crystal diameter of I mm and a rms slope of 0.07. Reasonable agreement between theory and experiment is obtained except for vv polarization of 140 GHz; we again attribute the difference to the backscattering enhancement effect, although no strong evidence exists to support this contention. 3. EFFECT OF SNOW SURFACE ROUGHNESS According to the model results presented in the previous paper, snow surface roughness should 10. 0. -- VVTherI -- vv........ Theo f = 35GHz...,...,..,... i.. I... t...,0. -20. 0. 10. 20. 30. 40. 50. 60. Incident Angle 0, ( degrees ) 10...... --..... i =3 - -= 1: u co 'a ic 23 0. -10. -20.......... HV Th r X VV DuA a f9.GN f = 95 GHz... I.-......-..... 0. 10. 20. 30. 40..50. 60. 70. 0 0 'E,.i 8 u 11.2 -4 I Incident Angle 9, ( degrees ) 10. I --- —- I ---I- l ------ ' --,.........-1.. O...W... - --- Iv'' VI~ -10. -....... rV IIby * VV Dwa * KVDM f= 140 GHz 0. 10. 20. 30. 40. 50. 60. 70. Incident AngIe e ( degrees ) Fig. 1. Measured and calculated backscattering coefficient for a dry snowpack with the following parameters: depth, 12 cm; snow density, 0.2 g/cm3; mean crystal diameter, I mm; and rms surface slope, 0.07. have a minor effect on the level of Tr~ (except at normal incidence) when the snow is dry. For wet snow, however, o-~ should increase by as much as 5; dB at 35 GHz if the surface is made rough relative to the wavelength, but the increase should not be significant at 95 GHz or higher frequencies. There is

207 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 345 10. m 0 O Z5.2 u E 8 II o4 M0.a se C" m 7) ap 0 E ao C: u co F9 Q u r.r. r I; U C's m. 0. 10. 20. 30. 40. 50. 60. n 10. r 0. *...... T *** I *.*.*.* T.... *. *.... 95 GHz ---- WThory mv=5%.-......HV Theno * VV -m ~* W Da * HVD ma * -a... 0. 10. 20. 30. 40. 50. 60. -10. t -20. I 10. 0. 140 GHz _. w-V r n=3%........- HYVlo VV WDNa \ * HVYDV - r.~ —....,.................................... *.. snow at 95 GHz. The data includes measurements for a 25-cm-deep snowpack made before and after artificially roughening the surface. The measured rms height was found to be 0.7 cm for the undisturbed surface and 2.0 cm for the roughened surface. Considering that A = 0.3 cm at 95 GHz, both surface conditions are electromagnetically rough, so it is not surprising that the angular curves shown in Figure 3 for the two surface conditions are within I dB of each other for both vv and hv polarizations. A more detailed examination of the effect of surface roughness on the radar response from snow was conducted by observing the radar backscatter as a function of time for three sections of a 71-cmdeep snowpack with different surface roughnesses. Twelve sets of observations were made, all at 00 = 40~, for each surface roughness, starting at 0630 LT and ending at 2200 LT; each set consisted of measurements at 35 and 95 GHz for all linear polarization combinations. The measured rms heights of the three surfaces were s = 0.49 cm, s2 = 0.89 cm, and 53 = 1.98 cm, and the liquid water content of the 5-cm surface layer exhibited a Gaussianlike variation with a peak value of 4.8% (Figure 4). Figure 4 shows a sample of the measured data, specifically the cross-polarized diurnal responses of all three roughnesses at both 35 and 95 GHz. According to the 35 GHz data, increasing the rms height from 0.49 to 0.88 cm (which corresponds to increasing s/A from 0.57 (slightly rough) to 1.02 (rough)) causes a~ to increase by 1-3 dB, but increasing the roughness further to s3 = 1.98 cm does not seem to have much of an impact on or~. At 95 GHz, even the least rough surface (with sl/A = 1.5) is electromagnetically very rough. Hence cr~ exhibits approximately the same diurnal pattern for all three surface roughnesses. A similar behavior was observed for vv polarization. 4. DIURNAL RESPONSE The plot shown in Figure 5 depicts the temporal variation of the volumetric liquid water content of the top 5-cm snow layer over a period of 14 hours, starting at 0800 LT. The snowpack was dry until 1100 LT, then m, increased rapidly to a peak value of 7% at 1400 LT, and then returned to the dry condition by 1900 LT. This description applies to only the top 5-cm layer, treated as a single layer; it does not provide any information on lower layers or on the depth profile of m, within the 5-cm layer. -10. -20. 0. 10. 20. 30. 40. 50. 60. 70. Incident Angle O (degrees ) Fig. 2. Measured and calculated backscattering coefficient for a wet snowpack with the following parameters: depth, 27 cm; snow density, 0.4 g/cm3; mean crystal diameter, 1 mm; and rms surface slope, 0.07. plenty of experimental evidence to support the model expectations with regard to dry snow, both at 35 GHz [Stiles and Ulaby, 1980] and 95 GHz [Williams et al., 1988]. Additional support is provided by the data shown in Figure 3 which shows plots of oa~ versus the incidence angle 00 for dry

208 346 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 6 - 7_ hv Polanzanon 4 4 - J -2 - 0 * - Moderately Rough Surface (s = 0.7cm) * - Very Rough Surface (s= 2cm) m, =0 -4 Density =0.35 g/cnm] 95 GHz 0 i b Li 0 10 20 30 40 50 50,0 Incidence Angle (degrees) Fig. 3. Measured backscattenng coefficient of a dry snowpack at 95 GHz under two different surface roughness conditions. The significance of this statement will become apparent later when we compare measured values of m, with those predicted by theory. Ideally, it would be desirable to measure the entire depth profile of m, with a vertical resolution on the order of millimeters. In practice, however, it is very difficult to make accurate measurements of m, for samples thinner than 5 cm. Furthermore, because it takes about 30 min to process each snow sample (using the freezing calorimeter technique) and because m, varies rapidly with time, it is impractical to sample more than one depth layer, unless additional trained manpower is used (two people, working full time, were required to generate the data in Figure 5). The radar response to the observed variation in liquid water content is shown in Figure 6 for 35, 95, and 140 GHz, measured at an incidence angle of 40~. We observe that: 1. The shapes of the diurnal patterns of a7~ are qualitatively mirror images of the my pattern for all polarizations and frequencies. 2. The dynamic range of the variation is greater at 35 GHz and smallest at 140 GHz. 3. The dynamic range of the variation is slightly greater for cross polarization than for like polarization. All three observations are consistent with the predictions of the theoretical model given in the preceding paper. However, as we will see later, quantitative agreement between theory and experiment for the melting part of the cycle is quite different from the agreement for the freezing part of the cycle. To examine this question, we computed o~ as a function of m, using the numerical solution outlined in the preceding paper for the snowpack parameters measured in the field (depth, density, crystal size, etc.). It was assumed that the value of m, measured for the top 5-cm layer is valid for the entire depth of the snow layer. The results are shown in Figure 7 for hh polarization at 00 = 40~. (Also shown in Figure 7 is the computed response of r~ based on the first-order solution of the radiative transfer equation, which will be discussed later in section 4). Allowing for a possible level shift between data and theory, we used the measured ar( data shown in Figure 6 (relative to cr0 of dry snow) to compute mv according to the numerical solution curves given in Figure 7. The plots shown in Figure 8 depict the diurnal variation of mv as measured in the field and as predicted by the o-~ data and the model calculations at 35, 95, and 140 GHz. We observe good overall agreement between the mea.sured values of ma, and those predicted by the 35 -and 95-GHz data for the period between 0800 and 1300 LT, and although the values predicted on the basis of the 140-GHz data are somewhat lower in level than those predicted by the lower-frequency t

209 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 347 Zs E w 0 C: ~o u - 00 u m Co 0t 5. 4. 3. 2. 1. 0. 0. -2. -4. -6. -8. -10. -12. Snow Wetness 6. 8. 10. 12. 14. 16. 18. 20. 2: -- Smooh — 4. —. Siituhy Rough A A 0 —. — v u o — o --- ~ ~ ~ ~ - 0 * cYo~ -/. - 35GHz HV X. 6. 8. 10. 12. 14. 16. 18. 20. 2 10. -, ---,-) ---,-.-,-, -, —,-I-..- Liquid Water Content (%) - Snow Tan'x -v./, -......... Air Tempeaure ------ iquid WarCom m, - ' I I I - I -. i - I 2. -10. 2. I.., I I I I I I. 0.- Aa -- I AV.~.2. - A - ' I95GHz VH 4 -4. — I — I-I — I 6. 8. 10. 12. 14. 16. 18. 20. 22 Tune of Day (hour) Fig. 4. Observed diurnal variation of m.,, the liquid water content of the top 5-cm layer, and o~ at 35 and 95 GHz for cross polarization. The snowpack was 28 cm deep, had a mean crystal diameter of 0.6 mm, and a density of 0.4 g/cm3. The measured rms heights of the surfaces were: sl = 0.49 cm (smooth), s2 = 0.89 cm (slightly rough), and s3 = 1.98 cm (very rough). channels, the temporal variation is similar to that exhibited by the other plots. A totally different situation is observed in Figure 8 for the period from 1300 to 1900 LT, during which the snowpack undergoes a refreezing process. In fact, when the measured mrn is at its peak value of 7% (at approximately 1400 LT), all three radarpredicted values of mv are in the 0.8-1.5% range. Clearly, the models, as used thus far, are inappropriate for modeling the backscatter from wet snow under these "refreezing" conditions. This lack of correspondence between theory and measurement is caused by the fundamental assumption that mv 8. 10. 12. 14. 16. 18. 20. 22. 24. Time of Day ( hour ) Fig. 5. Measured air temperature, snow surface temperature, and liquid water content of the top 5-cm layer on March 1, 1990. has a uniform depth profile between the snow and ground surfaces. During the melting phase of the diurnal cycle the melting process starts at the airsnow interface because the source of thermal energy is the warmer air mass above the snow pack or direct solar radiation (under cloud-free conditions). For the experiment under discussion the air temperature rose from -8~C at 0800 LT to a high of +6~C at 1145 LT, and then decreased down to -3~C at 2000 LT. As the snow surface layer starts to melt, its thermal conductivity increases, thereby allowing the transfer of thermal energy down to lower layers. During the melting phase, mv is highest at the surface and decreases in an exponentiallike manner with depth. At millimeter wavelengths the penetration depth in dry snow is of the order of 1 m at 35 GHz, 5 cm at 95 GHz, and 1 cm at 140 GHz; and for mv = 3%, the penetration depth is only 3 cm at 35 GHz and less than I cm at the higher frequencies. Hence, when the top snow layer is wet, it alone governs the radar backscatter, and the lower layers exercise an inconsequential influence on o~. A different situation occurs during the refreezing phase. As the air temperature starts to decrease, we have a reversal in energy balance with the snowpack becoming the source of excess energy. Freezing starts at the surface and proceeds slowly downward to lower layers. This results in an inverted profile with m, = 0 at the surface, increasing slowly to a maximum value at some depth below the surface, and then decreasing again as the depth is

'210 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW. 2 348 a 8 -.I 22 0. -10. -20. 0., a 10. 12. 14. 16. 18. 20. 22. 24. 0. -VV Pou1.a 95 GHz.....V or* ----- Poimma U U 0 8 U.5 22 0. - 10. - 20. -30. 0. 9 5 G Hz Fm Oirde - 95 0Hz NuozncA1 -10. 10. S. 10. 12. 14. 16. 1. 20. 22. 24... I I O. -. -....... VN Pod 14G -— H - im 140 GHz..... HH P,,,,, -20. 0. -10. -20. 0. 2. 4. 6. 8. 10. 1; - 140 GHz Fir Order, 140GHz Numecal....................... I -10. A I! 0. 2. 4. 6. 8. 10. 12. -20. -1 - 8. 10. 112. 14. 16. 18. 20. 22. 24. Tune of Day ( hour ) Fig. 6. Backscattering variation measured at 35, 95, and 140 GHz on March 1, 1990. The incidence angle was 40~. The snowpack density was 0.32 g/cm3, depth was 12 cm, and other parameters are given in Figure 5. increased further. Figure 9 shows a family of profiles generated on the basis of this simple logic, constrained by the condition that the value of m,, averaged over the 5-cm layer has to be equal to the value measured by the freezing calorimeter. With the uppermost surface layer being dry in the inverted-profile case, the radar response will be governed not only by the dry surface layer, but by the wet layer immediately underneath it as well. It is not surprising, therefore, that the model-predicted val Volumetric Liquid Water Content me ( % ) Fig. 7. Calculated backscattering coefficient for a deep snowpack versus liquid water content for hh polarization at an incidence angle of 40~. The calculations are based on the firstorder and the numerical solutions of the radiative transfer equation for a snowpack with the following properties: mean crystal diameter. 1 mm; and snow density, 0.32 g/cm3. ues of m,, are quite different from the 5-cm average measured in the field. To solve this problem, we propose to use a hybrid first-order numerical model for characterizing the radar response of snow. 4.1. Hybrid first-order-numerical model Let us consider the advantages and limitations of the first-order solution and the numerical solution of the radiative transfer equation as we apply them to the snow problem. The first-order solution takes

211 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 349 Estimations of Liquid Water E a A2 c u To.3....U 0r 8. 7. 6. 5. 4. 3. 2. 1. 0. 8. 10. 12. 14. 16. 18. Time of Day ( hour) Fig. 8. Diurnal variation of liquid water c( labeled "measured" refers to the value of m, layer as measured by the freezing calorimeter. plots refer to m, as predicted by the model observations at the noted frequencies (with hh 00 = 40~). into account only single scattering volume, whereas the numerical soluti account all orders of multiple scatte quently, the first-order solution may u the magnitude of a~, particularly whe ing albedo is large. The scattering albe with increasing liquid water content b increased absorption by the backgrc ("wet" air) and because the ratio of Content refraction of the ice crystal to that of the back-,, ground decreases with increasing snow wetness. To - compare the first-order solution with the numerical 3GHz solution, we refer to the plots in Figure 7 which 95GHz show o0hh versus my for an incidence angle of 40~. I40GH For dry snow (m = 0) the first-order solution underestimates 0o~ by about 3 dB at 35 GHz, 8 dB at 95 GHz, and 11 dB at 140 GHz. The differences between the two solutions decrease with increasing mv, becoming insignificant at 35 GHz for m,, - 2%, but not so at the higher frequencies. Another major advantage of the numerical solution is that it pro-, 4 vides a result for the cross-polarized scattering coefficient o-~, whereas the first-order solution does not predict cross polarization because the scattering ontent. The plot particles are assumed to be spherical in shape. of the top 5-cm Having stated the advantages of the numerical The other three solution over the first-order solution, let us now using the radar consider the converse. The numerical solution repolarization and quires a great deal of computation, in spite of the fact that the snow medium is assumed to have uniform properties throughout. If we are to formuin the snow late the solution for a multilayer structure in order ion takes into to accommodate a nonuniform depth profile for mv, enrng. Conse- the complexity of the numerical approach would nderestimate make the solution computationally impractical. On:n the scattern the scatter- the other hand, the first-order solution is perfectly:do decreases do decreases amenable to computing the backscatter from a iecause of the o t medium with nonuniform properties in the depth )und material nthe mndex f dimension. Hence, we propose to use a hybrid model that takes advantage of the more accurate feature of the numerical solution and the easier structure of the first-order solution. The procedure -,.. is as follows. - 1:30 For a snowpack with specified density and crystal 00 size distribution, we compute ar~ for the semi0 infinitely deep case (at the angle of incidence, - 16:30 - frequency, and polarization of interest) as a func1.7.:30 tion of m, (assuming a uniform profile with depth) ---- 18:15 using the numerical solution outlined in the preced19:15 ing paper. We shall denote this backscattering co20.30 - efficient acr(0o, f, pq), where p and q denote the -- - polarization configuration of the receive and transmit antennas, respectively. Figure 10 shows the.'-... —._ results of such a calculation at 35, 95, and 140 GHz 7. 80 9. 10 for 0o = 40~ and pq = hh. The form of the variation can be fitted to the functional form: E C 0 u C!q m:, Moisture Profile 10. 1 I I I 9. 8. -. 7. - '. 6. - 7., 4. 3..,,7 2.. -.. o. 1. 2. 3. 4. 5. 6. Depth ( cm) Fig. 9. Proposed depth profiles for m,. The pack is assumed to be totally dry both prior to 1100 LT and after 2200 LT. 10 log [o'~(m,)] = ao + alm, 2 + a2m, (1)

212 350 p o ~0O E._ 8 c to r u 10. ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW. 2 tional techniques is in r1, (or its "equivalent" in the HH Polarization case of the numerical solution). We propose to use e = 400 the results of the numerical technique to obtain a more exact estimate of 7q, and then use it in (2) to. ', compute r~ for any my(z) profile. We do this by evaluating (2) for a semi-infinite medium with uni'"" - form liquid water content m,,: i,,(m) cos 00 95- 3GH1 oAm) T,(00o m,)Tq(0o m) -, (3) *-e — 95GHz - 2Ke(mt,) —.,- 140 GiH _, ________ ___ and then equating the result to (1). This process o. 2. 4. 6. 8. to. 12. leads to -30. Liquid Water Content nm ( % ) Fig. 10. Numencal fits using the form given by (1). where a0, al, and a2 are constants (for a given set of radar parameters (80, f, pq)). For a snow medium with wetness profile m,(z) observed at an incidence angle 80, the first-order solution of the radiative transfer equation leads to the following expression for the pq-polarized backscattering coefficient. al = Tp[0o, m,(0)]Tq[o, mv(0)] x, exp -2 z K sec O' dz dz' (2) o Jo ] where Tp and Tq are the transmissivities of the air-snow boundary for backscattered polarization p and incident polarization q, mu(0) is the snow liquid-water-content at z = 0 (air-snow surface), d is the depth of the snow layer, O' is the refraction angle in the snow medium, K' is the extinction coefficient, and 7r, is the volume backscattering coefficient in m2/m3. Both K' and i,, are functions of z if mW(z) is nonuniform. The form of (2) is very convenient for computing the backscattering coefficient when m,(z) is not uniform with depth, but as was mentioned earlier, the first-order solution underestimates the level of a~ for like polarization, and it does not predict a cross-polarized component because 7T, = 0 if p # q for spherical particles. The transmissivity coefficient TpTq is identical for both the numerical and first-order solution, and the same approach used for computing the extinction coefficient Ke can be used in both cases. The key difference between the two computa 2a,(m,.)K;(m,,) sec 80 TV (m,,) = Tp(0o, m,,)Tq(0o, rm.) 2K;(m,,) sec 60 [ao+almn+am, = m T ----antilo [ 1 Tp(o, m,)Tq(0o, m)antilog 10 [aLamr+im, (4) (5) Thus 77,(m), as given by the above expression, represents the "effective" backscattering coefficient per unit volume of the snow medium at the specific frequency f and polarization combination pq used in computing c~. To summarize, the proposed procedure consists of the following steps: 1. For a given set of snow parameters (density and particle size distribution), the QCA method (adjusted with respect to experimental data as discussed in section 2.3 of the preceding paper) is used to compute K.(m,), and the numerical solution is used to compute cr(mv) versus mn, (for a medium with uniform wetness profile) at the specified wave parameters (0, f, pq). 2. For convenience, cr(mv) is fitted to a function of the form given by (1). 3. Equation (5) is used to compute 7rv(mv). 4. Equation (2) is then used to compute oa~ for any specified profile mu(z). 4.2. Results This procedure was applied to the wetness profiles shown in Figure 9, and the results are shown in Figure 11 for hh polarization at 35, 95, and 140 GHz. The plots compare the measured diurnal variation of cr~ (relative to its value for dry-snow conditions) with the variation computed using the hybrid solution outlined above. Very good overall agreement is observed at all three frequencies, and

213 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 351 u 8 13 o i-= y u 04 3 em cs z l u -0.O CE 10. O. O0. 5. 0. -5. 10. -15. -20. 10. 5. 0. -5. - 10. -15. similar agreement was obtained for the other polarizations and for other diurnal experiments. 5. CONCLUSIONS The radiative transfer model with the quasi-crystalline approximation presented in the preceding paper, and the hybrid first-order numerical solution proposed in this paper, together provide an excellent tool for examining the radar response of snow at millimeter wavelengths. This conclusion is supported by comparisons of theoretical predictions with experimental observations made at 35, 95, and 140 GHz. Acknowiedwment. This work was supported by U.S. Army Research otfice contract DAAL03-89-K-0056. REFERENCES Haddock, T. F., and F. T. Ulaby, 140-GHz scatterometer system and measurements of terrain, IEEE Trans. Geosci. Remote Sens., 28, 492-499, 1990. Kuga, Y., F. T. Ulaby. T. F. Haddock, and R. DeRoo, Millimeter-wave radar scattering from snow, 1, Radiative transfer model, Radio Sci.. this issue. Stiles, W. H., and F. T. Ulaby, Microwave remote sensing of snowpacks, Rep. 3263, contract NAS5-23777, NASA, Goddard Space Flight Cent., Greenbelt, Md., 1980. Ulaby, F. T., T. F. Haddock, J. R. East. and M. W. Whitt, A millimeterwave network analyzer based scatterometer, IEEE Trans. Geosci. Remote Sens., 26, 75-81, 1988a. Ulaby, F. T., T. F. Haddock, and R. T. Austin, Fluctuation statistics of millimeter-wave scattering from distributed targets, IEEE Trans. Geosci. Remote Sens., 26, 268-281. 1988b. Williams, L. D., J. G. Gallagher, D. E. Sugden, and R. V. Birnie, Surface properties effects on millimeter-wave backscatter, IEEE Trans. Geosci. Remote Sens., 26, 300-306, 1988. R. T. Austin, T. F. Haddock, Y. Kuga, and F. T. Ulaby, Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109. -20. L8. io_. 10. 12. 14. 16. 18. 20. 22. 24. t 2. 4. 16...... -.! 5. 0. -5.. I I I, I, I i 140 GHz * t. - D~ I.. * D -.o. -15. -20. 8. 10. 12. 14. 16. 18. 20. 22. 24. T'ue of Day ( hour ) Fig. 11. Comparison of the diurnal variation of the measured backscattering coefficient (relative to dry snow) with that using the hybrid first-order-numerical model.

BIBLIOGRAPHY 214

215 BIBLIOGRAPHY [1] R. J. Adler, The Geometry of Random Fields. new York: Wiley, 1981. [2] Richard T. Austin and Anthony W. England, "Multi-Scale Roughness Spectra of Mount St. Helens Debris Flows," Geophysical Research Letters, vol. 20, pp. 1603-1606, 1993. [3] Richard T. Austin. Anthony W. England. and Gregory H. NWakefield. "Special Problems in the Estimation of Power-Law Spectra as Applied to Topographical Modeling," IEEE Transactions on Geoscience and Remote Sensing, vol. 32, pp. 928-939, 1994. [4] Bateman Manuscript Project, Tables of Integral Transforms. vol. 2, A. Erdelyi, Ed. New York: McGraw-Hill, 1954. [5] P. Beckmann and A. Spizzichino, The Scattering of Electromagnetic Waves from Rough Surfaces. Oxford: Pergamon Press, 1963. [6] T. H. Bell, Jr., "Statistical features of sea-floor topography," Deep-Sea Research, vol. 22, pp. 883-892, 1975. [7] T. H. Bell, Jr., "Mesoscale sea floor roughness," Deep-Sea Research, vol. 26A, pp. 65-76, 1979. [8] M. V. Berry and J. H. Hannay, "Topography of Random Surfaces," Nature, vol. 273, p. 573, 1978. [9] R. B. Blackman and J. W. Tukey, The Measurement of Power Spectra from the Point of View of Communications Engineering. New York: Dover, 1958. [10] S. R. Brown and C. H. Scholz, "Broad bandwidth study of the topography of natural rock surfaces," Journal of Geophysical Research, vol. 90, pp. 12 575 -12 582, 1985. [11] Fred M. Bullard, Volcanoes of the Earth, 2nd edition. Austin: University of Texas Press, 1984. [12] R. H. Cameron and W. T. Martin, "The Orthogonal Development of Non-linear Functionals in Series of Fourier-Hermite Functionals," Annals of Mathematics, vol. 48, pp. 385-392, 1947.

216 [13] J. Capon, "High-Resolution Frequency-Wavenumber Spectrum Analysis," Proceedings of the IEEE, vol. 57, pp. 1408-1418, August 1969. [14] M. F. Chen and A. K. Fung, "A numerical study of the regions of validity of the Kirchhoff and small-perturbation rough surface scattering models," Radio Science, vol. 23, pp. 163-170, 1988. [15] M. F. Chen, S. C. Wu, and A. K. Fung, "Regions of Validity of Common Scattering Models Applied to non-Gaussian Rough Surfaces," Journal of WaveMaterial Interaction, vol. 2, pp. 9-25, 1987. [16] Cornel Eftimiu, "Scattering by a rough dielectric interface: a modified WienerHermite expansion approach," Journal of the Optical Society of America A, vol. 7, pp. 875 —884, May 1990. [17] Cornel Eftimiu, "First-order Wiener-Hermite expansion in the electromagnetic scattering by conducting rough surfaces," Radio Science, vol. 23, pp. 769-779, September/October 1988. [18] Cornel Eftimiu, "Second-order Wiener-Hermite expansion in the electromagnetic scattering by conducting rough surfaces," Radio Science, vol. 23, pp. 1075-1084, November/December 1988. [19] Cornel Eftimiu, "Modified Wiener-Hermite expansion in rough- surface scattering," Journal of the Optical Society of America A, vol. 6, pp. 1584-1594, October 1989. [20] Cornel Eftimiu and G. W. Pan, "First-order Wiener-Hermite expansion in the electromagnetic scattering by dielectric rough interfaces," Radio Science, vol. 25, pp. 1-8, January/February 1990. [21] A. W. England. "The fractal dimension of diverse topographies and the effect of spatial windowing," Geological Survey of Canada Paper 90-4: Ground Penetrating Radar, J. A. Pilon, Ed., pp. 57-61, 1992. [22] C. G. Fox and D. E. Hayes, "Quantitative Methods for Analyzing the Roughness of the Seafloor," Reviews of Geophysics, vol. 23, pp. 1-48, 'February 1985. [23] Bruce L. Foxworthy and Mary Hill, Volcanic eruptions of 19809 at Mount St. Helens: The First 100 Days. U. S. Geological Survey Professional Paper 1249, 1982. [24] Adrian K. Fung, Microwave Scattering and Emission Models and Their Applications. Boston: Artech House, 1994. [25] A. K. Fung and K. K. Lee, "A semi-empirical sea-spectrum model for scattering coefficient estimation," IEEE Journal of Oceanic Engineering, vol. OE-7, pp. 166-176, 1982.

217 [26] L. E. Gilbert and A. Malinverno, "A Characterization of the Spectral Density of Residual Ocean Floor Topography," Geophysical Research Letters, vol. 15, pp. 1401-1404, November 1988. [27] Harry Glicken, Rockslide-debris avalanche of May 18, 1980, Mount St. Helens Volcano, Washington. U. S. Geological Survey Professional Paper 1488, 304 pp., 7 plates, 1989. [28] J. A. Goff, "Comment on 'Fractal Mapping of Digitized Images: Application to the Topography of Arizona and Comparison With Synthetic Images' by J. Huang and D. L. Turcotte," Journal of Geophysical Research, vol. 95, p. 5159, 10 April 1990. [29] Roger F. Harrington, Time-Harmonic Electromagnetic Fields. New York: McGraw-Hill, 1961. [30] J. Huang and D. L. Turcotte, "Fractal Mapping of Digitized Images: Application to the Topography of Arizona and Comparisons With Synthetic Images," Journal of Geophysical Research, vol. 94, pp. 7491-7495, 10 June 1989. [31] J. Huang and D. L. Turcotte, "Fractal image analysis: application to the topography of Oregon and synthetic images," Journal of the Optical Society of America: A, vol. 7, pp. 1124-1130, June 1990. [32] S. M. Kay, Modern Spectral Estimation. Englewood Cliffs, NJ: Prentice Hall, 1988. [33] Z. Li and A. K. Fung, "Scattering from a Finitely Conducting Random Surface," in Proceedings of the Progress in Electromagnetics Research Symposium (PIERS), Boston, 1989, p. 144. [34] B. B. Mandelbrot, The Fractal Geometry of Nature. New York: NV. H. Freeman and Company, 1983. [35] B. B. Mandelbrot and J. W. Van Ness, "Fractional Brownian motions, fractional noises, and applications," SIAM Review, vol. 10, pp. 422-437, 1968. [36] S. L. Marple, Jr., Digital Spectral Analysis. Englewood Cliffs, NJ: Prentice Hall, 1987. [37] W. C. Meecham and W. C. Clever, "Use of C-M-W Representations for Nonlinear Random Process Applications," in Statistical Models and Turbulence, J. Ehlers, K. Hepp, and H. A. Weidenmiiller, Eds. Berlin: Springer-Verlag, 1972. [38] W. C. Meecham, and W. W. Lin, "Reflection from random rough surfaces," Journal of Wave-Material Interactions, vol. 1, pp. 4-15, 1986. [39] William C. Meecham and Wuu-Wen Lin, "Radiation Reflected from Random Rough Surfaces," IEEE Journal of Oceanic Engineering, vol. OE-12, pp. 357 -361, 1987.

218 [40] William C. Meecham and Armand Siegel, "Wiener-Hermite Expansion in Model Turbulence at Large Reynolds Numbers," Physics of Fluids, vol.. 7, pp. 1178-1190, 1964. [41] J. Nakayama, H. Ogura, and B. Matsumoto, "A probabilistic theory of scattering from a random rough surface," Radio Science, vol. 15, pp. 1049-1057, 1980. [42] J. Nakayama, H. Ogura, and M. Sakata, "A probabilistic theory of electromagnetic wave scattering from a slightly random surface, 1, Horizontal polarization," Radio Science, vol. 16, pp. 831-845, 1981. [43].J. Nakayama, M. Sakata, and H. Ogura, "A probabilistic theory of electromagnetic wave scattering from a slightly random surface, 2, Vertical polarization," Radio Science, vol. 16, pp. 847-853, 1981. [44] C. Eric Nance, "Scattering and Image Analysis of Conducting Rough Surfaces," Doctoral Dissertation, University of Texas at Arlington, May 1992. [45] C. E. Nance, A. K. Fung, and J. W. Bredow, "Comparison of Integral Equation Predictions and Experimental Backscatter Measurements from Random Conducting Rough Surfaces," Proceedings of the International Geoscience and Remote Sensing Symposium (IGARSS '90), College Park, Maryland, 1990, pp. 477-480. [46] A. Nashashibi, K. Sarabandi, and F. T. Ulaby, "A Calibration Technique for Polarimetric Coherent-on-Receive Radar Systems," submitted to IEEE Transactions on Antennas and Propagation, 1994. [47] A. V. Oppenheim and R. W. Schafer, Digital Signal Processing. Englewood Cliffs, NJ: Prentice Hall, 1975. [48] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (FORTRAN Version). Cambridge: Cambridge University Press, 1989. [49] S. O. Rice, "Reflection of Electromagnetic Waves from Slightly Rough Surfaces," Communications in Pure and Applied Mathematics, vol. 4, pp. 361-378, 1951. [50] George T. Ruck, Donald E. Barrick, William D. Stuart, and Clarence K. Krichbaum, Radar Cross Section Handbook, 2 volumes. New York: Plenum Press, 1970. [51] D. Saupe, "Algorithms for random fractals," in The Science of Fractal Images, H.-O. Peitgen and D. Saupe, Eds. New York: Springer-Verlag, 1988, pp. 71 -136. [52] R. S. Sayles and T. R. Thomas, "Surface topography as a nonstationary random process," Nature, vol. 271, pp. 431-434, 1978.

219 [53] Thomas B. A. Senior, Mathematical Methods in Electrical Engineering. Cambridge: Cambridge University Press, 1986. [54] Alexander Slingeland. "Creating Surfaces with Mastercam Mill & the Techno Computer Controlled Milling Machine," EECS 499 Directed Study, University of Michigan, May 1993. [55] Steven R. Thiel, 'Determination of a Lossy Dielectric Material That Can Be Used in the Casting/Molding Process," EECS 499 Directed Study, University of Michigan, December 1991. [56] F. T. Ulaby and C. Elachi, Editors, Radar Polarimetry for Geoscience Applications. Norwood, MA: Artech House, 1990. [57] Fawwaz T. Ulaby, Thomas F. Haddock, and Richard T. Austin, "Flucuation Statistics of Millimeter-Wave Scattering From Distributed Targets," IEEE Transactions on Geoscience and Remote Sensing, vol. 26, pp. 268-281, May 1988. [58] Fawwaz T. Ulaby, Thomas F. Haddock, Richard T. Austin, and Yasuo Kuga, "Millimeter-wave radar scattering from snow: 2. Comparison of theory with experimental observations," Radio Science, vol. 26, pp. 343-351, 1991. [59] F. T. Ulaby, R. K. Moore, and A. K. Fung, Microwave Remote Sensing: Active and Passive, vol. II. Reading, MA: Addison-Wesley, 1982. [60] R. F. Voss, "Fractals in Nature: From characterization to simulation," in The Science of Fractal Images, H.-O. Peitgen and D. Saupe, Eds. New York: Springer-Verlag, 1988, pp. 21-70. [61] P. D. Welch, "The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms," IEEE Transactions on Audio and Electroacoustics, vol. AU-15, pp. 70-73, June 1967. [62] Michael W. Whitt, Fawwaz T. Ulaby, and Thomas F. Haddock, "The Development of a Millimeter-wave Network Analyzer Based Scatterometer," Radiation Laboratory Report 022872-1-T, University of Michigan, January 1987. [63] N. Wiener, Nonlinear Problems in Random Theory. Cambridge: MIT Press, 1958. [64] R. E. Ziemer and W. H. Tranter, Principles of Communications: Systems, Modulation, and Noise, 2nd ed. Boston: Houghton Mifflin Company, 1985.