The problem of diagnosing a grid of small (in terms of the probing wavelength) dielectric scatterers is considered. The aim is to detect and locate possible defects occurring within a known grid when one (or more) scatterer is removed/missing (fault). The study is developed for the canonical case of a TM scalar two-dimensional geometry with the scatterers consisting of dielectric cylinders of small circular cross section. The scattering by a fault is modeled by relaying only to a priori information about the complete grid which leads to a numerically effective inversion procedures as the bulk of the numerical effort is to be done only once. Inversion is achieved by a truncated singular value decomposition scheme and results are provided in terms of closed form expressions for the probability of detection and of false alarm. This allows us to foreseen the achievable performance and to highlight the role of scattering configuration parameters. Numerical examples are also enclosed to corroborate theoretical outcomes. The case of two or more faults is considered as well. For such a case it is numerically shown that detection method still works well even though multiple scattering (occurring between faults) is neglected.
© 2015 Optical Society of America
Grid structures arise in different applicative contexts in optics and electromagnetics literature. A periodic lattice can be designed in order to force wave propagation within a given range of directions and frequencies . Photonic bandgap crystals [2, 3] and metamaterials  are examples that find application in high gain antennas [5, 6] and waveguides  design. The planned electromagnetic response of the structure depends on precise positioning and dielectric properties of the grid elements. The absence of one or more elements of the grid leads to changing the behavior upon which the structure has been devised. Even though literature has paid attention on developing numerical methods in order to predict the field scattered by the grid when such a circumstance occurs [8–10], the need to test such structures in a nondestructive way arises naturally.
In this paper the problem of diagnosing a grid of dielectric small scatterers is addressed. In particular, the study focuses on the canonical TM scalar case where the scatterers consist of small (in terms of wavelength) circular cross section dielectric cylinders.
The problem is cast as an inverse electromagnetic scattering problem where the goal is to detect possible defects occurring in the grid from scattered field measurements taken outside the scattering structure. Such defects can be due to damage or even to wrongly manufacture. In particular, as defect the lack of some scatterer that should have crowed the grid is considered. Such missing scatterers are herein termed as ’faults’.
In order to attack the problem two different strategies can be followed: one can search for all the scatterers that populate the grid (determining where faults occur by comparison with the full ideal grid) or one can straightly seek for the missing scatterers . It is clear that the latter approach appears more convenient. Indeed, by searching for the faults poses less severe problems due to mutual scattering (as the fault number will be certainly much lower than the actual scatterers) and hence to the corresponding non-linearity of the inverse scattering problem. Moreover, this approach allows to tackle the inverse problem under a more favourable ratio between independent data - unknowns, which can be a severe limitation when searching for the scatterers (instead of faults) due the finite number of degrees of freedom of the scattered field . By contrary, searching for the missing scatterers requires adopting a suitable scattering model that must take into account that the fault scatter within an inhomogeneous background medium .
In this paper, the grid diagnostics is pursued by detecting and localizing the faults occurring within a known grid.
In this regard, the main contribution of this manuscript is the development of a scattering model which allows to reduce the numerical complexity of the inversion procedure by a factor of N, N being the number of scatterers populating the whole ideal grid. In particular, in the case of a single fault this model is exact whereas when multiple faults occur it neglects mutual scattering between them.
The second aim pursued in this paper is to foreseen in closed form the achievable performance in the detection problem. To this end, we find convenient to adopt a hybrid two steps based strategy. In the first step the scattering model is inverted by a truncated singular value decomposition scheme (TSVD)  once the unknown fault position has been represented as the support of a delta function. In the second step, a classical detection scheme is applied. The proposed method allows to foreseen the achievable performance in terms of probability of detection (PD) and of false alarm (PFA) when data are corrupted by an additive zero mean white Gaussian noise. What is more, the role of the scattering configuration and of the noise is explicitly highlighted.
The case of two or more faults is addressed as well. For such a case the mutual scattering (between the faults) is neglected so that the derived result turns to be not exact. Nonetheless, in view of the adopted problem formulation model error is not too high so that the analytical estimation still works very well.
In  the same geometry as here was adopted but the Born approximation (with respect to the full scattering grid) has been exploited to linearize the problem, even in the case of a single fault. Moreover, even though noisy data has been used to check the proposed inversion method the role of the available SNR (Signal to Noise Ratio) on the achievable performance was not clearly pointed out.
In , the authors considered the same problem with the grid consisting of metallic scatterers. There, a suitable scattering model was introduced (in terms of proper equivalent magnetic currents) and formulas for PD and PFA were derived as well.
In this regards, this contribution can be considered as a generalization of such previous papers.
2. Scattering by a fault occurring in a dielectric grid
Consider a grid scatterer Σ consisting of N identical dielectric cylinders embedded within a homogeneous host medium. The scatterers’ axes are directed along z and their cross section is circular with radius a ≪ λ, λ being the hosting medium wavelength. Denote with εr their relative dielectric permittivity with respect to the host medium whereas their centers are located at the points rn which are numbered according to the progressive index n, as shown in Fig. 1. The illumination is provided by a z-directed filamentary electric current of amplitude I located at rs. Therefore, we tackle the fault detection problem for a two dimensional scalar TMz configuration where only the z directed electric field component is relevant.
Assume that one of the scatterers occupying the cross section Ωm is removed. In this case we say that a fault has occurred (in correspondence to the removed scatterer) in the ensemble Σ and the background field is perturbed by a scattering contribution given by
The field Es can be though as being radiated by the equivalent polarization current just supported over Ωm, so thatEq. (3). However, as we are going to detail in sequel, such a model is ”somehow uncomfortable” in view of numerical implementation when N is very large. To cope with this drawback, in next sections, model in Eq. (3) is cast in alternative more convenient way.
2.1. Green’s function computation
According to the addressed geometry, Green’s function can be conveniently expressed in terms of cylindrical harmonics (which fit well the geometry of each scatterer). In this subsection we just recall some basic facts and notation concerning Green’s function computation which are useful in the sequel.
As the scatterers’ cross sections are small, their scattered fields are well approximated by only the leading term of the corresponding cylindrical harmonic expansions . Accordingly, the Green’s function of the full grid is given byEqs. (4) and (5) is the Hankel function of zero order and second kind and J0(⋅) is the Bessel function of zero order. Note that reciprocity arguments, that is GN(r,r′) = GN(r′,r), have been exploited in Eq. (4) and Eq. (5) derivation.
The expansion coefficients an,bn can be found once the appropriate continuity conditions are forced on each cylinder’s surface. This requires the solution of the following system of linear equations 
The elements of the N × N matrix are given by:
2.2. The Scattered Field
When the m-th cylinder is removed, the harmonic expansion coefficients of the scattered field due to the remaining N −1 scatterers change with respect to the full grid situation, and the total field at rm due to a source located in rs is given by:Eq. (6), the first two terms in the square brackets can be rearranged so that:
2.3. The scattering model
By using the same approach as followed in  for PEC scatterers, here it is shown that the second term in Eq. (13) is proportional to the first one by a factor depending only on the fault position and not on source or observation positions, so that the scattered field can be put into a form more suitable for the inversion.
Looking at the boundary conditions in Eq. (6), it can be easily shown that they can be rearranged as:Eq. (16) eliminating the m-th row and column from matrix . The solution of the new equation system, presenting N − 1 equations and N − 1 unknowns, provides the correct scattering coefficients. However, the same result can be obtained by simply forcing to zero the elements Aij,i = m or j = m and i ≠ j and leaving unchanged the number N of rows and columns of matrix .
To this end, the ”modified ” matrix is introduced, where:
In this way the actual coefficients in presence of the fault can be calculated as:
Equation (20) returns the actual (relevant to the N − 1 grid) scattering coefficient but the m-th one. The latter is incorrect as it is instead of 0. However, it can be ignored, because as shown below, it does not affect the final result.
Since the only elements of matrix different from zero belong the m-th row and column, it is straightforward to see that its product for returns a vector whose elements are all proportional to am, except for the m-th that shows a more complicated expression:
Now, the summation appearing in right-hand side of Eq. (13) can be recast under matrix notation as:
In order to cast the expression above in a more suitable way we exploit once again boundary conditions in Eq. (6) that provide the following relationship between am and bm:Eq.(24) in Eq. (13):
Eventually, Eq. (26) represents that scattering model upon which the inversion/detection procedure will be based.
3. Inversion and detection
Equation (26) can be rewritten as:
Let us denote with NO and NS the number of observation and of source positions, respectively. Accordingly, Eq. (29) gives rise to the following discrete linear operator:
Now we come back to the claim we have anticipated above about the numerical convenience of the model provided by Eq. (30).
Then, by performing similar rearrangements as in Eq. (30), it is obtained:
Even though Eqs. (30) and (32) look similar, actually they are profoundly different from the numerical (complexity) point of view. Indeed, while only requires the knowledge of the full grid coefficients bm, relies also on the fault location as the latter varies over all the grid. Accordingly, in order to build up the forward problem has to be solved only once (of course for each source/observation position). Instead, also depends on EN−1(rn,rs) and hence N2 forward problems must be solved. It is clear that all these efforts can be accomplished off-line. Nonetheless, as N grows up the computation can become prohibitive and the formulation in Eq. (30) with no doubts gives great advantages.
Let us turn to the detection problem. It is clear that the problem at hand is a classification problem where from the observation one has to choose the ”right” belonging to the state space , with having only the i-th entry different from zero. Assuming an additive white Gaussian noise and also that prior about S are not informative, then Bayes criterion (under the hit or miss or even the quadratic loss function) can be stated as
Note that the solution space is now and hence it is broader than the state space S. Also, the classification is now cast as an estimation or equivalently as an imaging problem. It is well known that the inversion of matrix is ill-posed. Therefore, some regularization scheme is required. This entails that, while under ideal situation Eqs. (33) and (34) would return exactly the same , due to noise and to necessity to adopt a regularization scheme solving for eq. (34 will in general retun a solution which is spread out over all the grid. Therefore, the need for a detection stage arises naturally.
In order to regularize the reconstruction a Truncated Singular Values Decomposition (TSVD) inversion scheme is employed. To this end, first the matrix is diagonalized by means of its singular spectrum, that is , where the upper-script H denotes conjugate transpose. Then, the reconstruction is obtained by truncating the SVD representation of the unknown by retaining the first NT singular values.
More in detail, say the projections of the unknown vector onto the i-th singular vector (i.e., the i-th column of ), then under ideal case the reconstruction is provided by the series:
When data are affected by noise , the actual data becomes
As can be seen, two terms contribute to the reconstruction: the first one accounts for the truncation filtering effect, and depends on the actual fault position (i.e. on the ci values for i ≤ NT); the second one is due to the additive noise present on the data. As the truncation index NT increases, the first term tends to the ideal value but the second term increases as well. Accordingly, truncation index NT should be chosen properly in order to compromise accuracy and stability.
Of course different regularization methods can be adopted as well. They basically result in a different filtering of the singular spectrum . Therefore, while the following discussion holds true provided to replace the truncation index NT with a suitable regularization parameters,TSVD is preferred for its simplicity.
Let us define the so-called spatial content  and side-lobe content. The spatial content is defined asEq. (37) calculated where there is not a fault, i.e. the reconstruction that it would be obtained from noiseless data at locations different from the fault itself.
If the noise is a zero mean Gaussian process whose variance is , it results in the reconstruction as a zero mean Gaussian error with variance . As a consequence, the reconstructed function is still a Gaussian distributed random process whose mean value is provided by SC(n) or SL(m,n) depending on the actual fault position, and its amplitude distribution is provided by a Rician probability density function:
Accordingly, in order to detect an locate the fault a further stage consisting of a classical detection procedure can be run. Therefore, once the reconstructed image is obtained, a second threshold Ath, the detection threshold, is introduced and the position of the fault is recovered as the locus of points where .
The probability of detection is defined as the probability that the amplitude of the reconstructed image has value above the detection threshold where the actual profile shows a fault. The probability of false alarm is the probability that the amplitude of the reconstructed image has values above the detection threshold where there is not a fault. It can be shown that, in the case at hand, their expression are, respectively:
4. Numerical analysis
In this section, we check the theory developed above by means of some numerical examples. In particular, we analyze the effects of noise, truncation index and detection threshold level on PD and PFA.
We consider a grid of N = 11×11 identical dielectric cylinders of εr = 4, radius a = 0.034λ, λ being the relevant wavelength, arranged uniformly at a step of Δ = λ/2 over a square of side 5λ. The grid is illuminated by three sources positioned around it and the observation points are equally spaced over a circle as it is shown in Fig. 2.
The singular value behavior is reported in Fig. 3. As it can be appreciated they decrease faster and faster as their index grows up. So, a truncation index NT = 60 corresponds to cutting out singular values whose amplitude is 10dB lower than the maximum one, whereas NT = 87 and NT = 112 correspond to 20dB and 50dB respectively. The truncation index affects both SL and SC, that give rough information (say qualitative) about the chance of succeeding in the detection and the occurrence of false alarms. Indeed, where SC is higher it is more likely to succeed in the detection. On the other hand, false positives are more likely to occur where SL is high. In Fig. 4 SC is shown for different values of the truncation index: for the case at hand, as it was expected, it assumes higher values on the external part of the grid and tends to increase everywhere as the truncation index increases. As for the side-lobe content SL, in Fig. 5 it is shown for a fixed truncation index and two different positions of the fault: central and external. As can be appreciated, the side-lobe content remains lower than the corresponding spatial content value for most of the locations.
Next, let us study the behavior of PD(n) and PFA(m,n) with respect to the detection threshold level Ath for different values of the data Signal to Noise Ratio (defined as ) and of the truncation index NT. In particular, noise is added to the data before any processing takes place. For exemplification we fix n = 61 (corresponding to the central position inside the grid) and choose the value of m where PFA achieves its maximum, so that the worst case is considered.
In Fig. 6 PD and PFA are plotted as function of Ath for different values of NT. As can be appreciated, when a large number of singular functions are retained in the reconstruction (i.e. the truncation level is low), it is difficult to obtain a high PD with a low PFA. Moreover, increasing the detection threshold entails a reduction of PFA but PD decreases as well. This of course is due to the ill-posedness of the inverse problem. Indeed, as the SVD truncation threshold lowers smaller singular values are retained and the variance of the noise in the reconstruction grows up. On the contrary, when the the singular values are truncated at a higher level a favorable situation occurs as an optimal interval for Ath exists where PD approaches unity and PFA is nearly zero. We call this interval as the detection interval. Finally, if an excessive limitation is applied on the number of retained singular values, there is not margin for conveniently choosing Ath. In addition, the actual value of SNR affects performances only when a low truncation level is applied, as can be appreciated by comparing top and bottom lines in Fig. 6.
Eventually, Eqs. (41) and (42), or equivalently the graphical counterparts reported above, as involves all the key ingredients od the detection problem (measurement configuration parameters, grid structure, SNR, TSVD and detection thresholds) provide a tool to a priori evaluate the achievable performance and to choose the more suitable detection setup in order to succeed in the detection.
In order to check this statement, some reconstruction examples are considered below for a fixed SNR = 10dB and different detection threshold values. First, we consider a fault located at n = 61, i.e. in the central position. In Figs. 7 and 8 results concerning the case NT = 60 and NT = 112 are shown respectively. As can be seen, as long as the detection interval exists and the detection threshold is chosen within it (look at Fig. 6) detection is successful. On the other hand, false positives appear or the relevant scatterer is missed. The same results are obtained for defects located in different positions of the grid. Note that, coherently with the theory, in case of low truncation level (Fig. 8) many false alarms affect the result for each detection threshold value.
We pass to consider the case when two or more faults occur jointly. For such cases, mutual interactions (between the faults) take place so the adopted model is no more exact. The aim is to check the ”robustness” of the presented detection scheme when it is forced for situations beyond the model. As mutual scattering depends on the fault number and their relative positions we consider different fault arrangements (see Figs. 9 and 10. As can be seen, also in these cases theory works well. Basically, the reason why the proposed method still works under mutual scattering is that the high order effects for the considered cases, remarkably even for very close faults, are low. Indeed, computing the relative mean quadratic model error one finds that it ranges between 0.2% and 0.3%. However, this could be expected, since, as we have anticipated in the introduction, this is a characteristic of casting the problem of grid diagnostics by looking for faults instead of detecting the scatterers which actually remained in the grid. Finally, it is noted that this is a further point in favor of our methods (wrt classification) as the complexity of procedure is the same also in the multi fault case.
In this paper we have considered the problem of detecting and localizing faults occurring within a known grid of dielectric cylinders in the framework of a scalar two-dimensional geometry.
To this end, first an exact scattering model for one single fault has been introduced. This, on one side, led to representing the scattering in terms of the full grid Greens function avoiding thus the need to recompute it as the fault position is changed. On the other side, this suggested to search directly for the fault instead of the remaining (within the grid) scatterers. In particular, as usual, the number of faults is a small fraction of the elements forming the grid, this allows to mitigate the nonlinearity of the problem and to better exploit information given by scattered field data. Moreover, the proposed model reveled computationally effective while building up the scattering operator to be inverted as it requires solving the forward problem only once (of course for each source position).
In the case of a single fault, the detection problem turned to be linear. Accordingly, the discrete scattering operator has been inverted by means of a TSVD scheme. In particular, analytical expressions for the PD and PFA have been found which accounts for the role of all the relevant parameters. Therefore, they can be used to foresee the achievable performance or to suitably set up the scattering experiment in order to succeed in the detection. Numerical examples confirmed the goodness of our results.
We also considered numerical examples with two or more faults. In this case linearity assumption is broken as mutual scattering takes place between the faults. However, once again the previous PD estimation worked well, not only for well separated faults (where mutual scattering can be considered negligible) but also for two adjacent faults. A general theory in order to establish the effect on multiple scattering on the performance achievable by the presented method can derived by following the approach adopted in . However, it is clear that this point deserves its own study and it is therefore beyond the scope of the present contribution and will be addressed in a future work.
This work was supported by Italian Ministry of University and Research through the FIRB initiative under the project MICENEA ( RBFR12A7CD).
References and links
1. C. Elachi, “Waves in active and passive periodic structures: a review,” IEEE Trans. Microwave Theory Tech. 24(12), 1666–1698 (1976).
2. J.D. Joanopulos, R.D. Meade, and J.N. Winn, Photonic Crystals: Modeling the Flow of Light (Princeton U. Press, 1995).
3. K. Yasumoto, Electromagnetic Theory and Applications for Photonic Crystals (Optical Engineering, CRC Press, 2005). [CrossRef]
4. C. Caloz and T. Itoh, Electromagnetic Metamaterials: Transmission Line Theory and Microwave Applications (Wiley-IEEE Press, 2005).
5. S. Enoch, B. Gralak, and G. Tayeb, “Enhanced emission with angular confinement from photonic crystals,” Appl. Phys. Lett. 81, 1588–1590 (2002). [CrossRef]
6. E.R. Brown, C.D. Parker, and E. Yablonovitch, “Radiation properties of a planar antenna on a photonic-crystal substrate,” J. Opt. Soc. Am. B 10, 404–407 (1993). [CrossRef]
8. G. Tayeb and D. Maystre, “Rigorous theoretical study of finite-size two-dimensional photonic crystals doped by microcavities,” J. Opt. Soc. Am. A 14(12), 3323–3332 (1997). [CrossRef]
9. A. Ludwig and Y. Leviatan, “Analysis of bandgap characteristics of two-dimensional periodic structures by using the source-model technique,” J. Opt. Soc. Am. A 20(8), 1553–1562 (2003). [CrossRef]
10. A. Ludwig and Y. Leviatan, “Analysis of arbitrary defects in photonic crystals by use of the source-model technique,” J. Opt. Soc. Am. A 21(7), 1334–1343 (2004). [CrossRef]
11. R. Solimene and G. Leone, “MUSIC algorithms for grid diasgnostics,” IEEE Geosc. Rem. Sens. Lett. 10(2), 226–230 (2013). [CrossRef]
12. A. Brancaccio, G. Leone, and R. Pierri, “Information content of Born scattered fields: results in the circular cylindrical case,” J. Opt. Soc. Am. A 15(7), 1909–1917 (1998). [CrossRef]
13. A. Brancaccio, G. Leone, and R. Solimene, “Fault detection in metallic grid scattering,” J. Opt. Soc. Am. A 28(12), 2588–2599 (1998). [CrossRef]
14. M. Bertero and P. Boccacci, Introduction to Inverse Problems in Imaging (Institute of Physics Publishing, 1998). [CrossRef]
15. J.P. Groby and D. Lesselier, “Localization and characterization of simple defects in finite-sized photonic crystals,” J. Opt. Soc. Am. A 25(1), 146–152 (2008). [CrossRef]
16. R.F. Harrington, Time Harmonic Electromagnetic Fields (Wiley and Sons, 2001). [CrossRef]
17. C. Balanis, Advanced Engineering Electromagnetics (Wiley and Sons, 1989).
18. R. Pierri, R. Solimene, A. Liseno, and J. Romano, “Linear distribution imaging of thin metallic cylinders under mutual scattering,” IEEE Trans. Antenn. Propag. 53 (9), 3019–3029 (2005). [CrossRef]