A method for non-linear reconstruction of high-energy flash x-ray images based on MCMC

Abstract

The application discloses a high-energy flash X-ray image nonlinear reconstruction method based on MCMC. According to the high-energy flash X-ray imaging principle, a discrete form of a nonlinear forward model is constructed, and a corresponding Jacobian matrix form is derived, the solution and uncertainty quantification of the inverse problem are considered in combination with the Bayesian theory, a hyperparameter based on weak information prior is introduced to construct a nonlinear hierarchical Bayesian model. By accelerating the solution of the optimization problem of random disturbance to sample the conditional distribution, the solution of the optimization problem is constrained in combination with the Jacobian matrix projection, and the proposal distribution of the target parameter is designed to reduce the sample statistical bias. Under the minimum variance criterion, the sample values of the linear and nonlinear Bayesian models are fused to obtain the final reconstructed image. The application can improve the sample estimation efficiency while ensuring that the reconstructed result presents clear edges and high precision.

Description

Technical Field

[0001] This invention relates to a nonlinear reconstruction method for high-energy flash X-ray images, belonging to the field of image processing technology. Background Technology

[0002] High-energy flash X-ray radiography, as a crucial tool for diagnosing the internal structure and density information of equipment under the context of nuclear testing bans, can solve a series of problems from equipment design to weaponization and engineering through laboratory simulation and numerical simulation, thereby ensuring the reliability and safety of the equipment. It has received high attention from the defense sectors of various countries. Unlike common medical X-ray imaging, high-energy flash X-rays irradiate targets with high material density and equivalent thickness, requiring high-intensity X-ray sources, typically in the MeV (million electron volts) range. Furthermore, under high-temperature and high-pressure detonation conditions, the radiographic system must be able to capture the internal structural information of targets with complex morphologies and rapidly changing speeds of up to several kilometers per second with a sufficiently narrow pulse width, exhibiting characteristics of high light emission speed and high safety. How to reconstruct high-precision target internal density information from low signal-to-noise ratio radiographic data is a key research focus in the quantitative diagnosis and analysis of high-energy flash X-ray radiography.

[0003] Due to the degradation effects of system ambiguity, noise, and scattering during high-energy flash X-ray radiography, its forward projection model exhibits nonlinear characteristics. Projected images often suffer from noise interference and detail blurring, affecting the accuracy of target internal density inversion and interface detection. To reduce the computational burden of reconstruction algorithms for high-dimensional data, most current high-energy flash X-ray image reconstruction work is conducted on simplified imaging models, such as linear approximation models. Linear reconstruction methods for X-ray images mainly include the Filtered Back Projection (FBP) algorithm, the Algebraic Reconstruction (ART) algorithm, and the Constrained Conjugate Gradient (CCG) algorithm. While these methods can solve the linear reconstruction problem in high-energy flash X-ray radiography, they cannot analyze the uncertainty of the reconstruction results and cannot handle nonlinear reconstruction problems involving system ambiguity in high-energy flash X-ray radiography. For the sake of convenience in model solving, linear reconstruction methods neglect the influence of system ambiguity on the reconstruction results. Their linear approximation reconstruction results are easily affected by ambiguity in the target interface region, impacting reconstruction accuracy and the accuracy of subsequent target interface detection. Therefore, it is necessary to construct a Bayesian reconstruction model based on a nonlinear model that better reflects the actual imaging process. Due to the need for more complex reasoning and higher computational costs, research on nonlinear reconstruction is relatively limited. Considering the nonlinear reconstruction of high-energy flash X-ray images from a Bayesian theoretical perspective, thereby estimating the target's linear absorption coefficient distribution and quantifying the uncertainty of target parameters, presents a significant challenge. Based on the nonlinear imaging model, this invention studies a Markov chain Monte Carlo (MCMC) reconstruction method for nonlinear reconstruction of high-energy flash X-ray images. This method effectively improves the efficiency and accuracy of nonlinear model sample estimation, ensuring reconstruction accuracy while also increasing the nonlinear reconstruction speed. Summary of the Invention

[0004] To address the issues that existing linear reconstruction results of high-energy flash X-ray images are susceptible to system fuzziness, and that nonlinear reconstruction algorithms require more complex reasoning and higher computational costs, this invention proposes a nonlinear reconstruction method for high-energy flash X-ray images based on MCMC. This method aims to improve the accuracy and efficiency of high-energy flash X-ray image reconstruction while providing uncertainty in the reconstruction results, thus ensuring the status assessment and operational maintenance of combat-ready equipment.

[0005] To address the aforementioned technical problems, this invention provides a nonlinear reconstruction method for high-energy flash X-ray images based on MCMC, comprising the following steps:

[0006] 1) Based on the principle of high-energy flash X-ray imaging, a discrete nonlinear forward model is constructed;

[0007] 2) Derive the Jacobian matrix form of the residual matrix of the nonlinear positive model in step (1), and optimize the solution by constructing a least squares model to obtain the deterministic reconstruction result;

[0008] 3) Combining Bayesian theory to consider the solution and uncertainty quantification of this nonlinear inverse problem, a nonlinear hierarchical Bayesian model is constructed by introducing hyperparameters based on weak information priors.

[0009] 4) By accelerating the solution of the optimization problem with random perturbation, the conditional distribution is sampled, and the solution of the optimization problem is constrained by the Jacobian matrix projection. The proposed distribution of the objective parameter is designed to reduce the statistical bias of the sample. The random perturbation strategy can ensure that the sampled samples are randomly generated in the high probability region of the posterior distribution, and then the confidence interval of the sample can be calculated.

[0010] 5) A multi-model fusion strategy is proposed, which fuses sample values ​​from linear and nonlinear Bayesian models under the minimum variance criterion to obtain the final reconstructed image.

[0011] The beneficial effects achieved by this invention are as follows: This invention solves the problem that the linear reconstruction results of high-energy flash X-ray images are affected by system ambiguity and noise, effectively improves the efficiency of sample estimation while ensuring that the reconstruction results present clear edges and high accuracy, and can achieve high-precision nonlinear reconstruction of high-energy flash X-ray images. At the same time, it provides uncertainty analysis of the reconstruction results, and better quantifies the internal information of weapons. Attached Figure Description

[0012] Figure 1 This is a flowchart of the present invention.

[0013] Figure 2 This is a schematic diagram of a high-energy flash X-ray imaging system.

[0014] Where 1 represents a pulse-driven injector, 2 represents injected electrons, 3 represents a light-sensing accelerator, 4 represents an electron beam, 5 represents a metal target, 6 represents a collimator, 7 represents X-rays, 8 represents a front protective device, 9 represents an imaging object, 10 represents a rear protective device, 11 represents a scintillator, 12 represents visible light, and 13 represents a CCD camera. Detailed Implementation

[0015] Schematic diagram of the present invention as follows Figure 1 As shown, the nonlinear reconstruction method for high-energy flash X-ray images based on MCMC of the present invention includes the following specific steps:

[0016] 1) Based on the principle of high-energy flash X-ray imaging, a discrete nonlinear forward model is constructed, and the preferred implementation process includes:

[0017] 11) To construct the discrete forward model for the required solution, the high-energy flash X-ray radiography system and imaging theory are first explained. The high-energy flash X-ray radiography system is as follows: Figure 2 As shown, it mainly consists of a light source, a rotationally symmetric target, a collimator, a scintillator, and an imaging detector.

[0018] High-energy flash X-ray sources are regional surface sources with continuous spectra, thus exhibiting spectral effects and source blurring. Source blurring and blurring caused by the imaging detector are collectively referred to as system blurring. During the interaction between high-energy X-rays and the target, scattered X-rays still pass through the scintillator, adversely affecting image quality. Furthermore, noise caused by the low energy conversion efficiency of the imaging detector also affects the quality of the projected image to some extent. Considering all these factors, the high-energy flash X-ray imaging equation can be expressed as follows: [Equation omitted for brevity]

[0019] g(x,y)=g(x0,y0)+f{(i0(x,y)T(x,y))*P s (x,y)+X s (x,y)}*P d (x,y)+n(x,y)(1)

[0020] In the formula, g(x,y) and g(x0,y0) are the static image and image background of the forward projection, respectively. f is the GX curve of the CCD system, representing the response function of the imaging detector. i0(x,y) represents the intensity of the rays received on the imaging plane in the empty field condition, indicating that there is no target in the optical path. P s (x,y) and P d (x, y) represent the point spread functions for light source blur and detector blur, respectively. s (x,y) represents the scattering amount on the imaging plane, n(x,y) represents the statistical noise term, and * indicates convolution operation. The product of i0(x,y) and transmittance T(x,y) is the direct irradiation amount.

[0021] 12) Subtract the image background from the static image and then divide it point-to-point by the empty field image to obtain the transmittance image. According to Beer's Law, the transmittance T(x,y) of high-energy flash X-rays penetrating the target is expressed as:

[0022]

[0023] Among them, S x (E) represents the source photon spectrum, μ ρ Let ρ be the mass absorption coefficient of the target and ρ be the target density. The product of the mass absorption coefficient and the density is equal to the linear absorption coefficient of the target.

[0024] Collimators can narrow the dynamic range of X-ray irradiation and effectively subtract the scattering intensity. Therefore, under monoenergetic irradiation, no scattering, and no noise conditions, the projected image received by the imaging detector is represented as:

[0025]

[0026] Here, B represents the system blur kernel composed of detector blur and light source blur. The transmittance image y is obtained by subtracting the image background from the forward-projected static image and then dividing it point-to-point by the empty field image.

[0027]

[0028] As can be seen, high-energy flash X-ray radiography is a typical nonlinear process. Converting the transmittance image y into vector form y, we obtain the discretized transmittance equation considering ambiguity and noise:

[0029]

[0030] Where, y∈R m Let G be the vector form of the transmittance image, where G∈R m×n Let x ∈ R be the forward projection matrix. n Let n be the linear absorption coefficient of the target to be reconstructed, n∈R m For noise, R es (·) represents the operation of converting a one-dimensional vector into an image matrix. This represents the transformation from an image matrix to a one-dimensional vector.

[0031] Optical path image R es The convolution of (Gx) with the fuzzy kernel is converted into matrix multiplication, yielding a discrete-form nonlinear forward model:

[0032] y = B·exp(-Gx) + n (6)

[0033] Where, y∈R m Transmittance data in vector form, B∈R m×m This is the fuzzy kernel extension form of matrix convolution after matrix multiplication, where G∈R m×n It is a positive matrix, x∈R n Let n be the linear absorption coefficient of the target to be reconstructed, n∈R m This is the noise term.

[0034] The form of the blur matrix B in a real photographic system is usually unknown. In this invention, the blur kernel of the transmittance image is estimated using existing image restoration methods.

[0035] 2) Derive the Jacobian matrix form of the residual matrix of the nonlinear positive model in step (1), and optimize the solution by constructing a least squares model to obtain the deterministic reconstruction result. The preferred implementation process includes:

[0036] The nonlinear positive model in step (1) is optimized by constructing a least squares model, as shown in equation (7), where r(x) represents the residual vector of the objective function. Non-negative constraints are applied to the optimized solution of the model to reduce the reconstruction error of the target vacuum region and low-density region.

[0037]

[0038] To solve this nonlinear model, this invention employs the non-negative constrained Levenberg-Marquardt (LM) algorithm to iteratively solve the objective function. The LM algorithm is based on the trust region algorithm and adds a damping term to the Gauss-Newton algorithm for regularization constraints. The descent direction Δx of each iteration is calculated using equation (8), and the damping coefficient is adjusted by using the gain ratio ρ to reflect the similarity between the objective function and its first-order Taylor expansion approximation.

[0039] (H+μI)△x=-g (8)

[0040]

[0041] In equation (8), H = J T J, g = J T r, J is the Jacobian matrix of the residual function.

[0042] Since the LM algorithm involves the calculation of the residuals of the objective function and the Jacobian matrix, it is necessary to derive the Jacobian matrix of the nonlinear model of high-energy flash X-ray radiography. The discrete form of the residuals of the objective function is expressed as follows:

[0043]

[0044] Where, x j This represents the value of the j-th element in the vector. B j,k Let $\frac{j}{k}$ represent the value in the $j-th row and $k-th column of matrix B. The corresponding Jacobian vector in the $j-th row is:

[0045]

[0046] Then the complete Jacobian matrix of the nonlinear positive model can be derived:

[0047]

[0048] In the formula, i and j represent the row and column of the matrix, respectively. [·] denotes the matrix form.

[0049] 3) Combining Bayesian theory to consider the solution and uncertainty quantification of this nonlinear inverse problem, a nonlinear hierarchical Bayesian model is constructed by introducing hyperparameters based on weak prior information. The preferred implementation steps include:

[0050] 31) First, under the definition of the nonlinear forward model, it is assumed that the noise n has a mean of zero and a covariance of λ. -1 Let I be a Gaussian random variable, and λ be the noise precision parameter. Then define the Gaussian likelihood function p(y|x,λ):

[0051]

[0052] Assume that the prior probability density function of the target parameter x follows a Gaussian distribution, and that it is a vector with zero mean and covariance δ. -1 Γ pr The precision matrix L is a discrete differential matrix with Neumann boundary conditions.

[0053]

[0054] The Jeffreys prior, based on weak information, is unaffected by changes in parameter form and yields more accurate parameter estimation results compared to the conditional conjugate Gamma prior. Therefore, this invention defines hyperparameters δ and λ based on the Jeffreys prior. λ is defined as 1 / u. 2 ,δ=1 / σ 2 And ζ=σ 2 / u 2 Jeffreys' priors satisfy:

[0055]

[0056] p(σ 2 |u 2 )=u -2 (1+σ 2 / u 2 ) -2 (16)

[0057] Combining the above equation, we obtain the form of the joint posterior probability density function:

[0058]

[0059] The fully conditional probability form of this joint posterior distribution is obtained as follows:

[0060]

[0061]

[0062]

[0063] 32) For the hyperparameter u in equation (18) above 2The sampling, whose distribution satisfies the inverse Gamma distribution, can be sampled using equation (21). Equation (19) can be transformed into the form of equation (22), where the first term ζ -(n / 2+1)-1 The exponential term of the third term also follows an inverse Gamma distribution, from the distribution... The sample value ζ generated by sampling is processed by the independent MH algorithm. 2 / (1+ζ) 2 The dominant acceptance rate calculation is used to determine whether to accept the current sample value.

[0064]

[0065]

[0066] Considering the nonlinear characteristics of f(x), it is not possible to derive the mean and variance of the target parameters in equation (20) as in a linear model and then directly sample them. The solution method will be explained in detail in step (4).

[0067] 4) By accelerating the solution of the optimization problem with random perturbation, the conditional distribution is sampled, and the solution of the optimization problem is constrained by the Jacobian matrix projection. A proposed distribution of the objective parameters is designed to reduce the statistical bias of the samples. The random perturbation strategy can ensure that samples are randomly generated in the high-probability region of the posterior distribution, thereby enabling the calculation of the confidence interval of the samples. The preferred implementation steps include:

[0068] 41) First, the conditional posterior probability density function shown in equation (20) is transformed into the following form so that it has the least squares form.

[0069]

[0070]

[0071] This invention constructs the following stochastic optimization problem to accelerate the generation of sample values ​​of the target parameter from p(x|λ,δ,y).

[0072]

[0073] To solve this stochastic optimization problem, the LM algorithm with non-negativity constraints from step (2) is used iteratively. The Jacobian matrix is ​​then:

[0074] J M =[λ 1 / 2 J,δ 1 / 2 L] T (26)

[0075] To prevent the iterative solution of equation (25) from failing to converge effectively and thus affecting the reconstruction accuracy, this invention proposes to use the Jacobian matrix J of the next sampling state (t+1)M In the residual matrix r M The method involves comparing the projection onto the given state with the projection values ​​from previous empirical states. If the adaptive constraints shown in the following equation are not satisfied, the solution is considered invalid.

[0076]

[0077] 42) Then, the MH algorithm is introduced to calculate the acceptance rate by designing the probability density function of the sample, thereby reducing the sample statistical bias obtained by equation (25) and realizing the sample debiasing estimation. The random perturbation equation of the observed data and target parameters can be expressed as:

[0078]

[0079] In the formula, ε y ~N(0,I m ),ε x ~N(0,I n ).

[0080] Combining equations (25) and (28), the MAP estimate can also be expressed as:

[0081]

[0082] In the formula,

[0083] To calculate the probability density function of the sample, we take the first-order partial derivative of the above equation, and obtain:

[0084]

[0085] Variables can be constructed based on the above formula. and Mapping function between And assume that the function is locally invertible.

[0086]

[0087] because The distribution form is known, and its probability density function is defined as follows: probability density function Then the determinant of the Jacobian matrix of the mapping function and The product of is shown in equation (33).

[0088]

[0089]

[0090] Among them, |J x | represents the mapping function The determinant of the Jacobian matrix, which contains f(x). T The second derivative of the equation is quite complicated in both theoretical derivation and computation, so only the first derivative term is retained, as shown in the following equation.

[0091]

[0092] 43) The above formula |J x The calculation of the determinant in high-dimensional models involves inverting the covariance matrix and multiplying first derivatives. Since the forward projection process sequentially scans each row of the linear absorption coefficient matrix to generate projected data, adjacent rows are independent in linear models, but in nonlinear models, any row of data is affected by the fuzzy kernel, influenced by the data in the rows above and below it within an interval defined by the radius of the fuzzy kernel. It is a symmetric matrix with values ​​concentrated in the region near the diagonal. (Linear and nonlinear models) The non-zero regions are distributed in blocks along the diagonal, called "effective matrix blocks". Assuming the linear absorption coefficient matrix has M rows and M columns, and the fuzzy kernel size is k×k, in the linear model, any effective matrix block is (M / 2)×(M / 2) and corresponds to any row of linear absorption coefficient values. In the nonlinear model, any row of linear absorption coefficients corresponds to an effective matrix of size ((2k-1)M / 2)×((2k-1)M / 2). Therefore, this invention uses the effective matrix corresponding to any row of linear absorption coefficients to approximate the calculation of |J| instead of the complete first-order derivative product. x The number of rows and columns in the precision matrix is ​​reduced from (M / 2)×M to (M / 2)×(2k-1), effectively reducing the computational cost of acceptance rate calculations for high-dimensional data. This invention selects the effective matrix corresponding to the M / 2th row of the line absorption coefficient matrix, i.e. The central region is used to calculate J. x |

[0093] probability density function The calculation is quite complex, so this invention approximates the sample probability density function as:

[0094]

[0095] Based on the calculated candidate samples and probability density function In the MH algorithm, the sample x of the next state * The probability of being accepted is α λ,δ The method of subtracting the logarithms first and then taking the exponent is used to calculate α. λ,δ The calculation is performed to prevent instability in numerical computation at high dimensions.

[0096]

[0097] 5) The sample values ​​of linear and nonlinear Bayesian models are fused under the variance criterion to obtain the final reconstructed image. The preferred implementation steps include:

[0098] The aforementioned random perturbation algorithm can effectively solve the sampling problem of nonlinear models. However, in Bayesian hierarchical models, due to the sampling of hyperparameters, a large number of samplings are required to ensure that inter-chain and intra-chain states reach a stationary state, which is difficult to satisfy when solving high-dimensional nonlinear models. This invention studies a fusion strategy that introduces multi-fidelity technology into the nonlinear reconstruction of high-energy flash X-ray images based on a simplified linear substitution model, in order to improve the efficiency and accuracy of sample estimation for high-dimensional nonlinear models.

[0099] 51) Without loss of generality, the nonlinear model of the present invention defines the input domain Y and the output domain X respectively, and establishes the input domain X according to equation (25). and output The mapping function between them is H:Y→X, where N s The total number of samples counted for the model. And define l models H. (1) H (2) ,…,H (l) H (1) The first model is the nonlinear model studied in this invention; the others are alternative models. It is worth noting that for these models, the same input domain Y is used, but the total amount of input samples N differs. s =[N s1 N s2 ,…,N sl ] T We obtain sample estimates under each model, thereby shifting the computational load from nonlinear models to alternative models and improving the performance of sample statistics.

[0100] Consider l models H (1) H (2) ,…,H (l) And set the corresponding total number of samples vector N. s =[N s1 N s2 ,…,N sl ] T Satisfying N s1 ≤N s2 ≤…≤N sl And assume that the first N have all been discarded. b A non-stationary sample. For nonlinear or alternative models N sl There are several input random variables. For i = 1, ..., l, model H needs to be modified. (i) Sample N using the MCMC algorithm si Next, we get N. siIndividual sample values:

[0101]

[0102] Sample mean E[H (i) (Y)] can be estimated unbiasedly. get:

[0103]

[0104] The posterior expectation of the nonlinear model can be calculated by equation (39). In order to obtain a more accurate sample estimate and less computational burden, and to explore the influence of alternative models on the posterior expectation of the nonlinear model, this invention constructs the Monte Carlo estimation fusion model shown in equation (40) to obtain the final reconstructed image.

[0105]

[0106]

[0107] In equation (40), the first term is the nonlinear model term, the second term is the substitution model term, and the coefficient ω i =[ω i,1 ,ω i,2 ,…,ω i,n ] T Used to balance the weights between different models.

[0108] 52) The fusion model defined by equation (41) depends on the coefficients ω1, ω2, ..., ω l and the total number of model evaluations N s Based on the Monte Carlo estimation fusion model that minimizes the mean squared error (MSE), an optimized model with minimum variance for the model coefficients and the total number of evaluations can be constructed.

[0109]

[0110] Where e(·) represents the model’s MSE and Var[·] represents the variance.

[0111] Define H (i) (Y) The variance Var[H] at any point q in x (i) [Y,q] and the Pearson product-moment correlation coefficient γ between different models i,j,q :

[0112]

[0113] γ i,j,q =Cov[H (i) (Y,q),H (j) (Y,q)] / (σ i,q σj,q (43)

[0114] In the formula, Cov[·] represents the covariance. For different models, the unbiased estimate at point q... and Its covariance is derived as follows:

[0115]

[0116] Based on the fact that the variance of multiple random variables equals the sum of their covariances, the variance of the fusion model can be derived:

[0117]

[0118] Therefore, according to equation (45), we can construct about O(ω) 1,q ,…,ω l,q N s The optimization problem is to minimize the total number of fusion models by assigning appropriate coefficients and model evaluations. The variance.

[0119]

[0120] For the alternative model in this nonlinear problem, this invention selects the simplified linear model as the only alternative model, explores its impact on the reconstruction results in model fusion, and sets N... s Assuming N is a constant that is pre-set, s1 =γ N ×N s2 , 0≤γ N ≤1 represents the scaling factor. In this case, the variance of the fusion model can be expressed as:

[0121]

[0122] For the above equation ω 2,q Taking the derivative, we get O(ω) 1,q ,ω 2,q The pixel-wise ω corresponding to the minimum value 2,q Value, and ω 2,q Substituting into equation (47), we get O(ω) 1,q ,ω 2,q Minimum value O(ω) 1,q ,ω 2,q ) min .

