Two-step proximal gradient descent algorithm for photoacoustic signal unmixing.

Photoacoustic microscopy uses multiple wavelengths to measure concentrations of different absorbers. The speed of sound limits the shortest wavelength switching time to sub-microseconds, which is a bottleneck for high-speed broad-spectrum imaging. Via computational separation of overlapped signals, we can break the sound-speed limit on the wavelength switching time. This paper presents a new signal unmixing algorithm named two-step proximal gradient descent. It is advantageous in separating multiple wavelengths with long overlapping and high noise. In the simulation, we can unmix up to nine overlapped signals and successfully separate three overlapped signals with 12-ns delay and 15.9-dB signal-to-noise ratio. We apply this technique to separate three-wavelength photoacoustic images in microvessels. In vivo results show that the algorithm can successfully unmix overlapped multi-wavelength photoacoustic signals, and the unmixed data can improve accuracy in oxygen saturation imaging.

Introduction

Optical-resolution photoacoustic microscopy (OR-PAM) is a non-invasive biomedical imaging technique that provides subcellular resolution based on optical-absorption contrast [1], [2]. Multi-wavelength photoacoustic excitations of the same target can probe various molecules [3], [4], [5], [6], [7] and quantify hemoglobin concentration [8], [9], [10], [11], [12], oxygen saturation (sO2) [13], [14], [15], and blood flow [16], [17], [18]. Multiple excitations usually sacrifice the scanning speed, which is unfavorable for acquiring dynamic information and mitigating motion artifacts. To implement fast broad-spectrum OR-PAM, the time delay between adjacent excitations must be sufficiently short [19], [20], [21], [22], [23], [24].

Several methods have been developed to achieve this goal. Sequentially triggering multiple lasers can flexibly control the time interval between excitations [25], [26]. Another method uses electro-optical modulators (EOMs) [27], [28] to switch among different light paths. Nevertheless, both methods increase costs and complicate the laser systems. The stimulated Raman scattering (SRS) effect in optical fibers can produce wavelength shift and pulse delay at a relatively low cost and, thus, become a practical choice for multi-wavelength OR-PAM [2], [29], [30], [31], [32]. However, the shortest time to switch among wavelengths or pulses in all these methods must be longer than the time duration of an individual signal. Otherwise, different photoacoustic signals would be overlapped.

Two approaches have been developed to unmix the overlapped signals. One is a frequency-domain approach [33] that uses partial data to solve an optimization problem. Experimental results show that the method can unmix two overlapped signals but becomes inaccurate when the two signals are less than 30-ns apart or the signal-to-noise ratio (SNR) is less than 31 dB [33]. Another method uses a deep learning model to separate overlapped signals [34] and has demonstrated signal unmixing with a 30-ns delay. However, it is costly to acquire a large amount of training data.

Both methods only implement separating two signals. It remains an open problem to separate three or more overlapped signals with short overlapping and high noises.

We present a two-step proximal gradient descent (TS-PGD) algorithm to separate overlapped signals. It does not require training data and thus avoids extra hardware costs. The algorithm uses the entire photoacoustic signal for fitting and directly outputs accurate results even in short-delay and low-SNR conditions. TS-PGD can separate multi-wave mixed signals with high accuracy.

Using TS-PGD, we simulate photoacoustic signal unmixing in various conditions. Then, we use a three-wavelength photoacoustic microscope for experimental validation. The time delays between adjacent laser pulses are ~30 ns. In phantom validation, we can separate three overlapped signals when SNR is as low as 19.1 dB. In in-vivo experiments, we acquire microvascular images of partially overlapped signals at three wavelengths. We can successfully unmix the three overlapped signals and use them to calculate sO2.

Method

We use a laser pulse train to excite one absorber and generate a series of 1D photoacoustic A-line signals. If the time delay between adjacent pulses is too short, these A-line signals will be overlapped. We can assume these A-line signals have the same waveform within the linear range, and the only difference is the amplitude. Fig. 1(a) and (b) show individual and overlapped signals. We can model an overlapped photoacoustic signal as:whereis the unknown amplitude factor of the i A-line signal,represents the waveform of an individual A-line signal,is the time delay and can be measured, and n is the number of non-overlapped signals. The goal is to estimateandfrom the overlapped signal. Because bothandare unknown, it becomes a nonlinear optimization problem. We develop an algorithm named two-step proximal gradient descent to solve the nonlinear model with fast convergence, high accuracy, and great tolerance to noises.

Fig. 1

Schematic of TS-PGD algorithm: (a) Ground-truth PA signals excited with single pulses; (b) Overlapped signal excited with a short-delayed pulse train; (c) TS-PGD process; (d) Intermediate fitting result; (e) Fitted result; (f) Separated signals.

Two-step proximal gradient descent algorithm

As shown in Fig. 1(c), TS-PGD algorithm solves the nonlinear optimization problem in two steps. In the first step, we linearize Eq. (1) to obtain as:

For sampled data, is usually not an integer multiple of the sampling interval , which may cause errors in the unmixing results. We propose a dual-k mode based on linear interpolation to address this problem. Let and are integer times of ,, ≤ < , and , then the discrete form of can be written as:

The dual-k mode uses twoto solve the problem of inaccurate delay estimation. Thereby, Eq. (3) can be linearized as:

The second step is to seek an optimal solution using a proximal gradient descent method. Because thesignal has no delay, we set , , and constantly. Eq. (4) can be written in a matrix format aswhere = [, ,…, ,…, , ], = [, ,…, ,…, , ,,,…,,,…,,], l is the length of a mixed-signal, m is the length of . is a matrix composed of the current , , and , and its detailed information is shown in supplementary Fig. S1.

In iteration, the initial value of is zero, is the difference between the mixed PA signal and the current , and the initial values of and are 1. can be solved via least-squares fitting:

Because Eq. (6) meets the condition of Lipschitz continuity, the PGD algorithm [35], [36] is suitable for searching for the optimal solution. We use the following algorithm to search for a new point :where is the step size. We use a fixed step size without backtracking, i.e., , and is the maximum absolute value of the second-order derivative of . To speed up the optimization, we efficiently select the next starting point referring to the fast iterative shrinkage-thresholding algorithm (FISTA) [37].

After converging, we can obtain new , and based on the , , and. Because the solution of the linearized model (Eq. (2)) does not have to be the optimal solution of the nonlinear model (Eq. (1)), we iterate this process until it converges. Representative results, including the estimated overlapped signal and separated signals, are shown in Fig. 1(d), (e), and (f).

We implement the algorithm in MATLAB (2021a, MathWorks, USA) on a computer (Inter Core i5@1.60 GHz, 16 GB of RAM).

Results and discussion

Simulation

We simulate signal separation under different numbers of overlapped signals, delay times, and SNRs. The simulation data is measured PA signals excited from the black tape by single laser pulses. We shift and sum several simulated signals to produce overlapped signals with different parameters. In simulations, zeros are appended to the signals to change the data length, and different levels of Gaussian noises are added to the signals.

First, we test unmixing different numbers of overlapped signals. After multiplying the signal excited from the black tape by 8 different random numbers within 0.5–2, we acquired 9 individual A-line signals with the same waveform but different intensities. Then we mix them with a 26-ns delay between adjacent signals. The number of signals is set to 2–9. As shown in Supplementary Fig. S2, all overlapped signals are separated successfully. The SNR range between 2-wave to 9-wave mixed signals is 27.5–34.2 dB. Representative separation results are shown in Fig. 2(a)-(c). Fig. 2(d) plots the mean absolute error of the normalized A-line (NA-MAE) between the separated and the original signals, and detailed information about calculating NA-MAE is shown in Supplementary Eqs. (S1) and (S2). The result indicates that the average error between the separated and the original signals is less than . We repeat the simulation with different signals from the black tape by ten times.

Fig. 2

Simulation on multi-wave separation: (a) Original, (b) overlapped, and (c) separated signals in nine-wave separation; (d) NA-MAEs between the separated and original signals versus different numbers of mixed waves, the error bar is the STD of calculating NA-MAEs, and the sample size is ten.

The error bars in Fig. 2(d) give the standard deviation (STD) of the calculated NA-MAEs under different numbers of overlapped signals, and all STDs are less than . Thus, TS-PGD can separate at least nine overlapped signals.

Next, we test unmixing signals with various delays. Because the sampling interval of the data is 2 ns, three waves are mixed with different delay times from 8 ns to 30 ns with a 2-ns interval. The SNR of mixed signals changes within 22.8–30.0 dB. Representative separated results with 12-ns delay are presented in Fig. 3(a)-(c), and all separation results are shown in Supplementary Fig. S3. As shown in Fig. 3(d), for the delays not less than 12 ns, the NA-MAEs are less than The STDs of the unmixing results are less than , indicating that the proposed method is stable.

Fig. 3

Simulations of signal unmixing with different delays: (a)-(c) are the original, overlapped, and separated signals with 12-ns delay; (d) NA-MAEs in delay simulation, the error bar is the STD of calculating NA-MAEs, and the sample size is ten.

In the last simulation, we test signal unmixing with different SNRs. By adding noise with different amplitudes, the SNRs of the mixed signals are set to 13.1, 15.9, 18.9, 21.9, 25.1, 27.9, 31.0, and 33.9 dB. Fig. 4(a)-(c) shows a representative set of separation results when the SNR is 15.9 dB. All separation results are shown in Supplementary Fig. S4. Fig. 4(d) shows the NA-MAEs between the separated and true signals. When SNR is not less than 15.9 dB, the NA-MAEs are less than , and STDs are less than .

Fig. 4

Simulation of signal unmixing at different SNRs: (a)-(c) corresponds to the original, overlapped, and separated signals, the SNR of the overlapped signal is 15.9 dB; (d) NA-MAEs at different SNRs, the error bar is the STD of calculating NA-MAEs, and the sample size is ten.

According to the simulation results, TS-PGD accurately unmixes signals with various numbers of overlapped signals, different delays, and low SNRs. Even when the overlapped number increases to nine, there is still no significant degradation in unmixing accuracy. A short delay and low SNR have a great impact on unmixing quality. Nevertheless, TS-PGD remains effective in many challenging scenarios. Furthermore, for the cases with different delays between adjacent excitations or a single excitation with multiple dominant peaks, TS-PGD can also unmix them with NA-MAE less than . The detailed results are shown in Figs. S5-S8 of the supplementary materials.

Phantom results

To validate the proposed method, we separate three-wavelength overlapped signals generated from black tape. The schematic of the setup is shown in Fig. 5. A pulsed laser beam (532-nm, laser model VPFL-G-20, Spectra-Physics) is divided into three paths, two of which generate 545-nm and 558-nm wavelengths via the SRS effect in polarization-maintaining fiber (PMF1–2, PM-S405-XP, Nufern). Light from three paths goes through three polarization beam splitters (PBS1-3, PBS051, Thorlabs Inc) and a long pass dichroic mirror (DM, T550lpxr, Chroma) to merge into a single-mode fiber. Light from the fiber is delivered to a photoacoustic probe for imaging. Two half-wave plates (HWP1 −2, WPH10E-532, Thorlabs Inc.) are used to adjust the laser polarization and the SRS efficiency. Neutral density filters (NDF1-3, NDC-50C-2, Thorlabs Inc.) are applied in all three paths to adjust the energy of each wavelength. A data acquisition card (ATS9360, Alazar Tech Inc) digitizes the signal with 12-bit resolution at 500 MHz.

Fig. 5

Schematic of OR-PAM with ultrafast three-wavelength excitation. DM, 550-nm long-pass dichroic mirror; FC1-6, fiber coupler; PMF, polarization-maintaining fiber; HWP1-2, half-wave plate; M1-2, mirror; NDF1–3: neutral density filter; PBS1-3, polarizing beam splitter; PMF1-2, polarization-maintaining fiber; SMF, single-mode fiber; UL, ultrasound lens; UT, ultrasonic transducer; WT, water tank.

We acquire overlapped signals and single-pulse-excited signals from the black tape. The time delays between adjacent laser pulses are 30 ns and 36 ns in the overlapped signals, which corresponds to 6-m PMF1 and 13.5-m PMF2. NDF1-3 adjust the pulse energy of each wavelength so that the overlapped signals have different SNRs of 19.1, 22.7, 25.9, 27.8, and 30.2 dB. The overlapped signals and the single-pulse-excited signals are collected ten times. We apply TS-PGD to separate the overlapped signal and compare them with the single-pulse signals. The results are shown in Fig. 6. A set of representative separation results (19.1-dB SNR) are presented in Fig. 6(a)-(c), and all other results are shown in Supplementary Fig. S9.

Fig. 6

Separation of mixed black-tape signals at different SNRs: (a)-(c) corresponds to the original, overlapped, and separated signals, the SNR of the overlapped signal is 19.1 dB; (d) NA-MAE at different SNRs, the error bar is the STD of calculating NA-MAEs, and the sample size is ten.

The separation errors at different SNRs are shown in Fig. 6(d). All NA-MAEs between the separated signals and the original signals are less than , the NA-MAEs are less than when SNRs are higher than 25.9 dB. The standard deviations are less than for all results. The errors of phantom results are higher than the simulation errors. The reasons are that the waveforms excited from different wavelengths are not identical, and the pulse widths from the three paths are also a little different, both of which are adverse to the TS-PGD separation.

In vivo experiment

Nude mice (BALB/c mouse, 27–30 g) were used for in vivo PA imaging. The mice were anesthetized with 1.2% (v/v) vaporized isoflurane and placed on a heating pad (37 °C). A membrane-sealed water tank was placed on top of the ear. US gel was applied between the ear and the membrane for acoustic coupling. In order to collect the overlapped signal for fast imaging, the mouse ear was sequentially excited by 532-nm, 545-nm, and 558-nm pulses at 0, 34, and 68 ns, and the PMF1-2 are 7 and 14 m. The laser pulse repetition rate is 4 kHz,

and the scanning step size is 2.5 µm. The animal ethical committee of the City University of Hong Kong approved all procedures in animal experiments.

We apply the TS-PGD method to separate overlapped signals. One high-resolution image consists of 1000 × 1000 overlapped A-lines. The SNRs of the overlapped signals range from 0.02 to 33.7 dB. TS-PGD can separate these overlapped A-lines in about 10 min. Maximum amplitude projections (MAP) of the separated images are shown in Fig. 7(a)-(c).

Fig. 7

In vivo functional imaging of the mouse ear with three pulses at 532 545, and 558 nm: (a)-(c) Separated PA images of the mouse ear. (d)-(f) sO2 of the mouse ear calculated from the separated data.

One application of separating the multi-wavelength image is to compensate for the fluence error in sO2 imaging [38], [39]. We can compute sO2 images using two or three wavelengths from the separated A-lines, i.e., 532&558 545&558, and 532&545&558 nm. Fig. 7(d)-(f) shows the three sO2 images. Due to scattering or other wavelength-dependent attenuations, the sO2 results are different, particularly in the small vessels. Because 545&558 nm are closer to each other than 532&558 nm, wavelength-dependent errors in the 545&558-nm result are less than that in the 532&558-nm one. We can see the 545&558-nm sO2 values in small vessels are higher and closer to normal physiological values [14], [40], [41] than the 532&558-nm ones. Using 532&545&558 nm, we can further compensate for wavelength-dependent errors using a linear model [38]. Fig. 7(f) shows the sO2 values in both the vein and small vessels are significantly higher than the two-wavelength results.

We compare the three-wavelength sO2 values with the two-wavelength results pixel-by-pixel in comparative scatter plots and histograms. From the comparative scatter plots (Fig. 8(a) and (b)), we can see that the 532&545&558 results have an enhancement in most of the sO2 values compared to the two-wavelength results. From the histograms in Fig. 8(c), the 532&545&558 image has more pixels in sO2 values between 0.65 and 0.75 than the other two results and is closer to the normal physiological conditions [25], [26]. This trend reflects that wavelength unmixing can effectively separate three mixed signals for fluence compensation in sO2 imaging.

Fig. 8

Comparison of sO2 calculation using different wavelengths: (a)-(b) Comparative scatter plots of three sO2 values; (c) Histograms of the distribution of the sO2 values.

Conclusion

In summary, we present an efficient signal separation method for unmixing short-delayed PA signals. Via linearization and PGD optimization, TS-PGD can efficiently and accurately solve the nonlinear signal unmixing problem.

In the numerical simulation, we test the method in different numbers of overlapped signals, delay times, and SNRs. TS-PGD presents a robust performance in unmixing signals with multiple waves (2–9), short-delay (12 ns), and low-SNR (15.9 dB). Phantom experiments further verify the separation accuracy at different noise levels. Three overlapped signals with SNR greater than 19.1 dB can be separated with NA-MAE less than . In in vivo experiments, we separate three-wavelength photoacoustic signals. The unmixed results are used to calculate sO2 and compensate for fluence error. The in vivo results further demonstrate the effectiveness of the signal unmixing method. TS-PGD can accurately unmix overlapped photoacoustic signals in various situations, offering a new approach for high-speed multi-wavelength photoacoustic imaging.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.