CN102920457A - 磁共振弹性成像精确度检测方法 - Google Patents

Abstract

一种磁共振弹性成像精确度检测方法，包括以下步骤：获取在预定外力作用下，通过磁共振弹性成像所得的预定成像区域的成像位移值；获取所述预定外力及所述预定成像区域的弹性模量的大小；根据所述预定外力大小，得到所述预定成像区域的应力大小；由所述预定成像区域的应力及弹性模量的大小，得到所述预定成像区域的理论位移值；将所述预定成像区域的理论位移值与通过磁共振弹性成像所得的成像位移值进行比较，检测出磁共振弹性成像的精确度。上述磁共振弹性成像精确度检测方法采用有限体积元法，由预定成像区域的应力及弹性模量的大小，得到理论位移值，为在磁共振弹性成像获取位移图后，逆计算求解成像区域的弹性模量提供了参考。

Description

磁共振弹性成像精确度检测方法

技术领域

本发明涉及磁共振技术，特别是涉及一种磁共振弹性成像精确度检测方法。

背景技术

弹性（stiffness）是人体组织物理性质中一种重要的机械力学参数。生物组织的弹性变化常与病理现象紧密相关，病变组织和正常组织往往存在弹性模量或硬度的差异，这种差异为临床上疾病的诊断提供了重要的参考信息。磁共振弹性成像（Magnetic Resonance Elastography，MRE）作为一种新型的无创成像方法，能直观显示和量化人体内部组织弹性，并对组织的弹性成像，使“影像触诊”成为了可能，在乳腺癌检测、肝硬化分期，动脉粥样硬化斑块、肌肉损伤、大脑疾病检测、射频消融等治疗和监控方面具有重要意义。

假设已知物体的弹性模量和施加在物体上的应力张量，在一定条件下，可将一在外力作用下处于平衡状态的弹性体视作平面弹性问题。磁共振弹性成像技术的基本原理是利用磁共振成像技术（Magnetic Resonance Imaging，MRI），检测体内组织在某种外力作用下产生的质点位移，以此为基础通过对弹性力学的逆求解，得出组织内各点的弹性系数的分布图。传统的磁共振弹性成像技术中，缺乏对其精确度进行检验的手段，导致成像精确度不高。

发明内容

基于此，有必要提供一种能检测磁共振弹性成像精确度的磁共振弹性成像精确度检测方法。

一种磁共振弹性成像精确度检测方法，包括以下步骤：

获取在预定外力作用下，通过磁共振弹性成像所得的预定成像区域的成像位移值；

获取所述预定外力及所述预定成像区域的弹性模量的大小；

根据所述预定外力大小，得到所述预定成像区域的应力大小；

由所述预定成像区域的应力及弹性模量的大小，得到所述预定成像区域的理论位移值；

将所述预定成像区域的理论位移值与通过磁共振弹性成像所得的成像位移值进行比较，检测出磁共振弹性成像的精确度。

在其中一个实施例中，所述由所述预定成像区域的应力及弹性模量的大小，得到所述预定成像区域的理论位移值的步骤具体为：由所述预定成像区域的应力及弹性模量的大小，经有限体积元法，并通过弹性力学方程得到所述预定成像区域的理论位移值。

在其中一个实施例中，所述由所述预定成像区域的应力及弹性模量的大小，经有限体积元法，并通过弹性力学方程得到所述预定成像区域的理论位移值的步骤具体包括以下步骤：

将所述预定成像区域进行三角剖分，得到原始单元；

对三角剖分后的所述预定成像区域进行对偶剖分，得到多个对偶单元；

对三角剖分后的所述预定成像区域建立试探函数，对对偶剖分后的所述预定成像区域建立检验函数；

由试探函数及检验函数得到有限体积元格式的弹性力学方程；

对每一个对偶单元求解所述有限体积元格式的弹性力学方程，得到所述预定成像区域的理论位移值。

在其中一个实施例中，所述对偶剖分为外心对偶剖分。

在其中一个实施例中，所述试探函数为：

其中，u_h为原始单元上的试探函数，u_i是u_i(x)在第i个节点x_i的值。

在其中一个实施例中，所述检验函数为：

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

其中，P为原始单元的节点，

为P₀上的特征函数。

在其中一个实施例中，所述有限体积元格式的弹性力学方程为：

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

其中，f为应力，ψ_j为单元上的检验函数，u_h为原始单元上的试探函数。

上述磁共振弹性成像精确度检测方法，将其所得到的预定成像区域的理论位移值与通过对预定成像区域进行磁共振弹性成像所得的成像位移值进行比较，如二者之间存在明显误差，则表明磁共振弹性成像系统的精确度还需提高，通过对磁共振弹性成像系统的调节，并通过上述磁共振弹性成像精确度检测方法来检测其误差信息，可以逐步提高磁共振弹性成像系统成像的精确度。同时，上述磁共振弹性成像精确度检测方法采用有限体积元法，由预定成像区域的应力及弹性模量的大小，得到理论位移值，为在磁共振弹性成像获取位移图后，逆计算求解成像区域的弹性模量提供了参考。

附图说明

图1为一实施例的磁共振弹性成像精确度检测方法的流程图；

图2为磁共振弹性成像精确度检测方法的具体流程图；

图3为对预定成像区域进行三角剖分的示意图；

图4为对三角剖分后的预定成像区域进行的对偶剖分的示意图。

具体实施方式

为了解决传统的磁共振弹性成像技术中，缺乏对其精确度进行检验的手段，导致成像精确度不高的问题，提出一种能检测磁共振弹性成像精确度的磁共振弹性成像精确度检测方法。

在磁共振弹性成像时，在被测物体的预定成像区域上施加的预定外力进行激励，由于预定外力较小，可视该预定成像区域仍处于平衡状态，故磁共振弹性成像中应力与应变力之间的关系可视为平面弹性问题。

请参阅图1，本实施例的磁共振弹性成像精确度检测方法，包括以下步骤：

步骤S100，获取在预定外力作用下，通过磁共振弹性成像所得的预定成像区域的成像位移值。

对弹性模量大小已知的预定成像区域上施加的预定外力进行激励，通过磁共振弹性成像，得到预定成像区域上的各点的成像位移值，进而可以通过该成像位移值及所施加的预定外力获得预定成像区域内各点的弹性系数分布图。获取上述通过磁共振弹性成像所得的预定成像区域的成像位移值。

步骤S200，获取预定外力及预定成像区域的弹性模量的大小。获取预定成像区域的弹性模量的大小，及施加在预定成像区域上的预定外力大小。

步骤S300，根据预定外力大小，得到预定成像区域的应力大小。由于在外力作用下，处于平衡状态的弹性体可视作平面弹性问题，预定外力与预定成像区域的应力相平衡，即可由预定外力大小得到预定成像区域的应力大小。

步骤S400，由预定成像区域的应力及弹性模量的大小，得到预定成像区域的理论位移值。由预定成像区域的应力及弹性模量的大小，即可通过弹性力学方程，得到预定成像区域的理论位移值。由于在外力作用下，处于平衡状态的弹性体可视作平面弹性问题，弹性力学方程为偏微分方程的一种，可以通过有限差分法、有限元法等方法对其进行求解。

假设预定成像区域的表面是一个平面域Ω，是其的边界，刻画弹性体平衡的状态变量有三组，即应力张量σ=(σ₁₁,σ₂₂,σ₁₂)^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{\νE}}{(1+\mathrm{\ν})(1-2\mathrm{\ν})}$

$\mathrm{\μ}=\frac{E}{2+2\mathrm{\ν}}$

这里v是泊松系数（Poisson’s Ratio），为指在固体力学中，材料的横向变形系数。E是杨氏模数（Young’s Modulus），即材料在弹性变形范围内，正应力与正应变的比值。λ和μ即可表征预定成像区域的弹性模量。σ、ε、u满足如下三组方程：

$\mathrm{\ϵ}=B^T(\▿)u---\left(1a\right)$

$B(\▿)\mathrm{\σ}+f=0---\left(1b\right)$

σ=Aε（1c）

其中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是

的单位外法向量。

假定Γ分为两段不重叠的线段Γ₀和Γ₁，在线段Γ₀上的位移边界条件：

在Γ₁上给定力条件

其中

为表面力。实际求解时，可利用公式（1a）与公式（1c）消去σ和ε，得出位移向量u满足的二阶椭圆微分方程组，即弹性力学方程：

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

具体的，可通过有限差分法、有限元法或有限体积元法等方法，由预定成像区域的应力f及弹性模量的大小λ及μ，求解公式（2）中的位移向量u，求解得的位移向量u即为预定成像区域的理论位移值。在本实施例中，具体采用有限体积元法。

有限体积元法又将它称为广义差分法（generalized difference method，GDM），它具有可以处理复杂的边值条件及不规则区域及离散化得到稀疏的线性方程组，在不降低数值解的收敛阶的前提下求解方程的计算量小、计算速度快的优点，可以很好的处理磁共振弹性成像中的平衡方程的复杂边值条件和组织的不规则性等情况。

是Ω上的Hilbert空间，定义在其上的函数在Ω的边界上的值为零。以

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

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

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

在边界Γ₁上满足：

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

故方程（2）的变分形式是：

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

其中，u∈(H¹(Ω))²，H¹(Ω)是Ω上的Hilbert空间，

要将公式（5）构造成有限体积元格式，分别需要构造原始剖分、定义于原始剖分上的试探函数空间及对偶剖分、定义于对偶剖分上的检验函数空间。

请参阅图2，步骤S400具体包括以下步骤：

步骤S410，将预定成像区域进行三角剖分，得到原始单元。

请一并参阅图3，对预定成像区域，即平面域Ω进行三角剖分，将平面域Ω划分为多个的三角形单元之和，使得多个的三角形单元之间互相没有重叠且任意一个三角形单元的顶点都不在其它任何一个三角单元的边上，另外边界Γ上的每一个顶点都是三角形的顶点，这样可以得到一个由多个三角形单元组成的原始剖分T_h，其中h是所有三角形单元边的最长边。记K_Q为原始剖分T_h内的三角形单元，记为原始单元。

步骤S420，对三角剖分后的预定成像区域进行对偶剖分，得到多个对偶单元。对偶剖分分为外心对偶剖分及重心对偶剖分等，具体在本实施例中，采用外心对偶剖分的方式，对三角剖分后的预定成像区域进行对偶剖分。请一并参阅图4，Q₁~Q₆是原始剖分T_h中以P₀为顶点的所有三角单元的外心，且原始剖分T_h中任意一个三角形单元的内角都不大于90°，取三角单元△P₀P_iP_i+1(i=1,2,…6,P₇=P₁)的外心Q_i为对偶剖分的节点，即依次连结Q₁~Q₆可得到阴影部分即为以P₀为顶点的外心对偶剖分单元，此时

是

的中垂线，并分别过各边的中点M_i。

步骤S430，对三角剖分后的预定成像区域建立试探函数，对对偶剖分后的预定成像区域建立检验函数。

设

是定义在三角剖分后的预定成像区域，即原始剖分T_h上的试探函数空间。是分片线性函数空间，即：

$U_h=\left\{u_h\right|u_h\∈C\left(\mathrm{\Ω}\right),u_h{|}_{K_Q},\∀K_Q\∈T_h,u_h{|}_{{\mathrm{\Γ}}_{0h}}=u_0\}$

上式为一次多项式，C(Ω)为Ω上的连续函数空间完全由三角形单元K_Q三个顶点上的值所确定。在三角形单元K_Q中，设其位于原始剖分T_h的内部的节点编号为1，2，…，N₀。而在原始剖分T_h上的界点分两类，给定力条件的界点编号为N₀+1，…，N₁，给定位移条件的界点编号为N₁+1，…，N。用

表示内节点i∈{1,2,…,N₁}的基函数，Γ_0h是Γ₀的近似。则对u_h∈U_h，

可表示为：

上式即三角剖分后的预定成像区域的试探函数，其中，u_h为原始单元上的试探函数，u_i是u_h(x)在第i个节点x_i的值。

设V_h为定义在对偶剖分后的预定成像区域

的检验函数空间，它为分片常数空间。对任意的内节点，如图3中的任意的一个内部节点P₀，其相应的基函数为

的特征函数，即检验函数为：

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

其中，P为原始单元的节点，

为P₀上的特征函数。

步骤S440，由试探函数及检验函数得到有限体积元格式的弹性力学方程。

由上述试探函数及检验函数构造并得到有限体积元格式的弹性力学方程。具体的，设节点基函数为ψ_j=(ψ_1j,ψ_2j)^T，线段Γ_1h是线段Γ₁的近似。则基于位移法的有限体积元法方程为：

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

其中，u_h∈U_h，

f为应力，ψ_j为单元上的检验函数，u_h为原始单元上的试探函数。

步骤S450，对每一个对偶单元求解有限体积元格式的弹性力学方程，得到预定成像区域的理论位移值。

