Anne-Lise Chesnel
email: anne-lise.chesnel@cea.fr

Renaud Binet
email: renaud.binet@cea.fr

Lucien Wald
email: lucien.wald@ensmp.fr

An object oriented approach for quantitative assessment of building damage in urban areas using very high resolution images

Very High Resolution images are particularly well adapted to damage assessment in urban area because on one hand it allows an analysis focused on the buildings solely through an object-oriented analysis, and on the other hand it permits a quantitative evaluation of this damage assessment using a visually established ground truth. We proposed in this paper a method of damage assessment that uses these two benefits. First an original object oriented approach to register the images is presented. Then a simple and fast damage assessment method based on correlation is proposed and tested on several test-cases. Each buildings of a test-area is classified using Support Vector Machines. The performances of the method in each case is evaluated thanks to a manually constructed reference database that uses the European Macroseismic Scale. Our method shows good results, with up to 75% of buildings well classified among four different EMS damage grades.

I. INTRODUCTION

Natural disasters impact an increasing number of persons. A large part of the human and economical losses are due to the destruction of constructions, where people live and work. Damage assessment in an urban context is of interest because the concentration of buildings makes these implications more important. Very High Resolution (VHR) images, i.e. with spatial resolution better than 1m, offer a strong potential to achieve this evaluation because it allows identification of changes on buildings. However, this assessment is timeconsuming because of the number of objects to be analysed. An automated data processing could help [START_REF] Yésou | Personnal communication[END_REF]. Our present work focuses on the roof of buildings to assess damage. The perception of roofs is less affected by changes in viewing angle than other parts of a building. Moreover, roofs are visible by means of remote sensing and it is expected that they are less concerned with other changes. We estimate that, most of the time, damage corresponding to grades 4 and 5 on the European Macroseismic Scale (EMS) [START_REF] Grünthal | L'Echelle Macrosismique Europenne[END_REF] can be detected in VHR images by analysing roofs, because the roofs are often affected by heavy structural damage. However a collapse of a soft storey (e.g. ground floor) is for example undetectable on the roof itself. Grade 3 can be detected at times but in a less reliable way [START_REF] Yamazaki | Damage detection from high-resolution satellite images for the 2003 boumerdes, algeria earthquake[END_REF], [START_REF] Yamazaki | Visual damage interpretation of buildings in bam city using quickbird images following the 2003 bam, iran, earthquake[END_REF]. Damage inside the buildings and cracks on the walls cannot be detected, notably regarding for the spatial resolution of the actual sensors; buildings with a EMS damage grade lower that 3 will be considered as intact (grade 0) in our work. Thus, 4 distinct damage grades will be considered in this study. Methods based on VHR images face new challenges [START_REF] Adams | Application of high-resolution optical satellite imagery for post-earthquake damage assessment: The 2003 boumerdes (algeria) and bam (iran) earthquakes[END_REF]- [START_REF] Sakamoto | Automatic detection of damaged area of iran earthquake by high-resolution satellite imagery[END_REF]. The natural changes are more numerous as the resolution increases, leading to new difficulties: shadows changes, apparitions and disappearances of objects due to the human activity (e.g. cars, refugees tents). In addition the relative influence of errors in registration of building roofs increases [START_REF] Bitelli | Image change detection on urban area: the earthquake case[END_REF]. Pixel-based methods are more sensitive to these problems and show higher false-alarm rates [START_REF] Bitelli | Image change detection on urban area: the earthquake case[END_REF], [START_REF] Kouchi | Damage detection based on object-based segmentation and classification from high-resolution satellite images for the 2003 boumerdes, algeria earthquake[END_REF]. They may not allow to distinguish damage from other changes. Recently, more attention is given to region-based methods with VHR images, reducing a part of the problems quoted above. We have chosen such an approach. However, the task of building extraction is not dealt within this work, since it represents a problem in itself [START_REF] Girard | Building detection from high-resolution color images[END_REF]- [START_REF] Meinel | A comparison of segmentation programs for high-resolution remote sensing data[END_REF]. Instead, we suppose that we are in the ideal situation where we have an a priori knowledge of the buildings footprint in the reference image. Our approach is divided in four steps. Firstly, the footprints of the roofs of buildings are segmented in the reference image. The roofs in the two images, the reference and crisis images, are then registered by a correlation method. In a third step, pixel-based tools for change detection are used to assess the damage on the roofs alone. The buildings are classified in the fourth step. A reference database has been built independently and is used as a reference for registration of roofs (second step), training for classification (fourth step) and assessment of the performance of the approach. This approach is applied to several cases: the Bam (Iran) earthquake of December 26, 2003; the Boumerdes (Algeria) earthquake of May 21, 2003; the Muzzafarabad (Pakistan) earthquake of October 8, 2005; and the Ryongchon (North Korea) train explosion of April 22, 2004. Two VHR sensors are used: QuickBird and Ikonos, with respectively a spatial resolution of 60cm and 1m. The evaluation of the performance of our method is made by comparison to the reference database. The effectiveness of the method is quantitatively evaluated as a function of the test case, of the spatial resolution of these images, and of the damage degree. Thus the application domain and the limits of the registration method are investigated.

II. IMAGES AND REFERENCE DATABASES

The images used in this study are panchromatic VHR images from QuickBird [16] and Ikonos [17] sensors (table I). These sensors respectively enable spatial resolution equal to 60cm and 1m at nadir. The images of each cases of study have different acquisition parameters (table I). The elevation and azimuth of the sun are also different due to the season change. Different combinations of images from the two sensors acquired before and after each events are used (table II). In the case of a multi-sensor images (e.g. QuickBird/Ikonos), the image with higher resolution is downsampled to the resolution of the other. For that purpose, an apodised sinc, with a window size equal to 21 pixels is used. The image before the event is the reference image, and the image after is the crisis image. The registration method is applied to several cases of urban disaster: the Bam (Iran) earthquake of December 26, 2003; the Boumerdes (Algeria) earthquake of May 21, 2003; the Muzzafarabad (Pakistan) earthquake of October 8, 2005; and the Ryongchon (North Korea) train explosion of April 22, 2004. In each test case, the buildings and the damage on them differ by some aspects. We investigate the robustness of our registration method with these varied cases. We have constructed one database per case to serve as a reference to assess the performances of our approach. It contains the footprints of the roofs of the buildings prior to the disaster. The footprints were manually extracted from the reference images. For example, in the case of the Bam earthquake, 2168 buildings are contained in this database, and they are located in the eastern part of the town. An extract of the Bam QuickBird reference image on which the footprints of the buildings have been overlaid is presented in figure 1. The other databases contain respectively 610 buildings in Boumerdes, 725 buildings in Ryongchon and 937 buildings in Muzzafarabad.

III. REGISTRATION

Standard methods for image registration do not allow registration of buildings of unknown heights when the off-nadir viewing angles differ from one image to another [START_REF] Bitelli | Image change detection on urban area: the earthquake case[END_REF], [START_REF] Thönnessen | Change detection in satellite images[END_REF]. Thus the footprint of the buildings roofs have to be registered prior to the analysis of changes, by computing individually the offset for each building. First, the two images are roughly registered; the ground registration error is estimated to be at most three pixels. Then the registration of the buildings roof is conducted with a subpixel accuracy. The buildings in the crisis image do not need to be segmented to extract the buildings; the registration process uses solely the footprint of the buildings of the reference database.

A. Image registration

The QuickBird images are orthorectified by the mean of the SRTM Digital Terrain Model (3 arcsec) to compensate for the ground elevation, the Rational Polynomial Coefficients (RPC) provided with the QuickBird data to model the sensor and a set of four ground control points in order to compensate the pointing error. The Ikonos images were already orthorectified. The crisis image is registered to the reference image using a cubic convolution interpolation method. For both type of images, the remaining registration error of the ground is due to the imprecision of the DTM and/or RPC modelling error. It can be considered as locally constant. It is constant in the image for an area with smooth relief; steep relief can make it more varying if not accurately corrected. In addition, the images having different viewing angles, roofs of buildings are not registered precisely; considering each pair of images, as the offset is proportional to the height of the building, it can differ for each roof. We have devised a method for assessing and correcting those offsets for each roof.

B. Roof registration

Our registration method estimates the maximum a posteriori of a similarity measurement. This measurement is the correlation (equation 1) computed on the group of pixels held in the roof footprints on both images. In the reference image, the footprint is specified by the manual segmentation. This footprint is translated in the crisis image by a quantity o estimated by the algorithm of roofs registration (equation 1). Let I 1 and I 2 be respectively the reference and the crisis images. Let c and l be respectively the column and the line coordinates, and k c and k l be respectively the shift value of the footprint of the roof in the crisis image along the columns and along the lines.

R i (c, l) = I i (c, l).H(c, l) - c,l I i (c, l).H(c, l) c,l H(c, l) , for i = {1, 2}, o = arg max kc,k l            c,l R 1 (c, l).R 2 (c + k c , l + k l ) c,l R 2 1 (c, l). c,l R 2 2 (c + k c , l + k l )            , for {c, l|H(c, l) = 1} (1)
where The function H(c, l) may be defined by the means of the segmentation of buildings footprint as discussed earlier.

H(c, l) = 1 if I 1 (c, l) ∈ roof footprint 0 otherwise o r
In our case, its values are assigned by the means of the reference database of buildings. The position (k * c , k * l ) of the maximum correlation yields the offset o. Equation ( 1) is applied to a limited search area. Defining this search area is the first challenge. If it is too small, one may not find the true maximum correlation and underestimate the offset. Conversely if it is too large, it is likely that the correlation peak does not correspond to the right result, especially in urban area where edge features are numerous and have often the same orientation. The second challenge is that when considering a damaged building, the correlation has no sense because there is no similarity at all between the roof in the reference image and the corresponding place in the crisis image. One should find a way to retrieve the correct location of the building remains. In the following we propose a method to find out an adapted search area. The offset o can be written as the sum of a ground offset o g due to the ground misregistration and a height offset o h due to the unknown height of the building:

o = o g + o h . Let o ⊥ be such that o g ∧ o ⊥ = 0.
The ground offset o g has no expectable orientation. On the contrary, the height offset o h is solely dependent on the viewing angles for both images (direction) and on the height of the building (amplitude). Except for pancaked buildings (soft story collapse), this offset is observed in an epipolar-like direction θ given by the equation (2). Figure 2 illustrates how θ is defined as a function of the viewing angles of the two images.

θ = arctan tan e 1 sin a 2 -tan e 2 sin a 2 tan e 1 cos a 2 -tan e 2 cos a 2 (2)
where a 1 , e 1 , a 2 , e 2 are respectively the satellite azimuth and elevation angle above horizon of the reference and crisis images (table I).

The height offset is due to the height of the constructions, so its amplitude differs for each building. On the contrary, we will suppose that the ground offset does not vary much. These differences are exploited to define the search area. Note that a 2D search is necessary due to the unpredictability of the ground offset, reducing the relevance of an epipolar image projection.

The search area is the same for all buildings. It is defined by the angle θ and by two search intervals defined along the axes o h and o ⊥ (figure 3). These search intervals are different along these directions. The roofs registration is divided in two steps.

A first approximation of the search area (dotted rectangle in figure 3) is defined along the o ⊥ axis by setting the search interval to [-s, +s], where s is twice the expected registering error of the ground between the images (3 pixels in our case).

Along the o h axis, the interval is set to [-s, s + s h ], where s h is twice the roof offset associated to the highest building and is defined by:

s h = 2h max sin a 2 tan e 2 - sin a 1 tan e 1 2 - cos a 2 tan e 2 - cos a 1 tan e 1 2
(3) where h max is the estimated maximum height of all the buildings. An exhaustive computation over this whole first search area is conducted, and the correlation peak corresponds to the approximate roof offset of each building. A 2D histogram of the collection of offsets is built (figure 4). One may note, in the search area, points having very small values of occurrence which are outliers. It corresponds as a majority to damage buildings that do not exhibit correct offsets. In this case, they would corrupt the distribution of the offsets and should be filtered out. For that purpose, the histogram is fitted by a 2D Gaussian function, whose mean (µ h , µ ⊥ ) and standard deviation (σ h , σ ⊥ ) are estimated with a least square fit. The choice of a Gaussian function instead of another one is not very important since it serves only to compute a first approximation. The approximation of the roofs offset is then refined for each building. The offsets are computed on the new search area (plain rectangle in figure 3) centred on (µ h , µ ⊥ ) with intervals on each axis set to [-8σ h , +8σ h ] and [-4σ ⊥ , +4σ ⊥ ] (figures 3 and 4). Note that if there is no outlier, this computation would have no other action than centred the search area and refined the related tolerances in a better way.

Practically, the correlation values are computed by steps of 0.25 pixel in each orientation in the search region previously defined. A subpixel step is used for correlation, as the offset is not supposed to be a multiple of a pixel size. Thus it is necessary to compute a resampled version of the crisis image. For the first step of the registration, the bilinear interpolation is chosen as a good trade-off between efficiency and computational cost. We estimate that, except for the outlier values, the correlation peak is estimated with a 0.5 pixel accuracy. This is visible with the histogram in figure 4: offsets with an integer value seem more likely. This is an artefact generated by the bilinear interpolation. The same computation can be done using a apodised sinc instead of the bilinear polynomial (figure 5): the artefact disappears but the computation for this step is very expensive (at least ten times more in our case). For the second step of the registration, an apodised sinc is used to resample the crisis image. The correlation peak value is estimated using a truncated-Newton optimization method to speed up the computation. The location of the peak, i.e. the roof offset, is thus estimated with a subpixel accuracy.

To improve the effectiveness of the correlation, the images are filtered by a Canny filter [START_REF] Canny | A computational approach to edge detection[END_REF]. It emphasizes the important features (edges, elements on the roof) that characterize the buildings. This filtering especially allows a reduction of the radiometric changes in low frequencies and high frequencies (noise, aliasing) that are present in the images.

C. Registration results

For the buildings with EMS damage grade 0 and 3, the standard deviation (SD) of the registration error is lower than 1 pixel for the test-case of Bam. It is between 1 and 2 pixel for the intact building for test-case of Ryongchon. For the Boumerdes and Muzzafarabad test-cases, the SD of the registration error for grade 0 and 3 is higher; if most of the buildings are precisely registered, some exhibits large errors due to radiometric natural changes. The more the damage, the less precise is the registration. The buildings with a collapse of a soft storey in Boumerdes are shown to be difficult to register, because although their roofs are intact, their position and orientation are often unexpected. Finally, another difficulty is pointed out: the imprecision of the SRTM is more sensitive in an area with steep relief (Muzzafarabad), leading to geometric deformations uncorrected by the orthorectification (building deformation, location errors).

IV. DAMAGE ASSESSMENT

The roof of the buildings being registered, one can compare the roofs and thus assess the damage on them between the two dates of acquisition. The correlation coefficient that corresponds to the best offset for each building yields information about the amount of similarity. We use this coefficient to classify the buildings among the different detectable damage classes. We have chosen a supervised classification for the following reasons: we consider that we have no a priori knowledge about the number and the different degree of damage classes that are present for a given disaster; moreover, for a given damage degree, the damage feature (here the correlation coefficient) can differ from one test-case to another, according to the type of buildings, to the atmospheric conditions, to the acquisition parameters, etc. However, the required training set is constrained to be small, with five examples per damage grade. In our case, these examples are randomly selected among the reference database. Finally, because the classification result depends more or less on the chosen training patterns, 100 classification tests are conducted for each case with different training sets, to ensure the robustness of the results.

A. Damage assessment as a function of the test-case

First we investigate the quality of the damage assessment as a regard to the test-case. The rate of buildings that are well classified is computed respectively for every test-case over the whole corresponding reference database, without taking into account the damage degree. Moreover we only consider in this section the pair of QuickBird images so as to compare fairly the results. For the test-case of Bam, the mean classification performance is 75%, the best performance being 78%. For the Boumerdes test-case, the performance is 70%, with a maximum equal to 72%. For the Muzzafarabad test-case, only 36% of the buildings are well classified, with a maximum performance equal to 39%. Finally, the buildings of the database of the Ryongchon test-case are classified with an mean accuracy of 78% and a maximum performance of 79%. First, we note the stability of the results as a function of the selected training examples: the mean and the maximum classification performance differ only of 2 to 3 points. It shows that even if the size of the training set is small, the classification method is robust. We note that the performance are notably lower for the testcase of Muzzafarabad. The damage assessment in this case is unsatisfying. We investigated more precisely this case and noticed important geometric deformations on the buildings between the reference and the crisis images. This is due mainly due to the steep relief that is ineffectively corrected during the orthorectification with the SRTM; the resolution of the DTM is too low. By consequence, the pixel-based similarity measurement used for the registration and for the damage assessment is less efficient. For the other test-cases, the results are quite close and seem satisfactory.

B. Damage assessment as a function of the sensor

The different images pair that have been investigated per test-case are reported in table II. For the test-case of Bam, the classification performances are equal to 75% for the QB/QB case, to 72% for the QB 1m /QB 1m case and to 67% for the QB/IK case. In this case, the better the resolution, the better the classification performance. We also see that, for the same spatial resolution, the mono-sensor case gives better results than the multisensor case. For the test-case of Boumerdes, the original QuickBird pair also given better results (70% of buildings correctly classified) than the downsampled ones (64% of buildings correctly classified). For the test-case of Ryongchon, the difference of performance between the QB/QB pair and the downsampled one is lower (from 78% to 76%). It could be due to the type of damage that is considered for this case of explosion: damage are either highly visible, or undetectable on the images, whatever the resolution considered. Finally, for the test-case of Muzzafarabad, the classification performance is equal to 36% for the QB/QB case; regarding this low performance, the downsampled case has not been tested. The IK/IK case gives a result of 32% of buildings well classified. Thus, even if the performances are low, the trend of decrease of performance along with the resolution is confirmed.

C. Damage assessment as a function of the visually estimated damage degree

In this section, the performance of classification is evaluated per damage grade, and the global performance is computed giving the same weight to each class (instead of giving the same weight to each building independently of its damage grade). Because, for each test-case, every damage class has not the same number of buildings than the others, we first note that the global performances differ when evaluating the classification accuracy per damage class. Moreover, classes corresponding to EMS grade 0 and 5 contain more buildings than the intermediate grades. The classification results per class are reported in table III. For Bam and boumerdes test-cases, the buildings with damage grade 0 and 5 are the most accurately classified. The damage grade 4 is detected with a notably best accuracy than the damage grade 3. Indeed, in the case of Bam, grade 3 corresponds to slight damage on the roof, that do not impact strongly the correlation comparing to the natural changes (shadows for example) that appends on most of the intact roofs. In the case of Boumerdes, grade 3 corresponds to buildings with no visible damage on their roof, and with rubbles/dust at their base; hence they are not detected with our method. For the test-case of Ryongchon, the intact and totally destroyed buildings are detected with the best accuracy. The intermediate class presents more confusion with the other damage classes. For the Bam, Boumerdes and Ryongchon test-cases, the classification results per class are stable; conclusion can be drawn without a loss of generality. Conversely, for the Muzzafarabad test-case, there is an important confusion between the different damage classes. The correlation coefficients are low, whatever the grade considered. The trend is to classify the buildings alternatively between only two classes, the two others being almost ignored. One example of result is shown in table III.

V. CONCLUSION

We propose in this paper a method to assess damage on buildings in an urban area using a pair of very high resolution images. Its performances are evaluated on several cases, for four different disasters and using images from two sensors. We use as a reference the European macroseismic Scale, and decide to distinguish 4 different damage grades for the buildings in the images. First, a solution to the problem of registration of building roofs is presented. A method that takes into account the conditions of the image acquisition is proposed. It is based on correlation. The less the damage on the buildings, the better the registration. The roof registration is particularly efficient for the buildings of EMS damage grade 0 and 3. This preliminary important step of registration being accomplished, we show how the correlation coefficient can be an appropriate feature to efficiently and rapidely classify the buildings in several damage grades using SVM along with a small training set. For three of the four test-cases, the classification performances are between 70% and 80%. The last test-case show the limitations of the method for damage assessment. The reasons of the difficulties encountered are investigated, and we propose some possible solutions. One could improve the damage assessment if some of the natural changes between the images could be managed. This is notably the case of the shadows created by some elements on the top of the roof, of pixels with specular reflexion and of the geometric deformations uncorrected by the image preprocessing.

These results shows that having a database of the footprint of the buildings before a disaster can lead to an accurate and reliable damage assessment with very high resolution images. GIS of cities tend to be built more and more in order to achieve for example better urban development and environmental management. Such damage assessment using available GIS has yet to be tested.

Fig. 1 .

 1 Fig. 1. QuickBird reference image (before the earthquake) of Bam. The buildings footprint of the corresponding reference database is overlaid. (a) : entire scene; (b) : image extract. Colour legend of the database: Green = damage grade 0; Yellow = damage grade 3; Orange = damage grade 4; Red = damage grade 5.Images: copyright DigitalGlobe 2003, 2004.

Fig. 2 .

 2 Fig.2. Illustration of the direction θ, along which the roof projection of the buildings moves between the images, as a function of the angles of acquisition of each image.

Fig. 3 .

 3 Fig. 3. Definition of the referential used for the registration computation. The ellipse depicts the 2D histogram of the offsets.

Fig. 4 .

 4 Fig. 4. A 2D histogram of the roof offsets computed in the first search area. The cross stands for the mean of the Gaussian fitted distribution; the dotted box will be the second search area obtained from the estimated distribution. The bilinear interpolation method used generates artefacts corresponding to integer values of offset in the og and o ⊥ directions.

Fig. 5 .

 5 Fig. 5. A 2D histogram of the roof offsets computed in the first search area. The cross stands for the mean of the Gaussian fitted distribution; the dotted box will be the second search area obtained from the estimated distribution. An apodised sinc interpolation method used.

TABLE I IMAGES

 I USED. QB=QUICKBIRD AND IK=IKONOS

	Area	Sensor	Date	Sat. elevation Sat. azimuth
				angle above	(clockwise
				the horizon	due North)
	Bam	QB	2003-03-30	78.8 o	191.9 o
	(Iran)	QB	2004-01-03	63.3 o	233.8 o
		IK	2004-01-04	62.9 o	132.1 o
	Boumerdes	QB	2002-04-22	78.6 o	352.8 o
	(Algeria)	QB	2003-05-23	64.2 o	276.5 o
		QB	2004-08-13	79.4 o	151.4 o
	Muzzafarabad	QB	2005-10-22	64.8 o	118.6 o
	(Pakistan)	IK	2002-09-15	71.7 o	227.8 o
		IK	2005-10-09	79.7 o	306.7 o
	Ryongchon	QB	2003-05-13	59.9 o	274.1 o
	(North Korea)	QB	2004-04-27	72.6 o	282.8 o

TABLE III CLASSIFICATION

 III PERFORMANCES PER DAMAGE GRADE USING THE PAIRS OF ORIGINAL QUICKBIRD IMAGES.

				Events
		Bam	Boumerdes Muzzafarabad	Ryongchon
	EMS grade 0 / intact	81%	79%	22%	82%
	EMS grade 3	20%	0%	52%
	EMS grade 4	68%	31%	4%
	EMS grade 5	91%	91%	87%
	Partially damaged				45%
	Totally destroyed				86%

1-4241-4244-0712-5/07/$20.00 c 2007 IEEE

ACKNOWLEDGMENT

The authors would like to thank...


L Wald

In-flight interband calibration of the AVHRR data by a cloud-viewing technique

A significant degradation in the responsivity of the AVHRR radiometers aboard the NOAA satellite series, affects the index vegetation (NDVI), which is an important source of information for monitoring vegetation conditions on regional and global scales. Many studies have been carried out which use the viewing Earth calibration approach in order to provide accurate calibration correction coefficients for the computation of the vegetation index using the visible and near-infrared spectral channels 1 and 2 of AVHRR. This study deals with the interband calibration of AVHRR visible and near-infrared data by means of an cloudviewing technique. This technique is simple to implement and can be used in real-time. It is also well-suited to the processing of large time-series of data. Results are presented for various NOAA satellites and are in full agreement with the calibration degradation model proposed by NOAA and various authors.

INTRODUCTION

The visible and near-infrared data (channels 1 and 2) of the Advanced Very High Resolution Radiometer (AVHRR) aboard the NOAA satellites are widely used in quantitative applications in vegetation studies by means of the Normalised Difference Vegetation Index (NDVI). The AVHRR instrument is not equipped with on-board calibration of channels 1 and 2. Degradation of the responsivity of these channels cannot be accounted for and comparison of data between years cannot be done. Successive AVHRR radiometers may even induce discontinuities in timeseries of NDVI. It prevents the generation of timeseries of accurate vegetation index, which are of the highest interest in the monitoring of the vegetation for global change, for crop survey or disease warning. Many studies have been carried out which use the Earth calibration approach in order to provide accurate calibration correction coefficients for the computation of the vegetation index. Variations in time of these coefficients may be important and should be taken into account. This communication presents an accurate, easy-to-use method to calibrate vegetation index derived from AVHRR data.

The Earth viewing calibration approach has been recently developed as a backup solution to the possible failure or unreliability of on-board calibration devices. It is based on the knowledge of physical characteristics of some Earth phenomena as well as upon the processing of the digital imagery flowing down from the sensor itself (see e.g. Abel 1990).

If  i is the calibration coefficient, C 0i the deep space count for band i,  s the solar zenith angle and taking into account the solar extraterrestrial flux F 0i , the relation between the reflectance  i derived from the spectral radiance detected by the radiometer and the integer count value on a computer tape C i is:

 i =  i (C i -C 0i ) / cos ( s ) with  i =  i  / F 0i (1)
The change in the sensor calibration r i between the true calibration  i and the pre-flight calibration coefficient  i * given on the computer tapes supplied by the NOAA is:

r i =  i * /  i (2)
All published studies (see e.g., Abel et al. 1990;[START_REF] Rao | Degradation of the visible and near-infrared channels of the Advanced Very High Resolution Radiometer on the NOAA-9 spacecraft: assessment and recommendations for corrections[END_REF][START_REF] Rao | Post-launch calibration of the visible and near-infrared channels of the Advanced Very High Resolution Radiometer on NOAA-7, -9, and -11 spacecraft[END_REF] found that the deep space values did not change substantially with time and were close to the pre-flight values. Hence, the procedure for calibration correction only consists from the derivation of the coefficients  i or r i . If  i * denotes the reflectance computed using the pre-flight calibration coefficient:

Hence from here on the calibration procedure for vegetation index consists in the accurate evaluation of the interband calibration by the means of the coefficient r 12 .

CLOUDS AS TARGETS FOR CALIBRATION

Some Earth's bodies present spectrally flat reflectances in the solar spectrum and the interband calibration is done by adjusting the coefficient r 12 so that the ratio of the spectral reflectances observed over these targets is equal to 1, if other atmospheric effects can be corrected or do not change this ratio. Clouds are such bodies. The very reflectance of the cloud is spectrally constant in most of the solar spectrum. Such a property is often used in procedures for cloud detection (see e.g., [START_REF] Saunders | An automated scheme for the removal of cloud contamination from AVHRR radiances over western Europe[END_REF][START_REF] Saunders | An improved method for detecting clear sky and cloudy radiances from AVHRR data[END_REF][START_REF] Wald | Automatic cloud screening in NOAA-AVHRR day-time imagery[END_REF]. Moreover clouds are very frequent in satellite imagery and this renders an operational procedure possible.

At satellite level, the reflectance observed over a cloud is a function of the very reflectance of the cloud, the albedo of the underlying surface, the water vapour distribution of the environment in which the cloud is located, the clear-sky layer above the clouds in which molecular and particle scattering take place, and gaseous transmittance. If the reflectance of the ground is large enough with respect to the optical depth of the cloud, some of the spectral variations in reflectance observed at satellite level are due to the spectral variations of the ground reflectance. Also multiple reflections between the ground surface and a highly reflective cloud base may become important. Such cases must be avoided and only clouds over ocean areas should be examined, out of the sunglint area.

Asmami, Wald (1993) examined the potentials of the cloud-viewing technique for the interband calibration and for various instruments, with special emphasis on the calculation of the NDVI vegetation index. They analysed NOAA-9 data over several parts of the world ocean, and used the radiative transfer model of [START_REF] Paris | A cloudy-sky radiative transfer model suitable for calibration of satellite sensors[END_REF] for a cloudy atmosphere in order to guide the processing and help in understanding the results. It appears that with respect to the scope of the study, the most important among the various parameters are the optical depth and the altitudes of the base and top of the cloud. Under reasonable conditions, some clouds may serve as targets for interband calibration. In particular the clouds must have a large optical depth but too much reflective clouds (i.e. very large optical depths) are not suitable. Reasonable changes in the granulometry of the cloud have little effects on the interband calibration for such clouds. Attention has been paid to the influence of the viewing angle on the ratio; it was found that the influence of the viewing angle on the interband calibration is negligible. From their findings, it appears feasible to use reflectances observed at satellite level over clouds for interband calibration, at least for the spectral bands not located in oxygen or water vapour absorption windows.

They proposed a method for the particular case of the AVHRR sensor, which they calibrated against the calibration degradation model proposed by [START_REF] Kaufman | Calibration of the AVHRR visible and near-IR bands by atmospheric scattering, ocean glint and desert reflection[END_REF], for eight NOAA-9 images, most of them acquired in summer 1988. While they should have found a value of 1.0, they found 1.01 with a standard-deviation of 0.01. The accuracy was similar to that claimed by [START_REF] Kaufman | Calibration of the AVHRR visible and near-IR bands by atmospheric scattering, ocean glint and desert reflection[END_REF].

The effects of viewing angle upon the correction factor r 12 have been examined by processing portions of images and have been found negligible if any. This confirms the theoretical findings (see also [START_REF] King | Determination of the spectral absorption of solar radiation by marine stratocumulus clouds from airborne measurements within clouds[END_REF]).

AN ALGORITHM FOR INTERBAND CALIBRATION

The Asmami, Wald method is the core of the proposed algorithm. It has been slightly improved and constraints have been added to increase reliability and robustness for a better operational implementation. The algorithm selects within an image the clouds which are appropriate for interband calibration. It is based upon the results of the numerical simulations and is made up of a set of criteria dealing with the density of probability of the visible reflectances of clouds. They express the fact that only clouds of medium reflectivity are an ideal target. They have been completely defined by the analysis of the NOAA-9 images. Optical thickness of appropriate clouds is comprised between 10 and 70. The suitable clouds have been identified as being a very widespread layer of low clouds, stratus and stratocumulus. If temperature observed in AVHRR4 is related to standard vertical profiles of temperature and pressure, the altitudes of the top of these clouds are mostly comprised between 2 and 3 km, and range from 1 to 6 km. Rain may be observed under these clouds.

However within a cloud field optical properties may fluctuate strongly and it is necessary to investigate the density of probability of the reflectance of clouds. For example, very bright clouds may have portions of lower reflectance, which will be considered as clouds of medium reflectivity but which are still inappropriate for the calibration procedure. Thus the probability of having very bright clouds must not be too large. Similar reasoning holds with low reflectance clouds which exhibit bright portions. If too much such clouds are present, the sample of pixels representing the medium reflectivity clouds is corrupted. Indeed the density of probability should be fairly large for medium reflectances. In the following, the density of probability is approximated by the means of the histogram of the cloudy pixels.

Therefore the algorithm puts conditions on the population of ocean pixels of medium reflectances. It applies to cloudy pixels, and is the following: -define five classes of reflectances: 0.4-0.5; 0.5-0.6; 0.6-0.7; 0.7-0.8; 0.8-0.9. Typical cloud optical thickness is respectively 10, 15, 20, 30, 70 and 200 (see e.g. [START_REF] King | Determination of the scaled optical thickness of clouds from reflected solar radiation measurements[END_REF]). Define a sixth class which lower limit is 0.4 and the higher limit, the highest reflectance a cloud can take. This class contains all the cloudy pixels, having a reflectance greater than 0.4. This constitutes the population of cloudy pixels; -compute the histogram of the channel 1 reflectances of the cloudy pixels by scanning the image; -the first four classes must contain at least 250 pixels when summed up, an arbitrary number ensuring some kind of statistical significance. Otherwise reject the scene; -each of the second, third and fourth classes must contain at least 10 % of the population of cloudy pixels. This ensures that medium reflective clouds are significantly represented. Otherwise reject the scene; -the first class must not contain more than 20 % of the population of cloudy pixels, that is not too much low reflective clouds. Otherwise reject the scene; -the mean channel 1 reflectance of the cloudy pixels must not be greater than 0.7 and the lower limit of the most populated class must be less than 0.7, that is not too much bright clouds. Otherwise reject the scene; -compute at each pixel entering one of the four first classes, i. e. for which the channel reflectance is comprised between 0.4 and 0.8, the ratio of the reflectance in channel 1 to reflectance in channel 2; -the mean temperature for these same pixels should be less than 260 Kelvin. This prevents non-cloudy high reflective ocean pixels to corrupt the results. Otherwise reject the scene;

-compute the mean value of these ratios. This algorithm contains a number of limits of classes and some other numbers, which have been defined from the analyses of Asmami, Wald. It has been checked that a relative change of 10 % of these values does not induce noticeable change on the assessment of the correction factor r 12 .

RESULTS FOR VARIOUS SATELLITES

The above algorithm was tested in a number of cases, covering various areas (Northern Pacific ocean, tropical Atlantic, Gulf of Biscay, Western and Eastern parts of the Mediterranean sea, Barents sea). It was implemented in routine operations, some of them outside Ecole des Mines de Paris. These operations provide correction factors which were then sent to me. Many thanks are due to the company GM-Images which supplied many of the results. It should be noted that these routine operations were all aiming at the mapping of the sea surface temperature, and that they mostly used cloud-free scenes. This impeded a continuous assessment of the correction factor.

For NOAA-7, -9, and -11, when several computations of the correction factor r 12 have been made over a period of 50 days, they are averaged over this period before reporting in Figures 1 to 3. Otherwise, the single value is reported. For the other NOAA satellites, since the number of assessment is very limited (3 or 4), the values are reported in Tables 1 and2.

In the Figures, our approach is called the 'cloud method'. The assessments made by this method are compared to the models proposed by [START_REF] Kaufman | Calibration of the AVHRR visible and near-IR bands by atmospheric scattering, ocean glint and desert reflection[END_REF][START_REF] Rao | Post-launch calibration of the visible and near-infrared channels of the Advanced Very High Resolution Radiometer on NOAA-7, -9, and -11 spacecraft[END_REF], and to other assessments made by [START_REF] Teillet | Three methods for the absolute calibration of the NOAA AVHRR sensors in-flight[END_REF]), Justus (1989), [START_REF] Santer | Inflight calibration of channels 1 and 2 of AVHRR using desertic sites and clouds[END_REF][START_REF] Vermote | Absolute calibration of AVHRR visible and near-infrared channels using ocean and cloud views[END_REF]. It should be noted that the last three methods are also using clouds, but highly reflective ones for the last two methods.

The 'cloud method' provides assessments which are in very good agreement with the others. The correction factors are most often less than those of Vermote, Kaufman (which are most often the highest). They are very close to the model of Rao, Chen, and close to that of Kaufman, Holben, despite the fact that the latter was used by Asmami, Wald to establish that method.

The deviation of the correction factor r 12 from 1.0 is large and varies in time, often decreasing as the time since launch increases. This means on the one hand that time-series of NDVI should be corrected for the degradation of the calibration of the AVHRR instruments, and on the other hand, that periodic assessments of this correction factor r 12 should be made.

PRACTICAL IMPLEMENTATION

This method is composed of a few elementary operations and can be easily implemented. If a cloud detection procedure is applied, then these clouds can be automatically screened and when suitable the interband calibration is computed. But it will also work by using a simple threshold to detect cloud (reflectance should be greater than 0.4 over the ocean for a pixel to be declared as cloudy).

Given an AVHRR image containing five channels of raw data, the practical implementation can be the following:  convert raw digital numbers of channels 1 and 2 into albedo by means of the standard NOAA procedure using the pre-flight calibration values:  i and C 0i ;  convert albedo into reflectance  i * , by dividing the albedo by cos ( s ), where  s is the solar zenith angle. Actually, one should take into account the daily correction of the distance from the earth to the sun, but it is not necessary here;  reject non-ocean pixels;  avoid sunglint (see e.g., [START_REF] Monget | Remote sensing of the seastate using the 0.8-1.1 microns channel[END_REF] or simply use only pixels for which the viewing direction (azimuth) is close to the sun direction. In the case of an afternoon overpass of the satellite, keep only pixels located in the eastern part of the swath. This step is not mandatory, because the test on the mean temperature avoids such cases. However it may corrupt the statistics of the population of medium to high reflectances, and a scene may be rejected while it contains suitable clouds;  reject pixels which reflectance is lower than 0.4 in channel 1. The remaining pixels are cloudy pixels;  define five classes of reflectances: 0.4-0.5; 0.5-0.6; 0.6-0.7; 0.7-0.8; 0.8-0.9;  define a sixth class which lower limit is 0.4 and the higher limit is arbitrarily set to 2.0 for sake of simplicity. This class contains all the cloudy pixels and hence provides the population;  compute the histogram of the channel 1 reflectances of the cloudy pixels by scanning the image;  the first four classes must contain at least 250 pixels when summed up. Otherwise reject the scene;  each of the second, third and fourth classes must contain at least 10 % of the population of cloudy pixels. Otherwise reject the scene;  the first class must not contain more than 20 % of the population of cloudy pixels. Otherwise reject the scene;  the mean channel 1 reflectance of the cloudy pixels must not be greater than 0.7 and the lower limit of the most populated class must be less than 0.7. Otherwise reject the scene;  compute at each pixel entering one of the four first classes, the ratio of the reflectance in channel 1 to reflectance in channel 2;  compute for these same pixels the mean temperature. This temperature should be less than 260 Kelvin. Otherwise reject the scene;  compute the mean value of these ratios. This provides the quantity r 12 ;  compute the NDVI according to Equation 4.

CONCLUSION

From these results, it can be concluded that the interband calibration of AVHRR data is possible using clouds as targets. Following the procedure described above, it has been demonstrated that an accuracy (rms.) better than 5 % can be achieved in a simple way. The method is made up of elementary procedures, and has been proven to be robust. It has been implemented in routine operations, and has given reliable and accurate results without the help of an operator. The method was developed using a few images of NOAA-9 in summer 1988, and has been applied to various satellites, from TIROS-N to NOAA-11.

The deviation of the correction factor r 12 from 1.0 is large and varies in time, often decreasing as the time since launch increases. This means on the one hand that time-series of NDVI should be corrected for the degradation of the calibration of the AVHRR instruments, and on the other hand, that periodic assessments of this correction factor r 12 should be made.

Figure 1 .

 1 Figure1. Correction to the ratio of AVHRR channel 1 to channel 2 for NOAA-7. When computed with the pre-flight values, this ratio should be multiplied by this correction factor r 12 . The results of the cloud method are compared to the models of Kaufman, Holben, and Rao, Chen, and also to other results.

FigureFigure 3 .

 3 Figure 2. As Figure 1, but for NOAA-9

Table 1 .

 1 Correction values to the ratio of AVHRR channel 1 to channel 2 found for TIROS-N, and NOAA-6, for several days, using the cloud method.

	TIROS-N

Table 2 .

 2 As Table1, but for NOAA-8, and NOAA-10.