2.1.1.

Model

The temperature-change estimation method described in this paper is based on the thermal dependence of the ultrasound echo that accounts for two different physical phenomena: (1) local change in speed of sound due to changes in temperature and (2) thermal expansion of the propagating medium. The former produces an apparent shift in scatterer location and the latter leads to a physical shift. Along an A-line, however, the two effects lead to echo time-shifts that can be estimated and are shown to be related to local change in temperature in the medium. In many situations of pratical interest, these effects are typically small, so that a linearized approach can be used in the analysis that follows. It should be mentioned, however, that for large temperature rise and/or large size of the heated region, the linearized analysis breaks down. Different formulation is needed to handle these cases. This formulation will be introduced in a future report.

The baseline temperature in the sample (tissue) medium is assumed to be constant and equal to θ₀ at the initial time T₀. We start by considering the time delay of the echo from a scatterer at axial depth z as being exclusively a function of the temperature-dependent speed of sound in the propagating medium:

where θ(ξ) = θ₀ + δθ(ξ) is the temperature at depth ξ, c(ξ,θ(ξ)) represents the speed of sound at depth ξ and temperature θ(ξ), and the subscript c indicates that only the thermal dependence of speed of sound has been considered. In order to account for thermal expansion in the propagating medium we substitute dξ by (1 + α(ξ)δθ(ξ))dξ in Eq. (1), where α(ξ) is the linear coefficient of thermal expansion of the medium at axial depth ξ. We get:

At the initial time T₀ = 0, the echo time delay is . Once the temperature in the medium changes, a time-shift will be observed on the echo from scatterer at axial depth z:

Differentiating Eq. (3) with respect to the scatterer location z, we get:

The thermal dependence of the speed of sound in tissue can be assumed to be approximately linear in the temperature range of interest⁹:

where , and c₀(z) = c(z, θ₀). From Eqs. (4), (5), and using the fact that |β(z)δθ(z)| ≪ 1, solving for δθ(z) = (θ(z) – θ₀) we get:

The term is a medium (material) dependent parameter. For a homogeneous medium this parameter can be experimentally determined. For the remainder of this paper we assume k(z) and c₀(z) to be invartiant with respect to axial depth z, i.e., we assume k(z) = k and c₀(z) = c₀. In practice, these parameters should be allowed to vary in space and as a function with temperature. For the purposes of this paper, however, we will limit our attention to homogeneous media and small changes in temperature so that these assumptions are approximately valid. It should be mentioned that the linearized approach is valid for temperature changes on the order of 10 - 15 °, even in inhomogeneous media. Therefore, the temperature estimation algorithm described below will be of practical for temperature estimation in hyperthermia and, more importantly, for image guidance of thermal therapy devices.

2.3.1.

Image Guidance for Noninvasive Surgery

Figure 1 shows an experimental setup used in image guidance based on the temperature imaging method described in this paper. In this experiment, a thick piece of bovine muscle was subject to a sequence of heating patterns due to the therapeutic array. Eight 2-second duration foci were used for simulating a coagulative thermal necrosis in tissue. The individual foci were spaced by two mm and deivered to the target in sequence with the foci approaching the imaging transducer. Ample wait time between shots was given to allow temperature to return to baseline before the next shot is delivered. Upon completion of the experiment, the tissue sample was cut and it was confirmed that tissue destruction (as indicated by visible change of color) at the desired location. In order to demonstrate the capability of the temperature estimation algorithm for guiding the therapy, the sequence of shots were initially delivered at low power (causing only 2 °C rise at the focal point. That is, the therapy array was operated well below the threshold for tissue damage and was run in a sequence that resembles the planned therapy session. While this was being done, temperature images were collected and processed to determine if the locations and quality of the foci are consistent with the treatment objectives. This allows for characterizing the quality of the heating foci in-situ at subthreshold (nondestructive) levels just before the theraputic pulses are to be delivered. Figure 2 shows the spatio-temporal temperature distribution along one image line going through all eight heating spots. The grayscale temperature image is obtained from an M-mode image consisting of 128 A-lines. For each shot, 16 A-lines were acquired and processed for temperature estimation, 2 before, 4 during, 5 after (fast rate), and 5 after (slow rate). Please note that the time axis is simply labeled A-line Number. For the temperature estimation experiment shown, each 16 lines represent 1 minute, but the A-lines were not acquired uniformly in time. Figure 2 shows clearly that the temperature estimation algorithm has the spatial contrast and sensitivity needed to characterize the quality of the therapeutic transducer. The peaks of the temperature rise profiles are spaced by 2 mm as planned. Also, the extent of each heated region is between 3-5 mm axially which also corresponds to the expected temperature profile of a 2 mm wide beam heating for two seconds. While the image clearly demonstrate the potential of the temperature estimation algorithm, it should be noted that some level of ripple artifact exists. However, this ripple is not a limiting factor in this case. It should also be emphasized that the maximum temperature rise needed to produce this result is approximately 2 °C. This is a major advantage for ultrasound in that the quality assurance of focusing and energy deposition can be performed energy levels well below the therapeutic thresholds. Furthermore, since these images can be acquired and processed within a fraction of a second, the initial heating rate can be evaluated two dimensionally. Therefore, it is possible to visualize the specific absorption rate (SAR) of the energy source in situ. These are perhaps the major strengths of using ultrasound in image guidance of energy sources for noninvasive or minimally invasive thermal surgery.

Figure 1.

Experimental setup used in the image guidance experiment.

Figure 2.

A spatio-temporal map of tissue temperatures along the axis of the imaging transducer.

3.1.1.

Backscattered Signal Model

A widely accepted model for ultrasound backscattered signal from tissue media expresses the received signal as a superposition of reflections from N scatterers within the resolution cell of the transducer:

where A_i is the strength of the ith scatterer assumed to be independent of frequency, p_z(t)wo-way propagation del is the basic ultrasound wavelet emitted from the transducer with z (range) dependence to account for attenuation, diffraction, etc., and τ_i is the two-way propagation delay to the ith scatterer . A simplified version of the spectrum of the bacscattered signal can be given by:

where we have explicitly accounted for the attenuation (a is the attenuation coefficient in dB/cm).

3.1.2.

Autoregressive Modeling

Autoregressive (AR) modeling of ultrasound backscattered data has proved to be extremely beneficial in many applications, e.g., Doppler, mean scatterer spacing estimation, attenuation measurements, etc. The pth order AR model for the received signal given in Equation 7 is given by:

where x(n) is the sampled version of x(t) at t = nT, a_i are the AR filter coefficients to be determined from the time series, and ϵ(n) is white Gaussian noise sequence with zero mean and variance .

A typical use of the AR model calls for a model order p such that ϵ is truly a white noise sequence. However, for many applications, a low-order AR model is preferred for computational simplicity and robustness. Furthermore, for values of p equal to 2 and 3 (AR2 and AR3), the coefiicients of the AR model take on physical significance with intuitive appeal. Regardless, once the model order has been determined the filter coefficients can be determined by many techniques based on second oder or higher order statistics. Furthermore, algorithms for determining the optimal filter coefficients have been developed for cases when the underlying model is stationary or nonstationary. The interested reader is referred to several textbooks on the subject.

3.1.3.

AR Parameter Estimation

In an AR2 model, the coefficients a₁ and a₂ are direcly related to the reflection coefficients as well as the center frequency of the spectrum of the backscatter. The latter is an important parameter in determining the attenuation coefficients. The attenuation coefficient is known to change dramatically when the tissue is subjected to coagulative necrosis.²⁹ Assuming that α is a linear function of frequency, α = α₀f, then it can be shown that

where σ_s is related to the transducer bandwidth and f_o is the center frequency of the backscattered signal. It can also be shown that

where f_s is the sampling frequency.

4.3.1.

Data collection system

Imaging data was obtained by driving the transducer in pulse-echo mode. The transmitted pulses were generated by a pulser/receiver (Panametrics 5052-PR). The received echoes were first amplified by 30dB using an ultrasound signal amplifier (NTR pfs-005a). The transmit and receive circuitry were connected to the array elements through a VXI programmable 64×4 relay matrix (Tektronix VX4380). Two different data acquisition systems were used to digitize the data (Fig. 7).

Figure 7.

Schematics of the circuitry used in the imaging experiments. When the custom made gated digitizer is used, the digital oscilloscope is replaced by the blocks in dashed lines.

The first digitizer consists of a digital signal analyzer (Tektronix DSA-601), the captured data being transfered to the controlling workstation through a GPIB interface. The duration of the data collection in this case is limited by the speed of the digitizer, a complete set of 64×64 A-lines taking about 3 hours.

In order to reduce the total duration of the data collection, a gated digitizer was implemented. The data is digitized using a Burr Brown ADS-801 12-bit fast-A/D operated at 16 MHz (data sampling frequency). Gating is implemented using conventional fast TTL logic. The digitized data is stored in a dual port 128 Mb memory buffer, and transfered at slower rate to the controlling workstation. The duration of the acquisition is limited in this case by the dead time of the relay matrix after each switching (approximately 40 ms), a complete set of 64×64 A-lines being acquired in about 3 minutes. The current prototype implementation of this circuit provides data with lower SNR than the data obtained using the Tektronix digital signal analyzer. This increases the level of the background noise in the reconstructed image, ultimately limiting the image dynamic range. Nevertheless, the system is still capable of providing good image quality, as seen in Fig. 12.

4.3.2.

Synthetic aperture imaging

Images were obtained using a full synthetic aperture technique.³¹ In this case the system is in focus at every pixel of the image, in transmit and receive. Furthermore, in order to reduce side-lobes, apodization of the transmit and receive apertures was applied by weighting the echo responses with a Kaiser window³⁴ (β = 3.25). The image pixel at coordinates (x_p, z_p) was therefore computed by (Fig. 8):

Figure 8.

Coordinate system used in the synthetic aperture imaging system.

where i is the transmit element index, j is the receive element index, A_i is the transmit apodization weight at element i, B_j is the receive apodization weight at element j, R_ip and R_jp are the distances from the transmit and receive elements, respectively, from the image pixel, c is the speed of sound in the medium being imaged, and S_i,j(t) is the base-band converted echo acquired when transmitting with element i and receiving with element j. The images were first sampled in the R × sin(θ) domain, and then scan converted into a Cartesian spatial sampling.

4.4.1.

Multiple-wire phantom

A multiple-wire phantom consisting of thirteen nylon wires spaced by 10 mm along a horizontal line was used to characterize the point spread function of the imaging system described above (Fig. 9). The phantom was aligned so that the middle wire was located at the geometric center of the spherically-shaped therapeutic transducer. The phantom and transducer were immersed in water to provide acoustic coupling. A full 64×64 data-set was collected using the setup represented in Fig. 7, the data being digitized by the Tektronix digital signal analyzer.

Figure 9.

Cross section of the multiple wire phantom. The 13 wires are indicated as well as the Plexiglas frame that holds the wires. Dimensions are indicated in millimeters.

The pulse-echo synthetic aperture reconstructed image is shown in Fig. 10 at dynamic ranges of 20, 30, 40, and 50 dB, where all images are referenced to the peak value in the reconstructed region. At 20 dB seven central wires are seen. As this transducer is highly focused, and no compensation for element directivity was applied during reconstruction, objects closer to the geometric center of the transducer appear to be brighter than objects farther away. Furthermore, as the aperture is undersampled grating lobes become visible when imaging off-axis (image at 40 and 50dB). At 50 dB eleven of the thirteen wires are seen.

Figure 10.

Image of the multiple wire phantom shown at different dynamic ranges. The vertical direction corresponds to axial depth in the interval from 80 mm to 120 mm. The horizontal direction corresponds to lateral from -60 mm to +60 mm.

The results in Fig. 10 indicate that this system provides good image reconstruction of objects close to the transducer geometric center. However, the quality of the image degrades as objects off-center are considered. Although this image degradation is certainly not appropriate for diagnostic systems, in the case of therapeutic guidance the most relevant information is contained in the region undergoing treatment. Typically, this corresponds to the region close to the center of the therapeutic transducer, where the image quality is better.

4.4.2.

Speckle-generating phantom

In order to characterize the performance of this system with speckle media, a tissue-mimicking phantom was designed, its dimensions being given in Fig. 11. The phantom is made of rubber material (plastic for making plastic worms, M-F Manufacturing Co., Fort Worth, TX), with embedded polystyrene scatterers (Amberlite IR-120plus, Sigma Chemical Co., St. Louis, MO). The scatterers are homogeneously distributed within the phantom, except for a 12.7 mm diameter cyst region, that was filled with water. The phantom and the transducer were immersed in water to provide acoustic coupling. The phantom was aligned so that the cyst region corresponded to the geometric center of the therapeutic transducer, as indicated in Fig. 11. A full 64×64 data set was collected using the experimental setup shown in Fig. 7, the custom gated digitizer being used in this case.

Figure 11.

Dimensions of the speckle generating phantom and its position relative to the transducer. Dimensions are indicated in millimeters.

Figure 12.

Image of the speckle generating phantom shown at different dynamic ranges. The vertical direction corresponds to axial depth in the interval from 70 mm to 120 mm. The horizontal direction corresponds to lateral from -50 mm to +50 mm.

The synthetic aperture image of the phantom is shown at four different dynamic ranges, namely 20, 30, 40 and 50 dB, referenced to the peak value detected within the imaged region. The image shows the cyst region, illustrating the performance of the system in the vicinity of the focal region. The strong reflection at the top and bottom edges of the cyst are caused by the acoustic impedance mismatch between the rubber material and water.