1'04)()())43-!4 HF-UHF Propagation Prediction Over Rough Terrain Submitted to: Army Research Office Attn:Dr. James F. Harvey P.O. Box 12211 Research Triangle Park, NC 27709-2211 Prepared by: Kamal Sarabandi Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, MI 48109-2122 Phone: (313) 936-1575 Fax: (313) 647-2106 Email: saraband@eecs.umich.edu July 26, 1999 37227-1-T = RL-2512

CONTENTS Contents 1 Introduction 1 2 SUMMARY OF ACCOMPLISHMENTS 4 2.1 Small Dipole Over a Random Rough Dielectric Surface........................ 4 2.2 Numerical Approaches................................................... 5 2.2.1 Wavelet-Based Method of Moments Approach..................... 5 2.2.2 Physical Optics Iterative Approach.............................6 3 Students Supported 7 4 Publications 7 Appendix A Appendix B Appendix C HF-UHF Propagation Prediction Over Rough Terrain

1 HF-UHF Propagation Prediction Over Rough Terrain Abstract In this final report a summary of our activities with regard to HF-UHF propagation prediction over rough terrain is provided. Our research activities in the HF-UHF Propagation Prediction can be categorized into three categories: (1) development of a theoretical method to predict the field propagation from a small dipole over a half-space dielectric, (2) Numerical Wavelet-based Method of Moments approach to predict the scattering from random rough surfaces, and (3) Numerical Iterative Physical Optics technique to predict the scattering from random rough surfaces that have low to moderate rms slopes. The accurate prediction of radio wave coverage at HF-UHF frequencies over irregular terrain features is of importance in the design and development of low-cost, low power communications system. The terrain effects consist of multipath, diffraction, scattering and depolarization of the electromagnetic wave. Current methods of prediction are overly simplistic and tend to neglect phenomena which have a significant effect on radio wave propagation modeling. We investigated a number of techniques in order to predict accurately the propagation over rough terrain in the HF-UHF range. 1 Introduction With the rapid expansion of technology for mobile and wireless systems an accurate method for prediction of radio wave propagation in various environments has become essential in the design and development of efficient, low-cost, low-power communications systems which can operate in these rapidly varying environments. A problem of particular interest is the prediction of High Frequency, Ultra High Frequency (HF, UHF) radio wave propagation over irregular terrain, where there may be no line of sight (LOS) path and the received signal may consist of components of multipath, diffraction, scattering, and depolarization effects from various, and usually numerous terrain features. Realistically all terrain effects cannot be accounted for, either statistically or deterministically, with the current level of computer technology, however assumptions can be made which produce a computational model of realistic size while maintaining an acceptable degree of accuracy. Several methods are currently used in propagation prediction over irregular terrain. These consist of heuristic models, based on measurements, and electromagnetic models. The irregular terrain models can be broken down into two types, those that predict scattering from small-scale roughness, i.e., surfaces whose irregularities are small compared to electrical wavelength, and those that predict scattering and diffraction from large scale irregularities such as mountains or hills. The heuristic models are limited to very specific physical conditions at the time the measurement was made and the measurement system attributes such HF-UHF Propagation Prediction Over Rough Terrain

as frequency, bandwidth. and polarization. Moreover the cost and time to perforn these measurements is prohibitive. Current electromagnetic models assume the propagation between transmitter and receiver is limited to a very narrow Fresnel zone and the probleml can be reduced to a 2-D problem, in a plane between transmitter and receiver, of propagation over intervening obstacles. Effects of surface roughness are handled in a strictly empirical manner. One of the more widely used models, applied to large scale irregularities, assumes all obstacles consist of knife edges, essentially a screen between transmitter and receiver and Physical Optics (PO) diffraction methods are applied. [4] This method is referred to in the literature as the knife edge diffraction model. While this method is simple to implement and computationally fast, PO methods are accurate only at high frequencies, in the very far zone of the obstacle, and are invalid in or near the shadow boundary and near the surface of the obstacle which obviously does not reflect many realistic scenarios. Other high-frequency techniques have also been applied such as Geometrical Theory of Diffraction (GTD) for wedge diffraction, but these techniques tend to be cumbersome for many obstacles. Due to the limitations of both heuristic and PO models more rigorous electromagnetic models have been investigated in recent year. Several methods have been proposed for modeling of large terrain features, including parabolic equation techniques which assume cylindrical wave scattering [5], and integral equation (IE) techniques [6] [7]. IE techniques have been an area of major research in electromagnetics in recent years. An integral equation is formulated by enforcing boundary conditions on the surface of the scattering object. The unknown quantities in the IE are the surface currents. These are solved for either by numerical or iterative techniques. This method accounts for all electromagnetic phenomena and interactions and is accurate to the degree of accuracy of the solution technique. When a numerical solution is sought, this technique is limited by the electrical size of the problem. While advances in computer technology have made the solution of larger problems possible when using this technique, the size of radio wave propagation problems is still prohibitive if all electromagnetic interactions are to be accounted for. With this in mind, techniques are sought using IE methods to reduce the size of the radio wave propagation problem to a manageable size while maintaining an acceptable degree of accuracy. Two techniques to achieve this for large scale irregularities are proposed, a numerical technique using a MoM approach with wavelet basis functions, and an Iterative PO approach. In addition a theoretical technique is proposed, using a perturbation method, to predict scattering and diffraction from a surface with small scale irregularities when both transmitter and receiver are near the surface. A widely used numerical technique to solve IE equations is the Method of Moments (MoM) technique. This technique expands the unknown surface currents, the expansion consisting of basis functions with unknown coefficients. The expansion can be over the entire domain of the scattering object or discrete sub-domains, where the scatterer is subdivided. These sub-domains are at most 1/10A in size to produced the required degree of accuracy, where A is the wavelength. The most widely used is the sub-domain method. Weighting functions are applied to reduce the average residual error. Direct application of MoM produces matrices of order N2 and number of computations of order N3, where N is the number of unknowns in the matrix formulation. Both the order of the matrix size and number of computations are unacceptable for a problem the size of the radio wave propagation problem. HF-UHF Propagation Prediction Over Rough Terrain

Current research is directed towards reducing the problem size while maintaining accuracy using both acceptable approximations and computational techniques. Several methods have been proposed to achieve the goal of reducing the MoM problemi size in the area of radio wave propagation over hilly terrain. As stated previously most mnethods assume a 2-D problem, calculating propagation loss along a 2-D path, and all methods cited assume this unless otherwise noted. Hviid et.al [6] proposes a method which is valid for vertical polarization only. In this method the IE is formulated using the Magnetic Field Integral Equation (MFIE) and solves for magnetic currents on the surface of the scatterer. It assumes a smooth, perfectly reflecting smooth surface and to simplify the problem does not account for backscatter. Brennan et.al, [7] extends this method by applying the Fast Far Field Approximation (FAFFA), proposed by Chew et.al. [8] to reduce the computational requirements of the problem. In this technique local groups are defined in the problem where the near field interaction between sub-domains is significant and the IE is fully enforced. It is then assumed that the distance between groups is far enough to assume plane wave interaction between them and the full IE need not be enforced thus reducing the computational size of the problem. Both methods described have limitations, including smooth surface assumptions and no backscatter in the first, and in the second the need to define near and far zone. A method is sought to more fully account for the electromagnetic interactions while maintaining accuracy and a computationally efficient problem. With this in mind two techniques are proposed to solve the IE, a numerical method using wavelet basis functions to reduce the problem size while maintaining the accuracy requirement and an iterative PO technique, applied to the magnetic field integral equation, which also produces accurate results while significantly reducing run time. These techniques are valid when the fine details of the terrain are not significant in the model as is the case at frequencies up to UHF. The MoM techniques described previously use standard basis functions such as pulse basis functions and point matching weighting which essentially enforces the MoM formulation at the center of each sub-domain. Alternative basis functions have been investigated in recent years which produce acceptable accuracy while significantly reducing the problem size and runtime. Wavelets are an adaptive basis function which has seen wide application in recent years in both communications and signal processing and electromagnetic modeling applications. In recent years it has been shown that the orthogonal properties of the wavelet basis functions in addition to having vanishing moments produce a sparse MoM matrix which can be solved using iterative solver techniques such as Conjugate Gradient (CG), thus creating a significant reduction in problem size and a significant speedup in the run time of the problem. [9, 10, 11, 12] This method produces accurate results without the need to define near and far-field domains and produces a full bistatic pattern, while maintaining the accountability for all electromagnetic interactions inherent in MoM techniques. This method has recently been applied to the problem of scattering from rough surfaces very successfully. In the iterative PO technique the MFIE may be enforced repeatedly to obtain the solution, with the previous solution substituted for each successive iteration. The number of iterations is dependent on the accuracy desired. A significant advantage to the iterative PO solution is the memory requirements are of order N as opposed to order N2 for standard MoM. HF-UHF Propagation Prediction Over Rough Terrain

4 2 SUMMARY OF ACCOMPLISHMENTS Because no single method is comprehensive in its results (quick computations, exact, accurate over all ranges of physical parameters, etc.) our techniques have focused on providing a number of desired qualities including accuracy of results and computational speed. To this end we have focused on three different techniques: (1) Theoretical results of scattering iII a dielectric medium, (2) semi-exact scattering results from random rough surfaces using a wavelet-based technique, and (3) fast scattering results using an iterative Physical Optics approach for random rough surfaces with low to moderate rms slopes. A summary of the work accomplished in each area is given next. 2.1 Small Dipole Over a Random Rough Dielectric Surface With recent progress in wireless technology and ever increasing demand for reliable low power wireless systems, the need for predicting system performance has become an extremely important step in the design of such systems. At HF through UHF, the channel characteristics such as attenuation, and multipath fading statistics significantly affect the performance of wireless systems. Prediction of the channel characteristics can be accomplished using physics based propagation models. For this purpose precise diffraction models are developed and incorporated with accurate terrain data base to construct a realistic propagation model. In this study the problem of electromagnetic wave propagation, excited by a short dipole, above a dielectric ground plane with an arbitrary dielectric profile and an irregular interface is studied. This investigation is a natural extension of the classical Sommerfeld problem with the exception of the random surface irregularities at the interface between the two dielectric media. Assuming that the interface profile height variations are small compared to the wavelength, the problem is formulated as follows. First, the bistatic scattering of a plane wave illuminating the rough surface is solved using a perturbation solution of an integral equation for the induced polarization current. Analytical expressions for the coherent field (mean-field) and incoherent scattered power at an arbitrary observation point, including points near the interface, are obtained. Then the solutions for the mean-field and incoherent scattered power generated by a small dipole of arbitrary orientation and position are derived by expanding the field of the dipole in terms of a continuous spectrum of plane waves and using superposition. The effect of rough interface on the surface waves and the phenomenon of depolarization caused by the rough interface are studied. In reality, the interface between forest canopies and air is not flat, hence it is not clear whether the lateral wave can be excited and if it can how the surface roughness affects it. In this study, the effect of roughness of interface between canopy and air on the wave propagation in forested area is investigated. Also, an expression for the mean-field of an infinitesimal dipole of arbitrary orientation is derived by obtaining a partial second order solution of the Born approximation and a sensitivity analysis is carried out to demonstrate the variations of the mean-field to physical parameters such as a effective permittivity, location of the dipole and observation points, and surface roughness. More details may be found in Appendix A. HF-UHF Propagation Prediction Over Rough Terrain

2.2 Numerical Approaches 25 2.2 Numerical Approaches 2.2.1 Wavelet-Based Method of Moments Approach The problem of electromagnetic scattering from rough surfaces has been the subject of intensive investigation over the past several decades for its application in a number of important remote sensing problems. Radar remote sensing of the oceans, soil moisture, and mine detection using wideband radars are such examples. For these problems, where the rough surface is either the primary target or the clutter, the understanding of interaction of electromagnetic waves with the rough surface is essential for developing inversion or detection algorithms. An exact analytical solution for random rough surfaces does not exist. An alternative approach for evaluating of the scattered field and its statistics for rough surfaces is Monte Carlo simulation. In order to use Monte Monte Carlo simulation for evaluating the scattering statistics of rough surfaces more routinely, computationally more efficient scattering codes must be developed. In this paper the application of wavelets as a basis function for the expansion of induced surface currentss is considered. Traditional method of moments (MoM) in conjunction with Galerkin's method would require matrix fill computation time of the order of N2 and matrix inversion computation time of the order of N3 (using Gaussian elimination). It is well known that the solution of linear system of equations can be obtained far more efficiently using search routines, such as the Conjugate Gradient method, if the matrix of the coefficients is a sparse matrix. In MoM, the application of conventional pulse or rooftop basis and testing functions would usually produce full impedance matrices. Although the diagonal elements are usually larger than the rest of the elements, the smaller elements cannot be arbitrarily thresholded without drastically altering the resulting scattering pattern. The success of wavelet expansion function in generating sparse matrices have been demonstrated for many circuits and antenna problems [9, 10, 11, 12]. In the Monte Carlo simulation of scattering from rough surfaces the quantities of interest are the statistical parameters, such as the mean and variance of the scattered field, and therefore it is expected that the overall accuracy be less sensitive to the threshold level. A comparative study is carried out to demonstrate the application of wavelets for improving the computation time and reducing computational memory required for evaluating the statistics of the scattered field from rough surfaces using the method of moments in conjunction with a Monte Carlo simulation. In specific, Haar and the first-order B-spline wavelet basis functions are applied to the MoM formulation of two-dimensional rough surfaces in order to compare the computation time and sparsity for wavelets in the same family but of higher order. Since the scattering coefficient (the second moment of the backscatter field per unit area) is a gentle function of the surface parameters and the radar attributes, it is demonstrated that a relatively high thresholding level can be applied to the impedance matrix which leads to a sparser impedance matrix and faster computation time. It is also shown that applying a high threshold level the coefficients of the high order wavelets would increase out of proportion, however the effect of these current components averages out when computing the scattering coefficients. The resulting sparse impedance matrices are solved efficiently using fast search routines such as the conjugate gradient method. A systematic study is carried out to investigate the HF-UHF Propagation Prediction Over Rough Terrain

2.2 -Numerical Approaches 6 effect of different threshold levels on the accuracy versus computing speed criterion. The computed scattering coefficients are compared to previous results comnuted using a Coilven tional pulse basis function as well as the existing theoretical solutions for rough surfaces. It is shown that wavelet basis functions provide substantial reductions in botilh lemory requirements and computation time. This procedure is further discussed in Appendix B. 2.2.2 Physical Optics Iterative Approach Development of numerically efficient Monte Carlo-based models for simulations of electromagnetic scattering from random surfaces has attained significant prominence over the past decade. A major stumbling block in this endeavor has been the large memory and computation time requirement. This is because the size of the scatter is large compared to the wavelength and the Monte Carlo simulations demand computation of the scattering problem many times. Iterative methods offer an alternative approach when exact solutions are not available and have been used in different electromagnetic problems. By construct evaluation of iterative solutions are rather straight forward especially when the perturbation parameter is relatively small. Physical Optics (PO) approximation is known to provide accurate approximation for the induced surface currents provided that the local radii of curvature at each point on the surface of scatterer is large and the surface is convex. For concave surfaces and surfaces with many adjacent humps multiple scattering drastically alter the standard PO current. However these surface current variations can be estimated through an iterative process. Unlike Method of Moments (MoM) which requires matrices on the order of N2 to find the surface current of the sample surface with N elements, the iterative Physical Optics method only requires memory size of the order N. Thus, substantial memory savings are realized. Also, since no solver routine is necessary in order to solve for the surface currents, as in the MoM, substantial time savings are realized as well. Iterative methods, such as iterative PO offer an alternative approach when exact solutions are not available. Evaluation of iterative solutions are rather straight forward especially when the perturbation parameter is relatively small. PO approximation is known to provide accurate approximation for the induced surface currents provided that the local radii of curvature at each point on the surface of the scatter is large and the surface is convex. For concave surfaces and irregular surfaces with many adjacent obstacles multiple scattering drastically alter the standard PO current. These surface current variations can be estimated through an iterative process producing a problem size of order N as previously stated. The computation time of this problem then becomes of order N which is significantly less than the MoM technique which is typically of order N3. The application of iterative Physical Optics (PO) in conjunction with a Monte Carlo simulation for characterizing the bistatic scattering coefficient of random rough surfaces is examined. The iterative PO method offers decreased memory and computation time restrictions compared to the standard numerical methods such as the Method of Moments (MoM). Results from the iterative PO method are compared to the standard electric field integral equation (EFIE), the magnetic field integral equation (MFIE) as well as the existing theoretical solutions for rough surfaces. It is demonstrated that memory requirements and computation time is significantly decreased, even compared to traditional MoM techniques HF-UHF Propagation Prediction Over Rough Terrain

-7 while providing fairly accurate results for surfaces with moderate to low rnis slope. In addition, this technique is extended to the 3-D scattering problem, whliclh is prohibitively large for lMoM techniques on all but the largest and most powerful computers. Good agreement is found to occur between analytical methods and the 3-D iterative Physical Optics approach. Also. the iterative Physical Optics approach has been further simplified by using a large argument approximation for the kernel. For these results, even more time reduction is observed. We investigated the use of the iterative Physical Optics method upon a variety of surfaces in order to find the approximate region of validity for such a method. The results obtained using the iterative PO method are also compared to the results found using the electric field integral equation (EFIE) with tapered resistive sheets at the ends of the surface samples as well as the magnetic field integral equation (MFIE) for a horizontally polarized wave. This procedure is detailed in Appendix C. 3 Students Supported Throughout the course of this investigation the effort of the following Ph.D. students were supported fully or partially by this project. 1. Tsenchieh Chiu 2. Il-Suek Koh 3. Daniel Zahn 4 Publications Journal Publications: 1. Zahn, D., K. Sarabandi, and K.F. Sabet, "Numerical simulation of scattering from rough surfaces: A wavelet-based approach," IEEE Trans. Antennas Propagat., submitted for publication (Jan. 1998). 2. Chiu, T.C., and K. Sarabandi, "Field of a short dipole above a dielectric half-space with rough interface," IEEE Trans. Antennas Pr opagat., submitted for publication (August 1998). 3. K. Sarabandi, and I.S. Koh, "Effect of canopy-air interface roughness on HF-UHF propagation in forest," IEEE Trans. Antennas Propagat., submitted for publication (April 1999). Conference Publications: 1. K. Sarabandi, and D. Zahn, "Numerical Simulation of Scattering from Rough Surfaces Using an Iterative Physical Optics Approach," Proc. IEEE Trans. Antennas Propagat.& URSI Symp., Atlanta, GA, July 1999. HF-UHF Propagation Prediction Over Rough Terrain

8 2. D. Zahn and K. Sarabandi. "Numerical Simulation of Scattering from tRougll Surfaces Using a Fast Far-Field Iterative Physical Optics Approach." Proc. IEEE Trans. Antennas Propagat.& URSI Syrnp., Orlando. FL, July 1999. 3. D. Zahn, K.F. Sabet, and K. Sarabandi, 'Numerical simulation of scattering from rough surfaces: A wavelet-based approach," Proc. IEEE Trans. Antennas Propagat.&. URSI Symp., Montreal. Canada, July 1997. 4. Chiu, T.C., and K. Sarabandi, "Electromagnetic scattering from targets above rough surfaces,"Proc. IEEE Trans. Antennas Propagat.& URSI Symp., Montreal, Canada, July 1997. 5. K. Sarabandi and I. Koh, "Effect of Canopy-Air Interface Roughness on HF-UHF Wave Propagation in Forest,"Proc. IEEE Trans. Antennas Propagat.& URSI Symp., Orlando, FL, July 1999. 6. K. Sarabandi, and I. Koh, "HF-VHF Wave Propagation in Forest: Effect of CanopyAir Interface Roughness on Lateral Waves," XXVI-th General Assembly of Int. Union of Radio Sci. (URSI), Toronto, Canada, August 1999. 7. T.C. Chiu, and K. Sarabandi "Propagation of Electromagnetic Waves Near the Surface of a Half-Space Dielectric with Rough Interface," Proc. IEEE Trans. Antennas Propagat.& URSI Symp., Orlando, Florida, July 1999. HF-UHF Propagation Prediction Over Rough Terrain

REFE RENCES 9 References [1] Sommerfeld, A., Ann. Physik. Vol. 28, pp. 665. 1909. [2] Sarabandi,K., and T.C. Chiu "Electromagnetic scattering from slightly rough surfaces with inhomogeneous dielectric profile, IEEE Trans. Antennas Propagat.. vol. 45. no. 9, Sept. 1997. [3] Zahn, D., K. Sarabandi, and K.F. Sabet, "Numerical simulation of scattering from rough surfaces: A wavelet-based approach,"IEEE Trans. Antennas Propagat., submitted for publication (Jan. 1998). [4] Epstein, J., and D.W. Peterson, "An experimental study of wave propagation at 850 MC", Proc. Inst. Radio Eng., vol. 41, pp. 595-611, Jan. 1953. [5] Craig, K.H., and M.F. Levy, "Parabolic equation modeling of the effects of multipath and ducting on radar systems", IE Proc. F, vol. 138, pp. 153-162, April 1991. [6] Hviid, J.T., J.B. Andersen, J. Toftgard, and J. Bojer, "Terrain-Based propagation model for rural area-An integral equation approach", IEEE Trans. Antennas Propagat., vol. 43, pp. 41-46, Jan. 1995. [7] Brennan, C., and P.J. Cullen, "Tabulated Interaction for UHF terrain propagation problems", IEEE Trans. Antennas Propagat., vol. 46, pp. 738-739, May 1998. [8] Lu, C.C., and W.C. Chew, "Fast far field approximation for calculating the RCS of large objects", Microwave Opt. Technol. Lett., vol. 8, pp. 238-241, May 1995. [9] B.Z. Steinberg and Y. Leviatan, "On the use of wavelet expansions in the method of moments" IEEE Trans. Antennas Propagat.,vol. 41, pp. 610-619, May 1993. [10] K.F. Sabet and P.B. Katehi, "Analysis of integrated millimeter-wave and submillimeterwave waveguides using orthonormal wavelet expansions," IEEE Trans. Microwave Theory Tech., vol. MTT-42, Dec. 1994. [11] K.F. Sabet, Novel efficient integral-based techniques for characterization of planar microwave structures, Ph.D. dissertation, The University of Michigan, Apr. 1995. [12] K.F. Sabet and P.B. Katehi, "A study of dielectric resonators based on two-dimensional fast wavelet algorithm," IEEE Microwave Guided Wave Lett., vol. 6, pp. 19-21, Jan. 1996. [13] Rioul, 0. and P. Duhamel, "Fast Algorithms for Discrete and Continuous Wavelet Transforms", IEEE Trans. on Info. Theory, Vol. 38, No. 2, March 1992. HF-UHF Propagation Prediction Over Rough Terrain

APPENDIX A Small Dipole over a Random Rough Dielectric Surface

Field of a Short Dipole Above a Dielectric Iialf-SpIace with Rough Interface '.sen(cliieh Cliiu aind Kamal Sarabarldi Radiation Iaboratory Department of Electrical Engineering and Computer Science The University of Michigan, Ann Arbor, MI 48109-2122 Tel:(313) 936-1575 Abstract - The accurate prediction of radio wave coverage at lF-UHF frequencies over irregular terrain features is of importance in the design and development of low-cost and low power communications systems. In this paper the problem of electromagnetic wave propagation, excited by a short dipole, above a dielectric ground plane with an arbitrary dielectric profile and an irregular interface is studied. This investigation is a natural extension of the classical Sommerfeld problem with the exception of the random surface irregularities at the interface between the two dielectric media. Assuming that the interface profile height variations are small compared to the wavelength, the problem is formulated as follows. First, the bistatic scattering of a plane wave illuminating the rough surface is solved using a perturbation solution of an integral equation for the induced polarization current. Analytical expressions for the coherent field (mean-field) and incoherent scattered power at an arbitrary observation point, including points near the interface, are obtained. Then the solutions for the mean-field and incoherent scattered power generated by a small dipole of arbitrary orientation and position are derived by expanding the field of the dipole in terms of a continuous spectrum of plane waves and using superposition. The effect of rough interface on the surface.waves-and the phenomenon of depolarization caused by the rough interface are studied.

1 Introduction With the rapid expansion of technology for mobile and wireless systems, an accurate method for prediction of radio wave propagation has become essential in the design and development of efficient, low-cost, low-power communications systems. In many communication scenarios where both the transmitter and receiver are near the ground, shadowing and multipath significantly affect the signal strength and the coherent bandwidth at the receiver. This is specifically the case for the propagation over irregular terrain. Terrain irregularity so far as propagation is concerned can be categorized into two groups: 1) large-scale roughness, and 2) small-scale roughness. Large-scale terrain irregularities are generally referred to terrain irregularities large compared to the wavelength such as mountains and hills. Small-scale terrain irregularities, on the other hand, refers to surface roughnesses where the rms height and slope are small compared to the wavelength (at HF-UHF). These affect the wave propagation differently; for example, while large-scale terrain irregularities are the sources of shadowing and multipath, small-scale irregularities reduce the ground reflectivity and produces an incoherent field component due to surfaces scattering. Small-scale irregularities also affect the surface waves which are essential when both the transmit and receive antennas are close to the air-ground interface. Determination of field of a small dipole over a half-space dielectric is a classical problem with a well-known solution [1]. It is shown that when both the transmitter and receiver are near the surface, the contribution from surface waves is dominant. In practice, the transmit and receive antennas of mobile units are usually very close (relative to the wavelength) to the ground. Existence of surface roughness may alter the contribution of surface waves drastically. In this case the azimuthal symmetry of the problem may no longer be exploited, and the Sommerfeld solution must be modified significantly. The surface roughness generates incoherent scattered field which is the source of depolarization. In this paper, the effect of slightly rough surfaces on the radiation of a short dipole is studied. In what follows, first, a solution for the scattered field (including the near field) from a slightly rough surface illuminated by plane waves is formulated in Sections 2 and 3. To investigate the effect of small-scale surface roughness on surface waves, ground reflectivity, and the significance of the incoherent scatter fields, an analytical solution based on a perturbation theory is proposed. In this formulation the perturbation theory is applied to a volumetric integral equation for the induced polarization current in the top rough layer of the dielectric interface. The perturbation parameter is the normalized rms height of the rough surface and an iterative solution starting from the unperturbed problem (dielectric half-space with smooth interface) is obtained. Basically, the formulation is similar to what has recently been applied to evaluate far-field scattering from rough surfaces with inhomogeneous profile when illuminated by a plane wave [2]. In Section 4 the solution for dipole excitation is obtained by expanding the radiated

