Density-adapted k-space trajectory design method for sodium element magnetic resonance imaging

Abstract

The application discloses a density adaptive k-space trajectory design method for sodium element magnetic resonance imaging. Firstly, a k-space trajectory curve is constructed, the k-space trajectory curve comprises a radial trajectory, and the starting point of the radial trajectory is at the origin of the polar coordinate system of the k-space; corresponding to the radial trajectory, a spiral trajectory is adopted at the periphery of the k-space, and the radial trajectory and the corresponding spiral trajectory are connected through an arc trajectory; then on the k-space trajectory curve, sampling points on the radial trajectory and sampling points on the spiral trajectory are selected as the target to minimize the sampling time, and sampling points on the spiral trajectory are selected as the target to realize uniform sampling density, so that the density adaptive k-space trajectory is obtained. The density adaptive k-space trajectory designed by the application effectively improves the filling efficiency of the k-space trajectory during sampling, reduces the sampling times, simultaneously realizes uniform sampling at the periphery of the k-space, and improves the signal-to-noise ratio of the sodium element magnetic resonance imaging.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of multi-nuclear magnetic resonance imaging (MRI) methods, and particularly relates to a density-adaptive k-space trajectory design method for sodium element magnetic resonance imaging. BACKGROUND

[0002] There are more than 60 elements in the human body, among which many elements play a key role in life activities. Magnetic resonance imaging is a non-invasive and non-ionizing radiation imaging method, which can detect elements with non-zero nuclear spin, including nuclei with spin 1 / 2 (such as 1 H and 31 P), and nuclei with spin greater than 1 / 2 (such as 14 N and 23 Na). Traditional magnetic resonance imaging technology is mainly based on proton ( 1 H) signal, which is used to obtain structural and functional information of tissues. With the development of magnetic resonance technology, multi-nuclear magnetic resonance imaging methods based on non-proton nuclei have gradually emerged, providing rich multi-element information for life activity detection. Among them, sodium ( 23 Na) is the highest endogenous element that can be detected by magnetic resonance in the human body except for 1 H. Its concentration in the body is closely related to ion balance and energy metabolism, and can provide functional metabolic information that 1 H nuclear magnetic resonance imaging cannot obtain. 23 Na nuclear magnetic resonance imaging has become a unique visualization tool for studying ion balance in the body due to its non-invasive detection capability, and has shown broad application prospects in tissue metabolism, pathology and functional research.

[0003] 23 Na is a quadrupolar nucleus, and its T2 relaxation time is short, showing significant double relaxation characteristics in tissues. Under the magnetic field strength of magnetic resonance imaging, 60% of 23 Na signal in the tissue is fast-relaxation signal (T2 relaxation time is in the order of milliseconds), and 40% of 23 Na signal is slow-relaxation signal (T2 relaxation time is in the order of tens of milliseconds). In order to reduce the influence of T2 relaxation on 23 Na signal, 23 Na magnetic resonance imaging mainly uses short echo time (TE) non-Cartesian sampling method, using radial or spiral k-space trajectory. Its typical feature is that the sampling density is larger in the center of the spatial frequency domain (i.e. k-space) and smaller in the periphery, and the non-uniform sampling introduces additional noise. In order to reduce the influence of non-uniform sampling on 23 Na signal,Na signal-to-noise ratio, Konstandin et al. proposed a density adaptive trajectory based on radial k-space trajectory, the main feature of which is to achieve uniform sampling density at the periphery of k-space, thereby improving the signal-to-noise ratio of the image. However, to obtain a sampling matrix of MxM (M is the number of rows or columns of the sampling matrix), the k-space filling efficiency of the density adaptive trajectory based on radial sampling is low, and more than Mxπ samplings are required, and the total sampling time is long. 23 Na magnetic resonance image, the k-space filling efficiency of the density adaptive trajectory based on radial sampling is low, and more than Mxπ samplings are required, and the total sampling time is long. SUMMARY

[0004] In view of the defects of the prior art, the present application provides a density adaptive k-space trajectory design method for sodium element magnetic resonance imaging.

[0005] The above-mentioned object of the present application is realized by the following technical means:

[0006] The density adaptive k-space trajectory design method for sodium element magnetic resonance imaging comprises the following steps:

[0007] Step 1, constructing a k-space trajectory curve, the k-space trajectory curve comprising a radial trajectory, the starting point of the radial trajectory being at the origin of the polar coordinate system of the k-space, corresponding to the radial trajectory, a spiral trajectory is adopted at the periphery of the k-space, the radial trajectory and the corresponding spiral trajectory are connected by a circular arc trajectory, and the two ends of the circular arc trajectory are tangent to the radial trajectory and the spiral trajectory, respectively;

[0008] Step 2, selecting sampling points on the k-space trajectory curve, the trajectory of the sampling points being a density adaptive k-space trajectory, and the k value of the sampling points being the intensity value G i of the encoding gradient;

[0009] The sampling points on the radial trajectory and the sampling points on the circular arc trajectory in the k-space trajectory curve are selected according to the following rules: from the starting time, the intensity value G i of the encoding gradient changes at a climbing rate upper limit value SR max , if the intensity value G i of the encoding gradient reaches the encoding gradient intensity upper limit value G max , the intensity value G i of the encoding gradient remains the encoding gradient intensity upper limit value G max .

[0010] The sampling points on the spiral trajectory in the k-space trajectory curve are selected according to the following rules: the k-space distance of adjacent sampling points remains constant.

[0011] The density-adaptive k-space trajectory design method for sodium element magnetic resonance imaging as described above further comprises step 3: obtaining a complete density-adaptive k-space trajectory of the entire k-space by rotating the density-adaptive k-space trajectory determined in step 2 multiple times.

[0012] Step 1 as described above comprises the following steps:

[0013] Step 1.1, calculating the radial component maximum value of the k value of the sampling point in the k-space:

[0014]

[0015] k max The radial component maximum value is FOV1, which represents the size of one dimension component of the field of view, and Matrix1 represents the size of one dimension component of the sampling matrix;

[0016] Step 1.2, determining the radial component of the k value of the critical point of the switching from the circular arc trajectory to the spiral trajectory as:

[0017]

[0018] Ns represents the preset number of excitations, and k swi represents the radial component of the k value of the critical point of the switching from the circular arc trajectory to the spiral trajectory, which is denoted as the critical point radial component k swi ;

[0019] Step 1.3, setting the spiral trajectory as:

[0020] k(τ)=k max τe jωτ , k swi / k max ≤τ≤1

[0021] Wherein, ω is the angular frequency of the spiral trajectory, ω=k max / k swi ; τ is a time coefficient, k(τ) represents the k value corresponding to the time coefficient τ, and the k value includes the radial component and the angular component;

[0022] Step 1.4, based on the critical point radial component k swi and the spiral trajectory, obtaining the k value of the critical point and the derivative of the spiral trajectory at the critical point; according to the k value of the critical point, the derivative at the critical point and the k-space radius of the circular arc trajectory, obtaining the k value of the center of the circular arc trajectory;

[0023] Based on the fact that the other end of the circular arc trajectory is tangent to the radial trajectory, the k value of the tangent point of the circular arc trajectory and the radial trajectory is obtained according to the k value of the center of the circular arc trajectory and the k-space radius of the circular arc trajectory.

[0024] As described above, step 2 includes the following steps:

[0025] The intensity value G of the encoded gradient is calculated from the start time, at each time interval Δt. i Size based on the upper limit of the climb rate SR max Incrementing, where i represents the index of the gradient intensity value encoded in chronological order; simultaneously, at each time interval Δt, based on the radial and circular trajectories in the k-space trajectory curve, and using the intensity values ​​G of the first N encoded gradients... i Integrating over time yields the corresponding value of k, denoted as k. N Value; k N The radial component of the value is denoted as k. Nr ;

[0026] When the radial component k Nr Radial component k not reaching the critical point swi If the strength value G of the encoded gradient i The size reaches the upper limit of the encoding gradient strength G. max Then the strength value G of the encoded gradient i The size remains the upper limit of the encoded gradient strength G. max The corresponding climb rate SR = 0;

[0027] When the radial component k Nr Reaching or exceeding the critical point radial component k swi At that time, the k-space distance between adjacent sampling points remains a constant preset spatial spacing. Based on the spiral trajectory in the k-space trajectory curve and the preset spatial spacing, sampling points are selected sequentially from the first sampling point on the spiral trajectory until the radial component k corresponding to the sampling point is reached. Nr Reaching or exceeding the maximum value k of the radial component max Stop selecting sampling points.

[0028] As described above, in step 2, the intensity value G of the first N encoded gradients i Integrating over time yields the corresponding value of k, based on the following formula:

[0029]

[0030] N represents the number of intensity values ​​of the encoded gradients participating in the integration. N increases with the length of time Δt, where N∈{1, 2, 3, ..., N′}, and N′ represents the total number of intensity values ​​of the encoded gradients corresponding to the entire density-fit k-space trajectory; γ represents the gyromagnetic ratio of the nuclear spin; G i Let G represent the intensity value of the encoding gradient, denoted as the intensity value of the encoding gradient. i ;k N Let k represent the integral of the strength values ​​of the first N encoded gradients with respect to time. NValue; k N The radial component of the value is denoted as k Nr .

[0031] The k-space distance between adjacent sampling points on the spiral trajectory is constant as 1 / FOV1, and the k value corresponding to the first sampling point on the spiral trajectory is denoted as k n The k-space distance between the first sampling point on the spiral trajectory and the critical point is equal to wherein is the k-space distance between the last sampling point on the circular arc trajectory and the critical point, and k n-1 is the k value corresponding to the last sampling point on the circular arc trajectory, and the last sampling point on the circular arc trajectory is the sampling point before the first sampling point on the spiral trajectory.

[0032] The sampling bandwidth corresponding to the sampling point is as follows:

[0033] SW = max(G i )·γ·FOV1

[0034] SW represents the sampling bandwidth; max(G i ) represents the maximum intensity value in the encoding gradient.

[0035] The step 2 is further comprised of the following steps:

[0036] The intensity values of the respective encoding gradients are calculated:

[0037] G N = (k N -k N-1 ) / (γΔt)

[0038] G N represents the Nth intensity value of the encoding gradient; k N-1 represents the k value corresponding to the integral of the intensity values of the previous N-1 encoding gradients with respect to time.

[0039] The angle of each k-space trajectory rotation in the step 3 is 2π / Ns, and Ns represents the preset number of excitations.

[0040] The k-space radius of the circular arc trajectory is greater than 0 and less than

[0041] Compared with the prior art, the present application has the following beneficial effects:

[0042] 1. The density-adaptive k-space trajectory proposed in the method can efficiently fill the k-space, and the number of sampling points is reduced 23The number of excitation and sampling required by Na magnetic resonance imaging is reduced, and the sampling time is shortened; meanwhile, uniform sampling is realized in the outer periphery of k-space, noise introduced by non-uniform sampling is reduced, and the signal-to-noise ratio of magnetic resonance imaging is improved; therefore, the density adaptive k-space trajectory can be used for designing the k-space trajectory of sodium element magnetic resonance imaging.

[0043] 2, The density adaptive k-space trajectory proposed by the method satisfies the hardware limitations of maximum gradient strength and maximum gradient climbing rate. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 is a flowchart of the method of the present application;

[0045] Figure 2 is a k-space trajectory curve constructed by the present application;

[0046] Figure 3 is a schematic diagram of the sampling point distribution in the density adaptive k-space trajectory obtained by the method of the present application;

[0047] Figure 4 is the encoding gradient distribution obtained by the method of the present application; in the legend, Gx is the component of the encoding gradient on the x-axis, Gy is the component of the encoding gradient on the y-axis, and G is the absolute value of the encoding gradient strength.

[0048] Figure 5 is the complete density adaptive k-space trajectory obtained by the method of the present application;

[0049] Figure 6 is the Venn diagram of the complete density adaptive k-space trajectory obtained by the method of the present application;

[0050] Figure 7 is the density compensation coefficient of the complete density adaptive k-space trajectory obtained by the method of the present application;

[0051] Figure 8 is a comparison of reconstructed images and original images under different k-space trajectories; (a) is the original image, (b) is the reconstructed image corresponding to the radial trajectory, and (c) is the reconstructed image corresponding to the complete density adaptive k-space trajectory of the present application. DETAILED DESCRIPTION

[0052] In order to facilitate those skilled in the art to understand and implement the present application, the present application will be further described in detail below in conjunction with examples, and the examples described herein are only used to illustrate and explain the present application, and are not a limitation on the present application.

[0053] Example 1:

[0054] The density adaptive k-space trajectory design method for sodium element magnetic resonance imaging, as shown in Figure 1 includes the following steps:

[0055] Step 1, constructing k-space trajectory curve

[0056] The target of the density-adaptive k-space trajectory generation of the method of the present application is to complete the k-space center filling as soon as possible and to achieve uniform sampling at the k-space periphery. Based on the above requirements, the k-space trajectory curve in the present embodiment comprises a radial trajectory, the starting point of the radial trajectory is at the origin of the polar coordinate system at the k-space center, corresponding to the radial trajectory, a spiral trajectory is adopted at the k-space periphery, the radial trajectory and the corresponding spiral trajectory are connected through an arc trajectory, and the two ends of the arc trajectory are tangent to the radial trajectory and the spiral trajectory, respectively.

[0057] Step 1.1, determining the maximum value of the radial component of the k value in k-space.

