Estimation of aerosol particle distributions with Kalman Filtering – Part 1 : Theory , general aspects and statistical validity

Introduction Conclusions References


Introduction
Atmospheric aerosol particles have significant effects on visibility (Hand and Malm, 2007), cloud formation (McFiggans et al., 2006), atmospheric radiative transfer (Myhre, 2009), and public health (Pope and Dockery, 2006;Gurjar et al., 2010).According to the Intergovermental Panel on Climate Change (IPCC; Forster et al., 2007), uncertainties related to the direct and indirect climate effects of aerosols are a significant uncertainty factor in the climate change assessment.Both particle size and chemical composition largely determine their climatic impacts, and are thus important to be accurately characterized.
Atmospheric particles can cover several orders of magnitude in size, and contain various chemical species in internal and external mixtures.Consequently, there is no single instrument capable of measuring the entire range of aerosol quantities (McMurry, 2000).Rather, a number of instruments are needed, each providing information on measurable quantities, such as particle electrical mobility or light scattering intensity.In addition, mathematical techniques are often needed to invert the raw observations into physical and/or chemical particle properties (e.g., Kandlikar and Ramachandran, 1999;Fiebig et al., 2005).Therefore, even though obtained particle properties are related to each other, combining them into a unique conceptual framework which accurately describes the aerosol state has proven to be very challenging.For example, several algorithms have been developed to determine the aerosol size distribution state starting from independently measured quantities from optical, aerodynamic and electrical mobility detectors (Hand and Kreidenweis, 2002;Shen et al., 2002;Khlystov et al., 2010).However, assumptions on aerosol particle shape, density or chemical composition are needed to finally obtain the best possible closure between the measurements.In this approach, the errors introduced by each instrument are not easily accounted for and does not guarantee a physically justified continuity of the obtained solution.
Here we will introduce data assimilation as an option in aerosol physics to obtain a consistent state estimate based on measured aerosol properties.Data assimilation is a mathematical framework in which information from different sources is blended to obtain a maximum-likelihood, or minimum-variance, state estimate.It is widely used in geosciences and engineering, for example in numerical weather prediction (Lorenc, 1986;Rabier et al., 2000).In atmospheric chemistry and aerosol modeling, data assimilation has been applied in the context of air quality (Elbern et al., 2001), and environmental monitoring (Dubovik et al., 2008;Chung et al., 2010).Lately, it has also been used to improve aerosol mass, aerosol optical depth or extinction profile measurements (e.g.Lin et al., 2008;Tombette et al., 2009;Sekiyama et al., 2010;Schutgens et al., 2010).In this article, data assimilation is, as far as we know, for the first time introduced as a potential unifying framework to estimate the particle number size distribution from multiple in-situ measurements.
Data assimilation always incorporates time-evolution models which carry the state estimate forward in time.Consequently, the obtained state tends to maintain a temporal continuity unlike a pure instrumental closure and is thereby physically more justified.In the references above, these models were 3-dimensional chemical transport or general circulation models.Focus in this article is on in-situ multiinstrument aerosol measurements.The modeling framework is thus a size segregated local 0-dimensional aerosol microphysical process model (a so-called "Eulerian box model").This enables very detailed time evolution of aerosol processes but lacks the spatial aspects, similar to point-wise atmospheric aerosol measurements.
Here, in Part 1, Extended Kalman Filter (EKF) was only applied with information from similar type of detectors with different measurement ranges that slightly overlap.This limited approach is chosen for three reasons: it allows (1) to detect how the general aspects of aerosol size distributions affect the EKF implementation, (2) experimenting on how the EKF implementation handles drastic changes in the aerosol size distributions, and (3) examination of the statistical validity of the EKF implementation.The motivation for this article is to demonstrate the applications of data assimilation in interpretation of in-situ observations and to discuss the strengths and weaknesses of this technique relative to inversion of the instrument kernel matrix technique.
In Part 2 (Viskari et al., 2012), the EKF implementation are extended to include information also from different instrument types simultaneously measuring different particle size ranges.

The Kalman Filter
Kalman Filter (KF; Kalman, 1960) is a sequential state estimation method for linear dynamical systems.Extended Kalman Filter (EKF; for text-book treatment, e.g., Kaipio and Somersalo, 2004) is a standard extended version of KF for non-linear systems.In both KF and EKF it is assumed that observation and model errors are Gaussian with zero mean.KF and EKF both proceed from one observation time to the next in two alternating steps: time evolution and observation updating.The difference between KF and EKF is that in EKF the error propagation must be assumed to be dominated by tangent-linear terms while in KF it is known to be linear.
The notation used here was defined in Ide et al. (1997).At the time evolution updating, the prior state estimate, also referred to as the background state estimate, x k at time k is obtained by propagating in time the preceding state estimate x a,k−1 at time k − 1 using Here M is the non-linear time-evolution model.Note, that the prior state estimate x k is the best estimate at time k before observations are available.The state estimate x a,k−1 at time k − 1, on the other hand, is a posterior estimate after the observations are available.
The error covariance B k of the prior state estimate x k is obtained using The term MB a,k−1 M T corresponds to the tangent-linear time-evolution of the error covariance.Note, that M is tangent-linear with respect to M. The model error term Q represents the errors due to model discrepancies and system noise, excluded from M. It should be noted that for EKF, Eq. ( 2) is the approximated B k as the non-linear evolution of the B k−1 is not accounted for.At the observation updating, the posterior state estimate x a,k at time k is obtained by optimally combining the prior state estimate x k and the observations y k at time k using Here H is a possibly non-linear observation operator, which produces a simulated observation corresponding to the prior state.The observation minus model counterpart term is called innovation.The Kalman gain K k is defined as Here O k is an observation error covariance matrix and H is a tangent-linear version of H .The Kalman gain optimally weights observations and prior state based on their respective accuracies, and spreads the innovation to the different elements of the vector x a,k .The error covariance of the state estimate x a,k is The error variances of the posterior state estimate are smaller than in the prior due to introduction of new information contained in the observations.In essence, EKF estimates the new state by correcting the time-evolved state with the recent observations.The relative accuracies of the observations and prior state determine the impact of the observation.More weight is laid upon the more accurate information sources.In EKF, there are two implicit information sources, on top of the observations: the prior state, and the model.The prior state carries the information of all the past observations forward in time.The evolution model encapsulates our physical and chemical insight of the system in a compact way, and the observation operator characterizes the measurement event.
3 Information sources

Observations
Particle number mobility distributions were measured using a Differential Mobility Particle Sizer (DMPS; Hoppel, 1978).The DMPS classifies particles according to their electrical mobility as a function of voltage, which can further be related to the characteristic particle size.A Differential Mobility Analyser (DMA) is used for particle classification and a Condensation Particle Counter (CPC) for subsequent particle counting.The smallest, 3 to 40 nm, particles were measured using a DMPS consisting of a 10.9 cm long Vienna type DMA (Winklmayr et al., 1991) and a CPC model 3025 manufactured by the TSI company (Stolzenburg and Mc-Murry, 1991).The bigger, 10-1000 nm particles, were measured using a DMPS consisting of a 28 cm long Vienna type DMA and a TSI CPC model 3010 (e.g.Quant et al., 1992).The two DMPS systems were operated in parallel with an overlapping size range from 10 to 40 nm.By changing the voltage step-by-step in the DMAs, a full particle distribution was obtained every 10 min.Prior to the DMPS instruments, a radioactive bipolar neutralizer was used in order to obtain a steady-state charge distribution.The charging probability of the particles at different sizes was calculated according to Wiedensohler (1988).Details of the DMPS measurement system used are presented in Aalto et al. (2001).
The mathematical DMPS data inversion from raw observations to the geophysical quantities is performed as follows.The particle number size distribution x is obtained as a solution of the equation In the equation y r is the vector of raw measurements and R is an instrument specific kernel matrix.Essentially, R takes into account the particle charging probabilities (Wiedensohler, 1988), transfer functions of each DMPS channel (Stolzenburg, 1988) and particle size dependent losses, such as the calibration based CPC detection efficiencies and diffusion losses (Aalto et al., 2001).Different approaches (e.g.Stratmann and Wiedensohler, 1996;Stolzenburg and McMurry, 2008;Stolzenburg, 1988) for modeling and handling of the DMA transfer functions can be used in the DMPS inversion, but they lead to qualitatively similar results.Here, the theoretical transfer functions taking into account the diffusional broadening in the DMA (Stolzenburg, 1988) were utilized.
The inversion of the instrument kernel matrix solution was validated in Wiedensohler et al. (2012).
The solution for x in Eq. ( 6) is determined through the following steps: (i) the size distribution diameters are calculated as transfer function peak diameters for DMPS 1 and 2. The diameters for the overlapping size range are taken from DMPS 1 only.(ii) Transfer functions for both DMPSs are integrated separately for this diameter grid and combined as a transfer matrix R. (iii) The best solution for concentrations is found using a least-square nonnegative pseudoinverse method where solutions in the overlapping region are averaged from the two information sources.
The DMPS measurements deployed here are a part of the EUSAAR network, thereby following the network quality standards (Wiedensohler et al., 2012).It is important to note that the DMPS measures each channel at a separate time instant over a 10 min measurement period.The particle size distribution, however, evolves during this period and we assume here that the observations at different channels are made simultaneously at the nominal observation time.

The time-evolution model
The University of Helsinki Multi-component Aerosol model (UHMA) is a size-segregating aerosol dynamical box model of the state evolution of a particle size distribution.The model is discussed in more detail in Korhonen et al. (2004).Essentially UHMA contains the major microphysical atmospheric aerosol processes.The activation scheme (Sihto et al., 2006;Kulmala et al., 2006) and the adjusted Fuchs-Sutugin method (Fuchs and Sutugin, 1971;Lehtinen and Kulmala, 2003) were used to determine nucleation and condensation rates, respectively.The coagulation kernel is calculated with the Fuchs equation (Fuchs, 1964).The dry deposition velocity is determined according to Rannik et al. (2003) for 10-500 nm and extrapolated to particles smaller than 10 nm.Here, the model is discretized such that the size distribution is presented with 50 size bins, evenly spaced according to the logarithm of particle diameter.The smallest (largest) size bin is chosen to be 1.5 nm (2 µm) in diameter.For this study, the UHMA model variables are the particle number and volume concentrations in each size bin.The new particle formation is due to the nucleation process.Sudden changes in the particle size distribution due to external The model is initialized with a measured particle number size distribution.The initial particle composition is approximated as 40 % organic compound and 60 % sulphuric acid.After each observation update step, the model re-initialized with an estimated particle number size distribution.Ambient vapour concentrations used in the model are based on measured vapour concentrations.
A tangent-linear version of the UHMA model required by Eq. ( 2) was presented and validated in Viskari et al. (2008).

The observation operator for the DMPS instrument
The observation operator H computes the observed quantity corresponding to the model state, that is, the particle number concentration within a size bin.In the model output, the number concentration is given for discrete particle diameters according to the model discretization.The DMPS measures electrical mobility for a defined voltage which can be converted to a characteristic diameter corresponding to geometrical diameter for a spherical particle.The observation operator for DMPS H DMPS is therefore where an interpolation matrix P interpolates the model output to the characteristic diameters of the DMPS, i.e., interpolates from the model grid to the instrument grid.The instrument kernel matrix R, same in Eq. ( 6), converts the interpolated number concentrations at the characteristic diameters to number concentrations at the electrical mobilities of corresponding voltages.The kernel matrix used here is constant in time.
The actual interpolation was performed with a cumulative distribution function (CDF) of the number concentration.The CDF for each model size bin is obtained by a simple summation of number concentrations in the model grid.Then, 4th order Lagrange polynomials are used to interpolate the CDF values to the characteristic diameters of the instrument.Finally, particle number concentrations in the instrument grid are computed by differentiating the CDF.We note that both matrices in Eq. ( 7) are linear, and thus H is also a linear.
The observation operator was tested and validated by comparing raw observations to the values computed by H from a size distribution obtained with the inversion of the instrument kernel matrix.

