027735-1-F Final Report for: NASA GRANT NAGW-2199 THEORETICAL AND EXPERIMENTAL MODELS OF THE DIFFUSE RADAR BACKSCATTER FROM MARS Work performed during July 1, 1990 through June 30, 1994 (plus a 6 month, no-cost extension) A.W. eJgld" Principal'Investigator Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Ph: (313) 936-1340 FAX (313) 763-1503 england@eecs.umich.edu May 19, 1995

2 Table of Contents I Statement of Objectives........................................................3 II Approach........................................................3 III Assumed Surface Statistics........................................................3 Multi-Scale Roughness Spectra of Mount St. Helens Debris Flows.....................5 IV Estimation of Power-Law Spectra........................................................9 Special Problems in the Estimation of Power-Law Spectra as Applied to Topographical Modeling........................................................ 10 V Em pirical M odel.............................................................................................................22 Design and Manufacture of Artificial Power-Law Surfaces..................................23 Radar Measurements..............4................... 48 VI Theoretical Models.........................................64 Rough Surface Scattering Models and Comparisons to Experiment.................65 VII Continuing W ork........................................79 VIII Acknowledgment 79........................................ IX References.79................

3 I Statement of Objectives Our general objective for this work was to develop a theoretically and experimentally consistent explanation for the diffuse component of radar backscatter from Mars. The strength, variability, and wavelength independence of Mars' diffuse backscatter are unique among our Moon and the terrestrial planets. This diffuse backscatter is generally attributed to wavelength-scale surface roughness and to rock clasts within the Martian regolith. Through the combination of theory and experiment, we attempted to bound the range of surface characteristics that could produce the observed diffuse backscatter. Through these bounds we gained a limited capability for data inversion. Within this umbrella, specific objectives were: (a) To better define the statistical roughness parameters of Mars' surface so that they are consistent with observed radar backscatter data, and with the physical and chemical characteristics of Mars' surface as inferred from Mariner 9, the Viking probes, and Earth-based spectroscopy. (b) To better understand the partitioning between surface and volume scattering in the Mars regolith. (c) To develop computational models of Mars' radio emission that incorporate frequency dependent, surface and volume scattering. II Approach There is never enough radar scattering data to invert Mars radar backscatter profiles to obtain unambiguous estimates of surface roughness - even if we possessed the scattering theory that would allow us to do so. Our approach to the problem has been to assume something about the statistical roughness characteristics of Mars' surface, and, based upon this roughness model, to develop a limited scattering theory or empirical law that could be used to relate observed backscatter to roughness. We are currently testing this approach by developing a classification algorithm based upon the empirical law developed in this investigation. Applying this algorithm to SIR-C data, we will attempt to classify Andean Volcanics according to their roughness. This application is funded by the SIR-C project and is not strictly a part of the Mars investigation. III Assumed Surface Statistics For roughness characteristics, we assumed a power-law surface. Several investigations have suggested that roughness profiles of many natural surfaces obey a power spectrum of the type,

Sz(f)= cf-P (1) [e.g., 7-8, 10, 22, 26, 28, 30-31, 34-35, 52, 60] where Sz is the power spectral density, f is spatial frequency, c is a scaling constant, and f is a constant that is related to the fractal dimension of the surface [e.g., 21]. We tested this hypothesis on several flows of the Mount St. Helens volcano and found that fresh flows in that environment do, in fact, exhibit a power-law spectrum over several decades of spatial frequency both below and above the range of spatial frequencies of observing radars. In contrast, water deposited sediments of Mount St. Helens ejecta do not obey a powerlaw. While aeolian deposition on Mars may be like fluvial deposition at Mount St. Helens and not obey a power-law, the assumption that Mars volcanic deposition does obey a power law gives us a place to begin. We have presented various of these findings in several forums [65-66] and published the work in Geophysical Research Letters [2]. The GRL paper is included on the following 4 pages.

GEOPHYSICAL RESEARCH LETTERS, VOL. 20, NO. 15. PAGES 1603-1606. AUGUST t. 1993 M1I'LTI-S('CA.E R()I'(;IIN.ES. SPI'('TIR.\ (F lMOUNT' SXT. tIELENS DE31BRIS FLOWS licllarl I 1.\ lstill andl.ltlionv \V. EnIglarnd IH(itiitoi L.(,, at (r. 1) )t. of Elc 1*t I icHI Frl 21 I 1 a(l 1I!ri)llt er Science. I[ni\versit of MIichIliTar..b..t,'act.\ r:; prie tec.ltri1nl allows surface struc- ext remelv olllr t errains. hlie sulrfaces t hat we examnie(ld tlure to h)e iterlrlt (l..1 as a <ni of siniusoidal components were located in the (ebl)ris avalanche wvest-northlwest of the faith (lifferring WaelengSths. lnowledge of the roughness volcano. along the North Fork Tolutle River valley (Figslpectrulm gives inisi2lit into t he mechanisms responsible for ure 2). MIost of the d(el)ris was dleposited (ulring the erupelectronlagnetic scattering at a given Nwavelength.. lea- tion of Mlount St. Helens on IS Mlav 1980. Since then. tile suredl spectra from 10-vear-oll p)rimary debris flowsv sur- deposits have undergone erosion bty windl and water. The faces at M\ount St. tIelens conform to a power-law spectral sites chosen for measurement where visually representative model. -iggesting that these surfaces are scaling over the of principal surface textures, accessible to logging roads. measul. range of spatial frequencies. MIeasured spectra and located in large homogeneous areas which should be from water-deposited surfaces deviate from this model. distinguishable in a radar image. Figure 3 shows one of the rougher surfaces. A geologic description of the debris Introduction avalanche is given by Glicken [1989]. Our objectives were to characterize the surface roughRecent investigations [Evans et al.. 19SS: Gaddis et al.. ness of the debris flows at scales smaller than, on the or1989] demonstrate that radar images are a useful tool in der of, and larger than the radar wavelength of common the study of volcanic terrain. Although there have been remote sensing radars. Two measurement techniques were experimental investigations of radar scattering by such used. A computer-driven. 2D laser profilometer recorded terrains, there have been few rigorous attempts to link surface height profiles of square grids with sides between the scattering cross section to quantitative surface rough- 8 cm and 1 m in length. The grids were sampled at interness. Scattering from many volcanic terrains is difficult vals ranging from 2 mm to 2 cm. We used surveying into model due to their extreme roughness and short cor- struments to measure elevation at grid points within larger relation lengths. Rough surface scattering theories of- square areas measuring 32 m on a side. Sampling intervals ten assume that surfaces possess a Gaussian autocorrela- within the larger squares were 1-4 m. tion function (and corresponding Gaussian roighness spec- The 2D laser profilometer system was designed speciftrum) [Chen and Fung, 1988; Thorsos and Jackson, 1989]. ically for this experiment. Its main component is a surIn contrast, mans natural surfaces are better described by veying electronic distancemeter (EDM), which uses an inpower-law spectra [Berry and Hannay, 1978; Sayles and frared laser to measure distance. The EDM is mounted on Thomas, 1978: Huang and Turcotte, 1989, 1990; England, an XY table, allowing the EDM to scan a surface area of 1992]. (These spectral forms are compared in Figure 1.) up to 1 m2. The EDNI laser has a spot diameter of approxRough surfaces that have structure over a wide range of imately 1.5 mm. In its most precise mode, the standard spatial scales may be described using the concepts of frac- deviation of the measured surface height is 3 mm. Surface tal geometry introduced by Mandelbrot [1983]. Random height profiles from a typical grid are shown in Figure 4. rough fractal surfaces have power-law spectra; the slope of While the two large rocks are quite prominent, a larger such a spectrum is related to the surface's fractal dimen- sample of the surface would show even larger "lumps." sion (Austin et al.. Special Problems in the Estimation of Larger-scale topography was surveyed. Square grids Power-Law Spectra as Applied to Topographical Model- with sides of 32 m were delineated by cables with markers ing, submitted to IEEE Transactions on Geoscience and at 1 m intervals. A self-leveling surveying level and staRemote Sensing, 1993. hereinafter referred to as Austin et dia rod were then used to measure the surface height at al.. submitted manuscript. 1993). In this report, we exam- intervals of 1-4 m, depending on the surface roughness. ine estimates of the roughness spectra of volcanic debris flows at several spatial scales. _ Debris Flow Measurements la 10'4 Surface roughness measurements were performed on sev- X eral debris flows near NMount St. Helens (Cascade Range, northwestern United States) during September 1990. Mt. St. Helens was chosen because of its ease of access and its I II:i~ Copyright 1993 by the American Geophysical Union......,. 0.1 1 10 100 Paper Number 93GL01875 Spatial frequency (m-) 0094-8534/93/93GL-01875$03.00 Fig. 1. A comparison of Gaussian and power-law spectra. 1603

1604 As.-\tin and(l LEnglarnl: [Rollilles, S.pectra of l)eblris ilows 0.9 ERS Elk 0.8 Rock _ 0.7 Flow JRS Debris Flow JRS~ohOstovr, _ 0.5 - ~ % %,0.4 - B.C., -'"!ash. A I SH 0 1 2 3 4 5 km, 0.0 0.1 0.2 0.3 0.4 0.5 0.6 ______ HDistance (m) Oregon Mount St. Helens North Fig. -1..-\ typlical lrofilorleter salan. slowiil o sulrface heighlt profiles. The ri(l nleasture's 10 C x 10 cIll. wNit i Fil. I. Iile.Jollnstonl Ridge site (JRS) is located about lard points spacedl 1 cmn atpart. l'otal slrfa(e iltight.ari10 krn rlort}liwest of thlie crater. Tile Elk Rock site (ERS) atior is 10.8 cut. is located farther west. abolit 1!) knm west-northwest of the crater. we assume in this study that the surfaces have isotropic Spectral Estimation statistics as a first-order approximation. Power-law spectra introduce unique difficulties in the Estimates of the roughness spectrum were calculated spectral estimation process. In a related paper (Austin et from both types of elevation data. Profilometer data suf- al., submitted manuscript. 1993). we describe how Capon's fered from errors which prevented the direct application estimator provides an estimate of the roughness spectrum of spectral estimation techniques. Since the measurements that is better than that given by Fourier-based estimators could not be repeated, it was necessary to apply corrective (such as the periodogram) in the power-law case. Capon's filtering. The procedure is fully described in Austin and estimator (described in Kay [198S]) essentially customizes England [1991]. a filter at each frequency of interest to minimize the total As previously noted, investigators have found that many power output subject to the constraint that the gain at natural surfaces are scaling. Linear profiles of such surfaces the frequency of interest is unity. In the power-law case have (over the range of measurable spatial frequencies) the low-frequency side lobes are reduced, thus Capon's esa surface roughness spectrum (power spectral density of timator suffers less from spectral leakage. which can make surface height) Sz that has the form of a power law: Fourier estimates insensitive to the power-law slope 3. Capon's estimator also has less variance than comparable Sz(f) = clf1-, ( 1 ) Fourier-based estimators. Capon's estimator, PCAP(f), is obtained by first calwhere f is the spatial frequency and c and d are constants, culating an estimate of the autocorrelation matrix Rzz, with 1 < 3 < 3. /3 is sometimes called the spectral slope; a whose elements are defined as power-law spectrum is linear with slope (-3) when plotted on a log-log scale. [Rzz]kl = E [Z[n]Z[n + k - 1]], (3) Power-law surfaces with isotropic statistics have a twodimensional roughness spectrum SZ2D that has a similar where Z[n] is the surface height and E[.] indicates an exform: pected value. Capon's estimator is then given by Sz2D(fr) = af~", (2) where f,. is the radial spatial frequency satisfying f2 = PCAP(f) - eHf e' (4) f2 + fy2, and 2 <' < 4. Because estimation of directional spectra requires a much greater quantity of data, where p is the maximum dimension of the autocorrelation matrix, A is the sampling interval, R-1 is the inverse of an estimate of the autocorrelation matrix, and 4y1i~1"' i wh [h 1 ej2nfA eturfA... j2w(p-.l)f A 1 (5);/l~r[' "'~':':'/'' *iin which f is the spatial frequency. eH is the Hermitian JCI 1.~-i.-,~, transpose of (5). 7=~!~''i;.;;?J:~ "- The bias of Capon's estimator is independent of N (the _~4 -number of points in a sampled row) but is dependent on _- -2.X ~..,,~~ p, the dimension of the autocorrelation matrix. For a data ar~ set of length N, an increase in p will result in reduced bias _ -j~ H a at the cost of increased variance. We set p 0.3N and calculated a Capon's spectral estimate for each row of each profilometer or survey grid. (The mean row height was Fig. 3. A sample surface scanned by the EDMI. The sur- subtracted before calculating the estimate.) The spectral face has been spray-paintedl to increase its reflectivity. estimates were then averaged over the rows of each grid.

.-\ustill aI(l [:lilarlld: I\('ln,,litrIe - ";)('(tra, t lht,'iiD I:''wl, 1t05 10 v.lmAalre( [t),11m re;..l tess'ct ra A-.1 m A o10 T- hile first site tilat we ex.r anailedl i.JRS; wais l(,oateI ill 10 tile (lel)ris flow area a.lja('telt t(, alnl ({le' wcst of l.,}llnston 10 Ridge Isee l ilir, 2). L'I'l i area of t le lft eiris flow \,h'xili,10 5_ =2 cr itedl lairk xvriat in-I in relit if.. ai ill m I I X:2 i 1 10 w - Eq. (6) mm tis ar(e wereX t )lie ialllar 1 hiri- tilow siirface~. i.e.. (Iieris 10 -8 A=5mm flow sutilfa(ces w}ho,( p)rew.ilt stats i s (1lle to tie r, tl(o'al 10' of material }) win(i an1 water. Spectra fioni tie ro= li0.1 1 10 100 est of t hese are shownT in Fiolure.-\A.'Ihic survey s)pectrumnl E 10_ is the curve extending lr roll 0.0 to i.25 n-i thc eigher E 10' frequllency spectra were compultel from roirolonieter scans. 2 102l Tile linear charact(er of these sp(ct ral estimates suiggests'A 3 " that a i)owver-law slectrunl is appropriate for this surfac-e. 10'-4 The equlation of tlie power-law fit to these estitnates is i - 1 cm'z(.f) = (:3.'1 x 10-4 )fl-2.34 (6) o. - - Eq. (8) m 10'7.10. -' (The hat denotes an estimate.) This fit is shown as a It10 V dashe(l line in Figure 5A. \\hile we do not know that the 0.1 1 10 100 spectrum follows this power law at the intermediate (un10 0 2m measured) spatial frequencies, we do think that it is sig10, A=2m C nificant that both the meter- and centimeter-scale spectra'210-2' - are well fit by the same power law. 3 We can convert (6) to an estimate of the surface (21D):.0'4 - spectrum, obtaining the following power law: 10~106 _ A=5mm SZ2D(f,) = (1.77 X 10-4)f334. (7) 1-0 - Eq. (10) 10 - In Figure 5B. the same survey spectrum is shown to10. gether with a profilometer spectrum collected on a similar 0.1 1 10 100 Spatial frequency (m') Fig. 5. (A) Survey and profilometer spectra from a JRS 10 a2m A primary debris flow. (B) The same survey spectrum with 2 10.~' a profilometer spectrum from an adjacent debris flow. (C) Spectra from a mixed terrain survey and a primary debris 4 Eq. (10) 10-4 flow scan at ERS. E 04. EDM noise floor The high-frequency half of each averaged estimate was dis- C 10-8 cmcarded to reduce the effect of aliasing. Averaged estimates 0.1 1 10 1 10 at a single scale or at multiple scales were then fit with a E 100o power-law function (1), where appropriate, using a min- A 1 m B imum absolute deviation fit. Profile (1D) spectra fitting a 10'.2 the power-law model were then converted to surface (2D) 10 spectra using the formulation given in (Austin et al., sub- Eq (6) mitted manuscript. 1993), in which the 1D and 2D roughness spectra are related through their corresponding auto- 1 EDM noisefloor correlation functions and the 2D autocorrelation function is assumed isotropic. 10.8 A5 Although this method of determining the surface spec- trum from linear profiles was made necessary by the pro- 0.1 1 10 100 filometer problems, it is not without benefit. One grid of Spatial frequency (ma) measured surface heights yields many linear profiles. The Fig. 6. (A) A profilometer spectrum (A = 1 cm) of an profile rows are not completely independent, but the row ERS sedimented surface, shown with the survey specspectral estimates can still be averaged to obtain a signifi- trum of the local mixed terrain, the power-law fit (10), cant decrease in variance, allowing us to set p (in Capon's and the EDM noise floor. (B) A profilometer spectrum estimator) higher to decrease the bias. This averaging (A = 5 mm) of a JRS sedimented surface, shown with the wxvould not he possible with a two-dimensional estimator, survey spectrum of the surrounding primary debris flow, given the same measured data. the power-law fit (6). and the EDM! noise floor.

ltV,.,, ia,'l t,, a }1. IIIe-.,rv,aw,..- A fit t,) cited (ietnlozisttrat, - hf. t'id r!tov,. 1'4-law..clt terir:,t r11()dt }~,..~.,..~t'it,,,7'.,., v1.-' W,,'.~t,"' l ra: r'lt ili tte rOftl,,,-> z ll i' ot',l,'atli, (if iris t,.-. a II kt)I ( -r Ilat tlral surface(,. \\W ar- tlr,.llirlL stU(}c a nii tIel hlltroil."'> ~~;' = b-),. t 1() -ovo! ttlol I t-i. ld. (l.\ip r illoill i'tIii Ild!:l 1X,,.I`) l',k I lit) (' V\1 ll.v lite w'est'f Elk l1)ok. Iletit tir tI r,,,,)'rdt i,)1. 1[\li.\l(')otiall alt1 lil) ii- l;.; f e tll,f }lirtro(kv debris mo(ids.,-\,Ist i for assist ai' }il t(he s, frl'f e ii('asuirvtie'lt s. 1a(1 v ((lfte1'e acros..1 ia11..,.{i)(t'd( plain. A\ profilonieter aolltoittis revi(wirs frl t,.i'r }lml>htftill ()llll(t(lilti. dll.(all tlolil()Ill }o t tic ()(lf islisillolill.l~xi(}l~l(}(ltile}Fr~e~t - o hlg estolls. I- 1i<- _t1lIl))( 4ti 1oller. (;Xclilt 5C'all'4)111 Otie ot' i Is,eYiris tiottI(ls yielded th~e spectral es- sitg(-e.! otis. Ii, work is sit puort,.i t tiller N.\S.\ ( rat I italite sliowi' at th1e ligli-fre(llencv etd of Figure C. The N A( \V-21 99 1 l a sre.\.\ Si1-(' )rO.jec('t tio.lire also shlows an estimate from a sirvey grid containitli )ot Ih terraitI tvl)es. ( ih (leIbris mounds were typically Ief ere'l(es sltiall(er thanl tlie 32 1n x:32 m survey area.) A power-law Auistin. R. T.. ad.\. \W. Enalan(d. Surface (liaracterizatio fit tlav be al)prol)priale liere as well and is shown in the of Volcanic Debaris Flows at.\1llt iple Scales. in Proct(:ding. figuire. The equations for this fit are of the 1991 Inrnational Gos ( notc nsng = x -4 2.51 (10) Symposium, pp. 1675-167~,. IEEE. New York. 1991. Sz(.f) = (4.98 x 104I.fI[-2't (10) Berry. SI. V.. and J. II. Ilannay. Topography of random surSZ2D(.fr) = (2.S6 X 10-4)fJ-3-51 (11) faces. Nature. 273. 57:3. 1978. Chen.,I. F.. and A. K. Fung. A numerical study of the regions Surfaces formed of water-deposited sediments within the of validity of the Kirchhoff and small-perturbation rough hummockyl debris had spectra that were markedly differ- surface scattering models. Radio Science. 2., 163-170. 1988. ent from those of the primary surfaces. For example, in England. A. XW.. The fractal dimension of diverse topographies Figure 6A. we see that the spectral estimate from a scan and the effect of spatial windowing, in Geol. Soc. of Callda of such a surface lies on the noise floor due to the lim- Papr 90-4: Ground Penetrting Radar. ed. J. A. Pilon. 57iteci precision of the EDMI. The true spectrum of this sur- 61. 1992. Evans, D. L.. T. G. Farr. J. J. van Zvl. and H. A. Zebker. face must lie below this level and therefore cannot form Radar Polarimetrv: Analysis Tools and Applications, IEEE a power-law trend with the survey spectrum of the local si Tralns. Geosci. Remote Sensing, 26. 774-789. 1988. mixed terrain. This difference in spectral behavior is con- Gaddis, L., P. Mouginis-MIark, R. Singer, and V. Kaupp. Gesistent with the visible character of the sedimented sur- ologic analyses of Shuttle Imaging Radar (SIR-B) data of face, which appeared quite flat. Similar remarks apply to Iilauea Volcano, Hawaii, Geol. Soc. Am. Bull., 101, 317the profilometer spectrum of a water-deposited sediment 332. 1989. at the Johnston Ridge site, which is shown in Figure 6B Glicken, H., Rockslide-debris avalanche of.lMay 18, 1980, with the survey spectrum of the surrounding primary de- Mlount St. Helens Volcano. Washington (USGS Prof. Paper bris flow. 1488), 304 pp., 7 plates, 1989. Huang, J., and D. L. Turcotte. Fractal Mapping of Digitized Conclutsions Images: Application to the Topography of Arizona and Comparisons With Synthetic Images, J. Geophys. Res., 94, Spectral estimates derived from multi-scale topographi- 7491-7495,1989. cal measurements of several debris flow surfaces at Mount Huang, J., and D. L. Turcotte, Fractal image analysis: application to the topography of Oregon and synthetic images, J. St. Helens show differences that correspond to observed Opt. Soc. Am. A. 7. 1124-1130,1990. surface characteristics. The extremely rough surfaces affil- Malldebrot. B. B.. The Fmctl Geometry of Nature, 468 pp., iated with 10-vear-old primary debris flows possess roughW. H. Freeman and Company, New York, 1983). ness spectra that are well modeled by a power-law func- Sayles, R. S., and T. R. Thomas, Surface topography as a tion at certain length scales, implying that these surfaces nonstationary random process, Nature, 271, 431-434, 1978. are scaling over the measured range of spatial frequencies. Thorsos, E. I., and D. R. Jackson, The validity of the perSurfaces created through sedimentation have spectra that turbation approximation for rough surface scattering using lie below the noise floor at higher spatial frequencies and a Gaussian roughness spectrum, J. Acoust. Soc. Am., 86, thus do not fit the power-law model. 261-277, 1989. While we are not yet able to develop a classifier of volcanic terrains based on surface roughness spectra, this R. T. Austin and A. W. England, Radiation Laboratory, study shows that a power-law spectrum is a reasonable 3228 EECS Bldg., University of Michigan, Ann Arbor, MI descriptor of some volcanic surfaces. The scattering prop- 48019-2122. erties of such surfaces are not well understood. There has been some application of standard scattering models to surfaces of this type, but routgher volcanic terrains often (Received February 11, 1993; violate assumptions made in these models (e.g., rms height revised May 14, 1993; < radar wavelength). The present study and the others accepted June 3, 1993.)

9 IV Estimation of Power-Law Spectra Having measured topography at scales from millimeters to 32 meters on several volcanic flows, our next task was to estimate the roughness spectra of these topographies. We immediately found that simple Fourier-based estimators, like the periodogram, produced greatly distorted spectra. The low frequency leakage of a power-law surface produced an artificial 3 in equation (1) that tended to lay near 2. Many of the early papers on the spectra of natural surfaces [e.g., 1] have suggested that the f for natural surfaces always lay near 2. We believe that this conclusion may have been the result of the leakage phenomenon that we observed. In fact, while the Mount St. Helens debris flows obeyed a power-lay, their f always lay between 2.31 and 2.51 - which are significantly different from 2. An artificial surface produced with a f of 2 would appear too rough at the higher spatial frequencies when compared to the Mount St. Helen surface profiles. We presented these findings at the Lunar Science Conference [67] and published the results in the Transactions of Geoscience and Remote Sensing [3]. The TGARS paper is included on the following 12 pages.

IEEE FRANSACTIO)NS ON GEOSCIENCE AN[) REMOTF I\, L;' t t i, -. Special Problems in the Estimation of Power-Law Spectra as Applied to Topographical Modeling Richard T. Austin, Student Member, IEEE. Anthony W. England. Senior Member, IEEE, and Gregory H. Wakefield. Member, IEEE Abstract-An increasing number of topographical studies hld of spatial scales may also be described using the concepts that natural surfaces possess power-law roughness spectra. of fractal geometry introduced by Mandelbrot [14]. FracPower-law spectra introduce unique difficulties in the spectral estimation process. We describe how an improper window tal objects are continuous but not differentiable: the fracestimation process. WAe describe how an improper windowm' choice allows leakage that yields a spectral estimate that is in- tal dimension Df is a real-valued measure of how a line sensitive to the spectral slope. In addition, the commonly used (surface) fills a plane (space). For a one-dimensional Fourier-based spectral estimates have higher variances than (l-D) surface profile. Df takes on values between 1.0 other available estimators. Higher variance is particularly (smooth and differentiable) and 2.0 (plane-filling). Twoproblematic when data records are short, as is often the case in remote-sensing studies. We show that Capon's spectral es- dimensional (2-D) surfaces have Df between 2 timator has less variance than Fourier-based estimators and Adler [15] shows that the surface profile created by the measures the spectral slope more accurately. We also show how intersection of a plane and a 2-D fractal surface is itself estimates of a 2-D roughness spectrum can be obtained from fractal with a fractal dimension equal to that of the 2-D estimates of the I-D spectrum for the isotropic power-law case. surface decreased by one [2]. Random rough fractal surfaces have power-law spectra; Mandelbrot and Van Ness I. INTRODUCTION [16] and Voss [17] derive the relation between Df and the S PECTRAL analysis is a tool of increasing importance spectral slope B of a 1-D profile of a rough fractal surface: in the characterization of natural surfaces. One- and S- ( two-dimensional spectral estimates are useful because they Df = (2) allow surface structure to be interpreted as a sum of sinusoidal components with differing wavelengths. A num- Thus, the spectral exponent satisfies 3 > / > 1 for 1 < ber of studies have found that spectra computed from one- Df < 2. dimensional surface profiles may be reasonably modeled Although synthetic topographies appear most realistic using a power-law spectrum of the form with Df 5 2.2 [14], [17], yielding linear profiles of Df = 1.2 and / = 2.6, several studies [1], [4], [7], [8] have Sz(f) = C lfl- (1) reported measured values of A (the hat denotes an estiover some range of spatial frequencies. Sz(f) is the power mate) near 2.0, leading some investigators to suggest that spectral density of the surface height random process; it spectral slopes of 2.0 are charactenstic of natural topoghas units of meters2/meters-' or meters2 per unit spatial raphy [1]. One explanation for the discrepancy between frequency. The spectral exponent B indicates the relative measured A and inferred values of B is that spectral analcontribution of different wavelength components. Higher ysis of nonstationary regions tends to result in A near 2.0 values of / indicate that the surface roughness is mostly [9], [18]. We show in this paper that leakage in spectral due to low-frequency components, while lower values de- estimators may also yield A that cluster near 2.0. note significant high-frequency contributions. B is also While leakage is a well-known problem of classical known as the spectral slope; a power-law spectrum is lin- spectral estimators such as the periodogram, common ear with slope -,/ when plotted on a log-log scale. The Fourier-based spectral estimators also suffer from high roughness amplitude c is a multiplicative factor scaling variance. The usual remedy is to average spectral estithe roughness at all spatial frequencies. Applications of mates derived from multiple data sets or multiple segpower-law spectra include studies of regional topography ments of a long data set. This is not always possible: in [1]-[5], rock surfaces [6], seafloor morphology [7]-[11], remote-sensing studies, surface profiles are often only one and other surfaces [12], [13]. of many "ground-truth" parameters to be collected. As a Rough surfaces that have structure over a wide range result, investigators may be faced with the problem of estimating the surface spectrum using a limited data set.' In Manuscript received May 29, 1993; revised March 18, 1994. This work was supported in part by the National Aeronautics and Space Administra-'In the case of I-D surface profiles, we find that a typical record consists tion under Grant NAGW-2199 and by a shared SIR-C project. of at most 500 points; 2-D profiles are even shorter, rarely exceeding 50 The authors are with the Department of Electrical Engineering and Comrn- x 50 points. Numbers of this magnitude are considered small from a staputer Science, University of Michigan, Ann Arbor, MI 48109. tistical inference viewpoint, where samples of thousands or hundreds of IEEE Log Number 9402672. thousands of points are more common. 0196-2892/94$04.00 ~ 1994 IEEE

this paper. we describe a spectral estimator that reduces 102 the leakage problem and predicts c and 3 with reduced E... WT variance in the case of short data runs. 101' -\ - pwer-law. -2.8 II. SPECTRAL ESTIMATES FROM LINEAR PROFILES 10 Uz 10'L i. Fourier-based spectral estimators (e.g.. the periodo- j 102, gram) render an estimate of the power spectral density ~ 3 lI whose expectation or mean value is a convolution of the i 1, spectrum of the sampled profile and the Founer transform 10 "- i I of a window function introduced by the finite extent of 0.001 0.01 0.1 the sample profile. Various window functions are used; normalized spatial frequency, dimensionless common windows include rectangular. triangular. and Fig. 1. Spectral domain window functions', 10() tor the periodogram Hanning windows. The selection of a particular window and WH (13 for the modified periodogram sith Hanning windowv. both with N = 256. The consolution of the Aindow tunction with the true specis based on a trade-off between spectral resolution and trum (e.g.. the power-lav spectrum shown as the dashed line) man result spectral leakage. in a spectral bias (an inaccuracy in the mean salue ort the estimate). Spectral leakage refers to the inaccuracy of a spectral estimate at a given frequency due to convolution with the are spaced at intervals A. The periodogram spectral estitransformed window; i.e., spectral power "leaks" from mator PPER (fa) for the sampled surface height is given by nearby frequencies according to the shape of the trans- Kay [20]: formed window. We show spectral domain window func- 1 - I tions corresponding to a periodogram and a modified pe- PPER(fA) = - Z Z[xr] exp (-j 2rfan) (4) riodogram with Hanning window in Fig. 1. At each point N n = or frequency at which the spectrum is to be estimated, the The frequency is written f as a reminder that (4) is given true spectrum (the dashed line in Fig. 1) is weighted by in terms of a normalized frequency f, where the window function centered on that frequency and then f 5 averaged over frequency. Contributions from frequencies f ( ) other than the frequency of interest constitute spectral and f satisfies - f <. (Most spectral estimation leakage. We can reduce spectral leakage by choosing an texts use normalized frequency; we will follow this conestimator whose window has sidelobes that decay quickly vention until a later section when we compare spectra ob(such as WH in Fig. 1), but such a window will have a tained from sampled profiles with differing A.) broader main lobe, decreasing the spectral resolution. (See While the periodogram may be evaluated for any f.l [19] for a discussion.) In the power-law case, leakage less than (the normalized Nyquist frequency). it its often from the large peak at low spatial frequencies affects the calculated using a fast Fourier transform (FFT) and conestimate at upper spatial frequencies (making the slope sequently is evaluated only at discrete frequencies fai = shallower) unless a window function with very low side- ilN, where i = -N12, ~ ~, -1, 0, 1, 1,~ ~ ~, N/2. To lobe levels is selected. We now show this result for the evaluate the periodogram using an FFT at additional frecase of 1-D spectral estimates based on 1-D surface pro- quencies, we can "zero-pad" the original data series to files. length M by adding M - N zeros to the end. The periodSpectral estimates derived from sampled data suffer ogram will then be evaluated at frequencies fi = ilM. from aliasing if the sampled process has spectral com- Most FFT algorithms require that N be an integer power ponents at frequencies greater than the Nyquist frequency of two. Series for which N is not an integer power of two fc, must be zero-padded up to the next higher power of two. 1 The periodogram is an unreliable estimator of the power fc = 2 (3) spectral density because it has a variance that is equal to 2A ~the square of its expected value, independent of N [20]. where A is the sampling interval. In the sections that fol- The usual practice in the field of spectral estimation is to low, we will assume that aliasing is not present. Simu- average multiple periodograms to reduce the variance, but lated data sets will be bandlimited to prevent aliasing. this may not be possible if the quantity of data is limited.2 Modifications due to aliasing will be presented in a later We show later in this section that PPER is a poor estimator section. in the power-law case even if the variance problem were absent. A. Fourier-Based Estimators 2For example, suppose the data record is 1000 points long. To reduce the variance by a factor of eight, we can divide the data record into eight Consider a one-dimensional surface Z(x), where Z, the segments of length 125, calculate eight spectral estimates, and average the surface height, is a single-valued function (i.e., there are 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 no overhangs). Suppose that the surface height has been than eightfold if the segments are not uncorrelated. If the entire data record measured at N locations xi (i = 0, ~ ~ ~, N - 1), which has only 125 points, this procedure is not useful.

4'1 e I,~T~R.A\,+( 1(EE'RA )''C.tI,\, R)\ t'lt:, \I) K[ M il - -\l\c.' \ J L': A modified penriodogram PHAN (fa) may be obtained b - — 2.4 multiplying the data senies by a Hanning window before E 1 2.0.oz calculating the spectral estimate:. C 1021 U ) ] 10' io'' nHltl] =, I - cos ( I l | (7), 0.001 0.01 0.1 replaces the rectangular window (WRR[n] = 1) implicit in normalized spatial frequency, dimensionless (4), and an extra normalization factor is introduced: Fig. E xpected value of the penodogram spectral estimator for powerlaw surfaces with 2 = 2.0 and 2.4 and c., = 1.0 10W'. The estimates are 1 I - I based on 256-point profiles that were zero-padded to 16 384 points so that u = - Z W[n]. (8) the penodogram would be evaluated at many frequencies. Solid dots IndiN n=0o cate the expected values of the penodogram when evaluated without zeropadding. (Equation (6) follows the formulation of Welch [21] using a single data segment.) The transform of the Hanning where window wH[n] has sidelobes that are lower and a main lobe that is roughly twice as wide as that of the transform WH(fA) =-1 (b + - b2 + - b3i (13) of a rectangular window WR [n] of the same width. Con- 4NU 2 2 volution with a transformed Hanning window results in less leakage from low spatial frequencies (because the b, = N) (14) sidelobes are smaller) at the cost of decreased spectral res- sin (Trf) olution (smoothing due to the wider main lobe). We must Cin 1 1 ) consider the smoothness of the spectrum and how quickly sin lo r i - N - N it rolls off and then decide whether the leakage reduction is worth the loss in resolution. In the power-law case, ( F I reducing the leakage seems more important. sin 7 fA NWe now examine the expected values of these two es- zLI timators for a power-law surface spectrum. The expecta-. ( 1 tion of the periodogram is given by [20] + N 1 N c 1/2 b -(16) E[PpER(f)] = Wr(fa - )Sz () dt (9) sin r f+N 1j} where W7Af) is the transform of a triangular window: and U is given by (8). 1 / )2 To show the performance of these Fourier-based estiWT(f) = } sin RfN (10) mators, we calculate their expected values using (9) and N sin 7rf (12). We avoid the singularity at the low end of a pure power-law spectrum by modeling the exact surface specand SZ,(f,) is the normalized form of the true power spec- trum as a Rayleigh function at low spatial frequencies: tral density of Z(x) chosen such that I Sza(fa) dfA = i Sz(f) df: (luI2)jfkl exp (-f A IfA ( 0o.0001 Sz,(fA) = SAf) Si(fa) = ctfIB - -1 0.0001 I fJ S A -'c(fA) - where the parameters a and a are chosen by enforcing the =Sa (11)) continuity of Sza and its first derivative at i faa = 0.0001. zt(fA) CAf. () Fig. 2 shows the expected value of the periodogram The expected value of the modified periodogram esti- spectral estimator calculated from theoretical spectra with mator PHAN has a form similar to (9) [23, p. 553]: spectral exponents of 2.0 and 2.4. The oscillatory behavI/2 ior is due to the sidelobes of WT(fA). These oscillations E[PHAN(fA)] = W(A - )S() d (12) are not visible if the periodogram is evaluated only at fre-t I/2 ~~~~quenciesf = iN (i.e., without zero-padding) as indicated

Al.STIN.'. PROBLEM\b iN ESTIM-ATI()N ()F PO)%ER-LA, SPECIRA -. E 101.. I — 1_-.2 ---- 24i "- 10 E 100 _ - — _ ---. __ ~~~ 20 ~~~~~~~ 10 N 120 N.. 10........ 1. -- -2.8 1, 10~]..... 1.... I.... 1=.2.81EO~,12 12. — ~~~~ ~ ~~~~~~~~~~~~~~~~.0.1~. i~~~~~~~104<-~~~~,104 10~~~~~~~i'-...'.....'-.:2 10 1u0.01 0.1 0.01 0.1 normalized spatial frequency, dimensionless normalized spatal frequency, dimensionless Fig. 3. Expected value of JPER for five 256-point power-law profiles eval- Fig. 4. Exact values of the power-lab spectra used in calculations of the uated at frequencies corresponding to the case of no zero-padding. expected values of the spectral estimates. in the figure. Since it is difficult to fit power-law functions TABLE I to such oscillatory estimates, we sampled the expectations POWER-LAW' PARAMETERS DERIVED FROM FITS TO THE EXPECTATIONS OF of periodograms of surfaces with five different spectral PERIODOGRAM (1/N < f. < 0.2). MODIFIED PERIODOGRAM (4/N < f < 0.5), PREWHITENED (I/N < fA < 0.2). AND CAPON'S (l/(2p) < fA exponents at frequencies corresponding to the case of no < 0.5) ESTIMATORS zero-padding (Figure 3) 3 These curves may be compared T tic EPE EPAI E[pw] EICAP 4 to the exact spectra in Figure 4.4 While the exact spectra, =.10 6to-. ~,7 = -] c -9.9 t. o-' | e., = 8.42 10'- c,, 10.10-6 vary in slope, the periodogram estimates are insensitive 3 =1.2 =1.231 = 1.201 = 1.255 = 1.198 ca = 1.0.-10- c, 1.04.10' a 9.99.10-:a 8.63.10' =1.01 ~ 10-' to the slope for several values of 3. Fitting power-law 8 = 1.6 = 1.656 = 1.601 = 1.643 = 1.598 functions to the linear portions of these estimate expec- c0 -.6 1 -.Oi 0 = 2.0 /0 = 1.982 0 = 2.001 6 = 2.042 0=1.998 tations (fa < 0.2), we see (Table I) that the estimated cA=. l10' C=1.87 10' 6& 9.98 10- CA=865 10' e=1.02 10' spectral exponent; clusters near 2.0 for 2 < (3 < 3. C.-24' 2.032 - =-2.e&99640-2 =2.43 - c = 1.02.3D0=: — 10 383-10-4 /:a =9.96.iC 10c: =8.43-10-* /:a - 1.02'10-' These values of f correspond to fractal dimensions that 0 = 2.8 = 1.986 = 2. = 2.859 0=2.796 are typical for natural terrains (1.0-1.5 for profiles, 2.02.5 for surfaces). Slope insensitivity is not a problem with PHAN(fA), lobe and first two sidelobes do not overlap the low-frewhose expected value is shown for various 3 in Fig. 5. quency peak (fA, > 0.015625). (Power-law fits to the upThe expected values of these estimates closely approxi- per regions of these estimates are listed in Table I). While mate the exact values for frequencies at which the main PHAN has a mean value that closely approximates the exact 3.. ^ ~~~~~~~~~~~spectrum this estimator suffers from a deficiency that is 3We emphasize that Fig. 3 shows expected values, E[5PER(f)], (i.e., spectrum, this estimator suffers from a deficiency that is mean values) of the periodogram estimator, calculated using (9) at fre- common to Fourier-based estimators: undesirably high quencies corresponding to the case of no zero-padding for random profiles variance. We examine the variance of PHAN(A) in a later Z(x) having power-law spectra given by (17) with five different values of section.6 t. Plots of the estimator PPER itself would show wide variance and departures from linearity (similar to the estimates shown in Fig. 8) and would vary for different realizations of the random process, i.e., for different sam- B. Prewhiening ples of a given surface. B. Prewhtenng The expected values of the estimator appear linear at the low-frequency The periodogram (Fig. 2) and, to a lesser degree, the end because we have assumed in (17) that the true surface spectrum is modified periodogram (Fig. 5) suffer from spectral leakpower-law (with linear slope) down to a frequency of fa = 0.0001-a fre modified peodogram (Fg. 5) suffer from spectral leakquency too low to be accurately measured by the sampled profile segment. age in the power-law case, especially for spectra with We used this assumption because surface profiles are often too short to 2 < ( < 3. Fox and Hayes [9] and Gilbert and Malinshow very long wavelength roughness components. (This was certainly true for our field work.) verno [10] avoid the leakage problem by using prewhiten4The expected values of the periodogram estimator in Fig. 3 have a higher ing procedures in which they modify the surface height total energy level than the corresponding exact spectra in Fig. 4 due to the spectral leakage that is also responsible for the slope insensitivity. Remem- 6We have also exmined the Blackman-Tukey spectral estimator T as ber that E[PPER] is the convolution of the true spectrum Sza with the trans- described in Blackman and Tukey 22]. For proc es other than white >.... ~~~~~~~~~~~~~described in Blackman and Tukey [221. For processes other than white form of a window function introduced by the finite surface sample [Wr in noise the BT estimator decreases variance but increases bias as compared T noise, the BT estimator decreases variance but increases bias as compared (9)]. If the estimator were perfect, WT would be a delta function, and the to the periodogram (Kay [20, p. 80). We did not include in our coestimate at frequency f would reflect only power in the true spectrum at prioda ( [0p 8e i o icl e n o rthat frequency. In the periodogram case, however, WT is given by (10) (a parison because it is even more susceptible to spectral leakage (and therethat frequency. In the periodogram case, however, WT IS given by (10) (a fore, slope Insensitivity) than PPER. The expected value of PBT is equal to sinc2 function), so power from many frequencies above and below f affect fore, slope insensitivity) thanp The expected value of PT is equal to (f). Leakage from lower frequencies is particularly significant in the the expected value of P.ER convolved with another window, as derived in PPER f Leakageyrom lower prequencies is particularly signifcant. 98]: power-law case. y p SFor a given single periodogram, ( may indeed take values between 2.0 / 112 and 3.0, due to the variance of the periodogram. Long data records do not E[PBT(fa)] = 3 W(fA - ~)E[/Ips(D)] dl. reduce this variance; the variance is independent of N (see [20, p. 422]). One may reduce the variance by segmenting the long data record as de- Since the BT estimator has an expected value that is the result of the conscribed previously, but this process will lead to the expected values shown volution of the true spectrum with two window functions, its value at a in Fig. 3. Sampling PPER at frequencies other than fa, = i/N could lead to given frequency can be greatly corrupted by leakage from lower frequena variety of incorrect slopes, as can be seen in Fig. 2. cies.

iEFE. TRAsACTIO(, Nb I, G;E(St.CIENCE AD KREI()TF SECNI.; \,OL 3_' (t; JIL:. E,E~ 2 ~~~E -. —.-'=-0.4' —0.8 1 -2 -1.2 -— =2.4 2-i.. 10 - o.' ------.''...-1.6.-2.8... 100 —.- ~'2.0 V._~ - -. - g 10 ui 10- - -10 L0 -6_ 0.01 o,1 O.01 0.1 normalized spatial frequency, dimensionless normalized spatial frequency, dimensionless Fiz. 5. Expected valueso POf' (ta for five different 256-point power- Fig. 6. Expected values of P\, (f) for five different 257-point pow\er-law Fi.-. Ep.oiec.te salu ies ot r PHA\(f4 o 638 oit. profiles. The profiles were zero-padded to 16 384 points. law profiles. The profiles were zero-padded to 16 384 points. dalu to have a relatively flat spectrum. perform spectral terms of the first differences of the sampled data: estimation, and then correct the estimate according to the N — -I original modification. ^I _ = Z We will examine here the simplest prewhitening ap- N n=o proach similar to that used by Gilbert and Malinverno. If the power spectral density of the surface height is defined (21) as [24, p. 270] 1 T/2 2 The expected value of the prewhitened estimator may be S7f) = AA E [ 2Z(x1) exp (- -21fr) d | derived by a procedure similar to the derivation of (9): T-.o T -T2 1/2 (18) E[Ppw(f4)] = 2 - 2 WT(fa - O a(0 then the derivative of the surface height dZ(x)Idx will have [1 - cos (27rt)] dt. (22) a related spectrum: r ratieT/2 2d 2 Using (22) and the modified power-law spectrum defined Sz (f ) =dlim - TjE j | in (17), the expected values of the prewhitened estimator S T cm T I -T2 d Iare calculated and are shown in Fig. 6. These estimates 1 rI S T/2 2 are based on 257-point profiles. lim -E j2irf Z(x)exp (-j2wf) dX Power-law functions were fit to the linear portions of T -a T - T/2 these estimate expectations (f, < 0.2), and the parame= 4r 2f2SZ (f ters ci' and t' were converted to e and, using (20) and 2 = C'/(4r2). The results, given in Table I, show that the = 47c f - 2 expected value of the spectral slope is fairly accurate beSzp(f) = C' If (-19) cause the spectral leakage has been largely eliminated. A () Hanning window may still be useful in practice (since we and we see that the spectral exponent y' of the derivative do not know a priori that a spectrum is power law), but process takes on values between -1 and 1 for 1 < ( < we see that leakage is not a problem for the prewhitened 3. The spectrum of the derivative process is more nearly estimator in the power-law case. The variance of PP (fi) flat (it is white noise for (' = 0). Since there is no pro- will be examined in a later section. nounced peak at low frequencies, spectral leakage is less of a problem. C. Capon's Estimator Gilbert and Malinverno approximate the derivative by An alternative procedure that is more direct than the taking first differences of the sampled data: Z' [xi] = prewhitened periodogram is the use of Capon's estimator Z [x + - Z [xi. They then apply a Hanning window and [25] (the so-called "minimum variance spectral estiuse a periodogram to obtain a spectral estimate of Sz'I(fA). mator") described in Kay [20, chap. 11]. This estimator After averaging multiple periodograms, an estimate of the was developed originally for geological signal processing spectral exponent 0 is obtained from an estimate of,': problems in which the number of sensors (and therefore = (' + 2. (20) the number of spatial samples) were extremely limited and *posed restrictions on what could be inferred from standard We now examine this procedure for the limited data case spectral analysis. Capon's estimator essentially customin which multiple periodograms are not available for av- izes a filter at each frequency of interest to minimize the eraging. We omit the Hanning window for simplicity. total power output, subject to the constraint that the gain We define a prewhitened spectral estimate Ppw(fa) in at the frequency of interest is unity. Therefore, the filter

ALSrI\ e: ai PROBLEM N1\ ESTIN1MATION OF POW(ER-LAA SPECTRA may be asymmetric in sidelobe level according to the shape of the signal spectrum: in the power-law case, the E 10 O -.12 - -- P2.4 sidelobes are adjusted to reduce the leakage from low- i I - - |-... 0 10'1'' -2.0 frequency components.. 10'2 Capon's estimator PCAP(fA) is obtained by first calcu- 1 -. lating an estimate of the autocorrelation matrix Rzz. whose 10..o elements are defined as, [Rzz], = E[Z*tIn]Zln + ill. (23) a ] We can estimate the elements of the autocorrelation ma- 001 0.01 0.1 trix using the modified covariance method, described in normalized spatial frequency, dimensionless [20]: Fig. 7. Expected values of PC,,p(f) for synthetic profiles with vanous I r N - X using an autocorrelation matnx of dimension p = 70. Estimates are eval=[R(N -p) L Z[n - i]Z[n -j] uated at the same 16 384 frequencies used in Figs. 2, 5. and 6. 2 (N- p) = p N- I -p Zln - D. Variance Comparison + Z[n + i]Z[n +1j] (24) Capon's estimator, the modified periodogram with Hanning window, and the prewhitened periodogram all Capon's estimator is then given by produce parameter estimates whose expected values closely approximate the exact values for a power-law PcAp(fa) = (25) Pspectrum. We now compare the estimators in terms of e Rzz e variance. The exact expression for the variance of the modified periodogram involves higher order moments and where p is the dimension of the autocorrelation matrix and is usually evaluated only for a Gaussian white-noise case e = [I ei2Trf e 4TwfA... ed(p - l1g] T. (26) [21]. The statistical properties of Capon's estimator for time series data are not known [20], although some results Since (24) is an unbiased7 and consistent estimator of have been derived for array data. To avoid these problems the covariance matrix, we can evaluate the expected value and test these estimators under actual-use conditions, we of Capon's estimator by substituting the exact covariance generated short synthetic topographic profiles of known matrix (obtained from the exact theoretical spectrum) into statistics using a spectral synthesis algorithm. We then (25). This procedure was used to calculate the expected compare parameters obtained from Capon's, the modified values of PCAP (f) for five values of 3 using a covariance periodogram, and the prewhitened periodogram estimates matrix of dimension p = 70. The expected values of these with those of the spectra used to create the profiles. estimates (Fig. 7) have good agreement with the exact Our spectral synthesis algorithm closely follows [26] spectra over a wide range of spectral slopes. Parameters and consists of the following steps. 1) Generate a set of of power-law functions fit to the expected values of discrete Fourier amplitudes that, when squared and mulPCAP (f) for f, >-.1I(2p) are compared with the exact tiplied by 1/N, satisfy the desired power law. While [26] spectral parameters in Table I. used a pure power-law spectrum, we use the modified The bias of Capon's estimator is independent of N be- spectrum (17) to assure a zero-mean surface height. 2) cause (24) is an unbiased estimator of lAzz. However, the Multiply the amplitudes by a Gaussian random variable bias does depend on p ( pA is the longest lag for which such that the mean value of the amplitude satisfies the the covariance is estimated). For a given data set, an in- power law. 3) Generate a set of complex discrete Fourier crease in p will result in reduced bias at the cost of in- coefficients (i.e., the FFT of a real surface) using the rancreased variance. Thus, N indirectly influences the bias domized amplitudes and a uniformly distributed random because the range of p is constrained by N. We generally phase, enforcing symmetry conditions such that the infound p = 0.3N to be a good compromise between vari- verse FFT will be real-valued. 4) Calculate the inverse ance and bias. FFT of the coefficients, resulting in a synthetic surface In the power-law case, the bias of Capon's estimator profile. was found to increase noticeably for spatial frequencies Using this algorithm, we generated 10 synthetic 64corresponding to wavelengths much longer than pA. We point surface profiles with (3 = 2.4 and ca, a, and a chotherefore discard calculated values of PCAP(fa) for spatial sen as before. We then calculated 1) Capon's estimates frequencies below 1/(2p). using a covariance matrix of dimension p = 20, 2) modified periodogram estimates using a full-width Hanning 7Equation (24) is approximately but not strictly unbiased because all lags window, and 3) periodogram estimates based on the first are not weighted equally and overlapping lags are used. A variety of co- differences of the surface profiles. We see in Fig. 8 that variance estimators have been examined, and (24) was found to be the. preferred estimator for short data records. See Kay [20] or other texts for the modified periodogram and prewhitened periodogram more information. estimates have noticeably greater variance.

EIE TR.ANsA-c('T I(\ i) GEOSCI.\NCF -+\XD RIF:,-!\ N. -,\ Jt Ni:~, 4 -nsestimator P 20 a measured surface profile. In remote sensing studies. esE 10' timates of the surface spectrum in terms of real spatial t 102 - -..'.....:-.... frequenc f(in meters- ) are necessar} so that 1) surtace a 10 ~ [Hanningl features may be compared to the electromagnetic wave- 1.2':. " length and 2) spectral estimates derived from profile,, with 10: different sampling intervals.A ma! be compared or comVP 10whi bined into a composite spectrum. ~ 104l * t w + WijR" The following expressions correspond to (4). (9). (10). a,, (6). (12). (13). (14). (15). (16). (21). (22). (25). and (26). 2 0i. 1 i; respectively, for the case of real (nonnormalized) fre0.1 normalized spatial frequency, dimensionless quenc f. where satisfies -'( A) < J ~ I (2A): Fig. 8. Capon'. modified penodogram. and prewhitened periodogram A spectral estimates for 10 different power-law profiles. The profiles were 64 PPER(f) = - Z [, exp -) (27) points long. The modified periodogram and Capon's estimates have been N,, =() offset for clarit,. The solid lines represent the exact spectrum. E[PPER(f)] = i WT(f - ) Sz(t) dt (28) TABLE II -1 2a SPECTRAL STATISTICS DERIVED FROM FITS TO MODIFIED PERIODOGRAMS, PERIODOGRAMS BASED ON PREWHITENED DATA. AND CAPON'S ESTIMATORS 1 sin rA N OF TEN SYNTHETIC 64-POINT PROFILES WITH C, = 1.0 10-6 AND WT(f N sin rf (29) = 2.4 sin / Estimator mean C, std. dev. Ca | mean i std. dev. N- I PHAN 0.724- 10-6 0.819. 10-6 2.767 0.631 PHAN(f) = NU Z Z[x]wH[n] exp (-j2 rfAn) PPW 1.310.10-6 1.690.10-6 2.469 0.567 PCAP 0.891 10-6 0.433 10-6 2.371 0.313 (30) Estimation of power-law parameters from the spectral E[PHAN(f)] = A /2 WH(f - ) ) (31) estimates is an independent problem. In this analysis, we 2 estimate c and 3 by performing a minimum absolute de- WH(f) = 1b + - + b3 (32) viation fit of a power-law function to the estimated values 4NU \ 2 2 / of the power density spectrum. The minimum absolute sin (fAN) deviation fit, which is performed in log-log space, is more b, = (33) robust than a least-squares fit. To make the power-law fit sin (irfA) independent of the frequency sampling of the spectral es- f 1 N1 timates, we evaluated the estimates at frequencies spaced (N - AJ NA very closely together so that the estimator would be a b2 = (34) smooth curve between the sampled points. While it may 1 [ 1 be possible to design a parameter estimate based on fewer i l (N- 1)Aj samples of the spectral estimate, such a method would require additional assumptions in the spectral model which sin tr f + A I N] are outside the scope of this study. (N - 1)j ) Values of the spectral estimates outside their regions of b3 = (35) high accuracy (4/N < fA < 0.5 for PHAN, 1IN <f a < sin r f+ N AJS 0.2 for PpW, and 1/(2p) _< for PCAP) were dis- (N - 1) carded as before. The mean and standard deviation of the estimated roughness amplitude cA and spectral slope / areA { Z[Xn + - Z[xn given in Table II. We see that the spectral parameters pre- N n = O dicted using Capon's estimator have roughly half the 2 variance of estimates derived from either the modified pe- exp (-j 2rfA n) (36) riodogram with the Hanning window or the periodogram based on prewhitened data. We therefore select Capon's 1L2A estimator as the preferred estimation algorithm for use E[Ppw(f)] = 2 A WT(f-) Sz(() E. Formulas for Real Frequency p-A Expressions for the spectral estimators listed in the pre- PCAP(f) eHZ e (38) ceding sections were given as functions of the normalized frequency fa that is dependent on the sampling interval of e = [1 ei'rfa ej4TfA... e~(- lafT (39)

AUSTIN etr al PROBLEMS IN ESTIMATION OF POWER-LAW SPECTRA IE seta fnoestimators, our comparison also illustrates the relative -. __ 4 \ |-S~rperformance of these estimators for spectra that have E - d hh o t variations in slope with frequency or local perturbations in level that violate monotonicity. ist tre:~ncy 10.10 III. 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. { 4 4 I 4 I ~ t 6 ~ ~ Such a spectrum might show hidden anisotropy or may *0.1 spatial requency, m' 1 indicate that the surface is isotropic in the statistical sense. In either case, the full 2-D spectrum provides valuable Fig. 9. Aliasing in the power-law spectrum case. Points marked (x) are In either case, the full 2-D spectrum provides valuable discarded to reduce errors due to aliasing. additional information. Two-dimensional spectral analyses are precluded when 2-D data are not available. Spectral studies of seafloor F. Aliasing morphology, for example, usually utilize profiles colThe estimates presented in the previous sections as- lected by bathymetric instruments that are towed by a ship, sumed that aliasing was not present, and the synthetic sur- producing a seafloor profile along the ship's path. Some face profiles were constructed to be bandlimited; i.e., they directional insight can be obtained by obtaining profiles had no spatial frequency components above the Nyquist in orthogonal directions. frequency f, = 1/(2A). Real surfaces will generally have On land, 2-D data are somewhat more accessible. frequency components at spatial frequencies above fc. As Huang and Turcotte [1], [4] and England [5] use digital a result of a necessarily insufficient sampling rate, power elevation models to estimate surface spectra, but their at frequencies higher thanf, will be aliased into the range methods average over azimuth angle in the spectral do-f <- f < fc. We can derive an equation showing this main. Two-dimensional data at scales of centimeters to result in terms of real frequency by extending the deri- tens of meters can be tedious and costly to obtain, necesvation of Oppenheim and Schafer [23, pp. 26-28], ob- sitating smaller dimensions of a sampling grid, i.e., fewer taining than 100 x 100 measured points per 2-D profile. Sz.(f) = z s (rf +. (40) A. Estimator Selection nl= -o A) Two-dimensional Fourier-based spectral estimators In the power-law case, the power spectral density falls such as the 2-D periodogram suffer from the same proboff quickly with increasing frequency. The n = 0 (exact) lems as the corresponding 1-D estimators: leakage in the and n = 1 (first alias) terms and their sum are plotted for power-law case and high variance. A 2-D Capon's estia sample power-law case in Fig. 9. We see that the effects mator is one solution to these problems. The 1-D Capon's of aliasing are most significant at the high-frequency end estimator is readily extended to the two-dimensional case of the spectral estimate. As a first-order correction to avoid [20, sect. 15.8], at the cost of increased complexity. Inthe aliasing problem, we discard values of the spectral version of the covariance matrix may become problematic estimate computed for frequencies greater than fd/2, as because the matrix is now of dimensions p2 x p2 rather indicated in the figure. than p x p. For a similar level of resolution, the 2-D Capon's estimator will require a much greater volume of G. Applicability to Real Data data (as will the other 2-D estimators). If such high quanReal geological data possess a variety of irregularities tities of data are not available, as is often the case for 2-D and surface statistics. Spectra of natural topography data sets, the variance of the Capon's estimate will be clearly cannot conform to a power-law model at all spatial correspondingly higher. frequencies-roughness features on the scale of millions of meters are unphysical, as are those at subatomic scales. B. Two-Dimensional Spectral Estimation of sotropic However, in the studies cited in Section I as well as in Surfaces our own investigations, measured spectra of topography If a surface area is known to have or can be assumed are well modeled by a power-law spectrum in the form of to have isotropic statistics (from knowledge of the origin (1) over some range of spatial frequencies. We have based of the surface, or perhaps from measured l-D spectra in our comparison of spectral estimators on an idealized several directions), then an estimate of the 2-D roughness spectrum with constant slope (except for the very-low- spectrum can be obtained from estimates of the 1-D specfrequency roll-off introduced in (17) to avoid an unphys- trum. These 1-D spectral estimates may be obtained from ical singularity) in order to have an easy reference for individual rows of the 2-D surface height grid; we can comparison. Because both the Fourier estimators with calculate a spectral estimate for several or all rows and windowing and Capon's estimator are local-neighborhood then average the estimates together. The resulting esti

q It, IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL 32 N L' mate will have a variance that is much lower than that of Let SZID(f., fV) be given by (41). and let f, be restricted a single row estimate in spite of the fact that the grid rows to the power-law region, i.e.. f, > fo,,: are not completely independent. M Suppose that an isotropic surface has a 2-D power spec- SZD(fq) = 2a [(f,) + (f)r] df. (47) tral density Sz(f, f,!) that has the form of a power law ~ down to some low spatial frequency fo,, and some un- Changing variables, we have specified finite form below that: Sz(e fr.f 0) = SZ( fr) == 2af fr' a fow SzlD(fJ) = 2a df, (48) finite fr < fow which can be integrated analytically [27, p. 201]: where the radial spatial frequency fr = (f + f2)"2. We now need to relate the parameters a and -y to the param- a - _r ( eters of the 1-D spectrum of surface heights measured 2 / along a straight line on the same surface. This 1-D spec- SID() = r( (49) trum will also have the form of a power law forf > j-owand can therefore be described by the parameters c and 1. where r (x) is the gamma function. To obtain this relation, we begin with the general Therefore, if SzlD(fx) = cf.- forf >L Jiow, the 1-D and expression for the Wiener-Khintchine relation in terms of 2-D spectral parameters are related by the Fourier transform and write the correlation E[Z(x, = + 1 (50) y)Z(x + A6, y + by)] as Rz(6x, 60): Rz (6, by) = i i_ SZ(f I fay) r c. (51) exp (j 2 r[fx x + fy by]) dfx dfy. (42) Equation (50) agrees with the result derived by Voss; he Changing to radial coordinates (6x = 6, cos 60, by = 6, sin does not give an expression comparable to (51). 60) and noting that Sz(fr, Jo) is independent of fJ, we ob- As previously discussed, 1-D profiles of a fractal surtain face with 1 < Df < 2 have power-law spectra with ex0o 2a' ponent l given by i = 5 - 2Df, so that 3 > j > 1. From RZ(6r, ae) S | (Jfr)Jfr exp (j2Jrf6, cos?) dr dr, (50), the 2-D power-law exponent'y satisfies 4 > ey > 2, o 0 corresponding to 2-D surfaces with 2 < Df < 3, where Df = (8 - y)/2. We next investigate the performance of the linear, avwhere' = eo - 60, leading to eraged Capon's estimator in the determination of the 2-D co power-law parameters a and "y. A variation of the 2-D RZ(6r) = 27 1 Sz(r)Jo(2rfr r)fr d (44) spectral synthesis algorithm described by [26] was used Jo to generate 2048 x 2048-point synthetic topographic surfaces with spectral exponents of 3.0, 3.4, and 3.6 with where JO is the zero-order Bessel function. This expres- sampling intervals A of 1 cm. The exact surface spectrum sion differs from a similar (but incorrect) expression given was modeled as isotropic, with a form similar to (17): by Voss [17, eq. (1.52)]. The 1-D and 2-D roughness spectra of an isotropic ran- (a/a2) I rl exp (-f2/[2a2]) Jr flow dom field are related by an Abel transform [2]. This re- S (f- a lation is more easily visualized using the surface autocor- f relation functions (the inverse Fourier transforms of the 0 fr > fc l-D and 2-D surface spectra). The 1-D autocorrelation function is a slice of the 2-D autocorrelation function: co cm wherefr is the radial spatial frequency, a = 2.384' 10-", SZlD(fX) = 3 3 Rz(68, by)6(6y) flow = 0.01 m-l, and ac and o are chosen as before. The I-0 I cut-off frequency fc = 1/2 A is equal to the Nyquist freexp (-j2 r[fJ6: + fyy]) d6 d6y. (45) quency along the fx and fy axes. As noted previously, arrays of data containing over a Using Fourier transform properties: million points are often unavailable. In practice, in situ rX1D~f,> Sfsurface height measurements may yield data sets that are SzlD(JX) = Sz(fJ, Jy) * 6(Ji) = 3 Sz(fJ, fJ) dfy. much smaller, say, 40 x 40 points. Since such a data set ~-a,~ ~is useful over a relatively narrow band of frequencies, (46) several sparse grids may be collected at different scales,

-t[STIN e: J PROBLEMS IN ESTIMATION OF POW'ER-LAW SPECTRA E 3.0 2....... _10 w 10 -..,.. E 2.5 Z E 5- -1. CI 1.0. — 10 U= 0.5 0.0 _.__ _ _ _ __ 0 5 10 15 20 25 30 35 distance, m. 1S.....~ PW: Eq. (56) —. PER: Eq. (55) Fig. 10. 41 x 41-point sample of a 2048 x 2048-point synthetic power- ---—..... HAN: Eq. (54) A= 1 cm law surface with a = 2.384 ~ 10- and y = 3.0. The sampling interval -104 - CAP: Eq. (53) is 50 cm.. l I I X I I Cill I I i I 0.1 1 10 100 spatial frequency, m-' E -.. A=50 cm Fig. 12. Spectral estimates calculated with (from bottom to top) PCAP, E 10 [ ----- power4aw fit PHAN, PPER (with no zero-padding), and PPw from surface profiles of a debris flow near Johnston Ridge at Mount St. Helens. (Spectral estimates c -10 - derived from techniques other than Capon's estimator have been offset for 0A = 10 an 10cm clarity.) The dashed lines are power-law fits to the estimates; the equations of these fits are (53), (54), (55), and (56), respectively. (This and the other 10.11 measured spectra from Mount St. Helens debris flows are presented in [28].).,1012 I,~........ I. I.... I I * C. Application to Measured Data 0.1, 10 Surface roughness measurements were performed on spaa frequency, m several debris flows near Mount St. Helens in support of Fig. 11. Averaged Capon's estimates from 2-D synthetic surface height o e l n data at three scales with power-law fit. a study of electromagnetic scattering by volcanic terrains. The surface profiles, which were collected using several TABLE III different sampling rates and profile lengths, were used to POWER-LAW PARAMETERS DERIVED FROM FITS TO AVERAGED CAPON'S obtain estimates of the roughness spectra of the measured ESTIMATES AT THREE SCALES debris flows using the techniques described in the present Surface (exact) Estimated parameters report. The resulting spectra are reported elsewhere [28]. a 7 mean a std. dev. a mean j std. dev. ~ One of the spectral fits calculated from measurements of 2.384.10-11 3.0 2.34. 10-11 2.83.10-'2 2.971 0.041 2.384 10-I' 3.4 2.28. 10- 1 1.55 - 10-'1 3.386 0.041 a primary debris flow near Johnston Ridge is shown in 2.384.10-" 3.6 2.17- 10-' 1.14. 10-'1 3.555 0.061 Fig. 12. This composite spectrum was derived from two surface profile grids, one with A = 1 m, the other with A = 1 cm. The equation of the power-law fit to the two a spectral estimate calculated for each, and a function fit Capon's estimates is to the composite of the individual spectra. Scale factors are particularly important in the calculation of these Sz(f) = (3.57 - 10-4)lfI-2'31 (53) multiscale spectra. Spectral estimates based on PHAN, PPER (with no zero-padIn our investigation, we sample each 2048 x 2048- ding), and Ppw were calculated from the same data for point synthetic surface at sampling intervals of 50, 10, comparison with (53) and are shown in Fig. 12. Equations and 4 cm, obtaining three 41 x 41-point data grids. (A of the power-law fits derived from these estimators are sample data grid is shown in Fig. 10.) At each scale, a -2.36 l-D Capon's estimate (with p = 12) is calculated for each Z, HAN(f) = (6.28 10-4)fl (4) grid row at frequencies 1/(2pA) < f < 1/(4A), and the SZ PER(f) = (7.50 ~ 10-4) lf-2.36 (55) estimates are averaged over the rows. A power-law function is then fit to the averaged spectral estimates at the Sz,pw(f) = (4.20 ~ 10-4)I fl-240 (56) three scales. (A sample fit is shown in Fig. 11.) The upper The spectral exponents in (53)-(56) seem acceptably simfrequency limit of 1/(4A) was chosen to reduce the effect ilar, but this is principally due to the combination of specof aliasing. Estimates of 1-D power-law parameters c and tral estimates at different scales into a single composite /5 are then converted into estimates of a and y using (50) spectrum. If we fit power-law functions to spectral estimates based on the A = 1 cm data alone, we obtain the This process was performed on 10 synthetic surfaces at follo each of three spectral exponents. The mean and standard wing deviation of a and fj are listed in Table III. ]CAP(f) = (1.15. 10-3) fl-273 (57)

Usa IEEE TRA.NSACrTI(IUNS ON 3Et()CiIE\jt. A\I[) Ri-l) e ti —Ul\u T\AL so'4t L L'i I14: PHAN (f) = (3.95 ~ 10-3) f -3 (58) [7] T. H Bell. Jr. Statistical features of sea-floor topograph, " DeepSea Res. vol. 22. pp. 883-892. 1975 PPER(f) = (2. 12 10 4) Ifl-1.85 81(59) 8 -. "Mesoscale sea floor roughness." Deep-Sea Res.. vol 26A. pp. 65-76. 1979. P (992 10) -75 [9] C. G. Fox and D. E Haves. "Quantitative methods for analszing the ~Ppw(f) -= (9.92 ~ 10-) I]f-, (60) roughness ot the seafloor." Rev. Geoph.s.. sol. 23. pp. 1-48. Feb he Capon's, modified peiodogram with Hannin win- 1985. Fhe Capon's, modified periodogram with Hanning win- [1101 L. E. Gilbert and A. Nialinvemo. "A characterization of the spectral low, and prewhitened spectral estimators give spectral density ot residual ocean floor topography." Geophys. Res. Lett..;opes between 2.73 and 3.0. but the periodogram (with sol. 15. pp 1401-14(. No. 1988I [1 11 A. Malinvemo. "Segmentation of topographic profiles of the seafloor io zero-padding) gives a significantly different slope i3 = based on a selt-affine model." IEEE J. Ocean Engineer.. vol. 14..85, illustrating the effects of spectral leakage. pp. 348-359. Oct. 1989 1121 R. S. Savles and T. R. Thomas. "Surface topography as a nonstationars random process.'' Nature. vol. 271. pp. 431-434, Feb. 2. IV. CONCLUSIONS 1978. [131 M. V. Berm and J. H. Hannay. "Topography of random surfaces." Spectral analysis is a tool of increasing importance in Nature, vol. 273. p. 573. June 15. 1978. he characterization of topography and other natural sur- 1141 B. B. Mandelbrot. The Fractal Geometry of Nature. New York: Freeman. 1983. aces. Numerous studies indicate that such surfaces have [151 R. J. Adler. The Geometry of Random Fields. New York: Wilev.,ower-law spectra in the form of (1) over some range of 1981. patial frequencies. This spectral form introduces unique [16] B. B. Mandelbrot and J. W. Van Ness. "Fractional Brownian motion process.s fractional noises, and applications," SIAM Rev., vol. 10, pp. ifficulties in the spectral estimation process. 422-437. 1968. In the present work, we have shown how leakage causes [17] R. F. Voss. "Fractals in nature: From characterization to simulaome Fourier-based estimators to yield an estimated spec- tion," in The Science of Fractal Images, H.-O. Peitgen and D. Saupe, Eds. New York: Springer-Verlag, 1988, pp. 21-70. 1al exponent of 2.0 for surface profiles having spectral [18] S. E. Hough, "On the use of spectral methods for the determination xponents between 2 and 3. This may explain the number of fractal dimension." Geophys. Res. Lett.. vol. 16, pp. 673-676, f studies reporting power-law exponents of 2.0 for nat- July 1989.,studiesgowe aw eponet o2.0for [19] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, ra~l terrains. We show how Capon's estimator may be Numerical Recipes (FORTRAN Version). Cambridge: Cambridgesed to avoid the leakage problem and measure the sur- Univ., 1989. ice spectrum more accurately. We also show that [20] S. M. Kay. Modern Spectral Estimation. Englewood Cliffs, NJ: I,~~~~~~~~~~~~ ~~~~Prentice Hall, 1988. lapon's estimator has reduced variance, which is useful 121] P. D. Welch, "The use of fast Fourier transform for the estimation,hen surface data records are short, as is often the case of power spectra: A method based on time averaging over short, modI remote-sensing studies. ified periodograms," IEEE Trans. Audio Electroacoust., vol. AUThe two-dimensional Capon's estimator may be em- 15, pp. 70-73, June 1967. [22] R. B. Blackman and J. W. Tukey, The Measurement of Power Specloyed to compute 2-D spectra of rough surfaces. How- tra from the Point of View of Communications Engineering. New ver, two-dimensional data are often unavailable or too York: Dover, 1958. )arse. For surfaces that are known or can be assumed to [23] A. V. Oppenheim and R. W. Schafer, Digital Signal Processing. >arse. For surfaces that are known or can be assumed to Englewood Cliffs, NJ: Prentice Hall, 1975. Englewood Cliffs, NJ: Prentice Hall, 1975. > isotropic, we show how a linear, averaged Capon's [24] R. E. Ziemer and W. H. Tranter, Principles of Communications: Sysitimator can be used to estimate the 2-D power-law pa- tems, Modulation, and Noise, 2nd ed. Boston, MA: Houghton Mifflin, 1985. tmeters. [25] J. Capon, "High-resolution frequency-wavenumber spectrum analysis," Proc. IEEE, vol. 57, pp. 1408-1418, Aug. 1969. [26] D. Saupe, "Algorithms for random fractals," in The Science of FracACKNOWLEDGMENT tal Images, H.-O. Peitgen and D. Saupe, Eds. New York: SpringerThe authors wish to thank several anonymous review- Verlag, 1988, pp. 71-136. [27] Bateman Manuscript Project, Tables of Integral Transfornns, vol. 2,'s for their helpful comments and suggestions. A. Erddlyi, Ed. New York: McGraw-Hill, 1954. [28] R. T. Austin and A. W. England, "Multi-scale roughness spectra of Mount St. Helens debris flows," Geophys. Res. Lett., vol. 20, pp. REFERENCES 1603-1606, Aug. 6, 1993. 1] J. Huang and D. L. Turcotte, "Fractal mapping of digitized images: Application to the topography of Arizona and comparisons with synthetic images," J. Geophys. Res., vol. 94, pp. 7491-7495, June 10, 1989. Richard T. Austin (S'84) was born in Maryville, 2] i. A. Goff, "Comment on'Fractal mapping of digitized images: Ap- TN, on September 2, 1964. He received the B.S. plication to the topography of Arizona and comparison with synthetic degree in electrical engineering from the Univerimages,"' J. Geophys. Res., vol. 95, p. 5159, Apr. 10, 1990. sity of Kentucky, Lexington, in 1986 and the 3] J. Huang and D. L. Turcotte, "Reply," J. Geophys. Res., vol. 95, M.S.E. degree in electrical engineering from the p. 5161, Apr. 10, 1990. University of Michigan, Ann Arbor, in 1987. t] -, "Fractal image analysis: Application to the topography of He is presently pursuing the Ph.D. degree in Oregon and synthetic images," J. Opt. Soc. Amer.: A, vol. 7, pp. electrical engineering at the Radiation Laboratory 1124-1130, June 1990.. at the University of Michigan. During his studies, 5] A. W. England, "The fractal dimension of diverse topographies and - he has held a National Science Foundation Gradthe effect of spatial windowing," Geolog. Surv. Can., Paper 904: uate Fellowship and a Michigan Space Grant FelGround Penetrating Radar, J. A. Pilon, Ed., pp. 57-61, 1992. lowship. His research interests including scattering and emission from nati] S. R. Brown and C. H. Scholz, "Broad bandwidth study of the to- ural surfaces and radiometric studies of the atmosphere. pography of natural rock surfaces," J. Geophys. Res., vol. 90, pp. Mr. Austin is a member of Tau Beta Pi and the American Geophysical 12 575-12 582, 1985. Union.

-ALSTI\. PRO(BLEMIS IN ESTIM,T~1()N, OF POW'ER-L.AV SPECTRA Anthon, W'. England (M'87-SM'89) received Gregory H. Wakefield (S'84-M'84) receled the the B.S. degree in earth sciences (1965). the M.S. B.A. degree summa cum laude in mathematics and degree In geology and geophysics (1965), and the psychology. the M.S. degree in electrical engiPh.D. degree In geophysics (1970) from the Mas- neering. the Ph.D degree in electrical engineersachusetts Institute of Technology. Cambridge, ing. and the Ph.D degree in psycholog,. all from M A. the Universit, of Minnesota. Minneapolis. in -- i _His research Interests have included terrestrial 1978. 1982. 1985. and 1988. respectivel, heat flow; geomagnetic and gravimetric studies in From 1978 to 1985. he was with the Pschothe northern Rockies and in Antarctica: radar acoustics Laborator'. IUniversitv of Minnesota. studies of temperate and polar glac:. ->. and the From 1983 to 1985. he was with the Communimicrokwave radiometric signatures ot snow, ice, cations Laborator\ and the St. Anthony Falls H\frozen and thawed soils. and planetarD regoliths. He was Deputy Chief of draulics Laboratory. Universits of Minnesota. From 1982 to 1983. he was the Ottice of Geochemistrn and Geophysics of the U.S. Geological Survey, a consultant to the Technology Strategy Center. Honeywell. Minneapolis. and was an Associate Editor for the Journal of Geophysical Research. He From 1982 to 1985. he was consultant to the Bioscience Laborator. Life served on several federal committees concerned with Antarctic policy, nu- Science Sector, 3M. St. Paul. Since 1986. he has been on the faculty of clear waste containment. and federal science and technology. He was ap- the University of Michigan. Ann Arbor. where he is currently an Associate pointed as a NASA Scientist-Astronaut in 1967. acted as Mission Scientist Professor with the Department of Electrical Engineering and Computer for Apollo's 13 and 16. and fle`w as a Mission Specialist on Space Shuttle Science. His research interests are in statistical signal processing. signal Challenger's Spacelab 2 Mission in 1985-a solar astronomy and plasma analysis and synthesis, image processing, spectral estimation, sensory sysphysics mission. He has been a Program Scientist for NASA's Space Sta- tems modeling, psychoacoustics. sensory prosthetics, speech processing, tion. He is currently a member of the NAS/NRC's Space Studies Board and music processing. He has been involved in teaching and research in and a member of several NASA advisory committees. He has been an Ad- these areas and has led in the development of the educational laboratories junct Professor at William Marsh Rice University, Houston, TX, and is in signal processing (SPEL) and in image processing (IPEL) at the Uninow a Professor of Electrical Engineering and Computer Science at the versity of Michigan. University of Michigan., Ann Arbor, MI. Dr. Wakefield received the Presidential Young Investigator Award Dr. England is a member of the American Geophysical Union. (1987), the University of Minnesota Doctoral Dissertation Fellowship (1982). and the USPHS Traineeship with the Center for Research in Human Learning (1979-1982). He is a member of Phi Beta Kappa and Sigma Xi and his professional affiliations include the Acoustical Society of America.

22 V Empirical Model There are no reliable scattering models for surfaces whose roughness spans the spatial frequencies of the radar's wavenumber. In fact, a comparison between an empirical model and the predictions of several candidate theories is one of the tasks of the study. The empirical model was based upon scattering measurements of milled surfaces that were carefully produced to exhibit power-law roughness for 3 values of f (2.29, 2.55, and 2.80). The process for manufacturing the scattering surfaces became Chapter V in R.T. Austin's Ph.D. dissertation [69]. The process of measuring the radar backscatter from these surfaces became Chapter VI of the dissertation. His results were presented at an International Geoscience and Remote Sensing Seminar (IGARSS'94) in Pasadena, CA [68]. They will also appear in a TGARS paper. Chapters V and VI of Dr. Austin's thesis appears on the following 41 pages.

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 difflicult 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 measu ment 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 subrface inhomogeneities is usually present. Both of these factors yield eects that are difficult to decouple from those due to urace rghna. J —,mst, intrpretation of scattering measurments is hindered by the dific ly of obtinin 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-laws 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 manufactured 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 [441, 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 f to avoid volume scattering) * high permittivity i 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 155]. 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 ZA 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 6inch and one 2-inch) and the radar bandwidth set to 3 GBz 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.333 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: c, = 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(f,) = afr, (5.2) where f,. = (f2 + f2)1/ is the radial spatial frequency. First, a set of 2048 x 2048 discrete Fourier amplitudes Zl(f3, f,) is generated such that al, 1A22,(ft, f)l12 = Sz(f,,f,) (5.3) for all points within the circle ft < f, where f., the maximum spatial frequency, was (2)-1 = 406.5 m1'. Point spacing in the frequency domain was = (N )1 =

86 0.39698 m-'. In other words, -Z (ff i)'= U(f ) (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(fs, fA) = Zl(fl, fA)ge', (5.5) where g is a Gaussian random variable with zero mean and unity variance, and 4 is a phase random variable uniformly distributed over [0, 2r]. 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(f, f,) = Z'(-f,, -f,), (5.6) The randomized coefficien s Z(f5, fs,4) 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: 1 N-1 N-1 Z(Xk, yl) = N2 Z2(f., f,,f) exp(-i2r[km + ln]/N). (5.7) The resultant 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 a in (5.2), were selected to be similar to those measured at Mount St. Helens and to be physically realizable n UHMW slabs using the milling machine. At the time the first artificial surface was designed,

30 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; = 3.29. (The surface spectral estimate was subsequently changed to a = 1.7672 x 10-4 and' = 3.345.) The UHMN 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): J (a/a'2)fr exp(-fr/[2a2]), f, < 0.01 m1Sz(f,) = af,?', 0.01 m-1 f <_ f, (5.8) 0, fr _ fc where a = 2.17143 x 10-5 = 3.29 a = 1.6438 a= 0.00482805 = 406.5 mThe 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 10Y 1 3.29 3.8 3.545 i a 1.6438 19.8523 5.70356 I a f 0.00482805 0.00456435 0.00469065 Df 2.355 2.1 2.2275 Table 5.1: Spectral parameters (a, r, ca, 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, y, 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 thremaxis gantry positioning table, a Kress variable-speed electronic grinder, and a Techao MAC 100 programmable position controller, all purchased from Tedhno, 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 serial port drive E motor z drive -tz-axis motor 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 zyz coordinates. The table and z- 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 minll a irl thick slab measuring one meter square in a single piece. One problem that became apparent soon after starting work with the table w that the table sagged somewhat under its own weight due to support bears that were

90 33 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 Masterm Mill package is an extensive program for desinin complex objects, producing machine toolpaths for the fabrication of these objects, and controlling the machining equipment while executing the machine toolpath instructions. (However, the Tchno controller must be operated by a TeDhno postprocessor program and not from ittin M tca) Becue the srycs in hi study were extremely ar and Cople, adbs. they we POdd by p external program on another computing plathem, many of the Itue EdAbCHai

91 tercam package were not used. The principal use of the Mastercamrn 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-aiis 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 1l.5inch tension spring and support rod

92 11.5 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 suficiently 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 I 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/4inch endmill had zi tIr ogh the slab and through part of the table itself. The simultanous failures of two uas led

93 3) PDS 486/33 Kress Electronic Techno computer Grinder MAC 100 AT6400 indexer ~18, 26 VAC (in slot) transformer AT6400 Compumotor SC30 3-axis milling table connector box Drive RaPd x drive x-axis motor encoder y drive y-axis motor encoder z Av z-axis motor encoder 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 Comnpumotor, was used with Microsoft Visual Basic to write new control sote. 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 z 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 3SS 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 ds2amcinv.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 zyz 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 edguaz.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. Mastram is then used to convert the surface elevation files into toolpath fl at two resolutions. Roughing files contain toolpaths for a 1/4inch ball-end eudmill. This tool is used to

39 96 reference surface 1/4 inch radius Figure 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. 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 I T 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 z 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 millimtrs in places where the surface was low near its edge. A Visual Basic program called aill 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. Rofghing 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 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, aillscan, 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 it 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. 1000 _m 800 600 400.. E c 200 -0.3 -0.2 -0.1 0.0 0.1 0.2 surface error (inchu) Figure 5.8: A histogram of differences between the specified surface elevations and the scanned surface height.

104 of the endmill position on the metal block position to within 3-6 steps (0.03-0.06 mm) after numerous milling seions 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 postprocessing 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 GIIz. 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 mixn with the signal from a second local oscillator, resulting in a 24 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.80 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 RECEIVER MODULE IF OUT V 2-4 GHz IF OUT H l 2-4F UP 24 o n 9wI Figure 6.1: Block diagram of the 35 GHz radar. 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~ 9 9 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 9 9~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ S~ ~ ~ ~ ~ ~ ~ l~~T! MOUL 9~~~~Fgr 6. " B o k dar m o Sh 5Gzrdr

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 6inch 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, G = 2.16 - j0.0022 at 10 GHz [29]) was used to fill the air gaps. Plots of

109 Gateway Polarization switch control computer servo control 0 servo control I transmit receive HP-IB HP 8720 IF UP Network Analyzer H l HP 3488A | Switch/Control Unit Aerotech Unidex 11 Turntable Controller Figure 6.2: Components of the radar msurement ystem.

110 radar absorber wall 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.7w5-inch sphere was about 38 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, T1, and T2: * 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 cl, C2, C3, r1, 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 — 27.6 _-2,7.6 -. Theoretical E.......... -W measured -278 8' - -...- HH measured.27.8. 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). o31 4 -31.X 3.1.. -32.2! U 3U I'' I' 1 " ~' 33.6 34.2 34.8 35.4 36.0 fretquency (GHz) Figure 6.5: Test measurmt of 1.188-inch diameter spere. Figure 6.5: Test measuremet of 1.188.tach 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 1800 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 150 to 600 at 5~ increments. The measurement sequence at each incidence angle consisted of the following steps: * with radar lowered, set incidence angle, 8 * 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 transformng a fequency-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 transrming the result back to the frequency domain. The network analyzer was set to its standard gate shape; the gate width was st to 5 us in order to include reflectio froqp 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 250, 200, 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 at(fi) = (2VfiSrt(fi)l2), (6.1) where r and t indicate the receive and transmit polarizations and may be either V or H, and St 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 a0(fi) are obtained by applying an illumination integral correction to a(fi). 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,(e) = P G(G), (6.2)

117 and the distributed target form is P, () = (4 3 I(h. O)ao~(), (6.3) where I(h. 0) is called the illumination integral: I(h, 9) = | gt(Ot, I'O)gr) R4(O I) dA. (6.4) 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 P,(O)/Pt equal in (6.2) and (6.3) and solving for a0 in terms of a: ao~(fi) = 1R4l(h ) a (6.5) The program measillum.f was used to convert values of a(f,) to values of ao0(fi) 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:' = lip, (6.6)

118 a' = pa, (conductivity) (6.7) f' = pf, (6.8) A' = /p, (6.9)' = {, (6.10) is = /s, (6.11) a' = alp2, (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, ao, is obtained by noting that the surface area, A, scales as length squared: 1 2 A A' = (11)2= (. (6.13) Writing the scattering coefficient as the average radar cross section normalized by area, we obtain () = (') = (e/p2) = (a) = a (6.14) A' A/p2 A a Because ao 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 50 between measurements. A few measurement runs were executed using smaller angular increments, Ad, 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.10 for 600 incidence and 4.50 for 150 incidence (based on 60-sample traces with AO = 0.5~). For UHMW1, the decorrelation angle is approximately 1.50 for 600 incidence (based on a 90-sample trace with AO = 2'). A data set with fine angular increments was not taken on UHMW1 at 150 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 N., therefore varied from 300'/7.5' = 40 (estimate for UHMW1 at 150) to 60 (the nwmber of samples collected) for UHMW1 at 60' and UHMW3 at all incidence angles. Values of a 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 a N., 1 fNvg Ntonal UHMW1 15' 40 (est.) 1.3 52 600 60 5.1 306 UHMW3 150 60 1.3 78 60~ 60 5.1 306 Table 6.1: Estimated number of independent samples B: B B sin a cr~~ 2 (6.15) 2fB a si 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, Nf1,9 can be shown to vary from 1.3 at 150 incidence to 5.1 at 600 incidence. The total number of independent samples, Ntota, is calculated by multiplying the numbers gained by spatial sampling and frequency averaging: Ntotal = NoNfv,,, (6.16) The resulting numbers are listed in Table 6.1. Values of Ntota, 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.

64 VI A Comparison with Various Scattering Models There are no theoretical scattering models for random surfaces which have significant surface roughness at spatial scales that are near to the wavenumber of the scattered signal. The Phased Wiener-Hermite (PWH) model had never been tested and there was reason to believe that it might yield better approximations to scattering by power-law surfaces. In fact, the PWH model did not do well. Of the common theories that were tested under this investigation, the Small Perturbation Model performed best. The comparison among theories and experimental scattering results became Chapter VII in R.T. Austin's Ph.D. dissertation [69]. These results will also appear in a TGARS paper. Chapter VII of Dr. Austin's thesis appears on the following 14 pages.

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 ao resulting from several rough surface scattering models. 7.1 Experimental Results Measured values of uw and aO 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 & values for both polarizations at 5.5~ for the two smoother surfaces. There is little difference between the O, and aOu for a given surface —a, 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 o ~ Surmce A (UHMW1) A A Surface B (UHMW3) t O U* Sualmce C (UHMW2) -1 HE t8 o i -15-' 0 f-20- r o I D e~A A ~0 D A -25 A A IA A Io ~ -30-., I I I I 20 30 40 50 60 incidence angle (degrees) Figure 7.1: Measured co-polarized backscattering coefficients, ao and aOh, 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, ao and acv are shown in Figure 7.2. The measured aoh, and ~ho 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 a~h and ah,, 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 0ao and Ar. 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 - 25 0 o 0 Surface A (UHMW1) A A Surface B (UHMW3) o 0 U Surface C (UHMW2) HV VH -30mQ 2.45 -40'L'i * ii 20 30 40 50 60 incidence angle (degrees) Figure 7.2: Measured cross-polarized backscattering coefficients, ao' and oa,, of the three artificial surfaces: Surface A (UHMW1), D, = 2.355; Surface B (UHMW3), D, = 2.2275; and Surface C (UHMW2), D1 = 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 sifce 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, I 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] 00 — = c0 +-0 + ao (7.4) Pq pqc pqn P9a 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 k2 IR.ol2 | 2 COS 4k2cos2 9ff2 E' ( -4k= c s2 0f2)" Iz)(-2k sin 0, 0) (7.5) 0n = o (7.6) where Rlo and R110 are the Fresnel reflection coefficients for horizontal and vertical polarization, respectively, and nc00 00 IF)(q,l)= s (RZ(wscu ))l i(qu+qmV) dd (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 1I().

125 The slope term a, was omitted for the power-law case because the non-analytical correlation function prevented the direct evaluation of the slope terms. The a0qs is much smaller than on in the Gaussian case, so it is assumed that the term is of little consequence. The physical optics results are compared with the measured ao 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 450 or so. The VV predictions start about 5 dB low at 150, 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 fic = (10A)1', where A is the radar wavelength. The low-frequency cutoff was imposed due to programming constraints in the calculation of the I(") terms. Letting fic vary from (5A)'1 to (25A)-1 resulted in a change in a' 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 [591. 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

70 126 O Surface A (UHMW1) 1 I 0 -_ AA Surface B (UHMW3) o O O Surface C (UHMW2) -e- Phys. Optcs, UHMW1 D-15-~o -.U- Phys. Optics, UHMW2 0 _ -20 LD 0 ~-25Q= Q E 0 -35 20 30 40 50 60 incidence angle (degrees) (a) ~-10-5. * C "CW H 3 0 0 3~~-25~~~ su A (UHMW) 5 0 Surface B (UHIW3) - 0 o- 0 Surface C (UHMW2) 0 -e- Phys. Opt*, UHMW1 = I ~c- Phys. Opbcs, UHMW2 -20- - 3 cm -25- AA -30-35-40 vv 20 30 40 50 60 incidence angle (degrees) (b) Figure 7.3: Physical optics results for (a) HH and (b) VV backscatter compared to the measured values for the three artificial surfaces.

7i 127 criteria for this solution are [56] ki > 6 (7.8) R, > A (7.9) koa > (7.10) cos 0, - cos 0, where a is the standard deviation of surface height, Oi and 9f are the incident and scattered angles, and k, 1, R,, and A are defined as before. The geometrical optics solution is given by Ulaby et al. [59]: o _ IR(O)12 dPP =2cos42Ip"(0)I (7.11) 2 cos' -'I0'p"(O)l 0_O (7.12) where R is the Fresnel reflection coefficient and p"(O) is the second derivative of the normalized correlation function, p(r), Rz(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): r - ro 1P2 - PI P - Po A A -- (p2 - 2Pi + po) (7.14) w here the first three correlation values po, p, and Ph are numerically caculated using the same bandlimited spectrum used for the physical optics case. For the three

128 ""-'. 0o Surface A (UHMW1).-,1ol\ A. IA Surfce B (UHMW3) O Surface C (UHMW2) ~ -1 5. o E-205 2 0 _-\ " "_ A ~ 4D \ -30 - -35 -40! I I I 20 30 40 50 60 incidence angle (degrees) Figure 7.4: Geometrical optics results for backscattering compared to the measured values for the three artificial surfaces. artificial surfaces, A = 1.23 mm, and estimates of the rms slope, m = JI (0), 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 kc < 0.3, (7.15)

73 129 m < 0.3, (7.16) where a. the standard deviation of surface elevation, is computed "from frequency components of the surface responsible for scattering at a given electromagnetic wavelength" [391. The small perturbation model is usually used for larger incidence angles (9. > 200) [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 o =- 8k4 a cos4 la1q l2W(2k sin e, 0), (7.17) where ahh = R1 (), sin2 8 - 4(1 + sin2 0) )[e, cos 0 + ((, - sin2 0)1/212 Clh = Ckhv =, and W(k., kv) is the normalized roughness spectrum, 1 co coo W(k:, k) = 2 f / p(u, v) exp[-j(kzu + kv)]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(rT, T,), is defined as Rz(T,) )= E Sz(fg,I f) p2r(fr, + frT)]dfdf,. (7.19) The normalized spectrum W(k, k,) with k. and k, in rad/m is therefore W(~, i,)=2= 2sz 2 w',2' (7.20)

74 130 and we see that W(2k sin 8,O) = r ksinz o (7.21) and therefore 0a = 8 OS48eaq[21 [ksin 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 c were tried. Changing, from 2.05 to 2.3 raised a~ by about 1 dB. Errors in the c, measurement of more than 0.3 are considered extremely unlikely. Adding an imaginary component of, = 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 0 O Surface A (UHMW1) -10 A B Surface B (UHMW3) o0 Surface C (UHMW2) o - -5\ I-O |............... SPM -35 - 0 ~ -. -... 0 G) ca-350-11 -40.... 20 30 40 50 60 incidence angle (degrees) (a) 0 Surhm A (UHMWI) G ~ope o tSurm C('UW2)I c~-20 ~~P cm 0 j-25 ~ ~ -30 -35 VY 20 30 40 50 60 incidence angle (degrees) (b) Fi'gure7.5: Small perturbation model results for (a) HR and (b) VV ba catte compared to the vaues.

76 132 current in a series of orthogonal Wiener-Hernite 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, C(s,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 ao is rather lengthy and involves a number of intermediate terms, and is therefore left in the Appendix. The PWH solutions for ao of the three surface analogues are shown for HH and VV polarizations in Figure 7.6. The modeled auh, which were calculated using the bandlimited roughness spectrum truncated at fic = (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 o values that were closer together in magnitude, but the values overestimated the measured values worse than in the narrower spectrum case. Values of a~ were closer, differing from the measured

133 values by 6 dB or less. The trends of oa vs. 0 are closer than those of Aoh, but still not as close as those predicted by the small perturbation method. W'hile 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 0~~~~~O Surface A (UHMW1) O- A Surface B (UHMW3) t ~'"~ ~o0 Surface C (UHMW2) ~O ~-~ ~~~ —- PWH, UHMW1 -.-II- PWH, UHMW2 co -10 0 A 3.) 0 A O A o A S -20 0 o O 0 O =0 20 30 40 50 60 incidence angle (degrees) (b) Surface A (UHMW1)pared to the measured. A Surface B (UHMW3) E 3 Surface C (UHMW2) - I 0' PWH, UHMWI O I- PWH, UHMW2._ -20- 0 20 30 4O 5O 60 incidence angle (degrees) (b) Figure 7.6' Phased Wiener-Hermite model results for (a) HH and (b) VV backscatvalues.

79 VII Continuing Work Empirical laws or, to a limited extent, small perturbation theory can be used to relate D and c in equation (1) to scattering cross section. We are currently testing this approach by developing a classification algorithm from the empirical laws discussed in Section V. Using this algorithm, we will attempt to classify Andean Volcanics according to their roughness. This application is funded by a SIR-C project and was not a part of our Mars investigation. VIII Acknowledgement Dr. Austin and I wish to thank the NASA Planetary Science Program for its support of our work. While the issues are not exhausted, we believe that we have made significant progress toward a quantitative examination of the diffuse component of Mars radar backscatter. IX References. [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. Wakefield, "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 seal-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.

80 [12] R.H. Cameron and W.T. Martin, "The Orthogonal Development of Nonlinear Functionals in Series of Fourier-Hermite Functionals," Annals of Mathematics, vol. 48, pp. 385-392, 1947. [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 Wave-Material Interaction, vol. 2, pp. 9-25, 1987. [16] Cornel Eftimiu, "Scattering by a rough dielectric interface: A modified Wiener-Hermite expansion approach," Journal of the Optical Society of American A, vol. 7, pp. 875-884, May 1990. [17] Cornel Eftimiu, "First-order Wiener-Hermite expansion in the electromagnetic scattering by conduction 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 1980 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. [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.

81 [28] J.A. Goff, "Comment of'Fractal Mapping of Digital 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 American: 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: W.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, "Reflections 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-161, 1987. [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. 10491057, 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

82 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-Receiver 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. 361378, 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. [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, "Fluctuation Statistics of Millimeter-Wave Scattering From Distributed Targets," IEEE Transactions on Geoscience and Remote Sensing, vol. 26, pp. 268-181, 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 stimulation," in The Science of Fractal Images, H.-O. Peitgen and D. Saupe, Eds. New York: Springer-Verlag, 1988, pp. 21-70.

83 [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. [65] Austin, R.T., and A.W. England, Radar scattering from volcanic debris flows, Progress in Electromagnetic Research Symposium, 1991, Cambridge, MA, July 1-5, 1991. [66] Austin, R.T., and A.W. England, Surface characterization of volcanic debris flows at multiple scales, Proc. of IGARSS'91, pp. 1675-1678, Espoo, Finland, June 3-6, 1991. [671 Austin, R.T., and A.W. England, Multi-scale roughness spectra of volcanic debris flows, Lunar and Planetary Science Conf. XXIII, Houston, Texas, March 16-20, 1992. [68] Austin, R.T., and A.W. England, Radar scattering by volcanic debris flow surface analogues, Proc. of IGARSS'94, Pasadena, CA, August 8-12, 1994. [69] Austin, R.T., Electromagnetic Wave Scattering by Power-Law Surfaces, a dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Electrical Engineering at The University of Michigan, 1994.