A quantitative susceptibility mapping method and device based on radial sampling

Abstract

The application discloses a kind of quantitative magnetic susceptibility imaging method and device based on radial sampling, comprising: by magnetic resonance scanner executes three-dimensional gradient echo multi-echo sequence, and the three-dimensional k space multi-echo data of tissue is collected in radial sampling mode;The three-dimensional k space multi-echo data collected is corrected and density is compensated to obtain the three-dimensional k space multi-echo data and corrected track RT after track correction and density compensation;The three-dimensional k space multi-echo data and corrected track RT after track correction and density compensation are input into non-uniform fast Fourier transform reconstruction complex image, to obtain corresponding amplitude diagram and phase diagram;Phase diagram is unwound and background field is removed, to obtain corresponding local field;By improved sparse linear equation and least square method ILSQR algorithm inversion local field, corresponding magnetic susceptibility distribution diagram is generated.The method and device of the application reduce the influence of motion artifact, improve imaging quality and diagnostic accuracy.

Description

Technical Field

[0001] This invention relates to the field of magnetic resonance imaging technology, and in particular to a quantitative magnetic susceptibility imaging method and apparatus based on radial sampling. Background Technology

[0002] Quantitative magnetic susceptibility imaging (QSM) is a technique that uses MRI phase information to infer local magnetic susceptibility. It can be used to detect the distribution of specific elements (such as iron) in the human body, thereby aiding in the diagnosis of specific diseases. Existing QSM methods mostly employ Cartesian sampling, reconstructing images after uniform sampling in K-space. However, Cartesian sampling is sensitive to motion, easily producing motion artifacts when imaging areas easily affected by respiratory movements, such as the liver; furthermore, the scan time is relatively long, making it less suitable for special populations who cannot hold their breath for extended periods. Summary of the Invention

[0003] To address the technical problems existing in the background art, the present invention proposes a quantitative magnetic susceptibility imaging method and device based on radial sampling.

[0004] This invention proposes a quantitative magnetic susceptibility imaging method based on radial sampling, comprising the following steps:

[0005] S1. A three-dimensional gradient echo multi-echo sequence is executed using a magnetic resonance scanner to acquire three-dimensional k-space multi-echo data of the tissue in a radial sampling manner.

[0006] S2. Perform eddy current delay correction and density compensation on the acquired three-dimensional k-space multi-echo data to obtain the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT.

[0007] S3. Input the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT into the non-uniform fast Fourier transform to reconstruct the complex image to obtain the corresponding amplitude map and phase map.

[0008] S4. Unwrap the phase diagram and remove the background field to obtain the corresponding local field;

[0009] S5. The local field is inverted using the improved sparse linear equation and the least squares ILSQR algorithm to generate the corresponding magnetic susceptibility distribution map.

[0010] Preferably, step S2 specifically includes:

[0011] The acquired 3D k-space multi-echo data is input into the SAGE algorithm for simultaneous automatic calibration and gradient delay estimation to obtain gradient delay correction parameters.

[0012] The trajectory T during the reconstruction process is obtained, and the trajectory T is corrected according to the gradient delay correction parameter to obtain the corrected trajectory RT;

[0013] The corrected trajectory RT is input into the Voronoi algorithm of the Thiessen polygon to obtain the corresponding density compensation weight W. Based on the density compensation weight W, the collected three-dimensional k-space multi-echo data is corrected to obtain the three-dimensional k-space multi-echo data after trajectory correction and density compensation.

[0014] Preferably, step S4 specifically includes:

[0015] The phase diagram is processed using the Laplace unwrapping algorithm to obtain the corresponding unwrapped phase diagram;

[0016] Divide the unwrapped phase diagram by the echo time TE to obtain the corresponding field diagram;

[0017] After removing the background field from the field map, it is input into the V-SHARP algorithm for complex harmonic artifact suppression with variable convolution kernel radius to obtain the local field.

[0018] Preferably, the execution process of the SAGE algorithm for simultaneous automatic calibration and gradient delay estimation in step S2 includes:

[0019] Construct the objective function:

[0020]

