Patent Publication Number: US-10327739-B2

Title: Measuring tissue shear wave properties using one excitation pulse and excitation pulses having different widths and shapes

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
     This application claims the benefit of U.S. Provisional Patent Application No. 62/143,017, filed Apr. 3, 2015. The disclosure of U.S. Provisional Patent Application No. 62/143,017 is incorporated herein by reference. 
    
    
     BACKGROUND 
     Quantitative ultrasound elastography measures tissue viscoelasticity by applying ultrasound radiation force to a tissue region, measuring the propagation properties of induced shear waves, such as displacements and propagation velocities (group velocity or phase velocity), and calculating viscoelasticity of the tissue. These technologies include ultrasound vibrometry (U.S. Pat. Nos. 7,753,847, 7,785,259, 8,659,975, 8,602,994), ARFI (US 20050215899), Supersonic shear Imaging (SSI) and Shear wave Spectroscopy (SWS) (U.S. Pat. No. 8,150,128), and other similar techniques such as Shear Wave Elasticity Imaging (SWEI), Spatially Modulated Ultrasound Radiation Force (SMURF), Crawling Wave Sonography (CWS). 
     Above said methods assume that tissue response to excitation pulses are linear, so that the measurements of tissue property such as group velocity and phase velocity of shear wave are independent to excitation pulse widths and shapes, as the said velocities are used to calculate the tissue viscoelasticity for charactering tissue shear property. 
     Above said methods do not measure tissue nonlinearity, nor consider the impact of the nonlinearity to the accuracy of the estimates of the tissue viscoelasticity. 
     Above said methods require the detection of shear wave at multiple locations to estimate shear wave velocities in order to calculate the viscoelasticity of tissue. 
     SUMMARY 
     This work finds that measured propagation velocities of shear waves are not independent to excitation pulse widths and shapes in biological tissues such as swine livers. It is partially due to the complexity of biological structures of tissue including nonlinearity, anisotropy, and non-uniformity, etc. For a tissue region, this disclosure finds that the measured tissue shear property of the induced shear waves including group velocities and phase velocities are different for different excitation pulses. Thus, the calculated viscoelasticity of the tissue are different for different vibration pulses in the same tissue region. Therefore, the estimated tissue viscoelasticity using the said prior art of elastography is pulse dependent. 
     This disclosure discloses the method and apparatus to characterize tissue shear property using different vibration pulses having difference widths and shapes. 
     According to one aspect of this disclosure, the detection location can be at one point at a time for imaging two-dimensional or three-dimensional tissue shear wave property investigated by different excitation pulses. The property includes nonlinear magnitude variations, nonlinear phase variations, and combination parameters based on the responses excited by different pulses, which are related to the tissue shear property and pathological statues. 
     According to another aspect of this disclosure, the induced shear wave can be also detected at multiple locations in the tissue region for calculating different group velocities and phase velocities using different excitation pulses. For example, excitation pulses having pulse widths of 50 μs and 100 μs are used to induce shear waves, respectively. This disclosure estimated the group velocities and phase velocities of shear wave induced by the said two different excitation pulses. The measured velocities and its dispersion with different pulses are used to characterize tissue shear properties, which are used to characterize different types of tissue and pathological statues of the tissues. 
     Another aspect of this disclosure discloses a method to measure tissue nonlinearity by using a single pulse or a sequence of the same pulses for a measurement, while as the measurement is repeated for different pules having different widths and shapes. 
     Another aspect of this disclosure is to correct the viscoelasticity measured by the said priori elastography technologies, by providing a standard of specified pulse widths. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
         FIG. 1  shows a radiation force pi(t) generated by one of pulses having different widths is applied to a tissue region to induce a shear wave yi(t, x), i=1, 2, . . . M, and x is a location according to an aspect of the disclosure. The process is repeated for another pulse till all shear waves induced by pulses having different widths are measured for estimating the tissue shear property including nonlinearity. 
         FIG. 2  shows a system diagram of tissue shear property characterization using ultrasound radiation force induced by excitation pulses having different widths and shapes according to an aspect of the disclosure. 
         FIG. 3  shows phase velocities of shear waves in a swine liver phantom according to an aspect of the disclosure. It shows that the phase velocity of the induced shear wave at a given frequency is not a constant, but varies with the pulse width. This dispersive property of phase velocity with different pulse widths represents tissue pathological statues. 
         FIG. 4  shows group velocities of shear waves in a swine liver phantom according to an aspect of the disclosure. It shows that the group velocity of the induced shear wave is not a constant, but varies with the pulse widths. This dispersive property of group velocity with different pulse widths represents tissue pathological statues. 
         FIG. 5  show estimated elasticities of a swine liver phantom. It shows that the estimated elasticity of a swine liver is not a constant, but varies with the pulse width according to an aspect of the disclosure. This dispersive property of estimated elasticity with different pulse widths represents tissue pathological statues. 
         FIG. 6  show estimated phase velocities of a swine liver phantom versus frequency according to an aspect of the disclosure. It shows that the estimated elasticity of a swine liver is not a constant at a frequency, but varies with the pulse widths. 
         FIG. 7  shows a procedure to measure tissue nonlinearity, viscoelasticity, and shear velocities by using a single pulse or different excitation pulses at one or more locations according to an aspect of the disclosure. 
         FIG. 8  shows the vibration beams and detection beams aimed at different locations according to an aspect of the disclosure. 
         FIG. 9  shows the vibration beams and detection beams aimed at the same location according to an aspect of the disclosure. 
         FIG. 10  shows the vibration beams and detection beams are aimed at different locations according to an aspect of the disclosure; detections are done at multiple locations in multiple dimensions in a region. 
     
    
    
     DETAILED DESCRIPTION 
     When ultrasound radiation force is applied to viscoelastic tissue, the induced tissue shear motion in a tissue region can be described by Newton&#39;s second law of motion: 
                       ∑     k   =   1     3     ⁢           ⁢       ∂     σ   jk         ∂     x   k           =     ρ   ⁢         ∂   2     ⁢     d   j         ∂     t   2         ⁢     (       j   =   1     ,   2   ,   3     )               (   1   )               
where (x 1 , x 2 , x 3 ) are Cartesian coordinates, (σ 11 , σ 12 , σ 13 ) are x 1 , x 2 , x 3 -components of the stress per unit cross the x 2 -x 3  plane, analogous to other components, d j  is the displacement of the motion in x j  direction, ρ is density. This equation states that the net internal forces per unit area are derivatives of the stress tensors in space, and the net forces induce the accelerated motions at location (x 1 , x 2 , x 3 ).
 
     The wave propagation speeds of the wave components described by (1) are useful for studying the mechanical properties of the tissue, such as shear elasticity and viscosity of tissue. 
     Considering a harmonic motion of the shear wave of which the tissue displacement is parallel to the compressional wave front and the propagation direction is perpendicular to the compressional wave front in an isotropic viscoelastic medium, (1) is reduced to: 
                             ∂   2     ⁢     d   1         ∂     x   k   2         +       h   2     ⁢     d   1         =   0     ,       ∀   k     =     2   ⁢           ⁢   or   ⁢           ⁢   3               (   2   )               
where d 1  is the displacement of the harmonic motion, x k  is in the direction that is perpendicular to the wave front of the compressional wave, h is the wave number of the shear wave in the x k  direction. The solution of (2) for a motion having only one single harmonic is:
 
 d   1 ( x,t )=Re{ De   j(ωt−hx)   }De   −h     I     x  cos(ω t−h   r   x )  (3)
 
where h r  and h I  are real and imaginary parts of h, respectively, D is the initial amplitude of the tissue displacement, ω is the angular frequency of the harmonic motion, and x is the travel distance. The velocity of the shear wave can be obtained by taking time derivative of (ωt−h r x), which is:
 
     
       
         
           
             
               
                 
                   
                     
                       v 
                       s 
                     
                     ⁡ 
                     
                       ( 
                       t 
                       ) 
                     
                   
                   = 
                   
                     ω 
                     
                       h 
                       r 
                     
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     The said priori methods measure the velocity v s (ω) and solve the tissue viscoelasticity that is related to h r , which is dependent to the tissue viscoelasticity model. Tissue viscoelasticity can be modelled by many different models. The simplest one is:
 
σ( t )=μ 1 ε( t )  (5)
 
where σ is stress, ε is strain, μ 1  is called elasticity. Thus,
 
                       v   s     ⁡     (   ω   )       =       ω     h   r       =         μ   1     ρ                 (   6   )               
where ρ is the density of tissue. Another popular model with one additional element for viscoelastic tissue is called Voigt model:
 
                     σ   ⁡     (   t   )       =         μ   1     ⁢     ɛ   ⁡     (   t   )         +       μ   2     ⁢       d   ⁢           ⁢     ɛ   ⁡     (   t   )         dt                 (   7   )               
where μ 2  is called viscosity. For a harmonic motion:
 
σ/ε=μ 1   +jωμ   2   (8)
 
For the Voigt medium:
 
 h =√{square root over (ρω 2 /(μ 1   +iωμ   2 ))}  (9)
 
Thus, the phase velocity of the shear wave in Voigt tissue:
 
     
       
         
           
             
               
                 
                   
                     
                       v 
                       s 
                     
                     ⁡ 
                     
                       ( 
                       ω 
                       ) 
                     
                   
                   = 
                   
                     
                       ω 
                       
                         h 
                         r 
                       
                     
                     = 
                     
                       
                         
                           2 
                           ⁢ 
                           
                             ( 
                             
                               
                                 μ 
                                 1 
                                 2 
                               
                               + 
                               
                                 
                                   ω 
                                   2 
                                 
                                 ⁢ 
                                 
                                   μ 
                                   2 
                                   2 
                                 
                               
                             
                             ) 
                           
                         
                         
                           ρ 
                           ⁡ 
                           
                             ( 
                             
                               
                                 μ 
                                 1 
                               
                               + 
                               
                                 
                                   
                                     μ 
                                     1 
                                     2 
                                   
                                   + 
                                   
                                     
                                       ω 
                                       2 
                                     
                                     ⁢ 
                                     
                                       μ 
                                       2 
                                       2 
                                     
                                   
                                 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   10 
                   ) 
                 
               
             
           
         
       
     
     It can be showed that the amplitude of the shear wave in distance is also a function of tissue viscoelasticity, and the similar equations to (6) and (10) can be found that relate amplitude attenuations to the viscoelasticity. 
     For other tissue viscoelastic models, the velocity related to wave number h and viscoelasticity can be found by using (4). They can be found in literatures and are not listed here. 
     The priori technologies calculate viscoelasticity based on the velocities or amplitudes of the induced shear waves over a distance. Thus, the detection of shear wave at multiple locations are required. This limits its applications for tissue points and tissue imaging. 
     In theory, the phase velocities and group velocity of shear wave in the tissue region represent tissue property, and they should not be dependent to the excitation pulses, which is true if the tissue response to the radiation force is linear. In this case, for Voigt tissue, the velocity is given by (10); for tissue having only elasticity, the velocity is given by (6). The shape of the applied radiation force may be changed over a distance due to attenuation, but this change has no impact to propagation phase velocity. 
     This work finds that the propagation velocities (both group and phase velocities) and the shapes of the induced shear waves may be changed in the tissue. These changes are functions of both the tissue property and the shapes of applied radiation forces. When excitation pulses are used to generate the ultrasound radiation force in a complex tissue region which may be nonlinear, the pulse width has impact to the measured propagation velocities. Thus, the calculated viscoelasticity is not unique as described by (6) and (10) and other equations. This situation is also true for the radiation force generated by continue waves or periodic pulse sequences. 
     For an example, assume that the induced shear wave has two harmonics (ω 2 =2ω 1 ) at location x:
 
 d ( x,t )= D (ω 1 ) e   −h     I     (ω     1     )x  cos(ω 1   t−h   r (ω 1 ) x )+ D (ω 1 ) e   −h     I     (ω     1     )x  cos(ω 1   t−h   r (ω 1 ) x )  (11)
 
Further assume that the tissue is nonlinear and its transfer function is modeled by a square term at distance x+Δx:
 
[ D (ω 1 ) e   −h     I     (ω     1     )(x+Δx)  cos(ω t−h   r (ω 1 )( x+Δx ))+ D (ω 2 ) e   −h     I     (ω     1     )(x−Δx)  cos(ω t−h   r (ω 2 )( x+Δx ))] 2   (12)
 
This quadrature term produces new frequency components of (ω 2 −ω 1 )=ω A  and (ω 2 +ω 1 )=ω B . Thus, the harmonic of ω A  at x+Δx includes phase information of harmonics ω 1  and ω 2  at x, but ω A −ω 1 . Similarly, the harmonic of ω B  at x+Δx introduce a new harmonic component at 3ω 1 . If such harmonic ω B  exists at the location, it will be modified. Thus, the phase difference used to calculate the velocity does not only depend on the tissue viscoelasticity, but also the nonlinearity of the tissue. As different excitation pulses have different harmonic distributions, the impact the nonlinearity is different to different pulses; thus, the measured velocity does not only depend on the tissue viscoelasticity, but also pulse shapes. This tissue nonlinearity can be measured by the method disclosed in this disclosure.
 
     The general concept of this disclosure to measure the tissue shear nonlinearity and shear viscoelasticity is shown in  FIG. 1 . One ultrasound excitation pulse (a tone burst), or a pulse sequence, is applied to a tissue region  20  to generate radiation force p i (t)  10 . The pulse or the pulses in a sequence have a prescribed pulse width. The radiation force induced shear wave in the tissue region  20 . The shear wave y i (t) induced by p i (t) is detected by pulse echo ultrasound. The procedure is repeated for another pulse, till all shear waves  30 , y 1 (t), y 2 (t) y 3 (t), . . . y M (t) induced by p 1 (t), p 2 (t), p 3 (t), . . . , p M (t) (20) are obtained. The measured shear waves  30  are processed for tissue shear property and nonlinearity estimation. 
     The block diagram of ultrasound system to implement this disclosure is shown in  FIG. 2 . The sinusoidal digital signal is synthesized and generated by a direct digital synthesis (DDS)  212  implemented in a FPGA  210 . The frequency of the signal is prescribed and controlled by configuring the DDS  212  using Control Unit  11  and Computer  90 . The digital sinusoidal signal is formatted and buffered in 13 and is converted to an analog signal by a digital to analog converter (DAC)  231 . DAC  231  outputs the analog signal that has a prescribed center frequency to the pulse with control unit  232 . Pulses  221 ,  222 ,  223 ,  224  show the envelopes of the signal at the output of pulse width control unit  232 . They are called excitation pulses. Only one excitation pulse or a sequence of excitations pulse having the same pulse width are generated to induce shear wave in a tissue region. Beamforming circuit  240  distributes the excitation pulse to multiple channels having different time delay for transmission focus. The excitation pulses are amplified by the power amplifiers  250  and applied to ultrasound transducer via T/R switches  261 . Shear wave is induced in a tissue region  264  by ultrasound radiation force generated by the excitation pulses. 
     Using the pulse echo ultrasound method which has been reported in many literatures, the shear wave can be measured. With proper pulse widths, amplifier gains, and timing control, the same transmitting circuits shown  FIG. 2 , which is discussed above, can be used to generate detection pulses. Periodic detection pulses are transmitted to a detection location in tissue region  264 . The echoes of the detection pulses are received by multiple channels of the ultrasound transducer  262 , passed through the T/R switch  261 , amplified by time-gain control amplifiers, and digitized by the ADC  280 . The digitized echoes are processed by beamforming unit  215  in FPGA  210 . The detected RF signal representing the shear wave is received by shear wave detector  217 , which demodulates the RF signal and extracts the motion representing the shear wave. The shear wave detection is documented by many literatures including the prior art of U.S. Pat. No. 7,785,259. The detected shear wave including amplitudes, group velocity, and phase velocity are transmitted to computer  290 . 
     Above transmission and detection processes are repeated for another pulse having another pulse width, till all shear waves induced by pulses having different pulse width are measured for tissue shear property characterization. 
     Above ideal has been experimentally illustrated by investigating shear wave speeds of a swine liver embedded in a gelatin phantom using an ultrasound experiment system. The ultrasound excitation pulses had a center frequency of 4.1 MHz with pulse widths of 10 μs, 20 μs, . . . , 100 μs. The receiving transducer had a center of 6.25 MHz. The group velocities and phase velocities for each excitation pulse were measured. The experiments were repeated ten times for each pulse in two regions of the liver to find the averages and standard deviations (SD) of the group velocities and phase velocities. 
       FIG. 3  shows the phase velocities of five harmonics of shear waves in the liver. The shear waves are induced by the pulses having different pulse widths from 10 μs to 100 μs for 100 Hz and 200 Hz, and 500 Hz. It shows that the phase velocity of the induced shear wave at a given frequency is pulse width dependent, not a constant. The different phase velocities produce different viscoelasticity estimates, as indicated by equation (10). On other hand, the dispersive property of the phase velocities with different pulse widths represents the tissue pathological property.  FIG. 3  also shows that the differences of the phase shifts due to difference of the pulse widths is larger in the lower frequency range, which is the important frequency range for ultrasound induced shear waves. 
       FIG. 4  shows the group velocities of shear waves in the swine liver phantom. It shows that the group velocity of the induced shear wave is not a constant, but varies with the pulse width. The different group velocities lead to different tissue elasticity estimates as shown in  FIG. 5 , which is obviously not a constant for different excitation pulses. This dispersive properties of group velocities and estimated elasticities with different pulse width represent tissue pathological statues. The phase velocities measured with different excitation pulses are also shown in  FIG. 6 , which shows the velocities versus frequencies. In the prior art, the velocity has a single value at a frequency. In this case, the velocities of shear wave induced by the wider pulses have lower values in lower frequencies than those with narrower pulses, the velocities of shear wave induced by the wider pulses have higher values in higher frequencies than those with narrower pulses. 
     According to one aspect of this disclosures is to measure tissue shear properties using different pulses. The shear group velocities, shear phase velocities, and viscoelasticity measured with different excitation pulses are used to characterize the tissue shear property and tissue pathological statues. The distribution of the phase velocities versus different pulse widths is illustrated in  FIG. 3 . The distribution of the group velocities versus different pulse widths is illustrated in  FIG. 4 . The distribution of the elasticity calculated from the measured group velocities versus different pulse widths is illustrated in  FIG. 5 . The phase velocity versus frequencies using different excitation pulses are shown in  FIG. 6 . 
     According to another aspect of this disclosure is to use the averaged value of all measurements with different pulse widths to characterize the tissue pathological statues. 
     According to another aspect of this disclosure is to use the measurements with a selected pulse width to characterize the tissue pathological statues with a labeled pulse width. 
     According to another aspect of this disclosure is to model and measure the nonlinear tissue responses using different excitation pulses. If tissue responses to different excitation pulses are linear, assumed in the prior art of elastography, shear wave velocities are independent to different excitation pulses. However, as illustrated in  FIGS. 3, 4, and 5 , the shear velocities are pulse dependent. This phenomenon can be modeled by nonlinear models. 
     For the Voigt tissue that has a linear response to an applied force, 
                           Y   i     ⁡     (   ω   )       =           P   i     ⁡     (   ω   )           A   0     ⁡     (       μ   1     +     j   ⁢           ⁢     ωμ   2         )         =       H   ⁡     (   ω   )       ⁢         P   i     ⁡     (   ω   )       /     A   0             ,           ⁢     
     ⁢           ⁢   for     ⁢     
     ⁢       i   =   1     ,   2   ,   …   ⁢           ,   M   ,             (   13   )               
where P i (ω) is Fourier transform of the radiation force that is proportional to the excitation pulses, Y i (ω) is Fourier transform of shear wave y i (t), H(ω) is a transfer function, A 0  represents a transfer efficiency of the applied radiation force to the detection location for the shear wave. (13) can be rearranged as:
 
     
       
         
           
             
               
                 
                   
                     
                       
                         
                           P 
                           i 
                         
                         ⁡ 
                         
                           ( 
                           ω 
                           ) 
                         
                       
                       
                         
                           Y 
                           i 
                         
                         ⁡ 
                         
                           ( 
                           ω 
                           ) 
                         
                       
                     
                     = 
                     
                       
                         
                           A 
                           0 
                         
                         ⁡ 
                         
                           ( 
                           
                             
                               μ 
                               1 
                             
                             + 
                             
                               j 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               ω 
                               ⁢ 
                               
                                   
                               
                               ⁢ 
                               
                                 μ 
                                 2 
                               
                             
                           
                           ) 
                         
                       
                       = 
                       
                         
                           
                             R 
                             i 
                           
                           ⁡ 
                           
                             ( 
                             ω 
                             ) 
                           
                         
                         + 
                         
                           
                             jI 
                             i 
                           
                           ⁡ 
                           
                             ( 
                             ω 
                             ) 
                           
                         
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       for 
                       ⁢ 
                       
                         
                             
                         
                         ⁢ 
                         
                             
                         
                       
                       ⁢ 
                       i 
                     
                     = 
                     1 
                   
                   , 
                   2 
                   , 
                   … 
                   ⁢ 
                   
                       
                   
                   , 
                   
                     M 
                     . 
                   
                 
               
               
                 
                   ( 
                   14 
                   ) 
                 
               
             
             
               
                 
                   
                     
                       
                         
                           I 
                           i 
                         
                         ⁡ 
                         
                           ( 
                           ω 
                           ) 
                         
                       
                       
                         
                           R 
                           i 
                         
                         ⁡ 
                         
                           ( 
                           ω 
                           ) 
                         
                       
                     
                     = 
                     
                       
                         ω 
                         ⁢ 
                         
                             
                         
                         ⁢ 
                         
                           μ 
                           2 
                         
                       
                       
                         μ 
                         1 
                       
                     
                   
                   , 
                   
                     
 
                   
                   ⁢ 
                   
                     
                       for 
                       ⁢ 
                       
                         
                             
                         
                         ⁢ 
                         
                             
                         
                       
                       ⁢ 
                       i 
                     
                     = 
                     1 
                   
                   , 
                   2 
                   , 
                   … 
                   ⁢ 
                   
                       
                   
                   , 
                   
                     M 
                     . 
                   
                 
               
               
                 
                   ( 
                   15 
                   ) 
                 
               
             
           
         
       
     
     For the linear tissue, values obtained by (14) and (15) are the same for all different radiation force p i (t) A 0  is unknown, because it dependents on locations in tissue and transmitting transducer and power etc. Thus, the μ 1  and μ 2  in (14) and (15) cannot be calculated with the measurement of y i (t) at a tissue location. The prior art of elastography requires measurements from two locations to estimate the group velocity or phase velocities so that (6) or (10) can be used to calculate μ 1  and μ 2 . 
     For the nonlinear response, the Fourier transforms of tissues responses are: 
                         Y   i     ⁡     (   ω   )       =         1     A   0       ⁡     [       H   ⁡     (   ω   )       ⁢       P   i     ⁡     (   ω   )         ]       α       ,     
     ⁢       for   ⁢             ⁢             ⁢   i     =   1     ,   2   ,   …   ⁢           ,   M   ,           (   16   )               
where α is defined as the shear nonlinear exponent coefficient (SNEC). For the Voigt tissue,
 
                         Y   i     ⁡     (   ω   )       =         P   i   α     ⁡     (   ω   )             A   0     ⁡     (       μ   1     +     j   ⁢           ⁢   ω   ⁢           ⁢     μ   2         )       α         ,     
     ⁢       for   ⁢             ⁢             ⁢   i     =   1     ,   2   ,   …   ⁢           ,   M           (   17   )               
As shown in (17) the linear relationship between A 0  and (μ 1 +jωμ 2 ) is no longer valid. Thus, the four unknown (A 0 , μ 1 , μ 2 , and α) at the right hand side of (17) are independent each other and can be numerically solved by fitting with the Y i (ω) and P i (ω) over a frequency range. Y i (ω) is given by the Fourier transform of the shear wave measurements. Recognizing the excitation pulse is rectangular, the relationship between the radiation force and the excitation pulse is a scaling factor. Thus, P i (ω) can be given by the Fourier transform of the excitation pulse, as the scaling factor is included in A 0 .
 
     For example, a tone burst (M=1) with a certain width is transmitted to induce a shear wave in a tissue region. The shear wave is detected at a location that is near or at the location where the radiation force is applied. The shear wave at this location is detected by pulse echo ultrasound and the Fourier analysis is applied to the detected shear wave. Let ω increases from 1 Hz to 800 Hz with a step size of 1 Hz, there are 2×800=1600 complex data of Y i (ω) and P i (ω) for fitting (17) and solving the four unknown variables: A 0 , μ 1 , μ 2 , and α. In this case, a single transmission of vibration pulse is transmitted to generate the shear wave and viscoelasticity is estimated at a selected location. 
     When the number M of different pulses increases, the numerical fitting is enhanced and the reliability of the estimates is increased. For example, select M=8, so that the pulse widths are 25, 50, 75, . . . , 200 μs, with a step size of 25 μs. Thus, there are 2×800×8=12,800 data available to fit (17) for unknown variable of A 0 , μ 1 , μ 2 , and α. The frequency range and the pulse width in this example are for illustration purpose, which are selected and prescribed for different applications. 
     According to another aspect of the disclosure, the viscoelasticity parameters (μ 1 , μ 2 )) of tissue can be estimated by the measurements of shear wave at one location, which is different than the prior art of elastography. 
     According to another aspect of the disclosure, the measured viscoelasticity of every location in a tissue region can be used to form a two-dimensional image or three-dimensional image, which presents the viscoelasticity of the tissue region. 
     According to another aspect of this disclosure, the accuracy and reliability of the estimates of viscoelasticity (μ 1 , μ 2 )) should be improved as the number of different pulses increases from M=1 to a large number. 
     According to another aspect of this disclosure, tissue is modeled by a general form that includes linear and nonlinear components of tissue:
 
 Y   i (ω)= A   1   {H (ω) P   i (ω)+β[ H (ω) P   i (ω)] α }, for  i= 1,2, . . . , M.   (18)
 
where A 1  is an inverse of A 0 . For the Voigt tissue,
 
 Y   i (ω)= A   1   {P   i (ω)/(μ 1   +jωμ   2 )+β P   i   α (ω)/(μ 1   +jωμ   2 ) α }, for  i= 1,2, . . . , M,   (19)
 
where β is defined as shear nonlinear content coefficient (SNCC) to represent the content of the nonlinearity of the tissue. Because the tissue responses are not completely linearly changed in the frequency range, five parameters A 0 , μ 1 , μ 2 , α, and β are mutually independent. Therefore, the Fourier transforms of the pulses and detected shear wave induced by either one single pulse or multiple pulses or different pulses can be applied to fit (19) for solving unknown A 0 , μ 1 , μ 2 , α, and β.
 
       FIG. 7  shows the key procedure to measure tissue nonlinearity, viscoelasticity, shear motion, and shear velocities by using a single pulse or different excitation pulses at one or more locations. 
     Certain methods and devices disclosed in this disclosure have several important advantages:
         1. The estimation of tissue viscoelasticity requires only the shear wave measured from one location in a tissue region.   2. The estimate of tissue viscoelasticity can be done with a single transmission of an excitation pulse.   3. The multiple transmission of the excitation pulses enhance the reliability of the estimates.   4. The transmission of the multiple excitation pulses enhance the accuracy and reliability of the estimates.   5. Quantitative 2D and 3D images for viscoelasticity of each point in a tissue region issue can be obtained.   6. The nonlinear parameter α and β are calculated to characterize tissue shear wave property and pathology statues.   7. Parameter β indicates the linearity level. When β=0.0, only the ratio between elasticity and viscosity can be obtained and measurements at two locations are required for estimating tissue viscoelasticity.   8. The phase shift distribution of shear waves induced by different excitation pulses is a unique tissue mechanical property, which directly represent tissue shear property including the nonlinearity.       

     As the measurement of shear waves is pulse dependent, one aspect of this disclosure to label viscoelasticity measurements with used pulse width and shapes. 
     While this disclosure is benefited by the prior art of elastography for the transmitting the ultrasound radiation force and detection of the shear wave using pulse-echo ultrasound, the following aspects of this disclosure are unique. 
     One of the aspects of this disclosure is to transmit one excitation pulse to a location in a tissue region to induce shear wave, detect the induced shear wave, take Fourier analysis and fit the results with (17) or (19) to estimate shear elasticity, shear viscosity, shear nonlinear coefficient of the tissue at the location. This allows the estimates of viscoelastic parameters without using shear velocities which require measurements of shear waves at two or more locations. 
     Another aspect of this disclosure is to transmit different excitation pulses having different pulse widths to a location in a tissue region to induce shear wave, detect the induced shear wave, take Fourier analysis and fit the results with (17) or (19) to estimate elasticity, viscosity, nonlinear coefficient α and β of the tissue at the location, shown in  FIG. 8 . The locations of excitation and detections can be the same location as shown in  FIG. 9 . This allows the estimates of viscoelastic parameters without using shear velocities which require measurements of shear waves at two or more locations. 
     Another aspect of this disclosure is to form a quantitative viscoelastic and nonlinear coefficients α and β 2D or 3D image of a tissue region as the viscoelasticity of each location in the region is truly independent measured as shown in  FIG. 10 . 
     Another aspect of the disclosure is to calculate the difference of phase shifts of shear waves induced by two or more different excitation pulses, and use the difference to characterize tissue shear mechanical property. 
     Another aspect of this disclosure is to transmit different excitation pulses having different pulse widths to two or more locations in a tissue region to induce shear wave, detect the induced shear wave, take Fourier analysis and fit the results with (17) or (19) to estimate elasticity, viscosity, nonlinear coefficient of the tissue in the region, as shown in  FIG. 9 . 
     Another aspect of this disclosure is to transmit different excitation pulses having different pulse widths to two or more locations in a tissue region to induce shear wave, detect the induced shear wave, and use the phase shifts distributions and differences with pulse widths to characterize tissue shear property.