# Application of Multistep Inversion Method for Online Monitoring Aerosol Particle Size Distribution and Aerosol Concentration

Authors: Zhenzong He; Liang Xu; Junkui Mao; Xingsi Han; Biao Zhang
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

Aerosol concentration in the flow is usually time varying, and aerosol particle size distribution (PSD) is considered to be unchanged, which increases the difficulty of the measurement of aerosol PSD and concentration online. To solve these problems, a kind of multistep inversion method based on the angular light-scattering (ALS) signals is proposed. First, the aerosol PSD is estimated using shuffled frog-leaping algorithms (SFLAs) from relative ALS signals. Then, with aerosol PSD as priori information, the aerosol concentration is obtained by the Kalman filter (KF) algorithm, widely used in the real-time control system of industrial facilities for its ability of fast predictions. The result reveals that the performance of the improved SFLA is better than that of the original SFLA in solving the aerosol PSD. Moreover, in studying the aerosol concentration, more accurate results can be obtained with larger standard deviation of process noise or smaller standard deviation of measurement noise, while decreasing sampling time interval can improve the accuracy of retrieval results and reduce time delay to a certain degree. So, to improve retrieval accuracy, the noise should be controlled, and appropriate sampling time interval should be selected. All the numerical simulations confirm that the methodology provides effective and reliable results in real-time estimating.


---

## Body

## 1. Introduction

Aerosols usually play a pivotal role in determining the properties of the atmosphere, e.g., by reducing visibility or affecting the net radiative fluxes and temperature [1–3]. Considerable researches have been carried out on studying the aerosol particle size, optical constant, and concentration in the particle dispersed medium, e.g., Aerosol Robotic Network (AERONET) and Moderate Resolution Imaging Spectroradiometer (MODIS), which are the global ground-based aerosol observation networks established to study the properties of the atmospheric aerosols [4, 5]. Usually, in the atmosphere, the aerosol affects the radiative transfer, meteorological phenomena, and climate trends. The aerosol particle size distribution (PSD) and concentration are regarded as the important indicators of industrial emissions. The aerosol PSDs are usually treated as constants, while the aerosol concentration is usually changing over time, which increases the difficulty of the online monitoring of the PSDs and time-varying particle concentration [6–9].The online particle monitoring has been widely studied [10–12], and the common measurement techniques for studying the properties of the particles contain optical measurement methods, acoustic measurement methods, and electric induction methods [10, 13]. In the past few years, the optical measurement method (e.g., spectral extinction method, angular light-scattering method, diffraction light-scattering method, dynamic light-scattering method, laser-induced incandescence method, and laser-induced breakdown spectroscopy method) has been widely utilized to study the properties of particle system owing to offering a useful and effective approach to characterize a large number of industrial production processes. Their measurements vary from nanometer to millimeter [14]. Among these optical measurement methods, the angular light-scattering (ALS) method is regarded as one of the most viable. For example, the ALS method was applied to study the fractal dimension and size distribution of noncompact soot aggregates [15], the optical constants of the particle [16], the particle size distribution [17, 18], and the particle concentration [19–22]. In the author’s previous work [23, 24], the ALS method combined with inverse radiative problem technique is successfully applied to study the PSD and optical constants. Moreover, the convergence accuracy and robustness of the ALS method was proved to be better than that of the spectral extinction method.In the present study, to make the problem mathematically trackable, the aerosol particle is assumed to be spherical, and the optical properties of aerosol particles are simulated by the Mie theory. The ALS method is introduced to study the aerosol PSD and time-varying aerosol concentration. Moreover, considering the properties of the aerosol PSD and aerosol concentration mentioned above, a kind of multistep inversion method is proposed. First, the aerosol PSD is retrieved by shuffled frog-leaping algorithms (SFLAs) [25] from the relative ALS signals, which is independent of the aerosol concentration. Then, the solved aerosol PSD is regarded as known quantity, and the aerosol concentration is monitored by the Kalman filter (KF) algorithm, which performs very well in real-time control systems of industrial facilities for its ability of fast predictions, such as the particle monitoring [26, 27], species monitoring in combustion applications [28, 29], and others [26, 28, 30–32]. Finally, main conclusions and prospects for further research are provided.


## 2. Methodology

### 2.1. Angular Light-Scattering (ALS) Method

When a collimated light with intensityI0 impinges on a particle system, the transmitted light is scattered and absorbed by the particles. If the aerosol is optically thin and the independent scattering dominates, the real-time angular light-scattering (Figure 1) intensity I (θ, t) at the angle θ can be expressed as follows [33–36]:(1)Iθ,t=πNDtλ∫DminDmaxiθ,mλ,D,λfDdD,where D is the diameter of the particle; Dmin and Dmax denote the lower and upper integration limits of particle size; I (θ, t) denotes the angular light-scattering intensity at angle θ and time t, which can be measured by the optical sensor; ND (t) is the unknown number concentration of suspended particle system which is a function of measurement time t; f (D) is the unknown particle size distribution (PSD) which is usually regarded as unchanged in the particle system; and i (θ, mλ, D, λ) denotes the light-scattering intensity of particles with diameter D at the angle θ from the incident direction, which can be derived in terms of the Mie scattering functions i1 (θ, mλ, D, λ) and i2 (θ, mλ, D, λ) [18, 37, 38]:(2)iθ,mλ,D,λ=i1θ,mλ,D,λ+i2θ,mλ,D,λ8π2λ2I0,where mλ, an optical constant, refers to the interaction between particles and incident light and λ denotes the incident wavelength (μm).Figure 1
Geometry and physical model of the angular light-scattering method.


### 2.2. Shuffled Frog-Leaping Algorithms (SFLAs)

#### 2.2.1. Original SFLA

The original SFLA, based on the evolution of memes carried by interactive individuals and global exchange of information among the interactive population, was developed by Eusuff and Lansey [39] to seek a global optimal solution by performing an informed heuristic search using a heuristic function. The concept of memetic comes from the word “meme” which can be considered as the unit of culture evolution and the sample of virtual frogs which stands for a possible solution constitutes a population.First, the total number of the frogsS in the swarm and the number of unsolved variables N are fixed and the initial virtual frog population PX = [X1, X2, …, XS] is generated randomly in the original SFLA. Each frog is represented by a 1-D vector Xi = [xi1, xi2, …, xiN], which denotes a potential solution of the inverse problem. The frogs are evaluated and sorted in a descending order according to the value of fitness function, and the frog with the global best fitness is recorded and identified as Xg. Then, the entire swarm is divided into m memeplexes each holding n frogs such that N = m × n. The memeplexes can be regarded as the parallel frog cultures for the same goal. Each frog culture proceeds towards their goal independently. The kth memeplexes Mk are defined as [25, 39, 40](3)Mk=Xk+ml−1∈PXl=1,2,…,n,k=1,2,…,m.The frogs with the best and worst fitness values of each memeplex are identified asXb and Xw. After the previous preparation, the evolution process is conducted to search the solution of the inverse problem. However, not all the frogs change their positions in the searching process, but only the ones with the worst fitness in each memeplex. The corresponding positions of the worst frogs are adjusted as follows:(4)d=Rand1⋅Xb−Xw,(5)Xw′=Xw+d,where d = [d1, d2, … , dN] with dmin≤dj≤dmax, dmax and dmin denote the maximum and minimum allowed limits for a frog’s position, respectively; Rand1 is the uniformly distributed random numbers in the range of 0,1; and Xw′ denotes the new position of the worst frog. The new generated position will be valuated after updating, and the worst frog will be replaced if the fitness function value of the new position is superior to that of the original one. Otherwise, the execution of equations (4) and (5) is repeated with respect to the global best frog. In the calculation of the step size d, the best frog in this memeplex Xb is replaced by the global best frog Xg. This operation is going to be repeated for Emax times before taking the next step. A new random frog will be generated to replace the worst frog if there is no improvement of the new frog. All the frogs are mixed together and redivided into m memeplexes when this process is done. The information of each memeplex is passed between each other through the shuffling process. The whole process is repeated unless the iteration number reaches the user-defined iteration limit, or the fitness function value of global best position is less than the tolerance. The details of the SFLA are available in Ref. [25].


#### 2.2.2. Improved SFLA

In the SFLA, the initial positions are generated by random numbers. In this way, some of the initial positions are too close to each other. To solve these problems, the Chaos theory, the highly unstable motion of deterministic systems in finite phase space often existing in nonlinear systems, is applied to initialize the positions of the frogs. The Chaos theory studies the behavior of dynamical systems that are highly sensitive to initial conditions, an effect that is popularly referred to as the butterfly effect [41]. The typical logistic mapping to generate the Chaos signal is described as follows [41]:(6)r0=Rand2,rk+1=μrk1−rk,k=0,1,…,S−1,where Rand2 is a uniformly distributed random number in the range of [0, 1] and μ is the control number; when μ = 4.0, the logistic mapping is in a fully chaotic state.Moreover, although each frog culture evolves a certain iteration without the exchange of information with other cultures to improve the local searching ability, the information exchange within each culture only proceeds between the best and worst frogs which will restrict the convergence speed and the local search ability. Each memeplex is iterated plenty of times caused by the low convergence, which occupies large amount of calculation time. To remedy these problems, the particle swarm optimization (PSO) algorithm is introduced to improve the SFLA by accelerating the convergence speed and local searching ability in each memeplexes. The basic principle of the algorithm can be explained as follows. After local search process, i.e., updating the worst frog in each memeplex, the whole frog population is shuffled and divided into new memeplexes. Subsequently, the PSO algorithm is introduced to improve the local search process in each memeplex. The velocity and position evolution equations of the worst frog in each memeplex are expressed as [42](7)Vw′=C1⋅Rand3⋅Vw+C2⋅Rand4⋅Xg−Xw′,(8)Xw″=Xw′+Vw′,where Rand3 and Rand4 are two uniformly distributed random numbers in the range of [0, 1] and C1 and C2 are two positive constants called acceleration coefficients. The pseudocode for the improved SFLA can be summarized as follows:(i)
Step 1: input the system control parameters of the improved SFLA, i.e., the total number of frogs in the swarmS, the numbers of memeplexes m, the number of frogs in each memeplex n, the number of evolutionary iterations Emax, the maximum iteration limit Nc, the number of unsolved variables N, the searching space [lowi, highi] of each variable, and the tolerance for minimizing the fitness function value ε.(ii)
Step 2: initialize the frogs’ positions by mapping the chaotic sequence, which is generated by equation (6), to the search space.(iii)
Step 3: evaluate the frogs and calculate the corresponding fitness, and then rank them in the descending order.(iv)
Step 4: divide the frogs intom memeplexes according to equation (3).(v)
Step 5: generate the new position of the worst frog in each memeplex by using equations (7) and (8). If the position of the new frog is better, then use it to update the worst frog. Otherwise, repeat this step by replacing Xb with Xg.(vi)
Step 6: evaluate the new generated frog. If the position of the new frog is superior, then use it to update the worst frog and go to Step 7. Otherwise go to Step 5 until the maximum evolution iterationsEmax is reached.(vii)
Step 7: go to Step 4 until one of the stop criteria is reached. The stop criteria are as follows:(1)
The iteration number reaches the user-defined iteration limitNc, iter (t) > Nc(2)
The fitness function value of global best position is less than the toleranceε, Fitness (Xg) < ε


### 2.3. Kalman Filter Algorithm (KF)

In this study, the unknown time-varying particle number concentration in a particle dispersed system is predicted by measuring real-time ALS intensity, and the equation (1) can be transformed to a linear system and described as(9)Iθ,tI0=B⋅NDt,B=λ8π∫DminDmaxi1θ,mλ,D,λ+i2θ,mλ,D,λfDdD.According to the Kalman filter algorithm, when there is white Gaussian noise added to the measurement signals, the linear system can be expressed by the following state and measurement equations [43, 44]:(10)Xk=Φk−1Xk−1+wk−1,Yk=HkXk+vk,where Xk,Xk−1∈Rn denotes the n-dimensional state vector at instant k and k − 1, respectively, and Xk=lnτtk in the present study; Yk∈Rm represents the m-dimensional observation vector at instant k; Φk−1 represents the state transition matrix at instant k − 1; H (k) stands for the observation matrix at instant k; and w (k − 1) and v (k) denote the mutually uncorrelated process noise at instant k − 1 and measurement noise at instant k, respectively, with zero mean and covariance matrix written as [43, 44](11)Ewk=0,EwkwTj=Qk=diagσQ2,Evk=0,EvkvTj=Rk=diagσR2,where Q and R denote the covariance matrices of the process noise and the measurement noise, respectively; σQ and σR are the standard deviations of the process noise and measurement noise, respectively. To simplify the problem, the measurement noise is assumed to be additive noise and obey the standard Gaussian distribution. Usually, the Kalman filter algorithm contains two update steps. The time update mainly can be divided into state prediction and covariance prediction, while the measurement update is composed of gain calculation, state correction, and covariance correction. The time update equations and measurement update equations of the Kalman filter algorithm can be written as follows [44, 45]: (i) Predicting the present stage on the basis of the previous moment stage(12)X¯kk−1=Φk−1X¯k−1,Pkk−1=Φk−1Pk−1ΦTk−1+Qk−1.(ii) Updating the measurement stage from the corrected prediction stage based on the present measurement signals(13)sk=HkPkk−1HTk+Rk,Kk=Pkk−1HTks−1k,Y¯k=Yk−HkX¯kk−1,X¯k=X¯kk−1+KkY¯k,N¯Dtk=X¯kB.Moreover, the posterior covariance is updated as follows:(14)Pk=I−KkHPkk−1,where P denotes the error covariance matrix of the state estimate vector. The flowchart of numerical procedures for retrieving the real-time particle number concentration using the KF algorithm is shown in Figure 2. First, the exact value of particle concentration is used to calculate the ALS signals, which is treated as the measurement signals in the simulation experiment without noise. Then, the measurement noise and process noise are added to the measurement signals to simulate observation measurement signals with noise. Finally, the KF algorithm is applied to obtain the retrieval value of particle concentration from the obversion measurement signals with noise.Figure 2
Flowchart of the numerical procedure for studying the real-time particle number concentration using the KF algorithm.


### 2.4. Inverse Process

According to the theory of the inverse radiative problems [46], the inversion of multiple parameters usually requires the inverse model to have better convergence characteristics to avoid the multivalue characteristics of the retrieval results and make sure each parameter with good inverse accuracy. So, when the inverse model and algorithm are fixed, the dimension of the unsolved parameters should be reduced as much as possible, and the inversion strategies are very important.To solve the problem, a kind of multistep inversion strategy is proposed in this study. The strategy is to perform the inversion studies on some easily accessible parameters first. Then, these parameters are regarded as known quantity to reconstruct other unsolved parameters. Based on this idea, the aerosol PSDs are estimated from the relative ALS signals at first, and the mathematical expression ofIR (θ) is derived as(15)IRθ=Iθ,tIθ0,t=∫DminDmaxiθ,mλ,D,λfDdD∫DminDmaxiθ0,mλ,D,λfDdD,where θ0 denotes the angle of the reference signals which can be found in Figure 1. From equation (15), it can be found that the relative ALS signals are not involved with the number concentration of suspended particle system. So, the aerosol PSDs can be reconstructed without considering the effect of the time-varying particle number concentration. Subsequently, the time-varying particle number concentration ND(t) is estimated from the real-time light-scattering intensity with the aerosol PSDs as the priori information. In this study, the aerosol PSDs are estimated by the SFLA and improved SFLA from relative ALS signals, respectively, and the time-varying particle number concentration ND(t) is estimated by the KF algorithm from the real-time light-scattering intensity. The logical relationship of numerical procedures for online monitoring the aerosol PSDs and the real-time particle number concentration is shown in Figure 3.Figure 3
Logical relationship of multistep inversion method for online monitoring the aerosol PSDs and real-time particle number concentration.


## 2.1. Angular Light-Scattering (ALS) Method

When a collimated light with intensityI0 impinges on a particle system, the transmitted light is scattered and absorbed by the particles. If the aerosol is optically thin and the independent scattering dominates, the real-time angular light-scattering (Figure 1) intensity I (θ, t) at the angle θ can be expressed as follows [33–36]:(1)Iθ,t=πNDtλ∫DminDmaxiθ,mλ,D,λfDdD,where D is the diameter of the particle; Dmin and Dmax denote the lower and upper integration limits of particle size; I (θ, t) denotes the angular light-scattering intensity at angle θ and time t, which can be measured by the optical sensor; ND (t) is the unknown number concentration of suspended particle system which is a function of measurement time t; f (D) is the unknown particle size distribution (PSD) which is usually regarded as unchanged in the particle system; and i (θ, mλ, D, λ) denotes the light-scattering intensity of particles with diameter D at the angle θ from the incident direction, which can be derived in terms of the Mie scattering functions i1 (θ, mλ, D, λ) and i2 (θ, mλ, D, λ) [18, 37, 38]:(2)iθ,mλ,D,λ=i1θ,mλ,D,λ+i2θ,mλ,D,λ8π2λ2I0,where mλ, an optical constant, refers to the interaction between particles and incident light and λ denotes the incident wavelength (μm).Figure 1
Geometry and physical model of the angular light-scattering method.


## 2.2. Shuffled Frog-Leaping Algorithms (SFLAs)

### 2.2.1. Original SFLA

The original SFLA, based on the evolution of memes carried by interactive individuals and global exchange of information among the interactive population, was developed by Eusuff and Lansey [39] to seek a global optimal solution by performing an informed heuristic search using a heuristic function. The concept of memetic comes from the word “meme” which can be considered as the unit of culture evolution and the sample of virtual frogs which stands for a possible solution constitutes a population.First, the total number of the frogsS in the swarm and the number of unsolved variables N are fixed and the initial virtual frog population PX = [X1, X2, …, XS] is generated randomly in the original SFLA. Each frog is represented by a 1-D vector Xi = [xi1, xi2, …, xiN], which denotes a potential solution of the inverse problem. The frogs are evaluated and sorted in a descending order according to the value of fitness function, and the frog with the global best fitness is recorded and identified as Xg. Then, the entire swarm is divided into m memeplexes each holding n frogs such that N = m × n. The memeplexes can be regarded as the parallel frog cultures for the same goal. Each frog culture proceeds towards their goal independently. The kth memeplexes Mk are defined as [25, 39, 40](3)Mk=Xk+ml−1∈PXl=1,2,…,n,k=1,2,…,m.The frogs with the best and worst fitness values of each memeplex are identified asXb and Xw. After the previous preparation, the evolution process is conducted to search the solution of the inverse problem. However, not all the frogs change their positions in the searching process, but only the ones with the worst fitness in each memeplex. The corresponding positions of the worst frogs are adjusted as follows:(4)d=Rand1⋅Xb−Xw,(5)Xw′=Xw+d,where d = [d1, d2, … , dN] with dmin≤dj≤dmax, dmax and dmin denote the maximum and minimum allowed limits for a frog’s position, respectively; Rand1 is the uniformly distributed random numbers in the range of 0,1; and Xw′ denotes the new position of the worst frog. The new generated position will be valuated after updating, and the worst frog will be replaced if the fitness function value of the new position is superior to that of the original one. Otherwise, the execution of equations (4) and (5) is repeated with respect to the global best frog. In the calculation of the step size d, the best frog in this memeplex Xb is replaced by the global best frog Xg. This operation is going to be repeated for Emax times before taking the next step. A new random frog will be generated to replace the worst frog if there is no improvement of the new frog. All the frogs are mixed together and redivided into m memeplexes when this process is done. The information of each memeplex is passed between each other through the shuffling process. The whole process is repeated unless the iteration number reaches the user-defined iteration limit, or the fitness function value of global best position is less than the tolerance. The details of the SFLA are available in Ref. [25].


### 2.2.2. Improved SFLA

In the SFLA, the initial positions are generated by random numbers. In this way, some of the initial positions are too close to each other. To solve these problems, the Chaos theory, the highly unstable motion of deterministic systems in finite phase space often existing in nonlinear systems, is applied to initialize the positions of the frogs. The Chaos theory studies the behavior of dynamical systems that are highly sensitive to initial conditions, an effect that is popularly referred to as the butterfly effect [41]. The typical logistic mapping to generate the Chaos signal is described as follows [41]:(6)r0=Rand2,rk+1=μrk1−rk,k=0,1,…,S−1,where Rand2 is a uniformly distributed random number in the range of [0, 1] and μ is the control number; when μ = 4.0, the logistic mapping is in a fully chaotic state.Moreover, although each frog culture evolves a certain iteration without the exchange of information with other cultures to improve the local searching ability, the information exchange within each culture only proceeds between the best and worst frogs which will restrict the convergence speed and the local search ability. Each memeplex is iterated plenty of times caused by the low convergence, which occupies large amount of calculation time. To remedy these problems, the particle swarm optimization (PSO) algorithm is introduced to improve the SFLA by accelerating the convergence speed and local searching ability in each memeplexes. The basic principle of the algorithm can be explained as follows. After local search process, i.e., updating the worst frog in each memeplex, the whole frog population is shuffled and divided into new memeplexes. Subsequently, the PSO algorithm is introduced to improve the local search process in each memeplex. The velocity and position evolution equations of the worst frog in each memeplex are expressed as [42](7)Vw′=C1⋅Rand3⋅Vw+C2⋅Rand4⋅Xg−Xw′,(8)Xw″=Xw′+Vw′,where Rand3 and Rand4 are two uniformly distributed random numbers in the range of [0, 1] and C1 and C2 are two positive constants called acceleration coefficients. The pseudocode for the improved SFLA can be summarized as follows:(i)
Step 1: input the system control parameters of the improved SFLA, i.e., the total number of frogs in the swarmS, the numbers of memeplexes m, the number of frogs in each memeplex n, the number of evolutionary iterations Emax, the maximum iteration limit Nc, the number of unsolved variables N, the searching space [lowi, highi] of each variable, and the tolerance for minimizing the fitness function value ε.(ii)
Step 2: initialize the frogs’ positions by mapping the chaotic sequence, which is generated by equation (6), to the search space.(iii)
Step 3: evaluate the frogs and calculate the corresponding fitness, and then rank them in the descending order.(iv)
Step 4: divide the frogs intom memeplexes according to equation (3).(v)
Step 5: generate the new position of the worst frog in each memeplex by using equations (7) and (8). If the position of the new frog is better, then use it to update the worst frog. Otherwise, repeat this step by replacing Xb with Xg.(vi)
Step 6: evaluate the new generated frog. If the position of the new frog is superior, then use it to update the worst frog and go to Step 7. Otherwise go to Step 5 until the maximum evolution iterationsEmax is reached.(vii)
Step 7: go to Step 4 until one of the stop criteria is reached. The stop criteria are as follows:(1)
The iteration number reaches the user-defined iteration limitNc, iter (t) > Nc(2)
The fitness function value of global best position is less than the toleranceε, Fitness (Xg) < ε


## 2.2.1. Original SFLA

The original SFLA, based on the evolution of memes carried by interactive individuals and global exchange of information among the interactive population, was developed by Eusuff and Lansey [39] to seek a global optimal solution by performing an informed heuristic search using a heuristic function. The concept of memetic comes from the word “meme” which can be considered as the unit of culture evolution and the sample of virtual frogs which stands for a possible solution constitutes a population.First, the total number of the frogsS in the swarm and the number of unsolved variables N are fixed and the initial virtual frog population PX = [X1, X2, …, XS] is generated randomly in the original SFLA. Each frog is represented by a 1-D vector Xi = [xi1, xi2, …, xiN], which denotes a potential solution of the inverse problem. The frogs are evaluated and sorted in a descending order according to the value of fitness function, and the frog with the global best fitness is recorded and identified as Xg. Then, the entire swarm is divided into m memeplexes each holding n frogs such that N = m × n. The memeplexes can be regarded as the parallel frog cultures for the same goal. Each frog culture proceeds towards their goal independently. The kth memeplexes Mk are defined as [25, 39, 40](3)Mk=Xk+ml−1∈PXl=1,2,…,n,k=1,2,…,m.The frogs with the best and worst fitness values of each memeplex are identified asXb and Xw. After the previous preparation, the evolution process is conducted to search the solution of the inverse problem. However, not all the frogs change their positions in the searching process, but only the ones with the worst fitness in each memeplex. The corresponding positions of the worst frogs are adjusted as follows:(4)d=Rand1⋅Xb−Xw,(5)Xw′=Xw+d,where d = [d1, d2, … , dN] with dmin≤dj≤dmax, dmax and dmin denote the maximum and minimum allowed limits for a frog’s position, respectively; Rand1 is the uniformly distributed random numbers in the range of 0,1; and Xw′ denotes the new position of the worst frog. The new generated position will be valuated after updating, and the worst frog will be replaced if the fitness function value of the new position is superior to that of the original one. Otherwise, the execution of equations (4) and (5) is repeated with respect to the global best frog. In the calculation of the step size d, the best frog in this memeplex Xb is replaced by the global best frog Xg. This operation is going to be repeated for Emax times before taking the next step. A new random frog will be generated to replace the worst frog if there is no improvement of the new frog. All the frogs are mixed together and redivided into m memeplexes when this process is done. The information of each memeplex is passed between each other through the shuffling process. The whole process is repeated unless the iteration number reaches the user-defined iteration limit, or the fitness function value of global best position is less than the tolerance. The details of the SFLA are available in Ref. [25].


## 2.2.2. Improved SFLA

In the SFLA, the initial positions are generated by random numbers. In this way, some of the initial positions are too close to each other. To solve these problems, the Chaos theory, the highly unstable motion of deterministic systems in finite phase space often existing in nonlinear systems, is applied to initialize the positions of the frogs. The Chaos theory studies the behavior of dynamical systems that are highly sensitive to initial conditions, an effect that is popularly referred to as the butterfly effect [41]. The typical logistic mapping to generate the Chaos signal is described as follows [41]:(6)r0=Rand2,rk+1=μrk1−rk,k=0,1,…,S−1,where Rand2 is a uniformly distributed random number in the range of [0, 1] and μ is the control number; when μ = 4.0, the logistic mapping is in a fully chaotic state.Moreover, although each frog culture evolves a certain iteration without the exchange of information with other cultures to improve the local searching ability, the information exchange within each culture only proceeds between the best and worst frogs which will restrict the convergence speed and the local search ability. Each memeplex is iterated plenty of times caused by the low convergence, which occupies large amount of calculation time. To remedy these problems, the particle swarm optimization (PSO) algorithm is introduced to improve the SFLA by accelerating the convergence speed and local searching ability in each memeplexes. The basic principle of the algorithm can be explained as follows. After local search process, i.e., updating the worst frog in each memeplex, the whole frog population is shuffled and divided into new memeplexes. Subsequently, the PSO algorithm is introduced to improve the local search process in each memeplex. The velocity and position evolution equations of the worst frog in each memeplex are expressed as [42](7)Vw′=C1⋅Rand3⋅Vw+C2⋅Rand4⋅Xg−Xw′,(8)Xw″=Xw′+Vw′,where Rand3 and Rand4 are two uniformly distributed random numbers in the range of [0, 1] and C1 and C2 are two positive constants called acceleration coefficients. The pseudocode for the improved SFLA can be summarized as follows:(i)
Step 1: input the system control parameters of the improved SFLA, i.e., the total number of frogs in the swarmS, the numbers of memeplexes m, the number of frogs in each memeplex n, the number of evolutionary iterations Emax, the maximum iteration limit Nc, the number of unsolved variables N, the searching space [lowi, highi] of each variable, and the tolerance for minimizing the fitness function value ε.(ii)
Step 2: initialize the frogs’ positions by mapping the chaotic sequence, which is generated by equation (6), to the search space.(iii)
Step 3: evaluate the frogs and calculate the corresponding fitness, and then rank them in the descending order.(iv)
Step 4: divide the frogs intom memeplexes according to equation (3).(v)
Step 5: generate the new position of the worst frog in each memeplex by using equations (7) and (8). If the position of the new frog is better, then use it to update the worst frog. Otherwise, repeat this step by replacing Xb with Xg.(vi)
Step 6: evaluate the new generated frog. If the position of the new frog is superior, then use it to update the worst frog and go to Step 7. Otherwise go to Step 5 until the maximum evolution iterationsEmax is reached.(vii)
Step 7: go to Step 4 until one of the stop criteria is reached. The stop criteria are as follows:(1)
The iteration number reaches the user-defined iteration limitNc, iter (t) > Nc(2)
The fitness function value of global best position is less than the toleranceε, Fitness (Xg) < ε


## 2.3. Kalman Filter Algorithm (KF)

In this study, the unknown time-varying particle number concentration in a particle dispersed system is predicted by measuring real-time ALS intensity, and the equation (1) can be transformed to a linear system and described as(9)Iθ,tI0=B⋅NDt,B=λ8π∫DminDmaxi1θ,mλ,D,λ+i2θ,mλ,D,λfDdD.According to the Kalman filter algorithm, when there is white Gaussian noise added to the measurement signals, the linear system can be expressed by the following state and measurement equations [43, 44]:(10)Xk=Φk−1Xk−1+wk−1,Yk=HkXk+vk,where Xk,Xk−1∈Rn denotes the n-dimensional state vector at instant k and k − 1, respectively, and Xk=lnτtk in the present study; Yk∈Rm represents the m-dimensional observation vector at instant k; Φk−1 represents the state transition matrix at instant k − 1; H (k) stands for the observation matrix at instant k; and w (k − 1) and v (k) denote the mutually uncorrelated process noise at instant k − 1 and measurement noise at instant k, respectively, with zero mean and covariance matrix written as [43, 44](11)Ewk=0,EwkwTj=Qk=diagσQ2,Evk=0,EvkvTj=Rk=diagσR2,where Q and R denote the covariance matrices of the process noise and the measurement noise, respectively; σQ and σR are the standard deviations of the process noise and measurement noise, respectively. To simplify the problem, the measurement noise is assumed to be additive noise and obey the standard Gaussian distribution. Usually, the Kalman filter algorithm contains two update steps. The time update mainly can be divided into state prediction and covariance prediction, while the measurement update is composed of gain calculation, state correction, and covariance correction. The time update equations and measurement update equations of the Kalman filter algorithm can be written as follows [44, 45]: (i) Predicting the present stage on the basis of the previous moment stage(12)X¯kk−1=Φk−1X¯k−1,Pkk−1=Φk−1Pk−1ΦTk−1+Qk−1.(ii) Updating the measurement stage from the corrected prediction stage based on the present measurement signals(13)sk=HkPkk−1HTk+Rk,Kk=Pkk−1HTks−1k,Y¯k=Yk−HkX¯kk−1,X¯k=X¯kk−1+KkY¯k,N¯Dtk=X¯kB.Moreover, the posterior covariance is updated as follows:(14)Pk=I−KkHPkk−1,where P denotes the error covariance matrix of the state estimate vector. The flowchart of numerical procedures for retrieving the real-time particle number concentration using the KF algorithm is shown in Figure 2. First, the exact value of particle concentration is used to calculate the ALS signals, which is treated as the measurement signals in the simulation experiment without noise. Then, the measurement noise and process noise are added to the measurement signals to simulate observation measurement signals with noise. Finally, the KF algorithm is applied to obtain the retrieval value of particle concentration from the obversion measurement signals with noise.Figure 2
Flowchart of the numerical procedure for studying the real-time particle number concentration using the KF algorithm.


## 2.4. Inverse Process

According to the theory of the inverse radiative problems [46], the inversion of multiple parameters usually requires the inverse model to have better convergence characteristics to avoid the multivalue characteristics of the retrieval results and make sure each parameter with good inverse accuracy. So, when the inverse model and algorithm are fixed, the dimension of the unsolved parameters should be reduced as much as possible, and the inversion strategies are very important.To solve the problem, a kind of multistep inversion strategy is proposed in this study. The strategy is to perform the inversion studies on some easily accessible parameters first. Then, these parameters are regarded as known quantity to reconstruct other unsolved parameters. Based on this idea, the aerosol PSDs are estimated from the relative ALS signals at first, and the mathematical expression ofIR (θ) is derived as(15)IRθ=Iθ,tIθ0,t=∫DminDmaxiθ,mλ,D,λfDdD∫DminDmaxiθ0,mλ,D,λfDdD,where θ0 denotes the angle of the reference signals which can be found in Figure 1. From equation (15), it can be found that the relative ALS signals are not involved with the number concentration of suspended particle system. So, the aerosol PSDs can be reconstructed without considering the effect of the time-varying particle number concentration. Subsequently, the time-varying particle number concentration ND(t) is estimated from the real-time light-scattering intensity with the aerosol PSDs as the priori information. In this study, the aerosol PSDs are estimated by the SFLA and improved SFLA from relative ALS signals, respectively, and the time-varying particle number concentration ND(t) is estimated by the KF algorithm from the real-time light-scattering intensity. The logical relationship of numerical procedures for online monitoring the aerosol PSDs and the real-time particle number concentration is shown in Figure 3.Figure 3
Logical relationship of multistep inversion method for online monitoring the aerosol PSDs and real-time particle number concentration.


## 3. Numerical Simulation

According to our previous work [23], the optical constant of the aerosol particle in the present study is set as mλ = 1.53 + 0.01i, and the measurement wavelength is set as λ = 0.5 μm. All the simulations were performed on an Intel Core i7–6500 PC by using Matlab 2016. According to Ref. [46], the particle size range is set as [0.1 μm, 10 μm]. To make the problem mathematically trackable, the particle size distribution was assumed to be constant and obeys the common log-normal (L-N) and gamma function, respectively, and the corresponding mathematical expressions are expressed as(16)fL−ND=12πDlnσ×exp−lnD−lnD¯22lnσ2,fGammaD=Dα×exp−βDγ,where D¯ is the average diameter of L-N distribution; σ denotes the width of the distribution; and α, β, and γ are the characteristic parameters of the gamma distribution. Usually, in the modified form, γ = 1, so only parameters α and β in the gamma distribution need to be investigated, respectively. The retrieval of the aerosol PSDs is solved by minimizing the fitness function value Fitness, which is defined as the sum of the square residual between the estimated and measurement signals ratios. The mathematic expression of Fitness is derived as(17)Fitness=∑i=1NθIRθiest−IRθimeaIRθimea2,where Nθ is the number of the measurement angles. For the inverse algorithms are the stochastic optimization method and all optimizations have certain randomness, all the calculations are repeated 20 times. To ensure the online monitoring, the relative ALS signals for repeated calculations are obtained at 20 different measurement moments. The PSD is assumed to be unchanged during different measurement moments, and the average value of the 20 times retrieval results is regarded as the final retrieval results. For the purpose of investigating the reliability and feasibility of the inverse algorithms, the relative standard deviation ξ, which means the sum of the deviation between the probability distribution estimated from the inverse calculation and the true distribution of aerosol PSDs in every subinterval, is studied to evaluate the quality of inverse results, and its mathematic expression is described as(18)ξ=∑i=1N′fX,estD˜i−fX,trueD˜i21/2∑i=1N′fX,trueD˜i21/2,X=L−N or Gamma,where N′ denotes the number of subintervals which the particle size range [Dmin, Dmax] is divided into; D¯i is the midpoint of the ith subinterval [Di, Di+1]; fX,trueD¯i is the true aerosol PSDs in the ith subinterval; and fX,estD¯i is the estimated aerosol PSDs in the ith subinterval.The time-varying particle number concentrationND(t) in this study is assumed to obey two common time-varying distributions, e.g., sine wave-type equation (19) and linear-type equation (20), and the corresponding mathematical expressions are as follows:(19)NDt=5+sin0.2t,(20)NDt=0.4t−3,10≤t≤25,−0.4t+17,25≤t≤40,1,others.To investigate the reliability and feasibility of the Kalman filter algorithm, the effects of measurement noise and process noise on the estimated results are studied. The average deviation δ, which means the average value of the deviation between estimated time-varying particle number concentrationND,estti and the true time-varying particle number concentration ND,trueti in each sampling time interval, is used to evaluate the estimated results, and its mathematic expression is described as(21)δ=1Nt∑i=1NtND,estti−ND,trueti21/2,where T is the total sample time and Nt denotes the number of sampling time intervals, Nt=T/Δt. Moreover, to make the problem mathematically trackable, the distribution functions of the aerosol PSDs and aerosol concentration are assumed to be known beforehand.


### 3.1. Retrieval of Aerosol Particle Size Distributions

