CN108010098B - Double-energy-spectrum CT (computed tomography) base material image iterative reconstruction method - Google Patents

Abstract

The invention discloses a double-energy-spectrum CT (computed tomography) base material image iterative reconstruction method which is used for reconstructing a base material density image of a measured object. The method comprises the following steps: scanning a measured object by using a dual-energy-spectrum CT system to obtain dual-energy-spectrum projection data of the measured object; based on the high-energy spectrum projection data and the low-energy spectrum projection data, respectively reconstructing a high-energy spectrum image and a low-energy spectrum image of the measured object; calculating a linear combination image of the images of the measured object with high and low energy spectrums; and solving a ratio image of the basis material density image and the linear combination image, establishing a regularization constrained image reconstruction model by taking the total variation minimization of the ratio image as a constraint condition, and iteratively reconstructing the basis material density image of the measured object. The invention fully utilizes the structural consistency among various reconstructed images of the dual-energy CT, effectively reduces the reconstruction noise amplified due to the decomposition of the base material and improves the convergence rate of the iterative reconstruction of the base material image.

Description

An Iterative Reconstruction Method of Dual-energy Spectral CT-Based Material Image

technical field

The invention relates to the technical field of X-ray computed tomography (CT), in particular to a method for iterative reconstruction of images of dual-energy spectrum CT-based materials.

Background technique

X-ray CT technology is a technology that uses the principle of interaction between X-rays and matter to image the internal information of objects. It is not only used in medical testing, but also widely used in non-destructive testing of key components in various industries. Traditional single-energy-spectrum X-ray CT scans the measured object using one X-ray energy spectrum. The CT reconstruction images of different substances may have the same or similar CT values, and it is difficult to distinguish different substances. Dual-energy spectral CT scans the measured object using two X-ray energy spectra, which can obtain more information about the measured object than traditional single-energy spectral CT. Using this information, the density image or double-effect image of a specific base material can be reconstructed, and then the equivalent atomic number and electron density image of the measured object can be calculated, and substances can be distinguished based on this information. Dual-energy spectral CT not only has better material discrimination ability, but also effectively removes hardening artifacts in reconstructed images.

Since the conventional X-ray source emits polychromatic X-rays that obey a certain spectral distribution, rather than mono-energy X-rays, the dual-energy spectrum CT scans the measured object to collect polychromatic projections of high and low energy spectra, not ideal high- and low-energy monochrome projections. From a mathematical point of view, the CT image reconstruction problem based on dual-energy spectral polychromatic projection is a nonlinear inverse problem, which is characterized by ill-posedness and high dimensionality. At present, the decomposition methods of dual-energy spectral CT-based materials can be roughly divided into two categories: projection domain decomposition method and image domain decomposition method. The projection domain and image domain decomposition methods can be further divided into direct methods and iterative methods, respectively. Due to the huge difference in the attenuation coefficients of different base materials in the dual-energy decomposition process, the noise of the reconstructed image of dual-energy CT is amplified, especially the noise of the image of the base material with low equivalent atomic number is seriously amplified. Although the steps of the direct decomposition method are simple, since there is no feedback correction mechanism, there are strong artifacts in the reconstructed images. The iterative decomposition method uses numerical methods or optimization methods to construct an iterative structure, and obtains information such as multi-base material density images by gradually correcting the image reconstruction results. Although this method can reconstruct images of better quality, it has a slow convergence speed and a large amount of computation, which makes it difficult to use in actual dual-energy CT systems. Using prior knowledge as a regularization constraint is an important means for iterative methods to effectively reduce reconstruction noise, improve reconstruction quality, and speed up convergence. Therefore, how to discover more effective prior knowledge of image reconstruction, especially the prior knowledge with the characteristics of dual-energy spectral CT imaging, for iterative reconstruction of dual-energy spectral images, is an important research content.

SUMMARY OF THE INVENTION

The present invention provides an iterative reconstruction method of a dual-energy spectrum CT base material density image, which is used to solve the problem that the decomposition noise is amplified in the base material image in the prior art, so as to improve the reconstruction quality of the base material image. The method of the present invention comprises the following steps:

S1: Scan the measured object with the X-ray dual-energy spectrum CT system, and collect the high- and low-energy spectrum multicolor projection data of the measured object;

S2: According to the collected high- and low-energy spectral polychromatic projection data, the traditional single-energy CT reconstruction method is used to directly reconstruct the high- and low-energy spectral images;

S3: Combine the high- and low-energy spectral images into a reference image, and decompose it into a 0-value part and a positive-value part;

S4: Find the ratio image of the base material density image and the linear combination image according to the positive value part of the reference image, and use the minimum total variation of the ratio image as a constraint to establish a regularized constrained image reconstruction model, and iteratively reconstruct the measured object's image. Base material density image.

Further, wherein the step S4 includes:

S41: Assign an initial value of 0 to the density image of the dual-base material of the object to be reconstructed;

S42: Iteratively reconstruct the dual-base material density image based on dual-energy spectrum multicolor projection, and update only the pixel value corresponding to the positive value part of the reference image in each iteration;

S43: Calculate the ratio image from the positive value part of the reference image and the corresponding base material density image;

S44: Make regularization constraints on the current dual-base material density image based on the ratio image, and correct the dual-base material density image;

S45: Repeat steps S42-S44 until the iteration end condition is satisfied, and reconstruct the density image of the dual-base material.

Further, the calculation formula of the reference image u _ref described in step S3 is as follows:

Here R ^-1 represents the inverse Radon transform, and R ^-1 (q _k ) is the image of the measured object directly reconstructed by the traditional single-energy CT reconstruction method based on the polychromatic projection q _k of the k-th energy spectrum, abbreviated as "energy spectrum". image".

Further, the objective function of the regularization constraint-based dual-energy spectrum CT-based material image decomposition in step S4 can be expressed as follows:

是相应的由正投影模型计算得到的多色投影值；w_k,L是方程的权值，由该射线测得的多色投影数据的方差决定；Ω_k是对应于第k个能谱的射线路径的集合；

是比值图像；α_m是用于调节各项平衡的权重因子。where f _m is the m-th base material density image to be reconstructed; q _k,L is the polychromatic projection data of the k-th energy spectrum collected along the ray L;

is the corresponding polychromatic projection value calculated by the orthographic projection model; w _k,L is the weight of the equation, determined by the variance of the polychromatic projection data measured by the ray; Ω _k is corresponding to the kth energy spectrum a collection of ray paths;

is the ratio image; α _m is the weighting factor used to adjust the balance of each item.

Further, the discretization formula of the dual-energy spectral orthographic model is:

Among them, S _k,E is the effective energy range of the kth energy spectrum divided into N _k intervals at equal intervals, and the length of each interval is the sampling value in the Eth interval; ψ _m,E is the measured object The mass attenuation coefficient _{of the mth} base material _with respect to the X _- photon energy E in The contribution weight of the jth pixel of _fm to the projection along ray L.

的定义如下：Further, the ratio image

is defined as follows:

为步骤S3中得到的参考图像u_ref的正值部分；

为第m种基材料密度图像中与

对应的部分像素。in

is the positive value part of the reference image u _ref obtained in step S3;

is the mth base material density in the image with

corresponding partial pixels.

Further, the objective function of image decomposition of dual-energy spectral CT-based material is divided into the following two sub-problems to solve:

表示其在第n次迭代时得到的解。得到

后，计算第n次迭代时的比值图像

然后用软阈值等方法求解第二个子问题，再乘以参考图像恢复到基材料密度图像域，完成了对基材料密度图像的正则化，得到第n次迭代的基材料密度图像

在计算中，对重建的图像加像素值非负约束。The first sub-problem can be solved by weighted extended algebraic iteration method, penalized likelihood method, penalized weighted least squares method, etc.

represents its solution at the nth iteration. get

After that, calculate the ratio image at the nth iteration

Then, the second sub-problem is solved by methods such as soft threshold, and then multiplied by the reference image to restore to the base material density image domain, completes the regularization of the base material density image, and obtains the base material density image of the nth iteration.

In the calculation, a pixel value non-negative constraint is imposed on the reconstructed image.

Further, the reconstruction process of the double-base material density image of the measured object is calculated by a parallel method, and is accelerated based on a hardware parallel computing platform.

To sum up, compared with the prior art, the method of the present invention has the following beneficial effects:

(1) Make full use of the structural consistency between various reconstructed images of dual-energy CT, effectively reduce the reconstruction noise amplified by the decomposition of the base material, and improve the convergence speed of the iterative reconstruction of the base material image;

(2) In each iteration calculation, only the non-zero pixels in the base material image are updated, which reduces the amount of calculation and speeds up the calculation speed of each iteration;

(3) It has strong applicability and can be applied to multi-spectral CT systems.

Description of drawings

In order to explain the embodiments of the present invention or the technical solutions in the prior art more clearly, the following briefly introduces the accompanying drawings that need to be used in the description of the embodiments or the prior art. Obviously, the accompanying drawings in the following description are only These are some embodiments of the present invention. For those of ordinary skill in the art, other drawings can also be obtained according to these drawings without creative efforts.

1 is a flow chart of a method for decomposing a dual-energy spectrum CT-based material according to an embodiment of the present invention;

Figure 2 is a central slice image of a FORBILD head model used as a test phantom;

Figure 3 is a schematic diagram of the high and low energy spectra of X-rays used in the experiment;

4a is a high-energy spectrum multicolor projection data image of a FORBILD head model collected under a tube voltage of 140kV according to an embodiment of the present invention;

4b is a low-energy spectrum multicolor projection data image of a FORBILD head model collected under a tube voltage of 80kV according to an embodiment of the present invention;

Fig. 5a is a high-energy-spectrum phantom image directly reconstructed by high-energy-spectrum polychromatic projection according to an embodiment of the present invention;

Fig. 5b is a low-energy-spectrum phantom image directly reconstructed by low-energy-spectrum polychromatic projection according to an embodiment of the present invention;

Fig. 5c is a reference image formed by combining a high-energy-spectrum measured object image and a low-energy-spectrum measured object image according to an embodiment of the present invention;

Fig. 6a is a cross-sectional view of line 287 of the high-energy spectrum phantom image in Fig. 5a according to an embodiment of the present invention;

FIG. 6b is a cross-sectional view of the 287th row of the low-energy spectrum phantom image in FIG. 5b according to an embodiment of the present invention;

FIG. 6c is a cross-sectional view of line 287 of the reference image in FIG. 5c according to an embodiment of the present invention;

Figure 7a is a result image of a bone-based material density image reconstructed by using the E-ART+TV method for 100 iterations according to an embodiment of the present invention, and the displayed grayscale window is [0.9, 1.1];

Fig. 7b is a result image of a water-based material density image reconstructed by using the E-ART+TV method for 100 iterations according to an embodiment of the present invention, and the displayed grayscale window is [0.95, 1.15];

Fig. 7c is a linear attenuation coefficient image of photons with energy of 70keV synthesized from the bone-based image of Fig. 7a and the water-based image of Fig. 7b, and the displayed grayscale window is [-50HU, 150HU];

Fig. 8a is a result image of a bone-based material density image reconstructed by the method of the present invention iteratively 100 times according to an embodiment of the present invention, and the displayed grayscale window is [0.9, 1.1];

Fig. 8b is a result image of a water-based material density image reconstructed by the method of the present invention iteratively 100 times according to an embodiment of the present invention, and the displayed grayscale window is [0.95, 1.15];

Figure 8c is a linear attenuation coefficient image of photons with energy of 70 keV synthesized from the bone-based image of Figure 8a and the water-based image of Figure 8b, showing a grayscale window of [-50HU, 150HU].

Detailed ways

The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, but not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative efforts shall fall within the protection scope of the present invention.

Conventional X-ray sources emit polychromatic X-rays that obey a certain spectral distribution, not single-energy X-rays. Therefore, dual-energy spectrum CT scans the measured object and collects not ideal high- and low-energy monochromatic projections, but high-energy X-rays. , Polychromatic projection of low energy spectrum. In the case of ignoring the X-ray scattering effect, the base material decomposition problem of dual-energy spectral CT can be expressed as the following mathematical model:

表示沿着给定X射线路径L由多色投影方程计算得到的多色投影值；Ω_k是对应于第k个能谱的X射线路径的集合。双能谱CT的基材料分解问题就是利用测得的高、低能谱在各个角度下的多色投影q_k,L(k＝1,2；L∈Ω_k)，反求被测物体的双基材料密度图像f_m(x)(m＝1,2)。Among them, the parameter k represents the energy spectrum index; _Sk (E) represents the k-th normalized X-ray energy spectrum; ψ _m (E) is the mass of the m-th base material in the measured object with respect to the X-photon energy E Attenuation coefficient; f _m (x) represents the density spatial distribution of this base material, that is, the base material density image;

