The study on the inverse problem of applied current thermoacoustic imaging based on generative adversarial network

Applied Current Thermoacoustic Imaging (ACTAI) is a new imaging method which combines electromagnetic excitation with ultrasound imaging, and takes ultrasonic signal as medium and biological tissue conductivity as detection target. Taking the high contrast advantage of Electrical Impedance Tomography (EIT) and high resolution advantage of ultrasound imaging, ACTAI has broad application prospects in the field of biomedical imaging. Although ACTAI has high excitation efficiency and strong detectable Signal-to-Noise Ratio, yet while under low frequency electromagnetic excitation, it is still a big challenge to reconstruct a high-resolution image of target conductivity. This paper proposes a new method for reconstructing conductivity based on Generative Adversarial Network, and it consists of three main steps: firstly, use Wiener filtering deconvolution to restore the electrical signal output by the ultrasonic probe to a real acoustic signal. Then obtain the initial acoustic source image with filtered backprojection technology. Finally, match the conductivity image with the initial sound source image, which are used as training samples for generating the adversarial network to establish a deep learning model for conductivity reconstruction. After theoretical analysis and simulation research, it is found that by introducing machine learning, the new method can dig out the inverse problem solving model contained in the data, which further reconstruct a high-resolution conductivity image and has strong anti-interference characteristics. The new method provides a new way to solve the problem of conductivity reconstruction in Applied Current Thermoacoustic Imaging.

www.nature.com/scientificreports/ the Applied Current Thermoacoustic Imaging (ACTAI). Instead of using the above-mentioned induced electric field, they utilize the applied current signal with microsecond pulse width to generate an ultrasonic signal with a stronger Signal to Noise Ratio (SNR). Compared with the induced current, the excitation frequency of the applied current is easier to adjust and it has a higher penetration depth, which improves the energy conversion efficiency and detectable SNR. As a new type of thermoacoustic imaging technology, ACTAI still has many problems in the high-resolution reconstruction of conductivity, which need to be further solved. In terms of image reconstruction, the current research mainly focuses on the reconstruction algorithm of sound source image, such as time inversion method, filtered back projection 12 and so on. Although these algorithms can reflect the distribution characteristics of thermoacoustic sources, they cannot obtain the true conductivity distribution. The fundamental reason is that the Signal-to-Noise Ratio of the ultrasound signal is low, while the inverse problem of rebuilding the conductivity image from the ultrasound signal has a strong nonlinearity, which makes the reconstruction process greatly affected by noise, and the reconstruction results of conductivity cannot meet the requirements of functional imaging.
In recent years, inverse problem reconstruction algorithms based on deep learning have been widely studied in the field of medical imaging, such as magnetic resonance super-resolution imaging, low-dose CT imaging, and so on. From the existing research, we find that the deep learning method has such unique advantages as high degree of nonlinearity, fast solving speed, easy access to simulation training dataset, etc. in resolving large-scale nonlinear reconstruction problems. Generative Adversarial Network [13][14][15] (GAN) can reconstruct image based on limited data, remove the view artifact 16 and improve the image SNR 17 , improve the image similarity and get higher image quality. ACTAI, as a complex Multi-Physics Coupling Imaging mode, its inverse conductivity reconstruction problem precisely meets the above characteristics. Therefore, this paper proposes to apply the deep learning network to conductivity reconstruction in ACTAI in order to obtain the real conductivity distribution image.

