Reconstruction of the acoustical impedance profile of a multilayer medium

ULTRASONIC

IMaGING

14,

40-68

RECONSTRUCTION

(1992)

OF THE ACOUSTICAL IMPEDANCE MULTILAYERMEDIUM

PROFILE OF A

Ran Lifshitx. Peder C. Pedersenl and Peter A. Lewin2 Cardiometrics Inc. 645 Clyde Avenue Mountain View, CA 94043

The recovery of the acoustical reflectivity function and impedance profile of a layered medium from bandliited and noisy pulse-echo ultmsonic data is considered. The effects of the transducer and the noise am reduced using Wiener filtering. With further signal processing and a priori knowledge of the attenuation and the velocity profiles, the compensated coefficients of the reflectivity function and the impedance profile are recovered. The reconstruction techniques are presented analytically, and are also evaluated in an experimental setup composed of a conventional pulse echo system, a data acquisition and a data processing system. Under experimental conditions, the RMS errors in the estimation of the compensated acoustical discrete impedance profile were within 3% of the calculated values. 0 1992 Academic Press, Inc. Key words: Deconvolution, impedance profile, inverse methods, medical imaging.

I. INTRODUCTION The majority of acoustic imaging techniques available today utilizes the energy backscattered from the inhomogeneities in the medium. The procedure of using the scattered wave to reconstruct the physical properties of an unknown medium is usually referred to as the inverse problem. Acoustical inverse techniques are widely employed in many fields, such as geophysical signal processing [l] and medical imaging [2]. Inverse techniques ate also applied in the field of ultrasound tissue characterization. In particular, the implementation of these techniques to reconstruct the acoustic impedance protile of biological tissue is of interest. This is because the knowledge of the tissue acoustical impedance profile along with attenuation and velocity profiles provides a complete acoustic characterixation of the tested tissue medium. This technique, if sufficiently accurate, can be applicable in diffbrentiating between normal tissue and tissue with different pathologies.

1 Department of Electrical Engineering, Worcester Polytechnic Inst., Worcester, MA 01609. 2 Department of Electrical and Computer Engineering, Drexel University, Philadelphia, PA 19104. 0161-7346/92 $3.00 Copyright Q 1992 by Academic Press, inc. All rights of reproduction in any form reserved.

40

IMPEDANCE

PROFILE RECONSTRUCTION

In a pulse echo mode of ultrasound imaging, pulses are transmitted and received with the same transducer. The reflected signal, y(t), is used to determine various properties of the medium. A characterization in terms of velocity, impedance and absorption fluctuation [3 - 51 may be performed. In order to determine one specific acoustic parameter, it is generally necessaryto make apriori assumptions about several other acoustic parameters. Solutions to the inverse problem are given for specific geometries and measurement conditions. For acoustic impedance profiling, the assumption of a stratified medium allows a relatively simple formulation. Furthermore, this assumption is appropriate in several applications in medical and industrial ultrasound, and the basic mathematical treatment of the inverse acoustical problem may then be formulated [6] as a convolution performed in one dimensional space. In this case, y(t), the received signal from of the ultrasonic system is given by: y(t) = x(t) * h(t) * e(t) + n(t)

(1)

where x(t) is the electrical excitation function, h(t) is the combined transmitting and receiving transducer impulse response and n(t) corresponds to the measurement noise. The symbol * indicates convolution. The term e(t) is the reflection impulse response of the layered medium under investigation and is determined by the acoustic impedance profile (only interfaces have reflection coefficients) of the medium in the direction of propagation. If the medium consists of lossless, homogeneous layers, e(t) will be a sum of impulses whose amplitudes are determined by the relative change in acoustic impedance between adjacent layers. In most cases, a nonnegligible attenuation exists which modifies the amplitudes of the pulse in e(t). From e(t), with the attenuation effect removed, one can in principle determine the impedance profile. Thus, clearly, a first requirement for recovering the impedance profile is to estimate the medium impulse response. If the time durations of x(t) and h(t) am short and, furthermore, the signal-to-noise ratio is high, then y(t) - e(t). However, in practice, this is not the case. The duration of the signal emitted by the transducers cannot be ignored. Also, noise is always present in the real measurement system and cannot be neglected. Hence, there is an uncertainty regarding the location of interfaces and closely spaced interfaces may not be resolved. Nonetheless, in many cases, an appropriate processing of the received signal allows the interface echoes to be identified and quantified. Given exact knowledge of the bandlimited transducer impulse response h(t) and the statistical properties of the noise n(t), an approximate determination of e(t) can be obtained by means of various deconvolution techniques. One very common deconvolution method is the Wiener filtering technique [8 - lo] which gives the least mean square amplitude error over the transducer passband. This filter is an optimum noise filter which combines a minimum signal distortion and a maximum reduction of the noise. This restoration filter approaches a pure inverse filter l/X-I(o) if the noise power is small and performs as a matched filter when the signal-to-noise ratio is poor. The Wiener filter is relatively easy to implement in software and may also be hardware implemented by means of a storagecorrelator 111, 121. Another deconvolution technique utilizes the matched filter approach [7] which involves the correlation of h(t), the transducer impulse response, with the data y(t). This filter is optimized if the noise is white and Gaussian.

41

LIFSHITZ ET AL

Time domain deconvolution techniques have been investigated by a French group [ 16 181 and other investigators [13 - 151. While in most other available techniques the processing is done after recording the data, the majority of the time domain techniques can in principle be carried out in real time. The main disadvantage is still due to the complexity of the time domain algorithms when compared to the frequency domain techniques. The French group reported good experimental results such as collecting real signals in vitro in a rabbit eye and detecting the variation of the internal diameter and wall thickness of the humeral artery. Another important class of deconvolution techniques is referred to as the “homomorphic method”. This method combines both time and frequency approaches. Homomorphic processing [19], or more specifically, cepsaal analysis, is useful primarily because it allows the removal of a convolution term in the time domain or a product term in the frequency domain. This method has been found useful especially when e(t) has a periodic structure. Yet another class of deconvolution technique, recently suggested by Schmolke et al [20] is the prior inverse filtering method. In this method, the transducer is excited by a synthesized signal which has the spectrum of the inverse filter l/H(o) within specified frequency limits. While this technique essentially produces the same results as the conventional (posterior) inverse filtering, it has several advantages such as simpler implementation, real time capability, greater flexibility and better signal to noise ratio (SNR). A great deal of work has been done in the past few years in studying continuously varying and discrete impedance profiles [21 - 281. An algorithm that is accurate for lossless media with small variations in acoustical impedance, but not strictly applicable for layered media [26,27] is known as the impediography equation : z(t) = q, exp[ jute(r)dz] where zo usually denotes the acoustical impedance of a coupling medium and z(t) is the acoustic impedance at point x = 42 where t is the propagation time. Several papers [29, 301 have been published addressing specifically the layered medium reconstruction problem. In these papers, e(t) is estimated by extrapolating H(w), the transducer bandlimited frequency response beyond its frequency range. In this paper, the goal is to formulate, implement and experimentally test practical and easy-to-use algorithms to determine the one dimensional impedance profile of a multilayer medium of finite thickness. The emphasis is on recovering the medium impulse response (reflectivity function). Once the medium impulse response is obtained and attenuation compensation has been applied, the impedance profile can be derived. The Wiener filter was used extensively in this study due to its relative simplicity and accuracy. In addition to Wiener filtering, algorithms were developed to determine the best match of the received signal to the true reflectivity function. Two different methods for estimating the reflection coefficients of the layers are evaluated. These methods permit detection of very small variations in the impedance from one layer to the other.

42

IMPEDANCE

PROFILE RECONSTRUCTION

This paper is organized as follows: The theory of the recovery procedure is discussed in detail. First, a model of the measurement system is formulated. Then, the deconvolution approach to correct for the distortion of the actual medium reflection impulse response by the measurement system is discussed. The procedure of the recovery of the reflectivity coefficients from the deconvolved signal, or directly from the measured data, is then presented. By compensating for the transmission and attenuation losses, the reflection coefficients are obtained from the previously estimated reflectivity function. The experimental setup and the overall experimental results and data analysis are finally discussed. II. THEORY In this section, the theory describing the recovery of the acoustical impedance profile of a multilayer medium by means of reflection measurements is discussed. Below are given the general assumptions on which the theory is based: (i) The measurement process is based on modeling the backscattered echo waveform in terms of a layered (one dimensional) medium. (ii) Both the signal generating process and the interaction between the acoustical wave form and the medium are linear processes. Therefore, deconvolution techniques can be applied. (iii) The interrogating wave is longitudinal plane wave. (iv) The propagating ultrasonic wave is assumed to be only moderately reduced in amplitude as it travels through the medium due to attenuation effects. (v) The variations in acoustic impedance are assumed to be relatively small, so that transmission losses and multiple reflections can be ignored. (vi) The transducer is inherently bandlimited. (vii) Any type of noise in the system is additive. Assumption (i) implies that the unknown medium (e.g., biological tissue) is stratified in plane parallel layers, and that the axis of the acoustic beam is normal to these layers. Assumption (ii) states that if both processes are linear, then the reflected signal is a convolution of the transmitted signal with the reflectivity function representing the medium. It is also assumed that the conversion of the ultrasonic wave into analog electrical signal is linear. Assumption (iii) implies that the sound field in the region of interest can be approximated to a plane wave field; this assumption holds in focal region of focused transducers and - in limited regions - in the far field of unfocused transducers. Assumption (iv) states that over some (limited) range of the unknown medium, the attenuation effect is sufficiently small so that the acoustic pulse is changed in amplitude only, but not in spectral content. Therefore, only a scalar amplitude correction is required to generate e(t) for the corresponding lossless medium. In this paper, small variation of the impedance within the tested medium are of interest. Small variation of the impedance implies very small pressure reflection coefficients which justifies assumption (v). Assumptions (iv) and (v) will be further discussed later in this paper. The assumption of a bandlimited transducer (assumption (vi)), is realistic since in fact, all available ceramic piezoelectric transducers show resonant behavior. (In contrast, polymer (PVDF) transducers are wideband systems [31] but their sensitivity is significantly lower than that of PZT transducers). The last assumption (vii) does not impose any condition on the noise type. That means that the noise can be either white, or colored.

43

LIFSHITZ ET AL

Model of the reflection response of the medium and the measurement system Consider the following simple pulse echo experiment. A very short excitation pulse is applied to the transducer, resulting in an acoustic pressure h,(t) at a given field point. h,(t) can be considered a good approximation to the transmitting transducer impulse response for that field point (this will generaliy be the time derivative of the velocity response on the surface of the transducer). The acoustic pressure pulse is then modified by propagation through and reflection from the insonitied medium, i.e., h,,(t) is convolved with e(t), the medium reflection impulse response. (In this paper, the function e(t) is referred to as the ‘reflectivity function’). Corrupted further by additive noise, the reflected acoustic signal is detected by the receiving transducer. The detected signal is convolved with h,(t), the receiving transducer impulse response. Hence, the measured signal y(t), received from the insonified medium, is a distorted version of e(t). One can therefore write the output from the measurement system, y(t), as y(t) = [&t(t) * e(t)+ n&l ) * h,(t) J + n,(t) The term “a(t) is usually referred to as the acoustical noise, while n,Jt) usually corresponds to the electrical noise. Assuming: h,(t) = h,(t) and h,(t) * h,(t) = h(t), one can write y(t) as y(t) = h(t) * e(t) + h,(t) * n,(t) + n,(t) = h(t) * e(t) + n(t)

(3)

with n(t) = h,(t) * n,(t) + n,(t) being the overall noise in the system. In the frequency domain, the output model Y(w) is given by Y(o) = H(o) E(o) + N(o)

(4)

where Y(o) = Fly(t)); H(w) = F(h(t)); E(o) = F(e(t)); and N(o) = F(n(t)) and where F[) represents Fourier transform. The above derivations and analysis were made with no restriction regarding the attenuation in the medium. Deconvolution and Wiener filtering Deconvolution or inverse filtering attempts to restore the generally broadband spectrum of e(t), by compensating for the spectral distortion of the transmitting and receiving transducers. If the strict inverse filter is denoted by R(o) = l/H(o), then applying R(w) on Y(o) one gets the estimation of E(o), I?(o), as &(a) = Y(o)R(w)

=#

Difficulties arise for those frequencies where H(o) approaches zero. In order to compensate for the bandlimited character of H(o), the inverse filter amplifies low power frequency components relatively to the high power frequencies. Consequently, the frequency cxmponents of the noise outside the system passband will be amplified very strongly by the inverse filter, and the noise will dominate the deconvolved signals in the time domain completely. The optimum filter, which combines a minimum signal distortion and a maximum reduction of the

44

IMPEDANCE

PROFILE RECONSTRUCTION

noise, is the Wiener filter [lo]. This optimal filter is given by:

w(w) =A-[ ,H(;g(“)]

(5)

where H(o) is the transducer transfer function. The function G(o) in Eq. (5) is defined as: cN(w)N*(w)> E(@)E*(w)

G(o) =

where the terms N(w) and E(o) are defined in Eq. (4), c > denotes the expected value and the superscripted * means complex conjugate. In most practical applications, the statistical properties of the processes involved are unknown, approximated by a constant $ whereby Eq. (5) becomes:

w(0)

1 = H(w)

[

IH( o) I2 IH(&+ 4

1

and in those cases G(o)

= R(o)lWs(w)l

can be

(6)

where IW,(w)l = IH(o)12/ ( IH( + o). This version of the Wiener filter is known as the Modified Wiener filter. It should be pointed out that this last version is theoretically true only if both n(t) and e(t) are white processes, in which case G(o) becomes $, the noise to signal power density ratio (NSR). This restoration filter approaches to a pure inverse filter if the noise power is small compared to the signal power. On the other hand, W(o) turns out to be a matched filter when the SNR is poor. In this case, the filter suppresses unwanted signals outside the system frequency band and can induce an improvement of the signal to noise ratio. A model for the reflectivity function In this approach, the problem of reconstructing e(t), the medium impulse response (reflectivity function) is considered. The medium consists of finite number of plane parallel layers, as depicted in figure 1. The impedance profile of the medium, z(t), is estimated from e(t). Therefore, the reflectivity function e(t) is a sum of impulses and is given by: M

e(t) = c

(7)

rkmk6 ( t-tk)

k=O

where M is the total number of layers in the medium and is generally unknown. In this paper, only moderate variations in acoustic impedance from layer to layer are considered, and hence muItiple reflections can be ignored and are not included in Eq. (7). The time value tk corresponds to the round trip travel time between the transducer and the boundary between the kth and the (k+l)* layer. Similarly, rk is the true reflection coefficient between the kth and the (k+l)* COeffiCient

layer, which in general can be assumed to be frequency independent. The reflection rk iS modified by an attenuation faCtOr, mk, which is the nXX3n-frequenCy attenuation

45

LIFSHITZ ET AL

0th

(k-l)st in-

1st intelfacc

interface

TIIlWlitting/ Receiving

kst

TranSdllC~

II

zo

0th layer

z1

1stlayer

‘k

4-l

.

.

.

.

k-th layer

‘k+l

(k+l)st layer * 0 l

Fig. 1 Layered mediummodel.

of the kth layer. Defining ak = rkmk

(8)

gives M

e(t) = c

(9)

aks (t-tk)

k=O

Basedon this model the searchfor e(t) is reducedto determiningthe ak and t, valuesof Eq. (9). Measurementmodel The measureddata, y(t) is defined in Eq. (3). Making the substitutionof Eq. (9) in (3), onegets:

y(t) = h(t) * e(t) + n(t) = h(t) *

+ n(t) =

Thus, the output signaly(t) is a train of transducerimpulseresponses h(t), scaledby the factor ak, delayedby 4(and corruptedby the additive noise.

46

IMPEDANCE

PROFILE RECONSTRUCTION

Recovery of the reflectivity function Fourier transforming Eq. (lo), one gets: M

Y(O) = cakH(o)exp(-jotk)

+ N(o)

(11)

k=O

Applying the Modified Wiener filter with the white noise assumption on the data Y(W) one attains an estimate of the reflectivity function, E( 0):

[

^ 1 E(w) = ‘(@H(w)

IH( ,H@$+

@

1

Y(o) = -lws(o)l H(o)

ln the equation above, IW,(o)l represents the gain function of the Wiener filter which is nearly unity for frequencies where the signal to noise ratio is high and which drops to near zero for frequencies where the signal power is vanishing. Since the Wiener filtering is applied to the reflections from each layer interface, it is appropriate to define a new function for each interface which consists of IWs(cu)l scaled by the reflection coefficient for that interface and modified by the attenuation of the layer in front of it. Calling this new function Bk(o), follows:

it is defined as

Bk(o) = IWs(w)l ak = lWs(w)l rk mk Introducing Y(w) from Eq. (11) into the expression for E(o) and applying the definition for Bk(0), a new expression

for

E(0)

iS

obtained:

IWs(@> 1 [ 1No)

M

B(W) = Cak

exp{-j 0 tk} + Ho

IWs(w)l

k=O M =

By c k=O

where N&o) = [ IWs(o)l/H(w) response estimation is given by

eXP{-j

O.J tk}

+ &(a)

(12)

] N(o) is a noise term. The time domain medium impulse

qt)=F-l@(W)]=fbk(t-tk)

+n,(t)

(13)

k=O

where b(t) = F-1 (B&0)) and bk(t-tk) = F-l { Bk(CO)eXp {-jwtk) ). One might expect the reflection response of an interface to be represented by very narrow spikes scaled by the ak coefficients as in Eq. (9). However, the actual experiments reveal narrow, but finite duration signals resembling sine functions. This behavior becomes clear considering Eq. (12). The first term in Eq. (12) represents the coefficients ak, multiplied by a wide bandpass filter IW,(w)l,

47

LIFSHIT’Z ET AL

and a phase term corresponds to the delay in the time domain. Therefore, the signals bl&t-tk) = F-I{ IWs(o)l eXp(-jatkak)) will have a finite duration. Hence, the sine function behavior is due to the missing (low) frequencies caused by the transducer filtering effect. Nonetheless, deconvolution by means of Wiener filtering enhances the resolution drastically at the expense of lower signal to noise ratio (SNR). Estimation of the ak coefficients: Maximum Amplitude Method The next step is to get an approximation for the original layered medium impulse response as in Eq. (9). In other words, generate a series of impulses that provide the best match to the experimentally derived function e(t) given in Eq. (13). Such a function would take the form, given in Eq. (14).

(14)

Note that %(t) given in Eq. (14) is stated as an approximationbecausee(t) in Eq. (13) and in Eq. (14) are not identical. e(t) in (13) is a train of short pulses,resembbngsine functions, wherease(t) in Eq. (14) is a train of impulses.^akdenotesthe esrimaredcoefficient, in contrast to ak which denotesthe true coefficient. M+l is the numberof reflections correspondingto a mediumwith M layers.The arrival timestk andthe polarity of & coefficientsaredeterminedby applying a peak detection algorithm to e(t). The arrival times are defined as the times of maximumamplitudesof the reflection signalsin the output of the Wiener filter. It is important hereto note that the time location of the peak of a given pulsein the Wiener filter signal,e(t), correspondsclosely to the time location of the beginningof that pulsein the received signal, y(t) before Wiener filtering. Using Eqs. (13) and (12) one can write the following expression for the ratio of two consecutivereflections

maX[ bk( t-tk) 1

maX[

ws(t-tk>akl

max] bk-1(t-tk-111 = ImXIWs(t-tk-$k-Il

=

akmax[ wS( t-tk) 1 = ak-lmax[ws(t-tk-l)l

ak = ak_l

(1%

wheremax[ ] corresponds to the maximumpositive or negativeexcursion.The maximumvalue of the bk (t-tk) is chosensinceit is clear that the signalis leastcorruptedby the additive noise rid(t) at the maximumvalue. The last stepin Eq. (15) is true due to the fact that W,(t) is time invariant. For thoseregionswhere the individual reflectionsbl&t-tk) are separatedand much larger thanthe noisen&t). the following relationshold, asobtainedfrom Eq. (15):

ii1 = 20 &a= il

=[

WWI

-[bo(t-to -[ -h

)I

wt-t2)l (t-t113

m4

=’

Wt-t2)

max{be(t-tu)]

48

I

IMPEDANCE

PROFILE RECONSTRUCI’ION

This leads to an expression for the ^akcoefficients k= 1, 2,....M

(16)

Notice that & has to be determined prior to estimating the ^ai,coefficients. This can be done by acquiring a reference signal from a reflector with a known reflection coefficient and then comparing the first echo in the received signal y(t) with the reference echo. If the first echo in y(t) overlap with the next one, another technique such as the one described in the following section can be applied. Therefore, given ^&, the entire set (^ak) can be recovered. Estimation of the ak coefficients: Adaptive Method In this section, an alternative method for reconstructing the reflectivity function is presented. In this approach, the magnitude of the ak coefficients are determined directly from the received signal, y(t) whereas the arrival times and polarity of the ak coefficients am found from the Wiener filtered signal, e(t), just as with the Maximum Amplitude Method. Again, $ will denote the estimared coefficient, in contrast to ak which denotes the true coefficient. Basically, the Adaptive Method determines the ratio of the energy from a given pulse relative to the energy of a pulse, obtained from a perfect reflector. The energy from a given pulse is obtained by integrating the received signal, y(t) over an appropriate time segment. The system first determines the ?$kmagnitude for the reflection closest to the transducer, i.e., ^&; next i, magnitude is determined, etc. Since the energy of a given reflection can be expected to be determined more accurately when partially overlapping reflections have been removed, the contribution due to a given reflection is stripped from y(t) after the & value for this reflection is determined. Thus, the Adaptive Method sequentially stips y(t) of the contributions from the individual reflections, one at a time, such that the contribution from the reflection closest to the transducer is removed first. To formulate this mathematically, a new function, fk(t), is produced which consists of y(t) from which the k first reflections have been stripped. Conceptually, the Adaptive Method contains the following steps: I. Determine the arrival times, tk and polarity of reflections from interfaces, based on peak detection of C(t). II. Acquire an h(t) for an ideal reflector, i.e., the impulse response of the combined transmitting and receiving transducer measured with a perfect reflector. III. Determine the duration of h(t), called T. This value wiIl be used in setting the integration limits when determining the energy of a given pulse in y(t). IV. Find the magnitude of &,, the first ak coefficient, from energy considerations, as discussed below. V. Based on & and h(t), produce a new signal, called ye(t) = &h(t-t,,) represents the received signal due to reflection from the 0th interface.

49

9

which

LIFSHITZ ET AL

VI. Define a signal, fl(t) which is y(t) with the first reflection removed, i.e., f*(t) = y(t) -Y&h VII. Find the magnitude of ^a,, the second ^ak coefficient, considerations, by integration of an appropriate part of fl(t).

from energy

1

VIII. Create yltt) = c

4 h(t-ti)

which

represents

the received signal due to

i=O

reflections from the Oth and lSt interfaces. IX. Define f2(t) as y(t) with the two first reflections removed, i.e., fi(t) = y(t) - yl(t). X. Based on the relationships in VI to IX, the following general expressions can be defined: k fk+l

= Ytt)

- h(t);

f,(t)

yk(t)

= y(t);

= C%h( i=O

(17)

t-t;)

Note that fk+l(t) consists of y(t) from which the (k+l) first reflections, starting with yu(t), have been stripped. XI. The magnitude of subsequent ak coefficients are found by applying steps VII, VIII, and IX and incrementing k by 1. As described earlier, the arrival times tk and the polarity of the ^akcoefficients are determined by applying a peak detection algorithm to e(t). The arrival times are defined as the times of the maximum amplitudes of the reflection signals in the output of the Wiener filter which correspond to the beginning of the pulses in y(t). Thus, in absolute time, the k* pulse in y(t) lasts from tk to tk + T, and the energy contained in the pulse may be found by integrating over this time interval. If adjacent pulses are overlapping, only the non overlapping part d the received signal can be integrated without introducing an error. For overlapping pulses, tk+l - tk < T, and the integration limits must in this case be from tk to tk+I. The energy, received from a perfect reflector, is denoted En whereas the energy received from the k* interface is denoted Ek. When adjacent pulses do not overlap, E, and Ek are given as follows: E, = ;[h(t)12dt

;

E+T [fk(t)12dt

Ek =

(18)

‘k

When adjacent pulses do overlap, El, and E, take the following form: 4+1

E,=

tk+l-tk[ h( t)] 2dt; 0

E,=

[fk(t>12dt tk

50

(1%

IMPEDANCE

PROFILE RECONSTRUCIION

The magnitude of & is then calculated from EP and Ek: Ek + Ii&J= E

[1

(20)

P

Discussion of attenuation compensation One can define a general loss function Mj(O), characterizing layer j. This function may be written as a product of three independent mechanisms such that each one of them can be treat& separately : Mj(O) = Tj(O) Hj(W,x) exp[-2aj(w)xj]

(21)

where Tj(W) is the transmission loss factor, Hj(WJ) is the loss due to diffraction effects and exp[-2aj(U)xj] is the absorption and scattering factor. For the experimental data reported here, the diffraction effects are presumed to be negligible and Hj(O,x) therefore is assumed to be unity; in other words, plane wave insonification is postulated. It is also presumed that for layered structures the reflection coefficient set (rk) is frequency independent. Let P be the amplitude of the normal pressure wave insonifying the layered medium. If no absorption exists in the medium, the pressure wave coming from the kth boundary is simply k-l Pk = P rkn(

l-

rf>

i=l

In this paper, as well as for applications in soft tissue, a typical reflection coefficient is around 5% and will not exceed 20%. In this case, transmission loss effects can be neglected. To demonstrate it, consider a single layer, bounded by semi-infinite media, with pressure reflection coefficients IrOl= lrll = 0.2 (20%). The resulting pressure waves from the single layer (front and back wall echoes), are po = 0.2P and pl = 0.192P. The transmission loss is only 2Olog[pl/po] = 2Olog(O.96) = -.35dB. Therefore, Tj(W) is set equal to one. In general, the attenuation coefficient is a function of the frequency. In biological tissues, the attenuation coefficient a(x,o) can be approximated by a linear frequency dependence: a(x,o) = ma(x); and in many other materials by a quadrature frequency dependence: a(x,o) = &a(x). The attenuation due to absorption will have a low-pass filtering effect on the lossless waveform. As mentioned earlier, the assumption is that the attenuation effects are relatively small. Also, the ultrasound transducer is bandlimited and can in practice be chosen such that its center frequency is low enough (but not too low if a good resolution is still to be maintained), so that the absorption effect for each layer can be considered as a scalar attenuation coefficient that represents the average loss in the frequency band of the measuring system. One way to do this is to calculate an average attenuation factor for a given layer by

51

LIFSHI’I-Z ET AL

weighing the attenuation coefficient by the transducer frequency response, as shown in Eq. cm

cm

With this definition, approximated by

the loss due to absorption corresponding

to the kth layer can be

k mk=

rI j=l

eXP{-2&jXj}

;

rnc = 1 and k = 1, 2,....M

(23

where xj is the thickness of the j’s layer. Therefore, combining Eqs. (8) and (23) the reflection coefficients can be approximated by :

1 rk =

ma = 1, ac = ra and k = 1, 2,.... M k rI j=l

(24)

mk

eXP{-2&jXj}

Eq. (24) gives an iterative method to compensate for the attenuation loss, once the & coefficients are recovered. Reconstruction of the impedance profile Knowing the impulse response of the system, i.e., the coefficients of the reflectivity function, the impedance profile can be reconstructed in the following way: M z(t)

= zou(

t) +

c i=l

(25)

^zi U( t-t;)

where k = 1, 2.....M

(26)

and u(t) is the unit step function. In summary, to obtain the estimation of the impedance profile from a normal incidence pulse echo measurement one should apply the following steps assuming that the attenuation coefficients, ami the velocities ck = x&k - tk-l), k = 0,1,2,....M are known a priori. a. Measure the reflected echo, y(t).

52

IMPEDANCE

PROFILE RECONSTRUCTION

b. Apply deconvolution (Wiener filtering) and signal processing to get the & and the tk coefficients. c. Apply Eq. (24) to estimate the ?k coefficients. d. Determine the estimated impedance profile using Eq. (25) or alternatively Eq. (2). Finally, in the following example, it will be shown that the reconstruction of the impedance profile can be quite accurate even for a relatively large error in the estimation of the reflection coefficients. The above is true only as long as the reflection coefficients am relatively small (less than 0.1). For a medium, consisting of two layers with impedances zu and zl, the impedance of the second layer can be calculated from 20 and the reflection coefficient, r. l+r Zl=zOF

Next, a relative error dr is introduced to simulate the uncertainty in the estimation of the reflection coefficient r, such that f = r(l + dr). Thus: 2, = zo-l+i 1-i

=

l+r(l+dr) %l-r(l+dr)

The relative error in St, is termed e,, and with some algebra can be approximated by:

i, - Zl e=-= z Zl

--l+i y-i

l+r Z0l-r l+r Z0l-r

l+-

rdr l+r rdr

= l--i-q

-1

2r dr 5- l-r2

Notice that if r is small, say less than 0.2, then e, _= 2r dr. If the reflection coefficient is .05 and the relative error dr is as high as 20%, the corresponding error in z is only about 2%. III. EXPERIMENTAL

SETUP AND RESULTS

Instrumentation setup An instrumentation and system block diagram of the measurement system used for the experimental study is shown in figure 2. The pulser-receiver used to excite the transducer and to receive the reflected echoes was Panametrics model 5052 with 10 MHz bandwidth. The chosen transducer center frequencies ranged from 1.OMHz to 2.25 MHz. The digitizer was the L&my 9400 which acquires data at 100 MHz with 8 bit accuracy. Although the LeCroy 9400 can acquire up to 32k points, only 2048 points were actually transferred to the IBM to speed up the acquisition and processing time. This was sufficient for the thickness of our investigated media. In other words, reflections coming from any portion of the media were contained within the 2048 data points. The microcomputer used for controlling the measurement system, acquiring and storing the data, and for data analysis and signal processing was the IBM PC-AT. For both the data

53

LIFSHJZ

ET AL

WATERTANK

n e (0 I Electrical

1 noise

noise

h, (0

I Oy(t)

-+

0-t

Transmitting tI&IlSClUCCS

impulse response

MdiUm impulse response

Receiving transdu~r impulse response 0)

Fig. 2 (a) Instrumentationblock diagramof measurement system.(b) Systemblock diagram of measurement system.

acquisitionand signalprocessing,the ASYST softwarepackagewasused.Various programs were developedto perform data acquisition,simulations,analysisand signalprocessingof the data. Test objectdesign The medium is modeled here as a planar layered structure. The layers were implementedin the experiment using specially designedtest cells. The test cell is shown schematicallyin figure 3. The purposeof the testcell is to containvariousliquidsin the form of thin planarslabs,correspondingto the desiredmultilayer medium.The 8 cm square,fluid filled area was large enough to be consideredessentially of infinite extent when placed near the transmittingtransducer.The key to the designof the testcell wasthe useof shrink-wrapplastic (10 to 20 pm thick), which provided a thin yet strongmembrane.The plastic film washeId in placeusing3 mm thick metalframes,andwassealedto the stiff Plexiglasshellstructureusing rubber cement. Then the membraneswere carefully stretched taut using a heat gun. The

54

IMPEDANCE

PROFILE RECONSTRUCTION

Plexiglas structure

Fluid fill Tested fluid supported by two thin stnched plastic membranes

Metal fm

Fig. 3 Designof the testcell usedto modela planarlayeredmedium.

stretched membraneswere about 18-20 pm thick, basedon micrometer measurementsof disassembledtest cells. At 2.25 MHz, the wavelength of the soundtraveling through the membraneis approximately 1100pm (seereference[32] for the acousticpropertiesof Mylar). With a membranethicknessof 20 pm, and a wavelengthof 1100pm, the acoustictransmission coefficient is greater than 99.8%, and hence,the membranecan be consideredacoustically transparent. Once the membraneswere sealedand stretched,the test cell was carefully filled with liquid through the fill hole using a hypodermic needleand syringe. In order to equalize the hydrostatic pressureof the liquid againstthe membranesasthe cell was filled, the cells were slowly dippedinto the water tank during filling. Thiskept the membranes from bowing out; the fill hole was then sealedwith a small set screw. Several test cells were constructed with different thicknesses.The test cells were mounted on a two-axis gymbal manipulator and manuallead-screwslider. Ideally, the liquids usedin the test cells shouldhave acousticimpedancecloseto, but not the same,as the acousticimpedanceof water. Several different liquids were tested for propagationproperties,and thosethat exhibited acousticalimpedancewithin 20% of the water impedancewerechosen.The liquids includedmineraloil andPeg400 oil. The test cell, either a three layer or one layer configuration, was placed in the water tank. The insonifying nonfocusedtransducerwaslocated at the desired(axial) distanceaway from the test cell. Proper alignmentof the transduceris critical. This is becausethe impedance profile reconstructionprocedureis basedon the assumptionthat the interrogatingacousticpulse is at normalincidenceto the interfacesof the multilayer medium.This requirementis difficult to achieve in practice, due to the fact that it is virtually impossibleto built a multilayer test cell suchthat all interfacesareexactly parallel.

CT

LIPSHITZ ET AL

Table I Calculated parameters of the layered medium: mineral oil (3 mm) - water (3 mm) mineral oil (1 mm). Layer thicknesses were measured before injecting the fluids into the cells. Parameter

Description

Units

k=O

k=l

k=2

k=3

tk

Arrival time from interface k

w

8.8

12.94

17.0

18.37

ak

Attenuation coefficient

dB/cm(lMhz)

0

.5

0

.5

k-4

0

of layer k mk

Attenuation factor of interface k

scalar

1.0

9664

.9fA6

ck

Speed of sound at layer k [32]

m/s

1480

1450

1480

‘k

Reflection coefficient of interface k

scalar

-.102

.102

-.102

.102

=k

Acoustic impedance of layer k [321

x106 Rayls

1.48

1.21

1.48

1.21

1.48

Density of layer k [321

g/cm3

1.00

.834

1.00

.834

1.00

.9.55

Transducerandtestobjectparameters The overall reconstructionprocedureis demonstratedby meansof a three layer target medium(testcell) consistingof : mineraloil (3 mm) - water (3 mm) - mineralQil ( 1 mm) ; for this medium,the expectedtheoreticalvalues of the parametersa& tk, $, mk, ckr rk and zk are summarizedin table 1 (alsorefer to figure 1 and1321).The transducerwasa low Q, 1.0 MHz center frequency PZT transducer,manufacturedby K.B. Aerotech. It was unfocused,with a radius of 0.75 cm. This transducerwas placed 4 cm (l/2 of the far field distance)from the multilayer medium. Procedureof reconstructingthe impedanceprofile The final recorded waveform was a normalized sum of 500 single sweeps.The averagingpmcedumincreasedthe SNR from about 30 dB for a singlesweepup to 45-50 dB for the averagedsignal.The digitized averagedwaveform wasfirst storedin the digitizer and thentransfemdthroughaGPIBbustotheIBMPCrok:storedandanolyzad.

56

IMPEDANCE

PROFILE RECONSTRUCTION

A

B

1.81

)r

151

iy

0.6 1.2.. 0.9 m I

!g 0.3 m

Fig. 4 (a) The received trace, y(t), from a three layer medium: mineral oil (3 mm) - water (3 mm) - mineral oil (1 mm). The imaging transducer is a nonfocused 1.0 MHz center frequency. (b) The magnitude frequency response of y(t).

Processingof y(t)

Starting with the measured data y(t), the final goal was to restore the nflectivity function and the impedance profile. figure 4(a) shows the measured waveform y(t), and its magnitude frequency response IY(o)l = IF(y(t))l is shown in figure 4(b). The transducer impulse response h(t) was obtained previously using a wide Plexiglas block as a reflector. Since the reflection coefficient for the water-Plexiglas interface is known, the measured h(t) is scaled to correspond to what would have been received from a reflector with a reflection coefficient of one. The first reflection was separated from the received echo coming from the Plexiglas plate and used as the overall transducer impulse response. figure 5 presents both h(t) and IH(o)l, the magnitude frequency response of this transducer. Next, the Wiener filter W(o) was obtained as in Eq. (6). The first choice for 0 is somewhat arbitrary and can be adjusted later to yield better results. The initial ‘guess’ is an estimated noise to signal ratio (NSR), here chosen to be I$ = 0.06. figure 6a shows the magnitude of the Wiener filter, IW(w)l. As in Eq. (12), the output from the filter, &co) = Y(o>W(o), is shown in figure 6b. The rapid

57

LIFSHITZ ET AL

5i 5 1.8 E 1.5 2

1.2

2

0.9

E Ox5 LG 0.3

0.5

1.0

1.5

2.0

2.5

3.0

FREQUENCY

3.5 (MHz)

(b)

Fig. 5 (a) Impulse response, h(t), of 1.0 MHz transducer. response, IH(f)l, of 1.0 MHz transducer.

(b)

Magnitude frequency

oscillations observed in figure 6b are due to the standing waves phenomena caused by the interference between the different reflections coming from all interfaces. In figure 7a, e(t) = F-l{ 8( w) ), the estimate of the reflectivity function, is presented as given in Eq. (13). The interfaces can now be easily identified. The resolution is far better than that of y(t), although the SNR is degraded.

Reconstruction with the impediography equation Using the impediography equation given in Eq. (2), the impedance profile with no attenuation compensation was obtained and is shown in figure 7b. The reconstruction of f(t) with the impediography equation is appropriate when no information is available as to whether the medium consists of discrete layers or has a continuously varying impedance. Despite the errors, the interfaces are clearly visible in figure 7b. The lack of sufficient bandwidth, even after Wiener filtering, is clearly observable in 2(t) as the profile has the appearance of having been highpass filtered.

58

IMPEDANCE

PROFILE RECONSTRUCTION

Fig. 6 (a) Magnitude frequency response of the Wiener filter, W(f). (b) Magnitude frequency response, I& f)l, of the output through the Wiener filter.

Reconstruction with the Maximum Amplitude Algorithm For the impedance profile reconstruction under the assumption of a planar layered medium, the peak detection algorithm is applied to e(t). In order to emphasize the peaks in k?(t), the function [e(t)]2 is calculated. The array data is now manually divided up into four segments (sub-arrays) : [0, 10.01 , [lO.O, 15.01 , [15.0,20.5] and [20.5, 22.01 (all numbers correspond to time in ps). Each segment contained one of the peaks representing one of the interfaces. The square root of the peak values, the location of the peaks and the polarities corresponding to the peaks for each segment is detected to yield the sets of the arrival times and the estimated reflectivity coefficients. The peak location and polarities, can also be seen directly in figure 7(a). (Notice that once the set ( tk) is obtained, the polarity of each peak is determined from e(t) and assigned to the peak values): (tk)(in ps) = (8.8, 12.28, 20.16, 21.2); (&) = (--.I, .0885, -.0878, .084). The 8.8 t~s value of the first peak is an arbitraq time offset.

59

LIFSHITZ ET AL

0.95

8.0

12.0

16.0 00

20.0 WQ

Fig. 7 (a) The estimation of the reflectivity function, e(t), for a three layer medium. (b) The estimation of the impedance profile using the impediography equation.

Although not easily implemented, the peak detection could have been performed by an automatic peak detection algorithm. Several test criteria in the form of amplitude and peak separation thresholds must be included to preclude detection of erroneous peaks. For example, spurious spikes superimposed on the a true peak must be rejected on the basis of being too close to a true peak, and peaks must exceed a threshold value to be recognized. Reconstruction with the Adaptive Algorithm Alternatively, the *$kcoefficients (of the reflectivity function) can be reconstructe&using the Adaptive Algorithm described earlier. The experimental steps that produce the & coefficients one by one are shown in figure 8. As described earlier, fk(t) is defined as the received signal, y(t) from which the k fist reflections, starting with ye(t), have been stripped. Hence, fo(t) is identical to y(t), and in figure 8(a), fo(t) = y(t) is shown. In figure 8(b), fi(t) is displayed which represents y(t) with the first reflection removed. Figures 8(d), 8(f), and 8(h)

60

IMPEDANCE

PROFILE RECONSTRUCTION

to(t) ,880 ~489 (1) .w -.400 -.wa

9(t) 1.91 .LW I@ ,299 -.ZOD -.608

.l Be

c

*

4,BU

”

’

12.9

2e,u

28.8

’

4,62

1280

4,BO

12ru

29.8

20,e

28.9

28.0

Fig. 8 Application of the Adaptive Method reconstruction algorithm. In (a), (b), (d), (f) and (h); fn(t), fi(t), f2(t) and fs(t), respectively, are shown; in (c), (e), (g) and (i) the estimates & coefficients for k = 0, 1, 2 and 3, respectively, are shown. The final estimate of e(t) is shown in (i).

show fz(t), f3(t), and fq(t), respectively. Since the deletion of this first reflection is based on the estimated magnitude, polarity, and time of the first coefficient k(t). the reflection is not removed perfectly. The magnitude of *&,(t) is determined from energy considerations, as given in Eq. (20), and the polarity and time of arrival are found by the peak detection scheme, applied to the Wiener filtered signal, e(t). h(t) is depicted in figure 8(c). In figures 8(e), 8(g), and 8(i), it is seen how the subsequent coefficients are added, such that 8(i) contains the full set of

61

LIFSHITZ ET AL

reflection coefficients, denoted e(t) where

q(t) = gi$(t-tt) k=O

In the case illustrated in figure 8, M = 3. The coefficients ^ak must of course still be compensated for attenuation effects before the true reflection coefficients are available for the impedance profile reconstruction. Examples of the numerical values for the determination of the ak coefficients are given in the next paragraph. The duration of h(t), T is measured to be 4.681s. The peak detection has located the fiit reflection at to = 8.8~s, with a negative polarity. From Eq. (201, the magnitude of&(t) is found to be 0.1 and therefore h(t) = - 0.1 6(t - 8.8~s). The second reflection is located at tt = 12.28~s, with a positive polarity. Since to + T < tl, the integration to determine it(t) uses to and tl - tu as limits. The magnitude of at(t) is determined to be 0.0943, and hence iI = 0.0943 6(t - 12.28~~). Following the same procedure as demonstrated above, the full set (4) = -.l, .0943, -.0883, .0858, is obtained.

Discussion of the arrival time measw-ements As mentioned under “Transducer and Test Object Parameters”, the three layered medium consists of a 3mm layer of mineral oil, a 3 mm layer of water and a 1 mm layer of mineral oil. The sound velocities for mineral oil and water are 1450 m/s and 1480 m/s, respectively, giving rise to estimated round trip travel times for the three layers of 4.05 ps, 4.14 ps, and 1.38 ps, respectively. When the 8.8 ps offset was removed from the measurement data, the ultrasonically measured round trip travel time are 3.48 ps, 7.88 ps, and 1.04 ps, corresponding to layer thicknesses of 2.52 mm, 5.83 O.mm, and 0.75 mm. The differences between the physically measured and the acoustically measured thicknesses stem most likely from errors in assessing the physical dimensions rather than from the acoustical measurement errors. Because the fluids used in the experiments had densities different from that of water, there was a tendency for the fluids to distort the test cells into a wedge shape due to the difference in hydrostatic pressure from the top to the bottom of the cell. In this example, water is mom dense than mineral oil, causing the test cells to sag near the bottom. The thicker the test cell, the greater the hydrostatic effect and the greater the wedge shaped distortion. The predicted test cell thicknesses were determined from measuring the distance between the membrane solid metal supports before injecting the fluids, not taking into account the wedge shaped distortion.

Experimental data with attenuation compensation The values of the attenuation coefficients is known apriori and given as : a mineraloil = .5 dB/cm @ 1 MHz; Eqs. (23) and (24) are next applied to yield the estimated

62

IMPEDANCE

PROFILE RECONSTRUCTION

reflection coefficients (calculated only for the Maximum Amplitude Method in this example): mo= 1; ?,=-0.1;

% -00878 ?2=;I;;=~=-0.0913 0.966

m,=exp[0]*ml=0.966,

m3= exp[-2%]*

exp[O]* m,=0.989 * 0.966=0.955;

a, i$= -= m3

0084 A=O.O88 0.955

The calculated values for ro, rt, r2. and r3 are -0.102, 0.102, -0.102 and 0.102, respectively as shown in table 1. By applying Eq. (25), the estimatedimpedanceprofile is finally obtained and comparedto the calculatedprofile. The resultsare given in table 2 and figure 9. Figure 9 showsthe calculatedimpedanceprofile (heavy solid line), the reconstructed impedanceprofile usingthe Maximum Amplitude Method (the dashedline) andthe impedance profile obtainedfrom reconstructingwith the Adaptive Method (the solidline). Overall the following mediawere tested:(all layers thicknesseswere measuredbefore injecting the fluids into the testcells) . . .

Mineral Oil ( 1 mm ) - Water - Mineral Oil ( 3 mm ) Mineral Oil ( 3 mm ) - Water(3 mm) - Mineral Oil ( 1 mm) MineralOil( lmm)-Water(3mm)-Peg400( lmm)

Table 3 presentsthe resultsfor the reconstructionof the coefficients:&, ?k and& of the fit medium.For this medium(aswell asthe others)the data were analyzedasdemonstrated aboveusingboth the Maximum Amplitude Method andthe Adaptive Method. Two transducers with center frequencies 1.0 MHz and 2.25 MHz were used. In table 3, the errors in the

Table II Calculatedand reconstructedimpedanceprofile of the medium:mineraloil (3 mm) water (3 mm) - mineraloil (1 mm); ( z valuesarex106 Rayls)

Method

z water

Calculated 1.48 Maximum Amplitude 1.48 Adaptive 1.48

z mineraloil 1.21 1.21 1.21

z water 1.48 1.455 1.473

63

z mineraloil 1.21 1.213 1.226

z water 1.48 1.447 1.468

LIFSHIXZ ET Al

1.51 -) 1 A8 d “2 8 .Q

--------

1.45 1.42 1.39 1.36 1.33 1.30 1.27 1.24 1.21 -

I

I

I

8.0

I

12.0

I

I

16.0

I

I

I

20.0

*

24.0

Fig. 9 Calculated and reconstructed acoustic impedance profiles: Calculated impedance profile (heavy solid line); impedance profile from Maximum Amplitude Method (dashed line); impedance profile from Adaptive Method (solid line).

reconstruction of the reflection coefficients and the impedance profile are defied by 1 T

1 3

Error in r =

c i=l

-z

4

(ii-r;)2

Error in z =

3 c i=l

r:

c i=2

(i

- Zi)2 4

(27)

< c i=2

where ri and Zi are the expected reflection coefficients and impedance values (the coefficients ro, ~0 and ZI are assumed to be known a priori therefore not included in Eq. (27)). IV. CONCLUSIONS In this paper, two algorithms were developed and tested to restore the acoustic impedance profile of a multilayer medium: The Maximum Amplitude Method and the Adaptive Method. These algorithms were presented analytically, tested by computer simulation and experimentally verified under controlled conditions. Both the Adaptive Method and the Maximum Amplitude Method provide the capability to recover the one dimensional acoustical impedance profile of a multilayer medium.

64

IMPEDANCE

PROFILE RECONSTRUCTION

Table III Summary of the reconstruction results of the medium: mineral oil (1 mm) - water (3 mm) - mineral oil (3 mm). All z values are x106 Rayls; 01“ineral oil = 0.5 dB/cm (l.OMHz), 2dB/cm (2.25 MHz). The calculated set of the pressure reflection coefficients and impedances is given by: (rk)= t-0.102, 0.102, -0.102, 0.102) and {zJ= (1.48, 1.21 , 1.48, 1.21, 1.48).

1.0 MHz transducer Max. Amplitude Method

Adaptive Method

ik = -.l, .087, -.086, .081 f, = -.l, .088, -.087, .085 $ = 1.48, 1.21, 1.44, 1.21, 1.44 Error (in r) = 13.4% Error (in z) = 2.3%

& = -.l, .094, -.086, .085 & = -.l, .095, -.087, .089 4 = 1.48, 1.21, 1.47, 1.23, 1.47 Error (in r) = 10.2% Error (in z) = 1.0%

2.25 MHz transducer Max. Amplitude Method

Adaptive Method

ik = -.l, .078, -.069, .052 i, = -.l, .082, -.072, .059 G = 1.48, 1.21, 1.43, 1.24, 1.39 Error (in r) = 30.0% Error (in z) = 4.4%

ik = -.l, .081, -.072, .056 i, = -.l, .085, -.075, .064 & = 1.48, 1.21, 1.44, 1.24, 1.4 Error (in r) = 26.7% Error (in z) = 3.9%

The Maximum Amplitude Method is less complex. Only Wiener filtering and peak detection algorithm is required. Its main disadvantage is low spatial resolution. In other words, for about 75% overlap between two succeeding wall echoes (meaning that the second arrival time tk+l = tk + .25T), this method will fail. The Adaptive Method will yield good results even for overlap larger than 75%. Moreover, this method behaves better under conditions such as low bandwidth and poor SNR. Its main disadvantage is relatively long execution time. In practice, both methods performed with good accuracy. For the different multilayer combinations, the estimated value for the acoustic impedance profile are within 5% of the theoretical values. This small error was obtained despite the larger errors (within 35%) in the recovery of the reflectivity function before attenuation compensation. As predicted by the theory and later confirmed experimentally, even when the impulse response was not reconstructed with a high degree of accuracy, the impedance profile could still be recovered with only small errors. It is assumed that the main source of errors is deviation from normal incidence.

65

LIFSHJTZ ET AL

The choice of the optimal transducer parameters for a given medium is not immediately obvious. Wide bandwidth alone is not sufficient to ensure good reconstruction, and selecting the best center frequency of the transducer is important. &XpiEe narrower bandwidth, imaging with the 1.O MHz center frequency transducer provided better results than those obtained with the 2.25 MHz center frequency transducer. In addition, the physical dimensions of the transducer are of importance. In real situation such as in medical application, the multilayer medium will not in general contain perfectly parallel layers. In this case, transducer angle misalignment relative to any one interface can create significant error. Here, the 1 MHz, 0.75 cm radius transducer produced little distortion up to 3 degrees angle of misalignment. Oti the other hand, with the same amount of angle misalignment, using a 2.25 MHz, 0.5 cm radius transducer, no reasonable recovery could be obtained. Attenuation loss compensation is necessary in order to achieve reliable reconstruction. In this work both the attenuation profile and the velocity profile were assumed to be known a priori. It is desirable to further develop and implement inverse techniques to be able to restore the entire set of acoustical parameters (impedance, attenuation and velocity) simultaneously for such a medium. ACKNOWLEDGMENT The assistance of people at Cardiometrics Inc. is greatly appreciated. This work was partly supported by the NSF grant no ECE-8504602.

REFERENCES

Ill

Robinson, E.A., Durrani, T.S., and Peardon, L.D., Geophysical Signal Processing, Chap. 6, pp. 120-211 (Prentice-Hall Inc., London, UK, 1986).

VI

Kak, A.C., Fry, F.J., and Jones, J.P., Acoustic Impedance Profiling, in Ultrasound: Its Applications in Medicine and Biology, Part 1, F.J. Fry, Ed., pp. 511-516 (Elsevier, New York, 1978).

[31

Jones, J.P., and Leeman, S.L., Ultrasonic Tissue Characteriation and Quantitative Ultrasound Scatter Imaging: Method and Approaches, in Proc. Znt. Workshop on Physics and Eng. in Medical Imaging, pp. 247-258, IEEE Publication CH1751-7, (IEEE, New York, 1982).

141

Lee, Q.. Wave propagation and profile inversion in lossy inhomogeneous media, Proc. IEEE 74219-228 (1982)

151

Leeman, S., Ferari, L., Jones, J.P., and Fink, M., Perspectives on attenuation estimation from pulse-echo signals, IEEE Transaction on Sonic and Ultrasonics SU31, 352-360 (1984).

[61

Hayward, G., and Lewis, J., Computer Aided Design of Inverse Filter Algorithms for Ultrasonic System, in Ultrasonic International Conference Proceedings, pp. 716-721 (London, UK, 1987).

66

IMPEDANCE

171

@I

PROFILE RECONSIRUCTION

Fatemi, M., and Kak, A.C., Ultrasound B-scan imaging: theory of image formation and technique for restoration, Ultrasonic imaging 2,1-47 (1980). Kino, G.S., Acoustic Waves Devices, Imaging & Analog Signal Processing, Chap. 4,

pp. 488-498 (Prentice-Hall Inc., Englewood, NJ, 1987). 191

Levin, M.J., Optimum estimation of impulse response in the presence of noise, IRE Trans. on Circuit Theory, 7,50&O (1960).

[101

Rosenfeld, A., and Kak, A.C., Digital Picture Processing, (Academic Press, San Diego, CA, 1982).

1, Chap. 7, pp. 268-347

1111 Widrow, B., Adaptive Filters, in Aspects of Network and System Theory, pp. 563587 (Halt, Rinehart and Winston, New York, 1970). WI

Behar, D., Kino, G.S., Bowers, J.E., and Olaisen, H., Storage correlator as an adaptive inverse filter, Electron. Left. 16, 130-131 (1980).

u31

Mahalanabis, A.K., Prasd, 8, and Mohandas, K.P., A fast optimal deconvolution algorithm for real seismic data using Kalman predictor model, IEEE Trans. Geoscience Remote Sensing GE-19,216-221 (1981).

[I41

Mendel, J.M., and Kormylo, J., New fast optimal white noise estimators for deconvolution, IEEE Trans. Geoscience Remote Sensing GE-15,32-41 (1977).

WI

Crump, N.D., A Kalman filter approach to the deconvolution of sesmic signals, Geophysics 39, 1-13 (1974).

P61 Demoment, G., and Reynaud, R., Fast minimum variance deconvolution, IEEE Trans. on Acoustic, Speech, and Signal Processing ASSP-33, 1324-1326 (1985). H71

Herment, A., Demoment, G., and Vaysse, M., Algorithm for On Line Deconvolution of Echographic Signals, in Acoustical Imaging, Vol. 10, P.Alais and A.F Metherell, Ed., pp. 325-344 (Plenum Press, New York, 1982).

WI

Demoment, G., Reynaud, R., and Herment. A., Range Resolution Improvement by a Fast Deconvolution Method, in Acoustical Imaging. Vol. 6, N. Booth, Ed., pp. 435451 (Plenum Press, New York, 1984).

[I91

Oppenheim, A.V., and Shafer, R.D., Digital Signal Processing, pp. 480-532 (Prentice-Hall Inc., Englewood, NJ, 1975).

PO1

Schmolke, J., Hiller, D., Ermert, H., Schaefer, J.O., and Grabner, G., Generation of Optimal Input Signals for Ultrasound Pulse-echo System, in IEEE Ultrasonic Symposium Proceedings, pp. 929-934, IEEE Publication CH 1832-4 SU, (IEEE, New York, 1982).

67

Chap. 10,

LIFSHITZ ET AL

PI

Pedersen, P.C., Tretiak, O.J., and He, P., Numerical techniques for the inverse acoustical scattering problem in layered media, Acoustical Imaging 12, 443-457 (1982).

1221 Pedersen, P.C., Tretiak, O.J., and He, P., Transmission and Reflection Properties of Inhomogeneous Layers and their Application to Medical Ultrasound, in Proc. 4th Annual conf. IEEE Eng. Med. Biol., pp. 551-556(IEEE, New York, 1982).

WI

He, P., Direct and inverse acoustical scattering problem for a medium with continuousimpedanceprofile, M.S Thesis,Drexel University (1981).

[241

Bai, J., One dimensional acoustic scattering: The direct and inverse problem underrealistic conditions,PhD. Thesis,Drexel University( 1985).

WI

Pedersen,PC., Tretiak, O.J., and Bai, J., The Inverse Acoustical ScatteringProblem for a Layered Media in the Presenceof BroadbandAcoustic Noise and with Limited TransducerBandwidth, in Acoustic Imaging, Vol. 13, M. Kavek, R.K. Mueller, and J.F. Greenleaf, Ed., pp. 57-73 (PlenumPress,New York, 1983).

VI

Jones,J.P., Impediography: A New Ultrasonic Techniquefor Diagnostic Medicine, in Proc 1974 Meeting of the AIUM, Vol. 1, pp. 499-508 (ALUM, Bethesda,MD, 1975).

WI

Ware, J.A., and Aki, K., Continuousand discrete inverse scatteringproblemsin a stratified elastic medium,J. Acoust. Sot. Amer. 45,911-921 (1969).

WI

Pedersen,P.C., Tretiak, O.J., and He, P., Impedancematching properties of an inhomogeneousmatching layer with continuously changing acoustic impedance, J. Acoust. Sot. Amer. 72, 327-336 (1982).

WI

Schwetlick, H., andKessel,W., Deconvolution of bandlimited signalsusinga-priori knowledge of their impulseresponse,Arch. Elektron. and Uebertragungstech41, 62-64 (Germany, 1987).

[301

Papoulis, A., and Chamzas, C., Improvement of range resolution by spectral extrapolation, UlrrasonicImaging 1,121-135 (1979).

[311

Lewin, P.A., and Schafer, M., Wide-band piezoelectric polymer acoustic sources, IEEE Trans.on UitrasonFerroelec Freq Control 35,175183 (1988).

~321

Selfridge, A.A., Approximate material properties in isotropic materials, IEEE Trans. SonicsUltrasonicsW-32,381-394 (1985).

68