[0058] The k value in the polar coordinate system comprises a radial component and an angular component, the k-space trajectory curve corresponding to the k-space trajectory proposed by the present application belongs to a non-Cartesian sampling trajectory, and a square field of view (FOV) and isotropic resolution are adopted. In the present embodiment, the sampling parameters are set as follows: the size of the field of view FOV is 32mm*32mm, the size of the sampling matrix Matrix is 64*64, and the preset number of excitations (Ns) is 64. Based on the above sampling parameters, the maximum value of the radial component of the k value of the sampling point in k-space is calculated as follows:

[0059]

[0060] wherein, k max represents the maximum value of the radial component, FOV1 represents the size of one dimensional component of the field of view, and Matrix1 represents the size of one dimensional component of the sampling matrix. Because a square field of view FOV and isotropic resolution are adopted, the same result can be obtained by taking the value of any one dimensional component of the field of view FOV and the sampling matrix Matrix.

[0061] Step 1.2, determining the defining condition of the arc trajectory in the k-space center region and the spiral trajectory at the k-space periphery for the k-space trajectory curve.

[0062] In the method of the present application, the defining condition is that the radial trajectory, the arc trajectory and the spiral trajectory under different excitations all satisfy the Nyquist sampling law, that is, the maximum distance of the adjacent sampling points on the k-space trajectory curve in k-space should not exceed 1 / FOV1. The included angle between the two adjacent k-space trajectory curves is 2π / Ns (Ns is the preset number of excitations), and the distance between the critical points of the arc trajectory and the spiral trajectory on the two adjacent k-space trajectory curves does not exceed 1 / FOV1, therefore, the radial component of the k value of the critical point of the arc trajectory and the spiral trajectory satisfies the following condition:

[0063]

[0064] wherein Ns represents a preset number of excitations, k swi represents a radial component of k value at the critical point where the circular arc trajectory switches to the spiral trajectory, denoted as critical point radial component k swi . When the radial component of k value is less than the critical point radial component k swi , the radial trajectory and the circular arc trajectory are adopted. Wherein the k-space radius of the circular arc trajectory is greater than 0 and less than In this embodiment, the k-space radius of the circular arc trajectory is set as k swi / 2. When the radial component of k value is greater than k swi , the spiral trajectory is adopted.

[0065] Step 1.3, setting the spiral trajectory.

[0066] The spiral trajectory is expressed as:

[0067] k(τ)=k max τe jωτ , k swi / k max ≤τ≤1………………3

[0068] wherein ω is the angular frequency of the spiral trajectory, τ is the time coefficient, k(τ) represents the k value corresponding to the time coefficient τ, and the spiral trajectory is equivalent to a circular trajectory with continuously increasing k-space radius. The spiral trajectory also needs to satisfy the Nyquist sampling law. Therefore, when the number of excitations is the preset number of excitations Ns, the k-space radius of the spiral trajectory needs to increase by 1 / FOV1 every time the arc increases by 2π / Ns. Therefore, it can be obtained that:

[0069] ω=k max / k swi …………………………………4

[0070] Under the angular frequency ω of the spiral trajectory defined in formula 4, the maximum distance between the spiral trajectories in multiple excitations is exactly 1 / FOV1.

[0071] Step 1.4, determining the center of the circular arc trajectory, and the point of intersection between the circular arc trajectory and the radial trajectory.

[0072] Based on the critical point radial component k swi and the spiral trajectory (including formula 3 and formula 4), the k value (equal to the polar coordinate value) of the critical point where the circular arc trajectory switches to the spiral trajectory is k swi ·e 1j , and the derivative of the spiral trajectory at the critical point is k max (1+j)e 1jSince one end of the circular arc trajectory is tangent to the spiral trajectory at the critical point, according to the k value of the critical point, the derivative at the critical point and the k space radius of the circular arc trajectory, the k value of the center of the circular arc trajectory can be calculated;

[0073] Based on the other end of the circular arc trajectory being tangent to the radial trajectory, the k value of the tangent point of the circular arc trajectory and the radial trajectory can be calculated according to the k value of the center of the circular arc trajectory and the k space radius of the circular arc trajectory.

[0074] Therefore, the k space trajectory curve designed by the scheme is: starting from the center of the k space polar coordinate system, switching from the radial trajectory to the circular arc trajectory at the tangent point of the radial trajectory and the circular arc trajectory; then, switching to the spiral trajectory at the tangent point of the circular arc trajectory and the spiral trajectory on the circular arc trajectory, and continuing to change on the spiral trajectory until the radial component of the k value reaches the maximum radial component k max , that is, the end point of the k space trajectory curve.

[0075] The corresponding k space trajectory curve under one excitation designed by the method is as shown in Figure 2 .

[0076] Step 2, obtaining a density-adaptive k space trajectory and an encoding gradient

[0077] The sampling points on the k space trajectory curve are selected, and the trajectory of the sampling points is the density-adaptive k space trajectory. The k space trajectory curve obtained in step 1 ensures that the distances between the corresponding k space trajectories under different excitations satisfy the Nyquist sampling theorem. In one excitation, the distances between the sampling points on the same k space trajectory curve also need to satisfy the Nyquist sampling theorem. In magnetic resonance imaging, the sampling points on the k space trajectory curve are discrete, and the sampling points are selected, that is, the k values of the sampling points are selected. The k value corresponding to each sampling point is the integral of the intensity value of the encoding gradient with time, and when the interval length of the change of the intensity value of the encoding gradient is Δt, the k value can be expressed as:

[0078]

[0079] N represents the number of the intensity values of the encoding gradient participating in the integration, and N increases correspondingly with the increase of the number of the time length Δt, N ∈ {1, 2, 3, … N′}, N′ represents the total number of the intensity values of the encoding gradient corresponding to the entire density-adaptive k space trajectory, since the k value is limited by the maximum radial component k max , the total number N′ of the intensity values of the encoding gradient is limited by the maximum radial component k max ; γ represents the gyromagnetic ratio of the nuclear spin; G i represents the i-th intensity value of the encoding gradient, which is denoted as the intensity value G i of the encoding gradient, and the intensity value G iIn a Cartesian coordinate system, including x-axis and y-axis components, the intensity value G of the encoded gradient is... i The magnitude is determined by both the x-axis and y-axis components, which are two orthogonal coordinate axes in a Cartesian coordinate system; i represents the index of the gradient intensity value encoded in chronological order; k N Let k represent the integral of the strength values ​​of the first N encoded gradients with respect to time. N Value; k N The radial component of the value is represented as the radial component k. Nr Due to limitations in magnetic resonance imaging hardware, the maximum intensity value and maximum gradient ascent rate that the system can achieve are finite. The intensity value G of the encoded gradient... i The upper limit of the encoded gradient strength G cannot be exceeded. max The strength value G of the encoded gradient i The rate of change per unit time is denoted as the climb rate SR, and the climb rate SR does not exceed the upper limit value of the climb rate SR. max .

[0080] With the goal of minimizing the sampling time, sampling points on the radial trajectory and the circular trajectory of the k-space trajectory curve are selected. The rule for selecting the corresponding sampling points is: from the starting time, the intensity value G of the encoded gradient is... i The size of the gradient varies with the maximum ascent rate if the strength value of the encoded gradient G is... i The size reaches the upper limit of the encoding gradient strength G. max Then the strength value G of the encoded gradient i The size remains the upper limit of the encoded gradient strength G. max At this point, the intensity value G of the radial trajectory, the circular trajectory, and the encoded gradient is... i The size of the value determines the k value of the corresponding sampling point.

[0081] To achieve uniform sampling density, sampling points are selected along a spiral trajectory within the k-space trajectory curve. The rule for selecting these sampling points is that the k-space distance between adjacent sampling points remains constant at a preset spatial interval. The k-value of the sampling point is then determined by the spiral trajectory and the preset spatial interval, which in turn determines the intensity value G of the encoding gradient. i .

[0082] Selecting sampling points on the k-space trajectory curve includes the following steps:

[0083] The intensity value G of the encoded gradient is calculated from the start time, at each time interval Δt. i Size based on the upper limit of the climb rate SR maxIncrementing, i represents the index of the intensity value of the encoded gradient in chronological order; simultaneously, at each time interval Δt, the intensity value G of the first N encoded gradients is calculated sequentially based on the radial and circular trajectories in the k-space trajectory curve, as well as Equation 5. i Integrating over time yields the corresponding value of k, denoted as k. N Value; k N The radial component of the value is denoted as k. Nr ;

[0084] As time increases, when the radial component k Nr Radial component k not reached the critical point swi (Corresponding to radial and circular trajectories), if the intensity value G of the encoded gradient... i The size reaches the upper limit of the encoding gradient strength G. max Then the strength value G of the encoded gradient i The size remains the upper limit of the encoded gradient strength G. max (The intensity value G of the encoding gradient corresponding to the circular trajectory) i The x-axis and y-axis components change over time, corresponding to a ramp rate SR = 0; thus ensuring the strength value G of the encoding gradient. i Not exceeding the upper limit of the coding gradient strength G max And the climb rate SR does not exceed the upper limit of the climb rate SR. max At the same time, the goal is to minimize the sampling time and complete the sampling of the k-space center as quickly as possible;

