Beamformed nearfield imaging of a simulated coronary artery containing a stenosis

This paper is concerned with the potential for the detection and location of an artery containing a partial blockage by exploiting the space-time properties of the shear wave field in the surrounding elastic soft tissue. As a demonstration of feasibility, an array of piezoelectric film vibration sensors is placed on the free surface of a urethane mold that contains a surgical tube. Inside the surgical tube is a nylon constriction that inhibits the water flowing through the tube. A turbulent field develops in and downstream from the blockage, creating a randomly fluctuating pressure on the inner wall of the tube. This force produces shear and compressional wave energy in the urethane. After the array is used to sample the dominant shear wave space-time energy field at low frequencies, a nearfield (i.e., focused) beamforming process then images the energy distribution in the three-dimensional solid. Experiments and numerical simulations are included to demonstrate the potential of this noninvasive procedure for the early identification of vascular blockages-the typical precursor of serious arterial disease in the human heart.


I. INTRODUCTION
C ORONARY artery disease (CAD) is a primary precursor of heart attacks-the leading cause of death in the United States. CAD is the buildup of plaque in the coronary arteries resulting in a condition referred to as stenosis, in which blood flow is restricted and the oxygen supply to the heart muscle is decreased [1]. Presymptomatic detection of coronary stenosis is extremely important because, at this stage, the disease can be treated by diet and medication rather than by such costly and invasive procedures as angioplasty or coronary bypass surgery. On a wide scale, early CAD diagnosis would not only preserve the quality of life for millions of patients by allowing the medical community to intercede before serious heart damage occurs, but it would also save the United States health care system billions of dollars through the avoidance of extensive medical intervention.
The most common technique for the identification of constrictions in coronary arteries is a procedure called an angiogram. During this process, a catheter injects a radiopaque Manuscript received March 2, 1998;revised October 26, 1998. This work was supported by the Office of Naval Research. The Associate Editor responsible for coordinating the review of this paper and recommending its publication was M. Insana. Asterisk indicates corresponding author.
contrast agent into the coronary arteries, while an X-ray screen simultaneously images the arterial system. However, an angiogram is costly and its invasiveness places the patient at risk.
Attempts to noninvasively diagnose CAD [2] and intercranial vascular lesions and aneurysms [3] by passive listening (auscultation) for the effects of blood flow turbulence, which is frequently distal to a stenosis [4], have been reported in the recent literature. The potential for these techniques to result in low-cost, mass-screening examinations is considerable. However, the feasibility of such noninvasive detection for other than a near-surface artery depends on the capability to discriminate signals with a very low-level coronary-lesioninduced blood flow turbulence simultaneously in the spatial and temporal domains from any masking noise source. Examples of such noise interference would include blood flow turbulence caused by diastolic interval valve flow and coronary artery branch point flow velocity gradients.
The space-time signals of interest result from turbulentblood-induced forces normal and parallel to a diseased artery internal wall, which, in turn, couple and propagate as wave energy to the body surface as a time-varying spatial vibration field. The research presented here demonstrates the principle that this space-time field can be used to image the energy distribution within a viscoelastic medium by placing an array of noninvasive vibroacoustic sensors on a model of the human body. The spatial discrimination properties of the sensor array are then used as a generalized focusing antenna to enhance the spatial localization and signal-to-background interfering noise ratio relative to that which could be achieved with a single sensor [5]- [7]. It shall also be demonstrated here that shear wave energy in the frequency decade from 100-1000 Hz [8] is an exploitable energy transport mechanism.
A murmur in a coronary artery induced by blood flow turbulence can be discriminated acoustically from masking background sound in the temporal, spectral, and spatial domains. For example, temporal discrimination against background noise in listening for CAD is based on the simple fact that sounds characterizing occlusive coronary artery stenosis occur during diastole [9]. Energy from stenosed coronary arteries must therefore be extracted from a background noise environment that occurs during the filling of the ventricle chambers. Additionally, it has been suggested that resonance effects in arteries excited by blood flow produce characterizable, albeit low-level, sonic frequency band signatures [10]. Double resonant frequencies below 400 Hz are predicted, whereas measured CAD-related resonant energy as high as U.S. Government work not protected by U.S. copyright 1100 Hz due to turbulent blood flow caused by lesions on the inner wall of a partially occluded artery is reported [2]. In a related study of the stenotic carotid, quantitative blood flow rate may be measurable in the absence of highly turbulent flow-induced resonance by estimating the level and cutoff frequencies of the broadband energy produced by low-level turbulence, which is, in turn, related to the diameter of the stenotic carotid segment [11]. In the study presented here, the key to reliable detection, identification, and location of arterial lesions and aneurysms may very well lie in the spatial domain. It is thus proposed that the detection and location of vibrationinduced acoustic frequency band energy from blood flow can potentially be accomplished through the simultaneous filtering and combining of the outputs of a spatially diverse array of vibroacoustic sensors.
In diagnosing CAD, it must be recognized that the signal levels produced by the effects of relatively small coronary artery blood flow rates are substantially below the levels produced by nearby valve flow activity, such as ventricle refilling during diastole. Moreover, the ambient environmental noise produced by noncardiac functions is a recognized inhibitor to low-level sound detection, localization, and diagnosis. It is also well known that the audible sonic field produced by the cardiac cycle exhibits a high degree of spatial variability [12]. The traditional auscultation in phonocardiography involves the placement of a single vibroacoustic sensor (namely, a stethoscope) in the area associated with high-quality energy transmission from a particular region of the heart. In contrast, multiple auscultation point (MAP) space-time phonocardiography proposes to receive energy simultaneously at an array of sensors. Through the application of space-time signal processing technology, the MAP approach could provide the following.
• An improved signal-to-noise ratio (SNR) estimate of the waveform produced by the turbulent blood flow activity in a specific region of the body. This is accomplished by coherently combining the spatially diverse samples of the vibration field to enhance the desired signal and discriminate against interfering noise in the measured space-time field. • A spatially and spectrally resolved image of the turbulentflow-induced, spatially extended vibratory energy source for either a resonant phenomenon or a spectrally wideband phenomenon, such as produced by valve clicks and valvular flow. Both enhanced SNR waveform estimation, which is critical to the development of low-error-rate detection algorithms for low-level signal activity, and the high-resolution MAP imaging of energy fields due to turbulent fluid jets would be invaluable for the prognosis and treatment of lesions and aneurysms. A technique combining these two capabilities could potentially provide a relatively low-cost alternative to invasive and expensive imaging procedures.
The new method proposed to detect and locate arterial stenoses uses phased array sensor technology to image the shear wave space-time energy field produced in the volume under the array and in the acoustic frequency band [13]. This approach is based on the turbulent normal wall forces in an artery that are present in and downstream from a stenosis [4]. These wall forces, while extremely small, exert a pressure on the wall of the vessel, which causes a displacement in the surrounding soft tissue. This displacement occurs in the form of two waves in the medium: a compressional wave and a shear wave [14]. The shear wave energy that dominates the response at low frequencies (100-500 Hz) is measured at the surface of the solid by a spatially diverse array of displacement sensors, and a MAP nearfield (focused) beamforming imaging processor allows detection and location of the turbulent source. An in vitro experiment is described to validate the physical principles that enable this method, and numerical simulations are presented to confirm these experimental model results.
In the next section of this paper, the principles of a nearfield focused beamformer are presented and related to the volumetric imaging process. This is followed in Section III by a description of the numerical simulation of the focused beamformer imaging process that includes a discussion of the effects of the spatial extent of shear wave energy sources due to turbulent flow in a tube. Section IV details the in vitro urethane model experimental setup, and Section V presents the comparison of numerically predicted and experimentally measured results for urethane models with and without embedded ribs. In Section VI, both the conclusions and intended future efforts are discussed.

II. FOCUSED BEAMFORMING
The MAP beamformer is defined by first assuming that is a zero-mean, stationary, complex vector of sensor random outputs that are Fourier transformed at frequency Next, the ensemble-averaged image scalar pixel value of the so-called conventional focused beamformer (CFB) at Cartesian coordinates and is expressed by [15] (1) In (1), is the cross-spectral density matrix (CSDM), is the statistical expectation operator, the superscript denotes the complex conjugate transpose matrix operation, and is the -by-1 focused MAP imaging column vector, with the nth element given by (2) In (2), is the geometric spreading loss factor, is the average wave speed, and is the Euclidean distance from the th image point to the th sensor located at and given by (3) All CSDM's used throughout this paper to produce beamformed images have been normalized, as discussed in the Appendix.
By scanning the image point throughout a volume, the MAP beamformer produces an image of an underlying volume based on the energy intensity sensed at each focal point. The vectors are typically called the steering vectors. The numerator of (1) provides the two-or three-dimensional (3-D) image, which is unscaled with respect to geometric spreading, and the denominator normalizes the image so that it is correctly scaled at each image point location. When this process is applied to the data on a time-averaged basis, the high-energy spatial source locations produce proportionately larger image scalar values.

III. SIMULATION
Numerical simulation of this problem consists of analytically formulating the CSDM, inserting it into (1), and then calculating the simulated image that would be produced by the focused beamformer. The analytical model of the turbulent flow energy field, based on previous laboratory wall pressure measurements downstream from a blockage [4], is modified to include energy in the region of the blockage. Additionally, a spatial correlation term is added to the model to account for the spatial coherence of the turbulent flow [16]. The expression for the streamwise wall pressure autospectrum is given as (4) where is the amplitude of the wall pressure shown in Fig. 1, is the streamwise location of the flow measured from a reference point (m), is the spatial decay constant (dimensionless), and is the fluid convection speed (m/s). In Fig. 1, the ordinate is the root mean square (rms) wall pressure nondimensionalized by the fluid density , the square of the constriction velocity , the constricted diameter , and the unconstricted diameter ; the abscissa is the spatial location For this study, the spatial decay constant was 0.125 and the  fluid convection speed was the measured flow velocity through the constricted portion of the tube. Fig. 2 is a plot of the absolute value of (4), without the amplitude term , versus spatial location at a frequency of 300 Hz. Fig. 3 shows a plot of the absolute value of (4), without the amplitude term versus frequency when mm (0.125 in). Figs. 2 and 3 illustrate the spatial decay of the wall pressure amplitude term versus spatial distance and frequency, respectively.
Equation (4) is used to calculate the CSDM component produced by the source (blockage and downstream energy), which has the element in the th row and th column expressed as shown in (5) at the bottom of the page.
For the distributed multiple source model used here, the indexes and range from 1-21 and correspond to point sources at the circled data points in Fig. 1, and is the maximum value of the pressure with respect to spatial Once this matrix is known, the simulated CSDM can be computed from (6) where is an matrix giving the response of the th sensor to the th unit amplitude level source as represented by the th element (7) In (7), is the Euclidean distance from the th sensor to the th source, given as (8) From (6), is the autospectral density (ASD) of the assumed zero-mean, spatially uncorrelated additive noise at each sensor and is the average signal ASD at the peak of the spatial spectrum at a sensor output.

IV. EXPERIMENT
A homogeneous urethane mold containing a latex surgical tube was manufactured for the in vitro tests, as shown in Fig. 4. The urethane, which is Hexcel 195-RE, represents the soft tissue surrounding the simulated artery and has a shear wave  speed of 7 m/s at 100 Hz and 12 m/s at 1000 Hz. A Metravib Viscoanalyzeur was used to measure the shear modulus of the urethane from which this shear wave speed was calculated. The 100-1000-Hz shear wave speed is approximated by using a linear dispersion relationship with respect to frequency. It is noted here that the slowest shear wave speed available in a commercial grade urethane is somewhat faster than soft human tissue shear wave speed would be at low frequencies [8]. The surgical tube, intended to simulate the mid-left anterior descending coronary artery, has an inner diameter of 3.18 mm (0.125 in) and is embedded at a depth of 31.8 mm (1.25 in) in the urethane. It contains a flow constriction (nylon), which represents a coronary (partial) occlusion having a length of 12.7 mm (0.5 in) and a hole with an inner diameter of 1.59 mm (0.0625 in), as shown in Fig. 5. This blockage was chosen because it has a 50% diameter reduction (or a 75% area reduction), which is similar to the reduction used by other researchers who have measured break frequencies [17], velocity [18]- [20], and pressure [4], [21], [22] caused by a constriction in a cylindrical tube. Moreover, this degree of blockage would be typical of presymptomatic CAD. The inflow end of the nylon occlusion contains a 45 taper to provide a relatively smooth entrance region.
Water, with its kinematic viscosity of 1.01 10 m /s, is gravity fed into the surgical tube from a holding tank and captured below the urethane mold in a calibrated flask, as shown in Fig. 6. Although this is a steady-state flow experiment, it has been shown that the turbulent wall force levels are nearly identical to a corresponding pulsatile flow experiment [22]. The laboratory configuration permits measurement of the flow volume from which corresponding velocities and Reynolds numbers, based on the obstructed and unobstructed diameters, can be computed. The moving fluid turbulence transmits a normal force into the embedded tube wall, which, in turn, transmits energy into the urethane. This energy transfer is due mainly to the compliance of the tube wall immediately downstream of the stenosis. Some energy is also transferred from the region directly inside the stenosis, although this is at a lower level because the nylon constriction is stiffer than the surrounding surgical tube. The energy then propagates to the surface of the urethane, primarily in the form of shear waves, where it produces a local displacement field that fluctuates in space and time. After a polyvinylidene fluoride (PVDF) piezoelectric sensor array measures this field, the MAP beamformer images the energy distribution in the 3-D solid. The experiment was conducted in an electromagnetically (EM) shielded anechoic chamber to minimize ambient noise contamination.
All the measurements in this paper were taken with a twodimensional (2-D) array that consists of 15 sensors (elements) on two orthogonal axes, as shown in Fig. 7. A common middle sensor is shared by the nine sensors (equally spaced) on the primary axis and the seven sensors (not equally spaced) on the secondary axis. The element-to-element midpoint spacing on all of the primary axis sensors is 9.53 mm (0.375 in). The element-to-element midpoint spacing of the outer three sensors on each secondary axis is also 9.53 mm (0.375 in), and the distance from the common middle sensor to its adjacent sensor is 28.6 mm (1.125 in). Each group of outboard sensors on the secondary orthogonal axis is built into a separate aluminum box that is attached to the main box by hinges so that the array can conform to nonplanar surfaces. A coordinate system is chosen such that the primary axis is the -axis and the secondary axis is the -axis. The -axis corresponds to the depth of the image point into the urethane block. The coordinate system origin is the intersection of the primary and the secondary array axes. The main axis of this array was designed to fit in an intercostal region between the ribs of a human, and the secondary axis was designed to fit in the adjacent upper and lower intercostal regions.
The sensing elements of the array are PVDF film strips that are 38.1 mm (1.50 in) long and approximately 9.27 mm (0.365 in) wide. Attached to the middle of each strip is a nearly circular piece of plastic, which has a 9.27-mm (0.365in) approximate diameter and a 1.59-mm (0.0625-in) thickness. This circular plastic contact piece is intended to provide an omnidirectional response to propagating shear waves with wavelengths that are greater than the spacing between the adjacent surface contacts. From the spatial sampling perspective, the spacing between contact centers (9.53 mm) is appropriate for assumed wavelengths of 37.8 mm at 200 Hz and 18.4 mm at 500 Hz. The sensitivity of the PVDF film along its major axis is 22 176 V/unit strain for a standard thickness of 28 10 mm. The sensor outputs of the array were connected to a multichannel low-noise preamplifier with lightweight multistrand coaxial cables (Cooner Corporation). The received signals were passed from this preamplifier to a bank of amplifiers (Ithaco Corporation Model 451) with a 6 dB/octave rolloff below 1000 Hz to equalize the spectrum level and a multichannel programmable filter (Precision Filter Corporation) with a 48 dB/octave rolloff above 800 Hz to prevent spectral aliasing in the digitizer sample-and-hold operation. Finally, the signals were passed to a data acquisition system (LabView) residing in a computer (Apple Power Macintosh 8100/110), where the time series data were stored for postprocessing.
Two typical beam patterns of the array are shown in Figs. 8 and 9. They are computed using the steering vectors according to the relationship (9) where (10) and the subscript denotes the focal point of the array. Fig. 8 is computed at 200 Hz and Fig. 9 is computed at 500 Hz.

