CN101530333A - Ultrasonographic system, distortion distribution display method, and elastic modulus distribution display method - Google Patents

Abstract

A measurement point is set at frame data consisting of an envelope signal obtained by orthogonal detection of an RF signal output from an ultrasonic probe before and after compression of an examinee. The invention provides an ultrasonographic system, a sidtortion distribution display method and elastic modulus distribution display method, making the envelope move relatively in order in the ultrasonic beam direction -first dimensional direction, second dimensional direction or third dimensional direction so as to obtain a position where the envelope signal correlation coefficient becomes maximum and the phase difference in the position before and after compression, according to this, the displacement of the measurement point accompanying the compression is obtained and delayed, making one of a pair of signals, obtained by orthogonal detection of a pair of echo signals, move relatively in order in the ultrasonic beam direction -first dimensional direction and set intervals, thus, it is possible to reduce the calculation time, and cope with a horizontal direction displacement.

Description

Compuscan and stress distribution and elastic modelling quantity distribution display packing

The application is that application number is " 03820755.9 (PCT/JP03/09731) ", and the applying date is on July 31st, 2003, and denomination of invention is divided an application for the application of " compuscan and distortion distribution display method ".

Technical field

The present invention relates to compuscan and stress distribution (strain distribution) display packing, this compuscan and distortion distribution display method allow the user to utilize ultrasonic diagnostic equipment, and the hardness of tissue is carried out quantitative measurement.

Background technology

Ultrasonic diagnosis not only is applied to the observation to organizational structure, and is applied to following field (ultrasonic tissue sign): to measuring such as physical quantity in the tissue of the velocity of sound, damped coefficient etc., and according to the physical quantity of such measurement, generate diagnostic image.As the part in this field, known such field, wherein to the hardness of tissue, promptly elastic characteristic is measured.Because the elastic characteristic of tissue and pathological condition have relation closely, therefore above-mentioned field is by primary study.For example, the known tissue that is subjected to following influence shows the hardness bigger than normal structure: the hardening tumor, such as breast carcinoma, thyroid carcinoma etc.; Liver cirrhosis; Or the like.By convention, by touching the hardness that detects tissue.Yet, have the shortcoming that is difficult to carry out objective analysis by the detection that touches, need skilled surgeon, and have following restriction: can only detect and have certain size or large scale and be positioned near the body surface affected tissue more.

On the other hand, known a kind of like this method, wherein apply static pressure to tissue, so that make tissue compression and distortion, and utilize ultrasound wave to measure and organize internal strain (strain) with applied pressure is corresponding, so that estimate the elastic characteristic (J.Ophir of tissue, I.Cespedes, H.Ponnekanati, Y.Yazdi, and X.Li, " Elastography:A quantitative methodfor imaging the elasticity of biological tissue ", Ultrasonic Imaging, Vol.13, pp.111-134,1991).Routine techniques is based on that the following fact is developed: have being organized in of big hardness and show the little strain of organizing under the pressure, and on the other hand, have being organized in of little hardness and show the big strain of organizing under the pressure.Just,, apply static pressure, and, estimate the elastic characteristic of organizing according to the in-house stress distribution under such applied pressure to tissue about above-mentioned conventional method.

Particularly, not by ultrasound probe under the situation that tissue is exerted pressure, utilize the ultrasonic diagnostic equipment have ultrasound probe, carry out the normal measurement of ultrasound echo signal (radio frequency under the situation of not exerting pressure (RF) signal).Subsequently, after this surgeon measures the ultrasound echo signal (the RF signal under the pressure) that passes the tissue that is pressed by ultrasound probe, exert pressure (approximately a few percent) to tissue slightly.Then, utilize the space correlation method, according to the RF signal under the situation of exerting pressure and not exerting pressure to tissue, estimate Displacements Distribution, this Displacements Distribution is represented the displacement of the every bit of the tissue that causes owing to such applied pressure.

The space correlation method has following mechanism: wherein, by the template matching (mask matching) of utilizing two-dimensional correlation function, according to the RF signal under the situation of exerting pressure and not exerting pressure (or envelope signal of RF signal), estimate the in-house Displacements Distribution under applied pressure to tissue.Promptly, the two-dimensional correlation window (template) that will have certain size is applied to and does not have a corresponding RF signal data of tomography data under the situation that tissue is exerted pressure, so that, estimate the displacement of the expectation measurement point on the two-dimensional surface by utilizing auto-correlation processing, detecting at the RF signal data of having used correlation window and the maximal correlation between the RF signal data under the situation that tissue is exerted pressure.For example,, all carry out above-mentioned auto-correlation processing, estimate stress distribution thus each measurement point that the form with grid is provided with.In general, because the sampling of horizontal direction is more coarse than axial sampling, make for the processing that utilizes the space correlation method the axial displacement detecting low precision of displacement detecting ratio of precision of its horizontal direction (scanning direction of ultrasonic beam).As mentioned above, space correlation method has the advantage that makes the user can estimate the two-dimension displacement vector.And, though above-mentioned space correlation method has the shortcoming that the Displacement Estimation precision is subjected to the restriction in sampling interval, but the space correlation method has the following advantages: (for example, about 5% under) the situation, also make the user can estimate Displacements Distribution even at tissue big distortion takes place.Yet different with conventional ultrasonic diagnosis, the space correlation method has following shortcoming: space correlation is handled and is expended a large amount of calculating, has caused being difficult to handle in real time.

Therefore, the purpose of this invention is to provide a kind of method that is used for obtaining in real time Displacements Distribution, stress distribution and elastic modelling quantity distribution (elastic modulus distribution).

Summary of the invention

In order to address the above problem, first compuscan of the present invention comprises: ultrasound probe, and carry out the transmitting-receiving of ultrasonic signal between the object; Correlation calculation apparatus, allow the envelope of a pair of reflection echo signal before and after the compression of the described object that received by described ultrasound probe relatively move successively on the direction of principal axis of ultrasonic beam and set at interval, the correlation coefficient of obtaining the envelope before and after the compression becomes maximum relevant maximum position and the phase contrast under this relevant maximum position; With the strain accountant,, obtain the stress distribution of the tissue of the displacement of tissue of the described object of following described compression and described object according to described relevant maximum position of obtaining by described correlation calculation apparatus and described phase contrast; Described correlation calculation apparatus has: first orthogonal demodulation circuit, carry out orthogonal detection to the reflection echo signal before compressing; Second orthogonal demodulation circuit carries out orthogonal detection to the reflection echo signal after the compression; With a plurality of delay circuits, the output of described second orthogonal demodulation circuit is postponed setting-up time successively; According to the output of described first orthogonal demodulation circuit and described a plurality of delay circuits, obtain described correlation coefficient and described phase contrast.

In order to address the above problem, second compuscan of the present invention, according under the situation of exerting pressure and not exerting pressure, the reflection echo signal (RF signal) that utilizes ultrasound probe measuring object (subsject) to obtain, obtain the displacement of the tissue of object, this compuscan comprises: storage device, be used for that the described ultrasound probe of memory by using detects, such as the characteristics of signals of the envelope signal of signal; Correlation calculation apparatus, be used for described characteristic according to described storage device is stored under the situation of exerting pressure and not exerting pressure to object, the correlation coefficient of calculating between the described characteristic under the situation of exerting pressure and not exerting pressure, and the phase contrast between the received signal under the situation of exerting pressure and not exerting pressure to described object; And the displacement accountant, be used for according to by correlation calculation apparatus like this that obtain, between the RF signal under the situation of exerting pressure and do not exert pressure correlation coefficient and phase contrast, calculate displacement owing to described each measurement point that causes of exerting pressure.In addition, compuscan according to the present invention comprises: the strain accountant, be used for carrying out space differentiation by displacement to each measurement point, and come the stress distribution of the tissue of calculating object; And display device, be used to show the stress distribution of such acquisition.

As mentioned above, about compuscan according to the present invention, according to exert pressure with the situation of not exerting pressure under relevant such as between the characteristic of envelope signal, the displacement of calculating each measurement point can be estimated Displacements Distribution thus in real time.In addition, about compuscan according to the present invention, can be by on described ultrasonic beam direction, change described measurement point with the spacing of 1/2nd wavelength of described ultrasonic signal, obtain between the envelope signal under the situation of exerting pressure to measurement point and do not exert pressure, to show the position of each measurement point of maximum correlation coefficient, solve aliasing (aliasing) problem thus as the shortcoming of Doppler (Doppler) method.

Note, for such correlation computations, wherein by change the envelope signal under the situation of exerting pressure along time shaft, to obtain the auto-correlation function of the envelope signal of each measurement point under the situation of exerting pressure, calculate the position that between the envelope signal under the situation of exerting pressure and not exerting pressure, shows each measurement point of maximum correlation coefficient to measurement point, this correlation computations has caused long computation time, often causes the problem that can not calculate in real time.

Therefore, about compuscan according to the present invention, at first, preferably calculate the auto-correlation function of the envelope signal under the situation of exerting pressure and do not exert pressure.Then,, for example change phase contrast between the auto-correlation function, obtain correlation coefficient by spacing with 1/2nd wavelength of ultrasonic signal by changing the phase contrast between the auto-correlation function that changes corresponding to predetermined measurement point, obtain like this.This has reduced the computation time that is used for displacement calculating, can carry out high speed processing thus.In addition, compuscan according to the present invention may further include: storage device is used to store the envelope signal of the RF signal of orthogonal detection; Correlation calculation apparatus is used to calculate the position that shows each measurement point of maximum correlation coefficient between the envelope signal under the situation of exerting pressure and not exerting pressure to measurement point, and described measurement point is surrounded by the two-dimensional correlation window; And the displacement accountant, be used for calculating the two-dimension displacement at least of each measurement point that causes owing to exerting pressure according to each the measurement point position that obtains like this by correlation calculation apparatus, show correlation coefficient and phase contrast.

Just, proposed to be called as the method for " combination autocorrelation method (CA method) " in this manual.As mentioned above, the combination autocorrelation method has the following advantages: the two and three dimensions displacement measurement in utilizing the space correlation method of correlation window; And the real-time and high precision computation in Doppler (Doppler) method.The combination autocorrelation method allows the user to estimate Displacements Distribution, and irrelevant with the displacement of horizontal direction to a certain extent.In this case, two-dimensional directional can comprise: the ultrasonic beam direction, and ultrasound probe receives ultrasonic signal in this ultrasonic beam direction; And ultrasonic beam scanning direction.In addition, about compuscan according to the present invention, preferably, by on the ultrasonic beam direction with the spacing of 1/2nd wavelength of ultrasonic signal and on the ultrasonic beam scanning direction with the ultrasonic beam spacing, change measurement point, detect the position of each measurement point that shows maximum correlation coefficient.Notice that though according to the present invention, the change spacing of the measurement point on the ultrasonic beam direction is not limited to 1/2nd wavelength of ultrasonic signal, preferably adopts the spacing less than 1/2nd wavelength of ultrasonic signal.

In order further to improve the computational speed that displacement is calculated, at first, preferably calculate the auto-correlation function of the envelope signal under the situation of exerting pressure and not exerting pressure.Then,, change the phase contrast between the auto-correlation function, obtain to show the position of each measurement point of maximum correlation coefficient by in the scope that changes corresponding to predetermined measurement point.

In addition, displacement calculating according to the present invention not only can be applied to two dimension and calculate, and can be applied to three-dimensional computations.Have the three-dimensional configuration of the ultrasound probe of one-dimensional array structure about employing, the frame data of storing in the storage device comprise the volume data of being made up of a plurality of frame data collection, and each of a plurality of frame data collection is as slice of data.On the other hand, have the three-dimensional configuration of the ultrasound probe of two-dimensional array structure about employing, data comprise the envelope signal that obtains by scanning on slice direction.Correlation calculation apparatus is by on three-dimensional, about volume data, the measurement point that change is surrounded by three-dimensional correlation window, detect the position that shows each measurement point of maximum correlation coefficient between the envelope signal under the situation of exerting pressure to measurement point and do not exert pressure, described measurement point is surrounded by three-dimensional correlation window.Simultaneously, correlation calculation apparatus calculates the RF phase difference between signals under the situation of exerting pressure and not exerting pressure.In this case, three-dimensional can comprise: the ultrasonic beam direction, and ultrasound probe receives ultrasonic signal in this ultrasonic beam direction; The ultrasonic beam scanning direction; And with the orthogonal slice direction of above-mentioned both direction.In addition, correlation calculation apparatus preferably calculates, ultrasonic beam direction, ultrasonic beam scanning direction and with the orthogonal slice direction of above-mentioned both direction on RF phase difference between signals under the situation of exerting pressure and not exerting pressure.In addition, above-mentioned high speed processing method can be applied to three-dimensional configuration.In addition, can be unit by slice distance, on slice direction, change measurement point, carry out aforementioned calculation with ultrasonic beam.

In addition, be used to obtain the method that elastic modelling quantity distributes according to the present invention and can comprise the elastic modelling quantity accountant, this elastic modelling quantity accountant is used for: by object being divided into the unit (element) of limited quantity, set up two dimension or three-dimensional FEM (finite element) model at least; And, come the calculating elastic modulus to distribute according to the information and the such stress distribution that obtains that are used to set up model.In addition, can utilize display device to show that the elastic modelling quantity of such acquisition distributes.In this case, the elastic modelling quantity accountant preferably passes through, tissue at object shows under isotropic elasticity and the supposition near Incoercibility, and the tissue of object is divided into the rectangular parallelepiped protrusion part unit of limited quantity, sets up three-dimensional finite element model.In addition, the elastic modelling quantity accountant preferably all shows under the supposition of uniform elastic modelling quantity, the homogeneous state of stress and homogeneous strain in each unit of the tissue of object, utilize elastic equation, according to about the information of above-mentioned stress distribution, come the calculating elastic modulus to distribute.

As mentioned above, show under isotropic elastic supposition, carry out according to calculating of the present invention at tissue.Its reason is, is putting under the outside static pressure of tissue, and the relation between the stress and strain shows linear relationship.Therefore, can carry out approximate calculation at tissue as under the supposition of elastic model.In addition, tissue shows isotropic characteristic, and therefore, can carry out according to calculating of the present invention at tissue as under the supposition of isotropic elastic model.On the other hand, about the present invention, under tissue shows supposition near Incoercibility, calculate.Its reason is, if show at tissue under the supposition of complete Incoercibility, calculate, promptly utilize 0.5 the interior constant Poisson of tissue to calculate than (Poisson's ratio), then elastic equation becomes special circumstances, has caused the problem that can not utilize Finite Element Method to calculate.Note,, utilize the interior Poisson's ratio of uniform tissue to calculate about the present invention.In this case, only should estimate Young modulus (Young's modulus), distribute, simplify inversion problem thus to estimate elastic modelling quantity.Notice that Poisson's ratio shows uniformity in enough tissues, and therefore,, preferably utilize 0.49 Poisson's ratio to calculate about the present invention.Elastic modelling quantity distribution according to the present invention is calculated and can only be distributed according to the axial strain that can calculate accurately, come the reconstruct elastic modelling quantity to distribute, thus calculating elastic modulus distribution stably.

Distortion distribution display method of the present invention comprises: first step, and input is by contacting obtained ultrasound echo signal with ultrasound probe with the tissue of object; Second step, allow the envelope of a pair of reflection echo signal before and after the compression of the described object that received by described ultrasound probe relatively move successively on the direction of principal axis of ultrasonic beam and set at interval, the correlation coefficient of obtaining the envelope before and after the compression becomes maximum relevant maximum position and the phase contrast under this relevant maximum position; Third step according to described relevant maximum position of obtaining in described second step and described phase contrast, is obtained the stress distribution of the tissue of the displacement of tissue of the described object of following described compression and described object; With the 4th step, show described stress distribution.Described second step has: the first orthogonal detection step, the reflection echo signal before compressing is carried out orthogonal detection; The second orthogonal detection step is carried out orthogonal detection to the reflection echo signal after the compression; With postpone step, the output of the described second orthogonal detection step is postponed setting-up time successively; According to the output of described first orthogonal detection step and described delay step, obtain described correlation coefficient and described phase contrast.

Elastic modelling quantity distribution display packing of the present invention comprises: first step, and input is by contacting obtained ultrasound echo signal with ultrasound probe with the tissue of object; Second step, allow the envelope of a pair of reflection echo signal before and after the compression of the described object that received by described ultrasound probe relatively move successively on the direction of principal axis of ultrasonic beam and set at interval, the correlation coefficient of obtaining the envelope before and after the compression becomes maximum relevant maximum position and the phase contrast under this relevant maximum position; Third step according to described relevant maximum position of obtaining in described second step and described phase contrast, is obtained the stress distribution of the tissue of the displacement of tissue of the described object of following described compression and described object; The 4th step is divided into the unit of limited quantity with the tissue of described object, sets up two-dimensional finite element model at least, according to the information and the described stress distribution that are used to set up described model, comes the calculating elastic modulus to distribute; With the 5th step, show that described elastic modelling quantity distributes.Described second step has: the first orthogonal detection step, the reflection echo signal before compressing is carried out orthogonal detection; The second orthogonal detection step is carried out orthogonal detection to the reflection echo signal after the compression; With postpone step, the output of the described second orthogonal detection step is postponed setting-up time successively; According to the output of described first orthogonal detection step and described delay step, obtain described correlation coefficient and described phase contrast.

Description of drawings

Fig. 1 has described the mechanism of ultrasonic diagnostic equipment.

Fig. 2 has shown the specific examples by the tissue elasticity measurement method that applies static pressure, and the mechanism of the tissue elasticity measurement method by applying static pressure.

Fig. 3 has shown the mechanism of space correlation method.

Fig. 4 has shown the mechanism of Doppler (Doppler) method.

Fig. 5 has shown the mechanism of combination autocorrelation method.

Fig. 6 has shown the block diagram of the circuit configuration of the rudimentary algorithm that is used to carry out the combination autocorrelation method.

Fig. 7 has shown the block diagram according to the example arrangement of the compuscan of the embodiment of the invention.

Fig. 8 has shown the flow chart of the rudimentary algorithm of three-dimensional arrangement autocorrelation method.

Figure 9 shows that, according to the flow chart of the rudimentary algorithm of the three-dimensional arrangement autocorrelation method that is adopted in the compuscan of the present invention, and the detail flowchart of the part of processing shown in Figure 7.

Figure 10 has shown, is used to describe in detail the flow chart of the combination autocorrelation method of the computational speed with raising, and this combination autocorrelation method is performed in step S15 shown in Figure 9.

Figure 11 shown, is used for carrying out the block diagram of circuit configuration of the rudimentary algorithm of the three-dimensional arrangement autocorrelation method that adopts according to compuscan of the present invention.

Figure 12 has shown the sketch map of simulation process.

Figure 13 shown, the example of the point spread function of every bit under the situation of the ultrasound wave mid frequency of 5.0MHz.

Figure 14 has shown the sketch map of organize models.

Figure 15 shown, utilize that each method of estimation is estimated, because the estimation strain error that the horizontal direction displacement causes.

Figure 16 shown, the stress distribution of utilizing each method of estimation (combination autocorrelation method, expansion combination autocorrelation method and space correlation method) to estimate under the situation of the horizontal direction displacement of 0.0mm.

Figure 17 shown, the stress distribution of utilizing each method of estimation (combination autocorrelation method, expansion combination autocorrelation method and space correlation method) to estimate under the situation of the horizontal direction displacement of 0.4mm.

Figure 18 shown, is used to simulate the organize models of the compression of incline direction.

Figure 19 has shown, by the simple compression of axial organize models being simulated the estimated result of the stress distribution that obtains.

Figure 20 has shown, by the estimated result that the stress distribution that obtains is simulated in organize models's compression of incline direction.

Figure 21 has shown two kinds of organize models's examples, and each of these two kinds of organize models's examples all has three dimensional structure.

Figure 22 has shown first figure that utilizes inclusions to comprise the model results estimated.

Figure 23 has shown second figure that utilizes inclusions to comprise the model results estimated.

Figure 24 has shown first figure that utilizes layer structure model results estimated.

Figure 25 has shown second figure that utilizes layer structure model results estimated.

Figure 26 has shown the basic configuration of compuscan.

Figure 27 has shown, the customized configuration of the ultrasound scanner that adopts in the compuscan.

The specific embodiment

Below with reference to accompanying drawing, the compuscan according to the embodiment of the invention is described.Adopt a kind of being called as the method for " expansion combination autocorrelation method " according to compuscan of the present invention, the one dimension detection of wherein utilize the combination autocorrelation method, carrying out correlation computations by the one dimension window, according to envelope signal handle be expanded into, by the two-dimensional correlation window, utilize two-dimentional detection, come the displacement of processing horizontal direction.In addition, utilization is according to the expansion of embodiment combination autocorrelation method, only to axial spacing be the ultrasound wave wavelength half, horizontal-direction pitch is the mesh point of bunch spacing (beam-line pitch), carries out the envelope correlation computations, so that reduce amount of calculation, make it possible to carry out supercomputing thus.Note, make up autocorrelation method as the combination autocorrelation method, adopt phase information to improve the estimated accuracy of axial displacement according to the expansion of present embodiment.Yet,, thereby do not utilize phase information to estimate the displacement of horizontal direction owing to the signal that lacks as carrier (carrer).Therefore, as the space correlation method, the estimated accuracy of horizontal direction displacement is subject to the sampling interval (ultrasonic beam distance between centers of tracks).Note, do not utilize the improve the standard estimated accuracy of direction displacement of special mechanism according to the expansion of present embodiment combination autocorrelation method, because can utilize the elastic modelling quantity distribution reconstructing method described later, only distribute, estimate that elastic modelling quantity distributes according to axial strain (displacement).

Before the customized configuration of description, will the combination autocorrelation method that make up the basis of autocorrelation method as expansion be described referring to figs. 1 to Fig. 6 according to the expansion combination autocorrelation method of present embodiment.Fig. 1 has described the principle of ultrasonic diagnostic equipment.Be can clearly be seen that to have as the ultrasound probe of ultrasonic detector and convert the electrical signal to ultrasound wave and convert ultrasound wave the function of the signal of telecommunication to by Fig. 1, this allows user to tissue 11 projection ultrasonic pulse signals.A part of passing the ultrasonic pulse signal of tissue 11 is reflected on first interface 12 that has between the zone of mutually different acoustic impedance.The reflection supersonic wave that will be called as " reflection echo signal 12a " passes tissue 11, invests ultrasound probe 10, and remaining ultrasound wave passes first interface 12.A part of passing the ultrasonic pulse signal at first interface 12 is reflected on the second contact surface 13 having between the zone of mutually different acoustic impedance.Similarly, reflectance ultrasound wave pulse signal reflection, that will be called as " reflection echo signal 13a " passes tissue and invests ultrasound probe 10 on second contact surface 13, and on the other hand, remaining ultrasonic pulse signal passes second contact surface 13.The reflection supersonic wave echo-signal is received by ultrasound probe 10, so that be converted into reflection echo signal, this reflection echo signal is the signal of telecommunication.In this case, represent with following formula (1), from ultrasound probe 10 projection ultrasonic pulse signals, up to receive from and the echo-signal of the reverberation 14 of ultrasound probe 10 distance L (having the interface between the zone of mutually different acoustic impedance) reflection till period t.

$t=\frac{2L}{c}\·\·\·\left(1\right)$

At this, c represents the in-house velocity of sound, and when passing soft tissue, c can be confirmed as the constant of about 1500 meter per seconds.Therefore, according to from throwing ultrasound wave, calculate the distance L between probe and the reverberation up to the time t that receives reflection echo signal.In addition, reflection echo signal has the information of acoustic characteristic about tissue, and therefore can be according to reflection echo signal, will be such as the organizational information pictorial display of B pattern tomographic map on monitor.

For example, known a kind of like this method wherein utilizes ultrasonic diagnostic equipment to measure the elastic characteristic that hardness is organized in representative.The above-mentioned method that is used to measure elastic characteristic has a kind of mechanism, wherein mechanical vibration are acted on tissue, and pass sclerous tissues and pass this fact of soft tissue with high spread speed with low propagation speed based on shear wave (transverse wave), according to the spread speed of the shear wave of such generation, estimate hardness information.Say that strictly the spread speed v that propagates the shear wave that passes tissue depends on density, the modulus of shearing μ that organizes ρ ₁, shear viscosity μ ₂With the angular frequency of vibration, represented as shown in the formula (2).

$v=\sqrt{\frac{2({{\mathrm{\μ}}_1}^2+{\mathrm{\ω}}^2{{\mathrm{\μ}}_2}^2)}{\mathrm{\ρ}({\mathrm{\μ}}_1+\sqrt{{{\mathrm{\μ}}_1}^2+{\mathrm{\ω}}^2{{\mathrm{\μ}}_2}^2})}}\·\·\·\left(2\right)$

On the other hand, as shown in Figure 2, proposed a kind of like this method, wherein static pressure has been put on tissue, and, estimate the elastic characteristic of organizing according to the stress distribution of organizing under the static pressure that is applied.Said method is based on the following fact and designs: this static pressure causes small strain in sclerous tissues, and causes big strain in soft tissue.Fig. 2 (A) has shown the specific examples by the tissue elasticity measurement method that applies static pressure.Fig. 2 (B) has shown the mechanism by the tissue elasticity measurement method that applies static pressure.Can clearly be seen that by figure, about said method, utilize have the conventional ultrasonic diagnostic equipment of ultrasound probe 10, under the situation of not exerting pressure, normally measure the ultrasound echo signal (the RF signal under the situation of not exerting pressure) of tissue 11.Subsequently, by ultrasound probe 10 slightly (approximately a few percent) push tissue 11, and the ultrasound echo signal (the RF signal under the situation of exerting pressure) of the tissue 11 of measurement under exerting pressure.Then, according to the RF signal of under the situation of exerting pressure and not exerting pressure, measuring, the Displacements Distribution of the displacement of the in-house every bit when estimating to represent pressurized to tissue.The main example of Displacements Distribution method of estimation comprises the method for utilizing space correlation and utilizes the method for Doppler (Doppler) effect.

Fig. 3 has shown the mechanism of space correlation method.About this method, carry out template matching by utilizing two-dimensional correlation function, according to the RF signal of under the situation of exerting pressure and not exerting pressure, measuring (or envelope signal of RF signal), estimate the in-house stress distribution under exerting pressure to tissue.The particular procedure of this method below will be described.At first, respectively with i ₁(t, x) and i ₂(t x) is illustrated in the RF signal of measuring (or envelope signal of RF signal), cross correlation coefficient C (t, x between then above-mentioned two signals under the situation that tissue is exerted pressure and do not exerted pressure; N m) represents with following formula (3).

$C(t,x;n,m)=\frac{\underset{v={-\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{\overset{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{\mathrm{\&Sigma;}}}\underset{w={-x}_0/2}{\overset{x_0/2}{\mathrm{\&Sigma;}}}{i}_1(t+v,x+w)i_2(t+v+{\mathrm{nL}}_t,x+w+{\mathrm{mL}}_x)}{\sqrt{\underset{v={-\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{\overset{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{\mathrm{\&Sigma;}}}\underset{w={-x}_0/2}{\overset{x_0/2}{\mathrm{\&Sigma;}}}{i_1}^2(t,x)}\&CenterDot;\sqrt{\underset{v={-\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{\overset{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{\mathrm{\&Sigma;}}}\underset{w={-x}_0/2}{\overset{x_0/2}{\mathrm{\&Sigma;}}}{i_2}^2(t+v+{\mathrm{nL}}_t,x+w+{\mathrm{mL}}_x)}}\&CenterDot;\&CenterDot;\&CenterDot;\left(3\right)$

At this, t represents the ultrasonic beam direction coordinate points of (axially), the coordinate points of x representative and ultrasonic beam direction quadrature (horizontal direction), t ₀Represent the axial dimension of correlation window, x ₀Represent the horizontal direction size of correlation window, L _tRepresent the axial sampling interval, and L _xRepresent the sampling interval of horizontal direction.In addition, n and m are integer.In this case, if will demonstrate peaked cross correlation coefficient combination (n, m) be expressed as (k, l), then at measurement point (t, the axial displacement u that x) locates _yWith the horizontal direction displacement components u _xRepresent with following formula respectively.

u _y＝kL _t

u _x＝lL _x

Note sampling interval L in the horizontal direction _xThan axial sampling interval L _tCarry out data sampling under the coarse situation, caused the displacement detecting ratio of precision axial displacement accuracy of detection of horizontal direction poor thus.Each measurement point is carried out above-mentioned processing, estimate Displacements Distribution thus.The advantage of space correlation method is, makes the user can estimate the two-dimension displacement component of a vector.In addition, even occur big strain (about 5%) in the tissue, the space correlation method also allows the user to estimate Displacements Distribution.Yet, being different from conventional ultrasound measurement system, said method causes a large amount of calculating, causes the difficulty of real-time calculating.In addition, the Displacement Estimation precision is subjected to the restriction in sampling interval, has caused thus and the problem that Doppler (Doppler) method is compared, precision is relatively poor, will be described later Doppler (Doppler) method.

Fig. 4 has shown the mechanism of Doppler (Doppler) method.About this method, according to the RF signal under the situation of exerting pressure and not exerting pressure to tissue, utilize Doppler (Doppler) effect to estimate in-house Displacements Distribution under exerting pressure, Doppler (Doppler) effect also is used for blood flow measurement.The particular procedure of this method below will be described.At first, utilize, be illustrated in the RF signal under the situation that tissue is exerted pressure and do not exerted pressure as shown in the formula (4) represented model.

At this, i ₁(t) the RF signal of representative under the situation of not exerting pressure, i ₂(t) the RF signal of representative under the situation of exerting pressure, A (t) represents envelope signal, ω ₀Represent hyperacoustic central angle frequency, and τ represents time shift (time shift).In case each the execution orthogonal detection to two RF signals has just obtained baseband signal, is shown below.

s ₁(t)＝A(t)e ^jθ

$s_2\left(t\right)=A(t-\mathrm{\&tau;})e^{j({\mathrm{\&omega;}}_0\mathrm{\&tau;}+\mathrm{\&theta;})}\&CenterDot;\&CenterDot;\&CenterDot;\left(5\right)$

Equally, represent compound correlative function R between above-mentioned two signals with following formula ₁₂(t).

$R_{12}\left(t\right)={\&Integral;}_{{-\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}^{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{s}_1(t+v)s_2{(t+v)}^{*}\mathrm{dv}=R_A\left(t\right)e^{-j{\mathrm{\&omega;}}_0\mathrm{\&tau;}}\&CenterDot;\&CenterDot;\&CenterDot;\left(6\right)$

At this, R _A(t) represent the auto-correlation function of envelope signal, and t ₀Represent the axial correlation window size of ultrasonic beam.In addition, complex conjugation operator represented in asterisk " * ".Therefore, from correlation function R ₁₂(t) phase (t) has obtained the time shift τ and the axial displacement u that cause owing to exerting pressure _y, shown in (7).

$\mathrm{\&tau;}=-\frac{\mathrm{\&phi;}\left(t\right)}{{\mathrm{\&omega;}}_0}$

$u_y=\frac{\mathrm{c\&tau;}}{2}\&CenterDot;\&CenterDot;\&CenterDot;\left(7\right)$

Notice that c represents the in-house velocity of sound, and supposition c in tissue is a constant.

Utilize Doppler (Doppler) method, each measurement point carried out above-mentioned processing, and by with the same mode of blood flow measurement according to Doppler (Doppler) effect exploitation, estimate Displacements Distribution.Therefore, Doppler (Doppler) method has the advantage of real-time measurement.In addition, Doppler (Doppler) method utilizes phase information to calculate, and compares with the space correlation method thus, can estimate displacement more accurately.Yet the shortcoming of Doppler (Doppler) method is: the measurement of big displacement, and for example 1/4th of the wavelength of ultrasound wave mid frequency or the measurement of bigger displacement, caused aliasing (aliasing), caused estimating the problem of correct displacement.In addition, the shortcoming of Doppler (Doppler) method also is, can not estimate the displacement except that axially as inferring from following formula.

In order to address the above problem, the inventor has proposed " combination autocorrelation method (CA method) ", and this method has the advantage of above-mentioned two kinds of methods simultaneously.Fig. 5 has shown the mechanism of the combination autocorrelation method that is proposed by the inventor.About the combination autocorrelation method, utilize the relevant of envelope signal of RF signal to calculate, solved aliasing problem thus as the shortcoming of Doppler method.The particular procedure of this method below will be described.

At first,, utilize, be illustrated in the RF signal under the situation that tissue is exerted pressure and do not exerted pressure as shown in the formula represented model as the situation of Doppler method.

At this, i ₁(t) the RF signal of representative under the situation of not exerting pressure, i ₂(t) the RF signal of representative under the situation of exerting pressure, A (t) represents envelope signal, ω ₀Represent hyperacoustic central angle frequency, and τ represents time shift.In case each the execution orthogonal detection to two RF signals has just obtained baseband signal, is shown below.

s ₁(t)＝A(t)e ^jθ

$s_2\left(t\right)=A(t-\mathrm{\&tau;})e^{j({\mathrm{\&omega;}}_0\mathrm{\&tau;}+\mathrm{\&theta;})}\&CenterDot;\&CenterDot;\&CenterDot;\left(9\right)$

Then, utilize following formula to represent compound correlative function R between above-mentioned two signals ₁₂(t; N).

$R_{12}(t;n)={\&Integral;}_{{-\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}^{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{s}_1(t+v)s_2{(t+\mathrm{<figure-callout\; id="2"\; label="n\; T"\; filenames="A200910006835E00521.png,A200910006835E00531.png"\; state="\{\{state\}\}">n</figure-callout>}\frac{T}{}\frac{}{2}+v)}^{*}\mathrm{dv}=R_A(t;\mathrm{\&tau;}-n\frac{T}{2})e^{-{\mathrm{j\&omega;}}_0(\mathrm{\&tau;}-n\frac{T}{2})}\&CenterDot;\&CenterDot;\&CenterDot;\left(10\right)$

(n＝...，-2，-1，0，1，2...)

At this, T represents hyperacoustic cycle, R _A(t; τ) represent the auto-correlation function of envelope signal, and t ₀Represent the correlation window size.In addition, complex conjugation operator represented in asterisk " * ".At this, n represents variable number, and different n is calculated at every turn.Around the point number definite by variable, the displacement of calculating each measurement point.Notice that under the situation of n=0, the auto-correlation function in combination auto-correlation function and Doppler as the formula (6) (Doppler) method equates.Just, under the situation that 1/4th or bigger displacement of ultrasound wave wavelength are measured, n equals 0 combination auto-correlation function and has caused the aliasing problem.In order to address the above problem, about the combination autocorrelation method, with envelope correlation coefficient C (t; N) be defined as shown in the formula (11) represented.

$C(t;n)=\frac{\left|R_{12}(t;n)\right|}{\sqrt{|R_{11}(t;0)|\&CenterDot;|R_{22}(t;n)|}}\&CenterDot;\&CenterDot;\&CenterDot;\left(11\right)$

Note R ₁₁(t; 0) represents s ₁(t) auto-correlation function, and R ₂₂(t; N) represent s ₂(t+nT/2) auto-correlation function.If will demonstrate peaked envelope correlation coefficient C (t; N) n is expressed as k, then utilizes the R of n=k ₁₂(t; K) phase is calculated.In this case, no aliasing ground displacement calculating.Its reason is, carries out the envelope correlation computations under the spacing of 1/2nd wavelength.Notice that the calculating spacing of 1/2nd wavelength is to be used for displacement calculating, to prevent the maximum spacing of aliasing simultaneously.Thereby, utilize φ (t; K) calculate time shift T and axial displacement u _y, be shown below.

$\mathrm{\&tau;}=-\frac{\mathrm{\&phi;}(t;k)}{{\mathrm{\&omega;}}_0}+\mathrm{<figure-callout\; id="2"\; label="k\; T"\; filenames="A200910006835E00521.png,A200910006835E00531.png"\; state="\{\{state\}\}">k</figure-callout>}\frac{T}{}\frac{}{2}$

$u_y=\frac{\mathrm{c\&tau;}}{2}\&CenterDot;\&CenterDot;\&CenterDot;\left(12\right)$

Notice that c represents the in-house velocity of sound, and supposition c in tissue is a constant.

Utilize the combination autocorrelation method, each measurement point is all carried out above-mentioned processing, estimate Displacements Distribution thus, it is the extended method of Doppler (Doppler) method.Therefore, the combination autocorrelation method has the advantage of the measurement implemented.In addition, different with Doppler (Doppler) method, the combination autocorrelation method has the following advantages: make the user can utilize envelope to be correlated with to estimate the Displacements Distribution that comprises big displacement (that is, 1/4th of the ultrasound wave wavelength or bigger displacement).

Fig. 6 has shown the block diagram of the circuit configuration of the rudimentary algorithm that is used to carry out the combination autocorrelation method.In Fig. 6, the echo-signal x that obtains under the situation of not exerting pressure (t) is imported into not pressurized orthogonal demodulation circuit (QD) 131, be used for orthogonal detection, and the orthogonal detection signal Ix (t) that detects like this and Qx (t) are imported into the first correlation computations circuit 133 and the first correlation coefficient counting circuit 1350 to 135N.On the other hand, the echo-signal y that obtains under the situation of exerting pressure (t) is imported into the first pressurized orthogonal demodulation circuit (QD) 1320, be used for orthogonal detection, and orthogonal detection signal Y (the t)=Iy+jQy of Jian Ceing (Iy (t), Qy (t)) is imported into the first correlation coefficient counting circuit 1340 and the second correlation computations circuit 1350 like this.First delay circuit 134 makes echo-signal y (t) postpone hyperacoustic period T, and the echo-signal y that postpones ₁=y (t-T) is imported into second pressurized orthogonal detection (QD) circuit 1321.Equally, second delay circuit 135 makes the echo-signal y that has been postponed by first delay circuit 134 ₁=y (t-T) postpones hyperacoustic period T, and the echo-signal y that postpones ₂=y (t-2T) is imported into next second pressurized orthogonal detection (QD) circuit, 1322 (not shown)s.Notice that this circuit has N delay circuit, echo-signal is postponed continuously, and the echo-signal that has been delayed the integral multiple period T is imported into corresponding pressurized orthogonal demodulation circuit in the same way.

The first correlation computations circuit 133 calculates correlation Rxx according to signal Ix and Qx, and correlation Rxx is outputed to the second correlation coefficient counting circuit 1380 to 138N.The second correlation computations circuit 1340 receives orthogonal detection signal Iy (t) and Qy (t) from pressurized orthogonal demodulation circuit 1320, according to signal I _yCalculate correlation Ryy with Qy, and correlation Ryy is outputed to the second correlation computations circuit 1380.The first correlation coefficient counting circuit 1350 never pressurized orthogonal demodulation circuit 131 receives orthogonal detection signal Ix (t) and Qx (t), and receives orthogonal detection signal Iy (t) and Qy (t) from the first pressurized orthogonal demodulation circuit 1320, calculates complex baseband signal S in view of the above _RAnd S _I, and with baseband signal S _RAnd S _IOutput to third phase and close counting circuit 1360 and phase difference calculating circuit 1370.Third phase closes counting circuit 1360 and receives complex baseband signal S from the first correlation coefficient counting circuit 1350 _RAnd S _I, calculate correlation in view of the above | Rxy|, and with the correlation that calculates | Rxy| outputs to the second correlation coefficient counting circuit 1380.Phase difference calculating circuit 1370 receives complex baseband signal S from the first correlation coefficient counting circuit 1350 _RAnd S _I, and according to complex baseband signal S _RAnd S _ICalculate phase difference ₀(t).The second correlation coefficient counting circuit 1380 receives correlation Rxx from the first correlation computations circuit 133, close counting circuit 1360 from third phase and receive correlation | Rxy|, and from the second correlation computations circuit, 1340 reception correlation Ryy, according to above-mentioned correlation value calculation correlation coefficient C ₀And the correlation coefficient C that calculates of output (t), ₀(t).

The second pressurized orthogonal demodulation circuit (QD) 1321 receives the echo-signal y that is postponed by first delay circuit 134 ₁=y (t-T), and with orthogonal detection signal Y ₁(t)=Iy ₁+ jQy ₁(Iy ₁(t), Qy ₁(t)) output to the first correlation coefficient counting circuit 1341 and the second correlation computations circuit 1351.The second correlation computations circuit 1341 receives orthogonal detection signal Iy from the second pressurized orthogonal demodulation circuit (QD) 1321 ₁(t) and Qy ₁(t), according to signal Iy ₁(t) and Qy ₁(t) calculate correlation Ry ₁y ₁, and with correlation Ry ₁y ₁Output to the second correlation computations circuit 1381.The first correlation coefficient counting circuit 1351 never pressurized orthogonal demodulation circuit 131 receives orthogonal detection signal Ix (t) and Qx (t), and receives orthogonal detection signal Iy from the second pressurized orthogonal demodulation circuit (QD) 1321 ₁(t) and Qy ₁(t), calculate complex baseband signal S in view of the above _R1And S _I1, and with complex baseband signal S _R1And S _I1Output to third phase and close counting circuit 1351 and phase difference calculating circuit 1371.Third phase closes counting circuit 1361 and receives complex baseband signal S from the first correlation coefficient counting circuit 1351 _R1And S _I1, according to complex baseband signal S _R1And S _I1Calculate correlation | Rxy ₁|, and with the correlation that calculates | Rxy ₁| output to the second correlation coefficient counting circuit 1381.Phase difference calculating circuit 1371 receives complex baseband signal S from the first correlation coefficient counting circuit 1351 _R1And S _I1, and according to complex baseband signal S _R1And S _I1Calculate phase difference ₁(t).The second correlation coefficient counting circuit 1381 receives correlation Rxx from the first correlation computations circuit 133, closes counting circuit 1361 from third phase and receives correlation | Rxy ₁|, and from the second correlation computations circuit, 1341 reception correlation Ry ₁y ₁, according to above-mentioned correlation value calculation correlation coefficient C ₁And the correlation coefficient C that calculates of output (t), ₁(t).

Similarly, from from each of the second pressurized orthogonal demodulation circuit (QD) 1322 to 132N of the downward phase delay circuit received signal of first delay circuit 135, each of the second correlation computations circuit 1342 to 134N, each of the first correlation coefficient counting circuit 1352 to 135N, third phase closes each of circuit 1362 to 136N, each of phase difference calculating circuit 1372 to 137N, and second each of correlation coefficient counting circuit 1382 to 138N, all carry out processing, export correlation coefficient C thus as the aforesaid first order and second level component ₂(t) to CN (t) and phase value φ ₂(t) to φ _N(t).As mentioned above, the circuit that is used to carry out the rudimentary algorithm of combination autocorrelation method shown in Figure 6 has following configuration: wherein for each of delay circuit 134 to 13N, all make the echo-signal y (t) under exerting pressure postpone the period of hyperacoustic two/one-period T/2 (1/2nd wavelength), and corresponding orthogonal demodulation circuit (QD) 1320 to 132N carry out orthogonal detection (quadrature-detection) to each echo-signal of such delay.

As mentioned above, the space differentiation by the estimation Displacements Distribution under the situation of exerting pressure to tissue has obtained stress distribution.The relative resilient characteristic of stress distribution representative tissue, and therefore shows and the similar effect of diagnosis based on the elastic modelling quantity distribution based on the diagnosis of stress distribution.Yet, under the situation of the hardened liver cirrhosis that causes whole affected tissue, being difficult to carry out the organizational diagnosis that distributes as elastic characteristic, the permission surgeon that distributes quantitatively estimates elastic characteristic.Therefore, in the last few years, the reconstructing method that is used for the distribution of tissue elasticity modulus was studied.Notice that all these reconstructing methods all are in conceptual phase, and do not set up standard method.

On the other hand, can obtain the tissue elasticity modulus and distribute based on aforesaid, in-house stress distribution and stress distribution.Yet, be difficult to utilize prior art directly to measure stress distribution.Therefore, especially distribute,, promptly need to solve inversion problem to satisfy the boundary condition under the situation of exerting pressure to tissue based on the stress distribution elastic modelling quantity of resetting.In general, be difficult to solve inversion problem, and only proposed elastic modelling quantity reconstructing method seldom.Conventional elastic modelling quantity reconstructing method below will be described.

At first, known a kind of like this method, this method are based on that the supposition representing to organize by one-dimensional model (one dimension elastic model) proposes.Just, known a kind of like this method wherein based on the supposition of representing by the one dimension elastic model to organize, comes the calculating elastic modulus to distribute.Based on above-mentioned supposition, elastic modelling quantity is defined as strained anti-number (inverse number).Say that strictly said method is not the elastic modelling quantity reconstructing method.About said method, only obtained strained inverse, and therefore, as the situation of stress distribution, the relative resilient characteristic that can only obtain to organize.

Secondly, known a kind of like this method, wherein elastic equation is reduced to does not have the stress item (the supposition tissue demonstrates isotropic elasticity, Incoercibility and plane strain).About said method, formed, be reduced to the elastic equation of representing plane stress state and do not have the stress item, and utilize the simplification elastic equation do not have the stress item, come reconstruct tissue elasticity modulus to distribute based on stress distribution (institute of strain tensor that comprises the shear strain component is important), to satisfy boundary condition (distribution of exerting pressure on the body surface, the perhaps Displacements Distribution on the body surface).Notice that this method needs such zone (reference area), has obtained elastic modelling quantity in this zone.

The 3rd, known a kind of like this method wherein combines the elasticity differential equation (the supposition tissue demonstrates isotropic elasticity, Incoercibility and plane strain).About said method, formed, be reduced to the elastic equation of representing plane stress state and do not have the stress item, and by continuously combine about near elastic modelling quantity unstressed the simplification differential equation of elastic modelling quantity and as a reference the body surface, based on stress distribution (comprise the shear strain component strain tensor important), come the distribution of reconstruct tissue elasticity modulus.Notice that this method needs such body surface near zone (reference area), has obtained the elastic modelling quantity distribution in advance at this body surface near zone.In addition, said method has following problem: owing to as a reference body surface near elastic modelling quantity combine and caused the accumulation of error, make the table that exsomatizes far away more, the error of calculation is just big more.

The 4th, known a kind of method (the supposition tissue demonstrates isotropic elasticity, Incoercibility and plane strain) of utilizing method of perturbation.About said method, utilize the method for perturbation that has iterative method, exerting pressure based on the body surface of ultrasonic beam direction (axially) distributes and the body surface stress distribution, is formed to represent the elastic equation of plane stress state by separating, and comes reconstruct tissue elasticity modulus to distribute.

The particular problem of having described basic mechanism and will having solved.Below will describe, so that address the above problem according to embodiments of the invention.Fig. 7 has shown the concise and to the point configuration block diagram according to the compuscan of the embodiment of the invention.Ultrasound probe 91 comprises conventional sector scanning probe

(fan-shaped phase modulation array probe), linear scanning probe (linear array probe) or convex scanheads (convex array probe) have to the tissue projection ultrasound wave that will observe and the hyperacoustic function that receives reflection.

The RF signal that obtains under the situation of exerting pressure and not exerting pressure to tissue outputs to quadrature detector 92 from ultrasound probe 91.Quadrature detector 92 is the RF conversion of signals under the situation of exerting pressure and not exerting pressure to tissue the complex envelope signal (IQ signal) of exerting pressure and not exerting pressure, and the IQ signal is outputed to two-dimentional multiple correlation computing unit 93.Two dimension multiple correlation computing unit 93 calculates at the RF signal under the situation of exerting pressure to tissue and less than the two-dimensional correlation between the RF signal under the situation that tissue is exerted pressure, the position that demonstrates maximal correlation is outputed to horizontal direction displacement computing unit 94 and axial displacement computing unit 95, and corresponding correlation function phase place is outputed to axial computing unit 95.Note, axially, calculate relevantly with the spacing of 1/2nd ultrasound wave mid frequencyes, this spacing is to obtain phase place, prevent the maximum spacing of aliasing simultaneously.Calculate with this spacing relevant so that allow the real-time demonstration of compuscan.Therefore, the invention is not restricted to calculate with the spacing of 1/2nd wavelength, but, can construct such configuration, wherein calculate relevant accurately.

Horizontal direction displacement computing unit 94 bases and the corresponding position of horizontal direction maximal correlation of receiving, calculated level direction displacement components u from two-dimentional multiple correlation computing unit 93 _x, and this displacement outputed to horizontal direction strain computing unit 96.On the other hand, axial displacement computing unit 95 bases and axial maximal correlation of receiving from two-dimentional multiple correlation computing unit 93 and the corresponding position of phase place, reference axis is to displacement components u _y, and this displacement outputed to axial strain computing unit 97.96 pairs of horizontal direction displacement components u of receiving from horizontal direction computing unit 94 of horizontal direction strain computing unit _xCarry out space differentiation, so that the stress distribution ε of calculated level direction _x, and this stress distribution outputed to quantifying unit 98.On the other hand, 97 couples of axial displacement u that receive from axial computing unit 95 of axial strain computing unit _yCarry out space differentiation, so that reference axis is to stress distribution ε _y, and this stress distribution outputed to quantifying unit 98.98 pairs of horizontal direction stress distribution of quantifying unit ε _xWith axial strain distribution ε _yQuantize, show (or colored demonstration) so that these stress distribution are carried out GTG, and on display unit 99 demonstration information.Display unit 99 shows quantized like this each stress distribution.

Next, will the operation of the expansion combination autocorrelation method that adopts in the compuscan shown in Figure 7 be described.At first, consider such a case, wherein tissue is by compression slightly (be a few percent or still less).In this case, from local viewpoint, assumed stress has caused the translation of in-house every bit.Just, by the model shown in the following formula, represent the RF signal under the situation that tissue is exerted pressure and do not exerted pressure.

At this, i ₁(t, x) the representative RF signal under the situation of not exerting pressure, i ₂((t x) represents envelope signal, ω to A for t, x) the representative RF signal under the situation of exerting pressure ₀Represent hyperacoustic central angle frequency, τ represents time shift, and this time shift is as the time parameter of representing axial displacement, u _xRepresent the horizontal direction displacement, and θ represents initial phase.Notice that different with the combination autocorrelation method with Doppler (Doppler) method, about this method, the RF signal under the situation of exerting pressure and do not exert pressure is represented by the model of considering the horizontal direction displacement.In the end the parameter of stage acquisition is axial displacement u _y=c τ/2 (being time shift τ) and horizontal direction displacement components u _xNotice that c represents the in-house velocity of sound, and supposition c in tissue is a constant.

Then, the RF signal under the quadrature detector 92 subtend tissues situation of exerting pressure and not exerting pressure carries out orthogonal detection.Just, the sine wave and the cosine wave that will have the frequency identical with the ultrasound wave mid frequency are applied to the RF signal, further the RF signal are carried out low-pass filtering subsequently, obtain complex baseband signal s thus ₁And s ₂, shown in (14).

s ₁(t，x)＝A(t，x)e ^jθ

$s_2(t,x)=A(t-\mathrm{\&tau;},x-u_x)e^{j({\mathrm{\&omega;}}_0\mathrm{\&tau;}+\mathrm{\&theta;})}\&CenterDot;\&CenterDot;\&CenterDot;\left(14\right)$

Then, with s ₁(t, x) and s ₂(t+nT/2, x+mL) the two-dimentional compound correlative function R between ₁₂(t, x; N m) is defined as shown in the formula shown in (15).

$R_{12}(t,x;n,m)={\&Integral;}_{{-\mathrm{<figure-callout\; id="0"\; label="x"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">x</figure-callout>}}_0/2}^{x_0/2}{\&Integral;}_{{-\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}^{{\mathrm{<figure-callout\; id="0"\; label="t"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">t</figure-callout>}}_0/2}{s}_1(t+v,x+w)s_2{(t+\mathrm{<figure-callout\; id="2"\; label="n\; T"\; filenames="A200910006835E00521.png,A200910006835E00531.png"\; state="\{\{state\}\}">n</figure-callout>}\frac{T}{}\frac{}{2}+v,x+\mathrm{mL}+w)}^{*}\mathrm{dvdw}$

$=R_A(t,x;t-\mathrm{<figure-callout\; id="2"\; label="n\; T"\; filenames="A200910006835E00521.png,A200910006835E00531.png"\; state="\{\{state\}\}">n</figure-callout>}\frac{T}{}\frac{}{2},u_x-\mathrm{mL})e^{-{\mathrm{j\&omega;}}_0(\mathrm{\&tau;}-n\frac{T}{2})}\&CenterDot;\&CenterDot;\&CenterDot;\left(15\right)$

(n＝-N _min，...，-2，-1，0，1，2，...，N _max)

(m＝-M _min，...，-2，-1，0，1，2，...，M _max)

At this, T represents hyperacoustic cycle, and L represents sampling interval (bunch spacing), R _A(t, x; τ, u _x) represent the auto-correlation function of envelope signal, t ₀Represent the axial length of correlation window, x ₀Represent the cross-directional length of correlation window.On the other hand, the v representative is used for time (τ) the axial variate-value of integration, and the w representative is used for the variate-value of the bunch direction of integration, and complex conjugation operator represented in asterisk " * ".Then, utilize, define two-dimentional envelope correlation coefficient C (t, x as shown in the formula the two-dimentional compound correlative function shown in (16); N, m).

$C(t,x;n,m)=\frac{\left|R_{12}(t,x;n,m)\right|}{\sqrt{|R_{11}(t,x;\mathrm{0,0})|\&CenterDot;|R_{22}(t,x;n,m)|}}\&CenterDot;\&CenterDot;\&CenterDot;\left(16\right)$

Note R ₁₁(t, x; 0,0) represents S ₁(t, auto-correlation function x), and R ₂₂(t, x; N m) represents S ₂(t+nT/2, auto-correlation function x+mL).As the combination autocorrelation method, envelope correlation coefficient is used to solve the aliasing problem.Just, (t x) locates, and for all variable n and m, all obtains (t, x by C at each measurement point; N, m) and φ (t, x; N, the m) combination of Gou Chenging { C (t, x; N, m), φ (t, x; N, m) }.About this method, suppose in enough scopes, promptly in being enough to carry out the relevant scope of envelope, determine variable number to (n, m).In this case, with demonstrate the largest enveloping correlation coefficient (n, m)=(k, l) corresponding phase (t, x; K, l) phase place of the no aliasing of coupling.Its reason is, if the variable number that will demonstrate the largest enveloping correlation coefficient to (n, m) be expressed as (k, l), s then ₁(t, k) and s ₂(t+kT/2, x+lL) time shift between | τ-kT/2| is less than T/2.That is, | φ (t, x; K, l) |=ω ₀| t-kT/2| is less than π.Just, about this method, utilize phase (t, the x of no aliasing; K l) calculates, and (t x) locates to obtain correct time shift τ, correct axial displacement u at each measurement point thus _yAnd correct horizontal direction displacement components u _x, shown in (17).

$\mathrm{\&tau;}=-\frac{\mathrm{\&phi;}(t,x;k,l)}{{\mathrm{\&omega;}}_0}+\mathrm{<figure-callout\; id="2"\; label="k\; T"\; filenames="A200910006835E00521.png,A200910006835E00531.png"\; state="\{\{state\}\}">k</figure-callout>}\frac{T}{}\frac{}{2}$

$u_y=\frac{\mathrm{c\&tau;}}{2}$

u _x＝lL ...(17)

Notice that c represents the in-house velocity of sound (be assumed to the constant of 1500m/s, 1500m/s is the normal velocity of sound in the soft tissue).About this method, all measurement points are carried out aforementioned calculation, obtain axial displacement distribution u thus _y(x is y) with horizontal direction Displacements Distribution u _x(x, y).

In addition,, each of above-mentioned Displacements Distribution is carried out space differentiation, obtain axial strain distribution ε thus about this method _y(x is y) with horizontal direction stress distribution ε _x(x, y), shown in (18).

${\mathrm{\&epsiv;}}_y(x,y)=\frac{{\&PartialD;u}_y(x,y)}{\&PartialD;y}$ ...(18)

${\mathrm{\&epsiv;}}_x(x,y)=\frac{{\&PartialD;u}_x(x,y)}{\&PartialD;y}$

As mentioned above, about this method,, estimate axially and the displacement (strain) of horizontal direction distribution according to the RF signal under the situation of exerting pressure and not exerting pressure to tissue.Note, as can be from following formula u _x=lL finds out that the displacement detecting precision of horizontal direction is subjected to the restriction in the sampling interval (bunch spacing) of horizontal direction, and therefore, this method has the lower shortcoming of horizontal direction precision.Yet this method has the advantage of Real Time Observation, has improved actual performance thus.

Above-mentioned expansion combination autocorrelation method has following function: for calculating each time, the two-dimensional correlation window that utilization is used in preset range, analyze relative displacement, handle the horizontal direction displacement of tissue thus in the horizontal direction under the situation of exerting pressure, between tissue and the ultrasound probe to tissue.Yet, this two-dimensional expansion combination autocorrelation method with this function can not estimate, under the situation of exerting pressure to tissue, comprise with axial and all orthogonal direction of horizontal direction (with the orthogonal direction of the two-dimensional ultrasound plane of scanning motion (slice direction)) on the stress distribution of displacement.In order to address the above problem, at an easy rate above two-dimensional expansion combination autocorrelation method is expanded to the three-dimensional extended combination autocorrelation method that utilizes the three dimensional window that is applied to three-dimensional scope, make the system can be more stable thus.

Fig. 9 and Figure 10 are the flow charts of describing the rudimentary algorithm of three-dimensional arrangement autocorrelation method.Notice that Figure 10 is the flow chart that the part to processing shown in Figure 9 is described in detail.

In step S11,1 storage " 1 " of scanning line depositor, scanning line depositor 1 usefulness acts on the enumerator of execution for the identical processing of first scanning line to the L scanning line.The value of storage is examined in definite processing of step S23 in the scanning line depositor 1.

At step S12,, the variable on the thickness direction (frame direction) is increased progressively between-U to U for each circular treatment.Note, in definite processing of step S18, check counting.

At step S13,, the variable on the scanning direction (scan-line direction) is increased progressively between-V to V for each circular treatment.Note, in definite processing of step S17, check counting.

At step S14, for each circular treatment, make on the range direction (axially) variable 0 and M between increase progressively.Note, in definite processing of step S16, check counting.

At step S15, utilize the combination autocorrelation method, correlation coefficient C (l, the t of the envelope signal of computed range direction (axially); U, v, n).Notice that conventional combination autocorrelation method can not show the computational speed that is enough to carry out three-dimensional computations, and therefore, utilize high speed to make up autocorrelation method and calculate correlation coefficient C (l, t; U, v, n).The high speed autocorrelation method will be described in the back.

At step S16, the variable of determining in above-mentioned steps S14 is handled.That is, determine whether the variable n that stores in the range direction depositor has reached predetermined maximum M.If determine that variable n has reached M, then flow process advances to step S17.Otherwise flow process turns back to step S14, and the variable n that stores in the range direction depositor is increased progressively.

At step S17, the variable of determining in above-mentioned steps S13 is handled.That is, determine whether the variable v that stores in the depositor of scanning direction has reached predetermined maximum V.If determine that variable v has reached V, then flow process advances to step S18.Otherwise flow process turns back to step S13, and the variable v that stores in the depositor of scanning direction is increased progressively.

At step S18, the variable of determining in above-mentioned steps S12 is handled.That is, determine whether the variable u that stores in the thickness direction depositor has reached predetermined maximum U.If determine that variable u has reached U, then flow process advances to step S19.Otherwise flow process turns back to step S12, and the variable u that stores in the depositor of scanning direction is increased progressively.

At step S19, system u=(U ..., 0 ..., U), v=(V ..., 0 ..., V) and n=(0,1 ..., in scope N), search shows maximum correlation coefficient C (l, the t that calculates in step S12 to S18; U, v, variable combination n) (u, v, n).Note,, and adopted this range direction scope n=(0 though the present embodiment supposition has only normal pressure to put on tissue, 1, ..., N), but much less, normal pressure and negative pressure all are applied in the structural configuration therein, should adopt range direction scope n=(M ..., 0, ..., N).

At step S20, calculate the correlation coefficient C (l, the t that in step S19, obtain; U, v, phase difference n) (l, t; u ₀, v ₀, n ₀).

At step S21, phase contrast n ₀π+φ (l, t; u ₀, v ₀, n ₀) calculated, as the phase contrast of terminal stage.

At step S22, utilize point (u ₀, v ₀, n ₀) near point (u, v, the correlation coefficient C (l, the t that n) locate; U, v n), calculates the displacement v=v of scanning direction ₀Displacement components u=the u of+Δ v and thickness direction ₀+ Δ u.

At step S23, check the variable of in above-mentioned steps S11, storing in the scanning line depositor 1.That is, determine whether variable 1 has reached L.If determine that 1 has reached L, then flow process advances to step S24.Otherwise flow process turns back to step S11, and makes the variable increment of storage in the scanning line depositor 1.

At step S24,, calculate stress distribution under the situation of exerting pressure to tissue by the Displacements Distribution of estimating is carried out space differentiation.

Figure 10 is the flow chart that the above-mentioned high speed combination autocorrelation method among the step S15 shown in Fig. 9 is described in detail.

At step S31, the RF signal under the situation of not exerting pressure and exert pressure is carried out orthogonal detection, obtain I signal and Q signal thus, as described below.

X (t) → Ix, Qx (X (t)=Ix+jQx) wherein

Y (t) → Iy, Qy (Y (t)=Iy+jQy) wherein

At step S32, calculate relevant Rxy, Rxx and Ryy according to following formula.

Rxy＝∫X(t+v)·Y*(t+v)dv

Rxx＝∫X(t+v)·X*(t+v)dv

Ryy＝∫Y(t+v)·Y*(t+v)dv

At step S33, utilize following formula, calculate correlation coefficient C according to relevant Rxy, Rxx and the Ryy of such acquisition ₀

${\mathrm{<figure-callout\; id="0"\; label="C"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">C</figure-callout>}}_0=\left|\mathrm{Rxy}\right|/\sqrt{\mathrm{Rxx}}\&CenterDot;\sqrt{\mathrm{Ryy}}$

At step S34, variable n is set to 1.

At step S35, calculate Y _n(t)=Y (t-nT) e ^{J ω nT}

At step S36, calculate Rxy according to following formula _nAnd Ry _ny _n

Rxy _n＝∫X(t+v)·Y _n*(t+v)dv

＝∫X(t+v)·Y*(t-nT+v)e ^jωnTdv

Ry _ny _n＝∫Y _n(t+v)·Y _n*(t+v)dv

＝∫Y(t-nT+v)·Y*(t-nT+v)dv

At step S37, utilize following formula, according to Rxy _nAnd Ry _ny _nCalculate correlation coefficient C _n

$C_n=\left|{\mathrm{Rxy}}_n\right|/\sqrt{\mathrm{Rxx}}\&CenterDot;\sqrt{{\mathrm{Ry}}_n{y}_n}$

At step S38, variable n is increased progressively.

At step S39, determine whether variable n has reached predetermined maximum M.If determine that n has reached M, then processing finishes.Otherwise flow process turns back to step S35, and the complex phase of laying equal stress on is with handling.

Utilize the method shown in the flow chart among Figure 10, in step S35, calculate Y from Y _n, so that calculate Rxy _nAnd Ry _ny _nThis allows simple circuit configuration.Below will describe from Y and calculate Y _nMethod.

At first, represent echo-signal x (t) under the situation of not exerting pressure with following formula.

x(t)＝u(t)cos(ωt+θ)

On the other hand, be illustrated in echo-signal y (t) under the situation of axially exerting pressure with following formula.

y(t)＝x(t+τ)＝u(t+τ)cos(ω(t+τ)+θ)。

Then, represent x, y and y with following formula _nThe orthogonal detection signal.

x(t)→

Ix＝0.5u(t)cosθ

Qx＝-0.5u(t)sinθ

(X(t)＝Ix+jQx＝0.5u(t)e ^-jθ)

y(t)→

Iy＝0.5u(t+τ)cos(ωτ+θ)

Qy＝-0.5u(t+τ)sin(ωτ+θ)

(Y(t)＝Iy+jQy＝0.5u(t+τ)e ^-j(ωτ+θ))

yn(t)＝y(t-nT)

＝u(t+τ-nT)cos(ω(t+τ-nT)+θ)

＝u(t+τ-nT)cos(ωt+ω(τ-nT)+θ)

At this, T represents for two/one-period.Then, obtain the Iy be shown below _nAnd Qy _n

Iy _n＝0.5u(t+τ-nT)cos(ω(τ-nT)+θ)

Qy _n＝-0.5u(t+τ-nT)sin(ω(τ-nT)+θ)

(Y _n＝Iy _n+Qy _n＝0.5u(t+τ-nT)e ^{-j(ω(τ-nT)+θ)})

Therefore, obtained following relational expression according to following formula.

Y _n(t)＝Iy _n+Qy _n

＝0.5u(t+τ-nT)e ^{-j(ω(τ-nT)+θ)}

＝Y(t-nT)e ^jωnT

As calculating Y from Y (t)=Iy+jQy as seen from the above equation _n(t).

Therefore, calculate Rxy from X and Y _nAnd Ry _ny _n, be shown below.

Rxy _n＝4∫X(t+v)·Y _n*(t+v)dv

＝4∫X(t+v)·Y*(t-nT+v)e ^jωnTdv

|Rxy _n|＝Ru _n

＝∫u(t+v)u(t+τ-nT+v)dv

＝4∫|X(t+v)·Y _n*(t+v)|dv

＝4∫|X(t+v)·Y*(t-nT+v)e ^jωnT|dv

＝4∫|X(t+v)·Y*(t-nT+v)|dv

Ry _ny _n＝∫u(t+τ-nT+v)u(t+τ-nT+v)dv

＝4∫|Y _n(t+v)·Y _n*(t+v)|dv

＝4∫Y(t-nT+v)·Y*(t-nT+v)dv

At this, complex conjugation operator represented in asterisk " * ".

Figure 11 has shown the block diagram of the circuit configuration of the rudimentary algorithm that is used to carry out the three-dimensional arrangement autocorrelation method.About being used to carry out the custom circuit configuration of combination autocorrelation method as shown in Figure 7, be provided with a large amount of orthogonal demodulation circuit 1320 to 132N, and the processing in these orthogonal demodulation circuits 1320 to 132N needs a large amount of processing times, caused to be difficult to carry out supercomputing, caused being difficult to carrying out real time imaging and shown.Therefore, present embodiment is used to carry out the circuit configuration of above-mentioned rudimentary algorithm as shown in figure 11, allows operation to be used to carry out the high speed circuit of combination autocorrelation method thus.

Pressurized orthogonal demodulation circuit (QD) 131 does not receive the echo-signal x (t) under the situation of not exerting pressure, the echo-signal of receiving is carried out orthogonal detection, and the signal Ix (t) of orthogonal detection and Qx (t) are outputed to the first correlation computations circuit 133 and the first correlation coefficient counting circuit 1350 to 135N.The echo-signal y (t) that pressurized orthogonal demodulation circuit (QD) 132 receives under the situation of exerting pressure, the echo-signal of receiving is carried out orthogonal detection, and signal Y (the t)=Iy+jQy (Iy (t), Qy (t)) of orthogonal detection outputed to the first correlation coefficient counting circuit 1350, the second correlation computations circuit 1340, first delay circuit 134 and second delay circuit 135.First delay circuit 134 and second delay circuit 135 make orthogonal detection signal Y (t) postpone hyperacoustic period T, and the orthogonal detection signal Y (t-T) that postpones is outputed to the first correlation computations circuit 1351, the 3rd delay circuit 136 and the 4th delay circuit 137.Each of the 3rd delay circuit 136 and the 4th delay circuit 137 all makes orthogonal detection signal Y (t-T) postpone hyperacoustic period T, and the orthogonal detection signal Y (t-2T) that postpones is outputed to subsequently the first correlation coefficient counting circuit and delay circuit (not shown).In the same way, utilize each of N delay circuit, make orthogonal detection signal delay period T continuously, and the signal that postpones is outputed to the corresponding first correlation coefficient counting circuit.

The first correlation computations circuit 133 calculates correlation Rxx according to signal Ix and Qx, and correlation Rxx is outputed to the second correlation coefficient counting circuit 1380 to 138N.The second correlation computations circuit 1340 receives orthogonal detection signal Iy (t) and Qy (t) from pressurized orthogonal demodulation circuit 132, calculates correlation Ryy according to signal Iy and Qy, and correlation Ryy is outputed to the second correlation coefficient counting circuit 1380.The first correlation coefficient counting circuit 1350 never pressurized orthogonal demodulation circuit 131 receives orthogonal detection signal Ix (t) and Qx (t), and receives orthogonal detection signal Iy (t) and Qy (t) from pressurized orthogonal demodulation circuit 132, calculates complex baseband signal S _RAnd S _I, and with complex baseband signal S _RAnd S _IOutput to third phase and close counting circuit 1360 and phase difference calculating circuit 1370.Third phase closes counting circuit 1360 and receives complex baseband signal S from the first correlation coefficient counting circuit 1350 _RAnd S _I, according to complex baseband signal S _RAnd S _ICalculate correlation | Rxy|, and with correlation | Rxy| outputs to the second correlation coefficient counting circuit 1380.Phase difference calculating circuit 1370 receives complex baseband signal S from the first correlation coefficient counting circuit 1350 _RAnd S _I, and according to complex baseband signal S _RAnd S _ICalculate phase difference ₀(t).The second correlation coefficient counting circuit 1380 receives correlation Rxx from the first correlation computations circuit 133, close counting circuit 1360 from third phase and receive correlation | Rxy|, and from the second correlation computations circuit, 1340 reception correlation Ryy, according to these correlation value calculation correlation coefficienies C ₀And output correlation coefficient C (t), ₀(t).

The second correlation computations circuit 1341 calculates correlation Ry from the orthogonal detection signal Iy (t-T) and the Qy (t-T) of first delay circuit 134 and second delay circuit, 135 receive delays according to signal Iy (t-T) and Qy (t-T) ₁y ₁, and with correlation Ry ₁y ₁Output to the second correlation coefficient counting circuit 1381.The first correlation coefficient counting circuit 1351 never pressurized orthogonal demodulation circuit 131 receives orthogonal detection signal Ix (t) and Qx (t), and, obtain complex baseband signal S from the orthogonal detection signal Iy (t-T) and the Qy (t-T) of first delay circuit 134 and second delay circuit, 135 receive delays _R1And S _I1, and with complex baseband signal S _R1And S _I1Output to third phase and close counting circuit 1361 and phase difference calculating circuit 1371.Third phase closes counting circuit 1361 and receives complex baseband signal S from the first correlation coefficient counting circuit 1351 _R1And S _I1, according to complex baseband signal S _R1And S _I1Obtain correlation | Rxy ₁|, and with correlation | Rxy ₁| output to the second correlation coefficient counting circuit 1381.Phase difference calculating circuit 1371 receives complex baseband signal S from the first correlation coefficient counting circuit 1351 _R1And S _I1, and according to complex baseband signal S _R1And S _I1Calculate phase difference ₁(t).The second correlation coefficient counting circuit 1381 receives correlation Rxx from the first correlation computations circuit 133, closes counting circuit 1361 from third phase and receives correlation | Rxy ₁|, and from the second correlation computations circuit, 1341 reception correlation Ry ₁y ₁, according to these correlation value calculation correlation coefficienies C ₁And output correlation coefficient C (t), ₁(t).

In the same way, utilization is closed counting circuit 1362 to 136N, phase difference calculating circuit 1372 to 137N and the second correlation coefficient counting circuit 1382 to 138N from the second correlation computations circuit 1342 to 134N, the first correlation coefficient counting circuit 1352 to 135N, the third phase of the 3rd delay circuit 135 and the 4th delay circuit 136 arrangement downwards, to the orthogonal detection signal Iy (t-2T) of continuous delay ..., Iy (t-NT) and Qy (t-2T) ..., Qy (t-NT) each carry out aforesaid same treatment, export correlation coefficient C thus ₂(t) to C _N(t) and phase difference ₂(t) to φ _N(t).

Next, description is utilized the elastic modelling quantity distribution reconstructing method of dimensional finite element method.To be used for the inversion problem that the reconstruct elastic modelling quantity distributes in order simplifying,, to represent tissue with model about present embodiment.The elastic modelling quantity distribution reconstructing method that this allows user to utilize Finite Element Method to carry out to propose in the present embodiment.Particularly, about present embodiment, represent tissue by model based on following supposition.

At first, suppose that tissue shows isotropic elasticity.On the other hand, to being in the tissue under the outside static pressure, estimate to organize strained distribution.Note, tissue is applied static pressure,, allow the user to utilize the relation between the RF signal under the RF signal under tissue is exerted pressure situation and the situation of not exerting pressure to calculate thus to impel the mild compression of tissue.Therefore, in this case, the relation between the stress and strain demonstrates linear relationship.Just, can suppose and represent tissue, carry out proximate calculating with elastic model.In addition, in this case, suppose that tissue is isotropic, and therefore, suppose that tissue shows isotropic elasticity.

In addition, suppose that tissue shows near Incoercibility.In general, well-known, tissue shows the compressibility near Incoercibility (Poisson is than (Poisson's ratio) v=0.5).About present embodiment, utilize in 0.49 the tissue constant Poisson's ratio to calculate.Noting, about present embodiment, is not to calculate under tissue shows the supposition of complete Incoercibility.Its reason is, calculate if show at tissue under the supposition of complete Incoercibility, promptly utilize in 0.5 the tissue constant Poisson's ratio to calculate, then elastic equation will become special circumstances, the problem that has caused utilizing the Finite Element Method that proposes in the present embodiment to be calculated.In addition,, suppose that Poisson's ratio is constant in tissue, and therefore, only should be elastic modelling quantity distribution estimating Young modulus (Young's modulus), simplified inversion problem thus about present embodiment.Note, in general,, utilize in 0.49 the tissue constant Poisson's ratio to calculate about present embodiment.

Next, represent tissue with three-dimensional finite element model.Elastic modelling quantity distribution reconstructing method about according to present embodiment utilizes Finite Element Method to calculate, and therefore, tissue is divided into the rectangular parallelepiped protrusion part unit of limited quantity.Note, suppose that each unit shows constant elastic modelling quantity, constant stress and constant strain therein.In general, in order to find out the method that solves inversion problem, find out with the corresponding forward problem of inversion problem be important.About present embodiment, inversion problem is to estimate that according to stress distribution and boundary condition elastic modelling quantity distributes.Therefore, be to calculate stress distribution according to elastic modelling quantity distribution and boundary condition with the corresponding forward problem of above-mentioned inversion problem.About present embodiment, utilize Finite Element Method (FEM: Finite Element Method), solve the problems referred to above as a kind of known numerical analysis method.

About Finite Element Method, by the model that the unit by limited quantity constitutes, represent to serve as the tissue of the continuum that will estimate, and utilize numerical analysis to find the solution to represent the simultaneous linear equations of each unitary characteristic.Notice that the specific calculations of utilizing Finite Element Method will be described in the back.In brief,, distribute and boundary condition, obtain to distribute and stress distribution as the strain (displacement) of output valve according to tissue elasticity modulus as input value about Finite Element Method.

About present embodiment, show at tissue under the supposition of isotropic elasticity and carry out approximate calculation, and therefore, elastic equation (relation between equilibrium equation, strain and the displacement, the relation between the stress and strain) keeps under the condition in the tissue that is shown below.

Represent equilibrium equation with following formula (19).

$\begin{array}{cc}\frac{{\&PartialD;\mathrm{\&sigma;}}_{\mathrm{ij}}}{{\&PartialD;x}_j}+f_i=0& (i,j=\mathrm{1,2,3})\end{array}\&CenterDot;\&CenterDot;\&CenterDot;\left(19\right)$

At this, serve as the reference number 1,2 and 3 of i and j and represent x, y and z respectively.On the other hand, represent relation between strain and the displacement with following formula (20).

${\mathrm{\&epsiv;}}_{\mathrm{ij}}=\frac{1}{2}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})\&CenterDot;\&CenterDot;\&CenterDot;\left(20\right)$

Represent relation (unitized Hooke rule) between the stress and strain with following formula (21)

${\mathrm{\&sigma;}}_{\mathrm{ij}}=\frac{\mathrm{<figure-callout\; id="1"\; label="E"\; filenames="A200910006835E00531.png,A200910006835E00561.png"\; state="\{\{state\}\}">E</figure-callout>}}{1+v}({\mathrm{\&epsiv;}}_{\mathrm{ij}}+\frac{v}{1-2v}{\mathrm{\&delta;}}_{\mathrm{ij}}{\mathrm{\&epsiv;}}_{\mathrm{nn}})\&CenterDot;\&CenterDot;\&CenterDot;\left(21\right)$

Utilize the above-mentioned elastic equation of tensor representation.Therefore, in fact three equilibrium equations are arranged, six strain-displacement relation formulas, and six stress-strain relations.Note coordinate (x ₁, x ₂, x ₃) representative (x, y, z).The following characteristic of other reference character representative.

E: Young's modulus (being also referred to as " elastic modelling quantity ")

V: Poisson's ratio

ε _Ij: strain tensor

(ε nn=ε ₁₁+ ε ₂₂+ ε ₃₃: bulk strain)

s _Ij: stress tensor

δ _Ij: Kronecker delta (kronecker delta)

u _i: motion vector

f _i: the body force vector

(notice that it is negligible that gravity is considered to, and therefore, suppose f in this case _iBe 0).

Then, be ε _IjFind the solution the relational expression between the stress and strain.As a result, the relational expression between strain and the stress is changed into as shown in the formula shown in (22).

${\mathrm{\&epsiv;}}_{\mathrm{ij}}=\frac{1+v}{E}({\mathrm{\&sigma;}}_{\mathrm{ij}}-\frac{\mathrm{<figure-callout\; id="1"\; label="v"\; filenames="A200910006835E00531.png,A200910006835E00561.png"\; state="\{\{state\}\}">v</figure-callout>}}{1+v}{\mathrm{\&delta;}}_{\mathrm{ij}}{\mathrm{\&sigma;}}_{\mathrm{nn}})\&CenterDot;\&CenterDot;\&CenterDot;\left(22\right)$

At this, s _Nn=s ₁₁+ s ₂₂+ s ₃₃So, with i=j=2 substitution following formula (22), and following formula found the solution Young's modulus E, obtain following formula (23) thus.

$E=\frac{{\mathrm{\&sigma;}}_{22}-v({\mathrm{\&sigma;}}_{11}+{\mathrm{\&sigma;}}_{33})}{{\mathrm{\&epsiv;}}_{22}}$ ...(23)

$=\frac{{\mathrm{\&sigma;}}_y-v({\mathrm{\&sigma;}}_x+{\mathrm{\&sigma;}}_z)}{{\mathrm{\&epsiv;}}_y}$

As finding out, can calculate Young's modulus according to the components of strain that axially (are the ultrasonic beam direction in the present embodiment) and the components of stress of all directions, i.e. elastic modelling quantity from following formula (23).Yet, be difficult to utilize current techniques that the stress distribution that is used for following formula is directly measured.Therefore,, alternately estimate and upgrade stress distribution and elastic modelling quantity to distribute, approach actual distribution so that estimated elastic modelling quantity distributes to become about present embodiment.As the execution of getting off is used for the particular procedure that the reconstruct elastic modelling quantity distributes.

At first, suppose uniform distribution with acting on the initial distribution { E^ that estimates that elastic modelling quantity distributes ⁰.The second, utilize dimensional finite element method, calculate { E^ because the initial elasticity modulus distributes ⁰Caused stress distribution { S^ ⁰.Particularly, in strain-displacement relation formula and the equilibrium equation more than the stress-strain relation substitution, form new equilibrium equation thus, shown in (24).The equilibrium equation that this is new is applied to each unit in the organize models.

$\frac{\&PartialD;}{{\&PartialD;x}_j}[\frac{E}{2(1+v)}(\frac{{\&PartialD;u}_i}{{\&PartialD;x}_j}+\frac{{\&PartialD;u}_j}{{\&PartialD;x}_i})+\frac{\mathrm{Ev}}{(1+v)(1-2v)}{\mathrm{\&delta;}}_{\mathrm{ij}}\frac{{\&PartialD;u}_n}{{\&PartialD;x}_n}]+f_i=0$

At this,

$\frac{{\&PartialD;u}_n}{{\&PartialD;x}_n}=\frac{{\&PartialD;u}_1}{{\&PartialD;x}_1}+\frac{{\&PartialD;\mathrm{<figure-callout\; id="2"\; label="u"\; filenames="A200910006835E00521.png,A200910006835E00531.png"\; state="\{\{state\}\}">u</figure-callout>}}_2}{{\&PartialD;x}_2}+\frac{{\&PartialD;\mathrm{<figure-callout\; id="3"\; label="u"\; filenames="A200910006835E00592.png,A200910006835E00593.png"\; state="\{\{state\}\}">u</figure-callout>}}_3}{{\&PartialD;x}_3}\&CenterDot;\&CenterDot;\&CenterDot;\left(25\right)$

Then, with under the downstream condition, utilize Gaussian (Gauss) elimination approach that simultaneous equations are found the solution displacement, obtain thus and the elastic modelling quantity { E^ that distributes ⁰Corresponding Displacements Distribution { u^ ⁰.

u _i| _{The y=bottom}=0

σ _i| _{The y=top}=p _i... (26)

σ _n| _{X, the z=side}=0

In following formula, p _iRepresent the external pressure vector on the body surface, and s _nThe representative and the orthogonal components of stress in side.First formula is represented the fixed boundary condition in bottom, and second formula is represented the boundary condition that the stress distribution coupling external pressure of body surface distributes, and the 3rd formula represents that each side has free-ended boundary condition.Next, with Displacements Distribution { u^ ⁰In the substitution strain-displacement relation formula, obtain thus and the elastic modelling quantity { E^ that distributes ⁰Corresponding stress distribution { ε ^ ⁰.Then, with stress distribution { ε ^ ⁰In the substitution stress-strain relation, obtain thus and the elastic modelling quantity { E^ that distributes ⁰Corresponding stress distribution { S^ ⁰.

The 3rd, utilize following formula (27), according to axial (y direction) stress distribution { ε of the stress distribution of utilizing three dimension finite element method and the estimation of utilization expansion combination autocorrelation method _y, upgrade elastic modelling quantity distribution { E^ ^k.

${\hat{E}}^{k+1}=\frac{{{\widehat{\mathrm{\&sigma;}}}_y}^k-v({{\widehat{\mathrm{\&sigma;}}}_x}^k+{{\widehat{\mathrm{\&sigma;}}}_z}^k)}{{\mathrm{\&epsiv;}}_y}\&CenterDot;\&CenterDot;\&CenterDot;\left(27\right)$

Note,, obtain following formula by under the situation of i=j=2, above-mentioned stress-strain relation is found the solution Young's modulus E.In following formula, the k representative is code name repeatedly.

The 4th, elastic modelling quantity distribution and above-mentioned boundary condition according to such renewal carry out three-dimensional finite element analysis continuously, upgrade stress distribution thus.

Then, carry out third and fourth continuously and handle, elastic modelling quantity distributes and distributes near actual elastic modelling quantity thus.Note, satisfy following formula (28), determine that then estimated elastic modelling quantity distribution has reached convergence, and finish to estimate to handle if elastic modelling quantity distributes.

$\frac{1}{N}\underset{l}{\overset{N}{\mathrm{\&Sigma;}}}|{{\hat{E}}_l}^{k+1}-{{\hat{E}}_l}^k|<\mathrm{\&Gamma;}\&CenterDot;\&CenterDot;\&CenterDot;\left(28\right)$

At this, 1 representative unit number, the sum of N representative unit, and Γ represents threshold value.

The elastic modelling quantity distribution reconstructing method that utilizes three-dimensional finite element model that proposes in the present embodiment has been described.About this method, estimate that according to the three-dimensional balancing equation elastic modelling quantity distributes.Therefore,, under the supposition that more approaches actual tissue than conventional method, calculate, make it possible to thus estimate more accurately that elastic modelling quantity distributes about this method.In addition, about this method, only distribute according to the axial strain that can estimate accurately, come the reconstruct elastic modelling quantity to distribute, making it possible to reliably thus, the reconstruct elastic modelling quantity distributes.Note,, estimate three-dimensional elastic modulus for tissue and distribute about this method, and therefore, need to adopt two-dimensional array ultrasonic probe, perhaps need the one-dimensional array ultrasound probe of slice direction is carried out mechanical scanning, so that the tissue that will analyze is carried out 3-D scanning.

Below, with describe according to present embodiment, based on the advantage of the elastic modelling quantity distribution reconstructing method of expansion combination autocorrelation method and three-dimensional finite element model, utilize simulation to confirm the advantage of this method.Figure 12 has shown the sketch map of simulation process.

At first, create and to have the organize models that is used to estimate that the elastic modelling quantity tested distributes.Note, comprise scattering point in this organize models, these scattering point reflection supersonic wave echo-signals.The second, apply external pressure to organize models, so that in computer simulation, cause the compression of organize models.Then, utilize Finite Element Method etc., calculate the displacement of each scattering point that causes owing to compression.The 3rd, according to the scattering point Displacements Distribution under the situation of exerting pressure and not exerting pressure, simulate the RF signal under the situation of exerting pressure and do not exert pressure to organize models.The 4th, expansion combination autocorrelation method is applied to analog rf signal under the situation of exerting pressure and do not exert pressure, estimate to organize stress distribution thus.The 5th, by utilizing the elastic modelling quantity distribution reconstructing method of three-dimensional finite element model, according to utilizing expansion to make up the stress distribution of autocorrelation method estimation and being the definite boundary condition of the compression of simulated tissue model (external pressure distributes or the like), estimate that the tissue elasticity modulus distributes.

Though in this simulation, different elastic modelling quantity distributes and is used for organize models, suppose that various as used herein elastic modelling quantity distributions all show isotropic elasticity.Note general coupling of the elastic modelling quantity that uses in this simulation and elastic modelling quantity according to the corresponding breast tissue of main uses of the tissue elasticity modulus distribution measurement system of present embodiment.On the other hand, all comprise scattering point in each organize models, these scattering points are used to simulate the reflection RF signal under the situation of exerting pressure and not exerting pressure to tissue.Particularly, each organize models all to comprise average density be 500 points/cm ³Scattering point.In addition, utilizing average is 1 .0, standard deviation is 0.3 normal state pseudo random number, determines the scattering point position under the situation of not exerting pressure.Then, by according to the displacement of calculating each scattering point under the situation of not exerting pressure, obtain the scattering point position under the situation of exerting pressure from the result of finite element analysis.At this, though be unknown, determine each parameter of scattering point like this, so that the B mode image that generates according to mimic RF signal is similar to the B mode image of actual tissue about the information of the scattering point in the actual tissue.

About this method, carry out convolution by making the scattering point function under the situation of exerting pressure and do not exert pressure and the point spread function of ultrasonic system, calculate the analog rf signal under the situation of exerting pressure and not exerting pressure, shown in (29) to tissue.

i ₁(x，y，z)＝∫∫∫h(x-x′，y-y′，z-z′)t ₁(x′，y′，z′)dx′dy′dz′

...(29)

i ₂(x，y，z)＝∫∫∫h(x-x′，y-y′，z-z′)t ₂(x′，y′，z′)dx′dy′dz′

At this, i ₁(z) representative does not have the RF signal under the situation that tissue is exerted pressure, i for x, y ₂(x, y, z) the RF signal of representative under the situation that tissue is exerted pressure, (x, y z) represent the point spread function (impulse response) of ultrasonic system, t to h ₁(z) representative does not have the scattering point function under the situation that tissue is exerted pressure, and t for x, y ₂(x, y, z) the scattering point function of representative under the situation that tissue is exerted pressure.Notice that the scattering point function is represented predetermined like this scattering coefficient at each scattering point place, and represent scattering coefficient 0 in the position except that scattering point.On the other hand, according to because the displacement of each scattering point that the strain of organize models causes, never the scattering point function t under tissue is exerted pressure situation ₁(x, y z), calculate the scattering point function t under tissue is exerted pressure situation ₂(x, y, z).Note, carry out linear interpolation, calculate because the displacement of each scattering point that the compression of organizing causes by motion vector to the cell node place that utilizes finite element analysis computation.

About simulation, suppose the non-focusing ultrasonic system that adopts no ultrasound wave damping according to present embodiment.Just, and assumed position spread function h (x, y, z) all constant for every bit.In addition, suppose point spread function can be divided into as shown in the formula shown in (30), with respect to the function of all directions.

h(x，y，z)＝h _x(x)hy(y)h _z(z) ...(30)

At this, h _y(y) represent the point spread function of ultrasonic beam direction.On the other hand, h _x(x) and h _z(z) each representative and the point spread function on the orthogonal direction of ultrasonic beam direction.Particularly, suppose that the x direction represents direction in the plane of ultrasonic tomographic image (horizontal direction), and the representative of z direction is perpendicular to the direction (slice direction) of ultrasonic tomographic image.Note,, set up the point spread function of each direction according to utilizing ultrasonic device, distributing to carry out the echo-signal that actual measurement obtains from the echo-signal of line target (extend the water, diameter be the line (wire) of 0.13mm) reflection.Figure 13 has described and has been used to utilize the ultrasound wave mid frequency of 5.0MHz to carry out mimic each point spread function.Figure 13 (A) has shown axial point spread function h _y(y).By Gaussian (Gauss) function and sine wave are multiplied each other, obtain point spread function h _yAnd h (y), _y(y) APPROXIMATE DISTRIBUTION that distributes as actual reflection echo from the line target reflection.On the other hand, Figure 13 (B) has shown the point spread function h of horizontal direction _x(x), Figure 13 (C) has shown the point spread function h of slice direction _z(z).Note utilizing Gaussian (Gauss) function to set up each of above-mentioned point spread function, so that be used as the APPROXIMATE DISTRIBUTION that distributes from the actual reflection echo of line target reflection in the same way.Each function parameters becomes with mid frequency, will be in the back, in explanation each when simulation this is described.

Next, will describe, confirmed the advantage of this method by simulation as the advantage that makes up autocorrelation method according to the expansion of displacement (strain) distribution estimation method of present embodiment.At first, will describe the advantage of the displacement of tissue of estimating horizontal direction, this is the advantage that is better than making up autocorrelation method.

Figure 14 has shown organize models.Organize models is made of external dimensions and the even elastic modelling quantity distribution of 60mm * 60mm (two-dimensional).Then, the compression of simulated tissue is to cause axial 3% homogeneous strain.Suppose in this simulation, adopt simple one dimension elastic model, only estimate the organize models that the advantage of autocorrelation method is made up in expansion as being used for.In addition, utilize the axial compression of the shift simulation tissue of horizontal direction 0.0mm to 1.4mm, to estimate the advantage of handling horizontal direction displacement (about horizontal direction relative displacement ultrasound probe, tissue).In addition, suppose the simple parallel displacement of the tissue on the dummy level direction, the situation that this displacement is slided with respect to tissue corresponding to ultrasound probe.

Then, have the strain of organizing and do not organizing under the strained situation, be the RF of enterprise modeling signal.Being used for the mimic parameter of ultrasonic system comprises: the mid frequency of 5.0MHz, the pulse width of 0.5mm; 1.0mm the ultrasound wave beamwidth; 0.5mm trace interval; And the sample frequency of 30MHz.Then, utilize the expansion combination autocorrelation method that proposes in the present embodiment,, estimate stress distribution according to the analog rf signal under the situation of exerting pressure and do not exert pressure.In addition, prepare to utilize combination autocorrelation method stress distribution of estimating and the stress distribution of utilizing the space correlation method to estimate, to compare.Note, the correlation window of same size is applied to have each method of same search scope, allow the user between them, precision etc. simply to be compared thus.Particularly, expansion combination autocorrelation method and space correlation method adopt in the two-dimensional search scope that is applied to 5.6mm (axially) * 7.5mm (horizontal direction) and two-dimentional window that be of a size of 1.6mm (axially) * 2.5mm (horizontal direction).On the other hand, about being used to analyze the combination autocorrelation method of one dimension displacement, the one-dimensional correlation window application that will have identical 1.6mm axial length is in identical 5.6mm axial range.

Figure 15 to Figure 17 has shown the stress distribution estimated result that utilizes each method.Figure 15 has shown the strain error of the horizontal direction of utilizing each method estimation.At this, Reference numeral " ◇ " representative utilizes the error of combination autocorrelation method, "

" representative utilizes the error of expansion combination autocorrelation method, and " Δ " representative utilizes the error of spatial autocorrelation method.Note, represent the strained error e of estimating with following formula (31).

$e=\sqrt{\frac{\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{({\widehat{\mathrm{\&epsiv;}}}_i-{\mathrm{\&epsiv;}}_i)}^2}{\underset{i}{\overset{N}{\mathrm{\&Sigma;}}}{{\mathrm{\&epsiv;}}_i}^2}}\&CenterDot;\&CenterDot;\&CenterDot;\left(31\right)$

At this, ε ^ _iRepresent estimated strain, ε _iRepresent actual strain (ideal value), i representative unit number, and the sum of N representative unit.On the other hand, Figure 16 has shown, utilizes each method (combination autocorrelation method, expansion combination autocorrelation method and space correlation method) stress distribution that estimate, that comprise the tissue of 0.0mm horizontal direction displacement.Figure 17 has shown, utilizes each method (combination autocorrelation method, expansion combination autocorrelation method and space correlation method) stress distribution that estimate, that comprise the tissue of 0.4mm horizontal direction displacement.Notice that Figure 16 and Figure 17 have shown strained meansigma methods of estimation and the standard deviation with respect to vertically each degree of depth.

As finding out from analog result, combination autocorrelation method (CA method) about routine, the relative displacement of organizing that surpasses ultrasonic beam size (being 1/2nd beamwidths, i.e. 0.5mm in this case) on the horizontal direction has caused the rapid deterioration of strain estimated quality.On the other hand, as can be seen, expansion combination autocorrelation method can stably be estimated strain, and is irrelevant with the size of horizontal direction displacement.In addition, as can be seen, though the space correlation method also can stably be estimated strain, and irrelevant with the size of horizontal direction displacement, estimated result demonstrates relatively poor precision, promptly compares with expansion combination autocorrelation method, demonstrates twice or bigger error.In addition, processing time to said method compares, though the processing time that expansion combination autocorrelation method needs is 3.6 times of the combination autocorrelation method, the processing time that expansion combination autocorrelation method needs only is 1/ (7.7) times of space correlation method, and is as shown in the table.As finding out from The above results, expansion combination autocorrelation method can be realized real-time calculating to a certain extent.

Method	Processing time	The processing time of normalization
			The combination autocorrelation method	26 seconds	1/(3.6)
Expansion combination autocorrelation method	1 minute 34 seconds	1.0
			The space correlation method	12 minutes 5 seconds	7.7

Next, with the estimated result of describing under the situation of in an inclined direction exerting pressure.Note, different with the above-mentioned simulation that utilizes simple tissue model with two-dimentional homogeneous texture, estimate that about this organize models that utilization has as the three dimensional structure of actual tissue estimates.Note, in ideal measurement, should utilize ultrasound probe to go up and exert pressure to tissue in ultrasonic beam direction (axially).Now, description is subjected to the estimated result of the influence of the tissue compression on the incline direction.

Figure 18 has shown the organize models of the influence that is used to estimate the tissue compression on the incline direction.Shown in Figure 18 (A), this organize models has the three dimensional structure that external dimensions is 60mm * 60mm * 60mm, and to comprise diameter be that 15mm, length are the cylinder inclusions of 60mm, high rigidity.Suppose that the material that surrounds inclusions has the elastic modelling quantity (Young's modulus) of 10Kpa, and the elastic modelling quantity of inclusions is 3 times of above-mentioned material, promptly has the elastic modelling quantity of 30Kpa.Note,, determine that above-mentioned elastic modelling quantity, the diagnosis of breast carcinoma are main purposes of the present invention according to the elastic modelling quantity of breast tissue and the elastic modelling quantity of breast carcinoma.Then, under following two kinds of situations, the compression of simulated tissue model.A kind of situation corresponding to, 2% the axial organize models compression that applies to organize models that the even external pressure of 200Pa causes owing to vertically, from end face is shown in Figure 18 (B).Another kind of situation corresponding to, owing to along inclined direction, compress to the organize models that organize models applies the incline direction that even external pressure (axially 200Pa and horizontal direction 30Pa) causes from end face, shown in Figure 18 (C).

Then, for above-mentioned two kinds of situations, the RF signal of simulation under the situation of exerting pressure and do not exert pressure then utilizes expansion combination autocorrelation method to estimate stress distribution.Note, also utilize combination autocorrelation method and space correlation method to estimate stress distribution, to compare.At this, for each method, the correlation window that will have same size is applied to the interior calculating of same search scope, allows relatively simple thus.Notice that as above-mentioned simulation, correlation window has same size.In addition, other parameter that will be used for analog rf signal is defined as the value as above-mentioned mimic situation.Just, other parameter comprises: the mid frequency of 5.0MHz, the pulse width of 0.5mm; 1.0mm horizontal direction ultrasound wave beamwidth; 2.0mm slice direction ultrasound wave beamwidth; 0.5mm trace interval; And the sample frequency of 30MHz.

Figure 19 and Figure 20 have shown above-mentioned mimic analog result.Figure 19 shown, wherein the estimated result of the stress distribution of organize models under axial compressed simple case.Figure 20 shown, wherein the estimated result of the stress distribution of organize models under the compressed situation of incline direction.Notice that for each situation, the axial strain that obtains by three-dimensional finite element analysis is distributed as the stress distribution of utilizing Perfected process to estimate, this axial strain distributes and serves as actual stress distribution.Notice that Figure 19 and Figure 20 are that Figure 19 and Figure 20 have shown estimated result along organize models's viewgraph of cross-section that the center line of organize models is got.In Figure 20, the stress distribution of utilizing Perfected process to estimate shows left and right sides unsymmetry.Its reason is that pressure applies at incline direction.Particularly, in this case, pressure is that in the accompanying drawings lower right upwards applies, and has caused the big displacement on the horizontal direction of upper right portion in the accompanying drawing.

At first, confirmed that as the combination autocorrelation method, for the stress distribution that comprises axial compression, expansion combination autocorrelation method demonstrates general identical estimated result.On the other hand, in an inclined direction under the situation that tissue is exerted pressure, also confirmed, though because the big displacement of horizontal direction, can not utilize the combination autocorrelation method to estimate some zone and cause, but expansion combination autocorrelation method can stably be estimated stress distribution, and irrelevant with the displacement size of horizontal direction, and is described as above above-mentioned simulation.The space correlation method on the other hand, also confirmed, though can estimate stably that and irrelevant with the displacement size of horizontal direction, the estimated accuracy of space correlation method is more very different than expansion combination autocorrelation method.Except that above-mentioned analog result, confirmed that expansion combination autocorrelation method has the advantage that stress distribution is estimated.

As mentioned above, confirmed, can utilize above-mentioned expansion combination autocorrelation method, estimate to organize stress distribution at a high speed, accurately.To describe now, by the elastic modelling quantity distribution reconstructing method that proposes in the present embodiment, utilize three-dimensional finite element model is simulated the estimated result that obtains.Notice that elastic modelling quantity distribution reconstructing method is to estimate the method that elastic modelling quantity distributes according to stress distribution, this method is as the method for carrying out in the second level of tissue elasticity imaging system.

The major function of the elastic modelling quantity distribution reconstructing method that proposes in the present embodiment is to estimate that according to the three-dimensional dynamics equilibrium equation elastic modelling quantity distributes.Now, by elastic modelling quantity distribution reconstructing method and two-dimentional reconstructing method according to present embodiment are compared, confirm the advantage of the elastic modelling quantity distribution reconstructing method that proposes in the present embodiment, this two dimension reconstructing method has same processing except utilizing the two-dimensional dynamical equilibrium equation calculates.Note, suppose plane strain in tissue, to occur, utilize two-dimentional reconstructing method to calculate.In this simulation, adopt two kinds of models as organize models, each of these two kinds of models has the three dimensional structure as actual tissue, as shown in figure 21.Just, Figure 21 has shown two kinds of organize models's examples with three dimensional structure.Figure 21 (a) has shown that the inclusions that comprises the inclusions of serving as the mammary neoplasms model comprises model.Particularly, this inclusions comprises the external dimensions that model has 100mm * 100mm * 100mm, and comprises the hard inclusions that diameter is 20mm.Suppose that this inclusions has the elastic modelling quantity of 30kPa, and the material that surrounds this inclusions has the elastic modelling quantity of 10kPa.In mode, determine above-mentioned elastic modelling quantity according to the elastic modelling quantity of actual breast tissue as above-mentioned simulation.On the other hand, the material of inclusions and encirclement inclusions all shows near Incoercibility, and therefore, all is assumed that to have 0.49 identical Poisson's ratio.Figure 21 (b) has shown the layer structure model that is used to simulate such as the layer structure tissue of muscle.The layer structure model has the external dimensions of 100mm * 100mm * 100mm, and comprises the hard formation that thickness is 20mm.Suppose that hard formation has the elastic modelling quantity of 30kPa, and the material that surrounds this hard formation has the elastic modelling quantity of 10kPa.Notice that the layer structure model also has 0.49 even Poisson's ratio.

Then, comprise under the situation of model simulation on computers, the compression under the even external pressure of 100Pa that applies vertically, from the end face of model in the inclusions shown in Figure 21 (a).On the other hand, under the situation of the layer structure model shown in Figure 21 (b), simulation on computers, the compression under the even external pressure of 150Pa that applies vertically, from the end face of model.Note, to the compression under the different external pressures of above-mentioned two kinds of modelings, so that to the same strain of every kind of modeling about 1%.Then, for every kind of organize models, the RF signal of simulation under the situation of exerting pressure and not exerting pressure, and utilize expansion combination autocorrelation method to estimate axial stress distribution.Subsequently, utilize three-dimensional elastic modulus distribution reconstructing method, simulate the boundary condition of determining, estimate that elastic modelling quantity distributes according to the axial strain distribution of estimating and for compression.Equally, utilize two-dimentional reconstructing method,, estimate that elastic modelling quantity distributes, to compare according to identical axial strain distribution and same boundary conditions.At this, the parameter that is used for analog rf signal comprises: the mid frequency of 3.75MHz, the pulse width of 0.75mm; 2.0mm horizontal direction ultrasound wave beamwidth; 2.0mm slice direction ultrasound wave beamwidth; And the trace interval of 2.0mm.On the other hand, being used for utilizing expansion combination autocorrelation method to carry out parameters calculated comprises: the correlation window size of 3.2mm (axially) * 4.0mm (horizontal direction), and the hunting zone of 11.2mm (axially) * 14.0mm (horizontal direction).On the other hand, about utilizing the elastic modelling quantity distribution reconstructing method of three-dimensional finite element model, each organize models is made of 50000 rectangular parallelepiped protrusion part unit, and each rectangular parallelepiped protrusion part unit all has the size of 2mm (axially) * 2mm (horizontal direction) * 5mm (slice direction).

Figure 22 to Figure 25 has shown analog result.Figure 22 and Figure 23 have shown that inclusions comprises the estimated result of model.On the other hand, Figure 24 and Figure 25 have shown the estimated result of layer structure model.Notice that though can utilize three-dimensional reconstruction method to estimate that three-dimensional elastic modulus distributes, each width of cloth figure of Figure 24 and Figure 25 has shown the model viewgraph of cross-section of obtaining along the center line of model.Its reason is, utilizes two-dimentional reconstructing method can only estimate that the two-dimension elastic modulus distributes.Therefore, Figure 24 and Figure 25 have shown the viewgraph of cross-section of the estimated result of obtaining along the model center line, allow the user that they are compared thus.On the other hand, in following table, listed and utilized the three-dimensional reconstruction method that each organize models obtains and the estimated value of two-dimentional reconstructing method.

Estimated parameter comprises as used herein: the elastic modelling quantity error e in the zone of encirclement model core _SStandard deviation S D in the zone of encirclement model core _SAnd the contrast error e in the model nucleus _C, they are defined by following formula.

$e_S=\frac{\frac{1}{N_S}\underset{i}{\overset{N_S}{\mathrm{\&Sigma;}}}|{\hat{E}}_i-E_i|}{{\stackrel{\&OverBar;}{E}}_S}$

${\mathrm{SD}}_S=\frac{\sqrt{\frac{1}{N_S}\underset{i}{\overset{N_S}{\mathrm{\&Sigma;}}}{({\hat{E}}_i-{\stackrel{\&OverBar;}{E}}_S)}^2}}{{\stackrel{\&OverBar;}{E}}_S}$

$e_C=\frac{|({\hat{E}}_C-{\stackrel{\&OverBar;}{E}}_S)-(E_C-E_S)|}{E_C-E_S}\&CenterDot;\&CenterDot;\&CenterDot;\left(32\right)$

Note, in following formula, the S representative surround in the zone of core and, the elastic modelling quantity that the E^ representative is estimated, E represents actual elastic modulus, N _SUnit sum in the zone of representative encirclement core, E ^- _SThe meansigma methods of the elastic modelling quantity in the zone of representative encirclement model core, E^ _CRepresent the elastic modelling quantity in the estimated model nucleus, E _CActual elastic modulus in the representative model nucleus, and E _SActual elastic modulus in the zone of representative encirclement model core.

As finding out from above-mentioned analog result, confirmed, compare with the two-dimentional reconstructing method that forms based on the supposition of plane stress occurring in the tissue, the elastic modelling quantity distribution reconstructing method that proposes in the present embodiment, utilize three-dimensional finite element model has the advantage of estimating that more accurately elastic modelling quantity distributes., three-dimensional reconstruction method estimates accurately that the elastic modelling quantity that utilizes two-dimentional reconstructing method to estimate distributes and demonstrates than the little value of actual elastic modulus distribution though can distributing to elastic modelling quantity.Much less, its reason is, about two-dimentional reconstructing method, is restricted in the calculating of calculating on the direction of planar quadrature with two dimension.In this is estimated, clearly confirmed, utilize the elastic modelling quantity distribution reconstructing method of the three-dimensional computations corresponding to be applicable to actual tissue is analyzed with actual tissue.

Next, will the customized configuration of the compuscan of the above-mentioned elastic modelling quantity distribution reconstructing method that adopts above-mentioned expansion combination autocorrelation method and utilize three-dimensional finite element model be described.Figure 26 has shown the basic configuration of compuscan.Compuscan comprises: three-dimensional ultrasound scanner 281, and ultrasonic diagnostic equipment 282, personal computer 283, pulse electric machine controller 284, pulse motor 285, and piezometer 286 etc.Three-dimensional ultrasound scanner 281 has to tissue projection ultrasonic pulsative signal and receives from the function of the ultrasound echo signal of tissue reflection.Notice that this system adopts the elastic modelling quantity reconstructing method that utilizes three-dimensional finite element model, and therefore, this system needs in-house three-dimensional data.Therefore, this compuscan has three-dimensional ultrasound scanner utilized 281 and carries out the configuration of 3-D scanning.Ultrasonic diagnostic equipment 282 has the function of control three-dimensional ultrasound scanner 281 and demonstration real-time ultrasound B mode image, allows the position in the definite zone that will measure of user thus.Note, can be under situation about not changing, with the ultrasonic diagnostic equipment of routine as ultrasonic diagnostic equipment 282.This ultrasonic diagnostic equipment comprises digital device, and this digital device has the frame memory of the RF signal that is used for interim storage measurement.Personal computer 283 receives the RF signal of being measured by ultrasonic diagnostic equipment 282, carries out and handles (processing that proposes in the above-mentioned present embodiment) with estimation tissue elasticity characteristic, and show estimated result.Pulse motor 285 has the function that control puts on structural pressure.Pulse motor 285 is fixed on the top of fixed arm, and three-dimensional ultrasound scanner 281 is installed on the movable part of pulse motor 285.This system has following configuration, pulse electric machine controller 284 control impuls motor 285 wherein, put on structural pressure to regulate the vertical direction position of ultrasound scanner 281, to allow the user to regulate thus, this pressure impels tissue that the mild compression of a few percent takes place accurately.Piezometer 286 has the function of measuring the pressure on the body surface; Pressure is used as the required boundary condition of reconstruct elastic modelling quantity distribution.Notice that piezometer 286 is between ultrasound scanner 281 and body surface.At this, suppose that the tissue compression by ultrasound scanner 281 has caused the uniform pressure on the body surface to distribute, the pressure that then will utilize piezometer 286 to measure is used as the pressure on the body surface.

Figure 27 has shown the customized configuration of the ultrasound scanner 281 that comprises in the compuscan.Three-dimensional ultrasound scanner 281 does not have the two-dimensional array configuration of wherein having disposed ultrasonic transducer two-dimensionally, but has the 3-D scanning configuration, in this 3-D scanning configuration, carries out the mechanical scanning of two-dimentional convex scanheads in water, on slice direction.

Compuscan shown in Figure 26 mainly is designed to the diagnosis of breast carcinoma, and therefore, determines the characterisitic parameter of ultrasound scanner at the diagnosis of breast carcinoma.Particularly, the key property parameter of ultrasound scanner comprises as used herein: the ultrasound wave mid frequency of 7.5MHz, the sample frequency of 30MHz, 71 number of scanning lines, 44 frame number, 30 ° transducer scan angle, and the persistent period in 0.5 second 3-D scanning cycle.Note,, on slice direction, carry out the scanning of convex probe by above-mentioned transducer scan angle.In addition, above-mentioned frame number is illustrated in the quantity of the scanogram (frame) that obtains in the scan period of convex probe.On the other hand, comprise by the characteristic of utilizing line (wire) in the water to carry out the ultrasonic pulse that actual measurement obtains: the about pulse width of 0.5mm; The horizontal direction beamwidth of about 1.5mm; And the slice direction beamwidth of about 2.6mm.

Next, carry out the operation that elastic characteristic is measured with describing compuscan shown in Figure 26.At first, the user moves the spatial digitizer 281 that is installed on the arm, so that ultrasound scanner 281 is positioned at the expectation part that will measure, and the real-time B mode image that monitoring simultaneously obtains by ultrasonic diagnostic equipment 282.Note, when ultrasound scanner 281 is positioned, do not carry out the 3-D scanning (that is, not carrying out the mechanical scanning of convex probe) of ultrasound scanner 281, be displayed on the ultrasonic diagnostic equipment 282 and have only with the corresponding B mode image of the median plane of scanner.Its reason is, about ultrasonic diagnostic equipment 282 as used herein, and can not be when carrying out 3-D scanning, show real-time B mode image.Much less, a kind of like this ultrasonic diagnostic equipment of can sampling, this equipment has can be when carrying out 3-D scanning, show the function of real-time B mode image, allows the user when carrying out 3-D scanning, three-dimensional ultrasound scanner 281 is positioned thus.After ultrasound scanner 281 was positioned at the desired locations that will measure, the movable part of user's fixed arm was so that fixing ultrasound scanner 281.Then, system is not measured there being the three-dimensional RF signal under the situation that tissue is exerted pressure.Notice that in case the user pushes the button that is used for 3-D scanning, 3-D scanning just automatically carries out.At this, the persistent period in a 3-D scanning cycle only is 0.5 second.Measure R F data under the situation of not exerting pressure are stored in the interior frame memory of ultrasonic device.Subsequently, in case the user push be used for operating impulse motor controller 284 button once, with the compression of control tissue, the pulse motor 285 that is fixed on the arm just makes ultrasound scanner 281 move preset distance, thus by ultrasound scanner 281 compress tissues.Subsequently, pulse motor 285 stops, and keeps the compression of tissue simultaneously.Under this state, the user pushes the button that is used for 3-D scanning once more, measures the RF data under the situation of exerting pressure to tissue thus.Note, in mode, in the frame memory that the RF data storage under the situation of exerting pressure to tissue is comprised in ultrasonic device 282 as the RF data under the situation of not exerting pressure.In addition, utilize the piezometer of the end that is installed in ultrasound scanner 281, measure and put on structural pressure.Then, measure to finish, and discharge and put on structural pressure, after this discharge object.

After discharging object, system visits the frame memory that comprises in the ultrasonic diagnostic equipments 282 by personal computer 283, and is stored on the hard disk that comprises in the personal computer 283 in the RF data under the situation of exerting pressure and not exerting pressure to tissue.Carry out this processing and be because, the frame memory that comprises in the ultrasonic device has the function of temporary storaging data, promptly this frame memory only is useful on the capacity of the data of a measuring period of storage.Personal computer 283 is carried out such program, the elastic modelling quantity distribution reconstructing method of this program utilization expansion combination autocorrelation method and employing three-dimensional finite element model calculates, so that according to the RF data under the situation of exerting pressure and not exerting pressure, estimate that stress distribution and elastic modelling quantity distribute to tissue.Then, personal computer 283 shows B mode image, strain pattern and elastic modelling quantity image by carrying out display routine on monitor.

According to of the present invention, utilize the compuscan of distortion distribution display method and elastic modelling quantity distribution display packing to have the following advantages: can be only according to the ultrasonic beam direction stress distribution of (axially), come the reconstruct elastic modelling quantity to distribute; And can with the horizontal direction displacement irrespectively, estimate Displacements Distribution.

Note,, the invention is not restricted to above-mentioned configuration, but represent any parameter of the relation between the undulatory property that comprises amplitude, wave height and wave number to be used though described the configuration of adopting envelope signal.

Claims (8)

1, a kind of compuscan comprises:

Ultrasound probe, and carry out the transmitting-receiving of ultrasonic signal between the object;

Correlation calculation apparatus, allow the envelope of a pair of reflection echo signal before and after the compression of the described object that received by described ultrasound probe relatively move successively on the direction of principal axis of ultrasonic beam and set at interval, the correlation coefficient of obtaining the envelope before and after the compression becomes maximum relevant maximum position and the phase contrast under this relevant maximum position; With

The strain accountant according to described relevant maximum position of being obtained by described correlation calculation apparatus and described phase contrast, is obtained the stress distribution of the tissue of the displacement of tissue of the described object of following described compression and described object;

Described correlation calculation apparatus has:

First orthogonal demodulation circuit carries out orthogonal detection to the reflection echo signal before compressing;

Second orthogonal demodulation circuit carries out orthogonal detection to the reflection echo signal after the compression; With

A plurality of delay circuits postpone setting-up time successively to the output of described second orthogonal demodulation circuit;

According to the output of described first orthogonal demodulation circuit and described a plurality of delay circuits, obtain described correlation coefficient and described phase contrast.

2, compuscan according to claim 1 is characterized in that,

Described correlation calculation apparatus has:

The first correlation computations circuit according to the output of described first orthogonal demodulation circuit, is obtained the autocorrelation value of the envelope of the preceding reflection echo signal of described compression;

The second correlation computations circuit according to the output of described second orthogonal demodulation circuit, is obtained the autocorrelation value of the envelope of the reflection echo signal after the described compression;

Third phase closes counting circuit, according to the output of described first orthogonal demodulation circuit and described second orthogonal demodulation circuit, obtains the cross correlation value of the envelope of a pair of reflection echo signal before and after the described compression;

The correlation coefficient counting circuit, close the correlation that counting circuit is obtained respectively according to described first orthogonal demodulation circuit and described second orthogonal demodulation circuit and described third phase, obtain the described correlation coefficient of the envelope of a pair of reflection echo signal before and after the described compression; With

The phase difference calculating circuit according to the output of described first orthogonal demodulation circuit and described second orthogonal demodulation circuit, is obtained the described phase contrast of a pair of reflection echo signal of described compression front and back.

3, compuscan according to claim 1 is characterized in that,

Described a plurality of delay circuit postpones the output maximum of described second orthogonal demodulation circuit interval of 1/2nd wavelength of described ultrasonic signal successively.

4, compuscan according to claim 1 is characterized in that,

Also comprise the elastic modelling quantity accountant, the tissue of described object is divided into the unit of limited quantity, set up two-dimensional finite element model at least; According to the information and the described stress distribution that are used to set up described model, come the calculating elastic modulus to distribute.

5, compuscan according to claim 1 is characterized in that,

Also comprise the elastic modelling quantity accountant, suppose the isotropic elastomer of being organized as of described object and near Incoercibility, the tissue of described object is divided into the Nogata body unit of limited quantity, set up three-dimensional finite element model, and supposition elastic modelling quantity, stress and strain in each described unit is all even, described stress distribution is used for the elastic equation formula, comes the calculating elastic modulus to distribute.

6, a kind of distortion distribution display method comprises:

First step, input is by contacting obtained ultrasound echo signal with ultrasound probe with the tissue of object;

Second step, allow the envelope of a pair of reflection echo signal before and after the compression of the described object that received by described ultrasound probe relatively move successively on the direction of principal axis of ultrasonic beam and set at interval, the correlation coefficient of obtaining the envelope before and after the compression becomes maximum relevant maximum position and the phase contrast under this relevant maximum position;

Third step according to described relevant maximum position of obtaining in described second step and described phase contrast, is obtained the stress distribution of the tissue of the displacement of tissue of the described object of following described compression and described object; With

The 4th step shows described stress distribution;

Described second step has:

The first orthogonal detection step is carried out orthogonal detection to the reflection echo signal before compressing;

The second orthogonal detection step is carried out orthogonal detection to the reflection echo signal after the compression; With

Postpone step, the output of the described second orthogonal detection step is postponed setting-up time successively;

According to the output of described first orthogonal detection step and described delay step, obtain described correlation coefficient and described phase contrast.

7, a kind of elastic modelling quantity distribution display packing comprises:

First step, input is by contacting obtained ultrasound echo signal with ultrasound probe with the tissue of object;

Second step, allow the envelope of a pair of reflection echo signal before and after the compression of the described object that received by described ultrasound probe relatively move successively on the direction of principal axis of ultrasonic beam and set at interval, the correlation coefficient of obtaining the envelope before and after the compression becomes maximum relevant maximum position and the phase contrast under this relevant maximum position;

Third step according to described relevant maximum position of obtaining in described second step and described phase contrast, is obtained the stress distribution of the tissue of the displacement of tissue of the described object of following described compression and described object;

The 4th step is divided into the unit of limited quantity with the tissue of described object, sets up two-dimensional finite element model at least, according to the information and the described stress distribution that are used to set up described model, comes the calculating elastic modulus to distribute; With

The 5th step shows that described elastic modelling quantity distributes;

Described second step has:

The first orthogonal detection step is carried out orthogonal detection to the reflection echo signal before compressing;

The second orthogonal detection step is carried out orthogonal detection to the reflection echo signal after the compression; With

Postpone step, the output of the described second orthogonal detection step is postponed setting-up time successively;

According to the output of described first orthogonal detection step and described delay step, obtain described correlation coefficient and described phase contrast.

8, elastic modelling quantity distribution display packing according to claim 7 is characterized in that,

Described the 4th step, suppose the isotropic elastomer of being organized as of described object and near Incoercibility, the tissue of described object is divided into the Nogata body unit of limited quantity, set up three-dimensional finite element model, and supposition elastic modelling quantity, stress and strain in each described unit is all even, described stress distribution is used for the elastic equation formula, comes the calculating elastic modulus to distribute.

Images