[0085] When the radial component k Nr Reaching or exceeding the critical point radial component k swi At time (corresponding to the spiral trajectory), the k-space distance between adjacent sampling points is constant at a preset spatial interval (the preset spatial interval is 1 / FOV1). Based on the spiral trajectory in the k-space trajectory curve and the preset spatial interval, sampling points are selected sequentially from the first sampling point on the spiral trajectory until the radial component k Nr Reaching or exceeding the maximum value k of the radial component max The selection of sampling points is stopped, and the corresponding density-adapted k-space trajectory ends, thus achieving the goal of uniform sampling density. Here, the k-value corresponding to the first sampling point on the spiral trajectory is denoted as k. n Then the k-space distance between the first sampling point on the spiral trajectory and the critical point is equal to k is the spatial distance from the last sampling point on the circular trajectory (i.e., the sampling point preceding the first sampling point on the spiral trajectory) to the critical point. n-1 This represents the k value corresponding to the last sampling point on the circular trajectory.

[0086] To meet the Nyquist sampling law, according to formula 5, the sampling bandwidth SW corresponding to the sampling point satisfies the following condition:

[0087] SW = max(G i ) · γ · FOV1 …………………… 6

[0088] Wherein, max(G i ) represents the maximum intensity value of the encoding gradient. When switching to the spiral trajectory, there is a case that the intensity value G i of the encoding gradient has not reached the upper limit value G max of the encoding gradient intensity, so it is possible that the maximum intensity value max(G i ) of the encoding gradient is less than the upper limit value G max of the encoding gradient intensity.

[0089] After sampling, the sampling point distribution of the k-space corresponding to a single excitation is shown in Figure 3 . The density adaptive k-space trajectory completes the discretization of the k-space trajectory and the density adaptive design.

[0090] According to the sampling point distribution in the density adaptive k-space trajectory, the intensity value of each encoding gradient can be calculated:

[0091] G N = (k N -k N-1 ) / (γΔt) …………………… 7

[0092] G N represents the Nth intensity value of the encoding gradient, which is denoted as the intensity value G N of the encoding gradient; k N-1 represents the k value corresponding to the integral of the intensity values of the previous N-1 encoding gradients with respect to time.

[0093] The intensity value of the encoding gradient corresponding to the density adaptive k-space trajectory of the method of the present application is shown in Figure 4 .

[0094] Step 3, obtain the complete density adaptive k-space trajectory of the entire k-space by k-space trajectory rotation

[0095] The sampling of one excitation only fills part of the k-space, the sampling density cannot be analyzed and the actual sampling needs to cover the whole k-space. In the embodiment, to analyze whether the sampling density of the density-adaptive k-space trajectory is consistent with the design, the complete density-adaptive k-space trajectory of the whole k-space after multiple excitations is obtained, and the complete density-adaptive k-space trajectory of the whole k-space is obtained by rotating the density-adaptive k-space trajectory determined in step 2 multiple times. One excitation corresponds to one density-adaptive k-space trajectory, and the preset number of excitations Ns in the embodiment is 64 set in step 1. The angle of each k-space trajectory rotation is 2π / Ns, and the preset number of excitations Ns is completed by k-space trajectory rotation, and the filling of the whole k-space is realized (as shown in Figure 5 ).

[0096] Step 4, analyze the sampling density of the density-adaptive k-space trajectory

[0097] Based on the complete density-adaptive k-space trajectory of the preset number of excitations Ns, the area occupied by each sampling point in the complete density-adaptive k-space trajectory of the whole k-space is analyzed using a Voronoi diagram (as shown in Figure 6 ). The sampling density at different k values is the reciprocal of the area occupied by the sampling point at different k values, and the sampling density distribution is calculated.

[0098] In non-Cartesian sampling, the sampling density is often corrected using a density compensation coefficient, which is the reciprocal of the sampling density at different k values and is proportional to the area occupied by the sampling point in the Voronoi diagram. The density compensation coefficient of the calculated density-adaptive k-space trajectory is as shown in Figure 7 It can be seen that when the radial component of the k value is less than the radial component of the critical point k swi , the sampling density is large and the density compensation coefficient is small. When the radial component of the k value is greater than the radial component of the critical point k swi , the sampling density remains basically unchanged, and uniform sampling density is achieved. The density-adaptive k-space trajectory of the method of the application realizes fast sampling in the center of the k-space and uniform sampling in the periphery of the k-space.

[0099] Step 5. Evaluate the imaging effect of the density-adaptive k-space trajectory

[0100] According to the sampling parameters of step 1, the imaging effect of the full density adaptive k-space trajectory is simulated using MATLAB, and is compared with the radial sampling trajectory. This embodiment evaluates the imaging effect by converting the original image into sampling data under the k-space trajectory, and then reconstructing the sampling data into an image. Under the condition of satisfying the Nyquist sampling theorem, the number of excitations using the full density adaptive k-space trajectory is 64, and the number of excitations using the radial sampling trajectory is 202. First, a sample original image with a spatial resolution of 64x64 is generated using the BART (Berkeley Advanced Reconstruction Toolbox) toolbox. Then, using the non-uniform fast Fourier transform (NUFFT) algorithm, the original image, the k-space trajectory (the k-space trajectory is the full density adaptive k-space trajectory and the radial sampling trajectory, respectively) is inputted to perform non-uniform fast Fourier transform (NUFFT) to obtain the corresponding sampling data. According to the full density adaptive k-space trajectory and the radial sampling trajectory, the density compensation coefficients of the density adaptive k-space trajectory and the density compensation coefficients of the radial trajectory are calculated. Then, the sampling data under the full density adaptive k-space trajectory, the full density adaptive k-space trajectory and the corresponding density compensation coefficients are inputted to perform inverse non-uniform fast Fourier transform (NUFFT) to obtain the reconstructed image corresponding to the full density adaptive k-space trajectory; the sampling data under the radial sampling trajectory, the radial sampling trajectory and the corresponding density compensation coefficients are inputted to perform non-uniform fast Fourier transform (NUFFT) to obtain the reconstructed image corresponding to the radial sampling trajectory. The sample original image and the reconstructed images corresponding to the two k-space trajectories are shown in FIG. 8. Figure 8 The number of excitations of the full density adaptive k-space trajectory of the method of the present application is 32% of the radial sampling, and the reconstructed image obtained is more uniform than the reconstructed image obtained by the radial sampling trajectory.

[0101] It should be noted that the embodiments described in the present application are only examples illustrating the spirit of the present application. Those skilled in the art of the present application can make various modifications or supplements to the described embodiments or use similar ways to replace them, without deviating from the spirit of the present application or exceeding the scope defined by the appended claims.

Claims

1. A method for density-adapted k-space trajectory design for sodium element magnetic resonance imaging, characterized in that, The method comprises the following steps: Step 1, constructing a k-space trajectory curve, the k-space trajectory curve comprising a radial trajectory, a starting point of the radial trajectory being at an origin of a polar coordinate system of the k-space, for the radial trajectory, a spiral trajectory is adopted at an outer periphery of the k-space, the radial trajectory and the corresponding spiral trajectory are connected through an arc trajectory, two ends of the arc trajectory being tangent to the radial trajectory and the spiral trajectory respectively; Step 2, selecting sampling points on the k-space trajectory curve, the trajectory of the sampling points is the density-adaptive k-space trajectory, and the k value of the sampling points is the intensity value G of the encoding gradient i Integral over time; The sampling points on the radial trajectory and the sampling points on the circular arc trajectory in the k-space trajectory curve are selected according to the following rules: from the starting time, the size of the intensity value G i of the encoding gradient varies at a rate of the upper limit value SR max of the climbing rate, if the size of the intensity value G i of the encoding gradient reaches the upper limit value G max of the encoding gradient intensity, the size of the intensity value G i of the encoding gradient remains the upper limit value G max of the encoding gradient intensity. Sampling points on the spiral trajectory in the k-space trajectory curve are selected according to the following rule: the k-space distance of adjacent sampling points remains constant.

2. The method for designing density-adapted k-space trajectories for sodium MRI of claim 1, wherein, Further comprising step 3: obtaining a complete density-adaptive k-space trajectory of the entire k-space through multiple k-space trajectory rotations on a density-adaptive k-space trajectory determined in step 2.

3. The method for designing density-adapted k-space trajectories for sodium MRI of claim 1, wherein, The step 1 comprises the following steps: Step 1.1, calculating a maximum value of a radial component of a k value of a sampling point in the k-space: k max FOV1 is the size of one dimension component of the field of view, and Matrix1 is the size of one dimension component of the sampling matrix. Step 1.2, determining that a radial component of a k value of a critical point at which the arc trajectory switches to the spiral trajectory is: Ns represents a preset number of excitations, k swi A radial component of the k value of the critical point at which the circular arc trajectory switches to the spiral trajectory is denoted as a critical point radial component k swi ; Step 1.3, setting the spiral trajectory as: k(τ) = k max τe jωt , k swi / k max ≤ τ ≤ 1 where ω is the angular frequency of the spiral trajectory, ω = k max / k swi ; τ is a time coefficient, k(τ) represents the k value corresponding to the time coefficient τ, and the k value includes a radial component and an angular component; Step 1.4, based on the radial component k of the critical point swi and the derivative of the spiral trajectory at the critical point; based on the k value of the critical point, the derivative at the critical point, and the k-space radius of the circular arc trajectory, the k value of the center of the circular arc trajectory is obtained; Based on the arc trajectory being tangent to the radial trajectory at the other end, the k value of the tangent point of the arc trajectory and the radial trajectory is obtained according to a k value of a center of the arc trajectory and a k-space radius of the arc trajectory.

4. The method for designing density-adapted k-space trajectories for sodium MRI of claim 3, wherein, The step 2 comprises the following steps: The intensity value G of the encoded gradient is calculated from the start time, at each time interval Δt. i Size based on the upper limit of the climb rate SR max Incrementing, where i represents the index of the gradient intensity value encoded in chronological order; simultaneously, at each time interval Δt, based on the radial and circular trajectories in the k-space trajectory curve, and using the intensity values ​​G of the first N encoded gradients... i Integrating over time yields the corresponding value of k, denoted as k. N Value; k N The radial component of the value is denoted as k. Nr ; when the radial component k Nr does not reach the critical point swi , if the magnitude of the intensity value G i of the encoding gradient reaches the upper limit value G max of the encoding gradient intensity, the magnitude of the intensity value G i of the encoding gradient is maintained as the upper limit value G max of the encoding gradient intensity, and the corresponding climbing rate SR = 0. When the radial component k Nr reaches or exceeds the critical point radial component k swi , the k-space distance between adjacent sampling points remains constant at a preset spatial interval, and sampling points are sequentially selected from the first sampling point on the spiral trajectory based on the spiral trajectory in the k-space trajectory curve and the preset spatial interval until the radial component k Nr of the sampling point reaches or exceeds the maximum radial component k max .

5. The method for designing density-adapted k-space trajectories for sodium MRI of claim 4, wherein, The step 2 by the intensity value G of the first N encoding gradients i With respect to the integration of time, the corresponding k value is obtained, based on the following formula: N denotes the number of the intensity values of the encoding gradient participating in the integration, N is increased correspondingly with the increase of the number of the length of time Δt, N ∈ {1, 2, 3, …N'}, N' denotes the total number of the intensity values of the encoding gradient corresponding to the whole density-adapted k-space trajectory; γ denotes the gyromagnetic ratio of the nuclear spin; G i denotes the i-th intensity value of the encoding gradient, denoted as the intensity value of the encoding gradient G i ; k N denotes the k value corresponding to the integration of the intensity values of the first N encoding gradients with respect to time, denoted as the k value; the radial component of the k N value is denoted as k N ; and k Nr .

6. The method for designing density-adapted k-space trajectories for sodium MRI of claim 4, wherein, The k-space distance between adjacent sampling points on the spiral trajectory is constant as 1 / FOV1, and the k value corresponding to the first sampling point on the spiral trajectory is denoted as k n The k-space distance between the first sampling point on the spiral trajectory and the critical point is equal to wherein is the k-space distance between the last sampling point on the circular arc trajectory and the critical point, k n-1 is the k value corresponding to the last sampling point on the circular arc trajectory, and the last sampling point on the circular arc trajectory is the sampling point before the first sampling point on the spiral trajectory.

7. The method for designing density-adapted k-space trajectories for sodium MRI of claim 5, wherein, A sampling bandwidth corresponding to the sampling point is: SW = max(G i ) · γ · FOV1 SW denotes the sampling bandwidth; max(G i ) denotes the maximum intensity value in the encoding gradient.

8. The method for designing density-adapted k-space trajectories for sodium MRI of claim 4, wherein, The step 2 further comprises the following steps: Calculating intensity values of various encoding gradients: G N = (k N -k N-1 ) / (γΔt) G N represents the Nth intensity value encoding the gradient; k N-1 represents the k value corresponding to the integration over time of the intensity values of the previous N-1 encoding gradients.

9. The method for designing density-adapted k-space trajectories for sodium MRI of claim 2, wherein, An angle of each k-space trajectory rotation in the step 3 is 2π / Ns, Ns representing a preset number of excitations.

10. The method for designing density-adapted k-space trajectories for sodium MRI of claim 4, wherein, The k-space radius of the circular arc trajectory is greater than 0 and less than

Images