A Mathematical Framework for Deep Learning in Elastic Source Imaging

An inverse elastic source problem with sparse measurements is of concern. A generic mathematical framework is proposed which incorporates a low- dimensional manifold regularization in the conventional source reconstruction algorithms thereby enhancing their performance with sparse datasets. It is rigorously established that the proposed framework is equivalent to the so-called \emph{deep convolutional framelet expansion} in machine learning literature for inverse problems. Apposite numerical examples are furnished to substantiate the efficacy of the proposed framework.


1.
Introduction. An abundance of real-world inverse problems, for instance in biomedical imaging, non-destructive testing, geological exploration, and sensing of seismic events, is concerned with the spatial and/or temporal support localization of sources generating wave fields in acoustic, electromagnetic, or elastic media (see, e.g., [8,10,20,32,37,41,45,47,50] and references therein). Numerous applicationspecific algorithms have been proposed in the recent past to procure solutions of diverse inverse source problems from time-series or time-harmonic measurements of the generated waves (see, e.g., [2,6,7,19,36,51,54,55,58,64,65]). The inverse elastic source problems are of particular interest in this paper due to their relevence in elastography [3,10,20,50]. Another potential application is the localization of the background noise source distribution of earth, which contains significant information about the regional geology, time-dependent crustal changes and earthquakes [21,32,34,37].
Most of the conventional algorithms are suited to continuous measurements, in other words, to experimental setups allowing to measure wave fields at each point inside a region of interest or on a substantial part of its boundary. In practice, this requires mechanical systems that furnish discrete data sampled on a very fine grid confirming to the Nyquist sampling rate. Unfortunately, this is not practically feasible due to mechanical, computational, and financial constraints. In this article, we are therefore interested in the problem of elastic source imaging with very sparse data, both in space and time, for which the image resolution furnished by the conventional algorithms degenerates. In order to explain the idea of the proposed framework, we will restrict ourselves to the time-reversal technique for elastic source localization conventional reconstruction frameworks. As it will be explained later on, the resulting algorithm can be extended to the deep convolutional framelet expansion in order to achieve an image resolution comparable to that furnished by the continuous/dense measurements [62].
The paper is organized as follows. The inverse elastic source problem, both in discrete and continuous settings, is introduced in section 2 and a brief review of the time-reversal algorithm is also provided. The mathematical foundations of the proposed deep learning approach are furnished in section 3. Section 4 is dedicated to the design and training of the deep neural network. A few numerical examples are furnished in section 5. The article ends with a brief summary in section 6.
2. Problem formulation. Let us first mathematically formulate the inverse elastic source problems with continuous and discrete measurements. Then, we will briefly review the time-reversal technique for elastic source imaging (with continuous data) as discussed in [5] in order to make the paper self-contained.
2.1. Inverse elastic source problem with continuous measurements. Let S : R d × R → R d be a compactly supported function. Then, the wave propagation in a linear isotropic elastic medium loaded in R d (d = 2, 3) is governed by the Lamé system, where u = (u 1 , · · · , u d ) : R d × R → R d is the elastic wave field generated by the source S, operator L λ,µ u = µ∆u + (λ + µ)∇(∇ · u) is the linear isotropic elasticity operator with Lamé parameters of the medium (λ, µ), and superscript indicates the transpose operation. Here, it is assumed for simplicity that the volume density of the medium is unit, i.e., λ, µ, and S are density normalized. Moreover, the source is punctual in time, i.e., S(x, t) = F(x)dδ 0 (t)/dt, where δ 0 denotes the Dirac mass at 0 and its derivative is defined in the sense of distributions.
Let Ω ⊂ R d be an open bounded smooth imaging domain with C 2 −boundary ∂Ω, compactly containing the spatial support of F(x) = (F 1 , · · · , F d ) ∈ R d , denoted by supp{F}, i.e., there exists a compact set Ω * ⊂ R d strictly contained in Ω such that supp{F} ⊂ Ω * ⊂ Ω. Then, the inverse elastic source problem with continuous measurement data is to recover F given the measurements It is precised that F and u can be decomposed in terms of irrotational components (or pressure components polarizing along the direction of propagation) and solenoidal components (or shear components polarizing orthogonal to the direction of propagation). In particular, in a two-dimensional (2D) frame-of-reference wherein xand y-axes are aligned with and orthogonal to the direction of propagation, respectively, the respective components of F are its pressure and shear components (see, e.g., Figure 1 for the imaging setup and source configuration).

2.2.
Inverse elastic source problem with discrete measurements. Most of the conventional algorithms require the measurement domain ∂Ω × (0, t max ) to be sampled at the Nyquist rate so that a numerical reconstruction of the spatial support is achieved at a high resolution. Specifically, the distance between consecutive receivers is taken to be less than half of the wavelength corresponding to the smallest frequency in the bandwidth and the temporal scanning is done at a fine rate so that the relative difference between consecutive scanning times is very small.
In practice, it is not feasible to place a large number of receivers at the boundary of the imaging domain and most often the measurements are available only at a few detectors (relative to the number of those required at the Nyquist sampling rate). As a result, one can not expect well-resolved reconstructed images from the conventional algorithms requiring continuous or dense measurements.
In the rest of this subsection, the mathematical formulation of the discrete inverse elastic source problem is provided. Towards this end, some notation is fixed upfront. For any sufficiently smooth function v : R → R, its temporal Fourier transform is defined byv where ω ∈ R is the temporal frequency. Similarly, the spatial Fourier transform of an arbitrary smooth function w : R d → R is defined bŷ with spatial frequency k ∈ R d . Let the function G ω be the Kupradze matrix of fundamental solutions associated to the time-harmonic elastic wave equation, i.e., where I d ∈ R d×d is the identity matrix. For later use, we decompose G into its shear and pressure parts as Here, H (1) 0 denotes the first-kind Hankel function of order zero, and c P = √ λ + 2µ and c S = √ µ are the pressure and shear wave speeds, respectively. Let y 1 , · · · , y M ∈ ∂Ω be the locations of M ∈ N point receivers measuring the time-series of the outgoing elastic wave u at instances 0 < t 1 < · · · < t N < t max for some N ∈ N. Then, the inverse elastic source problem with discrete data is to recover F given the discrete measurement set In this article, we are interested in the discrete inverse source problem with sparse data, i.e., when M and N are small relative to the Nyquist sampling rate. In order to facilitate the ensuing discussion, let us introduce the discrete measurement vector G ∈ R dM N by Here and throughout this investigation, notation [·] i indicates the i-th component of a vector and [·] ij indicates the ij-th component of a matrix. Thus, G n i , for 1 ≤ i ≤ d, denotes the vector formed by the i-th components of the waves recorded at points y 1 , · · · , y M at a fixed time instance t n . Let us also introduce the forward operator, D dis : . . .
Then, the inverse elastic source problem with discrete data is to recover F from the relationship 3. Time-reversal for elastic source imaging: A review. The idea of the time-reversal algorithm is based on a very simple observation that the wave operator in loss-less (non-attenuating) media is self-adjoint and that the corresponding Green's function possesses the reciprocity property [19]. In other words, the wave operator is invariant under time transformation t → −t and the positions of the sources and receivers can be swapped. Therefore, it is possible to theoretically revert a wave from the recording positions and different control times to the source locations and the initial time in chronology thereby converging to the source density. Practically, this is done by back-propagating the measured data, after transformation t → t max − t, through the adjoint waves v τ (for each time instance t = τ ) and adding the contributions v τ for all τ ∈ (0, t max ) after evaluating them at the final time t = t max . Precisely, the adjoint wave v τ , for each τ ∈ (0, t max ), is constructed as the solution to where δ ∂Ω is the surface Dirac mass on ∂Ω. Then, the time-reversal imaging function is defined by By the definition of the adjoint field v τ and the Green's theorem, Therefore, the time-reversal function can be explicitly expressed as The time-reversal function I TR in (2.5) is usually adopted to reconstruct the source distribution in an elastic medium. However, it does not provide a good reconstruction due to a non-linear coupling between the shear and pressure parts of the elastic field u at the boundary, especially when the sources are extended [34,8,5]. In fact, these components propagate at different wave-speeds and polarization directions, and cannot be separated at the surface of the imaging domain. If we simply backpropagate the measured data then the time-reversal operation mixes the components of the recovered support of the density F. Specifcally, it has been established in [5] that, by time reversing and back-propagating the elastic wave field signals as in (2.5), only a blurry image can be reconstructed together with an additive term introducing the coupling artifacts.
As a simple remedy for the coupling artifacts, a surgical procedure is proposed in [5] taking the leverage of a Helmholtz decomposition of I TR , (regarded as an initial guess). A weighted time-reversal imaging function (denoted by I WTR hereinafter) is constructed by separating the shear and pressure components of I TR as and then taking their weighted sum wherein the weights are respective wave speeds and the functions ψ ITR and φ ITR are obtained by solving a weak Neumann problem. Precisely, I WTR is defined by In fact, thanks to the Parseval's theorem and the fact that F is compactly supported inside Ω ⊂ R d , it can be established that for a large final control time t max with After tedious manipulations, using the elastic Helmholtz-Kirchhoff identities (see, e.g., [5, Proposition 2.5]), and assuming Ω to be a ball with radius R → +∞, one finds out that which comes from the integration of the time-dependent version of Eq. (2.1) between t = 0 − and t = 0 + , the following result holds (see, e.g., [5, Theorem 2.6]).
Let Ω be a ball in R d with large radius R. Let x ∈ Ω be sufficiently far from the boundary ∂Ω with respect to the wavelength and I WTR be defined by (2.6). Then, We conclude this section with the following remarks. Let D be the source-tomeasurement operator, defined in (2.2). Then, it is easy to infer from Theorem 2.1 that its inverse (or the measurement-to-source) operator is given by when imaging domain Ω is a ball with large radius R. However, there are a few technical limitations. Firstly, if Ω is not sufficiently large as compared to the characteristic size of the support of F, which in turn should be sufficiently localized at the center of the imaging domain (i.e., located far away from the boundary ∂Ω), one can only get an approximation of F which may not be very well-resolved. Moreover, I WTR may not be able to effectively rectify the coupling artifacts in that case as it has been observed for extended sources in [5]. Secondly, like most of the contemporary conventional techniques, time-reversal algorithm requires continuous measurements (or dense measurements at the Nyquist sampling rate). Therefore, as will be highlighted later on in the subsequent sections, very strong streaking artifacts appear when the time-reversal algorithm is applied with sparse measurements. In order to overcome these issues, a deep learning approach is discussed in the next section.
3. Deep learning approach for inverse elastic source problem. Let us consider the inverse elastic source problem with sparse measurements. Our aim is to recover F from the relationship (2.4). Unfortunately, (2.4) is not uniquely solvable due to sub-sampling. In fact, the null space, N (D dis ), of the forward operator D dis is non-empty, i.e., there exist non-zero functions, Moreover, the existence of the non-radiating parts of the source also makes the solution non-unique. This suggests that there are infinite many feasible solutions to the discrete problem (2.4). Hence, the application of the time-reversal algorithm requiring the availability of continuous or dense measurements results in strong imaging artifacts severely affecting the resolution of the reconstruction.
A typical way to avoid the non-uniqueness of the solution from sparse measurements is the use of regularization. Accordingly, many regularization techniques have been proposed over the past few decades. Among various penalties for regularization, here our discussion begins with a low-dimensional manifold constraint using a structured low-rank penalty [63], which is closely related to the deep learning approach proposed in this investigation.

Generic inversion formula under structured low-rank constraint.
Let {z q } Q q=1 ⊂ Ω, for some integer Q ∈ N, be a collection of finite number of sampling points of the region of interest Ω confirming to the Nyquist sampling rate. In this section, a (discrete) approximation of the density F is sought using piece-wise constants or splines ansatz where ϑ i (·, z q ) is the basis function for the i-th coordinate, associated with z q . Accordingly, the discretized source density to be sought is introduced by where superposed n and m indicate the dependence on n-th time instance for 1 ≤ n ≤ N and m-th boundary point y m for 1 ≤ m ≤ M , respectively. The subscripts 1 ≤ i, j ≤ d indicate that the (i, j)-th component of the Kupradze matrix is invoked and the index 1 ≤ q ≤ Q indicates that the basis function associated with the internal mesh point z q for the j-th coordinate is used. Accordingly, the sensing matrix Λ ∈ R dN M ×dQ is defined by Then, the discrete version of the relationship (2.4) is given by In order to facilitate the ensuing discussion, we define the wrap-around structured Hankel matrix associated to f i ∈ R Q , for i = 1, · · · , d, by where p i < Q is the so-called matrix-pencil size. As shown in [25,27,62,63] and reproduced in Appendix for self-containment, if the coordinate function [F] i corresponds to a smoothly varying perturbation or it has either edges or patterns, then the corresponding Fourier spectrumf (k) is mostly concentrated in a small number of coefficients. Thus, if f i is a discretization of [F] i at the Nyquist sampling rate, then according to the sampling theory of the signals with the finite rate of innovations (FRI) [53], there exists an annihilating filter whose convolution with the image f i vanishes. Furthermore, the annihilating filter size is determined by the sparsity level in the Fourier domain, so the associated Hankel structured matrix H pi (f i ) ∈ R Q×pi in the image domain is low-rank if the matrix-pencil size is chosen larger than the annihilating filter size. The interested readers are referred to Appendix or the references [25,27,62,63] for further details.
In the same way, it is expected that the block Hankel structured matrix of the discrete source vector f , constructed as is low-rank, where p = d i=1 p i . Let r i := rank(H pi (f i )) and r := d i=1 r i where rank(·) denotes the rank of a matrix. Then, a generic form of the low-rank Hankel structured constrained inverse problem can be formulated as It is clear that, for a feasible solution f = (f 1 , · · · , f d ) of the regularization problem (3.2), the Hankel structured matrix H pi (f i ), for i = 1, · · · , d, admits the singular value decomposition , v i ri ) ∈ R pi×ri denote the left and the right singular vector basis matrices, respectively, and Σ i = (Σ i kl ) ri k,l=1 ∈ R ri×ri refers to the diagonal matrix with singular values as elements. If there exist two pairs of matrices Φ i , Φ i ∈ R Q×S and Ψ i , Ψ i ∈ R pi×ri , for each i = 1, · · · , d and S ≥ Q, satisfying the conditions with the transformation C i : R Q → R S×ri given by which is often called the convolutional framelet coefficient [62]. In Eq. (3.4), where φ ik and ψ il denote the k-th and the l-th columns of Φ i and Ψ i , respectively. This implies that the Hankel matrix can be decomposed using the basis matrices B i kl .Here, the first condition in (3.3) is the so-called frame condition, R(V i ) denotes the range space of V i and P R(V i ) represents a projection onto R(V i ) [62]. In addition, the pair (Φ i , Φ i ) is non-local in the sense that these matrices interact with all the components of the vector f i . On the other hand, the pair (Ψ i , Ψ i ) is local since these matrices interact with only p i components of f i . Precisely, (3.4) is equivalent to the paired encoder-decoder convolution structure when it is un-lifted to the original signal space [62] C which is illustrated in Figure 2. The convolutions in (3.7) correspond to the multichannel convolutions (as used in standard CNN) with the associated filters, Here, the superposed prime over ψ ik ∈ R pi , for fixed i = 1, · · · , d, and k = 1, · · · , r i , indicates its flipped version, i.e., the indices of ψ ik are reversed [63]. Let us introduce the block matrices Then, thanks to conditions in (3.3), the pairs (Φ, Φ) and (Ψ, Ψ) satisfy the conditions Consequently, with the matrix transformation C : R dQ → R dS×r given by Let H i , for i = 1, · · · , d, refer to the space of signals admissible in the form (3.7), i.e., Then, the problem (3.2) can be converted to where the sub-matrices G i and Λ i are defined in (2.3) and (3.1), respectively.

Extension to Deep Neural Network.
One of the most important discoveries in [62] is that an encoder-decoder network architecture in the convolutional neural network (CNN) is emerged from Eqs. (3.4),(3.5), and (3.7). In particular, the non-local bases matrices Φ i and Φ i play the role of user-specified pooling and unpooling operations, respectively (see, subsection 3.3), whereas the local-bases Ψ i and Ψ i correspond to the encoder and decoder layer convolutional filters that have to be learned from the data [62].
Specifically, our goal is to learn (Ψ i , Ψ i ) in a data-driven fashion so that the optimization problem (3.8) (or equivalently (3.9)) can be simplified. Toward this, we first define Λ † (resp. Λ † i ) as a right pseudo-inverse of Λ (resp. Λ i ), i.e., ΛΛ † G = G (resp. Λ i Λ † i G i = G i ) for all G ∈ R dM N so that the cost in (3.8) (resp (3.9)) can be automatically minimized with the right pseudo-inverse solution. However, the solution leads to where f * denotes the true solution and f 0 ∈ N (Λ). Therefore, one looks for the matrices (Ψ i , Ψ i ) such that (3.10) can be engineered so that its output f belongs to R(Λ † ). This can be achieved by selecting the filters Ψ i 's that annihilate the null-space components f 0 i 's, i.e., In other words, the block filter Ψ should span the orthogonal complement of N (Λ). Therefore, the local bases learning problem becomes , denotes the training data-set composed of the input and ground-truth pairs. This is equivalent to saying that the proposed neural network is for learning the local basis from the training data assuming that the Hankel matrices associated with discrete source densities f i are of rank r i [62]. Still, the convolutional framelet expansion is linear, so we restricted the space so that the framelet coefficient matrices C i (f i ) are restricted to have positive elements only, i.e., the signal lives in the conic hull of the basis to enable part-by-part representation similar to non-negative matrix factorization (NMF) [38,39,40]: for i = 1, · · · , d. This positivity constraint can be implemented using the rectified linear unit (ReLU) [46] during training. Accordingly, the local basis learning problem (3.11) can equivalently be expressed as Here, the operator K i [Ψ i , Ψ i ] : R Q → R Q is defined analogously as in (3.10) but in terms of the mapping C i : R Q → R S×ri given by where denotes ReLU, i.e., for arbitrary matrix A ∈ R S×ri , we have [ (A)] kl ≥ 0, for all k and l.
The geometric implication of this representation is illustrated in Figure 3. Specifically, the original image f i is first lifted to higher dimensional space via Hankel matrix, H pi (f i ), which is then decomposed into positive (conic) combination using the matrix bases B kl i in (3.6). During this procedure, the outlier signals (black color) are placed outside of the conic hull of the bases, so that they can be removed during the decomposition. When this high conic decomposition procedure is observed in the original signal space, it becomes one level encoder-decoder neural network with ReLU. Therefore, an encoder-decoder network can be understood as a signal space manifestation of the conic decomposition of the signal being lifted to a higher-dimensional space. Fig. 3. Geometry of single layer encoder decoder network for denoising. The signal is first lifted into higher dimensional space, which is then decomposed into the positive combination of bases. During this procedure, the outlier signals (black color) are placed outside of the conic hull of the bases, so that they can be removed during the decomposition. When this high conic decomposition procedure is observed in the original signal space, it becomes one level encoder-decoder neural network with ReLU.
The idea can be further extended to the multi-layer deep neural network. Specially, suppose that the encoder and decoder convolution filter Ψ i and ζ i ( Ψ i ) can be represented in a cascaded convolution of small length filters: then the signal space is recursively defined as where, for all  = 1, · · · , J − 1 ∈ N, Here, the -th layer encoder and decoder filters, Ψ are the filter lengths, the number of input channels, and the number of output channels, respectively. This is equivalent to recursively applying high dimensional conic decomposition procedure to the next level convolutional framelet coefficients as illustrated in Figure 4(a). The resulting signal space manifestation is a deep neural network shown in Figure 4(b).

3.3.
Dual-Frame U-Net. As discussed before, the non-local bases Φ i and Φ i correspond to the generalized pooling and unpooling operations, which can be designed by the users for specific inverse problems. Here, the key requisite is the frame condition in (3.3), i.e., Φ i Φ i = I Q . As the artifacts of the time-reversal recovery from sparse measurements are distributed globally, a network architecture with large receptive fields is needed. Thus, in order to learn the optimal local basis from the minimization problem (3.11), we adopt the commonly used CNN architecture known as U-Net [48] and its deep convolutional framelets based variant, coined as Dual-Frame U-Net [22] (see Figure 5). These networks have pooling layers with down sampling, resulting exponentially large receptive fields.
As shown in [22], one of the main limitation of the standard U-Net is that it does not satisfy the frame condition in (3.3). Specifically, by considering both skipped connection and the pooling Φ i in Figure 5(a), the non-local basis for the standard U-Net is given by where Φ i,avg denotes an average pooling operator given by  Moreover, the unpooling layer in the standard U-Net is given by Φ i . Therefore, Consequently, the frame condition in (3.3) is not satisfied. As shown in [62], this results in the duplication of the low-frequency components, making the final results blurry.
In order to address the aforementioned problem, in the Dual-Frame U-Net [22], the dual frame is directly implemented. More specifically, the dual frame Φ i for the specific frame operator (3.13) is given by where the matrix inversion lemma and the orthogonality Φ i,avg Φ i,avg = I Q are invoked to arrive at the last equality. It was shown in [22] that the corresponding generalized unpooling operation is given by where B i denotes the skipped component. Equation (3.14) suggests a network structure for the Dual-Frame U-Net. More specifically, unlike the U-Net, the residual signal at the low resolution should be upsampled through the unpooling layer and subtracted from the by-pass signal to eliminate the duplicate contribution of the low-frequency components. This can be easily implemented using additional bypass connection for the low-resolution signal as shown in Figure 5(b). This simple fix allows the proposed network to satisfy the frame condition (3.3). The interested readers are suggested to consult [22] for further details.

Network design and training.
Let us now design the U-Net and Dual-Frame U-Net neural networks for the elastic source imaging based on the analysis performed in the previous section. For simplicity, consider a 2D case (i.e., d = 2) for the recovery of the x-component (i.e., [F] 1 ) of the unknown source. The y-component (i.e., [F] 2 ) of the source can be obtained in exactly the same fashion using the network architectures discussed above.

Description of the forward solver and time-reversal algorithm.
For numerical illustrations and generation of the training data, the region of interest Ω is considered as a unit disk centered at origin. Each solution of the elastic wave equation is computed over the box B = − β/2, β/2 2 so that Ω ⊂ B, i.e., (x, t) ∈ − β/2, β/2 2 × 0, t max with β = 4 and t max = 2. The temporal and spatial discretization steps are, respectively, chosen to be h t = 2 −6 t max and h x = 2 −7 β. The Lamé parameters are chosen in such a way that the pressure and the shear wave speeds in the medium are, respectively, c P = √ 3m.s −1 and c S = 1m.s −1 . The Lamé system u(x, 0) = 0 and ∂u ∂t (x, 0) = 0, is numerically solved over the box B with periodic boundary conditions. A splitting spectral Fourier approach [13] is used together with a perfectly matched layer (PML) technique [24] to simulate a free outgoing interface on ∂B. The weighted time-reversal function I WTR (x) also requires a Helmholtz decomposition algorithm. Since the support of the function I WTR (x) is included in Ω ⊂ B, a Neumann boundary condition is used on ∂B and a weak Neumann problem is solved in order to derive the Helmholtz decomposition. This decomposition is numerically obtained with a fast algorithm proposed in [59] based on a symmetry principle and a Fourier Helmholtz decomposition algorithm. The interested readers are suggested to consult [5, Sect. 2.2.1] for more details on the numerical algorithm.

Network architectures.
The original and Dual-Frame U-Nets consist of convolution layer, ReLU, and contracting path connection with concatenation (Figure 6). Specifically, each stage contains four sequential layers composed of convolution with 3 × 3 kernels and ReLU layers. Finally, the last stage has two sequential layers and the last layer contains only a single convolution layer with 1 × 1 kernel. The number of channels for each convolution layer is illustrated in Figure 6. Note that the number of channels is doubled after each max pooling layer. The differences between the original and Dual-Frame U-Nets are from additional residual paths illustrated in Figure 5.

Network training.
In order to deal with the imbalanced distribution of non-zero values in the label phantom images f * ( ) i and to prevent the proposed network to learning a trivial mapping (rendering all zero values), the non-zero values are weighted by multiplying a constant according to the ratio of the total number of voxel over the non-zero voxels. All the convolutional layers were preceded by appropriate zero-padding to preserve the size of the input. The mean squared error (MSE) is used as a loss function and the network is implemented using Keras library [18]. The weights for all the convolutional layers were initialized using Xavier initialization. The generated data is divided into 4000 training and 1000 validation data-sets. For training, the batch size of 64 and Adam optimizer [33] with the default parameters as mentioned in the original paper are used, i.e., the learning rate = 0.0001, β 1 = 0.9, and

5.
Numerical experiments and discussion. In this section, some numerical realizations of the proposed algorithm are presented for the resolution of the inverse elastic source problem and the performances of the proposed deep learning frameworks are debated. The examples of sparse targets with binary intensities and extended targets with variable intensities are discussed. The sparse targets are modeled by an elongated tubular shape and an ellipse. The extended targets are modeled by the Shepp-Logan phantom. The performance of the proposed framework is compared with the results rendered by the weighted time-reversal algorithm with sub-sampled sparse data and total variation (TV) based regularization approach applied on the low-resolution images provided by the time-reversal algorithm. The reconstructed images are compared under both clean and noisy measurement conditions with the TV-regulatization using fast iterative shrinkage threshholding algorithm (FISTA) of Beck and Teboulle [12]. For comparison, the peak-signal-to-noise ratio (PSNR) and the structural similarity index (SSIM) are used as metrics, where PSNR := 20 log 10 Here, M and N are the number of pixels in the rows and columns,f 1 and f * 1 are the reconstructed image and ground truth, µf 1 and µ f * 1 are the expectations, σ 2 f1 and σ 2 f * 1 are the variances, and σ 2 f1f * 1 is the covariance off 1 and f * 1 , respectively. Here, c 1 and c 2 are stabilization parameters and are chosen as c 1 = (0.01ξ) 2 and c 1 = (0.03ξ) 2 with ξ being the dynamic range of the pixel intensity. Figure 7 shows the variations of PSNR and SSIM values of the reconstructed test images using the standard and Dual-Frame U-Net. By increasing the number of recorders (N R ) or scanning rate (N K ), it is observed that the PSNR and SSIM values of the images show a monotonically increasing trend except for the SSIM value of the image from U-Net with N K = 128 (Figure 7(b)). On the other hand, the performance of the Dual-Frame U-Net always improved with more measurement data. In addition, the PSNR and SSIM values of the images from the Dual-Frame U-Net are always higher than the ones from the standard U-Net. This suggests that the Dual-Frame U-Net, which satisfies the frame condition, is a robust and predictable reconstruction scheme.  Figures 10 and 11 show the reconstruction results from the test data-set and Shepp-Logan phantom data. In particular, Figures 8 and 10 correspond to the noiseless measurements, while Figures 9 and 11 are from the noisy measurements. In a noisy condition, a white Gaussian noise with SNR= 5dB is added to the measurement and the images are reconstructed using the time-reversal algorithm. In both conditions, for the total variation algorithm, a regularization parameter γ = 0.02 is chosen without any other constraint and the FISTA algorithm is used. Here, we could not find a significant improvement in the quality of the images by varying the hyperparameter γ. These results are compared with the results by the neural networks which are trained on the images from the clean measurements only. Note that the network has seen neither the images from the noisy measurements nor the Shepp-Logan phantom during the training phase.

5.2.
Discussion. The denoising methods using the neural networks showed a superior performance over the total variation algorithm. Among those, the Dual-Frame U-Net showed the best results in both PSNR and SSIM. Though the standard U-Net recovered the overall shapes of the inclusions, it failed to find an accurate outfit and lacks the fine details of the inclusions (see Figures 12 and 13). For example, the recovered shapes of the ellipses using Dual-Frame U-Net in thin and sparse inclusions case have sharper ends than standard U-Net relative to that of the ground truth as highlighted in Figure 12. In addition, the standard U-Net failed to remove artifacts around the inclusions and bias in the background (Figures 12 and 13). On the other hand, the Dual-Frame U-Net recovered the oval shapes of the inclusions and their pixel values more accurately in both sparse and extended targets (pointed out by white arrows in Figures 12 and 13). These differences come from the overly emphasized low frequency components in the U-Net configuration that does not meet the frame 6. Conclusion. In this article, we showed that the problem of elastic source imaging with very sparse data, both in space and time, can be successfully dealt with our proposed deep learning framework. While the conventional denoising algorithm using TV regularization gives an unsatisfying reconstruction quality, deep learning approaches showed more robust reconstruction with better peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM). We showed that the network performance can be further improved by using the Dual-Frame U-Net architecture, which satisfies a frame condition. Appendix. To make this paper self-contained, here we briefly review the origin of the low-rank Hankel matrix as extensively studied in [62,63].
Note that many types of image patches have sparsely distributed Fourier spectra. For example, as shown in Figure 14(a), a smoothly varying patch usually has spectrum content in the low-frequency regions. For the case of an edge as shown in Figure 14(b), the spectral components are mostly localized along the k x -axis. Similar spectral domain sparsity can be observed in the texture patch shown in Figure 14(c), where the spectral components of the patch are distributed at the harmonics of the texture. In these cases, if we construct a Hankel matrix using the corresponding image patch, the resulting Hankel matrix is low-ranked [63].
In order to understand this intriguing relationship, consider a 1-D signal, whose spectrum in the Fourier domain is sparse and can be modeled as the sum of Dirac Suppose that we have a (r + 1)-length filter h that has the z-transform representation   Then, it is easy to see that [53] f h = 0. (6.4) Thus, the filter h annihilates the signal f and is accordingly referred to as the anni- hilating filter. Moreover, since Eq. (6.4) can be represented as the Hankel matrix H p (f ) is rank-deficient. In fact, the rank of the Hankel matrix can be explicitly determined by the size of the minimum-size annihilating filter [63]. Therefore, if the matrix pencil size p is chosen bigger than the minimum annihilating filter size, the Hankel matrix is low-ranked.