Method of simultaneous reconstruction of acoustic pressure and acoustic velocity using multi-grid iteration in photoacoustic tomography

Abstract

The application discloses a sound pressure-sound velocity synchronous reconstruction method of multi-grid iteration in photoacoustic tomography, and relates to the technical field of medical image reconstruction, and comprises the following steps: constructing a joint reconstruction least square problem by minimizing the two norms of a real detector detection signal and a simulation detection signal, and dividing the joint reconstruction least square problem into two problems for solving; problem one is to solve a sound pressure distribution by a known sound velocity distribution, and problem two is to calculate a sound velocity distribution by taking the sound pressure distribution obtained by problem one as a known sound pressure distribution, and the obtained sound pressure distribution and sound velocity distribution are taken as initial values and cycled to problem one for solving; multi-grid iteration reconstruction: joint reconstruction iteration is carried out on coarse and fine grids, and finally, the output of the fine grid is taken as the reconstructed sound pressure distribution and the reconstructed sound velocity distribution; by joint reconstruction under different scale grids, the method can simultaneously reconstruct the initial sound pressure distribution and the sound velocity distribution from photoacoustic signals, improve image quality, and accelerate the joint reconstruction speed.

Description

Technical Field

[0001] The present invention relates to the technical field of medical image reconstruction, and in particular to a multi-grid iterative sound pressure-speed synchronous reconstruction method in photoacoustic tomography. Background Art

[0002] Photoacoustic Computed Tomography (PACT) is an emerging biomedical imaging technology in which the object to be imaged is irradiated with short laser pulses. The biological tissue within the object absorbs the light energy, causing local heating, which in turn induces thermoelastic expansion and produces an initial acoustic pressure distribution. The pressure wave propagates outward and is received by ultrasonic transducers placed around the object. The received signal is processed using photoacoustic tomography image reconstruction methods, ultimately determining the initial acoustic pressure distribution of the tissue.

[0003] Existing photoacoustic image reconstruction methods typically assume that the object to be imaged has the same sound velocity as the surrounding background. However, the sound velocity distribution in biological tissue is often heterogeneous. Traditional methods, such as filtered back projection algorithms, require adjusting a single sound velocity to improve imaging quality in a specific area, or simply perform dual sound velocity estimation, estimating the sound velocity of the object to be imaged and the background separately to improve imaging quality. However, when the sound velocity distribution in biological tissue is highly heterogeneous, these methods have limited improvement in image quality. Another class of methods, such as time reversal and model-based iterative methods, can obtain higher-quality images using the known object sound velocity distribution. However, obtaining the object sound velocity distribution requires combining it with ultrasonic tomography methods, which requires the detector to have additional emission capabilities on the equipment and requires additional computational time, increasing both economic and time costs.

[0004] Synchronously reconstructing the initial sound pressure distribution and sound speed distribution from the detection signal is a problem that needs to be broken through. This problem is called the joint reconstruction problem. The main drawback of the current joint reconstruction method is that the reconstruction problem is highly ill-posed and the iterative process will fall into a local optimal solution. Some researchers have proposed using a small amount of prior knowledge of sound speed to segment the sound speed distribution and then reconstruct it, but this still cannot get rid of the dependence on prior knowledge. A method called feature coupling has been proposed to solve this type of problem. It divides the sound speed area through prior segmentation, adjusts the sound speed in different areas, and compares the coherence of the images reconstructed by the two semi-ring detectors to evaluate the correctness of the current sound speed. However, this method is only applicable to cases with small heterogeneity. Summary of the Invention

[0005] Based on the technical problems existing in the background technology, the present invention proposes a multi-grid iterative sound pressure-sound velocity synchronous reconstruction method in photoacoustic tomography. By performing joint reconstruction on grids of different scales, the initial sound pressure distribution and sound velocity distribution can be simultaneously reconstructed from only the photoacoustic signal, thereby improving image quality and accelerating the joint reconstruction speed.

[0006] The method for simultaneously reconstructing sound pressure and sound velocity in photoacoustic tomography by using multi-grid iteration, comprises the following steps:

[0007] Step one, constructing a joint reconstruction least square problem by minimizing the two-norm of the real detection signal of the detector and the simulated detection signal, and solving the problem in two sub-problems, wherein the simulated detection signal is a signal obtained by solving a wave equation by using the sound pressure distribution and the sound velocity distribution as initial values;

[0008] Step two, solving the sound pressure distribution by using the conjugate gradient method in the first sub-problem, and using the sound pressure distribution obtained in the first sub-problem as the known sound pressure distribution to construct the two-stage mass source term in the two-stage adjoint field calculation in the second sub-problem, so as to obtain the sound velocity distribution, and taking the reconstructed sound pressure distribution and the reconstructed sound velocity distribution as initial values to cycle to the wave equation solving in step one;

[0009] Step three, multi-grid iteration reconstruction: performing the joint reconstruction iteration of steps one to two in a coarse grid until the coarse grid iteration convergence condition is reached, interpolating the sound velocity distribution reconstructed in the coarse grid to a fine grid, and performing the joint reconstruction iteration of steps one to two in the fine grid scale by taking the interpolated sound velocity distribution as the initial value until the fine grid iteration convergence condition is reached, and finally taking the output of the fine grid as the reconstructed sound pressure distribution and the reconstructed sound velocity distribution.

[0010] Further, the joint reconstruction least square problem is specifically:

[0011] ;

[0012] wherein, represents the sound pressure distribution, represents the sound velocity distribution, represents the real detection signal received by the detector, represents a system matrix related to , i.e. represents the simulated detection signal , represents the square of the two-norm, and represent the regularization terms of and , and are regularization coefficients, is a sampling matrix, and the sound field obtained by calculation is interpolated to the detector position to obtain the sound field signal collected by the detector.

[0013] Further, the joint reconstruction least square problem is specifically: the first sub-problem is specifically:

[0014] ;

[0015] Question 2 is specifically:

[0016] ;

[0017] in, and They are and The sound pressure function and the sound speed function.

[0018] Furthermore, in step 2, the sound velocity distribution is known, and the sound pressure distribution is solved using the conjugate gradient method, which is:

[0019] ;

[0020] in, When the sound velocity distribution is known, the optimal sound pressure distribution estimate is obtained by solving the regularized least squares problem using the conjugate gradient method. is the Laplace operator, is the real detection signal, is the sound pressure distribution, is the regularization parameter, is the system matrix, is the transpose of the system matrix.

[0021] Furthermore, the The effect of the operator is obtained through the one-stage adjoint method, specifically:

[0022] Set up the one-stage state equation, add the one-stage mass source term to the one-stage adjoint state equation, and use the solved one-stage adjoint field as , thus obtaining Operator pair Calculation of

[0023] The quality source item of the first stage for:

[0024] ;

[0025] in, represents the sound pressure distribution, is the spatial position, is the detector position, is the Dirac function, is the total time for the sound wave to propagate, is the positive time, To reverse time, is the time-reversed residual signal.

[0026] Furthermore, in step 2, the second-stage accompanying field is specifically:

[0027] A two-stage adjoint state equation is set up, and the difference between the simulated detection signal and the real detection signal is time-reversed and taken as a negative value. This is added as a two-stage mass source term to the one-stage adjoint state equation, and the two-stage adjoint field is obtained by solving it.

[0028] Furthermore, the second-stage mass source term Specifically:

[0029] ;

[0030] in, is the spatial position, is the detector position, is the real detection signal, To simulate the detection signal, is the Dirac function, is the total time for the sound wave to propagate, is the positive time, To reverse time.

[0031] Furthermore, in step 2, the sound velocity gradient is calculated using the two-stage adjoint field, specifically:

[0032] Calculate the sound pressure distribution versus time After the first-order difference of The dimension is integrated and multiplied by the coefficients of the sound speed distribution generated in the previous iteration to obtain the gradient of the sound speed function with respect to the sound speed, that is, the sound speed gradient.

[0033] Furthermore, the sound velocity gradient is specifically:

[0034] ;

[0035] ;

[0036] in, is a function of the speed of sound Speed ​​of sound The sound velocity gradient, is the medium density, is the real detection signal, is the sound speed distribution generated in the previous iteration, The speed of sound The relevant system matrix, To simulate the detection signal, For location Time The sound pressure distribution, is the time reversal result of the two-stage adjoint field, is the regularization coefficient, The speed of sound Regularization term.

[0037] Furthermore, in the multi-grid iterative reconstruction, the root mean square error of the sound pressure distribution of two adjacent reconstructions is used as a convergence condition, and the iteration is stopped when the root mean square error is less than a set error value.

[0038] The advantages of the multi-grid iterative sound pressure-sound velocity synchronous reconstruction method in photoacoustic tomography provided by the present invention are: by constructing a joint reconstruction least squares problem and dividing it into two problems for solution, no additional measurement is required, and the relationship between the detection signal and the sound pressure distribution and sound velocity distribution is established through a mathematical model. Through an alternating cycle iteration strategy, high-quality sound pressure distribution and sound velocity distribution are reconstructed in a single measured photoacoustic signal, thereby improving image quality and accelerating the joint reconstruction speed; through multiple grids for joint reconstruction, it is verified in a two-dimensional simulation that the initial sound pressure distribution and sound velocity distribution can be simultaneously reconstructed using only a single detected photoacoustic signal without any prior conditions. Compared with the traditional algorithm that assumes that the tissue sound velocity is consistent, this reconstruction method improves image quality. Compared with the method combined with ultrasonic tomography, no prior sound velocity knowledge is required, which reduces equipment cost. The application of the multi-grid strategy also accelerates the speed of the joint reconstruction algorithm and improves reconstruction accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 It is a schematic diagram of the process of the present invention;

[0040] Figure 2 are the initial conditions for the simulation settings, (a) is the schematic diagram of the simulation settings of the detector, (b) is the schematic diagram of the initial sound pressure distribution, and (b) is the schematic diagram of the initial sound velocity distribution;

[0041] Figure 3 Schematic diagram of the iterative results of the joint reconstruction of the two problems using only a fine grid. (a) is the schematic diagram of the reconstructed sound pressure distribution, (b) is the difference between the reconstructed sound pressure distribution and the true sound pressure distribution; (c) is the sound pressure profile intensity diagram;

[0042] Figure 4 Schematic diagram of the iterative results of the joint reconstruction of the two problems only on the fine grid. (d) is the schematic diagram of the reconstructed sound speed distribution, (e) is the difference between the reconstructed sound speed distribution and the true sound speed distribution, and (f) is the sound speed profile intensity diagram;

[0043] Figure 5 Schematic diagram of the interpolation of the iterative results of the joint reconstruction of the two problems on a coarse grid onto a fine grid. (a) is a schematic diagram of the reconstructed sound velocity distribution, (b) is a schematic diagram of the difference between the reconstructed sound velocity distribution and the true sound velocity distribution, and (c) is an intensity diagram at the sound velocity distribution profile.

[0044] Figure 6 For Figure 5 Based on the results of the joint reconstruction iteration on the coarse grid, the joint reconstruction iteration is continued on the fine grid. Schematic diagram of the results, (a) is a schematic diagram of the reconstructed sound pressure distribution, (b) is a schematic diagram of the difference between the reconstructed sound pressure distribution and the true sound pressure distribution, and (c) is a sound pressure profile intensity diagram;

[0045] Figure 7 For Figure 5 Based on the results of the joint reconstruction cycle iteration on the coarse grid, the joint reconstruction cycle iteration results are continued on the fine grid. (d) is the reconstructed sound speed distribution, (e) is the difference diagram between the reconstructed sound speed distribution and the true sound speed distribution, and (f) is the sound speed profile intensity diagram. DETAILED DESCRIPTION

[0046] The technical solutions of the present invention are described in detail below through specific embodiments. Numerous specific details are set forth in the following description to facilitate a full understanding of the present invention. However, the present invention can be implemented in many other ways than those described herein, and those skilled in the art may make similar modifications without departing from the scope of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below.

[0047] like Figures 1 to 7 As shown, the method for synchronous reconstruction of sound pressure and sound velocity in photoacoustic tomography proposed by the present invention includes:

[0048] Step 1: construct a joint reconstruction least squares problem by minimizing the two norms of the detector's real detection signal and the simulated detection signal, and divide the problem into two problems to solve. The simulated detection signal is a signal obtained by solving the wave equation using the sound pressure distribution and the sound velocity distribution as initial values;

[0049] Step 2: Problem 1 is to solve the sound pressure distribution given the known sound velocity distribution. Problem 2 is to use the sound pressure distribution obtained in Problem 1 as the known sound pressure distribution to construct the second-stage mass source term in the two-stage adjoint field calculation, and then obtain the sound velocity distribution. The reconstructed sound pressure distribution and sound velocity distribution are used as initial values ​​to loop back to solve the wave equation in Step 1.

[0050] Step 3, multi-grid iterative reconstruction: perform the joint reconstruction iteration of step 1 to step 2 on the coarse grid until the coarse grid iterative convergence condition is reached, interpolate the sound velocity distribution reconstructed by the coarse grid to the fine grid, and use the interpolated sound velocity distribution as the initial value at the fine grid scale to perform the joint reconstruction iteration of step 1 to step 2 until the fine grid iterative convergence condition is reached. Finally, the output of the fine grid is used as the reconstructed sound pressure distribution and reconstructed sound velocity distribution.

[0051] This embodiment addresses the main drawbacks of current joint reconstruction methods, namely, their high pathological nature, the tendency for the iterative process to fall into local optimal solutions, the slowness of reconstruction due to the iterative nature of the joint reconstruction, and the need for prior knowledge of the sound velocity distribution for initial segmentation and estimation. A reconstruction method is provided that constructs a joint reconstruction least squares problem and solves it using these two problems, eliminating the need for additional measurements. A mathematical model is used to establish the relationship between the detection signal and the sound pressure and sound velocity distributions. Using an alternating cyclic iteration strategy, high-quality sound pressure and sound velocity distributions are reconstructed from a single measured photoacoustic signal, improving image quality and accelerating joint reconstruction.

[0052] That is, by performing joint reconstruction at different scales of grids, the initial sound pressure distribution and sound velocity distribution can be simultaneously reconstructed from the photoacoustic signal alone, improving image quality and accelerating the speed of joint reconstruction. By performing joint reconstruction using multiple grids, it has been verified in a two-dimensional simulation that the initial sound pressure distribution and sound velocity distribution can be simultaneously reconstructed from the photoacoustic signal alone, without any prior conditions. Compared to traditional algorithms that assume that tissue sound velocity is consistent, this embodiment improves image quality. Compared to methods combined with ultrasonic tomography, this embodiment does not require prior knowledge of sound velocity, reducing equipment costs. The application of a multi-grid strategy also accelerates the speed of the joint reconstruction algorithm and improves reconstruction accuracy.