Method
The implementation of the new method consists of the following steps: first, the electrical signal measured by the ultrasonic probe is preprocessed by Wiener filter deconvolution 18,19 to obtain the original acoustic signal emitted by the measured sample. Then the filtered back projection (FBP) is used to reconstruct the sound source image, match the reconstructed sound source image with the conductivity image of the measured sample, and use it as a training sample to train a Generative Adversarial Network. Finally, the trained network is used to process the sound source image to obtain high-resolution conductivity images. This method comprehensively considers the characteristics of the physical model as well as the data-driven model, and constructs a deep neural network specifically for the ACTAI inverse problem, which can effectively improve the imaging resolution and imaging speed.
In ACTAI, the duration of the pulse current is much shorter than the heat conduction time inside the sample, so only the thermoacoustic signal caused by the adiabatic expansion of biological tissues needs to be considered. The thermoacoustic pressure p(r, t) generated at field point r can be described as below: Refer to (1), c s represents the speed of sound, C p and β represent the specific heat capacity and volume expansion coefficient of the target body respectively, while r and r ′ respectively indicate the position of the ultrasound probe and the sound source. H(r ′ , t) is the Joule heat absorbed by the biological tissue per unit time and volume, and it can be expressed as H(r ′ , t) = H(r ′ )I(t) , among which, H(r ′ ) indicates the spatial heat absorption distribution, and I(t) is the time domain distribution of pulse excitation intensity. Using the integral solution of Green function, p(r, t) can be further expressed as below in (2): Refer to (2), represents the integral domain containing all sound sources, δ(t−|r ′ −r|/c s ) is the Green function of the sound field. According to (2), the thermoacoustic signal p(r 0 , t) detected at probe r 0 at time t = |r 0 − r|/c s is the integral of all points along the projection arc, as shown in Fig. 1. The electrical signal p 0 (r 0 , t) output by the ultrasonic probe is the result of convolution of the thermoacoustic signal p(r 0 , t) and the probe impulse response, which is expressed as below in (3): Refer to (3), * represents convolution, p(r 0 , t) means an ideal thermoacoustic signal without noise, i(t) is the time domain impulse response of the ultrasonic probe, and γ (t) represents random noise. This formula can be expressed in the frequency domain as below p 0 (r 0 , t): P 0 (r 0 , ω) , P(r 0 , ω) , I(ω) and γ (ω) are the frequency spectrum of p 0 (r 0 , t) , p(r 0 , t) , I(t) and γ (t) , respectively. The original acoustic signal at the ultrasound probe can be restored by the Wiener filter deconvolution method, namely: (4) P 0 (r 0 , ω) = P(r 0 , ω)I(ω) + γ (ω) www.nature.com/scientificreports/ Refer to (5), G(ω) is the Fourier transform spectrum of the Wiener filter, C = 1/(α|I(ω)|) is the regularization factor, and α is a mature factor to match the ratio of the power spectrum of the input signal and the noise signal in the entire period. FFT −1 represents the inverse fast Fourier transform. Using formula (6), the acoustic signal p 0 (r 0 , t) collected by the ultrasonic probe can be converted into the original acoustic signal p(r 0 , t).
According to the above original sound signal p(r 0 , t) , the thermal sound source H(r ′ ) can be reconstructed by the Filtered backprojection algorithm, as shown below in (7): Refer to (7), η = β/C p , and φ 0 is the circumferential angle of the location where the ultrasound probe is located.
The inverse problem of acoustic field is to reconstruct the thermoacoustic source from the collected thermoacoustic signals, while the inverse problem of electromagnetic field is to reconstruct the conductivity of the measured sample from the thermoacoustic source. Assuming that the target body is an acoustically homogeneous medium with electrical conductivity of σ (r ′ ) , and that the ultrasonic probe and the target body are both immersed in the insulating medium, with low-frequency electromagnetic excitation and electrical quasi-static approximation, the boundary conditions and governing equations that the target area meets are described as below in (8): where φ represents the electric scalar potential, 1 is the plane where the high voltage plate is located, 2 is the ground plate plane, and 3−6 are the other four boundaries of a cube imaging region except for the electrode. According to the governing equation given in (8), the electric field intensity in the target body can be calculated as below: The relationship between the thermal sound source term and the conductivity of the target body can be expressed as follows in (10): Refer to (10), the scalar potential φ is also a function of the target body's conductivity σ , therefore it is a nonlinear optimization problem to reconstruct the conductivity with known thermoacoustic source H(r ′ ). A two-dimensional ACTAI area is established for the imaging tomographic plane, which is discretized into rectangular cells with M rows and N columns, and converted into M × N-dimensional column vectors. Let n = M × N, and assume that the conductivity inside the j-th rectangular area is uniform, where j = 1,2,…,n. www.nature.com/scientificreports/ According to the ACTAI forward problem model, the thermal function in the imaging area can be calculated as H(r ′ ) = f ∈ R n×1 , and the thermal sound source reconstructed by filtered back projection is y ∈ R n×1 . The least square method is adopted to find the optimal solution, and the optimization problem is established as follows: Then the problem of reconstructing the conductivity of the target body σ is transformed into finding the optimal combination of σ to minimize the objective function S.
The conductivity reconstruction algorithm based on least squares solves the nonlinear problem through approximate linearization. However there are problems of low image reconstruction accuracy and poor antinoise performance. In addition, multiple iterations are often required in the calculation process, and the calculation cost is very high. Therefore, it is necessary to explore a more efficient, accurate and stable conductivity reconstruction method.
The Generative Adversarial Network draws on the idea of a two-player game, and is composed of two independent units: a Generator Network (GN) and a Discriminator Network (DN). As shown in Fig. 2, in the training process, the output of GN, that is, the generated sample, is usually used as the input of DN to complete the forward propagation process of the network. For the applied current thermoacoustic imaging problem, we use the thermoacoustic source image as the input of GN, and pair the conductivity image output by GN with the real conductivity image to generate samples, which together serve as the input of DN. In the back propagation process, GN and DN networks need to be separately trained to optimize the weight and bias of each network. When training DN, the pairing of the "fake image" generated by GN and the real conductivity image is marked as FALSE; while training GN, the pairing of the "fake image" generated by GN and the real conductivity image is marked as TRUE. In other words, the ultimate goal of training the GN network is to make the "fake images" generated by the GN deceive the DN network as much as possible, so that the DN network thinks this is a real pairing. While the ultimate goal of training the DN network is to make the DN network be able to identify as much as possible the "Fake image" generated by GN network. Through batch and repeated training of a large amount of data, the probability of FALSE and TRUE in the output result of DN is each 50%, that is, the DN network cannot judge the authenticity of the image generated by the GN network. At this time, the conductivity image generated by the GN network can "mix the false with the genuine" to satisfy the matching relationship with the real conductivity image.
The Generator Network of GAN is demonstrated in Fig. 3a. We use an improved U-shaped convolutional neural network 20 structure based on the residual network (ResNet). Firstly, the addition of the ResNet network avoids the degradation phenomenon that occurs with the increase in the number of network layers. Secondly, we also introduced a new convolutional layer and activation function Relu 21,22 in the front end of ResNet 23 , and the back end uses the Dropout layer to prevent over-fitting. The U-shaped symmetric structure and jump connection mode of the generating network can further improve the problem of gradient disappearance in the training process of the deep learning network. All convolutional and deconvolutional layers used in this paper have a convolution kernel size of 3 × 3, and the sliding step size is 2. The overall structure of the number of feature maps and the size of feature maps is shown in Fig. 3a.
The Discriminator Network of GAN is demonstrated in Fig. 3b. It adopts Patch-GAN (Markov Discriminator) structure and consists of 4 fully convolutional layers. First, the image is divided into 30 N × N patches, then each patch is judged, and finally the input is mapped to a matrix X. The average value of each element in X is the final output of the discriminator.
The loss function of a standard GAN network is: In addition to the loss function of standard GAN network, the loss function of pix2pix also introduces the loss function of L 1 : where x is the input image, y is the true image, G is the generator, and D is the discriminator. The loss of the standard GAN network is responsible for capturing the image, and L 1 is responsible for capturing low-frequency features such that the generated result is both true and clear. After the GAN training, the parameters of the Generator Network are fixed and saved for the reconstruction of the conductivity image. We can then use the simulation data to test the above-mentioned GAN network. The specific process is as follows: first, for each new conductivity sample, calculate the sound field distribution in the COMSOL software according to the multi-physics coupling positive problem model, and then the measurement signal is obtained after convolution probe characteristics to solve the positive problem; secondly, the original sound signal is obtained by wiener filtering deconvolution and the sound source distribution is reconstructed through filtered back projection; finally use the reconstructed sound source as the input of the generator in the GAN network, the output of which is the conductivity image to be solved.

