Next Article in Journal
Assessing Storm Response of Multiple Intertidal Bars Using an Open-Source Automatic Processing Toolbox
Next Article in Special Issue
Characteristics of Regions with High-Density Initiation of Flashes in Mesoscale Convective Systems
Previous Article in Journal
Damaged Building Extraction Using Modified Mask R-CNN Model Using Post-Event Aerial Images of the 2016 Kumamoto Earthquake
Previous Article in Special Issue
Close Observation of the Evolution Process during Initial Stage of Triggered Lightning Based on Continuous Interferometer
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Technical Note

3D Lightning Location Method Based on Range Difference Space Projection

1
School of Computer Science, Chengdu Normal University, Chengdu 611130, China
2
Key Laboratory of Interior Layout Optimization and Security, Institutions of Higher Education of Sichuan Province, Chengdu Normal University, Chengdu 611130, China
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(4), 1003; https://doi.org/10.3390/rs14041003
Submission received: 11 January 2022 / Revised: 13 February 2022 / Accepted: 16 February 2022 / Published: 18 February 2022

Abstract

:
Most lightning location networks obtain the position results by optimizing the goodness of fit to determine that all combinatorial time differences of arrivals (TDOAs) are due to a common discharge. This paper proposes a three-dimensional (3D) lightning location method based on range difference (RD) space projection. The proposed method projects all the measurements into the RD space, which has the space-invariant feature of the equivalence cell and can be partitioned soundly. Aiming at the problem of computational cost of the procedure of the projection, the hierarchical strategy is proposed to improve computational efficiency. The performance of the RD space projection based on the hierarchical strategy is analyzed via Monte-Carlo simulations. The results show that the proposed method can locate lightning sources in real time with high accuracy. The results also show that the location accuracy is limited by the level of the inherent time uncertainty, the layout, and the size of the receiver network. Under the fixed layout and size of the receiver network, and the fixed measurement noise uncertainty, the positioning precision cannot be improved more even if the grid step is small enough or the number of receivers is large enough.

Graphical Abstract

1. Introduction

Lightning is a kind of strong discharge phenomenon occurring in convective weather, which is one of the most common disastrous weather phenomena in nature. With the wide application of microelectronic equipment and electrical equipment, the direct and indirect disasters caused by lightning are becoming more and more serious. Therefore, the detection and early warning of lightning have become urgent tasks to be solved.
The three-dimensional (3D) location technology of lightning provides key information for the research of the development mechanism of the lightning and early warning of the lightning because it not only can monitor the cloud-to-ground lightning (CG) well, but it can also observe the development process by using early lightning electromagnetic radiation during intra-cloud lightning (IC) [1,2,3,4]. Therefore, various studies have recently been carried on the fine location technology of the total lightning. The representative location systems are Lightning Mapping Array (LMA) system and Los Alamos Sferic Array (LASA) system. The LMA system adopted the time difference of arrival (TDOA) technique that is based on GPS clock synchronization and can detail the lightning channel structure [5,6,7,8]. The LASA system is a fast electric field location network in the VLF/LF band and is capable of high-precision, three-dimensional locations of total lightning discharge using the TOA method [9,10,11]. Later, a series of similar detection systems have been developed around the world, including the Huntsville Alabama Marx Meter Array, a Broadband Observation network, an LF near-field interferometric-TOA 3D Lightning Mapping Array, etc. [12,13,14,15,16,17,18,19]. Based on the above work, Fan et al. introduced a signal processing method (empirical mode decomposition (EMD) technology) to filter the original waveform to improve the accuracy of the peak time extraction, which further improves the location accuracy [18]. Furthermore, the time reversal (TR) technology can also be used in the 3D location of lightning discharges [20,21,22]. For example, Chen et al. used the TR method to find the optimal solution in the space limited by the linear initial solution of the TOA method [21]. Li et al. combined the TDOA technique and electromagnetic TR technique to image lightning channels [22]. In addition, the deep-learning technology has been introduced in the location of total lightning [23]. The ABA (A, B double time base) data acquisition system based on first input first output memory is adopted in a low frequency 3D lightning mapping network [24].
The lightning location based on the TOA technique can be roughly divided into two categories. The one works by solving a group of non-linear equations by some iteration [23,25,26] or analytical algorithm in principle [27,28], and the other adopts the numerical grid traverse algorithm (GTA), which builds a grid database first and then finds the grid point that matches the time differences [15,29]. Different from the iteration algorithm, the GTA does not need an initial value and is inherently accurate but has computing burdensome. Qin et al. proposed a GPU-based GTA algorithm, which has a fast location speed, although the matching ability is not improved [30].
Recently, multi-target positioning for a passive sensor network via bistatic range (BR) space projection is proposed [31,32]. This method uses the projection strategy to solve multi-target positioning problems in the viewpoint of the imaging technique. Projection strategy is the traversal method and similar to the GTA, which partitions the surveillance region into a group of grids and substitutes the grid point into equations. However, the difference between the projection strategy (PS) and GTA is that GTA obtains the positioning results by optimizing the objective function, and PS accumulates the vote for each grid and picks up the grid that has the maximum vote value as the positioning results. In addition, the BR space projection method is proposed to eliminate the space-variant feature in the geographic space in the face of a vast surveillance region.
Based on the work of [31,32], in this study, we propose a 3D lightning location method based on the range difference (RD) space projection. We present a detailed derivation procedure of the space-variant feature in the geographic space and project the TDOAs into the RD space, which is space-invariant. Considering the computational cost, a hierarchical search strategy from coarse to fine is proposed. This method has fine location ability based on the analysis of the Monte Carlo simulation.

2. Model of TDOA Lightning Detection Sensor Network

The lightning location receiver network, based on TDOA, consists of N + 1 ( N   3) spatially separated receiver stations, which is illustrated in Figure 1. The sampled time series from all the receivers are synchronized using a GPS clock.
Since a TDOA-based positioning system does not measure absolute time, but instead measures the time difference that a signal arrives at the TDOA receivers with respect to the reference TDOA receiver, the N TDOAs are expressed as
Δ t i = t i t 0 ,   i = 1 , 2 , , N ,
where t 0 is the absolute time of arrival to the reference receiver, t i is the absolute time that the lightning source arrives at the i th receiver, and N is the number of receivers, excluding the reference receiver (receiver 0). The TDOAs are converted to RDs by multiplying by c , the speed of light in air:
d i = c Δ t i = c ( t i t 0 ) ,   i = 1 , 2 , , N ,
where d i is the RDs of the lightning source to the i th receiver and the reference receiver.
Assuming that the receivers’ positions are given, we have
x r i 2 x r 0 2 = r i r 0 = d i ,
where, r i = ( x i , y i , z i ) T denotes the positions of the i th receiver, r 0 = ( 0 , 0 , 0 ) T denotes the reference receiver that is located at the origin of the coordinates, x = ( x , y , z ) T denotes the position of the lightning source, r i denotes the range of the i th receiver and the lightning source, r 0   denotes the range of the reference receiver and the lightning source, and · 2 means 2-norm.
Considering the reference receiver is located at the origin, one can create the RD equations by combining all the receivers’ information as:
{ x r 1 2 x 2 = d 1 x r 2 2 x 2 = d 2 x r N 2 x 2 = d N ,
Equation (4) is a group of hyperbolic equations that can be solved by some iteration [23,25,26] or analytical algorithm [27,28]. Unfortunately, we cannot determine that all the observed TDOAs are due to a common discharge because there might have been flipping and superposition in the different receivers. Flipping is caused by the change of relative positions of the receivers and the lightning sources. Superposition is caused by the coinciding of the different lightning flashes in some receivers. Flipping and superposition might cause some TDOAs that should be used more than once when creating the RD equation, which indicates that the relative position cannot be used as a rule to allocate TDOAs. Besides these, there might be some false TDOAs in the echoes because of the existence of noise.
Since a TDOA may participate in multiple combinations in the calculation at the same time, leading to multiple location results, the goodness of fit χ 2 is compared in all combinations to select the best result [13,28]. However, flipping, superposition, and false echoes lead to the idea that there might be dozens of candidate TDOAs in practice, and the combination of the possible RD equations is considerable (assume that there are ten receivers and five candidates for every receiver, the number of the possible RD equations is 510). Thus, we must find some new ways to solve the lightning location problem.

3. Geographic Space Projection

3.1. Geographic Space Projection

Differing from the iteration algorithm to solve the non-linear equations, the traversing algorithm discretizes the network region into a series of grids, finds the grid point that matches the TDOAs, and then votes the corresponding grid. After traversing all of the grids and equations, the grids whose votes are equal to the number of equations is the solution of the equations.
In fact, the traversing strategy is popular in the microwave imaging processing field, such as the back projection algorithm [33]. Thus, the multi-target location problem can be solved by some classical imaging algorithm, such as the back projection algorithm [34,35], which is illustrated in Figure 2.
Based on the back projection algorithm, a geographic projection is proposed, and the steps are listed as follows:
Step 1.
Partition the solution region into a group of grids by the grid step and assign a representative for each of them.
Step 2.
Select one representative, calculate its’ TDOA Δ t ˜ i according to Equation (4). Compare the calculated TDOA Δ t ˜ i with the observed TDOA Δ t i ; if the absolute difference is less than the designated threshold η 0 , | Δ t ˜ i Δ t i | < η 0 , the representative gets a vote, if else do nothing.
Step 3.
Change the receiver one by one and repeat step 2 until all receivers are processed.
Step 4.
Change the representative one by one and repeat steps 2 and 3 until the whole solution domain is processed.
Step 5.
Select the representatives whose value of vote is greater than the designated threshold η 1 . The coordinates of those representatives are the estimated lightning sources’ position.
A discretized 2D region is illustrated in Figure 3.

3.2. Equivalence Cell

Because of the discretization of the surveillance region, the RD equation should be rewritten as a group of inequalities in practice:
{ | x + δ r 1 2 x + δ 2 d 1 | ρ / 2 | x + δ r 2 2 x + δ 2 d 2 | ρ / 2 | x + δ r N 2 x + δ 2 d N | ρ / 2 ,
where δ is the variable, and ρ denotes the grid step. The collection of those δ satisfying the inequalities is named as ambiguity region in the field of imaging processing. Obviously, given a representative x , all those positions δ in the same ambiguity region are indistinguishable and should be partitioned into the same grid, which is named as equivalence cell in this paper.
By introducing a group of intermediate variables y i x + δ r i 2 x + δ 2 d i , the equivalence cell ω ( x ) can be considered as a mapping from a high-dimensional cube into the 3D space, i.e.,:
ω ( x ) = { δ | Ψ ( δ ; x ) = y ; y [ ρ 2 , ρ 2 ] N } ,
Ψ ( δ ; x ) = [ | x + δ r 1 2 x + δ 2 d 1 | | x + δ r N 2 x + δ 2 d N | ] ,
where, y [ y 1 , y 2 , , y N ] , which is a N-dimensional cube.
Given a representative x , by the multivariate Taylor expansion theorem, Ψ ( δ ; x ) can be approximated as:
ω ( x ) = { δ | Ψ ( x ) δ = y ; y [ ρ 2 , ρ 2 ] N } ,
Ψ ( x ) = [ α 1 T α 2 T α N T ] ,
α i = x r i x r i 2 x x 2 ,
Equations (8) to (10) indicate that the equivalence cell of the geographic space projection is approximately a polygon in 3D space.
To facilitate the analysis, we present the following property (seeing Appendix A for details) at first.
Property 1.
Given three receivers whose positions are r 1 , r 2 , and r 3 , and assume that there is a pair of x and δ satisfying: | α 1 T δ | ρ 2 , | α 2 T δ | ρ 2 and | α 3 T δ | ρ 2 , then | α T δ | ρ 2 (approximately), where α denotes a point in the triangle with vertices r 1 , r 2 , and r 3 .
Property 1 indicates that the shape of the equivalence cell is dominant by the three receivers covering the largest area. In other words, one cannot improve the positioning precision via increasing the number of receivers in the triangle (compatible with the discussion on the position dilution of precision (PDOP) [36] of GPS system).
Under this property, we can continue our discussion on the equivalence cell by examining only four receivers.
The geometry of α i is shown in Figure 4. Given a representative x , receiver r i , and reference point 0, the α i associated with Equation (10) is the lower side of the isosceles triangle that is determined by the unit vectors of x r i and x . When the three sensors are not collinear, α 1 , α 2 , and α 3 are not co-planar. Thus, the equivalence cell is the linear transform of a cube, which is a parallel hexahedron and is shown in Figure 5. The eight vertices of the hexahedron can be calculated as:
δ V = Ψ 1 ( x ) y V ,
where y V denotes the vertices of the cube.
Equation (11) indicates that the size of the equivalence cell is proportional to the grid step; the smaller the grid step is, the smaller the equivalence cell is. On the other hand, the equivalence cell is space-variant with respect to the representative   x (seeing Section 5 for details). This space-variant characteristic makes it difficult to partition the solution domain in the geographic space effectively. A fine partition means that there might be several representatives for some equivalence cells (for those regions far away from the receiver network); a coarse partition means that there might be no representative for some equivalence cells (for those regions near the receiver network). Thus, a novel projection algorithm is proposed in the next section.

4. RD Space Projection

According to the discussion in Section 3.2, it is difficult to partition the solution domain soundly, as the equivalence cell in the geographic space is irregular and space-variant. On the contrary, the equivalence cells of the intermediate variable y i (corresponding to the RD of different receivers) are all regular and space-invariant (N-dimensional cube). Therefore, we can project the data into the RD space, rather than the geographic space, which can be defined as:
Ω R n g = { y | y = Ψ ( x ) , x Ω G e o } ,
{ y 1 = x + δ r 1 2 x + δ 2 y 2 = x + δ r 2 2 x + δ 2 y 3 = x + δ r 3 2 x + δ 2 ,
where Ω R n g denotes the solution domain in the RD space, Ω G e o denotes the solution domain in the geographic space, and Equation (13) gives the mapping between the RD space and the geographic space. Obviously, given an x or y , there is only one y or one x (only one is reasonable) corresponding to it. Thus, we can traverse Ω R n g one by one to substitute traversing Ω G e o .
In the RD space, we can simply find a partition Ξ as:
Ξ { ω i , j , k , i , j , k } ,
ω i , j , k { y y = y i , j , k + δ y , δ y [ ρ 2 , ρ 2 ] } ,
y i , j , k [ ρ i , ρ j , ρ k   ] T ,
where y i , j , k is the representative of ω i , j , k .
For the N + 1 ( N   3) receivers network, the three whose triangle can cover as many other receivers as possible are selected, and their positions are used to create the mapping by Equation (13). According to Property 1, if an equivalence cell’s counterpart satisfies the equations in Equation (5) associated with the three sensors, it satisfies the other equations approximately too.
In addition, the geographic projection is computing costly and its overall performance, such as the speed and location precision, relies heavily on the grid step and the surveillance region. Therefore, we proposed the hierarchical search strategy to decrease the grid step and shrink the solution region gradually. At the first, we partition the solution domain by using the large grid step and obtain the coarse location. Then, we shrink the solution domain into some specific regions. We partition the new solution domain by using another smaller grid step and obtain the new location until the accuracy of the position result is satisfied the system requirement. The computational cost can be reduced geometrically by the hierarchical strategy iteratively.
According to the above discussion, the steps of the RD space projection based on hierarchical strategy are listed as follows:
Step 1.
Select three receivers whose triangle covers as many receivers as possible.
Step 2.
Calculate the representatives in the RD space according to Equation (16) and solve the non-linear Equation (13) to obtain representatives’ counterparts in the geographic space.
Step 3.
Calculate the RDs using the representatives’ counterparts via Equation (4) and perform the geographic space projection discussed in Section 3.1.
Step 4.
If the votes of the representatives are greater than the threshold, store the positions of them. Set a new solution domain with some specific size of a 3D cube and center at those positions. Change the grid step to a smaller one. Repeat step 1 to step 4 of the geographic space projection by using the new solution domain and the new grid step until the accuracy of the position result satisfies the system requirement. Then, the representatives’ positions whose votes are greater than the threshold are the estimated lightning sources’ positions.

5. Performance Analysis

5.1. Space-Invariant Feature of the RD Space

This subsection demonstrates the space-invariant feature of the RD space. Assume that there are 10 receivers. Figure 6, Figure 7, and Figure 8a plot geometries of the receiver network with the lightning flash, one of the receivers (the red circle) is placed at the origin, three of the receivers (the cyan circles) are located at [−50, 0, 0] km, [50, 0, 0] km, and [0, 87, 0] km. the others (the black circles) are distributed uniformly in the triangle determined by the three receivers, and the lightning flash (the blue star) is located at [0, 0, 6] km, [100, 120, 6] km, and [300, 200, 6] km, respectively. The equivalence cells (ambiguity region whose values are 0.7 times larger than the maximum value) in the RD space and the geographic space are shown in Figure 6, Figure 7, and Figure 8b,c respectively.
From Figure 6, Figure 7, and Figure 8, we find that the equivalence cells in both the RD space and the geographic space are points when the lightning source is located at [0, 0, 6] km. The equivalence cell in the geographic space diffuses significantly, but the one in the RD space is still a point when the source is located at [100, 120, 6] km. The equivalence cell in the geographic space diffuses in a mass region, but the one in the RD space is still a point when the source is located at [300, 200, 6] km.
Assume that the surveillance region ranges from [0, 0, 0] km to [300, 200, 6] km; if a coarse grid step is selected, some lightning sources located in those regions near to the receiver network might be lost. On the contrary, if a fine grid step is selected, a lot of projection operations are performed for an equivalence cell in those regions far from the network in the geographic space, which is useless for positioning. In the RD space, the equivalence cell can easily avoid the redundant operations and lightning source loss as it is space-invariant.
It is worth mentioning that, since one needs to solve the geographic positions from the coordinates in the RD space via Equation (13), the space-invariant feature of the RD space does not mean that its positioning precision is superior to that in the geographic space. The theoretical analysis on the positioning precision of the receiver network is similar to that of the GPS system, i.e., geometric dilution of precision (GDOP), position dilution of precision (PDOP), etc., which can be found in the literature [36]. Further discussion of the positioning precision is presented in Section 5.2 based on some numerical experiments.

5.2. Positioning Performance Analysis

This subsection discusses the influences of the additive timing noise, number of receivers, and lightning source positions on the positioning performance. The surveillance regions are from [−200, −200, 0] km to [200, 200, 10] km. Assume that there are 5, 10, 20, and 30 receivers, respectively; one of the receivers is placed at the origin, three of the receivers are located at [−50, −50, 0] km, [50, −50, 0] km, and [0, 50, 0] km, and the others are distributed uniformly in the triangle determined by the three receivers. The baseline of the network is 100 km. Use four-level hierarchical grids. The first-level grid is 41 × 41 × 10 with a grid step of 10 × 10 × 1km to cover the whole surveillance region. The second-level grid is 100 × 100 × 10 with a grid step of 1000 × 1000 × 1000 m, the third-level grid is 100 × 100 × 100 with a grid step of 100 × 100 × 100 m, and the fourth-level grid is 100 × 100 × 100 with a grid step of 10 × 10 × 10 m.
Firstly, we discuss that there must be just one lightning source in the surveillance region. For analysis, the different lightning source positions on the positioning performance partition the surveillance region into 41 × 41 × 10 grids by interval 10 × 10 × 1 km. Assume that the lightning flash occurred in each grid one by one, and each time the source’s height is fixed at 6 km in the z dimension.
The simulation steps are listed as follows:
Step 1:
Select one grid, record the geographic position, and substitute the position in the Equation (4) to calculate the TDOA.
Step 2:
Add some gauss noise with standard deviation 100 ns in TDOA calculated in Step 1 to simulate the real observed measurement Δ t i (according to the GPS measurement error).
Step 3:
Use the RD space projection algorithm based on the hierarchical strategy discussed in Section 3.1 and Section 4.
Step 4:
Extract the cells whose voting values are greater than the threshold. Store those cells’ position as the estimated lightning source’s positions and the corresponding onset time.
Step 5:
Repeat step 2 to step 4 1000 times, calculate the root mean square error (RMSE) in x-y and z between the simulated lightning source position and the estimated source position. The x-y RMSE and the z RMSE were calculated according to
RMSE k x y = 1 J j = 1 J ( ( x k j x ^ k j ) 2 + ( y k j y ^ k j ) 2 ) ,
RMSE k z = 1 J j = 1 J ( z k j z ^ k j ) 2 ,
RMSE k = ( RMSE k x y ) 2 + ( RMSE k z ) 2 ,
where x k j , y k j , and z k j are the x , y , and z   positions of the simulated lightning source; x ^ k j , y ^ k j , and z ^ k j are the x, y, and z positions of the estimated lightning source by the RD space projection at j th and k th simulation, respectively; j = { 1 , 2 , , J } is the number of times of Monte-Carlo simulation of step 2 to step 4; J = 1000 , k = { 1 , 2 , , K } is the cell’s number of the partitioned surveillance region; K = 41 × 41 × 10 = 16 , 810 is the whole grids’ number.
Step 6:
Repeat step 1 to step 5 K times, until all surveillance region grids are processed, store all RMSE k x y , RMSE k z , and RMSE k .
The maximum RMSE max x y , RMSE max z , and RMSE max with the number of receivers 5, 10, 20, and 30 outside and inside the receiver network are listed in Table 1 and Table 2, respectively.
From Table 1 and Table 2, we find that, as the grid step decreases, the RMSE max is smaller and the location accuracy is higher. When the grid step is 10 × 10 × 10 m, the position results inside the network with 5, 10, 20, and 30 receivers are less than 38 m location accuracy as the Gaussian timing error is 100 ns, and the position results outside the network with 5, 10, 20, and 30 receivers are less than 38 m for the horizontal and less than 500 m for the vertical. To the authors’ knowledge, the results are sound because the location accuracy is ~100 m in the coverage region of the network and ~1000 m out of the coverage region of the network by the existing algorithms. In addition, due to the measurement error of the arrival time, we cannot improve the location accuracy by increasing the receivers or decreasing the grid step further. Besides, we find that when the grid steps are 10,000 × 10,000 × 1000 m, 1000 × 1000 × 1000 m, and 100 × 100 × 100 m iteratively, the location accuracy becomes significantly worse as the receiver decreased. The location errors outside the network are much larger than the location errors inside the network under the same receivers and the same grid step, and the vertical position accuracy is less than the horizontal position accuracy under 5 to 30 receivers and different grid steps no matter inside or outside the network, except for the results inside the network under the 10×10×10m grid step. Especially, the RMSE max x y and RMSE max z outside the network are larger than the current grid interval. It means that when we center at the last estimated points, shrink the resolve domain into 100 × 100 × 100 with the last grid step, and start the next searching by using smaller interval, the true location of the flash may be out of the new resolve domain. Then, we must risk losing the lightning source at the next hierarchical searching.
For discussing the location errors further, the location errors under 100 × 100 × 100 m grid searching step and 10 × 10 × 10 m grid searching step with 5, 10, and 30 receivers are plotted in Figure 9, and Figure 10, respectively. Figure 9, and Figure 10a,c,e are the horizontal location errors, and Figure 9, and Figure 10b,d,f are the vertical location errors. From Figure 9, and Figure 10, we find that, as the distance from the lighting source to the central station increases, the location accuracy gradually decreases. The RMSE max x y and RMSE max z of the whole detection domain appear at the left and right corners. Figure 9, and Figure 10a,c,e show that the horizontal location error within a radius of 150 km is less than 50 m, and Figure 9, and Figure 10b,d,f show that the vertical location error within a radius of 150 km is less than 1000 m. This region is three times the baseline of the receiver network. The result can guide us to further subdivide the searching domain. For example, the searching domain in the radius of three times of the baseline can be set as 10 × 10 × 10 with the current grid step, and the searching domain out of the radius of three times of the baseline can be set as 100 × 100 × 100 with the current grid step, in order to not lose the source and balance the computational cost.
Computational cost is a big issue of the grid traverse algorithm. To balance the accuracy and computation, combining the results of Table 1 and Table 2, and Figure 9, and Figure 10, a 20 × 20 × 20 searching grid is accepted. Table 3 gives the processing time of one receiver for a one-time Monte-Carlo simulation. Because each receiver has the same computation procedure, it means N receivers can be processed parallelly. Therefore, under the parallel processing, the total processing time of the proposed algorithm is 0.8203 s. Moreover, for the first searching level, because the searching domain is fixed and determined by the user, we can calculate the results of the grid traverse and store them in advance. Then, the total processing time is reduced to 0.3536 s. It satisfies the requirement of the real-time lightning monitoring. Furthermore, we can divide the searching domain into two parts according to the results of Figure 9, and Figure 10 to reduce the computational cost further.
Finally, we discuss that there are two lightning sources in the surveillance region simultaneously. The 1000-time Monte-Carlo simulations are conducted with the standard deviation of timing noise 100 ns and the number of receivers 10. Two lightning flashes are distributed uniformly in a 3D cube with size 400 × 400 × 400 m3 and centered at [100, 100, 6] km. Figure 11 plots the positioning result of the RD space algorithm. From them, we find that the algorithm can position the two lightning sources correctly.

6. Conclusions

In this paper, a 3D lightning location method based on RD space projection, which combines with the hierarchical strategy, is proposed. The main conclusions of this paper are presented as follows:
  • Based on the GPS principle, the receiver network can position the lightning flash in 3D space. Different from most lightning location networks, it obtained the position results by optimizing the goodness of fit, and the projection strategy projects the TDOAs into the RD space to accumulate the echoes’ energy and overcomes the problem of how to determine all combinatorial TDOAs due to a common discharge.
  • It is difficult to partition the geographic space soundly because of the space-variant feature of the geographic space in face of a vast surveillance region. On the contrary, it is an easy task in the RD space due to its’ space-invariant feature. Thus, the proposed algorithm adopts the grid traverse algorithm to traverse the geographic space by traversing the RD space.
  • Increasing the number of receivers can promote the anti-noise performance but cannot improve the positioning precision. The location accuracy is limited by the level of the inherent time uncertainty, the layout, and the size of the receiver network. Due to the measurement time uncertainty, we cannot improve the location accuracy by increasing the receivers or decreasing the grid step further under the fixed size of the receiver network.

Author Contributions

Conceptualization, C.Z. and L.F.; methodology, L.F. and C.Z.; software, L.F. and C.Z.; validation, L.F. and C.Z.; formal analysis, L.F. and C.Z.; investigation, L.F. and C.Z.; data curation, L.F. and C.Z.; writing—original draft preparation, L.F.; writing—review and editing, C.Z.; funding acquisition, L.F. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the High-Level Talent Project of Chengdu Normal University (YJRC2020-12), the National Natural Science Foundation of China (51978089), the Innovation Team of Chengdu Normal University (CSCXTD2020B09), Department of Science and Technology key research and development project of Sichuan Provincial (22ZDYF2726).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing is not applicable to this article.

Acknowledgments

The authors are thankful to all members of the Key Laboratory of Interior Layout Optimization and Security of Institutions of Higher Education of Sichuan Province for their challenging conversation and comments.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

This appendix will discuss the Property 1.
Defining κ 1 x r 2 / x r 1 2 , we have:
| ( x r 1 x r 1 2 x x 2 ) T δ | = | ( κ 1 ( x r 1 ) x r 2 x x 2 ) T δ | ρ 2 ,
Defining δ x x T δ / x 2 , we have:
ρ 2 + δ x κ 1 ( x r 1 ) T x r 2 δ ρ 2 + δ x ,
i.e.,
ρ 2 + δ x κ 1 ( x r 1 ) T x r 2 δ ρ 2 + δ x κ 1 ,
ρ 2 + δ x κ 1 δ x ( x r 1 x r 2 x x 2 ) T δ ρ 2 + δ x κ 1 δ x ,
| ( x r 1 x r 2 x x 2 ) T δ | ρ 2 κ 1 + | ( 1 κ 1 ) δ x κ 1 | ,
Similarly, we have:
| ( x r 2 x r 2 x x 2 ) T δ | ρ 2 κ 2 + | ( 1 κ 2 ) δ x κ 2 | ,
| ( x r 3 x r 2 x x 2 ) T δ | ρ 2 κ 3 + | ( 1 κ 3 ) δ x κ 3 | ,  
Since r is in the triangle with vertices r 1 , r 2 , and r 3 , i.e.,
r = β ( α r 1 + ( 1 α ) r 2 ) ( 1 β ) r 3 ; α [ 0 ,   1 ] , β [ 0 ,   1 ] , we have:
α β | ( x r 1 x r 2 x x 2 ) T δ | α β ρ max 2 ,
β ( 1 α ) | ( x r 1 x r 2 x x 2 ) T δ | β ( 1 α ) ρ max 2 ,
( 1 β ) | ( x r 1 x r 2 x x 2 ) T δ | ( 1 β ) ρ max 2 ,
Thus,
| ( x r x r 2 x x 2 ) T δ | ρ max 2 ,
where ρ max 2 max [ { ρ 2 κ q + | ( 1 κ q ) δ x κ q | } ] , q = 1 , 2 , 3 .
Equation (A11) provides an accurate bound but is too complex to be used. When the range from the lightning strike to the receiver is far larger than the size of the receiver network, we approximate that κ q 1 ,   q = 1 , 2 , 3 , then ρ max ρ , and we have Property 1.

References

  1. Wu, T.; Wang, D.; Takagi, N. Locating preliminary breakdown pulses in positive cloud-to-ground lightning. J. Geophys. Res. Atmos. 2018, 123, 7989–7998. [Google Scholar] [CrossRef]
  2. Karunarathne, S.; Marshall, T.C.; Stolzenburg, M.; Karunarathna, N.; Vickers, L.E.; Warner, T.A.; Orville, R.E. Locating initialbreakdown pulses using electric field change network. J. Geophys. Res. Atmos. 2013, 118, 7129–7141. [Google Scholar] [CrossRef]
  3. Zheng, D.; Shi, D.; Zhang, Y.; Zhang, Y.; Lyu, W.; Meng, Q. Initial leader properties during the preliminary breakdown processes of lightning flashes and their associations with initiation positions. J. Geophys. Res. Atmos. 2019, 124, 8025–8042. [Google Scholar] [CrossRef]
  4. Zheng, D.; MacGorman, D.R. Characteristics of flash initiations in a super cell cluster with tornadoes. Atmos. Res. 2016, 167, 249–264. [Google Scholar] [CrossRef] [Green Version]
  5. Wu, P. The Lightning Direction-Finding Location System and Time Difference Location System. High Volt. Eng. 1995, 21, 3–7. [Google Scholar]
  6. Chen, J.H.; Zhang, Q.; Feng, W.X.; Fang, Y.H. Lightning location system and lightning detection network of China power grid. High Volt. Eng. 2008, 34, 425–431. [Google Scholar]
  7. Rison, W.; Thomas, R.J.; Krehbiel, P.R.; Hamlin, T.; Harlin, J. A GPS-based three-dimensional lightning mapping system: Initial observations in central New Mexico. Geophys. Res Lett. 1999, 26, 3573–3576. [Google Scholar] [CrossRef] [Green Version]
  8. Thomas, R.J.; Krehbiel, P.R.; Rison, W.; Hunyady, S.J.; Winn, W.P.; Hamlin, T.; Harlin, J. Accuracy of the Lightning Mapping Array. J. Geophys. Res. 2004, 109, D14207. [Google Scholar] [CrossRef] [Green Version]
  9. Shao, X.M.; Stanley, M.; Regan, A.; Harlin, J.; Pongratz, M.; Stock, M. Total lightning observations with the new and improvedLos Alamos Sferic Array (LASA). J. Atmos. Ocean. Technol. 2006, 23, 1273–1288. [Google Scholar] [CrossRef]
  10. Smith, D.A.; Eack, K.B.; Harlin, J.; Heavner, M.J.; Jacobson, A.R.; Massey, R.S.; Shao, X.M.; Wiens, K.C. The Los Alamos Sferic Array: A research tool for lightning investigations. J. Geophys. Res. 2002, 107, ACL 5-1–ACL 5-14. [Google Scholar] [CrossRef] [Green Version]
  11. Smith, D.A.; Heavner, M.J.; Jacobson, A.R.; Shao, X.M.; Massey, R.S.; Sheldon, R.J.; Wiens, K.C. A method for determining intracloud lightning and ionospheric heights from VLF/LF electric field records. Radio Sci. 2004, 39, RS1010. [Google Scholar] [CrossRef]
  12. Bitzer, P.M.; Christian, H.J.; Stewart, M.; Burchfield, J.; Podgorny, S.; Corredor, D.; Hall, J.; Kuznetsov, E.; Franklin, V. Characterization and applications of VLF/LF source locations from lightning using the Huntsville Alabama Marx Meter Array. J. Geophys. Res. Atmos. 2013, 118, 3120–3138. [Google Scholar] [CrossRef]
  13. Yoshida, S.; Wu, T.; Ushio, T.; Kusunoki, K.; Nakamura, Y. Initial results of LF sensor network for lightning observation and characteristics of lightning emission in LF band. J. Geophys. Res. Atmos. 2014, 119, 12034–12051. [Google Scholar] [CrossRef]
  14. Lyu, F.C.; Cummer, S.A.; Lu, G.P.; Zhou, X.; Weinert, J. Imaging lightning intracloud initial stepped leaders by low-frequency interferometric lightning mapping array. Geophys. Res. Lett. 2016, 43, 5516–5523. [Google Scholar] [CrossRef]
  15. Lyu, F.C.; Cummer, S.A.; Solanki, R.; Weinert, J.; McTague, L.; Katko, A.; Barrett, J.; Zigoneanu, L.; Xie, Y.; Wang, W. A low-frequency near-field interferometric-TOA 3-D lightning mapping array. Geophys. Res. Lett. 2014, 41, 7777–7784. [Google Scholar] [CrossRef]
  16. Wu, T.; Wang, D.H.; Takagi, N. Lightning mapping with an array of fast antennas. Geophys. Res. Lett. 2018, 45, 3698–3705. [Google Scholar] [CrossRef]
  17. Fan, X.P.; Zhang, Y.J.; Zheng, D.; Zhang, Y.; Lyu, W.T.; Liu, H.Y.; Xu, L.T. A new method of three-dimensional location for low frequency electric field detection array. J. Geophys. Res. Atmos. 2018, 123, 8792–8812. [Google Scholar] [CrossRef]
  18. Shi, D.D.; Zheng, D.; Zhang, Y.; Zhang, Y.J.; Huang, Z.G.; Lu, W.T.; Chen, S.D.; Yan, X. Low-frequency E-field Detection Array (LFEDA)— Construction and preliminary results. Sci. China Earth Sci. 2017, 60, 1896–1908. [Google Scholar] [CrossRef]
  19. Zhu, Y.; Bitzer, P.; Stewart, M.; Podgorny, S.; Corredor, D.; Burchfield, J.; Carey, L.; Medina, B.; Stock, M. Huntsville Alabama Marx Meter Array 2: Upgrade and Capability. Earth Space Sci. 2020, 7, e2020EA001111. [Google Scholar] [CrossRef] [Green Version]
  20. Liu, B.; Shi, L.H.; Qiu, S.; Liu, H.Y.; Sun, Z.; Guo, Y.F. Three-dimensional lightning positioning in low-frequency band using time reversal in frequency domain. IEEE Trans. Electromagn. Compat. 2020, 62, 774–784. [Google Scholar] [CrossRef]
  21. Chen, Z.; Zhang, Y.; Zheng, D.; Zhang, Y.; Fan, X.; Fan, Y.; Xu, L.; Lyu, W. A method of three-dimensional location for LFEDA combining the time of arrival method and the time reversal technique. J. Geophys. Res. Atmos. 2019, 124, 6484–6500. [Google Scholar] [CrossRef]
  22. Li, F.; Sun, Z.; Liu, M.; Yuan, S.; Wei, L.; Sun, C.; Lyu, H.; Zhu, K.; Tang, G. A new hybrid algorithm to image lightning channels combining the time difference of arrival technique and electromagnetic time reversal technique. Remote Sens. 2021, 13, 4658. [Google Scholar] [CrossRef]
  23. Wang, J.; Zhang, Y.; Tan, Y.; Chen, Z.; Zheng, D.; Zhang, Y.; Fan, Y. Fast and fine location of total lightning from low frequency signals based on deep-learning enconding features. Remote Sens. 2021, 13, 2212. [Google Scholar] [CrossRef]
  24. Ma, Z.; Jiang, R.; Qie, X.; Xing, H.; Liu, M.; Sun, Z.; Qin, Z.; Zhang, H.; Li, X. A low frequency 3D lightning mapping network in north China. Atmos. Res. 2021, 249, 105314. [Google Scholar] [CrossRef]
  25. Koshak, W.J.; Blakeslee, R.J.; Bailey, J.C. Data retrieval algorithms for validating the optical transient detector and the lightning imaging sensor. J. Atmos. Ocean. Technol. 2000, 17, 279–297. [Google Scholar] [CrossRef] [Green Version]
  26. Koshak, W.J.; Blakeslee, R.J. TOA lightning location retrieval on spherical and oblate spheroidal earth geometries. J. Atmos. Ocean. Technol. 2001, 18, 187–199. [Google Scholar] [CrossRef]
  27. Chan, Y.T.; Ho, K.C. A simple and efficient estimator for hyperbolic location. IEEE Trans. Signal Process. 1994, 42, 1905–1915. [Google Scholar] [CrossRef] [Green Version]
  28. Wang, Y.; Qie, X.; Wang, D.; Liu, M.; Su, D.; Wang, Z.; Tian, Y. Beijing Lightning Network (BLNET) and the observation on preliminary breakdown processes. Atmos. Res. 2016, 171, 121–132. [Google Scholar] [CrossRef]
  29. Qin, Z. The Observation of Ionospheric D-Layer Based on the Multi-Station Lightning Detection System. Master’s Degree Dissertation, Earth and Space Science, University of Science and Technology of China, Hefei, China, 2014. [Google Scholar]
  30. Qin, Z.; Chen, M.; Lyu, F.; Cummer, S.A.; Zhu, B.; Du, Y. A GPU-Based Grid Traverse Algorithm for Accelerating Lightning Geolocation Process. IEEE Trans. Electromagn. Compat. 2020, 62, 489–497. [Google Scholar] [CrossRef]
  31. Shi, J.; Fan, L.; Zhang, X.L.; Shi, T.Y. Multi-target positioning for passive sensor network via bistatic range space projection. Sci. China Inf. Sci. Lett. 2016, 59, 1–3. [Google Scholar] [CrossRef]
  32. Shi, J.; Zhao, W.; Zhang, X.L. Space-invariant projection method for multi-target positioning. In Proceedings of the TENCON IEEE Region 10 Conference, Xi’an, China, 23 January 2014. [Google Scholar]
  33. Hong, I.K.; Chung, S.T.; Kim, H.K.; Kim, Y.B.; Son, Y.D.; Cho, Z.H. Fast forward projection and backward projection algorithm using SIMD. In Proceedings of the 2006 IEEE Nuclear Science Symposium Conference Record, San Diego, CA, USA, 29 October–1 November 2006; Volume 6, pp. 3361–3368. [Google Scholar]
  34. Hong, I.K.; Chung, S.T.; Kim, H.K.; Kim, Y.B.; Son, Y.D.; Cho, Z.H. Ultra-Fast Symmetry and SIMD-Based Projection-Back projection (SSP) Algorithm for 3-D PET Image Reconstruction. IEEE Trans. Med. Imaging 2007, 26, 789–803. [Google Scholar] [CrossRef] [PubMed]
  35. Bai, X.; Xing, M.; Zhou, F.; Bao, Z. High-resolution three-dimensional imaging of spinning space debris. IEEE Trans. Geosci. Remote Sens. 2009, 47, 2352–2362. [Google Scholar]
  36. Elliott, D. Kaplan. In Understanding GPS: Principles and Applications, 2nd ed.; Artech House: Boston, MA, USA, November 2005. [Google Scholar]
Figure 1. Illustration of the receiver network.
Figure 1. Illustration of the receiver network.
Remotesensing 14 01003 g001
Figure 2. Illustration of the back projection algorithm. A and B are echoes of the lightning flash, F is the false alarm.
Figure 2. Illustration of the back projection algorithm. A and B are echoes of the lightning flash, F is the false alarm.
Remotesensing 14 01003 g002
Figure 3. Illustration of the discretized 2D region in a lightning location network.
Figure 3. Illustration of the discretized 2D region in a lightning location network.
Remotesensing 14 01003 g003
Figure 4. Geometry of α i .
Figure 4. Geometry of α i .
Remotesensing 14 01003 g004
Figure 5. Shape of the equivalence cell with N = 3.
Figure 5. Shape of the equivalence cell with N = 3.
Remotesensing 14 01003 g005
Figure 6. (a) Geometry of the receiver network with the lightning flash at [0, 0, 6] km; (b) in the RD space, the equivalence cell is a point; (c) in the geographic space, the equivalence cell is almost a point.
Figure 6. (a) Geometry of the receiver network with the lightning flash at [0, 0, 6] km; (b) in the RD space, the equivalence cell is a point; (c) in the geographic space, the equivalence cell is almost a point.
Remotesensing 14 01003 g006
Figure 7. (a) Geometry of the receiver network with the lightning flash at [100, 120, 6] km; (b) in the RD space, the equivalence cell is a point; (c) in the geographic space, the equivalence cell diffuses significantly.
Figure 7. (a) Geometry of the receiver network with the lightning flash at [100, 120, 6] km; (b) in the RD space, the equivalence cell is a point; (c) in the geographic space, the equivalence cell diffuses significantly.
Remotesensing 14 01003 g007
Figure 8. (a) Geometry of the receiver network with the lightning flash at [300, 200, 6] km; (b) in the RD space, the equivalence cell is a point; (c) in the geographic space, the equivalence cell diffuses in a mass region.
Figure 8. (a) Geometry of the receiver network with the lightning flash at [300, 200, 6] km; (b) in the RD space, the equivalence cell is a point; (c) in the geographic space, the equivalence cell diffuses in a mass region.
Remotesensing 14 01003 g008
Figure 9. (a) RMSE k x y under five receivers; (b) RMSE k z under five receivers; (c) RMSE k x y under 10 receivers; (d) RMSE k z under 10 receivers; (e) RMSE k x y under 30 receivers; (f) RMSE k z under 30 receivers. The location errors under the 100 × 100 × 100 m grid step. The triangles indicate the receivers’ locations. The white triangle shows locations at the origin, and the three red triangles form a triangle that covers the others, whose color is green.
Figure 9. (a) RMSE k x y under five receivers; (b) RMSE k z under five receivers; (c) RMSE k x y under 10 receivers; (d) RMSE k z under 10 receivers; (e) RMSE k x y under 30 receivers; (f) RMSE k z under 30 receivers. The location errors under the 100 × 100 × 100 m grid step. The triangles indicate the receivers’ locations. The white triangle shows locations at the origin, and the three red triangles form a triangle that covers the others, whose color is green.
Remotesensing 14 01003 g009
Figure 10. (a) RMSE k x y under five receivers; (b) RMSE k z under five receivers; (c) RMSE k x y under 10 receivers; (d) RMSE k z under 10 receivers; (e) RMSE k x y under 30 receivers; (f) RMSE k z under 30 receivers. The location errors under the 10 × 10 × 10 m grid step. The white triangle shows locations at the origin, and the three red triangles form a triangle that covers the others, whose color is green.
Figure 10. (a) RMSE k x y under five receivers; (b) RMSE k z under five receivers; (c) RMSE k x y under 10 receivers; (d) RMSE k z under 10 receivers; (e) RMSE k x y under 30 receivers; (f) RMSE k z under 30 receivers. The location errors under the 10 × 10 × 10 m grid step. The white triangle shows locations at the origin, and the three red triangles form a triangle that covers the others, whose color is green.
Remotesensing 14 01003 g010
Figure 11. The positioning results of two lightning discharges.
Figure 11. The positioning results of two lightning discharges.
Remotesensing 14 01003 g011
Table 1. Maximum RMSE under different grid step and number of receivers inside the receiver network.
Table 1. Maximum RMSE under different grid step and number of receivers inside the receiver network.
Grid Step(m)10,000 × 10,000 × 10001000 × 1000 × 1000100 × 100 × 10010 × 10 × 10
RMSE(m)
RMSE max x y 543826657038
1038934106336
2030672985432
3027821564331
RMSE max z 5867483610939
1076325409238
2067563417533
3060142636531
Table 2. Maximum RMSE under different grid step and number of receivers outside the receiver network.
Table 2. Maximum RMSE under different grid step and number of receivers outside the receiver network.
Grid Step(m)10,000 × 10,000 × 10001000 × 1000 × 1000100 × 100 × 10010 × 10 × 10
RMSE(m)
RMSE max x y 57774243120852
107725186015255
207177133813558
3069589278356
RMSE max z 516,45413,4324336443
1016,23089722998487
2015,88646731931474
3015,58227821817471
Table 3. The processing time of one receiver.
Table 3. The processing time of one receiver.
Searching LevelSearching GridsGrid Step (km)Searching Domain (km)Processing Time (s)
First level41 × 41 × 1010 × 10 × 1410 × 410 × 100.4667
Second level20 × 20 × 201 × 1 × 120 × 20 × 200.1098
Third level20 × 20 × 200.1 × 0.1 × 0.12 × 2 × 20.1153
Fourth level20 × 20 × 200.01 × 0.01 × 0.010.2 × 0.2 × 0.20.1285
Total processing time 0.8203
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Fan, L.; Zhou, C. 3D Lightning Location Method Based on Range Difference Space Projection. Remote Sens. 2022, 14, 1003. https://0-doi-org.brum.beds.ac.uk/10.3390/rs14041003

AMA Style

Fan L, Zhou C. 3D Lightning Location Method Based on Range Difference Space Projection. Remote Sensing. 2022; 14(4):1003. https://0-doi-org.brum.beds.ac.uk/10.3390/rs14041003

Chicago/Turabian Style

Fan, Ling, and Changhai Zhou. 2022. "3D Lightning Location Method Based on Range Difference Space Projection" Remote Sensing 14, no. 4: 1003. https://0-doi-org.brum.beds.ac.uk/10.3390/rs14041003

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop