Repetitive Genetic Inversion of Optical Extinction Data

 

B. R. Lienert, J.N. Porter and S.K. Sharma

 

Hawaii Institute of Geophysics and Planetology,

2525 Correa Rd, Honolulu HI 96822, U.S.A.

 

submitted to Applied Optics, 10-1-00

 

Abstract

 

We describe a genetic method of deriving aerosol size distributions from multi-wavelength extinction measurements. The genetic inversion searches for lognormal size distribution parameters whose calculated extinctions best fit the data.  By repetitively applying the genetic inversion using different random number seeds, we are able to generate multiple solutions that fit the data equally well. When these solutions are similar, they lend confidence to an interpretation, while when they vary widely, they demonstrate the non-uniqueness. In this way, we show that even in the case of a single lognormal distribution, many different distributions can fit the same set of extinction data unless the misfit is reduced below typical measurement error levels. In the case of a bimodal distribution, we find many dissimilar size distributions that fit the data to within 1% at six wavelengths. In order to satisfactorily recover the original bimodal distribution, we found extinctions at ten wavelengths must be fitted to within 0.5%.  Our results imply that many size distributions recovered from existing extinction measurements can be highly non-unique and should be treated with caution.

 

Introduction

 

Analytical (calculus-based) inversions of optical extinction measurements for aerosol size distributions are ill-posed due to the oscillatory nature of the Mie kernel1, as well as the kernel's rapid decrease when the aerosol diameters become less than a few wavelengths in size. Such inversions therefore require a priori assumptions about the nature of the aerosol size distribution (regularization) that can take many forms2-4. Even when a solution is successfully obtained, its uniqueness is still questionable due to the non-linearity of the problem and the local nature of linearized solutions5. One widely used analytical method6 iteratively applies corrections to a Junge Distribution7 having the form

where dN/dD is the differential concentration per unit volume of non-absorbing spherical particles in the diameter range D to D+dD and A and u are constants.  Constraints are then placed on the smoothness of the resulting size distribution in order to obtain an iterative solution. Equation 1 is physically unrealistic in that dN/dD continues to increase without limit as D à 0. A recent modification to this method8 has addressed some of these problems by using a more realistic starting distribution, along with additional measurements of the angular dependence of the optical extinction. However, the method8 still assumes that the problem is uniquely formulated after regularization.

 

We decided to address the limitations of analytical methods by using random search techniques, which make no assumptions about the problem's uniqueness. Their only drawback is the amount of time taken to explore the parameter search space. The ideal approach would be to explore the data fit of all possible combinations of the size distribution parameters within realistic ranges. However, the computer time to explore even a six-parameter problem in this way would be prohibitive, even using only 256 values for each parameter.  A practical approach is the Monte Carlo method9, which explores random selections of the parameters. It has been used to explore many non-unique problems10-12. Even so, the computer time required for Monte Carlo inversion has resulted in methods for further limiting the parameter search space such as genetic inversion13-15 and simulated annealing16,17. We chose genetic inversion mainly because an algorithm was publicly available5 that was easily incorporated into our existing Fortran code.

 

Our goal in this paper is to test the genetic inversion technique with extinction values generated from realistic aerosol size distributions, as well as with measured data. We then attempt to recover size distributions using genetic inversions of both these generated and measured data.  By comparing the recovered distributions for the generated data with their starting distributions, we can test for non-uniqueness.

 

Method

 

Measurements of aerosol size distributions18-20 using optical particle counters (0.1<D<10m) and electrostatic classifiers (0.01<D<0.1m) have indicated that real aerosol distributions are reasonably well represented by a combination of three lognormal distributions (nucleation, accumulation, and coarse modes), i.e.,

where Ni is the integrated number concentration,  is the median diameter and si is the geometric standard deviation in each distribution. Furthermore, only the two largest aerosol modes (accumulation and coarse) significantly affect the optical scattering in the range 0.35-1.06m20.

We will assume that such a multi-mode lognormal distribution adequately describes common atmospheric aerosol sizes.  Although we only treat the case of lognormal distributions, the approach described here could be used with any type of mathematical distribution function.

 

The wavelength dependence of the aerosol optical extinction on aerosol size distributions is calculated using the integral

aa(l,m) is the aerosol extinction, Qs(D,l,m) is the Mie extinction cross section, l is the wavelength and m is the aerosol refractive index. The integration in eq. 3 was performed using numerical quadrature and a Mie Fortran routine21 to calculate the extinction cross sections, Qs(D,l,m). To ensure adequate representation of observed distributions, a diameter range 0.002m < r < 50m was used for the calculations. The quadrature integration was tested by calculating results for a gamma distribution and comparing them with tabulated values22.

The number of quadrature points was determined empirically by increasing it until five figures of accuracy were achieved in the tabulated gamma function values.  It was found that 1500 equally spaced logarithmic points were sufficient to achieve this nominal accuracy. In order to speed up the calculations, the Mie cross-sections at each diameter value and wavelength were stored in an array on the first pass of the program. Subsequent integrations could then be performed using the stored Mie extinctions, as the integration ordinates remained constant regardless of the lognormal parameters. This latter technique was found to increase the calculation speed by a factor of about 30.

 

Although a graphical method has been given23 to determine one possible solution of a unimodal lognormal distribution, our own attempts to find a solution to the bimodal lognormal problem using linearized inversion24 were unsuccessful due to the high degree of non-uniqueness25. We hoped to find a method that would successfully obtain families of both unimodal and bimodal lognormal solutions that would fit the data.  Genetic inversion finds a best-fit solution to an inverse problem by randomly selecting n sets of parameter values (individuals) to provide a population of solutions. Each individual is formed by concatenating the binary numbers representing each of its parameter values (alleles) into a single string of ones and zeros (a chromosome) having length L. The squared sum of data residuals for each individual is then evaluated to provide a measure of fitness. We use the equation

to calculate the fitness, F, where sd is the average experimental error expected in N observations and Ddi is the ith data residual. A fitness of one then implies that the root mean square (r.m.s.) data residual is equal to the experimental error. Solutions are then selected from the initial population based on a fitness rating, e.g., a solution having a fitness of 0.5 is five times more likely to be selected than one having a fitness of 0.1. The chromosomes of the solutions with the highest fitness (parents. are then crossed and mutated to produce new solutions (children). Once a new child population is generated, the procedure is iterated using new generations of solutions until the desired (or maximum. fitness level is attained. We adapted the genetic inversion routine of Carroll5 to invert for the aerosol size distributions. Use of the program requires setting a number of options and parameters that are briefly described below.

 

1. Micro GA. This option starts with very small random populations that converge in a few generations of genetic inversion. The results of each of these “micro" inversions are then used for the main genetic inversion.

2. Population Size. The total number of solutions used in each iteration of the genetic inversion.

3. Mutation Fraction. The probability that each binary element of a child chromosome will be inverted.

4. Crossover Probability. The probability that the second parent's chromosome will be crossed with the first's in the formation of a child chromosome.

5. Tournament Selection. If this is one, random pairs are selected from the population and the most fit of each pair is used as a parent.

6. Elitism. If this is one, the fittest individual is always reproduced in a new generation by giving its chromosome to a randomly selected child.

7. Parameter Creep. The probability that each parameter value will be incremented by one parameter change unit (i.e., the total parameter range divided by the number of parameter values).

 

8. Uniform Selection. If this is zero, a random dividing point, p, is chosen in the two parent's chromosomes. If the crossover probability is satisfied, the two chromosomes are crossed by concatenating bits 1-p of the first chromosome with bits p+1-L of the second. If uniform selection is set to one, each of the L bits of the first parent's chromosome is separately replaced with the same bit of the other parent's chromosome according to the crossover probability, which is evaluated separately for each bit.

9. Niche Selection. This option implements a multidimensional phenotypic sharing scheme13 whereby similar individuals fitnesses are reduced according to a triangular sharing function of parameter distance.

10. Number of Children. The number of children produced by each pair of parents.

 

The settings we used for all the inversions described here are given in Table 1.

We found that the most critical settings for the genetic inversions were the total numbers and ranges of the lognormal parameter values. We experimented with different population sizes and empirically determined that 30 provided the most rapid convergence for the initial unimodal inversion, whereas 20 performed more efficiently for the bimodal inversions. We also empirically determined that searching for 256 (28) parameter values produced the most rapid convergence. We should emphasize that 256 values were sufficient to obtain convergence in all cases, implying that searching the parameter space in smaller increments was unwarranted by the constraints provided by the data. In order to keep their relative precision constant, the parameters were spaced at equal logarithmic, rather than linear intervals8. We also found it useful to repeat the genetic inversion multiple times using different seeds for the random number generator used to make the parameter selections. As well as overcoming non-convergence in a number of cases, this technique proved extremely useful for exploring the uniqueness of the solutions, as it provided specific examples of different solutions that fit the data. Although such examples can sometimes be obtained by examining the results for different generations having close to the best fitness, we found that there was a tendency of the fitness to increase by a factor of 2-3 in a single generation, preventing similar fitness solutions from being compared.

 

Results

Text Box:   Figure 2. Extinctions calculated for the models shown in Fig. 1. The black dots are the input data used for the genetic inversions. The r.m.s. extinction misfits are all <=1x10-5 m-1.

Text Box:   
Figure 1. Size distribution model (heavy black curve. and four models (light curves. obtained by genetic inversions of the extinctions predicted by Mie theory at four wavelengths for this model. Units for the size distribution concentrations are area/unit volume (m-1).

The first test case is a unimodal distribution (heavy solid curve in Fig. 1) that approximately represents the accumulation mode of city pollution20 corrected for hygroscopic growth at a relative humidity of 75%. We generated extinction coefficients at the four harmonic wavelengths of an Nd:YAG laser (0.266, 0.355, 0.532 and 1.064m) using a real refractive index of 1.38, which is a reasonable value for hygroscopic aerosols20.  Genetic inversions of this data using four different random number seeds resulted in the four size distributions shown in Fig. 1. The fitness for each inversion corresponds to an r.m.s. misfit of 1x10-5 m-1 (~ 10% misfit). The parameter ranges used for the unimodal inversion were: N1=100-100,000/cc; D1=0.1-2mum and s1=3-12m. The extinctions calculated for the four inverted distributions are shown in Fig. 2, along with the four initial data points. Fig. 3 shows the change in fitness with generation for the first of the inversions in Fig. 1 for different population sizes. It is clear that a population size of 30 provides the best performance, particularly when the calculation time (proportional to generation number times population size) is taken into account. It is apparent that although the larger diameters of the inverted distributions in Fig. 2 appear to be reasonably well constrained, the models exhibit a large degree of ambiguity at diameters smaller than 0.3m. We then repeated the inversions to Text Box:   
Figure 4. Size distribution models obtained using the same data as those in Fig. 1, but with the r.m.s. extinction misfit reduced to  2x10-6 m-1

Text Box:   
Figure 3. Changes in fitness with generation for different population sizes obtained using genetic inversion of the data in Fig. 2.

achieve an extinction misfit of less than 2x10-6 m-1 (~ 2% misfit). The results (Fig. 4) indicate that at this level of misfit, the unimodal distribution is well constrained by the genetic inversions.

 

Text Box:  
Figure 5. Bimodal size distribution models(light curves).obtained by genetic inversion of synthetic data at six wavelengths generated using the heavy black curve model.

Text Box:   Figure 6. Extinctions calculated using the models in Fig. 5. The black dots are the input data used for the genetic inversions.

Size distribution models obtained using the same data as those in Fig. 1, but with the r.m.s. extinction misfit reduced to 2x10-6 m-1.We then generated extinction coefficients for the bimodal distribution shown in Fig. 5. This bimodal distribution represents a marine environment under light wind conditions20, corrected for hygroscopic growth at a relative humidity of 75%. In order to constrain the six lognormal parameters, we used extinctions at six wavelengths. These were 0.355, 0.532 and 1.064m (the first three harmonics of an Nd:YAG laser), and 0.6, 0.8 and 0.9 m (three that could be obtained using variable-wavelength lasers). The parameter ranges used for the coarse mode peak in the bimodal inversions were:  N2=0.1-100/cc; D2=0.7-15m and s_2=3-12m. The four genetically inverted size distributions (with a misfit of 2x10-6 m-1) are also shown in Fig 5. Fig. 6 shows the extinctions for the initial and inverted size distributions. The inverted distributions show a large amount of variation, dramatically demonstrating the non-uniqueness inherent retrieving a bimodal size distribution from extinction measurements in the wavelength range 0.35-1.06m.

Text Box:   Figure 8. Measured sun photometer optical thicknesses converted to extinction assuming a 2 km thick aerosol layer (black squares). The curves are theoretical extinctions for the models in Fig. 9.

Text Box:   Figure 7. Size distribution models (light curves) obtained by genetic inversions of synthetic data generated using the solid black curve model at 10 wavelengths from 0.26 to 1.54m. The r.m.s. data misfits for these models are all < 5x10-7 m-1.

 

Text Box:  Figure 10.  Unimodal size distribution models obtained by repeated genetic inversion of the extinction data in Fig. 8. The data misfits are all <2.5x10-6m-1.

Text Box:   Figure 9. Unimodal size distribution models obtained by repeated genetic inversion of the extinction data in Fig. 8. The data misfits are all <5x10-6 m-1.How good would the data need to be to reasonably recover the bimodal distribution in Fig. 5? It is apparent that data at both lower and higher wavelengths are necessary. We therefore generated synthetic data at four additional wavelengths, 0.266, 0.45, 0.8 and 1.54m and repeated the inversions, still without success.  The data misfit was then reduced by another factor of four (to 5x10-7 m-1 or 0.5%). In order to achieve the four low misfits in Fig. 7, twenty genetic inversions (using different random numbers seeds) of up to 500 generations were necessary. The resulting models are shown in Fig. 7 and can be considered reasonable recoveries of the original bimodal distribution.  These studies illustrate the importance of reducing the error in the extinction measurements to the maximum degree possible.

 

As a final test, we applied the genetic inversion to the mean of four sets of sun photometer measurements collected on the Hawaii Ocean Time-series cruise26 (HOT89).  Fig. 8 shows the measured sun photometer optical thicknesses, which have been converted to extinctions by assuming that all the aerosols are in a 2000 m thick uniform layer.  These sun photometer measurements were collected on 1-12-98, 100 km north of Oahu. During this time period SW (Kona. winds brought volcanic aerosol up from the Kilauea volcano. This distribution pattern of volcanic aerosol was consistent with the wind fields as well as aerosol optical depths derived from the SeaWifs satellite. The volcanic aerosols are somewhat unique in that they are composed of sulfuric acid and water and have negligible absorption based on volatility tests. In situ measurements have shown that the volcanic aerosol size distribution diameter peaks near 0.5m at ambient relative humidities20. 

 

 Fig. 9 shows the size distributions (dashed curves) obtained by the genetic inversion.  The corresponding dashed curves in Fig. 8 are the extinctions calculated from these size distributions. Based on the ambient relative humidity, we assumed a value of 1.35 for the refractive index20. The extinctions were fitted to within 5x10-6 m-1, r.m.s., which corresponds to the error typically claimed for sun photometer optical depths of +/-0.0127.  In our own cross comparisons between well calibrated sun photometers we have found agreement can be better26 (<+/-0.005). We therefore tried improving the r.m.s. fit of the genetic inversion to 2.5x10-6 m-1 (Fig. 10). It is clear that improving the fit has a dramatic effect on the range of size distributions, particularly at the smaller sizes.  The family of size distributions in Fig. 10 are more physically reasonable, based on in situ measurements of marine aerosols20.

 

 

Discussion

 

We have demonstrated that it is possible to invert for a unimodal lognormal size distribution model using data at four wavelengths, provided that the total misfit in the measured extinctions does not exceed  2x10-6 m-1. This implies that the measurement errors in real extinction measurements would have to be of the same range. For a measured set of six optical thicknesses, we were able to show (Fig. 4. that reasonable constraints can be placed on a unimodal size distribution provided the data is fitted to within +/-0.005. Successful inversion of a bimodal lognormal distribution requires data in a wider wavelength range and of significantly greater accuracy. The repetitive genetic inversion technique presented here provides specific examples of the high degree of non-uniqueness that is present when the errors in the measured extinctions exceed 2x10-6 m-1. This lack of uniqueness, which is well-known in geophysical problems26,27 has been pointed out previously28.  As mentioned initially, the non-uniqueness is in large part due to the oscillatory variation of the Mie scattering cross section with particle diameter, which results in particles of widely varying sizes having the same scattering cross section at a given wavelength. It is also due, to a more limited extent, to the poor constraints on the smallest size aerosol using presently available wavelengths.

 

We have not addressed the problem of additional unknowns, such as refractive index/absorption, non-spherical particles, multiple scattering and lidar ratio (for lidar measurements). It is relatively trivial to solve for these using genetic inversion providing that they can be mathematically modeled. However, since we have demonstrated that the inversions are non-unique even for ``perfect" data, where these unknowns are well constrained, it is evident that their introduction can only make this problem even more non-unique. Additional constraints are essential in order to determine these additional unknowns8.

 

Acknowledgements

 

We thank the three anonymous reviewers for their suggestions. This work was supported by the Office of Naval Research under grant #N00014-96-1-0317. This is School of Ocean and Earth Science and Technology contribution no. 5347 and HIGP contribution no. 1137.

 

References

 

1. D. Muller, U. Wandinger and A. Ansmann, ``Microphysical particle parameters from extinction and backscatter data by inversion with regularization: theory", Appl. Optics, 38, 2346-2357 (1999).

2. G.P. Box, K.M. Sealy and M.A. Box, ``Inversion of Mie extinction measurements using analytic eigenfunction theory", J. Atmos. Sci., 49, 2074-2081 (1989).

3. D.P. Donovan and A.I. Carswell, ``Principal component analysis applied to multiwavelength lidar aerosol backscatter and extinction measurements", Appl. Opt. 36, 9406-9424 (1997).

4. E.R. Westwater and A. Cohen, ``Application of Backus-Gilbert inversion technique to determination of aerosol distributions from optical scattering measurements", Appl. Opt. 12, 1341-1348 (1973).

5. D.J. Carroll,  ``Chemical laser modeling with genetic algorithms", Amer. J. Aeron. Astron., 34, 338-346 (1996).

6. M.D. King, D.M. Byrne, B.M. Herman and J.A. Reagan, ``Aerosol size distributions obtained by inversion of spectral optical depth measurements", J. Atmos. Sci., 35, 2153-2167 (1978).

7. C.E. Junge,  ``The size distribution and aging of natural aerosols as determined from electrical and optical data on the atmosphere", J. Meteorol., 12, 13-25 (1961).

8. O. Dubovik and M.D. King, ``A flexible algorithm for retrieval of aerosol optical properties from sun and sky radiance measurements", J. Geophys. Res. 105,20, 673-20,696 (2000).

9. M. Kalos and P. Whitlock, {it Monte Carlo Methods, Wiley-Interscience, New York (1986).

10. M. Sambridge, ``Geophysical inversion with a neighbourhood algorithm --I. Searching a parameter space", Geophys. J. Int., 138, 479-494 (1999).

11. Besag, J. and Green, P.J. and Higdon, D. and Mengersen, K., ``Bayesian Computation and Stochastic Systems", Stat. Sci., 10, 3-66 (1995).

12. Gentile, J.E., ``Random Number Generation and Monte Carlo Methods", Springer Verlag, (1998).

13. D.E. Goldberg, in Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley, Reading MA (1989).

14. J.H. Holland, in Adaptation in Artificial and Natural Systems, Univ. Michigan Press, Ann Arbor MI, (1975).

15. D.E. Goldberg, and M.P. Samtani, in ``Engineering optimization via genetic algorithm", in Proceedings of the Ninth Conference on Electronic Computation ASCE, Birmingham, Al, 471-482 (1986).

16. S. Kirkpatrick, Jr C. D. Gelatt, and M. P. Vecchi, ``Optimization by simulated annealing", Science, 220, 37-119 (1983). 

17. A. Corana, M. Marchesi, C. Martini, S. Ridella, ``Minimizing multimodel functions of continuous variables with the simulated annealing algorithm", ACM Trans. Math. Soft., 13, 262-280, (1987).

18. K.T. Whitby, ``The physical characteristics of sulfur aerosols", Atmos. Environ., 12, 135-159, (1978).

19. P.K. Quinn, D.S. Covert, T.S. Bates, V.N. Kapustin, D.C. Ramsey-Bell and L.M. Melnnes,  ``Dimethylsulfide/cloud condensation nuclei/climate system: Relevant size-resolved measurements of the chemical and physical properties of atmospheric aerosol particles", J. Geophys. Res., 98, 10411-10427, (1993).

20. J.N. Porter and A.D. Clarke, ``Aerosol size distribution models based on in situ measurements", J. Geophys. Res., 102, 6035-6045, (1997).

21.  C.F. Bohren and D.R. Huffman, in Absorption and Scattering of Light by Small Particles, Wiley and sons, (1983).

22. D. Deirmendjian, in Electromagnetic Scattering on Spherical Polydispersions, Elsevier, New York, (1969).

23. Post, M.J., ``A graphical technique for retrieving size distributions parameters from multiple measurements:Visualization and error analysis", J. Atmos. Ocean. Tech., 13, 865-873, (1996).

24. B.R. Lienert, E. Berg and L.N. Frazer, ``HYPOCENTER: An earthquake location method using centered, scaled and adaptively damped least squares", Bull. Seis. Am. 76, 771-783, (1986).

25. M.J. Post, ``Limitations of cloud droplet size distribution by Backus-Gilbert inversion of optical scattering data", J. Opt. Soc. Am., 66, 483-486, (1976).

26. J.N. Porter, M. Miller, C. Pietras and C. Motell, ``Ship-based sun photometer measurements using Microtops sun photometers", J. Atmos. Oceanic Technol,. (in press, 2001).

27. G.E. Shaw, ``Sun Photometry", Bull. Amer. Meteor. Soc., 64, 4-11, (1983).

28. G. Backus and F. Gilbert, ``The resolving power of gross earth data'',  Geophys. J. R. Astr. Soc., 266, 169-205, (1969).

29. D.D. Jackson, ``Interpretation of inaccurate, insufficient and inconsistent data", Geophys. J. Roy. Astr. Soc., 28, 97-109, (1972).

30. H.G. Jorge and J.A. Ogren, ``Sensitivity of retrieved aerosol properties to assumptions in the inversion of spectral optical depths", J. Atm. Sci., 53, 3669-3683, ( 1996).