Photoacoustic Doppler measurement of flow using tone burst excitation

In this paper a novel technique for flow measurement which is based on the photoacoustic (PA) Doppler effect is described. A significant feature of the proposed approach is that it can be implemented using tone burst optical excitation thus enabling simultaneous measurement of both velocity and position. The technique, which is based on external modulation and heterodyne detection, was experimentally demonstrated by measurement of the flow of a suspension of carbon particles in a silicon tube and successfully determined the particles mean velocity up to values of 130 mm/sec, which is about 10 times higher than previously reported PA Doppler setups. In the theoretical part a rigorous derivation of the PA response of a flowing medium is described and some important simplifying approximations are highlighted. ©2010 Optical Society of America OCIS codes: (110.4153) Motion estimation and optical flow; (280.2490) Flow diagnostics; (110.5125) Photoacoustics; (170.5120) Photoacoustic imaging. References and links 1. L. V. Wang, Photoacoustic imaging and spectroscopy (CRC Press, 2009) 2. M. Xu, and L. V. Wang, “Photoacoutsic imaging in biomedicine,” R. Sci. Inst. 77(4), 041101 (2006). 3. H. Fang, K. Maslov, and L. V. Wang, “Photoacoustic Doppler effect from flowing small light-absorbing particles,” Phys. Rev. Lett. 99(18), 184501 (2007). 4. H. Fang, K. Maslov, and L. V. Wang, “Photoacoustic Doppler flow measurement in optically scattering media,” Appl. Phys. Lett. 91(26), 264103 (2007). 5. H. Fang, and L. V. Wang, “M-mode photoacoustic particle flow imaging,” Opt. Lett. 34(5), 671–673 (2009). 6. Y. Wang, D. Xing, Y. Zeng, and Q. Chen, “Photoacoustic imaging with deconvolution algorithm,” Phys. Med. Biol. 49(14), 3117–3124 (2004). 7. S. M. Blinder, “Delta functions in spherical coordinates and how to avoid losing them: fields of point charges and dipoles,” Am. J. Phys. 71(8), 816–818 (2003). 8. W. R. Brody, and J. D. Meindl, “Theoretical analysis of the CW doppler ultrasonic flowmeter,” IEEE Trans. Biomed. Eng. 21(3), 183–192 (1974). 9. G. Guidi, C. Licciardello, and S. Falteri, “Intrinsic spectral broadening (ISB) in ultrasound Doppler as a combination of transit time and local geometrical broadening,” Ult. Med. Biol. 26(5), 853–862 (2000). 10. A. Sheinfeld, E. Bergman, S. Gilead, and A. Eyal, “The use of pulse synthesis for optimization of photoacoustic measurements,” Opt. Express 17(9), 7328–7338 (2009).


Introduction
Photoacoustic (PA) imaging has attracted much attention in the last decade, especially in the context of biomedical applications, where advantages of both optical and ultrasound imaging can be exploited to obtain improved performance and new functionalities [1,2].PA techniques are based on measuring the acoustical signal which is generated due to the absorption of modulated light in a tested medium.They offer high contrast and functional information based on differences in the optical absorption properties of the imaged areas while enjoying the high spatial resolution and weak acoustical scattering which characterize ultrasound bioimaging.
A few recent studies [3][4][5] suggested the application of the PA effect for blood flow measurement, either via determination of the Doppler shift [3,4] or through implementation of the ultrasound M-mode approach [5].However, both of these techniques demonstrated limited maximum measurable velocity, which was reported to be 10 mm/sec in the first method and 1 mm/sec in the second.Also, the PA Doppler demonstrations so far were based on continuouswave (CW) modulation and lock-in detection, which cannot resolve the axial position of the detected moving particles and feature inherent tradeoff between sensitivity and maximum measurable Doppler frequency shift.In this paper we present, for the first time, a PA Doppler setup which uses burst excitation, and demonstrate the detection of velocities as high as 130 mm/sec.These characteristics are enabled via the use of a particularly flexible external modulation scheme and highly sensitive heterodyne detection.The use of heterodyne detection alleviates the tradeoff between sensitivity and maximum measurable Doppler shift since, in contrast with a lock-in measurement scheme, the receiver center frequency can be varied and is not limited to DC.It is also shown that the velocity induced Doppler shift is accompanied with a proportional Doppler broadening.This broadening, following proper calibration, can be used as an alternative method for measurement of the particles mean velocity.
The theory of flow measurement via the Doppler effect has been extensively treated in the context of ultrasound imaging theory.Several inherent attributes make the PA Doppler flow measurement different than the ultrasound Doppler flow measurement.While the latter is based on scattering, typically with the emitter and receiver collocated, the former can be thought of as based on the emission of a non-stationary distributed source.In the first paragraph we introduce a rigorous derivation of the PA response of a flowing medium to sinusoidal excitation and highlight several important simplifying approximations.The derivation yields expressions for slowly varying time-dependent amplitude and frequency modulations which contribute to the homogeneous broadening of the PA response.It is also shown that the magnitude of the PA frequency response depends linearly on the spatial Fourier component of the PA efficiency distribution at the wavenumber of the corresponding acoustical wave.

Theory
Consider a fluid, with inhomogeneous photoacoustic properties, flowing through a cylindrical volume in space (Fig. 1).The fluid is optically excited by temporally modulated light with uniform spatial distribution.Generalizing the well known Green-integral solution to the PA wave equation [6] we express the pressure produced due to the photoacoustic effect, as detected by an ideal ultrasonic detector at the origin as: where r = r and we assume that ( , ) H t r ɶ describes both the optical excitation and the physical properties of the excited medium, namely its thermal volume expansion, specific heat and absorption properties.For simplification, it is also assumed that the speed of sound in the fluid, in the frame of reference of the ultrasonic detector, is negligibly affected by the motion of the fluid and is equal to the speed of sound c in the surrounding medium.At zero flow rate ( , ) H t r ɶ can be represented as where ( ) I t is the timedependent optical excitation and ( ) A r is a function of position describing the PA efficiency distribution.Note that according to our definition of ( ) A r it is zero in the surrounding medium and it includes three physical parameters that contribute to the PA effect, namely, the absorption, thermal volume expansion and specific heat.In the case where the fluid is moving, the PA efficiency distribution becomes also a function of time and the integrand in Eq. ( 1) can be separated into two terms as follows: The first term describes the PA response to the modulated optical excitation.The second term represents a PA effect which stems from the temporal dependence of the absorption medium.Interestingly, this effect would lead to PA response even with a temporally constant (unmodulated) illumination.It is shown in the appendix, however, that for the parameters of a typical flow measurement this term can be neglected.
To see the effect of flow on the PA response we start by considering a single point absorber moving at a constant velocity v .The absorber is represented by a 3-dimentional delta function, ( ) , where a is the PA efficiency coefficient of the particle, v is the absorber's velocity vector and 0 r is its position at 0 t = .The PA response to this point absorber will be used later to construct the response to a general distribution.In the retardedtime formalism of Eq. ( 2), the 3-dimentional delta function becomes ( )  ) ( ) It can be shown through exact calculation that use of the approximated form amounts to neglecting terms of the order of 2 2  v c in the expression for the Doppler shift and the PA amplitude.Adopting a spherical coordinate system and using the appropriate delta-function representation [7], Eq. (2) becomes: where r v describes the radial motion of the point absorber, and φ are its time-dependent azimuth and elevation respectively.The explicit dependence of these functions on time is not important since the other terms in the integral are independent of θ and φ .Equation ( 3) can be easily integrated yielding: The complex PA response for continuous sinusoidal excitation at frequency 0 f can then be found: where 0 I is the peak-to-peak amplitude of the optical instantaneous power and the actual response is given by ( )

{ }
Re p t .Clearly this response is no longer a "pure" sinusoidal signal.
First, its amplitude is weighted by the varying function of time, ( ) This term leads to broadening of the response in the frequency domain of approximately (which was limited by 15 Hz in our experiment), however, given the targeted velocities, the typical distance between the particle and the detector and the other broadening mechanisms, this term has a very minor effect on the measurement.Second, the instantaneous frequency is also varying with time and is Doppler shifted by: where v = v and ( ) and v .This result, which reproduces the expression for the PA Doppler frequency given by Fang et al. [3], manifests another broadening mechanism of the Doppler peak that occurs if the motion is not strictly radial.Typically and the detection of the Doppler shifted tone can be implemented by heterodyne spectral analysis at the vicinity of 0 f .In the case of a multitude of particles in a fluid or a continuous flowing substance it is required to integrate the response in Eq. ( 5) over all contributions.Assuming uniform rectilinear flow, whose only effect on the initial distribution is simple translation, the spatial absorption distribution can be expressed as: where the integration is effectively bounded to the cylindrical volume where ( ) . The total PA response to CW excitation then becomes: Some physical insight regarding this result can be attained when considering the case of radial motion.In this case the flow is directed towards the detector so that ( ) 0 p r t r v t = − .The PA response in this case becomes: where is the Fourier transform of ( ) 9) is a manifestation of ( ) p t being a sum of coherent contributions.For an efficient effect it is required that the spatial spectra of ( ) A r would have a strong component at the wavenumber of the acoustical wave.This is similar to the condition for scattering of a Bragg grating but the required "pitch" in photoacoustics is doubled since no round trip is involved.In realistic setups the spectral response near 0 PAD f f + is further broadened due to several other mechanisms, similar to those known from Doppler ultrasound [8,9].These mechanisms can be classified into two main groups: homogeneous and inhomogeneous broadening.The finite dimensions and spatial inhomogeneities of the excitation beam is a major cause for homogeneous broadening.As a particle traverses the illuminated region its PA response is modulated by the relatively slow variations of the local intensity of the light.The induced broadening due to this effect can be estimated to be / where d is the characteristic length of the spatial inhomogeneities of the optical beam or the dimension of the beam in the direction of the flow if the beam is homogenous.Another homogeneous broadening factor stems from the finite size of the ultrasound detector.This leads to a range of angles ( ) t γ that contribute to the PA response, each with a different Doppler shift, as described by Eq. ( 6).It is helpful to define the mean angle γ which is a weighted average of ( ) t γ over all the PA-medium and over all detector points.Inhomogeneous broadening mechanisms result from differences in the Doppler shifts of different particles.As the particles have different sizes and shapes and as they flow at different radial distances from the axis of the tube their velocities are distributed over a wide range.As an estimate to the distribution of particles velocities we recall that in a laminar flow regime, which applies for the velocities measured in our experiment, the velocity of a flowing substance follows a parabolic shape: where r is the radial distance from the tube axis, R is the tube radius, F is the flow rate, A is the tube cross-section area and avg v is the mean velocity.It is interesting to note that all these broadening mechanisms show a linear dependence on the mean velocity.This observation, which was also experimentally observed, can be used as an alternative method for estimation of the particles velocity.Deviation from laminar regime may induce broadening that is nonlinearly dependent on the velocity, however our experimental results did not show such dependence.

Experimental setup
To enable extraction of axial position information in addition to the Doppler frequency shift, the continuous sinusoidal modulation is replaced by sinusoidal bursts with rectangular envelopes as expressed by: Here c T is the Pulse Repetition Interval (PRI) and r T is the burst duration.The main effect of the burst operation on the spectrum of the PA response is the appearance of spectral replicas of the CW-spectrum at periodic frequency points separated by 1 c T .Consequently, the PRI must be chosen to be smaller than 1 ( 2) PAD f to avoid aliasing.An equivalent consideration in the time-domain leads to the constraint max / PRI r c > , where max r is the distance between the detector and the farthest point in the measured volume, to allow enough time for the PA signal to reach the detector before the next excitation burst starts.
The optical part of the experimental setup (Fig. 2) comprised a fiber-coupled Tunable Laser Source (TLS) in the range of 1510-1620nm, a polarization controller, a Lithium-Niobate Electro-Optic Modulator (EOM) driven by an Arbitrary Waveform Generator (AWG) and an Optical Erbium Doped Fiber Amplifier (EDFA) with a cleaved fiber output.This particular setup which was used to demonstrate the proposed method operated in the fiber optic C-band (around 1550nm), and took advantage of components available for fiber-optic communications.This spectral range, however, has the disadvantage of relatively high optical absorption of water, which limits the penetration depth in biological tissues.Nevertheless, the growing interest in emerging non-telecom photonic applications has led to the development of comparable components in other wavelength ranges, in particular in the spectral range of 980-1150nm, which is more suitable for biological applications.The external modulation of the optical excitation allowed a precise burst modulation at any desired pulse rate and duty-cycle and in frequencies up to 100 MHz.As will be described below, a judicious biasing approach was required to ensure high extinction of the light between bursts.A similar setup was previously proposed for the application of sensitivity and responsivity enhancement of PA measurements [10].The modulated light, at mean power of 100mW, exited the cleaved fiber and was directed, without any collimation, to a transparent silicon tube (inner diameter: 3mm) which was coiled to form a loop as described in Fig. 2. The illuminated region was set to a length of 4mm via a circular aperture and was centered on one of the vertical sections of the loop.This configuration enabled a relatively small angle between the direction of flow and the transducer's acoustical axis, thus guarantying a relatively large Doppler shift, according to Eq. ( 6).The acoustical detection was implemented using a commercial ultrasound transducer with a center frequency of 1MHz (Olympus IR-0108-S-RU) followed by a low-noise preamplifier (Olympus 5660C).The detector was placed below the illuminated part with its symmetry axis slightly tilted with respect to the vertical axis.The gap between the transducer active aperture and the tube was around 15mm wide and it was filled with ultrasound gel for ensuring good acoustical coupling.The output of the preamplifier was split into two branches where one branch was connected to a RF spectrum analyzer and the other to a digital storage oscilloscope (DSO).The spectrum analyzer was used to measure the Doppler shift and was set to its maximum resolution bandwidth (1Hz) and a span of 100-500Hz.The DSO enabled measurement of the axial position of the PA generating region.
A syringe pump (Fusion-200, Chemyx) with a 50cc syringe was used for generating flow of up to 55 ml/min (either withdrawal or infusion) in the tube, giving a maximum mean velocity of 130 mm/sec for the given tube diameter.The fluid in the tube was made from an untreated graphite powder with typical particles dimensions of 37-149µm (Activated charcoal C3345, Sigma-Aldrich) suspended in deuterium oxide (D 2 O) with volume fraction φ ≈ 5%.
The choice of D 2 O was due to its relatively low absorption (around 1 cm −1 ) in the relevant spectral range.The absorption of the suspension was estimated to be around 5 cm −1 by measuring its PA response and using the PA response of double distilled water for calibration.
To produce the optical excitation the EOM was biased, through the DC port of the bias-T, at its minimum transmission point and a train of sinusoidal bursts with rectangular envelope was fed to the AC port of the bias-T.The bursts width and the PRI were chosen to be r T = 5µsec and c T = 20µsec respectively, to avoid aliasing in the time and frequency domain as described earlier.The frequency of the sinusoidal modulation was set to 0.5MHz which yielded at the output of the EOM optical pulses with modulated instantaneous power centered at 1MHz.The PA response to the burst excitation was measured in both the time-domain (using the DSO) and the frequency domain (using the spectrum analyzer) for flows between 1 to 55 ml/min in the infusion and withdrawal modes of the pump, yielding a wide range of velocities between 2.5 and 130 mm/sec.The signals were averaged in both instruments over the entire infusion/withdrawal time for a syringe volume of 50ml.

Experimental results
The excitation bursts and an example of the response in the time-domain are shown in Fig. 3.The temporal width of the generated PA bursts is 4.5µsec (FWHM) representing a spatial resolution of approximately 7mm.It can be seen from this figure that the PA response arrives 12µsec after the excitation.Assuming the speed of sound in the ultrasonic gel is approximately 1500m/s this delay corresponds to a distance of ~1.8cm and agrees well with the measurement geometry.f are due to stationary absorbers such as the tube itself or particles that were deposited on the innr walls of the tube.The spectra of moving particles appear as wings at the sides of the dominant peak at the modulation frequency.The dependence of both the center frequency and the spectral width of these wings on the flow rate and its direction is evident.The plots of the spectrum at higher rate are noisier due to the shorter averaging time.In order to estimate the particles velocity, the detected spectra, excluding a portion of 7 Hz around 0 f , was fitted to a shifted Gaussian, using least mean square (LMS) minimization.Figure 5 shows the measured Doppler shift for the different flow rates, compared to the theoretical shift as was calculated by Eq. ( 6).The sound velocity in the fluid was taken as A c = 1500 m/s and γ was found from calibration to be γ =60° which is in agreement with the geometry of our setup.The strong peak at 0 f , whose width is determined by the resolution bandwidth of the spectrum analyzer, limited the minimum measurable velocity to roughly 2.5 mm/sec.This limitation may become less restrictive by choosing an optical wavelength in which the particles has significantly higher absorption than the background.Due to the limitations of the syringe pump that was used, the maximum measurable velocity was not characterized.The maximum mean velocity that was measured was 130 mm/sec, but it is expected that higher velocities are achievable by the proposed method.
In addition to the mean velocity, the LMS fitting yielded the spectral width of the Doppler-shifted signal, calculated as the FWHM of the fitted Gaussian (Fig. 6).It can be seen that the spectral width is linearly dependent on the particles mean velocity and the proportion coefficient is around 2.5 Hz/(mm/sec).Using Eq. ( 10) it can be shown, via numeric simulation, that the expected broadening proportion coefficient calculated from the distribution of the velocities in laminar flow and transit-time broadening for a homogeneous optical beam is smaller by a factor of 4 from the measured coefficient.Similarly, the homogeneous broadening due to the finite size of the ultrasound detector and the rotation of the relative velocity vector between the particle and the detector cannot account for the observed proportion coefficient.The increased measured spectral broadening is attributed to spatial inhomogeneities of the excitation beam.Such inhomogeneities are expected due to the particles aggregates that were deposited on the inner walls of the tube and partially obstructed and scattered the excitation beam.The characteristic length of these spatial inhomogeneities can be estimated via / and the experimental data to be around 0.5mm for γ =60° (red dotted).Positive rates were measured for infusion (where the flow was towards the transducer) and negative for withdrawal.Inset includes zoom-in on the low velocities section.

Conclusions
In this paper the PA response of a flowing medium to sinusoidal excitation was derived and a new setup for photoacoustic Doppler flow measurement was presented and experimentally characterized.Besides reproducing the expression for the Doppler frequency shift from the basic wave equation, the theoretical analysis yielded also expressions for two inherent broadening mechanisms and showed the dependence of the PA response on the spatial spectra of the absorbing medium.Experimentally, the use of heterodyne detection yielded measurements of velocities up to 130 mm/sec which is 10 times higher than previously reported PAD setups.Moreover, for the first time in Doppler PA measurements, burst excitation was used.It was implemented via a particularly flexible external modulation scheme and enabled extraction of information about the position of the PA sources in addition to their velocities.Complete two-dimensional characterization in the time and frequency domains can be achieved through time-gating the input to the spectrum analyzer or via using a DSO with a higher dynamic range.In addition, implementing the setup using optical excitation at another wavelength range may reduce the background signal, allowing measurements of lower velocities and increasing the penetration depth.Using a detector at higher frequencies would improve spatial resolution, which was around 7mm in the presented setup.These improvements are currently under investigation.The advantages of PA imaging can be harnessed to in-vivo flow measurements, where the combination of relatively low scattering of the ultrasound signal with the high optical contrast may give superiority over purely optical or purely ultrasonic flow measurements techniques.

Fig. 1 .
Fig. 1.Illustration of PA flow measurement of inhomogeneous fluid flowing in a cylindrical volume

Fig. 3 .
Fig. 3.The optical burst excitation (black) and the resulting PA response (pink) in the time domain Typical spectra of the PA-response for 2 different rates in both infusion and withdrawal are shown in Fig. 4. The strong peaks at 0f are due to stationary absorbers such as the tube itself or particles that were deposited on the innr walls of the tube.The spectra of moving particles appear as wings at the sides of the dominant peak at the modulation frequency.The dependence of both the center frequency and the spectral width of these wings on the flow rate and its direction is evident.The plots of the spectrum at higher rate are noisier due to the shorter averaging time.

Fig. 4 .
Fig. 4. Example of spectra for infusion (blue) and withdrawal (pink) for rates of 20 ml/min (left plot) and 40 ml/min (right plot)

Fig. 6 .
Fig. 6.Spectral width (FWHM) vs. average particles velocity.Inset includes zoom-in on low velocities section ) and therefore expected to be generally negligible.