f1(ld I of il' (ipolc i, (I I I I1 I 1,(lO il i IS I,1S) tilll ol'| pIil ire W aves ridLl' l ridi g I ie s I,('!ll', io (M or ' l pl, l W\',a\; (),,l,'','-(i v.. lSt( i Is( af i I u-I I is i; (.(,l\ i r d, i(' I o I i,i\'f Il fi;,ll for,lI,1, lfl. i ri reI. (;1olI (rI( ll (,1 I III ) I rid,r olf ()I it (I i,lll l I io(, ) l fields. 'I e (.sll lt s,rar( (*oiip ared wit I tI e Sor.i i,erfJeid soluilon af 11 lh (cp)ollari/, ioll 1'feels ar(' ir v sligal.edl. 2 Polarization Current in a Slightly Rough Surface As mentioned earlier the first step towards evaluating the field generated by an arbitrary dipole above a ground plane with rough interface is to consider the plane wave illumination. To obtain tile scattered field, a perturbation solution to a volumretric integral equation for thle induced polarization current over tihe top rougli layer of HtIe surface is derived using a procedure similar to wliat is presented in [2]. Figure 1 shows the geometry of tlie scattering prob)lemn wliere a dielectric lhalf-space with an arbitrary dielectric profile and rough interface is illuminated by a plane wave from the upper medium. Suppose the surface height variation is small compared to the wavelength (A) of the incident wave. The incident wave with an arbitrary polarization Q can be written as Ei(r) = Q etkok'r where ko = 2 is the free space propagation constant, and kI is the unit vector along the direction of propagation, given by kt = sin Oi cos Oixx + sin Oi sin iY- - cos Oiz = k' - zk' To make the solution tractable, the permittivity of the top layer down to a depth of d is considered to be uniform, where -d < min{surface profile}. Denote the surface height profile by function z = Af (x, y), where f (x, y) is a zero-mean stationary random process with a known autocorrelation function and variance 1, and A << A is a small constant known as the perturbation parameter. In the following derivation, it is assumed that the medium below the top layer is stratified, that is, the relative permittivity is only a function of z. In the absence of the top homogeneous rough layer, the incident wave would be reflected at the smooth interface between the free space and the stratified half space soil medium. This reflected wave can be expressed by Er(r) = Er(O) eikokr- r where kr is the direction of propagation of the reflected wave, given by k = - -2(i. - ')z = k' + kZ andl E-(O) denllotes the magnitude and )olarization vector of tlie reflected %wave, wliic can be obtained from Er() = [^.vr +r ^ I.] Q

IHre r,., and ro, are the nresnel reflection coefficients, and the horizontal and vertical unit vectors are given by = k x i |k.1 X Z I s, = h, x 5s (1) where the subscript s can be i or r for the incident and reflected waves, respectively. In presence of the homogeneous rough layer, the incident and reflected waves induce a polarization current within the top dielectric layer which is the source of the scattered field. The polarization current in terms of the total field and the permittivity of the layer is J(r) = -ikoYo( - 1)Et, (2) where % o= - is the characteristic admittance of the free space, and E' E = + Er + ES The scattered field E can in turn be expressed in terms of the polarization current and is given by E = ikoZo G(r,r') J(r')dv', Vslab (3) where G(r, r') is the dyadic Green's function of the half-space stratified medium (in the absence of the top rough layer), and is given by [3] G(r, r') =- - k ) d d2k 8i 2, ek(p-p ') 87 2 d kz {[rhh(kz)ekzz + h(-,kz)e-ikz] h(-k,) + [rvv-(kz)eiz + vi(-kz)e-ikz] V(-k) }eikzz' {I(kz) [rhh(-kz)eikz' + h(kz)e-i'z'] +v(kz) [rvv(-kz)eikzz' + v(kz)e'-k'] }eikz', if z < z' if z > z' (4) In (4), k, = /k2- k2- k k2 x+ky, and h(+~k) and v(~-k) can be obtained from (1) with k, = (kix + kyy ~ kz)/ko. Substituting (3) into (2), the following integral equation for the polarization current can be obtained: 1(r) = (E + E + ~J(r) = -ikoYo(E' + E) + k2 ~ —1 00 d+ f (x',y') -C. 0 G(r, r')- J(r') dv'. (5) 4

Al aijpprox Il L(.soliltioll for tlie il tegral ('(p1iti011 Ca1 be obtailf(l ied si a p )ertr-;atLion f.(hrliqll<'. 'I 'l total ll iiri,, (i'11 t''f' is ('.Xparlldel inl l.terrls of,a p) jrl)tr atiotl s,('ie(s.. g ivei l),y J(r) = J,,(r)LA, (6) n=0 with the expectation that lim Jn(r) = 0. The most inner integral in (5) is expanded into n-oo a power series in terms of A and then a recursive set of equations for Jn arc obtained. These currents are expressed in term of their two-dimensional Fourier transforms defined by Jn(r) =-2 d2k J,,(kl, z)e P (7) After much algebraic manipulation, analytical solution for the induced polarization current to any desired order is obtained [2]. The expression for Jo(kj, z) is given by Jo(kL, z) =(27r)2 (kj - k) [Joh(z)h, + Jot(z)i, + Joz(Z)Z], (8) where ti = z x hi, and Joh(Z) = - Yo (c-1) Co(k z) [Q hi] e-ik z ' z kpz + kl, ) Jo(z) =- kktkZ -Q Jot(z) = - k (e + (-1) C (k, z) [Q. ] e (9) kl = ko, V=- sin2O, k = k0 sin, RR = k R =' ^ k I ^U cki + ktI Ch(k, z) (-1)n (Rh - rh) eiklz + (Rhrh - 1) e-iklz n p Rh (Rh - rh) eizd + (Rhrh- 1) e-izd C(k,) (-1)" (rv - Rv) ei'kz + (R,, - 1) e-iz-' R/ (/R - rv) elkild + (Rrv - 1) e-ikid As before 7h and rv denotes tlhe Fresnel reflection coefficients of tlhe half-space medium. If the half space dielectric is homogeneous (r/, = Rh and rt = /I?), the values of C7' and C are one. Tle expressions for the first-order currents are similar to those of to those of te zeroth order, and arc given in tlhe Appendix.

3 Evaluation of Scattered Fields Sul)stituting thlc expressions for polarization currents into (3), and following a similar procedure which resulted in (6) of [2], the scattered field to the Nth order in A is obtained: d E'(r, k) -=~ f /2 L {/ G(kl; z, z') J(kl, z') dz' o (27)2J "-1 n (n)n+.l (m ~ n-m n+, mG (n + 1)] z'tm G (k, z d) [z0 JF(k, ) n=O m=O (n + 1)! /-m (10) where G z, z') h(k,) [rhih(-kz)e'.z' + h(kz)e'z' ] + v(k) [r^(-k)eiz + v(kz)e-ik'] }e', if z > z' (11) which is the Fourier transformation of the dyadic Green's function. Equation (10) can be divided into two parts: E'(r, ks, k) =Esf(r, k, kt) + EsT(r, kS, k). (12) Esf(r, k, kic) is the scattered field to the zeroth order in A. Substituting the zeroth-order currents (9) into (10), the zeroth-order field is given by Ei (r, k, ki)) = eiokrr {( )-)~)(k') + (c)(-)R(~)(k)} Q (13) where Q and P are the polarization vectors of the incident and scattered fields, respectively, and R (k) [(1 + R)C(kd)- 1] e-2ikd, R( )(k) [(1 + R)C(kd)-1] e-. (14) The zeroth-order solution is equivalent to the reflected field from the original multilayer medium with a flat interface (f(x,y) = 0), and the expressions in (14) give the total reflection coefficients at the air-medium interface. It should be noted that (,(-k}),h(-k')) = (vihi) and (i'(k'),ih(kt)) = (irhr). The superscript "f" in (13) denotes the flat interface. 6

E"r(r, /",/' ) iS tli( IIigl,,i ',,,,ord sc( tl t(e(d iild, wlinci only cxists wJic' ItJll e surrfalc is roIIgII. SIII)Stiitgui11g (lie JoI li J / Ii( il Ii ii s in1o ( 1)).,,,,1 eor Isomrei fIg(l'el i r;irup,li olt 11. s,.lie' followillg ('XIpI(ssiO. i,,- f1 t,,1 st (fle( l ll i('l is oilitaliled k. -. ',. I n - Q/,',) = 8w2 c. L (N- 7)! N=I " n=O m=n= P (. h( t)(k,) [ h + ( )m]J Cm(k, d) + (kl)t(kz) [R' - (-1 )] C (kd) V -m +k )n [r/)+ -n-m-l, ~I o+~(k)i [l + (-1)z] (+ -1 (k-; d), i} 0-n(k,)) ~ (k j ) (15) Note that (15) is valid for all points z > d, including observation points very close to the surface. In a special case where the observation point is far from the surface (z >> d), the stationary phase approximation can be used for evaluating the integral in (15). The far-field expansion for the scattered field can be written in terms of a scattering matrix which depends on the direction of observation. Comparing (15) with (14) and (15) in [2], E'Q can be expressed in terms of scattering matrix elements: Et(r, k, k) = 2, d ik P {h(k,)h(-k S(v)(k,kN ) +h(k)(-k) 27r N1 J kz sh(k, ki) + v(kz)h(-kz) S (k, k) + (k)v(-k ) S (k, k,). Q, (16) where.the scattering matrix elements in (16) are given in [2]. The near-field expression in (16) can be interpreted as the superposition of the scattered fields from all different directions denoted by k. The integrand corresponding to Ikl < ko can be interpreted as the upward propagating waves emanated from the surface. When Ikll > ko, the corresponding waves are non-propagating which are known as the surface wave whose contributions are confined in the vicinity of the interface. It should be emphasized that the quantities of interest are the statistical mean and the standard derivation of the scattered field. These surface waves are caused by the rough-surface scattering and does not exist when the surface is flat. Performing the ensemble averaging of (12), it is found that, up to the second order, (E(r, kT, )) m E(r, kT k) + (E (2)(r, k, i)). (17) lere (E's(2)(r, r, k')), like ES/(r, k, ki') in (13), can be expressed as (Es(2)(r 1)) io'r p. {( )(-) i) + ()(2)(}. Q (IS) /

where J()(, k) 2 anl /I((/, ( ) are also given in the Appendix. Since f(x, y) is a zero mean GCaussianl process, it calr e)(' shown' that tlie average of tlie odd-order fields vanish. Th'e next term of the coherent field is the fourth-order ESr, which will be ignored due to the assumption of the slight roughness. For the evaluation of the incoherent scattering power (variance of the field), only the first-order scattered field is retained. Re-arranging (16), we have EPQ(r, kk) = d2kIpQ(k, k )AF(k - kl), (19) where IpQ(k_, kkL) is given by IPQ(kk) =- k P {h(kz)h(-kz) S()(k, ki+ h(k)v(-k) S(k k) kzAF(k~ - k). + (k)h(-k') S() (k,k) + v(k)(-kz) S()(k,ki)} Q. (20) Noting that A2(F(k~)F*(k')) = A2(F(kl)F(-kI)) = (27r)2 2(kl - k')W(kl), (21) the incoherent scattering power, up to the second order in A, is given by (IE - (E )12) ( Esr( 2) = A2 Jd2k IpQ(kk'lk)2 W(kl - k). (22) As mentioned previously, (16) is expressed as a continuous spectrum of scattered plane waves. What is expressed mathematically by (21) indicated that these plane waves are mutually uncorrelated. Therefore, (22) is simply the integration of the power carried by each plane wave. For observation points near the surface, (22) must be carried out numerically and cannot be simplified any further. The convergence of the integral can be examined noting that W(kl) decreases as Ikll increases and the fact that for [kjl > ko, kz becomes pure imaginary which causes the integrand (JIpQ(kl, kl) 2) to decay rapidly. To demonstrate the effect of the surface roughness on the surface reflectivity, a numerical example is considered. Both coherent reflectivity and incoherent reflectivity (( PQ )/ IE'I ) as a function of observation point height are calculated. These plane wave illumination examples simulate situation where the transmitter (receiver) is airborne and the receiver (transmitter) is near the rough interface. Consider a rough soil surface with rms height of 0.016m, correlation length 0.16m, and dielectric constant E = 8 + ii illuminated by a plane wave generated by a source operating at f = 890 MHz. At this frequency, the normalized rms height and correlation length are, respectively, ks = 0.3 and ki = 3.0. Figure 2 shows the magnitude of the zeroth order and complete second order mean-field in (14) versus incidence angle. In this simulation, the correlation function S

o()1 Llly Io(,l g1).sl1r fa( ('(. i. S..ljll ('(j (;Jll.'.,iLll. 'I '}(; (()IIIl)J('t(' s('('() flI ()I'(j('r.()lJiOl) d(']'lO lfor thie mui suac I Is a-Ssiiumed (iauiss5. '& I li- ohe soormp ete scmi (1 der solutiton deiimo01 -st.r;f.tes ( e. f c.lfe f[l t of.sir;face(( roughnes I (.s ol.li.( sirl;f(c('e ri(elle. ivitl.. I i.,sir(; ill, tf.h.surface oDIi li(ss r Iu s.l il suilface r l( Cti vil y a-1(1d (use;s('S A sliilght shift iII 1-i 1t 1 w t.(' aigle(' It sh(iouIll I(e rioted th.ii forr il;ation of the se('cond or(l(der coherent. re'flection (co('lici('elt. ill (18) does (not c(onve(Irg(e for( silrfaces wiith exponlential (correlatioll fulrictioln. j'hlis mlay 1be due to tihe fact that ligler ordcr terms are excluded. 11owcver, tliis problem is not observed in tihe formulation for the incoherent wave. Figure 3a and 3b show conmparisons between the zeroth order coherent and incoherent reflectivities as the function of the height of the observation point for an exponential correlation function. It is shown that for vertical polarization and observation point heights less than 0.1A, the incoherent reflectivity is significant and dominant near the lBrewster angle. However, for horizontal polarization, independent of the incidence angle, the incoherent reflectivity is much smaller than the coherent reflectivity. Generally, the incoherent reflectivities decrease as incidence angle decreases. This could be qualitatively explained by Rayleigh criterion [4]. The criterion is stated as follows: For a surface characterized by a distribution of irregularities of height h, if h satisfies A h < ---- (23) - 16cos ' ( where 0 is the incidence angle, the surface can be considered smooth. As the incidence angle increase, the surface appears "more flat". Therefore, the incoherent scattering decreases. Same is true for the coherent field as shown in Fig. 2 where the coherent reflectivity approaches unity when 0 is increased to 90~. It is noticed that the incoherent reflectivities vary as the height of the observation point changes. As mentioned previously, the scattered field can be decomposed into two components: upward propagating waves and surface waves. When the height increases, while the surface wave components attenuate, the propagating waves remains unattenuated. This phenomenon is demonstrated in Fig. 4 where the integrand of (19) is plotted in k1 space. The normalized magnitude of the integrand is shown in gray scale over an area with the radius of 2ko in the spectral domain. The propagating waves are confined in a circle of radius ko, while the surface waves are outside the k,-circle. Figure 4(a) and 4(d) show the incoherent vv- and hh-polarized power spectral densities for an observation point 0.01A above the rough surface when the incidence angle is 200. The integrand is normalized with respect to the value at (kr., k) = (ko sin 18.80, 0.0) for vv-polarization and (ko sin 21.2~, 0.0) for hh-polarization. In this case, most of the power is in the ko-circle, which justifies the lack of sensitivity of the incoherent reflectivity to the height variation at 20~ shown in Figs. 5(a) and 5(b). Figures 4(b) and 4(e) respectively show the incoherent vv- and hh-polarized power spectral densities when the observation point is 0.01A above the rough surface at incidence angle S0~. The integrand is normalized with respect to the value al (kxky) = (-ko,0.0) for vv-polarization and (ko,0.0) for hh-polarization. A significant. component of incohllrent scattering is from the contribution of the suralface ()

waves (thi area outside tile k,-circle. It is also noti(ced thlat for vv-polarizationl incohlierent scattering is mostly from t hie waves aroud(l thie backscattering (lirectionl. As the lheight of the observation point increases to 0.5A, the contribution from thle surface waves almost vanish, as shown in Figs. 5(c) and 5(f). In near-field region, all field components are present in general. Decomposing the field components into v, h, and k components, the behavior of the cross-polarized incoherent scattered waves are demonstrated next. Figures 5(a)-(d) show the incoherent cross-polarized scattering, including vh, hv, kv, and kh polarizations. Like the copolarized scattering, the cross-polarized incoherent scattered power is stronger for observation points close to the surface. Also for v-polarized incident field, the cross-polarized scattering powers are stronger than those of the polarization of h-polarized incident field. 4 Evaluation of Field of a Short Dipole above a Rough Surface Another problem of practical importance is the characterization of the field of a short dipole above a rough surface. Consider an infinitesimal current element given by Q 6(r - r') where Q denotes the polarization of the dipole antenna, and r' represents the location of the dipole. At the observation point, r, the direct radiated field from the dipole is given by [5] Ed(r,r) iZ [Q1 + (qk- R) R2] [3R (24) where R = r - r' and R = |RI. For z < z', (24) can be expressed in terms of integral of plane waves given by -koZo 2 ki z -1 Ed(r, r') = 8 Jd2k k [h(-kz)h(-kz) + v(-kz)v(-kz) Q e-k(-c') (25) The mean field can be obtained by evaluating the integral of the mean fields corresponding to each individual plane wave and is given by -Eor) / d2ke -iK E r K) (26) 10

wli.c. k k - /,-.',,,- I.l K kj - /:. S l.S' if,4 ti, (1:) into (2';), lf ' foIi,,l ti, io( I(/,'(r r')) - - (r, (- d - ')-! r') ~ 8' 1 IIdIi,, Z0 {(I - h)) + (I + I?,)( 'o(kp, d)C Q(L P)} (27) Equation (27) is obtained by noting that the reflection coefficients arc azimuthally symImetric and therefore the integration with respect to ) is carried out analytically and are expressed in terms of PQ(kp) and Pk(kp) which are given in (A.3) and (A.4), respectively. The integral in (27) is a Sommerfeld type. When botll r and r' are close to the ground and far from each other, the first term on the right-hand side, which can be viewed as the negative image source, dominates. This results in a destructive interference with the direct wave, and hence the surface waves, which is accounted for in the integral of (27), become dominant. The integral in (27) can also be written in terms of asymptotic expressions available in literatures [6]. The numerical technique for the evaluation of the Sommerfeld integral is not discussed here. Interested readers are referred to [7]. Here the objective is to investigate the significance of the rough surface which are included in the integral in (27) and in the incoherent scattered field. The first-order incoherent scattered field is written in terms of superposition of incoherent scattered fields: -koZo f, iK-r' ~ ) - 872 _I(28) EPQ(r, r=' 82'k Es r k K), (28) from which the incoherent scattering power is obtained and given by 2 Zk Z f2 42[ g2 ed(K")* -K'.r' (I pQ (r, r) 6d2k d / d2klpQ(kl, ki)IPQ(k + kl - k, kl)W(kl - k). (29) Note that the integral in (29) is six-fold which is extremely difficult to evaluate numerically. Practically, the distance between the dipole and the observation point is large, which can be used to simplify (29). Suppose r = z7 and Ir'I >> A. Using the stationary phIase approximation to evaluate the integrals onI k and kV, the incoherent scattering powerl can be obtained from (Eg(r, Z )2) O d2 k /p(k1 z-o') I4/(kl + ko) (30) o7r2 Jr 12|Irl Irl II

Nlrnerical sirnulatioris have been performe(d to dlemonstrate the effect of the surface rouglhness. Consider the rough surface of previous example witlh parameters ( = 8 + il, ks = 0.3, kl = 3.0, and exponential correlation function. Suppose the infinitesimal current source is placed at 0.2A above the surface at a location r' = x'x' + y'y' + 0.2A', and the observation point is at (0, 0, 0.2A). Figures 6(a) (Q = P = z2) and (b) (Q = P = y) show different components comprising the received power versus the radial distance (vIx2 + y'2) between the current source and the observation point. It is shown that, as the distance increases, the reflection coefficients approach to -1, and the coherent ground contribution cancel the direct wave. Note that as the incidence angle approaches 900, the effect of the surface roughness on the coherent field becomes insignificant, as explained previously when stating the Rayleigh criterion. Under this circumstance, the dominant propagation mechanism is the coherent and incoherent surface waves. H-polarized surface waves attenuate rapidly [8]. Therefore, in Fig. 6(b), the surface wave is less significant than that in Fig. 6(a), and the total field shows obvious interference phenomenon between the direct wave and reflected wave. The incoherent rough surface scattering is found to be somewhat insignificant in both cases, which was also predicted in Fig. 3 at high incidence angles. In the next simulation, the distance between source and observation point is fixed at 20A, but the source point is moved on a circle of radius 20A in the x-z plane, as shown in Fig. 7. As before, the observation point is at (0,0,0.2A). However the source is at (20A sin 0, 0, 0.2A + 20A cos 0) with 0 e [0, 90~]. Figure 8(a) shows the coherent and incoherent powers with Q = P = h = y. When 0 approaches 90~, the direct field and ground reflected field interfere with each other destructively, but the total coherent field is still about 10 dB higher than the incoherent field. In Figs. 8(b) and (c), choosing Q = v = (-sin0,0, cos0) whereas P = x in Fig. 8(b) and P = z in Fig. 8(c). and Fig. (c)( z) show the components of the received power. Note that in both simulations, the polarizations at the observation point are not suitable to receive the ground reflected waves, which should be (sin 0,0, cos 0). Thus, the direct field dominates. When Q and P become perpendicular to each other, the coherent field diminishes, and the incoherent field becomes significant. 5 Conclusions The radiation of a short electric dipole above a slightly rough surface is studied. This investigation is a natural extension of the classical Sommerfeld problem with the exception of the random surface irregularities at the interface between the two dielectric media. In this paper, the formulation for the near scattered field from a slightly rough surface when illuminated by plane waves is developed first. A perturbation technique is applied to solve the integral equation for the induced polarization current. Analytical expressions for the coherent field (mean-field) and incoherent scattered power at an ar 12

iitr.r'y o( S jIvti(oli )oiiit., i'l Iifig ploirlt.s ie;ar the iitL face 'L(, arc( oltlainled. SilmJilatiosf. show' th1a. while tlhe coellt s(tet.t('Red fi(ld generally (dIaiiitfs, (.,Slir te i (icoric-ct s(-l..r('(l fi('ll iS oily, e'sse( l.ill.rodlli lf I iw r(' wslt. r a('. le' fo[ r Vv - l(l,.iiztl.iO. 'I lie p)h(i oiio('ii ' of (depolarizatiorn caUtised( by t.1e incoherenii, rotigli-su rface S(catering ar( also studied. (enerally, lie incohlieretL scattered field becorime more signiificant as tlie observation )oint app)roaclhes the intcrfa(c. f'lieni lic solutions for the meanri-ficl( arid incoliherent scattered power generated by a small dipole of arbitrary orientation and position arc derived by expanding the field of the dipole in terms of a continuous spectrum of plane waves and using superposition. Although it is found that the direct and coherent reflected (reflected + surface waves) fields are dominant in most cases, the incoherent scattering could be inmportant, when the pathl along the line of sighlt is obscured. Acknowledgment: This research was supported by Army Research Office under contract DAAG5S-98-10-0458. 13:

Akppendix The expression for the Nth~ order polarization currents are given by JNh(kL, z) = 0k (c- - 1) Ch)(k, z) [VN, h (k,)1 kz + ki z Z'koklz (c - 1)Cj(k)[V i(k ] JNt(k~L,z) = ck, + kj, Cvk1z ~N- — z I J JN.(kiL,z) = -kkP lCv(k, ) ckz + kiz [VN 'kz where VNN-i N-n-I n=O m=O (N-nlI1) (i'kz)m [ N-n-m-1 1 N-n *O(F(k~L)I where * is the convolution operator, F(k~L) is the Fourier transform of f (x', y'), and0 n n represents n-fold self-convolution (O&F = F * F ****F). The second-order expressions for reflection coefficients used in (18) are given by R(2) =k 2ik S2d) + oh [c(~)k-k )] 2 Rh kz -2 khkk,%-Co(kkddCjsin, () - q7k)]2 Jd2k' W(k' - ki) IC'(' k~o2q + kz k~kD) i (A. 1) 2- ck2) R$J2) =k -2ikzd{ 2 ( - 1)(1 - +v (i ~ p "C(k, d)Ci'(k, d) + 2 r(C - 1)2{(1 _RV) 2Cov(k, d) I d2kI W(kI - k ~) [ kCoh(k', d)Cov(k, d) sin 2(qf - sbk) k'k ~ k' )~k, d co2(qk-~ 0 7 - k + z k,Co k',kd)Cov (k', d) c s(Ok, d) cos 'k - k + k+ k(ke~k +)kJ __ __ __ __ __ _ __ __ __ __ _ k k C Ik' dCo (k', d) C os2(k 4kd osO - O + (2Rk Cv(k, d) ckz + ki I I d2kI W(kI - k-L) [k'-'yk1C v (k', d) Cov (k, d) cos (q$ -kCvkq dCvkk) +S(O O'kJJ j(A.'.), 1.1

The' filtegritdI1(In (27) itrc giveii b~y 0 = (PrQ1r + PvV J10(k.,(jp- -,YJ1)) +7r (PrQr, - P~Q,)cos(2q0') + (PZQ, + PyQ,) sin(24/)]J2(k,(i,~I)(A.-3) 27r 0 - [2kP Q - k (P Q. +PQY)] Jo (ko (I, -'I)) + k2 [PQ P1 o S + ( Q - 0 + wr 2o [(PXQX - PyQy) cos(2qY) + (PxQY, + PyQZ,) sin(2qY)J J2(kp( j, - f') (A.4) di=Cos1

References [1] Sonirmerfeld A, lirtial lDifftr'ntiail Equations in liysics, New Yorik, Academic 1)Press, 1985. [2] Sarabandi 1K. and T. Chiu, "Electromagnetic Scattering from Slightly Rough Surfaces with Inhomrogeneous Dielectric Profiles," IEEE Trans. Antennas Propagat.., vol. 45, No. 9, pp. 1419-1430, Sep. 1997. [3] Tsang L., J. Kong, and R.T. Shin, Theory of Microwave Remote Sensing, John Wiley and Sons, New York, 1985. [4] Brown R.G., R.A. Sharpe, W.L. Hughes, and R.E. Post, Lines, Waves, and Antennas: the Transmission of Electric Energy, the Ronald Press Company, New York, 1992. [5] Balanis C.A., Advanced Engineering Electromagnetics, John Wiley and Sons, New York, 1989. [6] King R. W. P. and S. S. Sandler, "The Electromagnetic Field of a Vertical Electric Dipole over the Earth or Sea," IEEE Trans. Antennas Propagat., vol. 42, No. 3, pp. 382-389, Mar. 1994. f7] James J.R., and P.S. Hall, Handbook of Microstrip Antenna, Peter Peregrinus Ltd., London, United Kingdom, 1989. [8] Hall M.P.M, L.W. Barclay, and M.T. Hewitt, Propagation of Radiowaves, the Institution of Electrical Engineerings, London, United Kingdom, 1996. 16

List of Figures I An inhomogeneous hlalf-space medium with a rough interface. Left side of this figure shows the dielectric profile..................... 17 2 The magnitude of the reflection coefficients in (14) versus incidence angle. The underlying ground is homogeneous, and the dielectric constant is c = 8 + l........................................ 18 3 The comparison between the coherent and incoherent reflectivities versus incidence angle. The underlying ground is homogeneous, and e = 8 + il, and the rough surface has exponential correlation function, ks = 0.3, and kl = 3.0. The incoherent reflectivities are plotted for different observation point heights..................................................... 18 4 The spectral distribution of the incoherent scattered power generated from the rough surface. (a) and (d) are respectively for vv and hh polarization at incidence angle 20~, and the height is 0.01A. (b) and (e) respectively for vv and hh polarization at incidence angle 80~, and the height is 0.01A. (c) and (f) are the same as (b) and (e), but the heigh is 0.5A.......... 19 5 The cross-polarized incoherent scattering fields versus incidence angle. (a) kv, (b) kh, (c) hv, and (d) vh. The incoherent reflectivities are plotted for different observation point heights..................... 20 6 The electric field versus the distance between current and observation point. The underlying ground is homogeneous (e = 8 + il). The rough surface has exponential correlation function, ks = 0.3, and kl = 3.0. In (a),Q-= P =. In(b)Q=P-............................ 21 7 The distance between source and observation point is fixed at 20A, but the source point is moved on a circle of radius 20A in the x-z plane. The observation point is at (0, 0, 0.2A), and the source is at (20A sin 0, 0, 0.2A + 20A cos 0) with 0 e [0, 90~]....................... 21 8 Components of received power as a function of transmitter height for a fixed distance (20A) between the transmitter and receiver. The observation point is at (0, 0, 0.2A), and the source is at (20A sin 0, 0, 0.2A + 20A cos 0). In (a) Q = P = h = y, and in (b) and (c), Q = = (-sin 0, 0, cos 0), and P = x and P = z, respectively......................... 22 I

z Incidence Field 9i ',,: 1 \ \ v \%\ ~(z) dA --- —-- 0 x ~; *;. -.. Figure 1: An inhomogeneous half-space medium with a rough interface. Left side of this figure shows the dielectric profile. 17

0. I ' - I - l - - l - l - I - f0 I — x C7) (U 04 u -5. 10. -v -V -15. -20. -25. -30. Rvv, zeroth-order Rhh, zeroth-order Rvv, complete secon - - v - - Rhh, complete ~ - I i i i I i I _ -35 I...... 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. Incidence Angle (degree) Figure 2: The magnitude of the reflection coefficients in (14) versus incidence angle. The underlying ground is homogeneous, and the dielectric constant is e = 8 + il. 0. I *. s.......... I-I-I-I-I. I i I I 4-. - "0 C) -10. -20. -30. -40. -50. coherent i-' - -— ~ — inco, z=0.01 A — I --- inco, z=0.1 X - - v - - inco, z=0.2 X - -O- - - inco, z=0.5 X ----- inco, z=1.0.. i. i. I. I _ i. i PQ '0 C) X z 4-. '. 4 --- a) I nr\ I i i -oU. -... -.. 7.. 0. 10. 20. 30. 40. 50. 60. 70. 80. < -60..-..-.. - J... -. -. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 90. Incidence Angle (degree) Incidence Angle (degree) (a) (b) Figure 3: The comparison between the coherent and incoherent reflectivities versus incidence angle. The underlying ground is homogeneous, and c = 8 + il, and the rough surface has exponential correlation function, ks = 0.3, and kl = 3.0. The incoherent reflectivities are plotted for different observation point heights. 18

KvJ 2 Ko - I IKo Ky 2 Ko 1, %. I I I I I I I i I I i 1 Kx I I I I I I 1. I I Kx I I I 0.0 0. 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 (a) VV, 0- = 20',h=O.O1A 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 (b) VV, 0, = 80',h=O.O1A 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 (c) VV, Oi = 80',h=O.5A KYJ 2 Ko Ko. I I I I 'I.4 I I I I I I N I I I I I I i / — 14 I / I/K Ko - - - —! 7,.. -. — - % I..... 10 % 'I.-- I I.... -: V.,. ;, - -, '! -;.'- ".:... % ol %.1 /1 0.0 01 0.2 0.3 0.4050.07080910 (d) IIH, 01 = 20',h=O.O1A 0.0 0.1 020.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 (e) HH, 0i = 80',h=O.O1A 0.0 0. 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 (f) 111, 0i = 80',h=0.5A Figure 4: The spectral distribution of the incoherent scattered power generated from the rough surface. (a) and (d) are respectively for vv and hh polarization at incidence angle 200, and the height is O.O1A. (b) and (e) respectively for vv and hh polarization at I nci'dence angle 800, and the height 'is O.O1A. (c) and (f) are the same as (b) and (e), but the heigh is O.5A 1 9

~ -0 >I A@ A-a.- c) - 0. I -30. ---- z=O.Il e -40. -40 — ^ — z= 0.52 A — v — z=O.5X -50. - --- z= 1.0 O ---- z= 2.0X -60..........111 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. Incidence Angle (degree) la C.-) 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. Incidence Angle (degree) (b) (a) OE Z A/ >-.4 -C) '4-i 0 QU 04 10. -20. -30....... -40. -50. -60. -1....... 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. Incidence Angle (degree) C-) 04 -10. -20. -30. -40. " -50. -60. ''0. 10. 20. 30. 40. 50. 60. 70. 80. 90. Incidence Angle (degree) (c) (d) Figure 5: The cross-polarized incoherent scattering fields versus incidence angle. (a) kv, (b) kh, (c) hv, and (d) vh. The incoherent reflectivities are plotted for different observation point heights. 20

-l.e - 70. \-. X-'5 -5 ~0 V.- -70. { u u - — 70. ~ V. -90 ~ 'total. coherent -90.... -- dire ct XlO. \ - ground coherent - 30. - o. 10. - o 30.- o 20. 30. Distance () 10. Dsac Distance (X) (b) (a) t aL observation point. u tnebet-ween curren sraceha exponential field versus the distance betTeen rcuren t 20X observat io n oin. Figue 6undeerlyi ggon is homogeneous ( -8+ il.Teruhsrac a xoeta The underlying ~c= 0.3, andKcorrelation function, ks \ ^\ z~ t i fiedat 90N, but the' source a dobservattion point The, observation point betw~e So)' h - lane.t e o,'ue7 '~distance betw ee —:s 0inte-2,osO witl,06[,O1 xe on a, cirl u u 0, mo.,.t / s,s^ source potn[ Is i ouace isia is at (0,, 0.2X), an:d hr,' '2l

f I -30. -40. -50. — v- -— s -60. 7C 7 - -o- total, coherent ', -— 40 --- direct — I --- ground, coherent -90. - -, - - ground, incoherent -100. - -'.,.......,.. 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 0 (degree) (a) -30. 1 -30. 7 -a-1-40. 7O -40o., 50. 7 -50. ':, v-v-v-v-v..., V.-60. -v -60. -70. ---- total, coherent ' L -70. t total, coherent -80. - ' ' direct 4 -80. " "- direct J -9| —& --- ground, coherent -- ound.cohent -90. -7 0. - - --- -ground, incoherent- ground, incoherent -100. - I - -- I --- - I. * *.... 0. 10. 20. 30. 40. 50. 60. 70. 80. 90. 0. 10. 30. 40. 50. 60. 70. 80. 90. 0 (degree) 0 (degree) (b) (c) Figure 8: Components of received power as a function of transmitter height for a fixed distance (20A) between the transmitter and receiver. The observation point is at (0, 0, 0.2A), and the source is at (20A sin 0, O, 00.2A + 20A cos 0). In (a) Q = P = h = y, and in (b) and (c), Q = v = (-sin ~, 0, cos 0), and P = x and P = -, respectively.

Effect of Canopy-Air Interface Roughness on HF-UHF Wave Propagation in Forest Kamal Sarabandi, Il-Suek 1.olo Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan, Ann Arbor, S1I 48109-2122 Tel:(313) 936-1575, Fax:(313) 647-2106 email: saraband eecs.umich.edu Abstract In this paper, the problem of electromagnetic wave propagation in a realistic forest environment at HF through UHF bands is considered. In particular, the effect of the non-planar interface between the air and the canopy, which has been ignored in the previous models, is examined. An analytical formulation is obtained for the mean field when both the transmitter and receiver are within the foliage. This formulation is based on distorted Born approximation and accounts for the surface roughness that exists between the canopy and air interface. It is shown that the surface roughness attenuates the so called lateral wave which is the dominated source of the field at the receiver locations far from the transmitter. It is also shown that this attenuation rate increases when the rms height of the surface roughness is increased. 1 Introduction With recent progress in wireless technology and ever increasing demand for reliable low power wireless systems, the need for predicting system performance has been become an extremely important step in the design of such systems. At HF through UHF, the channel characteristics such as attenuation, and multi-path fading statistics significantly affect the performance of wireless systems. Prediction of the channel characteristics can be accomplished using physics based propagation models. For this purpose precise diffraction models must be developed and incorporated with accurate terrain data base to construct a realistic propagation model. Since a large portion of earth's surface is covered with vegetation, understanding of electromagnetic wave propagation through foliage is of great importance in devel — opment of a comprehensive physics based propagation model. It is well known that 1

at HF through UHF frequencies the direct propagation betweeni two distant points withinr a forest experiences significantly less attenuation than that predict(ed by models. That is the predicted attenuation of the direct wave far exceeds the attenuation values obtained from experimental results. A very interesting model that describes the phenomena and attributes the wave propagation in forest to lateral wave was first developed by Tamir [2]. In this approach, the forest is modeled by a homogeneous half-space dielectric medium with a planar interface and a permittivity equal to the effective dielectric constant of the foliage medium. This effective dielectric constant can be obtained using a dielectric mixing formula [3] or the formulation for calculation of propagation constant in random media [4, 5, 6]. The former is appropriate for low frequencies when typical dimensions of the constitutive particles are small compared to the wavelength whereas the latter is useful for sparse media and accounts for scattering losses. Tamir's formulation gives the field of dipole within the effective medium at an observation point near the interface using an asymptotic integral evaluation [11, 12]. This solution shows that the field at the observation point is dominated by a ray emanated from the source point and traveled in the direction of the critical angle toward the interface and along the interface before leaving it in the direction of critical angle, in order to reach the observation point. Due to renewed interest in wireless communications at HF through UHF, recently the problem of wave propagation in forested areas has gained some prominence. The origin Tamir's model is extended in [7, 8, 9] by representing the forest by a two-layer anisotropic dielectric medium. These models accounts for the anisotropy that exists in the trunk layer and to a much lesser extent in the canopy layer. It is shown that upto UHF frequencies the effect of anisotropy can be ignored [7]. In reality, the interface between forest canopies and air is not flat, hence it is not clear whether the lateral wave can be excited and of it can how the surface roughness affects it. In this paper, the effect of roughness of interface between canopy and air on the wave propagation in forested area is investigated. In section 2, analytical formulation using the volumetric integral equation in conjunction with the distorted Born approximation is presented. An expression for the mean-field of an infinitesimal dipole of arbitrary orientation is derived by obtaining a partial second order solution of the Born approximation. In section3, a sensitivity analysis is carried out to demonstrate the variations of the mean-field to physical parameters such as a effective permittivity, location of the dipole and observation points, and surface roughness. 2

2 Analytical Formulation In this section. an analytical formulation for the evaluation of mean field generated by a small dipole in a homogeneous half-space dielectric medium with rough interface is derived. The effective permittivity of forest canopies is slightly higher than that of free space since the vegetation particle density is very small [2]. In this case, distorted Born approximation may be used to evaluate the scattered fields. The geometry of the diffraction problem is shown in Fig.l where the dipole and the observation point, are respectively located at h and h' inside a canopy with effective dielectric constant E1. The envelope of the canopy-air interface is denoted by a random process, d(x, y). The permittivity of the upper medium(air) is denoted by E2. We modify the problem by extending the canopy to z = 0 plane and assume that in this region, there exists a volumetric polarization current J = iklYe (2 - 1 )Et, where Et = E E + E E. The unknown total internal field can be obtained from the following integral equation roo roo r0 Et(r) =E()+E(r)+ ko(2 - E) / G(r, il)Et(')dx'dy'dz' oo J -oo J -d(x,y) where El is the incident field, Er is the reflected wave from the planar interface, z = 0 and, Es is the scattered field generated by J itself and G is the dyadic Green's function of the half space dielectric medium(see(2) [1, 10]). To the first order in (62 - 61), the Born approximation can be used to the scattered field roo roo rO E = ko2(62 -E) G(, r') [EI(r) + Er(r )]dx'dy'dz' (1) J-oo 0 o J-d(x,y) In (1), d(x, y) is a two dimensional random process describing the interface between the canopy and air and is assumed to be Gaussian with a mean value of m(a positive number) and a standard deviation of a. Noting that the z-dependence of Ei and Er a f te fm o e) where of, the form of exp(integration with respect to z' can be carried out analytically. This integration would result in an algebraic expression in terms of d(x, y) which is amenable for the calculation of statistical average of E5. It is shown that distorted Born approximation provides a more accurate solution for Es. In this approximation, a phase correction term is incorporated into the approximate expressions for the internal field to account for the difference in the propagation constant between the air and the canopy. The explicit expressions for the incident and reflected fields, including the appropriate 3

phase correction terms for a plane wave illumination. are given by i = P( )( i[rk' l +k y'+kL-' ]. -i(k -k ) (,) -' =.i -k)'y'-k '] ei(k-k)m (3) where P(kl:) is a unit vector indicating polarization of the incident wave (P = v or h) and Rv,h is the Fresnel reflection coefficient for vertical or horizontal polarization. The unit horizontal and vertical polarization vectors are defined with respect to the plane of incident and given by h(klz) = X- h(-k1z) = h(k1z) jk' X v(kz) = h(k1z) x k1. V(-kz) = h(kz) x k where k[ = k - 2(kl x Z)z denotes the direction of the reflected wave. Since the layer thickness is not uniform, the phase correction terms in (2) and (3) are chosen for a layer with a uniform thickness m. The dyadic Green's function for the half-space dielectric medium when z < z'(observation point inside the foliage) is given by [1, 10] 6( ) + i j j k1 _{h (F-Fk) [Rh/(kl/). G = - zzx) ++ k (y (-( 2z2 + 8TT2 JV dkzdkye[k (z-z)+(Y-y)]k ^ Ih( e-'ikjZ + h(-klz)ek1zz] + v(-k) [Rv(klz)ei + vk(-l,)eik"zz'] } e-ik (4) After substituting (4) into (1), the integration with respect to z' is carried out analytically. After performing this integration, evaluation of (Es) is attempted which requires computation of the term like (eJsd). Assuming a Gaussian p.d.f. for d, it can easily be shown that ess ~'2s2 OIQ2S2 (eisd) = e [e erfc(-jc) + e 2 erfc(j ) (5) (e-js) = e2 -erfc(-j ) + e2 erfc(j i)] (6) where complex variable, s can be either s1 = kz + klz or s2 = k2 - k1, and erfc(.) is the complementary error function defined by [13]: erfc(z) = 1 - [j e-dt + je-2 et2 tdt V/0 / \] 4

where z = a + j,. Noting that erfc(jx) = 1 - j Sf_ 2dt for real x. it can readily be showll that (e-jsd) reduces to the more familiar expression, (e-jsd) = (-J'Sn-"'2 The integration with respect to x' and y' can now be carried out analytically which results in S(kx-k). (ky-k'). This, in turn, simplifies the integration with respect to k, and ky. The final result is a plane wave propagating along kr. Hence the total field in medium 2 is the sum of two plane waves: 1) the reflected wave at the hyNpothetical planar interface at z = 0, and 2) the mean scattered field. Hence to the first order in ko2(E2 - El), E = E ( (Rr + R(E) (RIref+ ))e[kx+kY-'' (7) where Rref is the Fresnel reflection coefficient for the canopy-air boundary. After some algebraic manipulations, the first order reflection coefficients (referenced to z = 0 plane) for the mean field are obtained and given by R (2(1) _(2- l) 1 + R(eis) - ( -isld)) kP - kl Born(v) - 2k~ v + ((ei2d)- (e-iS2d))ko2 eis2m S2 ( ko 2(2 - l) -1 R2 2 2(eis-d * — l$d Born(h) - 2k 1 - + sRh(e ( )) (9) +. (S2d) (e-i2d))} eis2m S2 Close examination of R(1) reveals that Rref + RBorn is accurate enough for horizontal polarization, but it cannot accurately predict the Brewster angle for the vertical polarization. To rectify this difficulty, higher order solutions must be obtained. This is accomplished by obtaining a partial second order solution for the vertical polarization. The second order solution for Es is given by s(2)() = k(52 - )2 G( i ). { G(r, f') [Ei(') + r(f)]} dv'dv (10) Using only -zz k2) term of the dyadic Green's function (see[4]) in the integrand of the inner integral, a partial second order solution is obtained. - k_(.5) 2 2 (( Er())] d' (11) /: Jv I 5

The cornplexity of this integral is of the same order as that of (1) aInd the eInsemI)le average of E$(1"5) can be obtained in a manner similar to what was used in computation of (Es(l)). Including the partial second-order term, the reflection coefficient for tlhe vertical polarization is modified as follows. Rt..5) 2 I R I2 +-1( - R2 + R (eisd) -(e-isd)) (-k - k) }eis S1 E2 \ 2 i2 where kp2 = kI + ky. Now it can easily be shown that at Brewster angle where R= 0, R(1o() vanishes also. Having found an approximate analytical solution for a plane wave excitation, the solution for an infinitesimal dipole embedded within the foliage can be obtained by expanding the field of the dipole in terms of a continuous angular spectrum of plane waves. Without loss of generality, consider a dipole, whose orientation is denoted by a unit vector 1, located at ro = -(h + m)z. The field of this dipole within a uniform medium with permittivity cl can be expressed by [10] Ed() 8r2 - I ok 0kZj_[_ k[ ]ki(m e} k dkd (13) Ed(l 7 r / z 1 J- 1 - E = 0^ 00 /00 Z (3 where I is the amplitude of the sinusoidal current carried by the dipole. The integrand of (13) can be regarded as Eoeik'r where o= Ik~Zo 1 - (k 'I)ki] is the amplitude of each plane wave propagating along k = kl/kl. Using superposition, the mean scattered field in the presence of the upper free-space medium with rough interface, can be computed from the coherent sum of all reflected plane waves. That is t 7 f ~~ o 1 i i A5)h 81 0o o k00- [(1 vi(kl))vi(klz)R +I (1. (-ki ))hi(kz ) 1 Bh (E') Eref - 8w- 8 2 _Z l zBv + z lz- BhRj 0 J-ooJ0-000 - e-ikl(z -m-h)ei(k'x+ky)dk d (14) 6

where Eref is the field of the dipole in the presence of the uipper free-space mediium1 with a planar interface at z = 0. The expression for Erjf is similar to the integraild III (14) with the exception that RB is replacedl with Rref. Using the change of variables x = p cos, y = psinp k = kp cos 7, k'y = kp sin l and making the use of the following identities 2 cos mveikppcos(~-)d -- = 2i m cos 7nJm, ( kpp) /27r j sin mveikPp cos(-')djv = 27im sin m Jm(kpp), o the integral in (14) can be simplified to IkZ Z k - eiklz (z-m-h) tR(1) E' = ~~ --- dkpP {(Jo(kpp) + J2(kpp)cos 2p)lr 47r k i lz 2 R(i.5) k2 + J2(kpp)sin 2yly}+ - {- z(Jo(kpp) - J2(kpp) cos 2p)l, /0 2 (1) + z J2(kpp) sin 2oly - ikpklz cos Ji (kpp)lz} + Bh {J2(kpp) (1-5) k2 sin 2Plx + (Jo(kpp) - J2(kpp) cos 2)ly} + 2J(kpp) k1 2 + 2 [ikpki(cos ly + sinoly)J (kpp) + kJo(kpp)lz]z} k2 J (15) Obviously a similar expression can be obtained for Eref. As will be shown later, the accuracy of the distorted Born approximation degrades as the mean height m and/or the normalized permittivity difference " `2 increase. To rectify this problem to some extant, we use the fact that an exact solution is available for a planar dielectric interface. In fact for a plane interface(7 = 0) at z = -m, for the reflection coefficient 7

of each plane we have Rexact(o = 0) = Rref eI2k Hence, for a rough interface, we may write Reract({) - Reract(o = 0) - RBorn(UT = 0) + RBorn(a) (16) and substituting (16) into (15), the modified solution is given by: -* -* -S '(Elct()) Eeat( = 0) - Eorn(a = 0) + (Esorn(O)) (17) 3 Asymptotic Evaluation In the previous section, an closed form solution for the mean-field of a dipole inside a dielectric medium with rough interface was obtained. When the radial distance between the observation point and the source point(p) is large, the integrand becomes highly oscillatory and therefore accurate numerical evaluation of the integral becomes extremely inefficient. This is especially true when both points are near the interface. The standard approach is to change the contour of the integral of integration by first extending the limit of the integral over the entire real axis(-00,00) and then using a change of variable kp = k1 sin w. An approximate analytical expansion can be obtained by applying the standard technique of the steepest descent [11, 12]. The expressions for the reflection coefficients, in the w-plane, take the following forms cos w - V sin2 w K cos w- - sin2 w Rh-=, R= cos w + V/ - sin w n cos w + V/'-sin w where K = ~- < 1. These introduce a branch cut (with the branch point at Wb = sin-1 V/) in the w-plane which may be encountered by the steepest descent path. In this case, diffracted field is dominated by the saddle point and the branch cut contribution. Hence, in general, each component of the diffracted field may be written as Ei Is + U(Os- Oc)IBc (18) where Is and IBC represent the saddle point and the branch cut contributions respectively and U(.) is the Heaviside step function, Os is the saddle point, and Oc is given by 0c = Re{wb} - cos-1 [sech(Im{wb})] 8

The saddle point conitributioll call be obtained rather easily. The saddle poiiits coIrresponding to the integrand of (14) is wt, = 01,2 where sil 01,2 =, c 01 h+= h cos02 = h+'+2 = /p2 + (h + h')2, and 7-2 = v/p2 + (h + h'+ 2w))2. The sa(dle point contribution is found to be -,sad _j koZo [ e-j3k 2 E. -. 0 [ Rha-ih-,r + R, cos20laiv+ cos i sin 201R/,z I/ + - 7 r {(19),-jkl r2 {RBhaih + RBV COS2 02ai+ cos sin 202RB,,Iz}] r2 for subscript i = x or y, and ah = (1 - cos2c)l/ - sin2ly, a = -( + sin 22)y, a - sin 2 lly, ayh = - sin 21 I + (1 + cos 2o)ly, ayv =- sin 2p1' - (1 - cos 2c)ly. The z component of the field is given by r (20) e-jkr2,0 RBv sin 02- {cos Plx + sin Cly + sin 021z}] /2 In (18) and (19), the Fresnel (Rv,Rh) and Born (BBh,R(,5 ) reflection coefficients are evaluated at the saddle point. Evaluation of the branch cut is far more complex. In this case, the integral is expanded near the branch cut and the contour is deformed to the steepest descent at the branch cut using a change of variables [11] cos(w - 0) = cos(Wb - 0) + js2 After much algebra, the branch cut contribution of (14) is given by Eibr IZ1 /2j sin tVb ra n [U(01 - )fl { Rah + RVcos2wba v+ cos cos 2wbR'l1} + 167r p U(02 - Oc)f2 {R'haih + R (5)COS2baiv+ cos Cos 2wbR(,5) lZ} (21) for the x, and y components, and Eba IZ /sin3/2wb [U(01 - O)f1R'+ U(02 - ) Of2R( 5) { b (s l + sin )-sin w(22) { cos Wb (cOS plx + sin Cpl.) - sin Wblz} 9

for the z coInponent. II (21) and (29). f;= - k ' C'-) anl II ' thl, and. sin(ub-, )r h ' RIh are coefficients of the linear term in the Taylor expansion of the each reflection coefficient at the branch point given by R' 2o R' - 2a, k cos wb ' h cos wb 1' 5)' _ 2 —1 '- (9J +, I )h 2 CO - ) B - -(gu + gv2 R' 2- 1v2) By 21 ki cosu,'b Bh - 2kl cos Wb (ghl + gh2) where ADo 9vi = kl cos Wb AEo gv2 = 1 cos tUb Aa ghl = k1 cos Wb Ac gh2 = k1 cos Wb ( f - 2)C - 2kme-j2km coswb CB K COS Wb "[k 4-KC CB] 2kim ---, K cos Wb Aa LBC -j2k1me-j2kl mcos t-b Cwb A2 - COS Wb A [j2kim 4+C BC] L cos Wb Aa J and 2 A = Re{h(wb)-1}, B=Re{h(wb) 2 k cos wb - } 2 27 C =1 - e-2klmcoswb D=kcS2 Wb + 1 sin2Wb, = - o s2w+ in2 W, h(wb)e-2 c2 wb/2erfc(_- a j sin 2Wb / sin(wb - 0) The above formulation is valid when the observation point 0 is away from the complex critical angle. It can easily be shown that this formulation reduces to the flat boundary case by letting a go to zero. In this case, h(wb) = 1 which reduces A and consequently R(t'5) and R' to zero. The remaining terms in (21) and (22) would be the branch cut contribution from a flat boundary. To extend the valid region to observation point near the critical point, (18) may be modified as Ei -- IS + u(O - O,)IBc * F(3) (23) 10

where F(3) is a correction factor. This function cani be obtained by retaiiiIlg hligllhe order terms of the Taylor series expansion of the integranid near the branch cut. For a dielectric with flat interface. F(3) is represented by the \eber or parabolic cylinllder function as F(/3) = eJ(32/2+3r/8) (232)3/4D._3/24 ( + j/3) where 3 = 29/kir1 sin 0 -'. Here D-3/2(3 + ji3) is Weber function of order -3/2 and is defined as [11, 12]. ( 2e/4 0t 2n-1 -(x+t2)2/2dt As i3 increases, F(3) approaches unity and (23) reduces to (18). It is, however. extremely difficult to modify the formulation for a rough boundary in such simple manner as the higher order derivative of the special functions in the integrand are difficult to evaluate. The asymptotic expressions clearly indicate that the diffracted field due to the branch cut contribution (see(21) and (22)) decays as 1/p2 whereas the saddle point contribution decays exponentially (eJkl/r). Hence for far observation points the lateral wave is the dominant source of the diffracted field. As will be shown later, the lateral wave for rough boundaries are weaker (depending on kou) than the lateral wave for a flat boundary. 4 Numerical Simulation In this section, numerical examples are considered to examine the validity of the solution based on distorted Born approximation and to demonstrate the sensitivity of the field intensity at the observation point to the canopy parameters such as effective dielectric constant, canopy-air interface roughness, and transmitter and receiver positions. In the following simulations, the transmitter dipole is assumed to be operating at f=30MHz having II = 1. To examine the effect of transmitter polarization, a vertical dipole (1 =- ^) and a horizontal dipole (1 = y) are considered in a canopy with e = 1.03+jO.006, and flat interface at h=3m (0.3A). Figure 2(a) and 2(b) show the three components of the diffracted electric field (excluding the direct field contribution from the dipole) as a function of distance for a observation point 3m below the interface (h' = 3m) for the vertical and horizontal dipole respectively. Also shown in these figures are the results obtained from the asymptotic evaluation(see (18-22) for koo = 0). From these figures, it is obvious that the field of a vertical dipole experiences far less propagation loss than the field of a horizontal dipole. 11

Figure:3 showrs the vertical component of the total field (diffracted plus thle direct field) of a xvertical dipole as a function of distance wahen Ih = 2m (0.2A) and I/' = 1ir (0.1A) for three different values of E1. In this simulation, the imaginary part of -- is kept constant and the real part is changed from 1.01 to 1.05. Figure 4 shows a simulation similar to that of Fig.3 with the exception that here the real part of the effective permittivity is kept constant and the imaginary part is increased. The total field calculation for Figs. 3 and 4 is accomplished using a numerical integration for the observation points as far as 2km and for the further point the asymptotic method is used. As expected, the propagation loss increases with increasing the imaginary part of e1. Next we examine the effect of location of transmitter and receiver in the canopy. It should be noted that the dependence of the field on h and h' is of the form. h + h'. Figure 5 show the dependence of the field of a vertical dipole as a function of ko(h + h') for three different permittivity values at fixed lateral distance (p = 2.8km). Unless the observation and source points are very close to the surface, the field decays exponentially as a function of ko(h + h') with an exponential factor proportional to Im[/-t~]. To examine the accuracy of the distorted Born solution, a vegetation layer with smooth boundary is considered (a = 0, and a non-zero m). For this case an exact solution exists. First the Fresnel reflection coefficients are considered. Figure 6(a) and 6(b) compare the magnitude and phase of the phase compensated Fresnel reflection coefficient(Rhe2jklzm) of a dielectric interface with 1 = 1.03 + j0.005 at m = 0.6A with those computed by the approximate distorted Born solution for perpendicular polarization. The computations are performed at 30MHz for two different values of 1 and two different layer thickness values as a function of sin Oi = kp/ko. Similar results for parallel polarization are shown in Fig.7(a) and 7(b). Their errors in magnitude and phase do not exceed 5% and 50 respectively. The accuracy of the distorted Born approximation degrades as the dielectric constant and layer thickness(m) increases. To determine the region of validity of this solution, the exact solution was compared with the distorted Born approximate solution for a wide range of e1 and m. It was found that the accuracy of the approximate solution improves as the imaginary part of the dielectric constant increases. To specialize the region of validity to the problem at hand and lower the number of independent variable, we considered a canopy with vegetation particle permittivity e, = 50 + j25 and used Polder Van Santer dielectric mixing formula [3] to compute the effective dielectric constant e1 from E - 1 e - 1 f = eV + 2e1 3e1 where f is the volume fraction. Tolerating 5% error in magnitude and 50 error in 12

phase of the reflection coefficient, the region under curve showin in Fig.S denotes values of 7n and -1 where the distorted Born approximattion produces valid results for the reflection coefficients. Wtith a confidence in the formulation of the distorted Born approximation, the effect of surface roughness on the reflection coefficient is considered next. Figure 9(a)(perpendicular polarization) and 9(b)(parallel polarization) show the magnitude of reflection coefficients as a function of kl/ko = sin 0i and two different values of surface rms height a = 7,/10 and a = 7r/5. In these simulations, e = 1.03 + j0.015 and rn was chosen to be 2a. It is shown that the reflection coefficient is drastically reduced by the surface roughness for low incident angle and for large value of klp/ko. The reduction in reflection coefficient is less prominent near the critical angle. This property is very important so far as the field of a dipole is concerned since a significant portion of the contribution of the integrand to the integral comes from this point. Figure 10 shows the magnitude of the reflection coefficient for both polarizations versus normalized surface rms height for a canopy with E1 = 1.03 + jO.015 at 30MHz. Figure 11 compares the field of a vertical dipole (h = 0.2A, h' = 0.1A) obtained using exact solution and the approximate distorted Born solution in a canopy with es = 1.01 + j0.006 and two different values, m = A/4 and rn = A/2 at 30MHz where an excellent agreement is shown. A similar comparison is shown in Fig.12 where the effect of imaginary part of the dielectric constant is examined. The effect of surface roughness is shown in Figures 13 and 14. Vertical and horizontal(l = y) dipole in a canopy with el = 1.01 + j0.006 and es = 1.03 + j0.015 and three different values of surface rms height, a = 0, oa = A/2, and a = A are considered (observation point at 4 = 450). It is shown that the surface roughness reduces the field intensity which is a function of the surface rms height. It is very interesting to note that despite relatively significant surface rms height, lateral wave is still the dominant source of the field in a forested area. This can be attributed to the fact that most contribution of the mean field comes from the integrand of the equation(15) for values of k,/ko near the critical angle. As shown before the value of the reflection coefficient near the critical angle does not experience a significant reduction due to the surface rms height. 5 Conclusion The effect of canopy-air interface roughness on the propagation of electomagnetic waves in forested environment was investigated in this paper. An analytic formulations that takes effect of the rough boundary into account, are obtained for both cases of plane wave and an arbitrary oriented dipole excitation. The solution is obtained 13

usilng distorted Born approximation to a volumnetric integral equation for the inIuced( polarization current in a hypothetical layer above the canopy. Thiis form-ulation was validated by comparing the approximate results with the exact results in the special case of smooth interface. It is shown that the canopy-air interface roughness reduces the mean field surface reflectivity drastically for plane wave illumination at incident angles below the critical angle. A significant result of plane wave simulations was the discovery of the fact that the mean surface reflectivitvy near the critical angle is not drastically affected by the surface roughness, thereby allowing the propagation of the lateral wave despite significant dielectric interface roughness. Direct simulations of the field of an arbitrary dipole in a forest with different effective permittivity and surface roughness show that the field at the observation point is anywhere between 0 to 10dB lower than that for a smooth air-canopy interface dependent on the value of the rms height surface roughness(koa). 6 Acknowledgment This work was supported by the U.S. Army Research Office under contacts DAA655 - 98 - 10458 and DAAH04 - 96 - 1 - 0377. 14

References [1] Sarabandi K., Tsenchieh Chiu, "Electromagnetic Scattering from Slightly Rough Surfaces with Inhomogeneous Dielectric Profiles," IEEE Trans. Antennas Prop — agat., vol. 45 pp.1419-1430, Sep. 1997. [2] Theodor Tamir, "On Radio-Wave Propagation in Forest Environments," IEEE Trans. Antennas Propagat., vol. AP-15 pp.806-817, Nov. 1967. [3] Polder, D, and J.H. van Santen, "The Effective Permeability of Mixture of Solids," Physica, 12, pp.257. [4] Foldy,L.L., "The Multiple Scattering of Waves," Phys. Rev., vol. 67, pp.107-119, 1945. [5] Tavakoli, A., K. Sarabandi, and F.T. Ulaby, "Microwave Propagation Constant for a Vegetation Canopy at x-band," Radio Sci., vol. 28, no. 4, pp.549-588, 1993. [6] E. Mougin, "Microwave Coherent Propagation in Cylindrical-Shape Forest Components: Interpretation of Attenuation Observation," IEEE Trans. Geosci. Remote Sensing, vol. 28, no. 3, pp. 315-324, 1990. [7] Le-Wei Li, Tat-Soon Yeo, Pang-Shyan Kooi, Mook-Seng Leong "Radio Wave Propagation Along Mixed Paths Through a Four-Layered Model of Rain Forest: An Analytic Approach," IEEE Trans. Antennas Propagat., vol. 46, NO. 7, pp.1098-1110, Nov. 1998. [8] G.P.S. Cavalcante, and A.J. Giardala, "Optimization of Radio Communication in Media with Three Layers," IEEE Trans. Antennas Propagat., vol. AP-31, pp.141-145, 1983. [9] S.S. Seker "Radio Pulse Transmission along Mixed Paths in a Stratified Forest," IEE Proceeding Microwave & Propagat., vol. 136, pp.13-18, 1989. [10] L.Tsang, J.Kong, and R.T.Shin, Theory of Microwave Remote Sensing. New York: Wiley, 1985. [11] L.B.Felsen, and N.Marcuvitz, Radiation and Scattering of Waves, New Jersey: Prentice-Hall, 1973. [12] L. Brekhovskikh Waves in Layered Media. New York: Academic Press, 1960. 15

[13] Abramowitz. M. and Stegul. I.A. Handbook of mathelmatical functionsr Dover. 1965. 16

z E2 is replaced with E1 plus a polarization current p Figure 1: A short dipole embedded in a forest canopy with rough interface modeled by an effective dielectric constant el. A layer above the canopy is, equivalently, replaced with a dielectric slab with planar interface(z = 0 plane) and permittivity e1 in addition to a polarization current, J = iklYe - E1)Et. in addition to a polarization current, J = Z-kjYj(F2 - E,)Et. 17

a -60 --. -... Z. it -80 -90 -100 -110 -120. 10 104 Distance[m] (a) I=: 30 ___ E Asym 0, E Asym -4 E Asym o -.:.. E xact 0 E Exact::.i:: ' ' ' 1 E Exact - o0 -.........................|..................,. A...... 4,, -90 -100. -110 103 104 Distance[m] (b) =y Figure 2: The diffracted field intensity components (see (15))for a vertical(a) and horizontal(b) dipole in a canopy with e = 1.03 + jO.006 with a smooth boundary (a = 0), h = 3m, h' = 3m at 30 MHz. The results are obtained with the approximation solution obtained from the asymptotic integral evaluation. 18

19n -.... I... F I I r........ I - I. I I I,.. - I. I -1 - -. I I - I I . -1 I I I I. E1 =1.01 + jO.006 0- I I: \. =.... -- - - - =1.03 + jO.006 - E =1.05 + jO.006 -2 0............::..::::::::::: ' i::: \::::::::::::::: -40 --60 -80 -100 -80 -. 2. -::......:..-'.......:.. '. i-........................................... -1 0 101 102 103 104 105 Distance(m) Figure 3: Electric field of a vertical dipole as a function of distance in a canopy with effective permittivity e. Three different values for e. are used, 1.01 + jO.006, 1.03 + jO.006, 1.05 + jO.006, and the location of the dipole and observation points are respectively, h = 2m (0.2A), h' = Im (0.1A). 19

- \!! I I, i i i I a 1i11...::::::::::::: -80......................::::::.:::.:.:-:::-: -:::::\:.:.:-::.::::.....:.::.:..: \. -100 i::i:::::: -120 -...............:::::::::............ -140-....... 10~ 101 102 103 104 105 Distance(m) Figure 4: Electric field of a vertical dipole as a function of distance in a canopy with effective permittivity e1. Three different values for e1 are used, 1.03 + jO.006, 1.03 + jO.03, 1.03 + jO.06, and the location of the dipole and observation points are respectively, h = 2m (0.2A), h' = im (0.1A). 20

-fn( I I I -64 - 6 6......................................,-68 -J. ~ = 1.01+jO.006............. ~1 = 1.03+jO.006 70 - - - = 1.05+jO.006 - 7........................................... -74 --------—........................ ---------------.-74 - 7 6.............................................................................................. -78 0 2 4 6 8 10 12 14 k0 (h+h') Figure 5: Electric field of a vertical dipole(f=30MHz) in a canopy with effective dielectric constant ^El as function of ko(h +h') at a fixed observation point(p = 2.8km). 21

so, - I - I I I I I I IT 60 40 20 0 -20 -40 I - Rej] = 1.03, m=6 / 6 Re{E,J = 1.01. m 10 - e~lr!. 'i -60 II I I, I I I 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 kp/ (a).a (b) Figure 6: Comparison of the magnitude(a) and phase(b) of exact reflection coefficient(Rhe2jklzm) with the approximate distorted Born solution(see (8)) for perpendicular polarization. Two permittivity values and two mean layer thickness values are shown as a function of sin O = -kp / ko. El = 1.03+jO.005, and m = 6m (0.6A) and the other is for e1 = 1.01+jO.005, and m = 10m (A).: frequency is 30MHz. 22

70. - 1 - I I- - - I - r I 60 50 - / 40 - RlE, ]= 1 01. m =10 30 - / Re(1 1.03, m =6 0 1 Exact -10 - - Bom -20 kpko (a) 150s 100F 50 so 7 O -50 -100 -150.Rer,=l].01, m 10.................. /.. i// i: | Exact I --- Bom: i...................................!/ /I if ---. -z t I I.... 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 kko p0 (b) Figure 7: Comparison of the magnitude(a) and phase(b) of exact reflection coefficient(Rve2jk'z) with the approximate distorted Born solution(see (9)) for parallel polarization. Two permittivity values and two mean layer thickness values are shown as a function of sin o = kp / ko. e = 1.03+jO.005, and m = 6m (0.6A) and the other is for el = l.01+jO.005, and m = 10m (A).: frequency is 30MHz. 23

Re[E1] E 0 0 0.5 1 1.5 2 2.5 3 Volume Fraction(%) Figure 8: The region of validity of the distorted Born approximation (points under the curve) is shown with regard to the reflection coefficient for 5% error criterion. 24

(a) - 2 ' '! i \ x. - 2 0..................... -30 - 0 0.2 0.4 0.6 0.8 1 1.2 14 16 1.8 2 k/o (b) Figure 9: The magnitude of reflection coefficient for perpendicular(a) and parallel(b) polarizations as a function of kp/ko = sinOi for a canopy with 1 =- 1.03+jO.015 assuming flat surface (a = 0), a- = A/10 (m = A/5), and cr = A/5 (m = 2A/5). It is shown that surface roughness significantly decrease the magnitude of the reflection coefficient 25

_1 - --.. Re[] = 1.01 -1.......Re[ 1.03.'.-....... -2 --. -3 -. m::,. Re[c]= 1.03 " - 24............................................................ -.:.:. 7 7 0 1 2 3 4 5 6 7 8 9 10 perpendicular and arallel polarization. L — v-pol " 0 1 2 3 4 5 6 7 8 9 10 the critical angle for a canopy with 61 = 1.03 + jO.015, and ~1 = 1.01 + jO.005 for perpendicular and parallel polarization. 26

-20 - -- - -4 0 -................ - -650 -............................,:V..:. -40 -80 -90 IU..:;......... '... I 102 103 104 Distance(m) Figure 11: Field of a vertical dipole in a canopy with 1 = 1.01+j0.006 and smooth interface obtained from the exact solution is compared with distorted Born approximation for two values of m. An operation frequency, 30MHz, and h = 0.2A and h' = 0.1A are assumed. 27

103 104 Distance(m) Figure 12: Field of a vertical dipole in a canopy with1 = 1.03+jO.006 or 1.03+jO.012 and smooth interface obtained from the exact solution is compared with distorted Born approximation for fixed value of kom = ir/5. An operation frequency, 30MHz, and h = 0.2A and h' = 0.1A are assumed. 28

E I T~ r -50 o=0...... S=~. -60 - -. '\.......... \. -70 -:,\.. -80 - 0D * N N",l = 1.01+jO.006,............... ~.... \:: os. -90 -. 1-,.. \: o 1.03+jO.01 5...... - -100 F -110 - IU I. 103 104 Distance(m) Figure 13: Diffracted field of a vertical dipole in a forest canopy (h = 0.2A, h' = 0.1A) for two different values of e1 = 1.01 + jO.006 and e1 = 1.03 + jO.015 and for three different values of surface rms height, c = 0, C = A/2, and a = A. 29

E y _fi - -:.. N:. -- - c= 0 2. ':,.-... -7 0................. -8 0............ I', _ _. '. ' '":'...................................... 0:.:9:0 -1 10.-...,................. "" <.:: ' ' L:,,.I: -110.. -120 103 1.. Distance(m) Figure 14: Diffracted field of a horizontal dipole( = in a forest canopy (h =.2, h' = 0.1A) for two different values of E = 1.01 + j.006 and E1 = 1.03 + j0.015 and for three different values of surface rms height, a = 0, a = A/2, and a = A. for~~~~ 'he ifretvle.o.ufc rshih,a 0,, /,ada. 30

Effect of Canopy-Air Interface Roughness on HF-VHF Wave Propagation in Forest Kamal. Sarabandi. and Il-Suek. Ioh Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan Ann Arbor, Michigan 48109-2122. Submitted to the IEEE AP-S Microwave Conference 1999, Orlando, FL Abstract The problem of wave propagation in forest is revisited. In particular, the effect of the non-planar interface between the air and the canopy on lateral waves is examined. An analytical formulation is obtained for the mean field when both the transmitter and receiver are within the foliage. This formulation is based on distorted Born approximation and is shown that compared to a planar interface, the field of a dipole in a canopy with rough interface is significantly reduced. 1 Introduction A generally adopted model of a forest at HF through VHF, that attributes wave propagation in forest to a lateral wave, was first developed by Tamir [1]. In this approach, the forest is modeled by a homogeneous halfspace dielectric medium with a planar interface. The field of a dipole within this medium was then evaluated at an observation point near the interface using the asymptotic form of the integral involved. This solution shows that the field at the observation point is dominated by the so called lateral wave that travels along the flat canopy-air interface. In reality, however, the interface between forest canopies and air is not flat, hence it is not clear as to what happens to the lateral wave or whether it can even be excited or not. In this paper, the effect of roughness of interface between canopy and air on the wave propagation in forest areas is investigated. An analytical solution is obtained using the volumetric integral equation in conjunction with the distorted Born approximation. 2 Analytical Formulation Geometry of the diffraction problem is shown in Fig.1 where the dipole located at heights h and h' inside a canopy with effective dielectric constant el. The envelope of the canopy-air interface is denoted by d(x, y). The permittivity of the upper medium(air) is denoted by E2. We modify the problem by extending the canopy to z = 0 plane and assume that there exists a volumetric polarization current J = ik1Yl (E2 - E1)E, where Et = Ei + EW + Es. Here E' is the incident field, E is the reflected wave from the planar interface and ES is the scattered field generated by J itself. To the first order in (E2 - el) the Born approximation can be used to the scattered field ES - k02(c2 - l)/ [E(r) E( ES = k -(E2 o),y) G((r, F). [Ei() + Er( )]dx'dy'dz' (1) J-oo J -oo J-d(x,y) whereG is the dyadic Green's function of the half space dielectric medium [2, 3]. In (1) d(x, y) is a two dimensional random process describing the interface between the canopy and air and is assumed to be Gaussian with a mean value of m(a positive number) and standard deviation of r. Distorted Born approximation provides a more accurate solution for Es. In this approximation, a phase correction term is 1

incorporated into the expressions for t he internal field to accouiint for the difference in t ie prol)agatiol coislltailt between the air and the canopy as seen in tlhe following expression of incident anid reflected waves. Ei = P(k )ci[k;x'+kY'+ k' ']. (i( k')m Er p k' ) RCi[k '+k'y' z'] ei(k2 k-' ) where P(kl) is a unit vector (P = v or h) as is defined in [2, 3] and R,,h is the Fresnel reflection coefficient for vertical or horizontal polarization. Since the layer thickness is not uniform, the phase correction terms in (2) are chosen for a layer with a uniform thickness m. Substituting 2-D Fourier transform ofG in (1) and after performing the integration \with respect to z', evaluation of (ES) would require computation of the term like (ei~Sd) which for a Gaussian d are found to be Csm a2s2 02,S,. 2 s2 2* (ed) = 2 [e2 erfc(-j ) + erfc( ) -sm 2*2 2 2 2 2(3) (e- = [e- 2 erfc(-j X ) + e erfc( ) where s can be either si = k'2 + k[, or s2 = k= - kz. The integration with respect to x' and y' can now be carried out analytically which results in 6(kx-k).'6(ky-ky). This in turn simplifies the integration with respect to kx and ky. Thus the final result is readily obtained as, Et = E ( + (Es) (Rref + Rrn)ei[ki +ky-kz (4) where Rref is the Fresnel reflection coefficient for the canopy-air boundary at z = 0. Rref + R^lrn is accurate enough for horizontal polarization, but it cannot accurately predict the Brewster angle for vertical polarization. To rectify this deficiency, higher order solutions must be obtained but it is sufficient to use only -zzv, ) term of the dyadic Green's function. Thus the partial second order solution can be k22 obtained and is given by s(1.5) __k02 (E2 _- 1)2 = E (1.5) -k02 ) E2 G [Z.i (Ei + Er)]dv (5) s2 Jv The ensemble average of Es(l 5) can be obtained in a manner similar to what was used in computation of (Es(1)). After some algebraic manipulations, the reflection coefficients for the mean field are obtained and given by R(1) k02(E2- El) [R2(id 1) 1 -Born(h) = 2klz -1) Rh desm + ((ei2d) - (e-is2d)) }ei2 S2:(1.5) (l2 El, R -i d El 2 k2 Born() - cl z. ((e i2d)-(e -2d)). (E 2 + k Born(v) 2kteEl s2 2 + -(1 - R2 + Rv(eid) - (e-s)). (1kk2 k 2- k2 e2m Sl V2 where kp2 = k 2 + ki 2 Using the above result, the solution for an infinitesimal dipole embedded within the foliage can be obtained by expanding the field of the dipole in terms of a continuous spectrum of plane waves. Assuming that a dipole, whose orientation is denoted by a unit vector 1, is located at ro = -(h + m)i, and using superposition, the mean scattered field can be computed from the coherent sum of all reflected plane waves. That is (EB n) = o k [( i)R + ( hi)R ]e-ikz(-m- h)e i(k +k Y)dki dk' I ( J-oJ -i k1 y B 2

For more accurate comrnputation. the meain field for a dielectric nedium with a rougl surface is calculated by. Eeract () o E.ract(or = 0) + EBorn(U) - EBorn(o = 0) where Eexact(i = 0) is the field of the dipole in the presence of the upper free-space medium with a planar interface at z = -m. 3 Numerical Simulation In order to verify the accuracy of the distorted Born approximation, first, the approximate analytical solution for reflection coefficient of a planar boundary is compared with the exact Fresnel reflection coefficient when a plane wave is incident at the boundary. We consider a canopy with effective permittivity E = 1.03 + i0.001 and the interface between the canopy and air is assumed planar at a distance m = 6[m] from the x-y plane of the reference coordinate system. Figure 2 and 3 show the comparison between the reflection coefficients as predicted by the distorted Born approximation and the exact solution for both parallel and perpendicular polarizations. A very good agreement is obtained for this example. Further sensitivity analysis show that the accuracy of the distorted Born approximation degrades. Similar behavior is obtained when m is kept fixed and dielectric constant of the dielectric layer is increased. Next we considered the field of dipole inside a dielectric layer with E = 1.03 + iO.001 at a depth of 2[m] below the interface. The field is observed at a depth in l[m] below the interface as a function of radial distance. The exact solution is compared with the distorted Born approximation for a chosen value of m = 2[m] at 30[MHz] in Fig.4. An excellent agreement is obtained. Close examination of these results indicates a maximum relative error of 2 % between the two solutions. As for the plane wave illumination simulations, the discrepancy between the distorted Born and the exact solution increases with increasing m. With confidence in the distorted Born solution, the effect of the interface roughness on the field can now be examined. Fig.5 shows the variation of the field as function of radial distance between the transmitter and the receiver (h = 2[m], h' = 1[m]) for two cases. In the first case, el = 1.01 + iO.6 and koa = 3 (roughness parameter) and the second case, Ei = 1.03 + iO.6 and ka = 2. The results are also compared with those had kar = 0 (flat interface). It is shown that these surface roughnesses reduce the field by a factor of 3-5dB. For these simulations, we used a p.d.f. for d of the following form fd(d) = 2 e-d2 /2o2 where d can only assume positive numbers. 4 Conclusions Analytical formulation for the mean-field of a short dipole embedded in a forest is computed. In this formulation, the effect of the roughness of the air-canopy interface is taken into account. Distorted Born approximation is shown to provide a very accurate results for the limiting case when the interface roughness disappears. Simulated results indicate that the roughness of the interface reduces the contribution of the lateral waves significantly. References [1] Theodor Tamir, "On Radio-Wave Propagation in Forest Environments" IEEE Trans. Antennas Propagat., vol. AP-15 pp.806-817, Nov. 1967 [2] Sarabandi K., Tsenchieh C1hiu, "Electromagnetic Scattering from Slightly Rough Surfaces with Inhomogeneous Dielectric Profiles", IEEE Trans. Antennas Propagat., vol. 45 pp.1419-1430, Sep. 1997 [3] L.Tsang, J.Kong, and R.T.Shin, Theory of Microwave Remote Sensing. New York: Wiley, 1985 3

Magnitude ~2 k, m r -?- - -- ----- m *o E. El i i,;h I k h --- t I Receiver I 0 20 40 60 80 Incident Angle [Degrees] 100 Figure 1. Geometry of a dipole with a canopy with rough interface Phase Figure 2. Comparison of the magnitude of reflection coefficients for a planar interface using exact and Born approximation 0n, 101 a) a) c, o> 0 -20 -30 -40 -:n........................... Exact -- Born........................................................................................................................ ~~~~..........................::..........*... -61( I, 20 40 60 80 100 Incident Angle [Degrees] 0 -uv O 500 1000 1500 Distance Figure 3. Comparison of the phase of reflection coefficients for a planar interface using exact and Born approximation Figure 4. Comparison of exact versus approximate Born solution for the field in a dielectric with flat interface -..-.. -... -...... -60 -80 -100 m -120 -140 -160 -180 -Ann.....,.........::......... Flat:...::...-Flat - - Rough -........................ Re[ - 1.03, ka3:.:;.:.::::Rele1] = 1.01, k;o= 2...................N.. -I. -|, —... 102 Distance Figure 5. Field intensity under two different roughness and effective dielectric const. 4

APPENDIX B Wavelet-Based Method of Moments Approach

Numerical Simulation of Scattering from Rough Surfaces:A Wavelet-Based Approach Daniel Zahn, Kamal Sarabandi, Kazem F. Sabet, James Harveyt Radiation Laboratory Department of Electrical Engineering and Computer Science The University of Michigan, Ann Arbor, MI 48109-2122 tArmy Research Office Research Triangle Park, NC 27709-2211 Abstract In this paper a preliminary study is carried out to demonstrate the application of wavelets for improving the computation time and reducing computational memory required for evaluating the statistics of the scattered field from rough surfaces using the method of moments in conjunction with a Monte Carlo simulation. In specific, Haar and the first-order B-spline wavelet basis functions are applied to the MoM formulation of one-dimensional rough surfaces in order to compare the computation time and sparsity for wavelets in the same family but of higher order. Since the scattering coefficient (the second moment of the backscatter field per unit area) is a gentle function of the surface parameters and the radar attributes, it is demonstrated that a relatively high thresholding level can be applied to the impedance matrix which leads to a sparser impedance matrix and faster computation time. It is also shown that applying a high threshold level the coefficients of the high order wavelets would increase out of proportion, however the effect of these current components averages out when computing the scattering coefficients. The resulting sparse impedance matrices are solved efficiently using fast search routines such as the conjugate gradient method. A systematic study is carried out to investigate the effect of different threshold levels on the accuracy versus computing speed criterion. The computed scattering coefficients are compared to previous results computed using a conventional pulse basis function as well as the existing theoretical solutions 1

for rough surfaces. It is shown that %axvelet basis functions provide substantial reductionis in both memory requirements and computation time. I Introduction The problem of electromagnetic scattering from rough surfaces has been the subject of intensive investigation over the past several decades for its application in a number of important remote sensing problems. Radar remote sensing of the oceans, soil moisture, and mine detection using wideband radars are such examples. For these problems, where the rough surface is either the primary target or the clutter, the understanding of interaction of electromagnetic waves with the rough surface is essential for developing inversion or detection algorithms. An exact analytical solution for random rough surfaces does not exist. However, approximate analytical solutions exists for rough surfaces with specific types of surface roughness conditions. For surfaces with small root mean square (rms) height and slope, the small perturbation method (SPM) is the most commonly used formalism. Formulations based on SPM exist for perfectly conducting [10], homogeneous dielectric [4], and inhomogeneous dielectric [12] rough surfaces. Another classical solution which is valid for surfaces with large radii of curvature is based on the tangent plane approximation [15]. The region of validity of these classical approaches are rather limited. In recent years much effort has been devoted to extend the region of validity of these models [3, 9], however, the improved techniques still have the basic limitations of the original models. An alternative approach for evaluating of the scattered field and its statistics for rough surfaces is Monte Carlo simulation. In this approach many sample surfaces with the desired roughness statistics are generated, and then the scattering solution for each sample surface is obtained using a numerical method. Monte Carlo simulation have primarily been considered for evaluating performance of and characterizing the region of validity of approximate analytical models [1, 2, 3, 9]. In general, the limitations of Monte 2

Carlo simulation of scattering from rough surfaces are the computation time and tile required memory as the typical size of the scatterers (sample surfaces) must be chosen to be much larger than the wavelength. Another issue is that the rough surfaces are the targets of infinite extent which must be truncated appropriately before the numerical scattering solution can be obtained. This can be done either using a tapered illumination [9], or padding the sample surfaces with a tapered resistive sheet [11]. It has been shown that with the tapered illumination, larger sample surfaces must be used as a considerable portion of the induced currents on the surface do not contribute significantly to the total scattered field. Application of the tapered resistive sheet is advantageous in that a relatively small portion of the sample surfaces is actually used to suppress the edge currents. This improves the computation time and reduces the required memory. In order to use Monte Carlo simulation for evaluating the scattering statistics of rough surfaces more routinely, computationally more efficient scattering codes must be developed. In this paper the application of wavelets as a basis function for the expansion of induced surface currents is considered. Traditional method of moments (MoM) in conjunction with Galerkin's method would require matrix fill computation time of the order of N2 and matrix inversion computation time of the order of N3 (using Gaussian elimination). It is well known that the solution of linear system of equations can be obtained far more efficiently using search routines, such as the Conjugate Gradient method, if the matrix of the coefficients is a sparse matrix. In MoM, the application of conventional pulse or rooftop basis and testing functions would usually produce full impedance matrices. Although the diagonal elements are usually larger than the rest of the elements, the smaller elements cannot be arbitrarily thresholded without drastically altering the resulting scattering pattern. The success of wavelet expansion function in generating sparse matrices have been demonstrated for many circuits and antenna problems [5, 6, 7, 8]. In the Monte Carlo simulation of scattering from rough surfaces the quantities of interest are the statistical parameters, such as the mean and variance of the scattered field, and therefore it is expected that the overall accuracy be less sensitive to 3

the threshold level. An investigation is conducted on the use of two different types of wavelets with compact support, Haar and B-spline wavelets with edge wavelets and the effect of different threshold level with regard to the overall accuracy and the computation time. The method is applied to one-dimensional perfectly conducting random rough surfaces to demonstrate the improvements achieved. In the Monte Carlo analysis presented here the tapered resistive sheet approach is used to suppress the edge current for plane wave illumination. The numerical results are also compared with the approximate analytical solutions. II Integral Formulation of Scattering Problem In order to characterize both the backscattering coefficient and the bistatic scattering coefficient using the MoM and a Monte Carlo simulation, a large number of random surfaces with known statistical parameters (kfs, k) must be generated. Then it is desired to find the surface current density Je induced by a plane wave from which the scattered field can be computed. For a horizontally polarized plane wave excitation the scattered electric field is given by: ES(P) = - koZo j Je(p)Hl)(ko, P-P')df', (1) here ko is the wave number, ZO is the intrinsic impedance of free space, H') is the zerothorder Hankel function of the first kind, and p and p' are the position vectors of observation and source points, respectively. The sample surface is discretized into sufficiently small cells, as shown in Figure 1, and equation (1) is cast into a matrix equation. The expression for the plane wave propagating along ki = sin Oix - cos OjY is given by, Ei (m,Y ) - eiko(xm sin -ymcosi) (2) E{xm Iym) = (2 \) 4

As mentioned before, because of the singular behavior of the current near the edges of the surface when excited by a horizontally polarized incidence wave, tapered resistive sheets must be added at the edges of the sample surface in order to suppress the edge currents. The induced surface current on a resistive sheet is proportional to the tangential electric field, or mathematically: n x (n x E) = -RJ. (3) where R is the surface resistivity (for a perfect conductor, R = 0). Another boundary condition for resistive sheets mandates continuity of the tangential electric field across the resistive sheet, that is, [n x E]+ = 0. Therefore the electric field integral equation for the surface current is given by E ) = R()Je(()) + J( )) (ko(e) -('))d'. (4) III Tapered Resistive Sheets Tapered resistive sheets, as introduced in section II, are added on the edges of the surface samples to suppress edge effects on the induced surface current and the scattered fields. An optimum tapered function for resistivity (found by trial and error) is given by [11]: 0 ~ Ll < D DR, 27- 2~+D where D is the width of the sample surface and DR is the width of the resistive section. This taper function was found to significantly decrease diffraction due to the edge discontinuity. Figure 2 shows a normalized resistivity profile of the tapered resistive sheets and 5

the placement of the resistive sheets on a surface described by a Gaussian huImp. 1Figures 3a and 3b show the induced current distribution on a flat surface with and without thin tapered resistive sheets, respectively. Expanding the current in terms of the basis functions Je(p) = E aOnn(p) and applying Galerkin's method to equation (4) we have: jfqm(p)E'(P)dt = n=1 S n=l{ janm (P)n( p)R( p)de + 4 + k je jan4m (P))n(p p')H( )(kol5 - p'l)d'dI},(6) where q)m(p) is the testing function, and qn(p') is the basis function. Equation (6) can be solved using numerical integration, and cast into a matrix equation, given by [Z] [I] = [V]. This matrix equation can easily be solved to find the surface current Je. IV Applications of Compact-Support Wavelets As stated above, the EFIE used in conjunction with Galerkin's method can be cast into a matrix equation using appropriate expansion and testing functions. The restriction on the expansion and testing functions is that they have to be in the domain of the integral equation operation. To expand the induced current in terms of a multiresolution expansion, first, the current must be projected onto the x-axis or the surface must be arclength parameterized. For natural rough surfaces with moderate rms slope it is more convenient to project the current on the x-axis since the domain will be identical for all sample surfaces. Applying a multiresolution expansion to the projected current distribution (f(x)) on the x-axis we have: mh-1 f(x) = ECmhnmhn(X) = ECmenqmn(X) + E E dm,n)m,n(X) (7) n n m=me n 6

which consists of the scaling functions at the lowest resolution, plus wavelets at the lowest resolution and subsequent higher resolutions. This expansion is equivalent to one consisting of only scaling functions at the highest resolution (4mhn(x)). Using the Fast Wavelet Transform(FWT)[18], the program needs only to compute the EFIE impedance matrix at the highest level of resolution. From this computed impedance matrix at the highest resolution, lower resolutions can be found in terms of the original computed impedance matrix. This reduces computation time drastically from having to perform the integration for each wavelet at various resolution levels. Sparse matrices arise due to the fact that the scaling functions and wavelets are orthogonal and wavelets have zero moments. For example, Haar and linear B-spline wavelets both have zero mean and the first moment of the linear B-spline wavelet is also zero. Higher order wavelets also have vanishing higher moments, that is J xO(x) dx = 0 n E {O....N}, (8) where N depends on the order of the B-spline wavelet. Figures 4 and 5 show the Haar and linear B-spline scaling functions and wavelets. The significance of vanishing moments of the wavelets stems from the fact that the integration of the kernel over the domain of the wavelet would produce a small quantity. This quantity is smaller for wavelets with higher vanishing moments. Then, when integrating over the kernel, very small matrix element values are usually calculated for cells that are relatively far from one another with the kernel being fairly regular. For adjacent cells and self-cells the integration of the kernel produce element values that are significant. In fact, the self-cell contribution is the largest element of the impedance matrix. Once the impedance matrix for the highest resolution is computed, the FWT algorithm is applied to the matrix to find an equivalent impedance matrix for the multiresolution expansion. Using the FWT it is expected that a sparse matrix will be generated. 7

At this stage a user-defined threshold level. usually of the order of.01'7 to 1'7, is inplosed on the matrix and only significant elements of the matrix will be preserved. For a preliminary examination of the wavelet-based Method of Moment, a Gaussian hump described by equation y(x) = A -( ~o)2 (9) is used. The effects of the multiresolution expansion and application of different threshold levels on the moment matrix on the bistatic scattering is investigated. Initially, pulse basis functions with the highest resolution is used to generate a reference solution to the problem of a plane wave incident at normal incidence upon a Gaussian hump surface and the solution is found with no threshold applied. This solution became the standard of comparison for subsequent tests on the Gaussian hump surface involving both Haar wavelets and first order B-spline wavelets for various levels of resolution and threshold levels. Using Haar wavelets, described by equation h(X) ~2 (10) along with pulse basis functions, the moment matrix produced had a maximum sparsity level of over 89% when a threshold level of 0.3% was applied and 5 levels of resolution were used. The bistatic scattering pattern produced from the reduced moment matrix shows no significant differences in the scattering levels when compared with the results obtained from the original full moment matrix. The sparsity level increases with the number of resolution levels, because a larger number of localized wavelets contribute to the cancellation effect at slowly varying regions of the induced current distribution. Sharply varying components of the current which are mostly localized are captured by a small number of wavelets. Keeping the same threshold level for the moment matrix, the 8

sparsity achieved for 2. 3. and 4 levels of resolution are respectively 72.5%(. 80.7T 7. and 84.8%. Next, rooftop basis functions along with linear B-spline wavelets wrere investigated. The results from the linear B-spline wavelet case agreed well with the Haar case. and a sparsity of over 99% was observed when a threshold level of 3% was applied and 5 levels of resolution were used, and the overall bistatic scattering pattern shows no significant differences. As expected, higher sparsity is achieved with linear B-spline wavelets. For 2. 3, and 4 levels of resolution with the same threshold level applied to the moment matrix, sparsities of 93.2%, 97.3%, and 97.8%, respectively, are achieved for the linear B-spline wavelet expansion. Once a sparse matrix is obtained, a search routine, such as conjugate gradient method, can be used to find the solution. Since the number of non-zero matrix elements is far less than the original matrix, the computation time for obtaining the solution is reduced drastically. The comparison between the full matrix solution, Haar, and linear B-spline wavelets with 5 levels of resolution and the threshold levels stated above for the respective basis functions is shown in Figure 6. V Monte Carlo Simulations of Random Rough Surfaces Encouraged by results obtained for a simple Gaussian hump surface, the multiresolution expansion method is then applied to random rough surfaces of known statistical parameters. Near normal incidence, sample surfaces for numerical analysis must be at least 40 correlation lengths (f) long in order to accurately characterize the bistatic pattern [13]. This requirement becomes more stringent for near grazing incidence. Since the Monte Carlo simulation of the scattering problem requires the numerical solution of the problem many times, the computation speed achieved by thresholding the moment matrix of a multi-resolution expansion in conjunction with a search routine 9

linear system solver becomes very significant with regard to thie overall computation timlle necessary for evaluating the statistics of the scattered field. V.1 Random Surface Generation Monte Carlo simulations require a large number of sample surfaces of a random process with prescribed surface height statistics. To generate the sample surfaces, the procedure in [10, 15] is used. First, a long string of numbers is generated using a random number generator having the same pdf as the height distribution of the surface (for example, a zero-mean Gaussian pdf). Next, a subset of the numbers of the string are correlated with a weight vector related to the Fourier transform of the desired autocorrelation function [9]. For the simulations presented in this paper, surfaces with Gaussian correlation function and Gaussian height distribution are considered. Hence, the surface statistics are uniquely specified by the surface height standard deviation (rms height), s, and by the surface correlation length, e. Figure 7 shows a sample surface of a Gaussian process with ks = 0.3, k1 = 3.0 and the corresponding histogram of height and calculated correlation function generated from 60 independent sample surfaces. The calculated correlation function gives ke = 3.17, and the calculated standard deviation of height distribution gives ks = 0.3020. These agree closely with the desired surface roughness parameters. V.2 Validation and Results Numerical simulation of rough surface scattering is performed for three different surfaces denoted by Si, S2, and S3. The roughness parameters (ks, kt) for each of these surfaces are respectively ksl = 0.3, ktl = 3.0; ks2 = 0.5, k~2 = 6.13; and ks3 = 2.0, ke3 = 2.5. The first two surfaces fall within the region of validity of small perturbation and physical optics models respectively, and hence the two numerical results can be compared with the analytical models. Comparisons are also made on the threshold level imposed 10

on the moment matrix for both Haar and B-spline xwavelet basis functions and on the matrix solving times using a fast Conjugate Gradient solver for sparse matrices. Another test that was run was the effect of the MoNM and wavelet technique on backscattering enhancement. Finally, the effect of the number of resolution levels on the scattering pattern is investigated. SPM is known to be valid when ks < 0.3, kf < 3.0, and m < 0.3, where m is the rms slope and is given by m = V/s/t for a surface with a Gaussian correlation function. The analytical bistatic pattern for the SPM is derived in [17] for a 3-dimensional surface with a Gaussian correlation function. For a 2-dimensional surface, the bistatic scattering coefficient (oad) is related to the 3-dimensional bistatic scattering coefficient (a'd) via 2d = k 3d Then, '2d is given by: a2d(90, 0s) = 4k3cos2(90)cos2(0i)fhh W(Il - kli) (11) where W(k) is the power spectral density. For a Gaussian surface correlation, W(j k1 - kcj) is given by W(Ikl - kl)= /es2e -(4 [sin()-sin()]) (12) W(\k. - k1\) = \/fs' e(12) Furthermore, for a perfect electrical conductor (PEC), fhh = cos2(4, - /) which equals 1 for a 2-dimensional problem. The incoherent bistatic scattering coefficient for 61 = 30~ using the Monte Carlo simulation with 40 independent samples and threshold applied is compared to the analytical bistatic SPM from equation (11) in Figure 8, and a good agreement is observed. The Physical optics (PO) method using the tangent-plane technique for approximating the fields on a surface, S, is next investigated. Under the tangent-plane technique, the fields present at any point, P, on S are approximated by the fields that would be present on a plane tangent to P. This is a valid approximation if every point on S has a large radius of curvature. The three-dimensional bistatic scattering coefficient is derived 11

in [19] for the PO approximation. Using this formulation a' is calculated and is given by: o s k (1 +cos(Oi)cos(Os) -sin(O9)sin(Os) 2 012d( ----J (13) ( =cos(0,) + cos(0s) ) where -X -e J = 2 dx cos [kx(sin(O0) - sin(0,))] e L J - X-X (14) where X = ks[cos(Oi) + cos(0,)]. For surface S2 which falls into the PO region of validity, a comparison is made in Figure 9 between the PO model in equation (13) and results obtained using the wavelet based moment method. The bistatic scattering results from S2 agree well with the theoretical Physical Optics solution. A comparison of different threshold levels imposed on the moment matrix for surface Si using B-spline basis functions is shown in Figure 10. The parameters for the simulation on Si are 0, = 30.0~, length = 32A, resistive tapered ends length = IA and Ax = 0.1A. For this simulation, only a single surface (N = 1) is used for comparison. As is shown in Figure 10, the scattering pattern varies very slightly and only at angles near grazing observation. This figure indicates that a sparsity of more than 90% can be achieved without substantial compromise in the accuracy of the bistatic pattern using B-spline wavelets. The bistatic scattering coefficient obtained using the Haar and B-spline wavelet of surface S3 are shown in Figure 11. For this simulation, ks = 2.0, kf = 2.5, number of independent surfaces N = 40, sample length = 32A, resistive tapered ends length = 1A, Ax = 0.04A and Oi = 30.0~. Since this surface does not fall into the region of validity for any analytical models, no comparison may be made. This figure shows again that by thresholding the moment matrix, a high sparsity can be achieved without compromising 12

the accuracy. It is also shown that a higher sparsitv is achieved using the linear Bspline wavelet. An average sparsity of 97.3%7 is achieved by the linear B-spline wavelet whereas using Haar wavelets only S2.2% average sparsity is obtained for a threshold level of 0.15%. Figure 12 shows a comparison of the exact current distribution and one where the impedance matrix has a sparsity of about 97%o. A portion of a surface with characteristics given above (ks = 2.0, kf = 2.5, etc.) was used and the sparse matrix has a threshold of 0.5%. Even though this sample surface (which is indicitive of the current distribution for every sample surface) has a current distribution that varies significantly from the exact solution, because the far-field scattering pattern is an averaging process the bistatic scattering pattern does not vary significantly from Figure 11. The threshold level imposed on the moment matrix and the number of multiplications needed to solve the matrix using a fast Conjugate Gradient solver routine is studied next. As stated previously, by imposing a threshold level the scattering pattern remains relatively unchanged, yet the moment matrix could be made quite sparse. It was also found that linear B-spline wavelets would produce a sparser matrix for a given threshold level before the scattering pattern were to deviate significantly from the exact solution. This is due to the fact that the linear B-spline wavelet has a vanishing first moment. Table 1 provides the average number of multiplications necessary to solve for the surface current of S3 for both Haar and linear B-spline wavelets and for different values of threshold levels. It is found that the number of multiplications, or equivalently the computation time, decreases significantly with the first order B-spline scaling function and its corresponding wavelets. For the Haar wavelet-based MoM, a slight improvement was observed (approximately 20%), yet for the linear B-spline wavelet-based MoM, a factor of about 20 improvement was observed. The question arises if the solution produced using a wavelet based MoM technique show the effect of backscattering enhancement. The backscattering enhancement is produces primarily due to multi-path and surface-wave effects on a rough surface. A surface was chosen that has physical parameters where backscattering enhancement is known 13

to occur. and these parameters are ks = 10.636. kC = 19.472. Figure 13 was produces using the aforementioned physical parameters, as well as number of independent surfaces AN = 400. sample length = 30A, resistive tapered ends length = 1A, Ax = 0.0625A and 0i = 10.00. Also, 5 levels of resolution were used, and the matrices had an average sparsitv of 97.5%. It can be seen that around 0s = -10~ there is an enhanced backscattering effect. This figure is in linear scale to display the enhanced backscattering effect, and agrees very well with Figure 5 of reference [20]. The effect of the number of resolution levels on the scattering pattern is next investigated using surface S3. Since it has already been determined that both Haar and B-spline wavelet-based MoM produce similar results, the effect of the number of resolution levels is demonstrated with Haar wavelets only. Monte Carlo simulations are run on random rough surfaces for 5, 4, and 2 levels of resolution, for threshold levels of 0.1% and 0.3%, except for the Haar with 5-levels of resolution in which 0.1% and 0.15% threshold level is used. This is because 0.3% threshold level is too high for 5 levels of resolution and the conjugate gradient solver does not converge. The results based on the full matrix is also included at 5 levels of resolution for comparison. As is shown in Figure 14 the scattering pattern starts deviating from the exact pattern for 5 levels of resolution and 0.1% threshold level imposed on the moment matrix. The average sparsity, obtained from 40 independent samples, for each case in Figure 14 is summarized in Table 2. The scattering pattern agrees quite well with the exact solution for almost all levels of resolutions shown with only slight deviation at near grazing observations. One notable exception is for 4 levels of resolution and an imposed tolerance level of 0.3%, which deviates noticeably from the exact bistatic pattern at certain angles of observation. As shown in Table 2, the matrix is less sparse for fewer levels of resolution at a single stated threshold level. Therefore, by increasing the number of resolution levels, while holding a constant threshold level, the sparsity of the matrix will increase, as is shown in Table 2. 14

VI Conclusion It has been found that using wavelet basis functions with lMoMN and Galerkin s method along with a fast solver routine such as conjugate gradient can drastically reduce both the memory requirements of a system and the time necessary to solve the MoMi matrix. This leads to solutions for rough surface scattering that are quite accurate when compared to other basis functions, yet take significantly less memory and time to solve. Matrices can be made over 97% sparse yet still produce accurate bistatic scattering coefficients in scattering problems. Thus, it becomes possible to generate statistics for the scattering from surfaces of different roughness in a relatively short period of time. The number of resolution levels were shown to play a significant role in determining the sparsity of the matrix and the accuracy of the solution. It was shown that the higher the number of resolutions, the more sparse the matrix could be made without compromising the bistatic scattering pattern. Also, the higher order the B-spline was made, the higher the sparsity achieved in the matrix, and thus, the faster the computation time to solve the matrix. Acknowledgments: This study was supported jointly by the Jet Propulsion Laboratory (JPL) under grant JPL-958749 and by the U.S. Army Research Office under contract DAAH 04-96-0377. References [1] Thorsos, E.I., "The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum", J. Acoust. Soc. Am., vol. 83, pp. 78-92, Jan. 1988. [2] Thorsos, E.I., and D.R. Jackson, " The validity of the perturbation approximation for rough surface scattering using a Gaussian roughness spectrum", J. Acoust. Soc. Am., 86, pp. 261-276, July 1989. 15

[3] Rodriguez, E., Y. Kim. and S.L. Durden. "A numerical assessment of rough surface scattering theories: horizontal polarization". Radio Sci., vol. 27,pp. 497-513, JulyAug. 1992. [4] Nieto-Vesperinas, M. (1982), Depolarization of electromagnetic waves scattered from slightly rough random surfaces: A study by means of the extinction theorem, J. Opt. Soc. Am., 72, 539-547. [5] B.Z. Steinberg and Y. Leviatan, "On the use of wavelet expansions in the method of moments" IEEE Trans. Antennas Propagat.,vol. 41, pp. 610-619, May 1993. [6] K.F. Sabet and P.B. Katehi, "Analysis of integrated millimeter-wave and submillimeter-wave waveguides using orthonormal wavelet expansions," IEEE Trans. Microwave Theory Tech., vol. MTT-42, Dec. 1994. [7] K.F. Sabet, Novel efficient integral-based techniques for characterization of planar microwave structures, Ph.D. dissertation, The University of Michigan, Apr. 1995. [8] K.F. Sabet and P.B. Katehi, "A study of dielectric resonators based on twodimensional fast wavelet algorithm," IEEE Microwave Guided Wave Lett., vol. 6, pp. 19-21, Jan. 1996. [9] Axline, R.M. and A.K. Fung, "Numerical computation of scattering from a perfectly conducting random surface," IEEE Trans. Antennas Propagat., Vol. AP-26, pp. 482 -488, May 1978. [10] S.O. Rice, "Reflection of electromagnetic waves by slightly rough surfaces," Commun. Pure Appl. Math., vol. 4, pp. 351-378, 1951. [11] Y. Oh, K. Sarabandi "An Improved Numerical Simulation of Electromagnetic Scattering from Perfectly Conducting Random Surfaces," IEE Proceedings on Microwave, Antenna and Propagation, To Be Published. 16

[12] 1K. Sarabandi, Y. Oh., F.W. Ulabv "A Numerical Simulation of Scattering from OneDimensional Inhomogeneous Dielectric Random Surfaces". IEEE Trans. oni Geoscience and Remote Sensing, vol. 34, no. 2, pp. 425-432. March 1996. [13] Sarabandi, K., and Y. Oh, "Effect of antenna footprint on the statistics of radar backscattering,"Proc. IEEE Trans. Geosci. Remote Sensing Symp., Firenze, Italy, July 1995. [14] T.B.A. Senior "Scattering By Resistive Strips," Radio Sci., vol. 14, pp. 911-924, 1979. [15] P. Beckmann and A. Spizzichino, The Scattering of Electromagnetic Waves from Rough Surfaces. Norwood, MA: Artech House, 1987. [16] Oh, Y., K. Sarabandi, F.T. Ulaby, "Re-examination of the Kirchoff approximation for scattering from a rough surface", 1993 URSI Radio Science Meeting, Digest p. 406, Ann Arbor, MI, June28-July2, 1993. [17] Tsang, L., J.A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing. John Wiley & Sons, New York, 1985. [18] Rioul, 0. and P. Duhamel, "Fast Algorithms for Discrete and Continuous Wavelet Transforms", IEEE Trans. on Info. Theory, Vol. 38, No. 2, March 1992. [19] Ishimaru, Akira, Wave Propagation and Scattering in Random Media: Volume 2., Acedemic Press, New York, 1978 [20] Li, L., C.H. Chan, L. Tsang, K. Pak, P. Phu, S.H. Lou "Monte Carlo simulations and backscattering enhancement of random metallic rough surfaces at optical frequencies," Journal of Electromagnetic Waves and Apps., vol. 8, no. 3, pp. 277-293, 1994 17

[21] Phu, Phillip, A. Ishimaru, dV. IKuga. -Controlled millimeter-wave experiments and numerical simulations on the enhanced backscattering from one-dimensional very rough surfaces" Radio Science, vol. 28. no. 4, pp. 533-548, July-August 1993 18

V 0 1 2 3 4 5 Position (in X) Figure 1: A typical discretized sample surface. Tonnrn,1 D Ck-t-, #o O WI adpert s nlbiilv o inee t1 Fx 20.8 -., 0.6-x x 0.4 5 Surface 2 0.8-,0.6 -' 0.4 - c, 0.2 03 0 1 2 3 4 5 Figure 2: Normalized resistivity profile and tapered resistive sheets placed on a surface described by a Gaussian hump 19

Surface current - tapered resistive sheet 2.5 -2 7 --- —--- l 0.5 -0 5 10 15 20 25 30 35 Surface current - without tapered resistive sheet 5 4 3 2 1 — 0 5 10 15 20 25 30 35 Figure 3: Surface currents on a 32-A flat plate when a)resistive tapered sheets are used and b) resistive tapered sheets are not used. The surface (not shown) is a PEC flat sheet with normal incidence. 1.5 - - 1.5 1 1 0.5 0.5 -0.5 -0.5 -1 -1 -1.5 - -1.5 -1 0 1 2 -1 0 1 2 Figure 4: Pulse basis function and the corresponding Haar wavelet 20

~ E!I 1 1 ' 0.5 0.5 X o -0.5 - -0.5 -1 -1 -2 0 2 -2 0 2 Figure 5: Linear basis function and its corresponding wavelet 25 Exact Solution 20 - - Haar Solution -- - B-spline Solution 15 -c -10 01 5 - C Angle (in Degrees) Figure 6: Comparison between the exact solution, Haar wavelets with 5 levels of resrs-5 olution, and B-spline wavelets with 5 levels of resolution. Normal incidence is used on a surface described by a Gaussian hump. Threshold levels are 0.3% and 3.0% for Haar and B-spline moment matrix with sparsity levels of 89.6% and 99.1% respectively. 21

Mean surface height (a) /' p(z) IS.......... P(z) -I e (b) f (Ic) Figure 7: (a) Generated random surface, with mean surface height shown (b) Theoretical and Monte Carlo probability density function of the surface height distribution, ks (c) Theoretical and Monte Carlo correlation function -TOO -50 0 50 100 Angle (in Degrees) Figure 8: Comparison of SPM with Monte Carlo simulation. Generated surface parameters are: ks = 0.3, ki = 3.0, N = 40, length = 32A, resistive tapered ends length = 1A, Ax = 0.1A, and 0- = 30.0~. 99

-80 -60 -40 -20 0 20 40 60 80 Angle (in Degrees) Figure 9: Comparison of the Physical Optics model and a random surface with ks = 0.5, k = 6.13, N = 50, length = 100A, resistive tapered ends length = IA, Ax = 0.2A and 0i = 30.0~. 5 * - Threshold = 0.0%, Sparsity = 0.0% -80 -0 40 -20 0 20 40 60 80 Angle (in Degrees) Figure 10: Bistatic Scattering from a single random surface with ks = 0.3, k = 3.0, N = 1, length = 32A, resistive tapered ends length = 1A, Ax = 0.1A and 0i = 30.0~. The various threshold levels imposed on the moment matrix and its corresponding sparsity level are given in the figure. 23

10r C 5 CO o Cr) C, m —, B-spline, 0.0% Tol., 0.0% Sparse ~ - B-spline, 0.1% Tol., 93.4% Spars 0 * * B-spline, 0.5% Tol., 97.3% Sparse 8- - - Haar, 0.0% Tol., 0.0% Sparse.. - Haar, 0.15% Tol., 82.2% Sparse ~^0 -50 0 50 100 Angle (in Degrees) Figure 11: Comparison of threshold levels on scattering pattern for a random surface with ks = 2.0, ki = 2.5, N = 40, length = 32A, resistive tapered ends length = 1A, Ax = 0.04A, 9? = 30.0~, and 5 levels of resolution. 4.5 -- 0.0% Thresholding | 4 — 4 — L 0.5% Thresholding \ 3.5 i 3 10 12 14 16 18 20 Position (in Dr) 1.5of 0.5%). A portions of a random surface with ks = 2.0, k 2.5, length 32A, resistive 0 12 14 16 1 20 tapered ends length = 1X. Ax = 0.04A, Oi = 30.00, and 5 levels of resolution is plotted. 24

-80 -60 -40 -20 0 20 40 60 80 Angle (in Degrees) Figure 13: Test for backscatter enhancement effect. The surfaces under test have ks = 10.636, k~ = 19.472, N = 400, length = 30A, resistive tapered ends length = 1A, Ax = 0.0625A and 0. = 10.00. As can be seen, a significant backscatter enhancement is shown at 0s = -10~. 10. 1 5 W0 0 -50 0 50 100 () -5 /M Haar, 0.0% Tol., 5 Levels / -.. Haar, 0.15% Tol., 5 Levels -- Haar, 0.3% Tol., 4 Levels 0 Haar, 0.3% Tol., 2 Levels -5 0 50 100 Angle (in Degrees) Figure 14: Comparison of resolution levels used on scattering pattern for a random surface with ks = 2.0, ke == 2.5, N = 40, length = 32A, resistive tapered ends length = 1A, Ax = 0.04A and Oi = 30.0~. 25

Threshold Sparsity Number of Mults Bspline 0.0(%7 0.0'(% 5.3 x 108 Bspline 0.05%7 90.7% 8.4 x 10' Bspline 0.1% 93.4%7. 6.4 x 10' Bspline 0.5% 97.3% 2.6 x 10' Haar 0.0% 0.0%' 7.0 x 108 Haar 0.05% 66.1%0 9.5 x 108 Haar 0.1% 76.9% 7.0 x 108 Haar 0.15% 82.2% 5.6 x 108 Table 1: Number of multiplications needed to solve the matrix and sparsity vs. imposed threshold level for a random surface with ks = 2.0, ke = 2.5, N = 40, length = 32A, resistive tapered ends length = 1A, Ax = 0.04A, 0I = 30.0~, and 5 levels of resolution. Threshold Sparsity Haar 0.0% 0.0% Haar-5 levels res. 0.1% 76.9% Haar-5 levels res. 0.15% 82.2% Haar-4 levels res. 0.1% 68.9% Haar-4 levels res. 0.3% 85.5% Haar-2 levels res. 0.1% 51.1% Haar-2 levels res. 0.3% 74.4% Table 2: Sparsity vs. threshold for a number of different resolution levels used for a random surface with ks = 2.0, kf = 2.5, N = 40, length = 32A, resistive tapered ends length = 1A, Ax = 0.04A and Oi = 30.0~. 26

APPENDIX C Physical Optics Iterative Approach

Submitted to the IEEE AP-S Microwave Conference 1998. Atlanta. (GA NUMERICAL SIMULATION OF SCATTERING FROM ROUGH SURFACES USING AN ITERATIVE PHYSICAL OPTICS APPROACH Kamal Sarabandi and Daniel Zahn Radiation Laboratory, Department of Electrical Engineering and Computer Science University of Michigan, Ann Arbor, MI 48109-2122 Abstract The application of iterative Physical Optics (PO) in conjunction with a Monte Carlo simulation for characterizing the bistatic scattering coefficient of random rough surfaces is examined in this paper. The iterative PO method offers decreased memory and computation time restrictions compared to the standard numerical methods such as the Method of Moments (MoM). Results from the iterative PO method are compared to the standard electric field integral equation (EFIE), the magnetic field integral equation (MFIE) as well as the existing theoretical solutions for rough surfaces. It is demonstrated that memory requirements and computation time is significantly decreased while providing fairly accurate results for surfaces with moderate to low rms slope. I Introduction Development of numerically efficient Monte Carlo-based models for simulations of electromagnetic scattering from random surfaces has attained significant prominence over the past decade [4, 2, 1]. A major stumbling block in this endeavor has been the large memory and computation time requirement. This is because the size of the scatter is large compared to the wavelength and the Monte Carlo simulations demand computation of the scattering problem many times. Iterative methods offer an alternative approach when exact solutions are not available and have been used in different electromagnetic problems. By construct evaluation of iterative solutions are rather straight forward especially when the perturbation parameter is relatively small. Physical Optics (PO) approximation is known to provide accurate approximation for the induced surface currents provided that the local radii of curvature at each point on the surface of scatterer is large and the surface is convex. For concave surfaces and surfaces with many adjacent humps multiple scattering drastically alter the standard PO current. However these surface current variations can be estimated through an iterative process. Unlike Method of Moments (MoM) which requires matrices on the order of N2 to find the surface current of the sample surface with N elements, the iterative Physical Optics method only requires memory size of the order N. Thus, substantial memory savings are realized. Also, since no solver routine is necessary in order to solve for the surface currents, as in the MoM, substantial time savings are realized as well. We investigated the use of the iterative Physical Optics method upon a variety of surfaces in order to find the approximate region of validity for such a method. The results obtained using the iterative PO method are also compared to the results found using the electric field integral equation (EFIE) with tapered resistive sheets at the ends of the surface samples as well as the magnetic field integral equation (MFIE) for a horizontally polarized wave. II Formulation Using a Monte Carlo simulation, first many sample surfaces of a stochastic random process representing the desired random rough surface is generated using a standard procedure [4, 2]. As mentioned earlier it is far more convenient to solve for the induced polarization current using an iterative method as opposed to brute force numerical methods. Magnetic field integral equation (MFIE) can be used as the basis for the iterative PO solution using the surface curvature as the perturbation parameters. The MFIE for the induced surface current (J,) over a perfect electric conductor is given by [6] 1 ( 1 eik(r-r') J,(r) = 2 (n x Hi) + 2 fJs)(n. (r'-r)) ik - jr - 1 ds' (1) where Hi is the incident magnetic field, n is the unit normal at the observation point, and f represents the principle value integral. It is obvious from (1) that for a flat surface (r' - r) - n = 0 which renders 1

J,(r) = 2(n x Hi) which is the standard Physical Optics approxiimation for the induced surface current. For undulating surfaces with large radii of curvature tlhe contribution from the integral has a secondary effect and therefore (1) may be solved in an in iterative fashion. To examine the accuracy of the iterative PO approacll tile results must be compared with an exact solution obtained from a numerical method. Numerical solutions for one-dimensional rough surfaces are tractable and thus the accuracy of the iterative PO approach is examined for one-dimensional rough surfaces. The MFIE for two-dimensional problems (one-dimensional roughness) can easily be obtained. For 2-D scattering problems the transverse (Jt) and longitudinal (J ) components of the induced current are decoupled and the integral equations for a TM and TE incidence are respectively given by:,ik n.(p - p') J,(p) = 2(n x H'). Z + 2 fJ (p/)Hl)(kolp- p'l ) d' (2) ik 0W - P).(1 (p' - p) d~' (3) -J(p) = 2H + 2 fJt(')Hi (koIp- p'l) d (3) where t = z x in is the tangent unit vector. Again we recognize that for a flat surface (n = h ) the contribution from the integrals in (2) and (3) are zero and the PO currents are retrieved. Consider a perfectly conducting rough surface illuminated by an E-polarized (TM) plane wave. Points on the surface can be grouped into two categories: 1) lit points and 2) shadowed points. As a first-order approximation the induced current over the lit region is given by J(1) = 2(i x H') and over the shadow regions J(1) = 0. This current produces a scattered magnetic field whose tangent on the surface induces a secondary PO current. It can be shown that the expression for this first order scattered field is given by x Hi) = J(')(p')H(')(lko) p - p- l ) n (- de' (4) (1)- 4 J1 IP - P'I and hence the second order PO current can be obtained from J(2) = 2(n x Hl)) which is exactly the same as the integral in (2). From this argument it is evident that starting with J1), the integral in (2) can be used in a recursive manner to find a higher order solution. There is a subtlety in the computation of the second order current over the shadow regions. The reason for shadow is that the scattered field is out of phase with the incident wave. Basically we have to add the contribution from the incident field to the first-order scattered fields over the shadow regions. Therefore the corrected second order solution is J = J(2) + 2 (n x Hi) over shadow region (5) Once J(2)corrected is characterized, higher order currents can be obtained from (n) = io J(n ) (p) H (l)(ko lp-p (I d-' (6) I2 p-p' and Jz = En Jj(n). The convergence is examined by monitoring the normalized difference in the successive solutions (| j)III MFIE and the Method of Moments Using Galerkin's Method, equation (2) is re-written as O/im2[n x H'(p)]. = Om JZ(P)- n Jz )H (ko Ip - p-l)" - (P ). (7) n=l 'P - P'l n=' 2

which can be solved using numerical integration, and cast into a matrix equation. In equation (7). o,, anll O,, are the basis function and testing function respectively. This is also compared with the EFIE which, using Galerkin's nlethod is given by m E ( [p)d = f,,n R(p)Je(p) + 4 m j OnkJe(p)H((kop -p )d'd (8) where R is the surface resistivity (for a perfect conductor, R = 0). For the EFIE, tapered resistive sheets are used to reduce the edge effects. An optimum tapered function is reported in [2]. IV Results To verify the accuracy of the iterative PO method, a perfectly conducting surface with two Gaussian humps, as shown in Figure 1(a), was studied. The total bistatic scattering pattern for normal incidence excitation is shown in Figure l(b) and excellent agreement is observed between the electric field integral equation (EFIE), magnetic field integral equation (MFIE) and the iterative Physical Optics (PO) method. After verifying that the results from the iterative PO method, EFIE and MFIE gave similar results, numerical Monte Carlo simulations of rough surface scattering is performed for four differ ent surfaces denoted by S1, S2, S3, and S4. The roughness parameters (ks, ki) for each of these surfaces, created using [4], are respectively kso = 0.3, kt1 = 3.0; ks2 = 0.5,k(2 = 6.13; ks3 = 1.0,kf3 = 3.0; ks4 = 0.5,k 4 = 2.0. Surfaces S1 and 5S2 were chosen because they fall into the regions of validity for the Small Perturbation Method (SPM) and the Physical Optics (PO) method, respectively. For surfaces Si and S2, whose bistatic scattering patterns (a0) are shown in Figures 2(a) and 2(b) respectively, good agreement is shown between the EFIE, MFIE, iterative PO method and the analytical solutions. Si sample surfaces are 30 A in length (except for the EFIE where 1 A resistive tapered sheets are added on each end) with Ax = 0.2A and Oi = 30.0~. 40 sample surfaces are used to find the estimate of 0. Good agreement is shown in the figure until about -60~. S2 sample surfaces are 46 A in length with Ax = 0.1A and Oi = 30.0~. Again, good agreement is shown in the figure until about -60~. For surfaces S3 and S4, good agreement is also shown for the bistatic scattering pattern. For surface S3, N = 40, length = 18A, Oi = 30.0~, rms slope = 0.471 and Ax = 0.1A, and the bistatic scattering pattern is shown in Figure 3(a). For surface S4, N = 40,length = 18A = 300, 30.0~, rms slope = 0.354 and Ax = 0.1A, and the bistatic scattering pattern is shown in Figure 3(b). V Conclusion It has been found that using an iterative PO method for characterizing the bistatic scattering pattern for random rough surfaces with low to moderate rms slope (< 0.6), relatively good agreement with the exact solution using both the EFIE and the MFIE is demonstrated. Computation time is significantly decreased as well as memory size necessary to perform the numerical solution for the iterative PO method. Since for 3-D problems the kernel of the integral equation decays faster with distance than that for 2-D problems, it is expected that iterative P.O. to perform even better for 3-D problems. References [1] K. Sarabandi, Y. Oh, F.W. Ulaby "A Numerical Simulation of Scattering from One-Dimensional Inhomogeneous Dielectric Random Surfaces", IEEE Trans. on Geoscience and Remote Sensing, vol. 34, no. 2, pp. 425-432, March 1996. [2] Y. Oh, K. Sarabandi "An improved Numerical Simulation of Electromagnetic Scattering from Perfectly Conducting Random Surfaces," IEE Proceedings on Microwave, Antenna and Propagation, To Be Published. [3] T.B.A. Senior "Scattering By Resistive Strips," Radio Sci., vol. 14, pp. 911-924, 1979. [4] Axline, R.M. and A.K. Fung, "Numerical computation of scattering from a perfectly conducting random surface," IEEE Trans. Antennas Propagat., vol. AP-26, pp. 482-488, May 1978. [5] Oh, Y., K Sarabandi, F.T. Ulaby, "Re-examination of the Kirchoff approximation for scattering from a rough surface", 1993 URSI Radio Science Meeting, Digest p. 406, Ann Arbor, MI, June28-July2, 1993. [6] Poggio, A.J. and E.K. Miller, Computer Techniques for Electromagnetics, Editor R. Mittra, Summa Book, New York 1987. [7] Thorsos, E.I., "The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum", J. Acoust. Soc. Am., vol. 83, pp. 78-92, Jan. 1988. 3

[8] i'sang. L.. I A k'rn 0.9 08 07 306 005 M04 0.3 0.2 0.1 0 ~o!-;.... 4 1 0 ( Z,;, 7';i,, --, If A f.1111 91, - Q, - oo'- - 1-61, 'fe l, t- N.- N,)r,-. 1985. - --- -- I' \ /I J. i I I - 01 3 4 5 6 7 2 3 4 5 6 7 1 Position (in A) 0 Angle (in Degrees) 100 (a) Figure 1: (a) Surface profile with (b) two Gaussian humps. (b) Total bistatic scattering 5 - 101 0 m -5.: -10 -15 8 20 S-25. -30 -35 -40 I -- EFIE - - - MFIE - - Iterative PO method Analytical SPM - 0 m -o 8 --40 |-20 m '-4O - EFIE - -- Iterative PO method - - - Kirchhoff Approximation -Ar I An I -4;, - - -100 -50 0 Angle (in Degrees) 50 100 -do0 -50 0 Ange(in degrees) 50 100 (a) (b) Figure 2: (a) Comparison of bistatic scattering from a random surface with ks = 0.3 ke = 3.0, N = 40, length = 30A, resistive tapered ends length for EFIE = 1A, Ax = 0.2A and Oi = 300. Comparison shows difference between EFIE, MFIE, iterative PO method, and the analytical solution. (b) Comparison of bistatic scattering from a random surface ks = 0.5 ke = 6.13, N = 40, length = 46A, resistive tapered ends length for EFIE = 1A, Ax = 0.1A and Oi = 300. Comparison shows difference between EFIE, iterative PO method, and the Kirchoff Approximation. 5r 10.. 0.v i -5:9 0-10._ m -15 -Y.o -20 m 5 0 c 1o-10 -.S / -- EFIE / — M- FIE -. — Iterative PO method - EFIE -... MFIE - - - Iterative PO method -2Zb &. I I -roo -50 0 Angle (in Degrees) 50 100 -roo -50 0 Angle (in Degrees) 50 100 (a) (b) Figure 3: (a) Comparison of bistatic scattering from a random surface with ks = 1.0 ki = 3.0, N = 40, length = 18A, resistive tapered ends length for EFIE = 1A, Ax = 0.1A and Oi = 30.0~. Comparison shows difference between EFIE, MFIE, and the iterative PO method. (b) Comparison of bistatic scattering from a random surface with ks = 0.5 ke = 2.0, N = 40, length = 18A, resistive tapered ends length for EFIE = 1A, Ax = 0.1A and Oi = 30.00. Comparison shows difference between EFIE, MFIE, and the iterative PO method. 4

Submitted to the IEEE AP-S Microwave Conference 1999. Orlando. FL NUMERICAL SIMULATION OF SCATTERING FROM ROUGH SURFACES USING A FAST FAR-FIELD ITERATIVE PHYSICAL OPTICS APPROACH Daniel Zahn and KIamal Sarabandi Radiation Laboratory, Department of Electrical Engineering and Computer Science University of Michigan, Ann Arbor, MI 48109-2122 Abstract The application of a fast far-field iterative Physical Optics (FIPO) method in conjunction with a Monte Carlo simulation for characterizing the bistatic scattering coefficient of random rough surfaces is examined in this paper. The FIPO method offers decreased memory and computation time restrictions compared to the standard numerical methods such as the Method of Moments (MoM), and decreased computation time compared to an exact iterative PO method. Results from the FIPO method are compared to the standard electric field integral equation (EFIE), the magnetic field integral equation (MFIE), a complete iterative PO (IPO), as well as the existing theoretical solutions for rough surfaces. It is demonstrated that memory requirements and computation time is significantly decreased while providing fairly accurate results for surfaces with moderate to low rms slope. I Introduction Numerically efficient Monte Carlo-based models for simulations of electromagnetic scattering from random surfaces has attained significant prominence over the past decade [4, 2, 1]. A major stumbling block in this endeavor has been the large memory and computation time requirement, but as previously reported, the computation time and memory size has been reduced using a standard iterative PO method which produces fairly accurate results [9]. For a surface consisting of N elements, the iterative Physical Optics method only requires memory size of the order N. Thus, substantial time and memory savings are realized compared to the Method of Moments (MoM). Previously, we investigated the use of the iterative Physical Optics method upon a variety of surfaces in order to find the approximate region of validity for such a method [9]. The difference between a complete iterative PO approach and the new fast far-field IPO (FIPO) approach is that the original IPO was order N2 in computation time, while the new FIPO approach is only order N computation time. Both routines are order N in memory requirements. The results obtained using the FIPO method are also compared to the results found using the electric field integral equation (EFIE) with tapered resistive sheets at the ends of the surface samples, the magnetic field integral equation (MFIE), and the complete IPO method for a horizontally polarized wave. II Formulation For a Monte Carlo simulation, many sample surfaces representing the desired random rough surface characteristics are generated using a standard surface generation routine [4, 2]. As previously reported, the Magnetic field integral equation (MFIE) is used as the basis for the iterative PO method [9] and the currents over a perfect electric conductor is given by [6] I 1 e1k(1r-1) J,(r) = 2 (h x Hi) + J(r')(i(. (r' - r)) (ik - Ir ) ds' (1) where Hi is the incident magnetic field, n is the unit normal at the observation point, and f represents the principle value integral. For the MFIE for two-dimensional scattering problems (one-dimensional roughness), 1

thie transverse (Jt) and longitudinal (J. ) components of thle induced current are decoupled and tlhe integral equations for a TM and TE incidence are respectively given by: ik, p (P - p') d J (p) = 2(n x H') + - f J(n-)()H (k) (op- p'I) d'( 2 p, - d~' (3) -Jt(p) = 2H + Jt,(p)H(1)(kolp - pI)) ' (3) I Jz 2IP - PI where t =, x n is the tangent unit vector. Consider a perfectly conducting rough surface illuminated by an E-polarized (TM) plane wave. Points on the surface are grouped into two categories: 1) lit points and 2) shadowed points. As a first-order approximation the induced current over the lit region is given by J(') = 2(ni x Hi) and over the shadow regions J () = 0. This current produces a scattered magnetic field whose tangent on the surface induces a secondary PO current. The expression for this first order scattered field is given by nx Hl) = fJ ')(p')H()(kolp-') np- p') df' (4) and hence the second order PO current can be obtained from J(2) = 2(n x Hi1)) which is exactly the same as the integral in (2). For the fast far-field iterative PO (FIPO), instead of a complete solution as given in (2) or (3), the surface is divided into a near-field region and far-field regions. Then, (2) and (3) can be written as iko fj(l )(p)H (Pl(kolp-p'l)n p de' p )J(1) ( \ fpl) koIp ) di' J:n)(pn-1) (PI) _ (-' i p- p;jl)) d. - i( 2 Jnearfield P- P I mm 2 Jfarfield,, IP- P I (5) iko (j(nl) (pt)H l)(koIp-pt') (Pl - () M = iko fJ(nEi Jn -). () ( j (-) (kol p)-P )I - p d+ io j("l)(pt))Hl)(ko p-pt[) _ p de' Jnearfieldm IP - PI m=l: mim 2 Jfarfield, 1P -PI (6) Then, the FIPO can be found by summing the contributions from the near-field plus the far-field contributions. The far-field contributions, or last terms, from (5) and (6) may be approximated respectively using the large argument expansion of Hankel functions and noting that in the far field P 'l u where u = x if p' > pc or Ip-p'l wr -i > pa, or = -x if p, < p, iko J (n-1)-ik~(pi OH(1)(kop -p Pi)4 J- l)(p)e-iko(Pd-P)(n'. )d i (7) f arf ield i koH)(kop-p i)(n' u) f (n-1)(p)eiko(P ')ud' (8) Jfarfield where Pmid is the midpoint of the far field section. For the far field, since the fields fall off quickly due to the kernel, large sections of surface are integrated over to give us the far-field contribution to the integral. As shown in (7), for any near-field section, the far-field section only has to be calculated once. Then instead of an N2 calculation the order of calculation is NM, that is, a computational savings factor of (N) is achieved, where M is the total number of far-field sections. As an example, for a 40A surface consisting of 400 elements divided every 5A which leads to m = 8 distinct sections, FIPO is approximately 50 times faster than the traditional IPO. In the actual implementation, the adjacent sections should also be excluded from the far-field 2

approximation. that is. m': {nm - 1n. n +7 1}. Then, for thie example above, for each section m, there are 6 far-field sections if m is equal to I or 8, and 5 far-field sections if 7 is in thle range {2... 7}. For any section. mn, there are 100 or 150 near-field elements which have to be iterated over. and 5 or 6 far-field sections, which are calculated once. The above argument allows one to start with JP) and recursively find a higher order solution. As reported previously there exists a subtlety in the computation of the second order current over the shadow regions [9]. The corrected second order solution is j(2)corrected = J72 over lit region o J ') + 2 (h x Hi) over shadow region From Jz()correc te higher order currents can be obtained. The convergence still is examined by monitoring the normalized difference in the successive solutions ( jE a,) as in [9]. III Results Previously, it was determined that the iterative PO method produced very favorable results for surfaces with small RMS slope. The FIPO method was then compared to the IPO method for a number of surfaces and found to also produce very favorable result. These results are compared to the complete IPO, the Electric Field Integral Equation (EFIE) as well as the Magnetic Field Integral Equation (MFIE). First, the FIPO method was compared to the EFIE, MFIE, and IPO method for a single Gaussian hump, as shown in Figure la. The bistatic scattering pattern for a single Gaussian hump is shown in Figure lb and excellent agreement between the EFIE, MFIE, IPO, and FIPO method are shown. After verifying these results, two surfaces, S1 and S2 were compared for the IPO method as well as the FIPO method. Surface Si is a surface that was measured in [9] and produced excellent agreement with the EFIE and MFIE method. Si has surface characteristics ks = 0.5 and ki = 2.0. In addition, a Monte Carlo simulation of surfaces that have characteristics ks = 0.65, ki = 2.0 and where each surface is 800A in length (or 8000 elements) is compared. For surface S1, the total bistatic scattering pattern (o-o) is shown in Figure 2a. Si sample surfaces are 18A in length, with Ax = 0.1A, the number of surfaces averaged over, N, is 40 and Oi=30.0~. For surface S2, the total bistatic scattering pattern is shown Figure 2b, and the IPO and FIPO methods are compared to result obtained by surfaces with the same characteristics (ks, kr). It is obvious the both the IPO and FIPO methods work very well compared to the EFIE and MFIE methods. For surface 5S2, ks = 0.65, kf = 2.0, N = 30,length = 800A for FIPO and IPO, while length = 20A for EFIE, and 9i = 30.0~. IV Conclusion It has been found that when using a fast far-field iterative PO method for characterizing the bistatic scattering pattern for random rough surfaces with low to moderate rms slope (< 0.6), relatively good agreement with the complete iterative PO method, as well as with the exact solutions using both the EFIE and the MFIE is demonstrated. Computation time is significantly decreased as well as memory size necessary to perform the numerical solution for the iterative PO method. Since for 3-D problems the kernel of the integral equation decays faster with distance than that for 2-D problems, it is expected that iterative PO to perform even better for 3-D problems. References [1] K. Sarabandi, Y. Oh, F.W. Ulaby "A Numerical Simulation of Scattering from One-Dimensional Inhomogeneous Dielectric Random Surfaces", IEEE Trans. on Geoscience and Remote Sensing, vol. 34, no. 2, pp. 425-432, March 1996. [2] Y. Oh, K. Sarabandi "An improved Numerical Simulation of Electromagnetic Scattering from Perfectly Conducting Random Surfaces," IEE Proceedings on Microwave, Antenna and Propagation, To Be Published. 3

[3] T.B.A. Senior "Scattering By Resistive Strips," Radio Sci., vol. 14, pp. 911-924. 1979. [4] Axline, R.M. and A.K. Fung, "Numerical computation of scattering from a perfectly conducting random surface." IEEE Trans. Antennas Propagat., vol. AP-26. pp. 482-488, May 1978. [5] Oh, Y., K Sarabandi, F.T. Ulaby, "Re-examination of the Kirchoff approximation for scattering from a rough surface", 1993 URSI Radio Science Meeting, Digest p. 406, Ann Arbor, MI, June28-July2. 1993. [6] Poggio, A.J. and E.K. Miller, Computer Techniques for Electromagnetics, Editor R. Mittra. Summa Book. New York 1987. [7] Thorsos, E.I., "The validity of the Kirchhoff approximation for rough surface scattering using a Gaussian roughness spectrum", J. Acoust. Soc. Am., vol. 83, pp. 78-92, Jan. 1988. [8] Tsang. L., J.A. Kong, and R. T. Shin, Theory of Microwave Remote Sensing. John Wiley &t Sons, New York, 1985. [9] Sarabandi, K, and D. Zahn "Numerical Simulations of Scattering using an Iterative Physical Optics Approach" Ant. and Prop. Symposium, June 1998 1..... 0.9 0.8 0.7 -0.6 c 0.5 I0.4 0.3 0.2 0.1 J U ) 5 10 15 20 25 Position (in X) 30 35 4 3 (a) (b) Figure 1: (a) Surface profile of a Gaussian humps. (b) Total bistatic scattering Angle (in Degrees) Angle (in Degrees) (a) (b) Figure 2: (a) Comparison of bistatic scattering from a random surface with ks = 1.0 ki = 3.0, N = 40, length = 18A, resistive tapered ends length for EFIE = 1A, Ax/ = 0.1A and Oi = 30.00. Comparison shows difference between IPO, FIPO, and EFIE method. (b) Comparison of bistatic scattering from a random surface using FIPO and IPO ks = 0.65,ke = 2.0, N = 20, length = 800A, Ax = 0.1A and Oi = 30.00 vs. the EFIE for the same physical characteristics, but N = 40, length = 20A with tapered resistive sheets of 1A. Comparison shows difference between IPO, FIPO, and EFIE method. 4