确定血管壁的弹性特征的方法、计算设备和介质

摘要

本公开提供了一种确定血管壁的弹性特征的方法、计算设备和计算机可读存储介质。该方法包括：基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相；基于所述多维动态血管壁模型和所述目标对象的多维血管影像的位移场信息，确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量；基于所述多维动态血管壁模型和所述目标对象的血压参数，确定所述血管壁的膨胀阶段的每个时相的血管壁应力张量；基于所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量分析所述血管壁的拟合曲线；以及基于所述血管壁的拟合曲线和预定应变幅值，确定所述血管壁的弹性值。

权利要求

1.一种确定血管壁的弹性特征的方法，包括：

基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相；

基于所述多维动态血管壁模型和所述目标对象的多维血管影像的位移场信息，确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量，其中所述位移场信息用于对所述目标对象的所述多维血管影像进行配准；

基于所述多维动态血管壁模型和所述目标对象的血压参数，确定所述血管壁的膨胀阶段的每个时相的血管壁应力张量；

基于所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量获取所述血管壁的拟合曲线；以及基于所述血管壁的拟合曲线和预定应变幅值，确定所述血管壁的弹性值，其中确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量包括：基于所述血管壁的膨胀阶段的初始时相的模型坐标和所述目标对象的多维血管影像的位移场信息确定所述膨胀阶段的每个时相的模型坐标；

确定所述膨胀阶段的每个时相的模型坐标相对于所述膨胀阶段的初始时相的模型坐标的变形梯度；

基于所述变形梯度确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量，其中确定所述血管壁的膨胀阶段的每个时相的血管壁应力张量包括：为血管壁模型选择材料模型，并设置所述血管壁模型的材料参数以使得模拟出的血管壁接近刚性壁面；

基于所述血管壁模型确定所述血管壁的能量方程；

将所述血管壁的膨胀阶段的每个时相对应的血压数据施加到所述血管壁模型中，并且利用反向有限元分析对所述血管壁的能量方程进行求解以确定所述血管壁应力张量，所述方法还包括基于所述目标对象的多维血管影像确定所述目标对象的多维动态血管壁模型，其包括：采集所述目标对象的多个时相的多维血管影像；

从所述多个时相中选择一个参考时相；

基于所述参考时相的多维血管影像生成所述目标对象的参考模型；以及基于对所述多维血管影像进行图像配准得到的位移场信息和所述参考模型来确定所述目标对象的多维动态模型。

2.如权利要求1所述的方法，其中基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相包括：基于在拍摄所述目标对象的多维血管影像时同步记录的心电门控心电图信号确定所述血管壁的膨胀阶段。

3.如权利要求1所述的方法，其中基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相包括：基于所述目标对象的多维动态血管壁模型确定在不同时相下所述血管壁的体积容量变化曲线；

确定所述血管壁的体积容量变化曲线中体积单调递增的阶段作为所述膨胀阶段。

4.如权利要求1所述的方法，其中所述血管壁应变张量包括Cauchy‑Green应变张量或Green‑Lagrange应变张量。

5. 如权利要求1所述的方法，其中所述材料模型包括超弹性材料模型，并且其中所述血管壁应力张量为Cauchy应力或Second PK应力。

6.如权利要求1所述的方法，其中对所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量进行拟合以获取所述血管壁的拟合曲线包括：基于所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量确定用于拟合的血管壁等效应变和血管壁等效应力；

为所述血管壁选择符合生物超弹性材料特性的拟合公式；以及基于所述血管壁等效应变和所述血管壁等效应力对所述拟合公式进行求解以确定拟合系数，从而确定所述血管壁的拟合曲线。

7.如权利要求1所述的方法，其中基于所述血管壁的拟合曲线和预定应变幅值，确定所述血管壁的弹性值包括：基于所述血管壁的拟合曲线求取所述血管壁等效应力相对于所述血管壁等效应变的一阶偏导数；以及将所述一阶偏导数与所述预定应变幅值之商确定为所述血管壁的弹性值。

8.如权利要求7所述的方法，还包括：

基于所述血管壁的拟合曲线求取所述血管壁等效应力相对于所述血管壁等效应变的二阶偏导数；

将所述一阶偏导数与所述预定应变幅值之商以及将所述二阶偏导数与所述预定应变幅值之商确定为所述血管壁的弹性值。

9.如权利要求1所述的方法，还包括：

基于所述目标对象的多维血管影像确定所述目标对象的位移场信息。

10. 一种计算设备，包括：

至少一个处理器；以及

至少一个存储器，所述至少一个存储器被耦合到所述至少一个处理器并且存储用于由所述至少一个处理器执行的指令，所述指令当由所述至少一个处理器执行时，使得所述计算设备执行根据权利要求1至9中任一项所述的方法的步骤。

11.一种计算机可读存储介质，其上存储有计算机程序代码，所述计算机程序代码在被运行时执行如权利要求1至9中任一项所述的方法。

说明书

确定血管壁的弹性特征的方法、计算设备和介质

技术领域

[0001] 本发明概括而言涉及医学图像处理领域，更具体地，涉及一种确定血管壁的弹性特征的方法、计算设备和计算机可读存储介质。

背景技术

[0002] 主动脉瘤是指主动脉病理性的扩张，超过正常血管直径的50%。当前，主动脉瘤的治疗方式包括开放手术和介入治疗。其中，开放手术是指通过开腹或开胸进行动脉瘤治疗，属于传统的治疗方法，存在着手术创伤大、风险高、对患者身体条件要求高的缺点。介入治疗是沿外周动脉将支架推送至主动脉瘤病灶并放置，从而将主动脉瘤与血流进行隔绝的治疗方法。介入治疗不需要打开胸腔或腹腔，具有微创、安全、操作简便及近期疗效确切等优点。

[0003] 当前，主动脉瘤的介入治疗的时机通常根据瘤径这一形态学参数做为临床指标。然而，瘤径无法反应患者个性化的差异，例如拥有相同瘤径的不同患者，该指标无法区分其病情的严重程度。从力学角度理解，动脉瘤的破裂属于一种结构失效行为，即血压引起的血管壁应力超出血管壁的强度极限进而造成组织破裂。因此，确定血管壁的力学特征参数对于辅助临床诊断和制定治疗方案具有重要意义。

[0004] 早期获取血管壁力学特征参数的技术手段主要为体外拉伸试验。尽管体外拉伸试验获取的力学特征参数（如拉伸模量等）有助于理解病变血管壁力学状态的变化，但由于其需要通过开放手术获取活体组织样本，因此无法辅助术前的病情诊断和治疗方案制定。

[0005] 随着动态医学影像技术的应用，如多维动态CT（计算机断层扫描，Computed Tomography）、4D核磁成像等，为进行血管壁在体应力或应变分析提供了可能。血管壁应变的在体分析技术主要基于医学影像获取的血管壁形态开展分析，然而血管壁应变与患者采集图像时的血压高度相关，且单一的应变所提供的信息有限。

[0006] 当前血管壁应力的分析技术，主要利用假设的血管壁本构模型，结合患者的血压信息，借助有限元仿真技术，对血管壁应力进行仿真分析。但是，血管壁组织具有高度各项异性、各项异质的特点，且不同患者之间的差异巨大，使用统一的本构模型假设会使应力分析结果引入不可预料的误差。

发明内容

[0007] 针对上述问题中的至少一个，本发明提供了一种确定血管壁的弹性特征的方法，其通过引入新的血管壁弹性指数来反应血管壁的硬化病变程度，在不需要假设无应力状态，也不需要假设材料本构模型的情况下，构建了一套基于多维动态血管壁影像的在体无创血管壁力学特征参数识别方法。

[0008] 根据本发明的一个方面，提供了一种确定血管壁的弹性特征的方法。该方法包括：基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相；基于所述多维动态血管壁模型和所述目标对象的多维血管影像的位移场信息，确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量；基于所述多维动态血管壁模型和所述目标对象的血压参数，确定所述血管壁的膨胀阶段的每个时相的血管壁应力张量；基于所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量获取所述血管壁的拟合曲线；以及基于所述血管壁的拟合曲线和预定应变幅值，确定所述血管壁的弹性值。

[0009] 在一些实现中，基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相包括：基于在拍摄所述目标对象的多维血管影像时同步记录的心电门控心电图信号确定所述血管壁的膨胀阶段。

[0010] 在一些实现中，基于目标对象的多维动态血管壁模型确定所述血管壁的膨胀阶段的多个时相包括：基于所述目标对象的多维动态血管壁模型确定在不同时相下所述血管壁的体积容量变化曲线；确定所述血管壁的体积容量变化曲线中体积单调递增的阶段作为所述膨胀阶段。

[0011] 在一些实现中，确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量包括：基于所述血管壁的膨胀阶段的初始时相的模型坐标和所述目标对象的多维血管影像的位移场信息确定所述膨胀阶段的每个时相的模型坐标；确定所述膨胀阶段的每个时相的模型坐标相对于所述膨胀阶段的初始时相的模型坐标的变形梯度；基于所述变形梯度确定所述血管壁的膨胀阶段的每个时相的血管壁应变张量。

[0012] 在一些实现中，所述血管壁应变张量包括Cauchy‑Green应变张量或Green‑Lagrange应变张量。

[0013] 在一些实现中，确定所述血管壁的膨胀阶段的每个时相的血管壁应力张量包括：为血管壁模型选择材料模型，并设置所述血管壁模型的材料参数以使得模拟出的血管壁接近刚性壁面；基于所述血管壁模型确定所述血管壁的能量方程；将所述血管壁的膨胀阶段的每个时相的血压数据施加到所述血管壁模型中，并且利用反向有限元分析对所述血管壁的能量方程进行求解以确定所述血管壁应力张量。

[0014] 在一些实现中，所述材料模型包括超弹性材料模型，并且其中所述血管壁应力张量为Cauchy应力或Second PK应力。

[0015] 在一些实现中，对所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量进行拟合以获取所述血管壁的拟合曲线包括：基于所述血管壁的膨胀阶段的多个时相的血管壁应变张量和血管壁应力张量确定用于拟合的血管壁等效应变和血管壁等效应力；为所述血管壁选择符合生物超弹性材料特性的拟合公式；以及基于所述血管壁等效应变和所述血管壁等效应力对所述拟合公式进行求解以确定拟合系数，从而确定所述血管壁的拟合曲线。

[0016] 在一些实现中，基于所述血管壁的拟合曲线和预定应变幅值，确定所述血管壁的弹性值包括：基于所述血管壁的拟合曲线求取所述血管壁等效应力相对于所述血管壁等效应变的一阶偏导数；以及将所述一阶偏导数与所述预定应变幅值之商确定为所述血管壁的弹性值。

[0017] 在一些实现中，该方法还包括：基于所述血管壁的拟合曲线求取所述血管壁等效应力相对于所述血管壁等效应变的二阶偏导数；将所述一阶偏导数与所述预定应变幅值之商以及将所述二阶偏导数与所述预定应变幅值之商确定为所述血管壁的弹性值。

[0018] 在一些实现中，该方法还包括：基于所述目标对象的多维血管影像确定所述目标对象的多维动态血管壁模型。

[0019] 在一些实现中，该方法还包括：基于所述目标对象的多维血管影像确定所述目标对象的位移场信息。

[0020] 根据本发明的另一个方面，提供了一种计算设备。该计算设备包括：至少一个处理器；以及至少一个存储器，该至少一个存储器被耦合到该至少一个处理器并且存储用于由该至少一个处理器执行的指令，该指令当由该至少一个处理器执行时，使得该计算设备执行根据上述方法的步骤。

[0021] 根据本发明的再一个方面，提供了一种计算机可读存储介质，其上存储有计算机程序代码，该计算机程序代码在被运行时执行如上所述的方法。

附图说明

[0022] 通过参考下列附图所给出的本发明的具体实施方式的描述，将更好地理解本发明，并且本发明的其他目的、细节、特点和优点将变得更加显而易见。

[0023] 图1示出了用于实现根据本发明的实施例的确定血管壁的弹性特征的方法的系统的示意图。

[0024] 图2示出了根据本发明的一些实施例的用于确定血管壁的弹性特征的方法的流程图。

[0025] 图3示出了根据本发明一些实施例的血管体积在一个心动周期的变化曲线的示意图。

[0026] 图4示出了根据本发明一些实施例的确定血管壁应变张量的过程的示意性流程。

[0027] 图5示出了根据本发明一些实施例的确定血管壁应力张量的过程的示意性流程。

[0028] 图6示出了根据本发明一些实施例的获取血管壁的拟合曲线的过程的示意性流程。

[0029] 图7示出了根据本发明的血管壁拟合曲线的示意图。

[0030] 图8示出了根据本发明一些实施例的用于确定目标对象的位移场信息以进行图像配准的方法的示例性流程图。

[0031] 图9示出了根据本发明一些实施例的用于确定目标对象的多维动态血管壁模型的方法的示意性流程图。

[0032] 图10示出了适合实现本发明的实施例的计算设备的结构方框图。

具体实施方式

[0033] 下面将参照附图更详细地描述本发明的优选实施方式。虽然附图中显示了本发明的优选实施方式，然而应该理解，可以以各种形式实现本发明而不应被这里阐述的实施方式所限制。相反，提供这些实施方式是为了使本发明更加透彻和完整，并且能够将本发明的范围完整的传达给本领域的技术人员。

[0034] 在下文的描述中，出于说明各种发明的实施例的目的阐述了某些具体细节以提供对各种发明实施例的透彻理解。但是，相关领域技术人员将认识到可在无这些具体细节中的一个或多个细节的情况来实践实施例。在其它情形下，与本申请相关联的熟知的装置、结构和技术可能并未详细地示出或描述从而避免不必要地混淆实施例的描述。

[0035] 除非语境有其它需要，在整个说明书和权利要求中，词语“包括”和其变型，诸如“包含”和“具有”应被理解为开放的、包含的含义，即应解释为“包括，但不限于”。

[0036] 在整个说明书中对“一个实施例”或“一些实施例”的提及表示结合实施例所描述的特定特点、结构或特征包括于至少一个实施例中。因此，在整个说明书的各个位置“在一个实施例中”或“在一些实施例”中的出现不一定全都指相同实施例。另外，特定特点、结构或特征可在一个或多个实施例中以任何方式组合。

[0037] 此外，说明书和权利要求中所用的第一、第二、第三、第四等术语，仅仅出于描述清楚起见来区分各个对象，而并不限定其所描述的对象的大小或其他顺序等。

[0038] 图1示出了用于实现根据本发明的实施例的确定血管壁的弹性特征的方法的系统1的示意图。如图1中所示，系统1可以包括操作台10、扫描床20和射线发生器30，其例如可以是一个CT系统。在系统1工作时，患者可以躺在扫描床20上，医生或操作员可以通过操作台
10控制扫描床20移动，以使得射线发生器30发出的射线对患者特定部位进行扫描，并将扫描产生的多维血管影像返回给操作台10。在本文中，多维血管影像不局限于CT影像，还可以包括MRI（核磁共振成像，Nuclear Magnetic Resonance Imaging）影像、血管超声成像影像等，取决于多维血管影像的类型不同，系统1可以有不同的结构和形式，而不局限于图1中所示的具体结构和形式。

