Acoustic transmission model-based algebraic iteration photoacoustic image imaging method

Abstract

The invention discloses an algebraic iteration photoacoustic image imaging method based on an acoustic transmission model, and belongs to the technical field of photoacoustic / ultrasonic imaging. Comprising the following steps: S1, designing a photoacoustic imaging system based on photoacoustic and ultrasonic bimodal; S2, discretizing an ultrasonic signal forward propagation equation; s3, calculating a forward model; and S4, constructing a forward propagation equation and iteratively solving the initial sound pressure. A traditional photoacoustic reconstruction method is influenced by limited bandwidth and limited visual angle of a transducer, so that reconstructed images are often influenced by artifacts. According to the method, the iterative reconstruction method based on the sound transmission model is introduced, so that the phenomena of high artifacts and multiple side lobes existing in a traditional image reconstruction algorithm are effectively inhibited. The invention provides a new method for reconstructing the photoacoustic image, and the resolution of the reconstructed image based on the photoacoustic effect is remarkably improved. In the process of iteratively solving the forward propagation equation, the regularization method is introduced, the convergence direction of the iterative solution is controlled by the L2 norm information of the initial sound pressure, and the robustness of the reconstructed image influenced by noise is effectively enhanced.

Description

Technical Field

[0001] This invention belongs to the field of photoacoustic / ultrasound imaging technology, specifically relating to a photoacoustic image imaging method based on an algebraic iteration of an acoustic transmission model. Background Technology

[0002] Photoacoustic imaging, as an emerging method for imaging biological functional information, combines the high contrast of optical imaging with the high penetration depth of ultrasound imaging, attracting widespread attention in the field of biological tissue therapy. However, due to the inherent limited viewing angle and bandwidth of transducers in photoacoustic imaging, its resolution is affected by negative value artifacts, fringe artifacts, and splitting artifacts. This can lead to inaccurate imaging of the target area during the imaging process, contaminating the reconstructed biological modal information. Fortunately, the rapid development of model-based photoacoustic image reconstruction technology, especially its applications in image reconstruction and signal processing, has provided new solutions to this problem. Summary of the Invention

[0003] To overcome the limitations of transducer-limited viewing angle and bandwidth in existing photoacoustic imaging methods, this invention provides an algebraic iterative photoacoustic imaging method based on an acoustic transmission model. This method starts with the initial sound pressure from light excitation and establishes a forward propagation model based on the partial differential equation of sound wave propagation in a medium. This model then constructs a sound pressure propagation matrix equation and iteratively solves the matrix equation using an improved conjugate gradient method and regularization techniques. This reconstructs a high-quality photoacoustic image that largely overcomes the limitations of transducer-limited viewing angle and bandwidth, and exhibits robustness against certain noise effects, making it suitable for clinical use in biological photoacoustic imaging.

[0004] The technical solution adopted in this invention is:

[0005] The photoacoustic image imaging method based on the algebraic iteration of the acoustic transmission model includes the following steps:

[0006] S1. Design a photoacoustic imaging system based on photoacoustic and ultrasonic dual-modal imaging;

[0007] S2. Discretization of the forward propagation equation of ultrasonic signals;

[0008] S3. Calculate the forward model;

[0009] S4. Construct the forward propagation equation and iteratively solve for the initial sound pressure;

[0010] S5. System Verification and Application.

[0011] Compared with the prior art, the present invention has the following advantages:

[0012] 1. Traditional photoacoustic reconstruction methods are often affected by artifacts due to the limited bandwidth and viewing angle of the transducer. This invention introduces an iterative reconstruction method based on the acoustic transmission model, which effectively suppresses the high artifacts and multi-sidelobe phenomena present in traditional image reconstruction algorithms.

[0013] 2. This invention proposes a novel method for photoacoustic image reconstruction, which significantly improves the resolution of reconstructed images based on photoacoustic effects.

[0014] 3. In the process of iteratively solving the forward propagation equation, this invention introduces a regularization method to control the convergence direction of the iterative solution with the L2 norm information of the initial sound pressure, which effectively enhances the robustness of the reconstructed image to noise. Attached Figure Description

[0015] Figure 1 This is a flowchart of the invention;

[0016] Figure 2 This is a comparison of photoacoustic images reconstructed by the iterative reconstruction method and the FBP method. Detailed Implementation

[0017] To better understand the purpose, structure, and function of this invention, the invention will be described in further detail below with reference to the accompanying drawings.

[0018] like Figure 1 As shown, this invention provides an algebraic iterative photoacoustic imaging method based on an acoustic transmission model, comprising the following steps:

[0019] Step 1: Hardware composition and configuration of the photoacoustic imaging system based on photoacoustic and ultrasonic dual-modality:

[0020] This invention provides a photoacoustic imaging system based on photoacoustic and ultrasound dual-modal imaging, aiming to provide medical personnel with basic functional information needed within biological tissues. The photoacoustic imaging system based on photoacoustic and ultrasound dual-modal imaging consists of the following main components: a pulsed laser, a beam shaping system, an arbitrary-shaped ultrasound transducer, a multi-channel data acquisition card, and a central computer. Specifically: the pulsed laser emits nanosecond-width pulsed light to excite light absorbers within the target area, generating photoacoustic signals; the beam shaping system adjusts the shape of the laser beam, including the beam size and divergence angle, to ensure the laser beam can be transmitted to the optical system with sufficient coupling efficiency; the ultrasound transducer acquires ultrasound and photoacoustic dual-modal signals; the multi-channel data acquisition card digitizes the ultrasound and photoacoustic dual-modal signals received by the ultrasound transducer and transmits the processed data to the central computer for further analysis and processing; the central computer provides a human-computer interface for interaction between the user and the system. Users can issue commands such as start imaging, end imaging, and save data on the central computer. The central computer processes user commands and automatically matches the working sequence of each system component to ensure the normal operation of the system. During operation, the user controls the pulsed laser by issuing commands through the central computer. The central computer automatically matches the timing of the ultrasonic transducer and the multi-channel data acquisition unit under different operating modes to complete the image reconstruction. Ultrasonic and photoacoustic signals are essentially sound pressure signals, and the ultrasonic transducer can acquire information in both ultrasonic and photoacoustic modes. Through preset timing design, this system can achieve real-time ultrasonic / photoacoustic dual-modal image reconstruction of the target area, thereby providing rich information on biological tissue function.

[0021] Step 2: Discretization of the forward propagation equation of the ultrasonic signal:

[0022] Based on the partial differential linear equations describing sound wave propagation, the integral equation satisfied by the ultrasonic signal propagating in the medium can be derived as follows:

[0023]

[0024] in, For the first The initial sound pressure signal received by each transducer; For Grunison parameters; Sensitivity factor; Speed ​​of sound; This is a coefficient related to sound attenuation; It is a directional function; It is a constant related to the direction angle; For located The initial sound pressure level at that location; The pulse response of the transducer; coordinate point Place and No. Transducer position The square of the Euclidean distance between them; Indicates time information.

[0025] The above equation can be reduced to a surface integral of the first kind within a spherical cap with the product of time and sound speed as the radius and the coordinates of the ultrasonic transducer as the center. Since computers cannot process continuous integrals, it is necessary to convert the continuous integral expression into a discrete form:

[0026] The continuous sound propagation space is discretized into a Cartesian grid, and the integral of the first type of curve in the above formula is transformed into a discrete form. By dividing the continuous integral arc segment into a finite number of subdivided arc segments, and starting the loop with each sound source point as the sound pressure excitation point, the time series of each transducer receiving the sound source point is calculated.

[0027] Includes the following steps:

[0028] Step 2.1: Division of the region of interest and the Cartesian grid:

[0029] The coordinates of each spatial point in the optical / ultrasound imaging region are represented using Cartesian coordinates. During the imaging phase, it is not necessary to reconstruct all elements of the imaging region. The coordinates of the center point of the imaging region are specified as the origin (0, 0). A regular quadrilateral is selected as the region of interest, and the side length of the quadrilateral is specified as needed. The coordinates of the region of interest are discretized into a finite number of Cartesian grids. For example, if the side length of the quadrilateral is L, the side length of each Cartesian grid is... Then you can use and Two grid matrices representing the x and y coordinates represent the coordinates of each Cartesian grid point, where Each column is:

[0030]

[0031] in, Represents all coordinates of the grid points in the x-axis direction, and similarly, Each action:

[0032]

[0033] Simply provide the matrix index of the corresponding discrete grid point, and you can access it. and The coordinates of the corresponding matrix elements are determined, and at the same time, based on the relative position of the ultrasonic transducer to the origin of the Cartesian coordinate system, a corresponding Cartesian coordinate is assigned to each ultrasonic transducer unit. Complete step 2.1.

[0034] Step 2.2, coarse division of the integral arc segment:

[0035] The surface integral of the first kind degenerates into the line integral of the first kind in a two-dimensional plane. Since the region of interest is finite, the integration arc segments will intersect with the boundary of the region of interest, requiring coarse partitioning of the integration arc segments. The partitioning method is as follows: First, using the product of the current time point and the speed of sound as the arc segment radius, solve for the intersection points of the arc segment and the region boundary according to the angular direction, and select the arc segment within the region, denoted as . ,in, It is represented by Cartesian coordinates of its starting and ending points.

[0036] Step 2.3, Refining the integral arc segment:

[0037] After completing the coarse division of the integral arc segment, for each coarse arc segment... Further subdivision of the arc segments is performed by selecting each coarse arc segment. Find the intersection points of the line with each Cartesian grid line, and save each intersection point to... middle, It is a 4-row matrix. Each column of the matrix represents an intersection point. The first row stores the global x-coordinate of the intersection point, the second row stores the global y-coordinate of the intersection point, the third row stores the polar coordinate angle of the intersection point relative to the center of the circle, and the fourth row indicates whether these points are intersections of perpendicular grid points. If they are, they are stored as 1, and if they are not, they are stored as 0. Every two intersection points determine a subdivision arc segment, and then the corresponding integration can be performed on each subdivision arc segment.

[0038] Step 3, Calculation of the forward model:

[0039] The calculation of each subdivided arc point corresponds to a part of a column of the forward model matrix. Specifically, the first... The integral value of the i-th transducer at time point t corresponds to the i-th column idp in the model matrix. arrive The elements of the row, where idp represents the integer pixel index of the reconstructed image pixel in column-major order. This represents the total time received by each transducer. Since the integral formula is continuous, and the previous steps discretized the integration arc into a finite number of sub-segments, the integral on each sub-segment remains continuous. Therefore, a bilinear interpolation method is needed to solve for the integral formula on each sub-segment, as follows:

[0040]

[0041] in, For the first The sound pressure signal values ​​generated by each transducer at the corresponding 4 pixel coordinate points at the current calculation time t; It is a constant factor that is related to the relative position of the Cartesian grid endpoints where the subdivided arc segments are located and the ultrasonic transducer. To subdivide the arc segment between its two endpoints The difference in coordinates, i.e. ,same, To subdivide the arc segment between its two endpoints Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the coordinates of the transducer element are... Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the coordinates of the transducer element are... Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The coordinates between the lower left grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The coordinates between the lower left grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The angle subtended by the two endpoints of the subdivided arc segment with the center of the arc segment as the center.

[0042] Step 4: Constructing the forward propagation equation and iteratively solving for the initial sound pressure:

[0043] Previously, in step 3, the integral calculation of the photoacoustic signal generated by the coordinates of each image pixel at each time t was completed. Step 4 requires storing the corresponding calculated sound pressure values ​​in the matrix according to the matrix arrangement order, constructing the forward propagation equation, and after processing the original sound pressure signal obtained by the transducer in actual work, solving the initial sound pressure using the algebraic iteration method.

[0044] Step 4.1, Construction of forward propagation equations:

[0045] In step 3, two cycles are required: transducer unit The loop for the current time point t and the loop for each transducer unit, in each calculation, simultaneously saves the pixel index value idp, the time point t, and the transducer number. Calculated sound pressure value ,correspond The values ​​are filled into the index coordinates of the forward propagation matrix A as follows:

[0046]

[0047] in, This represents the total number of transducer units; The total sampling time for each transducer; idp is the index value of the pixel calculated according to the column priority rule; t is the current time. The calculation of each element of the forward propagation matrix can only be completed after both the inner and outer loops have finished.

[0048] However, only the forward propagation matrix elements have been calculated so far; it still needs to be converted into matrix form. This will involve calculating the time series from each sound source. according to The indexes are arranged in a column from top to bottom. Following the order of sound source points, the time indices in each column are then arranged row-by-row from left to right to form the forward propagation matrix. .

