Lk( cJ- j C Fu, 0e. ma-;t 7-t~ - -56 ts ( Ic MILLIMETER WAVE MEASUREMENT AND MODELING OF TERRAIN SCATTERING FINAL REPORT U.S. ARMY RESEARCH OFFICE CONTRACT DAAL03-89-K-0056 February, 1992 THE VIEW, OPINIONS, AND/OR FINDINGS CONTAINED IN THIS REPORT ARE THOSE OF THE AUTHOR(S) AND SHOULD NOT BE CONSTRUED AS AN OFFICIAL DEPARTMENT OF THE ARMY POSITION, POLICY, OR DECISION, UNLESS SO DESIGNATED BY OTHER DOCUMENTATION.

F.Qi2 V$Y

MASTER COPY KEEP THIS COPY FOR REPRODUCTION PURPOSES i m REPORT DOCUMENTATION PAGE Form Approved OMB No. 0704-0188 Public reporting burden for this collection of information is estmated to average 1 hour e rsponse, icluding the time for reviewing instructions. searching existlng data sources. gathernng and maintaining the data needed. and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of nformation, including suggestions for reducing thi burden. to Washington Headquarters Services. Directorate for Information Operations and Reports. 1215 Jefferson Davis Highway. Suite 1204, Arlington. VA 222024302. and to the Office of Management and Budget. Paperwork Reduction Project (0704-0188). Washington, DC 20503. 1. AGENCY USE ONLY (Leave blank). REPORT DATE 3. REPORT TYPE AND DATES COVERED 1I92, February 15 Final 2-1-89 - 1-31-92 4. TITLE AND SUBTITLE S. FUNDING NUMBERS MILLIMETER WAVE MEASUREMENT AND MODELING OF TERRAIN SCATTERING 6. AUTHOR(S) Fawwaz T. Ulaby 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION REPORT NUMBER Radiation Laboratory, University of Michigan Ann Arbor, Michigan 48109 026247-1-F 9. SPONSORING /MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSORING /MONITORING AGENCY REPORT NUMBER U. S. Army Research Office P. O. Box 12211 DAAL03-89-K-0056 Research Triangle Park, NC 27709-2211 11. SUPPLEMENTARY NOTES The view, opinions and/or findings contained in this report are those of the author(s) and should not be construed as an official Department of the Army osition, policy, or decision, unless so designated by other documentation. 12a. DISTRIBUTION/ AVAILABILITY STATEMENT 12b. DISTRIBUTION CODE Approved for public release; distribution unlimited. 13. ABSTRACT (Maximum 200 words) This final report provides a summary of the results realized over the past three years with regard to improved characterization of radar backscatter from terrain. The papers and reports published under this program cover three types of contributions: (a) techniques for measuring the polarimetric scattering response of distributed targets using coherent-onreceive systems, (b) theoretical models for polarimetric scattering from tree canopies and snow-covered terrain, and (c) extensive experimental measurements at 35, 94, and 140 GHz for model targets (such as small trees) under laboratory conditions as measurements for natural surfaces from truck-mounted platforms. 14. SUBJECT TERMS 15. NUMBER OF PAGES Millimeter waves, radar scattering, terrain clutter 104 16. PRICE CODE 17. SECUR1TY CLASSIFICATION 18. SECURITY CLASSIFICATION 19. SECURITY CLASSIFICATION 20. LIMITATION OF ABSTRAC OF REPORT, OF THIS PAGE OF ABSTRACT UNCLASSIFIED UNCLASSIFIED UNCLASSIFIED UL NSN /754U 1-280U-5500 Standard Form 298 (Rev. 2-89) Prescriel- by ANSI Std Z39-18

Contents 1 INTRODUCTION 2 2 SUMMARY OF RESULTS 2.1 MM-Wave Polarimetric Radar Systems..... 2.2 Modeling Polarimetric Backscatter From Terrain 2.2.1 Tree Canopies................ 2.2.2 Snow-Covered Terrain........ 2.2.3 Polarization Phase Statistics.... 2.3 Millimeter-Wave Radar Data Base....... 2 ~.. ~.. 2 ~..... 3...... 4 *.. - 0 e 5 3 LIST OF PUBLICATIONS 4 PARTICIPATING SCIENTIFIC PERSONNEL 5 CONCLUSIONS APPENDIX A: REPRINTS OF SELECTED PAPERS 5 9 10 11 1

1 INTRODUCTION An intensive program, comprised of modeling efforts and intensive experimental measurements, was conducted over the past three years aimed at improving our understanding of the radar response of terrain at millimeterwave frequencies. The experimental measurements, which were conducted primarily at 35, 94, and 140 GHz, included laboratory investigations of specific targets-such as small trees-as well as field observations of trees and snow-covered terrain. Of particular interest was to model and characterize the polarimetric radar backscatter of statistically homogeneous targets. To this end, polarimetric radar systems were constructed and appropriate calibration techniques were developed. This report provides a summary of the results realized under this program, with the details provided in Appendix A in the form of reprints of articles published in scientific journals and symposia proceedings. Not included in the Appendix are two long reports [24, 25] documenting MM-Wave radar observations acquired under this program, as well as data reported by other institutions. These reports have been requested by and provided to numerous U.S. Army and other DOD groups and laboratories. The references cited in the body of this report refer to the publication list given in Section 3. 2 SUMMARY OF RESULTS 2.1 MM-Wave Polarimetric Radar Systems Under a previous ARO contract, we had developed radar scatterometer systems at 35, 94, and 140 GHz. In order to investigate the polarimetric response of terrain, we needed to modify the system and develop appropriate calibration techniques to allow us to make accurate measurements of the 2

backscattered signal in a polarimetric mode. For indoor measurements in an anechoic chamber, the fully coherent polarimetric technique can be used to measure the scattering matrix C the target, from which both the magnitude and phase statistics can bc Jbtained directly. Under outdoor conditions, however, wind-caused fluctuating movements of the radar antenna or movements of the illuminated targets (such as leaves and branches on trees) can lead to large errors in the polarization phase measurements because the spatial variations are comparable to the wavelength at millimeter-wave frequencies. Hence, it was necessary to develop an alternate measurement technique that is not vulnerable to spatial movements of the radar antenna relative to the illuminated target. The coherent-on-receive polarimetric measurement technique satisfies this condition. To evaluate the accuracy of this technique, the 35-GHz radar system was configured to operate in both the fully coherent mode and the coherent-on-receive mode, and then it was used to measure a target placed on a rotating platform in both modes. Having verified that the two techniques indeed provide the same Mueller matrix, modifications were then made to convert the 94 and 140 GHz systems into coherent-on- receive polarimetric configurations. These results are discussed in [19, 21, and 22]. 2.2 Modeling Polarimetric Backscatter From Terrain 2.2.1 Tree Canopies A radiative transfer model was developed for characterizing the polarimetric radar response of tree canopies at millimeter wavelengths [1]. The model uses a phase matrix whose functional form is based on experimental bistatic measurements obtained under laboratory conditions. The model was found to be in very good agreement with experimental data acquired under field conditions. One of the major conclusions of this study is that the first-order solution of the radiative transfer equation provides reasonable agreement 3

with data for the like-polarized scattering coefficient, but not for the crosspolarized component. However, the second order solution provides very good agreement for both. Related contributions on scattering from foliage are given in references [6, 14, 21]. 2.2.2 Snow-Covered Terrain A radiative transfer model with tthe quasi-crystalline approximation was developed for snow-covered terrain. The model provided very good agreement with experimental observations of dry snow at 35, 94, and 140 GHz. To incorporate the dependence on the wetness profile of snow under wet snow conditions, a hybrid first-order/numerical model was developed and used to predict the temporal variation of the backscatter from snow. The results are available in references [3, 4, 15, 23]. 2.2.3 Polarization Phase Statistics The co-polarized and cross-polarized phase distributions can be obtained from multiple measurements of the scattering matrix of the terrain under observation. Such measurements are provided directly by fully coherent polarimetric systems. However, with coherent-on-receive systems, the quantity measured is the Mueller matrix. Also, theoretical models cannot predict the scattering matrix but can provide results for the Mueller matrix. Until recently, no process was available by which the probability density functions of the co-polarized and cross-polarized phases could be derived from the Mueller matrix. References [7] and [8] contain derivations of new relationships that allow the calculation of the probability density functions in terms of the elements of the Mueller matrix. The new relationships provide a critical link between theory and experimental observations of the polarized phase statistics, making it possible to use theoretical models to relate the phase statistics to physical parameters of the terrain. 4

2.3 Millimeter-Wave Radar Data Base The data handbooks [24, 25] of millimeter-wave radar scattering from terrain contain a compilation of data measured at 35, 94, 140, and 215 GHz, as reported in the open literature. The unique features of these handbooks are: (a) They are the only handbooks available for millimeter- wave radar data of terrain. (b) All data is presented using a uniform format throughout. (c) Only calibrated data that has been augmented with "ground- truth" information is included. 3 LIST OF PUBLICATIONS The following list contains papers published in scientific journals or symposia proceedings, or presented at technical symposia and workshops. All of these papers were generated in support of this project. Journal Papers 1. Ulaby, F.T., T.H. Haddock, and Y. Kuga, "Measurement and Modeling of Millimeter-Wave Scattering From Tree Foliage," Radio Science, Vol. 25, pp. 193- 203, 1990. 2. Haddock, T.F. and F.T. Ulaby, "140 GHz Scatterometer System and Measurements of Terrain," IEEE Trans. on Geoscience and Remote Sensing, Vol. 28, pp. 492-499, July, 1990. 3. Kuga, Y., F.T. Ulaby, T.F. Haddock and R. DeRoo, "Millimeter-Wave Radar Scattering From Snow: Part I-Radiative Transfer Model with 5

Quasi-Crystalline Approximation" Radio Science,.Vol. 26, pp. 329 -341, March 1991. 4. Ulaby, F.T., T. Haddock, R. Austin and Y. Kuga, "Millimeter-Wave Radar Scattering From Snow: Part II-Comparison of Theory with Experimental Observations" Radio Science,.Vol. 26, pp. 343-351, March 1991. 5. Sarabandi, K. and F.T. Ulaby, "High Frequency Scattering from Corrugated Stratified Cylinders", IEEE Trans. Antennas and Propagation, Vol. 39, pp. 512-520, April, 1991. 6. Sarabandi, K., F.T. Ulaby, and T.B.A. Senior, "Millimeter Wave Scattering Model for a Leaf', Radio Science, V. 25, N. 1, pp. 9-18, Jan.-Feb. 1990. 7. Sarabandi, K., "Derivation of Phase Statistics From the Mueller Matrix," accepted for publication in Radio Science, 1992. 8. Ulaby, F.T., K. Sarabandi, and A. Nashashibi "Statistical Properties of the Mueller Matrix of Distributed Targets", Special Issue of IEE Proceedings Part F: Remote Sensing Radar, 1992. 9. Sarabandi, K., Y. Oh, and F.T. Ulaby, "Measurement and Calibration of Differential Mueller Matrix of Distributed Targets", submitted for publication in IEEE Transactions on Antennas and Propagation, February, 1992 Symposia Papers 10. Kuga, Y., Austin, R. T., Haddock, T. F., and F. T. Ulaby, "Calculation of the Diurnal Backscattering Characteristics of Snow at 35 and 95 6

GHz," Progress in Electromagnetics Research Symposium, July 25-26, 1989, Cambridge, Massachusetts. 11. Haddock, T. F. and F. T. Ulaby, "140 GHz Scatterometer Measurements," 1989 International Geoscience and Remote Sensing Symposium (IGARSS '89), July 10-14, 1989, Vancouver, Canada. 12. Sarabandi, K., F. T. Ulaby, T. B. A. Senior, "Millimeter Wave Scattering Model for a Leaf," 1989 International Geoscience and Remote Sensing Symposium (IGARSS '89), July 10-14, 1989, Vancouver, Canada. 13. Ulaby, F. T., T. F. Haddock, and Y. Kuga, "Measurements and Modeling of Millimeter-Wave Scattering from Tree Canopies," 1989 International Geoscience and Remote Sensing Symposium (IGARSS '89), July 10-14, 1989, Vancouver, Canada. 14. Ulaby, F. T., "Millimeter-Wave Bistatic Scattering by Terrain," AGARD Conference Proceedings, Kopenhagen, Denmark, Sept. 9-13, 1989. 15. Kuga, Y., R. T. Austin, T. F. Haddock, and F. T. Ulaby, "Calculation of the Diurnal Backscattering Characteristics of Snow at 35 and 95 GHz," 1989 Progress in Electromagnetic Research Symposium (PIERS), July 25-26, 1989, Massachusetts Institute of Technology, Cambridge, Massachusetts. 16. Ulaby, F. T., Y. Kuga, and T. Haddock, "Measurement and Modeling Millimeter-Wave Scattering from Snow, URSI Specialist Meeting on Signature Problems in Microwave Remote Sensing, 16-18 May, 1990, Hyannis, Massachusetts. 17. Ulaby, F., T. Haddock, and Y. Kuga, "Radar Clutter Database at Millimeter-Wave Frequencies", Workshop on Detection, Discrimination and Classification of Targets in Clutter, 13-15 November, 1990. 7

18. Kuga, Y., F.T. Ulaby, and T.F. Haddock, "Millimeter-Wave Radar Scattering from Snow", Workshop on Detection, Discrimination and Classification of Targets in Clutter, 13-15 November, 1990. 19. Kuga, Y., K. Sarabandi, A. Nashashibi, and F. Ulaby, "Millimeter Wave Polarimetric Scatterometer Systems: Measurement and Calibration Techniques", AGARD, 6-10 May, 1991, Montreal, Canada. 20. Ulaby, F., and Y. Kuga, "Modeling and Measuring Millimeter-Wave Scattering from Snow-Covered Terrain", AGARD, 6-10 May, 1991, Montreal, Canada. 21. Nashashibi, A., Y. Kuga, and F.T. Ulaby, "Polarimetric Observations of Trees at 35 and 94 GHZ", 1991 North American Radio Science Meeting and International IEEE/AP-S Symposium, June 24-28, 1991, The University of Western Ontario, London, Ontario. 22. Kuga, Y., A. Nashashibi and F.T. Ulaby, "Clutter Measurements by Millimeter- wave Radars", National Telesystems Conference, Atlanta, GA, 26-27 March, 1991. 23. Kuga, Y., A. Nashashibi and F.T. Ulaby, "Millimeter-wave Polarimetric Radar Scattering from Snow", IGARSS '91, 3-7 June, 1991, Finland. Technical Reports 24. Haddock, T.F., and F.T. Ulaby, "Millimeter - Wave Radar Scattering from Terrain: Data Handbook I", Radiation Laboratory Technical Report 026247-T-1, The University of Michigan, Ann Arbor, Michigan, January, 1990. 8

25. Haddock, T.F., and F.T. Ulaby, "Millimeter - Wave Radar Scattering from Terrain: Data Handbook II", Radiation Laboratory Technical Report 026247-T-2, The University of Michigan, Ann Arbor, Michigan, September, 1990. 4 PARTICIPATING SCIENTIFIC PERSONNEL The following people participated in the MMW radar program: Faculty and Research Scientists Professor Fawwaz T. Ulaby Professor Yasuo Kuga Dr. Tom Haddock Dr. Kamal Sarabandi Graduate Students Dr. Kamal Sarabandi, received Ph.D (1989) Mr. Adib Nashashibi, expected PhD completion in 1993 Mr. Roger DeRoo, received MS (1990), expected PhD completion in 1993 Mr. Richard Austin, expected PhD completion in 1992 Mr. Steve Ciccarelli, expected MS completion in 1992 Mr. Paul Siqueira, expected PhD completion in 1994 Research Engineer Mr. Ron Hartikka 9

5 CONCLUSIONS The work conducted under this program has led to significant contributions in terms of characterizing the radar backscatter from terrain at millimeter wavelengths, particularly with regard to the new field of radar polarimetry. A major strength of the program has been the strong connection between theoretical modeling and experimental observations. A key area of future research is to explore the relationships between the statistics of the copolarized and cross-polarized phase angles of the backscattered signal and the physical properties of the terrain. The phase statistics provide a new vector of observational parameters which may convey significant information about the observed scene. This area of research is in its infancy and deserves serious exploration. 10

APPENDIX A: REPRINTS OF SELECTED PAPERS This appendix contains reprints of a select set of papers generated under this program. These are: A.1 Ulaby, F.T., T.H. Haddock, and Y. Kuga, "Measurement and Modeling of Millimeter-Wave Scattering From Tree Foliage," Radio Science, Vol. 25, pp. 193- 203, 1990. A.2 Haddock, T.F. and F.T. Ulaby, "140 GHz Scatterometer System and Measurements of Terrain," IEEE Trans. on Geoscience and Remote Sensing, Vol. 28, pp. 492-499, July, 1990. A.3 Kuga, Y., F.T. Ulaby, T.F. Haddock and R. DeRoo, "Millimeter-Wave Radar Scattering From Snow: Part I-Radiative Transfer Model with Quasi-Crystalline Approximation" Radio Science,.Vol. 26, pp. 329 -341, March 1991. A.4 Ulaby, F.T., T. Haddock, R. Austin and Y. Kuga, "Millimeter-Wave Radar Scattering From Snow: Part II-Comparison of Theory with Experimental Observations" Radio Science,.Vol. 26, pp. 343-351, March 1991. A.5 Sarabandi, K., F.T. Ulaby, and T.B.A. Senior, "Millimeter Wave Scattering Model for a Leaf', Radio Science, V. 25, N. 1, pp. 9-18, Jan.-Feb. 1990. A.6 Kuga, Y., K. Sarabandi, A. Nashashibi, and F. Ulaby, "Millimeter Wave Polarimetric Scatterometer Systems: Measurement and Calibration Techniques", AGARD, 6-10 May, 1991, Montreal, Canada. A.7 Nashashibi, A., Y. Kuga, and F.T. Ulaby, "Polarimetric Observations of Trees at 35 and 94 GHZ", 1991 North American Radio Science 11

Meeting and International IEEE/AP-S Symposium, June 24-28, 1991, The University of Western Ontario, London, Ontario. A.8 Kuga, Y., A. Nashashibi and F.T. Ulaby, "Clutter Measurements by Millimeter- wave Radars", National Telesystems Conference, Atlanta, GA, 26-27 March, 1991. A.9 Sarabandi, K., "Derivation of Phase Statistics From The Mueller Matrix," accepted for publication in Radio Science, 1992. 12

Radio Science, Volume 25, Number 3, Pages 193-203, May-June 1990 Measurement and modeling of millimeter-wave scattering from tree foliage F. T. Ulaby, T. H. Haddock, and Y. Kuga Radiation Laboratory, University of Michigan, Ann Arbor (Received August 21, 1989; revised December 18, 1989; accepted December 20, 1989.) Because the constituent elements of a tree canopy, namely the leaves, needles, branches, and trunks, have complex geometries with curvatures and surface roughness scales that are comparable to or larger than the wavelength at millimeter wavelengths, the traditional approach used to compute the phase function of the vegetation volume is totally impractical. In this paper we propose a relatively simple model for characterizing the phase function on the basis of direct experimental data. The model is used in conjunction with a solution of the radiative transfer equation to predict the backscattering behavior of tree canopies. The model is found to provide very good agreement with radar observations made at 35, 94, and 140 GHz. 1. INTRODUCTION The leaves, needles, branches, and trunks comprising a vegetation canopy are lossy dielectric structures with complex geometries. Whereas it may be acceptable to approximate a leaf as a thin, flat disc at centimeter and longer wavelengths, such a treatment is invalid at millimeter wavelengths because the leaf curvature and its thickness are comparable to or larger than A. Similar statements can be made with regard to the size and surface roughness of branches and other components of a vegetation plant or tree. Consequently, it is very difficult, if not impossible, to use numerical quadrature techniques for solving the vector radiative transfer equation [Ulaby et al., 1986; Tsang et al., 1985] to compute the radar backscattering coefficient of vegetation at millimeter wavelengths. The problems encountered are that (1) it is not possible to compute the scattering phase function of the vegetation volume because accurate models for the scattering matrices of the scattering elements (curved leaves, rough-surface branches, etc.) are not available at millimeter wavelengths, and even if such models were available, the numerical computations that would have to be performed to obtain the phase function (which involve integration over size and orientation parameters) would be extremely expensive, and (2) when the phase function has a complicated dependence on the bistatic scatCopyright 1990 by the American Geophysical Union. Paper number 89RS03764. 0048-6604/90/89RS-03764$08.00 tering angles, it is very difficult to compute the solution of the radiative transfer problem beyond the first order. Hence an alternate approach is needed for computing the radar backscatter from vegetation at millimeter wavelengths. In this paper we shall propose a relatively simple model for characterizing the phase function of vegetation canopies at millimeter wavelengths, and then use it in a second-order solution of the radiative transfer equation to compute radar backscattering from trees. Model results are compared with radar backscatter measurements for tree canopies at 35, 94, and 140 GHz. 2. PHASE MATRIX Except for the main trunk, tree foliage consists primarily of randomly distributed leaves (or needles) and branches, most of which are much larger than the wavelength in size (at millimeter wavelengths), have complex shapes, and are characterized by quasi-uniform orientation distributions. These properties suggest that whereas individual scattering elements may exhibit highly complex and polarization-dependent scattering patterns, an elemental volume dV containing many of these elements is likely to exhibit propagation and scattering properties that are weakly polarization-dependent and characterized by a relatively simple scattering pattern. This prediction is supported by experimental observations made by Ulaby et al. [1988] which show that bistatic scattering from trees exhibits comparable results for horizontal and vertical polarizations. 193

194 ULABY ET AL.: SCATTERING FROM TREE FOLIAGE The random nature of the tree foliage supports the use of radiative transfer theory [Ishimaru, 1978] for modeling millimeter wave propagation in the canopy [Ulaby et al., 1988; Schwering et al., 1988]. In the radiative transfer model the formulation is given in terms of the phase function P(Os, ki; Ai, hi) relating the specific intensity scattered by a unit volume of the scattering medium into the direction (Os, Cs) to the specific intensity incident upon the unit volume from the direction (0i, 4i), with both sets of orientation angles being defined with respect to a reference coordinate (x, y, z). The phase matrix represents the average Stokes matrix of the particles constituting the unit volume. To relate P to the properties of the medium, we start by considering scattering by a single particle. For a plane wave with electric field vector E' incident upon the particle in the direction ki = (8i, /i), the far-field wave scattered by the particle in the direction ks = (Os, s5) is a spherical wave with field vector Es. The vertical and horizontal polarization components of Es at a range r from the scatterer are related to the components of the incident field through the scattering matrix S(Os, 4s; Oi, hi) of the particle [Ulaby et al., 1986], a,,h(ks, ki) = [4irr2 ~ [Es12/E I 2] = 4,rISv,hl2 (3) When considering an elementary volume containing N randomly distributed particles per unit volume, we can characterize bistatic scattering by the volume in terms of the bistatic scattering cross section per unit volume (or bistatic scattering coefficient) Kvh = N(COVh) = 4rN(ISvh I2) (4) where angle brackets denotes ensemble average. The vector radiative transfer equation is formulated in terms of the specific intensity vector I defined through the modified Stokes parameters I,,, Ih, U, and V as follows: Ih fl~(lEv 12) IVU '2Re (IEh I') I&I Al = (5) 2I1m (E,ES,) For an elementary volume of length ds illuminated, in the general case, from all directions ki by incident intensity I'(ki), the intensity scattered in the direction ks is given by Es E 1 eikor, =h "-r Eh r LE LE El IS(ks) = f I ds P(ks, k)Ii(ki) dfli (6) (1) The matrix S is given by four scattering amplitudes, St,v Svh Shv Shh where P is the phase matrix given by P11 P12 P13 P14 P(21 P P23 P24 P31 P32 P33 P34 P41 P42 P43 P44 whose elements are related to those of S by (2) For a specified scattered/incident polarization com (IS1512) P~ii, ii) = N i (ISh, 12) P(kk)=N 2(Re (SvSh*v)) 2(Im (S&ShI, )) (ISvh 12) (IShh 12) 2(Re (S hS h)) 2(Im (SnS h )) (Re (SwS*v)) (Re (ShvS h)) (Re (S, S * + SvhS v )) (Im (Svs h + SvhShv)) -(Im (SS *i)) -(Im (ShvS Ih)) -(Im (S S h - SvhS? )) (Re (SIvS h - SvhS v)) (7) bination, the bistatic scattering cross section of the particle is defined in terms of the ratio of the scattered to incident power densities. Fot vh polarization, for example, 3. PROPOSED PHASE MATRIX The 16 elements of P can be readily computed provided we know (1) the number density N, (2) the probability density functions for the sizes, shapes,

ULABY ET AL.: SCATTERING FROM TREE FOLIAGE 195 and orientations of the particles, and (3) the dielectric properties of the particles, and additionally, we have available appropriate models for computing the scattering matrices of the particles. In most cases, this information is not available for terrain surfaces and volumes, which forces investigators to estimate the physical parameters of the canopy and to treat the canopy constituents as spheres, cylinders, and discs. These approximations lead to errors, and because the number of parameters involved is large, it is difficult to assess the sources of error. 3.1. Phase matrix in the scattering plane Instead of using the first approach described above to compute the elements of the phase matrix P, we propose to use a semiempirical approach based on experimental measurements. Ulaby et al. [1988] used a 35-GHz bistatic radar system to examine the scattering patterns of small trees under laboratory conditions. Two types of experiments were conducted: (1) transmission measurements to determine the extinction coefficient Ke for horizontal and vertical polarizations and (2) bistatic scattering measurements in the plane of scattering (de 1.6 m 1.4 m Ficus 5 cm 4~~~~~~~ Arbor Vitae T Fig. 2. Sketches of the tree architectures and photography of a Ficus leaf and an Arbor Vitae branch. fined to be the plane containing the incident and scattered directions and orthogonal to the polarization planes of the waves [Chandrasekhar, 1960]) to evaluate the angular variations of the like- and cross-polarized bistatic scattering cross sections per unit volume K,,,(l/), Kvh( q), Kh,((li), and Khh(ql), where tr is the angle shown in Figure 1. While the transmitter remained in one location with the beam pointing at the crown section of the tree, the receiver was moved in discrete steps around a circle in the horizontal plane with the tree at its center. At each receiver location, defined by the angle i, the average received power was measured and then used to compute KV,(I/), Kvh(/), Khv,(I), and Khh(OlI). The averaging process was realized by placing the tree on a rotating platform and measuring the received power (for a given receiver location) as the tree was rotated over 360~. Two distinctly different types of trees were selected for examination: Ficus and Arbor Vitae (Figure 2). The Ficus tree had small, flat, simple leaves approximately 10 cm2 in area, whereas the Arbor Vitae tree had a branching trunk arrangement with branches supporting needles approximately 1.5-3 I R w = 90 6m Fig. 1. Configuration used for measuring bistatic scattering from tree foliage. The tree was placed on a rotating platform, the transmitter was in a fixed location, and the receiver could be set at any angle dq.

196 ULABY ET AL.: SCATTERING FROM TREE FOLIAGE FICUS TREE ARBOR VITAE TREE E - 3 w 20.4 17.6 14.8 12.0 9.2 6.4 3.6 0.8 -2.0 K =, 4.7 Np.m ia, 0. 1395 w,(Y) p, ~~~~~ 0.1515 Like Polarization C = 0.932 / ~~~~Data Model a2 = 0.0548:2 (W)'"/:":,. Cross Polarization 0.2489 O,r,,,' Data~ra ~~~~~~~~.,,............ --- —- ----.~ '' —.,..'~~~,....~~~~~ I.....'~I..I..;~~l.I.. 17.9 15.0 12.2 9.4 E V D 6.5 3.7 0.9 -1.9 ~ ~~~~~~~,~~~~K 4.7 Np-m'l /, ' r,~ a, 0.0985 Like Polarization 0 6 C 0.9586 Model Data a = 0.0223 K 2(W/)2~ Cross Polarization: 2 =. 1 573 Data / Model....'. I'...,...,.,..-.........'::..;....-....,.... -4.8 -7.6 -4.8 -7.6 -10.4 -80 -64 -48 -32 -16 0 16 32 48 64 80 -80 -64 -48 -32 -16 16 32 48 64 80 11 (Degrees) Fig. 3. Comparison of measured bistatic scattering cross section per unit volume for Ficus tree foliage with calculations based on the model functions given by (14) and (15). mm in length. More detailed information about these test trees and the measurement procedure is given in Ulaby et al. [1988]. The major conclusions derived from the experimental observations that pertain to the present study are: 1. For both types of trees, the like-polarized scattering patterns, K,v(li) and Khh(I), were approximately the same, and a similar result was observed for Kh(tI) and Khv(t/). Thus w (L;'4vees3) Fig. 4. Comparison of measured bistatic scattering cross section per unit volume for Arbor Vitae tree foliage with calculations based on the model functions given by (14) and (15). assume azimuthal symmetry with respect to the forward scattering direction (t$ = 0), the like- and cross-polarized scattering coefficients can be expressed as KI(4/) = KsgIl (1) (11) K2(i/) = Ksg2(WI) (12) and to satisfy (10), the sum ofgI(4i) and g2(2tf) has to satisfy the relation A Khh(4I) ' Kz,v,(I/) = KI(I/) A Khi,(t1) == Kh( 0) = K2(41) (8) (9) -I '[gIl(4i) + g,(ti)] sin ch dl = I 2f0's (13) 2. In spite of the fact that the two trees were markedly different in terms of the shapes and sizes of their scattering elements (leaves, needles, branches), both exhibited similar scattering patterns. Figures 3 and 4 shows plots of the measured values of Kl(i) and K2(4) for the two types of trees. Also shown are plots calculated using the expressions discussed below. The scattering coefficient Ksh for a h-polarized incident intensity is given by K= =~ff [,.h(Iiks) + K1kh() d+i. (10) and a similar expression can be defined for K'. In view of (8) and (9) we shall set K h = Ks' = Ks. If we treat foliage as an isotropic medium and if we In view of the shapes of the measured patterns (Figures 3 and 4), g 1(0) and g9(4J) can each be described as the sum of a relatively weak isotropic component and a Gaussian-shaped, strong and narrow forward-scattering lobe f(4), g!I() = [a If(l() + (I - ai)]C = al exp - +(l- Ca,) C I- IS I A I/ r/ 414) g2(4) = [a f.2(4) + (l - a_)](l- C) = a,(2 2 exp +- 2 + (l - a) (I -C) (15)

ULABY ET AL.: SCATTERING FROM TREE FOLIAGE 197 where pI and I2 are the effective beamwidths of the like- and cross-polarized forward-scattering lobes. These expressions have the following properties: J2 gl(4) sin i di = C (16) 2:92() sin q dS = 1 - C (17) J aft () sin 4 dS * [j (1 - a) sin i d = ca/(l - a) (18) The sum of properties (16) and (17) satisfies (13), the ratio (1 - C)/C represents the ratio of total scattered cross-polarized energy to total scattered likepolarized energy, and the ratio (1 - a)/a represents the ratio contained in the isotropic component to the energy contained in the main lobe. The "calculated" plots shown in the Figures 3 and 4 are based on (11), (12), (14), and (15), with the values of the parameters selected to provide good agreement between the measured and calculated plots. Now let us return to the phase matrix given by (7). The element P l is given by P (,) = N(IS1,,I 2) N 1 = (ap, - Kt() (19) 41' t) ' 47 ( Ks 4 91( 4ir P13 = N(Re (S,,S*S h)) = N(Re (IS,,,le/JllST,h e -j"*Ah)) - N(|S,,,IS,,hIcos (4,. -,vh)) - N(IS,.ll S,.hl)(o s (4,,1, - bvh)) (21) where 04 is the phase of the scattering amplitude S, (and similar definitions apply for the other scattering amplitudes). In the last step of (21) it was assumed that the magnitude IS,,v ISvhl and the phase difference (&, - kAh) are independent random variables. According to 35-GHz radar measurements of the backscattering from rocks [Whitt and Ulaby, 1988] and 1.25-GHz polarimatric data extracted from airborne radar images of forested areas, the phase difference (4, - kvh) is uniformly distributed over [0, 21r]. Hence the average value of cos (4v, - k4h) is zero, and therefore P13 = 0. Similarly, all terms in P involving the product of a like-polarized scattering amplitude and cross-polarized scattering amplitude may be set equal to zero. It was also observed (in the same investigations) cited above that the phase difference (04 - Mhh) corresponding to the product of the like-polarized scattering amplitudes has a Gaussian-like distribution centered at 0~. We shall therefore adopt the approximations (cos (4, - Ohh)) 1 and (sin (4h, - Mhh)) - 0. Furthermore, in view of (8) we shall assume that (lSI IShhl) -- (ISw12). Hence for the terms involving (SvS^h) we have N(Re (S,,Sh)) = N(lSv,,llShhlcos (c - Ohh)) = N(lSz,lIlShh)<COS (os, - 4hh)) = N(IlS,,,12) - P (22) and Similarly, it is easy to show that N(Im (S,,,S th)) = 0 (23) P22(0) = P!!() Upon incorporating the preceding results in (7), we obtain the simplified matrix and K= P PI2(41) = P21 (4) = 4 g2(4t) (20) 9g g2 Ks 92 91 =4 0 0 00 0 0 91 +92 0 0 0 0 g1 -g2J o_ O i 91 - 92 (24) Next, we shall make certain assumptions to simplify the remaining terms. Let us consider the term P13 in P, with gl and g2 as given by (14) and (15), respectively.

198 ULABY ET AL.: SCATTERING FROM TREE FOLIAGE P(Os, &s; Oi, 4i) = L(Tr - Y2)P(41)L(-YI) (25) where L is a linear transformation given by cos2 y sin2 y -sin 2y 0 sin2 e cos2 y sin 2y 0 I sin 2y -I sin 2y cos 2y 0 0 0 0 1 (26) x The angles Yl and Y2 are defined in Figure 5. Fig. 5. Coordinate system showing the incident and scattered Introducing the abbreviations, angles, (6i, Xi) and (8s,:b), and their relation to r,. (l, I) = cos 7l cos 72 - sin yl sin Y2 (27a) 3.2. Phase matrix for any incident and scattered directions The phase matrix P(fi) in (24) is defined in terms of the scattering angle /i as shown in Figure 1. If the incident direction is (0i, ci) and the scattered direction is (Os, Os), denoted by points P1(6i, Xi) and (1, r)= -cos 7y1 sin 72 + sin 71 cos 72 (r, 1) = sin y1 cos Y2 + cos Y1 sin Y2 (r, r) = cos lI cos Y2 + sin yl sin Y2 (27b) (27c) (27d) we can write the phase matrix as P(Os,,s; O, bi) = K4 91(4k)(1, 1)2 + 92(41)(1, r)2 g1(4i)(r, 1)2 + g2(4i)(r, r)2 gl(/i)(r, 1)2 + 92(4i)(r, r)2 gl(4)(l, 1)2 + g2(J)(1, r)2 2[g1(4,)(1, 1)(r, 1) + g2(0)(r, r)(1, r)] - 2[g1(I)(r, 1)(1, 1) + 92(4,)(r, r)(1, r)] 0 0 -gl(O)(r, 1)(1, 1) + 92(O)(r, r)(l, r) gl(1)(r, 1)(1, l) - g2(0)(r, r)(1, r) gl( )[(l, 1)2 _ (r, 1)2] + g2(qi)[(r, r)2 - (1, r)2] 0 0 0 0 91(') - 92(') (28) P2(0s, As) in the polar coordinate system shown in Figure 5, the plane of scattering contains the triangle OP1P2 and the scattering angle I, is given by the angle P10P2. The radiative transfer equation, on the other hand, is written in terms of the polar angles 0 and t. We need to obtain a new phase matrix in terms of 0 and, in order to use it in the radiative transfer equation. The details of the transformation from P(u,) to P(Os, 4s; Oi, Oi), which involves linear transformations through angles yl and Tr - Y2, are described in Chandrasekhar [1960]. The transformed phase matrix is given by Using the cosine and sine laws of a spherical triangle, we can write (1, 1), (r, r), (r, 1), and (1, r) in terms of Os, Asp, Oi and Xi, sin2 (si - i ) (1, l) =. 2 -sin Os sin si sin 41 ~ (cos i - 1) - cos (Oi - bs) (29a) (r sin2 (oi - As) s (r, r) = 2 --- - sin2 i, ~ sin Oi (cos 4, + 1) - cos (0i - 9s) (29b)

ULABY ET AL.: SCATTERING FROM TREE FOLIAGE 199 (r, ) = {(sec Oi + sec Os)[sin (Oi - s) sin O sin Oi cos i cos (ki - (s))] sin2 O - sin. i ) [sin2 Oi sec Oi + sin2 0s sec Os] sin' 2 (29c) (l,r) = {(sec i - sec Os)[sin (Oi - 4s) sin Os sin Oi cos lb cos ('i - OS) ] ( sin ( -sin 2) -sin (0i Os) [sin2 Oi sec i - sin2 Os sec s]} sin 4, I' I Air Z-O Forest Canopy;/' I. Ground Fig. 6. Geometry of the scattering problem. (29d) where cos 4, = cos Os cos Oi + sin 0s sin Oi cos (4i - Os) (30) 4. RADIATIVE TRANSFER MODEL At millimeter wavelengths the penetration depth of foliage rarely exceeds 1 m. Hence it would be reasonable to neglect the backscatter contribution of the underlying ground surface, and in the case of most vegetation canopies it may also be possible to treat the canopy as semi-infinite in depth. In this section we seek an expression for the vector specific intensity Is scattered from a forest canopy characterized by a phase matrix of the form given by (28). To this end, we shall develop a first-order solution and a second-order solution of the radiative transfer equation and then compare the results with the exact solution (based on numerical computations using the quadrature method) and with experimental data. The vegetation canopy is modeled as a continuous, statistically homogeneous, horizontal layer of vertical extent d. The layer has diffuse upper and lower boundaries (Figure 6) at z = 0 and z = -d, respectively, and it is illuminated by an intensity I = Io S(cos 0 - cos o0)8(0 - 4o) (31) incident upon the upper boundary in the direction (r - 00, 0o). Upon solving the radiative transfer equation to obtain an expression for the intensity IS(Os, Os) scattered in any direction (O8, )s), we can set Os = 00 and Os = 7r + 00, which corresponds to scattering in the backward direction, to compute the backscattering coefficient from the equation o0q(00) = 4Ir cos 80oIp(0, 7r + 40). (32) I,(rr'- 6o, 0o) where p, q = v or h polarization. 4.1. Radiative transfer equations When formulating the radiative transfer problem for bounded media, the standard practice is to split the intensity vector into upward-going (I +(Os, Os, z)) and downward-going (I- (ir - Os, Os, z)) components, noting that Os varies between 0 and rr/2 [Ulaby et al., 1986]. In the vegetation layer the intensity I + (Os, s, z) traveling in the upward direction (Os, Os) and the intensity I-(r - Os, (s, z) traveling in the downward direction (Ir - Os, Os) must satisfy the coupled radiative transfer equations d Ke dzI + (AS, O s, z) =- I + (S, Os, z) dz As + F + (s, 4,s, z) (33a) d Ke - I- (-As, Os, z) = — I (-Ls, Os, z) + F (-ps, Os, z) (33b)

200 ULABY ET AL.: SCATTERING FROM TREE FOLIAGE where Ke is the extinction coefficient of the vegetation medium, /,5 = cos Os, and -As = cos (Or - 0s). The source functions F+(s,, Os, z) and F-(-/ts, 4s, z) account for directing the energy incident upon an elemental volume from all directions into the direction (0s, bs) and (Tr - Os, Os), respectively, and are given by F -(,s, s, z) I I t 2 P(As, Os; /i, oi)I+(.ti, 4bi, z) dli f+ |0 P( s u$ P s; -pAi, ki)I (- i, Hi, Z) dAi Jo Jo (34a) F -(-As, Os, z) f21I. +02 | P(-s, s; — is, O)I - ( —i, 4i, z) d fi] (34b) where dfli = diiid4i = sin OidOidOi, and P(xs5, Os; /xi, Oi) is the phase matrix of the vegetation layer relating the intensity incident (upon a unit volume in the medium) in the direction (6i, 4i) to the intensity scattered in the direction (Os, Os). The phase matrix is defined by (28). The solution to differential equations (33a) and (33b) can formally be expressed as I + (Als, >s, z) = e -K(z + d)/I + (Ls, Os, -d) + e -,( - z')/t"F + (A s, &s, z') dz' (35a) I - (-us, Ois, z) = eK/t'lsI - (-LAs, OS, 0) + t e"" (z - z')e/IAF - (-A.s, 4s, z) dz' (35b) Because there is no reflection at the (diffuse) air-vegetation boundaries at z = 0 and z = -d, the following boundary conditions must be satisfied: I &(-/zs, Os, 0)= Io$(/Xs - /Lo)8(4s - (o) I + (5s, s, -d) = 0 (36a) (36b) and because F+ and F- are themselves integral functions of I+ and I-, we have to use numerical techniques involving segmentation in z and (0, 4) in order to obtain an exact solution for Is = I + (6s,,S z = 0). While this may be useful, particularly for comparing with results based on approximate solutions, the numerical technique does not provide much insight with regard to the relative importance of various scattering contributions. Hence we shall use the iterative technique to develop expressions for the first-order and second-order solutions of (35a) and (35b) and then compare their results with the exact results of the numerical solution. The assumptions underlying the iterative technique is that the medium is weakly scattering, i.e., the scattering albedo o = KsIKe << 1. At millimeter wavelength, w = 0.6-0.9 for vegetation [3], and therefore the condition is not satisfied. Nonetheless, we shall now proceed with the iterative technique and then evaluate its usefulness in a later section. 4.2. First-order solution We start with the zeroth-order solutions, which are obtained by setting P = 0 in (34a) and (34b), which renders F o = Fo = 0 in (35a) and (35b), where the zero subscript denotes zero order. Using the boundary conditions given by (34) and (35), the zeroth-order specific intensities are given by Io+ (As', c, z) = 0 (37a) Io (-As, bs, Z) = emrz&MsIos(p.5 - go)8(os - 00) (37b) The zeroth-order solution corresponds to propagation of the coherent wave through the medium with scattering ignored, except for its contribution to extinction. To obtain the first-order solution, we first need to insert (37a) and (37b) into (34a) and (34b) to compute the first-order source function F and F- and then insert the results in (35a) and (35b). This process leads to e -K,.Z//. II (ast, &s, z) = /LsK I ~ [eKIZ - e - 'd]P(s, Cs' k - O, - O)Io (38)

ULABY ET AL.: SCATTERING FROM TREE FOLIAGE 201 II- (-/Ss, s z) = eK" ~~Io0(/Xs -/ o)&(dbs - 4o0) For a thick canopy such that (2Ked sec 00) >> 1, the term in the second square bracket in (42) and (43) reduces to 1. e K, Z/z, + K [- e Kz]P(- Ms, 4S; -Cto, ko)Io wesK2 where (39) 4.3. Second-order solution A K! = Ke(I/ILO + l/AsL) (40) K2 = Ke(1/0L - Il/tLs) (41) The first-order solution for the backscattering coefficient can be obtained by setting z = 0, xs, = AxO, and bs = rr + Oo in (38) and then inserting the result in (32). These steps lead to the expression r '(0 ) 4= cos _,-2K,d sec 0o] 2Ke ' [P(o, ' + o0; -0o, 0)]11 The second-order solution for the backscattered intensity at the surface I+2(0o, ir + ko, 0), can be obtained by (1) replacing (O8, bs, z) in (38) and (39) with (8i, Ni, z'), (2) inserting the resultant expressions in (34a) and (34b) to obtain F + (as, Os, z') and Fl(-As, Os, z'), (3) inserting those expressions in (35a), and finally (4) replacing (/zs, Os, z) with (AxO, Ir + 40, 0). This process leads to 12 (Oo, Vr + 4o, O) eKez'/O FI+ (,uo, + o0, z) dz' d (44) Ks Cos 0 [I - e-2Ked sec O~]g1 (7T) 2Ke (42) where (P) l is the 1 1 element of P and gl (IT) is given by (14) with O = tr. Similarly, the other principal-polarization backscattering coefficients are given by ahh(0o) = a~v,(6o) oa%,(80) = a~,h(o) with F I+ (Io, ir + 0o, z') = -[e KCZ/~op(lao, IT + b0; -/01, O0)IO /Zo r2r r1 e- K'i + J J [e Kez' - e- Kd]P(/xo, IT Jo o LiK3 + Oo; Ai, Oi)P(/i, 4i; -Mo, 1o)Io dfli - Ks cOs o [1- e-2K,d sec 0]g92( r) 2Ke (43) J27r IP e KZ /i iL + AI K4 [1 - e K4z']P(Mo, Ir o MiK4 m 0 t0 o o) o 0 0 c o o Un 0 m 0 SO lNum.10. Co-Pol. Second Order,.- First Order -20.,., f /.,,.-'"~',, Numeric - - - Second Order....... Firsn Order X-Pol. 0e, -30 + Oo; -Mi, 4i)P(-Mi, 'i; -Mo, (o)Io dfli] (45) where A /1 1\ K3 = Ke +) AC i A O A /1 1\ K4 = Ke - \0 Mi (46) (47) I 0.0 0.2 0.4 0.6 0.8 1.0 Albedo o Fig. 7. Backscattering coefficient as a function of albedo at normal incidence, computed with the first-order, second-order, and numerical solutions. Using (14), (15), and (28) to define P, the integrals can be evaluated numerically, and the computed intensity I2 can be inserted in (32) to compute Or~q for any p, q, = v or h polarizations.

202 ULABY ET AL.: SCATTERING FROM TREE FOLIAGE m co a 4, u,.c 4) O. -10 -20 --30 Numencaj --- Second Order First Order.0.. 0.7 8,olo 0 la 0 c a 120 - -10 a o aJ i. c -30 m Co-Pol... *, X-Pol. - - -= 35 GHz ', 0W=0.6 -- Model (Second Order), Co-Pol. ----- Model (Second Order), X-Pol. a Data VV * Data VH -- 0 1 2 3 4 Optical Distance t 5 { 20 40 60 Incident Angle Oo (degrees) (a) 35 GHz 80 Fig. 8. Variation of backscattering coefficient with optical thickness r. 5. MODEL BEHAVIOR AND EXPERIMENTAL OBSERVATIONS Using the phase matrix given by (28) with the parameters measured for the Ficus tree, the copolarized (co-pol) and cross-polarized (x-pol) backscattering coefficients were computed for a variety of canopy conditions in accordance with (1) the first-order solution of section 4.2, (a) the secondorder solution of section 4.3, and (3) the exact numerical si:tion using the quadrature gradient technique [LU,,y et al., 1986]. Figure 7 shows the variation of or~ with the albedo w = Ks/IKe, for a canopy with an optical thickness T = Ked = 1 Np. For an error within 1 dB of the numerical solution the first-order solution is useful up to Co w 0.4 for the co-pol component but is not at all useful for the x-pol component, and the second-order solution is useful up to w - 0.85 for the co-pol component but only useful up to co v 0.5 for the x-pol component. If we relax the error margin to 2 dB for the x-pol component, the useful range of wa may be extended up to 0.85 for the second-order solution. The dependence on optical thickness is illustrated in Figure 8 for all three solutions. For all intents and purposes, a~ is independent of T for - 1I Np. This condition is almost always satisfied for tree canopies at millimeter wavelengths. Comparison of the model behavior with experimental data is provided in Figure 9 which shows measurements of cr~ as a function of incidence angle for a canopy of Spruce trees at 35 GHz and a canopy of Bur Oak trees at 94 and 140 GHz. The model calculations, based on the second-order solution, are in good agreement with the experimental m C) la 00 c 0 Ot C U (a U C 00 0 -co 0 (n1 a0 -10 --20 - Co-Pol. X-Pol. * '. 94 GHz 0)0.8 Model (Second Order), Co-Pol. -... -Model (Second Order), X-Pol. a Data VV * Data VH -30 20 40 60 Incident Angle eo (degrees) 80 (b) 94 GHz Q0 m 00 cl 0 4) z (a U 0 Q ce m -10 -20 Co-Pol. Is 140 GHz c=0.95 --- Model (Second Order), Co-Pol...-..Model (Second Order), X-Pol. ~ Data VV * Data VH 20 40 60 Incident Angle 00 (degrees) 80 (c) 140 GHz Fig. 9. Comparison of theory with experimental observations (a) 35 GHz for a canopy of Spruce trees and at (b) 94 GHz and (c) 140 GHz for a canopy of Bur Oak trees.

ULABY ET AL.: SCATTERING FROM TREE FOLIAGE 203 data at all incidence angles and frequencies except for the 70~ case at 140 GHz (Figure 9c). No obvious explanation is available except for the fact that ao = 0.95, which exceeds the region of validity of the second-order solution. The canopies had continuous crown sections, the trees were about 10 m in height, and the leaves had a moisture content of 53% in the case of the Spruce trees and 27% for the Bur Oak trees. The computations are based on the second-order solution using the Ficus phase-function model shown in Figure 3. The only free parameter used in attempting to match the model results with the data is the albedo w, which was chosen to be equal to 0.6 at 35 GHz, 0.8 at 94 GHz, and 0.95 at 140 GHz. Similar results were obtained in attempting to match the model to experimental observations for (horizontally) continuous canopies comprised of other types of trees. This observation is not surprising in view of the strong similarity noted earlier between the bistatic scattering patterns shown in Figures 3 and 4 for two trees with very dissimilar tree architectures. In other words, the proposed model appears to apply to a wide range of tree types of continuous-crown canopies, with the only major parameter controling the levels of the co-pol and x-pol backscattering responses being the albedo w. In turn, w is strongly dependent on the wavelength and probably dependent on leaf moisture content. Further study is needed to establish the dependence of w on these two parameters. 6. CONCLUSIONS Using the phase matrix model proposed in this study, radiative transfer theory appears to provide excellent agreement with experimental observations of the backscatter from tree canopies at 35, 94, and 140 GHz. The only free parameter used in matching the model to data is the scattering albedo w which appears to depend on only two parameters, the wave frequency and the leaf moisture content. The roles of shape and size of the tree leaves or needles and the tree branch architecture appear to be secondary in importance. Further study is needed to establish the exact dependence of co on moisture content and frequency. Acknowledgment. This work was supported by the U.S. Army Research Office contract DAAG29-85-K-0020. REFERENCES Chandrasekhar, S., Radiative Transfer, pp. 34-35, Dover, New York, 1960. Ishimaru, A., Wave Propagation and Scattering in Random Media, vol. 1, chap. 7, Academic, San Diego, Calif., 1978. Schwering, F. K., E. J. Violette, and R. H. Espeland, Millimeter-wave propagation in vegetation: Experiments and theory, IEEE Trans. Geosci. Remote Sens., 26, 355-367, 1988. Tsang, L., J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing, chap. 3, John Wiley, New York, 1985. Ulaby, F. T., R. K. Moore, and A. K. Fung, Microwave Remote Sensing, vol. III, chap. 13, Artech House, Dedham, Mass., 1986. Ulaby, F. T., T. E. van Deventer, J. R. East, T. F. Haddock, and M. E. Coluzzi, Millimeter-wave bistatic scattering from ground and vegetation targets, IEEE Trans. Geosci. Remote Sens., 26, 229-243, 1988. Whitt, M. W., and F. T. Ulaby, Millimeter-wave polarimetric measurements of artificial and natural targets, IEEE Trans. Geosci. Remote Sens., 26, 563-574, 1988. T. H. Haddock, Y. Kuga, and F. T. Ulaby, Radiation Laboratory, University of Michigan, Ann Arbor, MI 48109.

T-GE/28/4//35831 140-GHz Scatterometer System and Measurements of Terrain Thomas F. Haddock Fawwaz T. Ulaby Reprinted from IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING Vol. 28, No. 4, July 1990

492 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 28. NO. 4. JULY 1990 140-GHz Scatterometer System and Measurements of Terrain THOMAS F. HADDOCK, MEMBER, IEEE, AND FAWWAZ T. ULABY, FELLOW, IEEE Abstract-The goal of the University of Michigan millimeter-wave radar program is to characterize terrain scattering at 35, 94, and 140 GHz. The 140-GHz channel of a truck-mounted scatterometer system has recently been added to give the full desired operating capability. Two injection-locked 45.33-GHz Gunn oscillators use triplers to supply the up- and down-converters. Full polarization capability is obtained through the use of rotatable quarter-wave plates. Real-time signal processing and data reduction takes place in an HP 8510A automatic network analyzer on the truck-mounted platform. Sample measurements of millimeter-wave radar backscattering from vegetation and snow are given. I. INTRODUCTION ILLIMETER-WAVE systems offer the inherent advantages of high resolution, large bandwidth, and small antenna size. In recent years, significant advances have been realized in the development of millimeter-wave components for the construction and operation of imaging airborne radar systems at the atmospheric window frequencies of 35, 94, 140, and 215 GHz. Hence, there is much interest in measuring terrain scattering at these frequencies and in the associated development of theoretical and empirical scattering models. While recent work has been carried out at millimeter wavelengths on trees [1] and snow [2], [3], such data are still sparse, particularly at 140 GHz. The University of Michigan 140-GHz scatterometer system is the latest addition to the network-analyzer-based millimeter-wave scatterometer system, a truck-mounted full-polarization scatterometer developed in support of a program to characterize radar scattering from terrain at 35, 94, and 140 GHz. The basic operation of the scatterometer system has been described in detail elsewhere [4], [5]. Conversion from a swept 2-4-GHz intermediate frequency (IF) to the millimeter-wave RF frequencies is made in the front end, allowing flexible real-time signal processing by the remotely located HP 8510A automatic network analyzer. An outline of the system is given in Fig. 1. An HP 8350B sweep oscillator is swept from 2 to 4 GHz by the HP 8510A network analyzer. After a portion of this IF signal is taken off and fed to the a, referManuscript received October 8, 1989; revised February 2, 1990. This work was supported by ARO Contract No. DAAG 29-85-K-0220. The authors are with the Radiation Laboratory, 3228 Electrical Engineering and Computer Science Building, University of Michigan, Ann Arbor. MI 48109-2122. IEEE Log Number 903583 1. ence port of the HP 8511 frequency converter, it is upconverted within the 140-GHz radar and transmitted to the target. The returned signal is down-converted to the 2-4-GHz range and fed to the b, port of the frequency converter. Signal processing of the return and reference signals takes place within the network analyzer and is sent on the HPIB bus to the HP 9920S computer where final data reduction takes place, and the results are printed out and saved on disk. The HP 8510 computer-control system allows vector error correction of system imperfections through its calibration algorithms. The system has previously operated in this mode at 35 and 94 GHz, and the 140-GHz channel is its latest extension in frequency capability. II. 140-GHz SCATTEROMETER DESIGN A block diagram of the 140-GHz front end is shown in Fig. 2. The transmit portion across the top and the receive portion across the bottom are driven by a common local oscillator (LO) chain. The LO consists of a 45.33-GHz free-running Gunn oscillator, two circulator-coupled 45.33-GHz injection-locked Gunn oscillators acting as amplifiers and two third-harmonic frequency multipliers. This combination provides a nominal output power of 10 dB (1 mW) from each multiplier to power the up- and down-converters. Other combinations of fundamental oscillators, amplifiers or frequency multipliers are possible. This particular combination provided the best combination of performance and cost. Wave-polarization is controlled by a fixed quarter-wave plate followed by a rotatable quarter-wave plate. A 90~ rotation of the movable wave-plate moves the electric field vector through 90~ to give either vertical or horizontal polarization. The polarized signal is transmitted though a conical standard-gain horn with a half-power beamwidth of 1 1.8~. The received RF signal passes though a 3.0-in diameter lens-corrected horn antenna with a half-power beamwidth of 2.2. Hence the antennas' product pattern is essentially controlled by the receiver-antenna pattern, resulting in an effective beamwidth of 2. 1. Receive polarization is determined by movable and fixed quarter-wave plates, in the same manner as the transmit section. The RF signal is down-converted using a tripled 45.33-GHz LO. Since the conversion processes must be phase-coherent, the up- and downconverter LO's are each injection-locked from a central dual-ended Gunn oscillator running at 45.33 GHz. This 0196-2892/90/0700-0492$01.00 ~ 1990 IEEE

HADDOCK AND ULABY: 140-GHz SCATTEROMETER SYSTEM 493 TARG. HP 851i 1A bl FREQUENCY RETURN CONVERTER h2!F TEST SET INTERCONNECT HP. HP 8510A COMI NFTWORK ANALYZER D _ UIDISC DRIVE Fig. 1. Millimeter-wE arrangement gives a phase-coherent LO of sufficient power to supply both up- and down-converters. III. CALIBRATION AND PERFORMANCE For each data set, measurement of a sphere of known size and range is used to generate the 401 VV and HH calibration constants for each of the 401 frequencies in the 2-4-GHz IF band. A calibration target with known cross-polarization response is used for VH and HV calibration. For an incident signal consisting of either pure vertical (or pure horizontal) polarization, a return signal oriented at 450 to vertical is generated by a calibrator consisting of a rectangular standard-gain horn followed by a 38. l-cm-long section of WR6 waveguide with a short on the end. This calibrator is placed in the far field of the 140-GHz radar, and pointed toward the radar with the rectangular aperture of the horn oriented at 45w to horizontal. While a portion of the incident radiation is reflected from the horn, another portion passes from the horn into the waveguide and propagates in the TEV0 mode with the electric field in line with the short axis of the waveguide, which is oriented at 45~ to the horizontal. This signal, comprising equal amplitude vertical and horizontal components, is reflected by the short, and returns to the radar. It can be distinguished from the return from the return from the horn aperture by its longer time delay. The waveguide and short are encased in a metal cylinder to prevent return from the outsid the the guide and flange at the range of the short. Fig. 3 illustrates the W, and HH, VH, and HV responses of the cross-polarization calibration target. While the return from the horn is complex, the return from the short gives a known cross-polarization response. At the range of the horn in Fig. 3. the like- and cross-polarized responses differ by approximately 5 dB. At the range of the short, where for a-perfect radar all four responses would be the same, all returns fall within a + 1 dB range (within experimental uncertainties). Measurement of this signal is used to generate the 401 cross-polarization calibration constants. ave scatterometer system. Transmit Receive i Mixaer r en. Fig. 2. 140-GHz scatterometer front end. Sphere calibration is made on a daily basis, but the standard-gain horn cross-polarization calibration is more cumbersome and is made less frequently. Cross-polarization isolation of the system is typically about 15 dB, and this is checked at each use of the system by making cross-polarization measurements of the sphere. For most natural targets, the cross-polarized return at 140 GHz is only 3-6 dB below the like-polarized return. Hence the cross-polarization isolation of the system is quite adequate at 140 GHz. Noise performance of the system is checked after each calibration by making measurements of the sky at typical target ranges. Table I lists the measured system performance parameters. Fig. 4 illustrates the combined effects of stability and repeatability of calibration of the 140-GHz system over a diurnal cycle. Repeatability of sphere measurements due to pointing only is typically within + 0.5 dB. Variations are the cumulative result of system gain variations and sphere pointing errors. Installation of a controlled heater

494 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 28. NO. 4, JULY 1990 INCIDENT -- WAVEGUDE DELAY UNE I RER RCTED - STANDA DGREEN OERN ORIENTED AT 4S DEGREES TO VERTICAL O> AC m o u 3 o:a o o Cx. Adu m RT,...__ -65 -85 - [ L75 / k II,_ / v -95 41.5 ns 49.5 ns time (ns) Fig. 3. Measured backscatter response of 140-GHz system when observing pyramidal horn antenna connected to short through 38.1-cm-long waveguide section. 0 CD s U, c 0 0 c.li c0 E a) V) 45 40 35, 30 | 140HH 25.. I.. I 10 11 12 13 14 15 1i 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 Time (Hours) Fig. 4. System calibration constant over diurnal cycle, February 27, 28, 1989 (starting 10 A.M.). TABLE I PARAMETERS OF TRUCK-MOUNTED 140-GHz SCATTEROMETER RF frequencies Transmit power RF bandwidth Sweep rate Polarization Product-gain beamwidth Incidence angles Platform height Noise equivalent o~ Stability Measurement repeatability Near-field distance Footprint Signal processing Output products 138-140 GHz -4 dB (1 mW) 0-2.0 GHz I ms/freq., 51, 101, 201, 401 freq/sweep VV, VH, HV, HH 2.10 00-700 2.7 m minimum to 18 m, maximum - -20 dB -0.2 dB/h -I dB 2.7 m minimum: 0.013 m2 maximum: 16.8 m2 HP 8510A/8511A based received power versus range (AR = c/2B) received power versus frequency (at fixed R)

HADDOCK AND ULABY: 140-GHz SCATTEROMETER SYSTEM 495 10 S 5 = 0 a -10 - m -1 I. -15 -20 Amaranthus Retroflexus (Pigweed), about 50 cm tall. Stellana Media (Chickweed) is a low ground-hugging weed. Gravimetric Water Content - 37.6% * * 1. - I... I.. -. __..,.. 140HV *-0 —O 140HV - - * - 140HH 0 10 20 30 Incidence Angle (Degrees) 40 50 60 Fig. 5. Measured backscatter angular response of pigweed (Amaranthus retroflexus) grass over Stellaria media at 140 GHz. on the triple LO unit was required to achieve the excellent system stability observed in Fig. 4 ( + 0.8 dB). For most terrain measurements, a data set consists of measurements of the backscattering coefficient a0 as a function of incidence angle for VV, HH, and VH (or HV) polarizations. The incidence angle is set by an elevation positioner located at the top of the truck-mounted boom. The target is scanned in azimuth to obtain spatially independent samples. For each polarization configuration the number of independent samples, including bandwidth averaging, is at least 50, which corresponds to a measurement precision of about +0.66 dB [6]. Data are tabulated as they are recorded and examined in real time. IV. SAMPLE RESULTS Several types of terrain surfaces and covers were observed by the 140-GHz scatterometer in 1988 and 1989. Sample results are shown next for grasses, trees, and snow. A. Backscatter from Grasses The backscatter plots shown in Figs. 5 and 6 correspond to a field of Amaranthus retroflexus, a spiny weed about 50 cm tall, commonly known as pigweed, over ground cover of Stellaria media, a low ground-hugging weed, commonly known as chick weed. Fig. 5 shows the 140-GHz backscatter response as a function of incidence angle (measured relative to normal incidence) for VV, HH, and HV. Throughout this paper the "receive-transmit"' convention is used. As expected for such a medium, volume scattering effects predominate, and the VV and HH returns are comparable to one another at all incidence angles. The cross-polarized return is approximately 6 dB lower than the like-polarized returns. Fig. 6 shows the HH backscatter response from the same target at all three of the scatterometer operating frequencies: 35, 94, and 140 GHz. The target shows a weak sensitivity to frequency, exhibiting a maximum spread of 5 dB between the three curves. B. Backscatter from Trees Fig. 7 shows 140-GHz measurements of the backscattering coefficient, plotted as a function of incidence angle, for a uniform tree canopy of Thuja occidentalis, commonly known as white cedar. The trees were approximately 10 m in height, and the average water content of the needles was measured to be 56.3%. The like-polarization components (HH and VV) are essentially identical in level and exhibit an approximately cos 0 dependence between 20~ and 70~. The HV component, on the other hand, increases with increasing incidence angle, and its level approaches those of the like-polarization components at 70~. In a separate investigation, the backscattering coefficient at 35 GHz was observed as a function of time over a two-week period for a canopy of deciduous trees (bur oaks). The observation period covered the autumn senescence stage during which the moisture content of the trees decreased. The temporal response of the backscattering coefficient (Fig. 8) exhibited a 3-dB change in level between October 2 and 4 as the leaves underwent a rapid change in moisture content. C. Backscatter from Snow In February and March of 1989, the University of Michigan millimeter-wave system was used to measure the backscatter from snow at a site near Ann Arbor, MI. Figs. 9 and 10 illustrate the angular variation of a0 for wet and dry snow. Fig. 9 shows the response of dry metamorphosed snow with a crystal size of approximately 2.2 mm, and Fig. 10 corresponds to fresh wet unmetamorphosed snow, with crystal size of approximately 1.0 mm and a gravimetric water content of approximately 1.9%. While the angular dependence of the two plots is similar, there is a 4-5-dB level shift between the like-polarized responses of the dry and wet snow targets. At these frequencies, snow is predominantly a volume-scattering medium, and the presence of liquid water in the snow me

496 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 28. NO. 4. JULY 1990 15 10 a o S - 5 0 Cm -5 2 = -10 c I,._ - 35HH -— o — 94HH — _ — 140HH Gravimetric Water Content = 37.6% 881 202.15 -20 -25 0 I 10 20 30 40 Incidence Angle (Degrees) 50 60 Fig. 6. Mcasured angular response of pigweed grass at 35, 94, and 140 GHz for HH polarization. 10 m m 5 -5 ca - z -o R -15 u 10 u, -15 _ _ — - - ---- 140W..... --- 140HV - - I- - 140HH Gravimetric Water Content-56.3% 881116.20 c._ - 0 -oc a u, m 0 10 20 30 40 50 60 7 Incidence Angle (Degrees) Fig. 7. Measured angular response of white cedar trees at 140 GHz. 5 I.. *. ~ * * i - i - i * i - i ' i - i -: 0 -5 10 Full Gre en 15 Leaves Leaves Sparce Dry Drying Leaves 20 25 Z- 35VV(70deg.) 30 ------ 35HV(70deg.) - - - - 35HHt(70deg.) 35 Incidence angle is 70 degrees for all measurements. 880930 40. O -4 ~ 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 October 1988 Fig. 8. Two-week temporal response of 35-GHz backscatter from bur oak (Querces macrocarpa) trees at 70~.

HADDOCK AND ULABY: 140-GHz SCATTEROMETER SYSTEM 497 15 10 I m s 'E.! Ao O g g CD -C 0 0 c9 CD 5 1 0 -5 -10 Snow depth is 10 cm Crystal size approximately 2.2 mm rms surface roughness approximately 0.15 cm - 140VV -— 0 --- 140HV -- --- 140HH 890302., I.... 0 10 20 30 40 50 60 Incidence Angle (Degrees) Fig. 9. Angular response of backscatter from dry metamorphosed snow at 140 GHz. 15 10o E 8 0 0 0 -C Z) 9.X lu 5 Snow depth is 6.5 cm Crystal size approximately 1.0 mm Gravimetric Water Content = 1.9 % IW~ ----* --- —--- L ~ ~ ~ ~ ~ ~~~~~~~~~,._ 0 -54 I - 140W -10 -— o — 140HV I..... —"I — 140HH..... 890220 1I. ) 10 20 30 Incidence Angle (Degrees) 40 50 60 Fig. 10. Angular response of backscatter from fresh wet unmetamorphosed snow at 140 GHz. dium leads to increased attenuation and decreased albedo. For the dry snow, the cross-polarized response is lower than the like-polarized response by about 4 dB, while for the wet snow the difference in level is about 6 dB. Fig. 11 shows 140-GHz backscatter measurements made at an incidence angle of 40~ as a function of time over a 12-h interval extending from noon to midnight on February 27, 1989. The liquid water content measured with a freezing calorimeter for the top 5-cm snow layer and the air temperature are also shown. The backscattering coefficient is observed to exhibit a 3-dB change in level at around 14: 30 in response to the decrease in temperature and liquid water content. At 140 GHz, the penetration depth is on the order of 1 cm, particularly when the snow is wet. As the air temperature drops below 0~C, the snow layer starts to freeze from the top surface downward. Hence, although the liquid water content of the top 5-cm layer may still be greater than zero, the radar responds only to the top 1-2-cm layer and therefore exhibits a time response that leads the temporal variation exhibited by the measured liquid water content. This dependence on penetration depth is illustrated further by the data in Fig. 12 which were measured by the 35-GHz channel. Because of the greater penetration depth, the 35-GHz system exhibits a much more gradual change in level between 1400 and 2400 h. Also, the magnitude of the change in level is 12 dB at 35 GHz, compared to only 3 dB at 140 GHz. At 94 GHz, the measured diurnal pattern (not shown) exhibited a response similar to the 35-GHz data, but with a total change in level of 8 dB. This observed decrease in sensitivity (of the backscattering coefficient to liquid rate content) with increasing frequency is in agreement with earlier observations reported at 35 GHz and lower frequencies [7], [8].

498 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING. VOL. 28, NO. 4, JULY 1990 8 6 c o 2.4 a -6 ____-x................. Indence Angie -40 deg. ~{~*'0'...... -— 140VW.... — 140HV ----- 14HH 8902278......~~~ -1 -10 i r I I I - i Ia I. 12 13 14 15 16 17 18 19 20 21 22 23 24 Time (hours) (a) 6 - _ o A -- I ii -l -------— o "O.4.8 e -12.4 0 12 13 14 15 16 17 18 19 Time (hours) (b) 20 21 22 23 - LWC (%) - -4. -- Temp. (deg. C) 2 Fig. 11. Diurnal measured temporal variation: (a) 140-GHz radar backscatter, and (b) air temperature and snow gravimetric liquid water content (top 5-cm layer). 0 - 5 0.10 - 15 -.20. v..to a' bO~ " 0 a SP" "a d a a.a a..a. a' * - - - -I —a Incience Angie. 40 deg -- 35W.... *.35VH — _ —3SHH 890227/8;o 12 15 18 21 24 27 30 33 Time (houn) (a) 8.,.. ~..,. 6 - ~ 4 i v._ __-4 2 -8 O... - - - — ' - * -12 a o I. 12 15 18 21 24 27 30 33 Time (hours) LWC (%) Time (houm) -- -,4 — Temp. (deg. C) (b) Fig. 12. Diurnal measured temporal variation: (a) 35-GHz radar backscatter, and (b) air temperature and snow gravimetric liquid water content (top 5-cm layer).

HADDOCK AND ULABY: 140-GHz SCATTEROMETER SYSTEM 499 V. CONCLUSION This paper describes the operation of a 140-GHz scatterometer system with a measured accuracy of 1 dB. Sample measurements of terrain backscatter at 35, 94, and 140 GHz are shown for grasses, trees, and snow. In all cases, the angular dependence is approximately as cos 0. The two like-polarized components (HH and VV) exhibit essentially identical levels, and the cross-polarized response is anywhere from 1 to 6 dB below the like-polarized responses, depending on target type. At 140 GHz, the backscatter from snow exhibits a dynamic range of about 3 dB, in response to change in liquid water content. REFERENCES [11 R. M. Narayanan, C. C. Borel, and R. E. McIntosh, "Radar backscatter characteristics of trees at 215 GHz," IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 217-228, May 1988. [2] E. P. Baars and H. Essen, "Millimeter-wave backscatter measurements on snow-covered terrain," IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 282-299, May 1988. [3] L. D. Williams, J. G. Gallagher, D. E. Sugden, and R. V. Birnie, "Surface snow properties effects on millimeter-wave backscatter," IEEE Trans. Geosci. Remote Sensing, vol. 26, no. 3, pp. 300-306, May 1988. [4] F. T. Ulaby, T. F. Haddock, J. R. East, and M. W. Whitt, "A millimeterwave network analyzer based scatterometer," IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 75-81, Jan. 1988. [5] M. W. Whitt and F. T. Ulaby, "Millimeter-wave polarimetric measurements of artificial and natural targets," IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 562-573, Sept. 1988. [6] F. T. Ulaby, T. F. Haddock, and R. T. Austin, "Fluctuation statistics of millimeter-wave scattering from distributed targets," IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 268-281, May 1988. [7] W. H. Stiles and F. T. Ulaby, "The active and passive microwave response to snow parameters: Part I —Wetness," J. Geophys Res., vol. 85, no. C2, pp. 1037-1044, Feb. 1980. [8] F. T. Ulaby and W. H. Stiles, "The active and passive microwave response to snow parameters: Part 2-Water equivalent of dry snow," J. Geophys Res., vol. 85, no. C2, pp. 1045-1049, Feb. 1980. Thomas F. Haddock (M'86) was born in Washington, DC, on November 2, 1949. He received the B.A. degree in mathematics and the M.S. and Ph.D. degrees in physics from the University of Michigan, Ann Arbor, in 1972, 1977, and 1984, respectively. From 1984 to 1985 he was Manager of Development Projects at Applied Intelligent Systems, a machine vision firm involved in real-time optical, infrared, and X-ray vision systems. He is currently with the Radiation Laboratory and the Department of Electrical Engineering and Computer Science, University of Michigan. He has conducted research in the fast flux density variations of quasi-stellar objects at a wavelength of 12.5 mm. Other research has included development of real-time alphanumeric character recognition algorithms and ultrasonic weld inspection algorithms. Prior to receiving the Ph.D. degree, he worked as Applications Engineer for Sarns/3M. a manufacturer of heart-lung machines and cardiac assist devices, where he developed electrodes for manufacturing applications. His current research interests are millimeter-wave scattering and emission from natural targets. Dr. Haddock is a member of the American Astronomical Society. Fawwaz T. Ulaby (M'68-SM'74-F'80) received the B.S. degree in physics from the American University of Beirut, Lebanon, and the M.S.E.E. and Ph.D. degrees in electrical engineering from the University of Texas, Austin, in 1964, 1966, and 1968, respectively. He is currently Professor of Electrical Engineering and Computer Science of the University of Michigan, Ann Arbor, and Director at the NASA Center for Space Terahertz Technology. His current interests include microwave and millimeter wave remote sensing, radar systems, and ratio wave propagation. He has authored several books and published many papers and reports on these subjects. Dr. Ulaby is the recipient of numerous awards including the IEEE Geoscience and Remote Sensing Distinguished Achievement award in 1983, the IEEE Centennial Medal in 1984, and the Kuwait Prize in applied science in 1986.

Radio Science, Volume 26, Number 2, Pages 329-341, March-April 1991 Millimeter-wave radar scattering from snow 1. Radiative transfer model Yasuo Kuga, Fawwaz T. Ulaby, Thomas F. Haddock, and Roger D. DeRoo Radiation Laboratory, Department of Electrical Engineering and Computer Science University of Michigan, Ann Arbor (Received May 17, 1990; revised October 30, 1990; accepted November 6, 1990.) Millimeter-wave (MMW) remote sensing of ground snow has attracted considerable interest in recent years. Because the size of the snow ice particle is comparable to the wavelength in the millimeter-wave region, we can no longer use a simple Rayleigh phase function or the small particle approximation usually used at microwave frequencies for calculating the extinction coefficient. In this paper we present a model for MMW scattering from snow using the vector radiative transfer theory and a Mie phase function. Assuming snow to consist of randomly distributed spherical particles embedded in a mixture of air and water, the vector radiative transfer theory is solved using the discrete ordinate method. The values of the extinction coefficient used in the calculations are based on a combination of experimental data and calculations using the quasi-crystalline approximation. The backscattering coefficient is calculated for different liquid water contents at 35, 95, and 140 GHz. We show that the backscattering coefficient is sensitive to liquid water content at all three frequencies, with 35 GHz being the most sensitive. Except for normal incidence, the effect of snow surface roughness is negligibly small for dry snow, and it is somewhat significant for wet snow at 35 GHz, but not at the higher frequencies. 1. INTRODUCTION Although several papers have appeared recently documenting the results of millimeter-wave (MMW) radar observations of snow-covered terrain [Stiles and Ulaby, 1980; Ulaby and Stiles, 1980; Hayes et al., 1984; Hallikainen, 1985; Baars and Essen, 1988; Williams et al., 1988; Ulaby et al., 1988; Currie et al., 1988; Narayanan and Mcintosh, 1990], the interaction mechanisms responsible for the observed radar response are not well understood at the present time. At millimeter wavelengths, snow is a highly lossy medium, particularly when wet; consequently, the penetration depth is only of the order of a few centimeters [Ulaby et al., 1986, p. 1608; Hallikainen, 1985; Hallikainen et al., 1986]. For dry snow the attenuation is dominated by scattering because the ice particles are comparable to the wavelength in size, and for wet snow both absorption and scattering are important. The physical parameters that exhibit the strongest importance on the radar backscatter from snow are snow surface roughness, crystal size, snow depth, and the liquid-water profile with depth. This Copyright 1991 by the American Geophysical Union. Paper number 90RS02560. 0048-6604/91/90RS-02560$08.00 paper, which is part 1 of a two-paper sequence, provides a radiative transfer model for characterizing MMW scattering from snow, using experimental data in conjunction with the quasi-crystalline approximation [Tsang et al., 1985, pp. 461-469] to compute the extinction coefficient of the snow medium. Part 2 [Ulaby et al., this issue] describes the results of experiments conducted at 35, 95, and 140 GHz, and includes comparison between theory and experiment for certain cases. 2. SNOW MODEL In our snow model we assume ground snow to consist of spherical ice particles embedded in a background medium. Liquid water, when present, is included as part of the background medium. The size of the water inclusion is usually much smaller than the wavelength for millimeter-wave remote sensing. Therefore it is reasonable to assume that the water is uniformly distributed in the snow and the dielectric constant of water can be included as a part of the background. Thus the wet snow medium is modeled in terms of ice crystals in a "wet air" background. It is also possible to consider the snow as lossy particles with a thin film of water surrounded by air, in which case the dielectric constant 329

330 KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I TABLE 1. Dielectric and Propagation Properties of a Snow Medium With Ice Volume Fraction f = 0.4, Average Ice Particle Diameter of 1 mm, and Number Density N = 7.64 x 108 eb K(QCA) m,,, % Real Imaginary Kag, m-1 Real Imaginary Re[n(QCA)], m-I Ke(QCA), m35 GHz Eice = 3.15 + i0.003 ar, = 1.2 x 10-8 Ca = 5 x 10-10 Ke(EFA) = 8.24 m-l 0 1 0 0 930 0.269 1.268 0.538 1 1.0427 0.0129 5.56 942 0.267 1.285 0.534 2 1.104 0.0532 22.2 958 0.264 1.307 0.528 3 1.1776 0.0869 35.2 977 0.262 1.332 0.524 4 1.2405 0.1268 50.0 992 0.259 1.353 0.518 5 1.3046 0.1723 66.2 1007 0.257 1.373 0.514 95 GHz eic, = 3.15 + iO.0085 a, = 5.89 x 10-7 aa = 4.28 x 10-5 Ke(EFA) = 402 m-1 0 1 0 0 2562 11.35 1.2884 22.69 1 1.0296 0.0138 16.2 2582 10.91 1.298 21.82 2 1.0626 0.0397 45.9 2604 10.43 1.3095 20.86 3 1.1055 0.0697 78.7 2632 9.84 1.323 19.68 4 1.141 0.0961 107.2 2654 9.39 1.335 18.78 5 1.1765 0.1244 136.6 2676 8.95 1.346 17.9 140 GHz Rice = 3.15 + iO.012 a, = 1.95 x 10-6 aa = 2.33 x 10-8 Ke(EFA) = 1331 m-l 0 1 0 0 3876 56.58 1.32 113.15 1 1.0254 0.0115 20 3891 54.95 1.326 109.9 2 1.0534 0.0307 52.67 3909 53.94 1.331 107.9 3 1.0883 0.0540 91.15 3933 53.36 1.339 106.7 4 1.1181 0.0741 123.38 3955 53.2 1.347 106.4 5 1.1482 0.0953 156.54 3977 53.1 1.354 106.3 Here mv: liquid water content by volume, Eb: background dielectric constant, Kag: background absorption coefficient, K(QCA): propagation constant in the snow layer obtained by QCA, Re[n(QCA)]: real part of the effective index of refraction obtained by QCA, Ke(QCA): extinction coefficient of ice particles obtained by QCA, K.(EFA): extinction coefficient of ice particles obtained by EFA, o, and ca: total and absorption cross sections of ice particles obtained by EFA. of the water inclusions is part of the dielectric constant of the lossy particles. This approach is inappropriate, however, because liquid water in snow usually occupies the spaces between adjacent ice crystals rather than coat the crystals [Ulaby et al., 1986, p. 2067]. From ground-truth data we know that the ice particles have an average diameter of the order of 0.1-2 mm and their shapes are round but nonspherical. In our model we treat the snow ice particles as spheres with a normal-size distribution with an average diameter of 1 mm and a standard deviation of 0.2 mm. The values of the dielectric constant of the ice particles needed in this study were obtained from Mtiizler and Wegmiiller [1987] and are listed in Table 1. Also from the ground-truth data measured in conjunction with the experimental observations reported in part 2, the volume concentration of snow is approximately 40% and the snow depth is 0.45 m. These values are used in our model calculations. 2.1. Background absorption by water inclusions The imaginary part of the background dielectric constant is directly related to the background absorption coefficient. We use a dielectric mixing formula to calculate the background dielectric constant, and we assume that scattering by the water particles to be much smaller than absorption. The absorption coefficient Kag of the background medium is given by K ag = 2 Im {kb(l -f)} Kag = 2ko(l -f) Im {I } (1) (2) where ko is the free space wavenumber, Eb is the background dielectric constant, and f is the volume fraction of ice particles. To obtain the dielectric constant of the background, which is a mixture of air and water, we use the Van Santen mixing formula [Ulaby et al., 1986, p. 1609]:

KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW. I 331 Lb 1 + — -1~ 1 F 3b I +A. —1) 3,=cI +A — I l where eq, is the dielectric constant of water, m., is the volumetric snow wetness, and Au is the depolarization factor which depends on the shape of the water droplet. If the water droplet is spherical, Au is constant and given by Au = 1/3, but the water particles in snow are usually nonspherical and they change shape with the amount of liquid water present [Hallikainen et al., 1986]. If the liquid water content is low, known as the pendular regime, the values of Au are close to those for a needle. On the other hand, if the liquid water level is high, known as the funicular regime, the values of Au are close to those for a disk. Hallikainen et al. [1986] provide plots for Aa, Ab, and Ac as a function of my, which show that the transition from the pendular to funicular regime occurs at around ml, = 2.5. These plots are used in all calculations of Eb in this and the following papers (see Table 1). 2.2. Scattering characteristics of ice particles The scattering characteristics of ice particles are calculated using the Mie solution [Ishimaru, 1978, p. 27]. The background dielectric constant is assumed to be the real part of eb listed in Table 1. The imaginary part of Eb is not used because of the difficulty in calculating the Mie solution when the imaginary part of the normalized dielectric constant is negative. The average total and absorption cross sections are shown in Table 1. 2.3. Extinction coefficient of dry snow When snow is dry, the attenuation at millimeter wavelengths is due mainly to scattering by the snow particles. In a sparsely distributed medium in which the correlation between particles can be neglected, the effective field approximation (EFA) can be applied and the extinction coefficient is linearly proportional to the concentration of particles [Tsang et al., 1985, p. 458]. Ground snow, however, has an ice volume fraction of 10 to 40%, and the dielectric constant of ice is much larger than that of the background medium. Hence the independentscatterers assumption is inappropriate for snow, which means that correlation between adjacent ice particles should be considered [Tsang et al., 1985, chaps. 5 and 6]. Two additional important phenomena must be considered when treating propagation in a dense medium like snow, namely the backscattering-enhancement effect and the decrease of the extinction coefficient when the density is high [Ishimaru and Kuga, 1982; Kuga and Ishimaru, 1984; Tsang and Ishimaru, 1984]. Backscattering enhancement is caused by the constructive interference of two waves propagating in opposite directions and is important only in the backscattering direction when the phase difference of the two waves is zero. The backscattering enhancement effect has been observed for both randomly distributed discrete particles and very rough' surfaces. The importance of the backscattering enhancement effect that has been recognized in optics and solid state physics, but its significance in microwave and millimeter-wave remote sensing has not yet been evaluated. Since the angular width of the backscattering enhancement pattern is much less than 1~ for discrete particles, the backscattering enhancement effect is not observable if the detector's field of view is large. In most microwave and millimeter-wave systems, the receiving cone of the antenna is much larger than 1~ and the observation configuration is not truly monostatic. Hence the backscattering-enhancement effect is probably not significant and may be ignored, at least to first order, although we suspect that it may be responsible for the difference in level between the calculated and observed backscatter at 95 GHz and 140 (antenna beamwidth about 2~) reported in part 2 of this sequence [Ulaby et al., this issue]. When the volume fraction of particles is more than 1%, the extinction coefficient is no longer linearly proportional to the number density. The deviation from the linear relationship applicable at low densities is related to the size parameter, the dielectric constant and volume fraction of the particles. Extensive experimental and theoretical studies on the extinction coefficient in a dense medium have been conducted in recent years [Ishimaru and Kuga, 1982; Tsang et al., 1985]. The theoretical models include Twersky's model, the perturbation solution with hole-correction, the quasi-crystalline approximation (QCA) with Percus-Yevick pair correlation function, and the quasi-crystalline approximation with coherent potential (QCA-CP). Twersky's model is simple but it is applicable only for

332 KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 1 1000 9 GHz QCA C 24 5 GHz EFA 35 GHz QCA.1 0 1 2 3 4 5 6 Water Content mv(%) Fig. 1. Optical distance T versus liquid water content m,, at 35 and 95 GHz. Snow depth is 0.45 m, the ice volume fraction is 0.4, and average ice particle diameter is I mm. EFA is the effective field approximation and QCA is the quasi-crystalline approximation. small particles. The formula based on the holecorrection is valid if the volume fraction is much less than 10%. For higher concentrations, QCA and QCA-CP with the Percus-Yevick pair correlation function have been shown to be effective [Tsang et al., 1985 p. 479]. In thte millimeter-wave region where the size parameter is close to 1 we cannot use a smallparticle approximation. We need to solve the QCA numerically. We calculated the extinction coefficient of snow using the QCA with the PercusYevick pair distribution function at 35, 95, and 140 GHz. The results are listed in Table 1. When the volume fraction f is 0.4, the extinction coefficients given by QCA are only 6.5% of those calculated according to the EFA at 35 GHz and similar percentages apply at 95 and 140 GHz. The optical distance, which is defined as T = (Ke + Kag)d, is shown in Figure 1 as a function of liquid water content for both EFA and QCA. According to recent extinction measurements conducted for dry snow at 35 and 95 GHz [Hallikainen, 1985], the extinction coefficient was found to exhibit a strong dependence on snow type. The reported values of the extinction coefficient covered the range between 0.96 and 15.4 m - at 35 GHz and between 1.9 and 30.7 m-l at 95 GHz, with esti mated median values of about 3.7 m-l at 35 GHz and 19 m - at 95 GHz. The median value at 95 GHz is close to that computed using the QCA method, but the median extinction coefficient at 35 GHz is much higher than that computed according to the QCA method. If the optical distance is greater than 5, the backscattering coefficient of dry snow becomes essentially independent of the optical distance (in the millimeter-wave region). For a snow thickness of 0.45 m and extinction coefficient of 19 m-I the optical distance is already more than 8 at 95 GHz. Therefore the exact value of extinction coefficient is not important at 95 and 140 GHz. However, an accurate estimate of the extinction coefficient is important at 35 GHz because - is smaller. In our model calculations with a volume fraction f = 0.4 we used extinction coefficients of 3.7 m-1 at 35 GHz and 19 m - I at 95 GHz, which seem to offer the best fit to the experimental data. For snow with other values off the extinction coefficient is scaled in proportion to the change in the QCA model. This method is chosen because the extinction coefficient is a nonlinear function off and the density dependence curve given by QCA has been confirmed experimentally [Ishimaru and Kuga, 1982]. With this correction, for example, the extinction coefficient at 35 GHz becomes 5.18, 4.45, and 3.7 m-I when the snow density f is 20, 30, and 40%, respectively. Because of the lack of measured extinction-coefficient data at 140 GHz, we will use the value calculated by the QCA model. 3. RADIATIVE TRANSFER THEORY Microwave remote sensing of ground snow has been studied by many researchers in the past [UIaby et al., 1986, chap. 13; Tsang et al., 1985, chap. 4]. The model is usually based on the radiative transfer theory, and the Rayleigh phase function is used for modeling the ice particles. This is a good approximation for microwave remote sensing because the ice-particle size is much smaller than the wavelength. For millimeter-wave remote sensing, where the ice particle size is comparable to the wavelength, the Mie phase function must be used instead. In a dense medium like snow we need to modify the conventional radiative transfer theory to take into account the correlation between particles. A dense-medium radiative transfer theory was re

KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I 333 cently developed using the Dyson equation with QCA-CP and the Bethe-Salpeter equation under the ladder approximation of correlated scatterers [Tsang and Ishimaru, 1987; Wen et al., 1990]. The form of the dense-medium radiative transfer theory is the same as the conventional radiative transfer theory and, therefore the same numerical techniques can be used for its solution [Cheung and Ishimaru, 1982; Ishimaru et al., 1982; Wen et al., 1990]. The difference between the conventional radiative transfer theory and the dense-medium radiative transfer theory is the extinction rate, which can be obtained by the QCA-CP and the new form of the albedo. The phase function is still the single particle phase function, which can be approximated by the Rayleigh phase function when the particles are small in size compared to the wavelength. In this paper we shall use a modified version of the conventional radiative transfer model, but the formulation is not as rigid as the one given by Wen et al. [1990]. We will use the extinction coefficient values given at the end of the previous section rather than those based on the effective field approximation. This seems to give the best fit to the experimental data. The effective dielectric constant of the snow layer will be calculated using the QCA and the real part will be used for calculating the reflectivity and transmissivity of the air-snow boundary. The phase function will be calculated using the Mie theory. The background absorption will be calculated using a mixing formula and included in the total extinction coefficient. This will reduce the effective value of the albedo when the liquid water content increases in the background. 3.1. Problem formulation Air el- 1 'OO OQ d O snow 0 rI e2 x= O 0 ground e3- 3.17 Fig. 2. Geometry showing the snow model. previously in terms of the Fourier-series expansion for the Stokes vectors [Ishimaru et al., 1982]. If the incident wave is normally incident and linearly polarized, only two terms in the Fourier series are necessary, but the oblique-incidence case requires all the components of the Fourier series. In this paper we will briefly describe the formulation. The details can be found elsewhere [Cheung and Ishimaru, 1982; Ishimaru et al., 1982; Ma et al., 1990]. For convenience we will use the modified Stokes parameters (II, 12, U, V), II = (E1ET) (4a) 12 = (E2EN) (4b) U = 2Re(EIEl) (4c) V = 2Im(EIEl) (4d) where El and E2 are the electric field in the 6 and 4 directions. The equation of transfer for the incoherent specific intensity I propagating in the s direction (Figure 2) is given by We consider a plane-parallel medium containing spherical particles as shown in Figure 2. A linearly polarized wave, which can be either vertically or horizontally polarized, is obliquely incident with incident angles 80 and 0o. The dielectric constants of media 1, 2, and 3 may be different. In our model medium 1 is air, medium 2 is snow, and medium 3 is the ground. We assume that E1 and E3 do not change with temperature, but E2 (the background dielectric constant of medium 2) varies in response to changes in liquid water content of the snow layer. The formulation of the radiative transfer theory for an oblique-incidence case has been derived ds t(F, ))])] + Ke d [S] [s(-Kel, S)] -,[( s] + Ke da'[S] ds ~~~~~~~~~~~~~~~~~~~~~~r [I(T, S')] + Ke[Ii] where II [I]= = 4 x lI incoherent specific intensity matrix, (5)

334 KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I [S] = [S0] = 4 x 4 Mueller matrix, [Ii] = 4 x 1 incident specific intensity matrix (defined later by (14) and (22)). The Mueller matrix (which also is called the scattering matrix) is expressed in terms of the scattering amplitudes f 1, f12, and f22 [Cheung and Ishimaru, 1982]: In a sparse medium with number density N and Ke = N(a,) the above equation becomes Ke (a,) Ke (ort) + (Kag/N) (10) This form is similar to the albedo of a single particle which is defined as O's U-S wo =- = ra, a + U a (11) / (Ifll>2 IS] = Ia (Ifi)2l [S] = ~(- ) 2 Re (f, I fl,) 2 Im (fil ftl) (If 21)2 (f221)>2 2 Re (fl2ft2) 2 Im (fi2J'2) where a, and oa are the particle's scattering and absorption cross sections. If we define the effective albedo as Re (fil fT2) Re (f21 f12) Re (f, I f2 + fi2ftl) Im (fi X f2 + f2 fl) -Im (fl ff2) -Im (f21 fA2) -Im (fil f2 - fi2ftl) Re (fi, f2 - fi2ft1l) ars aot + (K ag/N) (12) (6) where ( ) denotes ensemble averaging over particle size distribution and (oAt) is the extinction cross section of the particles. In ia the optical distance r and the extinction Ke, equation of transfer becomes d[I(r, S)] Ke d = -[I(i', S)] + -- I d'[S][I(i, s')] do~ ~ ~ Ke J4r Ke + - [PI] Ke with we observe that it decreases as the background absorption increases. For example, the effective albedo of dry snow at 95 GHz is Wo = 0.99. If the -age background is wet and the snow wetness is 5%, the s of effective albedo becomes Wo = 0.69. A similar the 0 calculation at 35 GHz shows that Wo decreases from 0.95 for dry snow to Wo = 0.09 for a snow wetness of 5%. In our model we assume Ke/Ke to be independent of the density and given by (9). For calculating the optical distance and the reflection coefficients, however, we used QCA to obtain the effective propagation constant of the snow. We chose this approach because if QCA is used for calculating Ke/K, the backscattering coefficient becomes very small for slightly wet snow and it does not agree with the experimental data. Our approach (8a) seems to give a reasonable fit to experimental data, as discussed in Part 2. One problem with our (8b) approach is that when r becomes small, the intensity of dry snow may become less than that for wet (8c) snow. For millimeter-wave remote sensing, how(8d) ever, r is large and this deficiency is not critical. (8d) / = cos 0 dw' = dA' d+' dr = LKe' ds Ke = Ke + Kag Absorption by the background medium (wet air) is 3.2. Incident specific intensity now included in Ke in (7). If the background absorp- for vertical polarization tion is zero, (7) reduces to (2) in the work by Because of reflection at the boundaries, we have Ishimaru et al. [1982]. However, if Kag is not zero, upward and downward traveling incident waves, Kag acts like an additional loss term in the particle's extinction cross section, [Iil = [I/+] + [I1-] (13) Ke Ke where [1/+] is for i > 0 and [1/] is for /t < 0. For a K; Ke + Kag vertically polarized incident wave 1I[1 0 0 0]T inci

KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 1 335 dent along the direction (A0o, 0), [1:+] is given in terms of the reduced incident intensity Ir+ by i+ ] = f [S][,r+ ] dw' = oC,, f [S ]( I)'- -8o )6(')e -e d,4' d4' 0 ISi0 (AW)e [li- ] = f [S[lr] dw' = Io CpR23 f [S]( 5(A' +,o)48(4') dL' do' exp ( 2ro - r) /Fl F21 = Io CpF F31 F41 / = 7 e 4Mo (14) Fit =I0C R23 F21 ( \ F31 e 4 — Lo 2r - r (22) with (IfI,, 12) F11 = (CT,) 2 Re (fl f5l) F31 = -xo = Cos 00, 0o = 0; evaluated at / = go; (lf2,12) F2, = (,I) (15) 3.3. Incident specific intensity for horizontal polarization For a horizontally polarized incident wave given by Io[O 1 0 0]T, [I+] is given by 2 Im (fil f.l) F41 = Fll, F21, F31, F41 are F12 i+ = o F42 exp,F42 A-' AO K AOo (23) Cp I - R23R21 exp 2To (16) (17) (18) with n2 COS 02 it1 nt2 = c 6 12 nR l cos 01 R2, = lr2'v12 (lf1212) F12 -- (a,) 2 Re (fi2fl2) F32 = (ort) (lf2212> F22 = (24) 2 Im (f12ft2) F42 = (or,) R23= Ir'312 (19) where the superscript v denotes vertical polarization, nj is the index of refraction of mediumj, and ryi and ty are the Fresnel reflection and transmission coefficients at boundary ij, The reflected specific intensity for horizontal polarization is given by FI2 F22 ( ~ -) [li-] = IoCR23 F32 exp - F42, = — O (25) nj cos Oi - ni cos Oj r = r nj cos 6i + ni cos Oj 2nj cos Oi i nj cos Oi + ni cos Oj (20) The quantities Cp and R23 in (23) and (24) are given by the same expressions given previously for the vertical polarization case, except for replacing the (21) superscript v with h everywhere and using The reflected incident specific intensity [I- ] is given by ni cos 0i - nj cos Oj r cos ' ni cos Oi + nj cos Oj (26)

336 KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I 2ni cos Oi where h t. U (27) ni cos Oi + nj cos Oj [I]o = [lo2 [L]o = [S]o 3.4. Outline of the solution The upper and lower boundary conditions for the incoherent specific intensity [I] inside medium 2 are given by (33) F +] = F210, S =,o -F21 =A IF[]= [F-]o = 1 = [l(T = 0, -A) = [R21O(-)][I(T = o, A)] (28) -1 < / < 0 For m > 0 d[I]m KKe K -[i]m + [L]m[l]m d/' + Ke JI e [I(r = r0, -A.)] = [R23(A)][I(r = To, a)] (29) 0< 1<l where [R21 (x)]) and [R23(,/)] are the reflectivity matrices given by the generic form Ke * Cp[F + ]me -(/ o.) + - CpR23[-F]me-(2 - t)/ -Ke (34) Ir/jl2 0 li ~1 12 0 pr, l2 [Rij(A)l = o 1O r O O 0 0 0 0 Re (r~vr.*) Im (rgr.*) where It is possible to eliminate the 4 dependence fror and obtain a new equation in terms of r an [Cheung and Ishimaru, 1982; Ishimaru et al., 19 First we expand [I], [S], and [F] in Fourier serie 4. For the plane wave case it has been shown all Fourier components of the equation of tran are independent of each other, and for the vertic and horizontally polarized incident waves, the two terms of the stokes vector are even function 4 and the last two terms are odd functions of. ' incoherent specific intensity is given by the sun all Fourier terms, (30) i (7) d I '82]..s in that isfer:ally first.- _ZO Iml m]m Vm /Fl, F21 F31, \F41 A -=-o Fl F21 IF +ira = F31 F41 / Ao [L]ira ( [SI3 ] [s4]a [S4 ]m/ (35) [I] = E [I]m cos m4 + 2 [lim sin m4 m = o m = I We can obtain the equation of transfer for e Fourier component. For the Fourier componen = 0 and a vertically polarized incident wave, d[l]o= -[Io +- fI [L]0[]0 d/' + eI dr [Io t [L]o[l]o d~ + — Ke J -e Ke ' Cp[F+]oe -('i~) + - CpR23[F-]oe-(2'7 - r)l/" Ke I1s01 The equation of transfer given by (32) and (34) The together with the boundary conditions given by (28) n of and (29), constitute the complete mathematical formulation of the problem. An analytical solution of the integrodifferential equation given by (32) and (34) is not available, but it is possible to solve the equation of transfer numerically by several techniques. In this study we use the discrete ordinate:ach method. The details of this technique are given by It m Ishimaru [1978, chap. 11]. Once [1]o and [I]m have been found, the total incoherent specific intensity emerging from the layer at a point above the airsnow boundary can be computed from (32) [I] =[T2,] 2 cos mO tm = 0 K 0

KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I 337 m (i sin m4 r = where [T21] is the transmissivity matrix of the ul boundary. If the incidence angle is larger than critical angle, [Tij] = 0; otherwise, it is given b [Tir] on the wave transmitted into the snow layer because, as we have shown in Table 1, the effective (36) index of refraction of the snow layer is around 1.3. However, we have to account for the direct backscatter from the rough surface, which we can add pper incoherently to the volume backscattering contribut the tion derived in the preceding section. iY Two widely used techniques for analyzing surface scattering are the Kirchhoff approximation and the small perturbation method [Tsang et al., 1985, chap 2]. If a surface has a large rms slope, we can obtain a simple expression using the Kirchhoff stationary-phase approximation. Assuming that the surface characteristics can be modeled by a Gaussian correlation function, we can write the copolar(37) ized backscattering coefficient of the snow surface nj cos ej n 3 cos i tl O i (It0l2 0 0 O 0 Itei2 0 0 O Oth 0 0 Re (t1t.*) Im (t.th*) lU' )i 0 -Im (titi~* Re (t itih*) x' tJ t where tij is given by (21) for vertical polarization and by (27) for horizontal polarization. Once the specific intensity has been computed for a vertically polarized incident wave, the backscattering coefficients of the snow volume for vv and hv polarizations can be obtained from the first and second elements of [I], 01rv = 47r cos 0joI Orv cs61(38) hv = 4ir cos o12 A similar solution is used for computing %hh and vh oTsv 4. SURFACE SCATTERING AT BOUNDARIES In the millimeter-wave region, most natural surfaces, such as ground snow, are rough compared to the wavelength. Surface scattering at the snowground interface is not important for millimeterwave remote sensing because of the large amount of attenuation in the snow layer, however, the scattering at the air-snow interface becomes important when the snow is wet and the volume fraction of ice particles is high. If the top surface is rough, the incident wave transmitted into the snow layer is no longer a plane wave. Hence we need to consider the interaction between surface scattering and volume scattering. This interaction has been examined in the literature [Ulaby et al., 1986, chap 13; Tsang et al., 1985, p. 203] for a Rayleigh layer bounded by a Kirchhoff rough-surface interface. On the basis of these studies we can ignore the influence of the rough-surface aa Ir(0)12 exp [-(tan2 0,/2m2)] aT ss= 2m2 cos4 0o (39) where 00 is the incidence angle, m is the rms slope, and r(0) is the Fresnel reflection coefficient evaluated at normal incidence. The Kirchhoff approximation does not produce depolarization in the backscatter direction. Therefore the total-snow backscattering coefficient is Otr = r Osv + Oass (40a) hh hh 0ts = rsv + Orss hv = orvh = rvh = hv ts ts sv sv (40b) (40c) 5. RESULTS In this section we shall discuss the variation of the backscattering coefficient with snow wetness (liquid water content), snow thickness, and angle of incidence, and we shall evaluate the contribution of surface scattering to the total backscattering coefficient. All calculations were performed for a snow layer with an ice volume fraction of 0.4, containing ice particles characterized by a normal size distribution with a mean diameter of 1 mm and a standard deviation of 0.2 mm. The snow surface is assumed to have a rms slope m = 0.5. Evaluation of the Fourier components of (36) led to the conclusion that it is necessary to include only the first four components in the computation of the scattered intensity I because the contributions of

338 KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 1 0P 35 GHz * Co-pol. + Surface Scattering 10 - Q. Co-pol. - -20 $is _ Cross-pol. o 1 2 3 4 s 6 Liquid Water Content mv (%) Fig. 3. Backscattering coefficients a,, and crts versus m,, at 35 GHz. Symbols are the calculated values based on the radiative transfer theory. Solid lines are third-order polynomial curve fits to calculated results. The incidence angle is 8o = 400 and the incident wave is vertically polarized. The rms slope of the snow surface is m = 0.5. the higher-order components are negligibly small. Moreover, except when the snow layer is dry and its thickness is only a few centimeters, the effects of reflections by the underlying ground surface may be ignored. 5.1. Surface scattering contribution The plots shown in Figure 3 depict the 35-GHz variation of the snow-volume backscattering coefficient a-s with liquid water content for a 0.45-m thick snow layer at an incidence angle 0o = 40~. An additional curve is shown for the copolarized case representing the total-snow backscattering coefficient aots, which includes os and the surface contribution ss,. We observe that os decreases rapidly with increasing liquid water content due to the corresponding increase in background absorption. The rate of decrease of the cross-polarized component is much higher than the rate for the copolarized component, indicating less multiple scattering as the medium becomes highly absorptive. Because we use spherically shaped ice particles in the model, single scattering does not produce depolarization; the depolarized return is caused exclusively by multiple scattering. For dry snow (mv = 0) the contribution of surface scattering crs, is negligibly small in comparison with osv. However, oas increases with increasing m~, and becomes comparable in magnitude with os, at mv = 5%. These observations pertain to 35 GHz. Similar calculations made at 95 and 140 GHz indicate that css << oasv even at mv = 5%. Thus, surface scattering need not be included at these higher frequencies. The variation of the backscattering coefficient with liquid water content at 95 and 140 GHz is shw.vn in Figure 4. Compared with the 35-GHz results, a,, decreases with liquid water content at a much slower rate, and the cross-polarized component remains significant in magnitude even at high liquid water contents, indicating strong volume scattering within the snow layer. Although the background absorption coefficient Kag is approximately twice as large at 95 GHz than its value at 35 GHz (see Table 1), the effect of the background absorption on the total backscattered signal is less important at 95 GHz than at 35 GHz because the scattering cross section of the ice particles is much larger at 95 GHz, as a result of which multiple scattering becomes dominant in comparison to absorption. As was noted earlier, at 95 GHz the effective albedo of the snow volume decreases from 0.99 for dry snow to 0.69 for wet snow with mv = 5%, compared to a corresponding decrease from 0.95 to 0.09 at 35 GHz. At 140 GHz the effective albedo decreases from 0.99 for dry snow to 0.90 for mv = 5%. 5.2. Angular variation The angular variation shown in Figure 5 is very close to a cos O-dependence for both the copolarized and cross-polarized scattering coefficients. The copolarized response represents both hh and vv polarizations because according to the model calculations, there is little difference between the two responses. This cos -like variation was also observed at 95 and 140 GHz. 5.3. Variation with snow thickness Figure 6 shows the response of the backscattering coefficient to snow thickness at 35 and 95 GHz for 80 = 40~. Two wetness conditions are shown. In all cases the minimum thickness considered in the model calculations was 5 cm. For dry snow, av, increases rapidly with increasing snow thickness, particularly for cross polarization, until it ap

KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I 339 4).0 l) a) 0 z W) U ft a) m cr 80 4" 0.0 CO 62 cn.w 95 GHz Co-pol. 'a- A_~~~~~~~~~~ to v i.w 0 C S U co U S a U 0 J9 C) a 0 -5 -10 -15' -20 35 GHz Co-pol. Cross-pol. m 0%v Co-pol. V Cross-pol. *h age mvml% -101 Cross-pol -20..... 0 1 2 3 4 5 Liquid Water Content mv (%) 10.. 0 10 20 30 40 50 60 70 80 Incident Angle (degrees) Fig. 5. Backscattering coefficient cr, versus incidence angle at 35 GHz for a 0.45-m-thick layer. the wetness need not be uniform in depth because the radar response is essentially controlled by the wetness of the very surface layer (except when the surface layer is dry and the snow contains a subsurface layer with nonzero wetness). C o c 0 O.a C. U at an 140 GHz Co-pol. Cross-pol. -10 co C c 0 a o o c CI Uy 0 --10 --20 35 GHz Co-pol. o f rol$-pol. _ ) my VO%,I Co-pol. 4 my.2% Cross-pol. V e _ + ~~~~. 0 1 2 3 4 a 0 1 2 3 4 5 Liquid Water Content mv (%) -1 I..... 0 20 40 60 80 100 12 Snow Thickness (cm) 0 Fig. 4. Backscattering coefficient a,, versus my at 95 and 140 GHz. Symbols are the calculated values based on the radiative transfer theory. Solid lines are third-order polynomial curve fits to calculated results. The incident angle is 0o = 40~ and the incident wave is horizontally polarized... -u C 0 'a U 8 U a C U op 00 a S U U U J9 LI a 0 95 GHz I Co-pol. Cross-po.m 0% Cpo-pol..5 -10 proaches a saturation level beyond which the increase becomes very gradual. 5.4. Simulation of diurnal response Let us consider a diurnal cycle during which the liquid water content of a 0.45-m thick snow layer exhibits the variation shown in Figure 7a, which is Gaussian shaped with a peak value of 2%. In fact, 0 20 40 60 80 1 Snow Thickness (cm) 00 120 Fig. 6. Backscattering coefficient o,, versus snow thickness. The incident angle is 0o = 400 and the incident wave is vertically polarized.

340 KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I 3 E E 0 0 0 c 0 m 2 1 6 8 10 12 14 16 18 Time (hours) (a) 20 22 co Q) 4-0 c 0 0 cn o 0 O m 0 0 o 0 0 10 0 -20 --20. 6 8 10 12 14 16 18 Time (hours) (b) 20 22 35 GHz Co-pol. Cross-pol *f O,,, *" 0' 4" a Ar asf is,4 A A to C c 0 0 0 CD C Q) 0 10 12 14 1s Time (hours) (c) 10 12 14 16 18 Time (hours) (d) Fig. 7. (a) Model for the liquid water content variation in time. The coefficient m,, = Wpeak exp [-(t - 13.85)2/a2 where Wpeak = 2% is the peak value, t is time in hours and ao = 2 is the standard deviation. (b), (c), and (d) Calculated diurnal variations of the backscattering coefficient o,,s at 35, 95, and 140 GHz. Figures 7b, 7c, and 7d show the corresponding diurnal radar response at 35, 95, and 140 GHz. All cases include both the copolarized and cross-polarized responses. As expected, the radar diurnal responses are approximately mirror images of the liquid-water diurnal variation, with 35 GHz exhibiting the greatest dynamic range and 140 GHz exhibiting the smallest. 6. CONCLUSIONS Using radiative transfer theory, we developed a reasonably and computationally efficient, uncomplicated, model for relating the backscattering coefficient of snow at millimeter wavelengths to the physical properties of the snow layer. According to the model, snow-surface roughness is unimportant when the snow is dry, but when the snow is wet

KUGA ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, I 341 snow-surface roughness is important at 35 GHz but not at the higher frequencies examined in this paper (95 and 140 GHz). Within an error of about 1 dB the backscattering coefficient or, varies with incidence angle as cos 0o for both the copolarized and cross-polarized configurations at all millimeter-wave frequencies and liquid water contents considered in this study. The effective penetration depth of the snow medium is strongly dependent on frequency and the liquid water content. At 35 GHz the effective penetration depth decreases from about 30 cm for dry snow down to about 1.5 cm for m,a = 5%. At 95 and 140 GHz the penetration depth of dry snow is 5 and I cm, respectively. The backscattering coefficient decreases with increasing liquid water content, with 35 GHz exhibiting the strongest sensitivity to liquid water variations and 140 GHz exhibiting the weakest. Acknowledgment. This work was supported by U.S. Army Research office contract DAAL03-89-K-0056. REFERENCES Baars, E. P., and H. Essen, Millimeter-wave backscatter measurements on snow-covered terrain, IEEE Trans. Geosci. Remote Sens., 26, 282-299, May 1988. Cheung, R. L.-T. and A. Ishimaru, Transmission, backscattering, and depolarization of waves in randomly distributed spherical particles, Appl. Opt., 21, 3792-3798, 1982. Currie, N. C., J. D. Echard, M. J. Gary, A. H. Green, T. L. Lane, and J. M. Trostel, Millimeter-wave measurements and analysis of snow-covered ground, IEEE Trans. Geosci. Remote Sens., 26, 307-318, 1988. Hallikainen, M., Scattering properties of snow in the 10 to 90 GHz range, Proc. ISAP, 667-670, 1985. Hallikainen, M., F. T. Ulaby, and M. Abdelrazik, Dielectric properties of snow in the 3 to 37 GHz range, IEEE Trans. Antennas Propag., 34, 1329-1340, 1986. Hayes, D. T., U. H. W. Lammers, and R. A. Marr, Scattering from snow backgrounds at 35, 98, and 140 GHz, Tech. Rep. RADX-TR-84-69, Rome Air Dev. Cent., N. Y., 1984. Ishimaru, A., Wave Propagation and Scattering in Random Media, vols. I and II, Academic, San Diego, Calif., 1978. Ishimaru, A., and Y. Kuga, Attenuation constant of coherent field in a dense distribution of particles, J. Opt. Soc. Am., 72, 1317-1320, 1982. lshimaru, A.. R. Woo, J. W. Armstrong, and D. C. Blackman, Multiple scattering calculations of rain effects, Radio Sci., 17. 1425-1433, 1982. Kuga, Y., and A. Ishimaru, Retroreflectance from a dense distribution of spherical particles, J. Opt. Soc. Am., A Opt. Image Sci., 1, 831-835, 1984. Ma, R., A. Ishimaru, and Y. Kuga, Scattering and depolarization of waves incident upon a slab of random medium with refractive index different from that of the surrounding media, Radio Sci., 25, 419-426, 1990. Matzler, C., and U. Wegmuller. Dielectric properties of fresh water ice at microwave frequencies, J. Phvs. D Appl. Phys., 20, 623-630, 1987. Narayanan, R. M., and R. E. McIntosh, Millimeter-wave backscatter characteristics of multilayered snow surfaces, IEEE Trans. Antennas Propag., 38, 693-703, 1990. Stiles, W. H., and F. T. Ulaby, The active and passive microwave response to snow parameters, I, Wetness, J. Geophys. Res., 85, 1037-1044, 1980. Tsang, L., and A. Ishimaru, Backscattering enhancement of random discrete scatterers, J. Opt. Soc. Am. A Opt. Image Sci., 1, 836-839, 1984. Tsang, L., and A. Ishimaru, Radiative wave equations for vector electromagnetic propagation in dense nontenuous media, J. Electromagn. Waves Appl., 1, 52-72, 1987. Tsang, L., J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing, John Wiley, New York, 1985. Ulaby, F. T., and W. H. Stiles, The active and passive microwave response to snow parameters, II, Water equivalent of dry snow, J. Geophys. Res., 85, 1045-1049, 1980. Ulaby, F. T., R. K. Moore, and A. K. Fung, Microwave Remote Sensing, vol. 11I, Artech House, Norwood, Mass., 1986. Ulaby, F. T., T. F. Haddock, and R. J. Austin, Fluctuation statistics of millimeter-wave scattering from distributed targets, IEEE Trans. Geosci. Remote Sens., 26, 268-281, 1988. Ulaby, F. T., R. Austin, T. Haddock, and Y. Kuga, Millimeterwave radar scattering from snow, II, Comparison of theory with experiment, Radio Sci., this issue. Wen, B., L. Tsang, D. P. Winebrenner, and A. Ishimaru, Dense medium radiative transfer theory: Comparison with experiment and application to microwave remote sensing and polarimetry, IEEE Trans. Geosci. Remote Sens., 28, 46-59, 1990. Williams, L. D., J. G. Gallagher, D. E. Sugden, and R. V. Birnie, Surface snow properties effects on millimeter-wave backscatter, IEEE Trans. Geosci. Remote Sens., 26, 300-306, May 1988. R. D. DeRoo, T. F. Haddock, Y. Kuga, and F. T. Ulaby, Radiation Laboratory, Department of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109.

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

344 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 measured using the freezing calorimeter technique, and (7) microscope photographs of thin snow samples, from which the ice crystal size is estimated. The standard measurement procedure consisted of measuring the backscattered power for hh, hv, vh, and vv polarizations at each of 30 or more spatial locations. The backscattering coefficient was obtained by averaging the measurements for the different spatial locations. In addition to spatial averaging, frequency averaging was used to further improve measurement precision [Ulaby et al., 1988b]. The estimated uncertainty associated with the values reported in this paper is +0.5 dB. Analysis of t,. data shows that the hv and vh measurements are essentially identical (within a fraction of 1 dB), which is expected from the reciptwcity relation. In almost all cases the copolarized esponses, 0h and,, were within 1-2 dB of each other. 0Uhh an07wrewti -2d feahohr 2. ANGULAR RESPONSE The data shown in Figure 1 were measured for a 12-cm-thick layer of dry, freshly fallen, unmetamorphosed snow composed of ice crystals with diameters on the order of 1 mm. The measured rms height was 1.4 mm. Only hv- and vv-polarized data are shown because the difference between Oh and ory is I dB or less across the entire angular range at all three frequencies. The curves shown in the figure were calculated according to the theoretical model described in the preceding paper. For vv polarization, theory and experiment are in good agreement at 35 and 95 GHz. At 140 GHz, however, the level predicted by theory is lower than the experimental observations for vv polarization by about 4 dB. We attribute the difference to the backscattering enhancement effect, which the model does not take into account. A typical example of the angular dependence for wet snow is shown in Figure 2. The theoretical curves were computed assuming a mean ice crystal diameter of 1 mm and a rms slope of 0.07. Reasonable agreement between theory and experiment is obtained except for vv polarization of 140 GHz; we again attribute the difference to the backscattering enhancement effect, although no strong evidence exists to support this contention. 3. EFFECT OF SNOW SURFACE ROUGHNESS According to the model results presented in the previous paper, snow surface roughness should 10. 10 0*0 U la - L ci 'E Q (A.W u:2 0..10. -VV Thwsa.HV Tlwov vv DM ~ VVDWa a KV DM o a a 0~~~~ 0 0..................................................... ~~~ -. "........ f = 35 GHz -. 1 --- I.... - --......... i.... -20. 0. 10. 20. 30. 40. 50. 60. 70. Incident Angle e0 ( degrees ) 10. co la 0 ~1 -@10 C I u rt 2 0. vv Thw7......... Thm/ VV VDm 0 HV Dw "'~~~~~~I.~~~~~~~~~~~. f= 95 GHz -10. -20. 10. 0. 10. 20. 30. 40. 50. 60. 70. Incident Angle Oo (degrees ) A I 100 4.I jO '! mi 0. 0 0 * 0 4~~0 I~~~~~~~~~ ~ ~~~~~~............lw e*...... -..... HV Thwy ~WDI fa 140GHz ~ HV D.,ii t...i.. -10. 0. 10. 20. 30. 40. 50. 60. 70. Incident Angle S( degrees ) Fig. 1. Measured and calculated backscattering coefficient for a dry snowpack with the following parameters: depth, 12 cm; snow density, 0.2 g/cm3; mean crystal diameter, 1 mm; and rms surface slope, 0.07. have a minor effect on the level of ar (except at normal incidence) when the snow is dry. For wet snow, however, o- should increase by as much as 5 dB at 35 GHz if the surface is made rough relative to the wavelength, but the increase should not be significant at 95 GHz or higher frequencies. There is

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

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

ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 347 5. 4. 3. 2. 1. 0. on._ 04 u.E as m % t;= u mn 00 0. -2. -4. -6. Snow Wetness 6. 8. 10. 12. 14. 16. 18. 20. 2 o Smooth... -. Sligtly Rough M-, -- VerRoh,, Ao~ ~ ~~~'. a ~. 8... 12. 14..,.. 2. 6. 8. 10. 12. 14. 16. 18. 20. 22. 10. 0. b Liquid Water Content (o),", '. -.., \ Air Temerature ('C) 1 1...........Tempa ------ Liquid WsarConmu ni 8. 10. 12. 14. 16. 18. 20. 22. 24. -10. Time of Day ( hour ) Fig. 5. Measured air temperature, snow surface temperature, and liquid water content of the top 5-cm layer on March 1, 1990. -8. -10. -12. 22. 2. C 0. -2. I -4. Fig. 4. content of t polarizatior diameter of rms heights 0.89 cm (sli has a uniform depth profile between the snow and ground surfaces. During the melting phase of the em< --- —_- diurnal cycle the melting process starts at the air-A;. ~snow interface because the source of thermal energy is the warmer air mass above the snow pack or 95GHz VH * direct solar radiation (under cloud-free conditions). [?O, I,"-: — -_ L,.,.,For the experiment under discussion the air temper6. 8. 10. 12. 14. 16. 18. 20. 22. ature rose from -8~C at 0800 LT to a high of +6~C Time of Day (hour) at 1145 LT, and then decreased down to -3~C at 2000 LT. As the snow surface layer starts to melt, Observed diurnal variation of in, the liquid water its thermal conductivity increases, thereby allowing the top 5-cm layer, and oO at 35 and 95 GHz for cross the transfer of thermal energy down to lower layers. a. The snowpack was 28 cm deep, had a mean crystal During the melting phase, mv is highest at the f 0.6 mm, and a density of 0.4 g/cm3. The measured surface and decreases in an exponentiallike manner s of the surfaces were: sl = 0.49 cm (smooth), s2 with depth. At millimeter wavelengths the penetraightly rough), and s3 = 1.98 cm (very rough). tion depth in dry snow is of the order of 1 m at 35 GHz, 5 cm at 95 GHz, and 1 cm at 140 GHz; and for mv = 3%, the penetration depth is only 3 cm at 35, the temporal variation is similar to that GHz and less than 1 cm at the higher frequencies. I by the other plots. Hence, when the top snow layer is wet, it alone Ily different situation is observed in Figure governs the radar backscatter, and the lower layers period from 1300 to 1900 LT, during which exercise an inconsequential influence on a~. 'pack undergoes a refreezing process. In A different situation occurs during the refreezing.n the measured mv is at its peak value of phase. As the air temperature starts to decrease, we pproximately 1400 LT), all three radar- have a reversal in energy balance with the snowI values of my are in the 0.8-1.5% range. pack becoming the source of excess energy. Freezthe models, as used thus far, are inappro- ing starts at the surface and proceeds slowly downr modeling the backscatter from wet snow ward to lower layers. This results in an inverted ese "refreezing" conditions. This lack of profile with mv = 0 at the surface, increasing slowly ndence between theory and measurement to a maximum value at some depth below the I by the fundamental assumption that mv surface, and then decreasing again as the depth is channels, exhibited A total 8 for the I the snow fact, whe 7% (at a predicted Clearly, 1 priate for under the correspoi is caused

348 ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 10. 35 GHz od 10 00.. o Zq U to. _ 2 a CP m 0. I, WP.I,, ~,,, p '-,.,. o.........HHolbd - 0 Q 00 I "3 C-) zo. _-.: Bo. _ 13 -10. an -20. 8. 10. 12. 14. 16. 18. 20. 22. 24. 10. - I ' I - I I.. -I, — _ 0 c 0.o Im QP 0., V Polizam 95 GHzHpoim. i. i. I. I.!. |i.... -10. F c0 00 D. () E._ (0 me o0 8 Go -i I.) c tO c (0 0. -10. -20. -30. 0. -10. -20. 0. 0. 2. 4. 6. 8. 10. 12. i I I " I I I......... 95 GHz FwrtOrder 95 GH Numerica 0. 2. 4. 6. 8. 10. 12. 4.u. -20. _ 8 10. r. 10. 12. 14. 16. 18. 20. 22. 24. I I " I I I ' I.I I E a l1 CS 0..... VV Pohmimim........ VH Poemm 140GHz - HH j -10. ' ' ' I ' I ', ' 140 GHz Fits Orde 140 GHz Numericl..................,.................... 0. 2. 4. 6. 8. 10. 12. -lo0. -20. Volumetric Liquid Water Content m, ( % ) 8. 10. 1 14. 16. 18. 20. 22. 24. Time of Day ( hour ) Fig. 6. Backscattering variation measured at 35, 95, and 140 GHz on March 1, 1990. The incidence angle was 40~. The snowpack density was 0.32 g/cm3, depth was 12 cm, and other parameters are given in Figure 5. increased further. Figure 9 shows a family of profiles generated on the basis of this simple logic, constrained by the condition that the value of m,, averaged over the 5-cm layer has to be equal to the value measured by the freezing calorimeter. With the uppermost surface layer being dry in the inverted-profile case, the radar response will be governed not only by the dry surface layer, but by the wet layer immediately underneath it as well. It is not surprising, therefore, that the model-predicted val Fig. 7. Calculated backscattering coefficient for a deep snowpack versus liquid water content for hh polarization at an incidence angle of 40~. The calculations are based on the firstorder and the numerical solutions of the radiative transfer equation for a snowpack with the following properties: mean crystal diameter, 1 mm; and snow density, 0.32 g/cm3. ues of mv are quite different from the 5-cm average measured in the field. To solve this problem, we propose to use a hybrid first-order numerical model for characterizing the radar response of snow. 4.1. Hybrid first-order-numerical model Let us consider the advantages and limitations of the first-order solution and the numerical solution of the radiative transfer equation as we apply them to the snow problem. The first-order solution takes

ULABY ET AL.: MILLIMETER-WAVE RADAR SCATTERING FROM SNOW, 2 349 8. 7. e 6. S 5. 3 2. a. 1. 0. -1. Fig. 8. labeled "m layer as me plots refer observatior 0o = 40~). into acc( volume, account quently, the magn ing albed with incr increasec ("wet" ' 10. 9. E 8. 7. C_ 6. 5. 4. *3~ 3. 2. 0. Estimations of Liquid Water Content refraction of the ice crystal to that of the back-,.,.,,,,,, ground decreases with increasing snow wetness. To compare the first-order solution with the numerical /5.G, s --- 39asolution, we refer to the plots in Figure 7 which 95 show 0hh versus m0 for an incidence angle of 40~. \. I4QGH, For dry snow (mv = 0) the first-order solution underestimates cr0 by about 3 dB at 35 GHz, 8 dB at 95 GHz, and 11 dB at 140 GHz. The differences between the two solutions decrease with increasing; I N.~ \1 \ 4 mv, becoming insignificant at 35 GHz for m. - 2%, but not so at the higher frequencies. Another major 'i ~ I.,.,.,.,.,. advantage of the numerical solution is that it pro08. l. 12. 14. 16. 18. 20. 22. 24. vides a result for the cross-polarized scattering coefficient aOrv, whereas the first-order solution does Time of Day ( hour) not predict cross polarization because the scattering Diurnal variation of liquid water content. The plot particles are assumed to be spherical in shape. leasured" refers to the value of m,, of the top 5-cm Having stated the advantages of the numerical easured by the freezing calorimeter. The other three solution over the first-order solution, let us now to mV as predicted by the model using the radar consider the converse. The numerical solution reIs at the noted frequencies (with hh polarization and quires a great deal of computation, in spite of the fact that the snow medium is assumed to have uniform properties throughout. If we are to formuount only single scattering in the snow late the solution for a multilayer structure in order whereas the numerical solution takes into to accommodate a nonuniform depth profile for m, all orders of multiple scattering. Conse- the complexity of the numerical approach would the first-order solution may underestimate make the solution computationally impractical. On uitude of 0.0, particularly when the scatterlitude of, particularl wheno is large. The scattering albedo decreases - the other hand, the first-order solution is perfectly lo is large. The scattering albedo decreases easing liquid water content because of the amenable to computing the backscatter from a ~easing liquid water content because of the medium with nonuniform properties in the depth i absorption by the background material dimension. Hence, we propose to use a hybrid iir) and because the ratio of the index of dimension. Hence, we propose to use a hybrid model that takes advantage of the more accurate feature of the numerical solution and the easier Moisture Profile structure of the first-order solution. The procedure....,....,....,...,..,..,,..,..,. is as follow s. -1:1:30 For a snowpack with specified density and crystal.. ----.12:00 size distribution, we compute a~ for the semi------ 13:00......14:30 l infinitely deep case (at the angle of incidence, -: *'-.1 6:30 frequency, and polarization of interest) as a func"-. '....7:30 tion of mv (assuming a uniform profile with depth) --—. 18:15 using the numerical solution outlined in the preced-... IRI 5 ing paper. We shall denote this backscattering co-;.v..- 230 - efficient O.~(0o, f, pq), where p and q denote the polarization configuration of the receive and trans4..... - mit antennas, respectively. Figure 10 shows the results of such a calculation at 35, 95, and 140 GHz 0. 1. 2- 3- 4- 5. 6. 7. 8. 9. 10. for 00 = 40~ and pq = hh. The form of the variation Depth ( cm) can be fitted to the functional form: Fig. 9. Proposed depth profiles for my,. The pack is assumed to be totally dry both prior to 1100 LT and after 2200 LT. 10 log [a (my)] = ao + alm1/2 + a2m, (1) n v~~~~~~~~~~~~~

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

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

Radio Science, Volume 25, Number 1, Pages 9-18, January-February 1990 Millimeter wave scattering model for a leaf K. Sarabandi, F. T. Ulaby, and T. B. A. Senior Radiation Laboratory, Electrical Engineering and Computer Science Department, University of Michigan, Ann Arbor (Received February 23, 1989; revised November 6, 1989; accepted November 7, 1989.) At millimeter wave frequencies a typical leaf is a significant fraction of a wavelength in thickness, and its nonuniform dielectric profile now affects the scattering. To provide a simple and efficient method for predicting the scattering, two types of physical optics approximations are examined. The first approximates the volume polarization current by the current which would exist in an infinite dielectric slab with the same profile, while the second (and simpler) one employs the surface current which, on the infinite slab, produces the known reflected field. It is shown that the first method is superior, and provided the actual dielectric profile is used, it predicts the scattered field to an accuracy which is adequate for most practical purposes. 1. INTRODUCTION Leaves are a key feature of any vegetation canopy, and in order to model the scattering from vegetation-covered land, it is necessary to develop an efficient and effective technique for predicting the scattering from a single leaf. At microwave frequencies where a typical leaf is electrically thin with lateral dimensions at least comparable to the free space wavelength A0, several methods have been proposed [e.g., Le Vine et al., 1985; Willis et al., 1988] all based on the physical optics approximation applied to a uniform dielectric slab. In particular, if the leaf thickness is no more than about A0/50, physical optics in conjunction with a resistive sheet model predicts the scattering at most angles of incidence [Senior et al., 1987] and can also handle curved leaves [Sarabandi et al., 1988]. On the other hand, at millimeter wavelengths the thickness can be a significant fraction of a wavelength, and it is also necessary to take into account the internal structure of a leaf. At least two different types of cell can be distinguished, and their differing water content affects the dielectric constant, leading to a nonuniform dielectric profile. To compute the scattering at these higher frequencies, two different physical optics approximations are examined. The first of these employs the polarization current which would exist in an infinite slab consisting of one, two, or more layers simulating the dielectric profile of the leaf, and this is referred to as Copyright 1990 by the American Geophysical Union. Paper number 89RS03459. 0048-6604/90/89RS-03459$08.00 the volume integral physical optics (VIPO) approximation. When there are many layers, a convenient method of implementation is described in the appendix. In the second (and simpler) approach the equivalent surface currents are approximated by the electric and magnetic surface currents of an infinite slab. New expressions for the physical optics currents, which are more convenient than the standard ones [Beckmann, 1965], are introduced, and this technique is referred to as the surface current physical optics (SCPO) approximation. For an electrically thin leaf or plate the two approximations are indistinguishable, but as the thickness (or frequency) increases, the predicted scattering differs in most directions, and by comparison with the results of a moment method solution of the volume integral equation it is shown that VIPO is superior. In addition, for a two-layer material it is no longer adequate to treat the plate as a homogeneous one having an average dielectric constant. Provided the actual dielectric profile of a leaf is simulated, it appears that VIPO can predict the scattering behavior of a leaf to an accuracy that is sufficient for most practical purposes at millimeter wavelengths. 2. STRUCTURE OF A LEAF The structure of a typical vegetation leaf is shown in Figure 1. The type and number density of cells may vary as a function of depth into the leaf which, in turn, results in a nonuniform dielectric profile. The effect of this nonuniformity becomes observable at higher frequencies where the thickness of the leaf is comparable to the wavelength. 9

10 SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF Upper Cuice Palisade Layer Spongy 0 0 0 0 00 t Layer Oo,, o Lr COicle Fig. 1. The structure of a typical vegetation leaf. Leaves contain two types of photosynthetic cells: palisade parenchyma, consisting of column-shaped cells in which most photosynthesis takes place, and spongy parenchyma, consisting of irregularly shaped cells with large spaces between them [Curtis and Barnes, 1985]. Because a large part of the vegetation material is water, its dielectric constant is strongly influenced by the dielectric constant of water and the water content. For most leaves the water content is higher in the upper layer (palisade region) than in the undersurface (spongy region). The sensitivity of the dielectric constant to water content is much greater in the lower part of the millimeter wave spectrum than in the upper, but this is more than counterbalanced by the thickness to wavelength ratio. The net result is that the sensitivity to dielectric variations is greater at the higher frequencies. To examine the effect of the nonuniform dielectric profile on the scattering properties of the leaf at millimeter wavelengths, we computed the normal incidence reflection coefficient ro of a two-layer dielectric slab and compared it with the reflection coefficient of a uniform dielectric slab whose dielectric constant is the average. The computation was performed for a leaf thickness of 0.5 mm, and the water content ratio of the two layers was chosen to be 4 to 1, representing a marked variation between the upper and lower surfaces of the leaf. From the data in Table 1 it is seen that when the two-layer slab is approximated by a uniform slab the error in the reflection coefficient increases with increasing frequency and is as large as 4 dB at 140 GHz. 3. PHYSICAL OPTICS APPROXIMATIONS At microwave frequencies where a typical leaf is no more than about AO/50 in thickness with lateral dimensions comparable to or larger than the wavelength, the scattering properties can be accurately predicted using the physical optics approximation applied to a resistive sheet model of a leaf [Sarabandi et al., 1988]. In effect, the leaf is modeled as an infinitesimally thin layer, but as the frequency increases, it is necessary to take the leaf thickness into account. There are now two types of physical optics approximation that can be employed. One is the surface current (SCPO) approach in which an infinite dielectric slab is replaced by an equivalent sheet current that produces a plane wave identical to the reflected wave of the slab. This current is then used as an approximation to the equivalent surface current over the upper surface of a finite dielectric plate. Alternatively, the induced (volume) polarization current in the plate can be approximated by the current in the infinite dielectric slab, and we shall refer to this as the volume integral physical optics (VIPO) method. It is more accurate than the SCPO method, although the latter is more convenient to use for evaluating the scattered field. To illustrate the two procedures, consider a dielectric plate consisting of a homogeneous dielectric of thickness dl and relative permittivity E1 atop a second material of thickness d2 - dl and relative permittivity E2. The plate occupies the region -a/2 s x - a/2, -b/2 < y - b/2, and -d2 - z - 0, as shown in Figure 2, and is illuminated by an Epolarized plane wave whose electric vector is Ei = yeiko(xsinGo - zcos0o) (1) where ko is the propagation constant in the free space medium above and below the plate. When the TABLE 1. Voltage Reflection Coefficient for a Two-Layer and Average Dielectric Slab El + E2 If, GHz El EE(avg) 2 Fo Fo(avg) 35 20 + i21 6 + t3 13 + i12 0.74L6 0.78L-0.16 94 6 + i5 ' 2 + il 4 + t3 0.59L 12 0.48L27 140 5 + i4 2 + il 3 + t2.5 0.50L20 0.34L26.1

SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF 11 aRoo 0- 8 Region 0 T The corresponding results for a single layer of thickness dl and relative dielectric constant E1 can be obtained by putting d2 = dl and k2z = kLz, implying B2 = B1 and A2 = A 1. Given a volume distribution of electric current J in free space, the corresponding Hertz vector is a 2 x K e ma i' X. XX1X X NNN\\\NX\\\\\NIyI JO —-,r-....\\\\\\\\\\\\\\~ ' s\\\\\\\\\\\\\~ l~,ll. q Regio 3 II() = 4rko J(r) Ir-Ft dv v r- (6) Fig. 2. The geometry of the scattering of a plane wave from a two-layer dielectric slab. plate is treated as an infinitely extended slab, the electric field can be written as where ZO(= 1/Yo) is the free space impedance, and the resulting field is E(r)=V x V x H(r) H() = -ikoYoV x H(r) In the far zone of the current distribution, Ey = (e -ikoz + reikozZ)eikosin0sx O < z Ey = (Ble-iklz +A eiketz)eikOs'nofx -dl sz < 0 Ey = (B2e-ik2zz + A2eik2z)eikOsin0ox -d2 <z < -d (2) eikor iZ0 r' -ikorFC dv I(rD) - - J(QF)-/kor dv kor 4ir (7) Ey = B3e -ikoZeikosin0Ox z < -d2 E(r) =-ko2f x f x n(r) (8) where koz = ko cos 0o kjz = ko(Ej - sin2 0) 1/2 for j = 1, 2. If R1 and R2 are the reflection coefficients at the upper and lower surfaces where koz - ktz koz - k2z R1=k+k R2=k...+k.. and C.+ = 1 + {1 + R2e2i12'(d: -di} {1 - R2e2ik2:(d2 -d,1-1 application of the boundary condition at the three interfaces gives In the dielectric slab the volume current J is the polarization current J = -ikoYo(Ej - 1)Eyy (9) where Ey has the value appropriate to each layer (j = 1, 2), and when this is inserted into (6) and the integration carried out over the volume occupied by the plate, we obtain the VIPO approximation. For scattering in the direction Os indicated in Figure 2 the expression for the Hertz vector is eikor koab sin X = kor 4:r X F (10) where C + (1 +R1) B1 =C + + C R le2 td' A ei(ku - k2z)dl 2 B2 = 1 _ R2e2k2(d2 -d,) C + _ 2ikdB1 A2 = -R2e2ik2zd2 1 - e -i(kz - kocosO,)d A (3) F=(E — 1) i(kkz-ko cos Os) A 1 ei(khz + kocsOw)di } B2 -i(klz +ko cos Os) (4) e -i(k2z - kocoss )dl e -i(kz - kocoss )d2 + (E2 - 1) i(k2 - ko cos ) A2 i~2z -k o s B3 = ei(k2z - kozW)2(1 - R2)B2 and C+R1 + C _ e2ikd' '.C+ + CRleZ likzd ei(k2z + kocsOs)dl _ ei(k2z + kocosO,)d2 (5) - i(k2z + ko cos 8) B2 (11)

12 SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF koa X= - (sin es + sin 80) (12) 2 The far zone scattered field can then be obtained from (7) and written as e ikor Es= - SE(Os, 00) (13) kor where SE(O8, 00) is the far-field amplitude, and for the VIPO approximation the result is k3ab sin X SEIP~(, 00) = - F (14) In terms of the far-field amplitude the bistatic scattering cross section is o(Os, 00) =- IS(O, 00)12 (15) The more conventional SCPO approximation can be obtained by noting that the electric current sheet J = -2Yo cos 0oFeikOsinsoxS(Z)9 (16) produces a plane wave identical to the field reflected from the dielectric slab. As evident from the impulse function 6(z) in (16), the current is located at the upper surface of the slab, and when (16) is inserted into (7) we find where E = (ikoej)-lZoaHy,/x and Ez = -(ikoej)- ZoaHyIa x have the values appropriate to each layer (j = 1, 2). The Hertz vector can be computed using (6), and for the scattered field Hs the far-field amplitude is found to be k3ab sin X SVIP~(os, 80) = 4 r X-(cos OsF' - sin OsFj) (21) where F klz(e - 1) 1- e -i(ki, - kccos8,)d, Fj= =IAl ko el i(klZ - ko cos Es) 1 - ei(klz + kocosO,)dl ] k2z(E2 - 1) i(klz + ko cos Es) + ke 2 e -i(k2z - kocoss,)dl _ e-i(ki, - kocoses)d2 A- k cos i(k2z - ko cos 0s) ei(k2 + kocosO,)dl _ ei(k2 + kocosOs)d2 + i(k2z + ko cos Os) B2 l F2= sin 00 {( - 1) E I (22) [1- e -i(kl - kocos6O)dl. i(klz - ko cos 8s) 1 - ei(klz + kocosO,)dl i(klz + ko cos 08) Bf eikor -i sin X HnscPo~ y T 2 1 cos OoFab X and the far-field amplitude is then SCPO( = ~-ik o sin X SE Os, 00) = cos OoFab X (17) (E2 - 1) e -i(k2z -kocos, )dl -_ e -i(kit -kocos6,)d2 + i(k -ko cos E2 i(k2 -ko cos 86) 2 ei(k2z + kocossO)di - ei(k2 + kocosO,)d2 (18) i(k2z + ko cos 0s) BjI} (23) In the specular (08 = - 80) and backscattering (0s = 00) directions it can be verified that (14) and (18) are identical, but in the other directions the two approximations differ. In the case of H polarization for which Hi = geiko(xsin"o - zcOSBo) (19) the analysis is similar. With Hy represented as shown in (2) the various coefficients (now indicated by primes) differ from those for E polarization in having klz replaced by klz/el and k2z replaced by k2z/E2 everywhere except in the exponents. The induced polarization current then has two components and is given by The SCPO approximation can also be obtained by noting that a magnetic current sheet of the form J* = -2Zo cos oor'eikosinOoxs(z)y (24) generates a plane wave identical to the reflected wave. When we use this as the equivalent surface current on the dielectric plate, the magnetic far-field amplitude becomes -ik2 sin X sSCPO (0s, 00) = - cos 00F'ab SH 2~rr x (25) As in the case of E polarization, the two approximations are identical in the specular direction, but (21) and (25) differ in all other directions, including backscattering (0s = 80) unless 00 = 0. (20) J = -ikoZo(ej - 1)(Exi + EI)

SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF 13 0.............................. -90 -70 58 -38 -18 8 38 so 78 90 artt elng Anglo (Do.e') Fig. 3. Amplitude of the ratio of the bistatic far-field amplitude of VIPO to SCPO for E polarization of a dielectric plate with d2 = Ao/4 and El = e2 = 3 + iO.1 at 00 = 300. 4. NUMERICAL RESULTS To illustrate the difference between the VIPO and SCPO approximation, we consider a homogeneous (single layer) plate of thickness d2 = Ao/4 with E2 = 61 = 3 + i0.1. For an E-polarized plane wave incident at 300 the amplitude and phase of SVIPO/ SSCPO are given in Figures 3 and 4, and these show that the difference increases away from the specular and backscattering directions. At a fixed scattering angle the difference increases with the electrical thickness of the plate up to the first resonance and then decreases. To test their accuracy, the two approximations have been compared with the results of a moment method solution of the volume integral equation. The particular code used is a two-dimensional one which was extended to three dimensions by assuming that the induced currents are independent of the y coordinate. Since the dielectric constant of most vegetation materials is high, it is necessary to have the cell sizes very small, and one consequence of this is the need to compute the matrix elements extremely accurately, especially for H polarization [Sarabandi, 1989]. For a 2Ao square plate formed from the above mentioned layer and illuminated by an E-polarized plane wave at normal incidence, the two approximations are compared with the moment method solution in Figure 5, and the superiority of VIPO is clear. In the case of a thin plate the two approximations are indistinguishable. This is illustrated in Figure 6 showing the VIPO expression (14) and the moment method solution for a 2A0 square plate of thickness -39 -1s le 3X noattring Angle (Doroee) Fig. 4. Phase of the ratio of the bistatic far-field amplitude of VIPO to SCPO for E polarization of a dielectric plate with d2 = A0/4 and EI = E2 = 3 + iO.1 at 00 = 300. d2 = A0/50 for E polarization. The plate is a homogeneous one having E = 13 + i12 corresponding to the average permittivity at 35 GHz in Table 1. The SCPO expression (18) yields the same results, as does a two-layer model having the permittivities listed in Table 1. The analogous data for H polarization are given in Figure 7, and over a wide range of scattering angles the approximate and moment method solutions are in excellent agreement. As the frequency and hence the electrical thickness of the plate increase, the superiority of the VIPO approximation becomes apparent and, in addition, it becomes necessary to take the layering -48.M............................................. -58 -90 -7 -3 - s 3 s W ii -Gu Satten Angl (Deree.) Fig. 5. The bistatic cross section of a 2Ao x 2A0 plate for E polarization with d2 = Ao0/4 and E 1 = 2= 3 + iO.l at normal incidence: moment method solution is denoted by the solid line, VIPO by the short-dashed line, and SCPO by the long-dashed line.

14 SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MObEL FOR A LEAF -10 -SO / i i i i.........: s90 -70 -60 -3 -18 1s 3 so0 70 90 a ettaw Angle (Doure) Fig. 6, The bistatic cross-section area of a 2Ao x 2AO plate for E polarization with d2 = AO/50 and Eavs = 13 + i12 at normal incidence: moment method solution is denoted by the solid line. VIPO or SCPO is denoted by the dashed line. of the plate into account. In Figures 8 and 9 the simulated frequency is 140 GHz, but to keep the moment method calculations tractable, the plate has been reduced in size to 1.4A0 by 2Ao. The curves shown are for a two-layer plate having d2 = 2dl = 0.5 mm with El = 5 + i4 and E2 = 2 + il, and for a single layer having the average permittivity earsvg = 3.5 + i2.5 (see Table 1). Since the accuracy of the physical optics approximation increases with the plate size, the agreement between the two-layer VIPO approximation and the moment method solution is remarkably good, and significantly better than if a single layer had been used. -20 -30 i -50 -7i0 -90 -m -w -38 -X L 38 m m w Sattl" Anle (Dew*") Fig. 7. The bistatic cross section of a 2A0 x 2A0 plate for H polarization with d2 = A0/50 and Eavs = 13 + i12 at normal incidence: moment method solution is denoted by the solid line. VIPO or SCPO is denoted by the dashed line. Fig. 8. The bistatic cross section of a 1.4Ao x 2A0 plate for E polarization with d2 = 2dl = 0.5 mm andf= 140 GHz at normal incidence: The solid line denotes moment method solution with El = 5 + i4, e2 - 2 + il; the short-dashed line denotes VIPO with El = 5 + i4, E2 = 2 + il; and the long-dashed line denotes VIPO with E2 = 61 = 3.5 + i2.5. 5. CONCLUSIONS A typical leaf has at least two dielectric layers whose cells have differing water content, and this produces a nonuniform dielectric profile which can now affect the scattering. At microwave frequencies where the leaf is no more than (about) A0/50 in thickness, the nonuniformity is not important, and as shown by Senior et al. [1987] the leaf can be modeled as a resistive sheet using an average value for the permittivity. If the physical optics approxi s8 tluem Age (oeene) Fig. 9. The bistatic cross-section area of a 1.4Ao x 2A0 plate for H polarization with d2 = 2dl = 0.5 mm andf = 140 GHz at normal incidence: The solid line denotes moment method solution with el = 5 + i4, e2 = 2 + il; the short-dashed line denotes VIPO with E1 = 5 + i4, 62 = 2 + il; and the long-dashed line denotes VIPO with e2 = 61 = 3.5 + i2.5.

SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF 15 mation is then applied, the resulting scattering is attributed to a surface current, and this method is equivalent to the SCPO approximation. At higher frequencies, however, the thickness and structure of a leaf are more significant. At 100 GHz and above, a leaf is a considerable fraction of a wavelength in thickness, and in spite of the reduced sensitivity to water content, the nonuniformity affects the scattering. For a two-layer model of a leaf the SCPO approximation has been compared with the volume integral (VIPO) approximation. When the leaf is thin, the two approximations are identical and in good agreement with data obtained from a moment method solution of the integral equation, but as the electrical thickness increases, the two approximations diverge in all directions except the specular and (for E polarization) backscattering ones. Although the VIPO approximation is more complicated, its accuracy is greater, and the agreement with the moment method data is better using a two-layer model than when a single layer of average permittivity is employed. For most practical purposes it would appear that VIPO in conjunction with an accurate dielectric profile of a leaf provides an adequate approximation to the scattering at millimeter wavelengths. As our knowledge of the profile increases, it may be desirable to use a multilayer model which could even simulate a continuous, nonuniform profile, and a convenient way of doing this is described in the appendix. We also note that at frequencies for which the leaf thickness is comparable to Am/2, where Am is the (average) wavelength in the leaf, the scattering is greatly reduced at some angle of incidence, and because the permittivity is complex, there is actually a range of angles for which this is true. Since the reduction is accompanied by an increase in the field transmitted through the leaf, this could provide a means for penetration through a vegetation canopy. APPENDIX Al. Combined sheets model When using the VIPO approximation, an efficient way to take into account the effect of any nonuniformity in the dielectric profile is to model the leaf as a stack of N combined current sheets. Each sheet simulates a very thin dielectric layer whose thickness is less than A/15 where A is the wavelength in the material. A combined sheet consists of coincident resistive and modified conductive sheets that support electric and magnetic currents, respectively, with the conductive sheet accounting for the electric currents flowing perpendicular to the dielectric layer. The mth layer sheets are characterized by a complex resistivity and conductivity Rm and R*,, respectively, where iZ0 Rm = i- kOAm (e m - 1) R*= - koAm(Em - 1) (Al) Here em and Am are the relative dielectric constant and thickness of the mth layer, and - = IN=l1 Am is the total thickness of the dielectric slab. The boundary conditions at the mth combined sheet are as follows [Senior and Volakis, 1987]: fix {fix [E + +E-I} = -2RmJm Jm = a X [H+ - H-] (A2) (A3) where Jm is the total electric current supported by the resistive sheet, and t x {a x [H+ + H}-] iYo a -i x- [E+ + E-]= -2RJ -k0 J* an = -a x [E+ - E-] (A4) (A5) where J* is the total magnetic current supported by the conductive sheet. The superscripts plus and minus refer to the upper (+) and lower (-) sides of the sheet, and fi is the unit vector outward normal to the upper side. A2. The induced currents of N planar sheets Consider a stack of N infinite planar combined sheets all parallel to the xy plane of a Cartesian coordinate system (x, y, z) as depicted in Figure Al. The top sheet is in the z = 0 plane, and the mth sheet is located at z = -dm, where d = O. The space between the mth and (m + I)th sheets is referred to as region m, and we note that region 0 (z > 0) and region N (z < -dn) are semi-infinite free space. A plane wave whose plane of incidence is parallel to the xz plane impinges on the stack of

16 SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF the following relations are obtained: Agion O E m -1 - 1 + (2 Y0 cos 80 Rm - l )e 2 iocoso(dm + - dm)rE 1 + 2Yo cos 80 Rm + e2ikocos~o(da+~ dm)rF (A7) dX-b X PAGrW1 Fisgon 2 -d3 Cm = + e2 ikocosO(d+ i - d,)rE Cm - 1 m The induced electric current in the mth sheet can be found from (A3) and expressed as (excluding the phase factor eiksinfOx) te dielecJm = 2Y0 cos 0oeikocossodm PAleon N Fig. Al. Layer of N combined sheets simulating infinil tric slab. sheets from above. From the symmetry of the problem, all the field vectors are independent of y (i.e., /alay = 0), as a result of which the field components in each region can be separated into E- and H-polarized waves which are the dual of each other. The analysis is similar to that employed in studying reflection and transmission by a layered halfspace [Kong, 1985]. In the case of E polarization the incident field is given by (1), and the field components in region m can be expressed as Emy = [Ce -ikocosOoz + Cr eikocosOoz]eikosin0ox Hmx = Yo cos Oo[Cie-ikocosOoz - C eikocosOoeeikosin0ox (A6) Hmz = Yo sin 80[C e-ikOcoso~z + C eikocos-0zleikosineox The coefficients C' and Cm are the amplitudes of the waves traveling in the -z and +z directions, respectively, in region m. In region 0, C' = 1 and C0 = rE (the total reflection coefficient), and in region N, Cr = 0 and C' = TE (the total transmission coefficient). Hence, using the boundary conditions (A2)-(A5), there are 2N unknowns and 2N equations that can be solved simultaneously. On substitution of (A6) into (A4) the left-hand side vanishes, showing J* = 0. As expected, the conductive sheet is not excited with this polarization since there is no current in the z direction in the dielectric slab, and in the absence of a magnetic source the tangential component of the electric field must be continuous as given by (A5). On inserting the expressions (A6) into (A5) and (A2) and defining the reflection coefficient in region m as rE = r e-2ikocosOodm+ Cm 1+ rE1 ~.1-m 1 + e2ikocoso~(dm + I - d,)rE m- e I Ik + r - I 1 + e2ikocoseo(d*+ + - dl)rE / =1I I+ (A8) The total reflection coefficient in region 0 (rE(6) = E) can be evaluated from the recursive recursive relation (A7) by noting that NF = 0 (the region N is semi-infinite). The total transmission coefficient can also be obtained from (A7) as follows: i' N I + rE CO mSCN + e Mikocos -1~d,+ l -d,)r~ T =C i -[ 1 + e2ikoc~se~(d'"+ i - dF)rFE rE 0 M= cI LI- +e (A9) Unlike the E-polarized case where the magnetic current is zero, an H-polarized wave excites a magnetic current in the y direction, and the tangential electric and magnetic fields are both discontinuous across the combined sheets. For H polarization the tangential field vectors in region m can be obtained by applying the duality relationships to (A6). In this case the amplitudes of the waves traveling in -z and +z directions are denoted by B'm and Br, respectively. By applying the boundary conditions (A2)-(A5) at the mth sheet and denoting the reflection coefficient in region m by Br rH _ e B-2ikocosOodm + B4 Bm after some algebraic manipulation we obtain rHm-I (QmPm - )-(l P,)(Qm - 1)FrHee2ikoCsOo(dm+l -dm) (1 +Pm)(l +Qm)+(l -QmPm)rFHe2ikocso(dm+l - di,) (A10)

SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF 17 Bi i(Qm-l)+(l+Qm)rmH l -(1 +Qm)+(-Qm)FrHe2ikocos0o(dm+ i - d,) B - I (All) where the parameters Qm and Pm are sin2 80 2R* Zo cos 00 2Rm sec 80 Zo (A12) Since FrN = O the recursive relation (A11) gives the total reflection coefficient in region 0, and using (A12) the total transmission coefficient is given by TH(O) Fig. A2. The geometry of the scattering of a plane wave from a finite N-layer combined sheet. N [ (Qm l)+(i+Qm)Frm m= -(1 +Qm)+(1-Qm)IrmHe2ikocosso(dm + -d ) (A13) The induced electric and magnetic currents can be expressed in terms of the reflection coefficients as JH = _ieikocoseodm- [(1 + rH I) + (1 + Fe2ikocoseo(dm+I - d,)) (Qm -1) + (1 + Qm)rml1 1 (1 + Qm) + (1 - Qm)rHeikOcOso(dm+l -dm) A3. Scattering by a rectangular stack Consider a portion of the N-layered stack of combined sheets in the form of a rectangle occupying the region -a/2 - x - a/2, -b/2 -< y - b/2 as depicted in Figure A2. In the far zone the approximation I - TmI ~ r + sin 6sx' + cos Osdm leads to () ikor iZo aI2 b/2/ N kor 41r J n Jm ( I )eikocose eikosin,'x' dx' dy' (A16) eikor iYo_2b/2 ( 9i r mje ikocosOdm\ ko' 4 J Ian -b/2 m - ) eikosin&sx' dx' dy' rm- I[1+ r>(Qt-l)+(l+Qt)rf-l t- Q)fe a:O(d, -d-) I (Il+Ql) +(I1-lr~ kco d (A14) JH* - Zo cos O eikocosOod - I [(1 rm- 1) + (1 - He2ikocoseo(d + - dm)) (Qm - 1) + (1 + Qm)rfm (1 + Qm) + (1 - Qm)rmHneakocoso(dhl -dm) m [1 (Ql - (Q-)+(l+Qi)r_ I tl _ ( 1 l) + (1-Q,)FfreikocosOo(di + - di) 1=1 -(iQ,)+(iQl~rl] Using the physical optics approximation, the currents obtained for the infinite sheets are substituted into (A16) to find the scattered fields. For E and H polarizations the far-field amplitudes are i N sin X SE(S, 0o)= 9i- k2abZo ( JmEeikcssdsn X (A17) (AI5) SH(O, 80) = 9 i k2ab 4ir where JH(x) = JHeikosin0ox JH*(x) = JH*eikosin0ox m m IN [ COS NrH*H1.ikocosO,d. sin X (COS sJm OJm,)ec (A18) m=l L. c-sH Y 1H*a*Y i1sin]

18 SARABANDI ET AL.: MILLIMETER WAVE SCATTERING MODEL FOR A LEAF where, as before, X = (koa/2) (sin Os + sin 80). In the backscattering (6s = 00) and specular (6s = -80) directions the summation term in (A17) reduces to a telescopic series resulting in N JEeikocoseod = 2Yo cos 0o CO = 2Yo cos 0oFE(O) m=l (A19) and the backscattering cross section is then (ab)2 ~aE(80, 00)=4r k2 Cos2 001rE(80)12 sin2 (koa sin 00) 0 V (A20) (koa sin o)2 A0) Also, for H polarization, N (COS 0oJH + yoJ*)eikocosodm m=N =- cos 00 Z (Bm 1-Br) m=l = -2 cos OoBO= -2 cos 0o FH(0o) (A21) which leads to (ab) 2 "H(80, 00)=47r a2 Cos2 olFrH( 0o)l sin2 (ka sin 0o) 0 - 00)2 (A22) (ka sin 60)2 These results are identical with the ones obtained by the SCPO technique. The extinction cross sec tion can also be obtained from the far-field amplitude using A2 ext= - Im [S(80, 00 + 7r)] 7r (A23) from which we obtain aoext = 2ab cos 80 Re [1 - TE(O)] CHrxt = 2 ab cos 80 Re [1 - TH(0O)] (A24) Acknowledgment. This work was supported by the U.S. Army Research Office under contract DAAG 29-85-k-0220. REFERENCES Beckmann, P., The Depolarization of Electromagnetic Waves, Golem Press, Boulder, Colo., 1965. Curtis, H., and N. S. Barnes, Invitation to Biology, Worth Publishers, New York, 1985. Kong, J. A., Electromagnetic Wave Theory, John Wiley, New York, 1985. Le Vine, D. M., A. Snyder, R. H. Lang, and H. G. Garter, Scattering from thin dielectric disks, IEEE Trans. Antennas Propag., 33, 1410-1413, 1985. Sarabandi, K., Electromagnetic scattering from vegetation canopies, Ph.D. dissertation, Univ. of Mich., Ann Arbor, 1989. Sarabandi, K., T. B. A. Senior, and F. T. Ulaby, Effect of curvature on the backscattering from a leaf, J. Electromagn. Waves Appl., 2, 653-670, 1988. Senior, T. B. A., and J. L. Volakis, Sheet simulation of a thin dielectric layer, Radio Sci., 22, 1261-1272, 1987. Senior, T. B. A., K. Sarabandi, and F. T. Ulaby, Measuring and modeling the backscattering cross section of a leaf, Radio Sci., 22, 1109-1116, 1987. Willis, T. M., H. Weil, and D. M. Le Vine, Applicability of physical optics thin plate scattering formulas for remote sensing, IEEE Trans. Geosci. Remote Sens., 26, 153-160, 1988. K. Sarabandi, T. B. A. Senior, and F. T. Ulaby, Radiation Laboratory, Electrical Engineering and Computer Science Department, University of Michigan, Ann Arbor, MI 48109.

I ADVISORY GROUP FOR AEROSPACE RESEARCH & DEVELOPMENT 7 RUE ANCELLE 92200 NEUILLY SUR SEINE FRANCE Paper Reprinted from AGARD Conference Proceedings 501 Target and Clutter Scattering and their Effects on Military Radar Performance (Diffraction par les Cibles et le Fouillis et ses Effets sur les Performances des Radars Militaires) * i - NORTH ATLANTIC TREATY ORGANIZATION I.I

28-1 MILLIMETER WAVE POLARIMETRIC SCATTEROMETER SYSTEMS: MEASUREMENT AND CALIBRATION TECHNIQUES Y. Kuga, K. Sarabandi, A. Nasbashibi, F. T. Ulaby and R. Austin University of Michigan The Radiation Laboratory Department of Electrical Engineering and Computer Science 3228 EECS Building Ann Arbor, Michigan 48109-2122 SUMMARY The target and system phase-stability during the time to measure the scattering matrix is a major problem for millimeter wave polarimetric radars. This is particularly true for network analyzer-based systems. To circumvent this phase-stability problem, we have developed new fully polarimetric radars at 35 and 94 GHz. The system is based on a relatively inexpensive network analyzer and is capable of operating in either the coherent or the incoherent polarimetric measurement mode. In the coherent mode, the scattering matrix can be measured within 2 ms. In the incoherent mode, the average Mueller matrix is measured directly by transmitting four different polarizations and measuring the Stokes vector of the backscattered signal. To compare the performance of the true measurement modes, the average Mueller matrix and the statistics of the phase difference of the two co-polarized signals were measured for a rhododendron tree and for a metallic tree. The average Mueller matrices obtained from the coherent and incoherent polarimetric measurement modes were similar. The target motion during the data acquisition period did not change the average Mueller matrix in the incoherent measurement mode. The probability density function of the phase difference of the two co-polarized signals computed from the average Mueller matrix is essentially the same as the one measured with the coherent polarimetric measurement mode. 1 INTRODUCTION Increasing interest has been expressed in recent years for understanding the statistical properties of data obtained with fully polarimetric radars for remote sensing applications [Ulaby and Elachi, 1990]. At centimeter wavelepgths, polarimetric data has been found to be useful for land-use classification [Van Zyl et. al., 1987] and for measuring the biophysical properties of forest canopies [McDonald et. al., 1990]. For the MMW region, however, it is still not clear what type of information can be extracted from polarimetric radar, over and above the magnitude information provided by conventional radar systems. Unlike the microwave region, the complexity and the cost of building a fully polarimetric radar at millimeter-wave frequencies is still very expensive, and progress has been rather slow, which is due, in part, to the limited availability of experimental data. At microwave frequencies the traditional approach used for measuring the polarimetric radar response of a given target is based on the direct measurement of the target's scattering matrix, S. For distributed targets, such as terrain surfaces, multiple measurements of S are made, corresponding to statistically independent samples, each measurement is used to compute its corresponding Mueller matrix L, and then an ensemble average is performed to obtain an estimate of the average Mueller matrix, < C >. Whereas the scattering matrix measurement technique is appropriate at microwave frequencies, it is difficult to implement at millimeter wavelengths because it requires that both the system and target phases remain stable during the time it takes to measure S. This is particularly true for network analyzer-based polarimetric radars [Ulaby et. al., 1990]. To circumvent this phase-stability problem, we have developed new fully polarimetric radars at 35 and 94 GHs. The system is based on relatively inexpensive network analyzer and is capable of operating in either the coherent or the incoherent polarimetric measurement mode. In the coherent mode, the scattering matrix can be measured withinr 2 ms. In the incoherent mode the average Mueller matrix is measured directly by transmitting four different polarizations and recording the horizontally polarized and vertically polarized' components of the backscattered field. This paper includes a detailed analysis of the two measurement modes, and provides comparisons of data measured using the two modes for a rhododendron tree and an artificially made metallic tree. 2 NWA BASED POLARIMETRIC RADARS The fully polarimetric radar configuration based on the vector network analyzer (NWA) is easy to construct and is widely used for remote sensing investigations [Ulaby et. al., 1990]. These systems usually are operated in the swept frequency mode over a given bandwidth. The minimum sweep time, which depends on the number of frequency points and the type of NWA, is typically between 100 to 400 ms. The decorrelation time of the MMW wave scattered from trees, on the other hand, can be shorter than 10 ms [Narayanan et. al., 1988]. Hence, when using the fully coherent measurement configuration, it is necessary that all four components of the scattering matrix be measured within a few milliseconds in order to obtain accurate data. If the V- and Hpolarized signals are transmitted sequentially in the swept frequency mode, it will take at least 0.5 to 1 second to get a complete scattering matrix, including the data transfer time between the NWA and the computer. Obviously, the NWA-based MMW radar used in the swept frequency mode is not suited for coherent polarimetric measurements. There are two ways to overcome the shortcoming of the traditional swept-frequency NWA based polarimetric radar. The first approach is the incoherent polarimetric measurement technique. With this technique the swept frequency mode can still be used for the NWA operation but the radar transmitter must be modified to transmit four independent polarizations. The data processing and calibration are substantially more complicated than those associated with the coherent polarimetric technique. The second approach is the coherent polarimetric measurement technique using Coupled/Chop mode and point by point external triggering of the NWA. We have developed both coherent and incoherent polarimetric radars based on these techniques at 35 and 94 GHz. The radar front end and data acquisition system are the same for both systems. The only difference is the operating mode of the NWA and the data processing. It is, therefore, possible to obtain polarimetric data of the same targets coherently and incoherently. The block diagram of the MMW/ radar system and the 35 GHz front end are shown in Figs. 1 and 2. The block diagram of the 94 GHz is essentially the same as that of the 35 GHz system. In the following section the details of the coherent and incoherent systems will be discussed. 2.1 Coherent Polarimetric Radar The coherent polarimetric radar has many advantages over the incoherent polarimetric radar. For example, with the coherent polarimetric radar the statistical data including the phase difference between the two copolarized channels, can be easily obtained. Another advantage is the significantly simpler signal processing and calibration processes compared to those of the incoherent polarimetric radar. As discussed in the previous section, the NWA-based radar operated in the swept frequency mode is not suited for coherent polarimetric measurements. In this section, we will describe a new technique which utilizes a relatively inexpensive NWA (HP8753C) that allows the acquisition of coherent polarimetric data at a much faster rate. With this system it is possible to measure the scattering matrix within 2 ms at 35 and 94 GHz. The Hewlett-Packard network analyzer, HP8753C, has two independent receiving channels which can be used in the Coupled/Chop mode. It also has a point by point external triggering capability in the swept frequency mode. These functions are ideally suited for the coherent polarimetric radar. For example, the simultaneous acquisition of V and H channels can be done by operating A and B inputs in the Coupled/Chop mode. The point by point external triggering can be used for transmitting V and H sequentially and synchronizing a polarization control circuit to create different polarizations. At present, HP8753C does not support the external point by point triggering in the CW mode but a near CW mode can be created

28-2 in the swept frequency mode by choosing the output frequency bandwidth to be 1Hz. The minimum time to get a complete scattering matrix is approximately 2 ms in the present system. The polarization of the transmitted MMW signal is controlled by a Faraday rotator whose switching time is less than 5 us. Using the maximum number of points provided by the HP8753C, it is possible to obtain 800 scattering matrices within 3.2 a without transferring data into a computer. The separation of signal from unwanted noise, such as antenna coupling, is accomplished by the hardware gating circuit in the IF path as shown in Fig. 2. The transmitted pulse length is 20 ns and the pulse-repetition-rate is 5 MHz. Although it is not necessary to scan the RF frequency band in the coherent polarimetric mode, additional independent samples can be realized by averaging the backscattering coefficient over the RF bandwidth [Ulaby et. al., 1988]. A bandwidth of 1 GHz at 35 GHz, for example, offers 5 to 10 independent samples per spatial observation for the tree measurements. The calibration of the coherent system is straightforward. Because the system has more than 23 dB of isolation between the V and H channels, a simple calibration technique that requires a sphere and a depolarizing target is used [Sarabandi et. al., 1990]. 2.2 Incoherent Polarimetric Radar In the incoherent polarimetric radar technique, the Mueller matrix of the target is measured directly by transmitting four independent polarizations and receiving the Stokes vector of the scattered signal. Because the correlation between the V- and H- polarized signals is inherently included in the received Stokes vector, the measurement time between the different incident polarizations can be much longer than the decorrelation time of the target. The incoherent polarimetric technique also permits the use of MMW sources that do not have good phase-stability in the transmitter section [Mead, 1990]. A desired polarization can be created by placing two quarter-wave plates in front of the transmitting antenna and by adjusting the orientation angle of each wave plate relative to the incident polarization. The received Stokes vector for a given incident polarization is usually obtained by two different approaches, incoherent and coherent-on-receive techniques. The incoherent receive technique, which often is employed in optics measures the intensity of six different receive polarizations, but the phase measurement is not required. The Stokes vector is obtained by taking appropriate ratios of the receive intensities, as shown in Appendix A. The receiver of the coherent-on-receive technique is similar to that of the coherent polarimetric radar. The coherent-on-receive method requires the measurement of the magnitudes of the Vand H- polarized receive signals and the phase difference between them, but it does not have to measure the phase angle relative to the transmitted signal, as is the case with the coherent polarimetric radar. The Stokes vector can be computed from the magnitudes of the V and H components of the received signal and the phase difference between them as shown in Appendix A. Because it is relatively easy to measure the phase difference between the V and H channels, our system is based on the coherent-on-receive technique. Calibration of incoherent polarimetric radar systems involves two steps [Mead, 1990]. In the first step, the receiver distortion matrix is obtained by placing a wire grid polarizer in front of the receiving antenna at three different positions. In the second step, the exact polarization properties of the transmitter are determined by measuring the backscatter from a point target with known scattering matrix using the calibrated receiver. 3 EXPERIMENTAL DATA To demonstrate that the coherent and incoherent polarimetric measurement techniques do indeed provide identical information for distributed targets, experiments were conducted using a rhododendron tree and a metallic structure resembling a short tree. Photographs of these targets are shown in Fig. 3. The metallic structure is used for creating a target return with strong correlation between the S,, and SA,, components. To create many independent samples and also to show that the incoherent polarimetric technique can provide accurate results even if the data acquisition time is much longer than the target decorrelation time, the trees were rotated at slow (0.67 rpm) and fast (1.33 rpm) speeds during the data-collection process. Table 1 shows the average Mueller matrices of the rhododendron and metallic trees obtained by the coherent and incoherent polarimetric radar techniques. The Mueller matrices are normalized with respect to the L11 component to show the relative magnitude of the matrix elements. The average Mueller matrix of the coherent polarimetric radar technique was computed from the 8000 scattering matrices obtained over the 34-35 GHz band. Because of the slow data-acquisition speed in the incoherent polarimetric radar, the average Mueller matrix is obtained from only 500 samples, including those due to frequency averaging over the 1-GHz RF bandwidth. The sum of the Mueller matrix elements L33 and L44, which is a function of the correlation between S,, and Shh, is higher for the metallic tree than for the rhododendron tree. Although the polarimetric signature computed from the average Mueller matrix is useful for showing the characteristics of the target, it is not easy to directly relate the target characteristics to the values of the Mueller matrix elements. Figure 4 shows the probability density function of the phase difference between the two co-polarized channels (0b = phase of S,, - phase of Shi) obtained with the coherent polarimetric radar. As expected, p (~) of the metallic tree is much narrower than that of the rhododendron tree, showing strong correlation between S,, and Shh. Unlike the coherent polarimetric radar, the information obtained with the incoherent polarimetric radar is limited to the average Mueller matrix and it is not possible to measure the probability density function p (,) directly. Due to a recent theoretical derivation, however, the phase statistics of X, can be estimated from the average Mueller matrix [Sarabandi, 1991]. The probability density function p (0,) is given by P 13All>"-A 13-4A1?+4+ D [+ t I where B where Lll All = 2' 2 ' L33 + L44 4 L22 2 ' Ls4- L43 4 D = A13 cosc + IAI4 sin 4c, B = [LLA33 - D2] t The function p (0) is completely specified in terms of the elements of the average Mueller matrix m,. The average Mueller matrix given in Table 1 and the probability density function of the phase difference shown in Fig. 4 are obtained from the same target by two different polarimetric measurement techniques. If the probability density function given by Eq. 1 is correct, p (,) estimated from the average Mueller matrix must be similar to the one shown in Fig. 4. Figure 5 shows the probability density function computed from the average Mueller matrix obtained by the incoherent polarimetric radar. The agreement between Figs. 4 and 5 is excellent for both trees. 4 CONCLUSION The work described in this paper has demonstrated that the average Mueller matrices obtained using the coherent and incoherent polarimetric measurement techniques are essentially identical. The advantage of the coherent polarimetric radar over the incoherent polarimetric radar is its ability to measure the statistical distributions of the magnitudes and relative phases of the scattering matrix elements. The incoherent polarimetric radar, however, is particularly useful if the target decorrelation time is much faster than the data acquisition time. REFERENCES (11 Ulaby, F.T., and C. Elachi, Radar Polarimetry for Geoscience Applications, Artech House Inc., Massachusetts, 1990. [21 Van Zyl, J.J., H.A. Zebker, and C. Elachi, 'Imaging Radar Polarization signatures: Theory and Observation," Radio Science, 1987, Vol. 22, pp. 529-543. [3) Sarabandi, K., F.T. Ulaby, and M.A. Tassoudji, "Calibration of Polarimetric Radar Systems with Good Polarization Isolation," IEEE Tins. on Geosci. Remote Sensing, 1990, GE-28.

28-3 [41 Ulaby, F.T., T.F. Haddock, and R.T.. Austin, "Fluctuation Statistics of Millimeter-Wave Scattering from Distributed Target," IEEE TLamU. on Geosci. and Remote Sensing, 1988, Vol. 26, pp. 268-281. [5] McDonald, K.D., M.C. Dobson, and F.T. Ulaby, "Using MIMICS to Model L- Band Multiangle and Multitemporal Backscatter from a Walnut Orchard," IEEE Trans. on Geosci. and Remote Sensing, 1990, Vol. 28, pp. 477-491. [6] Mead, J.B., "Polarimetric Measurements of Foliage and Terrain at 225 GHz," (PH.D. Thesis, Univ. of Mass., 1989). [7] Ulaby, F.T., M.W. Whitt, and K. Sarabandi, "AVNA-based Polarimetric Scatterometers," IEEE Antennas Propagation Magazine, 1990. [8] Sarabandi. K., "Derivation of Phase Statistics of Distributed Targets from the Averaged Mueller Matrix," Technical Report 026511-T, 1991, Radiation Laboratory, University of Michigan, Ann Arbor, MI. [9] Narayanan, R.M., C.C. Borel, R.E. McIntosh, "Radar Backscattering Characteristics of Trees at 215 GHz," IEEE Trans. on Geosci. and Remote Sensing, 1988, Vol. 26, pp. 217-218. APPENDIX A COHERENT-ON-RECEIVE TECHNIQUE Complete polarimetric characterization of the scattering properties of a distributed target can be obtained by measuring either the scattering matrix S or the Mueller matrix 4m,. Measurement of the scattering matrix requires accurate phase measurements. Also 4 elements of S must be obtained within the decorrelation time of the target which is in the order of milliseconds at MMW frequencies. The scattered electric field E', in terms of the scattering matrix S and the incident electric field Et, is given by where X, is called the Mueller matrix. The totally incoherent method does not require phase measurements. With this method, the 4 elements of Stokes vector are obtained by receiving 6 polarizations (V, H, 45,135, LHC, RHC). For example, the third element of Stokes vector U is given as a ratio of intensities at 45 linear to 135 linear. For a given incident polarization, therefore, we can obtain a column of the Mueller matrix. To get the complete Mueller matrix, we need to repeat this process for 4 independent incident polarizations. Altogether, at least 24 magnitude only measurements are required to obtain the complete Mueller matrix. Although the phase measurement is not required with the incoherent method, it is necessary to receive all 6 polarizations. The elements of the Stokes vector, in terms of 6 polarizations and a set of 4 independent incident polarizations, given by: I= (A.6) WI Wh + Wv Wh Is = "'a+W, (A.7) Wh + W. U = W 5 Wl- (A.8) w45 + wIN WLHC - WRtHC WLHC + WRHC I. = [0 Ih' = 0[] Is5 = [1. ILHC= 1/21 (A.10) where W is the received intensity of polarization. If a receiver is able to measure the phase between the V and H channels, it is possible to do the incoherent method without measuring 6 polarizations. This method is called the coherent-on-receive (COR) technique. The elements of the Stokes vector can be expressed as 1I =1 E. 12 U = 21 E, II E, I cos6 V =2 E, I1 E I sin6 (A.11) (A.12) (A.13) (A.14) E = SE' r Eh =Eh S[$hS, Sh (A.1) (A.2) where 6 is the phase difference between V and H channels. To obtain S, we need to send [E,, 0]' and [0, EA]', and measure E, and Eh simultaneously. The polarized wave can also be expressed in terms of the Stoke's vector Fm which is defined as 'I ' I E, 12 Fr = = - RIEt 1] mlfl I 2Re[E.EA)] V t2Im [E, El] then (A.1) in terms of Stokes vector becomes (A.3) F-1 CmF r = (A.4) ISb, I IISh21 2R e(SS, ) 2Re(ShS) * *I) L 2Im(S.Sh,,) 2Im(S,AS*,) Re(SS,,S,) -Im(S*,s,,S.) Rm(Sn, Sa + ) -Im(S. - SaA) 'a' KRc(S..S,, + Sqa) -Im(S$,, — S"AS,) Im(SS;,h + S.,AS,,) Re(S,.Sa -SiS.) (A.5)

28-4 Radm Figure 3 Photographs of the Metallic Tree and Rhododendron Tree. Figure 1 Block Diagram of the MMW Polarimetric Radar. 35 GHz Radar (Fully Polarlmetric) Transmitter IF bandwidth up to 2 GHz Power +23 dBm Antenna 6" Lens (beamwidth 4.2 degrees) Polarization Any polarization Receiver Dual Channel V and H Mixers Fundamental mixing Antenna 6' Lens (beamwidth 4.2 degrees) Polarimetrle data Incoherent (coherent-on-receive) Mueller matrix Coherent Scattering amplitude matrix Rhododendron Tree (Target in motion) Incohaent Polarimetric Measurment Mode (500 samples) Fast Motion 1 0..16 -0.006 0.007 0.186 0.828 -0.017 0.017 0.04 0.059 0.735 0.056 -0.059 -0.023 -0.019 0.472 Slow Motion 1 0.231 -0.052 0.126 1.035 -0.021 -0.038 0.082 0.697 -0.013 -0.015 -0.064 0.012 ] 0.049 0.003 0.619 POL Coherent Polarimetric Measurment Mode (8000 samples) 1 0.159 -0.002 -0.006 0.179 0.823 -0.003 -0.018 0.0 -0.01 0.683 -0.023 -0.033 -0.001 0.003 0.596 Metallic Tree (Target in motion) Incoherent Polarimetric Measurement Mode (500 samples) Fast Motion 1 0.089 0.02 0.03 0.094 0.74 0.02 0.001 0.026 0.011 0.888 0.105 -0.028 -0.007 -0.126 0.619 Coherent Polarimetric Measurement Mode (8000 samples) 1 0.06 0.0 0.0 0.06 1.16 0.009 0.004 0.011 0.002 0.973 -0.053 0.008 0.002 0.049 0.89 Table 1. Average Mueller matrices of rhododendron and metallic trees measured with coherent and incoherent polrimtric measurement ndes. Figure 2 Block Diagram of the 35 GHs Radar Frontend.

28-5 1.0 * 1.0 - 0.8 0.8 M i., a 0.6 - 11$~~~~~~~~~~I 0.6 -Rhododendron tree.L Rhododendron tree Rhoddendro tree (Slow motion) 0.e Rhododendron tree (Past motion). A Ir.150 -100 -50 0 50 100 150 -150 -100 -50 0 50 100 150 Degree (0c) Degree (#c) Figure 4 Probability Density Function Measured by the Coherent Figure 5 Probability Density Function Computed from the Average Polarimetric System. Mueller Matrix. DISCUSSION What is the realized cross polarization isolation? Author's Reply I am not sure but I believe that it is listed in the preprint article. E. Schweicher, BE You mentioned (in your second talk) a solid state design. Does that mean a solid-state transmitter and in that case what kind of transistors are used at 36 GHz and 94 GHz? Author's Reply The presenter did not know. G. Brown, US The histograms of the phase differences between the co- and cross polarization returns for data (from real and metallic trees) were relatively narrow in spread. What range bin was used to generate these data and how do you explain the rather small phase difference variance? Author's Reply I do not know the range bin used in the phase distribution plots since I did not work on those data. However, I do know that a very large number of samples were used in the preparation of the phase statistics, so I imagine that the bin was rather small.

1r7 POLARIMETRIC OBSERVATIONS OF TREES AT 35 AND 94 GHZ Adib Nashashibi*, Yasuo Kuga, and Fawwaz T. Ulaby Radiation Laboratory Department of Electrical Engineering and Computer Science University of Michigan, Ann Arbor, MI 48109 I. SUMMARY Two different kinds of trees were measured with a newly developed polarimetric radars at 35 and 94 GHz. The system is based on the COR technique and is capable of measuring the Mueller matrix directly. The degree of polarization, phase statistics, and polarization signatures are presented in this paper. Our results show that the degree of polarization is sensitive to the radar frequencies and tree types. II. SYSTEM DESCRIPTION The block diagram of the University of Michigan MMW polarimetric radar is shown in Fig. 1. The 35 and 94 GHz radars are fully polarimetric and are capable of measuring the Mueller matrix using the coherent-on-receive (COR ) technique [2,3]. With the COR technique, 6 different polarizations (V, H, 45, 135, LHC, and RHC) are transmitted sequentially, and the phase and magnitude of V and H channels are received simultaneously. Unlike the coherent polarimetric radar which measures the scattering amplitude matrix S, the COR radar measures the Mueller matrix Lm. The original MMW system was based on the HP8510A and it did not have the capability of obtaining data from two different channels simultaneously. To overcome this problem, a delay line technique was developed. With this technique, one channel, usually H, is delayed by a long delay line (200 nsec delay) and combined with the non-delayed channel. The received data in the time domain shows V and H channels separated by the electrical length of the delay line. Using the time gating function of the network analyzer, V and H channels can be processed separately. The calibration of the COR system is performed in two steps. First, the receiver is calibrated using an odd bounce reflector and a polarizing grid placed in front of the receiving antenna. The receiver distortion matrix is obtained by positioning the polarizer at three different angles. Second, the actual transmitted polarizations are obtained using the calibrated receiver and the odd bounce reflector [3]. III MUELLER MATRIX AND DEGREE OF POLARIZATION The polarization of the EM wave can be expressed by the modified Stokes vector, F,, which is defined in terms of the vertical, Ev, and horizontal, Eh, electric fields [1]. Fm=[[l, I2, U, yVT, Il=Ev 212.= I 1Eb 12, U=2Re(EvEh), and V=2Im(EvEh*) (1) The transmitted and received modified Stokes vectors are related to each other by F,r = (1/r2) L, Fmt (2)

I S 12 IjS I 2 Re(S S S) -Im(S S ) L,, = | S.2 IS,,1 Re(S;S,.) -[ S,). | Lxg SkV S kv Id ' ' (3 2Re (S S. ) Re (S (S V-S + SS) m(S~ S' - S So 21m(S$S$) 21m(S,,,S,,) Im(S,,S, + S,,S) Re(SvS; - SSv,,)! where r is the distance between the target and the receiving antenna, Lm is the modified Mueller matrix of the target, and SjL (l,k=v or h) are the four elements of the scattering matrix S. The COR radar obtains Lm directly by sending 4 to 6 independent polarizations. Using the elements of Lm, the phase difference of the co- and -crosspolarized components of S can be written as vv -hh = tan-1[(-43-L34)/(L33+L4)] (4) wvv -hv = tan-1(L41/L31) (5) 4vh -hv = tan'- [(L43+L34)/(L33-L44)], (6) The degree of polarization, m, which is directly obtained from the received Fm, is given by m= \q((<IIP> - <12p>)2+<Up>2+<Vp>2) / (<IP>+<12p>) (7) whereo is an ensemble average, and "P" denotes the transmitted polarization. IV. RESULTS Measurements were conducted using leafy evergreen and coniferous trees. The first was a rhododendron (rhodo) with a planophil leaf orientation distribution and an average leaf area of 40 cm2. The second was a spruce blue shiner (spruce) with a uniform leaf orientation distribution and an average needle length of 1.5 cm. To obtain the statistical data, more than 100 independent samples were measured for each incidence angle. Table 1 shows the degree of polarization for 6 different incident polarizations at 35 and 94 GHz. The degree of polarization is quite sensitive to frequencies and tree types. The MMW return from spruce tree, for example, shows a significant change in m between 35 and 94 GHz The dependence on the incidence angle is also clear for the planophil type tee. Figure 2 shows the phase distribution of the received signals. In all cases, Ov -hv had a uniform distribution showing the decorrelation of co- and cross-pol. components. The most interesting case is the phase difference between the two copolarized returns shown in Fig. 2. The distribution of v -hh seems to depend on the tree type and on the incidence angle. *w -hh of rhododendron has a much narrower spread than that of spruce. The phase *h -hv of a truly monostatic system should be a delta function centered at zero degrees. Our data, however, shows a wide distribution. We think this is caused by the bistatic configuration of the radar system.

Figure 3 is the co-polarized polarization signatures of rhododendron and spruce. The pedestal which is related to depolarization is much larger for spruce tree showing a strong depolarization. IV. REFERENCES 1. F.T. Ulaby and C. Elachi, Radar Polarimetry for Geoscience Applications, Artech, 1990. 2. F.T. Ulaby, M. Whitt and K. Sarabandi,"AVNA-Based Polarimetric Scatterometers",IEEE AP Magazine, Vol 32, 1990. 3. J. B. Mead, Polarimetric Measurements of Foliage and Terrain at 225 GHz", Ph. D. thesis, University of Massachusetts, Amherst, MA, 1990. TABLE 1. Degree of polarization m of rhodo and spruce trees at 35 and 94 GHz. PL..35 GHz 94 GHz RHODO SPRUCE RHODO SPRUCE 300 500 900 30~ 900 50~ 900 900 V 0.86 0.89 0.80 0.47 0.60 0.92 0.82 0.64 45 0.80 0.85 0.66 0.48 0.43 0.91 0.87 0.71 H 0.76 0.83 0.69 0.37 0.57 0.87 0.87 0.71 135 0.83 0.87 0.69 0.49 0.49 0.82 0.77 0.58 LHC 0.69 0.77 0.61 0.1 0.11 0.81 0.73 0.39 RHC 0.61 0.79 0.51 0.06 0.21 0.77 0.73 0.39 J4 phase plates tmum - —, m Delay line rigure 1: Block digageu of the CO. systm

20 15 10 - 5 - 0 - 20 - 15 - 10 5. 0 - a 1 I o * rhodo., vv-.hh, 94 GHZ, 50~1 0 rhodo., )vv-hh, 35 GHz, 50" 5 0 2 1n4 iI L! nn a O O O a O O O0 q N O 4 Co q N Phase Difference 2: Distributions of rvv-hh O o O N q O O O O 0 40 O N a a _ O O O qD a m VW v. 0vv.hh (Degrees) at 35 and 94 GHz Figure: Spruce Tre at 3 GaHz, W Rhodo. Tree at 35 Hz, 900 t Figure 3: Copolarized signatures at 35 GHz and 900 incidence angle.

U., UI' tKt MEASUREMENTS BY MILLIMETER-WAVE RADARS Yasno Kuga, Adib Nuashibi, and Fawwaz T. Ulby D ticnt of ecrical Engineering and Computer Science University of Michigan, Ann Arbor, M 48109 SUMMARY An extensive radar clutter database was generated by the University of Michigan's millimeter-wave mobile polarimetric radar system during the past few years [11. The data base includes millimeter-wave observations of snow, trees, vegetation, and soil and road surfaces at 35, 94, 140 and 215 0Hz. The radar measurements were often augmented with close-up observations of the targets including such measurement as water contents and surface roughness, when appropriate. For each data set, a summary of these observations and photographs of the target scene are provided. The millimeter-wave system consists of truck-mounted radars capable of making observations from a 20-m high platform at incidence angles between 0 and 70 degrces The 35 and 94 oHz radars are fully polarimctric and capable of obtaining the Mueller matrix in sit. MEASURIMENT SYSTEMS Radar System The block diagram of the MMW radar system at the University of Michigan is shown in Figure 1. The system consists of MMW radar front-ends at 35, 94, 140, and 215 GHz, a receiver, and a system control unit. The receiver is based on the HP8753 vector network analyzer (NWA). The IF signal (2 to 3 GHz) from the NWA is upconverted to the MMW frequency in each radar front-end and the reccived MMW signal is down-converted to the IF signal and fed into the NWA. Both software and hardware gatiags are utilized to eliminate the unwanted noise from the returned signal. Because the system has one GHz bandwidth, the software gating using the time domain capability of NWA gives an approximately 15 cm spatial resolution. Unfortunately, it takes 400 millisecond to sweep 401 frequency points on HP8753. If a target is moving, the sweep time can be much longer than the target decorrelation time at the MMW frequency. The software gating, therefore, i mainly used for calibration and indoor measurements in which the target motion is negligibly small. For outdoor tree measurements, a hardware gating circuit, which creates 20 nsec pulses with a pulse-repetition-rate of 5 MHz, is used. Because the pulse-repetition-rate is much greater than the NWA IF bandwidth (3KHz), the system is essentially operating in the CW mode. The original MMW system was based on the HP8510A and it did not have the capability of obtaining data from two different channels simultaneously. To overcome this problem, a delay line technique was developed. With this technique, one channel, usually H, is delayed by a long delay line (200 nsec delay) and combined with a non-delayed channel. The received data in the time domain shows V and H channels separated by the electrical length of the delay line. Using the time gating function on the NWA, V and H channels can be processed separately. The data acquisition and the system control arc performed by a HP computer with a control program written in HP BASIC. An IBM AT-compatible computer controls the polarization control unit which can be controlled by HP computer in the automatic mode through a RS232 interface. 35 and 94 Radar Front-ends The 35 and 94 OHz radars are fully polarimetric and capable of measuring the Mueller matrix using the coherent-on-receive (COR) technique [2,3]. With the COR technique, 6 diffeent polarizations (V, H. 45, 135, LHC, and RHC) are ransmitted sequentially, and the phase and maanitude of V and H channels ar received simultaneously. Unlike the coherent polarimctic radar which measures the scattering amplitude matrix S the COR radar measures the Mueller marix Lm. The description of COR technique is given in Appendix I. The block diagram of the 94 GHz radar is shown in Figure 2. The ansmitted wave polarization is controlled by two in tl rttable quaer-wa plte A precise control of wave plate motion is provided by stepping motors which are controlled by an IBM AT-compatible computer. To stabilize the MMW LO, a stable X-band signal is multiplied and injected into the Gunn oscillator. In our original design, the upconvered MMW signal was transmitted without amplificatio resulting m poor S/N ratio. To improve the S/N ratio and to obtain stable phase da, MMW amplifiers were added to boost the output power from +4 dBm to +23 dBm at 35 GHz and from +0 dBm to +13 dBm at 94 GHz........_ -. - r.. ---- -0 -I I w0 ILO.. _... Fig. 2. Block diagram of 94 GHz radar. 35 GHz radar has a similar design. Rdm ilo Fig. 1. Block diagram of MMW radar systecm

10 140 and 215 $Hz Radar Front-ends The present 140 and 215 GHz radars are limited to co- and cross-pol. measurements only. The block diagram is shown in Figure 3. The polarization of transmitted and received signals are controlled by wave plates and the receiver is limited to measuring either V or H channel at a time. Because Gunn oscillators are not available at these frequencies, the MMW signal is obtained by multiplying a lower frequency signal by 2 at 215 GHz and by 3 at 140 GHz. The transmitted power is -4 dBm at 140 GHz and -10 dBm at 215 GHz. I I I I I I I I 35 GHz k II t, v e.f.1. - - wv,/......... M l * %',~............... Itl IIulums --- I I I II I I Is ~ J ~ II ~ g ~ -10. II a ssa MIXER.20. in SWEPT RF SOURCE - I 0. L I I I I I 0 I I _ I I I, I, L r' `~~r- ~~ I I~~~~~~ I~~~~ I 10. 12. 14. 16. 1. 20. 22. 24. I - vVM bim VNHP-I m- 94 GHz *-,.- mm Idmtdm ].20. a I!... I..t.... & 10. lQ 1 l 16 IL 20..L U. 2,A. f *I bi Fig. 3. Block diagram of 140 and 215 GHz radars. RXPERIZMEN AL DATA In the following section we will describe two data sets obtained with our MMW system The first is a snow diurnal data set obtained in March 1990 and the second data set shows the polarimetric observations of two kinds of trees conducted recently. 0. I Q0 I 0 I 0 I I I - - - I -, >= ~~~~~~~~~~ r k 140OHz _ I _ I..- vH I-I - i I w....m.20 L 10. 1L 14. L I. 20. 2. 24. Tin of Day (hour) Because of absorption by water, MMW backscattering from snow is quite sensitive to liquid water content (LWC), and this sensitivity depends on the radar frequency. When snow is frozen, the scattering by the snow ice particles is significant and consequently, the backcattering cross-section is large. As snow melts during the warum daytime hours and LWC increases, the radar backscattering cross-section decreases dnmadcaiy. Figure 4 shows the backscattering coefficient at 35, 94, san 140 GHz as a function of time. Figure 5 is the air temperature and LWC measured by a freezing calorimeter a the same time. The sharp increase of LWC in snow is clearly detemed by MMW radar Among the three channels, the 140 GHzdaa is least sensitive to LWC. A detailed theoretical analysis and more extensive experimental results uc given in References 4 and 5. Fig. 4. Backscattering variation measured at 35, 94, and 140 GHz on 1 March, 1990. The incident angle was 40'. U LqU WtW Comte (%al.! S. T. -. 'IC)........ Air Tq mlpr ( -..... l w CmA m,. 10. 12. 14. 16. It 20. 22. 24. Thm of Day ( hour ) Fig. 5. Measured air temperature, snow surface temperature, and liquid water content of the top 5-cm layer on I March, 1990.

Pol'imneni Observations of'ftr MMW scattering by tree is usually dominated by the tree canopy during the spring and summer time. The polarization of the scatteed wave may contain information which can be used for extracting tree characteristics such as leaf size and density. The 35 and 94 GHz radars at the University of Michigan are capable of obtaining fully polarimetric information including the Mueller matrix, phase statistics, and degree of polarization. Measurements were conducted for leafy evergreen and coniferous trees. The first was a rhododendron (rhodo) with a planophil leaf orientation distribution and an average leaf area of 40 cm2. The second was a spruce blue shiner (spruce) with a uniform leaf orientation distribution and an average needle length of 1.5 cm To obtain the statistical data, more than 100 independent samples were measured for each incidence angle. Table 1 shows the degree of polarization for 6 different incident polarizations at 35 and 94 GHz. The degree of polarization is quite sensitive to frequency and tree type. The MMW return from the spruce tree, for example, shows a significant change in m between 35 and 94 GHz. The dependence on incidence angle is also clear for the planophil type tree. 2 tO 0* ehodo., Ov... 3s OGH, SO0 O rhodo.,,.... 35 OHt:. O 1$ 11 Is ' I0. S * spruce, 0,".r. 35 GH:, 30' O fhodo.,,,.*. 35 oHz, 30' I n I. f L nLhiE L Ii IIJI l JI...L.. a I I a I 1 11 1 1 1 n J {n t1ohl i II- U. m..| I * I la 0 rhodo., wiM. 34 G"z, 50 0 rhode.. Ov.w 315 GoI, S0O 0 - -II- - - - -1 - WILD I~, a TABLE 1. Degree of polarization m of rhod and spruce trees at 35 and 94 GHI. - - "= = z * a S T A O: $ 9 2.:: X Z;; Phase DitItroncO Ovv*hh (Degoru) Fig. 6. Distributions of vv -hh at 35 and 94 GHz. POL 35 GHz RHODO 30' 50~ 90~ SPRUCE 30~ 90' V 0.86 0.89 0.80 0.47 0.60 45 0.80 0.85 0.66 0.48. 0.43 H 0.76 0.83 0.69 0.37 0.57 135 0.83 0.87 0.69 0.49 0.49 LHC 0.69 0.77 0.61 0.1 0.11 RHC 0.61 0.79 0.51 0.06 0.21 94 GHz RIIODO 50 90" 0.92 0.82 0.91 0.87 0.87 0.87 0.82 0.77 0.81 0.73 0.77 0.73 SPRUCE 90" 0.64 0.71 0.71 0.58 0.39 0.39 Figure 7 is the co-polarized polarization signatures of rhododendron and spruce. The pedestal, which is related to depolarization, is much larger for the spruce tree than for the rhododendron tree. Figure 6 shows the phase distribution of the received signals. In all cases, *w -hv had a uniform distribution, suggesting that the co- and cross-pol. components are decorrelated. The most interesting case is the phase difference between the two co-polarized returns shown in Fig. 6. The distribution of vv -hh seems to depend on tree type and on the incidence angle. *vv -hh of rhododendron has a much narrower spread than that of spruce. The phase ~vh -hv of a truly monostatic system should be a delta function centered at zero degrees. Our data, however, shows a wide distribution. We think this is caused by the quasi-bistatic configuration of the radar system (small separation between transmit and receive antennas). Fig. 7. Co-polarized signatures at 35 GHz and 90' incidence angle. CONCLUSION MMW technology has progressed significantly in recent years and the cost of making MMW radars has been

decreasing. To utilize MMW radars for target detection and identification, it is important to understand how millimeterwaves interact with geophysical media. We believe the MMW Clutter Handbook published at the University of Michigan will be a useful source of data for engineers and scientists working with MMW radars (I[. Fm=[li, 12, U, VIT, I=I-Ev i2 12=lEh 12, U=2Rc(EvEh*), and V=2Im(EvEh*) A-3 Then (A-1) in terms of Stoke's vector becomes Fm, = (1i/r2) Lm Fmt A4 [1] F. T. Ulaby and T. F. Haddock, "Millimeter-wave Radar Scattering from Trrain:Data Handbook Version 2," Technical Report 026247-3-T, University of Michigan, 1990. [2] F.T. Ulaby, M. Whit and K. Sarabandi,"AVNA-Based Polarimetric Scatterometers", IEEE AP Mag., Vol 32, 1990. [3] J. B. Mead, Polarimetric Measurements of Foliage and Terrain at 225 Gllz", Ph. D. thesis, University of Massachusetts, Amherst, MA, 1990. (4] Y. Kuga, F. T. Ulaby, T. F. Haddock, and R. DeRoo, "MMW Radar Scattering from Snow: Part l-Radiative Transfer Model," to be published in Radio Science, 1991. [5] F. T. Ulaby, T. F. Haddock, R. Austin, and Y. Kuga, "MMW Radar Scattering from Snow: Part 2-Radiative Transfer Model," to be published in Radio Science, 1991. (6] F.T. Ulaby and C. Elachi, Radar Polarimetryfor Geoscience Applications, Artech, 1990. APPENDIX I COHERENT.ON-RECEIVE TECHNIOUE The complete polarimetric data can be obtained by measuring either the scattering amplitude matrix S or the Mueller matrix Lm. The ment of the cattering amplitude matrix requires curate phase measurement. Also 4 elements of S must be obtained within the decorrelation time of the target which is in the order of millisecond at the MMW frequency. The scattered electric field Er in tems of the scattering amplitude matrix S and the incident electric field Et is given by ZR e (S$.S) 2R e(SS,,) 21 m(S.,S,) 21 m(S,S,) Rc(S bs,) Re(S,,S + Re(SV.S, 4 S,,S,) Im(S,,S,, 4 SSSW).Im(s,$,S,9.Lm(S bSb,) * I M(SV, * S.01: A-5 where Lm is called a Mueller matrix. The totally incoherent method does not require the phase measurement. With this method the 4 elements of Stoke's vector are obtained by receiving 6 polarizations (V, H, 45,135, LHC, RHC) simultaneously. For example, the third element of Stoke's vector U is given as a ratio of intensities at 45 linear to 135 linear. For a given incident polarization, therefore, we can obtain a column of the Mueller matrix. To get the complete Mueller matrix, we need to repeat this process for 4 independent incident polarizations. Altogether at least 24 magnitude only measurements are required to obtain the complete Mueller matrix. Altllough the phase measurement is not required with the incul;rent method, it is necessary to receive all 6 polarizations simulutneously. The Stoke's vector in terms of 6 polarizations and a set of 4 independent incident polarizations are shown below. W~ I1..~ Wh + Wv Wh 12.... Wh + W; U. ws W45 -W3 W4s + W,35 A-6 A-7 A-8 A-9 WuHC WRHC WU4. + WnRC r t I. Shy SM E '-Ehl t:]t [S vSh A-1 A-2 I 0 1/2 1/2 IVwe 0 I {1 145 '/2 ILUI 1'2 A-l0 whet W* is a received intensity of * polarization. If a receiver is able to measure phase between V and H channels, it is possible to do the incoherent method without measuring 6 polarizations. This method is called the coherent-on-receive (COR) technique 12,3]. The elements of the Stoke's vector can be expressed as A1 1 To obtain S. we need to send [Ev,,0] and [OEh], and measure Ev and Eh simultaneously. Because a long time is required to sweep the frequency point the scatterometer based on the NWA is not suited for the fast acquisition of scattering matrix. For example, sweeping 401 points takes 400 msec and if we measure VV, HV, VH, and HH sequentially, it takes at least 2 seconds to acquire 4 elements of S. Obviously this technique does not work at the MMW frequency. The polarized wave can also be xpressed in terms of the Stoke's vector Fm which is defined as [61 -%1 I 12. I E,12 12.IEI U 21 EvlIIEhlICOS V 2IEvI Eh|sin & A-12 A-13 A-14

where 6 is a phase difference between V and H. The requirements of the COR method are (1) the transmitter must be able to transmit at least 4 independent polarizations and (2) the receiver must be able to get V and H channels simultaneously. The measurement time between the four incident polarizations can be much longer than the decorrelation time of target because correlation between V and H is obtained directory. Since it is relatively easy to obtain the phase difference between V and H at the MMW frequencies, we chose the COR method for our radars. Also COR is suited for a radar system based on NWA as a receiver. Using the elements of L,, the phase difference of the co- and cross-polarized components of S can be written as vw -hh = tan1l[(L43-L34)/(L33+L44)] A-15 wvv-hv = tan-'(L41/31) A-16 4vh -hv = tan-11 [(43+L34)/(L33-L44)], A-17 The degree of polarization, m, which is directly obtained from the received Fm, is given by m= 4((<Ijp> - <2p)2+<Upl'2+<Vp!2) / (<IIp>+<12P>) A-18 whereo is an ensemble average, and "P" denotes the transmitted polarization.

DERIVATION OF PHASE STATISTICS FROM THE MUELLER MATRIX K. Sarabandi Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan AbstractTo answer the question of what radar polarimetry has to offer to the remote sensing of random media, statistics of the phase difference of the scattering matrix elements must be studied. Recent polarimetric measurements of rough surfaces have indicated that the statistical parameters of the phase difference (mean, standard deviation, etc.) are very sensitive to some of the physical parameters. In this paper the probability density function of the phase differences is derived from the Mueller matrix assuming that the elements of the scattering matrix are jointly Gaussian. It is shown that the probability density functions of the co-polarized and cross-polarized phase differences are similar in form and each can be determined by two parameters (a and () completely. The expressions for the probability density functions are verified by comparing the histograms, the mean, and the standard deviations of phase differences derived directly from polarimetric measurements of variety of rough surfaces to the probability density function and its mean and standard deviation derived from the Mueller matrices of the same data. The expressions for the probability density functions are of special interest for non-coherent polarimetric radars and non-coherent polarimetric models for random media such as vector radiative transfer. 1

1 Introduction In the past decade substantial effort within the microwave remote sensing community has been devoted to the development and improvement of polarimetry science. Polarimetric radars are capable of synthesizing the radar response of a target to any combination of the receive and transmit polarizations from coherent measurements of the target with two orthogonal channels. Polarimetric radars have demonstrated their abilities in improving point-target detection and classification [Ioannidis and Hammers, 1979]. That is, for a point target in a clutter background the transmit and receive polarizations can be chosen such that the target to clutter response is maximum. Also, different point targets in the radar scene can be classified according to their optimum polarization. Although radar polarimeters have shown a great potential in point-target detection and classification, their capabilities in remote sensing of distributed targets is not completely understood yet. Considering the complexity involved in designing, manufacturing, and processing the data of an imaging polarimeter as opposed to a conventional imaging radar, it is necessary to examine the advantages that the imaging polarimeter provides about the targets of interest. For example, in retrieving the biophysical parameters from the polarimetric radar data one should ask whether there exists a dependency between the parameters and the measured phase of the scattering matrix components. If the answer is negative, obviously gathering polarimetric data for inversion of that parameter is a waste of effort. One way of confirming this question is by collecting data over a range of the desired parameter while keeping other influential parameters constant. This procedure, 2

if not impossible, is very difficult to conduct because of problems in repeatability of the experiment and difficulties in controlling the environmental conditions. Moreover at high frequencies (millimeter-wave frequencies and higher) coherent measurement of the scattering matrix is impossible because of instabilities of local oscillators and relative movements of the target and the radar platform [Meads and McIntosh, 1991]. At these frequencies non-coherent radars are employed which provide the Mueller matrix of the target. Another approach to examine the dependency of the radar response to the desired parameters of the targets is the application of theoretical models. One of the most successful polarimetric models for random media is the vector radiative transfer theory [Tsang et al., 1985]. This model is based on conservation of energy and the single scattering properties of the constituent particles. The solution of the radiative transfer equation relates the scattered-wave Stokes vector to the incident-wave Stokes vector via the Mueller matrix. The Mueller matrix, as computed by this method, is an ensembleaveraged quantity because of the inherent nature of the radiative transfer theory. Since the Mueller matrix is related to the scattering matrix through a nonlinear process and the components of the scattering matrix are statistically dependent, the information about the phase difference of the scattering matrix components cannot be obtained from the Mueller matrix directly. To achieve information about the phase statistics, one may resort to the Monte Carlo-type models which are computationally inefficient and in general inaccurate. Experimental observations of phase difference statistics from a polarimetric SAR at L-band [Ulaby et al., 1987; Zebker et al., 1987] over agricultural terrain and bare soil 3

surfaces indicate that the statistics of the co-polarized phase difference depends on the target type and its conditions. Recent measurements of bare soil surfaces by polarimetric scatterometers show that the variance of the co-polarized phase difference is a function of the roughness parameters and incidence angle but is less sensitive to moisture content [Sarabandi et al., 1991]. In view of difficulties in measuring the scattering matrix at high frequencies and performing controlled experiments, it is necessary to establish a relationship between the Mueller matrix and the statistics of the phase differences of the scattering matrix elements. In the next section we derive the probability density function of the co- and cross-polarized phase difference in terms of the Mueller matrix elements assuming that the scattering matrix elements are jointly Gaussian. Then the assumptions and final results are compared with the experimental data acquired by polarimetric scatterometers in Section 3. 2 Theoretical Derivation of Phase Difference Statistics The polarimetric response of a point or distributed target can be obtained by simultaneously measuring both the amplitude and phase of the scattered field using two orthogonal channels. If the incident and scattered field vectors are decomposed into their horizontal and vertical components, the polarimetric response can be represented by the scattering matrix S, which, for plane wave illumination we can write eikr SV, Svh ~Es=_ S h E] (1) r Shv Shh Eb~~ 4

where r is the distance from the radar to the center of the distributed target. It should be noted that in the backscattering case reciprocity implies Svh = Shve Each element of the scattering matrix, in general is a complex quantity characterized by an amplitude and a phase. When the radar illuminates a volume of a random medium or an area of a random surface, many point scatterers contribute to the total scattered energy received by the radar and therefore each element of the scattering matrix may be represented by N Spq = jSpei p,= v, h (2) n=1 Here N is the total number of scatterers each having scattering amplitude ISqn I and phase npq. It should be mentioned that the phase of each scatterer, as given in (2), includes a phase delay according to the location of the scatterer with respect to the center of the distributed target. Without loss of generality all multiple scattering over the surface or in the medium can be included in (2). Since the location of the scatterers within the illuminated area (volume) is random, the process describing the phasor Spq is a Wiener process (random walk) [Davenport, 1970]. If N is large enough, application of the central limit theorem shows that the real and imaginary parts of the scattering matrix element Sp, are independent identically distributed zero-mean Gaussian random variables. Equivalently it can also be shown that ISpq and kpq are, respectively, Rayleigh and uniform independent random variables. The three elements of the scattering matrix, in general, can be viewed as a six-element random vector and it is again reasonable to assume that the six components are jointly Gaussian. Observation of polarimetric data for a variety of distributed targets such as bare soil surfaces and different kinds of vegetation-covered terrain all indicate that the cross 5

polarized component of the scattering matrix (Shv) is statistically independent of the copolarized terms (S,, and Shh). Therefore the statistical behavior of Shy can be obtained from a single parameter, namely the variance (a2) of the real or imaginary part of Shy = X5 + iX6, that is 1 x52 + x6 fxs,x6 (5, 6) = 27rcrc2exp[- 2c2 or equivalently the joint density function ISvhl and q,h is 1 ISuiI2 fisvhl,O.h(lsvhIvh) = 2 Ira2Ishlex[- 2-2- ] ' (3) which indicates that Oh is uniformly distributed between (-xr, +i). Since measurement of the absolute phase of the scattering matrix elements is very difficult, it is customary to factor out the phase of one of the co-polarized terms, for example S,,, and therefore the phase difference statistics are of concern as opposed to the absolute phases. Since Sh, is assumed to be independent of Sv,, (not a necessary assumption) and both 4h, and q,, are uniformly distributed, it can be easily shown that the cross-polarized phase difference Xc = kvh - Ov is also uniformly distributed between (-r, +7r). The co-polarized elements of the scattering matrix, however, are dependent random variables which can be denoted by a four-component jointly Gaussian random vector X. Let us define Svv = X1 + iX2 Shh = X3 + iX4 and since X1, — X4 are Gaussian their joint probability density function can be fully determined by a 4 x 4 symmetric positive definite matrix known as covariance matrix 6

(A) whose entries axe given by [2] Aij = Aji =< XiXj > i,j E {1,...,4} The joint probability density function in terms of the covariance matrix takes the following form: fX(XI,,4) = 4 lA exp - A-X (4) where X is transpose of the column vector X. To characterize the covariance matrix the following observations are in order. First, it was shown that the real and imaginary parts of the scattering matrix elements are mutually independent and identically distributed zero-mean random variables, therefore X11 = X22 =< X1 >=< X2 >, (5) A12 =< X1X2 >= 0, (6) A33 = A44 =<X2>=< X2 >, (7) X34 =< X3X4 >= 0. (8) Second, it was shown that the absolute phase pp, is uniformly distributed and is independent of ISppI. Thus the random variable 0,, + khh is also uniformly distributed and is independent of ISvvIlShhl from which it can be concluded that < SvIllShhl CoS(vv, + khh) >= 0 (9) < Isv llShh Sin(v,, + 4hh) >= 0 In fact, the complex random variable S,,Shh is obtained from a similar Wiener process 7

which led to the random variables S,,v and Shh- On the other hand X1X3 - X2X4 = ISvvllshhl Cos(O+ + hh), (10) X1X4 + X2X3 = ISvvllShhlsin(qvv + qhh) In view of (9) and (10) it can easily be seen that 13 = 24, (11) A14 = -A23 ~ (12) The properties derived for the entries of the covariance matrix, as given by (5)-(8) and (11)-(12), indicates that there are only four unknowns left in the covariance matrix. The unknowns, namely A11, A13, A14, and A33, can be obtained directly from the Mueller matrix of the target as will be shown next. The Mueller matrix relates the scattered-wave Stokes vector to the incident-wave Stokes vector by [van Zyl and Ulaby, 1990] 1 F = MFi, r where F'i3 are the modified incident- and scattered-wave Stokes vector defined by IEVl2 IEhl2 F = 2R[EvEZ] T2s[ErEh e ] e The Mueller matrix can be expressed in terms of the elements of the scattering matrix 8

as follows [6] IStvl2 IS|l Shl2 S[hSvvh - Z[SVohStv] ISh l2 ISahhl2?[ShhShv] -Q[ShhShv] M2s[SvvShv] 2R[ShvShh] R[SvvShh + ShvS I] -Z[SvvSh -ShvSh] 2Qs[SvvSZ ] 2Z[ShvS*h] S[SVVShh + ShvS*h] R[SvvShh - ShvSvh] In the case of a random medium we are dealing with a partially polarized scattered wave and the quantity of interest is the ensemble averaged Mueller matrix. Using the assumption that the co- and cross-polarized terms of the scattering matrix are independent and employing the properties given by (5)-(8) and (11)-(12), the Mueller matrix in terms of the entries of the covariance matrix is given by 2A11 2a2 0 0 2a2 2A)33 0 0 M =< M >=. (13) O 0 2A13 + 2c2 2A14 O 0 -2A14 2A13 - 2a2 Equation (13) provides enough equations to determine the unknown elements of the covariance matrix and variance of the cross-polarized component, i.e. a2 = M.12 C11 2 All- = il, 1\33 = 77 2 2 A13 _= L.I+M44 A14 M,4-M.4a 4 A With the covariance matrix, the joint density function of X1,.., X4 can be obtained as given by (4). Using a rectangular to polar transformation, i.e. X1 = Pi cos 4vv, x2 = pl sin 4,, X3 = P2 COS (qhh, X4 = P2 sin4 hh 9

the joint probability density function of the amplitudes and phases takes the following form 1 ri fPlP2,Ov,,vhh (Pi P2,7 vv, Ohh)= 4r2apl12eXp t —[alp2 + a2p2 - 2a3plp2]j, (14) where A= ill = (A11i-A3_ - A 24)2 al= A33/V. a2 = A1/VI ' a3 = [A13 coS(khh - kvv) + A14 sin(4hh - kvv)] /VA To obtain the co-polarized phase difference statistics the joint density function of Ov, and khh is needed which can be obtained from roo roo /fvvOhh (4vv,4 hh) = I fPll2,Pvv,,hh(Pl,P2, vv,4 hh)dpldP2 * (15) Noting that al is a positive real number, the integration with respect to Pi can be carried out which results in fv., hh(vv.,hh) = 4 a r- p2e 2 P2dp2 + a3 J P22 [1erf( Ja- )] e- 2- al 2 dP2 (16) where erf(-) is the error function and the plus or minus sign is used according to the sign of a3. To evaluate the integrals in (16), we need to show that both a2 and a1a2 - a3 are positive numbers. By definition, a2 is positive and to show ala2 - a3 is positive we note that A is a symmetric positive definite matrix, therefore its eigenvalues must be positive. It can be shown that A has two distinct eigenvalues 7y and 72 each with multiplicity 2 and their product is given by 7172 = A1iA33 - A23 - A4 > 0. 10

Thus ala2 - a3 = y172 + [A13 coS(qhh - 0,,)- A14 sin(4hh - vv)12 is positive. After integrating the first integral and the first term of the second integral in (16) directly and using integration by parts on the second term of the second integral, (16) becomes fl',, ~hh ('vv,,hh) = 1 + 4rt / aa1 a2 ala2(ala2 - a3) x/ii~a3I ~00 Jasl -a2-p2 + -rlaa [ erf( [a2l )e-2-l(, 2-al)p } + 8v"~'~(ala2 - a2) o ' 8/d By expanding the error function in terms of its Taylor series, interchanging the order of summation and integration, then using the definition of the Gamma function it can be shown that j erf( P2)e 21 (a2-a)Pdp2 = 2r(al a2) tan-i ( la31 J7r(ala2 -a3) 2aa2-a The joint density function of kv, and qhh is a periodic function of 0 = qhh - 0,v and therefore the random variable X, after some algebraic manipulation, can be shown to have the following probability density function over the interval (-ir, +ir) f*(2) = "2 (A)- AD) {12 - D2 + tDan- "X D ] (17) where we recall that D = Al3 Cos + A14 sin q and the elements of the covariance matrix in terms of the Mueller matrix elements are 11

given by -1 1 -22 2 ' 2 ' 1=M33 + M44 M34 - M43 A13 = ~14 = 4 414 Some limiting cases can be considered in order to check the validity of (17). For example, when S,, and Shh are uncorrelated, then both A13 and A14 are zero for which fD(q) = 1/(2r), as expected. Also, for the case of completely polarized scattered wave where S,, and Shh are completely correlated, the determinant of A is zero and so f<(q) is a delta function. It is interesting to note that the p.d.f. of the phase difference is only a function of two parameters defined by a= + =t2 14 A11 A33 A13 where a and C can vary from O to 1 and -r to ir respectively. In fact if the wave were completely polarized, ( would have been the phase difference between the co-polarized terms. The parameter ( will, henceforth, be referred to as the polarized-phase-difference. In terms of these parameters (17) can be written as 1 - a2 2w [1 - a2 cos2(0 - C)] a cos( ) 7 + ta- a cos(-() (18) 1 ~1-a~ 2(q5- ) [2 V tan-' JI~ —;E~ 2~`;;p1 - a2 os2(CJ-J (') It can be shown that the maximum of the probability density function occurs at + = ( independent of a. However, the width of the p.d.f. (e.g. the 3 dB angular width) is only a function of a which will be referred to as the degree of correlation. The probability distribution function given by (18) is the analog of Gaussian distribution for periodic random variables where ( and ca are, respectively, the counterparts of the mean and 12

variance for Gaussian random variables. Figure 1 shows the p.d.f. for different values of ~ while keeping a constant, and Fig. 2 shows the p.d.f. for a fixed value of a while changing C as a parameter. The calculated mean and standard deviation of the phase difference as a function of both the polarized-phase-difference and the degree of correlation are depicted in Figs. 3 and 4 respectively. Lastly, it is necessary to point out that the formulation of the co-polarized phase difference p.d.f., as given in (17), is not restricted to the backscattering case or to the co- and cross-polarized components being uncorrelated. In fact we can derive the crosspolarized phase difference statistics in a similar manner and the p.d.f. in this case for the backscattering case can be obtained from (17) upon the following substitution for the elements of the cross-polarized covariance matrix A11 =- 11 A = Ml2 A13 -- A14=,14 2 14 2 3 Comparison with Measurements Using the polarimetric data gathered by scatterometers from a variety of natural targets, the assumptions leading to the probability density function of phase differences as derived in the previous section are examined. Also by generating the histograms, means, and standard deviations of the phase differences from the data and comparing them with the results based on the p.d.f. derived from the measured Mueller matrices validity of the model is also examined. The polarimetric radar measurements of bare soil surfaces were performed at L-, C-, and X-band frequencies for a total of eight different soil surface conditions (four roughness and two moisture conditions). For this experiment we tried to 13

preserve the absolute phase of the measured scattering matrix by calibrating the surface data with a metallic sphere located at the same distance from the radar as the center of the surface target. For each frequency, surface condition, and incidence angle a minimum of 700 independent samples were collected. The detailed procedure of the data collection and calibration is given in reference [4]. By generating the histograms of the real and imaginary parts of the elements of the scattering matrix for all surfaces, it was found that they have a zero-mean Gaussian distribution as we assumed. Figure 5 represents a typical case where the histogram of the real and imaginary parts of S,, and Shh of a dry surface with rms height 0.32 cm and correlation length 9.9 cm at C-band have a bell-shaped distribution. The properties of the covariance matrix as given by (5)-(8) and (11)-(12) are verified by calculating the covariance matrices of the data for all cases. Table 1 represents a typical situation where the covariance matrix for the same rough surface at C-band possesses the mentioned properties approximately, that is A11; A22, A12 t A34 " 0, A33 A44, A13 A24, and A14 e -A23. The small discrepancies are due to the fact that the measurement of the scattering matrix with absolute phase has an uncertainty of ~30 degrees. Table 2 gives the Mueller matrix of the typical surface (Table 1) at C-band from which the co- and cross-polarized phase difference probability density functions are calculated using (17) and are compared with the measured phase histograms in Figs. 6 and 7 respectively. Similar comparisons were also made for the rest of surfaces, frequencies, and incidence angles and it was found that the expression (17) predicts the density functions very accurately. Some example of these comparisons are shown in Figs. 8 and 9. Figures 8 and 9 compare the mean and standard deviation of the co-polarized phase 14

difference versus incidence angle at L- and X-band for a surface with rms height 0.4 cm and correlation length 8.4 cm in dry condition using the results based on the direct calculation and the results derived from (17). 4 Conclusions Prompted by the experimental observations which show strong dependence of phase differences of the scattering matrix elements on the physical parameters of random media, the statistical behavior of the phase differences for distributed targets is studied. The probability density functions of the phase differences are derived from the Mueller matrix of the target. In derivation of the density functions it is assumed that the real and imaginary parts of the co- and cross-polarized terms of the scattering matrix are jointly Gaussian and their covariance matrices are found in terms of the Mueller matrix elements. The functional form of the co- and cross-polarized density functions are similar and are obtained independently. It is shown that the density function of the phase difference is completely determined in terms of only two parameters. The assumptions and final expressions are verified by using a set of polarimetric data acquired by scatterometers from rough surfaces. Acknowledgement This research was supported by NASA Grant NAGW-2151 and ARO contract DAAL 03-91-G0202. 15

References [1] Davenport, W. B., Probability and random processes, New York: McGraw-Hill, 1970. [2] Ioannidis, G.A. and D.E. Hammers, "Optimum antenna polarization for target discrimination in clutter," IEEE Trans. Antennas Propagat., vol. 27, no. 5, 1979. [3] Meads, J. B., and R. E. McIntosh, UPolarimetric backscatter measurements of deciduous and coniferous trees at 225 GHz," IEEE Trans. Geosci. Remote Sensing, vol. 29, no. 1, 1991. [4] Sarabandi, K., Y. Oh, F.T. Ulaby, "Polarimetric radar measurement of bare soil surfaces at microwave frequencies," Proc. of IEEE Geosci. Remote Sens. Symp., Espoo, June 1991. [5] Tsang, L., J.A. Kong, and R.T. Shin, Theory of Microwave Remote Sensing, New York: John Wiley and Sons, 1985. [6] Ulaby, F. T., D. Held, M.C. Dobson, K.C. McDonald, and T.B.A. Senior, "Relating polarization phase difference of SAR signals to scene properties," IEEE Trans. Geosci. Remote Sensing, vol.25, no. 1, 1987. [7] van Zyl, J.J., and F.T. Ulaby, Radar Polarimetry for Geoscience Applications, Chapter 2, ed. F.T. Ulaby and C. Elachi, Norwood, MA: Artech house, 1989. [8] Zebker, H.A., J.J. van Zyl, and D.N. Held, "Imaging radar polarimetry from wave synthesis," J. Geophys. Res., vol. 92, no. B1, 1987. 16

1.00 0.03 0.75 -0.12 0.03 0.90 0.08 0.68 0.75 0.08 0.77 0.05 -0.12 0.68 0.05 0.69 Table 1: Normalized covariance matrix of co-polarized terms of scattering matrix for a surface with rms height 0.3 cm and correlation length 9 cm at C-band and at 30 degrees incidence angle. 1.000 0.030 0.000 0.000 0.028 0.767 0.000 0.000 0.000 0.000 0.770 -0.11 0.000 0.000 0.110 0.711 Table 2: Normalized Mueller matrix for a surface with rms height 0.3 cm and correlation length 9 cm at C-band and at 30 degrees incidence angle. 17

List of Figures 1 The probability density function of the co-polarized phase difference for a fixed value of a (degree of correlation) and five values of ( (coherentphase-difference)................................. 18 2 The probability density function of the co-polarized phase difference for a fixed value of ( (coherent-phase-difference) and four values of a (degree of correlation)............................ 19 3 The mean value of the co-polarized phase difference as a function of a (degree of correlation) and ( (coherent-phase-difference)........... 20 4 The standard deviation of the co-polarized phase difference as a function of a (degree of correlation) and C (coherent-phase-difference)........ 20 5 The histogram of the real and imaginary parts of S,, and Shh for a rough surface with rms height 0.32 cm and correlation length 9.9 cm at C-band and 30~ incidence angle........................... 21 6 The histogram and p.d.f. of the co-polarized phase difference for a rough surface with rms height 0.32 cm and correlation length 9.9 cm at C-band and 30~ incidence angle........................... 22 7 The histogram and p.d.f. of the cross-polarized phase difference for a rough surface with rms height 0.32 cm and correlation length 9.9 cm at C-band and 30~ incidence angle........................ 23

8 Angular dependency of the mean of the co-polarized phase difference for a dry rough surface with rms height 0.4 cm and correlation length 8.4 cm at L- and X-band............................... 24 9 Angular dependency of the standard deviation of the co-polarized phase difference for a dry rough surface with rms height 0.4 cm and correlation length 8.4 cm at L- and X-band........................ 25

0.020 a).-7, d) a, P._;,m CO to q) Oz + 0.015 0.010 0.005 0.000................ -180. -135. -90. -45. 0. 45. khh -,vv(Degrees) 90. 135. 180. Figure 1: The probability density function of the co-polarized phase difference for a fixed value of a (degree of correlation) and five values of C (coherent-phase-difference). 18

61 (uoI.elalo o Jo oalaip) o Jo soUnlA snoj pus (aouaa.Ujp-aseqd-ualaqo:) ) jo anreA paxi e ioJ aualaJJ!p aised pazTelod-o- all Jo uo!punnJ Xl!suap XIn!qeqold aqL:Z ainj (s'aaoj) 'n0 -'1 -St.'0 gSv- '06 - '081 'g1 '06 *Sg1- '08... ' '' -i" " I... I... IL.-. -:-.:-.- _',_,~' — *,X\ ~/,, -- \ ~-0 / - 'I I\ ",.'1 I I / I I I I I I I 11 \ I \I.I. 0000 0000 01O'O 010,0 OZO'O

C -50 0 0.25 0.5 Figure 4: The mestandard deviation of the co-polarized phase difference as a function of a (degree of correlation) a nd d (coh erent-ph hase-diff erence). 180 -9 -ZO. 7 (degree of correlation) and C (coherent-phase-difference). 20

0.500 -.0 10 "0 Z, N 0.400 0.300 0.200 0.100 0.000 m-K-r-. I.I.. -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 Normalized Magnitude Distribution 3.0 4.0 Figure 5: The histogram of the real and imaginary parts of Sl, and Shh for a rough surface with rms height 0.32 cm and correlation length 9.9 cm at C-band and 300 incidence angle. 21

0.020 0.018 L | -— Measured(Mean=-5.00,Std. Dev=47.7) --------- Calculated (Mean=-7.8', Std. Dev=47.20) 0.016 0.014 1 0.012 0.010 0.008 0.006 0.004n 0.002 0.000 -180. -120. -60. 0. 60. 120. 180. Ohh - vv Figure 6: The histogram and p.d.f. ofthe co-polarized phase difference for a rough surface with rms height 0.32 cm and correlation length 9.9 cm at C-band and 300 incidence angle. 22

0.020.. 0.018, - Measured (Mean=-1.00, Std. Dev=102.5~) ---------- Calculated (Mean=0.8~, Std. Dev=103.7) 0.016.... 'm 0.014 0.012 = 0.010 -, 0.008:g 0.006 0.004 0.002 0.000 -180. -120. -60. 0. 60. 120. 180. hh - 'kvv Figure 7: The histogram and p.d.f. of the cross-polarized phase difference for a rough surface with rms height 0.32 cm and correlation length 9.9 cm at C-band and 300 incidence angle. 23

60. 50.1 - X-band (Calculated) O X-band (Measured).-. ---... L-band (Calculated) o L-band (Measured) CA o L as I a) e O E lTo 40. 30. 20. 10. 0. B.f O M Ul -10. - -20.1 -30. 0. 10. 20. 30. 40. 50. 60. 70. 80. Incidence Angle (Degrees) Figure 8: Angular dependency of the mean of the co-polarized phase difference for a dry rough surface with rms height 0.4 cm and correlation length 8.4 cm at L- and X-band. 24

90. 80. 70. 0 "I co e t) 4.) 3 05 O 0 CC m~ 4-i 60.. 50. l 40.1 - 30.. 20.1 X-band (Calculated) 0 X-band (Measured) --------- L-band (Calculated) o L-band (Measured) 10. O.0. 10. 20. 30. 40. 50. 60. Incidence Angle (Degrees) 70. 80. Figure 9: Angular dependency of the standard deviation of the co-polarized phase difference for a dry rough surface with rms height 0.4 cm and correlation length 8.4 cm at L- and X-band. 25

UNIVERSITY OF MICHIGAN 3 1111111015 02527 80971111111 3 9015 02527 8097