Multifrequency Linear Sampling Method on Experimental Datasets

We investigate the use of the linear sampling method (LSM) for determining the shape of a scatterer from multifrequency experimental data. We study three multifrequency indicators for two 2-D datasets available online: one is provided by the Institut Fresnel, and another by the Electromagnetic Imaging Laboratory of the University of Manitoba. We show that the multifrequency LSM works exceptionally well on the 2-D Fresnel database, and also acceptably well on the Manitoba one. In particular, a new multifrequency indicator is tested, and data completion for the Fresnel dataset is studied. We also test an adaptive technique to cut down on the number of evaluations of the indicator function for well-resolved scatterers.


I. INTRODUCTION
T HE linear sampling method (LSM) introduced in [1] is a technique for solving inverse scattering problems that uses multistatic time-harmonic scattering data to reconstruct the boundary of penetrable dielectric or impenetrable metallic objects.By solving a sequence of linear ill-posed problems associated with sampling points in a region containing the scatterer, an indicator function is computed that is expected to be small outside the scatterer and larger inside.A key feature is that the coefficient matrix for the linear system does not depend on the sampling point, allowing for a fast solution algorithm.Since the original article [1], the LSM has been studied for many different equations (e.g., Helmholtz, Maxwell's, and linear elasticity) in different geometries (for example, waveguides [2]), and extended to time domain electromagnetism [3].Of particular relevance to this article, numerical aspects of this problem using regularization via Morozov's principle were studied in [4].For an overview of the LSM and related methods in the context of acoustic or electromagnetic scattering, we refer the reader to [5] and [6].
An issue related to the practical use of the LSM in nondestructive testing is the question of what quality of reconstruction can be obtained from electromagnetic scattering data available from experimental antenna-transmitting systems.There has already been work on testing the LSM with real data, including the early study [7] using single-frequency data, as well as [8] concerning the use of the LSM on data from the Institut Fresnel that we shall discuss shortly.These studies did not consider combining results for multifrequency data, although multifrequency data are often available.
The first theoretical work to consider the LSM for multifrequency data is [9], where two different ways (termed parallel or serial indicators) were proposed to combine the indicator functions from data at each frequency.We discuss these two multifrequency indicators in more detail in Section IV-B.
More recently, to handle multifrequency data and enhance reconstructions, two modifications of the LSM have been suggested in [10] and [11].The first [10] suggests a method for combining multifrequency data by constraining the results for nearby frequencies using changes in path length, while the second modification [11] uses extra constraints related to the boundary conditions to help sharpen the image.These modified LSM variants show a great deal of promise for data that are available only in a small aperture or from sparse measurements.However, both algorithms complicate the LSM implementation: the multifrequency modification in [10] results in a different coefficient matrix for the linear system at each sampling point, requiring much greater computational effort than the standard LSM.We return to this point in Section IV-B.The second modification in [11] uses boundary information and results in a nonlinear and possibly nonconvex optimization problem.Because these issues will likely slow down the LSM, we instead propose an alternative way to combine single-frequency data (see Section IV-B) that is aimed at preserving the speed of the original LSM; this is our so-called product indicator.To further enhance the speed of the reconstruction, we also test an adaptive approach similar to that in [12] on the experimental data.This is described in more detail in Section IV-C.
The articles most relevant to ours are [13] and [14].In these articles, a special metamaterial-based antenna system is used to probe metallic cylinders in a 2-D configuration and spheres in 3-D; for the latter, the antenna system has limited aperture.
In addition, the authors suggest a new indicator based on the geometric mean of the individual indicators at each frequency [14].We will study the use of a closely related indicator for penetrable dielectric and impenetrable metallic targets.
We shall apply the LSM to multifrequency experimental data from the two repositories that exhibit contrasting issues for the LSM.Both repositories provide data for elongated targets that can be considered as 2-D scattering problems, but each dataset has different challenges.
The first dataset is from the Electromagnetic Imaging Laboratory of the University of Manitoba [15], and is discussed in more detail in Section III-A.Measurements are available for several targets including pipes and more exotic phantoms and combinations of materials.From the point of view of LSM application, for the Manitoba dataset the measurement array is sufficiently dense and surrounds the target so that the standard LSM can be applied (after data completion only in the backscattering direction).Moreover, multifrequency data can be combined using the serial or parallel approach in [9] since the reconstructions for each frequency alone carry information about the respective targets (although we shall show that our new indicator provides sharper reconstructions).However, the mutual interaction between nearby antennas which we do not model results in a noisy dataset that challenges the LSM.
The second dataset we use to test the LSM is from the Institut Fresnel [16] and is described in more detail in Section III-B.In this setup, the interaction of nearby antennas is removed by having only two of them present: the transmitter which is fixed, and the receiver which is allowed to rotate around the scatterer.To measure different angles of incidence the scatterer is then rotated.This, however, allows for error in the positioning of the scatterer that is manifested as a solid rotation, as we will see in Section V-B.Due to physical restrictions, both antennas cannot be closer than 60 • , hence there is a lack of scattering data for angles surrounding the transmitting direction.
As we shall see a simple completion strategy first suggested in [8] and [17] works well in this case.We shall offer an explanation of why this is the case in Section IV-A.
For both datasets, it is desirable to improve on the parallel and serial indicators of [9], particularly in the case of the nonconvex U-shaped scatterer probed in the Fresnel dataset, in which the reconstruction tends to resemble the convex hull of the U.This is the motivation for the product indicator in this article.
The main contributions of this article are: 1) tests of the basic LSM on multifrequency real data; 2) tests of the product indicator for the multifrequency problem; 3) a partial explanation of why the completion strategy of [17] works; and 4) tests of an adaptive algorithm on real data.
The outline of the article is as follows.In Section II, we describe the basic single-frequency LSM in the near-field for the underlying scattering model.Next, in Section III, we introduce the experimental settings used by the Electromagnetic Imaging Laboratory of the University of Manitoba, and by the Institut Fresnel, and explain the challenges of each dataset for the LSM.In Section IV, we describe the usual implementation of the LSM for an individual frequency.We also derive a justification of the simple data completion proposed in [8] for the Fresnel database, propose three different multifrequency indicators, and describe an adaptive approach.In Section V, we provide LSM reconstructions for some of the objects inspected in the Manitoba and Fresnel experimental setups.We conclude in Section VI with some observations about our results and a few comments about possible future work.

II. LSM IN THE NEAR-FIELD
We now describe the LSM in the ideal case of infinitely many measurements and sources.In this work, we consider experimental setups in which an elongated dielectric or metallic scatterer is irradiated by microwaves from some linearly polarized antennas oriented in the same direction as the aforementioned scatterer.This allows modeling the scattering problem using the 2-D scalar Helmholtz equation.More precisely, we assume that each transmitter can be modeled as a 2-D point source giving rise to an incident field.If the point source is located at x 0 , for the wavenumber κ > 0 the incident field is where k is the fundamental solution of the Helmholtz equation, H (1)  0 is the first kind Hankel function of zero order, and C κ is a complex amplitude.This complex amplitude does not play a significant role in the single-frequency LSM, however, we will comment on this point more after we have defined the multifrequency indicators.
When the target is metallic, the corresponding total and scattered fields, u(•) ≡ u(• ; x 0 , κ) and u sc (• ; where D denotes the cross section of the scatterer.If the scatterer is a penetrable dielectric, then these fields solve instead Here, the superscripts + and − denote the limit values when approaching the boundary ∂ D from its outside and inside, respectively, and ε r is the relative electrical permittivity of the scatterer, which we assume to be equal to unity outside D. As we will describe in Section III, both experimental setups give rise to near-field measurements. More precisely, we suppose that the scattered field is measured on a simple closed curve R enclosing the target D. We also suppose that the source transmitters are located on a simple closed curve T which contains D (possibly R = T ).Thus, in the ideal case, we know u sc (x; y, κ) for all x ∈ R and y ∈ T , and for several choices of the wavenumber κ.From these data, we want to determine the boundary ∂ D.
For the previously mentioned full continuous data, the LSM uses the near-field operator where L 2 ( ) is the vector space of complex-valued square integrable functions on with norm ∥•∥ L 2 ( ) (where = R or = T ).
The LSM then proceeds by solving the following ill-posed linear near field equation for each sampling point z under study: Having approximated g κ,z using the Tikhonov regularization, we expect that 1/∥g κ,z ∥ L 2 ( T ) considered as a function of z is an indicator function for D (see [5,Sections 5.6 and 11.5] for an analysis of the behavior of the norm of g with respect to the sampling point).In Section V, we describe the discrete version of the LSM, which we implement in this work since we only have access to the scattered field at a finite set of directions.For more information on the discrete version of the LSM see [6].

III. EXPERIMENTAL DATABASES
In this section, we present the two open access databases of experimental measurements used in this work.

A. Manitoba Database
A detailed description of the setup can be found in [15] and [18].As depicted in the right panel of Fig. 1, 24 doublelayer Vivaldi antennas are placed in a Plexiglas cylinder with an inner radius of 22.2 cm.Each of them is 7 cm long and held 3 cm apart from the wall, so that its tip is at 10 cm from the inner edge of the wall.All the targets investigated are inside a 12 × 12 cm square, so this shall be the inspection region in our computations.For each scatterer, experiments at eight frequencies linearly spaced between 3 and 10 GHz are performed.
Note that, in this dataset, full aperture measurements (except for backscattering) are provided.However, these measurements have a low signal-to-noise ratio: the 24 antennas are present at the same time, allowing for self-induction (a phenomenon that is not handled in the current model).Furthermore, as is mentioned in [15], the lengths of the wires connecting the antennas to the spectrum analyzer differ, hence having different impedances.The authors recognize this in [15] and attempt to ameliorate this issue by postprocessing of the raw data.

B. Fresnel Database
A thorough description of the experimental setup can be found in [16].As depicted on the left panel of Fig. 1, a 2-D bistatic measurement system is used: an emitter is placed at a fixed position on a circular rail with a radius of 720 ± 3 mm, and a receiver is rotating with the arm around the target (at a total distance 760 ± 3 mm from the center).In this configuration, the angular range of the receivers is outside a 60 • angle on either side of the transmitter: the target is rotated from 0 • to 350 • in steps of 10 • , whereas the receiver rotates from 60 • to 300 • in steps of 5 • (see Fig. 1).We complete the data by setting the unknown measurements to zero; this was the choice in [17], and we provide a justification in Section IV-A.
The operating frequencies are chosen differently for each target under study: from 1 to 8 GHz in steps of 1 GHz or from 4 to 16 GHz in steps of 4 GHz for either one or two rounded dielectrics; or from 2 to 16 GHz in steps of 2 GHz for a metallic rectangular or U-shaped target.

IV. NUMERICAL IMPLEMENTATION
In this section, we shall describe a general implementation of the LSM for finitely many transmitters and receivers.Because we use data completion (see Section IV-A), we can assume that full multistatic data are available.
For this purpose, let us denote by {x T t } N T t=1 and {x R r } N R r =1 the locations of the transmitting and receiving antennas, respectively.In both experimental datasets, the transmitting antennas are linearly spaced in angle on a circumference around the target for t = 1, . . ., N T , and where R T = 0.124 m and N T = 24 for the Manitoba database, and R T = 0.720 m and N T = 36 for the Fresnel one.Similarly, the receiving antennas are located on a circumference for r = 1, . . ., N R , where now R R = 0.124 m and N R = 24 for the Manitoba database, and R R = 0.760 m and N R = 72 for the Fresnel one.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
For each of the scatterers, measurements are provided for several frequencies: the corresponding wavenumbers are denoted as {κ k } N f k=1 .Accordingly, we organize the measurements in a where we recall that u sc (•; x T t , κ k ) denotes the measured the scattered field when the target is irradiated from a point source situated at x T t and for the wavenumber κ k .We denote the set of sampling points by {z s } N S s=1 .These points are usually arranged into a grid over the region to be probed for the scatterer but can form an unstructured mesh, which, as we will see in Section IV-C can be exploited in an adaptive manner.
The near-field equation ( 3) at the sampling point z s and for the wavenumber κ k is rewritten at the discrete level after numerical integration by the trapezoidal rule as follows: where .
This equation is ill-conditioned, so it is solved approximately via Tikhonov regularization: we define g k,s as the minimizer of for some regularization parameter α > 0, where ∥•∥ denotes the Euclidean norm.
For an efficient implementation, it is important to note that, although A k depends on the wavenumber κ k , it does not depend on the sampling point z s .Let us denote by {σ k,n } N n=1 , the singular values of the matrix A k , and by {u k,n } N n=1 and {v k,n } N n=1 , the corresponding left and right (normalized) eigenvectors, respectively: Then, the minimizer g k,s can be computed as where N = min{N R , N T }.For each sampling point z s and wavenumber κ k , we compute the minimizer g k,s = (g tks ) N T t=1 as described above, and the standard single-frequency LSM indicator for the wavenumber κ k at the point z s is defined by It is suggested in [4] to choose α using the Morozov discrepancy principle, and we will follow that approach here using a Morozov parameter δ > 0.

A. Data Completion
Here, we investigate data completion for the Fresnel database (the same method is used for the missing backscattering data in the Manitoba database).
Let us denote by u sc m : [0, 2π] → C the restriction of the scattered field u sc (•; x T t , κ k ) to the circumference r(θ ) = R R (cos(θ )i + sin(θ )j).If we knew u sc m exactly on the subinterval [(π/3), 2π − (π/3)] ⊂ [0, 2π ], we could invoke analytic continuation to extrapolate its values to the whole interval.However, the use of analytic continuation is an ill-posed problem, and hence unstable to data error.
Instead, we can see this as a data-fitting problem since, in practice, we only have measurements of u sc m at a finite set of points: let us denote by , the scattered field for a given angle of incidence and frequency measured for the points x R r = r(θ r ) (r = 1, . . ., N R ).Taking into account that u sc m admits a Fourier expansion, we can try to extrapolate it to the full aperture [0, 2π ] by using a linear least square method with trigonometric polynomials.That is, given a fixed N = 2M +1, we look for c ∈ C N that minimizes the distance and F r n = e inθ r .
This least squares problem is numerically unstable, rendering it infeasible to approximate the first N Fourier coefficients of u sc m in this way.This behavior is clearly shown in the top row of Fig. 2. Alternatively, a common option (c.f.[19]) consists in regularizing the minimization problem via Tikhonov penalization of the solution, i.e., looking for the minimum of where µ > 0 is the regularization parameter.However, recalling Parseval's equality 2π 0 we see that minimizing the norm of c also minimizes the L 2 norm of the extrapolation.This means that, for N fixed and µ increasing, such extrapolation will tend toward 0 everywhere; and for fixed µ > 0 and increasing N , the extrapolation will tend to zero everywhere except at the measuring angles: this is shown in Fig. 2. For N = 2M + 1 fixed, the trigonometric fit drops rapidly to zero away from the data (particularly for larger N ).Of course, in cases where there is significant noise in the data, a larger value of the regularization parameter µ is needed.The above discussion suggests why using zero for unknown data, as done in [8], is a reasonable choice and the approach we have finally decided to use.However, we have no explanation why this choice works so well in practice, as we shall see in Section V-B.In the top row it can be seen that, with no regularization (µ = 0), the extrapolated part diverges as the number of modes M increases.If M is fixed and the regularization parameter µ increases, the solution tends to 0 away from the measurements.Finally, for small but nonzero µ fixed and a large number of modes, the solution looks as if the data were completed with zeroes.

B. Multifrequency LSM
When considering multifrequency data, it is potentially useful to couple information from different frequencies.This is done in [10] by adding a term in (5) which promotes coherence in the phase of the solution g κ,z of (2) across different values of κ.This approach is very interesting, however, we do not use it here.The penalty term added to the Tikhonov functional in ( 5) is dependent on each sampling point z and prevents us from solving the minimization problem for every sampling point in a fast way since the relevant singular value decomposition needs to be recomputed for each sampling point.
To maintain the efficiency of the LSM, we opt for the more classical multifrequency indicators from [9].In particular, we use the so-called serial and parallel indicators given, respectively, by In [9], it is shown that the serial indicator may not distinguish the interior and exterior of the scatterer if one of the frequencies is a transmission eigenvalue for the domain, whereas the parallel indicator has no such issue.We shall show results for both indicators.Note that changes in the amplitude of the incident wave [C κ in ( 1)] and changes due to differences in scattering for different frequencies will cause the terms in the single-frequency indicators to have different magnitude.When we present the single-frequency plots, this different magnitude can be seen in the color bars (e.g.Fig. 3).We return to this point in Section V-B where the effect is strongly apparent.
During our study, we noticed that a good indicator could be obtained by multiplying the single-frequency indicators We call this the product indicator and note that it is related to the logarithm of the parallel indicator This suggests that this indicator is going to have sharper boundaries than the parallel indicator: the rationale for the single-frequency indicator is that it attains values close to zero outside the scatterer and nonzero inside it.With this in mind, the product indicator promotes values that are close to zero whenever at least one of the single-frequency indicators is near zero.Moreover, the product indicator is found to be less sensitive to the normalization C κ of the individual single-frequency incident fields (indeed, for α = 0 different only add an offset to the product indicator).An indicator closely related to ours was previously proposed in [14] but only tested for metallic targets; here, we also consider dielectric targets and discuss its limitations.The product indicator in [14] differs from ours in that it uses the geometric mean so possibly smoothing the result.

C. Adaptive Method
In the computation of the single-frequency indicator, the value of ||g k,s || at each sampling point z s is computed independently of the others.This allows for easy parallelization of  coarse triangular mesh for the inspection zone and gradually refine in places where the relevant indicator varies the most.This approach was suggested in [12] and we test a modified version of that method on real data.The algorithm is as follows: 1) define a coarse triangular mesh on the inspection zone; 2) evaluate ||g k,s || at the vertices of the mesh and compute the chosen indicator; 3) mark the triangles to be refined depending on the variation of the relevant indicator and some given threshold; the procedure stops when there are no more triangles to be refined, and otherwise it proceeds to step 4); 4) refine the mesh; 5) compute ||g k,s || at the newer vertices of the mesh; and 6) go back to step 3).In this iterative way, we are choosing efficiently the points {z s } N S s=1 .One possible indicator for the refinement of a triangle abc in step 3) is max{I a , I b , I c } − min{I a , I b , I c } > t where {I a , I b , I c } are the values of the chosen indicator (singlefrequency, serial, parallel, or product) on the three vertices of a triangle and t is some threshold.Another possible criterion for the refinement of a triangle could be based on the gradient of the indicator on the given triangle, but limited testing suggests the aforementioned criterion results in a grid that is more concentrated on a neighborhood of the boundary of the scatterer compared to a gradient-based rule.

V. NUMERICAL RESULTS
We start by presenting results for the Manitoba database, and then proceed to the Fresnel database.
For each database, we reconstruct two targets that illustrate some of the issues with the multifrequency LSM indicators in a noisy environment.We start by performing a standard single-frequency reconstruction of each scatterer using the Morozov discrepancy principle to choose the regularization parameter [4].In this article, the Morozov noise parameter was chosen by trial and error; in practice, we would need to estimate the noise in the data or to calibrate the measurement system by choosing δ based on a known target.Except for the adaptive results, we always use a 100 × 100 uniform grid of sampling points {z s } N S s=1 .We then present multifrequency results for the serial, parallel, and product indicators.

A. Manitoba Data
We start by reconstructing a single hollow steel cylinder.The individual single-frequency results are shown in Fig. 3 and the multifrequency indicators are shown in Fig. 4. It is tempting to conclude that the LSM is imaging both the outer and inner boundaries of the pipe (we do not have information about the inner diameter), see, for example, the results at 4.5 GHz and also with the multifrequency product indicator.However, there is no theoretical justification for this observation.
Any of the three multifrequency indicators delivers a satisfactory reconstruction, although the product indicator is the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.sharpest and gives an excellent reconstruction of the outer boundary.
We next consider a single large penetrable dielectric Nylon-66 cylinder.In Fig. 5, for some of the single-frequency we again see a ring structure although we believe this object is a solid dielectric.Moreover, the product indicator shown in Fig. 6 fails in this case.The reason is that unfortunately the yellow regions shown for single-frequencies in Fig. 5 (for example at 3.5 and 4 GHz) do not overlap near the boundary of the scatterer and therefore the product is small near the boundary.Instead, the product indicator emphasizes the bright spot at the center.An explanation for this behavior could be that we have transmission eigenvalues close to the frequencies used in the experiment.Indeed, the reported value of the electrical permittivity for Nylon-66 is ε r = 3.00-0.03i.Ignoring the imaginary part, we have calculated transmission eigenvalues at 3.91, 4.03, 4.31, 4.64, and 5.06 GHz (since the imaginary part of ε r is small, continuity of the transmission eigenvalues with respect to changes in ε r suggests that the true eigenvalues are close by).However, attempts to modify the LSM along the lines of [20] did not succeed.
The proximity of transmission eigenvalues to the measurement frequencies may indicate why the parallel indicator in Fig. 6 provides the best reconstruction.Given the wide discrepancy between the results for the serial, parallel, and product indicators, we might suspect that the reconstructions are unreliable.Greatly increasing the regularization parameter to δ = 0.45 provides more consistent reconstructions between the parallel and product indicators, see Fig. 7, and a much better reconstruction of the outer boundary.Indeed, we suggest that mismatch between the point source antenna model and the true device is larger for the Nylon-66 cylinder than for the metal cylinder because the boundary of the Nylon-66 cylinder is much closer to the antennas than the metal cylinder.
The lesson from reconstructing the Nylon-66 cylinder example is that a comparison of both the parallel and product indicators may be necessary to establish the reliability of the reconstruction.

B. Fresnel Data
Next, we turn to the Fresnel database which contains several different targets, each having measurements at several frequencies (not the same in each experiment).We have tried data completion and inversion on all the 2-D targets but elect only to show results for two of the more complex targets.
First, we show results for a U-shaped metallic scatterer in Figs. 8 and 9.As is to be expected, for low-frequency measurements the scatterer is reconstructed without good definition.Moreover, the single-frequency LSM can have difficulties clearly resolving that the scatterer is an open U, and attempts to improve the multifrequency serial and parallel reconstructions lead us to our new product indicator.At each frequency, there are always some artifacts visible inside the scatterer.However, the position of the artifacts varies with frequency (whereas the U-boundary does not) and this is one motivation for the product indicator.
Remarkably, the higher frequency reconstructions start to show artifacts outside the scatterer.These artifacts are unlikely to be due to the measurement array being too coarse.Indeed, Fig. 12. From left to right we show the initial coarse mesh, and then five adaptively generated meshes together with the resulting product indicator.The final reconstruction is comparable to the product indicator result in Fig. 9. Fig. 13.Adaptive procedure for the Two-diel target.The reconstruction after five iterations is comparable to the product indicator in Fig. 11.
in [21] it is recommended for a wavenumber κ to take an angular separation of at most π/(κa) radians, where a = 0.047 m is the radius of the smallest circumscribing circle around the object.For the largest frequency ν = 16.0GHz this evaluates to 0.199 rad whereas the measuring receivers are separated by 5 • = 0.087 rad, so the artifacts are unlikely to be due to a coarse sensor array.
In Fig. 9, we show the results of using the serial, parallel, and product indicators defined in Section IV-B.The product indicator helps to clear up the artifacts inside the U compared to those shown in the serial and parallel results (and it removes completely those outside the U in the parallel result).
The second example from the Fresnel database consists of two disconnected dielectric cylinders with a measured relative electric permittivity of ε r = 3 ± 0.3.This example is chosen to illustrate the behavior of the LSM for disconnected targets.
In Fig. 10, we show individual reconstructions at each frequency.As is to be expected, at low frequency it is not possible to distinguish the two objects, but as the frequency increases the separation of the two scatterers becomes clear.The three multifrequency indicators under study are shown in Fig. 11.It is interesting to note that the reconstructed scatterers are slightly rotated from the published exact case, and this has been noticed too in [16] and [22].It is also interesting to emphasize that the product indicator again sharpens the reconstruction.
As we have already mentioned the amplitude C κ of the incident field is not the same among all the frequencies as can be seen in Table I.Taking into account that the scattered field is linear with respect to the incident field, we expect its amplitude to differ between frequencies due to a combination of weaker/stronger scattering as well as the variation in emitted power.This will also affect the single-frequency indicators.
In general, the higher the incident amplitude the higher will be the single-frequency indicator values.This is shown by the differing range of the color maps in the monofrequency plots.For example, in Fig. 10, the indicator for frequency 1.0 GHz is much smaller than for the other frequencies perhaps due to weak scattering at low frequency.Although the single-frequency colormaps are the ones one would use when performing single-frequency calculations, it is important to take into consideration the relative amplitudes when computing multifrequency indicators.Fig. 10 explains why the parallel indicator performs better than the serial indicator for this target: on one hand, the parallel indicator depends on the sum of the monofrequency indicators.While the first two monofrequency indicators yield the worst reconstructions, they are also the ones with the smallest Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
amplitudes, so they do not have a great influence on the parallel indicator.On the other hand, the serial indicator depends on the sum of the inverse of the monofrequency indicators, hence these two first reconstructions play a big role in determining the outcome.Finally, the product indicator is almost immune to different rescaling of the monofrequency indicators.
Our final results test the adaptive procedure on the Fresnel examples.We chose the product indicator in the adaptive algorithm outlined earlier.We use t = 0.3 for the U-shaped scatterer and t = 0.5 for the two dielectrics.The initial coarse mesh, and subsequent refinements are generated using Netgen [23].For the U-shaped scatterer (see Fig. 12) and the two dielectric cylinders (see Fig. 13) we show five refinement steps.In the case of the U-shaped scatterer the reconstruction improves with each step.For the two dielectrics, the reconstruction does not improve after the third refinement step.
The adaptive approach can reduce the number of evaluations of the indicator function for each frequency significantly and, hence, reduce the computational time.For example, for the two dielectrics and using ten iterations of the Newton solver for the Morozov criterion, the use of a uniform 100 × 100 grid (see Fig. 9) requires 10 000 evaluations of the indicator for each frequency and the entire multifrequency reconstruction takes 42 s CPU time running vectorized in Python on a single-core AMD Ryzen 7. By contrast, the adaptive approach, after five refinement steps, uses 250 sampling points and an overall CPU time of 6 s running on the same computer.The latter includes time for mesh initialization and refinement, and we emphasize that the adaptive scheme is not fully optimized since the indicator is computed at each mesh vertex for every iteration (i.e., without taking into consideration the already known values).

VI. CONCLUSION
We have provided results for the multifrequency LSM using serial, parallel, and product indicators on real data, and used an adaptive approach to improve CPU time for the associated reconstructions.Generally, the product indicator provides a sharper reconstruction as long as the individual single-frequency reconstructions conform to the LSM theory in that the indicator does not vanish near the boundary of the scatterer.In our examples, the parallel indicator always provides a better reconstruction than the serial indicator.For a more reliable reconstruction of the scatterer, both the product and parallel indicators should be compared before assessing a reconstruction.
The adaptive approach shows promise for speeding up the LSM by reducing the number of indicator evaluations needed, at least in cases where the scatterer is well resolved.Note that both the uniform grid case and the adaptive case could benefit from parallelization since the evaluation of the indicators for each sampling point is embarrassingly parallel.But this is beyond the scope of our work.
Multifrequency data can greatly improve the LSM reconstruction using relatively simple and easily computed indicators.

Fig. 1 .
Fig. 1.Experimental setup for both databases (left: Fresnel, right: Manitoba).The location of the transmitting antennas is shown with red crosses, whereas for the receiving ones it is shown in black.It is important to realize that the whole Manitoba setup fits inside the inspection zone of the Fresnel database (shown in lavender color in both cases).

Fig. 2 .
Fig.2.Several extrapolations of the same measurements in the Fresnel dataset.In the top row it can be seen that, with no regularization (µ = 0), the extrapolated part diverges as the number of modes M increases.If M is fixed and the regularization parameter µ increases, the solution tends to 0 away from the measurements.Finally, for small but nonzero µ fixed and a large number of modes, the solution looks as if the data were completed with zeroes.

Fig. 3 .
Fig. 3. Single frequency LSM results for the hollow steel cylinder, using the Morozov discrepancy principle to choose the regularization parameter α; the Morozov noise parameter is fixed to δ = 0.00001.The bottom right panel shows the true target.The colormap is individual to each frame.

Fig.
Fig. Multifrequency LSM results for the hollow steel cylinder.For each individual frequency, the Morozov discrepancy principle is used to choose the regularization parameter, and the Morozov parameter is fixed to δ = 0.00001.In the bottom right panel, we show the exact scatterer together with bars indicating the maximum and minimum wavelength of the probing field.

Fig. 5 .
Fig. 5. Single frequency LSM results for the Nylon-66 cylinder, using the Morozov discrepancy principle to choose the regularization parameter α; the Morozov noise parameter is fixed at the same value as for the steel cylinder, δ = 0.00001.

Fig. 6 .
Fig.6.LSM results for the Nylon-66 cylinder.For each individual frequency, the Morozov discrepancy principle is used to choose the regularization parameter (for the noise parameter fixed to δ = 0.00001).In the bottom right panel, we show the exact scatterer together with bars indicating the maximum and minimum wavelength of the probing field.

Fig.
Fig. Multifrequency indicators for the U-shaped target.These combine the results from Fig. 8 and show the power of the product indicator to reduce artifacts.

Fig. 10 .
Fig. 10.Single frequency LSM indicator for the two dielectric cylinders using the Morozov discrepancy with the noise parameter fixed to δ = 0.00001.

TABLE I INCIDENT
FIELD AMPLITUDE MEASURED AT THE ANTENNA IN FRONT OF THE EMITTER AT DIFFERENT FREQUENCIES FOR THE FRESNEL DATA