CN101634618A - Three-dimensional elastic modulus imaging method - Google Patents

Abstract

The invention discloses a three-dimensional elastic modulus imaging method which comprises the following steps: (A) establishing a three-dimensional finite element model of an initial state of an object; (B) applying a system load to the object; (C) measuring the three-dimensional motion information of a finite element node of the object, which is applied with the system load, wherein the three-dimensional motion information comprises displacement and speed; (D) obtaining a continuous state estimation equation by an H infinite filtering estimation method based on a state space, estimating the three-dimensional distribution of the elastic modulus of the object by the continuous state estimation equation and realizing the three-dimensional elastic modulus imaging of the object. The invention has the advantage that the material characteristic of the object can be robustly acquired on the basis of the H infinite filtering estimation method of the state space by combining a finite element model under the condition of knowing the three-dimensional motion information of the node of the object.

Description

A kind of three-dimensional elastic modulus imaging method

Technical field

The present invention relates to a kind of three-dimensional elastic modulus imaging technology, specifically a kind of deformation inverting elastic modulus that under the effect of power, produces according to object, thus realize imaging method.

Background technology

Modern imaging technique can make people observe naked eyes the internal structure of body that can't see and the situation of motion change.Early stage X-ray can only simply be expressed the transmissivity of different material, nowadays, because the continuous development of computer graphic image technology, ultrasound wave, MRI emerging imaging technique means such as (magnetic resonance imagings) not only can show the inner structure of object more clearly, but also can obtain the deformational displacement and the velocity information of interior of articles particle dynamically.

Many times, we observe the structure of interior of articles and distorted movement information is in order to understand the variation of its character, such as the physical construction fatigue strength survey, structural integrity analyzes or the like.Yet this is not a directly reflection.The variation of physical property has at first caused the variation of material parameter, and when being subjected to the outside time spent of doing, this variation just can show by the unusual of internal structure of body and motion.If these information of Direct observation can not be understood the variable condition of physical property very intuitively.Therefore, the material parameter that obtains interior of articles distributes and also directly shows, and all there is very important meaning in many fields.

In each parameter of material properties, most importantly elastic modulus is called Young modulus again, and the ratio of the stress and strain that produces behind the expression object suffered external action has reflected the rigidity of material, and it is the material parameter that will estimate in this method.

Because the restriction of technical merit also lacks the directly method of the elastic modulus of non-invasive measurement at present.Present research method generally is that the distorted movement ten-four of suppose object (is implanted the seed acquisition or measured distorted movement information by the computer graphics method such as relying on, be that shift value or velocity amplitude are known), estimate material parameter such as elastic modulus from these known information according to constitutive relation then.Actual imaging detection means can be subjected to The noise inevitably in observation process.The reflection of changes in material on motion state was just not obvious originally, added interference of noise, and it is just difficult more to obtain material parameter from these information, can only estimate by certain mathematical method.

The object of reality all is three-dimensional, at present full 3D imaging and motion detection technology, image three-dimensional analytical technology are also in continuous development, utilize the more approaching reality of elastic modulus of the whole object of three-dimensional motion information direct estimation, and be familiar with for the people more intuitively, judge thereby make more accurately.

Summary of the invention

The object of the present invention is to provide a kind of three-dimensional elastic modulus imaging method, may further comprise the steps:

(A) set up the three-dimensional finite element model of object original state;

The foundation of three-dimensional finite element model can realize according to existing dimensional finite element method.

(B) to object application system load;

Described system load can be the acting force of known dimensions, also can be measure be applied to variation of force value on the object;

(C) three-dimensional motion information of finite element node after application system load of measuring object, described three-dimensional motion information comprises displacement, speed;

Because the intrinsic property of different measuring method of the prior art make the finite element node three-dimensional motion information of some object to measure, and some finite element node three-dimensional motion information can't be measured, and need estimate;

(D) use obtains following continuous state estimate equation (20)～(23) based on the H ∞ Filtering Estimation method of state space, utilizes equation (20)～(23) to estimate the distributed in three dimensions of object elastic modulus, realizes the three-dimensional elastic modulus imaging of object.

$\left[\begin{array}{c}\widehat{\stackrel{\·}{x}}\\ \widehat{\stackrel{\·}{\mathrm{\θ}}}\end{array}\right]=\left[\begin{array}{cc}F& A_C\left(x\left(t\right)\right)\\ 0& 0\end{array}\right]\left[\begin{array}{c}\hat{x}\left(t\right)\\ \widehat{\mathrm{\θ}}\end{array}\right]+\left[\begin{array}{c}{B}_{C}w\left(t\right)\\ 0\end{array}\right]+\mathrm{\Σ}{\left(t\right)}^{-1}{D}^T(y-D\hat{x}\left(t\right))---\left(20\right)$

$\stackrel{\·}{\mathrm{\Σ}}=-\mathrm{\Σ}\left(t\right)\left[\begin{array}{cc}F& A_C\left(x\left(t\right)\right)\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}{F}^T& 0\\ A_C^T\left(x\left(t\right)\right)& 0\end{array}\right]\mathrm{\Σ}\left(t\right)+\left[\begin{array}{cc}I& 0\\ 0& -{\mathrm{\γ}}^{-2}Q\end{array}\right]-\mathrm{\Σ}\left(t\right)\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]\mathrm{\Σ}\left(t\right)---\left(21\right)$

$\left[\begin{array}{c}\hat{x}(t+\mathrm{\Δt})\\ \widehat{\mathrm{\θ}}\end{array}\right]=\left[\begin{array}{c}\hat{x}\left(t\right)\\ \widehat{\mathrm{\θ}}\end{array}\right]+\left[\begin{array}{c}\widehat{\stackrel{\·}{x}}\\ \widehat{\stackrel{\·}{\mathrm{\θ}}}\end{array}\right]\mathrm{\Δt}---\left(22\right)$

$\mathrm{\Σ}(t+\mathrm{\Δt})=\mathrm{\Σ}\left(t\right)+\stackrel{\·}{\mathrm{\Σ}}\mathrm{\Δt}---\left(23\right)$

Wherein:

$F=\left[\begin{array}{cc}0& I\\ 0& -\mathrm{\αI}\end{array}\right];$

$A_C=\left[\begin{array}{cc}0& 0\\ -M^{-1}{G}_1& -\mathrm{\β}{M}^{-1}{G}_2\end{array}\right];$

$B_C=\left[\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right];$

$\hat{x}\left(t\right)={\left[\begin{array}{cc}\hat{U}\left(t\right)& \widehat{\stackrel{\·}{U}}\left(t\right)\end{array}\right]}^T;$

U is the motion vector of the finite element node of object, its single order differential

And second-order differential Represent velocity vector and vector acceleration respectively;

I is and the corresponding unit matrix of U dimension;

α and β are the Rayleigh damping constant, desirable 0.01～0.05;

M is the finite element mass matrix of object, is the function of subject material density;

G ₁And G ₂Be the transition matrix in the computation process, building method such as step (a)～(e):

(a) two N * N of initialization _eEmpty matrix G ₁And G ₂, N is the system node variable number, in three dimensions, each node has the degree of freedom of three directions, i.e. the total node number of system * 3, N _eFor the total unit number of system (referring to the model among Fig. 1, the unit is exactly one of them tetrahedron, node is exactly tetrahedral four end points, a node may be by a plurality of units shareds in system, the system unit number is exactly the number that Object Segmentation becomes tetrahedron element, determines when setting up model meshes); The empty matrix K of N * N of initialization (line number of representing matrix and columns) _g', as the interim stiffness matrix of the system that does not contain elastic modulus E;

(b) for a unit that is numbered j, structure does not comprise this unitary elasticity modulus E _jStiffness matrix K _j';

The coding rule of unit and node during (c) according to system's generating mesh is K _j' in the interim stiffness matrix K of element insertion system _g' in the corresponding position, and make subscript in its unit into corresponding system subscript;

(d) _KG ' and the vector that U multiplies each other and obtains promptly are matrix G ₁J row; K _g' with The vector that multiplies each other and obtain promptly is matrix G ₂J row;

(e) turn back to (b) step up to all node traversals, finish G ₁And G ₂Structure;

T is a time variable;

Δ t is interval estimated time;

Be the elastic modulus vector estimated value of object, the θ initial value preestablishes according to the experimental knowledge of object, and general requirement is approaching with the true elastic modulus of object;

W (t) is a control vector, w (t)=[0 R] ^T, R is the system load vector;

The three-dimensional motion information that y (t) obtains for measuring object in the step (C);

D is for measuring matrix, the line number of D is the finite element node number that three-dimensional motion information can measure, the columns of D is the dimension of x (t), in every row of D, if the three-dimensional motion information of the finite element node of index position correspondence can measure, the element of this index position is made as 1, otherwise element is made as 0;

Be the derivative of x (t), be illustrated in the variable quantity of concrete t x (t) constantly;

Modified value for θ;

γ influences the estimated value degree of accuracy for the binding occurrence of the quality of estimation, generally can get 1;

$\mathrm{\Σ}=\left[\begin{array}{cc}{\mathrm{\Σ}}_1& {\mathrm{\Σ}}_2\\ {\mathrm{\Σ}}_2^T& {\mathrm{\Σ}}_3\end{array}\right],$ ∑ wherein ₁, ∑ ₂, ∑ ₃The initial value ∑ ₁(0), ∑ ₂(0), ∑ ₃(0) is respectively ∑ ₁(0)=Q ₁, ∑ ₂(0)=0, ∑ ₃(0)=Q ₀

Be the derivative of ∑, be illustrated in the variable quantity of concrete t ∑ constantly;

Q, Q ₀, Q ₁Be weight matrix, its value is set as required, characterizes the influence degree of different factors to the distributed in three dimensions of object elastic modulus;

Subscript ^ in each parameter is the estimated value of this parameter, subscript ^TRepresent this transpose of a matrix matrix.

Advantage of the present invention is based on the H ∞ Filtering Estimation method of state space, in conjunction with finite element model, is implemented under the situation of known object node three-dimensional motion information, obtains the material behavior of object robust.

Description of drawings

Fig. 1 is the realistic model that the present invention is used for verification algorithm;

Fig. 2 is the elastic modulus distribution schematic diagram that the realistic model among Fig. 1 uses this method to reconstruct on displacement field basis shown in Figure 2.

Embodiment

Below in conjunction with accompanying drawing the specific embodiment of the present invention is elaborated, elastic modulus imaging method of the present invention may further comprise the steps:

(A) set up the three-dimensional finite element model of object original state;

(B) apply the acting force of known dimensions to object, perhaps measure the variation of force value that has been applied on the object, as system load;

(C) the measuring object node is after application system load, or the three-dimensional motion information of existing system load after changing, and described three-dimensional motion information comprises displacement, speed;

(D) according to the experimental knowledge about object, the value of the priori that the true elastic modulus of setting and imaging object is approaching uses Finite Element Method to set up the mechanical model of object, sets up the system dynamics equation; Set up the state equation group according to the system dynamics equation, use H ∞ Filtering Estimation method, estimate the distributed in three dimensions of object elastic modulus, be embodied as picture based on state space.

In three-dimensional problem, common cell type is tetrahedron element and hexahedral element.Because tetrahedron element has better adaptability for the gridding of complicated shape, the grid dividing algorithm is comparative maturity also, and calculated amount is little, therefore adopts this cell configuration here.

Displacement components u in the tetrahedron element can be only with the displacement components u of its four nodes _eRepresent (subscript e all represents a variable in the unit):

u＝Nu _e (1)

Wherein N is the tetrahedron element shape function matrix, concrete form:

$N=\left[\begin{array}{cccccccccccc}{N}_i& 0& 0& N_j& 0& 0& N_k& 0& 0& N_m& 0& 0\\ 0& N_i& 0& 0& N_j& 0& 0& N_k& 0& 0& N_m& 0\\ 0& 0& N_i& 0& 0& N_j& 0& 0& N_k& 0& 0& N_m\end{array}\right]---\left(2\right)$

N wherein _i, N _j, N _k, N _mBe respectively four node corresponding shape of tetrahedron element function, can calculate by following formula:

$N_p=\frac{1}{6V^e}(a_{p,1}+a_{p,2}x+a_{p,3}y+a_{p,4}z)(p=i,j,k,m)---\left(3\right)$

The matrix that note is made of four apex coordinates of tetrahedron element:

$\mathrm{\Λ}=\left[\begin{array}{cccc}1& x_i& y_i& z_i\\ 1& x_j& y_j& z_j\\ 1& x_k& y_k& z_k\\ 1& x_m& y_m& z_m\end{array}\right]---\left(4\right)$

