Three-Dimensional Mapping on Lightning Discharge Processes Using Two VHF Broadband Interferometers

: Lightning Very-high-frequency (VHF) broadband interferometer has become an effective approach to map lightning channels in two dimensions with high time resolution. This paper reports an approach to mapping lightning channels in three dimensions (3D) using two simultaneous interferometers separated by about 10 km. A 3D mapping algorithm was developed based on the triangular intersection method considering the location accuracy of both interferometers and the arrival time of lightning VHF radiation. Simulation results reveal that the horizontal and vertical location errors within 10 km of the center of the two stations are less than 500 m and 700 m, respectively. The 3D development of an intra-cloud (IC) lightning ﬂash and a negative cloud-to-ground (-CG) lightning ﬂash with two different ground terminations in the same thunderstorm are reconstructed, and the extension direction and speed of lightning channels are estimated consequently. Both IC and CG ﬂash discharges showed a two-layer structure in the cloud with discharges occurring in the upper positive charge region and the lower negative charge region, and two horizontally separated positive charge regions were involved in the two ﬂashes. The average distance of the CG ground terminations between the interferometer results and the CG location system was about 448 m. Although disadvantages may still exist in 3D real-time location compared with the lightning mapping array system working with the principle of the time of arrival, interferometry with two or more stations has the advantage of lower station number and is feasible in regions with poor installation conditions, such as heavy-radio-frequency-noise regions or regions that are difﬁcult for the long-baseline location system.

Radiation in the very-high-frequency (VHF) range widely occur in various lightning processes and usually burst from the tip of the lightning breakdown [13][14][15].It can be used to map the lightning channel development with high resolution and is suitable for studying small-scale lightning-discharge processes.There are mainly two mapping methods to determine the position of a lightning VHF radiation source: the time-of-arrival (TOA) method [16][17][18] and the interferometry method [19,20].
The TOA method considers the VHF signal from a single lightning radiation source that arrives at radio sensors ranging from a few kilometers to tens of kilometers as a spherical wave.Based on the difference in arrival times at different sensors, the occurrence time and location of the radiation source can be estimated.This method is generally applied in a mapping network with five or more stations to obtain three-dimensional location results, e.g., the Lightning Mapping Array (LMA) [21,22].Thousands of radiation sources can be located for a single lightning flash with a temporal resolution of about tens of microseconds.The location accuracy mainly depends on the TOA estimation accuracy and the number and layout of stations [23].Usually, more stations can be fitted to obtain more accurate results due to the uncertainty in TOA estimation, and cover a wider detection area.However, it will also increase the complexity of the system and the requirement for time consistency of stations.
For the interferometry method, one station with three or more closely spaced sensors can provide mapping results.The distance between two sensors is typically a few meters to tens of meters [24].It can be regarded that the lightning radiation signal from about tens of kilometers arrives at an interferometer as a plane wave.The elevation and azimuth of the lightning radiation source relative to the interferometer can be determined as a two-dimensional location result by the cross-correlation analysis of VHF signals detected by the interferometer sensor array [20,25].Many interferometers have been designed and developed gradually.Compared with the TOA mapping system, the interferometer coherently detects the VHF signals of all sensors, allowing both isolated pulses and continuous radiation bursts to be mapped with higher time resolution and fewer stations.Commonly, the time resolution can reach a few microseconds.Breakthroughs in lightning science have also been achieved by using interferometry [26][27][28].
Due to the lack of distance information in the 2D mapping result, it is difficult to judge the development direction of lightning channels.Based on the previous understanding that the positive breakdown emits much weaker VHF radiation than the negative one and positive channels appear to have a more scattered development [29], the direction of the charge transfer can be determined in combination with the pulse polarity of the simultaneous electric field change.However, recent observations have shown that some positive breakdowns in virgin air can also produce intense VHF radiation, such as fast positive breakdown during lightning initiation [27,30], positive breakdown immediately after the negative return stroke [26], etc.Thus, there is uncertainty in judging the polarity of breakdown and the direction of the channel progressing.Moreover, uncertainty also exists in estimating the extension length and velocity of the lightning channel.Therefore, 3D mapping of the VHF interferometry is needed, and many methods have been proposed at present, such as synchronous observation of dual station interferometers [31][32][33][34], synchronous observation of an interferometer with a high-speed camera [35] or acoustic array [36], and the combination of 3D structures provided by the LMA and higher-timeresolution 2D mapping of an interferometer [27].The combination of different detection objects or different mapping technologies may have the problem of inconsistent temporal resolution and thus reduce the temporal and spatial resolution of the final 3D mapping results.Further attempts are still needed to deal with high-precision 3D mapping using only interferometers and the corresponding lightning discharge process.
A lightning VHF broadband interferometer with high time resolution has been developed since 2009 and has been applied in the SHAndong Triggering Lightning Experiment (SHATLE) in China [37][38][39].Synchronous observation of a two-station VHF broadband interferometer has been made since 2017, combined with simultaneous observations of fast/slow electric field changes.In this paper, the 3D lightning mapping method and the location system are first introduced, and the 3D mapping results of an intra-cloud (IC) lightning flash and a cloud-to-ground (CG) lightning flash with multiple-ground terminations are presented.

Lightning 3D Mapping Algorithm
The basic principle of 2D interferometry mapping is to determine the time differences (or phase differences) of the lightning radiation waveform arriving at the interferometer antenna array and obtain the incident wave vector of the radiation source relative to the interferometer, that is, the 2D position (elevation and azimuth) of the radiation source.Many 2D mapping algorithms for interferometers have been developed.Here, the algorithm presented in Stock and Krehbiel [40] is used for an interferometer with an arbitrary configuration.
The three-dimensional position can be resolved at the point of the intersection of two interferometer incident wave vectors from the same radiation source.Due to the location error of each interferometer, atmosphere refraction, and the asynchronous time error, two incident wave vectors may be on two different surfaces.A double theodolite measurement algorithm is used to calculate the 3D position of the radiation source.The radiation source is estimated on the common vertical line of two wave vectors.
Figure 1 shows the schematic diagram of the three-dimensional location method for two stations of interferometer O 1 (x 1 ,y 1 ,z 1 ) and O 2 (x 2 ,y 2 ,z 2 ).The unit incident vectors of the lightning radiation source point are l 1 and l 2 for two stations, respectively.
where cos α and cosβ represent the cosine of the incident angle related to the north-south baseline and east-west baseline, respectively.For simplicity, the three-dimensional spatial vector is defined as (u, v, w), where µ 2 + ν 2 + ω 2 = 1.Points M and N are intersection points of the common vertical line with two incident wave vectors.As the common vertical line MN is perpendicular to the l 1 and l 2 exactly, their dot products are zero.We introduce , ,  to express  ⃑ .From Equations ( 3), (4), and ( 7), x, y, z can be calculated by Equation (8), thus the common vertical line  ⃑ and its positio The distances between the intersection point and the corresponding interferometer station are l 1 and l 2 , respectively.The O 1 M, O 2 N and common vertical line MN can be written as We introduce (x, y, z) to express MN.From Equations ( 3), ( 4) and ( 7), x, y, z, l 1 , l 2 can be calculated by Equation (8), thus the common vertical line MN and its position relative to the coordinate origin can be obtained.
The position of the radiation source can be estimated to be located at point P on the common vertical line MN.In this paper, the weight coefficient ρ = is used to locate the position of point P in Equation ( 9), which considers the 2D location error of each interferometer and the distance between the interferometer and the corresponding intersection point with the common vertical line.Point P is close to the interferometer station with a low 2D location error and a short distance.σ 1 and σ 2 are 2D location errors of two stations, determined by the baseline length and the estimation error of the time difference of arrival.Here, in this case, we only consider the difference in baseline length for each station with the same time estimation error as a form of simplification.
Ideally, each interferometer can collect the lightning radiation source from the same lightning discharge process.Because of the different attenuation of the lightning radiation signal in the propagation path to each station, the time and geographical position error of the GPS, and the signal recording mode of each station, etc., not all the times of radiation sources of one station can correspond to that of another station.Considering the triangle inequality, that is, the time difference of the signal arriving at two stations cannot be greater than the propagation time of light through the distance between two stations, radiation sources within a certain time window are used to identify all possible matches of the same radiation source at both stations.The 3D location result is calculated for each pair of possible radiation sources, and the most reasonable result is the one with the closest time differences of arrival between two stations compared with the actual time difference between two stations.
The algorithm is given as follows: (1) Select one station as the main station and obtain all of the possible radiation source matches of each source of the main station within a certain time window considering the triangle inequality.The time delay of two 2D results of each station is estimated by the cross-correlation waveform.(2) Calculate the 3D location of all matches using Equations (8) and ( 9).
(3) Quality control is performed using the following four metrics for each possible 3D location: (a) Length of the common vertical line | MN| is less than a certain threshold.Here it is set to less than half of the minimum value of l 1 and l 2 .
(b) ∠MO 1 P and ∠NO 2 P are less than an empirical value of 10 • , that is, the difference between the incident vector of the 3D location point relative to each station and the incident vector provided by the 2D location of each station should be small.(c) Difference between the time delay of the calculated position to two stations and the actual time delay of two 2D results (denoted by DT) is less than a threshold, which is set to be 5 µs.

(d)
Since the antenna gain of each station is the same, the signal strength received by the station near the 3D location point should not be smaller than that of the farther one.
(4) Considering that there may be a radiation source satisfying multiple sets of matches, the matching pair with the smallest DT is selected as the most reasonable match in the time window.The corresponding 3D radiation source location is included in the final 3D mapping results.If none of the above quality controls are met, no 3D radiation source will be involved.( 5) Remove all 2D result matches involving the 3D radiation source obtained in Step 4,and repeat Step 4 to obtain another 3D radiation source location until the end of all matches.Because the detection efficiency varied for a source with different distances to two stations, there will be sources detected only by one station and only included in the 2D result.

Experiment and Data
Since 2017, a multi-station synchronous observation experiment of lightning VHF broadband interferometry was carried out in SHATLE.The distance between the two stations used in this paper was about 9.6 km.There were two VHF broadband interferometers with different bands of 140-300 MHz and 35-70 MHz in each station.Each interferometer had at least three identical broadband discone antennas to receive lightning VHF radiation signals and recorded signals in the continuous or sequential recording mode.Baselines ranged from 8 to 20 m.The maximum baseline of Station 2 was larger than that of Station 1.For the signal in the range of 140-300 MHz, the data sampling rate was 1 GS/s and the data vertical resolution was 12 bits.For the 35-70 MHz one, the data sampling rate and the data vertical resolution were 250 MS/s and 12 bits, respectively.The system could determine the location of lightning radiation sources in two dimensions (elevation and azimuth) by using a generalized cross-correlation time delay estimation algorithm and a parabolic interpolation algorithm [20,25].In this paper, 2D location results at 140-300 MHz are used for each lightning flash.Meanwhile, the fast antenna and slow antenna were used to measure the lightning electric field changes simultaneously.Their decay time constants are 0.1 ms and 0.22 s, and their bandwidths are 1.5 kHz-2 MHz and 10 Hz-1 MHz, respectively.The trigger time of each interferometer was provided by GPS with a high time accuracy of 25 ns for a cooperating analysis with other observations.

Error Analysis
The 3D locating error of two interferometers mainly depends on the incident vector error of the radiation source relative to each interferometer.The incident vector of the lightning radiation source is l = (µ, ν, ω) and has a relationship with the baseline l • d cτ, where d is the baseline of two antennas separated by a distance of d, τ is the time delay between two antennas, and c is the speed of light in air.To simplify the analysis, two non-collinear baselines in the same plane are considered to resolve the l .
Taking the coordinates µ in the east-west direction as an example, where α is the direction cosine between the incident plane wave and the east-west baseline of the interferometer antenna array, d 1x and d 1y are the projected distance of the baseline d 1 in the east-west and north-south directions, respectively, and the same for d 2x and d 2y for d 2 .Further, the coefficient of the time delay τ can be simplified to D µi and D vi , and i is the baseline index.Thus, as the time delay errors σ τi of each baseline are independent, the incident vector errors σ µ and σ µ are: The Monte Carlo method is used to study the error distribution of the 3D mapping on lightning using two separated VHF broadband interferometers in SHATLE.The simulation area is set in the range of 30 km × 30 km near the center of the two stations, and the grid size is 100 m × 100 m.The radiation sources are located at each grid point in three altitude layers of 1 km, 5 km, and 10 km.The ideal incident vectors of their radiated electromagnetic waves reaching each station are calculated, then Gaussian random noise with zero mean and standard deviation of 0.5 ns is superimposed as a time delay error.The time delay errors for each station and each baseline are independent of each other.The positions of the radiation sources are estimated using algorithms in Section 2.1, and finally, the mapping error is obtained by comparing the position of each grid point and the position estimated.The arithmetic mean values of mapping error at each altitude and each grid are obtained by repeating the above steps 1000 times.
Figure 2 shows the error distribution and the corresponding length of the common vertical line for sources at various altitudes.Both the horizontal and vertical error distributions are almost symmetrically distributed along the baseline direction, which is the line connecting the two stations.For sources at the same altitude, both the horizontal and vertical mapping errors are smallest in the direction perpendicular to the baseline and increase with the distance increasing to the center of the two stations.Also, the error increases as the source approaches the direction of the baseline.In the direction of the baseline, the error behind each station is the largest.
Within 10 km of the center of both stations, the horizontal mapping error at the same grid point is less than the vertical error.Except for radiation sources behind each station, the horizontal and vertical errors at different heights are less than 500 m and 700 m, respectively.Both horizontal and vertical errors decrease and then increase as the height of the radiation source increases, and have the smallest error distribution at a height of 5 km.For radiation sources behind each station, vertical errors show a similar tendency, while horizontal errors decrease with increasing height.The error distribution is different 10 km outside the center of both stations.The horizontal error is close to the vertical error or even greater than the vertical error for sources at the same grid point.Both the horizontal and vertical errors decrease with the increase in altitude.Therefore, the two-station interferometer observation has great advantages in the study of the close lightning process.
Similar to the distribution of horizontal and vertical mapping errors, the average lengths of the common vertical lines at different altitudes are symmetrically distributed along the baseline direction, while the lengths of the same horizontal grid points decrease with the increase in height.Contrary to the mapping error, the common vertical line of points close to the baseline direction is shorter than that of points perpendicular to the baseline direction, which suggests that it is not enough to judge the mapping quality only based on the length of the common vertical line, but the time difference between the actual time and the calculated time should also be considered.Within 10 km of the center of both stations, the horizontal mapping error at the sa grid point is less than the vertical error.Except for radiation sources behind each stati the horizontal and vertical errors at different heights are less than 500 m and 700 m, spectively.Both horizontal and vertical errors decrease and then increase as the heigh the radiation source increases, and have the smallest error distribution at a height of 5 k For radiation sources behind each station, vertical errors show a similar tendency, wh horizontal errors decrease with increasing height.The error distribution is different 10 outside the center of both stations.The horizontal error is close to the vertical error or ev greater than the vertical error for sources at the same grid point.Both the horizontal a

An IC Lightning Flash
Figure 3 shows the 2D mapping results of an IC lightning flash observed by two interferometers.The color of VHF radiation sources varies with time.This IC flash occurred on 30 July 2017 at about 18:24:00 UT.Due to the triggering setting, no simultaneous fast and slow electric field changes were recorded for this flash.As seen from 2D mapping results, the duration of the whole flash was about 1.1 s.
Figure 3 shows the 2D mapping results of an IC lightning flash observed by two interferometers.The color of VHF radiation sources varies with time.This IC flash occurred on 30 July 2017 at about 18:24:00 UT.Due to the triggering setting, no simultaneous fast and slow electric field changes were recorded for this flash.As seen from 2D mapping results, the duration of the whole flash was about 1.1 s.There were 9154 and 8322 points mapped for the IC flash at Station 1 and Station 2, respectively.This flash was distributed throughout the sky at Station 1, while it was mainly concentrated in the southeast at Station 2. The lightning discharge was closer to Station 1 and almost over it.Accordingly, there were more 2D mapping points at Station 1. From 2D mapping results, the overlap of discharge channels leads to difficulties in distinguishing the structure of discharge channels and thus in analyzing charge distribution involved in the discharge.
Figure 4 shows the 3D mapping result of the IC lightning flash.A total of 6104 3D points were obtained.This flash was mapped within 20 km of Station 1, which is consistent with the 2D mapping results in Figure 3.The 3D lightning structure and leader progression were well reconstructed.The IC flash showed a bilevel structure with vertical breakdown channels connecting the two layers.Considering that the negative breakdown radiation is stronger than the positive breakdown, and the associated negative channels tend to converge more, it can be judged that the upper layer was the positive charge region where negative breakdowns occurred, while the lower layer was the negative charge region where positive breakdowns occurred.The upper and lower layers were located at altitudes of 8-11 km and 4-7 km, respectively.There were 9154 and 8322 points mapped for the IC flash at Station 1 and Station 2, respectively.This flash was distributed throughout the sky at Station 1, while it was mainly concentrated in the southeast at Station 2. The lightning discharge was closer to Station 1 and almost over it.Accordingly, there were more 2D mapping points at Station 1. From 2D mapping results, the overlap of discharge channels leads to difficulties in distinguishing the structure of discharge channels and thus in analyzing charge distribution involved in the discharge.
Figure 4 shows the 3D mapping result of the IC lightning flash.A total of 6104 3D points were obtained.This flash was mapped within 20 km of Station 1, which is consistent with the 2D mapping results in Figure 3.The 3D lightning structure and leader progression were well reconstructed.The IC flash showed a bilevel structure with vertical breakdown channels connecting the two layers.Considering that the negative breakdown radiation is stronger than the positive breakdown, and the associated negative channels tend to converge more, it can be judged that the upper layer was the positive charge region where negative breakdowns occurred, while the lower layer was the negative charge region where positive breakdowns occurred.The upper and lower layers were located at altitudes of 8-11 km and 4-7 km, respectively.Comparing 3D and 2D mapping results, most discharge channels could be w mapped, such as the negative leader in the northeast direction and recoil leaders.Ho ever, for the positive channels in the dashed rectangle area in Figures 3 and 4, only a f points were obtained in the 3D result, whereas there was a fan-like channel extension the head of the positive leader in the 2D result of Station 1 in Figure 3a.Two factors mi related to the poor mapping result.First, this positive leader discharge occurred ab 20 km from Station 2 and near the baseline direction of two stations, corresponding t Comparing 3D and 2D mapping results, most discharge channels could be well mapped, such as the negative leader in the northeast direction and recoil leaders.However, for the positive channels in the dashed rectangle area in Figures 3 and 4, only a few points were obtained in the 3D result, whereas there was a fan-like channel extension at the head of the positive leader in the 2D result of Station 1 in Figure 3a.Two factors might be related to the poor mapping result.First, this positive leader discharge occurred about 20 km from Station 2 and near the baseline direction of two stations, corresponding to a low incident vector angle of the radiation source close to the horizon.According to the simulation results in Section 2.3, the mapping error in this area is indeed larger than in the other areas.The corresponding 3D results that could not meet the quality controls were excluded.Secondly, both interferometers used the triggering acquisition mode.As the positive leader was weak in the VHF radiation and far from Station 2, Station 2 recorded fewer data segments, resulting in fewer 2D mapping points.The mismatch in the number of 2D mapping points between two stations for the same discharge also reduced the ability to reconstruct the discharge channel.
The flash initiated at an altitude of 8.4 km and progressed vertically to an altitude of 10 km with a speed of 1.6 × 10 6 m/s.After that, the discharge channel changed to develop towards the southeast horizontally for about 20 km, gradually increasing in height, shown as negative leader1 in Figure 4d.About 0.1 ms after the flash initiation, the first sources associated with positive breakdown were located at 7.3 km altitude, and then a few radiation sources occurred in the negative charge region.During the preliminary stage of IC flash, the development of positive and negative leaders was mainly concentrated in the southeast direction of the origin of the flash, and positive leader discharges were basically in the area directly below the negative leader discharges.
As shown in Figure 5, during the late stage of the IC flash, a recoil leader initiated from the lower negative charge region at about 1.5 s and propagated along the previous discharge channel to the upper positive charge region with an average speed of 8.6 × 10 6 m/s.About 1 ms later, the recoil leader reached the tip of a previously negative leader branch, and developed into a new negative leader in the virgin air.It first extended horizontally towards the northeast at 10-11 km altitude.Then, it suddenly dropped to a height of 8 km in 0.1 ms and kept developing horizontally towards the northeast with plenty of branches.Also, there was an increase in the extension speed from 4.6 × 10 5 m/s to 8.6 × 10 5 m/s compared with that at the higher altitude.
Figure 6 shows the IC flash in a time-distance representation, with the initiation point of the flash as a reference.The color changes with the altitude.The slope can reflect the average radial speed of different leaders.The difference in speed and altitude between positive and negative leaders can be seen directly.The negative leader discharge of the late stage developed at an average rate of about 10 5 m/s, which was closer to that of the negative leader during the early stage, but with more branches, a longer extension, and a lower altitude.It indicates that two regions of positive polarity charge at different altitudes might be involved in this leader discharge successively.
With the development of the negative leader, many scattered and transient recoil leaders were also detected in the lower negative charge region, where there was no discharge before.They appeared gradually further away from the initiation position with time at an average extension speed closer to the reference of 2 × 10 4 m/s, indicating that positive leaders were developing.It suggests that the discharge of the upper negative leader activated the lower negative charge region and influenced the positive leader activity, while this correlation is difficult to be clarified using only the 2D mapping results due to the lack of vertical distribution information.

A CG Lightning Flash
A negative CG lightning flash with multiple ground terminations occurred at about 18:25:41 UT, about 1.6 s after the above IC flash.Figure 7 shows the 2D mapping results observed by two interferometers.This negative flash lasted about 0.8 s and consisted of four negative leader/return stroke sequences.The time intervals between two adjacent return strokes were 28.8 ms, 65.3 ms, and 37.9 ms, respectively.The CG flash contacted the ground by two ground terminations, with the first three return strokes having the same ground termination.Like the IC flash, this CG flash occurred above the interferometer Station 1 and southeast of Station 2. In total, 13,717 and 7961 points were mapped by Station 1 and Station 2, respectively.At Station 2, mapping results below the lightning origin point were compressed to about 25 degrees from the horizon, bringing a poorer resolution of the downward negative leader than the others.m/s.About 1 ms later, the recoil leader reached the tip of a previously negative lead branch, and developed into a new negative leader in the virgin air.It first extended ho zontally towards the northeast at 10-11 km altitude.Then, it suddenly dropped to a heig of 8 km in 0.1 ms and kept developing horizontally towards the northeast with plenty branches.Also, there was an increase in the extension speed from 4.6 × 10 5 m/s to 8.6 × 1 m/s compared with that at the higher altitude.average radial speed of different leaders.The difference in speed and altitude between positive and negative leaders can be seen directly.The negative leader discharge of the late stage developed at an average rate of about 10 5 m/s, which was closer to that of the negative leader during the early stage, but with more branches, a longer extension, and a lower altitude.It indicates that two regions of positive polarity charge at different altitudes might be involved in this leader discharge successively.With the development of the negative leader, many scattered and transient recoil leaders were also detected in the lower negative charge region, where there was no discharge before.They appeared gradually further away from the initiation position with time at an average extension speed closer to the reference of 2 × 10 4 m/s, indicating that positive leaders were developing.It suggests that the discharge of the upper negative leader activated the lower negative charge region and influenced the positive leader activity, while this correlation is difficult to be clarified using only the 2D mapping results due to the lack of vertical distribution information.

A CG Lightning Flash
A negative CG lightning flash with multiple ground terminations occurred at about 18:25:41 UT, about 1.6 s after the above IC flash.Figure 7 shows the 2D mapping results observed by two interferometers.This negative flash lasted about 0.8 s and consisted of four negative leader/return stroke sequences.The time intervals between two adjacent return strokes were 28.8 ms, 65.3 ms, and 37.9 ms, respectively.The CG flash contacted the ground by two ground terminations, with the first three return strokes having the same ground termination.Like the IC flash, this CG flash occurred above the interferometer Station 1 and southeast of Station 2. In total, 13,717 and 7961 points were mapped by Station 1 and Station 2, respectively.At Station 2, mapping results below the lightning origin point were compressed to about 25 degrees from the horizon, bringing a poorer resolution of the downward negative leader than the others.Figure 8 shows the 3D mapping result of the CG flash.A total of 5993 3D points were obtained.The 3D lightning structure and leader progression were also well reconstructed.Compared with the IC flash 1.6 s ago, the overall spatial scale of the CG flash was smaller and more concentrated in a 10 km area around Station 1.The polarity of the discharge can be determined according to the discharge characteristics of different polarity radiation, the development direction of the 3D channel, and the polarity of the corresponding elec- Figure 8 shows the 3D mapping result of the CG flash.A total of 5993 3D points were obtained.The 3D lightning structure and leader progression were also well reconstructed.
Compared with the IC flash 1.6 s ago, the overall spatial scale of the CG flash was smaller and more concentrated in a 10 km area around Station 1.The polarity of the discharge can be determined according to the discharge characteristics of different polarity radiation, the development direction of the 3D channel, and the polarity of the corresponding electric field change [24].Similar to the previous IC flash, the in-cloud structure of the CG flash also showed a double-layer structure with negative discharge channels in the upper positive charge region and positive discharge channels in the lower negative charge region.Altitudes of upper and lower layers were 8-11 km and 4-6 km, respectively.As shown in Figure 9, the flash initiated at about 4.3 km with obvious preliminar breakdown pulses in the fast electrical field.Different from the IC flash, the initiation poin was closer to the lower negative charge region.The leader developed upward obliquely to 10 km altitude with several branches and continued to propagate horizontally in thre directions for about 140 ms.Multiple recoil leaders propagated along the preexisting neg ative preliminary breakdown channels in the positive charge region before the negativ downward stepped leader occurred.A few radiation sources were located in the negativ charge region at a 4 km altitude several milliseconds after the lightning initiation.Mean time, some recoil leaders in the negative charge region extended the positive charge chan nels in a retrograde manner, and then many sources scattered at the starting points o As shown in Figure 9, the flash initiated at about 4.3 km with obvious preliminary breakdown pulses in the fast electrical field.Different from the IC flash, the initiation point was closer to the lower negative charge region.The leader developed upward obliquely to 10 km altitude with several branches and continued to propagate horizontally in three directions for about 140 ms.Multiple recoil leaders propagated along the preexisting negative preliminary breakdown channels in the positive charge region before the negative downward stepped leader occurred.A few radiation sources were located in the negative charge region at a 4 km altitude several milliseconds after the lightning initiation.Meantime, some recoil leaders in the negative charge region extended the positive charge channels in a retrograde manner, and then many sources scattered at the starting points of recoil leaders.The average speed of recoil leaders was about 10 7 m/s.Figure 10 shows the 3D mapping result and electric field of the leader/return strok sequences and inside-cloud discharge after the first return stroke.About 2 ms after fre quent recoil leaders along the preliminary breakdown channels, the negative downwar stepped leader channel started at about 2.8 km altitude and approached the ground wit three main branches extending discharge channels simultaneously, which could only b seen from the near interferometer.About 197 ms later, one leader finally connected to th ground and induced the return stroke, and several recoil leaders occurred following th return stroke.Two types of recoil leaders could be recognized: some carried negativ charges propagating from the preexisting positive leader channel to the flash initiatio point, and some carried positive charges propagating from the charge region of the nega tive leader to space.Figure 10 shows the 3D mapping result and electric field of the leader/return stroke sequences and inside-cloud discharge after the first return stroke.About 2 ms after frequent recoil leaders along the preliminary breakdown channels, the negative downward stepped leader channel started at about 2.8 km altitude and approached the ground with three main branches extending discharge channels simultaneously, which could only be seen from the near interferometer.About 197 ms later, one leader finally connected to the ground and induced the return stroke, and several recoil leaders occurred following the return stroke.Two types of recoil leaders could be recognized: some carried negative charges propagating from the preexisting positive leader channel to the flash initiation point, and some carried positive charges propagating from the charge region of the negative leader to space.The dart leader causing the second return stroke initiated not from the starting poin of the first leader, but from a previous channel in the negative charge region.It propagated down to the flash initiation region and then traveled along the same channel of the firs leader-stroke sequence.The leader exhibited chaotic pulse trains (CPTs) with significan slow electric field variations accompanying strong VHF radiation signals.The result indi cated intense breakdown discharges during the downward propagation of the CPT leader.The peak amplitude of the second return stroke was 6.8 times smaller than that o the first.
The leader preceding the last leader-return stroke process initiated from the startin region of the flash and propagated downward to the ground as a stepped leader.Durin the development of the stepped leader, there were sources detected along the previou return stroke channel in the 2D location result toward the ground.Then a weak return stroke pulse in the electric field waveform indicated a third return stroke occurring using The dart leader causing the second return stroke initiated not from the starting point of the first leader, but from a previous channel in the negative charge region.It propagated down to the flash initiation region and then traveled along the same channel of the first leader-stroke sequence.The leader exhibited chaotic pulse trains (CPTs) with significant slow electric field variations accompanying strong VHF radiation signals.The result indicated intense breakdown discharges during the downward propagation of the CPTs leader.The peak amplitude of the second return stroke was 6.8 times smaller than that of the first.
The leader preceding the last leader-return stroke process initiated from the starting region of the flash and propagated downward to the ground as a stepped leader.During the development of the stepped leader, there were sources detected along the previous return stroke channel in the 2D location result toward the ground.Then a weak return stroke pulse in the electric field waveform indicated a third return stroke occurring using the same termination as the first two strokes.In the meantime, the stepped leader kept propagating with plenty of small branches and finally induced the return stroke.The field peaks of the last two strokes were 14.9 and 10.2 times smaller than that of the first stroke.
The terminations in the first three leader-return stroke sequences and the fourth one were different due to the newly established stepped leader path preceding the fourth return stroke.The CGLLS (Cloud-to-Ground Lightning Location System) located three return strokes of the CG except the weakest third one.The average position of locations during 1 ms before the return stroke was chosen as the ground termination given by the interferometer, and distances between the CGLLS results and the interferometer results for three ground terminations were about 341 m, 590 m, and 413 m, with an average of 448 m.Although there is a location error of about 710 m in CGLLS [41], it can reflect the accuracy of the interferometer 3D results to some extent.
After the last leader-return stroke process, there was a long period of inside-cloud discharge.A new negative leader was detected to start from the channel of the last downward stepped leader and extended upward obliquely to 7.9 km altitude with several branches.The corresponding positive charge region was different from that of the preliminary breakdown process.Two horizontally separated positive charge regions were involved in the discharge of the CG flash.

Summary
Two stations' synchronous observation of the VHF broadband interferometer has been conducted.Three-dimensional radiation source mapping results of an intra-cloud lightning flash and a negative cloud-to-ground flash indicated the location system could effectively reconstruct the lightning discharge channel and visually describe the spatiotemporal development characteristics of the lightning discharge process.In this paper, the source position in the common vertical line depended on the distance from the source to two stations and the baseline length of each station.As the location accuracies of elevation and azimuth are differently determined by the elevation angle, in the future the elevation angle will be considered in the weight coefficient ρ to calculate the source position near the common vertical line.
The development of both IC and CG flash discharges inside the cloud appeared to be a bilevel structure in the upper positive and lower negative charge regions, and two horizontally separated positive charge regions were involved in both flashes.The discharge channel of the preliminary breakdown process sloped upward and connected the bilevel structure.Recoil leaders occurred in almost all the stages of the whole flash.The speeds of recoil leaders were about 10 7 m/s.They propagated along the preceding electrified channels not only in the upper positive charge region but also in the lower negative charge region.For the CG flash, multiple terminations were produced by two different stepped leaders with different negative charge regions.
According to the Monte Carlo simulation results and the three-dimensional mapping results of the two lightning flashes, the mapping accuracy is affected by the 2D mapping error of each interferometer and the position of the discharge channel relative to two interferometers, while the mapping efficiency can also be limited by the strength and distance of the discharge signal and the simultaneous observation efficiency of two interferometers.It can be concluded that for the lightning discharge relatively close to the interferometers, 3D mapping using two VHF broadband interferometers has the better mapping effect.Currently, the VHF broadband interferometer uses full waveform information for hightemporal-resolution mapping and therefore applies post hoc rather than real-time location.Although disadvantages may still exist in the wide-area and real-time 3D location compared with the long-baseline lightning-mapping array system working with the principle of the time of arrival, 3D mapping interferometry with two or more stations has the advantages of fewer station numbers and is feasible in regions with poor installation conditions, such as heavy-radio-frequency-noise regions or mountainous regions.
Though the efficiency of 3D mapping for two interferometers could not be comparable to that of 2D mapping, 3D mapping can supply additional information on the temporal and spatial development of discharge channels, such as 3D morphology, extension direction, and speed.Moreover, the overall structure of the lightning can be reproduced in fine detail.It is possible to improve the capability of high-temporal-resolution 2D mapping interferometry to analyze the details of lightning discharge and distribution of the thunderstorm charges involved in the discharge.As the number of interferometers increases, the time-of-arrival information of each VHF signal can also be used to improve the efficiency of 3D mapping.

Figure 1 .
Figure 1.Schematic diagram of the three-dimensional location method.

Figure 3 .
Figure 3. Two-dimensional mapping results of an intra-cloud lightning flash.Cosine projection (a) and elevation with time (c) observed by Station 1. (b,d) are the same as (a,c), but for Station 2. The white + shows the zenith of each interferometer, and the circle is the horizon.Dashed lines intersect at the origin of the flash.

Figure 3 .
Figure 3. Two-dimensional mapping results of an intra-cloud lightning flash.Cosine projection (a) and elevation with time (c) observed by Station 1. (b,d) are the same as (a,c), but for Station 2. The white + shows the zenith of each interferometer, and the circle is the horizon.Dashed lines intersect at the origin of the flash.

Figure 4 .
Figure 4. Three-dimensional location results of an intra-cloud lightning flash.(a) height with ti (b) north-south vertical projection, (c) height distribution of the number of radiation sources, plane projection, and (e) east-west vertical projection of the three-dimensional location result.terisks are positions of the interferometer.Dashed lines intersect at the initiation point of the fla

Figure 4 .
Figure 4. Three-dimensional location results of an intra-cloud lightning flash.(a) height with time, (b) north-south vertical projection, (c) height distribution of the number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the three-dimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect at the initiation point of the flash.

Figure 5 .
Figure 5. Three-dimensional location results of the late stage of the IC flash in Figure 4. (a) heig with time, (b) north-south vertical projection, (c) height distribution of the number of radiat sources, (d) plane projection, and (e) east-west vertical projection of the three-dimensional locat

Figure 5 .
Figure 5. Three-dimensional location results of the late stage of the IC flash in Figure 4. (a) height with time, (b) north-south vertical projection, (c) height distribution of the number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the three-dimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect at the initiation point of the flash.

Figure 6 .
Figure 6.Time-distance graphs of sources of the IC flash.The solid line indicates the reference location for the distance, which is the initiation point of the IC flash.The sign of distance is consistent with the sign of the difference in height between the radiation source and the lightning initiation point.The dashed reference lines indicate slopes corresponding to speeds of 2 × 10 4 m/s, 1 × 10 5 m/s, and 1 × 10 6 m/s.

Figure 6 .
Figure 6.Time-distance graphs of sources of the IC flash.The solid line indicates the reference location for the distance, which is the initiation point of the IC flash.The sign of distance is consistent with the sign of the difference in height between the radiation source and the lightning initiation point.The dashed reference lines indicate slopes corresponding to speeds of 2 × 10 4 m/s, 1 × 10 5 m/s, and 1 × 10 6 m/s.Remote Sens. 2022, 14, x FOR PEER REVIEW 12 of 19

Figure 7 .
Figure 7. Two-dimensional mapping results of a cloud-to-ground lightning flash.Cosine projection (a) and elevation with time (c) observed by Station 1. (b,d) are the same as (a,c), but for Station 2. The white + shows the zenith of each interferometer, and the circle is the horizon.Dashed lines intersect at the origin of the flash.

Figure 7 .
Figure 7. Two-dimensional mapping results of a cloud-to-ground lightning flash.Cosine projection (a) and elevation with time (c) observed by Station 1. (b,d) are the same as (a,c), but for Station 2. The white + shows the zenith of each interferometer, and the circle is the horizon.Dashed lines intersect at the origin of the flash.

1 Figure 8 .
Figure 8. Three-dimensional location results of the CG lightning flash.(a) height with time, (b north-south vertical projection, (c) height distribution of the number of radiation sources, (d) plan projection, and (e) east-west vertical projection of the three-dimensional location result.Asterisk are positions of the interferometer.Dashed lines intersect at the initiation point of the flash.

Figure 8 .
Figure 8. Three-dimensional location results of the CG lightning flash.(a) height with time, (b) northsouth vertical projection, (c) height distribution of the number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the three-dimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect at the initiation point of the flash.

1 Figure 9 .
Figure 9. Three-dimensional location results of the preliminary breakdown process of the CG flas in Figure 8.(a) height with time, (b) north-south vertical projection, (c) height distribution of th number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the three dimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect a the initiation point of the flash.

Figure 9 .
Figure 9. Three-dimensional location results of the preliminary breakdown process of the CG flash in Figure 8.(a) height with time, (b) north-south vertical projection, (c) height distribution of the number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the threedimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect at the initiation point of the flash.

1 Figure 10 .
Figure 10.Three-dimensional location results of the leader-return stroke sequences of the CG flas in Figure 8.(a) height with time, (b) north-south vertical projection, (c) height distribution of th number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the three dimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect a the initiation point of the flash.

Figure 10 .
Figure 10.Three-dimensional location results of the leader-return stroke sequences of the CG flash in Figure 8.(a) height with time, (b) north-south vertical projection, (c) height distribution of the number of radiation sources, (d) plane projection, and (e) east-west vertical projection of the threedimensional location result.Asterisks are positions of the interferometer.Dashed lines intersect at the initiation point of the flash.