A FRACTAL-BASED COHERENT SCATTERING AND PROPAGATION MODEL FOR FOREST CANOPIES by Yi-Cheng Lin A dissertation submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy (Electrical Engineering) in The University of Michigan 1997 Doctoral Committee: Associate Professor Kamal Sarabandi, Chair Professor Fawwaz T. Ulaby Professor Anthony W. England Professor Thomas B. A. Senior Professor William R. Kuhn RL-954 = RL-954

" Let the land produce vegetation: seed bearing plants and trees on the land that bear fruit with seed in it, according to their various kinds. " Genesis 1:1 1

ABSTRACT A FRACTAL-BASED COHERENT SCATTERING AND PROPAGATION MODEL FOR FOREST CANOPIES by Yi-Cheng Lin Chair: Kamal Sarabandi Dependency of tree structures on SAR backscatter and modeling of radar interferometry on forestry are two important topics in the field of microwave remote sensing. However, traditional scattering models for forest canopies based on radiative transfer theory do not preserve the coherent effect caused by the relative position of scatterers within a tree structure and cannot provide the phase information of the backscattered fields. To address these deficiencies, a fractal-based coherent scattering and propagation model is developed in this dissertation. The model comprises three components: 1) generation of realistic tree structures, 2) modeling the interaction of the electromagnetic wave with vegetation components, and 3) Monte Carlo simulation. Generation of realistic tree architectures is implemented by a fractal-based stochastic Lindenmayer system. The electromagnetic problem is formulated by considering the generated tree structure as a cluster of scatterers composed of dielectric cylinders and disks with specified size, position, and

1 orientation. Using the single scattering approximation, the total scattered field is obtained from the coherent addition of the individual scattering from each scatterer illuminated by a mean field. Foldy's approximation is invoked to model the coherent wave propagation within the forest canopy where the mean field at a given point accounts for the accumulated attenuation and phase change caused by vegetation components. Finally, the desired statistics of the scattered field are acquired using a Monte Carlo simulation over a large number of realizations. Combined with the recently developed Ak-radar equivalence algorithm, the model can numerically predict the scattering phase center height of forest canopies, the very quantity measured by an interferometric SAR (INSAR). The accuracy of the model is verified using SIR-C (SAR) and TOPSAR (INSAR) measurements over test stands of coniferous and deciduous types near Raco, Michigan. A sensitivity analysis is performed to characterize the response of the forest canopy with controlled physical parameters to a multi-parameter (frequency, polarization, and incidence angle) SAR/INSAR. To improve the model utility and efficiency, a procedure for developing simplified empirical formulation for scattering from complex forest structures is presented. Finally, an inverse scattering model based on genetic algorithms is developed for retrieval of forest parameters that may take multi-parameter SAR/INSAR data as its input.

o Yi-Cheng Lin 1997 All Rights Reserved

To my parents for their support and encouragement and my wife for her endurance and love. 11

ACKNOWLEDGEMENTS I gratefully acknowledge my dissertation committee for their help and advice. Special thanks are due to Professor Kamal Sarabandi for his guidance and financial support throughout the course of this work. I would also like to thank the following friends and colleagues in Radiation Laboratory for their friendship and encouragement during my graduate studies: Professor Chen-To Tai, Drs. Valdis Liepa, Leland Pierce, Craig Dobson, Adib Nashashibi, DeRoo Roger, James Stiles, Paul Siqueira, Chen-Yu Chi, Jian Gong, John-Gwan Yook, Kazem Sabet, Messrs. Ron Hartikka, Jui-Ching Cheng, Tsen-Chieh Chiu, and Eric Li. My grateful heart goes to God, the amazing creator, who prepared people to help me to know, reconcile with, and live within Him during my graduate studies in Ann Arbor. I wish to thank Mr. Chris Hardy at International Student Institute for his time and patience in Bible study/discussion and Pastor Fu at Ann Arbor Chinese Christian Church for his help and care. I also appreciate for the sharing and encouragement that I received from my spiritual brothers and sisters at AACCC, especially from Drs. Ming-Jen Tseng, Song-Ping Cheng, Chiu-Fong Lin, I-Ming Tso, Jih-Tsung Kuo, Jih-Liang Tseng, Joe Juan, and their families. Most importantly, I wish to express my deep gratitude to my family members: my wife, Chwen-Huei Liou, for her endurance and inspiration, my parents, Shan-Jin Lin and Su-Yen Hsu, for their selfless support, and my parents-in-law as well as my elder brothers and sister for their constant encouragement. iii

TABLE OF CONTENTS DEDICATION.................................. ii ACKNOWLEDGEMENTS.......................... iii LIST OF FIGURES............................... vii LIST OF TABLES................................ xii CHAPTER I. INTRODUCTION......................... 1 1.1 Background........................... 1 1.1.1 Forestry and Remote Sensing.............. 1 1.1.2 SAR/INSAR Technologies.............. 2 1.1.3 SAR/INSAR in Forest Remote Sensing....... 4 1.1.4 Existing Scattering Models for Forest Canopies... 5 1.2 Goals............................... 8 1.3 Dissertation Overview...................... 9 II. FRACTAL MODEL FOR GENERATION OF TREE STRUCT U R E S................................. 11 2.1 Fractal Theory and Lindenmayer Systems........... 11 2.2 Modeling of Tree Structures................... 14 2.3 Computer Implementation.................... 21 2.4 Application to Remote Sensing of Vegetation......... 22 III. ELECTROMAGNETIC SCATTERING MODEL FOR A TREE TRUNK ABOVE A TILTED GROUND PLANE....... 28 3.1 Introduction........................... 28 3.2 Global Coordinate System and Notations........... 30 3.3 Polarimetric Reflection from A Tilted Ground Plane..... 34 3.4 Semi-Exact Solution for Electromagnetic Scattering from A Finite Stratified Dielectric Cylinder............... 36 iv

3.5 Physical Optics Approximation................. 41 3.6 Modeling of a Corrugated Bark Layer................ 43 3.7 Simulation Results........................45 3.8 Conclusions............................ 53 IV. COHERENT SCATTERING AND PROPAGATION MODEL FOR CONIFEROUS AND DECIDUOUS CANOPIES.... 59 4.1 Introduction.......................... 59 4.2 Electromagnetic Scattering from A Tree Structure...... 61 4.3 Coherent Wave Propagation within Tree Canopies...... 65 4.3.1 Deciduous Canopies.................... 68 4.3.2 Coniferous Canopies.................... 71 4.4 Monte Carlo Simulation....................... 73 4.5 Model Validation......................... 74 4.5.1 SIR-C/X-SAR....................... 75 4.5.2 Ground Truth..................... 75 4.5.3 Simulation Results................... 76 4.6 Conclusions............................. 91 V. SIMULATION OF INTERFEROMETRIC SAR RESPONSE FOR CHARACTERIZING THE SCATTERING PHASE CENTER STATISTICS OF FOREST CANOPIES......... 92 5.1 Introduction.......................... 92 5.2 Model Description.......................... 95 5.3 Verification of Identity Between Ak-radar and INSAR System 97 5.3.1 INSAR Case...................... 98 5.3.2 Ak-Radar Case...................... 101 5.4 Scattering Mechanisms vs. Scattering Phase Centers...... 103 5.5 Model Validation for INSAR Applications.......... 06 5.6 Sensitivity Study................... 10 5.7 Conclusions............................ 1.22 VI. RETRIEVAL OF FOREST PARAMETERS USING A FRACTALBASED COHERENT SCATTERING MODEL AND GENETIC ALGORITHMS......................... 124 6.1 Introduction........................... 124 6.2 Empirical Model Development.................... 128 6.3 Inversion Algorithms........................ 136 6.3.1 Least-Square Approach................ 139 6.3.2 Genetic Algorithms................... 140 6.4 Conclusions............................. 146 v

VII. CONCLUSIONS AND RECOMMENDATIONS........ 148 7.1 Summary of Achievements................... 148 7.2 Recommendations for Future Work................. L50 BIBLIOGRAPHY................................. 152 vi

LIST OF FIGURES Figure 2.1 An example of the growing process of a fractal tree................13 2.2 3-D branching and branch tapering at a splitting node.............15 2.3 A fractal tree structure of deciduous type...................... 17 2.4 A fractal tree structure of coniferous type............................ 18 2.5 Leaf attachment schemes for the red pine: the end cluster, the stem cluster, and the leaflet on the leaf bud along three concentric spirals. 19 2.6 A red pine with attached needles............................... 20 2.7 The top and side views of a deciduous tree with attached leaves.. 25 2.8 The photograph of a red pine stand (Stand 22), and the simulated tree structure using the fractal model........................... 26 2.9 The photograph of a red maple stand (Stand 31), and the simulated tree structure using the fractal model....................... 27 3.1 Global coordinate system and notations......................... 31 3.2 The first-order scattering components for scattering from an object above a ground plane.................................... 33 3.3 Geometry and composition of a finite stratified cylinder (top-view and side-view).......................................... 37 3.4 A corrugated layer (a) and its equivalent anisotropic layer (b).... 44 vii

3.5 Comparison of the PO approximation with the semi-exact solution. The ratio of the interior radius a2 and the exterior radius al(= a) is kept constant (a2/al = 0.9). Other parameters are: b = 20Ao, 1 = 4 + il, E2 = eg = 10 + i5, Oi = 120~, 'i = 180~, 0, = 60~, s = 0~, g = 0c = 0~................................ 47 3.6 The normalized backscattering RCS as a function of the incidence angle 0 = 7- 0i for a two-layer cylinder with and without the corrugation. Other parameters are: b = 20Ao, al = 2o, a2 = 1.8A,t = 0.1o, d/L = 0.7, = 1800, q = 0~, = Oc = 0............ 48 3.7 The normalized bistatic a,, (a) and Chh (b) as a function of the scattering elevation angle 0, in X-Z plane (0, > 0 when Os = 0~; 0, < 0 when 4s = 180~). The backscattered and the specular directions are shown at 0s = 60~ and O, = -60~ respectively. Other parameters are: Oi = 120~, i = 180~, g = O 00, b = 10Xo, al = 2Ao, a2 = 1.8Ao. 49 3.8 The normalized bistatic av (a) and ghh (b) as a function of the scattering azimuth angle q,. The backscattered and forward-scattered directions are shown at,s = 0~ and 4, = 180~ respectively. Other parameters are: 0i = 120~,4i = 180~,Os = 60~,Og = 0 = 0~,b = 200o, al = 4Ao, a2 = 3.6Ao.......................... 51 3.9 The normalized backscattering RCS as a function of the incidence angle r - Oi for a two-layer cylinder above a tilted ground plane. Other parameters are: g = 20~, g = 90~, = 0~, = 180~, = 0~, 0 = - Oi, b = 20Aoal = 2Ao,a2 = 1.8Ao............... 55 3.10 The normalized backscattering RCS as a function of the ground azimuth angle g. Other parameters are: g = 20~,0c = 0~,0 = 135~, 4i = 180~, 0s = 45~, Oqs = 0~, b = 20Xo, al = 2A, a2 = 1.8Ao.... 56 3.11 The normalized backscattering RCS as a function of the ground elevation angle 0g (Og > 0 when cg = 0~;O0g < 0 when 14 = 180~). Other parameters are: 0i = 135~, i = 180~, 0, = 45~,,s = 0~, b = 10o a = 2o, a2 = 1.8A0............................ 57 3.12 The geometry of the scattering configurations for a cylinder over a tilted ground where a strong backscatter can be observed...... 58 4.1 Scattering from a cluster of scatterers above a tilted ground plane.. 62 4.2 Four scattering components from an object above a tilted ground plane.................................. 63 viii * * Vlll

4.3 Extinction of a coherent wave in random media.......... 67 4.4 The extinction profile calculated for Stand 31................. 70 4.5 The position dependence of the transmissivity for coniferous trees and the shadow effect caused by neighboring trees. 72 4.6 Convergence behavior of the Monte Carlo simulation for Stand 31 at L-band(a) and C-band (b) at incidence angle i = 30~. 78 4.7 Comparison between the model predictions (lines) and SIR-C data (symbols) at (a) L-band and (b) C-band..................... 80 4.8 Contribution of different scattering components (of~ and crt) to the overall backscattering coefficient (r^0otal) as a function of incidence angle.................................. 81 4.9 Comparison of the 11-layer and 2-layer extinction models simulated at C-band for hh polarization (a) and vv polarization (b)....... 83 4.10 Histogram of the backscattering coefficients for Stand 31 at Oi = 43.6~. 84 4.11 One-way transmissivity as a function of leaf area index (LAI).... 85 4.12 Histogram of branch orientation resulted from two branching angles: (a) A0b = 22~ ~ 5~ (Stand 31); (b)AOb = 15~ ~ 3~ (Case 1). 86 4.13 Tree structures generated for (a) Stand 31, (b) Case 1, and (c) Case 2.................................... 8 7 4.14 The effect of the ground tilt angle on radar backscatter: (a) backscattering coefficients for Stand 31 over a tilted ground with 0g = 10~ simulated at L-band with incidence angle Oi = 25.4~, (b) a vertical cylinder on a tilted dielectric plane and the associated coordinates, (c) RCS simulation for (b).................................... 89 4.15 Comparison of the polarimetric backscattering coefficients at different frequency bands (L, C, and X)......................... 90 5.1 Multipath scattering mechanisms of an object (at C with projection at 0 and image at Ci) above a ground plane observed by an interferometer with two antennas at Al and A2..................98 ix

5.2 An object above a ground plane illuminated by a Ak-radar......:101 5.3 Four configurations for a cylinder above a ground plane with orientation angles: (a) Oc = 60~, c = 180~, (b) Oc = 00~,c = 0~, (c) 0C = 60~, c = 0~, (d) Oc = 45~, 4c = 150~, and their principal scattering mechanisms respectively. The center of the scatterer is denoted by (,) and equivalent scattering phase center by (D).............104 5.4 A portion of a TOPSAR C-band image (aO ), indicating Stand 22 (red pine) at an airport near Raco, Michigan (April 1995)...... 107 5.5 The estimated height of scattering phase center of Stand 22 (red pine), compared with the data extracted from two TOPSAR images (0 = 39~ and 53~) of the same stand................... 109 5.6 The simulated backscattering coefficient of Stand 22 (red pine), compared with the measured 4o extracted from two TOPSAR images of the same stand.........................111 5.7 The estimated scattering phase center height of Stand 22 (red pine) as a function of soil moisture, simulated at Oi = 45~........ 113 5.8 The estimated scattering phase center height of Stand 22 (red pine) as a function of tree density, simulated at Qi = 45~.......... 114 5.9 The estimated scattering phase center height of Stand 31 (red maple) as a function of incidence angle, with fully-polarimetric L- and C-and response.................................. 16 5.10 The estimated scattering phase center height of Stand 31 (red maple) over a tilted ground with tilt angle 0g = 10~, at Oi = 25.40...... 117 5.11 PDF of the scattering phase center height of Stand 22 (red pine) at 6i = 45~ as a function of frequency and polarization.......... 119 5.12 The correlation coefficient as a function of the frequency shift, simulated from Stand 22 (red pine)..................... 120 5.13 The correlation coefficient as a function of the frequency shift, simulated from Stand 31 (red maple)...................121 x

6.1 The sensitivity analysis of the C-band polarimetric scattering phase center height as a function of the physical parameters: (a) trunk diameter, (b) tree height, (c) tree density, (d) branching angle, (e) soil moisture, and (f) wood moisture, simulated at incidence angle 0 = 25.................................... 1131 6.2 The sensitivity analysis of the C-band polarimetric backscattering coefficient as a function of the physical parameters: (a) trunk diameter, (b) tree height, (c) tree density, (d) branching angle, (e) soil moisture, and (f) wood moisture, simulated at incidence angle 0 = 25~. 132 6.3 The angular dependence of the polarimetric backscatter in terms of: (a) the scattering phase center height Ze and (b) the backscattering coefficient a~. The simulation results are fitted with polynomials of degree 3..............................137 6.4 A comparison between the empirical model and the Monte Carlo coherent scattering model for a red pine stand at vv polarization configuration............................... 138 6.5 A comparison between the empirical model and the Monte Carlo coherent scattering model for a red pine stand at vh polarization configuration................................ 138 6.6 A comparison between the empirical model and the Monte Carlo coherent scattering model for a red pine stand at hh polarization configuration.............................. 138 6.7 A flow chart of a genetic algorithm................... 1141 6.8 Comparison of the input parameters (x) and the output of the inversion algorithm (x') using the synthetic data obtained from the Monte Carlo simulation, for (a) trunk diameter d,, (b) tree height Ht, (c) tree density Dt, (d) branching angle 0b, (e) soil moisture m, and (f) vegetation moisture m,. Here r7 is a measure of the average error in the inversion process defined by equation(6.15).......... 144 6.9 Comparison of the input parameters (x) and the output of the inversion algorithm (x') using the synthetic data obtained from the empirical model, for (a) trunk diameter d,, (b) tree height Ht, (c) tree density Dt, (d) branching angle 0b, (e) soil moisture m,, and (f) vegetation moisture mo. Here r is a measure of the average error in the inversion process defined by equation (6.15)............ 145 xi

LIST OF TABLES Table 4.1 Ground Truth of Stand 31.............................. 76 4.2 Dielectric properties of Stand 31........................ 77 4.3 Effect of tree structures on backscattering coefficients, simulated at 0i = 43.6~................................... 87 5.1 The normalized height of the scattering phase center, the normalized scattering components, and the radar cross section for four different scattering configurations as shown in Figure 5.3.......105 5.2 Ground truth data of Stand 22 (red pine) collected on April 26, 1995.108 6.1 The ranges of the selected ground truth parameters and the corresponding percentage variations to the centroid.............. 130 6.2 The inversion results using the interferometric TOPSAR data as the measured channels M. Here x is the actual ground truth data and x' is the output of the inversion process................. 146 xii

CHAPTER I INTRODUCTION 1.1 Background 1.1.1 Forestry and Remote Sensing Forests, constituting about 90% of the standing biomass in the biosphere, play a crucial role in the Earth's systems including the hydrologic cycle, the carbon dioxide circulation, and the energy balance [1, 2]. Since carbon dioxide absorbs radiant heat (infra-red energy) retained in the atmosphere, forest ecosystems may indirectly affect global climate change. Furthermore, forests are valuable natural resources to human beings and provide a habitat for a wide variety of species. However, forests worldwide are being destroyed or reduced by logging, fires, floods, and pollution. It is a global vision that the preservation of these important resources is of interest and requires more endeavors in terms of comprehensive investigation, timely monitoring, and effective management. To achieve these goals, the first task is to develop a remotesensing technology that can probe the forest canopies continuously and globally and provide accurate, reliable, and complete information about the physical parameters. Our knowledge of the vegetation biomass in the biosphere is quite limited clue to the difficulty in obtaining sufficiently high quality observations that are representative of a region or a ecosystem type. Measurements on the ground are very 1

2 time-consuming, labor-intensive, and often constrained by lack of access. Satellite remote sensing, however, provides a means to collect the geophysical information of the Earth's surface repeatedly on a regional or global scale. Intense research efforts have been conducted mainly with optical and microwave remote sensing instruments [3,4]. The disadvantages of optical sensors are that cloud-free conditions and sunlight illumination are required. It turns out that optical sensors are restricted to day-time use only and may be ineffectual in the tropics where clouds and rain often occur all year. In addition, optical remotely sensed data cannot work for most terrestrial biomass densities, because there is a saturation effect at very low levels of biomass. Microwave remote sensing systems can perform measurements independent of both weather conditions and the time of day. There are two measuring schemes for microwave remote sensing technology: active and passive [5,6]. Radiometers, known as passive sensors, measure the brightness temperature through the thermal radiation from the observed scenes. Imaging radars, on the other hand, actively assess the backscattering properties of the distributed target by transmitting radiowave radiation from an antenna and receiving the returned signal from the target. The disadvantage of this real-aperture radar system is its poor resolution which is proportional to the measuring distance and the antenna directivity. This deficiency can be solved by using the synthetic aperture technique discussed next. 1.1.2 SAR/INSAR Technologies Synthetic Aperture Radar (SAR) is a remote-sensing technology which uses the motion of the aircraft or spacecraft carrying the radar to synthesize an antenna aperture much longer than the physical size of the antenna to yield high-resolution

3 imaging capability [4, 6, 7]. This microwave imaging radar is known as an active illumination system, in contrast to passive imaging systems which require the Sun's illumination (like an optical camera) or thermal radiation (like a radiometer). The illumination configuration of SAR systems is usually side-looking where the antenna boresight is perpendicular to the vehicle's direction of flight. The returned signal (echo) originated from the transmitting antenna, mounted on an aircraft or spacecraft, reflects from the objects within the illuminated area (scene), and is received at the same antenna, which is then measured, recorded, and processed to construct the image. This active feature together with the high penetration through clouds and rain enable SAR systems to perform all-weather and round-the-clock geophysical measurements of the Earth over wide surface areas. Since the launch of the Seasat spacecraft in 1978, the technology of obtaining global-scale SAR imagery of the Earth's land and ocean surfaces has been established and refined. Following Seasat, Space Shuttle Imaging Radars (SIR-A in 1981, SIR-B in 1984, and SIR-C/X-SAR in 1994) and operational satellites such as -the European ERS-1 in 1991 and ERS-2 in 1995, the Japanese JERS-1 in 1992, and -the Canadian RADARSAT in 1995 have collected multi-parameter (frequency, polarization, and incidence angle) SAR data. The gathered data show that this technology offers significant and unique features needed to the advancement of scientific research in understanding of environmental change and the impact of various human activities on that. The scientific community has also collected sufficient evidence suggesting that SAR may prove to be extremely valuable in operational applications particularly in monitoring of agriculture/forestry, mapping, resource management, and environmental observations in support of warnings and predictions. Interferometric SAR (INSAR) [8,9] is an evolutionary steep from SAR technol

4 ogy, which can provide topographic information via backscatter phase comparison of successive SAR images obtained from repeat-pass or two-antenna SAR systems. Even though typical imaging radars display only the magnitude of the backscatter, an important attribute of SAR is its coherent imaging capability which allows simultaneous measurement of the radar echo during data acquisition and subsequent processing. INSAR exploits this coherence capability to infer differential range and range change from two or more SAR images of the same scene. In practice, INSAR can generate large-scale high-resolution topographic maps with a horizontal resolution of 30 m and height accuracies of several meters [10]. Over longer time scales of several years or more, high-resolution topographic data can also be used for change detection with accuracy near the millimeter level by comparing elevations of the same scenes at different times. These unprecedented features have made INSAR a valuable technique that has broad applications including 3-D topographic mapping, estimation of vegetation height and ocean wind speed, and studies of catastrophic topographic change due to volcanic eruption, landslides, and major floods. 1.1.3 SAR/INSAR in Forest Remote Sensing SAR/INSAR applications in forest ecology include land-cover classification and estimation of biophysical parameters. For classification applications, various approaches have been investigated, covering the hierarchical knowledge-based approach [11], the dynamic neural network [12], and the interferometrically derived algorithms [13]. For estimation of forest biomass, the multi-frequency (P-, L-, C-band) multi-polarization (hh,vv,vh) JPL AIRSAR system was employed to examine the dynamic range of the SAR backscatter from stands of different ages [14]. The results showed that the P-band SAR is more sensitive for forest parameter retrieval because

5 at low frequencies saturation level of backscatter to forest biomass occurs at high biomass levels. In the studies of SAR dependence on forest parameters, the multiparameter SIR-C/X-SAR data were employed, for the first time, to show that forest structures must be taken into account prior to developing retrieval algorithms [15]. In recent years some experimental and theoretical studies have been carried out to demonstrate the potential of INSARs in retrieving forest parameters. For example in [16] and [17] experimental data using ERS-1 SAR repeat-pass are employed to show the applications of SAR interferometry for classification of forest types and retrieval of tree heights. 1.1.4 Existing Scattering Models for Forest Canopies Interpretation of SAR images is not as straightforward as for optical images because of the complex interaction of the electromagnetic wave with the scatterers. Hence, parallel to the experimental studies, development of theoretical scattering models for vegetation is required in order to interpret the collected SAR data. A tree structure can be used to describe the classification of these models. First, all the scattering models may be divided into two types: phenomenological and physical models. The phenomenological models are based on an intuitive understanding of the relative importance of the individual backscattering components. In this type of scattering model, the total backscattered intensity is calculated by adding up the significant components believed to be important [18, 19]. On the other hand, the physical models are based on the physical interaction of an electromagnetic wave with the forest. This type of model may again be divided into two categories by considering the forest as a continuous medium or a medium with uniform discrete scatterers. For continuous medium modeling, the canopy is treated as a continuous

6 medium with a fluctuating permittivity function [20]. For a discrete model, the canopy is treated as a collection of randomly distributed scatterers with prescribed size and orientation distribution functions for various tree components. In this type of scattering model, the first step is to obtain the general bistatic scattering matrices for various tree components such as trunks, branches, and leaves [21-24]. The solutions of the discrete models may be divided into the distorted Born approximation (DBA) [25,26] and radiative transfer (RT) theory [27-29]. DBA, based on wave theory, starts with Maxwell equations, takes into account the extinction properties of the medium, and calculates the average backscattering coefficients. In general, DBA is complicated and approximations are required for the numerical calculation. On the other hand, RT starts with the radiative transfer equations that govern the transport of radar wave energy through the scattering medium. Its development is quite heuristic, and lacks mathematical rigor. RT theory is the most widely used approach in the radar remote sensing community today. There are three major deficiencies in the incoherent discrete models. First, the incoherent models cannot provide the absolute phase of the backscattered fields, therefore they are inappropriate for INSAR applications. Secondly, the assumption of modeling a forest canopy as a two-layer medium containing uniform distributed scatterers is oversimplified. Actually, the tree structures in the real world have apparently inhomogeneous size distributions in the vertical direction because of the natural branch tapering. Thirdly, the coherent effect caused by the interferences among the relative position of the scatterers and the multipath scattering components is not taken into account. This effect may be an important factor in investigating tree canopies in low-frequency radar remote sensing. Modeling vegetation using coherent approaches has attained significant proini

7 nence over the past five years. For short vegetation, the effect of the soybean plant structure on radar backscatter using a two-scale branching model was first examined in [30]. Similarly, a coherent model for cultural grass canopies, where the dimensions of the vegetation components such as grass blades and stems are comparable to the medium height dimension, was also developed in [31]. In these models the plant structure is considered from a statistical point of view and therefore only the second moments of the scattered fields are provided, that is, the absolute phase information is lost. Further investigations have explored the coherent scattering from a 3-D tree structure. In [32] the radar backscatter was simulated for various deciduous tree types using fractal theory [33] for the tree structure. In a more recent paper [34] Lindenmayer Systems (L-systems) [35], useful tools for implementation of fractal patterns or structures, were employed to develop simple 3-D tree structures of the order of few wavelengths to examine the importance of coherent and multiple scattering. The tree structures in [34] might be unrealistic, since the emphasis was focused on the comparison of various scattering models performed on the same tree structures. In other words, the tree structures can be arbitrarily virtual. A straightforward approach in constructing the tree structure was carried out in [36] where an accurate description of particle positions was characterized for a red pine tree using surveying tools. In addition to the time-consuming and labor-intensive work, this approach excluded the statistical properties of the stand because the average backscatter-ing coefficient was obtained by rotating the same surveyed tree. In this model the tree structure is divided into cylinders whose backscattered fields are added coherently via the distorted Born approximation. Theoretical models have also been developed to establish relationships between

8 both the interferogram phase and the correlation coefficient and the physical parameters of the vegetation and the underlying soil surface [37,38]. Although these models give qualitative explanations for the measured data and provide a basic understanding of the problem, they are not accurate enough for most practical applications because of the oversimplified assumptions in the description of vegetation structure. For example the shape, size, number density, and orientation distributions of vegetation in forest stands are nonuniform along the vertical direction. The nonuniform distributions of the physical parameters of the vegetation structures give rise to inhomogeneous scattering and extinction which significantly affects the correlation coefficient and the location of the vegetation scattering phase center [39]. 1.2 Goals The purpose of this dissertation is to establish a robust scattering model for forest canopies that resolves the aforementioned deficiencies of current models and simulates the polarimetric and interferometric radar backscattering from forest canopies. More specifically, the goals to be achieved in this study include: 1. High fidelity of tree structures. Realistic tree architectures are systematically generated, including both coniferous and deciduous trees. The dependence of tree structures on radar backscatter is investigated. 2. Fully coherent approach. The coherent interference due to the relative position of the tree structures and multipath scattering components are considered. Backscattered fields with magnitudes and phases are predicted, and their statistics are analyzed. 3. Utility for multi-parameter SAR/INSAR systems. The model is capable of

9 simulating polarimetric and interferometric radar responses under various frequencies, polarizations, incidence angles, and baselines. 4. Validation with measurements. The model must be verified using SAR/INSAR measured data with reliable ground truth data. 1.3 Dissertation Overview The above goals have been achieved in this dissertation and the developed approaches and results are depicted in the following chapters. In chapter 2, fractal theory and Lindenmayer systems are employed to generate tree structures, including both coniferous and deciduous types, which are extensively used as the simulated targets in this study. Botanical modeling and 3-D visualization are incorporated to develop realistic tree structures. Electromagnetic modeling of trees and their constituents are addressed in chapter 3 and 4. Once a fractal tree is generated, it is considered as a cluster of scatterers representing the tree components (trunk, branch, and leaves). In chapter 3, the bistatic scattering from individual tree components is investigated. Because of its important role in the total backscatter, scattering from the tree trunk is especially addressed taking into account the effect due to the stratified, corrugated, and tilted features. In chapter 4, a coherent scattering and propagation model for the entire forest canopy is constructed. The backscattered field is obtained from the coherent addition of the individual scattering from each scatterer illuminated by a mean field. Foldy's approximation is invoked to model the coherent wave propagation within the forest canopy where the mean field at a given point accounts for the accumulated attenuation and phase change caused by vegetation particles. A sensitivity analysis is conducted to characterize the dependency on the tree structures, moisture content,

10 and the ground tilt angle. The validation of the model is made by comparing the simulation results with SIR-C data. In chapter 5, the model is applied to the numerical simulation of INSAR response for characterizing the scattering phase center statistics. A Ak radar equivalence algorithm is introduced first to simulate the INSAR response. The concept of scattering phase center is discussed, and its relation with the multipath scattering components is demonstrated. Again, a sensitivity analysis and verification, using TOPSAR data, are also included. In chapter 6, the empirical model and inverse algorithms are derived. The empirical model is constructed based on the simulation results of the Monte Carlo coherent model conducted on a certain stand. The inverse model is implemented for retrieval of forest parameters through genetic algorithms. Lastly, the contributions and recommendation for the future work are described in chapter 7.

CHAPTER II FRACTAL MODEL FOR GENERATION OF TREE STRUCTURES Tree structures are complex and random in nature. A mathematical description of these structures seems to require a large number of independent parameters. However, it has been shown that geometrical features of most botanical structures can be described by only a few parameters using fractal theory [33]. In this study, the fractal concept is employed in the development of the tree generation model which creates realistic tree structures according to the desired forest stand and characterizes the components of the generated trees with specified size, location, and orientation. The geometric data of the generated tree structure are then used in the later electromagnetic modeling of this study. The emphasis is put on the practical aspects of the computerized simulation as well as the fidelity of the tree structure which is a crucial factor when applied to remote sensing applications. 2.1 Fractal Theory and Lindenmayer Systems The mathematical concept of fractals was originated by Mandelbrot [33] in the early seventies. Currently fractal theory is the most popular mathematical model used for relating natural structures to abstract geometries. Mandelbrot defined a 11

12 fractal as a set where Hausdorff-Besicovitch dimension strictly exceeded the topological dimension. In other words, the notion of fractal is defined only in the limit. However, in order to apply the fractal concept to practical problems, a finite curve is usually considered as an approximation of an infinite fractal so long as the significant properties of both are closely related. A distinctive feature of a finite fractal is the self-similarity which is kept through the derivation process. To implement fractal theory, Lindenmayer systems (L-systems) have been wellknown tools for the construction of fractal patterns or structures where the selfsimilarity is preserved through a so-called rewriting process [35]. L-systems were originally proposed by A. Lindenmayer [40], who applied it to the development of lower forms of plant life, such as red algae. L-systems, also called developmental systems, have since been applied in many fields, including formal language theory and biomathematics. The features of L-systems include the structural grammar rules and recursive processes which are very suitable to be implemented by modern computers. In this study, the notation and concept of L-systems are extensively utilized. In L-systems, a tree structure G is specified by three components: (1) a set of edge labels V, (2) an initiator w, called axiom, with labels from V, and (3) a set of tree growth productions P. In compact notation this tree growing process is symbolized by G =< V, w, P >. Given a tree structure G, a tree T2 is directly derived from a tree T1 (T1 => T2) if T2 is obtained from T1 by simultaneously replacing each edge in T1 by its successor according to the production set P. A tree T is generated by G in a derivation of length n if there exists a sequence of trees To, T1,..., Tn such that To = w, and To =* T1 =0... =0 Tn = T. These notations and processing rules are illustrated in the following example.

13 Figure 2.1: An example of the growing process of a fractal tree. Figure 2.1 shows an example of a simple two-dimensional fractal tree of length 4, where the self-similarity can be easily observed through each successive process. In this example, the grammars of growing process are given by G = G(V, w, P) w=X w = X pl: X -- FF{-X}F{+ + X}F{+X}{-X} p2: F -FF where the meanings of labels and operators are described as follows X a seed, no drawing F move one step forward, drawing a line at current angle + increment the current angle (turn to the right) - decrement the current angle (turn to the left) { save the current position and angle on the stack

14 } restore the current position and angle from the stack 2.2 Modeling of Tree Structures In order to simulate realistic tree structures, fundamental botanical properties must be incorporated into the tree generation model. Tree structures in nature grow randomly, but they follow some specific rules according to their various kinds. For example, most deciduous trees have forking branches while coniferous trees have an apparent central trunk growing all the way up to the tip. In this section, several botanical features that have been considered in the development of the tree generation model in this study are described, including three-dimensional (3-D) branching, branch tapering, tree type development, and leaf attachment. 3-D Branching At a node where an originating branch splits into two or more subsequent branches, the branching angle for a subsequent branch is defined by the angle between the originating and subsequent branches. For 3-D branching, the branching angle is specified by two parameters: the tilt angle 0 and the rotation angle 0, as shown in Figure 2.2. Note that in Figure 2.1 only the tilt angle 0 is assigned to the branching angle, therefore the generated tree structure is two-dimensional. In developing algorithms for generating tree structures, the parameters of branching angle are usually the decisive factors to determine the resultant tree structure, which may affect the radar backscatter significantly, as will be shown later in this dissertation. The branching angles are dependent on the tree species and the environment to some extent, therefore the in-situ measurements or observation are sometimes required, especially in the applications of microwave remote sensing, to generate faithful tree structures. Branch Tapering The cross section and the length of subsequent branches decrease

15 I c Figure 2.2: 3-D branching and branch tapering at a splitting node. as the branching processes progress. When a branch splits, as shown in Figure 2.2, a common practice for determining the relationship between the radius of the originating branch (ra) and the radius of the subsequent branches (rb and rc) is the application of the conservation law of the cross sectional area. Take an example for a branch splitting into two branches, it can be mathematically described by: ra = r rc. (2.1) The relationship between the radius of the new branches (i.e. rb and rc) is specific to the tree type, and should be prescribed according to the in-situ measurements. Similarly, the lengths of the originating and subsequent branches

16 (la and lb, respectively) can be related by a tapering factor g, expressed by lb = la/g. (2.2) Again, g is a user-defined parameter which may have different values depending on the current position in a tree structure. This parameter may affect the tree height and crown coverage. Tree Type Development Based on the architectural characteristics, tree structures can be categorized into three primary classes: columnar, decurrent, and excurrent, which may be represented by coconut, maple, and pine trees respectively [41]. In the majority of deciduous trees, the lateral branches grow as fast as, or faster than, the terminal shoot, giving rise to the deliquescent growth habit where the central stem eventually disappears from repeated forking to form a large spreading crown. This branching pattern is termed decurrent. On the other hand, most coniferous species belong to the excurrent class. In the branching pattern of excurrent trees, the main stem outgrows the lateral branches giving rise to cone-shaped crowns and a clearly defined bole. In this study, a library of typical tree structures is constructed that can easily be fine-tuned to simulate the desired tree stands. Figures 2.3 and 2.4 show, respectively, the fractal trees of deciduous and coniferous types. In L-systerns, the developing concept for generating deciduous trees is quite different from that for coniferous trees. For example, only 4 rewriting processes are used to generate the deciduous tree shown in Figure 2.3, while 12 rewriting processes are required to perform the coniferous tree shown in Figure 2.4. Leaf Attachment Efficient algorithms for leaf attachment may greatly reduce the modeling tasks. Once the fractal tree is created, the numerical data specifying

17 '." 9.......... Figure 2.3: A fractal tree structure of deciduous type. the location, size, and orientation of each tree component (trunk, branch or leaf) are stored. The number of leaves for a coniferous (or prosperous deciduous) tree are so large, that positioning the needles on branches and specifying their location and orientation would require significant computation time and

18 Figure 2.4: A fractal tree structure of coniferous type. storage memory. To solve this problem, a concept of cluster leaflets is used. For example, a close examination of red pine trees revealed that the position of needles on branches follows precise mathematical patterns. Leaf buds of red pine are equally spaced on three concentric spirals drawn on a branch. Regu

19 End Cluster | Concentric Spirals Stem Cluster Figure 2.5: Leaf attachment schemes for the red pine: the end cluster, the stem cluster, and the leaflet on the leaf bud along three concentric spirals. larly, one spiral period contains 7 leaf-buds of which each grows a two-needle leaflet, as shown in Figure 2.5. Also shown in this figure are the end cluster which is found at the tip of each branch and the stem cluster which is composed of a thin stem/branch and surrounding needles. It can be observed that the needle density and orientation of an end cluster are quite different from those of a stem cluster. Nevertheless, the end clusters at individual branches are very much identical, and can be treated as a single scatterer whose orientation is defined by that of the connecting branch. Figure 2.6 shows a fractal-generated red pine with attached needles, which consists of 792 branch segments, 391 end clusters, and 747 stem clusters. The tree generation algorithms developed in

20 Figure 2.6: A red pine with attached needles. this study allow the user to specify the number of leaves/needles per stem as well as a local orientation angle distribution for the attached leaves.

21 2.3 Computer Implementation The computer work in the development of the tree structure consists of three main components: the encoding, decoding, and visualization. In L-systems, the encoding is accomplished by iterating the labels with prescribed productions and length. A long label string, like DNA in biology, is obtained at the end of the processes, holding embedded information about the tree structure. Then this long label string is decoded (or translated) into a tree structure through a so called turtle graph interpreter [35]. Numerical calculation is performed in this stage to quantify the geometries of the entire tree structure. A fast Fortran code has been developed for this end which allows for user-defined parameters. Some special features of the tree structure must be treated carefully using user-defined parameters. For example, to model the upward property observed at the end of each branch of red pine, as shown in Figure 2.6, the branch angle with the position dependence must be incorporated into the branch angle parameter. It is important to visualize the structures of trees during the development process. Once the fractal tree is created, the tree data file usually contains a large number of tree components and it is difficult to examine the accuracy by manual inspection of the numerical data. Visual inspection of the tree image is a better way at this point. In addition, real-time visualization of the tree structure during the developing stage can also assist the user in learning the sensitivity of the fractal parameters to the tree structure. In this study, a visualization program is developed using the PostScript language where real-time display and printout can be easily performed without any extra software. This program is capable of projecting a 3-D fractal tree structure into a

22 2-D image with the functions of arbitrary scaling and perspective view. The red pine shown in Figure 2.6 is viewed at 200 measured from the horizontal plane. Moreover, Figure 2.7 shows top- and side-view images of a coniferous tree of which broad leaves (modeled as circular disks) are attached to the tree skeleton shown in Figure 2.3. 2.4 Application to Remote Sensing of Vegetation Not only qualitative but quantitative fidelity of the tree structure is required when the application of remote sensing is concerned. Although there are many computer graphics available for generating tree-like structures, most are not appropriate for the purpose of remote sensing modeling. At microwave frequencies, the radar return and its statistics strongly depend on tree structures, which necessitates the application of a realistic model for generating accurate tree structures, rather than just 'tree-like' structures. Therefore, the final and most important step in the fractal model is to incorporate the information about the tree structure and its statistics obtained from the in-situ ground truth measurements. Figures 2.8 and 2.9 compare the photographs of the actual forest stands with the simulated trees produced by the fractal model developed in this study. These simulated structures are generated according to the in-situ measurements collected from two test stands denoted by Stand 22 (red pine) and Stand 31 (red maple) near Raco, Michigan. The parameters for L-systems for generating the maple tree in Figure 2.9 is given below. Note that some user-defined symbols other than those used by L-systems [35] are rendered to account for the sophisticated features of the tree structures. For example, the branch tapering symbols (),[], and {} are used to denote small, medium, and large branch tapering respectively.

23 * Fractal Coding Parameters: Length (number of iteration) n: 4 Axiom/initiator w: FFF!(+A){!FF(+A){!FF(++A){!F(+B)!(++B){F!(+A)!(++A) [!F[++B]![++B]![+B]]}}}} Production pl: A: — ff(+A)!f(+A){!(++A){!f[+A]!f[+B]}}} Production p2: B: -- f(+A)[!f(++B)[!f[++B]!f[+B][-B]]] Production p3: F: -FF Production p4: f: - ff where F and f are respectively vertical and horizontal forward steps. * Geometrical parameters include tree diameter at breast height (DBH), branching angle (0b) and rotation angle (0b), trunk tilt angle (0t), leaf orientation angle (Os), number of leaflet (Nleafiet), leaf radius (al), leaf thickness (ti), stem radius (as), and stem length (1s). mean(DBH), std(DBH), mean(F), std(F), mean(f), std(f) (cm) 14. 3. 9. 1. 7.5 1. mean(0b), std(Ob), mean(4b),std(bb) (deg) 22. 5. 137.5 10. mean(Ot), std(Ot), mean(Ol), std(0i) (deg) 0. 3. 5. 10. Nleaflet, al, tl, as, Is (cm)

24 5 4.0 0.02 0.1 8. where mean(-) and std(.) refers to the statistical mean and the standard deviation of the parameters respectively. * Display postscript parameters: Xstart, Ystart, Zstart, Oview, qkview, Scale 300 300 80 90 0 25

25 f Figure 2.7: The top and side views of a deciduous tree with attached leaves.

26 (b) Fractal Pine (a) Stand 22 Figure 2.8: The photograph of a red pine stand (Stand 22), and the simulated tree structure using the fractal model.

27 ( (a) Stand 31 (b) Fractal Maple Figure 2.9: The photograph of a red maple stand (Stand 31), and the simulated tree structure using the fractal model.

CHAPTER III ELECTROMAGNETIC SCATTERING MODEL FOR A TREE TRUNK ABOVE A TILTED GROUND PLANE 3.1 Introduction Once a tree is generated, it is treated as a cluster of scatterers composed of cylinders (trunks and branches) and disks (leaves), each with specified size, position, and orientation. It is very complicated to simulate the radar backscatter from this highly chaotic random media. The entire problem, however, can be decomposed into two sub-problems -the scattering from individual tree components and the total backscattering from the whole forest canopies taking into account the propagation extinction. Of all the tree components, the tree trunk is the largest and therefore most important. In this chapter the general bistatic scattering from a tree trunk is investigated. The total backscatter of the tree canopy will be discussed in the next chapter. In the conventional scattering models for forest canopies, a tree trunk is simply modeled by a vertical, homogeneous, finite-length dielectric cylinder. The scattering solutions for a finite-length cylinder, reported in the literature, are either based on the eigen-function expansion solution for an infinite cylinder [42-44], or low fre 28

29 quency approximation where all dimensions of the cylinder are small compared to the wavelength [23,45]. When the cylinder radius is large compared to the wavelength the eigen-function solution becomes, numerically, inefficient due to the poor rate of convergence of the series involved in the solution. This is the case in the microwave region where the radius of tree trunks in a forest stand can be significantly larger than the wavelength. An inefficient solution for the calculation of scattering properties of a canopy's constituent particles makes the canopy model numerically intractable because the scattering solution for individual particles must be evaluated many times to account for the particle variability in size and orientation. Moreover, in modeling a tree trunk as a dielectric cylinder, the stratified property and the barklayer roughness have been overlooked. In the real world, the dielectric profile of a tree trunk cross section is inhomogeneous with multiple concentric rings. Furthermore, for many trees the bark layer is rough and can be represented by longitudinal grooves on the surface of a dielectric cylinder. The effect of the bark layer on the radar cross section (RCS) of a tree trunk was demonstrated recently by representing the bark layer with a corrugated dielectric layer [46]. Using a hybrid scattering model based on the method of moments and physical optics it was shown that the RCS of a tree trunk is significantly reduced when the effect of the bark layer is taken into account. However this model is not numerically efficient enough to be used in conjunction with the scattering model for a forest canopy. In this study, a realistic and efficient scattering model for a tree trunk above a ground plane is developed. The coordinate systems are first defined in both local and global terms. The local coordinate system is for the calculation of the scattering matrix of the individual tree components, while the global coordinate system is for the whole environment configuration including the slope of the ground plane, the

30 position and orientation of the individual tree components, and the radar perspective. In the scattering model for a tree trunk, the effect of the radial inhomogeneity as well as the rough bark layer are taken into account. For finite-length cylinders having radii comparable to the wavelength, the eigen-function expansion in conjunction with the field equivalence principle is used. At high frequencies where the radius of curvature is large compared to the wavelength, the physical optics (PO) approximation is employed, based on the fact that the dielectric constants of tree trunks are highly lossy. The bark layer is represented by a periodic corrugated layer and equivalently replaced by an anisotropic layer as suggested in [47]. The ground plane is considered to be a half-space dielectric medium having a smooth surface. Both the cylinder and the ground plane are allowed to have arbitrary orientation described by the global coordinates. The direction and polarization of the incident wave are also arbitrary. At the end of this chapter, thorough numerical simulations are conducted to demonstrate the region of validity of the PO approximation, the effects of a bark layer and the ground tilt angle on the backscattering RCS. 3.2 Global Coordinate System and Notations In this study, the problem of scattering from an object above a ground plane in the most general configuration is considered. A global coordinate system (X, Y, Z) is constructed to describe the two-parameter directional unit vector u, given by u(0, &) = Xsin0cos+Ysin0sin + +ZcosO. (3.1) This unit vector can represent the incidence direction ki(Oi,oi), the scattering direction k,(O8,s ), the orientation of the object ic(c, c), or the unit normal vector ng{(0g, fg) of the ground plane, as shown in Figure 3.1. In this coordinate system, the horizontal and vertical polarization of the incident and the scattered waves are

31 A hi A k A A 5. S / " ^: Figure 3.1: Global coordinate system and notations. defined by hp = Z x kp/\Z x kp\ (3.2) A ' l vp = hp x kp (3.3) where subscript p can be i or s representing the incident and scattered waves respectively. In this study, the forward scattering alignment convention [45] is assumed. The components of the scattered field Es and the incident field E' in the global coordinate system can be related to each other by the scattering matrix S. defined

32 by eikor E = S - E', (3.4) r or in the matrix form (EV eikor ShS Svh) Ev (35) Eh VSh Shh VEh Since a tree is modeled as a class of discrete scatterers above a ground plane, the most fundamental problem is the calculation of the coherent scattered field frorm a single object above the ground plane. First, we model the ground plane by a halfspace dielectric medium having a smooth interface and an arbitrary tilt angle in the global coordinate system. In this case the effect of the ground plane on scattering can simply be taken into account by including the mirror image contributions. Neglecting the multiple scattering between the scatterer and its mirror image, the total bistatic scattering matrix S, as shown in Figure 3.2, consists of four main terms: 1) the direct scattering St, 2) the ground-target scattering Stg, 3) the target-ground scattering Sgt, and 4) the ground-target-ground scattering Sgtg. Coherent addition of these first order scattering terms gives rise to the total scattering matrix S, S = St + Sgt + Stg + Sgtg, (3.6) where S = S ~(kki), (3.7) Sgt = eiTsr(ksg, kgS). S~(s, k), (3.8) St, = eTiS(ksI, ki) r(ki,,n, ki), (3.9) Sgtg = ei(7+rs)r(ks, ng, kgs) S~(kgs kgi) (kg ni, ki), (3.10) and kgki= ki-2ng(ng. ki), (3.11)

33 A / A nl A g / SOt /g / A ks / A A k' k Ag ks kgs g Figure 3.2: The first-order scattering components for scattering from an object above a ground plane. kgs = ks-2ng(n ks), (3.12) Ti = -2koh(ng ki), (3.13) Tr = 2koh(hg ks). (3.1.4) In the above expressions, h is the height measured from the geometric center of the target to its projection on the ground plane surface. Here, the optical lengths Ti and Ts account for the extra path lengths of the image excitation and the image scattered waves respectively. S0(ks, ki) is the scattering matrix of the isolated target in free space, and is called the intrinsic scattering matrix. r(kr, nig, ki) is the reflection mnatrix of the ground plane accounting for the reflection and polarization transformation. In order to provide a physical insight for the four main terms, subscripts t and

34 g, representing the target and the ground plane respectively, indicate the sequence of scattering from right to left. The unit vectors (ks, ks, kgi and kgs) indicating the propagating directions are also expressed in the arguments of S~ and F in the same manner. In the following sections, a general polarimetric reflection matrix F for a tilted ground plane with arbitrary slope is first obtained, and then the bistatic expressions for the scattering matrix S~ of a stratified finite cylinder in free space based on the eigen-function expansion and the PO approximation are derived. 3.3 Polarimetric Reflection from A Tilted Ground Plane In most scattering models for forest stands, the effect of the surface topography is ignored. In a recent paper [48], the importance of the ground tilt angle in the radar backscatter was emphasized through the SAR measured data. However, it lacks a rigorous and general formulation to account for the effect of an arbitrarily tilted ground plane in a polarimetric radar system. In this section an expression for the reflection matrix of a tilted ground plane is derived. Consider a tilted ground plane with a unit normal ig(0g(, 4g) that is illuminated by a plane wave propagating in the direction ki(0, i). The direction of the reflected wave is determined by kr= ki-2ng(hgki), (3.15) which is normal to the polarization of the reflected field Er having Er and Er as its vertical and horizontal components in the global coordinate system. Defining the reflection coefficient matrix F by A = r.16 E = r(kr I ki) -El (3.16)

35 the objective is to express the elements of F in terms of the Fresnel reflection coefficients of the ground plane. In the local coordinate system of the ground plane, tIhe vertical and horizontal polarizations of a wave are defined by hp = x9g x kp/llng X kpl (3.17) AI A"I p = hp x kp (3.18) where the subscript p can be i or r. By representing both the incident and reflected field vectors in the local coordinate system (Vp, hp, kp) and noting that =- (F i):) (3.:19) V Eh, V O 0 J V Eh, J the elements of the reflection coefficient matrix can be obtained from Fpq = (pr' vi)F,(I (vqi) + (pr' h)F)rh(h' I) (3.20) where p and q can be v or h, and Fv, and Fh, are, respectively, the vertical and horizontal Fresnel reflection coefficients of the ground plane. The inner products in the above expression in terms of the global coordinate parameters are given by A'.A A ' n- Z-(k3Z)(ng.kj) j. 3 - 3lgxkjJhxkjA^ ^ A I ng. hj Vj hj -hj v. = lg 3 3xgkjI where j can be i or r. As will be shown in the later numerical simulation, a tilted ground plane has two significant effects on the backscatter: (1) the ground-trunk component in the backscattering direction may decrease drastically off the peak of strong specular backscattering due to the slightly tilted angle, and (2) a significant cross-polarized component can be generated by a slanted ground plane and/or a cylinder object which is not oriented in the principal plane.

36 3.4 Semi-Exact Solution for Electromagnetic Scattering from A Finite Stratified Dielectric Cylinder Scattered fields of an infinite stratified cylinder can be obtained by the standard eigen-function expansion method [42]. However, for finite-length cylinders, no exact solution exists. In the microwave region where the length of a tree trunk is much larger than the wavelength and the dielectric constant has a significant imaginary part, the effect of the longitudinal traveling waves on a finite cylinder can be ignored. Therefore, the internal fields of a finite cylinder may be approximated by those of an infinite cylinder having the same radial characteristics. In [22], the bistatic scattering from a finite stratified cylinder is addressed by invoking the field equivalence principle. However, the formulae given in [22] are quite brief and symbolic. In this section, a detailed alternative expression ending up with an explicit form of the scattering matrix is given in terms of the global parameters. Consider an infinite stratified cylinder with M layers oriented along z axis. The radius of the m-th layer is denoted by pm and its dielectric constant by cm, as shown in Figure 3.3. At an oblique incidence, the z-component of the incident fields can be expanded as [45] 00 Ez(p,, z) = S ezJn(kopsin ),, ((,z), (3.21) n=-oo oo ZoHz(p,, z) = S hzJn(kopsin /)Fn(q,z), (3.22) n=-oo where e, = Ei, z (3.23) h, = ZoH', (3.24) Fn(, z) = (-i)neineikozcosp3 (3.25)

37 \ P \ \ / x \, A P\\ \, ki Figure 3.3: Geometry and composition of a finite stratified cylinder (top-view and side-view). Here, Fn is the mode function associated with the z-dependent phase, and Jn is the Bessel function of the first kind. Since the cylinder is infinite, the z-dependent phase in F, is kept in the scattered fields, which can be written as 00 Es(p,,z) = 3 AnH(p) (kopsinp) Fn(q, z), (3.26) n — n=-oo ZoHS(p,,z) = E BBnH(1)(kopsin/3)Fn(b )z), (3.27) where Hn() is the Hankel function of the first kind accounting for the outgoing wave ( eiwt is assumed and suppressed). The coefficients An and Bn are functions of the

38 polarization of the incident wave and can be expressed by An Cez + Chz (3.28) Bn =Cn- + Chz. (3.29) Note that the cross-polarized coefficients Cn and Cnh are included and it can be shown that C = -Ch if the cylinder is homogeneous [42]. Expressing the q$-component of the internal fields, Ed and Ho, in terms of Ez and Hz and matching these tangential fields along each layer interface [44], a closed-form solution of the coefficients can be obtained from = ( 1 - a 1), (3.30) 2 0308 - a4a7 C 2a3 - 4 ), (3.31) 2 a3a8 - Of4C7 C'h = 1(c~c8-(6cX7) (3.32) 2 a3a8 - a4a7 h = 1 a3C6 - c4a5 (333) 2 a308 - a4a7 where a1 = d11 + d12, a2 = d13+ d14, a3 = d21 + d22, 04 = d23+ d24, ca5 = d31+d32, a6 = d33+d34 c7 = d41 +d42, a8 = d43+d44, and dij is the element of matrix D4x4 which can be calculated from L()1 L(pi)... L-1(pm)1 L L(pM). (3.34) Here, counts for the distribution of the matrix accounts for the distribution of the internal fields of the n-th mode confined within the m-th layer (pm+l < p < pm), and is given by

39 HHl)(k(p) H2)(km p) 0 0 0 0 Hnl)(k p) Hn 2(kpp) -nkocosH)(krp) -nkocos H (k) -ikoHk (km)2p (km)2p kk -ikoEm H(l) kmp) ko m (2)' km -nko Cos H() -nkoc Cos H(2)(k mp) kmp P kk( p -kkp c n (P (kp)2p n P where k = -ko -os2, (3.35) and cm is the dielectric constant of the m-th layer, as shown in Figure 3.3. By invoking the field equivalence principle, the scattered field outside the cylinder region can be attributed to fictitious electric and magnetic currents which are related to the total fields over the surface of the cylinder [21]. To find the scattered fields for a cylinder of finite length b, we can assume the currents on a finite cylinder are the same as those of the infinite cylinder. After some mathematic manipulation similar to [21], the scattering amplitude of the finite stratified cylinder is given by S = osin V [(Peez + x x (Qhh+ h Qkhez)k x i], (3.36) where Pe = J (-1) e '(xo)J(yo)- -in(xo )(yo ) (3.37) B (X0) 00 n=-oo +ce [H(i't(xo)Jn(yo) - sin H)(x)(Y)] (3.38) in cos (3 xo sin /3 k. - xo yoB cos fiC +Qh1- = (1)e0 5 X )JY(xo) (339) Qh_ = (-l) e ( B J( (o)J (yo) (3.40) n= —o- I

40 B Pe =+ ( —[H) n' C H ')(o) J- (o)S sin H(1)(xo)J(y) (3.413) in cos/3 xo sin3 k' z + (1 )Jn(yo) [Jn(xo) CH() () (3.44) xo yoB cos3p J with V = (k-k), (3.L7) o = kopl sin /, (3.48) YO = kopl B, (3.49) n=-oo L J iB cos ( 3.s i 0), = tan-l ( L:). (3.51) Note that the above expression of scattering amplitude S is in the global coordinates. By definition, the elements of the scattering matrix are calculated from Spq = y + S(Ci), (3.52) where p and q can be v or h. Note that ez - i, hz = i., for v-pol incidence; e = hi, hz = -vi i, for h-pol incidence; and. ( x k x ) = -S. sz, (3.53)

41 hs (ks x ks x ) = -hs i, (3.54) Vs (ks x z) = -hs z, (3.55) hs. (ks xz) = Vs-z. (3.56) Therefore, the intrinsic scattering matrix in the global coordinate can be expressed in the matrix form: kopsin V V[s hs Pe Pe zvl z h 2sin VA A A J A A A h' -s Qh Qh h i -'i 3.5 Physical Optics Approximation The semi-exact solution described in the previous section becomes inefficient at high frequencies where the radius of the cylinder is large compared to the wavelength and fails when the cross section of the cylinder is not circular. These deficiencies can be removed at high frequencies by employing the PO approximation. This approximation is valid when the radius of curvature of the cylinder is large compared to the wavelength and the permittivity of the cylinder has a relatively large imaginary part so that the effect of the glory rays and the creeping waves can be ignored. As before, the cylinder is replaced by fictitious electric and magnetic currents; however, the currents in this case are approximated by those of the local tangential plane which are proportional to the sum of the incident and reflected fields. To simplify the integration of the currents over the lit surface, the stationary phase (SP) approximation may be used. This approximation is valid so long as the stationary point falls over the lit region. For convenience, a local coordinate (n, t. I) is established at the SP point. The local tangential directions are defined by t = fi x ki/n x ki( (3.58)

42 I = n x t (3.59) where n is a unit vector normal to the cylinder surface at the SP point. For the general case of an anisotropic medium (the bark layer may exhibit anisotropic properties) a dyadic reflection coefficient R is introduced to relate the polarization coupling between the incident and reflected waves, i.e. ET = R.ET. (3.60) Combining the incident and reflected fields, the total fields E(= E'r + E') and H(= HT + H') on the surface of the cylinder can be obtained from El 1 + RVV Rvh ( E 31 I= (3.61) Et \ Rhv 1 + Rhh E and ( Hi (1 - Rhh Rhv ( H) (3 \Ht \ Rh 1 - Rv V Ht Applying the stationary phase approximation, it can be shown that the Hertz vector potentials are given by = Zo eikor IIe o e QJ (3.63) iYo eikor II = ~ QM. (3.64) where J and M are the fictitious currents evaluated from the total fields E and H at the SP point ('/ = ), and Q j jz e-ikoaBc. (3.65) 47Q -b/2 J-T/2 ib sin V -ikoBa a {F kBa v7r F [koBa l)]} = r e ~ 2 +F 2 -+ (2-+) 2wr V 2B [ 2 2 J2 2

43 with B = {[(k - k) ]2 )[(k -ki) Y]2}1/ - tan-l((k i) ' and F(-) is the Fresnel Integral. This approximation is valid provided koaB >> 1 and q is away from the shadow boundary. Using a similar procedure as in the previous section, the scattering matrix elements are found to be V = Q[(i )ZoJI, + (i. V)ZoJty + (I h,)Mv + (t. hs)Mt] (3.66a) Sh = Q[(i. s)ZoJlh + (t. v)ZoJth + (i hs)Mh + (. hs)Mth] (3.66b) S = Q[(l h,)ZoJ1, + (t. h,)ZoJ - (i vs)Ml - (t. s)Mt] (3.66c) Sh = Q[(- h,)ZoJlh + (t. hs)ZoJth - (. Vs)Mh - (t. s)Mth] (3.66d) where Jpq and Mpq are the currents along the p direction induced by a q polarized incident wave (p can be t or I and q can be v or h ). The inner products of the vectors in the above expressions can easily be calculated in terms of the global coordinates. The above results fail in the case of forward scattering for which B = 0. However, in directions close to the forward direction, an alternative approximation for the scattered field is possible [21], given by -2iab A sin V sin W eikoi E = (ki.)xi) rE, (3.67) where W = koa(k y') and V is given in (3.47). 3.6 Modeling of a Corrugated Bark Layer For some tree species, the bark layer is corrugated with grooves along the longitudinal direction. In this study, the bark is simply modeled as a periodic corrugated

44 d e3 L (a) e3 % I t ~ e, (b) Figure 3.4: A corrugated layer (a) and its equivalent anisotropic layer (b). layer with period L and width d as shown in Figure 3.4. It is shown in [47] that,when L < Ao/2 (single Bragg mode), the corrugated layer can be equivalently replaced by an anisotropic layer (see Figure 3.4) with the same thickness whose permittivity tensor is given by 611 0 0 e = 0 E22 0. (3.68) 0 0 0 33 The entries of the tensor in terms of the permittivity, period, and width of the corrugated layer, when L < 0.2Ao, are approximated by e6l = (1-dL d(3.6;9) Cr(1 - d/L) + d/~

45 622 = 33 = 1 + (r - l)d/L. (3.70) Assuming that the radius of the cylinder is much larger than the wavelength, the permittivity of the bark layer can be represented by c(4, z', n) where Ie = 1 and Czz = Enn = C22. To employ the PO approximation, a coordinate transformation from the local (j, z', n) to (t, 1, n) at the SP point is needed. The resultant permittivity tensor in coordinate (t, 1, n) is ezz sin2 qz + eCo cos2 z (e+ - e,,) sin 4, cos qz 0 E = (e - Ezz) sin z cos 4z ezz cos2 z + 6e, sin2 z 0 (3.71) 0 0 Enn where 1- - z( 0z = cos- () (3.72) /1-(n k. k)2 The reflected field from a stratified anisotropic dielectric half space is computed using the method described in [49]. 3.7 Simulation Results In this section a number of numerical examples for the scattering from a finite cylinder above a ground plane are presented. In all the considered examples the normalized RCS, defined by 47tISpq\12 apq kab2 (3.73) is displayed for a two-layer cylinder with height b, exterior radius a, and interior radius a2. The permittivity of the exterior and interior layers are chosen to be

46 ie = 4 + ii and 62 = 10 + i5 respectively. The cylinder is positioned on a tilted ground with permittivity cg = 10 + i5. First, the validity region of the PO approximation in the backscatter direction is examined. Figure 3.5 compares a,, and ahh using the PO and semi-exact solutions. It is found that the PO solution agrees well with the semi-exact solution when koa > 10. For small values of koa the resonance behavior of the backscatter is shown by the semi-exact solution. Figures 3.6, 3.7, and 3.8 show the monostatic and bistatic polarimetric scattering patterns which are simulated for a two-layer cylinder with and without corrugation. The thickness of the corrugated layer and its filling factor are respectively chosen to be t = 0.1A0 and d/L = 0.7 (see Figure 3.4). Figure 3.6 shows the backscattering pattern as a function of incidence angle. At small angles of incidence, the PO approximation differs slightly from the semi-exact solution because the radial component of the propagation constant (kp = k0 sin 0) is small in this region and the condition kpa > 10 is not satisfied. The vv-polarized backscattering RCS has two minima corresponding to the two Brewster angles, one occurring on the surface of the cylinder (0 ~ 25~) and the other occurring on the ground plane (0 - 75~). The ripples on the curves are due to the coherent interference between the components St and Sgtg, of which the oscillation amplitudes become significant at incidence angles close to 90~ (grazing angles) where the normal incidence of the cylinder occurs. Moreover, the oscillation rate of the ripples is proportional to the res s ppa thight of the cylinder. Figure 3.6 also shows the effect of the bark layer on the backscattering RCS. Depending on the incidence angle, the RCS of the cylinder may be reduced as high as 10 dB. The reduction in the RCS is a function of the cylinder length and the corrugation parameters. Basically the corrugated layer behaves as an impedance transformer

47 "0 -5 --10 - -15 C -20/ vv (Semi-Exact) - - - - Chh (Semi-Exact) -25 (P.O. Approx.) -30 L.. j 0 5 10 15 20 25 30 35 40 k.a Figure 3.5: Comparison of the PO approximation with the semi-exact solution. The ratio of the interior radius a2 and the exterior radius a1(= a) is kept constant (a2/al = 0.9). Other parameters are: b = 20Ao, 61 = 4+il, e2 = g = 10 + i5,i - 120~, i = 180~, Os =- 60~,,s = 0~, Og 0 = 0~.

48 I I I I I I I I " -10 -20 x! \ O \ hh (Semi-Exact) o / - - Ohh (P.O.) -\ / o (Semi-Exact) -30 - - vv (P.O.) -Corrugated, ohh C orrugated, Ovv -40 I I I I I i 0 10 20 30 40 50 60 70 80 90 Incidence Angle n-0i (Degrees) Figure 3.6: The normalized backscattering RCS as a function of the incidence angle o = 7r- Oi for a two-layer cylinder with and without the corrugation. Other parameters are: b = 20A0, al = 20, a2 = 1.8A,t = O.1A0,d/L - 0.7,? = 180~, q = 0~, 0g = 0c = 0~. between the air and the vegetation material. Figure 3.7 shows the bistatic scattering pattern as a function of elevation angle when 0i = 120~, = 180~ and the observation point is moving in the XZ-plane. Figure 3.8 shows the bistatic scattering pattern as a function of azimuth angle ($s) with 0i = 120~, 5 = 180~, and 0, = 60~. The discontinuities found on the PO solution near the forward directions are because the scattering formula in (3.67) is used instead of that in (3.66).

49 (a) "0 Iz) 10 0 -10 -20 -30 -40 -50 -60 '-90 -60 -30 0 30 60 Scattering Elevation Angle Os (Degrees) 90 Figure 3.7: The normalized bistatic a,, (a) and ahh (b) as a function of the scattering elevation angle 0, in X-Z plane (O, > 0 when;, = 0~; 0, < 0 when Os =- 180~). The backscattered and the specular directions are shown at 0s == 60~ and 0, = -60~ respectively. Other parameters are: i = 120~, qi = 180~, 09 = 0c = 0~ b = 1OAo, al = 2Ao, a2 = 1.8o.

50 (b) 20... 1 0 Semi-Exact 10 P.O. I I - Corrugated 0 - 1.iiI- i iiiiI -20 - II a II 11 'I. II II/ i -40 '!1\ v i 1' I I!t -50. ',' -90 -60 -30 0 30 60 90 Scattering Elevation Angle 9s (Degrees)

51 (a) 10 5 Semi-Exact 0 P.O. - Corrugated -5 --O -l - -. 15 -20 -30 0 30 60 90 120 150 180 Scattering Azimuth Angle Os (Degrees) Figure 3.8: The normalized bistatic orr, (a) and Crhh (b) as a function of the scattering azimuth angle qs. The backscattered and forward-scattered directions are shown at 4s = 0~ and 4s = 180~ respectively. Other parameters are: Oi = 120~, i= 180~, 0 = 60~, 09 = 0 = 0~, b = 20Ao, a, = 4Ao, a2 = 3.6Ao.

52 (b) 20 Semi-Exact P.O. 10 - -- - Corrugated 0 -20 ' \ -30 0 30 60 90 120 150 Scattering Azimuth Angle <s (Degrees) 180

53 Figures 3.9, 3.10, and 3.11 show the effects of the tilted ground plane on the backscattering RCS. All the parameters in Figure 3.9 are the same as those given in Figure 3.6 except for the tilt angle of the ground at Og = 20~ and qSg = 90~. Comparing Figure 3.6 with Figure 3.9, it can be seen that a significant cross-polarized backscattered signal is generated due to the slope of the ground plane. Figure 3.10 shows the variation of backscattering RCS as a function of the ground azimuth angle Og where 0, = 135~, qi = 180~ and Og = 20~. One can observe that the peak of the backscattering RCS occurs at qg = 70~. Figure 3.11 shows the backscattering RCS as a function of the ground elevation angle 0g where 0Q = 135~, b = 180~ and qg = 0~ and 1800. The regions in the positive and the negative 0g represent the ascending and descending sides of a mountain respectively. In this case no cross-polarized signal is generated because the cylinder is in the principal plane (X-Z plane). Note that two maxima take place at 0g = 0~ and 0g = -22.5~. The first maximum corresponds to the dihedral-like ground -trunk interaction. The second maximum corresponds to a reflection from the ground plane which illuminates the cylinder at normal incidence. The backscatter from the cylinder bounces off the ground plane and returns toward the radar (see Figure 3.12). This strong backscatter component can be observed where ki, nig and Z, are in the same plane and satisfy the condition: 7' 0i = 20g + 2. (3.74) 3.8 Conclusions An efficient and realistic electromagnetic scattering model for a tree trunk above a ground plane is presented in this chapter. The trunk is modeled as a finite-length stratified dielectric cylinder with a corrugated bark layer arbitrarily positioned and

54 oriented above a ground plane. The ground is considered to be a smooth homogeneous dielectric with an arbitrary slope. An asymptotic solution based on the PO approximation for high frequencies is derived. This solution provides a fast algorithm with excellent accuracy when the radii of tree trunks are large compared to the wave — length. The effect of the bark layer is also taken into account by simply replacing the bark layer with an anisotropic layer. It is shown that the corrugated layer acts as an impedance transformer which may significantly decrease the backscattering RCS. The RCS reduction depends on the corrugation parameters. It is also shown that for a tilted ground plane a significant cross-polarized backscattered signal is generated while the co-polarized backscattered signal is reduced.

55 o -10 hh v3 -- -- o.....Ovh,, - -20 U -30 L 10 20 30 40 50 60 70 80 90 Incidence Elevation Angle n-0i (Degrees) Figure 3.9: The normalized backscattering RCS as a function of the incidence angle 7r - Oi for a two-layer cylinder above a tilted ground plane. Other parameters are: 0g = 20~', = 90~, c0~, = 180~, = 0~, 0, = r - -, b = 22 = So 20A0, al -- 2,~o, a2 -- 1.8)to.

56 HH-Pol -10 -....... VV-Pol.. VH-pol -20 -30 X - Z -40 -50 / ' "..:: ) i: il. '.I ' ' x ' 0 30 60 90 120 150 180 Scattering Azimuth Angle tg (degrees) Figure 3.10: The normalized backscattering RCS as a function of the ground azimuth angle qbg. Other parameters are: 0g = 20~, c =- 0~, i = 135~, i = 180~, 0s = 45~, Os = 0~, b = 20Ao, a, -= 2Ao a2 = 1.8Ao.

57 ---- VV-Pol -10 -HH-Pol R:: A I aI -20 I/) I I^:; a Th,,,,,,,a as a..o o '~...... 1 *; i i';; I ^,, ~ q<" ', ',,';. ',i' ' " \ a are. ", =1,, a 18)~ g... ', a, ' i= a aI, a, a a.a,,' I) II J _ i | -40 - 1:1: 1 -50 I -50 -40 -30 -20 -10 0 10 20 30 40 50 Scattering Azimuth Angle 0g (degrees) Figure 3.11: The normalized backscattering RCS as a function of the ground elevation angle 0g (Og > 0 when 4g = 0~; Og < 0 when 4g = 1800). Other parameters are: I0 = 35, = 180,O = 450, = 00,b = 10Aoa = 2Ao,a2 = 1.8Ao.

58 A Figure 3.12 The geometry of the trscatteing configurations for a cylinder over a tilted ground where a strong backscatter can be observed.

CHAPTER IV COHERENT SCATTERING AND PROPAGATION MODEL FOR CONIFEROUS AND DECIDUOUS CANOPIES 4.1 Introduction Among the existing scattering models for forest canopies [18,19,26-29], radiative transfer (RT) theory [50] is the most widely used theory for characterization of radar response to a forest canopy. When the medium consists of sparse scatterers that are small and far separated compared to the field correlation length within the random medium, RT theory can accurately predict the second moments of the radar backscatter statistics. RT theory is an incoherent approach based on the energy conservation and therefore cannot provide the phase information of the backscattered field. However, the phase quantity is of importance required for investigating the response of a forest to an interferometric SAR. The other shortcoming of RT theory is its inability to account for the coherent effects caused by the relative position of the discrete scatterers within a tree structure as well as the multi-path scattering components. Recent investigations on radar backscatter from forest canopies have shown that both backscatter and attenuation are significantly influenced by tree architectures [51]. Therefore, it is crucial to develop a coherent scattering model 59

60 that simulates the phase statistics of the backscattered field and takes into account; the coherent effect due to the tree structures. In this study, an accurate and comprehensive coherent scattering model for forest canopies is developed, which can preserve the coherent interference among the tree structure but also provide information about the absolute phase of the backscattered field from forest canopies. The proposed model is composed of three major components: (1) accurate generation of tree structures based on only a few fractal parameters, (2) evaluation of scattered fields from a single tree structure, and (3) Monte Carlo simulation of the entire canopies taking into account the coherent wave propagation and extinction. In the tree structure modeling, stochastic L-systems [35] based on fractal theory are employed to construct a realistic tree structure according to a user-defined input parameter set, as described in Chapter II. As will be shown in this chapter, the spatial and angular distribution of branches strongly influences the behavior of radar backscatter, indicating that the tree structure should be taken into account in modeling the radar backscatter and estimating vegetation biomass from the forest canopy. In the scattering model, individual tree components located above a tilted ground plane are illuminated by the mean field, and the scattered fields are computed and then added coherently. The branches and tree trunks are modeled by stratified dielectric cylinders, as described in Chapter III, and leaves are modeled by dielectric disks and needles of arbitrary cross sections. The mean field at a given point within a tree structure, which accounts for the accumulated phase change and attenuation due to the scattering and absorption losses of vegetation particles, is calculated using Foldy's approximation [50]. Finally, a Monte Carlo simulation is performed on a large number of fractal generated trees to characterize the statistics of the backscattered fields. The developed coherent model is able to account for

61 the effect of an inhomogeneous scattering and extinction profile of a forest canopy. The accuracy of the model is compared with the measured backscattering coefficient acquired by SIR-C SAR over a forest test stand near Raco, Michigan. In what follows, the first order coherent scattering approximation for calculating the backscatter from a tree structure is described first. Next the coherent wave propagation and extinction within the tree canopies of both coniferous and deciduous types are discussed. Then the model results based on the Monte Carlo simulation are examined against L- and C-band backscatter measurements over a well characterized forest stand. At last, a sensitivity study is carried out to demonstrate the dependency of individual forest parameters on the multi-parameter (frequency, polarization, and incidence angle) radar response. 4.2 Electromagnetic Scattering from A Tree Structure In this study, the first order coherent scattering approximation is employed to calculate the polarimetric scattered field of a tree structure. Once a tree is created, it is treated as a cluster of scatterers composed of cylinders (trunks and branches) and disks (leaves) with specified position, orientation, and geometric shape and size. It is assumed that the entire tree is illuminated by a plane wave, whose direction of propagation is denoted by a unit vector ki and is given by E'(r) = E e.kokir (4.1) The scattered field in the far zone is to be calculated for individual trees, as shown in Figure 4.1. Since the uncertainty in the relative position of trees with respect to each other is usually of the order of many wavelengths, the total scattered power from a forest canopy can simply be determined by the incoherent addition of scattered power from individual trees.

62 Es Ei / A ng IL Figure 4.1: Scattering from a cluster of scatterers above a tilted ground plane. To the first order of approximation, the scattering from a tree is approximated by the superposition of the scattered field from each scatterer within the tree structure. Hence, neglecting the effect of multiple scattering among the scatterers, the total scattered field from a single tree can be evaluated from ikr N Es= e einSn~ Et, (4.2) r n=l where N is the total number of the scatterers within a tree structure, Sn is the scattering matrix of the n-th scatterer above a dielectric ground plane and On is a phase compensation term accounting for the shift of the phase reference from the local coordinate system of the n-th scatterer to the global coordinate phase reference, which is assumed at the joint of the tree and the ground plane, as shown in Figure 4.1. Denoting the position of the n-th scatterer in the global coordinate system by

63 4t S n t n sgt 'S t D A / / Ags 9} Figure 4.2: Four scattering components from an object above a tilted ground plane.. rn1, 4 is given by n = (k - ks) ~ rn, (4.3) where ks is the unit vector representing the propagation direction of the scattered field. In order to compute the local scattering matrix Sn, let us consider a single scatterer above a ground plane. Neglecting the multiple scattering between the scatterer and its mirror image, each scatterer mainly contributes four scattering components, denoted by: (1) direct scattering denoted by St, (2) ground-scatterer scattering denoted by S'9, (3) scatterer-ground scattering denoted by Sat, and (4) groundscatterer-ground scattering denoted by Sn, as shown in Figure 4.2. The scattering matrix Sn can be written in terms of its components as Snt + Sgt + st + S949 (4.4)

64 where St = S~(k:i), (4.51) S9t eisR(kskg S~ (kgs,ki), (4.6) S - eiS~(ks,kgj) R(kgi,, k ), (4.7) S -g = ei('+)R(ks, kgs) S~(kgs, kg) R(kgi ki), (4.8) and kg = ki- 2ng(nAg ki), (4.9) kgs ks,-2g(ng ks), (4.10) Ti = -2ko(rn ns)( ki)), (4.11) Ts = 2ko(rn nig)(ng' ks). (4.12) In the above expressions, S~ is the bistatic scattering matrix of the n-th scatterer in free space. The direction of incidence and scattering are denoted by unit vectors in the arguments of So. In the above expressions, hg is the unit vector normal to the tilted ground plane. The phase terms r; and Ts account for the extra path lengths of the image excitation and the image scattered waves respectively. R is the reflection matrix of the ground plane whose elements are derived in terms of the Fresnel reflection coefficients and the polarization transformation due to the ground tilt angle. The explicit expressions of the reflection matrix of a tilted ground plane (R) with an arbitrary slope and the expressions for the bistatic scattering matrices (SO) of large scatterers like trunks and primary branches are given in Chapter III, where the semi-exact solution as well as the physical optics approximation are derived for the calculation of scattering from a stratified dielectric cylinder above a tilted ground plane. The formulae for the scattering matrices of small scatterers like needles

65 and broad leaves are constructed based on the expressions given in [23,24]. 4.3 Coherent Wave Propagation within Tree Canopies The analysis in the previous section is not quite complete since in the calculation of scattering from the n-th scatterer the other scatterers are assumed to be transparent. The second or higher-order analysis, which takes into account the multiple scattering among the tree structures, is fairly complicated and is beyond the scope of this study. However, the effect of attenuation and phase change of the coherent wave propagating in the random media can be readily modeled by calculating the mean field within the random medium. Consider a coherent radar wave propagating in a statistically uniform random medium. Based on Foldy's approximation [50], the variation of the mean field E with respect to the distance s along the direction k is generally governed by dE -= iK E, (4.13) ds where ko + Mv Mvhh K =, (4.14) Mhv ko + Mhh and Mpq = k < S~q(kk) >' (4.15) Here ko is the wave number of the free space; no is the volume density of the scatterer; and < S (k, k) > is the ensemble average of the forward scattering matrix, (p and q can be v or h). Using the standard eigen-analysis, the above differential equation (4.13) can easily be solved and the solution is given by E(s) = eiksT(s k). EO, (4.16)

66 where E~ is the field at s = 0 and T is the transmissivity matrix accounting for the extinction due to scattering and absorption. In most natural tree structures, azimuthal symmetry can be assumed where M h = Mhv = 0 and the transmissivity matrix is reduced to Ci'Mvvs 0 T =. (4.17) 0 eiMhhs Note that the transmissivity matrix in (4.16) excludes the phase terms due to free space path lengths, and merely accounts for the perturbation in propagation caused by the vegetation. To include the effect of wave extinction in the scattering model, consider a situation when the entire tree structure is embedded in an effective medium with an effective propagation constant given by (4.15). Under the aforementioned approximations the expressions for the components of the n-th scattering matrix in the backscatter direction should be modified as S = - T S~(-ki, ki) T, (4.18) Sgt = ei"Tt R T. S~(-i, ki) T, (4.19) Stg - einTt S~( —ki] kr)' T. R Tt, (4.20) Sgtg = ei2nTt. R. T. S(-krkr) T. R Tt (4.21) with kr = ki - 2ng(s~g ki), (4.22) rn = 2ko(rn ng)(ing kr), (4.23) where Tf,Te, and Tt are the transmissivity matrices, respectively, for the direct, reflected, and total traveling path as shown in Figure 4.3. In the above equations the reciprocal property of wave propagation is employed, i.e., T(s, k) = T(s, -k),

67 Figure 4.3: Extinction of a coherent wave in random media. which results in the expected reciprocal scattering relation, Sgt = (Stg)-t. Here the superscript (.)-t denotes the operation of matrix transposition followed by negation of the cross-polarized elements in order to be consistent with the the forward scattering alignment convention [45]. Due to lack of knowledge, the aforementioned existing scattering models assumed the distribution of vegetation particle size and type to be uniform. However, in the real world, this distributions of vegetation particle type and size is non-uniform along the vertical extent of most forest stands. Based on the fractal model developed in this study where the exact description of particle distributions is readily available, propagation and scattering of the mean field within the forest medium can be characterized quite accurately. To calculate the transmissivity matrices Tn. T, and Tt, different propagation models are developed for deciduous and coniferous forest canopies.

68 4.3.1 Deciduous Canopies For deciduous forest stands it can be assumed that the canopy is continuous and homogeneous in horizontal direction but inhomogeneous in vertical direction. To account for the vertical inhomogeneity, consider an M-layered random media representing the tree canopy which is illuminated by a plane wave, as shown in Figure 4.3. Each layer, with thickness dm(m - 1,2,..., M), is assumed to be parallel to the ground plane, and its boundaries are diffuse such that no reflection or refraction can take place. In the backscatter case, only the incident directions ~cki and the reflected directions ~kr are of interest. Therefore for each layer (say the m-th layer), the layered transmissivity matrix is computed by "i;'~mLm 0 T' (Lm), (4.24) 1riMi Lm 0 CMhh m where Lm = drn/(tg kr) is the slant thickness, and M = kdmt, E/ S~(ki,r, >, (4.25) n is the averaged effective propagation constant of the m-th layer. In the above expression Dt is the tree density and Nm is the number of particles of a single tree in the m-th layer.. Here < * > is implemented by a sufficient number of independent realizations of fractal trees through the Monte Carlo simulation. Suppose the n-th scatterer is within the m-th layer, the final expressions for the transmissivity matrices can be written as Tt = T (L1)T'(L2)...T'(LM), (4.26) T = Tm(Ln)Tm+l (Lm+l)... TM(LM), (4.27) Tn = TM(Lmn)Ti-(Lm_)....Ti(Li), (4.28) M 1Um

69 where Ln and LJ are the slant distances from the n-th scatter to the top and bottom of the m-th layer boundary along the ki and kr directions respectively, as shown in Figure 4.3. These distances are given by Lmn = (Hm - rn)/(ng. kr) (4.29) Lrn = (rn, g -— H)/(ng A kr) (4.30) where m Hm = dk (4.31) k=l represents the height of the upper interface of the m-th layer. Figure 4.4 shows the calculated extinction profile of the simulated maple stand (Stand 31) where the entire tree canopies are divided into 11 uneven layers and the extinction coefficient is defined as the imaginary part of the effective propagation constant Mpq for each layer. It is noted that the wave attenuation at C-band is much greater than L-band, and the extinction coefficient for vertical polarization is slightly greater than that for horizontal polarization at both frequencies. This extinction coefficient profile is shaped according to the tree architecture and composition, which plays an important role in radar backscatter parameters including the position of the scattering phase center [37]. In this example, the entire tree canopy is divided into eleven layers, and the extinction coefficient is calculated as described in the previous section. It should be pointed out that the number of layers can be determined by imposing a step discontinuity threshold. Basically the algorithm starts with a moderate number of layers, calculates the extinction coefficient for each layer, and examines the step discontinuity. If the discontinuity between any two layers is larger than the prescribed threshold, these layers are divided into finer layers.

70 16 14 Lhh 12 - Cvv — s — Chh o I 0..... ' ----.. — Lhh 0.00 0.10 0.20 0.30 Extinction Coefficient( Np/m ) Figure 4.4: The extinction profile calculated for Stand 31.

71 4.3.2 Coniferous Canopies For coniferous canopies, the mean field within the canopy depends on both the vertical and the horizontal location of observation point. In general the canopy can be divided into three main regions: (I) conical crown, (II) continuous layer of overlapped branches, and (III) trunk layer, as shown in Figure 4.5(a). The trunk and the overlapped branches layers are treated as homogeneous layers similar to the deciduous type. However, the conical crown layer is discontinuous. Based on numerical simulation it was found that within the conical region the effective propagation constant is almost constant. This is due to the fact that at the top region there are more needles and thinner/shorter branches while at the lower region there are less needles and thicker/longer branches. Therefore, the transmissivity is merely determined by the path length traveled by the mean field within the conical region. This path lengths can be obtained easily from the intersection of the propagation line, X - Xn y — Yn Z- Zn432) XX __ - _ ZZ, (4.32) kx y z and the cone surface, (x - xt)2 + (y - yt)2 = (zt- z) tan(a/2), Zb < z < t, (4.33) where (xni Yn, zn) is the position of the n-th scatterer, (Xt, Yt, Zt) is the position of the cone tip, and a is the cone angle, as shown in Figure 4.5(a). Note that for incidence angle 0 > a/2 the adjacent trees may shadow each other. To calculate the path length across the neighboring trees, the coordinated cone tip (Xt, Yt, Zt) of the neighboring tree should be used in (4.33). The positions of the neighboring trees are chosen randomly according to the tree density. Figure 4.5(b) shows the position-dependent transmissivity of the incident wave at three different locations

72 A k/ /e, r I,.. i (a) 0 P?..., c> C/O E -5 -10 -15 0 10 20 30 40 50 60 70 80 90 0 (deg) (b) Figure 4.5: The position dependence of the transmissivity for coniferous trees and the shadow effect caused by neighboring trees.

73 shown in Figure 4.5(a) as a function of incidence angle. The path lengths for the other scattering components can be obtained using a similar ray-tracing approach. 4.4 Monte Carlo Simulation In our model the radar backscatter is expressed in terms of the scattering matrix. However, for distributed targets, the radar backscattering coefficients and phase difference statistics are usually the quantities of interest. These quantities can be derived from the second moments of the backscattered field components [52]. The statistics of the scattered field are approximated from a Monte Carlo simulation where a large number of tree structures are generated using stochastic L-systems and then the scattering matrix of all generated trees are computed. Computation of the scattering matrices is accomplished in the following manner. First the canopy height is discretized into M layers and the extinction coefficients of individual layers and the integrated transmissivity matrices are computed as outlined previously. Then these quantities are used in (4.18)-(4.21) for calculating the scattering matrix of individual trees. The computation involved in the calculation of the scattering matrices of individual leaves for many trees is too extensive to be carried out even with the fastest available computers. To solve this problem, the 47r solid angle covering the entire vector space representing the orientation direction of a leaf is discretized into a finite number, and a look-up table for scattering matrices of a leaf oriented along all the discrete directions is generated for the three principal directions: backscattering (So(-k,kki) and So(-kr,kr)), forward scattering (S~(ki,ki) and So(k, kr)), and bistatic scattering (So(-k,, ki) and S(-ki, kr)). The number of discrete orientation directions is determined from the ratio of a typical leaf dimension to the wavelength

74 (a/A). According; to this scheme the number of the discrete points should increase with increasing a/A. A similar scheme may be used for branches: however, we found that this may unnecessarily increase the CPU time due to a large variability in diameter and length of the branches. In order to calculate the desired backscatter statistics, the differential covariance matrix of the backscattered field must be evaluated. As described earlier, the backscattered fields of adjacent trees in a forest are uncorrelated at microwave frequencies and above. Therefore the backscattered power from individual trees can be added and the covariance matrix elements are proportional to the tree density Dt and are given by W~qs Dt< SPqSst >, (4.34) where p, q, s, t E {v, h}. According to this definition for the differential covariance matrix, the backscattering coefficient can be obtained from P= 4Wp7r (4.35) pq ppqpq 4.5 Model Validation In this section, the accuracy and validity of the developed model is examined using a set of measured data acquired by the Shuttle Imaging Radar-C/X-Band Synthetic Aperture Radar (SIR-C/X-SAR). The collected ground truth and the radar parameters, such as frequency and incidence angle, are used as model input. In this section we also present some examples to demonstrate the sensitivity of radar backscatter to some important forest parameters.

75 4.5.1 SIR-C/X-SAR The SIR-C/XI-SAR radar system [53] was flown aboard the shuttle Endeavor in the spring (SRL-1) and fall (SRL-2) of 1994. This mission was the first of its kind where a beam-steerable, multi-frequency, and multi-polarization space-borne synthetic aperture radar was deployed. The SIR-C/X-SAR system operated at L(1.25 GHz), C- (5.3 GHz), and X-band (9.6 GIIz). The L- and C-band SARs were configured to collect polarimetric data whereas the X-band SAR was a single channel radar and collected the backscatter data at vv polarization. The look angle of the system was varied from 15~ to 60~. In this study, the polarimetric SIR-C data (Land C-band) during the SRL-2 is compared with the results predicted by the model developed in this study. 4.5.2 Ground Truth Raco Supersite, located in the eastern part of Michigan's Upper Peninsula, was designated by NASA as a calibration and ecological Supersite for the external calibration of SIR-C and has been a test site for our radar remote sensing activities since 1991 [54,55j. Great efforts have been devoted towards characterizing ground inventories, and the site has been imaged by ERS-1, JERS-1, SIR-C/X-SAR, and JPL AIRSAR. The main research objective at this site has been relating the measured SAR backscatter data to the forest ecological/biophysical parameters, which are essential input parameters for the ecological models used for the study of land and atmosphere processes. The Raco Supersite contains most boreal forest species and many of the temperate species. The SIR-C/X-SAR overflight occurred in the fall, a time of some seasonal change where trees begin to dry and the deciduous leaves begin to undergo their

76 Tree Density: 1700/Hectare Tree Height: 16.8 m Trunk Diameter (DBH): 14 cm Leaf Density: 382 #/m3 Leaf Area: 50 cm2/# Leaf Thickness: 0.2 mm Leaf Moisture (mg): 0.51 Wood Moisture (mg): 0.60 Soil Moisture (ma): 0.18 Table 4.1: Ground Truth of Stand 31 fall color change. During the SIR-C overflight (October 1994), the leaves were still predominantly green. Color change happened towards the end of the mission. In this study, a deciduous forest stand, denoted in the existing reports as Stand 31 [54], is selected as a test stand. This stand consists of a large number of red maple as well as a few sugar maple, uniformly covering an area about 300 m by 300 m on flat terrain. The ground truth of this stand has been collected since 1991, and a summary of its pertinent parameters is reported in Table 4.1. The vegetation and soil dielectric constants during the SIR-C overflights are reported in Table 4.2, and are derived from the measured moisture values using the empirical models described in [56,57]. 4.5.3 Simulation Results The first step in obtaining the model prediction is to generate fractal trees faithful to the real tree structure of the desired forest stand. There are two phases for determining the input parameters for the tree generating code. The first phase is

77 L-band C-band Leaf 17.9 + i6.0 14.7 + i4.7 Wood 32.1 + ilO.O 27.7 + i8.4 Soil 9.7 + il.6 9.4 + il.5 Table 4.2: Dielectric properties of Stand 31 to characterize the coarse parameters such as the branching nature of the trees, the growth factors, and the finite fractal order. In the second phase, some of the fine input parameters, such as the branch tilt angle and its distribution, are slightly tuned in order to minimize the difference between the simulated and measured backscattering coefficients a~. In general accomplishing the second phase is much more difficult than the first phase because there is no apparent rule for adjusting the parameters. To establish a set of rules of thumb for fine-tuning the parameters of tree structures, we performed a sensitivity analysis. The gradient of the desired radar backscatter parameters with respect to the desired tree structure parameters was determined and used for determining the fine tuning procedure. In this procedure, we allowed the fine tree parameters to be adjusted to within 10% of the measured ground truth parameter to account for the uncertainty in the ground truth measurements. The model (including tree generation and coherent scattering) is assumed to be verified if the simulation results can simultaneously match with the polarimetric SIR-C data for both frequencies and different incidence angles. In performing Monte Carlo simulations, one should be careful of the convergence properties of the simulation. In all simulation results reported in this study, conver —

78 L-band C-band 0 0 co3 m -5 hh -5 w 0 5 - - ---------------- -- --------- V -20 -20 '- ' vh -25 -25; -30 -30 0 50 100 0 50 100 Realization Number Realization Number (a) (b) Figure 4.6: Convergence behavior of the Monte Carlo simulation for Stand 31 at L-band(a) and C-band (b) at incidence angle QO = 30~. gence was achieved to within ~0.5 dB of the estimated mean values for less than 100 tree realizations. Figure 4.6 shows, respectively, the convergence behavior of the backscattering coefficients at L- and C-band for forest Stand 31. Figure 4.7 shows the comparison between the model prediction and the measured backscattering coefficients for three consecutive SIR-C overflights as a function of incidence angle at L- and C-band respectively. It is shown that an excellent agreement is achieved for all incidence angles and polarizations except for the C-band cross-polarized backscattering coefficient. The lack of accuracy for this polarization can be attributed to the effect of multiple scattering between branches or branches and leaves in the canopy crown. It shows that the measured C-band data is consis

79 tently higher than the simulated results by 1.3 dB which can be attributed to the overestimation of radiometric calibration constants. The computation time for each incidence angle point is about 35 minutes at L-band and 65 minutes at C-band on a Sun Sparc-20 workstation. As mentioned previously, the total backscatter is comprised of different scattering components. Simulation results show that in all cases except for L-band hh polarization the backscattering coefficients are dominated by the direct backscatter component (or). The hh-polarized backscattering coefficient chh at L-band, depending on the incidence angle, is mostly dominated by the direct backscatter or the ground bounce term (4ot). The double ground bounce component (~tg) is negligible for all cases because of the low transmissivity for this canopy (see Figure 4.11 for LAI=12). Figure 4.8 shows the scattering components of crh as a function of incidence angle. The analysis for characterizing the contribution of each scattering component is essential in determining the position of the scattering phase center of the forest. It is also important to examine the effect of the inhomogeneity of the extinction profile (see Figure 4.4). Figure 4.9 compares the co-polarized backscattering coefficients of Stand 31 where the forest is both modeled by a 2-layer medium and by an 11-layer medium. In the 2-layer model the tree canopy is composed of a trunk layer extending from 0-5 m and a crown layer which extends from 5-17 m. It can be ob — served that the 2-layer model overestimates the backscatter at lower incidence angles and underestimates at higher incidence angles. The discrepancy in this example is as high as 2.5 dB, and can be even higher for stands with higher leaf density. It is also found that the discrepancy increases with increasing frequency. For example, the discrepancy at L-band is less than 0.3 dB. It should be mentioned that the CPU

80 ~i -5 5F -1 - --- - ~ --— A --- —-- Lvv (SIR-C) - Lhh (SIR-C) 0 A Thh (SIR-C) em2 -15 - -- Lvv (Model) -- Lvh (Model) ------—. Lh (Model) -20...... Lh (Model) = -25 10 20 30 40 50 60 70 Incidence Angle (degrees) (a) m O *. I. i.... -5 -o A = —. o. Cv (SIR-C) t -10 0 --- —------ - E CVh (SIR-C) U 0 A Chh (SIR-C) eo -15 - - - - -. - ---- - - - Cv (Model) - -- -- Cvh (Model) gC -20 20 ----..Chh (Model) g -25 - -25.-L.-.... I.... a....I...... 10 20 30 40 50 60 70 Incidence Angle (degrees) (b) Figure 4.7: Comparison between the model predictions (lines) and SIR-C data (symbols) at (a) L-band and (b) C-band.

81 0.... E -10 -' --- o o 's -5 -- -------------------------------------------------- 0 N -15 ------- t.... ' gt, -20 - 0___ Ct 20 total > F..................... 'I\...1 P -25 10 20 30 40 50 60 70 Incidence Angle (degrees) Figure 4.8: Contribution of different scattering components (or~ and cr~) to the overall backscattering coefficient (Ctotal) as a function of incidence angle.

82 time for calculation of the backscattering coefficient for a 2-layer and an 11-layer forest is almost the same, because the mean field profile of the canopy is calculated before the Monte Carlo simulation is carried out. The statistical behavior of the backscatter can also be obtained from the present model. Through the Monte Carlo simulations the desired histograms can be constructed by recording the backscatter results for each realization. Figure 4.10 shows the estimated probability density function (pdf) of backscattering coefficients in dB3 at incidence angle 43.6~. The pdf can provide additional information about the distributed target if the backscatter statistics are non-Gaussian. For instance, although the mean values of hh and oa at L-band are nearly identical, their pdfs are somewhat different. The transmissivity is another quantity used to characterize a stand. Based on the extinction profile of the forest canopies, the transmissivity can be computed by integrating the attenuation of each layer. In Figure 4.11, the one-way transmissivity (from top to bottom) is calculated as a function of leaf area index (LAI), defined as the total leaf area (single side) per unit area of forest. It is shown that the horizontally polarized wave can more easily penetrate the canopies than the vertically polarized wave. This phenomenon results from the fact that the tree trunk and branches are oriented mostly along the vertical direction. To demonstrate the effect of tree structures on the radar backscatter two examples are considered. In the first example denoted as Case 1, we changed the branching angle AO from 22~ ~ 5~ (used in Stand 31) to 15~ ~ 3~ while keeping the other parameters the same. Figure 4.12 compares the orientation distribution of the branches for Stand 31 and Case 1. The pdfs of branch orientation are obtained by counting the number of branches in small increments of orientation angle for all the

83 "Ic O 0 U d)) 0 0.4 j 4-)) M uC cnP -5 -10 ' I ' ' ' I..I ' '. I. ' ' I ' '. I ' r. - - 11-layer (Chh), -- 2-layer (Chh) i -15 1 (. ~, I ~ ~. ~ ~ ~ " ~ ~ ~ ~ ~, ~ I. ~ " ) 20 30 40 50 60 70 80 Incidence Angle 0i (Degrees) (a) -5........................... "0 0 m 0 o o*, C, CZ cm c, -10 h ' ' I ' ' ' I ' ' ' I ' ' ' I ' ' ' I ' ' ' I ' ' ' ll-layer (Cvv) ~~ ---3.. ---- 2-layer (Cvv).. I. ~.. I... I.,.. I..,,,,. I -15 L 10 20 30 40 50 60 70 80 Incidence Angle Oi (Degrees) (b) Figure 4.9: Comparison of the 11-layer and 2-layer extinction models simulated at C-band for hh polarization (a) and vv polarization (b).

84 0. L-band.1II 0. C-band.1II vv-pol )5m LL: 6i 0.0 0.0 -40 -20 0 20 -40 -20 0 20 0. LI: ci0. 0 0:.1 I - vh-pol '5 - 0, -40 -20 0 20 I.1 vh-pol -40 -20 2C n lank hh-pol 0. UL: 6i0. 0.1I hh-pol '5 - 0, 0.c W -40 -20 0 aY0 (dB) 20 4 -20 0 G 0 (dB) 20 Figure 4.10: Histogram of the backscattering coefficients for Stand 31 at Oi = 43.6'.

85 -5......... -........"-... -- A -------- --- -10 -15 -i ] ---;' ' -. Lvv. -20 -- M -3 — C, -20 -- - c -25'.. ' I.. '. I I, -25. 0.0 5.0 0 10.0 15.0 20.0 LAI Figure 4.11: One-way transmissivity as a function of leaf area index (LAI).

86 (1 f25 r- I I I I I I I I 0.02,.,0.015 '0 Q' 0.01 0.005 - (a) rn frf f I F Stand 31 1[ri IrI............................................. 0 10 20 30 40 50 60 70 80 90 Ob 90 Ob Figure 4.12: Histogram of branch orientation resulted from two branching angles (a) AOb = 22~ ~ 5~ (Stand 31); (b)AOb = 15~ ~ 3~ (Case 1). branch segments of a fractal tree, which included about 7500 branch segments in more than 20 classes of diameters and lengths. These two tree structures are shown in Figure 4.13 (a) and (b). In the second example, referred to as Case 2, we changed the tree height and trunk diameter while keeping the dry biomass unchanged (14.4 kg/m2). A tree structure of Case 2 having height 8.6 m and trunk diameter 17cm is shown in Figure 4.13(c). The backscattering coefficients calculated for these tree structures are given in Table 4.3. Simulation results indicate a significant variability in the backscattering coefficients among the three simulated forest stands of different geometrical structures having identical biomass.

87 (a) (b) (c) Figure 4.13: Tree structures generated for (a) Stand 31, (b) Case 1, and (c) Case 2. Tree L-band C-band Structures au (dB) () ( dB) Chh ( dB) <T~ (dB) hh (dB) Stand 31 -8.8 -14.6 -8.2 -9.3 -16.4 -10.1 Case 1 -8.4 -16.1 -7.1 -7.9 -16.2 -8.7 Case 2 -13.2 -19.6 -9.1 -11.4 -21.6 -10.3 Table 4.3: Effect of tree structures on backscattering coefficients, simulated at Oi = 43.6~.

88 The effect of the ground tilt angle on the radar backscatter from the forest canopy is also investigated in this study. Using the same parameters of Stand 31 except for changing the ground tilt angle 0g from 0~ to 10~, the backscattering coefficients are computed as a function of the azimuthal look angle Xi for incidence angle Oi at 25.40. As expected, simulation results show that the radar backscatter is not sensitive to the ground plane at C-band since most of the backscattered field emanates from the crown layer. However, at L-band the radar backscatters, especially aUh and Uvh, show sensitivity to the ground tilt angle (see Figure 4.14(a)). Near the azimuthal look angle /i = 70~, there are noticeable increases in orhh and aCh which can be attributed to a ground-trunk interaction. For a cylinder oriented along the zc direction standing on a tilted ground plane with a unit normal ng, the ground-trunk scattered ray is parallel to the incident ray when the following relationship is satisfied: ki Zc = (ki -ng)(9g zc). (4.36) Assuming Zc is along the vertical direction Z, the above equation can be readily reduced to an explicit expression given by 2 sin2 g cos -. (4.37) oitan 0i sin 20g ( Using this expression a maximum backscatter is expected at Xi = 68~ when 0g = 10~ and 0i = 25.4~. For this particular look angle a significant cross-polarized backscatter is generated. Figure 4.14(c) shows the backscattering cross section of a cylinder with radius a = 7.2cm, length b = 7.2m, and dielectric constant c = 32.1 + il10.0 on a tilted ground plane with 0g = 10~ and cg = 9.7 + il.6 illuminated by a plane wave with incidence angle 0i = 25.4~, as shown in Figure 4.14(b). It is well known that the backscatter from forest canopies as a function of frequency is governed by various scattering mechanisms. For instance, at lower fre

89 1-o 0 Q 4.). 0 -5. ~ ~~ ~ ~ ~.......... I I I I (a) hh - vh I. ---- Wr^ - ' - - - -10 F -15 - ') Ifl L. I.... I - 20 40 60 80 100 Azimuthal Look Angle, {i (Degrees) 120 (b) z A, A A ki Fc I y I -WI -— r - - I - -,l.....0.0j: - - r-. -. x 20 - 10 U) r 0 rA, U1; -10 CZ -20 - 20 40 60 80 100 120 Azimuthal Look Angle, o, (Degrees) Figure 4.14: The effect of the ground tilt angle on radar backscatter: (a) backscattering coefficients for Stand 31 over a tilted ground with 0g = 10~ simulated at L-band with incidence angle Oi = 25.4~, (b) a vertical cylinder on a tilted dielectric plane and the associated coordinates, (c) RCS simulation for (b).

90 5.... I................ i......... o VV real structure (seven-layer) o HH real structure (seven-layer) 0 < VV random dist. (two-layer) * VHH random dist. (two-layer) -10 - 0 — -. -.-15 -15........................................ 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 Frequency (GHz) L/C/X - Band Figure 4.15: Comparison of the polarimetric backscattering coefficients at different frequency bands (L, C, and X). quencies the dominant backscatter term is the trunk-ground bounce, while at high frequencies the backscatter is mainly influenced by the direct scattering term from the crown layer. The difference between the backscatter from a forest stand using a multilayer model and a two-layer model is examined. In the two-layer model the forest is divided into a trunk layer (0-8 m) and a crown layer (8-17 m), and extinction parameters are averaged overte the whole layer. Figure 4.15 shows the results of this study where the backscattering coefficients are computed at different frequencies, L-band (1.25 GHz), C-band (5.3 GHz), and X-band (9,6 GHz). It is shown that at high frequencies the discrepancy between the two models increases.

91 4.6 Conclusions A coherent scattering model for forest canopies based on a Monte Carlo simulation of fractal generated trees is developed in this study. A coherent model offers three major advantages over the existing incoherent scattering models: (1) the model preserves the effect of the aritectural structure of the trees which manifests itself in the extinction and scattering profiles; (2) the model provides complete statistics of the scattered field instead of just its second moments; and (3) the model is capable of simulating the scattering from forest stands on tilted ground planes. In general the coherent scattering model is composed of two main components: (1) a tree structure generating model which is developed based on stochastic L-systems, and (2) a firstorder scattering model which can handle radially stratified cylinders and dielectric disks and needles of arbitrary cross section. The validity and accuracy of the model was demonstrated by comparing the results based on the model simulation with the backscattering coefficients measured by polarimetric L- and C-band SIR-C at three different incidence angles. A very good agreement between the measured quantities and the model predictions is obtained with the exception of aoh for C-band. This disagreement may result from the existence of multiple scattering in the crown layer. A sensitivity study was also carried out to demonstrate the effects of tilted ground planes and tree structures on the radar backscattering coefficients.

CHAPTER V SIMULATION OF INTERFEROMETRIC SAR RESPONSE FOR CHARACTERIZING THE SCATTERING PHASE CENTER STATISTICS OF FOREST CANOPIES 5.1 Introduction As mentioned previously, one of the features of the present model is its ability to predict the absolute phase of the radar backscatter from forest canopies. This feature makes possible the numerical simulation of the interferometric synthetic aperture radar (INSAR) over a forest canopy. In this chapter, the present coherent scattering model is applied to simulate INSAR response to forest canopies in terms of the scattering phase center height and the correlation coefficient. Recent advancements in the field of radar interferometry have opened a new door to the radar remote sensing of vegetation. In addition to the backscattering coefficient, radar interferometers measure the scene of interest with two additional quantities, namely the correlation coefficient and interferogram phase [58]. To interpret these quantities and to characterize their dependency on the physical parameters of the target, a thorough understanding of the coherent interaction of electromagnetic waves with vegetation particles is required. 92

93 The premise of this investigation with regard to retrieving vegetation parameters from INSAR data stems from the fact that the location of the scattering phase center of a target is a strong function of the target structure. For example the scattering phase centers of non-vegetated terrain are located at or slightly below the surface depending upon the wavelength and the dielectric properties of the surface media, whereas for vegetated terrain, these scattering phase centers lie at or above the surface depending upon the wavelength of the SAR and the vegetation attributes. It also must be recognized that the vegetation cover adds noise in other interferometric SAR applications where the geological field mapping or surface change monitoring, instead of vegetation itself, is the primary interest. In these cases it is also important to characterize the effect of vegetation on the topographic information obtained from the interferometric SAR and to develop correcting schemes for those vegetated areas. In recent years some experimental and theoretical studies have been carried out to demonstrate the potential application of INSARs in retrieving forest parameters. For example in [16] and [17] experimental data using ERS-1 SAR repeat-pass are employed to show the applications of SAR interferometry for classification of forest types and retrieval of tree heights. Also theoretical models have been developed to establish relationships between the interferogramrn phase and correlation coefficient to the physical parameters of vegetation and the underlying soil surface [37,38]. Although these models give qualitative explanation for the measured data and provide a basic understanding of the problem, due to the oversimplified assumptions in the description of vegetation structure, they are not accurate enough for most practical applications. For example the shape, size, number density, and orientation distributions of vegetation in forest stands are nonuniform along the vertical direction. The nonuniform distributions of physical parameters of vegetation particles

94 give rise to inhomogeneous scattering and extinction which significantly affects the correlation coefficient and the location of the vegetation scattering phase center. On the other hand, there are a number of electromagnetic scattering models available for vegetation canopies [27-29], but most are of little use with regard to INSAR applications due to their inability to predict the absolute phase of the scattered field. The absolute phase of the scattered field is the fundamental quantity from which the interferogram images are constructed. In this chapter, the present model is for the first time applied to the INSAR simulation on forestry. This investigation wil] demonstrate the capability of the present model in predicting the scattering phase center statistics, provide a physical explanation for the effect of scattering mechanism on the location of the scattering phase center, and discuss the optimum design of INSAR system configuration for remote sensing on vegetation parameters. In the next section, the proposed model for INSAR applications is described. The fractal-based coherent scattering model is briefly reviewed and a Ak-radar equivalence algorithm is employed for the extraction of the scattering phase center. Also in this chapter, the effect of scattering mechanism on the location of the scattering phase center is discussed. The validation of the proposed model is verified by comparing the simulation results with those measured by JPL TOPSAR [59]. Finally a sensitivity study is conducted to demonstrate the variations of the scattering phase center of a forest stand in terms of target parameters such as tree density, soil moisture, tree type, and ground tilt angle, as well as INSAR parameters such as polarization, frequency and incidence angle.

95 5.2 Model Description In this section an overview is given of the approaches which are employed to extract statistics of the scattering phase center of forest canopies. Three tasks must be undertaken for the calculation of the correlation coefficient and the location of the scattering phase center. These include: 1) accurate simulation of tree structures, 2) development of the scattering model, and 3) development of an algorithm for evaluation of the location of the scattering phase center. The first two have been described in detail in the previous chapters. In this chapter, the reemphasis is on the algorithm for evaluation of the location of the scattering phase center. As mentioned earlier, the overall objective of this investigation is to study the relationship between the coherent phase and correlation coefficient of an INSAR interferogram and the physical parameters of a forest stand. An INSAR system measures the backscatter of a scene at two slightly different look angles, and the phase difference between the two backscattered fields is used to derive the elevation information. In a recent study [37] it has been established that similar information can be obtained by measuring the backscatter of the scene at two slightly different frequencies, provided that the look angle is known. For an INSAR system with known baseline distance (B) and angle a operating at frequency fo, the frequency shift (Af) of an equivalent Ak-radar is given by Af - foB sin(a - 0) (5.1) mr where 0 is the looking angle, m = 1,2 for repeat-pass, and two-antenna INSAR con — figurations respectively, and r is the distance between the antenna and the scatterer. This equivalence relationship is particularly useful for numerical simulations and controlled experiments using stepped-frequency scatterometer systems. In Monte Carlo

96 simulations, once the tree structure and the scattering configuration are determined, the backscatter signals are calculated twice at two slightly different frequencies. The backscatter at i = fo and f2 = fo + Af are represented by E1 and E2 respectively, which are computed from N El= e2ikokirnS(ko) E, (5.2) n=1 N E = e2i(ko+k)'rnSr(ko + Ak) E. (5.3) n=l It is also shown that the height of the equivalent scatterer above the x-y plane of the global coordinate system can be determined from -Aq e 2Akcos (5.4) where Ak = 27rAf/c, and A<> = Z(E1E2) represents the phase difference between E1 and E2. Note that the equivalent frequency shift for most practical INSAR configurations is only a small fraction of the center frequency (Af/fo < 0.1%) and therefore the far field amplitudes of individual isolated scatterers (S ) do not change when the frequency is changed from fc to fo + Af, that is, So(ko) ^ S~(ko + Ak). This approximation speeds up the Monte Carlo simulation without compromising the overall accuracy of scattering phase center height estimation. For a random medium like a forest stand, the scattering phase center height (Ze) is a random variable whose statistics are of interest. Usually the mean value and the second moment of this random variable are sought. Based on a rigorous statistical analysis [37] it is shown that the statistics of A<D can be obtained from the frequency correlation function of the target by computing < E*E2 > ceC = 1' (5.5) H IeEll2 >~< E1 212

97 where a is the correlation coefficient, and ( is the coherent phase difference. In (5.5) < * > denotes the ensemble averaging which is evaluated approximately using a sufficiently large number of realizations through the Monte Carlo simulation. When the backscatter statistics are Gaussian, a and ( provide a complete description of the statistics of A1. The apparent height of the scattering phase center of a forest stand is proportional to ( and can be obtained from -C5 Ze = (5.6) 2Ak cos (5.6) Note that ( is not the statistical mean of AD, but rather the phase value at which the probability density function of At> assumes its maximum. In fact, using the mean value may result in a significant error for calculating the apparent height Ze. To demonstrate this, two cases may be considered where in one case ( = 0 and in the other case ( = 180~. In both cases the mean value of AD is zero, whereas the apparent heights calculated from (5.6) are obviously different. 5.3 Verification of Identity Between Ak-radar and INSAR System To test the above Ak-radar algorithms we use a simple configuration, an object above a ground plane, to examine if the scattering phase center derived from a Akradar is identical to that from an INSAR system. However, it is extremely difficult to analytically determine the scattering phase center of this even simple configuration where the multipath scattering components take place due to the ground plane. Here we use an alternative approach in which the comparison between INSAR and Ak-radar is conducted on the basis of each scattering mechanism, whose scattering phase center is distinct according to the path length. In addition, we present a numerical example to demonstrate how the scattering mechanisms influence the

98 A1 A2 \ \ D \\ Di Figure 5.1: Multipath scattering mechanisms of an object (at C with projection at O and image at Ci) above a ground plane observed by an interferometer with two antennas at A1 and A2. scattering phase center. 5.3.1 INSAR Case Consider an object at C above a ground plane measured by an interferometric system with two antennas at A1 and A2, as shown in Figure 5.1. Radar waves are transmitted from antenna Al, scattered and reflected by the object at C and the ground plane respectively, and received by both antennas at A1 and A2. Based on the first order scattering approximation that the object is hit only once during the traveling path, the scattered field is comprised of three scattering mechanisms known as (1) direct backscatter St, (2) ground-bounce backscatter Sg6, and (3) double

99 ground-bounce backscatter Sgtg, that is, S = Steit + Sabei$gb + StYeitgt, (5.7) where Qt, gb, and Ogtg are the associated phases accounting for the different path lengths. As follows are the description of individual scattering mechanisms and the derivation of the correspondent scattering phase centers. 1. S: In this direct backscattering component, the incident field is scattered by the object and immediately returned to the receivers. As shown in Figure 5.1, the traveling path of the scattered field received at antenna A1/A2 is denoted as A1 -- C.- A1/A2, and the scattered fields received are written as S1 = Steiko(2A1C) (5.8) S2 = Steiko(AC+A2C). (5.9) The inner product of these two signals is readily given by S*S2 = S| tl2eko(A2C-A1C) (5.10) The phase term, AI =- ko(A2C - AIC), representing the phase difference between the two antennas, is then processed along with the slant range to derive the 3D digital elevation model. This phase difference implies that the scattering phase center of the direct backscattering component is at C. 2. Sgb:In this ground-bounce scattering component, a reciprocal pair of traveling paths take place simultaneously denoted as Al -* D1 -+ C -+ A1/A2 and A1 -+ C -- (D1 -> A1)/(D2 -~ A2), representing the ground-target and targetground path respectively. The scattered fields received at the two antennas are

100 written as S1 = 2gb iko(AIC+CDI+A1 D), (5.11) S2 = Sgb(eiko(Ac+cD2+A2D2) + eiko(A1D1+CD1+A2C)) (5.12) Noting CD112 = CiD112, the inner product of these two signals can be reduced to SrS2 = 21Sgb12(eiko(A2C-Alc) +eiko(A2Ci-A1Ci)) (5.13) (5.14) Hence, the phase term Ad = ko[(A2C + A2C) - (A1C + ACi)]/2, (5.15) - ko[(A20) -(A10)] (5.16) where O is the projection of the object on the ground plane. Here we use a reasonable approximation, A112C + A1/2Cdi - 2A1/20 since the two antennas are in the far zone. This resultant phase term of the above equation indicates that the scattering phase center of this ground-bounce scattering component is located at 0. 3. Sgt9: In this double ground-bounce scattering component, the incident field is first reflected by the ground plane, then scattered from the object, and again reflected by the ground plane towards the two antennas. The traveling path of the scattered field received at antenna 1/2 is denoted as Al -> D1 -D C -> (D1 -+ A:I)/(D2 -- A2). The scattered fields are written as S1 = 5gtgeiko2(A1D1+CD1) (5.17) S2 = Sgtg iko(A1DI+CD1+CD2+A2D2) (5.18)

101 Figure 5.2: An object above a ground plane illuminated by a Ak-radar. Again, with the identity CD112 = CiD112, the inner product of these two fields can be reduced to 15S2 = IS9tl2ei'k(A2Ci-ACi) (5.19) The phase term of the above expression indicates that the scattering phase center of the double ground-bounce scattering component is located at Ci. 5.3.2 Ak-Radar Case Consider an object above a ground plane with height h illuminated by a plane A wave in the direction of ki, as shown in Figure 5.2. The total backscattered field can be expressed as eikor Es = (St + Sgb + S9tg)e-2ikohcos (5.20) with multipath scattering components: with multipath scattering components: St - So(-k, i,) (5.21)

102 S = [R. So- kr, ki) + So( —ki, kr) R]e2~cos (5.22) Sgtg - [R -So(-kr, kr)' R]e4ikohcos (5.23) where 5o is the bistatic scattering amplitude of the object in free space, and R is the reflection coefficient of the ground plane. In general So and R are in the form of matrices for the polarimetric analysis, and they are assumed invariant while the frequency is slightly changed (e.g. Ak/ko == 10-4 for JPL TOPSAR). Now we examine the scattering phase center for each component: 1. S': In this case the scattering amplitude from (5.2) is E1 = So(-ki, ki)-2ihc (5.24) Using the frequency decorrelation process, the phase difference is obtained as A<1 = phase(E*E2) = -2A\kh cos 0. Therefore, from (5.4) the scattering phase center is obtained as Ze = h. 2. Sgb: In this case the scattering amplitude from (5.2) is E1 =R So(-kr, ki) +- So(-ki, kr) R (5.25) The phase difference is AD = 0 since there is no path-length dependent phase term associated. Therefore the scattering phase center in this case is on the ground, z,e = 0. 3. Sgt9: In this case the scattering amplitude from (5.2) is written as E1 [[R S0(-kr, kr). R]e2ikohcos ~ (5.26) The phase difference can be obtained as Az(I = 2Akh cos 0. Therefore, from (5.4) the scattering phase center is obtained as Ze = -h, the image of the object.

103 Based on the above analysis, Ak-radar and INSAR are shown to preserve the same scattering phase center for each scattering component of the multipath mechanisms. Since the total backscattered field is the coherent addition of these components, it can be inferred that the equivalent scattering phase center of the whole system derived from zAk-radar and INSAR is identical. 5.4 Scattering Mechanisms vs. Scattering Phase Centers In order to develop some intuition about the scattering phase center and to examine the validity of the above equivalence algorithms, let us consider a simple case where a single scatterer target above a ground plane is considered. Through this illustrative case the relationship between the location of the scattering phase center and the scattering mechanisms can be demonstrated. Consider a dielectric cylinder of radius a = 5cm, length b = 3m, and dielectric constant ct = 22 + 10, which is located at height h = 6m above a ground plane having a complex permittivity cg = 9.7 + il.6. Suppose the target is illuminated by a plane wave whose direction of propagation is determined by the incident angles Oi = 30~, 4i = 180~, as shown in Figure 5.3. As mentioned previously, the backscattered field is mainly composed of four scattering components with different path lengths. In general it is quite difficult to characterize the location of the scattering phase center of a scatterer analytically when multi-path scattering mechanisms are involved. However, in cases where a single scattering mechanism is dominant it is found that the location of the scattering center is strongly dependent upon the path length of the dominant scattering component. Here we illustrate this fact through an experimental study where the orientation of the cylinder is properly arranged in four configurations, as shown in (a)-(d) of

104 (a) Z St A 4 I h Ze z (b) + I: A h X (C) Z ~~~ --- —S gtg I / r Figure 5.3: Four configurations for a cylinder above a ground plane with orientation angles: (a) 0, = 60~,+ = 180~, (b) 0, = 0~, = 0~, (c) 0, = 60~, c = 0~, (d) Oc = 45~, c = 150~, and their principal scattering mechanisms respectively. The center of the scatterer is denoted by (e) and equivalent scattering phase center by (D').

105 Configuration (a) (b) (c) (d) Polarization vv hh vv hh vv hh vv hh ze/h 1.00 0.99 0.01 -0.01 -1.06 -0.96 0.49 0.43 IS/SI 0.99 0.97 0.02 0.01 0.05 0.05 0.62 0.42 ISb1SI 0.03 0.02 0.99 1.01 0.13 0.08 0.62 0.59 iSt_ / S 0.00 0.00 0.01 0.00 1.10 1.03 0.00 0.00 RCS (dBsm) 8.06 5.19 -0.46 6.16 -6.05 -5.23 -17.2 -15.1 Table 5.1: The normalized height of the scattering phase center, the normalized scattering components, and the radar cross section for four different scattering configurations as shown in Figure 5.3. Figure 5.3, such that the total backscatter is dominated by (a) St, (b) Sgb, (c) Sgtg, and (d)St + S:b respectively. Note that Sgb is the combination of the reciprocal pair Sgt and Stg. The simulation results at fo =1.25 GHz are shown in Table 5.1, which includes the scattering phase center height normalized to the physical height Ze/h, the ratio of the amplitude of individual scattering components to the total backscattered field IS()/SI, and the overall radar cross section (RCS) of the target for each orientation configuration and for both polarizations. It is obvious from the results reported in Table 5.1 that in scattering configuration (a) where the backscatter is dominated by the direct component (St/S = 0.99) the location of the scattering phase center appears at the physical location of the scatterer above the ground (ze/h _- 1). Similarly in scattering configurations (b) and (c) where the backscatter is dominated, respectively, by the single ground bounce component and the double ground bounce component, the locations of the scattering phase center appear on the ground surface and at the mirror image point as shown in

106 Figure 5.3. In scattering configuration (d), the direct and the single ground bounce components of the backscatter are comparable in magnitude and as shown in Table 5.1 the location of the scattering phase center in this case appears at a point between the physical location of the scatterer and the ground surface. When the number of scatterers is large the location of the scattering phase center is a convoluted function of physical locations of the constituent scatterers and the relative magnitudes and phases of the scattering components. Note that \Ze/hI and ISH)/S\ in Table 5.1 may exceed 1 since the total backscattered field S is the superposition of four scattering components which are not necessarily in phase. 5.5 Model Validation for INSAR Applications In this section full simulations of forest stands are carried out. As a first step, the model predictions are compared with the JPL TOPSAR measurements over a selected pine stand, denoted as Stand 22. Then a sensitivity study is conducted to characterize the variations of the scattering phase center height and correlation coefficient as a function of both forest and INSAR parameters. Stand 22 is a statistically uniform red pine forest located within Raco Airport, Raco, Michigan. This scene was selected for this study because the stand is over a large flat terrain which reduces the errors in the measured tree height due to possible surface topographic effects. In addition, the nearby runway provides a reference target at the ground level. Ground truth data for this stand have been collected since 1991 [54] and careful in-situ measurements were conducted by the authors during the overflights of TOPSAR in late April, 1995. The relevant physical parameters of this stand are summarized in Table 5.2. The vegetation and soil dielectric constants are derived from the measured moisture contents using the empirical models described

107 Runway Stand 22 Figure 5.4: A portion of a TOPSAR C-band image (cr), indicating Stand 22 (red pine) at an airport near Raco, Michigan (April 1995). in [56,57]. The JPL TOPSAR is an airborne two-antenna interferometer, operating at Cband (5.3 GHz) with vv polarization configuration [59]. During this experiment, Stand 22 was imaged twice at two different incidence angles 39~ and 53~. Figure 5.4 shows a portion of the 39~ radar image which includes the test stand. Each side of the dark triangle in this image is a runway of about 2 miles long. The measured height of the stand is obtained from the elevation difference between the stand and the nearby runway. Using the ground truth reported in Table 5.2, the backscattering coefficient and the location of the scattering phase center as a function of the incidence angle were simulated at 5.3 GHz. As shown in Figures 5.5 and 5.6, excellent agreement be

108 Tree Density: Tree Height: Trunk Diameter (DBH) Needle Length: Needle Diameter: Needle Moisture (mg): Wood Moisture (mg): Soil Moisture (mV) Table 5.2: Ground truth data of Stand 22 1142#/Hectare 8.9 m: 14.6 cm 10 cm 1.2 mm 0.62 0.42 0.18 (red pine) collected on April 26, 1995. = tween the model predictions and TOPSAR measurements is achieved. The simulated height of the scattering phase center of the same forest for an hh-polarized INSAR, having the same antenna configuration and operating at the same frequency is also shown in Figure 5.5. It is shown that the estimated height at the hh-polarization configuration is lower than that obtained from the vv-polarization configuration. This result is usually true for most forest stands since the ground-trunk backscatter for hh-polarization is higher than that for vv-polarization because hh-polarized wave has higher transparency within the canopy and stronger reflection from the dielectric plane. Also noting that the location of the scattering phase center for a ground-trunk backscatter component is at the air-ground interface, the location of the scattering phase center of trees for hh-polarization is lower than that for vv-polarization. The comparison between the simulated aco and the measured ao acquired by TOPSAR as a function of the incidence angle is shown in Figure 5.6. Also shown in this figure is the contribution of each scattering component (the direct backscatter at and the ground-bounce backscatter cab) to the overall backscattering coefficient.

109 10...... I... I... I.. * TOPSAR (vv-pol) o8. Model (vv-pol) N[ 6 -— e -- Model (hh-pol) 6 i4 ~-i C4 —. 4 --. - 10 20 30 40 50 60 70 Incidence Angle Oi (Degrees) Figure 5.5: The estimated height of scattering phase center of Stand 22 (red pine), compared with the data extracted from two TOPSAR images (0 = 39~ and 53~) of the same stand.

110 It was found that the contribution of the double ground-bounce component agtg was relatively small and for most practical cases can be ignored. In this case, at low incidence angles (0i < 30~) the ground-bounce backscatter is the dominant component, whereas at higher incidence angles the direct backscatter becomes the dominant factor. This trend is the cause of the increasing behavior of the scattering phase center height as a function of the incidence angle found in Figure 5.5. It is worth mentioning that the contribution of pine needles to the overall backscattering coefficient is found to be negligible compared to the contribution from the branches and tree trunks. However, inclusion of the needles in the scattering simulation is necessary because of their significant effect on the extinction. 5.6 Sensitivity Study With some confidence in the scattering model and the algorithm for evaluation of the scattering phase center height, further simulation can be performed to characterize the dependence of the scattering phase center height of a forest stand on the system parameters such as frequency, polarization, and incidence angle, and the forest parameters such as tree density, soil moisture, and tree types. In addition, we demonstrate the capability of the present model as a tool for determining an optimum system configuration for retrieving physical parameters of forest canopies. Figure 5.7 shows the estimated height of Stand 22 for two principal polarizations at C-band (5.3 GHz) and L-band (1.25 GHz) as a function of the ground soil moisture, simulated at 0i = 45~. As the soil moisture increases, the ground plane reflection will also increase, which in turn causes the ground bounce scattering component to increase. As a result of this phenomenon, the scattering phase center height decreases with soil moisture as shown in Figure 5.7. This effect is more pronounced for L-band

111 Illi o -10 U -20 - * TOPSAR *) ---- Model (total)\ u -30 - -— a --- Component o C — A — Component ogb -40.... a........ i.... I.... i.... -40.... 10 20 30 40 50 60 70 Incidence Angle 0i (Degrees) Figure 5.6: The simulated backscattering coefficient of Stand 22 (red pine), compared with the measured ao extracted from two TOPSAR images of the same stand.

112 vv-polarization than other INSAR configurations, suggesting a practical method for monitoring the soil moisture using the apparent height of the forest stand. This high sensitivity at Lvv is achieved because of the existence of competitive scattering components. Basically, at low soil moisture the direct backscatter component is comparable with the ground bounce component and the scattering phase center lies amidst the canopy. As the soil moisture increases, the ground bounce scattering component becomes more dominant, which results in lowering the apparent height of the stand. On the other hand, the least sensitive configuration is Lhh since the dominant scattering component, independent of the soil moisture, is the ground bounce component. Figure 5.8 shows the effect of the tree density on the estimated height of a red pine stand having a similar structure as that of Stand 22 at Oi = 45~. As the tree density increases, the extinction within the canopy increases, which reduces the ground-bounce component. Increasing the tree density would also increase the direct backscatter component. As a result of these two processes, the apparent height of the canopy increases with increasing tree density as demonstrated in Figure 5.8. As before, the apparent height for Lhh configuration does not show any sensitivity to the tree density indicating that the ground bounce component remains dominant over the entire simulation range of 700-1200 trees/Hectare. This lack of sensitivity to the apparent height of coniferous stands for Lhh suggests that this configuration is most suitable for mapping the surface height of coniferous forest stands. Now let us examine the response of INSAR when mapping deciduous forest stands. For this study a red maple stand, denoted by Stand 31, is selected whose structure and scatterers are different from the previous example. A fractal generated red maple tree and a picture of the stand are shown in Figure 2.9. This stand

113 6-,\ 5 - N to 3 d0 a 2 0-0 co. I o 0.0 ---— O --- - -o I I I I I ---— c ----- L ----- Lhh hh — I --- Lv -— v --- Lh 0.1 0.2 0.3 0.4 0.5 0.6 Soil Moisture, m, Figure 5.7: The estimated scattering phase center height of Stand 22 (red pine) as a function of soil moisture, simulated at Oi = 450~

114 6 5 a) N 4 —) b' ct E r.~ *4-U CA^ w 4 3 2 1 0 I ' I I I I A - -- - - -: ^.- ------— J... --- —-----— "..... -....... -- -" V...... V-... — ^.........V --- — -- ^ -0* ----v ----- L 600 700 800 900 1000 1100 1200 1300 Tree Density, (#/Hectare) Figure 5.8: The estimated scattering phase center height of Stand 22 (red pine) as a function of tree density, simulated at 0, = 45~.

115 was selected as a test stand to validate the previously developed coherent scattering model [60], using the SIR-C data. The average tree height and tree number density were measured to be 16.8m and 1700 trees/Hectare respectively. Table 4.1 provides the detailed ground truth data from Stand 31. The simulations for estimating the scattering phase center height are performed fully-polarimetrically at L-band and Cband. Figure 5.9 shows the variation of the apparent height of Stand 31 as a function of the incidence angles for co- and cross-polarized L- and C-band INSAR configurations. Simulation results at C-band show that except at very low angles of incidence, the scattering phase center is near the top of the canopy. In this case the backscatter in all three polarizations is dominated by the direct backscatter components of particles near the canopy top. The same is true for L,, and Lsh configurations; however, since penetration depth at L-band is higher than C-band, the location of the scattering phase center appears about 1-3 m below the apparent height at Cband. The scattering phase center height for L;h configuration, on the other hand, is a strong function of the incidence angle where it appears near the ground surface at low incidence angles and increases to a saturation point near grazing angles. Close examination of this figure indicates that a pair of C,, and Lhh INSAR data at low incidence angles can be used to estimate the tree height of deciduous forest stands with closed canopies. A C-band foliated canopy behaves as a semi-infinite medium and as shown in [37] the knowledge of extinction would reveal the distance between the location of the scattering phase center and the canopy top (Ad) using Ad = cos O/(2K). If an average extinction coefficient (K) of 0.2Np/m is used in the above equation, a distance Ad = 1.77mn is obtained at 0 = 45~. However, a simple relation for evaluating the apparent height for Lhh does not exist yet. Figure 5.10 shows the effect of the ground tilt angle on the estimated scattering

116 20............ ~S 15 |- r-' — -tA - -- - = 10 F- -> l —a — C 1- - L Lvv L ---.C. v -I --- -C, 0 Chh 10 20 30 40 50 60 70 Incidence Angle Oi (Degrees) Figure 5.9: The estimated scattering phase center height of Stand 31 (red maple) as a function of incidence angle, with fully-polarimetric L- and C-and response.

117 20 ' ' 15 - -—? — - - E — - -. - --- - \- -- - - L 'U; CVV 5. - - -E. c --— e ---* Chh.)....... hh vv 20 40 60 80 100 120 Azimuthal Look Angle, Xi (Degrees) Figure 5.10: The estimated scattering phase center height of Stand 31 (red maple) over a tilted ground with tilt angle 0g = 10~, at,O = 25.4~. phase center height. This simulation is obtained by setting a ground tilt angle Og = 10~ for a forest stand similar to Stand 31 (red maple) and calculating the estimated scattering phase center height as a function of the azimuthal incidence angle qi at O- = 25.40. As already discussed in Figure 4.14, there is a strong groundtrunk backscatter around Xi = 700, particularly for Lhh and Lvh. Consistently, this accounts for the dip in the apparent height simulations for Lhh and Lvh configurations at oi = 700, as shown in Figure 5.10, noting that the scattering phase center of the ground-target scattering component is on the ground plane. So far only the behavior of the mean value of the scattering phase center height

118 has been investigated; however, the model has the ability to provide an approximate probability distribution function of the scattering phase center height. The histograms of the scattering phase center height can be constructed by recording the simulated results for each scattering simulation. Figure 5.11 shows the simulated probability density function (PDF) of the scattering phase center height of Stand 22 at 0 = 45~ for the three principal polarizations and for both L- and C-band. At Lband the scattering phase center height has a narrow distribution for hh-polarization, indicating that a relatively small number of independent samples are sufficient for estimating the apparent height. At C-band the scattering phase center height of the cross-polarized backscatter exhibits a narrower PDF. As mentioned earlier, the correlation coefficient (a) is an independent parameter provided by INSAIRs which, in principle, may be used for inversion and classification processes. The measured correlation coefficient is a function of INSAR parameters such as look angle, baseline distance and angle, radar range to target and target parameters. To examine the behavior of a as a function of target parameters, the Ak-radar equivalence relationship given by (5.5) is used where the dependence on INSAR parameters are lumped into one parameter, namely, the frequency shift. Figures 5.12 and 5.13 respectively show the calculated correlation coefficients (a) as a function of the normalized frequency shift (Af/fo) (corresponding to the baseline distance in an INSAR), simulated for Stand 22 (red pine) and Stand 31 (red maple) at Q- = 45~. As shown in [37] the correlation coefficient is inversely proportional to the width of the PDF, that is, a high value of a indicates a narrow distribution. A comparison between the histograms shown in Figure 5.11 and the values of a shown in Figure 5.12 demonstrates this relationship. It is interesting to note that simulated a for Stand 31 at Lhh is significantly smaller than the correlation coefficients at other

119 L-band ~ 0.2.- 0.1 Q. VL ~ 0.2 1 - 0.1 ru C-band n II fL IIrrrfln I O-10 0 10 20 U-1 — 10 0 10 20 0. 0.2 I c-,-0.1 VL I I n 0.O..c.~ 0. - O. a 0L.2.1 ('. n m ( -10 0 10 20 — 10 0 10 20 0 0.0.2 I -C C0.1,-0.1 V. -10 0 10 20 Ze (m) Ze (m) Figure 5.11: PDF of the scattering phase center height of Stand 22 (red pine) at Oi == 45~ as a function of frequency and polarization.

120 Stand 22 1.0 0.8 0.6 0.4 0.2 -0 — --- E ---* Lvv Lvh Lhh Cvv Cvh Chh 0.0 - 0.0 2.0 4.0 6.0 8.0 10.0 Frequency Shift Af/f0 (x 10-3 ) Figure 5.12: The correlation coefficient as a function of the frequency shift, simulated from Stand 22 (red pine). polarizations (see Figure 5.13). This behavior is a result of the fact that the direct backscatter and ground-bounce backscatter components are comparable. It is shown that for the same baseline to distance ratio (B/r) which corresponds to a constant Af/fo,? a at C-band is smaller than a at L-band independent of polarization. It should be mentioned here that for most practical situations Af/fo is of the order of 10-4 or smaller which renders a value for a near unity (a > 0.99). That is, for practical INSAR configurations, the effect of forest parameters on the correlation coefficient appears on the third digit after the decimal point. This puts a

121 Stand 31 1.0, - - 0.6.0. ----.vh. |Qj 0.4 —.,-. —Lh 0. 2.0 4.0 6.0 8.0 10.0 0.2 Cvh ------- Chh 0.0 ' ' ' ' 0.0 2.0 4.0 6.0 8.0 10.0 Frequency Shift Af/fo (x 103 ) Figure 5.13: The correlation coefficient as a function of the frequency shift, simulated from Stand 31 (red maple).

122 serious limitation on the applicability of a for inversion and classification algorithms, since accurate measurement of a with three significant digits is not practical even with two antenna INSARs. For repeat-pass interferometry, the a values reported for forest stands are below 0.7 which is caused mainly by the temporal decorrelation of the target. Therefore, it does not seem logical to use a as a parameter for classifying forest types. The TOPSAR measured as for Stand 22 at incidence angles 39~ and 530 are, respectively, 0.935 and 0.943 which are below the calculated values of 0.998 and 0.999. This discrepancy can be attributed to processing errors and thermal noise. 5.7 Conclusions In this chapter, a scattering model capable of predicting the response of interferometric SARs when mapping forest stands is described. The model is constructed by combining a first-order scattering model applied to fractal generated tree structures and a recently developed equivalence relation between an INSAR and a Ak-radar. Using this model, for the first time accurate statistics of the scattering phase center height and the correlation coefficient of forest stands are calculated numerically. The validity and accuracy of the model are demonstrated by comparing the measured backscattering coefficient and the scattering phase center height of a test stand with those calculated by the model. Then an extensive sensitivity analysis is carried out to characterize the dependence of the scattering phase center height on forest physical parameters, such as soil moisture, tree density, and tree types, and INSAR. parameters such as frequency, polarization, and incidence angle. The ability of the model to predict the PDF of the scattering phase center height and the correlation coefficient is also demonstrated. It is shown that for practical INSAR configurations, the correlation coefficient of forest stands is near unity, much larger than what can

123 be measured by existing INSAR systems.

CHAPTER VI RETRIEVAL OF FOREST PARAMETERS USING A FRACTAL-BASED COHERENT SCATTERING MODEL AND GENETIC ALGORITHMS 6.1 Introduction Retrieval of gross biophysical parameters of forest stands, such as basal area, tree height, and leaf area index (LAI), is of great importance in many environmental research programs. Radar remote sensing at lower microwave frequencies has been proposed as a sensitive instrument for such applications [14,15]. In support of pro — grams pertaining to radar remote sensing of vegetation, many advanced polarimetric (SIR-C, AIRSAR) [53] and interferometric (TOPSAR) [59] radar instruments have been developed. The study of the inversion problems in geophysical science and engineering has been of great importance from the onset of the remote sensing science [61,62]. For example, in microwave remote sensing of vegetation the inverse problem is defined as the application of the measured quantities such as the polarimetric backscattering coefficients (from a SAR) [45] and/or the scattering phase center heights (from an interferometric SAR) [37,38] in an algorithm in order to retrieve forest parameters 124

125 such as tree type, tree density and height, and moisture content of vegetation and soil. Over the past two decades significant effort has been devoted towards the development of scattering models for vegetation canopies [21,26-29,60] as well as inversion models to retrieve forest parameters from the measured data [63-65]. So far the emphasis of the scattering model development has been on the construction of simplified models with as few input parameters as possible so that the inversion problem becomes tractable. In this process the importance of structural features of the canopy (particle arrangement), coherence effects, and multiple scattering were ignored. Even with these simplifications, the inversion process is rather complex. In [63,65] neural network approaches are suggested for the inversion process where extensive computer simulations or experimental results are used to train a neural network in a reverse order (the model outputs are fed as the input to the program). This method is computationally extensive and its success depends on the fidelity and the extent of the training data. In [64] a gradient-based search routine is applied to a nested linearized model. This model is comnputationally efficient; however, its applicability is limited to models with small dimensionality and its success depends on the fidelity of the forward model. This chapter describes the application of a high fidelity scattering model in an inversion process based on a stochastic global search method. Basically, the aforementioned coherent scattering model that preserves the structural features of tree canopies using fractal models is employed to generate simplified empirical models (for different tree species) that can predict the polarimetric and interferometric radar response of a forest stand efficiently and accurately. The premise for the successful development of such empirical models stems from the fact that the model outputs

126 are averaged quantities, such as backscattering coefficients or the mean height of the scattering phase center, and therefore are very gentle functions of model inputs. As demonstrated in the previous chapters, the fractal-based coherent scattering model (FCSM) offers two advantages over the traditional scattering models, namely, FCSM is more versatile and accurate. Basically, FCSM is capable of simulating the fully polarimetric (including the phase statistics) and polarimetric-interferometric (scattering phase centers and correlation coefficients for any polarization configuration [39]) radar responses of coniferous and deciduous forest stands. High accuracy is achieved by FCSM through incorporating the coherent effects among the individual scatterers and scattering components and by accounting for the accurate position of scatterers which is manifested in inhomogeneous scattering and extinction profiles. However, this versatility and accuracy has been achieved at the expense of the model complexity which demands extensive computational power. For example, the number of input parameters needed to accurately characterize the tree structures and the environment may easily exceed 30 (it should be noted that once a tree type is chosen much fewer free parameters are needed to model the natural variabilities). On the other hand, to obtain a solution with a reasonable accuracy in the Monte Carlo simulation, a sufficiently large number (> 100) of realizations are required. The required computational time for each simulation limits the model's utility in inverse processes which may demand the calculation of the forward problem many times. To circumvent the aforementioned problem, development of empirical models based on FCSM is proposed. Construction of an empirical model can be achieved using a standard procedure such as curve-fitting and regression method. Unlike physical models, empirical models are simple mathematical expressions formed from a set

127 of data acquired from measurements or a physical model prediction. Once empirical formulae are obtained, they are easy to use and require minimal computation time. However., it should be noted that an empirical model is usually valid only for a specific case within a certain range of the parameter space over which the model is constructed. For the development of the empirical model used in this study, first a sensitivity analysis is conducted in order to determine the significant parameters, the number of which determines the dimensionality of the input vector space. A red pine stand is chosen in this study and six parameters are selected as the input parameters. Each selected parameter is allowed to have about 30% variation with respect to a centroid. Using the Monte Carlo simulation results obtained from FCSM a database is constructed by varying the individual parameters over a prescribed range of the input vector space around the centroid. The parameters at the centroid are obtained from the ground truth data of a red pine test stand (Stand 22) in Raco, Michigan. For the inversion process, first a least-square estimator is used and is shown to work properly when the number of measured channels is equal to or larger than the dimension of the input vector space. But since this may not be the case in general situations, a genetic algorithm (GA) [66] is developed and employed as a search routine for the nonlinear optimization problem. GAs are known to be very successful when the dimension of the input vector space is large and/or when the objective function is nonlinear.

128 6.2 Empirical Model Development In general, the output of the Monte Carlo coherent scattering model can be expressed as M = (f, p, 0; da, Ht, Dt. Ob, ms, mw), (6.1) where C is a complex operator relating the input and output of the model and the output M is a vector which may contain the backscattering coefficient (a 0, Crh, Chh), scattering matrix phase difference statistics, the scattering phase center height Ze, or the interferogram correlation coefficient. The input parameters are divided into two categories: 1) radar system parameters, and 2) target parameters. Radar parameters include the radar frequency f, the polarization configuration p, and the incidence angle 0. The number of target parameters can be very large, consisting of the tree structural parameters and the dielectric properties of the constituent components. However, the number of these parameters is reduced drastically once a tree type is chosen. In this case only a few structural parameters are sufficient to allow for natural variabilities observed for that type of tree. The rest of the structural parameters are embedded in the fractal code of the tree. In this study we demonstrate development of an empirical model for a red pine tree where only six free parameters are sufficient to describe the stand. These include the trunk diameter da, tree height Ht, tree density Dt, branching angle Ob, soil moisture ms and wood moisture m". It should be noted that these parameters themselves are statistical in the coherent model with prescribed distribution functions and here we are referring to their mean value. Multi-frequency polarimetric SAR systems operate at discrete frequencies, usually at P-, L-, C-, and X-band, and the polarization configuration p are vv, vh, and hh. In this study, we demonstrate a model with three fundamental backscattering

129 coefficients and the associated mean scattering phase center height as the model output and fix the frequency at C-band (5.3 GHz). The empirical model is developed to operate over the angular range 25~- 70~. Therefore the output and input vectors M and x are defined as: zev da Zvh Ht L =, and x= 0 0 'vh hh m As mentioned earlier since no resonance behavior is expected, the output vector M is a gentle function of the input vector x and the incidence angle 0 which may be related to each other via a simple empirical relationship M = L(0; x), (6.2) where L is the simple empirical operator and M is the output of the empirical model. It is expected that M be as close to M as possible. In general, the output parameters are non-linear functions of the incidence angle and other input parameters. In order to establish these relationships the coherent model must be run by varying the incidence angles and other input parameters. Through an extensive sensitivity study it was found that over a finite domain of the input vector space a logarithmic relationship between the backscattering coefficient (linear in dB scale) and a linear relation between the scattering phase center height and the input parameters exist. The dependence on the incidence angle was found to be nonlinear.

130 Parameter Range Variation Trunk Diameter 10.9 14.9(cm) 31 % Tree Height 8.0 9.9(m) 21 % Tree Density 807 - 1027 (trees/Hectare) 24 % Branching Angle 13.8 - 15.8 (deg) 14 % Soil Moisture 0.36 -0.56 (g/g) 43 % Wood Moisture 0.28 - 0.48 (g/g) 52 % Table 6.1: The ranges of the selected ground truth parameters and the corresponding percentage variations to the centroid. The first step in the construction of the empirical model is to choose the domain of the input vector space. In this investigation we chose the structural parameters of a young red pine stand, a test forest stand in Raco Michigan (Stand 22), and the seasonal average of soil and vegetation moisture as the centroid of the input domain. These parameters and their range of variation used in the model development are shown in Table 6.1. The Monte Carlo simulation was then carried out for specific incidence angles by varying the six free parameters within the prescribed ranges. The average scattering phase center heights (Ze) for each polarization configuration and backscattering coefficients (a~ in dB) are shown in Figures 6.1 and 6.2 as a function of each parameter respectively. These figures clearly demonstrate the linear relationship previously described. Hence the output vector can be readily approximated by the Taylor series expansion of the exact model to the first order, and is given by ~(x) = L(xo) +A. (x - xo) (6.3) where xO denotes the input vector at the centroid and A is the matrix of partial

131 (a) (b) (c).. I 5 - - - -*- -. — E 0) N 1 -E a) N N N) 3 2 1 0 n — K -- -— T — ek-. -i *- * 8 8.5 9 9.5 Ht (m) -1 0.085 0.09 0.095 0.1 Dt (m-2) 11 12 13 14 dbh (cm) (d) (e) F -..I (f) ~,, 5 - 5 4 4 - -*- -*- - — * — - * - *- - -*- - 3 3F E - a) N - a) N r ---IN- |2. 1 0 -1 - -' — ~- -- --.- ~ - ~ 0.4 0.45 0.5 0.55 ms (g/g) 1 0 -1 -.-. - -. - ~* - ' — 0.3 0.35 0.4 0.45 m (9/9) 14 14.5 15 0b (deg) 15.5 Figure 6.1: The sensitivity analysis of the C-band polarimetric scattering phase center height as a function of the physical parameters: (a) trunk diameter, (b) tree height, (c) tree density, (d) branching angle, (e) soil moisture, and (f) wood moisture, simulated at incidence angle 0 = 25~.

132 (a) -.* - * -. - *- -i HH VV (b) (c) -10 -c I 15 0-1C s — -15. -2C ) m-10 0t t-1 -15 * * - tI - -* — - *- - VH -St- - i -20 11 12 13 14 dbh (cm) — 4 -.. "' -, 8 8.5 9 9.5 Ht (m) -20 - - --- 0.085 0.09 0.095 0.1 Dt (m-2) (f) M (d) (e) -. -* - - - w -, m - -- C) t) *10 m 8 -o0 -10- *- * --- -15 -201 - ---- - -- - 14 14.5 15 15.5 0b (deg) -15 -20 m -10 -0 0 -15 -20 ~ ---- - *- - -*- -. 0.4 0.45 0.5 0.55 m (g/g) S * * * I~t~ - - -* — * - --- - - 0.3 0.35 0.4 0.45 mw (g/g) Figure 6.2: The sensitivity analysis of the C-band polarimetric backscattering coefficient as a function of the physical parameters: (a) trunk diameter, (b) tree height, (c) tree density, (d) branching angle, (e) soil moisture, and (f) wood moisture, simulated at incidence angle 0 = 25~.

133 derivatives whose ij-th element is given by 9C; aij= - Ix=xo. (6.4) aij simply represents the derivative of the i-th output channel LC with respect to the j-th input parameter xj, evaluated at the centroid xo. In this matrix, each element was evaluated by calculating the slope of a fitting line over 5 sample points based on a least square method. In Figures 6.1 and 6.2 the symbols (*) are the simulation results and the lines are the best linear estimation. It should be pointed out that each point in each figure represents an ensemble average of 200 realizations of the Monte Carlo simulation. This indicates that the initial task of generating a matrix of coefficients is very tedious and time-consuming. However, once the empirical model is obtained, it can provide a highly accurate solution to an arbitrary input in almost real-time. This property of the empirical model is especially important in the inversion processes. Results in Figures 6.1 and 6.2 are for a fixed incidence angle 0 = 250. However, the simulations at other incidence angles show that the general form of (6.3) is valid for all incidence angles with the exception that ~(xo) and A are functions of the incidence angle., i.e., ~(0; x) = ~~(0) + A(0). (x - xo). (6.5) It is found that OC(0) and A(0) are non-linear, but gentle, functions of the incidence angle 0 over the range of interest (250 to 70~). In order to obtain the functional form of L~ and A on 0, the aforementioned Monte Carlo simulation was repeated at several different incidence angles, and the corresponding values L~ and A were evaluated. Polynomial functions are used to capture the angular variations of L~ and

134 A. It was found that ~o and A can be accurately expressed by 20(0) = 2O + ~10 + 2202 + 4303, (6.6) and A(0) = Ao + A10 + A202 + A303 + A404, (6.7) where Li and Ai are 6 x 1 and 6 x 6 matrices whose values are given by 4C 0 1 12 C 3 ] -3.1051 3.6061 3.6264 -12.2843 -50.1483 3.5502 0.2642 -0.0659 -0.3188 0.2713 2.1623 -0.3062 -0.0020 0.0029 0.0094 -0.0079 -0.0431 -0.0005 0.0000 0.0000 -0.0001 0.0000 0.0003 0.0000 and Ao = 1.5294 — 0.6424 7.1314 -13.5348 — 1.0359 — 4.7396 -0.2244 0.2687 0.1875 -1.8434 -1.0341 -2.3518 49.0879 -0.5762 10.3344 -45.9425 22.6264 -224.1260 10.8663 0.7727 — 7.0854 1_6.2122 5.2748 — 6.9444 5.2289 -2.4209 -3.6540 2.2269 6.6677 -4.8994 -59.6487 -7.0918 -2.5728 -108.8944 -73.7306 114.7993

135 -0.1943 0.0320 -4.2781 -0.9772 -0.7657 5.5036 0.0581 -0.0244 0.0609 -0.0693 0.1572 0.5872 -0.7710 -0.0071 -0.5252 0.6810 0.3879 0.3659 1.3292 0.1763 3.9880 -1.3883 0.3188 11.9278 0.0508 0.1137 -2.1445 -0.2531 -0.5110 6.0220 0.6771 0.1919 19.5566 0.7884 0.9928 -10.1014 0.0086 -0.0013 0.1302 0.0316 0.0312 -0.1844 -0.0019 0.0010 -0.0025 0.0029 -0.0038 -0.0198 0.0297 -0.0002 0.0005 -0.0233 -0.0161 -0.0121 -0.0469 -0.0057 -0.1266 0.0433 -0.0218 -0.4499 -0.0010 -0.0040 0.0673 0.0035 0.0144 -0.1697 -0.0305 -0.0052 -0.6039 -0.0325 -0.0354 0.3303 -0.0002 0.0000 -0.0017 -0.0004 -0.0005 0.0027 0.0000 -0.0000 0.0000 -0.0001 0.0000 0.0003 -0.0005 0.0000 0.0002 0.0003 0.0003 0.0001 0.0007 0.0001 0.0017 -0.0006 0.0004 0.0072 0.0000 0.0001 -0.0009 -0.0000 -0.0002 0.0021 0.0005 0.0001 0.0079 0.0006 0.0005 -0.0047 Z0.0097 -0.0013 0.0795 0.0222 0.0306 -0.1417 -0.0014 0.0009 -0.0025 0.0033 -0.0016 -0.0148 0.0278 -0.0008 -0.0218 -0.0176 -0.0142 -0.0044 -0.0390 -0.0038 -0.0831 0.0302 -0.0284 -0.4111 -0.0003 -0.0033 0.0469 -0.0002 0.0079 -0.0906 -0.0339 -0.0023 -0.3794 -0.0357 -0.0218 0.2407 10-4

136 Figure 6.3 compares the results of the empirical model given by (6.5) with those of the Monte Carlo simulation at the centroid (x = xo). It should be noted that the choice of the output parameters are arbitrary and depends on the available set of input data. For example, an empirical model for a two-frequency system with three backscattering coefficients could be developed using the same procedure. Equation (6.5) represents the overall empirical model whose accuracy can be evaluated through a comparison with the Monte Carlo simulations. For this comparison a large number Monte Carlo simulations with independent input vectors were carried out. Figures 6.4, 6.5, and 6.6 show the comparison between the results of the empirical model and those of the Monte Carlo coherent model using 200 independent input data sets randomly selected within the aforementioned domain of the empirical model. The figures show excellent agreement between the empirical model and the Monte Carlo coherent model, noting that the convergence criteria for the Monte Carlo model is ~0.5 dB. Having confidence on a fast and accurate empirical model, the inversion processes can be attempted which is the subject of the next section. 6.3 Inversion Algorithms Consider a physical system whose input-output relation is expressed by M = ~(x) where in general M and x are multi-dimensional vectors of arbitrary length. The inverse problem is mathematically defined as x = ~-1 (M) subject to certain physical constraints. Although the inverse problem may be well-defined mathematically, in practice the inverse solution may not exist for two reasons: 1) mathematical construction of the model may not be exact, and 2) the measured vector M may not be exact because of measurement errors. Hence, instead of casting the problem in terms of an inverse problem, the problem of finding x is usually cast in terms of a

137 v- - C(4 N 0 -3 - 4 2 I, 1 -F 20 30 40 50 60 70 0 (degrees) (a) VV +- 0 VHX 0 s + * a. -15 N,~-'+ HH -1o0 -15,g -20 -20..... ------- * -- *20 30 40 50 60 70 0 (degrees) (b) Figure 6.3: The angular dependence of the polarimetric backscatter in terms of: (a) the scattering phase center height Z, and (b) the backscattering coefficient or~. The simulation results are fitted with polynomials of degree 3.

138 4 6 Empirical Model (m) -20 -18 -16 -14 -12 -10 -8 Empirical Model (dB) Figure 6.4: A comparison between the empirical model and the Monte Carlo coherent scattering model for a red pine stand at vv polarization configuration. (a) (b), I a 3 0 2 3 4 5 6 7 Empirical Model (m) -20 -18 -16 -14 Empirical Model (dB) Figure 6.5: A comparison between the empirical model and the Monte Carlo coherent scattering model for a red pine stand at vh polarization configuration. (a) (b) 0 ii1 E co 0 0 2 0 0 2 4 6 — 15 -10 -5 Empirical Model (m) Empirical Model (dB) Figure 6.6: A comparison between the empirical model and the Monte Carlo coherent scattering model for a red pine stand at hh polarization configuration.

139 constraint minimization problem. Suppose there exists a set of measurements M, the problem is defined as characterization of x so that the objective function (or error function), defined by 6(x)= LC(x)- Mll2. (6.8) is minimized over a pre-defined domain for x. Here, 11 denotes the norm of the argument. As mentioned earlier there are a number of inversion processes available in the literature; however, in this study by constructing a simple empirical model a traditional least-square minimization approach and a stochastic global minimization method are examined. 6.3.1 Least-Square Approach As it was shown in Section 2, the scattering problem can be cast in terms of a linear system of equations of the form fC(u) = Au where A is an m x n matrix and u is an n-dimensional vector in D D C Rn. For a given m-dimensional vector G, (6.8) can be expanded as m n ~ = ( aijuj — g)2.6.9) i=l j=l A solution that minimizes ~ must satisfy O9 =0, j=1,2,...,n (6.10) and is referred to as the least-square solution. It is shown that the solution of (6.10) (um) can be obtained from the solution of the following matrix equation [67]: (A*A) ~ um = A* G. (6.11) Here, A* is the transpose of A. It is also demonstrated that the solution Um = (A*A)-lA*G exists if rank(A) = n. This requirement states that the number of independent equations should exceed the number of unknowns.

140 To apply (6.11) to our empirical model using (6.3), it is noted that C(x) -, ~ = A. (x - x0), (6.12) thus we use the substitution u = x - xo and G = M - ~~. Here C~ and A are evaluated from (6.6) and the solution is given by Xm = X0 + (A*A)-1A*(M - ~). (6.13) The least-square solution may not be suitable for the inverse problem at hand for two reasons. First, the number of output channels m is usually less than that of unknown parameters n. In this case, rank(A) < n, and A*A is not invertible. Even when the number of channels is larger than the unknowns, the solution provided by (6.13) may not be accurate. This happens when A*A is ill-conditioned. Basically, some elements of (A*A)-1 become very large which amplify the errors in M [68]. 6.3.2 Genetic Algorithms In recent years, applications of genetic algorithms to a variety of optimization problems in electromagnetics have been successfully demonstrated [69,70]. The fundamental concept of genetic algorithms (GAs) is based on the concept of natural selection in the evolutionary process which is accomplished by genetic recombination and mutation. The algorithms are based on a number of ad hoc steps including: 1) discretization of the parameter space, 2) development of an arbitrary encoding algorithm to establish a one-to-one relationship between each code and the discrete points of the parameter space, 3) random generation of a trial set known as the initial population, 4) selection of high performance parameters according to the objective function known as natural selection, 5) mating and mutation, 6) recursion of steps 4 and 5 until a convergence is reached. Figure 6.7 shows the flow chart of GAs. Note

141 Discretization of Parameter Space and Binary Transformation Population Size Decodingf____ 1 Initialization of -Decoding of Genes Population Random Number Generator -- No Evaluation of Cost Function Convergence Check - Natural Selection es Mating KStop Decoding of Genes Mutation Figure 6.7: A flow chart of a genetic algorithm. that the population size is provided by the user and an initial population of the given size is generated randomly. GAs, as the revolutionary-type search routines, are different from the traditional gradient-based methods and have the following features: 1) GAs work with a coding of the parameters, not the parameters themselves. 2) GAs search for a population of points, not a single point. 3) GAs use information based on the objective function, not derivatives or other auxiliary knowledge. 4) GAs use probabilistic transition rules, not deterministic rules. In this study, since we have as many as six input ground truth parameters and six output channels, it is expected that the objective function is complex and highly non-linear containing many local minima. In this case, the traditional gradient

142 based optimization methods usually converge to a local minimum and fail to locate the inverted data. One interesting feature of GAs is that the method would provide a list of optimal solutions instead of a solution. This is important in a sense that a solution that best meets the physical constraints (not included in the objective function) may be selected from the list of optimal solutions. For this problem, each of the input parameters was discretized and encoded into a 4-bit binary code, creating a discrete input vector space with 224 members. A population of 240 members was used for each generation and the objective function was defined by ~(x) = lw.[M - ~ A. (x-xo)] 112, (6.14) where w is a user-defined weighting function assigned to individual output channels. To examine the performance of this GA-based inversion algorithm, many arbitrary points within the domain of the input vector space were selected and then the Monte Carlo simulation was used to evaluate the polarimetric backscattering coefficients and the scattering phase center heights at 5.3 GHz. The output of FCSM for these simulations were used as a synthetic measured data set M for the inversion algorithm. Figures 6.8(a)-6.8(f) show the performance of the inversion algorithm through comparisons of the input parameters x and the inverted parameters x'. Also shown in each of the figures is the calculated average error ir, defined by EN [ _Z -/ -NAx — ' (6.15) where N is the number of points (N=10 in this case), and Ax is the range of validity of the parameter according to the empirical model. It should be noted here that the quantization error for 4-bit quantization (~3%) is also included in the results. To examine the importance of the quantization error and the stochastic nature of the

143 solution in the inversion process, the inversion process was applied to another set of synthetic measurement data generated by the empirical model. Figures 6.9(a) — 6.9(f) show the comparison between the actual input x and the inverted solution x'. It is noticed that the error in Figures 6.9(a)-6.9(f) are slightly smaller than those obtained from Figures 6.8(a)-6.8(f). This indicates that the quantization error and the stochastic nature of the solution are considerable factors on the overall error. Increasing the quantization level to 5 bits increases the members of the input vector space by factor of 26. This slows down the inversion process since the population in each generation must also be increased. However, this does not improve the overall accuracy drastically as the errors inherent in the empirical model and those caused by the stochastic nature of the GA solution are independent of quantization error. At last, the developed inversion algorithm is tested using the real measured data acquired by the JPL TOPSAR over a test stand of red pine forest in Raco, Michigan. Although only -four data points (C-band v-polarized backscattering coefficients and scattering phase center heights at incidence angles 0 = 39~ and 53~) are available, the inversion algorithm can be easily modified via the objective function of the GA. In this case, the objective function is given by ~(x)-=^(x) + ~2(x), (6.16) where ~i(x) = |w. [Mi-r;-Ai.(x-Xo)] ||2" (6.17) S2(x) = ||w [M2-,Co-Al (x-Xo)] || 2 (6.18) Here the subscript 1 and 2 denotes, respectively, the case for the incidence angle 0 = 39~ and 53~. Note that the weighting function and the measured vectors in this

144 III r ~ -~ 10 Q. 14 a "0 -o 13 > 12 _C 4 —' I() _c 9.5 9 Tj = 3.1% 8.5F 8 7.5 -' 7.5 8 8.5 9 9.5 10 10.5 Input Ht (m) (b) 16if ", (a) 0.105 0.1 a 0.095 ~' 0.09 a) > 0.085 C 0.08 =6.2% 0 0I0 0 0 o7v.8 a.8.9.9. C 15.5 - 15 1 14.5 14 n = 7.4% 0O 0 0/,60 II 7lli v —, i i I I 0.075 0.08 0.085 0.09 0.995 0.1 Input Dt (m ) (c) 13.5 "' 13.5 14 14.5 15 15.5 Input 0b (deg) (d) 0.5. 16 to 0.5 U 0.45 'c 0.4 0.45 E 0.4 "O ' 0.35 0.3 T1 = 6.3% 0/ 00 0 0v 0.3 v 0.3 0.35..9 I. 0.4 0.45 0.5 0.55 Input ms 0.25 0.3 0.35 0.4 0.45 Input mw 0.5 (e) (f) Figure 6.8: Comparison of the input parameters (x) and the output of the inversion algorithm (x') using the synthetic data obtained from the Monte Carlo simulation, for (a) trunk diameter da, (b) tree height Ht, (c) tree density Dt, (d) branching angle Tb, (e) soil moisture mi, and (f) vegetation moisture mr,. Here r is a measure of the average error in the inversion process defined by equation(6.15).

145 15 H-4 a 14 -0 13 ) > 12 (a) (b) 16, -. p 0.095 D 0.09 0) 0.085 > 0.085 c: 15.5 CY) UD 15 'a a) ' 14.5 14 c 14 T1=6.4% X 0O /0 13.5 13 9.! I 5 14 14.5 15 15.5 Input b (deg) 16 (c) (d) 0.55 o 0.5 '0. 0.45 -E 0.4 0.35 = 5.3% O I T _l L I 0.3 0.35 0.4 0.45 0.5 0.55 Input ms 0.35 0.4 0.45 0.5 Input mw (e) (f) Figure 6.9: Comparison of the input parameters (x) and the output of the inversion algorithm (x') using the synthetic data obtained from the empirical model, for (a) trunk diameter da, (b) tree height Ht, (c) tree density Dt, (d) branching angle 0b, (e) soil moisture m,, and (f) vegetation moisture mw. Here r1 is a measure of the average error in the inversion process defined by equation (6.15).

146 Measured Z V(0 = 39~) = 3.6m a (0 = 39~)= -10.26 dB Channels M1 Zv (O = 53~) = 6.1m ]r( = 53~)= -13.07 dB Parameters da(cm) Ht(m) Dt(#/lHectare) Ob ms mw Input x 12.9 9.00 907 14.8~ 0.46 0.38 Inverted x' 12.4 9.37 945 13.9~ 0.44 0.37 Table 6.2: The inversion results using the interferometric TOPSAR data as the measured channels M. Here x is the actual ground truth data and x' is the output of the inversion process. case are written as w= 1 0 0 0 0] (6.19) M = Zvv (0=390) 0 0 ao(0=O390) 0 0], (6.20) -- t M2 Z= vv(0=53) 0 0 7~,(0=53~) 0 0] (6.21) The simulation results are shown in Table 6.2 where a very good agreement is achieved. 6.4 Conclusions In this study, a simplified empirical model was developed using a high fidelity Monte Carlo coherent scattering model to be incorporated in an efficient inversion algorithm. The empirical model was specifically developed for a red pine forest stand which provides simple expressions for the polarimetric backscattering coefficients and scattering phase center heights at C-band as a function of the incidence angle. The accuracy of the empirical model was examined by comparing its output with that of the Monte Carlo fractal-based coherent scattering model. The empirical model in conjunction with a stochastic search algorithm (genetic algorithm) were

147 used to construct an inversion algorithm. The accuracy of the inversion algorithm was demonstrated by first using synthetic measured data generated from the empirical model and the Monte Carlo FCSM. It was shown that the inversion algorithm can accurately estimate the input parameters where synthetic data were used. Next we applied the inversion algorithm to an actual data set, obtained from TOPSAR, composed of vv-polarized backscattering coefficient and scattering phase center height at C-band and at two incidence angles. Excellent agreement was obtained between the ground truth data and the output of the inversion algorithm. Using the combination of a simplified empirical model and a GA-based inversion process, the investigation of this study provides a potential application to forest monitoring.

CHAPTER VII CONCLUSIONS AND RECOMMENDATIONS 7.1 Summary of Achievements The main theme of this dissertation is the development of a realistic and compre.hensive electromagnetic scattering model for characterizing radar response to forest canopies. Currently, the most widely used scattering model for forest canopies is based on radiative transfer (RT) theory which treats the vegetation as a random medium composed of uniformly distributed scatters. However, the natural forest canopy has an apparent inhomogeneous architecture in both vertical and horizontal directions. In addition, RT is an incoherent approach that restricts its application to radar interferometry. These two deficiencies have been addressed in this study through the development of a fractal-based coherent scattering model for forest canopies. Based on fractal theory, realistic 3-D tree structures have been first generated as the simulated targets for the simulation of radar backscatter. Using L-systems, the useful tools for constructing fractal patterns, versatile generating algorithms have been developed in this study covering coniferous and deciduous trees with fine botanical features (like needle-leaflet clusters) incorporated. Two real forest stands (red pine and red maple) near Raco, Michigan have been simulated and high fidelity 148

149 has been achieved. Robust computer codes have been developed to generate and visualize the desired tree structures. A comprehensive electromagnetic scattering model has been constructed in a most general configuration. The ground plane was modeled as a half-space dielectric medium with an arbitrary tilt angle. The tree trunk was modeled as a stratified lossy dielectric cylinder with corrugated bark layer on the surface, from which the bistatic scattering has been investigated by deriving the theoretical semi-exact solution and PO approximation. The broad leaf and needle leaf were modeled as homogeneous dielectric thin disks and cylinders respectively, where the generalized Rayleigh-Gan's approximations were used to characterize their scattering properties. To model the backscattering from a forest canopy, the fractal-generated tree structure was considered as a cluster of scatterers composed of dielectric cylinders and disks. Using the single scattering approximation, the total scattered field was obtained from the coherent addition of the individual scattering from each scatterer illuminated by a mean field. Foldy's approximation has been invoked to model the coherent wave propagation within the forest canopy where the mean field at a given point accounts for the accumulated attenuation and phase change caused by vegetation components. Finally, the desired statistics of the scattered field were acquired using a Monte Carlo simulation over a large number of realizations. Combined with the recently developed Ak-radar equivalence algorithm, the developed model, for the first time, can simulate the correlation coefficient and the scattering phase center height of a forest canopy, the very quantity measured by an interferometric SAR (INSAR). The accuracy of the model has been verified using SIR-C (SAR) and TOPSAR (INSAR) measurements over test stands of coniferous and deciduous types near Raco, Michigan.

150 A thorough sensitivity analysis has been performed to characterize the response of the forest canopy with controlled physical parameters (geometry, density, and composition) to a multi-parameter (frequency, polarization, incidence angle, and baseline) SAR/INSAR. To improve the model utility and efficiency, a procedure for developing simplified empirical formulation for scattering from complex forest structures was presented in this study. Finally, an inverse scattering model based on genetic algorithms has been developed for retrieval of forest parameters that may take multi-parameter SAR/INSAR data as its input. 7.2 Recommendations for Future Work Several possible directions could be taken as future research: Database Construction Since most of the efforts of this study were put into the exploration of a new model, only two forest stands were investigated. In the future, to expand the application of the present model, more fractal generating algorithms for other tree species or stand should be developed. Associated with the ground truth inventory, a updated database is required to have accurate model input for the scattering model. Model Improvements Surface roughness and multiple scattering are two important factors to be taken into account in the future improvement work. For short tree canopies, the direct backscattering for the rough surface may not be neglected. To accurately model the cross-polarized backscatter, the multiple scattering among the scatterers may play a, crucial role. Application to Forest Management In this study, a procedure for developing an empirical model and inverse algorithms was proposed and some preliminary

151 results have achieved. However, only six parameters were selected with a limited range of validity centering the calibrated ground truth. In the future work, more parameters and wider range of validity should be added and expanded respectively to meet the real needs for the application in forest monitoring or management.

BIBLIOGRAPHY 152

153 BIBLIOGRAPHY [1] G. M. Woodwell, The Role of Terrestrial Vegetation in the Global Carbon Cycle. New York: Wiley, 1984. [2] A. M. Solomon and H. H. Shugart, Vegetation Dynamics and Global Change. New York: Chapman and Hall, 1992. [3] C. Elachi, Introduction to The Physics and Techniques of Remote Sensing. New York: Wiley, 1987. [4] C. Elachi, Spaceborne Radar Remote Sensing: Applications and Techniques. New York: IEEE Press, 1988. [5] F. T. Ulaby, R. K. More, and A. K. Fung, Microwave Remote Sensing: Active and Passive, Vol. I - Fundamental and Radiometry. Artech House, 1986. [6] F. T. Ulaby, R. K. More, and A. K. Fung, Microwave Remote Sensing: Active and Passive, Vol. II - Radar Remote Sensing and Surface Scattering and Emission Theory. Artech House, 1986. [7] J. C. Curlander, Synthetic Aperture Radar: systems and signal processing. New York: Wiley, 1991. [8] L. C. Graham, "Synthetic interferometric radar for topographic mapping," Proc. IEEE, vol. 62, pp. 763-768, 1972. [9] H. A. Zebker, T. G. Farr, R. P. Slazar, and T. Dixon, "Mapping the world's topography using radar interferometry: the topsat mission," Proc. IEEE, vol. 82, pp. 1774-1786, 1994. [10] H. A. Zebker, C. Werner, P. A. Rosen, and S. Hensley, "Accuracy of topographic maps derived from ers-1 interferometric radar," IEEE Trans. Geosci. Remote Sensing, vol. 32, pp. 823-836, 1994. [11] M. C. Dobson, L. E. Pierce, and F. T.T. Ulaby, "Knowledge-based land-cover classification using ers-1/jers-1 sar composites," IEEE Trans. Geosci. Remote Sensing, vol1. 34, pp. 83-99, 1996.

154 [12] K. S. Chen, W. P. Huang, D. H. Tsay, and F. Amar, "Classification of multifrequency polarimetric sar imagery using a dynamic learning neural network," IEEE Trans. Geosci. Remote Sensing, vol. 34, pp. 814-819, 1996. [13] U. Wegmuller and C. Werner, "Retrieval of vegetation parameters with sar interferometry," IEEE Trans. Geosci. Remote Sensing, vol. 35, pp. 18-24, 1997. [14] M. C. Dobson, F. T. Ulaby, T. L. Toan, A. Beaudoin, and E. S. Kasischke, "Dependence of radar backscatter on conifer forest biomass," IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 402-415, 1992. [15] M. C. Dobson, F. T. Ulaby, L. E. Pierce, T. L. Sharik, K. M. Bergen, J. Kellndorfer, E. L. J. R. Kendra, Y. C. Lin, A. Nashashibi, K. Sarabandi, and P. Siqueira, "Estimation of forest biophysical characteristics in Northern Michigan with SIRC/X-SAR," IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 877-895, 1995. [16] J. 0. Hagberg, L. M. H. Ulander, and J. Askne, "Repeat-pass sar interferometry over forested terrain," IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 331 — 340, 1995. [17] U. Wegmuller and C. L. Werner, "Sar interferometry signature of forest," IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 1153-1161, 1995. [18] J. A. Richards, G. Q. Sun, and D. S. Simonett, "L-band radar backscatter modeling of forest stands," IEEE Trans. Geosci. Remote Sensing, vol. 23, pp. 487-498, 1987. [19] S. L. Durden, J. J. van Zyl, and H. A. Zebker, "Modeling and observation of the radar polarization signature of forested areas," IEEE Trans. Geosci. Remote Sensing, vol. 27, pp. 290-301, 1989. [20] A. K. Fung and H. S. Fung, "Application of first order renormalization method to scattering from a vegetation-like half space," IEEE Trans. Geosci. Remote Sensing, vol. 15, pp. 189-195, 1977. [21] K. Sarabandi, Electromagnetic Scattering from Vegetation Canopies. PhD thesis, University of Michigan, 1989. [22] Y. C. Lin and K. Sarabandi, "Electromagnetic scattering model for a tree trunk above a ground plane," IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 1063 -1070, July, 1995. [23] 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. [24] M. A. Karam, A. K. Fung, and Y. M. M. Antar, "Electromagnetic wave scattering from some vegetation samples," IEEE Trans. Geosci. Remote Sensing, vol. 26, pp. 799-808, 1988.

155 [25] R. H. Lang and J. H. Sidhu, "Electromagnetic backscattering from a layer of vegetation: a discrete approach," IEEE Trans. Geosci. Remote Sensing, vol. 21, pp. 62-71, 1983. [26] 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. [27] F. T. Ulaby, K. Sarabandi, K. MacDonald, M. Whitt, and M. C. Dobson, "Michigan microwave canopy scattering model," Int. J. Remote Sensing, vol. 11, pp. 1223-1253, 1990. [28] L. Tsang, C. H. Chan, J. A. Kong, and J. Joseph, "Polarimetric signature of a canopy of dielectric cylinders based on first and second order vector radiative transfer theory," J. Electromag. Waves and Appli., vol. 6, pp. 19-51, 1992. [29] M. A. Karam, A. K. Fung, R. H. Lang, and N. H. Chauhan, "A microwave scattering model for layered vegetation," IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 767-784, 1992. [30] S. H. Yueh, J. A. Kong, J. K. Jao, R. T. Shin, and T. L. Toan, "Branching model for vegetation," IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 390 — 402, 1992. [31] J. M. Stiles and K. Sarabandi, "Scattering from cultural grass canopies: a phase coherent model," in Int. Geoscience and Remote Sensing Symposium (IGARSS) Digest, pp. 720-722, 1996. [32] C. C. Borel and R. E. McIntosh, "A backscattering model for various foliated deciduous tree types at millimeter wavelengths," in Int. Geoscience and Remote Sensing Symposium (IGARSS) Digest, pp. 867-872, 1986. [33] B. Mandelbro, Fractal Geometry of Nature. W. H. Freeman Company, 1983. [34] G. Zhang, L. Tsang, and Z. Chen, "Collective scattering effects of trees generated by stochastic Lindenmayer systems," Microwave and Opt. Tech. Lett., vol. 11, pp. 107-111, 1996. [35] P. Prusinkiewicz and A. Lindenmayer, The Algorithmic Beauty of Plants. Spring-Verla, 1990. [36] R. H. Lang, R. Landry, O. Kavakhioglu, and J. C. Deguise, "Simulation of microwave backscatter from a red pine stand," in Multispectral and Microwave Sensing of Forestry, Hydrology, and Natural Resources, SPIE, vol. 2314, (Rome, Italy), pp. 538-548, 1994. [37] K. Sarabandi, "Ak-radar equivalent of interferometric sars: a theoretical study for determination of vegetation height," IEEE Trans. Geosci. Remote Sensing, vol. 35, pp. 1267-1276, 1997.

156 [38] R. N. Treuhaft, S. N. Madsen, M. Moghaddam, and J. J. van Zyl, "Vegetation characteristics and underlying topography from interferometric radar," Radio Science, vol. 31, pp. 1449-1485, 1996. [39] K. Sarabandi and Y. C. Lin, "Simulation of interferometric SAR response for characterizing the scattering phase center statistics of forest canopies," IEEE Trans. Geosci. Remote Sensing, submitted March 1997. [40] A. Lindenmayer, "Developmental algorithms for multicellular organisms: a survey of L-systems," Journal of Theoretical Biology, vol. 54, pp. 3-22, 1975. [41] M. H. Zimmermann and C. L. Brown, Tree Structures and Function. SpringerVerlag, 1971. [42] G. T. Ruck, D. E. Barrick, W. D. Staurt, and C. K. Krichbaum, Radar Cross Section Handbook. Plenum Press, 1970. [43] S. S. Seker and A. Schneider, "Electromagnetic scattering from a dielectric cylinder of finite length," IEEE Trans. Antennas Propag., vol. 36, pp. 303-307, 1988. [44] M. 0. Kolawole, "Scattering from dielectric cylinders having radially layered permittivity," J. Electromag. Waves and Appli., vol. 6, pp. 235-259, 1992. [45] F. T. Ulaby and C. Elachi, Radar Polarimetry for Geoscience Applications. Artech House, 1990. [46] K. Sarabandi and F. T. Ulaby, 'High frequency scattering from corrugated stratified cylinders," IEEE Trans. Antennas Propag., vol. 39, pp. 512-520, 1991. [47] K. Sarabandi, "Simulation of a periodic dielectric corrugation with an equivalent anisotropic layer," International Journal of Infrared and Millimeter Waves, vol. 11, pp. 1303-1321, 1990. [48] J. J. van Zyl, "The effect of topography on radar scattering from vegetated areas," IEEE Trans. Geosci. Remote Sensing, vol. 31, pp. 153-160, 1993. [49] M. A. Morgan, "Electromagnetic scattering by stratified inhomogeneous anisotropic media," IEEE Trans. Antennas Propag., vol. 35, pp. 191-197, 1987. [50] L. Tsang, J. A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing. Wiley Interscience, 1985. [51] M. L. Imhoff, "A theoretical analysis of the effect of forest structure on SAR backscatter and the remote sensing of biomass," IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 341-352, 1995. [52] K. Sarabandi, "Derivation of phase statistics from the mueller matrix," Radio Science, vol. 27, pp. 553-560, 1992.

157 [53] R. L. Jordan, B. L. Huneycutt, and M. Werner, "The SIR-C/X-SAR synthetic aperture radar system," IEEE Trans. Geosci. Remote Sensing, vol. 33, pp. 829 -839, 1996. [54] K. M. Bergen, M. C. Dobson, T. L. Sharik, and I. Brodie, "Structure, composition, and above-ground biomass of SIR-C/X-SAR and ERS-1 forest test stands 1991-1994, Raco Michigan Site," Tech. Rep. 026511-7-T, The University of Michigan, Radiation Laboratory, 1995. [55] K. M. Bergen, M. C. Dobson, L. E. Pierce, J. Kellndorfer, and P. Siqueira, "October 1994 SIR-C/X-SAR mission: Ancillary data report, Raco Michigan Site," Tech. Rep. 026511-6-T, The University of Michigan, Radiation Laboratory, 1995. [56] 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. [57] F. T. Ulaby and M. A. Elrayes, "Microwave dielectric spectrum of vegetation, Part II: Dual-dispersion model," IEEE Trans. Geosci. Remote Sensing, vol. 25, pp. 50-557, 1987. [58] E. Rodriguez and J. Martin, "Theory and design of interferometric synthetic aperture radars," IEE Proceedings, vol. F139, pp. 147-159, 1992. [59] H. A. Zebker, S. N. Madsen, J. Martin, K. B. Wheeler, T. Miller, Y. Lou, G. Alberti, S. Vetrella, and A. Cucci, "The topsar interferometric radar topographic mapping instrument," IEEE Trans. Geosci. Remote Sensing, vol. 30, pp. 933-940, 1992. [60] Y. C. Lin and K. Sarabandi, "Monte Carlo coherent scattering model for forest canopies using fractal generated trees," IEEE Trans. Geosci. Remote Sensing, revised, September 1997. [61] V. Dimri, Deconvolution and Inverse Theory: application to geophysical problems. Elsevier, 1992. [62] M. K. Sen, Global Optimization Methods in Geophysical Inversion. Elsevier, 1995. [63] L. Pierce, K. Sarabandi, and F. Ulaby, "Application of an artificial neural network in canopy scattering inversion," Int. J. Remote Sensing, vol. 15, pp. 3263 -3270, 1994. [64] P. F. Polatin, K. Sarabandi, and F. T. Ulaby, "An iterative inversion algorithm with application to the polarimetric radar response of vegetation canopies," IEEE Trans. Geosci. Remote Sensing, vol. 32, pp. 62-71, 1994.

158 [65] F. Amar, M. S. Dawson, and A. K. Fung, "Inversion of the relevant forest and vegetation parameters using neural networks," in Proc. of Progress in Electromagnetic Research Symposium (PIERS), 1993. [66] D. E. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning. Addison-Wesley, 1989. [67] S. H. Friedberg, A. J. Insel, and L. E. Spence, Linear Algebra. Prentice Hall, Inc., 1979. [68] A. Ishimaru, Wave propagation and scattering in random media. Academic Press, 1978. [69] K. Sarabandi and E. S. Li, "Characterization of optimum polarization for multiple target discrimination using genetic algorithms," IEEE Trans. Antennas Propag., p. accepted for publication, 1997. [70] R. L.. Haupt, "An introduction to genetic algorithms for electromagnetics, IEEE Ant. and Propagat. Magazine, vol. 37, pp. 7-15, 1995.