A then _{P, n}For the p of matrix Λ is capable, the algebraic complement of n column element, unit volume V _e=| Λ |/6.

Consider that from the angle of calculating suppose object is the linear elastomer of isotropic, its stress vector σ with answer variable vector ε linear:

σ＝Dε (5)

Wherein, material matrix D:

$D=\frac{E}{(1+v)(1-2v)}\left[\begin{array}{cccccc}1-v& v& v& 0& 0& 0\\ v& 1-v& v& 0& 0& 0\\ v& v& 1-v& 0& 0& 0\\ 0& 0& 0& 1-2v& 0& 0\\ 0& 0& 0& 0& 1-2v& 0\\ 0& 0& 0& 0& 0& 1-2v\end{array}\right]={\mathrm{ED}}^{\′}---\left(6\right)$

E is an elastic modulus, and the expression object produces the difficulty of deformation; V is a Poisson ratio, the ratio of expression transverse strain and longitudinal strain.These two material characters that all belong to object, in the reality variation of Poisson ratio less, can think constant.Elastic modulus is the material parameter that will estimate herein.

The strain in the space and the relation of displacement are by the three-dimensional geometry The Representation Equation:

$\mathrm{\ϵ}=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂y}& 0\\ 0& 0& \frac{\∂}{\∂z}\\ \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\\ 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\end{array}\right]u=\mathrm{Lu}---\left(7\right)$

Can obtain by formula (1) and (7):

ε＝LNu _e＝Bu _e (8)

The strain that this shows each point in the tetrahedron element is the same, is normal strain unit, and the strain size determines that by the displacement of four nodes B is an element strain matrix.The formula that is calculated as follows of tetrahedron element stiffness matrix:

K _e＝B ^TDBV _e (9)

The mass matrix of tetrahedron element adopts the form of lumped mass matrix, thinks that promptly mass concentration is on the node of unit:

$M_e=\frac{\mathrm{\ρ}{V}_e}{4}I---\left(10\right)$

Wherein ρ is the material average density of object, and I is 12 * 12 unit matrix.

According to physics principle, derive the matrix form of finite element model nodal displacement system dynamics equation:

$M\stackrel{\·\·}{U}+C\stackrel{\·}{U}+\mathrm{KU}=R---\left(11\right)$

All nodes of whole object are formed a system, and M represents the mass matrix of system in the formula (4), is the function of subject material density; C is the damping matrix of system; K is the stiffness matrix of system, is the strain tensor of object and the relation function between the stress tensor; R is the load vector (being system load) of all nodes of system, the acting force that reflection object is suffered; U is the motion vector of all nodes of system, its single order differential

And second-order differential

Represent velocity vector and vector acceleration respectively; C is relevant with frequency, is assumed to be Rayleigh damping, and then C satisfies formula C=α M+ β K (α, β get 0.01-0.05); C expression formula substitution (11) is obtained formula (12):

$M\stackrel{\·\·}{U}+\mathrm{\αM}\stackrel{\·}{U}+\mathrm{\βK}\stackrel{\·}{U}+\mathrm{KU}=R---\left(12\right)$

Rebuild elastic modulus, promptly, utilize existing displacement information and boundary force information to try to achieve the distribution of elastic modulus E by recombination system dynamic equation (12).Therefore the left side of equation (12) need be transformed into the function of E vector, i.e. KU=G ₁E, $K\stackrel{\·}{U}=G_{2}E,$ By the so as can be known conversion of equation (6) is feasible, obtains thus:

$M\stackrel{\·\·}{U}+\mathrm{\αM}\stackrel{\·\·}{U}+\mathrm{\β}{G}_{2}E+G_{1}E=R---\left(13\right)$

According to element stiffness matrix K in the Finite Element Method _eDefinition, K _eCan be expressed as containing unitary elasticity modulus E _e(subscript e represents the unit to the form of item, K _eAnd E _eBe respectively the stiffness matrix and the elastic modulus of a unit, do not descend target just to represent the stiffness matrix and the elastic modulus of total system):

$K_e=E_e\underset{\mathrm{\Ωe}}{\∫}{B}_e^T{D_e}^{\′}{B}_{e}d{\mathrm{\Ω}}_e=E_e{K_e}^{\′}---\left(14\right)$

