Kalman Filter and Wavelet Cross-correlation for VHF Broadband Interferometer Lightning Mapping

A lightning mapping system based on perpendicular crossed baseline interferometer (ITF) technology has been developed rapidly in recent years. Several processing methods have been proposed to estimate the temporal location and spatial map of lightning strikes. In this paper, a single very high frequency (VHF) interferometer is used to simulate and augment the lightning maps. We perform a comparative study of using different processing techniques and procedures to enhance the localization and mapping of lightning VHF radiation. The benchmark environment involves the use of different noise reduction and cross-correlation methods. Moreover, interpolation techniques are introduced to smoothen the correlation peaks for more accurate lightning localization. A positive narrow bipolar event (NBE) lightning discharge is analyzed and the mapping procedure is confirmed using both simulated and measured lightning signals. The results indicate that a good estimation of lightning radiation sources is achieved when using wavelet denoising and cross-correlations in wavelet-domain (CCWD) with a minimal error of 3.46 (degrees). The investigations carried out in this study confirm that the ITF mapping system could effectively map the lightning VHF radiation source.


Introduction
Lightning is a natural phenomenon that occurs in the atmosphere. When electrical discharges are generated, electromagnetic radiations over different ranges of frequencies are produced, usually extending from the Ultra-Low Frequency (ULF) through the Ultra-High Frequency (UHF) [1]. Electromagnetic fields are emitted either from cloud-to-ground (CG) strikes, intra-cloud (IC) lightning, and some in-cloud processes [2]- [6]. The physics behind lightning initiation suggested that the lightning is produced when a considerable amount of warm air is up-drafted to a height where an extremely low temperature can cause particles of frozen ice to produce electric charge disconnections [7].
Lightning strikes can cause harm to human beings and extend effects to systems operating in the surrounding environment. These effects can be minimized if reliable lightning mapping and localization systems are deployed to enable accurate prediction and detection of lightning strikes [4], [8]- [11]. Lightning mapping provides sufficient information on detected CG flashes at a particular 2 of 13 range of frequencies. The very high frequency (VHF) band has been important for lightning mapping and for studying the lightning process over the past four decades [12].
The VHF lightning mapping is carried out in two different ways: time of arrival (TOA) and direction of arrival (DOA) or what is generally known as interferometry (ITF). The VHF lightning mapping has been traditionally performed using the TOA approach and narrow-band interferometer systems to locate the radiation sources [13]- [15]. The TOA uses the time difference of arrival (TDOA) technique which requires at least four to five antennas [6]. On the other hand, the narrowband ITF system uses a closely spaced multiple baseline antennas (in order of n wavelength, nλ) to resolve the radiation measurement ambiguities and obtain the desired angular resolution [9]. The use of broadband ITF for lightning localization was first introduced by [2] and has been developed rapidly during the past few decades [16]- [20]. Unlike the narrowband ITF, the broadband ITF uses all frequency components which are then used to produce a better resolution of lightning discharges. In the latest development of ITF systems, an array of closely spaced antennas (three perpendicular antennas) are used to capture and map the lightning radiation into the two-dimensional mapping of sources [9]. Despite the significant improvements made to map lightning using interferometer systems, the mapping accuracy remains a challenge.
In this paper, we exploit the effects of different parameters of well-established methods to accurately estimate the lightning maps. We propose a modified cross-correlation method to solve problems in the DOA calculation and to improve the lightning location accuracy. We examine various noise filtering to reduce the low and high noise frequency contents. Noisy lightning signals impact the cross-correlation performance and lead to incorrect estimation of phase difference. Different generalized cross-correlation techniques are used to estimate the time difference of arrival in a segment-wise manner. These techniques include conventional time and frequency domain crosscorrelation. However, these generalized cross-correlation methods are more sensitive to the embedded noise in the signal which results in wrong time difference of arrival estimations. Therefore, we extend the proposed framework to include the wavelet-based cross-correlation which reduces the noise and improves the system location accuracy. We then compare the proposed technique to the conventional cross-correlation techniques. We also investigate the effect of interpolating the crosscorrelation function in determining the time difference. To assess the performance of the proposed lightning mapping approach, we introduce a simulation technique to generate a broadband ITF lightning. In the simulation part, we examine the effectiveness of the proposed methods in a noisy environment, where two Gaussian noises with different SNRs (in dB) are added to the simulated lightning signals.
This paper is organized as follows. In Section 2, we discuss the background of the broadband ITF, the preprocessing of the original data, the related steps of lightning mapping estimation, and the simulation procedure of lightning data. Section 3 presents the obtained lightning mapping from the measured VHF radiation sources and also a comparison with the simulation result. This paper is concluded in Section 4.

Materials and Methods
This section introduces the broadband interferometer structure and the related methods used for two-dimensional lightning mapping. Then, we present the proposed procedure for lightning data simulations. Figure 1 shows the proposed lightning mapping general procedures consisting of three stages: (1.) Pre-processing to assess the VHF signal quality, by removing noise components and performing amplitude normalization. (2.) A piecewise procedure is used to segment the lightning signals into overlapped sliding windows to perform the cross-correlations and calculate the phase difference from overlapping VHF signals time difference at each pair of antennas. Different cross-correlation techniques in time-domain (TD), frequency-domain (FD) and wavelet-domain (WD) were applied.

Pre-processing
Prior to capturing the lightning electric field (E-field) signal at the VHF interferometer system, the lightning electric field signal travels from its origin in the clouds and usually propagates through open space. The signal is affected by reflection, refraction, absorption, and other interference means. To remove the unwanted noise from the recorded lightning, signals are usually filtered using bandpass filters with different cut-off frequencies [9], [21], [22]. In this paper, we investigate three types of noise filtering approaches including bandpass filter (BPF), Kalman filter (KF) [23] and wavelet transform (WT) denoising [24].

Time Difference of Arrival
To calculate the time difference of arrival and determine the direction of the lightning VHF radiation source, the cross-correlation has been traditionally used for estimating the similarity of two crossed baselines (consisting of three antennas). This process is achieved by using sliding overlapped windows. Besides, the generalized cross-correlation technique is one of the most classical time delay estimation methods. Figure 2 illustrates the basic geometry of an orthogonal baseline for two-dimensional lightning mapping. From Figure 2, the given two baselines separated by distance d, (15 m in this study), the time difference of arrivals can be calculated from, is the phase difference of a distant signal arriving at two antennas.   can be found implicitly by performing the cross-correlation procedure. For sources located in two spatial dimensions, two angles are needed to specify the direction and height of the radiation source, namely the azimuth (Az) and elevation (El) angles. From the spherical coordinate in Figure 2, the Az and El angles can be found by, . This constrains the instances between the brackets of Eqn. (3) and (4) to have values in the range of ±1. However, any values exceed the ±1 range are considered as wrong estimations due to existing noise, and the corresponding Az and El angles can be neglected from the cross-correlation calculations.

Hardware
Figure3 ,shows the configuration diagram of the proposed ITF system which consists of four sets of antennas, band-pass filters, low noise amplifiers (LNA), coaxial cables, and oscilloscope. Two types of parallel-plate antennas of size A4 and A3 are used to capture the VHF E-field and Fast E-field (10 Hz up to 3 MHz) of lightning flashes respectively. For the VHF E-field antenna, the gap between the two plates is about h = 1 cm and 3 cm for Fast E-field to avoid the fringing effects for the VHF sensor [25,26]. The Fast E-field antenna is used to determine the type of lightning flash, and the VHF E-field antenna is used to capture the lightning radiation emissions. Figure 3, illustrates the implemented antennas B, C, and D were arranged in perpendicular baselines with two sides BC and BD to measure the phase difference. The broadband ITF system was installed in Klebang Beach, Malacca Strait, Melaka, Malaysia (2.216304° N, 102.192288° E).
The received VHF signal at each antenna is sequentially passed through preamplifier and BPF before being digitized and stored in external memory. The LNA amplifies the low-power signals preserving the signal-to-noise-ratio of the received signal. The bandwidth of the BPF is limited to the range 40-80 MHz with a center frequency of 60 MHz [17]. The filtered analog signal is then passed to the oscilloscope (PicoScope 3000 series) through three coaxial cables of 15 m (3λ) length. The output of the antenna was digitized in the oscilloscope at a sampling frequency of 250 MS/s with a vertical resolution of 8 bits. Data records were event-triggered and stored in segments of 200 ms length. Then, the retrieved data from the oscilloscope is transferred to a personal computer with MATLAB software for lightning localization analysis. A summary of the interferometer's attributes is given in Table 1. Figure 3. VHF broadband digital interferometer data acquisition system.

Simulation
In this sub-section, numerical simulations are conducted to assess the performance of the proposed mapping methods using simulated VHF interferometer signals. An actual measured signal (from antenna B) is used as a reference to simulate the other perpendicular antennas (C and D). Given handcrafted azimuth and elevation incident angles, signals of antennas C and D are a copy of antenna B with phase shift for each temporal sliding window. This simulation can provide important insight into the best methodological approach for estimating lightning mapping. Figure 4 illustrates the proposed steps to simulate the lightning mapping signals and to evaluate the performance of the estimation methods.

. Generation of Azimuth and Elevation Angles
This part presents some basic but powerful data augmentation techniques. Data augmentation is a strategy that enables practitioners to significantly increase the diversity of data available for analysis, without actually collecting new data. Various data augmentation techniques were implemented in this study such as adding Gaussian noise, scaling, and flipping and their combinations (see Figure 4). The generated Az and El angles from the hand-crafted lightning signal are assumed to be real-positive or negative values. The Gaussian noise  is an independent and identically distributed random variables, ~(0,1) with zero-mean and unit-variance, which is added to the El angles and the output ˆi For scaling augmentation, since the simulated azimuth ( Az ) and elevation ( El ) angles fall in a relatively narrow range, the outward scaling is used as it expands the azimuth and elevation range, forcing the proposed methods to make assumptions about what lies beyond the simulated values boundary. The lightning representations can be flipped horizontally and vertically. However, in practice, lightning mapping is not preferred to be flipped vertically. A horizontal flip is equivalent to rotating the lightning curve by 360 degrees. c. d.  Consequently, substituting Eqn. (7) into Eqn. (5) we can obtain the time difference between the two antennas,

Simulated ITF Signals
Given a captured lightning signal of antenna B (located at Cartesian coordinate 0 of Figure 2), Figure 5 illustrates the proposed simulation steps including, windowing, calculating the amount of temporal shift for each particular window, and using Eqn. (7) and Eqn. (8) to produce corresponding simulated step-wise signals for antennas C and D respectively. The simulated sliding windows are then concatenated to form a full track of VHF interferometer lightning signals.

Results and discussion
In this paper, we consider the use of the broadband ITF technique for the analysis of VHF lightning mapping. A comparison of various processing methods of the VHF lightning signals is performed. The simulation was conducted to investigate the impact of noise and data resolution on lightning mapping accuracy. Different filtering and cross-correlation techniques were implemented introducing new processing methods such as Kalman filter and wavelet-based cross-correlation. The best performing set of methods are validated using collected real lightning signals.
As above mentioned, the antenna B VHF signal was recorded with a sampling rate of 250 MS/s for 200 ms duration segments. Also, a Butterworth BPF filter of order 4 and cut-off frequencies of  was used along with the Kalman filter and Wavelet denoising technique. Then, the amplitudes of the signals are normalized within the range of (±1). Three cross-correlation methods are used namely, Cross-correlation in the time domain (CCTD), Cross-correlation in the frequency domain (CCFD), and Cross-correlation in the Wavelet domain (CCWD). Figure 6 shows an example of cross-correlation steps to determine the time difference of arrival using two antennas sliding windows.

Simulation results
The performance of the lightning mapping estimation was measured based on the number of sliding windows over the time course of the analyzed ITF lightning signals, where each pair of − angles are correctly estimated relative to the ground-truth simulated lightning map. The Euclidean distance between the estimated and the actual-simulated lightning map was used to measure the estimation error, where is the total number of sliding windows in the processed lightning data. To summarize the performance overall the simulated dataset, we compute the average of Euclidean distances. Table 2 summarizes the performance comparisons of the suggested lightning mapping estimation methods. We compare three main filtering approaches listed in the first column of Table  2. It can be seen that, when using the CCTD method, the KF and BPF show the lowest Euclidean distances of 3.55º and 3.63º respectively. The linear-interpolation factor used to achieve these results was 8. The WT denoising also shows relatively low distances of 3.94º when using Fejer-Korovkin function with cubic-interpolation of factor 8.
For the CCFD method, the best performance was achieved using the Kalman filter without interpolation with a distance of 3.59 º while the BPF and Wavelet denoising performed comparably to that of the CCTD method showing minimum distances of 3.76º and 3.94º, respectively.
On the other hand, using CCWD, the WT denoising shows the best performance with a distance of 3.46º using the Symlet function without any interpolations. BPF and KF also provide comparable results to that of the WT with 3.76º and 3.47º.
Generally, from Table 2, the filtering methods do not show significant enhancements in the performance when increasing the cross-correlation interpolation factor due to the relatively small window size adapted from Stock et al. [9]. The introduced wavelet filtering shows better performance than that of the BPF. However, the variation of wavelet functions and the thresholding methods show almost similar performances. On the other hand, the KF outperforms both the BPF and the WT in general. This is because KF can recursively estimate the noise regardless of the frequency contents.
In comparing the cross-correlation techniques, the CCWD method outperforms the other methods in all cases. This is due to the inherent capability of wavelet transforms to further enhance the signals and reduce the noise effects, providing a better time-domain resolution and reliable estimation of lighting maps. Moreover, comparing several wavelet functions in WT denoising helps to identify the best combination of WT parameters that fits the analyzed signals. A noticeable enhancement of cross-correlation methods was found introducing cross-correlation in the Wavelet domain for estimating the lightning maps.

Real data results
The main challenge of source localization is to identify the radiation emitted by the source out of noisy signals received by the measuring ITF system. Figure 7 shows an example of a positive narrow bipolar event (NBE). Figure 8 shows measured VHF signals related to a lightning flash triggered by the discharge process of the NBE as observed in Figure 7. This measurement represents a time expanded window of 10 μs (2500 samples). The observed signals of antenna B and D show significant similarity with possible phase shift or a time lag. However, the signal of antenna C shows slightly different noise characteristics.
In the simulation results, shown before, it is noted that the wavelet-based CCWD method produced the best performance of lightning mapping using the Symlet function with linear interpolation of level 8. Therefore, the same setup is implemented for real lightning data mapping. . shows the constructed map of the lightning source location for positive NBE. The starting position of radiation sources was in the higher altitude about 66° elevation. However, these detected sources are isolated from the well-defined lightning channel and considered outliers. The colour changes in Figure 9 (a) with time from blue to green is corresponding to the lightning progression time in Figure 8. It can be seen that the channel from the bottom to the 19° elevation was approximately straight with small changes in azimuth angles range, corresponding to the trajectory of the lightning flash. Up to the 32° elevation, the shape and the height of the channel showed certain consistency with the observation of positive NBE. In Figure 9 (b), the radiation elevation angles and fast electric field change observations (black waveform) superimposed on the VHF waveform (gray), showing the downward propagation of the VHF source. The negative electric field waveform is indicative of an upward-directed current, consistent with a positive NBE. Overall, all three filtering methods implemented in this paper show comparable accuracy for this lightning mapping with some errors  that could be resulted from the limited windowing approach conducted in this study or due to the existing in-band hardware-interference noises when collecting the data. To avoid such noise effects, it is better to properly filter the analog signal from the antennas before the signal is digitized, such that no portion of the bandwidth which contributes to the cross-correlation is contaminated by extreme noises.

CONCLUSION
Broadband interferometry for two-dimensional mapping of lightning progression was used to examine the effects of different factors on accurately estimate the lightning maps. In this paper, various lightning signal proceeding approaches are compared, including different filtration methods, cross-correlation techniques, introducing interpolation to effectively estimate the fractional time of arrival delay realizing the location of VHF lightning radiation sources in two dimensions. Simulations were conducted to provide enough analysis when comparing the performance of implemented approaches for estimating lightning maps. Aggregated results show that the CCWD outperforms the other correlation methods in estimating the time difference of arrivals regardless of the interpolations. This is because the wavelet transform not only can reduce the effective noise but also weakens the effects of the peak resolution of cross-correlation. The filtering methods also show comparable performance with Wavelet denoising which shows slightly better performance when using Symlet function. Several improvements of the current study could be achieved by involving other parameters, such as tuning the sliding window size, the windowing function, the amount of overlap.