CN103049663B - 磁共振弹性成像中的弹性模量重建方法和系统 - Google Patents

Abstract

本发明提供了一种磁共振弹性成像中的弹性模量重建方法和系统。所述方法包括：将成像组织的表面假设为平面域，通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数；通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代；将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。采用本发明能降低计算量。

Description

磁共振弹性成像中的弹性模量重建方法和系统

技术领域

本发明涉及磁共振成像技术，特别是涉及一种磁共振弹性成像中的弹性模量重建方法和系统。

背景技术

弹性是人体组织物理性质中一种重要的机械力学参数，生物组织的弹性变化通常是与一定的病理现象紧密相关的，也就是说，病变组织和正常组织往往存在着弹性模量的差异，这一差异为临床上疾病的诊断提供了重要的参考信息。

磁共振弹性成像（MagneticResonanceElastography，简称MRE）作为一种无创成像方法，能够直观地显示和量化人体内部组织弹性，实现对人体内部组织的弹性成像，使得“影像触诊”成为了可能，在乳腺癌检测、肝硬化分期、动脉粥样硬化斑块、肌肉损伤、大脑疾病检测和射频消融等治疗和监控方面具有重要意义。

磁共振弹性成像中弹性模量重建的方法是一个由质点位移图反推弹性分布的逆问题求解，因此，其本质上是不稳定的。为了避免该问题的病态性，弹性弹性模块重建方法将根据应用范围进行假设和简化。目前提出的弹性模量重建方法包括：（1）局部频率估计（LocalFrequencyEstimation，简称LFE）算法及其变种，该算法将假设介质是均匀的和不可压缩的，并忽略波动中的衰减，机械波在介质中的传播方程因而简化为亥姆霍兹方程，以该方程为模型进行直接逆问题代数求解，但是局部频率估计算法存在着分辨率低、精度有限的缺陷，对尖锐的边界无法估计出精确的弹性系数，其假设也不适用于某些临床中；（2）基于有限元分析的弹性模量重建算法，计算出一幅质点位移图，通过最小化该质点位移图和磁共振质点位移图得到弹性系数分布图，与局部频率估计算法及其变种相比较，该方法对介质等没有做特定假设，对噪声的干扰不敏感，可产生较高分辨率的图像，但计算量非常庞大。

发明内容

基于此，提供一种能降低计算量的磁共振弹性成像中的弹性模量重建方法。

此外，还有必要提供一种能降低计算量的磁共振弹性成像中的弹性模量重建系统。

一种磁共振弹性成像中的弹性模量重建方法，包括如下步骤：

将成像组织的表面假设为平面域，通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数；

通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代；

将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。

在其中一个实施例中，所述通过有限体积元算法求解得到所述平面域中的位移初值的步骤为：

根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。

在其中一个实施例中，所述通过有限体积元算法求解得到所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤包括：

将所述平面域划分为若干个子区域，并划分所述子区域为若干个单元，所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上，平面域边界的顶点为单元的顶点；

在所述子区域中，以构成单元的顶点作为所述子区域的节点，并构建所述节点对应的有限体积元方程，并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移；

通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数，所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。

在其中一个实施例中，所述通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的步骤包括：

通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值；

根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。

在其中一个实施例中，所述根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值的步骤之后还包括：

获取通过对所述成像组织进行磁共振成像所得到的位移图；

根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差；

判断所述最优化问题平方差是否小于预设的容忍误差，若是，则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中，若是，则

判断所述单元中的最小单元对应的迭代次数是否达至预设的迭代次数，若是，则

停止进行牛顿迭代。

在其中一个实施例中，还包括：

若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数，则

返回所述通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤。

一种磁共振弹性成像中的弹性模量重建系统，包括：

有限体积元运算模块，用于将成像组织的表面假设为平面域，通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数；

迭代模块，用于通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代；

分布形成模块，用于将所述最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布。

在其中一个实施例中，所述有限体积元运算模块还用于根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。

在其中一个实施例中，所述有限体积元运算模块包括：

划分单元，用于将所述平面域划分为若干个子区域，并划分所述子区域为若干个单元，所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上，平面域边界的顶点为单元的顶点；

子区域位移求解单元，用于在所述子区域中，以构成单元的顶点作为所述子区域的节点，并构建所述节点对应的有限体积元方程，并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移；

子区域导数求解单元，用于通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数，所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。

在其中一个实施例中，所述迭代模块包括：

改进值运算单元，用于通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值；

弹性模量迭代单元，用于根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。

在其中一个实施例中，所述迭代模块还包括：

位移图处理单元，用于获取通过对所述成像组织进行磁共振成像所得到的位移图，并根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差；

判断单元，用于判断所述最优化问题平方差是否小于预设的容忍误差，若是，则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中，若是，则判断所述单元中的最小单元对应的迭代次数是否达到预设的迭代次数，若是，则停止进行牛顿迭代。

在其中一个实施例中，所述判断单元还用于若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数，则通知所述有限体积元运算模块。

上述磁共振弹性成像中的弹性模量重建方法和系统，将成像组织的表面假设为平面域，引入有限体积元算法以求解得到位移和该位移对未知弹性模量的导数，进而通过求解得到的位移和导数进行牛顿迭代得到弹性模量值，并在弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代，进而将最终迭代得到的弹性模量值组成成像组的弹性模量分布，通过引入有限体积元算法，使得弹性模量的重建在保证了计算精度的前提下降低了计算量，提高了计算速度。

附图说明

图1为一个实施例中磁共振弹性成像中的弹性模量重建方法的流程图；

图2为图1中通过有限体积元算法求解得到平面域中子区域的位移、位移对未知弹性模量的导数的方法流程图；

图3为一个实施例中子区域的原始剖分的示意图；

图4为一个实施例中对子区域的某一节点的外心对偶剖分单元的示意图；

图5为一个实施例中通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的方法流程图；

图6为另一个实施例中通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的方法流程图；

图7为一个实施例中平面域的示意图；

图8为一个实施例中磁共振弹性成像中的弹性模量重建系统的结构示意图；

图9为图8中有限体积元运行模块的结构示意图；

图10为一个实施例中迭代模块的结构示意图；

图11为另一个实施例中迭代模块的结构示意图。

具体实施方式

如图1所示，一种磁共振弹性成像中的弹性模量重建方法，包括如下步骤：

步骤S10，将成像组织的表面假设为平面域，通过有限体积元算法求解得到平面域中的位移初值，以及平面域中子区域的位移、位移对未知弹性模量的导数。

本实施例中，在一定条件下，可将一个在外力作用下处于平衡状态的弹性体视为平面弹性问题，而在磁共振弹性成像中，成像组织被施加的外力是非常小的，可将该成像组织视为处于平衡状态的，因此，可将磁共振弹性成像中应力和应变之间的关系视为平面弹性问题。

将成像组织的表面视为平面域Ω，获取通过平面域对应的平衡方程所构建得到的二阶椭圆微分方程组，并在构建的二阶椭圆微分方程组中对未知弹性模量求导得到包含了位移对未知弹性模量的导数的方程。

进一步的，对于平面域Ω，其边界为描述成像组织平衡的状态变量包括三组，即，应力张量σ＝(σ₁₁,σ₂₂,σ₁₂)^T、应变张量ε＝(ε₁₁,ε₂₂,ε₁₂)^T和位移向量u＝(u₁,u₂)^T。假设Ω是均匀各向同性的弹性体，令

$B(\▿)=\left(\begin{array}{ccc}\frac{\∂}{\∂x_1}& 0& \frac{\∂}{\∂x_2}\\ 0& \frac{\∂}{\∂x_2}& \frac{\∂}{\∂x_1}\end{array}\right)$

$A=\left(\begin{array}{ccc}\mathrm{\λ}+2\mathrm{\μ}& \mathrm{\λ}& 0\\ \mathrm{\λ}& \mathrm{\λ}+2\mathrm{\μ}& 0\\ 0& 0& \mathrm{\μ}\end{array}\right)$

其中，正数λ和μ是Lamé常数：

$\mathrm{\λ}=\frac{\mathrm{vE}}{(1+v)(1-2v)},$ $\mathrm{\μ}=\frac{E}{2+2v}$

其中，ν是Poisson系数（Poisson’sRatio，是已知的常数），E是弹性模量（Young’sModulus）。σ，ε，u满足如下三组方程：

（应变-位移关系）（1）

（平衡方程）（2）

σ＝Aε（应力-应变关系）（3）

其中，f是体积力。

由Green公式可推出：

${\∫}_{\mathrm{\Ω}}{\mathrm{\σ}}^T{B}^T(\▿)\mathrm{udx}+{\∫}_{\mathrm{\Ω}}{\left(B(\▿)\mathrm{\σ}\right)}^T\mathrm{udx}={\∫}_{\∂\mathrm{\Ω}}{\left(B\left(\mathrm{\υ}\right)\mathrm{\σ}\right)}^T\mathrm{udx}$

其中，υ＝(υ₁,υ₂)^T是的单位外法向量。

假定Γ分为两段不重叠的线段Γ₀和Γ₁，在Γ₀上给定位移边界条件：在Γ₁上给定力条件（表面力）。

在求解时，利用公式（1）-（3）消去部分未知量，留下的是待求量，即弹性模量E，得到包含了位移向量u的二阶椭圆微分方程组：

$-\mathrm{\μ}\▿u-(\mathrm{\λ}+\mathrm{\μ})\mathrm{grad}\mathrm{div}u=f---\left(4\right).$

其中，以分别乘（4）式两端，关于x∈Ω积分，然后利用Green公式，得积分形式：

$a(u,v)-{\∫}_{\mathrm{\Ω}}[\mathrm{\μ}\frac{\∂u}{\∂\mathrm{\υ}}+(\mathrm{\λ}+\mathrm{\μ})\left(\mathrm{div}u\right)v]\mathrm{vds}=(f,v)---\left(5\right)$

其中

$a(u,v)={\∫}_{\mathrm{\Ω}}[\mathrm{\μ}\▿u\▿v+(\mathrm{\λ}+\mathrm{\μ})\mathrm{div}u\·\mathrm{div}v]\mathrm{dx}---\left(6\right)$

是边界上的法向导数，是Ω上的Hilbert空间，定义在其上的函数在平面域Ω的边界上的值为零.

在边界Γ₁上有

$\mathrm{\μ}\frac{\∂u}{\∂\mathrm{\υ}}+(\mathrm{\λ}+\mathrm{\μ})\mathrm{div}u={\left(B\left(\mathrm{\υ}\right)\mathrm{\σ}\right)}^T=\stackrel{\‾}{P}$

故方程（4）的变分形式是：求u∈(H¹(Ω))²，使得

$a(u,v)=(f,v)+{\∫}_{{\mathrm{\Γ}}_1}\stackrel{\‾}{P}\mathrm{vds},\∀v\∈{\left(H_0^1\left(\mathrm{\Ω}\right)\right)}^2---\left(7\right)$

其中，H¹(Ω)是平面域Ω上的Hilbert空间.

进一步的，为求解出位移对未知弹性模量的导数，将对（4）式两端对未知弹性模量求导，以得到包含了位移对未知弹性模量的导数的方程。

在构建得到了二阶椭圆微分方程组和包含了位移对未知弹性模量的导数的方程之后，将应用有限体积元算法求解位移和位移对未知弹性模量的导数。

有限体积元法是吸取了有限差分法和有限元法的优点发展起来的，又称广义差分法，可以处理复杂的边值条件及不规则区域，并且由于离散化得到的线性方程组是稀疏的，在不降低数值解的收敛阶的前提下求解方程的计算量小，计算速度快，因此，将有限体积元法应用于磁共振弹性成像中将大大地降低了弹性模量重建的计算量。

在一个实施例中，上述通过有限体积元算法求解得到平面域中的位移初值的具体过程为：根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。

本实施例中，将首先进行位移初值的求解，以方便后续的位移求解过程。在位移初值的求解过程中，将进行弹性模量初值的假设，并根据这一假设通过有限体积元算法得到平面域所对应的位移初值，其中，该位移初值计算时所使用的边界Γ上的初值取为对成像组织进行磁共振成像所得到的位移图中对应的位移值。

如图2所示，在一个实施例中，上述通过有限体积元算法求解得到平面域中子区域的位移、位移对未知弹性模量的导数的步骤包括：

步骤S110，将平面域划分为若干个子区域，并划分子区域为若干个单元，子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上，平面域边界的顶点为单元的顶点。

本实施例中，将平面域Ω划分为若干个子区域，然后分别划分每一子区域为若干个单元之和，在优选的实施例中，划分得到的单元为三角形的形状，进而使得子区域之间、单元之间互相不存在重叠区域，并且任一单元的顶点不会在其它单元的边上，边界Γ的每一顶点均为单元的顶点，从而每一子区域将对应了一个原始剖分T_h，如图3所示，其中h是所有三角形单元边的最长边。记K_Q为T_h内的三角形单元，Q是三角形单元的外心。

进一步的，如图4所示，记是原始剖分T_h的对偶剖分，包括是重心对偶剖分和外心对偶剖分。以下设是外心对偶剖分。图4表示的是原始剖分T_h中以P₀为顶点的所有三角单元和围绕P₀的对偶单元（图中阴影剖分）。外心对偶剖分做法：设T_h的任一单元的内角不大于90°，取三角单元□P₀P_iP_i+1(i＝1,2,…6,P₇＝P₁)的外心Q_i为对偶剖分的节点，即依次连结各个外心点可得到外心对偶剖分，此时，是的中垂线，分别过边中点M_i。

步骤S130，在子区域中，以构成单元的顶点作为子区域的节点，并构建节点对应的有限体积元方程，并以位移初值作为有限体积元方程中的初值进行求解得到子区域对应的位移。

本实施例中，对于每一子区域，设是定义在原始分剖T_h上的试探函数空间，它是分片线性函数空间，即为一次多项式，完全由三角形单元中三个顶点上的值所确定，将T_h的内节点编号为1,2,…,N₀。其中，节点分两类，给定力条件的节点编号为N₀+1,…,N₁，给定位移条件的节点编号为N₁+1,…，N。用表示节点i∈{1,2,…,N₁}的基函数，则对u_h∈U_h，可表示为：（Γ_0h是Γ₀的近似）

u_i是u_h(x)在第i个节点x_i的值。又假设V_h是定义在的检验函数空间，它为分片常数空间，即对任意的内节点P₀，相应的基函数为的特征函数：

${\mathrm{\ψ}}_{P_0}\left(P\right)=\left\{\begin{array}{c}1,P\∈K_{P_0}^{*}\\ 0,P\∉K_{P_0}^{*}\end{array}\right.$

设节点基函数为ψ_j＝(ψ_1j，ψ_2j)^T。则得到有限体积元法方程是：

求u_h∈U_h，使得

$a(u_h,{\mathrm{\ψ}}_j)=(f,{\mathrm{\ψ}}_j)+{\∫}_{{\mathrm{\Γ}}_{1h}}\stackrel{\‾}{P}{\mathrm{\ψ}}_j\mathrm{ds},j=\mathrm{1,2},...,N_1---\left(8\right)$

其中，Γ_1h是Γ₁的近似。

将式（4）两端在上积分，应用Green公式并用u_h代替u，得

$-{\∫}_{\∂K_{P_0}^{*}}[\mathrm{\μ}\frac{\∂u_h}{\∂\mathrm{\υ}}+(\mathrm{\λ}+\mathrm{\μ})\left(\mathrm{div}{u}_h\right)\mathrm{\υ}]\mathrm{ds}={\∫}_{K_{P_0}^{*}}\mathrm{fdx}---\left(9\right)$

此即为在节点P₀的有限体积元方程。下面给出（9）式左端积分的计算。第一项积分分解为中垂线段上的积分和。例如沿的积分为

$-\mathrm{\μ}{\∫}_{\stackrel{\‾}{Q_1{Q}_2}}\frac{\∂u_h}{\∂\mathrm{\υ}}\mathrm{ds}=-\mathrm{\μ}\stackrel{\‾}{Q_1{Q}_2}\frac{u_{P_2}-u_{P_0}}{\stackrel{\‾}{P_0{P}_2}}---\left(10\right)$

其余中垂线段上的积分类推。第二项积分分解为沿折线 上的积分。例如，沿上的积分

$-(\mathrm{\λ}+\mathrm{\μ}){\∫}_{\stackrel{\‾}{M_1{Q}_1{M}_2}}\left(\mathrm{div}{u}_h\right)\mathrm{\υds}$

$=-(\mathrm{\λ}+\mathrm{\μ})\mathrm{div}{u}_h\left(Q_1\right)(\frac{\stackrel{\→}{P_0{P}_1}}{\stackrel{\‾}{P_0{P}_1}}\stackrel{\‾}{M_1{Q}_1}+\frac{\stackrel{\→}{P_0{P}_2}}{\stackrel{\‾}{P_0{P}_2}}\stackrel{\‾}{Q_1{M}_2})---\left(11\right)$

其中

$\mathrm{div}{u}_h\left(Q_1\right)=\frac{\∂u_{1h}}{\∂x_1}\left(Q_1\right)+\frac{\∂u_{2h}}{\∂x_2}\left(Q_1\right)$

$=\frac{1}{2S_{Q_1}}[(x_2\left(P_1\right)-x_2\left(P_2\right))u_1\left(P_0\right)+(x_2\left(P_2\right)-x_2\left(P_0\right))u_1\left(P_2\right)+(x_2\left(P_0\right)-x_2\left(P_1\right))u_1\left(P_2\right)]$

$+\frac{1}{2S_{Q_1}}[(x_1\left(P_2\right)-x_1\left(P_1\right))u_2\left(P_0\right)+(x_1\left(P_0\right)-x_1\left(P_2\right))u_2\left(P_1\right)+(x_1\left(P_1\right)-x_1\left(P_0\right))u_2\left(P_2\right)]---\left(12\right)$

其中，(x₁(P_i),x₂(P_i))是点P_i的坐标，是含外心Q₁的单元面积。其余折线段上的积分类推。

方程（10）是内点有限体积元法方程。在边界Γ_0h上给定位移值u₀。例如假设节点此时在处仍有方程（9），左端的线积分仍按公式（10）和（11）计算。

总之可列出2N₁个形如（9）式的有限体积元方程和Γ_0h上的位移边界条件。也就是说，假设已知给定的表面力，和正数λ、μLamé常数，就可以计算出位移所对应的数值。

步骤S150，通过求解得到的位移对包含了位移对未知弹性模量的导数的方程进行求解得到位移对未知弹性模量的导数，位移对未知弹性模量的导数是与位移所在的子区域相对应的。

本实施例中，在式（4）中对每一子区域的未知弹性模量求导数，即可得到慎勿使形谍了位移对未知弹性模量的导数的方程，即：

$-\mathrm{\μ}\▿u_j^{\′}-(\mathrm{\λ}+u)\mathrm{grad}\mathrm{div}{u}_j^{\′}=\frac{\∂\mathrm{\μ}}{\∂E_{z_j}}\▿u+\frac{\∂(\mathrm{\λ}+\mathrm{\μ})}{\∂E_{z_j}}\mathrm{grad}\mathrm{div}u---\left(13\right)$

其中z用于子区域z，（13）式的形式如同（4），这里是关于u′_j的方程，方程右端可看为是（4）式中的f，因此可以用解（4）式的有限体积元法解（13）式，得到u′_j的解，在此不再赘述。

步骤S30，通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代。

本实施例中，应用牛顿迭代法进行弹性模量值的，假设迭代第n步的解为则第n+1步的解为

$E_z^{(n+1)}=E_z^{\left(n\right)}+\mathrm{\Δ}{E}_z---\left(14\right)$

其中ΔE_z是弹性模量E_z的改进值，它是如下正则化方程组的解

$\left[(H_z^{\left(n\right)}+\mathrm{\αI})\right]\mathrm{\Δ}{E}_z=-f_z^{\left(n\right)}---\left(15\right)$

是不计二阶导数项的海森（Hessian）矩阵，元素为

$h_{\mathrm{ij}}^{\left(n\right)}=\frac{\∂f_i^{\left(n\right)}}{\∂E_{z_j}^{\left(n\right)}}=-\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_i}^{\left(n\right)}}\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_j}^{\left(n\right)}}+\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_i}^{\left(n\right)}}\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_j}^{\left(n\right)}})---\left(16\right)$

α是正则化对角矩阵的参数，使可逆，I是单位矩阵。

因此，由上述推导可知，在进行牛顿迭代计算弹性模量值时，将通过有限体积元算法求解得到的位移和导数经由式（15）和式（16）计算得到改进值，进而方可进行下一步的迭代。

如图5所示，在一个实施例中，上述通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的步骤包括：

步骤S310，通过每一子区域对应的位移和导数进行运算得到子区域对应的弹性模量改进值。

步骤S330，根据弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。

本实施例中，通过计算得到的弹性模量改进值和上一次迭代所得到的弹性模量值通过式（14）进行计算得到当前迭代次数所对应的弹性模量值。

如图6所示，在一个实施例中，上述步骤S330之后还包括如下步骤：

步骤S301，获取通过对成像组织进行磁共振成像所得到的位移图。

步骤S302，根据当前迭代得到的弹性模量值和位移图得到最优化问题平方差。

本实施例中，由于磁共振成像中的弹性模量重建问题被归结为一个带约束条件的最优化问题，其目标函数即为求解得到磁共振成像所测得的位移与计算得到的位移之间的最小平方和差。

该最优化问题为：minF(E)

其中

$F\left(E\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[(u_i^{\mathrm{meas}}-u_i^{\mathrm{comp}})+(v_i^{\mathrm{meas}}-v_i^{\mathrm{comp}})]$

这里和分别是在成像组织位置i处通过磁共振测得的x方向和y方向的位移，和在成像组织位置i处用有限体积元法计算出的位移，一共有N个不同的位置。E是一个M维的弹性参数向量，它由一个连续基底φ张成，用它来定义整个感兴趣区域的弹性组织分布。

如图7所示，假设整个感兴趣区域为Ω，将Ω分为若干个子区域Ω_sub，可将最优化问题重写为子区域上的最优化问题。假设一共有Q个子区域，有

$F\left(E\right)=\underset{z=1}{\overset{Q}{\mathrm{\Σ}}}{F}_z\left(E_z\right)$

其中，F_z(E_z)是第z个子区域上的最优化函数，对所有的子区域求和可得最优化问题为：

$\min F\left(E\right)=\min \left[\underset{z=1}{\overset{Q}{\mathrm{\Σ}}}{F}_z\left(E_z\right)\right]\⇒\underset{z=1}{\overset{Q}{\mathrm{\Σ}}}\min F_z\left(E_z\right)$

这里

$F_z\left(E_z\right)=\underset{i_Z}{\overset{N_z}{\mathrm{\Σ}}}[{(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})}^2+{(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})}^2]$

这里和分别是在成像组织第z个子区域上i_z处磁共振测得的x方向和y方向的位移，和在成像组织z个子区域上i_z处用有限体积元法计算出的位移，一共有N_z个不同的位置。E_z是M_z的向量，有N_z＜N，M_z＜M。对于在方程两边分别对Ez求一阶导数，并令导数等于0，可得如下非线性方程组：

$\begin{array}{c}{f}_1=\frac{\∂F_z}{\∂E_{z_1}}=\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_1}}+\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_1}}=0\\ f_2=\frac{\∂F_z}{\∂E_{z_2}}=\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_2}}+\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_2}}=0\\ \begin{array}{ccccc}.& & .& & .\\ .& & .& & .\\ .& & .& & .\end{array}\\ f_{M_z}=\frac{\∂F_z}{\∂E_{z_{M_z}}}=\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_{M_z}}}+\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_{M_z}}}=0\\ \end{array}$

进而通过如上所述的牛顿迭代法求解上述非线性方程组。

步骤S303，判断最优化问题平方差是否小于预设的容忍误差，若是，则进入步骤S304，若否，则返回步骤S310。

本实施例中，预设的容忍误差是根据实际需要进行灵活设定的。

步骤S304，进一步判断平面域中的每一单元是否均包含于至少一个子区域中，若是，则进入步骤S305，若否，则返回步骤S110。

步骤S305，判断单元中的最小单元对应的迭代次数是否达至预设的迭代次数，若是，则进入步骤S306，若否，则返回步骤S10。

步骤S306，停止进行牛顿迭代。

步骤S50，将最终迭代得到的弹性模量值组成成像组织的弹性模量分布。

