Improved method of estimating temperatures at meteor peak heights

For 2 decades, meteor radars have been routinely used to monitor atmospheric temperature around 90 km altitude. A common method, based on a temperature gradient model, is to use the height dependence of meteor decay time to obtain a height-averaged temperature in the peak meteor region. Traditionally this is done by fitting a linear regression model in the scattered plot of log10(1/τ) and height, where τ is the half-amplitude decay time of the received signal. However, this method was found to be consistently biasing the slope estimate. The consequence of such a bias is that it produces a systematic offset in the estimated temperature, thus requiring calibration with other co-located measurements. The main reason for such a biasing effect is thought to be due to the failure of the classical regression model to take into account the measurement error in τ and the observed height. This is further complicated by the presence of various geophysical effects in the data, as well as observational limitation in the measuring instruments. To incorporate various error terms in the statistical model, an appropriate regression analysis for these data is the errors-in-variables model. An initial estimate of the slope parameter is obtained by assuming symmetric error variances in normalised height and log10(1/τ). This solution is found to be a good prior estimate for the core of this bivariate distribution. Further improvement is achieved by defining density contours of this bivariate distribution and restricting the data selection process within higher contour levels. With this solution, meteor radar temperatures can be obtained independently without needing any external calibration procedure. When compared with colocated lidar measurements, the systematic offset in the estimated temperature is shown to have reduced to 5 % or better on average.


Introduction and background
As meteoroids enter the Earth's atmosphere, they produce ionized trails which can be detected as back-scattered radio signals by interferometric radars. After the trail has been formed, the ionisation begins to dissipate by various processes, such as, ambipolar diffusion, eddy diffusion, or electron loss due to recombination and attachment depending on the height of ablation. The rate at which the echo power decreases is also determined by the combined effect of electron line density of the trail, ambient pressure and temperature.
If the electron line density of the trail is less than 2.4×10 14 electrons m −1 the trail is called 'underdense', meaning each electron in the trail scatters independently. The decay of underdense trails is thought to be mainly due to ambipolar diffusion at a height range of 85-95 km, where the majority of the meteors ablate (Jones, 1975). In the weak scattering limit the backscattered 25 amplitude of the radio signal from an underdense trail decays with time (t) as where λ r is the radar wavelength and D a is the ambipolar diffusion coefficient (Kaiser, 1953). This coefficient depends on the ambient pressure (P ) and temperature (T ) of the neutral gas (Chilson et al., 1996) and can be estimated from the half-amplitude decay time (τ ) as 30 D a = K amb T 2 P = λ 2 r ln2 16π 2 τ (2) K amb in Eq.
(2) is a constant related to the ionic constituent of the plasma in the trail (Hocking et al., 1997). The pressure at a given height (h) is where m is the mass of a typical atmospheric molecule, g is the acceleration of gravity, k is the Boltzmann constant, z is an axis along the vertical. Substituting the equation for pressure in Eq.
(2), and differentiating provides the height profile of the decay time: log 10 D a (h) = 2log 10 T (h) + log 10 e mg k where Ψ is a constant. Equation (5) states that the height profile of decay time is a function of both temperature and temperature 40 gradient under the assumption of ambipolar diffusion for underdense meteor trails. In practice, most trail echoes are received at a small altitude range referred to as the region of peak meteor occurrence (Hocking, 1999). Hence a height-averaged temperature gradient near the peak height can be used to estimate the mean temperature (T 0 ) at the peak height by fitting a linear function (Hocking, 1999). A linear approximation of Eq. (5) is, 45 where β is the slope of the scattered plot of log 10 (1/τ ) and height, and T 0 is the average temperature of the atmosphere at the height of peak meteor occurrence. A typical scattered plot of height and log 10 (1/τ ) shows significant variation in the measured data along both abscissa and ordinate (Fig. 1). Traditionally, the slope (β) is estimated using the ordinary leastsquares (OLS) method with log 10 (1/τ ) as the independent variable. The justification of using log 10 (1/τ ) as independent variable is that the measurement errors in τ are smaller than those in heights (Hocking et al., 1997). While the pulse length, 50 angular resolution etc., of the radar introduces intrinsic measurement errors in heights, much of the variation in decay time is due to various geophysical effects that persist at all altitudes. At higher altitude the collision frequency with neutrals is reduced, and the diffusion is inhibited in a direction orthogonal to the geomagnetic field (Jones, 1991;Robson, 2001). This anisotropic diffusion causes an increase in the duration of meteor radar echoes as compared to ambipolar diffusion. Whereas at lower altitude decay time tends to decrease due to additional effect of electron loss by recombination and attachment (Younger et al., 55 2008). In addition, other geophysical factors, such as meteor fragmentation, turbulence within trail, chemical composition of the meteors, or the temperature variation due to passage of tides and gravity waves, can contribute to the measurements of decay time and heights at all altitudes (Hocking, 2004).
Temperature estimation from meteor radar (MR) data requires obtaining the best-fit regression line in the scattered plot of log 10 (1/τ ) and height. However, the pioneer work done by Hocking (1999) to implement this method using ordinary least-60 squares fitting showed a clear systematic offset between the MR temperature and colocated lidar measurements, indicating that the estimated slope was not determined correctly. To correct for this offset, a common practice is to calibrate the meteor radar temperatures using temperatures from lidar, OH spectrometer, satellite or rocket climatology. Hocking et al. (2001b) provided a statistical comparison technique (SCT) to calibrate the biased slope estimate as, Or, where s δ and s ε are the error variances of the (log of) diffusion coefficient (or decay time) and height respectively, s d and s h are data variances of log 10 (1/τ ) and of height respectively, β d OLS and β h OLS are slope coefficients when diffusion coefficient or height is treated as independent variable respectively in the OLS regression analysis. The calibrated slope, β SCT , is traditionally 70 obtained by arbitrarily choosing s δ or s ε that gives the best temperature estimates of MR as compared to optical or satellite data (e.g., Holdsworth et al., 2006;Hocking et al., 2007;Kim et al., 2012). A severe shortcoming of such calibration procedure is that the calibration factors (s δ and s ε ) in the parenthesis of Eq. (7) or Eq. (8) are dependent on the data selection criteria (e.g., limiting heights, decay time or zenith angle to a certain range). Moreover, the outcome of any such calibration routine will depend on the location of the MR and the choice of the calibration instrument. From a pure statistical context, the arbitrary 75 choice of calibration parameters makes the estimated temperature also an arbitrary quantity thereby making it impossible to draw any reasonable statistical inferences.
In practice, the ordinary least-squares method will not be valid for MR data since neither the height nor the decay time can be predetermined as independent variable, and both variables are subjected to intrinsic measurement errors and various geophysical effects. The reasons for such bias, and thus needing calibration, is discussed on theoretical and experimental 80 grounds in Sect. 3.1. In addition, a statistical procedure to estimate s δ and s ε using SCT calibration is formulated and presented in Sect. 3.1. An alternative method that includes measurement errors in the regression model is introduced in Sect. 3.2.1. This analysis does not require an absolute knowledge of s δ or s ε , but only the relative ratio is needed. As an alternative to SCT calibration, in Sect. 3.2.2, we have provided an analytic solution for estimating the bias-corrected slope coefficient. Using this solution and Eq. (6), an independent and absolute value of MR temperature is obtained. A comparison study of the estimated 85 MR temperatures with colocated lidar temperatures is discussed in Sect. 4.

Instrumentation and data
The All-Sky Interferometric Meteor Radar (SKiYMET) at Sodankylä Geophysical Observatory (SGO,67 • 22' N,26 • 38' E,Finland) has been routinely monitoring daily meteor-height averaged temperatures and wind velocity since December 2008 (Kozlovsky et al., 2016). The radar operates at a transmission power of 15 kW and frequency of 36.9 MHz, with a transmitting 90 antenna which has a broad radiation pattern designed to illuminate a large expanse of the sky. The meteor trails are detected within a circle of 300 km diameter around SGO. The phase differences in the five-antennae receiving array allow the determination of the azimuth, elevation, range, and line-of-sight Doppler velocity of the meteor trails. The 2144 Hz pulse repetition frequency of MR transmission introduces a range ambiguity of 70 km, and the built-in analysis software therefore assumes meteor trails are within the height of 70 to 110 km for unambiguous detections. Uncertainty in the height is ±1 km (or better 95 for large zenith angle), which is determined by the 2 km range resolution. In addition, the half-time (τ ) of the received signal is calculated from the width of the autocorrelation function. A detailed description of the algorithm of the SKiYMET signal processing software is outlined in Hocking et al. (2001a).  SGO is located at the corrected geomagnetic latitude of 64.1 • , which is statistically a region of the auroral oval. Hence the radar frequently detects nonmeteor targets during substorms associated with ionospheric plasma waves generated due to 100 Farley-Buneman instability (Kelley, 2009). The Doppler velocity of such echoes can be more than a few 100 ms −1 , which are mostly detected at low elevation (Lukianova et al., 2018). MR also detects ground echoes modulated by the ionosphere during pulsating auroras (Kozlovsky and Lester, 2015). These targets have near-zero Doppler velocities, and are also observed at low elevation. Our data selection criteria is kept to bare minimum, such that, all heights and decay times are included as long as they are unambiguous detections above 30 • elevation angle and Doppler radial velocity in the range ±100 ms −1 . The Mass-Spectrometer-Incoherent-Scatter or MSIS90 (Hedin, 1991) model temperatures are used to generate a temperature gradient model near the peak heights. Model temperatures are computed for each date at intervals of 6 h between 85 km to 95 km. A third-degree polynomial fit is carried out to obtain the height profiles at 6 UT, 12 UT, 18 UT and 24 UT. For each time interval, the gradient at the respective meteor peak height, as well as at 1 km above and below from peak height are estimated.

115
These 12 values are then used to obtain the mean and standard deviation of the temperature gradient for each day near the peak height, which varies in the range 89±1 km. The daily meteor detections, the peak heights, and the corresponding temperature Our reasons for choosing the MSIS temperature gradient are twofold. Firstly, MSIS data are easily accessible from the online 120 version, which guarantees reproducibility of this work independent of location. Secondly, even if the temperature gradient term in Eq. (6) is ignored, the resulting offset in the estimated temperature can be no more than 10% (Hocking, 1999). Hence an approximate estimate is sufficient for the main objective of this paper. However, the actual temperature gradient in the atmosphere may be slightly different from these model values, which can contribute to the biasing effect in the estimated temperatures. Any such possibility and its effect on the estimated temperatures is addressed in the subsequent section. In the following text we use the notation and formulation in Gillard and Iles (2005). The observables, log 10 (1/τ ) and height, are represented as d i and h i respectively, and the corresponding unobserved true values as ξ i and η i respectively, where the index i represents the i'th meteor detection. For consistency we also assume that log 10 (1/τ ) is presented in abscissa and height 130 is in the ordinate in the respective scattered plot (as shown in Fig. 1). Suppose, we are assuming a linear relation in variables ξ i and η i as Due to measurement errors and various geophysical processes, the true values ξ i and η i will be subjected to random errors and 135 hence the observable d i and h i will have scatter around the linear model in Eq. (9): And, where δ i and ε i are errors in the measured log 10 (1/τ ) and height respectively, and are assumed to be mutually uncorrelated, 140 have zero mean and independent of the suffix i. This implies that the measurement error variances, s δ and s ε , are constant with respect to the suffix i. Classical regression analysis or ordinary least-squares (OLS) treats d i as an independent variable without intrinsic errors (or, d i = ξ i ), and then minimises the sum of squared residuals along the ordinate. The slope from the OLS method can be represented in terms of the covariance (Cov) and variance (V ar) as (e.g., Smith, 2009;Keles, 2018), where β d OLS is the OLS slope estimated by considering log 10 (1/τ ) as independent variable, and r is the Pearson productmoment correlation coefficient between d and h. Likewise, by reversing the arguments above, it's trivial to show that the reciprocal value of OLS slope estimate with height as independent variable is (e.g., Smith, 2009), To see the effect of errors in the independent variable on the OLS slope estimate, e.g., for β d OLS , Eq. (10) and Eq. (11) are used 150 in Eq. (12), Since δ i , ε i and ξ i are mutually independent, Eq. (14) simplifies to , where ζ is known as the attenuation or regression dilution bias. Since variances are always positive by definition, Eq. (15) 155 shows that in the presence of measurement error in the so-called independent variable (in abscissa), the OLS slope estimate (β d OLS ) will always be smaller than the unbiased slope β. Likewise, β h OLS is greater than β if there is error in the measured height (specific example presented in Fig. 1a). By substituting d i = ξ i + δ i , we can rearrange Eq. (15) as . Statistical distribution of (a) log10(1/τ ) and (b) height obtained using SCT method with colocated lidar data. The variance ratios of these parameters in the normalised coordinate system is shown in (c). For comparison, (d) presents these variance ratios obtained using EIV method (Sect. 4) for the same data set. The properties of these histograms are presented in Table 1.
Equation (16) is a well known identity (e.g., Carroll and Ruppert, 1996;Frost and Thompson, 2000) in statistical literature, 160 that was re-derived by Hocking et al. (2001b) in Eq. (7) in the context of SCT correction. Equation (16) reveals that an absolute knowledge of error variances, s δ ( or s ε ), is required to obtain the bias-corrected slope (β) if we choose OLS fitting for the slope estimate. To this date, no such attempt has been made to assess these error variances in MR data. Instead, biased temperatures were calibrated with optical or satellite data by arbitrarily choosing a value of s δ or s ε for estimating the so-called SCT slope, β SCT (e.g., Holdsworth et al., 2006;Hocking et al., 2007;Kim et al., 2012). In the remaining part of this section, we 165 demonstrate how to obtain a representative value of the error variances in this data following a revised calibration procedure.
For each 24 h of the data set, we performed two OLS fittings to estimate β d OLS and β h OLS by using Eq. (12) and Eq. (13). The corresponding biased temperatures, T d M R and T h M R respectively, are estimated using Eq. (6). Experimental values of the parameters s δ and s ε can be obtained by comparing these estimated biased temperatures with the colocated lidar temperatures (T lidar ) as the reference values. Using Eq. (7) and Eq. (8), and noting that the slope is proportional to the estimated temperatures 170 from Eq. (6) we obtain: where T d M R and T h M R are MR temperatures estimated using OLS fitting with log 10 (1/τ ) and height as independent variable respectively. Furthermore, if the measurements are normalised with the mean and standard deviation (STD) as (18)): where the error variances, s ε and s δ , are in the dimension free system defined by Eq. (18). In essence, λ OLS ef f is a measure of all sources of errors in the normalised heights and decay times that cause the real data to deviate from the idealized physical 180 model of Eq. (6), thereby producing a typical scatter as seen in Fig. 1.
λ OLS ef f , s δ , and s ε were estimated by Eq. (17) and Eq. (19) using 24 h of MR data and colocated lidar temperatures at 88 km, 89 km and 90 km for the dates for which lidar data was available. The biasing effect on the OLS estimate of MR temperatures with log 10 (1/τ ) and height as independent variable respectively are presented in Fig. 3 (green and blue lines and histogram).
As expected from Eq. (7) or Eq. (8), a mean offset of -78 K and +367 K occurs depending on whether log 10 (1/τ ) or height 185 respectively is considered as independent variable. In practice, the magnitude of these biases is related to the total errors in heights and log 10 (1/τ ) from Eq. (17), which was not taken into account by the OLS regression model in Eq. (12) and Eq.   (17), we have used the statistical mean of differences for calibration. This is because lidar data is not available for all days during the 6 months of data used in this work. Moreover, both MR and lidar data have their own intrinsic 195 errors and technical differences in the observation time and volume of sky. MR temperatures are daily averages over 24 h of observation, whereas lidar data is just the nightly mean profile. The lidar probes a small volume limited to the diameter of the lidar beam, while the radar illuminates a large part of the sky. For a single observation, lidar may see the phase structure of large-scale gravity waves while the MR averages over the gravity wave structure due to different spatial resolutions. As a result, the radar averages over gravity waves with horizontal wavelengths smaller than few hundred km. On the other hand, 200 the lidar may resolve these gravity waves if the run time is shorter than the period of these waves. As gravity wave amplitudes are on the order of 10-15 K at these altitudes, we cannot expect perfect agreement between radar and lidar temperature due to geophysical variations which show up differently in the two data sets as a result of the different observational volumes.
While such calibration routine may prevent large offsets in the estimated temperatures, the day-to-day variation in these error variances (as shown in Fig. 4) due to natural geophysical processes will persistently introduce artefacts in the estimated 205 temperatures. Moreover, due to the continually changing atmospheric dynamic, these calibration parameters need to be updated at time intervals. This in turn requires availability of optical or satellite data throughout the year. In the following section, we provide a robust method of estimating MR temperatures that do not require any external calibration. Table 1. Representative value of the (square root of) error variances and their normalised ratio for SGO's meteor radar obtained by SCT method (Fig. 4)

Errors-in-Variables (EIV) model
For fitting a straight line model, such as, to a set of N data points (x i , y i ) measured with errors, the corresponding χ 2 merit function is (Press et al., 1992, p 660) where σ yi and σ xi are the standard deviation of the i'th data point, and the weighted sum in the denominator of Eq. (21) can be interpreted as the weighted-error of the i'th data point. The regression coefficients, a and b, can be found by minimising the 215 merit function with respect to these coefficients following any suitable numerical root-finding routine. However, under certain assumptions on the error variances, it is possible to derive an analytic solution for the regression coefficients. This analytic solution is introduced and discussed in Sect. 3.2.1. And in Sect. 3.2.2, we have looked more closely into the nature of various error variances in the data to construct a closed-form solution of the bias-corrected slope coefficient.

First estimate of slope (β)
220 Application of Eq. (21) to physics data requires that all measured variables are dimensionally consistent so that χ 2 is dimension free. Moreover, the analysis in this section requires that the measurements are presented in an appropriate dimension free system. This facilitate the direct comparison between different parameters, such as, the measured variables or the associated error variances. By applying the coordinate transformation introduced in Eq. (18) to Eq. (9), we therefore intend to solve the simplified bivariate linear system of equations, where η i , ξ i and β W are dimension free. For the specific choice of normalisation by Eq. (18), the intercept (α W ) is always zero in the transformed coordinate system. The merit function in Eq. (21) can be further simplified by invoking a homoscedastic standard weighting model (Macdonald and Thompson, 1992). This error model assumes that the error variances are independent of data point, thereby simplifying the merit function as (Macdonald and Thompson, 1992;Lolli and Gasperini, 2012), where s δ and s ε are constant error variances of the measured log 10 (1/τ ) and heights respectively in the normalised (or, dimension free) coordinate system, and λ is the ratio The χ 2 minimisation of Eq. (23) with respect to β W leads to the well-known (Carroll and Ruppert, 1996;Smith, 2009;Lolli and Gasperini, 2012) analytic expression for the EIV slope parameter in terms of the variances (s d , s h ) and covariances (s h d ) of the measured variables, Or, Since s d = s h = 1 from Eq. (18). And the covariance (s h d ) is computed using the standard definition, Equation (26) can be solved if a priori knowledge of λ is available, which in turn requires a precise estimate of all sources of errors in the measured data. In the more practical case for unknown λ, we need to initiate a good starting estimate. Using since β W = 1 when λ = 1 from Eq. (26). Following Eq. (22), we therefore have our first estimate of β in the scattered plot of log 10 (1/τ ) and heights, this is just the geometric mean (GM) of the two OLS estimates, β d OLS and β h OLS , as can be seen by combining Eq. (12) and Eq. (13), The GM solution in Eq. (30) has the unique feature that this is the only case of EIV estimate which is both scale-invariant and symmetric in the variables (e.g., Ricker, 1984;Smith, 2009). While these properties do not necessarily imply that the 260 GM solution is the correct solution (discussed below), but this first estimate is essential to solve the problem of biasing in the estimated slope (Sect. 3.2.2). A specific example of GM fitting is presented in Fig. 1a.
The systematic offset between MR temperature and colocated lidar temperature as a result of using Eq. (30) is shown in Fig. 3a (red curve) and Fig. 3b (red histogram). The MR temperature is estimated with 1 day of radar data. Each of these temperatures are then compared with lidar temperatures at 88 km, 89 km and 90 km for the dates when lidar data are available.

265
The intrinsic noise in the lidar temperature is about 5-10 K (Reichert et al., 2019), which implies no temperature gradient is observed between 88-91 km in these data. Figure 3b (red histogram) reveals that the MR temperatures are over-estimated by a mean value of +63 K for the case of GM solution .
When compared to the temperature-gradient model derived from optical, satellite and rocket climatology (e.g., Holdsworth et al., 2006), it can be easily argued that our MSIS derived gradient model values (Fig. 2a) is more negative than expected.

270
If these values are shifted by a constant positive offset of +1 Kkm −1 , the absolute value of the estimated temperatures will increase by 10 K (Singer et al., 2004). This will further increase the offset between lidar and MR temperature, thereby shifting the histogram (in red) in Fig. 3b further to the right.
On the other hand, lidar temperatures are usually obtained during the night time, which can lead to a systematic offset due to day-night differences or tidal variations. As discussed by Hocking et al. (2004), the day-night temperature difference 275 at these altitudes is of the order of 3-4 K. This is significantly less than the standard errors in these temperatures, which is on average 19 K. Hence we expect the day-night difference in SGO's MR temperatures to be insignificant during the winter period. Moreover, any attempt to estimate MR temperatures using only night-time data has the adverse effect of reducing the accuracy in the estimated temperatures due to data loss. While no specific studies of tidal variation have been made for this location, the data from other sites (e.g., Hocking and Hocking, 2002;Stober et al., 2008) show that the temperature variation 280 due to tidal activity is typically less than 10 K. We can therefore rule out the possibility of an offset in the MSIS gradient model or tidal effects as the primary cause for the +63 K offset seen in Fig. 3 (red curve and histogram). Ricker (1984) emphasised that the biasing in GM solution is conditional upon the value of the correlation coefficient (r) between the variables, while Kimura (1992) demonstrated that this solution will be an overestimate in the case of low r value.
For the data set used in this work, we found that the correlation between log 10 (1/τ ) and height is typically 0.50 ± 0.05, 285 thereby indicating the presence of significant natural variation in the measurements. Furthermore, Jolicoeur (1990) used error modelling to conclude that r must be more than 0.6 for the the GM solution to be acceptable. In fact, we have observed that if we restrict our data selection process by excluding all data beyond the density contour 0.2 (Fig. 1b), this artificially increases the r to be typically around 0.66 ± 0.06 with the consequence of reduced biasing in the GM solution. However, temperatures estimated with any such arbitrary choice of data rejection criteria will lack consistency.

290
From purely mathematical perspective, Eq. (25) provides two equivalent explanations for the overestimate in GM solution.
Either the assumption of λ = 1 is not correct, and thus a higher value of λ is required to correct for the bias. This was validated in Sect. 3.1 using SCT calibration, where we have obtained λ OLS ef f = 1.70 ± 0.30 (Table 1). Or, the assumption of λ = 1 is correct, but the effective variance of the normalised heights, s h , is less than 1. This can be achieved by introducing a correction term in the statistical model of Eq. (22). In Sect. 3.2.2, we follow the second approach. And in Sect. 4, we have combined both 295 approaches to obtain a numerical estimate of the effective value of λ.

Reweighted second estimate of slope (β)
Equation (26) indicates that the overestimate of the GM solution can be corrected by using an increased value of λ, although it's an unknown parameter. Carroll and Ruppert (1996) attributed the parameter λ as purely due to measurement errors alone.
Following such argument, Carroll and Ruppert (1996) criticised Eq. (22) as an incomplete model for real data, since it implies 300 that in absence of measurement error only, "the data would fall exactly on a straight line". In practice, meteor echoes are subjected to random variations due to natural geophysical processes at all heights. Hence an additional term, known as equation error (µ ), is necessary to account for such geophysical or natural variation in the data. Therefore, a correction term (µ i ) is introduced in the statistical model (Carroll and Ruppert, 1996, Sect. 3), The adjusted value of the slope (β adj W ) in presence of equation error is always less than β W in the statistical model of Eq. (22). This can be seen by expressing Eq. (22) and Eq. (31) in terms of the variances, For the specific case of GM fitting or for λ = 1, Eq. (32) implies β adj W < 1. This solution requires that the observed variance of the measured heights in the data need to be weighted down by an appropriate scaling factor to compensate for the presence of 310 natural variation in the data. Alternatively, the effective value of λ in presence of equation error can be written as (Carroll and Ruppert, 1996), The quantity s µ essentially incorporates additional variability in the data that doesn't conform to the original construct or the characteristics being measured according to the underlying assumptions of Eq. (6). For example, meteor trails can get 315 modified by wind effects, ion composition, meteor fragmentation, strong ionospheric currents as well as temperature and pressure fluctuations on various spatial and temporal scale (Hocking, 2004;Younger et al., 2014). Hence the practical effect of the extra term in Eq. (31) is to increase the effective variance of h i by an additional factor of s µ (Gillard and Iles, 2005).
Assessing the sampling distribution of s µ would require carefully designed replicates of observations as well as long term comparison of MR temperatures with other colocated instruments. This line of reasoning limits the practical use of Eq. (33).

320
Because an absolute knowledge of measurement error or various sources of natural variation is not available for this data on a day-to-day basis.
Smith (2009)  In contrast to Carroll and Ruppert (1996), we therefore no longer stringently distinguish between different types of errors in the data (e.g., Sprent, 1990, p 13). In other words, the error variances, s δ and s ε , now consist of both measurement errors 330 and partly the natural geophysical variation. Hence the equation error variance, s µ in Eq. (33), is simply a representative of the asymmetric component of the natural variation in the data. In MR data, such asymmetric effect of natural variation occurs mainly at altitudes approximately above 95 km due to geomagnetic effect, and approximately below 85 km due to electron recombination and attachment (Sect. 1 and the references therein). To compensate for this asymmetric component of natural variation in the data, the observed variance of the normalised heights (h i ) must be replaced by its reweighted value (h * i ) as, such that the scaling factor, ν, is less than 1. In the Appendix section we have derived an analytic solution of ν in terms of s µ : For λ = 1, Eq. (25) can be expressed in terms of the reweighted normalised heights as , Using Eq. (34) in Eq. (36), the bias-corrected slope is, where ν can be estimated from Eq. (35). Equation (37) is therefore the reweighted second estimate of the normalised slope.
It is to be noted that λ has been essentially eliminated in Eq. (37), and hence an absolute knowledge of the ratio of error variances is no longer necessary as a result of data normalisation. Equation (35)  (37) to obtain the normalised slope coefficient (β adj W ). In the final step, the slope and the intercept in the original coordinate system is calculated from, And, This value of β and the temperature gradient (Fig. 2)

Bootstrap analysis of the standard error in β
Monte Carlo simulation is carried out to assess the standard error in the estimated β (Press et al., 1992, Sect. 15.6). The actual daily data set, with its N data points, are used to generate 20000 synthetic data sets using random resampling. Each of these synthetic data set is made by randomly drawing N data points with replacement from the original data set, while allowing repeated draw of the same data. Hence for large N , approximately e −1 th fraction of the actual measurements are replaced by 365 duplicated original points. As mentioned above we have repeated the resampling procedure 20000 times and then assessed the standard error (SE) of the parameters, V ar(β adj W h i ) and V ar(d i ), from the 95% confidence interval (CI) of the distribution as, Equation (38) can be used to compute the standard error in β from the standard error in β adj W √ s h , √ s d and the algebra of error 370 propagation. Figure 5 presents the error estimates in β for the data sets used in this paper. As expected, the estimated standard error was found to be inversely related to √ N . A linear fit revealed the following relation, where ∆β represents the standard error as a percentage of the estimated β.

Results and Discussion
We have applied the method formulated in Sect. 3.2.2 to the MR data from the period October 2015 to March 2016. This is presented in Fig. 6, along with the data from colocated lidar observation. The percentage differences between MR temperatures 380 and lidar are shown in Fig. 7a for altitudes near peak meteor counts. The mean difference between lidar and MR temperatures are less than ±1 K (or about 0.2%) which is negligible within error consideration. This means we have achieved obtaining an unbiased temperature estimation. The root-mean-square differences (RMS) is about 7%. This is expected since a 0.7 Kkm −1 error in MSIS gradient model will introduce a 7-10 K variation in the data for a typical winter temperature of 200 K in the peak meteor region. Meek et al., 2013). In addition, lidar temperatures have intrinsic errors of 5-10 K. The 385 10 percentile and the 90 percentile value of daily number of meteor detections in this data set is 1815 and 4057 respectively.
Using Eq. (41), this corresponds to an error in the estimated slope to be in the range of 5-7%. Hence the RMS difference of 7% between MR temperature and lidar is not statistically significant within the error limits. In other words, we have managed to estimate MR temperatures correctly without requiring any external calibration.
The standard error in temperatures estimated using EIV analysis is on average 15 K. As can be seen in Fig. 6, around 390 middle of Feb 2016 there is an abrupt change in temperature that is much larger than 2-sigma error bars. This temperature drop corresponds to a minor sudden stratospheric warming (SSW). Then, it was a "major final warming" starting on 5-6 Mar 2016 (e.g., Manney and Lawrence, 2016), which also coincides with a temperature decrease seen in this data. This demonstrates that there are day-to-day geophysical processes that can be revealed by the method developed in this work.
For direct comparison with the results above, we have also estimated the MR temperatures using the revised SCT calibration 395 procedure described in Sect. 3.1. For this we have estimated the OLS slope (β d OLS ) and used √ s δ = 0.18 (Table 1)  between the MR temperatures and that of lidar is again less than 2%, thereby showing that the biasing effect has been properly 400 corrected by this calibration procedure. However, the RMS difference has increased to 10%, thereby clearly indicating the presence of artefacts in these calibrated MR temperatures.
As emphasised in Sect. 3.2, the statistical models in Eq. (22)  slope estimated showed overcorrection and thus implying a higher value of λ is required. Or equivalently, we have continued 405 to assume λ = 1, and instead changed the statistical model to Eq. (31). This model effectively reduced the variance of the measured heights to allow a bias-corrected GM fitting. This solution is given by Eq. (37), which has been validated by the excellent agreement between lidar and MR temperatures (Fig. 6). Alternatively, we can show numerically that Eq. (37) is equivalent to Eq. (26) as follows: First, the parameter λ is allowed to run as free variable from 0.5 to 2.5 in steps of 0.01. Each value of λ is then used in Eq.

410
(26) to obtain a value of β W . Second, the same data set is now used to obtain the bias-corrected slope using Eq. (37) for which we get β adj W = 0.75. The slope estimates from these two independent procedures are shown in Fig. 8 for one specific date. The value of β adj W = 0.75 corresponds to λ EIV ef f = 1.7 when λ is introduced as free variable in Eq. (26). Hence, the solution given by Eq. (37) for λ = 1 is essentially equivalent to increasing the effective value of λ in Eq. (26).
The effective value of λ (or, λ EIV ef f ) for each of the 24 h data set can estimated numerically by equating Eq. (26) with Eq.

415
(37), The distribution of λ EIV ef f obtained from the full data set shows a mean and standard deviation of 1.76 and 0.06 respectively ( Fig. 4d and Table 1). This is well within the limit of the experimental value of λ OLS ef f =1.7±0.3 which we have obtained using the SCT method ( Fig. 4c and Table 1). Another implication of Eq. (37) can be seen in Fig. 3 (red line and the red histogram).

420
The temperature estimated using the GM fitting (the first estimate) shows a systematic positive offset of +63±16 K. For a typical winter temperature of 200 K in the meteor regions Hall et al., 2006), this would correspond to a positive biasing of 23-40%. In the example provided for the sample data (Fig. 8), we see that use of Eq. (37) reduced the GM slope from 1 to 0.75 or by 25%. Thus, application of Eq. (37) has effectively reduced the offset by shifting the red histogram in Fig. 3b leftwards without requiring any external calibration. This is presented in Fig. 7a where the mean offset is nearly zero.
Despite the excellent agreement in the statistical measure of ratios of error variances using two independent methods, λ OLS ef f and λ EIV ef f , the estimated values of log 10 (1/τ ) and (height) are clearly higher than what has been reported in literature previously (Table 1). Hocking (2004) estimated (height) = 3.25 km from pure measurement error perspective, using a numerical model for a pulse length equivalent to 2 km, and a meteor at an altitude of 90 km and zenith angle 50 • . Likewise, their estimation of log 10 (1/τ )=0.14 was based on simulation studies for meteors below 95 km and confirmed that decay times 430 variability arises due to 27% variability in K amb and 8% variability in temperatures over meteor region. The argument for choosing 95 km as the maximum height is that above this altitude meteor decay rates are substantially affected by processes other than ambipolar diffusion. Holdsworth et al. (2006) and Kim et al. (2012) also restricted their data selection process below the altitude of 95 km (or within 86-96 km in latter), and concluded that (height) ≈ 1.0 km is the best value to calibrate their MR temperatures with optical, satellite or rocket climatology using the SCT correction. In addition, Holdsworth et al. (2006, 435 p 5) used various other data rejection criteria and assumed log 10 (1/τ ) = 0.14 for calibrating the slope. In principle, these literature values of the error variances ( Table 1) were not derived following any statistical procedure, but rather assumed on the basis of speculation and convenience.
Our estimates of log 10 (1/τ ) and (height) are higher, because we chose to include meteors at all altitudes and all decay times in the data selection criteria. Due to geomagnetic effect the decay times are higher than expected from the simple 440 ambipolar case at altitudes roughly above 95 km, whereas below 85 km the effect is the opposite (Sect. 1 and the references therein). More than 30% of the meteor echoes are from above 95 km and below 85 km. This means the value of s δ and s ε will be higher due to increased geophysical variation in the data if no restriction is applied to the height range and decay times in the data selection process. By including meteor echoes from all altitudes we have avoided sampling bias, and hence we believe our estimates of s δ and s ε incorporate all sources of variation in the data, both instrumental and geophysical, thereby 445 probing the true meteor population. Moreover, the altitude range where most meteor decays deviate from the assumption of ambipolar diffusion is a dynamical process, showing seasonal variation and can be location-dependent. In essence, there is no way to determine what are outlier in the data, since no error estimates are provided in the measured heights and decay times by the radar software. By keeping meteor echoes from all altitudes, we have retained most of the natural variation in the measurements, and this ensures that the robustness of the proposed EIV method can be assessed to the fullest extent. Moreover, 450 any arbitrary choice of data rejection criteria will also change the slope coefficient in an arbitrary way, and thus the estimated MR temperatures will lack consistency. To address this issue, in Sect. 3.2.2, we have incorporated the statistical effects of natural geophysical processes at all altitudes in the regression model. Due to lack of colocated lidar measurements at other times of the year we have restricted our temperature estimation only for the winter times. Complications might arise during certain days in the summer time when the turnaround of the diffusion 455 coefficient profile occurs at higher altitudes close to the peak meteor regions. Hence the decay of meteor radar echos deviate significantly from diffusion-only evolution (Lee et al., 2013). This would require an improvement in the original physical model in Eq. (6), or some other alternative method of temperature estimation technique as discussed below.
In addition to the temperature-gradient technique used in this work, another common method for estimating MR temperatures is to directly use some model pressures in Eq. (2). This method doesn't require any line-fitting analysis, and in theory, can be used to determine temperatures at any heights within the meteor height distribution. The main drawback of this technique is that the pressure is a sensitive function of heights. Hence, even a small uncertainty in these model pressure values can lead to a large bias in the estimated temperatures thereby requiring further recalibration. Such biasing effect can be seen in Meek et al. (2013Meek et al. ( , p 1273. Using the pressure values from two different models, Holdsworth et al. (2006)  for calibrating the MR temperatures against optical data. Later on, Dyrland et al. (2010) used the same data set as Hall et al. (2006), but calibrated their MR temperatures using satellite data and consequently derived a different calibration equation.
As an alternative to pressure and temperature-gradient methods, Hocking (2004) claimed that the correlation coefficient (r) between log 10 (1/τ ) and height is linearly related to the MR temperatures, and obtained T 0 = 360r − 42 for T 0 < 190 K.
However, this method is based on the flawed assumption that a large value of correlation coefficient between r and T 0 implies 470 a linear relation between these two variables. Any such linear relationship between r and T 0 must be validated with physics theory, which was not provided in Hocking (2004).
More recently, Lee et al. (2016) provided a good alternative technique for estimating MR temperatures near peak heights.
Using both theory and experimental data, Lee et al. (2016) demonstrated a linear relation between the MR temperatures and the width of the meteor height distribution. Once the constant of proportionality of this linear correspondence is found, the 475 MR temperatures can be directly estimated from the meteor height distribution. Lee et al. (2018) derived this proportionality constant empirically by using temperatures from satellite data, which in turn, is essentially a form of calibration process.
In all generality, it is desirable to avoid any kind of calibration process and instead, formulate an independent and unbiased estimate of temperatures using MR data alone. The statistical method presented in this paper pave the way towards achieving that goal.

Summary
The biasing effect in MR temperature has been a pressing issue for last two decades. Attempts have been made in the past to correct the slope in the scattered plot of log 10 (1/τ ) and height, usually either by direct calibration with optical or satellite data or by an arbitrary choice of data rejection criteria to exclude parts of measurements. This paper has addressed the underlying reasons for such biasing effect, which is mainly due to the presence of various error terms in the physical model. We have 485 reviewed the conventional calibration procedure (1), and then provided an alternative but more robust method (2) for estimating MR temperature that doesn't require any calibration. We have applied both of these methods to the MR data from winter 2015-2016, and assessed the quality of the estimated MR temperature using colocated lidar measurements. The key points from each of these two aspects of the paper are given below.
1. This paper has reviewed the statistical comparison technique (SCT), originally proposed by Hocking et al. (2001b), 490 within the context of MR temperature calibration. We have extended the theoretical basis of the SCT method to obtain an estimate of error variances of log 10 (1/τ ) and height using colocated lidar measurements. No significant offset was seen in the calibrated MR temperature, even without applying any outlier rejection criteria. But artefacts introduced due to the difference in measurement techniques between MR and lidar was clearly visible in the estimated temperatures in the form of increased RMS differences. 495 2. As an alternative method, we have applied the Errors-in-Variables (EIV) regression model to estimate the slope in the scattered plot of log 10 (1/τ ) and height. This model takes into account the measurement errors in both abscissa and ordinate. Following Carroll and Ruppert (1996), an important source of variation, known as equation error was introduced in the EIV method. The term equation error was shown to play a critical role in determining the MR temperature in the most unbiased way. The main cause of this error was identified as due to a number of geophysical effects in the mea-500 sured data that was not taken into account in the original physical model. We have used the theoretical considerations by Carroll and Ruppert (1996) and Smith (2009) to formulate a closed-form analytic solution of the unbiased slope coefficient. This solution allows an independent estimate of atmospheric temperatures using meteor radar at the altitude of peak meteor occurrences. The temperatures estimated using this method show very good agreement with colocated lidar measurement, with 7% or better accuracy and without any systematic offset. The quality of the estimated temperatures 505 using the EIV method was significantly better than the conventional SCT calibration method.
Data availability. The meteor radar data were collected at SGO (https://www.sgo.fi/Projects/SLICE/). Data used in the paper are available at https://www.sgo.fi/pub/AMT_2020_ESarkar/. The MSIS model data were obtained from https://ccmc.gsfc.nasa.gov/modelweb/. The lidar data are available at https://halo-db.pa.op.dlr.de/mission/109. The revised statistical model of Eq. (31) implies that the variable h i is now redefined in its adjusted form (h * i ), such that the expectation value of µ , or < µ > = 0, to account for the biasing effect due to asymmetric natural variation in the measured normalised heights. The variance of the reweighted heights are related to the original measurements as, By defining a parameter, ν, such that, Equation (A2) can be expressed as, V ar(h * i ) = V ar(νh i ) = ν 2 V ar(h i ) = ν 2 (A4) (A1) simplifies to, In essence, the central trend line defined by the revised statistical model of Eq. (31) for the reweighted normalised heights (h * i ) is equivalent to GM fitting corresponding to Eq. (22) for the original normalised heights (h i ). The χ 2 -minimisation for the two models are related as, where β W = 1 for λ = 1. Using Eq. (A1) and Eq. (A5), the merit function inside the parenthesis of the left-side of Eq. (A6) approximates to, where again we have assumed that all the error terms are mutually independent and uncorrelated to the measured variables.