MODELING AND INVERSION OF THE RADAR RESPONSE OF VEGETATION CANOPIES by Paul F. Polatin.'~..... - o....:..1, ~. A: dissert'ation' suibiitted in partial fulfillment of the requir;eents for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1993 Doctoral Committee: Assistant Professor Kamal Sarabandi, Co-Chairman Professor Fawwaz T. Ulaby, Co-Chairman Professor Anthony W. England Professor William R. Kuhn Associate Professor Valdis V. Liepa

- G /'4tSr /2l

This thesis is dedicated to my Father, my Elder Brother and my Wife 11

ACKN OWLED GEMENTS The author would like to express his gratitude to the members of his committee and to the Radiation Laboratory for supporting his research. Special thanks are due to Professor Kamal Sarabandi for his constant advice and encouragement and to Professor Fawwaz T. Ulaby for his forbearance and generous financial support throughout the course of this work. He would also like to thank the following friends and colleagues for their help, friendship and encouragement during his graduate studies: Dr. Leland Pierce, Professor C.T. Tai, Dr. Norman Vandenberg, Dr. Kyle McDonald, Dr. Hasnain Syed, Dr. Ahad Tavakoli, Roger DeRoo, Leo Kempel and Dan Ross.

TABLE OF CONTENTS DEDICATION.................................. ii ACKNOWLEDGEMENTS.......................... iii LIST OF FIGURES............................... vii LIST OF TABLES................................ xix CHAPTER I. INTRODUCTION........................ 1 II. THE SECOND-ORDER RADIATIVE TRANSFER MODEL OF THE CANOPY TRUNK LAYER.............. 9 2.1 Introduction......................... 9 2.2 The Radiative Transfer Model................. 13 2.2.1 Basic Equations and Definitions........... 13 2.2.2 Zeroth-Order Solution................. 20 2.2.3 First-Order Solution.................. 21 2.2.4 Second-Order Solution................ 24 2.3 Results and Comparison With Measurements......... 30 2.4 Conclusions............................ 38 III. MONTE CARLO SIMULATION OF SCATTERING FROM A LAYER OF VERTICAL CYLINDERS............ 40 3.1 Introduction........................... 40 3.2 Scattering from Two Adjacent Cylinders at Oblique Incidence 41 3.2.1 Second Order Interaction............... 46 3.2.2 Effect of Ground Plane............... 48 3.3 Monte Carlo Simulation............. 53 3.4 Comparison with Radiative Transfer.............. 57 3.5 Conclusions............................ 67 iv

IV. ELECTROMAGNETIC SCATTERING FROM TWO ADJACENT OBJECTS......................... 70 4.1 Introduction........................... 70 4.2 Secondary Scattered Field from Reciprocity......... 72 4.3 Electromagnetic Scattering by a Cylinder-Sphere Pair.... 76 4.4 Numerical Results........................ 85 4.5 Conclusions............................ 92 V. SIMULATION OF ELECTROMAGNETIC SCATTERING FROM A HETEROGENEOUS MEDIUM...... 93 5.1 Introduction........................... 93 5.2 Construction of Scattering Terms................ 95 5.2.1 The Sphere-Cylinder Interaction........... 96 5.2.2 Effect of Ground Plane................ 100 5.3 Validation of the Scattering Terms............... 106 5.4 Heterogeneous Canopy Simulation............... 110 5.4.1 Comparison with Radiative Transfer......... 115 5.5 Conclusions............................ 127 VI. A HYBRID MODEL FOR ELECTROMAGNETIC SCATTERING FROM FOREST CANOPIES............. 129 6.1 Introduction.................. 129 6.2 The Transmission Matrix for the Trunk Layer......... 134 6.3 Construction of the Hybrid Model............... 140 6.4 Results and Discussion...................... 146 6.5 Conclusions............................ 166 VII. AN ITERATIVE INVERSION ALGORITHM WITH APPLICATION TO THE RADAR RESPONSE OF VEGETATION CANOPIES.............171 7.1 Introduction........................... 171 7.2 Iterative Algorithm....................... 173 7.3 Implementation of the Algorithm................ 182 7.4 Application to a Vegetation Canopy.............. 184 7.5 Inversion Results......................... 188 7.6 Error Analysis........................ 195 7.7 Conclusions............................ 199 VIII. CONCLUSIONS AND RECOMMENDATIONS....... 200 8.1 Summary............................. 200 v

8.2 Recommendations for Future Work............... 205 APPENDICES.................................. 207 BIBLIOGRAPHY.2..............................210 vi

LIST OF FIGURES Figure 2.1 Scattering into the Specular Cone................... 12 2.2 First-Order Effects......................... 13 2.3 Second-Order Effects.......................... 14 2.4 Closed Canopy Model............... 15 2.5 Trunk Region Problem Geometry................... 18 2.6 Data/Model Comparison L-Band hh Polarization (Mv=0.05).... 31 2.7 Data/Model Comparison L-Band vh/hv Polarization (Mv=0.05). 31 2.8 Data/Model Comparison L-Band vv Polarization (Mv=0.05).... 32 2.9 Data/Model Comparison L-Band hh Polarization (Mv=0.16).... 32 2.10 Data/Model Comparison L-Band vh/hv Polarization (Mv=0.16). 33 2.11 Data/Model Comparison L-Band vv Polarization (Mv=0.16).... 33 2.12 Data/Model Comparison L-Band hh Polarization (Mv=0.32).... 34 2.13 Data/Model Comparison L-Band vh/hv Polarization (Mv=0.32). 34 2.14 Data/Model Comparison L-Band vv Polarization (Mv=0.32).... 35 3.1 Geometry of a long cylinder with equivalent electric and magnetic currents on the surface...................... 44 3.2 Geometry of dual cylinder configuration and second-order interaction. 45 3.3 Integration contour in complex 7-plane................ 46

3.4 First order interactions of a cylinder with the ground plane.... 49 3.5 Second order interactions of two cylinders with the ground plane.. 50 3.6 Broadside and endfire cylinder configurations. Wave incident in X-Z plane with incidence angle fi...................... 52 3.7 RCS of two rods broadside on metal plate at 9.5 GHz: 2cm separation................................... 53 3.8 RCS of two rods endfire on metal plate (vv-pol.) at 9.5 GHz: 4cm separation................................ 54 3.9 RCS of two rods endfire on metal plate (hh-pol.) at 9.5 GHz: 4cm separation................................. 55 3.10 The measurement of a random ensemble of vertical cylinders..... 57 3.11 Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 70cyls/m2......................... 58 3.12 Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 100cyls/m2......................... 59 3.13 Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 140cyls/m2......................... 59 3.14 Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 180cyls/m2...................... 60 3.15 Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 70cyls/m2............... 60 3.16 Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 100cyls/m2................... 61 3.17 Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 140cyls/m2................... 61 viii

3.18 Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 180cyls/m2................... 62 3.19 Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 70cyls/m2................... 62 3.20 Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 100cyls/m2................... 63 3.21 Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 140cyls/m2................... 63 3.22 Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 180cyls/m2................... 64 3.23 Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 140cyls/m2....... 65 3.24 Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 140cyls/m2...... 66 3.25 Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 70cyls/m2........ 67 3.26 Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 70cyls/m2........ 68 4.1 The current induced on particle #1 by the source produces a scattered field from particle #2........................ 73 4.2 Particle #1 is removed and an elementary current source is placed at the observation point........................ 73 4.3 Geometry and coordinates of the cylinder-sphere pair........ 78 4.4 Scattering geometry and angles for forward specular cone scattering. 86

4.5 The vv-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at 13=l8cm, q = 1800 and i=9cm relative to the cylinder base. The incidence and scattering angles are 1430 relative to vertical. The frequency is 9.25 GHz................ 87 4.6 The vv-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, - = 450 and,=9cm relative to the cylinder base. The incidence and scattering angles are 1450 relative to vertical. The frequency is 9.25 GHz................ 88 4.7 The vv-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at 13=2.Scm, q = 450 and,=9cm relative to the cylinder base. The incident wave azimuth angle is 1800 and the scattered wave azimuth angle is 1300. The frequency is 9.25 GHz.. 88 4.8 The vh-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, b = 45~ and,=9cm relative to the cylinder base. The incident wave azimuth angle is 1800 and the scattered wave azimuth angle is 3500. The frequency is 9.25 GHz.. 90 4.9 The hh-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at /=2.0cm, q = 450 and,=9cm relative to the cylinder base. The incidence and scattering angles are 1450 relative to vertical. The frequency is 9.25 GHz................ 90 4.10 The hh-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, q = 450 and,=9cm relative to the cylinder base. The incident wave azimuth angle is 1800 and the scattered wave azimuth angle is 900. The frequency is 9.25 GHz.. 91 5.1 Geometry and coordinates for the sphere-cylinder interaction..... 96 5.2 Image theory geometry for construction of the cylinder-ground-sphere interaction............................... 101 5.3 Image theory geometry for construction of the ground-sphere-groundcylinder interaction.......................... 104 5.4 Scattering geometry and reference angles.............. 108 5.5 Backscatter pattern for a perfectly conducting sphere above a conducting ground plane. The sphere has ka = 1.69, where a is the sphere radius, and is located 9cm above the ground plane. The frequency is 9.25 GHz................................. 109

5.6 Backscatter pattern for a pair of vertical cylinders on a conducting ground plane. The cylinders are 18cm in length, 0.1cm in diameter and separated by 2cm in the endfire configuration. The Conductivity of the cylinders is 100 S/m. vv-polarized return at 9.25 GHz.... 109 5.7 Backscatter pattern for a sphere-cylinder pair above a perfectly conducting ground plane with p=l1.5cm, q = 180~ and i=9cm. vvpolarized return at 9.25GHz................... 111 5.8 Bistatic scattering pattern for a sphere-cylinder pair above a perfectly conducting ground plane with 1=1.5cm, b = 180~ and i=9cm. vvpolarized return at 9.25GHz. The elevation angle is 370....... 111 5.9 Bistatic scattering pattern for a sphere-cylinder pair above a perfectly conducting ground plane with )=2.0cm, X = 45~ and i=18cm. The azimuthal scattering angle is 350~. vh-polarized return at 9.25 GHz. 112 5.10 Bistatic scattering pattern for a sphere-cylinder pair above a perfectly conducting ground plane with p=1.5cm, q = 450 and z=18cm. The elevation angle is 520. vh-polarized return at 9.25 GHz....... 112 5.11 Backscatter pattern for a sphere-cylinder pair above a perfectly conducting ground plane with p=2.0cm, q = 45~ and Z=9cm. hhpolarized return at 9.25 GHz................... 113 5.12 Bistatic scattering pattern for sphere-cylinder pair above a perfectly conducting ground plane with p=2.0cm, - = 450 and i=9cm. The elevation angle is 35~. hh-polarized return at 9.25 GHz....... 113 5.13 Configurations for the pure two-layer and mixed two-layer canopy models.................. 116 5.14 Typical convergence behavior for the Monte Carlo simulation of the backscattering coefficient. The number density is 106.1 cylinders per square meter. vv-polarized return at 9.25 GHz............ 117 5.15 Typical convergence behavior for the Monte Carlo simulation of the phase statistic, a. The number density is 106.1 cylinders per square meter. The frequency is 9.25 GHz............. 117 5.16 Homogeneous layer of 14,147 conducting spheres per cubic meter distributed in 10cm layer above conducting ground. Spheres have ka - 0.62. vv-polarized backscatter cross section.............. 118

5.17 Homogeneous layer of 106.1 vertical dielectric cylinders per square meter above conducting ground. Cylinders are 5.6 wavelengths long and 0.17 wavelengths in diameter. Relative dielectric,- = 35 + ill. vv-polarized backscatter cross section................. 119 5.18 Homogeneous layer of 106.1 vertical dielectric cylinders per square meter above conducting ground. Cylinders are 5.6 wavelengths long and 0.17 wavelengths in diameter. Relative dielectric er = 35 + ill. hh-polarized backscatter cross section................. 119 5.19 Heterogeneous canopy consisting of 35.4 cylinders per square meter in lower layer and 14,147 spheres per cubic meter in an upper layer 10cm thick. The cylinders are 18cm high and 0.55cm in diameter. The spheres have ka = 0.62. vv-polarized backscatter coefficient at 9.25 GHz............................... 121 5.20 Heterogeneous canopy consisting of 106.1 cylinders per square meter in lower layer and 14,147 spheres per cubic meter in an upper layer 10cm thick. The cylinders are 18cm high and 0.55cm in diameter. The spheres have ka = 0.62. vv-polarized backscatter coefficient at 9.25 GHz............................... 122 5.21 Comparison of Monte Carlo simulation of vv-polarized backscatter coefficient for a homogeneous canopy and a heterogeneous canopy. The homogeneous canopy has 106.1 cylinders per square meter. The heterogeneous canopy has the same number density of cylinders but has 14,147 spheres per cubic meter in a 10cm layer on top...... 122 5.22 Heterogeneous canopy consisting of 106.1 cylinders per square meter in lower layer and 14,147 spheres per cubic meter in an upper layer 10cm thick. The cylinders are 18cm high and 0.55cm in diameter. The spheres have ka = 0.62. hh-polarized backscatter coefficient at 9.25 GHz............................. 123 5.23 Comparison of Monte Carlo simulation of vh-polarized backscatter coefficient for a homogeneous canopy and a heterogeneous canopy. The homogeneous canopy has 106.1 cylinders per square meter. The heterogeneous canopy has the same number density of cylinders but has 14,147 spheres per cubic meter in a 10cm layer on top...... 124 5.24 Homogeneous canopy consisting of 106.1 cylinders per square meter. The degree of copolarized phase correlation at 9.25 GHz....... 124 xii

5.25 Effect on the degree of copolarized phase correlation of adding a layer of 14,147 small conducting spheres per cubic meter in a layer 10cm thick on top of the canopy of 106.1 cylinders per square meter... 125 5.26 Comparison of pure and mixed layer R.T. models for a canopy consisting of 106.1 cylinders per square meter and 14,147 spheres per cubic meter. The crown layer is 10cm in height in both cases. vvpolarized return at 9.25 GHz. The pure and mixed copolarized Monte Carlo simulations are almost identical................. 125 5.27 Comparison of pure and mixed layer Monte Carlo simulations for a canopy consisting of 106.1 cylinders per square meter and 14,147 spheres per cubic meter. The crown layer is 10cm in height in both cases. vh-polarized return at 9.25 GHz................ 126 6.1 The vv-polarized backscattering coefficient of 106.1 vertical dielectric cylinders per square meter above a conducting ground plane at 9.25 GHz. The cylinders are 18cm long and 0.55cm in diameter and have a dielectric constant e7 = 35 + ill.................... 130 6.2 The vv-polarized backscattering coefficient of 35.4 vertical dielectric cylinders per square meter above a conducting ground plane at 9.25 GHz. The cylinders are 18cm long and 0.55cm in diameter and have a dielectric constant e, = 35 + ill............... 131 6.3 The vv-polarized backscattering coefficient of 14,147 metal spheres per cubic meter in a layer 10cm high above a conducting ground plane at 9.25 GHz. The spheres have a diameter of 0.6350cm.... 132 6.4 Geometry of the cylinder layer for transmission matrix derivation.. 133 6.5 Tree canopy layout for the hybrid model showing intensity designations................................... 141 6.6 The effective vv-polarized scattering matrix of a distribution of cylinders as a function of number density. The cylinders are 18cm long and 0.1cm in diameter and have a conductivity of 100 mhos per meter. The angle of incidence is 80~ from normal and the frequency is 9.25 GHz............................... 148 6.7 The effective vv-polarized scattering matrix of a distribution of cylinders as a function of number density. The cylinders are 18cm long and O.lcm in diameter and have a conductivity of 100 mhos per meter. The angle of incidence is 20~ from normal and the frequency is 9.25 GHz................................ 149 xiii

6.8 The magnitude of the correlation function for the z-component of the near electric field within a medium consisting of 177 thin conducting cylinders per square meter as computed by the method of moments. The cylinders are 5.5 wavelengths long and 0.03 wavelength in diameter at 9.25 GHz. The excitation is vertically polarized and incident at 30~ from the vertical axis................... 151 6.9 The vv-polarized forward scattering coefficient of a layer of 177 thin cylinders per square meter. The cylinders are 5.5 wavelengths long, 0.03 wavelength in diameter and have a conductivity of 100 mhos per meter at 9.25 GHz. The cylinders were positioned with a random number generator. The single scattering approximation provides a good estimate at this fairly high density of scatterers........ 152 6.10 The vv-polarized forward scattering coefficient of a layer of 264 thin cylinders per square meter. The cylinders are 5.5 wavelengths long, 0.03 wavelength in diameter and have a conductivity of 100 mhos per meter at 9.25 GHz. The cylinders were positioned with a random number generator. The second-order scattering approximation is needed here to provide an acceptable estimate of the scattering coefficient................................. 153 6.11 The vv-polarized intensity transmission matrix element as a function of cylinder height for a layer consisting of 30 cylinders per square meter. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants er = 35+ill at 9.25 GHz. The angle of elevation is 60~ from vertical for the incident wave......... 154 6.12 The vv-polarized intensity transmission matrix element as a function of cylinder height for a layer consisting of 90 cylinders per square meter. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants e, = 35+ill at 9.25 GHz. The angle of elevation is 60~ from vertical for the incident wave.......... 155 6.13 The hh-polarized intensity transmission matrix element as a function of cylinder height for a layer consisting of 90 cylinders per square meter. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants c, = 35+ill at 9.25 GHz. The angle of elevation is 60~ from vertical for the incident wave......... 156 xiv

6.14 The vv-polarized intensity transmission matrix for a layer of cylinders 5.5 wavelengths in height and 0.17 wavelengths in diameter. The cylinders have a relative dielectric constant of er = 35+ill at 9.25 GHz and the transmissivity of the layer has been computed using a first-order Monte Carlo simulation................... 157 6.15 The difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of a layer of cylinders above a conducting ground plane as a function of the real part of the relative dielectric constant of the cylinders. The imaginary part of the dielectric constant is zero in this case. The cylinders are 5.5 wavelengths long and 0.17 wavelengths in diameter......... 158 6.16 The difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of a layer of cylinders above a conducting ground plane as a function of the imaginary part of the relative dielectric constant of the cylinders. The real part of the dielectric constant is 20 in this case. The cylinders are 5.5 wavelengths long and 0.17 wavelengths in diameter........ 159 6.17 The difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of a layer of cylinders above a conducting ground plane as a function of the layer height. The relative dielectric constant of the cylinders is 6, = 35+ill in this case and the cylinder diameter is 0.17 wavelengths.......... 160 6.18 Area convergence behavior of the vv-polarized transmissivity for a layer of 130 thin cylinders per square meter. The cylinders are 5.5 wavelengths long, 0.03 wavelengths in diameter and have a relative dielectric constant 6, = 1+i1943 at 9.25 GHz. The distribution is circular............. 161 6.19 Comparison of models for the vv-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 3540 spheres per cubic meter. 162 6.20 Comparison of models for the hh-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 3540 spheres per cubic meter. 163 XV

6.21 Comparison of models for the vv-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 14,147 spheres per cubic meter. 163 6.22 Comparison of models for the hh-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 14,147 spheres per cubic meter. 164 6.23 Comparison of models for the vv-polarized backscattering coefficient of a layer of vertical cylinders 5.5 wavelengths high above a conducting ground plane. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants e, = 35+ill at 9.25 GHz. The number density is 35 cylinders per square meter........... 164 6.24 Comparison of models for the hh-polarized backscattering coefficient of a layer of vertical cylinders 5.5 wavelengths high above a conducting ground plane. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants c, = 35+ill at 9.25 GHz. The number density is 35 cylinders per square meter........... 165 6.25 Comparison of models for the vv-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 35 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter... 166 6.26 Comparison of models for the hh-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 35 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter... 167 6.27 Comparison of models for the vv-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 106 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter... 167 6.28 Comparison of models for the hh-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 106 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter... 168 xvi

6.29 Comparison of models for the vv-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. All the second-order interactions between spheres and cylinders and between cylinders alone have been included in the simulation while the hybrid model includes only the second-order coupling between cylinders................................. 168 6.30 Comparison of the hybrid model and a Monte Carlo simulation for the vh-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. All the second-order interactions between spheres and cylinders and between cylinders alone have been included in the simulation while the hybrid model includes only the second-order coupling between cylinders........... 169 7.1 The first four vv-polarized Fourier coefficients of a simulated vegetation canopy as a function of volumetric soil moisture at 1.5 GHz.. 177 7.2 The first four vh-polarized Fourier coefficients of a simulated vegetation canopy as a function of volumetric soil moisture at 1.5 GHz.. 177 7.3 Comparison of radiative transfer model with Fourier representations for vv polarization at 1.5 GHz.................. 178 7.4 Comparison of radiative transfer model with Fourier representations for vh polarization at 1.5 GHz.................. 178 7.5 2-D parameter space discretization showing subrange centroids as solid marks.......................... 180 7.6 Layered structure of the radiative transfer model for a vegetation canopy.................................. 185 7.7 The determination of soil moisture over the entire parameter range near the primary centroid..................... 190 7.8 The determination of trunk canopy density over the entire parameter range near the primary centroid.................... 190 7.9 The rmination of leaf canopy detnsity over the entire parameter range near the primary centroid................... 191

7.10 The determination of trunk height over the entire parameter range near the primary centroid......................... 191 7.11 The determination of leaf canopy density over a second level parameter sub-range.............................. 192 7.12 The determination of trunk height over a second level parameter subrange.................................. 192 7.13 Uniform convergence of the iterative inversion algorithm....... 193 7.14 Non-uniform convergence of the iterative inversion algorithm.... 194 7.15 Anomalous convergence of the iterative inversion algorithm. 194 7.16 Oscillatory non-convergence of the iterative inversion algorithm.. 195 7.17 Typical inversion error bounds for variable soil moisture with other parameters held fixed........................ 198 7.18 Typical inversion error bounds for fixed leaf density with variation of the soil moisture parameter.................... 198

LIST OF TABLES Table 7.1 Variable parameter set for the vegetation canopy used to test the inversion algorithm......................... 186 7.2 Fixed parameter set for the vegetation canopy used to test the inversion algorithm............................ 187 xix

CHAPTER I INTRODUCTION The objective of microwave remote sensing observations of the earth is to obtain information about the terrestrial environment from airborne and spaceborne platforms. The use of synthetic aperture radars as active sensors in this type of work has become increasingly important in the last few decades. The power and flexibility of this tool has increased greatly with the development of radar polarimeters which allow simultaneous measurement of amplitude and phase at four polarizations and several frequencies for distributed targets. Along with this increase in ability to measure such targets comes the need to be able to utilize the information so obtained. For this reason it is essential to understand both the quantitative and qualitative aspects of electromagnetic scattering from terrain. Because a large portion of the earth's surface is covered by vegetation having many different types of canopy configurations it is necessary to understand the detailed aspects of microwave interaction with such media. The overall purpose of modeling radar scatter from vegetation canopies is to be better able to interpret and invert signals from the remote radar sensors used in terrain mapping operations. With this objective in mind, the goal of this thesis is to examine both the modeling and inversion aspects of the interaction of electromagnetic waves with general

2 random media having structural properties closely related to those found in vegetation canopies. It is hoped that this investigation will make clear the areas in which methods currently in use for the analysis of radar scattering data from vegetation canopies have their limits of applicability. It is also hoped that this work will point the way toward a solution to these problems. Due to limitations of space and time it has been necessary to devote the greater part of the volume of this thesis to the modeling aspects of the problem. In fact, five chapters have been devoted to modeling the radar scattering properties of vegetation. Only Chapter 7 has been dedicated to the inversion of remotely sensed data. In Chapter 7 no attempt has been made to investigate the most popular technique for inverting sensor data which is based on the application of artificial neural networks. Instead, the focus has been on developing a simple, model-based, iterative inversion algorithm that permits the investigator to develop insight into the fundamental nature of both the forward and inverse aspects of the problem under consideration. Chapter 7 also attempts to provide a framework into which the modeling work fits by considering the sensor dependent nature of the remote sensing problem in terms of both sensor configuration and the usual experimental errors that are an integral part of any measurement process. The remaining chapters of this work are devoted to the various aspects of modeling radar scattering from vegetation canopies. In recent years considerable effort has been invested in the development of theoretical scattering models for forest canopies. This is because a large fraction of the global vegetation cover is in the form of forests and other trunk dominated canopies. A forest canopy is considered to be an inhomogeneous medium comprised of scattering elements with different sizes, shapes, and electrical properties [95, 96, 40, 39]. Except for the analytical wave approach

3 [81, 24, 44], which accounts for the particles in the medium through a fluctuating permittivity function, all existing methods account for the particles through the single scattering properties of the individual scatterers [45, 72, 42, 82]. The analytical wave approach is appropriate for media in which the ratio of the fluctuating dielectric function to the mean dielectric constant is small. Therefore, at microwave and millimeter wave frequencies, where the dielectric constant of leaves and branches is much larger than that of air, the analytical wave approach may not be appropriate. The solution of the scattering problem for a canopy consisting of discrete scatterers can be tackled in two ways: (1) the field approach [45] and (2) the intensity approach [96]. Due to the complexity of the methodology the solution based on the field approach is limited to sparse distributions of weakly scattering particles. The intensity approach, or the radiative transfer (RT) method, is very general, easily formulated and mathematically convenient. These are precisely the characteristics that make it a useful tool for application with inversion algorithms. Because of its many favorable attributes, RT theory has gained wide acceptance and usage in modeling electromagnetic scattering from vegetation. Radiative transfer theory was initially developed by physicists for the analysis of electromagnetic wave propagation in the interstellar medium and planetary atmospheres [9]. Because of its usefulness in this capacity, RT was subsequently applied to other types of particulate random media such as vegetation cover. When the medium consists of sparsely distributed scatterers that are small in comparison with the wavelength of the incident radiation, the underlying assumptions on which RT theory is based are valid, and the model can be expected to perform reasonably well. However, in constructing RT models for forests and stalk dominated agricultural canopies some basic conditions necessary to the validity of the method have been overlooked. The RT approach is

formulated in terms of the single scattering properties of the particles in the medium; i.e., it is assumed that particles are in the far-field of each other and are illuminated, locally, by plane waves. A tree canopy usually contains particles, such as trunks and branches, which are much larger in dimension than a wavelength. In addition, the forest medium consists of strong scatterers that may be quite densely distributed. In this case both the far-field and local plane wave conditions are violated. The purpose of the modeling chapters of this thesis is to investigate how the limitations imposed on RT theory by its underlying assumptions affect the result of its application in the analysis of radar scattering by vegetation canopies. In the initial stages of this work it was discovered that the application of first-order RT theory to the modeling of dense, stalk dominated agricultural canopies, such as corn, produced results that were in sharp disagreement with the body of previously measured data. Much of the previously acquired experimental data was obtained using the earlier generation of radar scatterometers, and it was felt that with better equipment and improved calibration techniques it would be possible to resolve the modeling problems. An investigation of radar scattering from mature cornfields was then initiated that used the improved equipment and methodology coupled with extensive and careful collection of information about the physical characteristics of the canopy. It was found that the newly acquired data confirmed the existence of a discrepancy between the predictions of first-order RT theory and measurements. As a result of these investigations it was felt that perhaps higher orders of iteration of the radiative transfer equations, corresponding to higher levels of multiple scattering, would provide the solution to modeling the radar return from this type of dense canopy. Chapter 2 of this thesis presents the development of the second-order RT model for the trunk layer of a general vegetation canopy and compares the results

5 of application of the first and second-order models with radar measurements of the full corn canopy. This chapter shows that while second-order RT provides a dramatic overall improvement in the estimate of the backscattering coefficient for both the co-polarized and cross-polarized radar returns from this type of vegetation, it does not reproduce the measured angular trend for either the vv or cross-polarized canopy response. The conclusion drawn from this part of the work is that while multiple scattering is indeed important for dense canopies of this sort, RT theory cannot correct the modeling problem by higher levels of iteration since it cannot account for the true nature of the local illuminating fields inside the canopy. Multiple scattering in random media consisting of large scatterers causes non-uniformity of illumination of the scatterers with respect to amplitude, phase and polarization of the local excitation. It is apparent that the true phase and extinction matrices for the constituent particles in the medium, which control the angular response of the radar return, are no longer accurately represented by their simple single scattering characteristics based on the assumption of local plane-wave illumination. In particular, the correlation distance for local fields inside the medium may be smaller than the physical size of the scatterers which causes overestimation of any quantity, such as the extinction matrix, that depends on the phase matrix of these particles. This situation provides the motivation for Chapter 3 of this work, which develops a Monte Carlo scattering model for the trunk layer of a forest canopy. The model accounts for near-field interactions between the cylinders composing the canopy trunk layer up to second-order and provides a benchmark against which to evaluate the performance of RT theory in a regime for which no other analytical models are satisfactory. The interaction terms between pairs of cylinders are validated with measured data. Monte Carlo simulations based on this model are presented for various repre

sentative cylinder number densities, and the simulated results are compared with radar measurements made on prepared ensembles of randomly positioned vertical rods above a ground plane. First and second-order RT solutions for the same media are also compared with the simulated data. This work verifies that the RT model gives incorrect results for media in which there is a preponderance of scatterers that have dimensions that are large compared to the excitation wavelength. The natural extension of this investigation is to canopies having more than one type of particle. The importance of the heterogeneous canopy problem is dependent on the nature and significance of smaller particles distributed either above or within the canopy volume occupied by the large vertical trunks. What is the effect of smaller scatterers such as leaves and twigs on the radar response of the canopy as a whole? Do these components couple significantly with the trunks, and in what ways does their presence alter the performance of the RT model for similar canopies? Chapter 4 of this thesis details a general technique based on the reciprocity theorem for deriving the secondary scattered field from a pair of adjacent objects. The general formulation is then applied to obtain analytical expressions for the secondary scattered field from a cylinder-sphere pair. The expressions are validated using the method of moments. This development provides the basis for the construction of Monte Carlo simulations for a heterogeneous canopy structure consisting of vertical cylinders, representing trunks, and smaller spheres, representing components of the canopy crown. This is the subject of Chapter 5 which applies the formulation of the previous chapter to the case of spheres and cylinders distributed above a ground plane, and which takes into account mutual coupling between cylinders and coupling between spheres and cylinders in the ensemble. The results of these simulations are compared with those obtained from RT theory for the same canopy configurations.

7 It is shown that the radiative transfer models do not properly predict the scattering behavior of media composed of densely distributed scatterers having dimensions large compared with the excitation wavelength. This is partly a result of the fact that the extinction matrix for this type of medium as computed by RT theory is overestimated, especially at angles of incidence far away from vertical. It is also partly due to the incorrect treatment of the RT source function for volume scattering in media consisting of particles having dimensions that are not small compared to a bounding dimension of the medium. Chapter 5 also shows that RT models that attempt to divide the canopy into an upper layer consisting of trunks and smaller particles and a lower layer consisting of trunks alone can produce results that are seriously in error. One further purpose of this chapter is to demonstrate that while second-order interactions between cylinders in the trunk layer can cause significant changes in the level of backscattered co-polarized radiation, the predominant effect of the interaction of trunks with smaller particles is to produce a change in the angular trend of the cross-polarized canopy radar response. Second-order interactions between spheres and cylinders in the heterogeneous canopy structure are also shown to result in distinct changes in the co-polarized phase statistics of the backscattered wave. In Chapter 6 of this work a hybrid model is presented for computing radar scattering from layered vegetation media consisting of vertical trunks and small, weakly scattering particles above a dielectric ground plane. The crown layer is modeled using radiative transfer theory, and the trunk layer is simulated using the Monte Carlo method. The analytical derivation of the transmissivity matrix for the trunk layer is presented, and some of its important features are investigated. The concept of the effective scattering matrix for an average cylinder in the trunk medium is developed

and evidence is presented that it assumes limiting behavior as the effect of multiple scattering in the medium becomes increasingly important. Evidence is also presented to confirm that the exponential extinction model used in RT theory works fairly well for sparse distributions of small particles and even for sparse distributions of extended scatterers. However, it is shown that the RT extinction model breaks down for high densities of strongly scattering particles that are large compared with the excitation wavelength. It is demonstrated that the source function integration used in RT theory to account for volume scattering leads to results that are severely in error for layers of long cylinders, and it is concluded that this is true for distributions of other large scatterers. Finally, several types of test canopies consisting of dielectric cylinders and small metallic spheres are investigated. It is shown that the hybrid model developed in this chapter can be an effective way to compute scattering from vegetation media having this general structure.

CHAPTER II THE SECOND-ORDER RADIATIVE TRANSFER MODEL OF THE CANOPY TRUNK LAYER 2.1 Introduction Because most canopy cover is at least semi-random in character and because much of it is statistically homogeneous and area extensive, the radiative transfer (RT) method has been widely applied in modeling electromagnetic scattering from vegetation. In formulating such models the most important features involve scattering by canopy constituents such as leaves, branches, trunks and the ground. Inherent in the radiative transfer approach is the ability, by iterative solution of the equations, to elucidate the individual scattering interactions between the electromagnetic wave and the canopy constituents. Such interactions may consist of direct scattering of the incident wave by the individual leaves and branches in the crown layer or by the underlying rough ground surface. In addition, depending on the degree of iteration, RT theory accounts for mutual coupling between the canopy constituents by using the single scattering properties of the constituent particles and the local plane-wave approximation. The electromagnetic scattering behavior of the canopy constituents is clearly a

function of their various geometries and electrical properties. In this regard, one would like to use relatively simple approximations to the constituent geometries in order to keep the model tractable. Cylinders may be used to model trunks or stalks and branches. Circular or elliptical discs, with or without curvature, may be used in modeling leaves. Simple rough surface scattering models such as the Kirchoff model or the small perturbation model may be appropriate for the ground layer or, alternatively, an empirical may be used. Various models have been developed based on a first-order solution of the radiative transfer equation and simple canopy geometries [16, 42, 45, 96]. Such models have had some degree of success in predicting the co-polarized radar return from canopies with relatively low densities of strong scatterers or in cases where there is limited penetration by the wave into the medium (as exists, for example, in leaf dominated canopies at high frequencies) [53]. In many circumstances, as for instance in agricultural and dense forest canopies, the low-frequency microwave behavior is dominated by strongly scattering stalks packed in high density. In a mature corn canopy, the wavelength of L-band microwave radiation in the dielectric stalk medium is close to the diameter of the stalk. Under these conditions resonant scattering from the stalks is exceptionally strong; and it has been found that the leaves do not seem to play an important role in determining the co-polarized backscatter response of the canopy [94, 104]. Measurements made by this investigator [61], do not compare well with the firstorder RT model of corn canopies at low microwave frequencies. In addition, the first-order RT backscatter model for a canopy consisting of primarily vertical trunks does not generate any cross-polarized return at all, whereas experimental data on scattering from corn shows that there is a significant level of depolarization that

may not be accounted for by the natural orientational diversity of the stalks. It is also observed that the first-order RT model severely underestimates the co-polarized radar return from trunk dominated canopies and gives an incorrect angular trend for the vv-polarized backscattering coefficient, falling off much more sharply at angles of incidence away from vertical than is seen in measured data. These observations have provided the motivation for an investigation of the second-order RT scattering mechanisms in trunk dominated vegetation which is the subject of this chapter. When plant stalk heights are large compared with the wavelength of the incident radiation it can be shown for vertical trunks that scattering occurs primarily within a narrow cone centered on the specular direction as illustrated in Figure 2.1. Within this specular cone there may be other stalks present as well as the ground itself. If reflections from the ground are considered to be mainly specular in nature, the cone of scattered radiation from a particular stalk when inverted on reflection from the ground surface may intersect other stalks depending on the density of the medium. The firstorder model consists then of wave interactions between a single stalk and the ground as well as diffuse backscatter from the rough ground surface. These terms are referred to as the ground-trunk, trunk-ground and direct ground interactions and are illustrated in Figure 2.2. The second-order model includes an additional wave-stalk interaction between pairs of stalks in the presence of a specular ground surface. The terms generated by the second-order model are referred to as trunk-trunk-ground, groundtrunk-trunk and trunk-ground-trunk interactions. Another kind of second-order effect involves the interaction between single stalks and diffuse bistatic scatter from the rough ground surface. These terms are called diffuse ground-trunk and diffuse trunkground interactions. All the second-order terms are illustrated in Figure 2.3. In the remainder of this chapter the radiative transfer model will be derived for

kmc Fi 2 S i t u Figure 2.1: Scattering into the Specular Cone a layer of vertical trunks above a rough ground surface up to, and including, the second-order terms. The first and second-order RT models will then be compared with measured data to determine whether second-order RT theory provides a satisfactory solution to the problem of reproducing the electromagnetic backscattering behavior of trunk dominated vegetation canopies.

13 Trunk-Ground Interaction Ground-Trunk Interaction Direct-Ground Interaction Trunk Specular Ground (Rough Ground Specular Ground / Surface Surface Surface Figure 2.2: First-Order Effects 2.2 The Radiative Transfer Model 2.2.1 Basic Equations and Definitions We now consider the salient features of a general closed canopy. By closed canopy it is meant a canopy that is continuous and statistically homogeneous in the horizontal plane but that has significant variation from top to bottom. The top of the canopy consists of a crown layer made up of leaves and branches. The leaves are usually modeled as flat dielectric discs and the branches as dielectric cylinders. The leaves and branches in the crown layer are described in terms of a number density (number per unit volume). The dimensions and orientations of these components are

14 Diffuse Ground-Trunk Trunk-Ground Rough Ground Surface Trunk-Trunk-Ground Trunk-Ground-Trunk Ground-Trunk-Trunk Interaction Interaction Interaction Smooth Ground Surface Figure 2.3: Second-Order Effects usually specified either as fixed values or in terms of probability density functions. Underneath the crown layer is the trunk layer which is composed of dielectric cylinders. Again these components are described in terms of the number density (number per unit area), the diameter, the height and the orientation. The lower level of the canopy is a rough dielectric surface that is used to represent the ground. It is usually characterized by rms height and correlation length of the surface roughness scale. One of the rough surface scattering models such as Kirchoff's scalar approximation

TCrown Layer\ Trunk Layer,-Ground Layer —Figure 2.4: Closed Canopy Model to physical optics, the small perturbation model, or, in the case of very large-scale roughness, geometrical optics is often used depending on the roughness scale of the ground surface under consideration [92, 95]. The general canopy geometry is shown in Figure 2.4. In a crown region made up of scatterers with small albedo the first-order radiative transfer solution may be a sufficient approximation to the canopy behavior. However, in a trunk region composed of strongly scattering stalks at high number densities the first-order approximation may be a poor representation of reality and higherorder terms should be included in the model. In this section the vector radiative including the second-order terms.

For an incident plane-wave, the electric field vector may be decomposed into a vertical and horizontal linear polarization basis: E = (Eivvi + Ehhi) ik okr (2.1) where ko is the free-space wavenumber, and vi = cos Oicos qsix + cos Oisin qi4 - sin Oiz (2.2) hi = -sin id + cos (2.3) ki = sin Oicos qix + sin Oisin qY - cos i (2.4) In (2.1) a time dependence of the form e-iwt is assumed and suppressed. The vector specific intensity (W nm-2 sr-'1 Hz-1) for the incident coherent wave is defined through the modified Stokes parameters I,, Ih, U and V as follows: IV, IEEIE Io - 1 (2.5) Ui? o 2Re (Ev E.*) V' 2S2m (Ev Ef*) where T0 is the intrinsic impedance of free-space. The scattered wave from the distributed target is of spherical character and is partially coherent. If A is the illuminated target area and 08 is the angle between the outward normal to A and the scattered wave-vector, then the scattered intensity must be normalized by the solid angle A, where r is the distance from the target to the observation point. Thus the intensity of the scattered wave may be written in

17 terms of the Stokes parameters as: VI (lEVI 2) tIs =_____r2 (IE ) (2.6) (2.6) uswA cos 0 2Re ((Es E1*)) 2Qjm ((Ev El,*)) and ( ) denotes the ensemble average. The intensity of the plane wave incident on the upper diffuse boundary may be written as: Ii Io6(cos 0 - cos 0o)6(q - o0) (2.7) Inside the trunk (stalk) layer the intensity It is separated into two components, It(y, q, z) and It (I, 4, z) which represent the upward-going and downward-going intensities respectively. The subscript t denotes the trunk layer, and P = cos 0 where 0 is defined with respect to the positive z-axis. The geometry of the trunk region problem is shown in Figure 2.5. The radiative transfer equations inside the trunk layer are: 2i9 Ot I- (,,z) + t(2.8) ~ —Rt(~ t Z) t(2.8) -dIt(-p,, tz) = -I'(-p, q, z)+ Ft(-, q, z) (2.9) where we have defined: Ft(p,, z) P= jJ Q;t(t',,') It (I', I', z)ddQ' + Pt(y, 0;-p', O') It(-,l',' z)dQ' (2.10) and pFj(-p qi, z) = jo JPt(-1 ~;', sb') It (,u, q', z)d'

z-axis A A:0:0'y Io( go Free Space Diffuse Boundary z=-d Trunk Layer!o/ +'S It(- O' Rough Ground Surface z=d Ground Layer E g Figure 2.5: Trunk Region Problem Geometry dQ' is the differential element of solid angle, df' = -dp'dq' = sin O'dO'dq'. The quantity Pt (u, 0; yl', q') is the phase-matrix of the trunk layer which accounts for scattering by the trunks of radiation incident from the direction (#', q') into the direction (1t, q). In (2.8) and (2.9), Act is the extinction matrix for up-going or down-going radiation in the trunk layer. Therefore, the first term in the radiative transfer equations accounts for extinction of radiation as it travels through the medium, while the source terms Ft~ account for scattering by the trunks of radiation from all directions at depth z into (p, q). It has been assumed that the trunk layer is isotropic with

respect to the azimuthal angle q. The solutions to the coupled radiative transfer equations are: It+([, X, z) - e-t(z+d)/"It(, (, X,-d) + L e-t(ZZ')/'Ft (l,, z') dz' (2.12) and It(-i, /, z) = e-7tZ/iI(-I, _, 0) + j e-t(Z-Z')/"Fi(-tu,, z') dz' (2.13) For the diffuse boundary at the air-trunk interface we have: It-(-, 4, 0) = - IoS(t- _o)5(5 - ) (2.14) while at the lower solid boundary with the ground we write: r2-r rl It' (,, -d) = jJ(,, q; -', O')I7(-/s', 4',-d)dp'dO' (2.15) where g is the scattering phase-matrix for the ground surface. Substituting these into (2.12) and (2.13) we have: It(1 of z) = e4K(z+d)/Il Ljrg1(,, 4; -pl.X (-,., -d)dp'do' +-ie- e t(z-z')/lFt (I,, z') dz' (2.16) d and It(-L, X, z) = e-z/iIo6(, - ~o)b(o - q) + j e _(zz_)/FT( z') dz (2.17) For the stalk dominated field under consideration we will use the scattering matrix based on a modification of the matrix for infinite vertical cylinders. This is a good approximation for cases of practical interest since the length of stalks for most of

20 the growing season is equal to many wavelengths, even at L-band. The solution of the radiative transfer equations is constructed in an iterative fashion. The zerothorder solution is obtained by neglecting the source terms. The zeroth-order source terms are then formed from the zeroth-order solutions. The first-order solution may then be constructed using these source terms. Proceeding in this way solutions can be generated to as high an order as is desired. As a side benefit of this process, the solution is cast in a form that elucidates the individual scattering mechanisms present in the medium. It is now demonstrated how solutions up to second-order are generated. 2.2.2 Zeroth-Order Solution For the zeroth order solution scattering by the trunks is ignored. Setting the phase-matrix for the trunks equal to zero, we effectively turn off the source terms. Then: p,(, ~;',l')- =0 0 = Ft(u,q,Z)-=0, V (u,J',z,q',z) Thus: It(~)-(-i, q, z) = e tI'/zIoS(/l-lu o)6(~-uo) (2.18) and Ito)+(y, q, z) = e-4tg(z+d)//g(,, c;g-o, Co)e-It d/~OIo (2.19) from which, following substitution, the zeroth-order source functions are found to be: Fr?)+(1, I, z)-117 o'(. O;,' )e- (z+d)/,( t,. O';-t, O)e-'td/MoI~dQ'

21 +P-1,;H X~(2.21) It is evident from (2.18) and (2.19) that the zeroth-order solution contains the direct-ground term alone. The first-order solutions can now be derived. 2.2.3 First-Order Solution Substituting these source terms back into the original equations we find: F'-(_, X, z) eKt-ZI/ 6(L - )6( - 0) e-+ e (Z )/"F~-(-,, z') dz' (2.22) t - ( - p% d o (2.22)!~TM (p, ~, z) = e-l(z+d)/] ] 5(p, (;1-',,;.)I.(-1(', q,-d)d/u'dqY + I e- Ct(z-z')l)/F(~)+()p,, z') dz' (2.23) d To compute the first-order solution explicitly, several assumptions are made. We simplify Fo)' + and Ft() - using the fact that the phase-function for the vertical cylinders used to represent the trunks is such that there is no scatter from them back into the hemisphere of incidence. Thus:'Pt (-A, 1; P', S') ='Pt(p, q; -ii', q') = 0 (2.24) Therefore there is also no direct backscatter from the stalks themselves, except at normal incidence. Since it has been assumed that the stalks are much longer than a wavelength, radiation scattered by the stalks is scattered into the forward cone with half-angle

22 equal to the specular angle as was shown in Figure 2.1. This is a property of the infinite cylinder scattering matrix. Then:'Pt(p1, 6; i"', q') = P7::t(,'P; iA',') &k(I - 1u') (2.25) where Sk(# - Y') is the Kroneker delta. After utilizing these relations the source functions become: Ft()-(-, I, z)= -Pt (-#i, P; - o, Oo)e<tz/'~Io (2.26) and Ft()+( -1,, Z) = - Pt(',; i',' (z+d)/(1 "I'; o )e-d/I'~Iodq' (2.27) Then we can write for the region above the stalk layer: It(')+(o,,0o + r, 0) = Ig + It + I t (2.28) where the superscript "b" means that the term is evaluated in the backscatter direction. The three terms correspond to direct-ground, trunk-ground and ground-trunk interactions as shown in Figure 2.2 and are given by: I = e-tl'/~o(oa0, ro + 7r; -lo, qo)e-id/l'~Io (2.29) I = -e 1(/o, ~0 + 7r;-o, -). JIt jPoe/o e-dt (- _/,Io; -, o Oo)e~?z'/l odz' dr'Io (2.30) d~.t. txt..

23 The extinction matrix is given by e-CQ/t' = Qt(p, ) Dt(,,; -z/t) Q-'(#, q) where Qt(p, q) is a matrix whose columns are eigenvectors of the extinction matrix Ct(/t, q) and E)t(u, b; -z/u) is a diagonal matrix with elements: [EVt(l, q; -zl/)]ii =- e (")z/' (2.32) with Ai(/p, q) the ith eigenvalue of ICt(tt, 0). For trunks we have t4 = t(t, b) = Kt(-p, 0) = ai due to vertical symmetry of the trunks and specular forward scattering in the trunk layer. We can then write: Ibg =-_et/o (o, 0 + 7r; -Po, ) A1(;-(-,o,; — o, q') Iod1' (2.33) where: A1(-1Po, q;-Po, 0o) = Qt(-jo, 0) Ai(-o, q;- -o, 0o) Q;-'(-po, Oo) (2.34) and [ (_ k *; _ ~ =)~ l [eAi(-o,,,)d/po _ e-A(-juo, d)d/ ] [- o(-#o, q) + Aj(-po, 0o) *[Q.7(po, ) )Pt(-/o,0;-Po, o) Qt(-po, o)]j (2.35) Similarly: - I fI 2r gt= - I 2(#o, o + i; oi, q') 9(po, q'; -/o Oo) do' e-td/uo" Io (2.36) where: A2(Puo, qo + 7r; Po, q) = Qt(Po, qo + 7r).A~(/Po, ko + 7r; io, I4) Q 1 (o, I ) (2.37) and ~o [e-Aji("o~~o+r)d/"o _ e-Aj(~o,')d/~o] [A4 (po, Oo + 7r; Po, qI)],= -xi(io-, 0 -r) + Xj(,o, 4) *[ -A' (go, o + 7r) Pt(o, +o + r; io, d) Q,(io, )] i (2.38) These are the first-order terms.

24 2.2.4 Second-Order Solution The procedure for derivation of the second-order terms is the same as for the first-order terms. The first-order source functions are obtained from the first-order solutions and these source functions are used to compute the second-order solutions. The second-order terms consist of the single stalk interactions with diffuse reflectance from the ground and double stalk interactions with specular reflectance from the ground. The terms are: diffuse ground-trunk, diffuse trunk-ground, trunk-trunkground, trunk-ground-trunk and ground-trunk-trunk as shown in Figure 2.3. For the diffuse ground-trunk term we find: gt2 = ~t et'/ P(o,,(, o + 7r; o, ). ent(/'+d)/,oa(lo, ('1; —Po, 0o)d4'dz'e- td/,oI o (2.39) where Qd is the diffuse-ground phase matrix. Using the approach of the previous section we have: gt2 = 1 A3(joi, o + ir; o,') d(,o,'; -/-1o, $o) d' - "etd/' Io 1 (2.40) where A.3(0o, q$O + O r; 0 o, >) = Qt(oo, 6o + rr) (/uoo, 6o + 7r; +ro;, ) Qt (to, 6) (2.41) For vertical trunks Qt is independent of (0, q), so we will write simply Qt from now on without including the angular dependence. Then: vA(Lo, qo + 7r; Po, ) = Ot(lo, o0 + ir; z'/Ho) ~Qt-l t(lo, Co + A; pox, ) Qt ont(Lo, 4'; -(z' + d)/lo) dz' (2.42)

25 Since q0 is arbitrary it may be set equal to zero. Upon performing the z'-integration above we find: /~ [e-A'(4"o~)d/to _ e-Aj(o,~')d/o] [w(Po, r; Po, )]j = -(o ) + =(o ) * [Q;- pt(po,; o, qB) Qt] (2.43) Similarly, for the diffuse trunk-ground term we obtain: 1 _Ktd/.o/2i Itg2 = -e d/ d(o, r;-P-o, ) Qt A4 A(-po,'; -po, 0) Q-' d' I, (2.44) lPo Jo with: o [e-"'(-~'~')d/~ _ e-j(-o~O)d/"~] [.A4(-lo, /);-,o, O)],ij = -x.:(-, o ( ) + Ai-g,, 0) * [Q-1'Pt(-Yog, q0;-/o, 0) Qt] i (2.45) The trunk-ground-trunk term is computed as follows:.el(t(z'+d)/opt(_o,; -o, o)eKtaz codz/dzd) Io = i Qt5(o); -=Qo,01 Qc) dtIo (2.46) where 5c is the specular (coherent) ground reflectivity matrix [84]. We also have: /J_.d I (yto, 7, Po i; z") 7"(Po).2A(-po, 0,-po, 0; z') dz'dz" (2.47) with 7E"(po)- Qt; gc(po) Qt (2.48)

and 41 (yLo, 7r,o,; Z") = t(pto, Ir; Zl"/to) S1 (/o, 7r; o,,'1) )Vt(pol,'; -(Z" + d)lapo) (2.49) 4M2(-Lo, q, -Po, 0; z') = Vt(-po, +;-(z' + d)lpo) S2(-Po, +;-po, O) VD(-Lo, O; z'lpo) (2.50) The definitions of Si and S2 are: Sl (Po, 7r; Po, O) = Qt 1 Pt(Po, 7r; Po, O) Qt S2(-Po,; -Po, 0) -= Qt1 Pt(l-o, O; -Po, 0) Qt We now evaluate the kernel of (2.47): [Ai,1'1"?A2]ij = E (MI)ik (7Z")kl (A42)/j (2.51) k'l It may be easily verified that'R" is diagonal. So we can write: (-R")kl= rkl lk1 = [ Qt-1 c(Po) Qt]kl b5kl where 6kl is the Kroneker delta. Thus: [A41 RI -t 2] ij - (l41)ik rkk (M42)kj (2.52) k Using the definitions (2.49) and (2.50): (.A4I )ik - eAiz/o e-k(z" +d)/po (S )k (2.53) and (.A2)kj - eAjZ'/~oe-Xk(z'+d)/MO(S2)kj (2.54)

27 where the directionality (p, q) of the eigenvalues A has been assumed and suppressed for the sake of notational convenience. Combining (2.54) and (2.53) in equation (2.52) gives: [A1 Z" A42]ij _= > e-2Akd/lo e(Ai-Ak) z"/o e(Aj-Ak) z'/l0 ( Sl)ik (S2)kj rkk (2.55) k Producing the final result: (45)ij = E e-2kd/~ (Sl)ik (S2)kj rkk e(AiAk)z"/h e(Aj-Ak)z'O dz'dz" ~(e-id/ _- Aid/o~) A/(e-Akd/o~ - e-Ajd/o) (Sl)ik (S2)kj rkk k (A - Ak) (Aj - (.6 (2.56) which is substituted back into equation (2.46) to give the trunk-ground-trunk term. The trunk-trunk-ground term is derived as follows: 1 o2'r [~ tdluoo (z"d p I t91= -- l -JtdJLIe og (,go)e-t(-z + )/ ( -,oo,,)e t(,-z,)/r; o -Pt(-Io, 0'; -Po, O)eCtZ'l"odz' dz" dd' Io 1 _Ktd/M() (27! e_;td/.ogl(/Zo). Qt, 46 Q dO' Io (2.57) ]o o (2.5) After rewriting the exponentials as the similarity transformation of the eigenmatrix At we have: 46 = (- 7r; -(z" + d)/a) A (- u- 0 $'; z') D.(-, )0; z'/lto) dz' dz" (2.58)

28 where: S4(-o., r; -Po, phi)Dt (-po, phi; (z" - z')/Po) S3(-Po, 0; -no, O) (2.59) and S3(-Po, 4; -Io, 0) = Q7t' pt(-1o,; -Po, O) Qt S4(-Po, X; -/o, 4) = Qi''Pt(-Poi, r; -o, q) Qt In addition, since Vt is diagonal: [Vt A4 Vt]ij = e-Xi(z"+d)/O eAjZI'/l (A)ij (2.60) Using these results in equation (2.58) we find: (A4) = J~ e-X (z"+d)/o eAjZ'/Io [(-o,'; zz")]ij dz' dz" (2.61) Expressing (2.59) in element form: (M)ij = E (S4)ik (Vt)kl (S3)/j (2.62) kl but (Dt)kl = (Vt)kk 6kl therefore: (-)ij = (S4)ik (t)kk (S3)kj k = (S4)ik ($ 3)kj eAk (~"-Z')/"o (2.63) k We can now express (2.61) as the sum: (.A6)ij = Y (S4)ik ($3)kj (T)ijk (2.64) k

29 where (7T)ijk = j eAi(z"+d)/,o eAXjZ'//O eAk (z"-z')/,o dz' dz" d z - -- e-/ (2.65) (Aj - Ak) (A3 - Ak) (Ai - A) J(2) The elements of T may be evaluated in the degenerate case by using L'Hopital's rule. For the ground-trunk-trunk term, the same method is applied and yields: Itt =- (j2 Qt A7 Q-' d+') c(po) e-'d/~ Io (2.66) where.A7 is given by: A47= I d V t(Po, r; Zl/Po).A12(1o, ~); z'i Z"o) t(Po, 0; -(z' + d)/#o) dz' dz"(2.67) with A-4 22(Po, 0; z', z") = S6 V)t(o, phi; (z' - z")/lo) S5 (2.68) and s5 = Q-1'Pt((,o, 7r; Po,') Qt S6 = Q'Pt(L (o, (o;,o, 0) Qt Following a procedure exactly analogous to that for the trunk-trunk-ground term it is found after lengthy but straightforward manipulation: (c7)ij = (S6)ik (Ss)k, (T)ijk (2.69) k where (Y)ijk is the same as in equation (2.65).

The derivation of the extinction matrix; its eigenvectors and eigenvalues as well as the phase-matrix for the trunk layer is given in reference [95] and will not be repeated here. The first and second-order model results are now compared with experimental data taken on mature cornfields at the end of the growing season. At this time of year the dielectric constant of the leaves is much lower than that of the stalks, therefore the effect of scattering by leaves in the canopy may be ignored relative to scattering by the stalks. 2.3 Results and Comparison With Measurements Figures 2.6 through 2.14 show the comparison between radar scatterometer data taken at L-band frequencies on mature corn canopies and the first and second-order radiative transfer models. The comparison is shown as a function of incidence angle for several values of soil volumetric moisture (Mv) within the typical moisture range. It can be seen that, in general, there is a 5 to 10 dB discrepancy between the first-order model and the hh-polarized data. The first-order model also predicts no depolarization for vertical trunks. The difference between the first-order model and the vv-polarized data is, at best, around 10 dB at low angles of incidence and deteriorates rapidly with incidence angle. The agreement between first-order RT theory and these experimental results is quite poor. Inclusion of the second-order terms provides a marked improvement in the results for all polarizations, atleast as far as the mean signal level is concerned. The predicted hh-polarized return shows a 5-10 dB improvement in overall signal level as well as improvement in the angular trend of the data. The vv-polarized RT results are in much better general agreement with the measured data, however there are some

31 10.0 o 9/05/90 corn data 5.0 _ Total I st + 2nd Order......... Total I st Order 0.0 "0CP ~~~o -5.0 -10.0 -15.0. I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.6: Data/Model Comparison L-Band hh Polarization (Mv=0.05) 10.0 I 5.0 o 9/05/90 corn data (VH) o 9/05/90 corn data (HV) Total Second Order 0.0'0 -5.0 -10.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.7: Data/Model Comparison L-Band vh/hv Polarization (Mv=0.05)

32 10.0 o 9/05/90 corn data 5.0 _ Total l st + 2nd Order......... Total I st Order 0.0 Ott -5.0 -10.0 -15.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.8: Data/Model Comparison L-Band vv Polarization (Mv=0.05) 10.0 0 8/28/90 corn data 20.0 30.0 40.0 50.0 60.0 Figure 2.9: Data/Model Comparison L-Band hh Polarization (Mvde 0.1) Figure 2.9: Data/Model Comparison b-Band hh Polarization (My=O.16)

33 10.0 I 5.0 8/28/90 corn data (VH) 0 8/28/90 corn data (HV) Total Second Order 0.0 o -5.0 B -10.0 -15.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.10: Data/Model Comparison L-Band vh/hv Polarization (Mv=0.16)!0.0 o 8/28/90 corn data Total 1st + 2nd Order 5.0 5.0 - -- Total 1st Order 0.0 Ot) -5.0 -10.0' -15.0.. i. 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.11: Data/Model Comparison L-Band vv Polarization (Mv=0.16)

34 10.0 5.0 0.0 -5.0 o 9/07/90 corn data -10.0 Total I st + 2nd Order........ Total 1st Order -15.0 I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.12: Data/Model Comparison L-Band hh Polarization (Mv=0.32) 10.0 5.0 o 9/07/90 corn data (VH) 9/07/90 corn data (HV) Total Second Order 0.0 o -5.0 -10.0 -15.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.13: Data/Model Comparison L-Band vh/hv Polarization (Mv=0.32)

35 10.0 5.0 0.0 -01 -5.0 o 9/07/90 corn data -10.0 Total Ist + 2nd Order......... Total 1st Order -15.0. 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 2.14: Data/Model Comparison L-Band vv Polarization (Mv=0.32) serious problems with the angular behavior of the model. The model also shows resonant behavior for angles of incidence between 200 and 400 which is completely absent from the measured data. This resonance behavior is actually a result of two competing processes occurring within the RT model for the stalk medium. As the angle of incidence increases away from the vertical direction, the vv-polarized phase function of the vertical cylinders also increases rapidly. At the same time, the diagonal extinction matrix element for the vv-polarized intensity is also becoming larger. Thus, while the overall level of scattering in the medium is increasing with incidence angle, the transmissivity of the cylinder layer is decreasing exponentially. Above 40~ the exponential decrease in the vv-polarized transmissivity is the dominant stalk related effect to be observed in the RT model response. Also, coupling between the incident wave and a single cylinder is a sensitive function of the cylinder diameter especially

within the resonant scattering regime. This effect is not observed in measured radar data because there is a random distribution of stalk diameters in real cornfields that is not accounted for in this RT model. Another problem with the vv-polarized response of the RT model is the presence of the Brewster angle effect in the rough ground phase matrix. This effect seems to be completely absent from the scatterometer data for agricultural fields and is also conspicuously absent in the SAR data for all types of canopy covered terrain. The model depends to a large extent on specular forward scatter from the ground which provides the basis for the many ground-stalk type interactions. This specular forward scatter is always coupled to the reflection coefficient of the ground surface which produces a null in the model response at the Brewster angle. If physical optics is used to compute the diffuse reflectance from a rough ground surface, it is also found that the reflection coefficient dependent scattering terms disappear at the Brewster angle. Actual experimental vv-polarized data show a gentle decreasing trend with angle of incidence and is almost identical with the trend seen in the hh-polarized data. The vv data are, however, a few dB lower in magnitude. This means that the effective reflection coefficient for both polarizations is approximately the same in a real canopy ground layer. Actual ground is composed of layers of dirt, stones, twigs and other debris compacted in regions of varying density separated perhaps by pockets of air and moisture in the volume just below the surface. Obviously, different soils make ground layers with differing macroscopic properties. The point is, as far as electromagnetic scattering is concerned, it is hard to imagine that real ground may be characterized in terms of a simple, homogeneous medium with a well defined surface layer. It is more reasonable to expect that volume scattering in the soil sub-surface region may play a

significant role in scattering from the canopy ground layer; especially approaching the Brewster angle where transmission into the underlying medium is of greater importance than reflection from the surface. It is consistent with this view to suppose that the true scattering mechanism in real-world soil media is a mixture of surface and volume scattering. The relative predominance of one mechanism over another would be a function of polarization, frequency, angle of incidence, soil volumetric moisture and the complexities of the ground structure itself in as far as both its surface and sub-surface properties are concerned. The investigation of this aspect of the problem is, however, outside the scope of this work and will not be treated in this thesis. One important feature of the angular response exhibited by the RT model may be seen in a comparison of Figures 2.8 and 2.14 for the vv-polarized backscattering coefficient. In Figure 2.8 the volumetric moisture of the soil surface is 0.05. At this moisture level, the model of El-Rayes and Ulaby [91] predicts a value of = 3+i0 for the relative dielectric constant of this soil type. For this dielectric constant, the Brewster angle is at 60 and a null in the predicted vv-polarized backscattering coefficient is quite apparent in the Figure. However, Figure 2.14 illustrates the case for a soil volumetric moisture value of 0.32. At this moisture level, the Brewster angle would be at around 80~ for a specular surface. In addition, the soil now has a loss factor near to unity which prevents the sharp null that would otherwise be present at the Brewster angle. In this case the predicted angular response for the backscattering coefficient is being produced almost entirely by the RT extinction model. The measured data shows a decrease of approximately 5dB over the angular range from 200 to 600 while the computed RT data falls off by over 10dB within the same range. It is evident that even though we have ignored the effect of leaves in the RT model which would, if anything, be expected to increase the attenuation of the

incident wave, the RT model severely overestimates the extinction of the vv-polarized wave by the canopy. Examination of Figures 2.7, 2.10 and 2.13 shows that, with inclusion of secondorder effects, the overall level of depolarization predicted by the RT model is commensurate with the level indicated in the measured data. The model, however, exhibits a much more pronounced decrease in the backscattering coefficient with angle of incidence than does the data. This decrease is on the order of 10dB over the entire angular range and is present for all values of the soil moisture. The measured data also shows a decrease in crosspol level over the angular range, but the magnitude of this variation is only about 3dB for all soil moisture conditions. It is likely that there would be less extinction of the predicted canopy crosspol response if orders of scattering higher than 2 were included in the model. This is because the amount of depolarization produced by the canopy is directly related to the level of multiple scattering present within it. As the angle of incidence increases, the level of multiple scattering should also increase. However, it is by no means certain that the RT solution will converge as the iteration level increases without bound since RT is an incoherent approach. 2.4 Conclusions In this chapter the second-order radiative transfer model for a layer of vertical trunks above a dielectric ground surface has been presented. The first and secondorder models have been compared with experimental measurements made at L-band frequencies on mature corn canopies. It has been shown that second-order RT theory provides an overall improvement in the predicted backscattering coefficient for both the co-polarized and cross-polarized radar returns. It has also been demonstrated that

RT theory does not reproduce the correct angular trend for the vv-polarized or crosspolarized canopy response. This may be attributed in part to the Brewster angle effect which is absent from the measured data, but it is also partly a result of overestimation of the extinction matrix by the RT canopy model for the vv-polarized coherent wave. In the case of the cross-polarized backscattering coefficient, it is believed that the inclusion of higher order scattering terms might improve the angular response of the model, but there is no guarantee that this is correct because the addition of more terms would increase the mean signal level that is already too high at the lower angles of incidence.

CHAPTER III MONTE CARLO SIMULATION OF SCATTERING FROM A LAYER OF VERTICAL CYLINDERS 3.1 Introduction In applying radiative transfer theory to the modeling of a forest medium, some basic conditions necessary to the validity of the method have been overlooked. This model is based on the single scattering properties of the particles in the medium; i.e., it is assumed that particles are in the far-field of each other and are illuminated, locally, by plane waves. A tree canopy usually contains particles, such as trunks and branches, which are much larger in dimension than a wavelength; therefore, the farfield condition is not satisfied. Moreover, since these large particles are embedded in a random medium, the magnitude and phase of the field distribution illuminating the particles are non-uniform; thus, the plane wave illumination condition is violated. It is the purpose of this chapter to demonstrate the shortcomings of the radiative transfer technique for a medium containing particles that are large compared to the wavelength of radiation in the medium. The trunk layer of a forest canopy consisting of vertical dielectric cylinders over a dielectric surface is considered. A Monte Carlo simulation of the scattering problem, which includes multiple scattering up to second

order, is developed to provide a realistic solution for comparison. The second-order scattering term in this solution allows the long cylinders to be in the near-field of each other. The validity of the Monte Carlo simulation is verified with experimental backscatter data collected from a random distribution of metallic cylinders over a perfectly conducting ground plane. 3.2 Scattering from Two Adjacent Cylinders at Oblique Incidence Exact analytical solutions to the electromagnetic scattering problem exist for only a very limited number of geometries, including infinite cylinders. It has been shown that for cylinders that are very long relative to the wavelength, an approximate solution can be obtained based on the solution for the infinite length case provided d/L < 1, where d and L are the cylinder diameter and length, respectively [100, 41]. Although more accurate solutions for the finite-length cylinder can be obtained using numerical techniques, the solutions obtained in this way are not desirable since they become very inefficient when the dimensions of the cylinder are large compared to the wavelength. Similarly, for two finite-length cylinders adjacent to one another, an exact solution does not exist and numerical solutions are even more inefficient. In this chapter we resort to an approximate iterative scattering solution. We assume the cylinders are much longer than the excitation wavelength, and that they are mutually in the near-field with respect to their longitudinal dimensions but are in the far-field with respect to their diameters. The approach taken is to find the scattered field from the first cylinder as an isolated body given a primary plane-wave excitation. The response of a second cylinder to the cylindrical-wave excitation from the first cylinder is then found. The effect of the second cylinder on the first is obtained by reciprocity. In this way an ap

proximate analytical solution may be obtained which accounts for multiple scattering to second-order. In principle, this procedure can be continued to any desired order of approximation; however, a price is paid in terms of the complexity of the final solution. By invoking the field equivalence principle, the cylinders can be replaced by equivalent electric and magnetic surface current densities given by Ja(, z) = n x H (3.1) and Jm(c/) z) = -i x E (3.2) where E and H are the total electric and magnetic fields on the cylinder surface and n is the outward surface normal. If the surface currents are known, the scattered electric field may be obtained from E8(r) = V x V x Ie(r) + ikoZoV x Im(r) (3.3) where He and Im are the electric and magnetic Hertz vector potentials respectively. The electric Hertz vector potential is given by = iZo eiko lr-r' I II = i J(r') dS' (3.4) 4irkoI r - r I In this chapter the e-iwt time convention has been assumed and will be suppressed throughout. The magnetic Hertz vector potential has a similar form with Zo replaced with Zo1 and Je replaced with Jm. The surface currents on a long cylinder are approximated by the surface currents on a corresponding infinite cylinder of the same diameter. These surface currents have also been separated into a travelling wave component along the cylinder's axial

43 direction and a circumferential component [70]. The expressions for these currents are given by J e (, Z) = J c ()e-ikOcOS3z The circumferential components of surface current density are Je(q5) = Yo(sin qx - cos y) +m=~~_(-i)m{hzJm(xo) + BmHml)(xo)}eimk ko i E+oo2 + (-i) {ko sin/ [ezJm(xo) + AmH()'(xo)] (35) iko sin2 0 - oo +imcos~ [hzJm(xo)+ BmH$I(x)] }(im and J() = -(sin x-cos y) =_0o(-i)m{ezJm(xo) + AmHm)(xo)}Ceim k o2Z / m=-o (-i)m {ko sin/3 [hzJm(xo) + BmH,() (Xo)] (3.6) -imcose [ezJm(xo)+ AmHm() (xo)]}eim Pi where 3 is the angle of incidence (see Figure 3.1) and xo = koa sin 3, with a being the cylinder radius. In this case, the expression for the electric Hertz vector potential simplifies to H=iZo r2l2r' rL eiko(lr-r'l-cosazI) 4'e o J (q$) e adz'dO' (3.7) e 4rIko Jo [r - r' I If the observation point (x,y) satisfies the condition ko sin O/3 T+/jy~ > 1, the z' integration can be evaluated using the stationary phase approximation. The condition for the stationary phase is d ( r- I-cos z') = 0 which, in this case, implies,cos/3x -x')2 + (y- y,)2 sin'

44 L >>. z p= a Figure 3.1: Geometry of a long cylinder with equivalent electric and magnetic currents on the surface. and the stationary point is on the surface of a cone of half-angle / which contains the observation point. Performing the integration with respect to z' yields f~L,eko(I-r' -r'os)dZ = /I i( iko sin/-7r/4) e-iko COSz e-ikoasinOcos('-)'- (3.8) Applying the far-field condition in the x-y plane, the curl operator reduces to V x ikok, x where k,, = sin f(cos Ox' + sin q)) - cos Oz and the expression for the scattered field simplifies to E= k2(k8 x ks x Ie- Zok x Hm) (3.9) Substituting (3.5) and (3.8) into (3.7) the remaining integration with respect to 4' is accomplished with the aid of the following integral relations [90] fo2' e-ikop Bcos(~'-g)eim9' = 2 (-i)mJm (yo)eim

45 Stationry Points Primary Scatteri\ Cone X Figure 3.2: Geometry of dual cylinder configuration and second-order interaction. -r{~os~' c-ikopB s('o(,_)eim,,- - 2~.(-i)m[i{ Csin }J(Yo) + { s}J}(y oJ (Y)] with B = -[( k-1i) * ]2 + [(k -I) *..]2}1/2 and yo = koaB. The definitions of. and p are given in Figure 3.2. After lengthy algebraic manipulations the scattered field from the first cylinder is found to be E= - F(q)H(l)(ko sin PP)e-ikOzCoS (3.10) in which -1 +_ sF() =;2 j1 (-l)m[Am(ks x c, x )+ Bm( x B)]eim (3.11) The coefficients Am and Bm are given in terms of components that are TEz and TMz: Am - CM Ei. + iCm Hi. Bm= CE H' *- iCm Ei * z (3.12) and the expressions for aTE, C0M and Cm are given in [64].

46 Figure 3.3: Integration contour in complex -y-plane. 3.2.1 Second Order Interaction Referring to (3.10), the Hankel function of the scattered cylindrical wave may now be expressed as a continuous spectrum of plane waves H()'(ko sin /3p) = 1 e ikn in+sinll)dy (3.13) where the contour of integration r is shown in Figure 3.3. Each scattered wave from the second cylinder as a result of the plane wave spectrum excitation from the first cylinder can be summed by superposition to obtain the total scattered field from the second cylinder. In regard to the stationary phase approximation, only a subset of stationary points on the surface of the first cylinder act as sources of illumination for the second cylinder. Stationary points not located on the cylinder surface give no contribution to the secondary scattered wave. We have assumed that since the length L of the first cylinder is much greater than the excitation wavelength, most of the primary scattered fields are contained in the forward scattering cone. We now also assume that (L - pcos f) ~ A0, so the scattered field from each incident plane wave in the expansion is confined to the forward scattering

47 cone. The total scattered field from the second cylinder is then given by E = er (7)d- (3.14) E2 7 in which the direction of each incident plane wave is ki = sin /3(cos yx + sin y) - cos #z and the contour r has been described previously. Again each individual secondary scattered field can be described in terms of TE, and TMz components. Then, in the far-zone, these fields may be written as e,(7) - eikoki2 eikoR -i(L - cot)] sin V 2zll R [ r sin2 V +00 y (-1)m [A'(k8 xks x z) + Bm(ks x (3.15) m= -00 with Ao(L - - cot +] V k(L-cot/3) (cos~, - cos$) (3.16) b= sin/,(cos,~ + sin,.) - cos, and +00 A = 1 (-1)"[AnCTM + iBCm ]e',n nl=-oo00 +00 B = (-1)"[BnCE - iACm]eti (3.17) n= - oo where An and Bn are as given in (3.12) and cosI, f, =. Substituting (3.15) into (3.14) and using the change of variable y' = q + ry, an analytical expression for the secondary scattered wave in the far-zone is found to

48 be e _ eiko sin cos(r,- ) [-i(L - pcot#P) sin V iv r 7rsin2 ] V +oo Z (-1)m [A' (k x k' x 2) + B( k x -)] (3.18) H()(ko sin,3) e'm(~'-i) 3.2.2 Effect of Ground Plane The next step towards obtaining the scattered field from a layer of cylinders above a dielectric half-space is to consider the problem of two adjacent cylinders above a dielectric half-space. The level of difficulty involved in obtaining a solution for this problem is greater than for the previous problem due to the complexity of the Green's function in this case. The Green's function for this type of problem has an integral form and obtaining an analytical expression for this solution is impossible. A numerical solution with a great deal of complexity can be obtained by applying exact image theory [50, 67]. Here again we are approximating the form of the Green's function by assuming that the image of the source point is a source point located on the opposite side of the half-space interface (mirror-image point) and modified by the appropriate reflection coefficient. This approximation is very accurate when the source or observation point is not close to the interface and becomes exact when the ground-plane is perfectly conducting. To apply this simple approximation we decompose the incident and scattered fields into TE and TM components. The effect of the ground-plane will then be accounted for by modifying these field components by their appropriate reflection coefficients. The first-order interactions are shown in Figure 3.4. In the trunk-ground (TG) interaction, the primary scattered field is modified by interaction with the ground plane. The reflected field for this interaction is obtained by multiplying the TE and

49 TG Figure 3.4: First order interactions of a cylinder with the ground plane. TM components of (3.10) by the appropriate reflection coefficient to produce eikor -iL sin U Eg - eik~ (cosp-cosp)[ -ir sin2] U +oo00 Z (-1)m[AmRji(ks x k, x k z) + BmRi(k, x z)]eimp' (3.19) m= -oo00 where Am and Bm are defined in (3.12), o,3 and As represent the direction of the scattered wave and U = koL(cos /3 -cos $)/2. The other type of first-order interaction is ground-trunk (GT) and consists of modification of the primary scattered wave from a cylinder by the ground plane reflection coefficients. It is found in the same manner as above and is given by eE ik= r -iL sin U gt r _r sin2/9 U +00 Z (-1 )m[Am( x k x z) + Bm(k x z)] eim8 (3.20) m=-oo where AL = CTM Rl, Ei z + iCm RL Hi. z

50 Figure 3.5: Second order interactions of two cylinders with the ground plane. The second-order interactions are trunk-trunk-ground (TTG), trunk-ground-trunk (TGT) and ground-trunk-trunk (GTT). The mechanisms are illustrated in Figure 3.5. The TTG interaction may be obtained simply by taking the two cylinder scattering equation (3.18) and modifying the appropriate components by their respective ground plane reflection coefficients to obtain E iko eincos() i(L - cot/) sin V ei ttg r ir sin2 2 V +00 Z (i)m[A'$ R11 (kA x k x z) + BT' R. (kis x z) (3.22) m=-oo H(1) (ko sin /P) eim(~-$k) where A' and B' are given in (3.17) and V in (3.16). The TGT scattered wave is obtained by exciting a single cylinder with the expanded field of the TG interaction over its lower cot / portion in a manner completely analagous to what has already been done. The result is ikor ~ Es =,. ei Etgt eikr -ik sinpcos(s-) [-ipcot/] sinV' v r Lr sin2/3 V' +00 ~ (i)m[A" (k, x k, x ) + B" (k8 x.)] (3.23) m=-0o H1) (ko sin /p) eim(4f-d)

51 where V' k0o cot/3 (c os/3 os) (3.24) 2 and +oo00 A = (-1)n [A, RII CTM + iB, R1 Cm,]e " n=-oo +oo B1 = E (-1) [Bn R CTE - iAn R1 Cm]eio (3.25) n=-oo Finally, the GTT field is obtained using the same approach Ett e-ikos((L - cot sinV _iw ik0T rev'osinO co s(3Bc t/)] gtt r ir sin2 / V +00 (i)m[AI (k x x k ) + B"d (k: x Z )] (3.26) Hm )(ko sin /p) eim(ks4) with W ko(L + cot (os - cos (cos) (3.27) 2 and +00 A'" =. (-l)n [ATCM + iBr m]ein$ n=-oo +oo Bm= z C (-1)n[BC.E CT iAoCm]ei (3.28) n=-00 where A, and Bn are as given in (3.21) and V as in (3.16). The total far-zone scattered field to second-order is then the sum of the terms given above: E' = E- + E,, + Etg + Ett + Ett (3.29) The corresponding scattering matrix elements for each term are given in the appendix.

Z,Z~52 ki ki A A Y Y BROADSIDE CONFIGURATION ENDFIRE CONFIGURATION Figure 3.6: Broadside and endfire cylinder configurations. Wave incident in X-Z plane with incidence angle P. Measurements were performed at X-band (9.5 GHz) using a polarimetric scatterometer to validate the expressions derived for two adjacent cylinders above a ground plane. A pair of metallic rods 18cm in length and 0.56cm in diameter were arranged above a large metallic ground plane with various separations and orientations relative to the illuminating beam. The cylinder positions relative to the incident wave direction and the incidence angle for the cases of broadside and endfire illumination are shown in Figure 3.6. Figure 3.7 compares the first and second order theoretical predictions of the radar cross section for the two cylinders to the measured values in the case of 2cm separation and broadside confguration. Figures 3.8 and 3.9 compare the measured and predicted values of RCS for a separation of 4cm in the endfire configuration. The experimental data agrees well with the second-order results. The first-order approximation does not provide an adequate estimation of RCS for cylinder separations within this general range. It is seen that the first and second order

53 0.0 o, -5.0 —,.. ——. -- -10.0..m9 -15.0 __" "~ El _,,o'E O W Measred El,-'~ t )El HH Measured -20.0 -,, VV(2nd Order)?/ ---- - HH (2nd Order) -25.0 --- v (VVst Order)...... HH (1st Order) -30.0 I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.7: RCS of two rods broadside on metal plate at 9.5 GHz: 2cm separation. results do not differ significantly for ahh because the cylinders are relatively thin and therefore the horizontally polarized incident wave is only weakly scattered. Discrepancies between the measured and computed values of RCS are observed at the lower angles of incidence due to the effect of scattering by the cylinder end caps which is ignored in the theoretical formulation. 3.3 Monte Carlo Simulation Having validated the expressions derived for adjacent pairs of cylinders, we now attempt to obtain the scattering properties of random collections of such cylinders above a ground plane. For a given arrangement consisting of many cylinders, the solution of the scattering problem can be obtained to second-order by computing the single and pairwise interactions for every cylinder in the ensemble. The statistical properties of a random medium comprised of such scatterers are simulated by

54 0.0 -10.0,f -.... -20.0 It~,~I -30.0 _:!-40.0 O Measured 4-0._... - 2nd Order ------ 1st Order -50.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.8: RCS of two rods endfire on metal plate (vv-pol.) at 9.5 GHz: 4cm separation. application of the Monte Carlo method. The principle of the Monte Carlo simulation based on the second-order algorithm is as follows: 1 An ensemble of randomly positioned cylinders is generated using a random number generator. In this case the cylinder positions are uniformly distributed within a circular area. The number of cylinders used is dependent on the specified number per unit area and the area of the circular region. 2 The scattering is computed for all cylinders and between all pairs of cylinders within the ensemble up to second-order. 3 The ensemble is re-randomized and the scattering recomputed as discussed

55 0.0 -10.0 -.,r " - " -20.0. 3 -30.0:1, ~~~~40.0 _O~~~ 4~ 3 Measured -40.0 4 Rii -------- 2nd Order s ------ Ist Order -50.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.9: RCS of two rods endfire on metal plate (hh-pol.) at 9.5 GHz: 4cm separation. above. The number of independent samples is chosen so as to make the variance as small as possible within limits depending on the computing time. For the cases analyzed in this work the sample number is greater than one hundred. 4 The values of the scattering coefficients (oaV,, ah) are found from the ensemble average. The same is true for the co-polarized phase difference ( and the degree of co-polarized phase correlation a [66]. Measurements of radar backscatter were also made on random collections of cylinders distributed within a circular area as described for the Monte Carlo simulation as an experimental verification of the simulation. Identical metal cylinders 18cm in length and 0.56cm in diameter were uniformly distributed within a 60cm diameter

circular area with densities of 70, 100, 140 and 180 per square meter. Radar measurements were made with an X-band scatterometer at incidence (elevation) angles relative to the vertical axis ranging from 20 to 60 degrees in 5 degree increments. Samples were generated by arranging the cylinders on four separate thin plastic disks with independently generated holes for supporting the cylinders. The edge diffraction contribution of the ground plane to the overall RCS of the cylinder/ground plane system was removed by measuring the ground plane without the cylinders present and coherently subtracting this measurement from that of the combined system. Each disk was rotated in 10 degree increments to create distinct aspects for the radar. In this way 144 independent samples were generated at each density level. The final measurements were tested for correlation to ensure that the sampling was independent. The experimental setup is illustrated in Figure 3.10. The correspondence of the simulation to the measurements of a& as a function of incidence angle is shown in Figures 3.11 through 3.14. The first and second-order Monte Carlo simulations agree very well with the experimental data. Apparently, for the type of random medium considered in this study, the effect of multiple scattering between cylinders is averaged out as far as the magnitude of the radar backscatter is concerned. This is why the first-order results agree well with the co-polarized radar cross section measurements. Figures 3.15 through 3.22 illustrate the agreement between measurement and simulation for the co-polarized phase statistics. It should be noted that for both oa and C, inclusion of the second-order terms provides the correct phase statistics while the first-order scattering theory is significantly in error. Since a is a sensitive function of the degree of multiple scattering within a medium, first-order theory incorrectly predicts a value of unity independent of the number density of particles.

57 Radar Field of Metallic Cylinders Ground Plane Figure 3.10: The measurement of a random ensemble of vertical cylinders. 3.4 Comparison with Radiative Transfer In a medium in which the individual particle scattering albedos are small and the particles are in the far-zone with respect to each another, one would expect that there would be little difference between the Monte Carlo simulation results and those of radiative transfer. Because the scattered fields are added coherently in the Monte Carlo simulation, the computed values of a" by this method should be 3 dB above those computed by the incoherent addition of power as in radiative transfer theory for reciprocal scattering mechanisms such as GT and TG. The medium described above in the experimental section of this article consists of metallic cylinders that are long compared with the wavelength of radiation in the medium, however the diameters of the cylinders are fairly small compared with the excitation wavelength. In this case the

58 15.0 5.0 00/ O a, Measured 0 E a Measured ----- da* (2nd Order) -15.0, ~~~~~~-15.0 * ~ —- a, (lst Order) -.-.- o, (2nd Order) b —- a, (st Order) -25.0. 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.11: Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 70cyls/m2. far-field condition is satisfied for the wave polarized in the dimension transverse to the cylinder length (hh Polarization) and propagation for this polarization is dominated by single scattering. This effect can be seen in Figure 3.23 which represents hhpolarized backscatter from a fairly dense medium with 140 rods per square meter. A vertically polarized wave travelling in this cylinder medium encounters particles that are both strong scatterers and that are in the near-field with respect to each other. Because of these conditions the correlation distance for the vertically polarized field in the medium is significantly smaller than the length of a cylinder, yet it is larger than the distance between particles. Therefore the local plane-wave approximation is no longer valid in this regime and, in addition, there is a significant degree of coupling between particles. This would be expected to affect radiative transfer results in

59 15.0 0m 0 5.0 --- --- (2nd Ord" ) 1-5.0 -. (I stOrder) - - - -.?bb (2nd Order) e., (stOrde) -25.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.12: Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 100cyls/m2. 15.0I v_, -< iL____ —-"d, —. <gx <!>_______' 5.0 _-. a:, -5.0 -15.0......?,, (It Order).... a;e (2nd Ordea).ba (Ist Order) -25.0 I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.13: Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 140cyls/m2.

60 15.0 ~~5.0 - 0r o r_,,'I 5.0 L. to O co,, Measured El a-~, Measured ---- a, (2nd Ordr) -15.0 ------ ad, (IstOrder) - -- - op (2nd Order) - b (1stOrder) -25.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.14: Backscattering coefficient of a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 180cyls/m2. 1.0... 0.- 0 0.8 o-, - 0.6 0.4 0.2o Measured 0.2-------— 2nd Order 0.0 - - 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.15: Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 70cyls/m2.

61 1.0 -, 0 0 0.8 - 0.6 0 0 <0.4 0.2 0 Measured 0.2 -------- 2nd Order 0.0 I I. I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.16: Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 100cyls/m2. 1.0 0.8 20.0 30.0 40.0 50.0 60.0 0.4 o Measured 0.2......... 2nd Order 0.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.17: Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 140cyls/m2.

62 1.0 - 0.8, 0 0 0.6, - <0.4 - 0.2 - 0O Measured - 2nd Order 0.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.18: Co-polarized degree of correlation for a layer of uniformly distributed metal cylinders above perfectly conducting ground plane. The density of scatterers is 180cyls/m2. 360.0 0 Measured..-...... 2nd Order 270.0 - st Order 0 180.0 90.0 a 0.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.19: Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 70cyls/m2.

63 360.0 0 Measured -—. 2nd Order 270.0 ------ 1st Order 180.0 ~~90.0 i-:,0 0 — 1 0.0... 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.20: Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 100cyls/m2. 360.0 0 Measured -. ——... 2nd Order 270.0 ------ 1st Order 180.0 90.0. -..-0-'"'.-... 0.0. " 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.21: Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 140cyls/m2.

64 360.0 0 Measured -. 2nd Order 270.0 ------ 1st Order 180.0 90.0 90.-"o -.-. 0.0' 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.22: Co-polarized phase difference for a layer of uniformly distributed metallic cylinders above a perfectly conducting ground plane. The density of scatterers is 180cyls/m2. two ways. The phase matrix is utilized by radiative transfer to generate the source functions which drive the coupled set of differential equations for the upward and downward travelling radiation intensities. This phase matrix is linearly dependent on the height squared of the cylinders. If the mechanism of cylinder excitation is not by a plane wave then field decorrelation in the medium can make the effective height of the cylinders to be smaller than their actual length. This means that radiative transfer theory would overestimate the phase matrix of the cylinder medium. On the other hand, the computation of the extinction matrix, which accounts for attenuation of the wave intensity as it propagates in the medium, is linearly dependent on the height of the individual cylinders. If the effective scattering matrix for cylinders in the medium is actually smaller than would be expected on the basis of the local planewave approximation, the extinction matrix for the medium would also be smaller than

65 15.0 5.0 o~o..-'" " J~ ) o, O Measured -. Monte-Carlo (2nd Order) -15.0 ----- Monte-Calo(lst Order) - Rad. Tras. (2nd Order) Rad. Trans. (1st Order) -25.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.23: Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 140cyls/m2. that normally used in radiative transfer. Because extinction is an exponential process it has a large effect on the computation of the vv-polarized radar cross section. Figure 3.24 shows that while first and second-order Monte Carlo simulations agree fairly well with measured data in the computation of ao, first-order radiative transfer gives an estimate that is significantly low. This is also shown by the fact that the first-order radiative transfer solution becomes worse as the angle of incidence increases which is consistent with overestimation of the extinction and which more than compensates for overestimation of the phase matrix. Even though second-order radiative transfer provides a solution that is more in harmony with the measured results, it would be expected that if the second-order theory were entirely correct this solution should maintain a level consistently about 3 dB lower than the measurements and Monte Carlo results. In addition, in light of the trend illustrated, it seems unlikely

66 15.0.....~_f... ~- ~ ---.-............ 5.0 -5.0 0 0 Measured. —------ Mont-Carlo (2nd Order) -15.0 ------ Monte-Carlo (IstOrder)...... Rad. Trans. (2nd Order) - Rad. Trans. (Ist Order) -25.0 I I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.24: Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 140cyls/m2. that the radiative transfer solution will converge to the correct level as successively higher order terms are added since all further contributions to the net scattered intensity are guaranteed to be positive and the second-order result is already too high. It is also to be noted that in all cases evaluated by these authors the first-order Monte Carlo simulation provides an excellent estimate of radar cross section. If the radiative transfer approach is valid under these conditions one would wonder why the first-order R.T. theory does not give results consistent with this finding. The degree of correlation a for the phase distribution PDF is sensitive to multiple scattering effects because it is the multiple scattering that produces phase decorrelation in random media. Figure 3.25 shows that the second-order Monte Carlo simulation gives good agreement with experimental measurements indicating that multiple scattering to second-order is significant in the medium. The second-order radiative

67 1.0 0.8 - 0.6 0.4 o Measured 0.2 ------—. Monte-Carlo.....- Rad. Trans. 0.0 I. I I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.25: Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 70cyls/m2. transfer theory gives an erroneously high degree of copolarized phase decorrelation. This demonstrates the overestimation of the phase matrix previously mentioned. In Figure 3.26 both the Monte Carlo simulation and radiative transfer are in agreement with the measured data. It is evident that the co-polarized phase difference is not significantly affected by parameters that reflect the differences between these two methods of modeling the canopy. 3.5 Conclusions In this chapter a Monte Carlo scattering model for the trunk layer of a forest canopy has been developed which takes into account scattering effects up to secondorder. Experimental data have been presented for the purpose of validating the two-cylinder scattering solution, and the results of a Monte Carlo simulation based on this solution have also been presented and compared with measured results. First

68 360.0 270.0 _ Measured.. Monte-Carlo ------ Rad. Trans. 180.0 90.0 0.0. I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 3.26: Comparison of RT model with Monte Carlo simulation and measured data for the cylinder layer with a density of 70cyls/m2. and second-order radiative transfer model solutions for the same medium have been given and compared with those for the Monte Carlo simulation. It is verified that the radiative transfer model provides incorrect results under conditions such that its basic assumptions are violated. This occurs in media for which the size of the particles is large compared to the wavelength in the medium which causes the illumination to be non-uniform and/or the medium is dense and therefore the near-field interaction becomes significant. The source of this problem is attributable to several factors: (1) the extinction matrix computed by the radiative transfer model is overestimated because multiple scattering reduces the actual coherence length for fields within the medium making it less that the size of the scatterers, (2) the cross-coherence terms that are accounted for in the Monte Carlo simulations are absent from the RT model, (3) the source functions of RT theory do not properly describe the nature of the

volume scattering that occurs within this type of medium.

CHAPTER IV ELECTROMAGNETIC SCATTERING FROM TWO ADJACENT OBJECTS 4.1 Introduction Typical vegetation canopies consist of objects such as trunks, stems, branches and leaves or needles and, in general, vegetation tends to have some complex structural features. However, analytical solutions for the problem of scattering of electromagnetic (EM) plane-waves by objects exist for only a limited number of canonical geometries. If the scattering body is inhomogeneous, or the polarization and phase front of the illuminating field is non-uniform, analytical solutions of the vector wave equation do not exist even for canonical geometries. Almost all models that are currently being used for the analysis of EM scattering from collections of discrete scatterers rely on the single scattering properties of the constituent particles [45, 72, 42, 82]. However, when the sizes and/or number densities of the scatterers become large enough that they are in the near-field of each other, solutions based on their single scattering properties are no longer valid. Certain types of vegetation canopies such as forest stands and some agricultural fields have high number densities of strong scatterers [12, 105, 56]. In addition, the particle sizes in these canopies are large enough so that adjacent scatterers are not in each others 70~~~~~~~~~h

far-field zone, especially in the microwave region. To model vegetation of this type, it is necessary to be able to treat electromagnetic interactions between particles that are not only non-plane-wave in character but have non-uniformities in amplitude, phase and polarization. As has been stated, it is impossible to find exact analytical solutions for this type of problem except in a small number of cases [49, 5]. In some few circumstances it may be possible to obtain an analytical solution by employing a plane-wave expansion technique [69] or some other specialized approach. Even so, such solutions may yield results that are difficult to evaluate or have tedious multiple integrations that must be done numerically. This is obviously a distinct disadvantage when the desired end result is to simulate EM scattering from a dense random medium. EM modeling of vegetation canopies usually involves the construction of simplified geometrical representations for the constituent scattering elements [95, 96, 40, 39]. When this is the case, it is a fairly simple matter to obtain expressions for the firstorder scattered field using the single-scattering properties of the isolated particles when the primary excitation is a plane-wave. To obtain the secondary scattered field from interacting particles it is necessary to account for illumination of the secondary scatterer by the scattered field from the primary scatterer. This chapter presents a technique for obtaining the secondary scattered field analytically by employing the reciprocity theorem. The technique is then applied to obtain an analytical solution for bistatic scattering from a cylinder-sphere pair. This cylinder-sphere interaction has some importance because, along with the electromagnetic coupling between pairs of cylinders, it provides a basic building block from which the EM scattering properties of a heterogeneous two-component forest canopy may be simulated. The results of analytical field calculations for the cylinder-sphere pairs are then compared with

numerical computations based on the method of moments. 4.2 Secondary Scattered Field from Reciprocity In this section, a procedure utilizing reciprocity for evaluation of the secondary scattered field from a particle when illuminated by the primary scattered field of another adjacent particle is outlined. The reciprocity theorem simplifies this evaluation significantly when the observation point is in the far-field zone of both particles. Derivation of the expression for the secondary scattered field of perfectly conducting particles is slightly different from that of dielectric particles. First we consider perfectly conducting particles. Suppose the incident field induces a surface current density J1 on the surface of particle #1 in absence of particle #2. The objective is evaluation of the scattered field from particle #2 with J1 as the excitation source. Suppose the electric and magnetic fields produced by J1 (as an impressed source) in the presence of particle #2 are denoted by E1 and H1 respectively as shown in Figure 4.1. Obviously E1 and H1 satisfy Maxwell's equations everywhere in the medium. Now, let us consider another situation where the source J1 is removed and an infinitesimal current source J, given by: Je- = _ (r- r,) is placed at observation point P as in Figure 4.2. The electric and magnetic fields produced by the elementary current source are denoted by E,2 and H,2 and also satisfy Maxwell's equations. Applying the reaction theorem [30] to fields E1, Ee2 H1, and He2 over the entire medium, it can be shown that Js+s (ElxH, Ex H)2-E2x ds=-A J1 Ee2 ds + E1 (41)

73 Observation S Point E,11 Particle #1 S2 Particle #2 Figure4.1: The current induced on particle #1 by the source produces a scattered field from particle #2. Observation Infinitesimal Point Current Source A A t/S Ee2 S2 local plane wave Particle #2 Figure 4.2: Particle #1 is removed and an elementary current source is placed at the observation point.

where S. represents a closed sphere at infinity. Since at distant points Ej = ZoHj x ii, where j = 1, e2, the surface integral over SO vanishes. Also the surface integral over the surface of particle #2 vanishes since in x E = O over a perfectly conducting surface. Since the elementary current source Je is in the far field zone of particle #2, Ee2 can be approximated by the scattered field of particle #2 when illuminated by a plane-wave. Noting that J1 is a function of the incident wave polarization (vii, hi) and the polarization of the elementary source p can also be chosen to be either Vi or h,, the expression for the scattered field is given by EP s EE= js Ee2(is, I,) J(i, hi) ds (4.2) The scattered field from particle #1 when illuminated by the field generated by the induced current on particle #2 can be obtained in a similar manner and is given by EP. Ee = E,( h). J2(Oi, h,) ds. (4.3) Here, Ee, (v,, h,) is the sum of the radiated field from the elementary current source and the scattered field of particle #1 when illuminated by this elementary source. The source is characterized as possessing either a vertical or an horizontal polarization state and is located at the point of observation. The primary induced currents J1 and J2 are functions of the incident wave polarization, thus, by selecting the incident polarization to be vi or hi, the scattering matrix elements can be computed. It should be noted that E' includes the first-order scattered field of particle #1 and secondary scattered field of particle #2. Identically, E' includes the first-order scattered field of particle #2 plus the secondary scattered field of particle #1. Therefore Es + E2 is the total scattered field of both particles up to second-order.

Now let us consider the case in which both particles are dielectric objects. The incident wave induces a polarization current J1 in particle #1 in the absence of particle #2. When particle #2 is placed at its location with no incident wave present, the volumetric current J1 induces a volumetric current J12 in particle #2. The total fields in this case will be denoted by E1 and H1. In a second experiment, we place the elementary current source Je at the observation point as before and remove the current source J1 (particle #1). The elementary current then induces a volumetric current Je2 in particle #2 for which the electric and magnetic fields will be denoted as Ee2 and He2. The currents induced in particle #2 may be expressed in terms of the total electric field and relative dielectric constant (62) of particle #2. They are given by Jl2(r) = -ikoYo (e2 - 1) El(r), r E V2 (4.4) Je2(r) = -ikoYo (e2 -1) Ee2(r), r e V2 (4.5) where ko and Yo are the wave number and characteristic admittance of free space respectively, and V2 is the region occupied by particle #2. Application of the reaction theorem over the entire medium results in j (E E xe2-e2 nds= -- J. Ee2 dv- J12.e2dv + J J2e E1 dv + p- E1 The integral on the left-hand side vanishes as before, and, by substituting (4.4) and (4.5) into the second and third integrals on the right-hand side, it can be shown that the last two integrals in the expression given above cancel each other. Thus, the sum of the primary scattered field of particle #1 and the second-order scattered field

of particle #2 is given by pi.E1- J. Ee2 dv (4.6) Similarly, the sum of the scattered field of particle #2 and the second-order scattered field of particle #1 is given by p' E2 = J2' Eel dv. (4.7) In (4.7), J2 is the volumetric current induced in particle #2 by the incident wave in the absence of particle #1 and Eei is the total field of the elementary current located at the observation point in the presence of particle #1. When one of the particles is dielectric and the other one is metallic, the expressions for the scattered fields can be obtained in a similar way. Let us say that particle #1 is metallic and particle #2 dielectric, then the expressions for the scattered fields are given by (4.2) and (4.7) respectively. Conversely, if particle #1 is dielectric and particle #2 metallic, expressions (4.6) and (4.3) give the scattered fields. 4.3 Electromagnetic Scattering by a Cylinder-Sphere Pair The expressions for the scattered fields from the two particles as derived in the previous section are very general and can be applied to any particle pair with known geometries and dielectric properties. In this calculation only the scattered fields and induced currents of isolated particles when illuminated by a plane-wave are required. The sphere and cylinder are among the few geometries for which an exact analytical scattering solution is known. Additionally, as mentioned previously, a collection of randomly positioned spheres and vertical cylinders above a ground plane can be used to simulate a heterogeneous forest medium which includes the effect of multiple scattering between canopy components.

Analytical evaluation of the integrals (4.2), (4.3), (4.6) and (4.7), even for the cylinder-sphere pair, is not possible, and one must resort to numerical methods. However, under some physical conditions, approximate analytical expressions may be derived. In reality, a tree trunk is very long (L>> A) and, therefore, other tree trunks and leaves will tend to fall within its near field region. In this section an approximate analytical solution will be derived for a cylinder-sphere pair. The assumption is made that the cylinder is in the far-field region of the sphere. However, the sphere is assumed to be in the near-field of the cylinder with respect to the cylinder's longitudinal dimension and in the far-field region of the cylinder's transverse dimension. If the radius of the cylinder and sphere are denoted by a, and a, respectively, and the cylinder length is represented by L, then the conditions previously specified may be stated mathematically as: 2a2 L2 2a2 -> C > (4.8) where - is the distance between the cylinder axis and the sphere center. Suppose a plane-wave, whose direction of propagation is denoted by ki, is incident on a cylinder-sphere pair and is given by: Ei =i eikk'r. The cylinder axis coincides with the z-axis of the Cartesian coordinate system and the sphere center is located at r = pcos qx+p- sin qy+zz'. The observation point is located at a distance ro in the direction k,. The geometry of the problem is given in Figure 4.3. The current induced on the sphere when illuminated by a plane-wave can be easily

Y Figure 4.3: Geometry and coordinates of the cylinder-sphere pair. computed by the standard method of separation of variables [4]. Since the cylinder is assumed to be much longer than the excitation wavelength, the current induced in this finite length cylinder can be approximated by that in an infinitely long cylinder of the same radius and electrical properties [100, 41]. The field generated by the elementary current source located at the observation point, over the volume (surface) of each scatterer when the other one is absent is calculated as follows. Let us first consider the case where the sphere is absent. The field generated by the elementary current J, = ip(r- ro) is composed of two components. The first component is the direct contribution of the current source and is given in the far zone by Eed(r) =-ikoZo eikOO e-ik0okr x k, x. (4.9) The second component is the scattered field from the cylinder when illuminated by

79 the radiated field of the elementary current source. This illuminating field can be approximated locally by a plane-wave propagating along the -k, direction. Since the point r is in the near field of the cylinder which satisfies condition (4.8), the scattered field is given by [69] E (r) = —ik~Zeik~~ F(q - q$,)H(1) (ko sin O,p) e-iko cosOz (4.10) 4irr0 where (p, 0, z) is the cylindrical coordinate of position vector r, H(1) is the Hankel function of the first kind and zeroth order, and 8, and 4o are the spherical angles specifying the unit vector k8, that is k, = sin 0 cos qx + sin 9, sin q,y + cos Oz (4.11) The expression for the vector quantity F(q - Os) in (4.10) is given by F(- n2 -} )1 (+1o)m [A'(P' x k' x z) + B'(k' x z)] eim(68) where k' = sin 0, [cos X * + sin ~] - cos 0, z and Am C= M ks x (k, x ) + iC(k, x ) BA, - CmEk,x iCk, kx (k~ x *z The expressions for CTM, CTE, and Cm are given in [64]. The direct scattered field from the sphere and the secondary scattered field from the cylinder can be obtained from (4.2) or (4.6) depending on whether the sphere is perfectly conducting or dielectric. The first component of this scattered field is due

80 to the excitation Eed(r) which yields the direct contribution from the sphere. This is specified by p. E = J J,(r). Eed(r) dv sph. Using the far-field expression for Eed(r), it can be shown that ~ikoro = eko=o-(-ko2)[ksx(k x Z0x k Js(r)e-~ko dv)].ip ek e-ikokrs.,(kik, s). i (4.12) where S,(ki, ks) is the bistatic scattering amplitude of the isolated sphere. The farfield amplitude of a sphere is given by S,(ki, ks) = (4.13) I k, x 1 [1i k, S,(O;)(k, x ki) x k, + e.l (k x ki)2(os) (k, x ki)] where Si(09) = -iZ(-1)nn(n + TAn ncsin98) +B 16 Pn()(cosit)} S2(,) = iy..(_1)n 2n + 1 fA 0P( )(csO:)+B P+ (cosO ) } n(n + 1) a00n sin f In these expressions Pn1) is the associated Legendre polynomial, An and Bn are given in [90] and cos 0' = -ki. k,. The second component of the scattered field is referred to here as the sphere-cylinder interaction and is given by I'. Elsc = J J(r) Eec(r)dv. (4.14) sph.

81 Integral (4.14) can be evaluated analytically, keeping in mind the conditions on the dimensions of, and distance between, the particles as specified in (4.8). Under these conditions the angle subtended by the sphere is small, therefore the following approximation is made F(q$-q0s) -F(q$-qOs) which is a constant vector and comes out of the integral. Also the Hankel function is approximated by its large argument expansion, and we have H(1) (ko sin 0,p) e-ik ~SZ os,k 2 e-ir4 exp[-iko(- sin 0s + cos Oz) r] /7rko sinS (4.15) Again, by approximating p and p^ in the amplitude and argument of (4.13) by p and p respectively, (4.5) can be written as Ho1) (ko sin OBp) e-ik~o cosGz Ho)(ko sin Os,)e-ikoc~~8z~ ex p [-ik. r'] where lk = -sinO, + cos 08, r= r —r Thus, the sphere-cylinder interaction term is obtained from P. -lc 1 e k~Hk)(ko sinO,P)e-ikocos Oa{[ ( As i (E, _.4)] El EHPc (c- si x z

82 where E'(ki, k) is the bistatic scattered field amplitude of the isolated sphere. Now we consider the case for which the cylinder is absent and evaluate the field generated by the elementary current source. As previously, the total field generated by the elementary current consists of two components. The first component is the direct field given by (4.9), and the second component is the scattered field from the sphere when illuminated by the direct field of the elementary current source. Noting that this current source is in the far-field of the cylinder-sphere pair, the sphere illumination appears locally as a plane-wave. Since the cylinder is in the far-field of the sphere, the second component, in the vicinity of the cylinder axis, is given by ) -ikoZo eikoo e.-ikk.r 4ir ro r S(' where r' is the distance between the sphere center and the point at r, that is r' =I r-r I and P' = The direct scattered field from the cylinder and the secondary scattered field from the sphere can be obtained from (4.3) or (4.7) depending on whether the cylinder is perfectly conducting or dielectric. The first component of this field is the direct contribution from the cylinder alone and is given by P. E2c J Jc(r). Eed(r) dv cyl. As in (4.12) it can easily be shown that eikoro p = E2c =- Sc(k,, ik,). (4.16) where S,(ki, i:) is given by

83 SC -iL sin U s A {At ^ x) + Bm (ks x }im(O,0) -iL sinU i U +O with koL U= ~- (cos Oi-cos 0). In this expression Oi is defined by cos i- ki = z. Also, the expressions for At and B/ are given by Am Cm ei Z + iCm(ki x ej) z Bm Cm (ki X ei Z-iCm(ei Z For analytical evaluation of the cylinder-sphere interaction, it is noted that the current induced in the cylinder has a progressive phase factor along the cylinder axis, that is Jc(r) = j(p, 4) eiko cos iz The second term in (4.7) is given by -ikoZo eikoro _oeiko Ir.i 4ir r0 Jujo rE -2 c.9Zo ei9F~ -_ikok.r i(+ p)eiko cosOil * S,(-kS7, rm) dzds (4.17) where the first integral is over the cylinder cross-section. By rewriting r- | in cylindrical coordinates, i.e.

84 r- = /I P-,I + (z- z)2 and noting that ko I p - j j> 1, the stationary phase approximation can be used to evaluate the z integration. The stationary point of the function F(z) =I r - | I + cos Oiz is at Z = z - cos O I p - I and (4.17) simplifies to -ikOZO eikCr~ -ko0, E2c= ikZ e-k~o r H(')(ko sin0i I p-P I)Jc(r).S,(-k, r',)pdpdo (4.18) p -o4ar ro 0 where ri is the unit vector r' evaluated at the stationary point. Noting that I p I is much smaller than I p F S - in Oj- cos - sin Oi 0 - cos oi Z Under this approximation S,(-k,, r,) is not a function of the integration variables. Therefore (4.18) can be written as pi E2C = iko e-ikok0 r H(l)(ko sin Oij) So(-/k, r'). fr Jc(r)e-iko~nOiP pdpdck Note that with values of z for which Z, > L or Z, < O (stationary point outside the cylinder surface) the scattered field E2,, is negligible and can be ignored. Using (4.13), the far field amplitude of the sphere takes the following form SJ(-k,,')= I k 1 [x * r' S(cos (k, *')) r x (rS x Ic) -p (kIc, x r) S2 (cos-1(k~s* r.)) r' x /q] (4.19)

85 Substituting (4.19) into (4.18), after some algebraic manipulation it is found that p' E2cs = — e~~ e~; I OH()0o sin 0i ). [p;S (cost-(ks r',))E~(k,,-r' ) =ie ikoro e ikokA where E,(ki,-r,') is the bistatic scattered field amplitude of an infinitely long cylinder with the same radius and dielectric constant as those of the finite cylinder. That is, 1 +00 ~k E~(k,-%) - sin 0j m-E- 1A(- cos i i + sin Ol) + B eim('-k) where = oz x p. The derivation of the scattered field for a cylinder-sphere pair can be easily generalized to a cylinder and an arbitrary scatterer so long as the expression for the scattered field of the isolated scatterer is known and the dimensions of the scatterer satisfy the conditions specified in (4.8). 4.4 Numerical Results The theoretical development presented in the previous section for second-order scattering from the cylinder-sphere pair has been validated using the Numerical Electromagnetics Code (NEC) [6] which is a computational package based on the method of moments. This approach was chosen because we were interested in the bistatic scattering behavior of the pair, particularly in the forward specular cone, which is quite difficult to obtain experimentally. The forward specular cone is referred to as the set of azimuthal angles for which the scattered wave-vector lies on the conical surface of revolution generated by rotating the incident wave-vector around the z-axis as

86 z Elevation Angle X I Azimuthal Angle A Figure 4.4: Scattering geometry and angles for forward specular cone scattering. shown in Figure 4.4. This sub-domain of the scattering pattern is particularly important for simulating the interaction of electromagnetic waves with vertical structures above smooth surfaces. In this case, the radar return is dominated by the scattered field in the specular cone. The validation for the vv and vh polarization states were made with a model consisting of a cylinder 18.0 cm in length and 0.1 cm in diameter with a finite conductivity of 100 mhos/m, and a perfectly conducting sphere with ka = 1.69 at an excitation frequency of 9.25 GHz. The cylinder was chosen to be of finite conductivity because this damps the axial standing wave pattern that exists on a finite length perfectly conducting cylinder. Our finite length cylinder model does not need to account for this standing wave behavior because in all real vegetation, cylindrical structures are composed of lossy dielectric material and do not support standing waves of signifi

87 -20.0 -E 0-25.0 -30.0 Second-Order Theory -------- First-Order Theory A MoM Computation -35.0....I 0.0 90.0 180.0 270.0 360.0 Azimuth Angle (deg) Figure 4.5: The vv-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p/=18cm, q = 1800 and z=9cm relative to the cylinder base. The incidence and scattering angles are 143~ relative to vertical. The frequency is 9.25 GHz. cant amplitude. A cylinder having a small diameter as compared with the excitation wavelength was used because the version of NEC we have only provides for finite conductivity in thin wire structures. The number of unknowns for the thin cylinders was on the order of 10 per wavelength or a total of about 60 for the 18cm length. The sphere was composed of variable segmented perfectly conducting rectangular patches as described in [7]. In general, the relative configurations of the cylindersphere pair and the scattering patterns were chosen so as to present as great a contrast as possible between the first and second-order scattering behaviors. The angle of elevation and the azimuthal angle are defined in Figure 4.4. The plane of incidence is the x-z plane and the azimuthal incidence angle is 180 degrees. The cylinder is always located at the origin and the relative cylindrical coordinates of the sphere are

88 -20.0.. -25.0 -30.0 Second-Order Theory -......... First-Order Theory A MoM Caomputation -35.0... 0.0 90.0 180.0 270.0 360.0 Azimuth Angle (deg) Figure 4.6: The vv-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, q = 45~ and z=9cm relative to the cylinder base. The incidence and scattering angles are 145~ relative to vertical. The frequency is 9.25 GHz. -20.0.................ider Theory.... 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 4.7: The vv-polarized bistatic scattering cross section of a cylinder-sphere pair. The s-polhere is at ati2.5cm, s = 450 and 9cm relative to the cylinder base. The incident wave azimuth angle is 180~ and the scattered wave azimuth angle is 130~. The frequency is 9.25 GHz.

presented. Figure 4.5 shows the vv-polarized azimuthal pattern for a cylinder-sphere pair. In this case the cylinder and sphere are located far enough apart to be effectively isolated. This figure demonstrates that the single scattering models for the sphere and cylinder are in agreement with the moment method computation to within less than about 0.2 dB and establishes a baseline for comparison. Figures 4.6 and 4.7 illustrate the vv scattering behavior for the case where the sphere is close to the cylinder. In Figure 4.6 the relative cylindrical coordinates of the sphere are (i, ~, z) = (2cm,45~,9cm), while in Figure 4.7 they are (b, ~, z) = (2.5cm,450,9cm). It may be seen from these two figures that the vv second-order result provides a reasonable approximation and is in agreement with the moments method data to within about 0.4 dB over the angular range. The maximum difference between the first and second-order cross-polarized (vh) response occurs in regions close to 0 and 180 degrees in azimuthal angle. At these two points the first-order cross-polarized response disappears while the second-order response is low but non-zero. Figure 4.8 illustrates the difference between the first and second-order scattering behavior for the case having relative cylindrical coordinates of (i, X, z) = (2cm,450,9cm) and an azimuthal scattering angle of 350 degrees. This provides good contrast between the scattering orders, and the scattering amplitude is strong enough that the accuracy of the numerical computation is sufficient for comparison. The agreement of the second-order analytical result with the moment method computation has a mean deviation of about 1 dB or so over most of the angular range. For verification of the hh-polarized response, it was necessary to use a thicker cylinder since the contrast between the first and second-order terms is insufficient

90 -40.0' -50.0 m -60.0o "-. Second-Orde Theory -- First-Order Theory A MoM Computation -70.0............ I........ -70.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 4.8: The vh-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, q = 450 and,=9cm relative to the cylinder base. The incident wave azimuth angle is 180~ and the scattered wave azimuth angle is 3500. The frequency is 9.25 GHz. -25.0 -35.0'"-, -- Second-Orde Theory....... Firt-OrderThry A MoM Computation -55.0 0.0 90.0 180.0 270.0 360.0 Azimuth Angle (deg) Figure 4.9: The hh-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, b = 45~ and z=9cm relative to the cylinder base. The incidence and scattering angles are 145~ relative to vertical. The frequency is 9.25 GHz.

91 -20.0 -30.0 0 -40.0.Y..CP I4; a** Ift- Second-Order Theory -50.0......... First-Order Tey - A MoM Computation -60.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 4.10: The hh-polarized bistatic scattering cross section of a cylinder-sphere pair. The sphere is at p=2.0cm, q = 450 and i=9cm relative to the cylinder base. The incident wave azimuth angle is 180~ and the scattered wave azimuth angle is 90~. The frequency is 9.25 GHz. with very thin cylinders. For cylinders with larger diameter to wavelength ratios than the one used for verification of the vv polarization case, NEC requires the use of a patch model. A patch model of a cylinder 18.0 cm in length and having a diameter of 0.55 cm was constructed using perfectly conducting patches. The model had 15 sides and consisted of over 800 rectangular patches. Because the patches were perfectly conducting, the axial standing wave made the model inappropriate for verification of the vv-polarized case. However, the standing wave seems to have a much smaller effect as far as the hh response is concerned, becoming significant only for elevation angles less than about 30~. Figures 4.9 and 4.10 present azimuth and elevation patterns for the cylinder-sphere pair having relative cylindrical coordinates (p,, z) = (2.0cm,45',9cm).

4.5 Conclusions In this chapter a general technique based on the reciprocity theorem has been developed for deriving the secondary scattered fields from a pair of objects. The general formulation has been applied to obtain approximate analytical expressions for the secondary scattered fields from a cylinder-sphere pair. The validity of the analytical results were verified by comparison with method of moments computations, and good agreement was obtained for both the co-polarized and cross-polarized bistatic scattering cross sections in the forward specular scattering cone. This chapter provides the basis for the construction of computational simulations of electromagnetic wave scattering from heterogeneous two-component vegetation canopies that include the effect of multiple scattering up to second-order.

CHAPTER V SIMULATION OF ELECTROMAGNETIC SCATTERING FROM A HETEROGENEOUS MEDIUM 5.1 Introduction In the application of radiative transfer (RT) theory to the modeling of forests and other dense vegetation cover, several fundamental conditions necessary to the validity of the model have been ignored. The RT model is based on the single scattering properties of the particles that constitute the medium. That is, the assumption has been made that particles are in the far field of each other and that they are illuminated locally by plane-waves. A tree canopy usually contains particles such as trunks and branches which are much larger in dimension than the excitation wavelength at microwave frequencies and above. These large scatterers are usually in each other's near zones. Smaller particles such as leaves, needles and twigs are also in the near zone of the larger canopy components. These smaller canopy constituents may be present in large numbers and can have a significant effect on the electromagnetic properties of the medium both as isolated scatterers and through multiple scattering interactions. In addition, it is to be expected that scatterers with dimensions large compared with the dimensions of the medium will be illuminated with a fairly high

degree of amplitude, phase and polarization non-uniformity as a result of the finite coherence distance for the mean field in this type of random medium. Heterogeneous forest canopies consisting of predominantly vertical, cylindrical trunks and large numbers of diverse smaller scatterers have traditionally been divided into two main regions for modeling purposes. The upper region in which there is usually a high density distribution of the smaller particles is referred to as the crown layer, while the lower portion of the canopy is referred to as the trunk layer. In general, RT modeling of forests is constrained by the constructs of RT theory to treat the trunk and crown regions of the canopy as separate layers. This is not an accurate representation of most canopies of this sort which have larger and smaller scatterers co-existing to varying extents within the same volume. Some attempt has been made to overlap the crown and trunk layers in RT modeling by creating a two layer medium in which the lower story contains only trunks while the upper story contains both trunks and smaller scatterers [43]. In this chapter a Monte Carlo simulation is constructed for a heterogeneous twocomponent random medium consisting of large vertical cylinders and small spheres above a ground plane. The second-order interactions between the cylinders and between spheres and cylinders have been taken into account in formulating the scattering model, and all scattering terms used in the model have been validated using the method of moments technique. The backscattering properties of the medium have been obtained using the simulation and the results are compared with those of the corresponding radiative transfer model. Second-order interaction effects in the medium are examined, and the validity of the layered RT model is also investigated.

5.2 Construction of Scattering Terms The heterogeneous two-component medium under consideration consists of vertical dielectric cylinders and small perfectly conducting spheres. Dielectric cylinders are employed to represent the strongly scattering stalks and trunks that exist in a typical vegetation canopy, and the spheres are representative of the smaller scatterers that comprise the upper (crown) layer in such media. Spheres are used in this study because their scattering matrix is independent of orientation, and one purpose of this study is to determine what effect a distribution of such scatterers might have on the overall response of the medium without regard to their orientational characteristics. The spheres are chosen to be perfectly conducting because this reduces the complexity of their scattering matrix and any such reduction will translate into less overall computation time in the Monte Carlo simulations. Since spheres have relatively low scattering cross-sections, it is also desirable to boost their albedo by making them perfectly conducting in order to bring them on par with typical small scatterers that might be found in the canopy. The electromagnetic interaction between adjacent cylinders and between cylinders and spheres is computed to second-order, while the interaction between spheres is ignored. It is assumed that the coupling between spheres is significantly small in comparison with the other terms so as to be negligible. Numerical studies indicate that this should certainly be the case for the number densities of small particles considered in this chapter. In any event, the cylinders are by far the stronger scatterers in the medium and, in most circumstances, are the dominant component of the return. The cylinder-cylinder interaction has been discussed in detail in Chapter 3 and will not be considered here. It suffices to say that the scattering matrix for this interaction has been derived using the plane-wave expansion technique and is appli

96 Z Y X Figure 5.1: Geometry and coordinates for the sphere-cylinder interaction. cable to cylinders that are electromagnetically close with respect to their longitudinal dimension and, at the same time, are in each other's far-zone with respect to the transverse dimension. The sphere-cylinder interaction has been derived using an analytical technique based on the reciprocity principle. The details of this novel technique are presented in Chapter 4 and will not be duplicated here. However, the application of these results to the problem at hand requires extension due to the presence of a dielectric half-space. In this section the results for the sphere-cylinder interaction in the presence of a smooth dielectric ground are presented. 5.2.1 The Sphere-Cylinder Interaction

Figure 5.1 shows the coordinate system for the sphere-cylinder interaction. It has been assumed that the cylinders are electrically long. That is to say, their length is much greater than the excitation wavelength. It has also been assumed that the cylinder is in the far-field region of the sphere, but the sphere is in the near-field region of the cylinder with respect to the cylinder's longitudinal dimension and in the far-field region of the cylinder's transverse dimension. If we consider the analogy with a tree canopy, it is reasonable to assume that the trunks are very long electrically at microwave frequencies, therefore other trunks, branches and leaves will fall within the near-zone of their longitudinal dimension. This is also the rationale for considering adjacent cylinders to be within each other's near-zone. If the radii of the cylinder and sphere are denoted by ac and a, respectively, and the cylinder length is represented by L, then the conditions previously specified may be stated mathematically as 2a2 p > L2 2a2 - > > -2a (5.1) A, A As has been mentioned, it can be shown that the secondary scattered field from a sphere-cylinder pair may be obtained by application of the reciprocity theorem. The expression for the scattering matrix of the pair in the absence of the ground-plane is given in terms of a sum of scattering matrices for the individual particles as if they were isolated plus the second-order interaction terms. This may be expressed as stot = sdc + 8ds + ssc + Scs In the above equation,'dc' refers to the direct wave from the isolated cylinder,'ds' refers to direct scattering from the sphere,'sc' refers to the secondary scattered field from the cylinder when illuminated by the primary scattered field from the sphere

98 and'cs' refers to the secondary scattered field from the sphere when illuminated by the primary scattered field from the cylinder. The scattering matrices for the direct terms are easily found and are not reproduced here. The scattering matrix for the cylinder-sphere (cs) interaction is given by Scs(k1i, k) = -e-ik(cosoizk-,r ) H(1) (ko sin 0Oi) Sv~(-ks, s) * Ev(ki,,-) Ssv(-ks, ). Ec (ki(,-,) with ki = sini cos x + sin;sin iy + cos Oiz ks = sin O, cos Os$ + sin Os sin qy + cos 8,z Fr =- sin Oi cos - sin Oi sin y - cos Oiz and S,v(-k5, I) = 2 [ s S1() xs x (kx )- (k. x ) (o ) S2(]0) k. x rs] ko Ik, x F2 (5.3) Ssh(-k, ri) = 1 [h * f S1(O') r. x (k x F)-h, (k, x F') S2(OI) k x i] ko I kc, x r~ 12'X _ 12 In the last two expressions we have =ucos O cos Ocos+ + cos O~ sin qy - sin 0,z h, = sin q,$-cos qg

99 and +oo 2n + 1 ( Sl(O;): -i n(-l)' (,~ + ~).a. + {L na P(+)(cos8)~ +- 0 P=1-on(ns )+ ) PI( cos~in ) S2(9) = i n+(-1 r ) P2()( COS ) +oo E2(i-~') -C C: Mr(i) eAm(ngiZ < -P (cos0o)-+!3sin +sin) E(i,-T) - cTE(i) eim(-i) (_ sin O + cos y) n=9-oo where is the associated Legendre polynomial,, and m are given in [64]. The scatter symbol, E, given in (5.2) is the scattered field ffor the sphere-cylinder (sc) interaction is given by $'~(ci, k,) - e-iko ~S~@~J H(1) (ko sin 0,p). [rlsv&(ki, k')-r& (k, k')] [r.,(ki, k') - rs,(k, k')] | -[r25v(k,, k') + rS~(k,, k')] -[rShh(k,, k') + rSvh(k,, k')] for which +0000 m= -00 +oo h2- Y~ (-_)~cTr(w _,) eim(g-",) +oo m= -00 m=) = e-oo

100 where the Cm's are as defined previously. The Spq's (p and q equal v or h) are the scattering matrix elements for a sphere with an incident wave-vector kI and a scattered wave-vector k,, and k' = sin 0, cos(O + 7r)i + sin 0, sin(q + ir)y + cos 0~z 5.2.2 Effect of Ground Plane To obtain the solution for a sphere-cylinder pair above a smooth dielectric halfspace we use exact image theory [50, 67]. The form of the Green's function for the problem is approximated by assuming that the image of a source point is the mirror image of that point on the opposite side of the half-space interface and modified by the appropriate reflection coefficient. This approximation is very accurate when the source or field point is not close to the interface and becomes exact when the ground plane is perfectly conducting. When this approach is applied, ten interaction terms are obtained to account for the second-order scattered field from the sphere-cylinder pair above the ground plane. The terms are as follows: i cylinder-sphere-ground ii ground-cylinder-sphere iii cylinder-ground-sphere iv ground-cylinder-sphere-ground v cylinder-ground-sphere-ground vi sphere-cylinder-ground vii ground-sphere-cylinder

101 z -,, Sphere,,I L AO c / /,,,tI Image Excitationit viii sp here- ground- cylinI / Ground Planehere- cylinder-ground t erac' -Lri ion/. viii sphere-ground-cylinder ix ground-sphere-ground-cylinder x ground-sphere-cylinder-ground scattering approximation is valid.

102 In the discussion that follows, we shall obtain the scattering matrices for the cylinder-ground-sphere and ground-sphere-ground-cylinder interactions. The expressions for the remaining terms can be easily derived following the same approach and will not be presented here. Figure 5.2 shows the image geometry for the cylinderground-sphere interaction. In this case, the image cylinder is illuminated by the image excitation, kc. The scattered wave from the image cylinder then illuminates the sphere as shown in the figure. The expression for the image excitation is given by ki = sin 0i cos ixx + sin Oi sin iy' - cos;z The figure also indicates that, under the prevailing assumptions, such an interaction will exist only if the sphere lies somewhere in the specular scattering cone of the cylinder. It has been assumed that because the cylinders are quite long, most of the energy is scattered into the conical region about the direction of incidence. The mathematical statement of the condition for existence of the cylinder-ground-sphere interaction is -L < i + p tan(7r/2 - Oi) < 0 (5.5) For the image excitation, the expressions for the scattered field from the image cylinder are E"(Vc, r V- - I - hl Ec(k, -C) -2 h -r vF (5.6) for vertically and horizontally polarized fields respectively. The polarization vectors, VC and ht, are defined by

103 h x _ _= sin = x + cos z y sin Oi = ~, x = - - cos Oicoos i - cos Oi sin q - sin zO and', = -sin i cos - sin Oi sin ~ y + cos Oi z Modifying the appropriate terms by their corresponding reflection coefficients, R,(0) and Rh(O), we obtain Ec(kr, -sc) rF Rv,(i) C - RFh(i) h' EhrC(kI -Fc) = -r2 Rh(Oi) h' - r RIv(i) V- (5.7) which is the field scattered from the cylinder, reflected from the ground, and incident on the sphere. To obtain the proper far-zone scattering amplitudes for the sphere, the equations (5.3) are modified to account for the image excitation. When these modifications are made, the expression (5.2) becomes SCs9(ki, k) _-e-ikO~(CO~sGz+ks r~) H(1) (k0 sin Oilp) [ S8(-kS, ~) ECvc(k, -JCc) SSV(-kS TJC) j E1hC(kjC ) S~h(-l, r PS C) cr t Cr S (A c

104 z A Cylinder Ground Plane Image Sphere Figure 5.3: Image theory geometry for construction of the ground-sphere-groundcylinder interaction. The ground-sphere-ground-cylinder interaction geometry is shown in Figure 5.3. In terms of the image theory approach, the image sphere is illuminated by the primary plane-wave excitation after modification by the ground reflectivity. The secondary scattered wave from the sphere then illuminates the cylinder after being modified by the ground reflectivity a second time. The condition for existence of this second-order term under the specular cone approximation will be 0 < - o[ + 3tan(08 - r/2)] <L (5.9) If the igroumande sphere and the cylinder were confi geomured as shown in the figure and no ground plane was present, the scattering matrix would have the same form as (5.2). However, the presence of the ground plane alters the incident wave through multipli

105 cation by the surface reflection coefficient corresponding to the transmit polarization state for a particular scattering matrix element. The fraction of the modified incident wave scattered by the sphere that is not depolarized is again modified by the ground reflection coefficient for the incident wave polarization. This wave scatters from the cylinder with the polarization specified by the receive state of the scattering matrix element under consideration. The fraction of the modified incident wave depolarized by the sphere is multiplied by the reflection coefficient having the opposite polarization state from that of the incident wave. This wave is reflected from the cylinder and is either depolarized or not depending on the final polarization of the receive state. The scattering matrix elements for the ground-sphere-ground-cylinder interaction are then given by Sg9gc(ki, k ) = eiko(cos~~z~ Ho1) (ko sin OP) (5.10) [nl&sv(ki, k')- Q2&Sh(ki, k')]Rv(O,) [QlSvh(ki, k') - ~2Shh(ki, k')]Rh(0i) -[Q3 hv(ki, k') + n14Svv(ki, k')]Rv(Oi) -[Q3Shh(ki, kI) + 14S(vh(k,, k')]Rh(Oi) J where Q1 = r,(-ks) Rv(O) 2 = r(-k) Rh (0os) Q3 = r2(-ks)Rh(s3) Q4= r(-ks) R(o,) and all other symbols are as previously defined.

106 5.3 Validation of the Scattering Terms In order to construct the Monte Carlo simulation of the two component medium described previously, it is first necessary to validate the scattering model. In this section it is demonstrated that all the fundamental pieces of the model are correct. The Numerical Electromagnetics Code (NEC) [6], which is a software package based on the method of moments (MoM), has been chosen for this purpose because of its flexibility and accuracy. As discussed previously, the two component medium being simulated consists of a mixture of vertical cylinders and perfectly conducting spheres above a perfectly conducting ground plane. In this model it has been considered that the spheres do not themselves interact electromagnetically,.however the interactions between the cylinders and also between spheres and cylinders have been taken into account. It has been attempted to validate each of these interactions using the MoM code. For this validation work a perfectly conducting sphere, with ka=1.7, has been employed as well as two different types of cylinder models. The sphere model was composed of variably segmented, perfectly conducting rectangular patches as described in [7]. One of the cylinders used was 18 cm in length, having a diameter of 0.1 cm and a conductivity of 100 mhos/m. The other cylinder was also 18 cm in length but had a diameter of 0.55 cm and was composed of perfectly conducting patches. At the operational frequency of 9.25 GHz chosen for this study, the cylinders are between 5 and 6 wavelengths long which is sufficient length to provide a valid realization of the specular cone approximation. Thin cylinders were chosen so that the extended thin wire kernel and the finite conductivity feature provided in NEC for wire models could be used. In this case only 60 unknowns were necessary to achieve computational accuracies of 1 percent or

107 better. The thin cylinders were chosen to be of finite conductivity for two reasons. First, finite length perfectly conducting cylinders support strong axial travelling waves that are not accounted for by our analytical model. These travelling waves have a significant effect on the vv and vh-polarized scattering response of the cylinder. This effect is much less pronounced for hh polarized scattering. Dielectric and finitely conducting cylinders with lengths of greater than several wavelengths do not strongly support axial travelling waves. Our two-component medium is supposed to simulate effects that might be observed in vegetation canopies, and it is not anticipated that travelling waves play a significant role in scattering from branches and stems. The second reason for choosing finitely conducting cylinders is that the finite imaginary part of their effective dielectric constants requires the use of the same complete set of coefficients as is used in computing the scattering from dielectric cylinders, thus allowing us to thoroughly test these routines. MoM data for this extended thin wire model has an associated uncertainty of approximately ~0.25 dB. It was necessary to use the larger diameter cylinders to obtain useful results for hhpolarized scattering. Larger diameter cylinders require the use of perfectly conducting patches because the extended thin wire approximation used in NEC is not valid for structures of this size, and the version of the code presently available only allows for perfectly conducting patches. These cylinder models were constructed with in excess of 800 patches and required a fairly large amount of computation time on an IBM RS6000 workstation. The overall uncertainty in the MoM data for these cylinders is about ~0.5 dB or roughly twice that for the thin cylinders. For the scattering patterns used in this validation, the angle of elevation and the azimuthal angle is defined in Figure 5.4. The plane of incidence is the x-z plane and the azimuthal incidence angle is 180 degrees. The cylinder is always located at the

108 Z Elevation Angle A k LA!! I I Y Azimuthal Angle X Figure 5.4: Scattering geometry and reference angles. origin and the relative cylindrical coordinates of the sphere are presented. Figure 5.5 provides a comparison between the backscatter elevation patterns for the sphere single scattering model used in the simulation and the numerical computation. It is seen that the model is in agreement with the MoM computation to within better than ~0.2 dB over the entire angular range considered. In this case, the sphere is located 9 cm above the ground. It has been found that the sphere scattering model is good for such spheres to within about 3d2/A of the ground plane, where d is the diameter of the sphere. This indicates that the single scattering model for the spheres should be good for number densities providing average separations on this same order of magnitude. These number densities are well within the bounds employed in the Monte Carlo simulations.

109 -20.0 -25.0 -30.0 Ei -35.0 -40.0 - Theoretical (VV) -45.0 ------— Theoretical (HH) A MoM Computation (VV) -50.0 v MoM Computation (HH) -55.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.5: Backscatter pattern for a perfectly conducting sphere above a conducting ground plane. The sphere has ka = 1.69, where a is the sphere radius, and is located 9cm above the ground plane. The frequency is 9.25 GHz. -10.0 -20.0 - -30.0 -40.0::.- Scod-OrdrThry -...... First-Order Theory -50.0 A MoM Computation -60.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.6: Backscatter pattern for a pair of vertical cylinders on a conducting ground plane. The cylinders are 18cm in length, O.1cm in diameter and separated by 2cm in the endfire configuration. The Conductivity of the cylinders is 100 S/m. vv-polarized return at 9.25 GHz.

110 Figure 5.6 shows typical backscatter results for a pair of thin vertical cylinders above a conducting ground plane. The cylinders are separated by 2 cm and are positioned so as to lie along the direction of incidence (endfire configuration). A complete experimental verification of the second-order cylinder-cylinder interaction terms above a conducting ground may be found in [69]. Figures 5.7 through 5.12 illustrate that the agreement between the sphere-cylinder interaction model described in the previous section and the MoM computation is quite good. Elevation and azimuth scattering patterns have been provided for all polarizations. The elevation patterns shown are for the backscatter direction except in the case of the cross-polarized return for which the pattern was computed at a bistatic azimuthal angle of 350 degrees to provide the best comparison between the first and second-order results. The co-polarized second-order results agree with the MoM computation to within about ~0.5 dB over most of the angular range, however, the hh-polarized MoM results have a somewhat higher level of inaccuracy due to the details of the larger diameter cylinder model. The cross-polarized second-order results have an overall agreement with the numerical data to within ~1 dB or so for the elevation pattern and ~0.5 dB for the azimuth pattern. 5.4 Heterogeneous Canopy Simulation Having validated the scattering terms for the interaction of spheres and cylinders above a ground-plane, we will now attempt to obtain the backscattering properties of collections of such objects that comprise the heterogeneous two-component canopy. For a given arrangement consisting of many cylinders and spheres, the solution of the scattering problem can be obtained to second-order by computing the single and pairwise interactions for every particle in the ensemble, except that in this case

111 -15.0.... -20.0 -25.0 ~~~~~~-25.0 Socond-Order Theory -------- First-Orda Theory a MoM Canomputatiol -30.0.... 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.7: Backscatter pattern for a sphere-cylinder pair above a perfectly conducting ground plane with p=1.5cm, 6 = 1800 and z=9cm. vv-polarized return at 9.25GHz -15.0................, - ---- -20.0 -25.0 ----- Scaond-Order Theory......... Frst-Order Theory A MoM Coanputation -30.0..... 0.0 90.0 180.0 270.0 360.0 Azimuth Angle (deg) Figure 5.8: Bistatic scattering pattern for a sphere-cylinder pair above a perfectly conducting ground plane with p=1.Scm, rq = 180~ and z=9cm. vvpolarized return at 9.25GHz. The elevation angle is 37~.

112 -35.0........,.......... -45.0 -65.0 - Second-Order Theory -------- First-Order Theory A MoM Computation -75-0..... I.... I.... I.... I.... 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.9: Bistatic scattering pattern for a sphere-cylinder pair above a perfectly conducting ground plane with 1=2.0cm, q = 45~ and i=18cm. The azimuthal scattering angle is 350~. vh-polarized return at 9.25 GHz. -25.0... -35.0 - -45.0 -5.50 Second-Ordr Thy -------- First-Order Theory, MoM Comnputation -65.0 0.0 90.0 10.0 270.0 360.0 Azimuth Angle (deg) Figure 5.10: Bistatic scattering pattern for a sphere-cylinder pair above a perfectly conducting ground plane with 1=1.5cm, q = 45~ and, 18cm. The elevation angle is 52~. vh-polarized return at 9.25 GHz.

113 -12.0 ro*-17.0 M -22.0 -27.0 \,'- Second-Order Theory -------- Frst-Order Theoy \ MoM Computation -32.0.... I, A.. 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.11: Backscatter pattern for a sphere-cylinder pair above a perfectly conducting ground plane with p=2.0cm, b = 450 and ~=9cm. hh-polarized return at 9.25 GHz. -15.0 -20.0 -25.0 -. MoM Compuution -40.0 0.0 90.0 180.0 270.0 360.0 Azimuth Angle (deg) Figure 5.12: Bistatic scattering pattern for sphere-cylinder pair above a perfectly conducting ground plane with p=2.0cm, q = 45~ and z=9cm. The elevation angle is 35~. hh-polarized return at 9.25 GHz.

114 the sphere-sphere interactions have been excluded. The statistical properties of this random medium are simulated by application of the Monte Carlo method. The principle of the Monte Carlo simulation based on second-order interactions is as follows: 1 An ensemble of randomly positioned cylinders and spheres is generated using a random number generator. In this case the cylinder positions are uniformly distributed within a circular area, and the sphere positions are uniformly distributed within a cylindrical volume. The number of cylinders and spheres used is dependent on the specified numbers per unit area or volume and the dimensions of the cylindrical region that constitutes the medium. 2 The scattering is computed for all cylinders and spheres and between all pairs of particles within the ensemble up to second-order excluding the interactions between spheres. 3 The ensemble is re-randomized and the scattering recomputed as discussed above. The number of independent samples is chosen so as to make the variance as small as possible within limits depending on the computing time. For the cases analyzed in this article the sample number was chosen to be two hundred except in the highest density case for which the sample number was one hundred. 4 The values of the scattering coefficients (av, cahh) are found from the ensemble average. The same is true for the co-polarized phase difference (, and the degree of co-polarized phase correlation a [66]. The simulations examined in this article utilized perfectly conducting spheres

115 having diameters of 0.635 cm and characterized by number densities of 14,147 per cubic meter. The spheres were generated inside a cylindrical volume 0.1 meter in height and 0.6 meter in diameter. Identical cylinders 18 cm in height and 0.55 cm in diameter with a dielectric constant cr of 35+ill were uniformly distributed in the same 0.6 meter diameter circle used to generate the cylindrical volume for the spheres, but the cylinder bases were constrained to rest on the ground plane whereas the sphere layer could begin at any height above the plane. Simulations were performed for spheres alone, cylinders alone and combinations of cylinders and spheres together. Two basic heterogeneous canopy configurations were examined. In one case the sphere layer overlapped the upper 10 cm of the cylinder layer, while in the other case the sphere layer was distinct and existed separately above the cylinder layer. These configurations are shown in Figure 5.13. Figures 5.14 and 5.15 illustrate typical convergence properties of the simulations for the scattering cross section, ao, and the copolarized degree of correlation, a. In all cases tested, the simulations converged to within about 0.5 dB for cr~ and 5% for ac in one hundred samples or less. The operational frequency is 9.25 GHz for all the examples considered in the following discussion. 5.4.1 Comparison with Radiative Transfer Electromagnetic modeling of heterogeneous tree canopies has traditionally been accomplished using two-layer radiative transport (RT) theory. In this type of model the trunk and crown layers are treated as being either separate [96] or, if it is desired to have part of the crown layer in the same region as the trunk layer, a mixed two-layer geometry is formulated [43]. The mixed layer formulation treats the lower portion of the trunks as if they were a separate layer, while the crown layer of the canopy consists

116 0 000 0 0 00 00 0O Crown Spheres 0 O 0 0 0 000000 0 0 0 0 000 O Layer ~ 00 0 00 0 0 O O O 0 0 0 Trunk Cylinders 1 I I I Layer Ground Layered Canopy / Mixed Canopy els. 0 0 00 00 0 0 0 0 00 O 00 Crown -0 0 O 0 0 0 Trunks Ground els. RT theory. layer consisting of 14,147 spheres per cubic meter. In this case, the Monte Carlo that the RT results or a layer of firstonde saterr houdb bu el

117 10.0.,. 9.0 8.0 7.0 6.0 _ Second-Order Simulation -------- Fst-Order Simulation 5.0 ~ I..., I.. 0. 50. 100. 150. 200. Number of Samples Figure 5.14: Typical convergence behavior for the Monte Carlo simulation of the backscattering coefficient. The number density is 106.1 cylinders per square meter. vv-polarized return at 9.25 GHz. 1.0 0.8 0.6 0.4 Second-Order Simulation 0.2 0.0... 0. 50. 100. 150. 200. Number of Samples Figure5.15: Typical convergence behavior for the Monte Carlo simulation of the phase statistic, ac. The number density is 106.1 cylinders per square meter. The frequency is 9.25 GHz.

118 10.0 Radiative Transfer Monte Carlo Simulation 0.0 0,00 _ -10.0 -20.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.16: Homogeneous layer of 14,147 conducting spheres per cubic meter distributed in 10cm layer above conducting ground. Spheres have ka = 0.62. vv-polarized backscatter cross section. incoherent addition of power while the Monte Carlo simulation is based on coherent addition of the fields. However, for a collection of small discrete scatterers whose phase centers are distributed randomly within a volume above a smooth ground, the net effect of cross coherence is not significant over most of the angular range. The cross coherence terms do produce a minor but noticeable effect in certain sub-regions of the angular backscatter spectrum as can be seen in the Figure. Figures 5.17 and 5.18 show the typical response of a homogeneous layer of cylinders above a perfectly conducting ground plane. The number density for this layer is 106.1 cylinders per square meter which is sufficient to produce a significant amount of cylinder-cylinder coupling for vv polarization and somewhat less coupling for hh polarization. The coupling between cylinders also produces a fair amount of crosspol which is not provided by first-order RT in backscatter. At low angles of elevation (near nadir) the RT result for ai, is about 4dB below those for the simulation. This

119 10.0 o.o - "' -10.0' -— E- Second-Ordr Monte Carlo — s —- Fiust-OrdcrMonte Carlo...... Fust-Order R.T. -20.0.... 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.17: Homogeneous layer of 106.1 vertical dielectric cylinders per square meter above conducting ground. Cylinders are 5.6 wavelengths long and 0.17 wavelengths in diameter. Relative dielectric e, = 35 + ill. vv-polarized backscatter cross section. 10.0 9.0 0,, -10.0 -'. Second-Order Moie Carlo ---- First-Order Monte Carlo -..... Radiative Transfer -20.0...... 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.18: Homogeneous layer of 106.1 vertical dielectric cylinders per square meter above conducting ground. Cylinders are 5.6 wavelengths long and 0.17 wavelengths in diameter. Relative dielectric Er = 35 + ill. hh-polarized backscatter cross section.

120 difference increases with elevation angle and becomes quite large at angles near 800. This effect is partly due to overestimation of the extinction matrix in the RT formulation for large scatterers. The RT model computes extinction based on the dimensions of a scatterer. In the case of a layer of cylindrical scatterers, the extinction is related to the length of the cylinders. However, as multiple scattering becomes more significant at higher angles of incidence, the correlation length for the coherent field in the medium decreases. Thus, the effective scattering length of cylinders in a dense medium is actually shorter than the physical length of the cylinders, making the true extinction less than what RT theory would predict. In the plot of auh, we observe a 3dB difference to be present at low angles of incidence. This effect is due both to the significant cross coherence present in the simulation and also to the integration of the source functions in RT theory which is not actually valid for layers of extended scatterers. In the case of hh polarization, the extinction is much less than that for vv polarization due to the relatively small diameter of the cylinders. In this regime, the aforementioned anomalous behavior resulting from integration of the source functions drives the RT angular response. Figures 5.19 through 5.22 show the result of adding a volume consisting of 14,147 spheres per cubic meter to the canopy. The sphere layer begins at the top of the cylinder layer and extends for 10cm above it. In Figure 5.19 the cylinder layer is characterized by a number density of 35.4 cylinders per square meter. At this number density, neither the mutual coupling between cylinders nor the coupling between cylinders and spheres produces much second-order interaction in the simulation results. The RT result is seen to have a fairly constant offset from the simulation since extinction in both layers is minimal. Figure 5.20 is for the same canopy configuration as in Figure 5.19 except the cylinders have a number density of 106.1 cylinders

121 10.0 0.0 -10.0 -*- _ Second-Ordaer Monte Carlo — b —- Firat-OrderMonte Carlo - - - -...... First-Order R.T. -20.0 I., I, I I. I. I -20.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.19: Heterogeneous canopy consisting of 35.4 cylinders per square meter in lower layer and 14,147 spheres per cubic meter in an upper layer 10cm thick. The cylinders are 18cm high and 0.55cm in diameter. The spheres have ka = 0.62. vv-polarized backscatter coefficient at 9.25 GHz. per square meter. Once again, at this cylinder density level, the coupling between cylinders is not insignificant. It is evident from comparison with Figure 5.17 that the addition of the sphere layer does not change the vv response of the Monte Carlo simulation appreciably since the level of scattering from the sphere layer is approximately 10dB down from that of the cylinder layer. This can be seen more clearly in Figure 5.21 which gives a direct comparison between the simulation results for the 106.1 cylinders per square meter canopy with and without the sphere layer present. The small difference at angles of incidence closer to vertical is present because the cylinder response is significantly lower in that region. Comparison of Figures 5.22 and 5.18 shows that the addition of the sphere layer produces the same effect at lower angles of incidence for the hh-polarized response. At higher angles of incidence, extinction in the sphere layer reduces the hh backscattering response of the two layer canopy over that of the cylinder canopy. The most notable difference between the

122 10.0 ~, I, ~ 0.0 0, -10.0 _ —— ~ Second-Order Monte Carlo — = —- First-Order Monte Carlo - - -..... First-Order R.T. -20.0 I. I. I.. I. 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.20: Heterogeneous canopy consisting of 106.1 cylinders per square meter in lower layer and 14,147 spheres per cubic meter in an upper layer 10cm thick. The cylinders are 18cm high and 0.55cm in diameter. The spheres have ka = 0.62. vv-polarized backscatter coefficient at 9.25 GHz. 10.0 0.0' -10.0 Homnogencous Canopy * —-—.. Heterogeneous Canopy -20.0 ~ I. I I I I I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.21: Comparison of Monte Carlo simulation of vv-polarized backscatter coefficient for a homogeneous canopy and a heterogeneous canopy. The homogeneous canopy has 106.1 cylinders per square meter. The heterogeneous canopy has the same number density of cylinders but has 14,147 spheres per cubic meter in a 10cm layer on top.

123 10.0 -10.0 - _ Second-Order Monte Carlo — a —- First-Order Monte Carlo... —. First-Order R.T. -20.0 I I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.22: Heterogeneous canopy consisting of 106.1 cylinders per square meter in lower layer and 14,147 spheres per cubic meter in an upper layer 10cm thick. The cylinders are 18cm high and 0.55cm in diameter. The spheres have ka = 0.62. hh-polarized backscatter coefficient at 9.25 GHz. homogeneous and heterogeneous canopies can be seen in Figures 5.23 through 5.25. Figure 5.23 illustrates the dramatic increase in depolarization that results solely from coupling between the sphere and cylinder layers, and Figures 5.24 and 5.25 show that addition of the sphere layer markedly changes the behavior of the co-polarized degree of phase correlation, ca. This is not surprising since this phase statistic is a sensitive measure of the level of multiple scattering in the canopy. It should be noted that both the first-order Monte Carlo simulation and first-order RT theory yield a result for ac of unity for a homogeneous layer of vertical cylinders. Finally, Figures 5.26 and 5.27 show the effect of using a mixed two-layer canopy model versus a pure two-layer canopy model as was illustrated in Figure 5.13. In this example, the pure two-layer canopy has a 10cm layer of small conducting spheres beginning at the top of the cylinder layer, while in the mixed canopy the sphere layer begins at 8cm above the ground plane and extends to the top of the cylinder

124 -20.0 ~'" " " [i..., -.- -E -30.0-40.0 0 -50.0 -0- Hmnogencows Canopy -60.0. ---- Heterogeneous Canopy -70.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.23: Comparison of Monte Carlo simulation of vh-polarized backscatter coefficient for a homogeneous canopy and a heterogeneous canopy. The homogeneous canopy has 106.1 cylinders per square meter. The heterogeneous canopy has the same number density of cylinders but has 14,147 spheres per cubic meter in a 10cm layer on top. 1.0 0.8 0.6 0.4 -=- Secmd-Order Monte Carlo 0.2 - ---.... First-Order Monte Carlo ------ Radiative Transfer 0.0 I I I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.24: Homogeneous canopy consisting of 106.1 cylinders per square meter. The degree of copolarized phase correlation at 9.25 GHz.

125 1.0 --- 0.8 Carlo Radiative Transfer 0.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.25: Effect on the degree of copolarized phase correlation of adding a layer of 14,147 small conducting spheres per cubic meter in a layer 10cm thick on top of the canopy of 106.1 cylinders per square meter. 10.0 0.0 - 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.26: Comparison of pure and mixed layer R.T..models for a canopy consisting of 106.1 cylinders per square meter and 14,147 spheres per cubic meter. The crown layer is 10cm in height in both cases. vv-polarized return at 9.25 GHz. The pure and mixed copolarized Monte Carlo simulations are almost identical.

126 -20.0 -30.0 - -40.0 0 -50.0 Pure Layer Canopy -60.0 I-".... Mixed Layer Canopy -70.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 5.27: Comparison of pure and mixed layer Monte.Carlo simulations for a canopy consisting of 106.1 cylinders per square meter and 14,147 spheres per cubic meter. The crown layer is 10cm in height in both cases. vhpolarized return at 9.25 GHz. layer. Figure 5.26 shows that using the mixed two-layer model actually lowers the RT co-polarized backscatter response severely over that of the pure two-layer model. This decrease occurs because the phase matrix of the cylinders is proportional to the square of the cylinder length, and this length is divided into smaller subsections in the mixed layer model. For example, division of the trunks into two equal portions would reduce the phase matrix by a factor of 3dB. There is almost no difference between the co-polarized backscatter response of the Monte Carlo simulations for the two different canopy types. Figure 5.27 shows that the pure and mixed layer canopy simulations do have significantly different crosspol responses. This is because the spheres are, on the average, closer to the cylinders in the mixed layer model and this elevates the degree of coupling present in the medium.

127 5.5 Conclusions In this chapter a Monte Carlo simulation has been constructed for a heterogeneous two-component medium consisting of large vertical cylinders and small spheres above a ground plane. The second-order interactions between cylinders and between spheres and cylinders have been included in the scattering model, and all scattering terms have been validated using the method of moments. The results of the Monte Carlo simulation have been compared with those of corresponding radiative transfer models for similar media. It has been found that the radiative transfer models do not properly predict the scattering behavior of media which are composed of scatterers having dimensions large compared to the dimensions of the medium. This is a result of the fact that the extinction matrix for the medium as computed by radiative transfer is overestimated for large particles in any plane of polarization corresponding to a large dimension and produces excessive attenuation of the scattered wave for that polarization especially at higher angles of incidence. In addition, even though the extinction matrix is not incorrectly estimated for polarizations that are not associated with large dimensions of the scatterer, the source function for the medium will still be improperly estimated in that dimension and will not reproduce the proper angular trend for the scattered wave. It has also been determined that the mixed-layer radiative transfer model which attempts to divide the canopy into an upper layer consisting of trunks and smaller scatterers and a lower layer consisting of trunks alone can produce results that are seriously in error. Another conclusion is that while second-order interactions between the cylinders can yield significant changes in the level of backscattered copolarized radiation, the predominant effect of the interaction of the cylinders with the smaller spheres is to produce a change in the angular trend of the cross polarized response. The second

128 order interactions between spheres and cylinders have also been shown to produce a distinct change in the copolarized phase statistics of the backscattered wave.

CHAPTER VI A HYBRID MODEL FOR ELECTROMAGNETIC SCATTERING FROM FOREST CANOPIES 6.1 Introduction In previous chapters Monte Carlo simulations have been developed and validated for vegetation-like random media consisting of long vertical cylinders above a dielectric half-space representing tree trunks or plant stalks and including distributions of smaller spherical particles representing leaves, branches and other components such as may be found in the crown layer of a forest canopy [62]. The interactions between the cylinders and between spheres and cylinders up to second order were included in the simulation algorithm for the purpose of determining the effect of multiple scattering in the canopy. The motivation for developing the Monte Carlo canopy simulations was to provide a benchmark with which to evaluate the performance of the RT model in reproducing the electromagnetic scattering behavior of tree canopies. It has been found that the RT model does not accurately reproduce the scattering behavior of layers of electrically long vertical cylinders above a dielectric half-space. Figure 6.1 shows the discrepancy between the vv-polarized backscatter response of the RT model and the Monte Carlo simulation for a layer consisting of 106 cylinders per 129

130 10.0 0.0' -10.0 -0 —- Second-Order Monte Carlo —, —- First-Order Monte Carlo...... Radiative Transfer -20.0 I I... I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure6.1: The vv-polarized backscattering coefficient of 106.1 vertical dielectric cylinders per square meter above a conducting ground plane at 9.25 GHz. The cylinders are 18cm long and 0.55cm in diameter and have a dielectric constant e, = 35 + i11. square meter above a perfectly conducting ground. The cylinders are 5.5 wavelengths long and 0.17 wavelength in diameter and have a dielectric constant, Or = 35+ill1. The figure illustrates the typical overestimation of the extinction behavior of the canopy trunk layer that is obtained using the RT approach. In addition, there is an offset between the RT results and those of the simulation due to the improper treatment of coherent volume scattering found in RT theory and resulting from integration of the source functions and the inability of RT to account for cross coherence effects. The extinction problem occurs because RT uses the scattering matrix of the long cylinder to compute the extinction matrix of the trunk medium. Figure 6.2 shows the vv-polarized backscatter response of a layer of cylinders identical in all respects to those of the previous figure except that the layer has a number density of 35 cylinders per square meter. In this case, the medium is fairly sparse with all cylinders being

131 10.0 0.0 -10.0 --- Second-Order Monte Carlo.' ----- First-Order Monte Carlo...... First-Order R.T. -20.0. I I. 1 I I. I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.2: The vv-polarized backscattering coefficient of 35.4 vertical dielectric cylinders per square meter above a conducting ground plane at 9.25 GHz. The cylinders are 18cm long and 0.55cm in diameter and have a dielectric constant c, = 35 + ill. uniformly illuminated, and the RT extinction approximation is valid. However, when the population of cylinders becomes larger, interactions between the scatterers cause decorrelation of the coherent field as it passes through the medium. This decorrelation produces non-uniform illumination of the cylinders and makes them scatter as if they were shorter than their actual physical length. The end result is to cause the RT model to overestimate both the extinction matrix and the phase matrix for the layer. In contrast to their distinct behaviors for a layer of vertical trunks, the RT model and Monte Carlo simulation agree well for a medium consisting of small scatterers distributed throughout a volume above a dielectric half-space. Figure 6.3 shows the comparison between the vv-polarized backscatter response of the RT model and the Monte Carlo simulation for a layer of small spheres 10 cm deep above a perfectly conducting ground. The spheres have a diameter of about 0.2 wavelength and are

132 10.0 Radiative Transfer 0 Monte Carlo Simulation 0.0 -10.0 -20.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.3: The vv-polarized backscattering coefficient of 14,147 metal spheres per cubic meter in a layer 10cm high above a conducting ground plane at 9.25 GHz. The spheres have a diameter of 0.6350cm. perfectly conducting. The density is 14,147 spheres per cubic meter. In this case there is no offset between the methods because the cross coherence terms are small over most of the angular range, becoming slightly larger at elevation angles approaching horizontal. The extinction is properly estimated by RT since all the particles are small and uniformly illuminated by the incident plane-wave and, since the eigenvalues of the extinction matrix are small in this case, integration of the source functions does not produce serious errors. The detailed second-order Monte Carlo simulation is quite computationally intensive and requires a fairly large amount of CPU time as compared with the RT model. However, such second-order Monte Carlo simulations have been shown to give accurate estimates of scattering cross-section and phase statistics for layers of vertical scatterers [69]. First order simulations give excellent estimates of the co-polarized scattering cross-section for cylinder layers while at the same time being computation

133 z ki position of nth cylinder Rn Figure 6.4: Geometry of the cylinder layer for transmission matrix derivation. ally efficient. For the crown layer, the RT model is accurate, efficient and easy to formulate, thus providing a useful tool for sensor inversions. In this chapter a hybrid model has been developed that uses radiative transfer to compute the scattering from the canopy crown while employing a first-order Monte Carlo simulation to obtain the scattering behavior of the trunk layer. In section 2 of this chapter an analytical expression for the transmission matrix of the trunk layer is derived which may then be used to provide the proper extinction behavior for the RT crown model. Section 2 also introduces the concept of the effective scattering matrix of an average scatterer in the medium. Section 3 describes the development of the hybrid scattering model algorithm. Section 4 presents evidence to support the conceptual foundations of the hybrid model and gives some numerical results and examples of its application in comparison with RT theory and pure Monte Carlo simulations.

134 6.2 The Transmission Matrix for the Trunk Layer Consider an infinite layer of finite length vertical cylinders whose phase centers are randomly distributed in the xy-plane as shown in Figure 6.4. Suppose that a plane-wave whose direction of propagation is denoted by k- is incident on the layer and is given by: Ei = ei eikk where ei is the polarization vector of the incident field. The expression for the scattered electric field from the nth cylinder is: ikok.p eikolr-P II ES(r)-=ekk, e kkp Sn(Rnj ki) (6.1) Ir-p -'J where Pn = XnX+ YnY r = zz r -Pn Ir-pnl Rn= and ki = cos qi sin Oix + sin qi sin 9$p - cos Ojz is the unit propagation vector for the incident wave. In the expressions given above, n is the location of the phase center of the nth cylinder, r' is the location of the observation point which, in this case, lies along the negative z-axis and Rn is the unit vector specifying the direction between the location of the nth cylinder and the observation point. In equation (6.1), Sd(R-, k.) is the far-field scattering amplitude of the nth cylinder and is given by:

135 Sn [Svv(Rv, ki) ev + Svh(R, k,) e] V(Rn) +[Shv(Rn k + (, ) e + Shh( i)] h(R) (6.2) where ev = i Vi eh = ei * hi and /ki xZ hi = ki x ZI hi x ki Vi = -.11 Ihi x kil jRn xI h.(Rn)= h(R) x R v(Rn) = -I Ihs(Rn)x Rn Spq in equation (6.2) is the effective pq scattering matrix element of the nth cylinder, where p and q may be either v or h. At this point it should be noted that the effective scattering matrix of a particular cylinder is not necessarily the same as the scattering matrix of an isolated cylinder with the same dimensions and electrical properties. If the medium is not sparse, then each cylinder will experience a different local field depending on its location in the ensemble. If the level of multiple scattering in the medium is significant, then any particular cylinder may, in general, be subject to an excitation that is non-uniform in amplitude, polarization and phase due to its interaction with other cylinders. Thus, each cylinder is represented by its effective scattering matrix which

136 reflects the unique current distribution that exists on that cylinder as a result of its individual electromagnetic environment. Therefore, the total scattered field at the observation point is obtained by summing the effect due to all cylinders in the layer and is expressed as oo ikaki eikolr-P.I ^I^ ES(r) = E e= ~Pn Sn(Rn, k,) (6.3) >=l Ir- Pn'Writing (6.3) explicitly in terms of the cylindrical coordinates yields - eikoV/pT2+z2 E"(r) = E eikopnsincos(c- i) Sn(Rn, ki) (6.4) n=1 2 + In (6.4),,n is the cylindrical angular coordinate of the nth scatterer, /,i is the elevation angle of the incident wave as measured from the positive z-axis and Oi is the azimuthal angle of the incident wave in the usual sense. If the observation point is located very far away from the cylinder layer, the subtended angle between cylinders in the layer as seen by the observer is quite small, and it is possible to pass from the sum of equation (6.4) to an area extensive integration over the xy-plane. The expression for the scattered field becomes ES(r) M J eikOPSiniVPcos(-'i)e S(R, ki) pdpdq (6.5) in which M is the number of cylinders per unit area and R = r-p is the unit vector from the observation point to the integration point. Changing variables to p' = p/jzj we have 2r oo eikolzl{p' sinfli cos(0-0i)+V/p+ 121 ES(r) _ M S(R, ki)zlp'dp'do (6.6) Since kolzl > 1 the integral can be evaluated analytically using the stationary phase approximation. The stationary function is given by f(p', qS) = p' sin p/3i cos (i - Xi) + p' + 1 (6.7)

137 Taking the first partial derivatives with respect to p' and q we have of pi sin p, cos (O - O) + Of -p' sin,p i sin (O - 0i) (6.8) Setting the partial derivatives equal to zero provides the stationary phase point at 0 = qi +,r and p' = tan Pi which corresponds to the forward direction of propagation through the layer of cylinders. Using these values in the stationary phase integration formula along with the second partial derivatives, the expression for the scattered mean field becomes 27riM ikOk,.r (E (r)) k cos (S(ki, ki))ei~r (6.9) for which we have taken the mean value of the scattered field and the ensemble average of the far-zone effective scattering amplitude of the cylinder layer after completing the integration of (6.6). Taking the ensemble average of the effective scattering matrix over all cylinders in the medium results in the average effective cylinder scattering matrix for the medium. This might be thought of as the scattering matrix of a typical cylinder in the electromagnetic environment representative of the average medium. Using this average effective scattering matrix in equation (6.9) gives the expression for the mean scattered field from the cylinder layer. To obtain the total coherent field at the observation point, the incident field is added to the expression in (6.9) which yields (E'(r)) - E'(r) + (E"(r)) = [ei+ koc, (S(kiki)) eko (6.10) The transmission matrix T for the coherent field is defined by (Et(r)) = TE' (6.11)

138 The transmission matrix elements are found by first setting the incident field polarization to es = v and taking the dot product with v and h successively. The same procedure is repeated with ei = h and we find 2~riM TVV = 1 + - (svv(ki, ki)) ko cos /3i 27riM Tvh = (Sh(ki, Iki)) Th = (Sh,(ki, ki)) (6.12) ko cos Oi Thh = 1 + 2 M (Shh(ki, ki)) ko cos/3i where (Spq(ki, ki)) is the ensemble average effective cylinder pq scattering matrix element for the forward scattering direction. (p and q may be either v or h) For the purpose of use in combination with radiative transfer theory we will need to develop an expression for the intensity transmission matrix for the coherent field. Let us define the coherent intensity transmission matrix r as follows It = YIi (6.13) where Ii and It are the modified Stokes vectors for the incident and transmitted intensities respectively. They are defined by IEv12 = Eh 12 (6.14) t7 2,Re(EvE'*) 2Qm(E E-*) and (I E 12) It = (Ehl 2) (6.15) 27 2Re(Ef Ej*) 2Sm(Et Ej*)

139 where r7 is the intrinsic impedance of free space. From equation (6.11) we find (Ev) = Tv,,E + TvhE' (6.16) and from which we may derive (I E'2) = ITvI2IE'I2 + 2Ree[(TvvTv'h)E'Eh*] + ITvh 12EhI2 (6.17) Equation (6.13) can be also expressed as (IEl 2) = TIEv 12 + T121LEhL2 + 2T1,Re(EVEEh) + 2T14.m(EvEh ) (6.18) where Tij is the ijth matrix element of r. Equating the two expressions we find rll = ITvvl2 T22 = ITvhl2 T13:= Re(TvvTvh) T14 = -O-m(TvTv'h) The rest of the matrix elements are found in a similar manner, and the expression for the intensity transmission matrix is found to be Y = ITvv 12 ITvh 2 Re(TvvTvh) - am(TvvTvh) IThvl2 IThhl2 1e(ThvTh/h) -Sm(ThvT'h) 23e(TvVThv) 2JRe(TvhThTh) Re(TVvTh'h + TvhThv) -Sm(TVVT'h- TvhT~hv) 2Sm(T/vTr) 2Qm(TvhThh) Q~m(TvvTvh + TvhTj) 3e(TVVTh*h - TvhThv) (6.19)

141 Diffuse Boundary \ z=0O 0 0 0 000 0 O0 0 C rownO +I 0Layer 0 \j 0 0 0 0 0 \ z=d0 0 0 0 0 z = -d Trunk I, 7/ / /Layer z = -d' Specular Ground Figure 6.5: Tree canopy layout for the hybrid model showing intensity designations. where p = cos 0 and K, is the extinction matrix of the crown layer. The explicit expression for ic, can be found in [84] and will not be replicated here. The source functions F+ and Fc account for the scattering of energy that is incident on an elemental volume of the crown medium from all directions into the directions (p, q) and (-/u, I) respectively and are given by F+ (/, q, z) 1[?c(y, z) =,':)I ( (',;', z) do' + jj'P((, O -Pc',;' )I(-', I', z) d'] (6.22) F(- o, z) Pc (-p,;; #'7, Y')I )I' (-, O', z) ddg] + jherc( i, =; -i', 5)1 n(-i', z1, z) dn'] relates the intensity of radiation incident on an elemental volume of the crown from

142 the direction (I', q') to the intensity of radiation scattered by the elemental crown volume into the direction (j, t). Once again, the definition of the phase matrix can be found in [84] and is not reproduced here. The phase matrix of the crown is the sum of the phase matrices of the individual components of which the crown is composed and may be written Pc(7 8, 4,; i, i) = E'Pk(Ls, O.,; i, i) k where?k is the phase matrix of particle species k in the crown. The formal solutions to the radiative transfer equations (6.21) are expressed by I+(l, i, z) = e-'c(z+d)1/ I+(P, X,-d) + L e-C(z-z')// F+ (, X, z') dz' IC (-,, z) = e / IC(-,, 0) + f ec(z-z')/ F-(-p, X, z') dz' (6.23) for the upward-going and downward-going intensities inside the crown. If we consider the boundary condition at the diffuse boundary interface between the crown and trunk layers we have I+(p, 4, -d) = IT (H., -d) (6.24) where It is the upward-going vector specific intensity in the trunk medium. The hybrid Monte Carlo-radiative transfer formulation presented here obtains the scattering behavior of the trunk layer by simulation. In this case, the scattered field from a particular layer of cylindrical trunks is given by eikor oo E,(r) = E iko(ki-k).Pn Sn(ks, ki) (6.25) r - n=l where ki and k, are the unit propagation vectors for the incident and scattered waves respectively and pn is the position vector in the xy-plane of the nth cylinder. S, is the effective far-zone scattering amplitude for the nth cylinder, as discussed in the

143 previous section, and represents the current distribution on the cylinder as a result of its coupling to all the other cylinders in the medium. Then the specific intensity incident on the trunk layer may be related to the specific intensity scattered from the trunk layer through the relation It (O,, s) = 2-ct(Os, I; 9i, Xi) It(0i, Xi) (6.26) in which It and It are the modified Stokes vectors for the waves incident on, and scattered from the trunk layer. Ct is the Stokes matrix for the trunk layer and is defined by (I S 12) ( ISvh 2) e (Svv Sh,) - m ( S,, Svh) (IShv 12) (IShh 12) Re(ShvShh) - Sm(ShvShh) 2ReJe(SvvShv) 2~e(SvhSh) W Re(SvvShh + SvhSv) -Sm(SVVSh - SvhSv) 2Qm(SvvShv) 2Tm(SvhSWh)!m(SvvShh + SvhShv)'Re(SvvShh - SvhShv) (6.27) where Spq is the pq scattering matrix element of a particular distribution of cylinders in the trunk medium and < > represents the process of ensemble averaging accomplished by the Monte Carlo simulation. Then the upward-going intensity at the boundary between the trunk and crown layers is considered here to consist of two parts. The first part is the coherent intensity that has been reflected from the specular bottom surface and is subsequently attenuated by the trunk medium. The second part is the component that has been attenuated by the crown medium and is then scattered by the trunk layer. Therefore, the upwardgoing vector specific intensity at the bottom of the crown layer may be expressed in

144 the following form I+(,u, ), -d) = Y(1) It (, 1, -d') + Ct(u, 4; og, so) e-Itd/~o Io (6.28) where It(i, (, -d') is the upward-going vector specific intensity at the ground surface,,o = cos 9j, I, is the incident vector specific intensity and r is the intensity transmission matrix for the coherent wave in the trunk layer as defined in (6.19). The boundary condition at the canopy floor is It (, (, -d') =- R() It(-i, X, -d') (6.29) where I is the downward-going vector specific intensity in the trunk medium and' is the specular reflectivity matrix for the ground as given in [84]. Again, It(-A, q, -d') is considered to consist of two components. The first component is the coherent wave which is first attenuated in the crown medium and then attenuates further as it travels through the trunk medium to the canopy floor. The second component is the incoherent wave that is generated by scattering in the crown medium and attenuates in traversing the trunk medium. These terms can be written formally as I-(-A,1,.-d') - Y(-p) e-cdl/oo 6( - go) 6( - o) +r(-a) d e- (d+,)/ F-(-p,, z') dz' (6.30) Substituting this expression into (6.29) and substituting (6.28) into the first of equations (6.23), we obtain I+(, z) = e-,c(z+d)/lp IZ(p) e-Kcd/Jo 1 S(A - yo)6(q - o) +e-%(z+d)/ e'(p) e-Kc(d+)1 F-(-p, q, z') dz' (6.31) + 12-&c(~-z')/, F+ (,,, z') dz' +e-c(z+d)/ c~t(/L,; /o, qo) e-cd/o Io

145 where IV'(p) = Y(p) 7Z(y) YT(o). To obtain the zeroth-order approximations for the source functions F+ and Fthe expressions for the upward and downward travelling coherent waves in the crown medium, Ic(o, X, z) = eIKc(z+d)/M 7l(p) e- Ccd/M o 6( - _ )5(~ -,) Ico(-, O, z) = elcZ/' Io S(P -,o)6( - 0) (6.32) are used in equations (6.22). The resulting expressions for the source functions are CF (o.., z) = -[Pc(.,; tlo, 0o) e-(c(z+d)/" 7Z'(o,) e-Cd//Io +pc(,,; -,o, 4o) eKCcz/o] Io (6.33) F(j(-., 0, z) = [jpc(-p, 0; o, o,) e-,Kc(z+d)/lo'Z'(1o) e-Kcd/.o +'Pc(-P.,; -Po,,o) e Cz/"o] Io Using these source functions in equations (6.23), the scattered intensity above the crown layer for the hybrid model is found to be to be I.(t., ) = e-Ccd/d Xl'( ) e-cdl/ Io 6(u - po)S(q - _0) + 1 e;t,/ vR'(z~tI[e-lcc(Z2+d)l pc(-At,, a; dtom ZO) e- pc(zz+d)/jojZi(1lo eKCCd/o P -d + e- K(z'+d)/Ipc(-(, P; -p, ko) eCz'/1"~] dz'} Io (6.34) + { [f e c'/"~Pc(#. r; Oo, rOo) e-%c(z'+d)/I~'(Zoi() e-'cd/t" L -d + e Kcz'//pc(P, 0; -,ot, o) eKcCZ'/4o] dz'} Io + e- tcd/lpt(,I, k; go, Oo) e-;cd/'~o Io The integrals in expression (6.34) above can be easily integrated. The results of the integrations are found in [96] and will not be duplicated here.

146 6.4 Results and Discussion In order to evaluate the actual behavior of our expressions (6.12) for the field transmission matrix of the cylinder layer and the related quantity (6.19), the intensity transmission matrix, it is first necessary to sharpen our definition of the effective scattering matrix (S). It will now be defined exactly what is meant by the effective scattering matrix of an average cylinder in a random medium consisting of vertical cylinders. Then the typical behavior for (S) will be shown as the number density of scatterers in the medium increases. The scattered field from a typical cylinder in the medium may be expressed as Es =-kI2 x rx (6.35) where 1le is the electric Hertz vector given by He __ 4 Z J e(r') e-ik~r'cos dS' (6.36) In the expression above, Je is the electric surface current density, Zo is the intrinsic impedance of free-space and ~ is the angle between the vector to the source point and the vector to the field-point. The integration is over the cylinder surface. If we have a z-directed current density that is a function of z' alone on the surface of a thin vertical cylinder located at the origin, the above expression for the Hertz vector becomes He = 2k r I Jz(z') e-ikoz cosG dz' (6.37) where J, is the z-directed surface current density, O, is the elevation angle for the scattered wave, a is the cylinder radius and h is the cylinder height. If this integral is approximated by a sum over n equal segments of the cylinder length, making sure

147 that n is sufficiently large to properly sample the current, we obtain for the scattered field E ik0Zoab eikor ikoznCoSAa (6.38) 2 r in which 6 is the sample segment length, Jn is the z-directed current density sampled at the center of the nth segment and Zn is the coordinate of the nth segment. After taking the dot product of the scattered field with the polarization vector of an incident vertically polarized field we can then extract the vv-polarized scattering matrix element for a particular cylinder. Taking an ensemble average of the current density over all cylinders in the medium, the effective vv-polarized scattering matrix element is obtained for an average cylinder, (S) - 2 sin (J(n)) e-incos (6.39) The scattering matrix element (6.39) is, ideally speaking, what would be used in the first of equations (6.12) to compute the vv-polarized field transmission matrix of the cylinder medium. It is clear from the form of (6.12) that if the field transmission matrix is to remain finite as the number density of scatterers in the medium becomes large, the effective scattering matrix must decrease in amplitude as the effect of multiple scattering increases. In addition, as the elevation angle of the incident wave approaches 900, the cosine in the denominator of equations (6.12) approaches zero. This means that at angles of incidence approaching 90~, the effect of coupling between cylinders becomes quite large and the result should be to reduce the amplitude of (S) to preserve the finite amplitude of (6.12). This hypothesis was tested using the standard NEC method of moments code [6]. Thin, finitely conducting vertical cylinders were distributed within a circular area using a random number generator, and a test cylinder was placed with its base at the

148 0.0200 ~............................ 0.0150 0.0100 A 0.0050 Magnitude of Real Part......... Magnitude of Imag. Part 0.0000. I. I... I. I 0.0 100.0 200.0 300.0 400.0 500.0 Number Density Cylinders (cyls/m2) Figure 6.6: The effective vv-polarized scattering matrix of a distribution of cylinders as a function of number density. The cylinders are 18cm long and 0.1cm in diameter and have a conductivity of 100 mhos per meter. The angle of incidence is 80~ from normal and the frequency is 9.25 GHz. origin of the coordinate system. The number of cylinders to be placed within the area was determined by the desired number density. The cylinders were 5.5 wavelengths long and 0.03 wavelengths in diameter and had a conductivity of 100 mhos per meter. Weakly scattering finite conductivity cylinders were used to eliminate interference caused by the travelling waves that are supported by perfectly conducting cylinders. After a specified distribution was set up by the random number generator, the method of moments was used to determine the current on each segment of the test cylinder. This process was repeated a large number of times to obtain an ensemble average value for the current density on the test cylinder. Equation (6.39) was then used to obtain the effective scattering matrix for this test cylinder. Figure 6.6 shows the result of this procedure for number densities in the range from 0 to 500 cylinders per square meter at an incidence angle of 800 from the normal. At this angle of incidence

149 0.0200 0.0150 Z - ---- Magnitude of Real Part --—......... Magnitude of Imag. Part E 0.0100 A 0.0050 0.0000.. 0.0 100.0 200.0 300.0 400.0 500.0 Number Density Cylinders (cyls/m2) Figure 6.7: The effective vv-polarized scattering matrix of a distribution of cylinders as a function of number density. The cylinders are 18cm long and O.lcm in diameter and have a conductivity of 100 mhos per meter. The angle of incidence is 200 from normal and the frequency is 9.25 GHz. the relatively high level of multiple scattering produces a notable decrease in both the real and imaginary parts of (S) with increasing number density. Figure 6.7 is for the same distribution but at an angle of incidence of 200. At the 20 angle of incidence, the small degree of mutual coupling in the distribution produces little change in the effective scattering matrix with an increase in the number density of scatterers. It should also be noted that the 200 angle of incidence data has a much lower overall magnitude which is consistent with the lower level of interaction. It is expected that stronger scatterers would show a steeper decrease in (S) with packing density and would eventually produce some kind of limiting behavior in equation (6.12), but it is not possible to check this using the NEC code. It is interesting to note that multiple scattering within the cylinder medium produces a decorrelation of the local fields in the medium and this is the field that

150 actually excites the individual cylinders. This field decorrelation property is illustrated in Figure 6.8 which shows the magnitude of the field correlation function for the z-component of the electric field as a function of vertical distance into the medium from the top of the cylinder layer. The field correlation function was generated by creating a random distribution thin conducting cylinders on top of a metallic ground plane. The method of moments code was used to sample the near electric field as a function of vertical distance into the medium at many separate locations. The result was then ensemble averaged and the correlation function determined. It should be realized that the finite correlation length for the electric field in the medium will affect both the scattering and extinction behavior of the medium. It is not actually correct to describe scattering by the individual particles in the medium in terms of their single scattering properties when the medium becomes dense enough to produce significant multiple scattering. In this case the length of the scatterer, as far as the excitation is concerned, is controlled by the coherence length of the fields in the medium. Clearly, for the purpose of simulation of the transmission matrix, it is not possible to obtain all orders of scattering in the medium since an analytical approximation is used to compute the coupling between cylinders. The simulation utilized in this chapter contains the interaction terms for scattering between cylinders up to secondorder. It has been found in previous work [69] that the second-order approximation is valid up to fairly high densities, even for strongly scattering media. In most cases of interest for vegetation, the first-order approximation should provide a reasonable estimate of the co-polarized backscattering coefficients. However, to obtain the proper cross-polarized backscattering coefficients and phase statistics [66], it is necessary to take the second-order terms into account. Figures 6.9 and 6.10 show the vv-polarized forward scattering coefficient for two

151 1.000. 0.750 0.500 U 0.250 0.000 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Vertical Displacement (wavelengths) Figure 6.8: The magnitude of the correlation function for the z-component of the near electric field within a medium consisting of 177 thin conducting cylinders per square meter as computed by the method of moments. The cylinders are 5.5 wavelengths long and 0.03 wavelength in diameter at 9.25 GHz. The excitation is vertically polarized and incident at 300 from the vertical axis. different densities of thin finitely conducting cylinders. The cylinders have been positioned with a random number generator. In the lower density case of Figure 6.9 it can be seen that the single scattering approximation yields a satisfactory estimate of the coefficient while at the higher density level shown in Figure 6.10 the secondorder terms are required to obtain accurate results. The level at which orders of scattering greater than one become necessary is determined by the dielectric constant of the cylindrical scatterers, their heights and the density of packing of the medium. It is a rule of thumb for cylinders with dielectric constants in the region found in tree trunks and vegetation stalks, that when the cone of half-angle equal to the angle of incidence and height equal to the cylinder height encompasses more than about twenty other cylinders it will be necessary to include second-order effects to obtain

152 5.0 0.0 -5.0 ____ SeSond-Order Theory -.. — Frst-Order Theory A MoM Computation -10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.9: The vv-polarized forward scattering coefficient of a layer of 177 thin cylinders per square meter. The cylinders are 5.5 wavelengths long, 0.03 wavelength in diameter and have a conductivity of 100 mhos per meter at 9.25 GHz. The cylinders were positioned with a random number generator. The single scattering approximation provides a good estimate at this fairly high density of scatterers. satisfactory estimates of the co-polarized scattering coefficients. Any satisfactory estimate of phase statistics or cross polarized scattering coefficients will require the use of higher order couplings for all but the most sparse canopies. Figures 6.11 and 6.12 show the vv-polarized intensity transmissivity as a function of layer height for a layer of dielectric cylinders at two different number densities. The cylinders are 0.17 wavelengths in diameter and have a dielectric constant, Cr -35+ill. The frequency of operation is 9.25 GHz and the angle of elevation for the incident wave is 60~ from vertical. In Figure 6.11 the number density is 30 cylinders per square meter which is within the single scattering regime. It is seen that the first-order transmissivity computed using equation (6.19) and a first-order Monte Carlo simulation agrees quite well with the transmissivity computed using

153 5.0 0.0 Socond-Order heory,.-. Fst-Order Theory A MoM Computation -10.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.10: The vv-polarized forward scattering coefficient of a layer of 264 thin cylinders per square meter. The cylinders are 5.5 wavelengths long, 0.03 wavelength in diameter and have a conductivity of 100 mhos per meter at 9.25 GHz. The cylinders were positioned with a random number generator. The second-order scattering approximation is needed here to provide an acceptable estimate of the scattering coefficient. the exponential extinction model of radiative transfer. In Figure 6.12 the secondorder transmissivity computed using equation (6.19) and a second-order Monte Carlo simulation diverges significantly from the result obtained with the radiative transfer extinction model. It is significant that the two models diverge only at relatively high number densities of scatterers. In the first-order scattering regime which one might reasonably expect to exist in most tree canopies, the result of using equation (6.19) and the radiative transfer extinction model should be practically indistinguishable. In dense, stalk dominated, vegetation such as agricultural canopies, one would expect to find that the exponential extinction model overestimates the attenuation of the vertically polarized wave, especially at angles of incidence that are far away from the vertical direction.

154 1.0 0.8 > 0.6 0.4 0.2 ---- Simulation......... Radiative Transfer 0.0 I I.,I I I A I 0.10 0.20 0.30 0.40 0.50 0.60 0.70 Cylinder Height (m) Figure 6.11: The vv-polarized intensity transmission matrix element as a function of cylinder height for a layer consisting of 30 cylinders per square meter. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants E, = 35+ill at 9.25 GHz. The angle of elevation is 60~ from vertical for the incident wave. Figure 6.13 shows the hh-polarized transmissivity for the 90 cylinders per square meter case. At all number densities examined in this study, the hh-polarized transmissivities computed by the method of equation (6.19) and the exponential radiative transfer model are in exact agreement. This is because the approximations used in formulating the extinction matrix of radiative transfer theory are valid for the weak scattering regime presented by thin vertical cylinders to a horizontally polarized wave. Finally, Figure 6.14 shows the general trend of the vv-polarized transmissivity computed using a first-order simulation, for cylinders that are 5.5 wavelengths long. The extinction is seen to increase with increasing number density and the slope becomes steeper with increasing angle of incidence. Both trends are as would be expected from equations (6.12). The most dramatic difference between the Monte Carlo simulation and radiative

155 1.0 0.8 > 0.6 0.4 0.2 - Simulation Radiative Transfer 0.0 I I 0.10 0.20 0.30 0.40 0.50 Cylinder Height (m) Figure 6.12: The vv-polarized intensity transmission matrix element as a function of cylinder height for a layer consisting of 90 cylinders per square meter. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants cr = 35+ill at 9.25 GHz. The angle of elevation is 600 from vertical for the incident wave. transfer results for layers of cylinders is to be found in the treatment of the phase matrix. In radiative transfer theory, the source functions are integrated with respect to the cylindrical axis coordinate to account for volume scattering within the medium. This integration procedure is illustrated in [96]. In Chapters 3 and 5, it was demonstrated that radiative transfer theory produces results that are in sharp disagreement with measured data for layers of vertical cylinders. On the other hand, Monte Carlo simulations have been used to reproduce the scattering behavior of such layers with a fairly high degree of fidelity. To a certain degree, this difference is due to cross coherence between scattering mechanisms such as ground-cylinder and cylinder-ground, and this would account for a baseline offset of about 3 dB if everything else were the same. This is a result of the incoherent nature of radiative transfer theory which does not take into account constructive interference from the various multiple bounce

156 1.0 0.8 > 0.6. 0.4 0.2 - - Simulation -. —----- Radiative Transfer 0.0 I I I I 0.10 0.20 0.30 0.40 0.50 Cylinder Height (m) Figure 6.13: The hh-polarized intensity transmission matrix element as a function of cylinder height for a layer consisting of 90 cylinders per square meter. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants er = 35+ill at 9.25 GHz. The angle of elevation is 600 from vertical for the incident wave. paths and which will result in the radiative transfer results being 3 dB below coherent models. However, the error produced by integration of the source functions for layers of extended scatterers are even more serious. Figures 6.15 and 6.16 illustrate the difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of layers of vertical dielectric cylinders. In Figure 6.15 the dielectric constant is pure real and varies between 3 and 100. The data presented is for an angle of incidence of 20~ from the vertical. At this angle extinction plays a negligible role in the scattering behavior of the layer. The tremendous variability of the difference spectrum is due to the integration of the source functions that occurs in the radiative transfer model. Figure 6.16 shows the results for the same cylinder layer but, in this case, the real part of the relative dielectric constant is fixed at 20 while the imaginary part varies

157 1.0 0.8 > 0.6 - 0.4 30 Degrees Incidence 0.2 -........ 50 Degrees Incidence ------ 70 Degrees Incidence 0.0.. I.... 15.0 30.0 45.0 60.0 75.0 90.0 105.0 Number Density (cyls per sq. m) Figure 6.14: The vv-polarized intensity transmission matrix for a layer of cylinders 5.5 wavelengths in height and 0.17 wavelengths in diameter. The cylinders have a relative dielectric constant of er = 35+ill at 9.25 GHz and the transmissivity of the layer has been computed using a first-order Monte Carlo simulation. between 0 and 20. Figure 6.17 shows the difference between the simulated result and radiative transfer theory for a layer of cylinders with dielectric constants c, = 35+ill as a function of cylinder height. Convergence is always an issue with Monte Carlo simulations. In Chapter 5 it was shown that for simulations of scattering from cylinder layers such as are discussed here, 200 samples is more than sufficient to achieve convergence to within several percent. However, the issue of area convergence has not been discussed up to this point. That is, how big an area is required to achieve convergence of the second-order simulation terms for typical distributions of cylinders investigated in this study? If the area distribution of cylinders is not large enough in extent, the test cylinder used in determining the effective scattering matrix for an average cylinder in the layer will not be truely representative of the ensemble. In other words, the test cylinder will

158 10.0 8.0 VV-Polarized Difference 6.0 4.0 2.0 0.0 0.0 20.0 40.0 60.0 80.0 100.0 Real Part Dielectric Constant Figure 6.15: The difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of a layer of cylinders above a conducting ground plane as a function of the real part of the relative dielectric constant of the cylinders. The imaginary part of the dielectric constant is zero in this case. The cylinders are 5.5 wavelengths long and 0.17 wavelengths in diameter. not be representative of a typical cylinder in a medium of infinite transverse extent. Figure 6.18 shows the result of typical tests for area convergence of the vv-polarized transmissivity. The method of moments has been used to compute the convergence data using equation (6.12) as a basis. The cylinder density is 130 cylinders per square meter, and the relative dielectric constant is er = 1+i1943. This combination provides a level of multiple scattering that is within the second-order regime for angles of incidence below about 700. At angles of incidence below 750 or so, the convergence is to within about a percent for distribution radii below 20cm. However, at 80~, the level of multiple scattering is so high that convergence cannot be achieved to better than 5 percent for distribution radii below 35cm. Figures 6.19 and 6.20 show the comparison between the Monte Carlo simulation,

159 10.0. VV-Polarized Difference 8.0......... HH-Polarized Difference 6.0 4.0 2.0 0.0.... 0.0 5.0 10.0 15.0 20.0 Imaginary Part Dielectric Constant Figure 6.16: The difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of a layer of cylinders above a conducting ground plane as a function of the imaginary part of the relative dielectric constant of the cylinders. The real part of the dielectric constant is 20 in this case. The cylinders are 5.5 wavelengths long and 0.17 wavelengths in diameter. radiative transfer and the hybrid model for a layer of perfectly conducting spheres above a perfectly conducting ground at the moderately low density of 3540 spheres per cubic meter. The spheres have a diameter of 0.2 wavelength and are distributed in a layer 3.1 wavelengths thick at 9.25 GHz. The agreement between all three methods is excellent at this level of number density. Figures 6.21 and 6.22 show the same comparison for a sphere density of 14,147 spheres per cubic meter. Again, the agreement between methods is quite good even at this relatively high distribution density. It is expected that the agreement between radiative transfer and the simulation will break down for extremely high sphere densities due to amplification of the coherence effects. It should also be realized that at high number densities of scatterers, the single scattering representation used here for the sphere layer will no longer be valid. For

160 10.0 VV-Polarized Difference 8.0......... HI-IH-Polarized Diffaerence 6.0 0_ 4.0 2.0 0.0 L * I..I., 0.1 0.2 0.3 0.4 Cylinder Height (m) Figure 6.17: The difference between Monte Carlo simulation results and radiative transfer results for the backscattering coefficient of a layer of cylinders above a conducting ground plane as a function of the layer height. The relative dielectric constant of the cylinders is c, = 35+ill in this case and the cylinder diameter is 0.17 wavelengths. most types of vegetation canopies that have crown layers consisting of small weakly scattering particles in moderate densities, radiative transfer should work quite well for modeling the co-polarized backscattering coefficients. Figures 6.23 and 6.24 show the comparison between methods for a layer of dielectric cylinders on a conducting ground plane. The cylinder number density is 35 cylinders per square meter and is within the single scattering regime of the models. The cylinders are 5.5 wavelengths long and 0.17 wavelengths in diameter and have relative dielectric constants of t, = 35+ill. The hybrid model gives exactly the same result as the simulation but the radiative transfer model differs significantly. At this density level the hybrid model will reproduce the correct scattering and extinction behavior of the canopy with a first-order simulation alone. This is a very fast

161 1.0 0.8 > 0.6''''I —-........................... 0.4 0.2 - 70 Degrees......... 80 Degrees 0.0 I.. I. I. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Distribution Radius (m) Figure 6.18: Area convergence behavior of the vv-polarized transmissivity for a layer of 130 thin cylinders per square meter. The cylinders are 5.5 wavelengths long, 0.03 wavelengths in diameter and have a relative dielectric constant Er = 1+i1943 at 9.25 GHz. The distribution is circular. computation and actually requires less time than the radiative transfer model. Figures 6.25 and 6.26 show a comparison of the models for a low density of cylinders and moderately high density of spheres, while Figures 6.27 and 6.28 give the comparison when the density of cylinders has been increased to a level that is within the second-order scattering regime. In these and subsequent figures, the cylinders and spheres used as canopy components are the same as in Figures 6.19 through 6.24 above. The hybrid model presented in these figures is based on a first-order simulation since it is desirable to keep the computation time down. Because no second-order effects are included in the model, no crosspol information is obtainable. Figure 6.29 is for the same canopy configuration as Figure 6.27 but, in this case, all the secondorder interactions are included in the simulation. The hybrid model contains all the second-order interactions between cylinders but excludes the interaction between

162 5.0 o Monte Carlo Simulation ------—. Hybrid Method - ----- Radiative Transfer -5.0 0 -15.0- o -25.0 -25.0. I. I. I. I. I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.19: Comparison of models for the vv-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 3540 spheres per cubic meter. cylinders and spheres that have been included in the pure simulation. There is not much difference between the hybrid model and the simulation even for fairly high number densities of strongly scattering particles in the crown layer which suggests that coupling between the crown and trunk layers is not an important factor for the co-polarized backscattering coefficient. It therefore seems justifiable to neglect this coupling in the hybrid model. Comparison of Figure 6.29 and Figure 6.27 also reveals that the first-order hybrid scattering model should give a reasonable estimate of the co-polarized backscattering coefficient even when the cylinder densities are within the second-order scattering regime. This is especially true for angles of incidence below around 60~ from the vertical. Finally, Figure 6.30 presents a comparison between between the fully second-order simulation and the second-order hybrid model for the cross-polarized backscatter co

163 5.0 o Monte Carlo Simulation ------—. Hybrid Method ------ Radiative Transfer -5.0 cO 0 -15.0o o -25.0 I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.20: Comparison of models for the hh-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 3540 spheres per cubic meter. 10.0 0 Monte Carlo Simulation......... Hybrid Method -....- Radiative Transfer 0.0 0O -20.0. I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.21: Comparison of models for the vv-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 14,147 spheres per cubic meter.

164 10.0 0 Monte Carlo Simulation. —---—. Hybrid Method ------ Radiative Transfer 0.0 0 0 -10.0 _ o 0 -20.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.22: Comparison of models for the hh-polarized backscattering coefficient of a 3.1 wavelength thick layer of small metal spheres above a perfectly conducting ground plane. The spheres have a diameter of 0.2 wavelength and the number density is 14,147 spheres per cubic meter. 10.0 o Monte Carlo Simulation......... Hybrid Method -....- Radiative Transfer 0.0 G...... -10.0,, -20.0 -20.0. I I I.I I. I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.23: Comparison of models for the vv-polarized backscattering coefficient of a layer of vertical cylinders 5.5 wavelengths high above a conducting ground plane. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants Er = 35+ill at 9.25 GHz. The number density is 35 cylinders per square meter.

165 10.0 0 Monte Carlo Simulation......... Hybrid Method.....- Radiative Transfer 0.0 -10.0 -20.0. 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.24: Comparison of models for the hh-polarized backscattering coefficient of a layer of vertical cylinders 5.5 wavelengths high above a conducting ground plane. The cylinders are 0.17 wavelengths in diameter and have relative dielectric constants,: = 35+ill at 9.25 GHz. The number density is 35 cylinders per square meter. efficient. The effect of coupling between the crown and trunk layers is seen to be significant as far as depolarization is concerned, but the second-order hybrid model still provides a rough estimate of the cross-polarized backscatter coefficient while the radiative transfer model provides no depolarization information at all. It is to be expected that the second-order hybrid model will improve relative to the simulation as far as depolarization is concerned as the number densities and single scattering albedos of particles in the crown layer decrease. In any event, it is not very advantageous to employ a second-order hybrid model because the CPU time it requires tends to be rather large. It takes about 12 hours of run time on an IBM RS6000 workstation to compute the second-order hybrid model of Figures 6.29 and 6.30. If crosspol information is not essential, the first-order hybrid model will provide fairly good accuracy for

166 10.0 o Monte Carlo Simulation ------—. Hybrid Method ------ Radiative Transfer 0.0 m e —— e ------------------— ~ — -10.0 - -20.0 I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.25: Comparison of models for the vv-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 35 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter. most conditions found in typical vegetation canopies and the computation time will be of the same order of magnitude as that required for the radiative transfer model. 6.5 Conclusions This chapter has presented a hybrid model for the computation of scattering from layered vegetation media consisting of vertical trunks above a dielectric ground plane and a crown layer consisting of small weakly scattering particles. The crown layer has been modeled using radiative transfer theory and the trunk layer has been simulated with a Monte Carlo simulation. The derivation of the transmissivity matrix for the trunk layer has been presented and some of its important features investigated. The concept of the effective scattering matrix for an average cylinder in the trunk medium has been developed and evidence has been presented that it assumes limiting

167 10.0 o Mcnte Carlo b Auation Hybrid Method ------ Radiative Transfer 0.0 o o o -....... -10.0 -20.0 I. I. I, 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.26: Comparison of models for the hh-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 35 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter. 10.0. 0~ o - _9 —-''O Lj —-10.0 0 Monte Carlo Snimultion -------- Hybrid Method - ----- Radiative Transfer -20.0. I. I. I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.27: Comparison of models for the vv-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 106 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter.

168 10.0 0.0 - -10.0 o Monte Carlo Simulation -......... Hybrid Method ------ Radiative Transfer -20.0 I I. I. 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.28: Comparison of models for the hh-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. The density of the lower layer is 106 cylinders per square meter and the density of the upper layer is 14,147 spheres per cubic meter. 10.0. -, z.0 -10.0 _ 0 Monte Carlo Simulation ------—. Hybrid Method ----- Radiative Transfer -20.0. I - I I I 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.29: Comparison of models for the vv-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. All the secondorder interactions between spheres and cylinders and between cylinders alone have been included in the simulation while the hybrid model includes only the second-order coupling between cylinders.

169 -20.0 O O G 0 -30.0 - o -40.0 - "-, ~0 -50.0 " o Monte Calo Simulation ---—.........- Hybrid Mcthod -60.0 _, -70.0 20.0 30.0 40.0 50.0 60.0 70.0 80.0 Elevation Angle (deg) Figure 6.30: Comparison of the hybrid model and a Monte Carlo simulation for the vh-polarized backscattering coefficient of a two-layer canopy consisting of vertical cylinders on a conducting ground plane and a crown layer composed of small metal spheres. All the second-order interactions between spheres and cylinders and between cylinders alone have been included in the simulation while the hybrid model includes only the second-order coupling between cylinders. behavior as the effect of multiple scattering in the medium becomes more important. Evidence has also been presented to confirm that the exponential extinction model used in radiative transfer theory works fairly well for sparse distributions of particles and even for sparse distributions of extended scatterers such as cylinders, however the radiative transfer extinction model breaks down for high densities of strongly scattering particles that are large compared with the wavelength of the excitation. In addition, it has been demonstrated that the source function integration used in radiative transfer theory to account for volume scattering in random media leads to results that are severely in error for layers of long cylinders, and it may be inferred that this is also true for other types of large scatterers. Finally, several varieties

170 of test canopies consisting of dielectric cylinders and small metallic spheres over a conducting ground plane have been investigated. It has been shown that the hybrid model developed in this chapter can be an effective way to compute scattering from vegetation media with this general structure.

CHAPTER VII AN ITERATIVE INVERSION ALGORITHM WITH APPLICATION TO THE RADAR RESPONSE OF VEGETATION CANOPIES 7.1 Introduction In recent years a great deal of emphasis has been placed on the retrieval of information from synthetic aperture radars and radar polarimeters [93, 94, 101]. Some of the applications have been in remote sensing of soil moisture for bare soil [57, 75] and vegetation-covered soil [102, 28]. Other applications include the determination of vegetation canopy parameters [53, 48, 10, 35], sea surface characteristics [34, 11], and snow parameters [51, 99]. Within the general problem of classification of remotely sensed data, there exists the sub-problem of inversion of radar data to obtain parameters of interest for the scene under observation. The vast majority of the literature in this area has been concerned with two major approaches to the problem of inversion of radar data. The first approach involves the construction of an empirical scattering model specific to the type of problem being studied [57, 1, 76]. In this technique the scattering characteristics of a particular type of terrain or vegetation canopy are determined experimentally, usually at several frequencies using a polarimetric radar scatterometer. 171

172 The results are then fit to a fairly simple equation or a set of equations that describe the scattering behavior as a function of polarization and frequency over a specified range of parameters (region of validity) for the type of terrain or canopy being characterized. The empirical model obtained in this way is designed to be invertible over its region of validity. This type of technique can give accurate results for the case it has been designed to treat but is completely specific to that case and provides little physical insight. The second approach, which has received a large amount of attention lately, is based on the use of artificial neural networks [59, 14, 65, 52]. The transfer function of the network is determined by using training sequences of known input and output data. In the present case this would consist of polarimetric radar data as inputs and scene parameters as outputs. The network characteristics relating the radar data to the scene parameters are determined using a back propagation algorithm to successively refine the transfer function based on the set of training sequences. After the network has been trained it can be used to estimate unknown scene parameters given an input set of radar data. The neural network has great flexibility and, depending on the set of training data, the results may be excellent [59, 33]. However, the neural network is, in most cases, used essentially as a black box. There is no way currently known of discerning the underlying physical processes that give rise to the network behavior, which means that there is no systematic way of selecting the optimum set of data channels for use in a neural network-based inversion given the network response alone. To do this, one must ultimately rely on information provided by theoretical models and/or measured data. There is also no way of determining if the decision path taken by the neural network in arriving at its result is a reasonable one. In this chapter an iterative algorithm for the inversion of polarimetric radar data

173 is presented. The algorithm is completely general and may be applied, in principle, to any type of radar scattering problem for which a model exists. The behavior of this algorithm is derived entirely from the physical scattering mechanisms existing within the system being studied and represents the summation of knowledge gained in applying the scattering model across a representative range of states of the system. The amount of computation time required in applying the algorithm to any given case is small since all the information necessary for inversion of a set of input data has been precomputed. It is also possible to monitor the decision path taken by the algorithm in arriving at a result and therefore some measure of control over the reliability of such results is achievable. The second and third sections of this chapter present the theoretical development of the iterative inversion algorithm and some important considerations involved in its implementation. Sections 4 and 5 discuss the application of the algorithm to the particular case of retrieval of vegetation canopy parameters from radar data and give the results of an application based on inversion of the radiative transfer model for a simplified representative canopy. The final section presents the results of an error analysis for the algorithm in terms of both systematic and measured quantities and discusses how sensor configuration, algorithm implementation and data uncertainty influence the inversion accuracy. 7.2 Iterative Algorithm Let us assume that backscatter data is provided for a particular target at a given frequency and for a range of incidence angles 0 E [a, b]. This data may be represented in terms of a Fourier series in the restricted angular range. The Fourier coefficients obtained in this way provide an equivalent representation of the system response for

174 this type of target. The behavior of the Fourier coefficients as a function of the target parameters provides a convenient measure of the system response that is independent of incidence angle. If a subset of the Fourier coefficients can be found that represents the angular response of the system to sufficient accuracy and their behavior is known for all values of the target parameters, then in effect we have an empirical model that describes the target behavior in detail. In addition, if the functional forms of these Fourier coefficients with respect to the target parameters are known, it should be possible to construct subsets of the parameter space over which piecewise linear representations of the coefficients are obtainable. The resulting sets of linear equations should be directly invertible and, by means of an iterative process that successively converges on smaller domains, a solution set is found. Because it is impossible to determine experimentally the behavior of a complicated physical system under all conditions, a model that represents its behavior as a function of all its important parameters is used in the construction of the empirical Fourier representation. This leaves the determination of the optimum sensor configuration to the discretion of the system designer. Let us consider then a model M that operates on parameters [cl,..., am] and the angle 8 to produce an estimate of a' for a particular scene: a(a1,,..Ima, 8)= M(at... am, a ) (7.1) The range of validity of the model is now restricted to a subset of the range [0, ]. o e [a,b] - 9' E [oj] (7.2) with 0' = 3(s - a) such that 0' = O when 0 = a and 0' = when 8 = b.

175 Thus, we find that / = 2(b-a) and 0 = a + 0'/P. We now construct the function g(O') = ~o(a + 0'//3). Since the domain of g is [0, 2], we can expand g as a Fourier cosine series in 0': 1 00oo g(O') = co + E C, cos(2nO') (7.3) n=l where Cn - g (0')cos(2n') dO' (7.4) We can then write: 1 00 a(0) = CO + Can cos[2n3(0 - a)] (7.5) n=l with C, j 1 0(80) cos[2n#(0 - a)] dO (7.6) These are the coefficients for the restricted range Fourier series description of 6~(8). We describe the general polarimetric model response as a function of angle in the angular range [a, b] by its restricted range Fourier series: 1 N (l,..., m; f, 0) = CP + E Cpq cos[2n/3(0 -a)] (7.7) n=l with cP-Cq (cpq(l,..., am; f) where cal,..., am are the model parameters, f is the frequency and p,q is each either v or h for vertical or horizontal polarization. The series has been truncated at the (N + l)th term which is assumed to give a satisfactory approximation to the system response. This representation generally consists of six or seven terms in the case of

176 vegetation since the angular behavior of the radar cross section of such canopies away from normal incidence is a gentle function of 9. Thus, the unique set of Fourier coefficients has been determined for the model response as a function of angle within the restricted angular range for any given set of input parameters. By way of an example, the behavior of the first four Fourier coefficients as a function of the volumetric soil moisture beneath a vegetation canopy simulated using vector radiative transfer theory at 1.5 GHz for vv and vh polarizations is shown in Figures 7.1 and 7.2. It is seen that the functional behavior of the Fourier coefficients is fairly linear over the entire operational range of volumetric soil moisture and that the convergence of the Fourier series is rapid. This has been found to be true in most vegetation canopy applications. Figures 7.3 and 7.4 show the magnitude angular response of the radiative transfer model as compared with a Fourier series representation utilizing many terms and synthetic data created by using the four coefficient series of Figures 7.1 and 7.2. In this case the four term approximation agrees with the model to better than 0.25 dB over the entire angular range. The behavior of the CP(al,..., am; f)'s is now approximated as linear functions of the cas over the initial range of these parameters. That is, it is known initially, from experience, that each ck falls into a range kin < ak < ak; which is not unreasonable in, for example, typical vegetation canopies of this kind at this time within the growing season. Then we say that each initial range has a centroid or central value of that parameter range. The centroid for that parameter range is denoted by ac and is defined by: = (c4'ak + crn)/2 (7.8) To the first-order of approximation, the Taylor series expansion of the Fourier

177 3.0 Coef. 1 2.0 -- Coef. 2 ------ Coef. 3 V -Coef. 4 1.0 i......... 0.0 - - - - - - - - - - - - - - - -1.0 0.0 0.1 0.2 0.3 0.4 Soil Volumetric Moisture Figure 7.1: The first four vv-polarized Fourier coefficients of a simulated vegetation canopy as a function of volumetric soil moisture at 1.5 GHz. 0.6 0.5 Coef. 1 -— Coef. 2 0.4 ------ Coef. 3 -.... Coef. 4 0.3 1 0.2 0.2, ----........ - 0.1 0.0 -- i 0.0 0.1 0.2 0.3 0.4 Soil Volumetric Moisture Figure 7.2: The first four vh-polarized Fourier coefficients of a simulated vegetation canopy as a function of volumetric soil moisture at 1.5 GHz.

178 2.0 — 0 Fourier Series-, 0.0 -2.0 30.0 40.0 50.0 60.0 -4.0 Model Data -6.0 — 50 F e S Fourier Series ------ - - Synthetic Data -7.0 -10.0 -13.0. I, I, I 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg) Figure 7.4: Comparison of radiative transfer model with Fourier representations for vv polarization at 1.5 GHz. -3.0 Model Data -5.0 Fourier Series Synthetic Data -7.0 -9.0 -11.0 -13.0 -15.0 20.0 30.0 40.0 50.0 60.0 Incidence Angle (deg)

179 coefficients about the range centroids is given by: Cnp(~)- (C p+) E m (Cnpq- ) (7.9) k=lk1Oak) 0 where ca= (ca,..., m), c~ (c1,... ao) and (CP)o =C Pq(a). The coefficients, ( )k _=,o are determined by varying each ak over its range while holding all other a's fixed at their centroid values. At each value of ak in its range, the model is evaluated as a function of 0 and the Fourier coefficients are extracted. In this way all the partial derivatives may be computed from a linear least-squares fit of the coefficients for each parameter to be estimated. From this discussion it is seen that'q has, in effect, been linearized as a function of the model parameters for any arbitrary value of incidence angle in the range [a, b]. Thus we write: apq(0) -_ (q(0)) + E(r q(O))(ak- ao) (7.10) k=1 where (A 1 N 2 Cq)o + (q)o cos[2n3(O - a)] (7.11) n=l and pq (0)))> 1 Cgopq) N (On (,k(o)) 2 (\ ak) +E( k cos[2n3(a - a)] (7.12) It is possible to divide all the initial parameter ranges ak e [tm ina Otkax] into two equal subranges for each parameter. These subranges are denoted by {ak(l, 1)} = [~ak, ai] and {cak(1,2)} = [aot, ara]. Each of these new subranges has associated with it a centroid. These centroids are denoted by a?)(l) and aO1)(2), respectively. It is now possible to construct 2m new centroid vectors consisting of

180 2I''. I a) i.: m-tuples of centroids of the parameter subranges. According to the procedure deIscribed previously, the Fourier coefficients and gradients about the 2 new cen0: a two p arameter space. Now, let us assume we have m measurements of a; a.11 m... 12(2.............. a 2plpe *,,., p,,(mf) h: __ i_ m-tuples o centro of the paramte isubranges According to te roereo de- troidsuare 7.5: () 2- p...,r4m(tp) T.e.rest.showing (2p..... centroid. vector s may

181 by: 0a = (~) + (/7) (& - c~) (7.13) where (8-) - ((af'>,..., (ao))T (7.14) from (7.11), and (R)ij is the value of (rP') (0j,fj)) in the matrix of Fourier expansions of the coefficient slopes as in (7.12). a - (a&,..., &m) is the estimate vector for the parameter set and ca is the centroid vector for the set of parameter ranges. Then the zeroth-order estimate for the parameter set may be found from: &(0) = (Z)l (or - ()r~o) + _o~(O) (7.15) where ac(O) is the centroid vector for the initial set of parameter ranges, (t,)o is the vector (7.14) for these initial ranges, and (rIo is the matrix whose elements are the original values of (7.12) formed as described above. This estimate will have errors due to the fact that the Fourier coefficients computed using the initial parameter ranges are not actually linear functions of the parameters. However, it may be expected that for parameter ranges sufficiently small this first-order Taylor expansion approximation will become increasingly accurate. This forms the basis for an iterative algorithm. The parameter space is discretized to produce successively higher-order centroid vectors about subrange spaces of successively smaller volume as described previously. Within these successively smaller parameter subspaces, the first-order Taylor representation of the Fourier coefficients becomes an increasingly better approximation. Then let us suppose that we have a performance measure P(Ia - a'I) that decreases monotonically with argument and that satsfies P(O) = Pmin. The first-order estimate of the parameter set is then

182 computed as follows: A(1) = (7R)'1 (a,~ - (~)1) + ac(opt) (7.16) where car(opt) is a first-order centroid vector that satisfies: P(&t(0) - aC(i)l) > P(I&(0) - al(opt)j) > Pmj, V i E {1,...,2m} (7.17) and the other quantities with subscript "1" are the same type as in (7.15) with values computed at cxl(opt). This process may be repeated as many times as desired up to the highest order of discretization of the model parameter space, or until the rms difference between successive solution vectors stabilizes to within some arbitrary percentage. 7.3 Implementation of the Algorithm In practice, it may not be necessary to compute the Fourier representation of the model behavior. In many cases, one is constrained to the use of a particular sensor configuration. In such cases there is a limited set of channels available to users of the data. For example, given a polarimetric SAR with one look angle and two frequencies of operation, and ruling out the use of phase information, there would be a maximum of six channels (four copolarized and two cross polarized amplitudes) available for use in inversion. If there are more data channels available than there are parameters of interest for the system under observation, then it would be useful to perform a sensitivity analysis on the model to determine which of the available channels will provide an optimum set. In any event, if the set of data channels available for use is pre-determined then the piecewise linear representation of the model must be based on information provided by those channels and equation (7.13) may be inverted directly.

183 This reduces the number of values that need to be pre-computed by a factor equal to the number of coefficients required in the Fourier representation. Whether the Fourier representation is used or not, the matrix of slopes and vector of intercepts as indicated in (7.13) must be generated from the pre-computed model outputs for each centroid subrange in the parameter space, the number of such subranges being determined by the maximum level of discretization of the space. The intercept vector for a particular centroid is a fixed quantity which depends on the set of model outputs obtained for the inputs determined by the position of that centroid in parameter space. The slope matrix on the other hand is dependent on the set of data points used in sampling the parameter subranges connected with a particular centroid. The functional behavior of the model within a subrange of one of the parameters is considered to be linear, and the partial derivative with respect to that parameter is determined by a weighted linear least squares fit of the model outputs holding all other parameters fixed at the centroid value and varying the parameter in question over its subrange. The point that represents the centroid value is strongly weighted relative to the other points because that is the point about which the derivative must be computed. There are then at least two basic types of optimization that may be utilized in implementing the iterative inversion scheme. The first type will be referred to as intra-centroid optimization. In intra-centroid optimization the data points used in computing the linear fit (derivatives) about a particular centroid are initially equally weighted except that the centroid itself is strongly weighted. In succesive iterations, points that are found to be most distant from the refined estimate are assigned increasingly lower weights until the estimated value of the parameter vector no longer changes by more than a specified amount. The second type of optimization, referred

184 to here as inter-centroid optimization, is utilized after applying the intra-centroid scheme. If the refined estimate obtained using intra-centroid optimization represents a state of the system closer in parameter space to another centroid on either the same or a different level as was shown in Figure 7.5 the algorithm uses this new centroid as the basis for further refinement. It is evident that the algorithm may "wander" from centroid to centroid until it finds a point of local stability on which to converge. The convergence behavior of the algorithm using these optimization schemes will be discussed in a later section. Convergence of the algorithm itself is fairly rapid since most of the CPU intensive computations are pre-evaluated using the scattering model for the system under investigation and are stored on disk as data arrays to be loaded into memory at runtime. In this way systems of functions that are fairly non-linear may be inverted in a direct way, although the more the model behavior deviates from linearity in a given parameter, the less likely it will be to obtain a unique solution for that parameter. 7.4 Application to a Vegetation Canopy Modeling of vegetation canopies using the radiative transfer approach has become increasingly popular in recent years [39, 96, 53, 54, 40]. The vector radiative transfer equations take into account the individual scattering properties of the vegetation canopy components through a phase matrix which relates the incident to scattered intensities as a function of the wave directions, scatterer geometries, and the material composition of the scatterers. This formulation has the advantage of being general, mathematically tractable, and computationally fairly non-intensive which makes it convenient from the standpoint of inversion. In this section the iterative algorithm developed above is used to invert the radia

185 Crown Layer Trunk - Ground Layer Figure 7.6: Layered structure of the radiative transfer model for a vegetation canopy. tive transfer model of a simplified vegetation canopy. The canopy consists of leaves and vertical trunks above a rough ground layer. The layered structure of the canopy model is illustrated in Figure 7.6. The leaves are considered to be uniformly distributed in orientation and are modeled as thin dielectric sheets using physical optics to generate the scattering matrix [74, 71]. The scattered field from the rough ground surface is generated using the Kirchoff and scalar approximations [92]. Trunks are treated as finite length circular dielectric cylinders for which it has been shown that an approximate solution for the scattered field can be obtained based on the solution for the infinite length case provided the diameters of the cylinders are much smaller than their lengths [100, 64, 41]. The dielectric function of all canopy components is calculated using the dual-dispersion model of EI-Rayes and Ulaby [91, 29]. In this simplified vegetation canopy, five parameters are considered as variables: volu metric

186 Variable Parameter Parameter Range Volumetric Soil'Moisture 0.05 - 0.33 Ca.nopy Density: Trunks 0.05 - 0.25 trilnks/rn2 Trunk Height 2./10 - 3.68 m Trunk Diameter 11.0 - 1/1.2 cm Canopy Density: Leaves 100 - 300 leaves/n7, 3 Table 7.1: Variable parameter set for the vegetation canopy used to test the inversion algorithm. soil moisture, trunk canopy density, trunk height, trunk diameter and leaf canopy density. Actually, there are only four independent variables since the trunk height is usually related to the trunk diameter through a function. The ranges of these variable parameters and the values of the parameters considered to be fixed are listed in Tables 7.1 and 7.2. While the ranges of some of these parameters may not be exceptionally representative of all typical real canopies, they do model the important general features of electromagnetic scattering from vegetation layers and serve to demonstrate the invertibility of the radiative transfer equations using the iterative technique developed in this chapter. The range of volumetric soil moistures given is representative of the entire wetness scale for typical soils, from that just following a rain storm to drought conditions [2]. The trunk canopy density range is for conditions from fairly dense woods to sparsely wooded areas [12]. The number density of leaves in the canopy, considered together with the size of the leaves and the height of the crown layer combine to give a range of leaf area index (LAI) from 2.0-5.9 which covers almost the entire range of this parameter [106]. The trunk height and diam

187 Fixed Parameter Parameter Vallle Trunk Gravimeetric Moistlire 0.7 Leaf Gravimnetric Moistulre 0.7 Dry Density of Plant Ma~teria~ls 0.33 gmi/cc T,ea.f Diameter 10.0 cm T,eaf Thickness 0.02,5 cm Crown Height 2.5 m Slurface RMS Rolighness 1./1 cm Sllrface Correlation length 20.0 cm Soil Composition sand: 23.6 % clay: 33.7 % Table 7.2: Fixed parameter set for the vegetation canopy used to test the inversion algorithm. eter ranges are not extremely representative of an average forest canopy but would be more in accord with what one would expect to find in an orchard or grove [53]. Values of 0.7 for the gravimetric moisture of leaves and trunks would give an L-band dielectric constant of about 28-j9 and a C-band dielectric constant of around 25-j9 for both canopy constituents using the model of El-Rayes and Ulaby. This is consistent with typical measured data for trees [53]. Soil surface characteristics are within the validity region for the physical optics model [92]. To test the inversion algorithm it was first necessary to select four radar channels as data sources to be used as model outputs and inversion algorithm inputs. These four channels provide the correct number of equations to invert the linear system (7.13) for the four variable model parameters. The assumption was made that a pre-determined

188 system configuration consisting of an L-band radar operating at 1.5 GHz, and a Cband radar operating at 5.0 GHz were available for use. Further, an analysis taking into account the overall backscatter signal levels and sensitivity to variation of the inversion parameters indicated that an incidence angle of 45 degrees would be useful in this application. It was decided to use data from the co-polarized channels because they produce the highest signal levels and they can, in practice, be calibrated with the greatest accuracy. The four-parameter space was discretized up to the second level producing sixteen secondary range centroids. The number of centroids at each level is given by n = 2m(l-1), where m is the number of parameters and I is the level of discretization. The model was then run in the forward direction five times within each full parameter range for the primary centroid and three times within each subrange for every secondary centroid to provide the pre-computed data used in determining the slope matrices and intercept vectors. After the necessary information was produced to provide a mapping between the parameter space and the model output space, the model was then used as a data simulator for use in testing the inversion. 7.5 Inversion Results The region in parameter space for which the least inversion accuracy is to be expected is the region near the primary range centroid. In this region the estimate of the parameter vector is based on only a single level of refinement and intra-centroid optimization. Because the primary centroid is the nearest centroid to the estimated parameter vector in this case, the algorithm does not take advantage of the higher degree of discretization available on the second level, and inter-centroid optimization is not used. When the estimated parameter vector is closer in parameter space to

189 any higher level centroid the inter-centroid optimization may then be used to obtain further refinement of the estimate. Typical results for inversions obtained near the primary range centroid are shown in Figures 7.7 through 7.10. It is seen that primary level optimization alone provides excellent inversion accuracy for most of the parameters over almost the entire parameter range. However, there is a significant degree of instability in the trunk height determination for the region below about 2.7 meters (Figure 7.10). This is due to the small range of heights considered in this study and the fact that the radar response is not a sensitive function of either trunk height or diameter within the selected ranges of these two parameters for the data channels used in the inversion. This problem could be remedied by using a finer division of the parameter space on the second level so as to provide a convenient secondary range centroid corresponding to the position of the primary range centroid but with a restricted parameter sub-range. A larger range of trunk heights and diameters such as one might find in an actual forest, or a set of data channels more sensitive to these parameters, would also improve the relative accuracy of the trunk height estimate. The excellent results obtained for regions close in parameter space to the primary centroid are actually representative of the worst case estimates produced by the algorithm. Typically, when the algorithm is able to take advantage of a higher level of optimization by using a convenient second level centroid, excellent convergence stability is achieved. This is illustrated in Figures 7.11 and 7.12 which demonstrate that the stability problems encountered in using the entire parameter range of the primary centroid are greatly reduced when the parameter sub-ranges of the second level are utilized. The iterative inversion algorithm exhibits the four general types of convergence behavior shown in Figures 7.13 through 7.16. In about 90 to 95 percent of the cases studied, the algorithm converges

190 0.35 0.30 0.25.0 0.20 ~ - 0.15 Th 0.10 / Model Response El Algorithm Response O 0.05 _ 0.00.. I. I. I 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Volumetric Soil Moisture (Actual) Figure 7.7: The determination of soil moisture over the entire parameter range near the primary centroid. 0.275.. 1, 0.225 *.,, 0.175'~ 0.125 C 0.075 -- Model Response E 0.075 ~ /.a Algorithm Response 0.025 0.025 0.075 0.125 0.175 0.225 0.275 Actual Trunk Density (trunks/m2) Figure 7.8: The determination of trunk canopy density over the entire parameter range near the primary centroid.

191 350.0 E 300.0 250.0 a 200.0 150.0 a -150.0 t Model Response IE a Algorithm Response 100.0 - 50.0 50.0 100.0 150.0 200.0 250.0 300.0 350.0 Actual Leaf Density (leaves/m3) Figure 7.9: The determination of leaf canopy density over the entire parameter range near the primary centroid. 3.840... 3.520 3.200 o 2.880 E uS D — Model Response t; 2.560 a Algorithm Response 2.240 2.240 2.560 2.880 3.200 3.520 3.840 Actual Trunk Height (m) Figure 7.10: The determination of trunk height over the entire parameter range near the primary centroid.

192 0.35 I I 0.30. 0.025 0.20'~ 0.15 0.10 Model Response o Algorithm Response 0.05 0.00 0.00. I,, I. I ~ I 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 Actual Soil Moisture Figure 7.11: The determination of leaf canopy density over a second level parameter sub-range. 3.84 ____ Model Response a Algorithm Response 2.88. I 2.88 3.20 3.52 3.84 Actual Trunk Height (m) Figure 7.12: The determination of trunk height over a second level parameter subrange.

193 0.75 0.60'6E 0.45 8 " CT b —---— 0 —--— a 5' 0.30 I -;>. — Actual Value 0.15 — o-. —. Estimated Value 0.00. I. I. I 0 2 4 6 8 10 Iteration Number Figure 7.13: Uniform convergence of the iterative inversion algorithm. to the correct value as shown in Figures 7.13 and 7.14. In the remaining cases, the algorithm either does not converge at all as shown in Figure 7.16 or converges to the wrong value as in Figure 7.15. Of these last two situations it is definitely preferable to have non-convergence of the algorithm. With non-convergence it is clear that the inversion must be rejected and thus an erroneous result is easily avoided. Anomalous convergence represents less than 5 percent of all convergence behavior exhibited in this implementation of the iterative inversion algorithm. With anomalous convergence it is not always possible to distinguish an incorrect result, therefore this type of convergence behavior represents the most difficult problem to correct in applying the algorithm. Because the convergence behavior of the algorithm is to some degree implementation specific there are also some steps that can be taken to reduce the risk of obtaining incorrect results. It should be possible to utilize other types of algorithms to monitor the decision processes that the inversion algorithm uses in

194 0.75 1.0 - A,, Actual Value a', -— 0.-.- Estimated Value 0.45 E!. *= 0.30.. Q~ t o.', 0.0015 0.00. I. I, I, I. I. I, 0 2 4 6 8 10 12 14 Iteration Number Figure 7.14: Non-uniform convergence of the iterative inversion algorithm. 5.0 I 4.0 I~P ___ ------ 3.0 Actual Value -—. —-. Estimated Value 2.0 0 2 4 6 8 10 Iteration Number Figure 7.15: Anomalous convergence of the iterative inversion algorithm.

195 600.0 500.0 400.0': 300.0 Actual Value 200.0 -o-e —-. Estimated Value 100.0 I I I. I 0 2 4 6 8 10 12 Iteration Number Figure 7.16: Oscillatory non-convergence of the iterative inversion algorithm. achieving convergence. In this way it would be possible to minimize the probability of taking an incorrect decision path. However, the investigation of this approach is outside the scope of this work. 7.6 Error Analysis A complete analysis of errors for the inversion algorithm is lengthy and the derivation will not be given here. The rms error in the ith parameter, ai will be denoted E (ti). This error originates from two major sources. The first source is the overall systematic error which depends on the combined data channel characteristics, the accuracy of the scattering model and the degree of discretization utilized in the algorithm implementation. The systematic error in parameter ai is directly proportional to the offsets of the individual data channels from their respective range centroids for a given level of discretization. This offset is defined as t= (A- ()) for the

196 jth data channel as in equation (7.15). The second source of error in the inversion algorithm derives from the measurement uncertainties in the data channels. The measurement uncertainty in ac for the jth data channel, em(<j) produces an error in the inversion dependent on the gain characteristics of that channel. Both sources of error are inversely proportional to the determinant A of the slope matrix r as given in equation (7.13). The expression for the rms error in the parameter ai for an inversion using m data channels is given by =1 m E (as) = { Z [i?. (/7)2 + cofL2(i') Emi(O.;)] }1/2 (7.18) Where cofji(r) is the ji cofactor of?, and the factor,uij represents the overall uncertainty in parameter ai produced by an error in estimating the response of data channel j to that parameter. It is defined by #ij = [ E2(cofji(7Z)) + cof i (' ) ( (A\ )/A)2 ]1/2 (7.19) where (cofji('7)) is the error in cof, (TZ) and e (A) is the uncertainty in the determinant of the slope matrix. Examination of equations (7.18) and (7.19) reveal some important factors that affect the accuracy of the inversion. The overall rms system error will be minimized when A is maximized, the E (cofji(IZ))'s are minimized and the individual slope matrix element errors that contribute to e (A) are minimized simultaneously. The magnitude of A is a function of the data channel gain characteristics and the structure of the linear system of equations provided by the combination of channels used in the inversion. The slope matrix determinant is maximized when the overall channel gains are as large as possible and the system of equations is as close to being linearly independent as is practical for the set of sensors available. The slope matrix errors

197 are minimized as the linear regression slopes approach the local derivatives in the proximity of the parameter estimate vector. This occurs as the level of optimization increases. The upper limit on the reduction of systematic error as a result of increasing the level of discretization will ultimately be determined by the contribution to the c (cofji(R))'s originating from the model itself. If the model is inaccurate, no amount of optimization will permit accurate inversion of the measured data. However, if the data channels are chosen appropriately and the model is sufficiently accurate, the systematic error can be minimized by providing a suitable degree of discretization of the parameter space. In the case of systematic errors it is to be expected that the best error performance will occur when the actual parameter vector lies close to a range centroid. The second contribution to (7.18) is dependent on the measurement uncertainty in the data channels. If the data channels have been properly selected and the algorithm has been implemented so as to provide sufficient discretization of the parameter space, then measurement errors can be expected to give the largest contribution to the overall inversion error. These errors have a slight dependence on algorithm implementation but are largely influenced by the ratio of the measurement uncertainty to the sensitivity of the data channels. This type of error response is illustrated in Figures 7.17 and 7.18 which show typical error bounds for the inversion if there is a uniform ~0.5 dB measurement uncertainty in all of the data channels. If soil moisture were the only variable parameter in the inversion, measurement accuracy would control the error response to yield a worst case average uncertainty of about ~3.5%. The inversion accuracy for the leaf canopy density, with the density fixed at midrange, is within about ~2.0%. In this case the systematic error and the measurement uncertainty error are of the same order of magnitude.

198 0.40 Lower Bound (0.5dB Accuracy) -......- Ideal Case (No Measurement Error),/ ----- UpperBound(0.5dB Accuracy) v i 0.30 ()o 0.20 7 rii 0.10 0.00 0.00 0.10 0.20 0.30 0.40 Actual Soil Moisture Figure 7.17: Typical inversion error bounds for variable soil moisture with other parameters held fixed. 300.0 sE --- Upper Bound (O.5dB Accuracy) > 250.0 --- Ideal Case (No Measurement Error) ----- Lower Bound (0.5dB Accuracy) 200.0 3I 150.0 100.0. I. 0.00 0.10 0.20 0.30 0.40 Actual Soil Moisture Figure 7.18: Typical inversion error bounds for fixed leaf density with variation of the soil moisture parameter.

199 7.7 Conclusions In this chapter a general, model-based, iterative algorithm has been presented for use in the inversion of polarimetric radar data. The algorithm was implemented using two levels of discretization and two types of optimization and was applied to the case of a general representative vegetation canopy. This simplified test canopy was modeled using vector radiative transfer theory and consisted of vertical trunks, leaves and a rough ground surface. Four canopy parameters, soil moisture, trunk canopy density, trunk height and leaf canopy density were estimated over their ranges utilizing four radar data channels, which were simulated with the radiative transfer model, as inputs to the algorithm. Excellent inversion accuracy was obtained over almost the entire range of all four parameters. Successful convergence was achieved in 90 to 95% of the cases tested for the algorithm as implemented in this chapter, and it is concluded that further improvement in the convergence characteristics should be obtainable in future implementations. The results of the error analysis show that if the radar data channels are appropriately chosen so as to maximize the sensitivity of each channel to a separate inversion parameter, and the degree of discretization of the parameter space is sufficient, the overall inversion error will be minimized for a given level of uncertainty in the measured input data.

CHAPTER VIII CONCLUSIONS AND RECOMMENDATIONS 8.1 Summary This thesis provides an analysis of the application of radiative transfer (RT) theory in modeling radar scattering from vegetation canopies. Except for Chapter 7, which is concerned with model-based inversion of radar data, the specific purpose of this work has been to examine the limitations of applicability of radiative transfer theory in cases for which its fundamental underlying assumptions are not valid. A major component of this work has involved the development and validation of Monte Carlo simulations for general representative random media having structural properties that are closely related to those found in vegetation canopies. In these simulations the second-order couplings between constituent elements have been included to provide a conceptual foundation for understanding the effects of multiple scattering in shaping the radar response of dense vegetation. However, the primary purpose of constructing the Monte Carlo simulations has been to provide a benchmark with which to evaluate the performance of radiative transfer theory in domains where no other analytical theory is available for comparison. In this section we will summarize on a chapter by chapter basis the results of these investigations. In Chapter 2 the second-order radiative transfer model for a layer consisting of 200

201 vertical trunks above a dielectric ground has been developed. The first and secondorder RT models have been compared with experimental measurements made on mature corn canopies. It has been shown that second-order RT theory provides a dramatic overall improvement in the estimate of the backscattering coefficient for both the co-polarized and cross-polarized radar returns from this type of vegetation. It has also been demonstrated that the RT model does not reproduce the measured angular trend for either the vv-polarized or the cross-polarized canopy response. This is partly attributable to the presence of the Brewster angle in the RT model, which is apparently absent from the measured data, and partly to overestimation of the extinction and phase matrices by RT theory. This occurs because RT theory is unable to properly account for the decorrelation of the local fields in layers of extended scatterers and which results from multiple scattering in the canopy. Chapter 3 is concerned with the construction of a Monte Carlo scattering model for the trunk layer of a forest canopy which takes into account multiple scattering effects up to second-order. Experimental data has been presented for the purpose of validating the scattering solution for pairs of cylinders. Monte Carlo simulations based on this solution have also been presented for various representative cylinder number densities, and these results have been compared with measurements made on prepared ensembles of randomly distributed vertical cylinders. First and second-order RT solutions for the same media were also compared with the simulation results, and it has been verified that the RT model gives incorrect results for this type of medium. The failure of RT theory to predict the correct scattering behavior of media in which there is a predominance of scatterers that have dimensions large compared to the wavelength of the incident radiation is evident. In Chapter 4 a general technique based on the reciprocity theorem has been devel

202 oped for deriving the secondary scattered fields from a pair of objects. This general formulation has been applied to obtain approximate analytical expressions for the secondary scattered fields from a cylinder-sphere pair. The validity of the analytical results were verified by comparison with method of moments computations, and good agreement was obtained for both the co-polarized and cross-polarized bistatic scattering cross sections in the forward specular scattering cone. This chapter provides the basis for the construction of computational simulations of electromagnetic wave scattering from heterogeneous two-component vegetation canopies that include the effects of multiple scattering up to second-order. Chapter 5 is concerned with the construction of Monte Carlo simulations of heterogeneous two-component media consisting of large vertical cylinders and small spheres above a ground plane. The second-order interactions between cylinders and between spheres and cylinders have been included in the scattering model, and all scattering terms have been validated using the method of moments. The results of the Monte Carlo simulation have been compared with those of corresponding RT models for similar media. It has been found that the radiative transfer models do not properly predict the scattering behavior of media which are composed of scatterers having dimensions large compared to the dimensions of the medium. This is a result of the fact that the extinction matrix for the medium as computed by radiative transfer is overestimated for large particles in any plane of polarization corresponding to a large dimension and produces excessive attenuation of the scattered wave for that polarization especially at higher angles of incidence. In addition, even though the extinction matrix is not incorrectly estimated for polarizations that are not associated with large dimensions of the scatterer, the source function for the medium will still be improperly estimated in that dimension and

203 will not reproduce the proper angular trend for the scattered wave. It has also been determined that the RT model that attempts to divide the canopy into an upper layer consisting of trunks and smaller scatterers and a lower layer consisting of trunks alone can produce results that are seriously in error. Another conclusion is that while second-order interactions between cylinders can yield significant changes in the level of backscattered co-polarized radiation, the predominant effect of the interactions with the smaller spheres in the medium is to produce a change in the angular trend of the cross-polarized response. The second-order interactions between spheres and cylinders have also been shown to produce a distinct change in the co-polarized phase statistics of the backscattered wave. Chapter 6 presents a hybrid model for the computation of scattering from layered vegetation media consisting of vertical trunks above a dielectric ground plane and a crown layer consisting of small weakly scattering particles. The crown layer has been modeled using radiative transfer theory and the trunk layer has been simulated with a Monte Carlo simulation. The derivation of the transmissivity matrix for the trunk layer has been presented and some of its important features investigated. The concept of the effective scattering matrix for an average cylinder in the trunk medium has been developed and evidence has been presented that it assumes limiting behavior as the effect of multiple scattering in the medium becomes more important. Evidence has also been presented to confirm that the exponential extinction model used in radiative transfer theory works fairly well for sparse distributions of particles and even for sparse distributions of extended scatterers such as cylinders, however the radiative transfer extinction model breaks down for high densities of strongly scattering particles that are large compared with the wavelength of the excitation. In addition, it has been demonstrated that the source function integration used in radiative transfer theory

204 to account for volume scattering in random media leads to results that are severely in error for layers of long cylinders, and it may be inferred that this is also true for other types of large scatterers. Finally, several varieties of test canopies consisting of dielectric cylinders and small metallic spheres over a conducting ground plane have been investigated. It has been shown that the hybrid model developed in this work can be an effective way to compute scattering from vegetation media with this general structure. In Chapter 7 a general, model-based, iterative algorithm has been developed for use in the inversion of polarimetric radar data. The algorithm was implemented using two levels of discretization and two types of optimization and was applied to the case of a general representative vegetation canopy. This simplified test canopy was modeled using vector radiative transfer theory and consisted of vertical trunks, leaves and a rough ground surface. Four canopy parameters, soil moisture, trunk canopy density, trunk height and leaf canopy density were estimated over their ranges utilizing four radar data channels, which were simulated with the radiative transfer model, as inputs to the algorithm. Excellent inversion accuracy was obtained over almost the entire range of all four parameters. Successful convergence was achieved in 90 to 95% of the cases tested for the algorithm as implemented in this work, and it is concluded that further improvement in the convergence characteristics should be obtainable in future implementations. The results of the error analysis show that if the radar data channels are appropriately chosen so as to maximize the sensitivity of each channel to a separate inversion parameter, and the degree of discretization of the parameter space is sufficient, the overall inversion error will be minimized for a given level of uncertainty in the measured input data.

205 8.2 Recommendations for Future Work The heavy reliance on Mlonte Carlo simulations in this work has been a necessary evil. The simulations have been used for the purpose of evaluating the performance of RT theory in areas where no other analytical theories exist. Monte Carlo simulations may produce accurate results, but they can also be extremely inefficient with respect to computation time when the coupling between constituent elements of the simulated medium must be taken into account. While it has been demonstrated that first-order simulations are able, in some circumstances, to give useful results for the co-polarized backscatter coefficients, this reliance on first-order simulations imposes some considerable limitations on a model such as the hybrid model developed in this work. Future work should concentrate on developing analytical techniques for determining the ensemble average effective scattering matrix for media in which the coupling between particles is significant. This may be achieved by analytical integration of the second and higher order scattering matrices for such media over the assumed probability density functions of the constituent particles. These higher order effective scattering matrices may then be used in fast first-order type simulations to produce accurate results for dense vegetation canopies such as forests and some agricultural fields. If it is desired to be able to reproduce the scattering behavior of such vegetation, especially as far as crosspol is concerned, more work must be done on obtaining realistic analytical solutions for multiple scattering problems. Additionally, it would be useful to make a thorough experimental investigation of the extinction characteristics of distributions of different types of particles in order to formulate an accurate analytical extinction model that takes the interparticle couplings and orientational dependence of the particles into account. As far as the inversion aspect of this work is concerned, it is evident that the use of neural networks is rapidly

206 becoming the method of choice in this field. This author maintains that the best neural network available for use in inverting radar data exists between the ears of human beings and that model-based inversions can provide an additional dimension that is not available in packaged artificial neural network software. It would make an interesting study to compare the performance of the iterative inversion algorithm developed in this work with neural network based algorithms and also to see if it is possible to configure a hybrid algorithm that uses a neural network to improve the convergence properties of the model-based iterative algorithm.

APPENDICES 207

208 APPENDIX A CYLINDER MEDIUM SECOND-ORDER SCATTERING MATRIX ELEMENTS In this appendix the second order scattering matrix elements for two adjacent vertical cylinders above a ground plane are given. The elements are obtained by setting Et equal to Vi3 and Ah and taking the scalar product with 0s and Ah. First defining +o00 rTM E= (-1)nCTMeind n=-oo +00oo rTE _= (-1)nC-TEeinf n=-oo and +oo "= Z (_1)nCne n n= -oo the matrix elements for each type of interaction are found to be: L +00 i CTMRll C1 Rll 8TG =X -- (-1)m m emi Cm:m=-oo CmRi.iC1Rj

209 L +oo i TMRII Om RSGT = - E (-1) ] eim m=-oo Cm RIi -i CTERL STTG = e-iko.in/6cos(O.-$) (L - )cot 3) 7r +00 i(rTMcZM + rCm)R11 (rcTM + rTEm)RI ] m=-(oo (rTMCm + PCZrE)RI -i(rTEcrTE + JOm)Ri *H()(ko sin s3,p) eim(('-$) STGT = e-iko Sinf3pcos(4-~,),(L- cot r +0 [i(rTMCrTM + rCOm)R (rC~ZM + rTECm)RI m=-oo (rTM(mRi + rCTE)RII i(rTECER + CmRii Hm()(k0 sin p3) e im(Of-~) STGT _ -ikosinJcos(0,-d) PCOtc MM CZ~rTMCTMRII +rCmR,) (rCTMRjj +rTEcmRll -H,(n)(ko sin P-) e'm(03-~)

BIBLIOGRAPHY 210

211 BIBLIOGRAPHY [1] E.P.W. Attema and F.T. Ulaby, "Vegetation modeled as a water cloud," Radio Science, vol. 13, pp. 357-364, 1978. [2] P.P. Batlivala and F.T. Ulaby, "Feasibility of monitoring soil moisture using active microwave remote sensing," University of Kansas Remote Sensing Laboratory, Report 264-12, 1977. [3] C.F. Bohren and D.R. Huffman, Absorption and Scattering of Light by Small Particles, Wiley, New York, 1983. [4] J.J. Bowman, T.B.A. Senior and P.L.E. Uslenghi, Electromagnetic and Acoustic Scattering by Simple Shapes, New York: Hemisphere, 1987. [5] J.H. Bruning and Y.T. Lo, "Multiple scattering of EM waves by spheres, Part I - multipole expansion and ray-optical solutions," IEEE Trans. Antennas Propagat., vol. 18, pp. 378-390, 1971. [6] G.J. Burke, A.J. Poggio, Numerical Electromagnetics Code (NEC) - Method of Moments, Lawrence Livermore National Laboratory, 1981. [7] G.J. Burke, A.J. Poggio, Numerical Electromagnetics Code (NEC) - Method of Moments, Naval Ocean Systems Center, Technical Document 116, Vol. 2, 1981. [8] T.F. Bush and F.T. Ulaby, "Radar return from a continuous vegetation canopy," IEEE Trans. Antennas Propagation, vol. 24, pp. 269-276, 1976. [9] S. Chandrasekhar, Radiative Transfer, Dover, New York, 1960. [10] N.S. Chauhan, R.H. Lang and K.J. Ranson, "Radar modeling of a boreal forest," IEEE Trans. Geosci. Remote Sensing, vol. 29, pp. 627-638, 1991. [11] K.S. Chen and A.K. Fung, "An empirical bispectrum model for sea surface scattering," IGARSS Symp., vol. II, pp. 1454-1455, 1992. [12] M.C. Dobson, K. McDonald and F.T. Ulaby, "Modeling of forest canopies and analysis of polarimetric SAR data," University of Michigan Radiation Laboratory, Report 026143-1-F, 1988. [13] M.C. Dobson and F.T. Ulaby, "Preliminary evaluation of the SIR-B response to soil moisture, surface roughness, and crop canopy cover," IEEE Trans. Geosci. Remote Sensing, vol. 24, pp. 517-526, 1986.

212 [14] S.E. Decatur, "Application of neural networks to terrain classification," Proc. IJCNN, vol. 1, pp. 283-285, 1989. [15] S.L. Durden, J.J. Van Zyl and H.A. Zebker, "Modeling and observation of radar polarization signatures of forested areas," IEEE Trans. Geosci. Remote Sensing, vol. 27, pp. 290-301, 1989. [16] N. Engheta and C. Elachi, "Radar scattering from a diffuse vegetation layer over a smooth surface," IEEE Trans. Geosci. Remote Sensing, vol. 20, pp. 212-216, 1982. [17] A.W. England, "Thermal microwave emission from a half-space containing scatterers," Radio Science, vol. 9, pp. 447-454, 1974. [18] A.W. England, "Thermal microwave emission from a scattering layer," J. Geophys. Res., vol. 80, pp. 4484-4496, 1975. [19] H.J. Eom and A.K. Fung, "A scatter model for vegetation up to KU-band," Remote Sensing of Environment, vol. 15, pp. 185-200, 1984. [20] L. Foldy, "The multiple scattering of waves," Phys. Rev., vol. 67, pp. 107-119, 1945. [21] U. Frisch, "Wave propagation in random media," in Probabilistic Methods in Applied Mathematics, vol. 1, Academic Press, New York, 1968. [22] A.K. Fung and H.S. Fung, "Application of the first-order renormalization method to scattering from a vegetated-like half-space," IEEE Trans. Geosci. Electr., vol. 15, pp. 189-195, 1977. [23] A.K. Fung, "A review of volume scattering theories for modeling applications," Radio Science, vol. 17, pp. 1007-1017, 1982. [24] A.K. Fung and F.T. Ulaby, "A scatter model for leafy vegetation," IEEE Trans. Geosci. Remote Sensing, vol. 16, pp. 281-286, 1978. [25] N.S. Goel, "Modeling canopy reflectance and microwave backscattering coefficient," Remote Sensing of Environment, vol. 18, pp. 235-253, 1985. [26] N.S. Goel and T. Grier, "Estimation of canopy parameters for inhomogeneous vegetation canopies from reflectance data. I. Two-dimensional row canopy," Int. J. Remote Sensing, vol. 7, pp. 665-681, 1986. [27] N.S. Goel and R.L. Thompson, "Inversion of vegetation canopy reflectance models for estimating agronomic variables. V. Estimate of leaf area index and average leaf angle using measured canopy reflectances," Remote Sensing of Environment, vol. 16, pp. 69-85, 1984. [28] S. Gogineni, A.K. Fung,.S. Chen and J. Wang, "Comparison of measurements and theory for backscatter from vegetation-covered soil on the konza prairie," Proceedings IGARSS'92, Houston, Texas, vol. II, pp. 920-922, 1992.

213 [29] M.T. Hallikainen, F.T. Ulaby, M.C. Dobson, M.A. El-Rayes and L.K. Wu, "Microwave dielectric behavior of wet soil - Part I: Empirical models and experimental observations," IEEE Trans. Geosci. Remote Sensing, vol. 23, pp. 25-34, 1985. [30] R.F. Harrington, Time-Harmonic Electromagnetic Fields, New York: McGrawHill, 1961. [31] H. Hirosawa, H. Ishida, T. Ochi and Y. Matsuzaka, "Measurement of microwave backscatter from trees," Presented at the 3rd International Colloquium on Special Signatures of Objects in Remote Sensing held in Les Arcs, France on 16-20 December, 1985. [32] D.H. Hoekman, "Radar backscatter from forest stands," International Journal of Remote Sensing, vol. 6, pp. 325-343, 1985. [33] S. Hong, K. Fukue, H. Shimoda and T. Sakata, "Non-parametric texture extraction using neural network," Proceedings IGARSS'92, Houston, Texas, vol. II, pp. 1084-1086, 1992. [34] D. Holliday, G. St-Cyr and N.E. Woods, "Comparison of a new radar ocean imaging model with sarsex internal wave data," Int. J. Remote Sensing, vol. 8, pp. 1423-1430, 1987. [35] Y. Hussin, R. Reich and R. Hoffer, "Estimating slash pine biomass using radar backscatter," IEEE Trans. Geosci. Remote Sensing, vol. 29, pp. 427-431, 1991. [36] A. Ishimaru, Wave Propagation and Scattering in Random Media, vol. 1, Academic Press, New York, 1978. [37] A. Ishimaru, Wave Propagation and Scattering in Random Media, vol. 2, Academic Press, New York, 1978. [38] A. Ishimaru, R.J. Marks II, L. Tsang, C.M. Lam and D.C. Park, "Particle-size distribution using optical sensing and neural networks," Opt. Lett., vol. 15, no. 21, 1990. [39] M.A. Karam and A.K. Fung, "Electromagnetic scattering from a layer of finite length, randomly oriented, dielectric circular cylinders over a rough interface with application to vegetation," International Journal of Remote Sensing, vol. 9, pp. 1109-1134, 1988. [40] M.A. Karam and A.K. Fung, "Scattering from randomly oriented circular discs with application to vegetation," Radio Science, vol. 18, pp. 557-565, 1983. [41] M.A. Karam and A.K. Fung, "EM Scattering from a randomly oriented circular dielectric, finite-length cylinder," International Union Radio Science Commission F: Wave Propagation and Remote Sensing, University of New Hampshire, Durham, New Hampshire, 4.1.1-4.1.3.

214 [42] M.A. IKaram, A.K. Fung and Y.M. Antar, "Scattering models for vegetation samples," Proceeding of IEEE Geosci. Remote Sensing Symposium, vol. 2, pp. 1013-1018, 1987. [43] M.A. Karam, A.K. Fung, R.H. Lang and N.S. Chauhan, "A microwave scattering model for layered vegetation," IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 767-784, 1992. [44] R.H. Lang, "Electromagnetic backscattering from a sparse distribution of lossy dielectric scatterers," Radio Science, vol. 16, pp. 15-30, 1981. [45] R.H. Lang and J.S. Sidhu, "Electromagnetic backscattering from a layer of vegetation: a discrete approach," IEEE Trans. Geosci. Remote Sensing, vol. 21, pp. 62-71, 1983. [46] J.K. Lee and J.A. Kong, "Active microwave remote sensing of an anisotropic random medium layer," IEEE Trans. Geosci. Remote Sensing, vol. 23, pp. 910923, 1985. [47] T. Le Toan, A. Beaudoin, J. Riom and D. Guyon, "Relating forest biomass to SAR data," IEEE Trans. Geosci. Remote Sensing, vol. 30, March 1992. [48] T. LeToan, H. Laur, E. Mougin and A. Lopes, "Multitemporal and dualpolarization observations of agricultural vegetation covers by X-band SAR images," IEEE Trans. Geosci. Remote Sensing, vol. 27, pp. 709-718, 1989. [49] C. Liang and Y.T. Lo, "Scattering by two spheres," Radio Science, vol. 2, pp. 1481-1495, 1967. [50] I.V. Lindell and E. Alanen, "Exact image theory for the Sommerfeld half-space problem, Part III: General formulation," IEEE Trans. Antennas Propagat., vol. 32, pp. 1027-1032, 1984. [51] C. Matzler and E. Schanda, "Snow mapping with active microwave sensors," Int. J. Remote Sensing, vol. 5, pp. 409-422, 1984. [52] G.E. McClelland, R.N. DeWitt, T.H. Hemmer, M.L. Matheson and G.O. Moe, "Multispectral image processing with a three-layer backpropagation network," Proc. IJCNN, vol. 1, pp. 151-153, 1989. [53] K.C. McDonald, M.C. Dobson and F.T. Ulaby, "Using MIMICS to model Lband multi-angle and multi-temporal backscatter from a walnut orchard," IEEE Trans. Geosci. Remote Sensing, vol. 28, pp. 477-491, 1990. [54] K.C. McDonald and F.T. Ulaby, "MIMICS II: Radiative Transfer Modeling of Discontinuous Tree Canopies at Microwave Frequencies," IGARSS'90, pp. 20-24, 1990, College Park, Maryland. [55] K.C. McDonald and F.T. Ulaby, "Modeling microwave backscatter from discontinuous tree canopies," University of Michigan Radiation Laboratory, Technical Report 026511-2-T, 1991.

215 [56] K.C. McDonald and F.T. Ulaby, "Modeling microwave backscatter from discontinuous tree canopies," University of Michigan Radiation Laboratory, Report 026511-2-T, 1991. [57] Y. Oh, K. Sarabandi and F.T. Ulaby, "An empirical model and an inversion technique for radar scattering from bare soil surfaces," IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 370-381, 1992. [58] W.H. Peake, "Interaction of electromagnetic waves with some natural surfaces," IRE Trans. Antennas Propagation, S324-S329, 1959. [59] L. Pierce, K. Sarabandi and F.T. Ulaby, "Application of an artificial neural network in canopy scattering inversion," Proceedings IGARSS'92, Houston, Texas, vol. II, pp. 1067-1069, 1992. [60] D.A. Pierre, Optimization Theory with Applications, New York: Dover, 1986. [61] P.F. Polatin, Unpublished Data. [62] P.F. Polatin, K. Sarabandi and F.T. Ulaby, "Monte Carlo simulation of electromagnetic scattering from a heterogeneous two-component medium," PIERS 1993 Proceedings, Pasadena, CA, pp. 870, July 12-16, 1993. [63] J.A. Richards, G.Q. Sun and D.S. Simonett, "L-band radar backscatter modeling of forest stands," IEEE Trans. Geosci. Remote Sensing, vol. 25, pp. 487-498, 1987. [64] G.T. Ruck, D.E. Barrick, W.D. Stuart and C.K. Krichbaum, Radar CrossSection Handbook. New York: Plenum, 1970. [65] D.E. Rumelhart, J.L. McClelland and the PDP Research Group, Parallel Distributed Processing, vol. I - Foundations, Cambridge: MIT Press, 1986. [66] K. Sarabandi, "Derivation of phase statistics of distributed targets from the averaged mueller matrix," to be published in Radio Science, 1991. [67] K. Sarabandi, "Scattering from dielectric structures above impedance surfaces and resistive sheets," IEEE Trans. Antennas Propagat., vol. 40, pp. 67-78, 1992. [68] K. Sarabandi, P.F. Polatin, "Electromagnetic scattering from two adjacent objects," submitted for publication, IEEE Trans. Antennas Propagat.. [69] K. Sarabandi, P.F. Polatin, F.T. Ulaby, "Monte carlo simulation of scattering from a layer of vertical cylinders," IEEE Trans. Antennas Propagat., vol. 41, April 1993. [70] K. Sarabandi and T.B.A. Senior, "Low frequency scattering from cylindrical structures at oblique incidence," IEEE Trans. Geosci. Remote Sensing, vol. 28, pp. 879-885, 1990.

216 [71] K. Sarabandi, T.B.A. Senior and F.T. Ulaby, "Effect of curvature on the backscattering from leaves," J. Electromagnetic Waves and Applications, vol. 2, pp. 653-670, 1988. [72] S. Seker, "Microwave backscattering from a layer of randomly oriented discs with application to scattering from vegetation," I.E.E. Proceedings, Microwave, Antennas and Propagation, H, pp. 497-502, 1986. [73] S. Seker and A. Schneider, "Electromagnetic scattering from finite-length dielectric cylinders," Digest of International Geoscience and Remote Sensing Symposium, vol. 1, pp. 523-526, 1985. [74] T.B.A. Senior, K. Sarabandi and F.T. Ulaby, "Measuring and modeling the backscattering cross section of a leaf," Radio Science, vol. 22, pp. 1109-1116, 1987. [75] J. Shi, J.V. Soares, L. Hess, E.T. Engman and J. van Zyl, "SAR-derived soil moisture measurements for bare fields," Proceedings IGARSS'91, Espoo, Finland, vol. I, pp. 393-396, 1991. [76] H. Shinohara, T. Homma, H. Nohmi, H. Hirosawa and T. Tagawa, "Relation between L-band microwave penetration/ backscattering characteristics and state of trees," Proceedings IGARSS'92, Houston, Texas, vol. I, pp. 539-541, 1992. [77] H.S. Tan and A.K. Fung, "A first-order theory on wave depolarization by a geometrically anisotropic random medium," Radio Science, vol. 14, pp. 377386, 1979. [78] H.S. Tan, A.K. Fung and H.J. Eom, "A second order renormalization theory for cross-polarized backscatter from a half-space random medium," Radio Science, vol. 15, pp. 1059-1065, 1980. [79] L. Tsang, Z. Chen, S. Oh and R.J. Marks II, "Inversion of snow parameters from passive microwave remote sensing measurements by a neural network trained with a multiple scattering model," IGARSS'91, Espoo, Finland. [80] L. Tsang and J.A. Kong, "Thermal microwave emission from half-space medium," Radio Science, vol. 11, pp. 599-609, 1976. [81] L. Tsang and J. Kong, "Application of strong fluctuation random medium theory to scattering from a vegetation-like half space," IEEE Trans. Geosci. Remote Sensing, vol. 19, pp. 62-69, 1981. [82] L. Tsang, M.C. Kubasci and J.A. Kong, "Radiative transfer theory for active remote sensing of a layer of small ellipsoidal scatterers," Radio Science, vol. 16, pp. 321-329, 1981. [83] L. Tsang and J.A. Kong, "Radiative transfer theory for active remote sensing of half-space random media," Radio Science, vol. 13, pp. 763-773, 1978.

217 [84] L. Tsang, J.A. Kong and R.T. Shin, Theory of Microwave Remote Sensing. New York: Wiley-Interscience, 1985. [85] V. Twersky, "On propagation in random media of discrete scatterers," Proc. Symp. Appl. Math., vol. 16, pp. 84-116, 1964. [86] V. Twersky, "Interference effects in multiple scattering by large, low-refracting, absorbing particles," J. Opt. Soc. Am., vol. 60, pp. 908-914, 1970. [87] S. Twomey, H. Jacobowitz and H.B. Howell, "Matrix methods for multiple scattering problems," J. Atmos. Sci., vol. 23, pp. 289-296, 1966. [88] F.T. Ulaby and T.F. Bush, "Monitoring wheat growth with radar," Photogramm. Eng. Remote Sensing, vol. 42, pp. 557-568, 1976. [89] F.T. Ulaby and T.F. Bush, "Corn growth as monitored by radar," IEEE Trans. Antennas Propagation, vol. 24, pp. 819-828, 1976. [90] Fawwaz T. Ulaby and Charles Elachi, Radar Polarimetry for Geoscience Applications, Norwood, Massachusetts: Artech House, 1990, pp. 95. [91] F.T. Ulaby and M.A. El-Rayes, "Microwave dielectric spectrum of vegetation, Part II: Dual-dispersion model," IEEE Trans. Geosci. Remote Sensing, vol. 25, pp. 550-557, 1987. [92] F.T Ulaby, R.K. Moore and A.K. Fung, Microwave Remote Sensing: Active and Passive, Vol. II - Radar Remote Sensing and Surface Scattering and Emission Theory., Reading, Massachusetts: Addison-Wesley, 1982. [93] F.T Ulaby, R.K. Moore and A.K. Fung, Microwave Remote Sensing: Active and Passive, Vol. III - From Theory to Applications., Dedham, Massachusetts: Artech House, 1986. [94] F.T. Ulaby, D. Held, M.C. Dobson, K. McDonald and T.B.A. Senior, "Relating polarization phase difference of SAR signals to scene properties," IEEE Trans. Geosci. Remote Sensing, vol. 25, pp. 83-92, 1987. [95] F.T. Ulaby, K. Sarabandi, K. McDonald, M. Whitt and M.C. Dobson, "Michigan Microwave Canopy Scattering Model (MIMICS)," University of Michigan Radiation Laboratory Report 022486-3-T, 1988. [96] F.T. Ulaby, K. Sarabandi, K. McDonald, M. Whitt and M.C. Dobson, "Michigan microwave canopy scattering model," Int J. Remote Sensing, vol. 11, pp. 1223-1253, 1990. [97] F.T. Ulaby, A. Tavakoli and T.B.A. Senior, "Microwave propagation constant for a vegetation canopy with vertical stalks," IEEE Trans. Geocsi. Remote Sensing, vol. 25, pp. 714-725, 1987.

218 [98] F.T. Ulaby, M.W. Whitt and M.C. Dobson, "Measuring the propagation properties of a forest canopy using a polarimetric scatterometer," IEEE Trans. Antennas and Propagation, vol. 38, pp. 251-258, 1990. [99] F.T. Ulaby and W.H. Stiles, "The active and passive response to snow parameters, part II: Water equivalent of dry snow," Journal of Geophysical Research, vol. 85, pp. 1045-1049, 1980. [100] H.C. Van de Hulst, Light Scattering by Small Particles, New York: John Wiley, 1957. [101] J.J. Van Zyl, H.A. Zebker and C. Elachi, "Imaging radar polarization signatures: theory and observations," Radio Science, vol. 22, pp. 529-543, 1987. [102] J.R. Wang, E.T. Engman, T.J. Schmugge, T. Mo and J.C. Shiue, "The effects of soil moisture, surface roughness, and vegetation on L-band emission and backscatter," IEEE Trans. Geosci. Remote Sensing, vol. 25, pp. 825-833, 1987. [103] C. Wessman, J. Aber and D. Peterson, "An evaluation of imaging spectrometry for estimating forest canopy chemistry," Int J. Remote Sensing, vol. 10, pp. 1293-1316, 1989. [104] M.W. Whitt, Unpublished Data. [105] M.W. Whitt and F.T. Ulaby, "Microwave scattering from periodic rowstructured vegetation," University of Michigan Radiation Laboratory, Report 026511-3-T, 1991. [106] E. Wilcox and M.C. Dobson, "Leaf area index data for the Michigan forest test sites," University of Michigan Radiation Laboratory, Report 025761-1-T, 1992. [107] H.A. Zebker, J.J Van Zyl and D.N. Held, "Imaging radar polarimetry from wave synthesis," Journal of Geophysical Research, vol. 92, pp. 683-701, 1987. [108] M. Zuniga, T. Habashy and J. Kong, "Active remote sensing of layered random media," IEEE Trans. Geosci. Electron., vol. 17, pp. 296-302, 1979.

UNIVERSITY OF MICHIGAN Ii 90111[111115 02827 289911111[1 3 9015 02827 2899