A 3D Interferometer-Type Lightning Mapping Array for Observation of Winter Lightning in Japan

: We have developed and deployed a 3D Interferometer-type Lightning Mapping Array (InLMA) for observing winter lightning in Japan. InLMA consists of three broadband interferometers installed at three stations with a distance from 3 to 5 km. At each interferometer station, three discone antennas were installed, forming a right triangle with a separation of 75 m along their two orthogonal baselines. The output of each InLMA antenna is passed through a 400 MHz low-pass ﬁlter and then recorded at 1 GS/s with 16-bit accuracy. A new method has been proposed for ﬁnding 3D solutions of a lightning mapping system that consists of multiple interferometers. Using the InLMA, we have succeeded in mapping a positive cloud-to-ground (CG) lightning ﬂash in winter, particularly its preliminary breakdown (PB) process. A study on individual PB pulse processes allows us to infer that each PB pulse process contains many small-scale discharges scattering in a height range of about 150 m. These small-scale discharges in a series of PB pulses appear to be continuous in space, though discontinuous in time. We have also examined the positive return stroke in the CG ﬂash and found a 3D average return stroke speed of 7.5 × 10 7 m/s.


Introduction
As demonstrated by Wang et al. [1] and Zheng et al. [2], winter thunderstorms in Japan could exhibit various charge structures featured with low altitudes and small separations between charge regions of opposite polarities.Under such charge structures, many lightning flashes tend to have very complicated characteristics, as shown in a case study by Huang et al. [3].In order to understand such complicated lightning, it is necessary to track their progressions in both time and space.Recently, we have developed a 3D lightning mapping system called DALMA (Discone Antenna Lightning Mapping Array) [4] for observing the 3D progressions of winter lightning in Japan.DALMA is based on the TOA technique and works at MF-HF bands.As shown by Wang et al. [4], for each of the two examples of lightning flashes, more than three thousand sources are successfully located with the DALMA.However, as is expected, DALMA can only locate pulse discharges with sufficient MF-HF radiations.This significantly limits the time resolution of DALMA.Therefore, many lightning discharge processes, such as preliminary breakdowns and return strokes, cannot be studied in detail with DALMA.To overcome these problems, we have developed and deployed an interferometer-type lightning mapping array in the area covered by DALMA, and we have operated a supplementary system to DALMA for studying detailed lightning discharge processes in a small area.
As reviewed by Stock et al. [5], interferometers have been used in lightning mapping for more than 40 years.Historically, interferometers used narrowband radio signals of lightning through hardware interferometry techniques [6,7], but nowadays almost all Remote Sens. 2023, 15,1923 2 of 16 interferometers utilize broadband signals through software interferometry approaches pioneered by Shao et al. [8], Ushio et al. [9], Akita et al. [10], and Stock et al. [5].Recently, 3D interferometers with two stations have been developed and tested [10][11][12][13].In such systems, 3D mapping is usually achieved using the triangulation method by combining azimuth and elevation angles from the two stations, as detailed by Thyer [14] and Liu et al. [15].A similar method was used by Mardiana et al. [16].Recently, Jensen et al. [13] have included timing information for better matching of the azimuth and elevation angles from two stations.The 3D interferometer system developed in this study consists of three stations.For the 3D mapping, instead of using the triangulation method, we have proposed a new method that can be easily used in a system with multiple stations, which, of course, includes the two station systems already reported in the literature [10][11][12][13].Furthermore, as discussed by Stock et al. [5], the performance of an interferometer is related to many factors, such as noise level, sampling rate, bit number of a digitizer, and baseline length.In order to achieve high performance, we have used a digitizer that operates at 16 bits and has a maximum sampling rate of 1 GHz.The baseline used in this study is 75 m, which is more than twice longer than those used in previous 3D interferometer systems.This paper aims to report the instrumentation, the method, and the preliminary results of our new system.

Instrumentation and Methods
Figure 1 shows the schematic of the interferometer hardware at one station in this study.The outputs of three discone antennas are connected to three channels of a Gage PCI Express CompuScope digitizer through three 400 MHz low-pass filters, respectively.The 1 PPS output of a GPS receiver named FURY-10M is connected to the 4th channel of the digitizer.The rising edge of this 1 PPS signal is used to synchronize different stations.The digitizer is disciplined using the GPS 10 MHz output, so the synchronization accuracy should be better than 100 ns.We use a compact industrial computer to control the digitizer.The digitizer has a 4 GS/8 GB onboard sample memory.For each trigger, the 4-channel data, each sampled at 1 GHz for 0.8 s, are first recorded on this memory and are then transferred to the hard disk of the computer with a transit time of about 10 s.The digitizer can be operated in pre-trigger mode.However, for every channel, only 8192 sampling points are available for pre-trigger memory.The discone antennas that we used are called Super Discone Antenna D1300AM, which is the same for DALMA, and receives electromagnetic radiation in the frequency range from less than 1 MHz to over 1300 MHz.Three discone antennas were separated at a distance of 75 m, forming a right triangle.
As reviewed by Stock et al. [5], interferometers have been used in lightning mapping for more than 40 years.Historically, interferometers used narrowband radio signals of lightning through hardware interferometry techniques [6,7], but nowadays almost all interferometers utilize broadband signals through software interferometry approaches pioneered by Shao et al. [8], Ushio et al. [9], Akita et al. [10], and Stock et al. [5].Recently, 3D interferometers with two stations have been developed and tested [10][11][12][13].In such systems, 3D mapping is usually achieved using the triangulation method by combining azimuth and elevation angles from the two stations, as detailed by Thyer [14] and Liu et al. [15].A similar method was used by Mardiana et al. [16].Recently, Jensen et al. [13] have included timing information for better matching of the azimuth and elevation angles from two stations.The 3D interferometer system developed in this study consists of three stations.For the 3D mapping, instead of using the triangulation method, we have proposed a new method that can be easily used in a system with multiple stations, which, of course, includes the two station systems already reported in the literature [10][11][12][13].Furthermore, as discussed by Stock et al. [5], the performance of an interferometer is related to many factors, such as noise level, sampling rate, bit number of a digitizer, and baseline length.
In order to achieve high performance, we have used a digitizer that operates at 16 bits and has a maximum sampling rate of 1 GHz.The baseline used in this study is 75 m, which is more than twice longer than those used in previous 3D interferometer systems.This paper aims to report the instrumentation, the method, and the preliminary results of our new system.

Instrumentation and Methods
Figure 1 shows the schematic of the interferometer hardware at one station in this study.The outputs of three discone antennas are connected to three channels of a Gage PCI Express CompuScope digitizer through three 400 MHz low-pass filters, respectively.The 1 PPS output of a GPS receiver named FURY-10M is connected to the 4th channel of the digitizer.The rising edge of this 1 PPS signal is used to synchronize different stations.The digitizer is disciplined using the GPS 10 MHz output, so the synchronization accuracy should be better than 100 ns.We use a compact industrial computer to control the digitizer.The digitizer has a 4 GS/8 GB onboard sample memory.For each trigger, the 4-channel data, each sampled at 1 GHz for 0.8 s, are first recorded on this memory and are then transferred to the hard disk of the computer with a transit time of about 10 s.The digitizer can be operated in pre-trigger mode.However, for every channel, only 8192 sampling points are available for pre-trigger memory.The discone antennas that we used are called Super Discone Antenna D1300AM, which is the same for DALMA, and receives electromagnetic radiation in the frequency range from less than 1 MHz to over 1300 MHz.Three discone antennas were separated at a distance of 75 m, forming a right triangle.We have equipped 3 identical interferometers at 3 stations separated with a distance of 3-5 km.As a system, these 3 interferometers are called InLMA, short for Interferometer-type Lightning Mapping Array. Figure 2 shows the site map of InLMA.All 3 InLMA stations are set up near a windmill in the Uchinada town of Ishikawa, Japan, where we have been carrying out various lightning observation experiments for more than 17 years [17].In our recent observations, we also deployed FALMA [18,19] (fast antenna lightning mapping array) and DALMA.FALMA and DALMA sites are also included in Figure 2. FALMA covers a wider area than DALMA, while DALMA covers a wider area than InLMA.
We have equipped 3 identical interferometers at 3 stations separated with a distance of 3-5 km.As a system, these 3 interferometers are called InLMA, short for Interferometertype Lightning Mapping Array. Figure 2 shows the site map of InLMA.All 3 InLMA stations are set up near a windmill in the Uchinada town of Ishikawa, Japan, where we have been carrying out various lightning observation experiments for more than 17 years [17].In our recent observations, we also deployed FALMA [18,19] (fast antenna lightning mapping array) and DALMA.FALMA and DALMA sites are also included in Figure 2. FALMA covers a wider area than DALMA, while DALMA covers a wider area than InLMA.For a signal source whose geometric size is much smaller than its distance to the antennas of InLMA, it can be regarded as a point with its coordinates of , ,  as shown in Figure 3.We assume that this signal arrives at 3 stations with elevations  ,  , and  and azimuths  ,  , and  , respectively.At each station, the azimuth and the elevation can be calculated using the method outlined by M. G. Stock et al. [5], among others [20][21][22].Let us take station 1 as the example.Assume that the two distances between the antennas (ANT1, ANT2, and ANT3) along the two baselines of station 1 are  and  , the geometry angles between the signal incident wave and the two baselines are  and , and the corresponding time arrival differences of the signal are  and  .In this study,  and  are obtained using cross correlation between the signals from 2 corresponding antennas, according to the geometric relationship in Figure 3, For a signal source whose geometric size is much smaller than its distance to the antennas of InLMA, it can be regarded as a point with its coordinates of (x, y, z) as shown in Figure 3.We assume that this signal arrives at 3 stations with elevations el 1 , el 2 , and el 3 and azimuths az 1 , az 2 , and az 3 , respectively.At each station, the azimuth and the elevation can be calculated using the method outlined by M. G. Stock et al. [5], among others [20][21][22].Let us take station 1 as the example.Assume that the two distances between the antennas (ANT1, ANT2, and ANT3) along the two baselines of station 1 are d 12 and d 13 , the geometry angles between the signal incident wave and the two baselines are α and β, and the corresponding time arrival differences of the signal are τ 12 and τ 13 .In this study, τ 12 and τ 13 are obtained using cross correlation between the signals from 2 corresponding antennas, according to the geometric relationship in Figure 3, where c is the speed of light.According to geometry, cosine α and β can be substituted with the azimuth az 1 and elevation el 1 of the incident source.cosα = cos(az 1 )cos(el 1 ), ,  , ∆ , ∆ described above.
where  is the total number of stations ( 3 in this study),  ,  ,  are random measurement errors of elevation, azimuth, and time, respectively.Combining ( 1), ( 2), ( 3), ( 4), the azimuth and elevation can be obtained as follows.
Due to that arctan(x) ∈ − π 2 , π 2 and elevation is limited to positive after being inverted, Equations ( 5) and ( 6) cannot include all the possibilities of direction.Consequently, according to the value of cosβ, azimuth should be adjusted sometimes.If cosβ < 0 • , azimuth should be 180 • larger.In this way, for 2D mapping, azimuth is limited to − π 2 , 3π 2 and elevation is limited to 0, π 2 .Similarly, we can also get el 2 , az 2 , and el 3 and az 3 for stations 2 and 3.Meanwhile, since all 3 stations are equipped with GPS, we can also get the time differences ∆t 1j between that for a source to arrive at station 1 and the remaining two stations (j = 2, 3).
For a given source (x, y, z) discharge occurring at time t, according to geometry, we can find its azimuths az i (x, y, z) and elevations el i (x, y, z) of all stations, as well as the time differences ∆t 1j (x, y, z) between that for a source to arrive at station 1 and the remaining 2 stations (j = 2, 3).In this study, we propose the use of the following χ 2 equation to find the best solution of (x, y, z) for the set of the corresponding az 1 , el 1 , az 2 , el 2 , az 3 , el 3 , ∆t 12 , ∆t 13 described above.
where N is the total number of stations (N = 3 in this study), σ i el , σ i az , σ i t are random measurement errors of elevation, azimuth, and time, respectively.
The flow chart for realizing 3D location is shown in Figure 4a.First, we divided the data recorded at each channel of the 3 stations into a series of uniform windows based on the GPS PPS signal as the time reference.Secondly, using the window data and cross correlation, we found the time differences between the antennas of each station.Thirdly, from these time differences, we got all the azimuths and elevations.Fourthly, by matching windows of different stations and cross correlation, we found the time differences between different stations.The method for matching the time windows and then finding the time differences between the 2 stations is shown in Figure 4b.The red triangles are used to indicate the same discharge event observed at station 1 and station 2. Apparently, the time difference consists of two parts: ∆t w 12 and ∆t c 12 .The first part is the time difference, which can be measured in time windows (window 1 and window 3 in this case), and the second part is the time difference bounded by the time window and can be found by the cross correlation method.Finally, the Least-square method using the χ 2 equation in ( 7) is used not only for finding the best solution of (x, y, z), but also for matching the windows.If there are any windows used multiple times in previous solutions, only the solution with the smallest χ 2 is reserved.
radiates electromagnetic waves to the three stations.The corresponding angles are denoted and used to find the 3D location of the source , ,  .
The flow chart for realizing 3D location is shown in Figure 4a.First, we divided the data recorded at each channel of the 3 stations into a series of uniform windows based on the GPS PPS signal as the time reference.Secondly, using the window data and cross correlation, we found the time differences between the antennas of each station.Thirdly, from these time differences, we got all the azimuths and elevations.Fourthly, by matching windows of different stations and cross correlation, we found the time differences between different stations.The method for matching the time windows and then finding the time differences between the 2 stations is shown in Figure 4b.The red triangles are used to indicate the same discharge event observed at station 1 and station 2. Apparently, the time difference consists of two parts: ∆ and ∆ .The first part is the time difference, which can be measured in time windows (window 1 and window 3 in this case), and the second part is the time difference bounded by the time window and can be found by the cross correlation method.Finally, the Least-square method using the  equation in (7) is used not only for finding the best solution of , ,  , but also for matching the windows.If there are any windows used multiple times in previous solutions, only the solution with the smallest  is reserved.

Observation and Data
We finished deployment of the InLMA and started the observations of winter lightning in the Hokuriku area of Japan at the end of November of 2021.Since then, more than 100 lightning flashes have been recorded by InLMA.However, during December of 2021, one station of the InLMA was damaged by the voltage surge caused by nearby lightning.Only one of the lightning flashes, which occurred at 20:27:52 on 21 December 2021, was recorded by three stations of InLMA.The electric field change of the lightning flash recorded at one station of FALMA is shown in Figure 5.The "atmospheric electricity" sign convention is used.From the expanded waveforms of the E-field change, we know that it is a positive cloud-to-ground (CG) flash, which consists of two stages of the characteristic PB (preliminary breakdown) pulses (Figure 5b,d,e) and a single positive return stroke (Figure 5c,f).Based on Figure 5f, the peak current of the return stroke is estimated to be about 38 kA [23].Figure 6 shows its 3D mapping through DALMA.A total of 881 sources were located by DALMA.According to this mapping, the lightning flash started from a height of about 2.4 km MSL at a distance of about 15 km away from the central position of InLMA.The main sources mapped by DALMA occurred in the first 150 ms, which agrees with the waveform shown in Figure 5.An initial negative leader was observed to propagate upward and then divide into two main parts.The return stroke happened about 62.6 ms after lightning initiation.As shown in Figure 6, the lightning discharge sources spread into an area with a horizontal scale of several tens of kilometers, while only a few sources appeared above InLMA.
An example waveform of the corresponding raw data of InLMA along with the spectrums of the main lightning signal and background noise, separately denoted, is shown in Figure 7.The lightning signal mainly appears in the first 150 ms with a frequency bandwidth from about 1 MHz to more than 300 MHz.Previously, we have tried various frequency bands for the mapping and found some small differences in the mapping results [24], but in this study, we only use the frequency band from 80 MHz to 230 MHz.

Observation and Data
We finished deployment of the InLMA and started the observations of winter lightning in the Hokuriku area of Japan at the end of November of 2021.Since then, more than 100 lightning flashes have been recorded by InLMA.However, during December of 2021, one station of the InLMA was damaged by the voltage surge caused by nearby lightning.Only one of the lightning flashes, which occurred at 20:27:52 on 21 December 2021, was recorded by three stations of InLMA.The electric field change of the lightning flash recorded at one station of FALMA is shown in Figure 5.The "atmospheric electricity" sign convention is used.From the expanded waveforms of the E-field change, we know that it is a positive cloud-to-ground (CG) flash, which consists of two stages of the characteristic PB (preliminary breakdown) pulses (Figures 5b,d,e) and a single positive return stroke (Figures 5c,f).Based on Figure 5f, the peak current of the return stroke is estimated to be about 38 kA [23].Figure 6 shows its 3D mapping through DALMA.A total of 881 sources were located by DALMA.According to this mapping, the lightning flash started from a height of about 2.4 km MSL at a distance of about 15 km away from the central position of InLMA.The main sources mapped by DALMA occurred in the first 150 ms, which agrees with the waveform shown in Figure 5.An initial negative leader was observed to propagate upward and then divide into two main parts.The return stroke happened about 62.6 ms after lightning initiation.As shown in Figure 6, the lightning discharge sources spread into an area with a horizontal scale of several tens of kilometers, while only a few sources appeared above InLMA.An example waveform of the corresponding raw data of InLMA along with the spectrums of the main lightning signal and background noise, separately denoted, is shown in Figure 7.The lightning signal mainly appears in the first 150 ms with a frequency bandwidth from about 1 MHz to more than 300 MHz.Previously, we have tried various frequency bands for the mapping and found some small differences in the mapping results [24], but in this study, we only use the frequency band from 80 MHz to 230 MHz.  in Figure 7.The lightning signal mainly appears in the first 150 ms with a frequency bandwidth from about 1 MHz to more than 300 MHz.Previously, we have tried various frequency bands for the mapping and found some small differences in the mapping results [24], but in this study, we only use the frequency band from 80 MHz to 230 MHz.

Preliminary Results
When performing cross correlation to find the time difference between the different antennas of an InLMA station or different stations, there is a need to choose the proper length of the window and the sliding step.In this study, 1024 sampling points were chosen for the window and 256 sampling points were chosen for the sliding step after trial and error.When solving the χ 2 equation, the random measurement errors of elevation and azimuth σ i el , σ i az are chosen as one degree, based on our previous study using upward lightning data [24].The error of σ i t is chosen as 100 ns limited by the time accuracy of the GPS 1PPS signal.

2D Mapping Using Individual Stations
In order to demonstrate the basic functions of our system, we start by showing 2D mapping results from one station of InLMA in Figure 8 with (a) azimuth versus time, (b) elevation versus time, and (c) azimuth versus elevation for the initial 20 ms processes of the flash.The time is color-coded.In addition, Figure 8a shows the waveform of the corresponding VHF signal recorded by InLMA after being filtered into bandwidth from 80 MHz to 230 MHz.This bandwidth is used for all the mapping of InLMA in this study, as described above.The simultaneous E-field change waveform recorded by FALMA is shown in Figure 8b.As shown in Figure 8c, this flash started from an elevation of about 10 • with a negative leader moving upward.After propagating about 0.4 ms, the leader divided into two branches.Branch 1 continues to move upward and branch 2 propagates horizontally to the left.The first preliminary breakdown apparently occurred along Branch 1 and the second preliminary breakdown occurred along branch 2. The preliminary breakdown process should be an indicator of a strong electric field region.Therefore, this lightning flash initially appears to involve two regions of strong electric fields with their vectors in completely different directions.This result again indicates the complicated nature of the charge structure of winter thunderstorms.

Preliminary Results
When performing cross correlation to find the time difference between the different antennas of an InLMA station or different stations, there is a need to choose the proper length of the window and the sliding step.In this study, 1024 sampling points were chosen for the window and 256 sampling points were chosen for the sliding step after trial and error.When solving the  equation, the random measurement errors of elevation and azimuth  ,  are chosen as one degree, based on our previous study using upward lightning data [24].The error of  is chosen as 100 ns limited by the time accuracy of the GPS 1PPS signal.

2D Mapping Using Individual Stations
In order to demonstrate the basic functions of our system, we start by showing 2D mapping results from one station of InLMA in Figure 8 with (a) azimuth versus time, (b) elevation versus time, and (c) azimuth versus elevation for the initial 20 ms processes of the flash.The time is color-coded.In addition, Figure 8a shows the waveform of the corresponding VHF signal recorded by InLMA after being filtered into bandwidth from 80 MHz to 230 MHz.This bandwidth is used for all the mapping of InLMA in this study, as described above.The simultaneous E-field change waveform recorded by FALMA is shown in Figure 8b.As shown in Figure 8c, this flash started from an elevation of about 10° with a negative leader moving upward.After propagating about 0.4 ms, the leader divided into two branches.Branch 1 continues to move upward and branch 2 propagates horizontally to the left.The first preliminary breakdown apparently occurred along Branch 1 and the second preliminary breakdown occurred along branch 2. The preliminary breakdown process should be an indicator of a strong electric field region.Therefore, this lightning flash initially appears to involve two regions of strong electric fields with their vectors in completely different directions.This result again indicates the complicated nature of the charge structure of winter thunderstorms.InLMA and DALMA show consistent mapping results.InLMA mapped 1266 sources, while DALMA only mapped 119 sources.Clearly, as expected, InLMA has shown a more detailed progression of processes.The average speeds of branch 1 and branch 2 are about 1.5 × 10 6 m/s and 1.4 × 10 6 m/s, respectively.Figure 10 compares the source heights mapped by InLMA and DALMA.It appears that the DALMA-mapped sources of the first preliminary breakdown tend to have a lower height than those mapped by InLMA.Since the first PB is moving in an upward direction, this result may indicate that a PB discharge event involves many small-scale discharge events, as detected by InLMA, and a big scale event at the bottom that radiates lower frequency radio waves, as detected by DALMA.Additionally, the discharge channel mapped by DALMA is much thinner than those mapped by InLMA.It appears that a PB process may involve many small branches that have been mapped by InLMA, but not by DALMA.Interestingly, during the very initial stage, there is a small branch that has been mapped by DALMA, which is circled by a dotted line in both Figures 9b and 10, but InLMA could not map any sources from this small branch.

3D Mapping by Using 3 Stations
(b) Elevation result with E-field change versus time; (c) 2D view of azimuth versus elevation with two branches, 1 and 2, identified.

3D Mapping by Using 3 Stations
Figure 9 compares the 3D mappings of InLMA and DALMA for the initial 10 ms of the flash with the (a) InLMA mapping result and (b) DALMA mapping result.Time is color-coded.Three red pentagrams represent the three stations of InLMA.Generally, InLMA and DALMA show consistent mapping results.InLMA mapped 1266 sources, while DALMA only mapped 119 sources.Clearly, as expected, InLMA has shown a more detailed progression of processes.The average speeds of branch 1 and branch 2 are about 1.5 × 10 6 m/s and 1.4 × 10 6 m/s, respectively.Figure 10 compares the source heights mapped by InLMA and DALMA.It appears that the DALMA-mapped sources of the first preliminary breakdown tend to have a lower height than those mapped by InLMA.Since the first PB is moving in an upward direction, this result may indicate that a PB discharge event involves many small-scale discharge events, as detected by InLMA, and a big scale event at the bottom that radiates lower frequency radio waves, as detected by DALMA.Additionally, the discharge channel mapped by DALMA is much thinner than those mapped by InLMA.It appears that a PB process may involve many small branches that have been mapped by InLMA, but not by DALMA.Interestingly, during the very initial stage, there is a small branch that has been mapped by DALMA, which is circled by a dotted line in both Figures 9b and 10, but InLMA could not map any sources from this small branch.

3D Mapping by Using 3 Stations
Figure 9 compares the 3D mappings of InLMA and DALMA for the initial 10 ms of the flash with the (a) InLMA mapping result and (b) DALMA mapping result.Time is color-coded.Three red pentagrams represent the three stations of InLMA.Generally, InLMA and DALMA show consistent mapping results.InLMA mapped 1266 sources, while DALMA only mapped 119 sources.Clearly, as expected, InLMA has shown a more detailed progression of processes.The average speeds of branch 1 and branch 2 are about 1.5 × 10 6 m/s and 1.4 × 10 6 m/s, respectively.Figure 10 compares the source heights mapped by InLMA and DALMA.It appears that the DALMA-mapped sources of the first preliminary breakdown tend to have a lower height than those mapped by InLMA.Since the first PB is moving in an upward direction, this result may indicate that a PB discharge event involves many small-scale discharge events, as detected by InLMA, and a big scale event at the bottom that radiates lower frequency radio waves, as detected by DALMA.Additionally, the discharge channel mapped by DALMA is much thinner than those mapped by InLMA.It appears that a PB process may involve many small branches that have been mapped by InLMA, but not by DALMA.Interestingly, during the very initial stage, there is a small branch that has been mapped by DALMA, which is circled by a dotted line in both Figures 9b and 10, but InLMA could not map any sources from this small branch.

An Attempt to Better Characterize Preliminary Breakdowns
Preliminary breakdown processes have been studied by many authors [25][26][27][28].However, concerning each individual PB pulse, little is known about how and what happens.In this section, by taking advantage of InLMA, we try to reveal more detailed features for individual PB pulse processes.
Figure 11 shows the mapped height along with the three stations' azimuths and elevations of PB pulse discharges during the initial 1 ms with (a) for height, (b) for azimuths, and (c) for elevations.For comparison, the corresponding E-field change is also included in Figure 11a.As evident in Figure 11b, the located points by station 1 are apparently much fewer than those by station 2 and station 3 due to the larger background noise at station 1.Since 3D located points are basically limited by the station with the fewest points, the height points in Figure 11a are also much fewer than the azimuth and elevation points of station 2 and station 3. Additionally, as seen in Figure 11, for all the individual PB discharges, particularly those numbered from one to seven with simple structures, there is consistently no considerable change in the corresponding azimuths, but the elevations vary significantly.The height variation in Figure 11a can be approximately treated as a 3D distance change.It appears that each PB pulse process has a spatial scale of around 150 m along its progression direction.Furthermore, we note that the sub-discharges of two adjacent pulse discharges for the pulses numbered from one to seven in Figure 11a consistently appear progressively in connection to each other, though they are discontinuous in time.The progression speed measured from pulse one to pulse seven is about 1.4 × 10 6 m/s.
Preliminary breakdown processes have been studied by many authors [25][26][27][28].How ever, concerning each individual PB pulse, little is known about how and what happen In this section, by taking advantage of InLMA, we try to reveal more detailed features f individual PB pulse processes.
Figure 11 shows the mapped height along with the three stations' azimuths and e vations of PB pulse discharges during the initial 1 ms with (a) for height, (b) for azimuth and (c) for elevations.For comparison, the corresponding E-field change is also includ in Figure 11a.As evident in Figure 11b, the located points by station 1 are apparently mu fewer than those by station 2 and station 3 due to the larger background noise at statio 1.Since 3D located points are basically limited by the station with the fewest points, t height points in Figure 11a are also much fewer than the azimuth and elevation points station 2 and station 3. Additionally, as seen in Figure 11, for all the individual PB d charges, particularly those numbered from one to seven with simple structures, there consistently no considerable change in the corresponding azimuths, but the elevatio vary significantly.The height variation in Figure 11a can be approximately treated as a 3 distance change.It appears that each PB pulse process has a spatial scale of around 150 along its progression direction.Furthermore, we note that the sub-discharges of two a jacent pulse discharges for the pulses numbered from one to seven in Figure 11a consi ently appear progressively in connection to each other, though they are discontinuous time.The progression speed measured from pulse one to pulse seven is about 1.4 × 1 m/s.ground by optical means [35,36].Since return strokes tend to have a decreasing speed when they go up [37], their average speeds along their entire channels should be smaller.Indeed, for one negative return stroke, Shao et al. reported an average speed of 3-5 × 10 7 m/s [33].We believe that the average speed that we measured is reasonable.

Discussion and Concluding Remarks
To the best of our knowledge, the lightning mapping system reported in this paper is the first lightning mapping array that is formed by three broadband interferometers.The method proposed in this paper for conducting 3D mapping with multiple interferometers is also new.Although the system demonstrated in this paper only consists of three stations, it can be readily expanded to have more stations.More stations can not only cover a wider area, but also reduce uncertainty in the mapping, similar to the lightning mapping array using TOA techniques.In order to observe winter lightning in Japan, whose discharge sources usually occur with lower heights, we intentionally limited the distances between stations.In the cases of observing summer lightning, whose discharge sources usually occur at higher heights, we think the station distances can be increased to more than 10 km.In such cases, the cross correlation between stations may become difficult.As a matter of fact, in the method we propose, even without the cross correlation between stations, 3D mapping can still be achieved only by using the azimuth and elevation information of each station.
Despite the lightning analyzed in this study being more than 10 km away from our InLMA stations, we can still perform some detailed mappings, particularly for the PB processes, which exhibited two stages.Interestingly, these two stages of PB processes occurred, respectively, along two branches with completely different orientations.The corresponding charge structure must be very complicated.Each individual PB pulse involves many small-scale discharges that usually scatter in a vertical scale of around 150 m.The smaller scale (1-4 m scale) discharges appear continuously in time but scatter in a wider range of height, while discharges with a scale of about 10 m tend to occur in a discrete manner.For some small-scale discharges, we can identify their progressions.The PB pulse itself should be a discharge with a scale of some hundred meters, which tends to be located at a lower height than its small-scale discharges.Our results on PB pulse discharges appear

Figure 1 .
Figure 1.Schematic of the interferometer hardware at one station.Figure 1.Schematic of the interferometer hardware at one station.

Figure 1 .
Figure 1.Schematic of the interferometer hardware at one station.Figure 1.Schematic of the interferometer hardware at one station.

Figure 2 .
Figure 2. Antenna site map of InLMA, DALMA, and FALMA with the filling color representing the height above sea level.

Figure 2 .
Figure 2. Antenna site map of InLMA, DALMA, and FALMA with the filling color representing the height above sea level.

Figure 3 .
Figure 3. Illustration of three station InLMA with an assumed point discharge source (x, y, z) which radiates electromagnetic waves to the three stations.The corresponding angles are denoted and used to find the 3D location of the source (x, y, z).

Figure 4 .
Figure 4. (a) Flow chart used in 3D mapping program; (b) Description of the method for finding the time arrival differences windows between two stations.Figure 4. (a) Flow chart used in 3D mapping program; (b) Description of the method for finding the time arrival differences windows between two stations.

Figure 4 .
Figure 4. (a) Flow chart used in 3D mapping program; (b) Description of the method for finding the time arrival differences windows between two stations.Figure 4. (a) Flow chart used in 3D mapping program; (b) Description of the method for finding the time arrival differences windows between two stations.

Figure 5 .
Figure 5. E-field change recorded at one station of FALMA.(a) 1 s of waveform; (b) waveform from 44 ms to 58 ms including two preliminary breakdowns; (c) waveform of the return stroke; (d-f) Expanded waveforms to indicate the 1st preliminary breakdown stages and the 2nd preliminary breakdown and return stroke in detail.The "atmospheric electricity" sign convention is used.

Figure 5 .
Figure 5. E-field change recorded at one station of FALMA.(a) 1 s of waveform; (b) waveform from 44 ms to 58 ms including two preliminary breakdowns; (c) waveform of the return stroke; (d-f) Expanded waveforms to indicate the 1st preliminary breakdown stages and the 2nd preliminary breakdown and return stroke in detail.The "atmospheric electricity" sign convention is used.

Figure 6 .
Figure 6.3D mapping result of the lightning flash by DALMA color-coded in time.Three pentagrams represent stations of InLMA.(a) 3D mapping; (b) Top view.Figure 6. 3D mapping result of the lightning flash by DALMA color-coded in time.Three pentagrams represent stations of InLMA.(a) 3D mapping; (b) Top view.

Figure 6 .
Figure 6.3D mapping result of the lightning flash by DALMA color-coded in time.Three pentagrams represent stations of InLMA.(a) 3D mapping; (b) Top view.Figure 6. 3D mapping result of the lightning flash by DALMA color-coded in time.Three pentagrams represent stations of InLMA.(a) 3D mapping; (b) Top view.

Figure 6 .
Figure 6.3D mapping result of the lightning flash by DALMA color-coded in time.Three pentagrams represent stations of InLMA.(a) 3D mapping; (b) Top view.

Figure 7 .
Figure 7. Example waveform and frequency spectrum of the raw lightning data recorded at one channel of InLMA.(a) Signal waveform and background noise; (b,d,f) Frequency spectrum of the flash signal; (c,e,g) Frequency spectrum of the noise.

Figure 7 .
Figure 7. Example waveform and frequency spectrum of the raw lightning data recorded at one channel of InLMA.(a) Signal waveform and background noise; (b,d,f) Frequency spectrum of the flash signal; (c,e,g) Frequency spectrum of the noise.

Figure 8 .
Figure 8. 2D result from station 3 of InLMA.(a) Azimuth result with InLMA waveform versus time; (b) Elevation result with E-field change versus time; (c) 2D view of azimuth versus elevation with two branches, 1 and 2, identified.

Figure 9
Figure 9 compares the 3D mappings of InLMA and DALMA for the initial 10 ms of the flash with the (a) InLMA mapping result and (b) DALMA mapping result.Time is color-coded.Three red pentagrams represent the three stations of InLMA.Generally,

Figure 9 .
Figure 9.A comparison between 3D mapping results of InLMA and DALMA for the initial 10 ms of the flash.(a) 3D mapping of InLMA with locations of three stations of InLMA.(b) 3D location result from DALMA.

Figure 9 .
Figure 9.A comparison between 3D mapping results of InLMA and DALMA for the initial 10 ms of the flash.(a) 3D mapping of InLMA with locations of three stations of InLMA.(b) 3D location result from DALMA.

Figure 9 .
Figure 9.A comparison between 3D mapping results of InLMA and DALMA for the initial 10 ms of the flash.(a) 3D mapping of InLMA with locations of three stations of InLMA.(b) 3D location result from DALMA.

Figure 10 .
Figure 10.A comparison between the 3D heights mapped by InLMA and DALMA along with the corresponding E-field change.

Figure 11 .
Figure 11.Mapping results for the PB pulses occurred during the initial 1 ms stage.(a) 3D mapping by InLMA along with the E-field change.The pulses numbered from one to seven have a simple structure and will be studied further next; (b) Azimuth versus time for all three InLMA stations; (c) Elevation versus time for all three InLMA stations.For these seven pulses, we have tried a more detailed analysis by comparing the waveforms of the E-field change, the raw InLMA data, the 3 MHz~30 MHz, 30 MHz~80 MHz, and 80 MHz~230 MHz bandpass data, and the corresponding azimuth and elevation mapped by station 3, as shown in Figure12.Figure12a-dhave the same format, except for different pulses, with the top panel for E-field change, the second to the fourth panels for the 3 MHz-30 MHz, 30 MHz-80 MHz, and 80 MHz-230 MHz data, and the fifth panel for the azimuth and elevation.

Figure 13 .
Figure 13.2D mapping result of the discharge processes around the positive return stroke.(a) Azimuth and InLMA raw waveform; (b) Elevation and FALMA E-field change (c) 2D vertical view of elevation versus azimuth.Black points are all discharge sources occurring during the period from lightning initiation until 62.2 ms later; color points coded in time are the discharge sources during the period around the positive return stroke.Two clusters are circled with red dotted lines, with cluster 1 occurring before the return stroke, likely due to the downward positive leader of the stroke, while cluster 2 occurs at the initial stage of the return stroke.

Figure 13 .
Figure 13.2D mapping result of the discharge processes around the positive return stroke.(a) Azimuth and InLMA raw waveform; (b) Elevation and FALMA E-field change (c) 2D vertical view of elevation versus azimuth.Black points are all discharge sources occurring during the period from lightning initiation until 62.2 ms later; color points coded in time are the discharge sources during the period around the positive return stroke.Two clusters are circled with red dotted lines, with cluster 1 occurring before the return stroke, likely due to the downward positive leader of the stroke, while cluster 2 occurs at the initial stage of the return stroke.