Result
The simulation model shown in Fig. 4 is established, and the conductivity samples adopt the MNIST handwritten digit set that has been widely used in machine learning. The size of the target volume is 4cm × 4cm × 4cm , and its xoy interface is evenly divided into 64 × 64 grids. The handwritten digital image is assigned to the target plane as the conductivity distribution data and its size is consistent in the Z direction, in which the digital region conductivity is set as 1 and the conductivity of other regions is set as 0. 1000 handwritten digital images are selected as conductivity samples, and the electric field intensity, thermal sound source and sound field distribution are calculated by formulas (1), (8) and (10) using the finite element method. The 1000 sets of data obtained are randomly divided into two groups, with one set of 900 pieces targeted for network training, and the other set of 100 pieces for network testing.
Assuming that the volume expansion coefficient of the tissue is β = 3.8 × 10 −4 K −1 , the heat capacity is C p ≈ 3.94mJ/(gmK) , the propagation velocity of the thermoacoustic signal in the tissue and the ultrasonic coupling medium is 1404m/s , and the voltage applied by the electrode plate is U(t) = 10 4 × e −(t−b) 2 2c 2 , where b = 1 × 10 −6 and c = 1 × 10 −7 , then the thermal sound source generated by the target body at time 0.8µs is shown in Fig. 5.
If 72-channel ultrasound probes are placed around the target, then the acoustic signal they receive contains the frequency response characteristics of each probe. Figure 6 below demonstrates the comparison between the reconstructed sound source directly using FBP and the reconstructed sound source after preprocessing the sound signal with Wiener filter deconvolution. It can be found that the sound source reconstructed after deconvolution using Wiener filtering can more accurately restore the sound source distribution inside the target body. Due to the limitation of calculation amount, this paper only divides the tomographic plane into 64 × 64 grids. For future research, a tighter grid division can be used to get higher-quality results of sound source reconstruction.
Input the sound source in Fig. 6b as a test sample into the GAN network, and the reconstruction result is shown in Fig. 7.  Refer to (15), µ x is the average value of image x , µ y is the average value of image y , σ x is the variance of x , σ y is the variance of y , and σ xy is the covariance of x and y . c 1 = (k 1 L) 2 , c 2 = (k 2 L) 2 is a constant used to maintain stability, usually taking the following values: k 1 = 0.01 , k 2 = 0.03 , L = 255 . The value range of SSIM is 0-1, and when two images are exactly the same, the value of SSIM is 1.
When solving the PSNR value, it is necessary to first calculate the Mean Squared Error (MSE) of each target image, as indicated in (16): Then calculate the Peak Signal to Noise Ratio of the reconstructed image according to the MSE, as shown in (17): In the above formula, MAX I is the maximum pixel value in image I. The larger the PSNR value, the smaller the distortion. The conductivity reconstruction quality in Fig. 7 is shown in Table 1. Compared with Fig. 6, GAN not only accurately reconstructs the conductivity image of the target body, but also effectively suppresses the ripple noise generated during the reconstruction process. It is also worth mentioning that a well-trained GAN network can reconstruct a conductivity image in less than 1 s, which is much more efficient than traditional algorithms.
For the network trained by digital conductivity samples, we continue to use digital conductivity samples for testing. Although the test samples have never appeared in the training set, they are still lack of persuasion. Therefore, we use the sample which is quite different from the digital conductivity to test. In Fig. 8, the image of the thermoacoustic source as an input to the GAN network is shown. At the same time, in order to better verify the superiority of conductivity distribution reconstruction by GAN network, the conductivity reconstruction results of traditional method are also given for comparison. The traditional method used in this paper is the least  www.nature.com/scientificreports/ squares iterative algorithm 24 . By giving an initial conductivity value, the solution of conductivity is transformed into the solution of the objective function of conductivity. On the premise of meeting the calculation accuracy, the optimal solution is the conductivity. The conductivity reconstruction quality in Fig. 8 is shown in Table 2. It can be found that compared with the traditional method, the conductivity image reconstructed by GAN has less artifacts and higher resolution. In addition, GAN takes less time to reconstruct the conductivity image, which meets the real-time requirements of medical imaging.
We also tested the generalization ability of the GAN network: after training the network with handwritten digital samples, other samples that are significantly different from the training samples (handwritten digital samples) are used for imaging tests. Firstly, establish the phantom model of conductivity distribution as shown in Fig. 9. Set the major axis of outer ellipsoids to be 18 mm and its conductivity to be 0.1 S/m, the conductivity of two internal centered big ellipsoids to be 0.5 S/m and conductivity of other internal ellipsoids to be 1 S/m. According to the test sample, during the GAN network training, we make a complex upgrade to the traditional MNIST data set: any two different handwritten numerals are merged to form an irregular graph, and the conductivity is inconsistent. Although the background of phantom model is not uniform, the different conductivity values can be found in the training set. The complex training samples are shown in Fig. 10.
Secondly, place the ultrasound probe with a dominant frequency of 1 MHz at 0.02 m from the center of the phantom. When then probe position is at (− 0.02, 0) m, the time curve of the received acoustic signal is shown in Fig. 11a, and the time curve of acoustic signal after Wiener filter deconvolution is shown in Fig. 11b.  www.nature.com/scientificreports/ Thirdly, Fig. 12 shows two results of two different reconstruction methods with the FBP algorithm: one is with the detected signal directly, and the other is with signal preprocessed by Wiener filter deconvolution. The next step is to take the sound source image of Fig. 12b as the input of GAN, and reconstruct the conductivity image. The results are shown in Fig. 13 and Table 3. It can be found that although some noisy points are generated at the boundary of the reconstruction result, its main distribution is basically consistent with the real model.
Finally, in order to verify the anti-noise performance of the GAN network, white noise of different intensities is introduced into the thermal sound source image. Figure 14 displays the conductivity reconstruction results after applying noise with different Signal to Noise Ratio (SNR). It can be found that the error between the reconstructed conductivity and the conductivity of the original model increases as the SNR decreases. GAN network can get better reconstruction results when the SNR is 20 dB and 10 dB.
In order to ensure the feasibility of the method, we use cooled sodium chloride gel solution to prepare the target, and the concentration of sodium chloride determines the conductivity of the target. We conducted       Fig. 15, in which the conductivity was for the central portion of the mimic, at the periphery of the mimic, the outer circle diameter of the mimic was about 6 cm, and the inner circle diameter of the mimic was about 3.5 cm. The sources were reconstructed using filtered back projection method. The reconstruction results are shown in Fig. 16a. It can be seen that for the measured acoustic signal received by the ultrasound probe using the filtered back projection method, the conductivity information at the target body contour can be reconstructed more completely, but how its inner specific conductivity distribution can not be reconstructed well and the reconstruction results have more serious artifacts.
After obtaining the preliminary sound source reconstruction results, input them into the built GAN network for conductivity imaging. The results are shown in Fig. 16b. It can be seen that the designed network can still reconstruct the conductivity distribution of the target according to the measured ultrasonic data.   www.nature.com/scientificreports/