[0039] 在操作台10处，或者，在独立于操作台10的另一计算设备（如医生的计算设备，图中未示出）处，可以对上述产生的多维血管影像进行处理和分析以获取所需要的结果。在这种情况下，操作台10或另一计算设备（本文中也统称为计算设备）可以包括至少一个处理器和与该至少一个处理器耦合的至少一个存储器，该存储器中存储有可由该至少一个处理器执行的指令，该指令在被该至少一个处理器执行时执行如下所述的方法的至少一部分。计算设备的具体结构例如可以如下结合图10所述。

[0040] 图2示出了根据本发明的一些实施例的用于确定血管壁的弹性特征的方法100的流程图。方法100例如可以由图1中所示的系统1中的操作台10或另一计算设备执行。以下以在操作台10中执行为例，结合图1至图10对方法100进行描述。

[0041] 如图2中所示，方法100可以包括方框110，其中可以基于目标对象的多维动态血管壁模型确定血管壁的膨胀阶段的多个时相。

[0042] 在一些实施例中，可以基于在拍摄目标对象的多维血管影像时同步记录的心电门控心电图（ECG）信号来确定血管壁的膨胀阶段。心电门控是指在采集心脏周期性节律性运动带来的影像数据时，开启心动周期上的一段“时间窗”，使得采集数据与周期性节律性的心电活动同步，从而相对制动心脏运动的磁共振生理同步采集技术。在拍摄目标对象的多维血管影像时，通常可以同步记录患者的心电图，在使用心电门控方式记录心电图时，可以根据心电图曲线容易地确定心脏舒张期血管壁处于膨胀阶段。

[0043] 在另一些实施例中，可以基于目标对象的多维动态血管壁模型确定在不同时相下血管壁的体积容量变化曲线，并且确定血管壁的体积容量变化曲线中体积单调递增的阶段作为血管壁的膨胀阶段，从而确定相应的时相。

[0044] 图3示出了根据本发明一些实施例的血管体积在一个心动周期的变化曲线的示意图。如图3中所示，假设血管体积在时相2开始单调递增，并且在时相9达到了峰值，此后开始递减，则时相2至时相9可以被确定为该目标对象的血管壁的膨胀阶段的时相。

[0045] 这里，该目标对象的多维动态血管壁模型可以基于目标对象的多维血管影像确定。例如，目标对象可以是指患者或者患者的特定部位，如头部、胸部、腹部等。基于目标对象的多维血管影像确定该目标对象的多维动态血管壁模型的方法可以是各种已知的或未来开发的方法，更具体地，其例如可以是如下结合图9所述的方法，但是本领域技术人员可以理解，本发明并不局限于图9所述的确定多维动态血管壁模型的方法。并且，可以在方法100开始之前独立地或者作为方法100的一部分来确定该目标对象的多维动态血管壁模型。

[0046] 此外，在一些实施例中，在方法100开始之前独立地或者作为方法100的一部分，还可以基于目标对象的多维血管影像确定该目标对象的位移场信息。这里，确定目标对象的位移场信息以用于对目标对象的多维血管影像进行配准，其可以是各种已知的或未来开发的方法，更具体地，其例如可以是如下结合图8所述的方法，但是本领域技术人员可以理解，本发明并不局限于图8所述的确定目标对象的位移场信息以进行图像配准的方法。

[0047] 在方框120，可以基于目标对象的多维动态血管壁模型（例如如下所述通过图9所述的方法确定的）和目标对象的多维血管影像的位移场信息（例如如下所述通过图8所述的方法确定的），确定血管壁的膨胀阶段的每个时相的血管壁应变张量。这里，可以利用正向有限元分析方法来确定血管壁应变张量。例如，可以利用常规有限元技术或商业有限元软件如ANSYS、Abaqus等。

[0048] 图4示出了根据本发明一些实施例的确定血管壁应变张量的过程（方框120）的示意性流程。

[0049] 如图4中所示，在方框122，可以基于血管壁的膨胀阶段的初始时相（例如如上所述在方框110确定的）的模型坐标和目标对象的多维血管影像的位移场信息确定膨胀阶段的每个时相的模型坐标。

[0050] 例如，血管壁的膨胀阶段的每个时相的模型坐标x可以表示为：

[0051]            (1)

[0052] 其中，X表示血管壁的膨胀阶段的初始时相的模型坐标， 表示目标对象的多维血管影像的位移场信息，其表示的是模型坐标和时间（时相）之间的关系。

[0053] 在方框124，可以确定膨胀阶段的每个时相的模型坐标x相对于膨胀阶段的初始时相的模型坐标X的变形梯度。例如，该变形梯度F可以表示为：

[0054]                                     (2)

[0055] 在方框126，可以基于该变形梯度F确定血管壁的膨胀阶段的每个时相的血管壁应变张量。

[0056] 取决于所采用的应变张量的类型，血管壁应变张量可以有不同的计算方式。这里，T可以基于变形梯度F和变形梯度F的转置F来确定血管壁应变张量，该血管壁应变张量可以是Cauchy‑Green应变张量或Green‑Lagrange应变张量。

[0057] 例如，Cauchy‑Green（柯西‑格林）应变张量C可表示为：