[0049] Thus, we construct a system of linear equations in matrix form for forward propagation:

[0050]

[0051] Where p is the time series of the signals received by each transducer, arranged by transducer index along the column direction. Let be the forward propagation matrix, and x be the sound pressure sequence in column order of arranging the initial sound pressure matrix into column vectors.

[0052] Step 4.2, Transducer raw data preprocessing:

[0053] The raw data received by the transducer is often affected by white noise. Noise outside the transducer's bandwidth is filtered out by applying a bandpass filter with a center frequency equal to the transducer's center frequency in the frequency domain. Signal processing techniques that enhance the signal-to-noise ratio, such as Wiener filters, can also be used in this step.

[0054] Step 4.3, Iterative solution of initial sound pressure:

[0055] Due to the matrix obtained in step 4.1 Since the matrix is ​​a large, sparse, non-square matrix, directly solving the forward propagation equation is inefficient. Therefore, an improved conjugate gradient method is used to iteratively solve the matrix. Simultaneously, to increase the reconstructed image's resistance to noise, a regularization technique is employed, incorporating the L2 norm information of the initial sound pressure level as part of the gradient calculation in the iterative solution process. Specifically, this involves solving the equation...

[0056]

[0057] The problem of finding a solution is transformed into the following optimization problem:

[0058]

[0059] in, This is the final estimated initial sound pressure level; The initial sound pressure level used in the calculation; For regularization parameters; This is the square of the L2 norm. The photoacoustic reconstructed image based on the acoustic transmission model can only be obtained through iterative optimization.

[0060] Step 5, System Verification and Application:

[0061] To verify the effectiveness of this invention, relevant experiments were conducted to verify its reliability and applicability in practical applications. Experimental results show that the algebraic iterative photoacoustic image reconstruction method based on the acoustic transmission model performs excellently in photoacoustic imaging tasks. In in vivo experiments, the system was successfully applied to photoacoustic image reconstruction of target areas in mouse biological tissues after local injection of photothermal nanoprobes. Using this reconstruction method, the system demonstrated excellent performance in local imaging resolution and imaging signal-to-noise ratio, verifying its reliability and broad applicability in practical applications.

[0062] Figure 2 This paper presents a comparison between the photoacoustic reconstructed image of a finger based on an algebraic iterative reconstruction using a sound transmission model and the photoacoustic temperature image reconstructed using traditional FBP. Figure 2 It can be clearly seen that the photoacoustic reconstruction imaging achieved by the algorithm of this invention can image smaller targets and produce less reconstruction artifacts. Compared with the image reconstructed by traditional FBP, this invention can remove artifacts more effectively and improve imaging resolution.

[0063] This invention demonstrates high-precision photoacoustic imaging capabilities, providing a more advanced non-invasive method for detecting functional information in biological tissues. The invention integrates model building and iterative reconstruction modules. By constructing a forward propagation model based on the physical essence of the acoustic wave propagation equation and introducing regularization parameters to participate in the solution, it significantly improves imaging resolution compared to the high artifacts and multiple sidelobes in FBP-reconstructed images, showing significant advantages and great application potential in photoacoustic imaging. Furthermore, the use of multimodal image fusion methods ensures that multimodal information of the target area can be acquired during treatment, providing staff with more intuitive and visualized feedback on the state of biological tissues and improving treatment outcomes.

[0064] It is understood that the present invention has been described through some embodiments, and those skilled in the art will recognize that various changes or equivalent substitutions can be made to these features and embodiments without departing from the spirit and scope of the invention. Furthermore, under the teachings of the present invention, these features and embodiments can be modified to adapt to specific situations and materials without departing from the spirit and scope of the invention. Therefore, the present invention is not limited to the specific embodiments disclosed herein, and all embodiments falling within the scope of the claims of this application are within the protection scope of the present invention.

Claims

1. A photoacoustic image imaging method based on algebraic iteration of an acoustic transmission model, characterized in that: Includes the following steps: S1. Design a photoacoustic imaging system based on photoacoustic and ultrasonic dual-modal imaging; S2. Discretization of the forward propagation equation of ultrasonic signals; S3. Calculate the forward model; S4. Construct the forward propagation equation and iteratively solve for the initial sound pressure; S5. System Verification and Application.

2. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 1, characterized in that: In S1, the photoacoustic imaging system based on photoacoustic and ultrasonic dual-modality includes A pulsed laser is used to emit pulses of light with nanosecond-level pulse widths to excite light absorbers in a target area and generate photoacoustic signals. An optical path shaping device is used to adjust the spot shape of the pulsed laser beam to ensure that the spot can be transmitted to the optical system with sufficient coupling efficiency. An ultrasonic transducer is used to acquire ultrasonic and photoacoustic dual-mode signals. The multi-channel data acquisition card is used to convert the ultrasonic and photoacoustic dual-mode signals received by the ultrasonic transducer into digital-to-analog signals and transmit the processed data to the central computer for further analysis and processing. The central computer provides the interface between the user and the system, waits for the user to issue image commands, and then matches the working sequence of each system component to ensure that the system works normally.

3. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 1, characterized in that: In the discretization of the ultrasonic signal forward propagation equation in S2, When an ultrasonic signal propagates in a medium, it satisfies the following integral equation: in, For the first The initial sound pressure signal received by each transducer; For Grunison parameters; Sensitivity factor; Speed ​​of sound; This is a coefficient related to sound attenuation; It is a directional function; It is a constant related to the direction angle; For located The initial sound pressure level at that location, The pulse response of the transducer. coordinate point Place and No. Transducer position The square of the Euclidean distance between them; Indicates time information.

4. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 3, characterized in that: The discretization of the forward propagation equation of the ultrasonic signal in S2 includes the following steps: S21. Divide the region of interest into a Cartesian grid; S22. Perform coarse division of the integral arc segment; S23. Perform fine division of the integral arc segment.

5. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 4, characterized in that: The specific process of dividing the region of interest and the Cartesian grid in S21 is as follows: the coordinates of each spatial point in the optical / ultrasonic imaging region are represented by Cartesian coordinates, the coordinates of the center point of the imaging region are specified as the origin (0, 0), a regular quadrilateral is selected as the region of interest, and the coordinates of the region of interest are discretized into a finite number of Cartesian grids.

6. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 5, characterized in that: In the forward model calculation in S3, the bilinear interpolation method is used to solve for the integral formula on each subdivided arc segment, as follows: in, For the first The sound pressure signal values ​​generated by each transducer at the corresponding 4 pixel coordinate points at the current calculation time t; It is a constant factor that is related to the relative position of the Cartesian grid endpoints where the subdivided arc segments are located and the ultrasonic transducer. To subdivide the arc segment between its two endpoints The difference in coordinates, i.e. , To subdivide the arc segment between its two endpoints Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the coordinates of the transducer element are... Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the coordinates of the transducer element are... Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The coordinates between the top right grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The coordinates between the lower left grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The coordinates between the lower left grid of the subdivided arc segment and the midpoint between the two endpoints of the subdivided arc segment. Difference in coordinates; The angle subtended by the two endpoints of the subdivided arc segment with the center of the arc segment as the center.

7. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 6, characterized in that: The construction of the forward propagation equation and the iterative solution of the initial sound pressure in S4 include the following steps: S41. Construct the forward propagation equation; S42. Preprocess the raw transducer data; S43. Iteratively solve for the initial sound pressure level.

8. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 7, characterized in that: In the construction of the forward propagation equation in S41, Construct a system of linear equations in matrix form for forward propagation: in, This is the time series of signals received by each transducer, arranged by transducer index along the column direction. Forward propagation matrix, This is to arrange the initial sound pressure matrix into a sequence of column vectors based on column priority.

9. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 8, characterized in that: In S42, a bandpass filter with a center frequency equal to the transducer's center frequency is applied in the frequency domain to filter out noise outside the transducer's bandwidth, or a Wiener filter is used to enhance the signal-to-noise ratio, thus preprocessing the original transducer data.

10. The photoacoustic image imaging method based on algebraic iteration of the acoustic transmission model according to claim 8, characterized in that: In step S43, the initial sound pressure is iteratively solved by using an improved conjugate gradient method to iteratively solve the matrix. At the same time, regularization technology is used to incorporate the L2 norm information of the initial sound pressure as part of the gradient calculation into the iterative solution process.

Images