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
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 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.
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.
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).