CN105785298A - High-precision three-dimensional chemical shift imaging method - Google Patents

Abstract

Disclosed in the invention is a high-precision three-dimensional chemical shift imaging method. A water image and a fat image that are separated from each other and in-phase graphs and anti-phase graphs of the water image and the fat image can be obtained by scanning once by using a three-dimensional single excited long echo string collection way, a radio frequency refocusing pulse string or frequency coding gradient technology, a pre scanning scheme having an echo amplitude and phase error precision correction feature, and a data post processing method. The method has the low requirements on magnetic field uniformity and gradient system performances; and the image quality and the diagnosis value can be improved obviously.

Description

A kind of high-precision three-dimensional chemical shift imaging method

Technical field

The present invention relates to medical nuclear magnetic resonance imaging technology and Image Diagnosis technical field, particularly relate to a kind of high-precision three-dimensional chemical shift imaging method.

Background technology

In medical magnetic resonance imaging (MRI), press fat imaging that pathological changes is shown and be substantially better than conventional imaging technique, for instance, be conducive to improving the Sensitivity and Specificity of diagnosis articular cartilage damage, mammary gland and the contrast of optic nerve image can be improved.In t1 weighted image, fat appears as high signal owing to T1 is short and is prone to obscure with the focus of contrast medium enhancing in enhanced ct scans, have impact on the clinical definite meaning of enhanced ct scans.Additionally, it it is 210 hertz at the esterified displacement study difference on the frequency of 1.5T water, RF bandwidth is mostly in 1～3 kilohertz range, fatty scanning aspect spatially can move the part of thickness, several pixel can be covered at reading gradient bandwidth 50～200 hertz/pixel esterified chemical shift artefacts of Water Under in frequency coding direction simultaneously, cause human anatomic structure image distortion.Just because of this, pressure fat technology has important using value in clinical diagnosis medical science, conventional pressure fat method selective fat magnetic saturation, and water selective excites and recovers (STIR) with short TI.But, selective pulse is to magnetic field B₀With radio-frequency field B₁Inhomogeneities very sensitive, in low field owing to water fat chemical shift difference is only tens hertz and inapplicable.STIR is to B₁The inhomogeneities of field is sensitive, also reduces the signal of other tissue, and extend sweep time while suppressing fat.On the other hand, based on the chemical shift f of different tissues composition, and carry out the technology of the independent imaging of each composition to radio-frequency field B₁Uniformity and magnetic field intensity do not have particular/special requirement not lose signal noise ratio (snr) of image yet, it is also possible to measure fat and the relative scale of water in tissue, more valuable in clinical diagnosis.It is one of chemical shift imaging technology that Dixon water fat separates, can every layer of k-space data gathering two width (or three width) water fat signal phase difference respectively 0 π and (or 0, π and 2 π), then pass through data and process and obtain water images and fat image.Relatively broad based on 2 of spin echo and 3 Dixon separate imaging of water and fat technology application clinically, water fat homophase echo and anti-phase echo gather respectively in twice or three independent scanning, sweep time is the twice or three times of conventional T1 weighted scanning, it reduces clinical scanning efficiency and adds sufferer cooperation difficulty.As fast spin echo or echo planar sequence, Dixon separate imaging of water and fat technology can also adopt the mode of the multiple echo of single-shot to shorten sweep time several times.But; when magnet and gradient system performance are undesirable, not only the phase place of echo-signal usually can exceed-π to+π scope, causes serious phase place to be wound around artifact; and the phase symbol of the tissue regions water fat coexisted at water fat is prone to mistake in computation, cause that water fat separates not exclusively；Especially since both forward and reverse directions frequency encoding gradient is uneven, the use of the single-shot acquisition mode carrying out frequency coding based on bipolarity (i.e. positive-negative polarity) gradient can substantially increase phase error.These problems have become single-shot mode and have realized the technical bottleneck of conventional application in a clinical setting, reduce the diagnostic value of water fat separate picture.It addition, 2 and 3 Dxion imaging modes constrain the further raising of signal noise ratio (snr) of image or scan efficiency, the two-dimensional scan mode that conventional imaging adopts also limit the image resolution ratio selecting layer direction.

Summary of the invention

For this, the invention provides a kind of high-precision three-dimensional chemical shift imaging method, the method adopts three-dimensional single-shot long echo string acquisition mode and can obtain the water images and fat image that are kept completely separate and the same phasor of the two and anti-phase figure in conjunction with specially designed radio frequency reunion train of pulse or frequency encoding gradient technology in single pass with the prescan scheme and data post processing method with echo amplitude and the accurate correction feature of phase error, this technology is low to magnetic field homogeneity and gradient system performance requirement, can significantly improve picture quality and diagnostic value.

A kind of high-precision three-dimensional chemical shift imaging method provided by the invention, comprises the following steps:

The first step: first 180 ° of reunion pulse is applied to the slice selective gradient centre position to first echo, i.e. TE/2 place, after phase encoding gradient is all applied to first 180 ° of reunion pulse, a series of 180 ° of reunion pulses subsequently are applied to Δ t/2 place, Δ t is set to minima and Δ τ=1/ Δ f/2；Selecting gradient to be applied behind an opposite polarity phase place reunion gradient to avoid chunk to select gradient to cause magnetization vector phase dispersion at chunk, each sequence repeat period end applies the damage phase gradient of an amplitude change at random；In sequential parameter table, arranging TE is minima, and to arrange sampled data points be Dim1, and phase code step number is Dim2, and selecting layer direction phase code step number is Dim3；Above-mentioned imaging sequence module is run, successively to G in single-shot mode_p1And G_P2Amplitude carries out phase code circulation, and often step circle collection N group homophase and anti-phase echo are until all phase code step numbers complete, and thus builds four-dimensional complex matrix (Dim1 × Dim2 × 2N × Dim3)；Wherein, G_p1For selecting layer direction phase encoding gradient, G_p2For phase encoding gradient, Δ t is time delay, and Δ f is water fat Chemical shift differences, and TE is the echo time.

Second step: to four-dimensional complex matrix (Dim1 × Dim2 × 2N × Dim3) along selecting layer direction to carry out the three-dimensional complex matrix (Dim1 × Dim2 × 2N) of one-dimensional discrete Fourier transformation acquisition Dim3, extract N group two-dimensional complex number matrix (Dim1 × Dim2) more successively and carry out two dimension discrete fourier transform, it is thus achieved that N group homophase and anti-phase image；For the signal that two dimension version imaging sequence gathers, Dim3=1 is set.

3rd step: by following formula to same phasor complex matrix S₀With anti-phase figure complex matrix S₁Carry out data process:

$S_0=\underset{m=1}{\overset{N}{\Σ}}{S}_0^m=\underset{m=1}{\overset{N}{\Σ}}(S_w{e}^{-\[TE+(2m-2)\·\mathrm{\Δ}t\]/T_2^w}+S_f{e}^{-i(2\mathrm{\π}\mathrm{\Δ}f\mathrm{\Δ}\mathrm{\τ}\·(2m-2))}{e}^{-\[TE+(2m-2)\·\mathrm{\Δ}t\]/T_2^f})e^{-(2m-2)\·\mathrm{\Δ}\mathrm{\τ}/T_2^{*}}{e}^{-i({\mathrm{\φ}}_0+\mathrm{\φ}\·(2m-2\left)\right)}---\left(1\right)$

$S_1=\underset{m=1}{\overset{N}{\Σ}}{S}_1^m=\underset{m=1}{\overset{N}{\Σ}}(S_w{e}^{-\[TE+(2m-1)\·\mathrm{\Δ}t\]/T_2^w}+S_f{e}^{-i(2\mathrm{\π}\mathrm{\Δ}f\mathrm{\Δ}\mathrm{\τ}\·(2m-1))}{e}^{-\[TE+(2m-1)\·\mathrm{\Δ}t\]/T_2^f})e^{-(2m-1)\·\mathrm{\Δ}\mathrm{\τ}/T_2^{*}}{e}^{-i({\mathrm{\φ}}_0+\mathrm{\φ}\·(2m-1))}---\left(2\right)$

Wherein S_wAnd S_fIt is expressed as water and fat components in region, water fat homophase that N gathers after representing disposable radio-frequency drive and the group number of anti-phase echo, and N can strengthen according to signal noise ratio (snr) of image and weighting scheme need to take various different natural number,WithThe respectively T2 constant of the tissue of human body enrichment water and fat,It is and Magnetic field inhomogeneity degree Δ B₀Relevant apparent T2 constant, Δ τ is the chemical potential transposition evolution time, φ₀Being the initial phase of proton magnetization vector, φ is the phase error that the effect such as magnetic field bump and body area's susceptibility produces during a Δ τ.

Consider that magnetic field homogeneity is insufficient to ideal on home equipment, in above formulaWithItem can be ignored, therefore above formula can be reduced to when Δ τ=1/ Δ f/2:

$S_0=\underset{m=1}{\overset{N}{\Σ}}{S}_0^m=\underset{m=1}{\overset{N}{\Σ}}(S_w+S_f)e^{-(2m-2)\·\mathrm{\Δ}\mathrm{\τ}/T_2^{*}}{e}^{-i({\mathrm{\φ}}_0+\mathrm{\φ}\·(2m-2))}---\left(3\right)$

$S_1=\underset{m=1}{\overset{N}{\Σ}}{S}_1^m=\underset{m=1}{\overset{N}{\Σ}}(S_w-S_f)e^{-(2m-1)\·\mathrm{\Δ}\mathrm{\τ}/T_2^{*}}{e}^{-i({\mathrm{\φ}}_0+\mathrm{\φ}\·(2m-1))}---\left(4\right)$

The free induction decay signal gathered by nonlinear fitting T2 cycle tests or a succession of water fat echo-signal, to single exponential function, are obtained respectivelyWithApproximation, and obtain the same phasor complex matrix after echo amplitude correction and anti-phase figure complex matrix based on formula (1)-(2) or formula (3)-(4).

4th step: same phasor complex matrix and anti-phase figure complex matrix are carried out phase unwrapping and phasing, described in comprising the following steps that:

By rightIt is multiplied by after taking complex conjugateAnd divided byMould, it is thus achieved that there is no start-phase φ₀The same phasor complex matrix of impact | S₀| withAnd anti-phase figure complex matrix

Based onCalculate phase placeHere atan2 () is four-quadrant arctan function, or preferential based on same phasorCalculate phase placeAnd adopt conventional Branch cut (branchcut) or region growth method (regiongrowth) φ to carry out phase unwrapping to obtain real phasing matrix φ `.

RightIt is multiplied by e^{iφ`</sup>,It is multiplied by e<sup>iTM2φ`}WithIt is multiplied by e^iTM3φ`Eliminate phase error respectively.

At N > 2, other each group carries out phasing with all available similar fashion of phasor and anti-phase figure, namelyIt is multiplied by e^{i·(2m-2)φ`</sup>,It is multiplied by e<sup>i·(2m-1)φ`}。

