Our study constitutes an attempt to evaluate the air quality in semiarid region, Tiaret city, from Algeria. The first aim of this work (case I) was to assess the most frequent atmospheric heavy metal, especially (Pb, Zn, Cu), accumulated in leaves of two urban plant species, plane and cypress, in the context of their usefulness in urban environment as bioindicator or in phytoremediation. Results show that the ratio fresh weight/dry weight (FM/DM) of trees leaves, which is an indicator of air quality in the area of investigation, is negatively correlated with the concentration of trace elements. It is extremely low in the contaminated sites, reflecting the effect of urban pollutants on trees vigor and their physiological activity. Furthermore, we found that the plane tree leaves accumulate more zinc than lead against very low concentrations of copper. Nevertheless, among the cypress, the degree of accumulation of lead is greater than that of zinc, furthermore, very low concentrations or absent of copper were observed in both species. In addition, these results show the interest of the use of leaves of woody plants as bioindicators of air contamination and the impact of

In order to identify the main source of lead and zinc accumulation in trees leaves then, in the intend of estimate concentrations of these pollutants and mapping their spatial repartition in the city (case II), lead and zinc was measured in 48-lichen *Xanthoria parietina* samples, transplanted for while into different sites located near major highways and in urban areas of Tiaret city. Automatic classification of data has shown three classes of lead-polluted sites and four classes of zinc-polluted sites. These results established that the degree of lead pollution including traffic and roads infrastructures is causing a major source of heavy metals in this city. All the observations of pollution have been mapped to reveal spatial pollution originated from road traffic and spot on the vulnerable sites. Hot points sites were locating on land with relatively steep slope and high traffic, so, this forces the vehicle to release a lot of smoke. Sites with low traffic are also open enough to promote the dispersion of atmospheric deposition. These maps could answer questions about the air pollution problem in Tiaret city. These biological methods contribute to a health risk assessment

because they assist to identify areas potentially exposed to urban air pollution.

These maps could answer questions about the air pollution problem in Tiaret city. These biological methods contribute to a health risk assessment because they assist to identify areas potentially exposed to air pollution. Thus, urban trees can form networks of plant bioindicator of air pollution; the use appears to be particularly simple, economical and efficient method. The second part of the study has demonstrated the perfect suitability of lichens bioaccumulators for mapping different elementary deposits and location of sources

these studies are mosses and lichens through, inter alia, their lack of root system.

accumulate substantially lead better than the cypress leaves.

**5. Conclusion**

pollution on urban plants.

On the same trees (the same sites S1 and S2), lichens *Xanthoria pareitina* were transplanted following the steps described in previous paragraph. Results of this comparison are illustrated in Figure 8.

Fig. 8. Comparison between contents of lead and zinc accumulated by trees (Plane and Cypress) to those accumulated by the lichen Xanthoria in Tiaret City.

We notice very clearly that lead storage capacity of Xanthoria is very high compared to that of the trees, the values are 1.04 ± 0.60 µg/g in Xanthoria against 2.86 ± 0.17 x 10-2 µg/g and 3.39±0.01 x 10-2 µg/g respectively in the plane and cypress.

By the same way, the Xanthoria accumulates 460.38±43.85 µg/g against 380.45±91.53 x 10-2 µg/g and 93.20±20.96 x 10-2 µg/g respectively in plane and cypress leaves (Maatoug, 2010).

Lichens absorb indiscriminately by their pseudocyphelles all the nutrients as toxic substances. Absorption capacity of metals by lichens is directly related to their morphological and anatomical structures. Lichens accumulate more metals when their surfaces increase (Garty et al., 1996).

Sharing some biological traits (Asta et al., 2003), lichens have a relatively quick response to the deteriorating of air quality and are extremely sensitive to other types of environmental changes such as climate change and eutrophication (Galun, 1988), (Garrec and Van Haluwyn, 2002).

In trees, the emissions of lead and zinc originated from road traffic are mainly in the form of fine particles that are then collected by leaves surfaces. It is an assimilation of external metal pollutants. The significant difference in concentrations of lead and zinc accumulated between the plane and cypress is certainly due to the nature of leaves surfaces (mainly on epicuticular waxes) of each tree, the sycamore leaves, which are smooth and wide, accumulate substantially lead better than the cypress leaves.

Generally, optimal use of plants as means of investigation of air pollution concerns cases where the route of contamination by the soil is either nil, negligible or known. The transfer of a trace element from the soil to a harvest organ depends on both parameters related to soil and plant-specific factors. It results in some sort of meeting between ground supply and plants demand that are not completely independent. Therefore, the plants most used for these studies are mosses and lichens through, inter alia, their lack of root system.

#### **5. Conclusion**

240 Air Pollution – Monitoring, Modelling and Health

On the same trees (the same sites S1 and S2), lichens *Xanthoria pareitina* were transplanted following the steps described in previous paragraph. Results of this comparison are

> Pb (ug/g) Outliers Extremes Zn (ug/g) Outliers Extremes

Plane Cypress Xanthoria Species

**x 10-2**

Fig. 8. Comparison between contents of lead and zinc accumulated by trees (Plane and

We notice very clearly that lead storage capacity of Xanthoria is very high compared to that of the trees, the values are 1.04 ± 0.60 µg/g in Xanthoria against 2.86 ± 0.17 x 10-2 µg/g and

By the same way, the Xanthoria accumulates 460.38±43.85 µg/g against 380.45±91.53 x 10-2 µg/g and 93.20±20.96 x 10-2 µg/g respectively in plane and cypress leaves (Maatoug, 2010). Lichens absorb indiscriminately by their pseudocyphelles all the nutrients as toxic substances. Absorption capacity of metals by lichens is directly related to their morphological and anatomical structures. Lichens accumulate more metals when their

Sharing some biological traits (Asta et al., 2003), lichens have a relatively quick response to the deteriorating of air quality and are extremely sensitive to other types of environmental changes such as climate change and eutrophication (Galun, 1988), (Garrec and Van

In trees, the emissions of lead and zinc originated from road traffic are mainly in the form of fine particles that are then collected by leaves surfaces. It is an assimilation of external metal pollutants. The significant difference in concentrations of lead and zinc accumulated between the plane and cypress is certainly due to the nature of leaves surfaces (mainly on

**x 10-2**

illustrated in Figure 8.

0

surfaces increase (Garty et al., 1996).

Haluwyn, 2002).

200

400

**x 10-2**

**x 10-2**

Cypress) to those accumulated by the lichen Xanthoria in Tiaret City.

3.39±0.01 x 10-2 µg/g respectively in the plane and cypress.

600

800

1000

Concentrations

1200

1400

1600

1800

Our study constitutes an attempt to evaluate the air quality in semiarid region, Tiaret city, from Algeria. The first aim of this work (case I) was to assess the most frequent atmospheric heavy metal, especially (Pb, Zn, Cu), accumulated in leaves of two urban plant species, plane and cypress, in the context of their usefulness in urban environment as bioindicator or in phytoremediation. Results show that the ratio fresh weight/dry weight (FM/DM) of trees leaves, which is an indicator of air quality in the area of investigation, is negatively correlated with the concentration of trace elements. It is extremely low in the contaminated sites, reflecting the effect of urban pollutants on trees vigor and their physiological activity. Furthermore, we found that the plane tree leaves accumulate more zinc than lead against very low concentrations of copper. Nevertheless, among the cypress, the degree of accumulation of lead is greater than that of zinc, furthermore, very low concentrations or absent of copper were observed in both species. In addition, these results show the interest of the use of leaves of woody plants as bioindicators of air contamination and the impact of pollution on urban plants.

In order to identify the main source of lead and zinc accumulation in trees leaves then, in the intend of estimate concentrations of these pollutants and mapping their spatial repartition in the city (case II), lead and zinc was measured in 48-lichen *Xanthoria parietina* samples, transplanted for while into different sites located near major highways and in urban areas of Tiaret city. Automatic classification of data has shown three classes of lead-polluted sites and four classes of zinc-polluted sites. These results established that the degree of lead pollution including traffic and roads infrastructures is causing a major source of heavy metals in this city. All the observations of pollution have been mapped to reveal spatial pollution originated from road traffic and spot on the vulnerable sites. Hot points sites were locating on land with relatively steep slope and high traffic, so, this forces the vehicle to release a lot of smoke. Sites with low traffic are also open enough to promote the dispersion of atmospheric deposition. These maps could answer questions about the air pollution problem in Tiaret city. These biological methods contribute to a health risk assessment because they assist to identify areas potentially exposed to urban air pollution.

These maps could answer questions about the air pollution problem in Tiaret city. These biological methods contribute to a health risk assessment because they assist to identify areas potentially exposed to air pollution. Thus, urban trees can form networks of plant bioindicator of air pollution; the use appears to be particularly simple, economical and efficient method. The second part of the study has demonstrated the perfect suitability of lichens bioaccumulators for mapping different elementary deposits and location of sources

Bio-Monitoring of Air Quality Using Leaves of Tree and Lichens in Urban Environments 243

[8] Brodo .I M. Transplant experiments with corticolous lichens using a new technique.

[9] Chocat, B.,. Eurydice. Encyclopédie de l'hydrologie urbaine et de l'assainissement. Paris :

[10] Colandini, V., Effets des structures réservoirs à revêtement poreux sur les eaux de

[11] Dahmani-Muller H., Phytoréhabilitation des sols pollués par des éléments métalliques :

[12] Deletraz G., Géographie des risques environnementaux lies aux transports routiers en

[13] Delmas – Gadras C., Influence des conditions physico-chimiques sur la mobilité du

[14] Deruelle .S. Effets de la pollution atmosphérique sur la végétation lichénique dans le

[15] Deruelle .S. L'utilisation des lichens pour la détection de la pollution par le plomb.

[16] Domergue–Abak. M.F., Etude de synthèse sur les causes de dépérissement de la

[18] Fernandez-cornudet,C., Devenir du Zn, Pb et Cd issus de retombées atmosphériques

[20] Garrec j.p. et Van Haluwyn c., Biosurveillance végétale de la qualité de l'air. Concepts, méthodes et applications. Editions Tec et Doc Lavoisier, Paris (France) 2002,118 p. [21] Garty, J., Kauppi, M., Kauppi, A., Accumulation of airborne elements from vehicles in

[22] Hébrard-Labit. C et Meffray.L, Comparaison de méthodes d'analyse des éléments traces

agglomérations. LEN n° 9515. Bron : INRETS (France) 1995, 125 p.

de l'université de Pau et des Pays de l'Adour (France) 2000, 191 p.

(Phaseolus vulgaris L.). Environ. Toxicol. Chem.. (5),659-65.

[19] Galun, M. Handbook of lichenology. Springer, 1988; 181p.

ruissellement pluviales : qualité des eaux et devenir des métaux lourds. Thèse de doctorat. Pau. Université de Paut des pays de l'Adour (France) 1997, 161 p. +

facteurs et mécanismes de prélèvement dans les sols et d'accumulation par les espèces métalliques. Thèse de l'Ecole Nationale de Génie rural, des eaux et des

montagne. Incidences des émissions d'oxydes d'azote en vallées d'Aspe et de Biriatou (Pyrénées). Thèse de Doctorat en Géographie - Aménagement. Université de Pau et des pays de L'Adour. Institut de Recherche sur les Sociétés et

plomb et du zinc dans un sol et un sédiment en domaine routier. Thèse de docteur

bassin Parisien. Convention de recherche n° 79-15, Ministère de l'Environnement et

végétation en milieu urbain, et notamment les arbres d'alignement. Rapport du Ministère de l'Urbanisme et du Logement (France) 1981, N° 81/47142/00, 52p. [17] Edwards N.T., 1986. Uptake, translocation and metabolism of anthracene in bush bean

dans les sols, à différentes échelles d'étude. Influence de l'usage des sols sur la distribution et la mobilité des métaux. Thèse de doctorat de l'INA-PG (France)

transplanted lichens in urban sites. Journal of Environmental Quality 1996 (25) 265–

métalliques (ETM) et des hydrocarbures aromatiques polycycliques (HAP) sur les sols et les végétaux. Guide technique. CETE Nord Picardie (France) 2004, 121 p. [23] Joumard R., Lamure C., Lambert J., Politiques de transport et qualité de l'air dans les

Ecology 1961 (42): 838-841.

annexes

Tec et Doc Lavoisier. 1997, 1136 p.

Forêts ENGREF (France) 2000, 151 p.

l'Aménagement (France) 2002, 564 p.

du Cadre de Vie 1981 : 91-112.

Bull.Eco 1984 :1-6.

2006, 232p.

272.

of metals. Lichens and urban trees can form networks of plant bioindicators of air pollution, the use appears to be also very simple, flexible, economical and efficient to build good repartition in space and time maps of pollution. The identification of pollution within sensitive organisms can also detect the degradation of air quality before it severely affects the biota or humans.

Today, the use of motor vehicles has increased substantially in Algeria. However, in recent decades, urban and industrial developments have accelerated the use of motor vehicles. Unfortunately, at present, air pollution generated by traffic is taken into account only through its effects on human health. On the other hand, the reduction of pollutant emissions from road traffic demands renewal of Automobile Park and improvement of combustion engines using less pollutant fuels.

In the end, it appears that pollution related to road traffic has an impact on urban environment. Plant contamination comes mainly from the aerial parts, which are often the beginning of food chains with all the problems that may cause, particularly on human health. This study confirms the need to control the automobile park and to reduce emissions from road traffic by the renewal of the automobile park, improvements in the accuracy of regulating combustion engines and the use less pollutant fuels.

#### **6. References**

of metals. Lichens and urban trees can form networks of plant bioindicators of air pollution, the use appears to be also very simple, flexible, economical and efficient to build good repartition in space and time maps of pollution. The identification of pollution within sensitive organisms can also detect the degradation of air quality before it severely affects

Today, the use of motor vehicles has increased substantially in Algeria. However, in recent decades, urban and industrial developments have accelerated the use of motor vehicles. Unfortunately, at present, air pollution generated by traffic is taken into account only through its effects on human health. On the other hand, the reduction of pollutant emissions from road traffic demands renewal of Automobile Park and improvement of combustion

In the end, it appears that pollution related to road traffic has an impact on urban environment. Plant contamination comes mainly from the aerial parts, which are often the beginning of food chains with all the problems that may cause, particularly on human health. This study confirms the need to control the automobile park and to reduce emissions from road traffic by the renewal of the automobile park, improvements in the accuracy of

[1] ADEME,. Agence de l'Environnement et de la Maîtrise de l'Energie. Dérogations

[2] Alfani, A. et Baldantoni. D. Temporal and spatial variation of C, N, S and trace element

[3] Asta, J., Erhardt, W., Ferreti, M., Fornassier, F., Kirschbaum, U., Nimis, P.L., Purvis,

[4] Baize D., 1997. Teneurs totales en éléments traces métalliques dans les sols (France).

[5] Bazzaz, F.A., Rolfe, G.L., et Carlson, R.W., Effect of cadmium on photosynthesis and

[6] Belouahem D., Détection de la pollution atmosphérique fluorée d'origine industrielle à

[7] Bhatti G.H., Iqbal M.Z, 1988. Investigations into the effect of automobile exhausts on the

relatives à la réglementation sur l'épandage des boues de stations d'épuration. Comment formuler une demande pour les sols à teneurs naturelles élevées en éléments traces métalliques ?. Guide technique. J. Béraud et A. Bispo (Coordinateurs). D. Baize, T. Sterckeman, A. Piquet, H. Ciesielski, J. Béraud et A.

contents in the leaves of Quercus ilex within the urban area of Naples. Environ.

O.W., Pirintsos, S., Scheidegger, C., van Haluwyn, C., Wirth, V. (2003). European guideline for mapping lichen diversity as an indicator of environmental stress. 20p.

transpiration of excised leaves of corn and sunflower. Physiologia Plantarum, 1974

l'aide de certaines espèces végétales bioaccumulatrices dans les régions de Annaba et Taraf. Thèse de Magister. Institut National Agronomique INA (Algérie) 1993.

phenology, periodicity and productivity of some roadside trees, Acta Societatis

regulating combustion engines and the use less pollutant fuels.

the biota or humans.

**6. References**

engines using less pollutant fuels.

Bispo (2005). 120p

(32), 373-377.

165p.

Pollut. 2000 (109): 119-129.

INRA de France. Editions, Paris, 408p.

Botanicorum Poloniae, 57 (3). 395-399.

**11**

*Greece*

**Developing Neural Networks**

Kyriaki Kitikidou and Lazaros Iliadis

*Democritus University of Thrace, Department of Forestry and Management of the Environment and Natural Resources, Orestiada*

**to Investigate Relationships Between**

**Air Quality and Quality of Life Indicators**

Quality of life (QOL) is an integral outcome measure in the management of diseases. It can be used to assess the results of different management methods, in relation to disease complications and in fine-tuning management methods (Koller & Lorenz, 2003). Quantitative analysis of quality of life across countries, and the construction of summary indices for such analyses have been of interest for some time (Slottje et al., 1991). Most early work focused on largely single dimensional analysis based on such indicators as per capita GDP, the literacy rate, and mortality rates. Maasoumi (1998) and others called for a multidimensional quantitative study of welfare and quality of life. The argument is that welfare is made up of several distinct dimensions, which cannot all be monetized, and heterogeneity complications are best accommodated in multidimensional analysis. Hirschberg et al. (1991) and Hirschberg et al. (1998) identified similar indicators, and collected them into distinct clusters which could represent the dimensions worthy of distinct

In this research effort we have considered the role of air quality indicators in the context of economic and welfare life quality indicators, using artificial neural networks (ANN). Therefore in this presentation we have obtained the key variables (life expectancy, healthy life years, infant mortality, Gross Domestic Product (GDP) and GDP growth rate) and developed a neural network model to predict the air quality outcomes (emissions of sulphur and nitrogen oxides). Sustainability and quality of life indicators have been proposed recently by Flynn et al. (2002) and life quality indices have been used to estimate willingness to pay (Pandey & Nathwani, 2004). The innovative part of this research effort lies in the use of a soft computing machine learning approach like the ANN to predict air quality. In this way, we introduce the reader to a technique that allows the comparison of various attributes

**1. Introduction**

treatment in multidimensional frameworks.

that impact the quality of life in a meaningful way.

### **Developing Neural Networks to Investigate Relationships Between Air Quality and Quality of Life Indicators**

Kyriaki Kitikidou and Lazaros Iliadis

*Democritus University of Thrace, Department of Forestry and Management of the Environment and Natural Resources, Orestiada Greece*

#### **1. Introduction**

244 Air Pollution – Monitoring, Modelling and Health

[24] Kupper H, Lombi E, Zhao F J, Wieshammer G, McGrath S P., Cellular

[25] Leblond. S. Etude pluridisciplinaire du transfert des métaux de l'atmosphère vers les

[26] Maatoug M. Détection de la pollution de l'air d'origine routière par certaines espèces

[27] Maatoug M. Cartographie de la pollution atmosphérique par le plomb d'origine

[28] Madany IM, Ali SM, Akhter MS. Assessment of lead in roadside vegetation in Bahrain.

[29] Ouzounidou G, Ciamporova M, Moustakas M, Karataglis S.,. Responses of Maize (Zea-

[30] Ouzounidou G, Moustakas M, Eleftheriou E P.,. Physiological and ultrastructural effects

[31] Peulon V., Le dépérissement des arbres en ville. Edition du STU, Ministère de

[32] Semadi .A et Deruelle .S Détection la pollution plombique à l'aide de transplants

[33] Vazquez M D, Poschenrieder C, Barcelo J., Ultrastructural Effects and Localization of LowCadmium Concentrations in Bean Roots. New Phytol; 1992 (120), 215-226.

Roots. Environmental and Experimental Botany; 1995 (35), 167-176.

de l'Environnement, Université Paris 7 - Denis Diderot 2004. 212p.

pollution atmosphérique 2007 ; 196 : 385-394.

Environment International 1990 (16): 123-126.

Contamination and Toxicology; 1997 (32),154-160.

l'Equipement (France) 1988, 61p.

Octobre-Décembre 1993 : 86-101.

(52), 2291-2300.

compartimentation of nickel in the hyperaccumulators Alyssum lesbiacum, Alyssum bertolonii and Thlaspi goesingense. Journal of Experimental Botany; 2001

mousses (Scleropodium purum (Hedw.) Limpr.) : Suivi sur un site rural (Vouzon, France). Thèse de doctorat en Chimie de la Pollution Atmosphérique et Physique

végétales bioaccumulatrices de quelques métaux lourds (Pb, Zn, Cu). Revue

routière à l'aide de transplantation d'un lichen bioaccumulateur *Xanthoria parietina* dans la ville de Tiaret (Algérie). Revue pollution atmosphérique 2010 ; 205 : 93-101.

Mays L) Plants to Copper Stress .1. Growth, Mineral-Content and Ultrastructure of

of cadmium on wheat (Triticum aestivum L) leaves. Archives of Environmental

lichéniques dans la région de Annaba (Algérie). Revue pollution atmosphérique

Quality of life (QOL) is an integral outcome measure in the management of diseases. It can be used to assess the results of different management methods, in relation to disease complications and in fine-tuning management methods (Koller & Lorenz, 2003). Quantitative analysis of quality of life across countries, and the construction of summary indices for such analyses have been of interest for some time (Slottje et al., 1991). Most early work focused on largely single dimensional analysis based on such indicators as per capita GDP, the literacy rate, and mortality rates. Maasoumi (1998) and others called for a multidimensional quantitative study of welfare and quality of life. The argument is that welfare is made up of several distinct dimensions, which cannot all be monetized, and heterogeneity complications are best accommodated in multidimensional analysis. Hirschberg et al. (1991) and Hirschberg et al. (1998) identified similar indicators, and collected them into distinct clusters which could represent the dimensions worthy of distinct treatment in multidimensional frameworks.

In this research effort we have considered the role of air quality indicators in the context of economic and welfare life quality indicators, using artificial neural networks (ANN). Therefore in this presentation we have obtained the key variables (life expectancy, healthy life years, infant mortality, Gross Domestic Product (GDP) and GDP growth rate) and developed a neural network model to predict the air quality outcomes (emissions of sulphur and nitrogen oxides). Sustainability and quality of life indicators have been proposed recently by Flynn et al. (2002) and life quality indices have been used to estimate willingness to pay (Pandey & Nathwani, 2004). The innovative part of this research effort lies in the use of a soft computing machine learning approach like the ANN to predict air quality. In this way, we introduce the reader to a technique that allows the comparison of various attributes that impact the quality of life in a meaningful way.

Developing Neural Networks to Investigate

Emissions

of sulphur oxides (million tones of SO2 equivalent)

Statistics

Std.

movements will not inflate the growth rate.

Emissions

of nitrogen oxides (million tones of NO2 equivalent)

\*Number of observations (countries) for each variable. \*\*Number of countries that didn't had available data.

Table 1. Descriptive statistics for all variables used in the analysis.

continuous, single valued with first derivative existence (Picton, 2000).

(Eurostat, 2010). Descriptive statistics for all variables are given in **Table 1**.

Infant Mortality

Valid N\* 34 34 34 33 33 33 27 Missing\*\* 0 0 0 1 1 1 7 Mean 0.503 0.372 5.721 95.921 4.206 77.535 60.448

Deviation 0.648 0.482 4.227 46.620 2.521 3.244 5.443 Min 0.00 0.00 2.30 28.50 0.70 70.94 50.10 Max 2.37 1.63 23.60 254.50 10.60 81.54 69.30

For the performance of the analyses, multi-layer perceptron (MLP) and radial-basis function (RBF) network models were developed under the SPSS v.19 statistical package (IBM, 2010). We specified that the relative number of cases assigned to the training:testing:holdout samples should be 6:2:1. This assigned 2/3 of the cases to training, 2/9 to testing, and 1/9 to holdout. For the MLP network we employed the back propagation (BP) optimization algorithm. As it is well known in BP the weighted sum of inputs and bias term are passed to the activation level through the transfer function to produce the output (Bishop, 1995; Fine, 1999; Haykin, 1998; Ripley, 1996). The sigmoid transfer function was employed (Callan, 1999; Kecman, 2001), due to the fact that the algorithm requires a response function with a

Before using the input data records to the ANN a normalization process took place so that the values with wide range do not prevail over the rest. The autoscaling approach was applied. This method outputs a zero mean and unit variance of any descriptor variable (Dogra, Shaillay, 2010). Thus, each feature's values were normalized based on the following equation:

Relationships Between Air Quality and Quality of Life Indicators 247

Data were extracted for 34 European countries, for the year 2005, from the Eurostat database

GDP (Purchasing Power Standards, PPS)

GDP Growth Rate

Life

Expectancy At Birth (years)

Healthy Life Years

 GDP growth rate: The calculation of the annual growth rate of GDP volume is intended to allow comparisons of the dynamics of economic development both over time and between economies of different sizes. For measuring the growth rate of GDP in terms of volumes, the GDP at current prices are valued in the prices of the previous year and the thus computed volume changes are imposed on the level of a reference year; this is called a chain-linked series. Accordingly, price

#### **2. Materials and methods**

It is well known that the quality of the air in a locale influences the health of the population and ultimately affects other dimensions of that population's welfare and its economy. As a simple example, in cities where pollution levels rise significantly in the summer, worker absenteeism rates rise commensurately and productivity is adversely impacted. Other dimensions of the economy are influenced on "high pollution days" as well. For example, when outdoor leisure activity is restricted this may have serious consequences for the service sector of the economy (Bresnahan et al., 1997). In this chapter, we have introduced two measures of environmental quality or air quality as quality of life factors. A feature of these indices is the fact that these types of pollution are created by some of the very activities that define economic development. The two factors under investigation here are sulfur oxides (SOx) and nitrogen oxides (NOx) (million tones of SO2 and NO2 equivalent, respectively). Sulphur oxides, including sulphur dioxide and sulphur trioxide, are reported as sulphur dioxide equivalent, while nitrogen oxides, including nitric oxide and nitrogen dioxide, are reported as nitrogen dioxide equivalent. They are both produced as byproducts of fuel consumption as in case of the generation of electricity. Vehicle engines also produce a large proportion of NOx. SOx is primarily produced when high sulphur content coal is burned which is usually in large-scale industrial processes and power generation. Thus, the ratio of these emissions to the population is an indication of pollution control.

The following attributes of QOL have been used:

	- Gross Domestic Product (GDP) per capita: GDP is a measure of the economic activity, defined as the value of all goods and services produced less the value of any goods or services used in their creation. These amounts are expressed in PPS (Purchasing Power Standards), i.e. a common currency that eliminates the differences in price levels between countries allowing meaningful volume comparisons of GDP between countries.

It is well known that the quality of the air in a locale influences the health of the population and ultimately affects other dimensions of that population's welfare and its economy. As a simple example, in cities where pollution levels rise significantly in the summer, worker absenteeism rates rise commensurately and productivity is adversely impacted. Other dimensions of the economy are influenced on "high pollution days" as well. For example, when outdoor leisure activity is restricted this may have serious consequences for the service sector of the economy (Bresnahan et al., 1997). In this chapter, we have introduced two measures of environmental quality or air quality as quality of life factors. A feature of these indices is the fact that these types of pollution are created by some of the very activities that define economic development. The two factors under investigation here are sulfur oxides (SOx) and nitrogen oxides (NOx) (million tones of SO2 and NO2 equivalent, respectively). Sulphur oxides, including sulphur dioxide and sulphur trioxide, are reported as sulphur dioxide equivalent, while nitrogen oxides, including nitric oxide and nitrogen dioxide, are reported as nitrogen dioxide equivalent. They are both produced as byproducts of fuel consumption as in case of the generation of electricity. Vehicle engines also produce a large proportion of NOx. SOx is primarily produced when high sulphur content coal is burned which is usually in large-scale industrial processes and power generation. Thus, the ratio of these emissions to the

 Life expectancy at birth: The mean number of years that a newborn child can expect to live if subjected throughout his life to the current mortality conditions (age specific

 Healthy life years: The indicator Healthy Life Years (HLY) at birth measures the number of years that a person at birth is still expected to live in a healthy condition. HLY is a health expectancy indicator which combines information on mortality and morbidity. The data required are the age-specific prevalence (proportions) of the population in healthy and unhealthy conditions and age-specific mortality information. A healthy condition is defined by the absence of limitations in functioning/disability. The indicator is also called disability-free life expectancy (DFLE). Life expectancy at birth is defined as the mean number of years still to be lived by a person at birth, if subjected throughout the rest of his or her life to the current mortality conditions

 Infant mortality: The ratio of the number of deaths of children under one year of age during the year to the number of live births in that year. The value is expressed per 1000

 Gross Domestic Product (GDP) per capita: GDP is a measure of the economic activity, defined as the value of all goods and services produced less the value of any goods or services used in their creation. These amounts are expressed in PPS (Purchasing Power Standards), i.e. a common currency that eliminates the differences in price levels between countries allowing meaningful volume

**2. Materials and methods**

population is an indication of pollution control. The following attributes of QOL have been used:

comparisons of GDP between countries.

probabilities of dying).

(WHO, 2010).

live births.

 GDP growth rate: The calculation of the annual growth rate of GDP volume is intended to allow comparisons of the dynamics of economic development both over time and between economies of different sizes. For measuring the growth rate of GDP in terms of volumes, the GDP at current prices are valued in the prices of the previous year and the thus computed volume changes are imposed on the level of a reference year; this is called a chain-linked series. Accordingly, price movements will not inflate the growth rate.

Data were extracted for 34 European countries, for the year 2005, from the Eurostat database (Eurostat, 2010). Descriptive statistics for all variables are given in **Table 1**.

\*Number of observations (countries) for each variable.

\*\*Number of countries that didn't had available data.

Table 1. Descriptive statistics for all variables used in the analysis.

For the performance of the analyses, multi-layer perceptron (MLP) and radial-basis function (RBF) network models were developed under the SPSS v.19 statistical package (IBM, 2010). We specified that the relative number of cases assigned to the training:testing:holdout samples should be 6:2:1. This assigned 2/3 of the cases to training, 2/9 to testing, and 1/9 to holdout. For the MLP network we employed the back propagation (BP) optimization algorithm. As it is well known in BP the weighted sum of inputs and bias term are passed to the activation level through the transfer function to produce the output (Bishop, 1995; Fine, 1999; Haykin, 1998; Ripley, 1996). The sigmoid transfer function was employed (Callan, 1999; Kecman, 2001), due to the fact that the algorithm requires a response function with a continuous, single valued with first derivative existence (Picton, 2000).

Before using the input data records to the ANN a normalization process took place so that the values with wide range do not prevail over the rest. The autoscaling approach was applied. This method outputs a zero mean and unit variance of any descriptor variable (Dogra, Shaillay, 2010). Thus, each feature's values were normalized based on the following equation:

Developing Neural Networks to Investigate

Fig. 2. Radial-basis function network structure.

Relationships Between Air Quality and Quality of Life Indicators 249

#### Zi=(Xi-μi)/σ<sup>i</sup>

where *Xi* was the *i*th parameter, *Zi* was the scaled variable following a normal distribution and *σi*, *μi* were the standard deviation and the mean value of the *i*th parameter.

These networks were trained in an iterative process. A single hidden sub layer architecture was followed in order to reduce the complexity of the network, and increase the computational efficiency (Haykin, 1998). Two units were chosen in the hidden layer. The schematic representation of the neural network is given in **Fig. 1**.

Fig. 1. Multi-layer perceptron network structure.

As regards the RBF network (Bishop, 1995; Haykin, 1998; Ripley, 1996; Tao, 1993; Uykan et al., 2000), the architecture that was developed included nine neurons in the hidden layer. The transfer functions (hidden layer activation functions and output function) determine the output by depicting the result of the distance function (Bors & Pitas, 2001; Iliadis, 2007). The schematic representation of the neural network with transfer functions is given in **Fig. 2.**

Zi=(Xi-μi)/σ<sup>i</sup> where *Xi* was the *i*th parameter, *Zi* was the scaled variable following a normal distribution

These networks were trained in an iterative process. A single hidden sub layer architecture was followed in order to reduce the complexity of the network, and increase the computational efficiency (Haykin, 1998). Two units were chosen in the hidden layer. The

As regards the RBF network (Bishop, 1995; Haykin, 1998; Ripley, 1996; Tao, 1993; Uykan et al., 2000), the architecture that was developed included nine neurons in the hidden layer. The transfer functions (hidden layer activation functions and output function) determine the output by depicting the result of the distance function (Bors & Pitas, 2001; Iliadis, 2007). The schematic representation of the neural network with transfer functions is given in **Fig. 2.**

and *σi*, *μi* were the standard deviation and the mean value of the *i*th parameter.

schematic representation of the neural network is given in **Fig. 1**.

Fig. 1. Multi-layer perceptron network structure.

Fig. 2. Radial-basis function network structure.

Developing Neural Networks to Investigate

Table 3. MLP Parameter Estimates.

NO 0.005 1.049NO 2 2 *error*

( 2 2 SO SO *a b error*

values of MLP.

Predictor Predicted

Input Layer (Bias) -0.119 -0.537

Life Expectancy

At Birth

, NO NO 2 2 *a b error*

errors of the ANN are overestimations of the target values.

Relationships Between Air Quality and Quality of Life Indicators 251

1.646 1.226

H(1:1) -0.518 0.116 H(1:2) 1.396 1.395

Infant Mortality -0.805 0.752 GDP 1.033 -3.377 GDP Growth Rate 0.318 -3.767

Healthy Life Years 0.567 0.358 Hidden Layer 1 (Bias) -0.635 -0.877

Linear regression between observed and predicted values

reasonably good job of predicting emissions of sulphur and nitrogen oxides. Ideally, linear regression parameters *a* and *b* should have values 0 and 1, respectively, while values of the observed-by-predicted chart should lie roughly along a straight line. Linear regression gave

Fig. 3. Linear regression of observed values for emissions of sulphur oxides by predicted

 (**Fig. 4**), respectively. There appears to be more error in the predictions of emissions of sulphur oxides than in emissions of nitrogen oxides, something that we also pointed out in Table 2. **Figs 3 and 4** actually seem to suggest that the largest

results for the two output variables 2 2 SO 0.114 0.918SO *error*

Hidden Layer 1 Output Layer H(1:1) H(1:2) SO2 NO2

) showed that the MLP network does a

(**Fig. 3**) and

#### **3. Results – Discussion**

From the MLP analysis, 19 cases (70.4%) were assigned to the training sample, 2 (7.4%) to the testing sample, and 6 (22.2%) to the holdout sample. The choice of the records was done in a random manner. The whole effort targeted in the development of an ANN that would have the ability to generalize as much as possible. The seven data records which were excluded from the analysis were countries that did not had available data on Healthy Life Years. Two units were chosen in the hidden layer.

**Table 2** displays information about the results of training and applying the MLP network to the holdout sample. Sum-of-squares error is displayed because the output layer has scaledependent variables. This is the error function that the network tries to minimize during training. One consecutive step with no decrease in error was used as stopping rule. The relative error for each scale-dependent variable is the ratio of the sum-of-squares error for the dependent variable to the sum-of-squares error for the "null" model, in which the mean value of the dependent variable is used as the predicted value for each case. There appears to be more error in the predictions of emissions of sulphur oxides than in emissions of nitrogen oxides, in the training and holdout samples.

The average overall relative errors are fairly constant across the training (0.779), testing (0.615), and holdout (0.584) samples, which give us some confidence that the model is not overtrained and that the error in future cases, scored by the network will be close to the error reported in this table

Table 2. MLP Model Summary.

In the following **Table 3** parameter estimates for input and output layer, with their corresponding biases, are given.

Table 3. MLP Parameter Estimates.

250 Air Pollution – Monitoring, Modelling and Health

From the MLP analysis, 19 cases (70.4%) were assigned to the training sample, 2 (7.4%) to the testing sample, and 6 (22.2%) to the holdout sample. The choice of the records was done in a random manner. The whole effort targeted in the development of an ANN that would have the ability to generalize as much as possible. The seven data records which were excluded from the analysis were countries that did not had available data on Healthy Life

**Table 2** displays information about the results of training and applying the MLP network to the holdout sample. Sum-of-squares error is displayed because the output layer has scaledependent variables. This is the error function that the network tries to minimize during training. One consecutive step with no decrease in error was used as stopping rule. The relative error for each scale-dependent variable is the ratio of the sum-of-squares error for the dependent variable to the sum-of-squares error for the "null" model, in which the mean value of the dependent variable is used as the predicted value for each case. There appears to be more error in the predictions of emissions of sulphur oxides than in emissions of

The average overall relative errors are fairly constant across the training (0.779), testing (0.615), and holdout (0.584) samples, which give us some confidence that the model is not overtrained and that the error in future cases, scored by the network will be close to the

> Sum of Squares Error 14.029 Average Overall Relative Error 0.779

Sum of Squares Error 0.009 Average Overall Relative Error 0.615

Average Overall Relative Error 0.584

In the following **Table 3** parameter estimates for input and output layer, with their

Emissions of sulphur oxides (million tones of SO2 equivalent) 0.821

Emissions of nitrogen oxides (million tones of NO2 equivalent) 0.738

Emissions of sulphur oxides (million tones of SO2 equivalent) 0.390

Emissions of nitrogen oxides (million tones of NO2 equivalent) 0.902

Emissions of sulphur oxides (million tones of SO2 equivalent) 0.603

Emissions of nitrogen oxides (million tones of NO2 equivalent) 0.568

**3. Results – Discussion**

error reported in this table

.

Training

Testing

Holdout

Years. Two units were chosen in the hidden layer.

nitrogen oxides, in the training and holdout samples.

Relative Error for Scale

Relative Error for Scale

Relative Error for Scale

Dependents

Dependents

Dependents

Table 2. MLP Model Summary.

corresponding biases, are given.

Linear regression between observed and predicted values ( 2 2 SO SO *a b error* , NO NO 2 2 *a b error* ) showed that the MLP network does a reasonably good job of predicting emissions of sulphur and nitrogen oxides. Ideally, linear regression parameters *a* and *b* should have values 0 and 1, respectively, while values of the observed-by-predicted chart should lie roughly along a straight line. Linear regression gave results for the two output variables 2 2 SO 0.114 0.918SO *error* (**Fig. 3**) and NO 0.005 1.049NO 2 2 *error* (**Fig. 4**), respectively. There appears to be more error in the predictions of emissions of sulphur oxides than in emissions of nitrogen oxides, something that we also pointed out in Table 2. **Figs 3 and 4** actually seem to suggest that the largest errors of the ANN are overestimations of the target values.

Fig. 3. Linear regression of observed values for emissions of sulphur oxides by predicted values of MLP.

Developing Neural Networks to Investigate

reason.

indicators), followed distantly by other predictors.

oxides, in the training and holdout samples.

Relationships Between Air Quality and Quality of Life Indicators 253

that the results are dominated by GDP growth rate and GDP (strictly economical QOL

From the RBF analysis, 19 cases (70.4%) were assigned to the training sample, 1 (3.7%) to the testing sample, and 7 (25.9%) to the holdout sample. The seven data records which were excluded from the MLP analysis were excluded from the RBF analysis also, for the same

**Table 4** displays the corresponding information from the RBF network. There appears to be more error in the predictions of emissions of sulphur oxides than in emissions of nitrogen

The difference between the average overall relative errors of the training (0.132), and holdout (1.325) samples, must be due to the small data set available, which naturally limits

Training Sum of Squares Error 2.372

Testing Sum of Squares Error 0.081

Holdout Average Overall Relative Error 1.325

In **Table 5** parameter estimates for input and output layer are given for the RBF network.

aCannot be computed. The dependent variable may be constant in the training sample.

Average Overall Relative Error a

Average Overall Relative Error 0.132

Emissions of sulphur oxides

Emissions of nitrogen oxides

Emissions of sulphur oxides

Emissions of nitrogen oxides

Emissions of sulphur oxides

Emissions of nitrogen oxides

0.161

0.103

a

a

1.347

1.267

(million tones of SO2

(million tones of NO2

(million tones of SO2

(million tones of NO2

(million tones of SO2

(million tones of NO2

equivalent)

equivalent)

equivalent)

equivalent)

equivalent)

equivalent)

the possible degree of complexity of the model (Dendek & Mańdziuk, 2008).

Relative Error for Scale

Relative Error for Scale

Relative Error for Scale

Dependents

Dependents

Dependents

Table 4. RBF Model Summary.

Fig. 4. Linear regression of observed values for emissions of nitrogen oxides by predicted values of MLP.

The importance of an independent variable is a measure of how much the network's modelpredicted value changes for different values of the independent variable. A sensitivity analysis to compute the importance of each predictor is applied. The importance chart (**Fig. 5**) shows

Fig. 5. MLP independent variable importance chart.

Fig. 4. Linear regression of observed values for emissions of nitrogen oxides by predicted

The importance of an independent variable is a measure of how much the network's modelpredicted value changes for different values of the independent variable. A sensitivity analysis to compute the importance of each predictor is applied. The importance chart (**Fig. 5**) shows

values of MLP.

Fig. 5. MLP independent variable importance chart.

that the results are dominated by GDP growth rate and GDP (strictly economical QOL indicators), followed distantly by other predictors.

From the RBF analysis, 19 cases (70.4%) were assigned to the training sample, 1 (3.7%) to the testing sample, and 7 (25.9%) to the holdout sample. The seven data records which were excluded from the MLP analysis were excluded from the RBF analysis also, for the same reason.

**Table 4** displays the corresponding information from the RBF network. There appears to be more error in the predictions of emissions of sulphur oxides than in emissions of nitrogen oxides, in the training and holdout samples.

The difference between the average overall relative errors of the training (0.132), and holdout (1.325) samples, must be due to the small data set available, which naturally limits the possible degree of complexity of the model (Dendek & Mańdziuk, 2008).

aCannot be computed. The dependent variable may be constant in the training sample. Table 4. RBF Model Summary.

In **Table 5** parameter estimates for input and output layer are given for the RBF network.

Developing Neural Networks to Investigate

0

0 0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8

values of RBF.

values of RBF.

**Predicted values**

0,5

1

1,5

**Predicted values**

2

2,5

Relationships Between Air Quality and Quality of Life Indicators 255

0 0,511,522,5 **Observed emissions of sulphur oxides**

0 0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 **Observed values of nitrogen oxides**

Fig. 7. Linear regression of observed values for emissions of nitrogen oxides by predicted

Fig. 6. Linear regression of observed values for emissions of sulphur oxides by predicted

Table 5. RBF Parameter Estimates.

Linear regression between observed and predicted values ( 2 2 SO SO *a b error* , NO NO 2 2 *a b error* ) showed that the RBF network does also a reasonably good job of predicting emissions of sulphur and nitrogen oxides. Linear regression gave results for the two output variables 2 2 SO 0.0114 0.8583SO *error* (**Fig. 6**) and NO 0.026 0.7932NO 2 2 *error* (**Fig. 7**), respectively. In this case, it is difficult to see if there is more error in the predictions of emissions of sulphur or nitrogen oxides.

Predicted

1.708 1.517 -1.064 -1.279 -0.562 -0.491 -0.276 -0.204 -0.132

1.572 -0.164 0.575 1.390 -0.448 0.924 -0.544 -0.720 -1.212

H(1) -0.552 -0.668 H(2) 0.463 -0.395 H(3) -0.813 -0.773 H(4) -0.833 -0.795 H(5) -0.617 -0.401 H(6) 0.970 -0.253 H(7) -0.718 -0.547 H(8) 3.116 3.429 H(9) 2.698 2.790

Linear regression between observed and predicted values

reasonably good job of predicting emissions of sulphur and nitrogen oxides. Linear regression gave results for the two output variables 2 2 SO 0.0114 0.8583SO *error*

difficult to see if there is more error in the predictions of emissions of sulphur or nitrogen

) showed that the RBF network does also a

(**Fig. 7**), respectively. In this case, it is

GDP -1.092 -0.986 3.098 0.451 0.667 -0.714 -0.101 -0.076 0.161

Hidden Unit Width 0,606 0.363 0.363 0.363 0.645 0.363 0.576 0.363 0.363

Hidden layer Output layer H(1) H(2) H(3) H(4) H(5) H(6) H(7) H(8) H(9) SO2 NO2

Predictor

Infant Mortality

GDP Growth Rate

Life Expectancy At Birth

Healthy Life Years

Table 5. RBF Parameter Estimates.

, NO NO 2 2 *a b error*

(**Fig. 6**) and NO 0.026 0.7932NO 2 2 *error*

( 2 2 SO SO *a b error*

oxides.

Input Layer

Hidden Layer

Fig. 6. Linear regression of observed values for emissions of sulphur oxides by predicted values of RBF.

Fig. 7. Linear regression of observed values for emissions of nitrogen oxides by predicted values of RBF.

Developing Neural Networks to Investigate

NY, Physica-Verlag.

*Economics* 73, 340–357.

*Health* 1(4), 19-39.

90(10), 1175-1177.

Cambridge.

Callan, R. (1999). *The Essence of Neural Networks*. Prentice Hall, UK .

Information Science, Warsaw, Poland.

Eurostat. (2010). http://epp.eurostat.ec.europa.eu.

Melbourne, Parkville, Australia. IBM. (2010). SPSS Neural Networks 19. SPSS Inc, USA.

Stamoulis editions, Thessaloniki, Greece.

Kecman, V. (2001). *Learning and Soft Computing*. MIT Press, London.

pay for safety. *Structural Safety* 26(2), 181-199. Picton, P. (2000). *Neural Networks, 2nd ed*. Palgrave, New York.

**5. References**

Oxford.

Relationships Between Air Quality and Quality of Life Indicators 257

Bishop, C. (1995). *Neural Networks for Pattern Recognition, 3rd ed*. Oxford University Press,

Bors, A., Pitas, I. (2001). Radial Basis function networks In: Howlett, R., Jain, L (eds.). *Recent*

Bresnahan, B., Mark, D., Shelby, G. (1997). Averting behavior and urban air pollution. *Land*

Dendek, C., Mańdziuk, J. (2008). *Improving Performance of a Binary Classifier by Training Set*

Dogra, Shaillay, K. (2010). *Autoscaling*. QSARWorld - A Strand Life Sciences Web Resource.

Fine, T. (1999). *Feedforward Neural Network Methodology, 3rd ed*. Springer-Verlag, New York. Flynn P., Berry D., Heintz T. (2002). Sustainability & Quality of life indicators: Towards the

Haykin, S. (1998). *Neural Networks: A Comprehensive Foundation, 2nd ed*. Prentice Hall, UK. Hirschberg, J., Esfandiar, M., Slottje, D. (1991). Cluster analysis for measuring welfare and

Hirschberg, J., Maasoumi, E., Slottje, D. (1998). *A cluster analysis the quality of life in the United*

Iliadis, L. (2007). *Intelligent Information Systems and Applications in Risk Management*.

Koller, M., Lorenz, W. (2003). Survival of the quality of life concept. *British Journal of Surgery*

Maasoumi, E. (1998). Multidimensional approaches to welfare. In: Silber, L. (ed.). *Income*

Pandey, M., Nathwani, J. (2004). Life quality index for the estimation of social willingness to

Narasinga Rao, M., Sridhar, G., Madhu, K., Appa Rao, A. (2010). A clinical decision support

diabetes. *Diabetes & Metabolic Syndrome: Clinical Research & Reviews* 4, 57–59. Ripley, B. (1996). *Pattern Recognition and Neural Networks*. Cambridge University Press,

Slottje, D., Scully, G., Hirschberg, J., Hayes, K. (1991). *Measuring the Quality of Life Across Countries: A Multidimensional Analysis*. Westview Press, Boulder, CO.

system using multi-layer perceptron neural network to predict quality of life in

*Inequality Measurement: From Theory to Practice*. Kluwer, New York.

quality of life across countries. *Journal of Econometrics* 50, 131–150.

http://www.qsarworld.com/qsar-statistics-autoscaling.php

*Developments in Theory and Applications in Robust RBF Networks*, 125-153 Heidelberg,

*Selection*. Warsaw University of Technology, Faculty of Mathematics and

Integration of Economic, Social and Environmental Measures. *The Journal of Social*

*States over time*. Department of Economics research paper #596, University of

Finally, the importance chart for the RBF network (**Fig. 8**) shows that, once again, GDP growth rate and GDP are the most important predictors of sulphur and nitrogen oxides emissions.

Fig. 8. RBF Independent variable importance chart.

#### **4. Conclusions**

The multi-layer perceptron and radial-basis function neural network models, that were trained to predict air quality indicators, using life quality and welfare indicators, appear to perform reasonably well. Unlike traditional statistical methods, the neural network models provide dynamic output as further data is fed to it, while they do not require performing and analyzing sophisticated statistical methods (Narasinga Rao et al., 2010).

Results showed that GDP growth rate and GDP influenced mainly air quality predictions, while life expectancy, infant mortality and healthy life years followed distantly. One possible way to ameliorate performance of the network would be to create multiple networks. One network would predict the country result, perhaps simply whether the country increased emissions or not, and then separate networks would predict emissions conditional on whether the country increased emissions. We could then combine the network results to likely obtain better predictions. Note also that neural network is open ended; as more data is given to the model, the prediction would become more reliable. Overall, we find that predictors that include economic indices may be employed by investigators to represent dimensions of air quality that include, as well as go beyond, these simple indices.

#### **5. References**

256 Air Pollution – Monitoring, Modelling and Health

Finally, the importance chart for the RBF network (**Fig. 8**) shows that, once again, GDP growth rate and GDP are the most important predictors of sulphur and nitrogen oxides

The multi-layer perceptron and radial-basis function neural network models, that were trained to predict air quality indicators, using life quality and welfare indicators, appear to perform reasonably well. Unlike traditional statistical methods, the neural network models provide dynamic output as further data is fed to it, while they do not require performing

Results showed that GDP growth rate and GDP influenced mainly air quality predictions, while life expectancy, infant mortality and healthy life years followed distantly. One possible way to ameliorate performance of the network would be to create multiple networks. One network would predict the country result, perhaps simply whether the country increased emissions or not, and then separate networks would predict emissions conditional on whether the country increased emissions. We could then combine the network results to likely obtain better predictions. Note also that neural network is open ended; as more data is given to the model, the prediction would become more reliable. Overall, we find that predictors that include economic indices may be employed by investigators to represent dimensions of air quality that include, as well as go beyond, these

and analyzing sophisticated statistical methods (Narasinga Rao et al., 2010).

emissions.

Fig. 8. RBF Independent variable importance chart.

**4. Conclusions**

simple indices.

**12**

*Brazil*

**Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere**

A relevant issue nowadays is the monitoring and identification of the concentration and rate flux of the gases from the greenhouse effect. Most of these minority gases belong to important bio-geochemical cycles between the planet surface and the atmosphere. Therefore, there is an intense research agenda on this topic. Here, we are going to describe the effort for addressing this challenging. This identification problem can be formulated as an inverse problem. The problem for identifying the minority gas emission rate for the system ground-atmosphere is an important issue for the bio-geochemical cycle, and it has being intensively investigated. This inverse problem has been solved using regularized solutions (Kasibhatla, 2000), Bayes estimation (Enting, 2002; Gimson & Uliasz, 2003), and variational methods (Elbern et al., 2007) – the latter approach coming from the data

The inverse solution could be computed by calculating the regularized solutions. Regularization is a general mathematical procedure for dealing with inverse problems, looking for the smoothest (*regular*) inverse solution. For this approach, the inverse problem can be formulated as an optimization problem with constrains (*a priori*  information). These constrains can be added to the objective function with the help of a regularization parameter. For adding a regularization operator, the ill-posed inverse problem becomes in a well-posed one. The maximum second order entropy principle was used as a regularization operator (Ramos et al., 1999), and the regularization parameter is found by the L-curve technique. A recent approach for solving inverse problems is provided by the use of artificial neural networks (ANN). Neural networks are non-linear mappings, and it is possible to design an ANN to be one inverse operator, with

We have applied two different strategies for addressing the inverse problem: formulated as an optimization problem, and neural network. The optimization problem has been solved

**1. Introduction**

assimilation studies.

robustness to deal on noisy data.

Haroldo F. de Campos Velho1, Débora R. Roberti2,

*1Laboratory for Computing and Applied Mathematics (LAC),*

Eduardo F.P. da Luz1 and Fabiana F. Paes1

*2Department of Physics, Federal University of Santa Maria (UFSM) Santa Maria (RS)*

*National Institute for Space Research (INPE) São José dos Campos (SP)*

### **Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere**

Haroldo F. de Campos Velho1, Débora R. Roberti2, Eduardo F.P. da Luz1 and Fabiana F. Paes1 *1Laboratory for Computing and Applied Mathematics (LAC), National Institute for Space Research (INPE) São José dos Campos (SP) 2Department of Physics, Federal University of Santa Maria (UFSM) Santa Maria (RS) Brazil*

#### **1. Introduction**

258 Air Pollution – Monitoring, Modelling and Health

Tao, K. (1993). A closer look at the radial basis function (RBF) networks. In: Singh, A. (ed.).

*Computers*. IEEE Computational Society Press, Los Alamitos, California . Uykan, Z., Guzelis, C., Celebi, M., Koivo, H. (2000). *Analysis of input-output clustering for determining centers of RBFN*. IEEE Transactions on Neural Networks 11, 851-858.

WHO: World Health Organization. (2010). http://www.who.int.

*Conference Record of the Twenty-Seventh Asilomar Conference on Signals, Systems, and*

A relevant issue nowadays is the monitoring and identification of the concentration and rate flux of the gases from the greenhouse effect. Most of these minority gases belong to important bio-geochemical cycles between the planet surface and the atmosphere. Therefore, there is an intense research agenda on this topic. Here, we are going to describe the effort for addressing this challenging. This identification problem can be formulated as an inverse problem. The problem for identifying the minority gas emission rate for the system ground-atmosphere is an important issue for the bio-geochemical cycle, and it has being intensively investigated. This inverse problem has been solved using regularized solutions (Kasibhatla, 2000), Bayes estimation (Enting, 2002; Gimson & Uliasz, 2003), and variational methods (Elbern et al., 2007) – the latter approach coming from the data assimilation studies.

The inverse solution could be computed by calculating the regularized solutions. Regularization is a general mathematical procedure for dealing with inverse problems, looking for the smoothest (*regular*) inverse solution. For this approach, the inverse problem can be formulated as an optimization problem with constrains (*a priori*  information). These constrains can be added to the objective function with the help of a regularization parameter. For adding a regularization operator, the ill-posed inverse problem becomes in a well-posed one. The maximum second order entropy principle was used as a regularization operator (Ramos et al., 1999), and the regularization parameter is found by the L-curve technique. A recent approach for solving inverse problems is provided by the use of artificial neural networks (ANN). Neural networks are non-linear mappings, and it is possible to design an ANN to be one inverse operator, with robustness to deal on noisy data.

We have applied two different strategies for addressing the inverse problem: formulated as an optimization problem, and neural network. The optimization problem has been solved

Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere 261

swarm intelligence, artificial immune systems and others. Different from other methodologies for computing inverse solutions, neural networks do not need the knowledge on forward problem. In other words, neural networks can be used as *inversion operator* without a mathematical model to describe the direct problem. Several architectures can be adopted to implement a neural network, here will be adopted the multi-layer perceptron with back-propagation algorithm for learning process. Other approaches to implement artificial neural network for inverse problem are described by Campos Velho

Regularization is recognized as the first general method for solving inverse problems, where the regularization operator is a constrain in an optimization problem. Such operator is added to the objective function with the help of a Lagrangian multiplier (also named the *regularization parameter*). The regularization represents an *a priori* knowledge from the physical problem on the unknown function to be estimated. The regularization operator moves the original ill-posed problem to a well posed one. Figure 1 gives an outline of the idea behind the regularization scheme, one of most powerful techniques for computing

Fig. 1. Outline for regularization methods, where additional information is applied to get a

*a priori*  information

Physical reality

The regularization procedure searches for solutions that display *global* regularity. In the mathematical formulation of the method, the inverse problem is expressed as optimization

*Ax f x*

(Tikhonov & Arsenin, 1977). The problem (1) can be written as an optimization problem without constrains with the help of a Lagrange multiplier (penalty or regularization

> <sup>2</sup> <sup>2</sup> min ( ) [ ] *x X Ax f x*

represents the forward problem, and [ ] *x* is the regularization operator

 (2)

(1)

Well-posed problem

2 <sup>2</sup> min ( ) , suject to [ ]

*x X*

+

(2011).

**2.1 Regularization theory**

inverse solutions.

Ill-posed problem

well-posed problem.

problem with constraint:

where *A x*( ) *f*

parameter):

by using a deterministic method (quasi-Newton), implementation used in the E04UCF routine (NAG, 1995) – see Roberti et al. (2007), and applying a stochastic scheme (Particle Swarm Optimization: PSO) (Luz et al., 2007). The PSO is a meta-heuristic based on the collaborative behaviour of biological populations. The algorithm is based on the swarm theory, proposed by Kennedy & Eberhart (1995) and is inspired in the flying pattern of birds which can be achieve by manipulating inter-individual distances to synchronize the swarm behaviour. The swarm theory is then combined with a social-cognitive theory. The multilayer perceptron neural network (MLP-NN) with 2 hidden layers (with 15 and 30 neurons) was applied to estimate the rate of surface emission of a greenhouse gas (Paes et al., 2010). The input for the ANN is the gas concentration measured on a set of points.

### **2. Techniques for solving inverse problems**

There is an aphorism among mathematicians that the theory of partial differential equation could be divided into two branches: an easy part, and a difficult one. The former is called the *forward* (or *direct*) *problem*, and the latter is called *inverse problem*. Of course, the anecdotal statement tries to translate the pathologies found to deal with inverse problem.

For a given mathematical problem to be characterized as a *well-posed problem*, Jacques Hadamard in his studies on differential equations (Hadamard 1902, 1952) has established the following definition:

*A problem is considered well-posed if it satisfies all of three conditions:*

However, inverse problems belong to the class of ill-posed problems, where we cannot guaranty existence, uniqueness, nor stability related to the input data. A mathematical theory for dealing with inverse problems is due to the Russian mathematician Andrei Nikolaevich Tikhonov during the early sixties, introducing the regularization method as a general procedure to solve such mathematical problems. The regularization method looks for the smoothest (*regular*) inverse solution, where the data model would have the best fitting related to the observation data, subjected to the constrains. The searching by the smoothest solution is an additional information (or *a priori* information), which becomes the ill-posed inverse problem in a well-posed one (Tikhonov & Arsenin, 1977).

Statistical strategies to compute inverse solution are those based on Bayesian methods. Instead to compute a set of parameter or function as before, the goal of the statistical methods is to calculate the probability density function (PDF) (Tarantola, 1987; Kaipio & Somersalo, 2005). Bayesian methodology will not be applied here.

With some initial controversy, artificial neural networks are today a well established research field in artificial intelligence (AI). The idea was to develop an artificial device able to emulate the ability of human brain. Sometimes, this view is part of program of so called strong AI. Some authors have used the expression *computational intelligence* (CI) to express a set of methodologies with inspiration on the Nature for addressing real world problems. Initially, the CI includes fuzzy logic systems, artificial neural networks, and evolutionary computation. Nowadays, other techniques are also considered to belong to the CI, such as swarm intelligence, artificial immune systems and others. Different from other methodologies for computing inverse solutions, neural networks do not need the knowledge on forward problem. In other words, neural networks can be used as *inversion operator* without a mathematical model to describe the direct problem. Several architectures can be adopted to implement a neural network, here will be adopted the multi-layer perceptron with back-propagation algorithm for learning process. Other approaches to implement artificial neural network for inverse problem are described by Campos Velho (2011).

#### **2.1 Regularization theory**

260 Air Pollution – Monitoring, Modelling and Health

by using a deterministic method (quasi-Newton), implementation used in the E04UCF routine (NAG, 1995) – see Roberti et al. (2007), and applying a stochastic scheme (Particle Swarm Optimization: PSO) (Luz et al., 2007). The PSO is a meta-heuristic based on the collaborative behaviour of biological populations. The algorithm is based on the swarm theory, proposed by Kennedy & Eberhart (1995) and is inspired in the flying pattern of birds which can be achieve by manipulating inter-individual distances to synchronize the swarm behaviour. The swarm theory is then combined with a social-cognitive theory. The multilayer perceptron neural network (MLP-NN) with 2 hidden layers (with 15 and 30 neurons) was applied to estimate the rate of surface emission of a greenhouse gas (Paes et al., 2010). The input for the ANN is the gas concentration measured on a set of points.

There is an aphorism among mathematicians that the theory of partial differential equation could be divided into two branches: an easy part, and a difficult one. The former is called the *forward* (or *direct*) *problem*, and the latter is called *inverse problem*. Of course, the anecdotal

For a given mathematical problem to be characterized as a *well-posed problem*, Jacques Hadamard in his studies on differential equations (Hadamard 1902, 1952) has established

However, inverse problems belong to the class of ill-posed problems, where we cannot guaranty existence, uniqueness, nor stability related to the input data. A mathematical theory for dealing with inverse problems is due to the Russian mathematician Andrei Nikolaevich Tikhonov during the early sixties, introducing the regularization method as a general procedure to solve such mathematical problems. The regularization method looks for the smoothest (*regular*) inverse solution, where the data model would have the best fitting related to the observation data, subjected to the constrains. The searching by the smoothest solution is an additional information (or *a priori* information), which becomes the

Statistical strategies to compute inverse solution are those based on Bayesian methods. Instead to compute a set of parameter or function as before, the goal of the statistical methods is to calculate the probability density function (PDF) (Tarantola, 1987; Kaipio &

With some initial controversy, artificial neural networks are today a well established research field in artificial intelligence (AI). The idea was to develop an artificial device able to emulate the ability of human brain. Sometimes, this view is part of program of so called strong AI. Some authors have used the expression *computational intelligence* (CI) to express a set of methodologies with inspiration on the Nature for addressing real world problems. Initially, the CI includes fuzzy logic systems, artificial neural networks, and evolutionary computation. Nowadays, other techniques are also considered to belong to the CI, such as

statement tries to translate the pathologies found to deal with inverse problem.

*A problem is considered well-posed if it satisfies all of three conditions:*

c. *The solution has continuous dependence on input data (stability).*

ill-posed inverse problem in a well-posed one (Tikhonov & Arsenin, 1977).

Somersalo, 2005). Bayesian methodology will not be applied here.

**2. Techniques for solving inverse problems**

the following definition:

a. *The solution exists (existence);*  b. *The solution is unique (uniqueness);*  Regularization is recognized as the first general method for solving inverse problems, where the regularization operator is a constrain in an optimization problem. Such operator is added to the objective function with the help of a Lagrangian multiplier (also named the *regularization parameter*). The regularization represents an *a priori* knowledge from the physical problem on the unknown function to be estimated. The regularization operator moves the original ill-posed problem to a well posed one. Figure 1 gives an outline of the idea behind the regularization scheme, one of most powerful techniques for computing inverse solutions.

Fig. 1. Outline for regularization methods, where additional information is applied to get a well-posed problem.

The regularization procedure searches for solutions that display *global* regularity. In the mathematical formulation of the method, the inverse problem is expressed as optimization problem with constraint:

$$\min\_{\mathbf{x}\in\mathcal{X}} \left\| A(\mathbf{x}) - f^{\boldsymbol{\delta}} \right\|\_{2}^{2} \quad \text{subject to } \; \Omega[\mathbf{x}] \le \rho \tag{1}$$

where *A x*( ) *f* represents the forward problem, and [ ] *x* is the regularization operator (Tikhonov & Arsenin, 1977). The problem (1) can be written as an optimization problem without constrains with the help of a Lagrange multiplier (penalty or regularization parameter):

$$\min\_{\mathbf{x}\in\mathcal{X}} \left\| A(\mathbf{x}) - f^{\delta} \right\|\_{2}^{2} + \alpha \Omega[\mathbf{x}] \right\|\tag{2}$$

Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere 263

allowing a greater flexibility to introduce information about the expected shape of the true

For defining higher order entropy functions, two approaches can be derived from the maximization of the entropy of the vector of *first*- and *second*-*differences x*. Under the assumption *x*min < *xi* < *x*max, (*i*=1, …, *N*), the methods of *maximum entropy of order 1* (MaxEnt-1)

(zeroth order)

The non-extensive formulation for the entropy proposed by Tsallis (Tsallis, 1988) can be used to unify the Tikhonov and entropic regularizations (Campos Velho et al., 2006, see also

> 1 1

*q i*

The non-extensive parameter *q* plays a central role in the Tsallis' thermostatistics. The Boltzmann-Gibbs-Shannon's entropy is recovered setting *q*=1. However, for *q*=2, the maximum non-extensive entropy principle is equivalent to the standard Tikhonov regularization. Of course, for *q* 1 and 2, we have a new regularization operator. In context

In order to have a complete theory, the regularization operator should be known, and it is also necessary to have a scheme to compute the regularization parameter. Several methods have been developed to identify this parameter: Morosov's discrepancy criterion, Hansen's method (L-curve method), and generalized cross validation are the most used schemes.

The Morosov's criterion is based on the difference between data of the mathematical model and observations. It should have the same magnitude as measurement errors (Morosov,

2

For determining *N* parameter in an inverse problem, and assuming that measurement errors

<sup>2</sup> *Ax f N* ( )

The discrepancy principle can also be applied to the higher entropy regularization (Muniz et al., 2000). However, if the probability density function of measurement errors follows a

<sup>2</sup> <sup>2</sup>

\*

<sup>2</sup> *Ax f* ( )

*<sup>k</sup> S (x) <sup>x</sup> q*

( ) (first order) 2 2( ) (second order)

1

*i*

*M q*

= 10-10). The argument in the Eq. (4) is to be considered with

(6)

is the error in the measure process,

(7)

2 variance, the discrepancy criterion for

(7.1)

(5)

and *maximum entropy of order 2* (MaxEnt-2) are defined below (Muniz et al., 2000):

1 1max min

1 max min

*x xx x x*

Campos Velho, 2011). The Tsallis' non-extensive entropy is expressed as

of inverse problem, the Boltzmann constant (*k*) is assumed to be equal to 1.

physical model, or its derivatives, for the inverse solution.

min

*i ii*

1984; see also Morosov & Stessin, 1992). Therefore, if

could be modeled by a Gaussian distribution with

independent measures could expressed as:

is the root of the following equation:

is a small parameter (

where

*pi*, instead of *xi*.

*i i ii*

*x x p x xx x*

being the regularization parameter. The first term in the objective function (2) is the fidelity of the model with the observation data, while the second term expresses the regularity (or smoothness) required from the unknown quantity. Note that for 0, the fidelity term is overestimated; on the other hand, for all information in the mathematical model is lost. A definition of a family of regularization operator is given below.

**Definition 1:** A family of continuous regularization operators *RF U* : is called a regularization scheme for the inverse operation of ( ) *A x f* , when

$$\lim\_{\alpha \to 0} \mathcal{R}\_{\alpha} \left\{ A(\mathbf{x}) \right\} = \mathbf{x} \tag{3}$$

where *x X*, and is regularization parameter.

The expression (2) is a practical implementation of the Definition 1. Several regularizations operators have been derived from the pioneer works. Here, only one class of these regularization operators will be described.

#### **2.1.1 Entropic regularization**

The maximum entropy (MaxEnt) principle was firstly proposed as a general inference procedure by Jaynes (1957), on the basis of Shannon's axiomatic characterization of the amount of information (Shannon & Weaver, 1949). It has emerged at the end of the 60's as a highly successful regularization technique. Similar to others regularization techniques, MaxEnt searches for solutions with *global* regularity. Employing a suitable choice of the penalty or regularization parameter, MaxEnt regularization identifies the smoothest reconstructions which are consistent with the measured data. The MaxEnt principle has successfully been applied to a variety of fields: pattern recognition (Fleisher et al., 1990), computerized tomography (Smith et al., 1991), crystallography (de Boissieu et al., 1991), non-destructive testing (Ramos & Giovannini, 1995), magnetotelluric inversion (Campos Velho & Ramos, 1997).

For the vector of parameters *xi* with nonnegative components, the discrete entropy function *S* of vector *u* is defined by

$$\mathbf{S(x)} = -\sum\_{q=1}^{N} \mathbf{s}\_q \log(\mathbf{s}\_q), \quad \text{with} \quad \begin{cases} \mathbf{x} = \begin{bmatrix} \mathbf{x}\_1 & \cdots & \mathbf{x}\_N \end{bmatrix}^T\\ \mathbf{s}\_q = \mathbf{x}\_q \Big/ \sum\_{q=1}^{N} \mathbf{x}\_q \end{cases} \tag{4}$$

The (nonnegative) entropy function *S* attains its global maximum when all *sq* are the same, which corresponds to a uniform distribution with a value of max *S N* log , while the lowest entropy level, min *S* 0 , is attained when all elements *sq* but one are set to zero.

Following the Tikhonov regularization, new entropic higher order regularization techniques have been introduced: MinEnt-1 – applied to identify the 2D electric conductivity from electro-magnetic data (Campos Velho & Ramos, 1997), MaxEnt-2 – employed to the retrieval of vertical profiles of temperature in the atmosphere from remote sensing data (Ramos et al., 1999). They represent a generalization of the standard MaxEnt regularization method,

fidelity of the model with the observation data, while the second term expresses the

mathematical model is lost. A definition of a family of regularization operator is given

lim ( )

The expression (2) is a practical implementation of the Definition 1. Several regularizations operators have been derived from the pioneer works. Here, only one class of these

The maximum entropy (MaxEnt) principle was firstly proposed as a general inference procedure by Jaynes (1957), on the basis of Shannon's axiomatic characterization of the amount of information (Shannon & Weaver, 1949). It has emerged at the end of the 60's as a highly successful regularization technique. Similar to others regularization techniques, MaxEnt searches for solutions with *global* regularity. Employing a suitable choice of the penalty or regularization parameter, MaxEnt regularization identifies the smoothest reconstructions which are consistent with the measured data. The MaxEnt principle has successfully been applied to a variety of fields: pattern recognition (Fleisher et al., 1990), computerized tomography (Smith et al., 1991), crystallography (de Boissieu et al., 1991), non-destructive testing (Ramos & Giovannini, 1995), magnetotelluric inversion (Campos

For the vector of parameters *xi* with nonnegative components, the discrete entropy function

*<sup>T</sup> <sup>N</sup>*

*q q N q q q qq*

1 1

The (nonnegative) entropy function *S* attains its global maximum when all *sq* are the same, which corresponds to a uniform distribution with a value of max *S N* log , while the lowest

Following the Tikhonov regularization, new entropic higher order regularization techniques have been introduced: MinEnt-1 – applied to identify the 2D electric conductivity from electro-magnetic data (Campos Velho & Ramos, 1997), MaxEnt-2 – employed to the retrieval of vertical profiles of temperature in the atmosphere from remote sensing data (Ramos et al., 1999). They represent a generalization of the standard MaxEnt regularization method,

regularity (or smoothness) required from the unknown quantity. Note that for

**Definition 1:** A family of continuous regularization operators *RF U* :

is regularization parameter.

*Sx s s*

( ) log( ), with

entropy level, min *S* 0 , is attained when all elements *sq* but one are set to zero.

fidelity term is overestimated; on the other hand, for

regularization scheme for the inverse operation of ( ) *A x f*

the regularization parameter. The first term in the objective function (2) is the

, when is called a

all information in the

<sup>0</sup> *R Ax x* (3)

*N*

(4)

*xx x*

*sx x* 1

0, the

being

below.

where *x X*, and

**2.1.1 Entropic regularization**

Velho & Ramos, 1997).

*S* of vector *u* is defined by

regularization operators will be described.

allowing a greater flexibility to introduce information about the expected shape of the true physical model, or its derivatives, for the inverse solution.

For defining higher order entropy functions, two approaches can be derived from the maximization of the entropy of the vector of *first*- and *second*-*differences x*. Under the assumption *x*min < *xi* < *x*max, (*i*=1, …, *N*), the methods of *maximum entropy of order 1* (MaxEnt-1) and *maximum entropy of order 2* (MaxEnt-2) are defined below (Muniz et al., 2000):

$$p\_i = \begin{cases} \mathbf{x}\_i - \mathbf{x}\_{\min} + \boldsymbol{\zeta} & \text{(zeroth order)}\\ \mathbf{x}\_{i+1} - \mathbf{x}\_i + (\mathbf{x}\_{\max} - \mathbf{x}\_{\min}) + \boldsymbol{\zeta} & \text{(first order)}\\ \mathbf{x}\_{i+1} - 2\mathbf{x}\_i + \mathbf{x}\_{i-1} + 2(\mathbf{x}\_{\max} - \mathbf{x}\_{\min}) + \boldsymbol{\zeta} & \text{(second order)} \end{cases} \tag{5}$$

where is a small parameter ( = 10-10). The argument in the Eq. (4) is to be considered with *pi*, instead of *xi*.

The non-extensive formulation for the entropy proposed by Tsallis (Tsallis, 1988) can be used to unify the Tikhonov and entropic regularizations (Campos Velho et al., 2006, see also Campos Velho, 2011). The Tsallis' non-extensive entropy is expressed as

$$S\_q(\mathbf{x}) = \frac{k}{q-1} \left[ 1 - \sum\_{i=1}^{M} \mathbf{x}\_i^q \right] \tag{6}$$

The non-extensive parameter *q* plays a central role in the Tsallis' thermostatistics. The Boltzmann-Gibbs-Shannon's entropy is recovered setting *q*=1. However, for *q*=2, the maximum non-extensive entropy principle is equivalent to the standard Tikhonov regularization. Of course, for *q* 1 and 2, we have a new regularization operator. In context of inverse problem, the Boltzmann constant (*k*) is assumed to be equal to 1.

In order to have a complete theory, the regularization operator should be known, and it is also necessary to have a scheme to compute the regularization parameter. Several methods have been developed to identify this parameter: Morosov's discrepancy criterion, Hansen's method (L-curve method), and generalized cross validation are the most used schemes.

The Morosov's criterion is based on the difference between data of the mathematical model and observations. It should have the same magnitude as measurement errors (Morosov, 1984; see also Morosov & Stessin, 1992). Therefore, if is the error in the measure process, is the root of the following equation:

$$\left\| A(\mathbf{x}) - f^{\delta} \right\|\_{\alpha^\*}^2 \approx \delta \tag{7}$$

For determining *N* parameter in an inverse problem, and assuming that measurement errors could be modeled by a Gaussian distribution with 2 variance, the discrepancy criterion for independent measures could expressed as:

$$\left\| A(\mathbf{x}) - f^{\mathcal{S}} \right\|\_{2}^{2} \approx N\sigma^{2} \tag{7.1}$$

The discrepancy principle can also be applied to the higher entropy regularization (Muniz et al., 2000). However, if the probability density function of measurement errors follows a

Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere 265

*nnn nn n bxx u Jx Jx*

4. Compute the search direction *dn* as a solution of the following quadratic programming

<sup>1</sup> <sup>1</sup> minimize: [ ( , )] ( ) <sup>2</sup>

under constrain: ( ) min max ( ) *n n*

6. Test of convergence: Stop, if *x* satisfies the first-order Kuhn-Tucker conditions and

The particle swarm optimization (PSO) is a heuristic search method the uses the behavior of biological collective behavior, like birds, fishes or bees to drive the artificial particles in a search pattern that equalizes a global search with a local search, in a way that an optimum, or near to optimum, result is found (Kennedy & Eberhart, 1995, 2001). There is a sociocognitive theory for supporting the particle swarm scheme. The cultural adaptive process encloses two components: high-level (pattern formation trough the individual and the ability of problem solving), and a low-level (individual behavior), with some skills for

The tendency to evaluate is, maybe, the behavioral feature more present in several living beings. This is intrinsic for the learning process. The comparison among individuals from a population is a behavior process used to identify patterns to find a way to improve each one quality. The emulation behavior gives the living beings the ability to learn through the observations of others actions. This behavior is not so much common in the nature. The emulation behavior is closer to the human society, but it can be used to improve the performance of the algorithms. These three concepts can be combined, even in simplified

PSO searches for the optimum in an infinite (∞) or N (*N*-dimensional space of real numbers,). Indeed, computational solutions can only obtain a projection of an infinite space into N space. The algorithm allows the particles to fly in the search space, looking for the optimization of a predetermined function. A special position *x*\* of a particle in the search space represents the solution of the problem. This position can be represented by an vector of real numbers *x*=(*x*1, *x*2, …, *x*N). In the PSO, the solution *x*\* is searched in a iterative procedure. At each time step, or iteration of the algorithm, the position of the each particle *xk* is updated, by the addition of a speed vector *vk*, in each dimension (space direction). The

( ) ( 1) ( ) *kk k xt xt vt* (10)

approximate the solution of the problem. Otherwise, return to step 1.

1 1 and ( , ) ( , )

*n T n nT n n Jx d dHd*

*qq q qq x xd x x*

*n* minimizes *J*(*xn* +

> *n dn*).

specifies the accuracy to which one wishes to

where

subproblem:

5. Set: *xn*+1 = *xn* +

*n*||*d*|| <

iteration could be given by

**<IP_ADDRESS> Particle swarm optimization**

*ndn*, where the step length

evaluating, comparing, and emulating (Kennedy & Eberhart, 2001).

computational rules, providing a scheme for solving optimization problems.

1/2(1+||*x*||), where

distribution where the second statistical moment is not defined (Lévy or Cauchy distributions, for example) the *generalized Morosov's discrepancy principle* can be used (Shiguemori et al., 2004).

If the statistics on the observational data is not available, the generalized cross-validation method can be applied (Bertero & Bocacci, 1998; Aster et al., 2005). Considering now a linear forward problem: *A*(*x*)*Ax*, being *A* a matrix, the goal of the cross-validation scheme is to minimize the *generalized cross-validation function* (Aster et al., 2005):

$$V(\alpha) = \frac{N \left\| A(\mathbf{x}\_{\alpha}) - f^{\mathcal{S}} \right\|\_{2}^{2}}{\left\lceil Tr\left\{ I - B(\alpha) \right\} \right\rceil^{2}} \tag{8}$$

being *Tr*{*C*} the trace of matrix *C*, and *B*( ) is the following matrix:

$$B(\alpha) \equiv AA^\* \left(AA^\* + \alpha I\right)^{-1} \tag{9}$$

where *A*\* is the *adjoint matrix*: \* ( ,) (, ) *Af g f A g* and *I* is the identity matrix.

Another scheme to compute the regularization parameter is the L-curve method. The Lcurve criterion is a geometrical approach suggested by Hansen (1992) (see also Bertero & Bocacci, 1998; Aster et al., 2005). The idea is to find the point of maximum curvature, on the corner of the plot: 2 <sup>2</sup> [ ] () *x Ax f* . In general, the plot smoothness × fidelity shows a Lshape curve type.

**2.1.2 Optimization methods**

This is a central issue in science and engineering, where several methods have been developed. They can classified as a deterministic and stochastic schemes. Here, only the methods used in the application treated will be described.

#### **<IP_ADDRESS> Quasi-Newton optimization**

The minimization of the objective function *J*(*x*) given by equation (2), subjected to simple bounds on *x*, is solved using a first-order optimization algorithm from the E04UCF routine, NAG Fortran Library (1993). This routine is designed to minimize an arbitrary smooth function subject to constraints (simple bounds, linear and nonlinear constraints), using a sequential programming method. For the *n*-th iteration, the calculation proceeds as follows:

$$H^n = H^{n-1} + \frac{b^n (b^n)^T}{(b^n)^T u^n} - \frac{H^{n-1} u^n (u^n)^T H^{n-1}}{(u^n)^T H^{n-1} u^n}$$

where

264 Air Pollution – Monitoring, Modelling and Health

distribution where the second statistical moment is not defined (Lévy or Cauchy distributions, for example) the *generalized Morosov's discrepancy principle* can be used

If the statistics on the observational data is not available, the generalized cross-validation method can be applied (Bertero & Bocacci, 1998; Aster et al., 2005). Considering now a linear forward problem: *A*(*x*)*Ax*, being *A* a matrix, the goal of the cross-validation scheme is to

( )

( )

*N Ax f*

*Tr I B*

\*\* 1 () ( )

Another scheme to compute the regularization parameter is the L-curve method. The Lcurve criterion is a geometrical approach suggested by Hansen (1992) (see also Bertero & Bocacci, 1998; Aster et al., 2005). The idea is to find the point of maximum curvature, on the

This is a central issue in science and engineering, where several methods have been developed. They can classified as a deterministic and stochastic schemes. Here, only the

The minimization of the objective function *J*(*x*) given by equation (2), subjected to simple bounds on *x*, is solved using a first-order optimization algorithm from the E04UCF routine, NAG Fortran Library (1993). This routine is designed to minimize an arbitrary smooth function subject to constraints (simple bounds, linear and nonlinear constraints), using a sequential programming method. For the *n*-th iteration, the calculation proceeds

2 2 2

(8)

,*x*),

is the following matrix:

*x Ax f* . In general, the plot smoothness × fidelity shows a L-

1 1

() ()

*nT n nT n n*

*bu uHu*

() ()

 *n nT n n nT n n n*

*bb H uu H H H*

1

*B AA AA I* (9)

minimize the *generalized cross-validation function* (Aster et al., 2005):

( )

where *A*\* is the *adjoint matrix*: \* ( ,) (, ) *Af g f A g* and *I* is the identity matrix.

2

1. Solve the forward problem to *x* and compute the objective function *J*(

1

,*x*) by finite difference, 3. Compute a positive-definite approximation (quasi-Newton) to the Hessian *Hn*:

<sup>2</sup> [ ] ()

methods used in the application treated will be described.

*V*

(Shiguemori et al., 2004).

corner of the plot:

shape curve type.

as follows:

**2.1.2 Optimization methods**

**<IP_ADDRESS> Quasi-Newton optimization**

2. Compute the gradient *J*(

being *Tr*{*C*} the trace of matrix *C*, and *B*( )

$$b^{n} = \mathbf{x}^{n} - \mathbf{x}^{n-1} \quad \text{and} \quad \boldsymbol{u}^{n} = \nabla I(\boldsymbol{\alpha}, \mathbf{x}^{n}) - \nabla I(\boldsymbol{\alpha}, \mathbf{x}^{n-1})$$

4. Compute the search direction *dn* as a solution of the following quadratic programming subproblem:

$$\text{minimize:} \left[\nabla f(\boldsymbol{\alpha}, \boldsymbol{\alpha}^{n})\right]^{T} \boldsymbol{d}^{n} + \frac{1}{2} (\boldsymbol{d}^{n})^{T} \boldsymbol{H}^{n-1} \boldsymbol{d}^{n}$$

under constrain: ( ) min max ( ) *n n qq q qq x xd x x*

#### **<IP_ADDRESS> Particle swarm optimization**

The particle swarm optimization (PSO) is a heuristic search method the uses the behavior of biological collective behavior, like birds, fishes or bees to drive the artificial particles in a search pattern that equalizes a global search with a local search, in a way that an optimum, or near to optimum, result is found (Kennedy & Eberhart, 1995, 2001). There is a sociocognitive theory for supporting the particle swarm scheme. The cultural adaptive process encloses two components: high-level (pattern formation trough the individual and the ability of problem solving), and a low-level (individual behavior), with some skills for evaluating, comparing, and emulating (Kennedy & Eberhart, 2001).

The tendency to evaluate is, maybe, the behavioral feature more present in several living beings. This is intrinsic for the learning process. The comparison among individuals from a population is a behavior process used to identify patterns to find a way to improve each one quality. The emulation behavior gives the living beings the ability to learn through the observations of others actions. This behavior is not so much common in the nature. The emulation behavior is closer to the human society, but it can be used to improve the performance of the algorithms. These three concepts can be combined, even in simplified computational rules, providing a scheme for solving optimization problems.

PSO searches for the optimum in an infinite (∞) or N (*N*-dimensional space of real numbers,). Indeed, computational solutions can only obtain a projection of an infinite space into N space. The algorithm allows the particles to fly in the search space, looking for the optimization of a predetermined function. A special position *x*\* of a particle in the search space represents the solution of the problem. This position can be represented by an vector of real numbers *x*=(*x*1, *x*2, …, *x*N). In the PSO, the solution *x*\* is searched in a iterative procedure. At each time step, or iteration of the algorithm, the position of the each particle *xk* is updated, by the addition of a speed vector *vk*, in each dimension (space direction). The iteration could be given by

$$\mathbf{x}^{k}(t) = \mathbf{x}^{k}(t-1) + \boldsymbol{\upsilon}^{k}(t) \tag{10}$$

Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere 267

weight adjustments are conducted by backpropagating such error to the network, governed by a change rule. The weights are changed by an amount proportional to the error at that unit, times the output of the unit feeding into the weight. Equation below shows the general

[ ( ) ( 1)] [ ( ) ( 1)] *ji i ij ij i i ij ij*

decides to what degree this previous adjustment is to be considered so as to prevent any sudden changes in the direction in which corrections are made. The learning rate

important procedure, indicating the performance of a NN. The effective test is defined using an unknown vector *x* that did not belong to the training set. This action is called the

The up date for the weight matrix of neural connections is obtained with the innovation

( ) ( 1) *wn wn w ij ij j <sup>i</sup>*

where the iterative process follows up to reach a maximum number of epochs (iterations) or the required precision to be caught. The delta rule can be justified if the diagonal matrix

(the learning parameter multiply by identity matrix) could be understood as an approximate inverse of Hessian matrix (a matrix *C* is defined as an approximate inverse if: ||*I* - *CA*|| <

were set to 0.1 and 0.5, respectively. Nevertheless, the activation test is an

the momentum coefficient, which

is the learning rate

and

> *I*

(13)

*w x wn wn x wn wn*

*j* is the local gradient, *xi* is the input signal of neuron *j*,

weight correction according to the so-called delta rule (Haykin, 1999)

Fig. 2. Outline of fully connected MLP-NN with 2 hidden layers.

*ji*

parameter that controls the strength of change, and

*w*

where,

momentum

*generalization* of the NN.

described in Equation (13):

where *xk*(*t*-1) is the previous position of the particle. The velocity that is applied to the particle is function of two basic behaviors: a cognitive component, and a social component. The velocity is up dated taken into account the behaviors mentioned, as following:

$$\upsilon^k(t) = w\upsilon^k(t-1) + c\_1 r n m d\_1(\upsilon^k - \mathfrak{x}^k) + c\_1 r n m d\_2(\upsilon^g - \mathfrak{x}^k) \tag{11}$$

where *c*1,2 are constants (real numbers) that weights the cognitive and social components, respectively, *rand*1,2 are random numbers from a uniform distribution in the interval [0,1], *pk* is the best position achieved by the particle and *pg* is the best position of the swarm, with *xk* actual position of the particle *k*. The parameter *w* is used to give more stability to the algorithm (Shi & Eberhart, 1998). The selection of the *w* value follows a heuristic rule:

#### **2.2 Artificial neural networks**

The artificial neuron is a simplified representation of a biological neuron, where a weighted combination of the inputs is the value for the non-linear activation function (in general, a sigmoid one). In artificial neural networks (ANN), the identification of the connection weights is called the learning (or training) phase.

Several architectures can be used for the ANN. Related to the learning process, ANNs can be described into two important groups: supervised, and unsupervised NNs. For the supervised scheme, the connection weights are selected to become the output from the ANN as close as possible to the target. For example: the weights could be found by minimizing a functional of the square differences between the ANN output and the target values. Some authors have suggested the use of a Tikhonov's functional (Poggio & Girosi, 1990). The regularization suggested by Poggio & Girosi (1990) is applied on the output values *x* of the NN:

$$f[w] = \sum\_{i=1}^{N} [\mathbf{x}\_i^{\text{Target}} - \mathbf{x}\_i(w)]^2 + \alpha \left\| P\mathbf{x} \right\|\tag{12}$$

where *x* is the NN output, *w* is the connection weight matrix, and *P* is a linear operator. Shiguemori et al. (2004) employed three different architectures for the NNs, with focus on determination of the initial profile for the heat conduction problem. The NNs used were: multilayer perceptron (MLP), RBF-NN, and cascade correlation. The backpropagation algorithm was used for the learning process (without regularization). Figure 2 shows a cartoon to represent a MLP-NN. Cybenko (1989) has proved that the standard multilayer NN with a single hidden layer with a sigmoid activation function, for a finite number of hidden neurons, are universal approximators on a compact subset of N (Cybenko theorem).

The backpropagation training is a supervised learning algorithm that requires both input and output (desired) data. Such pairs permit the calculation of the error of the network as the difference between the calculated output and the desired vector: 1/2 *x*target – *x*(*w*). The

where *xk*(*t*-1) is the previous position of the particle. The velocity that is applied to the particle is function of two basic behaviors: a cognitive component, and a social component.

11 12 ( ) ( 1) ( ) ( )

where *c*1,2 are constants (real numbers) that weights the cognitive and social components, respectively, *rand*1,2 are random numbers from a uniform distribution in the interval [0,1], *pk* is the best position achieved by the particle and *pg* is the best position of the swarm, with *xk* actual position of the particle *k*. The parameter *w* is used to give more stability to the

If *w* > 1.2: the algorithm works in a exploitation mode, where the global is more

The artificial neuron is a simplified representation of a biological neuron, where a weighted combination of the inputs is the value for the non-linear activation function (in general, a sigmoid one). In artificial neural networks (ANN), the identification of the connection

Several architectures can be used for the ANN. Related to the learning process, ANNs can be described into two important groups: supervised, and unsupervised NNs. For the supervised scheme, the connection weights are selected to become the output from the ANN as close as possible to the target. For example: the weights could be found by minimizing a functional of the square differences between the ANN output and the target values. Some authors have suggested the use of a Tikhonov's functional (Poggio & Girosi, 1990). The regularization suggested by Poggio & Girosi (1990) is applied on the output values *x* of the

Target 2

(12)

1/2 *x*target – *x*(*w*). The

*i i*

*J w x x w Px*

where *x* is the NN output, *w* is the connection weight matrix, and *P* is a linear operator. Shiguemori et al. (2004) employed three different architectures for the NNs, with focus on determination of the initial profile for the heat conduction problem. The NNs used were: multilayer perceptron (MLP), RBF-NN, and cascade correlation. The backpropagation algorithm was used for the learning process (without regularization). Figure 2 shows a cartoon to represent a MLP-NN. Cybenko (1989) has proved that the standard multilayer NN with a single hidden layer with a sigmoid activation function, for a finite number of hidden neurons,

The backpropagation training is a supervised learning algorithm that requires both input and output (desired) data. Such pairs permit the calculation of the error of the network as

1 [ ] [ ( )]

*i*

are universal approximators on a compact subset of N (Cybenko theorem).

the difference between the calculated output and the desired vector:

*N*

If 0.8 < *w* < 1.2: a balance point, where exploration and exploitation are well divided.

*k k k k <sup>g</sup> <sup>k</sup> v t wv t c rand p x c rand p x* (11)

The velocity is up dated taken into account the behaviors mentioned, as following:

algorithm (Shi & Eberhart, 1998). The selection of the *w* value follows a heuristic rule:

If *w* < 0.8: the algorithm works in a local search mode, in a exploration mode;

explored;

NN:

**2.2 Artificial neural networks**

weights is called the learning (or training) phase.

weight adjustments are conducted by backpropagating such error to the network, governed by a change rule. The weights are changed by an amount proportional to the error at that unit, times the output of the unit feeding into the weight. Equation below shows the general weight correction according to the so-called delta rule (Haykin, 1999)

$$
\Delta w\_{\vec{\eta}\vec{\iota}} = \eta \frac{\partial \mathcal{x}}{\partial w\_{\vec{\eta}\vec{\iota}}} \mathbf{x}\_i + \gamma \{w\_{\vec{\eta}}(\mathbf{n}) - w\_{\vec{\eta}}(\mathbf{n} - 1)\} = \eta \delta\_i \mathbf{x}\_i + \gamma \{w\_{\vec{\eta}}(\mathbf{n}) - w\_{\vec{\eta}}(\mathbf{n} - 1)\} \tag{13}
$$

where, *j* is the local gradient, *xi* is the input signal of neuron *j*, is the learning rate parameter that controls the strength of change, and the momentum coefficient, which decides to what degree this previous adjustment is to be considered so as to prevent any sudden changes in the direction in which corrections are made. The learning rate and momentum were set to 0.1 and 0.5, respectively. Nevertheless, the activation test is an important procedure, indicating the performance of a NN. The effective test is defined using an unknown vector *x* that did not belong to the training set. This action is called the *generalization* of the NN.

Fig. 2. Outline of fully connected MLP-NN with 2 hidden layers.

The up date for the weight matrix of neural connections is obtained with the innovation described in Equation (13):

$$w\_{ij}(n) = w\_{ij}(n-1) + \Delta w\_{ji}$$

where the iterative process follows up to reach a maximum number of epochs (iterations) or the required precision to be caught. The delta rule can be justified if the diagonal matrix *I* (the learning parameter multiply by identity matrix) could be understood as an approximate inverse of Hessian matrix (a matrix *C* is defined as an approximate inverse if: ||*I* - *CA*|| <

Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere 269

& Ferrero, 1997). The diffusion coefficients, ( , , ) *ij b t* **x u** , for both forward and backward

2 *<sup>i</sup>*

the Lagrangian time scale (Degrazia et al., 2000), respectively. With the coordinates and the

The inverse problem here is to identify the source term *S*(*t*). A source-receptor approach can be employed for reducing the computer time, instead of running the direct model, Eq. (14), for each iteration. This approach displays an explicit relation between the pollutant

1

where *VR i*, and *VS*, *<sup>j</sup>* are the volume for the *i*-th receptor and *j*-th source, respectively; *NS*, *<sup>j</sup>* and *NR i*, are the number of particle realised by the *j*-th source and *i*-th sensor, respectively; *NR i*, , *<sup>j</sup>* and *NS i*, , *<sup>j</sup>* are the number of particle released by the *j*-th source and detected by the

For testing the procedure to estimate the emission rate procedure, it is considered the area pollutant sources placed in a box volume, where the horizontal domain and vertical height are given by: (1500 m 1000 m) 1000 m. There are two embedded regions R1 and R2 into computational domain, with following horizontal domain (600 m 600 m) for each region, and 1 m of height, and they are realising contaminants with two different emission rates. Figure 1 shows the computational scenario in a two-dimensional projection ( , *x y* ): the six sensors are placed at 10 m height and spread horizontally with the coordinates presented in Table 1. The domain is divided into sub-domains with 200 m 200 m 1 m. The emission

R1 = *AAAAAA* <sup>234789</sup> *SSSSSS* = 10 gm-3s-1 (Region-1)

R2 = *AAAAAA* 12 13 14 17 18 19 *SSSSSS* = 20 gm-3s-1 (Region-2)

*Ns i ij j j C MS*

*Li*

*ij ij*

 *<sup>i</sup>* and *Li*

*b*

where the matrix *Mij* is the *transition matrix*, and each matrix entry given by

*V V tN N*

, ,,

**4. Determining gas flux between the ground and the atmosphere**

*Ri Si j*

, , , ,,

*S j Ri S j Ri j*

*tN N*

*ij*

rates for each sub-domain ( *Ak S* ) is as following:

*M*

mass of each particle, the concentration is computed (see below).

concentration of the *i*-th receptor related the *j*-th sources:

1/2 <sup>2</sup>

(16)

are velocity variance at each component and

(17)

(18)

integration is given by

is the Kronecker delta, <sup>2</sup>

where *ij*

*i*-th receptor.

1, see Conte and de Boor (1980)). Therefore, the free parameter should determined to be an approximation for the inverse of the Hessian matrix, producing similar results as the Newton method. However, the Newton method, or even the steepest descent scheme, is not good to deal if a shallow local minimum. The momentum term ( 0, in Eq. (13)) is added to circumvent some instabilities (Rumelhartet al., 1986) and avoid pathological behaviour with shallow error surface (Haykin, 1999). However, any other optimization method can be used instead of delta rule. Indeed, the PSO scheme described in Section <IP_ADDRESS> can be employed to determine the weights of the neural network (Li and Chen, 2006).

#### **3. The lagrangian stochastic model LAMBDA**

The Lagrangian particle model LAMBDA was developed to study the transport process and pollutants diffusion, starting from the Brownian random walk modeling (Anfossi & Ferrero, 1997). In the LAMBDA code, full-uncoupled particle movements are assumed. Therefore, each particle trajectory can be described by the generalized three dimensional form of the Langevin equation for velocity (Thomson, 1987):

$$du\_i = a\_i(\mathbf{x}, \mathbf{u}, t)dt + b\_{ij}(\mathbf{x}, \mathbf{u}, t)dW\_j(\mathbf{t})\tag{14a}$$

$$d\mathbf{x} = (\mathbf{U} + \mathbf{u})dt\tag{14b}$$

where *i j* , 1,2,3 , and **x** is the displacement vector, **U** is the mean wind velocity vector, **u** is the Lagrangian velocity vector, ( , ,) *<sup>i</sup> a t* **x u** is a deterministic term, and ( , ,) () *ij <sup>j</sup> b t dW t* **x u** is a stochastic term and the quantity ( ) *<sup>j</sup> dW t* is the incremental Wiener process.

The deterministic (drift) coefficient ( , ,) *<sup>i</sup> a t* **x u** is computed using a particular solution of the Fokker-Planck equation associated to the Langevin equation. The diffusion coefficient ( , ,) *ij b t* **x u** is obtained from the Lagrangian structure function in the inertial subrange, ( ) *K L t* , where *<sup>K</sup>* is the Kolmogorov time scale and *<sup>L</sup>* is the Lagrangian decorrelation time scale. These parameters can be obtained employing the Taylor statistical theory on turbulence (Degrazia et al., 2000).

The drift coefficient, ( , ,) *<sup>i</sup> a t* **x u** , for forward and backward integration is given by

$$
\rho\_i P\_E = \frac{\partial}{\partial u\_j} \left( B\_{i,j} P\_E \right) + \varphi\_i \left( \mathbf{x}\_i \mathbf{u}\_i t \right) \tag{15a}
$$

with,

$$\frac{\partial \varphi\_i}{\partial u\_i} = -\frac{\partial P\_E}{\partial t} - \frac{\partial}{\partial \mathbf{x}\_i} (\mathbf{u}\_i P\_E) \text{ and } \varphi\_i \to 0 \text{ when } \mathbf{u} \to \infty \tag{15b}$$

where 1 *<sup>v</sup> c* for forward integration and 1 *<sup>v</sup> c* for backward integration, *PP t <sup>E</sup>* **x u**, , is the non-conditional PDF of the Eulerian velocity fluctuations, and , , , (1 /2) *B bb <sup>i</sup> <sup>j</sup> i k <sup>j</sup> <sup>k</sup>* . Of course, for backward integration, the time considered is ' *t t* , and velocity ' **U U** , being **U** the mean wind speed. The horizontal PDFs are considered Gaussians, and for the vertical coordinate the truncated Gram-Charlier type-C of third order is employed (Anfossi

approximation for the inverse of the Hessian matrix, producing similar results as the Newton method. However, the Newton method, or even the steepest descent scheme, is not

circumvent some instabilities (Rumelhartet al., 1986) and avoid pathological behaviour with shallow error surface (Haykin, 1999). However, any other optimization method can be used instead of delta rule. Indeed, the PSO scheme described in Section <IP_ADDRESS> can be employed

The Lagrangian particle model LAMBDA was developed to study the transport process and pollutants diffusion, starting from the Brownian random walk modeling (Anfossi & Ferrero, 1997). In the LAMBDA code, full-uncoupled particle movements are assumed. Therefore, each particle trajectory can be described by the generalized three dimensional form of the

, ,) ( , ,) () *i i ij <sup>j</sup> du a t dt b t dW t* **(x u x u** (14a)

where *i j* , 1,2,3 , and **x** is the displacement vector, **U** is the mean wind velocity vector, **u** is the Lagrangian velocity vector, ( , ,) *<sup>i</sup> a t* **x u** is a deterministic term, and ( , ,) () *ij <sup>j</sup> b t dW t* **x u**

The deterministic (drift) coefficient ( , ,) *<sup>i</sup> a t* **x u** is computed using a particular solution of the Fokker-Planck equation associated to the Langevin equation. The diffusion coefficient ( , ,) *ij b t* **x u** is obtained from the Lagrangian structure function in the inertial subrange,

is the Kolmogorov time scale and *<sup>L</sup>*

, for forward and backward integration is given by

 and 0

correlation time scale. These parameters can be obtained employing the Taylor statistical

, , , *i E ij E i*

where 1 *<sup>v</sup> c* for forward integration and 1 *<sup>v</sup> c* for backward integration, *PP t <sup>E</sup>* **x u**, , is the non-conditional PDF of the Eulerian velocity fluctuations, and , , , (1 /2) *B bb <sup>i</sup> <sup>j</sup> i k <sup>j</sup> <sup>k</sup>* . Of course, for backward integration, the time considered is ' *t t* , and velocity ' **U U** , being **U** the mean wind speed. The horizontal PDFs are considered Gaussians, and for the vertical coordinate the truncated Gram-Charlier type-C of third order is employed (Anfossi

*aP B P t*

*j*

*i E*

*u*

*<sup>P</sup> u P*

 *i E*

*i i*

*u tx*

is a stochastic term and the quantity ( ) *<sup>j</sup> dW t* is the incremental Wiener process.

*d dt* **x Uu** ( ) (14b)

**x u** (15a)

*<sup>i</sup>* when **u** (15b)

is the Lagrangian de-

should determined to be an

0, in Eq. (13)) is added to

1, see Conte and de Boor (1980)). Therefore, the free parameter

good to deal if a shallow local minimum. The momentum term (

to determine the weights of the neural network (Li and Chen, 2006).

**3. The lagrangian stochastic model LAMBDA**

Langevin equation for velocity (Thomson, 1987):

( ) *K L*

The drift coefficient, ( , ,) *<sup>i</sup> a t* **x u**

, where *<sup>K</sup>*

theory on turbulence (Degrazia et al., 2000).

*t*

with,

& Ferrero, 1997). The diffusion coefficients, ( , , ) *ij b t* **x u** , for both forward and backward integration is given by

$$b\_{ij} = \delta\_{ij} \left[ 2 \frac{\sigma\_i^2}{\sigma\_{Li}} \right]^{1/2} \tag{16}$$

where *ij* is the Kronecker delta, <sup>2</sup> *<sup>i</sup>* and *Li* are velocity variance at each component and the Lagrangian time scale (Degrazia et al., 2000), respectively. With the coordinates and the mass of each particle, the concentration is computed (see below).

The inverse problem here is to identify the source term *S*(*t*). A source-receptor approach can be employed for reducing the computer time, instead of running the direct model, Eq. (14), for each iteration. This approach displays an explicit relation between the pollutant concentration of the *i*-th receptor related the *j*-th sources:

$$\mathbf{C}\_{i} = \sum\_{j=1}^{N\_{s}} \mathbf{M}\_{ij} \mathbf{S}\_{j} \tag{17}$$

where the matrix *Mij* is the *transition matrix*, and each matrix entry given by

$$\mathbf{M}\_{ij} = \begin{cases} \left(\mathbf{V}\_{\mathbf{S},j} / \mathbf{V}\_{\mathbf{R},i}\right) \left(\Delta t / \mathbf{N}\_{\mathbf{S},j}\right) \mathbf{N}\_{\mathbf{R},i,j} & \text{-forward model;}\\ \left(\Delta t / \mathbf{N}\_{\mathbf{R},i}\right) \mathbf{N}\_{\mathbf{S},i,j} & \text{-backward model.} \end{cases} \tag{18}$$

where *VR i*, and *VS*, *<sup>j</sup>* are the volume for the *i*-th receptor and *j*-th source, respectively; *NS*, *<sup>j</sup>* and *NR i*, are the number of particle realised by the *j*-th source and *i*-th sensor, respectively; *NR i*, , *<sup>j</sup>* and *NS i*, , *<sup>j</sup>* are the number of particle released by the *j*-th source and detected by the *i*-th receptor.

#### **4. Determining gas flux between the ground and the atmosphere**

For testing the procedure to estimate the emission rate procedure, it is considered the area pollutant sources placed in a box volume, where the horizontal domain and vertical height are given by: (1500 m 1000 m) 1000 m. There are two embedded regions R1 and R2 into computational domain, with following horizontal domain (600 m 600 m) for each region, and 1 m of height, and they are realising contaminants with two different emission rates. Figure 1 shows the computational scenario in a two-dimensional projection ( , *x y* ): the six sensors are placed at 10 m height and spread horizontally with the coordinates presented in Table 1. The domain is divided into sub-domains with 200 m 200 m 1 m. The emission rates for each sub-domain ( *Ak S* ) is as following:

$$\mathbf{R\_1} = \mathbf{S\_{A\_2}} = \mathbf{S\_{A\_3}} = \mathbf{S\_{A\_4}} = \mathbf{S\_{A\_7}} = \mathbf{S\_{A\_8}} = \mathbf{S\_{A\_9}} = \mathbf{1} \text{ 9 cm} \cdot \text{s} \text{s}^{-1} \text{ (Region-1)}$$

$$\mathbf{R\_2} = \mathbf{S\_{A\_{12}}} = \mathbf{S\_{A\_{13}}} = \mathbf{S\_{A\_{14}}} = \mathbf{S\_{A\_{17}}} = \mathbf{S\_{A\_{18}}} = \mathbf{S\_{A\_{19}}} = \mathbf{20 } \text{gm} \cdot \text{s} \text{s}^{-1} \text{ (Region-2)}$$

Strategies for Estimation of Gas Mass Flux Rate Between Surface and the Atmosphere 271

(h:min) 10 m 120 m 200 m 10 m 20 m 200 m 12:05 2.6 5.7 5.7 290 310 310 12:15 2.6 5.1 5.7 300 310 310 12:25 2.1 4.6 5.1 280 310 320 12:35 2.1 4.6 5.1 280 310 320 12:45 2.6 5.1 5.7 290 310 310 Table 2. Measured wind speed for the Copenhagen experiment (Gryning & Lyck, 1998).

The inverse approach is tested using synthetic measured concentration data The synthetic

random number associated to the Gaussian distribution with zero mean and unitary

2 6

where the vector 1 2 <sup>25</sup> [ , ,..., ]*<sup>T</sup>* **<sup>S</sup>** *SS S AA A* represent the emission rate from 25 sub-domains, presented in Figure 3. However, some sub-domains are considered releasing nothing. Therefore, the number of the unknown parameters to be estimated by inverse method decreases to 12 cells: 2 3 4 7 8 9 12 13 14 17 18 19 [,,,,,, , , , , , ] **S** *SSSSSSS S S S S S AAAAAAA A A A A A* . Another *a priori* information is smoothness of the inverse solution (regularization operator). As mentioned before, the second order maximum entropy is applied, Eq. (13c). The use of entropic regularization for this type of inverse problem was independently proposed by Roberti (2005) and Bocquet (2005a, 2005b) – see also Roberti et al. (2007) and Davoine &

There are several methods to compute the regularization parameter. Here, a numerical experimentation was employed to determine this parameter. Roberti (2005) had found the

The first inverse solution was obtained by applying the quasi-Newton optimization method (described in Section 2.1.1). In the Experiment-1, a global error around 2% is found for the inverse solution. For the region-1 (lower emission rate), the error in the estimation is 10%, and for region-2 the error is 3%. The biggest errors for estimation are verified for sub-

= 10-5 for Experiment-2.

, ( )

where *CMod*(*x*,*y*) is the gas concentration computed from the forward model (14),

experiments are performed, using two different maximum noise levels: 5% and 10%.

 *Exp Mod j j*

The optimization method is used to find the minimum for the functional:

*j F CC*

domains A9 and A14, with error of 55% and 42%, respectively.

1

is the level of the noise (Gaussian white noise). In our tests,

=0.1 (experiment-2) were used. In other words, two numerical

**S S <sup>S</sup>** (22)

(21)

is a

=0.05

Time Speed *U* (ms-1) Wind direction ()

( , ) ( , )[1 ] *Exp Mod C xy C xy jj jj*

= 10-6 for Experiment-1, and

observations are emulated as:

standard deviation, and

(Experiment-1) and

Bocquet (2007).

values

$$S\_{A\_1} = S\_{A\_5} = S\_{A\_6} = S\_{A\_{10}} = S\_{A\_{11}} = S\_{A\_{15}} = S\_{A\_{16}} = S\_{A\_{20}} = S\_{A\_{21}} = S\_{A\_{22}} = S\_{A\_{26}} = S\_{A\_{26}} = S\_{A\_{25}} = 0$$

The dispersion problem is modelled by the LAMBDA backward-time code, where data from the Copenhagen experiment are used, for the period (12:33 h – 12:53 h, 19/October/1978) (Gryning & Lyck, 1984, 1998). It is assumed a logarithmic vertical profile for the wind field

$$\text{LU}(z) = \frac{u\_\*}{\kappa} \ln \left( \frac{z}{z\_0} \right) \tag{19}$$

being *U*(*z*) the main stream, *u*\* the fiction velocity, the von Kármán constant, *z* the height above the ground, and *z*0=0.06 m the roughness. The wind speed was measured at 10 m, 120 m, and 200 m. The numerical value for *u*\* was obtained from the best fitting with the measured wind speed – see Table 2 – and the equation (19).

Table 1. Position of sensors in physical domain.

For this experiment the wind direction had an angle =180°, and the boundary layer height is *h*=1120 m (Gryning & Lyck, 1998). The turbulence parameterization follows Degrazia et al. (2000), for computing the wind variance ( <sup>2</sup> *<sup>i</sup>* ) and the Lagrangian decorrelation time scale ( *Li T* ). These two turbulent parameters are considered constant for the whole boundary layer, and their numerical values were calculate at *z*=10 m level. This characterises a stationary and vertically homogeneous atmosphere.

The sensors dimension, where the fictitious particles have arrived, was de 0,1 m 0,1 m 0,1 m, centred in the computational cells presented in Table 1. Next, the reverse trajectories are calculated for 1000 particles per sensor. The parameters for numerical simulations were 1800 time-steps with *t* = 1 s, meaning 30 min for the whole simulation. After 10 min, the concentration was computed at each 2 min, for the remaining 20 min of simulation. The mean concentration is found using the following expression:

$$\mathbf{C}^{\text{Mod}}(\mathbf{x}\_{j}, \mathbf{y}\_{j}) = \frac{1}{10} \sum\_{n=1}^{10} \left[ \sum\_{i=1}^{25} \mathbf{S}\_{i} \frac{\Delta t}{\mathbf{N}\_{\text{PES},j}} \mathbf{N}\_{\text{PVF},i,j,n} \right] \tag{20}$$

where ( , ) *j j x y* is position the each sensor (the number of sensors is 6), *NPES*, *<sup>j</sup>* =1000 is the number of particles arriving at each sensor. *NPVS i*,,,*<sup>j</sup> <sup>n</sup>* is the quantity of fictitious particles in the *i-*th volume source, at the *n-*th instant.

1 5 6 10 11 15 16 20 21 22 23 24 25 0 *SSSS S S S S S S S S S AAAA A A A A A A A A A*

The dispersion problem is modelled by the LAMBDA backward-time code, where data from the Copenhagen experiment are used, for the period (12:33 h – 12:53 h, 19/October/1978) (Gryning & Lyck, 1984, 1998). It is assumed a logarithmic vertical

\*

above the ground, and *z*0=0.06 m the roughness. The wind speed was measured at 10 m, 120 m, and 200 m. The numerical value for *u*\* was obtained from the best fitting with the

Sensor-number *x* (m) *y* (m)

is *h*=1120 m (Gryning & Lyck, 1998). The turbulence parameterization follows Degrazia et al.

( *Li T* ). These two turbulent parameters are considered constant for the whole boundary layer, and their numerical values were calculate at *z*=10 m level. This characterises a

The sensors dimension, where the fictitious particles have arrived, was de 0,1 m 0,1 m 0,1 m, centred in the computational cells presented in Table 1. Next, the reverse trajectories are calculated for 1000 particles per sensor. The parameters for numerical simulations were 1800 time-steps with *t* = 1 s, meaning 30 min for the whole simulation. After 10 min, the concentration was computed at each 2 min, for the remaining 20 min of simulation. The

10 25

*<sup>t</sup> C xy S N*

1 1 ,

*j j i PVF i j n n i PES j*

*N*

where ( , ) *j j x y* is position the each sensor (the number of sensors is 6), *NPES*, *<sup>j</sup>* =1000 is the number of particles arriving at each sensor. *NPVS i*,,,*<sup>j</sup> <sup>n</sup>* is the quantity of fictitious particles in

1 400 500 2 600 300 3 800 700 4 1000 500 5 1200 300 6 1400 700

( ) ln *u z U z z*

0

(19)

the von Kármán constant, *z* the height

=180°, and the boundary layer height

*<sup>i</sup>* ) and the Lagrangian decorrelation time scale

,,,

(20)

profile for the wind field

being *U*(*z*) the main stream, *u*\* the fiction velocity,

Table 1. Position of sensors in physical domain.

(2000), for computing the wind variance ( <sup>2</sup>

For this experiment the wind direction had an angle

stationary and vertically homogeneous atmosphere.

mean concentration is found using the following expression:

*Mod*

the *i-*th volume source, at the *n-*th instant.

<sup>1</sup> (,) <sup>10</sup>

measured wind speed – see Table 2 – and the equation (19).

Table 2. Measured wind speed for the Copenhagen experiment (Gryning & Lyck, 1998).

The inverse approach is tested using synthetic measured concentration data The synthetic observations are emulated as:

$$\mathbf{C}^{Exp}(\mathbf{x}\_{j'}y\_j) = \mathbf{C}^{Mod}(\mathbf{x}\_{j'}y\_j)[1 + \lambda\mu] \tag{21}$$

where *CMod*(*x*,*y*) is the gas concentration computed from the forward model (14), is a random number associated to the Gaussian distribution with zero mean and unitary standard deviation, and is the level of the noise (Gaussian white noise). In our tests, =0.05 (Experiment-1) and =0.1 (experiment-2) were used. In other words, two numerical experiments are performed, using two different maximum noise levels: 5% and 10%.

The optimization method is used to find the minimum for the functional:

$$F(\boldsymbol{\alpha}, \mathbf{S}) = \sum\_{j=1}^{6} \left[ \mathbf{C}\_{j}^{\text{Eip}} - \mathbf{C}\_{j}^{\text{Mod}} \left( \mathbf{S} \right) \right]^{2} + a \, \Omega(\mathbf{S}) \tag{22}$$

where the vector 1 2 <sup>25</sup> [ , ,..., ]*<sup>T</sup>* **<sup>S</sup>** *SS S AA A* represent the emission rate from 25 sub-domains, presented in Figure 3. However, some sub-domains are considered releasing nothing. Therefore, the number of the unknown parameters to be estimated by inverse method decreases to 12 cells: 2 3 4 7 8 9 12 13 14 17 18 19 [,,,,,, , , , , , ] **S** *SSSSSSS S S S S S AAAAAAA A A A A A* . Another *a priori* information is smoothness of the inverse solution (regularization operator). As mentioned before, the second order maximum entropy is applied, Eq. (13c). The use of entropic regularization for this type of inverse problem was independently proposed by Roberti (2005) and Bocquet (2005a, 2005b) – see also Roberti et al. (2007) and Davoine & Bocquet (2007).

There are several methods to compute the regularization parameter. Here, a numerical experimentation was employed to determine this parameter. Roberti (2005) had found the values = 10-6 for Experiment-1, and = 10-5 for Experiment-2.