[0053] It should be noted that in the process of solving the two problems in step 2, problem 1 is optimized through a small loop iteration using the conjugate gradient method to obtain the reconstructed sound pressure distribution, which is used as the known sound pressure distribution in solving problem 2; problem 2 is optimized through a small loop iteration using the L-BFGS method to obtain the reconstructed sound velocity distribution, and the reconstructed sound pressure distribution and sound velocity distribution are again used as the input of the wave equation in step 1. Steps 1 and 2 are used as steps to be executed in each iteration round of the large loop iteration. This embodiment uses a multi-grid reconstruction method with alternating coarse and fine grids to perform large loop iterations, and the reconstructed sound pressure distribution and sound velocity distribution output from the previous iteration round are used as the input of step 1 in the current iteration round.

[0054] In one embodiment, step 1 is to construct a joint reconstruction least squares problem by minimizing the two norms of the detector's real detection signal and the simulated detection signal and to solve the problem by dividing it into two problems, specifically:

[0055] The propagation model of sound waves in lossless media can be modeled by the following first-order wave equation:

[0056] , (1);

[0057] , (2);

[0058] , (3);

[0059] Its initial value conditions are:

[0060] , (4);

[0061] in, Represents spatial location, Indicates time, is the particle vibration velocity, It's location At the moment The sound pressure field, is the medium density, is the sound density distribution, is the speed of sound, is the pressure distribution at the initial moment, For the organization in position The speed of sound wave propagation at is the sound pressure gradient.

[0062] By numerically solving the acoustic wave equations of formulas (1) to (4) and obtaining the sound pressure distribution at the detector position, the simulated detection signal can be obtained. , the whole process of solving the acoustic wave equation can be written as ,in represents the initial sound pressure distribution, that is, here , Represents the system matrix, which is used to describe the forward propagation of sound waves in the tissue. is the sampling matrix, and the calculated sound field is interpolated to the detector position to obtain the sound field signal collected by the detector.

[0063] In this embodiment, the joint reconstruction problem can be described as a least squares problem, which is solved by minimizing the bi-norm value between the real detection signal and the simulated signal generated by the estimated pressure distribution and sound velocity. The mathematical formula is as follows:

[0064] , (5);

[0065] in, represents the sound pressure distribution, represents the sound velocity distribution, represents the real photoacoustic signal received by the detector, Represents The relevant system matrix, This represents the simulated detection signal , represents the square of the two-norm, and Representatives and The regularization term is used to alleviate the ill-posed nature of the problem. and The regularization coefficient of is the sampling matrix, and the calculated sound field is interpolated to the detector position to obtain the sound field signal collected by the detector.

[0066] The meaning of formula (5) is to find the solution that minimizes the right side of the equation. and , temporarily ignore the sampling matrix in the subsequent problem analysis , mainly explore and .

[0067] In one embodiment, in step 2, problem 1 is to solve the sound pressure distribution using the conjugate gradient method, given a known sound velocity distribution, as follows:

[0068] To solve the above problem, an alternating iterative optimization method is used to divide the joint reconstruction problem into two stages. The first stage is to find the sound pressure distribution given the sound velocity distribution:

[0069] , (6);

[0070] in, for The sound pressure function is .

[0071] The solution to problem one is achieved through the conjugate gradient method. During the solution process of problem one, problem one is solved through an internal small loop (implemented by the conjugate gradient method), thereby obtaining the sound pressure distribution after multi-grid iterative reconstruction corresponding to the current iteration round.

[0072] The sound speed distribution at this time As we know, the problem now becomes a linear least squares problem. Considering the smoothness of the initial sound pressure distribution, the regularization term Calculated using Laplace canonical The second-order differential of . For the solution of formula (6), considering the large system matrix of the forward process, it is difficult to directly invert it, so a gradient-based method is used to solve it. right The gradient of can be directly derived from (6) as:

[0073] , (7);

[0074] in, for right The gradient, is the system matrix.

[0075] Since the system matrix Large, find the system matrix It is difficult in itself, and this embodiment considers As an operator, rather than an explicit system matrix, then The solution of the operator requires the use of the adjoint state method, specifically for the forward operator The acoustic wave equations (1) to (3) and the initial value conditions (4), the first-stage adjoint state equation can be written as:

[0076] , (8);

[0077] , (9);

[0078] , (10);

[0079] The initial conditions are:

[0080] , (11);

[0081] in,

[0082] , (12);

[0083] in, is the mass source term of one stage, It is a stage accompanying field. For the organization in position The speed of sound wave propagation at is the detector position, is the total time for the sound wave to propagate, is the positive time, To reverse time, is the Dirac function, which means that only when hour There is a value, that is, there is only a one-stage mass source term at the detector position. is the time-reversed residual signal, which is used as the parameter of the mass source term. The residual signal is selected as the time-reversed residual signal, that is, the time-reversed difference between the real detection signal of the detector and the simulated detection signal. The first-stage adjoint field is obtained by numerically solving formulas (8) to (12) That is .

[0084] There are many gradient-based methods. This example uses the conjugate gradient (CG) method to iteratively solve Problem 1. The CG method is memory-efficient and highly efficient when solving linear least squares optimization problems for large sparse matrices. However, it is only applicable to positive definite systems of equations. This example rewrites the problem to be solved to facilitate the use of the CG method:

[0085] , (13);

[0086] in, When the sound velocity distribution is known, the optimal sound pressure distribution estimate is obtained by solving the regularized least squares problem using the conjugate gradient method. is the Laplace operator, is the real detection signal, is the sound pressure distribution, in formula (13) is the regularization parameter, is the system matrix, is the transpose of the system matrix.

[0087] The system matrix at this time can be regarded as , the conjugate gradient method is used to solve formula (13). This method does not need to calculate the specific gradient in formula (7), which involves The calculation of is performed using the one-stage adjoint state equation described in the above formulas (8) to (12). The conjugate gradient method is an existing method and will not be described in detail in this embodiment.

[0088] In step 2 of this embodiment, problem 2 uses the sound pressure distribution obtained in problem 1 as a known sound pressure distribution to construct a second-stage mass source term in the two-stage adjoint field calculation, thereby obtaining the sound velocity distribution, specifically:

[0089] Given the sound pressure distribution, the sound velocity distribution is calculated as follows:

[0090] , (14)

[0091] in, for The sound speed function is .

[0092] The solution to Problem 2 is achieved through the L-BFGS method. During the solution process of Problem 2, Problem 2 is solved through an internal small loop (implemented using the L-BFGS method), thereby obtaining the sound velocity distribution after multi-grid iterative reconstruction corresponding to the current iteration round.

[0093] The sound pressure distribution at this time It is known that this problem is a nonlinear optimization problem. right The gradient of is not as simple as formula (7). This embodiment refers to the full-wave inversion problem in geophysics and also adopts the adjoint state method, that is, setting a two-stage adjoint state equation, and taking the negative value of the difference between the simulated detection signal and the real detection signal after time reversal, and adding it as the second-stage mass source term to the first-stage adjoint state equation to obtain the second-stage adjoint field. Based on the two-stage adjoint field, we can get right The sound velocity gradient is expressed as:

[0094] , (15);

[0095] in, is a function of the speed of sound Speed ​​of sound The sound velocity gradient, is the medium density, To solve the sound velocity distribution generated in the previous iteration in the inner small loop of Problem 2, is the time-reversed result of the two-stage adjoint field. At this time, the two-stage mass source term calculated by the two-stage adjoint field is and the first-stage quality source term Slightly different:

[0096] , (16);

[0097] in, The calculation result of is the time reversal of the result of numerically solving equations (8) to (11) and equation (16), and For location Department The sound pressure at the moment and the simulated detection signal different, What is needed is the amount of the entire sound pressure distribution, and Only the sound pressure changes at the detector location are stored.

[0098] Based on the above calculations, the sound pressure field can be obtained and the two-stage adjoint field time reversal result , and then the sound velocity gradient is calculated by formula (15): .

[0099] For the solution of problem 2, this embodiment adopts the L-BFGS method (quasi-Newton method) for iterative solution. BFGS is the initials of the four founders of this method, and L represents limited memory. This method needs to calculate the value of right The gradient of , which has a better solution effect for nonlinear problems and is suitable for complex objective functions. Similarly, the LBFGS method is also an existing method, and this embodiment will not go into details.

[0100] In steps 1 and 2, the joint reconstruction performs alternating iterative solutions to the two-stage problem. Considering the strong pathological nature of problem 2, problem 1 is solved first during the iterative process. As one embodiment, the specific steps of the two-stage alternating iterative process are as follows:

[0101] (a1) Assume that the sound velocity distribution is that of background water and the initial sound pressure distribution is zero. When solving the problem, assume that the sound pressure distribution is solved iteratively in a homogeneous medium.

[0102] (a2) Solve Problem 2 based on the sound pressure distribution obtained in Problem 1 to obtain a new sound velocity distribution;

[0103] (a3) Using the sound velocity distribution obtained in (a2) and the sound pressure distribution obtained in (a1) as initial values, solve Problem 1 again. Repeat this cycle until the iteration convergence condition is met.

[0104] In step 3 of this embodiment, the multi-grid iterative reconstruction is specifically as follows:

[0105] During the two-stage alternating iteration process, errors with larger image frequencies are first rapidly attenuated, while errors with smaller frequencies are attenuated more slowly. The multi-grid iteration strategy considers iterating on a fine grid to attenuate errors with larger frequencies, then iterating on a coarse grid to rapidly attenuate errors with smaller frequencies. Finally, to ensure accuracy, the previously obtained result is used as the initial guess, and further iteration is performed on the fine grid. This fine-coarse-fine grid iteration is called a V-cycle. There are also W-cycles and full-grid cycles. The full-grid cycle begins with the coarsest grid, with V and W cycles in between. This strategy can find a good initial value.

[0106] Since the second problem in the joint reconstruction is highly pathological, a good initial value can greatly alleviate the pathological condition. The joint reconstruction based on the prior sound velocity segmentation provides a relatively good initial sound velocity. This embodiment considers the use of a multi-grid method, first performing joint reconstruction in a coarse grid (looping steps 1 to 2), and then performing joint reconstruction in a fine grid (looping steps 1 to 2). Since the joint reconstruction is relatively time-consuming, this embodiment only considers obtaining a good initial value to accelerate the joint reconstruction of the fine grid, that is, performing the joint reconstruction iteration of the two-stage problem on the coarse grid (looping steps 1 to 2). ), until the coarse grid iteration convergence condition is reached, the sound velocity distribution reconstructed by the coarse grid is interpolated to the fine grid, and the interpolated sound velocity distribution is used as the initial value at the fine grid scale to perform a joint reconstruction iteration of the two-stage problem (looping steps one to two) until the fine grid iteration convergence condition is reached, and finally the output of the fine grid is used as the reconstructed sound pressure distribution and the reconstructed sound velocity distribution. In the joint reconstruction loop iteration, the sound pressure distribution and the sound velocity distribution output by the previous iteration are used as the wave equation input of step one in the current iteration, and steps one to two are optimized through multi-grid iterative reconstruction.

[0107] It is understandable that this embodiment does not exclude the possibility of implementing joint reconstruction iterations of the two-stage problem by adding structures such as a V loop and a W loop.

[0108] As an embodiment;

[0109] Based on Matlab programming, version number is R2019b, the computer main configuration is: Intel Xeon Gold 6226RCPU, 128G memory, NVIDIA TITAN RTX graphics card. Figure 1 As shown, the joint reconstruction method in this embodiment includes numerical solution of the wave equation, CG solution of the sound pressure distribution of problem one, solution of the new pressure field and adjoint field based on the acoustic wave equation and the adjoint equation, solution of the sound velocity gradient of the sound velocity function with respect to the sound velocity in problem two, and LBFGS solution of the sound velocity distribution. The specific operations are as follows:

[0110] (b1) Numerical solution of wave equation;

[0111] Formulas (1) to (4) are calculated using the numerical simulation software Matlab. The k-Wave toolbox is used to solve the acoustic wave equation. The input sound pressure distribution for the first solution is zero, and the sound speed distribution is uniformly distributed with a size of 1480 m / s. The k-Wave toolbox is an open source acoustic simulation tool.

[0112] (b2) Use the conjugate gradient method (CG) to solve the sound pressure distribution:

[0113] The optimized sound pressure distribution of problem 1 is obtained by solving formula (13) using the conjugate gradient method, which involves The calculation of the operator is directly solved by using the formulas (1) to (4), involving The operator is calculated by one-stage adjoint field, i.e. solved by the formulas (8) to (12).

[0114] (b3) Solving a new pressure field and two-stage adjoint field based on the wave equation and two-stage adjoint state equation;

[0115] Taking the sound pressure distribution solved in (b2) as a new initial sound pressure distribution estimate, the sound pressure distribution is calculated according to the method in (b1) And the simulated detection signal The two-stage mass source term in the two-stage adjoint field calculation is formula (16), i.e. the simulated detection signal is subtracted from the real detection signal , and after time reversal, the negative value is taken, and the two-stage adjoint field is calculated.

[0116] (b4) Calculating the sound velocity function of problem two The sound velocity gradient of the sound velocity c:

[0117] The first-order difference of the sound pressure distribution with respect to time is calculated, multiplied by the time reversal result of the two-stage adjoint field, and integrated in the time dimension, multiplied by the coefficient formed by the sound velocity distribution generated in the previous iteration in the internal small loop of solving problem two, to obtain the gradient of the sound velocity function with respect to the sound velocity, i.e. the sound velocity gradient.

[0118] Specifically, the sound pressure distribution is calculated according to formula (15), and the first-order difference with respect to time is calculated to approximate the gradient calculation, which is multiplied by the time reversal result of the two-stage adjoint field, and integrated in the time dimension, and in the discrete case, it is summed, and then multiplied by the coefficient , where is the sound velocity distribution generated in the previous iteration in the internal small loop of solving problem two, and finally the sound velocity gradient of with respect to the sound velocity is calculated.

[0119] (b5) Solving the sound velocity distribution of problem two by LBFGS;

[0120] LBFGS is a quasi-Newton method, which has good solving effect on nonlinear problems, and its iteration format is

[0121] , (17);

[0122] where is the sound velocity distribution generated in the previous iteration in the internal small loop of solving problem two. The reconstructed sound velocity distribution at the second iteration is used as the input or final output of the next iteration in the inner small loop to solve problem 2. To solve the inner small loop of problem 2 The sound velocity distribution to be reconstructed at the iteration, To solve the problem 2, the sound speed function in the inner small loop is The sound velocity gradient, gradient Through the two-stage adjoint field calculation, is the search step size, which represents the amplitude of each iterative update in the inner small loop for solving problem 2. It requires a one-dimensional line search to calculate. is an approximate Hessian matrix, which represents second-order partial derivatives. Since the true Hessian matrix is ​​difficult to calculate, an approximate Hessian matrix is ​​necessary, which is the purpose of the quasi-Newton method. Since the fundamental purpose of joint reconstruction is to optimize the initial sound pressure distribution, the accuracy of the sound velocity is less demanding. Therefore, the convergence conditions can be appropriately relaxed to speed up the calculation.

[0123] (b6) Loop iteration:

[0124] The new pressure distribution obtained by (b5) and sound speed distribution As the new initial value (i.e., the new initial sound pressure distribution and initial sound velocity distribution), repeat the above steps (b1) to (b5) to perform large loop iteration. When the upper limit of the large loop iteration setting is reached or convergence is achieved, the large loop iteration is terminated and the new sound pressure is output. and sound speed distribution Since the essential purpose of joint reconstruction is to optimize the initial sound pressure distribution map, the overall convergence condition is set to stop iteration when the root mean square error (RMSE) of the two reconstructed sound pressures is less than 0.01. At this time, the multi-grid iterative reconstruction sound pressure has almost no change, and convergence is achieved. The RMSE is calculated as:

[0125] , (18);

[0126] in, is the sound pressure distribution obtained in the current iteration round, is the sound pressure distribution obtained in the last iteration round, is the number of grids. Since RMSE is used as the convergence condition in both coarse and fine grids, that is, in the coarse grid iteration, is the number of coarse grids, in the fine grid iteration, is the number of fine grids.

[0127] In order to verify the advantages of this embodiment for simultaneous reconstruction of sound pressure and sound velocity, the coarse grid size used in the multi-grid iteration in this embodiment is 128×128, with a grid spacing of 0.5 mm, and the fine grid size is 256×256, with a grid spacing of 0.25 mm. Figure 2 The simulation settings shown further illustrate the technical solution of this embodiment. Figure 2 As shown in (a), the outer black dots (outer circle) are the detector distribution, the dark red (branch-shaped) in the middle is the initial sound pressure distribution, which represents the blood vessel absorption density distribution, and the light red area (center circle) is the heterogeneous medium area, representing the high sound speed area. The setting value is 1530 m / s, and the background sound speed is 1480 m / s (the sound speed can be set as needed). Figure 2 (b) shows the sound pressure distribution, Figure 2 (c) shows the sound velocity distribution. Figure 2 The red lines in (b) and (c) represent the corresponding profile positions of the subsequent profile intensity maps. The actual detection signal was simulated using the Matlab k-Wave toolbox. The image grid size was set to 256×256, the grid spacing was 0.25 mm, and the sampling frequency was 40 MHz. Fifteen joint reconstruction iterations were performed on this fine grid, and the results are shown in Figure 2. Figure 3 and 4 As shown, through Figure 3 (b) and Figure 3 As can be seen from (c), when the iterative result of the sound pressure reaches the set convergence condition under the fine grid, the maximum difference between the reconstructed sound pressure distribution and the real sound pressure distribution is about 0.08. The sound pressure profile shows that the peak value is not very consistent. Figure 4 (e) and Figure 4 As can be seen from (f), the maximum difference between the reconstructed sound speed distribution and the true sound speed distribution is about 20 m / s. The sound speed profile shows that the reconstructed sound speed result has sound speed variation disturbance outside the actual high sound speed area.

[0128] The real detection signal is downsampled to a sampling frequency of 20 MHz to obtain a real signal corresponding to a coarse grid. The coarse grid size is set to 128×128 and the grid spacing is 0.5 mm. 10 iterations are performed on the coarse grid, and then the reconstructed sound velocity result is interpolated to a 256×256 grid. The interpolated result is as follows: Figure 5 As shown, through Figure 5 (b) and Figure 5As can be seen from (c), when the iterative result of the sound pressure reaches the set convergence condition under the coarse grid, the maximum difference between the reconstructed sound speed distribution and the actual sound speed distribution is about 30 m / s. However, the reconstructed sound speed result divides the actual high sound speed area well, and there is almost no sound speed change disturbance outside the actual high sound speed area. The sound speed distribution is used as the initial value to perform joint reconstruction on the fine grid, and 5 cycles of iteration are performed to obtain the following Figure 6 and 7 The multi-grid reconstruction results are shown. Figure 6 (b) and Figure 6 As can be seen from (c), when the iterative result of the sound pressure reaches the set convergence condition under the multi-grid method, the maximum difference between the reconstructed sound pressure distribution and the real sound pressure distribution is about 0.03, and the sound pressure profile shows a good match at the peak. Figure 7 (e) and Figure 7 (f) shows that the maximum difference between the reconstructed sound speed distribution and the actual sound speed distribution is only 10 m / s, and the reconstructed sound speed result can well divide the actual high sound speed area.

[0129] By comparison Figures 3 to 7 Compared with iterating directly on the fine grid, the strategy of first coarse grid and then fine grid in this embodiment significantly improves the image quality, the difference between the sound pressure distribution and the sound velocity distribution is smaller, the cross-section intensity overlap is higher, and the time required for 15 cycles of iteration directly on the fine grid is 3 hours, while this embodiment takes about 40 minutes to perform 10 cycles of iteration on the coarse grid and about 1 hour and 10 minutes to perform 5 cycles of iteration on the fine grid. Therefore, the overall time required for the multi-grid strategy is about 1 hour and 50 minutes, which saves about 1 / 3 of the time compared to directly calculating on the fine grid.

[0130] Therefore, this embodiment applies a multi-grid iterative reconstruction strategy to joint reconstruction, accelerating reconstruction speed and improving reconstruction accuracy. Furthermore, this embodiment does not require additional measurements. By establishing a mathematical model that relates the detection signal to the sound pressure and speed distributions, and using an alternating, iterative strategy, high-quality sound pressure and speed distributions can be reconstructed from a single photoacoustic signal.

[0131] The above description is only a preferred specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any technician familiar with the technical field, within the technical scope disclosed by the present invention, who makes equivalent replacements or changes based on the technical solution and inventive concept of the present invention, should be covered by the scope of protection of the present invention.

Claims

1. A method for simultaneous reconstruction of acoustic pressure and acoustic velocity in photoacoustic tomography, characterized in that: include: Step 1: construct a joint reconstruction least squares problem by minimizing the two norms of the detector's real detection signal and the simulated detection signal, and divide the problem into two problems to solve. The simulated detection signal is a signal obtained by solving the wave equation using the sound pressure distribution and the sound velocity distribution as initial values; Step 2: Problem 1 is to solve the sound pressure distribution given the known sound velocity distribution. Problem 2 is to use the sound pressure distribution obtained in Problem 1 as the known sound pressure distribution to construct the second-stage mass source term in the two-stage adjoint field calculation, and then obtain the sound velocity distribution. The reconstructed sound pressure distribution and sound velocity distribution are used as initial values ​​to loop back to solve the wave equation in Step 1. Step 3, multi-grid iterative reconstruction: perform the joint reconstruction iteration of steps 1 and 2 on the coarse grid until the coarse grid iterative convergence condition is reached, interpolate the sound velocity distribution reconstructed by the coarse grid to the fine grid, and use the interpolated sound velocity distribution as the initial value at the fine grid scale to perform the joint reconstruction iteration of steps 1 and 2 until the fine grid iterative convergence condition is reached. Finally, the output of the fine grid is used as the reconstructed sound pressure distribution and reconstructed sound velocity distribution; The joint reconstruction least squares problem is specifically: Problem 1 is specifically: ； Question 2 is specifically: ； in, represents the sound pressure distribution, represents the sound velocity distribution, Indicates the actual detection signal received by the detector, Represents The relevant system matrix, This represents the simulated detection signal , represents the square of the two-norm, and Representatives and The regularization term, and express and The regularization coefficient of and They are and The sound pressure function and the sound speed function.

2. The reconstruction method according to claim 1, characterized in that The joint reconstruction least squares problem is specifically: ； in, is the sampling matrix, and the calculated sound field is interpolated to the detector position to obtain the sound field signal collected by the detector.

3. The reconstruction method according to claim 1, wherein: In step 2, the sound velocity distribution is known, and the sound pressure distribution is solved using the conjugate gradient method. The conjugate gradient method solves the problem as follows: ； in, When the sound velocity distribution is known, the optimal sound pressure distribution estimate is obtained by solving the regularized least squares problem using the conjugate gradient method. is the Laplace operator, is the real detection signal, is the sound pressure distribution, is the regularization parameter, is the system matrix, is the transpose of the system matrix.

4. The reconstruction method according to claim 3, characterized in that: described The effect of the operator is obtained through the one-stage adjoint method, specifically: Set up the one-stage state equation, add the one-stage mass source term to the one-stage adjoint state equation, and use the solved one-stage adjoint field as , thus obtaining Operator pair Calculation of The quality source item of the first stage for: ； in, represents the sound pressure distribution, is the spatial position, is the detector position, is the Dirac function, is the total time for the sound wave to propagate, is the positive time, To reverse time, is the time-reversed residual signal.

5. The reconstruction method according to claim 1, characterized in that: In step 2, the two-stage accompanying field is specifically: A two-stage adjoint state equation is set up, and the difference between the simulated detection signal and the real detection signal is time-reversed and taken as a negative value. This is added as a two-stage mass source term to the one-stage adjoint state equation, and the two-stage adjoint field is obtained by solving it.

6. The reconstruction method according to claim 5, characterized in that: The second stage mass source term Specifically: ； in, is the spatial position, is the detector position, is the real detection signal, To simulate the detection signal, is the Dirac function, is the total time for the sound wave to propagate, is the positive time, To reverse time.

7. The reconstruction method according to claim 1, characterized in that: In step 2, the sound velocity gradient is calculated using the two-stage adjoint field, specifically: Calculate the sound pressure distribution versus time After the first-order difference of The dimension is integrated and multiplied by the coefficients of the sound speed distribution generated in the previous iteration to obtain the gradient of the sound speed function with respect to the sound speed, that is, the sound speed gradient.

8. The reconstruction method according to claim 7, characterized in that: The sound velocity gradient is specifically: ； ； in, is a function of the speed of sound Speed ​​of sound The sound velocity gradient, is the medium density, is the real detection signal, is the sound speed distribution generated in the previous iteration, The speed of sound The relevant system matrix, To simulate the detection signal, For location Time The sound pressure distribution, is the time reversal result of the two-stage adjoint field, is the regularization coefficient, The speed of sound Regularization term.

9. The reconstruction method according to claim 1, characterized in that: In the multi-grid iterative reconstruction, the root mean square error of the sound pressure distribution of two adjacent reconstructions is used as a convergence condition, and the iteration is stopped when the root mean square error is less than the set error value.

Images