Semi-Automatic Normalization of Multitemporal Remote Images Based on Vegetative Pseudo-Invariant Features

A procedure to achieve the semi-automatic relative image normalization of multitemporal remote images of an agricultural scene called ARIN was developed using the following procedures: 1) defining the same parcel of selected vegetative pseudo-invariant features (VPIFs) in each multitemporal image; 2) extracting data concerning the VPIF spectral bands from each image; 3) calculating the correction factors (CFs) for each image band to fit each image band to the average value of the image series; and 4) obtaining the normalized images by linear transformation of each original image band through the corresponding CF. ARIN software was developed to semi-automatically perform the ARIN procedure. We have validated ARIN using seven GeoEye-1 satellite images taken over the same location in Southern Spain from early April to October 2010 at an interval of approximately 3 to 4 weeks. The following three VPIFs were chosen: citrus orchards (CIT), olive orchards (OLI) and poplar groves (POP). In the ARIN-normalized images, the range, standard deviation (s. d.) and root mean square error (RMSE) of the spectral bands and vegetation indices were considerably reduced compared to the original images, regardless of the VPIF or the combination of VPIFs selected for normalization, which demonstrates the method’s efficacy. The correlation coefficients between the CFs among VPIFs for any spectral band (and all bands overall) were calculated to be at least 0.85 and were significant at P = 0.95, indicating that the normalization procedure was comparably performed regardless of the VPIF chosen. ARIN method was designed only for agricultural and forestry landscapes where VPIFs can be identified.


Introduction
Remote sensing observations are usually instantaneous and are affected by many factors, such as atmospheric conditions, sun angle, viewing angle, dynamic changes in the soil and plantatmosphere system, and changes in the sensor calibration over time [1,2]. The goal of radiometric corrections is to remove or compensate for all of the above effects. Exceptions to this procedure include corrections for actual changes in the ground target to retrieve surface reflectance (absolute correction) or to normalize the digital counts obtained under the different conditions and to establish them on a common scale (relative correction) [2].
Absolute radiometric corrections (ARC) make it possible to relate the digital counts in satellite image data to radiance at the surface of the Earth. This relation requires sensor calibration coefficients, an atmospheric correction algorithm and related input data among other corrections [2]. A considerable amount of research has been performed to address the problem of correcting images for atmospheric effects. Radiometric normalization of remote imagery requires all of the previously mentioned information at the time of image acquisition. For most historically remote scenes, these data are not available, and for planned acquisitions, the data may be difficult to obtain [3]. Consequently, absolute surface reflectance retrieval may not always be practical [2].
Relative radiometric normalization (RRN) based on the radiometric information intrinsic to the images themselves is an alternative whenever absolute surface radiances are not required. RRN of imagery is important for many applications, such as land cover change detection, mosaicking and tracking vegetation indices over time, and supervised and unsupervised land cover classification [2,3,4,5]. Several methods have been proposed for the RRN of multitemporal images collected under different conditions at different times [2,4]. All methods operate under the assumption that the relationship between the at-sensor radiances recorded at two different times from regions of constant reflectance is spatially homogeneous and can be approximated by linear functions. The most difficult and time-consuming aspect of these methods is the determination of suitable time-invariant features that will serve as the basis for normalization [3]. Ya'allah and Saradjian [5] developed an automatic normalization method based on regression applied to unchanged pixels in urban areas. The proposed method is based on efficient selection of unchanged pixels through image difference histogram modeling using available spectral bands and calculation of relevant coefficients for dark, gray and bright pixels in each band. Yang and Lo [6] studied five methods of RRN applied to Landsat images. This method includes pseudo-invariant features, radiometric control set, image regression, no-change set determined from scattergrams, and histogram matching, all of which require the use of a reference-subject image pair. Factors that affect the performance of RRN include land-use/land-cover distribution, water-land proportion and topographic relief [6]. Ground reference data are expensive and difficult to acquire for most remotely sensed (satellite) images, and the selection of PIF is generally subjective [2]. In practice, vegetative targets of absolutely constant reflectance do not exist. Therefore, the concept of PIF is adopted with the assumption that the reflectance is constant over time [2].
Several authors have developed powerful statistical approaches to determine invariant features for the atmospheric normalization of image pairs. Hall et al. (1991) [7] and Coppin and Bauwer (2004) [8] developed radiometric rectification techniques for land cover change detection through the use of landscape elements whose reflectance values are nearly constant over time. Hall et al. (1991) [7] selected PIFs with two sets of data, bright and dark. The two sets were selected in different images by visual inspection. Du et al. (2002) [2] developed a new procedure for radiometric normalization between multitemporal images of the same area. In this method, the selection of PIF is performed statistically, and quality control and principal component analysis (PCA) are used to find linear relationships between temporal images of the same area. Several authors have proposed a change-detection technique called multivariate alteration detection (MAD), which is invariant to linear and affine scaling [3,9]. Thus, if MAD was used for change detection applications, pre-processing by linear radiometric normalization is superfluous. An iteratively re-weighted  modification of the MAD transformation (IR-MAD) established a better background of no change upon which significant changes can be examined [10][11]. Some authors have used the IR-MAD transformation for relative radiometric normalization of multitemporal images and MAD transformation for change detection [12]. Others converted digital number values to reflectance directly by relative radiometric normalization using IR-MAD [13]. Baisantry et al. (2012;[14]) performed RRN using MAD transformation and selected PIFs automatically through the Bin-Division Method. Kim et al. (2012;[15]) developed a method designed to automatically extract pseudo-invariant features for the RRN of hyperion hyperspectral images and used band-to-band linear regression. Philpot and Ansty (2013; [16]) developed an analytical formula that relates pseudo-invariant features (PIFs) to the radiometric properties of the scenes. The formula is then inverted to yield an estimate of the ratio of the transmission spectra of the two images given the path radiance for each scene and a set of invariant features. Sadeghi et al. (2013) [17] proposed an automated RRN to adjust a non-linear based on artificial neural network and unchanged pixels.
QUAC (quick atmospheric correction) and FLAASH (fast lineof-sight atmospheric analysis of spectral hypercubes) are atmospheric correction modules used in ENVI image processing software (Exelis-Visual Information Solutions, Inc. 4990 Pearl East Circle Boulder, CO 80301 USA, http://www.exelisvis.com). QUAC is an on-the-fly method for use in real-time data processing that determines parameters directly from the information contained in the scene using the observed pixel spectra [18]. FLAASH is a physics-based correction method built on MOD-TRAN4 atmospheric correction software [19]. FLAASH allows the user to define all parameters that influence atmospheric absorption and scattering, such as relative solar position, atmospheric, aerosol, and scattering models, and visibility parameters, among others. The advantage of QUAC is that an in-scene approach is easily implemented, while FLAASH uses a very diverse atmospheric ancillary parameter, and the data are therefore highly tunable by the image expert. Hu et al. (2011, [20]) compared the FLAASH and MAD normalization methods on Landsat time-series images. Other authors applied FLAASH to change detection applications [21].
To our knowledge, no information is available on relative image normalization (ARIN) of multitemporal agricultural scenes based on VPIFs or on the development of software to achieve this semiautomatically. Our specific objectives were as follows: 1) to describe the ARIN procedure and implement it in the GeoEye-1 multitemporal image series; 2) to comparatively study the selected VPIFs in relation to the ARIN method efficacy; 3) to compare ARIN-transformed multitemporal images to the original (ORI) and to the FLAASH and QUAC calibrated images; and 4) to develop semi-automatic software to normalize any set of multitemporal remote imagery by identifying VPIFs. Materials and Methods

ARIN Procedure and Software
The procedure developed for the relative normalization of multitemporal images consists of the following steps: 1) selecting one or several VPIFs; 2) defining the same parcel or parcels for each selected VPIF in each multitemporal image; 3) extracting the VPIF spectral band data for each image; 4) calculating the correction factors (CFs) for each image band to fit each band value to the average value of the image series; and 5) obtaining the normalized images by transforming each band through CF linear functions. Further information of VPIF and of vegetative variant features (VVF) will be given later in this article. Basically they coincide with permanent orchards/mature tree plantations and annual/herbaceous crops, respectively. We select the VPIF parcels at random, among many parcels of very similar characteristics available in our agricultural scene. Main steps of ARIN procedure can be achieved by conventional image processing menus including the parcel definition, the spectral band data extraction, and the image liner transformation through the estimated correction factors (CF), as will be later defined. The CF can be calculate manually or in a excel sheet once the VPIF spectral band values has been extracted.
Environment for Visualizing Images (ENVI 4.8 and ENVI 5.0, Exelis-Visual Information Solutions) software was used to visualize and process the images. Generally, mature non-deciduous tree orchards and permanent lawn green cover are eligible to be VPIFs, and at least one must be present in the scene for normalization. A parcel of the selected VPIF needs to be drawn through the ROI/SHAPE menu in one image and then moved to the rest of the image series through the VECTOR/SHAPE menu (convert the ROI to a DXF vector). The B, G, R and NIR spectral bands of the selected VPIF can be extracted for each image through the ROI/SHAPE Tool Statistical menu.
The CFs of any VPIF spectral band, for example, the G band of the image i (CF Gi ), are defined as the ratio G m /G i , where G m is the average original value of the G band in the original image (ORI), and G i is the band value of image i. Then, each band of each ORI will be transformed by applying the corresponding linear CF through the Basic Tool-Band Math menu. The normalized image will be composed of the transformed bands through the Layer Stacking menu. To semi-automatically perform the previously described steps (steps 2 through 5), the so-called ARIN software and procedure were developed [22,23]. The ARIN flowchart is shown in Figure 1     are named V1 to V7, respectively. The panchromatic image was 0.50 m pixels 21 , and the multi-spectral-image spatial resolution was 2.00 m pixels 21 , providing information on blue (B, 450-510 nm), green (G, 510-580 nm), red (R, 655-690 nm) and nearinfrared (NIR, 780-920 nm) spectral bands. The swath width was 15.2 km. The ground was predominantly flat, with an average slope grade of 2.12%. The georeferencing accuracy of the GeoEye-1 images was improved by using ground control points (GCPs) and image-to-image co-registration [24].

Land uses and selected variant and invariant
vegetative features. The LaVentilla area was surveyed approximately every 3 weeks from April to October 2010 to identify the crop of each parcel, its stage of development, and any key agricultural features. A total of 23 land uses were identified in the Geo-Eyes-1 scenes. Vegetative systems, such as alfalfa, avena, broad beans, citrus orchards, chickpeas, corn, cotton, Mediterranean forest, olive orchards, potatoes, sunflower, rapeseed, poplar groves and winter wheat, among others, were identified. Additionally, non-vegetative land uses, such as rivers, water reservoirs, paved roads, bare soil roads, and civil buildings, were also found. The phenotypes of some herbaceous vegetation parcels, such as wheat, sunflower, corn and cotton, varied considerably throughout the growing season and can thus be designated variant vegetative features (VVFs). The phenotypes of high density adult tree plantation parcels, such as citrus orchards (CIT), olive orchards (OLI), and poplar groves (POP), demonstrate much less vegetation change throughout the cropping season. Thus, they can be designated pseudo-invariant vegetative features (VPIFs). To show the differences between VVF and VPIF parcels, the spectral bands and NDVI vegetative index evolution were determined in four parcels of approximately 0.3-0.5 ha for each selected VPIF (CIT, OLI and POP) and for wheat (WHT) and corn (CRN), as representatives of the winter (autumn-sown) and summer (springsown) VVFs. 2.3 Implementation of the ARIN procedure. CIT, OLI and POP were used as VPIFs. Four parcels of approximately 0.3-0.5 ha were drawn for each VPIF using the regions of interest (ROI-VECTOR)/SHAPE menu of ENVI. The ARIN procedure provided the B, G, R and NIR spectral bands, the NDVI and the G/B vegetative indices for the original and transformed/normalized images. Normalization was achieved using a single VPIF and by validating the other two VPIFs or by using two VPIFs consecutively and validating with the other.
2.4. Absolute corrections using QUAC and FLAASH. To allow comparison with ARIN, the original V1 to V7 images were transformed using QUAC and FLAASH software. QUAC software requires the presence of dark and bright pixels in the images to serve as a basis for the implemented corrections. FLAASH was implemented by fitting to GeoEye-1 sensor specifications and to the atmospheric mid-latitude summer geographic area where the satellite images were taken. This area was aerosol rural, the highest image ground elevation was 150 m, the water column retrieval parameter was 2.92, and the scene visibility parameter was 100 km for all images except for the V7 image (140 km).
2.5. Statistical parameters. VPIF and VVF parcels spectral band and vegetation index data were subjected to analysis of variance and means were separated at the 5% level of significance by the least significance difference (LSD) test with the use of SPSS Statistical-21 software (IBM North America, New York, NY, United States). For any original or transformed image, the VPIF mean, range, standard deviation (s. d.) and root mean standard error (RMSE) ofthe band spectral values, NDVI and B/G vegetation indices were determined. The root mean square error (RMSE) of the series of images was calculated by the following equation [25]: Where n is the number of images, X i is the values of each image and X m is the average of all images. The smaller range, s. d. and RMSE of data in a given series of images indicated more uniformity among images than in the others. The correlation coefficients and the level of significance between CFs were determined across VPIFs for each image band and for all images.

Evolution of Vegetative Variant and Pseudo-invariant Features
Generally, the evolution of the spectral bands and the NDVI throughout the growing season varied to a greater extent in VVFs than in VPIFs (Figures 6, 7, 8, and 9). For example, the NIR band average, range and s. d. for OLI were 466, 240, and 81 (Figure 6a), respectively, whereas for WHT, these values were 3494, 2683 and 936, respectively (Figure 6b), for the images series. Similarly, the NDVI evolution varied significantly in VVFs, while it was relatively stable in VPIFs. For example, the NDVI mean, range and s. d. were 0.63, 0.064 and 0.024 for CIT and 0.34, 0.70 and 0.26 for WHT, respectively (Figure 6c & d). These data confirmed that the phenotype, morphology, development, and observable physical characteristics of perennial plantation are very stable throughout the agricultural season for VPIFs (Figure 9), while the opposite is true for the herbaceous cropping systems (VVFs, Figures 7 and 8). This observation has been obvious to field workers and agronomists for many years. This simple, evident finding is very important for land use classification of agricultural scenes through remote sensing.

Relative Radiometric Normalization
2.1 Using a single VPIF. The spectral band values of the original GeoEye-1 image series for the given VPIFs varied considerably among the images (Table 1). For example, the CIT B band digital values for the images V2 and V6 were 241 and 322, respectively, and the POP NIR values for the same images were Figure 9. Canopy structure of citrus (A and B) and olive (C and D) orchards vary very little throughout the year and particularly in the main growing season (April to October). Photographs A and C were taken early May and B and D late September, roughly coinciding with the V2 and V7 satellite images of this study. In photograph D can be appreciated the orange fruit but not the olives in photograph D at the distance where the photograph was taken. Due to the little changes in the trees canopy and soil surface covered the NDVI values of citrus and olive are very similar throughout the growing season (Figure 6d), and can be used as pseudo-invariant features for radiometric normalization as shown in this study. doi:10.1371/journal.pone.0091275.g009 Table 1.         (Table 1). Consequently, the vegetation indices varied considerably among the original images at any given VPIF (Table 2). For example, the CIT NDVI values varied from 0.64 at V3 to 0.81 at V7, and the POP NDVI values ranged from 0.62 at V3 to 0.81 at V2. This wide variation in the VPIF band digital counts and vegetation indices among images indicates that the radiometric normalization process is highly recommended for a comparative follow-up of cropping systems and other environmental features in any multitemporal image series.
Generally, the range and s. d. of the spectral bands and vegetation indices in the ARIN-transformed images were considerably reduced compared to those of the original images (Table 1 and 2), regardless the VPIF considered. First, it should be noted that for a single transformation (using just one VPIF), the spectral bands and vegetation indices of the VPIF taken as reference produce exactly the same value for any transformed image. This value is the average of the image series, for example, 329 and 1051 for the CIT B and NIR spectral bands, respectively. Additionally, the corresponding range and s. d. are negligible.
The ARIN process was an efficient normalization process regardless of the single VPIF or the combination of VPIFs chosen. The selection of the single VPIF used for the ARIN normalization process only slightly affected the normalization results of other VPIFs. For example, the range, s. d. and RMSE of the CIT VPIF B and NIR bands at the POP transformed images series were 40, 15 and 14 and 288, 117 and 108, respectively, whereas these values were 304, 111 and 103 and 546, 218 and 202 in the original images (Table 1). Similarly, the range, s. d. and RMSE values for the R band of the CIT VPIF in the OLI transformed images were 55, 20 and 19, compared to 208, 70 and 65 for the original images, respectively (Table 1).
2.2 Using two VPIFs consecutively. For each selected VPIF, the normalization effect caused by the ARIN procedure was also determined using two VPIFs consecutively; thus, the potential slight stationary phenotypic variation could be balanced for each VPIF. Applying the ARIN process using two VPIFs consecutively is also an effective method to normalize the multitemporal series of images. For example, after normalization with the pseudoinvariants OLI+POP, the CIT NDVI range, s. d. and RMSE of the image series were 0.08, 0.03 and 0.02, respectively, in comparison to 0.20, 0.07 and 0.07for the original images (Table 2). Similarly, the consecutive implementation of ARIN CIT+POP resulted in OLI NDVI range, s. d. and RMSE values of 0.14, 0.05 and 0.04 in comparison to 0.21, 0.10 and 0.09 for the original images, respectively.
2.3 VPIF correlation factors. The VPIF spectral band CFs used to implement the ARIN linear normalization procedure varied greatly among the spectral bands for any given image ( Figure 4) and among images for any given spectral band. Furthermore, the correlation coefficients between the CFs among the VPIFs for any spectral band and for all bands were found to be at least 0.85 and were significant at P = 0.95 or higher (Table 3). This finding also demonstrates that the ARIN normalization process is efficient regardless of the VPIF selected.

ARIN vs. QUAC and FLAASH
Generally, for the series of GeoEye-1 images studied, ARIN was more efficient than QUAC and as efficient as FLAASH, varying slightly with the VPIF selected. For example, if we consider the VPIF CIT, the NDVI s. d. of the original, OLI+POP-transf., QUAC and FLAASH images were 0.07, 0.03, 0.05 and 0.08, respectively, and the same statistics for the B/G index were 0.29, 0.06, 0.11 and 0.07 ( Abbreviations: ORI; original images; CIT, citrus orchards; OLI, Olive orchards; POP, poplars grove; -transf., transformed images; B, blue; G, Green; R, red; NIR, near infra-red; S. d., standard deviation; RMSE, Root Mean Square Error. 3 For each VPIF, spectral band and image type the data of the multitemporal images followed by the same letter are not significantly different at P$0.05. 4 For each VPIF and spectral band statistical data of image types followed by a different letter are significantly different at P$0.05. doi:10.1371/journal.pone.0091275.t001 s. d. and RMSE statistical, the results are better for ARIN CIT+ OLI-transf. (0.03 and 0.03) and FLAASH (0.03 and 0.03), followed by QUAC (0.05 and 0.04) compared to the original images (0.08 and 0.07) ( Table 2). Considering the B/G index, the ARIN-transf. procedure was also more effective than QUAC and FLAASH in any VPIF studied.
ARC relates image digital counts to radiance at the surface of the Earth and requires sensor calibration coefficients, an atmospheric correction algorithm and related input data, among other corrections [2]. For most historical images, such data are not available, and for planned acquisitions, they may be difficult to obtain [3]. Consequently, ARC retrieval is not often a practical method [2]. RRN is based on the radiometric information intrinsic to the images themselves and is an alternative whenever absolute surface radiances are not required, as in change detection applications or for supervised land cover classification [3,4].
To our knowledge, no RRN methods for multitemporal remote images using VPIF as a reference have previously been developed. The concept of PIF is adopted with the assumption that the reflectance is constant over time [2]. Moreover, any individual plant, cropping system, or vegetative feature varies with time, and therefore, there is a unanimous agreement that the reflectance is not an absolute invariant. As we have shown, most agricultural areas have vegetative parcels that change drastically throughout the growing season, such as annual herbaceous crops (VVFs), while others features, such as dense forest, permanent lawn or dense non-deciduous orchard plantations, remain comparatively invariant throughout the growing season (VPIFs). We have shown drastic differences between the selected VVFs and VPIFs in the spectral bands and vegetation index evolution. The phenotypic or morphological aspects of VPIFs are well known, and therefore, the light reflectance changes very little throughout the annual growing season, which is a key factor for the land use classification of agricultural scenes through remote sensing. Additionally, the pseudo-invariability of VPIF is the characteristic used in the ARIN procedure to normalize a set of images of a common scene.
ARIN method was developed for the radiometric normalization of multitemporal images of agricultural and forestry scenes where vegetative pseudo-invariant features (VPIFs) can be identified. In our work, we have shown that selected mature CIT (citrus orchards), OLI (olive orchards) and POP (poplar groves) in a Mediterranean landscape can be chosen efficiently as VPIFs throughout spring and summer. Generally, the ARIN normalization process efficiently produced relatively uniform data for all images, regardless of the single VPIF or the combination of VPIFs chosen. This result is likely because the range and s. d. of the spectral bands and vegetation indices in the transformed images are considerably reduced compared to those of the original images.
Regardless of the single VPIF used to estimate the band CFs for the image series transformation, the results were relatively normalized when compared to those of the original images. The VPIF spectral band correction factors (CFs) used to implement the ARIN linear normalization procedure varied greatly among spectral bands for any given image and among images for any given spectral band. Furthermore, the high and statistically significant correlation coefficients between the CFs among the VPIFs for any spectral band and for all bands suggest that the ARIN normalization process was efficient regardless of the VPIF selected.
Implementing the ARIN procedure consecutively using two VPIF was also an efficient method of normalizing the multitemporal images series. Generally, for the multitemporal series of GeoEye-1 images studied, ARIN was more efficient than QUAC and as efficient as the FLAASH absolute calibration method. The advantage of the ARIN method is that weather calibration parameters are not necessary, whereas they are required for the highly tunable FLAASH methods.
VPIF size is clearly related to the image spatial resolution and should have a sufficient number of pixels to provide a solid average of the selected vegetative feature. With medium to high spatial resolution images from the satellite GeoEye-1 (i.e., ,5 m pixel), VPIF size normally coincides with uniform parcels of approximately 0.3 to 0.5 ha or larger. Moreover, ARIN can also be used with very high spatial resolution images (i. e. 3 to 5 cm pixel) such as those provided by unmanned aerial vehicles (UAV) [26]. In UAV images a VPIF size of 2 to 4 m 2 will cover a high number of pixels to provide a solid pseudo-invariant vegetative feature sample, and in practice, it may coincide with a non-deciduous tree of 2 to 4 m 2 . UAV images can be normalized through the ARIN procedure avoiding the use of the barium sulfate standard spectralon panel, which is placed in the middle of the field to calibrate data [26].
Step-by-step implementation of the VPIF-based ARIN normalization procedure through the available ENVI imageprocessing tools is time consuming, and therefore, it is not economically feasible. ARIN software [22] quickly and easily executes ARIN and generates the transformed images; consequently, its development is essential for the practical application of the ARIN method. However ARIN method can be applied only in any agricultural and forestry landscapes where a VPIF can be identified.

Conclusions
A novel method for the radiometric normalization of multitemporal images, named ARIN, was developed to be used in agricultural and forestry scenes where a vegetative pseudoinvariant features (VPIFs) can be identified. This new procedure identifies one common VPIF parcel in all scenes, extracts the spectral band values of each image and transforms them to common band values through linear transformation. We validated ARIN using a series of GeoEye-1 satellite images of one scene. ARIN worked correctly, regardless of the three VPIFs considered (citrus orchards, olive orchards and poplar groves). The ARIN method was slightly more efficient than the absolute calibration QUAC method and as efficient as the highly tunable FLAASH method, which uses solar position and weather calibration parameters. Implementing the ARIN procedure through conventional image processing is time consuming. The software ARIN executes ARIN semi-automatically in an economically feasible manner.  B and G are spectral bands. 3 For each VPIF, vegetation index and image type the data followed by the same letter are not significantly different at P$0.05. 4 For each VPIF and vegetation index statistical data of image types followed by a different letter are significantly different at P$0.05. doi:10.1371/journal.pone.0091275.t002 Table 3. Correlation coefficients between the spectral band CFs of VPIF CIT, OLI and POP.