represents the polychromatic projection value calculated by the polychromatic projection equation along a given X-ray path L; Ω _k is the set of X-ray paths corresponding to the kth energy spectrum. The base material decomposition problem of dual-energy spectral CT is to use the measured polychromatic projections q _k,L (k=1,2; L∈Ω _k ) of the measured high and low energy spectra at various angles to reversely find the dual Base material density image f _m (x) (m=1,2).

In numerical implementation, the reconstructed base material density image is discretized into pixels. Let f _m =(f _m,0 ,f _m,1 ,...,f _m,J-1 ) ^τ denotes the discrete form of the image f _m (x), where f _m,j is the jth in f _m (x) The sampled values on pixels, the parameter J is the total number of pixels in each image, and τ represents the transpose of the vector. Let A _L =(a _L,0 ,a _L,1 ,...,a _L,J-1 ) be the system projection vector, where a _L,j denotes the projection of the jth pixel pair of f _m along the ray L Contribution weight. The effective energy range of the k-th energy spectrum is divided into N _k intervals at equal intervals, and the length of each interval is δ. S _k (E) and ψ _m (E) are sampled in each interval to obtain their discrete The forms are:

Using the above notation definitions, the discrete form of the polychromatic projection can be expressed as follows:

In this way, the dual-energy spectral CT reconstruction problem is transformed into the measured polychromatic projection q _k,L (k=1,2; L∈Ω _k ). By solving the nonlinear equation system of (3), f _m (m= 1,2).

In order to obtain a high-quality base material density image, the method of the present invention reconstructs the base material density image of the measured object by using a regularization-constrained iterative reconstruction method. The specific method is as follows:

The above formula contains two terms, the first term is the data fidelity term, and the second term is the image regularization constraint term. The data fidelity term makes the reconstructed image satisfy the polychromatic projection transformation as much as possible by minimizing the square of the weighted polychromatic projection residual. Here w _k,L is the weight of the equation, determined by the variance of the multicolor projection data measured by this ray.

为比值图像，定义为：The parameter α _m in the second item is a weight factor used to adjust the balance of each item, which can be gradually reduced with iteration;

is a ratio image, defined as:

为参考图像u_ref中大于0的部分；

为第m种基材料密度图像中与

对应的部分像素；u_ref具体定义为：here

is the part greater than 0 in the reference image u _ref ;

is the mth base material density in the image with

The corresponding part of the pixel; u _ref is specifically defined as:

Here R ^-1 represents the inverse Radon transform, and R ^-1 (q _k ) is the image of the measured object directly reconstructed by the traditional single-energy CT reconstruction method based on the polychromatic projection q _k , referred to as "energy spectrum image" for short. The reconstructed image based on high-energy spectral polychromatic projection data is called "high-energy spectral image"; the reconstructed image based on low-energy spectral polychromatic projection data is called "low-energy spectral image". Since pleochroism is not considered, hardening artifacts are present in spectral images, especially in low-energy spectral images. The reference image u _ref is obtained by combining the high-energy spectrum image and the low-energy spectrum image. The structure of the measured object in the figure remains unchanged, not only the hardening artifacts are greatly reduced, but also there are no artifacts caused by dual-energy decomposition. Here, the combination coefficient β _k can be pre-calibrated.

对应的

由于减少了未知数，加快了计算速度。The reference image has the property of being non-negative, that is, the pixel value is ≥ 0. The method of the present invention is based on the fact that where the pixel value of the energy spectrum image or the reference image of the measured object is 0, the pixel value of the corresponding base material density image must also be 0. Therefore, when calculating the density image of the base material, it is only necessary to calculate and reference the image

corresponding

Speeds up computation due to fewer unknowns.

有很强的梯度稀疏性。本文充分利用了比值图像的梯度稀疏性，通过最小化两个比值图像

和

的梯度的L₁范数之和，对基材料图像进行正则化，去除图像中的伪影和噪声。Since the structure of the base material image of the tested object and the reference image are consistent, its ratio image

There is strong gradient sparsity. This paper takes full advantage of the gradient sparsity of ratio images, by minimizing the two ratio images

and

The sum of the L1 norm _of the gradients of , regularizes the base material image to remove artifacts and noise in the image.

Formula (4) is divided into two sub-problems to solve:

是子问题(7)在第n次迭代时得到的解。得到

后，基于预先计算的参考图像u_ref利用(5)式计算第n次迭代时比值图像

软阈值处理后，再乘以参考图像恢复到基材料密度图像域，完成了对基材料密度图像的正则化。注意对重建的图像加像素值非负的约束。Sub-problem (7) can be solved by weighted extended algebraic iterative method, penalized likelihood method (PL) and other algorithms; sub-problem (8) can be solved by soft threshold algorithm and so on. Here in sub-problem (8)

is the solution of subproblem (7) at the nth iteration. get

Then, based on the pre-calculated reference image u _ref , the ratio image at the n-th iteration is calculated by formula (5).

After soft thresholding, it is multiplied by the reference image and restored to the base material density image domain, completing the regularization of the base material density image. Note that the reconstructed image imposes a non-negative constraint on pixel values.

The formula for solving subproblem (7) using the weighted extended algebraic iteration method is as follows:

为第n-1次迭代时计算得到的第m个能谱的多色投影沿着射线L对第m个基材料图像的修正残差，定义如下：Here λ is the relaxation factor;

is the correction residual of the mth base material image along the ray L of the polychromatic projection of the mth energy spectrum calculated at the n-1th iteration, and is defined as follows:

的计算公式分别如下：this parameter

The calculation formulas are as follows:

As an embodiment, the main implementation steps of the iterative reconstruction of the density image of the dual-base material are summarized as follows:

和正值部分

1) The high and low energy spectral images are directly reconstructed from the measured high and low energy spectral polychromatic projection data, respectively, and the reference image u _ref is combined based on the pre-calibrated combination coefficients, and then the 0-value part in the reference image is obtained.

and the positive part

设迭代次数n＝0；2) Assign initial values to the density image of the dual-base material to be reconstructed

Set the number of iterations n = 0;

注意仅更新与

对应的像素值；3) Solve formula (7) to obtain the estimated value of the n-th iteration of the density image of the dual-base material

Note that only updates with

the corresponding pixel value;

和

计算出第n次迭代的比值图像

4) by

and

Calculate the ratio image for the nth iteration

5) Solve formula (8) to obtain the regularized correction value of the n-th iteration of the double-base material density image

6) Repeat steps 3) to 5) until the iteration end condition is satisfied, and the density image of the dual-base material is reconstructed.

FIG. 1 is a flowchart of a method for iterative reconstruction of a dual-energy spectrum CT image according to an embodiment of the present invention. In order to improve the reconstruction speed, the reconstruction process of the double-base material density image of the object under test in the embodiment of FIG. 1 can also be calculated by a parallel method, and the implementation is accelerated based on a hardware parallel computing platform.

