Quantitative photoacoustic tomography method and system based on PnP-ADMM

Abstract

The invention discloses a quantitative photoacoustic tomography method and system based on PnP-ADMM, and relates to the field of medical imaging, and the method comprises the steps: building a forward imaging physical model; determining a pre-trained deep denoising network; the deep denoising network is an improved DRUNet; reading a sound pressure signal measurement value; constructing a nonlinear least square optimization function; on the basis of a forward imaging physical model, a pre-trained deep denoising network and the sound pressure signal measurement value, PnP-ADMM is adopted to iteratively solve a nonlinear least square optimization function, and a light absorption coefficient distribution diagram of the biological tissue to be measured is obtained; and the PnP-ADMM adopts an L-BFGS algorithm to iteratively solve the light absorption coefficient distribution, and takes the pre-trained deep denoising network as a prior noise reduction module in the PnP to replace a near-end operator in the ADMM to iteratively solve the auxiliary variable. According to the method, high-quality and high-efficiency reconstruction of the biological tissue light absorption coefficient distribution diagram is realized.

Description

Technical Field

[0001] This application relates to the field of medical imaging technology, and in particular to a quantitative photoacoustic tomography method and system based on PnP-ADMM. Background Technology

[0002] Photoacoustic tomography (PAT) is an emerging biomedical imaging technique that combines the high absorption contrast of optical imaging with the high spatial resolution of ultrasound imaging. Its physical basis is the photoacoustic effect of biological tissues. This technique uses short-pulse lasers to irradiate biological tissues. Upon absorbing the light energy, the tissues expand elastically due to thermoelasticity, generating broadband ultrasound waves (i.e., photoacoustic waves). An ultrasound transducer collects the photoacoustic signals, and through image reconstruction, obtains a spatial distribution map of the initial sound pressure, light absorption coefficient, and chromaticity intensity within the biological tissue, reflecting its morphological structure and functional characteristics.

[0003] The inverse PAT problem comprises two aspects: first, the acoustic inverse problem, which involves reconstructing the initial sound pressure distribution or light absorption energy distribution within biological tissues based on photoacoustic signals acquired by an ultrasound detector using methods such as time reversal (TR), back projection (BP), or iterative reconstruction—essentially photoacoustic image reconstruction; and second, the optical inverse problem, which involves reconstructing functional parameters such as light absorption coefficient and chromophore concentration based on photoacoustic signal measurement data or the initial sound pressure distribution, thereby achieving quantitative imaging.

[0004] The classic method for solving the photoacoustic inverse problem is based on model-based reconstruction (MBR). The basic idea is as follows: First, an optical-acoustic coupled model is established. The optical transmission model describes the scattering and absorption of photons in biological tissue using radiative transfer equations or diffusion equations, while the acoustic propagation model characterizes the generation and propagation mechanism of photoacoustic waves based on photoacoustic wave equations. Then, an iterative algorithm is used to solve an optimization problem that minimizes the error between the measured and simulated values ​​of the photoacoustic signal, thus enabling the estimation of tissue characteristics and functional parameters. However, factors such as computational errors introduced by the discretization of the forward imaging model, ill-posedness of the inverse problem due to multi-parameter coupling, the rationality of the regularization strategy and its weight parameters, and the efficiency bottleneck of large-scale computation collectively limit the accuracy and efficiency of the MBR method in achieving quantitative photoacoustic imaging. Summary of the Invention

[0005] The purpose of this application is to provide a quantitative photoacoustic tomography method and system based on PnP-ADMM, which achieves high-quality and high-efficiency reconstruction of the light absorption coefficient distribution map of biological tissues.

[0006] To achieve the above objectives, this application provides the following solution.

[0007] Firstly, this application provides a quantitative photoacoustic tomography method based on PnP-ADMM. The method includes: establishing a forward imaging physical model; the forward imaging physical model simulating the physical process from photon absorption by biological tissue to sound pressure signal detection; determining a pre-trained deep denoising network; the deep denoising network being an improved DRUNet; the improved DRUNet having no bias terms in all convolutional layers, no activation functions after the first and last convolutional layers, as well as the SConv and TConv layers, and each residual block containing only one ReLU activation function; and reading sound pressure signal measurements; the sound pressure signal measurements are taken during short pulses. The detection process involves laser irradiation of the biological tissue under test; a nonlinear least squares optimization function is constructed; the nonlinear least squares optimization function characterizes the process of estimating the light absorption coefficient distribution of the biological tissue under test based on the measured sound pressure signal; based on the forward imaging physical model, the pre-trained deep denoising network, and the measured sound pressure signal, the nonlinear least squares optimization function is iteratively solved using PnP-ADMM to obtain the light absorption coefficient distribution map of the biological tissue under test; the PnP-ADMM uses the L-BFGS algorithm to iteratively solve the light absorption coefficient distribution, and uses the pre-trained deep denoising network as the prior denoising module in PnP, replacing the proximal operator in ADMM, to iteratively solve the auxiliary variables.

[0008] Secondly, this application also provides a quantitative photoacoustic tomography system based on PnP-ADMM, the PnP-ADMM-based quantitative photoacoustic tomography system comprising: a laser, an ultrasound detector, a processor, and a display; the ultrasound detector and the display are respectively connected to the processor; the laser is used to generate short-pulse laser to irradiate the biological tissue to be tested; the short-pulse laser is a laser with a pulse width on the order of nanoseconds or less; the processor is used to execute the PnP-ADMM-based quantitative photoacoustic tomography method described in the first aspect; the display is used to display the light absorption coefficient distribution map of the biological tissue to be tested.

[0009] Based on the specific embodiments provided in this application, the following technical effects are disclosed.

[0010] First, this application, based on a forward imaging physical model, employs a plug-and-play alternating direction method of multipliers (PnP-ADMM) framework to organically combine the numerical solution of the physical model with deep learning-based regularization prior processing. Second, in the iterative solution of the nonlinear least squares optimization function (i.e., the optical-acoustic coupling inverse problem), this application first utilizes the limited-memory BFGS (L-BFGS) algorithm to iteratively solve the light absorption coefficient distribution, overcoming the difficulty of solving nonlinear and ill-posed problems. The L-BFGS algorithm significantly reduces storage requirements and is particularly suitable for processing large-scale data. Third, this application also implicitly processes the regularization term through a pre-trained deep denoising network. Improvements in the deep denoising network structure enable it to effectively handle denoising tasks with different noise levels using a single model, providing powerful feature representation capabilities while maintaining the integrity of the image structure, making it particularly suitable as a regularization prior module for photoacoustic image reconstruction. Through the above-mentioned technical features, this application achieves high-quality and high-efficiency reconstruction of the light absorption coefficient distribution map of biological tissues. Attached Figure Description

[0011] To more clearly illustrate the technical solutions in the embodiments of this application or related technologies, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0012] Figure 1 This is a flowchart of a quantitative photoacoustic tomography method based on PnP-ADMM in one embodiment of this application.

[0013] Figure 2 This is a schematic diagram of the structure of a deep denoising network in one embodiment of this application.

[0014] Figure 3 This is a flowchart illustrating the training process of a deep denoising network in one embodiment of this application.

[0015] Figure 4 This is a flowchart illustrating the execution of a trained deep denoising network in one embodiment of this application.

[0016] Figure 5 This is a schematic diagram illustrating the solution of PnP-ADMM in one embodiment of this application.

[0017] Figure 6 This is a flowchart of the solution process for PnP-ADMM in one embodiment of this application.

[0018] Figure 7 This is a schematic diagram of the layout of a quantitative photoacoustic tomography system based on PnP-ADMM in another embodiment of this application.

[0019] Figure labels: Laser-1, Ultrasonic detector-2, Processor-3, Display-4, Biological tissue to be tested-5. Detailed Implementation

[0020] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0021] The classic method for solving the optical inverse problem (PAI) is the MBR method, which achieves quantitative reconstruction of tissue characteristics and functional parameters through numerical inversion of the forward imaging model. By embedding prior physical constraints into the objective function, the ill-conditioned nature of the problem can be alleviated, ensuring a globally optimal solution. However, the accuracy and efficiency of this method are mainly constrained by three factors: First, the discretization process of the forward imaging model introduces errors due to traditional interpolation methods, making it difficult to accurately describe the ultrasound propagation path. Second, the conventional gradient descent method used in the optimization process is prone to getting trapped in local minima and converges slowly, making it difficult to stably approach the global optimum when the solution space is large. Finally, in terms of regularization, both the Tikhonov strategy, which over-smooths high-frequency details, and the total variational method, which preserves edges but smooths weak gradient regions, have obvious limitations. While Bayesian framework-based methods are more flexible, their Gaussian priors do not match the discontinuous characteristics of biological tissues, and they suffer from problems such as the inability to adaptively preset hyperparameters and high computational costs. Therefore, solving the PAI inverse problem using the above-mentioned classic methods cannot achieve high-quality and high-efficiency reconstruction of quantitative optical imaging of biological tissues.

