Band-Based Best Model Selection for Topographic Normalization of Normalized Difference Vegetation Index Map

Topographic effect in remote sensing images is severe in high mountainous areas. Efficiently to reduce the effects, several topographic normalization models have been proposed. Since the performance of the models is largely dependent on the spectral band and land surface type, the best performance model can vary from image to image in an area as well as from band to band in an image. The normalized difference vegetation index (NDVI) map has been widely used for the vegetation monitoring and assessment. An efficient reduction of the topographic effect in the NDVI map must be required for the spatial analysis of the vegetation monitoring and assessment. In this paper, we propose an efficient method to select the best topographic normalization model in each band to reduce the topographic effect of NDVI maps. The histogram structural similarity (HSSIM) index was used for the model selection because the index allows to select the best model in each band of an image. Five topographic normalization models were used for the test, which include the sun-canopy-sensor (SCS), statistical-empirical, C-correction, Minnaert, and Minnaert + SCS. The performance of the proposed method was validated by using two different season Landsat-8 OLI images including the forest area of northern Malaysia. The standard deviations of the two NDVI maps generated from the test images were reduced by about 53.1% and 28.6% after correction in profile analysis. The coefficient of determination ( $R^{2}$ ) between the two different NDVI maps increased from 0.626 to 0.759. It indicates that the proposed method effectively reduced the topographic effect of the NDVI maps. This result implies that the proposed method can work well in the topographic normalization. Furthermore, the proposed method would be successfully applied to index maps including the normalized difference snow index (NDSI), normalized difference water index (NDWI), etc.


I. INTRODUCTION
Remotely sensed images have been widely used for the vegetation monitoring and assessment. However, it is not easy to monitor and assess the vegetation changes in mountainous areas due to severe topographic effect [1]- [6]. For example, the reflectance values in the sun-shaded slope area are usually lower than the sunlit slope area, even though the sun-shaded and sunlit areas are covered by the same vegetation type [2]. This effect may have a significant influence on the vegetation The associate editor coordinating the review of this manuscript and approving it for publication was Md. Moinul Hossain . cover identification results [5], [6]. Thus, the topographic normalization process is required to improve vegetation cover identification in mountainous area with high topographic variations [3], [4].
Topographic normalization is defined as the normalization of reflectance values on a slanted plane to a horizontal plane [7]- [9]. Several topographic normalization models have been developed to reduce the topographic effects. The models have been proposed by considering the surface reflection types (Lambertian or non-Lambertian), the slope and aspect of mountainous areas, the look angle, and the solar zenith angles [10]- [13]. Among the proposed models, the Sun-canopy-sensor (SCS) [10], Statistical-empirical [11], C-correction [11], Minnaert [12] and Minnaert + SCS [13] have been widely and successfully used for the topographic normalization.
Since the model performance largely depends on the spectral band and land surface type, the best performance model cannot be same in all images as well as it can be different even in each image band. Thus, several performance evaluation methods have been proposed to evaluate the model performance and find the best model among the topographic normalization models, but the methods have a disadvantage that they cannot determine the best model in each image band in most evaluation methods [14].
Recently, the evaluation methods using the structural similarity (SSIM) index [15] and the histogram structural similarity (HSSIM) index [14] have been proposed and they enable us to evaluate the model performance in each image band. To evaluate topographic normalization methods, we need to determine how similar the pixel values are in the sunlit and sun-shaded areas. To utilize the SSIM index, we must need to find pixels in a sunlit area that correspond to pixels in a sunshaded area. However, the corresponding pixels cannot be found in the sun-shaded area. Thus, the SSIM index requires a synthetic image that does not have topographic effect at all. It is almost difficult to meet the requirement for the SSIM evaluation, although the SSIM evaluation can provide a more precise statistical value. On the other hand, the HSSIM evaluation does not require the topography-free synthetic image, because it estimates from the histogram structural similarity between the sun-shaded and sunlit slopes in an image [15]. The HSSIM index has an advantage that it can also consider the topographic slope effect. However, the HSSIM evaluation needs an assumption which the sun-shaded and sunlit slopes have a similar condition. Since similar land covers can be found in both the sun-shaded and sunlit slopes, it would not be problematic in most cases.
The normalized difference vegetation index (NDVI) map has been widely used for the vegetation monitoring and assessment. Although the topographic effect of the NDVI map is much lower than the reflectance image, it needs be reduced for the forest type monitoring and classification. Several methods have been proposed to correct the topographic effects of the NDVI map and they has also applied to other ratio maps such as leaf area index (LAI), normalized burn ratio (NBR) and normalized difference snow index (NDSI) maps [2], [3], [16], [17]. In the researches, the ratio maps such as NDVI, LAI, NBR, and NDSI were generated by using one among the models such as the Cosine, C-correction and Minnaert. However, as aforementioned, since the model performance varies according to spectral band [1], [14], [15], [18], [19], the best model can be different in each band image. For example, Sola et al. [15] reported that the Minnaert + SCS model was suitable for the normalization of the nearinfrared (NIR) and short-wavelength infrared (SWIR) bands while the statistical-empirical model was suitable for the red, green and blue bands in their study area. Consequently, the best model needs to be differently selected according to spectral band for more precise topographic normalization of the ratio maps.
In this study, we propose an efficient band-based best model selection method for the topographic normalization of NDVI maps. A key of the proposed method is to determine the best normalization model in each band using the HSSIM evaluation approach. Two Landsat-8 operational land imager (OLI) images were used to validate the performance of the proposed method, which were obtained in the natural vegetation area of northern Malaysia. Since the two images were acquired in different two seasons, they had different solar zenith and azimuth angles. From the test data, the reduction of the standard deviation in the sun-shaded and sunlit areas was checked by using the original and normalized NDVI maps and the similarity of the two different season NDVI maps was compared before and after the topographic normalization as well.

II. STUDY AREA AND DATASET
As shown in Fig. 1a, Landsat-8 OLI images in different two seasons were used for this study. The images were obtained with different solar zenith and azimuth angles. The red and NIR images were selected for the NDVI map generation, and used for the performance evaluation of the proposed method. The atmospheric correction was first applied to the two images by using the cosine of the solar zenith angle (COST) model [20], and then the slope and aspect angles of the ground surface were calculated from the space shuttle radar VOLUME 8, 2020 topography mission (SRTM) digital elevation model (DEM) having the spatial resolution of 30 m. From the topographic slope and aspect angles, the local incidence angle were estimated [21], [22]. Fig. 1b shows the shaded relief topographic map created by the SRTM DEM.
The Landsat-8 images were obtained on December 12, 2014 (hereafter test image 1) and March 18, 2015 (hereafter test image 2). The temporal difference between the two images is 96 days. The sun zenith and azimuth angles of the test image 1 is about 37.4 • and 141.3 • , and the test image 2 is about 29.0 • and 102.9 • , respectively. The world geodetic system (WGS84) for geodetic datum and North-UTM47 for map projection were used in this study. The study area is located nearby Sik city in northern Malaysia, which is between 5 • 40'N to 5 • 60'N and 100 • 39'E to 100 • 55'E ( Fig. 1a). The mean altitude was approximately 234 m, ranging from 23 m to 1,006 m. The maximum and mean slope angles calculated from DEM were about 64.6 • and 13.7 • , respectively. The maximum and mean incidence angles were about 92.4 • and 38.8 • in the test image 1, respectively, and the maximum and mean incidence angles were about 82.3 • and 31.4 • in the test image 2, respectively. The images in the study area were cloud free. The study area consists mainly of hill dipterocarp forests, in which the trees of the family Dipterocarpaceae are dominant species [23]. The area consisted of a similar forest type, and hence it was suitable for the performance test of the best model selection method.
The sun-shaded slope areas can be easily distinguished because the reflectance values are lower than the sunlit slope areas due to topographic effects, as seen in Figs. 1c and 1d. Moreover, as seen in Figs. 1c and 1d, the contrast of the two slopes is more pronounced in the test image 1 due to higher incidence angle than in the test image 2. The test image 2 has a lower contrast and more emphasis on red color because the acquisition date was close to summer season. It also indicates that the study area is suitable to the performance test of the proposed method.
To evaluate the topographic normalization performance, the sample data of 99,000 pixels in the sun-shaded slopes and 87,000 pixels in the sunlit slopes were extracted from the test image 1. And the sample data of 106,000 pixels and 84,000 pixels were also extracted from the sun-shaded and sunlit slopes of the test image 2, respectively. The sample data were used to calculate the HSSIM index values in sun-shaded and sunlit slope areas in the test image 1 and 2.

III. METHODS
In this paper, an efficient approach is suggested to select band-based best model for the topographic normalization of NDVI maps. To validate the performance of the suggested model selection, the following steps are applied: 1) Red and NIR reflectance images are respectively normalized by using several topographic normalization models; 2) Best normalized red and NIR images are respectively selected by using the HSSIM evaluation; 3) Topography-normalized NDVI map is created from the best normalized red and NIR images. Five topographic normalization models in this study were SCS, Statistical-empirical, C-correction, Minnaert, and Minnaert + SCS.

A. TOPOGRAPHIC NORMALIZATION
From the early 1980s to the present, many topographic normalization models have been proposed to normalize topographic effects in mountainous areas. In this study, we just selected five of the developed topographic normalization models, as listed in Table 1, which are known to have a good performance. The models include SCS, Statistical-empirical, C-correction, Minnaert and Minnaert + SCS [10], [12], [13]. The selected models can be categorized according to surface reflection types and normalization approach as Lambertian, non-Lambertian and empirical approaches [24].
In the C-correction model, the constant C is calculated based on the fitting relationship between the reflectance values and the cosine of the solar incidence angle [11].
where a and b are the model parameters estimated by the linear regression fitting. The constant C is calculated as: In the Minnaert model, the constant k is calculated by using the log function as shown by Smith et al. [12]: where ln (cos (i) × cos (θ n )), ln (ρ cos (θ n )) and ln (ρ h ) can be replaced with x, y and m, respectively. Hence, this equation can be expressed as y = kx + m. Then, the constant k can be calculated by the linear regression method. The constant k for the Minnaert + SCS model can be calculated in the same way as constant k [13].

B. PERFORMANCE EVALUATION
The performance evaluation of normalization models can be carried out by the HSSIM index [15]. The HSSIM index can precisely evaluate the performance of the models. This index quantitatively calculates the histogram similarities between the sunlit and sun-shaded slopes for each band before and after normalization. The histograms similarities include the effect of the topographic slope gradients. The index is defined by [15]: where x and y are reflectance values in the sunlit and sunshaded slopes, respectively. V (x, y) and R(x, y) are the variation ratio and histogram structural similarity ratio between the original and normalized images, respectively, in the x and y data, as given by [15]: where σ x 0 , σ y 0 , σ x and σ y are the standard deviations in the x and y data of the original and normalized images, respectively. r H (x 0 )H (y 0 ) and r H (x)H (y) are the correlation coefficient between the x and y histogram data in the original and normalized images, respectively. If the topographic normalization is successful, 1) the standard deviation of the normalized image must be smaller than that of the original image, and 2) the correlation coefficient between the sunlit and sun-shaded slopes must be larger than the original image [15]. Thus, if the topographic effect is perfectly normalized, both, V (x, y) and R(x, y), are close to zero, and hence the HSSIM index is close to 0. On the other hand, if both, V (x, y) and R(x, y), show over-corrected values (V (x, y) > 1 and R(x, y) > 1), then the HSSIM index is unconditionally larger than 1 [15]. Thus, the HSSIM index allows to evaluate the normalization performance of each normalization model.

C. TOPOGRAPHIC NORMALIZED NDVI
As aforementioned, the best normalization model for each band can be different. The optimal normalization model can vary depending on the surface cover types and spectral wavelength bands [1], [14], [15], [18], [19]. Thus, for an example, the best topography-normalized NDVI image can be generated by using the red image corrected by the Statisticalempirical model and the NIR image by the Minnaert + SCS model. The best corrected NDVI image would be much more useful in identifying and classifying land cover type because the topographic effects can be minimized. The best topography-normalized NDVI map (NDVI TP ) can written as follows [2], [19]: whereρ NIR andρ RED are the NIR and RED reflectance images, which are topography-normalized by the best optimal models for NIR and RED images, respectively.  Fig. 1a. Figs. 2b to 2f show the topography-normalized images processed by using the SCS, statistical-empirical, C-correction, Minnaert, and Minnaert + SCS models. The model parameters of the topographic normalization models were estimated by using pixels in forest area. The result showed that most normalization models, except the SCS model, mitigated the topographic effects well. The Statistical-empirical, C-correction, Minnaert and Minnaert + SCS models showed that the topographic effects were normalized adequately. The results are very similar visually too (Fig. 2c-f). However, the SCS model showed severely over-corrected results due to the limitation of the Lambertian approach (Fig. 2b). Many research papers have reported 1) that the SCS model has a severe overcorrection problem [19], [25] and 2) that non-Lambertian and empirical approaches are better than the Lambertian approach [1], [18], [19]. It means that the result is in overall agreement with the researches [1], [2], [18], [19], [25]- [27]. Fig. 3 shows the histogram variations of the test image 1 before and after applying the normalization models for each band. The grey-and white-colored histograms represent the distribution of reflectance values in sun-shaded and sunlit slopes, respectively. The histograms have been successfully used to evaluate the normalized processing result [15]. In the whole bands of the original image, the sun-shaded and sunlit slopes are easily distinguishable because the reflectance values of sun-shaded slope are lower than those of the sunlit slope. If the performance of a normalization model is almost perfect, the difference between the two histograms should be close to zero. The result from the statistical-empirical model in the NIR band is representative in this study. It is difficult to find the difference between both the histograms so that they show a greater similarity than the other models. Most models with good performance in visual analysis made the sunlit and sun-shade histograms similar although the two histograms did not match perfectly. On the other hand, the SCS model showed overcorrected results in the RGB bands. The sun-shaded histogram had a higher reflectance value than the sunlit slope. It means that the similarity between the sunlit and sun-shaded histograms is in good agreement with the normalization performance. Table 2 summarizes the HSSIM index values calculated from before and after the topographic normalization in the test image 1 and 2. The HSSIM indexes of the SCS model had a value of more than 1 to indicate the over-corrected results in RGB bands. Most models have an index of less than 1 in each band wavelength. In particular, the normalization performance of the Statistical-empirical model in the NIR band had an index value of about 0, indicating almost perfect normalization. We can see that the normalization performance of all the models depends on the wavelength because the model with the smallest index is different for each  band. The models with the best normalization performance in the Blue to SWIR-2 bands are C-correction, Minnaert + SCS, C-correction, Statistical-empirical, C-correction and Minnaert + SCS, respectively, in the test image 1 and 2 as listed in Table 2.
From the evaluation result, the best topography-normalized image was generated by combining the best models in each band. For example, to create a best normalized image in case of the test image 1, the C-correction model can be used for the blue, red and SWIR-1 bands, the Minnaert + SCS model can be utilized to the green and SWIR-2 bands, and the Statistical-empirical model can be used for the NIR band. That is, the C-correction, Statistical-empirical and Minnaert + SCS models can be selected according to the performance evaluation to make the best topography-normalized image. Fig. 4 presents the NDVI maps of the test image 1 and 2. The original NDVI maps in the test image 1 and 2 are respectively shown in Figs. 4a and 4d, the topographynormalized NDVI maps with the statistical-empirical model VOLUME 8, 2020 are presented in Figs. 4b and 4e, respectively, and the topography-normalized NDVI maps with the proposed best selection model are presented in Figs. 4c and 4f, respectively. All the NDVI maps represent the characteristics of NDVI well, which has a high index value in the vegetation area and a low index value in the water and urban areas [6]. However, the NDVI maps calculated from the original image were different from the other topography-normalized NDVI maps. As seen in Figs. 4a and 4e, the original NDVI values in mountainous areas were severely affected by the topography. Even though the topographic effect was reduced by using the normalized difference index, the relationship between NDVI value and the topographic aspect and slope was very clear in Figs. 4a and 4e. The separation between the forest area and the water and urban areas can be easily applied to the original NDVI map regardless of the topographic effect, but the topographic effect makes it difficult to divide the forest types. This is the reason why the topographic normalization process is required in an NDVI map.
As shown in Figs. 4b and 4e, the topographic effect in the NDVI maps, which were topography-normalized by the Statistical-empirical model, was mitigated well, and hence the NDVI values in the sunlit and sun-shaded areas became similar. However, the normalized NDVI values still had the weak topographic effect. Especially for the test image 2, there was a problem of overcorrecting on the mountain ridges as shown in Fig. 4e. The NDVI maps topography-normalized by the proposed model selection method were remarkably improved (see Figs. 4c and 4f). The topographic effect was not severe in the NDVI maps, and the overcorrection problem was not found in the NDVI map of Fig. 4e. Fig. 5 shows the NDVI maps of test image 1 and 2 zoomed from the box in Fig. 1a to compare the topographic normalization performance between the original, Statisticalempirical normalized, and best model normalized NDVI maps. As seen in the original NDVI maps (Figs. 5a and 5e), most of the sun-shaded slope area had a relatively lower NDVI value than the sunlit slope area.
As presented in the normalized NDVI maps with the Statistical-empirical model, the topographic effect was remarkably reduced when compared to the original NDVI maps (Figs. 5b and 5e). However, some topographic effects on the sun-shaded slopes remained in the NDVI maps (see the blue arrow in Figs. 5b and 5e). This is because the Statisticalempirical model performed well in the NIR band while it was not carried out completely in the red band. The HSSIM index in the red band was two times larger than the result of the optimal model, i.e. the C-correction. Moreover, as seen in Fig. 5e, overcorrected NDVI values were found in sunlit areas.
The NDVI maps normalized by the proposed method were shown in Figs. 5c and 5f. When the maps were visually compared with the original and Statistical-empirical normalized maps, it is clear that the NDVI maps of Figs. 5c and 5f were the best topography-normalized NDVI maps. The overcorrection problem in the Statistical-empirical model was not found in the proposed method as well as the study area looked like a flat surface. This result indicates that the proposed method has a remarkable performance in normalizing the topographic effect of an NDVI map. Fig. 6a and 6b plot the NDVI values of the test image 1 and 2 in the profile shown in Fig. 5a, respectively. The profiles include both the sun-shaded and sunlit slope areas. As seen in Fig. 6a, the original NDVI profile of the test image 1 varied from 0.66 to 0.78 according to the slope angle, although the study area consists of the same vegetation cover. The NDVI values from the Statistical-empirical model showed a smaller range than the original NDVI. The NDVI values from the proposed method ranged from 0.70 to 0.75. The NDVI values of the test image 2 had the ranges of 0.71 to 0.78, 0.72 to 0.78, and 0.74 to 0.77 in the original, Statistical-empirical and proposed NDVI maps, respectively. Moreover, we could find the overcorrection problem in the Statistical-empirical NDVI map of the test image 2 (blue dot line in Fig. 6b). The standard deviations were calculated from the profiles of the test image 1 and 2. The standard deviations of the proposed NDVI maps were remarkably reduced from 0.032 to 0.015 and from 0.021 to 0.015 in the test image 1 and 2, respectively, which correspond to the reduction rates of about 53.1% and 28.6%, respectively. This result shows that the proposed method enables us to remarkably improve the topographic normalization performance of the NDVI maps. Fig. 7 shows the scattergrams of NDVI values between the test image 1 and 2. The data used for the scattergrams was obtained in mountainous areas with altitudes above 100m. The coefficient of determination (R 2 ) between the original NDVI maps of the test image 1 and 2 was about 0.626 and the root mean square difference (RMSD) was about 0.041 (Fig. 7a). As seen in Fig. 7a, the distribution of the NDVI values was larger than the other scattergrams of Figs. 7b and 7c. Since the two test images were obtained at different solar zenith and azimuth angles, the topographic effect appearing on each image was different. When the NDVI values were topography-normalized by using the Statisticalempirical model, R 2 and RMSD were largely improved to 0.711 and 0.036, respectively. This result shows that the topographic normalization processing must be applied to the NDVI maps. And it is also consistent with the reason why many researchers have applied the topographic normalization model to generate an improved index map [2], [3], [16], [17].
Moreover, when the proposed method was applied to the NDVI maps, the topographic normalization performance was more improved (Fig. 7c). Our result showed R 2 of about 0.759 and the RMSD of about 0.031. This result shows the highest correlation among all scattergrams. Most of all, the distribution of the scatter diagram from the proposed method is much closer to a perfect ellipse. It shows that the proposed method properly normalizes the topographic effect.
Further performance validation was performed by using the correlation coefficient between cosi and NDVI. The surface reflectance has a linear relationship with cosi, and hence the correlation coefficient (r) between cosi and surface reflectance enables us to estimate the performance of topographic normalization [28]. If the topographic normalization was well applied, r will be close to zero. In this study, the test was performed by using pixels acquired in mountainous forest areas. The r in the test image 1 were about 0.575, 0.063 and 0.032 in the original, Statistical-empirical, and proposed NDVI maps, respectively, and the r in the test image 2 were about 0.303, 0.095 and 0.031 in the original, Statistical-empirical, and proposed NDVI maps, respectively. The r values of the proposed NDVI maps were about 2 to 3 times closer to zero than the Statistical-empirical NDVI maps. Consequently, this result further validates that the VOLUME 8, 2020 proposed method normalizes the topographic effect of the NDVI maps well.
In this study, the HSSIM index was used for evaluating the topographic normalization performance. However, it could be replaced by another appropriate evaluation method. In addition, five topographic normalization models were used in this study, but other models can be used for the best model selection. If an existing DEM having a higher spatial resolution than the multi-spectral (MS) bands of Landsat OLI images can be used in the proposed method, the performance of topographic normalization would be improved [18], [29]. Finally, it should be noted that the proposed method can be applied to the topographic normalization of other index maps such as NDSI, LAI, NBR, etc.

IV. CONCLUSION
In this paper, we proposed an efficient method to generate a topography-normalized NDVI map. The key of the proposed method is to create the topography-normalized NDVI map after selecting and applying the best topographic normalization model for each band. For this, the HSSIM index is used to evaluate the performance of the five topographic normalization models such as the SCS, Statisticalempirical, C-correction, Minnaert and Minnaert + SCS. The HSSIM index enables us quantitatively to evaluate the topographic normalization performance based on the similarity between sun-shaded and sunlit slope areas. Two Landsat OLI images acquired from the forest area of northern Malaysia at two different seasons were used to validate the proposed method.
The models were evaluated by using the HSSIM index. Consequently, for the two test images, the C-correction, Minnaert + SCS, C-correction, Statistical-empirical, C-correction and Minnaert + SCS were selected as the best models in the Blue, Green, Red, NIR, SWIR-1, and SWIR-2 bands, respectively. To validate the performance of the proposed method, three NDVI maps were created from the original image, topography-normalized image by Statistical-empirical model and topography-normalized image by the proposed method. Although the topographic effect was remarkably reduced in the Statistical-empirical NDVI map, some topographic effects on the sun-shaded slopes remained in the NDVI maps as well as the overcorrected NDVI values were found in sunlit areas. The overcorrection problem was also found in the profile. The improvement by the Statistical-empirical model were clear.
Further improvement could be found in the proposed method. The undercorrection and overcorrection problems in the Statistical-empirical NDVI map were not shown in the proposed NDVI map. The scattergrams of the proposed NDVI maps created from the two test data showed that R 2 increased from 0.626 to 0.759 and RMSD reduced from 0.041 to 0.031. Most of all, the distribution of the scatterdiagram from the proposed method is much closer to a perfect ellipse. The results show that the proposed method has a higher performance.