A photoacoustic tomography image reconstruction method based on TV-CG
By combining total variational regularization and conjugate gradient method, the TV-CG method solves the problems of slow speed and blurred edges in photoacoustic tomography image reconstruction, achieving efficient and accurate image reconstruction while suppressing noise interference.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- SHANDONG NON METALLIC MATERIAL RESEARCH INSTITUTE
- Filing Date
- 2022-12-14
- Publication Date
- 2026-07-31
AI Technical Summary
Existing photoacoustic tomography image reconstruction methods suffer from problems such as slow reconstruction speed and blurred edges. Furthermore, due to limited boundary data and noise interference, it is difficult to generate accurate reconstructed images.
A TV-CG-based approach is adopted, combining total variation regularization and conjugate gradient method. By adding total variation perturbation during the iteration process, the iteration direction is optimized, and the combined total variation regularization and perturbation terms improve the reconstructed image effect.
It improves the speed and quality of image reconstruction, suppresses noise, maintains clear image edges, and provides higher performance reconstruction results.
Smart Images
Figure CN115953492B_ABST
Abstract
Description
Technical Field
[0001] This invention belongs to the field of nondestructive testing technology, specifically relating to a photoacoustic tomography image reconstruction method based on TV-CG. Background Technology
[0002] Photoacoustic tomography (PAT), as a novel non-destructive imaging technology, combines the advantages of pure optical imaging and pure ultrasonic imaging, and has the advantages of strong contrast, high sensitivity, and deep imaging depth.
[0003] Photoacoustic imaging uses modulated laser light to irradiate target tissue, causing the tissue to heat up and emit ultrasonic waves. The acquired sound waves are then used to reconstruct the light absorption distribution reflecting the internal structure of the target tissue. The reconstruction method affects the image quality, so its selection is crucial. Existing reconstruction methods include filtered back projection (FBP), time reversal, and Fourier transform-based methods. These algorithms are based on the spherical Radon transform model, which has inherent limitations, 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 actual measurements and the unavoidable noise in the measured photoacoustic signals, the recovery of the initial pressure rise distribution based on the model is an ill-posed problem.
[0004] Regularization, as a technique for solving ill-posed problems, has a wide range of applications and many classic regularization algorithms, such as TV regularization, which can reconstruct image edges well. The conjugate gradient method (CG) can quickly solve ill-posed problems, and the combination of the two has good reconstruction results in smooth image regions. However, due to the discontinuity of image signals, in some non-smooth regions such as edges, the L2 norm used can lead to excessive penalty for image edges, making it difficult to generate accurate reconstructed images. Summary of the Invention
[0005] To address the problems of slow reconstruction speed and blurred edges in photoacoustic tomography image reconstruction, this invention provides a high-performance, high-accuracy, and fast photoacoustic tomography image reconstruction method based on TV-CG.
[0006] The technical solution adopted in this invention is as follows: A photoacoustic tomography image reconstruction method based on TV-CG includes the following steps: S1. A forward model for photoacoustic imaging is established based on the photoacoustic wave equation, and the system matrix is constructed using the central difference method. Thus, measured acoustic data were obtained. Compared with the initial sound pressure distribution to be reconstructed Relationship ,in, This represents the initial sound pressure column vector corresponding to each pixel within the imaging region. This is a column vector of measured acoustic data on the detection boundary; S2. Initialize the reconstruction process, setting the initial solution as... and order Simultaneously, the initial gradient vector is calculated using the conjugate gradient method. Initial search direction and corresponding auxiliary quantities ,in , Then, based on the initial step size The solution is updated for the first time. ,get And set the iteration counter. ; S3, in the iteration counter No more than the preset maximum number of iterations And the current error function value is greater than the preset threshold. Under the given conditions, a TV-CG iterative update is performed on the current solution. The update process includes: first, updating the current solution... A total variational perturbation term is introduced for perturbation processing to obtain intermediate variables. The intermediate variable The current solution after incorporating the total variational perturbation is used as the input for subsequent conjugate gradient updates; then, the intermediate variable... Using the CG-PR-STEP gradient update rule as input, perform conjugate gradient update to calculate a new solution. The corresponding search direction and auxiliary quantities; finally, let Proceed to the next iteration for judgment; S4. When the number of iterations exceeds the preset maximum number of iterations. or the current error function value Not greater than the preset threshold When the time is up, terminate the iteration and output the reconstruction result. .
[0007] The system matrix in S1 The construction method is as follows: First, the basic equations of photoacoustic imaging are established based on the thermal equation, motion equation, and related wave relationship of photoacoustic imaging; in the process of photoacoustic image reconstruction, the key to the model-based reconstruction method lies in solving the system matrix. This invention uses the central difference method to solve the system matrix. The relationship between the acquired photoacoustic signals and the reconstructed images is represented in matrix form. .
[0008] To accurately establish the wave equation for the generation and transmission of photoacoustic signals, k-Wave simulation was performed. The process of acquiring photoacoustic signals at the sensor location can be represented as a time-varying causal system. The k-Wave toolbox was used to construct the system matrix of this system, and the impulse response (IR) of the imaging grid was stored pixel-by-pixel in the system matrix as columns of the considered geometry. Let the discrete regions of the reconstructed image be divided into... A grid of pixels, where the initial pressure increase at each pixel constitutes the vector to be reconstructed. Let the number of transducer elements in the imaging system be... The length of the signal collected at each array element position is The photoacoustic signals collected at the positions of each array element constitute the measurement vector. Therefore, the system matrix The dimension is OK, The forward model of photoacoustic imaging can be represented as follows: Each column represents the impulse response at the corresponding pixel location. .
[0009] The system matrix constructed in S1 It has a serious pathological condition, leading to It is an ill-conditioned underdetermined equation that requires a regularization method to solve. Based on the idea of photoacoustic tomography image reconstruction method of TV-CG, the iteration process is guided to the direction of the feasible point of optimization relative to the given function by adding TV perturbation in the original conjugate gradient iteration process. The conjugate gradient algorithm with added perturbation is combined with the total variation problem and applied to the photoacoustic reconstruction model.
[0010] The least squares (LS) problem of image reconstruction can be formulated as follows: The TV-CG algorithm is based on the formula Add TV perturbation item to the basis TV penalty items In the TV-CG algorithm, the coefficients of the total variation regularization term... and the coefficients of the total variational disturbance term Introduced into the iteration, reconstruction is transformed into an optimization problem. The reconstructed image is then expressed mathematically as: , The TV perturbation value is added; reconstruction using a nonlinear conjugate gradient algorithm with perturbation will produce a result containing total variation regularization and TV perturbation. Therefore, the reconstructed image... The following optimization problem needs to be solved: The The coefficients of the total variation regularization term are... Let be the coefficients of the total variational disturbance term, and and These are mutually independent adjustable parameters. The regularization penalty term added during image reconstruction aims to eliminate noise and other interference. For data fidelity items, For total variational regularization, The three factors, which are total variational perturbations, work together to improve the image reconstruction effect.
[0011] Based on the regularization idea of TV-CG, a TV perturbation is added to the original conjugate gradient iteration process, and the perturbation-added conjugate gradient algorithm is combined with the total variation problem and applied to the reconstruction model.
[0012] In processing photoacoustic data, the TV-CG algorithm processes data separately for each detection point when calculating the measurement matrix. This allows for updating the reconstructed image after each detection point iteration, avoiding the use of large matrices for computation and improving the efficiency of image updates during iteration, thus shortening image reconstruction time. In step S3, an iterative process involving total variational perturbation and conjugate gradient updates is performed on the current solution. Unlike common methods that improve the original optimization problem solely by adding penalty terms, the TV-CG algorithm guides the iterative process towards a feasible direction for optimizing the given objective function while maintaining the original optimization criteria; when the error function value is not greater than a preset threshold... Stop iteration when the time is reached. and When the preset value conditions are met respectively, the total variation regularization term and the total variation perturbation term work together to improve the reconstructed image effect while suppressing noise.
[0013] In photoacoustic image reconstruction methods, total variational regularization (TV) can effectively preserve image edges, while the conjugate gradient method (CG) can quickly solve ill-posed problems. This invention employs a TV-CG method combining total variational regularization and the conjugate gradient method. By introducing a TV perturbation into the original conjugate gradient iteration process, the iteration process is guided towards the optimization feasible direction of the given objective function. The perturbated conjugate gradient algorithm is then combined with the total variational problem and applied to the photoacoustic reconstruction model.
[0014] Compared with the prior art, the beneficial effects of this invention are: The TV-CG algorithm in this invention improves the quality of reconstructed images compared to the traditional standard TV regularization algorithm. Unlike common methods that improve the original optimization problem by adding penalty terms, the TV-CG algorithm guides the iterative process towards a feasible direction for optimizing the given objective function while maintaining the original optimization criteria; this is achieved when the error function value is no greater than a preset threshold. Stop iteration when the time is reached. and When the preset value conditions are met respectively, the total variation regularization term and the total variation perturbation term work together to improve the reconstructed image effect while suppressing noise.
[0015] The most significant feature of the conjugate gradient method used in this invention is that it uses a linear combination of the previous search direction and the negative gradient of the current starting point to generate the conjugate direction. This results in a smaller computational and storage requirement during the iteration process, and a significantly faster image reconstruction speed. Attached Figure Description
[0016] Figure 1 This is a map showing the pixel locations crossed by the integral curve; Figure 2 This is a flowchart of the photoacoustic tomography image reconstruction method based on the TV-CG algorithm. Detailed Implementation
[0017] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0018] The purpose of this invention is to propose a photoacoustic tomography method based on TV-CG, which can effectively improve the reconstruction effect of original images with arc-shaped edges and provide reconstructed images with higher performance than the traditional standard TV regularization scheme. To achieve the above objective, the technical solution adopted in this invention is a photoacoustic tomography image reconstruction method based on TV-CG.
[0019] Establishment of the wave equation:
[0020] When tissue is irradiated by a short-pulse laser, the change in the medium's temperature ultimately leads to a change in pressure; at this point, the photoacoustic pressure... It can be represented as: (1) in, It is the isobaric expansion coefficient. It is the speed of sound in the medium. The change in density The density of the medium, and These represent location information and time information, respectively. The temperature represents the increase in tissue temperature. The first term on the right side of the equation represents the contribution of the density change caused by the original sound field to the photoacoustic pressure, and the second term represents the effect of temperature change on the photoacoustic pressure.
[0021] Apply equation (1) to Differentiating, we get: (2) Will Substituting into equation (2), we have: (3) In the formula The velocity of the particle vibration, For relative position variables Vector differential operators.
[0022] Apply equation (3) After differentiation, substitute into Then we have: (4) Equation (4) shows the relationship between temperature change and photoacoustic pressure in the sound field.
[0023] When the pulse width of a pulsed laser is very short, such that the deposition time of laser energy on the medium is less than the diffusion time, the effect of volume change caused by temperature increase can be ignored. Therefore, temperature change is related to the deposition of light energy. The relationship between them is: (5) in, Specific heat capacity coefficient Indicates the temperature at which the tissue rises. The unit of heat source function representation Heat per unit volume over a given time.
[0024] Therefore, substituting equation (5) into equation (4), the relationship between photoacoustic pressure and photoelectric energy deposition is obtained as follows: (6) in, It is the Laplace operator relative to the position variable r, and equation (6) is the photoacoustic wave equation.
[0025] Assuming the medium through which the sound wave travels is inviscid and the density of biological tissue is uniform, the pressure distribution caused by the radial propagation of the sound wave generated by the photoacoustic effect at the speed of sound in the medium... With heat source function The photoacoustic wave equation of the relationship can be expressed as: (7) in The sound waves generated by the photoacoustic effect propagate radially at the speed of sound in the medium at the position. time The resulting pressure distribution is a function of time and location. The heat source function represents the amount of heat per unit volume per unit time.
[0026] As can be seen from equation (7), in order to generate sound waves, the energy of the irradiating laser must be variable, which is why pulsed laser excitation is chosen. Heat source function Space and time can be separated, that is... , It is the spatial electromagnetic absorption function. In most cases, the pulse duration of pulsed lasers is very short, typically a few nanoseconds, therefore... It can be approximated as a Function, i.e. Then equation (7) can be simplified to: (8) in, For position r, the duration is The initial sound pressure generated by pulsed laser excitation, and the sound pressure signal is expressed integrally. , For Grünersen parameters.
[0027] ultrasonic transducer in position The sound pressure signal detected at the location is also known as the photoacoustic signal. As a solution to formula (8), it is expressed as: (9) Equation (9) shows that the photoacoustic signal at position r and time t is the time derivative of the initial sound pressure distribution along the surface integral of a sphere with radius ct.
[0028] System Matrix Solving for: During the reconstruction process, the three-dimensional spherical problem needs to be simplified to a two-dimensional fault plane, and the integration object changes from the original sphere to a circular arc. Ignoring the constants in equation (9), equation (9) can be expressed as follows in the two-dimensional fault plane: (10) Given transducer position The polar coordinates of the point on the integral curve can then be expressed as: (11) In the formula Let be the polar angle of a point on the integral curve.
[0029] Equation (10) can be further simplified to: (12) The most crucial part of model-based reconstruction methods is solving for the system matrix. The central difference method will be used to solve for the system matrix below.
[0030] Using the central difference method, after simplifying equation (12), we get (13) In equation (13), Represented as: (14) In the formula This refers to the angular range covered by the transducer over the image reconstruction region.
[0031] Discretize the rectangular region containing the imaging object as follows: A grid, where the position of each pixel in the grid is represented as ,in, , The corresponding energy deposition is represented as Then equation (14) can be discretized as: (15) in, For sensor position, For a point in time, for The pixel position crossed by the integral curve at time step, pixel The distance between the sensor and the sensor meets the requirements. . Represents pixels The included angle corresponding to the intersecting arcs.
[0032] Combining equations (13) and (15), at time... ,Location Pressure measured at the point Represented as: (16) in, The elements in the system matrix represent the response relationship between the detection point and the pixel. pixel position Discrete values of the spatial absorption function at a given location.
[0033] The calculation of the system matrix involves the following steps: (1) Determine the angle The range, i.e., the calculation of the minimum angle and the maximum angle The and It relates to the first and last intersection points between the integral curve and the grid; (e.g.) Figure 1 As shown), they determined Scope; (2) Calculate the angles corresponding to each intersection point using polar coordinates. And make the obtained angle located at Within the range; (3) Obtaining various angles Then, calculate The system matrix is solved by combining equations (13), (15), and (16). .
[0034] Let the discrete regions of the reconstructed image be divided into The imaging system has a pixel grid and a transducer array element count of [number]. The length of the signal collected at each array element position is The photoacoustic signals collected at the positions of each array element and the pixel grayscale values of the reconstructed image Represented as a column vector: (17) Therefore, the relationship between the acquired photoacoustic signal and the reconstructed image can be expressed in matrix form: (18) Where, vector Let be the column vector of the image to be reconstructed. For length is column vectors, vectors For length is Column vectors. This is the forward model system matrix.
[0035] For a given photoacoustic tomography system, the system matrix depends only on the discrete image grid used and the sensor location, and is independent of the actual imaging object.
[0036] Due to the limited boundary data and the unavoidable noise in the measured photoacoustic signals during actual measurements, the linear equations cannot be solved directly. To reconstruct the initial sound pressure distribution, the TV regularization method can be used, and its functional expression is: (19) In the formula The system matrix contains the impulse responses of each pixel within the imaging region as columns. The initial pressure values for each pixel within the imaging region. This refers to the measured acoustic data at the boundary (detector location). It is a regularization parameter used to balance the remainder of the linear equation (the first term on the right) with the desired initial pressure distribution. Higher regularization tends to make images overly smooth, while lower regularization... The value will amplify noise in the image.
[0037] The total variational norm has two definitions in one-dimensional and two-dimensional spaces.
[0038] In one-dimensional space, suppose a function u is defined on the interval [0, 1]. On, then the function The total variation is defined as follows: (20) In two-dimensional space, that is For smooth functions The total variation in one-dimensional space can be defined as follows: (twenty one) In the formula, total variation , For time, It is a two-dimensional spatial position vector. If If a function is not smooth, then the total variation in two-dimensional space is defined as: (twenty two) in, , , Continuously differentiable, and (where Euclidean represents Euclidean space).
[0039] If u is a smooth function, then we have because It is not differentiable at the origin. To overcome this difficulty, consider... Approximate form Where parameters Therefore, TV(x) can be approximated as: (twenty three) The TV method does not impose any special restrictions on the smoothness of the solution, only requiring that the solution to the problem is bounded variation. This is also the advantage of the total variation method over other regularization methods in solving practical problems. However, the TV regularization method is relatively slow.
[0040] This invention proposes a photoacoustic tomography image reconstruction method based on TV-CG. This method utilizes the conjugate gradient method with added total variational perturbation, combined with TV regularization, for image reconstruction. Based on the least squares objective function, by introducing a total variational perturbation term, the perturbation objective function can be obtained: (twenty four) In the formula, This is a least-squares problem for image reconstruction. For the coefficients of the total variational disturbance term, This is the added TV perturbation value.
[0041] Based on the above perturbation objective function, a TV-CG reconstruction objective function is further formed by combining a TV regularization term, a data fidelity term, a total variation regularization term, and a total variation perturbation term. The corresponding TV-CG reconstruction objective function is expressed as follows: (25) Equation (25) is based on the TV regularization objective function of Equation (19), with a further perturbation term introduced. It was formed later. Among them, This is a total variation regularization term used to suppress noise and preserve edge features; This is a total variational perturbation term used to participate in the conjugate gradient iteration process.
[0042] This invention uses the perturbation-added conjugate gradient method combined with total variation regularization to reconstruct images. The specific algorithm of the perturbation-added conjugate gradient method is as follows.
[0043] The least squares (LS) problem of image reconstruction can be expressed as: (26) in, For the measurement data column vector, Let be the column vector of the image to be reconstructed. Let be the forward model system matrix. The LS problem is convex, and its solution satisfies a linear system: (27) Equation (27) is a widely understood method for solving the LS problem.
[0044] In image reconstruction, the LS problem is ill-conditioned and is usually addressed using regularization methods. This is achieved by adding a penalty term to the right-hand side of equation (26). (A common penalty is TV regularization) to promote normalization, so the image reconstruction after adding the penalty term is minimized as follows: (28) Equation (28) satisfies: (29) In the formula It is the gradient operator (total differential in all directions of space). These are the weighting coefficients for the total variation regularization term.
[0045] The conjugate gradient algorithm is one of the effective methods for solving equation (29). The conjugate gradient method is an algorithm based on the search direction, which uses the search direction from the previous iteration... Multiply by the conjugate search direction update coefficient and with the current negative gradient By performing linear combinations, a new search direction can be obtained. .
[0046] This invention optimizes the conjugate gradient algorithm by combining past gradients and current gradient information at a certain point, using their linear combination to construct a better search direction.
[0047] In the conjugate gradient algorithm In the parameter list, the input parameter is a matrix. sum vector Initial conjecture of the solution Allowable error Output parameters This is an approximate solution for reconstructing the image. Using... As an iteration counter exist Terminate immediately, among which According to the main standards Sure.
[0048] set up initial vector in and for: (30) (31) In each iteration, the provisional solution is: (32) In the formula, The step size is The search direction It is a positive integer starting from 0.
[0049] Rewrite equation (32) as follows In the formula A vector represents a function. , , Depending on , , .
[0050] To improve the conjugate gradient, this invention incorporates a perturbation into the CG algorithm (Algorithm 2). Although a single step is less computationally efficient than the CG algorithm, CG-PR-STEP can achieve a perturbed conjugate gradient solution for the LS problem through repeated iterations.
[0051] If for the current solution If a total variational perturbation term is added for perturbation processing, then an intermediate variable is used. This represents the current solution after perturbation processing.
[0052] If symbols are used Indicates the current solution The result after adding the total variation perturbation is:
[0053] algorithm
[0054] (1) Settings
[0055] (2) Settings
[0056] (3) Settings
[0057] (4) Settings
[0058] (5) Settings
[0059] (6) Settings
[0060] (7) Settings
[0061] (8) When
[0062] (9) Settings
[0063] (10) Call
[0064] (11) Settings
[0065] (12) Settings
[0066] algorithm
[0067] (1) Settings
[0068] (2) Settings
[0069] (3) Settings
[0070] (4) Settings
[0071] (5) Settings
[0072] (6) Settings
[0073] In Algorithm 2, Indicates the search direction of the previous iteration. Indicates the updated search direction. Update coefficients for conjugate search directions. This is the step size of the current iteration. Let be the solution for the point above the conjugate gradient direction. This is the auxiliary quantity corresponding to the previous iteration. For the updated auxiliary quantity, This is the updated solution.
[0074] The above description only illustrates the preferred embodiments of the present invention. However, the present invention is not limited to the above embodiments. Within the scope of knowledge possessed by those skilled in the art, various changes can be made without departing from the spirit of the present invention, and all such changes should be included within the protection scope of the present invention.
Claims
1. A method of photoacoustic tomographic image reconstruction based on TV-CG, characterized by, Includes the following steps: S1. A forward model for photoacoustic imaging is established based on the photoacoustic wave equation, and the system matrix is constructed using the central difference method. Thus, measured acoustic data were obtained. Compared with the initial sound pressure distribution to be reconstructed Relationship ,in, This represents the initial sound pressure column vector corresponding to each pixel within the imaging region. This is a column vector of measured acoustic data on the detection boundary; S2. Initialize the reconstruction process, setting the initial solution as... and order Simultaneously, the initial gradient vector is calculated using the conjugate gradient method. Initial search direction and corresponding auxiliary quantities ,in , Then, based on the initial step size The solution is updated for the first time. ,get And set the iteration counter. ; S3, in the iteration counter No more than the preset maximum number of iterations And the current error function value is greater than the preset threshold. Under the given conditions, a TV-CG iterative update is performed on the current solution. The update process includes: first, updating the current solution... A total variational perturbation term is introduced for perturbation processing to obtain intermediate variables. The intermediate variable The current solution after incorporating the total variational perturbation is used as the input for subsequent conjugate gradient updates; then, the intermediate variable... Using the CG-PR-STEP gradient update rule as input, perform conjugate gradient update to calculate a new solution. The corresponding search direction and auxiliary quantities; finally, let Proceed to the next iteration for judgment; S4. When the number of iterations exceeds the preset maximum number of iterations. or the current error function value Not greater than the preset threshold When the time is up, terminate the iteration and output the reconstruction result. .
2. The photoacoustic tomography image reconstruction method based on TV-CG according to claim 1, characterized in that, The TV-CG iteration in step S3 introduces a total variation regularization term into the least squares objective function. and total variational perturbation term Image reconstruction can be represented as an optimization problem: ,in, The coefficients of the total variation regularization term are... For the coefficients of the total variational disturbance term, Let be the total variational perturbation value added to the conjugate gradient iteration, and and These are mutually independent adjustable parameters.
3. The photoacoustic tomography image reconstruction method based on TV-CG according to claim 1, characterized in that, The invocation of CG-PR-STEP conjugate gradient update in step S3 is performed as follows: First, based on the current input solution... Calculate the gradient vector Combined with the search direction of the previous iteration and auxiliary amount Calculate the conjugate search direction update coefficients This leads to new search directions. And calculate the corresponding auxiliary quantities. Then, based on the step size Update the current solution to obtain Wherein, the conjugate search direction update coefficient satisfy The step size satisfy ;in, This refers to the search direction from the previous iteration. This is the auxiliary quantity corresponding to the previous iteration. For the updated search direction, For the updated auxiliary quantity, This is the updated solution.
4. The photoacoustic tomography image reconstruction method based on TV-CG according to claim 1, characterized in that, Constructing a system matrix The methods and steps are as follows: First, establish the photoacoustic wave equation, and the photoacoustic pressure. Represented as: (1) in, The coefficient of isobaric expansion is 1. For the speed of sound in the medium, The change in density For the density of the medium, and Representing spatial location and time variables respectively. Indicates temperature change; Apply equation (1) to Differentiating, we get: (2) Derivation of the sound field control equations, ； Substituting into equation (2), we get: (3) In the formula, The velocity of the particle vibration, For relative position variables The vector differential operator; applying equation (3) to Differentiate and combine get: (4) Equation (4) represents the relationship between temperature change and photoacoustic pressure in the sound field: Furthermore, temperature changes and light energy deposition The following conditions must be met: (5) in, Specific heat capacity coefficient Let be the heat source function, representing the heat per unit time and unit volume; substituting equation (5) into equation (4), we obtain the relationship between photoacoustic pressure and photoelectric energy deposition: (6) in, For relative position variables The Laplace operator, equation (6) is the photoacoustic wave equation; assuming the sound wave propagation medium is an inviscid medium and the biological tissue density is uniform, then the pressure distribution generated by the photoacoustic effect... With heat source function The relationship between them is represented as follows: Equation (6) is the photoacoustic wave equation; further, we have: (7) in, When a sound wave generated by the photoacoustic effect propagates radially at the speed of sound in the medium, at position ,time The pressure distribution formed at the location; when the heat source function satisfies: And the laser excitation satisfies: Then we have: (8) in, For position r, the duration is The initial sound pressure generated by pulsed laser excitation, and the sound pressure signal is expressed integrally. , For Grünersen parameters; ultrasonic transducer in position The sound pressure signal detected at that location, i.e., the photoacoustic signal. As a solution to equation (8), it is expressed as: (9) Location place, time The photoacoustic signal is along a radius of The time derivative of the initial sound pressure distribution of the spherical surface integral; During the reconstruction process, the three-dimensional spherical problem is simplified into a two-dimensional fault plane problem, and the integration object changes from a sphere to a circular arc. Ignoring the constants in equation (9), equation (9) is expressed in the two-dimensional fault plane as: (10) Given transducer position The position of a point on the integral curve can be represented as: (11) in, Let be the polar angle of a point on the integral curve; Equation (10) can be further written as: (12) Using the central difference method, after discretizing equation (12), we get: (13) in: (14) In the formula, The angular range covered by the transducer over the image reconstruction area; Discretize the rectangular region containing the imaging object as follows: A grid, where the position of each pixel in the grid is represented as ,in, , The corresponding energy deposition is represented as Then equation (14) can be discretized as: (15) in, For sensor position, For a point in time, for The pixel position crossed by the integral curve at time step, pixel The distance between the sensor and the sensor meets the requirements. ; Represents pixels The included angle corresponding to the intersecting arcs; combining equations (13) and (15), at time... ,Location Pressure measured at the point Represented as: (16) in, The elements in the system matrix represent the response relationship between the detection point and the pixel. pixel position Discrete values of the spatial absorption function at the location; The calculation of the system matrix includes the following steps: (1) Determine the angle The range, i.e., the calculation of the minimum angle and the maximum angle The and It relates to the first and last intersections between the integral curve and the grid; (2) Calculate the angles corresponding to each intersection point using polar coordinates. And make the obtained angle located at Within the range; (3) Obtaining various angles Then, calculate The system matrix is solved by combining equations (13), (15), and (16). ; Let the discrete regions of the reconstructed image be divided into The imaging system has a pixel grid and a transducer array element count of [number]. The length of the signal collected at each array element position is The photoacoustic signals collected at the positions of each array element and the pixel grayscale values of the reconstructed image Represented as a column vector: (17) Therefore, the relationship between the acquired photoacoustic signal and the reconstructed image can be expressed in matrix form: (18) Where, vector Let be the column vector of the image to be reconstructed. Let be the column vector of time-domain pressure signals at a given position of the ultrasonic detector. This is the forward model system matrix.