[0123]

[0124]

[0125] For any pixel q, when the variance of the fusion model as shown in equation (50) is less than the variance of the posterior expectation of the nonlinear model, a more accurate estimate of the sample mean can be obtained. At this time, the condition shown in equation (51) needs to be met.

[0126]

[0127]

[0128] In equation (51), N1 represents The total number of samples.

[0129] The above provides a detailed description of the nonlinear reconstruction method for high-energy flash X-ray images based on MCMC proposed in this invention, which can be widely applied in the field of internal density and interface inversion of high-energy flash X-ray radiographed targets.

[0130] The above embodiments are only used to illustrate the present invention. For those skilled in the art, several improvements and modifications can be made without departing from the basic principles of the present invention, and these improvements and modifications should be considered within the scope of protection of the present invention.

Claims

1. A nonlinear reconstruction method for high-energy flash X-ray images based on MCMC, characterized in that: 1) Based on the principle of high-energy flash X-ray imaging, a discrete nonlinear forward model is constructed; 2) Derive the Jacobian matrix form of the residual matrix of the nonlinear positive model in step (1), and optimize the solution by constructing a least squares model to obtain the deterministic reconstruction result; The discrete-form nonlinear forward model is shown below: y=B·exp(-Gx)+n (6) where, y∈R m Transmittance data in vector form, B∈R m×m Let G be the system fuzzy matrix, and G∈R m×n It is a positive matrix, x∈R n Let n be the target parameters to be reconstructed, where n ∈ R. m Noise term; The nonlinear forward model y=B·exp(-Gx)+n is optimized by constructing a least squares model, as shown in equation (7), where r(x) represents the residual vector of the objective function, x * The optimal solution for the linear absorption coefficient; Solving the nonlinear forward model involves calculating the residual vector of the objective function and the Jacobian matrix. The Jacobian matrix of the nonlinear forward model for high-energy flash X-ray radiography is derived, and the discrete form of the residual vector of the objective function is expressed as follows: Where, x i x j Let B represent the values ​​of the i-th and j-th elements in the vector, respectively. j,k Let $\matrix$ be the value in the $j$-th row and $k$-th column of matrix $B$. The corresponding Jacobian vector in the $j$-th row is: Then the complete Jacobian matrix J of the nonlinear positive model is derived: In the formula, i' and j' represent the rows and columns of the matrix, respectively, [·] represents the matrix form, B is the system fuzzy matrix, and G is the positive matrix; 3) Introduce hyperparameters based on weak information priors to construct a nonlinear hierarchical Bayesian model, and obtain the fully conditional probability form of the joint posterior distribution; 4) To accelerate the solution of optimization problems with random perturbations, the conditional distribution is sampled, and the solution of the optimization problem is constrained by the Jacobian matrix projection. A proposed distribution of the objective parameters is designed to reduce the statistical bias of the samples. 5) The sample values ​​of linear and nonlinear Bayesian models are fused under the minimum variance criterion to obtain the final reconstructed image.

2. The nonlinear reconstruction method for high-energy flash X-ray images based on MCMC according to claim 1, characterized in that: In step (3), the construction of the nonlinear hierarchical Bayesian model based on weak information prior includes: 21) Under the definition of the nonlinear forward model, assume that the noise n has a mean of zero and a covariance of λ. -1 I is a Gaussian random variable, and λ is the noise precision parameter; it is assumed that the prior probability density function of the target parameter x follows a Gaussian distribution, and satisfies the conditions of a zero-mean vector and a covariance of δ. -1 Γ pr The precision matrix L is a discrete differential matrix with von Neumann boundary conditions: Introduce hyperparameters δ and λ based on weak information priors, and define λ = 1 / u 2 ,δ=1 / σ 2 And ζ=σ 2 / u 2 We obtain the fully conditional probability form of the joint posterior distribution:

3. The nonlinear reconstruction method for high-energy flash X-ray images based on MCMC according to claim 2, characterized in that: Step (4) also includes: 31) Construct the following stochastic optimization problem, and use the LM algorithm to accelerate the generation of sample values ​​for the objective parameters. The Jacobian matrix is ​​J. M =[λ 1 / 2 J,δ 1 / 2 L] Τ : Through the Jacobian matrix J of the next sampled state (t+1) M In the residual matrix r M The solution to equation (25) is constrained by comparing the projection onto the previous empirical state with the projection value. If the adaptive constraint shown in the following equation is not satisfied, it is considered an invalid solution. 32) The MH algorithm is used to calculate the acceptance rate by designing the probability density function of the samples, thus removing statistical bias. The random perturbation equations for the observed data and target parameters are expressed as follows: In the formula, ε y ~N(0,I m ),ε x ~N(0,I n ); Combining equations (25) and (28), the MAP estimate is expressed as: In the formula, Calculate the probability density function of the sample, take the first-order partial derivative of the above equation, and construct the variables. and Mapping function between And assume that the function is locally invertible: x, probability density function The determinant of the Jacobian matrix of the mapping function and The product of is shown in equation (33): Among them, |J x | represents the mapping function The determinant of the Jacobian matrix, retaining its first derivative term, is shown in the following equation: 33) Assume the number of rows and columns of the line absorption coefficient matrix is ​​M, the fuzzy kernel size is k×k, the size of any effective matrix block in the linear model is (M / 2)×(M / 2) and corresponds to any row of line absorption coefficient values, and in the nonlinear model, any row of line absorption coefficients corresponds to an effective matrix of size ((2k-1)M / 2)×((2k-1)M / 2); select the effective matrix corresponding to the M / 2th row of the line absorption coefficient matrix, i.e., ▽f(x). Τ The central region of ▽f(x) is used to approximate the calculation of |J x |; Sample probability density The function is approximated as: Based on the calculated candidate samples and probability density function In the MH algorithm, the sample x of the next state * The probability of being accepted is α λ,δ The method of subtracting the logarithms first and then taking the exponent is used to calculate α. λ,δ Calculation:

4. The nonlinear reconstruction method for high-energy flash X-ray images based on MCMC according to claim 3, characterized in that: Step (5) also includes: 41) Define the input domain Y and output domain X according to the nonlinear model, and establish the input domain X according to equation (25). and output The mapping function between them is H:Y→X, where N s The total number of samples counted for the model; define l models H (1) H (2) ,…,H (l) H (1) The model is a nonlinear model, and the rest are alternative models of this nonlinear model; the Monte Carlo estimation fusion model shown in equation (40) is constructed to obtain the final reconstructed image; In equation (40), the first term is the nonlinear model term, the second term is the substitution model term, and the coefficient ω i =[ω i,1 ,ω i,2 ,…,ω i,n ] Τ Used to balance the weights between different models; 42) Construct an optimization model for the model coefficients and the total number of evaluations with minimum variance: Where e(·) represents the model's MSE, and Var[·] represents the variance; The variance of the fusion model is obtained: The simplified linear model is selected as the only alternative model, and N is... s Assuming N is a constant that is pre-set, s1 =γ N ×N s2 , 0≤γ N ≤1 represents the scaling factor, and the variance of the fusion model is expressed as: Calculate O(ω) 1,q ,ω 2,q Minimum value O(ω) 1,q ,ω 2,q ) min : For any pixel q, when the variance of the fusion model, as shown in equation (50), is less than the variance of the posterior expectation of the nonlinear model, a more accurate estimate of the sample mean is obtained, which requires satisfying the condition shown in equation (51): In equation (51), N1 represents The total number of samples.

Images