Time evolution updating
The non-linear UHMA model is used to propagate the state estimate x a,k−1 from time k − 1 to k.The tangent-linear UHMA model, used in the error covariance evolution, has been shown to be valid for time ranges up to about 30 min (Viskari et al., 2008).As the observation intervals are 10 min in this study and the simulations are re-initialized after each interval, the tangent-linear hypothesis is considered valid in this study.Auxiliary measurements of ambient sulphuric acid and non-volatile organic vapour concentrations were used to support the simulation (Petäjä et al., 2009;Paasonen et al., 2010).The ambient vapour concentrations are required by the UHMA model for the intermediate times between measurements and these were computed from the measurements using linear interpolation.A volatile organic vapour was included in the simulation with concentrations parametrized from the non-volatile organic vapour concentrations according to Vuollekoski et al. (2010).Ambient atmospheric conditions (e.g., temperature, pressure) were specified as constant values due to the model being not particularly sensitive to these values.
According to Eq. ( 2), the error covariance B k of the prior state is a sum of the time-evolved error covariance B a,k and the error model error Q.There are no tools available to properly estimate the model error term, which was therefore neglected here.This simplifying assumption, though, does have adverse effects.In EKF, the posterior error covariance decreases at each observation updating.If the error source is neglected, the decrease due to observations may exceed the increase due to time-evolution of the error covariance.As a consequence, the state error covariance may become gradually smaller and smaller, as was the case here.This is both unrealistic and undesirable: observations will have less and less weight, and the filter eventually diverges from the observations.This problem is common for both KF and EKF applications and is known as filter divergence (Whitaker and Hamill, 2002).Proper inclusion of the model error would increase the prior state error covariance.
In absence of a sound method to account for the model error, an ad hoc solution is developed.Here, we scale the prior state error variances before the observation updating.A minimum relative error was chosen for the background error.If the relative error value of the background error is smaller than this limit, the relative background error in that particle size is increased to the set limit, preventing the background error variance from becoming too small.By choosing this limit to be larger than the relative error of the observations, it is ensured that EKF will not trust the model information more than the observations.The downside of the approach is that the chosen limit will not accurately represent the relative errors of the background state.During the scaling, the background error covariances are adjusted so that the background error correlations remain unchanged.
Next, we will elaborate the definition of the observation error in the context of EKF.Let us consider the difference between observations and the model state in observation space (Eq.3).Obviously, where ε k is the perceived mismatch, or "error" at time k.In addition to the background state error, it is composed of three components: (i) the observations contain instrument noise, (ii) the x k is not representative of all the phenomena contained in the observations, and (iii) the observation operator cannot exactly reproduce the instrument response.These three components form the observation error.It is possible to estimate the instrument error before forming the EKF, but representativeness and observation operator errors need statistical material, i.e. innovation sequences produced with EKF (Dee and Da Silva, 1999).
In this study, we assume the observation error to be relative to the measured number concentration, specifically 15 % for DMPS I and 12 % for DMPS II (P.Aalto, personal communication, 2012).We also assume that the observations errors are uncorrelated.In actuality the observation operator does cause observation error correlation, but those covariances are generally small enough for the observation error to be considered uncorrelated.The prior state error is scaled such that the relative error for any particle size in the measurement range cannot be smaller than 20 %.This specification ensures that the observations will always be weighted more than the background state in the posterior state estimate.Refinements of this aspect are left for future work.

Observation updating
The observation updating is implemented without major simplifications.In order to compute the Kalman gain, Singular Value Decomposition (SVD) was applied for the matrix inversion.Additionally, in order to avoid numerical instabilities, both the observation error covariance and the prior state error covariance in the observation space were preconditioned so that their error standard deviations were divided by the corresponding observed number concentration, i.e., the errors are relative to the corresponding observed values.
The errors in the small particle sizes, where majority of particles are newly formed, and in particle sizes, where the total surface area is large and thus the influence on the ambient vapour concentration is the strongest, i.e., where the condensation sink values are high, are strongly correlated (Viskari et al., 2010).At the beginning of a nucleation event, there is usually a large difference between the background and observed state, which due to the strong correlation cause a large deviation in the estimated number concentration for the particles in the size range of the highest condensation sink value.This deviation, in turn, reinforces the initial difference in the smaller particle sizes during the observation updating via the background error covariance.The EKF implementation ultimately becomes unstable.To avoid this, the "length scale" of the prior state error covariance was artificially re- stricted at the observation update to 15 closest size bins in maximum.In other words, a large innovation can only affect the size distribution to a maximum distance of 15 size bins.This restriction is known as variable localization (Hamill et al., 2001).
No observation quality control was applied, except that negative observations were discarded as erroneous during the observation updating.The main reason for this omission is that the inversion of the instrument kernel matrix, our benchmark, has no built-in quality control either.It is important to note, though, that we assume that the observation error is always symmetric.However, due negative number concentrations are impossible, for very small measured number concentration values the error cannot be symmetric.We do not further focus on this here, as when the error becomes asymmetric is beyond the scope of this paper.

Results and analysis
We report here a case study utilizing measurements from 7 May 2007 from the SMEAR II-station in Hyytiälä, Finland (Hari and Kulmala, 2005).The measurements were a part of the EUCAARI project (Kulmala et al., 2009).The date was chosen because the measurements show different dynamical events, such as a strong nucleation event, and non-dynamical events, such as an apparent change in the air mass.The ambient sulphuric acid and non-volatile organic vapour concentrations used in the EKF implementation were also measured and are presented in Fig. 1.
Applicability of EKF in estimating particle size distributions was tested as follows.A size distribution obtained by interpolating the inversion of the instrument kernel matrix method for 00:00 local winter time (LT) to the model grid was used as an initial state for EKF.deviation of the initial state was set to 30 % of the initial number concentration and uncorrelated errors were assumed.The relative large initial error was chosen so that the observations weight relatively more in EKF at the very beginning.
The particle size distributions obtained with EKF (denoted as x EKF ) and direct inversion of the instrument specific kernel matrix (cf.Eq. 6; x INV ) are presented from 00:00 to 23:00 LT for 7 May 2007 in Fig. 2. Solutions for x EKF (Fig. 2a) and x INV (Fig. 2b) are qualitatively close to each other, although x EKF appears smoother compared to x INV .The total number concentrations for particles larger than 3 nm for both size distributions are presented in Fig. 2c.The differences in the total number concentrations are partially due to the diameters for x EKF and x INV not being the same, which makes it difficult to set the lower limit for x EKF when calculating the total number concentration as the lowest diameter of x INV .At approximately 15:00 and 17:00 LT, the total number concentration is significantly larger for x EKF than for x INV .
Validation based on independent observations was not considered feasible due to lack of reliable independent observations.Instead, in order to better understand the characteristics of the state estimate x EKF and to test self-consistency of the implementation, we compared it to the DMPS raw observations y.A residual r = y − H x EKF has to meet two conditions: (i) the bias and standard deviation of r has to be in the equal or better than the residual computed from the inversion of the instrument kernel matrix, (ii) large values of r are either due to measurement noise or special circumstances (e.g., precipitation, change of air mass).It is important to note that the focus in this article is not the inversion solution and we do not address how well the inversion solution estimates that state.Being the standard method to determine the particle number size distribution, the solution x INV obtained by inversion of the instrument kernel matrix was chosen as a reference for the comparison.
In order to quantitatively analyze the results, both x EKF and x INV were converted with the observation operator to the DMPS measurement channels at every measurement time (i.e., every 10 min) and compared to the raw DMPS measurements.The bias (the systematic difference between the estimated size distribution and the raw measurements; Fig. 3a) and standard deviation (the random difference between the estimated size distribution and the raw measurements; Fig. 3b) are shown for both solutions over the time window 00:00 to 23:00 LT.The bias for x EKF is largest in the size range of 10 to 40 nm, which is the area where DMPS I and II overlap (hereafter called the overlap region).The bias of the x EKF with respect to the DMPS I and II are roughly equal but of opposite sign in this region.Our interpretation is that the biases are a consequence of the instruments not completely agreeing in the overlap region.For the rest of the size range, there is a relatively small positive bias for DMPS I at particle sizes 8 to 11 nm and a small negative bias at 50 to 200 nm for DMPS II.At other particle sizes the bias is virtually zero.In comparison, the bias for x INV is generally larger than for x EKF .In the overlap region, the bias of x INV is positive for both DMPS I and II.In addition, for x INV there are large positive biases at particle sizes 60 to 80 nm and 600 to 800 nm and a large negative bias at 80 to 500 nm.
The standard deviation (Fig. 3b) is approximately as large for x EKF and x INV in the smallest (below 10 nm) and largest (above 300 nm) particle sizes.With respect to DMPS I, the two solutions are of about equal quality.With respect to DMPS II, x EKF has a smaller (larger) standard deviation than x INV for particles smaller (larger) than 40 nm.Especially, at 50 to 200 nm, x EKF has notably larger standard deviation than x INV .
To better understand these results, the data set is divided into four specific time windows: 00:00-09:00 LT (the aerosol system state is quasi-stationary), 09:00-17:00 LT (a nucleation event affects the size distribution), 17:00-19:00 LT (a sudden change in the size distribution), and 19:00-23:00 LT (a possible recovery phase).These correspond to time windows I, II, III and IV (Fig. 2).The time separation is only to facilitate the data analysis; it allows a more detailed data analysis, but due to shorter time windows, there are less data in each window which somewhat decreases the statistical reliability.

Common aspects for all time windows
The bias and standard deviation for each time window are presented in Fig. 4.Both quantities are generally large in the overlap region.This is related to the fact that the measure- ments of DMPS I and II are different in this region.In all the cases the bias of x EKF is smaller than the bias of x INV , especially in I and IV.It is our interpretation that the x INV is not the optimal solution in the overlap region.This indicates that x EKF performs very well in the overlap region, and this is very promising regarding the multi-instrument retrievals.
For x INV , the bias with respect to DMPS I changes between time windows, whereas for DMPS II it remains approximately the same.There are positive biases at 60 to 80 nm and 600 to 800 nm ranges and a negative bias at 80 to 500 nm.The standard deviation for x INV is very small outside the overlap region, except in time window I, where there are large values in 80 to 300 nm range.In contrast, the standard deviation of x EKF is significant for large particles outside the overlap region.This is mainly because of (1) large random errors in observations which are filtered out from x EKF but included in x INV , and (2) a sudden change in state for which x INV readily adjusts to, unlike x EKF which carries forward also the past information.These aspects will be discussed below.

Sudden changes in the system state
The bias of x EKF outside the overlap region is mainly due to the "memory" of EKF, which is caused by inclusion of the prior state estimate (and thereby past observations) into the solution.This implies a reduced sensitivity to sudden changes in observations.It tends to filter out measurement noise and produce dynamically consistent solutions from one observation time to another and from one particle size to the next.However, there may appear sudden changes in the measured size distribution by external causes (i.e., atmospheric processes not included in the 0-dimensional evolution model) and the system may continue to evolve from the new state.In this kind of a step-wise evolution, the "memory" of the previous observations delays the adjustment of x EKF to new observations.In such cases, x EKF is biased, and due to the gradual adjustment, this bias is a function of time.Additionally, this "drift" appears both in the bias and in the standard deviation as a time dependent bias cannot be statistically removed from the data when determining standard deviation.The results presented here, especially for time window III, include this phenomenon.Note how the bias of x EKF at 40 to 200 nm range (Fig. 4, time window III) coincides with the large standard deviation.The sudden change in the system state in time window III is possibly due to an air mass change.
This phenomenon also explains the large standard deviation of x EKF in time window IV. Figure 2b reveals how the number concentration at 50 to 300 nm range suddenly increases and then decreases again shortly after.Because of the two sudden changes of opposite sign within the time window IV, there are two short-term biases of opposite signs within the time window.They average each other out, which leads to the small bias but large standard deviation of x EKF .Better error estimation might alleviate this feature of our EKF implementation, especially in case we have over-estimated the observation errors and thereby increased the relative weight of the prior estimate.The drift is also partially responsible for the standard deviation for particles smaller than 20 nm in time window II, as the system adjusts to the nucleation event.We remind that our estimates of observation versus background errors are very preliminary.
The impact of sudden changes in the system state tends to be smaller in the overlap region.This is due to the fact that two independent instruments measure the same particle sizes, and provide mutually supporting evidence of the change in the number concentration.This accelerates the adjustment of EKF to the changes.

Measurement noise
Raw measurements contain random noise due to a variety of reasons.It is characteristic to EKF to filter out this noise whereas the inversion of the instrument kernel matrix technique aims to closely fit the solution to the measurements.We illustrate this point as follows.Figure 5 presents raw measurements from DMPS I and DMPS II for three consecutive measurement times (12:00, 12:10, and 12:20 LT).The measurements reveal large number concentration changes in particle sizes of about 10 nm and 100 nm, whereas in other particle sizes, the measurements remain far more stable.This implies that the air mass was broadly the same over this half an hour period.It is impossible to definitely partition these changes to measurement errors, temporary particle emissions, or other effects.In this particular case, x INV closely follows the measurements (not shown) while x EKF is, as expected, far more conservative and smoother.Therefore, the standard deviation of x EKF with respect to raw measurements would be larger than of x INV .
The measurement noise is present in all time windows.In time window II, at about 15:00 to 15:30 LT, the measurements for DMPS I and II in the overlap region disagree significantly due to random noise (not shown).In this case, x INV systematically underestimates both measured particle number concentrations, while x EKF is a compromise between the two measurements.This leads to larger values in the overlap region for x EKF than for x INV , with x EKF being closer to the observed state.There are several possible reasons for the large bias for x INV in the overlap region, but the exact reason is difficult to determine as the focus here was not the features of the inverse solution.Also in time window II, the significant standard deviation in x EKF in particle sizes above 30 nm, and some of the standard deviation in particles smaller than 20 nm, is due to the measurement noise.The same is also true in time window I for particles smaller than 20 nm and at 100 to 200 nm range.
We conclude by noting that our analysis covered about two months of EUCAARI measurements in April-May 2007.A simultaneous analysis of several days is difficult due to the bias being time-sensitive and the differences between the observed and estimated size distribution values depending on the absolute particle number concentration.To provide comparison with the results here, the bias and standard deviation for six days (2, 9, 13, 16, 18 and 28 April 2007) are presented in Fig. 6.These days were chosen because they contained noticeable events in the particle size distributions.Common to all days in Fig. 6 is: (i) in the overlapping measurement range of 10-40 nm the bias for x EKF is smaller than for x INV , (ii) in particle diameters of 50-200 nm the standard deviation for x EKF is noticeably larger than for x INV , and (iii) for particles smaller than 10 nm the statistics for x EKF and x INV are near each other despite the large changes in those particle sizes due to nucleation events.Detailed results presented here for 7 May 2007, though, are representative of the entire period regarding the general behaviour of measurement inversion techniques, including those results shown in Fig. 6.In particular, they high-light the merits and challenges of statistical state estimation in dynamical systems by the EKF technique.

Discussion
The initial results concerning the application of EKF in estimation of aerosol particle size distributions are promising, but also reveal limitations of the method and raise questions about their interpretation.First, the encouraging results and opportunities for EKF are summarized below.i. Enhanced accuracy: despite the approximations made in this implementation of EKF, the results are promising.Especially the treatment of the overlap region of DMPS I and II is very good, and promising regarding the multi-instrument retrievals.Improvements in microphysical modeling, for instance, will directly improve the accuracy of the solution.In the so-called controlled chamber experiments (Sipilä et al., 2010;Brus at al., 2010) the ambient variables and conditions are well known.Measurements from such experiments provide useful test data and opportunities to improve error estimation.
ii. Multi-instrument retrievals: there are no principal obstacles to include new measurements in the EKF inversion, as long as an observation operator and an estimate of observation error variance can be provided.The focus here has been on a single variable, the particle number concentration.However, it is possible in principle to extend EKF to simultaneous estimation of a multi-variate state, for example particle number concentration, particle composition, and ambient vapour concentration.vi.Additional information on the system: EKF provides a minimum-variance state estimate, and also the associated error covariance.This aspect is missing from the inversion of the instrument kernel matrix.In addition as the state estimate is based on the continuity of the system, it is possible to use the innovation term to determine when the system can be considered temporally continuous and when the system experiences large changes in comparison to the previous measurements.
vii.Insensitivity to the initial state: EKF does not appear to be very sensitive for the choice of initial state or the associated error covariance at the very beginning of the filtering.
The challenges concerning the current implementation of EKF are listed below.I. Discontinuities in measurements: in EKF, the background state is a prediction from the previous observation time.The limitation is that the state evolution only includes the dynamical processes of the model.When the system state suddenly changes due to external reasons, there will be a gradual adjustment toward the new state, lagged by the past information.Such situations are related to, for instance, air mass changes, precipitation, or particle influx from external sources.We hope to improve on this in the future by detecting the air mass change either by incorporating information from boundary conditions (e.g.wind and temperature) or by determining a statistically significant change in the aerosol number size distribution.
II. Model microphysics: the time evolution updating benefits of an accurate and universal forward model, in our case the UHMA model.There are numerous ways to describe microphysical processes to different conditions and environments, but some parametrizations are not necessary applicable to all situations.Accuracy of EKF is sensitive to the choices made regarding modelling of different processes, but it is very hard to quantify this sensitivity.
III. Model input data: external quantities (i.e., boundary conditions) influence the evolution of resolved state variables.Regarding particle number size distribution, for instance, the ambient vapour concentrations significantly affect nucleation and condensation processes.In this article, measured ambient vapour concentrations were used as boundary conditions.These are, however, not always available.
IV. Error estimates: EKF relies on estimates of the observation and background error covariances.In this article, the errors were not estimated but specified in a manner which is very likely to be sub-optimal.We intend to improve our uncertainty approximations.
V. Non-Gaussian error: EKF assumes that the associated error is Gaussian.However for very small number concentration values, the associated error cannot be considered Gaussian as the number concentrations cannot be negative.Further information on the error estimates is needed in order to determine when the error can no longer be considered Gaussian.

Conclusions
Extended Kalman Filter (EKF) was introduced to estimate particle number size distributions in a box-model context using observations from a DMPS instrument.Motivation for the research lies in the fact that it is generally hard to estimate size distributions from multiple observations, even of the same measured variable, using an inversion of the instrument kernel matrix technique.Here, a limited EKF implementation was applied to estimate the particle number size distribution by adjusting the time-evolved background state with observations from two DMPSs, which measure on different but partly overlapping particle size ranges.This allowed focus on how the general aspects of aerosol physics impact the EKF implementation.
The self-consistencies of the two solutions, EKF and an inversion of the instrument kernel matrix, were tested by calculating bias and standard deviation for the estimated size distributions with respect to the raw measurements.This was possible by applying the observation operators that are used in EKF to compute the observation counterpart of the model state vector.Despite the assumptions made in the EKF implementation, EKF was found to be more accurate than the inversion of the instrument kernel matrix method in terms of bias, and compatible in terms of standard deviation.The analysis covered about two months of EUCAARI measurements in April-May 2007.The detailed results were presented for 7 May 2007 which was selected as a representative example of the entire period regarding the overall behaviour of these inversion techniques.
Generally, the limited EKF implementation was found to be satisfactory, and justifies the more extensive EKF implementation examined in Part 2 of this work.Potential further improvements of the EKF implementation include more accurate estimation of error covariance of the measurements and the background state.

T.
Viskari et al.: Estimation of aerosol particle number distributions with Kalman Filtering causes, such as precipitation, air mass change or external influx of particles, are not included in the UHMA.

Fig. 1 .
Fig. 1.The ambient vapour concentrations for sulphuric acid (blue), non-volatile organic compound (red) and volatile organic compound (green) as applied in the EKF implementation for 7 May 2007.

Fig. 2b .
Fig. 2b.As in (a), but for inversion of the instrument kernel matrix (x INV ).

Fig. 2c .
Fig. 2c.The total number concentrations for particles larger than 3 nm for x EKF (solid; blue) and x INV (dashed; green).

Fig. 3a .
Fig. 3a.The bias of x EKF (blue) and x INV (red) as compared to the raw measurements from DMPS I (solid) and II (dashed) on 7 May 2007.The bias is on the y-axis [cm −3 ] and particle diameter on the x-axis [m].Note that the results are presented in the characteristic diameters.

Fig. 3b .
Fig. 3b.As in (a), but for the standard deviation.

Fig. 4 .
Fig. 4. The bias and standard deviation of x EKF (blue) and x INV (red) as compared to the raw measurements from DMPS I (solid) and II (dashes) for time windows I to IV.The bias and standard deviation are on the y-axis [cm −3 ] and particle diameter on the x-axisis [m].Note that the results are presented for the characteristic diameters.

Fig. 6 .
Fig. 6.The bias and standard deviation of x EKF (blue) and x INV (red) as compared to the raw measurements from DMPS I (solid) and II (dashes) for 6 days in April.The bias and standard deviation are on the y-axis [cm −3 ] and particle diameter on the x-axis [m].Note that the results are presented for the characteristic diameters.