A Flight Direction Design Method for Airborne Spectral Imaging Considering the Anisotropy Reflectance of the Target in Rugged Terrain

An excellent mission plan is the prerequisite for the acquisition of high quality airborne hyperspectral images which are useful for environmental research, mining etc. In order to minimize the radiance non-uniformity caused by the anisotropic reflectance of targets, the flight direction is mostly designed on the solar azimuth or 180° from it for whiskbroom and pushbroom imagers. However, the radiance to the observer is determined not only by the reflectance of the target, but also by the terrain, the illumination direction and the observation direction. So, the flight direction which is defined to minimize radiance non-uniformity might change with the terrain. In order to find the best flight direction for rugged terrain, we firstly analyze the causes of the effect of terrain on radiation non-uniformity based on the radiative transfer process. Then, the flight direction design method is proposed for composite sloping terrain. Tested by digital and physical simulation experiments, the radiance non-uniformity is minimized when the aircraft flies in the designated direction. Finally, a workflow for flight direction planning and optimizing is summarized, considering the flight mission planning techniques and the workflow of remote sensing missions.


Introduction
Due to the high spectral resolution of hyperspectral images, a large number of applications requiring fine identification of materials or estimation of physical parameters are enabled [1]. As the main equipment for acquiring hyperspectral images, airborne hyperspectral sensors have been widely used in recent decades. Whiskbroom and pushbroom hyperspectral imagers such as AVIRIS [2], HyMap [3], CASI/SASI [4,5], APEX [6] have been adopted for many operational routines. An excellent flight plan is the prerequisite of high imaging quality [7]. Firstly, to guarantee the proper ground sampling distance (GSD), the absolute flight altitude and speed should be calculated with the instantaneous field of view (IFOV) of the imager and the digital elevation model (DEM) of the area of interest (AOI) [8,9]. Secondly, as the AOI cannot be covered with a single flight stripe most of time, several parallel stripes are necessary for the entire flight routine to cover the entire AOI [10]. Then, the proper distance between adjacent stripes is calculated with the field of view (FOV) of imager, the flight height and the DEM of the AOI to guarantee specific overlaps between adjacent stripes [11]. Thirdly, to ensure illumination onto the scene and to limit the presence of shadows, the image is acquired on a sunny day when solar zenith angle is less than 60 • [9].
The effect of the bidirectional reflectance of the ground surface is also considered for some missions. Because of the anisotropy in the reflectance of the surface and the difference of observation directions  [24]. Table 1. Notations used for the description of radiance transfer process.

Symbols Explanations
( ) Figure 1. The geometric model for illumination and observation [24]. Table 1. Notations used for the description of radiance transfer process.

Symbols Explanations
→ T (θ T , φ T ) Normal vector of terrain in horizontal coordinate system → I (θ i , φ i ) Illumination vector (from sun to target) in horizontal coordinate system → E (θ e , φ e ) Observation vector (from observer to target) in horizontal coordinate system θ T The slope in horizontal coordinate system φ T The aspect in horizontal coordinate system θ i The solar zenith angle in horizontal coordinate system φ i The solar azimuth angle in horizontal coordinate system θ e The observation zenith angle in horizontal coordinate system φ e The observation azimuth angle in horizontal coordinate system θ i The incident zenith angle in local slope coordinate system φ i The incident azimuth angle in local slope coordinate system θ e The exit zenith angle in local slope coordinate system φ e The exit azimuth angle in local slope coordinate system The incident vector and exit vector in local slope coordinate system are acquired by rotating the coordinate system by θ T around Y-axis. Then, the new coordinate system is rotated by −φ T around the new Z-axis, as shown in Equations (1) and (2). R Z (φ T ) and R Y (θ T ) are the rotation matrixes shown as Equations (3) and (4).
Then the angles defining the incident direction and exit direction are shown as Equations (5) to (8).
θ i = arccos(sin θ T cos φ T sin θ i cos φ i + sin θ T sin φ T sin θ i sin φ i + cos θ T cos θ i ), (5) θ e = arccos(sin θ T cos φ T sin θ e cos φ e + sin θ T sin φ T sin θ e sin φ e + cos θ T cos θ e ), Sensors 2019, 19, 2715 4 of 20 φ e = arctan − sin θ e sin(φ T − φ e ) cos θ T sin θ e cos(φ T − φ e ) − sin θ T cos θ e , (8) The radiative transfer process is shown in Figure 2. For rugged terrain, the target is illuminated by the direct solar light (sunlight), downward scattered light (skylight) and background reflected light which are the three adding items in the bracket on the right side of Equation (9) [24]. Then the radiance reflected by the target transfers to the observer through the atmosphere, combining with the upwelled scattered radiance (path radiance), as shown in Equation (9) [24]. E sλ represents the exoatmospheric spectral solar irradiance coming from θ i , φ i . L dλ denotes the skylight radiance to the target from the direction σ, φ. L bλ denotes the radiance reflected from the background toward the target. F and 1-F are the shape factors denoting the fraction of the hemisphere above the target, which are sky and background respectively, as shown in Figure 2 [24]. L uλ denotes the path radiance toward the observer. τ 1 and τ 2 denote the transmittance from top of atmosphere (TOA) to the target and from target to the sensor. ρ(θ i , θ e , φ , λ) is the bidirectional reflectance factor (BRF) of the target in the local slope coordinate system.
Sensors 2017, 17, x FOR PEER REVIEW 5 of 21 The Hapke model is often adopted to simulate the BRDF of rock or desert [25]; as for vegetation, the SAIL model is often used [26]. The BRF of a surface can be expressed as its corresponding BRDF times π shown as Equation (11) [27]. Even though these models can simulate the BRF of targets precisely, many parameters are needed for a correct simulation. The fill factor, the phase function and the scattering albedo of sand particles are needed for Hapke model [28]. The leaf area, leaf inclination distribution, reflectance and transmittance of leaf, soil reflectance and illumination condition are needed for SAIL model [29]. However, these parameters are not always available. So, a semi-empirical kernel model for rock and vegetation is chosen to calculate the BRF of target [30], as shown in Equation (12 ( ) geo f λ Figure 2. The radiative transfer process [24].
The Hapke model is often adopted to simulate the BRDF of rock or desert [25]; as for vegetation, the SAIL model is often used [26]. The BRF of a surface can be expressed as its corresponding BRDF times π shown as Equation (11) [27]. Even though these models can simulate the BRF of targets precisely, many parameters are needed for a correct simulation. The fill factor, the phase function and the scattering albedo of sand particles are needed for Hapke model [28]. The leaf area, leaf inclination distribution, reflectance and transmittance of leaf, soil reflectance and illumination condition are needed for SAIL model [29]. However, these parameters are not always available. So, a semi-empirical kernel model for rock and vegetation is chosen to calculate the BRF of target [30], as shown in Equation (12). The volume kernel K vol (θ i , θ v , φ) and the geometric kernel K geo (θ i , θ v , φ) are functions of illumination direction and observation direction. f vol (λ) and f geo (λ) are the factors of kernels. ρ iso (λ) is the isotropic reflectance of the target. The RossThick and LiSparse kernels are used, as shown in Equation (13). The phase angle ζ is the angle between illumination direction and observation direction. t is a variable depending on particle shapes and the direction (θ i , θ e , φ ) [30].
According to the radiance transfer process, the radiance reaching the observer for each pixel can be simulated. The radiance gradient in the image and the radiance difference in the overlap area between adjacent stripes can be defined based on the radiance at nadir L λ (0, φ e , λ). Assuming the nadir direction of observer is normal to the earth, the maximum radiance difference in an image would appear at the edge of the image, with a view zenith angle of ± FOV 2 . For the flight direction φ F shown as Equation (14), the radiance gradient ∆L λ−intra (φ F , λ) is defined with the greatest radiance difference in the image as shown in Equation (15). As to the overlapping area, the observation zenith angles of adjacent stripes change slightly for different pixels. The angle between two observation directions comes to the maximum when the observation zenith angles are equal. For the extreme case, the radiance difference between the overlap area of two stripes ∆L λ−inter (φ F , λ) is defined for the condition that only the edge pixel overlaps, as shown in Equation (16). The BRF of sand is shown as Figure 3 [31]. In order to analyze the relationship between radiance difference caused by the BRF of target and the flight direction in the local slope coordinate system, the observation directions are firstly labeled in Figure 3 for specific flight direction. The flight direction is labeled as the blue arrow with a relative flight direction ∆φ F defined as Equation (17). φ F is the flight direction in the local slope coordinate system. Refering to Equation (2), φ F is calculated as Equation (18) by transforming the flight direction vector in horizontal coordinate system to the vector in slope local slope coordinate system. Then, the observation directions are red circles in Figure 3, c.f. Equation (14). It is obvious that the ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) come to minimum when ∆φ F is 0 • or 180 • . According to the radiance transfer process, the radiance reaching the observer for each pixel can be simulated. The radiance gradient in the image and the radiance difference in the overlap area between adjacent stripes can be defined based on the radiance at nadir . For the flight direction F φ shown as Equation (14), the radiance gradient is defined with the greatest radiance difference in the image as shown in Equation (15). As to the overlapping area, the observation zenith angles of adjacent stripes change slightly for different pixels. The angle between two observation directions comes to the maximum when the observation zenith angles are equal. For the extreme case, the radiance difference between the overlap area of two stripes defined for the condition that only the edge pixel overlaps, as shown in Equation (16). The BRF of sand is shown as Figure 3 [31]. In order to analyze the relationship between radiance difference caused by the BRF of target and the flight direction in the local slope coordinate system, the observation directions are firstly labeled in Figure 3 for specific flight direction. The flight direction is labeled as the blue arrow with a relative flight direction F φ Δ defined as Equation (17). F φ ′ is the flight direction in the local slope coordinate system. Refering to Equation (2), F φ ′ is calculated as Equation (18) by transforming the flight direction vector in horizontal coordinate system to the vector in slope local slope coordinate system. Then, the observation directions are red circles in Figure 3, c.f. Equation (14). It is obvious that the ( ) tan arctan cos   The best flight direction is defined as the φ F when ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) reaches a minimum. For horizontal surface, the X-Y-Z coordinate system is the same with the X'-Y'-Z' coordinate system. Then ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) come to minimum when the flight direction is parallel to the principle plane. This is consistent with the previous studies [15][16][17][18][19][20]. When the surface is not horizontal, the X-Y-Z coordinate system doesn't coincide with the X'-Y'-Z' coordinate system. Then ∆φ F changes with slope and aspect of the surface as well as the illumination direction. The best flight direction for different condition is calculated as Equation (19).
For more common conditions, there are several slopes in the AOI, whose slope and aspect angles are different from one another. Then, the best flight direction for each slope is not the same. In order to minimize the ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) for most pixels, the best flight direction for each pixel in AOI is firstly calculated, and then the probability density function (PDF) P(φ F ) is estimated with kernel density method [32]. At last, the best flight direction of AOI is defined as the expectation of the best flight direction for each pixel, shown as Equation (20).

The Solo Slope Digital Simulation Experiment
The relationship between the radiance gradient or the radiance difference in the overlap area with the slope of the surface was firstly analyzed with the solo slope digital simulation experiment. According to Equation (8), the radiance received by the sensor was simulated for different slopes and aspects of the surface, the illumination directions and observation directions. ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) were calculated with Equations (15) and (16), respectively. The sunlight irradiance onto the surface E sλ cos θ i τ 1 (λ) and the skylight radiance onto the surface L dλ (σ, ϕ) were calculated with MODTRAN 4.0 code [33]. Mid-Latitude Summer (MLS) and Rural were selected as the atmosphere model and aerosol model respectively. The other simulation parameters are listed in Table 2. For simplicity, the radiance reflected by background was ignored as the terrain doesn't change dramatically when the slope is less than the solar zenith angle. The target simulated is the gravel in the desert whose diameter is several millimeters. To simulate the reflectance from different illumination directions and observation directions, the kernel factors in Equation (11) were regressed with the BRF collected in a similar desert by using an ASD Fieldspec Pro spectroradiometer [34] between zenith of −40 • and +40 • with 5 • intervals in the principle plane and the cross-principle plane. Considering the commonly-used hyperspectral imager equipped on the airborne [1][2][3][4][5], the observation zenith was set to be ±30 • for the edge of FOV and 0 • for the nadir as equal to the FOV of HyMap imaging spectroradiometer [3].
For different slopes, ∆L λ−intra (φ F , λ) changes with solar azimuth and flight direction, as shown in Figure 4. The changing of rises with the increase in slope in general. For horizon surface whose θ T = 0 • , ∆L λ−intra (φ F , λ) comes to a minimum when φ F = φ i as shown by the yellow circles in Figure 4a,g. This phenomenon coincides with previous studies [15][16][17][18][19][20]. As θ i , θ e , φ changing with slope of the surface, the distribution of ∆L λ−intra (φ F , λ) for different φ F and φ i changes significantly. These phenomena show that the flight direction may not be the solar azimuth when surface is not horizontal. For example, as shown in Figure 4l, ∆L λ−intra (φ F , λ) is 23% when the flight direction is 90 • (or 270 • ) and solar azimuth is 94 • (or 266 • ). The best flight direction should be 200 • and 160 • when the solar azimuth angles are 94 • and 266 • respectively and the ∆L λ−intra (φ F , λ) would be less than 0.06. Calculated with Equation (19), the best fight directions for different azimuth are shown as the yellow dash lines and circles in Figure 4.       The above experiment didn't consider the changing of the solar zenith angle. The range of the solar zenith angle is affected by target location. The Huangshan copper mine in Xinjiang, China was chosen as the test area [35]. Huangshan copper mine is a typical many porphyry copper-gold mineralization subzone where there are many typical minerals. The latitude of Huangshan copper mine is 40 • , which leads to the solar zenith angle changing between 17.4 • and 46.9 • from summer morning to afternoon. In the meantime, the slope in this scene changes from 0 • to 36 • with an average of 23 • . Under this condition, the effect of shadow is avoided while the terrain is obviously rugged. The solar azimuth is listed in Table 2. ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) was calculated with the radiance simulated according to Equation (8). The ∆L λ−intra (φ F , λ) for different solar azimuth and flight direction is shown in Figure 5. The best flight directions were calculated with Equation (19) and are shown as the yellow circles in Figure 5. It is obvious that ∆L λ−intra (φ F , λ) comes to a minimum on the best directions. Therefore, for solo slope, the best flight directions can be designed with Equation (19).  . Radiance gradient in the image under different flight direction and solar azimuth: (a)-(f) the radiance gradient when aspect is 0° and slopes are 0° to 50° respectively, (g)-(l) the radiance gradient when aspect is 180°and slopes are 0° to 50° respectively.

Attribute Value
The above experiment didn't consider the changing of the solar zenith angle. The range of the solar zenith angle is affected by target location. The Huangshan copper mine in Xinjiang, China was chosen as the test area [35]. Huangshan copper mine is a typical many porphyry copper-gold mineralization subzone where there are many typical minerals. The latitude of Huangshan copper mine is 40°, which leads to the solar zenith angle changing between 17.4° and 46.9° from summer morning to afternoon. In the meantime, the slope in this scene changes from 0° to 36° with an average of 23°. Under this condition, the effect of shadow is avoided while the terrain is obviously rugged. The solar azimuth is listed in Table 2.  Figure 5. The best flight directions were calculated with Equation (19) and are shown as the yellow circles in Figure 5. It is obvious that comes to a minimum on the best directions. Therefore, for solo slope, the best flight directions can be designed with Equation (19).

The Composite Slope Digital Simulation Experiment
For more common conditions, an AOI is composed of many slopes; hence, there could be more than one slope in the image. The experiment was conducted to simulate the reflected radiance of the Huangshan copper mine for different illumination directions and observation directions. As the airborne hyperspectral image is affected by terrain and the anisotropic reflectance of surface, the reflectance retrieved from an airborne hyperspectral image contains complex non-uniformity. So, the simulation scene was built with the BRF of simulated rock samples shown in Figure 6 according to the classification map [36], instead of simply using the airborne image. The classification map was developed according to the lithological and mineral distribution in the area [37]. The simulated rock samples were made by mixing mineral powders. And the proportion of mineral powders was similar to the component of rock samples collected in the field. The BRFs of the samples were measured with an ASD FieldSpec Pro spectroradiometer in the laboratory. The DEM of the experiment scene was described by ASTER GDEM data of the AOI. The sunlight irradiance and skylight radiance were calculated with MODTRAN 4.0 model with the parameters listed in Table 3. The radiance reflected by background was ignored. The solar zenith and azimuth angles were set to be consistent with those used in solo slope experiment. The at-sensor radiance images were simulated with Equation (8), as shown in Figure 7.

The Composite Slope Digital Simulation Experiment
For more common conditions, an AOI is composed of many slopes; hence, there could be more than one slope in the image. The experiment was conducted to simulate the reflected radiance of the Huangshan copper mine for different illumination directions and observation directions. As the airborne hyperspectral image is affected by terrain and the anisotropic reflectance of surface, the reflectance retrieved from an airborne hyperspectral image contains complex non-uniformity. So, the simulation scene was built with the BRF of simulated rock samples shown in Figure 6 according to the classification map [36], instead of simply using the airborne image. The classification map was developed according to the lithological and mineral distribution in the area [37]. The simulated rock samples were made by mixing mineral powders. And the proportion of mineral powders was similar to the component of rock samples collected in the field. The BRFs of the samples were measured with an ASD FieldSpec Pro spectroradiometer in the laboratory. The DEM of the experiment scene was described by ASTER GDEM data of the AOI. The sunlight irradiance and skylight radiance were calculated with MODTRAN 4.0 model with the parameters listed in Table 3. The radiance reflected by background was ignored. The solar zenith and azimuth angles were set to be consistent with those used in solo slope experiment. The at-sensor radiance images were simulated with Equation (8), as shown in Figure 7.

The Composite Slope Digital Simulation Experiment
For more common conditions, an AOI is composed of many slopes; hence, there could be more than one slope in the image. The experiment was conducted to simulate the reflected radiance of the Huangshan copper mine for different illumination directions and observation directions. As the airborne hyperspectral image is affected by terrain and the anisotropic reflectance of surface, the reflectance retrieved from an airborne hyperspectral image contains complex non-uniformity. So, the simulation scene was built with the BRF of simulated rock samples shown in Figure 6 according to the classification map [36], instead of simply using the airborne image. The classification map was developed according to the lithological and mineral distribution in the area [37]. The simulated rock samples were made by mixing mineral powders. And the proportion of mineral powders was similar to the component of rock samples collected in the field. The BRFs of the samples were measured with an ASD FieldSpec Pro spectroradiometer in the laboratory. The DEM of the experiment scene was described by ASTER GDEM data of the AOI. The sunlight irradiance and skylight radiance were calculated with MODTRAN 4.0 model with the parameters listed in Table 3. The radiance reflected by background was ignored. The solar zenith and azimuth angles were set to be consistent with those used in solo slope experiment. The at-sensor radiance images were simulated with Equation (8), as shown in Figure 7.  Two areas marked as zone 'A' and 'B' respectively in Figure 7a were defined to test the proposed method of flight direction design. The PDFs of slope and aspect in zones A and B are quite different, as shown in Figure 8. The average slope of zone A is greater than that of zone B (12.5° for zone A and 6.9° for zone B). The slopes in zone A are generally facing west. And most of the slopes in zone B are facing north.   Two areas marked as zone 'A' and 'B' respectively in Figure 7a were defined to test the proposed method of flight direction design. The PDFs of slope and aspect in zones A and B are quite different, as shown in Figure 8. The average slope of zone A is greater than that of zone B (12.5 • for zone A and 6.9 • for zone B). The slopes in zone A are generally facing west. And most of the slopes in zone B are facing north. ∆L λ−inter (φ F , λ) for different solar azimuths and flight directions were calculated with Equation (16), as shown in Figure 9. The best flight direction calculated with Equation (20) is shown as the red dash line in Figure 8. Comparing with the distribution of ∆L λ−inter (φ F , λ), the best flight direction is very close to the direction where ∆L λ−inter (φ F , λ) coming to the minimum.

The Composite Slope Physical Simulation Experiment
The BRF of the mineral powder samples was different from the field targets' BRF in the previous digital simulation experiment. A physical simulation experiment was conducted to simulate the real imaging process one step further. As shown in Figure 10, a physically simulated scene under the illumination from a physical solar simulator and a physical skylight simulator [37] were imaged from several directions. As the rock sample in the Huangshan cooper mine was not enough for the scene simulation, the Huangshan Dong cooper mine was chosen as the experiment scene. The scaled scene model shown in Figure 11 was developed according to the ASTER GDEM of the scene, the lithological and the mineral distribution map [38], using the rock samples collected from the field. As the BRF changed with the diameter of smashed samples, as shown in Figure 12, the diameter of smashed sample was chosen by comparing the BRF of smashed sample with the BRF of rock sample. The RMSE of BRF simulated is less than 0.01. The illumination onto the scene was simulated with the illumination direction, atmospheric and aerosol modes setting as listed in Table 4. The solar azimuth

The Composite Slope Physical Simulation Experiment
The BRF of the mineral powder samples was different from the field targets' BRF in the previous digital simulation experiment. A physical simulation experiment was conducted to simulate the real imaging process one step further. As shown in Figure 10, a physically simulated scene under the illumination from a physical solar simulator and a physical skylight simulator [37] were imaged from several directions. As the rock sample in the Huangshan cooper mine was not enough for the scene simulation, the Huangshan Dong cooper mine was chosen as the experiment scene. The scaled scene model shown in Figure 11 was developed according to the ASTER GDEM of the scene, the lithological and the mineral distribution map [38], using the rock samples collected from the field. As the BRF changed with the diameter of smashed samples, as shown in Figure 12, the diameter of smashed sample was chosen by comparing the BRF of smashed sample with the BRF of rock sample. The RMSE of BRF simulated is less than 0.01. The illumination onto the scene was simulated with the illumination direction, atmospheric and aerosol modes setting as listed in Table 4. The solar azimuth

The Composite Slope Physical Simulation Experiment
The BRF of the mineral powder samples was different from the field targets' BRF in the previous digital simulation experiment. A physical simulation experiment was conducted to simulate the real imaging process one step further. As shown in Figure 10, a physically simulated scene under the illumination from a physical solar simulator and a physical skylight simulator [37] were imaged from several directions. As the rock sample in the Huangshan cooper mine was not enough for the scene simulation, the Huangshan Dong cooper mine was chosen as the experiment scene. The scaled scene model shown in Figure 11 was developed according to the ASTER GDEM of the scene, the lithological and the mineral distribution map [38], using the rock samples collected from the field. As the BRF changed with the diameter of smashed samples, as shown in Figure 12, the diameter of smashed sample was chosen by comparing the BRF of smashed sample with the BRF of rock sample. The RMSE of BRF simulated is less than 0.01. The illumination onto the scene was simulated with the illumination direction, atmospheric and aerosol modes setting as listed in Table 4. The solar azimuth was simulated by changing the orientation of experiment scene. Finally, the experimental scene was observed at the directions listed in Table 4. Two of the acquired images are shown in Figure 13. was simulated by changing the orientation of experiment scene. Finally, the experimental scene was observed at the directions listed in Table 4. Two of the acquired images are shown in Figure 13.   was simulated by changing the orientation of experiment scene. Finally, the experimental scene was observed at the directions listed in Table 4. Two of the acquired images are shown in Figure 13.   was simulated by changing the orientation of experiment scene. Finally, the experimental scene was observed at the directions listed in Table 4. Two of the acquired images are shown in Figure 13.     (16), the image observed from two view zeniths must be geometrically corrected with reference to the nadir image. In order to validate the flight direction for different scenes, the two zones whose terrain was rugged and different from each other were defined corresponding to the A and B in Figure 11. The terrain difference between zone A and B are shown as Figure 14. The slope in zone B is slightly greater than that in zone A (the average slope in zone A is 6.2°, the average slope in zone B is 8.5°). The slopes in zone A are generally facing west, and the slopes in zone B are generally facing south or north. The best flight direction for each zone was calculated with Equations (17) and (18), shown as the red circles in Figure 15. The   The geometric distortion in the image changes with the observation zenith. Hence, in order to calculate the radiance difference ∆L λ−inter (φ F , λ) between different view zenith with Equation (16), the image observed from two view zeniths must be geometrically corrected with reference to the nadir image. In order to validate the flight direction for different scenes, the two zones whose terrain was rugged and different from each other were defined corresponding to the A and B in Figure 11. The terrain difference between zone A and B are shown as Figure 14. The slope in zone B is slightly greater than that in zone A (the average slope in zone A is 6.2 • , the average slope in zone B is 8.5 • ). The slopes in zone A are generally facing west, and the slopes in zone B are generally facing south or north. The best flight direction for each zone was calculated with Equations (17) and (18), shown as the red circles in Figure 15. The ∆L λ−inter (φ F , λ) comes to the minimum in the best directions.

Discussion
The digital and physical simulation experiments showed that

Discussion
The digital and physical simulation experiments showed that in the digital simulation experiment, and were reduced from 0.45 to 0.02 in the physical simulation experiment. Although this radiance non-uniformity can also be corrected by BRDF correction algorithms, the accurate algorithms for rugged terrain scenes are still subjects of ongoing research Figure 15. Radiance difference in the overlap area between adjacent stripes: (a) the radiance difference in zone A (b) the radiance difference in zone B.

Discussion
The digital and physical simulation experiments showed that ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) for whickbroom and pushbroom airborne imaging are affected by rugged terrain and could be reduced by adjusting the flight direction. In the case of solo slope, the best flight direction dependent on the illumination direction can be calculated with Equation (19), where ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) come to the minimum. As to the composite slope scene which is composed of many facets, ∆L λ−intra (φ F , λ) and ∆L λ−inter (φ F , λ) of each facet are different. In order to minimize the effect of reflectance anisotropy and terrain in the whole image, the best flight direction is defined as the expectation of the best flight direction for each facet, which would be calculated with Equation (20). By design the flight direction with Equation (20), ∆L λ−inter (φ F , λ) were reduced from 0.40 to 0.04 in the digital simulation experiment, and were reduced from 0.45 to 0.02 in the physical simulation experiment. Although this radiance non-uniformity can also be corrected by BRDF correction algorithms, the accurate algorithms for rugged terrain scenes are still subjects of ongoing research [22]. It would be effective to reduce the radiance non-uniformity during data acquiring and improve the accuracy of correction by acquiring high performance image at first. A flight direction design method for whickbroom and pushbroom airborne imaging is summarized as in Figure 16. At the beginning of flight planning, the AOI is defined according to the goal of experiment. The solar zenith and azimuth angles are estimated with a software such as SolarBeam [39] and SolarEarthTools [40], according to the flight date and time planed. The flight time window is then designed to avoid the shadow effect and to guarantee sufficient solar irradiance onto the surface, according to the solar zenith and azimuth angles with a flight design software such as Leica Mission Pro [9]. The central time in the time window is chosen to design the best flight direction with Equation (20), accounting for the slope and aspect distribution in the AOI. Before conducting the flight mission, the actual time window is defined considering the weather of flight date. Then the actual imaging area is calculated with the flight design software according to the flight direction designed and time window [7][8][9]11,23]. Then the slope and aspect distribution in the actual imaging area, and the solar zenith and azimuth angles during the imaging process, can be recalculated to adjust the flight direction. Following this workflow, the best flight direction for the actual flight would be determined.  [22]. It would be effective to reduce the radiance non-uniformity during data acquiring and improve the accuracy of correction by acquiring high performance image at first. A flight direction design method for whickbroom and pushbroom airborne imaging is summarized as in Figure 15. At the beginning of flight planning, the AOI is defined according to the goal of experiment. The solar zenith and azimuth angles are estimated with a software such as SolarBeam [39] and SolarEarthTools [40], according to the flight date and time planed. The flight time window is then designed to avoid the shadow effect and to guarantee sufficient solar irradiance onto the surface, according to the solar zenith and azimuth angles with a flight design software such as Leica Mission Pro [9]. The central time in the time window is chosen to design the best flight direction with Equation (20), accounting for the slope and aspect distribution in the AOI. Before conducting the flight mission, the actual time window is defined considering the weather of flight date. Then the actual imaging area is calculated with the flight design software according to the flight direction designed and time window [7][8][9]11,23]. Then the slope and aspect distribution in the actual imaging area, and the solar zenith and azimuth angles during the imaging process, can be recalculated to adjust the flight direction. Following this workflow, the best flight direction for the actual flight would be determined.

Conclusions
The anisotropic reflectance of the surface significantly affects the performance of airborne hyperspectral remote sensing image. Current flight direction design methods cannot reduce the effects of rugged terrain as they do for flat terrain. Hence, this paper analyzes the effect of the anisotropic reflectance of the surface to the radiance reached the whickbroom and pushbroom airborne imagers. The radiance gradient in the image and the radiance difference in the overlapping area between adjacent stripes were analyzed for different flight directions and illumination directions. To minimize the radiance gradient and the radiance difference, a method is developed for the composite slope scene illuminated from different directions. Several digital and physical simulation experiments are conducted to demonstrate that the radiance gradient and radiance difference are reduced significantly on the designed flight direction. By cooperating the current flight

Conclusions
The anisotropic reflectance of the surface significantly affects the performance of airborne hyperspectral remote sensing image. Current flight direction design methods cannot reduce the effects of rugged terrain as they do for flat terrain. Hence, this paper analyzes the effect of the anisotropic reflectance of the surface to the radiance reached the whickbroom and pushbroom airborne imagers. The radiance gradient in the image and the radiance difference in the overlapping area between adjacent stripes were analyzed for different flight directions and illumination directions. To minimize the radiance gradient and the radiance difference, a method is developed for the composite slope scene illuminated from different directions. Several digital and physical simulation experiments are conducted to demonstrate that the radiance gradient and radiance difference are reduced significantly on the designed flight direction. By cooperating the current flight planning software with the proposed method, a basic workflow of flight direction design is summarized for flight planning. In this workflow, the flight direction is first calculated considering the DEM of AOI and the flight time. Then flight line is planned by the use of current flight planning software with the calculated flight direction. Finally, the flight direction is adjusted before the flight conducted considering the final flight time. With this workflow, the accuracy of airborne imaging can be improved with proper flight planning. As the scene in the physical simulation experiment is flat, the best flight direction in different subzones doesn't change significantly. More physical experiments will be conducted to simulate the hyperspectral imaging over scenes whose terrains are more rugged in order to analyze the influences of anisotropy of the surface on its retrieved reflectance characteristics.