[0058]        (3)

[0059] 其中， 表示应变张量C在三个维度的应变分量，当i=j表示该分量为正应变分量，当i≠j则表示该分量为切应变分量。

[0060] 另一方面，Green‑Lagrange（格林‑拉格朗日）应变张量E可以表示为：

[0061]      (4)

[0062] 其中，I是与变形梯度F大小相同的单位矩阵， 表示应变张量E在三个维度的应变分量，当i=j表示该分量为正应变分量，当i≠j则表示该分量为切应变分量。

[0063] 继续图2，在方框130，可以基于目标对象的多维动态血管壁模型和目标对象的血压参数，确定血管壁的膨胀阶段的每个时相的血管壁应力张量。这里，目标对象的血压参数例如可以在拍摄目标对象的多维血管影像时同步记录。这里，可以利用反向有限元分析方法来确定血管壁应力张量。

[0064] 图5示出了根据本发明一些实施例的确定血管壁应力张量的过程（方框130）的示意性流程。

[0065] 如图5中所示，在方框132，可以为目标对象的血管壁模型选择适当的材料模型，并且设置该血管壁模型的材料参数以使得模拟出的血管壁接近刚性壁面。例如，适合模拟血管壁的材料模型可以是超弹性材料模型，例如可以包括neo‑Hookean模型、Mooney‑Rivlin模型、Odgen模型等。

[0066] 在方框134，可以基于血管壁模型确定该血管壁的能量方程。例如，在选择的材料模型是专门针对腹主动脉瘤仿真提出的超弹性模型的情况下，其血管壁的变形能量方程可以表示为：

[0067]    (5)

[0068] 其中， ，即 的三角矩阵， ，表示变形梯度F的行列式， 是材料系数。

[0069] 在方框136，可以将血管壁的膨胀阶段的每个时相的血压数据（例如可以在获取目标对象的多维血管影像时同时获取）施加到上述血管壁模型，并且利用反向有限元分析对血管壁的能量方程（上述公式（5））进行求解以确定血管壁应力张量。

[0070] 所确定的血管壁应力张量例如可以为柯西（Cauchy）应力或Second PK应力。

[0071] 例如，在使用柯西应力的情况下，血管壁应力张量的准静态平衡方程的弱解形式可以表示为：

[0072]                       (6)

[0073] 其中， 是结构的当前状态， 是柯西应力分量， 血管壁的i方向的位移分量关于j方向的导数， 是血管壁的密度， 是血
管壁的体力， 是预定义的血压施加在血管壁上的面载荷，而 是血管壁流腔的内边界。
由当前形变状态反向推导初始零载荷状态时，需要将正向运动方程 转置成反向运动方程 ，同时将正向的形变梯度F逆向转换，得到 ，带
入弱公式计算得初始状态模型形态。

[0074] 根据公式（5）和公式（6），可以计算柯西应力σ为：

[0075]   (7)

[0076] 其中 是反向变形张量， ，表示形变梯度F的逆变换f的行列式。

[0077] 此外，在方框130中，可以进一步确定血管壁张力来代替血管壁应力张量。例如，在血管壁厚度已知的情况下，可以进一步确定血管壁张力 为：

[0078]               (8)

[0079] 其中，h表示血管壁厚度，σ表示血管壁应力张量。

[0080] 在方框120中获取了血管壁的膨胀阶段的多个时相的血管壁应变张量，以及在方框130中获取了对应的血管壁应力张量之后，在方框140，可以基于所获取的血管壁应变张量和血管壁应力张量来获取血管壁的拟合曲线。

[0081] 图6示出了根据本发明一些实施例的获取血管壁的拟合曲线的过程（方框140）的示意性流程。

[0082] 如图6中所示，在方框142，可以基于上述方框120确定的血管壁的膨胀阶段的多个时相的血管壁应变张量和上述方框130确定的血管壁应力张量确定用于拟合的血管壁等效应变和等效应力。

[0083] 如上所述，所确定的血管壁应变张量和血管壁应力张量存在六个分量，因此在进行拟合时可以选择适当的分量或其组合，例如最大主应力/应变、Von Mises应力/应变（即屈服准则下的当量应力/应变）、等效应力/应变等，作为血管壁的等效应变和等效应力。此外，所使用的应变张量和应力张量的类型也具有固定的对应关系。例如，在使用上述公式（3）所示的Cauchy‑Green应变时，对应的应力应当为Cauchy应力，而使用如上述公式（4）所示的Green‑Lagrange应变时，对应的应力应当为Second PK应力。

[0084] 在一种具体实例中，假设使用的应变类型为Green‑Lagrange应变，使用的应力类型为Second PK应力 ，并且使用等效应力/应变的形式，则用于拟合的血管壁等效应变和血管壁等效应力 可以如下计算：

[0085]        (9)

[0086]   (10)

[0087] 在方框144，可以为血管壁选择符合生物超弹性材料特性的拟合公式。这里，拟合公式的选择只需满足血管壁的材料模型（如上述方框132中所选择的超弹性材料模型）的力学特性即可，例如可以是二次曲线拟合、三次曲线拟合、幂函数曲线拟合等。

[0088] 在方框146，可以基于方框142确定的血管壁等效应变（如上述公式（9）所示的血管壁等效应变 ）和血管壁等效应力（如上述公式（10）所示的血管壁等效应力 ）对方框144选择的拟合公式进行求解以确定相应的拟合系数，从而确定血管壁的拟合曲线。图7示出了根据本发明的血管壁拟合曲线的示意图。如图7中所示，所确定的血管壁的拟合曲线指示了在不同时相下，该目标对象的血管壁等效应变和血管壁等效应力之间的关系。

[0089] 例如，假设选择如上述公式（9）所示的血管壁等效应变 和如上述公式（10）所示的血管壁等效应力 ，并且使用幂函数曲线拟合算法，则所构建的拟合公式可以表示为：

[0090]     (11)

[0091] 其中， 、 和 为拟合系数，可采用最小二乘或最小化损失函数等方法求解。

[0092] 继续图2，在方框150，可以基于方框140所确定的血管壁的拟合曲线和预定应变幅值，确定血管壁的弹性值。

[0093] 在本文中，设计了独特的血管壁弹性指数（SSI，Strain‑scaled stiffness index）作为血管壁的弹性值来指示血管壁的弹性特征。具体地，可以先针对上述公式（11）所示的血管壁拟合曲线求血管壁等效应力 相对于血管壁等效应变 的一阶偏导数：

[0094]

[0095] 然后，将该一阶偏导数 与预定应变幅值 之商确定为血管壁的弹性值，即

[0096]                  (12)

[0097] 其中， 和 是对上述公式（11）所示的拟合公式求解得到的拟合系数，N为在拟合曲线上进行均匀采样的采样点的数量。在本发明的一些实例中，通过数值试验表明，当N超过16后，继续增加采样点，SSI值的变化不会超过5%，在这种情况下，可以将N设置为16。

[0098] 在一些进一步的实施例中，除了SSI之外，血管壁的弹性特征还可以进一步通过SSI的偏导数来指示。在这种情况下，在方框160，还可以基于血管壁的拟合曲线求取血管壁等效应力 相对于血管壁等效应变 的二阶偏导数 。

[0099] 然后，将上述一阶偏导数 与预定应变幅值 之商以及将该二阶偏导数 与预定应变幅值 之商确定为血管壁的弹性值，即

[0100]                  (12)

[0101]         (13)

[0102] 可以看出，根据本发明的思想构造的血管壁弹性指数SSI和/或dSSI这两个参数主要由两个组成部分组成，第一个分量是收缩期应变幅度（即预定应变幅值 ）的倒数，第二个分量代表平均拉伸弹性刚度（或平均拉伸应变刚度的变化率）。由于血管壁（如主动脉壁）的应力‑应变曲线表现出高度非线性，弹性刚度及其变化率通常随着应变的增加而呈非线性上升趋势。单个时相下的值不能全面的反应血管壁的力学特性。从图7可以看出拟合曲线上的每个数据点对SSI和dSSI计算的贡献，因此使用平均拉伸弹性刚度计算了整个应变范围内的平均弹性刚度，更全面的反应了血管弹性特性，相应的其变化率也是如此。此外，当前已经有许多研究表明主动脉瘤破裂风险与应变幅值之间的关系，普遍认为主动脉瘤患者的应变幅值会明显降低，因此，应变幅值被引入SSI和dSSI指数中，以放大这些指数的变化，并更好地反映血管壁不同区域的弹性特性变化状态。

[0103] 利用本文提出的血管壁弹性指数SSI及其变化率dSSI来反应血管壁的硬化病变程度，在不需要假设无应力状态，也不需要假设材料本构模型的情况下，构建了一套基于多维动态影像的在体无创血管壁力学特征参数识别方法。

[0104] 图8示出了根据本发明一些实施例的用于确定目标对象的位移场信息以进行图像配准的方法800的示例性流程图。

[0105] 如图8中所示，方法800可以包括方框810，其中可以采集目标对象的多个时相的多维血管影像。这里，目标对象是指患者或者患者的特定部位，如头部、胸部、腹部等。多维血管影像例如可以通过CT血管造影、核磁共振成像或血管超声成像等医学成像技术，并结合心电门控（ECG）技术获取。

[0106] 在方框820，可以从方框810采集的多个时相的多维血管影像中确定一个参考时相的多维血管影像作为目标图像并且选择另一时相的多维血管影像作为浮动图像。

[0107] 如本领域所知的，图像配准是将不同时相拍摄的两个不同图像的对应点达到空间位置和解剖结构的一致。因此，可以以一个参考时相的多维血管影像（即目标图像）为基准，确定所有其他时相的多维血管影像（即浮动图像）相对于该目标图像的位移信息以构建所有多维血管影像的位移场信息。这里，参考时相可以是任一时相。例如，可以选择采集的起始时相作为参考时相。

[0108] 在方框830，可以确定方框820选择的目标图像和浮动图像之间的位移信息。例如，对于目标图像I和浮动图像J，该位移信息可以由目标图像I和浮动图像J之间的一个空间变换T来表示。利用该空间变换T能够将目标图像I变换到新的坐标系中，即配准的变换过程，，将空间变换T作用于一个图像（如目标‑1
图像I）可以得到一个新的图像，将空间变换T的逆变换T 作用于该新的图像则能够得到原来的图像。该变换T的参数包括时间t，空间坐标x，以及图像域上的速度场 ，其中速度场是平方可积的连续向量场，满足 。

[0109] 在本文中，空间变换T被构建为对称微分同胚变换，即将空间变换T拆分为第一变换T1和第二变换T2，并且使用第一变换T1从目标图像I到浮动图像J的路径与使用第二变换T2从浮动图像J到目标图像I的路径相等。

[0110] 对称微分同胚变换是将空间变换 分为两部分进行计算，分别为第一变换 和第二变换 ，两者分别作用于参考图像I和浮动图像J，并满足如下条件：

[0111]

[0112] 对称微分同胚变换保证了无论采用何种相似度矩阵和优化策略，从目标图像I到浮动图像J和从浮动图像J到目标图像I计算时的路径相等，从而使得浮动图像和目标图像在计算过程中的地位是对等的。

[0113] 在本文的实例中，可以基于目标图像I本身的信息、浮动图像J本身的信息以及目标图像I和浮动图像J之间的互相关信息来构建衡量相似度的指标并且寻找使得该相似度的指标达到最小的空间变换T。

[0114] 这里，可以使用迭代方式对第一变换T1和第二变换T2进行更新直至达到收敛。在此过程中，第一变换T1和第二变换T2被初始化为单位矩阵，如果将迭代次数设置得过低或者将收敛阈值设置得过大，则很可能在满足收敛条件时并不能得到最佳的空间变换，而如果将迭代次数设置得过高或者将收敛阈值设置得过小，则可能会降低运行效率。

[0115] 针对这种情况，为了提高迭代过程中的收敛效率，可以采用多尺度配准，即通过采样间隔依次降低的多次降采样来对目标图像I和浮动图像J进行降采样，并且将前一次降采样获得的第一变换T1和第二变换T2作为下一次降采样的第一变换T1和第二变换T2的初始值。

[0116] 在方框840，可以基于所获取的参考时相与多个时相中的所有其他时相的多维血管影像的位移信息确定用于对多个时相的多维血管影像进行配准的位移场。

[0117] 图9示出了根据本发明一些实施例的用于确定目标对象的多维动态血管壁模型的方法900的示意性流程图。

[0118] 如图9中所示，方法900可以包括方框910，其中可以采集目标对象的多个时相的多维血管影像。

[0119] 这里，例如可以采用心电门控技术，采集目标对象的多个时相的多维血管影像。心电门控是指在采集心脏周期性节律性运动带来的影像数据时，开启心动周期上的一段“时间窗”，使得采集数据与周期性节律性的心电活动同步，从而相对制动心脏运动的磁共振生理同步采集技术。这里，在每个心动周期采集的影像称为目标对象的一个时相的多维血管影像。

[0120] 在方框920，可以从方框910中所采集的多个时相中选择一个参考时相。

[0121] 这里，参考时相的选择用于产生目标对象的参考模型。在一些实施例中，可以选择图像质量最好的一个时相作为参考时相。具体地，可以分别确定方框910中所采集的多个时相的多维血管影像的图像质量，并且从中选择图像质量最高的多维血管影像所处的时相作为该参考时相。这里，图像质量是指图像的对比度、模糊、信噪比、伪像和失真中的一项或多项。本领域技术人员可以理解，本发明并不局限于此，也可以选择多个时相中的任一时相，如初始时相，作为参考时相。此外，这里的参考时相可以与上述结合图8所述的方法800中用于图像配准的参考时相相同或不同。

[0122] 在方框930，可以基于参考时相的多维血管影像生成该目标对象的参考模型。

[0123] 在一些实施例中，可以基于预定图像阈值范围对该参考时相的多维血管影像进行分割以形成初步分割图像。初步分割图像主要用于分割出感兴趣的组织区域（如血管区域）。

[0124] 在一些实施例中，可以先将该参考时相的多维血管影像转换为灰度图像，然后提取该灰度图像中处于该预定图像阈值范围的部分作为初步分割图像。

[0125] 然后，可以利用边缘检测算法提取该初步分割图像的边缘信息。在灰度图像中，器官边缘处会表现出灰度值局部突变的特征，因此可以通过边缘检测算法（如Canny边缘检测算法）来提取该初步分割图像的边缘信息，以用于进一步确定相应的组织。

[0126] 之后，可以基于所提取的边缘信息从该初步分割图像中分离不同组织的组织区域，并且提取感兴趣的组织对应的连通区域。

[0127] 在一些实施例中，可以利用形态学的膨胀算法来对边缘信息进行处理以更准确地确定感兴趣的组织对应的区域。具体地，膨胀算法可以把图像中物体周围的背景点合并到物体中，如果两个物体之间距离比较近，那么膨胀算法可以将这两个物体连接在一起，膨胀算法对于填补图像分割后物体中的空洞很有用。

[0128] 之后，可以将包含感兴趣的组织所处的连通区域的初步分割图像转换为网格形式，以产生目标对象的参考模型。

[0129] 在获得目标对象的参考模型之后，在方框940，可以基于对方框910采集的多维血管影像进行图像配准得到的位移场信息和上述参考模型来确定目标对象的多维动态模型。

