Magnetic particle reconstruction method based on improved magnetization model and prior information driving

Abstract

A magnetic particle reconstruction method based on an improved magnetization model and prior information driving accurately simulates the nonlinear magnetization behavior of magnetic nanoparticles through the improved magnetization model, designs an objective function through a double prior constraint strategy in combination with the inherent physical characteristics of the magnetic nanoparticles, and performs efficient iteration solution on the objective function to construct a depth expansion model. And accurate reduction of concentration distribution of the magnetic nanoparticles is realized. Through the magnetization model, a measurement voltage and a system matrix are obtained. And an objective function fusing the L1 regularization prior and the Laplacian regularization prior is constructed. An iterative solving process of the objective function is converted into a hierarchical structure of a deep network, so that a deep expansion model is constructed. And the high-fidelity reconstruction of the MNPs concentration distribution is realized. According to the method, the reconstruction quality of the MPI image is remarkably improved, and the method is of great significance in reducing system hardware requirements and promoting biological application.

Description

Technical Field

[0001] The invention belongs to the technical field of biomedical imaging, and in particular relates to a magnetic particle reconstruction method based on an improved magnetization model and prior information driving. Background Art

[0002] Magnetic nanoparticle imaging (MPI) is a high-contrast, high-sensitivity, and high-spatial-resolution molecular imaging technique that uses superparamagnetic iron oxide nanoparticles as tracers and uses the nonlinear magnetization characteristics of the particles to display molecular and protein abnormalities in the early stages of the disease, thereby achieving precise diagnosis and treatment. The spatial resolution of MPI directly determines the imaging quality and the accuracy of diagnosis and treatment, so it is necessary to propose new methods to improve the accuracy of spatial reconstruction. In MPI, a field-free region (FFR) is created using a gradient field generated by a permanent magnet or a powered coil. A sinusoidal oscillating drive field is then superimposed to allow the magnetic field-free region to scan the field of view (FOV) along a predetermined trajectory (such as a Lissajous trajectory). Within the FFR, the nanoparticles can adjust their magnetization according to the driving magnetic field, and the receiving coil captures this nonlinear response to generate a voltage signal. Then, the voltage signal is used to visualize the nanoconcentration distribution in the area to be measured, thereby directly imaging the interior of the object being measured.

[0003] The accuracy of magnetic nanoparticle imaging technology is closely related to the construction method of the system matrix, which consists of the Fourier components of the time signal generated by each magnetic nanoparticle at all possible spatial positions, representing the relationship between the measured signal and the nanoparticle distribution. The accuracy and quality of SM-based MPI reconstruction depends largely on how the SM is acquired. At present, there are two main methods for constructing SM: measurement-based methods and model-based methods. When the SM measurement records the corresponding signal, an incremental MNP sample must be repeatedly moved on each voxel within the FOV. In order to improve the quality of SM measurements, multiple averaging methods are usually used, which significantly prolongs the calibration time. In addition, when the properties of the tracer or the magnetic field environment change, the SM needs to be recalibrated, which brings a lot of labor and time costs. Model-based methods usually construct SM using mathematical models that can represent the physical properties of the imaging system and nanoparticles. These models take into account factors such as magnetic field distribution, nanoparticle characteristics and scanning procedures, without the need for a large number of measurements. The feasibility of obtaining SM based on the Langevin model has been confirmed in many studies. However, the Langevin model ignores the relaxation effect and interaction of SPIOs, which poses a challenge to accurately simulate the magnetization intensity of SPIOs in actual situations. Therefore, introducing a suitable model to accurately describe the dynamic magnetization behavior of magnetic nanoparticles in MPI is crucial for MPI voltage analysis and improving the resolution of MPI reconstruction results.

[0004] Under ideal conditions, SM-based MPI image reconstruction can be considered as a linear inverse problem, assuming that the complex interactions between particles in the concentration distribution are minimized. Various traditional regularization solutions, especially those related to iterative methods, such as the well-known Kaczmarz algorithm and the conjugate gradient method (CGS) algorithm, have been proposed for MPI reconstruction. In addition, advanced regularization models such as the non-negative fused LASSO technique, non-convex regularization techniques including minimum-maximum concavity, and TV penalty regularization techniques have been introduced to enhance image clarity and uniformity. Although these methods improve reconstruction accuracy, they usually require complex matrix operations and manual parameter adjustments, which may affect stability and image quality.

[0005] In recent years, it has become a trend to apply deep learning-based algorithms to MPI reconstruction. The Deep Image Prior (DIP) method enhances generalization by leveraging an untrained network learned during testing. Another approach, the Plug-and-Play MPI (PP-MPI) method, involves pre-training image denoising and combining it with the Split-Merge (SM) algorithm to achieve a balance between efficiency and generalization. However, these methods may face performance limitations when transferring priorities across different tasks. In addition, the source of large amounts of labeled data and poor model interpretability are also major obstacles. Summary of the invention

[0006] In order to overcome the shortcomings of the above-mentioned prior art, the present invention proposes a magnetic particle reconstruction method based on an improved magnetization model and prior information driven, the processing object is based on magnetic nanoparticle imaging technology, the improved Jiles-Atherton hysteresis function is used as the magnetic nanoparticle magnetization model, the gradient field is used as the selection field, the uniform magnetic field with sinusoidal variation is used as the driving field, the induced voltage generated by the magnetization is collected by the receiving coil, and then the prior information (saturation magnetic moment, particle concentration, vacuum permeability, etc.) in the magnetic nanoparticle imaging process is combined, and the magnetic nanoparticle concentration is reconstructed through the deep expansion model, so as to quantitatively visualize the magnetic particle concentration in the biological tissue. This method can be used to solve the inverse problem of MPI concentration reconstruction, solve the time-consuming and labor-intensive shortcomings of the iterative process and manual parameter adjustment in the traditional reconstruction algorithm, and perform efficient and stable MPI image reconstruction to obtain high-quality image reconstruction results.

[0007] In order to achieve the above object, the technical solution adopted by the present invention is:

[0008] A magnetic particle reconstruction method based on an improved magnetization model and driven by prior information comprises the following steps:

[0009] Step 1: Based on the basic theory of magnetic nanoparticle imaging, combined with the improved Jiles-Atherton magnetization model, the system matrix (a matrix representing the relationship between magnetic field and magnetization intensity) and the measured voltage of particle distribution (reflecting the induced voltage caused by the change of magnetization intensity) are obtained;

[0010] Step 2: Considering the unique properties of the magnetic nanoparticle solution space and the intrinsic correlation of magnetization between adjacent voxels in space, two types of prior information are designed. One is the L1 regularization prior, which aims to promote the sparsity of the solution, effectively suppress noise interference, and improve the quality of the reconstructed image; the other is the Laplacian regularization prior of the magnetization weight. This constraint further enhances the local coherence and spatial consistency of the reconstruction result by describing the relationship between the magnetization gradients in the space of adjacent particles.

[0011] Step 3: Based on the prior information proposed in step 2, a double prior constraint objective function is constructed. This objective function aims to balance the relationship between data fidelity and the sparsity and smoothness of the solution, laying a solid foundation for subsequent optimization solutions;

[0012] Step 4: In order to efficiently iteratively solve the objective function, a deep unfolding model is constructed based on non-convex optimization and iterative approximation theory to achieve high-fidelity reconstruction of the concentration distribution of magnetic nanoparticles;

[0013] Step 5: In order to verify the practical application effect of the proposed method, a reasonable loss function is designed to train the model and predict it using the test data set. At the same time, the predicted image and the real image are displayed.

[0014] Further, the specific steps of step 1 are as follows:

[0015] In the process of magnetic nanoparticle imaging (MPI), the system mainly consists of three parts: selection field, driving field and receiving coil. First, the selection field is a gradient field that produces a zero field region (FFR). The particles in this region have not reached magnetization saturation, so there is a prerequisite for generating a magnetization response. The region that generates a zero magnetic field is generally divided into a zero magnetic field point (FFP) and a zero magnetic field line (FFL). The present invention is based on FFP for research. Secondly, a time-varying driving field is superimposed to make the FFR move within the field of view (FOV). The driving field serves as an excitation field for the magnetic nanoparticles in the FFR, causing the particles to generate a magnetization response, and the receiving coil generates an induced voltage signal according to the change in magnetic flux density. Based on the basic principles of magnetism and the specific characteristics of magnetic nanoparticles, the present invention uses an improved Jiles-Atherton magnetization model to accurately model the magnetization behavior of magnetic nanoparticles in MPI, so as to obtain the measured voltage of the system matrix and particle distribution. First, considering the complex interaction between magnetic nanoparticles and the external field, the total effective magnetic field H is defined. e for:

[0016] H e =H+αM (1)

[0017] Where H is the external field, α is the coupling parameter, and M is the magnetization intensity of the magnetic nanoparticle. Secondly, the Langevin function is used as the hysteresis-free magnetization function to convert the hysteresis-free magnetization M an With effective magnetic field H e The relationship is expressed as:

[0018]

[0019] in, is the Langevin parameter, Ms is the saturation magnetization of the magnetic nanoparticle. rev Related to the domain wall bending, the irreversible magnetization M irr Related to domain wall displacement. Magnetic nanoparticles are aligned along the external magnetic field, resulting in domain wall displacement. Therefore, the irreversible magnetization intensity M irr Dominant, modified Jiles-Atherton magnetization model M JA It can be written as:

[0020]

[0021] In order to simplify the model, we can assume that the reversible parameter c is zero. In this case, M JA =M irr . M JA The model can be written as an ordinary differential equation describing the magnetization intensity of the particle under an external magnetic field H:

[0022]

[0023] in, The parameters α and k can be determined by the coercive point (H = H e , M=0) and remanence point (H=0, M=M r ) is determined. At the coercive point of δ = +1, M JA The model equation (4) is:

[0024]

[0025] Similarly, at the remanence point of δ = -1, (4) can be written as:

[0026]

[0027] By combining equations (5) and (6), the parameters α and k can be determined. Solving the ordinary differential equation, M JA The relationship between and H can be determined as:

[0028]

[0029] Based on formula (7), the induced voltage at r of a receiving coil with a sensitivity of p(r) can be expressed as:

[0030]

[0031] Where, M(r,t)=c(r)·M JA (r, t), c(r) is the particle concentration at position r. Define s(r, t) as the integral kernel describing the mapping between particle concentration and voltage:

[0032]

[0033] Then formula (8) can be expressed as:

[0034] u(t)=∫ Ω c(r)·s(r,t)dr (10)

[0035] In order to perform numerical calculations, the time and space of formula (10) are discretized:

[0036]

[0037] Among them, s(r n ,t k ) is calculated as:

[0038]

[0039] Where m represents the modulus of the magnetic moment of a single particle. Performing fast Fourier transform (FFT) on both sides of formula (11) in the time domain yields:

[0040]

[0041] Among them, S(f k ,r n ) is defined as the system matrix, is defined as the measured voltage signal spectrum vector of MNPs. k=1,2,3...,K is the number of frequency components, and n=1,2,3...,N is the number of voxel points on the imaging grid.

[0042] Further, the specific steps of step 2 are as follows:

[0043] Considering the unique properties of the magnetic nanoparticle solution space and the intrinsic correlation of magnetization between adjacent voxels in space, two prior information are designed. One is the L1 regularization prior, which aims to promote the sparsity of the solution, effectively suppress noise interference, and improve the quality of the reconstructed image; the other is the Laplacian regularization prior of the magnetization weight. This constraint further enhances the local coherence and spatial consistency of the reconstruction result by describing the relationship between the magnetization gradients in the space of adjacent particles.

[0044] First, formula (13) is described in matrix form as:

[0045]

[0046] in, is the system matrix, is the particle concentration distribution. is the measured signal spectrum vector of MNPs. K is the number of frequency components, and N is the number of voxel points on the imaging grid.

[0047] Second, since K is much smaller than N, the problem of solving the particle concentration distribution in formula (14) can be converted into solving an ill-posed inverse problem, which is generally solved by the least squares method plus regularization constraints, which can be expressed as:

[0048]

[0049] in, To ensure that the reconstruction result C is as close to the true value as possible, R(C) is the regularization constraint. By introducing additional constraints on the fidelity term of the inverse problem, the ill-posedness of the inverse problem can be overcome, thus obtaining a more stable and reliable solution.

[0050] In magnetic particle imaging technology, the generation of induced voltage mainly depends on the dynamic movement of the zero magnetic field area in the field of view (FOV), which stimulates the magnetization response of magnetic nanoparticles. It is worth noting that the actual distribution of magnetic nanoparticles in the field of view is often highly sparse, that is, they usually only occupy a very small number of voxel points in the FOV. This highly sparse characteristic places special requirements on the MPI image reconstruction algorithm. It is necessary to design a method that can effectively handle large-scale data sparsity to ensure that the spatial distribution of magnetic nanoparticles is accurately reconstructed from limited magnetization response data. Therefore, in view of the sparse distribution of magnetic nanoparticles in the field of view (FOV), the L1 regularization term is used as a priori constraint. The L1 regularization term helps to suppress the interference of noise during the reconstruction process by promoting the sparsity of the solution, that is, reducing the number of non-zero elements, thereby improving the quality and accuracy of the reconstructed image. For formula (15), R(C) can be expressed as follows:

[0051]

[0052] In addition, it must be pointed out that although the regularization term that promotes sparsity can effectively reduce non-zero elements and thus maintain the sparsity of the solution space, excessive sparsification may lead to the loss of the reconstructed target morphology, that is, the blurring of key details and structures. In order to balance sparsity and reconstruction quality, a strategy combining spatial prior information is further proposed to ensure that the reconstructed image has both morphological accuracy and rich details. In the spatial encoding mechanism of MPI, there is a significant spatiotemporal correlation between the magnetic nanoparticles in the zero magnetic field region (FFR) and the neighboring particles. In particular, at the moment when the central voxel point is spatially encoded, its surrounding voxel points will also produce weak magnetization responses, which are effectively captured during the signal acquisition process. Based on this high correlation, a magnetization weight prior constraint based on the Laplacian matrix is ​​adopted. This constraint term aims to accurately characterize the interaction and dependence between adjacent voxel points, so as to maximize the preservation and restoration of the morphological characteristics and detail information of the reconstructed image while ensuring sparsity.

[0053] The Laplacian matrix can be represented by L, where L = DW, W is the adjacency matrix, which represents the intensity relationship between each voxel point and each other point. D is the weighted degree matrix, where each element D ii Represents the total connection weight between voxel point i and other points. Through such a definition, the Laplacian matrix can effectively characterize the local connectivity and differences between voxel points in the image, and provides a mathematical basis for further using spatial prior information to constrain and optimize the morphology and details of the reconstructed image. First, in MPI, it is generally believed that if the applied external magnetic field makes the magnetization intensity of the particle reach 80% of the saturation magnetization intensity, the particle reaches magnetization saturation, and the modulus value of the external magnetic field intensity H is defined as the saturation field intensity H sat After the selection field and the driving field are superimposed, the particle magnetic field at the FFR is H(r, t) = H D (t)+H s (r) = 0, then the FFR surrounding particles r Δ The combined magnetic field can be expressed as:

[0054]

[0055] Where Δt is the displacement vector of the coordinate position around the FFR relative to the coordinate position of the FFR. Δ ,t)|| <H sat When , the spatial range of all components is defined as the δ neighborhood of FFR. For the relationship weight between each voxel point and FFR in the neighborhood, M JAThe derivative of the model is used to measure it, because formula (12) shows that the magnetization response generated at each coordinate position r in space is proportional to the rate of change of the magnetization intensity M(r,t) with respect to time, where Therefore, when the FFR position is determined, the spatial encoding time t FFR is also a fixed value, M JA The derivative of (β||H(r,t)||) with respect to H(r,t) reflects the intensity of the magnetization response of the entire space at that moment. The adjacency matrix W can be expressed as:

[0056]

[0057] Where i and j represent each voxel point. δ(i) represents the neighborhood of point i when it is the FFR point. Δr ij Represents the displacement vector from point j to i in the neighborhood. Indicates that H = G·Δr ij Time JA The derivative of (β||H||) with respect to H. Indicates that when point i is the FFR point, M JA The derivative of (β||H||) with respect to H is β / 3. For the determination of the neighborhood of δ(i), if the modulus of the applied magnetic field intensity at the voxel point i is less than H sat , then the voxel point is in the neighborhood δ(i) of FFR point i, that is, by solving 0.8M sat =cmM JA (β||H||) Solve for H sat Since the coordinates of voxel points are discrete during numerical calculation, to determine whether point j is in the neighborhood of i, it is only necessary to judge ||G·Δr based on the selection field gradient G. ij || Is it less than H sat For a selection field with good linearity, the number and distribution of voxel points in the magnetization space neighborhood of each voxel point are the same. Therefore, the calculation of the relative weight in the neighborhood of each voxel point can be simplified to only calculate the weight of each point in the neighborhood of a voxel point in space, and then match all the weights of the neighborhood of the point to other magnetization space neighborhood points to the voxel point according to the relative position in space, thereby completing the construction of the adjacency matrix W. The weighted degree matrix D can be expressed as:

[0058]

[0059] According to formula (18) and formula (19), the element L of the Laplace matrix L is ij It can be expressed as:

[0060] L ij =D ij -W ij (20)

[0061] According to the definition of D and W, L ij The value of can be divided into three cases:

[0062] 1) When i=j:

[0063]

[0064] Note that W ii =0, therefore,

[0065]

[0066] For clarity, the index in the summation was changed from j to k to avoid confusion with the external j, but in practice they are identical.

[0067] 2) When i≠j and i and j are adjacent (j∈δ(i)):

[0068] L ij =D ij -W ij =0-W ij =-W ij (twenty three)

[0069] therefore,

[0070]

[0071] 3) When i≠j and i and j are not adjacent

[0072] L ij =D ij -W ij =0-0=0 (25)

[0073] Based on the above, the expression of the Laplace matrix L is:

[0074]

[0075] For N individual voxel points, L is a symmetric matrix of size N×N symmetric about the main diagonal. So far, for formula (15), R(C) can be further expressed as follows after adding the Laplace prior constraint:

[0076]

[0077] Further, the step three comprises the following steps:

[0078] Based on the prior information established in step 2, a double prior constraint objective function is constructed. The core function of this objective function is to balance the delicate relationship between data fidelity (i.e., ensuring the consistency of the reconstruction results with the original measurement data) and the sparsity and smoothness of the solution. Therefore, based on the sparse prior constraint term and the spatial prior constraint term of formula (27), formula (15) can be written as follows:

[0079]

[0080] in, is the optimal approximate solution of the concentration vector C, λ1 and λ2 are regularization parameters, ||C||1 is the L1 regularization sparse constraint term, is the Laplace space constraint term. This objective function introduces sparse constraints according to the imaging characteristics of MPI to promote the sparsity of the solution and improve the anti-noise performance of the reconstructed image; introduces the Laplacian regularization term based on the magnetization weight to retain the neighborhood information of the pixel point and restore the morphological information of the image. By combining these regularization terms to balance the sparsity and morphology of the reconstructed image, the imaging quality and spatial resolution can be improved.

[0081] Further, the step 4 comprises the following steps:

[0082] In order to efficiently iteratively solve the objective function formula (28), based on non-convex optimization and iterative approximation theory, the present invention constructs a deep unfolding model to achieve high-fidelity reconstruction of the concentration distribution of magnetic nanoparticles. First, the gradient calculation is performed. For the first term of formula (28) We can use the chain rule to differentiate C. Expanding this term gives:

[0083]

[0084] The result of the derivation of C can be expressed as:

[0085]

[0086] The second term α||C||1, its partial derivative is non-smooth. The sub-gradient of the L1 norm can be expressed as:

[0087]

[0088] Where sign(C) is the sign function, which means that the derivative for each element is +1 (if the element is positive), -1 (if the element is negative), and 0 (if the element is 0).

[0089] Item 3 Expand this item to:

[0090] ω(LC)T (LC) = ωC T L T LC (32)

[0091] The result of the derivation of C can be expressed as:

[0092]

[0093] Combining the partial derivatives, the overall gradient formula can be expressed as:

[0094]

[0095] Based on non-convex optimization and iterative approximation theory, the optimal value of the particle information updated in the iterative process can be found after multiple iterations. The k+1th iteration can be expressed as:

[0096] C (k+1) =C (k) -[S T (SC (k) -U)+α·sign(C (k) )+2ωL T LC (k) ] (35);

[0097] Next, the core construction of the entire deep expansion reconstruction framework is based on the iterative formula (35), which transforms the iterative solution process into a hierarchical structure of the deep network. Each iteration is mapped to a layer in the deep network, and the layer-by-layer iterative approximation guides the advancement of the entire reconstruction process. The expanded deep neural network model can be formally expressed as:

[0098] C (k+1) =Ψ θ (C (k) ;U,S,L) (36)

[0099] Among them C (k) represents the concentration estimate of the kth iteration, and θ is the set of model parameters shared in all iteration steps. This model takes the Laplace matrix L, the complex measured voltage U and the system matrix S as input information, and the particle concentration C is initialized to all zeros in the initial stage to provide a benchmark starting point for the subsequent iteration process.

[0100] The key components of the model are as follows:

[0101] Gradient descent optimization strategy: In each iteration step, the concentration estimate C is updated using the gradient descent method. (k) Specifically, by calculating the negative gradient direction of the fidelity term, the result of each iteration is guided to move closer to the true solution. This process ensures the convergence of the iterative process and the accuracy of the reconstruction results.

[0102] Complex result processing mechanism: Since the result generated by formula (34) is in complex form, a channel and dimension conversion layer is used. This layer is responsible for converting the one-dimensional complex output generated by each iteration into a two-dimensional data format, and cleverly separates the real and imaginary parts, and allocates independent channels for processing. This design not only simplifies the complexity of complex operations, but also improves the model's ability to process complex data.

[0103] Feature extraction and enhancement: In order to effectively suppress noise and redundant information while ensuring the quality of reconstruction, a near-end projection layer consisting of an encoder, a decoder, and a spatial attention module is introduced. The encoder and decoder, as feature extractors, can effectively capture and reconstruct key feature representations in the image. The spatial attention module further enhances the model's sensitivity to core feature information. By dynamically adjusting the weights of different spatial positions, the model can focus on the areas that are most critical to the reconstruction results, thereby suppressing irrelevant information while improving reconstruction accuracy.

[0104] The calculation process of the model is as follows:

[0105] First, the iterative optimization process is explicitly expressed as the structure of a deep neural network, so that the basic idea of ​​gradient descent is used to optimize formula (35). In the kth iteration, the fidelity term is calculated relative to the current solution term (denoted as C (k) Gradient in Represents the fidelity loss function. Subsequently, the solution term is moved along the negative gradient direction of the fidelity term by a small step size η (learning rate), which can be expressed as:

[0106]

[0107] In this process, the results of the iteration are represented by complex numbers. Secondly, in the channel and dimension conversion layer, the one-dimensional data is converted into two-dimensional data, and the real part and the imaginary part are extracted into two different channels respectively. This provides a more intuitive and easier to process data format for subsequent network layers. The conversion process can be expressed as:

[0108]

[0109] Second, in order to effectively suppress noise and redundant information while ensuring the quality of reconstruction, the near-end projection layer integrates the encoder, decoder, and spatial attention module. The encoder part uses a batch normalization layer to accelerate convergence and stabilize gradients, and improves the complexity and depth of feature representation through stacked convolutional layers and ReLU activation function structures. The size of the convolution kernel is 3, and the number of output channels is 2, 4, 4, 8, 8, and 16, respectively. Assume that the input feature is The encoding process can be expressed as:

[0110]

[0111] Among them, h i is the feature representation after the i-th convolution layer. Next, the spatial attention module is introduced to add weights to the particle features at each position to further enhance the feature representation.

[0112] g6=F SAM (h6) (40)

[0113] Among them, F SAM (·) represents the spatial attention calculation. Then, the decoder gradually reduces the high-level enhanced feature representation through the six-layer convolution operation with the same structure as the encoder, setting the number of reverse mapping channels to 16, 8, 8, 4, 4, 2, which can be expressed as:

[0114] g i =ReLU(Conv(g i+1 )),i=5,4,…,0 (41)

[0115] Among them, g i is the feature representation after the i-th convolution layer. Fourth, in order to further improve the accuracy of the prediction and balance the tendency of the predicted value to approach the target, the output of the decoder is converted into a complex form through the back channel and dimension conversion layer, and the predicted value is updated using the residual connection strategy according to the iterative formula (35). This step helps the network retain the memory of the initial value during the iteration process and gradually improves the accuracy of the prediction. The number of iterations is preset to k = 5. Finally, after the entire iteration is completed, the internal covariate bias is removed through the BN layer to stabilize and accelerate the training process, and through two convolution operations, the nonlinear representation ability of the network is further improved to make the prediction result close to the target result.

[0116] In addition, it is supplemented that the spatial attention module F SAM (·) receives a feature map X as input, whose size is C×H×W, where C is the number of channels, H and W are the height and width of the feature map respectively. The feature map is The output is the result of The module re-weights different positions of the input feature map by calculating the weight map in space, thereby enhancing the focus on important areas. First, the spatial information of each channel is compressed into a global scalar using global average pooling and global maximum pooling.

[0117]

[0118] Among them, f avg and f maxThe two vectors are transformed by the convolutional layer to extract more complex features. The transformed vector is denoted as f avg and f avg .

[0119]

[0120] Among them, W avg and W max is the weight matrix, b avg and b max is the bias term, and σ is the ReLU nonlinear activation function. Then, the transformed vectors are concatenated element-wise and passed through an additional convolution and a Sigmoid activation function to generate a spatial attention weight map M with a size of 1×H×W.

[0121]

[0122] Finally, the attention weight map M is multiplied element-wise with the original input feature map X to obtain the weighted feature map X′.

[0123] X′(c,i,j)=X(c,i,j)×M(i,j) (45);

[0124] Further, the step five comprises the following steps:

[0125] The weighted sum of mean square error (MSE) and structural similarity index (SSIM) is used as the loss function for network model training. The total loss function can be expressed as:

[0126] Loss=λ·MSE+(1-λ)·(1-SSIM) (46)

[0127] Among them, λ = 0.9 is a weight parameter used to balance the relative importance between MSE and SSIM. MSE is used to evaluate the quality of the reconstructed image. SSIM is used to evaluate the similarity between the reconstructed image and the original image. MSE can be expressed as:

[0128]

[0129] Among them, Y and are the original and network predicted MPI images, respectively, and m is the number of images. SSIM can be expressed as:

[0130]

[0131] Among them, Y and are the original and network predicted MPI images respectively. Y and are image Y and The mean of . and They are image Y and The variance of . is the image Y and The covariance of . c1 and c2 are constants for stability, which are set to 0.01 and 0.03 for normalized images. Then, the trained deep neural network model is tested on the test dataset to verify the performance of the model in practical applications. At the same time, the error graph between the prediction results and the real image is shown.

[0132] Compared with the prior art, the present invention has the following advantages:

[0133] First, the present invention introduces an improved Jiles-Atherton hysteresis model to replace the Langevin function to more accurately simulate the magnetization behavior of SPIOs in the MPI system, thereby improving the simulation accuracy of the measured voltage and the system matrix.

[0134] Second, the present invention combines traditional iterative algorithms with deep learning to optimize algorithm performance by training implicit functions. Deep learning parameter training reduces the difficulty of manual parameter adjustment to adapt to different image reconstruction tasks and requirements.

[0135] Third, the present invention takes into account the spatial correlation between adjacent voxel points in the MPI magnetization space, which is ignored by traditional iterative algorithms. By introducing the Laplacian regularization term based on the magnetization weight, the correlation in the magnetization space can be more accurately described and simulated, thereby improving the accuracy and robustness of image reconstruction to a certain extent. BRIEF DESCRIPTION OF THE DRAWINGS

[0136] Figure 1 The schematic diagram of the MPI imaging system of the present invention, wherein (a) is an MPI imaging device system, and (b) is an MPI signal generation process and its spectrum diagram;

[0137] Figure 2 This is a flowchart illustrating an embodiment of the present invention;

[0138] Figure 3 This is the overall network diagram described in the embodiment of the present invention;

[0139] Figure 4 It is the spatial magnetization response generated by the center position of FOV and its adjacent voxel points during the driving field period of the present invention;

[0140] Figure 5This is a reconstruction result diagram based on the present invention. The first row in the figure is the original label result, and the second row is the result reconstructed by the present invention. DETAILED DESCRIPTION

[0141] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0142] like Figure 1 As shown in the figure, the MPI system schematic diagram is shown. Figure 1 As shown in (a), the MPI imaging system is mainly composed of a selection field, a driving field and a receiving coil. The selection field in the Z direction is generated by a pair of anti-Helmholtz coils, and a point with a magnetic field intensity of 0 is formed in the center, which is called FFP. Only the magnetic particles near the FFP will generate an induction signal, while the magnetic particles far away from the FFP are in a saturated state, so spatial positioning can be achieved through the FFP. The alternating magnetic fields (driving fields) in the X and Y directions are generated by a pair of Helmholtz coils, which are used to control the motion trajectory of the FFP and the magnetization response of the magnetic nanoparticles. Figure 1 (b) The curve in the upper right corner shows the magnetic response curve of magnetic nanoparticles under the excitation of an alternating magnetic field. The spectrum distribution of the signal is obtained by Fourier transform. Figure 1 (b) As shown in the lower right curve, the spectrum shows that M(t) contains not only the frequency component of the driving frequency but also a series of harmonic components. In MPI, the change in magnetization intensity is detected by measuring the voltage induced in the receiving coil.

[0143] like Figure 2 The process of the reconstruction method of the present invention is as follows:

[0144] S101: obtaining the measured induced voltage U and the system matrix S according to the improved magnetization model;

[0145] S102: Considering the unique properties of the magnetic nanoparticle solution space and the intrinsic correlation of magnetization between adjacent voxels in space, two types of prior information are designed respectively;

[0146] S103: Based on the prior information, construct an objective function with double prior constraints;

[0147] S104: Based on non-convex optimization and iterative approximation theory, a deep unfolding model is constructed, and U, S, and L are input into the model to iteratively solve the objective function;

[0148] S105: Design a reasonable loss function to train the model and use the test data set for prediction to obtain the reconstruction result.

[0149] The implementation process of step S101 is as follows:

[0150] According to the magnetization process of MPI, the system matrix S is constructed. The size of A is 2502×784. Save it as a mat file and set the key name to "S". Select 15,000 28×28 MNIST data sets and multiply them with the system matrix S to obtain the measured voltage signal U. Save it as a mat file and set the key name to "U". Save the used MNIST data set as label data as a mat file and set the key name to "labels". Among them, 10,000 model training data and 5,000 model test data are divided.

[0151] Step S102 specifically includes:

[0152] (1) Considering the unique properties of the solution space of magnetic nanoparticles, the L1 norm is designed as sparse prior information to ensure the sparsity of the solution space while taking into account the anti-noise performance of the model;

[0153] (2) Considering the intrinsic correlation of magnetization between adjacent voxels in space, the spatial prior information of Laplace magnetization weights is designed to establish the relationship between adjacent voxel points;

[0154] As mentioned above, the problem of solving the spatial concentration distribution of magnetic particles can be converted into solving an ill-posed inverse problem, which is generally solved by the least squares method plus regularization constraints, which can be expressed as:

[0155]

[0156] in, To ensure that the reconstruction result C is as close to the true value as possible, R(C) is the regularization constraint. By introducing additional constraints on the fidelity term of the inverse problem, the ill-posedness of the inverse problem can be overcome, thus obtaining a more stable and reliable solution.

[0157] In MPI, the distribution of magnetic nanoparticles usually occupies only a small number of voxel points in the field of view (FOV), and the number of non-zero voxels is much smaller than the total number of primes. Therefore, we use the L1 norm as a regularization constraint to ensure the sparsity of the solution space while taking into account the anti-noise performance of the reconstruction method. R(C) can be expressed as follows:

[0158] R(C)=||C||1

[0159] In addition, it must be pointed out that although the regularization term that promotes sparsity can effectively reduce non-zero elements and thus maintain the sparsity of the solution space, excessive sparsification may lead to the loss of the reconstructed target morphology, that is, the blurring of key details and structures. In order to balance sparsity and reconstruction quality, we further propose a strategy that combines spatial prior information to ensure that the reconstructed image has both morphological accuracy and rich details. In the spatial encoding mechanism of MPI, there is a significant spatiotemporal correlation between magnetic nanoparticles in the zero magnetic field region (FFR) and neighboring particles. In particular, at the moment when the central voxel point is spatially encoded, its surrounding voxel points will also produce weak magnetization responses, which are effectively captured during the signal acquisition process. Based on this high correlation, we innovatively designed a magnetization weight prior constraint based on the Laplacian matrix, which aims to accurately characterize the interaction and dependence between adjacent voxel points, so as to maximize the preservation and restoration of the morphological characteristics and detail information of the reconstructed image while ensuring sparsity.

[0160] R(C) can be further expressed as follows after adding the Laplace prior constraint:

[0161]

[0162] Where L = DW is the Laplacian matrix. W is the adjacency matrix, which represents the intensity relationship between each voxel point and every other point. D is the weighted degree matrix, where each element D ii Represents the total connection weight between voxel point i and other points.

[0163] The specific implementation of step S103 is:

[0164] Based on the sparse prior constraints and spatial prior constraints, we can construct the required objective function as follows:

[0165]

[0166] in, is the optimal approximate solution of the concentration vector C, α and ω are regularization parameters, ||C||1 is the L1 regularization sparse constraint term, is the Laplace space constraint term.

[0167] The specific implementation of step S104 is:

[0168] As mentioned above, based on non-convex optimization and iterative approximation theory, the objective function is iteratively optimized, and after multiple iterations, the optimal value of the particle information updated in the iterative process can be found. The k+1th iteration formula can be expressed as:

[0169] C (k+1) =C(k) -[S T (SC (k) -U)+α·sign(C (k) )+2ωL T LC (k) ]

[0170] Next, the core construction of the entire deep expansion reconstruction framework is based on an iterative formula, which guides the advancement of the entire reconstruction process. The expanded deep neural network model can be formally expressed as:

[0171] C (k+1) =Ψ θ (C (k) ; U, S, L)

[0172] Among them C (k) represents the concentration estimate of the kth iteration, and θ is the set of model parameters shared in all iteration steps. This model takes the Laplace matrix L, the complex measured voltage U and the system matrix s as input information, and the particle concentration C is initialized to all zeros in the initial stage to provide a benchmark starting point for the subsequent iteration process.

[0173] First, the iterative optimization process is explicitly expressed as the structure of a deep neural network, so that the iterative formula is optimized using the basic idea of ​​gradient descent. In the kth iteration, we calculate the fidelity term relative to the current solution term (denoted as C (k) ) in Represents the fidelity loss function. Subsequently, we move the solution term along the negative gradient direction of the fidelity term by a small step size η (learning rate), which can be expressed as:

[0174]

[0175] In this process, the results of the iteration are represented by complex numbers. Secondly, in the channel and dimension conversion layer, we convert the one-dimensional data into two-dimensional data and extract the real part and the imaginary part into two different channels respectively. This provides a more intuitive and easier to process data format for subsequent network layers. The conversion process can be expressed as:

[0176]

[0177] Third, in order to effectively suppress noise and redundant information while ensuring the quality of reconstruction, the near-end projection layer integrates the encoder, decoder, and spatial attention module. The encoder part uses a batch normalization layer to accelerate convergence and stabilize gradients, and improves the complexity and depth of feature representation through stacked convolutional layers and ReLU activation function structures. The size of the convolution kernel is 3, and the number of output channels is 2, 4, 4, 8, 8, and 16, respectively. Assume that the input feature is The encoding process can be expressed as:

[0178] h0=BN(h)

[0179] h i =ReLU(Conv(h i-1 )),i=1,2,…,6

[0180] Among them, h i is the feature representation after the i-th convolution layer. Next, we introduce the spatial attention module to add weights to the particle features at each position to further enhance the feature representation.

[0181] g6=F SAM (h6)

[0182] Among them, F SAM (·) represents the spatial attention calculation. Then, the decoder gradually reduces the high-level enhanced feature representation through the six-layer convolution operation with the same structure as the encoder, setting the number of reverse mapping channels to 16, 8, 8, 4, 4, 2, which can be expressed as:

[0183] g i =ReLU(Conv(g i+1 )),i=5,4,…,0

[0184] Among them, g i is the feature representation after the i-th convolution. Fourth, in order to further improve the accuracy of the prediction and balance the tendency of the predicted value to approach the target, we convert the output of the decoder into a complex form through the back channel and dimension conversion layer, and update the predicted value using the residual connection strategy according to the iterative formula. This step helps the network retain the memory of the initial value during the iteration process and gradually improve the accuracy of the prediction. The number of iterations is preset to k = 5. Finally, after the entire iteration is completed, the BN layer is used to remove the internal covariate bias to stabilize and accelerate the training process, and through two convolution operations, the nonlinear representation ability of the network is further improved to make the prediction result close to the target result.

[0185] Step S105 specifically includes:

[0186] (1) The weighted sum of mean squared error (MSE) and structural similarity index (SSIM) is used as the loss function for network model training;

[0187] (2) The Adam optimizer is used for optimization, and the learning rate is set to 1x10 -4, we trained the network 500 times with a batch size of 64. The model was trained using the PyTorch framework;

[0188] (3) Model testing and display of results;

[0189] We use the weighted sum of mean squared error (MSE) and structural similarity index (SSIM) as the loss function for network model training. The total loss function can be expressed as:

[0190] Loss=λ·MSE+(1-λ)·(1-SSIM)

[0191] Among them, λ = 0.9 is a weight parameter used to balance the relative importance between MSE and SSIM. MSE is used to evaluate the quality of the reconstructed image. SSIM is used to evaluate the similarity between the reconstructed image and the original image. MSE can be expressed as:

[0192]

[0193] Among them, Y and are the original and network predicted MPI images, respectively, and m is the number of images. SSIM can be expressed as:

[0194]

[0195] Among them, Y and are the original and network predicted MPI images respectively. Y and are image Y and The mean of . and are image Y and The variance of . is the image Y and The covariance of . c1 and c2 are constants for stability, which we set to 0.01 and 0.03 for normalized images. Then, the trained deep neural network model is tested on the test dataset to verify the performance of the model in practical applications. At the same time, the error graph between the prediction results and the real image is shown.

[0196] Figure 3 is a specific network diagram of the present invention. Figure (a) is the main framework of the network model, Figure (b) is the channel and dimension conversion layer, and Figure (c) is the proximal projection layer. For the initial concentration C (0)Perform full zero initialization and use S, U, and L as network inputs. S is the system matrix, U measures voltage, and L is the Laplace matrix, which is a symmetric matrix of 784×784 size that is symmetric about the main diagonal.

[0197] Figure 4 The figure in the figure shows the spatial magnetization response generated by the center of the FOV and its adjacent voxel points during the driving field cycle. During the spatial encoding process of MPI, the particles at the FFR and the surrounding particles will show a high temporal and spatial correlation. When the gradient is not very large, the particles around the FFR will also produce a small amount of magnetization response due to the low external magnetic field intensity. The spatial magnetization response generated by the center of the FOV and its adjacent voxel points in one driving field cycle.

[0198] Figure 5 In order to compare the reconstruction results of the present invention with the original label results, four groups of results were used to compare with the label values. The average MSE was calculated to be 0.055 and the SSIM was 0.96. The advanced nature of the invention can be proved from both the evaluation index and the visualization results. It can be seen from the comparison chart that the invention performs well in morphological restoration and has profound invention significance.

Claims