With the help of the original SFLA and improved SFLA, the aerosol PSDs are estimated from the relative ALS signals, and the real values are shown in Table1. According to our previous work [24], the measurement angles are set as θ = 5°, 10°, 15°, and 20°, and the reference angle θ0 = 30°. Table 1 also lists the retrieval results of aerosol PSDs using the SFLA and improved SFLA, respectively, and the corresponding inverse curves are depicted in Figure 4. It can be found that without random measurement noise, the value of relative standard deviation ξ obtained by the improved SFLA is smaller than that by the SFLA. When adding the random measurement noise to the relative ALS signals, the value of δ will increase, but it is still easy to find that the retrieval results obtained by the improved SFLA is acceptable. That is to say, the accuracy and robustness of the inverse results retrieved by the improved SFLA are higher than those by the SFLA, especially with random measurement noise.Table 1
Retrieval results of the aerosol PSDs by the SFLA and improved SFLA.
Aerosol PSDs
Noise,σR (%)
Improved SFLA
SFLA
D
¯
σ
ξ
D
¯
σ
ξ
L-ND¯,σ=3,2
0
3.00
2.00
0.00000
3.01
2.00
0.00474
3
2.95
1.98
0.01998
3.12
1.92
0.08235
5
3.08
1.96
0.04651
3.22
1.90
0.12335
Gammaβ,α=8,4
0
8.00
4.00
0.00000
8.01
4.00
0.00768
3
7.73
3.90
0.02543
8.14
4.09
0.08676
5
7.57
3.86
0.03971
8.26
4.16
0.14181
The mean inverse results and relative standard deviation of the 20 times calculations are shown in the tables.Figure 4
Reproducibility of the aerosol PSDs using the SFLA and improved SFLA, respectively.
(a)
(b)


### 3.2. Retrieval of Time-Varying Particle Number Concentration

Figure5 shows the real-time ALS signals. The black line denotes the signals without random noise, the gray line denotes the original observation signals, and the red line denotes the filtered results of the original observation signals using the KF algorithm. Figure 6 depicts the retrieval results with and without Kalman filter algorithm. The exact time-dependent particle number concentration is set as the sine wave type (19), and the standard deviation of measurement noise is assumed to be 0.05, i.e., σR = 0.05. It is obvious that if the measurement ALS signals with noise are used to reconstruct particle number concentration directly, the absolute value of maximal relative error reaches 11%. If the ALS signals are filtered by the Kalman filter algorithm beforehand, the absolute value of relative error will reduce quickly (less than 6%). The time-varying aerosol particle number concentration can be effectively reconstructed by employing the technique mentioned above.Figure 5
Different real-time ALS signals of the particle system.Figure 6
Reconstruction results of time-varying particle number concentration with and without the KF algorithm.Considering the applicability of the ALS method and Kalman filter algorithm, the following research should be done to study the effects of characteristic parameters in the measurement technique on the accuracy and stability of retrieval results to choose suitable parameter settings. Figure7 depicts estimated results obtained under different process noise. The sampling time interval Δt is set as 0.1 s. The standard deviation of measurement noise σR is set as 0.05. It can be found that reducing the standard deviation of process noise will result in larger deviation between real value and estimated result of aerosol particle number concentration and increase of time delay, i.e., requiring more time to response the mutation of particle number concentration. The amount of overshoot of peak value reduces with the decreasing standard deviation of process noise. These phenomena indicate that the dynamic performance and steady-state values of the system are reduced. The main reason is that decreasing standard deviation of process noise results in the reduction of filter gain matrix K (k), an important matrix that affects the correction weight, steady-state value, and dynamic performance of the system. Therefore, although the process noise and corresponding standard deviation are unknown beforehand in practical problem, a larger standard deviation of process noise should be selected to retrieve particle number concentration successfully.Figure 7
Effects of process noise on retrieval results: (a) linear-type distribution; (b) sine wave-type distribution.
(a)
(b)Figure8 shows the effects of measurement noise on the accuracy and stability of retrieval results. The sampling time interval Δt is set as 0.1 s. The standard deviation of process noise σQ is set as 0.01. From Figure 8, it can be found that increasing standard deviation of measurement noise will reduce the retrieval accuracy and increase time delay. The amount of overshoot of the peak value also shows increasing first and then reduces with the increasing of standard deviation of measurement noise. The reason is that according to the KF algorithm, the relationship between standard deviation of measurement noise and filtering gain K (k) is negative. Increasing the filtering gain K (k) will be conducive to reduce the time lag of retrieval results, and more valuable information obtained from measurement signals Y (k) rather than priori estimation vector X¯kk−1 and last estimated results N¯Dtk−1 will be used to improve the retrieval accuracy. Therefore, although the measurement noise is unknown in practice, we should better try our best to control it to improve the retrieval accuracy.Figure 8
Effects of measurement noise on retrieval results: (a) linear-type distribution; (b) sine wave-type distribution.
(a)
(b)During measurement experiments, except for these noises mentioned above affecting the measurement results, the sampling time intervalΔt also plays an important role in retrieving the time-varying particle number concentration. The standard deviations of process noise and measurement noise are set as 0.01 and 0.03, respectively. Figure 9 shows the effect of different sampling time intervals Δt on the retrieval results. It can be found that more accurate results and less time delay can be obtained, if the sampling time interval is reduced. The reason may be that more useful information about the particle system can be used to improve the retrieval accuracy when the sampling time interval Δt is smaller. Moreover, from Figure 9, it can also be found that although there are process noise and measurement noise existed in the measurement ALS signals, the retrieval results are still acceptable, if appropriate sampling time interval is selected. So, the ALS method combined with Kalman filter algorithm can be proposed as a kind of effective and reliable measurement technique to study the real-time estimation of aerosol concentration.Figure 9
Effects of sampling time intervalΔt on retrieval results: (a) linear-type distribution; (b) sine wave-type distribution.
(a)
(b)


## 3.1. Retrieval of Aerosol Particle Size Distributions

With the help of the original SFLA and improved SFLA, the aerosol PSDs are estimated from the relative ALS signals, and the real values are shown in Table1. According to our previous work [24], the measurement angles are set as θ = 5°, 10°, 15°, and 20°, and the reference angle θ0 = 30°. Table 1 also lists the retrieval results of aerosol PSDs using the SFLA and improved SFLA, respectively, and the corresponding inverse curves are depicted in Figure 4. It can be found that without random measurement noise, the value of relative standard deviation ξ obtained by the improved SFLA is smaller than that by the SFLA. When adding the random measurement noise to the relative ALS signals, the value of δ will increase, but it is still easy to find that the retrieval results obtained by the improved SFLA is acceptable. That is to say, the accuracy and robustness of the inverse results retrieved by the improved SFLA are higher than those by the SFLA, especially with random measurement noise.Table 1
Retrieval results of the aerosol PSDs by the SFLA and improved SFLA.
Aerosol PSDs
Noise,σR (%)
Improved SFLA
SFLA
D
¯
σ
ξ
D
¯
σ
ξ
L-ND¯,σ=3,2
0
3.00
2.00
0.00000
3.01
2.00
0.00474
3
2.95
1.98
0.01998
3.12
1.92
0.08235
5
3.08
1.96
0.04651
3.22
1.90
0.12335
Gammaβ,α=8,4
0
8.00
4.00
0.00000
8.01
4.00
0.00768
3
7.73
3.90
0.02543
8.14
4.09
0.08676
5
7.57
3.86
0.03971
8.26
4.16
0.14181
The mean inverse results and relative standard deviation of the 20 times calculations are shown in the tables.Figure 4
Reproducibility of the aerosol PSDs using the SFLA and improved SFLA, respectively.
(a)
(b)


## 3.2. Retrieval of Time-Varying Particle Number Concentration

Figure5 shows the real-time ALS signals. The black line denotes the signals without random noise, the gray line denotes the original observation signals, and the red line denotes the filtered results of the original observation signals using the KF algorithm. Figure 6 depicts the retrieval results with and without Kalman filter algorithm. The exact time-dependent particle number concentration is set as the sine wave type (19), and the standard deviation of measurement noise is assumed to be 0.05, i.e., σR = 0.05. It is obvious that if the measurement ALS signals with noise are used to reconstruct particle number concentration directly, the absolute value of maximal relative error reaches 11%. If the ALS signals are filtered by the Kalman filter algorithm beforehand, the absolute value of relative error will reduce quickly (less than 6%). The time-varying aerosol particle number concentration can be effectively reconstructed by employing the technique mentioned above.Figure 5
Different real-time ALS signals of the particle system.Figure 6
Reconstruction results of time-varying particle number concentration with and without the KF algorithm.Considering the applicability of the ALS method and Kalman filter algorithm, the following research should be done to study the effects of characteristic parameters in the measurement technique on the accuracy and stability of retrieval results to choose suitable parameter settings. Figure7 depicts estimated results obtained under different process noise. The sampling time interval Δt is set as 0.1 s. The standard deviation of measurement noise σR is set as 0.05. It can be found that reducing the standard deviation of process noise will result in larger deviation between real value and estimated result of aerosol particle number concentration and increase of time delay, i.e., requiring more time to response the mutation of particle number concentration. The amount of overshoot of peak value reduces with the decreasing standard deviation of process noise. These phenomena indicate that the dynamic performance and steady-state values of the system are reduced. The main reason is that decreasing standard deviation of process noise results in the reduction of filter gain matrix K (k), an important matrix that affects the correction weight, steady-state value, and dynamic performance of the system. Therefore, although the process noise and corresponding standard deviation are unknown beforehand in practical problem, a larger standard deviation of process noise should be selected to retrieve particle number concentration successfully.Figure 7
Effects of process noise on retrieval results: (a) linear-type distribution; (b) sine wave-type distribution.
(a)
(b)Figure8 shows the effects of measurement noise on the accuracy and stability of retrieval results. The sampling time interval Δt is set as 0.1 s. The standard deviation of process noise σQ is set as 0.01. From Figure 8, it can be found that increasing standard deviation of measurement noise will reduce the retrieval accuracy and increase time delay. The amount of overshoot of the peak value also shows increasing first and then reduces with the increasing of standard deviation of measurement noise. The reason is that according to the KF algorithm, the relationship between standard deviation of measurement noise and filtering gain K (k) is negative. Increasing the filtering gain K (k) will be conducive to reduce the time lag of retrieval results, and more valuable information obtained from measurement signals Y (k) rather than priori estimation vector X¯kk−1 and last estimated results N¯Dtk−1 will be used to improve the retrieval accuracy. Therefore, although the measurement noise is unknown in practice, we should better try our best to control it to improve the retrieval accuracy.Figure 8
Effects of measurement noise on retrieval results: (a) linear-type distribution; (b) sine wave-type distribution.
(a)
(b)During measurement experiments, except for these noises mentioned above affecting the measurement results, the sampling time intervalΔt also plays an important role in retrieving the time-varying particle number concentration. The standard deviations of process noise and measurement noise are set as 0.01 and 0.03, respectively. Figure 9 shows the effect of different sampling time intervals Δt on the retrieval results. It can be found that more accurate results and less time delay can be obtained, if the sampling time interval is reduced. The reason may be that more useful information about the particle system can be used to improve the retrieval accuracy when the sampling time interval Δt is smaller. Moreover, from Figure 9, it can also be found that although there are process noise and measurement noise existed in the measurement ALS signals, the retrieval results are still acceptable, if appropriate sampling time interval is selected. So, the ALS method combined with Kalman filter algorithm can be proposed as a kind of effective and reliable measurement technique to study the real-time estimation of aerosol concentration.Figure 9
Effects of sampling time intervalΔt on retrieval results: (a) linear-type distribution; (b) sine wave-type distribution.
(a)
(b)


## 4. Conclusions

A kind of multistep inversion method based on the angular light-scattering (ALS) measurement signals is proposed to solve the online monitoring of the aerosol PSDs and aerosol concentration. The aerosol PSD is retrieved by the shuffled frog-leaping algorithms (SFLAs) from relative ALS signals, and then the time-varying aerosol concentration is obtained by the Kalman filter algorithm from the real-time ALS signals. Moreover, the influences of characteristic parameters, i.e., sampling time interval, measurement noise, and process noise, on the estimated results have been studied. The conclusions drawn from the results are as follows:(1)
The convergence accuracy and robustness of the improved SFLA are higher than those of the original SFLA in studying the aerosol PSDs.(2)
With the aerosol PSDs known beforehand, the time-varying aerosol particle number concentration can be reconstructed effectively by the Kalman filter algorithm from the ALS signals even with process noise and measurement noise.(3)
The process noise and measurement noise have an important influence on the estimated results. Increasing standard deviation of process noise or reducing that of measurement noise benefit to obtaining more useful information to improve retrieval accuracy and reduce time delay. Moreover, the reliability and stability of retrieval results are also dependent on sampling time interval, and decreasing sampling time interval can improve the accuracy of retrieval results and reduce time delay to a certain degree. So, to improve retrieval accuracy and reduce time delay, the noise should be controlled, and suitable sampling time interval should be selected.The results obtained show that the proposed methodology is a promising approach for online monitoring of aerosol PSDs and time-varying aerosol particle number concentration. Further study will focus on performance improvement of the methodology as well as its applications in studying the nonspherical aerosol.



---

*Source: 1983460-2020-01-22.xml*