1. Introduction

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 [1]. Photonic bandgap crystals [2, 3] and metamaterials [4] are examples that find application in high gain antennas [5, 6] and waveguides [7] 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 [11]. 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 [12]. 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 [13].

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) [14] 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 (P_D) and of false alarm (P_FA) 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.

It is necessary to remark that a similar problem has been already addressed in [15] and [13].

In [15] 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 [13], 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 P_D and P_FA 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 r_n 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 r_s. Therefore, we tackle the fault detection problem for a two dimensional scalar TM_z configuration where only the z directed electric field component is relevant.

 

Fig. 1 Scatterers location and numbering.

Download Full Size | PDF

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 E_s can be though as being radiated by the equivalent polarization current $J({\mathbf{r}}^{\prime })=j\frac{\beta }{\zeta }(1-{\epsilon }_r)E^{N-1}({\mathbf{r}}_s,{\mathbf{r}}^{\prime })$ just supported over Ω_m, so that

Substituting the equivalent source expression in Eq. (2) and assuming that the field E^N−1 does not significantly vary on the small section Ω_m, Eq. (2) can be rewritten as

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 [16]. Accordingly, the Green’s function of the full grid is given by

The expansion coefficients a_n,b_n 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 [17]

The elements of the N × N matrix ${\underset{¯}{\underset{¯}{A}}}^i(i=0,1)$ are given by:

Accordingly, by using the expression of the Green’s function in Eq. (5) the integral in Eq. (3) can be calculated as:

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 r_m due to a source located in r_s is given by:

Substituting Eqs. (11) and (10) into Eq. (3) we obtain

2.3. The scattering model

By using the same approach as followed in [13] 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:

To this end, the ”modified ” matrix ${\underset{¯}{\underset{¯}{A}}}^{(m)}=\underset{¯}{\underset{¯}{A}}-\mathrm{\Delta }{\underset{¯}{\underset{¯}{A}}}^{(m)}$ 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 $a_m^{(m)}=y_m$ instead of 0. However, it can be ignored, because as shown below, it does not affect the final result.

By comparing Eqs. (20) and (16) the following expression is obtained:

Since the only elements of matrix $\mathrm{\Delta }{\underset{¯}{\underset{¯}{A}}}^{(m)}$ different from zero belong the m-th row and column, it is straightforward to see that its product for $\underset{¯}{a}$ returns a vector whose elements are all proportional to a_m, 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:

Therefore, accordingly to Eq. (21), Eq. (23) can be rearranged 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 a_m and b_m:

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 N_O and N_S 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).

First, we note that starting from Eq. (3) the detection problem can be formulated in a completely equivalent way. In particular, rewrite Eq. (3) as

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 $\underset{¯}{\underset{¯}{L}}$ only requires the knowledge of the full grid coefficients b_m, $\widetilde{\underset{¯}{\underset{¯}{L}}}$ relies also on the fault location as the latter varies over all the grid. Accordingly, in order to build up $\underset{¯}{\underset{¯}{L}}$ the forward problem has to be solved only once (of course for each source/observation position). Instead, $\widetilde{\underset{¯}{\underset{¯}{L}}}$ also depends on E^N−1(r_n,r_s) and hence N² forward problems must be solved. It is clear that all these efforts can be accomplished off-line. Nonetheless, as N grows up the $\widetilde{\underset{¯}{\underset{¯}{L}}}$ 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 ${\underset{¯}{E}}_s\in C^{N_O{N}_S}$ one has to choose the ”right” $\underset{¯}{\gamma }$ belonging to the state space $S=\{{\underset{¯}{\gamma }}_1,{\underset{¯}{\gamma }}_2,\cdots ,{\underset{¯}{\gamma }}_{N^2}\}$, with ${\underset{¯}{\gamma }}_i$ 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 $C^{N^2}$ 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 $\underset{¯}{\underset{¯}{L}}$ 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 $\underset{¯}{\gamma }$, due to noise and to necessity to adopt a regularization scheme solving for $\underset{¯}{\gamma }$ 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 $\underset{¯}{\underset{¯}{L}}=\underset{¯}{\underset{¯}{U}}\cdot \underset{¯}{\underset{¯}{S}}\cdot {\underset{¯}{\underset{¯}{V}}}^H$, 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 N_T singular values.

More in detail, say $c_i={\underset{¯}{v}}_i^H\cdot \gamma$ the projections of the unknown vector onto the i-th singular vector (i.e., the i-th column of $\underset{¯}{\underset{¯}{V}}$), then under ideal case the reconstruction is provided by the series:

When data are affected by noise $\underset{¯}{\eta }$, 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 c_i values for i ≤ N_T); the second one is due to the additive noise present on the data. As the truncation index N_T increases, the first term tends to the ideal value but the second term increases as well. Accordingly, truncation index N_T 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 [14]. Therefore, while the following discussion holds true provided to replace the truncation index N_T with a suitable regularization parameters,TSVD is preferred for its simplicity.

Let us define the so-called spatial content [13] and side-lobe content. The spatial content is defined as

If the noise $\underset{¯}{\eta }$ is a zero mean Gaussian process whose variance is $N_0^2$, it results in the reconstruction as a zero mean Gaussian error with variance ${\mathrm{var}}_{\eta }(n)=N_0^2{\displaystyle \sum _{i=1}^{N_T}\frac{{|v_i(n)|}^2}{{\sigma }_i^2}}$. As a consequence, the reconstructed function is still a Gaussian distributed random process whose mean value $mea{n}_{\widetilde{\underset{¯}{\gamma }}}(n)$ 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 $\tilde{\gamma }(n)$ is obtained, a second threshold A_th, the detection threshold, is introduced and the position of the fault is recovered as the locus of points where $|\tilde{\gamma }(n)|\ge A_{th}$.

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 P_D and P_FA.

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.

 

Fig. 2 Grid of N = 121 scatterers. The observation points are taken along the dot circle. Stars denote the source positions. Cylinders are numbered form the upper side from left to right. Red and green points indicate position m = 30 and m = 61 respectively.

Download Full Size | PDF

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 N_T = 60 corresponds to cutting out singular values whose amplitude is 10dB lower than the maximum one, whereas N_T = 87 and N_T = 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.

 

Fig. 3 Normalized singular values for the configuration of Fig. 2.

Download Full Size | PDF

 

Fig. 4 Spatial content for N_T = 60 (left) and N_T = 87 (right).

Download Full Size | PDF

 

Fig. 5 Side lobes content for N_T = 60: central fault at m = 61 (left); peripheral fault m = 6 (right). Note that the gray scale is enhanced with respect to Fig. (4) due to the lower value of side lobes content.

Download Full Size | PDF

Next, let us study the behavior of P_D(n) and P_FA(m,n) with respect to the detection threshold level A_th for different values of the data Signal to Noise Ratio (defined as $SNR=\frac{{\Vert {\underset{¯}{E}}_s\Vert }^2/N_{data}}{N_0^2}$) and of the truncation index N_T. 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 P_FA achieves its maximum, so that the worst case is considered.

In Fig. 6 P_D and P_FA are plotted as function of A_th for different values of N_T. 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 A_th 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 A_th. 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.

 

Fig. 6 P_D(61) (solid lines) and [P_FA(n,61)]_max(n) (dashed lines) versus A_th for different values of the truncation index: SNR = 10dB (top); SNR = 5dB (bottom).

Download Full Size | PDF

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 N_T = 60 and N_T = 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.

 

Fig. 7 Results for a fault whose actual position is denoted by the red circle. White circles denote detected faults for SNR = 10, N_T = 60 and different A_th.

Download Full Size | PDF

 

Fig. 8 Results for a fault whose actual position is denoted by the red circle. White circles denote detected faults for SNR = 10, N_T = 112 and different A_th.

Download Full Size | PDF

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.

 

Fig. 9 Results for two faults whose actual position is denoted by the red circles. White circles denote detected faults for SNR = 10, N_T = 60 and A_th = 0.3.

Download Full Size | PDF

 

Fig. 10 Results for three (left) and five (right) faults whose actual position is denoted by the red circles. White circles denote detected faults for SNR = 10, N_T = 60 and A_th = 0.3.

Download Full Size | PDF

5. Conclusion

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 [18]. 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.