如图8所示，在一个实施例中，一种磁共振弹性成像中的弹性模量重建系统，包括有限体积元运算模块10、迭代模块30和分布形成模块50。

有限体积元运算模块10，用于将成像组织的表面假设为平面域，通过有限体积元算法求解得到平面域中的位移初值，以及平面域中子区域的位移、位移对未知弹性模量的导数。

本实施例中，在一定条件下，可将一个在外力作用下处于平衡状态的弹性体视为平面弹性问题，而在磁共振弹性成像中，成像组织被施加的外力是非常小的，可将该成像组织视为处于平衡状态的，因此，有限体积元运算模块10可将磁共振弹性成像中应力和应变之间的关系视为平面弹性问题。

有限体积元运算模块10将成像组织的表面视为平面域Ω，获取通过平面域对应的平衡方程所构建得到的二阶椭圆微分方程组，并在构建的二阶椭圆微分方程组中对未知弹性模量求导得到包含了位移对未知弹性模量的导数的方程。

进一步的，有限体积元运算模块10对于平面域Ω，其边界为描述成像组织平衡的状态变量包括三组，即，应力张量σ＝(σ₁₁,σ₂₂,σ₁₂)^T、应变张量ε＝(ε₁₁,ε₂₂,ε₁₂)^T和位移向量u＝(u₁,u₂)^T。假设Ω是均匀各向同性的弹性体，令

$\▿=(\frac{\∂}{\∂x_1},\frac{\∂}{\∂x_2})$

$B\▿=\left(\begin{array}{ccc}\frac{\∂}{\∂x_1}& 0& \frac{\∂}{\∂x_2}\\ 0& \frac{\∂}{\∂x_2}& \frac{\∂}{\∂x_1}\end{array}\right)$

$A=\left(\begin{array}{ccc}\mathrm{\λ}+2\mathrm{\μ}& \mathrm{\λ}& 0\\ \mathrm{\λ}& \mathrm{\λ}+2\mathrm{\μ}& 0\\ 0& 0& \mathrm{\μ}\end{array}\right)$

其中，正数λ和μ是Lamé常数：

$\mathrm{\λ}=\frac{\mathrm{vE}}{(1+v)(1-2v)},$ $\mathrm{\μ}=\frac{E}{2+2v}$

其中，ν是Poisson系数（Poisson’sRatio，是已知的常数），E是弹性模量（Young’sModulus）。σ，ε，u满足如下三组方程：

（应变-位移关系）（1）

（平衡方程）（2）

σ＝Aε（应力-应变关系）（3）

其中，f是体积力。

由Green公式可推出：

${\∫}_{\mathrm{\Ω}}{\mathrm{\σ}}^T{B}^T(\▿)\mathrm{udx}+{\∫}_{\mathrm{\Ω}}{\left(B(\▿)\mathrm{\σ}\right)}^T\mathrm{udx}={\∫}_{\∂\mathrm{\Ω}}{\left(B\left(\mathrm{\υ}\right)\mathrm{\σ}\right)}^T\mathrm{udx}$

其中，υ＝(υ₁,υ₂)^T是的单位外法向量。

假定Γ分为两段不重叠的线段Γ₀和Γ₁，在Γ₀上给定位移边界条件：在Γ₁上给定力条件（表面力）。

在求解时，利用公式（1）-（3）消去部分未知量，留下的是待求量，即弹性模量E，得到包含了位移向量u的二阶椭圆微分方程组：

$-\mathrm{\μ}\▿u-(\mathrm{\λ}+\mathrm{\μ})\mathrm{grad}\mathrm{div}u=f---\left(4\right).$

其中，以分别乘（4）式两端，关于x∈Ω积分，然后利用Green公式，得积分形式：

$a(u,v)-{\∫}_{\mathrm{\Ω}}[\mathrm{\μ}\frac{\∂u}{\∂\mathrm{\υ}}+(\mathrm{\λ}+\mathrm{\μ})\left(\mathrm{div}u\right)v]\mathrm{vds}=(f,v)---\left(5\right)$

其中

$a(u,v)={\∫}_{\mathrm{\Ω}}[\mathrm{\μ}\▿u\▿v+(\mathrm{\λ}+\mathrm{\μ})\mathrm{div}u\·\mathrm{div}v]\mathrm{dx}---\left(6\right)$

是边界上的法向导数，是Ω上的Hilbert空间，定义在其上的函数在平面域Ω的边界上的值为零.

在边界Γ₁上有

$\mathrm{\μ}\frac{\∂u}{\∂\mathrm{\υ}}+(\mathrm{\λ}+\mathrm{\μ})\mathrm{div}u={\left(B\left(\mathrm{\υ}\right)\mathrm{\σ}\right)}^T=\stackrel{\‾}{P}$

故方程（4）的变分形式是：求u∈(H¹(Ω))²，使得

$a(u,v)=(f,v)+{\∫}_{{\mathrm{\Γ}}_1}\stackrel{\‾}{P}\mathrm{vds},\∀v\∈{\left(H_0^1\left(\mathrm{\Ω}\right)\right)}^2---\left(7\right)$

其中，H¹(Ω)是平面域Ω上的Hilbert空间.

进一步的，为求解出位移对未知弹性模量的导数，有限体积元运算模块10将对（4）式两端对未知弹性模量求导，以得到包含了位移对未知弹性模量的导数的方程。

在构建得到了二阶椭圆微分方程组和包含了位移对未知弹性模量的导数的方程之后，有限体积元运算模块10将应用有限体积元算法求解位移和位移对未知弹性模量的导数。

有限体积元法是吸取了有限差分法和有限元法的优点发展起来的，又称广义差分法，可以处理复杂的边值条件及不规则区域，并且由于离散化得到的线性方程组是稀疏的，在不降低数值解的收敛阶的前提下求解方程的计算量小，计算速度快，因此，将有限体积元法应用于磁共振弹性成像中将大大地降低了弹性模量重建的计算量。

在一个实施例中，上述有限体积元运算模块10还用于根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。

本实施例中，有限体积元运算模块10将首先进行位移初值的求解，以方便后续的位移求解过程。在位移初值的求解过程中，将进行弹性模量初值的假设，并根据这一假设通过有限体积元算法得到平面域所对应的位移初值，其中，该位移初值计算时所使用的边界Γ上的初值取为对成像组织进行磁共振成像所得到的位移图中对应的位移值。

如图9所示，在一个实施例中，上述有限体积元运算模块10包括划分单元110、子区域位移求解单元130和子区域导数求解单元150。

划分单元110，用于将平面域划分为若干个子区域，并划分子区域为若干个单元，该子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上，平面域边界的顶点为单元的顶点。

本实施例中，划分单元110将平面域Ω划分为若干个子区域，然后分别划分每一子区域为若干个单元之和，在优选的实施例中，划分单元110划分得到的单元为三角形的形状，进而使得子区域之间、单元之间互相不存在重叠区域，并且任一单元的顶点不会在其它单元的边上，边界Γ的每一顶点均为单元的顶点，从而每一子区域将对应了一个原始剖分T_h，其中h是所有三角形单元边的最长边。记K_Q为T_h内的三角形单元，Q是三角形单元的外心。

进一步的，划分单元110记是原始剖分T_h的对偶剖分，包括是重心对偶剖分和外心对偶剖分。以下设是外心对偶剖分。外心对偶剖分做法：设T_h的任一单元的内角不大于90°，取三角单元□P₀P_iP_i+1(i＝1,2,…6,P₇＝P₁)的外心Q_i为对偶剖分的节点，即依次连结各个外心点可得到外心对偶剖分，此时，是的中垂线，分别过边中点M_i。

子区域位移求解单元130，用于在子区域中，以构成单元的顶点作为子区域的节点，并构建节点对应的有限体积元方程，并以位移初值作为有限体积元方程中的初值进行求解得到子区域对应的位移。

本实施例中，对于每一子区域，子区域位移求解单元130设是定义在原始分剖T_h上的试探函数空间，它是分片线性函数空间，即

为一次多项式，完全由三角形单元中三个顶点上的值所确定， $\∀K_Q\∈T_h,u_h{|}_{\mathrm{\Γ}0h}=u_0\}.$

子区域位移求解单元130将T_h的内节点编号为1,2,…,N₀。其中，节点分两类，给定力条件的节点编号为N₀+1,…,N₁，给定位移条件的节点编号为N₁+1,…,N。用表示节点i∈{1,2,…,N₁}的基函数，则对u_h∈U_h，可表示为：（Γ_0h是Γ₀的近似）

u_i是u_h(x)在第i个节点x_i的值。又假设V_h是定义在的检验函数空间，它为分片常数空间，即对任意的内节点P₀，相应的基函数为的特征函数：

${\mathrm{\ψ}}_{P_0}\left(P\right)=\left\{\begin{array}{c}1,P\∈K_{P_0}^{*}\\ 0,P\∉K_{P_0}^{*}\end{array}\right.$

设节点基函数为ψ_j＝(ψ_1j，ψ_2j)^T。则得到有限体积元法方程是：

求u_h∈U_h，使得

$a(u_h,{\mathrm{\ψ}}_j)=(f,{\mathrm{\ψ}}_j)+{\∫}_{{\mathrm{\Γ}}_{1h}}\stackrel{\‾}{P}{\mathrm{\ψ}}_j\mathrm{ds},j=\mathrm{1,2},...,N_1---\left(8\right)$

其中，Γ_1h是Γ₁的近似。

将式（4）两端在上积分，应用Green公式并用u_h代替u，得

$-{\∫}_{\∂K_{P_0}^{*}}[\mathrm{\μ}\frac{\∂u_h}{\∂\mathrm{\υ}}+(\mathrm{\λ}+\mathrm{\μ})\left(\mathrm{div}{u}_h\right)\mathrm{\υ}]\mathrm{ds}={\∫}_{K_{P_0}^{*}}\mathrm{fdx}---\left(9\right)$

此即为在节点P₀的有限体积元方程。下面给出（9）式左端积分的计算。第一项积分分解为中垂线段上的积分和。例如沿的积分为

$-\mathrm{\μ}{\∫}_{\stackrel{\‾}{Q_1{Q}_2}}\frac{\∂u_h}{\∂\mathrm{\υ}}\mathrm{ds}=-\mathrm{\μ}\stackrel{\‾}{Q_1{Q}_2}\frac{u_{P_2}-u_{P_0}}{\stackrel{\‾}{P_0{P}_2}}---\left(10\right)$

其余中垂线段上的积分类推。第二项积分分解为沿折线 上的积分。例如，沿上的积分

$-(\mathrm{\λ}+\mathrm{\μ}){\∫}_{\stackrel{\‾}{M_1{Q}_1{M}_2}}\left(\mathrm{div}{u}_h\right)\mathrm{\υds}$

$=-(\mathrm{\λ}+\mathrm{\μ})\mathrm{div}{u}_h\left(Q_1\right)(\frac{\stackrel{\→}{P_0{P}_1}}{\stackrel{\‾}{P_0{P}_1}}\stackrel{\‾}{M_1{Q}_1}+\frac{\stackrel{\→}{P_0{P}_2}}{\stackrel{\‾}{P_0{P}_2}}\stackrel{\‾}{Q_1{M}_2})---\left(11\right)$

其中

$\mathrm{div}{u}_h\left(Q_1\right)=\frac{\∂u_{1h}}{\∂x_1}\left(Q_1\right)+\frac{\∂u_{2h}}{\∂x_2}\left(Q_1\right)$

$=\frac{1}{2S_{Q_1}}[(x_2\left(P_1\right)-x_2\left(P_2\right))u_1\left(P_0\right)+(x_2\left(P_2\right)-x_2\left(P_0\right))u_1\left(P_2\right)+(x_2\left(P_0\right)-x_2\left(P_1\right))u_1\left(P_2\right)]$

$+\frac{1}{2S_{Q_1}}[(x_1\left(P_2\right)-x_1\left(P_1\right))u_2\left(P_0\right)+(x_1\left(P_0\right)-x_1\left(P_2\right))u_2\left(P_1\right)+(x_1\left(P_1\right)-x_1\left(P_0\right))u_2\left(P_2\right)]---\left(12\right)$

其中，(x₁(P_i),x₂(P_i))是点P_i的坐标，是含外心Q₁的单元面积。其余折线段上的积分类推。

方程（10）是内点有限体积元法方程。在边界Γ_0h上给定位移值u₀。例如假设节点此时在处仍有方程（9），左端的线积分仍按公式（10）和（11）计算。

总之子区域位移求解单元130可列出2N₁个形如（9）式的有限体积元方程和Γ_0h上的位移边界条件。也就是说，假设已知给定的表面力，和正数λ、μLamé常数，就可以计算出位移所对应的数值。

子区域导数求解单元150，用于通过求解得到的位移对包含了位移对未知弹性模量的导数的方程进行求解得到位移对未知弹性模量的导数，位移对未知弹性模量的导数是与位移所在的子区域相对应的。

本实施例中，子区域导数求解单元150在式（4）中对每一子区域的未知弹性模量求导数，即可得到慎勿使形谍了位移对未知弹性模量的导数的方程，即：

$-\mathrm{\μ}\▿u_j^{\′}-(\mathrm{\λ}+u)\mathrm{grad}\mathrm{div}{u}_j^{\′}=\frac{\∂\mathrm{\μ}}{\∂E_{z_j}}\▿u+\frac{\∂(\mathrm{\λ}+\mathrm{\μ})}{\∂E_{z_j}}\mathrm{grad}\mathrm{div}u---\left(13\right)$

其中z用于子区域z，（13）式的形式如同（4），这里是关于u′_j的方程，方程右端可看为是（4）式中的f，因此可以用解（4）式的有限体积元法解（13）式，得到u′_j的解，在此不再赘述。

迭代模块30，用于通过有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代。

本实施例中，迭代模块30应用牛顿迭代法进行弹性模量值的，假设迭代第n步的解为则第n+1步的解为

$E_z^{(n+1)}=E_z^{\left(n\right)}+\mathrm{\Δ}{E}_z---\left(14\right)$

其中ΔE_z是弹性模量E_z的改进值，它是如下正则化方程组的解

$\left[(H_z^{\left(n\right)}+\mathrm{\αI})\right]\mathrm{\Δ}{E}_z=-f_z^{\left(n\right)}---\left(15\right)$

是不计二阶导数项的海森（Hessian）矩阵，元素为

$h_{\mathrm{ij}}^{\left(n\right)}=\frac{\∂f_i^{\left(n\right)}}{\∂E_{z_j}^{\left(n\right)}}=-\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_i}^{\left(n\right)}}\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_j}^{\left(n\right)}}+\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_i}^{\left(n\right)}}\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_j}^{\left(n\right)}})---\left(16\right)$

α是正则化对角矩阵的参数，使可逆，I是单位矩阵。

因此，由上述推导可知，迭代模块30在进行牛顿迭代计算弹性模量值时，将通过有限体积元算法求解得到的位移和导数经由式（15）和式（16）计算得到改进值，进而方可进行下一步的迭代。

如图10所示，在一个实施例中，上述迭代模块30包括改进值运算单元310和弹性模量迭代单元330。

改进值运算单元310，用于通过每一子区域所对应的位移和导数进行运算得到子区域对应的弹性模量改进值。

弹性模量迭代单元330，用于根据弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前失效次数对应的弹性模量值。

本实施例中，弹性模量迭代单元330通过计算得到的弹性模量改进值和上一次迭代所得到的弹性模量值通过式（14）进行计算得到当前迭代次数所对应的弹性模量值。

如图11所示，在一个实施例中，上述迭代模块30还包括位移图像处理单元301和判断单元303。

位移图处理单元301，用于获取通过对成像组织进行磁共振成像所得到的位移图，并根据当前迭代得到的弹性模量值和位移图得到最优化问题方差。

本实施例中，由于磁共振成像中的弹性模量重建问题被归结为一个带约束条件的最优化问题，其目标函数即为求解得到磁共振成像所测得的位移与计算得到的位移之间的最小平方和差。

该最优化问题为：minF(E)

其中

$F\left(E\right)=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}[(u_i^{\mathrm{meas}}-u_i^{\mathrm{comp}})+(v_i^{\mathrm{meas}}-v_i^{\mathrm{comp}})]$

这里和分别是在成像组织位置i处通过磁共振测得的x方向和y方向的位移，和在成像组织位置i处用有限体积元法计算出的位移，一共有N个不同的位置。E是一个M维的弹性参数向量，它由一个连续基底φ张成，用它来定义整个感兴趣区域的弹性组织分布。

位移图处理单元301假设整个感兴趣区域为Ω，将Ω分为若干个子区域Ω_sub，可将最优化问题重写为子区域上的最优化问题。假设一共有Q个子区域，有

$F\left(E\right)=\underset{z=1}{\overset{Q}{\mathrm{\Σ}}}{F}_z\left(E_z\right)$

其中，F_z(E_z)是第z个子区域上的最优化函数，对所有的子区域求和可得最优化问题为：

$\min F\left(E\right)=\min \left[\underset{z=1}{\overset{Q}{\mathrm{\Σ}}}{F}_z\left(E_z\right)\right]\⇒\underset{z=1}{\overset{Q}{\mathrm{\Σ}}}\min F_z\left(E_z\right)$

这里

$F_z\left(E_z\right)=\underset{i_Z}{\overset{N_z}{\mathrm{\Σ}}}[{(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})}^2+{(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})}^2]$

这里和分别是在成像组织第z个子区域上i_z处磁共振测得的x方向和y方向的位移，和在成像组织z个子区域上i_z处用有限体积元法计算出的位移，一共有N_z个不同的位置。E_z是M_z的向量，有N_z＜N，M_z＜M。对于在方程两边分别对E_z求一阶导数，并令导数等于0，可得如下非线性方程组：

$\begin{array}{c}{f}_1=\frac{\∂F_z}{\∂E_{z_1}}=\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_1}}+\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_1}}=0\\ f_2=\frac{\∂F_z}{\∂E_{z_2}}=\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_2}}+\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_2}}=0\\ \begin{array}{ccccc}.& & .& & .\\ .& & .& & .\\ .& & .& & .\end{array}\\ f_{M_z}=\frac{\∂F_z}{\∂E_{z_{M_z}}}=\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(u_{i_z}^{\mathrm{meas}}-u_{i_z}^{\mathrm{comp}})\frac{\∂u_{i_z}^{\mathrm{comp}}}{\∂E_{z_{M_z}}}+\underset{i_z=1}{\overset{N_z}{\mathrm{\Σ}}}(v_{i_z}^{\mathrm{meas}}-v_{i_z}^{\mathrm{comp}})\frac{\∂v_{i_z}^{\mathrm{comp}}}{\∂E_{z_{M_z}}}=0\\ \end{array}$

