A meteor head echo analysis algorithm for the lower VHF band

We have developed an automated analysis scheme for meteor head echo observations by the 46.5 MHz Middle and Upper atmosphere (MU) radar near Shigaraki, Japan (34.85 N, 136.10 E). The analysis procedure computes meteoroid range, velocity and deceleration as functions of time with unprecedented accuracy and precision. This is crucial for estimations of meteoroid mass and orbital parameters as well as investigations of the meteoroid-atmosphere interaction processes. In this paper we present this analysis procedure in detail. The algorithms use a combination of singlepulse-Doppler, time-of-flight and pulse-to-pulse phase correlation measurements to determine the radial velocity to within a few tens of metres per second with 3.12 ms time resolution. Equivalently, the precision improvement is at least a factor of 20 compared to previous single-pulse measurements. Such a precision reveals that the deceleration increases significantly during the intense part of a meteoroid’s ablation process in the atmosphere. From each received pulse, the target range is determined to within a few tens of meters, or the order of a few hundredths of the 900 m long range gates. This is achieved by transmitting a 13-bit Barker code oversampled by a factor of two at reception and using a novel range interpolation technique. The meteoroid velocity vector is determined from the estimated radial velocity by carefully taking the location of the meteor target and the angle from its trajectory to the radar beam into account. The latter is determined from target range and bore axis offset. We have identified and solved the signal processing issue giving rise to the peculiar signature in signal to noise ratio plots reported byGalindo et al.(2011), and show how to use the range interpolation technique to differentiate the effect of signal processing from physical processes.


Introduction
The flux of meteoroids onto Earth is the source of the neutral and ion metal layers in the middle atmosphere.The influx plays an important role in atmospheric dynamics and processes like the formation of high-altitude clouds, possibly through coagulation of meteoric smoke particles acting as condensation nuclei for water vapor (Summers and Siskind, 1999;Megner et al., 2006).Hunten et al. (1980) point out that estimating the deposition of mass in the atmosphere requires knowledge of not only the total mass influx of meteoroids, but also the size and velocity distributions and physical characteristics such as density and boiling point of the particles.
Meteor head echo observations with High-Power Large-Aperture (HPLA) radars are well suited for studying many aspects of the meteoroid influx in detail, as well as the atmosphere interaction processes (e.g.Pellinen-Wannberg, 2005).Meteor head echoes are radio waves scattered from the intense regions of plasma surrounding and co-moving with meteoroids during atmospheric flight.Head echoes were first reported by Hey et al. (1947), who observed the Giacobinid (now called Draconid) meteor storm of 1946 with a 150 kW VHF radar.HPLA radar systems, however, have a peak transmitter power of the order of 1 MW and array or dish antenna apertures in the range of about 800-7 × 10 4 m 2 (Pellinen- Wannberg, 2001), focusing their antenna gain pattern into a narrow main beam with a full-widthat-half-maximum (FWHM) of the order of 1 • at the VHF and/or UHF operating frequencies.This high power density permits numerous head echo detections from faint meteors.
Since the 1990s, head echo observations have been conducted with most HPLA radar facilities around the world (Pellinen- Wannberg and Wannberg, 1994;Mathews et al., 1997;Close et al., 2000;Sato et al., 2000;Chau and Woodman, 2004;Mathews et al., 2008;Malhotra and Mathews, 2011).These radar systems have diverse system characteristics in terms of operating frequency, dish or phased array antenna, aperture size etc.Characteristics of all but the Resolute Bay Incoherent Scatter Radar (RISR) are summarized in Table 1 of Janches et al. (2008).Methods of head echo analysis have been developed more or less independently at several of the facilities, and with emphasis on different aspects of meteor science and/or radio science issues.In Sect. 2 we give a brief review of references to previously developed meteor head echo analysis methods, and point out new features in our approach.
We have developed and implemented an automated analysis scheme for meteor head echo observations by the 46.5 MHz Middle and Upper atmosphere (MU) radar near Shigaraki, Japan (34.85 • N, 136.10 • E) (Fukao et al., 1985).Previous meteor head echo observations with the MU radar have been reported by Sato et al. (2000) and Nishimura et al. (2001).
The algorithms presented here use a combination of single-pulse-Doppler, time-of-flight and pulse-to-pulse phase correlation measurements, enabling the meteoroid radial velocity to be determined to within a few tens of m s −1 with 3.12 ms time resolution.Equivalently, the precision improvement of the determined line-of-sight velocity is at least a factor of 20 compared to previous single-pulse measurements.Furthermore, we have invented an interpolation scheme to find the target range within a small fraction of the 900 m long range gates.This, together with the upgrade of the MU radar receiver system from four analog to 25 digital channels (Hassenpflug et al., 2008), results in improved target position determination, crucial for accurately estimating meteoroid trajectory parameters and calculating true meteoroid velocities from the measured radial velocity component along the radar line-of-sight.
A block diagram of the analysis scheme is shown in Fig. 1.The paper is organized to describe the blocks as follows.
A brief description of the MU radar and our experimental settings is found in Sect.3. The initial search for meteor events (block A) and an overview of the decoding procedure (block B) is given in Sect. 4. The new range interpolation technique is presented in Sect. 5.It solves a major and systematic signal processing issue in meteor head echo data (Galindo et al., 2011), further discussed in Sect.5.1.The selection of data points constituting one meteor event (block C) is described in Sect.6.These data points are then subject to the pulse-to-pulse phase correlation technique (block D) reported in Sect.7.
The instantaneous position of a meteor target at each interpulse period (IPP) is determined by interferometry (block E) using the MUSIC method (Schmidt, 1986) explained in Sect.8. Meteoroid trajectories and radiant error estimations are determined by combining the interferometry data with the estimated range data and the radial velocity (block F) as detailed in Sect.9.In turn, the trajectories are used to redo the parts of the analysis given in blocks B-D, but with the re-ceiver beam post-steered towards the most probable location of the meteor target at each IPP.
In Sect. 10 we present meteor target radar cross section (RCS) calculations.Comparing the RCSs with and without post-steering the receiver beam gives an estimate of the validity of the position determination and the applied antenna gain pattern.To ensure as accurate position determination as possible, we have adopted an interchannel calibration routine, described in Sect.11.
The analysis output parameters are range, altitude, radial velocity, meteoroid velocity, instantaneous target position, RCS and meteor radiant.The parameter values calculated both with and without using post-beam steering are stored in a data base.
2 Meteor analysis methods at other HPLA radars Evans (1965Evans ( , 1966) ) describes the first head echo measurements with what today is termed a HPLA radar.He used the 440 MHz Millstone Hill radar, which has an operating frequency about an order of magnitude higher than classical specular meteor trail radar systems.Evans maximized the cross-beam detection area of the Geminid, Quadrantid and Perseid meteor showers by pointing the Millstone Hill radar towards the shower radiants at times when the radiants were located at very low elevations above the local horizon.This enabled velocity and deceleration determination for meteors belonging to the showers, for which the atmospheric trajectories were aligned with the radar beam.

EISCAT
Pellinen- Wannberg and Wannberg (1994) presented the first of the modern time HPLA meteor head echo observations.These were conducted using the radar systems of the European Incoherent SCATer (EISCAT) Scientific Association.Details of the analysis methods, focusing on how to identify, extract and analyse highly Doppler shifted meteor events in conventional Barker-coded power profile type incoherent scatter measurements, were reported by Wannberg et al. (1996).The earliest EISCAT observations were limited to time integrated data and therefore had time resolution of 2 s.
The first sets of tristatic meteor observations, using all three receiver stations of the EISCAT UHF radar, were conducted by Pellinen- Wannberg et al. (1999) and Janches et al. (2002).In contrast to monostatic observations, which only give the radial (line-of-sight) velocity component, multistatic (and also interferometric) observations enable calculation of the meteoroid velocity vector.However, the initially used antenna beam pointing geometry was such that the velocity vector suffered from large uncertainties (Wannberg et al., 2008).An improved geometry, where the linear dependence of the measured velocity components is minimised was therefore developed.Wannberg et al. (2008)  detailed description of this improvement and the signal processing development following the installation of the new digital signal processing and raw data recording systems in 2001, which enabled phase-coherent pulse-by-pulse analysis.Kero et al. (2008a) present a method for finding the position of a compact meteor target in the common volume monitored by the three UHF receivers, and how velocity, deceleration, RCS and meteoroid mass were estimated from the improved tristatic observations.The EISCAT UHF radar provided excellent precision and accuracy of meteors observed with all three widely separated receiver systems, but low rate of such events ( 10 h −1 ), mainly due to the small tristatic measurement volume (Szasz et al., 2008).Zhou et al. (1995) observed the first head echoes using the Arecibo Observatory (AO) 430 MHz UHF radar.The observations were limited to time integrated data and a time resolution of 11 s.Mathews et al. (1997) followed up the AO observations with an improved non-integrated data collection approach, enabling 1 ms time resolution.The Doppler technique for obtaining the instantaneous meteor Doppler velocity and deceleration is described in Janches et al. (2000aJanches et al. ( ,b, 2001)).Subsequent improvement of the signal processing techniques at AO has been particularized by Mathews et al. (2003); Wen et al. (2004Wen et al. ( , 2005a)); Briczinski et al. (2006); Wen et al. (2005bWen et al. ( , 2007)).The emphases of the analysis technique development have been to implement automated real-time analysis of meteor parameters (Wen et al., 2004), remove non-periodic bursty interference (Wen et al., 2005b), and separate incoherent scatter from meteor signals (Wen et al., 2005a(Wen et al., , 2007)).(Close et al., 2002).No conclusive evidence of shower meteor detections were found, in accordance with the (subsequently estimated) very low probability of detecting such meteors during the observations (Brown et al., 2001).Close et al. (2005) present a method for meteoroid mass estimation by converting the measured RCS to head echo plasma density utilizing a spherical electromagnetic scattering model.The ALTAIR radar has multi-frequency capability, and can transmit linear frequency modulated chirped pulses.This enables a variety of meteoroid range rate calculations, e.g. based on the difference in the measured ranges due to range-Doppler coupling (Loveland et al., 2011).Mathews et al. (2008) applied the analysis methods developed for the 430 MHz AO radar and described by Mathews et al. (2003) and Briczinski et al. (2006), to the 449.3 MHz 32 panel Advanced Modular Incoherent Scatter Radar at Poker Flat Alaska (PFISR-32), to the 1290 MHz Sondrestrom Radar Facility (SRF), and later also to the Resolute Bay Incoherent Scatted Radar (RISR) (Malhotra and Mathews, 2011).Mathews et al. (2008) estimated that AO is 77 times more sensitive than SRF and 2100 times more sensitive than PFISR.Yet, they found the lowest event rate at SRF (34 per hour) relative to PFISR (55 per hour) and AO (1000 per hour).Furthermore, the altitude distribution of SRF meteors was 10 km below that observed with AO/PFISR.These observations agree with a frequency dependent meteor head echo target RCS, further discussed in Sect.10, as well as the cut-off in the high-altitude end of the 930 MHz EISCAT UHF distribution as compared to the 224 MHz EISCAT VHF distribution (Westman et al., 2004).Sparks et al. (2009) report the results of concurrent PFISR observations using an independent but similar data analysis method.(Sparks et al., 2010) operated PFISR as a threechannel interferometer.They demonstrate that meteor radiants and orbits can be determined.Chau et al. (2009) describe an antenna compression approach to widen the PFISR beam width for meteor head echo observations, to about three times the width of the ordinary narrow beam.Chau et al. corrected the signal-to-noise ratio (SNR) depending on where in the beam meteors were detected, thus estimated a corrected relative RCS distribution, i.e. as if all meteors were detected within the narrow main lobe.Using a wider beam to detect a larger number of strong and/or long-duration meteor head echo events, which would not have been detected in the narrow beam, is an interesting and promising approach.However, it is not a necessary pro-cedure to enable beam shape correction for interferometric observations with the MU radar.Chau and Woodman (2004) and Chau et al. (2007) used the 50 MHz Jicamarca Radio Observatory (JRO) radar for meteor head echo observations.They utilize three-channel interferometry to calculate meteoroid trajectories and convert the radial velocity to vector velocity.13-bit Barker coded pulse sequences were transmitted to decrease interference from geophysical clutter, and pulse-to-pulse phase correlation was used to estimate radial deceleration.The sampling rate was equal to the subpulse (baud) rate (Chau et al., 2007, Table 1).Chau and Galindo (2008) report the first interferometric head echo observations of meteor shower particles.Galindo et al. (2011) describe a signal processing issue in JRO data that manifests itself as a peculiar signature in SNR plots.

Discussion
The outline of the analysis technique presented in this paper is similar to that presented by Chau and Woodman (2004) for interferometric JRO observations.The way target range and Doppler velocity are extracted from the raw data in a multi-step matched-filter procedure (Sects.3-4) largely follows the EISCAT analysis technique detailed by Wannberg et al. (2008) and Kero et al. (2008a).
The first main difference between the method at hand and published methods is that we have developed a range finding interpolation technique for BPSK (binary phase-shift key) coded pulse sequences (Sect.5).This technique solves the systematic signal processing issue causing ripples in the JRO data reported by Galindo et al. (2011).Also, Chau and Woodman (2004) report that there is a bias between the JRO time-of-flight velocity estimation and Doppler estimation.Our determined radial velocity component (Sect.3, Fig. 5) is unbiased, similarly to EISCAT observations (Wannberg et al., 2008).
Furthermore, we have developed an algorithm that combines Doppler velocity and velocity determined from pulseto-pulse phase correlation in such a way that we always find the most probable phase correlation velocity from a set of ambiguous possibilities (Sect.7).In retrospect, we have found that our approach is very similar to that developed by Elford (1999) and used to analyze occasional strong head echoes observed with the 30 kW Buckland Park 54.1 MHz narrow-beam VHF radar at the University of Adelaide (e.g.Cervera et al., 1997, and references therein).Our method is tailored to work as an automatic procedure on typical HPLA radar data containing numerous weak head echoes and using a transmitted pulse length much longer than the samples in the receiver data stream, but the general idea is the same.
The third main difference is that we have implemented interferometry utilizing all 25 channels of the MU radar receiver system (Sect.8), enabling unambiguous target localization.Three receiver channels were used for interferometric JRO observations (Chau and Woodman, 2004), as well as previous interferometric MU observations (Nishimura et al., 2001).Three channels are, in principle, enough to locate meteors inside the transmitter beam, but Chau et al. (2009) note that more than three antennas are required to remove angular ambiguities as a significant fraction of the meteors appear in sidelobes.
The fourth principal difference is the way we convert radial velocity to vector velocity (Sect.9).Equation (4) in Chau and Woodman (2004) is a good approximation, but does not utilize all information of a meteor event.
Improving the accuracy and precision of meteoroid velocity vector determination in head echo observations is important to provide useful data for the modelling of Solar System dust, e.g. for studying the evolution of meteoroid streams and predicting meteor shower outbursts (Jenniskens, 2006;Sato and Watanabe, 2007;Atreya et al., 2010).

The MU radar experimental setup
The present setup of the MU radar hardware comprises a 25 channel digital receiver system.It was upgraded from the original setup (Fukao et al., 1985) in 2004 and is described by Hassenpflug et al. (2008).After the upgrade, the MU radar always transmit right-handed circular (RC) polarization and receive left-handed circular (LC) polarization, with a phase accuracy of 2 • .The output of each digital channel is the sum of the received radio signal from a subgroup of 19 Yagi antennas.The whole array consists of 475 antennas, evenly distributed in a 103 m circular aperture.A schematic view of the array and the subgroups is given in Fig. 2. It is possible to combine the output from several subgroups into the same digital channel to reduce the total number of channels and hence decrease the data rate without decreasing the total aperture.We have, however, chosen to use all 25 channels to enable subgroup phase offset reduction and to optimize interferometric target position determination and post-steering of the receiver beam.The maximum continuous data rate is about 20 GB h −1 due to system limitations.
The range finding interpolation algorithm works best if the transmitted code is selected as to have a minimum value next to the central maximum in its autocorrelation function (ACF).The autocorrelation of a 13-bit Barker code has zeros next to the central peak.This property maximizes the precision of the range interpolation for a given code length.Other properties of the transmission schedule as the number of bits in the code, baud length, IPP, etc., are not restricted by the range finding interpolation algorithm and should be chosen according to hardware limitations and other constraints.(Fukao et al., 1985).The array is divided into 25 subgroups (A1-F5), each consisting of 19 antennas and connected to its own transmitter and receiver module (Hassenpflug et al., 2008).
However, the implementation described in Sect. 5 is designed for a radar setup where the transmitted pulse sequence is oversampled at reception.In the present MU observations it was oversampled by a factor of two, accomplished by using a sampling period of T s = 6 µs while transmitting the Barker code with a 12 µs baud length.The MU radar hardware does not allow receiver sampling period and transmitter subpulse length to differ.Each 12 µs baud of the 13-bit code is therefore defined as two 6 µs subpulses of equal phase in the radar experimental setup definition file.The transmitter and receiver bandwidths are defined by the 6 µs subpulse length and the 6 µs sampling period, and approximately equal to b w = 1/6 µs 167 kHz.In the decoding procedure we use an ideal, boxcar version of the transmitted code pattern, as exemplified in Fig. 3, and further described in Sect. 5 and by Eq. (3).
It is important to make sure that the receiver bandwidth is wide enough to accommodate both the modulation bandwidth and the target Doppler shift.Figure 4 shows the spectral width of the transmitted code, the receiver bandwidth, the received spectrum of a meteor with zero radial velocity, and the spectrum of a meteoroid with 70 km s −1 radial velocity, corresponding to a Doppler shift of 21.7 kHz.The energy loss from the most Doppler-shifted meteors due to finite receiver bandwidth, compared to non-Doppler shifted meteors, is less than 5 %.This loss is small enough to be negligible in the estimated target RCS that due to other reasons vary over several orders of magnitudes (Sect.10).However, the slight loss of received energy is asymmetric, as can be seen in Fig. 4. To confirm the validity of our Doppler estimates, Fig. 5 displays a comparison to independent time-of-flight estimates.They agree to within the order of one part in a thousand.This comparison demonstrates that the small but asymmetric loss of spectral energy does not bias the Doppler estimation of the velocity.Furthermore, the result agrees with the investigation by Wannberg et al. (2008), who found that no contribution from slipping plasma could be detected within the measurement accuracy of the EISCAT UHF meteor observations, and that the Doppler velocities were unbiased.Doppler and time-of-flight methods are further described in Sect. 5.The selection of a 156 µs pulse length and 6 µs sampling is a tradeoff between time resolution, range resolution, signalto-noise ratio (SNR) and the maximum possible data rate at the present MU radar system.A longer pulse length would indeed improve the SNR but also increase the time between consecutive pulses due to the 5 % transmitter duty cycle limitation.
A longer pause between pulses has two drawbacks; it increases the ambiguity of velocity data calculated from pulseto-pulse phase correlations (cf.Sect.7) and decreases the time resolution of the determined meteoroid parameters for one and the same meteor event.We have tried several different setups in search for a good tradeoff and found that increasing the IPP beyond ∼3 ms complicates the selection procedure of a velocity for the meteoroid among ambiguous possibilities determined by pulse-to-pulse phase correlation, described further in Sect.7. The IPP we finally decided for, T ipp = 3.12 ms, gives a separation v equal to between possible ambiguous velocities.Due to data rate limitations a 13-bit code is the longest code we can oversample by a factor of two at reception using all 25 channels of the MU radar and still monitor the most important part of the meteor zone, an altitude interval ∼70-130 km, where most meteor head echoes appear.Using a 156 µs long pulse together with a 6 µs sampling period gives a range interval of 73-127 km from where the echo of the whole transmitted pulse sequence is received.

Initial analysis: range and Doppler
The meteor head echo analysis procedure starts with a simple scan of the data.It is similar to the meteor data analysis performed on EISCAT VHF and UHF radar data described by Wannberg et al. (2008) and Kero et al. (2008a).The scanning procedure performs a search in the power domain by computing the boxcar function of 26 consecutive range gates (the length of a point target echo) and compares the result to the noise.If the boxcar function of seven consecutive IPPs exceed three noise standard deviations, the IPPs are flagged as a possible event.This choice of threshold keeps the number of false events at a reasonably low level without excluding analysable meteor head echoes.
An estimate of a meteoroid's radial (line-of-sight) velocity v r can in principle be deduced using the Doppler shift f D of one single received radar pulse as the Doppler shift depends on meteoroid velocity and operating frequency f 0 according to where c 0 is the speed of light.When a meteoric particle enters the atmosphere it will heat up in collisions with atmospheric constituents and generate a dense ionized plasma, generally detectable by radar along several kilometer of its trajectory, before the particle vanishes.To determine whether an enhanced signal in the data is due to a meteor target or not, we require the target timeof-flight velocity to agree with its Doppler shift.For this criterion to be applied, several IPPs worth of data needs to be recorded from each meteor.This is easily achievable with an IPP of 3.12 ms but demands precise range data.One range gate is defined by the sampling period, which in our observations equals an extent in range of about R s = T s c 0 /2 900 m.Our range finding interpolation algorithm enables a time-of-flight velocity calculation even when the meteor target is within one and the same range gate for the duration of the event.
Whenever there are several possible events with gaps smaller than 20 IPPs (62 ms) between them, we try treating them as one single meteor and analyze the whole set of IPPs together.Shorter gaps are, in case of MU meteor head echoes, most often caused by one meteor target of low SNR with an irregular ionization/RCS profile or moving through a minimum in the antenna radiation pattern.
Instead of limiting the temporal extent of a meteor event by a threshold on the signal and use all received radar pulses in between for determining meteor properties we have developed an automatic routine that looks for consistency in both velocity versus time and range versus time.Data points that do not fulfil the criteria are excluded.By looking for consistency we also try to include data points from before and after the initially flagged sequence of IPPs.Therefore, the interval that we analyze and look for consistency within include 20 IPPs before and after the marked sequence.
In the second step of the analysis procedure the 85 range gates from one transmission/reception is cross-correlated with a set of differently Doppler-shifted versions of the transmitted code in order to find an approximate Doppler shift and range of the echo.The unshifted code can be described as where the zero elements in both ends represent start and stop of transmission.These zero elements are necessary to define for the code sequence to be used successfully in the range finding interpolation technique described in Sect. 5.The ideal version of the code illustrated in Fig. 3 is Dopplershifted by the multiplication where the Doppler frequency f n = −30 000, −29 000, ..., 5000 Hz and T s = 6 µs.
The first guess of Doppler and range is done by selecting the Doppler frequency giving the highest peak decoded power and picking out the 26 + 1 range gates that correspond to the location of the highest peak in the cross-correlation for further analysis.

Range finding interpolation
We have invented a technique which uses the degree of asymmetry of the decoded signal to interpolate the code used in the decoding procedure and find the target range to within a few hundredths of a range gate.The twofold oversampling of the transmitted BPSK coded sequence means that the two range gates next to the peak in the cross-correlation sequence will in ideal cases have half the value of the peak, as was illustrated in Fig. 3.
The relation between sampling with a sampling period of T s and target range r is where r 0 is the target range at start of sampling (in the described observations r 0 ≈ 72 km), g r is the integer number of range gates (each of length T s c o /2 900 m) and is the remaining fraction of a range gate.The decoded signal s(r) will be symmetric with respect to a particular range gate (r g ), if and only if the target is located at a distance corresponding to an integer number of sampling periods from the radar, thus = 0.If the target location, however, is such that = 0 the signal s(r) becomes asymmetric.The ratio of the values next to the peak of the decoded power can be used to estimate according to where a and b are the differences between the peak and the adjacent range gates (illustrated in Fig. 7 described below).We use the value 1 from the first decoding attempt to select the 26 + 1 range gates containing the echo.Then we start an iterative procedure in which at each step a new value n is calculated at the same time as the Doppler shift used for decoding the signal is optimized.The optimization is accomplished by first increasing the Doppler shift with a given step size, and evaluate the cross-correlation until the decoded power is smaller than the previous value.At this point the step size is decreased and the search direction is reversed.The procedure is iterated until the step size is 5 Hz, corresponding to about 20 m s −1 .Both the Doppler frequency and the code interpolation used when calculating the cross-correlation affect the degree of symmetry and the value of the peak decoded power.Optimization of the two quantities are therefore searched for simultaneously.Each time a step is taken in Doppler frequency and a new cross correlation is computed, the new ratio n is added to form a sum of all evaluations according to = m n=1 n . (7) Generally, the absolute value | n | in each new iteration is smaller than the previous value.The sum thus forms a converging series, where m is the number of iterations.The sign of each n depends on which of a and b that are greater in that particular iteration.
The interpolated codes are found by adding zeroes to each end of the original 26 bit long code, as shown in Fig. 3, and thereafter interpolate adjacent bits of code.= 0 means that bit 1-26 of the code in Fig. 3 are used.If < 0, interpolation is performed towards left (bit 0) and if > 0 towards right (bit 27).It should be noted that = ±0.5 gives rise to zeros when the interpolation is performed on consecutive values of +1 and −1 (or −1 and +1).This is obvious also when looking at undecoded meteor data, as reception of an equal amount of signals with opposite phase cancels out to zero. Figure 6 shows an undecoded range-time and signal intensity plot of a meteor head echo detected 28 July 2009, 05:33:09 JST, and hereafter referred to as "meteor 1".Range gate 1 at the bottom of Fig. 6 is where the leading edge of an echo from a target at 73 km range appears in the data stream, while gate 60 corresponds to the leading edge of an echo from a target at 127 km range, as was described in Eq. ( 5).The occasions when the meteor target is located close to the middle of the range gates ( ±0.5) have weak signal power and are visible as dark bands in the plot.The reason for the weakened signal is that bauds transmitted with opposite phase are received in a subset of the range gates.When using the BPSK code given in Eq. ( 3), this subset consists of gates 11, 15, 19, 21, 23 and 25, given that the range gate where the leading edge of the echo appears is called number 1.It should be stressed that if = 0, the echo is spread out into 27 gates when a BPSK code consisting of 26 transmitted bauds is used.When = ±0.5, we have found that the signal in gates 11, 15, 19, 21, 23 and 25 drops completely below the noise floor, while the signal in gates 1 and 27 has a power level half that of the remaining gates 2-10, 12-14, 16-18, 20, 22, 24 and 26.The reason for this is that gate 1 and 27 contains reception of one half of a transmitted baud each.The efficient cancelling in gates 11, 15, 19, 21, 23 and 25 indicates that the meteor target is small when compared to a range gate (900 m).
The top row of Fig. 7 shows three examples of what interpolated codes look like.The columns correspond to IPP 72 (left column), 74 (middle column) and 76 (right column) of meteor 1 and are interpolated using = 0.287, = −0.477and = −0.248,respectively.The panels of the middle row illustrates the cross-correlation with a Doppler shifted but not interpolated version of the transmitted code, which clearly give asymmetric s(r).The bottom panels show the crosscorrelation with Doppler shifted interpolated codes.The interpolated codes give symmetric s(r).
To compensate for the signal power loss when = 0, the decoded signal must always be divided by the amplification of the interpolated code, which differs depending on the value of .This has been done for the bottom row in Fig. 7 showing decoded signal.In case of = 0, the amplification is 26.The weakest possible amplification occurs if = ±0.5, and is equal to 19.7 when using this code.Without compensation, periodic ripples will appear in SNR and RCS profiles of the meteor events.This is the cause of the signature reported by Galindo et al. (2011).We discuss this issue further in Sect.5.1.The compensation for the loss in signal power outlined above solves the problem of how to differentiate these signatures from actual physical processes, posed by Galindo et al. (2011).
Figure 8 shows the interpolated range data points of meteor 1 and their residuals compared to a quadratic fit.The residuals have a standard deviation smaller than 0.03 range gates or about 25 m.In the central part where SNR is above 15 dB, the standard deviation goes down towards a hundredth of a range gate or about 10 m.The simultaneously found Doppler data, SNR, RCS and position from interferometry of this meteor are summarized in Fig. 9.

The effect of signal processing on BPSK meteor head echo data
When using a 13-bit Barker code oversampled by a factor of 2, as given in Eq. ( 3), the worst loss of signal due to the cancelling of bauds with opposite phase is 25 %, or about −1.2 dB.If the baud length of a 13-bit Barker code instead is equal to the length of each sample, i.e. if no oversampling is performed, the loss is up to 50 %,or −3 dB.
Our initial analysis of MU meteor head echoes, before we developed this interpolation, resulted in ripples that were identical to those found in JRO head echo observations by Galindo et al. (2011), except for a difference in ripple amplitude due to our oversampling.Galindo et al. describe "a peculiar signature present in SNR plots from meteor-head radar returns".They explain that the signature has "... the following features: (1) strong correlation among fluctuations in SNR values and change in range of a meteor echo, and (2) the fluctuations exhibit periodic ripples with amplitude of 3 dB".Galindo et al. conclude that "... the understanding of this feature is critical to differentiate them from actual physical processes present in meteor returns.Failing to do so could lead to misinterpretation of meteor data." It is apparent from Fig. 1a and c in Galindo et al. (2011) that the systematic drop in SNR appears when the leading edge of the echo is in the middle of two range gates, i.e. when ±0.5.An additional investigation of the JRO decoded signal should show that it becomes asymmetric at the same time as SNR drops, in the manner we described for MU data in Sect. 5 and exemplified in Fig. 7. Galindo et al. (2011) suggest that a possible solution to avoid ripples is increasing the sampling rate with a factor of ∼60 above the transmitter subpulse rate, or from 1 to 60 MHz using their configuration (Chau et al., 2007, Table 1).Knowing the cause of the ripples enables a simple simulation, where we find that this would decrease the amplitude of the ripples to −0.04 dB.This shows that increasing the sampling rate indeed leads to a satisfactory result.However, the interpolation scheme outlined in this paper offers a "cheap" alternative to highly increased sampling, and is in any case advantageous to implement as a complement.It also provides a way to remedy the signal processing issues in already existing data.
The EISCAT meteor code described by Wannberg et al. (2008) and Kero et al. (2008a) is a 32-bit BPSK-coded sequence oversampled by a factor of 4 at reception.Our simulations show that ripples with an amplitude of 13 %, or −0.6 dB should be present in the data.However, the ripple amplitude is small compared to other SNR fluctuations caused by, e.g.fragmentation, quasi-continuous disintegration, etc. (Kero et al., 2008b).Also, the short sampling period of 0.6 µs, which corresponds to range gates with a length of ∼90 m, makes systematic appearance of these ripples rare.This feature has therefore passed unnoticed in the EISCAT UHF observations.
Figure 2 in Galindo et al. (2011) shows that ripples predominantly appear in events of long duration and high SNR.The main reason for this is simply that the ripples are much easier to spot in such events.Hence, it is likely that strong, long-duration events that apparently do not exhibit ripples contain other SNR fluctuations that conceal them, e.g.interference from several meteor targets.

Exclusion of data points
Due to deceleration and the geometry of the meteoroid trajectory, the radial meteoroid velocity component may change more than 10 % over the short time frame of a meteoroid's atmospheric interaction process.For this reason our automatic reduction algorithm must test whether the velocity and range values of consecutive echoes are consistent with a single target or not.This is accomplished as follows: each received radar pulse is analysed separately and its best-fitting Doppler shift, interpolated range gate, signal power and phase, as well as azimuth and elevation angle to the radar target (Sect.8) are stored as one row in a matrix hereafter called the event matrix.An iterative process performs linear least-squares fits on both range versus time and Doppler shift versus time.Residual values more than three standard deviations from either linear fit are excluded from the event matrix and the procedure repeated until there are no more such outliers.Deviating values in either best-fitting velocity, range, or both, caused by simultaneous signals other than the meteor head echo, e.g.echoes from an overdense and enduring meteor trail in a sidelobe, or volume scatter caused by mesospheric turbulence (Reid et al., 1989), or due to an enhanced noise level are thereby excluded.This also provides limits for the temporal extent of the event without having to specify a SNR threshold.
To be able to exclude false rows of data from the initial event matrix but keep those representative of the meteor, we first search for an initial set of data that is likely to represent the meteor.This is accomplished by computing the difference in range and velocity of consecutive rows.As range and velocity in case of a meteor event are estimates of continuous properties, for a row to be classified as representative we restrict the range values of neighbouring rows to be within one range gate and the Doppler velocity not to differ more than ±3 km s −1 .Linear least squares fits are performed on the selected range-time and velocity-time data.Next, the event matrix and these first least-squares fits are exhibited to an iterative procedure which excludes all rows with range values outside three residual standard deviations of the range-time fit, and velocity values outside ±3 km s −1 of the velocitytime fit.The rows remaining after exclusion of outliers are subject to new linear least-squares fits.Range and velocity is again compared to the respective fit and the procedure repeated until no further rows can be excluded.
It is expected that different pulse lengths give different best velocity limits.The velocity limit of ±3 km s −1 is empirically chosen with respect to the random spread of the Doppler data with the described MU radar experimental setup, and deviation of the meteoroid radial velocity from a linear fit of radial velocity due to its non-linear deceleration.The data points remaining after exclusion all have SNR exceeding about −3 dB, which may therefore be regarded as the detectability threshold of the analysis.

Pulse-to-pulse phase correlation
The fraction of a wavelength a target has moved during two adjacent transmissions can most often be determined very precisely by using pulse-to-pulse phase correlation.The main advantage of doing this is the possibility to determine the shape of the meteoroid velocity curve as a function of time (or altitude).This is necessary for dynamical meteoroid mass and atmospheric entry velocity estimations.
The peak of the convolution of the Doppler-shifted version of the transmitted code with the received signal containing a meteor echo (described in Sect.5) is a complex number.Its magnitude provides an estimate of the echo power, amplified from the SNR of each sample by a factor of 19.7-26, or 12.9-14.1 dB, depending on the offset ( ) between target and sampling.
The phase ( ) of the complex number is an estimate of the phase difference between the echo and the Doppler-shifted code.When the same phase is used as reference for analysing consecutive IPPs, their phase difference ( ) can be used to estimate how large fraction of a wavelength the target has moved during the IPP.A meteor head echo target will usually have moved several wavelengths when the IPP is of the order of 1 ms and the radar frequency is in the VHF band or higher (>30 MHz).Any integer number of wavelengths for which the target has moved cannot be revealed by pulse-topulse phase correlation.The velocity from the phase is therefore ambiguous with possible solutions separated according to Eq. ( 1).
The phase values ( ) of each decoded radar pulse of meteor 1 is plotted versus pulse number in the top panel of where p is an arbitrary integer value, λ is the wavelength and T ipp the length of an IPP.The correct (or at least the most probable) value of p is found by comparing the data points of the velocity from the phase (filled circles) with the Doppler velocity data points (open circles) in Fig. 11.For meteor 1 the value is p = 36.If the comparison of ambiguous phase data with Doppler data does not provide a clear distinction, range rate data can be used as a second alternative.The quality of the Doppler data in our analysis procedure is generally always good enough to give a clear distinction when a filtering procedure is used.The solid line in Fig. 11 is a linear least-squares fit to the Doppler velocity implemented for this purpose.
An example of an event with less precise Doppler data than meteor 1 is given in Fig. 12.This meteoroid's initial velocity and deceleration is not possible to determine accurately using the Doppler data alone.However, the Doppler data is good enough to discriminate which of the ambiguous but very precise sets of velocity from the phase that is the most likely one.The velocity determined from the phase reveals how the meteoroid's deceleration increases during the detection and enables dynamical modelling of it's mass loss.The integer number of wavelengths to add to the unwrapped phase difference is in this case p = 55 (yellow circles).
The standard deviation of the velocity determined from the phase as compared to a smooth curve (for this particular example a fourth degree polynomial gives quite random residuals) is 46 m s −1 .The standard deviation of the Doppler data is about 1 km s −1 .
Transmitting radar pulses with unequal IPPs would provide a robust way of unambiguously determining the velocity from phase-to-phase pulse correlation.

Complications in the calculations
A complication that has to be taken into account in order to acquire an accurate velocity from the phase is that the target will generally travel through several range gates.The fastest targets we detect have a radial velocity of about V r = 70 km s −1 .They traverse a range gate in R s /V r 13 ms, that is one every fourth IPP.Each time the leading edge of the echo appears in a different range gate than previously, the phase difference will not provide a correct velocity estimate, but an estimate that is biased by how much the phase changed during one or several sampling periods, depending on how many range gates the target crossed.We have chosen to compensate for this as follows: each IPP is analyzed independently as described in previous sections.To compare a phase value of an echo in IPP = i, where the echo appeared in range gates k to k +26, with a subsequent echo in IPP= i +1 and range gates l to l + 26, we need to estimate how much the phase has changed in the time (l − k)T s .
To do this we use the average Doppler shift ( fD ) from the single-pulse analysis according to where T s is the sampling period.The value of δ is added to the original phase difference when using Eq. ( 8) to estimate velocity.
For targets with a radial velocity of, e.g.V r = 70 km s −1 , the Doppler shift is f D 2f 0 V r /c 0 21.7 kHz when using an operating frequency of f 0 = 46.5 MHz.The phase compensation is in this case δ 0.82 rad or about 0.13 λ when the target passed from one range gate to another (l − k = 1).
For long-duration meteors with large total decelerations, the velocity at any given instant of time may differ from the average velocity by up to about ±5 km s −1 .Such a velocity difference equals a Doppler shift difference from the average of ±1.5 kHz at 46.5 MHz operating frequency.
The phase error (δ error ) introduced by using fD may therefore be up to about thus less than 0.01 λ.This is equivalent to introducing a velocity error of for this particular velocity estimate.An error of the order of V error 10 m s −1 is comparable to or smaller than the standard deviation of the data points of the velocity from the phase (as compared to a smooth curve).Thus, it is small enough not to tamper with further calculations.However, when the final radial velocity is estimated as a function of time, these data points can be recalculated to decrease errors if necessary.

Interferometry
Interferometry calculations are performed on all rows of the original event matrix before the exclusion of data described in Sect.6.We have for this purpose implemented the multiple emitter location and signal parameter estimation (MU-SIC) method developed by Schmidt (1986).It is based on a signal subspace approach suitable for point sources and where the data can be described by an additive noise model (Schmidt and Franks, 1986).Radar studies of meteor head echoes fulfill these criteria.MUSIC, therefore, allows rapid and precise estimations of the signal direction of arrival (DOA).When the criteria are fulfilled, Schmidt (1986) shows that MUSIC can be used to find asymptotically unbiased estimates of, e.g. the number of signals and their DOA for up to K < M multiple source directions, where M in the case of the MU radar is the number of subarrays M = 25.
A comparison of MUSIC with other methods as ordinary beamforming, maximum likelihood and maximum entropy is given by Schmidt (1986).
Guided by Manikas et al. (2001), we have defined an antenna manifold vector ϒ(θ,φ) as where θ is the azimuth (measured positive east of north), φ is the elevation, r = [r x ,r y ,r z ] T ∈ R 3×M are the antenna subgroup centre locations with respect to the geo-metric centre of the whole array expressed in radar wavelengths (and subgroup F5 is located at [0,0,0]), k = 2π[cosφ sinθ,cosφ cosθ,sinφ] T ∈ R 3×1 is the wavenumber vector, is the Hadamard product (elementwise multiplication of the matrices) and γ (θ,φ) ∈ C M×1 is a vector containing the directional gains of the subgroups.The one-way half power beam width of a single antenna subgroup is 18 • .We have in the calculations used unity directional gain for all subgroups, which works well for the purpose of direction finding of targets close to zenith.Furthermore, the MU radar antenna field being horizontally aligned gives r z = 0 and means that Eq. ( 12) can be simplified as ϒ(θ,φ) exp(−2πj (r x cosφ sinθ + r y cosφ cosθ )). (13) The displacements r x and r y of the subgroup centres are illustrated in Fig. 2. The MUSIC spectrum is calculated as where Q n contain the noise eigenvectors.We estimate Q n by first computing a spatial covariance matrix R from the M = 25 set of complex voltages, one from each receiver channel, and each one containing the N = 27 samples selected as containing the meteor echo (as described in Sect.5), according to R = XX /N, where X is a M × N matrix containing the received data.An eigendecomposition of the covariance matrix, [Q,D] = eig(R), gives a set of eigenvectors Q and associated eigenvalues D.
Each present source gives rise to a distinct nonzero eigenvalue D K .If the noise would be zero, there would only be as many nonzero eigenvalues as there are sources (Schmidt and Franks, 1986).Unfortunately noise is seldom zero in an experimental system.Source eigenvalues and noise eigenvalues must therefore be told apart, the former have larger magnitudes.
In our present implementation we are only searching for one point target.We assume that this target gives rise to the largest eigenvalue D max among D and that its associated complex vector Q max therefore defines the signal subspace.The exclusion of Q max and orthogonality of eigenvectors means that the remaining eigenvectors Q n (i.e.all eigenvectors of Q except the one associated with Q max ) now spans the orthogonal complement of the signal subspace, perhaps most appropriately called the signal nullspace (Schmidt and Franks, 1986).When evaluating Eq. ( 14) for different DOA, the denominator will approach zero in the vicinity of the signal DOA and there cause a narrow peak in the spectrum.
To find the DOA of the signal (and thus the direction to the meteor target) we initially evaluate the MUSIC(θ,φ) spectrum of Eq. ( 14) with 5 • steps in azimuth and 0.5 • steps in elevation from 75 • -89.5 • .This gives a densely spaced grid close to zenith where most meteor head echo targets appear.The area around the maximum of this first estimated spectrum is in two subsequent steps evaluated with finer grids, υ 1 Transmitter/Receiver υ 2 M e t e o r o i d t r a j e c t o r y Fig. 13.Exaggerated sketch of a meteoroid trajectory in the far-field of a radar beam (solid lines), where the phase fronts are spherical (dashed lines).The Doppler shift varies as cosυ along the trajectory.For an approaching meteoroid the angle υ 1 < υ 2 , which leads to a decreasing radial velocity component.and the DOA finally searched for within a small fraction of a degree.Each determined DOA is then stored in the event matrix to enable trajectory estimation.

From line-of-sight to vector velocity
It is very important to carefully take the geometrical considerations into account when estimating a meteoroid's velocity from measured radial quantities.This may sound as a superfluous comment, but underestimated flight parameter uncertainties may perhaps explain the anomalous acceleration reported by, e.g.Šimek et al. (1997) and Close (2004).Figure 13 shows an exaggerated sketch of a meteoroid trajectory in the far-field of a radar beam.The angle between the trajectory and the beam is υ, and the radial component of the velocity vector varies as cosυ.Kero et al. (2008a) showed that this gives rise to an apparent acceleration/deceleration term (depending on if the meteoroid is approaching or receding from the radar) that is of the same order of magnitude (and often larger) than the true meteoroid deceleration.
To calculate the meteoroid trajectory, we begin by searching for and including as many successful interferometry data points as possible from each meteor.We start by assuming that the detected trajectory is straight, i.e. the curvature of the trajectory (due to Earth gravity) within the radar beam is negligible.
Azimuth and elevation depend non-linearly on position along a straight trajectory.Therefore, we do all calculations in a Cartesian coordinate system; its origin located at the centre of the MU radar, the x-axis pointing east, the y-axis north South West 85º 86º 87º 88º 89º Fig. 14.Top view of the set of instantaneous target locations at each IPP (green circles) of meteor 1, the start of the event (red star) and a fit of the trajectory (black line).The contours of constant elevation, 85 • to 90 • , (blue) correspond to radial distances of about 1.7 km at 100 km range.and the z-axis completing the set by pointing towards local zenith.Figure 14 shows a top view of the set of DOA of meteor 1 (green circles).The event starts at the red star and a fit of the trajectory is drawn as a black line.
Occasionally, plasma in the trail left behind the meteoroid may constitute a target and interfere with the head echo position determination.To exclude these targets we fit both cartesian coordinates x and y versus time and exclude outliers rather than fitting y versus x.Such plots of meteor 1 is displayed in Fig. 9d.
We are ultimately interested in finding not only the azimuth of the radiant but also its zenith distance.For this reason, we estimate and compare the meteoroid's transversal velocity component obtained from interferometry to its radial velocity component.We proceed as follows.
First we make a linear fit of x versus time and y versus time, using the remaining data points after the iterative procedure described in Sect.6.Then we continue to exclude outliers (data points more than three standard deviation from any of the fits) in a iterative routine until no more points can be excluded.
We assume that the linear fits give reasonable estimates of the transversal velocity components at the central point (p c ) of the detection.For meteor 1, the radar pulse p c = 64 is the central point of the event.Thus, the slopes of the linear fits of x and y give us velocity components v x (p c ) and v y (p c ).If the values of these linear fits at p c are called X(p c ) and Y (p c ), they can together with the very precise range data of the same instant, r(p c ), be used to define a most probable meteoroid position The radial velocity component at the same instant (p c ) is best described by the velocity found using the phase correlation method explained in Sect.7. We use it to estimate the vertical velocity component v z (p c ) according to where θ (p c ) and φ(p c ) are the azimuth and elevation angles to the location of the target at p c as measured from the symmetry axis of the radar.The radiant of the meteor (i.e. the direction from which the meteoroid appears to originate) is expressed in terms of a different set of angles, the azimuth az radiant and the zenith distance zd radiant .These are in horizontal coordinates found by where az radiant is measured east of north, i.e. towards x from y) and arctan computes the arctangent within a range of [−π,+π], and where zd radiant is zero for a meteoroid originating from zenith.
Computing the velocity curve containing also the deceleration is more difficult than making a single vectorial velocity estimate.The trickiest part is converting the accurate radial meteoroid velocity to a reasonably accurate velocity along the trajectory.The reason is that the instantaneous position of the meteoroid (as well as a fit to the position data) has much lower precision than the radial velocity has.Furthermore, the error introduced by assuming, e.g. that the angle to the beam increases linearly (Chau and Woodman, 2004) leads in some cases to an acceleration at the beginning of the event and too fast deceleration at the end, compared to the true values, and in some cases to errors of opposite signs.
To circumvent these problems we use neither a linear assumption on the target angular velocity nor the instantaneous position of the target as a function of time found from interferometry and range, but propagate the target along the determined trajectory (assuming only it is straight) applying the radial velocity itself.For this we only need to use the already defined position P xyz (p c ) and the radiant.If the radial velocity at p c is v radial (p c ) then the meteoroid velocity at that point is where the angle between the trajectory and the line-of-sight vector from the radar to the target is α(p c ).This angle is given by where v xyz is the meteoroid velocity vector.As we have assumed that the trajectory is straight, only the magnitude v met =| v xyz | of the velocity vector changes whereas the direction v xyz = v xyz |v xyz | remains the same throughout the calculations.
The estimated v met (p c ) can now be used to propagate the meteoroid along the trajectory.Its location P xyz at an adjacent time of determined radial velocity is found by multiplying v met (p c ) with the time interval δt (where δt = T ipp if there is a velocity estimate available from the closest possible pair of received radar pulses) and therefore equal to The new position can be used to readily evaluate the new meteoroid velocity from the adjacent radial velocity estimate.We use Eqs.( 18) through ( 20) in an iterative procedure in both directions from the central point (p c ) and thus employ the full precision of the estimated radial velocity to find the meteoroid velocity curve, permitting deceleration and initial velocity to be deduced as accurately as possible.

Error estimation
The largest error in the velocity curve is introduced by the uncertainty of the angle between the trajectory and the beam.
To evaluate how this uncertainty affects the velocity curve we estimate confidence intervals (CI) for the linear fit coefficients of the interferometric data.Simultaneously, this also gives us radiant uncertainty regions.The CI are constructed by calculating the standard errors of the ordinary least squares solutions and multiplying them with the 95 % parameter of the student t distribution (e.g.Hamilton, 1992).Using so determined CI for both zenith distance and azimuth we construct an elliptical area (circular if the uncertainties in both directions are equal), which contains the true meteoroid radiant with 95 % certainty under the condition that the residuals of the interferometry data are random and normally distributed.To find boundaries for the velocity curve we apply the iterative process described above but with v xyz of Eq. ( 20) replaced by vectors corresponding to the smallest and largest angle to the beam within the radiant area.Meteoroid velocity curves for meteor 1 computed in this manner are plotted as dotted lines in Fig. 9c.Its initial velocity is 55.5 ± 0.2 km s −1 and the deceleration does not change significantly within the estimated uncertainty region.
Because a meteoroid's initial deceleration can be very small, an error as little as of the order of 1 • can indeed sometimes make a meteoroid appear to be accelerating along some part of, or the whole, detected trajectory.Exemplified in Fig. 15 is a meteoroid on a trajectory that crossed the beam at an angle greater than 80 • .The CI of the angle to the beam is 0.3 • .Yet, the amount of deceleration/acceleration exhibited during detection in the radar beam cannot be determined.Nevertheless, it is worth noting that its initial velocity can be determined to a reasonable degree of accuracy, to 29.2 ± 1 km s −1 .In fact, it can be limited even further, to 28.6 ± 0.4 km s −1 , if deceleration is presumed.The closer a meteoroid trajectory is to perpendicular to the beam, the more sensitive is the deceleration determination to errors.However, as the deceleration initially can be very small, an overestimated angle to the beam may cause meteoroids on all slant angles to appear to be accelerating.Conversely, if the angle to the beam is underestimated, a meteoroid will appear to decelerate faster than it does.The latter is an error less likely to be noticed, as meteoroids are expected to decelerate.Nevertheless, mass calculations based on the standard momentum equation (Bronshten, 1983, p. 12), using the velocity (v) and deceleration ( v) obtained from an event which angle to the beam is underestimated will result in an underestimated meteoroid mass.When the cross-sectional area of the meteoroid is rewritten using an arbitrarily chosen meteoroid shape factor (Bronshten, 1983, p. 14), it is easily seen that meteoroid mass is proportional to v 6 / v3 (Campbell- Brown and Koschny, 2004, Sect. 2.4 and Eq. 2).Small errors in v and v, therefore, quickly cause large errors in estimated mass (Kero et al., 2008a).

Radar cross section
Radar cross sections (RCS) of detected targets are evaluated by rewriting the classical radar equation (e.g.Skolnik, 1962) as where P r = received power, R = target range, G t = transmitter antenna gain, G r = receiver antenna gain, θ = azimuth of target (positive east of north), φ = elevation of target, λ = radar wavelength, and P t = transmitted power.
The received power is given by where SNR is the signal-to-noise ratio, T noise is the equivalent noise temperature, k B = 1.38 × 10 −23 J K −1 is the Stefan-Boltzmann constant, and b w 1/6 µs ≈ 167 kHz is the receiver bandwidth.T noise = T sys + T cosmic is the sum of the system noise (T sys ∼ 3000 K) and the cosmic background radio noise that varies from about T cosmic ∼ 5000-15 000 K throughout one diurnal cycle.T cosmic is dominated by the passage of two strong radio sources close to zenith, Taurus-A and Cygnus-A.Except for the receiver noise temperature and noise contribution due to losses in feed, T sys may in this context also include contribution from atmospheric emission and ground radiation (spillover and scattering).To proceed, we assume that T sys is constant throughout each diurnal cycle.This is not necessarily true, but as long as both T sys and its variance are small with respect to T cosmic , the assumption is not of major concern in the estimation procedure, as exemplified in Fig. 16 described below.24).Here, g = 7.5 × 10 −6 and T sys = 2500 K.
T cosmic towards zenith above the MU radar site is known, and varies by almost a factor of three during the course of one day.Given that the receiver response is linear and zero biased, we thus have a possibility to easily find the unknown T sys and the conversion factor g between the digital receiver output unit and power.We search for the best fit between a theoretical daily variation, P test , and the observed background noise level while varying g and T sys .Figure 16 displays an example of such a comparison.As the parameters g and T sys may change, appropriate parameters should be adapted for each observation.Another way to calibrate the parameters than described above, is to point the antenna towards known celestial sources and/or inject known amounts of noise at the antenna level of the receiver system.
The peak transmitted power is optimally P t = 1 MW, but is not continuously monitored and is therefore an estimated value.The attenuation of the signal as compared to a target at the boresight axis where G t = G r 34 dB is found by combining the radiation pattern of a single crossed-Yagi antenna with the array pattern illustrated in Fig. 2 (Fukao et al., 1985).Since the array is nearly circular the sidelobes are fairly symmetric in the azimuthal direction.
We have implemented a beam-steering algorithm where we utilize the linear fit of the target position to interpolate (and extrapolate) the expected target direction throughout the event matrix and calculate appropriate phase offsets which are added to each set of complex data from the receiver channels.The data from all channel are thereafter added coherently and the event reanalyzed using the same algorithms as already described (except interferometric calculations).This leads to extended durations of many events as compared to the original analysis, a natural consequence of the increased gain far from bore axis and the duration of many meteors being limited by the spatial extension of the detection volume rather than onset and/or end of ionization along the meteor trajectory.
Comparing the RCS estimated with and without postbeam steering of the receiver beam when the target is close to a minimum in the radiation pattern gives an indication of the validity: the difference between the two RCS estimates shows how well the position of the target is determined or, alteratively, how much the theoretical beam pattern used in the calculations deviates from the true radiation pattern close to the target.Near the centre of the main beam where the gain is close to the boresight axis value of 34 dB, the deviation is always small.At any rate, the RCS values where the two estimates differ should be discarded.
Investigations of signatures in SNR/RCS profiles (Kero et al., 2004;Kero et al., 2005;Kero et al., 2008b;Janches et al., 2009;Mathews and Malhotra, 2010;Mathews et al., 2010;Malhotra and Mathews, 2011) suggest that fragmentation and differential ablation are ubiquitous features of radar meteor results.Such signatures are present also in the MU radar SNR/RCS profiles and will be investigated in the future.Briczinski et al. (2009) report that the deceleration is not possible to determine for a large fraction of the meteor head echo events observed with the 430 MHz AO radar, and conclude that the reasons for this are primarily low SNR, fragmentation and short duration.In addition, we also note that even if the deceleration of a target seems well-defined, it cannot necessarily be interpreted using single body ablation theory (and easily converted to mass) without modelling the effects of fragmentation.
Typical RCS distributions of 10 000 sporadic meteors and 600 Orionid shower meteors observed by the MU radar using the methods described in the current paper are published in Fig. 12 of Kero et al. (2011).The distribution peaks at RCS 10 −3 m 2 .At beam centre, meteors down to RCS 3 × 10 −5 m 2 were detected.A handful of meteors had RCS > 1 m 2 .
As a comparison, Zhou et al. (1998) found that head echoes recorded with the 46.8 MHz Arecibo VHF radar typically had RCS of the order of 10 −3 m 2 , while the simultaneous 430 MHz UHF echoes were of the order of 10 −6 m 2 .Mathews et al. (1997) used the Arecibo 430 MHz UHF radar and found echoes with smaller scattering cross sections, RCS 10 −8 m 2 .Close et al. (2002) used the 160 MHz VHF and the 422 MHz UHF ALTAIR radars to study meteor head echoes.ALTAIR transmits RC polarization and can receive both RC and LC polarized signals.In order to compare with MU data, we focus on the LC polarized data.The LC RCS of AL-TAIR VHF meteors ranged from approximately 10 −4 m 2 to 10 −1 m 2 .The ALTAIR UHF LC RCS distribution was between approximately 10 −7 m 2 and 10 −3 m 2 .
These studies show that meteor head echo target RCS appears to be frequency dependent.Another conclusion that can be drawn is that the smallest RCS recorded by different systems largely depends on the sensitivity of the radar, i.e. the limit where the meteor signal is no longer strong enough to be analysable.
However, the largest RCS targets recorded with the MU radar are considerably larger than those of ALTAIR.This difference could partly be due to the different frequencies (AL-TAIR VHF 160 MHz versus MU 46.5 MHz), but also due to the difference in probabilities of larger than dust-sized meteoroids entering the radar beam during observation, i.e. the product of collection area times the observation time.The MU radar antenna gain pattern has relatively strong sidelobes, enabling detections of head echo targets with RCS > 1 m 2 .This detection volume has a horizontal cross-sectional area of the order of 1000 km 2 at 100 km altitude when the beam is pointed vertically (Kero et al., 2011).The observation period dedicated for RCS determination covered 33 h and resulted in ∼10 000 meteors.The ALTAIR VHF observation consisted of 734 meteors recorded during 29 min.This could explain why rarely occurring large RCS targets were not present in the ALTAIR data.
Moreover, the large RCS targets observed with MU agrees with the head echo target sizes that were recorded with the 32 MHz radar at the Springhill Meteor Observatory for meteors also observed visually (Jones and Webster, 1991).Jones et al. (1999) point out that such large targets are possible to detect even with a modest power VHF radar, but that they appear at a rate of less than five per day and are therefore generally ignored.

Channel phase offset compensation
The 25 channel setup of the MU radar enable ample opportunities for interferometric calculations as well as interchannel calibration when a point target is present in the beam.We have adapted a simple phase calibration algorithm for head echo targets.Since the head echo observations are run on campaign basis and not continuously, we save and keep all raw data and are consequently fortunate in having the possibility to rerun complete analyses on the whole data set, even after an observation has ended.The purpose of the first round of analysis is only to derive interchannel calibration coefficients, which are then used in the final analysis.From the results of the initial (uncalibrated) analysis, we select a set consisting of at least a few tens of well-behaved meteor events (high SNR, no apparent interference or contribution from trail plasma etc.) spread out in time during the usually 24 h long measurement.The selected events are examined one by one in an automated phase-error search routine.Using the estimated target direction at each IPP of an event, we calculate what the phase difference between the signal received at each subarray should be and compare with the actual differences in phase between the subarrays.This gives a set of phase offset values that we apply on the original data of the event, and reevaluate the target direction in an iterative procedure until the offset values found in an additional iteration are infinitesimal.The calculated phase offset values from the iteration procedure form a converging series and less than ten iterations are generally enough.This iterative procedure is applied once on each of the selected events.
If consistent phase offset values are found throughout the whole set of events, we calculate average offset values and input these to the final analysis of the complete data set of that observation period.At the MU radar, there is typically one channel (channel 9 = C1), which has an offset of 0.5-0.6 rad (i.e. 30 • ), while the rest are of the order of 0.1 rad.During one measurement occasion (January 2010), data from one channel contained glitches in complex amplitude and deviating phase values due to an unknown T/R hardware problem.This problem was easily identified and solved by assigning the phase offset value a weight of zero for that particular channel and data set for the reanalysis, effectively cancelling the contribution of the faulty channel on the final analysis results.
12 Validity control Many flagged events are not meteor head echoes but caused by, e.g.interference or echoes from meteor trail plasma.To ensure a high-quality data set only containing genuine meteor head echoes, we apply a qualitative data reduction technique.One criterion is agreement between the two independent estimates of the target radial velocity: the velocity derived from the Doppler shift of the received signal and the target range rate.As a first attempt, we demanded these to agree within ±2 km s −1 .The second criterion is the estimated trajectory uncertainty being within ±2 • and the estimated velocity vector uncertainty being smaller than ±2 km s −1 .Some echoes from drifting meteor trails, Earth-orbiting satellites and space debris fulfill the above criteria.However, these targets all have vector velocities smaller than the Earth escape velocity, ∼11.2 km s −1 , and are therefore easily removed.

Conclusions
We have developed an automated analysis method for meteor head echoes observed by the interferometric MU radar near Shigaraki in Japan.In this paper we focused on reporting the algorithms of the method in detail.
We have shown that the precision improvement of the radial velocity is at least a factor of 20 when using pulseto-pulse phase correlation in combination with single-pulse Doppler measurements compared to single-pulse Doppler measurements alone.The improved measurements reveal that deceleration increases significantly during the intense part of the meteoroid-atmospheric interaction process, which is equivalent to rapid mass loss.
When using an interpolation scheme in the decoding procedure of a transmitted 13-bit Barker code oversampled by a factor of two at reception, it is evident that the range of meteor targets can be determined to within a few tens of meters, or of the order of a few hundredths of the 900 m long range gates, at each received pulse.Also, we have identified and solved the signal processing issue giving rise to the peculiar signature in signal to noise ratio plots reported by Galindo et al. (2011), and show how to use the range interpolation technique to differentiate the effect of signal processing from actual physical processes.
We have developed a careful method of error estimation and exemplified that small errors in meteoroid trajectory parameters may lead to anomalous acceleration or overestimated deceleration.Underestimated flight parameter uncertainties may perhaps explain the anomalous accelerations reported by, e.g.Šimek et al. (1997) and Close (2004).Proper error estimation is crucial for the credibility of using meteor head echo measurements for meteoroid mass, velocity and orbit distributions and for developing models of the head echo radio wave scattering process.
Using the analysis methods described in the current paper, we have from June 2009 to December 2010 (except for August 2009) carried out monthly 24 h or longer meteor head echo observations using the MU radar.An overview of the complete set of data (>500 h) containing >100 000 highquality meteor events is given by Kero et al. (2012).A survey of 10 000 meteors collected during 33 h of observation in October 2009 is reported by Kero et al. (2011).In that paper we estimated an equivalent MU radar collection area by examining the detection probability of a meteor as a function of its maximum RCS and location within the MU radar detection volume.Furthermore, we used the collection area to convert the detection rates of the 2009 Orionid meteors to a meteoroid influx, comparing it to the Orionid flux estimated from visual and specular meteor radar detection rates.The MU radar Orionid radiant density distribution is as compact as the distribution of all precisely reduced Orionids (66 photographic and 19 video meteors) of the IAU Meteor Data Center presented by Lindblad and Porubcan (1999).This indicates that the meteor head echo radar method described here provides precision and accuracy comparable to the photographic reduction of much brighter meteors with longer detectable trajectories.

Fig. 3 .Fig. 4 .
Fig. 3.A representation of a 13-bit Barker code oversampled by a factor of two used in the described MU radar experimental setup (upper panel) and its ACF (lower panel).

Fig. 5 .
Fig. 5. Doppler velocity versus time-of-flight velocity for >100 000 MU radar meteors.The solid line is a linear least-squares-fit with a slope of 0.998.

Fig. 6 .
Fig. 6.Range-time and signal intensity plot of a meteor head echo event detected 28 July 2009, 05:33:09 JST, in subsequent figure captions and the text referred to as "meteor 1".

Fig. 8 .
Fig. 8. Upper panel: interpolated range data (open circles) of meteor 1 and a quadratic fit (solid line).Lower panel: the residuals have a standard deviation of less than 0.03 range gates or about 25 m.In the central part where SNR > 15 dB, the standard deviation goes down towards a hundredth of a range gate or about 10 m.

Fig. 9 .
Fig. 9. Overview of the analysis parameters of meteor 1: (a) range, (b) radial velocity, (c) geocentric meteoroid velocity, (d) transversal displacement from beam centre in west and south directions, (e) angle of trajectory to the beam, and (f) RCS and equivalent signal temperature (T signal = SNR • T noise , where T noise 10 4 K).Blue curves represent parameters obtained with, and red curves without, post-beam-steering throughout all panels.The green markers in (b) trace the velocity determined from pulse-to-pulse phase correlation, presented in greater detail in Figs. 10 and 11.The dotted lines in panels (c) and (e) show the estimated 95 % CI of the meteoroid velocity uncertainty margin and the angle to the beam, respectively.

Fig. 10 .
Fig. 10.From top to bottom: phase values ( ), phase difference of consecutive radar pulses (), and unwrapped phase difference, all versus radar pulse number of meteor 1.

Fig. 10 .
Fig. 10.The middle panel shows the difference in phase ( ) found by comparing consecutive IPPs.The bottom panel presents an unwrapped version of the phase difference.The unwrap procedure is a search for a smooth phase curve by adding or subtracting integer values of 2π to each value of shown in the middle panel.In Fig. 11 the calculated phase curve is converted to radial velocity according to

Fig. 12 .
Fig. 12. Unwrapped phase difference (filled circles) and singlepulse Doppler data (open circles) of a meteor detected 28 July 2009, 05:33:08 JST.The integer number of wavelengths to add to the unwrapped phase difference is in this case p = 55 (yellow circles).The standard deviation of the velocity determined from the phase as compared to a smooth curve (for this particular example a fourth degree polynomial gives quite random residuals) is 46 m s −1 .The standard deviation of the Doppler data is about 1 km s −1 .

Fig. 15 .
Fig. 15.Overview a meteor detected 14 December 2010, 00:01:29 JST: (a) range, (b) radial velocity, (c) geocentric meteoroid velocity, (d) transversal displacement from beam centre in west and south directions, (e) angle of trajectory to the beam, and (f) RCS and equivalent signal temperature (T signal = SNR • T noise , where T noise 10 4 K).Blue curves represent parameters obtained with, and red curves without, post-beam-steering throughout all panels.The green markers in (b) trace the velocity determined from pulse-to-pulse phase correlation, presented in greater detail in Fig. 11.The dotted lines in panels (c) and (e) show the estimated 95 % CI of the meteoroid velocity uncertainty margin and the angle to the beam, respectively.

Fig. 16 .
Fig.16.The average power level (gray) of each MU radar data block (512 IPPs 1.7 s of data) during an observation, given in receiver output units.Data blocks without meteor events (or other strong signals) trace out the bottom part of the gray curtain-like shape and are representative of T noise .Theoretical values (black dots), are calculated according to Eq. (24).Here, g = 7.5 × 10 −6 and T sys = 2500 K.