一种三维弹性模量成像方法技术

本发明专利技术公开了一种三维弹性模量成像方法，包括以下步骤：（Ａ）建立对象初始状态的三维有限元模型；（Ｂ）向对象施加系统载荷；（Ｃ）测量对象的有限元节点在施加系统载荷后的三维运动信息，所述的三维运动信息包括位移、速度；（Ｄ）使用基于状态空间的Ｈ∞滤波估计方法，得到连续状态估计方程，利用连续状态估计方程估计出对象弹性模量的三维分布，实现对象的三维弹性模量成像。本发明专利技术的优点是基于状态空间的Ｈ∞滤波估计方法，结合有限元模型，实现在已知对象节点三维运动信息的情况下，鲁棒地得到对象的材料特性。

【技术实现步骤摘要】

本专利技术涉及一种三维弹性模量成像技术，具体地说是一种根据对象在力的作用下产生的形变反演弹性模量，从而实现成像的方法。
技术介绍
现代成像技术已经能够使人们观察到肉眼所无法看到的物体内部结构和运动变化的情况。早期的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。本文档来自技高网...

【技术保护点】
一种三维弹性模量成像方法，其特征在于，包括以下步骤：　（Ａ）建立对象初始状态的三维有限元模型；　（Ｂ）向对象施加系统载荷；　（Ｃ）测量对象的有限元节点在施加系统载荷后的三维运动信息，所述的三维运动信息包括位移、速度；　 　（Ｄ）使用基于状态空间的Ｈ∞滤波估计方法，得到如下连续状态估计方程（２０）～（２３），利用方程（２０）～（２３）估计出对象弹性模量的三维分布，实现对象的三维弹性模量成像；　＊＊＊　（２０）　＊＊＊　（２１）　＊＊＊　（２２）　 Σ（ｔ－Δｔ）＝Σ（ｔ）＋＊Δｔ　（２３）　其中：　＊＊＊　Ｕ为对象的有限元节点的位移向量，其一阶微分＊和二阶微分＊分别表示速度向量和加速度向量；　Ｉ为与Ｕ维数相应的单位矩阵；　α和β为瑞利阻尼常数，可取０．０１～０．０ ５；　Ｍ为对象的有限元质量矩阵，是对象材料密度的函数；　Ｇ↓［１］和Ｇ↓［２］是计算过程中的过渡矩阵；　ｔ为时间变量；　Δｔ为估计时间间隔；　＊为对象的弹性模量向量估计值，θ初值根据对象的经验知识预先设定，一般要求与对象的真实 弹性模量接近；　ｗ（ｔ）为控制向量，ｗ（ｔ）＝［０　Ｒ］↑［Ｔ］，Ｒ为系统载荷向量；　ｙ（ｔ）为步骤（Ｃ）中测量对象得到的三维运动信息；　Ｄ为测量矩阵，Ｄ的行数为三维运动信息可以测量得到的有限元节点数，Ｄ的列数为ｘ（ｔ）的维数，Ｄ的每 行中，索引位置对应的有限元节点的三维运动信息如果可以测量得到，该索引位置的元素设为１，否则元素设为０；　＊为ｘ（ｔ）的导数，表示在具体的ｔ时刻的ｘ（ｔ）的变化量；　＊为θ的修正值；　γ为估计的质量的约束值，影响估计值精确度，一般可 以取１；　＊＊＊，其中∑↓［１］、∑↓［２］、∑↓［３］的初值∑↓［１］（０）、∑↓［２］（０）、∑↓［３］（０）分别为∑↓［１］（０）＝Ｑ↓［１］、∑↓［２］（０）＝０、∑↓［３］（０）＝Ｑ↓［０］；　＊为∑的导数，表示在具体的ｔ 时刻的∑的变化量；　Ｑ、Ｑ↓［０］、Ｑ↓［１］为权重矩阵。...

【技术特征摘要】

【专利技术属性】
技术研发人员：金杰锋，刘华锋，
申请(专利权)人：浙江大学，
类型：发明
国别省市：86[中国|杭州]