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 ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ 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 ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$. 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 co-located lidar measurements, the systematic offset in the estimated temperature is shown to have reduced to 5 % or better on average.

1 Introduction and background

As meteoroids enter the Earth's atmosphere, they produce ionised 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¹⁴ electrons m⁻¹, the trail is called “underdense”, meaning each electron in the trail scatters independently (e.g. Bronshten, 1983, p. 356). 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 amplitude of the radio signal from an underdense trail decays with time (t) as

$$\begin{array}{}\text{(1)}& A\left(t\right)=A\left(\mathrm{0}\right)e^{-\mathrm{16}{\mathit{\pi }}^{\mathrm{2}}{D}_{\mathrm{a}}t/{\mathit{\lambda }}_{\mathrm{r}}^{\mathrm{2}}},\end{array}$$

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

$$\begin{array}{}\text{(2)}& D_{\mathrm{a}}=K_{\mathrm{amb}}{\displaystyle \frac{T^{\mathrm{2}}}{P}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{r}}^{\mathrm{2}}\ln \mathrm{2}}{\mathrm{16}{\mathit{\pi }}^{\mathrm{2}}\mathit{\tau }}}.\end{array}$$

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

$$\begin{array}{}\text{(3)}& P\left(h\right)=P\left(\mathrm{0}\right)e^{-{\int }_{\mathrm{0}}^{\mathrm{h}}\frac{mg}{kT\left(z\right)}\mathrm{d}z},\end{array}$$

where m is the mass of a typical atmospheric molecule, g is the acceleration of gravity, k is the Boltzmann constant and z is an axis along the vertical. Substituting the equation for pressure in Eq. (2) and differentiating Eq. (4) provides the height profile of the decay time:

$$\begin{array}{}\text{(4)}& {\log }_{\mathrm{10}}{D}_{\mathrm{a}}\left(h\right)=\mathrm{2}{\log }_{\mathrm{10}}T\left(h\right)+{\log }_{\mathrm{10}}e{\displaystyle \frac{mg}{k}}\underset{\mathrm{0}}{\overset{\mathrm{h}}{\int }}{\displaystyle \frac{\mathrm{1}}{T\left(z\right)}}\mathrm{d}z+\mathrm{\Psi },\end{array}$$

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}h}}{\log }_{\mathrm{10}}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\tau }\left(h\right)}}\right)=\mathrm{2}{\log }_{\mathrm{10}}e{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}h}}{\displaystyle \frac{\mathrm{1}}{T\left(h\right)}}+{\log }_{\mathrm{10}}e{\displaystyle \frac{mg}{k}}{\displaystyle \frac{\mathrm{1}}{T\left(h\right)}},\end{array}$$

where Ψ is a constant. Equation (5) states that the height profile of decay time is a function of both temperature and temperature 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₀) at the peak height by fitting a linear function (Hocking, 1999). A linear approximation of Eq. (5) is

$$\begin{array}{}\text{(6)}& T_{\mathrm{0}}=\mathit{\beta }(\mathrm{2}<{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}h}}>+{\displaystyle \frac{mg}{k}}){\log }_{\mathrm{10}}e,\end{array}$$

where β is the slope of the scattered plot of ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height, and T₀ is the average temperature of the atmosphere at the height of peak meteor occurrence. A typical scattered plot of height and ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ shows significant variation in the measured data along both abscissa and ordinate (Fig. 1). Traditionally, the slope (β) is estimated using the ordinary least-squares (OLS) method with ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ as the independent variable. The justification of using ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ as an independent variable is that the measurement errors in τ are smaller than those in heights (Hocking et al., 1997). While the pulse length and 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., 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).

Figure 1 (a) Typical scatter plot of ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height. The lines correspond to best-fit models using different regression methods described in the text. The green and blue line corresponds to the ordinary least-squares method (OLS), with ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height as independent variables respectively. The red line corresponds to the geometric mean (GM) of ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$ and ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}$. (b) The bivariate distribution of the data. The measured height and ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ are converted to dimension-free coordinates using Eq. (18). The relative density contours are obtained by counting the number of detections in a circle of unit area relative to the detection density at the height of peak meteor occurrences at the centre.

Temperature estimation from meteor radar (MR) data requires obtaining the best-fit regression line in the scattered plot of ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height. However, the pioneer work done by Hocking (1999) to implement this method using ordinary least-squares fitting showed a clear systematic offset between the MR temperature and co-located 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

$$\begin{array}{}\text{(7)}& {\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}=(\mathrm{1}-{\displaystyle \frac{s_{\mathit{\delta }}}{s_{\mathrm{d}}}}){\mathit{\beta }}_{\mathrm{SCT}}\end{array}$$

or

