CN109799471B - 一种磁共振波谱成像模拟方法及系统、存储介质、电子设备 - Google Patents

Abstract

本发明公开了一种磁共振波谱成像模拟方法，包括如下步骤：根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；将脉冲序列中的模块的左右时间边界及其内部关键点投影到同一时间轴上，形成若干个长短不一且物理控制量恒定的序列片，对设置循环优先级的循环脉冲序列解析形成所有序列片；逐次对所有序列片进行核磁共振动力学演化模拟，得到模拟数据；该发明可通过直观统一的编程环境完成自旋体系构建、脉冲序列构建和解析以及优化、实验任务仿真和数据分析，本发明提出的可编程波谱成像模拟系统应用范围广泛。

Description

一种磁共振波谱成像模拟方法及系统、存储介质、电子设备

技术领域

本发明属于磁共振波谱和成像领域，具体涉及一种磁共振波谱成像模拟方法及系统。

背景技术

磁共振技术作为一种无损的探测手段已广泛应用于波谱学领域和医用成像诊断领域。前者可用于物质分子成分测定、活体小动物代谢定量，后者则广泛用于临床组织结构成像。磁共振模拟或仿真是基于磁共振自旋动力学理论的一种数值化仿真手段，除可用于演示教学目的，还大量用于新型脉冲序列开发、各类实验参数优化、理论与实验验证等研发工作，并与核磁仪器本身可形成优势互补，因而具有重要的应用价值。

本发明的提出主要面向以下两个现状问题：

一、现有磁共振模拟工具主要面向具体应用领域而设计，缺乏通用性。比如磁共振成像方面主要以Bloch方程作为计算模型，而波谱学方面则主要以Liouville方程作为计算模型。Bloch方程尽管简单，却无法处理复杂的耦合自旋体系和非1/2自旋体系，因而成像模拟的应用范围受限，比如无法进行Na成像模拟、无法进行磁共振波谱成像(MRSI)模拟等等。

二、无论是波谱学还是磁共振成像应用仿真，其相应的脉冲序列设计均缺乏兼具统一、直观、灵活的编写环境。脉冲序列一般由射频脉冲、梯度脉冲和信号采集构成，用于实现自旋体系的有效操控，从而产生预期的探测信号。因此，脉冲序列是磁共振的灵魂技术。尽管脉冲序列构件仅有三类，但其组合千变万化，实现的功能也各不相同。目前已开发出的脉冲序列已超过1000种，新的脉冲序列也伴随着应用的发展而不断被开发出来。现有的脉冲序列设计编写工具灵活性、直观性不足，比如布鲁克核磁谱仪和瓦里安核磁谱仪均采用C语言作为用户自定义接口，序列代码冗长、可读性较差。这在一定程度上阻碍了新型脉冲序列的创新研究。

发明内容

本发明的目的针对磁共振波谱(MRS)和成像(MRI)的不断融合发展的应用趋势，提出一种可同时兼容于磁共振波谱和成像的仿真模拟方法，可通过直观统一的编程环境完成自旋体系构建、脉冲序列构建和解析以及优化、实验任务仿真和数据分析。

本发明提供一种磁共振波谱成像模拟方法，包括如下步骤：

对象生成，根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；

脉冲序列生成，若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；

脉冲序列解析，将脉冲序列中的模块的左右时间边界及其内部关键点投影到同一时间轴上，形成若干个长短不一且物理控制量恒定的序列片，对设置循环优先级的循环脉冲序列解析形成所有序列片；

模拟实验，按照时间轴先后顺序，逐次对所有序列片进行核磁共振动力学演化模拟，得到模拟数据。

优选地，在步骤脉冲序列生成中，所述模块包括单一模块或高级模块，所述单一模块包括射频脉冲模块、梯度模块、信号采集模块，对所述单一模块中的任意模块通过串行和/或并行形成高级模块，对所述单一模块或高级模块设置循环优先级。

优选地，在步骤脉冲序列解析中还包括，所述循环优先级相同的多个单一模块或高级模块选取循环次数最大的单一模块或高级模块作为该优先级的循环次数；利用检索对循环脉冲序列进行循环次数和循环优先级统计，生成脉冲序列的全局循环矩阵；所述全局循环矩阵包括矩阵列和矩阵行，所述矩阵列为脉冲序列的优先级等级总数；所述矩阵行为脉冲序列各优先级的循环次数的乘积，即脉冲序列需要执行的总次数。

优选地，在步骤模拟实验中还包括：信号采集，利用重复多点采集将每一自旋体系在序列片结束时的状态投影到预期信号采集状态获取自旋体系对应的复数信号值和信号采集的原始数据；所述信号采集的原始数据用于后期数据分析处理。

优选地，在步骤模拟实验后还包括步骤：

数据处理和分析，模拟数据存储于本地或在线处理，利用函数对任意自旋体系状态的演化过程进行追踪和分析。

优选地，在步骤脉冲序列生成之前还包括如下步骤：

S11、数据流图创建，根据定义的自旋体系演化所需的数据流节点、自旋体系的初始态和期望态节点，定义自旋体系所需的目标函数节点，同时建立数据流连接；

S12、脉冲参数初始化，根据提供的自旋体系初始态和期望态分别为数据流图中的对应节点赋值；

S13、获得目标函数，利用对数据流图进行传播计算获得求解的最终目标函数，所述目标函数用于得出优化脉冲的梯度值；

S14、梯度迭代优化，利用步骤S13中的梯度值和当前的脉冲数据，对脉冲的各个节点数据进行更新；

S15、脉冲数据优化，重复步骤S13、S14直至目标函数取得极值，得到优化脉冲数据。

本发明还提供一种电子设备，包括：处理器，

存储器；以及程序，其中所述程序被存储在所述存储器中，并且被配置成由处理器执行，所述程序包括用于执行一种磁共振波谱成像模拟方法。

本发明还提供一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行一种磁共振波谱成像模拟方法。

本发明还提供一种磁共振波谱成像模拟系统，包括对象生成模块、脉冲序列生成模块、脉冲序列解析模块、模拟实验模块、数据流图创建模块、脉冲参数初始化模块、获得目标函数模块、梯度迭代优化模块、数据处理和分析模块、脉冲数据优化模块和信号采集模块；其中，

对象生成模块用于根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；

脉冲序列生成模块用于若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；

脉冲序列解析模块用于将脉冲序列中的模块的左右时间边界及其内部关键点投影到同一时间轴上，形成若干个长短不一且物理控制量恒定的序列片，对设置循环优先级的循环脉冲序列解析形成所有序列片；

模拟实验模块用于按照时间轴先后顺序，逐次对循环脉冲序列解析形成的所有序列片进行核磁共振动力学演化模拟，得到模拟数据；

数据流图创建模块用于根据定义的自旋体系演化所需的数据流节点、自旋体系的初始态和期望态节点，定义自旋体系所需的目标函数节点，同时建立数据流连接；

脉冲参数初始化模块用于根据提供的自旋体系初始态和期望态分别为数据流图中的对应节点赋值；

获得目标函数模块用于利用对数据流图进行传播计算从而获得求解的最终目标函数对于优化脉冲的梯度值；

梯度迭代优化模块用于利用获得目标函数模块中的梯度值和当前的脉冲数据，对脉冲的各个节点数据进行更新；

数据处理和分析模块用于模拟数据存储于本地或在线处理，利用函数对任意自旋体系状态的演化过程进行追踪和分析；

脉冲数据优化模块用于重复执行梯度迭代优化模块及数据处理和分析模块直至目标函数取得极值，得到的脉冲数据为优化脉冲数据；

信号采集模块用于利用重复多点采集将每一自旋体系在序列片结束时的状态投影到预期信号采集状态获取自旋体系对应的复数信号值和信号采集的原始数据；所述信号采集的原始数据用于后期数据分析处理。

优选地，所述脉冲序列解析模块包括全局循环矩阵单元，所述全局循环矩阵单元包括矩阵列和矩阵行，所述矩阵列为脉冲序列的优先级等级总数；所述矩阵行为脉冲序列各优先级的循环次数的乘积同时也为需要执行的总次数；所述脉冲序列生成模块包括单一模块或高级模块，所述单一模块包括射频脉冲模块、梯度模块、信号采集模块，对所述单一模块中的任意模块通过串行和/或并行形成高级模块，对所述单一模块或高级模块设置循环优先级。

相比与现有技术，本发明的有益效果在于：

本发明提出一种可同时兼容于磁共振波谱和成像的仿真模拟方法，可通过直观统一的编程环境完成自旋体系构建、脉冲序列构建和解析以及优化、实验任务仿真和数据分析，本发明提出的可编程波谱成像模拟系统应用范围广泛，包括新型脉冲序列研发、波谱成像模拟实验、基础教学演示等，本发明提出的脉冲序列编程语法还可移植应用于磁共振谱仪控制台，从而服务于核磁共振仪器系统。

本发明提供的一种磁共振动力学模拟方法和系统，以全面覆盖磁共振波谱学和磁共振成像的各类应用模拟需求，所有实验任务可编写、存储在脚本文件中，也可在模拟系统终端逐行执行。

本发明提供的一种磁共振波谱成像模拟系统，可工作于离线模式，即在计算机或工作站上运行，也可以部署于终端电子设备(如手机或平板)，实现在线的可编程设计，相关的模拟计算量可通过远程云计算完成，计算完成后传输模拟数据结果即可。

附图说明

图1为本发明的一种磁共振波谱成像模拟方法流程示意图；

图2为本发明的一种磁共振波谱成像模拟方法具体流程示意图；

图3为梯度回波序列的模块化分解示意图；

图4为脉冲序列分解示意图；

图5为本发明述及的一种最优控制的计算网络图模型；

图6(a)为本发明模拟一中的丙烯酸体系的模拟幅值谱图；

图6(b)为本发明模拟一中的500MHz布鲁克波谱仪采集到的试验图谱；

图7为运用本发明模拟二中的二维空间选择性MRI序列；

图8为运用本发明模拟二中的空间选择性脑部成像模拟图像；

图9为运用本发明对同一模体进行传统的自旋回波序列(SE)成像模拟获得的T1加权图像；

图10为运用本发明模拟三中的七步脉冲同步优化的脉冲形状结果图；

图11为运用本发明模拟三中的七步脉冲进行时频域的波谱分析结果图；

图12为脉冲序列构建步骤中的二维相关谱(COSY)试验过程示意图；

图13为本发明的脉冲序列构建的流程示意图；

图14为本发明的一种磁共振波谱成像模拟系统的示意图；

附图标记：10、射频脉冲模块，20、延时模块，30、梯度模块，40、并行模块(梯度模块+采集模块)，50、脉冲序列时间轴，60、物理控制量，70、成千上万的目标自旋体系，80、90度射频脉冲(rf90),90、180度反转脉冲(rf180)。

具体实施方式

下面通过实施例对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。

应当理解，本文所使用的诸如“具有”，“包含”以及“包括”术语并不排除一个或多个其它元件或其组合的存在或添加。

本发明提供一种可编程的磁共振动力学模拟方法和系统，首先基于量子自旋动力学理论构造目标自旋体系或成像模体模型，利用脉冲序列编程语法构建仿真实验序列，在此基础上，通过计算机运行仿真脚本进行数值计算完成实验任务，从而获取原始实验数据，继而进行数据分析如一维、二维、多维傅里叶变换(FFT)，得到波谱图或成像图；本发明示例的优选可编程环境为Lua脚本语言，本发明适应的编程语言不限于此。