[0021] in, Indicates the number of sampling repetitions; Indicates the number of coils used in the sampling; It is a low-resolution image that requires self-calibration, acquired by coil i; It is the k-space center data used for calibration, where i represents the i-th coil and j represents the j-th repeated sample; This indicates a non-uniform Fast Fourier Transform operator, which is subject to delay. The resulting k-space shift effect; F is the Fourier transform operator, H is a multidimensional linear operator used to construct a structured Hankel block correction matrix, i.e., correction matrix A, from the relevant Cartesian k-space data; y is the k-space data obtained in magnetic resonance sampling. and It is an unknown quantity;

[0022] Alternately optimize the coil image x and the gradient delay correction parameters, and output the converged gradient delay correction parameters.

[0023] Preferably, the V-SHARP algorithm for complex harmonic artifact suppression with variable convolution kernel radius uses a variable convolution kernel radius, employing a small-diameter spherical function near the edge of the imaging region and a large-diameter spherical function in the central region to remove the background field and preserve the local field.

[0024] Preferably, the radial sampling is based on a star-studded sampling model, employing a radial sampling trajectory in the Kx-Ky plane, with each spoke passing through the center of the K space, and uniform sampling in the Kz direction, with the sampling direction controlled by adjusting the rotation matrix.

[0025] Preferably, in step S1, the radial sampling is controlled by a rotation matrix to control the gradient direction. The rotation matrix is ​​dynamically generated based on the list of spoke rotation angles. Specifically, the rotation matrix is:

[0026] ;

[0027] in, Let be the rotation angle of the j-th spoke, and let the rotation matrix be based on the initial matrix. By spoke rotation angle Calculated.

[0028] Preferably, the tissue is one of brain tissue, liver tissue, muscle tissue, adipose tissue, or other soft tissue.

[0029] This invention proposes a quantitative magnetic susceptibility imaging device based on radial sampling, comprising:

[0030] The data acquisition module is used to acquire three-dimensional k-space multi-echo data of tissues by performing a three-dimensional gradient echo multi-echo sequence through a magnetic resonance scanner in a radial sampling manner.

[0031] The correction module is used to perform eddy current delay correction and density compensation on the acquired three-dimensional k-space multi-echo data to obtain the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT.

[0032] The reconstruction module is used to input the trajectory-corrected and density-compensated three-dimensional k-space multi-echo data and the corrected trajectory RT into a non-uniform fast Fourier transform to reconstruct a complex image, so as to obtain the corresponding amplitude map and phase map.

[0033] The field processing module is used to unwrap the phase map and remove the background field to obtain the corresponding local field.

[0034] The QSM reconstruction module is used to invert the local field using an improved sparse linear equation and the least squares ILSQR algorithm, and generate the corresponding magnetization distribution map.

[0035] This invention proposes a radial sampling-based quantitative magnetic susceptibility imaging method and apparatus that combines radial sampling with quantitative magnetic susceptibility imaging (QSM) reconstruction. It performs three-dimensional K-space radial sampling using a three-dimensional gradient echo multi-echo sequence. After data acquisition, the data undergoes delay correction, density compensation, image reconstruction, phase unwrapping, background field removal, and magnetic susceptibility inversion to obtain a magnetic susceptibility map. This invention employs radial sampling, reducing sensitivity to motion, minimizing respiratory motion artifacts in areas such as the liver, and improving image quality. It is particularly suitable for imaging areas susceptible to respiratory motion interference, such as the liver, while also shortening scan time, ensuring the accuracy of data correction and reconstruction, and improving the quantitative precision of magnetic susceptibility imaging. Attached Figure Description

[0036] Figure 1 This is a schematic diagram of the sampling trajectory structure of the radial sampling in the quantitative magnetic susceptibility imaging method based on radial sampling proposed in this invention.

[0037] Figure 2 This is a schematic diagram illustrating the workflow of a quantitative magnetic susceptibility imaging method based on radial sampling proposed in this invention.

[0038] Figure 3 This is a timing diagram of radial sampling in a quantitative magnetic susceptibility imaging method based on radial sampling proposed in this invention.

[0039] Figure 4 This is a Voronoi diagram of a quantitative magnetic susceptibility imaging method based on radial sampling proposed in this invention.

[0040] Figure 5 This is a flowchart of the SAGE algorithm for a quantitative magnetic susceptibility imaging method based on radial sampling proposed in this invention. Detailed Implementation

[0041] Reference Figures 1-5 The present invention proposes a quantitative magnetic susceptibility imaging method based on radial sampling, comprising the following steps:

[0042] S1. A three-dimensional gradient echo multi-echo sequence is executed using a magnetic resonance scanner to acquire three-dimensional k-space multi-echo data of the tissue in a radial sampling manner.

[0043] In this embodiment, radial sampling is based on the star-studded sampling model. A radial sampling trajectory is used in the Kx-Ky plane, with each spoke passing through the center of the K space. Uniform sampling is used in the Kz direction, and the sampling direction is controlled by adjusting the rotation matrix.

[0044] Specifically, based on the stack-of-stars magnetic resonance sampling model, radial sampling is performed using a three-dimensional gradient echo multi-echo sequence. The stack-of-stars sampling model is a sampling method for three-dimensional radial sampling, such as... Figure 1 As shown.

[0045] During 3D K-space sampling, sampling points are distributed on parallel sampling surfaces in K-space, and the same radial sampling structure is used on each sampling surface. That is, radial sampling is used on the Kx-Ky plane, while uniform sampling is used in the Kz direction. When implementing a 3D gradient echo multi-echo sequence, multi-echo sampling is achieved through regression. Regression means that after each echo acquisition, a dephasing gradient waveform is used to return the frequency direction in K-space to the initial sampling position before acquiring the next echo. In the radial sampling process, the direction of frequency encoding is modified and acquired each time by adjusting the gradient magnitudes in the output x and y directions. Encoding in the selected layer direction is achieved through phase encoding.

[0046] In this embodiment, a Fitt spectrometer is used to realize a three-dimensional gradient echo multi-echo sequence, such as... Figure 3 As shown, the three-dimensional gradient echo multi-echo sequence refers to the timing diagram applied by modules such as gradient waveform, radio frequency pulse, and receiver during the sampling process.

[0047] The specific implementation process is as follows:

[0048] In sequence design, the required field of view (FOV) for imaging is generally determined first based on the specific scanning conditions. Then, considering factors such as sampling time, signal-to-noise ratio, and resolution, parameters such as the number of spokes (noRadials) and the number of sampling points (noSamples) on each spoke are determined. The coordinates of the sampling points, i.e., the sampling trajectory, can be calculated based on noRadials, noSamples, and FOV. The x and y coordinates of the sampling points on each spoke are located within [-]. ]and[- A uniform distribution between ] . Among them, This is the rotation angle corresponding to each spoke. The set of sampling point coordinates can be called the sampling trajectory T. Then, considering the trade-off between signal-to-noise ratio and sampling time, the sampling period samplePeriod is used to determine the gradient magnitude Amp in the sequence waveform. Finally, the application time t of different waveforms and pulses is determined to avoid aliasing, thus initially determining the timing diagram S0 of the sequence.

[0049] Specifically, the effective field of view (FOV) of the reconstructed image is subject to the following limitations in terms of resolution and minimum number of spokes:

[0050]

[0051] in, This represents the sampling interval for each sampling point in k-space; This indicates the sampling interval between sampling points on each spoke; This indicates the sampling interval between sampling points between two spokes; This indicates the difference between the sampling angles of the spokes; Represents the maximum amplitude value of the coordinates of the sampling point in k-space; res represents the resolution in the imaging space; FOV represents the length of the effective (non-aliasing) FOV.

[0052] Common sampling strategies include:

[0053]

[0054] Radial in or centripetal radial sampling refers to K-space data being filled from the edge to the center, similar to scanning "from the outside in". Radial sampling, also known as centrifugal radial sampling, refers to filling K-space data from the center outwards, similar to a "from the inside out" scan. , These correspond to the rotation angles of the spokes when uniformly sampling the semicircle and the entire circle, respectively.

[0055] The mutual constraints between some imaging parameters are determined based on the Nyquist sampling theorem.

[0056] In this embodiment, in step S1, radial sampling controls the gradient direction through a rotation matrix. The rotation matrix is ​​dynamically generated based on the list of spoke rotation angles. Specifically, the rotation matrix is ​​as follows:

[0057] ;

[0058] in, Let be the rotation angle of the j-th spoke, and let the rotation matrix be based on the initial matrix. By spoke rotation angle The calculation shows that the z-axis is used as the layer selection axis at this time.

[0059] Specifically, around vectors The rotation matrix is ​​as follows:

[0060]

[0061] In actual sampling, rotation around the x, y, and z axes is used.

[0062] In this embodiment, the list of angles AL required for sampling is first determined, and then the table ML of the rotation matrix is ​​calculated based on AL using the formula above. The index of the table ML is used to control the sampling direction during the sampling process.

[0063] In this embodiment, the tissue is specifically one of brain tissue, liver tissue, muscle tissue, adipose tissue, or other soft tissue.

[0064] In this embodiment, a radial sampling mode is used for sampling, which can be applied to imaging various body parts. The radial sampling mode can reduce motion artifacts when imaging areas with motion artifacts, such as the liver. After obtaining the image domain data, the magnetic susceptibility distribution is obtained by combining it with the QSM reconstruction algorithm.

[0065] S2. Perform eddy current delay correction and density compensation on the acquired three-dimensional k-space multi-echo data to obtain the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT.

[0066] In this embodiment, step S2 specifically includes:

[0067] The acquired 3D k-space multi-echo data is input into the Simultaneous Auto-calibration and Gradient Delays Estimation (SAGE) algorithm to obtain gradient delay correction parameters.

[0068] The trajectory T during the reconstruction process is obtained, and the trajectory T is corrected according to the gradient delay correction parameter to obtain the corrected trajectory RT;

[0069] The corrected trajectory RT is input into the Voronoi algorithm of the Thiessen polygon to obtain the corresponding density compensation weight W. Based on the density compensation weight W, the collected three-dimensional k-space multi-echo data is corrected to obtain the three-dimensional k-space multi-echo data after trajectory correction and density compensation.

[0070] In this embodiment, the execution process of the SAGE algorithm for simultaneous automatic calibration and gradient delay estimation in step S2 includes:

[0071] Construct the objective function:

[0072]

[0073] in, Indicates the number of sampling repetitions; Indicates the number of coils used in the sampling; It is a low-resolution image that requires self-calibration, acquired by coil i; It is the k-space center data used for calibration, where i represents the i-th coil and j represents the j-th repeated sample; This indicates a non-uniform Fast Fourier Transform operator, which is subject to delay. The resulting k-space shift effect; F is the Fourier transform operator, H is a multidimensional linear operator used to construct a structured Hankel block correction matrix, i.e., correction matrix A, from the relevant Cartesian k-space data; y is the k-space data obtained in magnetic resonance sampling. and It is an unknown quantity;

[0074] Alternately optimize the coil image x and the gradient delay correction parameters, and output the converged gradient delay correction parameters.

[0075] Specifically, x and The solution is as follows:

[0076] First, ignore it. The influence of the pseudo-inverse matrix is ​​used to calculate A:

[0077]

[0078] Where A is the correction matrix, which satisfies The approximate solution is H, which is a multidimensional linear operator used to construct A. yes The pseudo-inverse matrix, It is the collected k-space data. For equations An approximate solution for x.

[0079] Then perform singular value decomposition on A:

[0080]

[0081] in, and It is a unitary matrix used for singular value decomposition. Indicates transpose. Describes the singular value matrix of A, with diagonal elements from or The square roots of the eigenvalues ​​are sorted in descending order.

[0082] Then, a hard threshold is applied to the singular values ​​to obtain the corrected values. :

[0083]

[0084] Then, by combining the pseudo-inverse matrix, we can obtain an approximate x:

[0085]

[0086] in, This represents the pseudo-inverse matrix.

[0087] After solving for x, an approximate estimate is then made using the Gauss-Newton method. For the formula:

[0088]

[0089] First of all middle The effect is approximately first-order, that is:

[0090]

[0091] The reciprocal can be calculated using the chain rule, and by combining this with the frequency shift property of the Fourier transform, we can obtain:

[0092]

[0093] Where L is a ramp filter. , The partial derivative of the delay is related to the sampling trajectory. In radial sampling, it corresponds to the cosine of the sampling spoke. sin .

[0094]

[0095] This is how we obtain it. It can be used for updates later. , Get new The approximate value is then used for the next iteration. This process continues until a certain number of iterations are reached or an approximate value is obtained. The iteration can stop when the amplitude value is less than a certain size.

[0096] S3. Input the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT into the non-uniform fast Fourier transform to reconstruct the complex image to obtain the corresponding amplitude map and phase map.

[0097] Specifically, during radial sampling, delays caused by eddies and other system issues may occur, resulting in unintended data trajectories. This invention employs the Simultaneous Automatic Calibration and Gradient Delay Estimation (SAGE) algorithm to correct these delays. This method utilizes the low-rank property of multi-channel acquired images to estimate and correct the effects of eddies. By inputting the acquired 3D k-space multi-echo data into the SAGE algorithm, gradient delay correction parameters can be obtained. These parameters can then be used to correct the trajectory T used in the reconstruction process, yielding RT for subsequent reconstruction. Regarding density compensation, since radial sampling is non-uniform, density compensation is necessary. This invention uses a Voronoi diagram to calculate density compensation weights, such as... Figure 4As shown, the corrected trajectory RT is input into the Voronoi algorithm of the Thiessen polygon to obtain the weight W. W is then used to correct the acquired 3D k-space multi-echo data, resulting in 3D k-space multi-echo data with trajectory correction and density compensation. Mesh reconstruction is then performed using the NUFFT algorithm. The 3D k-space multi-echo data with trajectory correction and density compensation, along with the trajectory RT, are input to obtain the reconstructed image I. Since it is a complex domain image, it includes an amplitude map IM and a phase map IP, which are subsequently processed for QSM reconstruction.

[0098] The obtained phase map IP is then de-wound. In this invention, the Laplace unwinding algorithm is used. The following formula illustrates the phase winding situation:

[0099]

[0100] in, Represents a local field; The Larmor precession frequency represents the influence of the magnetic field on the resonant frequency; TE is the echo time. The range is , Represents the actual phase; the Laplace unwinding is implemented as follows:

[0101]

[0102] That is what I wanted. The derivation is as follows:

[0103] First, construct a function:

[0104] Let's construct another formula: Then Substitute Replacement Later obtained By equipping the imaginary parts on both sides, we can obtain the untangled relation:

[0105] ;

[0106] This is the phase relationship before and after unwinding, which can be used as a basis for... and The relationship between them was obtained .

[0107] S4. Unwrap the phase diagram and remove the background field to obtain the corresponding local field.

[0108] In this embodiment, step S4 specifically includes:

[0109] The phase diagram is processed using the Laplace unwrapping algorithm to obtain the corresponding unwrapped phase diagram;

[0110] Divide the unwrapped phase diagram by the echo time TE to obtain the corresponding field diagram;

[0111] After removing the background field from the field map, it is input into the V-SHARP algorithm for complex harmonic artifact suppression with variable convolution kernel radius to obtain the local field.

[0112] In this embodiment, the V-SHARP (Variable-kernel Sophisticated Harmonic Artifact Reduction for Phase Dat) algorithm employs a variable convolution kernel radius, using a small-diameter spherical function near the edge of the imaging region and a large-diameter spherical function in the central region to remove the background field and preserve the local field.

[0113] Specifically, the background field refers to the magnetic field of the susceptibility source outside the imaging region, while the local field refers to the magnetic field generated by the susceptibility source within the imaging region. The V-SHARP method utilizes Maxwell's equations, because the background magnetic field... The magnetic susceptibility originates from a source outside the region of interest, therefore it satisfies the following internally. , and local field The equation is not satisfied. Therefore... It is a harmonic function. According to the mean-mean theorem for harmonic functions, if the function u is... The harmonic function on, then u in The value of u(x) on the sphere is equal to the value of u on the sphere. The average value over the circle is also equal to the function value at the center. That is, the harmonic function is the same as a non-negative, normalized radially symmetric function (the sum of all elements equals 1). The convolution of the global field is equal to itself. Therefore, the following operation can be performed on the global field:

[0114]

[0115] achievable ,in The unit impulse function, This represents the total magnetic field. It can be obtained using Fourier transform. . This represents the tissue mask matrix, and C is the convolution kernel. The Fourier transform form of the function is given. Considering that a small-diameter spherical function corresponds to a large deconvolution kernel, and a large-diameter spherical function corresponds to a small deconvolution kernel, using a small-diameter spherical function leads to a larger residual phase error. Using a large-diameter spherical function at the tissue edge introduces unwanted external signals, blurring the edge, which will be removed during tissue masking. The larger the radius of the spherical function used, the more edge areas are removed. Therefore, a small-diameter spherical function is used near the edge, and a large-diameter spherical function is used near the center. This is the complex harmonic artifact suppression method with variable convolution kernel radius.

[0116] S5. The local field is inverted using the improved sparse linear equation and least-squares algorithm (ILSQR) to generate the corresponding magnetic susceptibility distribution map.

[0117] Specifically, due to the long scan time of magnetic resonance imaging (MRI), motion artifacts can occur in the final image if there is movement during the scan, affecting observation and diagnosis. QSM reconstruction can obtain the magnetic susceptibility distribution in tissues, and for some areas such as the liver and brain, it can be used to detect and diagnose specific diseases. However, liver scans are easily affected by respiratory motion. Therefore, we can suppress respiratory motion artifacts by using radial sampling, thereby improving the accuracy of the final image reconstruction and diagnostic results.

[0118] Reference Figures 1-5 The present invention proposes a quantitative magnetic susceptibility imaging device based on radial sampling, comprising:

[0119] The data acquisition module is used to acquire three-dimensional k-space multi-echo data of tissues by performing a three-dimensional gradient echo multi-echo sequence through a magnetic resonance scanner in a radial sampling manner.

[0120] The correction module is used to perform eddy current delay correction and density compensation on the acquired three-dimensional k-space multi-echo data to obtain the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT.

[0121] The reconstruction module is used to input the trajectory-corrected and density-compensated three-dimensional k-space multi-echo data and the corrected trajectory RT into a non-uniform fast Fourier transform to reconstruct a complex image, so as to obtain the corresponding amplitude map and phase map.

[0122] The field processing module is used to unwrap the phase map and remove the background field to obtain the corresponding local field.

[0123] The QSM reconstruction module is used to invert the local field using an improved sparse linear equation and the least squares ILSQR algorithm, and generate the corresponding magnetization distribution map.

[0124] The above description is only a preferred embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any equivalent substitutions or modifications made by those skilled in the art within the scope of the technology disclosed in the present invention, based on the technical solution and inventive concept of the present invention, should be covered within the scope of protection of the present invention.

Claims

1. A quantitative magnetic susceptibility imaging method based on radial sampling, characterized by, The method comprises the following steps: S1, performing a three-dimensional gradient echo multi-echo sequence by a magnetic resonance scanner to acquire three-dimensional k-space multi-echo data of a tissue in a radial sampling manner; S2, performing eddy current delay correction and density compensation on the acquired three-dimensional k-space multi-echo data to obtain three-dimensional k-space multi-echo data after trajectory correction and density compensation and corrected trajectory RT; S3, inputting the three-dimensional k-space multi-echo data after trajectory correction and density compensation and the corrected trajectory RT into a non-uniform fast Fourier transform to reconstruct a complex image to obtain corresponding amplitude and phase images; S4, performing unwrapping and background field removal on the phase image to obtain a corresponding local field; S5, inverting the local field by an improved sparse linear equation and least square method ILSQR algorithm to generate a corresponding magnetic susceptibility distribution map; Step S2 specifically comprises: inputting the acquired three-dimensional k-space multi-echo data into a simultaneous auto-calibration and gradient delay estimation SAGE algorithm to obtain gradient delay correction parameters; acquiring a trajectory T in a reconstruction process, and correcting the trajectory T according to the gradient delay correction parameters to obtain a corrected trajectory RT; inputting the corrected trajectory RT into a Voronoi algorithm to obtain corresponding density compensation weights W, and correcting the acquired three-dimensional k-space multi-echo data based on the density compensation weights W to obtain three-dimensional k-space multi-echo data after trajectory correction and density compensation; The execution process of the simultaneous auto-calibration and gradient delay estimation SAGE algorithm in step S2 comprises: Constructing the objective function: where, denotes the number of repetitions of the sampling; denotes the number of coils used in the sampling; is the low-resolution image acquired by coil i that needs to be self-calibrated; is the k-space center data used for calibration, i denotes the i-th coil, and j denotes the j-th repetition of the sampling; denotes the non-uniform fast Fourier transform operator, which is affected by the k-space shift induced by the delay ; F is the Fourier transform operator, and H is a multi-dimensional linear operator used to construct a structured Hankel block correction matrix from the relevant Cartesian k-space data, i.e., the correction matrix A; y is the k-space data obtained in the magnetic resonance sampling, and are unknown quantities; alternately optimizing a coil image x and gradient delay correction parameters, and outputting converged gradient delay correction parameters.

2. The quantitative susceptibility mapping method based on radial sampling of claim 1, wherein, Step S4 specifically comprises: processing the phase image by a Laplacian unwrapping algorithm to obtain a corresponding unwrapped phase image; dividing the unwrapped phase image by an echo time TE to obtain a corresponding field image; inputting the field image after background field removal into a variable convolution kernel radius complex harmonic artifact suppression V-SHARP algorithm to obtain a local field.

3. The quantitative susceptibility mapping method based on radial sampling of claim 2, wherein, The variable convolution kernel radius complex harmonic artifact suppression V-SHARP algorithm uses a variable convolution kernel radius, uses a small-diameter spherical function near an edge of an imaging region, and uses a large-diameter spherical function in a central region to remove a background field and reserve a local field.

4. The quantitative susceptibility mapping method based on radial sampling of claim 1, wherein, The radial sampling is based on a star stack sampling model, adopts a radial sampling trajectory in a Kx-Ky plane, each spoke passes through a K-space center, adopts uniform sampling in a Kz direction, and adjusts a rotation matrix to control a sampling direction.

5. The quantitative susceptibility mapping method based on radial sampling of claim 4, wherein, In step S1, the radial sampling controls the gradient direction by a rotation matrix, which is dynamically generated according to a spoke rotation angle list, and the rotation matrix is specifically as follows: ; wherein, is the rotation angle of the jth spoke, the rotation matrix is based on the initial matrix is calculated by the spoke rotation angle .

6. The quantitative susceptibility mapping method based on radial sampling of claim 1, wherein, The tissue is specifically one of a brain tissue, a liver tissue, a muscle tissue, a fat tissue, or other soft tissues.

7. A quantitative magnetic susceptibility imaging apparatus based on radial sampling, characterized by, The method comprises: a data acquisition module configured to perform a three-dimensional gradient echo multi-echo sequence by a magnetic resonance scanner to acquire three-dimensional k-space multi-echo data of a tissue in a radial sampling manner; a correction module configured to perform eddy current delay correction and density compensation on the acquired three-dimensional k-space multi-echo data to obtain three-dimensional k-space multi-echo data after trajectory correction and density compensation and a corrected trajectory RT; The reconstruction module is configured to input the trajectory-corrected and density-compensated three-dimensional k-space multi-echo data and the corrected trajectory RT into a non-uniform fast Fourier transform to reconstruct a complex image, so as to obtain a corresponding amplitude map and a phase map; The field processing module is configured to unwrap and remove a background field from the phase map, so as to obtain a corresponding local field; The QSM reconstruction module is configured to invert the local field by using an improved sparse linear equation and least square method ILSQR algorithm, so as to generate a corresponding magnetic susceptibility distribution map; The correction module specifically comprises: The acquired three-dimensional k-space multi-echo data is input into a simultaneous automatic calibration and gradient delay estimation SAGE algorithm, so as to obtain gradient delay correction parameters; The trajectory T in the reconstruction process is acquired, and the trajectory T is corrected according to the gradient delay correction parameters, so as to obtain a corrected trajectory RT; The corrected trajectory RT is input into a Voronoi algorithm, so as to obtain corresponding density compensation weights W, and the acquired three-dimensional k-space multi-echo data is corrected based on the density compensation weights W, so as to obtain trajectory-corrected and density-compensated three-dimensional k-space multi-echo data; The execution process of the simultaneous automatic calibration and gradient delay estimation SAGE algorithm in the correction module comprises: Constructing the objective function: where, denotes the number of repetitions of the sampling; denotes the number of coils used in the sampling; is the low-resolution image acquired by coil i that needs to be self-calibrated; is the k-space center data used for calibration, i denotes the i-th coil, and j denotes the j-th repetition of the sampling; denotes the non-uniform fast Fourier transform operator, which is affected by the k-space shift induced by the delay ; F is the Fourier transform operator, and H is a multi-dimensional linear operator used to construct a structured Hankel block correction matrix from the relevant Cartesian k-space data, i.e., the correction matrix A; y is the k-space data obtained in the magnetic resonance sampling, and are unknown quantities; Alternately optimizing the coil image x and the gradient delay correction parameters, and outputting the converged gradient delay correction parameters.

Images