$$\begin{array}{}\text{(8)}& {\mathit{\beta }}_{\mathrm{SCT}}=(\mathrm{1}-{\displaystyle \frac{s_{\mathit{\epsilon }}}{s_{\mathrm{h}}}}){\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}},\end{array}$$

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 }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and of height respectively and ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$ and ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}$ are slope coefficients when diffusion coefficient or height is treated as an independent variable respectively in the OLS regression analysis. The calibrated slope, β_SCT, is traditionally 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 a calibration procedure is that the calibration factors (s_δ and s_ε) in the parentheses 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 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 an independent variable, and both variables are subject to intrinsic measurement errors and various geophysical effects. The reasons for such a bias, and thus the need for calibration, are discussed on theoretical and experimental 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. This analysis does not require an absolute knowledge of s_δ or s_ε, but only the relative value is needed. As an alternative to SCT calibration, in Sect. 4, an independent slope estimation is obtained using the errors-in-variables model. A comparison study of the estimated MR temperatures with co-located lidar temperatures is discussed to validate the method.

2 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 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 for 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 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 non-meteor targets during substorms associated with ionospheric plasma waves generated due to Farley–Buneman instability (Kelley, 2009). The Doppler velocity of such echoes can be more than a few hundreds of metres per second, which are mostly detected at low elevation (Lukianova et al., 2018). For the Sodankylä radar, Kozlovsky and Lester (2015) identified ground echoes modulated by the ionosphere during pulsating auroras. These targets have near-zero Doppler velocities and are also observed at low elevation. Furthermore, the SKiYMET system detects both underdense and overdense echoes as valid meteors. However, more than 95 % of detections are underdense (Hocking et al., 2001b). The percentage of overdense trails may be larger during some meteor trails such as Geminids or Quadrantids, and this leads to underestimation of temperatures (Kozlovsky et al., 2016). This artefact in temperature was found to be reduced for Sodankylä radar for a zenith angle of less than 50^∘. The initial data selection criteria is kept to a bare minimum, such that all heights and decay times are included, as long as they are unambiguous detections above a 40^∘ elevation angle with Doppler radial velocity in the range ±100 m s⁻¹. Subsequently, to improve the temperature estimation, we used a contour selection process; i.e. the data outside a certain contour in the normalised height–${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ distribution were rejected (Fig. 1b).

For temperature estimation we considered daily data for a 6-month period from October 2015 to March 2016 since simultaneous lidar measurements were available during this time. The Compact Rayleigh Autonomous Lidar (CORAL) provided vertical profiles of the atmospheric temperature at 27–98 km over Sodankylä as part of the GW-LCYCLE-II (Gravity Wave Life Cycle Experiment) campaign in winter 2015/2016 (Reichert et al., 2019). The median number of daily meteor detections in this data set is 1652, with a minimum of 410 meteor detections on 4 October 2015 and a maximum of 3456 detections during the Geminids meteor shower on 14 December 2015.

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 and 95 km. A third-degree polynomial fit is carried out to obtain the height profiles at 06:00, 12:00, 18:00 and 24:00 UT. For each time interval, the gradient at the respective meteor peak height, as well as at 1 km above and below the peak height, is estimated. 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 gradient model values are shown in Fig. 2. The standard deviation of the MSIS model values corresponds to roughly of the order of 0.7 K km⁻¹.

Figure 2 (a) Temperature gradient model derived from MSIS90. (b) Peak meteor heights for the data used in this work and (c) the daily meteor detection for zenith angle less than 50^∘ and velocity in the range ±100 m s⁻¹.

Our reasons for choosing the MSIS temperature gradient are twofold. Firstly, MSIS data are easily accessible from the online 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 is on average 10 % (Hocking, 1999) or less. 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.

3 Method: regression analysis

3.1 Estimation of error variances in decay time and height

In the following text we use the notation and formulation in Gillard and Iles (2005). The observables, ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ 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 ith meteor detection. For consistency we also assume that ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ is presented in abscissa, and height 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

$$\begin{array}{}\text{(9)}& {\mathit{\eta }}_i=\mathit{\alpha }+\mathit{\beta }{\mathit{\xi }}_i,i=\mathrm{1},\mathrm{2},\mathrm{\dots },N.\end{array}$$

Due to measurement errors and various geophysical processes, the true values ξ_i and η_i will be subject to random errors, and hence the observable d_i and h_i will have scatter around the linear model in Eq. (9):

$$\begin{array}{}\text{(10)}& d_i={\mathit{\xi }}_i+{\mathit{\delta }}_i\end{array}$$

and

$$\begin{array}{}\text{(11)}& h_i={\mathit{\eta }}_i+{\mathit{\epsilon }}_i=\mathit{\alpha }+\mathit{\beta }{\mathit{\xi }}_i+{\mathit{\epsilon }}_i,\end{array}$$

where δ_i and ε_i are errors in the measured ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height respectively and are assumed to be mutually uncorrelated, have zero mean and be 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 (Var) as (e.g. Smith, 2009; Keles, 2018)

$$\begin{array}{}\text{(12)}& {\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}={\displaystyle \frac{\text{Cov}(d_i,h_i)}{\text{Var}\left(d_i\right)}}\equiv r{\displaystyle \frac{\sqrt{s_{\mathrm{h}}}}{\sqrt{s_{\mathrm{d}}}}},\end{array}$$

where ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$ is the OLS slope estimated by considering ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ as an independent variable, and r is the Pearson product-moment correlation coefficient between d and h. Likewise, by reversing the arguments above, it is trivial to show that the reciprocal value of the OLS slope estimate with height as an independent variable is (e.g. Smith, 2009)

$$\begin{array}{}\text{(13)}& {\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}={\displaystyle \frac{\text{Var}\left(h_i\right)}{\text{Cov}(d_i,h_i)}}\equiv {\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\sqrt{s_{\mathrm{h}}}}{\sqrt{s_{\mathrm{d}}}}}.\end{array}$$

To see the effect of errors in the independent variable on the OLS slope estimate, e.g. for ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$, Eqs. (10) and (11) are used in Eq. (12):

$$\begin{array}{}\text{(14)}& {\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}={\displaystyle \frac{\text{Cov}({\mathit{\xi }}_i+{\mathit{\delta }}_i,\mathit{\alpha }+\mathit{\beta }{\mathit{\xi }}_i+{\mathit{\epsilon }}_i)}{\text{Var}({\mathit{\xi }}_i+{\mathit{\delta }}_i)}}.\end{array}$$

Since δ_i, ε_i and ξ_i are mutually independent, Eq. (14) simplifies to

$$\begin{array}{}\text{(15)}& {\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}={\displaystyle \frac{\mathit{\beta }\text{Cov}({\mathit{\xi }}_i,{\mathit{\xi }}_i)}{\text{Var}({\mathit{\xi }}_i+{\mathit{\delta }}_i)}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{\text{Var}\left({\mathit{\delta }}_i\right)}{\text{Var}\left({\mathit{\xi }}_i\right)}}}\mathit{\beta }=\mathit{\zeta }\mathit{\beta },\end{array}$$

where ζ is known as the attenuation or regression dilution bias. Since variances are always positive by definition, Eq. (15) shows that in the presence of measurement error in the so-called independent variable (in abscissa), the OLS slope estimate (${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$) will always be smaller than the unbiased slope β. Likewise, ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}$ is greater than β if there is error in the measured height (specific example presented in Fig. 1a). By substituting $d_i={\mathit{\xi }}_i+{\mathit{\delta }}_i$, Eq. (15) can be rearranged as

$$\begin{array}{}\text{(16)}& {\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}=(\mathrm{1}-{\displaystyle \frac{\text{Var}\left({\mathit{\delta }}_i\right)}{\text{Var}\left(d_i\right)}})\mathit{\beta }.\end{array}$$

Equation (16) is a well known identity in statistical literature (e.g. Carroll and Ruppert, 1996; Frost and Thompson, 2000), 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. A common practice is to calibrate the biased slope with optical or satellite data by arbitrarily choosing a value of s_δ or s_ε (e.g. Holdsworth et al., 2006; Hocking et al., 2007; Kim et al., 2012). In the remaining part of this section, we demonstrate how to obtain an average value of the error variances in these data following a revised calibration procedure.

For each 24 h of the data set, we performed two OLS fittings to estimate ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$ and ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}$ by using Eqs. (12) and (13). The corresponding biased temperatures, $T_{\mathrm{MR}}^{\mathrm{d}}$ and $T_{\mathrm{MR}}^{\mathrm{h}}$ 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 co-located lidar temperatures (T_lidar) as the reference values. Using Eqs. (7) and (8), and noting that the slope is proportional to the estimated temperatures from Eq. (6), we obtain

$$\begin{array}{}\text{(17)}& s_{\mathit{\delta }}\approx \left({\displaystyle \frac{T_{\mathrm{lidar}}-T_{\mathrm{MR}}^{\mathrm{d}}}{T_{\mathrm{lidar}}}}\right)s_{\mathrm{d}}\text{and}{s}_{\mathit{\epsilon }}\approx \left({\displaystyle \frac{T_{\mathrm{MR}}^{\mathrm{h}}-T_{\mathrm{lidar}}}{T_{\mathrm{MR}}^{\mathrm{h}}}}\right)s_{\mathrm{h}},\end{array}$$