一种磁共振波谱成像模拟方法，如图1-2所示，包括如下步骤：

S1、对象生成，根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；

本实施例中，自旋体系构造由spin_system{}函数完成，具体地，用户可根据实验样本的自旋体系参数提供：a、自旋体系组分，如spin＝"13C13C13C"指定该体系为三个碳13原子体系；b、化学位移参数，如zeeman＝"1scalar 15.74kHz"指定第一个碳13自旋的化学位移为15.74kHz，其他自旋化学位移可依此特点自行增加；c、耦合常数，如Jcoupling＝"12scalar54.2Hz2 3scalar 35.4Hz"指定自旋1和2之间的J耦合常数为54.2Hz，而自旋2与3之间的耦合常数则为35.4Hz；d、弛豫参数，如relaxation＝"T1 500ms T2 50ms"指定自旋体系的宏观T1和T2弛豫时间。

例如碳13标记的丙氨酸体系构建示例为：

local sys＝spin_system{

spin＝"13C 13C13C",

zeeman＝"1scalar 15.74kHz3scalar-4.3kHz",

Jcoupling＝"1 2scalar 54.2Hz2 3scalar 35.4Hz"}

基于上述自旋体系参数信息，仿真系统采用磁共振量子力学模型构建自旋体系模型，用户可利用一系列脚本函数获取必要的体系参量。

1、利用state{}获得自旋体系状态rho，例如sys:state("I1z+I2z+I3z")定义了体系的热平衡态I1z+I2z+I3z,注意自旋体系状态定义是由表达式完成的，这是最为接近用户的一种定义方式，可为使用者提供一种所见即所得的编程体验。自旋体系状态定义根据需要可经由数学表达式任意构造，比如组合态sys:state("2*I1xI3x-0.5*I2xI3y")。

2、利用sys.L0,sys.L和sys.R可分别提取自旋体系的自由演化哈密顿量，系统总体哈密顿量和体系弛豫矩阵。

成像模体则由成千上万个自旋体系构成，优选地，在磁共振成像领域，成像模体中自旋体系模型通常由单自旋氢核构成，成像模体构造需提供一组数据集，可包括但不限于模体物理尺寸、模体分辨率、模体空间偏移量、模体T1分布矩阵、T2分布矩阵、组织归属矩阵，上述参数矩阵大小可任意，但数据集内部矩阵大小相同，所述数据集可打包为单一MAT或HDF格式的文件。

S2、脉冲序列生成，若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；所述模块包括单一模块或高级模块，所述单一模块包括射频脉冲模块、梯度模块、信号采集模块，对所述单一模块中的任意模块通过串行和/或并行形成高级模块，对所述单一模块或高级模块设置循环优先级；本实施例中，具体地，如图3所示，脉冲序列的可编程设计是通过序列的模块化实现的，具体地，所有的脉冲序列均分解为两个等级的序列子块，等级一均为单一的不可再分模块，包括射频脉冲模块10、梯度模块30、延时模块20和信号采集模块等。等级二为高级模块，是由各模块串行(按照先后顺序执行)或者并行(同时同步执行)组合而成，如并行模块(梯度模块+采集模块)40；由此，所有脉冲序列最终均可由等级二的串行模块表示。

具体地，所述脉冲序列的构建包括如下步骤，如图13所示：

一、单一模块配置，所述单一模块配置包括射频脉冲模块配置、梯度模块配置和信号采集模块配置；

1、脉冲模块配置包括硬脉冲配置和形状脉冲配置；所述硬脉冲配置由hardRF{}函数完成,hardRF{beta＝60,phase＝"-y"}指定一个作用于氢核的60度翻转角，-y相位的硬脉冲,位参数不仅限于四个坐标轴(x/y/-x/-y)，也可以是任意的角度,于相位循环实验，可直接将相位循环列表(如"x-xy-y")添加到phase参数中，则模拟实验中该脉冲将逐次切换相应的相位参数；默认状态下，激发脉冲是作用在1H上的，作用在氢原子之外其他原子核时候需要定义channel，如hardRF{beta＝180,channel＝"1H|13C",phase＝"x|x"}定义为相位为x的作用在1H-13C两个不同的核自旋上的180度激发脉冲。

所述脉冲(Sinc,Gaussian,Rectangle等)配置由shapedRF{}函数完成，shapedRF{width＝5.12,step＝256,max_amp＝1000,pattern＝"sinc"}指定一个带宽为5.12ms，步长为256，最大幅度为1000的默认为5瓣的sinc脉冲；形状脉冲的pattern参数还支持外部文件格式和数学表达式；如shapedRF{channel＝"1H|13C",width＝5,pattern＝"shape.RF"}指定形状脉冲从shape.RF的.RF后缀文件中读取，带宽为5ms，作用在1H-13C两个不同的核自旋上。

2、所述梯度模块配置，脉冲场梯度最常用于MR脉冲序列中，或者用于成像中的空间编码，或者用于光谱学中特定磁化传递路径的重新相位(选择)和去相位(消除)。优选地，本发明示例提供trapGrad、shapedGrad和exprGrad，分别从典型的梯形形状、自定义波形文件和数学表达式创建梯度。

3、所述采集模块配置，数据采集模块acq用于信号观察，模拟磁共振扫描仪中接收器的正交检测。如acq{np＝1024,sw＝50000}指定采样点数为1024，采样率为50kHz。对于更复杂的采集需求，如观测自旋体系在特定感兴趣状态上的投影，可额外指定观测状态。如acq{np＝1024,sw＝50000,observer＝"2*I1zI2z"}指定采集的信号来源于自旋体系状态到预期的观测状态2I1zI2z的投影。

二、脉冲序列的构建，脉冲序列的构建装是指根据实验目的和实验任务要求进行组装的，将上述预先定义完成的单一模块进行串行(时间上先后执行)和并行(时间上同步执行)组合，从而进一步形成用于该实验的脉冲序列；优选地，本发明运用特定函数seq{}完成脉冲序列组装；序列子块之间用特定符号(逗号或分号)分隔，即为串行模块的组装；并行模块的组装可以通过特定运算符(如+)连接。

例如local s1＝seq{d1,rf,acq}表示一个简单的串行组装的脉冲序列(其中包含延迟块d1,射频块rf，采集块acq)，用于实现FID信号采集；再例如local s2＝seq{rf+gz,gxPre+gyPre+gzReph,gx+acq,delay}中由若干并行模块与单一模块进行混合式串行组装；由于脉冲序列本质上被抽象为等级二的串行模块，因此，脉冲序列本身也可以嵌套在任何其他序列中；此特性允许重复使用预先设计的脉冲序列模板，这可以进一步提高脚本编写效率。

三、循环脉冲序列构建，上述步骤二中的组装方式适应于单次实验，seq内部各串行模块运行完毕，即完成脉冲序列的模拟；而核磁共振波谱学如二维谱或核磁共振成像均依需通过阵列实验获得一系列编码信号，从而能够进行二维波谱或图像重建，因此，本发明提出一种简洁直观的用于阵列实验设计的脉冲序列语法，还包括全局循环和局部循环：

1、全局循环,全局阵列实验要求脉冲序列执行多次，每次执行需对脉冲序列中某一个或若干模块进行参数调整，本发明中该循环操作可通过在感兴趣模块后增加特定循环符号(如#)进行表征。

例如seq{rf90I,tau#,rf180,tau#,rf90,acq}为一种改进的INEPT阵列实验，其中tau为延迟模块，tau定义了N个不同的延时时间，该序列即设计了不同的延时时间下的阵列化实验，该序列执行完毕后，将获得N组原始信号。

对于三维或更高维实验任务，本发明通过在符号#后面分配不同的循环优先级数值，从而实现多维序列，数值越高，表示该层循环优先级越高。

如图12为相位循环功能的二维波谱COSY序列，其中，黑色矩形代表90度射频脉冲(rf90)80，白色矩形代表180度反转脉冲(rf180)90，不同延时t1条件下，该序列均需要按照相位循环表1对子模块rf90，rf180，acq同步进行相位循环处理，本例中循环次数为4次，因此，根据上述描述，该脉冲序列的延迟模块t1优先级需低于其他相位循环模块。

该脉冲序列可进一步组装为seq{rf90#2,t1#1,rf180#2,acq#2}，注意，相位循环周期中的步骤总数n(n＝4)由所有相关块的最大循环次数(本例中为rf180)决定，因此，如果延迟t1的循环次数为N，则该序列累计运行次数为N*n，所获得的原始信号数组数量也为N*n。

表1

此表格为表1自旋回波实验的四步相位周期，该脉冲序列执行的解释性过程如下：

seq{rf90(0),t1(0),rf180(0),acq(0)}

seq{rf90(1),t1(0),rf180(1),acq(1)}

seq{rf90(2),t1(0),rf180(2),acq(2)}

seq{rf90(3),t1(0),rf180(3),acq(3)}

seq{rf90(0),t1(1),rf180(0),acq(0)}

seq{rf90(1),t1(1),rf180(1),acq(1)}

seq{rf90(2),t1(1),rf180(2),acq(2)}

seq{rf90(3),t1(1),rf180(3),acq(3)}

......

2、局部循环，在脉冲序列全局循环过程中，在循环内往往还可能存在重复执行的子模块，为简化串行操作，本发明通过对目标模块追加特定符号(如“～”)来实现全局循环内部的重复性操作。

例如脉冲序列seq{rf,(gx+gy+acq)～}，假定gx和gy模块循环次数可定义为3次，则其与seq{rf,gx(0)+gy(0)+acq,gx(1)+gy(1)+acq,gx(2)+gy(2)+acq}等价，模块重复次数由该模块中的单一模块定义获得。

S3、脉冲序列解析，如图4所示，其中，成千上万的目标自旋体系70，将子序列的左右时间边界及其内部关键点投影到同一脉冲序列时间轴50上形成若干个长短不一且物理控制量60恒定的序列片；通过检索对含有循环标识的子序列进行循环次数和循环优先级统计，生成脉冲序列的全局循环矩阵；所述全局循环矩阵包括矩阵列和矩阵行，所述矩阵列为脉冲序列的优先级等级总数；所述矩阵行为脉冲序列各优先级的循环次数的乘积，即脉冲序列需要执行的总次数；具体地，脉冲序列编写完毕，经解析完成后可进一步执行模拟实验任务,解析过程如下：

例如，脉冲序列seq{A#1,B#2,C#3,D#2,E}，对序列内各子序列模块(如A，B，C，D)进行检索，对含有全局循环标识(#)的各子模块进行循环次数和循环优先级统计(例如A：定义循环次数3次，优先级1级；B循环次数4次，优先级2级；C循环次数2次，优先级3级；D循环次数5次，优先级2级)，从而获得优先级等级数量以及各优先级对应的循环次数，其中优先级相同的子模块选取最大的模块循环次数作为该优先级的循环次数。

根据上述循环次数和循环优先级的统计，生成全局循环矩阵C，其中矩阵列代表不同的优先等级，优选地，矩阵从左到右可按优先级升序或降序排列，矩阵列数即为该序列的优先级等级总数，比如本例中共有三级循环，列数即为3；矩阵行代表不同的全局循环索引，其行数为该序列累计需要执行的总次数，为各优先级的循环次数的乘积(如本例中优先级#1循环次数为3，#2循环次数为5，#3为2，即该序列累计执行次数为3*5*2)。

循环矩阵C的元素初始化的原则为：

1、每列的数值为自然数，用于可变序列子模块的索引定位；

2、每列按列方向逐行进行循环式赋值，从0开始，依次增加1，直至该列对应的优先级的循环次数；

3、相邻优先级的两列，优先级低的列在优先级高的列实施一次完整的循环赋值后方可进行增量(加1)赋值。

例：一种升序的循环矩阵示意(其中#3循环3次，#2循环2次，#1循环2次)，

#1	#2	#3
			0	0	0
0	0	1
			0	0	2
0	1	0
			0	1	1
0	1	2
			1	0	0
1	0	1
			1	0	2
1	1	0
			1	1	1
1	1	2

按照上述全局循环矩阵C，逐行切换到索引的脉冲序列模块参数后，将脉冲序列各子块左右时间边界及其内部关键时间点投影到同一时间轴，从而形成若干长短不一的序列片，各序列片的特点在于该序列片内的物理控制量(Ctrl Vector)均为恒定量，包括射频脉冲的幅值、相位和频率，梯度脉冲的轴向(X/Y/Z)、幅值、信号采集的标志位等；当全局循环矩阵所有行处理完毕后，即完成该脉冲序列的格式解析。

S4、模拟实验，按照时间轴先后顺序，逐次对循环脉冲序列解析形成的所有序列片进行核磁共振动力学演化模拟，得到模拟数据；具体地，在完成脉冲序列解析后，可进行核磁共振物理学过程模拟，以获取模拟信号数据及进行深入的后续处理等；

首先，按照时间轴先后顺序，逐次对所有序列片进行核磁共振动力学演化模拟，每序列片内外部控制量恒定，即保证自旋体系的总体哈密顿量不变，基于刘维尔(liouville)方程即可实现自旋体系从当前状态rho0(序列片开始)演化至下一状态rho1(序列片结束)；每一序列片的结束时状态均作为下一序列片的起始状态。对于成像模拟来说，自旋体系数量根据模体尺寸精度可达上万甚至上千万，因此，每一序列片的演化，所有自旋体系均需遵循上述过程。

自旋体系的动力学方程，其中rho为自旋体系状态，L为自旋体系的总体哈密顿量，由自旋体系的自由演化哈密顿量和外部射频场哈密顿量构成，L由该序列片内的射频场、梯度场和自旋体系本身特性决定；R为自旋体系的弛豫算符。

每一序列片的动力学计算可通过计算机并行计算实现加速，包括多核CPU加速和多核GPU加速。

其次，在序列片控制量出现信号采集标志时，即进行数据采集操作，具体地，将每一自旋体系在该序列片结束时的状态rho投影到预期信号采集状态rhoAcq，从而获取该自旋体系对应的复数信号值，如模拟对象为成像模体，则需将所有自旋体系的复数信号值进行累加，从而获得总体的复数信号值。

最后，根据实验需求和脉冲序列设计，信号采集一般可多次多点进行采集，进而获得一定大小的原始信号矩阵，可进一步用于后期数据分析处理(如二维谱或二维成像)；信号采集的原始数据，可重构为多维数据格式，从而对应于多维核磁共振实验。

在核磁共振中，自旋体系的状态的准确性影响着实验质量，因此，该领域所谓的自旋体系的最优控制，其主要目标是在特定的约束条件下，利用优化的脉冲形状，使得体系从初态迁移到末态后，该实际末态最大限度地保持与期望状态的一致性。

优化算法主要有数值优化方法和解析方法两类，脉冲优化所涉及的算法大多基于梯度方法，如GRAPE和Krotov，为了快速开发不同类型的控制策略，我们提出了一种基于计算图的优化控制框架，该框架可基于google的TensorFlow进行构造，并通过图的自动微分能力进行反向梯度求解，从而实现自动化的梯度优化计算。该方法无需人工求解目标函数的对应梯度表达式，因而极其适合于复杂多变的脉冲优化控制目的。

为了单一模块的形状射频脉冲提供优化波形，具体优化脉冲数据的方法包括如下步骤：

S11、数据流创建，定义自旋体系演化所需的数据流节点、自旋体系的初始态和期望态节点，进而定义自旋体系所需的目标函数节点，同时建立数据流连接；

主要包括两部分构建：一是自旋体系演化所需的数据流节点(包括张量数据和若干操作函数)定义；二是依据自旋动力学方程建立数据流连接；

具体地，如图5所示，将待优化脉冲离散化为N等步长段，由此形成N个数据输入节点，每一节点数据即为该步对应的射频脉冲波形参数。自旋体系的每一种核自旋的射频通道数据包括幅值A、相位phi和偏移频率df，因此，根据所定义的自旋体系的通道数量，上述各输入节点的数据格式uk即为：通道1(A，phi，df)、通道2(A，phi，df)....

依据自旋动力学理论,对于每一步k(k＝1,2,3…N)，分别建立自旋体系的当前态ρ_k-1和下一态ρ_k、自由演化哈密顿量H0、弛豫超算符R、射频场哈密顿量Hrf、自旋体系单步演化指数函数等关键节点，依据自旋动力学方程，建立上述所有节点的数据连接。

依据射频脉冲的优化目标任务，定义自旋体系的初态和期望态节点，在此基础上进一步定义所需的目标函数节点，同时建立数据流连接，即完成了自旋体系的演化数据流图创建。

S12、脉冲优化参数初始化,根据提供的自旋体系初始态和期望态分别为数据流图中的对应节点赋值；

根据提供的初始脉冲波形数据，为步骤S11中所述的N个射频脉冲节点进行初始化赋值，即(A1，phi1，df1)(A2，phi2，df2)…根据提供的自旋体系初始状态和期望状态分别为数据流图中的对应节点赋值。

S13、获得目标函数,利用对数据流图进行传播计算从而获得求解的最终目标函数对于优化脉冲的梯度值；

优选地，基于TensorFlow对上述初始化后的数据流图进行前向传播计算，可获得相应的目标函数值，同时，对其进行反向传播计算，可获得各节点对于待优化的脉冲数据的梯度值，从而自动获得求解的最终目标函数对于N步优化脉冲的梯度值；

S14、梯度迭代优化，利用步骤S13中的梯度值和当前的脉冲数据，对脉冲的各个节点数据进行更新；

基于上述N步梯度求解数值和当前的N步脉冲波形数据，进行脉冲波形更新，从而对N步脉冲的各个节点数据进行更新。

S15、脉冲数据优化,重复步骤S13、S14直至目标函数取得极值，此时得到的脉冲数据为优化脉冲数据；

重复实施S13、S14，直至目标函数取得极值，此时最新的脉冲波形数据即为优化脉冲数据；可用于进一步集成到自定义脉冲序列，开展仿真模拟实验，或直接应用于实际磁共振仪器实验中。

一种电子设备，包括：处理器，存储器；以及程序，其中所述程序被存储在所述存储器中，并且被配置成由处理器执行，所述程序包括用于执行一种磁共振波谱成像模拟方法。一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行一种磁共振波谱成像模拟方法。

一种磁共振波谱成像模拟系统，如图14所示，包括对象生成模块、脉冲序列生成模块、脉冲序列解析模块、模拟实验模块、数据流图创建模块、脉冲参数初始化模块、获得目标函数模块、梯度迭代优化模块、数据处理和分析模块、脉冲数据优化模块和信号采集模块；其中，

对象生成模块，根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；

脉冲序列生成模块用于若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；

脉冲序列解析模块用于将脉冲序列中的模块的左右时间边界及其内部关键点投影到同一时间轴上，形成若干个长短不一且物理控制量恒定的序列片，对设置循环优先级的循环脉冲序列解析形成所有序列片；

模拟实验模块用于按照时间轴先后顺序，逐次对循环脉冲序列解析形成的所有序列片进行核磁共振动力学演化模拟，得到模拟数据；

数据流图创建模块用于根据定义的自旋体系演化所需的数据流节点、自旋体系的初始态和期望态节点，定义自旋体系所需的目标函数节点，同时建立数据流连接；

脉冲参数初始化模块用于根据提供的自旋体系初始态和期望态分别为数据流图中的对应节点赋值；

获得目标函数模块用于利用对数据流图进行传播计算从而获得求解的最终目标函数对于优化脉冲的梯度值；

梯度迭代优化模块用于利用获得目标函数模块中的梯度值和当前的脉冲数据，对脉冲的各个节点数据进行更新；

数据处理和分析模块用于模拟数据存储于本地或在线处理，利用函数对任意自旋体系状态的演化过程进行追踪和分析；

脉冲数据优化模块用于重复执行梯度迭代优化模块及数据处理和分析模块直至目标函数取得极值，得到的脉冲数据为优化脉冲数据；

信号采集模块用于利用重复多点采集将每一自旋体系在序列片结束时的状态投影到预期信号采集状态获取自旋体系对应的复数信号值和信号采集的原始数据；所述信号采集的原始数据用于后期数据分析处理。

所述脉冲序列解析模块包括全局循环矩阵单元，所述全局循环矩阵单元包括矩阵列和矩阵行，所述矩阵列为脉冲序列的优先级等级总数；所述矩阵行为脉冲序列各优先级的循环次数的乘积同时也为需要执行的总次数；所述脉冲序列生成模块包括单一模块或高级模块，所述单一模块包括射频脉冲模块、梯度模块、信号采集模块，对所述单一模块中的任意模块通过串行和/或并行形成高级模块，对所述单一模块或高级模块设置循环优先级。

本发明的波谱成像模拟方法及系统的具体应用模拟：

模拟一、波谱学试验模拟，包括如下步骤：

1、试验环境确定，如主磁场强度、梯度场最大强度、梯度场切换率等；

2、自旋体系构建；

3、脉冲序列定义与组装，首先根据实验需求，定义序列各子块，如本示例实验为FID采集序列，包含翻转脉冲和采集模块两部分，序列组装将该两模块进行串联连接，先后执行，注意该步骤所需子块和组装方式均可根据实验需求自由定制；

4、模拟实验，通过模拟运行函数(如本发明的run{})，将上述自旋体系和脉冲序列绑定后，即可在计算机上完成模拟过程，并获取原始数据和傅里叶变换后的波谱数据。本发明示例中所有数据存储于result中，用户可通过result["raw"]["fid:re"]、result["raw"]["fid:im"]、result["raw"]["fid:abs"]分别获取原始信号的实部、虚部和幅值，也可通过result["raw"]["spec:re"]、result["raw"]["spec:im"]、result["raw"]["spec:abs"]获取信号处理后的波谱的实部、虚部和幅值。

丙烯酸体系的模拟波谱实验脚本示例,如图6(a)和图6(b),

B0{"500MHz"}//指定核磁实验的主磁场强度(单位MHz或T)

local acrylic＝spin_system{//丙烯酸自旋体系构建

spin＝"1H 1H1H",

zeeman＝"2scalar 88.42Hz

3scalar 214.9Hz",

jcoupling＝"1 2scalar 10.4Hz

1 3scalar 1.2Hz

2 3scalar 17.4Hz"

}

local rf45＝hardRF{beta＝45}//45度硬脉冲模块定义

local adc＝acq{np＝128,sw＝500}//采集模块定义

local fid＝seq{rf45,adc}//FID脉冲序列组装

result＝run{exp＝fid,spinsys＝acrylic}//实验任务模拟，指定脉冲序列和实验体系。

模拟二、磁共振成像实验仿真，包括如下步骤：

1、实验环境确定，如主磁场强度、梯度场最大强度、梯度场切换率、成像矩阵大小、成像视野(FOV)等。

2、脉冲序列构建，首先根据实验需求，定义序列各子块，如本示例实验为二维空间选择性MRI序列，包含空间选择性射频脉冲及其匹配的形状梯度、相位编码梯度、读出预散相梯度、读出梯度和延时等构成，序列组装可按图进行，如图7所示。

3、模拟实验，通过模拟运行函数(如本发明的run{})，将脉冲序列和成像模体(如本示例中的MNIBrain.h5)绑定后，即可在计算机上完成模拟过程，并获取原始数据和傅里叶变换后的成像数据。本发明示例中所有数据存储于result中，用户可通过result["raw"]["fid:re"]、result["raw"]["fid:im"]、result["raw"]["fid:abs"]分别获取原始信号的实部、虚部和幅值，也可通过result["raw"]["IMG:re"]、result["IMG"]["IMG:im"]、result["raw"]["IMG:abs"]获取信号处理后的图像的实部、虚部和幅值。如图8-9所示，空间选择性MRI成像序列的模拟实验脚本示例如下：

B0{"3T"}//指定核磁实验的主磁场强度(单位MHz或T)

peak_grad{40}--T/m//指定梯度磁场最大值

slew_rate{200}--T/m/s//指定梯度场最大切换率

seqParam{fov＝'240*240',matrix＝'128*128'}//指定成像视野和成像矩阵

local TR＝350//定义序列重复时间

local TE＝10//定义序列回波时间

//定义空间选择性梯度脉冲形状

local gxRf＝exprGrad{axis＝"X",width＝10,expr＝"D(t)1e6/gamma1H*20*(1-t/10)*Cos(2*Pi*8*t/10)"}--mT/m

local gyRf＝exprGrad{axis＝"Y",width＝10,expr＝"D(t)1e6/gamma1H*20*(1-t/10)*Sin(2*Pi*8*t/10)"}--mT/m

//定义对应射频脉冲形状

local rf＝shapedRF{width＝10,step＝100,beta＝20,pattern＝"40/10*Exp(-2^2*(1-t/10)^2)*Sqrt((2*Pi*8*(1-t/10))^2+1)"}

//定义其他所需的成像模块，如频率编码梯度、相位编码梯度、延时等。

local gy＝trapGrad{axis＝"Y",func＝"phase_encode",width＝2}

local adc＝acq{np＝128,sw＝32000}

local gx＝trapGrad{axis＝"X",func＝"read_out"}

local gxPre＝trapGrad{axis＝"X",area＝0.5*area(gx),width＝2}

local rf180＝hardRF{beta＝180,width＝0.002}

local d1＝delay{width＝TE/2-rf.tau/2-gy.tau-rf180.tau/2}

local d2＝delay{width＝TE/2-gx.tau/2-rf180.tau/2}

local d3＝delay{width＝TR-TE-gx.tau/2-rf.tau/2}

//空间选择性脉冲序列的组装

local spsel＝seq{rf+gxRf+gyRf,d1,gxPre+gy#,rf180,d2,gx+adc,d3}

//成像实验任务模拟，指定脉冲序列和实验的模体

result＝run{exp＝spsel,phantom＝"MNIBrain.h5"}

模拟三、磁共振射频脉冲的优化控制，包括如下步骤：

1、实验环境确定，如主磁场强度、梯度场最大强度、梯度场切换率等；

2、自旋体系构建；

3、优化控制器选择，首先根据实验需求，选取适当的优化控制工具，如本发明中单脉冲优化采用Opt{}，而多脉冲同步优化则采用multiOpt{}；

4、优化控制目标参数设定，根据需要，设置待优化脉冲的步数、时间步长、初始形状、自旋体系的初始状态和目标状态，以及必要时提供优化控制器所需的迭代次数、收敛条件等，即可利用计算机进行脉冲数据的迭代优化。

5、数据处理和分析，射频脉冲优化完毕，数据可存储于本地，或在线处理如图10所示，利用投影函数(如本发明的projection{})可实现对任意感兴趣自旋体系状态的演化过程追踪。进一步地，还可应用本发明的时频域分析方法对优化脉冲的数据进行时频域波谱图分析，如图11所示。

用于量子伪纯态制备的多步协作脉冲的优化实验脚本示例如下：

//丙烯酸自旋体系构建

local sys＝spin_system{

spin＝"13C 13C13C",

zeeman＝"1scalar 15.74kHz

3scalar-4.302kHz",

jcoupling＝"1 2scalar 52.2Hz

2 3scalar 35.4Hz"

}

//基于上述体系构建脉冲优化所需的数据流图

local oc＝multiOpt(sys)

//脉冲优化数据初始化、正反向传播计算及梯度迭代计算

local multi_rf＝oc:optimize{

ncoop＝7,width＝18.5,step＝1850,

init_state＝"I1z+I2z+I3z",

targ_state＝"I1z+I2z+I3z+2*I1zI2z+2*I1zI3z

+2*I2zI3z+4*I1zI2zI3z"}

如图11所示，对上述优化的7步脉冲进行时频域的波谱分析，可解析出优化射频脉冲包含了自旋体系三个碳13原子的相应化学位移信息，从而从时频域角度解释了优化脉冲在时域波形的高度调制特性。

尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用，它完全可以被适用于各种适合本发明的领域，对于熟悉本领域的人员而言，可容易地实现另外的修改，因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的实施例。

Claims (9)

1.一种磁共振波谱成像模拟方法，其特征在于，包括如下步骤：

对象生成，根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；

脉冲序列生成，若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；

脉冲序列解析，将脉冲序列中的模块的左右时间边界及其内部关键点投影到同一时间轴上，形成若干个长短不一且物理控制量恒定的序列片，对设置循环优先级的循环脉冲序列解析形成所有序列片；

模拟实验，按照时间轴先后顺序，逐次对所有序列片进行核磁共振动力学演化模拟，得到模拟数据；

在步骤脉冲序列生成中，所述模块包括单一模块或高级模块，所述单一模块包括射频脉冲模块、梯度模块、信号采集模块，对所述单一模块中的任意模块通过串行和/或并行形成高级模块，对所述单一模块或高级模块设置循环优先级。

2.如权利要求1所述的一种磁共振波谱成像模拟方法，其特征在于，在步骤脉冲序列解析中还包括，所述循环优先级相同的多个单一模块或高级模块选取循环次数最大的单一模块或高级模块作为该优先级的循环次数；利用检索对循环脉冲序列进行循环次数和循环优先级统计，生成脉冲序列的全局循环矩阵；所述全局循环矩阵包括矩阵列和矩阵行，所述矩阵列为脉冲序列的优先级等级总数；所述矩阵行为脉冲序列各优先级的循环次数的乘积，即脉冲序列需要执行的总次数。

3.如权利要求1所述的一种磁共振波谱成像模拟方法，其特征在于，在步骤模拟实验中还包括：信号采集，利用重复多点采集将每一自旋体系在序列片结束时的状态投影到预期信号采集状态获取自旋体系对应的复数信号值和信号采集的原始数据；所述信号采集的原始数据用于后期数据分析处理。

4.如权利要求1所述的一种磁共振波谱成像模拟方法，其特征在于，在步骤模拟实验后还包括步骤：

数据处理和分析，模拟数据存储于本地或在线处理，利用函数对任意自旋体系状态的演化过程进行追踪和分析。

5.如权利要求1所述的一种磁共振波谱成像模拟方法，其特征在于，在步骤脉冲序列生成之前还包括如下步骤：

S11、数据流图创建，根据定义的自旋体系演化所需的数据流节点、自旋体系的初始态和期望态节点，定义自旋体系所需的目标函数节点，同时建立数据流连接；

S12、脉冲参数初始化，根据提供的自旋体系初始态和期望态分别为数据流图中的对应节点赋值；

S13、获得目标函数，利用对数据流图进行传播计算获得求解的最终目标函数，所述目标函数用于得出优化脉冲的梯度值；

S14、梯度迭代优化，利用步骤S13中的梯度值和当前的脉冲数据，对脉冲的各个节点数据进行更新；

S15、脉冲数据优化，重复步骤S13、S14直至目标函数取得极值，得到优化脉冲数据。

6.一种电子设备，其特征在于包括：

处理器，存储器；以及程序，其中所述程序被存储在所述存储器中，并且被配置成由处理器执行，所述程序包括用于执行如权利要求1所述的方法。

7.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行如权利要求1所述的方法。

8.一种磁共振波谱成像模拟系统，其特征在于，包括对象生成模块、脉冲序列生成模块、脉冲序列解析模块、模拟实验模块、数据流图创建模块、脉冲参数初始化模块、获得目标函数模块、梯度迭代优化模块、数据处理和分析模块、脉冲数据优化模块和信号采集模块；其中，

对象生成模块用于根据样本分子结构或核磁序列扫描提取的成像模体模型，利用磁共振量子力学理论建立对应的自旋体系；

脉冲序列生成模块用于若干个模块串行组成脉冲序列，对脉冲序列中的模块设置循环优先级构建循环脉冲序列；

脉冲序列解析模块用于将脉冲序列中的模块的左右时间边界及其内部关键点投影到同一时间轴上，形成若干个长短不一且物理控制量恒定的序列片，对设置循环优先级的循环脉冲序列解析形成所有序列片；

模拟实验模块用于按照时间轴先后顺序，逐次对循环脉冲序列解析形成的所有序列片进行核磁共振动力学演化模拟，得到模拟数据；

数据流图创建模块用于根据定义的自旋体系演化所需的数据流节点、自旋体系的初始态和期望态节点，定义自旋体系所需的目标函数节点，同时建立数据流连接；

脉冲参数初始化模块用于根据提供的自旋体系初始态和期望态分别为数据流图中的对应节点赋值；

获得目标函数模块用于利用对数据流图进行传播计算从而获得求解的最终目标函数对于优化脉冲的梯度值；

梯度迭代优化模块用于利用获得目标函数模块中的梯度值和当前的脉冲数据，对脉冲的各个节点数据进行更新；

数据处理和分析模块用于模拟数据存储于本地或在线处理，利用函数对任意自旋体系状态的演化过程进行追踪和分析；

脉冲数据优化模块用于重复执行梯度迭代优化模块及数据处理和分析模块直至目标函数取得极值，得到的脉冲数据为优化脉冲数据；

信号采集模块用于利用重复多点采集将每一自旋体系在序列片结束时的状态投影到预期信号采集状态获取自旋体系对应的复数信号值和信号采集的原始数据；所述信号采集的原始数据用于后期数据分析处理。

9.如权利要求8所述的一种磁共振波谱成像模拟系统，其特征在于，所述脉冲序列解析模块包括全局循环矩阵单元，所述全局循环矩阵单元包括矩阵列和矩阵行，所述矩阵列为脉冲序列的优先级等级总数；所述矩阵行为脉冲序列各优先级的循环次数的乘积同时也为需要执行的总次数；所述脉冲序列生成模块包括单一模块或高级模块，所述单一模块包括射频脉冲模块、梯度模块、信号采集模块，对所述单一模块中的任意模块通过串行和/或并行形成高级模块，对所述单一模块或高级模块设置循环优先级。

Images