Discussion
Although the network we built is specifically trained on MNIST handwritten digit samples, satisfactory reconstruction results can still be obtained for other types of test samples. This indicates that in the Applied Current Thermoacoustic Imaging system, the GAN network has learned the mapping relationship between the input thermoacoustic source and the output conductivity. In other words, it can not only complete the image matching between input and output, but also has the ability to learn nonlinear problems. To conclude, this study focuses on Applied Current Thermoacoustic Imaging, mainly exploring the reconstruction method of inverse problem based on GAN network. Firstly, a physical model of electromagnetic and acoustic fields of Applied Current Thermoacoustic Imaging is established. Then the electric field in the model is calculated by finite element method, and the acoustic signal detected by the ultrasonic probe is simulated by the coupling relationship between electric field and acoustic field, and the thermoacoustic source distribution is reconstructed. In response to the problem of conductivity reconstruction, the research proposes an imaging method based on Generative Adversarial Network, and finally, with conductivity-based handwritten digit set training, it finally verifies that the network can learn the non-linear mapping relationship between the sound source and the conductivity, and thus realizes the deep learning reconstruction process of conductivity. Theoretical analysis and simulation results show that the method proposed in this paper can quickly and accurately reconstruct the conductivity image from the sound source image under the condition of low Signal to Noise Ratio, which verifies the applicability of the GAN network in the ACTAI problem.