将公式（2）两端在对偶单元上积分，应用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(7\right)$

请再次参阅图3及图4，公式（7）即为在节点P₀的有限体积元方程。假定P₀及其相邻节点的分布如图4所示。对公式（7）左端积分进行计算，左端第一项积分可分解为中垂线段

上的积分和。例如沿

的积分为：

$-\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(8\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(9\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)]$

这里(x₁(Pi),x₂(P_i))是点P_i的坐标，

是含外心Q₁的三角单元面积。剩余的

折线段上的积分类推。

式（9）是内部节点的有限体积元法方程。Γ_0h是Ω的边界Γ₀的近似，Γ_1h是Ω的边界Γ₁的近似，对于边界Γ_0h及Γ_1h上在边界上的界点，Γ_0h上给定位移值u₀。而在属于Γ_1h的界点需建立补充方程，例如假设界点

此时在

处仍适用于式（7），式（7）左端的线积分仍可按式（8）及式（9）进行计算。

综上，由预定成像区域的应力及正数λ和μLamé常数，即弹性模量的大小，就可以通过有限体积元法得到预定成像区域的理论位移值。

步骤S500，将预定成像区域的理论位移值与通过磁共振弹性成像所得的成像位移值进行比较，检测出磁共振弹性成像的精确度。

将计算得到的预定成像区域的理论位移值与通过对预定成像区域进行磁共振弹性成像所得的成像位移值进行比较，如二者之间存在明显误差，则表明磁共振弹性成像系统的精确度还需提高，通过对磁共振弹性成像系统的调节，并通过上述磁共振弹性成像精确度检测方法来检测其误差信息，可以逐步提高磁共振弹性成像系统成像的精确度。同时，上述磁共振弹性成像精确度检测方法采用有限体积元法，由预定成像区域的应力及弹性模量的大小，得到理论位移值，为在磁共振弹性成像获取位移图后，对弹性力学方程进行逆计算求解成像区域的弹性模量提供了参考。

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

Claims (7)

1.一种磁共振弹性成像精确度检测方法，其特征在于，包括以下步骤：

获取在预定外力作用下，通过磁共振弹性成像所得的预定成像区域的成像位移值；

获取所述预定外力及所述预定成像区域的弹性模量的大小；

根据所述预定外力大小，得到所述预定成像区域的应力大小；

由所述预定成像区域的应力及弹性模量的大小，得到所述预定成像区域的理论位移值；

将所述预定成像区域的理论位移值与通过磁共振弹性成像所得的成像位移值进行比较，检测出磁共振弹性成像的精确度。

2.根据权利要求1所述的磁共振弹性成像精确度检测方法，其特征在于，所述由所述预定成像区域的应力及弹性模量的大小，得到所述预定成像区域的理论位移值的步骤具体为：由所述预定成像区域的应力及弹性模量的大小，经有限体积元法，并通过弹性力学方程得到所述预定成像区域的理论位移值。

3.根据权利要求2所述的磁共振弹性成像精确度检测方法，其特征在于，所述由所述预定成像区域的应力及弹性模量的大小，经有限体积元法，并通过弹性力学方程得到所述预定成像区域的理论位移值的步骤具体包括以下步骤：

将所述预定成像区域进行三角剖分，得到原始单元；

对三角剖分后的所述预定成像区域进行对偶剖分，得到多个对偶单元；

对三角剖分后的所述预定成像区域建立试探函数，对对偶剖分后的所述预定成像区域建立检验函数；

由试探函数及检验函数得到有限体积元格式的弹性力学方程；

对每一个对偶单元求解所述有限体积元格式的弹性力学方程，得到所述预定成像区域的理论位移值。

4.根据权利要求3所述的磁共振弹性成像精确度检测方法，其特征在于，所述对偶剖分为外心对偶剖分。

5.根据权利要求3所述的磁共振弹性成像精确度检测方法，其特征在于，所述试探函数为：

其中，u_h为原始单元上的试探函数，u_i是u_i(x)在第i个节点x_i的值。

6.根据权利要求3所述的磁共振弹性成像精确度检测方法，其特征在于，所述检验函数为：

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

其中，P为原始单元的节点，

为P₀上的特征函数。

7.根据权利要求3所述的磁共振弹性成像精确度检测方法，其特征在于，所述有限体积元格式的弹性力学方程为：

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

其中，f为应力，ψ_j为单元上的检验函数，u_h为原始单元上的试探函数。

Images