What Can Multifractal Analysis Tell Us about Hyperspectral Imagery?

Hyperspectral images provide complex information about the Earth’s surface due to their very high spectral resolution (hundreds of spectral bands per pixel). Effective processing of such a large amount of data requires dedicated analysis methods. Therefore, this research applies, for the first time, the degree of multifractality to the global description of all spectral bands of Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) data. Subsets of four hyperspectral images, presenting four landscape types, are analysed. In particular, we verify whether multifractality can be detected in all spectral bands. Furthermore, we analyse variability in multifractality as a function of wavelength, for data before and after atmospheric correction. We try to identify absorption bands and discuss whether multifractal parameters provide additional value or can help in the problem of dimensionality reduction in hyperspectral data or landscape type classification.


Fractals
Mandelbrot ( [1,2]) first drew the attention of scientists to fractals, and they have since proven to be helpful in developing a numerical description of irregularities occurring in nature. He argued that fractals can be characterised by a quantitative parameter-the fractal dimension (D F ). The latter corresponds to measures used in Euclidean geometry (e.g., length or area), but in a scale of measurement function. Over time, applications have expanded from a simple description of linear features to surface feature analysis and research on information stored in the form of images [3]. Remote sensing researchers have applied the concept to aerial and satellite imagery analysis, including radar data processing for irrigation mapping [4], surface modelling [5], land cover mapping [6], estimation of building density [7] and segmentation [8], among others.
The most popular D F calculation methods used in remote sensing applications include the blanket method [9], the triangular prism [10], the isarithm method [11], Sevcik's method [12], the power spectrum method [13], the modified variogram method [14] and the adapted Hausdorff metric [15]. Besides the number of algorithms for D F calculation, there are a few applications where this parameter is associated with remote sensing data. For example, in a single spectral band, it can be calculated for each pixel (local description)-or for a square subset/patch of the image (global description) [16]. In some cases, the subset covers the whole image. 2 of 21 In the context of hyperspectral data, D F has been used in two ways. The first approach concerns an analysis of spectral signatures/profiles that refer to single pixels (e.g., [15,17,18]). Usually, the purpose of these analyses is to reduce the dimensionality of hyperspectral data. For example, Mukherjee et al. [19] showed that D F can be almost as accurate as conventional methods, but with lower computational requirements.
In a second approach, D F is calculated for spectral bands in the whole image or a subset ( [20][21][22]). Qiu et al. [20] compared a few calculation methods for each spectral band in two aerial images acquired by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) sensor. Krupiński et al. [22] analysed the same images with the Differential Box Counting (DBC) method [23] and additionally analysed four types of landscape. The latter work found significant variability in D F determined for radiance and reflectance data as a function of wavelength. In addition, in noisy bands, disturbance in D F values was observed as a result of strong radiation absorption by water molecules. The authors suggested that this parameter could be used, inter alia, to determine the quality of spectral bands or separability of landscape types.

Multifractals
A number of experiments have shown, however, that a fractal formalism may not be sufficient for satellite imagery description, and that multifractals (an extension of fractal theory) should be used instead (e.g., [24][25][26]). The multifractal approach is based on the assumption that it is necessary to use a number of non-trivially connected fractals, each with a different dimension of self-similarity, to describe data complexity. What is more, there are various functions (e.g., generalised dimensions, multifractal spectra [27,28]) and a number of quantitative parameters associated with multifractal description. Together, these open up a wide range of possible applications of multifractals in the Earth Observation domain (see Table 1). A review of existing implementations of a multifractal formalism in the analysis of remote sensing images indicates that they are mostly related to the segmentation and classification of radar imagery. Other applications include, for example, the detection of changes ( [16,31]), the edge-preserving smoothing of high-resolution satellite images [32], super-resolution image analysis [33], compression of remote sensing imagery [34] and usage as textural features for content-based image retrieval [35]. There are also combinations with methods such as granulometric analysis [36] or principal component analysis [37] (see Table 2). Our previous research comprehensively analysed the multifractal character of panchromatic very high resolution (VHR) satellite images acquired by the following satellites: EROS [26], WorldView-2 [25], GeoEye-1 and Pléiades 1A [45,46], using the box-counting based moment method. In particular, we proposed a quantitative parameter, named the degree of multifractality (∆), and demonstrated its usefulness in attempts to distinguish the basic forms of land cover (water, urban areas, forests and agriculture). In general, the method was effective in automatically assigning image subsets to certain classes. We also showed that the degree of multifractality can be a very effective descriptor of image content, compared to other textural features [47]. However, it should be noted that our global analysis was exclusively focused on a panchromatic band of satellite imagery. The exception was 30-m pixels from Landsat satellites, where we initially considered six spectral bands [48].
In the context of hyperspectral data analyses, we found only four examples in the literature where a multifractal formalism was applied. More precisely, Combrexelle et al. [30] calculated coefficients of the polynomial describing the multifractal spectrum for each spectral band. Li et al. in [42] showed that using multifractal parameters for spectral profiles' description may improve the average accuracy by 7-8%. Incorporation of multifractal parameters resulted in higher overall accuracy (almost 10%) in the algorithm proposed by Wan et al. [43]. Krupiński et al. showed [44] that overall classification accuracies of methods that use multifractal parameters to describe the spectral curve may differ by 7-9%, depending on the pre-processing method. A summary of the current state of knowledge regarding the usefulness of (multi)fractal formalism in the analysis of hyperspectral data is presented in Table 1 (global and local descriptions of spectral bands) and Table 2 (spectral curve description). It is worth noting that in all of the reviewed papers, the analysed datasets were limited to a small number of images or subsets (from 1 to 6), and mostly applied to airborne imagery. Our summary indicates that the multifractal analysis of hyperspectral imagery remains very limited, and we aim to fill this gap.
Therefore, the aim of this research is to consider, for the first time, degree of multifractality as a new global description of all spectral bands for AVIRIS data, for four landscape types. In particular, we verify whether multifractal character is present in all spectral bands. Additionally, we analyse the variability of multifractal features in the function of wavelength, both before and after atmospheric correction. Following this analysis, we try to discuss the areas of potential applicability for the presented methodology. We consider identification of absorption bands and whether multifractal parameters can help to overcome the problem of dimensionality reduction in hyperspectral data or with landscape type classification. Finally, a short comparison of determined multifractal features with statistical and fractal descriptors is performed.
This paper is organised as follows. In Section 2, we present our dataset. The methodology related to the concept of the multifractal is described in Section 3. The results of our analysis are presented and discussed in Section 4, the discussion is in Section 5 and the main conclusions are outlined in Section 6.

Data
The Jet Propulsion Laboratory (JPL), managed by NASA, designed and developed the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS). It was the first imaging spectrometer able to measure the solar spectrum from 400 to 2500 nm as 224 contiguous channels (spectral bands) [49].
Data acquired by AVIRIS over California (USA) are used in our analyses. More precisely, we considered four subsets of four images. Each subset measured 512 × 512 pixels. Each subset represented a different type of landscape, i.e., agriculture, mountains, urban and water (see Figure 1). Analyses were performed on two types of data-radiance (224 bands) and reflectance (197 bands). The reflectance dataset did not include bands with strong radiation absorption caused by water particles in the air, which are removed by the JPL during the atmospheric correction process. As noted in [49], the Atmosphere Removal Algorithm (ARTEM), developed by Gao et al. [50] and updated by Gao and Davis [51], is applied to generate Level-2 products. The spatial resolution of analysed images is about 15-16 m. A list of images used for sample selection is presented in Table A1 (Appendix A). The size of the subsets is given by the calculation and requires square images, where the number of rows (and columns) is a power of 2. The same dataset, among others, was previously analysed with the D F parameter by Krupiński et al. in [22]. According to the reviews presented in Tables 1 and 2, AVIRIS is the most commonly used sensor for fractal-related research. moments (mean, variance, skewness and kurtosis) for radiance (left column) and reflectance (right column) samples.

Methodology
The degree of multifractality (Δ) is a quantitative multifractal parameter that describes the nonhomogeneity of considered data and several methods can be used in its estimation (e.g., [25,44]). In this study, to calculate the degree of multifractality (Δ) for each spectral band separately, for both types of imagery (radiance and reflectance), we applied the box-counting based moment method [27] as described in, for example, [25]. More precisely, we considered the variability of the spectrum of the generalised dimension as a non-increasing function of the real index , −∞ < < +∞ (see, for example, [52]). In the first step, an image (in a given spectral band) of size × was divided into ( ) = ( / ) 2 square boxes of size × (starting with a single pixel of size = 1 and ending with a box of size = m); see [25] and Figure 1. For a given box, the normalised measure was calculated according to the formula: (1) where = 1, . . . , ( ) are labels given to individual boxes sized , and ( ) = ∑ ( , ) is the greyscale intensity at point ( , ), and is the set of all pixels in box . It is worth noting that in the context of multifractal analysis, different measures can be constructed that emphasise various effects and describe different physical processes in the considered data ( [47,53]). To facilitate the results' description, we divided the spectrum into six zones (VIS, NIR, A 1, SWIR 1, A 2 and SWIR 2). The bands of zones A 1 and A 2 were removed from reflectance as these bands were strongly affected by radiation absorption caused by water particles in the air (noisy bands) [49]. Exact ranges (bands and wavelengths) of the six zones are summarised in Table 3. More details about the analysed subsets may be found in Figure A1 (Appendix A), which presents four statistical moments (mean, variance, skewness and kurtosis) for radiance (left column) and reflectance (right column) samples.

Methodology
The degree of multifractality (∆) is a quantitative multifractal parameter that describes the non-homogeneity of considered data and several methods can be used in its estimation (e.g., [25,44]). In this study, to calculate the degree of multifractality (∆) for each spectral band separately, for both types of imagery (radiance and reflectance), we applied the box-counting based moment method [27] as described in, for example, [25]. More precisely, we considered the variability of the spectrum of the generalised dimension D q as a non-increasing function of the real index q, −∞ < q < +∞ (see, for example, [52]). In the first step, an image (in a given spectral band) of size m × m was divided into N (δ) = (m/δ) 2 square boxes of size δ × δ (starting with a single pixel of size δ = 1 and ending with a box of size δ = m); see [25] and Figure 1. For a given box, the normalised measure was calculated according to the formula: where i = 1, . . . , N (δ) are labels given to individual boxes sized δ, and p i (δ) = (k,l)∈Ω i g(k, l), where g(k, l) is the greyscale intensity at point (k, l), and Ω i is the set of all pixels in box δ. It is worth noting that in the context of multifractal analysis, different measures can be constructed that emphasise various effects and describe different physical processes in the considered data ( [47,53]). In the next step, the partition function χ(q, δ) for various values of δ and q was computed according to the formula: As q varies in (2), different subsets associated with different densities dominate. For multifractal measures, the partition function χ(q, δ) scales with the box size δ → 0 as: Based on (2) and (3), we obtained the spectrum of generalised dimensions D q defined, for the first time, by [28]: Additionally, the uncertainties of determining the corresponding slopes log (4)), obtained using unweighted least squares fitting, informed us about the accuracy of the calculation of the D q .
Next, the variability of the multifractal function D q , given as the difference between the maximum (D −∞ ) and minimum (D +∞ ) dimensions, defines the degree of multifractality ∆ [52,53] as follows: Theoretically, this generalised dimension function D q is defined for all real values of q [53]. In practice, the limited dataset means that we can only determine values of D q for a narrow number of moments q (see, for example, [45,47,52]). In this study, we performed an initial examination and verified what kind of q range each sample yields. Then, an optimum range −3 ≤ q ≤ 8 that overlaps the individual q ranges was chosen and applied in further steps. Finally, parameter ∆ was calculated as the difference between D −3 and D +8 , while the sum of errors of D −3 and D +8 gave the error obtained for each value of ∆.

Results
The multifractal analysis described in Section 3 was performed for each spectral band for four examples of different landscapes and, in each case, for radiance and reflectance data. In all considered cases, the scaling given by Equation (3) was revealed and the multifractal nature of the analysed data confirmed. As a consequence, the function (4) and final degree of multifractality (∆) could be determined and used as a quantitative data descriptor. Figure 2 presents the change in values of degree of multifractality ∆ (Y-axis) with wavelength/ band number (X-axis) for the four landscape types: water (blue), mountains (black), agriculture (green) and urban (red). The information about error estimation is also given in the form of vertical bars (∆ ± error). A general overview of the results shows that the values of ∆ determined for radiance data samples are within the range of 0 and 0.768. It is worth noting that both lower and upper extremes relate to water landscape sample (bands 4 and 111, respectively).

Figure 2.
Degree of multifractality (Δ) determined for radiance images of water (blue), mountain (black), agriculture (green) and urban (red) landscapes. The whole spectrum is divided according to ranges presented in Table 3.
In the case of agriculture, the minimum value of Δ is 0.006 and the maximum is 0.498. For mountains, Δ ranges from 0.001 to 0.581, and from 0.002 to 0.484 for urban areas. Water landscape has the lowest Δ values in all spectral bands, except the A 1 and A 2 zones. Figure 2 shows that for all landscape types, Δ is lowest within the VIS zone. With the exception of agriculture, the highest Δ value is found in the A 1 or A 2 zones. Besides the A 1 and A 2 zones, Δ values are the highest and the most variable in the SWIR 2 zone, for all four types of landscape. More detailed information about Δ values (minimum, mean and maximum) is presented in Table A2.

Degree of Multifractality for Reflectance
In order to study the multifractal character of reflectance data, we prepared a comparison between Δ values computed for radiance and reflectance. Figure 3 presents these values for each zone of the electromagnetic spectrum listed in Table 3 (except A 1 and A 2). Moreover, we calculated the difference |Δ − Δ | and it is presented in Table 4. The first thing we notice, analysing the difference, is the huge increment (exceeding 7000) in Δ values calculated for the first few bands of the mountain landscape subset. A closer examination of these bands revealed that they contain between 1.5% and 99.7% of pixels with a value of 0 (this issue is not observed in radiance data). Moreover, checks of the original source image indicated that most mountain areas are covered by 0-value pixels in the first few bands of reflectance data. As analogous situations were noted in other types of landscape, we decided to analyse only bands without this error. Consequently, all bands where the increment in number of pixels with a value of 0 (between radiance and reflectance) was higher than 1% (~2600 pixels) were excluded from further comparison. Bands 1-3 and 224 were removed from the agriculture reflectance sample; 1-7 from the mountain Figure 2. Degree of multifractality (∆) determined for radiance images of water (blue), mountain (black), agriculture (green) and urban (red) landscapes. The whole spectrum is divided according to ranges presented in Table 3.
In the case of agriculture, the minimum value of ∆ is 0.006 and the maximum is 0.498. For mountains, ∆ ranges from 0.001 to 0.581, and from 0.002 to 0.484 for urban areas. Water landscape has the lowest ∆ values in all spectral bands, except the A 1 and A 2 zones. Figure 2 shows that for all landscape types, ∆ is lowest within the VIS zone. With the exception of agriculture, the highest ∆ value is found in the A 1 or A 2 zones. Besides the A 1 and A 2 zones, ∆ values are the highest and the most variable in the SWIR 2 zone, for all four types of landscape. More detailed information about ∆ values (minimum, mean and maximum) is presented in Table A2.

Degree of Multifractality for Reflectance
In order to study the multifractal character of reflectance data, we prepared a comparison between ∆ values computed for radiance and reflectance. Figure 3 presents these values for each zone of the electromagnetic spectrum listed in Table 3 (except A 1 and A 2). Moreover, we calculated the difference ∆ re f lectance − ∆ radiance and it is presented in Table 4.
sample; 1-2 and 223-224 from the urban sample; and 223-224 from the water sample. Table 4 summarises the statistics after the removal of bands with this error. Figure 3. Degree of multifractality calculated for radiance and reflectance images of water (blue), mountain (black), agriculture (green) and urban (red) landscapes, divided into four zones according to Table 3.
Detailed analysis of Figure 3 and Table 4 reveals that the first zone, VIS, is the part of the electromagnetic spectrum where the biggest differences between radiance and reflectance results appear. More precisely, in the whole VIS zone, Δ reflectance values are higher than Δ radiance values for all four landscape types (almost four times in average). It means that the influence of atmospheric corrections is reflected in Δ values.
In the NIR zone, the average absolute change in Δ values is 9% (the highest in water-25%-and the lowest in the agriculture sample-4%). In the SWIR 1 zone, the average absolute change in Δ is 7% (the highest is for water (11%); the lowest is for agriculture (4%)). The biggest changes in Δ occur in bands 115-120 and 151-152 for mountains, observed as a smoothing of curve shape in Figure 3. These bands represent regions of the spectrum that are adjacent to zones of water absorption (A 1 and A 2).
In the SWIR 2 zone, there is the lowest average change in Δ in the whole analysed spectrum (3%), so influence of atmospheric correction on these bands is very small. Like SWIR 1, the highest relative change is found for water (6%) and the lowest for agriculture (1%). The most significant changes are also observed on the edges of the SWIR 2 zone for mountains (bands 169-170 and 221-224). . Degree of multifractality calculated for radiance and reflectance images of water (blue), mountain (black), agriculture (green) and urban (red) landscapes, divided into four zones according to Table 3. The first thing we notice, analysing the difference, is the huge increment (exceeding 7000) in ∆ values calculated for the first few bands of the mountain landscape subset. A closer examination of these bands revealed that they contain between 1.5% and 99.7% of pixels with a value of 0 (this issue is not observed in radiance data). Moreover, checks of the original source image indicated that most mountain areas are covered by 0-value pixels in the first few bands of reflectance data. As analogous situations were noted in other types of landscape, we decided to analyse only bands without this error. Consequently, all bands where the increment in number of pixels with a value of 0 (between radiance and reflectance) was higher than 1% (~2600 pixels) were excluded from further comparison. Bands 1-3 and 224 were removed from the agriculture reflectance sample; 1-7 from the mountain sample; 1-2 and 223-224 from the urban sample; and 223-224 from the water sample. Table 4 summarises the statistics after the removal of bands with this error.
Detailed analysis of Figure 3 and Table 4 reveals that the first zone, VIS, is the part of the electromagnetic spectrum where the biggest differences between radiance and reflectance results appear. More precisely, in the whole VIS zone, ∆ reflectance values are higher than ∆ radiance values for all four landscape types (almost four times in average). It means that the influence of atmospheric corrections is reflected in ∆ values.
In the NIR zone, the average absolute change in ∆ values is 9% (the highest in water-25%-and the lowest in the agriculture sample-4%). In the SWIR 1 zone, the average absolute change in ∆ is 7% (the highest is for water (11%); the lowest is for agriculture (4%)). The biggest changes in ∆ occur in bands 115-120 and 151-152 for mountains, observed as a smoothing of curve shape in Figure 3. These bands represent regions of the spectrum that are adjacent to zones of water absorption (A 1 and A 2).
In the SWIR 2 zone, there is the lowest average change in ∆ in the whole analysed spectrum (3%), so influence of atmospheric correction on these bands is very small. Like SWIR 1, the highest relative change is found for water (6%) and the lowest for agriculture (1%). The most significant changes are also observed on the edges of the SWIR 2 zone for mountains (bands 169-170 and 221-224).
Errors in ∆ estimation for reflectance data were analysed in comparison to errors for radiance data. Absolute and relative differences for all landscape types and the four spectrum zones are combined in Table A4. The biggest change in error values is observed in the VIS zone (0.005). This may be because in this zone, the change in ∆ values is the highest. In the three other zones (NIR, SWIR 1 and SWIR 2), average error is the same-0.001. For all zones, the absolute average change in error values is highest for urban landscape and lowest for water. In most cases, error values increase after atmospheric correction, notably for all bands in urban and water landscapes.

Discussion
To the best of our knowledge, in this paper, for the first time, global multifractal parameters are used for hyperspectral data analyses. For this reason, our results may be discussed only by reference to the global multifractal analysis of multispectral and panchromatic images or by reference to the global fractal analysis of hyperspectral data (listed in Table 1). In this section, we discuss the issues of atmospheric absorption, landscape type separability, comparison to fractal dimension and statistical moments and indicate potential directions for further analyses.

Influence of Atmospheric Absorption
The influence of atmosphere components on data in this experiment is analysed in two ways. Firstly, variability of ∆ curves may be compared to the transmittance spectrum of the atmosphere (Figure 2 in [49]). Secondly, comparing ∆ for radiance and for reflectance, we may observe the magnitude of changes and indirectly study how atmospheric correction modifies the data of specific spectral bands.
The biggest influence of atmosphere components on acquired radiance and, in consequence, on pixel values is visible in the A 1 and A 2 zones. An analysis of the shape of ∆ curves in Figure 2 highlights that the smoothness of the curves is disrupted in zones A 1 and A 2. For mountain and water landscapes, ∆ values in A 1 and A 2 zones exceed the values of neighbouring zones. In the case of urban landscape, we observe the opposite situation (∆ values are lower than in adjacent zones).
Values seem to preserve their continuity mostly in agricultural areas. Visual interpretation of these bands confirmed that images of other landscape types are noisier in these two zones. These differences may be related to different levels of humidity in the atmosphere over specific types of landscape and different image acquisition dates.
In the VIS zone, after band 35, ∆ decreases rapidly until band 40, which corresponds to O 2 absorption. In the NIR zone, three subzones can be clearly observed (bands 41-63, 64-84 and 85-103). This division may be caused by the absorbance of solar radiation by H 2 O vapour at the edges. The values for mountains reveal two fragments where there is a sudden increment (bands 63-65 and 81-86). These parts of the electromagnetic spectrum are where absorption caused by H 2 O is observed. Atmospheric correction causes smoothing of the ∆ curve for the mountain landscape in these bands. It could suggest that the ∆ parameter may be used to evaluate atmospheric correction quality. Moreover, in the bands adjacent to these zones (in SWIR 1 and SWIR 2), we observed increased values only for agriculture and mountain landscapes. It may be related to the fact that these two landscape samples include plants. After atmospheric correction, ∆ for these bands also seems to be smoothed.
It is worth noting that plots with mean and variance values of specific bands also reveal the influence of atmosphere on the data acquired (compare Figure 2 with Figure 2 in [49]). However, these two parameters describe only pixel values' distribution, while ∆ describes spatial and spectral complexity at the same time.

Landscape Types and Dimensionality Reduction
It is possible to identify parts of the electromagnetic spectrum where certain landscape samples have different ∆ values and differences between them are higher than values of errors. In such cases, we will note that ∆ values of certain bands could be used to clearly distinguish landscape types. Distinct ∆ values with differences smaller than error values could possibly distinguish landscape types. In the VIS zone, clear separation of all landscape types is observed in bands 23-38. In the NIR zone, it is possible to clearly distinguish landscape types using bands from the second and third subzones, although differences in ∆ are lower than in the VIS zone. Similar values of ∆ for neighbouring bands can be explained by the fact that these bands represent adjacent wavelengths. Low variability in ∆ within specific landscape types (except agriculture) indicates that information capacity, from the point of view of multifractal analysis, is similar. Therefore, ∆ could possibly be useful for data dimensionality reduction (especially in NIR) if agricultural land is not of interest.
In SWIR 1, no overlap between ∆ values (±errors) indicates that all four types of landscape could be separated using bands 115-117 or 121-125. In SWIR 2, all four types of landscape could be separated using almost any band from this zone. The exception are bands 182-187, where ∆ values for mountains are very close to those for urban landscape.
In this research, single samples of four landscape types have been used. To confirm if specific a landscape type is characterised by the same or similar ∆ curves, more samples should be analysed. To assess the usefulness of ∆ for dimensionality reduction, additional experiments should also be performed.

Size of Images and Spatial Resolution
The size of the analysed images is related to the methodology used to estimate the degree of multifractality and influences the range that is used for multifractal scaling (3). Our review of the literature found that sizes ranging from 17 × 17 to 512 × 512 pixels have been analysed (Table 1). However, only Myint et al. [21] investigated different sizes of hyperspectral subsets representing the same land cover class, and those subsets did not represent exactly the same area. Moreover, the data were analysed using a fractal, and not a multifractal, formalism.
In our previous work, we calculated the degree of multifractality for several subset sizes: 1024 × 1024 in [25], 512 × 512 in [22,26,47], 256 × 256 in [16,46] and 64 × 64 and 32 × 32 in [45]. These earlier studies indicated that the bigger the size, the more reliable the multifractal scaling and ∆ estimations. To verify this hypothesis, image subsets were resampled from their original size to 256 × 256, 128 × 128 and 64 × 64 pixels and analysed. In consequence, we obtained images with pixel sizes of 30, 60 and 120 m.
Stable results were found for the degree of multifractality determined for different image sizes, for all analysed landscape types (see Figure 4). However, estimation errors increased as image size decreased. A detailed analysis was performed with absolute and relative mean values of change in ∆ (and errors), and the results are presented in Table 5.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 21 0.019 (39%) for the smallest subsets (64 × 64). The biggest absolute differences are noted for urban landscape, 0.023 (256 × 256), 0.033 (128 × 128), 0.037 (64 × 64), and the smallest for water (0.005, 0.008 and 0.011, respectively). In most cases, resampling the agriculture subset to 256 × 256 or 128 × 128 caused an increment in Δ values. The exception was 64 × 64, where Δ tended to fall. In the analysis of mountain and urban landscapes, it was found that Δ fell almost in all bands. In the case of the water subset, each resampling increased Δ values. To analyse how errors in the Δ estimation change when images are resampled, we compared absolute relative values. In general, we observe that the smaller the image size, the bigger the change in error. The smallest average absolute change in errors was 0.003 (for water), and the highest was 0.009 (for urban landscape).   A smaller subset size changes ∆ values in all landscapes. Generally, the smaller the subset (and lower spatial resolution), the bigger the changes compared to the original image (512 × 512 pixels). Average absolute change increases evenly by 0.005, starting at 0.010 (15%) for 256 × 256, and rising to 0.019 (39%) for the smallest subsets (64 × 64). The biggest absolute differences are noted for urban landscape, 0.023 (256 × 256), 0.033 (128 × 128), 0.037 (64 × 64), and the smallest for water (0.005, 0.008 and 0.011, respectively). In most cases, resampling the agriculture subset to 256 × 256 or 128 × 128 caused an increment in ∆ values. The exception was 64 × 64, where ∆ tended to fall. In the analysis of mountain and urban landscapes, it was found that ∆ fell almost in all bands. In the case of the water subset, each resampling increased ∆ values.
To analyse how errors in the ∆ estimation change when images are resampled, we compared absolute relative values. In general, we observe that the smaller the image size, the bigger the change in error. The smallest average absolute change in errors was 0.003 (for water), and the highest was 0.009 (for urban landscape).

Correlation with Statistical Moments
Following Combrexelle et al. [30], we calculated the coefficient of determination (r 2 ) to evaluate the complementarity between ∆ and statistical moments. In the first assessment, values of ∆ for all bands for the four landscape types in reflectance images were combined. Low values of r 2 were obtained with four moments: mean (0.003), standard deviation (0.225), skewness (0.114) and kurtosis (0.034). In the next step, each landscape type was analysed separately (Table 6). In this case, we observed the highest coefficient values for urban landscape, notably for the mean (0.834) and skewness (0.858). This may be because the urban sample has the most complex structure and contains a dense mixture of dark and bright objects (see Figure 1c). Coefficient values were much lower for the three other landscape types. In general, the comparison performed here of statistical and multifractal description of hyperspectral data states the confirmation our previous analysis with panchromatic bands ( [26,35,47]) and indicates that ∆ provides complementary information to moments and may be used as an additional classification feature. This is in agreement with the fact that multifractal analysis considers both positive and negative higher-order moments (see Section 3), while in the case of statistical description, only positive or particular moments are used.

Comparison with Fractal Dimension
The data used in this experiment were previously analysed with a fractal parameter-fractal dimension (D F )-in [21], where the differential box-counting method was applied [3]. Values of D F (±error of estimation) from that paper are combined in Figure A2 (for radiance) and Figure A3 (for reflectance). Similarly to ∆ values' shape, D F values also increase with wavelength in VIS. There, however, water reached the highest values, then urban, agriculture and mountains. It may be confirmation of the fractal character of water, observed previously on panchromatic VHR images inter alia in [47]. When D F was calculated with other methods (isarithm and triangular prism) in [20], urban landscapes had also higher values than rural ones. Moreover, the few first bands had very high D F values. It is worth noting that the image samples used in [20] present more heterogeneous landscape samples than those in [22] and in this paper.
Unlike ∆, D F values in [22] did not allow for clear separation of landscape types in any band because of overlapping error bars; in [20], errors were not presented. It reveals an advantage of using the multifractal over the fractal parameter.
Both ∆ and D F show the increment in values caused by atmospheric correction in the VIS zone. In NIR, three subzones indicated by ∆ values were also visible when D F was used in [22]. Moreover, both parameters change at the edges of these subzones (compare our results to Figure 2 in [49]) when atmospheric correction is applied, especially in the case of mountain landscape. In the SWIR 1 zone, the shape of the ∆ curve for agriculture is very similar to the shape of the D F curve presented by [22]. However, in the latter, agriculture has the lowest values from all landscape types. The other three landscapes present almost flat shapes, with the highest values for water, then mountains, urban and agriculture.
Previous research on D F ( [20,22]) revealed disturbances in the A 1 and A 2 zones, similarly to ∆. Using ∆ or D F (estimated with the DBC method), the continuity of curve was disturbed. Values were mostly elevated, but with few exceptions. Using isarithm and triangular prism methods, D F values were clearly elevated for rural and urban samples. It may suggest that both parameters (∆ and D F ) may be used for absorption bands' evaluation.

Conclusions
The research presented here is one of the first attempts to describe hyperspectral aerial images using multifractal parameters. Our results confirmed the multifractal nature of hyperspectral data acquired by AVIRIS. In particular, we calculated, for the first time, the degree of multifractality (∆) for each spectral band of four image samples that present different landscape types. For each image, two levels of data processing were used (radiance and reflectance). The analysis of variability in ∆ highlighted a clear division of the spectrum into the following four ranges: VIS, NIR, SWIR 1 and SWIR 2. Our results are presented in detail for each range. Two zones of absorption were visible in radiance subsets (A 1 and A 2). Here, radiation absorption caused by water particles in the atmosphere leads to significant disturbances in the shape of the ∆ curve. This observation suggests that ∆ could be used as an indicator to discriminate noisy bands.
Based on our results, we can assume that variability in ∆ could possibly be used as a parameter for feature selection and identification of optimal bands to distinguish between different landscape types in the context of context-based image retrieval or image information mining. A comparison with results obtained from fractal analysis [22] indicates that our method identifies more bands where the degree of multifractality can distinguish landscape types. Moreover, the smaller ∆ error estimation confirms that the formalism is better-suited to aerial hyperspectral images analysis. In general, a comparison of the usability of ∆ to distinguish landscape types highlights that better performance is achieved using radiance data. The exception is the NIR zone, where reflectance data provide more bands with clear interclass separability than radiance. Moreover, in most cases, ∆ estimates increase after atmospheric correction.
In some bands, large differences in ∆ were found for data before and after atmospheric correction. The increase in ∆ after atmospheric correction shows that the parameter could be used to assess atmospheric correction quality, especially given that one value of ∆ may be calculated for the whole image. Both ∆ and its error can be used to accurately detect noisy bands. Like D F used in [22], both approaches (fractal and multifractal) indicate that atmospheric correction causes the biggest increase in parameter values in the visible part of the electromagnetic spectrum. Additionally, data acquired in absorption zones (A 1 and A 2) cause disturbances in both fractal and multifractal parameters, and smoothing effects are found in bands 62-64 and 81-85 for mountain landscape. However, unlike D F , the degree of multifractality can detect bands dominated by 0-value pixels, which should be excluded from the study.
Variability in ∆ does not change significantly with decreasing subset size. Experiments with images of four different sizes (512 × 512, 256 × 256, 128 × 128 and 64 × 64) revealed that variability in ∆ is stable, regardless of landscape type and size. Moreover, an analysis of change in ∆ (and errors) found that the smaller subset size, the bigger the differences in ∆ and its error. The biggest absolute differences in ∆ (and its error) were noted for urban landscape, and the lowest for water.
Another interesting observation is the generally low correlation coefficient values found between ∆ and statistical moments (mean, standard deviation, skewness and kurtosis). This finding suggests that ∆ may be used to complement the description of hyperspectral data.
Finally, it should be noted that our research introduces a new parameter for the global description of hyperspectral data. To sum up, we may conclude that the performed multifractal analysis resulted in several interesting findings and seems to be useful for a number of possible applications. Further studies on larger and more varied datasets are required to confirm our conclusions. These studies will be performed in the next stage of our work, during which we will analyse different land cover/use classes and images from other hyperspectral sensors. Another interesting avenue will be to apply a different method to estimate ∆ and different feature extraction methods. We believe that these studies may make a significant contribution to the problem of dimensionality reduction in hyperspectral data.  . Values of first four statistical moments (from top: mean, variance, skewness and kurtosis) determined for radiance (left column) and reflectance (right column) data for four landscape types: water (blue), mountains (black), agriculture (green) and urban (red). Figure A1. Values of first four statistical moments (from top: mean, variance, skewness and kurtosis) determined for radiance (left column) and reflectance (right column) data for four landscape types: water (blue), mountains (black), agriculture (green) and urban (red).  Table A3. Errors in the ∆ estimation for radiance data. Bold font indicates the highest mean value within a zone; green represents the lowest mean error for a specific landscape type; and red, the highest.