High-frequency photoacoustic computed tomography method based on sub-nyquist sampling

Abstract

The present invention relates to the technical field of biomedical photoacoustic imaging. Disclosed is a high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling. The method comprises: receiving photoacoustic signals from a target by means of an ultrasound transducer; passing the photoacoustic signals through a band-pass filter to filter out noise and frequency spectra of no interest outside the frequency spectrum of interest from the photoacoustic signals; sampling the photoacoustic signals that have passed through the band-pass filter, then flipping and copying the frequency spectrum of interest from a high-frequency range to a low-frequency range to obtain a flipped frequency spectrum; flipping and folding the flipped frequency spectrum back to the high-frequency range to obtain a recovered frequency spectrum that is identical to the frequency spectrum of interest; and performing photoacoustic computed tomography by means of the recovered frequency spectrum to obtain an image of the target. By introducing the high-frequency PACT method based on sub-Nyquist sampling, the present invention overcomes the limitations of existing PACT systems, thereby achieving high-frequency PACT with a relatively low sampling rate while obtaining relatively high spatial resolution.

Description

A high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling Technical Field

[0001] The present invention relates to the technical field of biomedical photoacoustic imaging, and in particular to a high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling. Background Art

[0002] Photoacoustic computed tomography (PACT) is a biomedical imaging modality that has rapidly developed over the past two decades. Based on the photoacoustic effect, PACT enables high-temporal and spatial resolution imaging of deep biological tissues, thus finding widespread application in a range of biomedical fields, including neurology, oncology, rheumatology, surgical navigation, and other medical imaging areas. High-resolution PACT can reveal the finer structures within biological tissues. The spatial resolution of a PACT system is affected by multiple system factors and is closely related to the characteristics of the ultrasound transducer, including center frequency, bandwidth, aperture size, and viewing angle. Improving spatial resolution, and thus improving PACT image quality, can be achieved simply by increasing the center frequency of the ultrasound transducer. According to the Nyquist sampling theorem, high-frequency PACT requires a data acquisition (DAQ) system with a high-sampling-rate analog-to-digital converter (ADC). In practice, the highest sampling rates supported by existing PACT systems typically range from 40 MHz to 60 MHz, resulting in a maximum supported Nyquist sampling frequency of no more than 30 MHz.

[0003] In the prior art, DAQ systems with sampling rates exceeding 60 MHz are typically used to implement high-frequency PACT with a bandwidth exceeding 30 MHz. Pan et al. used an 80 GHz sampling oscilloscope to sample photoacoustic signals detected by a micro-ring transducer array with a -6 dB bandwidth of 175 MHz. Vionnet et al. customized a DAQ system with a 125 MHz sampling rate to sample photoacoustic signals recorded by a 24 MHz transducer. Li et al. used a commercial imaging system with a 128 MHz sampling rate to sample photoacoustic signals received by a linear transducer array with a center frequency of 40 MHz (bandwidth: 33 MHz). However, DAQ systems with such high sampling rates are typically expensive. Inspired by a previously reported delayed-trigger DAQ method for high-frequency ultrasound imaging, Fu et al. proposed a new method called interleaved sampling photoacoustic imaging. This method enables high-frequency PACT at relatively low sampling rates, for example, sampling a 30 MHz transducer at a 41.67 MHz sampling rate. The proposed method requires two acquisitions, precisely offset by half the sampling period, and performed at the same sampling rate (e.g., 41.67 MHz). These two sampled signals are interleaved to form a virtual sampled signal equivalent to a signal sampled at twice the sampling rate (e.g., 83.33 MHz). Using this method, they achieved 63 μm axial resolution and 91 μm lateral resolution for high-frequency photoacoustic imaging.

[0004] Summary of the Invention

[0005] To solve the above technical problems, the present invention provides a high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling. By introducing a high-frequency PACT method based on sub-Nyquist sampling, the limitations of existing PACT systems are overcome, and high-frequency PACT is achieved at a relatively low sampling rate to obtain higher spatial resolution.

[0006] In order to solve the above technical problems, the present invention adopts the following technical solutions:

[0007] A high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling, comprising:

[0008] The photoacoustic signal of the target is received by the ultrasonic transducer; the center frequency of the photoacoustic signal is f c , bandwidth is B;

[0009] Passing the photoacoustic signal through a bandpass filter to filter out noise outside the spectrum of interest and spectrum of no interest in the photoacoustic signal;

[0010] Sub-Nyquist sampling is performed on the photoacoustic signal after passing through the bandpass filter, and the spectrum of interest will be shifted to 0 to Frequency interval, get the spectrum after translation; where m is the number of translations; at this time the sampling rate The shifted spectrum is restored back to the frequency interval where the spectrum of interest is located to obtain the restored spectrum; at this time, the equivalent sampling rate The recovered spectrum is exactly the same as the spectrum of interest. At this point, the actual sampling rate of the data acquisition system is But the recovered signal (recovered spectrum) is equivalent to using The signal is sampled at a sampling rate of

[0011] Photoacoustic computed tomography is performed by restoring the spectrum to obtain an image of the target.

[0012] Furthermore, the passband of the bandpass filter is arrive

[0013] Furthermore, the number of translations m = 1, the signal bandwidth Sampling rate The frequency range of the spectrum of interest is to So that the use Sub-Nyquist sampling can achieve the same sampling rate as Nyquist sampling effect on the sampling rate.

[0014] Further, in arrive The frequency spectrum of interest is between 0 and The frequency spectrum after the frequency interval is shifted Symmetrically, when restoring the shifted spectrum, the shifted spectrum is moved along Flip fold back to frequency range.

[0015] Furthermore, on a data acquisition system with a maximum sampling rate of 62.5 MHz, a photoacoustic signal with a center frequency of 30 MHz and a bandwidth of 20 MHz can be subjected to sub-Nyquist sampling at a sampling rate of 41.67 MHz, whereas a photoacoustic signal with a center frequency of 30 MHz and a bandwidth of 20 MHz cannot be subjected to Nyquist sampling at a sampling rate of 62.5 MHz. This is because a photoacoustic signal with a center frequency of 30 MHz and a bandwidth of 20 MHz has a maximum frequency of 40 MHz. Normal Nyquist sampling requires a sampling rate of 40*2=80 MHz, while sub-Nyquist sampling only requires a sampling rate of 30*4 / 3=40 MHz. Therefore, on a data acquisition system with a maximum sampling rate of 62.5 MHz, this signal can be subjected to sub-Nyquist sampling, but not Nyquist sampling.

[0016] Compared with the prior art, the beneficial technical effects of the present invention are:

[0017] The present invention proposes an imaging technology based on sub-Nyquist sampling, which can achieve high-frequency PACT at a relatively low sampling rate, thereby reducing the need for expensive high-sampling rate DAQ systems. Compared with traditional imaging methods, the imaging method of the present invention shows significant advantages in spatial resolution. In order to suppress the ringing artifacts caused by this technology, the present invention incorporates another signal that samples the low-frequency components into the signal collected by the sampling method. Under the condition of a relatively low sampling rate (i.e., 41.67MHz), the present invention uses this imaging technology to achieve an axial resolution of 22μm. It is understood that when the sampling rate does not exceed 60MHz, the above resolution is the highest axial resolution ever achieved by PACT. Therefore, the present invention provides a practical solution for high-frequency PACT under limited sampling rates. BRIEF DESCRIPTION OF THE DRAWINGS

[0018] The above and other objects, features and advantages of the present invention will become more apparent through the following description of the embodiments of the present invention with reference to the accompanying drawings, in which:

[0019] Figure 1 schematically illustrates the basic principles of the high-frequency PACT method based on sub-Nyquist sampling proposed in this invention. (a) is the flowchart of the invention; (b) is the original analog signal spectrum; (c) is the filtered analog signal spectrum, which is also the spectrum of interest; (d) is the sampled digital signal spectrum, where the high-frequency spectrum of interest is copied and flipped to a low-frequency range; and (e) is the high-frequency recovered spectrum obtained after flipping and recovery.

[0020] Figure 2 schematically illustrates the performance of the imaging method of the present invention verified by numerical simulation. (a) is a schematic diagram of the simulation setup, showing a 20 μm diameter microsphere placed 8 mm from the linear transducer array; (b) is the gold standard simulation obtained using sampling method 6 in Table 1; (c) to (g) are the simulation results obtained using sampling methods 1 to 5 in Table 1, respectively; (h) is the result of adding (c) and (d); (i) and (j) are comparisons of the axial and lateral resolutions of different sampling methods.

[0021] Figure 3 schematically illustrates the imaging results of a phantom experiment using fine hair in the present invention. (a) through (e) are the imaging results obtained using sampling methods 1 through 5 in Table 1, respectively; (f) is the sum of (a) and (b).

[0022] Figure 4 schematically illustrates the imaging results of an in vivo experiment using a mouse ear, as described in the present invention. (a) is a photograph of a mouse ear; (b) to (f) are images of blood vessels (indicated by arrows) obtained using sampling methods 1 to 5 in Table 1; (g) is the sum of (b) and (c). The inset within the dashed box is a magnified image of the blood vessels. (h) and (i) are the intensity distributions of the blood vessels along the z-axis and x-axis, respectively, representing the full width at half maximum. DETAILED DESCRIPTION

[0023] A preferred embodiment of the present invention will be described in detail below with reference to the accompanying drawings.

[0024] Sub-Nyquist sampling, also known as bandpass sampling, allows sampling a bandpass signal at a sampling rate less than twice the highest frequency component and is now widely used in ultrasound imaging to reduce computational complexity. The basic idea of ​​the present invention is based on the fact that the frequency response of high-frequency ultrasound transducers is usually band-limited.

[0025] Assume there is a center frequency f c , a bandpass signal with a bandwidth of B can be obtained with f without aliasing according to the sub-Nyquist sampling theorem. s The bandpass signal is sampled at a sampling rate of f. s The following conditions are met:

[0026] Where m is the value that satisfies f s Any positive integer ≥ 2B.

[0027] Formula (1) can be rewritten as:

[0028] The smaller the bandwidth, the larger the value of m can be, and the smaller the sampling rate f can be. s Since bandwidth is crucial for high-quality PACT, a smaller m should be chosen. For example, in the present invention, m = 1. In this case, the highest bandwidth supported by sub-Nyquist sampling is The corresponding f s yes The basic principle of the high-frequency PACT method based on sub-Nyquist sampling in the present invention is shown in FIG1 (a), and the sub-Nyquist sampling parameters used are: m=1, First, a bandpass filter is used to filter out the noise and uninteresting spectra outside the spectrum range of interest, as shown in Figure 1(b) and Figure 1(c). Then, the filtered photoacoustic signal (with a center frequency of f c , the bandwidth is ) for sampling, and the sampling rate is consistent with the highest frequency of the filtered photoacoustic signal

[0029] After sampling, it is in the high frequency range ( arrive ) will be flipped and copied to the low frequency range (0 to ), as shown in (d) in Figure 1. Then, the spectrum copied to the low-frequency interval is along Flip and fold back to the high frequency range, as shown in (e) in Figure 1. At this time, the sampling rate has been doubled, that is, In Figure 1(e), the recovered spectrum is exactly the same as the spectrum of interest shown in Figure 1(c). Finally, the recovered signal is input into the algorithm for PACT image reconstruction.

[0030] Figure 1 demonstrates that a signal with a bandwidth of 67% can be sampled at a sampling rate equal to the signal's highest frequency and then accurately recovered through post-processing. The recovered signal is identical to the signal obtained by directly sampling the original signal at a sampling rate equal to twice the signal's highest frequency (Nyquist sampling criterion). The key to the sampling method of this invention is ensuring that the original analog signal (photoacoustic signal) is first filtered by an analog bandpass filter before ADC sampling to eliminate noise and uninteresting spectra outside the spectrum of interest. Otherwise, these noise and uninteresting spectra will be aliased into the replicated spectrum after ADC sampling. In this case, the original spectrum of interest cannot be accurately recovered.

[0031] Table 1 shows different sampling methods used for performance comparison in the present invention.

[0032] Table 1

[0033] Example 1: Numerical simulation

[0034] The numerical simulation was performed by imaging a microsphere with a diameter of 20 μm under a uniform background (sound velocity: 1500 m / s). The photoacoustic signal of the microsphere was generated using an analytical model. The generated photoacoustic signal was received by a linear transducer array with 256 array elements and a 15 mm aperture. The center frequency of each array element was 30 MHz and the bandwidth was 20 MHz. The bandwidth characteristics of the transducer array were simulated by applying a Gaussian bandpass filter to the generated signal. The microsphere was placed at a distance of 8 mm from the transducer array. The received photoacoustic signal was first sampled at 1 GHz to be used as the original signal in the numerical simulation. The signals were individually filtered using different filter configurations in Table 1 and downsampled to different sampling rates to compare the performance of the present invention with other methods. In this numerical simulation, it was assumed that the highest sampling rate supported by the DAQ was 62.5 MHz. Sampling rates of 83.33 MHz and 125 MHz were used to simulate interleaved sampling based on 41.67 MHz and 62.5 MHz, respectively. A first-order Butterworth high-pass filter and a third-order Butterworth low-pass filter were combined as a bandpass filter. The PACT image was reconstructed using the filtered back-projection algorithm.

[0035] The numerical simulation setup is shown in Figure 2(a). Microspheres with a diameter of 20 μm are placed 8 mm from the linear transducer array. Figure 2(b) shows the imaging results of the microspheres obtained using sampling method 6 in Table 1. Although the reconstructed microspheres shown in Figure 2(b) are not circular due to the limited transducer bandwidth and viewing angle, Figure 2(b) is still considered the gold standard in this numerical simulation. Because the non-circularity of the reconstructed microspheres is caused by the hardware, the impact on all sampling methods is consistent and does not affect the comparison of the different sampling methods here. Figures 2(c) to (g) show the imaging results of the microspheres obtained using sampling methods 1 to 5 in Table 1. The negative artifacts shown in Figures 2(b) to (g) are caused by the limited transducer bandwidth and viewing angle. The intensity distribution along the dashed line shown in Figure 2(b) is shown in Figures 2(i) and (j), respectively. The full width at half maximum (FWHM) value represents the axial (z-axis) and lateral (x-axis) resolution. As can be seen, Figure 2(c) is clearer than Figure 2(e). This is because the signal obtained using sampling method 1 has higher frequency components than the signal obtained using sampling method 3. However, compared with Figure 2(b), Figure 2(c) shows more obvious ringing artifacts (indicated by arrows). This is because sampling method 1 removes low-frequency signals and only contains high-frequency signals. Combining Figures 2(d) and 2(c) yields a reconstruction result that contains both low-frequency (sampling method 2) and high-frequency components (sampling method 1), as shown in Figure 2(h). Compared with Figure 2(c), Figure 2(h) better suppresses ringing artifacts (indicated by arrows). It can be seen that the results obtained by the imaging method of the present invention (Figure 2(c)) are as clear as those obtained based on interleaved sampling (Figures 2(f) and (g)). The high similarity between Figures 2(f) and (g) indicates that the spectrum of the microspheres is primarily within 35 MHz. These findings are strongly supported by the axial and lateral resolution. As shown in (i) and (j) of Figure 2, the conventional sampling method (sampling method 3) produces the worst imaging results, with the microspheres being the widest. Numerical simulation results show that the proposed high-frequency PACT method based on sub-Nyquist sampling (sampling method 1) outperforms the conventional imaging method (sampling method 3) in spatial resolution and is comparable to the interleaved sampling method (sampling methods 4 and 5). To suppress ringing artifacts, it is necessary to add another sampling method (sampling method 2) to the proposed sampling method (sampling method 1) that only samples low-frequency components.

[0036] Example 2: Phantom Experiment

[0037] In order to further demonstrate the performance of the imaging method of the present invention, first, a hair with a diameter of 10 μm to 20 μm was removed from the human forearm and embedded in an agar model. Then, the imaging method of the present invention was used to image this model. The imaging system uses the Vantage 256 system (Verasonics Inc., WA, USA), which is designed for real-time multi-channel (up to 256) data acquisition and processing, and supports a maximum sampling rate of 62.5 MHz. The system is connected to a commercial 256-channel linear transducer array (MS400, Verasonics Inc., Canada) with a center frequency of 30 MHz, a bandwidth of 20 MHz, and an aperture of approximately 15 mm. The detected photoacoustic signals are filtered separately using a variety of filter configurations and sampled at different rates (sampling methods 1 to 5 in Table 1) to compare the performance of the imaging method of the present invention with other imaging methods. The other setting parameters are the same as those in the numerical simulation. Figures 3 (a) to (e) are the hair imaging results obtained using different sampling methods (1 to 5 in Table 1). As expected, the imaging method of the present invention (see (a) in FIG3 ) provides finer structures, but has more ringing artifacts (indicated by arrows) compared to the conventional imaging method (see (c) in FIG3 ). By fusing (b) in FIG3 with (a), the ringing artifact shown in (a) in FIG3 is well suppressed, as shown in (f) in FIG3 . The interleaved sampling method (see (d) and (e) in FIG3 ) also provides clearer images compared to the conventional imaging method (see (c) in FIG3 ). The ringing artifacts (indicated by arrows) shown in (e) in FIG3 are mainly caused by the limited transducer bandwidth and viewing angle. It can be seen that the combined imaging method (see (f) in FIG3 ) performs better in terms of spatial resolution and image quality.

[0038] Example 3: In vivo experiment

[0039] In order to verify the performance of the imaging method of the present invention, the present invention conducted an in vivo experiment on the ears of nude mice. The mouse ear blood vessels shown in Figure 4 (a) were imaged using different sampling methods (sampling methods 1 to 5 in Table 1), and the results are shown in Figures 4 (b) to (f). As expected, the imaging method of the present invention (see Figure 4 (b)) provides higher spatial resolution, but ringing artifacts will appear compared with the traditional imaging method (see Figure 4 (d)). By fusing Figures 4 (c) and (b), the ringing artifacts shown in Figure 4 (b) are well suppressed, as shown in Figure 4 (g). The performance of the combined imaging method (see Figure 4 (g)) is comparable to that of the interleaved sampling imaging method (see Figures 4 (e) and (f)). These results are further verified by the FWHM values ​​of the blood vessels shown in Figures 4 (h) and (i).

[0040] In practice, the rising and falling edges of analog bandpass filters typically have finite slopes, which means that noise and unwanted spectra outside the spectrum of interest cannot be completely eliminated. As a result, the recovered spectrum is contaminated to some extent. The ringing artifacts shown in Figure 2(h) and Figure 3(f) arise not only from the limited transducer bandwidth and viewing angle, but also from the finite slopes of the rising and falling edges of the analog bandpass filters. High-performance analog bandpass filters are expected to reduce ringing artifacts.

[0041] It will be apparent to those skilled in the art that the present invention is not limited to the details of the exemplary embodiments described above and that the invention can be embodied in other specific forms without departing from the spirit or essential characteristics of the invention. Therefore, the embodiments should be considered in all respects as illustrative and non-limiting, and the scope of the invention is defined by the appended claims rather than the foregoing description. It is intended that all variations within the meaning and range of equivalents of the claims be embraced herein, and any reference signs in the claims should not be construed as limiting the claims to which they relate.

[0042] In addition, it should be understood that although this specification is described in terms of implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each embodiment can also be appropriately combined to form other implementation methods that can be understood by those skilled in the art.

Claims

1. A high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling, comprising: The photoacoustic signal of the target is received by the ultrasonic transducer; the center frequency of the photoacoustic signal is f c , bandwidth is B; Passing the photoacoustic signal through a bandpass filter to filter out noise outside the spectrum of interest and spectrum of no interest in the photoacoustic signal; Sub-Nyquist sampling is performed on the photoacoustic signal after passing through the bandpass filter, and the spectrum of interest will be shifted to 0 to Frequency interval, get the spectrum after translation; where m is the number of translations; at this time the sampling rate The shifted spectrum is restored back to the frequency interval of the spectrum of interest to obtain the restored spectrum; at this time, the equivalent sampling rate And the recovered spectrum is exactly the same as the spectrum of interest; Photoacoustic computed tomography is performed by restoring the spectrum to obtain an image of the target.

2. The high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling according to claim 1, characterized in that: The passband of the bandpass filter is arrive 3. The high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling according to claim 1 or 2, characterized in that: The number of translations m = 1, the signal bandwidth Sampling rate The frequency range of the spectrum of interest is to So that the use Sub-Nyquist sampling can achieve the same sampling rate as Nyquist sampling effect on the sampling rate.

4. The high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling according to claim 3, characterized in that: in arrive The frequency spectrum of interest is between 0 and The frequency spectrum after the frequency interval is shifted Symmetrically, when restoring the shifted spectrum, the shifted spectrum is moved along Flip fold back to frequency range.

5. The high-frequency photoacoustic computed tomography method based on sub-Nyquist sampling according to claim 3, characterized in that: On a data acquisition system with a maximum sampling rate of 62.5 MHz, sub-Nyquist sampling can be performed on a photoacoustic signal with a center frequency of 30 MHz and a bandwidth of 20 MHz at a sampling rate of 41.67 MHz, while Nyquist sampling cannot be performed on a photoacoustic signal with a center frequency of 30 MHz and a bandwidth of 20 MHz at a sampling rate of 62.5 MHz.

Images