V. EXPERIMENTAL VALIDATION OF THE NUMERICAL MODEL
Three different geometrical configurations were numerically modeled, and corresponding experimental measurements were obtained to validate the numerical model predictions. The first  is an ideal case, where the main axis of the array is parallel to the surgical tube and the middle sensor is directly over the occlusion, as shown in Fig. 7. Figs. 10 and 11 are comparisons of the simulation and the experimental data for 300 and 500 Hz with the CFB. The surfaces are normalized such that the maximum value is 0 dB and all values in the beamformer response surface below 10 dB are mapped into 10 dB. The experimental data were computed with 48 fast Fourier transforms that were each 1024 data points long; the sampling frequency was 2016 Hz per channel, which corresponds to a frequency bin resolution of 1.969 Hz. The wave speed for computing all the images was linearly interpolated for each specific frequency, resulting in a wave speed of 8.11 m/s at 300 Hz and 9.22 m/s at 500 Hz. The spreading parameter was 1, which represents spherically symmetric propagation spreading attenuation at each source location. The Reynolds number in the surgical tube (Re was 1790, and the fluid velocity was 0.57 m/s. The Reynolds number through the constriction (Re was 3580, and the fluid velocity was 2.28 m/s. Although it should be noted that the effective constricted diameter (the vena contracta) is slightly smaller than the actual constricted diameter [23], numerical simulations with smaller constricted diameters produced the same results as those using the actual constricted diameter.
In Figs. 10 and 11, the -axis represents the coordinate that is parallel to the main array axis, with the middle sensor centered at , and the -axis is image depth into the urethane solid. The two parallel black lines centered at mm represent the surgical tube in the urethane. The solid black rectangle centered at and mm is the nylon constriction. The fluid is flowing from the right side of the image to the left side. The yellow rectangles at depict the individual sensor locations on the main axis of the array; the sensors on the secondary axis of the array are not shown. The color bar on the right-hand side of each image is the scale of the image in decibels. The second geometrical configuration illustrates the case where the main axis of the array is orthogonal to the surgical tube and the blockage is between two of the sensors of the array. Shown in Fig. 12, this would be representative of a human subject where an artery with a partial occlusion is at some unknown location. Figs. 13 and 14 are comparisons of the simulation and the experimental data for 200 and 400 Hz with the CFB. The Reynolds number in the surgical tube (Re was 1820, and the fluid velocity was 0.58 m/s. The Reynolds number through the constriction (Re was 3640, and the fluid velocity was 2.32 m/s. The simulated and experimental constants used in the first case were also used in this one.
In Figs. 13 and 14, the -axis represents the coordinate that is parallel to the main array axis, with the middle sensor centered at , and the -axis indicates the image depth into the urethane solid. The arrow-marked black circle, centered at mm and mm, represents the surgical tube in the urethane. The fluid is flowing up and out of the image and is orthogonal to the image plane. The yellow rectangles at depict the sensor locations on the main axis of the array; the sensors on the secondary axis of the array are not shown. The color bar on the right-hand side of the images is the scale of the image in decibels.
Modeling energy propagation in a homogeneous solid is a straightforward task. Energy travels in paths that can be predicted with simple straight-line propagation. Unfortunately, the human thorax is not homogeneous. A primary source of inhomogeneity in the medium between the left anterior descending coronary artery and the surface of the chest is the rib cage. Forward-propagating energy reaching the ribs tends to be both scattered and spectrally reflected from the ribs, altering the space-time field that is measured on the surface of a human when compared to the surface of a homogeneous solid. Additionally, the ribs reradiate energy by converting the wavenumber of incoming waves to create an extended space-time field.
An initial understanding of how shear wave energy propagates through the ribs required the third test, which was identical to the second experiment, with two exceptions. First, a set of dog ribs was embedded in the urethane between the surgical tube and the top surface of the mold, as shown in Fig. 15. The spacing (intercostal region) between each rib was approximately 9.53 mm (0.375 in). Second, the array was placed in an intercostal region of the dog ribs, requiring that it be located slightly downstream of the occlusion [approximately 6.4 mm (0.25 in)] rather than directly over it, as in the previous experiment. The 15-channel planar array was used to acquire the data, which was processed with the nearfield beamforming algorithms given by (1)-(3). No  At some frequencies, the ribs tend to make the source detection more accurate, as shown in Figs. 16 and 17. However, the images were not as consistent across frequency bins for the case with ribs as were the images from experiments one and two for the case without ribs. (Because of space considerations, all of the frequency images are not shown.) There were several frequencies where there was a second virtual source that imaged in a location that appeared to be a forward path reflection artifact. It is believed that this was a coherent reflection of the actual energy source from the ribs, although the exact cause has yet to be determined.
The conventional focused MAP beamformer produced images that variously detected and localized the blockage. As the frequency increases, the images have better resolution, which is the result of the shorter wavelength energy producing a beamformer with a more focused (smaller spot size) main  beam. However, as the distance to the focal point (into the imaged volume) increases, beams produced with this type of MAP imaging algorithm provide less resolution. Radial resolution decreases in proportion to the squared image focal point depth. The source depth of 31.8 mm (1.25 in) was chosen because this is the shortest distance from the chest wall to the midpoint of the left anterior descending coronary artery in the average adult human [24].
It is noted that the simulations and experimental data for the MAP beamformer are in good, although not exact, agreement. There are several possible reasons for the variance.
• The turbulent energy is modeled with a line-segment source that has no physical extent in the radial direction, whereas the experiment uses a tube that is 1/8 in in diameter. • The point-to-point correlation function of the source model was derived from a flat plate model, whereas the point-to-point correlation of the experimental tube may be slightly different. • The line-segment source model uses an rms pressure frequency-independent measurement rather than a frequency-dependent measurement that could have a different spatial shape at each frequency. • An ideal spherical propagation model was assumed, even at distances from the array that are short compared with a wavelength. In addition, the blockage used in the previous source measurements [4] had a blunt inlet region compared with the tapered inlet region used in this experiment. Also, the background noise assumed in the numerical simulations is spatially uncorrelated with zero mean and uniform variance from sensor to sensor, whereas the noise in the actual experiment could have angular dependence. This type of noise mismatch in the model would cause a bias error for energy location because the focusing process has inaccurately described the effective noise directional wavenumber spectrum.

VI. CONCLUSIONS
In summary, it has been demonstrated that the turbulent flow energy in and downstream from an occlusion creates a normal wall pressure on the tube that is converted to shear waves which, in turn, are dominant over compressional waves in an elastic solid at frequencies of 100-500 Hz. These shear waves propagate to the surface, where they can be measured by a vibration displacement sensor array. Signal processing is then accomplished by a system that includes a MAP nearfield beamforming algorithm to perform imaging of the spatial shear wave energy field intensity. Such imaging of the volumetric region under the array by focal point scanning produces a representation of the internal energy field. The use of a 3-dB contour image deflection criterion in the study allows the occlusion to be detected and localized to within 10 mm of depth and crossdepth at frequencies near 500 Hz.
Further research should investigate the capability of this method to more accurately locate the blockage by using a larger aperture array, by reducing the sources of focusing MAP beamformer mismatch, by using high spatial resolution beamformer signal processing, and by more fully characterizing the temporal spectrum of energy that is produced by the turbulent flow excitation. Additionally, complex molds that are more representative of the in vivo conditions that the experiment must replicate should be investigated so that the ultimate goal of this effort-the detection of early stage arterial disease-will be achieved.

APPENDIX CSDM NORMALIZATION
Before the MAP beamformer operates on the CSDM, the matrix is normalized to offset the effects of channel-to-channel amplitude discrepancies in the data. Although the array used to acquire the data can be calibrated in the laboratory, previous experience has shown that the calibrations change moderately when the array is placed in an operating environment. The normalization helps to moderate the effect of the changing calibration. It is based on a vector of the diagonal elements of the CSDM, (where is obtained by time averaging) and is expressed as . . . (11) where the bar denotes the unnormalized CSDM. With this vector, a matrix can be computed by using the outer product of the vector with itself as (12) Once is known, the entries of the CSDM can be individually normalized with the equation (13) where the subscripts and denote the individual entries in each row and column, respectively. This normalization was applied to all of the experimental data and numerical simulations presented in this paper.