气溶胶粒径分布.pdf

Hybrid regularization method for the ill posed inversion of multiwavelength lidar data in the retrieval of aerosol size distributions Christine Bo ckmann A specially developed method is proposed to retrieve the particle volume distribution the mean refractive index and other important physical parameters e g the effective radius volume surface area and number concentrations of tropospheric and stratospheric aerosols from optical data by use of multiple wavelengths This algorithm requires neither a priori knowledge of the analytical shape of the distri bution nor an initial guess of the distribution As a result even bimodal and multimodal distributions can be retrieved without any advance knowledge of the number of modes The nonlinear ill posed inversion is achieved by means of a hybrid method combining regularization by discretization variable higher order B spline functions and a truncated singular value decomposition The method can be used tohandledifferentlidardevicesthatworkwithvariousvaluesandnumbersofwavelengths Itisshown to my knowledge for the fi rst time that only one extinction and three backscatter coeffi cients are suffi cient for the solution Moreover measurement errors up to 20 are allowed This result could be achieved by a judicious fusion of different properties of three suitable regularization parameters Fi nally numerical results with an additional unknown refractive index show the possibility of successfully recovering both unknowns simultaneously from the lidar data the aerosol volume distribution and the refractive index 2001 Optical Society of America OCIS codes 010 0010 010 1110 100 0100 100 3190 280 0280 290 0290 1 Introduction Aerosol particle properties can be derived by mea surement of a variety of scattering properties These properties describe the particle s infl uence on the Earth s radiation budget on clouds and on pre cipitation as well as their role in chemical processes of the troposphere and the stratosphere This may include extinction or scattering information at mul tiple wavelengths scattering information at multiple angles or multiple scattering information Here the inversion of particle properties from lidar mea surements of backscatter and extinction at multiple wavelengths is discussed The inversion requires the solution of a Fredholm system of integral equations of the fi rst kind The mathematical model consists of two integral equations for the backscatter and extinction coeffi cients baerand aaer baer l z 5 rmin rmax k p r l m n r z dr 5 rmin rmax pr2Qp r l m n r z dr 1 aaer l z 5 rmin rmax k ext r l m n r z dr 5 rmin rmax pr2Qext r l m n r z dr 2 Here r denotes the particle radius m is the refractive index rminand rmaxrepresent suitable lower and up per limits respectively of realistic radii l is the wave length z is the height n is the aerosol size distribution we are looking for k pis the backscatter and k extis the extinction kernel The kernel functions refl ect the shape size and material composition of the particles WeassumeMieparticles Whenactuallidardataare The author bockmann rz uni potsdam de is with the Institute of Mathematics University of Potsdam PF 601553 D 14415 Pots dam Germany Received 4 April 2000 revised manuscript received 31 August 2000 0003 6935y01y091329 14 15 00y0 2001 Optical Society of America 20 March 2001 y Vol 40 No 9 y APPLIED OPTICS1329 used therefractiveindexisalsoanunknown There fore the problem is a nonlinear ill posed problem i e the correct kernel function is also unknown On the other hand all mathematical algorithms currently in use assume the kernel function to be known correctly with the following formulas1providing the extinction and the backscatter effi ciencies Qpand Qextof Eqs 1 and 2 Qp5 1 k 2r2U n51 2n 1 1 21 n an2 cn U 2 Qext5 2 k 2r2 n51 2n 1 1 Re an1 cn 3 where k is the wave number defi ned by k 5 2pyl and anand cn are the coeffi cients that one gets from the boundary conditions The standard formula1for computing Qpand Qextin the case of homogeneous spheres is an5 mcn mv cn9 v 2 cn v cn9 mv mcn mv zn9 v 2 zn v cn9 mv 4 cn5 cn mv cn9 v 2 mcn v cn9 mv cn mv zn9 v 2 mzn v cn9 mv 5 cn t 5 tjn t xn t 5 tyn t zn t 5 c t 1 x t i where jnand ynare spherical Bessel functions and v 5 k r The lidar setup2of the Institute for Tropospheric Research in Leipzig Germany typically uses six backscatter and two Raman channels see the fi rst setup case in Table 1 This number cannot how ever be achieved with standard currently available lidar systems which work with only three backscat terandoneRamanchannel seethefi fthsetupcasein Table 1 Moreover it is probable that during a mea surement campaign one or two channels may fail This inspired the consideration of the other three setup cases which are shown as cases 2 4 in Table 1 The resulting determination of the aerosol size dis tribution function n r from only a small number of backscatter and extinction measurements by means of a system of Fredholm integral equations of the fi rst kind is the most familiar and common example not only for a compact operator but also for a severely ill posed inverse problem Inthepastfewyearsdifferentinversiontechniques have been proposed for solving similar multiwave length ill posed inversions of more special operators The algorithms3 4are based on quadrature discreti zations that contain no regularization effect this will be considered in Section 2 The fi rst algorithm3 deals with the inversion of optical thickness data by means of the Tikhonov regularization TR at seven equally spaced wavelengths The second algorithm4 deals with scattering intensity data by means of trun cated singular value decomposition TSVD regular ization by use of between fi ve and eight wavelengths at a scattering angle of 176 References 5 and 6 propose a review of different regularization tech niques and compare some of these constraint least squares TSVD Landweber iterations7 by determi nation of the global aerosol distribution for which solar radiation data with known basic aerosol com ponents are used Other interesting works in re lated fi elds describe for example channel selection by use of optical depth data 8a purely graphic tech nique for parameter estimations of known distribu tion shapes 9or the aerosol mass concentration distribution 10Also discussed are analytic inver sions11 12of the special nonabsorbing anomalous dif fraction extinction operator12or the optical thickness data inversion by truncated analytic eigenfunction expansion 11with the Mie theory extinction kernel used as a product kernel Qext r l m 5 Qext kr m An overview and a critical review of different inver sion techniques for in situ devices can be found in Ref 13 References 14 20 contain the most related al gorithms with respect to lidar ill posed inversion The methods of Refs 14 17 work with projection discretization by means of fi rst order B splines and by means of TR16at eight wavelengths or a Bayesian based iterative regularization method 15respectively References 18 20 work with second order B splines by means of TR but still with eight wavelengths Of these only Ref 20 treats the case of an unknown refractive index In this paper a new hybrid regularization method is proposed that uses variable projection dimension and variable B spline order as well as TSVD simul taneously for that severely ill posed inversion Fi nally the method can handle the inversion with only oneextinctionandthreebackscattercoeffi cientsupto 10 noise and for two extinction and six backscatter coeffi cients up to 20 noise In addition the un known refractive index can be captured The paper is organized as follows After the mathematical background is described in Section 2 the specially developed hybrid regularization method is proposed in Section 3 and the projection method as a regular ization tool is discussed Moreover inversion re sults for simulated data are shown in order to fi nd a suitable B spline basis In Section 4 several mea surement situations with simulated noiseless and noisy data for numerous different distributions with known refractive indices are studied as well as sev eral examples with additional unknown refractive in dices In these latter examples we are able to retrieve simultaneously both the distribution and the refractive index Finally conclusions are given in Section 5 Table 1 Different Common Setup Cases of Multispectral Lidar Devices Setup CaseLp nm NLext nm M 1355 400 532 710 800 10646355 5322 2355 400 532 710 800 10646 0 3355 400 532 800 10645355 5322 4355 400 532 710 8005355 5322 5355 532 106435321 1330APPLIED OPTICS y Vol 40 No 9 y 20 March 2001 2 Mathematical Background Equations such as Eqs 1 and 2 i e with compact operators K are always ill posed generally on all three counts existence uniqueness and stability for which stability21means a solution that changes only slightly with a slight change in the problem We consider an equation of the form Kx 5 y where K X 3 Y is a compact linear but not necessarily self adjoint operator from a Hilbert space X into a Hilbert space Y The general theory of compact op erators evolved from the theory of integral operators of the form y l 5 rmin rmax k l r x r dr 6 Indeed if k z z is a square integrable over lmin lmax 3 rmin rmax then it is a well known classic result that K is a compact operator from L2 rmin rmax into L2 lmin lmax In addition if k z z is contin uous as in our case then K is a compact operator from C rmin rmax into C lmin lmax The operators K K X 3 X andKK Y 3 Y arecompactself adjoint linear operators where K is the adjoint operator of K The nonzero eigenvalues of K K or of KK they have the same eigenvalues can be enumerated as l1 l2 If we designate by v1 v2 an associated sequence of orthonormal eigenvectors of K K then vj is complete in the range R K K 5 N K orthogonal complement of the null space N K of K Let mj5 lj then Kvj5 mjujand K uj5 mjvj Moreover KK uj5 mjKvj5 mj2uj5 ljuj and it is not diffi cult to see that the orthonormal eigenvectors uj of KK form a complete orthogonal set for R KK 5 N K The system vj uj mj is called a singular system for K and the numbers mjare called singular valuesofK Inthenondegenerate infi nite dimensional case the values li i 5 1 2 of the singular value expansion SVE of the operator clus ter at zero Any function of the form x 5 j51 y uj see also the remark in Section 5 Second additional information about the noise d i e the measurement errors may be helpful e g deterministic or stochastic information about the noise Third one can deal with subjective a priori information by using a support functional V D X 3 R V x 5 min x D X 9 where x is compatible with the data This is the nondescriptive regularization for which the solution has certain desired features such as smoothness In general regularizations are families of operators Kgwith g R 23 Kg Y 3 X with limg30Kgy 5 K y for all y D K i e the convergence is pointwise on D K The pa rameter g is the regularization parameter In the case of noisy data ydwith iy 2 ydi d we determine a solution xgd5 Kgyd However the total error consists of two parts i e two summands xgd2 x 5 Kg yd2 y 1 Kg2 K y 10 The fi rst part is the data error and the second part is the approximation or regularization error If g 3 0 the approximation error tends to zero whereas the data error tends to infi nity Therefore the total er ror can never be zero We are looking for a trade off between these two terms i e for an optimal regular ization parameter g that minimizes the total error Many regularization methods can be found in Ref 24 First the most popular and well known regulariza tion is the classic TR or the iterative TR 25There are other examples such as truncated and general izedtruncatedsingular valuedecomposition23 26 SVD asymptotic regularization 24iterative meth ods e g linear Landweber iterations 7nonlinear conjugate gradient iterations 27 28Lardy method or Schulz method25 mollifi er methods i e fi rst smooth the function y by mollifi cation29 31and then approx imate the mollifi ed function or maximum entropy methods 32 Second there are also many different discretiza tion possibilities simple classic quadrature meth ods with different weights collocation points and nodes degenerate kernel approximations expan sions by eigenfunctions or by orthonormal systems approximations by Taylor or interpolation projec tion methods Galerkin moment collocation or least squares methods33 35 and multilevel algo rithms 36 The fi rst two discretization methods ap proximate the underlying equation whereas the projection methods approximate the solution On 20 March 2001 y Vol 40 No 9 y APPLIED OPTICS1331 the one hand it is now possible to combine any reg ularization method with any discretization method when one recognizes that each numerical algorithm has its advantages and its drawbacks On the other hand it is well known that the straightforward naive quadrature approach to integral equations of the fi rst kind is usually unsatisfactory and leads to disastrous results if the data contain any errorsor uncertainties The approximations become even worse for fi ner quadrature discretization schemes Fortunately in a projection method discretization and regulariza tion work hand in hand to produce the linear system this happens because pure projection into fi nite dimensional spaces itself acts as regularization in which the regularization parameter is the dimension n i e g 5 1yn see Section 3 This is in sharp contrast to the fi rst method in which the discretiza tion by quadrature takes place fi rst without any reg ularization and then an additional regularization is performed after the discretization is complete 37 Therefore the application of the quadrature methods is in principle confi ned to equations of the second kind and these are generally well posed Because my experience shows that such a standard method works well for only small measurement er rors I decided to develop a hybrid regularization method to amalgamate different advantages of vari ous methods 3 Hybrid Regularization Method Based on knowledge of the properties of different stan dard discretization and regularization techniques a specially developed hybrid regularization technique is proposed that can handle different lidar systems that work with various values and numbers of wavelengths and up to 20 noisy data The regularization process consists of a discretization part by means of projection and a regularization part the fi rst part itself has reg ularization properties To ensure that both regular ization stages work hand in hand we use three regularization parameters simultaneously The fi rst is the projection dimension the second is the order of the B spline basis used and the last which is impor tant during the solution of the resulting system of linear equations is the level of the TSVD A Regularization by Discretization by Means of Projection Methods Let X be a normed space and U X be a nontrivial subspace A bounded linear operator P X 3 U with the property Pt 5 t for all t U is called a projection operator from X onto U Two important examples for projection operators are given by the orthogonal pro jection e g Galerkin methods and by the interpola tion projection e g collocation methods We study thelastone LetXn XandYm Ybetwosubspaces with dim Xn5 n dim Ym5 m and with basis Xn5 span f1 fn Ym5 span r1 rm 11 Different possibilities of taking out a projection38 39 exist The dual projection method is preferred and the approximation xn X with xn5 b1f11 1 bnfn 12 is the solution of the projected equation PmKxn5 Pmy 13 Here Pm Y 3 Ymis an appropriate projection op erator let Km 5 PmK The projected equation yields a system of linear equations for all n unknown coordinates bi i 5 1 n The error of the pro jection method depends on how well the exact solu tion can be approximated by elements of the subspace Xn Increasing n and m will make the error smaller in relation to the discretization However because of ill conditioning the computation will be more con taminated by errors of the given data y Note that in actual numerical computations such errors will occur automatically because of round off effects on the data On the other hand if n and m are small then the approximation is robust against errors in y but will be inaccurate because of a large discretization error This problem is the same as in Eq 10 which calls for a compromise in the choice of discretization parameters n and m In general all discretizations have shown that they lead to matrices with a high condition number 40which grows with the increasing dimension of the matrix Besides the degree of ill conditioning of a matrix is measured by the condition number which is the ratio between the largest and the smallest singular value Let K y X be the solution of Kx 5 y and m 5 n By xn Xn we denote the unique solution of the equation PnKxn5 Pny We can represent the solutions xn Xnin the form xn5 Rny where Rn Y 3 Xn X is defi ned by Rn 5 PnKuXn Pn Y 3 Xn X The projection method is called convergent if the approximate solutions xn Xnconverge to the exact solution K y X for every y R K i e if Rny 5 PnKuXn Pny 3 K y n 3 14 for every y D K We observe that this defi nition of convergence coincides with the defi nition of the previous section i e a regularization strategy for the equation Kx 5 y with regularization parameter g d 5 1yn Therefore the projection method converges if and only if Rnis a regularization strategy for the equation Kx 5 y Hence pure discretization turns out to be a regularization method For convergence proofs of the special case m 5 n see Refs 24 and 33 In more detail we will now see that there is a hidden regularization parameter namely the smallest sin gular value nnof the operator Kn 5 PnK where Pnis the orthogonal projector onto Ynwith iPni 5 1 For a stability analysis we assume that noisy data ydare 1332APPLIED OPTICS y Vol 40 No 9 y 20 March 2001 available with iPn y 2 yd i d Because iKn i 5 1ynnwe obtain i