A specific implementation process of the dual-energy spectrum CT image iterative reconstruction method of the present invention will be described below through a specific embodiment. Figure 2 is a central slice image of the FORBILD head model used as a test phantom in this example, excluding the ears. Water and bone tissue were selected as the two base materials, the density of water was 1.0 g/cm ³ , and the density of bone tissue was 1.92 g/cm ³ . FIG. 3 shows the low and high energy spectrums of X-rays emitted by the dual-energy spectrum CT system in this embodiment. These two energy spectra are emitted when the tube voltages of the X-ray tube are 80kV and 140kV (1mm copper filter) respectively. In the numerical implementation, the sampling interval of the energy spectrum is set to 1keV. In dual energy spectral CT scanning, fan beam scanning is used. The scanning parameters of the system are: the distance from the ray source to the detector is 1200mm, the distance from the ray source to the center of the turntable is 1000mm, and the detector consists of 512 detection units of 0.6mm. Based on these parameters, the field of view is 256mm. The measured object is scanned 360° with high and low energy spectrum X-rays, and multi-color projection data is collected every 1°. The collected multicolor projection data of the test phantom is shown in Figure 4, wherein Figure 4a is the multicolor projection data collected under the high energy spectrum; Figure 4b is the multicolor projection data collected under the low energy spectrum. Based on these projection data, the algebraic iterative algorithm ART + total variation TV minimization algorithm is used to correct the high-energy spectral image and low-energy spectral image of the measured phantom directly after one iteration, as shown in Figure 5a and Figure 5b, respectively. The size of the image is 512×512. The reference image formed by combining these two images is shown in Fig. 5c, with the combination parameters β ₁ =8.62 and β ₂ =-1.50. Figures 6a, 6b, and 6c are cross-sectional views of line 287 of Figures 6a, 6b, and 6c, respectively. As can be seen from these figures, there are hardening artifacts in the images directly reconstructed by the polychromatic projection, especially in the low-energy spectrally reconstructed images. With the combination, hardening artifacts in the reference image are greatly reduced. In order to highlight the advantages of the method of the present invention, a comparison was made with the E-ART+TV method under the same conditions. Both the E-ART+TV method and the method of the present invention are implemented in parallel based on GPU. Figures 7a and 7b are respectively the result images of the bone-based material density image and the water-based material density image of the tested phantom reconstructed by the E-ART+TV method for 100 iterations. In order to clearly see the artifacts and noise in the reconstructed images, the grayscale windows displayed are set to [0.9, 1.1] and [0.95, 1.15], respectively. Fig. 7c is a linear attenuation coefficient image of photons with energy of 70 keV synthesized from the bone-based image of Fig. 7a and the water-based image of Fig. 7b, showing a grayscale window of [-50HU, 150HU]. It can be seen from these figures that although the E-ART+TV algorithm can effectively decompose the base material image, the convergence speed is slow, and there are still obvious artifacts and unsatisfactory decomposition in the images generated by 100 iterations. . Fig. 8a and Fig. 8b are respectively the result images of the bone-based material density image and the water-based material density image reconstructed by the method of the present invention after 100 iterations, and the displayed grayscale windows are also set to [0.9, 1.1] and [0.95, 1.15] . Figure 8c is a linear attenuation coefficient image of photons with energy of 70 keV synthesized from the bone-based image of Figure 8a and the water-based image of Figure 8b, showing a grayscale window of [-50HU, 150HU]. It can be seen from Figure 8 that after the same number of iterations as the E-ART+TV method, the artifacts in the base material density image reconstructed by the method of the present invention are effectively suppressed, and the reconstruction effect is significantly better than the E-ART+TV method. reconstruction effect.

Those of ordinary skill in the art can understand that the accompanying drawing is only a schematic diagram of an embodiment, and the modules or processes in the accompanying drawing are not necessarily necessary to implement the present invention.

Those skilled in the art may understand that: the modules in the apparatus in the embodiment may be distributed in the apparatus in the embodiment according to the description of the embodiment, and may also be located in one or more apparatuses different from this embodiment with corresponding changes. The modules in the foregoing embodiments may be combined into one module, or may be further split into multiple sub-modules.

Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, but not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that it can still be The technical solutions described in the foregoing embodiments are modified, or some technical features thereof are equivalently replaced; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the spirit and scope of the technical solutions in the embodiments of the present invention.

Claims (7)

1. a dual-energy spectrum CT base material image iterative reconstruction method, is characterized in that, comprises the following steps: S1: Scan the measured object with the X-ray dual-energy spectrum CT system, and collect the high- and low-energy spectrum multicolor projection data of the measured object; S2: According to the collected high- and low-energy spectral polychromatic projection data, the traditional single-energy CT reconstruction method is used to directly reconstruct the high- and low-energy spectral images of the measured object; S3: Combine the high and low energy spectrum images of the measured object into a reference image, and decompose it into a 0-value part and a positive-value part; S4: Find the ratio image of the base material density image and the linear combination image according to the positive value part of the reference image, and use the minimum total variation of the ratio image as a constraint to establish a regularized constrained image reconstruction model, and iteratively reconstruct the measured object's image. Base Material Density Image, The traditional single-energy CT reconstruction method described in step S2 includes an analytical algorithm and an iterative algorithm, the analytical algorithm includes a filtered back-projection algorithm and a back-projection filtering algorithm, and the iterative algorithm includes an additive iterative algorithm and a multiplicative iterative algorithm, Wherein step S4 includes: S41: Assign an initial value of 0 to the density image of the dual-base material of the object to be reconstructed; S42: Iteratively reconstruct the dual-base material density image based on dual-energy spectrum multicolor projection, and update only the pixel value corresponding to the positive value part of the reference image in each iteration; S43: Calculate the ratio image from the positive value part of the reference image and the corresponding base material density image; S44: Make regularization constraints on the current dual-base material density image based on the ratio image, and correct the dual-base material density image; S45: Repeat steps S42-S44 until the iteration end condition is satisfied, and reconstruct the density image of the dual-base material, The objective function of the image decomposition of the dual-energy spectrum CT-based material based on the regularization constraint in step S4 is expressed as follows:

where f _m is the m-th base material density image to be reconstructed; q _k,L is the polychromatic projection data of the k-th energy spectrum collected along the ray L;

is the corresponding polychromatic projection value calculated by the orthographic projection model; w _{k, L} is the weight of the equation, set as the reciprocal of the variance of the polychromatic projection data measured by the ray; Ω _k is corresponding to the kth energy a collection of spectral ray paths;

is the ratio image; α _m is the weight factor used to adjust the balance of each item, the specific value is determined by the image of the measured object, and it is set as the maximum reciprocal of the absolute value of the gradient of the base material image. 2 . The method for iterative reconstruction of an image of a dual-energy spectrum CT-based material according to claim 1 , wherein the X-ray dual-energy spectrum CT system described in step S1 comprises a high- and low-energy spectrum of a single-ray source and single-detector system. 3 . Any one of the two scanning modes, the simultaneous scanning mode of dual ray sources and dual detectors, the fast switching scanning mode of high and low energy spectrum, and the scanning mode of different energy spectrum segments based on photon detectors. 3. The iterative reconstruction method of dual-energy spectrum CT base material image according to claim 1, wherein the calculation formula of the reference image u _ref described in step S3 is as follows:

Here R ^-1 represents the inverse Radon transform, and R ^-1 (q _k ) is the image of the measured object directly reconstructed by the traditional single-energy CT reconstruction method based on the polychromatic projection q _k of the k-th energy spectrum, abbreviated as "energy spectrum". image"; β _k is a pre-scaled combination coefficient. 4. The iterative reconstruction method of dual-energy spectrum CT base material image according to claim 1, wherein, based on the objective function of regularization-constrained dual-energy spectrum CT base material image decomposition, the effective energy of the kth energy spectrum is The range is divided into N _k intervals at equal intervals, and the length of each interval is δ, then the discretization formula of the dual-energy spectrum orthographic projection model is:

Among them, E is the energy interval number. When sampling the energy spectrum at 1keV interval, E is a positive integer, 1≤E≤Nk; _{Sk , E} is the kth normalized X-ray energy spectrum at the Eth energy The average value of the energy in the interval; ψ _m,E is the average value of the mass attenuation coefficient of the mth base material in the measured object in the Eth energy interval; A _L =(a _L,0 ,a _L,1 , ...,a _L,J-1 ) is the system projection vector, a _L,j represents the contribution weight of the jth pixel of f _m to the projection along the ray L, the upper limit of the value of j is equal to the size of the current reconstructed image, The upper limit of the value of j is the square of 512 or the square of 1024. 5. The method for iterative reconstruction of dual-energy spectrum CT-based material images according to claim 1, wherein the objective function for image decomposition of dual-energy spectrum CT-based materials based on regularization constraints, the ratio image

is defined as follows:

in

is the positive value part of the reference image u _ref obtained in step S3;

is the mth base material density in the image with

corresponding partial pixels. 6. The method for iterative reconstruction of dual-energy spectrum CT base material images according to claim 1, wherein the objective function for image decomposition of dual-energy spectrum CT base materials based on regularization constraints is divided into the following two sub-problems Solve:

The first sub-problem can be solved by any one of the weighted extended algebraic iteration method, the penalized likelihood method, and the penalized weighted least squares method.

Indicates that the solution obtained at the nth iteration is obtained

After that, calculate the ratio image at the nth iteration

Then use the soft threshold method to solve the second sub-problem, and then multiply the reference image to restore to the base material density image domain, complete the regularization of the base material density image, and obtain the base material density image of the nth iteration

In the calculation, a pixel value non-negative constraint is imposed on the reconstructed image. 7 . The method for iterative reconstruction of dual-energy spectrum CT-based material images according to claim 1 , wherein the reconstruction process of the double-based material density image of the measured object is calculated by a parallel method, and is accelerated based on a hardware parallel computing platform. 8 .

Images