CN109091109B - Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator - Google Patents

Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator Download PDF

Info

Publication number
CN109091109B
CN109091109B CN201810706816.4A CN201810706816A CN109091109B CN 109091109 B CN109091109 B CN 109091109B CN 201810706816 A CN201810706816 A CN 201810706816A CN 109091109 B CN109091109 B CN 109091109B
Authority
CN
China
Prior art keywords
matrix
filtering
information
imaging
photoacoustic
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810706816.4A
Other languages
Chinese (zh)
Other versions
CN109091109A (en
Inventor
陶超
芮纬
刘晓峻
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University
Original Assignee
Nanjing University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University filed Critical Nanjing University
Priority to CN201810706816.4A priority Critical patent/CN109091109B/en
Publication of CN109091109A publication Critical patent/CN109091109A/en
Application granted granted Critical
Publication of CN109091109B publication Critical patent/CN109091109B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/0093—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy
    • A61B5/0095—Detecting, measuring or recording by applying one single type of energy and measuring its conversion into another type of energy by applying light and detecting acoustic waves, i.e. photoacoustic measurements
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7225—Details of analog processing, e.g. isolation amplifier, gain or sensitivity adjustment, filtering, baseline or drift compensation

Landscapes

  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Animal Behavior & Ethology (AREA)
  • General Health & Medical Sciences (AREA)
  • Biophysics (AREA)
  • Biomedical Technology (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Medical Informatics (AREA)
  • Molecular Biology (AREA)
  • Surgery (AREA)
  • Signal Processing (AREA)
  • Pathology (AREA)
  • Public Health (AREA)
  • Veterinary Medicine (AREA)
  • Power Engineering (AREA)
  • Acoustics & Sound (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physiology (AREA)
  • Psychiatry (AREA)
  • Ultra Sonic Daignosis Equipment (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses an image reconstruction method of photoacoustic tomography based on full-matrix filtering and time reversal operators, which comprises the steps of firstly, constructing an information matrix according to photoacoustic signals received by an ultrasonic transducer array, wherein the information matrix comprises information containing direct waves and information containing scattered waves causing interference; then, according to the characteristic inverse diagonal property of the direct wave part, through operations such as matrix 45-degree rotation, matrix filtering, matrix 45-degree reverse rotation reduction and the like, scattered wave components in the information matrix can be filtered; finally, a spatial distribution image of the intensity of the acoustic wave emitted by the optical absorber can be obtained from the filtered information matrix by using a time reversal operator. The method provided by the invention can effectively filter the interference of the scattered wave and completely reserve effective detection information, thereby effectively improving the reconstruction accuracy of the photoacoustic image in and behind the scattering medium, and having higher usability and wide applicability.

Description

Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator
Technical Field
The invention relates to an optimized photoacoustic tomography image reconstruction method based on full-matrix filtering and time reversal operators, and belongs to the photoacoustic tomography technology.
Background
Photoacoustic tomography has been rapidly developed in recent years as a non-invasive biomedical imaging technique. Photoacoustic tomography has the advantages of high resolution of deep tissues by an acoustic method and the advantages of an optical method in functional imaging, molecular imaging, imaging contrast and the like [ 1 ].
However, photoacoustic tomography generally applies to media with uniform acoustic characteristics, and imaging of scattering media remains a big problem faced by photoacoustic imaging technology. In the scattering medium, the received photoacoustic signal includes a scattered wave component after multiple scattering by the scattering layer in addition to a direct wave from the imaging target. Due to irregular distribution of scatterers, scattered wave components are random, image information is difficult to extract, and scattered waves can bring speckle noise and image distortion to a photoacoustic image. When the scattered wave component is dominant, strong scattered wave interference will seriously reduce the depth and imaging quality of photoacoustic imaging, even make it impossible to obtain a recognizable image.
To solve this problem, various photoacoustic imaging schemes suitable for scattering media have been proposed. For example, a method of using correlation matrix filtering is proposed to improve the photoacoustic imaging quality [ 2 ] in a scattering medium, but the method has the limitations of narrow imaging area, low resolution and the like because the detected photoacoustic data cannot be fully utilized.
Reference to the literature
【1】L.V.Wang and S.Hu,“Photoacoustic tomography:in vivo imaging from organelles to organs,”Science 335(6075),1458–1462(2012)
【2】Wei Rui,Chao Tao,Xiaojun Liu,“Photoacoustic imaging in scattering media by combining a correlation matrix filter with a time reversal operator,”Optics Express 25(19),22840.
Disclosure of Invention
The purpose of the invention is as follows: in order to improve the image quality of photoacoustic tomography in the acoustic scattering medium, the invention provides an optimized photoacoustic tomography image reconstruction method based on matrix filtering and time reversal operators, so as to effectively improve the image quality of photoacoustic imaging in the acoustic scattering medium. Specifically, according to the inverse diagonal characteristics of an information matrix corresponding to a direct sound wave, interference signals of a scattered wave are filtered through operations of 45-degree anticlockwise rotation, zero padding, matrix filtering, 45-degree clockwise rotation reduction of the matrix and the like, and then time reversal imaging is carried out, so that the interference of the scattered sound wave is reduced, and the image quality of photoacoustic imaging in a scattering medium is improved.
The technical scheme is as follows: in order to achieve the purpose, the invention adopts the technical scheme that:
the image reconstruction method of the photoacoustic tomography based on the full-matrix filtering and the time reversal operator comprises the following steps:
step 1: irradiating an imaged sample by using pulse laser, and after a light absorber in the sample absorbs the energy of the pulse laser, radiating ultrasonic waves to surrounding media due to a photoacoustic effect;
step 2: receiving an ultrasonic signal radiated by a sample by using an ultrasonic transducer array, wherein the ultrasonic transducer array is provided with N transducer units which are linearly arranged along the x-axis direction, the center distance between every two adjacent transducer units is w, and the transverse coordinate of the nth transducer unit is as follows:
Figure BDA0001715523810000021
and step 3: carrying out operations including time gating and short-time Fourier transform on the received ultrasonic signals to construct an information matrix K with the size of N multiplied by N;
and 4, step 4: sequentially performing operations including 45-degree anticlockwise rotation, zero filling, filtering and 45-degree clockwise rotation reduction on the information matrix K, filtering scattered wave components and reserving direct wave components;
and 5: singular value decomposition is carried out on the filtered information matrix, and a time reversal operator is applied to the eigenvectors corresponding to the largest M singular values to obtain the photoacoustic image with the corresponding depth;
step 6: obtaining correction parameters according to an acoustic transfer function in a uniform medium, and reducing image defects including image distortion and false contrast caused by limited visual angle ultrasonic signal detection;
and 7: and (3) changing the time gating area, repeating the steps 3-6, obtaining the photoacoustic images at different depths, and finally obtaining the photoacoustic image of the whole area.
Specifically, in step 3, the construction process of the information matrix K is as follows:
step 31: intercept data received by an ultrasonic transducer at
Figure BDA0001715523810000022
The signal in the time window is Fourier transformed to obtain the frequency domain expression P of the nth transducer unitn(T, f) the frequency domain represents Pn(T, f) including a direct wave signal PD(T, f) and scattered wave signal PS(T, f) two moieties; wherein, T represents the shortest time of acoustic wave propagation from the imaging plane to the transducer array, also called time gating; Δ t represents the width of the time window, f represents the frequency;
the direct wave is a signal which starts from an imaging area, is directly transmitted to a transducer without being refracted and scattered by any scatterer and is detected; the scattered wave refers to a signal which is transmitted to the transducer and detected after being refracted and scattered once or more times from the imaging area through the scatterer;
the information vector of the ultrasonic transducer array is denoted as P (T, f) ═ P1(T,f),…,Pn(T,f),…,PN(T,f)]The direct wave signal is recorded as
Figure BDA0001715523810000031
Scattered wave signal PS(T, f) is represented by
Figure BDA0001715523810000032
Direct wave signal of nth transducer unit
Figure BDA0001715523810000033
Satisfies the following conditions:
Figure BDA0001715523810000034
wherein A is0＝Z1/2K is 2 pi f/c, and c is medium sound velocity; z ═ cT denotes the vertical distance of the imaging plane to the transducer array, also referred to as the depth of the corresponding time gate T; x represents the X-axis coordinate value of an arbitrary point on the imaging plane, the arbitrary point being represented as (X, Z);
if L scattering paths exist between any point (X, Z) on the imaging plane and the nth transducer unit, the scattered wave signal of the nth transducer unit
Figure BDA0001715523810000035
Satisfies the following conditions:
Figure BDA0001715523810000036
wherein A islAnd slThe amplitude and phase parameters corresponding to the first scattering path; since the scatterers are randomly distributed, AlAnd slCorresponding randomness is also embodied;
step 32: and (3) constructing an information matrix K by using the information vector P:
Figure BDA0001715523810000037
the information matrix K is divided into two parts according to whether random parameters are contained: the first entry on the right of the equal sign is recorded as matrix KCThe remaining entries on the right side of the equal sign are recorded as a matrix KR；
Matrix KCDoes not contain any random parameter, namely is independent of the distribution position of the scatterer, and the matrix elements have the following correlation:
Figure BDA0001715523810000038
wherein m is 1,2, …, N, 1 is not less than m plus or minus q is not less than N;
Figure BDA0001715523810000039
representation matrix KCThe phase difference between the elements located anti-diagonals and at a distance q from the diagonals,
Figure BDA0001715523810000041
representation matrix KCThe element at the (m-q, m + q) position,
Figure BDA0001715523810000042
representation matrix KCAt the (m, m) position;
matrix KRAll elements in the matrix contain random parameters and show randomness, so that the matrix elements have no correlation.
Specifically, in the step 4, the operation process of filtering the scattered wave component in the information matrix K is as follows:
step (ii) of41: unlike the method in document [ 1 ], the invention rotates the information matrix K by 45 ° counterclockwise according to the following formula and fills the corresponding elements with zeros to ensure that the matrix dimensions are unchanged, resulting in two corresponding matrices a1 ═ a1u,v}、A2＝{a2u,v}：
If N is an odd number:
Figure BDA0001715523810000043
Figure BDA0001715523810000044
if N is an even number:
Figure BDA0001715523810000045
Figure BDA0001715523810000046
wherein k ism,nRepresenting the element of the information matrix K at the (m, n) position,
Figure BDA0001715523810000047
the dimension of the matrix A1 is N, and the dimension of the matrix A2 is N-1; matrix A1 and matrix A2 are collectively denoted by matrix A and by NARepresents the dimension of matrix a;
step 42: constructing a filter matrix
Figure BDA0001715523810000051
Filtering the matrix A, wherein the expression of a filtering vector C is as follows:
Figure BDA0001715523810000052
wherein the content of the first and second substances,
Figure BDA0001715523810000053
is a conjugate transpose of C; u-1, 2, …, NA，
Figure BDA0001715523810000054
Step 43: filtering the matrix A, and recording the filtered matrix A as the matrix AF：
AF＝FA＝F(AC+AR)＝FAC+FAR
Wherein, the matrix ACThe part of the matrix A without random parameters, the matrix ARThe matrix A is a part containing random parameters; first item FAC＝ACMatrix ACBefore and after filtering, and matrix ARAfter the filtering is used, the scattered wave is weakened powerfully, and the scattered wave part is effectively filtered;
step 44: will matrix AFClockwise rotating 45 degrees to obtain a filtered information matrix K, and recording the filtered information matrix K as a matrix KF：
If (m-n) is an odd number,
Figure BDA0001715523810000055
if (m-n) is an even number,
Figure BDA0001715523810000056
wherein the content of the first and second substances,
Figure BDA0001715523810000057
representation matrix KFOf the element at the (m, n) position,
Figure BDA0001715523810000058
representation matrix A1FAn element at the ((m-N)/2+ (N +3)/4, (m + N)/2) position,
Figure BDA0001715523810000059
representation matrix A2FIs located in ((m-N-1)/2+ (N +3)/4, (m + N-1)) /2) element at position, matrix A1FIs a filtered matrix A1, matrix A2FIs a filtered matrix a 2;
matrix KFThe size of (2) is still NxN, which is consistent with the information matrix K, and no information dimension is lost before and after filtering.
Specifically, in the step 5, the time reversal imaging process according to the filtered information matrix is as follows:
step 51: for matrix KFSingular value decomposition:
Figure BDA00017155238100000510
wherein, UFIs a unitary matrix of order NxN, the row vector of which is the matrix KFLeft singular vectors of (d); lambdaFDetermining an NxN order diagonal matrix for the semi-positive matrix;
Figure BDA00017155238100000511
is a unitary matrix of order N x N,
Figure BDA00017155238100000512
is a VFBy conjugate transposition of (V)FThe row vector of is the matrix KFRight singular vectors of (d); lambdaFThe diagonal elements are the matrix KFThe singular value of (a);
step 52: calculating a direct wave transmission matrix G between an arbitrary point (X, Z) on an imaging plane and an ultrasonic transducer array, wherein an element G at the position of (m, n) in the matrix Gm,nComprises the following steps:
Figure BDA0001715523810000061
wherein r ism,nIs the distance between the nth transducer unit and the mth transducer unit in the ultrasonic transducer array;
step 53: the eigenvectors corresponding to the largest M singular values and the direct wave transmission matrix G apply a time reversal operator, the ith singular value
Figure BDA0001715523810000062
The photoacoustic signal intensity at depth Z is:
Figure BDA0001715523810000063
wherein the content of the first and second substances,
Figure BDA0001715523810000064
represents VFThe ith row of singular vectors; superposing the strength of the photoacoustic signals correspondingly obtained by the maximum M singular values to obtain a photoacoustic image at the depth Z
Figure BDA0001715523810000065
Figure BDA0001715523810000066
Is IFThe mth element of the vector, which represents the preliminary imaging value for the location on the imaging plane corresponding to the mth transducer element.
Specifically, in step 6, the generation and utilization processes of the correction parameters are as follows:
step 61: imaging correction value corresponding to mth transducer unit at depth Z
Figure BDA0001715523810000067
Comprises the following steps:
Figure BDA0001715523810000068
wherein the content of the first and second substances,
Figure BDA0001715523810000069
is gm,nConjugation of (1);
step 62: using correction values for imaging
Figure BDA00017155238100000610
To the preliminary imaging value
Figure BDA00017155238100000611
Correction is made so as to obtain a final imaging value
Figure BDA00017155238100000612
Figure BDA00017155238100000613
I.e. the final imaging result is filtered by the imaging result IFAnd a correction matrix IRThe corresponding elements are obtained by dividing.
Has the advantages that: the image reconstruction method of photoacoustic tomography based on matrix filtering and time reversal operator provided by the invention has the following advantages: 1. according to the invention, the information matrix received by the ultrasonic transducer array is utilized, so that the filtering processing can be effectively carried out on the sound wave signals, the interference of a scattering medium on the imaging effect is reduced, and high-quality target imaging is obtained; 2. the invention only needs one excitation and receiving process, has simple and efficient operation, higher practicability and wide applicability; 3. the invention can not cause information loss before and after filtering, and can utilize the received whole effective information; 4. according to the invention, the correction parameters are obtained according to the acoustic transfer function in the uniform medium, and the image defects such as image distortion, false contrast and the like caused by limited visual angle ultrasonic signal detection are reduced.
Drawings
FIG. 1 is a schematic diagram of an image reconstruction principle of optimized photoacoustic tomography based on matrix filtering and a time reversal operator;
FIG. 2 is a schematic diagram of matrix rotation, zero padding and filtering operations;
FIG. 3 is a comparison diagram of characteristic space matrices before and after filtering, where 3(a) is a real part image of a signal matrix K received without a scattering layer, 3(b) is a real part image of the signal matrix K after a scattering layer is added to a model, and 3(c) is a matrix K filtered by a filtering method in document [ 2 ]FThe real part image 3(d) is filtered by adopting the filtering method provided by the inventionMatrix K of wave rearFThe real part image of (a);
fig. 4 is a comparison graph of imaging results of a conventional filtering imaging method and an optimized filtering imaging method, where 4(a) is an imaging result of a conventional beamforming imaging method, 4(b) is an imaging result of a time reversal method without filtering, 4(c) is an imaging result of a filtering imaging method in document [ 2 ], and 4(d) is an imaging result of an optimized filtering imaging method proposed by the present invention.
Detailed Description
The present invention will be further described with reference to the accompanying drawings.
Fig. 1 is a schematic diagram of a principle of an image reconstruction method of photoacoustic tomography based on matrix filtering and time reversal operator, taking imaging of a small steel ball behind a scattering layer as an example, the method includes the following steps.
Step 1: 40 small steel balls with the diameter of 0.8mm are randomly and irregularly arranged to serve as a random scattering layer and are arranged between the ultrasonic transducer and the target area. Due to irregular distribution of the scattering bodies, the received photoacoustic signal includes not only a direct wave from the imaging target but also a scattered wave component that is scattered multiple times by the scattering layer. The scattering layer is used for generating interference signals. Four small steel balls with the diameter of 0.8mm are placed on the scattering layer and then are used as imaging targets. The imaged target is irradiated by a pulse laser, and after the small steel balls in the imaged target absorb laser energy, ultrasonic waves are radiated to the surrounding medium due to the photoacoustic effect (the system is shown in figure 1).
Step 2: an ultrasonic transducer array with the sampling frequency of 20MHz is used for receiving ultrasonic signals radiated by a sample, the ultrasonic transducer array is provided with N101 units, the distance between every two adjacent units is w 0.5mm, and therefore the transverse coordinate of the nth unit is xn＝[n-(N+1)/2]w, where N is an odd number, the method of the present invention is still true when N is an even number.
And step 3: and carrying out time gating and short-time Fourier transform operation on the received ultrasonic signals to obtain signal components corresponding to the frequency f, and constructing an N multiplied by N information matrix K.
Step 31: truncated ultrasound transducer receptionTo be in the time range
Figure BDA0001715523810000071
The signal in the spectrum is subjected to short-time Fourier transform to obtain frequency domain expression, which is marked as Pn(T, f), wherein the subscript n denotes the corresponding nth transducer element. The frequency domain signals obtained by the N transducers form an information vector P (T, f) ([ P ]1(T,f),…,Pn(T,f),…,PN(T,f)]. The information vector P (T, f) contains the direct wave signal PD(T, f) and scattered wave signal PS(T, f) two parts (as shown in FIG. 1). The expression of the corresponding element of the direct wave signal satisfies:
Figure BDA0001715523810000081
wherein A is0＝Z1/2And k is 2 pi f/c, and c is medium sound velocity. Accordingly, the scattered wave signal corresponding element expression satisfies:
Figure BDA0001715523810000082
wherein A islAnd slThe amplitude and phase parameters corresponding to the first scattering path; since the scatterers are randomly distributed, AlAnd slCorresponding randomness is also exhibited.
Step 32: and (3) constructing an information matrix K by using the information vector P:
Figure BDA0001715523810000083
the matrix K is divided into two parts depending on whether random parameters are included. The first term on the right of the equal sign is denoted as KCThis term does not contain any random parameter, i.e. is independent of the location of the scatterer distribution, and the matrix elements have the following correlation:
Figure BDA0001715523810000084
in contrast, the remainder on the right of the equal sign contains random parameters, which are collectively denoted as KR。KRShows randomness without K between elementsCThe correlation embodied.
And 4, step 4: and carrying out operations such as 45-degree anticlockwise rotation, zero padding, filtering, 45-degree clockwise rotation reduction and the like on the information matrix K, filtering scattered wave components and reserving direct wave components.
Step 41: the information matrix is rotated 45 ° counter-clockwise, resulting in the corresponding matrices a1, a2, according to:
Figure BDA0001715523810000091
Figure BDA0001715523810000092
in the formula, m is u + v- (N +1)/2, and N is v-u + (N + 1)/2. The matrix A1 has a dimension of N and the matrix A2 has a dimension of N-1. For convenience of description, hereinafter, A represents A1 or A2 and N representsARepresenting the dimensions of matrix a.
The rotation and zero padding process is shown in fig. 2, and the original signal matrix is rotated by 45 ° counterclockwise. The coherence of the original matrix is shifted from the anti-diagonal direction to the column direction. The rotated matrix is then divided into two parts according to the parity of the number of columns, the light green circular parts corresponding to the odd columns and the grey diamond parts corresponding to the even columns. Finally, the two parts are separated and then each zero-filled (light blue square part) to obtain the matrices a1, a 2.
Step 42: constructing a filter matrix
Figure BDA0001715523810000093
The a matrix is filtered. The expression of the filter vector C is:
Figure BDA0001715523810000094
wherein u is 1,2, …, NA，
Figure BDA0001715523810000095
Step 43: filtering the matrix A to obtain a filtering matrix AF
AF＝FA＝F(AC+AR)＝FAC+FAR
First term FA in the above equationC＝ACCorrelation term ACRemains unchanged before and after filtering, while the random term ARAfter being acted on by the filter matrix, the scattered wave part is weakened strongly and effectively filtered out.
Step 44: will matrix AFClockwise rotating 45 degrees to obtain a filtered information matrix KF：
If (m-n) is an odd number,
Figure BDA0001715523810000096
if (m-n) is an even number,
Figure BDA0001715523810000097
information matrix K obtained after filteringFThe size of the filter is still NxN, the filter is consistent with the original signal matrix K, and no information dimension is lost before and after filtering.
When the selection time T is 58 mus (corresponding to the time at [ -5mm,87 mm)]Small steel ball), when the center frequency f is 2.0MHz, the comparison effect before and after matrix filtering is shown in fig. 3. Fig. 3(a) is a real part image of the received signal matrix K in the case of no scattering layer. At the moment, the matrix signal only contains direct wave components, and the image shows clear correlation. When we add scattering layer in the model, the scattered wave component in the received signal is dominant, and the real part image of the signal matrix K is shown in fig. 3(b), and the image shows strong randomness. FIG. 3(c) shows a matrix K after filtering by the filtering method of document 1FThe scattering wave in the central region of the image is disturbedTo a certain extent of filtering, a correlation pattern can be observed, however, a large amount of information of the surrounding area is missing. FIG. 3(d) is a diagram of a matrix K after filtering by the filtering method of the present inventionFThe real part image of (1). It can be observed that on the one hand the interference of the scattered wave components is more effectively filtered out and a clear correlation pattern can be observed, the result being highly close to 3 (a). On the other hand, the dimension of the matrix is kept unchanged, and no information loss occurs, so that better guarantee is provided for the subsequent imaging process.
And 5: and applying a time reversal method to the filtered matrix to obtain the imaging of the corresponding propagation interval.
Step 51: for filter matrix KFSingular value decomposition:
Figure BDA0001715523810000101
step 52: calculating a direct wave transmission matrix G between a point on a focal plane and an ultrasonic transducer array unit, wherein the element expression of the direct wave transmission matrix G is as follows:
Figure BDA0001715523810000102
in the formula rm,nIs the distance between the nth element of the transducer array and the mth point on the focal plane, m, N is (1,2, …, N).
Step 53: applying a time reversal operator to the eigenvector corresponding to the first several largest singular values and the transmission matrix G, the ith singular value
Figure BDA0001715523810000103
The photoacoustic signal intensity at depth Z is:
Figure BDA0001715523810000104
superposing the strength of the photoacoustic signals correspondingly calculated by the first M nonzero singular values to obtain a photoacoustic image with corresponding depth
Figure BDA0001715523810000105
Step 6: and obtaining correction parameters according to the acoustic transfer function in the uniform medium, and reducing image defects such as image distortion, false contrast and the like caused by limited visual angle ultrasonic signal detection.
Step 61: at a specified depth Z, the imaging correction value corresponding to the mth point on the focal plane is obtained by:
Figure BDA0001715523810000106
step 62: using the correction value
Figure BDA0001715523810000111
To the preliminary imaging value
Figure BDA0001715523810000112
Correction is made so as to obtain a final imaging value
Figure BDA0001715523810000113
The final imaging result is filtered by the imaging result IFAnd a correction matrix IRThe corresponding elements are divided to obtain:
Figure BDA0001715523810000114
and 7: and (3) changing the time gating area, repeating the steps 3-6, obtaining the photoacoustic images at different depths, and finally obtaining the photoacoustic image of the whole area.
Fig. 4 is a comparison graph of imaging results of a plurality of imaging methods, and fig. 4(a) and 4(b) are the imaging results of a conventional beamforming imaging method and a time reversal method without filtering, respectively. It can be seen that the imaging quality of the two methods is poor, a large amount of interference images and distortion exist in the background area, and the influence of scattered wave components on the imaging quality is obvious. Fig. 4(c) shows the imaging result of the filtering imaging method in document 1, in which the contrast between the light absorber and the background region is significantly improved, and the effect of the filtering process on the improvement of the imaging quality is reflected. However, this method only images the central region, and the two light absorbers in the upper and lower side regions are not imaged. Fig. 4(d) shows the imaging result of the optimized filtering imaging method proposed by the present invention. Because the utilization rate of effective information is improved, the contrast ratio of the light absorber and the background area is further improved, the imaging result of the complete area is obtained, and all the light absorbers are effectively imaged.
The above description is only of the preferred embodiments of the present invention, and it should be noted that: it will be apparent to those skilled in the art that various modifications and adaptations can be made without departing from the principles of the invention and these are intended to be within the scope of the invention.

