Abstract

The paper deals with the stochastic collocation analysis of a time domain response of a straight thin wire scatterer buried in a lossy half-space. The wire is excited either by a plane wave transmitted through the air-ground interface or by an equivalent current source representing direct lightning strike pulse. Transient current induced at the center of the wire, governed by corresponding Pocklington integrodifferential equation, is determined analytically. This antenna configuration suffers from uncertainties in various parameters, such as ground properties, wire dimensions, and position. The statistical processing of the results yields additional information, thus enabling more accurate and efficient analysis of buried wire configurations.

1. Introduction

The analysis of scattering properties of the perfectly conducting (PEC) objects buried in a lossy half-space has been of continuous interest in the various areas of application, such as ground penetrating radar (GPR) systems and related applications in civil engineering [1], as well as design and modeling of lightning protection systems (LPS) [2, 3]. Consequently, a deeper insight of scattering phenomena occurring in a lossy medium [4, 5] as well as the development of uncertainty analysis is required [6, 7].

Induced transient current distribution represents a fundamental quantity that can be used for determining other parameters of interest in the analysis of a scatterer (e.g., scattered voltage and transient input impedance) [8, 9].

The transient current induced along the wire scatterer, featuring a thin wire approximation, can be obtained using direct approach related to modeling of the electromagnetic (EM) wave coupling to a thin wire structure directly in time domain, providing better physical insight, accurate modeling of highly resonant structures, and easier implementation of nonlinearities [10]. This current is governed by integrodifferential equation of Pocklington type. An analytical solution of this equation can be obtained when dealing with canonical geometries, using carefully chosen set of approximations [11, 12]. The main advantage of analytical solution is the possibility of implementation for benchmark purposes, as well as for quick engineering estimation of the observed phenomenon [13]. Furthermore, analytical solution can be used within hybrid approaches for modeling complex structures, wherein computational time can be significantly reduced [9].

Variations in parameters can be explained by different causes, for example, by inability to obtain precise input parameters via measurement or changes in the obtained parameters due to variations in environmental conditions. The analysis of random variations regarding geometries and environment and/or materials of interest that may lead to significant misunderstanding or even errors in the analysis of buried objects is of interest [14]. Due to uncertain variations of parameters of interest, some techniques for an efficient integration of stochastic modeling have been developed [15]. The paper aims to demonstrate the ability of a precise analytical deterministic method to compute the currents induced on a straight buried wire, combined with an efficient and accurate stochastic method (SC) to integrate uncertainties with regard to parameters accuracy.

In Section 2, the fundamentals of the analytical solution of the space-time integrodifferential equation of Pocklington type are given for plane wave excitation as well as equivalent current source. The basis of the Stochastic Collocation method with the emphasis on its advantages in regard to other methods is given in Section 3. Computational results and discussion are outlined in Section 4, while the concluding remarks are given in Section 5.

2. Antenna Theory Formulation and Analytical Solution

2.1. Plane Wave Excitation

A horizontal thin wire scatterer of length and radius is buried in a lossy medium at depth . Properties of the medium are given with electrical permittivity and conductivity . The wire is illuminated by a transmitted part of a transient electromagnetic (EM) wave of normal incidence, as shown in Figure 1.

Time domain formulation for the transient analysis of horizontal straight wire buried in a lossy medium is based on the space-time Pocklington integrodifferential equation given with [16]where is the unknown space-time dependent current, is the tangential component of the transmitted field, and is the corresponding reflection coefficient arising from the Modified Image Theory (MIT) [17]. Detailed derivation of (1) can be found in [16].

The distance from the source point in the wire axis to the observation point located on the wire surface iswhile the distance from the source point on the image wire to the observation point on the original wire, according to the image theory, is

Time constant and propagation velocity in the lossy medium are given with

The influence of the earth-air interface is taken into account via the reflection coefficient arising from the MIT and is given with [17]where the corresponding time constants are

Note that the reflection coefficient (5) represents rather simple characterization of the earth-air interface, taking into account only medium properties. An accuracy of (5) has been discussed in [16].

Undertaking the analytical solution procedure documented in [16], the expression for the time dependent induced current for the case of impulse excitation can be written as follows:where function and coefficients and represent physical properties of the system, taking into account the dimensions of the wire and the distance from the interface:

Furthermore, other coefficients in (7) correspond to the properties of the medium and are given as follows:

Expression (7) represents the impulse response. Consequently, the response to an arbitrary excitation requires convolution. Only normal incidence is considered and the plane wave in the form of the double exponential pulse is assumed:

The transmitted electric field exciting the buried wire in the Laplace domain is given by [18]where represents Fresnel transmission coefficient defined by [16]

As the convolution, that is, the time domain counterpart of (11), would be too complex to be calculated analytically, the numerical convolution is carried out, as it is presented in [16].

2.2. Equivalent Current Source Excitation

A horizontal thin wire excited at one end with an equivalent current source (i.e., grounding electrode) of length and radius is buried in a lossy medium at depth . Properties of the medium are given as electrical permittivity and conductivity , as shown in Figure 2.

Governing equation for the unknown current induced along the wire is given in the form of time domain homogeneous Pocklington integrodifferential equation [19]:

The parameters in (13) are defined with expressions (2) through (6).

Integrodifferential equation (13) is solved analytically, using some carefully chosen approximations. As it has been shown in [19] time domain impulse response counterpart is obtained as follows:where corresponding coefficients and are the same as in (9), whereas is defined as

Equation (14) represents an analytical expression for the space-time distribution of the current induced along the thin wire excited by an equivalent current source in the form of the Dirac pulse.

To take into account realistic excitation of a lightning strike, one of the functions most frequently used is the double exponential pulse, given as [19]

Analytical convolution is undertaken with (14) and (16), to obtain the following expression for the current flowing along the electrode:

Equation (17) represents the expression for the space-time distribution of the current induced along the wire due to a double exponential pulse equivalent current source excitation.

3. Theoretical Principles and Statistical Methodology

First of all, the aim of this proposal is to focus on one particular aspect of uncertainty modeling. As detailed in [20], a huge difference exists between quantification of sources of uncertainty and Uncertainty Propagation (UP). In the following, we will lay emphasis on UP, assuming arbitrary variations of random inputs.

Among the huge diversity of statistical approaches available in the literature, the purpose of this paper is to focus on spectral stochastic techniques [6]. Relying on the physical problem under consideration, the statistics of the current induced in the center of the buried wire is expanded over an adapted function basis. The proposed methodology is relying on stochastic collocation methods as previously introduced in [21] for efficient computations of radar cross section from targets assuming uncertain geometrical inputs. A particular care was taken in [22] to optimization of SC methods.

In the following section, we will depict the proposed methodology according to an illustrative test case including one random variable (RV).

3.1. Theoretical Principles over 1-RV Polynomial Expansion

Without loss of generality, the theoretical principles are depicted considering only one random parameter (e.g., assuming uncertainty regarding soil conductivity and its impact on current ). In this framework, we may assume a given mapping through function and a RV given by an a priori probabilistic distribution ( standing for the relation between and ). In the following, we will establish the relations needed to access mean and variance of RV . First, an approximation of function is required over th order Lagrange polynomial basis:

Lagrange polynomials are given as follows:

One of the main advantages of using stochastic collocation method and polynomial expansion based upon these polynomials is to notice that , with standing for Kronecker symbol. Collocation points (so-called “sigma” points) correspond to the points given by Gauss quadrature rule attached to the probability distribution of random inputs (e.g., Legendre polynomials for a uniform distribution and Hermite polynomials for a Gaussian distribution). The choice of an adapted quadrature rule (different cubature points may be chosen [21]) allows for straightforward approach integral such as

In relation (19), is the probability density of random variable and a function with a sufficient regularity. The coefficients are called integration weights (or “collocation” weights). The weighted points are chosen to ensure an exact quadrature rule for polynomials with a degree lower or equal to . By definition, computation of the first statistical moment of a given random output considers integral terms. Thus, the expectation (mean) of (with a given RV) is given as follows:

Relying on relation (18), the output realization is replaced by its polynomial expression and expectation of RV may be written as

Then, the Gauss quadrature rule [23] provides the following relation:

Equation (24) is modified since every term (from sum) where does not match vanishes, and we obtain

By injecting relation (24) into (23), we deduce an intermediate expression involving expectation of RV :

Then, we substitute subscripts and as well as and in relation (18) to reach the crucial Kronecker property:

Finally, relation (25) is straightforward computed based upon weighted points :

As one may notice by decomposing the successive steps, computing the set of weighted points is quite simple because of the well-known property of Lagrange polynomials as the aforementioned. As it is the case for Galerkin-like methods [6], previous remarks entirely justify the choice of the same set of points for Gauss quadrature (computing efficiently the required integral) and collocation points (from Lagrange polynomials expansion of the given mapping). It is relatively easy to obtain the mean value from relation (27). Consequently, the same principle may be followed when computing variance of RV :

Similar to relation (18), we begin with the expansion of function over Lagrange polynomial basis as follows:

Based on Gauss quadrature rule from relation (20), the aforementioned relation may be written as

The replacement of relation (30) in (29) leads to

Similar to mean computing, most terms in this equation vanish (due to Kronecker property: i.e., nonzero terms are governed by ) which involves

By introducing (26), we may derive an approximated value of variance for RV , as follows:

In the following and in order to simplify the purpose, we define the statistical moments around one given output as follows, , where represents a given order of statistical moment (e.g., previous expectation of is defined as ). Thus, relations (27) and (33) may be rewritten using previous term as follows: , with such as and . In this framework, stands for a mapping (out of any consideration around polynomials used for expansion) dedicated to one statistical moment .

Thus, the generalization to multi-RVs cases may be obtained easily following the same strategy. We propose here to adapt previous notations by assuming random vector . Relying on previous notations for mapping , we may now compute th statistical moment of output such aswhere stands for weight of the expansion for random vector of size obtained from the generalization of relations (27) and (33) previously established. Relation (34) may imply different expansion orders ( representing the index of the random parameter ) among (see Section 4.3.3 for application). Some more details related to this approach are available in [23].

3.2. Multiple Independent Random Variables Principle

The fundamentals of SC technique [23] applied to the buried wire configuration taking into account three Random Variables (3-RVs) are outlined. The principle through which one operates with a random output (current) depending on random parameters (ground conductivity in mS/m, wire length , and depth in m) is illustrated. As depicted in relation (34), the SC technique is compatible with higher RV dimensions. The problem of interest requires one to model 3-RVs , , and (randomly modelled physical parameters), respectively , , and . The random variations related to independent , , and may be defined from the initial values , , and . Applying the same strategy as documented in [23], function is projected on a Lagrangian basis:where , , , and are the SC points required in corresponding random direction (i.e., accordingly to , , and parameters).

3.3. Computation of SC Statistical Moments: Advantages and Drawbacks

As previously stated, the computation of output statistics from relation (34) through a tensor product in each direction (i.e., for each RV) is rather simple. The SC technique gives the collocation sets of weighted points necessary to compute the needed statistics (e.g., mean and standard deviation). However, the technique also requires a particular attention regarding the cost/benefit ratio when increasing dimensions (tensor product of RV). Some methods for the improvement of stochastic techniques have been presented in [7, 14]. The part to follow proposes an alternative strategy to iteratively (i.e., increasing one RV at a time) and completely construct random model.

3.4. Statistical Convergence in relation to System Sensitivity

This part focuses on the SC requirements for variance convergence to improve the knowledge of the model sensitivity regarding each parameter independently. However, the interactions between RVs are not taken into account, but a view of the SC convergence and an idea of the RV global sensitivity are given.

First, it is compulsory to ensure the number of SC points necessary to assess convergence by computing statistics of the current.

There is no mathematical proof for SC convergence; despite all, based upon criterion proposed in [23], the SC convergence from SC method with weighted points for RV number is given bywhere is th order statistical moment computed from SC orders (i.e., involving +1 SC points according to relation (34), SC points are so-called “sigma” points). It can be noticed that the level of required convergence may be chosen arbitrarily (e.g., convergence criterion in relation (36) is lower than −1 dB). Criterion from relation (36) is defined considering, respectively, odd or even SC orders ( Indeed, even (resp., odd) number for implies (Gauss quadrature laws) taking into account (resp., or not taking into account) mean value of RV as a given “sigma” point. This may lead to a kind of “chessboard” effect when dealing with iterative convergence of the technique, and classically to distinguish SC convergence from even and odd number is well admitted. It is to be noted that the iterative process proposed in this work is quite similar to practice for the assessment of convergence and/or error Monte Carlo-like sampling methods. Indeed, assuming some given error/convergence threshold, it is possible to define stopping criterion (similar to previous relation) by computing confidence intervals (CIs). This often requires assumptions about the statistical distribution followed by the random output of interest, but some other techniques exist to infer trustworthy CIs. Some works were proposed in [24] to demonstrate the interest of bootstrapping techniques in this framework. An interesting idea may be to implement this method jointly with SC expansion to improve the assessment of CIs using, for instance, Bayesian Bootstrap [25].

SC method provides interesting global (i.e., relative to the entire bandwidth of study) information about the convergence and the number of SC points required. As a first approximation, the convergence of SC method informs about the continuity and/or discontinuity of the deterministic mapping under study (e.g., the current at the center of the buried electrode subject to random variations). Variance-based criteria are often used [23] to assess the sensitivity of random parameters in such modeling.

The variance computing to rank most influential parameters over the whole time simulation duration is proposed. Assuming a random model with three RVs (, , and ), the SC method may provide the experimental design (i.e., set of weighted points) required to compute the variance of a given output . In this framework, the criteria based upon cross-correlation between time variance obtained from the full tensorized SC model in relation (34) and the variances assuming one RV at a time (, , and ) may be helpful to provide global sensitivity indices of most influential random parameters given withwhere is the cross-correlation between the variances obtained from 1-RV model including RV and the full tensorized model including the whole RVs. As a consequence, the denominator of relation (37) stands for autocorrelation of the variance of the full tensorized model. Summing all terms over the observation time offers a global criterion regarding the entire time simulation. Practically, this criterion may be obtained without any additional computing costs since full tensorized model contains the results required to compute statistical moments regarding only one RV at a time.

4. Numerical Results from an Iterative Construction of the Random Model

Considering deterministic cases presented in Sections 2.1 and 2.2, the random model proposed in this paper will be validated in this section assuming different random parameters.

4.1. Definition of Test Cases

Different test cases that will be used are introduced. Table 1 summarizes three different configurations based upon the definition of parameters in Figures 1 and 2, including different kinds of uncertain parameters: materials (soil conductivity), geometries (burial depth and length of electrode), and sources parameters (lightning front times and time-to-half).

Data provided in Table 1 assumes uniform distribution of RVs. Without any loss of generality, it would be possible to apply SC strategy with different random distributions (e.g., normal, log-normal, and exponential).

The validation of the benefits that could be expected from SC modeling will be provided in the next section including different test cases, an overview of the SC convergence, and capability to predict the sensitivity of parameters. The advantage of SC will be highlighted in comparison with classical global sensitivity analysis techniques such as Sobol’s indices [26].

4.2. Equivalent Current Source Excitation

The purpose of this section is to focus on two kinds of uncertainties: environmental parameters such as the soil conductivity and the length of the electrode will be assumed as random; particular attention will be taken regarding random variations considering source (lightning) parameters.

4.2.1. SC Convergence and Stochastic Results for the Electrode and Soil Parameters (Test Case  (1), Table 1)

The random parameters ( and ) are modeled modifying the intensity of uncertainty with a uniform distribution (in mS/m) for soil conductivity and (in m) for electrode length.

First, it is necessary to define the number of SC points required to obtain converged results. The SC set of points (tensor product), regarding × , needed for soil conductivity and electrode length is constructed. Based upon previous work [27], far less sensitivity is expected from than from . A convergence criterion is defined in (36).

Figure 3 depicts the relative convergence gaps obtained considering different SC set of points. Except for (zero current), (red curve) is lower than −4 dB which validates the convergence of mean current obtained from 3 simulations. Increasing the SC approximation order from 5 to 9 points (resp., green, blue, and pink curves in Figure 3) shows that 5 points are necessary for soil conductivity to reach a converged result with lower than −2 dB. Similar to mean computation, Figure 4 offers a clear validation of the SC convergence regarding higher statistical moments (for instance, kurtosis) over the whole simulation duration: a set of 7 × 3 points is required to get a convergence lower than −2 dB.

With reference to Figures 3 and 4, a view of mean current and the statistical dispersion around it is proposed in Figure 5. The proposed “deterministic” results (black curves) are obtained from extreme values and show the time nonlinearity of the system. The converged results from SC technique (green curves) provide robust and trustworthy margins for computing mean current.

Figures 6 and 7 show the diversity of statistical results obtained from SC method (without any additional computing costs) with 5 × 3 points. The current variance gives a quick view of more sensitive time areas: for instance, there is no match between maxima of mean and variance, Figure 6. The converged data from Figure 7 is crucial to better assess statistical behavior of the system with only 7 × 3 points.

4.2.2. Stochastic Results for the Parameters of the Lightning (Test Case  (2), Table 1)

This part focuses on the importance of the equivalent current source parameters (representing lightning current and considering uncertain front time and time-to-half) as uncertain inputs [28].

Figure 8 was obtained from 32 + 52 + 72 = 83 simulations (results from SC order 2 not shown here with 3 × 3 points). Obviously, high rates of convergence are reached for current mean since results (green, Figure 8) from 5 × 5 points (i.e., 5 points for each RV) almost overlap data involving 7 × 7 SC points (blue, Figure 8). Current variance provides interesting information about the dispersion of data around mean values (via standard deviation in Figure 8). Thus, the standard deviation is about 20 mA at μs.

Figure 9 shows the statistical variations of average and mean + 2 standard deviation of the current in function of position along the electrode, at time instant μs. Similar to Figure 8, the convergence is reached with a restricted number of simulations. The statistical observation lays emphasis on the greater dispersion of the results in the first half of the electrode (zero-variance terminals). The standard deviation of the current is about 50 mA at the center of the electrode.

4.3. Plane Wave Excitation (Test Case  (3))

This section deals with some numerical results and statistics considering the current at the center of the wire excited by a transmitted plane wave. The entire stochastic modeling is based upon realistic values of environmental parameters (Table 1):(i)Soil conductivity : with and a zero-mean RV with a uniform distribution from 1 to 9 mS/m.(ii)Length of wire: with and a zero-mean RV with a uniform distribution from 9.5 to 10.5 m.(iii)Burial depth d: with and a zero-mean RV with a uniform distribution from 2.5 to 5.5 m.

Without loss of generality, the problem can be addressed following different assumptions about the statistical distribution laws [29].

4.3.1. Numerical Results for Entire Random Model: Fully Tensorized Modeling

Figure 10 shows mean (+ one standard deviation) of the current at the center of the wire under uncertain constraints fully tensorized (i.e., with 3 RVs). The main difficulty relies on the number of samples needed to assess converged statistics: 33 + 53 + 73 + 93 = 1224 points are required to ensure the asymmetrical convergence of 6th order (i.e., 7 points for each RV).

The sensitivity analysis provides relevant information needed to decrease the total number of SC points required for each RV and optimize the “full-tensor” random model to an “asymmetrical” one. To that end, and assuming the effects between random parameter are weak (which should be the case here regarding the physical nature of those parameters), SC may help to provide information about the relative importance of each parameter.

4.3.2. Sensitivity Analysis Assuming Low Levels of Interaction between Random Parameters

Figure 11 shows the results obtained by applying relation (37) for RV1, RV2, and RV3. As expected, fewer points are expected to ensure high convergence rate for RV2 () than for RV1 () and RV3 ().

Figure 12 emphasizes the importance of RV1 over other RVs: the relative dispersion from soil conductivity is larger than the one given by the length of the wire. An intermediate level is expected from RV2 (burial depth). In Figure 12, a quick overview of the 1-RV model sensitivity throughout the simulation time is given. This first-order trend is validated regarding the current variance computed from a full tensorized SC model ( points, black in Figure 12): in comparison with results obtained from 1-RV SC modeling and mostly over the whole simulation duration soil conductivity is the most influential parameter. From the beginning to 100 ns, burial depth seems to play a key role. Finally, as expected, the line length of the electrode does not play a major role. As previously explained, the variances (Figure 12) offer a first assessment of the influence of each random parameter. The SC method gives an overview of the potential relative influence of each random variable. This may be emphasized from the criterion given in relation (37) where , , and . This even more emphasizes the high similarity between variances given by (conductivity) relative to the full tensorized random model, following parameter (burial depth), and finally (length of the electrode) which is in accordance with previous statement for sensitivity analysis.

The advantage of SC comparatively to alternative classical global sensitivity analysis is to provide quick and accurate view of the influence of RVs at each time of the simulation. Thus, assuming low levels of interactions between random parameters, the first-order global sensitivity analyses from [26] are relevant to rank most influential parameters. In this framework, since it would be too demanding to compute previous criteria at each time of simulation, the focus is put on the maximum of current over the entire duration of simulation (max(I)).

Figure 13 gives an overview of the distribution of currents at the center of the electrode with random MC distributions of , , and from 0 to 1 μs. It is interesting to put the emphasis on the current amplitude since it plays a key role during the design of the lightning protection systems. In this framework and relying on Figure 13, Figure 14 shows the variability of max(I) throughout its empirical cumulative distribution function (CDF). The levels of variations are quite spread since max(I) is between 2 and 10 mA. Regarding time to maximum of current, Figure 15 offers also a synthetized view of its variability: it is between 80 ns and 220 ns. As depicted in Figures 10 and 12, the maximum statistical dispersion (i.e., maximum variance) is expected over this time interval [50 ns; 250 ns]. By computing classical first-order global sensitivity indices from [26] regarding max(I), it is possible to rank RVs from most to least influential: respectively, , , and . Thus, the results are in accordance with those previously obtained from relation (37) and SC simulations. It should be noticed that more than ten times speedup is obtained by SC approach relatively to classical calculation of Sobol’s indices. Indeed, simultaneously to the quantification of the statistical moments of current, the simulations with full tensorized SC models (including convergence testing for SC orders) were achieved in less than 1 hour, whereas global sensitivity analysis requires 8 hours of simulation (PC Intel Xeon 3.30 GHz, 16 Go RAM). Indeed, depending on the nature of the assumed random inputs (type: geometrical, electrical ones, for instance; king of uncertainty: distribution, range of variation…), the computing time required for one call to deterministic time domain model is varying between 5 and 7 seconds. It is to be noticed that SC method reveals to be competitive in comparison to MC simulations (required to compute Sobol’s indices) since a quasi-10x speedup is reached in that case. Moreover, from a sensitivity analysis point of view, the experimental design (so-called “sigma” points used for SC) provides not only a quantitative view of the impact of uncertainties around input parameters but also a trustworthy overview of the sensitivity of these parameters throughout variance-based method. Finally, the global sensitivity analysis as the one expected from MC-Sobol’s indices may be extended without any additional computing costs to a huge variety of outputs (not only restricted to max(I) here for Sobol’s indices).

4.3.3. 3-RV Full-Tensor Optimization (Asymmetrical SC)

In this section, the optimization of the symmetrical full tensorized SC experimental design is provided (in order to improve the efficiency of SC strategy).

Figure 16 provides convergence rates from the current variance including a complete random model: only 5 points are necessary to precisely describe the influence of random burying depth d (RV3). Nearly zero levels of the current (mean and variance) below 0.03 μs involve instability of the convergence criterion (and positive SC gaps). Finally, Figure 17 shows a good agreement between fully tensorized statistics of the current obtained with 343 points and the results given with 105 following previously depicted strategy.

5. Conclusion

The coupling of deterministic time domain analytical solutions of the Pocklington integrodifferential equation with stochastic collocation technique provides crucial information for the calculation of the response of wire configurations buried in lossy ground. The robustness, accuracy, and convergence of the two techniques (deterministic and stochastic) ensure useful statistics for designing various GPR or lightning protection systems by taking into account their intrinsic random characteristics (variations due to material parameters). The rapidity and the accuracy of the analytical expression offer interesting prospects to take benefit from SC method: robustness for time simulations (characterized by their high transient variations), quick convergence, and nonintrusiveness. Future work will be devoted to the analysis of benefit of the use of the space-time Pocklington equation coupled with proposed stochastic strategy in comparison with other costly sampling methods such as Monte Carlo.

Competing Interests

The authors declare that they have no competing interests.

Acknowledgments

This work benefited from networking activities carried out within the EU funded COST Action TU1208 “Civil Engineering Applications of Ground Penetrating Radar.”