F Kuznik

J Virgone

Experimental Investigation of Wallboard Containing Phase Change Material: Data for Validation of Numerical Modeling

Keywords: experimental data, PCM composite, phase change material, wallboard, thermal energy storage, hysteresis

In construction, the use of Phase Change Materials (PCM) allows the storage/release of energy from the solar radiation and/or internal loads. The application of such materials for lightweight construction (e.g. a wood house) makes it possible to improve thermal comfort in summer and reduce heating energy consumption in winter. The choice of a PCM depends deeply on the building structure, on the weather and on building use: numerical modeling is indispensable. In this paper, an experimental comparative study is described, using cubical test cells with and without PCM composite. A set of experimental data is detailed, concerning the air and wall temperatures. The results are compared with a numerical modeling and show that hysteresis must be taken into account to predict correctly the heat transfer.

Introduction

Nowadays, thermal energy storage systems are essential for reducing dependency on fossil fuels and then contributing to a more efficient environmentally benign energy use [START_REF] Dincer | Thermal Energy Storage -Systems and Applications[END_REF]. As demand in thermal comfort of buildings rises increasingly, the energy consumption is correspondingly increasing. For example, in France, the energy consumption of buildings has increased by 30% the last 30 years. Housing and tertiary buildings are responsible for the consumption of approximatively 46% of all energies and approximatively 19% of the total CO 2 emissions [START_REF]Climate Plan 2004: let's act together to challenge of climate change[END_REF].

Thermal energy storage can be accomplished either by using sensible heat storage or latent heat storage. Sensible heat storage has been used for centuries by builders to store/release passively thermal energy, but a much larger volume of material is required to store the same amount of energy in comparison to latent heat storage. The principle of the phase change material (PCM) use is simple. As the temperature increases, the material changes phase from solid to liquid. The reaction being endothermic, the PCM absorbs heat. Similarly, when the temperature decreases, the material changes phase from liquid to solid. The reaction being exothermic, the PCM desorbs heat.

The main disadvantage of light weight buildings concerning thermal comfort and energy consumption is their low thermal energy storage potential in walls.

Obviously, they tend to large temperature fluctuations due to external cooling or heating loads. Using PCM material in such building walls can decrease the temperature fluctuations, particularly in case of solar radiations loads. It is then a potential method for reducing energy consumption in passively designed buildings. This tendency is confirmed by numerous papers. For a review, see in [START_REF] Tyagi | PCM Thermal Storage in Buildings: A State of Art[END_REF].

When selecting a PCM, the average room temperature should be close to the melting/freezing range of the material. Moreover, the day fluctuations should allow the material phase change. Then many factors influence the choice of the PCM: weather, building structure and thermophysical properties,... Numerical simulation must then be used to achieve the practical application of this technology. In order to validate a PCM numerical model, experimental data are essential.

The table 1 summarizes, not exhaustively, some experimental studies concerning measurements held in a room with walls containing PCM wallboards.

Most of these experiments were carried out in outdoor conditions. Even if these conditions are more realistic, all of the weather parameters are not known (or measured). These studies cannot then be used to validate a numerical model requiring all of the experimental conditions. Concerning the number of the test cells, few studies use two identical cells with and without PCM.

Design of the cells

Two identical test cells are used in the study. Each test cell is a cubical enclosure with an internal dimension of 0.50m (figure 1). The skeleton of MICRO-BAT is made of polish aluminium and is not painted in order to minimize the shortwave radiative heat transfer with the surroundings.

The front face of MICROBAT is an aluminium plate of 2mm thickness. This face is called the active face because of its low thermal inertia and low thermal resistance. It enhances the heat transfer between the exterior and the interior of MICROBAT. The other faces of the test cells are of two types: =⇒ the normal face (figure 2), mainly composed of insulating material in order to have a low thermal inertia wall, =⇒ the PCM face (figure 2), which contains 5mm of a phase change material wallboard.

The witness box has five normal faces whereas the PCM box contains three PCM faces which are the back face, the right face and the left face.

The thermophysical properties of the materials are given in table 2. All of the internal faces are painted in black, the absorption coefficient in short-wave being measured at about 0.96.

Phase change material tested

The product tested, EN ERGAIN has been achieved by the Dupont de Nemours Society and is constituted of 60% of microencapsulated paraffin within a copolymer. The final form of the composite PCM material (see fig- The thickness of the PCM is the result of an optimization study that can be found in [START_REF] Kuznik | Optimization of a phase change material wallboard for building use[END_REF].

The composite PCM heat capacity has been measured using a differential scanning calorimeter; the heating and cooling rate being 0.05Kmn -1 . This rate corresponds to 3 • C/h which is an average heating rate in a light weight building during summer season and when solar gains are maximum. The thermal analysis is given for the range [1 • C; 34 • C]. Two curves are presented: the freezing curve (cooling from 34 • C to 1 • C) and the melting curve (heating from 1 • C to 34 • C). Figure 4 shows the measured values of heat capacity.

The nearly Gaussian curves correspond to an intermediate phase change for the material [START_REF] Darkwa | Simulation of Phase Change Drywalls in a Passive Solar Building[END_REF]. Paraffin PCMs display a wider melting range compared to inorganic PCMs (narrow phase change zone). The physical properties of the composite PCM are:

⋆ peak melting temperature = 22.3 • C, ⋆ peak freezing temperature = 17.8 • C,
⋆ specific heat at the melting peak = 13.4Jg -1 K -1 , ⋆ specific heat at the freezing peak = 12.9Jg

-1 K -1 , ⋆ solid specific heat = 2.4Jg -1 K -1 .
The latent heats of freezing and melting are respectively 71.0J/g and 72.4J/g. The difference between the melting temperature and the freezing temperature characterizes the hysteresis of the material, i.e. the fact that the mixture is not an eutectic. The figure 5 shows the measured enthalpy of the composite PCM function of temperature. This curve is typical of the hysteresis loop.

Concerning the numerical modeling of this phenomenon, there is a memory effect and then the order of previous events influences the order of subsequent events.

The thermal conductivity has been measured using guarded hot-plate apparatus [START_REF]Standard Test Method for Steady-State Heat Flux Measurements and Thermal Transmission Properties by Means of the Guarded-Hot-Plate Apparatus[END_REF]. The thermal conductivity is 0.22W.m -1 .K -1 in liquid phase (θ > 25 • C) and decreases to about 0.18W.m -1 .K -1 in solid phase (θ < 5 • C).

Instrumentation and measurements

All of the temperatures are measured using Pt100 sensors with a calculated resolution of ±0.25 • C. The exterior temperature is measured. In each test cell the set of probes is:

• the internal face temperature, measured at the center of the face, with a probe included in the wall surface,

• the active face external temperature, measured at the center of the face,

• the air temperature, measured at the center of the cubical air volume.

A total of 17 probes is used in the experiment. The acquisition is done by the means of a multiplexer-multimeter connected to a personal computer. The time step chosen between two series of measurement is 2mn and the duration of each test is about three days.

Experimental protocol

The two MICROBAT test cells are placed into a climatic chamber. The climatic chamber temperature can vary between -10 • C and 40 • C as a function of time. The active face allows to impose the same external temperature for each test cell. The air in the climatic chamber is mixed in order to avoid temperature gradients; the mean air velocity in this chamber is about 0.5m/s. Two types of external temperature evolutions are investigated:

• a temperature step, heating and cooling,

• a sinusoidal temperature evolution.

Concerning the heating/cooling temperature step, three cases are tested depending on the slope (SL) i.e. the time needed to reach the constant temperature value: 1hour, 2hours or 3hours. The figure 6 shows the temperature step parameters. A time delay of 20hours is needed for the two cells to reach the same stationary thermal state.

The sinusoidal temperature evolution is 24hours periodic and varies between 15 • C and 30 • C (figure 6).

Temperature step

The first experiment concerns the temperature steps (heating and cooling) and the present part deals with the results of the cases tested. The part 3.1 deals with general considerations concerning the measurements. The air temperatures is analyzed part 3.2 and the internal face temperature part 3.3.

General considerations

The box without PCM is strongly insulated but with a relatively low thermal inertia. The figure 7 shows the measured temperatures for the internal faces and air in the test cell without PCM for the case SL = 1hour. Except for the active face, the other faces temperatures and the air temperature are nearly identical. The temperature gradient in the box does not affect the heat transfer, and particularly the convection heat transfer. This conclusion is valid for all the cases tested and then only air temperature is presented in the following paragraphs of the article.

Air temperature analysis 3.2.1 Heating step

The air temperature measured for the three kinds of heating steps and for the two test cells is presented in figure 8. The air temperature for the box without PCM is close to the exterior temperature with a little time lag due to the low thermal inertia of the box. The storage effect of PCM affects much the air temperature.

The air temperature evolution of the witness box is nearly exponential. For the box with PCM, the temperature evolution is composed of three stages:

(1) for θ(t) 21 • C, the air temperature is increasing linearly with a slope SL1 (definition of SL h is given figure 6),

(2) for 21 • C θ(t) 30 • C, the air temperature is increasing with a slope lower than SL h due to the phase change,

(3) for θ(t) ≈ 30 • C, the temperature has reached a steady state.

In order to evaluate the time lag due to the PCM, the τ value is introduced and defined as:

θ (2τ ) = θ (0) + 1 -e -2 (θ (∞) -θ (0)) (1)
where θ (t) is the air temperature at time t. For an exponential temperature evolution, τ is the time needed to reach 63.2% of the infinite temperature value.

The table 3 summarizes the τ values for the heating step cases. The τ ratio between the case with PCM and without PCM varies between 3.9 and 2.4:

the more the slope of the external temperature step is important, the more the ratio is important. Consequently, the more the heating excitation is rapid, the more the PCM wallboard is efficient for slowing down the temperature increase.

Cooling step

The air temperature measured for the three kinds of cooling steps and for the two test cells is presented in figure 9. The conclusions from this figure are quite similar to those of the heating step i.e.: the air temperature for the box without PCM is close to the exterior temperature with a little time lag due to the low thermal inertia of the box and the storage effect of PCM affects much the air temperature.

The air temperature evolution of the witness box is nearly exponential. For the box with PCM, the temperature evolution is composed of three stages:

(1) for θ(t) 19 • C, the air temperature is decreasing linearly with a slope The main difference between the heating and the cooling is the flat part of the curve around 19 • C appearing only for the cooling step cases. This flat part of the curve lasts about 1h. The PCM composite is composed of microencapsulated paraffin spheres included in a copolymer matrix. During the solidification process, the solid paraffin is evenly formed through the sphere, starting from the outer surface and moving inward. As the solidification proceeds, the melt volume decreases with a simultaneous decrease in the magnitude of natural convection within the melt and the process is therefore much longer [START_REF] Ettouney | Heat transfer during phase change of paraffin wax stored in spherical shells[END_REF].

SL c , (2)
The table 3 summarizes the τ values for the cooling step cases. The τ ratio between the case with PCM and without PCM varies between 5.2 and 3.0: the more the slope of the external temperature step is important, the more the ratio is important. Then, the more the cooling excitation is rapid, the more the PCM wallboard is efficient for slowing down the temperature decrease.

Discussion

The response of MICROBAT to an external temperature step is interesting to characterize the time lag due to the PCM. Of course, the presence of PCM allows to delay the box air temperature increase or decrease (as long as the temperature varies within the phase change range).

The analysis of the time lag defined paragraph 3.2.1 shows that the more the thermal excitation is rapid, the more the PCM is efficient. In buildings, rapid thermal excitation can be solar spot or internal load.

For the PCM used in wallboards, hysteresis exists and has clearly been exhibited in paragraph 2.2. This phenomenon has been few studied in literature and never been taken into account in numerical modeling as well as in experimental analysis. In the present study, an optimal use of the PCM requires a wall temperature above 22. The response of the MICROBAT to an external temperature step allows to characterize the time lag due to PCM wallboards. In order to characterize the phase difference and the decrement factor due to PCM, an external sinusoidal temperature evolution is used, case which is closer to the real building configuration.

The figure 12 shows the air temperature curves for the two cells. The phase difference ζ is defined as the time difference between outdoor temperature maximum and indoor temperature maximum. For the cell with PCM, the phase difference ζ = 138min while this value is about 38min without PCM.

The decrement factor f is defined as the ratio between indoor temperature amplitude and outdoor temperature amplitude. For the cell with PCM, the decrement factor f = 0.74 while this value is about 0.89 without PCM.

The figure 13 presents the wall surface temperatures for the PCM box. For the walls with PCM, the temperature curve has two breaks of slope: one at about 19 • C and one at about 22 • C. This feature is a characteristic of the hysteresis phenomenon of the PCM.

Numerical modeling

In the previous section, experimental data are described so that in the present section numerical modeling can be held. The paragraph 5.1 is dedicated to a presentation of the numerical modeling. Then, the numerical results are presented in paragraph 5.2. The paragraph 5.3 is a discussion concerning the numerical modeling.

Presentation of the numerical modeling

The first modeling step consists in decomposing the test cell (or the building) in elementary objects: air zones, walls,.. The temperature evolution in the test cell is calculated using the energy balance of the air volume.

For each wall, the conduction heat transfer is supposed to be unidirectional.

The phase change is modeled using the equivalent heat capacity method: the exact heat capacity curve is approximated using splines. In order to solve numerically the problem, a finite-difference method is used: the continuous information contained in the exact solution of the differential equation are replaced by discrete temperature values. More information about such numerical scheme can be found in [START_REF] Kuznik | Energetic Efficiency of Room Wall Containing PCM Wallboard: A Full-Scale Experimental Investigation[END_REF].

Only the longwave radiation exchanges are considered in the present modeling.

The external exchanges occur between the external surface of the walls and the climatic chamber walls. The internal exchanges occur between the internal surfaces of the walls. All the longwave radiation exchanges are fully calculated, using the form factor of the various surfaces.

All of the temporal problem equations are discretized using an implicit scheme.

The time selected chose for the simulation is 60s. For the conduction problem, the spatial discretization results in 10 PCM nodes and 20 insulating material nodes. All the simulation parameters are optimized from a solution accuracy point of view. The software MATLAB is used for all the simulations.

Numerical results

In this paragraph, the experimental data are compared with the numerical results obtained from our modeling. For the sake of brevity, only the results concerning the air temperature are presented in this paper.

The figure 14 The figure 16 shows the comparisons between experiments and modeling for the sinusoidal temperature evolution. From a modeling point of view, this type of evolution is complicated because it is composed of cooling and heating periods. Neither the freezing curve, nor the melting curve allows to predict the correct air temperature evolution.

Discussion concerning the numerical modeling

The main problem of the PCM modeling is the way to introduce the phase change. In this paper, the equivalent specific heat has been tested. Unfortunately, the paraffin used is not an eutectic mixture, so hysteresis occurs during the phase change. When simulating heating step or cooling step separately, the numerical data are in good agreement with experiments only if the corresponding specific heat curve is used, e.g. melting curve and cooling curve.

For a sinusoidal temperature evolution, neither the melting curve nor the freezing curve allows to predict correctly the temperature evolution. The hysteresis phenomenon must the be taken into account correctly in order to predict the PCM composite thermal behavior.

Conclusions and outlook

The objective of this article is first, providing reliable experimental data that can be used for the validation of numerical modeling and then, studying some features related to the use of phase change material wallboard. A double test cell, MICROBAT, is used to measure the temperatures and characterize the phase change effects. The external temperature is the thermal excitation; a heating/cooling step and a sinusoidal evolution are tested.

The effects of PCM wallboard are to cause time lag between indoor and outdoor temperature evolutions and to reduce the external temperature amplitude in the cell. A comparison between the cases with and without PCM shows that the more the thermal excitation is rapid, the more the PCM is efficient.

This enables the use of the phase change material for solar spot in building rooms.

The experimental data are then compared with numerical simulation. The heating steps and cooling steps are correctly predicted by modeling when the corresponding specific heat curve is used. For a sinusoidal temperature evolution, neither the melting curve nor the freezing curve allows to predict correctly the temperature evolution.

The effect of hysteresis phenomenon has been clearly exhibited with the experimental data: the melting process arises at a temperature higher than for the solidification process. The numerical modeling shows that the hysteresis phenomenon must be taken into account correctly in order to predict the PCM composite thermal behavior. Further investigations are needed to have a better numerical description of this special phase change feature.

NOMENCLATURE f decrement factor [-] θ temperature [ • C] SL temperature step

Heating

Step 1h 48min 188min 3.9

Step 2h 70min 207min 3.0

Step 3h 93min 221min 2.4

Cooling

Step 1h 48min 249min 5.2

Step 2h 70min 262min 3.7

Step 3h 93min 280min 3.0 Table 3 Time lag values for for the various temperature steps -heating and cooling.

  From a numerical point of view, having experiments with and without PCM is useful to verify the basic modeling. The part 2 of this article presents the two identical test cells called MICRO-BAT and the PCM tested. Two kinds of external temperature evolutions are tested: heating and cooling steps with various slopes and sinusoidal temper-ature evolution with 24h period. The part 3 deals with the results of the temperature steps. The sinusoidal external temperature results are presented in the next part of the paper. The part 5 is dedicated to the numerical modeling of the experiments: the main assumptions of the model are presented, then the simulations are compared with experimental data and finally, the numerical modeling is discussed. 2 Description of the test cells MICROBAT Two identical test cells are used to investigate the effects of the PCM wallboards. The design of the test cells is presented part 2.1. The phase change material tested is described part 2.2. The part 2.3 deals with the measurement systems and the probes, the experimental protocol is given part 2.4.

ure 3 )

 3 is a flexible sheet of 5mm thickness which density is about 900kg.m -3 .

  for 15 • C θ(t) 19 • C, the air temperature is decreasing with a slope lower than SL c due to the phase change, (3) for θ(t) ≈ 15 • C, the temperature has reached a steady state.

3. 3

 3 figures 10 and 11 present the internal faces temperature , for the cell with PCM, and for the heating and cooling steps.

3 •

 3 C for energy storage process, and a wall temperature below 17.8 • C for energy release process. From an engineering point of view, this phenomenon must be taken into account! 4 External sinusoidal temperature evolution

  . The air zones are connected via walls. Each elementary object has homogeneous physical values: one surface temperature per wall, one temperature per air zone,...The problem treated in this article is composed of two air zones (inside and outside the test cell) and 13

  shows the comparisons between experiments and modeling for the heating step. The modelings (a), (b) and (c) use the equivalent specific heat obtained with the data of the figure 4 melting curve. For the three cases, the numerical data are in good agreement with experiment. The figure (d) shows that the freezing specific heat curve is not adapted for the heating step modeling. The figure 15 shows the comparisons between experiments and modeling for the cooling step. The modelings (a), (b) and (c) use the equivalent specific heat obtained with the data of the figure 4 freezing curve. For the three cases, the numerical data are in quite good agreement with experiment, excepted for the flat part of the curve around 19 • C. Due to the nature of the flat plate of the curve, this feature is not modeled. The figure (d) shows that the cooling specific heat curve is not adapted for the cooling step modeling.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Schematic view of a test cell.

Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 345678910 Fig. 3. Duoont de Nemours material.

Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .Fig. 15 .

 1112131415 Fig. 11. Temperatures for the boxes with PCM -external temperature decrease.

Fig. 16 .

 16 Fig. 16. Comparison between experimental data and numerical modeling -Sinusoidal temperature evolution.

  Box without PCM Box with PCM τ Ratio

			τ [min]
	Ref.		PCM composite	Size	Number of cells Conditions
	t	[4] [5] [6] [7]	slope [hour] butyl stearate / gypsum time [s] 2.82m, 2.22m, 2.24m τ polyethylene glycol / PVC panel 0.9m, 0.9m, 0.9m density specific heat thermal conductivity time lag [min] 1 2 paraffin / polymer 3.10m, 3.10m, 2.50m 1 [kg/m 3 ] [J/kgK] [W/mK] capriclauric acids / gypsum 5.00m, 3.30m, 2.80m 1 aluminium 2700 8800 230.00 insulating material 35 1210 0.04	outdoor outdoor laboratory outdoor
		[8]	paraffin / plaster cardboard 730	Subscripts 4.43m, 2.95m, 3.00m 1110 0.27	2	outdoor
	Greek Letters [9] Table 2	paraffin / gypsum	c 2.00m, 2.00m, 3.00m	cooling	2	outdoor
	ζ [10] Thermophysical properties of the materials. phase difference [min] paraffin / polyethylene 3.00m, 2.00m, 2.00m h	heating	1	outdoor
	Table 1
	Experimental studies involving PCM wallboards.

Acknowledgments

The authors thank Hébert Sallé from the Centre Scientifique et Technique du Bâtiment (CSTB) for the DSC analysis.


Dominique Béréziat
email: dominique.bereziat@upmc.fr

Isabelle Herlin
email: isabelle.herlin@inria.fr

SOLVING ILL-POSED PROBLEMS USING DATA ASSIMILATION. APPLICATION TO OPTICAL FLOW ESTIMATION

Keywords: Data Assimilation, Ill-posed Problem, Variational Formulation, Optical Flow, Missing Data

Data Assimilation is a mathematical framework used in environmental sciences to improve forecasts performed by meteorological, oceanographic or air quality simulation models. Data Assimilation techniques require the resolution of a system with three components: one describing the temporal evolution of a state vector, one coupling the observations and the state vector, and one defining the initial condition. In this article, we use this framework to study a class of ill-posed Image Processing problems, usually solved by spatial and temporal regularization techniques. A generic approach is defined to convert an ill-posed Image Processing problem in terms of a Data Assimilation system. This method is illustrated on the determination of optical flow from a sequence of images. The resulting software has two advantages: a quality criterion on input data is used for weighting their contribution in the computation of the solution and a dynamic model is proposed to ensure a significant temporal regularity on the solution.

Introduction

In the research field of Image Processing, most problems are ill-posed in the sense that it is not possible to provide a unique solution [START_REF] Hadamard | Lecture on Cauchy's Problem in Linear Partial Differential Equations[END_REF]. A common cause of ill-posedness is that the equations used to model image properties, called Image Model, are under-determined. An example is the famous "aperture problem" occurring in the estimation of optical flow.

A possible strategy to solve ill-posed problems is to provide the Image Model with additional information. Two options are considered: 1) Providing explicit information: additional images are used to enlarge the set of input data. However, this is generally not possible because other acquisitions having the requested properties are often not available. 2) Providing implicit information: either hypotheses on image properties or constraints on the solution can be used. A usual constraint is to restrict the dimension or the size of the space of admissible solutions: for instance the result is searched among the functions having bounded spatial variations [START_REF] Tikhonov | Regularization of incorrectly posed problems[END_REF]. Such methods are called "Tikhonov regularization methods". In the general case, these additional image properties and constraints are expressed as equations, which, added to the Image Model, lead to a new invertible model.

Assuming an image sequence is available, the Tikhonov regularization can be extended by considering a regularization acting both in spatial and temporal domains. This has been proposed by Weickert [START_REF] Weickert | Variational optic flow computation with a spatio-temporal smoothness constraint[END_REF] for determining optical flow. The interest is to produce a solution which is regular in time and in space. This approach has however several drawbacks. First, the temporal regularization imposes an arbitrary regularity which can not deal with complex dynamics. Second, if irrelevant input data values occur on large regions, they are however taken into account in the estimation process and introduce errors in the solution. These aberrant values may be due to a sensor failure or an acquisition noise. We call them "missing data" in the remain of this paper. Third, as the solution is looked for in the spatio-temporal domain, this leads to a high complexity in comparison with a pure spatial model, which is a limiting factor for operational purposes.

An alternative to the Weickert's approach is to model the temporal evolution of images. Such a modelling can be inferred, for instance, from a priori knowledge on image acquisition. The challenge becomes to define an effective and pertinent model, and to include it in the solution computation. Data quality has also to be evaluated in order to ignore missing data.

In this paper, we propose to use the Data Assimilation framework as a generic tool to solve ill-posed problems in Image Processing. A Data Assimilation method solves simultaneously an evolution equation, describing the evolution of a state vector over time, and an observation equation, modelling the links between the state vector and the observations. Each equation can be inaccurate and a description of the errors is stated in terms of a Gaussian noise characterized by a covariance matrix. As it will be shown in the paper, Data Assimilation is a suitable approach to solve the three drawbacks of Weickert's methods because: the evolution model describes the images' temporal dynamics; the covariance associated to the observation equation error overcomes the problem of missing data by weighting, at each pixel, the contribution of the observation equation in the computation of the solution; Data Assimilation is a frame-by-frame process allowing an affordable computational cost.

The Data Assimilation framework

The Data Assimilation framework aims to solve the system (1,2,3) with respect to a state vector X(x,t), depending on the spatial coordinate x and time t:

