Open AccessArticle

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

Huijie Zhao

Bolun Cui

and

Guorui Jia

School of Instrumentation & Optoelectronic Engineering, Beihang University, Key Laboratory of Precision Opto-Mechatronics Technology, Ministry of Education, 37# Xueyuan Road, Haidian District, Beijing 100191, China

Author to whom correspondence should be addressed.

Sensors 2019, 19(12), 2715; https://doi.org/10.3390/s19122715

Submission received: 16 May 2019 / Revised: 10 June 2019 / Accepted: 12 June 2019 / Published: 17 June 2019

(This article belongs to the Section Remote Sensors)

Download

Browse Figures

Versions Notes

Abstract

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.

Keywords:

hyperspectral imaging; airborne remote sensing; rugged terrain; flight mission planning

1. 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 among pixels in different columns, there are radiance gradients in the cross direction of an image and radiance differences in the overlap between adjacent stripes [12,13]. These radiance non-uniformities are more obvious in the data from whiskbroom and pushbroom imagers, as the FOVs of these imagers are often larger than 30°. The radiance gradient in the image and the radiance difference in the overlap between adjacent stripes are often affected by the flight direction. As mentioned by Pepe [10] and Myers [14], the sun glint and the ‘hotspot’ should be avoided when designing the flight direction. The ‘hotspot’ is the heavily reflecting sun light at ‘antisolar point’ [10]. In order to avoid the sun glint on the water, Myers [14] and Montes [15] mentioned that the flight line should be the solar azimuth or 180° from the solar azimuth. Downey [16], Beisl [17], Schiefer [18], Langhans [19], Kötz [20] and Zou [21] also designed the flight direction aligned with the sun-target-observer plane (also called principal plane) to minimize the radiance gradient in the image. Schiefer planned a flight direction which was aligned with the principal plane to get the reference image for the validation of the correction of the bidirectional reflectance–distribution function (BRDF) effects [18].

The radiance gradient and radiance difference reach a minimum while the flight direction is parallel to the principle plane only for the horizontal surface and the azimuthal symmetry BRDF condition. However, as the ground surface is not always horizontal, the reflected radiance also depends on the illumination direction, the observation direction, the facing direction of the surface and the bidirectional reflectance of the surface [22]. For rugged terrain, designing a flight line which is aligned with the principal plane might not be the best choice. In the meantime, as the AOI covered by many stripes, the illumination direction of the sun changes when the imaging process occurs [23]. So, the flight line cannot align with the principle plane all the time. Hence, the design method of flight direction still requires improvements.

In this paper, we analyze the at-sensor radiance from a rugged terrain to find the relationship between the effect of directional reflectance of the surface and the flight direction. Then, a method of flight direction design is developed based on this relationship. To test the performance of the method, two digital simulation experiments over a solo slope scene and a composite slope scene are conducted, respectively, for different illumination directions and observation directions. The physical simulation experiment over the composite slope scene is also conducted with more real targets. Finally, a basic workflow of flight direction design for planning flight mission is developed.

2. Methods

The radiance gradient in the image and the radiance difference between the overlap areas of adjacent stripes caused by the directional surface reflectance on the rugged terrain is firstly analyzed with radiance transfer process modeling. The illumination-observation geometry is shown in Figure 1. A horizontal Cartesian coordinate system was defined with the target as the origin, X-axis pointing to the north and Z-axis perpendicular to the earth. For a solo slope as the grey plane in Figure 1, the vector normal to terrain plane is defined as

\vec{T} (θ_{T}, ϕ_{T})

, where

θ_{T}

and

ϕ_{T}

are the slope and aspect of the terrain, respectively. The Z’-axis in the local slope coordinate system X’-Y’-Z’ is defined to be aligned with

\vec{T} (θ_{T}, ϕ_{T})

. The X’-axis and Y’-axis are rotated from X-axis and Y-axis according to

θ_{T}

and

ϕ_{T}

. The illumination direction, observation direction and the relevant angles in these two coordinate systems are listed in Table 1.

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).

\vec{I^{'}} (θ_{i}^{'}, ϕ_{i}^{'}) = R_{Z} (ϕ_{T}) R_{Y} (θ_{T}) \vec{I} (θ_{i}, ϕ_{i}),

(1)

\vec{E^{'}} (θ_{e}^{'}, ϕ_{e}^{'}) = R_{Z} (ϕ_{T}) R_{Y} (θ_{T}) \vec{E} (θ_{e}, ϕ_{e}),

(2)

R_{Z} (ϕ_{T}) = [\begin{matrix} \cos ϕ_{T} & \sin ϕ_{T} & 0 \\ - \sin ϕ_{T} & \cos ϕ_{T} & 0 \\ 0 & 0 & 1 \end{matrix}],

(3)

R_{Y} (θ_{T}) = [\begin{matrix} \cos θ_{T} & 0 & - \sin θ_{T} \\ 0 & 1 & 0 \\ \sin θ_{T} & 0 & \cos θ_{T} \end{matrix}],

(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}),

(6)

ϕ_{i}^{'} = \arctan \frac{- \sin θ_{i} \sin (ϕ_{T} - ϕ_{i})}{\cos θ_{T} \sin θ_{i} \cos (ϕ_{T} - ϕ_{i}) - \sin θ_{T} \cos θ_{i}},

(7)

ϕ_{e}^{'} = \arctan \frac{- \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.

L_{λ} (θ_{e}, ϕ_{e}) = τ_{2} (λ) (\begin{array}{l} {E^{'}}_{s λ} \cos θ_{i}^{'} τ_{1} (λ) \frac{ρ (θ_{i}^{'}, θ_{e}^{'}, ϕ^{'}, λ)}{π} + \int_{F} L_{d λ} (σ, φ) \cos σ \frac{ρ (σ, θ_{e}^{'}, φ, λ)}{π} d Ω \\ + \int_{1 - F} L_{b λ} (σ, φ) \cos σ \frac{ρ (σ, θ_{e}^{'}, φ, λ)}{π} d Ω \end{array}) + L_{u λ},

(9)

ϕ^{'} = ϕ_{e}^{'} - ϕ_{i}^{'},

(10)

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_{v o l} (θ_{i}, θ_{v}, ϕ)

and the geometric kernel

K_{g e o} (θ_{i}, θ_{v}, ϕ)

are functions of illumination direction and observation direction.

f_{v o l} (λ)

and

f_{g e o} (λ)

are the factors of kernels.

ρ_{i s o} (λ)

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].

r_{B R D F} = \frac{ρ}{π (s r)},

(11)

ρ_{B R F} (θ_{i}^{'}, θ_{e}^{'}, ϕ^{'}, λ) = ρ_{i s o} (λ) + f_{v o l} (λ) K_{v o l} (θ_{i}^{'}, θ_{e}^{'}, ϕ^{'}) + f_{g e o} (λ) K_{g e o} (θ_{i}^{'}, θ_{e}^{'}, ϕ^{'}),

(12)

{\begin{cases} K_{v o l} = \frac{4}{3 π} \frac{1}{c o s θ_{i}^{'} + c o s θ_{e}^{'}} [(\frac{π}{2} - ζ) c o s ζ + s i n ζ] - \frac{1}{3} \\ K_{g e o} = \frac{1}{π} (t - s i n t c o s t) (\frac{1}{c o s θ_{i}^{'}} + \frac{1}{c o s θ_{e}^{'}}) - (\frac{1}{c o s θ_{i}^{'}} + \frac{1}{c o s θ_{e}^{'}}) + \frac{1 + c o s ζ}{2 c o s θ_{i}^{'} c o s θ_{e}^{'}} \end{cases},

(13)

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

\pm \frac{F O V}{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°.

ϕ_{F} = ϕ_{e} + 90^{\circ},

(14)

Δ L_{λ - intra} (ϕ_{F}, λ) = \frac{\max (| L_{λ} (\frac{F O V}{2}, ϕ_{e}, λ) - L_{λ} (0, ϕ_{e}, λ) |, | L_{λ} (- \frac{F O V}{2}, ϕ_{e}, λ) - L_{λ} (0, ϕ_{e}, λ) |)}{L_{λ} (0, ϕ_{e}, λ)},

(15)

Δ L_{λ - inter} (ϕ_{F}, λ) = \frac{L_{λ} (\frac{F O V}{2}, ϕ_{e}, λ) - L_{λ} (- \frac{F O V}{2}, ϕ_{e}, λ)}{L_{λ} (0, ϕ_{e}, λ)},

(16)

Δ ϕ_{F} = ϕ_{F}^{'} - ϕ_{i}^{'},

(17)

ϕ_{F}^{'} = \arctan \frac{\tan (ϕ_{T} - ϕ_{F})}{\cos θ_{T}},

(18)

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).

ϕ_{F} = ϕ_{T} + \arctan \frac{\cos θ_{T} \sin θ_{i} \sin (ϕ_{T} - ϕ_{i})}{\sin θ_{T} \cos θ_{i} - \cos θ_{T} \sin θ_{i} \cos (ϕ_{T} - ϕ_{i})},

(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).

E (ϕ_{F}) = \int_{ϕ_{F}} ϕ_{F} \cdot P (ϕ_{F}) d ϕ_{F},

(20)

3. 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

Δ L_{λ - inter} (ϕ_{F}, λ)

is similar to that of

Δ L_{λ - intra} (ϕ_{F}, λ)

Δ L_{λ - intra} (ϕ_{F}, λ)

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).

4. 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.

Δ 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.

5. 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.

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.

6. 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.

7. 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.

Author Contributions

H.Z. and B.C. conceived and designed the experiments; B.C. and G.J. analyzed the theory of the experiments; B.C. performed the experiments; B.C. and G.J. analyzed the data; B.C. wrote the paper.

Funding

This work was partly supported by the National Key Research and Development Program of China (Grant No. 2016YFB0500502), the Civil Aerospace Technology Advance Research Project (Grant No. Y7K00100KJ), and the National Natural Science Foundation of China (Grant No. 61675012).

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

References

Bioucas-Dias, J.M.; Plaza, A.; Camps-Valls, G.; Scheunders, P.; Nasrabadi, N.; Chanussot, J. Hyperspectral Remote Sensing Data Analysis and Future Challenges. IEEE Geosci. Remote Sens. Mag. 2013, 1, 6–36. [Google Scholar] [CrossRef]
Foster, J.; Townsend, P. Mapping Forest Composition in the Central Appalachians Using AVIRIS: Effects of Topography and Phenology. Available online: https://www.researchgate.net/profile/Jane_Foster/publication/228386320_Mapping_forest_composition_in_the_central_appalachians_using_AVIRIS_Effects_of_topography_and_phenology/links/0912f50a2cf058d877000000.pdf (accessed on 14 June 2019).
Ong, C.; Cudahy, T. Deriving Quantitative Monitoring Data Related to Acid Drainage Using Multi-Temporal Hyperspectral Data. Available online: https://www.researchgate.net/profile/Cindy_Ong/publication/265264343_DERIVING_QUANTITATIVE_MONITORING_DATA_RELATED_TO_ACID_DRAINAGE_USING_MULTI-_TEMPORAL_HYPERSPECTRAL_DATA/links/548f6d230cf214269f263bad.pdf (accessed on 14 June 2019).
Ben-Dor, E.; Levin, N.; Singer, A.; Karnieli, A.; Braun, O.; Kidron, G. Quantitative mapping of the soil rubification process on sand dunes using an airborne CASI hyperspectral sensor. Geoderma 2006, 131, 1–21. [Google Scholar] [CrossRef]
Ben-Dor, E.; Chabrillat, S.; Demattê, J.A.M.; Taylor, G.R.; Hill, J.; Whiting, M.L.; Sommer, S. Using Imaging Spectroscopy to study soil properties. Remote Sens. Environ. 2009, 113, S38–S55. [Google Scholar] [CrossRef]
Vreys, K.; Iordache, M.D.; Bomans, B.; Meuleman, K. Data acquisition with the APEX hyperspectral sensor. Misc. Geogr. 2016, 20, 5–10. [Google Scholar] [CrossRef] [Green Version]
Cline, J.D.; Wilson, J.A.; Feher, S.H.; Ward, G.D. Airborne flight planning and information system. Canadian Vet. J. Rev. Vétérinaire Can. 1987, 47, 999–1002. [Google Scholar] [CrossRef]
Leica XPro Data Processing at the Speed of Flight. Available online: https://leica-geosystems.com/products/airborne-systems/software/leica-xpro (accessed on 10 January 2019).
Ip, A.W.; Mostafaa, M.M.; Huttona, J.; Barriere, J.P. An optimally integrated Direct georeferencing and flight management system for increased productivity of airborne mapping and remote sensing. Remote Sens. Spat. Inf. Sci. Int. Arch. Photogramm. 2008, XXXVII, 579–584. [Google Scholar]
Pepe, M.; Fregonese, L.; Scaioni, M. Planning airborne photogrammetry and remote sensing missions with modern platforms and sensors. Eur. J. Remote Sens. 2018, 51, 412–436. [Google Scholar] [CrossRef]
Möllney, M.; Kremer, J. Contour Flying for Airborne Data Acquisition. Photogramm. Week 2013, 13, 117–129. [Google Scholar]
Collings, S.; Caccetta, P.; Campbell, N.; Wu, X. Techniques for BRDF Correction of Hyperspectral Mosaics. IEEE Trans. Geosci. Remote Sens. 2010, 48, 3733–3746. [Google Scholar] [CrossRef]
Feingersh, T.; Ben-Dor, E.; Filin, S. Correction of reflectance anisotropy: A multi-sensor approach. Int. J. Remote Sens. 2010, 31, 49–74. [Google Scholar] [CrossRef]
Myers, J.S.; Miller, R.L. Optical airborne remote sensing. In Remote Sensing of Coastal Aquatic Environments; Remote Sensing and Digital Image Processing; Miller, R.L., Del Castillo, C.E., Mckee, B.A., Eds.; Springer: Dordrecht, The Netherlands, 2005; Volume 7, pp. 51–67. [Google Scholar] [CrossRef]
Montes, M.J.; Gao, B.-C.; Davis, C.O.; Moline, M. Analysis of AVIRIS Data From LEO-15 Using Tafkaa Atmospheric Correction. In Proceedings of the 12th AVIRIS/HYPERION Earth Science Workshop, Pasadena, CA, USA, 1 January 2003. [Google Scholar]
Downey, M.; Uebbing, R.; Gehrke, S.; Beisl, U. Radiometric processing of ads imagery: Using atmospheric and BRDF corrections in production. In Proceedings of the ASPRS 2010 Annual Conference, San Diego, CA, USA, 26–30 April 2010. [Google Scholar]
Beisl, U.; Adiguezel, M. Validation of the reflectance calibration of the ADS40 airborne sensor using ground reflectance measurements. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci. ISPRS Arch. 2010, 38, 80–85. [Google Scholar]
Schiefer, S.; Hostert, P.; Damm, A. Correcting brightness gradients in hyperspectral data from urban areas. Remote Sens. Environ. 2006, 101, 25–37. [Google Scholar] [CrossRef]
Langhans, M.; Van der Linden, S.; Damm, A.; Hostert, P. The influence of bidirectional reflectance in airborne hyperspectral data on spectral angle mapping and linear spectral mixture analysis. In Proceedings of the 5th EARSeL Workshop on Imaging Spectroscopy, Bruges, Belgium, 23–25 April 2007. [Google Scholar]
Kötz, B.; Schaepman, M.; Morsdorf, F.; Bowyer, P.; Itten, K.; Allgöwer, B. Radiative transfer modeling within a heterogeneous canopy for estimation of forest fire fuel properties. Remote Sens. Environ. 2004, 92, 332–344. [Google Scholar] [CrossRef]
Zou, X.; Hernández-Clemente, R.; Tammeorg, P.; Lizarazo Torres, C.; Stoddard, F.L.; Mäkelä, P.; Pellikka, P.; Mõttus, M. Retrieval of leaf chlorophyll content in field crops using narrow-band indices: Effects of leaf area index and leaf mean tilt angle. Int. J. Remote Sens. 2015, 36, 6031–6055. [Google Scholar] [CrossRef]
Wen, J.; Liu, Q.; Xiao, Q.; Liu, Q.; You, D.; Hao, D.; Wu, S.; Lin, X. Characterizing Land Surface Anisotropic Reflectance over Rugged Terrain: A Review of Concepts and Recent Developments. Remote Sens. 2018, 10, 370. [Google Scholar] [CrossRef]
Product Description Leica MissionPro. Available online: https://leica-geosystems.com/products/airborne-systems/software/leica-missionpro (accessed on 1 March 2019).
Schott, J.R. Remote Sensing: The Image Chain Approach, 2nd ed.; Oxford University Press: New York, NY, USA, 2007. [Google Scholar]
Hapke, B. Theory of Reflectance and Emittance Spectroscopy; Cambridge University Press: New York, NY, USA, 2005. [Google Scholar]
Jacquemoud, S.; Verhoef, W.; Baret, F.; Bacour, C.; Zarco-Tejada, P.J.; Asner, G.P.; François, C.; Ustin, S.L. PROSPECT+SAIL models: A review of use for vegetation characterization. Remote Sens. Environ. 2009, 113, S56–S66. [Google Scholar] [CrossRef]
Schaepman-Strub, G.; Schaepman, M.E.; Painter, T.; Dangel, S.; Martonchik, J.V. Reflectance quantities in optical remote sensing-definitions and case studies. Remote Sens. Environ. 2006, 103, 27–42. [Google Scholar] [CrossRef]
Bachmann, C.M.; Eon, R.S.; Ambeau, B.; Harms, J.; Badura, G.; Griffo, C. Modeling and intercomparison of field and laboratory hyperspectral goniometer measurements with G-LiHT imagery of the Algodones Dunes. J. Appl. Remote Sens. 2017, 12, 012005. [Google Scholar] [CrossRef]
Lunagaria, M.M.; Patel, H.R. Evaluation of PROSAIL inversion for retrieval of chlorophyll, leaf dry matter, leaf angle, and leaf area index of wheat using spectrodirectional measurements. Int. J. Remote Sens. 2018, 10. [Google Scholar] [CrossRef]
Koukal, T.; Atzberger, C.; Schneider, W. Evaluation of semi-empirical BRDF models inverted against multi-angle data from a digital airborne frame camera for enhancing forest type classification. Remote Sens. Environ. 2014, 151, 27–43. [Google Scholar] [CrossRef]
Eon, R.; Bachmann, C.M.; Gerace, A. Retrieval of Sediment Fill Factor by Inversion of a Modified Hapke Model Applied to Sampled HCRF from Airborne and Satellite Imagery. Remote Sens. 2018, 10, 1758. [Google Scholar] [CrossRef]
Botev, Z.I.; Grotowski, J.F.; Kroese, D.P. Kernel density estimation via diffusion. Ann. Stat. 2010, 38, 2916–2957. [Google Scholar] [CrossRef] [Green Version]
Berk, A.; Anderson, G.P.; Acharya, P.K.; Chetwynd, J.H.; Bernstein, L.S.; Shettle, E.P.; Matthew, M.W.; Adler-Golden, S.M. MODTRAN4 User’s Manual. Air Force Research Laboratory; Space Vehicles Directorate, Air Force Materiel Command: Hanscom AFB, MA, USA, 2000. [Google Scholar]
Hatchell, D.C. Analytical Spectral Device Technical Guide; Analytical Spectral Devices, Inc.: Boulder, CO, USA, 1999. [Google Scholar]
Zhao, J.; Wang, W.; Cheng, Q.; Agterberg, F. Mapping of Fe mineral potential by spatially weighted principal component analysis in the eastern Tianshan mineral district, China. J. Geochem. Explor. 2016, 164, 107–121. [Google Scholar] [CrossRef]
Deng, Y.F.; Song, X.Y.; Hollings, P.; Chen, L.M.; Zhou, T.; Yuan, F.; Xie, W.; Zhang, D.; Zhao, B. Lithological and geochemical constraints on the magma conduit systems of the Huangshan Ni-Cu sulfide deposit, NW China. Miner. Depos. 2017, 52, 845–862. [Google Scholar] [CrossRef]
Zhao, H.; Cui, B.; Jia, G.; Li, X.; Zhang, C.; Zhang, X. A “Skylight” Simulator for HWIL Simulation of Hyperspectral Remote Sensing. Sensors 2017, 17, 2829. [Google Scholar] [CrossRef]
Xia, M.; Jiang, C.Y.; Qian, Z.Z.; Xia, Z.D.; Wang, B.Y.; Sun, T. Geochemistry and petrogenesis of Huangshandong intrusion, East Tianshan, Xinjiang. Acta Petrol. Sin. 2010, 26, 2413–2430. [Google Scholar] [CrossRef]
Matusiak, B. The Graphical Tool for Sky Component, Solar Glare and Overheating Risk Prediction. Available online: https://www.researchgate.net/profile/Barbara_Matusiak/publication/268979885_THE_GRAPHICAL_TOOL_FOR_SKY_COMPONENT_SOLAR_GLARE_AND_OVERHEATING_RISK_PREDICTION/links/547c7eba0cf2cfe203bfc832/THE-GRAPHICAL-TOOL-FOR-SKY-COMPONENT-SOLAR-GLARE-AND-OVERHEATING-RISK-PREDICTION.pdf (accessed on 14 June 2019).
Sun Earth Tools-Sun Position. Available online: http://www.sunearthtools.com (accessed on 5 March 2019).

Figure 1. The geometric model for illumination and observation [24].

Figure 2. The radiative transfer process [24].

Figure 3. The BRF of sand* [31]. *The azimuth angle changing on circumferential direction align with the blue circle; the zenith angle changing on the radius direction; the blue arrow denotes the flight direction; for specific flight direction, the view directions of two stripes on the same target are the red circles.

Figure 4. 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.

Figure 5. Radiance gradient for different illumination directions and flight directions: (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.

Figure 6. The BRF of samples collected with illumination zenith 45° and view zenith 0°.

Figure 7. The simulated radiance of the overlap area, (a) the image acquired from 60° flight direction and −30° zenith; (b) the image acquired on 60° flight direction and +30° zenith.

Figure 8. The PDFs of slope and aspect in zone A and B: (a) the PDFs of slope, (b) the PDFs of aspect.

Figure 9. Radiance difference in the overlap area between adjacent stripes: (a) the radiance difference in zone A (b) the radiance difference in zone B.

Figure 10. Flowchart of physical simulation experiment.

Figure 11. Photo of the physically simulated scene.

Figure 12. BRF of smashed rock samples of different diameters.

Figure 13. The reflected radiance of scene, (a) on flight direction 90° and view zenith −30°; (b) on flight direction 90° and view zenith −30°.

Figure 14. PDFs of slope and aspect in zone A and B: (a) the PDFs of slope, (b) the PDFs of aspect.

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.

Figure 16. The workflow of the optimized flight direction design method.

Table 1. Notations used for the description of radiance transfer process.

Symbols	Explanations
$\vec{T} (θ_{T}, ϕ_{T})$	Normal vector of terrain in horizontal coordinate system
$\vec{I} (θ_{i}, ϕ_{i})$	Illumination vector (from sun to target) in horizontal coordinate system
$\vec{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

Table 2. Simulation Parameters for solo slope digital simulation experiment.

Attribute	Value
Solar Zenith Angle ( $θ_{i}$ )	20°
Solar Azimuth Angle ( $ϕ_{i}$ )	94°–266°
Observation Zenith Angle ( $θ_{e}$ )	−30°/0°/+30°
Flight direction ( $ϕ_{F}$ )	90°–270° (5° interval)
Slope	0°–50° (10° interval)
Aspect	0°/180°
Spectral Range	400 nm–2500 nm
Spectral Resolution	10 nm
Atmosphere Model	Mid-Latitude Summer
Aerosol Model	Rural
Visibility	40 km

Table 3. Simulation Parameters for composite slope digital simulation experiment.

Attribute	Value
Location	40.0°N, 93.9°E
Solar Direction	Zenith: 17.4°~46.9°, Azimuth: 94.3°~265.8°
Spectral Range	400 nm~2500 nm
Spectral Resolution	10 nm
Atmosphere Model	Mid-Latitude Summer
Aerosol Model	Rural
Visibility	40 km
View Zenith	−30.0°/0°/+30.0°
Flight direction	90.0°~−90.0° (30.0° interval)
Slope	0°~36° (mean: 23°)
Aspect	−180°~180° (mean: 1.4°)

Table 4. Simulation Parameters for composite slope physical simulation experiment.

Attribute	Value
Location	40.0°N, 93.9°E
Solar Direction	Zenith: 48°, Azimuth: 90°~270° (30.0° interval)
Atmosphere Model	Mid-Latitude Summer
Aerosol Model	Rural
Visibility	40 km
View Zenith	−30.0°/0°/+30.0°
Flight direction	−90.0°~+90.0°
Slope	0°~50.2° (mean: 8.6°)
Aspect	−180°~180° (mean: 7.5°)

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Zhao, H.; Cui, B.; Jia, G. A Flight Direction Design Method for Airborne Spectral Imaging Considering the Anisotropy Reflectance of the Target in Rugged Terrain. Sensors 2019, 19, 2715. https://doi.org/10.3390/s19122715

AMA Style

Zhao H, Cui B, Jia G. A Flight Direction Design Method for Airborne Spectral Imaging Considering the Anisotropy Reflectance of the Target in Rugged Terrain. Sensors. 2019; 19(12):2715. https://doi.org/10.3390/s19122715

Chicago/Turabian Style

Zhao, Huijie, Bolun Cui, and Guorui Jia. 2019. "A Flight Direction Design Method for Airborne Spectral Imaging Considering the Anisotropy Reflectance of the Target in Rugged Terrain" Sensors 19, no. 12: 2715. https://doi.org/10.3390/s19122715

APA Style

Zhao, H., Cui, B., & Jia, G. (2019). A Flight Direction Design Method for Airborne Spectral Imaging Considering the Anisotropy Reflectance of the Target in Rugged Terrain. Sensors, 19(12), 2715. https://doi.org/10.3390/s19122715

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu