Iterative regularization method for lidar remote sensing

Iterative regularization method for lidar remote sensing

Computer Physics Communications 174 (2006) 607–615 www.elsevier.com/locate/cpc Iterative regularization method for lidar remote sensing Christine Böc...

293KB Sizes 0 Downloads 113 Views

Computer Physics Communications 174 (2006) 607–615 www.elsevier.com/locate/cpc

Iterative regularization method for lidar remote sensing Christine Böckmann ∗ , Andreas Kirsche University of Potsdam, Institute of Mathematics, Am Neuen Palais 10, 14469 Potsdam, Germany Received 6 December 2004; received in revised form 29 June 2005; accepted 9 December 2005 Available online 15 February 2006

Abstract In this paper we present an inversion algorithm for ill-posed problems arising in atmospheric remote sensing. The proposed method is an iterative Runge–Kutta type regularization method. Those methods are better well known for solving differential equations. We adapted them for solving inverse ill-posed problems. The numerical performances of the algorithm are studied by means of simulations concerning the retrieval of aerosol particle size distributions from lidar observations.  2006 Elsevier B.V. All rights reserved. PACS: 02.60.Nm; 07.05.Kf; 42.68.Wt Keywords: Remote sensing; Laser-radar; Aerosol particle distribution; Inverse ill-posed problem; Iterative regularization

1. Introduction Inverse ill-posed problems are frequently encountered in practice in connection with engineering and physical problems. These range in the latter case from geophysical imaging to atmospheric remote sensing in meteorology and atmospheric physics. Such problems occur, e.g., in the retrieval of molecular concentrations from limb sounding observation [1] and in atmospheric temperature retrievals from a far infrared spectrum observed by an airborne uplooking heterodyne instrument [2] as well as in the ozone profile retrieval from backscattered ultraviolet radiances [3]. In this study we investigate in aerosol science, i.e. in the retrieval of aerosol size distributions from optical data measured by a lidar system (light detection and ranging, i.e. an optical radar system), see [4–6]. Tropospheric aerosols are a potentially important climate forcing factor. Due to their ability to scatter and absorb solar radiation (direct effect) and the modification of the optical properties and lifetime of clouds (indirect effect), they strongly influence the global radiation balance [7]. Despite the great attention atmospheric aerosols have found in the past, quantification of their effects on the climate is still associ* Corresponding author.

E-mail address: [email protected] (C. Böckmann). 0010-4655/$ – see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2005.12.019

ated with large uncertainties [8]. On a global scale, tropospheric aerosol particles seem to lead a negative radiative forcing [9]. Polluted air masses are generally assumed to lead to cooling due to high amounts of sulfates. However, as there are many other absorbing components as, for example, soot or iron oxides, the cooling caused by anthropogenic sulfates may be balanced [10]. In this regard, there is a need to determine microphysical aerosol particle properties by inversion using optical lidar data. The inversion problem in a mathematical sense is ill-posed. Its solution requires the application of appropriate mathematical regularization techniques. The performances of some regularization algorithms for lidar remote sensing were discussed by [11–14]. However, in some applications those regularization tools may lead to unrealistic solutions, therefore, we propose an iterative Runge–Kutta (RK) regularization method with discrepancy principle as parameter choice rule which is an improvement of [12] or [14], respectively. In Section 2, we describe the iterative RK type methods and show some simple test examples for illustrating the methods. In Section 3 we present a particular, very effective two-stage implicit iterative RK method by applying it to the above mentioned retrieving of particle size distribution. Finally, in Section 4 a conclusion is given. For the convenience of the general readers most of the mathematical details only given in Appendices A and B.

608

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

2. Runge–Kutta type iterative regularization methods 2.1. Preliminaries The mathematical ill-posed model which relates the optical and the physical particle parameters, consists of a Fredholm system of two integral equations of the first kind for the backscatter and extinction coefficients r1 n (λ, r; m, s)n(r) dr, Γ (λ) = Kπ/ext (1) r0

where r denotes the particle radius, m is the complex refractive index, assumed to be known and wavelength-independent, s is the shape of the particles, r0 and r1 represent suitable lower and upper limits, respectively, of realistic particle radii, λ is the measurement wavelength, n is the particle number distribution, n is the extinction number kernel. Kπn is the backscatter and Kext The kernel functions reflect shape, size, and material composition of the particles. Here homogeneous spherical particles, i.e. Mie scattering theory [15], are used under the assumption that small particles in a first approximation behave like spherical scatterers [16]. Γ stands for the backscatter coefficient β and/or the extinction coefficient α, respectively, depending on the measurement data. Before coming back to this mathematical model in Section 3 we have to consider some numerical regularization tools for solving such ill-posed problems. We consider the linear operator equation of the form Kx = y,

(2)

where K : X → Y is a linear, compact operator from a Hilbert space X into a Hilbert space Y . The general theory of compact operators evolved from the theory of integral operators of the form s1 y(t) = k(t, s)x(s) ds, t ∈ [t0 , t1 ], (3) s0

which is also the case in lidar remote sensing, see Eq. (1). The variables t, s, s0 , s1 , and functions x, y, k stand here for n , respectively. λ, r, r0 , r1 , and n, Γ, Kπ/ext We assume a square-integrable kernel function k(·, ·) over [t0 , t1 ] × [s0 , s1 ]. The operator (3) is a compact operator from L2 [s0 , s1 ] into L2 [t0 , t1 ]. Notice, in case that K has a nontrivial nullspace x does not need to be uniquely defined by Eq. (2). Therefore, one typically searches the unique solution x† = K †y

(4)

of minimal norm. This solution is equivalent to the solution of minimal norm of the normal equations K ∗ Kx = K ∗ y,

(5)

where K ∗ is the adjoint operator of K and K † is called the Moore–Penrose generalized inverse operator. In particular, we regard the equation Kx = y δ since measurements are made to a finite accuracy and in practice only

a noisy data vector y δ is available. We assume a deterministic error δ with y − y δ   δ. While K is a bounded operator between the Hilbert spaces X and Y , the linear operator K † is often unbounded, namely, if the space X is not compact which is in general the case. This means even small measurement errors δ can cause large uncertainties in the solution function x we are looking for. Therefore, one approximates the inverse operator K † by a linear bounded regularization operator Rγ with a regularization parameter γ > 0. This parameter is used in order to balance approximation error and propagated data error. Roughly spoken, if the parameter γ decreases the approximation error also decreases but the propagated data error increases and vice versa. For this reason, it is important to find the optimal regularization parameter γopt where the sum of both errors is smallest. For selecting an estimation for the optimal parameter we use an a posteriori parameter choice rule in the frame of this paper, namely, the discrepancy principle, e.g., see [17–19]. 2.2. Derivation of the methods Now we explain briefly the derivation of the iterative regularization methods of RK type. For this reason, we examine the initial value problem x  (t) + K ∗ Kx(t) = K ∗ y δ ,

x(0) = 0, t ∈ R.

(6)

We notice Eq. (6) can be seen as a “continuous” analogue of an already well-known iterative regularization method, e.g., if one solves (6) by the explicit Euler method, see Appendix A, with steplength ω, one gets the well-known Landweber iterative regularization method where the steplength ω takes over the role of the relaxation parameter, see [17]. We generalize the idea for the first time to our knowledge in solving Eq. (6) by any explicit or implicit RK method. This technique results in a wide range of new RK type iterative regularization methods for ill-posed problems, see [20–22]. One can derive this generalization for infinite dimensional problems but for the moment for our convenience we discretize the problem to a finite-dimensional one, which is necessary in any case during computer processing. We regard K as an injective matrix as well as x and y δ as vectors in the frame of the next two paragraphs. Eq. (6) is then a system of differential equations and the solution, see [23], is given by    x(t) = (K T K)−1 I − exp t −K T K K T y δ . (7) The adjoint operator is here the transpose matrix K T , I is the identity matrix and exp(t (−K T K)) is the matrix exponential. If t tends to ∞ the solution of the differential equation tends to the solution (4) with y = y δ and K † = (K T K)−1 K T , since K is assumed to be injective. The “regularization” is in integrating this initial value problem (6) not as far as possible (theoretically to +∞), but only up to an abscissa t = 1/γ , where γ serves as the regularization parameter. As γ becomes smaller, i.e. the abscissa up to which one integrates becomes larger, the approximation accuracy improves but in general the propagated data error worsens. In using any explicit RK method with s = p  4 to solve Eq. (6), see Appendix A, Eqs. (A.3) and (A.4), one gets the

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

iteration formula  s   1 i s,ex s,ex T −ωK K xk xk+1 = i! i=0  s−1    1  i +ω −ωK T K K Tyδ , (i + 1)!

(8)

i=0

with (k = 0, 1, 2, . . .) and x0s = 0; s stands for the stage and p is convergence order of the used explicit (ex) RK method, see Appendix A. In using implicit (im) RK methods we get, e.g., for the stages s = 1 or s = 2 with bii  0 (i = 1, 2)  −1 1,im xk+1 = I + ωb11 K T K    × I + ω(b11 − 1)K T K xk1,im + ωK T y δ (9) or 2,im xk+1 = xk2,im + ω(c1 k1 + c2 k2 )

0

for t < s, for s  t,

s ∈ [0, 1], t ∈ [0, 1],

and the right-hand side 1 (−4t 3 + 3t) for t < 1/2, y(t) = 24 1 3 2 24 (4t − 12t + 9t − 1) for t  1/2

with i = 1, 2 and i + 1, i = 1, j= i − 1, i = 2, respectively. We remark that there is no difference in using different explicit RK methods of the same stage as s = p  4. Obviously, the first-stage explicit iteration method, as already mentioned above, is the well-known Landweber method with the relaxation parameter 0 < ω < 2/K2 . In solving the differential equation (6) the relaxation parameter works as limited steplength. Any linear compact operator K between Hilbert spaces and its K † have singular value expansions as follows ∞  1 δ y , wi Y ui , σi i=1 i=1 (11) see Appendix B. Since the singular values σi tends to zero data errors δ in y of fixed size even very small ones can be amplified arbitrarily much, namely by the factors 1/σi , which increase without bound. Many regularization operators Rγ can be represented as filter methods

and

k(t, s)x(s) ds

with the kernel t (1 − s) k(t, s) = s(1 − t)

and

σi x, ui X wi

We consider first two mathematical examples to illustrate the iterative RK regularization techniques itself. It will be obviously that explicit and implicit RK methods of different stages differ only very slightly concerning the retrieval accuracy.

y(t) =

−1 T K K Gi = I + ωbii K K − ω2 bij bj i K T KBjj

∞ 

2.3. Numerical test examples

1

T

Kx =

with a suitable η > 1, occurs.

Example 1. The first kind Fredholm integral equation

 −1 T  2,i G1 · k1 = −K T K I − ωb12 B22 K K xk  −1  T δ T + I − ωb12 K KB22 K y ,  −1 T  2,i K K xk G2 · k2 = −K T K I − ωb21 B11  −1  T δ K y , + I − ωb21 K T KB11 Bii = I + ωbii K K

the number of iteration steps N , i.e. N plays the role of the regularization parameter by γ = 1/N(N ∈ N). The filters are able to damp the amplification of data noise δ during the solution process. We select the optimal number of iteration steps Nopt in using an a posteriori parameter choice rule called discrepancy principle, see [17]. The iteration is terminated by N = N (δ, y δ ) when for the first time δ  ηδ, y − Kx δ (13) δ N (δ,y )

(10)

with

T

609

K †yδ =

has the analytical solution s, for s < 1/2, x(s) = 1 − s, for s  1/2.

(14)

In a first step one has to discretized the integral equation for computer processing. We choose B-splines of second order as basis functions {ϕ1 , . . . , ϕn } and {ψ1 , . . . , ψ m } (n = m = 5) in using Galerkin discretization with xn = ni=1 αi ϕi and Mm,n (α1 , . . . , αn )T = gm where {M}ij := (ψi , Kϕj ) (i = 1, . . . , m; j = 1, . . . , n) and gm := ((y, ψ1 ), . . . , (y, ψm ))T assuming M is regular. In a second step the resulting equation system is solved by different iterative RK regularization methods, see Fig. 1. Whereas the first example is a polynomial one with a kernel function corresponding to boundary value problems of ordinary differential equations the second example is an exponential one.

(12)

Example 2. A second first kind Fredholm integral equation 1 y(t) = 0 k(t, s)x(s) ds with the kernel k(t, s) = exp(ts), s ∈ [0, 1], t ∈ [0, 1], and the right-hand side y(t) = (exp(t +1)−1)/ (t + 1) which has the exact solution x(s) = exp(s) is tested. Here we use first order B-splines (m = n = 50) and again Galerkin discretization. The results are shown in Fig. 2.

with filter functions Fi = F (γ , σi ); the developed RK methods, too, see Appendix B. We can interpret the RK methods as discrete filters whose error-sensitivity bandwidth is controlled by

Additionally, in both examples the integrals were computed by using the adaptive quadrature rule of Simpson from the MATLAB 6.5 library on a PC platform.

xγδ = Rγ y δ =

∞  i=1

F (γ , σi )

y δ , wi Y ui σi

610

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

Fig. 1. Numerical inversion results of noiseless and noisy data with different noise levels for Example 1 using Galerkin discretization and different Runge–Kutta iteration methods: (a) implicit (s = 1; Backward–Euler), (b) implicit (s = 2; Radau), (c) explicit (s = 1; Landweber), (d) explicit (s = 4; historical RK).

Fig. 2. Numerical inversion results of noiseless and noisy data with different noise levels for Example 2 using Galerkin discretization and different Runge–Kutta iteration methods: (a) implicit (s = 1; Backward–Euler), (b) implicit (s = 2; Radau), (c) explicit (s = 1; Landweber), (d) explicit (s = 4; historical RK).

Figs. 1 and 2 show the computed solution by the implicit Euler method (a), by the implicit Radau method (b), by the explicit Euler method (c) and by the explicit fourth-stage RK method (d), see Appendix A for the particular methods, with noiseless data as well as with noisy data of different noise levels δ (δ = 0.01; 0.05; 0.08; 0.1) for the first example and δ (δ = 0.01; 0.05; 0.2; 1.0) for the second one, respectively. The noisy right-hand side was determined by y δ (t) = y(t) + δ ∗ cos(100t), the parameter η was set to 1.02 and for the steplength was cho-

sen ω = 9/10 ∗ 1/K22 where, as usual, the subscript 2 denotes the spectral norm. Figs. 1 and 2 show quite good results also for noisy data and demonstrate that the developed iterative RK regularization methods work well. Concerning the reconstruction accuracy in using different RK methods one can see in the noiseless data case that the differences are negligible. In the next section one will discover that the Radau method is very favorable concerning the computer effectiveness.

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

3. Application to the retrieval of aerosol size distributions The mathematical model in Eq. (1) can be reformulated into a more specific form for homogeneous spheres and in using the volume distribution v(r) instead of the number distribution n(r) for stability reasons as r1 Γ (λ) =

v Kπ/ext (λ, r; m)v(r) dr r0

r1 =

3 Qπ/ext (λ, r; m)v(r) dr, 4r

(15)

r0

with 4πr 3 (16) n(r). 3 The backscatter and extinction efficiencies Qπ/ext (λ, r; m) for v homogeneous spheres can be computed from [16]. Kπ/ext (λ, r; m) are the volume kernel functions for backscatter and extinction, respectively, which are comparable with the kernel function k(t, s) of the previous section and v(r) stands for x(s).

v(r) =

611

In using a three-wavelength lidar with two Raman channels [4], i.e. we have data for the backscatter coefficient at three wavelengths and for the extinction coefficient at two wavelengths (case 3 + 2) namely β Aer (355 nm), β Aer (532 nm), β Aer (1064 nm) as well as α Aer (355 nm), α Aer (532 nm). Monomodal logarithmic-normal distributions are used to describe the particle number distributions   (ln r − ln rmed )2 1 Nt exp −0.5 n(r) = (17) , √ r 2π ln σ ln2 σ where Nt , σ and rmed denote the total number concentration, the geometric standard deviation and the median radius, respectively. We consider four examples (j = 1, . . . , 4), namely, rmed,j = 0.05, 0.1, 0.3, 0.1 µm, σj = 1.8, 1.8, 2.0, 1.6 for the distribution in Eq. (17) and the refractive index was mj = 1.7 + 0.01i, 1.7 + 0.05i, 1.4 + 0.01i, 1.5 + 0.01i. Discretization by collocation with {M}ij := (Kϕj )(ti ) (i = 1, . . . , m; j = 1, . . . , n) and gm := (y(t1 ), . . . , y(tm ))T at the collocation points ti (i = 1, . . . , m) needs less computing time as the Galerkin discretization. Therefore, discretization by collocation with 3 + 2 data points (coefficients above), B-spline basis functions of order 4 (n = 20) and exact right-hand side as well as a perturbed

(a)

(b)

(c)

(d)

Fig. 3. Retrieval results of the volume distribution of Example 1 without noise (a) with noise (b) and of Example 2 without noise (c) with noise (d) for the measurement case 3 + 2. In particular 1% normally distributed noise results in 2.22% relative deterministic error (b) and 2.17% (d).

612

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

(a)

(b)

(c)

(d)

Fig. 4. Retrieval results of the volume distribution of Example 3 without noise (a) with noise (b) and of Example 4 without noise (c) with noise (d) for the measurement case 3 + 2. In particular 1% normally distributed noise results in 2.07% relative deterministic error (b) and (d).

right-hand side (β(λ1 ), β(λ2 ), β(λ3 ), α(λ1 ), α(λ2 ))T by random numbers of normally distributed noise of 1% are used. Additionally, in applications the measurement data points are discrete ones, i.e. the function y(t) is not known analytically as in Examples 1 and 2 of the previous section which is necessary for Galerkin discretization. The measurement data, i.e. the wavelengths, are mostly taken as collocation points (m = 5). B-spline basis functions of order 4 are locally cubic polynomials. In retrieving a particle distribution we are in general looking for such a shaped function which can be approximated by locally cubic polynomials very well in contrast to Example 1 of Section 2. Therefore, in the frame of Example 1 we used locally linear polynomials to approximate the solution. The Bspline order choice, roughly spoken, depends on the shape of the expected solution. Besides in using a larger number of basis functions, see Example 2 of the previous section, the choice of the B-spline order does not seem to be so important. Moreover, in general a distribution grows steeper on the lefthand side as on the right-hand side. Therefore, we deal with a non-equidistant node grid for the B-splines by using the roots of the Chebyshev polynomials on the left-hand side in contrast to Examples 1 and 2. The expected shape serves as useful a-priori information in the retrieving process.

The results are shown in Figs. 3(a)–(d), 4(a)–(d). We compare two RK methods the Landweber iteration and the Radau iteration. As already remarked in Section 2 both methods show more or less the same very good retrieval quality if one has in mind that only 5 measurements are used. But the Landweber iteration needs a huge number of iteration steps between 37 and 14,038 steps, see the legends of Figs. 3(a)–(d), 4(a)–(d), whereas the Radau iteration is a very effective one it needs only between 3 and 7 iteration steps, see Appendix B or for more details [21,22]. Finally, we retrieved the effective radius (µm), reff = 3

vt at

(18)

with

 v(r) at = 3 dr (µm2 cm−3 ) and r  vt = v(r) dr (µm2 cm−3 )

(19)

for different noise levels and a wider number of examples. Additionally three examples (k = 5, 6, 7) were examined, namely, rmed,k = 0.1, 0.1, 0.3 µm, σk = 1.4, 1.4, 1.6 for the distribu-

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

(a)

613

(b)

Fig. 5. Comparison of the retrieval results for the effective radius concerning the measurement cases 6 + 2 (a) and 3 + 2 (b).

tion in Eq. (17) and the refractive index was mk = 1.4 + 0.5i, 1.55 + 0.5i, 1.4 + 0.005i. For the inversion the coefficients were imposed with 0 (exact case), 0.1, 0.5, 1.0, 5.0 and 10% of normally distributed noise. 100 runs per noise level and example were investigated in order to increase the statistical significance. The results are shown in Fig. 5(b). In using Landweber or Radau method there is again no remarkable difference concerning the retrieval errors. But the average of iteration steps is again higher for the Landweber method. Fig. 5(a) compares the same circumstances but for the case 6 + 2. Case 6 + 2 uses 3 additional backscatter coefficients at the wavelengths 400, 710 and 800 nm which can be provided, e.g., by [25]. Fig. 5(a), (b) shows that for larger noise levels 5 and 10%, which is a more realistic noise range for optical lidar data, the three additional used backscatter coefficients have no influence on the retrieved quality. This result was already noticed in [14] in using a different and more time consuming regularization method. Only for very small noise levels between 0 and 1% the three additional used backscatter coefficients are very important for the retrieval quality. The error to noise examinations are not directly comparable with those which were made in [14]. Here they are slightly larger because of two reasons. First, in the frame of this paper the noise was added more realistically to each lidar channel whereas in previous examinations a perturbation with respect to the norm of the data vector was added. Second, in using a priori information about the noise level an automatic a posteriori regularization parameter choice rule was adapted.

ing the necessary computer run time. The Radau method needs only a very small number of iteration steps until the discrepancy principle terminates the iteration since the steplength ω, roughly spoken, can be chosen theoretically unlimited. This paper is a first step towards the development of a fast algorithm without loosing accuracy in the retrieval process to determine microphysical aerosol particle properties from optical lidar measurements. In a second part we will investigate in the more extensive case of an a priori unknown complex refractive index and in adding a strong nonnegative constraint. Acknowledgements The authors are grateful to an anonymous referee for his or her very valuable and extensive suggestions to rewrite the paper for a more general reader without particular mathematical knowledge. This work was supported by the PEP project from the HGF “Impuls und Vernetzungsfond” under the grant VHVI-100. Appendix A. Brief summary of Runge–Kutta methods We give here a very short description of the well-known Runge–Kutta (RK) methods, see [24], as solvers for initial value problems only as necessary in the frame of this paper. We regard an ordinary differential equation system of first order   x  (t) = f t, x(t) with x(t0 ) = x0 . (A.1) One can get a discrete approximative solution for the function x in Eq. (A.1) by (k = 0, 1, 2, . . .)

4. Conclusion xk+1 = xk + ω Iteratively regularizing RK methods for atmospheric retrievals are presented. Our analysis is focused on the improvement of the computer run time. The iterative Radau regularization method, an implicit two stage (s = 2) RK method allows us to derive microphysical particle properties from multiwavelength lidar measurements with an automatic a posteriori regularization parameter choice rule, the discrepancy principle, in a very effective manner concern-



s 

c i ki

with

i=1

ki = f tk + ai ω, xk + ω

s 

 bil kl ,

(A.2)

l=1

where s is the stage of the RK method and ω is the steplength with tk+1 = tk + ω. Usually the parameters for different methods and different stages are collected in so-called Butcher-tableaus

614

a1 .. .

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

b11 .. .

... .. .

b1s .. .

as bs1 . . . bss c 1 . . . cs

xkδ,s,ex

or shorter as a triple (a, B, c).

Those parameters have to fulfill the usual order conditions to achieve a particular convergence order p. If the matrix B is a strongly lower tridiagonal one the method is a simple explicit RK method, otherwise it is an implicit one. In the frame of this paper we use two first-stage (s = 1) methods, namely, the explicit Euler method (0, 0, 1) (Landweber iteration) and the implicit Backward–Euler method (1, 1, 1). Furthermore, we regard “the” historical explicit fourth-stage (s = 4) RK method and the implicit two-stage (s = 2) Radau method 0 1/2 1/2 1

0 1/2 0 0 1/6

0 0 1/2 0 2/6

0 0 0 1 2/6

0 0 0 0 1/6

and

s = xks + ω xk+1

s  i=1

j =1



ci K T y δ − K T K xks + ω

s  j =1

(B.2)

i=0

(k = 1, 2, . . .) as well as in using the singular system by  s j k−1  ∞    1 δ,s,ex 2 i xk = −ωσl i! l=1 j =0 i=0  s−1      1  2 i ×ω (B.3) −ωσl σl y δ , wl Y ul . (i + 1)! i=0

For k = N iteration steps we have

1/3 5/12 −1/12 1 3/4 1/4 3/4 1/4

δ,s,ex s δ = RN y xN

=

respectively. In solving the initial value problem in Eq. (6), i.e. for the particular function f (t, x(t)) = K T y δ − K T Kx(t), where K is regarded again as injective matrix, by RK methods the terms ki (i = 1, . . . , s) and xk+1 (k = 1, 2, . . .) of Eq. (A.2) get the following form   s  T δ T s bij kj , ki = K y − K K xk + ω (A.3) 

 s j k−1   1 ∗ i = (−ωK K) i! j =0 i=0  s−1   1 ×ω (−ωK ∗ K)i K ∗ y δ (i + 1)!

 bij kj

. (A.4)

Appendix B. Runge–Kutta iterations as filter regularization methods For any linear compact operator K : X → Y , between Hilbert spaces X and Y a singular system (σi ; ui , wi )i∈N exists. The {σi2 }i∈N are the eigenvalues of the self-adjoint operators K ∗ K and KK ∗ , written down in decreasing order with multiplicity. The singular functions {ui }i∈N are a corresponding complete orthonormal system of eigenfunctions of K ∗ K and the singular functions {wi }i∈N are defined via wi := Kui /Kui X and form also a complete orthonormal system of eigenfunctions of

KK ∗ . It holds Kx = ∞ i=1 σi x, ui X wi with the inner product ·, ·X in X. In general, there are infinitely many singular values, they accumulate only at zero, i.e. limi→∞ σi = 0. As already remarked in Section 2 in Eq. (12) many methods are filter methods. A continuous filter has regularization properties if it fulfills the two conditions    (a) F (γ , σi )  CF with a constant CF > 0, (B.1) (b) limγ →0 F (γ , σi ) = 1. Those filters are able to damp the amplification of data noise δ during the solution process. We show in [22] that those particular Runge–Kutta methods are filter methods in the sense of (12) and (B.1). By induction we get that the iterates xkδ,s,ex may be expressed non-recursively through

∞  1−(

s

1 2 i N i=0 i! (−ωσl ) )

σl

l=1

(y δ , wl )Y ul

(B.4)

and the filter for any explicit RK iteration with s = p  4 arises as  s N  1  i F ex (N, σl ) = 1 − (B.5) −ωσl2 i! i=0

with N = 1/γ , N ∈ N. With the discrepancy principle (13) one has to choose an optimal finite number Nopt of iteration steps. Furthermore, as a general representation for explicit and implicit RK iterations we found for the filters F (N, σl ) = 1 − r(−ωσl2 )N with rational functions r(x) = p(x)/q(x) and p, q are polynomials. For all explicit methods we have q(x) ≡ 1. Our RK filters fulfill the conditions (B.1) with CF  2. For certain steplengths ω our explicit and implicit RK iterations are convergent and in using the discrepancy principle they are order optimal regularization methods with infinite qualification. For more details see [22]. Additionally, the number of iteration steps until the discrepancy principle terminates the iteration depends on the used steplength ω. The steplength is restricted by |r(−ωσ 2 )| < 1 for all σ ∈ (0, K]. Thus the steplength depends on the used method and is limited for many methods as for the Landweber iteration by 0 < ω < 2/K2 . Therefore, such methods are very slow. In [21] was shown that the steplength for the Radau method, see Appendix A, is unlimited. Therefore, this method is a very fast one. References [1] A. Doicu, F. Schreier, M. Hess, Iteratively regularized Gauss–Newton method for bound-constraint problems in atmospheric remote sensing, Comput. Phys. Comm. 153 (2003) 59–65. [2] A. Doicu, F. Schreier, M. Hess, Iteratively regularized Gauss–Newton method for atmospheric remote sensing, Comput. Phys. Comm. 148 (2002) 214–226. [3] O. Hasekamp, J. Landgraf, Ozone profile retrieval from backscattered ultraviolet radiances: the inverse problem solved by regularization, J. Geophys. Res. 106 (2001) 8077–8088.

C. Böckmann, A. Kirsche / Computer Physics Communications 174 (2006) 607–615

[4] V. Matthias, J. Bösenberg, V. Freudenthaler, A. Amodeo, I. Balin, D. Balis, A. Chaykovski, G. Chourdakis, A. Comeron, A. Delaval, F. de Tomasi, R. Eixmann, A. Hagard, L. Komguem, S. Kreipl, R. Matthey, V. Rizi, J.A. Rodriguez, U. Wandinger, X. Wang, Aerosol lidar intercomparison in the framework of the EARLINET project. 1. Instruments, Appl. Opt. 43 (2004) 961–976. [5] C. Böckmann, U. Wandinger, A. Ansmann, J. Bösenberg, V. Amiridis, A. Boselli, A. Delaval, F. de Tomasi, M. Frioud, I.V. Grigorov, A. Hagard, M. Horvat, M. Iarlori, L. Komguem, S. Kreipl, G. Larcheveque, V. Matthias, A. Papayannis, G. Pappalardo, F. Rocadenbosch, J.A. Rodriguez, J. Schneider, V. Shsherbakov, M. Wiegner, Aerosol lidar intercomparison in the framework of the EARLINET project. 2. Aerosol backscatter algorithms, Appl. Opt. 43 (2004) 977–989. [6] G. Pappalardo, A. Amodeo, M. Pandolfi, U. Wandinger, A. Ansmann, J. Bösenberg, V. Matthias, V. Amiridis, F. De Tomasi, M. Frioud, M. Iarlori, L. Komguem, A. Papayannis, F. Rocadenbosch, X. Wang, Aerosol lidar intercomparison in the framework of the EARLINET project. 3. Raman lidar algorithm for aerosol extinction, backscatter and lidar ratio, Appl. Opt. 43 (2004) 5370–5385. [7] R.J. Charlson, T.M.L. Wigley, Sulfate aerosol and climatic-change, Sci. Amer. 270 (1994) 48–57. [8] T.L. Anderson, R.J. Charlson, S.E. Schwartz, R. Knutti, O. Boucher, H. Rode, J. Heintzenberg, Climate forcing by aerosol—a hazy picture, Science 300 (2003) 1103–1104. [9] The Intergovernmental Panel on Climate Change (IPCC), IPCC Third Assessment Report—Climate Change 2001: The Scientific Basis, Cambridge University Press, 2001. [10] M.O. Andreae, The dark side of aerosol, Nature 409 (2001) 671–672. [11] D. Müller, U. Wandinger, A. Ansmann, Microphysical particle parameters from extinction and backscatter data by inversion with regularization: theory, Appl. Opt. 38 (1999) 2347–2357. [12] C. Böckmann, Hybrid regularization method for the ill-posed inversion of multi-wavelength lidar data in the retrieval of aerosol size distributions, Appl. Opt. 40 (2001) 1329–1342.

615

[13] I. Veselovskii, A. Kolgotin, V. Griaznov, D. Müller, U. Wandinger, D. Whitemann, Inversion with regularization for the retrieval of tropospheric aerosol parameters from multiwavelength lidar sounding, Appl. Opt. 41 (2002) 3685–3699. [14] C. Böckmann, I. Mironova, D. Müller, L. Schneidenbach, R. Nessler, Microphysical aerosol parameters from multiwavelength lidar, J. Opt. Soc. Amer. A 22 (2005) 518–528. [15] G. Mie, Beiträge zur Optik trüber Medien speziell kolloidaler Metallösungen, Ann. Phys. 25 (1908) 377–445. [16] G.F. Bohren, D.R. Huffman, Absorption and Scattering of Light by Small Particles, John Wiley and Sons, New York, 1983. [17] H.W. Engl, M. Hanke, A. Neubauer, Regularization of Inverse Problems, Mathematics and its Applications, vol. 375, Kluwer Academic Publishers, Dordrecht, 1996. [18] S.F. Gilyazov, N.L. Gol’dman, Regularization of Ill-Posed Problems by Iteration Methods, Kluwer Academic Publishers, Dordrecht, 2000. [19] C.R. Vogel, Computational Methods for Inverse Problems, SIAM, Philadelphia, PA, 2002. [20] C. Böckmann, Runge–Kutta type methods for ill-posed problems and the retrieval of aerosol size distributions, in: Proc. GAMM-Conference 2001, in: Proc. in Appl. Math. and Mech. (PAMM), vol. 1, Wiley–VCH Verlag, Weinheim, 2002, pp. 486–487. [21] A. Kirsche, C. Böckmann, Rational approximations for ill-conditioned equation systems, Appl. Math. Comput. 171 (2005) 385–397. [22] A. Kirsche, C. Böckmann, Pade Iteration methods for regularization, Appl. Math. Comput. 172 (2005) (accepted). [23] H. Heuser, Gewöhnliche Differentialgleichungen, B.G. Teubner, Stuttgart, 1991. [24] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II, Springer-Verlag, Berlin Heidelberg, 1991. [25] D. Althausen, D. Müller, A. Ansmann, U. Wandinger, H. Hube, E. Clauder, S. Zörner, Scanning six-wavelength eleven-channel aerosol lidar, J. Atmos. Ocean. Technol. 17 (2000) 1469–1482.