1. A magnetic particle reconstruction method based on an improved magnetization model and driven by prior information, characterized in that: The following steps are involved: Step 1: Based on the basic theory of magnetic nanoparticle imaging and combined with the improved Jiles-Atherton magnetization model, the magnetization behavior of magnetic nanoparticles in MPI is accurately modeled to obtain the system matrix and the measurement voltage of the particle distribution; the system matrix is ​​a matrix that characterizes the relationship between the magnetic field and the magnetization intensity, and the measurement voltage of the particle distribution is the induced voltage caused by the change of the magnetization intensity; According to the above process, first, considering the complex interaction between magnetic nanoparticles and the external magnetic field, the total effective magnetic field H is defined. e for: H e =H+αM (1); Where H is the external magnetic field, α is the coupling parameter, and M is the magnetization intensity of the magnetic nanoparticle. Secondly, the Langevin function is used as the hysteresis-free magnetization function to convert the hysteresis-free magnetization M an The relationship with the effective magnetic field He is expressed as: in, is the Langevin parameter, M s is the saturation magnetization of the magnetic nanoparticles. In MPI, the reversible magnetization M rev Related to the domain wall bending, the irreversible magnetization M irr It is related to the displacement of domain walls. The magnetic nanoparticles are arranged along the external magnetic field, which leads to the displacement of domain walls. Therefore, the irreversible magnetization intensity M irr Dominant, modified Jiles-Atherton magnetization model M JA It can be written as: M JA =M rev +M irr =c(M an -M irr )+M irr (3); In order to simplify the model, we can assume that the reversible parameter c is zero. In this case, M JA =M irr , M JA The model can be written as an ordinary differential equation describing the magnetization intensity of the particle under an external magnetic field H: in, The parameters α and k can be determined by the coercive point (H = H e , M=0) and remanence point (H=0, M=M r ) determines that at the coercive point of δ = +1, M JA The model equation (4) is: Similarly, at the remanence point of δ = -1, (4) can be written as: By combining equations (5) and (6), we can determine the parameters α and k, M JA The relationship between and H can be determined as: Based on formula (7), the induced voltage at r of a receiving coil with a sensitivity of p(r) can be expressed as: Where, M(r,t)=c(r)·M JA (r, t), c(r) is the particle concentration at position r, and s(r, t) is defined as the integral kernel describing the mapping between particle concentration and voltage: Then formula (8) can be expressed as: u(t)=∫ Ω c(r)·s(r,t)dr (10) In order to perform numerical calculations, the time and space of formula (10) are discretized: Among them, s(r n ,t k ) is calculated as: Where m represents the modulus of the magnetic moment of a single particle; Performing fast Fourier transform (FFT) on both sides of formula (11) in the time domain yields: Among them, S(f k ,r n ) is defined as the system matrix, It is defined as the measured voltage signal spectrum vector of MNPs, k = 1, 2, 3..., K is the number of frequency components, and n = 1, 2, 3..., n is the number of voxel points on the imaging grid; Step 2: Considering the unique properties of the magnetic nanoparticle solution space and the intrinsic correlation of magnetization between adjacent voxels in space, two prior information are used. One is the L1 regularization prior, which aims to promote the sparsity of the solution, effectively suppress noise interference, and improve the quality of the reconstructed image. The other is the Laplacian regularization prior of the magnetization weight, which further enhances the local coherence and spatial consistency of the reconstruction result by describing the magnetization gradient relationship in the space of adjacent particles. Based on the above content, formula (13) is described in matrix form as: Among them, S∈R K×N is the system matrix, C∈R N×1 is the particle concentration distribution, U∈R K×1 is the measured signal spectrum vector of MNPs, K is the number of frequency components, and N is the number of voxel points on the imaging grid; Since K is much smaller than N, the problem of solving the particle concentration distribution in formula (14) can be converted into solving an ill-posed inverse problem, which can be solved by the least squares method plus regularization constraints, and can be expressed as: in, In order to ensure that the reconstruction result C is as close to the true value as possible, R(C) is a regularization constraint term. By introducing additional constraints on the fidelity term of the inverse problem, the pathological nature of the inverse problem can be overcome, thereby obtaining a more stable and reliable solution. To ensure accurate reconstruction of the spatial distribution of magnetic nanoparticles from limited magnetization response data, the L1 regularization term is used as a priori constraint condition in view of the sparse distribution of magnetic nanoparticles in the field of view (FOV). The L1 regularization term promotes the sparsity of the solution, that is, reduces the number of non-zero elements, which helps to suppress the interference of noise during the reconstruction process, thereby improving the quality and accuracy of the reconstructed image. For formula (15), the regularization constraint term R(C) can be expressed as follows: In order to balance sparsity and reconstruction quality, a strategy combining spatial prior information is proposed. The magnetization weight prior constraint term based on the Laplacian matrix is ​​designed to accurately characterize the interaction and dependence between adjacent voxel points, so as to maximize the preservation and restoration of the morphological features and detail information of the reconstructed image while ensuring sparsity. The Laplacian matrix is ​​represented by L, where L = DW, W is the adjacency matrix, which represents the intensity relationship between each voxel point and each other point, and D is the weighted degree matrix, where each element D ii Represents the total connection weight between voxel point i and other points. Through this definition, the Laplacian matrix can effectively characterize the local connectivity and differences between voxel points in the image; Step 3: Based on the prior information described in step 2, a double prior constraint objective function is constructed. This objective function balances the relationship between data fidelity and the sparsity and smoothness of the solution. The data fidelity is to ensure the consistency between the reconstruction result and the original measurement data. The sparsity constraint aims to improve the effective information content of the reconstructed image by reducing the interference of unnecessary information; the smoothness constraint helps to eliminate the influence of noise; based on the sparse prior constraint term and the spatial prior constraint term of formula (27), we can write formula (15) as follows: in, is the optimal approximate solution of the concentration vector C, λ1 and λ2 are regularization parameters, ||C||1 is the L1 regularization sparse constraint term, is the Laplace space constraint term; Step 4: In order to obtain the objective function in step 3 Perform efficient iterative solution, build a deep unfolding model based on non-convex optimization and iterative approximation theory, and achieve high-fidelity reconstruction of the concentration distribution of magnetic nanoparticles. This is done in four steps: First, perform gradient calculation: For the first term of formula (28) We can use the chain rule to differentiate C and expand it to: The result of the derivation of C can be expressed as: The second term α||C||1, its partial derivative is non-smooth, and the sub-gradient of the L1 norm can be expressed as: Where sign(C) is the sign function, which means that the derivative of each element is +1 (if the element is positive), -1 (if the element is negative), and 0 (if the element is 0). Item 3 Expand this item to: ω(LC) T (LC)=ωC T L T LC (32); The result of the derivation of C can be expressed as: Combining the partial derivatives, the overall gradient formula can be expressed as: Based on non-convex optimization and iterative approximation theory, the optimal value of the particle information updated in the iterative process can be found after multiple iterations. The k+1th iteration can be expressed as: C (k+1) =C (k) -[S T (SC (k) -U)+α·sign(C (k) )+2ωL T LC (k) ] (35); Next, the core construction of the entire deep unfolding reconstruction framework is based on the iterative formula (35), which transforms the iterative solution process into a hierarchical structure of a deep network. Each iteration is mapped to a layer in the deep network. The layer-by-layer iterative approximation guides the advancement of the entire reconstruction process. The expanded deep neural network model can be formally expressed as: C(k+1)=Ψ θ (C(k);U,S,L)(36); Where C(k) represents the concentration estimate of the kth iteration, and θ is the set of model parameters shared in all iteration steps. This model takes the Laplace matrix L, the complex measured voltage U and the system matrix S as input information. The particle concentration C is initialized to all zeros in the initial stage to provide a benchmark starting point for the subsequent iteration process. Second, build a deep expansion model: In each iteration step, the gradient descent method is used to update the concentration estimate C(k). By calculating the negative gradient direction of the fidelity term, the result of each iteration is guided to move closer to the true solution. Given that the result generated by formula (34) is in complex form, a channel and dimension conversion layer is used. This layer is responsible for converting the one-dimensional complex output generated by each iteration into a two-dimensional data format, and cleverly separates the real and imaginary parts, and assigns independent channels for processing. Third, feature extraction and enhancement: The near-end projection layer consists of an encoder, a decoder, and a spatial attention module. The encoder and decoder, as feature extractors, can effectively capture and reconstruct key feature representations in the image. The spatial attention module further enhances the model's sensitivity to core feature information. By dynamically adjusting the weights of different spatial positions, the model can focus on the areas that are most critical to the reconstruction results, thereby suppressing irrelevant information while improving reconstruction accuracy. In order to effectively suppress noise and redundant information while ensuring the reconstruction quality, the proximal projection layer integrates the encoder, decoder and spatial attention module. The encoder part uses a batch normalization layer to accelerate convergence and stabilize the gradient, and improves the complexity and depth of feature representation through stacked convolutional layers and ReLU activation function structures. The size of the convolution kernel is 3, and the number of output channels is 2, 4, 4, 8, 8, and 16, respectively. Assume that the input feature is The encoding process can be expressed as: Among them, h i is the feature representation after the i-th convolution layer. Next, we introduce the spatial attention module to add weights to the particle features at each position to further enhance the feature representation; g6=F SAM (h6) (40); Among them, F SAM (·) represents the spatial attention calculation. Then, the decoder gradually reduces the high-level enhanced feature representation through the six-layer convolution operation with the same structure as the encoder. The number of reverse mapping channels is set to 16, 8, 8, 4, 4, 2, which can be expressed as: Among them, gi is the feature representation after the i-th convolution; Fourth, improve the accuracy of prediction: The tendency of the prediction value to approach the target is balanced. The output of the decoder is converted into a complex form through the back channel and the dimension conversion layer, and the prediction value is updated using the residual connection strategy according to the iterative formula (35). The number of iterations is preset to k = 5. After the entire iteration is completed, the internal covariate bias is removed through the BN layer to stabilize and accelerate the training process. Through two convolution operations, the nonlinear representation ability of the network is further improved to make the prediction result close to the target result. Step 5: In order to verify the practical application effect of the proposed method, the model is trained through the loss function and predicted using the test data set. At the same time, the predicted image and the real image are displayed; the weighted sum of the mean square error (MSE, Mean Squared Error) and the structural similarity index (SSIM, StructuralSimilarityIndex) is used as the loss function for network model training. The total loss function can be expressed as: Loss=λ·MSE+(1-λ)·(1-SSIM) (46); Among them, λ = 0.9 is a weight parameter used to balance the relative importance between MSE and SSIM. MSE is used to evaluate the quality of the reconstructed image, and SSIM is used to evaluate the similarity between the reconstructed image and the original image. MSE can be expressed as: Among them, Y and are the original and network predicted MPI images, respectively, and m is the number of images. SSIM can be expressed as: Among them, Y and are the original and network predicted MPI images respectively; μ Y and are image Y and The mean of and are image Y and The variance of is the image Y and covariance; c1 and c2 are constants for stability, which we set to 0.01 and 0.03 for the normalized image; then, the trained deep neural network model is used to test the test data set to verify the performance of the model in practical applications. At the same time, the error map between the prediction results and the real image is shown.

2. The magnetic particle reconstruction method based on an improved magnetization model and prior information driving according to claim 1, characterized in that: Construction of the Laplace matrix L = DW: First, the particles reach magnetization saturation when the external magnetic field applied in MPI makes the magnetization intensity of the particles reach 80% of the saturation magnetization intensity, and the modulus of the external magnetic field intensity H is defined as the saturation field intensity H. sat After the selection field and the driving field are superimposed, the particle magnetic field at the FFR is H(r, t) = H D (t)+H s (r) = 0, then the FFR surrounding particles r Δ The combined magnetic field can be expressed as: Where Δr is the displacement vector of the coordinate position around the FFR relative to the coordinate position of the FFR. Δ ,t)|| <H sat When , the spatial range of all components is defined as the δ neighborhood of FFR. For the relationship weight between each voxel point and FFR in the neighborhood, M JA The derivative of the model is used to measure the magnetization response at each coordinate position r in space, because formula (12) shows that the magnetization response generated at each coordinate position r in space is proportional to the rate of change of the magnetization intensity M(r, t) with respect to time, where M(r, t) = c(r)· Therefore, when the FFR position is determined, the spatial encoding time t FFR is also a fixed value, M JA The derivative of (β||H(r,t)||) with respect to H(r,t) reflects the intensity of the magnetization response of the entire space at that moment. The adjacency matrix W can be expressed as: Among them, i and j represent each voxel point, δ(i) represents the neighborhood of point i when it is the FFR point, and Δr ij represents the displacement vector from point j to i in the neighborhood, Indicates that H = G·Δr ij Time JA The derivative of (β||H||) with respect to H, Indicates that when point i is the FFR point, M JA The derivative of (β||H||) with respect to H is β / 3. For the determination of the neighborhood of δ(i), if the modulus of the applied magnetic field intensity at voxel point i is less than H sat , then the voxel point is in the neighborhood δ(i) of FFR point i, that is, by solving 0.8M sat =cmM JA (β||H||) Solve for H sat， Since the coordinates of voxel points are discrete during numerical calculation, to determine whether point j is in the neighborhood of i, it is only necessary to judge ||G·Δr based on the selection field gradient G ij || Is it less than H sat That’s it; for a selection field with good linearity, the number and distribution of voxel points in the magnetization space neighborhood of each voxel point are the same, so the calculation of the relative weight in the neighborhood of each voxel point can be simplified to only calculate the weight of each point in the neighborhood of a voxel point in space, and then match all the weights of the neighborhood of the point to other magnetization space neighborhood points to the voxel point according to the relative position in space, thereby completing the construction of the adjacency matrix W; The weighted degree matrix D can be expressed as: According to formula (18) and formula (19), the element L of the Laplace matrix L is ij It can be expressed as: L ij =D ij -W ij (20); According to the definition of D and W, L ij The value of can be divided into three cases: 1) When i=j: Note that W ii =0, therefore, For clarity, the index in the summation was changed from j to k to avoid confusion with the external j, but in practice they are identical. 2) When i≠j and i and j are adjacent (j∈δ(i)): L ij =D ij -W ij =0-W ij =-W ij (23); therefore, 3) When i≠j and i and j are not adjacent L ij =D ij -W ij =0-0=0 (25) Based on the above, the expression of the Laplace matrix L is: For n individual voxel points, L is a symmetric matrix of size n×n symmetric about the main diagonal. In summary, for formula (15), R(C) can be further expressed as follows after adding the Laplace prior constraint:

3. The magnetic particle reconstruction method based on an improved magnetization model and prior information driving according to claim 1, characterized in that: The calculation process of the model in step 4 is as follows: First, the iterative optimization process is explicitly expressed as the structure of a deep neural network, so that formula (35) is optimized using gradient descent. In the kth iteration, the gradient of the fidelity term relative to the current solution term is calculated. The previous solution is denoted as C(k), where Denotes the fidelity loss function. Subsequently, we move the solution term along the negative gradient direction of the fidelity term by a small step size η (learning rate), which can be expressed as: In this process, the results of the iteration are represented by complex numbers. Secondly, in the channel and dimension conversion layer, the one-dimensional data is converted into two-dimensional data, and the real part and the imaginary part are extracted into two different channels respectively. The conversion process can be expressed as:

4. The magnetic particle reconstruction method based on an improved magnetization model and prior information driving according to claim 1, characterized in that: The spatial attention module F SAM (·) receives a feature map X as input, whose size is C×H×W, where C is the number of channels, H and W are the height and width of the feature map respectively. The feature map is The output is the result of The module re-weights the different positions of the input feature map by calculating the weight map in space, thereby enhancing the focus on important areas, and uses global average pooling and global maximum pooling to compress the spatial information of each channel into a global scalar: Among them, f avg and f max The two vectors are transformed by the convolution layer to extract more complex features. The transformed vector is recorded as f avg and f avg ; Among them, W avg and W max is the weight matrix, b avg and b max is the bias term, and σ is the ReLU nonlinear activation function. Then, the transformed vectors are concatenated element-wise and passed through an additional convolution and a Sigmoid activation function to generate a spatial attention weight map M with a size of 1×H×W; Finally, the attention weight map M is element-wise multiplied with the original input feature map X to obtain the weighted feature map X′, X′(c,i,j)=X(c,i,j)×M(i,j)(45).

Images