- Email: [email protected]
Application of an improved artificial bee colony algorithm to inverse problem of aerosol optical constants from spectral measurement data
Accepted Manuscript Title: Application of an improved artificial bee colony algorithm to inverse problem of aerosol optical constants from spectral measurement data Authors: Zhenzong He, Junkui Mao, Xingsi Han PII: DOI: Reference:
S0030-4026(17)30707-6 http://dx.doi.org/doi:10.1016/j.ijleo.2017.06.038 IJLEO 59301
To appear in: Received date: Accepted date:
15-3-2017 11-6-2017
Please cite this article as: Zhenzong He, Junkui Mao, Xingsi Han, Application of an improved artificial bee colony algorithm to inverse problem of aerosol optical constants from spectral measurement data, Optik - International Journal for Light and Electron Opticshttp://dx.doi.org/10.1016/j.ijleo.2017.06.038 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Application of an improved artificial bee colony algorithm to inverse problem of aerosol optical constants from spectral measurement data Zhenzong He Junkui Mao* Xingsi Han College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, P. R. China, 210016
*Corresponding author: Junkui Mao Jiangsu Province Key Laboratory of Aerospace Power System, College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, 29 Yudao Street, Qinhuai District, Nanjing, P.R. China, 210016 Tel: +86-1395173772 Email: [email protected] (J. K. Mao)
Abstract: An improved Artificial Bee Colony (IABC) algorithm is introduced to obtain the optical constants of aerosols from the spectral measurement data of the aerosol dispersion medium. The direct problem solved using the Finite Volume Method (FVM) combined with Mie scattering theory is used to study the radiative transfer process in the aerosol dispersion medium. Compared with standard Artificial Bee Colony (ABC) algorithm, the IABC can avoid local optima and improve convergence accuracy. Based on the IABC, the optical constants of aerosols over Harbin, China are retrieved under different random measurement errors. Results indicate that there is acceptable retrieval accuracy without errors, and the retrieval accuracy reduces if errors increase. To improve retrieval accuracy, two improved inverse models, namely as double-concentration and double-layer inverse models, are proposed. The investigation reveals that the improved inverse models, especially double-concentration inverse model, can give more accurate predictions even with errors. Finally, the optical constants of aerosols over Beijing available on the website of AERONET are also reconstructed by IABC satisfactorily. All the results show that the
IABC is a potential and effective optimal algorithm to retrieve the optical constants of aerosols, especially under improved inverse models.
Key words: Artificial Bee Colony algorithm; Aerosol; Optical constant; Inverse radiative problem
1. Introduction Aerosols play a crucial role in the atmospheric radiative balance by absorbing and scattering solar radiation (direct effects), and by acting as cloud condensation nuclei (indirect effects) [1, 2]. On a global scale, the main direct effect is the backscattering of solar radiation into space which results in a negative climate forcing. On a regional or even local scale, the influence of anthropogenic greenhouse gases is dominated by the absorption properties of aerosol particles, which can also heat the atmosphere, cause cloud dissipation and affect atmospheric circulation [2-4]. So, accurate measurements of aerosol absorption and scattering properties are important for estimating Earth’s energy balance effectively. To obtain the scattering and absorption properties of aerosol accurately, the detail information about the aerosol size, morphology and optical constant should be known first. Especially, the optical constant, a function of the chemical, phase composition and their mixing state, describes the optical properties and the interaction of the aerosol with radiation [5-7]. Therefore, obtaining effective optical constants of aerosols will allow us to better characterize the effects of aerosols on the climate changes and atmospheric optics. Generally, the optical constant m is composed of refractive index n and absorption index k (i.e. m=n±ik), and may vary over a wide range for different chemical components [8, 9]. Various techniques have been proposed to obtain the optical constants of particles, which can be generally divided into three categories [10]: the optoacoustic method, the direct sampling method and the inverse method with the remote sensing data. Ruan et al. [11] compared those three techniques and pointed out that the inverse method with the remote sensing data showed more satisfactory accuracy and broader range of application than the other two methods. To date, many scholars have focused their attentions on
studying aerosol optical constants by the narrowband remote sensing or wideband solar radiation [12-14]. Ku et al. [12] applied the spectral extinction measurements to obtain the optical constant, size distribution, and number density of soot particles in flames. Dombrovsky et al. [15] used the Kramers-Kronig (K-K) relation to calculate the optical constants of several types of diesel fuel and studied the absorption properties of the diesel fuel. Ruan et al. [11, 16] employed the Mie theory, the K-K relation and the inverse radiation problem solving techniques to retrieve the optical constants of fly-ash particles and aerosols from the spectral transmittance distribution of a cloud fly-ash or aerosol particles in potassium bromide (KBr) pellets. Combined with FTIP and Abbe refractometer measurements, Zhao et al. [17] proposed an inverse model based on the shooting method, Mie theory and the improved K-K relation to predict the optical constants of various infrared opcifiers effectively. The potential inverse problem solving techniques combined optical measurement methods above are available to obtain the optical constants of the aerosols, but the corresponding experiments are complicated and may bring in inevitable errors. Moreover, the application of K-K relation usually requires the integral of all spectrums, which will increase the difficulty of the experiments and the requirement of equipment. To overcome the difficulties, a spectral transmittance and reflectance measurement technique combined with an improved artificial bee colony (IABC) algorithm is developed to retrieve the optical constants of aerosols directly. To improve the retrieval results, two improved inverse models, i.e. the double-concentration inverse model and the double-layer inverse model, are proposed. The remainder of the paper is organized as follows. First, the IABC algorithm is introduced and compared with the standard ABC algorithm in predicting the optical constants of aerosols. Then, based on the IABC algorithm, the optical constants of aerosols over Harbin, China at different wavelengths are retrieved under different random measurement errors. Moreover, two improved inverse models are proposed to improve the retrieval accuracy of the results. Finally, the optical constants of aerosols over Beijing, China available on AERONET are retrieved by the improved inverse models. 2. Direct problem
Fig. 1 Figure 1 depicts a beam of collimated monochromatic laser impinges on a slab filled with homogeneous, absorbing, and scattering aerosols under common temperature (emitting effect of aerosols ignored at room temperature). The light transfer process can be described by the one-dimensional steady-state radiative transfer equation (RTE) as follows [18, 19]:
s
I ( z, s ) I ( z, s ) z 4
I ( z, s )Φ (s , s )d 4
i
i
i
(1)
where I ( z, si ) is the radiative intensity in direction s at location z; denotes the extinction coefficient;
s, is the scattering coefficient; (si , s ) is the scattering phase function. The mathematical definitions of , s, and (si , s ) are derived as follows:
s, ( si , s )
1
n
4
D N
4
n
2 j
0
f ( D j )Qext, j
(2)
D N
f ( D j )Qsca, j
(3)
j 1 n
j 1
4D N j 1
2 j
2 j
0
0
f ( D j )Qsca, j p, j ( D j , si , s )
(4)
where N 0 is the total number concentration of the sample; f ( D j ) is the size distribution function when the aerosol diameter is D j ; Qext , Qsca and p denote the extinction efficiency factor, scattering efficiency factor and scattering phase function of single aerosol, respectively, which can be predicted by using the Mie theory [20]. The boundary conditions of the sample can be expressed as:
0 I , I (0, ) 0 0, 0 / 2 I ( L, ) 0, / 2
(5) (6)
where I 0 is the total incident light intensity; I (0, ) and I ( L, ) denote the light intensity incident to the internal medium from the left side and right side of the sample, respectively; denotes the angle between the m and positive x axis; L (in m) is the geometrical thickness of the slab. Solving the RTE by numerical methods
such as the Finite Volume Method (FVM) and Monte-Carlo method (MCM), the radiative intensities on the two borders of the sample I (0, ) and I ( L, ) can be obtained and the spectral reflectance R and transmittance can be derived as follows:
I (0, ) cos sin d /2 I0
R 2 T 2
/2
0
I ( L, ) cos sin d . I0
(7) (8)
Roughly, for most dry aerosols, the refractive index of the optical constant is about 1.53, with a small absorption index of 0.005-0.01. If the aerosols absorb water vapor and become hygroscopic particles, the value of the optical constant will be diluted and decrease towards the value of pure water with optical constant of 1.33+0.00 i. For instance, the oceanic aerosols, mainly composed of hygroscopic sea salt, almost have no absorption, and the refractive index of the optical constant is around 1.38, which is just slightly larger than that of pure water. Soot or black carbon has a very high absorption with a representative optical constant of 1.75+0.45 i, although black carbon is most often found internally mixed with other diluent particles, especially for other organic materials [21, 22]. In this study, the optical constants of aerosols measured over Harbin, China are studied. The corresponding optical constants and aerosol size distribution are shown in Fig. 2.
Fig. 2
Considering the RTE is a non-linear integral-differential equation, simulated transmittance and reflection signals, which contain useful information for retrieving the optical constants of aerosols, under different thickness of the slab L or different total number concentration N0 will be different and linear independence. So, how to chose L and N0 is of great significance to obtain more useful information to retrieve the optical constants of aerosol more accurately. To make this point, a detailed examination of the sensitivity analysis of transmittance or reflection to the optical constants at different L or different N0 is employed to find out the optimal thickness of the slab and different particle concentration. The sensitivity coefficients are one of the most important characteristic parameters in the sensitivity analysis, which is the first derivative of the transmittance or reflectance to optical constants in present study [21]:
a ( )
(a a, ) (a a, ) , N0 or L a 2a
(9)
where a denotes the independent variable which stands for refractive index n and absorption index k; represents a tiny change and is set 0.5% in present work; is the measurement radiative signals, i.e. reflection R or transmission T . Figure 3 depicts the sensitivity coefficients of n and k with respect to reflection and transmission signals under different N0 at different wavelengths. It is obvious that if the sensitivity coefficients are either small or correlated with one another, the estimation problem will be very sensitive to measurement errors and difficult to be solved. Therefore, the optimal N0 region for inverse problems is the region with larger absolute value of sensitivity coefficient generally, i.e. N0 2 109 , 2 1011 / m3 . Similarly, the optimal region of L for inverse problems can be obtained from Fig.4, i.e. L 0.01,0.2m . Fig. 3 Fig. 4
3. Inversion algorithm 3.1 The principle of standard ABC algorithm Bees are social insects and live in colonies. Their behavior is governed by the goal of colony survival rather than being focused on the survival of individuals. When looking for food, the bees can be classified into three categories: employed bees, onlooker bees and scout bees. The bees that are going to discover food sources are named employed bees. The ones that are waiting in the hive for making decision to choose a food source are called onlooker bees. The other kind of bees that carries out random search for finding new sources is named scout bees. The employed bees bring back the nectar information from the food source to the hive and share it with onlooker bees by dancing on the dance area in the hive. The longer the duration of dance is, the larger the nectar amount of the food source exploited is. Onlooker bees watch the dances of many employed bees, and finally tend to choose a food source according to the probability proportional to the quality of the food source.
Therefore, the food sources with abundant nectar will attract more onlooker bees. Whenever a food source is exhausted, this food source will be abandoned and the corresponding onlooker bee will become scout bee or onlooker bee. Scout bees can be treated as performing the job of exploration, while employed and onlooker bees can be seen as performing the job of exploitation [23-25]. Being inspired by foraged behavior of the bee colony, Karaboga proposed the ABC algorithm and employed it to solve the problem of multi-variable and multi-modal continuous function optimization firstly [26]. Then, the ABC algorithm has drawn widespread attention and been popularized to tackle some practical optimization applications all over the world for its satisfactory global optimization and fewer control parameters [23, 25, 26]. In standard ABC algorithm, each food source represents a possible solution to the concerned problem and the abundant degree of the corresponding nectar implies the quality of the solution represented by the fitness value. So, the bee colony usually divides to two equal parts. One half is the employed bees whose number is the same as the number of food source, which make sure exactly one employed bee for every food source. The other half is the onlooker bees. At first, the food source positions Xi [ xi ,1 , xi ,2 ,
, xi ,Nn ](i = 1, 2,
, SN ) are initialized randomly
in the search space, which contains the food source sites [27], xi , j lowj rand1 (high j low j )
where i 1, 2,
, SN and j 1, 2,
(10)
, Nn , SN denotes the number of food source and Nn denotes the
dimension of unsolved problem; xi , j denotes the i th food source position at the j th dimension problem; high j and low j is the upper and lower limits of the search space
for the j th dimension problem; rand1
is the uniformly distributed random number in the range of [0, 1]. Then, evaluate Xi and record current personal best food source position Pi and pbesti . After initialization, the employed bees are all sent out for food. The i th employed bee produces a candidate food source position Vi [vi ,1 , vi ,2 ,
, vi ,Nn ] from the old personal best
one Pi in its memory depending on local information to find a neighboring food source as follows [27]:
Vi Pi i (Pi Pk )
(11)
where i is a uniformly distributed real random number in the range of [0, 1]; k 1, 2,
, SN and k i .
Then, evaluate the quality of the candidate food source position and calculate the corresponding fitness value Fiti [28]: 1 / (1 Fi ) if Fi 0 Fiti 1 abs( Fi ) if Fi 0
(12)
where Fi is the cost value of the food source position Vi . For maximization problem, the cost function can be directly used as a fitness function. Comparing the fitness function value of Xi and that of Vi , the one with more satisfactory fitness value will be selected and the corresponding position will be memorized. If Xi cannot be improved, its counter holding the number of trials is incremented by 1, otherwise, the counter is reset to 0. When the employed bees bring back the nectar information of the food sources, the onlooker bees choose a food source depending on the probability Pi [28]:
Pi
Fiti SN
Fit i 1
.
(13)
i
In this probabilistic selection scheme, as the nectar amount increases, the number of onlooker bees visiting corresponding food source increases. If a position of food source cannot be improved further after a predetermined number limit of cycles, the food source will be abandoned, and the corresponding employed bee will be changed to onlooker bee or scout bee. Assume that the abandoned source is Xn , n 1, 2,
, SN , the
scout bee will discover a new food source as follow [23]: xn, j low j rand2 (high j low j )
(14)
where rand2 is the uniformly distributed random number in the range of [0, 1]. At the end of every generation, personal best position Pi and corresponding fitness value pbesti of each food source are recorded. Meanwhile, the global best position Pg and corresponding fitness value gbest of all the food sources are also recorded.
3.2. The theory of the IABC algorithm
In standard ABC algorithm, the fast information flow between bees seems to be the reason for clustering of bees, which will also lead to the diversity decline rapidly and leave the ABC algorithm with great difficulties of escaping local optima. Moreover, when solving those unimodal problems, the convergence speed of the ABC algorithm is typically slower than those of representative population-based algorithm (e.g. genetic algorithm and particle swarm optimization) [26]. Therefore, to remedy this problem, an improved ABC algorithm (IABC) applies a logistic mapping to generate a chaotic sequence in its initialization and introduces the Differential Evolution (DE) algorithm to improve the personal best food source positions. In the ABC algorithm, the initial food source positions are usually generated by random numbers that may allow some of the initial positions too close to each other which will increase the convergence time, see Eq. (9). The Chaos theory is proposed to overcome those difficulties in this paper. Chaos is the highly unstable motion of deterministic systems in finite phase space that often exists in nonlinear systems. The typical logistic mapping to generate the Chaos signal is described as follows [29, 30]: r (0) rand3 r (k 1) r (k ) 1 r (k ) , k 0,1,
, Nn 1
(15)
where rand3 is a uniformly distributed random number in the range of [0, 1]; is control number, when 4.0 , the logistic mapping is in a fully chaotic state. So, Eq. (9) can be expressed as: xi , j lowj r ( j ) (high j low j )
(16)
Moreover, the DE algorithm is applied to improve the personal best position Pi of each food source and avoid local optima. The DE algorithm, developed by Storn and Price [31], is an evolutionary optimization algorithm based on population cooperation and competition of individuals and has been successfully applied to solve optimization problems particularly involving non-smooth fitness function [31-33]. The basic optimization process in the DE algorithm consists of mutation, crossover and selection [33]. Specifically, the mutation operation generates a new parameter vectors by adding the weighted difference between two population vectors to a third vector. Then, the mutated vector parameters are mixed with the parameters of another predetermined vector,
named as the target vector, to yield the so-called trial vector. This operation is named as crossover. The last operation, called selection, is to judge whether the fitness value yielded by the trial vector is superior to that of the target vector. If so, the trial vector replaces the target vector in the following generation. Each population vector has to serve once as the target vector so that SN competitions take place in one generation [31]. The mutant vector Pi,mut can be generated according to [31]: Pi,mut Pk F P j Pm
(17)
where k , l , m are random integers uniformly selected from 1, SN , ik l m ; F , the mutant factor, is a real and constant value, F0, 2 , which controls the amplification of the differential variation P j Pm . In order to increase the diversity of the perturbed parameter vectors, crossover is introduced and the trial vector Pi,tri ( pi1,tri , pi 2,tri , pi3,tri ..., piNn ,tri ) can be derived from [31]:
pij ,mut pij ,tri (t ) pij
if rand 4 ( j ) CR or j rnbr (i ) if rand 4 ( j ) CR or j rnbr (i )
(18)
where rand4 ( j ) is the j th evaluation of a uniform random number generator, rand4 ( j )0, 1 ; CR is the crossover constant CR0, 1 ; rnbr (i) is a random integer uniformly selected from 1, SN . The search procedures of the IABC algorithm are shown below and the corresponding flowchart is illustrated in Fig. 5. Step 1.
Input system control parameters of the IABC algorithm, i.e. the total number of the food sources SN , the maximum number of iterations M , the number of the inversion parameters (dimensions) N n , the searching space
= [lowi , highi ](i = 1,2,
, Nn ) , the tolerance for minimizing the fitness value , the
predetermined number of cycles limit and the values of , F , CR . Step 2.
Initialize the positions of food sources in the N n dimensional problem by mapping the chaotic sequence, which is generated by Eq. (16), to the search space. Evaluate the fitness value of each food source position Fiti . Obtain the initial personal best food source position Pi and the global best food source position Pg , and record corresponding fitness values pbesti and gbest . Set the number of current iteration t 0 .
Step 3.
Determine whether the program matches one of the following two stop criterions. If so, end the calculation and record the estimated values; if not, proceed to Step 4. (i) The fitness value of global best food source position is less than the tolerance , gbest ; (ii) The number of the iteration reaches the user-defined iteration limit M , t M .
Step 4.
Send employed bees to food sources, every employed bee produces a modification near the position of former personal best food source position Pi to create a candidate food source position Vi by using Eq.(11). Evaluate Vi and compare the fitness values of Pi and Vi , and record the position with more satisfactory fitness value as the personal best food source position. Then, judge whether the position of the i th food source has been improved. If not, its counter holding the number of trials Tri(i) is incremented by 1; if so, the counter is reset to 0.
Step 5.
Employed bees bring back the nectar information of food sources, and share it with onlooker bees. Onlooker bees choose food sources with probabilities related to nectar amount of food sources according to Eq. (13).
Step 6.
Judge whether Tri(i) > limit , i = 1, 2,
, Nn . If so, abandon the i th food source position and
change the corresponding employed bee to scout one. The scout bee discovers new food source by using Eq. (14). Step 7.
Calculate the mutant vectors for current personal best positions of all food sources using Eq. (17). Operate the crossover and calculate the trial vector Pi ,tri using Eq. (18). Finally, evaluate the trial vector and compare the fitness value of the trial vector and that of the target vector (current personal best positions of food source). If the trial vector is superior, update the fitness value pbesti and corresponding personal best food source position Pi .
Step 8.
Determine whether the fitness values of the personal best food source positions pbesti are superior to those of the current global best food source position gbest . If so, update the value gbest and
corresponding location Pg . Loop to Step 3. Fig. 5 4. Inverse models and simulation The retrieval of the aerosol optical constants is solved by minimizing the fitness function value Fit , which is defined as the sum of the square residual between the estimated and measured signal ratios. The mathematical expression of Fit is derived as follows: 2 2 1 Rest Rmea Test Tmea Fit . 2 Rmea Tmea
(19)
where Rest and Test are the estimated transmittance and reflectance; Rmea and Tmea denote the real transmittance and reflectance. Considering that the algorithm is a stochastic optimization method and has certain randomness, all the calculations are repeated 50 times. The standard deviation of the inverse results is investigated to evaluate the reliability and feasibility of the inverse results, and the mathematical expressions are described as follows:
n
nest
1 50
50
nest nest, i , 2
1 50
k
i 1
1 50
50
i 1
nest, i ,
kest
1 50
50
k
est
kest, i
2
(20)
i 1
50
k
est, i
(i 1,
,50)
(21)
i 1
4.1 Comparison of standard ABC and IABC algorithms At first, the performance of the IABC is investigated through comparison with standard ABC. Figure 6 displays the comparison of the fitness function values between standard ABC and the IABC in retrieving the optical constants of aerosols over Harbin, China at wavelength =8.73 m (m=1.39655+0.14207 i). The geometric thickness of the measurement sample L is set as 0.01 m, and the total number concentration N0 is set as 2.0× 1010/m3. The aerosol size is assumed to vary from 0.001~10 m, which covers the main size range of aerosols [34, 35]. The corresponding system control parameters are listed in Table 1. Table 2 lists the corresponding retrieval results. The investigation shows that the fitness function value of the IABC converges much faster than the
standard ABC. Moreover, the IABC can arrive at the lowest best fitness value than the standard ABC within a smaller number of generations (see Fig. 6). All the results demonstrate that the IABC can avoid the phenomenon of local optima and low convergence accuracy which exit in the standard ABC.
Table 1 Table 2 Fig. 6
4.2 Retrieval of the optical constants of aerosols The optical constants of aerosols over Harbin, China at different wavelengths (between 6.0 m and 17.0 m) are estimated. The geometrical thickness of the slab L is set as 0.01 m, and the total number concentration N0 is set as 2.0×1010/m3. The system control parameters of the IABC are shown in Table 1. The searching ranges of the refractive and absorption indices are set as [1.0, 3.0] and [0.001, 1.0], respectively. Figures 7-9 depict the retrieval results of the optical constants with different random measurement errors added to the measured reflectance and transmittance. Figure 7 indicates that the inverse results of the optical constants show high accuracy without random measurement errors. The largest absolute errors of the refractive and absorption indices are at =6.68 m, which are less than 0.04 and 0.09, respectively. If the random measurement error is added, the retrieval accuracy will reduce and the standard deviation and absolute errors become larger (Figs. 8 and 9). When the random error is added up to 5%, the maximum absolute errors of the refractive and absorption indices increase to 0.87 and 0.17, respectively. Moreover, it is obvious that there is similar retrieval accuracy between the retrieval results of the refractive index and those of absorption index.
Fig. 7
Fig. 8
Fig. 9
4.3 The improved inverse models The light transmittance and reflectance measurement technique provides a potential method to retrieve the optical constants of aerosols. Without random errors, the estimated results are satisfactory. However, when there is a random error, the inverse results deteriorate (see Figs. 8 and 9). To remedy this situation, two improved inverse models, i.e. the double-concentration and the double-layer inverse models, are proposed to provide more useful radiative measurement information related to the optical constants to improve the retrieval accuracy. The original inverse model (single-concentration and single-layer inverse model) studied in Part 4.2 is named as Model 1.
Model 2: The double-concentration inverse model In this model, the geometrical thickness of the measurement sample remains unchanged, and the concentration of the aerosol dispersion medium is changed ( N0 =2.0 1010 / m3 and N0 =2.0 1011 / m3 ). The searching ranges of the real and imaginary parts are also set as [1.0, 3.0] and [0.001, 1.0], respectively. The fitness function Fit can be defined as follow: 2 2 2 2 Test,1 Tmea,1 Rest,2 Rmea,2 Test,2 Tmea,2 1 Rest， 1 Rmea,1 Fit . 4 Rmea,1 T R T mea,1 mea,2 mea,2
(22)
Model 3: The double-layer inverse model In this model, the concentration of the aerosol dispersion medium remains unchanged, and the geometrical thickness of the slab is changed ( L1 =0.1 m and L2 =0.2 m ). The searching ranges of the refractive index and are also set as [1.0, 3.0] and [0.001, 1.0], respectively. Thus, the fitness function Fit can also be defined as Eq. (22). Compared with Model 2, Model 3 needs two optical vessels with different geometrical thickness during actual experiement, which will increase the cost of experiment.
4.4 Comparison of the retrieval results under different inverse models Table 3 lists the estimated results of optical constants of aerosols at wavelengths =6.97 m and 11.07 m under the Model 1, Model 2 and Model 3, respectively. All calculations are repeated 50 times. It can be found that the average retrieval results of multiple calculations under Models 2 and 3 are more accurate than those under Model 1, and the standard errors under Models 2 and 3 are less than those under Model 1. The reason is that two improved inverse models provide more useful information about the optical properties of aerosols than Model 1, although it may be more time-consuming. Moreover, comparing the estimated results of two improved inverse models, it is easy to find that the accuracy of estimated results under Model 2 shows little better than that under Model 3, when there are the same random measurement errors.
Table 3
Figure 10 depicts the distribution of fitness function values at =6.97 m and 11.07 m under three different inverse models. From Fig. 10(a), it can be found that the minimum value region of fitness function is not a point but looks like a canyon, which indicates that the optimal solution is not unique and more than one appropriate value for absorption index are available to satisfy minimum fitness function under the same refractive index value. In other words, the multiple-value characteristics of the retrieval results exist and there are many couples of (n, k) values corresponding to the same experimental transmittance and reflectance. This is the difficulty of inverse estimated of aerosol optical constants, which is due to the complication of the scattering amplitude functions in the precise Mie theory. Moreover, comparing the distribution of the fitness function values under different inverse models, it is obvious that the minimum fitness function value regions under the improved inverse models are smaller than those under Model 1. This phenomenon indicates that corresponding to the same experimental transmittance and reflectance, less couples of (n, k) values may exist under the improved inverse models than
those under Model 1, which can explain why the retrieval results of the improved inverse models are better than those of Model 1. Fig. 10 4.5 Study of aerosols available on the AERONET Figure 11 illustrates the aerosol size distribution and the optical constants over Beijing, China available on the AERONET. The retrieval results under two improved inverse models are shown in Fig. 12. It can be found that the average estimated results of multiple calculations under the improved inverse models, especially Model 2, are more accurate than Model 1 in solving the problem of aerosol optical constants, even with 3% random measurement errors. Moreover, comparing the inverse results of refractive index and those of absorption index, it is obvious that the estimated accuracy of absorption index shows little better than that of refractive index. The phenomenon is attributed to the multi-value characteristics of the estimated results of refractive index, which can also be demonstrated by studying the distribution of fitness function values. For simplifying the manuscript, the study of fitness function values for the aerosols over Beijing is not available here, and interested readers can contact us for the details.
Fig. 11
Fig. 12 5. Conclusions Two improved inverse models, named double-concentration and double-layer inverse models, combined with an improved artificial bee colony (IABC) algorithm are proposed to improve the retrieval accuracy of aerosol optical constants. The results demonstrate that the improved inverse models combined with the IABC are suitable for retrieving aerosol optical constants. The following conclusions can be drawn: (1) Compared with standard ABC, the IABC can avoid the local optima and obtain higher convergence accuracy within a smaller number of iterations in solving the problem of retrieving aerosol optical constants.
(2) The retrieval results of two improved inverse models (i.e. the double-concentration inverse model and double-layer inverse model) are proved to be more accurate than the original inverse model, because more useful optical information about the aerosol particles can be obtained, and the effect of multi-value characteristics of the retrieval results can be weakened. Moreover, comparing the results of these two improved inverse models, the double-concentration inverse model shows more retrieval accurate than the double-layer inverse model. (3) Combined with the IABC, the optical constants of aerosols over Beijing, China obtained from AERONET can be estimated effectively under the improved inverse models, even with random measurement errors.
Acknowledgements The supports of this work by the Foundation of Nanjing University of Aeronautics and Astronautics of China for new teachers start work (No. 1022-YAH16051), the National Natural Science Foundation of China (No: 51606095) and Jiangsu Provincial Natural Science Foundation (No: BK20160794) are gratefully acknowledged. The author would also like to thank all the investigators of AERONET for their free provision of the aerosol properties over Beijing. A very special acknowledgement is made to the editors and referees who make important comments to improve this paper. References 1.
Dinar E, Riziq AA, Spindler C, Erlick C, Kiss G, Rudich Y. The complex refractive index of atmospheric and model humic-like substances (HULIS) retrieved by a cavity ring down aerosol spectrometer (CRD-AS). Faraday Discussions 2008; 137: 279-295.
2.
Ebert M, Weinbruch S, Hoffmann P, Ortner HM. The chemical composition and complex refractive index of rural and urban influenced aerosols determined by individual particle analysis. Atmospheric Environment 2004; 38: 6531–6545.
3.
Ebert M, Weinbruch S, Rausch A, Gorzawski G, Helas G, Hoffmann P, Wex H. Complex refractive index of aerosols during LACE 98#x2010; as derived from the analysis of individual particles. Journal of Geophysical Research Atmospheres 1999; 3(2): 275-280.
4.
Noia AD, Hasekamp OP, Van Harten G, Rietjens JHH, Smit JM, Snik F, Henzing JS, DeBoer J, Keller CU, Volten H. Use of neural networks in ground-based aerosol retrievals from multi-angle spectropolarimetric observations. Atmospheric Measurement Techniques 2015; 7: 281-299.
5.
Spindler C, Riziq AA, Rudich Y. Retrieval of Aerosol Complex Refractive Index by Combining Cavity Ring Down Aerosol Spectrometer Measurements with Full Size Distribution Information. Aerosol Science & Technology 2007; 41(11): 1011-1017.
6.
Flores JM, Washenfelder RA, Adler G, Lee HJ, Segev L, Laskin J, Laskin A, Nizkorodov SA, Brown SS, Rudich Y. Complex refractive indices in the near-ultraviolet spectral region of biogenic secondary organic aerosol aged with ammonia. Physical Chemistry Chemical Physics 2014; 16(22): 10629-10642.
7.
Guo LF, Shen JQ, Wan C. Particle analysis based on light scattering of particles illuminated by a divergent Gaussian beam. Optics and Lasers in Engineering 2013; 51: 826-831.
8.
Miles RE, Walker JS, Burnham DR, Reid JP. Retrieval of the complex refractive index of aerosol droplets from optical tweezers measurements. Physical Chemistry Chemical Physics Pccp 2012; 14(9): 3037-3047.
9.
Zarzana KJ, Cappa CD, Tolbert MA. Sensitivity of Aerosol Refractive Index Retrievals Using Optical Spectroscopy. Aerosol Science & Technology 2014; 48(11): 1133-1144.
10. Li JY, Dong SK, Tan HP. Effective Optical Constant of Alumina Particle Containing Carbon. Journal of Thermophysics & Heat Transfer 2009; 23(1): 216-219. 11. Ruan LM, Wang XY, Qi H, Wang SG. Experimental investigation on optical constants of aerosol particles. Journal of aerosol science 2011; 42(11): 759-770. 12. Ku JC, Felske JD. Determination of refractive indices of Mie scatterers from Kramers-Kronig analysis of spectral extinction data. Journal of the Optical Society of America A 1986; 3(5): 617-623. 13. Herman BM, Browning RS, De Luisi JJ. Determination of the Effective Imaginary Term of the Complex Refractive Index of Atmospheric Dust by Remote Sensing: The Diffuse-Direct Radiation Method. Journal of
the Atmospheric Sciences 1975; 32(5): 918-925. 14. Wei DJ, Qiu JH. Wide-Band Method to Retrieve the Imaginary Part of Complex Refractive Index of Atmospheric Aerosols. Part I: Theory, Scientia Atmospherica Sinica, 1998. 15. Dombrovsky LA, Sazhin SS, Mikhalovsky SV, Wood R, Heikal MR. Spectral properties of diesel fuel droplets Fuel 2003; 82(1): 15-22. 16. Ruan LM, Qi H, An W, Tan HP. Inverse radiation problem for determination of optical constants of fly-ash particles. International Journal of Thermophysics 2007; 28(4):1322-1341. 17. Zhao JJ, Duan YY, Wang XD, Zhang XR, Han YH, Gao YB, Lv ZH, Yu HT, Wang BX. Optical and radiative properties of infrared opacifier particles loaded in silica aerogels for high temperature thermal insulation. International Journal of Thermal Sciences 2013; 70(8): 54-64. 18. Modest MF. Radiative heat transfer, McGraw-Hill New York 2013. 19. Ren YT, Qi H, Chen Q, Ruan LM, Tan HP. Simultaneous retrieval of the complex refractive index and particle size distribution. Optics Express 2015; 23(15): 19328-19337. 20. Bohren CF, Huffman DR. Absorption and scattering of light by small particles. Wiley. Com 2008. 21. Chudzyński S, Czyżewski A, Ernst K, Karasiński G, Kolacz K, Pietruczuk A, Skubiszak W, Stacewicz T, Stelmaszczyk K, Szymański A. Multiwavelength lidar for measurements of atmospheric aerosol. Optics and Lasers in Engineering 2002; 37: 91-99. 22. Lenoble J, Remer L, Tanre D. Aerosol Remote Sensing. Springer 2013. 23. Gao W, Liu S, Huang L. A global best artificial bee colony algorithm for global optimization. Journal of Computational & Applied Mathematics 2012; 236(11): 2741–2753. 24. Akay B., Karaboga D. Artificial bee colony algorithm for large-scale problems and engineering design optimization. Journal of Intelligent Manufacturing 2012; 21(21): 1-14. 25. Singh A. An artificial bee colony algorithm for the leaf-constrained minimum spanning tree problem. Applied
Soft Computing 2009; 9(2): 625-631. 26. Karaboga D, Akay B. A comparative study of Artificial Bee Colony algorithm. Applied Mathematics & Computation 2009; 214(1): 108-132. 27. Akay B, Karaboga D. A modified Artificial Bee Colony algorithm for real-parameter optimization. Information Sciences 2012; 192(6): 120-142. 28. Maghali I, Oliveira SD, Schirru R, Medeiros JACCD. On the Performance of an Artificial Bee Colony Optimization Algorithm Applied to the Accident Diagnosis in a PWR Nuclear Power Plant. 2009 International Nuclear Atlantic Conference (INAC2009), Sep. 27th-Oct. 2nd, 2009, Rio de Janeiro, Brazil. 29. Liu F, Duan H, Deng Y. A chaotic quantum-behaved particle swarm optimization based on lateral inhibition for image matching. Optik-International Journal for Light and Electron Optics 2012; 123(21): 1955-1960. 30. Zhang B, Qi H, Sun SC, Ruan LM, Tan HP. Solving inverse problems of radiative heat transfer and phase change in semitransparent medium by using Improved Quantum Particle Swarm Optimization. International Journal of Heat and Mass Transfer 2015; 85: 300-310. 31. Storn R, Price K. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization 1997; 11(4): 341-359. 32. Wu Y, Lu J, Sun Y. An improved differential evolution for optimization of chemical process. Chinese Journal of Chemical Engineering 2008; 16(2): 228-234. 33. Lu S, Sun C, Lu Z. An improved quantum-behaved particle swarm optimization method for short-term combined economic emission hydrothermal scheduling. Energy Conversion and Management 2010; 51(3): 561-571. 34. Yuan Y, Yi HL, Shuai Y, Liu B, Tan HP. Inverse problem for aerosol particle size distribution using SPSO associated with multi-lognormal distribution model. Atmospheric Environment 2011; 45(28): 4892-4897. 35. Tang H, Lin JZ. Retrieval of spheroid particle size distribution from spectral extinction data in the
independent mode using PCA approach. Journal of Quantitative Spectroscopy and Radiative Transfer 2013; 115: 78-92. 36. Aerosol robotic network (AERONET) 2014, NASA,
Figure Captions Fig. 1
The schematic model of the parallel slab with collimated laser incident
Fig. 2
The size distribution and optical constants of aerosols over Harbin, China
Fig. 3
Sensitivity coefficients of n and k with respect to reflection or transmission signals under different N0 at different wavelengths
Fig. 4
Sensitivity coefficients of n and k with respect to reflection or transmission signals under different L at different wavelengths
Fig. 5
Flowchart of the IABC algorithm
Fig. 6
Comparison of fitness function values of standard ABC and IABC algorithms
Fig. 7
Reproducibility of optical constants of aerosols without random error.
Fig. 8
Reproducibility of optical constants of aerosols with 3% random error.
Fig. 9
Reproducibility of optical constants of aerosols with 5% random error.
Fig. 10
The distribution of the fitness function values under different inverse models with different optical constants when 1.0
Fig. 11
The size distribution and optical constants of aerosols over Beijing available on AERONET [36]
Fig. 12
The retrieval results of optical constants of aerosols over Beijing under different ramdon measurement errors
Aerosol sample Incident Laser
Reflection signals Fig. 1
Transmittance signals
10 Refractive index Absorption index
3.0 2.5 2.0 1.5 1.0
-1
1.5
10
1.4
10
-2
0.5 1.3
0.0 1
10
Diameter (m)
(a) Aerosol size distribution Fig. 2
100
-3
6
7
8
9
10 11 12 13 14 Wavelength, (m)
(b) Optical constants
15
16
10 17
Absorption index, k
3.5
0
1.6
Harbin, China, on 6th Mar. 2014 True distribution
Refractive index, n
Volume frequency distribution, f(D)
4.0
0.035
0.010
(b)
(a)
0.030 0.025 0.020
n T )
nR)
0.005
L=0.01m -3 =2.58m, n=1.40, k=1.14x10 =8.73m, n=1.39, k=0.142 =10.48m, n=1.47, k=0.137 =12.53m, n=1.46, k=0.139
0.015 0.010 0.005
0.000
0.000 -0.005 -0.005 8 10
10
9
10
10 3 N0(1/m )
10
11
10
12
0.005
(d)
(c) 0.000
10
10 3 N0(1/m )
11
10
11
10
10
12
0.05 0.04
0.01 k T )
kR)
9
0.02
-0.010 -0.015
Fig. 3
10
0.03
-0.005
-2 -4 -6 -8 -10 -12 -14 -16 -18 8 10
-0.010 8 10
0.00
10
9
10
10 3 N0(1/m )
10
11
10
12
-2 -4 -6 -8 -10 -12 8 10
10
9
10
10 3 N0(1/m )
10
12
0.010 (a)
0.020 (b) 0.015 0.010
10
3
=2.58m, n=1.40, k=1.14x10 =8.73m, n=1.39, k=0.142 =10.48m, n=1.47, k=0.137 =12.53m, n=1.46, k=0.139
-3
n T )
nR)
0.005
N0=1x10 /m
0.005 0.000
0.000
-0.005 -0.005 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m) 0.005 (c)
-0.010 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m) 0.035
0.030
0.000
0.025 0.020
-0.010
k T )
kR)
-0.005
0.015 0.010 0.005
-0.015 -10 -12 -14 -16 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m)
Fig. 4
0.000 -4 -8 -12 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m)
Start Input system datum and control parameters
Initialization t=0
Y
gbest ?
N t=t+1 i=0 i=i+1 Evaluate candidate food position Vi and calculate fitness value Fiti
N
Fiti pbesti ? Y Update the personal best position Pi Vi
Onlooker bees choose food sources N
Tri(i ) limit ? Y Abandon the ith position, discover new position
Evaluate Pi ,tri and calculate the fitness value Fiti ,tri N
Fiti ,tri pbesti ?
Y Update the personal best position Pi Pi ,tri i SN ? N
N
Y
pbesti gbest
Y Update the global best position Pg Pi tM ?
N Output result End
Fig. 5
Y
-1
10
ABC IABC
-3
10
-5
10
-7
Fit
10
-9
10
-11
10
-13
10
-15
10
-17
10
0
200
400 600 Interation
800
Fig. 6
28
1000
-3
10
-5
10
1.3
-7
10 1.2
-9
10 9
10 11 12 13 Wavelength, (m)
14
15
16
10
0
10 -1
-2
10
10
-4
10 -2
10
-6
10
-11
8
2
10
-3
10 17
10
-8
6
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
10 17
Fig. 7
10
1.4
10
1.2
10
1.0
10
0
-1
-2
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
Absolute error
0
10
3
-1
10
-2
10
-3
10
10
1
10
-1
10
-3
10
-3
6
10 Real value Inverse results, Err=3%
-4
10 17
10
Absolute error of inverse results
1
1.6
5
10
10
Absorption Index, k
Absolute error
Absolute error of inverse results
Refractive Index, n
Real values Inverse results, Err=3%
0.8
1
2
1.8
-5
6
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
10 17
Fig. 8
2
10
Refractive Index, n
1.8
1
10
1.5 0
10
1.2 -1
10
0.9
0.6
1
10
-2
10
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
10 Real value Inverse results, Err=5%
Absolute error
0
10
2
-1
10
-2
10
-3
10
10
0
10
-2
10
-4
10
-3
6
4
10
Absorption Index, k
Aboslute error
Aboslute error of inverse results
Real values Inverse results, Err=5%
-4
10 17
10
Fig. 9
29
-6
6
7
8
9
10 11 12 13 14 Wavelength, (m)
15
16
10 17
Absolute error of inverse results
3
2.1
Absolute error of inverse results
1.4
Absolute error
0
Absorption index, k
Refractive index, n
-1
10
7
10 Real value Inverse results, Err=0%
1
10
1.5
6
4
10
10
Absolute error
Absolute error of inverse results
Real values Inverse results, Err=0%
1.1
1
3
1.6
(a)
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
n
n
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
1.5
1 0
0.2
0.4
k
0.6
0.8
1 0
1
3
0.2
0.4
k
0.6
0.8
1
3
(c)
(d) Fitness
Fitness 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2.5
2
n
n
2.5
1.5
2
1.5
1 0
0.2
0.4
k
0.6
0.8
1 0 3
1
3
0.2
0.4
k
0.6
0.8
1
(f)
(e) Fitness 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
1.5
1.5
1 0
Fitness
2.5
2.5
n
1
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
Fitness
2.5
1.5
Fitness
(b)
Fitness
2.5
n
8
3
3
0.2
0.4
0.6 k =6.97m, n=1.39, k=0.08284
0.8
1
Fig. 10
30
1 0
0.2
0.4
0.6 k =11.07m, n=1.47, k=0.131
0.8
1
Volume size distribution, N0f(D)
x10-2 0.9
Beijing, China at 22:51:24 on 18 May 2014 k n
Beijing, China at 22:51:24 on 18 May 2014
Refractive index, n
0.20
0.15
0.10
1.45
0.6
1.40
0.3
0.05
1.35 0.4
0.00 0.1
1 Radius (m)
(a) aerosol size distribution Fig. 11
10
0.5
0.6
(b) optical constants
31
0.7 0.8 0.9 Wavelength, (m)
1.0
0.0 1.1
Absorption Index, k
1.50
0.25
1.52 Real values Inverse results, Model 2, Err=0% Inverse results, Model 3, Err=0%
x10-3 Real values Inverse results, Model 2, Err=0% Inverse results, Model 3, Err=0%
8.0 7.5
Absorption Index, k
Refractive index, n
1.48
8.5
1.44
1.40
7.0 6.5 6.0 5.5 5.0 4.5
1.36 0.4
0.5
0.6
0.7 0.8 0.9 Wavelength, (m)
1.0
4.0 0.4
1.1
16
2.4 Real values Inverse results, Model 2, Err=5% Inverse results, Model 3, Err=5%
0.6
0.7 0.8 0.9 Wavelength, (m)
1.0
1.1
x10-3 Real values Inverse results, Model 2, Err=5% Inverse results, Model 3, Err=5%
14 12
2.0
Absorption Index, k
Refractive index, n
2.2
0.5
10
1.8 1.6 1.4 1.2 1.0 0.4
8 6 4 2 0
0.5
0.6
0.7 0.8 0.9 Wavelength, (m)
1.0
1.1
Fig. 12
32
0.4
0.5
0.6 0.7 0.8 Wavelength, (m)
0.9
1.0
1.1
Table 1
The system control parameters of the standard ABC and IABC algorithms
Parameters ABC IABC
SN 50 50
Nn 2 2
M
0
limit
F
CR
1000
10-16
100
—
—
—
1000
10-16
100
4.0
0.5
0.4
33
Table 2
The retrieval results of the optical constants by standard ABC and IABC algorithms
Errors Algorithms
IABC
Parameters
0%
3%
5%
n
1.39655
1.36343
1.35979
n
4.2×10-8
0.01389
0.01043
k
0.14207
0.13602
0.11252
k
2.6×10-7
0.04484
0.05649
n
1.39654
1.35970
1.34258
n
4.9×10-6
0.01877
0.03852
k
0.14196
0.12181
0.10754
k
1.7×10-4
0.08285
0.08137
ABC
34
Table 3
The retrieval results of the optical constants by different inverse models
Models
Model 1
=6.97m, n=1.39, k=0.08284
=11.07m, n=1.47, k=0.131
Err=0%
Err=3%
Err=5%
Err=0%
Err=3%
Err=5%
n
1.40866
1.16026
1.56586
1.47345
1.34190
1.30214
n
4.8×10-7
0.19485
0.24662
1.8×10-7
0.16306
0.29876
k
0.11098
0.24095
0.35821
0.13100
0.19843
0.28753
k
1.4×10-6
0.11173
0.24378
2.8×10-7
0.12599
0.34621
n
1.39377
1.32328
1.30076
1.47344
1.41366
1.35234
n
2.4×10-7
0.15774
0.19371
3.5×10-8
0.06476
0.18762
k
0.08285
0.13264
0.14454
0.13100
0.16147
0.20629
k
2.0×10-7
0.06512
0.08754
2.5×10-9
0.11835
0.32768
n
1.39377
1.28401
1.23445
1.47344
1.36311
1.31030
n
2.3×10-7
0.15340
0.20131
1.5×10-7
0.21500
0.11703
k
0.08285
0.16095
0.21629
0.13100
0.18349
0.26876
k
1.8×10-7
0.26462
2.2×10-7
0.12411
0.15248
Parameters
Model 2
Model 3 0.11608
35
S0030-4026(17)30707-6 http://dx.doi.org/doi:10.1016/j.ijleo.2017.06.038 IJLEO 59301
To appear in: Received date: Accepted date:
15-3-2017 11-6-2017
Please cite this article as: Zhenzong He, Junkui Mao, Xingsi Han, Application of an improved artificial bee colony algorithm to inverse problem of aerosol optical constants from spectral measurement data, Optik - International Journal for Light and Electron Opticshttp://dx.doi.org/10.1016/j.ijleo.2017.06.038 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
Application of an improved artificial bee colony algorithm to inverse problem of aerosol optical constants from spectral measurement data Zhenzong He Junkui Mao* Xingsi Han College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing, P. R. China, 210016
*Corresponding author: Junkui Mao Jiangsu Province Key Laboratory of Aerospace Power System, College of Energy and Power Engineering, Nanjing University of Aeronautics and Astronautics, 29 Yudao Street, Qinhuai District, Nanjing, P.R. China, 210016 Tel: +86-1395173772 Email: [email protected] (J. K. Mao)
Abstract: An improved Artificial Bee Colony (IABC) algorithm is introduced to obtain the optical constants of aerosols from the spectral measurement data of the aerosol dispersion medium. The direct problem solved using the Finite Volume Method (FVM) combined with Mie scattering theory is used to study the radiative transfer process in the aerosol dispersion medium. Compared with standard Artificial Bee Colony (ABC) algorithm, the IABC can avoid local optima and improve convergence accuracy. Based on the IABC, the optical constants of aerosols over Harbin, China are retrieved under different random measurement errors. Results indicate that there is acceptable retrieval accuracy without errors, and the retrieval accuracy reduces if errors increase. To improve retrieval accuracy, two improved inverse models, namely as double-concentration and double-layer inverse models, are proposed. The investigation reveals that the improved inverse models, especially double-concentration inverse model, can give more accurate predictions even with errors. Finally, the optical constants of aerosols over Beijing available on the website of AERONET are also reconstructed by IABC satisfactorily. All the results show that the
IABC is a potential and effective optimal algorithm to retrieve the optical constants of aerosols, especially under improved inverse models.
Key words: Artificial Bee Colony algorithm; Aerosol; Optical constant; Inverse radiative problem
1. Introduction Aerosols play a crucial role in the atmospheric radiative balance by absorbing and scattering solar radiation (direct effects), and by acting as cloud condensation nuclei (indirect effects) [1, 2]. On a global scale, the main direct effect is the backscattering of solar radiation into space which results in a negative climate forcing. On a regional or even local scale, the influence of anthropogenic greenhouse gases is dominated by the absorption properties of aerosol particles, which can also heat the atmosphere, cause cloud dissipation and affect atmospheric circulation [2-4]. So, accurate measurements of aerosol absorption and scattering properties are important for estimating Earth’s energy balance effectively. To obtain the scattering and absorption properties of aerosol accurately, the detail information about the aerosol size, morphology and optical constant should be known first. Especially, the optical constant, a function of the chemical, phase composition and their mixing state, describes the optical properties and the interaction of the aerosol with radiation [5-7]. Therefore, obtaining effective optical constants of aerosols will allow us to better characterize the effects of aerosols on the climate changes and atmospheric optics. Generally, the optical constant m is composed of refractive index n and absorption index k (i.e. m=n±ik), and may vary over a wide range for different chemical components [8, 9]. Various techniques have been proposed to obtain the optical constants of particles, which can be generally divided into three categories [10]: the optoacoustic method, the direct sampling method and the inverse method with the remote sensing data. Ruan et al. [11] compared those three techniques and pointed out that the inverse method with the remote sensing data showed more satisfactory accuracy and broader range of application than the other two methods. To date, many scholars have focused their attentions on
studying aerosol optical constants by the narrowband remote sensing or wideband solar radiation [12-14]. Ku et al. [12] applied the spectral extinction measurements to obtain the optical constant, size distribution, and number density of soot particles in flames. Dombrovsky et al. [15] used the Kramers-Kronig (K-K) relation to calculate the optical constants of several types of diesel fuel and studied the absorption properties of the diesel fuel. Ruan et al. [11, 16] employed the Mie theory, the K-K relation and the inverse radiation problem solving techniques to retrieve the optical constants of fly-ash particles and aerosols from the spectral transmittance distribution of a cloud fly-ash or aerosol particles in potassium bromide (KBr) pellets. Combined with FTIP and Abbe refractometer measurements, Zhao et al. [17] proposed an inverse model based on the shooting method, Mie theory and the improved K-K relation to predict the optical constants of various infrared opcifiers effectively. The potential inverse problem solving techniques combined optical measurement methods above are available to obtain the optical constants of the aerosols, but the corresponding experiments are complicated and may bring in inevitable errors. Moreover, the application of K-K relation usually requires the integral of all spectrums, which will increase the difficulty of the experiments and the requirement of equipment. To overcome the difficulties, a spectral transmittance and reflectance measurement technique combined with an improved artificial bee colony (IABC) algorithm is developed to retrieve the optical constants of aerosols directly. To improve the retrieval results, two improved inverse models, i.e. the double-concentration inverse model and the double-layer inverse model, are proposed. The remainder of the paper is organized as follows. First, the IABC algorithm is introduced and compared with the standard ABC algorithm in predicting the optical constants of aerosols. Then, based on the IABC algorithm, the optical constants of aerosols over Harbin, China at different wavelengths are retrieved under different random measurement errors. Moreover, two improved inverse models are proposed to improve the retrieval accuracy of the results. Finally, the optical constants of aerosols over Beijing, China available on AERONET are retrieved by the improved inverse models. 2. Direct problem
Fig. 1 Figure 1 depicts a beam of collimated monochromatic laser impinges on a slab filled with homogeneous, absorbing, and scattering aerosols under common temperature (emitting effect of aerosols ignored at room temperature). The light transfer process can be described by the one-dimensional steady-state radiative transfer equation (RTE) as follows [18, 19]:
s
I ( z, s ) I ( z, s ) z 4
I ( z, s )Φ (s , s )d 4
i
i
i
(1)
where I ( z, si ) is the radiative intensity in direction s at location z; denotes the extinction coefficient;
s, is the scattering coefficient; (si , s ) is the scattering phase function. The mathematical definitions of , s, and (si , s ) are derived as follows:
s, ( si , s )
1
n
4
D N
4
n
2 j
0
f ( D j )Qext, j
(2)
D N
f ( D j )Qsca, j
(3)
j 1 n
j 1
4D N j 1
2 j
2 j
0
0
f ( D j )Qsca, j p, j ( D j , si , s )
(4)
where N 0 is the total number concentration of the sample; f ( D j ) is the size distribution function when the aerosol diameter is D j ; Qext , Qsca and p denote the extinction efficiency factor, scattering efficiency factor and scattering phase function of single aerosol, respectively, which can be predicted by using the Mie theory [20]. The boundary conditions of the sample can be expressed as:
0 I , I (0, ) 0 0, 0 / 2 I ( L, ) 0, / 2
(5) (6)
where I 0 is the total incident light intensity; I (0, ) and I ( L, ) denote the light intensity incident to the internal medium from the left side and right side of the sample, respectively; denotes the angle between the m and positive x axis; L (in m) is the geometrical thickness of the slab. Solving the RTE by numerical methods
such as the Finite Volume Method (FVM) and Monte-Carlo method (MCM), the radiative intensities on the two borders of the sample I (0, ) and I ( L, ) can be obtained and the spectral reflectance R and transmittance can be derived as follows:
I (0, ) cos sin d /2 I0
R 2 T 2
/2
0
I ( L, ) cos sin d . I0
(7) (8)
Roughly, for most dry aerosols, the refractive index of the optical constant is about 1.53, with a small absorption index of 0.005-0.01. If the aerosols absorb water vapor and become hygroscopic particles, the value of the optical constant will be diluted and decrease towards the value of pure water with optical constant of 1.33+0.00 i. For instance, the oceanic aerosols, mainly composed of hygroscopic sea salt, almost have no absorption, and the refractive index of the optical constant is around 1.38, which is just slightly larger than that of pure water. Soot or black carbon has a very high absorption with a representative optical constant of 1.75+0.45 i, although black carbon is most often found internally mixed with other diluent particles, especially for other organic materials [21, 22]. In this study, the optical constants of aerosols measured over Harbin, China are studied. The corresponding optical constants and aerosol size distribution are shown in Fig. 2.
Fig. 2
Considering the RTE is a non-linear integral-differential equation, simulated transmittance and reflection signals, which contain useful information for retrieving the optical constants of aerosols, under different thickness of the slab L or different total number concentration N0 will be different and linear independence. So, how to chose L and N0 is of great significance to obtain more useful information to retrieve the optical constants of aerosol more accurately. To make this point, a detailed examination of the sensitivity analysis of transmittance or reflection to the optical constants at different L or different N0 is employed to find out the optimal thickness of the slab and different particle concentration. The sensitivity coefficients are one of the most important characteristic parameters in the sensitivity analysis, which is the first derivative of the transmittance or reflectance to optical constants in present study [21]:
a ( )
(a a, ) (a a, ) , N0 or L a 2a
(9)
where a denotes the independent variable which stands for refractive index n and absorption index k; represents a tiny change and is set 0.5% in present work; is the measurement radiative signals, i.e. reflection R or transmission T . Figure 3 depicts the sensitivity coefficients of n and k with respect to reflection and transmission signals under different N0 at different wavelengths. It is obvious that if the sensitivity coefficients are either small or correlated with one another, the estimation problem will be very sensitive to measurement errors and difficult to be solved. Therefore, the optimal N0 region for inverse problems is the region with larger absolute value of sensitivity coefficient generally, i.e. N0 2 109 , 2 1011 / m3 . Similarly, the optimal region of L for inverse problems can be obtained from Fig.4, i.e. L 0.01,0.2m . Fig. 3 Fig. 4
3. Inversion algorithm 3.1 The principle of standard ABC algorithm Bees are social insects and live in colonies. Their behavior is governed by the goal of colony survival rather than being focused on the survival of individuals. When looking for food, the bees can be classified into three categories: employed bees, onlooker bees and scout bees. The bees that are going to discover food sources are named employed bees. The ones that are waiting in the hive for making decision to choose a food source are called onlooker bees. The other kind of bees that carries out random search for finding new sources is named scout bees. The employed bees bring back the nectar information from the food source to the hive and share it with onlooker bees by dancing on the dance area in the hive. The longer the duration of dance is, the larger the nectar amount of the food source exploited is. Onlooker bees watch the dances of many employed bees, and finally tend to choose a food source according to the probability proportional to the quality of the food source.
Therefore, the food sources with abundant nectar will attract more onlooker bees. Whenever a food source is exhausted, this food source will be abandoned and the corresponding onlooker bee will become scout bee or onlooker bee. Scout bees can be treated as performing the job of exploration, while employed and onlooker bees can be seen as performing the job of exploitation [23-25]. Being inspired by foraged behavior of the bee colony, Karaboga proposed the ABC algorithm and employed it to solve the problem of multi-variable and multi-modal continuous function optimization firstly [26]. Then, the ABC algorithm has drawn widespread attention and been popularized to tackle some practical optimization applications all over the world for its satisfactory global optimization and fewer control parameters [23, 25, 26]. In standard ABC algorithm, each food source represents a possible solution to the concerned problem and the abundant degree of the corresponding nectar implies the quality of the solution represented by the fitness value. So, the bee colony usually divides to two equal parts. One half is the employed bees whose number is the same as the number of food source, which make sure exactly one employed bee for every food source. The other half is the onlooker bees. At first, the food source positions Xi [ xi ,1 , xi ,2 ,
, xi ,Nn ](i = 1, 2,
, SN ) are initialized randomly
in the search space, which contains the food source sites [27], xi , j lowj rand1 (high j low j )
where i 1, 2,
, SN and j 1, 2,
(10)
, Nn , SN denotes the number of food source and Nn denotes the
dimension of unsolved problem; xi , j denotes the i th food source position at the j th dimension problem; high j and low j is the upper and lower limits of the search space
for the j th dimension problem; rand1
is the uniformly distributed random number in the range of [0, 1]. Then, evaluate Xi and record current personal best food source position Pi and pbesti . After initialization, the employed bees are all sent out for food. The i th employed bee produces a candidate food source position Vi [vi ,1 , vi ,2 ,
, vi ,Nn ] from the old personal best
one Pi in its memory depending on local information to find a neighboring food source as follows [27]:
Vi Pi i (Pi Pk )
(11)
where i is a uniformly distributed real random number in the range of [0, 1]; k 1, 2,
, SN and k i .
Then, evaluate the quality of the candidate food source position and calculate the corresponding fitness value Fiti [28]: 1 / (1 Fi ) if Fi 0 Fiti 1 abs( Fi ) if Fi 0
(12)
where Fi is the cost value of the food source position Vi . For maximization problem, the cost function can be directly used as a fitness function. Comparing the fitness function value of Xi and that of Vi , the one with more satisfactory fitness value will be selected and the corresponding position will be memorized. If Xi cannot be improved, its counter holding the number of trials is incremented by 1, otherwise, the counter is reset to 0. When the employed bees bring back the nectar information of the food sources, the onlooker bees choose a food source depending on the probability Pi [28]:
Pi
Fiti SN
Fit i 1
.
(13)
i
In this probabilistic selection scheme, as the nectar amount increases, the number of onlooker bees visiting corresponding food source increases. If a position of food source cannot be improved further after a predetermined number limit of cycles, the food source will be abandoned, and the corresponding employed bee will be changed to onlooker bee or scout bee. Assume that the abandoned source is Xn , n 1, 2,
, SN , the
scout bee will discover a new food source as follow [23]: xn, j low j rand2 (high j low j )
(14)
where rand2 is the uniformly distributed random number in the range of [0, 1]. At the end of every generation, personal best position Pi and corresponding fitness value pbesti of each food source are recorded. Meanwhile, the global best position Pg and corresponding fitness value gbest of all the food sources are also recorded.
3.2. The theory of the IABC algorithm
In standard ABC algorithm, the fast information flow between bees seems to be the reason for clustering of bees, which will also lead to the diversity decline rapidly and leave the ABC algorithm with great difficulties of escaping local optima. Moreover, when solving those unimodal problems, the convergence speed of the ABC algorithm is typically slower than those of representative population-based algorithm (e.g. genetic algorithm and particle swarm optimization) [26]. Therefore, to remedy this problem, an improved ABC algorithm (IABC) applies a logistic mapping to generate a chaotic sequence in its initialization and introduces the Differential Evolution (DE) algorithm to improve the personal best food source positions. In the ABC algorithm, the initial food source positions are usually generated by random numbers that may allow some of the initial positions too close to each other which will increase the convergence time, see Eq. (9). The Chaos theory is proposed to overcome those difficulties in this paper. Chaos is the highly unstable motion of deterministic systems in finite phase space that often exists in nonlinear systems. The typical logistic mapping to generate the Chaos signal is described as follows [29, 30]: r (0) rand3 r (k 1) r (k ) 1 r (k ) , k 0,1,
, Nn 1
(15)
where rand3 is a uniformly distributed random number in the range of [0, 1]; is control number, when 4.0 , the logistic mapping is in a fully chaotic state. So, Eq. (9) can be expressed as: xi , j lowj r ( j ) (high j low j )
(16)
Moreover, the DE algorithm is applied to improve the personal best position Pi of each food source and avoid local optima. The DE algorithm, developed by Storn and Price [31], is an evolutionary optimization algorithm based on population cooperation and competition of individuals and has been successfully applied to solve optimization problems particularly involving non-smooth fitness function [31-33]. The basic optimization process in the DE algorithm consists of mutation, crossover and selection [33]. Specifically, the mutation operation generates a new parameter vectors by adding the weighted difference between two population vectors to a third vector. Then, the mutated vector parameters are mixed with the parameters of another predetermined vector,
named as the target vector, to yield the so-called trial vector. This operation is named as crossover. The last operation, called selection, is to judge whether the fitness value yielded by the trial vector is superior to that of the target vector. If so, the trial vector replaces the target vector in the following generation. Each population vector has to serve once as the target vector so that SN competitions take place in one generation [31]. The mutant vector Pi,mut can be generated according to [31]: Pi,mut Pk F P j Pm
(17)
where k , l , m are random integers uniformly selected from 1, SN , ik l m ; F , the mutant factor, is a real and constant value, F0, 2 , which controls the amplification of the differential variation P j Pm . In order to increase the diversity of the perturbed parameter vectors, crossover is introduced and the trial vector Pi,tri ( pi1,tri , pi 2,tri , pi3,tri ..., piNn ,tri ) can be derived from [31]:
pij ,mut pij ,tri (t ) pij
if rand 4 ( j ) CR or j rnbr (i ) if rand 4 ( j ) CR or j rnbr (i )
(18)
where rand4 ( j ) is the j th evaluation of a uniform random number generator, rand4 ( j )0, 1 ; CR is the crossover constant CR0, 1 ; rnbr (i) is a random integer uniformly selected from 1, SN . The search procedures of the IABC algorithm are shown below and the corresponding flowchart is illustrated in Fig. 5. Step 1.
Input system control parameters of the IABC algorithm, i.e. the total number of the food sources SN , the maximum number of iterations M , the number of the inversion parameters (dimensions) N n , the searching space
= [lowi , highi ](i = 1,2,
, Nn ) , the tolerance for minimizing the fitness value , the
predetermined number of cycles limit and the values of , F , CR . Step 2.
Initialize the positions of food sources in the N n dimensional problem by mapping the chaotic sequence, which is generated by Eq. (16), to the search space. Evaluate the fitness value of each food source position Fiti . Obtain the initial personal best food source position Pi and the global best food source position Pg , and record corresponding fitness values pbesti and gbest . Set the number of current iteration t 0 .
Step 3.
Determine whether the program matches one of the following two stop criterions. If so, end the calculation and record the estimated values; if not, proceed to Step 4. (i) The fitness value of global best food source position is less than the tolerance , gbest ; (ii) The number of the iteration reaches the user-defined iteration limit M , t M .
Step 4.
Send employed bees to food sources, every employed bee produces a modification near the position of former personal best food source position Pi to create a candidate food source position Vi by using Eq.(11). Evaluate Vi and compare the fitness values of Pi and Vi , and record the position with more satisfactory fitness value as the personal best food source position. Then, judge whether the position of the i th food source has been improved. If not, its counter holding the number of trials Tri(i) is incremented by 1; if so, the counter is reset to 0.
Step 5.
Employed bees bring back the nectar information of food sources, and share it with onlooker bees. Onlooker bees choose food sources with probabilities related to nectar amount of food sources according to Eq. (13).
Step 6.
Judge whether Tri(i) > limit , i = 1, 2,
, Nn . If so, abandon the i th food source position and
change the corresponding employed bee to scout one. The scout bee discovers new food source by using Eq. (14). Step 7.
Calculate the mutant vectors for current personal best positions of all food sources using Eq. (17). Operate the crossover and calculate the trial vector Pi ,tri using Eq. (18). Finally, evaluate the trial vector and compare the fitness value of the trial vector and that of the target vector (current personal best positions of food source). If the trial vector is superior, update the fitness value pbesti and corresponding personal best food source position Pi .
Step 8.
Determine whether the fitness values of the personal best food source positions pbesti are superior to those of the current global best food source position gbest . If so, update the value gbest and
corresponding location Pg . Loop to Step 3. Fig. 5 4. Inverse models and simulation The retrieval of the aerosol optical constants is solved by minimizing the fitness function value Fit , which is defined as the sum of the square residual between the estimated and measured signal ratios. The mathematical expression of Fit is derived as follows: 2 2 1 Rest Rmea Test Tmea Fit . 2 Rmea Tmea
(19)
where Rest and Test are the estimated transmittance and reflectance; Rmea and Tmea denote the real transmittance and reflectance. Considering that the algorithm is a stochastic optimization method and has certain randomness, all the calculations are repeated 50 times. The standard deviation of the inverse results is investigated to evaluate the reliability and feasibility of the inverse results, and the mathematical expressions are described as follows:
n
nest
1 50
50
nest nest, i , 2
1 50
k
i 1
1 50
50
i 1
nest, i ,
kest
1 50
50
k
est
kest, i
2
(20)
i 1
50
k
est, i
(i 1,
,50)
(21)
i 1
4.1 Comparison of standard ABC and IABC algorithms At first, the performance of the IABC is investigated through comparison with standard ABC. Figure 6 displays the comparison of the fitness function values between standard ABC and the IABC in retrieving the optical constants of aerosols over Harbin, China at wavelength =8.73 m (m=1.39655+0.14207 i). The geometric thickness of the measurement sample L is set as 0.01 m, and the total number concentration N0 is set as 2.0× 1010/m3. The aerosol size is assumed to vary from 0.001~10 m, which covers the main size range of aerosols [34, 35]. The corresponding system control parameters are listed in Table 1. Table 2 lists the corresponding retrieval results. The investigation shows that the fitness function value of the IABC converges much faster than the
standard ABC. Moreover, the IABC can arrive at the lowest best fitness value than the standard ABC within a smaller number of generations (see Fig. 6). All the results demonstrate that the IABC can avoid the phenomenon of local optima and low convergence accuracy which exit in the standard ABC.
Table 1 Table 2 Fig. 6
4.2 Retrieval of the optical constants of aerosols The optical constants of aerosols over Harbin, China at different wavelengths (between 6.0 m and 17.0 m) are estimated. The geometrical thickness of the slab L is set as 0.01 m, and the total number concentration N0 is set as 2.0×1010/m3. The system control parameters of the IABC are shown in Table 1. The searching ranges of the refractive and absorption indices are set as [1.0, 3.0] and [0.001, 1.0], respectively. Figures 7-9 depict the retrieval results of the optical constants with different random measurement errors added to the measured reflectance and transmittance. Figure 7 indicates that the inverse results of the optical constants show high accuracy without random measurement errors. The largest absolute errors of the refractive and absorption indices are at =6.68 m, which are less than 0.04 and 0.09, respectively. If the random measurement error is added, the retrieval accuracy will reduce and the standard deviation and absolute errors become larger (Figs. 8 and 9). When the random error is added up to 5%, the maximum absolute errors of the refractive and absorption indices increase to 0.87 and 0.17, respectively. Moreover, it is obvious that there is similar retrieval accuracy between the retrieval results of the refractive index and those of absorption index.
Fig. 7
Fig. 8
Fig. 9
4.3 The improved inverse models The light transmittance and reflectance measurement technique provides a potential method to retrieve the optical constants of aerosols. Without random errors, the estimated results are satisfactory. However, when there is a random error, the inverse results deteriorate (see Figs. 8 and 9). To remedy this situation, two improved inverse models, i.e. the double-concentration and the double-layer inverse models, are proposed to provide more useful radiative measurement information related to the optical constants to improve the retrieval accuracy. The original inverse model (single-concentration and single-layer inverse model) studied in Part 4.2 is named as Model 1.
Model 2: The double-concentration inverse model In this model, the geometrical thickness of the measurement sample remains unchanged, and the concentration of the aerosol dispersion medium is changed ( N0 =2.0 1010 / m3 and N0 =2.0 1011 / m3 ). The searching ranges of the real and imaginary parts are also set as [1.0, 3.0] and [0.001, 1.0], respectively. The fitness function Fit can be defined as follow: 2 2 2 2 Test,1 Tmea,1 Rest,2 Rmea,2 Test,2 Tmea,2 1 Rest， 1 Rmea,1 Fit . 4 Rmea,1 T R T mea,1 mea,2 mea,2
(22)
Model 3: The double-layer inverse model In this model, the concentration of the aerosol dispersion medium remains unchanged, and the geometrical thickness of the slab is changed ( L1 =0.1 m and L2 =0.2 m ). The searching ranges of the refractive index and are also set as [1.0, 3.0] and [0.001, 1.0], respectively. Thus, the fitness function Fit can also be defined as Eq. (22). Compared with Model 2, Model 3 needs two optical vessels with different geometrical thickness during actual experiement, which will increase the cost of experiment.
4.4 Comparison of the retrieval results under different inverse models Table 3 lists the estimated results of optical constants of aerosols at wavelengths =6.97 m and 11.07 m under the Model 1, Model 2 and Model 3, respectively. All calculations are repeated 50 times. It can be found that the average retrieval results of multiple calculations under Models 2 and 3 are more accurate than those under Model 1, and the standard errors under Models 2 and 3 are less than those under Model 1. The reason is that two improved inverse models provide more useful information about the optical properties of aerosols than Model 1, although it may be more time-consuming. Moreover, comparing the estimated results of two improved inverse models, it is easy to find that the accuracy of estimated results under Model 2 shows little better than that under Model 3, when there are the same random measurement errors.
Table 3
Figure 10 depicts the distribution of fitness function values at =6.97 m and 11.07 m under three different inverse models. From Fig. 10(a), it can be found that the minimum value region of fitness function is not a point but looks like a canyon, which indicates that the optimal solution is not unique and more than one appropriate value for absorption index are available to satisfy minimum fitness function under the same refractive index value. In other words, the multiple-value characteristics of the retrieval results exist and there are many couples of (n, k) values corresponding to the same experimental transmittance and reflectance. This is the difficulty of inverse estimated of aerosol optical constants, which is due to the complication of the scattering amplitude functions in the precise Mie theory. Moreover, comparing the distribution of the fitness function values under different inverse models, it is obvious that the minimum fitness function value regions under the improved inverse models are smaller than those under Model 1. This phenomenon indicates that corresponding to the same experimental transmittance and reflectance, less couples of (n, k) values may exist under the improved inverse models than
those under Model 1, which can explain why the retrieval results of the improved inverse models are better than those of Model 1. Fig. 10 4.5 Study of aerosols available on the AERONET Figure 11 illustrates the aerosol size distribution and the optical constants over Beijing, China available on the AERONET. The retrieval results under two improved inverse models are shown in Fig. 12. It can be found that the average estimated results of multiple calculations under the improved inverse models, especially Model 2, are more accurate than Model 1 in solving the problem of aerosol optical constants, even with 3% random measurement errors. Moreover, comparing the inverse results of refractive index and those of absorption index, it is obvious that the estimated accuracy of absorption index shows little better than that of refractive index. The phenomenon is attributed to the multi-value characteristics of the estimated results of refractive index, which can also be demonstrated by studying the distribution of fitness function values. For simplifying the manuscript, the study of fitness function values for the aerosols over Beijing is not available here, and interested readers can contact us for the details.
Fig. 11
Fig. 12 5. Conclusions Two improved inverse models, named double-concentration and double-layer inverse models, combined with an improved artificial bee colony (IABC) algorithm are proposed to improve the retrieval accuracy of aerosol optical constants. The results demonstrate that the improved inverse models combined with the IABC are suitable for retrieving aerosol optical constants. The following conclusions can be drawn: (1) Compared with standard ABC, the IABC can avoid the local optima and obtain higher convergence accuracy within a smaller number of iterations in solving the problem of retrieving aerosol optical constants.
(2) The retrieval results of two improved inverse models (i.e. the double-concentration inverse model and double-layer inverse model) are proved to be more accurate than the original inverse model, because more useful optical information about the aerosol particles can be obtained, and the effect of multi-value characteristics of the retrieval results can be weakened. Moreover, comparing the results of these two improved inverse models, the double-concentration inverse model shows more retrieval accurate than the double-layer inverse model. (3) Combined with the IABC, the optical constants of aerosols over Beijing, China obtained from AERONET can be estimated effectively under the improved inverse models, even with random measurement errors.
Acknowledgements The supports of this work by the Foundation of Nanjing University of Aeronautics and Astronautics of China for new teachers start work (No. 1022-YAH16051), the National Natural Science Foundation of China (No: 51606095) and Jiangsu Provincial Natural Science Foundation (No: BK20160794) are gratefully acknowledged. The author would also like to thank all the investigators of AERONET for their free provision of the aerosol properties over Beijing. A very special acknowledgement is made to the editors and referees who make important comments to improve this paper. References 1.
Dinar E, Riziq AA, Spindler C, Erlick C, Kiss G, Rudich Y. The complex refractive index of atmospheric and model humic-like substances (HULIS) retrieved by a cavity ring down aerosol spectrometer (CRD-AS). Faraday Discussions 2008; 137: 279-295.
2.
Ebert M, Weinbruch S, Hoffmann P, Ortner HM. The chemical composition and complex refractive index of rural and urban influenced aerosols determined by individual particle analysis. Atmospheric Environment 2004; 38: 6531–6545.
3.
Ebert M, Weinbruch S, Rausch A, Gorzawski G, Helas G, Hoffmann P, Wex H. Complex refractive index of aerosols during LACE 98#x2010; as derived from the analysis of individual particles. Journal of Geophysical Research Atmospheres 1999; 3(2): 275-280.
4.
Noia AD, Hasekamp OP, Van Harten G, Rietjens JHH, Smit JM, Snik F, Henzing JS, DeBoer J, Keller CU, Volten H. Use of neural networks in ground-based aerosol retrievals from multi-angle spectropolarimetric observations. Atmospheric Measurement Techniques 2015; 7: 281-299.
5.
Spindler C, Riziq AA, Rudich Y. Retrieval of Aerosol Complex Refractive Index by Combining Cavity Ring Down Aerosol Spectrometer Measurements with Full Size Distribution Information. Aerosol Science & Technology 2007; 41(11): 1011-1017.
6.
Flores JM, Washenfelder RA, Adler G, Lee HJ, Segev L, Laskin J, Laskin A, Nizkorodov SA, Brown SS, Rudich Y. Complex refractive indices in the near-ultraviolet spectral region of biogenic secondary organic aerosol aged with ammonia. Physical Chemistry Chemical Physics 2014; 16(22): 10629-10642.
7.
Guo LF, Shen JQ, Wan C. Particle analysis based on light scattering of particles illuminated by a divergent Gaussian beam. Optics and Lasers in Engineering 2013; 51: 826-831.
8.
Miles RE, Walker JS, Burnham DR, Reid JP. Retrieval of the complex refractive index of aerosol droplets from optical tweezers measurements. Physical Chemistry Chemical Physics Pccp 2012; 14(9): 3037-3047.
9.
Zarzana KJ, Cappa CD, Tolbert MA. Sensitivity of Aerosol Refractive Index Retrievals Using Optical Spectroscopy. Aerosol Science & Technology 2014; 48(11): 1133-1144.
10. Li JY, Dong SK, Tan HP. Effective Optical Constant of Alumina Particle Containing Carbon. Journal of Thermophysics & Heat Transfer 2009; 23(1): 216-219. 11. Ruan LM, Wang XY, Qi H, Wang SG. Experimental investigation on optical constants of aerosol particles. Journal of aerosol science 2011; 42(11): 759-770. 12. Ku JC, Felske JD. Determination of refractive indices of Mie scatterers from Kramers-Kronig analysis of spectral extinction data. Journal of the Optical Society of America A 1986; 3(5): 617-623. 13. Herman BM, Browning RS, De Luisi JJ. Determination of the Effective Imaginary Term of the Complex Refractive Index of Atmospheric Dust by Remote Sensing: The Diffuse-Direct Radiation Method. Journal of
the Atmospheric Sciences 1975; 32(5): 918-925. 14. Wei DJ, Qiu JH. Wide-Band Method to Retrieve the Imaginary Part of Complex Refractive Index of Atmospheric Aerosols. Part I: Theory, Scientia Atmospherica Sinica, 1998. 15. Dombrovsky LA, Sazhin SS, Mikhalovsky SV, Wood R, Heikal MR. Spectral properties of diesel fuel droplets Fuel 2003; 82(1): 15-22. 16. Ruan LM, Qi H, An W, Tan HP. Inverse radiation problem for determination of optical constants of fly-ash particles. International Journal of Thermophysics 2007; 28(4):1322-1341. 17. Zhao JJ, Duan YY, Wang XD, Zhang XR, Han YH, Gao YB, Lv ZH, Yu HT, Wang BX. Optical and radiative properties of infrared opacifier particles loaded in silica aerogels for high temperature thermal insulation. International Journal of Thermal Sciences 2013; 70(8): 54-64. 18. Modest MF. Radiative heat transfer, McGraw-Hill New York 2013. 19. Ren YT, Qi H, Chen Q, Ruan LM, Tan HP. Simultaneous retrieval of the complex refractive index and particle size distribution. Optics Express 2015; 23(15): 19328-19337. 20. Bohren CF, Huffman DR. Absorption and scattering of light by small particles. Wiley. Com 2008. 21. Chudzyński S, Czyżewski A, Ernst K, Karasiński G, Kolacz K, Pietruczuk A, Skubiszak W, Stacewicz T, Stelmaszczyk K, Szymański A. Multiwavelength lidar for measurements of atmospheric aerosol. Optics and Lasers in Engineering 2002; 37: 91-99. 22. Lenoble J, Remer L, Tanre D. Aerosol Remote Sensing. Springer 2013. 23. Gao W, Liu S, Huang L. A global best artificial bee colony algorithm for global optimization. Journal of Computational & Applied Mathematics 2012; 236(11): 2741–2753. 24. Akay B., Karaboga D. Artificial bee colony algorithm for large-scale problems and engineering design optimization. Journal of Intelligent Manufacturing 2012; 21(21): 1-14. 25. Singh A. An artificial bee colony algorithm for the leaf-constrained minimum spanning tree problem. Applied
Soft Computing 2009; 9(2): 625-631. 26. Karaboga D, Akay B. A comparative study of Artificial Bee Colony algorithm. Applied Mathematics & Computation 2009; 214(1): 108-132. 27. Akay B, Karaboga D. A modified Artificial Bee Colony algorithm for real-parameter optimization. Information Sciences 2012; 192(6): 120-142. 28. Maghali I, Oliveira SD, Schirru R, Medeiros JACCD. On the Performance of an Artificial Bee Colony Optimization Algorithm Applied to the Accident Diagnosis in a PWR Nuclear Power Plant. 2009 International Nuclear Atlantic Conference (INAC2009), Sep. 27th-Oct. 2nd, 2009, Rio de Janeiro, Brazil. 29. Liu F, Duan H, Deng Y. A chaotic quantum-behaved particle swarm optimization based on lateral inhibition for image matching. Optik-International Journal for Light and Electron Optics 2012; 123(21): 1955-1960. 30. Zhang B, Qi H, Sun SC, Ruan LM, Tan HP. Solving inverse problems of radiative heat transfer and phase change in semitransparent medium by using Improved Quantum Particle Swarm Optimization. International Journal of Heat and Mass Transfer 2015; 85: 300-310. 31. Storn R, Price K. Differential evolution–a simple and efficient heuristic for global optimization over continuous spaces. Journal of global optimization 1997; 11(4): 341-359. 32. Wu Y, Lu J, Sun Y. An improved differential evolution for optimization of chemical process. Chinese Journal of Chemical Engineering 2008; 16(2): 228-234. 33. Lu S, Sun C, Lu Z. An improved quantum-behaved particle swarm optimization method for short-term combined economic emission hydrothermal scheduling. Energy Conversion and Management 2010; 51(3): 561-571. 34. Yuan Y, Yi HL, Shuai Y, Liu B, Tan HP. Inverse problem for aerosol particle size distribution using SPSO associated with multi-lognormal distribution model. Atmospheric Environment 2011; 45(28): 4892-4897. 35. Tang H, Lin JZ. Retrieval of spheroid particle size distribution from spectral extinction data in the
independent mode using PCA approach. Journal of Quantitative Spectroscopy and Radiative Transfer 2013; 115: 78-92. 36. Aerosol robotic network (AERONET) 2014, NASA,
Figure Captions Fig. 1
The schematic model of the parallel slab with collimated laser incident
Fig. 2
The size distribution and optical constants of aerosols over Harbin, China
Fig. 3
Sensitivity coefficients of n and k with respect to reflection or transmission signals under different N0 at different wavelengths
Fig. 4
Sensitivity coefficients of n and k with respect to reflection or transmission signals under different L at different wavelengths
Fig. 5
Flowchart of the IABC algorithm
Fig. 6
Comparison of fitness function values of standard ABC and IABC algorithms
Fig. 7
Reproducibility of optical constants of aerosols without random error.
Fig. 8
Reproducibility of optical constants of aerosols with 3% random error.
Fig. 9
Reproducibility of optical constants of aerosols with 5% random error.
Fig. 10
The distribution of the fitness function values under different inverse models with different optical constants when 1.0
Fig. 11
The size distribution and optical constants of aerosols over Beijing available on AERONET [36]
Fig. 12
The retrieval results of optical constants of aerosols over Beijing under different ramdon measurement errors
Aerosol sample Incident Laser
Reflection signals Fig. 1
Transmittance signals
10 Refractive index Absorption index
3.0 2.5 2.0 1.5 1.0
-1
1.5
10
1.4
10
-2
0.5 1.3
0.0 1
10
Diameter (m)
(a) Aerosol size distribution Fig. 2
100
-3
6
7
8
9
10 11 12 13 14 Wavelength, (m)
(b) Optical constants
15
16
10 17
Absorption index, k
3.5
0
1.6
Harbin, China, on 6th Mar. 2014 True distribution
Refractive index, n
Volume frequency distribution, f(D)
4.0
0.035
0.010
(b)
(a)
0.030 0.025 0.020
n T )
nR)
0.005
L=0.01m -3 =2.58m, n=1.40, k=1.14x10 =8.73m, n=1.39, k=0.142 =10.48m, n=1.47, k=0.137 =12.53m, n=1.46, k=0.139
0.015 0.010 0.005
0.000
0.000 -0.005 -0.005 8 10
10
9
10
10 3 N0(1/m )
10
11
10
12
0.005
(d)
(c) 0.000
10
10 3 N0(1/m )
11
10
11
10
10
12
0.05 0.04
0.01 k T )
kR)
9
0.02
-0.010 -0.015
Fig. 3
10
0.03
-0.005
-2 -4 -6 -8 -10 -12 -14 -16 -18 8 10
-0.010 8 10
0.00
10
9
10
10 3 N0(1/m )
10
11
10
12
-2 -4 -6 -8 -10 -12 8 10
10
9
10
10 3 N0(1/m )
10
12
0.010 (a)
0.020 (b) 0.015 0.010
10
3
=2.58m, n=1.40, k=1.14x10 =8.73m, n=1.39, k=0.142 =10.48m, n=1.47, k=0.137 =12.53m, n=1.46, k=0.139
-3
n T )
nR)
0.005
N0=1x10 /m
0.005 0.000
0.000
-0.005 -0.005 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m) 0.005 (c)
-0.010 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m) 0.035
0.030
0.000
0.025 0.020
-0.010
k T )
kR)
-0.005
0.015 0.010 0.005
-0.015 -10 -12 -14 -16 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m)
Fig. 4
0.000 -4 -8 -12 0.01 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 L(m)
Start Input system datum and control parameters
Initialization t=0
Y
gbest ?
N t=t+1 i=0 i=i+1 Evaluate candidate food position Vi and calculate fitness value Fiti
N
Fiti pbesti ? Y Update the personal best position Pi Vi
Onlooker bees choose food sources N
Tri(i ) limit ? Y Abandon the ith position, discover new position
Evaluate Pi ,tri and calculate the fitness value Fiti ,tri N
Fiti ,tri pbesti ?
Y Update the personal best position Pi Pi ,tri i SN ? N
N
Y
pbesti gbest
Y Update the global best position Pg Pi tM ?
N Output result End
Fig. 5
Y
-1
10
ABC IABC
-3
10
-5
10
-7
Fit
10
-9
10
-11
10
-13
10
-15
10
-17
10
0
200
400 600 Interation
800
Fig. 6
28
1000
-3
10
-5
10
1.3
-7
10 1.2
-9
10 9
10 11 12 13 Wavelength, (m)
14
15
16
10
0
10 -1
-2
10
10
-4
10 -2
10
-6
10
-11
8
2
10
-3
10 17
10
-8
6
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
10 17
Fig. 7
10
1.4
10
1.2
10
1.0
10
0
-1
-2
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
Absolute error
0
10
3
-1
10
-2
10
-3
10
10
1
10
-1
10
-3
10
-3
6
10 Real value Inverse results, Err=3%
-4
10 17
10
Absolute error of inverse results
1
1.6
5
10
10
Absorption Index, k
Absolute error
Absolute error of inverse results
Refractive Index, n
Real values Inverse results, Err=3%
0.8
1
2
1.8
-5
6
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
10 17
Fig. 8
2
10
Refractive Index, n
1.8
1
10
1.5 0
10
1.2 -1
10
0.9
0.6
1
10
-2
10
7
8
9
10 11 12 13 Wavelength, (m)
14
15
16
10 Real value Inverse results, Err=5%
Absolute error
0
10
2
-1
10
-2
10
-3
10
10
0
10
-2
10
-4
10
-3
6
4
10
Absorption Index, k
Aboslute error
Aboslute error of inverse results
Real values Inverse results, Err=5%
-4
10 17
10
Fig. 9
29
-6
6
7
8
9
10 11 12 13 14 Wavelength, (m)
15
16
10 17
Absolute error of inverse results
3
2.1
Absolute error of inverse results
1.4
Absolute error
0
Absorption index, k
Refractive index, n
-1
10
7
10 Real value Inverse results, Err=0%
1
10
1.5
6
4
10
10
Absolute error
Absolute error of inverse results
Real values Inverse results, Err=0%
1.1
1
3
1.6
(a)
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
n
n
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
1.5
1 0
0.2
0.4
k
0.6
0.8
1 0
1
3
0.2
0.4
k
0.6
0.8
1
3
(c)
(d) Fitness
Fitness 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2.5
2
n
n
2.5
1.5
2
1.5
1 0
0.2
0.4
k
0.6
0.8
1 0 3
1
3
0.2
0.4
k
0.6
0.8
1
(f)
(e) Fitness 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
2
1.5
1.5
1 0
Fitness
2.5
2.5
n
1
1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -3.5 -4 -4.5 -5 -5.5 -6 -6.5
Fitness
2.5
1.5
Fitness
(b)
Fitness
2.5
n
8
3
3
0.2
0.4
0.6 k =6.97m, n=1.39, k=0.08284
0.8
1
Fig. 10
30
1 0
0.2
0.4
0.6 k =11.07m, n=1.47, k=0.131
0.8
1
Volume size distribution, N0f(D)
x10-2 0.9
Beijing, China at 22:51:24 on 18 May 2014 k n
Beijing, China at 22:51:24 on 18 May 2014
Refractive index, n
0.20
0.15
0.10
1.45
0.6
1.40
0.3
0.05
1.35 0.4
0.00 0.1
1 Radius (m)
(a) aerosol size distribution Fig. 11
10
0.5
0.6
(b) optical constants
31
0.7 0.8 0.9 Wavelength, (m)
1.0
0.0 1.1
Absorption Index, k
1.50
0.25
1.52 Real values Inverse results, Model 2, Err=0% Inverse results, Model 3, Err=0%
x10-3 Real values Inverse results, Model 2, Err=0% Inverse results, Model 3, Err=0%
8.0 7.5
Absorption Index, k
Refractive index, n
1.48
8.5
1.44
1.40
7.0 6.5 6.0 5.5 5.0 4.5
1.36 0.4
0.5
0.6
0.7 0.8 0.9 Wavelength, (m)
1.0
4.0 0.4
1.1
16
2.4 Real values Inverse results, Model 2, Err=5% Inverse results, Model 3, Err=5%
0.6
0.7 0.8 0.9 Wavelength, (m)
1.0
1.1
x10-3 Real values Inverse results, Model 2, Err=5% Inverse results, Model 3, Err=5%
14 12
2.0
Absorption Index, k
Refractive index, n
2.2
0.5
10
1.8 1.6 1.4 1.2 1.0 0.4
8 6 4 2 0
0.5
0.6
0.7 0.8 0.9 Wavelength, (m)
1.0
1.1
Fig. 12
32
0.4
0.5
0.6 0.7 0.8 Wavelength, (m)
0.9
1.0
1.1
Table 1
The system control parameters of the standard ABC and IABC algorithms
Parameters ABC IABC
SN 50 50
Nn 2 2
M
0
limit
F
CR
1000
10-16
100
—
—
—
1000
10-16
100
4.0
0.5
0.4
33
Table 2
The retrieval results of the optical constants by standard ABC and IABC algorithms
Errors Algorithms
IABC
Parameters
0%
3%
5%
n
1.39655
1.36343
1.35979
n
4.2×10-8
0.01389
0.01043
k
0.14207
0.13602
0.11252
k
2.6×10-7
0.04484
0.05649
n
1.39654
1.35970
1.34258
n
4.9×10-6
0.01877
0.03852
k
0.14196
0.12181
0.10754
k
1.7×10-4
0.08285
0.08137
ABC
34
Table 3
The retrieval results of the optical constants by different inverse models
Models
Model 1
=6.97m, n=1.39, k=0.08284
=11.07m, n=1.47, k=0.131
Err=0%
Err=3%
Err=5%
Err=0%
Err=3%
Err=5%
n
1.40866
1.16026
1.56586
1.47345
1.34190
1.30214
n
4.8×10-7
0.19485
0.24662
1.8×10-7
0.16306
0.29876
k
0.11098
0.24095
0.35821
0.13100
0.19843
0.28753
k
1.4×10-6
0.11173
0.24378
2.8×10-7
0.12599
0.34621
n
1.39377
1.32328
1.30076
1.47344
1.41366
1.35234
n
2.4×10-7
0.15774
0.19371
3.5×10-8
0.06476
0.18762
k
0.08285
0.13264
0.14454
0.13100
0.16147
0.20629
k
2.0×10-7
0.06512
0.08754
2.5×10-9
0.11835
0.32768
n
1.39377
1.28401
1.23445
1.47344
1.36311
1.31030
n
2.3×10-7
0.15340
0.20131
1.5×10-7
0.21500
0.11703
k
0.08285
0.16095
0.21629
0.13100
0.18349
0.26876
k
1.8×10-7
0.26462
2.2×10-7
0.12411
0.15248
Parameters
Model 2
Model 3 0.11608
35