Thus, the method according to the stiffness matrix structure obtains constructing G ₁And G ₂Method:

(a) two N * N of initialization _eEmpty matrix G ₁And G ₂, N is the system node variable number, in three dimensions, each node has the degree of freedom of three directions, i.e. the total node number of system * 3, N _eFor the total unit number of system (is seen the model in the accompanying drawing, the unit is exactly one of them tetrahedron, node is exactly tetrahedral four end points, a node may be by a plurality of units shareds in system, the system unit number is exactly the number that Object Segmentation becomes tetrahedron element, determines when setting up model meshes); The empty matrix K of N * N of initialization (line number of representing matrix and columns) _g', as the interim stiffness matrix of the system that does not contain elastic modulus E;

(b), utilize equation (14) structure not comprise the stiffness matrix K of this unitary elasticity modulus Ej for a unit that is numbered j _j';

The coding rule of unit and node during (c) according to system's generating mesh is K _j' in the interim stiffness matrix K of element insertion system _g' in the corresponding position, and make subscript in its unit into corresponding system subscript.

(d) K _g' with the vector that U multiplies each other and obtains, promptly be matrix G ₁J row; K _g' with

The vector that multiplies each other and obtain promptly is matrix G ₂J row;

(e) turn back to (b) step up to all node traversals.

Construct N * N thus _eMatrix G ₁, G ₂

In order to solve the problem that has noise in the amount of exercise process of measuring, adopt and realize estimation procedure based on the noise disturbance total state information approach (NPFSI) of H ∞.

In order to use this Filtering Estimation method, at first set up system state equation.If state vector is $x\left(t\right)={\left[\begin{array}{cc}U\left(t\right)& \stackrel{\·}{U}\left(t\right)\end{array}\right]}^T,$ Control vector is w (t)=[0 R] ^T, with formula (12) distortion, set up linear stochaastic system state-space expression continuous time:

$\stackrel{\·}{x}\left(t\right)=\left[\begin{array}{cc}0& I\\ -M^{-1}K& -M^{-1}C\end{array}\right]x\left(t\right)+\left[\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right]w\left(t\right)+v\left(t\right)---\left(15\right)$

y(t)＝Dx(t)+e(t) (16)

D is for measuring matrix, the line number of D is the finite element node number that three-dimensional motion information can measure, the columns of D is the dimension of x (t), in every row of D, if the three-dimensional motion information of the finite element node of index position correspondence can measure, the element of this index position is made as 1, otherwise element is made as 0.V (t) represents process noise, and e (t) represents observation noise.

Because being included in the parameter matrix of elastic modulus E, for the purpose that realizes that we estimate, we establish the material vector and are θ=[E], again former movement-state x and material parameter θ are formed a new state vector, equation (15) (16) is out of shape, obtain the state equation new as next group, we suppose that material parameter θ does not change in time here:

$\left[\begin{array}{c}\widehat{\stackrel{\·}{x}}\\ \widehat{\stackrel{\·}{\mathrm{\θ}}}\end{array}\right]=\left[\begin{array}{c}\mathrm{Fx}\left(t\right)+A_C\left(x\left(t\right)\right)\mathrm{\θ}+B_{C}w\left(t\right)\\ 0\end{array}\right]+\left[\begin{array}{c}v\left(t\right)\\ 0\end{array}\right]---\left(17\right)$

$y=\left[\begin{array}{cc}D& 0\end{array}\right]\left[\begin{array}{c}x\\ \mathrm{\θ}\end{array}\right]+e\left(t\right)---\left(18\right)$

Wherein $F=\left[\begin{array}{cc}0& I\\ 0& -\mathrm{\αI}\end{array}\right],$ $A_C=\left[\begin{array}{cc}0& 0\\ -M^{-1}{G}_1& -\mathrm{\β}{M}^{-1}{G}_2\end{array}\right],$ $B_C=\left[\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right].$

H ∞ optimal estimation problem is by the design point estimator, makes process noise and observation noise to any situation, and evaluated error is all less, promptly realizes the optimal estimation under the unknown noise.The cost function of said process is as follows:

$\sup \frac{{\left|\right|\mathrm{\&theta;}-\widehat{\mathrm{\&theta;}}\left|\right|}_Q^2}{{\left|\right|v\left|\right|}^2+{\left|\right|e\left|\right|}^2+{|\mathrm{\&theta;}-\widehat{\mathrm{\&theta;}}\left(0\right)|}_{Q_0}^2+{|x-\hat{x}|}_{Q_1}^2}<{\mathrm{\&gamma;}}^2---\left(19\right)$

H ∞ wave filter is actually the problem of game between external disturbance and the estimator, i.e. a minimum maximization problems.Obtain as next group continuous state estimate equation by state equation:

$\left[\begin{array}{c}\widehat{\stackrel{\&CenterDot;}{x}}\\ \widehat{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}\end{array}\right]=\left[\begin{array}{cc}F& A_C\left(x\left(t\right)\right)\\ 0& 0\end{array}\right]\left[\begin{array}{c}\hat{x}\left(t\right)\\ \widehat{\mathrm{\&theta;}}\end{array}\right]+\left[\begin{array}{c}{B}_{C}w\left(t\right)\\ 0\end{array}\right]+\mathrm{\&Sigma;}{\left(t\right)}^{-1}{D}^T(y-D\hat{x}\left(t\right))---\left(20\right)$

$\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}=-\mathrm{\&Sigma;}\left(t\right)\left[\begin{array}{cc}F& A_C\left(x\left(t\right)\right)\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}{F}^T& 0\\ A_C^T\left(x\left(t\right)\right)& 0\end{array}\right]\mathrm{\&Sigma;}\left(t\right)+\left[\begin{array}{cc}I& 0\\ 0& -{\mathrm{\&gamma;}}^{-2}Q\end{array}\right]-\mathrm{\&Sigma;}\left(t\right)\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]\mathrm{\&Sigma;}\left(t\right)---\left(21\right)$

$\left[\begin{array}{c}\hat{x}(t+\mathrm{\&Delta;t})\\ \widehat{\mathrm{\&theta;}}\end{array}\right]=\left[\begin{array}{c}\hat{x}\left(t\right)\\ \widehat{\mathrm{\&theta;}}\end{array}\right]+\left[\begin{array}{c}\widehat{\stackrel{\&CenterDot;}{x}}\\ \widehat{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}\end{array}\right]\mathrm{\&Delta;t}---\left(22\right)$

$\mathrm{\&Sigma;}(t+\mathrm{\&Delta;t})=\mathrm{\&Sigma;}\left(t\right)+\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}\mathrm{\&Delta;t}---\left(23\right)$

For algorithm for estimating is got a desired effect, before beginning estimation, must set suitable noise Reduction Level, weight matrix Q, Q earlier ₀, Q ₁With the γ value.Usually be difficult to the optimum γ of decision ^*, it may be infinite, means the noise controlling level that can't reach expectation.Yet, if in fact Q is chosen, γ ^*Can obtain by analyzing.At first with the ∑ partitioning of matrix:

$\mathrm{\&Sigma;}=\left[\begin{array}{cc}{\mathrm{\&Sigma;}}_1& {\mathrm{\&Sigma;}}_2\\ {\mathrm{\&Sigma;}}_2^T& {\mathrm{\&Sigma;}}_4\end{array}\right]---\left(24\right)$

With (22) substitution (21) formula, obtain the time diffusion renewal equation of partitioned matrix:

${\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}}_1=I-{\mathrm{\&Sigma;}}_1^2-{\mathrm{\&Sigma;}}_{1}F-F^T{\mathrm{\&Sigma;}}_1,$ ∑ ₁(0)＝Q ₁ (25)

${\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}}_2=-{\mathrm{\&Sigma;}}_1{A}_c-F^T{\mathrm{\&Sigma;}}_2-{\mathrm{\&Sigma;}}_1{\mathrm{\&Sigma;}}_2,$ ∑ ₂(0)＝0 (26)

${\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}}_3=-{\mathrm{\&Sigma;}}_2^T{A}_C-A_C^T{\mathrm{\&Sigma;}}_2-{\mathrm{\&gamma;}}^{-2}Q-{\mathrm{\&Sigma;}}_2^T{\mathrm{\&Sigma;}}_2,$ ∑ ₃(0)＝Q ₀ (27)

According to the Schur principle, guarantee ∑＞0, have only and work as:

∑ ₁(t)＞0 (26)

$\mathrm{\&Pi;}={\mathrm{\&Sigma;}}_3-{\mathrm{\&Sigma;}}_2^T{\mathrm{\&Sigma;}}_1^{-1}{\mathrm{\&Sigma;}}_2>0---\left(28\right)$

∏ (t) is separating of following matrix differential equation:

$\stackrel{\&CenterDot;}{\mathrm{\&Pi;}}={\mathrm{\&Sigma;}}_2^T{\mathrm{\&Sigma;}}_1^{-1}{\mathrm{\&Sigma;}}_2-{\mathrm{\&gamma;}}^{2}Q,$ ∏(0)＝Q ₀ (29)

Because the parameter that This document assumes that will be estimated is a time constant, so order $Q={\mathrm{\&Sigma;}}_2^T{\mathrm{\&Sigma;}}_1^{-1}{\mathrm{\&Sigma;}}_2,$ γ＝1，∑ ₁(0)＝I。

In computation process, the initial value of elastic modulus E can be according to the experimental knowledge to cardiac muscle, the value of the some priori that the true elastic modulus of setting and imaging object is approaching is utilized equation (20), (21) (or (25-27)) then, (22), (23) default elastic modulus E is carried out recurrence and estimate to revise,, promptly obtain the distributed intelligence of elastic modulus E up to obtaining the convergent optimum solution, re-use computing machine 3D technology and show, realize the elasticity modulus of materials imaging.

The emulation experiment model is three-dimensional semielliptical shell, high 100mm, and internal diameter 35mm, external diameter 50mm, its Young modulus distributes as shown in Figure 1.Light-colored part is represented softer part, and Young modulus is 75kPa.Dark part representative is than hard portion, and Young modulus is 105kPa.The upper surface of model is fixed, and inwall applies the pressure in the 1000N, obtains three-dimensional motion information, add the certain noise influence therein, utilizing this method that elastic modulus is rebuild and imaging, obtain result as Fig. 2, can see substantially can the reflection bulk modulus true distribution.

Claims (2)

1, a kind of three-dimensional elastic modulus imaging method is characterized in that, may further comprise the steps:

(A) set up the three-dimensional finite element model of object original state;

(B) to object application system load;

(C) three-dimensional motion information of finite element node after application system load of measuring object, described three-dimensional motion information comprises displacement, speed;

(D) use obtains following continuous state estimate equation (20)～(23) based on the H ∞ Filtering Estimation method of state space, utilizes equation (20)～(23) to estimate the distributed in three dimensions of object elastic modulus, realizes the three-dimensional elastic modulus imaging of object;

$\left[\begin{array}{c}\widehat{\stackrel{\&CenterDot;}{x}}\\ \widehat{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}\end{array}\right]=\left[\begin{array}{cc}F& A_C\left(x\left(t\right)\right)\\ 0& 0\end{array}\right]\left[\begin{array}{c}\hat{x}\left(t\right)\\ \widehat{\mathrm{\&theta;}}\end{array}\right]+\left[\begin{array}{c}{B}_{C}w\left(t\right)\\ 0\end{array}\right]+\mathrm{\&Sigma;}{\left(t\right)}^{-1}{D}^T(y-D\hat{x}\left(t\right))---\left(20\right)$

$\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}=-\mathrm{\&Sigma;}\left(t\right)\left[\begin{array}{cc}F& A_C\left(x\left(t\right)\right)\\ 0& 0\end{array}\right]-\left[\begin{array}{cc}{F}^T& 0\\ A_C^T\left(x\left(t\right)\right)& 0\end{array}\right]\mathrm{\&Sigma;}\left(t\right)+\left[\begin{array}{cc}I& 0\\ 0& -{\mathrm{\&gamma;}}^{-2}Q\end{array}\right]-\mathrm{\&Sigma;}\left(t\right)\left[\begin{array}{cc}I& 0\\ 0& 0\end{array}\right]\mathrm{\&Sigma;}\left(t\right)---\left(21\right)$

$\left[\begin{array}{c}\hat{x}(t+\mathrm{\&Delta;t})\\ \widehat{\mathrm{\&theta;}}\end{array}\right]=\left[\begin{array}{c}\hat{x}\\ \widehat{\mathrm{\&theta;}}\end{array}\right]+\left[\begin{array}{c}\widehat{\stackrel{\&CenterDot;}{x}}\\ \widehat{\stackrel{\&CenterDot;}{\mathrm{\&theta;}}}\end{array}\right]\mathrm{\&Delta;t}---\left(22\right)$

$\mathrm{\&Sigma;}(t+\mathrm{\&Delta;t})=\mathrm{\&Sigma;}\left(t\right)+\stackrel{\&CenterDot;}{\mathrm{\&Sigma;}}\mathrm{\&Delta;t}---\left(23\right)$

Wherein:

$F=\left[\begin{array}{cc}0& I\\ 0& -\mathrm{\&alpha;I}\end{array}\right];$

$A_C=\left[\begin{array}{cc}0& 0\\ -M^{-1}{G}_1& -\mathrm{\&beta;}{M}^{-1}{G}_2\end{array}\right];$

$B_C=\left[\begin{array}{cc}0& 0\\ 0& M^{-1}\end{array}\right];$

$\hat{x}\left(t\right)={\left[\begin{array}{cc}\hat{U}\left(t\right)& \widehat{\stackrel{\&CenterDot;}{U}}\left(t\right)\end{array}\right]}^T;$

U is the motion vector of the finite element node of object, its single order differential

And second-order differential

Represent velocity vector and vector acceleration respectively;

I is and the corresponding unit matrix of U dimension;

α and β are the Rayleigh damping constant, desirable 0.01～0.05;

M is the finite element mass matrix of object, is the function of subject material density;

G ₁And G ₂It is the transition matrix in the computation process;

T is a time variable;

Δ t is interval estimated time;

Be the elastic modulus vector estimated value of object, the θ initial value preestablishes according to the experimental knowledge of object, and general requirement is approaching with the true elastic modulus of object;

W (t) is a control vector, $w\left(t\right)={\left[\begin{array}{cc}0& R\end{array}\right]}^T,$ R is the system load vector;

The three-dimensional motion information that y (t) obtains for measuring object in the step (C);

D is for measuring matrix, the line number of D is the finite element node number that three-dimensional motion information can measure, the columns of D is the dimension of x (t), in every row of D, if the three-dimensional motion information of the finite element node of index position correspondence can measure, the element of this index position is made as 1, otherwise element is made as 0;

Be the derivative of x (t), be illustrated in the variable quantity of concrete t x (t) constantly;

Modified value for θ;

γ influences the estimated value degree of accuracy for the binding occurrence of the quality of estimation, generally can get 1;

$\mathrm{\&Sigma;}=\left[\begin{array}{cc}{\mathrm{\&Sigma;}}_1& {\mathrm{\&Sigma;}}_2\\ {\mathrm{\&Sigma;}}_2^T& {\mathrm{\&Sigma;}}_3\end{array}\right],$ ∑ wherein ₁, ∑ ₂, ∑ ₃The initial value ∑ ₁(0), ∑ ₂(0), ∑ ₃(0) is respectively ∑ ₁(0)=Q ₁, ∑ ₂(0)=0, ∑ ₃(0)=Q ₀

Be the derivative of ∑, be illustrated in the variable quantity of concrete t ∑ constantly;

Q, Q ₀, Q ₁Be weight matrix.

2, three-dimensional elastic modulus imaging method as claimed in claim 1 is characterized in that, described transition matrix G ₁And G ₂Building method as follows:

(a) two N * N of initialization _eEmpty matrix G ₁And G ₂, N is the system node variable number, in three dimensions, each node has the degree of freedom of three directions, i.e. the total node number of system * 3, N _eBe the total unit number of system; The empty matrix K of a N * N of initialization _g', as the interim stiffness matrix of the system that does not contain elastic modulus E;

(b) for a unit that is numbered j, structure does not comprise this unitary elasticity modulus E _jStiffness matrix K _j';

The coding rule of unit and node during (c) according to system's generating mesh is K _j' in the interim stiffness matrix K of element insertion system _g' in the corresponding position, and make subscript in its unit into corresponding system subscript.

(d) K _g' with the vector that U multiplies each other and obtains, promptly be matrix G ₁J row; K _g' with

The vector that multiplies each other and obtain promptly is matrix G ₂J row;

(e) turn back to (b) step up to all node traversals, finish G ₁And G ₂Structure.

Images