where $T_{\mathrm{MR}}^{\mathrm{d}}$ and $T_{\mathrm{MR}}^{\mathrm{h}}$ are MR temperatures estimated using OLS fitting with ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height as an independent variable respectively. Furthermore, if the measurements are normalised with the mean and standard deviation (SD) as

$$\begin{array}{}\text{(18)}& d_i{}^{\prime }={\displaystyle \frac{d_i-\text{mean}\left(d_i\right)}{\sqrt{s_{\mathrm{d}}}}}\text{and}{h}_i{}^{\prime }={\displaystyle \frac{h_i-\text{mean}\left(h_i\right)}{\sqrt{s_{\mathrm{h}}}}},\end{array}$$

then $\text{Var}\left(d_i{}^{\prime }\right)=\text{Var}\left(h_i{}^{\prime }\right)=\mathrm{1}$, and the OLS estimate of the ratio of the measurement error variance is (from Eqs. 17 and 18)

$$\begin{array}{}\text{(19)}& {\mathit{\lambda }}_{\mathrm{eff}}^{\mathrm{OLS}}={\displaystyle \frac{s_{\mathit{\epsilon }{}^{\prime }}}{s_{\mathit{\delta }{}^{\prime }}}}\approx {\displaystyle \frac{T_{\mathrm{MR}}^{\mathrm{h}}-T_{\mathrm{lidar}}}{T_{\mathrm{lidar}}-T_{\mathrm{MR}}^{\mathrm{d}}}}\cdot {\displaystyle \frac{T_{\mathrm{lidar}}}{T_{\mathrm{MR}}^{\mathrm{h}}}},\end{array}$$

where the error variances, $s_{\mathit{\epsilon }{}^{\prime }}$ and $s_{\mathit{\delta }{}^{\prime }}$, are in the dimension-free system defined by Eq. (18). In essence, ${\mathit{\lambda }}_{\mathrm{eff}}^{\mathrm{OLS}}$ is a measure of all sources of errors in the normalised heights and decay times that cause the real data to deviate from the idealised physical model of Eq. (6), thereby producing a typical scatter as seen in Fig. 1.

Figure 3 (a) Temperature estimated in OLS method using ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ (green) and height (blue) as an independent variable. Also shown are (red) the temperatures obtained using the geometric mean (GM) fitting and the lidar temperatures (in black). (b) The offset between the lidar (T_lidar) temperatures and the estimated MR temperatures (T_MR) using OLS fitting and GM fitting (without contour selection).

${\mathit{\lambda }}_{\mathrm{eff}}^{\mathrm{OLS}}$, s_δ and s_ε were estimated by Eqs. (17) and (19) using 24 h of MR data and co-located lidar temperatures at 88, 89 and 90 km for the dates for which lidar data were available. The biasing effect on the OLS estimate of MR temperatures with ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height as an 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 −75 and +335 K occurs depending on whether ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ or height respectively is considered as an independent variable. In practice, the magnitude of these biases is related to the total errors in heights and ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ from Eq. (17), which was not taken into account by the OLS regression model in Eqs. (12) and (13). Furthermore, the relative density contour lines for each data point is obtained by counting the number of detections in a circle of unit area in the normalised height–${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ plane relative to the number of detections in unit area at the peak meteor region (Fig. 1b). In addition, all the error estimates from the SCT method are obtained at contour levels 0, 0.2 and 0.4, and the results are presented in Table 1. Several key features of these data are reflected in Table 1. Despite the data transformation via Eq. (18), the average error variances for the height data are more than those for the ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ data in these coordinates. This asymmetry in the error variances implies that the bivariate distribution is slightly skewed away from a perfectly normal distribution along the y direction. However, this effect of asymmetric error variances is less pronounced near the core of this distribution. This implies that the parameter ${\mathit{\lambda }}_{\mathrm{eff}}^{\mathrm{OLS}}$ is closer to 1 near the core of this distribution and further away from 1 near the tail of the distribution. The data near the outer contour area are subject to larger parameter estimation error due to observational limitation. For example, the zenith-angle-dependent error in height can easily skew the distribution in this direction. On the other hand, natural geophysical variability in these data contributes significantly at all altitudes. As seen in Table 1, the average error in height for contour levels 0, 0.2 and 0.4 are 5.0, 3.0 and 2.2 km respectively. This error is significantly higher than what is expected from purely parameter estimation error (for zenith angle less than 50^∘, the error in height is 1 km or less), indicating that geophysical variability dominates the total error variance in these data at all altitudes. The average error variances estimated at the contour level 0.2 (Table 1) correlate very well with the values reported in Hocking (2004) and Holdsworth et al. (2006). Hocking (2004) estimated △(height) =3.25 km 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 }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })=\mathrm{0.14}$ was based on simulation studies for meteors below 95 km and confirmed that “decay times' variability” arise due to 27 % variability in K_amb and 8 % variability in temperatures over the 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, p. 5) applied similar data rejection criteria and found out that $△{\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })=\mathrm{0.14}$ is required to calibrate the slope if ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ is used as an independent variable in the OLS regression model. Likewise, Thorsen et al. (1997) performed the comparison between the parameter estimation error and the geophysical variability for estimating the mean wind field in the middle atmosphere and found that the geophysical variability dominated at all heights.

It is worth noting that instead of directly using the individual observation between biased MR temperatures and lidar measurements from Eq. (17), we have used the statistical mean of differences for calibration. This is because lidar data are 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 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 are 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, the 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 kilometres. On the other hand, the lidar may resolve these gravity waves if the runtime is shorter than the period of these waves. As gravity wave amplitudes can be up to 10–15 K at these altitudes (Reichert et al., 2019), we cannot expect perfect agreement between radar and lidar temperature due to geophysical variation which shows 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 due to natural geophysical processes will persistently introduce artefacts in the estimated 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 for the given location. In all generality, it is desirable to avoid any kind of calibration process and instead formulate an independent estimate of temperatures using MR data alone. As an alternative to the OLS method, errors-in-variables (EIV) regression analysis provides a way to incorporate the error variances in both height and ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ data, thereby reducing the biasing effect in the estimated slope parameter.

Table 1 The average value of the (square root of) error variances in height and ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$; $s_{\mathit{\epsilon }{}^{\prime }}$ and $s_{\mathit{\delta }{}^{\prime }}$ are given along with the average value of λ from SCT calibration (${\mathit{\lambda }}_{\mathrm{eff}}^{\mathrm{OLS}}$) at contour levels 0, 0.2 and 0.4 for winter 2015–2016.

3.2 Errors-in-variables (EIV) model: GM solution

For fitting a straight line model, such as

$$\begin{array}{}\text{(20)}& y\left(x\right)=a+bx\end{array}$$

to a set of N data points (x_i, y_i) measured with errors, the corresponding χ² merit function is (Press et al., 1992, p. 660)

$$\begin{array}{}\text{(21)}& {\mathit{\chi }}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^N{\displaystyle \frac{(y_i-a-b{x}_i{)}^{\mathrm{2}}}{{\mathit{\sigma }}_{yi}^{\mathrm{2}}+b^{\mathrm{2}}{\mathit{\sigma }}_{xi}^{\mathrm{2}}}},\end{array}$$

where σ_yi and σ_xi are the standard deviation of the ith data point, and the weighted sum in the denominator of Eq. (21) can be interpreted as the weighted error of the ith data point. The regression coefficients, a and b, can be found by minimising the merit function with respect to these coefficients following any suitable numerical root-finding routine. However, under the assumption of symmetric error variances, it is possible to derive an analytic solution for the regression coefficients. This solution, when the data are appropriately normalised, leads to a slope estimate which is both scale-invariant and symmetric with respect to the data.

Application of Eq. (21) to physics data requires that all measured variables are dimensionally consistent so that χ² is dimension-free. Moreover, the analysis in this section requires that the measurements are presented in an appropriate dimension-free system. This facilitates the direct comparison between different parameters, such as the measured variables or the associated error variances. By applying the coordinate transformation introduced in Eqs. (18) to (9), we therefore intend to solve the simplified bivariate linear system of equations,

$$\begin{array}{}\text{(22)}& {\mathit{\eta }}_i{}^{\prime }=\mathit{\beta }{\displaystyle \frac{\sqrt{s_{\mathrm{d}}}}{\sqrt{s_{\mathrm{h}}}}}{\mathit{\xi }}_i{}^{\prime }\equiv {\mathit{\beta }}_{\mathrm{W}}{\mathit{\xi }}_i{}^{\prime },\end{array}$$

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)

$$\begin{array}{}\text{(23)}& {\mathit{\chi }}^{\mathrm{2}}({\mathit{\beta }}_{\mathrm{W}},{\mathit{\alpha }}_{\mathrm{W}}=\mathrm{0})\equiv \sum _{i=\mathrm{1}}^N{\displaystyle \frac{(h_i{}^{\prime }-{\mathit{\beta }}_{\mathrm{W}}{d}_i{}^{\prime }{)}^{\mathrm{2}}}{s_{\mathit{\delta }{}^{\prime }}(\mathit{\lambda }+{\mathit{\beta }}_{\mathrm{W}}^{\mathrm{2}})}},\end{array}$$

where $s_{\mathit{\delta }{}^{\prime }}$ and $s_{\mathit{\epsilon }{}^{\prime }}$ are constant error variances of the measured ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and heights respectively in the normalised (or, dimension-free) coordinate system, and λ is the ratio

$$\begin{array}{}\text{(24)}& \mathit{\lambda }={\displaystyle \frac{s_{\mathit{\epsilon }{}^{\prime }}}{s_{\mathit{\delta }{}^{\prime }}}}.\end{array}$$

The χ² minimisation of Eq. (23) with respect to β_W leads to the analytic expression (Carroll and Ruppert, 1996; Smith, 2009; Lolli and Gasperini, 2012) for the EIV slope parameter in terms of the variances ($s_{\mathrm{d}{}^{\prime }}$, $s_{\mathrm{h}{}^{\prime }}$) and covariances ($s_{\mathrm{h}{}^{\prime }\mathrm{d}{}^{\prime }}$) of the measured variables,

$$\begin{array}{}\text{(25)}& {\mathit{\beta }}_{\mathrm{W}}={\displaystyle \frac{s_{\mathrm{h}{}^{\prime }}-\mathit{\lambda }{s}_{\mathrm{d}{}^{\prime }}+\sqrt{(s_{\mathrm{h}{}^{\prime }}-\mathit{\lambda }{s}_{\mathrm{d}{}^{\prime }}{)}^{\mathrm{2}}+\mathrm{4}\mathit{\lambda }{s}_{{\mathrm{h}}^{\prime }{\mathrm{d}}^{\prime }}^{\mathrm{2}}}}{\mathrm{2}{s}_{\mathrm{h}{}^{\prime }\mathrm{d}{}^{\prime }}}}\end{array}$$

or

$$\begin{array}{}\text{(26)}& {\mathit{\beta }}_{\mathrm{W}}\equiv {\displaystyle \frac{\mathrm{1}-\mathit{\lambda }+\sqrt{(\mathrm{1}-\mathit{\lambda }{)}^{\mathrm{2}}+\mathrm{4}\mathit{\lambda }{s}_{{\mathrm{h}}^{\prime }{\mathrm{d}}^{\prime }}^{\mathrm{2}}}}{\mathrm{2}{s}_{\mathrm{h}{}^{\prime }\mathrm{d}{}^{\prime }}}}.\end{array}$$

since $s_{\mathrm{d}{}^{\prime }}=s_{\mathrm{h}{}^{\prime }}=\mathrm{1}$ from Eq. (18). And the covariance ($s_{\mathrm{h}{}^{\prime }\mathrm{d}{}^{\prime }}$) is computed using the standard definition,

$$\begin{array}{}\text{(27)}& s_{\mathrm{h}{}^{\prime }\mathrm{d}{}^{\prime }}={\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N\left(h_i{}^{\prime }{d}_i{}^{\prime }\right),\text{for}\text{large}N.\end{array}$$

Equation (26) can be solved if a prior 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 the calibration procedure described by Eqs. (17) and (18), the mean values of $s_{\mathit{\epsilon }{}^{\prime }}$ and $s_{\mathit{\delta }{}^{\prime }}$ are found to be 0.62±0.04 and 0.37±0.06 respectively (without any contour selection: Table 1). Since the EIV estimate of the slope requires only the ratio between $s_{\mathit{\epsilon }{}^{\prime }}$ and $s_{\mathit{\delta }{}^{\prime }}$, a good choice of this starting value is λ=1. Furthermore, for $\mathit{\lambda }=s_{\mathit{\epsilon }{}^{\prime }}=s_{\mathit{\delta }{}^{\prime }}=\mathrm{1}$, there is a simple geometric interpretation of the merit function in Eq. (23). This solution corresponds to minimising the Euclidean or orthogonal distance between the fitted line and the measured data. The residual function to be minimised with respect to the regression coefficients is

$$\begin{array}{}\text{(28)}& {\mathit{\chi }}^{\mathrm{2}}=\sum _{i=\mathrm{1}}^N{\displaystyle \frac{(h_i{}^{\prime }-{\mathit{\beta }}_{\mathrm{W}}{d}_i{}^{\prime }{)}^{\mathrm{2}}}{\mathrm{1}+{\mathit{\beta }}_{\mathrm{W}}^{\mathrm{2}}}}\equiv \sum _{i=\mathrm{1}}^N{\displaystyle \frac{(h_i{}^{\prime }-d_i{}^{\prime }{)}^{\mathrm{2}}}{\mathrm{2}}}\end{array}$$

since β_W=1 when λ=1 from Eq. (26). Following Eq. (22), we therefore have our first estimate of β in the scattered plot of ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and heights,

$$\begin{array}{}\text{(29)}& \mathit{\beta }={\displaystyle \frac{\sqrt{s_{\mathrm{h}}}}{\sqrt{s_{\mathrm{d}}}}}.\end{array}$$

Equation (29) is commonly referred to as a reduced major axis (RMA) solution in statistics literature (Smith, 2009). In practice, this is just the geometric mean (GM) of the two OLS estimates, ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$ and ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}$, as can be seen by combining Eqs. (12) and (13):

$$\begin{array}{}\text{(30)}& {\mathit{\beta }}_{\mathrm{GM}}=\sqrt{{\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}{\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{h}}}\equiv {\displaystyle \frac{\sqrt{s_{\mathrm{h}}}}{\sqrt{s_{\mathrm{d}}}}}.\end{array}$$

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 GM solution is the correct solution (discussed below), this first estimate reflects the nature of the biasing effect in the estimated slope in relation to the data selection process. A specific example of GM fitting is presented in Fig. 1a. The standard error in the OLS and GM slope estimate reported in Fig. 1a is from Vicente de Julián-Ortiz et al. (2010). Due to the relatively low detections with Sodankylä radar (Fig. 2c), the 2σ error in the estimated temperature using the GM solution is found to be significantly higher (13 K on average at contour level 0). This 13 K of noise level in the temperature can be reduced by a factor of $\sqrt{\mathrm{3}}$ or $\sqrt{\mathrm{5}}$ if a 3 or 5 d running mean of temperature is estimated with this radar. However, to test the robustness of the proposed method, this paper has estimated the daily averaged temperatures.

The systematic offset between MR temperature and co-located lidar temperature as a result of using Eq. (30) is reflected in Fig. 3a (red curve) and b (red histogram). Each of these temperatures are then compared with lidar temperatures at 88, 89 and 90 km for the dates when lidar data are available. 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 and 91 km in these data. Figure 3b (red histogram) reveals that the MR temperatures are overestimated by a mean value of +58 K for the case of the 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 (Fig. 2a) is more negative than expected. If these values are shifted by a constant positive offset of +1 K km⁻¹, 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 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 6 K at contour level 0. Hence we expect the day–night difference in Sodankylä 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 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 +58 K offset seen in Fig. 3b (red 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 values. For the data set used in this work, we found that the correlation between ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height is typically 0.50±0.05, 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 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 increases the r to be typically around 0.66±0.06, with the consequence of reduced biasing in the GM solution. Such a contour selection process essentially removes the erroneous data at higher and lower altitudes from the tail of the distribution, whereby the assumption of equal error variances (i.e. $s_{\mathit{\epsilon }{}^{\prime }}\approx s_{\mathit{\delta }{}^{\prime }}$) is achieved in the normalised coordinates (Table 1). In other words, the validity of GM solution is conditional upon how close the parameter λ is to 1 for a given data selection process.

4 Results and discussion

The convergence of λ towards 1 near the core of this bivariate distribution is evident from Table 1. This is demonstrated in Fig. 4a, where we have estimated the GM slope at various contour levels between 0 and 0.7. In addition, Fig. 4b reflects the asymptotic behaviour of GM solution at higher contour levels in normalised coordinates. For these data, beyond the contour level 0.4, any change in slope at higher contour level is within the 2σ error limit. This error in the slope corresponds to an average 2σ error of 11 K in the estimated temperatures at contour level 0.4. In addition, the estimated temperatures can be biased due to small variation of λ from 1 on a day-by-day basis. For example, If the true value of λ is 1.25, the GM slope will be overestimated by 4 % (from Eq. 26 for the date 14 November 2015). For a typical winter temperature of 200 K at 90 km, a 4 % offset translates to overestimation of temperatures by 8 K. In principle, this bias can be further reduced by selecting a higher contour level than 0.4, with the consequence of increased noise level in the estimated temperatures. For these data, the contour level 0.4 is found to provide an optimum condition such that a maximum of 25 % uncertainty in the parameter λ leads to 4 % bias in temperature, which, in turn, is comparable to the standard error in the temperature from regression analysis.

Figure 4 GM solution at different contour levels in (a) original coordinate and in (b) normalised coordinate for the date 14 November 2015. The vertical dashed lines in (a) correspond to the average value of λ obtained from SCT calibration at contour level 0.2 and 0.4 respectively for winter 2015–2016. The error bar in (a) corresponds to 2σ error .

We have applied the GM slope estimate at contour level 0.4 to the MR data from the period October 2015 to March 2016. This is presented in Fig. 5, along with the data from co-located lidar observations. The differences between MR temperatures and lidar are shown in Fig. 6a for altitudes near peak meteor counts. The mean difference between lidar and MR temperatures is +12 K, which is expected due to the variation in the true value of λ from 1. The root-mean-square (rms) difference is about 21 K. This difference is partly due the intrinsic errors of 5–10 K in lidar temperature and partly due to the statistical noise in the estimated MR temperatures.

Figure 5 Comparison of the bias-corrected MR temperatures with lidar data for the winter 2015–2016. The solid line corresponds to the temperature estimated using the GM solution at contour level 0.4. The dashed line corresponds to the SCT-calibrated temperatures using the co-located lidar measurements. The OLS estimates are obtained with ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ as an independent variable. The errors in lidar temperatures are 5–10 K, and the 2σ error (grey shade) of the temperature from EIV analysis is on average 11 K. The differences between lidar and MR temperatures are presented in Fig. 6.

Figure 6 Difference between MR temperatures and lidar data for (a) the GM solution at contour level 0.4 and (b) SCT calibration applied to OLS estimate of ${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$. MR temperatures are shown in Fig. 5.

For direct comparison with the results above, we have also estimated the MR temperatures using the revised SCT calibration procedure described in Sect. 3.1. For this we have estimated the OLS slope (${\mathit{\beta }}_{\mathrm{OLS}}^{\mathrm{d}}$) and used $\sqrt{s_{\mathit{\delta }}}=\mathrm{0.11}$ (Table 1) to obtain the calibrated temperatures from Eqs. (16) and (6). These calibrated temperatures are presented in Fig. 5. The histogram in Fig. 6b shows the differences between lidar data and the SCT-calibrated MR temperatures. The mean difference between the MR and lidar temperatures is again about −3 K, thereby showing that the biasing effect has been properly corrected by this calibration procedure. The rms difference is 15 K. Although the temperature estimated from the GM solution is slightly biased as compared to that estimated using SCT calibration, the presence of artefacts in the latter is clearly visible in Fig. 5. As evident from Fig. 7, the GM solution based on contour selection improves the temperature estimation significantly as compared to the traditional use of the OLS regression analysis.

Figure 7 (a) Improved temperature estimation using the GM solution (red) as compared to OLS estimates (blue and green) at contour level 0.4. (b) Reduced mean offset between MR and lidar temperature for GM slope estimate (red) as compared to OLS estimate (green and blue).

The EIV analysis does not distinguish between the measurement error and natural geophysical variability (e.g. Sprent, 1990, p. 13). In other words, the error variances, $s_{\mathit{\delta }{}^{\prime }}$ and $s_{\mathit{\epsilon }{}^{\prime }}$, consist of both measurement errors and the natural geophysical variation. The effective value of λ in a normalised coordinate may vary on a day-by-day basis and may be radar-dependent. For example, meteor trails can get modified by wind effects, ion composition, meteor fragmentation and strong ionospheric currents, as well as temperature and pressure fluctuations, on various spatial and temporal scales (Hocking, 2004; Younger et al., 2014). Despite the contour selection process, asymmetric effects of geophysical variation may have increased the effective variance of h_i^′, leading to overestimates in the GM solution (Gillard and Iles, 2005). Due to the high-latitude location of Sodankylä radar, the geomagnetic effect above 95 km can contribute to the systematic bias. Below 85 km, the decay of meteor radar echoes may deviate slightly from diffusion-only evolution (Lee et al., 2013), thereby requiring a better physical model that does not assume the linearity of Eq. (6). Assessing the contribution of geophysical variability at various altitudes would require carefully designed replicates of observations as well as long-term comparison of MR temperatures with other co-located instruments. In other words, the case for λ≠1 needs to be handled with careful modelling of errors by taking into account the dominant effect of geophysical variability in these data. This remains a subject of future research.

5 Summary

The biasing effect in MR temperature has been a pressing issue for the last 2 decades. Attempts have been made in the past to correct the slope in the scattered plot of ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ 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 a biasing effect, which is mainly due to the presence of various error terms. We have reviewed the conventional calibration procedure (1) and then provided an alternative method (2) for estimating MR temperature that does not 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 co-located 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), 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 }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height using co-located 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 were clearly visible in the estimated temperatures.

2. As an alternative method, we have applied the errors-in-variables (EIV) regression analysis to estimate the slope in the scattered plot of ${\log }_{\mathrm{10}}(\mathrm{1}/\mathit{\tau })$ and height. The model error in EIV analysis takes into account the total errors variances in both abscissa and ordinate. It is observed that the geophysical variability dominates at all altitudes as compared to measurement errors and is the key factor in addressing the biasing effect. Moreover, any asymmetry in the error variance is minimal near the meteor peak region. This allows for an independent estimate of weighted-averaged atmospheric temperatures at 90 km using a suitable contour selection procedure. The temperatures estimated using this method show very good agreement with co-located lidar measurement and with reduced systematic offset as compared to the traditional least-squares analysis.