进而位移图处理单元301通过如上所述的牛顿迭代法求解上述非线性方程组。

判断单元303，用于判断最优化问题平方差是否小于预设的容忍误差，若是，则进一步判断单元中的最小单元对应的迭代次数是否达到预设的迭代次数，若是，则停止进行牛顿迭代。

本实施例中，预设的容忍误差是根据实际需要进行灵活设定的。

在另一个实施例中，上述判断单元303还用于若判断到单元中的最小单元对应的迭代次数未达到预设的迭代次数，则通知上述有限体积元运算模块10。

分布形成模块50，用于将最终迭代得到的弹性模量值组成成像组织的弹性模量分布。

上述磁共振弹性成像中的弹性模量重建方法和系统，将成像组织的表面假设为平面域，引入有限体积元算法以求解得到位移和该位移对未知弹性模量的导数，进而通过求解得到的位移和导数进行牛顿迭代得到弹性模量值，并在弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代，进而将最终迭代得到的弹性模量值组成成像组的弹性模量分布，通过引入有限体积元算法，使得弹性模量的重建在保证了计算精度的前提下降低了计算量，提高了计算速度。

本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的程序可存储于一计算机可读取存储介质中，该程序在执行时，可包括如上述各方法的实施例的流程。其中，所述的存储介质可为磁碟、光盘、只读存储记忆体（Read-OnlyMemory，ROM）或随机存储记忆体（RandomAccessMemory，RAM）等。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (10)

1.一种磁共振弹性成像中的弹性模量重建方法，包括如下步骤：

将成像组织的表面假设为平面域，通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数；

通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代；

将最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布；其中，通过有限体积元算法求解得到所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤包括：

将所述平面域划分为若干个子区域，并划分所述子区域为若干个单元，所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上，平面域边界的顶点为单元的顶点；

在所述子区域中，以构成单元的顶点作为所述子区域的节点，并构建所述节点对应的有限体积元方程，并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移；

通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数，所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。

2.根据权利要求1所述的磁共振弹性成像中的弹性模量重建方法，其特征在于，所述通过有限体积元算法求解得到所述平面域中的位移初值的步骤为：

根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。

3.根据权利要求1所述的磁共振弹性成像中的弹性模量重建方法，其特征在于，所述通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值的步骤包括：

通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值；

根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。

4.根据权利要求3所述的磁共振弹性成像中的弹性模量重建方法，其特征在于，所述根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值的步骤之后还包括：

获取通过对所述成像组织进行磁共振成像所得到的位移图；

根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差；

判断所述最优化问题平方差是否小于预设的容忍误差，若是，则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中，若是，则

判断所述单元中的最小单元对应的迭代次数是否达至预设的迭代次数，若是，则

停止进行牛顿迭代。

5.根据权利要求4所述的磁共振弹性成像中的弹性模量重建方法，其特征在于，还包括：

若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数，则

返回所述通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数的步骤。

6.一种磁共振弹性成像中的弹性模量重建系统，其特征在于，包括：

有限体积元运算模块，用于将成像组织的表面假设为平面域，通过有限体积元算法求解得到所述平面域中的位移初值，以及所述平面域中子区域的位移、所述位移对未知弹性模量的导数；

迭代模块，用于通过所述有限体积元算法求解得到的位移和导数进行牛顿迭代得到弹性模量值，直至所述弹性模量值对应的最优化问题平方差小于预设的容忍误差，且达到预设的迭代次数时停止牛顿迭代；

分布形成模块，用于将最终迭代得到的弹性模量值组成所述成像组织的弹性模量分布；其中，

所述有限体积元运算模块包括：

划分单元，用于将所述平面域划分为若干个子区域，并划分所述子区域为若干个单元，所述子区域之间和单元之间不存在重叠且任一单元的顶点均不在其它单元的边上，平面域边界的顶点为单元的顶点；

子区域位移求解单元，用于在所述子区域中，以构成单元的顶点作为所述子区域的节点，并构建所述节点对应的有限体积元方程，并以所述位移初值作为有限体积元方程中的初值进行求解得到所述子区域对应的位移；

子区域导数求解单元，用于通过所述求解得到的位移对包含了所述位移对未知弹性模量的导数的方程进行求解得到所述位移对未知弹性模量的导数，所述位移对未知弹性模量的导数是与位移所在的子区域相对应的。

7.根据权利要求6所述的磁共振弹性成像中的弹性模量重建系统，其特征在于，所述有限体积元运算模块还用于根据设定的弹性模量初值通过有限体积元算法计算得到平面域中的位移初值。

8.根据权利要求6所述的磁共振弹性成像中的弹性模量重建系统，其特征在于，所述迭代模块包括：

改进值运算单元，用于通过每一子区域所对应的位移和导数进行运算得到所述子区域对应的弹性模量改进值；

弹性模量迭代单元，用于根据所述弹性模量改进值以设定的弹性模量初值为起始进行牛顿迭代得到与当前迭代次数对应的弹性模量值。

9.根据权利要求8所述的磁共振弹性成像中的弹性模量重建系统，其特征在于，所述迭代模块还包括：

位移图处理单元，用于获取通过对所述成像组织进行磁共振成像所得到的位移图，并根据当前迭代得到的弹性模量值和所述位移图得到最优化问题平方差；

判断单元，用于判断所述最优化问题平方差是否小于预设的容忍误差，若是，则进一步判断所述平面域中的每一单元是否均包含于至少一个子区域中，若是，则判断所述单元中的最小单元对应的迭代次数是否达到预设的迭代次数，若是，则停止进行牛顿迭代。

10.根据权利要求9所述的磁共振弹性成像中的弹性模量重建系统，其特征在于，所述判断单元还用于若判断到所述单元中的最小单元对应的迭代次数未达到预设的迭代次数，则通知所述有限体积元运算模块。