Bayesian Filtering in Incoherent Scatter Plasma Parameter Fits

Incoherent scatter (IS) radars are invaluable instruments for ionospheric physics, since they observe altitude profiles of electron density (Ne), electron temperature (Te), ion temperature (Ti) and line‐of‐sight plasma velocity (Vi) from ground. However, the temperatures can be fitted to the observed IS spectra only when the ion composition is known, and resolutions of the fitted plasma parameters are often insufficient for auroral electron precipitation, which requires high resolutions in both range and time. The problem of unknown ion composition has been addressed by means of the full‐profile analysis, which assumes that the plasma parameter profiles are smooth in altitude, or follow some predefined shape. In a similar manner, one could assume smooth time variations, but this option has not been used in IS analysis. We propose a plasma parameter fit technique based on Bayesian filtering, which we have implemented as an additional Bayesian Filtering Module (BAFIM) in the Grand Unified Incoherent Scatter Design and Analysis Package (GUISDAP). BAFIM allows us to control gradients in both time and range directions for each plasma parameter separately. With BAFIM we can fit F1 region ion composition together with Ne, Te, Ti and Vi, and we have reached 4 s/900 m time/range steps in four‐parameter fits of Ne, Te, Ti and Vi in E region observations of auroral electron precipitation.

primary energy spectra of precipitating electrons. However, N r equals N e only if T e = T i , which may be an unjustified assumption when the precipitation heats the electron gas.
We propose an IS analysis technique that combines Bayesian filtering (for example Särkkä, 2013) in time and correlation priors (Roininen et al., 2011) in range. The combination allows us to extend the idea of full-profile IS analysis (Holt et al., 1992;Hysell et al., 2008;, which assumes smoothness in range, to an assumption of smoothness in both time and range. With this approach we can fit ion compositions if both ion temperature and composition are assumed to vary smoothly with time and altitude, and we can include temperature fits in high-resolution electron density fits.
In Section 2 we give introduction to IS plasma parameter fits, Bayesian filtering and correlation priors. In Section 3 we explain how the prior models and Bayesian filtering are used in IS analysis and implemented as a "Bayesian Filtering Module" (BAFIM) in Grand Unified Incoherent Scatter Design and Analysis Package (GUISDAP; . In Section 4 we demonstrate BAFIM fits of N e , T e , T i , V i , and ion composition p = [O + ]/N e in the F 1 region, and high-resolution fits of N e , T e , T i , and V i in the E region.

Theoretical Background
Incoherent scatter signal from a small plasma volume is a zero-mean random process with autocorrelation function R(τ), where τ is time lag. IS radar data are discrete samples of the autocorrelation function at discrete ranges r i , times t j , and lags τ k . Power spectral density of the scattered signal, which is the Fourier transform of the autocorrelation function, is a known function of plasma parameters (for example Swartz & Farley, 1979, and references therein).
Typically, plasma parameters are extracted from the autocorrelation function samples by nonlinear least squares methods with optimization techniques such as Levenberg-Marquardt algorithm. Alternatively, Markov chain Monte Carlo methods can be used for parameter extraction (for example Virtanen et al., 2014), although optimization has remained as academic standard in IS analysis.

Gated Analysis and Full Profile Analysis
IS analysis techniques can be roughly divided into "gated" and "full-profile" techniques. In gated analysis one runs the fitting process for each range r i and time t j independently from the analysis of neighboring observational volumes. The EISCAT IS analysis tool GUISDAP  makes gated analysis. In full-profile analysis one fits range profiles of plasma parameters. Main benefit of the full-profile analysis is the possibility to include prior information of plasma parameter altitude profiles.
In its most general form the full-profile analysis performs also deconvolution of lag profiles (Holt et al., 1992;Hysell et al., 2008). A simpler approach is to use phase-coding, for example alternating codes (Lehtinen & Häggström, 1987), and to decode the autocorrelation function samples into high resolution before the plasma parameter fit . The two-stepped approach can be accomplished with arbitrary transmission modulations if the deconvolution is performed by statistical inversion (Virtanen et al., 2008(Virtanen et al., , 2009. It is technically possible to add prior information already in the lag profile inversion step, but expressing the prior in terms of the actual plasma parameters is difficult in this approach.

Bayesian filtering and Smoothing
Bayesian filtering (for example Särkkä, 2013) is a class of methods for estimating the state of a system from noisy indirect measurements. In IS analysis the state of the system reduces to point estimates of plasma parameter values and their standard deviations, while the indirect measurements are the observed autocorrelation function samples R.
The filtering consists of a sequence of prediction and update steps. The sequence starts from an initial set of parameters  1 x and its covariance matrix  1 P , which form our prior understanding of the unknown VIRTANEN ET AL. 10.1029/2020JA028700 2 of 13 parameters at time t 1 . Autocorrelation function samples R 1 are then used to update the prior model into our best estimates of the parameters and their covariance at time t 1 , x 1 and P 1 . The update step is accomplished using a measurement model M, The update step is followed by a prediction step, in which x 1 and P 1 are combined with our best understanding of dynamics of the system to create our best prediction of the parameters and their covariance at time t 2 ,  2 x and  2 P . The prediction step is accomplished using a dynamic model D, Measurements from time t 2 are then used to update the prediction into the final estimates x 2 and P 2 , etc.
The simplest "dynamic" model is to assume that the parameter values at subsequent time steps are close to each other, which reduces the prediction step into where Q is the system noise covariance matrix. The larger values Q has in its diagonal, the smaller is the correlation between subsequent state estimates and the larger is the filter gain.
Bayesian filtering allows one to recursively estimate unknowns using the whole time history of measurements. In Bayesian smoothing the idea is extended to use of also "future" measurements. Bayesian smoothing reduces variances of the unknown parameters and guarantees that equal amount of information from "past" and "future" measurements is included in each estimate of the unknowns. This removes a time shift that may be produced by a low-gain filter.
If the dynamic and measurement models are linear functions, Bayesian smoothing can be implemented as a recursive smoothing step called Rauch-Tung-Striebel (RTS) smoother (Rauch, 1963). The smoothing recursion runs backwards in time using equations where D j is the theory matrix of the linear dynamic model D and the superscript T denotes matrix transpose.
s j x and s j P form the Bayesian smoothing solution of the problem.

Correlation Priors
Correlation priors (Roininen et al., 2011) allow one to model mutual covariances of the unknowns of an inverse problem in a well controlled way. Assuming that our prior belief of the unknowns x is x p , the prior can be expressed as a linear inverse problem where Σ p is the error covariance of ɛ p , Ω p is the precision matrix, x p ′ is the final prior mean, and Σ p ′ is its error covariance matrix. It is important to notice that the initial profile is smoothed by the correlations and x p ′≠x p . In high-dimensional problems it is important that Ω p is a sparse matrix (Norberg et al., 2018).
The zeroth order part of the prior is where the diagonal error covariance matrix Σ p,0 contains the prior variances of x p . The first order terms are And the second order terms are The full prior covariance matrix Σ p is Variances of the zeroth, first and second order terms are (Roininen et al., 2011), where α i is correlation power in the i (th) range gate, Δh i is width of the i (th) range gate, and ℓ i is the corresponding correlation length. The correlation lengths define how smooth the profile is, and the correlation power defines width of the prior distribution. The constants c 0 , c 1 , c 2 define shape of the final covariance structure. For example, c 0 = 1, c 1 = 1/2, c 2 = 1/8 produces a Gaussian covariance. The model variances depend on the discretization and correlation length in a way that makes the model essentially grid-independent.

BAFIM Implementation
We have implemented an IS analysis tool based on Bayesian filtering in time and correlation priors in range as an additional Bayesian Filtering Module (BAFIM) to the GUISDAP IS analysis tool . We assume a five parameter fit of electron number density N, ion temperature T, ion-to-electron temperature ratio E, line-of-sight plasma velocity V and ion composition O in this section to simplify the equations. The vector of plasma parameters at time step j is where N j is the electron density profile in range gates i = 1, …, M, and the vectors of the other parameters are defined similarly.
The analysis starts from an initial guess of the plasma parameters at time t 1 ,  1 x , and their covariance  ( ) diag σ P are used as a prior in a normal GUISDAP fit to measurements R 1 . The GUISDAP fit is the update step of the Bayesian filter. The gated GUISDAP analysis does not produce a full error covariance matrix of x 1 , but the error covariance matrix P 1 contains mutual correlations of plasma parameters in each range gate.
After the first time step, priors for the following GUISDAP fits are not taken from the IRI model, but the fit results from t 1 are used to predict the parameters and their covariance at t 2 . The predicted values  2 x and diagonal of  2 P are used as prior mean and variance to fit x 2 and P 2 to measurements R 2 , x 2 and P 2 are used to calculate the predicted  3 x and  3 P , etc. The analysis steps are illustrated in Figure 1, whose first row shows predicted altitude profiles of The predicted values and variances form a Gaussian prior distribution in a GUISDAP fit, which produces the updated profiles on the second row. The profiles on the second row are used to predict the parameter profiles at time t j+1 (third row), the prediction is used as a prior when fitting the parameters at time t j+1 (fourth row), etc. Correlations in range are lost and reintroduced in each update and prediction step, correspondingly. This allows us to use the computationally light-weight gated analysis, and the approach is acceptable if the plasma parameters do not change much during a time step.
In the prediction step, a correlation prior is used to create smooth plasma parameter profiles. The measurements x 1 and their covariance P 1 are used as the zeroth order terms in (8), The first and second order differences in (8) The first order difference matrices (13) for each parameter are identical M × M − 1 matrices, A N,1 = A T,1 = A E,1 = A V,1 = A O,1 , and the full first order difference matrix is the block diagonal matrix  The second order differences and their variances are formed in a similar manner. As a result, we have a matrix equation of the from (8), from which parameter profiles smoothed in range, x′, and their covariance, Σ′, can be solved using (9) and (10).
The smoothed parameter profiles x′ are used as the prediction for time step t 2 , and the predicted covariance is the sum of the covariance of x′ and a process noise covariance Q, The process noise covariance is a diagonal 5M × 5M matrix with a different variance for each plasma parameter (35) in its diagonal, , , , , , , , , , , , , , .
The RTS smoother is implemented in BAFIM as a postprocessing step. Since only the first 5M elements of the vector m p are nonzero in (8) and (10), the matrix D in (5) consists of the first 5M columns of the The RTS smoother is only a linear approximation, but the approximation is reasonable if the time steps are short enough to keep changes in plasma parameters small in between subsequent time steps.
The correlation lengths ℓ i are proportional to the plasma scale height calculated from the IRI model parameters. Here k B is the Boltzmann constant, m i is the mean ion mass, g i is the acceleration of gravity, and the subscript i refers to the i (th) range gate. The correlation lengths of N are where h N s is a constant, and the correlation lengths of the other parameters are defined in a similar manner.
In the correlation prior, covariance of the zeroth order terms is the posterior covariance Σ p,0 = P 1 , and variances of the first and second order terms are proportional to   2 i and   4 i , respectively. Thus, at the limit of small correlation lengths ℓ i , the smoothed profile x′ approaches the fitted profile x 1 , and the covariance Σ′ approaches P 1 . BAFIM can thus be run without the smoothing in range if the correlation lengths ℓ i are small, i.e. the constants s h are small.

Plasma Parameter Fits with BAFIM
In this section we demonstrate plasma parameter fits with BAFIM in two use cases, ion composition fits in the F 1 region and high-resolution E region analysis during auroral electron precipitation. We use fieldaligned observations from the EISCAT Svalbard radar (ESR) and the EISCAT ultra high frequency (UHF) radar. We consider fits of electron density (N e ), ion temperature (T i ), ion-to-electron temperature ratio (T r ), line-of-sight plasma bulk velocity (V i ), and ion composition (p = [O + ]/N e ). In the results we show the electron temperature T e = T i ⋅ T r instead of T r . While the assumption of smoothness in range is necessary in the selected demonstrations, we emphasize that BAFIM can be used also without this assumption, for example to improve time resolution of four-parameter fits in low-elevation or bistatic observations. In this section, standard GUISDAP fits and GUISDAP fits with BAFIM are referred to as "GUISDAP" and "BAFIM", correspondingly.
Both ESR and UHF data are from experiments that use alternating codes (Lehtinen & Häggström, 1987). The ESR "ipy" experiment uses a 30-bit code sequence with 30 μs bit length and the data are decoded to 2.25 km resolution. The UHF "arc1" experiment uses a 64-bit code sequence with 6 μs bit length and the data are decoded to 900 m resolution. In high signal-to-noise conditions GUISDAP may underestimate plasma parameter variances because it neglects correlations between autocorrelation function samples . Both experiments use randomized (Lehtinen et al., 1997) codes to reduce the correlations. If highly correlated data were analyzed with BAFIM, smoothing in time and range would be reduced due to the underestimation of errors in the GUISDAP fits.

Ion Composition Fits
Ion frictional heating occurs when an electric field drives the ionospheric plasma through the neutral atmosphere and the ion gas is heated in collisions with neutral particles. The heating may affect F 1 region ion composition, because reaction rates of some important charge-exchange reactions depend on temperature, and expansion of the neutral atmosphere may increase neutral N 2 concentration in the F region (Kelly & Wickwar, 1981). Deviations from the IRI ion composition may bias F 1 region ion temperature estimates in four-parameter GUISDAP fits of N e , T e , T i and V i . An example of such an event is shown on the left in Figure 2, where four-parameter GUISDAP fit results with 60 s resolution are shown for 24 h of ESR data. Ion temperature (third panel on the left) has an artificial local maximum around 200 km altitude, where IRI predicts too much molecular ions (fourth panel).
In five-parameter BAFIM fit of the same data ( Figure 2, middle panels), also the ion composition p is fitted, and the analysis proceeds with 6 s time steps. Other BAFIM settings are listed in Table 1. The artificial ion temperature maximum, which is visible in the GUISDAP fit, is not produced in the BAFIM fit. The transition altitude, where number density of molecular ions is equal to O + density (p = 50%, black lines in the fourth panels), is generally lower than in the IRI model. Difference of the two fit results (GUISDAP -BAFIM) is shown on the right in Figure 2, where one can see how the difference in p affects also T i , T e and even N e profiles. While the artifact around 200 km altitude was removed by BAFIM, the true ion frictional heating events between 4 and 7 UT, as well as the weaker T i enhancements after 15 UT, are reproduced by BAFIM, demonstrating its ability to maintain true ion temperature maxima. We note that our results are very similar with those of Blelly et al. (2010), who used the same data to demonstrate a full-profile analysis technique based on ion energy equations.

High-Resolution Observations of Auroral Electron Precipitation
IS radars can detect impact ionization and electron heating caused by auroral electron precipitation. While existing high-latitude IS radars can typically reach a time resolution of some tens of seconds in the four-parameter fits of N e , T i , T r , and V i , optical observations show that the precipitation may change substantially in a few seconds and even below (for example Dahlgren et al., 2016). High-resolution E region observations often rely on raw electron densities (for example Dahlgren et al., 2011;Semeter & Kamalabadi, 2005;Virtanen et al., 2018), which are calculated assuming T e = T i . However, this assumption may not be justified, since the precipitation heats the electron gas.  Table 1. While plasma parameters from the GUISDAP fit are extremely noisy with the 4 s/900 m resolution, the BAFIM fit produces temperatures and velocities that match well with the standard coarse-resolution fit (for example, compare T i and V i in panels (a), (b), and (d)). Figure 3(c) which shows raw electron density N r , BAFIM-fitted N e , relative difference (N e − N r )/N r , and the temperature ratio T e /T i . The raw densities are VIRTANEN ET AL.

Importance of the temperature fit is demonstrated in
10.1029/2020JA028700 9 of 13   Virtanen et al. (2018), because the high-resolution four-parameter fits were practically impossible.

Discussion
BAFIM is the first implementation of Bayesian filtering to IS plasma parameter fits. In this section we discuss some important properties of BAFIM and potential future improvements.

Resolutions of BAFIM Fit Results
While the BAFIM analysis proceeds with short steps in range and time, each fit of plasma parameters (21) contains information from longer intervals because the steps are correlated. Exact "effective" resolutions cannot be easily calculated, since the correlation prior Equation 8 is nonstationary, the resolutions depend on measurement noise, and neglecting the error correlations in the GUISDAP implementation distorts VIRTANEN ET AL.
10.1029/2020JA028700 10 of 13 second moments of the posterior distribution. However, the correlation lengths in range (34) are known, and we can estimate the physical correlations in time from the fit itself. For example, in the ion composition fit in Section 3, the correlation lengths vary from 0.5 km (N e in the E region) to 30 km (T i in the F region), and random fluctuations in fitted N e , T i , T e , V i , and p are uncorrelated in time scales longer than, 6 s, 12 s, 24 s, 24 s, and 5 min, correspondingly.

Tuning and Validating BAFIM
Tuning the process noise variances and correlation lengths of BAFIM may be nontrivial, since the correlations in time allow part of the prior information introduced with the correlation priors to be passed from one time step to another. Any change in process noise variance must thus be compensated with a corresponding change in correlation length to keep the effective smoothing in range unchanged. In addition, changing the process noise and correlation length of one plasma parameter may affect the others due to error correlations.
In this paper, BAFIM was tuned to produce practically uncorrelated electron densities, while correlation lengths and process noise variances of the other parameters were selected in such a way that noise level of the fitted parameters roughly matched with the default GUISDAP fits with 60 s resolution. The only physics-based part of the model are the correlation lengths, which are proportional to the plasma scale heights. Physics-based, automatic ways to tune the filter will be topics of future works. Alternative ways to tune the filter would be to derive theoretical limits for gradients in space and time, or to extract information on the correlation structures from existing measurements. Correlation structures of mesospheric winds have been extracted from meteor radar observations by Vierinen et al. (2019), and a similar work for incoherent scatter radars could be possible.
Validation of BAFIM results, the ion composition fits in particular, is a challenging task due to lack of measurements from other instruments. Observations of F 1 region ion composition are mainly from rockets, and the rocket observation would need to be from vicinity of the radar beam to enable reasonable comparisons. Alternatively, one could analyze simulated radar data corresponding a realistic model ionosphere. Such simulations would be possible for example with the simISR tool (Swoboda et al., 2017).

Ion Composition Fits
In the ion composition fits a small process noise variance q O was used for the ion composition and a relatively large variance q T was used for the ion temperature, which is equivalent with the assumption that ion temperature varies much more rapidly than ion composition. Only slow variations in composition were allowed, because allowing rapid variations in both ion composition and temperature may lead to unrealistic oscillations due to the temperature-ion composition ambiguity. With the selected tuning BAFIM can follow the relatively slow ion composition variations associated with the large scale convection electric field, but rapid variations caused, for example, by small scale electric fields around auroral arcs are challenging.
Time resolution of the composition fits could be improved if physics-based models were included in the prediction step. One could either model the temperature profiles or include a chemistry model that solves temperature-dependent compositions. The temperature profiles could be modeled, for example, with the techniques of Blelly et al. (2010) and Zettergren et al. (2011) while chemistry modeling could be adopted for example from Richards and Voglozin (2011). Also D region ion composition and temperatures could be observed if a sufficient model, for example the Sodankylä Ion and Neutral Chemistry (SIC) model (Turunen et al., 2016) was used.

EISCAT_3D
EISCAT_3D (McCrea et al., 2015) is the next-generation geospace radar system currently being built in northern Norway, Sweden, and Finland. The radar will provide an order-of-magnitude improvement in measurement speed, and it will be the first multistatic, multibeam incoherent scatter radar system. EI-SCAT_3D will be able to conduct volumetric observations, including three dimensional observations of plasma flows.
If BAFIM-like analysis was applied to field-aligned EISCAT_3D measurements, the order-of-magnitude improvement would mean subsecond time steps in four-parameter fits, and resolutions sufficient for rapidly varying conditions in association with aurora in ion composition fits. The volumetric observations would allow one to implement three dimensional models of the ionosphere in the prediction step. An EISCAT_3D analysis tool could be designed for the volumetric observations and could make optimal use of the multistatic, multibeam data, following the idea of Virtanen et al. (2014).

Conclusions
We have introduced an incoherent scatter analysis technique that allows us to control plasma parameter gradients in both time and space using Bayesian filtering and correlation priors. The technique is implemented as a Bayesian Filtering Moculde (BAFIM) in the GUISDAP analysis package. BAFIM allows us to to fit F 1 region ion compositions and transition altitudes, and to include ion and electron temperatures in high resolution plasma parameter fits, in field-aligned incoherent scatter measurements. Improvements provided by the new analysis tool were demonstrated with EISCAT radar data, including fits of F 1 region ion composition and high-resolution E region plasma parameter fits during short-lived auroral precipitation events. The technique could be extended to volumetric, multistatic observations of the EISCAT_3D radar and supplemented with ion chemistry models.

Data Availability Statement
The EISCAT data and the GUISDAP software are available for download from the EISCAT web page (http:// www.eiscat.se). BAFIM is available at https://doi.org/10.5281/zenodo.4033904.