CN111833412B - A Tikhonov regularized image reconstruction method based on fractional filtering framework - Google Patents

Abstract

The invention belongs to the field of medical image processing, and specifically relates to a Tikhonov regularized photoacoustic imaging reconstruction method based on a fractional filtering framework, which includes the following steps: S1. Use Matlab simulation to build a system matrix to obtain a forward model of photoacoustic imaging; S2. Find the functional of Ax=b through Tikhonov regularization, and choose the fractional order in the (0,1) interval; S3, use the Morozov deviation principle to find λ; S4, find the reconstructed solution and filter function; S5, if S4 is obtained If the reconstruction solution and filter function can maximize the SNR or CNR value of the reconstructed image, then the most appropriate regularization parameter value and fractional order are obtained. Otherwise, repeat S2-S5 to further narrow the search interval. The fractional-order Tikhonov regularization of the present invention has higher reconstructed image quality than the traditional standard Tikhonov regularization. The superior performance of the fractional filtering framework is attributed to the inclusion of fractional order, rather than purely integer order regularization, which controls smoothness by increasing the norm of the reconstructed solution. The invention is used for the reconstruction of photoacoustic imaging.

Description

A Tikhonov regularized image reconstruction method based on fractional filtering framework

Technical field

The invention belongs to the field of medical image processing, and specifically relates to a Tikhonov regularized photoacoustic imaging reconstruction method based on a fractional filtering framework.

Background technique

As an auxiliary medical method, medical imaging can help humans detect, determine, understand and study diseases. Existing medical imaging methods are not enough to help humans explore and study diseases. Photoacoustic imaging (PAT), as a new type of non-destructive biomedical imaging technology, combines the advantages of pure optical imaging and pure ultrasound imaging, and has the advantages of strong contrast, high sensitivity, and deep imaging depth.

Photoacoustic imaging uses modulated laser light to irradiate target tissue, causing the tissue to heat and excite ultrasonic waves. The collected sound waves are then used to reconstruct the light absorption distribution that reflects the internal structure of the target tissue. The reconstruction algorithm will affect the imaging quality of the image, so its selection is very important. Existing reconstruction algorithms include filtered back projection (FBP), time reversal, and Fourier transform-based reconstruction algorithms. These algorithms are based on the spherical Radon transform model, which has the inherent limitation of requiring a large number of data points around the target object to accurately estimate the initial pressure distribution. Due to the limited boundary data in the actual measurement process and the inevitable mixing of noise in the measured photoacoustic signals, the recovery of the initial pressure rise distribution based on the model is an ill-posed problem.

As a technique for solving ill-conditioned inverse problems, regularization is widely used. There are also many classic regularization algorithms, such as Tikhonov regularization, which can make the reconstructed image have a small gradient, and Laplacian because of its isotropic nature. Evolution characteristics, better reconstruction effect in smooth areas of the image. However, due to the discontinuity of the image signal, in some non-smooth areas such as edges, the L2 norm used will cause excessive penalty on the image edges, making it difficult to generate accurate reconstructed images.

Contents of the invention

In view of the above technical problem that it is difficult to generate an accurate reconstructed image due to the discontinuity of the traditional Tikhonov regularized image signal, the present invention provides a Tikhonov regularized image reconstruction method based on a fractional filtering framework with high performance, high accuracy and small calculation amount.

In order to solve the above technical problems, the technical solution adopted by the present invention is:

A Tikhonov regularized image reconstruction method based on a fractional filtering framework, which is characterized by: including the following steps:

S1. Use Matlab simulation to construct the system matrix A, and obtain the forward model of photoacoustic imaging: Ax=b, where x is the initial value of the pressure of each pixel in the imaging area, and b is the measured acoustic data on the boundary;

S2. Find the functional of Ax=b through Tikhonov regularization, and select fractional order α in the (0,1) interval;

S3. Use the Morozov deviation principle to find the parameter λ, so that the residual norm is equal to the a priori upper limit δ;

S4. Obtain the reconstructed solution x _λ,α =VF(U ^T b) and filter function F=diag(σ ^α /(σ ^α+1 +λ));

S5. The reconstructed solution and filter function obtained from S4 can further establish the relationship between fractional α and noise level. The fractional α decreases as the noise increases, by changing the SNR or CNR value of the collected PA data. Change the fractional order α to maximize the SNR or CNR value of the reconstructed image, and obtain the most appropriate regularization parameter value λ _opt and fractional order α _opt . Otherwise, repeat S2-S5 to further narrow the search interval.

The method of constructing the system matrix A in S1 is to obtain the basic equation of photoacoustic imaging through the thermal equation, motion equation and diffusion equation of photoacoustic imaging: Said/> is the three-dimensional space position vector, the t is the time, the p is the sound pressure, the c is the sound speed, the β is the isobaric expansion coefficient, the C _p is the specific heat, let I be the delta function, and define The pulse emission time is the time zero point, that is, I(t)=δ(t). The solution of the basic equation of photoacoustic imaging is approximated by the Green function to obtain the photoacoustic signal:/> The photoacoustic signal is simulated through k-Wave, and the process of collecting power amplifier data at the sensor location is represented as a time-varying causal system. The toolbox of k-Wave is used to construct the system matrix A of the time-varying causal system. The system matrix A includes the impulse response of each pixel in the imaging area as a column. The imaging grid is n×n pixels. The imaging grid stacks all columns together and converts them into a high column vector of size n ² × 1. The initial pressure rise of the pixels of the imaging grid is x, and the system matrix A has a dimension of m×n ² , resulting in a forward model of photoacoustic imaging: Ax=b.

The functional function of Ax=b in S2 is: The weighted norm The W is a symmetric positive semi-definite matrix, and the /> The α is the fractional order α>0; when α=1, the matrix W=I, the I is the identity matrix, the weighted norm allows the selection of the fractional order α, and the fractional order α is in (0,1) Taking an appropriate value for the interval can reduce the shortcoming of the standard regular solution being too smooth, and the appropriate value is a value that maximizes the SNR or CNR value of the reconstructed image.

The method of calculating the parameter λ using the Morozov deviation principle in S3 is: let e be the observation error, and there is an upper bound ε, satisfying ||e||≤ε,

Define δ=ηε, so the solution x _λ,α satisfies ||b-Ax _λ,α ||=δ,

Let G(λ)=||b-Ax _λ,α || ² -δ ² =0,

When α is given, the G(λ)=||b-Ax _λ,α || ² -δ ² =0 is a nonlinear implicit equation about λ, and the Newton method is used to iteratively solve λ, that is

described is the iteration step size, said/> The expression of G'(λ) is:

described By minimizing x _Tikh = (A ^T A+λI) ^-1 A T b with respect to x by minimizing x _Tikh with respect to x, we can get the derivative of A ^T b with respect to λ, that is

Combined with the residual norm of the standard Tikhonov algorithm:

The expression that can be derived for the step size is

in

The residual residual sequence of the least squares solution is constructed from the residual norm of the standard Tikhonov algorithm, that is

When δ is given, based on min|ρ(i)-δ ² |, we can know the sequence number i when the expression takes the minimum value, and let the iteration initial value λ ₀ =σ _i .

The method for obtaining the reconstructed solution x _λ,α =VF( ^UT b) and filter function F=diag(σ ^α /(σ ^α+1 +λ)) in S4 is:

According to the deviation principle ||b-Ax _λ || ₂ = δ, we can get that λ is a function of δ, where F is used to replace A as an n×m matrix, and its diagonal elements are the filter factors of the fractional Tikhonov format, that is

get

Considering k=min(m,n ² ), then Then take the differential, that is, take δ(λ) as the inverse function of λ(δ) to get, Then the reconstructed solution x _λ,α =VF( ^UT b) and filter function F=diag(σ ^α /(σ ^α+1 +λ)) are obtained.

Compared with the prior art, the beneficial effects of the present invention are:

The fractional-order Tikhonov regularization of the present invention has higher reconstructed image quality than the traditional standard Tikhonov regularization. The superior performance of the fractional filtering framework is attributed to the inclusion of fractional order α, which controls smoothness by increasing the norm of the reconstructed solution, rather than purely integer order regularization. The regularization parameters are automatically selected by the deviation principle. Using this method for image reconstruction not only requires a small amount of calculation, but also provides a better balance between noise suppression and feature preservation, resulting in a more reasonable reconstructed image.

Description of the drawings

Figure 1 shows the process of automatically selecting fractional order and regularization parameters in this invention.

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 some of the embodiments of the present invention, rather than all 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 fall within the scope of protection of the present invention.

A Tikhonov regularized image reconstruction method based on the fractional filtering framework. When electromagnetic wave pulses are irradiated on the biological tissue to be measured, the tissue absorbs electromagnetic energy and generates heat. At this time, the thermal equation of photoacoustic imaging can be given:

in is the three-dimensional space position vector, t is time, function H represents the heat absorbed by the tissue, function T represents the temperature at which the temperature rises, ρ is the density of the tissue (assumed to be a constant here), C _p is the specific heat, and λ is the thermal conductivity . Usually, the duration of the electromagnetic wave pulse is much shorter than the heat conduction time, so the heat conduction can be ignored, so (1) can be simplified to:

The tissue expands due to heat, thus creating an initial sound field that changes over time (i.e., sound waves are generated). At this time, the motion equation and diffusion equation of the photoacoustic imaging sound field can be given:

where the vector function represents sound displacement, function p represents sound pressure, c is sound velocity, and β is the isobaric expansion coefficient. Simultaneous equations (2)-(4), writing the function H as the product of the light absorption distribution function A with respect to space and the irradiation pulse energy function I with respect to time can be obtained:

Equation (5) can be used as the basic equation of photoacoustic imaging, expressing the relationship between the photoacoustic signal and the light absorption characteristics of the tissue, and is a forward model describing photoacoustic imaging. Assume that I is a delta function, and define the pulse emission time as the time zero point, that is, I(t)=δ(t); the solution of equation (5) can be approximated by the Green function:

The linear equation obtained by discretizing equations (2) to (4) is Ax=b. Due to the limited boundary data in the actual measurement process and the inevitable mixing of noise in the measured photoacoustic signal, the linear equation cannot be solved directly. . In order to reconstruct the initial sound pressure distribution, as shown in Figure 1, the Tikhonov regularization method is commonly used, and its functional expression is

In the formula, A is the system matrix, including the impulse response of each pixel in the imaging area as a column, x is the initial value of the pressure of each pixel in the imaging area, and b is the measured acoustic data on the boundary (detector position). λ is a regularization parameter used to balance the remainder of the linear equation (first term on the right) and the desired initial pressure distribution (x). Higher regularization tends to over-smooth the image, while smaller λ values amplify noise in the image. It can be derived that the functional Γ _Tikh is minimized with respect to x as

x _Tikh =(A ^T A+λI) ^-1 A ^T b (8)

In the formula, I is the identity matrix. When λ=0, equation (8) is the traditional least squares solution. Substituting the singular value decomposition A=U∑V ^T of A into equation (8), the diagonal elements of ∑ are singular value sequences, then the standard Tikhonov solution x _Tikh can be expressed as

In the formula, the filter factor is

It can be seen from equation (9) that the solution norm and residual norm of the standard Tikhonov algorithm are respectively

Usually x _Tikh is too smooth, that is, the filter factor denoising process loses the data details contained in the exact solution. In response to this shortcoming, the present invention proposes a Tikhonov regularization method based on a fractional filtering framework. This method uses the fractional order of the Moore-Penrose pseudo-inverse of the AA ^T matrix as the semi-norm of the weighting matrix to measure the residual error in Tikhonov regularization. plan. The generic function in equation (7) is rewritten as

In the formula, the weighted norm Here W is a symmetric positive semidefinite matrix, and the matrix W is as follows

In the formula, α represents the fractional order of α>0. When α=1, that is, the matrix W=I, equation (13) is the standard Tikhonov regularization. Equation (13) has a unique solution for all positive values of the regularization parameter λ. The semi-norm _|| .|| quality. Differentiate equation (13) with x and equal it to zero, the result is

((A ^T A) ^α+1/2 +λI)x＝(A ^T A) ^α-1/2 A ^T b (15)

Equation (15) can be rewritten as,

(A ^T WA+λI)x=A ^T Wb (16)

SVD of A in Equation (16), we get,

(VSU ^T US ^α-1 U ^T USV ^T +λI)x＝VSU ^T US ^α-1 U ^T b (17)

(VS ^α+1 V ^T +λI)x＝VS ^α U ^T b (18)

The regular solution expression of fractional-order Tikhonov is

Different regularization methods contain different φ(σ _i ). When σ _i gradually tends to 0, φ(σ _i ) also tends to 0, and the too fast convergence of φ(σ _i ) is the reason why the solution vector is too smooth, that is, the small singular value associated vectors are lost, and these vectors reflect are the details of reconstructing the data. The filter factor φ _frac,α (σ) of fractional order Tikhonov is expressed as

In the present invention, the fractional order (α) is automatically selected based on maximizing the SNR/CNR of the reconstructed image. Taking an appropriate value for α slows down the convergence speed of φ _frac,α (σ). Therefore, compared with the standard Tikhonov method, the present invention is based on Tikhonov regularization of the fractional filtering framework results in higher accuracy solutions.

In the traditional Tikhonov regularization method, no matter how noisy the environment is, larger and smaller singular values are assigned the same weight and processed similarly. However, in noisier environments, to obtain good image reconstruction, smaller singular values require smaller weights, that is, in the weighted kernel norm, smaller singular values should be assigned smaller weights. Therefore, compared with the standard method based on Tikhonov, this project introduces fractional-order α into the weighting matrix. The fractional-order α can be changed according to the size of the noise, thereby controlling the damping and smoothness in the reconstruction solution to achieve a better Accurate reconstruction.

Singular vectors with smaller singular values usually contain high-frequency oscillations. Analyzing the filtering factors of the fractional Tikhonov method (given in Equation (20)) for small singular values, i.e. σ<<1, using Taylor series to expand Equation (20), ignoring higher order terms, we get,

Treat φ _t '(α) = σ ^α in equation (21) as a function of α and differentiate with respect to α,

φ _t '(α)＝σ ^α ln(σ) (22)

For this equation, consider σ<<1, where _φt (α) is a decreasing function, because when σ<<1, _φt '(α)<0. Thereafter, each term in equation (21) is a decreasing function of α. Therefore, reducing the fractional order from α = 1 (Tikhonov) will increase φ (α), which means that the high-frequency components in the fractional-order Tikhonov reconstructed solution increase, that is, the smoothness of the image is reduced.

For a large-scale inversion problem such as PAT image reconstruction, it is relatively difficult to choose an appropriate regularization parameter λ so that the algorithm has a better noise suppression effect while ensuring that it does not cause over-regularization problems. If the regularization parameter λ is chosen too small, the solution vector norm ||x _λ || of the equation is too large, and the solution is unstable; if λ is chosen too large, the denoising process is too strong, and _|| , the solution is too smooth. When the observation disturbance δ is known, the present invention adopts an automatic selection method of regularization parameters based on the Morozov deviation principle. Let e be the observation error, and there is an upper bound ε, satisfying

||e||≤ε (23)

Define δ = ηε, so the solution x _λ,α satisfies

||b-Ax _λ,α ||=δ (24)

make

G(λ)＝||b-Ax _λ,α || ² -δ ² =0 (25)

When α is given, equation (25) is a nonlinear implicit equation about λ, and the Newton method is used to iteratively solve λ, that is

In the formula, is the iteration step size,/> After calculation, the expression of G'(λ) is

In the formula, By deriving the derivative of λ from equation (8), we can get:

Combining Equation (12), the expression of the step length can be derived as

in

Whether the Newton method iteration sequence converges depends on the selection of the initial value λ _0. To address this feature, the following processing is adopted. The residual sequence of the least squares solution is constructed from Equation (12), that is

When δ is given, based on min|ρ(i)-δ ² |, we can know the sequence number i when the expression takes the minimum value, and let the iteration initial value λ ₀ =σ _i .

Only the preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the above-mentioned embodiments. Within the scope of knowledge possessed by those of ordinary skill in the art, various modifications can be made without departing from the spirit of the present invention. All changes should be included in the protection scope of the present invention.

Claims (1)

1. A Tikhonov regularized image reconstruction method based on fractional filtering framework, which is characterized by: including the following steps: S1. Use Matlab simulation to construct the system matrix A, and obtain the forward model of photoacoustic imaging: Ax=b, where x is the initial value of the pressure of each pixel in the imaging area, and b is the measured acoustic data on the boundary; S2. Find the functional of Ax=b through Tikhonov regularization, and select fractional order α in the (0,1) interval; S3. Use the Morozov deviation principle to find the parameter λ, so that the residual norm is equal to the a priori upper limit δ; S4. Obtain the reconstructed solution x _λ,α =VF(U ^T b) and filter function F=diag(σ ^α /(σ ^α+1 +λ)); S5. The reconstructed solution and filter function obtained from S4 can further establish the relationship between fractional α and noise level. The fractional α decreases as the noise increases, by changing the SNR or CNR value of the collected PA data. Change the fractional order α to maximize the SNR or CNR value of the reconstructed image, and obtain the most appropriate regularization parameter value λ _opt and fractional order α _opt . Otherwise, repeat S2-S5 to further narrow the search interval; The method of constructing the system matrix A in S1 is to obtain the basic equation of photoacoustic imaging through the thermal equation, motion equation and diffusion equation of photoacoustic imaging: Said/> is the three-dimensional space position vector, the t is the time, the p is the sound pressure, the c is the sound speed, the β is the isobaric expansion coefficient, the C _p is the specific heat, let I be the delta function, and define The pulse emission time is the time zero point, that is, I(t)=δ(t). The solution of the basic equation of photoacoustic imaging is approximated by the Green function to obtain the photoacoustic signal:/> The photoacoustic signal is simulated through k-Wave, and the process of collecting power amplifier data at the sensor location is represented as a time-varying causal system. The toolbox of k-Wave is used to construct the system matrix A of the time-varying causal system. The system matrix A includes the impulse response of each pixel in the imaging area as a column. The imaging grid is n×n pixels. The imaging grid stacks all columns together and converts them into a high column vector of size n ² × 1. The initial pressure rise of the pixels of the imaging grid is x, and the system matrix A has the dimension of m×n ² , resulting in the forward model of photoacoustic imaging: Ax=b; The functional function of Ax=b in S2 is: Weighted norm/> The W is a symmetric positive semi-definite matrix, and the /> The α is the fractional order α>0; when α=1, the matrix W=I, the I is the identity matrix, the weighted norm allows the selection of the fractional order α, and the fractional order α is in (0,1) Taking an appropriate value in the interval can reduce the shortcoming of the standard canonical solution being too smooth, and the appropriate value is a value that maximizes the SNR or CNR value of the reconstructed image; The method of calculating the parameter λ using the Morozov deviation principle in S3 is: let e be the observation error, and there is an upper bound ε, satisfying ||e||≤ε, Define δ=ηε, so the solution x _λ,α satisfies ||b-Ax _λ,α ||=δ, Let G(λ)=||b-Ax _λ,α || ² -δ ² =0, When α is given, the G(λ)=||b-Ax _λ,α || ² -δ ² =0 is a nonlinear implicit equation about λ, and the Newton method is used to iteratively solve λ, that is described is the iteration step size, said/> The expression of G'(λ) is: described By minimizing x _Tikh = (A ^T A+λI) ^-1 A T b with respect to x by minimizing x _Tikh with respect to x, we can get the derivative of A ^T b with respect to λ, that is Combined with the residual norm of the standard Tikhonov algorithm: The expression that can be derived for the step size is in The residual residual sequence of the least squares solution is constructed from the residual norm of the standard Tikhonov algorithm, that is When δ is given, based on min|ρ(i)-δ ² |, we can know the sequence number i when the expression takes the minimum value, and let the iteration initial value λ ₀ =σ _i ; The method for obtaining the reconstructed solution x _λ,α =VF( ^UT b) and filter function F=diag(σ ^α /(σ ^α+1 +λ)) in S4 is: based on the deviation principle ||b-Ax _λ || ₂ = δ can be obtained as λ is a function of δ, where F is used to replace A as an n×m matrix, and its diagonal elements are the filter factors of the fractional Tikhonov format, that is get Considering k=min(m,n ² ), then Then take the differential, that is, take δ(λ) as the inverse function of λ(δ) to get, Then the reconstructed solution x _λ,α =VF( ^UT b) and filter function F=diag(σ ^α /(σ ^α+1 +λ)) are obtained.