The estimation of the source pollutant strength is a relevant issue for atmospheric environment. This characterizes an inverse problem in the atmospheric

pollution dispersion studies. In the inverse analysis, a time-dependent pollutant source is considered, where the location of such source term is assumed known. The inverse problem is formulated as a non-linear optimization approach, whose objective function is given by the least-square difference between the measured and simulated by the mathematical model, pollutant concentration, associated with a regularization operator. The forward problem is addressed by a Lagrangian model, and a quasi-Newton method is employed for minimizing the objective function. The

second-order Tikhonov regularization is applied and the regularization parameter is computed by using the L-curve scheme. The inverse-problem methodology is verified

with data from the tracer Copenhagen experiment

Anfossi, D.

de Campos Velho, H. F.

Degrazia, G. A.

Roberti, D. R.

English

Scientific Open-access Literature Archive and Repository

DOI 10.1393/ncc/i2007-10003-yIL NUOVO CIMENTO Vol. 30 C, N. 2 Marzo-Aprile 2007Estimation of emission rate from experimental dataD. R. Roberti(1)(∗), D. Anfossi(2)(∗∗), H. F. de Campos Velho(3)(∗∗∗)and G. A. Degrazia(1)( ∗∗∗)(1) Department of Physics, Federal University of Santa Maria (UFSM)Santa Maria (RS), Brazil(2) ISAC, CNR - Torino, Italy(3) LAC, INPE - Sa˜o Jose´ dos Campos (SP), Brazil(ricevuto l’ 1 Marzo 2007; approvato il 30 Aprile 2007; pubblicato online il 16 Luglio 2007)Summary. — The estimation of the source pollutant strength is a relevant issue foratmospheric environment. This characterizes an inverse problem in the atmosphericpollution dispersion studies. In the inverse analysis, a time-dependent pollutantsource is considered, where the location of such source term is assumed known.The inverse problem is formulated as a non-linear optimization approach, whoseobjective function is given by the least-square difference between the measured andsimulated by the mathematical model, pollutant concentration, associated with aregularization operator. The forward problem is addressed by a Lagrangian model,and a quasi-Newton method is employed for minimizing the objective function. Thesecond-order Tikhonov regularization is applied and the regularization parameter iscomputed by using the L-curve scheme. The inverse-problem methodology is verifiedwith data from the tracer Copenhagen experiment.PACS 47.27.tb – Turbulent diffusion.1. – IntroductionThe pollutant dispersion in the atmosphere is a topic with high interest nowadays.Models for air monitoring are not only important for describing the pollutant impactover urban or rural areas, or for urban planning consideration, but other issues are alsorelevant, such as the pollutant source strength estimation, CO2 diurnal cycle, and to-tal ozone in the atmosphere. The last three issues are examples of inverse problems in(∗) E-mail: droberti@smail.ufsm.br(∗∗) E-mail: anfossi@isac.cnr.it(∗∗∗) E-mail: haroldo@lac.inpe.br( ∗∗∗) E-mail: degrazia@ccne.ufsm.brc© Societa` Italiana di Fisica 177178 D. R. ROBERTI, D. ANFOSSI, H. F. DE CAMPOS VELHO and G. A. DEGRAZIAatmospheric pollution modelling. Inverse problems are usually solved by an implicit tech-nique, where the inverse problem is formulated as a constrained non-linear optimizationproblem: the forward problem is iteratively solved for successive approximations of theunknown parameters. The associated forward problem is the solution of the dispersionmodel. As pointed out in [1], the major obstacle to a good inversion result is the accuracyof the dispersion model and the representativeness of the measurements.The purpose of this paper is to employ an inverse-problem methodology for estimatingthe emission rate from a measured dataset in Copenhagen experiment using a Lagrangianparticle model, named LAMBDA [2,3].2. – Forward modelThe Lagrangian particle model LAMBDA was developed to study the transport pro-cess and pollutants diffusion, starting from the Brownian random walk modelling. Thismodel is based on a three-dimensional form of the Langevin equation for the randomvelocity [4]. LAMBDA was tested against many data sets, and in particular with theCopenhagen experiment, giving satisfactory results [2, 3].LAMBDA uses a large number of fictitious particles to simulate the atmosphericdiffusion. Each particle can be marked for a mass, and the spatial distribution of particlesin the computational domain allows calculating the three-dimensional concentration field,CMod(r ) = C(x, y, z, t), through the calculation of how many of them lie in a cell orimaginary volume centred in x, y, z, as follows:(1) C(x, y, z, t) = mpNvVc,where mp is the mass of each particle, Nv is the particle number in the cell and Vc is thecell volume. The mass of each particle is determined from(2) mp =Q(t)ΔtNp,where Q(t) is the emission rate, Δt is the time step and Np is the number of particleemitted per time step.3. – Inverse modelIn order to set up the inverse analysis, it is assumed that the concentration obtainedwith the mathematic model is given by CMod(r,Q), where Q = Q(t) is the emissionrate and CExp(r ) are data from experimental data of concentration. The solution of theinverse problem is a function Q that minimizes the following objective function:(3) J(λ,Q) =∥∥CExp(r )− CMod(r,Q)∥∥22+ λΩ(Q) .The least-square difference between experimental data and calculated values is repre-sented by the first term in eq. (3); Ω(Q) is a regularization operator, a a priori piece ofinformation, that mathematically is a restriction about the function.ESTIMATION OF EMISSION RATE FROM EXPERIMENTAL DATA 179The regularization operator can be expressed by Tikhonov scheme [5]:(4) Ω(Q) =p∑m,j=0κm,j∥∥∥Q(m)∥∥∥22;here Q(m) denotes the m-th difference. In general the parameter κm,j is chosen asκm,j = δmj (Kronecker delta) and the regularization is named Tikhonov-j regularizationoperator, where j denotes the order of the regularization. Here, the second-order regu-larization (Tikhonov-2) is employed.3.1. Optimization algorithm. – The optimization problem is iteratively solved by thequasi-Newtonian optimiser routine E04UCF, from the NAG Fortran Library. This algo-rithm is designed to minimize an arbitrary smooth function subject to constraints (simplebound, linear or non-linear constraints), using a sequential programming method.This routine has been successfully used in several previous works: in geophysics,hydrological optics, heat transfer, and meteorology [6-12].4. – Simulating the Copenhagen experiments with forward problemThe choice of this experiment derives from the amount and quality of the measure-ments. For each experiment, meteorological information, like wind velocity standarddeviations at the emission height, mixing height and Obukhov length, were available foreach hour; and wind speed and wind direction were available at each 10min, and atthree levels along a tower. The ground level concentration was measured using an arrayof samplers.Copenhagen tracer dispersion experiments were carried out in the northern part ofCopenhagen [13,14]. In this experiment, the pollutant (SF6) was released without buoy-ancy from a tower at a height of 115m, and collected at ground level in three crosswindarcs located at 2–6 km from the releasing point. The site was mainly residential witha roughness length of 0.6m. The Copenhagen region lies on an island. Tracer concen-trations were measured and an average was computed at each 20min. We consideredthe exercise performed on 19 October 1978, and the concentrations evaluated in the sec-ond period (12:33h-13:53h). The Copenhagen tracer dispersion experiment used in thisanalysis was selected due to the amount and quality of the concentration and meteoro-logical data set.Meteorological data, available each 10min (see table I), were used to create the inputfor the simulation in the forward model (LAMBDA). The profiles of wind standard de-viations (σi, being i = u, v, w), and the Lagrangian decorrelation time scales (τLi) werecalculated according to the turbulence parameterisation scheme suggested by Degrazia etal. [15] (figs. 1a and b). The measured wind velocity variances σ2i (i = v, w) were used tofit the analytical model for the wind variance. It is worth mentioning that the ratios ofthese variances, (σ2v)Model/(σ2v)Exp|z=10 m and (σ2w)Model/(σ2w)Exp|z=10 m, were found tobe very close to one (1.08, 1.05, respectively), and this gives a further verification of thegood quality of the parameterisations used. The longitudinal wind standard deviationwas set equal, at each height, to the crosswind one. It was verified that the observedvertical profile of the wind speed U(z) was in a very good agreement with the classical180 D. R. ROBERTI, D. ANFOSSI, H. F. DE CAMPOS VELHO and G. A. DEGRAZIATable I. – Meteorological data.U (m/s) θ (◦)hour 10m 120m 200m 10m 120m 200m12:05 2.6 5.7 5.7 290 310 31012:15 2.6 5.1 5.7 300 310 31012:25 2.1 4.6 5.1 280 310 32012:35 2.1 4.6 5.1 280 310 32012:45 2.6 5.1 5.7 290 310 310    (a)     (b) Fig. 1. – (a) Wind standard deviations (σi); (b) Lagrangian decorrelation time scales (τLi) fromCopenhagen experiment (12:33h-12:53h, 19/10/1978).ESTIMATION OF EMISSION RATE FROM EXPERIMENTAL DATA 181logarithmic profile, namely(5) U(z) =u∗κln(zz0),where u∗ is the local friction velocity, κ = 0.4 is the von Karman constant, and z0 theroughness length. Equation (5) is used to compute u∗ at the three levels. Wind speed forall levels in the domain was also computed from eq. (5). The wind direction was linearlyinterpolated from experimental values. The mixing height zi, Obukhov length L, andwind standard deviations were only given as values averaged at each 1 hour and, therefore,were kept constant during the simulation: zi = 1120m, L = −71m, σv = 0.85ms−1 andσw = 0.68ms−1. In this experiment the emission rate was constant and equal to 3.2 g/s.The horizontal LAMBDA domain was determined according to maximum sampler dis-tance and the vertical domain was set equal to zi. The time step was maintained constant(Δt = 1 s). One hundred particles were released at each time step during 2880 time steps.The Gram-Charlier probability density function (PDF) truncated to the third order waschosen to represent the PDF of the vertical wind velocity [16,17].Since the meteorological information is available at every 10min, in the simulation allthe meteorological parameters, except zi and L, were allowed to vary each 10min. Thesimulation starts 28min before the second sampling period, in order to better mimethe experimental conditions (thus allowing the particles previously emitted to “fill” thecomputation domain as it occurred with the tracer in the experiment). To better estimatethe ground level concentration (g.l.c.) in this non-stationary case, g.l.c. is computed 10times (each 2min) during each 20min period, and then the average is computed.The model performance is shown in table II and figs. 2a and b. Table II presentsthe comparison between observed (CExp) and predicted (CMod) values of ground-levelconcentration in the samplers. Figure 2a shows the scatter diagram between observed andpredicted concentrations. Table III shows the result of a standard statistical analysismade with the data presented in table II. The statistical indices are the following:Normalized Mean Square Error, NMSE = (CExp − CMod)2/CExpCMod;Fractional Bias, FB = (CExp − CMod)/(0.5(CExp + CMod));Fractional Standard Deviation, FS = 2(σExp − σMod)/(σExp + σMod);Correlation Coefficient, COR = (CExp − CExp)(CMod − CMod)/σExpσMod;Factor 2, FA2 = 0.5 ≤ CExp/CMod ≤ 2,where σ are concentration standard deviations.From fig. 2b, the observed and computed g.l.c. at each sampler of the three arcsis shown, suggesting that the model was able to capture the general behaviour of thediffusion experiment. However, the simulated plume is wider than the observed one.Thus the maxima at each arc are lower than the observed ones, whereas at the furthersamplers the model computes more concentration than measured.5. – Results for the emission rate estimationIn the estimation of the emission rates it was assumed that the pollutant emission rateis changing with time, i.e. Q = Q(t). Actually, in the Copenhagen experiment, the emis-sion rate is constant and equal to 3.2 gs−1. The pollutant source term is estimated by182 D. R. ROBERTI, D. ANFOSSI, H. F. DE CAMPOS VELHO and G. A. DEGRAZIATable II. – Ground-level concentration (CExp) measured during the Copenhagen experiment(12:33h-12:53h, 19/10/1978) and predicted concentrations (CMod) by LAMBDA model usingturbulence parameterisation given by Degrazia et al. [15].Number X (m) Y (m) CExp (μg/s) CMod (μg/s)1 1398 −1312 0 0.7212 1404 −1214 0 1.3503 1492 −1131 0 2.2714 1516 −1044 0.614 3.2665 1582 −964 1.816 4.4786 1592 −884 5.455 5.1057 1602 −798 7.016 5.7278 1703 −767 6.770 5.5829 1766 −681 5.472 5.07710 1800 −593 3.806 4.06711 1877 −485 1.114 2.84212 1921 −405 0.919 1.96913 2067 −371 0.077 1.90214 2061 −284 0 1.00115 2055 −180 0 0.46316 2818 −2134 0 1.40017 2920 −1987 0.107 1.57918 3002 −1830 0.840 2.63519 3075 −1704 1.478 2.69120 3204 −1629 3.133 2.58021 3380 −1367 2.563 2.24922 3448 −1231 2.225 2.03323 3518 −1093 0.538 1.28224 3558 −919 0 1.05025 3729 −787 0 0.36126 3837 −550 0 0.12127 4027 −3616 0 0.30628 4283 −3447 0.018 0.78829 4390 −3277 0.021 1.08730 4459 −3010 0.085 1.117ESTIMATION OF EMISSION RATE FROM EXPERIMENTAL DATA 183Table II. – Continued.31 4572 −2795 0.800 1.55232 4668 −2514 1.502 1.23633 4824 −2260 2.035 1.29734 5029 −2108 1.112 1.06335 5286 −1939 0.434 0.77836 5378 −1570 0.053 0.12437 5395 −1399 0 0.04938 5375 −1139 0 0.10939 5323 −913 0 0.047minimizing the following functional:(6) F (λ,Q) =40∑i=1[CModi (Q)− CExpi]2+ λ5∑k=1(Qk+1 − 2Qk + Qk−1)2,CMod being the concentration computed from the Lagrangian model, and CExpi theconcentration measured by the i-th sensor. The time period for our estimation was50 minutes. The source term Q(t) is assumed to be constant into the period t ∈[t, t + 10min]. Therefore, the unknown source term can be represented by the vector:Q = [Q1, Q2, Q3, Q4, Q5]T, where Qi = Q(t0 + iΔt) with Δt = 10min.In order to complete the inverse analysis, the regularization parameter λ must bedetermined. This parameter was found by using the procedure proposed by Hansen [18]:some value on the corner of the L-curve, relating “Regularization operator term × SquareDifference” of the concentration— the second and first terms of the right-hand side ofeq. (6), respectively. Figure 3 displays the L-curve for different regularization parameters,using the Tikhonov-2 regularization scheme. From the plot, it can be seen that λ = 0.6is a good choice.The estimation results are presented in table IV and fig. 4, using the second-orderTikhonov regularization scheme. It can be realized that the regularization is the very im-portant issue in the inversion procedure. The estimation without regularization producedwrong estimation, indicating that the least-square scheme itself is not able to correctlyaddress the model to find good answers. However, a bad choice of the regularizationparameter (for example, λ = 40) will not produce a good inversion either. Therefore, forλ → 0 and λ →∝ implies a wrong inverse procedure. The appropriate inverse regular-Table III. – Statistical indices of the LAMBDA model performance for the Copenhagen experi-ment simulated.NMSE COR FA2 FB FS0.51 0.87 0.36 −0.38 0.19184 D. R. ROBERTI, D. ANFOSSI, H. F. DE CAMPOS VELHO and G. A. DEGRAZIA       (a)       (b) Fig. 2. – (a) Scatter diagram between observed and predicted samplers ground-level concentra-tions; (b) experimental and predicted samplers ground-level concentrations (the samplers at thethree arcs are sequentially plotted) from Copenhagen experiment (12:33h-12:53h, 19/10/1978).Table IV. – Emission rate estimation.Time Qtrue (gs−1)Qestimated (gs−1)(h:min) λ = 0 λ = 0.6 λ = 4012:05 3.20 3.57 3.185 0.65512:15 3.20 1.79 3.185 0.59512:25 3.20 4.82 3.187 0.50212:35 3.20 2.74 3.188 0.60712:45 3.20 0.25 3.188 0.704ESTIMATION OF EMISSION RATE FROM EXPERIMENTAL DATA 185Fig. 3. – L-curve for reconstruction of emission rate.ized solution is only obtained by using a correct value for the regularization parameter,indicating a good balance between the data adhesion (the square difference term in ex-pression (6)) and the smoothness (regularization) of the estimated function. Table IVshows that the inversion with regularization parameter λ = 0.6 presents good results.6. – ConclusionsAn ineluctable feature with dealing with experimental data is the noise associatedof such data. Inverse problems belong to the class of ill-posed problems, where they Fig. 4. – Pollutant source Q(t) estimation: effect of the regularization parameter. Without regu-larization (λ = 0): oscillatory behavior; strong regularization (λ = 40), only regularization oper-ator is focused in the minimization process; appropriate value (λ = 0.6), where Q(t) ≈ 3.2 gs−1.186 D. R. ROBERTI, D. ANFOSSI, H. F. DE CAMPOS VELHO and G. A. DEGRAZIAare unstable in the presence of noise (small variation in the input data implies a widevariation in the output data). The regularized inverse solution is a strategy to give agood answer [5]. The methodology used for estimating the emission rate was effectiveto produce correct reconstructions of the emission rate. In this paper, these results wereobtained using the second-order Tikhonov regularization.In this paper it was shown that the inverse procedure yielded good results, even dealingwith data from the real field (the Copenhagen tracer experiment), where some level ofnoise in the data is expected. However, one important issue is the determination of theregularization parameter. The L-curve scheme [18] permitted to find out an appropriatedvalue for regularization operator.∗ ∗ ∗The authors would like to thank to the Brazilians financial support agencies. The firstauthor acknowledges CNPq-Brasil (Proc. 20.0046/04-7) for the sandwich doctorate fel-lowship at ISAC-CNR—Turin—Italy.REFERENCES[1] Seibert P., in Inverse Methods in Global Biogeochemical Cycles, edited by KasibhatlaP. et al. (American Geophysical Union, Washington) 2000, pp. 147-154.[2] Ferrero E., Anfossi D., Brusasca G. and Tinarelli G., Int. J. Environ. Pollut., 5(1995) 360.[3] Carvalho J. C., Degrazia G. A., Anfossi D., Campos C. R. J de, Roberti D. R.and Kerr A., Meteorol. Z., 11 (2002) 89.[4] Thomson D. J., J. Fluid Mech., 180 (1987) 529.[5] Tikhonov A. N. and Arsenin V. I., Solutions of Ill-posed Problems (John Wiley & Sons)1977.[6] Ramos F. M. and Giovannini A., Int. J. Heat Mass Transfer, 38 (1995) 101.[7] Ramos F. M. and Campos Velho H. F., in 2nd International Conference on InverseProblems on Engineering, Le Croisic, France, 1996, pp. 199-206.[8] Stephany S., Ramos F. M., Campos Velho H. F. de and Mobley C. D., Comput.Modelling Simul. Eng., 3 (1998) 161.[9] Muniz W. B., Campos Velho H. F. de and Ramos F. M., J. Comput. Appl. Math.,110 (1999) 153.[10] Campos Velho H. F. de, Moraes M. R. de, Ramos F. M., Degrazia G. A. andAnfossi D., Nuovo Cimento C, 23 (2000) 65.[11] Muniz W. B., Ramos F. M. and Campos Velho H. F., Comput. Math. Appl., 40 (2000)1071.[12] Roberti D. R., Campos Velho H. F. and Degrazia G. A., Inverse Probl. Sci. Eng.,12 (2004) 329.[13] Gryning S. E. and Lyck E., J. Clim. Appl. Meteorol., 23 (1984) 651.[14] Gryning S. E. and Lyck E., The Copenhagen tracer experiments: reporting ofmeasurements. Available on line at: http://www.risoe.dk/rispubl/VEA/ris-r-1054.htm.[15] Degrazia G. A., Anfossi D., Carvalho J. C., Mangia C., Tirabassi T. and CamposVelho H. F., Atmos. Environ., 34 (2000) 3575.[16] Ferrero E. and Anfossi D., in Air Pollution Modelling and its Applications XII, editedby Gryning S. E. and Chaumerliac N. (Plenum Press, New York) 1998, pp. 673-680.[17] Ferrero E. and Anfossi D., Int. J. Environ. Pollut., 9 (1998) 384.[18] Hansen P. C., SIAM Rev., 34 (1992) 561.

Estimation of emission rate from experimental data

http://eprints.bice.rm.cnr.it/16230/1/ncc9229.pdf

Estimation of emission rate from experimental data

Abstract

Similar works

Full text

Available Versions

Scientific Open-access Literature Archive and Repository