Adaptive filter solution for processing lidar returns: optical parameter estimation

Adaptive filter solution for processing lidar returns: optical parameter estimation Francesc Rocadenbosch, Gregori Va´zquez, and Adolfo Comero´n Joint estimation of extinction and backscatter simulated profiles from elastic-backscatter lidar return signals is tackled by means of an extended Kalman filter ~EKF!. First, we introduced the issue from a theoretical point of view by using both an EKF formulation and an appropriate atmospheric stochastic model; second, it is tested through extensive simulation and under simplified conditions; and, finally, a first real application is discussed. An atmospheric model including both temporal and spatial correlation features is introduced to describe approximate fluctuation statistics in the sought-after atmospheric optical parameters and hence to include a priori information in the algorithm. Provided that reasonable models are given for the filter, inversion errors are shown to depend strongly on the atmospheric condition ~i.e., the visibility!and the signal-to-noise ratio along the exploration path in spite of modeling errors in the assumed statistical properties of the atmospheric optical parameters. This is of advantage in the performance of the Kalman filter because they are often the point of most concern in identification problems. In light of the adaptive behavior of the filter and the inversion results, the EKF approach promises a successful alternative to present-day nonmemory algorithms based on exponential-curve fitting or differential equation formulations such as Klett’s method. © 1998 Optical Society of America OCIS codes: 010.0010, 010.1290, 010.3640. 1. Introduction The superior qualities of laser radars, or lidars, with regard to collimation, spatial resolution, and interaction capability with atmospheric species, when compared with those of conventional microwave radar systems or passive visible instrumentation, hold promise that lidars will be long-lasting alternative observation systems. In an atmospheric lidar, the emission of a short laser pulse is followed by the reception of some radiation scattered from atmospheric constituents such as molecules, aerosols, and clouds. The interaction of the incident radiation with these constituents changes the intensity andyor the wavelength, depending on the strength of this optical interaction and the concentration of the interacting species. Consequently, it is possible to retrieve information about the physical state of the atmosphere along the exploration beam path. 1–3 In particular, estimation of the atmospheric optical parameters, namely, extinction and backscatter, based on pulsed elastic-backscatter lidars ~i.e., with no wavelength shift in reception!has been investigated in the literature. 4–6 The single-scattering range-return power for an elastic-backscatter lidar system can be expressed as 2 P~R!5A R2b~R!exp F 22 * 0 R a~r!dr G , (1) where P~R!is the range-received power ~W!,b~R!is the range-dependent volume backscatter coefficient of the atmosphere ~m 21 sr 21 !,a~R!is the rangedependent extinction coefficient ~m 21 !,Ris the range ~m!, and Ais the system constant ~Wm 3 !. Until now, the inversion of lidar signals has been tackled mainly by using classic procedures such as the slope method, 5 exponential-curve fitting, and Klett’s method. 6 Yet, all these methods assume simplifying andyor correlation hypotheses that limit the scope of the inversion results. These are discussed next. In the case of the slope-method algorithm, 5 the assumption of a homogeneous atmosphere is used to retrieve constant values ~a,b! as estimates of the The authors are with the Antennas, Microwaves, Radar and Optics Group, Department of Signal Theory and Communications, Universitat Polite`cnica de Catalunya, CySor Eulalia de Anzizu syn., 08034 Barcelona, Spain. Received 29 September 1997; revised manuscript received 8 April 1998. 0003-6935y98y307019-16$15.00y0 © 1998 Optical Society of America 20 October 1998 yVol. 37, No. 30 yAPPLIED OPTICS 7019 sought-after functions a~R!and b~R!. The key to the algorithm is a range-corrected function of the form S~R!5ln@R2P~R!#, (2) which enables us to find the extinction and backscatter coefficients from a linear regression of the form mina,biS~R!2@ln~Ab! 22aR#i2. (3) In the case of exponential-curve fitting, a similar approach is followed, but now the R 2 -corrected function is defined without the logarithm as F~R!5R2P~R!, (4) so that the norm minimization equivalent to expression ~3!takes the form mina,biF~R!2Abexp~22aR!i2. (5) Historically, this kind of fitting was introduced later because it is nonlinear in aand must be solved by using numerical methods. 7,8 In a different category, the inversion of the rangedependent function a~R!is first solved by Klett’s method. 6 The method assumes a power-law correlation between the extinction and backscatter atmospheric profiles as follows: b~R!5B0@a~R!#g, (6) and it requires a guess of the correlation constant g ~0.67 ,g,1!and a calibration at the far end of the inversion range interval in terms of S m 5S~R max ! and a m 5a~R max !. 9,10 Then the backward stable solution for a~R!becomes a~R!5exp@~S2Sm!yg# am2112 g * R Rmexp@~S2Sm!yg#dr . (7) In spite of the fact that this algorithm is significantly superior to the slope method and the exponential-curve fitting ~because the homogeneity approximation is not assumed!, the accuracy of the inverted profile a~R!is limited by that of the calibration a m and the correlation constant g. For this reason the algorithm retrieves a representative of the family a~R,g,a m !linked to the calibration pair ~g, a m !, which is thought to be close to the true extinction profile a~R!. In other words, the inversion of a range-dependent extinction profile from the return power is a many-to-one inversion problem ~see Appendix B in Ref. 11!, which can be solved only by adding appropriate a priori information ~e.g., calibrations along the observation path and physical constraints!. As we have seen, all these algorithms work with the present realization of the lidar return signal, so that correlation among past inverted returns remains unexplored. For example, in an elasticbackscatter pulsed lidar system, for each returnpower data stream received a new inversion, which is completely independent from those previously done, is performed. One of the things that distinguishes the Kalman filter 12–14 from nonmemory estimators such as those discussed above is the convenient way in which it accounts for any prior knowledge through a recursive process. As long as different power realizations are coming in, the filter updates itself, weighted by the imbalance between the a priori estimates of the optical parameters ~i.e., past inversions!and the new ones. Thus the new estimation of the optical parameters, or the project-ahead step ~a posteriori estimate!,isimproved based on a statistical minimum-variance criterion. In recent and pioneering work, Rye and Hardesty 15 and Lainiotis et al. 16 have found applications of the Kalman filter to the estimation of the return power and the logarithm of power for incoherent backscatter lidar with multiplicative noise. 17 Here we introduce an application of the filter to the solution of the inverse problem of joint estimation of the extinction and backscatter coefficients from the return power in an elastic-backscatter lidar. This paper is structured as follows: In Section 2 the problem is formulated from a theoretical point of view in terms of a first adaptive filter based on an extended Kalman filter ~EKF!; Section 3 describes the underlying statistics of the atmospheric model assumed in terms of the state-noise covariance matrix of the EKF; Section 4 discusses two examples of joint inversion of extinction and backscatter profiles from elastic-backscatter simulated lidar return signals; and Section 5 reviews some of the results presented by tackling a first real application of the filter to the inversion of power returns from a biaxial elastic-backscatter 1-J Nd:YAG lidar system. 2. Problem Formulation It is desirable now to study the feasibility of the derivation of the extinction and backscatter coefficients over the entire lidar inversion range. It is desired, then, to solve the functions a~R,t!and b~R,t!that, under a minimum-mean-square-error criterion, best fit the observable power P~R,t!at every time t. The term mean refers here to the ensemble average over time t.~Some revision of the EKF algorithm, along with the notation used below, is summarized in Appendix A.! A. State Vector Given the acquisition sampling rate of the system, f s , and considering the two-way path of the lidar signal, the power time samples P i correspond to a spatial sampling period DR5c 2fs . (8) Hence the spatial sampling points become Ri5Rmin 1~i21!DR,i51,...,N, (9) 7020 APPLIED OPTICS yVol. 37, No. 30 y20 October 1998 where R min is some predetermined minimum range of the system ~which is due to, for example, the minimum range of full overlap between the laser and the field of view of the receiving optics or some predefined minimum inversion range of interest!. The state vector to be estimated, x k , is a decimated version of the extinction and backscatter functions a~R!and b~R!over the whole lidar range. This is done in order to have more observables ~the power samples from each observation cell!than variables to estimate ~extinction and backscatter samples in the estimation cells!, and, as a result, it yields an overdetermined system with enhanced observability. 14,18 Furthermore, it can be shown 14 that if ~1!the system state vector is assumed to be a random constant, ~2! the measurement sequence z k yields an overdetermined set of linear equations, and ~3!the observation noise becomes negligible @i.e., the measurement noise covariance matrix of Eq. ~A4!,R k '0#, then the filter’s estimate depends more and more on current data ~which are rich in new information!and less and less on past inversions. Under these circumstances the filter behaves like a deterministic least-squares estimator ~which parallels the nonmemory approach!, and its estimate becomes 13 xˆk5Hk21zk, (10) where H k21 is the pseudoinverse matrix of H k @see relation ~A6!#. Historically, this bridges the gulf with past formulations of the problem in the form of expression ~5!. The model considers NyMobservation cells, where Mis the decimation ratio, so that the filter estimates NyMextinction samples and NyMbackscatter samples. Then the effective sampling period becomes MDR, which is Mtimes that of the return power. ~For simplicity, assume that Nis a multiple of M.! Mathematically, this can be expressed as ai5a~Rmin 1~i21!MDR!,i51,...,N M, bi5b~Rmin 1~i21!MDR!,i51,...,N M. (11) From these two halves of NyMelements, we form the state vector to be estimated: xk;~a1a2··· aNyMb1b2··· bNyM!T, (12) where the subscript kis a reminder of the discrete time t k . The nonstationarity or the dynamics of the state vector is described by the transition matrix F k @see Eqs. ~A1!and ~A19!# and the state-noise covariance matrix Q k @Eq. ~A2!#. The former represents how the state vector projects ahead from time t k to time t k11 , and the latter gives the filter key information about the underlying statistics of each component of the state vector ~the optical parameters under study!. Formulation of the system equations in terms of the transition matrix F k and the state-noise covariance matrix Q k gathers all the information the filter knows about the atmospheric model. This is tackled in Section 3. B. Measurement Equation If, according to Eq. ~9!, each power sample corresponds to a spatial increment DR, so that P i 5P~R i !, and a rectangle approximation is used to compute the transmittance term of Eq. ~1!, the observable power samples become P15A R12b1exp~22a1Rmin!, (13) · · · PM5A RM2b1exp$22a1@Rmin 1~M21!DR#%, (14) · · · PM115A RM112b2exp$22a1@Rmin 1~M21!DR# 22a2DR%, (15) · · · PN5A RN2bNyMexp H 22a1@Rmin 1~M21!DR# 22 ( i52 NyM aiMDR J . (16) The R 2 -corrected version @Eq. ~4!# of this set of N equations builds the measurement vector z k : zk;@F1~xk!F2~xk!··· FN~xk!#T. (17) By using F~R!rather than P~R!, we reduce the dynamic margin of z k and, consequently, numerical errors are reduced as we cycle through the Kalman loop @Eqs. ~A13!–~A17!#. At the far ranges, however, amplification of quantization noise generated during the analog-to-digital conversion might become significant and must be accounted for in a description of the statistics of the observation noise. Equations ~13!–~16!define the overdetermined set of equations discussed in Subsection 2.A that relates the measurement vector z k to the a priori estimate of the state vector, xˆ k 2 . From Eqs. ~12!and ~17!, the N3~2NyM!observation matrix H k can be computed by splitting it in two N3~NyM!submatrices of the form H5~H 1 H 2 !, where Hij ~1!5]Fi ]aj U x5xˆk 2 ,Hij ~2!5]Fi ]bj U x5xˆk 2. (18) 20 October 1998 yVol. 37, No. 30 yAPPLIED OPTICS 7021 This yields where H 1 and H 2 are evaluated at the a priori estimate xˆ k 2 . From the structure of H 1 and H 2 , it emerges that the formulation of the inversion problem involves, however, a trade-off between larger decimation ratios ~M! and model accuracy. From the point of view of M,if one compares the EKF formulation with the classical exponential-curve-fitting counterpart of expression ~5!, the exponential fitting algorithm works with a very large value of M, equal to the length of the inversion interval, which is, in turn, formed by a single inversion cell. Hence the larger the M, the more robust the system of Eq. ~10!against observation noise and the better the regression results @this is best seen by the products MDRF i in Eq. ~19!, which become larger for larger M#. Another advantage of increasing Mis the enhancement of the filter’s sensitivity to low return powers and large values of R min . Usually, large differences between the weight factors 22@R min 1~M2 1!DR#and 22MDRare not desirable, since then thetrajectory of the filter would be dominated by the estimation of the first cell. Yet, the most risky drawback for large Marises from the deterioration of the filter’s model @Eqs. ~13!–~16!#: Although the equivalent sampling period of the optical parameters in the filter’s model is MDR, the true atmospheric spacing DR9is differential in nature. For the time being, we assume the simplification DR95 DRin the atmospheric optical profile, so that such modeling errors are neglected. In other words, the atmosphere is assumed homogeneous inside any observation cell. Although much research is being done in this field, which is far from the scope of this study, a sensible solution might be achieved by combining the results of an array of Mcooperative filters or, perhaps, by using nonuniform spacing in the formulation of the problem, depending on the atmospheric situation at hand. H15 3 22RminF10 0 ··· 0 22~Rmin 1DR!F20 0 ··· 0 · · ·· · ·· · ···· · · · 22@Rmin 1~M21!DR#FM0 0 ··· 0 22@Rmin 1~M21!DR#FM1122DRFM110 ··· 0 · · ·· · ·· · ···· · · · 22@Rmin 1~M21!DR#FN22MDRFN22MDRFN··· 22MDRFN 4 N3~NyM! , (19) H25      F1 xNyM11 0 0 ··· 0 F2 xNyM11 0 0 ··· 0 · · ·· · ·· · ···· · · · FM xNyM11 0 0 ··· 0 0FM11 xNyM12 0 ··· 0 · · ·· · ·· · ···· · · · 0 0 0 ··· FN x2NyM      N3~NyM! , (20) 7022 APPLIED OPTICS yVol. 37, No. 30 y20 October 1998 Lidar measurements are corrupted mainly by Gaussian additive observation noise v k . Assuming that DR95DRand that the observation noise along the inversion range can be approximated by rangedependent stationary electronic thermal noise ~Gaussian additive noise with variance s r 2 !, the observation noise covariance matrix is computed as Rk5E~vkvkT!5 3 sr2~R1!R14··· 0 · · ····· · · 0 ··· sr2~RN!RN4 4 . (21) This assumes that electronic noise dominates observation noise. Following Refs. 17 and 19, the rangedependent noise variance can be written as sr2~R!5a@P~R!1Pback#1b, (22) where P~R!is the range return power defined in Eq. ~1!,P back is the background power from any other interfering source ~for example, the Sun!, and aand bare constants that depend only on specific parameters of the receiving system. The first term accounts for the contributions of the signal-induced shot noise to the total noise, and the second one merges into the variable bthe contributions of both dark-current shot noise and thermal noise. @Note that s r 2 ~R!has units of square volts or square watts, depending on whether equivalent noise is computed at the receiver’s output or input. In instances where other sources of measurement noise are present ~e.g., R 2 -amplified quantization noise!, one can increase pertinent terms along the main diagonal of R k to accommodate such an extra variance. Nonstationary noise can be tackled by recomputing R k at each succeeding step of the filter, and colored noise ~such as synchronized flash-lamp interferences!can be modeled by also using elements off the main diagonal of R k ~see also Ref. 14 for further insight!. 3. Atmospheric Model for the Extended Kalman Filter In their most general form, extinction and backscatter optical parameters are nonlinearly related. Unless microscale analysis is considered and plenty of boundary calibrations from balloon-borne instrumentation or other cooperative systems are given, the struggle to model physically the temporal and spatial evolution of the optical parameters leads to awkward and cumbersome results. A more convenient alternative is to try to model the macroscopic effects on them, rather than the underlying microphysical parameters. This is done by using the time–space stochastic correlation model sketched in Fig. 1. Each optical component ~extinction and backscatter!of an inversion cell is modeled as a stochastic process having both temporal and spatial correlation. Each output branch represents one component of the state vector x k , and the vector noise process w k is formed by spatially correlated components at the output of the linear system A. The atmospheric model is driven by an array of white-noise uncorrelated processes. A. Temporal Correlation Temporal correlation is perhaps the most attractive advantage of the EKF lidar inversion approach over the nonmemory solutions of expressions ~3!and ~5! and Eq. ~7!. This advantage comes from telling the filter that it should improve its projection steps based on the fact that the atmosphere usually has a long correlation time and that, consequently, swift changes in any optical parameter are not possible. Temporal correlation is achieved by modeling each component of the state vector x k ~with kthe discrete time!as a Gauss–Markov process. ~To simplify the notation, we define the Markovian process y k as the ith component of vector x k , so that y k 5x i,k .! The Gauss–Markov process, 14 which is often called Markovian noise, is zero-mean low-pass filtered Gaussian noise, whose autocorrelation function is given by Ry~t! 5sm2exp~2dutu!, (23) where s m 2 is the power of the process y k and dis the 3-dB cutoff frequency of the low-pass coloring filter ~IIR boxes in Fig. 1, where IIR stands for infinite impulse response!. The discrete-time equation of the process can be written in the form of an autoregressive movingaverage scalar process 20 as yk115exp~21yLc!yk1wk, (24) where y k and w k are the Markovian and white sequences, respectively, and L c is the temporal correlation length, defined as Lc51yd, (25) where L c has units of samples @the spatial period has already been defined in Eq. ~8!#. Finally, Eqs. ~24!and ~25!enable us to express the Fig. 1. Time–space EKF correlation model. 20 October 1998 yVol. 37, No. 30 yAPPLIED OPTICS 7023 state-vector transition matrix associated with Eq. ~A19!as Fk5exp~21yLc!I, (26) where Iisa~2NyM!3~2NyM!identity matrix and we have used the simplifications that F k is constant over time t k and that L c is the same for all the cells along the lidar exploration path. With F k a matrix, w k also becomes a vector, whose 2NyMcomponents represent white sequences at each time t5t k . In practice, Markovian noise is responsible for the time drift of the actual value of the optical parameters being estimated by the filter at each projection step. For this reason it is useful to define an intensity parameter pthat enables us to adjust the power of the Markovian process or, equivalently, to link the driving Gaussian noise power to the amplitude change caused in an optical parameter. From Eqs. ~23!and ~25!and Ref. 14, the Gaussian noise standard deviation s w and the Markovian one, s m , can be related as sw5sm Î 12exp~22yLc!. (27) A comprehensive collection of histograms for a given s m have shown that over 95% of the Markovian amplitudes distribute between 62.5s m . With that in mind, the white-noise strength s w 5s a i that is needed to cause a p-per-one change in the amplitude of a vector component of x k ~let it be a i !over a correlation length L c becomes sai5p 2.5 ai Î 12exp S 22 Lc D ,i51,...,N M. (28) B. Spatial Correlation Contrary to what happened with nonmemory algorithms, where analytical correlation relations were assumed @for example, homogeneity for the slope and exponential-curve-fitting algorithms of expressions ~3!and ~5!, respectively, or the power-law correlation of Eq. ~6!for Klett’s method#, the approach presented here is based on the correlation graph of Fig. 2. It paves the way for the introduction of loose stochastic relations among the sought-after optical parameters instead of tight analytical ones. It seems sensible to guess that, for example, any extinction change in a particular cell will, in turn, influence variations not only in the in-cell backscatter component but also in the extinction and backscatter components of its neighboring cells. The underlying physical phenomenon being the cause, the changes may well extend over several cells. From the correlation graph of Fig. 2, one can build the white-noise state-vector covariance matrix as follows: Cw5 F Caa Cab Cba Cbb G . (29) These block matrices can be developed as Caa 5 3 sa1 2rsa1sa2··· rn21sa1san ··· sa2 2··· rn22sa2san ··· ··· ··· ··· ··· ··· ··· san 2 4 , (30) Cab 5 3 r9sa1sb1r9rsa1sb2··· r9rn21sa1sbn ··· r9sa2sb2··· r9rn22sa2sbn ··· ··· ··· ··· ··· ··· ··· r9sansbn 4 , (31) where C ba 5C ab ,C bb is the same as C aa but with b and apermuted, ris the correlation coefficient between one cell and the next one along the beam path ~which is due to the physical continuity of the atmosphere!,r9 is the in-cell extinction-to-backscatter correlation coefficient, and s i has been defined above in Eq. ~28!. On the condition that uru,1, ur9u,1, it can be proved that the graph of Fig. 2 does represent a covariance matrix. Assuming that temporal and spatial correlation processes are independent, the Markovian noise state-vector covariance matrix ~i.e., the sought-after covariance matrix for the atmospheric model given to the filter, Q k !can be computed from Eqs. ~27!,~29!, ~30!, and ~31!as Qk5Cw 12exp~22yLc!, (32) where Q k and C w are basically the same except for a scaling factor, which could, in turn, be merged into an equivalent intensity parameter p9. Although the inversion of wind fields is far from the objective of this study, this independent hypothesis between temporal and spatial correlation is, however, doubtful in situations with a significant radial wind component ~i.e., the wind component along the exploration path!. In these instances radial wind strongly correlates both space and time fluctuations along the line of sight. Here one might consider only the spatial correlation of Fig. 1 ~i.e., Q k 5C w,k !, but this time a variant one, since the observation cells along the path become progressively affected by different correlation links as time goes on. In addition, boosting elements off the main diagonal would tell the filter of a significant increase in the correlation among neighboring cells. In any case the possibility of modeling nonstationary statistics in Q k by resetting it at each succeeding step of the filter offers a wide span of attractive possibilities yet to be investigated. Fig. 2. Spatial correlation graph of the state-vector components. 7024 APPLIED OPTICS yVol. 37, No. 30 y20 October 1998 Usually, the state-noise covariance matrix of the EKF, Q k , is the most difficult input to assess, since its atmospheric counterpart Q k,a ~the subscript arefers to atmospheric!is unknown. The problem of finding good models for Q k has sometimes been tackled by using a partitioned approach, 16,21,22 where the unknowns are merged into a vector of random variables Q5~u 1 ...u P !with known or assumed a priori probability density functions ~consider, for example, Ref. 16!. This approach yields a bank of EKF’s working in parallel, each matched to an appropriate value of u i , so that the overall vector Qspans the space of unknown parameters constrained by their related possible values. In theory, joint estimation of the extinction and backscatter parameters Qwould virtually apply to all the elements of Q k,a plus, possibly, the equivalent intensity parameter p9. In practice, this would involve a large array of filters that would probably exceed the framework of intelligent and selforganizing systems, and hence it would certainly prevent a straightforward formulation of the study. For this reason, at this first stage, estimation with a single EKF is preferred in this experimental work, even though this is done at the expense of larger modeling errors and, hence, worse performance. These model uncertainties justify a formulation of the a priori error covariance matrix as P0 25mQ0,m$1. (33) With regard to the atmospheric model, the simulations have used a set of parameters Q k,a ,L c,a , and p a different from those given to the EKF model, Q k , L c , and p, to test the performance of the filter under modeling errors. Eigenvalue decomposition is used to compute the linear correlator ~Ain Fig. 1!and the power of the white-noise uncorrelated sequences ~n 1 ...n NyM !, which are the driving inputs of the atmospheric simulator. 14 4. Simulation Results Through extensive simulation and simplified conditions, joint estimation of extinction and backscatter simulated profiles from elastic-backscatter lidar return signals have been inverted by using the formulation presented above. Next, two simulation sets are discussed; the first one ~Figs. 3–7!corresponds to a good-visibility scene, and the second one ~Figs. 8–10!corresponds to moderate-visibility conditions. Simulation parameters are summarized in Table 1. First, we give a brief outline of the choice of statistical parameters. The choice of r a and r9 a is based on cross-examined time–space plot sets of synthesized power return signals with nonwindy time–space real observations. Good agreement between typical real data sets and simulated ones has usually been achieved for large values of r9 a ~typically between 0.8 and 0.9!and medium values of r a ~typically between 0.3 and 0.7!. The former result is also in accordance with Eq. ~6!, where g51 is equivalent to r9 a 31. As for the latter, it has been found that r a values close to unity are not advisable because they yield stiff spatial profiles that are so correlated that it is difficult to accommodate even moderate heterogeneities along the lidar path. It has also been found that the intensity parameter p a is the most critical of all and that it must be adjusted to each particular scene. As a rule of thumb for low atmospheric extinctions, measurement of the fluctuations in the range-corrected power has yielded acceptable estimations of p.L c is usually determined from rough visual estimation. The first simulated set is related to a mean visibility of V M 539.12 km. Such visibility conditions are typical of standard clear to exceptionally clear air. From Refs. 23 and 24 and under the approximation of a homogeneous atmosphere, the visibility parameter can roughly be linked to the atmospheric optical parameters a50.1 km 21 and b54310 23 km 1 sr 21 or, equivalently, k a 525 sr and b54310 23 km 21 sr 21 , where k a is the extinction-to-backscatter ratio indicated in Table 1 and the subscript arefers to atmosphere. Hence one can without distinction talk about visibility or homogeneous atmospheric optical parameters ~a,b!. To simulate an inhomogeneous profile approximately related to the visibility V M , the simulator computes a range-dependent hump-shaped backscatter profile with mean b, such as the one shown in Fig. 3~a!. For other visibility margins, the profile is scaled accordingly. In this way it is ensured that the synthesized profile is always approximately related to the average visibility desired. We computed the Table 1. Simulation Parameters Basic parameters Optical parameters ~set 1!a50.1 km 21 ,b54310 23 km 21 sr 21 ,V M '39.12 km ~set 2!a51km 21 ,b53310 22 km 21 sr 21 ,V M '3.91 km Inversion range @Eq. ~9!# R min 5200 m, R max 55 km, DR5123.1 m @Eq. ~8!# Order parameters N540, M52@Eqs. ~11!#, iterations 5320 System constant @Eq. ~1!# A52.35 310 23 Wkm 3 Noise parameters @Eq. ~22!# a51.8 310 210 W, b55310 218 W 2 ,P back '2nW Model parameters @Eqs. ~28!–~32!# Atmosphere k a 5ayb,p a 540%, L c,a 550, r a 50.6, r9 a 50.9 EKF k50.9k a ,p550%, L c 5100, r50.3, r9 5 0.8 Initialization P 0 2 5Q 0 @Eq. ~33!#,xˆ 0 2 5~kb,...,kb,b,...,b! 20 October 1998 yVol. 37, No. 30 yAPPLIED OPTICS 7025 range-dependent extinction profile after the backscatter by reusing the extinction-to-backscatter ratio k a . From the extinction and backscatter initial profiles just computed, the lidar range-return power of Eq. ~1!follows as shown in Figs. 3~b!and 3~c!.To reduce the order of the filter, the physical problem has been discretized by using N540 and M52, which means N540 power samples and NyM520 observation cells. Since each cell is defined by both its extinction and backscatter parameters, there are 40 state-vector components, 20 for each optical parameter. To compute the observation noise, the simulator uses electrical and optical parameters from an elastic-backscatter lidar of the Polytechnic University of Catalonia in Barcelona, Spain ~system specifications are given in Section 5!to assess realistic noise parameters in Eq. ~22!. They are representative of a typical tropospheric lidar system ~see Table 1!and yield the range-dependent signal-to-noise ratio ~SNR!of Fig. 3~d!. The atmospheric behavior was simulated by using the simplified model described in Section 3 and the model parameters of Table 1. As for the EKF, a slightly mismatched model is input, so that, for example, the extinction-to-backscatter ratio is underestimated by 10%, the temporal correlation length is doubled, and the spatial correlation coefficients are changed as indicated in Table 1. These modeling errors translate into Q k,a ÞQ k , as suggested in Subsection 3.B. Since we are particularly concerned about the performance of the filter under different visibility conditions and atmospheric modeling errors, the simplification in which there are no mismatches in the model of R k , so that both the observables and the filter share the same covariance matrix, has been assumed. This can be justified because Q k,a is always the hidden parameter of the atmosphere, whereas R k can ultimately be measured from the lidar system. The initialization of the filter, xˆ 0 2 , may come from any of the methods discussed in Section 1; in particular, Eq. ~7!would yield the best approximation. Yet, to test the performance of the filter, it has been initialized in the simplest possible way by using a constant homogeneous profile for the extinction and backscatter components of the state vector, as indicated in Table 1. Since each simulation run takes 320 iterations, the filter depends more and more on the measurements and less and less on the initial state. As time goes on, the actual measurement data ~observables!received for any particular sample run change according to the atmospheric state model given by F k and Q k,a , so that slowly varying changes in both the extinction and backscatter profiles are accommodated. The filter keeps track of the timevarying nature of the observables from the beginning. Figures 4 and 5 illustrate the time evolution of the atmospheric model along with the EKF state-vector components. Recall that components 1–20 represent the extinction coefficient and components 21–40 represent the backscatter coefficient along with the observation cells, so that if one reads by cells, the first one comprises components 1 and 21, the second one Fig. 3. ~Set 1!initial state of the simulation: ~a!synthesized backscatter profile, ~b!range-corrected return power, ~c!return power as received by the lidar, ~d!associated SNR. Fig. 4. ~Set 1!time–space evolution of the extinction and backscatter profiles: ~a!synthesized atmospheric optical parameters ~extinction and backscatter!,~b!EKF inverted optical parameters. 7026 APPLIED OPTICS yVol. 37, No. 30 y20 October 1998 comprises components 2 and 22, and so on. The temporal evolution of the mountains in Fig. 4 is caused by the Markovian noise. Spatially, with r9 a 5 0.9 the in-cell extinction-to-backscatter correlation is so high that the two halves of each plot look virtually alike ~note that for illustrative purposes the backscatter half has been rescaled by k a !. Convergence of the EKF from the homogeneous initial profile to something close to a real profile can easily be tracked by monitoring the trace of the error covariance matrix P k @see Eq. ~A18!in Appendix A and Section 5 for further insight#. Since P k informs the filter about the expected error that it is committing at each time t k , a good convergence criterion is whether P k has reached a constant value. In the plots presented, as is always the case, the shape of the estimated profiles is retrieved fast, but their magnitudes ~especially the extinction one!take some more time to settle. In the simulations the trace of P k settles by iteration 150, although after iteration 50 most details from the true atmospheric profile at short ranges are recovered quite well. Figure 5 is a contour plot of Fig. 4 representing isoextinction and isobackscatter curves ~scaled by k a for illustrative purposes!along time for both the atmospheric and the estimated state vector. Both contours look virtually alike after the 50th iteration except for some slight deterioration in the far-range extinction components of the atmosphere ~components 10–20!, where the EKF performs more poorly. This, however, can easily be justified by the progressive reduction in the SNR of Fig. 3~d!for increasing ranges. From the point of view of the time–space correlation models introduced in Subsections 3.A and 3.B, Fig. 6 compares the time evolution of the EKF estimates in four observation cells successively farther along the lidar exploration range ~cells 5, 10, 15, and 20 located at 1307.7, 2538.5, 3769.2, and 5000 m, respectively!with their true atmospheric counterparts. The atmospheric backscatter evolution is denoted by solid curves, and the filter’s estimates are given by small circles. The filter follows the random drift of each cell fairly well in all the cases, but whereas Figs. 6~a!and 6~b!show the best-fitted cells, Figs. 6~c!and 6~d!show some slight underestimation of the atmospheric backscatter. Horizontal solid lines indicate the initial backscatter value in each cell before the atmospheric simulator starts. These values correspond to the 5th, 10th, 15th, and 20th components of Fig. 3~a!. As expected from the temporal correlation model formulation of Subsection 3.A, Markovian noise translates into a slow temporal drift of the backscatter figure. For example, the atmospheric temporal correlation length ~L c,a 550 samples!is best seen in Figs. 6~a!and 6~c!~solid curves!. Thus, in Fig. 6~c!, increasing and decreasing slopes last for approximately 50 samples on average, and the same happens in Fig. 6~a!except that now the plot includes some kind of horizontal interval. From the point of view of the spatial correlation, one has to compare all the plots. Since r a 50.6 and each plot Fig. 5. ~Set 1!contour plots of Fig. 4 showing very good correlation between the time–space evolution of the atmospheric optical parameters and the inverted ones: ~a!synthesized atmospheric optical parameters, ~b!EKF inverted optical parameters. Fig. 6. ~Set 1!temporal evolution of the backscatter profiles in four representative observation cells along the lidar beam path: ~horizontal lines!starting backscatter values for the atmospheric simulator, ~solid curves!atmospheric backscatter evolution, ~circles!EKF estimates. 20 October 1998 yVol. 37, No. 30 yAPPLIED OPTICS 7027 These expressions bridge the gulf with the classical linear filter if the equivalent matrices F k and H k are defined in the following way: Fk5]fk~x! ]x U x5xˆk , (A7) Hk5]hk~x! ]x U x5xˆk 2 . (A8) Identification with the first-order terms of approximations ~A5!and ~A6!yields xk11<fk~xˆk!1Fk~xk2xˆk!1wk, (A9) zk<hk~xˆk 2!1Hk~xk2xˆk 2!1vk. (A10) Approximations ~A9!and ~A10!represent the linearized version of the filter and resemble those of a linear Kalman filter except for the fact that rather than presenting total quantities to the filter, we consider incremental ones. In relation to approximations ~A9!and ~A10!, these are Dxk5xk112fk~xˆk!, (A11) Dzk5zk2hk~xˆk 2!. (A12) In summary the EKF’s recursive equation set becomes xˆk5xˆk 21Kk@zk2hk~xˆk 2!#, (A13) Pk5~I2KkHk!Pk 2, (A14) xˆk11 25fk~xˆk!, (A15) Pk11 25FkPkFkT1Qk, (A16) Kk5Pk 2HkT~HkPk 2HkT1Rk!21, (A17) where K k is the Kalman gain and P k 2 is the associated error covariance matrix, defined as Pk 25E~ek 2ek2T!5E@~xk2xˆk 2!~xk2xˆk 2!T#, (A18) where e k 2 is the a priori estimation error. Yet, careful attention should be drawn to the fact that use of the EKF is risky, as the linearization process takes places about the filter’s estimated trajectory of the state vector rather than about a precomputed nominal trajectory. That is, the partial derivatives are evaluated along a trajectory that has been updated with the filter’s estimates; thus it depends on the measurements. As a result, the filter is more likely to diverge. In the EKF problem formulated in this work, the system model is linear and the state-space representation of the atmospheric state vector is given by xk115Fkxk1wk, (A19) where F k is the transition state matrix from time t k to time t k11 . If both the system and the observation model are linear, Eqs. ~A13!–~A17!become the same after we replace F k by F k and h k by H k . We acknowledge the sponsorship of the CICYT ~Interministry Committee for Science and Technology! under grant AMB96-1144-C02-C01. References 1. R. T. H. Collis and P. B. Russell, “Laser measurement of particles and gases by elastic backscattering and differential absorption,” in Laser Monitoring of the Atmosphere,E.D. Hinkley, ed. ~Springer-Verlag, New York, 1976!, Chap. 4, pp. 91–102. 2. R. M. Measures, “Laser-remote-sensor equations,” in Laser Remote Sensing: Fundamentals and Applications ~Krieger, Malabar, Fla., 1992!, Chap. 7, pp. 237–280. 3. D. K. Killinger and N. Menyuk, “Laser sensing of the atmosphere,” Science 235, 37–45 ~1987!. 4. A. I. Carswell, “Lidar remote sensing of atmospheric aerosols,” in Propagation Engineering: Third in a Series, L. R. Bissonnette and W. B. Miller, eds., Proc. SPIE 1312, 206–220 ~1990!. 5. G. J. Kunz and G. de Leeuw, “Inversion of lidar signals with the slope method,” Appl. Opt. 32, 3249–3256 ~1993!. 6. J. D. Klett, “Stable analytical inversion solution for processing lidar returns,” Appl. Opt. 20, 211–220 ~1981!. 7. R. J. Barlow, Statistics ~Wiley, New York, 1989!. 8. J. J. More, “The Levenberg–Marquardt algorithm: implementation and theory,” in Numerical Analysis, Lecture Notes in Mathematics 630, G. A. Watson, ed. ~Springer-Verlag, New York, 1977!, pp. 105–116. 9. J. D. Klett, “Lidar calibration and extinction coefficients,” Appl. Opt. 22, 514–515 ~1983!. 10. J. D. Klett, “Lidar inversion with variable backscatteryextinction ratios,” Appl. Opt. 24, 1638–1643 ~1985!. 11. G. J. Kunz, “Probing of the atmosphere with lidar,” in Proceedings of Remote Sensing of the Propagation Environment ~AGARD-CP-502!,23, 1–11 ~1992!. 12. R. E. Kalman, “A new approach to linear filtering and prediction problems,” J. Basic Eng. 82, 35–46 ~1960!. 13. H. W. Sorenson, Kalman Filtering Techniques. Advances in Control Systems. Theory and Applications ~IEEE, New York, 1985!, Vol. 3. 14. R. G. Brown and P. Y. C. Hwang, Introduction to Random Signals and Applied Kalman Filtering ~Wiley, New York, 1992!. 15. B. J. Rye and R. M. Hardesty, “Nonlinear Kalman filtering techniques for incoherent backscatter lidar: return power and log power estimation,” Appl. Opt. 28, 3908–3917 ~1989!. 16. D. G. Lainiotis, P. Papaparaskeva, G. Kothapalli, and K. Plataniotis, “Adaptive filter applications to LIDAR: return power and log power estimation,” IEEE Trans. Geosci. Remote Sens. 34, 886–891 ~1996!. 17. R. J. McIntyre, “Multiplication noise in uniform avalanche photodiodes,” IEEE Trans. Electron Devices ED-13, 164–168 ~1966!. 18. P. S. Maybeck, Stochastic Models, Estimation and Control ~Academic, New York, 1977!, Vol. 1. 19. W. B. Jones, Introduction to Optical Fiber Communication Systems ~Holt, Rinehart & Winston, New York, 1988!, Chap. 7, 8. 20. A. Papoulis, Probability, Random Variables and Stochastic Processes ~McGraw-Hill, New York, 1991!. 21. D. G. Lainiotis, “Partitioned estimation algorithms. I: Nonlinear estimation,” J. Inf. Sci. 7, 203–255 ~1974!. 22. D. G. Lainiotis, “Partitioning: a unifying framework for adaptive systems. I: Estimation,” Proc. IEEE 64, 1126– 1143 ~1976!. 23. H. Koshmieder, “Theorie der Horizontalen Sichtweite,” Beitr. Phys. Freien Atmos. 12, 33–53 ~1924!. 24. P. W. Kruse, L. D. McGlauchlin, and R. B. McQuiston, Elements of Infrared Technology: Generation, Transmission and Detection ~Wiley, New York, 1962!. 7034 APPLIED OPTICS yVol. 37, No. 30 y20 October 1998