1 Introduction
In this paper we aim to numerically solve the inverse scattering problem for periodic media in . The periodic media of interest are unboundedly periodic in the horizontal direction and bounded in the vertical direction. These periodic media are motivated by one-dimensional photonic crystals and the inverse problem of interest is inspired by applications of nondestructive evaluations for photonic crystals. There have been an increasing amount of studies on numerical methods for shape reconstruction of periodic scattering media during the past years, see [Arens2005, Arens2003a, Bao2014, Elsch2012, Hadda2017, Lechl2013b, Nguye2014, Sandf2010, Yang2012, Nguye2016, Jiang2017, Cakon2019, Nguye2022b, Nguye2022]. Two major approaches that were studied in these papers are the factorization method and the near field imaging method. The latter method, which relies on a transformed field expansion, can provide super-resolved resolution. However, this method requires the periodic scattering structure to be a smooth periodic function multiplied by a small surface deformation parameter. The factorization method can essentially work for periodic scattering structures of arbitrary shape but it is not very robust against the noise in the scattering data. This method belongs to the class of sampling or qualitative methods that were introduced by D. Colton and A. Kirsch [Colto1996, Kirsc1998]. The factorization method aims to construct a necessary and sufficient characterization of the unknown scatterer from multi-static data. We refer to [Kirsc2008] for more details about the factorization method.
In this work we develop a sampling method with a new imaging functional to solve the inverse scattering problem for periodic media. This sampling method is inspired by the orthogonality sampling method [Potth2010] and the direct sampling method [Ito2012]. These two sampling methods share similar ideas and features and were studied independently. For simplicity we will refer to them as the orthogonality sampling method. The computation of the new imaging functional is simple and fast as one only needs to evaluate two finite sums involving the propagating modes of the scattered field data. Like the orthogonality sampling method, the proposed sampling method is also regularization free and very robust against noise in the data. The resolution of the new imaging functional is studied using Green’s identities and the Rayleigh expansion of the -quasiperiodic fields of the scattering problem. The stability of the functional is also established. The performance of the new imaging functional is studied in various contexts in the numerical study. The numerical study also shows that the proposed sampling method is more robust than the factorization method and more efficient than the orthogonality sampling method in reconstructing periodic scattering media. We also want to mention that although the orthogonality sampling method has been studied for inverse scattering from bounded objects [Potth2010, Gries2011, Ito2012, Ito2013, Kang2018, Harri2020, Le2022], its application to the inverse scattering problem for periodic media is still not known.
The paper is organized as follows. The basics of the scattering from periodic media and the inverse problem of interest are described in Section 2. The new imaging functional and its resolution and stability analysis are discussed in Section 3. Section 4 is dedicated to a numerical study of the new imaging functional and its comparison to the factorization method and the orthogonality sampling method.
2 Problem setup
We consider a two-dimensional medium which is unboundedly -periodic in -direction and bounded in -direction. Let be a bounded function which is -periodic with respect to . Suppose that the interior of the periodic medium is characterized by and that the exterior of the periodic medium is homogeneous which means in these areas. Note that the period can be any arbitrary value, but it is chosen to be for the convenience of the presentation. For , we define that a function is -quasiperiodic in if
From now on we will call functions with this property -quasiperiodic functions for short. A typical example of -quasiperiodic functions is a plane wave (e.g. with , ). Suppose that the periodic medium is illuminated by an -quasiperiodic incident field with wave number . Note that since the medium is unboundedly periodic in , we are only interested incident fields propagating downward or upward toward the medium. The scattering of this incident field by the periodic medium produces the scattered field described by
(1) |
where is the contrast given by
It is well known for this scattering problem that the scattered field must also be -quasiperiodic, and that the direct problem of finding the scattered field can be reduced to one period
Let . For such that
(2) |
the direct scattering problem is completed by the Rayleigh expansion condition for the scattered field
(3) |
where
and are the (complex-valued) Rayleigh sequences of the scattered field . The condition (3) means that the scattered field is an outgoing wave. Note that only a finite number of terms in (3) are propagating plane waves which are called propagating modes, the rest are evanescent modes which correspond to exponentially decaying terms. From now, we call a function satisfying (3) a radiating function. In addition, we also assume that is nonzero for all which means the Wood anomalies are excluded in our analysis.
Well-posedness of scattering problem (1)–(3) is well-known, see for instance [Bonne1994]. For define
Recall the constant in (2). For the inverse problem of interest we measure the scattered field on for some . From the Rayleigh expansion of , knowing on is equivalent to knowing the Rayleigh coefficients . Since the evanescent modes are associated with the exponentially decaying terms in the Rayleigh expansion (3), it is typically difficult to obtain these modes in practice unless one can measure extremely near (e.g. within one wavelength) the periodic scatterers. Thus, we consider only the propagating modes for the scattering data of our inverse problem as follows.
Inverse problem. Given the Rayleigh coefficients for such that , determine .
3 The imaging functional and its properties
Recall that the -quasiperiodic Green function of the direct problem is given by
(4) |
It is well known that the direct problem is equivalent to the Lippmann-Schwinger equation
(5) |
Let be the number of incident fields we use for the inverse problem. Let be the Rayleigh coefficients of the scattered field generated by incident field for . For , define the following imaging functional
where
This functional aims to determine in and plays the role of sampling points. We note that are also the Rayleigh coefficients of the -quasiperiodic Green function. In the proof of following theorem we will drop the dependence of on for the convenience of the presentation but we keep in the total field that is generated by incident field . The following theorem is the main theorem justifying the behavior of as an imaging functional.
Theorem 1.
For , the imaging functional satisfies
where
(6) |
where and are respectively the Bessel functions of first and second kinds.
Remark 2.
The series in (1) converge since is always nonzero for and . We know that and decay like as . Thus the terms with larger in (1) are small terms. The terms with smaller can be also small with a suitably larger wave number (e.g. the choice of works well in our numerical study). Therefore, the integral associated with is expected to be small compared with the integral with the kernel . Thanks to the behavior of the latter integral is expected to have small values as is outside and to have much larger values as is inside .
Proof.
For any , solves
Multiplying both sides by and integrating over gives
(7) |
Similarly, solves
thus by multiplying both sides by and integrating over we obtain
(8) |
Subtracting (8) from (7) yields
where
is an unit normal outward vector to
. Note that , but thanks to the -quasiperiodicity ofTherefore
(9) |
The left-hand side of (9) equals , now we calculate the right-hand side. Since satisfies the radiation condition, for a fixed ,
(10) |
and
Therefore
and since
we have
(11) |
Similarly
(12) |
Combining (11) and (12) yields
that is
(13) |
From the Rayleigh expansion (3) we have that
Thus, for , substituting (5) into the integral above gives
For , using the Rayleigh coefficients of above and (13) we obtain
Now the proof follows from using another representation of the -quasiperiodic Green function
where the series converges for , . ∎
In the next theorem we will establish a stability estimate for the imaging functional. Assume
for all .Theorem 3.
For , denote by and the noisy scattered wave and its Rayleigh sequences respectively, for which we have
Define
Then the following stability property holds
for every .
Proof.
For and such that , we have
thus, by Cauchy-Schwarz inequality,
Therefore
where
For brevity, set
then we have
for all . Hence, for all ,
This completes the proof. ∎
4 Numerical study
In this section, we test the performance of the proposed sampling method with respect to the number of incident sources, the levels of noise in the data, the wave numbers, the values of parameter (in (3)) and exponent , and the shape of the periodic scatterers. For the latter category we compare the performance of the proposed sampling method with those of the factorization method and the orthogonality sampling method.
In all of the numerical examples below we choose the following parameters:
The sampling domain is uniformly discretized in each dimension with sampling points. The boundary measurements are discretized uniformly with 64 points on each boundary. By using incident sources we mean to consider
where sources are uniformly located on and sources are uniformly located on . We generate the synthetic scattering data by solving the direct problem with the spectral Galerkin method studied in [Lechl2014]. With artificial noise added to the synthetic scattering data we implement the imaging functional
To compare with the orthogonality sampling method we implement following imaging functional
This imaging functional of the orthogonality sampling method can be rewritten in the modal form by using the Rayleigh expansion of the scattered field. In this modal form the evanescent modes or the exponentially decaying terms can be neglected and there remains only a finite number of the propagating modes. If we drop in , we will approximately obtain the modal form of . We refer to [Nguye2014] for the imaging functional of the factorization method implemented for the test in this section.
4.1 Reconstruction with one incident source (Figure 1)
In this section we test the performance of the sampling method for data generated by only one incident source. Here the parameters are chosen as , wave number , and exponent . We add 20 artificial noise to the scattered field data (). From Figure 1 we can see that the method is able to reconstruct small scatterers quite well. This capability is an advantage over classical sampling methods (e.g. linear sampling method, factorization method) in terms of computational efficiency since classical methods are only able to reconstruct targets with data generated by multiple incident fields. However, the proposed sampling method fails to reconstruct scatterers with extended shape, which is reasonable since only one incident source and one wave number are used for the data. We refer to Figure 2 for improved results when multiple incident sources are used for the reconstruction.
4.2 Reconstruction with multiple incident sources (Figure 2)
In this section we test the performance of the method for reconstructing extended scatterers with different number of incident sources. The parameters are the same as in Figure 1, meaning , wave number , exponent , and 20 artificial noise is added to the scattered field data. If in Figure 1 the method fails to reconstruct extended ellipses with one incident source, the reconstructions are improved with more sources, see Figure 2. The reconstructions already look reasonable with 32 incident sources and continue to improve when more sources are used. The results remain almost the same even if more than 128 incident sources are used to generate the scattering data.
![]() |
![]() |
![]() |
![]() |
4.3 Reconstruction with different levels of noise in the data (Figure 3)
In this section we test the performance of the method for reconstructing extended scatterers with different levels of artificial noise in the data ( and ). Here the parameters are chosen as , wave number , exponent , and 128 incident sources are used to generate the scattering data. It can be seen from Figure 3 that all reconstructions are not affected by the amounts of noise added to the data. Furthermore, the reconstructions will remain essentially the same even with much higher amounts of noise in the data. This is not a surprise since the evanescent modes are typically sensitive with noise but the sampling method only uses propagating modes. The great robustness against noise in the data was also seen in the orthogonality sampling methods studied in [Potth2010, Harri2020, Le2022].
![]() |
![]() |
![]() |
![]() |
4.4 Reconstruction with different wave numbers (Figure 4)
In this section we test the performance of the method for reconstructing extended scatterers with different wave numbers (). Here the parameters are chosen as , exponent , 128 incident sources are used to generate the scattering data, and noise added to the data. We can see from Figure 4 that the resolution of reconstructions improve as the wave number increases, which is reasonable. However, it is also known that as the wave number increases, the inverse problem will become more difficult to deal with. We can notice some small effect of larger wave numbers in the reconstruction for .
![]() |
![]() |
![]() |
![]() |
4.5 Reconstruction with different values of and (Figures 5–6)
In this section we test the performance of the method for different values of the parameters and . The other parameters are chosen as , 128 incident sources are used to generate the scattering data, and noise added to the data. From Figures 5–6, we can see that the sampling method works well for different values of . Although the reconstructions look reasonable for exponents and , the exponent seems to be an ideal exponent for the imaging functional. With larger the reconstructions will have cleaner background but lose some small details of the scatterers.
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
4.6 Reconstruction of different shapes and comparison with the factorization method and the orthogonality sampling method (Figures 7–10)
In this section we test the performance of the proposed sampling method for different shapes of periodic scatterers and compare its performance with those of the factorization method and the orthogonality sampling method. The parameters are chosen as , wave number , exponent , 128 incident sources are used to generate the scattering data, and noise added to the data. It is obvious from all reconstructions in Figures 7–10 that the factorization method suffers severely from the amount of noise added to the scattering data. Actually, the reconstructions of the factorization method can be greatly affected even with smaller amounts of noises (e.g. or ). This unstable behavior of the factorization method in reconstructing periodic media was also reported in [Arens2005, Nguye2016]. The reconstructions of the orthogonality sampling method are as stable as those of the proposed sampling method but it is also clear from the pictures that the accuracy in the reconstructions of the orthogonality sampling method is much worse than that of the proposed sampling method. The proposed sampling method may provide reasonable reconstructions for different shapes considered in the test. This indicates a high efficiency of this sampling method.
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
![]() |
Acknowledgement. This work was partially supported by NSF grant DMS-1812693.