[0130] 这里，目标对象的位移场信息例如可以通过上述图8所述的方法800或其他任何适当的方法获得。

[0131] 由于图像配准得到的位移信息处于体素格上，相对离散，而按照上述方法确定的目标对象的参考模型相对连续，难以正好在某个体素格上，因此需要通过适当的方式对二者进行叠加以得到目标对象的实际模型，即，多维动态模型。

[0132] 在一些实施例中，对于一个多维血管影像中的任一点（如点p），可以首先基于该点的空间坐标确定该点所处的体素。

[0133] 然后，可以基于图像配准的位移场信息确定点p周围的8个体素的位移信息。这里，以三线性插值为例，点p周围的8个体素包括点p所处体素周围的八个顶点。

[0134] 之后，可以基于点p周围的8个体素的位移信息，利用线性插值确定点p的位移以将点p转换为多维动态模型上的点p’（模型点）。

[0135] 在对方框910所采集的多个时相的多维血管影像中的每个点都执行上述转换以获得相应的模型点之后，对所有模型点进行叠加获得目标对象的三维动态模型。

[0136] 图10示出了适合实现本发明的实施例的计算设备1000的结构方框图。计算设备1000例如可以是如上所述的用于执行方法100的操作台10或另一计算设备。

[0137] 如图10中所示，计算设备1000可以包括一个或多个中央处理单元（CPU）1010（图中仅示意性地示出了一个），其可以根据存储在只读存储器（ROM）1020中的计算机程序指令或者从存储单元1080加载到随机访问存储器（RAM）1030中的计算机程序指令，来执行各种适当的动作和处理。在RAM 1030中，还可存储计算设备1000操作所需的各种程序和数据。CPU 1010、ROM 1020以及RAM 1030通过总线1040彼此相连。输入/输出（I/O）接口1050也连接至总线1040。

[0138] 计算设备1000中的多个部件连接至I/O接口1050，包括：输入单元1060，例如键盘、鼠标等；输出单元1070，例如各种类型的显示器、扬声器等；存储单元1080，例如磁盘、光盘等；以及通信单元1090，例如网卡、调制解调器、无线通信收发机等。通信单元1090允许计算设备1000通过诸如因特网的计算机网络和/或各种电信网络与其他设备交换信息/数据。

[0139] 上文所描述的方法100例如可由计算设备1000（如操作台10或另一计算设备）的CPU 1010执行。例如，在一些实施例中，方法100可被实现为计算机软件程序，其被有形地包括于机器可读介质，例如存储单元1080。在一些实施例中，计算机程序的部分或者全部可以经由ROM 1020和/或通信单元1090而被载入和/或安装到计算设备1000上。当计算机程序被加载到RAM 1030并由CPU 1010执行时，可以执行上文描述的方法100的一个或多个操作。此外，通信单元1090可以支持有线或无线通信功能。

[0140] 本领域技术人员可以理解，图10所示的计算设备1000仅是示意性的。在一些实施例中，计算设备1000可以包含更多或更少的部件。

[0141] 以上结合附图对根据本发明的用于确定血管壁的弹性特征的方法100以及可用作操作台10或另一计算设备的计算设备1000进行了描述。然而本领域技术人员可以理解，方法100的步骤及其子步骤的执行并不局限于图中所示和以上所述的顺序，而是可以以任何其他合理的顺序来执行。此外，计算设备1000也不必须包括图10中所示的所有组件，其可以仅仅包括执行本发明中所述的功能所必须的其中一些组件，并且这些组件的连接方式也不局限于图中所示的形式。

[0142] 本发明可以是方法、装置、系统和/或计算机程序产品。计算机程序产品可以包括计算机可读存储介质，其上载有用于执行本发明的各个方面的计算机可读程序指令。

[0143] 在一个或多个示例性设计中，可以用硬件、软件、固件或它们的任意组合来实现本发明所述的功能。例如，如果用软件来实现，则可以将所述功能作为一个或多个指令或代码存储在计算机可读介质上，或者作为计算机可读介质上的一个或多个指令或代码来传输。

[0144] 本文公开的装置的各个单元可以使用分立硬件组件来实现，也可以集成地实现在一个硬件组件，如处理器上。例如，可以用通用处理器、数字信号处理器（DSP）、专用集成电路（ASIC）、现场可编程门阵列（FPGA）或其它可编程逻辑器件、分立门或者晶体管逻辑、分立硬件组件或用于执行本文所述的功能的任意组合来实现或执行结合本发明所描述的各种示例性的逻辑块、模块和电路。

[0145] 本领域普通技术人员还应当理解，结合本发明的实施例描述的各种示例性的逻辑块、模块、电路和算法步骤可以实现成电子硬件、计算机软件或二者的组合。

[0146] 本发明的以上描述用于使本领域的任何普通技术人员能够实现或使用本发明。对于本领域普通技术人员来说，本发明的各种修改都是显而易见的，并且本文定义的一般性原理也可以在不脱离本发明的精神和保护范围的情况下应用于其它变形。因此，本发明并不限于本文所述的实例和设计，而是与本文公开的原理和新颖性特性的最广范围相一致。