The relationship between the Durbin-Watson statistics interval and the number of corresponding signals.
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.
- location finding
- mathematical modeling
- passive monitoring of the Earth’s electromagnetic field
- electric dipole
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 .
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 . The database of ISTC project #1822  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  give the expressions for the electromagnetic field of a dipole (moment p, direction n0, location r0, δ(*) is the Dirac delta function) in vacuum
here 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(P0), H(P0)) of the source P0, and the reflected field; the latter can be represented by the field (E(P1), H(P1)) of an imaginary dipole , which is the mirror image of the dipole P1 therefore we have (see Figure 1).
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 :
here is the inverse of the time of wave propagation from the source to the observer and the variables are determined from the equalities
Other parameters are given in Figure 1.
2.2. Inverse problem
Components , , and 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 for the dipole position and 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 . Exclusion of frames with interference of signals from different sources
Parameter φ, hereinafter called pseudo-bearing. For a vertical dipole (polar angle we have For an inclined dipole, the pseudo-bearing is differs from actual direction , this difference depends on the orientation angles and . For a horizontal dipole, we have and , i.e., the pseudo-bearing is independent of the actual direction and is determined by only the dipole orientation . 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 .
To calculate the pseudo-bearing we use the following formulas:
here , [T1,T2] 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
Really let be Nx, Ny are uncorrelated white noise in the antenna channels having the same intensity N. Then
It is easy to see that for function for any has a trend due to the difference in the functions of the radiation sources. On the contrary, for , the function has no such trend. Effect of interference make autocorrelation in function . Thus, the fact of autocorrelation in the error function 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 and , 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:
We have . The values of and correspond to the cases when there is a strict positive or negative linear relationship between the shifted series and , respectively. If , 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  and containing 796,112 atmospheric spheres suitable for computation . After calculating the value of the statistics 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:
If , then the hypothesis of the existence of first order autocorrelation in a series and is assumed. If , then this hypothesis is rejected. If the value of the statistics 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|
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
Further, we assume the parameter φ to be computed, and interference is absent in the recorded signals , , and . Let us further use the function .
4.1. Primal algorithm
An analytical solution of the problem to identify parameters is given in 
Calculated parameters , φ can be used to determine the distance to the dipole , 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 . 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
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  for determining parameters u, ν, and α is calculated as follows:
Step 1. Let be eigenvalues of matrix А, are corresponding eigenvectors.
Step 2. Let be the determinant of the matrix derived from the matrix through deleting lines having numbers m and n.
Step 3. For do
Step 4. Let ,
Step 5. Let
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 . Signals are presented as integral counts of rapid values in discrete time moments, step of time discretization is , number of quantification levels is . Input signals are preliminarily processed by a smoothing filter. Usual difference and quadrature formulas with the accuracy of 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 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 .
The spectrum of dipole moment 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 for a family of bandpass filters with an amplitude-frequency characteristic
is constructed. The elements of the set are considered as an implementation of a vector random variable, and for estimating the true position parameters and 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 . This method is based on the projection of the inverse image of the vector-valued function onto the linear manifold 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 . Therefore, the inverse image of the pair of signals belongs to , 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 
Function is continuous for . Algorithm for calculating is completely stable for signals represented of 4096 readouts with the sampling time 2–10−6 s, i.e., containing 2047 harmonics with the first-harmonic discretization frequency 122 Hz.
Set contains most reliable estimates. Statistical analysis of the set of estimates makes it possible to estimate 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 .
As before, we will consider signals represented as a fast Fourier transformation containing 2047 harmonics with the first-harmonic discretization frequency 122 Hz. Let us divide the set of all harmonics into subsets
Spectral statistical algorithm to find estimates of the parameters and their variance is calculated as follows:
Step 1. For each calculate
Step 2. Calculate the noise intensity
Step 3. For each calculate
Step 4. Calculate
Step 5. For each calculate
Step 6. Calculate
Step 7. For each calculate
Step 8. Calculate
The rationale for the effectiveness of the algorithm and computational experiments are given in .
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 , 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.
This equation includes known parameters: is distance from the dipole to the observation point, are the variables of the source dipole model (3)—and unknown Cartesian coordinates 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
They are equal
The coordinates of all probable points of a source location introduce on a plane a section of a straight line (10) between points and . The set of probable location of a dipole have one-parameter representation
A measure of the coordinates set equal
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
of the dipole location in space cell
here is cell size.
It is obvious that if . To find cells in which the dipole is located with a nonzero probability, it is sufficient to find the minimal covering of the projection of onto the plane XOY by flat cells
An example of coverage is shown in the Figure 2.
Cell in accordance with equality (12) uniquely defines a space cell 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 (see Figure 1). In this case the probability density for angle is 
Thus, for one can take
Let be the set of registered signals in the time interval The integral estimate grade of radiation locus membership to cell over period let be
It is suggested to display the stream of signals as a 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 be the Cartesian coordinates of the dipole location point. Let be Cartesian coordinates of -th observation post, Let be the dipole location parameters of -th observation post. Then the system of equations
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
The method uses for its calculations the values measured by -th autonomous lightning detection finder, and not the bearing estimation , 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.
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 , 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 of the -th point of observation, . Cartesian coordinates of the dipole EMR source and ranges from the -th observation points up to the dipole EMR source satisfy the equations system
reducible by the introduction of the variable
to a system of linear equations
here constant 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:
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 can be counted once for the implementation of the system.
For example, let observation points be
Then minimum value of the condition number equal 1 and is reached for The estimates of the dipole location coordinates are the following:
It is easy to determine the coordinates both the dipole orientation and its location with Cartesian coordinates (27) and the distances to the observation points.
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 .