A Kernel-Driven BRDF Approach to Correct Airborne Hyperspectral Imagery over Forested Areas with Rugged Topography

: Airborne hyper-spectral imaging has been proven to be an e ﬃ cient means to provide new insights for the retrieval of biophysical variables. However, quantitative estimates of unbiased information derived from airborne hyperspectral measurements primarily require a correction of the anisotropic scattering properties of the land surface depicted by the bidirectional reﬂectance distribution function (BRDF). Hitherto, angular BRDF correction methods rarely combined viewing-illumination geometry and topographic information to achieve a comprehensive understanding and quantiﬁcation of the BRDF e ﬀ ects. This is in particular the case for forested areas, frequently underlaid by rugged topography. This paper describes a method to correct the BRDF e ﬀ ects of airborne hyperspectral imagery over forested areas overlying rugged topography, referred in the reminder of the paper as rugged topography-BRDF (RT-BRDF) correction. The local viewing and illumination geometry are calculated for each pixel based on the characteristics of the airborne scanner and the local topography, and these two variables are used to adapt the Ross-Thick-Maignan and Li-Transit-Reciprocal kernels in the case of rugged topography. The new BRDF model is ﬁtted to the anisotropy of multi-line airborne hyperspectral data. The number of pixels is set at 35,000 in this study, based on a stratiﬁed random sampling method to ensure a comprehensive coverage of the viewing and illumination angles and to minimize the ﬁtting error of the BRDF model for all bands. Based on multi-line airborne hyperspectral data acquired with the Chinese Academy of Forestry’s LiDAR, CCD, and Hyperspectral system (CAF-LiCHy) in the Pu’er region (China), the results applying the RT-BRDF correction are compared with results from current empirical (C, and sun-canopy-sensor (SCS) adds C (SCS + C)) and semi-physical (SCS) topographic correction methods. Both quantitative assessment and visual inspection indicate that RT-BRDF, C, and SCS + C correction methods all reduce the topographic e ﬀ ects. However, the RT-BRDF method appears more e ﬃ cient in reducing the variability in reﬂectance of overlapping areas in multiple ﬂight-lines, with the advantage of reducing the BRDF e ﬀ ects caused by the combination of wide ﬁeld of view (FOV) airborne scanner, rugged topography, and varying solar illumination angle over long ﬂight time. Speciﬁcally, the average decrease in coe ﬃ cient of variation (CV) is 3% and 3.5% for coniferous forest and broadleaved forest, respectively. This improvement is particularly marked in the near infrared (NIR) region (i.e., > 750 nm). This ﬁnding opens new possible applications of airborne hyperspectral surveys over large areas. To


Introduction
Airborne hyperspectral remote sensing technology improved significantly during the past decades, permitting a wide range and increasing number of applications such as land-cover type classification, forest biodiversity assessment, and quantitative estimation of land surface parameters [1][2][3][4][5]. Airborne scanners equipped with a wide field of view (FOV) provide a sampling of the bidirectional reflectance distribution function (BRDF) for very high spatial resolution (VHR) imagery [6,7]. Multiple angular scans associated with various solar irradiance for the same survey, together with the rugged topography of the area, results in strong anisotropic scattering effects for the studied area even for the same land-cover type, which is represented by the BRDF [8][9][10][11][12][13][14]. BRDF features can either be treated as perturbing effects or as additional information, this depending on the targeted objectives and applications [15][16][17][18]. BRDF correction is advised when mosaicking multiple hyperspectral scenes of a region, resulting in a spectrally homogeneous analysis-ready mosaic [19,20]. It is worth mentioning that an angular normalization procedure is mandatory to reduce the confounding non-biophysical signals and to sustain the interpretation of the airborne imagery in respect to land surface properties. VHR airborne hyperspectral images acquired over a rugged topography requires disentanglement of soil, vegetation canopy, and relief BRDFs, which remains challenging [21][22][23]. The problem is even more complex with airborne data due to the strong influence of the local slope and aspect.
A correct understanding of the BRDF effects is crucial for reflectance-based quantitative remote sensing inversion models [24,25]. Previous studies developed various statistical correction methods to remove the BRDF effects from satellite and airborne images (e.g., [20,[26][27][28][29]). However, since these methods are not based on physical mechanisms of reflected light, they cannot be universally applied. Physical-based BRDF corrections usually require more detailed surface characteristics like the canopy structure and the ratio of direct and diffuse radiation for each pixel in order to simulate the radiative transfer mechanisms. This can be a time consuming procedure, which clearly limits its large scale use [30][31][32][33], while it remains still a challenging approach for even VHR airborne hyperspectral imagery. In this paper, the BRDF-RT model presented is based on a semi-empirical kernel-driven BRDF approach (e.g., [34]) utilizing the Ross-Li-Reciprocal model employed for Moderate Resolution Imaging Spectroradiometer (MODIS) BRDF/Albedo product suite (e.g., [35][36][37][38]) to normalize the off-nadir viewing directional reflectance to nadir BRDF adjusted reflectance. The MODIS BRDF algorithm and its physical basis have demonstrated to be a robust method and it has been widely used for BRDF correction of satellite images, for example to map land-cover types at the global scale [39][40][41].
The kernel-driven BRDF approach offers a synthetic vision of the complex scattering mechanisms of land surfaces via the linear combination of an isotropic kernel, a geometric-optical kernel, and a volumetric kernel. Various geometric-optical kernels (e.g., Li-Sparse-Reciprocal, Li-Dense-Reciprocal, Li-Transit-Reciprocal) and volumetric kernels (e.g., Ross-Thick, Ross-Thin, Ross-Thick- Maignan) having different mathematical expressions were proposed to optimally correct the observed scene. The prevalent combination of kernels in the literature is Ross-Thick and Li-Sparse-Reciprocal which describes the scene with a thick vegetation layer for the canopy and a sparse horizontal distribution of the stems. Numerous studies have attempted to extend the applicability of the MODIS satellite-specific BRDF algorithm to airborne hyperspectral images (e.g., [7,31,42]). Found issues were that different spatial scales of the image result in different BRDF effects, which limits the general applicability of the MODIS BRDF algorithm [43,44]. Hyperspectral images acquired using wide FOV airborne scanners typically have high spatial and spectral resolution. At the opposite to moderate resolution satellite images, pixels of a VHR image composed of less heterogeneous thematically different land-cover types is not common [45,46] but more physically as the limited number of vegetation attributes [15] infers precipitation per year. The vegetation in the region is dominated by Pinus kesiya var. langbianensis and Monsoon evergreen broadleaved tree species (e.g., Castanopsis echidnocarpa, Castanopsis hystrix, Schima wallichii, Betula alnoides Buch.-Ham. ex D. Don) [65]. The area surveyed for this study is characterized by mountainous topography, with altitude varying between 500 and 2200 m above sea level and an average slope of 27 • .

Remote Sensing Data
The hyperspectral data were acquired, respectively, under cloud-free conditions on 04 April and 06 April 2014 using the CAF-LiCHy system [59] at a nominal flight altitude approximately of 1500 m above mean ground elevation. Figure 1  Then we applied the approach to a larger area (04 April 2014 flight) to assess the repeatability of the RT-BRDF algorithm, with results shown in Section 4.6.
The CAF-LiCHy system comprises an Airborne Imaging Spectrometer for Applications (AISA) Eagle II hyperspectral sensor (SPECIM Spectral Imaging Ltd., Oulu, Finland), an IGI DigiCAM-60 digital camera (Integrated Geospatial Innovations, Kreuztal, Germany), and a RIEGL LMS-Q680i Light Detection And Ranging (LiDAR) scanner (Riegl GmbH, Horn, Austria), providing hyperspectral images simultaneously with high resolution CCD images and LiDAR point cloud data [59]. The CCD images are orthorectified with 0.18 m spatial resolution.
The AISA Eagle II is a push broom imaging system with FOV of 37.7 • . The hyperspectral images were obtained in 125/64 bands covering a spectral range spanning from 400 to 990 nm, with a spectral resolution of 4.6/9.2 nm. First, radiometric and geometric corrections of raw AISA Eagle II hyperspectral images were performed using CaligeoPRO v2.2.4 (SPECIM Limited, Oulu, Finland) which result in the radiance images with a ground sampling distance (GSD) of 1 m. Next, the hyperspectral images were atmospherically corrected using the ATCOR4 (ATmospheric CORrection 4) software (ReSe Applications LLC, Wil, Switzerland), to transform radiance into Bottom-Of-Atmosphere (BOA) apparent directional reflectance [7,66]. ATCOR4 calculates the database of atmospheric look-up table based on the MODTRAN5 (MODerate resolution atmospheric TRANsmission 5) code [67].
The LiDAR point cloud obtained has an average footprint density of 3 pts/m 2 . The LiDAR point cloud was separated into ground and non-ground return points data using the TerraScan software (Terrasolid, Helsinki, Finland) in order to generate a DEM and a digital surface model (DSM) with the resolution of 1 m. The DEM-derived slope and aspect were both calculated using a window size of 5 × 5 pixels in ENVI v5.3 (Exelis Visual Information Solutions, Boulder, CO, USA). The same GNSS and Inertial Navigation System (INS) module of the CAF-LiCHy leads to a high degree of geometric accuracy between the hyperspectral imagery and LiDAR-derived products.
Remote Sens. 2020, 12, 432 5 of 33   To assess the performance of the RT-BRDF model for the VHR airborne hyperspectral images (AISA Eagle II hyperspectral data in this paper), the MODIS BRDF/Albedo Model Parameters and Nadir BRDF-Adjusted reflectance product (MCD43A1 and MCD43A4 Collection 6), corresponding to the flight date and location, were chosen for further comparison.

Field Data
A field survey was conducted from November 2013 to January 2014 in the study area to collect 61 circular sample plots with a radius of 15 m. A total of 30 plots were located in Pinus kesiya var. langbianensis dominated stands, and the remaining 31 plots in broadleaved forest stands. Within each plot, tree species, tree height, under-branch height, diameter at breast high (DBH), and crown diameter were measured for all the individual trees with DBH greater than 5 cm. Plot center location was recorded using a Trimble R4 Global Navigation Satellite Systems (GNSS) receiver (Trimble Navigation Limited, Sunnyvale, USA), with final post-differential processing position error resulting within 1 m.

Methodology
The procedure comprises the following three main steps: (i) calculate the rugged topography kernel functions based on local viewing and illumination geometry; (ii) calculate the kernel coefficients of the RT-BRDF model for each spectral band and each forest type; and (iii) apply the RT-BRDF model to correct for the anisotropy effects of airborne hyperspectral data. The detailed processing chain is described in Figure 2.

Field Data
A field survey was conducted from November 2013 to January 2014 in the study area to collect 61 circular sample plots with a radius of 15 m. A total of 30 plots were located in Pinus kesiya var. langbianensis dominated stands, and the remaining 31 plots in broadleaved forest stands. Within each plot, tree species, tree height, under-branch height, diameter at breast high (DBH), and crown diameter were measured for all the individual trees with DBH greater than 5 cm. Plot center location was recorded using a Trimble R4 Global Navigation Satellite Systems (GNSS) receiver (Trimble Navigation Limited, Sunnyvale, USA), with final post-differential processing position error resulting within 1 m.

Methodology
The procedure comprises the following three main steps: (i) calculate the rugged topography kernel functions based on local viewing and illumination geometry; (ii) calculate the kernel coefficients of the RT-BRDF model for each spectral band and each forest type; and (iii) apply the RT-BRDF model to correct for the anisotropy effects of airborne hyperspectral data. The detailed processing chain is described in Figure 2.

Land-Cover Stratification (Pre-Classification)
Different land-cover types present a priori different anisotropic properties. Therefore, a priori knowledge of the type is needed first to extract the BRDF characteristics pre land-cover type. Moreover, different canopy structure can impact the anisotropic properties even for forested areas. Studies in the literature (e.g., [68][69][70]) indicate that the texture information of VHR airborne optical images can represent a proxy to describe forest structure variables (e.g., crown diameter, tree height, tree density or spacing). For these reasons we applied a Support Vector Machine (SVM) classification method [71] to our hyperspectral images, using as input the BRDF cover index [7], the three most prominent components of the hyperspectral reflectance images obtained with the Principal Components Analysis method [72], and the Gray Level Co-occurrence Matrix (GLCM) texture features [73]. The GLCM texture features (i.e., mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation; [74]) are extracted based on the first principal component data with a window size of 9 × 9 pixels. In this study, coniferous forest, broadleaved forest, bare soil, urban, and unclassified were classified, with most unclassified pixels in our dataset represented by the cast-shaded areas in the deep gullies and canopy gaps [71].

Local Viewing and Illumination Geometry
In aerial hyperspectral imaging, each pixel has a different viewing and illumination geometry due to different scanning angles, imaging time, and topographic features. The relative sun-surface-sensor geometry is fundamental for the RT-BRDF method. In order to obtain this information, we developed a two-step process to calculate the local viewing and illumination geometry for each pixel of the airborne hyperspectral images, based on the characteristics of the wide FOV scanner utilized and the topographic features of the area surveyed, i.e., (i) calculate viewing and illumination geometry based on flat topography assumption; and (ii) transfer the viewing and illumination geometry to local viewing and illumination geometry based on a realistic topography issued from a high resolution DEM.
In the first step, the viewing and illumination geometries are calculated under the assumption of a flat topography. Since each flight-line is acquired in a relatively short period of time (i.e., approximately 10 min), the solar position within this time period is assumed to be constant here. Thus, both solar zenith and azimuth angles are calculated based on date, average latitude and longitude, and median acquisition time in each flight-line. For each pixel, the image space coordinates are converted into object space coordinates using the co-linearity equations in Equations (1) and (2) [75]: where [u v w] T and [x y] T represent the pixel vector in the object space coordinates and image space coordinates, respectively; R is the rotation matrix calculated from the exterior orientation element ϕ, ω, and k ; f is the focal length of the optical camera; and [x 0 y 0 ] T is the principal point in image space coordinates. View zenith and azimuth angles are calculated for each pixel using Equations (3) and (4): Remote Sens. 2020, 12, 432 where θ v and ϕ v are the view zenith and azimuth angle, respectively.
Due to the fact that the majority of the species in our study area grows vertically independently of the slope of the terrain, and in order to apply the Li-Strahler geometric-optical model to sloping terrain, the ellipsoidal tree crown is replaced by spheres casting the same shadow area. The slope α is converted to the forest-adjusted slope α (Equation (5)). In the second step, the global coordinates obtained as described in Equation (1) are transformed into local coordinates based on converted slope and aspect angle of each pixel (Equation (6); [76]): where r/b represents the ratio of horizontal crown radius to vertical crown radius; [x y z ] T is the pixel vector in local coordinates; and R topo is the rotation matrix based on topographic features calculated from the converted slope α and aspect β angle of each pixel. Then, real local view angles are calculated for each pixel using Equations (7) and (8). Similarly, real local solar zenith angle θ s and solar azimuth angle ϕ s can be calculated according to each pixel using Equation (6): where θ v and ϕ v are the local view zenith and azimuth angle for each pixel, respectively. Finally, the cosine of the local solar zenith angle can be referred as illumination conditions, and it is calculated as in Equation (9): cos(θ s ) = cos(θ s )cos (α ) + sin(θ s )sin (α )cos(∆∅) (9) where θ s and ∆∅ are the solar zenith angle and relative azimuth angle between the aspect and sun azimuth based on flat topography, respectively.

Rugged Topography BRDF Model and Kernel Selection
We assume that each land-cover type in the images presents a homogeneous structure and does not vary with the slope and aspect of the underlying topography (i.e., relatively uniform anisotropic properties). The local viewing and illumination geometry result in a multi-angular data set for each land-cover type [77]. A kernel-driven BRDF approach (e.g., [34]), used in literature from the canopy to the satellite scale (e.g., [24,78,79]), was adopted in this study. Specifically, this model is based on the linear combination of volumetric, geometric, and isotropic scattering effects. However, despite having similar fitting capability to the Li-Sparse kernel, the Li-Transit kernel has demonstrated to perform better in the case of high solar and/or view zenith angles [80]. Therefore, in this study the Li-Transit-Reciprocal (LTR) kernel [81,82] was used as geometric scattering kernel, and the hotspot-revised Ross-Thick-Maignan (RTM) kernel [83] as volumetric scattering kernel. Using the local viewing and illumination angles in the RTM and LTR kernel functions, the RT-BRDF model can be expressed as in Equation (10): where: where θ v , θ s , ∆∅ , c, and λ are the local view zenith angle, local solar zenith angle, and local relative azimuth angle between the sun and the observer, land-cover type, and wavelength respectively; f iso , f vol , and f geo are coefficients of isotropic scattering, volume scattering and geometric scattering, respectively; k vol is the volume scattering kernel; and k geo is the geometric scattering kernel. A constant value of ξ 0 = 1.5 • is used to indicate the half-width of the hotspot for most of the targets [83]. The Li-Transit kernel requires to set the parameters for the tree crown geometry. In accordance with the field measurements, in this study we used h/b = 2 and b/r = 1 for coniferous forest; and h/b = 1.5 and b/r = 1 for broadleaved forest.

Determination of BRDF Model Coefficients
VHR airborne hyperspectral images generally contain a large number of pixels. For the purpose of an efficient processing and in order to obtain a comprehensive coverage of the viewing and illumination angles, a stratified random sampling method [84,85] was used to select a subset of pixels for the BRDF model fitting. For instance, under the assumption that the data set used in this study included a high number of pixels (N) in coniferous forest covered areas, a subset of n pixels were chosen from the N pixels (with n << N, with N ≈ 18,500,000 in this case) to build the BRDF model for coniferous forest in our study according to the following sampling rules ( Figure 3):

1.
n takes any values between 100 and 50,000, with step length of 100 pixels (i.e., 100 the root mean square errors (RMSE) between the observed band reflectance and the model fitted band reflectance is used to determine the optimal n (i.e., we define the optimal n here as the minimum n pixels with a relatively stable and minimum RMSE for the BRDF model fitting).
Remote Sens. 2020, 12, 432 10 of 33 Once the optimal number of pixels n is determined for each class type, the kernel coefficients for each band can be determined using the least square method, which solution is shown in Equation (23): in Equation (23), the , ,…, and , ,…, represent the volumetric and geometric scattering kernel for the 1 st , 2 nd ,…,n th pixel, respectively; and , ,.., is the observed band reflectance for the 1 st , 2 nd ,…, n th pixel. The Equation (23) can be written as a matrix form in Equation (24) and can be further expressed as a simplified form in Equation (25). Then the kernel coefficients can be calculated as = = .

BRDF Effects Correction for Image
The anisotropy factor (ANIF) is a ratio describing the relationship between directional reflectance and reference reflectance, and it is introduced to the BRDF correction model of the hyperspectral image. We followed the approach described in [31], and defined the ANIF with respect to the nadir reflectance as in Equation (26); therefore, the RT-BRDF corrected reflectance ( ) can be calculated using Equation (27): In order to obtain the optimal number of pixels for each spectral band (i.e., 125 bands, 400-990 nm), we developed an independent band normalized RMSE (Equation (22)) for each band: where RMSE min and RMSE max indicate the minimum and maximum RMSE for the different set n pixels (i.e., from 100 to 50,000 in this study), respectively. Once the optimal number of pixels n is determined for each class type, the kernel coefficients for each band can be determined using the least square method, which solution is shown in Equation (23): in Equation (23), the k 1 vol , k 2 vol , . . . , k n vol and k 1 geo , k 2 geo , . . . , k n geo represent the volumetric and geometric scattering kernel for the 1st, 2nd, . . . ,nth pixel, respectively; and ρ 1 , ρ 2 ,.., ρ n is the observed band reflectance for the 1st, 2nd, . . . , nth pixel. The Equation (23) can be written as a matrix form in Equation (24) and can be further expressed as a simplified form in Equation (25). Then the kernel coefficients can be calculated as f iso f vol f geo

BRDF Effects Correction for Image
The anisotropy factor (ANIF) is a ratio describing the relationship between directional reflectance and reference reflectance, and it is introduced to the BRDF correction model of the hyperspectral image. We followed the approach described in [31], and defined the ANIF with respect to the nadir reflectance as in Equation (26); therefore, the RT-BRDF corrected reflectance (ρ RT−BRDF ) can be calculated using Equation (27): nadir; θ s_ f ixed and ∆∅ f ixed represent the fixed solar position; and ρ BOA is the ATCOR-corrected apparent directional reflectance as mentioned in Section 2.2.

RT-BRDF Performance Assessment
In order to evaluate the performance of the RT-BRDF method, three assessment strategies were used qualitatively and quantitatively [7,53,86,87], which included (a) comparisons of the results with other topographic correction approaches; (b) comparisons with the MODIS products; (c) utilizing the fitted RT-BRDF correction coefficient to another larger area for the robustness verification. The specific operations are as follows.
(a) The results of three commonly used topographic correction approaches (i.e., C, SCS, and SCS+C, see Table 2) were compared to RT-BRDF. Specifically, we implemented the following goals: • First, since an effective correction method can produce a highly homogeneous mosaic imagery without topographic and BRDF (solar-viewing geometry) effects, a visual inspection was performed focusing on the topographic correction results and the variability of the overlapped area of adjacent flight-line images.

•
Previous studies [49,66,88] observed a positive linear relationship between surface reflectance and illumination conditions. An effective BRDF correction method can remove this relationship and result in the R2 and slope value of the linear function approaching zero. Table 2. Topographic correction methods for RT-BRDF performance assessment.

Term Equation 2 Description
Notes: 1 C is the ratio of p to q, and the it can be obtained from ρ BOA ≈pcos( θ s ) + q. 2 In these topographic correction methods, cos( θ s ) is calculated with the real slope α as input.
• The terrain-induced canopy BRDF can lead to reflectance variations with slope and aspect (i.e., topographic conditions). Practically, sun-facing pixels would potentially have higher reflectance than pixels lying on shaded slopes. An effective BRDF correction method can weaken the influence of topographic conditions, resulting in a more homogeneous spectral balance between sunlit and shaded areas.

•
Based on the assumption of homogeneous structure for each forest type, we followed the approach described in [87] focusing on the reflectance over the aspect of pixels perpendicular to the solar direction. An effective correction method can obtain a nearly perfect fit between uncorrected and corrected reflectance of these pixels.

•
Finally, this mitigation effect can be further confirmed by the analysis of the RMSE of the reflectance variation in the forested pixels (i.e., coniferous and broadleaved forest) from overlapping pixels of two adjacent hyperspectral images. In fact, a further effect of the BRDF correction is to decrease the standard deviation of the reflectance spectra for each land-cover type [53]. The coefficient of variation (CV) defined in Equation (28) is utilized to evaluate the results obtained with the RT-BRDF method: where σ and mean are the standard deviation and mean reflectance values of the land-cover type, respectively.
(b) The weighting parameters of the MODIS BRDF/Albedo Model Parameters product were identified from the relatively pure forest pixels. . (c) In order to assess the repeatability of the algorithm, we applied the RT-BRDF method to a larger area (04 April 2014 flight campaign) to verify whether the BRDF effects of airborne hyperspectral imagery over forested areas with rugged topography can be corrected.

Preprocessing Results
The training and test data for the forest type classification were obtained from field work (Section 2.3) and visual interpretation of high resolution CCD images. The overall classification accuracy and kappa coefficient of the airborne hyperspectral data are 94.9% and 0.93, respectively. Theoretical (i.e., flat topography hypothesis) and calculated viewing and illumination geometry from a selected region of the hyperspectral images are shown in Figures 4a-d and 4e-h, respectively. When compared to the hypothetical case, the slope and aspect of the surveyed topography lead to a wider range of local viewing and illumination geometry and consequently a multi-angular observation data set for the same land-cover type.
accuracy and kappa coefficient of the airborne hyperspectral data are 94.9% and 0.93, respectively. Theoretical (i.e., flat topography hypothesis) and calculated viewing and illumination geometry from a selected region of the hyperspectral images are shown in Figure 4a-d and Figure 4e-h, respectively. When compared to the hypothetical case, the slope and aspect of the surveyed topography lead to a wider range of local viewing and illumination geometry and consequently a multi-angular observation data set for the same land-cover type.

Optimal Pixel Numbers for BRDF Modeling
In this study, the optimal number of pixels for the BRDF modeling is defined as the minimum number of pixels, and was selected using a stratified random sampling method, which is necessary to achieve minimum RMSE between observed and model simulated reflectance. The calculated RMSE for a reference near infrared (NIR) band (i.e., band 94, 837.19 nm) is shown in Figure 5a,b for coniferous and broadleaved forest, respectively. The results in Figure 5a,b show that the RMSE declines dramatically with the number of pixels increasing from 100 to 10,000, and especially below 1000. In contrast, the RMSE remains low with increasing number of pixels from 10,000 to 50,000, and especially greater than 35,000. For both coniferous and broadleaved forest, for each band, the independent band normalized RMSE with number of pixels from 100 to 50,000 is shown in Figure 5c,d. For band 94, the trend of RMSE (Figure 5a,b) can therefore be represented as a consistent trend of independent band normalized RMSE (Figure 5c,d). The independent band normalized RMSE show similar overall trends for the majority of bands and forest types: (i) a rapid decline of RMSE with number of pixels increasing from 100 to 10,000; (ii) stabilization with number of pixels greater around 10,000; and (iii) consistently low RMSE value (i.e., <0.001) with number of pixels greater than 35,000. Ultimately, in this study the number of pixels n is set at 35,000 to ensure a comprehensive coverage of the viewing and illumination angles and higher fitting accuracy of the BRDF model for all bands. independent band normalized RMSE (Figure 5c,d). The independent band normalized RMSE show similar overall trends for the majority of bands and forest types: (i) a rapid decline of RMSE with number of pixels increasing from 100 to 10,000; (ii) stabilization with number of pixels greater around 10,000; and (iii) consistently low RMSE value (i.e., < 0.001) with number of pixels greater than 35,000. Ultimately, in this study the number of pixels n is set at 35,000 to ensure a comprehensive coverage of the viewing and illumination angles and higher fitting accuracy of the BRDF model for all bands.  Figure 6 shows a false-color composite (i.e., 837.19, 694.91, and 543.72 nm band for R, G, and B, respectively) close-up of the pre-(i.e., no terrain induced atmospheric effects corrected) and postcorrection hyperspectral images. The overlaps between adjacent hyperspectral images and the final mosaic outline is displayed in Figure 6a. At a first visual inspection, it can be observed that the reflectance in the BRDF uncorrected mosaic (Figure 6b) is significantly affected by topographic effects (e.g., shaded pixels darker than illuminated pixels). However, it can be noticed that the RT-BRDF, C, and SCS+C correction method effectively reduce these topographic effects as shown in Figure 6c,  Figure 6 shows a false-color composite (i.e., 837.19, 694.91, and 543.72 nm band for R, G, and B, respectively) close-up of the pre-(i.e., no terrain induced atmospheric effects corrected) and post-correction hyperspectral images. The overlaps between adjacent hyperspectral images and the final mosaic outline is displayed in Figure 6a. At a first visual inspection, it can be observed that the reflectance in the BRDF uncorrected mosaic (Figure 6b) is significantly affected by topographic effects (e.g., shaded pixels darker than illuminated pixels). However, it can be noticed that the RT-BRDF, C, and SCS+C correction method effectively reduce these topographic effects as shown in Figure 6c, Figure 6d,f, whereas the SCS correction in Figure 6e led to an obvious overcorrection in the weakly illuminated areas. Figure 6d,f also shows that both C and SCS+C corrected images still present a slight overcorrection in the backscattering regions (i.e., same observing direction and solar incidence direction [91] for each flight-line image). The full-sized RT-BRDF corrected and uncorrected images are displayed in the next with selected areas magnified to show significant differences (cfr. Figure 14).

BRDF Correction of Topographic Effect
To quantitatively analyze the topographic correction effects under varying illumination conditions, the relationship between reflectance of a selected NIR band (i.e., 837.19 nm) and the cosine of the local solar incidence angle (i.e., illumination condition) is shown in Figure 7, where 500 non-shadowed pixels were randomly selected from the hyperspectral mosaic for both coniferous and broadleaved forest. Figure 7a Figure 6e led to an obvious overcorrection in the weakly illuminated areas. Figure 6d,f also shows that both C and SCS+C corrected images still present a slight overcorrection in the backscattering regions (i.e., same observing direction and solar incidence direction [91] for each flight-line image). The full-sized RT-BRDF corrected and uncorrected images are displayed in the next with selected areas magnified to show significant differences (cfr. Figure  14).       Figure 8 shows the variation of reflectance at 837.19 nm with slope and aspect (i.e., topographic conditions). Figure 8a,f show an overall higher uncorrected reflectance of the broadleaved forest pixels in a great majority of topographic conditions compared to the uncorrected reflectance of coniferous forest pixels. Consistently with the results in Figure 8a,f, sun-facing aspect classes (i.e., high illumination conditions) display higher uncorrected reflectance values for both forest types. Conversely, the RT-BRDF, C, and SCS+C corrected reflectance of both coniferous and broadleaved forest pixels show a significant reduction of this positive relationship with topographic conditions. However, the SCS correction resulted in higher reflectance values for the shaded aspect classes (i.e., overcorrection; Figure 8d,i).
To assess the physical soundness of the results in Figure 8, we analyzed the relationship between uncorrected and corrected pixel reflectance at 837.19 nm and aspect perpendicular to the solar principle plane (i.e., approximately 10.8 • and 190.8 • for the first flight-line (06 April 2014)). Figure 9 confirms how pixels classified as broadleaved forest (Figure 9e-h) are generally characterized by higher reflectance than coniferous forest (Figure 9a-d). In addition, Figure 9 shows how the patterns of the relationships between corrected and uncorrected reflectance are similar in both forest types. The R 2 , RMSE, bias, and linear fitting function were analyzed to evaluate the physical soundness of correction methods. To assess the physical soundness of the results in Figure 8, we analyzed the relationship between uncorrected and corrected pixel reflectance at 837.19 nm and aspect perpendicular to the solar principle plane (i.e., approximately 10.8° and 190.8° for the first flight-line (06 April 2014)). Figure 9 confirms how pixels classified as broadleaved forest (Figure 9e-h) are generally characterized by higher reflectance than coniferous forest (Figure 9a-d). In addition, Figure 9 shows how the patterns of the relationships between corrected and uncorrected reflectance are similar in both forest types. The R 2 , RMSE, bias, and linear fitting function were analyzed to evaluate the physical soundness of correction methods. For the coniferous forest, RT-BRDF shows the highest R 2 (0.97), lowest RMSE

BRDF Correction Assessment for Multiple-Flights
Although the C and SCS+C BRDF correction methods demonstrated to be able to effectively remove the topographic effects (Figure 9), RT-BRDF achieves a better correction of the BRDF effects between multiple flight-lines. For example, Figure 10a,b shows the difference between RT-BRDF corrected and C corrected reflectance at 837.19 nm for the third and second flight-line (06 April 2014), respectively. The RT-BRDF correction significantly reduced the back-scattering effects (i.e., west of the flight-line), whereas forward-scattering (i.e., east of the flight-line) reflectance is increased or maintained. In addition, differences in reflectance values are noticeable in the mosaicked difference image (Figure 10c) for regions overlapping adjacent flight-lines. Figure 10 provides an example of false color composite close ups of the results obtained with the different correction methods. Although SCS shows a tendency to overcorrect (i.e., bright pixels), all other correction methods are capable of removing the terrain-induced canopy BRDF distortions with varying degrees of success. However, the edge of a hyperspectral image is clearly noticeable in the uncorrected, C and SCS+C

BRDF Correction Assessment for Multiple-Flights
Although the C and SCS+C BRDF correction methods demonstrated to be able to effectively remove the topographic effects (Figure 9), RT-BRDF achieves a better correction of the BRDF effects between multiple flight-lines. For example, Figure 10a,b shows the difference between RT-BRDF corrected and C corrected reflectance at 837.19 nm for the third and second flight-line (06 April 2014), respectively. The RT-BRDF correction significantly reduced the back-scattering effects (i.e., west of the flight-line), whereas forward-scattering (i.e., east of the flight-line) reflectance is increased or maintained. In addition, differences in reflectance values are noticeable in the mosaicked difference image (Figure 10c) for regions overlapping adjacent flight-lines. Figure 10 provides an example of false color composite close ups of the results obtained with the different correction methods. Although SCS shows a tendency to overcorrect (i.e., bright pixels), all other correction methods are capable of removing the terrain-induced canopy BRDF distortions with varying degrees of success. However, the edge of a hyperspectral image is clearly noticeable in the uncorrected, C and SCS+C corrected mosaics, indicating a less than optimal compensation for the BRDF effects in overlapping regions. In contrast, the RT-BRDF corrected mosaic image shows a higher spectral homogeneity. This mitigation effect is further confirmed by the analysis of the RMSE of the variation in reflectance in the forested pixels (i.e., coniferous and broadleaved forest) from the overlap of two adjacent hyperspectral images. It clearly appears that the RMSE decreased post-RT-BRDF ( Figure  11a,b), especially in the NIR region (i.e., >750 nm), implying that the selected forested areas shows more similar spectral values in the overlapped area of two images, and therefore an improvement in spectral homogeneity. The average decrease in RMSE of the spectra is 0.004 and 0.016 for coniferous This mitigation effect is further confirmed by the analysis of the RMSE of the variation in reflectance in the forested pixels (i.e., coniferous and broadleaved forest) from the overlap of two adjacent hyperspectral images. It clearly appears that the RMSE decreased post-RT-BRDF (Figure 11a,b), especially in the NIR region (i.e., >750 nm), implying that the selected forested areas shows more similar spectral values in the overlapped area of two images, and therefore an improvement in spectral homogeneity. The average decrease in RMSE of the spectra is 0.004 and 0.016 for coniferous forest and broadleaved forest, respectively.
Remote Sens. 2020, 12, 432 20 of 33 A post-RT-BRDF reduction of the CV can also be observed for both coniferous (Figure 11c) and broadleaved forest (Figure 11d) at every wavelength. Specifically, the average decrease in CV is 3% and 3.5% for coniferous forest and broadleaved forest, respectively. In addition, it can be noticed that this reduction is more marked at longer wavelengths (i.e., NIR), and in particular in correspondence and after the characteristic green reflectance peak of vegetation (i.e., from ≈490 nm), with the exception of the maximum chlorophyll absorption wavelength (i.e., ≈680 nm; [7]).  Figure 12a,b show the different sampling strategies of a multi-angle dataset for retrieving the kernel-based weighting parameters of AISA Eagle II and MODIS data, respectively. Then Figure 12ce shows the comparison of the VHR pixel-based three BRDF model weighting parameters ( , and ) of AISA Eagle II bands (400-990 nm) with the coarse scale (500 m) weighting parameters of MODIS bands 1-4, respectively. According to the strategies described in Section 3, the BRDF model weighting parameters of AISA Eagle II data are plotted for broadleaved forest and coniferous, respectively. Notably, all the BRDF model weighting parameters of broadleaved are slightly higher than those of coniferous in the visible (VIS) region but show much difference in the NIR region. The numerical ordering of the three weighting parameters is > > in each AISA Eagle II band for broadleaved forest and coniferous, respectively. A set of BRDF model weighting parameters for the corresponding MODIS bands 1-4 are generated, respectively, based on the relatively pure A post-RT-BRDF reduction of the CV can also be observed for both coniferous (Figure 11c) and broadleaved forest (Figure 11d) at every wavelength. Specifically, the average decrease in CV is 3% and 3.5% for coniferous forest and broadleaved forest, respectively. In addition, it can be noticed that this reduction is more marked at longer wavelengths (i.e., NIR), and in particular in correspondence and after the characteristic green reflectance peak of vegetation (i.e., from ≈490 nm), with the exception of the maximum chlorophyll absorption wavelength (i.e., ≈680 nm; [7]). Figure 12a,b show the different sampling strategies of a multi-angle dataset for retrieving the kernel-based weighting parameters of AISA Eagle II and MODIS data, respectively. Then Figure 12c-e shows the comparison of the VHR pixel-based three BRDF model weighting parameters ( f iso , f vol and f geo ) of AISA Eagle II bands (400-990 nm) with the coarse scale (500 m) weighting parameters of MODIS bands 1-4, respectively. According to the strategies described in Section 3, the BRDF model weighting parameters of AISA Eagle II data are plotted for broadleaved forest and coniferous, respectively. Notably, all the BRDF model weighting parameters of broadleaved are slightly higher than those of coniferous in the visible (VIS) region but show much difference in the NIR region. The numerical ordering of the three weighting parameters is f iso > f geo > f vol in each AISA Eagle II band for broadleaved forest and coniferous, respectively. A set of BRDF model weighting parameters for the corresponding MODIS bands 1-4 are generated, respectively, based on the relatively pure pixels with forest in the study area. The numerical ordering of the three weighting parameters is f iso > f vol > f geo in each MODIS band. Generally, the f iso of AISA Eagle II band is higher than the f iso of MODIS band.

MODIS BRDF Product Comparison
Compared with the f geo of AISA Eagle II band, the corresponding f geo of MODIS band is slightly lower in the VIS region and significantly drops in the NIR region to be close to zero. Meanwhile, f vol of AISA Eagle II band is nearly close to zero and much lower than the corresponding f vol of MODIS band. The three weighting parameters present significant bias between AISA Eagle II and MODIS data in the NIR region but these are interestingly corrected.
Remote Sens. 2020, 12, 432 21 of 33 pixels with forest in the study area. The numerical ordering of the three weighting parameters is > > in each MODIS band. Generally, the of AISA Eagle II band is higher than the of MODIS band. Compared with the of AISA Eagle II band, the corresponding of MODIS band is slightly lower in the VIS region and significantly drops in the NIR region to be close to zero. Meanwhile, of AISA Eagle II band is nearly close to zero and much lower than the corresponding of MODIS band. The three weighting parameters present significant bias between AISA Eagle II and MODIS data in the NIR region but these are interestingly corrected. As illustrated in Figure 13a,c,d, both uncorrected and RT-BRDF corrected AISA Eagle II aggregated reflectance show similar correlations with the MODIS nadir BRDF-adjusted reflectance. In addition, these figures show both the uncorrected and RT-BRDF corrected AISA Eagle II aggregated reflectance are generally higher or lower than the MODIS reflectance. This may be caused by the lower signal-to-noise ratio of airborne hyperspectral images in the blue band or other upscaling issues which make reconstructing the BRDF accurately more challenging. Notably, the AISA Eagle II RT-BRDF corrected aggregated reflectance (Figure 13b) are more consistent with the MODIS nadir BRDF-adjusted reflectance in the NIR (MODIS band2; R 2 = 0.83, RMSE = 0.007) showing a slightly higher consistency than the uncorrected scenario (R 2 = 0.69, RMSE = 0.007). As illustrated in Figure 13a,c,d, both uncorrected and RT-BRDF corrected AISA Eagle II aggregated reflectance show similar correlations with the MODIS nadir BRDF-adjusted reflectance. In addition, these figures show both the uncorrected and RT-BRDF corrected AISA Eagle II aggregated reflectance are generally higher or lower than the MODIS reflectance. This may be caused by the lower signal-to-noise ratio of airborne hyperspectral images in the blue band or other upscaling issues which make reconstructing the BRDF accurately more challenging. Notably, the AISA Eagle II RT-BRDF corrected aggregated reflectance (Figure 13b) are more consistent with the MODIS nadir BRDF-adjusted reflectance in the NIR (MODIS band2; R 2 = 0.83, RMSE = 0.007) showing a slightly higher consistency than the uncorrected scenario (R 2 = 0.69, RMSE = 0.007).

Figure 13.
Comparison between AISA Eagle II aggregated reflectance and MODIS nadir BRDFadjusted reflectance corresponding band1-4 (a-d), respectively . The best linear fit is in solid red/black, and the 1:1 line is in dashed black. Figure 14 shows the comparison of uncorrected and RT-BRDF corrected hyperspectral images of the 06 April 2014 flight campaign. Based on the full-sized observation of the image and the magnification of the selected areas, the RT-BRDF approach has the capability to reduce the BRDF effects caused by the combination of wide FOV airborne scanner and rugged topography in most forest covered areas. In addition, the RT-BRDF approach was utilized in another larger area hyperspectral images of the 04 April 2014 flight to verify the repeatability of the method. A crosstrack brightness gradient is clearly discernible in the uncorrected data (Figure 15a) but it is almost entirely absent in the RT-BRDF corrected data (Figure 15b). The close-up areas (site A, B, and C) in Figure 15 show the goodness of the correction in greater detail. Figure 13. Comparison between AISA Eagle II aggregated reflectance and MODIS nadir BRDF-adjusted reflectance corresponding band1-4 (a-d), respectively. The best linear fit is in solid red/black, and the 1:1 line is in dashed black. Figure 14 shows the comparison of uncorrected and RT-BRDF corrected hyperspectral images of the 06 April 2014 flight campaign. Based on the full-sized observation of the image and the magnification of the selected areas, the RT-BRDF approach has the capability to reduce the BRDF effects caused by the combination of wide FOV airborne scanner and rugged topography in most forest covered areas. In addition, the RT-BRDF approach was utilized in another larger area hyperspectral images of the 04 April 2014 flight to verify the repeatability of the method. A cross-track brightness gradient is clearly discernible in the uncorrected data (Figure 15a) but it is almost entirely absent in the RT-BRDF corrected data (Figure 15b). The close-up areas (site A, B, and C) in Figure 15 show the goodness of the correction in greater detail.

Discussion
The past two decades have seen significant advances in the use of airborne remote sensing data in the field of Earth system science and for ecosystem process modeling, from the canopy/stand level to more recently larger regions, increasing the demand for hyperspectral and LiDAR measurements. However, only few studies explored the applicability of BRDF correction methods to airborne hyperspectral imagery [6,7,[92][93][94], especially for the forested areas with rugged topography. A better understanding of the BRDF effects on airborne hyperspectral imagery will provide clearer insights into their influence on remote sensing-based quantitative inversion products (e.g., LAI, PRI, canopy biochemistry parameters). In this study, we proposed an improved correction method for airborne hyperspectral imagery over rugged terrain with forest.

Inter-Comparison of BRDF Correction Methods
Several studies demonstrated how airborne imagery with high spatial resolution (i.e., 0.6-4.5 m) are significantly affected by BRDF effects and can only partly follow the standardized BRDF processing methods of satellite images (e.g., [7,15,42,95]). The results of our study further confirmed that slope-based volumetric and geometric kernels, which correct for both incident-angle-dependent and observation-angle-dependent BRDF effects, can be suitable for BRDF correction [47,48,76,77,96]. Especially, the improved RTM kernels (see Equation (11)) and LTR (see Equation (13)) demonstrated to be successful in BRDF correction for high resolution airborne imagery of forested areas with both flat [7] and rugged topography. With the exception of SCS, the topographic effects are effectively removed by RT-BRDF, C, and SCS+C (see . In fact, the SCS, which is a semi-physical method, failed in topographic correction under diffuse light conditions with an overcorrection of the reflectance [48,87,97]. Although the empirical, parameter-based C and SCS+C methods proved their effectiveness for topographic correction of satellite image, it remains challenging to utilize a universal empirical parameter C to correct images with varying topography and viewing-illumination geometries [48,87,[98][99][100]. Both a visual inspection ( Figure 6) and quantitative analysis ( Figure 10) confirm that the empirical parameter-based methods are less than optimal for the BRDF correction of airborne images, especially for multiple flight-lines. In contrast, our results show that the physics-based RT-BRDF correction method can not only consistently compensate for terrain-induced but also for the viewing-illumination induced BRDF effects (Figure 11), which is supported by physical explanations (Figures 7 and 9). Moreover, despite previous research based on the flat terrain hypothesis (e.g., [20,101,102]), it is demonstrated here the potential to correct the BRDF effects between multiple flight-lines images. The RT-BRDF approach opens the possibility of failures due to insufficient observations and/or poor angular sampling capability (e.g., [95]). Accounting for topography in the BRDF correction method still remains a largely unexplored task [7,19,42], but by design the limitations above do not affect the RT-BRDF method (Figures 4 and 5).

Kernel-Driven BRDF Model Applied in Different Spatial Scales
Despite the spatial scaling assumptions underlying the semi-empirical kernel-driven BRDF approach for low spatial resolution [103], many studies demonstrated that it can be successful for reducing the BRDF effects of VHR airborne imagery [7,15,22,31,42]. However, only a few of these studies (i.e., [15]) give reasonable physical explanations for the applicability of the kernel-driven BRDF model with high spatial resolution imagery, which still needs further exploration. Results of Figure 12 show that for a coarse spatial resolution pixel such as 500 m of MODIS is dominated by the isotropic and volumetric scattering components. Likely, the geometric scattering component is inappropriate for the relief with a misrepresentation of shadowing effects in the observed scenes. It is reasonable that a MODIS pixel at this coarse-scale might contain many dense and continuous canopies in this study. Furthermore, a MODIS pixel contains both terrain-induced sunlit and shaded slope effects affecting the geometrical optics kernel in the kernel-driven BRDF modeling process (Figure 12b). Meanwhile, Figure 12 also reveals that a high spatial resolution pixel, such as for the VHR airborne imagery in this study (spatial resolution of 1 m), is more impacted by the isotropic and geometric components and less by the volumetric scattering component. A pixel of airborne image with ≈1 m resolution likely contains individual leaves and branches, and these physical attributes within the VHR pixels will lead to structure-induced mutual shadows within the canopy [15]. Moreover, topographic information, that differs from MODIS, was added to our RT-BRDF model which contributes to more geometric scattering effects of airborne imagery. It should be noticed that the cast-shaded pixels (deep gullies and canopy gaps) were ignored for the BRDF model reconstruction of airborne imagery in this study which could be the reason why f iso of AISA Eagle II band has a higher value than f iso of MODIS band. Ultimately, our results are consistent with other studies [43,104,105] suggesting that different scales can lead to different spatial and structural patterns containing different shadowing and clumping effects affecting the corresponding scattering component in the kernel-driven BRDF model.
Airborne hyperspectral images after RT-BRDF correction can reduce the reflectance in the back-scattering regions (or sunlit slope, if there is a mountain) and raise the reflectance in the forward-scattering regions (or shaded slope, if there is a mountain). While the RT-BRDF corrected airborne imagery performs better than the uncorrected images, there is no clear difference in Figure 13a,b. These results may be caused by the pixel aggregation process (500 m × 500 m) which averages the sunlit and shaded slope effects (swath width with ≈1000 m for one AISA Eagle II flight-line in this study). We therefore indicate that the airborne hyperspectral images under the acquisition conditions of this study have limited BRDF effects at course spatial resolutions (e.g., 500 m).

The Role of a Priori Knowledge on the BRDF Correction
Previous studies demonstrated that different land-cover types have different anisotropic properties and therefore different intrinsic BRDF characteristics [22,31,38,43]. Furthermore, the BRDF effects can vary with tree structure even within the same forest type [31,91]. While numerous studies use classification maps as a priori knowledge in the BRDF correction, this is by no means the only possible input [6,7,22]. A broader definition of a priori knowledge should contain external ancillary resources, such as field data, classification map, structure information, BRDF shapes from other products, topographic information, and ecological phenological season [38,86,106,107]. For example, for VHR airborne images, the structural information can be represented as texture based in algorithms such as GLCM, Sum Difference Histogram (SADH), and FOurier-based Textural Ordination (FOTO) [108,109]. Furthermore, the forest vertical and spatial structure can be derived from LiDAR data, such as canopy height, crown shape, as well as point cloud distribution of individual trees [110]. In this work, the general structure of the forest stand was investigated in the field, the classification map was produced using the SVM classification algorithm, and the high resolution DEM was derived directly from airborne LiDAR. Our results show that forest-type based reflectance is significantly different even for similar topographic conditions (Figure 8a,f), and as a consequence the BRDF models must be built separately. Due to the different structure of the coniferous and broadleaved forests in this study, it is possible to build accurate forest type-based BRDF models ( Figure 5), leading to improved forest type-based correction results (Figures 7 and 9). This further reinforces the fundamental role of a priori knowledge to describe the surface anisotropy based on multidimensional and temporal information in the BRDF correction.

Future Developments
Despite the RT-BRDF correction method for airborne hyperspectral imagery of forested areas with rugged topography was demonstrated to be successful across two study areas (Figures 14 and 15), more investigations are needed to analyze how each factor affects the performance. Furthermore, additional research is needed to evaluate for sparse canopies and other land-cover types for different ecosystems. Although a combination of the RTM and LTR fixed kernels was chosen in our study, recent studies demonstrate that an adaptive kernels selection approach can improve the accuracy of BRDF models for different land-cover/structural types compared to fixed kernels [23,31]. In addition, the fitting method for BRDF modeling might also affect the correction result as mentioned in [22], especially for non-forest cover types such as wetlands and snow covered surfaces.
A typical hyperspectral flight has duration of 4-5 h and is influenced by varying sky conditions (i.e., the ratio of direct and diffuse radiation). However, because of the challenge to accurately describe the contribution of direct and diffuse irradiance to total irradiance in the actual operating process [32,33], in the RT-BRDF model we simplified this contribution by neglecting the influence of diffuse and adjacent surrounding radiance in the correction process. Another critical aspect of our study is that the CAF-LiCHy system was simultaneously acquired for both the airborne hyperspectral imagery and the LiDAR data, enabling the creation of an accurate high resolution LiDAR-derived DEM to describe the canopy and the underlying topography. In addition, the sensitivity analysis of the topographic correction with the DEM at different spatial scales highlighted the need for highly detailed DEMs to describe the physically based BRDF mechanism over rugged terrain [49,54]. This can be achieved using airborne systems integrating hyperspectral and LiDAR sensors. However, although the extremely high detail level of the DEM can trigger the characterization of fine scale topographic elements such as deep gullies and crags, it remains challenging to physically remove the shadow/cast-shadow effects due to varying irradiance.
Quantitative research on the correction of remotely sensed reflectance over rugged terrain has recently gained increasing attention, with the creation of new models such as the Kernel-Driven reflectance model for Sloping Terrain (KDST; [47]) and the Geometric-Optical model for Sloping Terrains (GOST; [111][112][113]). In addition, recent studies have worked towards a more accurate description of direct and diffuse irradiance components at the pixel level by considering a physically based approach with the combination of fine resolution DEM/DSM and the Discrete Anisotropic Radiative Transfer (DART) model [33,114].

Conclusions
In this study, we presented an improved kernel based BRDF correction operational framework for airborne hyperspectral images of forested areas with rugged topography (RT-BRDF). The RT-BRDF method is based on the multi-angular nature of aerial hyperspectral observations, caused by the combination of sensor scanning angles and topographic effects. In particular, the RT-BRDF method takes into account both the BRDF effects derived from topographic conditions and observation/solar incident angles, with the advantage of providing a physical meaning for the observed BRDF effects, differently from empirical and semi-physical topographic correction models. Additionally, a stratified random sampling method for the VHR airborne image with a large number of pixels can be an efficient way to ensure good representative of viewing and illumination angles through the minimum fitting error of the BRDF model for all bands (35,000 in this study). We demonstrated that, compared to other correction methods, RT-BRDF reduced the BRDF effects from multi-line aerial hyperspectral images of forested areas with rugged topography, and we believe that land-cover based models obtained using RT-BRDF have the potential to be applicable to large areas comprising numerous flight-lines. Results of this study can be used for reference in follow-up research and would contribute to reduce the limitation of large-scale applications of airborne hyperspectral data.