Inverse Problem for an Electrical Dipole and the Lightning Location Passive Monitoring System

Abstract

We solve the problem of the locating parameters, identifying equivalent dipole electromagnetic radiation source through measured horizontal magnetic and vertical electric components at some point of the infinite conducting ground. Methods based on analysis of measured signals are suggested. The problem under consideration, like any inverse problem of mathematical physics, is ill-conditioned. The consequences of this are the high sensitivity of the algorithm to the errors in the source data and calculation errors. All these circumstances do not allow to estimate the accuracy and reliability of the results obtained with the help of single-scale algorithms. The considered problem is contained in a complex of mathematical models of the practically important problem of forecasting the development of thunderstorm foci. Lightning meteorology focuses on investigating the lightning activities in different types of convective weather systems and the relationship of lightning to the dynamic and microphysical processes in thunderstorms. With the development and application of advanced lightning detection and location technologies, lightning meteorology has been developed into an important interdisciplinary between atmospheric electricity and meteorology. This paper reviews (1) methods to identify the dipole location and (2) possibilities to analyze the pre-radiation of thunderstorm clouds by the passive methods.

1. Introduction

The problem of identifying the position parameters of an arbitrary-oriented electric dipole over a plane with infinite conductivity from its electromagnetic field induced at the observation point is considered. The considered problem is contained in a complex of mathematical models of the practically important problem of forecasting the development of thunderstorm foci. To estimate the location of the electromagnetic radiation source (EMR) based on the results of a single-point observation of the electromagnetic field induced by it, up to 1990, a number of devices were developed, based on the use of a vertical dipole as a model of an EMR source and physically realizable analog algorithms. In the next decade, the development of the use of an arbitrarily oriented dipole as a model of the EMR source and digital processing of observed signals was developed [1, 2, 3].

The effective algorithm for the single-point distance determination to a pulsed EMR source was proposed and investigated in the papers [4, 5], and pointed out to the emergence of irremovable uncertainty in the location of the EMR source by a one-point method caused by the difference in the orientation of the dipole from the vertical. The use of two or more observation points that do not belong to the same straight line makes it possible to determine the position parameters of the equivalent dipole [4, 6, 7].

A review of the status of passive storm monitoring systems by the end of 2003 and the demonstration of the use of lightning-position detection systems for the passive radar of hazardous meteorological phenomena are presented in [1, 2]. Modern methods for analyzing the field, allowing to determine the parameters of the source of EMP, characterizing its location and orientation are presented in the works [4, 5, 8] within the ISTC project #1822 was developed. As a result of the conducted field tests of this sample, during May–August 2004 more than 2.5 million atmospherics were recorded. Of these, not more than 10% were classified as radiation from a lightning discharge, the rest were classified as pre-threat radiation. Thus, the proportion of inter-cloud and intra-cloud discharges relative to the cloud-ground discharges turned out to be much higher than those noted for the work [9].

The registration of ominous cloud radiation (i.e., before the first lightning flash) by a single-point lightning protection system, for the purposes of forecasting thunderstorm development, was provided by the expansion of the dynamic range of receiving equipment and the further development of mathematical and software. At present, active radar facilities using a comparison of sounding and reflected signals are used to analyze the pre-threat state [9, 10, 11]. The obvious drawbacks of this approach are: (1) the high cost and the presence of the human factor; (2) the low probability of detecting individual discharges; (3) the lack of ecological compatibility due to the application of microwave radiation and its effect on the clouds. Passive methods of analysis of the pre-threat state, based on the analysis of the intrinsic radiation of clouds, exclude the noted shortcomings of active methods.

Previously, in passive thunderstorm monitoring systems, pre-threat radiation was either filtered out, or it was worked out incorrectly by single-point systems of thunderstorm location [12]. The database of ISTC project #1822 [8] field trials provides a large experimental material for testing the adequacy of thunderstorm models and testing of the developed software and devices. In general, the results of the project demonstrated the possibility and necessity of creating a new generation of thunderstorm detection systems and expanding the range of tasks they solve. This leads not only to a significant revision of the requirements for their technical characteristics, but also to the development of new mathematical models and algorithms for analyzing thunderstorm phenomena, their tracing and display, as well as archiving and their use by specialists from different subject areas.

2. Statement of the dipole location problem

2.1. Direct problem

Born and Wolf [13] give the expressions for the electromagnetic field of a dipole ${\mathbf{P}}_0=p\left(t\right){\mathbf{n}}_0\delta \left({\mathbf{r}}_0\right)$ (moment p, direction n₀, location r₀, δ(*) is the Dirac delta function) in vacuum

$$\mathit{E}\left({\mathit{P}}_0\right)=\frac{1}{4\pi {\epsilon }_0}\left[\frac{3p\left(t\right)}{r^3}+\frac{3p^{\prime }\left(t\right)}{c{r}^2}+\frac{p^{\prime \prime }\left(t\right)}{c^{2}r}\right]\left({\mathit{n}}_0,{\mathit{e}}_0\right){\mathit{e}}_0-\frac{1}{4\pi {\epsilon }_0}\left[\frac{p\left(t\right)}{r^3}+\frac{p^{\prime }\left(t\right)}{c{r}^2}+\frac{p^{\prime \prime }\left(t\right)}{c^{2}r}\right]{\mathit{n}}_0,$$ E1

$$\mathit{H}\left({\mathit{P}}_0\right)=\frac{1}{4\pi }\left[\frac{p^{\prime }\left(t\right)}{r^2}+\frac{p^{\prime \prime }\left(t\right)}{r}\right]\left({\mathit{n}}_0\times {\mathit{e}}_0\right)$$ E2

here ${\mathbf{e}}_0$ is the direction from the dipole to the observer and the single and double primes denote time derivatives. For distances up to 300 km, a conducting plane surface is used to model the ground. The dipole in the half-space bounded by the infinitely conducting plane results in a field (E, H) consisting of the field (E(P₀), H(P₀)) of the source P₀, and the reflected field; the latter can be represented by the field (E(P₁), H(P₁)) of an imaginary dipole ${\mathbf{P}}_1=p\left(t\right){\mathbf{n}}_1\delta \left({\mathbf{r}}_1\right)$ , which is the mirror image of the dipole P₁ therefore we have $\mathbf{E}=\mathbf{E}\left({\mathbf{P}}_0\right)+\mathbf{E}\left({\mathbf{P}}_1\right),\mathbf{H}=\mathbf{H}\left({\mathbf{P}}_{\mathbf{0}}\right)+\mathbf{H}\left(P_1\right)$ (see Figure 1).

Figure 1.

Geometrical parameters of the mathematical problem.

Using the Cartesian system of coordinates with origin at the observation point О and the Oz axis being the normal of the bounding plane, we present the electric (Е) and magnetic (Н) components of the field in the coordinate form [4]:

$$E_x\left(t\right)=E_y\left(t\right)=H_z\left(t\right)=0,E_z\left(t\right)=\left[q\left(t\right)+q^{\prime }\left(t\right)\right]\frac{3u-2v}{\alpha }+q^{\prime \prime }\left(t\right)\frac{u}{{\alpha }^2},$$ E3

$$H_x\left(t\right)=\text{sin}\phi \left(q^{\prime }\left(t\right)\frac{1}{\alpha }+q^{\prime \prime }\left(t\right)\frac{1}{{\alpha }^2}\right),H_y\left(t\right)=\text{cos}\phi \left(q^{\prime }\left(t\right)\frac{1}{\alpha }+q^{\prime \prime }\left(t\right)\frac{1}{{\alpha }^2}\right),$$ E4

here $\alpha =c/r$ is the inverse of the time of wave propagation from the source to the observer and the variables $u,v,q\left(t\right)$ are determined from the equalities

$$w\text{sin}\phi =\text{sin}\theta \text{cos}{\theta }_0\text{sin}\psi -\text{cos}\theta \text{sin}{\theta }_0\text{sin}{\psi }_0,$$ E5

$$w\text{cos}\phi =\text{sin}\theta \text{cos}{\theta }_0\text{cos}\psi -\text{cos}\theta \text{sin}{\theta }_0\text{cos}{\psi }_0,$$ E6

$$u=\text{sin}\theta \text{cos}\left(\phi -\psi \right),v=\frac{\text{sin}\left(\phi -{\psi }_0\right)}{\text{sin}\theta \text{sin}\left(\psi -{\psi }_0\right)},q\left(t\right)=w{\alpha }^{2}p\left(t\right)/\left(2\pi c^{2}r\right).$$ E7

Other parameters are given in Figure 1.

2.2. Inverse problem

Components $E_z$ , $H_x$ , and $H_y$ can be changed with the help of the antenna system consisting of a vertical electrical antenna and a pair of mutually orthogonal frame magnetic antennas, and hence used as initial data for the task of evaluating coordinates $\left(r,\mathrm{\psi },\mathrm{\theta }\right)$ for the dipole position and $\left({\psi }_0,{\theta }_0\right)$ for the dipole orientation.

The inverse problem has a few specific aspects that differs it from the direct problem. If the source model is known, one often solves the inverse problem by searching the values of the parameters of the model, which give the best explanation of the measured data. Such a search is usually performed by an optimization algorithm which varies the values of the parameters until a necessary agreement between the measured and the model predicted fields is achieved. The search algorithms may get stuck in a local minimum and the solution will not be determined correctly. The use of a global optimization technique usually diminishes the problem of local minima but increase computational complexity.

For our case, it is possible to directly invert the values of the measured field into the values of the parameters describing the source location. In the present article, we use an explicit identification to infer the parameters of an electric dipole source situated over a perfectly conducting plane.

3. Calculation of the parameter $\mathit{\phi }$ . Exclusion of frames with interference of signals from different sources

Parameter φ, hereinafter called pseudo-bearing. For a vertical dipole (polar angle ${\theta }_0=0$ we have $\phi =\psi ,w=\text{sin}\theta ).$ For an inclined dipole, the pseudo-bearing $\phi$ is differs from actual direction $\psi$ , this difference depends on the orientation angles ${\theta }_0$ and ${\psi }_0$ . For a horizontal dipole, we have $\phi ={\psi }_0+\pi$ and $w=\text{cos}\theta$ , i.e., the pseudo-bearing is independent of the actual direction $\psi$ and is determined by only the dipole orientation ${\psi }_0$ . That is why the well-known direction finders [1, 2, 3] correctly record only cloud-ground discharges and give false estimates of bearing on inter-cloud and intra-cloud discharges. In fact, these methods evaluate pseudo-bearing $\phi$ .

To calculate the pseudo-bearing we use the following formulas:

$$\text{cos}2\phi =\frac{h_{\mathit{xx}}-h_{\mathit{yy}}}{\sqrt{{\left(h_{\mathit{xx}}-h_{\mathit{yy}}\right)}^2+4h_{\mathit{xy}}^2}},\text{sin}2\phi =\frac{h_{\mathit{xx}}-h_{\mathit{yy}}}{\sqrt{{\left(h_{\mathit{xx}}-h_{\mathit{yy}}\right)}^2+4h_{\mathit{xy}}^2}},$$ E8

$$\overline{\phi }=\text{atan}\sqrt{\frac{1-\text{cos}2\phi }{1+\text{cos}2\phi }}\text{sign}{h}_{\mathit{xy}},$$ E9

here $h_{\mathit{xx}}={\int }_{T_1}^{T_2}{\left[H_x\left(t\right)\right]}^2\mathit{dt},h_{\mathit{yy}}={\int }_{T_1}^{T_2}{\left[H_y\left(t\right)\right]}^2\mathit{dt},h_{\mathit{xy}}={\int }_{T_1}^{T_2}{H}_x\left(t\right)H_y\left(t\right)\mathit{dt}$ , [T₁,T₂] is signal observation interval.

Obviously, the signals observed in the antenna system can represent interference of signals from several sources. For the correct work of the algorithms described in the work, it is necessary to exclude frames with the presence of interference. Let us consider a possible way of detecting interference.

Interference leads to the appearance of a trend in the array of instantaneous error values

$$∆\left(\overline{\phi },t_i\right)=H_x\left(t_i\right)\mathit{cos}\overline{\phi }-H_y\left(t_i\right)\mathit{sin}\overline{\phi },t_i\in \left[T_1,T_2\right],i=0,1,2,\dots ,\frac{T_2-T_1}{∆t}=N.$$ E10

Really let be $H_x\left(t\right)=\text{sin}{\phi }_1{S}_1\left(t\right)+\text{sin}{\phi }_2{S}_2\left(t\right)+N_x\left(t\right),H_y\left(t\right)=\text{cos}{\phi }_1{S}_1\left(t\right)+\text{cos}{\phi }_2{S}_2\left(t\right)+N_y\left(t\right),$ N_x, N_y are uncorrelated white noise in the antenna channels having the same intensity N. Then $∆\left(\overline{\phi },t\right)=S_1\left(t\right)\text{sin}\left(\overline{\phi }-{\phi }_1\right)+S_2\left(t\right)\text{sin}\left(\overline{\phi }-{\phi }_2\right)+N_x\left(t\right)\text{sin}\phi -N_y\left(t\right)\text{cos}\phi .$

It is easy to see that for ${\phi }_1\ne {\phi }_2$ function $∆\left(\overline{\phi },t\right)$ for any $\overline{\phi }$ has a trend due to the difference in the functions $S_1\left(t\right),S_2\left(t\right)$ of the radiation sources. On the contrary, for $\overline{\phi }={\phi }_1={\phi }_2$ , the function $∆\left(\overline{\phi },t\right)$ has no such trend. Effect of interference make autocorrelation in function $∆\left(\overline{\phi },t\right)$ . Thus, the fact of autocorrelation in the error function $∆\left(\overline{\phi },t\right)$ is criterion for the presence of interference.

Durbin-Watson test is often used to determine the presence of autocorrelation. This criterion is usually used to establish the fact of the presence of an autocorrelation dependence of the first order in the error series, i.e., between its neighboring values $∆\left(\overline{\phi },t_i\right)$ and $∆\left(\overline{\phi },t_{i+1}\right)$ , $t_i,t_{i+1}\in \left[T_1,T_2\right],i=0,1,2,\dots ,\frac{T_2-T_1}{∆t}=N,∆t$ is sampling step. Usually the neighboring error values are related by a stronger dependence than other values. Therefore, the absence of an autocorrelation dependence of the first order makes it possible to state quite confidently that there are no autocorrelation relationships in the errors.

Durbin-Watson statistics is calculated by the following formula:

$$d=\frac{{\sum }_{i=0}^{N-1}{\left(∆\left(\overline{\phi },t_{i+1}\right)-∆\left(\overline{\phi },t_i\right)\right)}^2}{{\sum }_{i=0}^{N-1}{\left(∆\left(\overline{\phi },t_i\right)\right)}^2}.$$ E11

We have $0\le d\le 4$ . The values of $d=0$ and $d=4$ correspond to the cases when there is a strict positive or negative linear relationship between the shifted series $\left\{∆\left(\overline{\phi },t_i\right)\right\}$ and $\left\{∆\left(\overline{\phi },t_{i+1}\right)\right\}$ , respectively. If $d=2$ , then the series is independent.

The method has been tested based on data generated as a result of three-month field testing of the lightning detection system in 2004 within the framework of the ISTC project #1822 [8] and containing 796,112 atmospheric spheres suitable for computation [8]. After calculating the value of the statistics $d$ for the entire set of signals from the signals (over the whole implementation), obtained because of the field tests of the lightning detection system and analysis of the results, the intermediate point values were chosen:

$${\mathrm{d}}_1=0.4,{\mathrm{d}}_2=0.7,{\mathrm{d}}_3=3.3,{\mathrm{d}}_4=3.6.$$

If $d\in \left[0,d_1\right]\cup \left[d_4,4\right]$ , then the hypothesis of the existence of first order autocorrelation in a series $\left\{∆\left(\overline{\phi },t_i\right)\right\}$ and $\left\{∆\left(\overline{\phi },t_{i+1}\right)\right\}$ is assumed. If $d\in \left[d_2,d_3\right]$ , then this hypothesis is rejected. If the value of the statistics $d$ is in other intervals, then an unambiguous answer about the existence of autocorrelation cannot be given.

The analysis of the detected signals using the Durbin-Watson test showed the following results (see Table 1).

Value Range of Durbin-Watson statistics	Number of signals
[0, 0.4]	572,005
[0.4, 0.7]	79,412
[0.7, 3.3]	144,695
[3.3, 4]	0

Table 1.

The relationship between the Durbin-Watson statistics interval and the number of corresponding signals.

Thus, we can confidently determine the coordinates of only 144,695 of the 796,112, based on the data of the entire implementation of each of the signals.

To eliminate this influence, you can select a part of the signal length (the implementation interval) at which interference appears less, or where there is none (there is no trend). Of course, the calculation of discharge characteristics over this implementation interval will be more accurate. In addition, the approach to isolating a “clean” signal can help in filtering the transients at the beginning of the signal reception.

4. Calculation of the parameters $\mathit{u}\mathit{,}\mathit{v}\mathit{,}\mathit{\alpha }$

Further, we assume the parameter φ to be computed, and interference is absent in the recorded signals $E_z$ , $H_x$ , and $H_y$ . Let us further use the function $H\left(t\right)=H_x\left(t\right)\text{sin}\phi +H_y\left(t\right)\text{cos}\phi$ .

4.1. Primal algorithm

An analytical solution of the problem to identify parameters $u,v,\alpha$ is given in [4]

$$\left\{\begin{array}{c}\alpha =\sqrt{\left(h_1{e}_2-h_2{e}_1\right)/\left(h_2{e}_0-h_1{e}_1\right)},\\ u=\left(e_0{e}_2-{e_1}^2\right)/\left(h_2{e}_0-h_1{e}_1\right),\\ v=\left(3h_{0}u\alpha -e_0\alpha -g\right)\left(2h_0\alpha \right),\end{array}\right.$$ E12

here

$$e_k={\int }_0^{\infty }{E}_z^{\left(k\right)}\left(t\right)H^{\left(k\right)}\left(t\right)\mathit{dt},h_k={\int }_0^{\infty }{\left(H^{\left(k\right)}\left(t\right)\right)}^2\mathit{dt},g={\int }_0^{\infty }{E}_z^{\prime }\left(t\right)H\left(t\right)\mathit{dt},k=0,1,2.$$ E13

Calculated parameters $u,v,\alpha$ , φ can be used to determine the distance to the dipole $r=c/\alpha$ , as well as guaranteed interval estimates and probabilistic characteristics of angular coordinates of the dipole position.

Primal algorithm calculates parameters using formulas (3). It is assumed that the signals are represented by counts of instantaneous values with a given sampling rate. To calculate the coefficients (4), we use the known quadrature formulas.

An analysis of the sensitivity of the direct algorithm to the noise in the measured signals is given in [5]. It is shown that to achieve acceptable accuracy, pre-processing of input signals is necessary. It is found the dependence of mathematical expectation and variance of the distance determination error on the source spectrum and transfer function of the input filter. The estimates for the case corresponding to the averaged source of lightning discharges indicate the possibility of using the algorithm under study in systems for locating the lightning centers and allow us to develop the requirements for hardware of such systems.

4.2. Extremal algorithm

The primal algorithm was obtained in the assumption that the radiation source is the dipole, and the underlying surface—a plane with the infinite conductivity. Due to non-ideal models and noises in observed signals, we are forced to allow discrepancies in the resultant equations. Hence, intuition shows that most preferable is the least-squares principle. According to this principle, parameters of the model are determined from

$$I\left(\mathit{x}\right)=\underset{T_1}{\overset{T_2}{\int }}{\left[x_1{E}_z^{\prime \prime }\left(t\right)+x_2{E_z}^{\prime }\left(t\right)+x_3{H}^{\prime \prime }\left(t\right)+x_4{H}^{\prime }\left(t\right)+x_{5}H\left(t\right)\right]}^2\mathit{dt}=\left(\mathit{x},\mathit{Ax}\right)\to {\text{max}}_{\mathit{x}}$$ E14

$$\mathbf{A}={\int }_{T_1}^{T_2}\mathbf{S}{\mathbf{S}}^T\mathit{dt},\mathbf{S}={\left(k{E}_z^{\prime \prime }\left(t\right){E_z}^{\prime }\left(t\right){\mathit{kH}}^{\prime \prime }\left(t\right)\begin{array}{ccc}{H}^{\prime }\left(t\right)& k^{-1}H\left(t\right)& \end{array}\right)}^T,$$

k is the time constant for differentiating and integrating sections. This constant is introduced to make scales of coordinates for vector x consistent.

The extreme algorithm [14] for determining parameters u, ν, and α is calculated as follows:

Step 1. Let ${\lambda }_0\le {\lambda }_1\le {\lambda }_2\le {\lambda }_3\le 4$ be eigenvalues of matrix А, ${\mathbf{x}}^{\left(i\right)},i=0,1,2,3,4$ are corresponding eigenvectors.

Step 2. Let ${∆}_{\mathit{mn}}$ be the determinant of the matrix derived from the matrix $\left(\begin{array}{ccc}{\mathbf{x}}^{\left(2\right)}& {\mathbf{x}}^{\left(3\right)}& {\mathbf{x}}^4\end{array}\right)$ through deleting lines having numbers m and n.

Step 3. For $i=1,2$ do

$${\beta }_i=\left\{\begin{array}{cc}\frac{2{∆}_{25}}{-{∆}_{24}-{∆}_{15}+{\left(-1\right)}^i\sqrt{{\left({∆}_{24}+{∆}_{15}\right)}^2-4{∆}_{14}{∆}_{25}}}& \text{if}\left|{∆}_{25}\right|\ge \left|{∆}_{14}\right|,\\ \frac{-{∆}_{24}-{∆}_{15}-{\left(-1\right)}^i\sqrt{{\left({∆}_{24}+{∆}_{15}\right)}^2-4{∆}_{14}{∆}_{25}}}{2{∆}_{14}}& \text{otherwise},\end{array}\right.$$

$$u_i=\frac{{\beta }_i{∆}_{13}+{∆}_{23}}{{∆}_{12}},w_i=\frac{-{∆}_{45}}{{\beta }_i{∆}_{24}+{∆}_{25}}$$

.

Step 4. Let $i^{\ast }=\text{arg}\underset{i=1,2}{\text{min}}\left[\left(1,{\beta }_i,-u_i,-{\beta }_i{w}_i,-{\beta }_i^2{w}_i^2\right)\mathbf{A}{\left(1,{\beta }_i,-u_i,-{\beta }_i{w}_i,-{\beta }_i^2{w}_i^2\right)}^T\right]$ ,

Step 5. Let $\alpha =\frac{{\beta }_{i^{\ast }}}{k},u=u_{i^{\ast }},v=\frac{w_{i^{\ast }}-3u}{2}.$

A calculation experiment was performed to analyze accuracy and stability of the extreme algorithm, as well as to compare its characteristics with those of the primal algorithm. Software of the algorithm is implemented with the observation of agreements for the instrumental environment developed under the ISTC project 1822 [8]. Signals $E_z\left(t\right),H\left(t\right)$ are presented as integral counts of rapid values in discrete time moments, step of time discretization is $∆t=2\mathrm{\mu s}$ , number of quantification levels is $2^{16}$ . Input signals $E_z\left(t\right),H\left(t\right)$ are preliminarily processed by a smoothing filter. Usual difference and quadrature formulas with the accuracy of $O\left(\mathrm{\Delta }{t}^2\right)$ were used for numerical differentiation and integration. The Givens rotation method is used to solve the complete symmetric problem of eigenvalues.

Extreme algorithm gives better accuracy if compared with the primal algorithms. Most significant difference in characteristics is observed for distances less than 50 km: error of the primal method sharply increases with the distance decrease, and the error of the extreme method is decreasing. Since with the decreasing distance, one has increased error in presentation of the real source as an electrical dipole, the experimental results demonstrate the extreme algorithm to be expedient at the distance $r\le 50$ km.

4.3. Parametrization of the algorithms

To solve the problem to identify the location parameters because of its poor condition is proposed to use a parameterized set of algorithms, and the final decision to accept the results of statistical analysis. Parameterizations of the primal and extremal methods are considered in paper [15].

The spectrum of dipole moment $p\left(t\right)$ lies in a rather narrow frequency range, so the use of a bandpass filter with suitable lower and upper cutoff frequencies will maximize the use of all useful information that the signal carries and weaken the effect of noise present in the signal. Since the position parameters and the spectrum of the radiation source are not known a set of estimates $\tilde{E}=\left\{\left({\alpha }_i,u_i,v_i\right):i=1,2,\dots ,N\right\}$ for a family of bandpass filters with an amplitude-frequency characteristic

$$F_i\left(\omega \right)=\left\{\begin{array}{cc}1,& \text{if}\omega \in \left[{\widehat{\Omega }}_i,{\stackrel{\mathrm{ˇ}}{\Omega }}_i\right],\\ 0,& \text{otherwise},\end{array}i=1,2,\dots ,N\right.$$ E15

is constructed. The elements of the set $\tilde{E}$ are considered as an implementation of a vector random variable, and for estimating the true position parameters $u,v,$ and $\alpha$ to apply methods of statistical robust estimation.

It would be well to exploit more methods for increasing of statistical significance of result estimation. Put computational experiments’ results confirms efficiency the approach, but а wide variety of the measurement results do not lets to form up the well-provided measuring estimation of the dipole location parameters adequately. It is suggested that the filtration of the measurement results based on the cleaner functionals and combination of them for improving of the measuring estimation of the dipole location parameters. For construction of the cleaner functionals is the projection method [16]. This method is based on the projection of the inverse image of the vector-valued function $\left(E_z\left(t\right),H\left(t\right)\right)$ onto the linear manifold $\mathcal{L}$ of solutions of the system of differential Eqs. (1) and (2). There is a solution q(t) satisfying the system of differential Eqs. (1) and (2) under error-free values of the estimated parameters $u,v,\alpha$ . Therefore, the inverse image of the pair of signals $\left(E_z\left(t\right),H\left(t\right)\right)$ belongs to $\mathcal{L}$ , and the length of the projection is maximal. This allows you for filtering to use the length of the projection as a utility function. In terms of Fourier transforms, the utility function has the form [16]

$$U\left(u,w,\alpha \right)=‖\frac{E\left(\omega \right)\left(-u{\omega }^2+\mathit{\text{jωwα}}+w{\alpha }^2\right)+H\left(\omega \right)\left(-{\omega }^2+\mathit{j\omega \alpha }\right)}{{\left|-u{\omega }^2+\mathit{\text{jωwα}}+w{\alpha }^2\right|}^2+{\left|-{\omega }^2+\mathit{j\omega \alpha }\right|}^2}‖.$$ E16

Function $U\left(u,w,\alpha \right)$ is continuous for $\left(u,w,\alpha \right)\ne 0$ . Algorithm for calculating $U\left(u,w,\alpha \right)$ is completely stable for signals represented of 4096 readouts with the sampling time 2–10⁻⁶ s, i.e., containing 2047 harmonics with the first-harmonic discretization frequency 122 Hz.

Set $E\left(U_0\right)=\left\{\left(u,w,\alpha \right)\in \tilde{E}:U\left(u,w,\alpha \right)\ge U_0\}\right)$ contains most reliable estimates. Statistical analysis of the set of estimates makes it possible to estimate $\left(u,v,\alpha \right)$ probability density.

4.4. Spectral statistical method for location parameters identifying of a dipole electromagnetic radiation source

The spectral statistical method makes it possible to obtain more stable solutions at a lower computation cost compared with the previously developed parametric extremum method. The spectral statistical method algorithm can be naturally parallelized. The proposed method based on the analysis of the measured-signal spectra, allows one to get many estimates of the source location, choose the final estimate of the results of analysis of the entire totality of these estimates, and therefore, reach stability in determining the source location. To reduce the computational complexity, it is preferable to analyze individual harmonics rather than the entire frequency band. In this case, the algorithms for estimating the dipole location are substantially simplified [17].

As before, we will consider signals represented as a fast Fourier transformation $E\left(\omega \right),H\left(\omega \right)$ containing 2047 harmonics with the first-harmonic discretization frequency 122 Hz. Let us divide the set $\mathrm{\Omega }=\left\{{\omega }_1,{\omega }_2,\dots ,{\omega }_{2047}\right\}$ of all harmonics into subsets

$${\mathrm{\Omega }}_1=\left\{{\omega }_1,{\omega }_2,\dots ,{\omega }_{131}\right\},{\mathrm{\Omega }}_2=\left\{{\omega }_{132},{\omega }_{133},\dots ,{\omega }_{400}\right\}$$ and $${\mathrm{\Omega }}_3=\left\{{\omega }_{401},{\omega }_{402},\dots ,{\omega }_{2047}\right\}.$$

Spectral statistical algorithm to find estimates $\widehat{u},\widehat{v},\widehat{\alpha }$ of the parameters $u,v,\alpha$ and their variance ${\sigma }_{\widehat{u}}^2,{\sigma }_{\widehat{v}}^2,{\sigma }_{\widehat{\alpha }}^2$ is calculated as follows:

Step 1. For each $\omega \in \mathrm{\Omega }$ calculate

$$f\left(\omega \right)=\text{Re}\left[E\left(\omega \right)\overline{H\left(\omega \right)}\right],g\left(\omega \right)=\text{Im}\left[E\left(\omega \right)\overline{H\left(\omega \right)}\right],h\left(\omega \right)={\left|H\left(\omega \right)\right|}^2.$$

Step 2. Calculate the noise intensity

$${\sigma }_N^2=\frac{\sum _{\omega \in {\mathrm{\Omega }}_3}{\left|H\left(\omega \right)\right|}^2}{\left|{\mathrm{\Omega }}_3\right|}.$$

Step 3. For each $\omega \in {\mathrm{\Omega }}_2$ calculate

$$\widehat{\mathrm{u}\left(\mathrm{\omega }\right)}=\left(1+\frac{{\mathrm{\sigma }}_{\mathrm{N}}^2}{\mathrm{h}\left(\mathrm{\omega }\right)}\right)\frac{\mathrm{f}\left(\mathrm{\omega }\right)}{\mathrm{h}\left(\mathrm{\omega }\right)}.$$

Step 4. Calculate

$$\widehat{u}=\frac{\sum _{\omega \in {\mathrm{\Omega }}_2}\left(h\left(\omega \right)\widehat{u\left(\omega \right)}\right)}{\sum _{\omega \in {\mathrm{\Omega }}_2}h\left(\omega \right)},{\sigma }_{\widehat{u}}^2=\frac{\sum _{\omega \in {\mathrm{\Omega }}_2}\left(f\left(\omega \right)-h\left(\omega \right)\widehat{u}\right)}{\sum _{\omega \in {\mathrm{\Omega }}_2}h\left(\omega \right)}.$$

Step 5. For each $\omega \in {\mathrm{\Omega }}_1$ calculate

$$A\left(\omega \right)=\omega \sqrt{\frac{h\left(\omega \right)\widehat{u}}{f\left(\omega \right)}-1},S\left(\omega \right)=\frac{{\sigma }_N^2}{h\left(\omega \right)},\widehat{\alpha \left(\omega \right)}=\sqrt{\frac{A^2\left(\omega \right)-S\left(\omega \right){\omega }^2}{S\left(\omega \right)+0.5+\sqrt{\frac{S\left(\omega \right)A^2\left(\omega \right)}{{\omega }^2}+0.25}}}.$$

Step 6. Calculate

$$\widehat{\alpha }=\frac{{\sum }_{\omega \in {\mathrm{\Omega }}_1}\left[\frac{\sqrt{h\left(\omega \right)}}{{\omega }^2}A\left(\omega \right)\widehat{\alpha \left(\omega \right)}\right]}{\sum _{\omega \in {\mathrm{\Omega }}_1}\left[\frac{\sqrt{h\left(\omega \right)}}{{\omega }^2}A\left(\omega \right)\right]},{\sigma }_{\widehat{\alpha }}^2=\frac{{\sum }_{\omega \in {\mathrm{\Omega }}_1}\left[\frac{\sqrt{h\left(\omega \right)}}{{\omega }^2}A\left(\omega \right){\left(\widehat{\alpha }-\widehat{\alpha \left(\omega \right)}\right)}^2\right]}{\sum _{\omega \in {\Omega }_1}\left[\frac{\sqrt{h\left(\omega \right)}}{{\omega }^2}A\left(\omega \right)\right]}.$$

Step 7. For each $\omega \in {\mathrm{\Omega }}_1$ calculate

$$W\left(\omega \right)=\frac{f\left(\omega \right)A\left(\omega \right)-\mathit{\omega g}\left(\omega \right)}{h\left(\omega \right)A\left(\omega \right)}.$$

Step 8. Calculate

$$\widehat{w}=\frac{{\sum }_{\omega \in {\mathrm{\Omega }}_1}\left[\sqrt{h\left(\omega \right)}W\left(\omega \right)\right]}{\sum _{\omega \in {\mathrm{\Omega }}_1}\sqrt{h\left(\omega \right)}},{\sigma }_{\widehat{w}}^2=\frac{{\sum }_{\omega \in {\mathrm{\Omega }}_1}\left[\sqrt{h\left(\omega \right)}{\left[W\left(\omega \right)-\widehat{v}\right]}^2\right]}{\sum _{\omega \in {\mathrm{\Omega }}_1}\sqrt{h\left(\omega \right)}}.$$

The rationale for the effectiveness of the algorithm and computational experiments are given in [17].

5. Secondary processing of electromagnetic field monitoring results

The most important problem in work of the single-point systems of storm location is ineradicable error in the bearing determination of each source of radiation taken separately. Amount of this error depends on orientation of the lightning equivalent dipole, and the “pseudo-bearing” determined by the device for horizontal dipoles does not depend from the real bearing at all and is defined only by projection of orientation of the dipoles to the ground plane.

The given problem can be solved only by the analysis of the information from all set of discharges of the storm. At International conference ETC’2006 [18], it was suggested to display the stream of signals as a map of lightning discharge hit probability density in this or that point of terrestrial surface. The representation proposed increases probability of definition of true thunderstorm discharges location, and hence, raises accuracy of their fixing.

Let us describe the technique of this approach realization, the difficulties arising, and the variants of their decision in detail.

Let values ( $r\mathit{,}\mathit{\psi }\mathit{,}\mathit{\theta })$ be the polar coordinates of the lightning equivalent dipole for the measuring device coordinate system. As it is known [4, 18], Cartesian coordinates $x,y$ of the dipole meet the condition

$$x\text{cos}\phi +y\text{sin}\phi =r\text{sin}\theta \text{cos}\left(\phi -\psi \right)=\mathit{ru}.$$ E17

This equation includes known parameters: $r$ is distance from the dipole to the observation point, $u,v$ are the variables of the source dipole model (3)—and unknown Cartesian coordinates $x,y$ of the dipole location that is a feature of an uncertainty. The coordinates of the points of possible location of a dipole most remote from origin are the solutions of a set of equations

$$x\text{cos}\phi +y\text{sin}\phi =\mathit{ru},x^2+y^2=r^2.$$ E18

They are equal

$$x_{1,2}=r\left(u\text{cos}\phi \mp \sqrt{1-u^2}\text{sin}\phi \right),y_{1,2}=r\left(u\text{sin}\phi \pm \sqrt{1-u^2}\text{cos}\phi \right).$$ E19

The coordinates $x,y$ of all probable points of a source location introduce on a plane $\mathit{XOY}$ a section of a straight line (10) between points $\left(x_1,y_1\right)$ and $\left(x_2,y_2\right)$ . The set of probable location of a dipole have one-parameter representation

$$\mathcal{L}\left(r,u,v\right)=\left\{\begin{array}{c}x\left(\lambda \right)=\mathit{ru}\text{cos}\phi +r\left(2\lambda -1\right)\sqrt{1-u^2}\text{sin}\phi \\ y\left(\lambda \right)=\mathit{ru}\text{sin}\phi +r\left(1-2\lambda \right)\sqrt{1-u^2}\text{cos}\phi ,\\ z\left(\lambda \right)=\sqrt{r^2-x^2\left(\lambda \right)-y^2\left(\lambda \right)}=2r\sqrt{\lambda \left(1-\lambda \right)\left(1-u^2\right)}\end{array}|0\le \lambda \le 1\right\}.$$ E20

A measure of the coordinates set $\mathcal{L}\left(r,u,v\right)$ equal

$$\mathrm{L}=\underset{0}{\overset{1}{\int }}\sqrt{{\left(\frac{\mathrm{dx}\left(\mathrm{\lambda }\right)}{\mathrm{d\lambda }}\right)}^2+{\left(\frac{\mathrm{dy}\left(\mathrm{\lambda }\right)}{\mathrm{d\lambda }}\right)}^2+{\left(\frac{\mathrm{dz}\left(\mathrm{\lambda }\right)}{\mathrm{d\lambda }}\right)}^2}\mathrm{d\lambda }=4.844\mathrm{r}\sqrt{1-{\mathrm{u}}^2}.$$

The value of L represents a measure of the uncertainty in the estimation of the Cartesian coordinates of the dipole location. Uncertainty is absent when L = 0. It is reached for a vertical dipole (u = ±1). It is possible value u = 0 for a horizontal dipole when L reaches the maximum value.

A more practical measure of uncertainty is probability

$$\mathbf{P}\left[i,j,k\right]\left\{\right[(x\left(\lambda \right),y\left(\lambda \right),z\left(\lambda \right)]\in C\left[i,j,k\right]\}$$

of the dipole location in space cell

$$C_{\left[i,j,k\right]}=\left\{\left(x,y,z\right)|\begin{array}{c}\mathit{i\Delta }\le x\le \left(i+1\right)\Delta \\ \mathit{j\Delta }\le y\le \left(j+1\right)\Delta \\ \mathit{k\Delta }\le z\le \left(k+1\right)\Delta \end{array}\right\},$$

here $\mathrm{\Delta }$ is cell size.

It is obvious that ${\mathbf{P}}_{\left[i,j,k\right]}\mathit{>}\mathbf{0}$ if $C_{\left[i,j,k\right]}\cap \mathcal{L}\left(r,u,v\right)\ne \varnothing$ . To find cells in which the dipole is located with a nonzero probability, it is sufficient to find the minimal covering $C_{\mathit{XOY}}\left(r,u,v\right)$ of the projection of $\mathcal{L}\left(r,u,v\right)$ onto the plane XOY by flat cells $C_{\left[i,j\right]}:$

$$⌊\text{min}\left\{\frac{ru\text{cos}\left(\phi -\text{arccos}u\right)}{\mathrm{\Delta }},\frac{ru\text{cos}\left(\phi +\text{arccos}u\right)}{\mathrm{\Delta }}\right\}⌋\le i\le ⌈\text{max}\left\{\frac{ru\text{cos}\left(\phi -\text{arccos}u\right)}{\mathrm{\Delta }},\frac{ru\text{cos}\left(\phi +\text{arccos}u\right)}{\mathrm{\Delta }}\right\}⌉,$$

$$⌊\text{min}\left\{\frac{ru\text{cos}\left(\phi -\text{arcsin}u\right)}{\mathrm{\Delta }},\frac{ru\text{cos}\left(\phi +\text{arcsin}u\right)}{\mathrm{\Delta }}\right\}⌋\le j\le ⌈\text{max}\left\{\frac{ru\text{cos}\left(\phi -\text{arcsin}u\right)}{\mathrm{\Delta }},\frac{ru\text{cos}\left(\phi +\text{arcsin}u\right)}{\mathrm{\Delta }}\right\}⌉.$$

An example of coverage $C_{\mathit{XOY}}\left(r,u,v\right)$ is shown in the Figure 2.

Figure 2.

An example of coverage $C_{\mathit{XOY}}\left(r,u,v.\right)$

Cell $C_{\left[i,j\right]}\in C_{\mathit{XOY}}\left(r,u,v\right)$ in accordance with equality (12) uniquely defines a space cell $C_{\left[i,j,k\right]},k=\left[\sqrt{\frac{r}{\mathrm{\Delta }}-i^2-j^2}\right]$ with a nonzero probability of the dipole location.

To find this probability, we accept a completely plausible hypothesis about the uniform distribution of the probabilities of all possible values of the angle ${\psi }_0$ (see Figure 1). In this case the probability density for angle $\psi$ is [4]

$$f_{\psi }\left(\eta \right)=\frac{s}{2\widehat{\psi }\left[s^2{\text{cos}}^2\left(\phi -\eta \right)+{\text{sin}}^2\left(\phi -\eta \right)\right]},-\text{arccos}u\le \eta -\phi \le \text{arccos}u$$ E21

here

$$s=\left|\frac{1-\mathit{uv}}{\mathit{uv}}\right|,\widehat{\psi }=\text{arctan}\left(\frac{1}{s}\text{tan}\text{arccos}u\right).$$

Thus, for $C_{\left[i,j\right]}\in C_{\mathit{XOY}}\left(r,u,v\right)$ one can take

$${\mathbf{P}}_{\left[i,j,k\right]}=f_{\psi }\left(\eta \right),k=\left[\sqrt{\frac{r}{\Delta }-i^2-j^2}\right],\eta =\text{arctan}\frac{j}{i}.$$

Let $R\left(T\right)$ be the set of registered signals in the time interval $T.$ The integral estimate grade of radiation locus membership to cell $C_{\left[i,j,k\right]}$ over period $T$ let be

$$\text{grade}\left[i,j,k\right]=\frac{\sum _{d\in R\left(T\right)}{\mathbf{P}}_{\left[i,j,k\right]}\left(d\right)}{\underset{\left[l,m,n\right]}{\text{max}}\sum _{d\in R\left(T\right)}{\mathbf{P}}_{\left[l,m,n\right]}\left(d\right)}.$$

It is suggested to display the stream of signals as a $\text{grade}$ map. The representation proposed increases in probability of definition of true thunderstorm discharges location, and hence, raises accuracy of their fixing.

6. Using a network of sensors

Currently, mainly for monitoring thunderstorm activity, multi-point systems based on the application of the difference-ranging method and its variations to the results of monitoring the Earth’s electromagnetic field in the VLF and VHF bands are used. Some lacks such an approach are: (1) insufficient reliability, which results from the need communication networks for the exchange of information between observation stations spaced over significant distances, and the need for time synchronization at the sub-microsecond level; (2) the inability to analyze the cloud radiation before the thunderstorm. Previously, this was the reason for the development of single-point systems.

Now, the rapid development of computer technologies and communication systems have made it possible to combine single-point lightning location finders into a system and to organize joint processing of signals from individual lightning range finder. Automating the collection of information from observation points helps different services to respond more exactly to changes in thunderstorm conditions.

The uncertainty in the dipole location is fundamentally unavoidable when single-point systems are using, however these errors can be eliminated by determining its generalized parameters at two or more points [4, 6, 7]. Thus, by integrating single-point systems of passive monitoring of thunderstorm activity into a single computer network, it is possible to increase the probability of detecting a lightning discharge and the accuracy of determining its coordinates. In addition, an increase in the degree of automation of information collection from observation posts makes it possible to receive data promptly and make adequate decisions in accordance with the current thunderstorm situation.

6.1. Direction-range method

Let $\left(x,y,z\right)$ be the Cartesian coordinates of the dipole location point. Let $\left(x_i,y_i\right)$ be Cartesian coordinates of $i$ -th observation post, $i\in J.$ Let $\left(r_i,{\phi }_i,u_i\right)$ be the dipole location parameters of $i$ -th observation post. Then the system of equations

$$\left(x_i-x\right)\text{cos}{\phi }_i+\left(y_i-y\right)\text{sin}{\phi }_i=r_i{u}_i,i\in J$$ E22

holds. If the observation points do not lie on one line, then system (22) has full rank and can be solved by the least squares method

$$\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}\frac{1}{\mathrm{\Delta }}\sum _{i\in J}\left(a_i\sum _{k\in J}\text{sin}{\phi }_k\text{sin}\left({\phi }_k-{\phi }_i\right)\right)\\ \frac{1}{\Delta }\sum _{i\in J}\left(a_i\sum _{k\in J}\text{cos}{\phi }_k\text{sin}\left({\phi }_k-{\phi }_i\right)\right)\end{array}\right),$$ E23

$$\mathrm{\Delta }=\sum _{1\le k<l\le \mid J\mid }{\text{sin}}^2\left({\phi }_k-{\phi }_i\right),a_i=x_i\text{cos}{\phi }_i+y_i\text{sin}{\phi }_i,i\in J.$$

The method uses for its calculations the values $\left(r_i=\frac{\mathbf{c}}{{\alpha }_i},u_i,{\phi }_i\right)$ measured by $i$ -th autonomous lightning detection finder, and not the bearing estimation $\psi$ , as a classical direction-finding method. This makes it insensitive to the presence of an anomalous component in the observed magnetic field. There is no requirement for high-precision synchronization for the direction and distance method, which is a characteristic of the difference-range method, and the possibility of estimating the position of the lightning discharge at each point eliminates collisions in identifying the correspondence of the recorded signals to specific lightning discharges. Finally, this method makes it possible to determine the three coordinates of the location of an equivalent dipole source.

The development of this method, which makes it possible to find the angular coordinates $\left({\theta }_0,{\psi }_0\right)$ of the EMR dipole source, is given in [6, 7].

It should be noted that the use of the direction-ranging method involves knowledge of the ratio of the effective heights of magnetic and electric antennas and high requirements for the accuracy of the orientation of magnetic antennas. The determination of the ratio of effective antenna heights is a complex technical problem requiring a reference meter or a long time for the collection and processing of statistical data. Therefore, it is of interest to develop other methods for determining the location of a lightning discharge, which do not require knowledge of this relationship.

6.2. Range method

Range method for determining the coordinates of a dipole EMR source [19], in contrast to a direction- range finder, does not require a prior knowledge of the ratio of the effective antenna heights, and high accuracy of the orientation of magnetic antennas. The essence of the method is that the Cartesian coordinates $\left(x_i,y_i\right)$ of the $i$ -th point of observation, $i\in J$ . Cartesian coordinates $\left(x,y,z\right)$ of the dipole EMR source and ranges $r_i$ from the $i$ -th observation points up to the dipole EMR source satisfy the equations system

$$r_i^2={\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2+z^2,i\in J$$ E24

reducible by the introduction of the variable

$$m=\frac{x_i^2+y_i^2+z_i^2}{C}$$

to a system of linear equations

$$\mathit{mC}-2x_{i}x-2y_{i}y=r_i^2-x_i^2-y_i^2,i\in J,$$ E25

here constant $C$ is the scaling factor for reducing the condition number of the equation system.

If the observation points do not lie on one line then the resulting system will have full rank, and its least-squares method, the solution is:

$$\left(\begin{array}{c}m\\ x\\ y\end{array}\right)={\left(A^{T}A\right)}^{-1}{A}^{T}B,\text{where}A=\left(\begin{array}{c}\begin{array}{c}\begin{array}{ccc}C& -2x_1& -2y_1\end{array}\\ \begin{array}{ccc}C& -2x_2& -2y_2\end{array}\end{array}\\ \begin{array}{c}\begin{array}{ccc}⋮& ⋮& ⋮\end{array}\\ \begin{array}{ccc}C& -2x_{\left|J\right|}& -2y_{\left|J\right|}\end{array}\end{array}\end{array}\right),B=\left(\begin{array}{c}\begin{array}{c}{r}_1^2-x_1^2-y_2^2\\ r_2^2-x_2^2-y_2^2\end{array}\\ \begin{array}{c}⋮\\ r_{\left|J\right|}^2-x_{\left|J\right|}^2-y_{\left|J\right|}^2\end{array}\end{array}\right).$$ E26

At fixed coordinates of all observation points, the matrix A is constant. Therefore, even at the design stage of the system it is easy to construct an optimal algorithm for calculating unknowns, since ${\left(A^{T}A\right)}^{-1}{A}^T$ can be counted once for the implementation of the system.

For example, let observation points be

$$\left(x_1,y_1\right)=\left(-R,O\right),\left(x_2,y_2\right)=\left(0,R\right),\left(x_3,y_3\right)=\left(R,O\right),\left(x_4,y_4\right)=\left(0,-R\right).$$

Then minimum value of the condition number equal 1 and is reached for $C=2R\sqrt{2}.$ The estimates of the dipole location coordinates are the following:

$$x=\frac{r_1^2-r_3^2}{2R},y=\frac{r_2^2-r_4^2}{2R}.$$ E27

It is easy to determine the coordinates both the dipole orientation and its location with Cartesian coordinates (27) and the distances $r_i$ to the observation points $i\in J$ .

7. Conclusion

New generation of thunderstorm passive monitoring systems expand the range of problems they solve: development of new mathematical models and algorithms for analyzing thunderstorm phenomena, their tracing and display [20].