Based on φ `=γ Δ B₀(2m-1) TM Δ τ obtains field pattern Δ B₀, and the cosine value defining the phasing matrix of the anti-phase figure of m-th is correction factor matrix κ_mFor determining that in anti-phase figure, the pixel of moisture fat signal should belong to water images or fat image.

5th step: water fat signal separation module and image reconstruction module produce, to all same phasor after phasing and anti-phase figure are cumulative respectively, the S that signal to noise ratio strengthens₀And S₁,WithIn the insignificant situation of item, produce water according to the following formula as S_wWith fat as S_f:

S_w=(| S₀|+κ_m·|S₁|/A)/2(5)

S_f=(| S₀|-κ_m·|S₁|/A)/2(6)

Or calculate the water producing to be sufficiently separated according to the following formula as S_wWith fat as S_f:

$S_w=\underset{m=1}{\overset{N}{\&Sigma;}}\left(\right|S_0^m|/A^{2m-2}+{\mathrm{\&kappa;}}_m\&CenterDot;|S_1^m|/A^{2m-1})/2---\left(7\right)$

$S_f=\underset{m=1}{\overset{N}{\&Sigma;}}\left(\right|S_0^m|/A^{2m-2}-{\mathrm{\&kappa;}}_m\&CenterDot;|S_1^m|/A^{2m-1})/2---\left(8\right)$

In above formula WithIt is the same phasor after eliminating phase error or anti-phase figure.

Wherein, take a step forward in the first step and include pre-scanning process: at the slice selective gradient of three orthogonal directions and position the test zone of human body under three corresponding soft pulse effects, then excite and gather free induction decay signal with pulse or by a succession of equally spaced 180 ° of pulse excitations and in a succession of water fat echo-signal of centric acquisition in 180 ° of pulse spacings.

Wherein, after the 5th step, the process of content and the distribution calculating fat is farther included: based on water as S_wWith fat as S_fThe meansigma methods of respective pixel or voxelWithCalculate content and the distribution of fat, it may be assumed that

High-precision three-dimensional chemical shift imaging method and the two dimension version thereof of the present invention provide general theoretical model and data analysing method, can effectively overcome the various range errors of the hardware imperfection generation of MRI machine and the impact of phase place artifact, and significantly improve signal noise ratio (snr) of image and scan efficiency, be conducive to improving the sensitivity of diagnosis or specificity and being applicable to the quantitative analysis of fat in the Differential Diagnosis of the pathological changes such as osteoarthrosis, mammary gland, optic nerve and adrenal gland and fatty liver.

Accompanying drawing explanation

Fig. 1 high-precision three-dimensional of the present invention chemical shift imaging method scan operation flow chart.

Fig. 2 high-precision three-dimensional of the present invention chemical shift imaging method vortex field cycle tests.

Wherein, gradient pulse amplitude and imaging sequence gradient are at the same order of magnitude, and list time delay is in 0.1ms and 50ms scope, and the data point of magnetic resonance signal FID is 256.

Fig. 3 high-precision three-dimensional of the present invention chemical shift imaging method local T2 constant (T₂) cycle tests.

Wherein, G_X、G_YAnd G_ZRepresenting the gradient of three orthogonal directions, its amplitude can be adjusted according to the bulk of water fat imaging region；Δ t represents the interval between 90 ° of soft pulse, and 2 τ represent the interval between adjacent 180 ° of pulsus durus punching, and the peak value of each echo of receiver collection is positioned at the midpoint during each 2 τ, and the data point of each echo is set to 128.

The apparent T2 constant (T of Fig. 4 high-precision three-dimensional of the present invention chemical shift imaging method local₂ ^*) cycle tests.

Wherein, G_X、G_YAnd G_ZRepresenting the gradient of three orthogonal directions, its amplitude and polarity can be adjusted by positioning picture according to the locus of water fat imaging region and size；Δ t represents the interval between 90 ° of soft pulse, and receiver starts to gather free induction decay signal at the 3rd Δ t end, and sampled data points is set to 512.

Fig. 5 high-precision three-dimensional of the present invention chemical shift imaging method three-dimensional single-shot long echo string chemical shift imaging sequence.

Wherein, 90 ° of excitation pulses are soft pulse (such as Sinc pulse), and 180 ° of reunion pulses are pulsus durus punching, and TE is the shortest echo time, and Δ t is the twice of the shortest interval time between echo summit and 180 ° of reunion pulses, Δ τ=1/ Δ f/2.Signals collecting is N group homophase and reverse phase gradient echo, is followed successively by in-phase signal, inversion signal, in-phase signal, inversion signal, and the rest may be inferred by analogy for it.Position in dashed box can need to repeat N-1 time by image parameter weighting.Gradient pre-emphasis waveform can switch in real time at Δ t/2 period.

Fig. 6 high-precision three-dimensional of the present invention chemical shift imaging method three-dimensional single-shot double echo chemical shift imaging sequence.

Wherein, 90 ° of excitation pulses are soft pulse (such as Sinc pulse), and 180 ° of reunion pulses are pulsus durus punching, and TE is the shortest echo time, Δ τ=1/ Δ f/2.The damage phase gradient in slice selective gradient direction is used for accelerating Spin System and returns to equilibrium state.Signals collecting gtadient echo is followed successively by water fat in-phase signal and water fat inversion signal.Gradient pre-emphasis waveform can switch in real time at Δ t/2 period.

Fig. 7 high-precision three-dimensional of the present invention chemical shift imaging method three-dimensional single-shot long echo string chemical shift imaging sequence.

Wherein, 90 ° of excitation pulses are soft pulse (such as Sinc pulse), and TE is the shortest echo time, and Δ τ is the interval time between echo summit, and Δ τ=1/ Δ f/2.Signals collecting is N group homophase and reverse phase gradient echo, is followed successively by in-phase signal, inversion signal, in-phase signal, inversion signal, and the rest may be inferred by analogy for it.Position in dashed box can need to repeat N-1 time by image parameter weighting.Gradient pre-emphasis waveform can switch in real time at Δ t/2 period.

Fig. 8 high-precision three-dimensional of the present invention chemical shift imaging method two dimension single-shot long echo string chemical shift imaging sequence.

Wherein, 90 ° of excitation pulses are soft pulse (such as Sinc pulse), and 180 ° of reunion pulses are Sinc pulse, and TE is the shortest echo time, and Δ t is the twice of the shortest interval time between echo summit and 180 ° of reunion pulses, Δ τ=1/ Δ f/2.Signals collecting is N group homophase and reverse phase gradient echo, is followed successively by in-phase signal, inversion signal, in-phase signal, inversion signal, and the rest may be inferred by analogy for it.Position in dashed box can need to repeat N-1 time by image parameter weighting.Gradient pre-emphasis waveform can switch in real time at Δ t/2 period.

Fig. 9 high-precision three-dimensional of the present invention chemical shift imaging method two dimension single-shot long echo string chemical shift imaging sequence.

Wherein, 90 ° of excitation pulses and 180 ° of reunion pulses are soft pulse (such as Sinc pulse), and TE is the shortest echo time, and Δ τ is the interval time between echo summit, and Δ τ=1/ Δ f/2.Signals collecting is N group homophase and reverse phase gradient echo, is followed successively by in-phase signal, inversion signal, in-phase signal, inversion signal, and the rest may be inferred by analogy for it.Position in dashed box can need to repeat N-1 time by image parameter weighting.Gradient pre-emphasis waveform can switch in real time at Δ t/2 period.

Detailed description of the invention

One, the feature of the present invention

Instant invention overcomes hardware imperfection and cause echo amplitude and phase error, and strengthen the three dimensional chemical displacement imaging of signal noise ratio (snr) of image, for quickly obtaining accurate magnetic resonance chemical shift image in MR imaging apparatus.

Single-shot long echo string sequence is adopted to realize a series of homophase and the Quick Acquisition of anti-phase echo-signal, and a general theoretical model is proposed, the accurately correction and the signal noise ratio (snr) of image that namely realize echo amplitude and phase place based on formula (1) and formula (2) strengthen.

Adopt three dimensions coded system to improve and select layer directional resolution, and the 180 ° of pulsus durus punchings adopting 2N pulse width very short reduce the loss of signal to try one's best, after each 180 ° of pulses and under the polarity maintenance uniform condition of frequency encoding gradient, gather N group homophase and anti-phase echo and build four-dimensional k-space matrix.

Three dimensions coded system adopts chunk to select gradient, thereafter a phase place reunion gradient is increased to make up the loss of signal, for this phase encoding gradient with after selecting layer direction phase encoding gradient to be placed in the punching of 180 ° of pulsus durus, and adopt 2N reading gradient to carry out frequency coding, wherein N may be configured as 1,2,3 or much higher value according to image weighting scheme.

Single-shot long echo string sequence is to may be used without 2N-1 opposite polarity frequency encoding gradient to replace 180 ° of pulsus durus punchings when N is bigger, and this gradient keeps integral area consistent with normal frequency encoding gradient but the amplitude significantly larger time shorter echo-signal that reduces to try one's best is lost and avoided introducing extra phase error.

N group homophase and anti-phase echo are that first echo time TE is set to the shortest time, then the interval between all adjacent echo summits is set to Δ τ=1/ Δ f/2 according to the chemical shift difference of water fat, then includes phase place reunion time Δ t under adopting 180 ° of reunion pulse situation.

Echo amplitude correction is to adopt two kinds of custom-designed local T2 cycle testss in a case of n=1, by the slice selective gradient of three orthogonal directions the test zone positioning human body under three corresponding soft pulse effects, excite and gather free induction decay signal again with pulse or by a succession of equally spaced 180 ° of pulse excitations and at the centric acquisition succession of echoes in 180 ° of pulse spacings, then nonlinear fitting is carried out by single exponential function, respectively obtain the decay for compensation water fat echo amplitude of T2 constant and apparent T2 constant.

Echo amplitude correction is to pass through correction factor when N=2Carry out amplitude correction, or select the k-space centrage (k that signal to noise ratio is the highest_y=0) and delivery, namelyWithAnd calculating obtains correction factorWithRespectively second to the 4th echo is carried out amplitude correction, and correction factor can be derived based on formula (1) and formula (2) by similar fashion when N takes much higher value.

Phase of echo correction is by obtaining a series of two dimension k-space data selecting layer direction that N group homophase and anti-phase echo-signal carry out one-dimensional discrete Fourier transform, then adopts Branch cut or region growth method to carry out two-dimensional phase solution to twine the phasing of the acquisition real phase ` scattergram water fat homophase for the same aspect of same position and anti-phase echo-signal.

Water fat homophase after phasing and anti-phase echo-signal by a general theoretical model, N group water that namely formula (5) and (6) obtain water and fat is sufficiently separated and the respective image of fat, and add up and realize signal noise ratio (snr) of image enhancing.

Adopt Sinc reunion pulse or 2N-1 opposite polarity frequency encoding gradient can realize accurately quickly two-dimensional chemical displacement imaging by same theoretical model and similar data processing method when replacing 180 ° of pulsus durus to rush.

It is applied in the MR imaging apparatus of various field intensity, according to particular field strength, sequential parameter Δ τ=1/ Δ f/2 is set, and based on general theoretical model, i.e. formula (1) formula (6), carry out similar data and process and obtain the image of water fat clean cut separation.

Two, principles of the invention

Medical magnetic resonance imaging instrument is mainly made up of hardware cells such as magnet, spectrometer, control station, gradient coil, radio-frequency coil, radio-frequency (RF) power amplification and gradient power amplifiers, and installs the scanning process work shown in Fig. 1 of three dimensional chemical displacement image-forming module and data processing module.First, according to prescan need on sequencer, to load the conventional gradients preemphasis sequence shown in Fig. 2 and debug vortex field compensating parameter carry out gradient waveform correction, then load custom-designed cycle tests here respectively and gather water fat signal for local water fat relaxation time mensuration, as shown in Figure 3 and Figure 4.Then, sequencer loads the three dimensional chemical displacement imaging sequence module shown in Fig. 5 (or Fig. 6 to Fig. 7) control each hardware cell and realize the exciting of homophase and anti-phase proton signal, space encoding and collection, and the data processing module that installation and operation is supporting on control station main frame, including echo amplitude correction module, phase of echo solution twines and correction module, water fat signal separation module and image reconstruction module, thus quickly accurately realizing chemical shift imaging function.

Here, cycle tests shown in Fig. 3 and Fig. 4 includes: at the slice selective gradient of three orthogonal directions and position the test zone of human body under three corresponding soft pulse effects, then excites and gather free induction decay signal with pulse or by a succession of equally spaced 180 ° of pulse excitations and in a succession of water fat echo-signal of centric acquisition in 180 ° of pulse spacings.

Three dimensional chemical displacement imaging sequence module includes chunk and selects gradient G_s, select layer direction phase encoding gradient G_p1, phase encoding gradient G_p2, 2N frequency encoding gradient G_rAnd the unit such as 2N echo acquirement, as shown in Figure 5.Wherein, first 180 ° of reunion pulse is applied to the slice selective gradient centre position to first echo, i.e. TE/2 place in Fig. 5, after phase encoding gradient is all applied to first 180 ° of reunion pulse, a series of 180 ° of reunion pulses subsequently apply Δ t/2 place in Figure 5, Δ t is set to minima and Δ τ=1/ Δ f/2；Further, 180 ° of reunion pulses here all adopt the very short pulsus durus punching of width to reduce the notable decay that echo amplitude causes due to transverse relaxation.Here the occupation mode of pulse of meeting again is different from common fast spin echo (FSE) situation, chemical shift effect and field inhomogeneous broadening effect and is not required to meet again, be used for realizing water fat exactly and separate and obtain Distribution of Magnetic Field figure when CPMG.Additionally, selecting gradient to be applied behind an opposite polarity phase place reunion gradient to avoid chunk to select gradient to cause magnetization vector phase dispersion at chunk, each sequence repeat period end applies the damage phase gradient of an amplitude change at random to avoid the interference to sequence follow-up operation of the residual magnetization vector.In sequential parameter table, arranging TE is minima, and to arrange sampled data points be Dim1, and phase code step number is Dim2, and selecting layer direction phase code step number is Dim3.Above-mentioned imaging sequence module is run, successively to G in single-shot mode_p1And G_P2Amplitude carries out phase code circulation, and often step circle collection N group homophase and anti-phase echo are until all phase code step numbers complete, and thus builds four-dimensional complex matrix (Dim1 × Dim2 × 2N × Dim3).As N=2, this sequence is exactly a kind of three-dimensional version 2 Dixon sequences based on spin echo, as shown in Figure 6.As N > 8 time, subsequent echoes string is it may happen that decay by a relatively large margin, when gradient intensity is big and switching rate is fast, available opposite polarity frequency encoding gradient replaces 180 ° of pulsus durus punchings in Fig. 5, this is equivalent to common echo wave plane (EPI) odd number or even number echo acquirement mode, but share the difference is that all echo of single-shot here and often walk phase code, and opposite polarity frequency encoding gradient used is keeping that integral area is consistent with normal frequency encoding gradient and its amplitude is significantly larger when the Δ f/2 of Δ τ=1/ decays with echo-signal during reducing Δ τ, as shown in Figure 7.Sequence shown in Fig. 8 and Fig. 9 is the two-dimentional version of sequence shown in Fig. 5 and Fig. 7, and all 180 ° of pulsus durus punchings are replaced by sinc soft pulse, can be selected for three lobes, single-lobe or half lobe type.

Data processing module is mainly made up of pretreatment module, amplitude correction module, phase unwrapping and correction module, water fat signal separation module and image reconstruction module, for the N group of sequence acquisition shown in Fig. 5 is carried out amplitude correction, phase unwrapping and correction and cumulative with phasor and anti-phase figure, it is derived from one group with phasor complex matrix S₀With anti-phase figure complex matrix S₁, specific implementation is as described below:

First, four-dimensional complex matrix edge is selected layer direction to carry out one-dimensional discrete Fourier transformation and is obtained Dim3 three-dimensional complex matrix (Dim1 × Dim2 × 2N) by pretreatment module, extract N group two-dimensional complex number matrix (Dim1 × Dim2) more successively and carry out two dimension discrete fourier transform, it is thus achieved that N group homophase and anti-phase image.For the signal that two dimension version imaging sequence gathers, Dim3=1 is set.

Secondly, amplitude correction module presses following formula to same phasor complex matrix S₀With anti-phase figure complex matrix S₁Carry out data process:

$S_0=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_0^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w{e}^{-\&lsqb;TE+(2m-2)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^w}+S_f{e}^{-i(2\mathrm{\&pi;}\mathrm{\&Delta;}f\mathrm{\&Delta;}\mathrm{\&tau;}\&CenterDot;(2m-2))}{e}^{-\&lsqb;TE+(2m-2)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^f})e^{-(2m-2)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-2\left)\right)}---\left(1\right)$

$S_1=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_1^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w{e}^{-\&lsqb;TE+(2m-1)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^w}+S_f{e}^{-i(2\mathrm{\&pi;}\mathrm{\&Delta;}f\mathrm{\&Delta;}\mathrm{\&tau;}\&CenterDot;(2m-1))}{e}^{-\&lsqb;TE+(2m-1)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^f})e^{-(2m-1)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-1))}---\left(2\right)$

Above formula provides the general mathematical model fully describing this chemical shift imaging signal, wherein a S_wAnd S_fIt is expressed as water and fat components in region, water fat homophase that N gathers after representing disposable radio-frequency drive and the group number of anti-phase echo, and N can strengthen according to signal noise ratio (snr) of image and weighting scheme need to take various different natural number,WithThe respectively T2 constant of the tissue of human body enrichment water and fat,It is and Magnetic field inhomogeneity degree Δ B₀Relevant apparent T2 constant, Δ τ is the chemical potential transposition evolution time, φ₀It is the initial phase of proton magnetization vector, by gradient, receives link group delay, vortex field, Maxwell field and radio-frequency field B₁The impact of phase place, φ is the phase error that the effect such as magnetic field bump and body area's susceptibility produces during a Δ τ.Consider that magnetic field homogeneity is insufficient to ideal on home equipment, in above formulaWithItem can be ignored, therefore above formula can be reduced to when Δ τ=1/ Δ f/2:

$S_0=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_0^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w+S_f)e^{-(2m-2)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-2))}---\left(3\right)$

$S_1=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_1^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w-S_f)e^{-(2m-1)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-1))}---\left(4\right)$

The data of same phasor and anti-phase figure that formula (3) and formula (4) are equally applicable to sequence acquisition shown in Fig. 6 to Fig. 9 process.

Free induction decay signal that amplitude correction module is gathered by local T2 cycle tests shown in nonlinear fitting Fig. 3 and Fig. 4 or a succession of water fat echo-signal, to single exponential function, obtain respectivelyWithApproximation, and obtain the same phasor complex matrix after echo amplitude correction and anti-phase figure complex matrix based on formula (1)-(2) or formula (3)-(4).

Then, phase unwrapping and correction module are for carrying out phase unwrapping and phasing to same phasor complex matrix and anti-phase figure complex matrix, and its implementation is as described below:

By rightIt is multiplied by after taking complex conjugateAnd divided byMould, it is thus achieved that there is no start-phase φ₀The same phasor complex matrix of impact | S₀| withAnd anti-phase figure complex matrix

Based onCalculate phase placeHere atan2 () is four-quadrant arctan function, or preferential based on same phasorCalculate phase place

And adopt conventional Branch cut (branchcut) or region growth method (regiongrowth) φ to carry out phase unwrapping to obtain real phasing matrix φ `.

RightIt is multiplied by e^{iφ`</sup>,It is multiplied by e<sup>i·2φ`}WithIt is multiplied by e^i·3φ`Eliminate phase error respectively.

At N > 2, other each group carries out phasing with all available similar fashion of phasor and anti-phase figure, namelyIt is multiplied by e^{i·(2m-2)φ`</sup>,It is multiplied by e<sup>i·(2m-1)φ`}。

Based on φ `=γ Δ B₀(2m-1) Δ τ obtains field pattern Δ B₀, and the cosine value defining the phasing matrix of the anti-phase figure of m-th is correction factor matrix κ_mFor determining that in anti-phase figure, the pixel of moisture fat signal should belong to water images or fat image.

Finally, water fat signal separation module and image reconstruction module produce, to all same phasor after phasing and anti-phase figure are cumulative respectively, the S that signal to noise ratio strengthens₀And S₁,WithWater is produced according to the following formula as S in the insignificant situation of item_wWith fat as S_f:

S_w=(| S₀|+κ_m·|S₁|/A)/2(5)

S_f=(| S₀|-κ_m·|S₁|/A)/2(6)

Or calculate the water producing to be sufficiently separated according to the following formula as S_wWith fat as S_f:

$S_w=\underset{m=1}{\overset{N}{\&Sigma;}}\left(\right|S_0^m|/A^{2m-2}+{\mathrm{\&kappa;}}_m\&CenterDot;|S_1^m|/A^{2m-1})/2---\left(7\right)$

$S_f=\underset{m=1}{\overset{N}{\&Sigma;}}\left(\right|S_0^m|/A^{2m-2}-{\mathrm{\&kappa;}}_m\&CenterDot;|S_1^m|/A^{2m-1})/2---\left(8\right)$

In above formula WithIt is the same phasor after eliminating phase error or anti-phase figure.Additionally, based on water as S_wWith fat as S_fThe meansigma methods of respective pixel or voxel (With) calculate fatty content and distribution, i.e.Qualitative assessment for this kind of disease of fatty liver and lipoma.

Three, embodiment

Embodiment 1

0.35T medical magnetic resonance imaging instrument is being edited imaging sequence shown in Fig. 6, the area arranging first frequency encoding gradient is the twice that gradient area is read in preparation, the interval (i.e. first echo time TE) setting radio-frequency drive soft pulse center and first chunk selection gradient center is minima, such as TE=6ms, set the Δ t width slightly larger than frequency encoding gradient, such as Δ t=2ms, the polarity arranging second frequency encoding gradient is identical with the polarity of first frequency encoding gradient, the two area equation, Δ τ=1/ Δ f/2=9.8ms is set according to the difference in resonance frequencies Δ f=51Hz of water and fat proton, and slice selective gradient and phase encoding gradient are set with reference to regular spin echo sequence, then saving sequence file.Sequential parameter table arranges acquisition matrix and is sized to 256 × 192 × 2 × 16, sequence repetition time TR=400ms, visual field FOV=250mm, block thickness THK=32mm, accumulative frequency NEX=1, other parameter debugging radio frequency power are set in the usual way to limit 180 ° of hard pulse widths within 100 μ s, then preserve parameter list file.

Then, arranging radio-frequency pulse mid frequency is water proton resonant frequency, in the following manner the reading gradient in frequency coding direction is carried out eddy current compensation on the spectrometer with preemphasis function:

Eddy current cycle tests shown in Fig. 5 arranges each reading gradient reading the many echo sequences of gradient shown in the amplitude of gradient pulse, width and polarity and Fig. 2 or Fig. 4Identical, run this sequence and debug vortex field timeconstantτ in following formula_nAnd amplitude alpha_nUntil free induction decay signal (FID) is completely the same when closing with gradient；

The reading gradient waveform of many gradin-echos is compensated respectively, namely according to gradient pre-emphasis waveform

$G_x^{\&prime;i}\left(t\right)=G_x^i\left(t\right)+\left(dG\right(t)/dt)\&CircleTimes;\underset{n}{\&Sigma;}{\mathrm{\&alpha;}}_n{e}^{-t/{\mathrm{\&tau;}}_n}---\left(9\right)$

Then, with cycle tests shown in Fig. 3 selected water fat imaging region collecting magnetic resonance signal and calculate in the manner aforesaidAgain with cycle tests shown in Fig. 4 selected water or fat rich region collecting magnetic resonance signal and calculate in the manner aforesaidWithThe chemical shift imaging sequence shown in service chart 6 when phase encoding gradient closedown, the phase place reunion gradient of debugging chunk choice direction is until echo amplitude is maximum, running pipeline displacement study imaging sequence gather four-dimensional complex matrix (256 × 192 × 2 × 16) under the normally-open condition of phase encoding gradient again, carry out one-dimensional discrete Fourier transformation and obtain 16 three-dimensional complex matrixs (256 × 192 × 2), extract one group of size successively and be the k-space matrix I of 256 × 192₀And I₁And carry out two dimension discrete fourier transform and obtain anti-phase figure matrix S₀With same phasor matrix S₁As follows:

$S_0=(S_w+S_f)e^{-TE/T_2}{e}^{-{\mathrm{i\&phi;}}_0}---\left(10\right)$

$S_1=(S_w-S_f)e^{-(TE+\mathrm{\&Delta;}t)/T_2-\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;})}---\left(11\right)$

T with reflection field inhomogeneities effect₂ ^*Relaxation process is compared, the T of water fat self₂The signal attenuation that relaxation process causes at short notice is generally less, and the difference between the transverse relaxation behavior of water fat can be ignored, therefore the T in above formula₂Can the signal attenuation that causes of the transverse relaxation behavior of compensation water fat self to a certain extent.

Due to Δ τ=1/ Δ f/2, formula (10) and (11) can be written as:

$S_0=(S_w+S_f)\&CenterDot;A_1\&CenterDot;e^{-{\mathrm{i\&phi;}}_0}---\left(14\right)$

$S_1=(S_w-S_f)\&CenterDot;A_2\&CenterDot;e^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;})}---\left(15\right)$

Here, we defineWithThe wherein T of imaging region₂Value can the sequence according to Fig. 3 be measured, T₂ ^*Value can the sequence according to Fig. 4 be measured.

To same phasor S₀Delivery, namely

|S₀|=(S_w+S_f)A₁(16)

Meanwhile, to same phasor S₀Anti-phase figure S it is multiplied by after taking complex conjugate₁And divided by S₀Mould, namely

$S_{1A}=S_1{S}_0^{*}/\left|S_0\right|=(S_w-S_f)A_2{e}^{-i\mathrm{\&phi;}}---\left(17\right)$

Then, we to obtain phase as follows:

φ=atan2 [Im (S_1A)/Re(S_1A)](18)

Often it is wound around beyond (-π ,+π) scope generation phase place due to φ and causes that mistake occurs in phasing, so adopting Such phase solution to twine algorithm carry out the PHASE DISTRIBUTION figure φ ` of the phase unwrapping each scanning slice of acquisition.Then, based on following formula to S_1ACarry out phase error correction:

S_1B=S_1Ae^iφ`=(S_w-S_f)A₂(19)

Finally, we obtain the water picture being sufficiently separated and fat as follows:

S_w=(| S₀|/A₁+κ·|S₁|/A₂)/2(20)

S_f=(| S₀|/A₁-κ·|S₁|/A₂)/2(21)

Here, correction factor κ=Re (S_1B)/|S_1B|, κ can in (-1 ,+1) scope consecutive variations.

Embodiment 2

nullImaging sequence shown in editor Fig. 5 on 1.5T medical magnetic resonance imaging instrument，The area arranging first frequency encoding gradient is the twice that gradient area is read in preparation，The interval (i.e. first echo time TE) setting radio-frequency drive soft pulse center and first chunk selection gradient center is minima，Such as TE=3ms，Set the Δ t width slightly larger than frequency encoding gradient，Such as Δ t=1.6ms，The polarity arranging second frequency encoding gradient is identical with the polarity of first frequency encoding gradient，The two area equation，By the same manner, the 3rd and the 4th frequency encoding gradient parameter are set，In sequence, Δ τ=1/ Δ f/2=2.2ms is set according to the chemical shift difference Δ f of water and fat proton，And slice selective gradient and phase encoding gradient are set with reference to regular spin echo sequence，Then saving sequence file.Sequential parameter table arranges acquisition matrix and is sized to 256 × 192 × 4 × 16, sequence repetition time TR=400ms, visual field FOV=250mm, block thickness THK=32mm, accumulative frequency NEX=1, other parameter debugging radio frequency power are set in the usual way to limit 180 ° of hard pulse widths within 100 μ s, then preserve parameter list file.

Then, the running pipeline displacement study imaging sequence when phase encoding gradient closedown, debugging chunk selects gradient until echo amplitude is maximum, running pipeline displacement study imaging sequence gather four-dimensional complex matrix (256 × 192 × 4 × 16) under the normally-open condition of phase encoding gradient again, and carry out one-dimensional Fourier transform 16 three-dimensional complex matrixs (256 × 192 × 4) of acquisition, the anti-phase figure matrix S of cumulative acquisition after extracting two groups of two-dimensional complex number matrixes successively and carrying out two dimension discrete fourier transform₀With same phasor matrix S₁Describe with following formula:

$S_0=\underset{m=1}{\overset{2}{\&Sigma;}}{S}_0^m=\underset{m=1}{\overset{2}{\&Sigma;}}(S_w+S_f)e^{-\left(\right(2m-2)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;})/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-2))}---\left(22\right)$

$S_1=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_1^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w-S_f)e^{-\left(\right(2m-1)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;})/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-1))}---\left(23\right)$

Consider that magnetic field homogeneity is insufficient to ideal on home equipment, in above formulaXiang HeItem is ignored；Further, in above formulaAll echoes are all identical, can be included in initial magnetization vector to simplify calculating.

Here, we defineTherefore formula (22) and formula (23) can be reduced to respectively:

$S_0^1=(S_w+S_f)\&CenterDot;e^{-{\mathrm{i\&phi;}}_0}---\left(24\right)$

$S_1^1=(S_w-S_f)\&CenterDot;A\&CenterDot;e^{-i({\mathrm{\&iota;}}_0+\mathrm{\&phi;})}---\left(25\right)$

$S_0^2=(S_w+S_f)\&CenterDot;A^2\&CenterDot;e^{-i({\mathrm{\&phi;}}_0+2\mathrm{\&phi;})}---\left(26\right)$

$S_1^2=(S_w-S_f)\&CenterDot;A^3\&CenterDot;e^{-i({\mathrm{\&phi;}}_0+3\mathrm{\&phi;})}---\left(27\right)$

Obviously, can try to achieve based on formula (24) and formula (25)ButMain by Δ B₀Rather than when being determined by local magnetic susceptibility effect, we can preferentially from complex matrixWithFourier transform pairWithK-space centrage (the k that middle selection signal to noise ratio is the highest_y=0) and delivery, namelyWithThus calculate and obtainTo reduce calculating error to greatest extent.

RightDelivery, namely

$\left|S_0^1\right|=S_w+S_f---\left(28\right)$

Thus obtainIt is multiplied by formula (24) to (27) again and eliminates initial phase, therefore the complex matrix stated of formula (24) to (27) is rewritable is:

$S_{0A}^1=S_w+S_f---\left(29\right)$

$S_{1A}^1=(S_w-S_f)\&CenterDot;\sqrt{{\left|I_0^2\right|}_{k_y=0}/{\left|I_0^1\right|}_{k_y=0}}\&CenterDot;e^{-\mathrm{i\&phi;}}---\left(30\right)$

$S_{0A}^2=(S_w+S_f)\&CenterDot;|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}\&CenterDot;e^{-i\&CenterDot;2\mathrm{\&phi;}}---\left(31\right)$

$S_{1A}^2=(S_w-S_f)\&CenterDot;{\left(\right|I_0^2{|}_{k_y=0}/\left|I_0^1{|}_{k_y=0}\right)}^{3/2}\&CenterDot;e^{-i\&CenterDot;3\mathrm{\&phi;}}---\left(32\right)$

Can calculate by formula (31) that to obtain phase error phi as follows:

$\mathrm{\&phi;}=atan2\&lsqb;\mathrm{Im}\left(S_{0A}^2\right)/\mathrm{Re}\left(S_{0A}^2\right)\&rsqb;/2---\left(33\right)$

Often it is wound around beyond (-π ,+π) scope generation phase place due to φ and causes that mistake occurs in phasing, so adopting Such phase solution to twine algorithm φ matrix carries out the PHASE DISTRIBUTION figure φ of the phase unwrapping each scanning slice of acquisition^{`</sup>.Then, to complex matrixWithIt is multiplied by e respectively<sup>iφ`}、e^{i·3φ`</sup>And e<sup>i·3φ`}Carry out phasing, then obtain following complex matrix:

$S_{1B}^1=(S_w-S_f)\&CenterDot;\sqrt{|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}}---\left(34\right)$

$S_{0B}^2=(S_w+S_f)\&CenterDot;|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}---\left(35\right)$

$S_{1B}^2=(S_w-S_f)\&CenterDot;{\left(\right|I_0^2{|}_{k_y=0}/\left|I_0^1{|}_{k_y=0}\right)}^{3/2}---\left(36\right)$

Then, these four groups obtain S with phasor and anti-phase figure are cumulative respectively₀And S₁And substitute into formula (5) and formula (6) respectively and obtain water fat separate picture, or be directly based upon following formula and produce two groups of water and fat image:

$S_w^1=\left(\right|S_{0A}^1|+{\mathrm{\&kappa;}}_1\&CenterDot;|S_{1B}^1/\left|\sqrt{|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}}\right)/2---\left(37\right)$

$S_f^1=\left(\right|S_{0A}^1|-{\mathrm{\&kappa;}}_1\&CenterDot;|S_{1B}^1|/\sqrt{|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}})/2---\left(38\right)$

$S_w^2=\left(\right|S_{0B}^2|/(|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0})+{\mathrm{\&kappa;}}_2\&CenterDot;|S_{1B}^2|/{\left(|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}\right)}^{3/2})/2---\left(39\right)$

$S_f^2=\left(\right|S_{0B}^2|/(|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0})-{\mathrm{\&kappa;}}_2\&CenterDot;|S_{1B}^2|/{\left(|I_0^2{|}_{k_y=0}/|I_0^1{|}_{k_y=0}\right)}^{3/2})/2---\left(40\right)$

Here, correction factorThe two consecutive variations in (-1 ,+1) scope.Finally, water and fat image carry out the cumulative water fat separate picture obtaining signal to noise ratio enhancing respectively.

Embodiment 3

nullImaging sequence shown in editor Fig. 8 (or Fig. 9) on 1.5T medical magnetic resonance imaging instrument，The area arranging first frequency encoding gradient is the twice that gradient integral area is read in preparation，The interval (i.e. first echo time TE) setting radio-frequency drive soft pulse center and first chunk selection gradient center is minima，Such as TE=4ms，Set the Δ t width slightly larger than frequency encoding gradient，Such as Δ t=2ms，The polarity arranging second frequency encoding gradient is identical with the polarity of first frequency encoding gradient，The two area equation，By the same manner, the 3rd and the 4th frequency encoding gradient parameter are set，In sequence, Δ τ=1/ Δ f/2=2.2ms is set according to the chemical shift difference Δ f of water and fat proton，And slice selective gradient and phase encoding gradient are set with reference to regular spin echo sequence.For sequence shown in Fig. 9, replace the frequency encoding gradient of 180 ° of reunion pulses to require opposite polarity, amplitude significantly larger when keeping integral area consistent with normal frequency encoding gradient and Δ τ=1/ Δ f/2 decays with echo-signal during reducing Δ τ.Sequential parameter table arranges acquisition matrix and is sized to 256 × 192 × 4, sequence repetition time TR=400ms, visual field FOV=250mm, thickness THK=3mm, accumulative frequency NEX=1, other parameter it is set in the usual way and limits 180 ° of soft pulse width within 2ms, then preserving parameter list file.

Then, the running pipeline displacement study imaging sequence when phase encoding gradient closedown, debugging slice selective gradient is until echo amplitude is maximum, running pipeline displacement study imaging sequence gather three-dimensional complex matrix (256 × 192 × 4) under the normally-open condition of phase encoding gradient again, after extracting two groups of two-dimensional complex number matrixes successively and carrying out two dimension discrete fourier transform, accumulation obtains anti-phase figure matrix S₀With same phasor matrix S₁, Theoretical Calculation and data processing method according still further to example two obtain the water fat separate picture that signal to noise ratio strengthens.

Claims (3)

1. a high-precision three-dimensional chemical shift imaging method, it is characterised in that comprise the following steps:

The first step: first 180 ° of reunion pulse is applied to the slice selective gradient centre position to first echo, i.e. TE/2 place, after phase encoding gradient is all applied to first 180 ° of reunion pulse, a series of 180 ° of reunion pulses subsequently are applied to Δ t/2 place, Δ t is set to minima and Δ τ=1/ Δ f/2；Selecting gradient to be applied behind an opposite polarity phase place reunion gradient to avoid chunk to select gradient to cause magnetization vector phase dispersion at chunk, each sequence repeat period end applies the damage phase gradient of an amplitude change at random；In sequential parameter table, arranging TE is minima, and to arrange sampled data points be Dim1, and phase code step number is Dim2, and selecting layer direction phase code step number is Dim3；Above-mentioned imaging sequence module is run, successively to G in single-shot mode_p1And G_P2Amplitude carries out phase code circulation, and often step circle collection N group homophase and anti-phase echo are until all phase code step numbers complete, and thus builds four-dimensional complex matrix (Dim1 × Dim2 × N × Dim3)；Wherein, G_p1For selecting layer direction phase encoding gradient, G_p2For phase encoding gradient, Δ t is time delay, and Δ f is water fat Chemical shift differences, and TE is the echo time；

Second step: to four-dimensional complex matrix (Dim1 × Dim2 × N × Dim3) along selecting layer direction to carry out the three-dimensional complex matrix (Dim1 × Dim2 × N) of one-dimensional discrete Fourier transformation acquisition Dim3, extract N group two-dimensional complex number matrix (Dim1 × Dim2) more successively and carry out two dimension discrete fourier transform, it is thus achieved that N group homophase and anti-phase image；For the signal that two dimension version imaging sequence gathers, Dim3=1 is set；

3rd step: by following formula to same phasor complex matrix S₀With anti-phase figure complex matrix S₁Carry out data process:

$S_0=\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}{S}_0^m=\underset{m=1}{\overset{N}{\mathrm{\&Sigma;}}}(S_w{e}^{-\&lsqb;TE+(2m-2)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^w}+S_f{e}^{-i\&lsqb;2\mathrm{\&pi;}\mathrm{\&Delta;}f\mathrm{\&Delta;}\mathrm{\&tau;}\&CenterDot;(2m-2))}{e}^{-\&lsqb;TE+(2m-2)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^f})e^{-(2m-2)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-2))}---\left(1\right)$

$S_1=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_1^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w{e}^{-\&lsqb;TE+(2m-1)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^w}+S_f{e}^{-i(2\mathrm{\&pi;}\mathrm{\&Delta;}f\mathrm{\&Delta;}\mathrm{\&tau;}\&CenterDot;(2m-1))}{e}^{-\&lsqb;TE+(2m-1)\&CenterDot;\mathrm{\&Delta;}t\&rsqb;/T_2^f})e^{-(2m-1)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-1))}---\left(2\right)$

Wherein S_wAnd S_fIt is expressed as water and fat components in region, water fat homophase that N gathers after representing disposable radio-frequency drive and the group number of anti-phase echo, and N can strengthen according to signal noise ratio (snr) of image and weighting scheme need to take various different natural number,WithThe respectively T2 constant of the tissue of human body enrichment water and fat,It is and Magnetic field inhomogeneity degree Δ B₀Relevant apparent T2 constant, Δ τ is the chemical potential transposition evolution time, φ₀Being the initial phase of proton magnetization vector, φ is the phase error that the effect such as magnetic field bump and body area's susceptibility produces during a Δ τ；

When magnetic field homogeneity is insufficient to ideal, in above formulaWithItem can be ignored, therefore above formula can be reduced to when Δ τ=1/ Δ f/2:

$S_0=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_0^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w+S_f)e^{-(2m-2)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-2))}---\left(3\right)$

$S_1=\underset{m=1}{\overset{N}{\&Sigma;}}{S}_1^m=\underset{m=1}{\overset{N}{\&Sigma;}}(S_w-S_f)e^{-(2m-1)\&CenterDot;\mathrm{\&Delta;}\mathrm{\&tau;}/T_2^{*}}{e}^{-i({\mathrm{\&phi;}}_0+\mathrm{\&phi;}\&CenterDot;(2m-1))}---\left(4\right)$

The free induction decay signal gathered by nonlinear fitting T2 cycle tests or a succession of water fat echo-signal, to single exponential function, are obtained respectivelyWithApproximation, and obtain the same phasor complex matrix after echo amplitude correction and anti-phase figure complex matrix based on formula (1)-(2) or formula (3)-(4)；

4th step: same phasor complex matrix and anti-phase figure complex matrix are carried out phase unwrapping and phasing, described in comprising the following steps that:

By rightIt is multiplied by after taking complex conjugateAnd divided byMould, it is thus achieved that there is no start-phase φ₀The same phasor complex matrix of impact | S₀| withAnd anti-phase figure complex matrix

Based onCalculate phase placeHere atan2 () is four-quadrant arctan function, or preferential based on same phasorCalculate phase placeAnd adopt conventional Branch cut (branchcut) or region growth method (regiongrowth) φ to carry out phase unwrapping to obtain real phasing matrix φ `；

RightIt is multiplied by e^{iφ`</sup>,It is multiplied by e<sup>i·2φ`}WithIt is multiplied by e^i·3φ`Eliminate phase error respectively；

At N > 2, other each group carries out phasing with all available similar fashion of phasor and anti-phase figure, namelyIt is multiplied by e^{i·(2m-2)φ`</sup>,It is multiplied by e<sup>i·(2m-1)φ`}；

Based on φ `=γ Δ B₀(2m-1) Δ τ obtains field pattern Δ B₀, and the cosine value defining the phasing matrix of the anti-phase figure of m-th is correction factor matrix κ_mFor determining that in anti-phase figure, the pixel of moisture fat signal should belong to water images or fat image；

5th step: water fat signal separation module and image reconstruction module produce, to all same phasor after phasing and anti-phase figure are cumulative respectively, the S that signal to noise ratio strengthens₀And S₁,WithIn the insignificant situation of item, produce water according to the following formula as S_wWith fat as S_f:

S_w=(| S₀|+κ_m·|S₁|/A)/2 (5)

S_f=(| S₀|-κ_m·|S₁|/A)/2(6)

Or calculate the water producing to be sufficiently separated according to the following formula as S_wWith fat as S_f:

$S_w=\underset{m=1}{\overset{N}{\&Sigma;}}\left(\right|S_0^m|/A^{2m-2}+{\mathrm{\&kappa;}}_m\&CenterDot;|S_1^m|/A^{2m-1})/2---\left(7\right)$

$S_f=\underset{m=1}{\overset{N}{\&Sigma;}}\left(\right|S_0^m|/A^{2m-2}-{\mathrm{\&kappa;}}_m\&CenterDot;|S_1^m|/A^{2m-1})/2---\left(8\right)$

In above formula WithIt is the same phasor after eliminating phase error or anti-phase figure.

2. high-precision three-dimensional chemical shift imaging method according to claim 1, it is characterised in that take a step forward in the first step and include pre-scanning process:

At the slice selective gradient of three orthogonal directions and position the test zone of human body under three corresponding soft pulse effects, then excite and gather free induction decay signal with pulse or by a succession of equally spaced 180 ° of pulse excitations and in a succession of water fat echo-signal of centric acquisition in 180 ° of pulse spacings.

3. high-precision three-dimensional chemical shift imaging method according to claim 1, it is characterised in that farther include the process of content and the distribution calculating fat after the 5th step:

Based on water as S_wWith fat as S_fThe meansigma methods of respective pixel or voxelWithCalculate content and the distribution of fat, it may be assumed that