Claims (5)

1. An image reconstruction method of photoacoustic tomography based on full matrix filtering and time reversal operators is characterized in that: the method comprises the following steps:
step 1: irradiating an imaged sample by using pulse laser, and after a light absorber in the sample absorbs the energy of the pulse laser, radiating ultrasonic waves to surrounding media due to a photoacoustic effect;
step 2: receiving an ultrasonic signal radiated by a sample by using an ultrasonic transducer array, wherein the ultrasonic transducer array is provided with N transducer units which are linearly arranged along the x-axis direction, the center distance between every two adjacent transducer units is w, and the transverse coordinate of the nth transducer unit is as follows:
Figure FDA0002926648550000011
and step 3: carrying out operations including time gating and short-time Fourier transform on the received ultrasonic signals to construct an information matrix K with the size of N multiplied by N;
and 4, step 4: sequentially performing operations including 45-degree anticlockwise rotation, zero filling, filtering and 45-degree clockwise rotation reduction on the information matrix K, filtering scattered wave components and reserving direct wave components;
and 5: singular value decomposition is carried out on the filtered information matrix, and a time reversal operator is applied to the eigenvectors corresponding to the largest M singular values to obtain the photoacoustic image with the corresponding depth;
step 6: obtaining correction parameters according to an acoustic transfer function in a uniform medium, and reducing image defects including image distortion and false contrast caused by limited visual angle ultrasonic signal detection;
and 7: and (3) changing the time gating area, repeating the steps 3-6, obtaining the photoacoustic images at different depths, and finally obtaining the photoacoustic image of the whole area.
2. The image reconstruction method of photoacoustic tomography based on full matrix filtering and time reversal operator as claimed in claim 1, characterized in that: in step 3, the construction process of the information matrix K is as follows:
step 31: intercept data received by an ultrasonic transducer at
Figure FDA0002926648550000012
The signal in the time window is Fourier transformed to obtain the frequency domain expression P of the nth transducer unitn(T, f) the frequency domain represents Pn(T, f) including a direct wave signal PD(T, f) and scattered wave signal PS(T, f) two moieties; wherein, T represents the shortest time of acoustic wave propagation from the imaging plane to the transducer array, also called time gating; Δ t represents the width of the time window, f represents the frequency;
the information vector of the ultrasonic transducer array is denoted as P (T, f) ═ P1(T,f),…,Pn(T,f),…,PN(T,f)]The direct wave signal is recorded as
Figure FDA0002926648550000021
Scattered wave signal PS(T, f) is represented by
Figure FDA0002926648550000022
Direct wave signal of nth transducer unit
Figure FDA0002926648550000023
Satisfies the following conditions:
Figure FDA0002926648550000024
wherein A is0＝Z1/2K is 2 pi f/c, and c is medium sound velocity; z ═ cT denotes the vertical distance of the imaging plane to the transducer array, also referred to as the depth of the corresponding time gate T; x represents the X-axis coordinate value of an arbitrary point on the imaging plane, the arbitrary point being represented as (X, Z);
if L scattering paths exist between any point (X, Z) on the imaging plane and the nth transducer unit, the scattered wave signal of the nth transducer unit
Figure FDA0002926648550000025
Satisfies the following conditions:
Figure FDA0002926648550000026
wherein A islAnd slThe amplitude and phase parameters corresponding to the first scattering path;
step 32: and (3) constructing an information matrix K by using the information vector P:
Figure FDA0002926648550000027
the information matrix K is divided into two parts according to whether random parameters are contained: the first entry on the right of the equal sign is recorded as matrix KCThe remaining entries on the right side of the equal sign are recorded as a matrix KR；
Matrix KCNot containing any random parameter, i.e. separation from scatterersThe positions of the cloth are independent, and the matrix elements have the following correlation:
Figure FDA0002926648550000028
wherein m is 1,2, …, N, 1 is not less than m plus or minus q is not less than N;
Figure FDA0002926648550000029
representation matrix KCThe phase difference between the elements located anti-diagonals and at a distance q from the diagonals,
Figure FDA00029266485500000210
representation matrix KCThe element at the (m-q, m + q) position,
Figure FDA00029266485500000211
representation matrix KCAt the (m, m) position;
matrix KRAll elements in the matrix contain random parameters and show randomness, so that the matrix elements have no correlation.
3. The image reconstruction method of photoacoustic tomography based on full matrix filtering and time reversal operator as claimed in claim 1, characterized in that: in step 4, the operation process of filtering the scattered wave component in the information matrix K is as follows:
step 41: rotating the information matrix K by 45 ° counterclockwise according to the following formula, and zero-filling the corresponding elements to ensure that the matrix dimension is not changed, resulting in two corresponding matrices a1 ═ a1u,v}、A2＝{a2u,v}：
If N is an odd number:
Figure FDA0002926648550000031
Figure FDA0002926648550000032
if N is an even number:
Figure FDA0002926648550000033
Figure FDA0002926648550000034
wherein k ism,nRepresenting the element of the information matrix K at the (m, n) position,
Figure FDA0002926648550000035
the dimension of the matrix A1 is N, and the dimension of the matrix A2 is N-1; matrix A1 and matrix A2 are collectively denoted by matrix A and by NARepresents the dimension of matrix a;
step 42: constructing a filter matrix
Figure FDA0002926648550000036
Filtering the matrix A, wherein the expression of a filtering vector C is as follows:
Figure FDA0002926648550000037
wherein the content of the first and second substances,
Figure FDA0002926648550000041
is a conjugate transpose of C; u-1, 2, …, NA，
Figure FDA0002926648550000042
Step 43: filtering the matrix A, and recording the filtered matrix A as the matrix AF：
AF＝FA＝F(AC+AR)＝FAC+FAR
Wherein, the matrix ACThe part of the matrix A without random parameters, the matrix ARThe matrix A is a part containing random parameters; first item FAC＝AC；
Step 44: will matrix AFClockwise rotating 45 degrees to obtain a filtered information matrix K, and recording the filtered information matrix K as a matrix KF：
If (m-n) is an odd number,
Figure FDA0002926648550000043
if (m-n) is an even number,
Figure FDA0002926648550000044
wherein the content of the first and second substances,
Figure FDA0002926648550000045
representation matrix KFOf the element at the (m, n) position,
Figure FDA0002926648550000046
representation matrix A1FAn element at the ((m-N)/2+ (N +3)/4, (m + N)/2) position,
Figure FDA0002926648550000047
representation matrix A2FOf the elements at the ((m-N-1)/2+ (N +3)/4, (m + N-1)/2) positions, matrix A1FIs a filtered matrix A1, matrix A2FIs a filtered matrix a 2;
matrix KFThe size of (2) is still NxN, which is consistent with the information matrix K, and no information dimension is lost before and after filtering.
4. The image reconstruction method of photoacoustic tomography based on full matrix filtering and time reversal operator as claimed in claim 1, characterized in that: in the step 5, the time reversal imaging process is performed according to the filtered information matrix as follows:
step 51: for matrix KFSingular value decomposition:
Figure FDA0002926648550000048
wherein, UFIs a unitary matrix of order NxN, the row vector of which is the matrix KFLeft singular vectors of (d); lambdaFDetermining an NxN order diagonal matrix for the semi-positive matrix;
Figure FDA0002926648550000049
is a unitary matrix of order N x N,
Figure FDA00029266485500000410
is a VFBy conjugate transposition of (V)FThe row vector of is the matrix KFRight singular vectors of (d); lambdaFThe diagonal elements are the matrix KFThe singular value of (a);
step 52: calculating a direct wave transmission matrix G between an arbitrary point (X, Z) on an imaging plane and an ultrasonic transducer array, wherein an element G at the position of (m, n) in the matrix Gm,nComprises the following steps:
Figure FDA0002926648550000051
wherein r ism,nIs the distance between the nth transducer unit and the mth transducer unit in the ultrasonic transducer array;
step 53: the eigenvectors corresponding to the largest M singular values and the direct wave transmission matrix G apply a time reversal operator, the ith singular value
Figure FDA0002926648550000052
The photoacoustic signal intensity at depth Z is:
Figure FDA0002926648550000053
wherein the content of the first and second substances,
Figure FDA0002926648550000054
represents VFThe ith row of singular vectors; superposing the strength of the photoacoustic signals correspondingly obtained by the maximum M singular values to obtain a photoacoustic image at the depth Z
Figure FDA0002926648550000055
Figure FDA0002926648550000056
Is IFThe mth element of the vector, which represents the preliminary imaging value for the location on the imaging plane corresponding to the mth transducer element.
5. The image reconstruction method of photoacoustic tomography based on full matrix filtering and time reversal operator as claimed in claim 4, characterized in that: in step 6, the generation and utilization process of the correction parameters is as follows:
step 61: imaging correction value corresponding to mth transducer unit at depth Z
Figure FDA0002926648550000057
Comprises the following steps:
Figure FDA0002926648550000058
wherein the content of the first and second substances,
Figure FDA0002926648550000059
is gm,nConjugation of (1);
step 62: using correction values for imaging
Figure FDA00029266485500000510
To the preliminary imaging value
Figure FDA00029266485500000511
Correction is made so as to obtain a final imaging value
Figure FDA00029266485500000512
Figure FDA00029266485500000513
I.e. the final imaging result is filtered by the imaging result IFAnd a correction matrix IRThe corresponding elements are obtained by dividing.
CN201810706816.4A 2018-07-02 2018-07-02 Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator Active CN109091109B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810706816.4A CN109091109B (en) 2018-07-02 2018-07-02 Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810706816.4A CN109091109B (en) 2018-07-02 2018-07-02 Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator

Publications (2)

Publication Number Publication Date
CN109091109A CN109091109A (en) 2018-12-28
CN109091109B true CN109091109B (en) 2021-04-20

Family

ID=64845298

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810706816.4A Active CN109091109B (en) 2018-07-02 2018-07-02 Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator

Country Status (1)

Country Link
CN (1) CN109091109B (en)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111435528A (en) * 2019-01-15 2020-07-21 中国科学院金属研究所 Laser ultrasonic visual image quality improvement processing method
CN109674490B (en) * 2019-01-17 2021-09-10 南京大学深圳研究院 Ultrasonic-guided photoacoustic microscope imaging method with low reflection artifact
CN109864707B (en) * 2019-01-17 2021-09-07 南京科技职业学院 Method for improving photoacoustic tomography resolution ratio under limited viewing angle condition
CN110881947B (en) * 2019-12-06 2022-07-05 北京信息科技大学 Optical coherence tomography imaging method

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101173869A (en) * 2006-05-05 2008-05-07 尤洛考普特公司 Method and apparatus for diagnosing a mechanism
US20150249493A1 (en) * 2012-09-24 2015-09-03 Telefonaktiebolaget L M Ericsson (Publ) Prefiltering in MIMO Receiver
CN105181805A (en) * 2015-09-30 2015-12-23 中国计量学院 Multi-filtering ultrasonic imaging method based on time-reversal MUSIC
CN107710014A (en) * 2015-05-12 2018-02-16 法国电力公司 The method and apparatus detected is propagated using ripple
CN107861517A (en) * 2017-11-01 2018-03-30 北京航空航天大学 The online trajectory planning method of guidance of great-jump-forward reentry vehicle based on linear pseudo- spectrum

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20040059265A1 (en) * 2002-09-12 2004-03-25 The Regents Of The University Of California Dynamic acoustic focusing utilizing time reversal
US20040128081A1 (en) * 2002-12-18 2004-07-01 Herschel Rabitz Quantum dynamic discriminator for molecular agents
RU2009108329A (en) * 2006-08-10 2010-09-20 Конинклейке Филипс Электроникс Н.В. (Nl) DEVICE AND METHOD FOR PROCESSING THE AUDIO SIGNAL
US8660328B2 (en) * 2008-06-27 2014-02-25 Wolfram R. JARISCH High efficiency computer tomography
CN102068277B (en) * 2010-12-14 2013-03-13 哈尔滨工业大学 Method and device for observing photoacoustic imaging in single-array element and multi-angle mode based on compressive sensing
CN102867294B (en) * 2012-05-28 2015-06-17 天津大学 Fourier-wavelet regularization-based coaxial phase contrast image restoration method
CN105785566B (en) * 2016-03-31 2017-12-19 南京大学 A kind of method that utilization space optical modulator improves photoacoustic imaging limited perspective

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101173869A (en) * 2006-05-05 2008-05-07 尤洛考普特公司 Method and apparatus for diagnosing a mechanism
US20150249493A1 (en) * 2012-09-24 2015-09-03 Telefonaktiebolaget L M Ericsson (Publ) Prefiltering in MIMO Receiver
CN107710014A (en) * 2015-05-12 2018-02-16 法国电力公司 The method and apparatus detected is propagated using ripple
CN105181805A (en) * 2015-09-30 2015-12-23 中国计量学院 Multi-filtering ultrasonic imaging method based on time-reversal MUSIC
CN107861517A (en) * 2017-11-01 2018-03-30 北京航空航天大学 The online trajectory planning method of guidance of great-jump-forward reentry vehicle based on linear pseudo- spectrum

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Photoacoustic imaging in scattering media by combining a correlation matrix filter with a time reversal operator;WEI RUI等;《Optics Express》;20170918;第22840-22850页 *

Also Published As

Publication number Publication date
CN109091109A (en) 2018-12-28

Similar Documents

Publication Publication Date Title
CN109091109B (en) Image reconstruction method of optimized photoacoustic tomography based on full-matrix filtering and time reversal operator
Matrone et al. High frame-rate, high resolution ultrasound imaging with multi-line transmission and filtered-delay multiply and sum beamforming
Song et al. Comb-push ultrasound shear elastography (CUSE) with various ultrasound push beams
Tong et al. Multi-transmit beam forming for fast cardiac imaging-a simulation study
JP5888833B2 (en) Apparatus and method for ultrasonic synthesis imaging
Sanchez et al. An ultrasonic imaging speckle-suppression and contrast-enhancement technique by means of frequency compounding and coded excitation
Nair et al. Robust short-lag spatial coherence imaging
US11684345B2 (en) System and method for visualization of tissue microvasculature using ultrasound
JP5905080B2 (en) Enhanced ultrasound imaging using qualified areas in overlapping transmit beams
KR20210042907A (en) Method and system for non-invasive characterization of heterogeneous media using ultrasound
Abadi et al. Frequency-sum beamforming for passive cavitation imaging
Hemmsen et al. Tissue harmonic synthetic aperture ultrasound imaging
Chen et al. Doppler-based motion compensation strategies for 3-D diverging wave compounding and multiplane-transmit beamforming: A simulation study
Kou et al. Grating lobe reduction in plane-wave imaging with angular compounding using subtraction of coherent signals
Toulemonde et al. Thomson’s multitaper approach combined with coherent plane-wave compounding to reduce speckle in ultrasound imaging
Lou et al. Zero-phase filtered delay multiply and sum in ultrasound computed tomography
Madhavanunni et al. Beam multiply and sum: A tradeoff between delay and sum and filtered delay multiply and sum beamforming for ultrafast ultrasound imaging
Ghoshal et al. Improved scatterer property estimates from ultrasound backscatter using gate-edge correction and a pseudo-Welch technique
Shin et al. Improving the quality of ultrasound images acquired using a therapeutic transducer
JP2006512129A (en) Small defect detection in medical ultrasound imaging
Jeong et al. High-spatial-resolution, instantaneous passive cavitation imaging with temporal resolution in histotripsy: a simulation study
Gong et al. Simulation studies of filtered spatial compounding (FSC) and filtered frequency compounding (FFC) in synthetic transmit aperture (STA) imaging
Islam et al. Evaluation of TR-MUSIC algorithm efficiency in detecting breast microcalcifications
Nili et al. Field of View and Resolution Improvement in Coprime Sparse Synthetic Aperture Ultrasound Imaging
EP3125770B1 (en) Ultrasonic contrast agent detection and imaging

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant