A Correction Method of NDVI Topographic Shadow Effect for Rugged Terrain

The influence of solar altitude and azimuth angles makes shadows prevalent in remote sensing images of rugged terrains. Consequently, the normalized difference vegetation index (NDVI) of the shadow areas is much lower than that of sunlit areas—a phenomenon known as the NDVI topographic shadow effect. In this study, we developed an NDVI topographic shadow effect correction (NTSEC) model. The NTSEC is based on the difference in solar radiation between the sunlit and shadow areas, and introduces a variable factor that indicates the intensity of shadows to simulate the reflectance for direct solar light not received in the shadow areas. The simulated reflectance for direct solar light in the shadow areas is used as a compensation value to be summed with the original reflectance for calculating the corrected NDVI. Landsat 8 Operational Land Imager images from two different regions were used to test the NTSEC method and multiple strategies were employed to provide objective evaluative results. The model performance was compared with four commonly used topographic correction (TC) models: C-correction, Minnaert, variable empirical coefficient algorithm, and statistical empirical. The results show that NTSEC can suppress the overcorrection of the self shadow areas that are difficult to eliminate in the TC methods, while simultaneously improving the undercorrection of the cast shadow areas. The corrected NDVI of the self shadow and cast shadow areas were almost the same. The NDVI differences between shady and sunny slopes were significantly reduced after correction. In addition, the NTSEC method did not produce outliers, and the NDVI of different land cover types retained interclass stability. In summary, the NTSEC model can be used as a simple and robust method for correcting the NDVI in shadow areas.


I. INTRODUCTION
R EMOTE sensing images in mountainous regions cause differences in target radiation characteristics due to terrain. Surfaces facing the sun usually receive more radiation than those facing away from the sun [1]. This condition is called Manuscript  the topographic effect and causes the sunlit and shadow areas to appear differently on remote sensing images, where the sunlit areas consist of the surface facing the sun, and the shadow areas consist of the surface turned away from the sun and the shadow of the mountain cast on the surface of the adjacent areas [2]. The presence of shadows leads to differences in the spectral response of the same land cover, which consequently leads to large differences in the vegetation index values of the same land cover in the same area and affects the inversion results of vegetation biophysical and biochemical parameters [3]- [5].
As the spatial resolution of remote sensing images increases, so do the vegetation index errors due to topographic effect [6]. The commonly used vegetation indexes include the ratio vegetation index (RVI), normalized difference vegetation index (NDVI), soil-adjusted vegetation index (SAVI), enhanced vegetation index (EVI), and perpendicular vegetation index (PVI) [7]- [11]. These indexes can be categorized as vegetation indexes in complete band ratio format (RVI and NDVI) and vegetation indexes in noncomplete band ratio format with a constant additive term in the formula (SAVI, EVI, and PVI) according to their construction form. The "ratio" properties enable the vegetation indexes in complete band ratio format to eliminate some of the influences of topographic shadows [12]. However, the elimination of these influences is incomplete, especially in mountainous areas with large terrain undulations [13]- [15]. Therefore, it is necessary to further eliminate the influence of shadows on vegetation indexes. Due to the advantage of the "ratio" property of RVI and NDVI, the correction of the topographic shadow effect for these two indexes is less difficult than that for SAVI, EVI, and PVI. In addition, NDVI was adopted by the moderate resolution imaging spectroradiometer (MODIS) land discipline group as the global-based vegetation index for monitoring the earth's terrestrial photosynthetic vegetation activity, and the audience of NDVI is higher than that of RVI [12], [16]. Therefore, in this article, we developed a topographic shadow effect correction model for NDVI.
The methods for eliminating the effect of shadows on NDVI can be divided into two categories. The first category involves methods wherein the shadow part is detected and compensated for through brightness information, such as digital number (DN) values [17], [18]. Liu and Yamazaki [19] used an object-oriented method to detect and compensate the DN values of the shadow areas. The principle of this method is based on repairing the spectral information of the shadow areas rather than directly correcting their NDVI; the restored information is subsequently used to calculate and correct the NDVI of the shadow areas This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ [20]. This correction is a relatively complicated process with a large calculation volume, and the calculation results of NDVI are unstable due to the different shadow restoration effects of different land cover types. The second category is the topographic correction model, which is the most commonly used method. Well-known topographic correction models include the C-correction (CC) [21], sun-canopy-sensor (SCS) [22], SCS with C (SCS+C) [23], statistical-empirical (SE) [21], variable empirical coefficient algorithm (VECA) [1], and Minnaert (MIN) models [24]. Under the correction method based on the topographic correction model the NDVI is indirectly corrected and the spectral information of each band needs to be corrected before NDVI calculation. This method has the advantage of clear physical meaning and a relatively simple process compared with the first category of methods [25]- [28]. However, when the solar altitude angle is low or the study area is located on steep mountainous terrain, the TC method, which is constructed by relying on the cosine of the solar incidence angle can easily cause overcorrection [29], [30].
In recent years, studies have used an approach where they start with the vegetation index and construct a series of vegetation indexes that consider the influence of shadows, such as the topography-adjusted vegetation index (TAVI) [31] and shadoweliminated vegetation index (SEVI) [32]. In this method, the vegetation index values in shadow areas are corrected by constructing a new vegetation index. However, the sensitivity of the new vegetation index to the vegetation information in each scene and the stability of the index need to be verified through a large number of experiments [33]. In addition, SEVI and TAVI require the selection of samples from both sunny and shady slopes in the study area to calculate the adjustment factor, which increases the complexity of practical application. Moreover, these are empirical indexes and their physical meaning is not clear. Therefore, this method should be based on a physical model to construct a vegetation index corrected for the topographic shadow effect so that it has a complete physical meaning. Avoiding the selection of samples to determine the appropriate index calculation formula is also an approach that can be considered.
Thus, the construction of a vegetation index correction model for shadow areas with clear physical meaning based on an existing vegetation index is a good solution to the vegetation index topographic shadow effect. In this article, we propose a correction method for the NDVI topographic shadow effect that can correct the NDVI distortion problem caused by shadows. This method and the topographic correction model based on the digital elevation model (DEM) were applied to Landsat 8 Operational Land Imager (OLI) images with severe shadow coverage and evaluated through visual inspection and statistical and correlation analysis to verify its accuracy for the NDVI correction of shadow areas.

A. Study Areas
Two areas were analyzed in this article. The first study area was located in the Gaoligong Mountains between Longyang District and Tengchong City, Baoshan City, Western Yunnan Province, China (see Fig. 1), extending between latitude 25°3 70 -25°18 43 and longitude 98°32 33 -98°52 54 . The terrain elevation of the study area ranged from 588-3578 m, with a mean elevation of 1788 m. The large differences in elevation mean that there was a large topographic slope in the study area. The maximum gradient reached 70.25°with a mean gradient of 19.79°. The land cover types in the study area were varied, with vegetation, bare land, impervious surface, and water being the main types.
The second study area was located in the Helan Mountains of the Ningxia Hui Autonomous Region, China, extending between latitude 38°40 2 -38°54 29 and longitude 105°48 35 -106°7 11 . The land cover type was mainly bare ground and vegetation, with a mean elevation exceeding 2000 m.
The Gaoligong Mountains were selected as one of the study areas because of their high vegetation cover and various types of land cover, which are conducive to analyzing the NDVI correction effect for various land cover types. In addition, the study area was located in a low-latitude mountainous area where the land cover categories show little seasonal variation. Therefore, images with a high solar altitude angle can be used as verification images. The Helan Mountains are located at a higher latitude, and the images we collected had a low solar altitude angle with a clear topographical shadow effect. There was also less vegetation cover in the area, which is crucial when evaluating the correction performance of the model for nonvegetation areas with extensive shadow cover.

B. Data
The Landsat 8 OLI multispectral images of the Gaoligong Mountain were acquired on January 22 and April 12, 2014, when the solar altitude and azimuth angles were 39.26°and 148.12°f or January and 62.60°and 121.53°for April, respectively. The Landsat 8 OLI image for Helan Mountain was acquired on November 21, 2013, when the solar altitude and azimuth angles were 29.21°and 162.39°, respectively. In addition to the Landsat 8 data, topographic data for the two study areas were obtained from the ASTER GDEM. These data were obtained from the Geospatial Data Cloud (http://www.gscloud.cn). In addition, the Landsat data were converted from DN values to top of atmosphere reflectance by radiometric calibration and atmospheric disturbances were removed by dark object subtraction before the NTSEC was applied.
The reference data for the Gaoligong Mountain study area also include the land-use/land cover information, which was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences. Based on this map, the land cover types of the study area were classified into eight broad classes. The percentage of surface covered for each land cover class is shown in Table I.

A. Theoretical Analysis
In rugged terrain, there are significant differences in the amount of solar radiation received by the ground due to the  different incidence angles of sunlight. The spatial distribution of solar radiation fields in rugged terrain has become extremely complicated due to the influence of factors such as mutual shelter between terrains, different slopes, different slope directions, and various types of underlying surfaces [34]. The total solar radiation received by the surface of rugged terrain is composed of three parts: 1) direct solar radiation, 2) sky diffuse radiation, and 3) terrain-reflected radiation [29], [30], [35]. It can be expressed as where E is the total solar radiation that a pixel receives, E d represents the direct solar radiation received by a pixel, E f is the sky isotropic diffuse radiation the pixel receives, and E a is the terrain-reflected radiation. For terrain surfaces, E d is calculated as follows [35]: where Θ is the binary coefficient, for pixels in shadow, Θ is equal to 0, otherwise Θ is equal to 1; i is the solar incidence angle; θ z is the solar zenith angle; and E h d is the direct solar radiation incident on a horizontal surface of the same height as the slope, calculated using the following formula: where E t d represents the direct solar radiation penetrating through the atmosphere, which can be computed by the second simulation of the satellite signal in the solar spectrum (6S) model. For the Gaoligong and Helan Mountain study areas, the specific parameters inputted into the 6S model were as follows: the atmospheric model was set in mid-latitude summer and mid-latitude winter, respectively; the aerosol type was continental aerosols; the aerosol concentration was 0.12 and 0.17 (from MODIS aerosol products), respectively; and the altitude of the target was 1.78 and 2 km (obtained from DEM data), respectively. cosi is the cosine of the solar incidence angle, calculated using the following equation: where θ s is the terrain slope, ϕ is the solar azimuth angle, and β is the aspect angle of the slope. Furthermore, if sky anisotropy scattering is disregarded, and we instead assume that the sky scattering is isotropic, the E f can be calculated using the following equation [36]: where E h f is the diffuse radiation on the horizontal surface, which can also be calculated using the 6S model. V d is the sky-view factor and can be approximated using the following equation [37]: The terrain-reflected radiation is complex and, thus, difficult to calculate; to simplify the complexity of the calculation, we used an approximate method. In the case of ignoring the multiple reflection between the earth's surface and the atmosphere, E a can be calculated using the following equation [35]: where E h is total radiation on a horizontal surface, C t is the terrain configuration factor and is equal to 1-V d , and ρ a is the average reflectance of the adjacent pixels.
In accordance with the variability of solar radiation received by the ground, the terrain surface can be divided into sunlit and shadow areas. Owing to the lack of direct solar radiation (E d ) in the shadow areas, the total solar radiation received in this area includes sky diffuse radiation (E f ) and terrain-reflected radiation (E a ). Thus, the solar radiation received in the sunlit areas E sun and the total solar radiation in shadow areas E shw can be expressed as We assumed that the land surface was Lambertian. Thus, it can be considered that the reflectance is an inherent property of the vegetation and does not vary with the received radiation energy. Based on this assumption, the radiance of vegetation in the sunlit and shadow areas can be calculated using the following equations: where L sun is the radiance of sunlit areas, L shw is the radiance of shadow areas, and ρ is the intrinsic reflectance.
The difference between the sunlit and shadow areas in the image because of the difference in the received solar radiation energy, which is not considered in the common image processing steps, can lead to significant differences when converting the pixel reflectance data. Hence, we simulated the reflectance of the vegetation in the sunlit areas (ρ sun ) and the reflectance of the vegetation in the shadow areas (ρ shw ), as follows: Under ideal conditions, the simulated reflectance for the same scene in the pixels should be the same in the sunlit and shadow areas. The difference in reflectance between the sunlit and shadow areas is that the reflectance for direct light is missing in the shadow areas. Therefore, we can calculate the difference between the simulated reflectance of the sunlit and shadow areas (ρ dd ) by using ρ sun and ρ shw , as shown the following equation: From (2), (3), and (8)-(11), we can obtain From (12), it can be seen that the difference between the simulated reflectance of the sunlit and shadow areas is largest when the sunlight is irradiated perpendicular to the terrain surface (cosi = 1). Therefore, the difference between the simulated reflectance of the sunlit and shadow areas when the sunlight is irradiated perpendicular to the terrain surface (ρ v dd ) is Therefore, the reflectance compensation for direct light in the shadow areas cannot exceed ρ v dd , otherwise overcorrection will occur. Therefore, we constructed the reflectance compensation formula for direct light of the shadow areas based on ρ v dd .
It is assumed that the reflectance compensation for direct light in the shadow areas varies with the intensity of the shadows, i.e., the greater the intensity of the shadows, the more direct light reflectance compensation can be obtained in the area. We can obtain the reflectance compensation for direct light in the shadow areas by adding a constraint variable α to (13) to control the amount of direct light reflectance compensation whereρ shw dc is the reflectance compensation for the direct light of the shadow areas.
The constraint variable α in (14) needs to meet the following requirements: 1) the stronger the shadows, the larger the value of α to ensure that the stronger shadows of the areas obtain more direct light reflectance compensation; 2) the range of α is 0-1 to ensure that the range of ρ shw dc is consistent with ρ dd . Due to the advantage of its band ratio format, the NDVI in the sunlit areas is negligibly influenced by the topography. Therefore, for sunlit areas in the image, compensating the reflectance for direct light is not required, i.e., α = 0. For shadow areas in the image, it is necessary to compensate the reflectance for direct light, i.e., α > 0. Combining (13) and (14), we can obtain where ρ dc is the surface reflectance compensation for direct light and ρ orig is the pixel reflectance before correction. The difference in the reflectance of land cover in the sunlit and shadow areas causes different NDVI values in the sunlit and shadow areas for the same type of land cover. The ρ NIR and ρ R used to calculate NDVI in the sunlit areas contain the surface reflectance for direct light and for diffuse light. However, the ρ NIR and ρ R used to calculate NDVI in the shadow areas only contain the surface reflectance for diffuse light. NTSEC aims to simulate the surface reflectance for direct light using the surface reflectance for diffuse light, and the simulated reflectance for direct solar light is used as a compensation value to be summed with the original reflectance for calculating the corrected NDVI.
Therefore, the NTSEC is calculated as follows:

B. Calculation of the Variable Coefficient (α)
In accordance with the principles of shadow composition, shadows can be categorized as self and cast shadows [38]. Self shadows are the those formed on the back side of mountains [39], [40], while cast shadows are formed on surfaces that do not receive direct solar illumination due to shading by other mountains [41], [42]. Unlike self shadows, cast shadows cannot be effectively identified by calculating the cosi. Thus, constructing a shadow index is a simple and effective way to extract both self shadows and cast shadows.
To analyze the differences in the spectral characteristics of different land covers in shadow and nonshadow areas, we selected six major land cover types within the two study areas and divided them into either shadow or nonshadow covered objects. Shadow covered objects include vegetation and bare land in shadow, while nonshadow covered objects include vegetation and bare land in nonshadow areas, impervious surfaces, and water bodies. A total of 200 samples were selected from the two study areas for each land cover type. Table II lists the number of samples for all six types. These samples were grouped according to whether they are shadow or nonshadow covered objects to plot the mean spectral signature curves (see Fig. 2). To ensure the purity of each sample class, we selected the central region of each object as the sample because they are unlikely to be confused with other objects [43]. Fig. 2 shows that the shadow covered objects have relatively higher reflectance in the Coastal band and the NIR band than in the other spectral bands. However, the reflectance of nonshadow covered objects is not the highest in these two bands. In addition, the Coastal band has a higher reflectance than the Green band for shadow covered objects, while the opposite is true for nonshadow covered objects. This finding is similar to that of Sun et al. [43]; however, their data were sourced from Sentinel-2A MSI images (with different spectral features than the Landsat 8 OLI). Sun et al. [43] demonstrated that the reflectance at B1 (the Coastal band) of shadow covered objects in Sentinel 2A MSI images is higher than that at B3 (the Green band) and that the reflectance at B8 (the NIR band) is lower than that at B9 (the Water vapour band); however, the opposite was observed for nonshadowed objects in their study. Thus, the subtraction of B3 and B8 from B1 and B9 will result in positive values for shadow covered objects and negative values for nonshadow covered objects, which form the expression of SEI. As B9 is not included in the Landsat 8 OLI data and considering that the wavelength of B9 is closest to that of the NIR band, we replaced B9 with the NIR band and modified the SEI to SI with the following expression: As described in Section III-A, α should be an SI that can distinguish shadows from nonshadows. The value of α for shadows is 0. For nonshadows, its value is between 0 and 1 and increases linearly with shadow intensity. Therefore, it is necessary to obtain the threshold value to distinguish between shadows and nonshadows using the thresholding method. Any areas greater than the threshold value are considered shadows, and areas less than the threshold value are considered nonshadows. In this article, the threshold value (c) of SI was calculated using Otsu's method [44]. Otsu's method is based on a histogram shape that assumes a distinctly bimodal image histogram and calculates the optimum threshold separating the two types, such that their intraclass variance is minimal and their interclass variance is maximal [45]. Since the variable coefficient α ranges from 0 to 1, it is necessary to stretch the SI values greater than or equal to the threshold value to 0-1. Meanwhile, an SI less than the threshold value is replaced by 0, as in where c is the threshold value for detecting shadow pixels using SI, and SI max is the maximum value of SI in the study area image.

A. Evaluation of the Shadow Detection Accuracy of SI
Numerous shadow pixels were marked as true shadow pixels in the image, and the recall ratio R r and precision ratio R p were calculated to reflect the detection accuracy of SI for shadows [2]. The calculation formula is shown as follows: where N tds is the number of true shadow pixels in the detected shadows, N ts is the number of true shadow pixels, and N ds is the number of detected shadow pixels.

B. Visual Inspection
Visual inspection is the easiest and most conventional method for inspecting the correction results. The main visual inspection criteria are whether: 1) the difference between the NDVI values of the shady slopes and the sunny slopes is reduced after correction; 2) the NDVI values in the shadow area are clearly overcorrected or undercorrected; 3) there are any null values in the corrected NDVI values; 4) the NDVI values in the sunlit areas are affected after correction; 5) the corrected NDVI images tend to be smoother compared with those before correction.

C. Correlation Analysis
The correlation between cosi and the band reflectance (or NDVI) can be used as a criterion to judge the elimination of the effect of topographic shadows [46]. Since cosi can indicate the degree of illumination, a successful NDVI topographic shadow effect correction model should eliminate or substantially reduce the correlation between NDVI and cosi.

D. Comparison of NDVI for Different Land Cover Types Between Sunny and Shady Slopes
In order to analyze the NDVI correction effect of the NTSEC method for different land cover types in shadow areas, the corrected NDVI values of the shady slopes (self shadows and cast shadows) and sunny slopes of eight different land cover types were compared. The relative error of mean NDVI (REMN) was calculated for each land cover type, and the correction effects of self shadow and cast shadow areas were quantified. The formula for this calculation is shown as follows: where NDVI sunlit is the NDVI values in adjacent sunny slopes without shadows and NDVI shadow is the NDVI values in self shadow or cast shadow areas. Compared with cast shadows, determining self shadows is simpler. Self shadows occur when the slope angle of the pixel in the solar incidence direction is greater than or equal to the solar elevation angle. The slope angle in the solar incidence direction (θ si ) can be calculated by using the following formula: When θ si ≥θ s , the pixel represents a self shadow. Once the self shadows are identified, identifying the cast shadows is easy, as any remaining shadows are cast shadows.

E. Comparison With Verification NDVI
The Gaoligong Mountains are located at low latitudes, and the land cover remains relatively constant throughout the year. In this article, the correction results of the Gaoligong study area were analyzed, with the April image used for verification. The solar altitude angle, at the time the image was acquired, was 62.60°, and there was less shadow covered area in the whole image at this time. Since NDVI can reduce some of the topographic shadow effect and April did not fully enter the rapid vegetation growth period, the NDVI changes caused by phenological changes were small [13]. Therefore, the NDVI values calculated from the April images were used as the verification NDVI. Note: δρ λ is the corrected reflectance for band λ; ρ λ is the original reflectance; θ z is the solar zenith angle (provided by the sensor metadata); i is the incidence angle with respect to surface normal; θ s is the terrain slope calculated from the DEM; k λ is called the Minnaert constant, which is the slope of the linear regression between x = log(cosicosθ s ) and y = log(ρ λ cosθ s ); b λ and m λ are the regression coefficients between cosi and ρ λ (ρ λ = b λ + m λ cosi); C λ is the quotient between b λ and m λ obtained by the previous equation in parentheses;ρ λ is the average reflectance of the image.

F. Outlier Analysis
Previously, most of the correction methods for the NDVI topographic shadow effect were based on topographic correction models, while the evaluation of TC algorithms was mostly carried out in favorable conditions. Overcorrections of self shadow areas and a few problematic pixels under low light conditions were usually excluded from the evaluation [47], [48]. Outliers can be generated when the correction model is unable to correct or overcorrect these pixels and are an important indicator for judging the adaptability of the model to extreme conditions. In this article, any pixels above the maximum and below the minimum original NDVI values were considered outliers, and their overall percentage in the image was calculated.

G. Sensitivity of NTSEC to NDVI of Different Land Cover Types
The NDVI of the sunlit areas both before and after correction was used as a reference because, according to the construction principle of the NTSEC method, NTSEC does not change the NDVI of the sunlit areas, meaning that the NDVI of sunlit areas is the same both before and after correction. For the same land cover type, if the difference between NDVI in the shadow areas and NDVI in the sunlit areas is significant after correction, it indicates that the method affects the inherent NDVI of the land cover, i.e., the method is insensitive to the inherent NDVI of the land cover. Conversely, if this difference is not significant, it means that the method is sensitive to the inherent NDVI of the land cover.

H. Comparison of NTSEC and Other TC Methods
Four topographic correction models (CC, MIN, VECA, and SE) were selected to correct the reflectance in the Red and NIR bands of the study area images and to calculate the NDVI. We compared NDVI corrected by the TC method with NDVI corrected by NTSEC. This process was performed to verify the effectiveness of the NTSEC model for NDVI correction in the shadow areas. Table III shows their expressions. The  four topographic correction models were selected because they perform excellently in most cases [46], [49], [50].

A. Result of Shadow Detection and the Variable Coefficient (α)
Taking the Landsat 8 OLI image of the Gaoligong Mountain study area as an example, shadow detection was performed using SI, where the threshold c of SI was obtained using Otsu's method (c = −0.01 in this article). The detected shadows are shown in Fig. 3(a). As can be seen from the figure, most of the shadows in this image have been detected, and the comparison with the shadow edges in the original image shows that the shadows detected by SI are consistent with the shadow areas of the image. To quantitatively evaluate the shadow detection accuracy of SI, 7989 shadow pixels and 10161 nonshadow pixels were marked as either "true shadow" or "true nonshadow" to calculate N tds , N ts , and N ds [see Fig. 3(b)]. For the shadow detection results [see Fig. 3(a)], N tds was 7425, N ts was 7989, and N ds was 7871, while R r and R p reached 92.94% and 94.33%, respectively. This demonstrates that the SI can correctly separate the shadow pixels and satisfy the requirement of the shadow index in the NTSEC model. The variable coefficient α in the NTSEC model can be calculated by substituting the SI and its threshold c into (18) [see Fig. 4(a)]. The self shadows and cast shadows, as determined by the process described in Section Ⅳ-D, are shown in Fig. 4(b).

B. Visual Inspection of the NDVI Correction Results
Fig . 5 shows the visual comparison results of NDVI both before and after correction using different methods. In order to ensure the consistency of the NDVI display, the corrected NDVI is 1 for pixels larger than 1 and −1 for pixels smaller than −1. The difference in NDVI between shady and sunny slopes before shadow effect correction is quite pronounced, particularly in rugged mountainous areas with steep slopes [see Fig. 5(b)]. The correction results for CC and VECA were mostly identical and all show overcorrection in a few pixels. The influence of shadows on NDVI was not eliminated in areas with significant terrain undulations [see The corrected NDVI using different methods for the Helan Mountain study area is shown in Fig. 6. The SE corrected results appeared to almost entirely eliminate the influence of topographic shadows on NDVI, particularly in nonvegetation areas [see Fig. 6(e)]. However, the NDVI correction for vegetation areas was obviously insufficient, and such results lead to distorted NDVI values of the corrected vegetation. In addition, similar to the Gaoligong study area, SE still had a few noise spots at the junction of the sunny slope and the shady slope. The voids in the correction results of MIN was greater than in the Gaoligong Mountain study area, and the CC and VECA did not compensate significantly for NDVI in the shadow areas compared with before the correction [see Fig. 6(b) and (d)]. The reason for this is that the terrain in the study area undulates significantly, and most of the shadows are formed by mountain shading (i.e., cast shadows). The TC methods are not sensitive to such shadows and, thus, they cannot be effectively corrected [2], [32]. Therefore, the correction results of three TC methods, CC, MIN, and VECA, had large differences between the NDVI of the sunny slopes and shady slopes, which means that the entire NDVI image was not smooth. In contrast, the correction results of NTSEC were smoother than those of CC, MIN, and VECA, and the correction effect on vegetation was better than that of SE [see Fig. 6 To further analyze the difference in NDVI between the corrected shady slopes (self shadows and cast shadows) and the sunny slopes, a subregion (see red rectangle, Fig. 5) was selected for detailed visual inspection, as shown in Fig. 7. When combining the analysis of the self shadow and cast shadow areas [see    We calculated the mean, standard deviation (Std), and coefficient of variation (CV) of NDVI both before and after correction for the sunny and shady slopes (self shadows and cast shadows) in the subregion (see Table IV). Among the correction results of the five methods, the corrected NDVI value exceeded the NDVI limit in all three methods except for the NTSEC and MIN methods. Therefore, we excluded the fraction exceeding the NDVI range values (i.e., the outliers) from the statistical calculations to ensure the accuracy of the statistical results. A closer mean NDVI of the corrected sunny slope to the mean of the original NDVI, a closer mean NDVI of the corrected self shadows and cast shadows area to the mean NDVI of the sunny slope, and smaller Std and CV values, are considered better [32], [46].
For the sunny slope, the mean, Std, and CV of the NTSEC correction results were the same as the original results, indicating that this method does not correct the NDVI of the sunny slope and retains the true NDVI information. For the self shadow areas, the means were closest to the sunny slope for CC and VECA, but they had high Std and CV values. Moreover, considering that we excluded the outliers from the calculations, this indicates that the NDVI correction of CC and VECA for the self shadow areas was highly unstable. This may be caused by overcorrection by these two methods. For the cast shadow areas, it is clear that the NDVI of the cast shadow areas was not as well corrected as the NDVI of the self shadow areas by all the methods (the mean NDVI of the cast shadow areas was smaller than that of the self shadow areas), except NTSEC. The NDVI of the cast shadow areas had higher stability after NTSEC correction.

C. Comparison With Verification NDVI
The density scatter plots between the original/corrected NDVI and the verification NDVI are shown in Fig. 8, and the coefficient of determination (R 2 ) and the root-mean-square error are shown in Table V. Since the spectral features of the water bodies in the images of the two periods differed significantly, we excluded the water bodies before analysis. Clearly, the scatter distribution of the original NDVI and the verification NDVI are concentrated and evenly distributed on both sides of the 1:1 line, with some symmetry [see Fig. 8(a)]; however, the original NDVI appears to be reduced at verification NDVI = 0.8, meaning the edge of the density distribution spills over. This indicates that the shadow area is undercorrected for NDVI. Although the density scatter plot results of MIN and NTSEC were similar to the original  NDVI results, the distribution of NTSEC had high symmetry and the highest R 2 among all the methods (see Table V), suggesting that NTSEC has a good correction effect on NDVI in shadow areas. The density scatter plots of CC, VECA, and SE were not concentrated and had a band distribution with no symmetry [see Fig. 8(b)-(d)]. In addition, the corrected NDVI at verification NDVI = 0.8 showed significant spillover, indicating that the corrected NDVI in the shadow areas was both undercorrected and overcorrected, which is why the R 2 of these three methods is low.
The NDVI profiles for the dashed lines in Fig. 7 are provided in Fig. 9. The dashed lines are the verification NDVI, and the solid lines are the NDVI after correction by the various methods. Fig. 9 shows that MIN [see Fig. 9(c)] was unable to correct the NDVI in the self shadow areas, leading to the creation of the vacant part of the figure, and that the CC [see Fig. 9(b)] and VECA [see Fig. 9(d)] methods overcorrect the NDVI in the self shadow areas. This result corroborates the findings in the visual inspection. Since the land cover on this NDVI profile was all vegetation, and almost all forest, the NDVI fluctuations should be small, and the corrected NDVI fold should be similar to the verification NDVI fold [32]. Therefore, we believe that NTSEC has better NDVI topographic shadow correction performance than the TC methods.

D. Analysis of Different Land Cover Types 1) Relative Error of Self and Cast Shadows:
The REMN values of the different land cover types are shown in Fig. 10. As CL and WB were almost all located in flat sunlit areas and were not affected by topographic shadow effect, we analyzed six types of features (i.e., all except for CL and WB). As shown in Fig. 10, the REMNs of almost all land cover types after MIN and SE correction were higher than before correction. Combined with the visual inspection results shown in Fig. 5, it can be observed that these two methods produce a clear undercorrection [see Fig. 5(d) and (f)]. For the self shadow areas, CC, VECA, and NTSEC all had lower REMNs, indicating that all three methods have good NDVI correction ability in self shadow areas. For the cast shadow areas, the correction result of NTSEC had the lowest REMN for the remaining five land cover types except for CUL, indicating that this method has a better correction ability of NDVI in the cast shadow areas than the TC method. Overall, the maximum and minimum REMN of NTSEC were 5.13% and 0.04%, respectively, with a lower range than the four TC methods.
2) Sensitivity of NTSEC to NDVI Values of Different Land Cover Types: In order to analyze whether there is a high sensitivity to different land cover types after NTSEC correction, we plotted the distribution of NDVI values both before and after correction for each land cover type (see Fig. 11). As shown in Fig. 11, the NDVI points before and after correction in the sunlit areas all overlapped, which indicated that NTSEC does not change the NDVI of the various land cover types in the sunlit areas. The corrected NDVI in shadow areas was slightly lower than the NDVI in sunlit areas for FL and GL, and the corrected NDVI in the shadow areas was slightly higher than the NDVI in the light areas for CUL. However, the difference between their NDVI in shadow areas after correction and the NDVI in sunlit areas was much smaller than their difference before correction, and there was no movement away from the NDVI point in sunlit areas. In addition, NTSEC does not apply any correction to COL and WB, as they are not influenced by topographic shadows. This indicates that NTSEC does not have a significant effect on the NDVI properties of different land cover types, meaning it has a high sensitivity to the NDVI inherent to This indicates that NTSEC does not have a significant effect on the NDVI properties of different land cover types, meaning it   has a high sensitivity to the NDVI inherent to different types of land cover.

E. Outlier Analysis
The outliers statistics of the NDVI before and after correction are shown in Fig. 12. Like the visual inspection results, the MIN corrected outliers showed null values due to the defect of the model expression and the percentage of outliers was close to 5%. For such results, we believe that the model does not apply to the correction of NDVI topographic shadow effect in areas with large terrain relief. The outliers for both CC and VECA were less than 0.5% and the difference in the number of positive and negative outliers was not significant. However,  there were many outliers in the correction results of SE, and these outliers corresponded exactly with the noise points found in the visual inspection. It is clear that there are no outliers after NTSEC correction, which means that the introduction of errors caused by topographic shadow correction is avoided using this method. In fact, for topographic correction models, the outliers increased as the solar altitude angle decreased and the slope angle increased [46]. Although these outliers may only become more pronounced under extreme conditions, a good NDVI topographic shadow correction model needs to be able to maintain a good correction even under extreme conditions.

F. Correlation Analysis
To observe the variation of NDVI with cosi for the sunny slopes, self shadow, and cast shadow areas, the scatter plots of cosi and NDVI were plotted for 1800 pixels selected in the Helan Mountain study area, as shown in Fig. 13. Table VI shows the correlation between NDVI and cosi.
As shown in Table VI, the best decorrelation ability was shown by CC and VECA, with a correlation coefficient of only 0.03. SE and NTSEC had similar correlation coefficients of less than 0.08, while MIN had the largest correlation coefficient, which exceeded even the original NDVI. However, Fig. 13 shows the same results as those of the outlier analysis (see Fig. 12), with outliers appearing for CC, MIN, VECA, and SE (i.e., the outlier distribution phenomenon). The outliers of CC, MIN, and VECA appeared in the self shadow areas, resulting in NDVI overcorrecting in these areas. The three methods also showed undercorrection in the cast shadow areas, which corresponds to the results of the visual inspection. The correction ability of NTSEC for self and cast shadows was almost identical. The corrected NDVI of both these types of shadow was significantly higher than before correction, with an even distribution of the NDVI points [see Fig. 13(f)]. In addition, the method had a good decorrelation ability (see Table VI) and the corrected NDVI had very little correlation with solar illumination. Therefore, the NTSEC method can effectively correct the NDVI topographic shadow effect in complex terrain.

VI. DISCUSSION
NDVI can quantify the growth of vegetation and is an important source of inversion data for vegetation biophysical and biochemical parameters such as vegetation cover, leaf area index, biomass, and photosynthetically active radiation absorption [51], [52]. Shadows are the major obstacle for accurate NDVI calculations in rugged terrain, and their presence greatly reduces the accuracy of applying NDVI to estimate vegetation parameters [32]. Therefore, the effect of eliminating shadows on NDVI is important for quantitative research. The analysis in this article clearly shows that the proposed NTSEC method has good potential for correcting NDVI in shadow areas for complex terrain.
The cosine of the solar incidence angle (cosi) is used as an important variable for the correction of the topographic shadow effect in most TC methods. However, for mountain surfaces with steep slopes, self shadows, and cast shadows from the surrounding terrain are common, especially in images with low solar altitude angles. Under such conditions, cosi cannot truly represent the actual illumination situation [53]. Specifically, the cosi of the self shadow areas is negative, and the cosi of the cast shadow areas is the same as the sunlit areas. Due to these issues, there is a significant difference in the correction ability of "cosi-dependent" topographic correction methods for self shadows and cast shadows, which are easily overcorrected in the self shadow areas and undercorrected in the cast shadow areas. Therefore, it is necessary to consider the influence of topographic correction methods on self shadows and cast shadows separately when investigating the correction of the topographic shadow effect in rugged mountainous areas.
From the results of NDVI correction in the shadow areas of the Landsat 8 OLI images of the two study areas, we found the appearance of several voids in the MIN [see Figs. 5(d) and 6(c)]. This is caused by the defective model expression of the MIN, which is not corrected if the solar incidence angle is greater than 85° [54]. In other words, the pixels with a cosi of less than 0 get invalid NDVI values after MIN correction [see Fig. 13(c)]. However, for complex mountainous areas with large terrain undulations, pixels with a cosi of less than 0 are not in the minority. The influence of a cosi less than 0 on the TC model does not only exist in MIN, but also exists in CC and VECA and is reflected as overcorrection [46]. Even though CC introduces a constant C to control the overcorrection, in mountainous areas with large terrain relief, more pixels have a cosi+C close to 0 and these pixels are overcorrected. Combining Figs. 7 and 9 shows that CC and VECA overcorrected for self shadows and undercorrected for cast shadows, while MIN had no self shadow correction capability. In fact, SE corrected reflectance better than most TC methods in complex mountainous areas with large terrain undulations [46], [55]. However, since its correction process reduces the reflectance of the sunlit areas while compensating for the reflectance of the shadow areas, the application of this method to the correction of the NDVI topographic shadow effect leads to a reduction in the corrected NDVI [see Figs. 6(e) and 13(e)], which is not helpful for quantitative analysis. In contrast, the NTSEC method introduces an index SI reflecting the shadow intensity to construct the constraint variable α to control the reflectance compensation for direct light. The range of α is 0-1; α is 0 for sunlit areas in which there is no need to compensate the reflectance for direct light, while for shadow areas, α gets closer to 1 as shadows are enhanced. This solves the problem of the lack of a valid variable to represent the illumination of the shadow areas, since α is constrained for full shadows and does not distinguish between self and cast shadows. It only represents the shadow intensity. Therefore, the NTSEC method does not overcorrect the NDVI in the self shadow areas and also improves the undercorrection in the cast shadow areas.
According to density scatter plots between the original/NTSEC corrected NDVI and the verification NDVI [see Fig. 8(a) and (f)], we can see that there is a high correlation between the original NDVI and the verification NDVI. However, the NDVI dynamic range "sags" in areas with high verification NDVI values, indicating that the NDVI in the shadow areas is low. NTSEC correction fixes the phenomenon and the NDVI dynamic range does not shrink, indicating that NTSEC does not significantly alter the inherent NDVI information of the land cover, a conclusion that is also confirmed by Fig. 11. In contrast, after the correction of CC, VECA, and SE, the vertical banding phenomenon and the NDVI dynamic range expanded significantly [see Fig. 8 , which led to a significant decrease in the correlation between the corrected NDVI and the verification NDVI (see Table V). This indicates that this method would change the inherent NDVI information of the land cover. The reason for this change in the intrinsic NDVI of the land cover caused by these three TC methods is the presence of outliers (see Fig. 12). For evaluating the NDVI correction capability of the different methods, we have limited those positive and negative outliers that are outside the theoretical range of NDVI values to be equal to 1 or −1. The purpose of this was to keep the corrected NDVI within a reasonable range. However, this does not mean that outliers due to methodological flaws do not exist. Sola et al. [46] stated that a good correction method for the topographic shadow effect should eliminate outliers. The NTSEC method does not produce any outliers, which is important, and therefore, NTSEC is a stable and effective correction method for the NDVI topographic shadow effect.
The construction process of NTSEC constrains the reflectance compensation for direct light by α. The α is not physically meaningful, but considering that the construction process of the NTSEC method is based on radiation transfer theory, NTSEC should be a semiempirical correction model. The semiempirical properties make the correction performance of NTSEC for NDVI somewhat dependent on the threshold c of α. Therefore, it is necessary to analyze the sensitivity of the correction results to c. The corrected NDVI values were calculated for the subregion of the Gaoligong Mountain study area at a threshold c from −0.03-0.02 in steps of 0.005, and the corrected NDVI in sunny and shady slopes were counted (see Fig. 14). When the threshold c is less than −0.015, the NDVI of the sunny slope becomes larger, i.e., NTSEC affected the NDVI of the sunlit areas, which is contrary to the core NTSEC concept. In addition, the corrected NDVI difference between shady and sunny slopes was the smallest when c was in the range of −0.015 to −0.01, suggesting that the best correction result was achieved in this range. Analysis of Fig. 14 suggests that the optimal c is the critical value at which the NDVI of the sunny slopes changes after correction. In Section Ⅴ-A, the threshold value calculated using Otsu's method lies in the interval, but it needs to be slightly lowered for c value to achieve the best correction. Therefore, NTSEC has a high sensitivity to threshold c. We suggest that, when applying the NTSEC method, after threshold c is calculated, a small range of c around the threshold can be selected for correction tests, and the optimal threshold can be judged by checking whether the NDVI of the sunny slopes is affected.
Although NDVI has good performance in most cases, it has the obvious disadvantage that it tends to be saturated in areas with high vegetation cover. EVI overcomes this problem, but has the disadvantage that it is much more anisotropic than NDVI and is strongly influenced by light, which is why the topographic shadow effect of EVI is more difficult to correct [12], [56]. We tried to introduce the proposed NTSEC method to the correction of EVI, but the correction results were not satisfactory. The area receiving more sky diffuse radiation and terrain-reflected radiation was able to get a good correction, but the area receiving less sky diffuse radiation and terrain-reflected radiation was hardly corrected (see Fig. 15). In fact, this problem also occurred in the correction of NDVI, as demonstrated by the poor correction of NDVI in nonvegetation shadow areas with very weak illumination intensity [see Fig. 6]. Because of the poor illumination in the shadow areas, the original reflectance before correction stays at a very low level, and the inner reflectance differences among different land cover types in the shadow areas cannot be fully displayed in the original image. Consequently, the type differences cannot be well represented in the corrected image [53]. NDVI as a vegetation index in complete band ratio format could reduce this effect to a certain extent, but EVI does not have this attribute. Therefore, the method proposed in this article is not recommended for application to EVI. However, it is worth mentioning that RVI is also a vegetation index in complete band ratio format. After our tests, the NTSEC method was also found to be applicable for the correction of RVI, and the stability and validity of the correction results were almost the same as NDVI (see Fig. 15). This result was expected because the formula for NDVI is mathematically derivable from RVI [12]. Considering the fact that NDVI is more sensitive to low vegetation cover areas and RVI is more sensitive to high vegetation cover areas, we suggest that the NDVI or RVI can be chosen according to the vegetation cover of the application area when applying the NTSEC method. Moreover, it is undeniable that EVI has other advantages, such as overcoming the supersaturation problem of NDVI, reducing atmospheric influence, and correcting the scattering of atmospheric aerosols and soil background. It is also imperative to study a method to correct for the topographic shadow effect of EVI. Therefore, in future research, we will focus on the development of new shadow indexes that are more sensitive to the illumination to construct α. The new α should have higher sensitivity, thereby improving the correction capability of NTSEC for EVI and NDVI in extreme conditions (nonvegetation areas with very weak light intensity).
In this article, we selected a DEM with 30 m resolution to ensure the consistency of the resolution with Landsat 8 OLI data. In this case, resampling during the data process could thus be simply avoided in order to retain the entire original spectral information [29], [57]. However, there is no consensus yet regarding the best combination of a DEM and original image resolutions for topographic correction. The advantage of higher resolution DEM data is that they can reduce errors in slope calculation and, thus, improve the accuracy of V d and C t calculations. In addition, we used a simple approximation method to obtain V d and C t . The advantage of this method is that the calculation process is simple and easy to implement as it does not consider the obstructing effect of the surrounding topography. However, there are also disadvantages to this method. For long slant surfaces or semirugged areas, this method does not cause large errors, but for rugged terrain, the errors caused by this method must be considered [58]. Dozier and Frew [37] and Li et al. [42] introduced horizon angles and proposed a more accurate method to calculate V d and C t , but its calculation process is complicated and difficult to apply in practice. This method was simplified by Mousivand et al. [58], who demonstrated that calculating V d in 64 directions for each pixel can achieve an accurate estimation, providing an important reference for improving the accuracy of the NTSEC model. Therefore, in future work, we will use a higher resolution DEM (e.g., ALOS 12.5 m DEM data) for further testing of the NTSEC model. Furthermore, we will obtain V d and C t using more advanced methods to improve the NTSEC model.

VI. CONCLUSION
To address the problem of NDVI outliers and distortions in shadow areas in rugged terrain, this article proposes an NTSEC method that can be applied to areas with large terrain relief. Landsat 8 OLI images of two regions were used to test the NTSEC and were evaluated using multiple strategies. The correction results were compared with four commonly used TC models (CC, MIN, VECA, and SE). In areas with large terrain relief, the CC, VECA, and SE methods were prone to overcorrection in self shadow areas and caused NDVI outliers, while MIN was unable to correct NDVI in self shadow areas. For cast shadow areas, all four TC methods suffered from undercorrection and NDVI distortion. NTSEC was able to suppress the overcorrection of self shadow areas, and the corrected NDVI results did not have any outliers. At the same time, NTSEC rectified the undercorrection of the cast shadow areas, so that the NDVI correction results of the cast shadow areas were almost the same as those of the self shadow areas. The NDVI differences between shady and sunny slopes were significantly reduced after correction. In addition, the NDVI of different land cover types still had interclass stability after NTSEC correction. For mountainous areas with large terrain relief, NTSEC provides a simple and effective method for correcting NDVI distortion caused by topography, which is difficult to achieve using traditional TC methods such as CC, MIN, VECA, and SE.
Wenbin Xie received the M.S. degree in surveying and mapping engineering from the Kunming University of Science and Technology, Kunming, China, in 2020.
He is currently with the Yunnan Institute of Geological Surveying and Mapping. His research interests include spatio-temporal big data, geographic information systems, and data processing.