Muti-frequency extended sampling method for the inverse acoustic source problem

1.   Introduction

Inverse source problems (ISP) have many important applications including pollution source detection, optical tomography, antenna synthesis, sound source localization ^{[1,2,3,4,5,6]}. Inverse source problem with a single frequency data does not have a unique solution due to the existence of non-radiating sources and thus is very challenging ^[4,7,8]. It is possible to impose additional constraints on the source to obtain uniqueness. For example, one can seek the solution with a minimum $L^2$-norm ^[9]. The uniqueness of the inverse source source problem can be restored using multiple frequency data if the frequencies coincide with the Dirichlet eigenvalues of $D$ ^[10]. It is shown in ^[11] that the uniqueness can be obtained using a set of frequencies with an accumulation point. A logarithmic stability is also obtained. Recently, Isakov showed that multiple frequency data can improve the stability for the inverse source problem ^[12]. Algorithms using multiple frequency data for the inverse source problem have attracted many researchers ^{[11,12,13,14,15,16,17,18,19,20]}.

The extended sampling method (ESM) was proposed for the inverse scattering problem using the incident wave of a fixed direction or a fixed point source ^[21]. Similar to the linear sampling method ^[22,23] but with much less data, the ESM uses the solutions of ill-posed linear integral equations based on the far field operator. Since the ESM uses a small set of data and is effective to reconstruct the location and approximate the size of the target, it is attractive for the inverse source problems. In our pervious work ^[24], we generalize the ESM for the inverse source problem. A two-step strategy is used. The location of the source is reconstructed first. Then one looks for a disc centered at the found location such that the source is contained in the disc. An iterative algorithm is devised to optimize the radius of disc. As a consequence, the ESM only provides a rough approximation of the compact support of the source.

In this paper, assuming multiple frequency far field data are given, we propose the multi-frequency ESM (MESM) for the inverse acoustic source problem by introducing a new indicator function. The indicator function is constructed based on some linear ill-posed integral equations. We show that the indicator is small if the sampling point is inside the compact support of the source function and is large otherwise. In contrast to the ESM, the MESM lifts the restriction on the choice of the radius of the disc and relies mainly on the indicator function to reconstruct the source. The performance improves significantly as expected since more data are used. Numerical examples are presented to validate the effectiveness of the MESM.

The proposed method is fast and easy to implement. It provides satisfactory reconstruction of the compact support of the source. The result can serve as the initial guess for other methods, e.g., the optimization method, which reconstruct further information of the source. The rest of this paper is arranged as follows. In Section 2, we introduce the acoustic source problem at a fixed frequency and state the multiple frequency inverse source problem. Section 3 contains a brief introduction of the ESM for the inverse source problem using the single frequency data. In Section 4 we propose a multiple frequency ESM and define an imaging function based on the solutions of some linear ill-posed integral equations. Section 5 contains numerical examples of various sources. We draw some conclusions and discuss future work in Section 6.

2.   Direct and inverse source problems

We first introduce the acoustic source problem in the frequency domain and then state the multiple frequency inverse source problem. Let $D \subset \mathbb R^2$ be a bounded simply-connected domain and $k > 0$ be the wave number. Note that the wave number is proportional to the frequency and thus we do not distinguish them in the rest of the paper. Assume the source function is given by $f(x) \in L^2(D)$, the set of the square integrable functions defined on $D$. Let $H_{loc}^{1}(\mathbb{R}^{2})$ be the set of locally integrable functions whose derivatives are also locally integrable. The direct problem is to find a function $u\in H_{loc}^{1}(\mathbb{R}^{2})$ satisfying the Helmholtz equation

and the Sommerfeld radiation condition

It is well-known that the above problem has a solution $u$ (see, e.g., ^[2]) given by

where $\Phi(x, y, k) = -\frac{i}{4}H^{(1)}_0(k|x-y|)$ is the fundamental solution of the Helmholtz equation and $H^{(1)}_0$ is the Hankel function of the first kind and order zero ^[2]. Note that $k > 0$ is a fixed parameter and we emphasize that $u$ depends on $k$.

The well-posedness of the above direct problem defines a solution operator. To be precise, let $B$ be a disc that is large enough such that $\overline{D} \subset B$. One defines an operator ${\mathcal S}: L^2(B) \to L^2(\partial B)$ such that it maps $f(x) \in L^2(B)$ to $u|_{\partial B} \in L^2(\partial B)$, where $u$ is the solution to the above problem. It is known that ${\mathcal S}$ is a compact operator. Consequently, the inverse problem of reconstructing $f(x)$ from $u|_{\partial B}$ is ill-posed. Furthermore, the kernel of ${\mathcal S}$ is non-empty, which relates to the non-uniqueness of the inverse source problem at a fixed frequency. In fact, the functions in the kernel of ${\mathcal S}$ are the so-called non-radiating sources ^[4,8].

Since the solution $u$ satisfies the radiation condition, it precesses an asymptotic behavior as the location $x$ moves away from the source

where $\hat{x} = \frac{x}{|x|}\in \mathbb{S}: = \{\hat{x} \in \mathbb{R}^2| |\hat{x}| = 1\}$ and $u_{\infty}(\hat{x}, k)$ is called the far field pattern of $u(x, k)$. Due to the asymptotic behavior of the fundamental solution $\Phi(x, y, k)$ and (2.3), it holds that

The multiple frequency inverse source problem (MISP) considered in this paper is to approximate the compact support $D$ of $f(x)$ given the far field pattern $u_\infty(\hat{x}, k)$ for $\hat{x} \in \mathbb{S}$ and $k \in [k_1, k_2]$, $0 < k_1 < k_2$.

3.   Extended sampling method

In this section, we give a brief introduction of the ESM for the inverse source problem proposed in ^[24] and then analyze the property of the regularized solutions to an ill-posed integral equation. We start with the radiating solution to the scattering of a plane wave by a sound soft obstacle, which is the main ingredient for the ESM.

Let $B\subset \mathbb{R}^2$ be a sound soft obstacle. Denote by $u^{i}(x, d): = {\rm e}^{i kx\cdot d}, x, d \in \mathbb{R}^2$ the incident plane wave, where $|d| = 1$ is the direction. The scattering problem is to find the total field $u = u^i+u^s$, where $u^s$ is the scattered field, such that

There exists a unique solution $u$ to the above problem ^[23]. Furthermore, similar to (2.4), the scattered field $u^s$ has the asymptotic expansion

uniformly in all directions $\hat{x} = x/|x|$.

When the obstacle is a disc centered at the origin with radius $R$, denoted by $B^R$, the far-field pattern for the scattered field can be written down in a closed form (Chp. 3 of ^[23]):

where $J_n$ is the Bessel function, $H_{n}^{(1)}$ is the Hankel function of the first kind of order $n$, $\theta = \angle (\hat{x}, d)$, the angle between $\hat{x}$ and $d$. If the center of the disc is shifted to $z$, i.e.,

the far field pattern is simply

Let $z \in \mathbb R^2$ and we define a far field operator $\mathcal{F}_z: L^2(\mathbb{S}) \to L^2(\mathbb{S})$ using the far field pattern of the disc $B_z$ due to the incident plane waves of all directions $d \in \mathbb{S}$:

Now we introduce the ESM for the inverse source problem for a fixed frequency proposed in ^[24]. Note that there is no incident wave any more, i.e., no dependence on $d$ for the radiating field. However, the kernel of the far field equation is the same, which depends on $d$. Let $u_\infty(\hat{x}, k)$ be the far field pattern for the inverse source problem defined in (2.5). The far field equation for the inverse source problem is

The above equation is ill-posed. One usually seek some regularized solution $g_z^\epsilon(d, k)$ for (3.5). The ESM uses $g_z^\epsilon(d, k)$ to decide the location and rough size of the source.

We check the property of the regularized solution $g_z^\epsilon(d, k)$. Let $B_z^R$ be a ball centered at $z$ with a large enough radius $R$. We shall consider two cases: (a) $\overline{D} \subset B_z^R$ and (b) $\overline{D}\cap \overline{B_z^R} = \emptyset$.

(a) $\overline{D} \subset B_z^R$. One can find a domain $D_0$ such that $D \subset D_0 \subset B_z^R$. Define

Then there exists a positive number $\sigma$ small enough such that $0 < d_2 < d_1 < \sigma$.

Consider an auxiliary scattering problem for a sound soft obstacle of finding the scattered field $u^s$ such that

where $u$ is the solution to the source problem (2.1) and (2.2). The above scattering problem has a solution $u^s$ such that $u^s$ coincides with $u$ in $\mathbb R^2 \setminus \overline{D}_0$ ^[23]. Define the Herglotz wave function operator $\mathcal{H}: L^2(\mathbb{S})\rightarrow H^{{1}/{2}}(\partial B_z)$:

We call $k^2$ a Dirichlet eigenvalue for the negative Laplacian for a domain $\Omega$ if there exists a nontrivial function $w$ satisfying

If the wavenumber $k^2$ is not a Dirichlet eigenvalue for the negative Laplacian for $B_z^R$, $\mathcal{H}$ is injective and has a dense range.

Define the near-to-far field operator ${\mathcal N}: H^{1/2}(\partial B_z^R)\rightarrow L^2(\mathbb{S})$, which maps the boundary value of $u\in H_{loc}^1(\mathbb{R}^2\setminus\overline{B_z^R})$ to its far field pattern $u_\infty$. Note that $u|_{\partial B_z^R}$ uniquely determines the far field pattern of $u$. Furthermore, the operator ${\mathcal N}$ is bounded, injective and has a dense range ^[23].

For any $\varepsilon > 0$, there exists a kernel function $g_z^\varepsilon$ such that

From the definition of the Herglotz wave function (3.7) and the far field pattern of the scattered field of a disc due to a plane incident wave (3.3), one has that

Due to the fact that $D_0\subset B_z^R$, it holds that

where $u^s$ is the solution to (3.6). Subtracting (3.9) from (3.10), and using (3.8), one obtains that

Taking $\varepsilon\rightarrow 0$ in (3.8), the Herglotz wave function $v_{g_z^\varepsilon}(x): = \int_{\mathbb{S}}e^{ikx\cdot d}g_z^\varepsilon (d)ds(d)$ converges to the unique solution $w\in H^1(B_z^R)$ of the following problem

(b) $\overline{D}\cap \overline{B_z^R} = \emptyset$. There exists $g_z^\varepsilon\in L^2(\mathbb{S})$ such that

Assume that there exists a sequence $v_{g_z^{\varepsilon_n}}$ such that $\|v_{g_z^{\varepsilon_n}}\|_{H^1(B_z^R)}$ remains bounded as $\varepsilon_n \to 0, n \to \infty$. Without loss of generality, we assume that $v_{g_z^{\varepsilon_n}}$ converges to $v_{g_z}\in H^1(B_z^R)$ weakly as $n\rightarrow \infty$, where $v_{g_z}(x) = \int_{\mathbb{S}}e^{ikx\cdot d}g_z (d)ds(d), x\in B_z^R.$ Let $v^s\in H_{loc}^1(\mathbb{R}^2\setminus \overline{B_z^R})$ be the unique solution of the following problem

The corresponding far field pattern is $v_\infty = \int_\mathbb{S} u^{B_z^R}_{\infty}(x, d)g_z(d)ds(d).$ From (3.12), as $\varepsilon\rightarrow 0$, we have $\int_\mathbb{S} u^{B_z^R}_{\infty}(x, d)g_z(d)ds(d) = u_\infty.$ Consequently, $v_\infty = u_\infty$. Then by Rellich's lemma (see Lemma 2.12 in ^[23]), the scattered fields coincide in $\mathbb{R}^2\setminus (\overline{B_z^R}\cup \overline{D}_0)$ and

We have that $v^s$ has an extension into $\mathbb{R}^2\setminus B_z$ and $u^s$ has an extension into $\mathbb{R}^2\setminus D_0$. Since $D_0\cap B_z^R = \emptyset$, $\tilde{u}^s$ can be extended from $\mathbb{R}^2\setminus (\overline{B_z^R}\cup \overline{D_0})$ into all of $\mathbb{R}^2$, that is, $\tilde{u}^s$ is an entire solution to the Helmholtz equation. Since $\tilde{u}^s$ also satisfies the radiation condition, it must vanish identically in all of $\mathbb{R}^2$. Since $\tilde{u}^s$ does not vanish identically, it leads to a contradiction and hence there does not exist a sequence $v_{g_z^{\varepsilon_n}}$ such that $\|v_{g_z^{\varepsilon_n}}\|_{H^1(B_z^R)}$ remains bounded as $\varepsilon_n \to 0, n \to \infty$.

The above derivation proves the following properties of the regularized solution for the integral equation (3.5).

Lemma 3.1. Let $B_z^R$ be a sound-soft disc centered at $z$ with a large enough radius $R$ such that $\overline{D} \subset B_z^R$. Furthermore, assume that $k^2$ is not a Dirichlet eigenvalue for $B_z^R$. Let $u_\infty(\hat{x}, k)$ be the far field pattern for the acoustic source problem defined in (2.5). Then the following properties hold.

(a) If $\overline{D}\subset B_z^R$, for a given $\varepsilon > 0$, there exists a function $g_z^\varepsilon\in L^2(\mathbb{S})$ to (3.12) and the Herglotz wave function $v_{g_z^\varepsilon}(x): = \int_{\mathbb{S}}e^{ikx\cdot d}g_z^\varepsilon (d)ds(d)$ converges to the solution $w\in H^1(B_z^R)$ of the Helmholtz equation with $w = u$ on $\partial B_z^R$ as $\varepsilon\rightarrow 0$.

(b) If $\overline{D}\cap \overline{B_z^R} = \emptyset$, every $g_z^\varepsilon\in L^2(\mathbb{S})$ that satisfies (3.12) for a given $\varepsilon > 0$ is such that $\lim_{\varepsilon\rightarrow 0}\|v_{g_z^\varepsilon}\|_{H^1(B_z)} = \infty.$

4.   Multiple frequency extended sampling method

The above method can reconstruct the location of the source. However, the size of the source is not always satisfactory. In fact, it is challenging to obtain a satisfactory reconstruction given a single frequency data. Recent studies have focused on the multiple frequency inverse source problems ^[10,11,15]. This motivates us to modify the ESM to process multiple frequency data to obtain a better reconstruction of the acoustic source. Numerical examples show that the method also works for frequency dependent sources, which makes it more attractive than some existing methods.

Recall the MISP and assume the far field pattern $u_\infty(\hat{x}, k)$ for $\hat{x} \in \mathbb{S}$ and $k \in [k_1, k_2]$, $0 < k_1 < k_2$ are given. In practice, one usually knows the far field pattern for a discrete set of frequencies, i.e., $u_\infty(\hat{x}, k_m), k_i \in [k_1, k_2], m = 1, \ldots, M$. Let $g_z^\epsilon(d, k_m)$ be a regularized solution to

Note that we write $g_z^\epsilon(d, k_m)$ instead of $g_z^\epsilon(d)$ to emphasize the dependence of $g_z^\epsilon$ on $k_m$ explicitly.

Based on the property of the regularized solution of the linear ill-posed equation, we define an image function

at a point $z \in S$, where $S$ is a region known as a priori and containing $D$. In general, $S$ is much larger than $D$. The following theorem shows the property of $I_z$.

Theorem 4.1. Assume that $R$ is large enough. If $z \in D$ such that $D \subset B^R_z$, a disc with radius $R$ and centered at $z$, then $I_z$ is bounded as $\epsilon \to 0$. If $z$ is such that $D\cap B^R_z = \emptyset$ and there exists a sequence of regularized solutions $g_z^\epsilon(d, k_m)$ satisfying (4.1), then $I_z \to \infty$ as $\epsilon \to 0$.

Proof. Let $B^R_z$ be such that $z \in D$ and $D \subset B^R_z$. In addition, let $k_m$ be fixed and $g_z^\epsilon(d, k_m)$ be a regularized solution to (4.1). According to (a) of Lemma 3.1, there exists a regularized solution in the form of the Herglotz wave function $v_{g_z^\varepsilon}(x, k_m): = \int_{\mathbb{S}}e^{ikx\cdot d}g_z^\varepsilon (d, k_m)ds(d), x\in B_z^R$. Furthermore, $v_{g_z^\varepsilon}(x, k_m)$ converges as $\epsilon \to 0$. Hence $\|g_z^\varepsilon (d, k_m)\|_{L^2(\mathbb S)}$ is bounded as $\epsilon \to 0$. Since the sum in (4.2) is over finitely many frequencies, we conclude that $I_z$ is bounded as $\epsilon \to 0$.

Next let $z$ be such that $D\cap B_z^R = \emptyset$. Again, for each fixed $k_m$, since $g_z^\epsilon(d, k_m)$ satisfies (4.1), it follows that $\|v_{g_z^\varepsilon}\|_{H^1(B_z^R)} \to \infty$ as $\epsilon \to 0$ again due to (b) of Lemma 3.1. Hence $\|g_z^\varepsilon (d, k_m)\|_{L^2(\mathbb S)} \to \infty$ as $\epsilon \to 0$. Consequently, $I_z \to \infty$ as $\epsilon \to 0$.

The above theorem justifies that $I_z$ can be viewed as an image function for the compact support $D$ of the source function. In particular, we expect $I_z$ is relatively small for $z \in D$ and relative large for $z \notin D$.

We devote the rest of this section to the implementation of the multiple frequency ESM (MESM) for the inverse source problem. The synthetic data of the direct source problems are computed using the integral formula (2.5). We first generate a triangular mesh ${\mathcal T}$ for $D$ with mesh size $h$. The far-field pattern is approximated by

where $T \in {\mathcal T}$ is a triangle, $y_T$ is the center of $T$, and $|T|$ is the area of $T$. For observation/measurement directions, the interval $[0, 2\pi]$ is uniformly divided into $N$ intervals and let $\theta_n = 2n\pi/N, n = 0, \ldots, N$. The discrete wave numbers are uniformly distributed on $[k_1, k_2]$, namely,

For all examples, we compute the far field data using (4.3)

Then we add 3% uniformly distributed random noise to the data for numerical simulations.

For each $k_i$, we discretize (3.5) and end up with a linear ill-posed system, which is solved using Tikhonov Regularization with the regularization parameter $\alpha = 10^{-2}$. We comment that the value does not affect the performance of the method significantly.

We denote the set of all sampling points by $S$ as well. For a fixed wavenumber $k_m$, we let $z \in S$ and, by (3.3), compute

Using the trapezoidal rule for (3.5), one obtains the following ill-posed linear system

where

Then we apply the Tikhonov regularization to compute a regularized solution ${\boldsymbol g}^\alpha_{k_m, z}$ of (4.4). The normalized discrete imaging function corresponding to (4.2) is

We summarize the algorithm of the MESM as follows. Assuming the sampling domain $S$ is known, the main program takes the noisy measurement data and outputs the image function $I_z, z \in S$. The subroutine $TSESM$ essentially computes the regularized solution for (4.4) given the wavenumber and radius of the disc. The subroutine $Radius$ chooses a suitable radius $R_m$ for the given data and wavenumber $k_m$.

$I_z =$MESM

  given $u_\infty(\hat{x}_n, k_m)$ and $S$.

  choose a regularization parameter $\alpha$.

  initialize the image function $I_z$.

${\rm{i}}$. for each $z \in S$,

  for $m = 1, \ldots, M$,

    * $R_m = Radius(u_\infty(\hat{x}_n, k_m))$.

    * $I_{z, m} = TSESM(u_\infty(\hat{x}_n, k_m), R_m, \alpha)$.

    * $I_z = I_z + I_{z, m}$.

${\rm{ii.}}$ return $I_z$.

Subroutine $R_m = Radius(u_\infty(\hat{x}_n, k_m))$

$1)$ choose a large enough radius $R$ for $B_z^R$.

$2)$ $I_{z, R}^m = TSESM(u_\infty(\hat{x}_n, k_m), R)$.

$3)$ take the global minimum point $z^*\in S$ for $\|g^\alpha_{z}\|_{l^2}$ as the location for $D$.

$4)$ generate a decreasing sequence of $R_j$ such that $R > R_1 > \ldots > R_J$.

$5)$ for $j = 1, 2, \cdots, J$

  ${\rm{a}}$ generate $S_j$ with the distance between sampling points being roughly $R_j$.

  ${\rm{b}}$ $I_{z, R}^m = TSESM(u_\infty(\hat{x}_n, k_m), R_j)$.

  ${\rm{c}}$ find the minimum point $z^*_j\in S_j$ and set $D_j: = B_{z_j}^{R_j}$. If $z_{j} \not\in D_{j-1}$, stop.

6) return $R_m = R_{j-1}$.

Subroutine $I_{z, R}^m = TSESM(u_\infty(\hat{x}_n, k_m), R, \alpha)$

  for each sampling point $z \in S$,

$1)$ compute $A^{k_m}$, ${\boldsymbol g}_{k_m}$, and ${\boldsymbol u}_{k_m}$.

$2)$ use the Tikhonov regularization to compute an approximate $g^\alpha_{z, k_m}$ to (4.4).

$3)$ $I_{z, R}^m = \|g^\alpha_{z, k_m}\|_{l^2}$.

5.   Numerical examples

We follow the discussion above to generate far field data. The algorithm is implemented using Matlab. We choose four different domains for $D$: a rectangle, a circle, an L-shaped domain, and a kite. The rectangle and the circle are defined as

The L-shaped domain is given by

The boundary of the kite is given by

In Figure 1, the domains and the associated meshes are plotted.

Four functions are used:

Note that the last two functions depend on the wavenumber $k$, which are not covered by the theory in Section 4. Note that some of the existing method do not consider the source functions depending on $k$ (see, for example, ^[10]). However, we shall see that the performance of the MESM is similar indicating the effectiveness of the proposed method for wavenumber depending sources.

For all examples, we set $k_1 = 1, k_2 = 6$ and uniformly divide the interval $[k_1, k_2]$ into $M = 51$ wavenumbers. We use $N = 32$ observation/measurement locations for the far field patterns. The regularization parameter is set to be $\alpha = 10^{-2}$. We notice that different regularization parameters do not affect the performance of MESM significantly. Assume that the search region is $S: = [-2, 2] \times [-2, 2]$. Then one uniformly divides $S$ into $40\times 40$ squares and ends up with $41 \times 41$ sampling points $z$'s.

To compute the far field data, we generate triangular meshes for these domains with mesh size $h \approx 0.01$. Note that since the mesh size is rather small, the numerical error using (4.3) can be ignored. Then 3% uniformly distributed noises are added to the far field data.The image functions are computed following the procedure described in Section 4.

We first consider the domain $D_1$. In Figure 2, we plot the image functions $I_z$ in the sampling domain for different source functions. The boundary of the exact support is the solid line. For different source functions, the image functions behave similarly with slight differences. It is clearly that the values of $I_z$ is small around the source and becomes larger when $z$ moves away from the source.

Next we consider the domain $D_2$, a disc with radius $0.6$. We plot the image functions $I_z$ in Figure 3. For all four source functions, the results are satisfactory.

In Figure 4, the results for $D_3$, an L-shped domain, are shown. The reconstructions again correctly provide the location and size of the source. Note that $D_3$ is a non convex domain. MESM cannot provide the details of the boundary, in particular, conners of the domain.

In Figure 5, the contour plots for $D_4$, a kite, are shown. Again $D_4$ is a non convex domain. MESM does not resolve the shape of the support well. Note that the source functions $f_3$ and $f_4$ are wave number dependent. The results indicate that MESM is effective for frequency-dependent sources.

Note that all the above examples use far field pattern computed by (4.3) with 3% uniformly distributed random noises. Sampling type methods are robust with respect to noises in general. It is also true for the proposed MESM, which is validated by the following example. Let the domain be the L-shaped domain and the source be $f_1$. In Figure 6, we show the contour plots using the far field pattern with 5% and 10% uniformly distributed random noises, respectively. It can be seem that the reconstructions are similar to the top-left picture in Figure 4.

6.   Conclusions

We consider the inverse source problem using multiple frequency data and propose a new version of the extend sample method to reconstruct the compact support of the source. The behavior of the image function is justified theoretically. Various numerical examples are presented to demonstrate the performance of the new method. Compared with the results using ESM of single frequency data in ^[24], MESM provides much better reconstructions than ESM. However, it is noted that MESM cannot provide the details of the boundary the domain. Nonetheless, the method is very fast and easy to implement.

There have been significant interests in the reconstruction of the actual source function $f$. The proposed MESM can provide a good estimate of the compact support of $f$ and thus can be integrated into other methods, e.g., the optimization method and Bayesian inversion, to reconstruct more information of the source. The current image function treats all frequency data in the same way. An interesting question is that shall we assign different weights to different frequencies? Finally, it is worthwhile to investigate how to decide a cutoff value for the image function if one requires a single contour reconstruction for the compact support. In fact, how to decide the cutoff value is pertinent to almost all sampling type methods such as the linear sample method and the direct sampling method and thus needs further investigations.

Conflict of interest

The authors declare there is no conflicts of interest.