∂X ∂t (x,t) + Å(X)(x,t) = ε m (x,t) (1) X(x, 0) = X b (x) + ε b (x) (2) À(Y,X)(x,t) = ε O (x,t) (3)
Equation (1) describes the temporal evolution of X.

Å, called evolution model, is supposed differentiable. As Å approximately describes the evolution of the state vector, a model error ε m is introduced to quantify inaccuracies. Equation (2) establishes the initial condition of the state vector and the error is expressed by the background error ε b . Equation (3), called observation equation, describes the links, possibly complex and non-linear, between the observation vector, Y(x,t), and the state vector. The observation error ε O represents the imperfection of À and the measurement errors. ε m , ε b and ε O are assumed to be Gaussian and fully characterized by their covariance matrices Q, B and R.

Let X denote a Gaussian stochastic vector depending on a space-time coordinate (x,t). X = X(x,t) and X ′ = X(x ′ ,t ′ ) are the values on two locations. The covariance matrix Σ, computed for X and X ′ , measures their dependency and is defined by:

Σ(x,t, x ′ ,t ′ ) = ZZ (X -X) T (X ′ -X ′ )dP X,X ′ P X,X ′
is the joint distribution of (X, X ′ ) and denotes the expectation. P X,X ′ is a Gaussian distribution whose covariance matrix depends on x, t, x ′ and t ′ . The inverse of a covariance matrix is formally and implicitly defined [START_REF] Oliver | Calculation of the inverse of the covariance[END_REF]

by: ZZ Σ -1 (x,t, x ′′ ,t ′′ )Σ(x ′′ ,t ′′ , x ′ ,t ′ )dx ′′ dt ′′ = δ(x -x ′ )δ(t -t ′ ) (4)
with δ denoting the Dirac function.

We note A = Ω × [0, T], Ω being the spatial domain and [0, T] the temporal domain. In order to solve equations (1), ( 2) and (3), we define the functional E, to be minimized, as stated by equation ( 5

(ii) introducing an auxiliary variable λ, called adjoint variable in the literature:

λ(x,t) = Z A Q -1 (x,t, x ′ ,t) ∂X ∂t + Å(X) (x ′ ,t ′ )dx ′ dt
(7) and (iii) writing the state vector X as X b +δX. Vectors X b and δX are respectively called background vector and incremental vector in the Data Assimilation literature. This leads to the following system of equations: λ(x, T) = 0 (8)

- ∂λ ∂t + ∂Å ∂X X b * λ = - Z A ∂À ∂X X b * R -1 (9) À(X b , Y) + ∂À ∂X X b (δX) dx ′ dt ′ X b (x, 0) = X b (x) (10) ∂X b ∂t + Å(X b ) = 0 (11) δX(x, 0) = Z Ω Bλ(x ′ , 0)dx ′ (12) ∂δX ∂t + ∂Å ∂X X b (δX) = Z A Qλ(x ′ ,t ′ )dx ′ dt ′ (13) E(X) = Z A Z A ∂X ∂t + Å(X) T (x,t)Q -1 (x,t, x ′ ,t ′ ) ∂X ∂t + Å(X) (x ′ ,t ′ )dxdtdx ′ dt ′ + Z A Z A À(X,Y) T (x,t)R -1 (x,t, x ′ ,t ′ )À(X, Y)(x ′ ,t ′ )dxdtdx ′ dt ′ + Z Ω Z Ω X(x, 0) -X b (x) T B -1 (x, x ′ ) X(x ′ , 0) -X b (x ′ ) dxdx ′ (5)
The reader is referred to [START_REF] Valur Hólm | Lectures notes on assimilation algorithms[END_REF] for details on the determination of the differential of E and Euler-Lagrange equations. The decomposition of the state vector into background and incremental vectors requires linearizing the operators Å and À using a first-order Taylor development:

Å(X) ≃ Å(X b ) + ∂Å ∂X X b (δX). The operator ∂• ∂X * is called adjoint
operator and represents a compact notation for integration by parts of the differential of the operators Å and À. It is formally defined by:

∂Å ∂X (η), λ = η, ∂Å ∂X * (λ) (14)
Riesz's theorem ensures the existence and uniqueness of the adjoint operator. To solve the system ( 8) to ( 13), equation ( 11) is first solved using initial condition ( 10) and provides X b . Equation ( 9) is solved retrogradely to provide λ using initial condition (8).

The incremental vector is finally provided by solving ( 13) with ( 12) as initial condition. If Å and À

are not linear, an iterative algorithm is used to sequentially improve the solution of (8,9,12,13): at each iteration, the value δX is computed from the background vector X b and the value δX + X b becomes the background vector for the next iteration.

Image Assimilation

This section explains how to solve ill-posed Image Processing problems using the framework of Data Assimilation and constitutes the core of the research and the main contribution of the paper. Using Data Assimilation to solve Image Processing Problem is an emerging domain. Studies have been done about curve tracking (Papadakis and Mémin, 2007) and determination of optical flow (Papadakis et al., 2007a;[START_REF] Papadakis | Image assimilation for motion estimation of atmospheric layers with shallow-water model[END_REF]. In [START_REF] Huot | Assimilation of sst satellite images for estimation of ocean circulation velocity[END_REF], a method is proposed to estimate the ocean surface circulation from SST images by assimilating image data in an image model for producing pseudo-observations for the oceanographic model. In this paper, we consider the general case of ill-posed problems, solved using Tikhonov regularization methods, and we define a method to convert them, in a generic way, from Tikhonov regularization to Data Assimilation framework.

First, state and observation vectors have to be defined. Obviously, the observations will be images or processed images, and the content of the state vector will strongly depends on the studied application. For example, segmentation, denoising and restoration use a state vector which is an image. Tracking, image registration and motion estimation are represented by a vector field. Active contours are modelled by curves. Sometimes, the links between the state vector and the observation are highly indirect. For instance, as it is not possible to deduce the ocean circulation from surface temperatures by a shallow-water model, an additional system has to be considered with relevant evolution and observation equations [START_REF] Huot | Assimilation of sst satellite images for estimation of ocean circulation velocity[END_REF].

To solve an Image Processing problem with the 4D-Var algorithm, as described in Section 2, a suitable evolution model describing the temporal evolution of the state vector has to be stated. Its associated error must also be defined. Next, an observation equation and its error should characterize the properties of observations. Last, the initial condition and its background error have to be established. This is discussed in the next Subsections.

The evolution model and error

One first choice of evolution model is to consider the simple transport equation dX dt = 0, or:

∂X ∂t + ∂X ∂x ∂x ∂t = 0 (15)
The evolution model is then obtained by identification: Å(X)(x,t) = ∂X ∂x ∂x ∂t . Section 4 describes the use of equation ( 15) for determining optical flow.

Another possible choice is to express the transport of the state vector as a diffusion, a physical law ruling the transport of chemical species or temperature.

The general formulation is: One property of this diffusion is to smooth the image while preserving the contours. In this case, identification gives: Å(X) = -∇ T (D∇X).

∂X ∂t = ∇ T (D∇X) with ∇ = ∂ ∂x , ∂ ∂y
The matrix Q of the model error ε m has to be specified. It is first important to understand the action of Q -1 in the functional (5). Let us describe two possible choices for Q and analyze their impact in a functional

E(X) = RR F T (X)Q -1 F(X)dxdx ′ which has to be minimized. First, let Q being the Dirac covariance defined by Q(x, x ′ ) = δ x -x ′ .
We have E(X) = R Ω F(X) 2 dx: a Dirac covariance is corresponding to a zero-order regularization. A second choice is the exponential covariance de-

fined by Q(x, x ′ ) = exp - x -x ′ σ . We have E(X) = 1 2σ Z F(X) 2 + σ 2 ∂F(X) ∂x 2 dx.
The exponential covariance is associated to a first-order regularization.

In the general case, the inversion of a covariance matrix is often non-trivial and usually inaccessible. Restrictive choices have to be made such as the previous ones.

For further details, the reader is referred to [START_REF] Oliver | Calculation of the inverse of the covariance[END_REF][START_REF] Tarantola | Inverse Problem Theory and Methods for Model Parameter Estimation[END_REF][START_REF] Béréziat | Solving ill-posed image processing problems using data assimilation[END_REF].

The observation equation and error

As previously pointed out in Section 2, the observation equation describes the links between the state vector and the observations. As the regularization part of the Tikhonov's method does not contain observation values, the observation equation will be exactly corresponding to the Image Model existing in the Tikhonov's approach.

The observation error has to be stated and its covariance matrix R specified. R has then to be inverted (see equation ( 9)) and, as explained, this is only possible for simple cases. The Dirac covariance can be used to express a null interaction between two spacetime locations. This allows observation errors to be located and quantified in space and time. R is written as:

R(x,t, x ′ ,t ′ ) = r(x,t)δ x -x ′ δ t -t ′ (16)
with r a real invertible matrix and:

R -1 (x,t, x ′ ,t ′ ) = δ x -x ′ δ t -t ′ r -1 (x,t) (17)
The matrix r -1 characterizes the quality of the observation: a high value indicates that the observation is relevant. Assuming the availability of a function f measuring the confidence in observation values ( f ∈ [0, 1], f = 0 for no confidence), r -1 can be defined by:

r -1 (x,t) = r 0 (1 -f (x,t) + r 1 f (x,t) (18)
r -1 (x,t) will be close to a "minimal value" r 0 on regions where confidence is low and r 1 otherwise. r 0 and r 1 are constant and invertible matrices. For efficiency, f is modeled as:

f (x,t) = f noise (x,t) f sensor (x,t) f À (x,t) (19)
f noise characterizes noise acquisition: it is close to 0 for large noise. f sensor is equal to 0 when data are not acquired by the sensor. f À measures the confidence in the observation model; it is close to 0 if the observation equation is not valid.

Observation values with a low confidence will then not be taken into account for the computation of the state vector and for the solution of the Image Processing problem.

Determining optical flow

Let I be a sequence of images on a spatial domain Ω and W(x,t) be the velocity vector of a point x ∈ Ω at time t. The velocity vector is expressed by the Optical Flow Constraint equation [START_REF] Horn | Determining optical flow[END_REF]):

∇I T (x,t)W(x,t) + ∂I ∂t (x,t) = 0 ∀x ∈ Ω (20)
This is an ill-posed problem: the velocity vector has two components and the optical flow equation is not sufficient to compute both. A solution can be obtained using Tikhonov's method by constraining spatially [START_REF] Horn | Determining optical flow[END_REF] or spatially and temporally [START_REF] Weickert | Variational optic flow computation with a spatio-temporal smoothness constraint[END_REF] the solution. Let us explain, advantages of the data assimilation methods on this simple example.

Observation and evolution models

W(x,t) is considered as the state vector X(x,t) and the image gradients (∇I(x,t), I t (x,t)) constitutes the observation vector Y(x,t). By identification in ( 3) and ( 20), the observation model is:

À(W,I) = ∇I(x,t) T W(x,t) + I t (x,t) (21)
To determine a suitable observation error covariance handling missing data, we use equations (17), ( 18) and ( 19) to define the inverse covariance. f noise is assumed to be equal to 1 without further information. f sensor is set to 0 where data are not acquired and to 1 otherwise. f À is chosen from the fol- lowing remark: on a sequence of images with regions of uniform grey level values, the spatio-temporal gradient is null on these regions and equation ( 20) is degenerated. For avoiding considering these points, f À is:

f À (x,t) = 1 -exp(-∇ 3 I(x,t) 2 ) (22)
where ∇ 3 denotes the spatio-temporal gradient. The observation model being scalar, r 0 and r 1 are scalar and respectively set to ε and 1 -ε with ε ∼ 10 -6 .

The transport of the velocity is chosen as evolution equation. The evolution model is:

Å(W) = Å 1 (W) Å 2 (W) = UU x + VU y UV x + VV y T with (U,V ) T = W.
Let us denote I 2 the 2×2 identity matrix. Q is chosen as:

Q(x,t, x ′ ,t ′ ) = I 2 exp(- 1 σ ( x -x ′ + |t -t ′ |)) (23)
to ensure a first-order spatial regularization. We choose the Horn and Schunk's algorithm to compute the velocity field on the first image of the sequence. This initial condition is supposed to be without errors: B(x, x ′ ) = δ(xx ′ ).

Adjoint operators

In order to determine the adjoint operators of Å and À, the directional derivatives must first be established. Using the definition (6), we obtain:

∂Å 1 ∂W (η) = ∂Å 1 ∂U (η 1 ) + ∂Å 1 ∂V (η 2 ) = Uη 1 x + U x η 1 + V η 1 y + U y η 2 ∂Å 2 ∂W (η) = ∂Å 2 ∂U (η 1 ) + ∂Å 2 ∂V (η 2 ) = Uη 2 x + V y η 2 + V η 2 y + V x η 1 with η = η 1 η 2 T
. By definition ( 14) and considering border terms equal to zero, the adjoint operator of Å is :

∂Å 1 ∂W * (λ) = -Uλ 1 x -V y λ 1 -V λ 1 y + U y λ 2 ∂Å 2 ∂W * (λ) = -U x λ 2 -Uλ 2 x -V λ 2 y + V x λ 1 with λ = λ 1 λ 2 T .
The directional derivative of the observation operator is:

∂À ∂W (η)(x,t) = ∇I T (x,t)η(x,t)
and the adjoint operator is:

∂À ∂W * (λ)(x,t) = ∇I(x,t)λ(x,t)

Discretization

Using the choices explained in Subsection 4.1, the three PDEs (9,11,13) become:

∂W ∂t + W T ∇W = 0 (24) - ∂λ ∂t -∇λ T W -λ T ∇ ⊥ W = A ⋆ L (25) ∂δW ∂t + W T ∇δW + ∇W T δW = Q ⋆ λ (26) with: A = - 1 2σ 2 ∇IR -1 , L = I t + ∇I T (W + δW), ∇ ⊥ W = ∇ ⊥ U ∇ ⊥ V and ∇ ⊥ U = U y -U x T .
These three equations are further approximated using the finite difference technique.

As W(x,t) is a vector of Ê 2 , equation ( 24) has two components. The first one combines a linear advection in direction y and a non linear in direction x.

A direct approximation leads to an unstable numerical scheme. In (Papadakis et al., 2007a), the scheme is stabilized by introducing a diffusive term with the drawback to smooth the solution. We propose an alternative to stabilize the numerical scheme associated to the original equation, as explain in the following. The first component of equation ( 24) is first expressed as a two-equation system using a splitting method (Verwer and Sportisse, 1998):

∂U ∂t + UU x = 0 (27) ∂U ∂t + VU y = 0 (28)
Equation ( 27) is rewritten using the Lax-Friedrich method [START_REF] Sethian | Level Set Methods[END_REF] in the following form

∂U ∂t + ∂F(U) ∂x = 0 with F(u) = 1 2 u 2
. This new equation is approximated by the explicit scheme:

U k+1 i, j = 1 2 (U k i+1, j + U k i-1, j ) - △t 2 (F k i+1, j -F k i-1, j ) with U k i, j = U(x i , y i ,t k ), F k i, j = F(U(x i , y i ,t k )). The term 1 2 (U k i+1, j + U k i-1, j
) stabilizes the numerical scheme (by adding a diffusive effect) when △t has a low value (the Courant-Friedrich-Levy condition). The linear advection (28) is approximated using a shock explicit scheme [START_REF] Sethian | Level Set Methods[END_REF]:

U k+1 i, j = U k i, j -△tS x (U k ,V k ) i, j
with S the shock operator, defined in the x-direction by: S

x (U,V ) i, j = max(V i, j , 0)(U i, j -U i-1, j ) + min(V i, j , 0)(U i+1, j -U i, j ).
It can be seen that the second component of ( 24) contains a linear advection term in direction x and a non linear advection term in direction y. The same strategy is then used to perform the discretization process.

Equation ( 25) combines a linear advection, a term of reaction and a forcing term. It also has two components. The first is -

∂λ 1 ∂t -Uλ 1 x -V y λ 1 -V λ 1 y + U y λ 2 = A 1 ⋆ L with A 1 = -1 2σ 2 I x R -1 .
It is split into two parts. The first one contains the linear advection in direction x and the reaction term -∂λ 1 ∂t -Uλ 1 x -V y λ 1 = 0 and is approximated in the same way as ( 28) with an explicit shock scheme. However, the equation is retrograde as the initial condition is given at time T :

(λ 1 ) k-1 i, j = 1 + 1 2 (V k i, j+1 -V k i, j-1 ) (λ 1 ) k i, j + △tS x ((λ 1 ) k ,U k )
The second part contains the linear advection term in direction y and the forcing term: -∂λ 1 ∂t -V λ 1 y + U y λ 2 = A 1 ⋆L. Again, a shock explicit scheme is used:

(λ 1 ) k-1 i, j = (λ 1 ) k i, j - △t 2 U k i, j+1 -U k i, j-1 (λ 2 ) k i, j - △t(A 1 ⋆ L) k i, j
-△tS y ((λ 1 ) k ,V ) Having the same structure as the first component, the second component of ( 25) is approximated using the same method

The last equation, ( 26), is similar to equation ( 25) and the same discretization technique is applied.

Results

The "taxi" sequence and a synthetic sequence have been chosen for discussing results. For both cases, the initial condition is provided by the Horn and Schunk method applied on the two-first frames of the sequence. Image gradients, providing observations, are computed with a convolution method and a derivative Gaussian kernel. The parameter σ of the matrix Q is set to 1. The number of iterations is set to 5.

The taxi sequence displays several cars with a slow and quasi uniform motion. In a first experiment, we compute the optical flow on the sequence using the Data Assimilation method (D.A.). The Horn & Schunk method (H.&S.) is also applied and results, on the fifth frame, are displayed for comparison purposes on Figure 1. Results of both methods are similar in this case because observation values are available on the whole sequence obviously. Having chosen H.&S. or another method does not really matter. The issue is not comparing optical flow method but proving the efficiency of D.A. for missing data or complex dynamics. Then, a second experiment is designed to prove the capability of D.A. to manage missing data. To simulate a sensor failure, a large region around the white car is set to zero in a frame and the function f sensor returns zero on these pixels. Image gradients are then computed. Figure 2 shows the results for D.A. and H.&S. methods. This latter obviously fails to provide acceptable velocity vectors over the missing data region while D.A. successes due to the eviction of aberrant values in the computation. Even a whole frame of observation can be missing: we force image gradient values to zero on the fifth frame of the taxi sequence and f À ≡ 0 on this frame. Figure 3 displays results obtained with D.A. with and without image gradients observation on frame 5. Results are similar, due to the fact that the evolution model correctly describes the temporal dynamics of the sequence and compensate missing data. H.&S. method is no more able to provide any result in that case, conducting to a discontinuity in the output.

A third experiment is dedicated to the demonstration that the evolution model improves the estimation of the image dynamics. For that purpose, we build a synthetic sequence with two squares, one moving horizontally and one vertically. At the end of the sequence, they meet each other. Figure 4 shows the result of both methods at the beginning and the end of the sequence. H.&S.'s method fails to estimate a correct velocity direction when the squares meet. This is a common drawback of image processing method due to over regularization. D.A. however provides accurate results: the evolution model correctly describes the squares dynamics and decreases the effect of spatial regularization.

Conclusion

In this paper we described a framework to solve ill-posed Image Processing problems using Data Assimilation algorithms. This approach is an alternative to the Weickert's method, which constrains the solu- tion's variations in space and time.

General principles to choose a suitable evolution model, an observation equation, and their matrices of covariance were discussed in the paper. The observation equation corresponds to the Image Model of the Tikhonov regularization method. The evolution equation constrains the solution in time and the model error constrains it in space. As an illustration, we proved that an exponential covariance model error is equivalent to a first order Tikhonov regularization. The three drawbacks of the Weickert's method, pointed out in the Introduction, are overcome. First, the state vector's dynamics is described by a suitable evolution model. Second, a covariance R built from the confidence on the observation values weights the influence of these latter: observations having a low quality, such as missing data, are not taken into account in the computation of the solution. And last, the algorithm processes a sequence of images frame by frame, which allows an implementation with low memory requirements. The choice of the evolution model remains a key point as it greatly depends on the applicative context. We described it in the case of determining optical flow: image's grey level values are assimilated in an evolution model representing a transport equation of the velocity and a stable numerical scheme, approximating this equation, has been established.

The main perspective of the work is to define a generic form of the evolution model which will then be specified depending of the Image Processing application. Diffusion laws, modelling the transport of the state vector using the physical conservative principles, are good candidates for a number of reasons. First, their approximation in a discrete space is robust as stable numerical schemes exist. Second, the adjoint of the model operator is known explicitly. Third, several types of diffusion can be considered depending on their impact on the spatial and temporal regularization of the solution. An isotropic diffusion is equivalent to a uniform regularization but has the drawback to smooth the solution while an adaptive diffusion, driven by the state vector, allows discontinuities to be preserved. As the spatial regularization may also be performed by the covariance matrix Q, it could be interesting to compare the two paradigms. Moreover, the diffusion can also be driven by the observation values or by their confidence. For instance, a strong diffusion on the state vector may operate when the observation confidence is low which is an alternative way to deal with missing data. Again, a comparison between such a diffusion and the covariance matrix R built from the observation confidence must be led.

T.

  D is a tensor characterizing the direction and intensity of the diffusion. It is possible to drive this diffusion according to the image characteristics. A well-known example is the Perona & Malik diffusion: the tensor matrix is equal to c( ∇X )Id with c a Gaussian function and Id the identity matrix.

Figure 1 :

 1 Figure 1: Result on frame 5: D.A. (top), H&S (bottom).

Figure 2 :

 2 Figure 2: Missing data: D.A. (top), H.&S. (bottom).

Figure 3 :

 3 Figure 3: Result with gradients set to zero on frame 5 (top) and with gradient available (bottom).

Figure 4 :

 4 Figure 4: Result with D.A. (left) and H.&S. (right).

  ), displayed at the top of next page. If ε m , ε b and ε O are assumed to be independent, E represents the log-likelihood of

	X. The minimization is carried out using a calcu-
	lus of variations: the differential of E with respect
	to X, denoted Euler-Lagrange equation is solved. ∂E , is calculated and the associated ∂X ∂E is obtained ∂X by (i) computing the derivative of E in the direction
	η:
	∂E ∂X	(η) = lim γ→0	d dγ	(E(X + γη))


Ubaid Abbasi
email: abbasi@labri.fr

Toufik Ahmed

SWOR: an Architecture for P2P Scalable Video Streaming using Small World Overlay

Keywords: P2P Network, Small-World, Overlay Organization, Push-Pull content delivery, Quality of Service (QoS). I

In this paper we propose SWOR, an architecture to deliver Scalable Video Coding (SVC) contents over P2P network. SWOR architecture is based on twofold mechanisms: (1) organization of peers in Small-World (SW) overlay networks, (2) delivering SVC content using push-pull mechanism. SWOR is evaluated and compared with CoolStreaming/DONet (CS) on QoS metrics using NS-2 simulator. Performance evaluation demonstrates that SWOR outperforms CS and provides a significant improvement in the received video quality in terms of lowering the packet loss and transmission delay.

INTRODUCTION

Over the past few years, P2P networks have emerged as a promising approach for distribution of multimedia contents over a large scale network. Recently, we have been witnessing the emergence of novel P2P applications such as P2P audio and video streaming. P2P streaming focuses on the efficient delivery of audio and video content under strict timing requirement. For real-time video streaming applications over P2P networks, a single sender peer cannot provide sufficient bandwidth to ensure smooth video delivery. This makes maintaining a smooth quality of service (QoS) an important research challenge. Our target objective is to design a P2P streaming architecture which achieves the following goals: -Low Delay: Real-time streaming applications are sensitive to delay. Thus, an overlay with low delay is required for these applications. We intend to construct a P2P overlay which minimizes the delay from source to receiver peer. -Peer Dynamics: In P2P networks, peers join or leave the network without prior notification. The desired overlay should be able to handle peer dynamicity.

-Scalable and Adaptive: The architecture should be scalable for large number of users. It should be adaptive in the sense that it continuously improves the overlay based on the existing user characteristics. -Quality Adaptation: Assigning different parts of the video contents to specific peer in order to ensure stable QoS. We address the above mentioned issues in our proposed video streaming architecture (SWOR) by combining the characteristics of push and pull content delivery and by incorporating peers organization into clustered overlays forming a Small World (SW) [START_REF] Watts | Collective dynamics of Small World Network[END_REF]. In this architecture different peers receive different parts of video from multiple sender peers using the push and pull mechanism as shown in Figure 1.

In this figure a single peer receives video from multiple sender peers using the PUSH mechanism and the missing parts are retrieved using the PULL mechanism. In SWOR a receiver side scheduler is implemented to orchestrate the delivery of the different parts of video from different sender peers. All the video packets received from the sender peers are combined in different buffers that are monitored regularly to identify when we need to implement pull-based mechanism to receive the missing parts of the video before their actual playback deadline. The proposed architecture is based on Scalable Video Coding (SVC) that is considered more suitable for heterogeneous networks and receivers due to its real-time content adaptations capabilities. In SVC encoding scheme, each video stream is encoded into multiple video quality tiers.

The rest of the paper is organized in different sections. A brief related work and motivation is presented in section 2. The proposed architecture is described in details in sections 3 addressing major components involved in video delivery over P2P networks. Section 4 illustrates the performance evaluation and section 5 presents a brief conclusion while highlighting some of the future perspectives.

II. RELATED WORKS

There have been tremendous efforts in the design and experimentation of video streaming systems in the past two decades, yet no single system has delivered the expected scalability and QoS [START_REF] Yin | TrustStream: A Secure and Scalable Architecture for Large-Scale Internet Media Streaming[END_REF]. In this section, we highlight some of the known solutions that exploit the push and pull video distribution. CoolStreaming/DONet (CS) [START_REF] Zhang | DONet: "A data-driven overlay network for client live media streaming[END_REF] presents a framework for live media streaming based on data-driven overlay networks where peers periodically exchange data availability information with each other and retrieves the unavailable data. This framework has shown great improvement in the video distribution with high scalability, however it has two main drawbacks: (1) long initial startup delay due to random peer selection process, and (2) higher failure rate in joining a program during flash crowd. GridMedia [START_REF] Zhang | A Peer-to-Peer Network for Live Media Streaming -Using a Push-Pull Approach[END_REF] is a well known system offering P2P-based IPTV services. It organizes the peers in unstructured overlay networks and implements a push-pull based approach to fetch the media contents from the neighbor peers. The iGridMedia [START_REF] Zhang | iGridMedia: Providing Delay-Guaranteed Peer-to-Peer Streaming Service on Internet[END_REF] focuses on providing delay-guaranteed services to support real-time applications for large number of users. There have been several recent works, which focused on effective content searching using the properties of small world. There are several fundamental questions that are unclear such as how an overlay evolves under a random peer selection scheme, in what sense is the system scalable and what are the limitations and trade-offs. The existing works didn't provide a comprehensive study on these crucial issues. In this paper, we proposed an architecture to organize the peers in overlay networks forming the small world. We attempt to provide an efficient P2P media streaming model based on chunk driven philosophy and a unique peering scheme to counter the shortcomings mentioned earlier.

III. SCALABLE VIDEO STREAMING ARCHITECTURE

Two major components are designed and implemented for SWOR: (1) overlay organization mechanism to organize participating peers and (2) end-to-end content delivery, utilizing push-pull mechanisms. The proposed P2P architecture comprises of four entities: (a) Source peers, (b) Tracker, (c) Super peers, and (d) Ordinary peers. The source peers are the seeders of the original video contents. Tracker plays an important role in the proposed architecture that is responsible for assisting participating peers to communicate with each others to find their bootstrapping peers. When a new peer enters the network, it provides its unique identity and bandwidth range [IP, BW] to tracker. The detailed streaming architecture is presented in the following sub-sections.

A. Overlay Organization

We organize the participating peers into overlay network in the form of a small world (SW) forming different clusters. The small world represents a class of random graph in which every peer is accessible from other peer in small number of hops [START_REF] Watts | Collective dynamics of Small World Network[END_REF]. Furthermore, in small world every node has an intense local clustering and it shares an edge with some far located nodes in the network. This ensures that once a packet is in the cluster, it can be easily obtained by other members of same cluster. The motivation of using small world is to ensure robustness and guaranteed bandwidth that can be achieved by decreasing the number of hops between source and receiver peers from application (overlay) level perspective. To mathematically define the shortest hop property of small world, let G = (N, E) denote a connected graph representing a small world network. Let D (i, j) represents the length (in hops) of the shortest path between two nodes i and j .

V ∈ Definition 1: The average shortest hop count of a graph G, denoted as F (G), is equal to

( ) ∑ ∈ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = V j i j i D N G F , , 2 1 ) (
Eq.1 F (G) is the ratio of the sum of all shortest paths between any two nodes in G and all possible pair wise connections of the connected graph. The participating peers in SWOR are organized in overlay network of two layers. First layer is comprised of super peers, whereas ordinary peers are organized in second layer forming different clusters of small world. All ordinary peers in the same cluster are connected with each other forming a mesh network (intense clustering) and also connected to a super peer. Super peers are connected to each other in a flat hierarchical structure as shown in Figure 2. This structure is similar to small world network, where nodes are connected to each other (ordinary peer) and also have some links elsewhere in the network (super peer). Super peers are selected among ordinary peers on the basis of their serving strength towards the overlay network. The availability of multiple super peers addresses the issues of single point of failure. An alternative super peer is selected in the case when a super peer is no more available and/or unable to contribute for the video content delivery. Super peers are not static entities but they are dynamically selected and updated by the tracker. The overlay organization is carried by the super peers selection and peers joining and leaving mechanisms that are presented in following subsections:

1) Mechanism for Super Peer Selection

Source peers are responsible for distributing video contents to a fixed number of super peers that depends on the uplink bandwidth capacity of the source peer and the streaming rate of the video content as presented in Eq. 2. Initially, super peers are selected on the basis of their uplink bandwidth as provided to tracker while joining the network. Once streaming process starts, these super peers are periodically updated according to their contribution towards actual content delivery. Eq. 2 is also used to calculate the maximum number of ordinary peers in a cluster.

Source-Super Super-Ordinary Super-Super Ordinary-Ordinary Source Peer Super Peer

Ordinary Peer Cluster

Tracker is responsible for calculating the peer's contribution and updates the list of super peers. Tracker periodically performs the active measurement to calculate the serving strength (SS) of every peer. Serving strength is defined as the ratio of upload capacity to download capacity. This "SS" is further used to determine the contribution ratio (CR) of each peer (P i ) present in the small world. The distance of a peer (P i ) from the source can be its delay time that is the time taken by packets to arrive at the peer after the source has streamed out. Eq. 3 shows that the selected super peer will exhibit the best serving capacity and more localized to source peer. The contribution ratio based selection of super peers also introduced the incentive mechanism. Peers providing greater contribution to the overlay network are selected as super peers and thus, they are able to receive contents directly from source peers. This mechanism encourages, ordinary peers to contribute more to the overlay network.

The "CR" estimation is carried out by every peer and obtained information is sent to tracker periodically. This "CR" estimation may not be feasible due to its excessive overhead when there is large number of peers present in the network because every peer performs such calculations. Thus, we consider the ordinary peers having "CR" value greater than a certain threshold to limit the overhead. Tracker regularly updates the source peers about new super peers present in the network. Each source peer also provides its "CR" estimates to tracker. All the source peers are organized on the basis of their respective offered QoS. The best one is assigned high quality layer "Base Tier" and the enhancements tiers are assigned to the source peers in descending orders of "CR" estimates. Furthermore, we assume that no malicious peer is part of our small world because peer trust & reputation is out of scope of this paper. Hence, the provided information by each peer is accurate and is used to determine the best super peer selection.

2) Peer Joining & Leaving

When a new peer requests the tracker to join the network, tracker in return provides a list of available super peers with in same bandwidth range. Alternatively, if there is no super peer available that can accommodate any new peer, a list of ordinary peers having similar bandwidth range is provided. In both cases, the new peer selects a super/ordinary peer from the list and sends the join request. The receiving super peer can accept the request or reject it (depending on cluster size). In case of rejection or no reply, requesting peer timeout after a short time interval and continues its search. Once, the new peer joins the network its parent peer (super/ordinary) updates the tracker that keeps track of number of peers joining any super peer. This information is used for the future reference because each super peer can support only a limited number of ordinary peers.

Ideally when a peer leaves the network, it initiates a departure message to its parent peer that releases the outgoing bandwidth and notifies the tracker for the graceful departure of this particular peer. The other case is not trivial, when a peer suddenly leaves the network without any prior notification. In SWOR, sudden departure is encountered by periodical "keep-alive" message between tracker and super peers and super peers and ordinary peers. If there is no response received from a peer with in significant amount of time, it is assumed that the particular peer is no more available in the network. In return, some alternative super peer is selected to ensure the smooth content delivery.

B. End-to-End Content Delivery

The main objective of the proposed streaming architecture is to ensure the End-to-End content delivery with improved QoS. In this sub-section, we present our proposed push-pull mechanism over P2P networks.

1) Proposed Push-Pull Mechanism

The proposed architecture is facilitated by incorporation of push and pull modes of the content delivery. In the first steps, video contents are pushed from source peers to receiver peers after passing through super peers. In the second phase receiver peer implements pull method to acquire the missing packets of the video layers to construct the original video. The objective of the pushing component is to quickly distribute a data block to a certain number of peers, in order to fuel the subsequent pull based exchanges. We implement a buffer map for the video contents available at different peers. The buffer map provides information for the respective sequence number and the playback deadline for a video packet. This information is used to pull the missing packets from other members present in the cluster. These buffer maps are exchanged periodically among ordinary peers within cluster. Buffer maps are also exchanged among super peer for content sharing. If cluster members are unable to provide those packets, missing packets are requested from super peer. The flat hierarchical structure of super peer layer, allows fetching the contents from other clusters.

In our architecture, a receiver side scheduler plays an important role in video content sharing mechanism. Generally, there are a large number of source peers available in the network for popular video contents. However, in real-time video sharing it is not acceptable to communicate with a large number of source peers simultaneously. The scheduler determines which part of the requested video will be pulled from which peer in the cluster.

2) Quality Adaptation

Quality adaptation becomes important when a receiver peer select multiple sender peers to receive any video. We focused to use H.264/SVC [START_REF] Reichel | Joint scalable video model JSVM-8[END_REF] for encoding the original video contents that are used to produce highly compressed bit streams to generate a wide variety of bit rates. In such video encoding schema, original video stream is truncated into many different video layers. The "base layer" provides a significant proportion of the video quality whereas the "enhancement layers" are used to enhance the video quality in different dimensions. "Base layer" is considered as the most important layer because all the "enhancement layers" are only decodable with reference to lower layers and "Base layer". Initially the source peer distributes the "Base Layer" among all super peers. The enhancement layers are exchanged among the clusters using the flat hierarchical structure of the super peers. The packet ordering based quality adaptation not only ensures the base video quality for streaming session but also adds resilience for the lower video layers. The SVC based video coding offers the receiver to select the video either with high SNR quality, temporal resolutions for varying frame rates, or spatial resolution for the different resolutions.

IV. PERFORMANCE EVALUATION

This section describes the performance evaluations of SWOR for different QoS parameters using NS-2 simulations.

A. Simulation Setup

We performed intensive simulations to study the performance of SWOR in order to quantify the overall QoS parameters against CoolStreaming (CS) mechanism that implements push-pull based content delivery mechanism on an organized overlay structure.

Performance Metrics: The performance evaluation is carried for different QoS metrics that include: packet delivery ratio, jitter and video throughput. These parameters have significant role in determining the overall QoS for the real-time streaming applications.

Video Coding: We used H.264/SVC as encoding scheme. We truncate the original SVC video in different layers (base layer and enhancement layers) and distribute them among different source peers. We perform simulations to receive the video quality of CIF/CGSO with 15 fps. We generated four different quality tiers using MPEG-4 trace files [START_REF]MPEG-4 Trace Files[END_REF]. These quality tiers are generated by fractions of DCT (Discrete Cosine Transform) matrix. Base Tier (T0,S0,B0) offers 40% throughput of original video, enhancement tier 1 (T1,S0,B0) offers 30% throughput of original video, enhancement tier 2 (T0,S1,B0) offers 20% throughput of original file and enhancement tier 3 (T1,S1,B0) offers 10% throughput of the original video file.

Network Topology: We used the BRITE universal topology generator [START_REF] Medina | BRITE: Universal Topology Generation from a Users Perspective[END_REF] in the top-down hierarchical mode to map the physical network. The network topology consists of autonomous system (AS) and fixed number of routers. All AS are assumed to be in the Transit-Stub manner. Each topology consists of 8 autonomous systems each of which has 625 routers. This gives us about 20000 links in the topology. The delay on inter-transit domains and intra-transit domains are assumed to be 90 ms and 40 ms respectively, while delay on stub-transit is assumed to be 30 ms and intra-stub transit links are randomly chosen between 5ms and 30ms. The incoming bandwidth of super and ordinary peers varies between 512 kbps to 5 Mbps and is uniformly distributed throughout the network. We deployed multiple source peers providing different layers of original video. The uplink bandwidth of source peers are also randomly selected from above mentioned range. The maximum numbers of super peer connected to source peers are calculated using equation 2, described earlier. Peers join the network at random intervals. We also introduced sudden peer departure at random intervals. The presented results are the average results of the multiple runs of these simulations.

B. Results and Discussion

Figure 3 shows the average delivery ratio for SWOR in comparison to "CS". Delivery ratio is defined as the number of packets delivered before deadline over total number of packets encoded. Packets are proactively pushed to ordinary peers in the small world. These missing packets are requested from the other members of small world. If a packet is pulled successfully within its playback deadline it is consumed otherwise dropped. Figure 3 depicts that both "CS" and SWOR perform better with increasing swarm size, however "CS" has lower packet delivery ratio due to random content delivery as compared to our proposed small world architecture where content flow in an organized way. Due to existence of super peers and push-pull delivery in SWOR, packets have higher chances to be delivered to the connected clients in a smaller period of time. Packet delivery ratio for SWOR improves with increasing number of nodes. The reason is that increase number of clusters causes more contents to flow from one cluster to other. The obtained results indicate that the performance of SWOR remains almost same with increasing number of peers where more than 85% packets are received before playback deadline. Figure 4 shows the average overlay hops from source to destination peer. The organization of peers into small world decreases the average hop counts for SWOR.

Moreover, we evaluated the jitter rate and the results are presented in Figure 5. The unavailability of a packet at its playback time causes a jitter. SWOR has comparatively lower jitter rate as compared to CS. Sometimes an intermediate piece (data packet) which can not be fetched may increase the jitter rate. In our proposed architecture chances of missing packets are relatively low due to organized content flow. If a packet is acquired by a cluster member, it implies that there are at least one or more sources for the contents. This flow of contents saves times as compared to CS.

The throughput of different video tiers at the receiver peer is obtained and plotted for both SWOR and CS. We noticed that SWOR improved the received throughput as compared to CS. Figure 6(a), 6(b), 6(c) and 6(d) presents the comparison of original base tiers, enh. tier 1, enh. tier 2, and enh. tier 3 respectively. By analyzing the obtained results, we can notice that our proposed architecture provides smooth video delivery with higher quality and lower loss as compared to CS.

We combined the different received SVC layers at receiver peers to construct the original video that is further used to evaluate QoS at application level through the peak signal to noise ratio (PSNR) and structural similarity (SSIM) [11] index measurements. PSNR is commonly used for measuring picture quality degradation for the received image quality compared to the original image, whereas SSIM index measures the structural similarity between the original and the received image. The comparison results for PSNR and SSIM of the video generated in both architectures with reference to the original video are shown in Figure 7 and Figure 8 respectively. PSNR and SSIM measurements show a remarkable improvement in the received video quality in case of SWOR.

V. CONCLUSION & FUTURE PERSPECTIVES

In this paper, we introduced a new architecture for scalable P2P media streaming namely SWOR. It utilizes small world peering organization on overlay to deal with the existing problems in chunk based streaming systems. SWOR evaluation against CS resulted in improved packet delivery and better video quality. Results also showed a significant reduction of average number of overlay hops which leads to lowering the transmission delay. The proposed SWOR architecture incorporates SVC to ensure smooth video quality delivery with guaranteed adapted QoS. Finally, we believe that our results are promising and could provide research insight towards development of newer and efficient peering strategies in P2P media streaming systems.

For the future perspective, we aim to perform real testbed evaluation for the more personalized VoD and IPTV services delivery over P2P network.

Figure 1 :

 1 Figure 1: Architecture for Receiver Centric Video Streaming

Figure 2 :

 2 Figure 2: SWOR Overlay Organization

Figure 3 :Figure 6 Figure 6

 366 Figure 3: Packet Delivery Ratio

Figure 4 :

 4 Figure 4: Average Hop Counts

Figure 7

 7 Figure 7: PSNR Measurements

Figure 5 :

 5 Figure 5: Average Jitter Rate

Figure 8 :

 8 Figure 8: SSIM Estimates