CN116430287A - 一种肺部磁共振动态成像的方法 - Google Patents
一种肺部磁共振动态成像的方法 Download PDFInfo
- Publication number
- CN116430287A CN116430287A CN202310308954.8A CN202310308954A CN116430287A CN 116430287 A CN116430287 A CN 116430287A CN 202310308954 A CN202310308954 A CN 202310308954A CN 116430287 A CN116430287 A CN 116430287A
- Authority
- CN
- China
- Prior art keywords
- motion
- magnetic resonance
- lung
- image
- pulmonary
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 210000004072 lung Anatomy 0.000 title claims abstract description 114
- 238000003384 imaging method Methods 0.000 title claims description 47
- 238000000034 method Methods 0.000 claims abstract description 71
- 230000002685 pulmonary effect Effects 0.000 claims abstract description 46
- 230000029058 respiratory gaseous exchange Effects 0.000 claims abstract description 44
- 238000009423 ventilation Methods 0.000 claims abstract description 35
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 30
- 238000012544 monitoring process Methods 0.000 claims abstract description 15
- 238000002595 magnetic resonance imaging Methods 0.000 claims abstract description 13
- 230000008569 process Effects 0.000 claims abstract description 6
- 230000000241 respiratory effect Effects 0.000 claims description 61
- 238000005070 sampling Methods 0.000 claims description 17
- 230000008859 change Effects 0.000 claims description 9
- 238000010586 diagram Methods 0.000 claims description 9
- 230000009466 transformation Effects 0.000 claims description 9
- 230000035945 sensitivity Effects 0.000 claims description 6
- 238000002939 conjugate gradient method Methods 0.000 claims description 4
- 230000003434 inspiratory effect Effects 0.000 claims description 2
- 230000006870 function Effects 0.000 description 12
- 230000007423 decrease Effects 0.000 description 7
- 238000012806 monitoring device Methods 0.000 description 6
- 238000005457 optimization Methods 0.000 description 5
- 238000012360 testing method Methods 0.000 description 5
- 229910015234 MoCo Inorganic materials 0.000 description 4
- 238000006073 displacement reaction Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 230000002123 temporal effect Effects 0.000 description 4
- 238000002591 computed tomography Methods 0.000 description 3
- 238000000338 in vitro Methods 0.000 description 3
- 230000003068 static effect Effects 0.000 description 3
- 208000019693 Lung disease Diseases 0.000 description 2
- 230000003187 abdominal effect Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000012937 correction Methods 0.000 description 2
- 238000012307 MRI technique Methods 0.000 description 1
- 206010035664 Pneumonia Diseases 0.000 description 1
- 208000010378 Pulmonary Embolism Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 238000010420 art technique Methods 0.000 description 1
- 230000000747 cardiac effect Effects 0.000 description 1
- 238000013184 cardiac magnetic resonance imaging Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009792 diffusion process Methods 0.000 description 1
- 238000013535 dynamic contrast enhanced MRI Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 230000005865 ionizing radiation Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 210000002345 respiratory system Anatomy 0.000 description 1
- 210000004872 soft tissue Anatomy 0.000 description 1
- 230000009469 supplementation Effects 0.000 description 1
- 230000001131 transforming effect Effects 0.000 description 1
- 238000012800 visualization Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves
- A61B5/055—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves involving electronic [EMR] or nuclear [NMR] magnetic resonance, e.g. magnetic resonance imaging
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/08—Detecting, measuring or recording devices for evaluating the respiratory organs
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
- G01R33/48—NMR imaging systems
- G01R33/4806—Functional imaging of brain activation
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
- G01R33/48—NMR imaging systems
- G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
- G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
- G01R33/5608—Data processing and visualization specially adapted for MR, e.g. for feature analysis and pattern recognition on the basis of measured MR data, segmentation of measured MR data, edge contour detection on the basis of measured MR data, for enhancing measured MR data in terms of signal-to-noise ratio by means of noise filtering or apodization, for enhancing measured MR data in terms of resolution by means for deblurring, windowing, zero filling, or generation of gray-scaled images, colour-coded images or images displaying vectors instead of pixels
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
- G01R33/48—NMR imaging systems
- G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
- G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
- G01R33/563—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution of moving material, e.g. flow contrast angiography
- G01R33/56308—Characterization of motion or flow; Dynamic imaging
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01R—MEASURING ELECTRIC VARIABLES; MEASURING MAGNETIC VARIABLES
- G01R33/00—Arrangements or instruments for measuring magnetic variables
- G01R33/20—Arrangements or instruments for measuring magnetic variables involving magnetic resonance
- G01R33/44—Arrangements or instruments for measuring magnetic variables involving magnetic resonance using nuclear magnetic resonance [NMR]
- G01R33/48—NMR imaging systems
- G01R33/54—Signal processing systems, e.g. using pulse sequences ; Generation or control of pulse sequences; Operator console
- G01R33/56—Image enhancement or correction, e.g. subtraction or averaging techniques, e.g. improvement of signal-to-noise ratio and resolution
- G01R33/565—Correction of image distortions, e.g. due to magnetic field inhomogeneities
- G01R33/56509—Correction of image distortions, e.g. due to magnetic field inhomogeneities due to motion, displacement or flow, e.g. gradient moment nulling
Landscapes
- Physics & Mathematics (AREA)
- Health & Medical Sciences (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- General Health & Medical Sciences (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- High Energy & Nuclear Physics (AREA)
- Radiology & Medical Imaging (AREA)
- General Physics & Mathematics (AREA)
- Condensed Matter Physics & Semiconductors (AREA)
- Signal Processing (AREA)
- Pathology (AREA)
- Heart & Thoracic Surgery (AREA)
- Surgery (AREA)
- Animal Behavior & Ethology (AREA)
- Medical Informatics (AREA)
- Public Health (AREA)
- Veterinary Medicine (AREA)
- Biophysics (AREA)
- Biomedical Technology (AREA)
- Molecular Biology (AREA)
- Pulmonology (AREA)
- Physiology (AREA)
- Vascular Medicine (AREA)
- Neurosurgery (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Magnetic Resonance Imaging Apparatus (AREA)
Abstract
一种肺部磁共振动态成像方法，包括：采用超短回波时间序列，采集自由呼吸状态下的肺部磁共振信号；在数据采集过程中，监测呼吸状态，得到呼吸曲线；利用所述呼吸曲线，将采集的肺部磁共振信数据进行运动分解重建，得到肺部运动分解图像；将所述肺部运动分解图像进行运动场估计；将运动场估计的结果进行运动状态加权的运动补偿重建，得到肺部磁共振动态图像；将所述肺部磁共振动态图像进行通气量图估计，得到肺部的通气量。
Description
技术领域
本发明涉及磁共振成像领域，尤其涉及一种肺部磁共振动态成像的方法。
背景技术
随着当代技术的发展，磁共振成像(Magnetic Resonance Imaging,MRI)技术已经展示出检查诸如肺栓塞、肺结节和新型冠状肺炎等肺部疾病的能力。与计算机断层扫描(Computed Tomography,CT)技术相比，肺部磁共振可以在无电离辐射的情况下提供良好的软组织对比度和肺部解剖结构信息，更加适用于儿童和其他需要长期随访检查的患者。然而，目前大部分肺部磁共振成像技术是借助呼吸门控来采集成像数据，这种方法扫描效率低，且只能获得单一呼吸状态的静态图像。动态磁共振成像是一项很有前景的技术，这种技术在进行肺部结构检查时扫描效率高，同时借助动态图像可以评估肺部通气功能。然而，由于肺部质子密度低、信号衰减快(T2*短)、呼吸运动以及呼吸过程中肺实质信号强度的改变等综合因素，使得动态肺部MRI具有一定难度。
为了解决肺部质子密度低、信号衰减快的问题，在肺部磁共振成像时一般选用超短回波时间(Ultrashort Time ofEcho,UTE)序列采集信号。UTE序列在射频激发后立即施加读出梯度并开始采集K空间信号，因此能够以亚毫秒级别的回波时间(echo time,TE)捕捉肺实质的短T2*信号。近期的一些研究表明，使用三维径向UTE序列可以检测一些肺部疾病，为CT扫描提供了一个可能的替代方案。尽管借助UTE序列较好解决了肺部磁共振信号采集的问题，肺部UTE成像依然需要准确的呼吸监测，达到减少图像中的呼吸伪影的问题。
在MR数据采集过程中，监测呼吸运动的方法有许多种。临床上最常用的方法是借助外部设备，如呼吸绑带，该方法通过间接测量腹部拉伸来估计呼吸运动。然而，这种方法操作过程繁琐，并且估计结果不太准确。另外一种方法是自导航(self navigation)，该方法是从采集到的K空间数据中直接估计呼吸信号，不需要额外的呼吸监测设备。在这种方法中，自导航信号可以是K空间中心点的数据，也可以是沿上下方向的导航回波(Superior-Inferior navigator，SI navigator)。此外，获取高时间分辨率、低空间分辨率的图像序列也可用于自导航。
利用检测到的呼吸信号，有多种运动校正或补偿重建方法可以用来重建出动态图像。目前最常用的方法是运动分解(motion resolved)重建，在该方法中，所有采集到的数据被分入不同的运动状态，并使用基于图像序列时空相关性的压缩感知(CompressedSensing,CS)算法重建出动态图像。这种方法已成功应用于腹部成像、心脏电影成像、动态对比度增强MRI和动态肺部UTE成像，代表性的有XD-GRASP(eXtra-Dimensional Golden-angle Radial Sparse Parallel)和XD-UTE(eXtra-Dimensional UTE)。然而，在运动位移较大的情况下，这种方法的重建结果容易产生模糊伪影。为了减少大的运动位移对运动分解重建的影响，一些研究提出了基于运动补偿(Motion Compensation,MoCo)的重建方法，并已成功应用于动态心脏MRI。近期，这些MoCo重建方法的一个变种，称为iMoCo(iterativeMoCo)UTE，已成功应用于肺部UTE MRI。在iMoCo中，使用配准得到的运动场(motion field)将参考状态下的单帧图像转换到其他运动状态来模拟呼吸运动，重建得到参考状态下的静态图像，而不是动态图像序列。然而，由于基于CS的运动分解重建和MoCo重建方法没有充分考虑肺实质信号强度与肺部磁共振动态图像的时间域相关性的变化关系，导致重建图像质量降低。
因此，本领域的技术人员致力于开发一种肺部磁共振动态成像的方法，利用运动补偿增加动态图像序列中图像间的时间相关性，减少重建伪影，并应用运动状态加权方法来适应肺实质信号强度随时间变化的先验信息，提高重建质量。
发明内容
与人体其他部位的动态成像相比，肺部磁共振动态成像具有如下特征：在呼吸气末端(End Expiration,EE)到吸气末端(End Inspiration,EI)的过程中，肺泡密度会降低，从而导致动态图像序列中肺实质区域的信号强度呈现下降趋势。因此，肺部磁共振动态图像的时间域相关性会随着肺实质信号强度的变化而降低。然而，现有已知的基于运动校正或补充实现肺部磁共振动态成像的方法(例如基于图像序列时空相关性的压缩感知(Compressed Sensing,CS)的重建方法，以及基于iMoCo(iterative MoCo)UTE的运动补偿重建方法)没有充分考虑到这种信号强度的变化，导致重建图像质量降低。
本发明提供了一种肺部磁共振动态成像方法，通过运动补偿重建算法来实现动态成像，包括三维径向UTE采集，数据采集过程中的呼吸状态监测，运动分解重建和运动场估计，以及一种新型动态MRI重建算法。该算法主要是利用运动补偿来增加动态图像序列中图像间的时间相关性，减少重建伪影；并应用运动状态加权方法来适应肺实质信号强度随时间变化的先验信息，提高重建质量。
为实现上述目的，本发明提供了一种肺部磁共振动态成像方法，包括：
步骤S10采用超短回波时间序列，采集自由呼吸状态下的肺部磁共振信号；
步骤S20获取呼吸曲线；
步骤S30利用所述呼吸曲线，将所述步骤S10中采集的数据进行运动分解重建，得到肺部运动分解图像；
步骤S40将所述肺部运动分解图像进行运动场估计；
步骤S50将所述步骤S40得到的结果进行运动状态加权的运动补偿重建，得到肺部磁共振动态图像；
步骤S60将所述肺部磁共振动态图像进行通气量图估计，得到肺部的通气量。
进一步地，在所述步骤S10中，所述超短回波时间序列的采样轨迹为交错式的螺旋叶序形式。
进一步地，所述步骤S10包括：在采样过程中，将每个叶序相较于前一个叶序的Z轴旋转一个黄金角度137.51°。
进一步地，所述步骤S10还包括：在数据采集过程中，监测呼吸状态；其中，对所述呼吸状态的监测可以使用自导航方式，在每个交错叶序采集之前，先采集一条沿着导航方向的回波信号用作导航信号；或者，对所述呼吸状态的监测也可采用外部呼吸监测设备，比如使用呼吸绑带，记录采集过程中的呼吸状态。
进一步地，所述步骤S20包括：
使用所述自导航方式监测呼吸状态时，对所述导航信号进行一维傅里叶反变换，得到所述肺部磁共振信号沿导航方向的各线圈通道的投影分布图；
从所述各线圈通道的所述投影分布图中估计出呼吸信号；
从所述呼吸信号中选出真实可信的呼吸信号；
将所述真实可信的呼吸信号的平均值作为最终估计出的呼吸曲线；或者
使用所述外部呼吸监测设备监测所述呼吸状态时，从所述外部呼吸监测设备直接获取记录到的所述呼吸曲线。
进一步地，所述步骤S30包括：
基于图像序列时空相关性的压缩感知算法重建出所述肺部运动分解图像。
进一步地，所述基于图像序列时空相关性的压缩感知算法的重建用如下公式表示：
其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，TVt和TVs分别代表时间域和空间域全变差稀疏变换，λt和λs分别代表时间域和空间域稀疏项的权重。
进一步地，利用共轭梯度法求解所述公式1。
进一步地，所述步骤S40包括：
根据所述肺部运动分解图像，将每个运动状态的图像配准到其他运动状态的图像，得到所述每个运动状态到其他所有运动状态的运动场。
进一步地，基于Demons的三维非刚性图像配准技术进行图像配准。
进一步地，所述步骤S50包括：
利用如下公式进行运动补偿重建：
其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，λs代表空间域稀疏项的权重，TVss代表空间域全变差稀疏变换，λt代表时间域稀疏项的权重，Ope()表示时间域稀疏变换，所述时间域稀疏变换将所述每个运动状态的图像与其它所有配准到这个运动状态图像的加权平均图像做差，来表示动态图像序列的时域稀疏性。
进一步地，所述Ope()的表达式如下：
其中，N是所述待求解的肺部运动分解图像X的总的运动状态数，xk是所述X中的第k个运动状态，Tik()代表将运动状态i的图像变换到状态k的运动变换算子，αik是运动状态加权系数。
进一步地，所述运动状态加权系数的表达式如下：
其中，c表示衰减系数。
进一步地，利用交错方向乘子法求解所述公式2。
进一步地，重复执行所述步骤S40和所述步骤S50。
进一步地，重复执行所述步骤S40和所述步骤S50包括：
将所述步骤S40得到的结果进行所述运动状态加权的运动补偿重建后，得到肺部磁共振动态图像序列；
借助质量更高的图像，将所述肺部磁共振动态图像序列重复执行所述步骤S40；
重复以上两个步骤直至达到迭代停止条件。
进一步地，所述步骤S60包括：
基于雅各比行列式估计所述通气量。
进一步地，所述步骤S60包括：
将所述步骤S50得到的所述肺部磁共振动态图像配准到中间运动状态，得到每个运动状态到所述中间运动状态的运动场；
计算每个所述运动场的雅各比行列式，以获得从所述中间运动状态到所述每个运动状态的肺部体积的百分比变化；
将呼吸周期内的所述雅各比行列式的曲线拟合到所述肺部体积的演化函数，然后估计所述通气量，其中，所述演化函数如下：
其中，v0为吸气末端的体积，Vr是所述肺部体积的正负峰间隔值，Tr是所述呼吸周期，n是一个决定呼吸曲线形式的正整数。
进一步地，所述步骤S60包括：
根据肺实质信号强度的时域变化估计所述通气量。
进一步地，所述步骤S60包括：
将所述步骤S50得到的所述肺部磁共振动态图像配准到中间运动状态，得到每个运动状态到所述中间运动状态的运动场；
将呼吸周期内每个像素的信号强度值随时间变化的值拟合到肺部密度的演化函数，然后估计所述通气量，其中，所述演化函数如下：
其中，i0为呼气末端的信号强度，Ir是肺部信号强度的正负峰间隔值，Tr是所述呼吸周期，n是一个决定呼吸曲线形式的正整数。
本发明提供的肺部磁共振动态成像方法，具有以下有益的技术效果：
1、利用SI导航的UTE序列采集肺部磁共振信号，无需外部呼吸检测设备的肺部自由呼吸扫描；
2、以及基于交错式螺旋叶序形式采集肺部磁共振信号，有助于减少动态图像重建时的伪影；
3、利用运动补偿来增加动态图像序列中图像间的时间相关性，减少了重建伪影；
4、应用运动状态加权方法来适应肺实质信号强度随时间变化的先验信息，提高了重建质量；
5、本申请在欠采样下能够获取高质量肺图像和通气量定量图像，缩短了扫描时间；且能提供肺部结构与通气功能的三维可视化。
以下将结合附图对本发明的构思、具体结构及产生的技术效果作进一步说明，以充分地了解本发明的目的、特征和效果。
附图说明
图1是本发明一个较佳实施例提供的肺部磁共振动态成像方法的流程图；
图2是本发明一个较佳实施例提供的带SI导航功能的UTE序列的数据采集方式图示；
图3是本发明一个较佳实施例提供的从SI导航信号中提取呼吸曲线的流程示意图；
图4是本发明一个较佳实施例提供的呼吸曲线图；
图5是本发明一个较佳实施例提供的运动分解重建示意图；
图6是本发明一个较佳实施例提供的运动场估计方式的图示；
图7是本发明一个较佳实施例提供的基于运动状态加权的运动补偿重建算法的框架图；
图8是图7的流程图；
图9是本发明一个较佳实施例提供的估计通气量的示意图；
图10是本发明一个较佳实施例提供的离体猪肺成像测试结果；
图11是本发明一个较佳实施例提供的离体猪肺通气量图估计结果；
图12是本发明一个较佳实施例提供的人体肺部自由呼吸动态成像测试结果；
图13是本发明一个较佳实施例提供的人肺通气量图估计结果。
具体实施方式
以下参考说明书附图介绍本发明的多个优选实施例，使其技术内容更加清楚和便于理解。本发明可以通过许多不同形式的实施例来得以体现，本发明的保护范围并非仅限于文中提到的实施例。
如图1所示，本发明提供的肺部磁共振动态成像方法，包括如下步骤：
S10、数据采集：可以利用带导航功能的超短回波时间(UTE)序列来进行采集。例如，利用采用带上下方向(SI)导航功能的UTE序列，在自由呼吸状态下采集肺部磁共振信号。其中UTE序列的采样轨迹为交错式的螺旋叶序形式，如图2所示，在采样过程中，采样方式为：将每个叶序相较于前一个叶序的Z轴旋转了一个黄金角度137.51°，黄金角度(GoldenAngle,GA)的定义为
其中1.618为黄金分割系数。这种采样方式保证了每个叶序采集的随机性，有助于减少动态图像重建时的伪影。此外，在每个交错叶序采集之前先采一条沿着SI方向的回波信号用作SI导航。在一些实施方式中，也可以采用外部呼吸监测设备(例如呼吸绑带)，记录采集过程中的呼吸状态。
S20、呼吸信号获取：借助步骤S10中的SI导航信号，估计采集过程中的呼吸曲线。如图3所示，对步骤S10所得各线圈通道SI导航信号(图3中的第一列)沿Kz方向进行一维傅里叶反变换，得到每个线圈通道中人体信号(即所采集的肺部磁共振信号)沿着SI导航方向的投影分布图(图3中的第二列)。采用主成分分析的算法从各线圈通道的投影分布图中估计出呼吸信号，具体做法是对各线圈通道的投影分布图进行特征值分解，将最大特征值对应的时域特征向量作为各线圈通道的呼吸信号(图3中的第三列)。在得到从各线圈通道中估计出的呼吸信号后，使用线圈聚类(coil clustering)的方法从所有的呼吸信号中选出真实可信的一部分，线圈聚类的具体做法是将所有线圈通道的呼吸曲线两两组合，计算两条曲线间的相关性，相关系数高于0.95的两条曲线就作为真实可信的呼吸曲线。如图3中第三列所示，线圈1和n的呼吸曲线相关性高，而线圈1和N的相关性低，线圈1和n的呼吸曲线就是真实可信的。最后，将这部分真实可信的呼吸曲线的平均值作为最终估计出的呼吸曲线(图3中的第四列)。图4为从某次扫描得到的SI导航信号中提取出的呼吸曲线。
在一些实施方式中，当使用外部呼吸监测设备监测呼吸状态时，从该外部呼吸监测设备直接获取记录的呼吸曲线。
S30、运动分解重建：借助步骤S20所得图4中的呼吸曲线，将步骤S10中采集到的成像数据按呼吸位移大小分入N个不同的运动状态，将同一个呼吸状态下采集到的数据分入同一个运动状态，这样共得到N个不同运动状态下的K空间数据，如图5第二列所示。此时，由于各个运动状态的K空间数据是高度欠采样的，直接对其进行非均匀快速傅里叶逆变换(inverse Non Uniform Fast Fourier Transform,iNUFFT)得到的图像有很严重的欠采样伪影，如图5第三列所示欠采样伪影图像。因此，这里使用基于图像序列时空相关性的压缩感知算法重建出运动分解图像。该重建可以用公式表示为如下优化问题：
其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道UTE序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换(NUFFT)算子。TVt和TVs分别代表时间域和空间域全变差(Total Variation,TV)稀疏变换，λt和λs分别代表时间域和空间域稀疏项的权重。该优化问题可以用共轭梯度法求解，得到图像质量一般的肺部运动分解图像，如图5第一列所示。共轭梯度法可以采用现有技术，例如参考文献“Lustig M,Donoho D,Pauly JM.Sparse MRI:The application ofcompressed sensing for rapid MR imaging.Magn Reson Med.2007；58:1182-95”所记载的。
S40、运动场估计：如图6所示，根据步骤S30所得运动分解图像，将每个运动状态的图像配准到其它(N-1)个运动状态的图像，得到每个运动状态到其它所有运动状态共计N×(N-1)个运动场，用于之后的图像重建。此步骤进行图像配准时可选用基于Demons的三维非刚性图像配准技术(基于Demons的三维非刚性图像配准技术可以采用现有技术，例如参考文献“Thirion JP.Image matching as a diffusion process:an analogy withMaxwell's demons.Med Image Anal.1998；2:243-60”)，也可选用其它常用的三维非刚性图像配准方法。
S50、运动状态加权的运动补偿重建，得到动态图像：如图7所示，使用步骤S30所得分好运动状态的成像数据与步骤S40所得运动场，进行运动补偿重建，得到肺部磁共振动态图像。运动补偿重建可以用公式表示为如下优化问题：
其中，以上公式的数据一致性项(第一项)中的X，d，D，S，F，空间域稀疏项(第二项)中的λs，TVs符号的含义与S30中所描述问题(公式1)的符号含义一致，即X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，λs代表空间域稀疏项的权重，TVs代表空间域全变差稀疏变换。时间域稀疏项(第三项)中的λt代表时间域稀疏项的权重，Ope()表示时间域稀疏变换，该变换将每个运动状态的图像与其它所有配准到这个运动状态图像的加权平均图像做差，来表示动态图像序列的时域稀疏性。Ope()的表达式(公式3)如下：
其中N是X总的运动状态数，xk是X中的第k个运动状态，Tik()代表将运动状态i的图像变换到状态k的运动变换算子。运动状态加权系数αik用下式(公式4)定义：
对于某一运动状态k，系数αik的大小随着xi和xk间时域距离的增大而减小，减小速度由衰减系数c决定。使用交错方向乘子法(Alternating Direction Method ofMultipliers,ADMM)求解运动补偿优化问题，以得到肺部磁共振动态图像。交错方向乘子法可以采用现有技术，例如参考文献“Boyd SP.Distributed optimization andstatistical learning via the alternating direction method of multipliers.NowPublishers Inc.；2011:126p”。
在一些较佳的实施方式中，可以通过重复迭代的方式提升图像重建质量。如图7和图8所示，重复执行步骤S40和步骤S50。具体地，步骤S50可以进一步包括：S501：运动状态加权的运动补偿重建，即利用公式2、3、4得到动态图像：S502；在得到动态图像序列后，借助质量更高的该图像重新估计运动场(执行步骤S40)，如S503所示。通过将该重建与运动场估计这两个步骤迭代执行，可以增加运动场估计精度并提升图像重建质量。迭代停止的条件可以根据实际需求设定，例如，设定具体的次数，或者设定图像质量达到预设条件。
S60、通气量图估计：根据步骤S50所得肺部动态图像，估计肺部通气量图。在呼吸过程中，肺泡的体积会发生改变，从动态图像可以估计出这种改变导致的运动场，运动场的雅各比行列式代表肺泡体积的百分比变化(参见图9)，用这种方式可以估计出肺部的通气量值；此外，呼吸过程中肺泡体积的改变也导致肺泡密度的改变，反映为磁共振图像中肺实质信号强度的时域变化，借助这种信号强度变化的程度大小也可以估计出肺部的通气量值，如图9所示。
基于雅各比行列式估计通气量的方法：如图9所示，将动态图像配准到一个选取的中间运动状态(如运动状态N/2)，得到每个运动状态到中间状态的运动场，其中，中间运动状态可以是第1个运动状态到第N个运动状态中的任意一个，较佳地，选取第N/2个运动状态。然后计算每个运动场的雅各比行列式(Jacobian determinant)，以获得从中间状态到每个运动状态肺部体积的百分比变化，雅各比行列式的计算方式如下(公式5)：
其中(ux,ux,ux)是配准得到的每个体素的位移大小。为了充分利用所有运动状态的信息并尽量减少噪声的影响，将呼吸周期内的雅各比行列式曲线拟合到肺部体积的演化函数(公式6)：
其中，v0为最大体积(吸气末端的体积)，Vr是肺部体积的正负峰间隔值，Tr是呼吸周期，n是一个决定呼吸曲线形式的正整数。根据该演化函数，拟合出最大体积v0以及正负峰间隔值Vr，便可以估计出肺部通气量Vr/v0。
基于肺部图像信号强度估计通气量的方法：与基于雅各比行列式的方法类似，如图9所示，将动态图像配准到一个中间运动状态，然后将呼吸周期内每个像素的信号强度值随时间变化的值拟合到肺部密度也即肺部磁共振信号强度的演化函数(公式7)：
其中，i0为最大信号强度(呼气末端的信号强度)，Ir是肺部信号强度的正负峰间隔值。根据该演化函数，拟合出最大信号强度i0以及正负峰间隔值Ir，便可以估计出肺部通气量Ir/i0。
本发明较佳实施例在3.0特斯拉的医用磁共振成像系统(uMR790,上海联影医疗)上进行了离体猪肺和人肺自由呼吸成像测试，证实了本发明前述较佳实施例进行肺部动态成像的可行性。
图10为离体猪肺成像测试结果。为了得到离体猪肺在不同呼吸状态下的全采样数据，用打气筒向猪肺中打入不同气量的空气：0毫升、200毫升、400毫升、600毫升、800毫升，然后对各个状态的猪肺做静态扫描，得到全采样数据。扫描参数如下：视场(field ofview)＝32×14×32cm3，各向同性分辨率＝0.8mm，翻转角(flip angle)＝8°，读出带宽(bandwidth)＝250Hz/pixel，回波时间(Time of Echo)＝60μs，重复时间(RepetitionTime)＝3ms。对于每一个运动状态，总共采集了220000根辐条，包含4400个交错叶序，每个叶序包含50根辐条，采集时间11分钟17秒。在得到全采样数据后，对每个运动状态的数据进行回顾式的欠采样，欠采样倍数为10倍。用不同重建方法对欠采样所得数据进行重建，重建结果如图10所示。
图10中，第一行是从全采样的K空间数据中重建出的作为参考的真实图像(GroundTruth)，第二行是第一行图像中方框区域的放大图。第三行是基于运动分解重建(XD-UTE)得到的作为参考的真实图像，第四行是第三行图像的误差图。第五行是基于iMoCo得到的作为参考的真实图像，第六行是第五行的误差图。第七行是基于本发明提供的方法得到的真实图像，第八行是第七行的误差图。从呼气末端到吸气末端，由于肺泡密度降低，可以看到肺实质的信号强度明显下降。在XD-UTE图像(第三行)中，一些气管结构和猪肺边缘存在图像模糊，在误差图(第四行)中用绿色箭头表示。在iMoCo的结果(第五行)中，这种图像模糊大大减少。然而，在iMoCo重建的图像中，肺实质信号强度的时域变化大大减少，在iMoCo的误差图(第六行)中可以看到由此导致的明显的重建误差。总的来说，与iMoCo相比，使用本发明的重建算法能够更好地保留了肺实质信号强度的时域变化，并且与XD-UTE相比，随机噪声减少，气管结构更清晰，肺部边缘更清晰。
图11是借助图10中动态猪肺图像，用基于雅各比行列式的方法估计出的猪肺通气量图。第一行分别显示了借助真实图像(Ground Truth)、XD-UTE、iMoCo和本发明重建算法图像估计出的通气量图。第二行显示由不同图像得出的通气量图与Ground Truth之间的误差。与另外两种方法相比，本发明重建算法的通气量图估计误差更小，如红色箭头所示。
图12为人肺自由呼吸成像测试结果。人体肺部成像扫描参数如下：视场(field ofview)＝40×28×30cm3，各向同性分辨率＝0.8mm，翻转角(flip angle)＝4°，读出带宽(bandwidth)＝250Hz/pixel，回波时间(Time of Echo)＝60μs，重复时间(RepetitionTime)＝3ms。对于人肺成像，总共采集了200000根辐条，包含4000个交错叶序，每个叶序包含50根辐条与1条SI导航信号，采集时间约10分钟。不同重建方法的重建结果如图12所示。
图12中，第一行是对高度欠采样的不同运动状态数据直接进行NUFFT变换得到的图像。其他几行是NUFFT、XD-UTE、iMoCo和使用本重建算法的结果放大图，放大区域如第一行中的方框所示。在NUFFT图像中，存在严重的欠采样伪影，这些伪影在XD-UTE图像中有所减少。然而，在XD-UTE图像中，一些肺部结构和膈肌区域存在图像模糊，如绿色箭头所示。在iMoCo图像中，这些残留噪声和图像模糊进一步减少，但iMoCo的图像看起来过于平滑，特别是位于红圈区域的一些肺实质结构模糊不清。相较于iMoCo，使用本发明重建算法的结果中肺实质的对比度更高，肺部结构和膈肌的清晰度与iMoCo的相当。
图13是借助图12中肺部动态图像，采用基于雅各比行列式的方法估计出的人肺通气量图。与图11中展示的猪肺通气量图评估结果一致，使用本发明重建算法的人肺通气量图要比XD-UTE更加真实。
以上详细描述了本发明的较佳具体实施例。应当理解，本领域的普通技术无需创造性劳动就可以根据本发明的构思作出诸多修改和变化。因此，凡本技术领域中技术人员依本发明的构思在现有技术的基础上通过逻辑分析、推理或者有限的实验可以得到的技术方案，皆应在由权利要求书所确定的保护范围内。
Claims (20)
1.一种肺部磁共振动态成像方法，包括：
步骤S10采用超短回波时间序列，采集自由呼吸状态下的肺部磁共振信号；
步骤S20获取呼吸曲线；
步骤S30利用所述呼吸曲线，将所述步骤S10中采集的数据进行运动分解重建，得到肺部运动分解图像；
步骤S40将所述肺部运动分解图像进行运动场估计；
步骤S50将所述步骤S40得到的结果进行运动状态加权的运动补偿重建，得到肺部磁共振动态图像；
步骤S60将所述肺部磁共振动态图像进行通气量图估计，得到肺部的通气量。
2.如权利要求1所述的肺部磁共振动态成像方法，其中，在所述步骤S10中，所述超短回波时间序列的采样轨迹为交错式的螺旋叶序形式。
3.如权利要求2所述的肺部磁共振动态成像方法，其中，所述步骤S10包括：
在采样过程中，将每个叶序相较于前一个叶序的Z轴旋转一个黄金角度137.51°。
4.如权利要求2所述的肺部磁共振动态成像方法，其中，所述步骤S10还包括：
在数据采集过程中，监测呼吸状态；其中，对所述呼吸状态的监测使用自导航方式，在每个交错叶序采集之前，先采集一条沿着导航方向的回波信号用作导航信号；或者，对所述呼吸状态的监测采用外部呼吸监测设备，记录采集过程中的所述呼吸状态。
5.如权利要求4所述的肺部磁共振动态成像方法，其中，所述步骤S20包括：
使用所述自导航方式监测所述呼吸状态时，对所述导航信号进行一维傅里叶反变换，得到所述肺部磁共振信号沿导航方向的各线圈通道的投影分布图；
从所述各线圈通道的所述投影分布图中估计出呼吸信号；
从所述呼吸信号中选出真实可信的呼吸信号；
将所述真实可信的呼吸信号的平均值作为最终估计出的呼吸曲线；或者
采用所述外部呼吸监测设备监测所述呼吸状态时，从所述外部呼吸监测设备直接获取记录到的所述呼吸曲线。
6.如权利要求1所述的肺部磁共振动态成像方法，其中，所述步骤S30包括：
基于图像序列时空相关性的压缩感知算法重建出所述肺部运动分解图像。
7.如权利要求6所述的肺部磁共振动态成像方法，其中，所述基于图像序列时空相关性的压缩感知算法的重建用如下公式表示：
其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，TVt和TVs分别代表时间域和空间域全变差稀疏变换，λt和λs分别代表时间域和空间域稀疏项的权重。
8.如权利要求7所述的肺部磁共振动态成像方法，其中，利用共轭梯度法求解所述公式1。
9.如权利要求1所述的肺部磁共振动态成像方法，其中，所述步骤S40包括：
根据所述肺部运动分解图像，将每个运动状态的图像配准到其他运动状态的图像，得到所述每个运动状态到其他所有运动状态的运动场。
10.如权利要求9所述的肺部磁共振动态成像方法，其中，基于Demons的三维非刚性图像配准技术进行图像配准。
11.如权利要求1所述的肺部磁共振动态成像方法，其中，所述步骤S50包括：
利用如下公式进行运动补偿重建：
其中，X是待求解的肺部运动分解图像，d是分好运动状态的多通道超短回波时间序列的K空间数据，D是非笛卡尔采样轨迹的密度补偿权重系数，S代表线圈灵敏度图，F是非均匀快速傅里叶变换算子，λs代表空间域稀疏项的权重，TVs代表空间域全变差稀疏变换，λt代表时间域稀疏项的权重，Ope()表示时间域稀疏变换，所述时间域稀疏变换将所述每个运动状态的图像与其它所有配准到这个运动状态图像的加权平均图像做差，来表示动态图像序列的时域稀疏性。
12.如权利要求11所述的肺部磁共振动态成像方法，其中，所述Ope()的表达式如下：
其中，N是所述待求解的肺部运动分解图像X的总的运动状态数，xk是所述X中的第k个运动状态，Tik()代表将运动状态i的图像变换到状态k的运动变换算子，αik是运动状态加权系数。
13.如权利要求12所述的肺部磁共振动态成像方法，其中，所述运动状态加权系数的表达式如下：
其中，c表示衰减系数。
14.如权利要求13所述的肺部磁共振动态成像方法，其中，利用交错方向乘子法求解所述公式2。
15.如权利要求1所述的肺部磁共振动态成像方法，其中，重复执行所述步骤S40和所述步骤S50。
16.如权利要求15所述的肺部磁共振动态成像方法，其中，重复执行所述步骤S40和所述步骤S50包括：
将所述步骤S40得到的结果进行所述运动状态加权的运动补偿重建后，得到肺部磁共振动态图像序列；
借助质量更高的图像，将所述肺部磁共振动态图像序列重复执行所述步骤S40；
重复以上两个步骤直至达到迭代停止条件。
17.如权利要求1所述的肺部磁共振动态成像方法，其中，所述步骤S60包括：
基于雅各比行列式估计所述通气量。
18.如权利要求17所述的肺部磁共振动态成像方法，其中，所述步骤S60包括：
将所述步骤S50得到的所述肺部磁共振动态图像配准到中间运动状态，得到每个运动状态到所述中间运动状态的运动场；
计算每个所述运动场的雅各比行列式，以获得从所述中间运动状态到所述每个运动状态的肺部体积的百分比变化；
将呼吸周期内的所述雅各比行列式的曲线拟合到所述肺部体积的演化函数，然后估计所述通气量，其中，所述演化函数如下：
其中，v0为吸气末端的体积，Vr是所述肺部体积的正负峰间隔值，Tr是所述呼吸周期，n是一个决定呼吸曲线形式的正整数。
19.如权利要求1所述的肺部磁共振动态成像方法，其中，所述步骤S60包括：
根据肺实质信号强度的时域变化估计所述通气量。
20.如权利要求19所述的肺部磁共振动态成像方法，其中，所述步骤S60包括：
将所述步骤S50得到的所述肺部磁共振动态图像配准到中间运动状态，得到每个运动状态到所述中间运动状态的运动场；
将呼吸周期内每个像素的信号强度值随时间变化的值拟合到肺部密度的演化函数，然后估计所述通气量，其中，所述演化函数如下：
其中，i0为呼气末端的信号强度，Ir是肺部信号强度的正负峰间隔值，Tr是所述呼吸周期，n是一个决定呼吸曲线形式的正整数。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310308954.8A CN116430287A (zh) | 2023-03-27 | 2023-03-27 | 一种肺部磁共振动态成像的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202310308954.8A CN116430287A (zh) | 2023-03-27 | 2023-03-27 | 一种肺部磁共振动态成像的方法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN116430287A true CN116430287A (zh) | 2023-07-14 |
Family
ID=87093697
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202310308954.8A Pending CN116430287A (zh) | 2023-03-27 | 2023-03-27 | 一种肺部磁共振动态成像的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN116430287A (zh) |
Cited By (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117078785A (zh) * | 2023-08-18 | 2023-11-17 | 厦门大学 | 一种快速非笛卡尔磁共振智能成像方法 |
-
2023
- 2023-03-27 CN CN202310308954.8A patent/CN116430287A/zh active Pending
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN117078785A (zh) * | 2023-08-18 | 2023-11-17 | 厦门大学 | 一种快速非笛卡尔磁共振智能成像方法 |
| CN117078785B (zh) * | 2023-08-18 | 2024-04-23 | 厦门大学 | 一种快速非笛卡尔磁共振智能成像方法 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US20200225307A1 (en) | System, method and computer-accessible medium for dynamic magnetic resonance imaging | |
| Scott et al. | Motion in cardiovascular MR imaging | |
| US9271661B2 (en) | Method for free-breathing magnetic resonance imaging using iterative image-based respiratory motion correction | |
| JP2018519909A (ja) | 動き検出を用いるｍｒ撮像 | |
| US12000918B2 (en) | Systems and methods of reconstructing magnetic resonance images using deep learning | |
| CN103502831A (zh) | 对运动中目标的磁共振成像 | |
| CN110652297A (zh) | 基于mri技术的肺功能成像处理方法 | |
| Moghari et al. | Accelerated whole‐heart MR angiography using a variable‐density Poisson‐disc undersampling pattern and compressed sensing reconstruction | |
| CN116430287A (zh) | 一种肺部磁共振动态成像的方法 | |
| He et al. | Low-rank and framelet based sparsity decomposition for interventional MRI reconstruction | |
| WO2022213666A1 (zh) | 结合k空间和图像空间重建的成像方法和装置 | |
| JP2004209084A (ja) | 核磁気共鳴を用いた検査装置 | |
| CN114236443B (zh) | 用于肺部动态通气功能快速定量评估的气体mri方法 | |
| Cao et al. | Motion-resolved and free-breathing liver MRF | |
| US8952693B2 (en) | Method for principal frequency magnetic resonance elastography inversion | |
| Qiu et al. | Sliding motion compensated low-rank plus sparse (SMC-LS) reconstruction for high spatiotemporal free-breathing liver 4D DCE-MRI | |
| WO2022212244A1 (en) | Distortion-free diffusion and quantitative magnetic resonance imaging with blip up-down acquisition of spin- and gradient-echoes | |
| Huang et al. | Optimization of PROPELLER reconstruction for free‐breathing T1‐weighted cardiac imaging | |
| Papp et al. | Deep learning for improving ZTE MRI images in free breathing | |
| US11294009B2 (en) | Cartesian sampling for dynamic magnetic resonance imaging (MRI) | |
| Aggarwal et al. | Spatio-temporal modeling and adaptive acquisition for cardiac MRI | |
| Wild et al. | Dynamic imaging of lung ventilation and gas flow with hyperpolarized gas MRI | |
| Balasch | Functional lung magnetic resonance imaging with ultra-short echo-time (UTE) | |
| EP4290265A1 (en) | Computer-implemented method, computer program and processing apparatus for reconstructing a dynamic series of magnetic resonance images | |
| CN118212317A (zh) | 用于肺动态通气功能定量评估的高分辨磁共振成像方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination |