专利名称：一种三维弹性模量成像方法
技术领域：
本发明涉及一种三维弹性模量成像技术，具体地说是一种根据对象在力的作用下产生的形变反演弹性模量，从而实现成像的方法。
背景技术：
现代成像技术已经能够使人们观察到肉眼所无法看到的物体内部结构和运动变化的情况。早期的X光还只能简单的表达不同物质的透射率，现如今，由于计算机图形图像技术的不断发展，超声波、MRI(磁共振成像)等新兴的成像技术手段不但能够更加清晰的显示物体的内部结构，而且还能够动态的获得物体内部质点的形变位移和速度信息。
很多时候，我们观察物体内部的结构和形变运动信息是为了了解其性质的变化，比如机械结构疲劳度探测、结构整体性分析等等。然而这并不是一个直接的反映。物质性质的变化首先引起了材料参数的变化，当受到外部的作用时，这种变化便会通过物体内部结构和运动的异常而表现出来。如果直接观察这些信息，并不能很直观地了解物质性质的变化状态。因此，获得物体内部的材料参数分布并直接的显示出来，对许多领域都有非常重要的意义。
在材料属性的各个参数中，最主要的是弹性模量，又称为杨氏模量，表示物体所受外部作用后产生的应力与应变之比，反映了材料的刚度，它是本方法中要估计的材料参数。
由于技术水平的限制，目前还缺少能够直接无创测量的弹性模量的方法。现在的研究方法一般是假设对象的形变运动信息已知(比如依靠植入种子获得或者通过计算机图形学方法测量形变运动信息，即位移值或者速度值为已知)，然后从这些已知的信息根据本构关系来估计材料参数如弹性模量。实际的成像探测手段在观测过程中会不可避免地受到噪声的影响。材料变化在运动状态上的反映本来就不明显，加上噪声的干扰，要从这些信息中获得材料参数就更加困难，只能通过一定的数学方法进行估计。
现实的物体都是三维的，目前全3D成像和运动探测技术、图像三维分析技术也在不断发展，利用三维运动信息直接估计整个物体的弹性模量更接近实际，而且能给人更直观的认识，从而做出更准确的判断。
发明内容
本发明的目的在于提供一种三维弹性模量成像方法，包括以下步骤 (A)建立对象初始状态的三维有限元模型； 三维有限元模型的建立可以根据现有三维有限元方法实现。
(B)向对象施加系统载荷； 所述的系统载荷可以是已知大小的作用力，也可以是测量得到的已经施加在对象上的作用力的变化值； (C)测量对象的有限元节点在施加系统载荷后的三维运动信息，所述的三维运动信息包括位移、速度； 由于现有技术中的不同测量方法的固有性质，使得有些对象的有限元节点三维运动信息可以测量得到，而有些有限元节点三维运动信息是无法测量的，需要进行估计； (D)使用基于状态空间的H∞滤波估计方法，得到如下连续状态估计方程(20)～(23)，利用方程(20)～(23)估计出对象弹性模量的三维分布，实现对象的三维弹性模量成像。
其中 U为对象的有限元节点的位移向量，其一阶微分
和二阶微分
分别表示速度向量和加速度向量； I为与U维数相应的单位矩阵； α和β为瑞利阻尼常数，可取0.01～0.05； M为对象的有限元质量矩阵，是对象材料密度的函数； G1和G2是计算过程中的过渡矩阵，构造方法如步骤(a)～(e) (a)初始化两个N×Ne的空矩阵G1和G2，N为系统节点变量数，在三维空间中，每个节点有三个方向的自由度，即系统总节点数×3，Ne为系统总单元数(参见图1中的模型，单元就是其中的一个四面体，节点就是四面体的四个端点，在系统中一个节点可能被多个单元共用，系统单元数就是对象分割成四面体单元的个数，建立模型网格的时候确定的)；初始化一个N×N(表示矩阵的行数和列数)的空矩阵Kg′，作为不含弹性模量E的系统临时刚度矩阵； (b)对于一个编号为j的单元，构造不包含该单元弹性模量Ej的刚度矩阵Kj′； (c)按照系统生成网格时单元和节点的编号规则，把Kj′中的元素插入系统临时刚度矩阵Kg′中相应的位置中，并将其单元内下标改成相应的系统下标； (d)Kg′与U相乘得到的向量，即是矩阵G1的第j列；Kg′与
相乘得到的向量，即是矩阵G2的第j列； (e)返回到第(b)步直到所有节点遍历，完成G1和G2的构造； t为时间变量； Δt为估计时间间隔；
为对象的弹性模量向量估计值，θ初值根据对象的经验知识预先设定，一般要求与对象的真实弹性模量接近； w(t)为控制向量，w(t)＝
T，R为系统载荷向量； y(t)为步骤(C)中测量对象得到的三维运动信息； D为测量矩阵，D的行数为三维运动信息可以测量得到的有限元节点数，D的列数为x(t)的维数，D的每行中，索引位置对应的有限元节点的三维运动信息如果可以测量得到，该索引位置的元素设为1，否则元素设为0；
为x(t)的导数，表示在具体的t时刻的x(t)的变化量；
为θ的修正值； γ为估计的质量的约束值，影响估计值精确度，一般可以取1； 其中∑1、∑2、∑3的初值∑1(0)、∑2(0)、∑3(0)分别为∑1(0)＝Q1、∑2(0)＝0、∑3(0)＝Q0；
为∑的导数，表示在具体的t时刻的∑的变化量； Q、Q0、Q1为权重矩阵，其取值根据需要设定，表征不同的因素对对象弹性模量的三维分布的影响程度； 各参数中的上标^为该参数的估计值，上标T表示该矩阵的转置矩阵。
本发明的优点是基于状态空间的H∞滤波估计方法，结合有限元模型，实现在已知对象节点三维运动信息的情况下，鲁棒地得到对象的材料特性。
图1是本发明用于验证算法的仿真模型； 图2是图1中的仿真模型在图2所示的位移场基础上使用本方法重建出的弹性模量分布示意图。
具体实施例方式 下面结合附图对本发明的具体实施方式
进行详细说明，本发明弹性模量成像方法，包括以下步骤 (A)建立对象初始状态的三维有限元模型； (B)向对象施加已知大小的作用力，或者测量已经施加在对象上的作用力的变化值，作为系统载荷； (C)测量对象节点在施加系统载荷后，或者是已有的系统载荷变化后的三维运动信息，所述的三维运动信息包括位移、速度； (D)根据关于对象的经验知识，设定与所成像对象的真实弹性模量接近的一个先验的值，使用有限元方法建立对象的力学模型，建立系统动态方程；根据系统动态方程建立状态方程组，使用基于状态空间的H∞滤波估计方法，估计出对象弹性模量的三维分布，实现成像。
在三维问题中，常见的单元类型是四面体单元和六面体单元。由于四面体单元对于复杂形状的网格化有更好的适应性，网格划分算法也比较成熟，而且计算量小，因此这里采用这种单元形状。
四面体单元内的位移u可以仅用其四个节点的位移ue来表示(下标e均表示一个单元中的变量) u＝Nue (1) 其中N为四面体单元形状函数矩阵，具体形式 其中Ni，Nj，Nk，Nm分别是四面体单元四个节点对应的形状函数，可以由下式计算得到 记由四面体单元四个顶点坐标构成的矩阵 则ap，n为矩阵Λ的第p行，第n列元素的代数余子式，单元体积Ve＝|Λ|/6。
从计算的角度考虑，假设对象是各项同性的线性弹性体，其应力向量σ与应变向量ε成线性关系 σ＝Dε(5) 其中，材料矩阵D E是弹性模量，表示物体产生形变的难易度；v是泊松比，表示横向应变与纵向应变之比。这两项都属于物体的材料性质，实际中泊松比的变化较小，可以认为是常量。弹性模量是本文中要估计的材料参数。
空间中的应变与位移的关系由三维几何方程表示 由式(1)和(7)可以得到 ε＝LNue＝Bue(8) 由此可见四面体单元内各点的应变是一样的，是常应变单元，应变大小由四个节点的位移确定，B是单元应变矩阵。四面体单元刚度矩阵的计算如下式 Ke＝BTDBVe(9) 四面体单元的质量矩阵采用集中质量矩阵的形式，即认为质量集中在单元的节点上 其中ρ是对象的材料平均密度，I是12×12的单位矩阵。
根据物理学原理，导出有限元模型节点位移系统动态方程的矩阵形式 整个对象的所有节点组成一个系统，式(4)中M表示系统的质量矩阵，是对象材料密度的函数；C是系统的阻尼矩阵；K是系统的刚度矩阵，是对象的应变张量和应力张量之间的关系函数；R是系统所有节点的载荷向量(即系统载荷)，反映对象所受的作用力；U是系统所有节点的位移向量，其一阶微分
和二阶微分
分别表示速度向量和加速度向量；C与频率有关，假设为瑞利阻尼，则C满足公式C＝αM+βK(α、β取0.01-0.05)；将C表达式代入(11)得到式(12) 重建弹性模量，即通过重组系统动态方程(12)，利用已有的位移信息和边界力信息来求得弹性模量E的分布。因此需要把方程(12)的左边变换成E向量的函数，即KU＝G1E，由方程(6)可知这样的变换是可行的，由此得到 根据有限元方法中单元刚度矩阵Ke的定义，Ke可以表示成含单元弹性模量Ee项的形式(下标e表示单元，Ke和Ee分别是一个单元的刚度矩阵和弹性模量，没有下标的就表示整个系统的刚度矩阵和弹性模量) 由此，根据刚度矩阵构造的方法，得到构造G1和G2的方法 (a)初始化两个N×Ne的空矩阵G1和G2，N为系统节点变量数，在三维空间中，每个节点有三个方向的自由度，即系统总节点数×3，Ne为系统总单元数(见附图中的模型，单元就是其中的一个四面体，节点就是四面体的四个端点，在系统中一个节点可能被多个单元共用，系统单元数就是对象分割成四面体单元的个数，建立模型网格的时候确定的)；初始化一个N×N(表示矩阵的行数和列数)的空矩阵Kg′，作为不含弹性模量E的系统临时刚度矩阵； (b)对于一个编号为j的单元，利用方程(14)构造不包含该单元弹性模量Ej的刚度矩阵Kj′； (c)按照系统生成网格时单元和节点的编号规则，把Kj′中的元素插入系统临时刚度矩阵Kg′中相应的位置中，并将其单元内下标改成相应的系统下标。
(d)Kg′与U相乘得到的向量，即是矩阵G1的第j列；Kg′与
相乘得到的向量，即是矩阵G2的第j列； (e)返回到第(b)步直到所有节点遍历。
由此构造出N×Ne的矩阵G1，G2。
为了解决测量运动量过程中存在噪声的问题，采用基于H∞的噪声扰动全状态信息方法(NPFSI)来实现估计过程。
为了应用该滤波估计方法，首先建立系统状态方程。设状态向量为控制向量为w(t)＝
T，将公式(12)变形，建立连续时间线性随机系统状态空间表达式 y(t)＝Dx(t)+e(t)(16) D为测量矩阵，D的行数为三维运动信息可以测量得到的有限元节点数，D的列数为x(t)的维数，D的每行中，索引位置对应的有限元节点的三维运动信息如果可以测量得到，该索引位置的元素设为1，否则元素设为0。v(t)代表过程噪声，e(t)代表观测噪声。
由于弹性模量E的是包含在参数矩阵中的，为了实现我们估计的目的，我们设材料向量为θ＝[E]，再将原运动状态量x和材料参数θ组成一个新的状态向量，将方程(15)(16)变形，得到如下一组新的状态方程，这里我们假设材料参数θ不随时间变化 其中 H∞最优估计问题，是通过设计状态估计器，使对任意情况的过程噪声和观测噪声，估计误差均较小，即实现未知噪声下的最优估计。上述过程的价值函数如下 H∞滤波器实际上为外部干扰与估计量之间的博弈问题，即一个最小最大化问题。由状态方程得到如下一组连续状态估计方程 为了使估计算法达到预期的效果，在开始估计之前，必须先设定合适的噪声干扰衰减水平，权重矩阵Q，Q0，Q1和γ值。通常很难决定最优的γ*，它可能是无穷的，意味着无法达到期望的噪声控制水平。然而，实际上如果Q被选定了，γ*能够通过分析获得。首先将∑矩阵分块 将(22)代入(21)式，得到分块矩阵的时间微分更新方程 ∑1(0)＝Q1(25) ∑2(0)＝0 (26) ∑3(0)＝Q0(27) 根据Schur原理，要保证∑＞0，只有当 ∑1(t)＞0(26) ∏(t)是下面矩阵微分方程的解 ∏(0)＝Q0(29) 由于本文假设要估计的参数是时间常数，所以令γ＝1，∑1(0)＝I。
在计算过程中，弹性模量E的初值可以根据对心肌的经验知识，设定与所成像对象的真实弹性模量接近的某一个先验的值，然后利用方程(20)，(21)(或(25-27))，(22)，(23)对预设的弹性模量E进行递归估计修正，直到获得收敛的最优解，即得到弹性模量E的分布信息，再使用计算机3D技术显示出来，实现材料弹性模量成像。
仿真实验模型为三维半椭球壳，高100mm，内径35mm，外径50mm，其杨氏模量分布如图1所示。浅色部分表示较软部分，杨氏模量为75kPa。深色部分代表较硬部分，杨氏模量为105kPa。将模型的上表面固定，内壁施加1000N内的压力，获得三维运动信息，在其中加入一定噪声影响，在利用本方法对弹性模量进行重建和成像，得到如图2的结果，可以看到基本能够反映物体弹性模量的真实分布。
权利要求
1、一种三维弹性模量成像方法，其特征在于，包括以下步骤
(A)建立对象初始状态的三维有限元模型；
(B)向对象施加系统载荷；
(C)测量对象的有限元节点在施加系统载荷后的三维运动信息，所述的三维运动信息包括位移、速度；
(D)使用基于状态空间的H∞滤波估计方法，得到如下连续状态估计方程(20)～(23)，利用方程(20)～(23)估计出对象弹性模量的三维分布，实现对象的三维弹性模量成像；
其中
U为对象的有限元节点的位移向量，其一阶微分
和二阶微分
分别表示速度向量和加速度向量；
I为与U维数相应的单位矩阵；
α和β为瑞利阻尼常数，可取0.01～0.05；
M为对象的有限元质量矩阵，是对象材料密度的函数；
G1和G2是计算过程中的过渡矩阵；
t为时间变量；
Δt为估计时间间隔；
为对象的弹性模量向量估计值，θ初值根据对象的经验知识预先设定，一般要求与对象的真实弹性模量接近；
w(t)为控制向量，R为系统载荷向量；
y(t)为步骤(C)中测量对象得到的三维运动信息；
D为测量矩阵，D的行数为三维运动信息可以测量得到的有限元节点数，D的列数为x(t)的维数，D的每行中，索引位置对应的有限元节点的三维运动信息如果可以测量得到，该索引位置的元素设为1，否则元素设为0；
为x(t)的导数，表示在具体的t时刻的x(t)的变化量；
为θ的修正值；
γ为估计的质量的约束值，影响估计值精确度，一般可以取1；
其中∑1、∑2、∑3的初值∑1(0)、∑2(0)、∑3(0)分别为∑1(0)＝Q1、∑2(0)＝0、∑3(0)＝Q0；
为∑的导数，表示在具体的t时刻的∑的变化量；
Q、Q0、Q1为权重矩阵。
2、如权利要求1所述的三维弹性模量成像方法，其特征在于，所述的过渡矩阵G1和G2的构造方法如下
(a)初始化两个N×Ne的空矩阵G1和G2，N为系统节点变量数，在三维空间中，每个节点有三个方向的自由度，即系统总节点数×3，Ne为系统总单元数；初始化一个N×N的空矩阵Kg′，作为不含弹性模量E的系统临时刚度矩阵；
(b)对于一个编号为j的单元，构造不包含该单元弹性模量Ej的刚度矩阵Kj′；
(c)按照系统生成网格时单元和节点的编号规则，把Kj′中的元素插入系统临时刚度矩阵Kg′中相应的位置中，并将其单元内下标改成相应的系统下标。
(d)Kg′与U相乘得到的向量，即是矩阵G1的第j列；Kg′与
相乘得到的向量，即是矩阵G2的第j列；
(e)返回到第(b)步直到所有节点遍历，完成G1和G2的构造。
全文摘要
本发明公开了一种三维弹性模量成像方法，包括以下步骤(A)建立对象初始状态的三维有限元模型；(B)向对象施加系统载荷；(C)测量对象的有限元节点在施加系统载荷后的三维运动信息，所述的三维运动信息包括位移、速度；(D)使用基于状态空间的H∞滤波估计方法，得到连续状态估计方程，利用连续状态估计方程估计出对象弹性模量的三维分布，实现对象的三维弹性模量成像。本发明的优点是基于状态空间的H∞滤波估计方法，结合有限元模型，实现在已知对象节点三维运动信息的情况下，鲁棒地得到对象的材料特性。
文档编号G01N3/00GK101634618SQ200910101968
公开日2010年1月27日 申请日期2009年8月27日 优先权日2009年8月27日
发明者金杰锋, 刘华锋 申请人:浙江大学