LIBRARY Michigan State University PLACE IN RETURN BOX to remove this checkout from your record. TO AVOID FINES return on or before date due. DATE DUE DATE DUE DATE DUE MSU Is An Affirmative ActiorVEqueI Opportunity Institution CWMI IDENTIFICATION OF MATERIAL CHARACTERISTICS OF LAYERED STRUCTURES BY ULTRASONIC INTERROGATION By Chih- Yeh King A DISSERTATION subnfiued u) hdkflfigan Sane Lhfivenfigy in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Electrical Engineering 1993 ABSTRACT IDENTIFICATION OF MATERIAL CHARAC- TERISTICS OF LAYERED STRUCTURES BY ULTRASONIC INTERROGATION By Chih-Yeh King In this research, the theory and approximations of scattered fields are studied, and mode retraction techniques that can be applied to improve the accuracy of ultrasonic iden- tification with the reflection waveform are developed. Since conventional interrogating algorithms are not valid when dealing with a thin layer, more sophisticated algorithms are needed. First, a wave transmission matrix uses a linearized version of the wave propagation equation and gives an especially simple reconstruction algorithm. This research reviews inverse problem methodology as applied to the backward scattering field, and investigates the quality of the reconstructions when the assumptions behind these approaches are vio- lated. The experimental results show that the evaluation procedures, using a spectral Prony method, are valid only when the wave number is linearly dependent on the fre- quency, and that the specified band-pass spectrum is valid when the incident wave is excited by a Gaussian-pulse shaped acoustic plane wave. ( - mil"... Better identification is based on a higher mode assumption about the reflection waveform. The advantage of this research is that the information obtained is based solely on the echo retum from the successive acoustic interfaces without using the signature of the incident wave which is difficult to capture experimentally. An additional feature of this research is that the results can be used to determine more than one material characteristic from a single trace of reflection waveform. ACKNOWLEDGEMENTS Foremost thanks are due Dr. Bong Ho, my academic adviser, for his guidance throughout this work Gratitude must also be expressed to Dr. H. Roland Zapp, Dr. John R. Deller and Dr. David Yen, for their sincere interest and encouraging advice. Furthermore, I must thank my wife, Wan - Li Ding, for her support and taking care of my parents and our children dur- ing these years of my graduate study. iv TABLE OF CONTENTS page LIST OF TABLES ' Vifi LIST OF FIGURE ix Chapter 1 : INTRODUCTION 1 Chapter 2 : PROPAGATION OF ULTRASONIC FIELDS 5 2.1 Introduction 2.2 Homogeneous Wave Equation 6 2.3 Inhomogeneous Wave Equation 11 2.4 Approximations for Solving the Scattered Field .................. 17 2.4.1 The Born Approximation 18 2.4.2 The Rytov Approximation 19 Chapter 3 : ULTRASONIC IDENTIFICATION OF THE LAYERED STRUCTURE 25 3.1 Introduction 25 3.2 Ultrasonic Identification Schemes 26 3.3 Backward Scattering Measurements ‘ 30 3.4 Angular- Spectrum Techniques 35 3.5 [PM Representation of the Layered Structure -------- 43 3.5.1 Determination of the Wave Transmission Matrix ...... 43 3.5.2 Integral Equation for Surface Reflection Force ~ --------- 52 3.5.3 Spectral Representation of the Reflected Field ......... - 58 Chapter 4 : EXTRACTION OF NATURAL MODE FACTORS FROM A MEASURED RESPONSE 4.1 Introduction 4.2 Least Squares Prony Method 4.3 Maximum Likelihood Method 4.4 MUSIC Method 4.5 Matrix Pencil Method 4.6 SPA Method 4.7 Comparison of Methods Chapter 5 : RECONSTRUCTION PROCEDURES USING THE SPA APPROACH 5.1 Introduction 5.2 SPA Reconstruction 5.2.1 Reflected Signal Perspective 5.2.2 PDF Representation 5.2.3 SPA Approach 5.3 Signal Processing Considerations 5.4 Implementation Considerations 5.4.1 Limitation of the Wave Propagation 5.4.2 Effect of a Finite Receiving Range Chapter 6 : ULTRASOUND FIELD WITHIN AN ARBITRARILY SHAPED AND INHOGENEOUS SCA'I'I‘ERER 6.1 Introduction 6.2 Integration Representation of the Scattered Field --------- 6.3 Evaluation Procedure - 6.4 Determination of Transformation Matrix Elements .......-....- 6.4.1 Off - Diagonal Elements 6.4.2 Diagonal Elements 6.5 Computer Simulations 6.5.1 Dimension of Pixels 6.5.2 Evaluation of Inhomogeneous Field Distribution . -------- vi 60 60 61 65 68 70 73 77 85 86 88 92 94 101 101 102 105 105 107 110 112 112 113 114 114 117 -.I‘-." Chapter 7 : EXPERIMENTS 123 7.1 Introduction 123 7.2 Experiment Facility and Measurement Equipments 124 7.3 Data Acquisition and Processing Procedure ---------------------------- 126 7.4 Experimental Investigation of IPM Validity WWW-"W 130 7.5 Experimental Verification of the SPA Reconstruction ------- 139 Chapter 8 : CONCLUSION 145 8.1 Summary 145 8.2 Suggestion of Future Studies 147 APPENDIX A : ORIGINAL PRONY METHOD 149 APPENDIX B : DERIVATION OF THE DYADIC GREEN’S FUNCTION IN A SOURCE REGION 152 APPENDIX C : EVALUATION OF THE PRINCIPAL INTEGRATION OF THE GREEN’S FUNCTION 155 BIBLIOGRAPHY 158 vii 4.1 : 4.2: 4.3: 6.1 : 6.2 : 7.1 : 7.2: 7.3 : 7.4 : LIST OF TABLES Evaluation of different mode-extraction methods in the presence of 30 & white noise . Evaluation of different mode-extraction methods in the presence of 20 & white noise . Evaluation of different mode-extraction methods in the presence of 10 (fl) white noise . The central displacement field of a (1.5x 1.5x 1.5 m3) homogeneous cube with various n( r) , insonified by a plane ultrasonic wave of difl'er- ent harmonic frequency . Specified values of tissue properties , under a harmonic frequency of 1 M12 used in the computational simulations . The reconstructed natural modes of a 1.475 mm thin plexiglass layer, evaluated via the SPA method . Acoustic properties of plexiglass layers reconstructed by the SPA method . The acoustic properties of different materials published by Kine [110]. Acoustic preperties of different materials evaluated by the SPA method . viii 80 80 81 115 117 134 138 144 144 .ve.’l':‘-.'r -F 2.1 : 2.2: 3.1 : 3.2: 3.3: 3.4 : 3.5 : 3.6 : 3.7 : 3.8 : 3.9 : LIST OF FIGURES A plane wave with direction cosines (‘1 1:02 - kyz , Icy) is propagating along the 7: direction . General scattering configuration . The 2-D backseattering problem to a planar-layered medium . 2-D map of the modulus of the x- and y- components of the scattered field Irr'(x, y)l versus (x, y) within the region of layer 1 . 2-D angular spectrum representation of the Helmholtz eqn. : (a) The medium function M(k), (o) The incident field Ui(k) (c) The spectral convolution mac) n UiUt), (d) The Green’s function G(k) , and (e) The scattered field U5(k) . Estimates of the two-dimensional transforrn of the scattering are avail- able along the solid arc for reflection propagation and the dashed arc for transmission propagation . The problem of the reflection and refraction of a plane wave 011‘ a lay- ered structure . Configuration for a two half-space layers . Configuration for an unbound layer with finite length . Cascaded system of a planar N-layered medium . The ultrasonic schematic diagram for a 2-layered medium in a fluid background . 3.10: The spectral distribution of a pulsed ultrasound wave. 4.1 : 4.2 : The quasi-impulse response of a 4-layered structure , constructed using the first three natural modes with 20 & noise added . Mode factors evaluated from noisy response using LS Prony method with three modes assumed (o) and six modes assumed (+) . (Stars denote the theoretical values) ix 10 14 32 37 38 42 45 48 51 51 83 4.3: 5.1 : 5.2: 5.3 : 5.4 : 5.5: 6.1 : 6.2: 6.3: 6.4 : 7.1 : 7.2: 7.3 : 7.4 : 7.5: Mode factors evaluated from noisy response using SPA method with three modes assumed (o) . (A star denotes the theoretical values) The signal processing steps in a typical SPA reconstruction are shown above . The steps that are needed are shown with a solid box while the optional steps are shown with dashed boxes . The material properties of reconstruction are shown with a 11.675 mm thin Plexiglass layer and the Hamming window added . (All reconstructions shown here are without band-pass filter.) ------- The material preperties of reconstruction are shown with a 11.675 mm thin Plexiglass layer and the Hamming window added . (All reconstructions shown here are with band-pass filter.) The spectrum of the field before (top graph) and after (bottom graph) multiplying by a Hamming window in the frequency domain are shown here . The field scattered by a layered structure is measured along a trans- ducer region with finite dimension . An inhomogeneous scatterer insonified by a plane ultrasonic wave . (a) The three different configuration of a homogeneous bar with n(r) = 10 , f = l Wz and (b) The ultrasonic field rflr) versus the z-axis. (a) Configuration of a fat-muscle plate , (b) The field distribution of the fat plate , and (c) The field distribution of the muscle plate , as a func- tion of spatial coordinates in y-z plane . (a) Configuration of a fat-muscle-fat bar , and (b) The x-component of the ultrasonic field as a function of the axial distance . Schematic representation of the ultrasonic interrogation system . ----- Typical layered structure trace of a reflection waveform insonified by a quasi-Gaussian pulse ultrasonic source . Flowchart for IBM-PC controlled data acquisition and processing . The measured reflection waveform of a 11.675mm thin plexiglass layer placed perpendicular to the transducer axis . One-mode best fit to the spectral ratio of the rear and front echoes , obtained from a 11.675 mm thin plexiglass layer . 84 95 98 100 104 108 116 119 121 125 128 129 131 131 ~4- 7.6 :The measured reflection waveform of a 1.475 mm thin plexiglass layer placed perpendicular to the transducer axis . 7.7 : Phase spectrum of the measured reflection waveform of a 1.475 mm thin plexiglass layer , obtained via the SPA method . 7.8: The measured reflection waveforms and single-mode best fit of plexi- glass layers of different thicknesses: (a) & (c) are 5.365 mm and (b) & (d) are 4.645mm. 7.9: (a) The measured reflection waveform of a plexiglass-aluminum layered model placed perpendicular to the transducer axis . ................... (b) The measured reflection waveform of a plexiglass-copper layered model placed perpendicular to the transducer axis . ................. 7.10: (a) Phase spectrum of the measured reflection waveform obtained from a plexiglass-aluminum layered model. (b) Phase spectrum of the measured reflection waveform obtained from a plexiglass-cepper layered model. 7.11 : (a) The Fourier spectrum of the measured reflection waveform obtained from a plexiglass-aluminum layered model. (b) The Fourier spectrum of the measured reflection waveform obtained from a plexiglass-copper layered model. 3.1: A general configuration of a source region problem . B.2: Evaluating the surface density r1 ( i) by using a differential spherical volume AV . 133 133 136 141 141 142 142 143 143 154 154 .1 s-4_-_'.' 3 1.7 5!: CHAPTER 1 INTRODUCTION The word interrogation generally refers to a type of verbal questioning. In phy- sics, ultrasonic interrogation refers to a procedure used to collect data about the inter- nal structure of an object and to then mathematically generate an image of some other- wise hidden properties of the object. Identification, on the other hand, extracts material properties from the echo retum which consists of incident pulse convolved with the transfer function of the object being investigated. Conventionally, only a single feature such as the interface reflection coefficient, acoustic impedance, attenuation constant, or velocity of prepagation is used for material discrimination. For a more complex structure, however, it is desirable to utilize as many features as possible so that a higher precision of identification can be achieved. Thus we require more sophisticated techniques than the one used for tradi- tional identification. These new techniques for advanced identification are the subject of this work. ' Techniques for the estimation of reflection coefficients at different interfaces have typically involved the manipulation of the echo amplitude in the time domain [1 - 2]. Spectral analyses are commonly used for the estimation of attenuation properties of materials [3 - 6]. However, these techniques suffer from the following shortcomings: (a) Amplitude detection processes have inherent limitations such as the loss of phase information, inability to cape with dispersive media, and low resolution. (b) Spectral resolution is restricted to the reciprocal of the sample spacing in the discrete Fourier transformation [7]. The reconstruction of the geometry and composition of material characteristics based on the measurement of the ultrasound field reflected from a layered structure is known as the inverse reflecting problem [8]. The concept of an inverse problem has also been used in other areas, such as aerodynamics, fluid mechanics, speech and image processing, and electromagnetics [9 - 11]. Since a true inverse problem in acoustics is extremely difficult to solve due to the scattering effect and the mode conversion at the interfaces, one needs to use a simplified model with reasonable assumptions to obtain meaningful results for identifying material characteristics. This work is presented in four parts: the derivation of the wave equation and the first order reconstruction algorithm, the development of the modified Spectral Prony Algorithm (SPA) technique which utilizes unique features in the frequency domain of the reflected response of .a layered structure, the evaluation of the ultrasonic field within an inhomogeneous and arbitrarily shaped object to complete the theoretical background, and finally an experiment verification of the proposed techniques. It is first hypothesized that the frequency domain scattered field response of a acoustic medium can be defined entirely as a function of medium geometry and pro- perties. Then, for a layered structure, the medium function approximately comprises a set of natural modes which uniquely determine the structure. The SPA is a specified phase detection technique constructed in such a way that, upon interaction with a cer- tain structure, results in a reflected waveform contain only a pro-specified component of the natural mode spectrum. The advantages of the SPA technique are twofold. First, it can be used to esti- mate several material characteristics in a single measurement. The information obtained from this technique is solely based on the echo return from the successive interfaces, rather than the signature of the incident pulse which is difficult to capture experimentally. Second, the diffraction effect can be reduced greatly when the reflections from the interfaces are used as references for successive estimation. In Chapter 2, the acoustic wave equation is deve10ped. This vector equation forms the basis of all work to be described. In addition, the Born and Rytov approximations are introduced and a linearized model for the scattered field as a function of the object Chapter 3 presents a brief overview of some of the more relevant identification schemes developed by other researchers, and underscores the need for the SPA tech- nique. A wave transmission matrix is used to find an expression for the reflected field of a layered structure. This forms the basis of the SPA technique which is useful to extract more than one property. The SPA concept is introduced in Chapter 4, and expanded upon in detail. Material in this chapter also includes other extraction schemes and provides numerical results to prove that the proposed method is preferred. Chapter 5 presents a theoretical analysis of SPA reconstruction and includes methods of SPA synthesis and interpreta- tion. Since the mathematical approximations and the experimental limitations contri- bute in different ways to the error in the final outcomes, an overview is also given in this chapter. The mathematical approximations are only valid for some regions of objects, as described in Chapter 2, and the experimental limitations are entirely caused by the limited availability of data. In addition, some signal processing issues will be discussed and experimental results presented. “J -;_. I‘m ,- In certain applications, it is desirable to gain knowledge of the field distribution inside a medium. For example, in hypertherrnia treatment of tumors, the induced tem- perature profile is directly related to the ultrasonic field distribution. Therefore, Chapter 6 follows the wave equations addressed in Chapter 2 and develops a specific algorithm to avoid the singularity of the integration of the Green’s function so that the unique- ness of the ultrasonic field is preserved within the medium. In order to complete this important topic, some inhomogeneous models for solving the given problem are also investigated. Finally, Chapter 7 introduces experimental validation of both the identification of material characteristics of the waveform reflected from thin plexiglass layers, and of the SPA concept itself. Verification of the SPA concept is provided by using the reflection waveform of some multi-layered structures. CHAPTER 2 PROPAGATION OF ULTRASONIC FIELDS 2.1 Introduction Characteristic identification from scattered energy cannot be modeled by a set of equations considering forward propagation wave only. Ultrasonic and electromagnetic waves do not travel along straight rays in an inhomogeneous medium with finite boun- daries. The backward flow of energy can be described with the wave equation under the assumption that the operating wavelengths are small compared with the physical dimensions of the system. It will be shown that the scattering identifications can be approximated by a non-scattering identification. First, let us consider the propagation of waves in homogeneous media. The wave equation is a second order linear difl'erential equation. The field intensity at any given ' location can be obtained by solving the wave equation. The problem is not to identify a homogeneous media but rather to identify one that is inhomogeneous in nature. To solve the inhomogeneous wave equation, one of two approximations, the Born [19] or the Rytov [20], needed to be used. With these two approximations, expressions for the inhomogeneous scattered field can be obtained. nu.‘ The theory to be developed will be applicable to both two- and three-dimensional structures. Even in a three-dimensional case, a two-dimensional model can often be used if the structure varies slowly in the third direction. This assumption, for example, is often made in the evaluation of electromagnetic fields distribution [12]. The theory of ultrasonic identification of materials will be presented almost entirely for the two- dimensional case. The reason is that the ideas behind the theory are often easier to reconstruct the material characteristics in two-dimensional case by using spatial transformation technique. In addition, technology has yet to make it practical to irnple- ment and to evaluate the results of threedimension methods. This limitation will cer- tainly be eliminated in the near future, and where the differences are significant, both the two- and three-dimensional solutions will be expected. 2.2 Homogeneous Wave Equation In a homogeneous medium the propagation of compressional ultrasonic waves can be modeled with the scalar Helmholtz equation [13]. For a temporal frequency of a) radians per second, a scalar potential (N?) satisfies the equation V2 cm + k} em = o . (2.1) For homogeneous media. the wave number k0 is a constant related to the wavelength 1. of the wave by 2 k, .. x . (2.2) The wavelength 3. is related to the temporal frequency of the wave by the propagation speed in the media tic or A = '23-“: . 1 (23) Furthermore, the compressional ultrasonic field 11’0") at position 7‘ is denoted as 4.. -_ H.-- 70") = - VdJC?) . (2.4) The technique developed in this research are based on harmonic ultrasonic fields, the time dependence of the fields will be suppressed in this work. Thus, all fields should be multiplied by e’j‘” to find the measured field as a function of time. For ultrasonic fields. 1??) is the displacement field at position 7’. Since the time depen- dence of the fields has been suppressed, 1??) represents the complex amplitude of the field. As a function of time and space, the field is given by itcr, t) = Re{- Verne-1“} . (2.5) The vector gradient operator, ”V", can be expanded into its two dimensional represen- tation and Eq. (2.1) becomes 2 2 a at? + a £2”) + tie?) = o . (2.6) As a trial solution let em = ei" ' 7’ (2.7) where the vector 2’ = (k, , Icy) is the two-dimensional propagation vector and (DC?) represents a two-dimensional plane wave of propagation factor IPI. This form of (X?) represents the basic function for the two-dimensional Fourier transform [14]. Using this equation, any two-dimensional function can be represented as a weighted sum of plane waves. Calculating the derivatives as indicated in Eq. (2.6), it can be seen that any plane wave that satisfies the condition m=g+g=g am is a valid solution to the wave equation. The homogeneous wave equation is a linear differential equation, so the general solution can be written as a weighted sum of each possible plane wave solution. In two dimensions, at a temporal frequency of to, the potential, M is given by 1 " x 1 " L x m: ElAaykflkx 4*,”de + fiimkflefl-k' +ky)')dky ’ (2.9) and by Eq. (2.8) k, = \lk} - k} . (2.10) The form of this equation might be surprising to the reader for two reasons. (1) The integral has been split into two parts. The coefficients of waves traveling to the right are represented by A (k,) and those traveling to the left by B (ky ). (2) In addition, the limits of the integrals are —oo and on. For It,2 greater than Ira2 the radical in Eq. (2.10) becomes imaginary and the plane wave becomes an evanescent wave. These are valid solutions to the wave equation, but, because It, is imaginary, the exponential has a real or attenuating component. This real component that causes the evanescent waves which decay rapidly within several wavelengths [13], and they can often be ignored. The limited range of valid solutions to the wave equation allows (under certain conditions) an expression to be written for the field in all of 2-D plane given the amplitude of the field along a line. The three-dimensional version of this idea gives the field in 3-D space, if the field is known at all points on a plane. Let us consider a plane acoustic wave travelling in the direction shown in Figure 2.1. By calculating the one-dimensional Fourier transform of the field along the direc- tion of propagation, the field can be decomposed into a number of one-dimensional components. Each of these one-dimensional components can then be attributed to one of the valid wave solutions to the homogeneous wave equation. Referring to Eq. (2.9), there will exist two solutions that satisfy the wave equation with a value of Icy. If the incident field has been defined to prepagate toward the right, then a one-dimensional Fourier component at a value of k, can be attributed to a two-dimensional wave with a propagation vector of ( W , k, ). This can be formulated on a more precise mathematical basis, if the one-dimensional Fourier transform of the field is compared with the general form of the wave equation. When only the forward travelling wave is considered, the general solut ion to the wave equation becomes 1 . cc?) = 31;- ]A(It,)e“"*’ ”malt, . (2.11) Now, if the origin of the wave propagation is chosen to coincide with x = 0, then the expression for the field becomes 1 . om, y) = ?1? [ A(k,)efl‘r’dk, . (2.12) This equation establishes the link between the one-dimensional Fourier transform of the field along the line to the two-dimensional field. The coefficients A (Icy) are given by the one-dimensional Fourier transform of the field A (to) = Fr {M0, n} . (2.13) where FT denotes the Fourier transform. 10 x=x0 x=x1 Figure 2.1: A plane wave with direction cosines (‘1 kg -k§ ,ky) is pr0pagating along the k direction . 11 If the propagation vector is known, then all the sources for the field potential d>(x=tx,,y) = O, e1 “Pu” ) can be decomposed into plane wave components. For a compressional plane wave, the ultrasonic field i? can be expressed as r't’(x=x,,y) = - V(x,, y). ' (2) Obtain the plane wave expression at x = 1:1 by multiplying by the phase factor ekur‘), where as before k, = W. (3) Obtain the field potential d>(x, y) by taking the inverse Fourier transform of the plane wave decomposition. 2.3 Inhomogeneous Wave Equation Considering an inclusion occupying a volume V, in space, embedded in a homo- geneous background medium (see Figure 2.2), the incident wave induces a scattered force [15]. This force, 7; is calculated for an observation point that arises from the spatial variations in density (p) and lame parameters 0. and 11) 12 7:05 = (02890???) + [810”) + 251105] V[v - rm] (2.16) — SuCFDVxVXH’C?) where to is the temporal frequency and 8p('r") , 8M7”) , 511(7) are the position- dependent incremental changes in elastic parameters with respect to the background. They are 59C?) = pt?) - pa . (2.17a) 57m = M?) - x, , (2.17b) and 5M?) = at?) - u. . (me with p, , 1,, and 11, being the density, bulk modulus, and shear modulus of the back- ground respectively. These incremental changes are proportional to the elastic parame- ter contrasts between the background and inhomogeneous media. A more general form of the compressional wave equation can be written as vzocr) + Itfocr) = ‘2'1 v 71(7) (2.18) a) 9. A + 2 where k, = :0)“ and t), = V—L—uf- is the propagation speed in the background. C 0 The term on the right-hand side of Eq. (2.18), is the forcing function of the differential equation [V2 + kg] OO’). Note that Eq. (2.18) is a scalar wave equation. As a result, all depolarization effects can be ignored. It has been shown [16] that the depolarization effects can be ignored when the compressional waves propagate through a viscous compressible fluid. If this condition is not met, the forcing function would be of a much more complicated form [17]. The potential field O0") contains two components -- the incident field O, (‘r") and the scattered field O, (7). The incident field satisfies the homogeneous wave equation 13 [V2 + k,2]O,(T’) = o , (2.19) and the scattered field O, (‘7’) is primarily based on medium inhomogeneities, «h (r) = W) - c. (r) . (2.20) The wave equation then becomes [V2+k,2]O,(i’)= (0;; V-ficr‘) . (2.21) The scalar Helmholtz Eq. (2.21) cannot be used to solve O, (7’) directly. However, a solution can be obtained through the use of the Green’s function technique [18]. The Green’s function, which is a solution of the differential equation [V2 + if] Gm? ') = - ace-7 ') (2.22) is written in free-space as IT” ”M I? ?’| OWL-41:1? ,R- - . (2.23) In two-dimensional problems, the solution of Eq. (2.22) can be written in terms of a zero-order Hankel function of the first kind and can be expressed as 6(7’1? ') = £11}th . (2.24) Notice that the Green’s function GCF’IT’ ’) is a function of the difference between the field point and source point 7’ - 7’ ’, so the function is often represented as G (P-‘P ’). Since the forcing function in Eq. (2.22) represents a point source, the Green’s function can be considered as a single point scatterer. 14 ( Transmitter / Receiver ) homogeneous Background Figure 2.2: General scattering configuration. 15 It is possible to represent the forcing function of the wave equation as an away of impulses suchthat V°fi®=IV’-fl(?’)6(?-?')dv’ . (2.25) In this expression, the forcing function of the inhomogeneous wave equation becomes a summation of impulses weighted by V - 7,0") and shifted by “P ’. Green’s function represents the solution of the wave equation for a single impulse response. The even- tual field solution can be obtained by summing the scattered fields due to each indivi- dual point scatterer. Following the approach, the total field due to the impulse V' - 7, (? var-r ') can be written as the summation of scalar and shifted versions of the impulse response G (P). This is a convolution process and the total field resulting from all sources on the right hand side of Eq. (2.21) is obtained from the following superposition: em = 0,1 [V’-}’,(?')G(?lr')dv’ . (2.25) 90 V Recall that of Eq. (2.4), the compressional ultrasonic scattered field 7'?) at a position 7’ can be expressed as for) = - vo,(r) (2.27) - 1 [VJ’V’TRP‘WCF’IT’W’] 0:29.. I? egir-r'l h c ' = —— w ere (7’ ) 4nl1’-r"l dinate, the unprimed gradient operation can be put inside the integral, . Since the integration is with respect to the source coor- r'cr)=-w21 [IV’-f;(f”)VG(?|?’)dv’]. (2.23) 90 V 16 If the observation point 7 is outside the inhomogeneous volume, V, , then the point 1’ ’ will not pass through ? when carrying out the integration. To carry out the integration, we need to use the dyadic identity [12] V’ - [710”) VG (PI? 3] = [V’ - fit? ')] VG (fl? ') (2.29) +720“)- V’vccrlr') and the symmetry property of the Green’s function V’G (7’1? ’) = — VG (Pl? ') . (2.30) If one substitutes equation Eq. (2.29) and Eq. (2.30) into Eq. (2.28), the scattered ultra- sonic field becomes lrm= .. 021 [j V’- [firivccrmyv' (2.31) a V - {Era-womanly] . l (,2 Defining 8m? ') = We (7|? '). Eq. (2.31) is rewritten as Po em: - [fl -flcr')]vc(rlrw - (2.32) S +£7,(7")- 8(T’l'i”)dv’ . In general. it will be assumed that the volume V, has a finite size; therefore fim is zero outside the volume. Thus, the first surface integration of Eq. (2.32) vanishes at a specified point 7’ outside V,, and Eq. (2.32) reduces to wm=1£m~ tiered)! . (2.33) From the conservation of momentum, 7(7) At = - p(?’) 4%], and when 17 At —>0,thescatteredforce?,(?)is 720') = - 890') [W] (2.34) = 0289557?) = m (mt-rm where m (7) = (0261)?) represents all inhomogeneities of the medium. With this scat- tered force expression, the scattered field becomes an?) = ”(7370”)- 8crlr'w . (2.35) At first glance it appears that this is the solution for the scattered field. However, notice that the scattered field i?’ C?) is expressed in terms of the total field, 11’0") = i“ (7') + i“ (3’). In order to solve the scattered field, some approximations must be made. 2.4 Approximations for Solving the Scattered Field In the last section, an inhomogeneous integral equation was derived for finding the scattered field WC?) as a function of the inhomogeneities of the medium. This equation cannot be solved in a straightforward manner. However, a solution can be obtained by using some approximations. These approximations, the Born [19] and the Rytov [20], are valid under different condition, but the forms of the resulting solutions are quite similar. These approximations are the foundation of solving the ultrasonic scattering problem. Mathematically speaking, Eq. (2.35) is a Fredholrn equation of the second kind [21]. A number of mathematicians have presented work describing the solution of scattering integrals [22, 23]. We will adopt these approximations to the ultrasonic 18 scattering problem with some modifications for reducing computation time. 2.4.1 The Born Approximation The Born approximation is the simpler of the two approaches. Recall that the total field E?) is expressed as the sum of the incident field 7‘ (i’) and a small scattered field ?' (P), rm=r‘(r)+r"(r) . (2.36) The integral of Eq. (2.35) is then fm=£m(?’)7‘(7")-G(?l?’)dv’ . (2.37) +1mcrir'tr3-8crlr'uv' . If the scattered field tr'cr) is small compared to a" (r) , then the effects of the second integral can be ignored and the following approximation can be made: fm=fism=fm(?’)fi"(?’)-6(?|?’)dl/ . (2.38) V As a result, under the Born approximation, the magnitude of the scattered field err) = rm - a“ (r) (2.39) is smaller than the magnitude of the incident field 75(7). If the medium is homogeneous, it is possible to express this condition as a func- tion of the size of the medium and its refractive index. Let the incident wave a" (r) be a plane wave propagating in the direction of the vector, 7:: = f, k, . For a large medium, the field inside the medium will not be adequately approximated by the incident field rm:r‘m=£,rr, e13" . (2.40) 19 In addition, the refractive index n 5 is a variable as well. The field inside the medium has the following variation arr) = IPA, Jam"; ' " . (2.41) The phase difference between the incident field and the field inside the medium can be obtained by integrating the refractive index over the length of the medium. For a homogeneous layer of thickness d , the total phase shift through the layer is (1 Phase Change = 21in 5? (2.42) where it. is the wavelength of the incident wave. For the Born approximation to be valid, a necessary condition is that the change in phase between the incident field and the wave propagating through the layer be less than it [24]. This condition can be expressed mathematically as dn, < 3i“- . (2.43) 2.4.2 The Rytov Approximation The Rytov approximation is another approximation for obtaining the scattered field and is valid under slightly different restrictions. The assumption is that the total field has a complex phase [20] m?) = I? A,e“"’ . (2,44) Recall that the wave equation for the ultrasonic field is Wm?) + k3 (arc?) = - rn (mt-rm . (2.45) In order to evaluate the complex phase, we substitute Eq. (2.44) into Eq. (2.45)and get the wave equation as 20 v - [V¢(F')e m] + k,2(?')em = - m (mm (2.46) or simply [vim]2 + We?) + k} = - m (r) . (2.47) The total phase function 4’ is the sum of the phase of the incident field it, and the scattered phase (1, , 00") = h, to + M?) . (2.48) Substituting Eq. (2.48) into Eq. (2.47), one has [V¢o(?)]2+2V¢o(7’)° Vac?) + [V¢t®]2+ Waco (2.49) +V2¢,(r')+k,2+m(r)=0. The incident field satisfies the homogeneous wave equation, [V , m]2 + v26, (7’) + k} = 0 . (2.50) With this consideration, the over-all wave equation becomes 2mm- vrm +V2¢rf7l = — [veer-mm . (2.51) It can be linearized by considering the following relationship, V2[u‘=‘1,mcrir‘cr')-8crrrw . (2.59) Thus, the complex phase of the scattered field becomes = L a" ' - 6’ Ir . 2.60 Substituting the expression for r” (r) in Eq. (2.38), the Rytov approximation can be written as M?) = i— . (2.61) (5’) 22 In the far-zone region, Eq. (2.58) is valid. Furthermore, "1(7) can be put in terms of the change of the refractive index m(?') = k,2[2n5(?’) + ngm] . (2.62) For small variation of n 5, mm = 2k,2n5(i’) , (2.63) and the Rytov approximation is valid, however, the following condition must be satisfied, 2 [Via (7’)] n 5(7) > ——-2—-— . (2.64) kc This can be considered by observing that to the assumption of Eq. (2.58), the scattered phase ¢,(i’) is linearly dependent on the refractive change; therefore, the second term in Eq. (2.62) above can be safely ignored for a small refractive change. The term V4), (7) is the change of the complex scattered phase per unit distance, while It, is the wave number, the alternate form of the Rytov approximation is W: 0")?» 2 "5(7) > [T] . (2.65) In this result, the change in scattered phase 0, (75 over one wavelength has a dominant effect on the solution of the scattered field. Evaluation of 1? a O") for the Rytov case is not as straightforward as in the case of Born approximation. Recall that the total field is 7(7) = r‘ (r) + W?) = dread-0"“ W’] . (2.66) after rearranging the exponential terms, one has arc?) = I? [A,c’-m[c"m - 1]] = r‘mLW’ - 1] . (2.67) 23 The scattered phase 41,0") now is hm=m Vim +1] . (2.68) 76’) From Eq. (2.61), the field a” (r) for the Rytov approximation becomes tr” =ir" Ingm+l]. 2.69 (7’) (7') [#(7’) ( ) Since the natural logarithm is a multiple-valued function, one must be careful in choosing its value over a given region. For continuous functions, there is no ambiguity about the solution since only one value will satisfy the continuity requirement. On the other hand, for sampled signals, the choice becomes more difficult and one must resort to a phase-wrapping algorithm for choosing the proper phase. Under the Rytov approximation, the total compressional field is expressed as 116’) =A,e’°m”'m . (2.70) Substituting the scattered phase Eq. (2.61), and the incident field Eq. (2.54), into this expression, u (r) = .4, J": "* “'m‘:"“"""3 "7 (2.71) = aim efl'm&"¢xp(-in ° 7) . The first exponential can be expanded into a power series. When u” (7’) is small, the total field can be linearized, u 0’) = u‘Cr')[l + 143(7) - life-fl: ' I] (2.72) = 14‘?) + 113(7) . When the magnitude of the scattered field is small, the Rytov solution approaches the Born solution in Eq. (2.38). 24 The similarity between the expressions for the Born and the Rytov solutions will form the basis of the reconstruction algorithm to be considered in this work. In the Born approximation, the complex amplitude of the scattered field is measured, and this measurement is used as an estimate of the function Pam, while in the Rytov case, i?‘B (i’) is estimated from the complex phase of the scattered field. Since the Rytov approximation is considered more complex [19] than the Born approximation, it might provide a better estimation of Va (7’). CHAPTER 3 ULTRASONIC IDENTIFICATION OF THE LAYERED STRUCTURE 3.1 Introduction The identification of material characteristics in its basic form is an inverse prob- lem - the reconstruction of the geometry and composition of a medium from measure- ments of scattered ultrasonic radiation [26]. The medium can then be identified or discriminated fi'orn other media or classes of media by comparison to known properties and compositions [27]. In a less rigid sense, identification of material characteristics involves the extrac- tion of features from the measured scattered ultrasonic field which correspond uniquely to an individual medium [28]. Identification is then accomplished by comparing these features. Ultrasonic material features for a layered structure medium are closely related to attenuation coefficients, wave velocities in the medium, and reflectivity at each inter- face. These features can take advantage of all the information contained in the scat- tered field such as the amplitude, phase and frequency content. 25 26 Whereas the solution of a true inverse problem requires an infinite amount of informa- tion (re. a measurement of the scattered field at all frequencies) [29], in practice, only a single feature can ordinarily be acquired for the identification of different materials [30]. In Section 3.2, a short overview of some of the more recent methods conceived for ultrasonic material identification is presented. This review is not exhaustive. The goal is to place the identification technique of this thesis in a proper perspective. The two schemes suggested in Sections 3.3 and 3.4 will be compared. It will also be shown that they are evaluations of the ultrasound field at a specified position. Finally, the wave transmission matrix and an ultrasound field integral equation of the ultrasound scattering field, resulting from a surface force induced by the incident field, will be described in Section 3.5. However, the main focus will be on overcoming two drawbacks in ultrasonic identification of materials: (1) The derivation is simply based on the surface reflection force at the planar sur- face of the medium rather than calculated by the complex integration of the 3- D dyadic Green function which is difficult to compute mathematically. (2) The parameter reconstruction can be reduced greatly when the scattered model of the medium is evaluated by using this approach. 3.2 Ultrasonic Identification Schemes The identification schemes reviewed in this section can be generally divided into two fundamental groups. The first involves the use of the amplitude of echo signals, the second utilizes the frequency response of the medium. [A] Reflection Techniques e Correlation Algorithm - In the case of time-harmonic excitation, it is con- venient to consider a layered structure as a "transmission line” [31]. The 27 ultrasonic wave backseattered from the perfectly homogeneous medium can be related to the incident pulse wave through it"(t) = h(t) * u‘(t) (3.1) where " "' " denotes the convolution operation and h (t) is the 1-D impulse response. This equation describes the delays in propagation between the incident and reflected waves, and is a function of the layer thicknesses and wave numbers. The impulse response h (t) of an echo return is in the form of hm=£aw-n) ea n=l where 1:, , properly scaled, corresponds to the locations of the reflecting inter- faces. Thus, the mapping is one-to-one and the reflection coefficients can be exactly evaluated by using the correlation theorem [32] Cueuicc) = $50,, Cuiui(‘tn) (3.3) n=l where " C " denotes the correlation function. e Deconvolution Algorithm -- The discrete component of the reflected signals in the opposite direction of the incident ultrasound wave can be written as [33] u’(n) = £a(k)u’(n - k) + r(n) (3.4) =1 where P is the order of the model and r (n) is an aperiodic pulse train, corresponding to the reflectivity r(n) = 213b(i)5(i - P5) (3.5) i=1 where b (i) is the it” coefficient and P, is related to the 1"“ interface. Then the common autoregressive model used in many digital signal processing [B] 28 applications [34] can be adopted to reconstruct the system r (n ). However, ultrasound identification by reconstructing the reflection coefficients is impracti- cal since measurements could be deteriorated due to the refraction effect. Rather, the reflection property of any interface can be characterized by its dispersion. Frequency Response Techniques Quite a large number of ultrasound identification schemes have been developed based on the frequency response of the medium. Two of these will be con- sidered in this section. e Spectral Difference Algorithm - A very interesting ultrasound identification scheme which simply describes the attenuation factor, a(f) = tr f , has linear frequency dependence [35]. Then the frequency domain response of the biologi— cal tissue, illuminated by an incident wave, can be denoted as [36] |H(f)|2 = e'W)‘ 2” = rm" (3.6) where 20 is the double thickness of the tissue. Using this relationship, the attenuation factor at can be evaluated by a: -10g£|fHD(f)|2] . (3.7) 0 Phase Spectrum Algorithm - An ultrasonic wave prepagating in an unbounded medium in a positive direction can be expressed as [4] u(x,t) =Aej(°"‘“"")e‘°"‘ (3.8) where A is amplitude, to is temporal frequency, 1: denotes the wave number, 9 is the related phase, and (1 denotes the attenuation. It can be shown [5] that any permissible propagating waveform combines with all reverberations can be expressed as 29 +0. '1'“ u(x,t)=% j j Marne-trio ej(°°""‘)e'°“dal . (3.9) For a referred point x = 0, Eq. (3.9) becomes «foo-foo u(0,t)=% j j Mary-tab ejm'do) . (3.10) CI... Applying the inverse Fourier transform technique, one obtains 4. +. U(0,t)= j A((0,t)e"d¢= j u(0,t)e'j°‘dt . (3.11) Similarly, we will get U(x,t) = mom-flue“ . (3.12) With the non-dispersive assumption, the difference in the phase spectrum of these near and far region waves is then alD “c A4) = t), - 9.. = to = (3.13) where D is the spatial distance, 4), and of are the corresponding phase to the near and far regions respectively, and t), is the wave velocity in the medium. As it is not possible to identify all materials with these assumptions, it is neces- sary to approximate the linear frequency dependence within a specified frequency band, and to consider the dispersive properties of the medium with some approxima- tions. These issues will be addressed in the following sections. 3.3 Backward Scattering Measurements In this section, we first review the salient features of the deviation of the 3-D dyadic Green function for a planar layered medium formerly proposed in Chapter 2. 30 When an elastic material is excited by a source, both compressional (C) and shear (S) waves could be launched. Shear wave can be further partitioned into shear vertical (SV) and shear horizontal (SH) polarizations [37]. As the medium is planar layered, C and SV waves are coupled by boundary conditions on horizontal planes (i.e., an incident C wave on the interface between two solids gives rise to both reflected and transmitted of C and SV waves). On the other hand, SH waves do not couple to C and SV waves [38]. Using this decomposition of waves and taking into account the cou- pling at boundaries and the multiple reflections within the layers, it is possible to derive explicit equations for a dyadic Green function for a perfectly layered structure. In particular, it can be shown that when the observation point 7 and the source point r" are in the same layer, 6’ is comprised of a free field ( singular at 7:7“) and a regular part resulting from the contributions of the layered structure. If the source and the observation points are not in the same layer, 6 is comprised of only a regular part. Once 8 is determined, the ultrasound field for an arbitrary source is obtained by applying this operator over the volume of the source weighted by an equivalent force equation (see Eq. (2.35)). For an identification problem, we can denote 17‘ (7) to represent the field at point 7 which was induced by the surface reflection force to compute the ultrasound field. Assume the medium is insonified by a time-harmonic transducer placed in the upper half-space labeled 1, as shown in Figure 3.1, and assume it is filled with fluid. The axis of the transducer is perpendicular to the plane of the figure, and the wave is trav- eling along the if direction as depicted on the figure. As a result, the displacement field anywhere in the space will depend only on the lateral coordinate y and the depth coor- dinate x. The origin is taken at the interface between layer 1 and layer 2. The incident wave is modeled by rim = (4.80! - 1:050” — an (3.14) 31 where 7, = (x, ,y,) locates the incident wave source in the rectangular coordinate sys- tem and u, is the wave magnitude at r,. A scattering force, as defined in Eq. (2.34), represented by a global reflection coefficient of the layered structure can be denoted as its”) = mcr')r(ril,.., (3.15) where mCi’,) is the global reflection coefficient at 7, . With the Born approximation, the scattered field, as defined in Eq. (2.35), becomes r‘crl= [Berri-Kiridv (3.16) V ?’=(x..r..z’) +6. =fo Igcflxo'Yorz’yfdz’ where f, = :71 (3)11, represents the force strength. Due to the spectral decomposition of the dyadic Green’s function as a superposition of plane waves, a 2-D ultrasonic field can be computed by setting the 2 component of the field to zero [39]. Thus, the analysis of the backseattered ultrasonic field can be straightforwardly reduced to a 2-D problem. 32 vi Figure 3.1 : The 2-D backseattering problem to a planar-layered medium. 33 Since layer 1 is a fluid which supports only compressional C waves, the x and y components of 7' become scalar quantities, and they can be written in a more compact forrnas .fO jk) (y-yo) 6.111 (I '4‘.)— u;=(r) —jT (0th In R dky, (3.17) u;m=— «fTo (021pr kye 3310'») e‘jks(x'xe)_ l _R_dky (3.18) where R is the Euclidean distance from (x, y) to (x,, y, ), and the x and y com- ponents of the compressional C wave vector obey k} + It} = k} . (3.19) It is also noted that the scattered field it" is only a function of the longitudinal dis- placement x - x, and the lateral displacement y — y, . A computer simulation gives a 2—D mapping of the modulus of the 3: component of the scattered field versus the coor- dinates x and y, is shown in Figure 3.2. The field mapping was sampled in the x direction at dx = 1m. The sampling in the y direction was dictated by an approxi- mated range of scattering effect, -25° 5 0 S +25°, resulting in dy = .5m; therefore the area shown has the dimensions 8x8 mm2 and was centered at (xT,y1-) = (~64, 0.0). The same computations were performed for the x component of the scattered field. These results provide more detailed information about the distribution of the backseattered filed and the optimal position for locating the transducer to gain a higher reflected energy. Using this approach, the backseattered field can be evaluated with a known global reflection coefficient. It is therefore a useful scheme for ultrasonic identification. 34 0.17 ”3;“ . ’7'}; \\\“‘:. \u is all/7] I,” ;.‘\ “x, 10‘“ ’."I’I I/"II," ’I“ o 0.16 0,“, ll'Il” 'O‘QIII’I7I”% '\\\‘ ‘32] '5 ' IH'I’I'I "I? 'I'. ‘VII’IIII’I I 'E o is .’ 2% "3’27?" 'ititoti’lilqullblf" “ u ”I: 3‘ DD . ’I’ I,” I I." ‘ “\“‘.’[’II III '0", 0‘ ““‘ 31'4”!” 1 v a II,I’ I “0.0““ [ III I], I, be“ ’l/ I u 8 13'3" '4’]? "1”, [III], I.” ‘05:“ ‘u: ’lllllllfll’llll".:“\“\“\:“ ’I [lg/”i ' ‘.;‘ 'Z.‘ 0:“‘0 “$3,: fl", ”0 “‘4‘“ '0 ’ l’ '0‘ o 0.14 \P ’3’"; 2‘ “e:““:::“;:;;;/ 471’”;I'"7;‘.’i‘§i:o:3’303’3'3'3‘gttm '5 I" {I}; ”ll '5' “‘5‘, “at: “ 3:0 ’lll'll’l’lo'g‘t“ :'I"I”'I”Ii'.. ‘: ,‘t‘ . jlll/ 9,”! ”’3'.“ I“ 9’” [I] I] I'.‘ W "I [I III ““\\““\ ' I N 9”], ”III ”I.” \\\ ‘5: III I I, “.33, '4 ’l, ’I,':::‘ u.“ “\r?‘ 3"], I 3 “3 WI! Ill/[I’ll’lll 0’0‘1t‘tt“.“\“::::,’”l,l/Ill ,’I,I,O,u\‘\\?t;I,,"I,f0""“‘\ fizz/(’3' m (,IIIIIIII [I’ll/0'2; “\\‘\“ .“é 0'53, 1,4,,” 10%.“ \, 0,9,3” I Ill/0:0, “\‘\\:‘.? III], I” life} 3“?“ \“‘ 9201 ”I”!!! {fl/'1" [120:0 |.\‘ ““0, ' ’III, ’0 ["34" ‘33:“ ‘3“ "W.” | uy‘ (x, y ) I I4 I], {In}, 711,43; 13:1.“ I H 1' //”/”I I” “‘ § '2 [51:06 'o III] I. s 'Il/I ”Ill/IIII’IIIO" \‘ ~"' '0 ”I [I] ’Iot“t:,;:'I, ’11” go I “WW ‘3‘ ’lllll’l [ill/1W; ‘S‘I‘."I”/lllllli;gill 00““ at u M l’.‘ Ill” 00!”; ”ill/’12“ 0.‘ “in E 6 ail/6’; 0‘ I’ll,” - 23,3 ”/1; ”(39" ’-\“ 9‘? I I .' g - 53:11 %‘,‘\\“ o I" ”II/III ”if I’M; ' Ill/If III/[”0], ’63“ _ 3:40,], ”I if}; '5 t ‘\\I "9’7”,” ”177”,, , ""fl 735,711” Ill/[anal] /' g “1‘ fiqlg.’ I], Ill/”'1'." , ”Wm , 20,049,39/1101/0” l/ 1qu o 2 "‘\\\\\“;o "’;”’///l//’/II ’lt'.‘~"o'l I ”'1’,” Ill/Ill” m ‘t\\“‘:';.'hl"llllllllII‘VHK1’01””0IIIIIII,” 0 " 11““ ”ff,“ "‘t‘t‘. I’I,I, I’ll/r _, 4 \‘l \wle I,’I it’ll; -60 \ V .31: C! -62 -60 0 x - axis _“ 4 y - axis (mm) (mm) Figure 3.2: 2-D map of the modulus of the x and y component of the scattered field lu’(x,y)| versus (x,y) within the region of layer 1. 35 3.4 Angular - Spectrum Techniques If the ultrasound scattered field were in a weakly scattering medium, it could be represented by the following convolution u8®=£m(?')u‘(7”)G(i’-?’)dv’ . (3.20) where Off-T”) is the Green’s function and u,-(?) is the incident field. Thus Eq. (2.35) can be considered entirely in the angular-spectrum [40]. The plots of Fig. 3.3 will be used to illustrate the various transformations that take place. Again, considering the effect of a single plane wave illuminating a perfectly lay- cred medium, the backward Scattered field can be obtained by using Eq. (3.20).The scattered field can simply be considered as a convolution of the Green function g? - 'P ’) with the product of the surface reflection coefficient mm and the incident field u‘ (7’ '). First, define the following spatial Fourier transform pairs mm <-.-> M (I?) . (3.21a) GC?-?')<-—)5(?) . (3.211)) and um (--) (10?) . (3.210) The integral solution to the wave equation can now be written in terms of these spatial Fourier transforms or Uta?) = 6(2‘)[M(I?) * Um] (3.22) where '* " represents convolution, and the propagation vector is I. = k,, Icy). An incident plane wave 14‘ (7‘) can be expressed as u‘ (I?) = 21:50? - 12;) , (3.23) 36 and thus the convolution of Eq. (3.22) becomes a shift in the spatial frequency domain or Ma?) * UR?) = 21: Ma? - 75,) . (3.24) This convolution is illustrated in Figures 3.3(a) - (c) for a plane wave propagating with the propagation vector, I; = (k, ,0). Figure 3.3a shows the spatial Fourier transform of a single cylinder of radius 1 A, and Figure 3.3b is the spatial Fourier transform of the incident field. The resulting multiplication in the space domain or convolution in the spatial frequency domain is shown in Figure 3.3c. To find the spatial Fourier transform of the Green’s function, the Fourier transform of the equation for a point scatterer [V2 + 1:02] G (Pl? ’) = - 8(7’ - 7’ ') (3.25) is taken to find [4:2 + k}]G(E’I? ') = — if? ' 7” . (3.26) Rearranging terms the following expression for the Fourier transform of the Green’s function is found .jf’. r t G PI? _ ‘ 3.27 < ') k2 _ k0, ( ) This expression has a singularity for all 5’ such that I?!2 = k} + k; = k} . (3.28) An approximation of GO?) is shown in Fig. 3.3d. 37 I \ o (a) InCIdent wave 2 0.3 x .‘é ' an 0.6 ~ 3 0.4 q '13 t6 3 0.2 ~ m 0: 20 I ~ (b) Medium function , o (I t) 0.8 I “a ‘ It) '8 Ml) DD u6~ , i ‘ . 2 M‘ .W \l“) ,, '3 IN ‘ \\\"" ~ 0.: M“ l“ g 4 ‘{;;;.;‘ .\\///I ;:.\1 III’ ;;:“ 3‘ I f h‘ “\“h l» (”to ‘3‘?! :02“ “L! ”1'6“: 0 .._ .. E15313”; - .o::.."‘f., If”! W" ’;‘l‘;;\.“ .““‘$’!" 5::“:-;‘ '3: ::v u ‘ ',, ,;';.':i:’3?;._o_3,.: (‘01::- ’. \\ \:- ftft. r" o- ,,-'.' :'.-. no . "15,’.?-?—’.'='§£.-'3 93472;. 15’”!!! 6“? w”! "1%“ “my " ‘ whit?“ .23... "6 3~'-:5:-. ’° '-‘4~:'.'“3-".v'//Il V“ ’1]! ’IO‘“ N I, . “1;",- . s\\\"~"fl ‘6‘3‘555-‘v ’93“ ‘ '.-'-’-'-'-‘«§‘\1*\-'!///,"o'o‘.‘c.I!) 0.4m (adv/o“ \“w \ttu“: - ,o‘ ’ u “w, I 0. .‘\\;._.’..‘.“\'\_;;v.:,‘,” '0. ““0, O .“ l"’1'l':'o::x'f:.‘-‘:".’-" ~. in ‘4’ \\\. 1. o, . I ‘.r. O. ‘,t I 9‘ ‘ .‘H ’I .Q ‘s..l.a '.. “. [II’ . ‘€;f-'.’:‘:":‘33"5:£3"?713‘33’33335124314937" ”In '0 "'- . ..~‘Z.T’-"r‘3”5?°"" ' o -10 . , . - axis .20 .20 kx - axis Figure 3.3 : 2-D angular-spectrum representation of the Helmholtz eqn. (a) The incident field U‘(k), (b) The medium function M(k) (c) The spectral convolution Mar) * U‘(k), (d) The Green function G(k), and (e) The scattered field Us(k) . U 9 P O I -l 1 1 1 p c. 1 Relative magnitude 9 H 3. i (c) Spectral convolution 38 Figure 3.3 (cont’d) '-'vc'l ..... .... ..,o' f“ ’03. 8» -’ .o' ..... .1 o. ".v :: “l 3"“ " p '9 O I I'- I I l p .p 1 Relative magnitude 9 to ‘3: \0 Relative magnitude Icy-axis ((1) Green ’s function I o’ .v‘. ‘4 ”'1' :'¢‘ .:\\q"’..: .: :v'...’ o9 - 30 3 -: 2332‘. 0' ' ‘ J : |\:- Icy-axis ..- a-.-A. ..... (e) Scattered field -"'.'_ ,.- . . r. ...... .“"o ‘.“9‘ . . ‘‘‘‘‘‘ I ' do: . K .. 3‘.- .- 35'. \ I. A; -. -. ,.. n b. 3 T .4-. 4‘: ’1‘...‘ ’0'” . - ’:,:;. 0": 1 (4:94"? 3“\‘:‘: , ,,,,,, 39 The spatial Fourier transfonn representation is somewhat misleading since it represents a point scatterer as both a sink and a source of waves. A single plane wave propagating downward can be considered in two different ways depending on the point of view. From the upper side of the scatterer, the point scatterer represents a sink of the wave while to the bottom of the scatterer the wave is spreading from a source point. Later, when our expression for the scattered field is reflected, it will be neces- sary to choose a solution that leads to backward-propagating waves only. The efl'ect of the convolution in Eq. (3.20) is a multiplication in the spatial fre- quency domain of the shifted medium function, Eq. (3.24), and the Green’s function, Eq. (3.27), evaluated at 7 ’ = (0,0). The scattered field can be written as Mw-a) U'a?) = 2n—— . (3.29 k2 k} ) Figure 3.3e shows the results for a plane wave propagating along the x-axis. Since the largest spatial frequency domain components of the Green’s function satisfies Eq. (3.28), the spatial Fourier transform of the scattered field is dominated by a shifted and sampled version of the medium’s Fourier transform. An expression for the field at the receiver will now be derived. For simplicity, we will continue to assume that the incident field is propagating along the positive 1 axis or. it; = (k, ,0). The scattered field along the receiver line (x = x, .y) is simply the inverse Fourier transform of the field in Eq. (3.29). That is u'(x=x..y)=4—‘, I I maritime, (330) R —.-.o Using Eq. ( 3.27), this can be expressed as 3 _ _ 1 - " MCY-kovky) , "(MW-rail .2..._,,.e e . 40 First we carry out the integral form with respect to 7. For a given Icy, the integral has a singularity at 71.2 = i Vic} - k): . (3.32) Using contour integration, the integral can be evaluated with respect to 7. By adding -1- of the residue at each pole, the scattered field becomes 21: u‘(x,y) = 31;] F1(x,ky)efl‘”dky + 311-:- j r,(x,k,)e""dk, (3.33) where w _ iM water. 1 .k,)- e (3.34a) NEE-T; and -°M «jki- i-Ir ,Ffl r2019): ’ ( ° '9 "Me-1’9””. (3.34h) ZVk} - k} It can be seen that 1‘1 represents the solution in terms of plane waves traveling along the positive 1: axis, while 1‘2 represents plane waves traveling along the negative 1: axis. As discussed earlier, the spatial Fourier transform of the Green’s function Eq. (3.27) represents the field due to both a point source and a point sink, but the two solutions are distinct for receiver lines that are outside the extent of the medium. First consider the case in which the scattered field along the line it = x, is less than that at the x-coordinate of all points in the medium . Since all scattered fields originate in the medium, reflected waves propagating along the negative 1: axis represent a backward propagating waves while waves propagating along the positive x axis represent waves from a point source. Thus for x 0 ( i.e., the receiver is above the medium ) the 41 reflected waves are represented by F2 or 1 u’(x,y)= E! r,(ch,)e”‘"dk, , x <0 . (3.35) Conversely, for a receiver along a line at = x, where x, is larger than the x-coordinate of any point in the object, the scattered field is represented by 1‘1 or 1 . when = 3;] Fr(ch,)e”‘”dk, . x > has... . (3.36) where 1M” represents the physical thickness of the layered structure. In general, the scattered field can be written as , 1 u (x,y) = a! I‘(Jt,lcy)e’7"ydky , (3.37) where it is understood that sign of the square root in the expression for 1‘ should be chosen so that only reflecting waves result. Taking the spatial Fourier transform of both sides of Eq. (3.37), the spatial Fourier transform of the scattered field at the receiver line is Iu‘(x =x,.y)e""‘v’dy = r(x,,k,) . (3.38) Judging from Eqs. (3.34a) and (3.34b), I“(x,,,k,) is a phase shifted version of the medium function. The Fourier transform of the scattered field along the line 3:: = x, is related to the Fourier transform of the medium along a circular arc. The use of the contour integration is further justified by noting that only those waves that satisfy the relationship k} + Ir,2 = k} (3.39) will be propagated, thus it is safe to ignore all waves not on the It), ~circle. 42 Aky reflection transmission; > O 7? Q o o o o o o-.. ..... Figure 3.4: Estimates of the two-dimensional transform of the scatterer are available along the solid are for reflection propagation and the dashed are for transmission propagation. 43 This result is depicted in Figure 3.4. The circular arc represents the locus of all points (kJr .15.) such that k, = :l: k, - k}. The solid line shows the incoming waves for a receiver line at x = are above the medium. This can be considered reflection identification. Conversely, the dashed line indicates the locus of solutions for the transmission case. 3.5 IPM Representation of the Layered Structure The inverse problem method (IPM) is an identification technique first formalized by T. K. Sarkar in 1981 for distinguishing the radar target from the complicated scat- tered field which corresponds uniquely to the individual target [41]. However, its basis is founded in the results of theoretical analysis (see [42 - 44] for examples) which revealed that the scattering field of such an approach is dependent on its geometry and material properties. Later developments in the area of IPM have generally taken one of two direc- tions. Many researchers have investigated the theoretical implications of the method [45 - 47], including the proper expansion of the system responses, while others have employed the method for the analysis of particular ultrasonic problems [48 - 49]. This section will present the theoretical approaches to the IPM technique. These will form the basis for an ultrasonic material identification method which will be investigated extensively in later chapters. 3.5.1 Determination of the Wave Transmission Matrix In general. the ultrasound wave scattered by a planar layered structure results from the surface force induced by the incident wave, is shown in Figure 3.5. Sys- tematic evidence suggests that the surface force can be represented as a wave reflection and transmission problem by a wave transmission matrix (WTM). Let us consider a plane ultrasound wave which is arriving at a planar interface, as shown in Figure 3.5, at some arbitrary angle. We denote the downward wave D‘ as the incident wave from above and the upward wave U ‘ as the back-scattered wave. In order to find the reflection and transmission coefficients at the interface, it suffices to assume the total resultant field is represented by U, =R,D, + T21U2 (3.40a) and D; = T12D1+ R2U2 (3.40b) where the corresponding reflection and transmission coefficients are defined as U D R, = —-1- , 7,, = —2 . (3.4la) D U R2: —3- , 1,, = -—‘ . (3.4lb) U: D.=o U2 D,=o After rearranging, Eq. (3.40b) becomes D e—l-D -£2—U =-l—[D -R U] (342a) l71227122732222. . With this relationship, Eq. (3.14a) can be rewritten as U, = £792,132 - 12,12,112] + 72,112 (3.42b) l = _7' [R102 + [712721- R1R2] U2] - 12 45 layer (1) Figure 3.5: The problem of the reflection and refraction of a plane wave off a layered structure. 46 Resulting fiom the reflection and transmission at the interface, two waves exist inside the upper layer with difl'erent directions of propagation. As a result, the expres- sion for the ultrasound wave, in the layer can be put in the following form Dz = W11 W12 U2 W21 W22 Dz Hz (3.43) [Dr] = __l__ 1 ‘R2 U‘ 712 R1 [leTzl‘RrRz] W W where [W] = [W11 W12] is the wave transmission matrix (WTM). Two cases of 21 22 WTM are of interest: Case A : Two half-space layers Figure 3.6 shows there are only two unbound half-spaces, namely the upper space 1 and the lower space 2, with their acoustic irnpedances Z, and 22 respectively. Assuming the downward wave D, and upward wave U 2 have unity strength, the reflection and transmission coefficients are related by l + R1 = T12 (3°443) 1 + R: = 721 . (3.441» The reflection coefficients are defined as R zz-z, d R 2"22 345 l-Zz-l-Zl a“ 2‘z,+z,' (') u whereZ = pi d istheirnpedance of the 1"“ medium, 0, isthe incident angleatthe ‘ cosfl, i "‘ interface, p, and 1),, are the density and wave velocity in the i "' medium respec- tively. Since R2 = - R1, the elements of the matrix [W] can then be expressed as (3.46a) 47 -R2 Zz-Z, Z2+Z1 Zz-Z, - = .. , .46b W12 T12 [22 + Zl][ 22 22 (3 ) R1 W21 '-'-' “i.— = W12 , (3.46C) 12 and ' = TIZTZI'RIRZ = [1+Rl][l+R2] -R1R2 =.1_ (346d) 22 712 Ta 712 Two half-space layers, the WTM can simply be represented as 1 1 R1 W] = — . .47 [ 1'12 [R1 1] (3 ) Case B : Finite length of an unbound medium For waves propagating through a finite length of medium, as depicted in Figure 3.7, the downward-going wave and upward-going wave are D, = D ,e-jh" (3.48a) y2 = 11,!“ (3.48b) x where lag, = 21:131. 3., is the wavelength of acoustic wave in the i "' medium and C x2, is the physical distance between point x, and x2. Thus one can define WTM, for a finite length of an unbound medium, as [w] = [‘15: 3' [gm] . (3.49) 48 Figure 3.6: Configuration for two half-space layers. Figure 3.7 : Configuration for an unbound layer with finite length. 49 Finally, by using the WTM technique, one can obtain the expression of the reflected wave which is excited by an incident wave as it hits a layered structure of a differing medium. For an N -layered structure with a planar interface as shown in Fig- ure 3.8, the downward wave D,- and upward wave U,- of the i 'h medium are related to the DM and U,“ waves of the i+1"’ layer as Di _1 lRi eflqx‘ 0 U,‘ --]-'i—Ri 1 0 e'jkixi ; 1_<.i .<.N (3.50) ‘ “79‘ 1' eJkixr‘ Rig 1 x _1__ Di+l Ti R- ejkl'xi e-jkixl‘ l Ui-t-l where k,- and x,- are the wave number and thickness of the 1"” medium, respectively. Using algebraic operations, the overall system becomes Dr N 1 ejk‘x" Rte—jk‘x‘ Di+l = — . . 3.51 U1 :1} Ti Rie’k‘x‘ e’mx‘ Ui+1 ( ) = W11 W12 DN+1 W21 W22 UN+1 is the overall wave transmission matrix. W11 W12 where [A] — [W21 W22 For example, for a two-layered medium sandwiched between two half-space back- grounds, as shown in Figure 3.9, the 4‘” medium is infinitely extended. As a result, there is no reflected wave fi'om the 4‘” medium (i.e., U 4 = 0 ). The incident and reflected waves in the upper half-space are in, U1 D4 W11 W12 - W21 W22 04:0 W11 W12 _ 3 1 ejk‘x‘ Rie—jk'x‘ _ r17“:— ROejkl‘xi e—jkixl' t b 50 Therefore, the surface reflection ratio m (7‘) at interface 1 (i.e., x = 0) is then = — (3.53) R,[l + R2R3e'ju’x’] + (”"7“[192 + R3e"2""’] [1 + R2R3e'j2k”’] + R,e""*°’i[R, + R3e'iui‘3] = R, + R2[l - Rae'ju’” + R3[l - R12][1- R22] -12[Iexz + km] e + 51 U1 Dr layer (1) Ui Ri 13' layer (i) f U” . l P” layer (N) Figure 3.8: Cascaded system of a planar N-lavered medium. Dr k1(0)) =k0 l x, W/T Transducer (Transmitter/Rad‘s? W k2(€°) =03: 41"062) 0: 7 U3 92 I A x Figure 3.9: The ultrasonic schematic diagram for a two-layered medium embedded in a fluid background. 52 3.5.2 Integral Equation for Surface Reflection Force The wave number has been assumed to have a linear frequency dependency, the complex wave number of the It“ layer can then be expressed as k, ((0) = [0,, - j or,]m (3.54) where B, and a, are the propagation and attenuation constants of the n" layer respec- tively. Therefore, the first three terms of Eq. (3.53) become m(o))=R,+R2[l-R12]e'jz(fii-ju")m"2 (3.55) + R3[l _ R12][1- R22]e"j2I(Bz ‘jazmz 4' (Br “133M131 = MIC-(01+le)” +M2e’(02+jvz)m + Mse'msi’jvs)” where M, = R,, M, =R,[l -R,2], M3 = R3[l -R,2][1 49,2], 0', +jv, = 0 +10, 02 + iv: = 2am + jzflzxs and 03 + jV3 = 2[Brxr+flszs] +12[asxz+asxs]. This result suggests that the surface reflection coefficient of an N -layered structure can be represented in the frequency domain by the sum of the natural mode functions as mass) = fiM,e'°-°’ (1” (3.56) n=l N = 21‘4““, u=l where N is the number of the interface, Mn is the coupling coefficient of the n“ mode, and 1,, = -o,, - j v,, is the n‘“ complex natural mode factor. Practically, a pulsed ultrasound field can be expressed as r E’Mejm" ; ltl S -A2—t um) =i A, (3.57) 0 ; ltl > 3- 53 where A: is the time duration of the pulse and to, is the harmonic frequency of the ultrasound wave. Figure 3.10 shows the spectral distribution of a pulsed ultrasound wave. The spectrum of this pulsed field can be written as Atl2 17(F‘,co)= j 7(7’,t)e’j°"dt = j i?(7’)e““‘°‘°""dt (3.58) -- -Atf2 . - _ At. - - A: _ [757:2 [6 1(9) «0&2 _eJ(¢o maz] (OJ-(0c 25in [(01) -- mc)-%’—] =31?) (ll—(Dc For a narrow pulse, the following approximation is valid, sin[(al - (nal—3i] 2 __ 21: Ar ~l , Ito-to, SE , (3.59) ((0 - (Dc)— 2 Finally, the surface reflection force becomes 7.0-1. .0» = m(?.m)?(?.w)5(? - e) (3.60) F . 2sin [(0) - 0),)9-5] _ N M ‘11” V? \ 2 - Lug "e 01 0) - me P N a = EMue‘” rice) Lfl=1 J N 1.0) Zn =2Mnfilm)‘ ’ lm-mc S-A—t- ’ 11:1 where 7.0;) = #01, )Ar is the spatial force distribution. 54 9 P b at I I Normalized magnitude 8 a r a P — I o 0.5 l 1.5 2 2.5 3 3.5 4 45 5 Frequency (MHz) Figure 3.10: The spectral ditribution of a pulsed ultrasound wave . 55 Evaluation of the steady-state scattered field excited by the induced surface reflection force can be simplified greatly by the use of a transform domain analysis. The z-transform of a discrete function h (n) is defined as [50] 11(2) = 2 {h (11)} = ill (n )z’" (3.61) n=0 and the corresponding inverse transform is l _ -l = n-l h (n) Z {H (z )} flnj H (z )2 dz (3.62) where the inversion integral is performed over the closed contour in the region of con- vergence of H (z ), and h (n) is taken to be a causal function. If the 1‘“ frequency sam- pleinthespecifiedband,m,=mc-%Stu,Smc+%=(oz,isdenotedas to,=to,+lAco, OSISL-l , (3.63) 02’01 where An) = L 1t . S 't—m—, Tmax = max{|t,|, ltzl, ° - ° , l‘tN'}, and L IS the total number of samples, then the discrete form of the surface reflection force can be denoted as N fire1>=ficit.cor)= zutficme‘”'*"°"- (3.64) n=l N = £Mn(zn)f:.(7:)z,’, . 051 SL -1 1181 Am. where 2,, = e and M', (z,,) = Mnem'“. Transforming these frequency samples into the z-domain, the surface reflection force becomes N fimsl=2m. (3.72) " s z - z, Replacing the normal unit vector i with If", '(F’) ds’ , where f", '(F’) is the conjugate 3 force distribution of the m"' mode, and utilizing the symmetry property of the dyadic Green’s function, then taking the limit 2,, -) z", and l’ Hospital’s rule to Eq. (3.72), the coupling coefficients become 58 M..= gl—jz? we're-Enema (3.73) "'3 where d,” is a normalization coefficient and is given by d... = I dsf‘m’cr’l - 180)?) - no») ds' . (3.74) s s Replacing 2.", by z, the coupling coefficients can be put in terms of the z-domain parameter. M'..lllv+l-'="llv+2='°' =TlL=52- (425) From the expression 11” > 11”,, =n~+2 = - ° - = m = 82 , we may obtain the number of reflections (N). (2) The eigenvectors corresponding to the minimum eigenvalue are orthogonal to the columns of the matrix Z. That is, they are orthogonal to the ”mode vector" of the signals: {VN+1’ ° ' ° , VL}.I. {C(Z 1)’ ° ' ’ , C(ZN)} (4.26) where c(z,,) = [1 2,, z,,2 --- 21'"? is the a“ column of the matrix L We define a VL to be the Lx(L—N) matrix whose columns are the L - N noise eigenvectors. Then, we can estimate the natural mode of each reflection by searching the peak position of the following function. ¢(Z )H C(Z ) (427) P - = . m“ (Z) e(z )” VL Vf’c(z) Properties (1) and (2), hold when the matrix hh" is nonsingular. However, the param- eters (11,, . - - , h”) are coherent in the case of measurements with a layered structure mechanism. This is because the parameters that are identical to the reflection coefficients do not change, but have constant values. Thus, the matrix E (hhH ) is singular, and the MUSIC method does not work properly [72-73]. The construction of E (hhh) has to employ some de-correlation preprocessing to destroy the parameter coherence which makes the problem more difficult. 70 4.5 Matrix Pencil Method Following the idea of the matrix pencil method, we consider the following set of "information” vectors: no, 11,, - - - , 11”, where u,- = [u(i), u(i+l), ~ - - , u(r°+L-M-1)]T . (4.28) The superscript " 7" denotes the transpose of a matrix. Based on these vectors, we define the matrices u, and U, as To look into the underlying structure of the two matrices, one can write where U1 = [“Ov “la ' ° ' . uM-l] and U2=[u1,u2,°-°,uu] . U1 = ZLHZR and U2 = ZLHZO ZR . 1 1 . z1 ZN z, - - _Z{‘-M-1 [Ii—M-l‘ F M-1. 1 Z! M-1 1 7-2 2 ZR - o o ’ i z), . . . “Ali-1 zo =diag[zl. ° ' ’ 92”] r and H=diag[h,, ~ ~ - 1,] . (4.29) (4.30) (4.3 1) (4.32) (4.33) (4.34) (4.35) (4.36) 71 Based on the above decomposition of U, and U2, one can show that if N SM SL -N, then the poles {an n = l,- - - ,N} arethe generalized eigenvalues of the matrix pencil Uz-zU,. Namely, ifN SM SL -N, then z =z,, is a rank- reducing eigenvalue of U, - zU,. To develop and illustrate the use of an algorithm for computing the generalized eigenvalues of the matrix pencil problem we can write U,*U, = zfirrlzzzmz, z, (4.37) =MAM where the superscript * denotes the (Moore-Penrose) pseudo-inverse [74], whereas we use “-1” for the regular inverse. It can be seen from Eq. (4.37) that there exist vectors want = l,- .. ,N suchthat UfU,w,, = w, (4.38) and U,*U2w,, = 2,, w, . (4.39) The w, are the generalized eigenvectors of U2 - zU,. To compute the pseudo-inverse U}, one can use the singular value decomposition (SVD) [75] of U, as follows: N Ul = znntnvrfl (440) n=1 = TDV” U,* = vn-l'r" (4.41) where T=[t,,-- - ,tN], V=[v,,--- ,vN], and D=diag[n,,-- - ,nN]. The superscript H denotes the conjugate transpose of a matrix. T and V are matrices of left and right singular vectors, respectively. Note that for noisy data 11 (1 ) one should 72 choose 11,, - - - , my to be the N largest singular values of U,, and the resulting U? is called the truncated pseudo-inverse of U,. Since U,*U, = W” and v” v = r. substitut- ing Eq. (4.41) into Eq. (4.39) and left multiplying Eq. (4.39) by v” yields [z - 2,, I]z,, = o (4.42) wheren =1, -- ~ ,N, z = 041‘” U,v , (4.43) and 2,, = VHwn . (4-44) Notice that Z is an N xN matrix and that 2,, and 2,, are, respectively, eigenvalues and eigenvectors of Z. It is important to realize that the number of poles, N, can be estimated from the singular values, 1112le2 ° ‘ ’ 2TIN 2 ' ‘ ‘ lenrin(L-M.M) . Since "N44 = ’ ' ' = nm(L_M. M) = O for DOiSClCSS data. IfM =N,theSVD ofU, isnotrequired, andz, ;n =1,--- ,N aretheeigen- -1 values of the NxN matrix [11:70,] U,”U, which is obtained by substituting -l U? = [1150,] U,” into Eq. (4.39). Furthermore, one can verify that with or without noise, -1 [U,"UJ Uf’u, = ° - ° (4.45) which is the companion matrix of the polynomial 73 N 1+ 2a,,z-l =0 . (4.46) n=l Where . 0N aN-l 1 = - [111,110,] Uf’up (4.47) . a.l .1 is the solution of the least square Prony method. So for M = N, the matrix pencil method is equivalent to the LS Prony method. However, in the matrix pencil method, it gives us different values ofM forN S M S L -N. 4.6 SPA Method A very interesting method for natural mode extraction which can be constructed using the frequency-sampling concept introduced in Chapter 3. This section will develop the SPA method such that the frequency response of the backseattered field will pass through a given response. It is a fact that the backseattered field is not a linear phase function, the sampled response must then contain both magnitude and phase. This method uses a criterion based on the SVD algorithm, which is helpful in increasing the accuracy of the mode extraction [76, 77], rather than rely upon the more commonly used LS error approximation between the actual and desired frequency response. Furthermore, it is a useful noncyclic design to avoid the coherence of the measured signals. First we define the z-transform of the ”causal" part of the spectral sequence 11(1) U(z)= iu(1)z" . z e c. (4.48) (=0 74 With Eq. (4.1), this series converges for lzl > 1 to a rational function in the z-plane: -- N U(z) = z[):h,a}.]a-' (4.49) (=0 n=l n=l = fihn [i424] [=0 N = 2])" Z u=l Z - 2,, N EbnzN-n _ 11=0 " N ZanzN-n 11:0 = B z A(z) ’ with a0 = l and b” = 0. Equation (4.49) can be rewritten as B(z)=U(z)A(z) . (4.50) which is the z-transform version of convolution. The convolution can be written as a matrix multiplication. Using the first L terms of the impulse response, we can write q b0 P11(0) 0 0 1 5: um um) . . . . 1 at 0 = 11(N) --°- 11(0) : . (4.51) a. _0 nal-1) - - - u(L-.N-1)_ d To compute the coefficients of a, and b, , we partition the matrices into the form of b U1 1 O “11U2 81 75 where b is the vector of the N +1 numerator coefficients, with b” = 0, of Eq. (4.49), a, is the vector of the N denominator coefficients, u, is the vector of the last L - N terms of the impulse response, U, is the (N+1)x(N +1) partition of Eq. (4.51), and U2 is the (L-N)xN remaining part. The lower L - N equations are written as 0 = u, + U,a, (4.53a) or u, = - U,a, , (4.53b) First, the vector 11,, the denominator coefficient in Eq. (4.49), has to be determined. The upper N + 1 equations of Eq. (4.52), with a = [ 1 a, ], are written as b = U,a . (4.54) Then the numerator coefficients vector b of the transfer function (4.49) can be evaluated. Now the natural mode factors and coupling coefficients are log(z,, ) ,, — A0) (4.55a) and - - _ B (z) . = . . . ha - 211-192(z z") m (Z) 9 n 19 a N 0 (45%) It is a fact that 11 (1) can contain only an unknown but finite number N of exponentials. A finite number L 2 2N of spectral samples 11(1) must be sufficient for the calculation of U (z ), thereby requiring the use of Eq. (4.55) to determine 1,, and 11,. According to Eq. (3.70), 1? ’(P’,o)) is known in a finite interval, the samples 11(1) = i? ’(i’,tu,); 1= 0,1, - - - , L, with LAO) S 0.12 — to, are also available. 76 These properties lead to an effective strategy for determining the natural mode factors of an ultrasound backseattering field out of band-limited measurement data. Starting with L spectral samples 11 (1 ), a rational function (7(2) of order M can be determined using SPA. In particular, the SVD algorithm, which was introduced in both the MUSIC and the matrix pencil methods, can be used for the accuracy estimation of the natural modes in the presence of additive noise. This yields the parameters of 17(2)= "’ ; 210:1 and 5,,=o. (4.56) During the procedure, the order M must be chosen large enough, i.e., M 2 N, to correspond to a first estimate of the expected number of reflections. One must calculate the N largest singular values of U, to evaluate poles 2,, of (7 (z ), and perform a noncy- clic convolution to get the information required for determining the estimates of natural mode factors and coupling coefficients by using Eq. (4.55). Because the preceding methods are an interpolation scheme to design a model that only produces the first L+1 terms of 9 the specified 11 (1 ), these processes say nothing about other points in the passband. To control 1‘1 (2) over the whole range of the passband. we pose an approximation problem where we define an equation-error vector for Eq. (4.52) as follows [i] + [e] = — :11— - -l— . (4.57) If U, has full rank, Eq.(4.57) can be solved for a,. This solution minimizes the lower part of e, and b = U,a gives zero error for the upper part. Then the natural mode factors and coupling coefficients are determined without any bound on 77 resolution. Therefore, it is not worthwhile to attempt to design a maximization scheme as in Sec- tions 4.3 and 4.4, since it is difficult to calculate the peaks of the norm of Eqs. (4.22) and (4.27) with respect to the natural modes. Rather, it is necessary to calculate these using an SVD algorithm, for which many computer routines are available. As with the maximization employed in Sections 4.3 and 4.4, initial guesses are required as a starting point for the iteration process. The main benefit of the SPA method over the ML method and MUSIC method is that no guesses for the complex mode factors are required since they are not used in the construction of the SPA method. The major problem in employing this scheme is, in reality, the quality or the measurement as the spectral samples 1‘1 (1) are perturbed by additive noise which causes variations in the estimation process. The resulting resolution is highly influenced by the choice of the modeling technique for the calculation of the parameters of (7 (z ). by the S IN -ratio in the measured data, and by the number L of spectral samples. For high S IN -ratios, however, some improvements can be achieved through additional echo estimation via parameter modeling. 4.7 Comparison of Methods Five methods for extracting the natural mode factors of the layered structure from a measured response have been presented in this chapter. The sensitivity to random noise and the computational load to numerical implementation have been noted as the motivations behind introducing the SPA method. In this section, a comparison will be made among the various methods, and justification for using the proposed technique will be presented 78 It has been shown that the new extraction routine SPA is successful in the simple case of a single-mode response, while other methods are also known to work well under certain circumstances. Without too much extra effort, we have shown that the spectral Prony method can be extended to cover multi-mode responses. There are con- ditions that the natural-mode extraction schemes must be satisfied. First, it must be able to perform in the presence of random noise, which will undoubtedly be introduced at points in the measuring process. Second, it must be able to discriminate nearly degenerate modes which are quite closely spaced compared with the separation between other modes. Third, it should work well when the number of modes in a measured response is underestimated. The last requirement leads to a particularly useful procedure. Since the number of modes in a measured response will not be known a priori , one would like to start from a small number of modes. For methods which need initial guesses, this small number of initial guesses will be very advantageous. This is only possible, of course, if the routine works well when the number of modes in measured response is underes- timated. To test the various requirements, consider a practical example with a three-mode response given by 11(1) = .6045 + e<--°57-1"-257>f - .6759 e<--‘21-11°-3°4V ; o s f s 5 MHz (4.58) Here, the complex mode factors used to construct 1’1‘(f ) are the first three natural mode factors fi'orn a 4-layer (water-Plexiglas-Aluminum-water) structure. The amplitudes have been chosen to accentuate the reflection coefficients. This response is then sam- pled at 38 equally spaced points between 1.5 MHz and 3.0 MHz. 79 The sensitivity of each of the methods to the presence of random noise can be tested by perturbing each point in the response by a random number, the value of which is defined by L I11 (1)|2 SNR = —— 4.59 and SNR (db) = 10 log,o(SNR) (4.60) where L is the total number of signal samples and 62 is the noise variance. Tables 4.1 - 4.3 show the results of using each method with various levels of noise. In the ML and MUSIC methods, only two phase factors have been evaluated. Using a more com- plex technique would have produced other factors, but would have required more com- putation time. A total of 200 Monte-Carlo runs were conducted in each of the pro- posed methods. For the case of Prony’s method, the DC component in Eq. (4.58) has been removed, since Prony’s method is unable to accommodate a DC component. Tables 4.1 - 4.3 reveal that, with the exception of the LS Prony method, the tech- niques work reasonably well in the presence of random noise, both the matrix pencil method and SPA method are among the best. It is reassuring to see that these two quite different approaches yield similar results. 80 Table 4.1 : Evaluation of different mode extraction methods in the presence of 30 db white noise. M099 01+1Vl 02+1V2 03 ”"3 Method 0 -057 - )7257 -.121 - j10.304 o, —— o, =- .079 i .016 03 =- .271 a .062 PRONY v, _ v, = 7275 a: .015 v3 = 10.441 x .079 a, _ 02 — o3 — MLM v, = -.031 a; .22 v, = 7.226 :t .22 v3 = 10.344 3; .22 o, — 02 — o3 — MUSIC II V, = -.031 1.22 v, = 7.226 a. .22 v3 = 10.344 x .22 o, = 0 02 =- .059 a .005 o, =-.121s.015 PENCIL ll v1 = 0 v2 = 7.258 a .006 v3 = 10.304 1.016 c, = 0 02 =- .057 a .006 03 =- .121 a .014 SPA II V, = 0 v2 = 7.257 :t .005 v3 = 10.305 i .012 20 db white noise . Table 4.2: Evaluation of different mode extraction methods in the presence of ’ M096 01+1Vl 02+1V2 03+1V3 Method 0 -.057 - 17257 -.121 - j10.304 o, -— oz =- .197 i .054 0'3 =- 1.355 1 .435 PRONY v, __ v2 = 7.384 i .053 v3 = 11.249 x .415 a, _ oz — o3 — MLM H v, = -.O31 i .22 v2 = 7.226 1.22 v, = 10.344 :t .22 o, — 02 — o3 _ MUSIC v, = -.031 i .22 v2 = 7.226 :1: .22 v3 = 10.344 1 .22 o, =-.0041.011 02 =- .069t .022 03 =-.188:I:.OS4 FENCE v, = 0 v2 = 7.388 i .053 v3 = 10.340 :1: .060 o, = -.001 3:014 02 =- .056 d: .017 03 =- .133 :t .043 SPA J v, = 0 v2 = 7.260 3: .017 v3 = 10.295 i .045 81 Table 4.3: Evaluation of different mode extraction methods in the presence of 10 db white noise. Mode ll Ol+lVl 02+1V2 03+IV3 Method J 0 -.057 - j7.257 -.121 - j10.304 o, — 62 =- .449 i .182 03 =- 4.560 :1: 1.732 PRONY v, _ v, = 7.549 x .131 v3 = 15.326 a 2.933 01 — 0’2 — 03 — MLM v, = -.053 i .22 v2 = 7.248 :t .22 v3 = 10.322 i .22 01 — 02 — O3 — I . MUS C v, =-.053:22 v2 =7204:.22 v3 =10.366:t.22 o, = -205 :t .511 0; =- .189 a: .156 03 =- .535 i .307 PEN IL C v, = .124: .324 v2 = 7.426: .761 v3 = 10.5462t .799 o, = -.006i .030 02 =- .089 a .106 03 =- .223 3: .151 PA S V1 = 0 V2 = 7.414 i .811 V3 = 10.346 21: .797 Note : 6,, + jv,, Theoretical modes used to construct impulse response PRONY Least square Prony’s Method MLM Maximum Likelihood Method MUSIC MUtiple Slgnal Classification PENCIL PENCIL-matrix Method SPA Spectral Prony’s Algorithm 82 Figure 4.1 shows a simulation of the quasi-impulse response of this 4-layered structure, with 20 1% random noise added. The first three natural modes were used. This represents a typical measured response. We intend to extract the three modes from the noisy response using the SPA method and Prony’s method. Figure 4.2 shows the results of a straightforward Prony method. The simulated results are deviated from the theoretical ones. Improvements can be made by assuming larger number of modes in the response. These results are also shown in Figure 4.2. The accuracy has improved for overestimated cases, but it is very difficult to determine which of the modes are actually present in the response. Without the a priori knowledge of the three mode values it would be difficult to identify the results. On the other hand, Figure 4.3 shows the modes extracted using the SPA method assuming only three modes are present. The results are quite reasonable. The results show that one does not gain much advantage by assuming more than three initial modes. 83 004 r I I I I I I I 0.3 at n o I i I Normalized amplitude b j I :7 1_ 412- 41.5- _o.‘ I J 1 l l I I o 2 4 e s 10 12 14 16 Time axis ( 11sec. ) Figure 4.1 : The quasi-impulse response of a 4-layered structure , constructed using the first three natural modes with 20 db noise added . 84 0- at“. .. +s. \,\ -o.r - - ‘5 -o.2- -\ "at . 3 '\ 4‘3 '\. g: -0.3~ \\ . .g ‘x. ‘3 ‘- 5 -o.4- "o": 3 modes assumed. ‘~ ‘ § "+": 6 modes assumed. \, ., -o.5- "*"z Theoretlcal case. ‘-\ f - 'x. 1' \. -o.a- 4 - l L l l 1 l l o 2 4 e a 10 12 Propagation factor Figure 4.2: Mode factors evaluated from noisy responses using LS Prony’s method with three modes asumed (o) and six modes assumed (+). (A stars denotes the theroretical values.) ‘°-3’ "o": 3 modes assumed. ‘ "*": Theoretical case. Attenuation factor h P l l o 2 4 e a 1 o Propagation factor b p -0.5 Figure 4.3: Mode factors evaluated from noisy responses using the SPA method with three modes (0). (A star denotes the theoretical values.) CHAPTER 5 RECONSTRUCTION PROCEDURES USING THE SPA APPROACH 5.1 Introduction Ultrasonic identification techniques based on natural resonance involve illuminat- ing an unknown layered structure with an arbitrary spectral waveform and analyzing the scattered ultrasonic field. No special restriction is given on the shape of the incident field other than demanding it to carry adequate energy contents to excite natural modes of the target for identification purpose. On the other hand, the angular- spectrum and backward scattering techniques investigate the probability of constructing a specified wavenumber-limited model which would result a wavenumber-limited scat— tered field. This chapter will introduce a method to enhance ultrasonic identification of layered structures by combining the SPA concept with the aspect independence of natural modes. An SPA reconstruction can be described as a phase-detection filter (PDF) which upon excitation of a particular layered structure extinguishes a pre-specified portion of the natural mode phase spectrum. Since the scattered-field response of a layered struc- ture is merely the convolution of the incident field waveform and the structure 85 86 scattered-field impulse response, it is easy to show that SPA reconstruction need not actually be implemented to affect identification. Indeed, two views of SPA reconstruc- tion will be considered. The first describes the effects of directly passing a specified PDF; it uses a z -transform approach. The second describes the more tractable scheme of convolving this PDF with the measured scattered field waveform; this strategy uses a discrete fiequency domain approach. Moreover, the effects of signal processing and implementation considerations will be addressed in Sections 5.3 and 5.4. Preliminary treatment of the SPA reconstruction technique has been carried out in [78]. 5.2 SPA Reconstruction The motivation for adopting a Specified PDF is revealed in the examination of the spectral representation, Eq. (3.79), of the scattered field. Notice that each of the com- plex terms in the natural mode series of the scattered field is multiplied by 11, (I... ). where 11,(z,, ) is the z-transform of the discrete incident field waveform 11, (1 ), and z. is the pole of 11‘“ natural resonance mode of the layered structure. If a specified PDF is defined as A(z,,)=0, forall 1Sn SN , (5.1) it will result in a null scattered-field waveform in the high-frequency samples. Conse- quently, it is easy to reconstruct ultrasonic identification using the SPA technique. Since the PDF is based entirely on the natural modes of a particular object. and since the natural modes of a layered structure are unique, each PDF will then correspond uniquely to one object. Thus, ultrasonic reconstruction using the SPA approach will be easier than analyzing the spectral response of any given scattered field. 87 As in the natural-resonance-based techniques described in Chapter 3, any con- venient incident field with adequate spectral content could be used to illuminate the unknown layered structure. Rather than analyzing the scattered field in terms of its natural mode composition, it would be easier to convolve the scattered field with a PDF defined by Eq. (5.1). Applying the z -transform of the discrete frequency samples of the scattered field to Eq. (3.79) yields fer-0’) N . u’(7’.z)= zvmiL— . (5.2) ”:1 z -z,, Convolution in the discrete frequency-domain is equivalent to multiplication in the z- transform domain, thus we have N 2 {am . u'm1)}=4(z)[z If..§ Padding =——>§ i—> F F T Data 5 Data Window : 1 3 Band 3 pa Interpolation —>§ 13.388 §-——> 5 Filter Reconstruction Figure 5.1 : The signal processing steps in a typical S P A reconstruction are shown above. The steps that are needed are shown with a solid box while the optional steps are shown with dashed boxes. 96 To evaluate the effects of each of these changes, Figures 5.2 and 5.3 Show the estimated natural mode factors using all eight possible combinations of options. The data were measured from a 11.675 mm Plexiglass layer under the ultrasonic wave with an operation frequency of 2.25 MHz . An important part of the reconstruction process is filtering the reflection data. The filter is implemented with an FFI‘ algorithm. These algorithms do not perform an aperiodic convolution like that used in linear system theory, but rather a circular convolution. If the data are padded with zeros so that the new data sequence is twice as long as the original, then the results produced by circu- lar and aperiodic convolution are the same. In addition, zero padding the original reflection increases the resolution of the data in the frequency domain and thus makes an interpolation scheme more accurate. Unfortunately, the extra data require more than double the computational time. For example: an FFI‘ takes N logzN operations so that when N is doubled, the computational expense goes up by a factor of (2N)1082(2N) _ 2 NlogzN -2[1 + logzN] ° (5.28) Based on the reconstruction results shown in Figure 5.2, it can be concluded that dou- bling the size of the reflection data only makes a small improvement in the quality of the reconstruction. Since the extra zeros cause more than doubling the computational expense involved in filtering the data, it is advisable not to zero pad. A second signal processing concern relates to data truncation. In a real-life exper- iment, it is only possible to collect and process a finite amount of data. Generally, this does not present a problem since the data eventually vanishes outside a specific range and the data can be truncated without loss of signal information. The data truncation error can then be mathematically modeled as a multiplication in the time domain by a rectangular window [79]. In the frequency domain this is equivalent to convolving the data with a sine function and thus leveling the frequency domain signal. Other 97 windows like the Hamming window have also been used to reduce the effects of data truncation. Figure 5.2 shows that a Hamming window does not have the same positive effect as the SPA reconstruction. In this case, most of the high frequency information is at either end of the reflected waveform; thus, the window only serves to attenuate the high frequency components. Figure 5.4 shows that the Fourier transform of the reflection is smoothed before (top graph) and after (bottom graph) a Hamming window is applied. The loss of high frequency data caused by the Hamming window leads to lower attenuation factors in reconstruction, as depicted in Figures 5.2 and 5.3. Finally, a very big improvement is observed by adding a BPF before SPA recon- struction of the data. Before the filtering, SPA implementation includes partial terms which also serve to enhance the low and high frequency noises. However, by adding a BPF, this effect is counteracted. Based on the information shown in Figures 5.2 and 5.3, we conclude that the best results are obtained if a BPF is used, but that zero padding the reflection data and applying a Hamming window to the reflection data do not significantly improve the results. 98 No zero padding Mth zero padding Attenuation 0 09 Attenuation .29 . 1 3 (db/ mm-MHz) (db/ mm-MHz) O 3 06 O 3 Acoustic Speed Acoustic Speed 3 2850.03 2848.45 g, (m/sec) (m/sec) E Reflection Coeff. ' 0.3713 Reflection Coeff. 0.3711 2 Acoustic Im . Acoustic Im . ped 3.272 *106 pad 3.270 *106 (Kg lmz—sec) (Kg / mz-sec) Attenuation 0 55 Attenuation 0 1484 .1 5 . (dbl mm-MHz) (db / mm-MHz) 3 .8 Acoustic Speed Acoustic Speed .5 2845.89 2841.15 3 (mlsec) (m/sec). 00 2 Reflection Coeff. 0.3707 Reflection Coeff. 0.3700 3 Acoustic Irn . Acoustic Irn . M 3.267 *106 pod 3.262 *106 (Kg lmz-sec) (Kg / mz-sec) Figure 5.2: The material properties of reconstruction are shown with a 11.675 mm thin Plexiglass layer and the Hamming window added. (All reconstruction shown here are without band-pass filtering) 99 No zero padding VVrth zero padding Attenuation Attenuation 1 .1 2 . 2 3 (db / cm-MHz) O 33 (db / mm-MHz) 0 3 8 O 3 Acoustic Speed Acoustic Speed 3 2754.41 2754.06 on (m/sec) (m/sec) g Reflection Coeff. 0.3565 Reflection Coeff. 0.3565 0 Z Acoustic Irnped. 6 Acoustic Irnped. 6 3.162 *10 3.162 *10 (Kg lmZ-sec) (Kg ”112-sec) Attenuation Attenuation 0.1123 0.1116 (db / rum-MHz) (db / mm-MHz) 3 0 Acoustic Speed Acoustic Speed 3 2756.01 2755.91 3 WM) (m/sec) DD g Reflection Coeff. 0.3567 Reflection Coeff. 0.3567 :1: Acoustic Irn . Acoustic Irn ed. pod 3.164 *106 p 3.164 *105 (Kg lmz-sec) (Kg / mz-sec) Figure 5.3: The material properties of reconstruction are shown with a 11.675 mm thin Plexiglass layer and the Hamming window added. (All reconstruction shown here are with band-pass filtering) 100 2 r r 1 I U E 1 5 _ Reflected Spectrum without Hamming window a | a 1- 3 E 0.6 52 o I l l I 0 1 2 3 4 Frequency axis (MHz) 2 I I I I 15b Reflected spectrum with Hamming window 1 Relative magnitude 0 1 2 8 4 6 Frequency axis (MHz) Figure 5.4: The spectrum of the field before (top graph) and after (bottom graph) multiplying by a Hamming window in the frequency domain are shown here. ' ””1 101 5.4 Implementation Considerations The quality of an SPA reconstruction is limited by both mathematical approxima- tions and experimental limitations. In the derivation of a model for the scattered field either the Born or the Rytov approximation is used to solve UFIE for the scattered field. These approximations are a source of error and limit the applications, as addressed in Chapter 2. The experimental considerations, on the other hand, are caused by a shortage of data. It is only possible to collect a finite amount of data about the backseattering field, and the experimental errors can be attributed to interpolation error and the finite aper- ture. Regardless of the limit in resolution caused by evanescent waves, the limit in quality due to the Born and the Rytov approximations still exist. 5.4.1 Limitation of the Wave Propagation The most fundamental limitation is that evanescent waves are ignored. Since these waves have a complex wave number, they are severely attenuated and phase delayed over a distance as short as a wavelength. This sets the upper limit of the wave number 10 km = 35:3 . ‘ (5.29) This also leads to a fundamental limit of the propagation process, and the range resolu- tion of the SPA reconstruction. These can be improved by conducting the ultrasonic identification at a higher operation frequency (or shorter wavelength). After the wave has been scattered by the target and received by the receiver, the signal must be detected and processed. Unfortunately, it is not possible to sample at every point, so a non-zero interval must be chosen. This introduces an interpolation error into the process. By the Nyquist theorem this can be modeled as a low-pass 102 filtering operation, where the highest evaluated propagation delay, as defined in Eq. (3.63), is given by rm, = 2% (5.30) and A0) is the sampling interval. Of course this analysis has ignored the non-linear pro- perties of wave number and the resulting frequency shifts that occur. 5.4.2 Efl'ect of a Finite Receiving Region Not only are there physical limitations on the finest sampling interval, there is usually a limitation on the region that can be covered. This generally means that sig- nals of the received waveform will be collected at only a portion of signal energy in the receiver region. Consider a single scatterer which is located at the far-zone field region. The wave propagating from this single scatterer is a cylindrical wave in two dimensions or a spherical wave in three dimensions. This efl‘ect is depicted in Figure 5.5. It is easier to analyze the effect by considering the expanding wave to be locally planar at any given point from the scatterer. This can be justified on a more analytical basis by evaluating the phase of the propagating wave. The received wave at a point (x = x, ,y,) due to a scatterer at the origin is given by crew 11 (x 0') = W This instantaneous spatial frequency along the receiver region (5.31) (x,-ASxSx,+A,y,-AySySy,-l-A) can be found by taking the spatial derivative of the phase with respect to y [80] phase = 1:0)le + y2 (5.32) 103 and key k... = _F‘_..2 + ,2 where (cm, is the spatial frequency received at the point (x,y ). From Figure 5.5, it is (5.33) easy to see that tan a = a} , (5.34) and we obtain a... i Jfl (5.35) .\/(x)z+l - ‘11an20+1 I Thus, km is a monotonically increasing function of the angle of view 0. It is easy to see that the maximum received energy can be increased by either moving the receiver closer to the focus line of the object or by increasing the region of the receiver. 104 Figure 5.5: The field scattered by a layered structure is measured along a transducer region with finite dimension. CHAPTER 6 ULTRASOUND FIELD WITHIN AN ARBITRARILY SHAPED AND INHOMOGENEOUS SCATTERER 6.1 Introduction In certain applications, it is desirable to gain the knowledge of the ultrasound field distribution inside a medium. For example, in hypertherrnia treatment of tumors [81 - 83], the induced temperature profile is directly related to the ultrasound field dis- tribution. Therefore, knowledge of the field distribution plays an extremely important role in the success of treatments. Much work has been reported on the evaluation of scattered acoustic fields [84 - 87]. Specifically, the displacement field has been obtained by computing the scattering integral in which the dyadic Green’s function has been used [88 - 90]. However, the field distribution inside the scatterer has not been investigated extensively due to the fact that the boundary conditions of an arbitrarily shaped scatterer are rather compli- cated. Without knowing the intemal field distribution, it is virtually impossible to determine the heat generation and, consequently, the temperature profile inside the medium structure. 105 106 Most of the work published in evaluating the internal field of a medium has been confined to the case of planar layered media with homogeneous and uniform attenua- tion properties [91 - 93], because of their simple in boundary conditions. In other areas, such as the case of underwater acoustics or seismology, one can ignore the boundary conditions because a large-scale problem is being considered [94 - 96]. Recent studies show that the scattering effect of a scatterer, embedded in an inhomo- geneous layered structure, can be solved by using the dyadic Green’s function with the Born approximation [97 - 99]. However, field distribution within the scatterer that remains unsolved could due to the existence of singularity in the integral equation. In this chapter, we seek to relax the restrictive assumptions commonly made on the medium and to design a novel approach to avoid the singularity of the integration. The outline of the chapter is. as follows: (1) In the next section, we develop a general 3-D forward scattering formulation [100] for the purpose of evaluating the scattered field of an arbitrarily shaped scatterer. Then, we decompose the displacement field into an incident field and a scattered field of interest. The incident field results from the ori- ginal excitation in the background medium on the absence of the scatterer whereas the scattered field results from the scatterer. (2) In Section 6.3, we consider an arbitrarily shaped medium embedded in a homogeneous background (such as water), and develop an expression for the 3-D displacement field induced by a plane ultrasonic wave. We then use the proposed technique to preset the coefficients of the transform matrix in Section 6.4. Finally, four simulations are performed. The results show that the fields obtained from the proposed technique are reasonable over a wide range of pixel dimensions. 107 6.2 Integral Representation of the Scattered Field Consider a finite scatterer of arbitrary dimension, with the density p,” (F) and the compressibility c," (7'), illuminated in a background medium (such as water) by a compressional plane ultrasonic wave as shown in Figure 6. l. The induced ultrasonic force in the scatterer gives rise to a scattered ultrasound displacement field 0"(7) , which might be accounted for by replacing an equivalent background force density function 7;, (7) [101] 740,) = - Lmzpmm - m2pa] VG.) (6-1) Cc... (aka?) - c. k3] rm - c'" (711430) —1 [c, 1}] m?) c0 k02 = - nm[c,t,2] at?) = - m (7’) EC?) . The wave number 11,, (7') = 0,, (‘1’) — j a", (7') is a complex quantity governed by the . . . . cm 0" 0) nature of the scatterer. The veloclty 1n the medrum rs u", (7’) = ——,A = —, Pm BM (7") a. (‘1’) is the attenuation coefficient, and p, and c, are the density and compressibility of the backgron medium respectively. Moreover, 11 (fl is the complex deviation ratio of the scatterer and the background, as defined in Eq. (6.1), and 3’0") is the total dis- placement field inside the scatterer. 108 [Vs Inhomogeneous Scatterer 3(a) (pawllcmtin _ , EGfimQFG’X '- _’ Tr’m=u‘(?) +1151?) Figure 6.1 : An inhomogeneous scatterer insonified by a plane ultrasonic wave. 109 In general, the scattered displacement field E" (7) inside the scatterer can be expressed in terms of K9?) by the use of the dyadic Green’s function 80’0”) as follows [102]: 11"(1’) = 174(72‘ 80%? ’)dv’ (6.2) V 1 arm-17c: - ,,_,.,, —3—vv , , to p, 4xl't’-? l where the dyadic Green’s function, GCPI? ’) = 45 represents the transfer function of the background medium, and IZI = m[&-] co denotes the wave number of the background medium. In this case in which the field point is inside the scatterer, fl" (1") should be evaluated with special care because of the existence of singularity and the requirement of uniqueness properties. The scattered field 7" (‘7) at an arbitrary point inside the scatterer can be written as ( see Appendix B) 11"(P')=PV £71,090- c’crlr'w -§LQ (6.3) 30120. where the symbol ”PV" denotes the principal value of the integral. We may now write the total ultrasonic displacement field 3‘?) inside the scatterer as the sum of the incident field 17‘ and the scattered field if as am=rm+rm . (6.4) Substituting Eq. (6.3) into Eq. (6.4) and rearranging terms gives the integral equation for 70’) : [l - m] fiG’H-PV £m(?')i?(?’)- 5(7’17’)dv’ =V‘0") . (6.5) 30129. In Eq. (6.5), if (7') is the incident ultrasonic field and is considered to be a known 110 quantity. 170’) is the unknown total ultrasonic field inside the scatterer, which can be discretized to yield a linear system of equations with complex coefficients. 6.3 Evaluation Procedure Assuming that 7(7) and m (7’) are smoothly-varying functions of space within the volume of the scatterer, we can partition the medium into N pixels, and treat 1??) and m(?) as constants in each pixel. Under these assumptions, the integral Eq. (6.5) can be transformed into a linear system problem. Denoting the m‘“ sub-volume by V,,, and the position of V", by 7:, , we can rewriteEq.(6.5) usingx,=x, x2=y, and x3=z ,as mm.) 1 .- I We. 3 N 114,07.“ 2: 2[m q=1 11:1 - PV 1'prer l?’)dv’:l 11,.(Z,)=u;'(fi,,) V. wherep,q=1,2,3and ,' 82 arm-1‘3; '(7111-?')) 0291: 3x9 81‘11 479$"? " 083.01,]?3: For simplicity, one can define 11:3. =m(7,,)‘ PV JG,',,(2,,I7')dt/ . Then Eq. (6.6) becomes 3 N 11101,) 2 Ag. ”no”, 1-—3m2 u,'(r;)= trim) , q=l 11:1 0 m=1,2,...,N and p,q=l,2,3. (6.6) (6.7) (6.8) (6.9) 111 Denoting [15"]“1‘, as a matrix, its elements are assigned by mm.) 30129. A” -A,',','§' +8Pq8m 1- w. " (6.10) where m, n = 1, 2, ..., N and p, q = 1, 2, 3. [11,5] and [11;] are N-dimensional vec- tors which are given by ”...,m)‘ In}, (7’1) [a5] = I and [14"] = I (6.11) _ux,(?N)j 11;, (TN) where p = l, 2, 3. For all possible values of m and p in Eq. (6.10), the matrix representation of Eq. (6.9) is obtained as . - ~ , . , [in] [4,1 (A...) [11,] 0111 [5,114,115,] 01,] = [11;] (612a) [A3] [A0] [A3] [“2] [”514 [11] = :u‘] (6.12b) where [A] is a 3N>GN matrix, while [11: and [11" have dimension 3N expressing the total displacement field and the incident displacement field at the centers of the N pix- els. Notice that in this way, a 3-D image matrix is represented by a 1-D vector stacked by rows. The inverse problem becomes that of determining the total ultrasound displacement field 170’) from the measurements 1? ‘(7). Therefore, one may evaluate the field distribution inside the medium by solving the Eq. (6.12) with available methods [103 - 105]. 112 6.4 Determination of Transformation Matrix Elements In order to evaluate the coefficients of the matrix [11,, a] , one must be careful in calculating the principal value of Eq. (6.8). Examining the off-diagonal elements of [74%] , we clearly see that r,. does not belong to V, , and prx'fi, I? ’) is continuous throughout V, . Therefore, we may omit the principal value operation from our evalua- tions. However, when we deal with the diagonal elements of [A ‘9‘] , r,,, now belongs to V,ll , and we need to make some approximations. 6.4.1 Off-Diagonal Elements First let us consider the off-diagonal elements of [19%], where p, q = l, 2, 3. From Eq. (6.9), it is clear that 11;"; ”3;. =nr(r:,)- J G,,(2,I7')dv' ; m at» . (6.13) With the same approximation described in Section 6.3, we have .43; as; = 11103,)Gx','(7},, Imav, ; m at n , (6.14) where AV, = I 111’. Using the definition of Eq. (6.7), we can rewrite Eq. (6.14) as v. A'”= 1 en: (-jk,, Rm,” (x"' — x”) (71!" - x1") (6.15) . [3 — (kgkm)2 +j3k0Rm] - (1 +jk,R,,,,,)5pq} ; m at n , whereRM = 117,, —T:,|,while7;, andl’), aredenoted m=rxr+yxg+£xg ,and 7,,=2x','+yx','+£x;. (6.16) 113 When N is chosen properly, reasonable results can be obtained by using Eqs. (6.14) and (6.15). For better accuracy, however, Eq. (6.13) should be integrated numer- ically by any conventional method. 6.4.2 Diagonal Elements For the diagonal elements of [4mm] in Eq. (6.9), we have "103.) 39129. 11;", =m(r;,) PV JG,,,(7;,1?')dv' +8W 1— (6.17) mfi’) = lim 1? dv’ + 8 1- " "1070;“ V. J “.0540: I) P9 360290 For better accuracy, the volume V, has been approximated as a sphere having the same volume as the cubic pixel which is centered at r, . Now, let a, be the radius of the sphere, after steps of manipulation (see Appendix C), Eq. (6.17) becomes mm.) . . mat.) r,., ,, 3m,” [11+1k.a.)exp<-Jk.a.)- 1] + 1- 30), . t J V (6.18) 173 3AV where a, =[ 411“] . If the actual shape of V, differs appreciably from that of a sphere, the approximation may lead to errors. In such a case, the integration throughout a small sphere surrounding r, should be performed by prowdures outlined in Appendix C; the integration throughout the remainder of V, must be carried out numerically. 114 6.5 Computer Simulation Numerical simulations of some simple composite models are performed to test the validity of the theory developed. Other models will also be used to explore the scatter- ing effects within an arbitrarily shaped inhomogeneous medium. 6.5.1 Dimension of Pixels When using a discretized pulse expansion in a linear system of equations, it is crucial to establish a suitable limitation on the dimensions of the sub-volumes. In order to arrive at an optimal choice, we consider a homogeneous bar with a deviation ratio of n (7') = 10, as shown in Figure 6.2a, illuminated by a 1 MHz plane wave with unit strength. Expressed in wavelengths, the dimensions of the this homogeneous bar is 6x%x%P. When the homogeneous medium is partitioned into cubical pixels, the induced field can be evaluated in each pixel. Models having different pixel numbers (N = 12, N = 96, and N = 324), with dimensions of $3., %1., and %h respectively, are used in this study. Computational results confirm with our assumption that the induced ultrasound field is dominantly the 1: component. Figure 6.2b shows the axial ultrasound fields along the f-axis for each model. For the %1 case, the field intensity distribution deviates noticeably from the other two cases. This is because too large a pixel dimension was chosen. For the 71-3. and %71. cases, the fields are symmetrical to the axial axis. A standing wave pattern is esta- blished along the plate due to the finite boundary condition. Since none of the pixels, as shown in the :14. case, lies along the f-axis, we plotted the average field intensity in its surroundings to facilitate a comparison with the cases of £4. and 715-1. The 115 results indicate that pixels with dimension of 7:— }. produce quite an accurate field dis- tribution; it is therefore not necessary to choose yet smaller pixel sizes. We have also investigated the effects of frequency and deviation ratio upon the field distributions. A small cube (1.5 x 1.5 x 1.5 mm 3) with the combinations of different harmonic frequencies and deviation ratios was used. Since the cube presents a symmetrical cross—section to the incident field, we only need to evaluate the induced field in one quadrant The results are summarized in Table 6.1. We observe that the field strength at the center of the cube is nearly independent of the harmonic frequency but is highly dependent on the deviation ratio. Table 6.1 : The central displacement field of a (1.5 x 1.5 x 1.5 m3) homoge- neous cube with various n( r), insonified by a plane ultrasonic wave of different harmonic frequencies. Displacement field a(r) at cubical center 1 1.5 0.5 2.25 11, = 0.1078 - j 0.0708 11y = 0.0179 - j 0.0021 u2 = uy 0.1 15 0.5 22.5 ux = 0.1065 - j 0.0853 uy = 0.0207 - 1 0.0053 ‘12 = uy 0.1 15 5.0 90.0 ux = -0.0468 + j 0.0102 11,. = -0.0021 +1 0.0013 u2 = uy 0.1 15 50.0 765 “x = -0.0036 + j 0.0003 uy = -0.0008 + J 0.0003 “2 = uy Note : Do=(C0/Po)1/2 and limo)=(c,,,(r)/P,,.(r))1/2 . lot. 1 I: I’ll” ¢\ ...... Q .lflll” Ill/If r ..4 O 9 I’ll]. Ill/II . .t. . r. O 9 \\\\\\ill\\\\\\ . .1 \\\\\\ \\\\\\ . \\\\\\ ’ ‘N‘Ni‘N‘N O J. o e . 9 ~ ~ A. \.\.¥ ’ c c 11.. emu 7... a ccc u+ _ 6 2 .t. I”, V\ 0 AAA .1. _ O. — w nun.n. V1.1 _ uni/\t. v in. _~ .11 o .\.\. p .. . r... 1. 9 1.1.4 . turtleseo S rue-n e .1.) \ .r..s+ . o 3.1. p - n p n P h p n 5 6 A 5 3 5 2 6 J 6 0 A o 3 o 2 o 1. 0 D o o 0 o 0 A e a... one 9:885 -‘I z-axisfl.) with ,n( r) = 10, f = 1 MHz, and (b) The ultrasonic field Figure 6.2 : (a) Three different configurations of a homogeneous bar u,( r) versus the z axis. 117 6.5.2 Evaluation of Inhomogeneous Field Distributions Two models composed of biological materials with their properties given in Table 6.2 [106] are used in our simulation. Table 6.2: Specified values of tissue properties, under a harmonic frequency of 1 MHz used in the computational simulations. Biological Velocity Attenuation density material (m/scc) (Np/m) (kg/m3) fat 1445 7 921 (1440 -1490) (5 - 9) muscle 1570 9 1138 (1500 - 1630) (4.5 - 15) (1) Model A : fat-muscle plate. (2) Model B : fat-muscle-fat bar. For the planar model, as shown in Figure 6.3a, we set the incident ultrasonic field to be perpendicular to the surface of the model (i.e., .t-axis). From the results of the previous section, no transverse components need to be considered since we are using a purely longitudinal ultrasonic wave. 118 The fat-muscle square plate has a dimensions 4.5 x 4.5 x .75 m3. The distribu- tions of the ultrasound field inside both layers, versus the f- and 2 - axes, are shown in Figure 6.3b. As expected, the field distributions at the fat-muscle interface are discon- tinuous due to the difference of their acoustic irnpedances. The field intensity is higher in the fat region than in the muscle region. Figure 6.4a shows the 13.5 mm long bar of fat with a 1.125 mm layer of muscle in its center. The axial field distribution is shown in Figure 6.4b. It can easily be seen that the field contains a forward propagation and a reflected component. These waves attenuate exponentially as they travel. As a result, a standing wave pattern is established between the boundaries. For a homogeneous medium, the field intensity is the lowest at the middle of the structure. This is expected, since no boundary exists inside the structure which could cause reflections to trap acoustic energy. On the other hand, structures composed of different materials could produce a local maximum in the field intensity distribution. In addition, if one wants to evaluate the acoustic power variation within a medium, one can relate power to pressure variation. The acoustic pressure can be expressed as a linear function of the phase constant B as pm=pma—“,€3 =jmpmu Ch. 1 b Ch. 2 _: Oscilloscope C3; Tektronix (466) 3% 3 3 Signal interface board 5 3 (Was II) {-.. CO 5 1 On a! water tank Data acquisition, processing, 7// fl and storage :5§:E-’:5§555§55553 Layered ObjCCt M" R\\\ Q (IBM-PC) Figure 7.1 : Schematic representation of the Ultrasonic Interrogation System. 126 7 .3 Data Acquisition and Processing Procedure As mentioned in the previous section, reception of the layered object reflected waveform is provided by means of a single transducer. The reason for using such a transducer is that its simplified configuration can easily be modeled. Moreover, there is no need to be concerned about the timing and amplitude differences between the transmitter and receiver. Generally, the measured signals are unavoidably corrupted by the presence of ran- dom noise. To overcome this problem, the data acquisition process involves making repeated measurement runs and averaging the results. This has been found to increase the signal to noise ratio considerably. However, inconsistencies in the measurements between runs due to triggering jitter, inaccurate reproduction of sampling density, and drift of the preamplifier characteristics and the external power supply voltage necessi- tate the implementation of fairly sophisticated data acquisition and processing software. Figure 7.2 shows a typical trace on the sampling reflected waveform, displaying the reflected waveform from the testing object. Because the AID converter will only accept voltage values between —5 and +5 volts, the first large positive going and nega- tive going peaks of the reflection response must be adjusted so that any portion of the trace to be measured falls within this range. The amplifier gain on the pulser is turned up to near maximum on the signal trace channel, to accentuate the steady-state portion of the reflection waveform. Frgure 7.3 shows a flowchart of the data acquisition and processing routines. The data acquisition phase consists of making multiple (usually 128) separate traces. The sampling density is adjusted such that each data run takes roughly 10 ms to complete. It is important that the sampling density be high enough for the layered structure’s response to be accurately reproduced. Data acquisition software (a C-program) allows the selection of an AID converter sampling rate of up to 40 Mlz . However, the 127 amount of available PC memory puts a limit on the number of points allowed in each data run. The data acquisition phase begins by running the C-program SCANS. This pro- gram controls the AID converter and interfaces with the experimenter via keyboard control. It then 100ps through the entire data acquisition process N titnes. After the data run is finished, the computer gives the operator the option of repeating the run (in case of problems incurred during the sweep), viewing the acquired data, resetting the time-domain window size, and/or writing the data to a floppy disk. It is important to keep the disk writing time to a minimum in order to make the entire data acquisition process as short as possible, because this minimizes the amount of transducer drifting between data runs. Data processing begins by averaging the subsequent data sets and then subtracting the DC offset. If satisfactory processing of all the data sets results in a high signal-to- noise ratio and differentiated response of the layered structure, this response is finally writtentoharddiskastheprocesseddata. ' At any point in the data analysis, the C-program DATAPLOT allows the plotting of a disk file data set on the CRT screen, and also allows a hard copy to be printed. Moreover, it contains a program employing the fast Fourier transform algorithm [108] to obtain the Fourier spectrum of a disk file data set. This allows a preliminary discrimination of the imaginary parts of the object’s natural modes, which can then be used as initial guesses in the more sophisticated natural mode extraction schemes of Chapter 4. One main goal of all the data processing routines discussed is to accomplish as much as possible in the laboratory using the IBM-PC computer. At this time, every- thing but the natural mode extraction schemes can be completed quite rapidly in the laboratory. 128 d d A reflected waveform from a multi-layered medium 0 b. 9 hr 0 L o» o h: Normalized amplitude 6 ho .03 l l l l 0 10 20 80 40 60 Time axis (usec.) Figure 7.2: Typical layered structure trace of a reflection waveform insonified by a quasi-Gaussian pulsed ultrasonic source. 129 Raw data Y Data acquisition — Acquisition of N Graphic plotting waveforms of up Hard disk T—-—-——J to 1024 pts. each _> Rough plots and Sequential storage visual inspection of of each waveform —> waveform sets Storage of the averaged data A f Data processing ‘ i — 1. Removal of DC bias 2. interpolation and averaging of data set on common scale ‘ ' Hard Disk Fully processed data set I Fourier transform — Application of FFI‘ for identification of material properties by SPA method Figure 7.3 : Flowchart for IBM-PC controlled data acquisition and processing. 130 7 .4 Experimental Investigation of IPM Validity The SPA concept discussed in this thesis is based entirely upon the hypothesis that the steady-state scattered field response of a layered structure can be written as in Eq. (3.79) - that is, it can be represented completely by a series of natural-resonance- based damped sinusoid functions. If this expansion is not correct, or if it is not com- plete, then the SPA scheme will fail. It thus becomes crucial to verify experimentally the natural mode behavior of a layered structure. Affirmation of the natural mode description of the frequency-domain response of a layered structure can be accomplished by comparing the extracted natural modes with those obtained from the theoretical IPM analysis in Chapter 3. To the extent that IPM involves implicitly the assumption of a pure natural mode description, close agreement between the two sets of modes would tend to substantiate the IPM recon- struction. This necessitates, however, the use of a layered structure which has been analyzed theoretically. The layered structure should also have clear boundaries and well-known properties, to make the mode extraction easier. Judging from these res- tricted criteria, we decided that the model should be made of thin Plexiglass layers. A visually striking confirmation of the natural mode resonance is shown in Figure 7.4. This is the measured reflection waveform of a thin plexiglass layer of a thickness of 11.675 mm , placed perpendicular to the ultrasonic transducer. Note that the indi- cated early-time and late-time regions are reflections from the front and rear surfaces. It is obvious that this structure response is dominated by a single natural mode which is determined by the product of the double layer thickness (D) and the wave number (I: (to) = fl(co) + j a(to)) of the specified plexiglass layer. Figure 7.5 shows a single mode best fit of the ratio of the rear and front responses from the SPA method. "J.. 131 Ultrasonic reflection from a 11.675 mm thinplexiglass layer. - o lo 0 a d e e N Normalized amplitude o o é 1b 115 2'0 25 Time axis (11860.) Figure 7.4: The measured reflection waveform of a 11.675 mm thin plexiglass layer placed perpendicular to the transducer axis. 2 " -- " :t Reflection spectrum 1.0- "- -": Single-mode curve fit . 1 .e - _ Relative magnitude d l 2 3 Frequency axis (MHz) Figure 7.5: One-mode best fit to the spectral ratio of the rear and the front echoes, obtained from a 11.675 mm thin plexiglass layer. 132 Based on analysis in Chapter 5, the first natural mode of this layer should occur at o, + jo, = -.0376 + j8.522 MHz". The single mode best fit from the SPA method ascribes a mode of 61 + jul = -.0358 + j 8.478 MHz'l. The agreement is extremely close. By thinning the layer it is possible to lower the mode factors and thus to stimu- late the existence of higher order modes caused by reverberation. Figure 7.6 shows the measured surface reflection waveform of a thin plate having a thickness of 1.475 mm placed perpendicular to the ultrasonic transducer. This reflection is apparently not dom- inated by any single mode. The question then becomes whether the reflection can be represented by a pure natural mode series. Figure 7.7 shows evidence for a positive answer. This figure displays the phase spectrum of the reflection waveform obtained via SPA. Three peaks very clearly dominate the spectrum and they can be used for ini- tial guesses in the SPA method. Due to scattering and absorption efl'ects, only the first three modes of this thin layer are absorbable. 133 9 ; Ultrasonic reflection from a 1.475 mm thin plexiglass layer ‘ O in O 0 d N I I —'—= l I .5 Normalized amplitude o b b H .' A 1 J_ l l e a 1B 1 2 1 4 1a 1 a Time axis (psec.) e lo O n. ,1. Figure 7.6 : Measured reflection waveform of a 1.475 mm thin plexiglass layer placed perpendicular to the transducer axis. 2nd echo 0 b 1st echo 0 b V f 1 O O b . 3rd echo ‘ O t- o b Normalized magnitude 0 hi o— l 1 J l l ‘ ‘ -2 -1 o 1 2 a 4 5 Delay time (11sec) , Figure 7.7 : Phase spectrum of the measured reflection waveform of a 1.475 mm thin plexiglas layer, obtained via SPA. 134 The theoretical and evaluated values of these natural modes are summarized in Table 7.1. Table 7.1 : The reconstructed natural modes of a l .475 mm thin plexiglass layer, evaluated via the SPA method. Theoretical Evaluated Evaluated Natural modes value value amplitude 0 .0041 + j0.0092 0.3872 1st mode 2nd mode -0.04 - jl.05 -.0213 - jl.0399 0.3451 3rd mode -0.08 - j2.10 -.0743 - j2.1004 0.0631 The imaginary parts of the natural modes compare well with the theoretical values. However, the damping coefficients are not quite as close. We suspect that the discrepancies are due to the scattering and attenuation effects. Confirmation of the IPM reconstruction can also be provided through reflected waveform measurements. For example, Figures 7.8a and 7.8b show the measured echo reflected by 5.365 mm and 4.645 mm thin plexiglass layers placed perpendicular to the transducer axis. It is obvious that the echo measurements are much noisier than that of 11.675 mm case. In fact, there is a noticeable amount of overlapping between the front and rear echo waveforms. Because of only two successive echoes, a single mode is . ‘ :nf“. our—i." l; l . .... 135 excited. Figures 7.8c and 7.8d show the single-mode best fit to the pass-band portion of the frequency response provided by the SPA method. These results for the thin plexi- glass layers provide a great amount of confidence in the natural mode extraction from the reflection waveform of layered structures. If the reconstruction is indeed complete, then the material properties evaluated from the layer reflection waveforms should be independent of layer thickness. This has been shown to be true in the case of the thin plexiglass layers above. A verification of the thickness independence of the material characteristics of this layered structure can be obtained by making measurements at different thicknesses, and evaluating the material characteristics from each reflection. Table 7.2 shows the experimental results. Material properties obtained from the experiment: such as velo- city, attenuation coefficients, reflection coefficients, and acoustic impedance, are almost the same for four plexiglass layers with different thicknesses. This completes the validation of natural modes reconstruction by using the IPM approach. 136 (a) Ultrasonic reflection from 05'9“” a 5.365 mm thin plexiglass layer ‘ 0.8...F - b .a D I r Normalized amplitude 66 hit as '2 b 11 8 o 5 4 .1 .l 1.0 1‘2 1‘4 116 1. Time axis (”SCC.) (b) 0.001» r . . . , 1 Ultrasonic reflection from - a 4.645 mm thin plexiglass layer WM l - o ‘l 3 U 1 o h 8 o Normalrzed amplitude O 8 O O 43.1... ~05... '- )- '°'7°°°o 2 4 a 3 1o 12 14 1e 10 Time axrs (11sec) Figure 7.8 : The measured reflection waveform and single-mode best fit of plexi- glass layers of different thicknesses: (a) & (c) are 5.365 mm and (b) & (d) are 4.645 mm. (to be cont’d) .‘i J.- Relative magnitude .. O 137 Figure 7.8 (cont’d) (C) " -- ": Reflection spectrum "- -": Single-mode best fit —h-~ 2‘ 3 Frequency axis (MHz) (d) 2.5 1.5 Relative magnitude 05> " -- ": Reflection spectrum "- -": Single-mode best fit 2 3 Frequency axis (MHz) 138 mo ow mo 238an8 a 3 Baotou 03 35.5898 2:. Am 4:2 mm.~ .«o 35:95 .883 a $83: 802.me a 8m eotozaoo mm covenants 05 .eomtafioo me 825:— 05 8m .2 ”802 Seaman 2 Son 2 82 S 83 corneas: e o a. e ... e ... convenes: EEO am. am. am. one .eooo =38qu . . . . 3.55 3 $28 a «3.8 cm $88 R mama £8; mafia. gonna. enoflz. 333. . genesis counseoeafl . ‘ohzaom SE mE : . macs—35. 6052: «mm 05 .3 325288 96.3— mmflwmxo—n mo 8qu89 oumsoo< ”NS 03¢... 139 7 .5 Experimental verification of the SPA Reconstruction With the natural mode representation validated by the results of the previous sec- tion, the next step is to verify the reconstruction of the layered structure experimen- tally. Two different layered su'uctures will be tested. Figure 7.9a shows the reflected waveform of a solid plexiglass-aluminum struc- ture. The geometry is shown in the same figure. Similarly, Figure 7.9b shows the reflected waveform of the plexiglass-copper structure. The Fourier spectrum of the plexiglas-aluminum response is shown in Frgure 7.10a. Three distinct modes are observed in the phase spectrum, as shown in Figure 7.11a. A similar analysis reveals a total of three dominant modes for the plexiglass-copper structure, as shown in Figures 7.10b and 7.11b. The extracted natural modes (U,,-l-j 1),.) corresponding to each struc- ture have been evaluated by the SPA method, and are indicated in Figures 7.9a and 7.9b. The experimental results suggest that as long as the echoes from the front and rear of the layer are not overlapping, the first order model is applicable. The transfer function of a given layer can be put strictly in terms of the properties of that layer, rather than the signature of the incident pulse and the characteristics of the measuring device. The transfer function of the n "' layer is U‘ (0))[A +lejmua+l+jcn+l)] A . Hnfi+l(m) = n _ _fleImW-auflauut) (7.1) vi ((0) [An ‘1 ”(Va +I0u)] - A" where U300) is the incident spectrum and 1),, fi+l+j 6,, n+1 = (1),,+1+j0',,+1) - (u,+jo,.) is the natural mode of that layer. For all cases, the SPA method is employed to evaluate material characteristics. The experimental results are summarized in Table 7.4. For the sake of comparison, published values [109, 110] of the material properties are listed in Table 7.3. It can 140 readily be seen that the velocity, reflection coefficient, and characteristic impedance compare favorable with the published data [109, 110]. However, it is noticed that the attenuation constant deviates from the published value by a significant margin. This is possibly due to the scattering of the incident wave at the interfaces and the absorption of the transmitted wave propagating through the layer. Reviewing the results, we observe that ( 1) The attenuation factors are approximately linear within the frequency band, 1.5 m1: (7’) as i’ " (7’) = —Vd>(?). If AV is small enough, the quasi-static approxima- tion can be used to evaluated (1) as _ l ,. , ¢XP('I7‘;'A7') M7)— (ozp, 5.1/v £40,) 41tlA?’| W (BA) whiz ° A7”) if, 410A? ’l +IM® AS where A?=T’-‘P’ is the distance between the source point and the field point, and AS is the surface of AV. It is evident that the first term of Eq. (B.4) goes to zero as AV approaches zero. However, th can be shown to converge to a finite value as AV approaches zero. This implies that V " C?) remains a finite value at the center of AV. For simplicity, we will evaluate the term -Vd> for a spherical volume AV as follows ?‘®-_£lim 91 I” 'k flds’ s - nMInU 22 exP(-Jo€)COS . (3.) -un. .. »_H‘-- ...u-——-————.-_—_ n a I 154 henna Region ‘ * ~ ’fj...(,p,,,~(r>), elm) .1 1 _ Figure B.l : A general configuration of a source region problem. AV = 4?1l5 a: an :Radius of sphere AV Figure B.2: Evaluating the displacement density 110) by using a differential spherical volume AV. APPENDIX C EVALUATION OF THE PRINCIPAL INTEGRATION OF THE GREEN ’8 FUNCTION Using a strategy similar to the one described in Appendix B, we modify V, as a sphere with the same volume centered at 7,, shown in Figure B.2. The radius of the sphere is denoted by 3AV, "3 a, = 4n . ((2.1) Considering the property of equality : V’OCPI‘P') = - VGW’) , Eq. (6.7) can be rewritten as (C.2) 2 -°? - —r G,” (3'7”): 21 ?'—, [exp( 1 a C? 3)] r ' (1) p, 3x, 3X, 41tIP-T’ ’l p, q =1, 2, 3 . Therefore, we are allowed to treat the point?=T, as the new origin with the transfor- mation of coordinates. Eq. (C.2) then becomes l a2 [apt-1K?) —— , C.3 (02% qu’ax,’ 47"? " ] ( ) G‘P‘Cm ’r')=erxt(r')= p,q=l,2,3 . 155 156 In the spherical coordinates, ———.::,l-~lAtari—tim]livelier-J. p,q=l,2,3 . The principal value of the integration becomes a. 2x x{xx PV G , I?’ = d d J’J- as it”??? W (Dz-”Poe “3,11% r2 () arm-1k: 7) l x lfwi ]*7[‘~’if] -73. . £_[exp(4)n: fl]}r2sin ode . Carrying out the integration by parts, it becomes O4 a d exp (’Jka ' 7") ' PV G ’)dv‘ = lim 2— . J. z'x'mlr (0po e-io {[r dr [ 41w H; (C 6) 21: xxx -jdpj-E-2-"—sinedo 0 0 1' 2: x xx +£d¢£ fipq-3-Ef-sin6d9 ‘- F’ d exp(-Jko°fl '£E[ 4.. ]”"} 157 arm-1'13 ° 7’) 41tr 0. O 0 d Althoughtheequatronlrgb i rdr[ 21: x xx jd¢j[5m-3J—2¢-]sinede=o ; p.q=l,2,3, o o r and 2: x xpxq - _4_fi . .. dd)! r2 srnOdO- 35,, , p,q-l,2,3 . Therefore, the Eq. (C.6) becomes -6 4 d apt-j? WH} = _H. - _E 2... 0 PV £0,541! (0po .131-r3{ 3 [r dr[ 4” a -5pq . . S 3029., [(1 + Jka ankxp we a..) - 1] _y where k, '7’ ,ggSkoan . Jdr —-) no, it is readily verified that (C.7) (C.8) (C.9) - C—AJ“- “1-9-7- , ‘. BIBLIOGRAPHY “‘1. mum-7 *lt'fi [l] [2] [3] [4] [5] [6] [7] [8] BIBLIOGRAPHY C. Gazanhes, J. P. Herault, and K. Stephankis, ”Reflection coefficient identification by means of correlation: Applications to layered medium," J. Acoust. Soc. Amer., vol. 69, No. 3, pp. 720-727, Mar. 1981. M. Cookey, H. J. Trussel, and T. J. Won, "Seismic deconvolution by multi-pulse methods," IEEE Trans. AS.S.P., vol. 38, No. 1, pp. 156-160, Jan. 1990. M. Fink, F. Hottier, and J. F. Cardoso, " Ultrasonic signal processing for in vivo attenuation measurement : short-time Fourier analysis,” Ultrasonic Imaging, vol. 5, pp. 117-135, Apr. 1983. R. Kuc, "Estimating acoustic attenuation from reflected ultrasound signals : Com- parison of spectral-shift and spectral-difference approaches," IEEE Trans. AS.S.P., vol. ASSP-32, No. 1, pp. 1-6, Feb. 1984. W. Sachse and Y. H. Pao, "On determination of phase and group velocity of dispersive wave in solids,” J. Appl. Phys, vol. 49, pp. 4320-4327, Dec. 1978. T. Pialucha, C. C. H. Guyott. and P. Cawley, ”Amplitude spectrum method for the measurement of phase velocity,“ Ultrasonic, vol. 27, pp. 270-279, Sep. 1989. R. Kumaresan and D. W. Tufts, ”Estimation of frequencies of multiple sinusoids making linear predication perform like maximum likelihood,“ Proc. IEEE, vol. 70, No. 9, pp. 975-989, Sep. 1982. S. W. Ross, "Review and examination of results on uniqueness in inverse prob- lems,” Int. Symp. on Electromagnetic Theory, , pp. 1026-1030, Aug. 1986. 158 159 [9] W. R. Smith, "Mathematical analysis, an inverse problem arising in convective- diffusive flow,” J. Appl. Math, vol. 45, No. 3, pp. 225-231, Mar. 1990. [10] Y. M. Ram and S. G. Braun, "Inverse problem associated with modification of incomplete dynamic systems," Trans. on ASME, vol. 58, No. 1, pp. 233-237, Mar. 1991. [l 1] J. A. Fawcett, ”Inversion of acoustic plane wave reflection coefficients for bottom elastic parameters,” Wave Motion, vol. 12, No. 5, pp. 451-460, Sep. 1990. [12] C. A. Balanis, Advanced engineering electromagnetics, New York: John Wiley & Sons, 1989. [13] R. A. lemons and C. F. Quate, ”Acoustic microscopy,” chap. 1 in Physical acoustics: principles and methods, vol. XIV, pp. 1-92, W. P. Mason and R. N. Thurston, eds., New York: Academic Press, 1979. [14] R. N. Bracewell, The Fourier transform and its applications, 2nd ed., New York: McGraw-Hill. 197 8. [15] w. M. Ewing, w. s. Jardetsky, and F. Press, Elastic waves in layered media, New York: McGraw-Hill, 1957. [16] A. Ishimaru, Wave propagation and scattering in random media, New York: Academic Press. 1978. [17] A. C. Kak, 'Tomographic imaging with diffracting and non-diffracting sources," in Array signal processing, 8. Haykin, ed., Englewood Cliffs, NJ: Prentice Hall, 1984. [18] P. M. Morse and H. Fesfbach, Methods of theoretical physics, New York: McGraw-Hill, 1966. [19] K. Iwata and R. Nagata, ”Calculation of refractive index distribution from inter- ferograms using the Born and Rytov’s approximations,” Jop. J. Appl. Phys, vol. 14, pp. 1921-1927, 1975. -au‘ a--. ”4.1.- -..‘Tu 160 [20] R. McGowan and R. Kuc, ”A direct relation between a signal time series and its unwrapped phase,” IEEE Trans. on A.S.S.P., vol. ASSP-30, pp. 719-726, Oct. 1982. [21] H. Hochstadt, Integral equations, New York: John Wiley & Sons, Inc., 1973. [22] D. Colton and R. Krees, Integral equation method in scattering theory, New York: John Wiley & Sons, 1983. [23] C. A. Balanis, Antenna theory: analysis and design, New York: John Wiley & Sons, 197 3. [24] R. G. Newton, Scattering theory of waves and particles, New York: McGraw- Hill, 1966. [25] M. Kareh, M. Soumekh and J. F. Greenleaf, "Signal processing for diffraction tomography,” IEEE Trans. on Sonics and ultrasonics, vol. 80-31, pp. 230-238, July 1984. ' [26] E. M. Kennaugh and D. I... Moffat, "Transient and impulse response approxima- tions,” Proc. IEEE, pp. 893-901, Aug. 1965. [27] L. Marin, ”Representation of transient scattered fields in terms of free oscillations of bodies,” Proc. IEEE, pp. 640-641, May 1972. [28] E. Mese, et. al., "Target identification by mean of Radar," Microwave Journal, pp. 85-102, Dec. 1984. [29] P. H. Johnston and E. I. Madaras, ”Modeling the pulse-echo response of a two- dimensional phase-insensitivity array for NDE of layered media,” in Proc. IEEE Ultrason. Symp., Chicago, IL, pp 1021-1025, 1988. [30] T. J. Cavicchi, S. A. Johnson, and W. D. O’Brien, ”Application of the sinc basis moment method to the reconstruction of infinite circular cylinders,” IEEE Trans. on U.F.F.C., vol. 35, no. 1, pp. 22-33, Jan. 1988. i ‘4' 161 [31] C. Gazanhes, J. P. Herault, and K. Stephanakis, "Reflection coefficients identification by mean of correlation: application to layered medium,” J. Acoust. Soc. Amen, vol. 69, no. 3, pp. 720-727, Mar. 1981. [32] Y. Lu, J. D. Achenbach, ”Effects of random deviations in interface properties on the propagation of ultrasound in thick composites," J. Acoust. Soc. Amen, vol. 90, no. 5, pp. 2576-2585, Nov. 1991. [33] M. Cookey, H. J. Trussell, and I. J. Won, "Seismic deconvolution by multiple methods,” IEEE Trans. on A.S.S.P., vol. 38, no. 1, pp. 156-160, Jan. 1990. [34] R. Yarlagadda, J. B. Bednar, and T. L. Watt, ”Fast algorithms for L, deconvolu- tion,” IEEE Trans. on A.S.S.P., vol. 33, no. 1, pp. 174-181, Jan. 1985. [35] S. A. Goss, R. L. Johnson, and F. Dunn, "Comprehensive compilation of empiri- cal ultrasonic properties of mammalian tissue,” J. Acoust. Soc. Amen, vol. 64, no. 2, pp. 423-457, Feb. 1978. [36] R. Kuc, ”Clinical application of an ultrasound attenuation coefficient estimation technique for liver pathology characterization," IEEE Trans. on Biomed. Eng, - vol. BME—27, pp 312-319, June 1980. [37] K. Sobezyk, Stochastic wave propagation, Warsaw, Poland: Elsevier, 1985. [38] K. Aki and P. Richards, Quantitative seismology: theory and methods, San Fran- cisco: W. H. Freeman, 1980. [39] F. R. DiNapoli and R. L. Deavenport, "Theoretical and numerical Green’s func- tion field solution in a plane multilayered medium," J. Acoust. Soc. Amen, vol. 67, no. 1, pp. 92-105, Jan. 1980. [40] J. W. Goodman, Introduction to Fourier optics, 3rd ed., San Francisco: McGraw- Hill, 1968. [41] T. K. Sarkar, D. D. Weiner, and V. K. Jain, ”Some mathematical considerations in dealing with the inverse problems," IEEE Trans. on Antennas Propagat, vol. 162 AP-19, pp. 373-379, Mar. 1981. [42] Y. Das and W. M. Boemer, 'On Radar target shape estimation using algorithms for reconstruction from projections," IEEE Trans. on Antennas Propagat, vol. AP-26, pp. 274-279, Mar 1978. [43] H. J. Schmitt, et al., ”Calculated and experimental response of thin cylindrical antennas to pulse excitation," IEEE Trans. on Antennas Propagat, vol. AP-14, pp. 120-127, Mar 1966. [44] C. 1.. Bennett, "Transient scattering from conducting cylinders,” IEEE Trans. on Antennas Propagat, vol. AP-l8, pp. 627-633, Sep. 1970. [45] L. B. Felsen, “Comments on early time SEM,” IEEE Trans. on Antennas Pro- pagat, vol. AP-33, pp. 118-119, Jan. 1985. [46] D. G. Dudley, ”Comments on SEM and parametric inverse problem," IEEE Trans. on Antennas Propagat, vol. AP-33, pp. 119-120, Jan. 1985. [47] M. A. Morgan, "Response to comments regarding SEM representations,” IEEE Trans. on Antennas Propagat, vol. AP-33, pp. 120, Jan. 1985. [48] K. M. Chen and D. Westrnoreland, "Impulse response of a conducting cylinder based on SEM," Proc. IEEE, vol. 69, pp. 747-750, Jan. 1981. [49] K. M. Chen, "Radar waveform synthesis for single-mode scattering by a thin cylinder and application for target discrimination,” IEEE Trans. on Antennas Pro- pagat, vol. AP-30, pp. 867-880, Sep. 1982. [50] R. Kuc, Introduction to digital signal processing, New York: McGraw-Hill, 1988. [51] H. P. Bucket, “Comparison of FFI' and Prony algorithm for hearing estimation of narrow-band signals in realistic ocean environment,” J. Acoust. Soc. Amen, vol. 61, pp. 756-762, Mar. 1977. [52] M. A. Rahmanand and K. B. Yu, ”Improved frequency estimation using total least square approach," in Proc. IEEE Int. Conf., Tokyo, Japan, pp. 1392-1400, 163 Apr. 1986. [53] M. P. Hurst and R. Mirna, ”Scattering center analysis via Prony’s method," IEEE Trans. on Antennas Propagat, vol. AP-35, pp. 986-988, Aug. 1987. [54] R. Can-idle and R. L. Moses, “High resolution Radar target modeling using a modified Prony estimator," IEEE Trans. on Antennas Propagat, vol. AP-40, no. 1, pp. 13-18, Jan. 1992. [55] R. Kumaresan and D. W. Tufts, "Estimating the parameters of exponential damped sinusoids and pole-zero modeling in noise,” IEEE Trans. on A.S.S.P., vol. ASSP-30, pp. 833-840, Dec. 1982. [56] I. Ziskind and M. Wax, ”Maximum likelihood localization of multiple source by alternating projection,” IEEE Trans. on A.S.S.P., vol. ASSP-36, no. 10, pp. 1553- 1560, Oct. 1988. [57] H. Yamada, M. Ohmiya, Y. Ogawa, and I. Itoh, "Superresoultion techniques for time-domain measurements with a network analyzer,” IEEE Trans. on A.S.S.P., vol. ASSP-39, no. 2, pp. 177-183, Feb. 1991. [58] Y. Hua and T. K. Sarkar, "Matrix pencil method for estimating parameter of exponentially damped/undamped sinusoids in noise,” IEEE Trans. on A.S.S.P., vol. ASSP-38, no. 5, pp. 814-824, May 1990. [59] Z. A. Mariéevié, T. K. Sarkar, Y. Hua, and A. R. Djordjevic, ”Time-domain measurements with the Hewlett-Packard network analyzer - HP 8510 using the matrix pencil method," Microwave Theory and Technique, vol. 39, no. 3, pp. 538-547, May 1991. [60] de Prony, Baron (Gaspard Riche), "Essai experimental et analytique, etc.," Paris J. of L'Ecole Polytech, Vol. 1, No. 2, pp. 24-76, 1795. [61] S. L. Marple, Jr., ”Spectral line analysis by Pisarenko and Prony methods,” in Proc. of 1979 IEEE Int. Confi, Washington, DC, pp. 155-161, 1979. 164 [62] R. N. McDonough and W. H. Haggins, ”Best least square representation of sig- nals by exponentials,” IEEE Trans. on Autom. Control, vol. AC-13, pp. 408-412, Aug. 1968. [63] S. I... Marple, Jr., Digital spectral analysis with application, Englewood Cliffs, NJ: Prentice-Hall, 1987. [64] T. J. Shari, M. Wax, and T. Kailath, "On spatial smoothing for direction-of-arrival estimation of coherent signals,” IEEE Trans. on A.S.S.P., vol. ASSP-33, no. 4, pp. 806-811, Apr. 1985. [65] D. W. Tufts and R. Kumaresan, ”Estimation of frequencies of multiple sinusoids making linear prediction perform like maximum likelihood,” in Prac. of 1982 IEEE Int. Conf, vol. ASSP-70, no. 9, pp. 975-989, Sep. 1982. [66] W. F. Gabriel, ”Using spectral estimation techniques in adaptive processing antenna systems," NRL-REP. 8920 Naval Res. Lab., Washington D.C., Oct. 1985. [67] R. O. Schmidt, ”Multiple emitter location and signal parameters estimation," IEEE Trans. on Antennas Propagat, vol. AP-34, pp. 276-280, Mar. 1986. . [68] D. W. Tufts and R. Kumeresan, “Improve spectral resolution,” in Prac. of IEEE, vol. 68, no. 3, pp. 419-420, Mar. 1980. [69] W. H. Press, et. al., Numerical recipes, the art of scientific computing, New York: Cambridge University Press, 1986. [70] R. K. Otnes and L. Enochson, Applied time series analysis, New York: John Wiley & Sons, 1978. [71] M. Zoltowski and F. Haber, ”A vector space approach to direction finding in a coherent multi-path environment,” IEEE Trans. on Antennas Propagat, vol. AP- 34. pp. 1069-1079, Sep. 1986. [72] R. T. Williams, S. Prasad, A. K. Mahaland, and L. H. Sibyl, "An improved spa- tial smoothing technique for hearing estimation in a multi-path environment," 165 IEEE Trans. on A.S.S.P., vol. 36, pp.425-432, Apr. 1988. . [73] G. E. Forsythe, M. A. Malcolm, and C. B. Molar, Computer methods for mathematical computations, Englewood Cliffs, NJ: Prentice-Hall, 1977. [74] A. W. Bojanczyk, M. Ewerbring, F. T. Luk, and P. V. Dooren, "An accurate pro- duct SVD algorithm," in SVD and Signal Processing II: algorithms, analysis and applications, R. Vaccaro, ed., Amsterdam: Elsevier Science Publisher, 1991. [75] R. Kumaresan and D. W. Tufts, ”Singular value decomposition and spectral analysis,” Prac. IEEE, pp. 6.4.1-6.4.12, Aug. 1981. [76] C. Lanczos, Applied numerical analysis, Englewood Cliffs, NJ: Prentice-Hall, 1956. [77] V. C. Klema and A. T. Lamb, "The singular value decomposition : its computa- tion and some applications,” IEEE Trans. Autom. Control, AC-25, pp. 164-176, Apr. 1980. [78] C. Y. King and B. Ho, "Simultaneous measurement of material characteristics of layered structures by a single acoustic interrogation,” IEEE Trans. an Instrument and Measurement, Oct. 1993. (forthcoming) [79] S. Haykin, ED., Advances in spectrum analysis and array processing, Vol. I and II, Englewood Cliffs, NJ : Prentice-Hall, 1991. [80] R. Gagliard, Introduction to communication engineering, New York: John Wiley & Sons, 197 8 [81] P. P. Lele and K. J. Parker, ”Temperature distributions in tissues during local hyperthermia by stationary or steered beams of unfocussed or focussed ultrasound,” Brit. J. Cancer Suppl. , vol. 45, pp. 108-121, 1982. [82] AIUM (The American Institute of Ultrasound in Medicine), “Bio-effects con- siderations for the safety of diagnostic ultrasound,” J. Ultrasound Med, vol. 7, no. 9, pp. 1-38, 1988. .t‘r... ‘ Am _n'am‘ - ...I. n .a-aF 166 [83] W. L. Lin, R. B. Roerner, and K. Hynynen, ”Theoretical and experimental evalua- tion of a temperature controller for scanned focussed ultrasound hyperthermia,” Med Phys, vol. 17, pp. 615-625, 1990. [84] B. A. Auld, Acoustic field and waves in solids: Vol. II, New York: John Wiley & Sons, 1973. [85] G. Beylkin and R. Burridge, ”Linearized inverse scattering problems in acoustic and elasticity,“ Wave Motion, vol. 12, pp. 15-52, 1990. [86] L. M. Brekhovskikh, Waves in layered media, translated by R. T. Beyer, New York: Academic press, 1980. [87] D. E. Bray and R. K. Stanley, Nondestructive evaluation: A tool for design, manufacturing, and service, New York: McGraw-Hill, 1989. [88] R. C. Wang, J. P. Astheimer, and G. W. Swartout, ”A characterization of wave- front distortion for analysis of ultrasound diffraction measurement made through an inhomogeneous medium,“ IEEE Trans. on Sonic & Ultrasan, vol. 32, no. 1, pp. 30-48, Jan. 1985. [89] S. M. Rao and B. S. Sridhara, ”Acoustic scattering from arbitrarily shaped multi- ple bodies in half space: method of moment solution,” J. Acoust. Soc. Amen, vol. 91, no. 2. DP. 652-657, Feb. 1992. [90] Y. L. Li, C. H. Liu, and S. J. Franke, ”Three dimensional Green’s function for wave pr0pagation in a linearly inhomogeneous medium - the exact analytic solu- tion," J. Acoust. Soc. Amen, vol. 87, no. 6, pp. 2285-2291, June, 1990. [91] E. A. Kraut, “Review of theories of scattering of elastic waves by cracks," IEEE Trans. Sonics & Ultrasan, vol. 23, no. 3, pp. 162-167, May 1976. [92] J. E. Gubernatis, E. Domany, J. A. Krumhansi, and M. Huberman, "Formal aspects of the theory of the scattering of ultrasound by flaws in elastic material,” J. Appl. Phys, vol. 48, no. 7, pp. 2804-2811, 1977. 167 O [93] J. E. Gubematis, E. Domany, J. A. Krumhansi, and M. Huberman, "The Born approximation in the theory of the scattering of elastic waves by flaws,” J. Appl. Phys, vol. 48, no. 7, pp. 2812-2819, 1977. [94] J. U. Tribolet. Seismic application of hamamarphic signal processing, Englewood ~ Cliffs, NJ: Prentice Hall, 1979. [95] G. S. Kine, Acoustic waves: devices, imaging and analog signal processing, Englewood Cliffs, NJ: Prentice Hall, 1987. [96] M. Bouchon, “Calculation of complete seismograrns for an explosive source in a layered medium," Geaphys, vol. 45, no. 2, pp. 197-203, Feb. 1980. [97] E. J. Ayme-Bellegarda and T. M. Habashy, “Forward ultrasonic scattering from multidimensional solid or fluid inclusions buried in multilayered elastic struc- tures," IEEE Trans. on U.F.F.C., vol. 39, no. 1, pp. 1-10, Jan. 1992. [98] E. J. Ayme-Bellegarda and T. M. Habashy, ”Ultrasonic inverse scattering mul- tidimensional object buried in multilayered elastic background structures,“ IEEE Trans. on U.F.F.C., vol. 39, no. 1, pp. 11-18, Jan. 1992. [99] E. J. Ayme-Bellegarda and T. M. Habashy, "Ultrasonic inverse scattering of 2-D objects buried in multilayered elastic background structures," Proc. IEEE I990 Ultrasan. Symp., pp. 845-848, Dec. 1990. [100] C. Y. King and B. Ho, "The ultrasonic field within an inhomogeneous and arbi- trarily shaped scatterer," IEEE Trans. on U.F.F.C., (under review) [101] B. A. Auld, Acoustic field and waves in solids: vol. I, New York: John Wiley & Sons, 1973. [102] D. Ensminger, Ultrasonics - fiendamentals, technology, and applications, New York: Marcel Dekker, 198 8. [103] E. J. Ayme-Bellegarda and T. M. Habashy, ”Constrained least-squares recon- struction of multidimensional objects buried in inhomogeneous elastic media,” in 168 Proc. IEEE 1991 Int. Conf. A.S.S.P., pp. 2502-2508, May 1991. [104] S. M. Rao and B. S. Sridhara, ”Application of the method of moments to acous- tic scattering from arbitrarily shaped rigid models coated with lossless, shear-less materials of arbitrary thickness,” J. Acoust. Soc. Amen, vol. 90, no. 3, pp. 1601- 1608, Mar. 1990. [105] Edited by R. Vaccaro, SVD and signal processing, II: algorithm, analysis and applications, Amsterdam: Elsevier Science publisher, 1991. [106] S. A. Goss, R. L. Johnson, and F. Dunn, ”Compilation of empirical ultrasonic properties of mammalian tissues, II," J. Acoust. Soc. Amen, vol. 68, pp. 93-108, 1980. [107] X. B. Fan and K. Hynynen, "The effect of wave reflection and refraction at soft tissue interfaces during ultrasound hypertherrnia treatments," J. Acoust. Soc. Amen, vol. 91, no. 3, pp. 1727-1736, Mar. 1992. [108] T. W. Parks and C. S. Burrus, Digital filter design, New York: John Wiley & Sons, 1981. [109] D. Ensminger, Ultrasonics: firndamental, technology, and applications, 2nd ed., New York: Marcel Dekker, 1988. [110] G. S. Kino, “The application of reciprocity theory to the scattering of acoustic waves by flaws,” J. Appl. Phys, vol. 49, no. 6, pp. 3190-3199, 1978. "‘uuuuuu