I. INTRODUCTION

One of the biggest problems of large and medium cities is the pollution of the surface atmospheric layer and soils by harmful emissions from vehicles. Large amount of harmful substances is concentrated in the vehicles exhaust gases. Among them, with high concentrations: Nitrogen oxides, Carbon oxides and Sulfur oxides. Motor transport is a source of atmospheric and soils pollution, which should be considered as distributed one. For reflecting and predicting of concentrations of harmful vehicle emissions, it is expedient to use mathematical models. They can be built based on the results of selective observations of their dynamics with known boundary errors of measurements. This approach was considered in [1].

The process of harmful emissions distribution and its dynamics is considered as a mass transfer process. For its description, the difference operators (schemes) are used. Their identification is carried out using the data of measurements of harmful emissions concentrations with known boundary errors. Such data are called interval data [2,3]. As is known, the methods of difference operators identification based on the interval data analysis require a uniform measurement grid that is impossible for the real city conditions. Mostly, the measurements of harmful emission concentrations is carried out in places with intensive traffic and accumulation of vehicles. This means that measurement grid is not uniform. Thus, for building of mentioned models, it is necessary to solve three tasks related to data preparation: execute cluster analysis for defining of homogeneous vehicular traffic intensity areas; identify the discrete values of the grid step; calculate estimates of harmful vehicle emissions concentrations in the nodes of the grid. The third task is solved by methods of interpolation [4,5]. The first one and second one are the subjects for research of this work.

II. STATEMENT OF THE PROBLEM

To solve the environmental monitoring tasks, it is necessary to build models of stationary and non-stationary fields of concentrations of harmful vehicle emissions [6]. The theoretical basis for solving this type of tasks are the mathematical models of objects with distributed parameters in the form of partial differential equations. Concentration of attention on the physical properties of environment requires to significantly complicate the mathematical model. Even despite the fact that, in practice, it is impossible to verify the results of modeling with real data obtained under conditions that meet the conditions of modeling. First of all, this is related to the complexity of the measurement experiment. For example, if a mathematical model in the form of a differential equation accurately enough describes the process of transferring of chemical substances in the atmosphere in case of wind gusts or other turbulent phenomena in the atmosphere, then an integrated value of the chemical substance concentration per volume unit is established in the process of measurement. In addition, the accuracy of such measurements is low, the relative measurement error may reach 50%. Consequently, it is enough to build a mathematical model with an accuracy that is equiualent to the accuracy of the measurement experiment. At this, it is expedient to represent the experimental data in the form of intervals of possible values of the modeled characteristic: , , , ,

1,..., , 1,..., , 1,..., , 1,...,

i j h k i j h k z z i I j J h H k K − + = = = =(1)

where , , , , , are the lower and upper bounds of the interval of possible values of measured concentration of harmful substances in the grid nodes with discretely given spatial coordinates

+ − k h j i k h j i z z , ,I i ,..., 1 = , J j ,..., 1 = , H h ,..., 1 = at the dicrete time value K k ,..., 1 =

, respectively. It is worth to note that, in the measurement experiment, we can set the lower and upper bounds based on the relative error of the measuring device: , , , , , , , , where

ε ⋅ − = − k h j i k h j i k h j i z z z , , ,, , , , , , and ε⋅ + = − k h j i k h j i k h j i z z z , ,k h j i z , , ,

is the measured value of the harmful substance concentration; ε is relative measurement error.

Under these conditions, macromodeling is the only way to reflect the distribution of harmful emissions concentrations. The building of such macromodels is convenient to carry out based on the obtained interval data in the form (1).

In the papers of O.G. Ivakhnenko [7], the inductive approach is described for choosing of acceptable way of mathematical description of these processes. Its essence consists in defining of some difference scheme in the way of its adjustment in accordance with the experimental data. The difference scheme itself, which actually converts the values of input variables into output values, is called a difference operator. The process of adjustment of this scheme is called structural identification [8,12].

In general case, the expression of linear in parameters difference operator (DO) has the following form [2]: , , , 0,0,0,0 0,0, 1,0 1,0,0,0 0, 1,0,0

T i j h k h i j i j h k i j h i j h k v f v v v v v u u g i d j d J h d H k d K − − − − − − − = ⋅ = = = =    (2)

where

) (• T f 

is the vector of basis functions (nonlinear, in general case) by which, the transformation of the modeled characteristic values, as well as the input variables in the spatial grid nodes for the certain discrete moments of time is carried out;

k h j i v , , ,

modeled concentration of harmful emissions in grid nodes with discretely-given spatial coordinates

I d i ,..., = , J d j ,..., = , H d h ,..., = at the moments of time K d k ,..., = ; k h j i h j i u u , ,, 0 , , , ,...,  

are the vectors of input variables (controls); d is the DO order; g  is the vector of unknown parameters of DO.

As a result of executing of structural identification procedure, we establish the difference computational scheme, in particular: the basis functions vector

) (• T f 

; sets and dimensionality of input variables (controls) vectors

k h j i h j i u u , ,, 0 , , , ,..

.,  

; order of the difference scheme d, which, as is known, is equivalent to the order of the differential equation (the analogue of the difference scheme). To realize the difference scheme, it is also necessary to set the initial conditions, that is, the value of each discrete element from the set ,..., ,..., , ,...,

0 , 0 , 1 , 0 0 , 0 , 0 , 1 0 , 1 , 0 , 0 0 , 0 , 0 , 0 − − − j i h v v v v , 1 , 1 , 1 , 1 − − − − k h j i v k h j i h j i u u , , , 0 , , , ,...,  

(as a rule, the initial ones) and to establish the values of the parameters vector g  components. If the structure of DO is known then, it remains the actual task of adjusting the parameters of DO (2) in such a way as to ensure the maximal consistency between the modeled characteristic and experimentally obtained values of this characteristic. Such a task is called the parametric identification task [9,13].

Based on the requirements of ensuring the mathematical model accuracy within the bounds of the measuring experiment accuracy, the conditions of consistency between experimental data, represented in the interval form (1), and data obtained based on the mathematical model in the form of DO (2), can be formulated in such form: , , , ,

,..., , 1,..., , ,..., , ,...,

i j h k i j h k i j h k v z z i d I j J h d H k d K − + ⊂ ∀ = ∀ = ∀ = ∀ = (3)

Based on the results of conducted analysis, we can state that for ensuring of conditions of the given accuracy (3) of the macromodel in the form of linear DO (2) during solving the task of its parametric identification, the application of interval data analysis methods [9] is substantiated.

Let's assume that the vector of parameters estimates g   in the DO ( 2) is obtained based on the interval data analysis. Substituting the vector of parameters estimates g   of DO instead of their true values g  in expression (2) together with given initial interval values of each element in the set ],...

, [ ],..., [ ], [ ],..., [ 0 , 0 , 1 , 0 0 , 0 , 0 , 1 0 , 1 , 0 , 0 0 , 0 , 0 , 0 − − − j i h v v v v     ] [ 1 , 1 , 1 , 1 − − − − k h j i v 

and given vectors of the input variables

k h j i h j i u u , ,, 0 , , , ,...

, T i j h k i j h k i j h k hi j i j h k i j h i j h k v v v f v v v v v u u g i I j J h H k K − + − − − − − − − = = ⋅ = = = =             (4)

Thus, the mathematical model of stationary and nonstationary fields of harmful emissions concentrations for the task of environmental control will be described by a DO in general form (4). Taking into account that all calculations in equation ( 4) are carried out using interval arithmetic rules [2], the difference operator (4) is called an interval difference operator (IDO).

The conditions of consistency of experimental data, represented in interval form (1), with the data obtained based on macromodel in the form of IDO (4) are formulated as follows:

, , , ,

1,..., , 1,..., , 1,..., , 1,...,

i j h k i j h k i j h k i j h k v v z z i I j J h H k K − + − + ⊂ ∀ = ∀ = ∀ = ∀ =  (5)

Let's substitute in expressions (5), instead of interval estimates of the harmful substance concentrations ] ; [ , , , , , ,

+ − k h j i k h j i v v  

, its interval values calculated using IDO (4), together with taking into account the given initial interval values of each element from the set 0,0,0,0 0,0,0,0 0,0, 1,0 0,0, 1,0 1,0,0,0 1,0,0,0 0, 1,0,0 0, 1,0,0

1, 1, 1, 1 1, 1, 1, 1 [ ] [ ],...,[ ] [ ], [ ] [ ],...,[ ] [ ],..., [ ] [ ] h h i i j j i j h k i j h k v z v z v z v z v z − − − − − − − − − − − − − − ⊆ ⊆ ⊆ ⊆ ⊆      (6)

and given vectors of input variables

k h j i h j i u u , , , 0 , , , ,...,  

. We obtain such interval system of non-linear algebraic equations (ISNAE) [3]: 0,0,0,0 0,0,0,0 0,0,0,0 0,0,0

, 1 0,0,0,0 0,0, 1,0 1,0,0,0

., i d j d h d k d i d j d h d k d i d j d h d k d i d j d h d k d T i j h kh i v v z z v v z z v f v v v − + − + − + − + − − − − − − − − − − − − − − − − − − − − − − ⊆ ⊆ =          0, 1,0,0 , , ,0 1

, , , 0,0,0,0 0,0, 1,0 1,0,0,0 0, 1,0,0 , , , 0 ,..., , 2,..

j i d j d h d k d k T i j h k h i j i d j d h d k d k i j h k v v u u g z f v v v v v u u g z i d I d − − − − − − − − − − + − − − − ⋅ ≤ ⋅ ≤ = =                 ., ,

,..., , ,..., .

J h d H k d K              = =   (7)

So, the ISNAE ( 7) is obtained by substituting the interval estimates of the output characteristic (given in the form of initial conditions and predicted using expression (4) in the remaining nodes of the grid) into conditions (5). Therefore, the task of parametric identification of IDO (4) under conditions ( 5) is the task of solving ISNAE in the form (7). Methods for estimation of solutions of the obtained ISNAE are described in [10].

The analysis of the proposed scheme for building of mathematical model of harmful vehicle emissions distribution showed that before its implementing, it is necessary to obtain a uniform grid of measured concentrations (3) in its nodes and vectors of influences on them

k h j i h j i u u , , , 0 , , , ,...,  

. The main among them, is the vehicular traffic intensity. This task is solved using modified method of subtractive clustering of data [11] on the traffic intensity.

III. MODIFIED SUBTRACTIVE CLUSTERING METHOD

As the basis for method of clustering of vehicular traffic distribution, it is expedient to use the "mountain" clustering method with subtractive algorithm of its implementation. This method does not require a large sample of experimental data and does not require to set a predetermined number of clusters that significantly reduces the time for its implementation. It is also worth to note that the number of clusters based on this method is regulated by the only parameter which is the cluster radius [11]. According to the clustering method, in the beginning, we form the potential cluster centers from the rows of data matrix for the clustering of input variables and calculate the potentials of identified cluster centers using the expression:

( ) ( ) 1 exp K h h h k k P c c x α = = − ⋅ − ∑   ,(8)

where ) ,..., , (

2 1 Kh h h h c c c c =  is a potential center of h-th cluster; а is a positive constant; k h x c   − is a distance between potential center of h-th cluster h c  and input experimental data k x  , k=1,…,K, H h ,..., 1 =

; H is a number of possible clusters.

In our case, if the only property of a cluster which is the number of vehicles The illustration of the potentials distribution is represented as a surface in the form of a mountainous relief (Fig. 1), whose peaks have the highest potential values and are pretenders to be the centers of the formed clusters. As we can see in the Fig. 1, one "mountain peak" is surrounded by other peaks that causes the problem of building of very similar data clusters with the corresponding centers. This does not provide the high quality clustering results. As centers of the clusters, we choose the coordinates of "mountain peaks", that is, the center of the cluster is the point on the city map with the highest value of potential:

, , exph h i j I J h h h x y k x y k i j P x y k u u α = = = − ⋅ − ∑ ∑ ,(9) where ( ) ,

( )

= . (10)

In order to avoid the building of similar clusters, we must recalculate the potential values for the remaining possible cluster centers: − is a distance between potential center of h+1 cluster and center of found h-th cluster. The procedure of cluster centers calculation is carried out until all the rows of the input variable matrix X, which is represented by the set (3), are excluded.

( ) ( ) ( ) 1 1 1 1 1 1 , ,, , , , , , , , e , ) 1 ( ,

The above procedure is based on the subtractive clustering algorithm, which is based on the following steps.

Step 1. Forming of potential cluster centers. They are all points of measured harmful emission concentrations and corresponding intensities of vehicular traffic.

Step 2. Calculation the potential of possible cluster centers based on (9).

Step 3. Selecting the data point with the maximal potential for representation of the cluster center based on (10).

Step 4. Excluding the influence of the found cluster center in the way of recalculating the potentials for other possible cluster centers by (11).

Step 5. Identifying the next cluster and the coordinates of its center. If the maximal value of the cluster center potential exceeds some predetermined threshold which is the cluster radius, that is

( ) , ,

h h h P x y k , then proceed to step 4, otherwise, the algorithm is completed. The iterative procedure for identification of cluster centers and the recalculation of potentials is repeated until all points in the space of input experimental data are located within the neighborhoods of the radius of sought cluster centers.

As a result of the clustering algorithm implementation, we obtain h clusters, h = 1,…,Н, with the corresponding centers. The next step is the identification of a uniform grid nodes for homogeneous parts of a cluster. The discrete values of the grid nodes coordinates are equal to the cluster diameter, and the value , ,

IV. EXAMPLE OF CLUSTERING METHOD APPLICATION

Let's consider the application of the developed clustering method for obtaining of uniform grid of nodes on an example of Ternopil city.

The fragment of map of central part of Ternopil city with the marked points of the measured vehicular traffic intensity for one discrete time (one hour) is shown in the Fig. 2. As we can see, the traffic is distributed not uniformly over the territory. Therefore, it is advisable to measure its intensity at some selected points, where this intensity is the highest, for example, as it is shown on the map of Ternopil city.

The points of measurement of traffic intensity are colored red on the map.

The application of cluster analysis for determination of areas with specific vehicular traffic intensities under condition of identification of cluster centers that are located on a certain uniform grid gives a possibility to define the discretization step for building a DO. In our case, the cluster is the set of points of a certain area of the city with similar values of current vehicles number. The result of the proposed method of cluster analysis application is schematically shown in Fig. 3. As we can see, during the clustering process, H clusters with different vehicular traffic intensity and radius r were defined. So, the discrete values of grid nodes coordinates are equal to the cluster diameter. Obtained grid for building of distribution model of harmful vehicle emission concentrations in the form of IDO ( 4) is schematically shown in Fig. 4. For building of mentioned model in the form of IDO (4), it is enough to execute interpolation and identify the pollution concentrations in the grid nodes.

V. CONCLUSIONS

The modified method of subtractive clustering and interval analysis for modeling of distribution of harmful vehicle emissions concentrations and vehicular traffic intensity under conditions of non-uniform measurement grid were proposed and substantiated.