[0022] In recent years, deep learning has demonstrated strong application potential in the field of medical imaging. By automatically learning multi-level feature representations through hierarchical networks, it overcomes the limitations of traditional methods that rely on manually designed feature extraction operators, significantly improving imaging quality. Combining deep learning with forward physics models, leveraging the powerful feature extraction capabilities of neural networks while incorporating the physical constraints of imaging, can effectively suppress artifacts and noise in reconstructed images. A typical method is the PnP framework, which decouples the inverse problem into data fidelity terms and deep learning priors, achieving excellent reconstruction performance while maintaining physical consistency, thus providing a new technical path for image reconstruction.

[0023] The purpose of this application is to provide a quantitative photoacoustic tomography method and system based on PnP-ADMM, which achieves high-quality and high-efficiency reconstruction of the light absorption coefficient distribution map of biological tissues.

[0024] To make the above-mentioned objectives, features and advantages of this application more apparent and understandable, the application will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0025] In one exemplary embodiment, a quantitative photoacoustic tomography method based on PnP-ADMM is provided, such as... Figure 1 As shown, the quantitative photoacoustic tomography method based on PnP-ADMM includes the following steps.

[0026] Step S1: Establish a forward imaging physical model.

[0027] In this embodiment, the forward imaging physical model consists of a coupled optical transmission model and an acoustic propagation model. The optical transmission model outputs the distribution of light absorption energy within the biological tissue, while the acoustic propagation model simulates the propagation of ultrasonic waves generated therefrom and obtains the corresponding sound pressure signal. The two models are coupled through the thermoelastic expansion effect, fully describing the physical process from photon absorption by biological tissue to sound pressure signal detection.

[0028] (1) Optical transmission model.

[0029] In optical imaging technology, the radiative transfer equation (RTE) is usually used to describe the transmission process of incident light in a turbid medium. The governing equation is as follows.

[0030] .

[0031] In the formula, Indicates position Along the direction The light radiation energy per unit area This is the set of all possible directions of photon propagation. and Representing positions respectively The light absorption coefficient and light scattering coefficient at that location; Let be the scattering phase function, which is a probability density function describing the scattering phase along the path. Photons propagating in a specific direction are scattered to The probability of direction; This refers to the light source item.

[0032] Because directly solving the radiative transport equation is extremely complex, diffusion approximation (DA) or Monte Carlo (MC) methods are generally used for approximation. DA is a deterministic model that is only effective in the diffusion domain (i.e., the region from the light source to the mean free path of several photons). However, DA performs poorly in regions near the light source and at the boundaries of the imaging region, which often constitute an important part of photoacoustic imaging and contain crucial information needed to assess the structure and functional components of biological tissues. In contrast, the MC method is a stochastic method that can accurately simulate the interaction between light and biological tissues. The MC method has the advantage of parallel computation and can more accurately describe the photon transport behavior in biological tissues. This embodiment uses the MC method to approximate the light transport equation, estimating the light radiant flux by tracking the propagation paths of a large number of individual photons in the biological tissue and counting the number of photons at each location.

[0033] .

[0034] In the formula, Indicates position The radiant flux at a given location. The absorbed energy density of light. With light radiation flux and light absorption coefficient The product is directly proportional.

[0035] .

[0036] (2) Sound propagation model.

[0037] When biological tissue absorbs light energy, it causes a local temperature increase, which in turn induces thermoelastic expansion and generates initial sound pressure. This process can be described by the thermoelastic effect equation.

[0038] .

[0039] In the formula, Indicates position The initial sound pressure at that location. Here are the Grüneisen parameters. The propagation of sound waves in a medium is described by the following wave equation.

[0040] .

[0041] In the formula, Indicates position Place t Sound pressure at any given moment; c For the speed of sound, The coefficient of thermal expansion is the isobaric coefficient. Specific heat capacity under a certain pressure; Let be the time-domain function of the incident laser pulse. The initial conditions for the above equation are as follows.

[0042] .

[0043] In this embodiment, the k-space pseudospectral method (PSM) is used to solve the above wave equation to obtain the sound pressure signal received by the ultrasonic detector.

[0044] (3) Operator representation of the forward imaging physical model.

[0045] Given the light scattering coefficient and sound velocity of a biological tissue, the forward imaging physical model can be described by nonlinear operators.

[0046] .

[0047] In the formula, This represents the spatial distribution of light absorption coefficients in the imaging plane. This represents the sound pressure signal matrix in the imaging plane. This represents the mapping from the light absorption coefficient to the sound pressure signal.

[0048] Step S2: Determine the pre-trained deep denoising network.

[0049] In this embodiment, the deep denoising network is an improved Deep Residual U-Net (DRUNet). Figure 2 As shown, the deep denoising network is a 15-layer network structure, based on the four-scale encoding and decoding structure of U-Net. Each scale includes 2×2 stride convolution (SConv) downsampling and 2×2 transposed convolution (TConv) upsampling operations, and features identity skip connections. From the first to the fourth scale, the number of network channels is 64, 128, 256, and 512, respectively. Each scale is configured with four consecutive residual blocks to enhance feature extraction capabilities. Regarding activation function design, referencing the super-resolution network architecture, no activation functions are set after the first and last convolutional layers, as well as after the SConv and TConv layers. Each residual block contains only one ReLU activation function. Furthermore, no bias terms are set in any convolutional layer (including regular Conv, SConv, and TConv). This design is based on two main considerations: First, the unbiased network combining ReLU activation and identity skip connections can enhance the scaling invariance of the image restoration task, i.e., satisfy... For any scalar a≥0 is true; secondly, in traditional biased networks, the magnitude of the bias term is often significantly larger than the filter weights, and this imbalance may weaken the network's generalization performance. Through the above network architecture design, the deep denoising network can effectively handle denoising tasks with different noise levels with a single model, providing powerful feature representation capabilities while maintaining the integrity of the image structure, making it particularly suitable as a regularization prior module for photoacoustic image reconstruction.

[0050] In this embodiment, step S2 is as follows.

[0051] Step S21: Construct a numerical biomimetic and slice the target part of the numerical biomimetic at equal intervals to obtain multiple slices.

[0052] As a preferred embodiment, a numerical mouse phantom was constructed, and its chest and abdomen were sliced ​​at equal intervals of 0.3 mm to ensure consistency with in vivo imaging parameters. Each slice was 256×256 in size and included different tissue components, such as the thoracic aorta, spinal muscles, and spinal cord.

[0053] Step S22: Set the optical and acoustic properties of different tissue components.

[0054] The tissue characteristic parameters that need to be set include optical characteristic parameters and acoustic characteristic parameters. Optical characteristic parameters include light absorption coefficient, light scattering coefficient, anisotropy factor, etc., while acoustic characteristic parameters include refractive index, sound velocity, density, etc.

[0055] Step S23: Use the true value distribution map of light absorption coefficients corresponding to different tissue components as a clean image.

[0056] Step S24: Add Gaussian noise with a mean of 0 and a standard deviation that is uniformly distributed in the interval [0, 50] to the clean image to generate noise images with different noise levels.

[0057] The above noise intensity range corresponds to a signal-to-noise ratio (SNR) of approximately 4.7 dB (the lowest calculated SNR) to 40.3 dB (the highest calculated SNR).

[0058] Step S25: Pair each noisy image with its corresponding clean image to construct a sample set.

[0059] In one preferred embodiment, the sample set contains at least 200 pairs of samples (i.e., image pairs).

[0060] Step S26: Divide the sample set into training set, validation set and test set according to the proportion.

[0061] As a preferred implementation, multiple pairs of samples in the sample set are randomly shuffled and divided into a training set, a validation set, and a test set in a ratio of 8:1:1.

[0062] Step S27: Use data augmentation techniques to expand the training set.

[0063] As a preferred implementation, in order to avoid overfitting, data augmentation techniques (including random flipping, rotation and translation) are used to expand the training set to 2000 samples.

[0064] Step S28: Based on the expanded training set, with noisy images as input and clean images as output, train the improved DRUNet to obtain a pre-trained deep denoising network.

[0065] As a preferred implementation, training the deep denoising network also requires the use of the Adam optimizer, which is used to adjust the network weights to minimize the loss function. For example... Figure 3 As shown, the specific training process is as follows.

[0066] (1) The initial number of times to traverse all mini-batch training sets is epoch=0.

[0067] (2) Initialize the index of the training set i =1.

[0068] (3) The first training set i In small batch datasets N Each sample is input into the network, and forward propagation is performed to obtain the network output.

[0069] (4) Compare the image output by the network with the clean image (target) and calculate the loss value.

[0070] Define a loss function. L 1. Loss is used to measure the difference between the denoised image and the clean image output by the network.

[0071] .

[0072] In the formula, The image output by the network. For a clean image.

[0073] (5) Set the gradient of each layer parameter in the network to 0.

[0074] (6) Calculate the gradient using the backpropagation algorithm and update the network weights.

[0075] (7) Optimize the network parameters using the Adam gradient descent algorithm, setting the learning rate to 1×10. ‒4 .

[0076] (8) Order i ← i +1. If i < N If , then return (3); otherwise, let epoch←epoch+1 and execute (9).

[0077] (9) If epoch < 500, return to (2); otherwise, stop traversing and the training process ends.

[0078] Step S3: Read the sound pressure signal measurement value.

[0079] In this embodiment, a short-pulse laser (i.e., a laser with a pulse width on the order of nanoseconds or less) is first used to irradiate the biological tissue to be tested. Simultaneously, an ultrasonic detector detects the corresponding sound pressure signal. Since this sound pressure signal is obtained based on the biological tissue to be tested, it is called the sound pressure signal measurement value. To facilitate the subsequent determination of the light absorption coefficient distribution, it is necessary to read this sound pressure signal measurement value from the ultrasonic detector.

[0080] Step S4: Construct a nonlinear least squares optimization function.

[0081] In this embodiment, the process of estimating the light absorption coefficient distribution of the biological tissue under test based on the sound pressure signal measurement can be expressed as the following nonlinear least squares optimization problem (function).

[0082] .

[0083] .

[0084] In the formula, and These are the estimated and optimized values ​​of the light absorption coefficient distribution, respectively. For regularization parameters; As a regularization term, a convex function is usually chosen to ensure the convergence and robustness of the solution; For data fidelity; This represents the mapping from the light absorption coefficient to the sound pressure signal; This is a matrix of sound pressure signal measurement values.

[0085] Step S5: Based on the forward imaging physical model, the pre-trained deep denoising network, and the measured sound pressure signal, the nonlinear least squares optimization function is solved iteratively using PnP-ADMM to obtain the light absorption coefficient distribution map of the biological tissue to be tested.

[0086] In this embodiment, the PnP-ADMM algorithm is used to solve the above optimization problem. This algorithm combines the flexibility of modern regularization techniques with the high efficiency of the ADMM algorithm, enabling it to obtain high-quality results in image reconstruction tasks. The specific implementation steps are as follows: First, auxiliary variables are introduced. The original optimization problem is reformulated as a constrained optimization problem (i.e., objective function) as follows.

[0087] .

[0088] Its corresponding augmented Lagrangian function is as follows.

[0089] .

[0090] In the formula, To augment the Lagrange function; These are dual variables, i.e., Lagrange multipliers; is the penalty coefficient. Then, the following subproblems are solved iteratively, alternating between them, until convergence.

[0091] .

[0092] .

[0093] .

[0094] In the formula, and The first k The second iteration and the first k +1 iterations of light absorption coefficient distribution; and The first k The second iteration and the first k Auxiliary variables for +1 iterations; and The first k The second iteration and the first k +1 iterations of Lagrange multiplier.

[0095] Because the aforementioned constrained optimization problem is nonlinear and ill-posed, the subproblem concerning the distribution of the light absorption coefficient is difficult to solve analytically. Therefore, this embodiment employs the L-BFGS algorithm for an approximate solution. L-BFGS is a widely used quasi-Newton algorithm in machine learning; it solves by only storing the most recent... mThe iterations significantly reduce storage requirements, making it particularly suitable for processing large-scale data. In the implementation, a line search strategy is used to adaptively determine the step size to accelerate convergence. A strong Wolfe condition is employed to ensure the convergence and efficiency of the step size, and a caching mechanism is implemented to reduce redundant computations of forward operators and gradients. The gradient of the data fidelity term in the objective function is solved using the adjoint method.

[0096] .

[0097] In the formula, for The adjoint operator.

[0098] The subproblem concerning the auxiliary variables is only related to the regularization term and can be considered as a regularized image denoising problem. Therefore, this embodiment uses a pre-trained deep denoising network as the prior denoising module in the PnP framework, replacing the proximal operator in the traditional ADMM algorithm. Figure 4 As shown, the specific execution steps of the pre-trained deep denoising network are as follows.

[0099] (1) Initialization parameters: initial number of feature channels Channels=2, number of convolutions Conv=0, number of stride convolutions SConv=0, number of transpose convolutions TConv=0.

[0100] (2) Remove the elements in the objective function that do not meet the stopping iteration requirement. As input, its size is S=H×W×1, where H is the length, W is the width, and the number of image channels is 1.

[0101] (3) Perform a convolution operation with a kernel size of 3×3, a stride of 1, and padding of 1 on the input of step (2) through 2 feature channels, and output feature channels = 64.

[0102] (4) Use residual blocks to construct skip connections in the network by directly adding the input and output. First, perform a convolution operation with channels, kernel size of 3×3, stride of 1, and padding of 1. Then, apply the ReLU activation function to achieve non-linear transformation, perform another convolution operation, and finally add the input feature map element by element to the output after convolution and activation.

[0103] (5) Perform stride convolution downsampling operation on the feature map obtained in step (4) through the channel(s) feature channels with a kernel size of 2×2, a stride of 2, and padding of 0. The output channel is channel(s)←2×channels. At the same time, the stride convolution and the transposed convolution at the same scale have an identity jump connection.

[0104] (6) Let SConv←SConv+1.

[0105] (7) If SConv < 3, return to step (4); otherwise, use residual connections to further extract feature maps.

[0106] (8) Perform transposed convolution upsampling operation on the feature map obtained in step (7) with a kernel size of 2×2, a stride of 2, and padding of 0 through the channel 'channels' feature channels. The output channel is channels←channels / 2.

[0107] (9) Let TConv ← TConv+1.

[0108] (10) Use residual blocks to construct the skip connections of the network by directly adding the input and output, where the operation details are the same as in step (4).

[0109] (11) If TConv < 3, return to step (8); otherwise, perform a convolution operation on the obtained feature map with a kernel size of 3×3, a stride of 1, and padding of 1.

[0110] (12) The number of channels becomes 1, and the output is a denoised image.

[0111] like Figure 5 As shown, step S5 in this embodiment can be summarized as follows.

[0112] Step S51: Determine the estimated value of the light absorption coefficient distribution and initialize the auxiliary and dual variables (Lagrange multipliers).

[0113] Step S52: Based on the estimated value of the light absorption coefficient distribution, the predicted value of the sound pressure signal is obtained using the forward imaging physical model.

[0114] Step S53: Based on the predicted sound pressure signal value, the measured sound pressure signal value, auxiliary variables and dual variables, update the light absorption coefficient distribution using the L-BFGS algorithm, and determine whether the updated light absorption coefficient distribution meets the stopping iteration requirements (i.e., the gradient is small enough or the number of iterations is greater than or equal to the set number of iterations).

[0115] If so, the updated light absorption coefficient distribution is output directly, and its visualization yields the light absorption coefficient distribution map of the biological tissue under test.

[0116] If not, the iterative steps are executed until the updated optical absorption coefficient distribution meets the stopping iteration requirement. The iterative steps include: updating the predicted sound pressure signal value using the forward imaging physical model based on the updated optical absorption coefficient distribution; updating the auxiliary variables using a pre-trained deep denoising network based on the updated optical absorption coefficient distribution; updating the dual variable based on the updated optical absorption coefficient distribution and the updated auxiliary variables; and updating the optical absorption coefficient distribution again using the L-BFGS algorithm based on the updated predicted sound pressure signal value, the measured sound pressure signal value, the updated auxiliary variables, and the updated dual variables.

[0117] like Figure 6 As shown, the specific process of reading the sound pressure signal measurement value and reconstructing the light absorption coefficient distribution map of the biological tissue under test in this embodiment is as follows.

[0118] (1) Input data: The guessed value of the light absorption coefficient distribution and the measured value of the collected sound pressure signal are used as input.

[0119] (2) Initialize variables: Initialize auxiliary variables v 0 and dual variables (Lagrange multipliers) w 0. These variables are used to break down complex problems into "data fidelity" problems and "regularization" problems.

[0120] (3) Calculate the gradient: based on the first k The light absorption coefficient distribution of the next iteration And its relationship with The gradient of the objective function is calculated from the residual.

[0121] (4) Double loop recursion to determine the search direction: The algorithm uses historical gradient information to quickly estimate the product of the inverse of the Hessian matrix and the gradient through double loop recursion, thereby obtaining the descent direction.

[0122] (5) Line search to determine step size: After determining the direction, a suitable step size is found through line search to ensure that the objective function value decreases. Strong Wolfe conditions are used to ensure the convergence and efficiency of the step size, and a caching mechanism is set to reduce the repeated calculation of forward operators and gradients.

[0123] (6) Update : Combine direction and step size to update .

[0124] (7) Convergence judgment: check If the iteration stops (gradient is small enough or number of iterations is greater than or equal to 20), the loop is exited and the output phase is entered; otherwise, the denoising / regularization phase is entered.

[0125] (8) A pre-trained deep denoising network is used to replace the proximal operator in the ADMM algorithm to solve the regularization problem, and image features learned from big data are used to suppress noise and artifacts. Combined with auxiliary variables, the input is fed into a pre-trained deep denoising network, and the output of this network becomes the updated auxiliary variables. v k .

[0126] (9) Utilization Update dual variables w k Used for coordination and v Consistency between them.

[0127] (10) Utilization w k and v k Start the next round Iteration.

[0128] (11) After the loop ends, the final This is the reconstructed light absorption coefficient distribution (matrix). Visualizing it yields the light absorption coefficient distribution map of the biological tissue under test, enabling quantitative imaging.

[0129] In one exemplary embodiment, a quantitative photoacoustic tomography system based on PnP-ADMM is provided, such as... Figure 7 As shown, the quantitative photoacoustic tomography system based on PnP-ADMM includes: a laser 1, an ultrasonic detector 2, a processor 3, and a display 4.

[0130] Specifically, the ultrasound detector 2 and the display 4 are respectively connected to the processor 3; the laser 1 is used to generate short-pulse laser to irradiate the biological tissue 5 to be tested; the processor 3 is used to execute the above-mentioned quantitative photoacoustic tomography method based on PnP-ADMM; and the display 4 is used to display the light absorption coefficient distribution map of the biological tissue to be tested.

[0131] In a preferred embodiment, laser 1 is a yttrium aluminum garnet (YAG) laser, capable of generating tunable laser light across the entire wavelength range (680nm~900nm), with a repetition rate of 10Hz, a pulse duration of 7ns, and a maximum incident pulse energy of 120mJ. The ultrasonic detector 2 consists of 256 focusing transducer elements with a center frequency of 5MHz and a bandwidth of 60%. The processor 3 can be a non-volatile computer-readable storage medium, and the computer program, when executed, can include the flow described in the above method embodiments.

[0132] In summary, this application constructs a complete photoacoustic coupling forward model (i.e., a forward imaging physical model) to accurately describe the generation and propagation of photoacoustic signals. During the ADMM iteration process, a deep learning denoiser (i.e., a pre-trained deep denoising network) is introduced as a regularization subproblem solver to enhance the algorithm's stability and noise resistance. Efficient convergence is achieved through alternating optimization. This application significantly outperforms mainstream comparative methods such as proximal gradient descent networks, deep gradient descent networks, extractor-attention-predictor networks, Tikhonov regularization, and total variational directional analysis in both visual quality and quantitative evaluation metrics. Therefore, this application achieves a good balance between accuracy and efficiency, providing a reliable approach for quantitative photoacoustic imaging. Furthermore, its framework integrating physical models and prior data offers valuable insights for other imaging techniques related to ill-conditioned inverse problems.

[0133] Any references to memory, databases, or other media used in the embodiments provided in this application may include at least one of non-volatile and volatile memory. Non-volatile memory may include read-only memory (ROM), magnetic tape, floppy disk, flash memory, optical memory, high-density embedded non-volatile memory, resistive random access memory (ReRAM), magnetic random access memory (MRAM), ferroelectric random access memory (FRAM), phase change memory (PCM), graphene memory, etc. Volatile memory may include random access memory (RAM) or external cache memory, etc. By way of illustration and not limitation, RAM may be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0134] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. The same or similar parts between the various embodiments can be referred to each other.

[0135] This document uses specific examples to illustrate the principles and implementation methods of this application. The descriptions of the above embodiments are only for the purpose of helping to understand the methods and core ideas of this application. Furthermore, those skilled in the art will recognize that, based on the ideas of this application, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of this application.

Claims

1. A PnP-ADMM based quantitative photoacoustic tomography method, characterized in that, The quantitative photoacoustic tomography method based on PnP-ADMM includes: A forward imaging physical model was established; the forward imaging physical model simulates the physical process from the absorption of photons by biological tissue to the detection of acoustic pressure signals; A pre-trained deep denoising network is determined; the deep denoising network is an improved DRUNet; no bias terms are set in all convolutional layers of the improved DRUNet, no activation functions are set after the first and last convolutional layers, as well as the SConv and TConv layers, and each residual block contains only one ReLU activation function. Read the sound pressure signal measurement value; the sound pressure signal measurement value is detected after the biological tissue to be tested is irradiated with a short pulse laser; A nonlinear least squares optimization function is constructed; the nonlinear least squares optimization function characterizes the process of estimating the light absorption coefficient distribution of the biological tissue under test based on the sound pressure signal measurement value. Based on the aforementioned forward imaging physical model, the pre-trained deep denoising network, and the measured acoustic pressure signal, the nonlinear least squares optimization function is iteratively solved using PnP-ADMM to obtain the light absorption coefficient distribution map of the biological tissue under test. The PnP-ADMM uses the L-BFGS algorithm to iteratively solve the light absorption coefficient distribution, and uses the pre-trained deep denoising network as the prior denoising module in PnP to replace the proximal operator in ADMM to iteratively solve the auxiliary variables.

2. The PnP-ADMM based quantitative photoacoustic tomography method of claim 1, wherein, The operator representation of the forward imaging physical model is as follows: ； wherein is the light absorption coefficient distribution; is the sound pressure signal matrix; denotes the mapping from the light absorption coefficient to the sound pressure signal.

3. The PnP-ADMM based quantitative photoacoustic tomography method of claim 1, wherein, The forward imaging physical model consists of a coupled optical transmission model and an acoustic propagation model. The optical transmission model uses the Monte Carlo method to determine the distribution of light absorption energy inside biological tissue. The acoustic propagation model uses the k-space pseudospectral method to determine the detected sound pressure signal. The optical transmission model and the acoustic propagation model are coupled through the thermoelastic expansion effect.

4. The PnP-ADMM based quantitative photoacoustic tomography method of claim 1, wherein, The determination of the pre-trained deep denoising network specifically includes: A numerical biomimetic is constructed, and the target part of the numerical biomimetic is sliced ​​at equal intervals to obtain multiple slices; each slice includes different tissue components. The optical and acoustic properties of different tissue components are set; the optical properties include at least the light absorption coefficient. The true value distribution map of the light absorption coefficient corresponding to different tissue components is used as a clean image; Gaussian noise with a mean of 0 and a standard deviation that is uniformly taken in the interval [0, 50] is added to the clean image to generate noise images with different noise levels. Each noisy image is paired with its corresponding clean image to form an image pair, thus constructing a sample set. The sample set is divided into a training set, a validation set, and a test set according to a certain ratio; Data augmentation techniques are used to augment the training set. Based on the expanded training set, the noisy image is used as input and the clean image is used as output to train the improved DRUNet, resulting in a pre-trained deep denoising network.

5. The PnP-ADMM based quantitative photoacoustic tomography method of claim 1, wherein, The improved DRUNet is based on the four-scale encoding and decoding structure of U-Net. Each scale includes 2×2 stride convolution downsampling and 2×2 transposed convolution upsampling operations, and has identity skip connections. From the first scale to the fourth scale, the number of network channels are 64, 128, 256 and 512 respectively. Each scale is configured with 4 consecutive residual blocks.

6. The quantitative photoacoustic tomography method based on PnP-ADMM according to claim 1, characterized in that, The nonlinear least squares optimization function is: ； ； In the formula, and These are the estimated and optimized values ​​of the light absorption coefficient distribution, respectively. For regularization parameters; For regularization terms; For data fidelity; This represents the mapping from the light absorption coefficient to the sound pressure signal; This is a matrix of sound pressure signal measurement values.

7. The quantitative photoacoustic tomography method based on PnP-ADMM according to claim 1, characterized in that, Based on the aforementioned forward imaging physical model, the pre-trained deep denoising network, and the measured acoustic pressure signal, the nonlinear least squares optimization function is solved iteratively using PnP-ADMM to obtain the light absorption coefficient distribution map of the biological tissue under test, specifically including: Determine the estimated value of the light absorption coefficient distribution and initialize the auxiliary and dual variables; Based on the estimated value of the light absorption coefficient distribution, the predicted value of the sound pressure signal is obtained using the forward imaging physical model; Based on the predicted sound pressure signal value, the measured sound pressure signal value, the auxiliary variable, and the dual variable, the light absorption coefficient distribution is updated using the L-BFGS algorithm, and it is determined whether the updated light absorption coefficient distribution meets the stopping iteration requirement. If so, the updated light absorption coefficient distribution will be output directly, and the light absorption coefficient distribution map of the biological tissue to be tested will be obtained after visualization. If not, then perform a loop iteration step until the updated light absorption coefficient distribution meets the stopping iteration requirement; the loop iteration step includes: updating the predicted sound pressure signal value using the forward imaging physical model based on the updated light absorption coefficient distribution; updating the auxiliary variable using a pre-trained deep denoising network based on the updated light absorption coefficient distribution; updating the dual variable based on the updated light absorption coefficient distribution and the updated auxiliary variable; and updating the light absorption coefficient distribution again using the L-BFGS algorithm based on the updated predicted sound pressure signal value, the measured sound pressure signal value, the updated auxiliary variable, and the updated dual variable.

8. The quantitative photoacoustic tomography method based on PnP-ADMM according to claim 6, characterized in that, The iterative solution for the light absorption coefficient distribution is expressed as follows: ； In the formula, and The first k The second iteration and the first k +1 iteration light absorption coefficient distribution; This is the penalty coefficient; For the first k Auxiliary variables for the next iteration; For the first k The dual variable of the next iteration; The iterative solution of the auxiliary variable is expressed as follows: ； In the formula, For the first k Auxiliary variables for +1 iterations.

9. A quantitative photoacoustic tomography system based on PnP-ADMM, characterized in that, The quantitative photoacoustic tomography system based on PnP-ADMM includes: a laser, an ultrasound detector, a processor, and a display; the ultrasound detector and the display are respectively connected to the processor; The laser is used to generate short-pulse laser light to irradiate the biological tissue to be tested. The processor is used to execute the quantitative photoacoustic tomography method based on PnP-ADMM as described in any one of claims 1-8; The display is used to show the light absorption coefficient distribution map of the biological tissue to be tested.

10. The quantitative photoacoustic tomography system based on PnP-ADMM according to claim 9, characterized in that, The laser is a YAG laser; the ultrasonic detector consists of 256 focusing transducer elements.

Images