Application of various approaches of multispectral and radar data fusion for modelling of aboveground forest biomass

The main source of data for ground-based estimation of Ukrainian forest biomass parameters is the forest inventory carried out by forest enterprises. Previous assessments of Ukrainian forest live biomass have been undertaken by various researchers (Lakyda 2002; Lakyda et al. 2011, 2012). In these cases, the aboveground forest biomass has not been measured directly during the forest inventory but computed using allometric models (Shvidenko et al. 1987, 2007; Lakyda 2002; Lakyda et al. 2012).

Algorithms for modelling aboveground forest biomass have been developed by Lakyda et al. (2011). The models for the assessment of live aboveground biomass and the carbon content of different components of the forest stand are based on the relations between forest productivity and the main inventory parameters of the forest stand (e.g. age, diameter and height) (Lakyda et al. 2012; Cortés et al. 2014; Lu et al. 2004; Singh et al. 2014). It has been shown that the relation between live biomass and the main biometric indicators of forest stands can be expressed as follows: 1 Mi=f(A,D,H,P) \[{{M}_{i}}=f(A,D,H,P)\] where:

M - aboveground biomass (t/ha),

For example, aboveground biomass of deciduous forests of the Ukrainian Polissya was estimated using the following equation (Lakyda et al. 2012): 2 Mi=ai×Dbi×Hci×Pdi \[{{M}_{i}}={{a}_{i}}\times {{D}^{{{b}_{i}}}}\times {{H}^{{{c}_{i}}}}\times {{P}^{{{d}_{i}}}}\] where:

a_i, b_i, c_i, d_i – regression coefficients calculated for each biomass fraction of forest (i).

The biomass for the main fractions of trees (stem (st), branches (br) and leaves (lv)) has been calculated using equation (2), while the biomass of the crown (M_cr) and the whole aboveground tree biomass (M_tr) is calculated using the following: 3 Mcr=Mbr+Mlv \[{{M}_{cr}}={{M}_{br}}+{{M}_{lv}}\] 4 Mtr=Mst+Mcr \[{{M}_{tr}}={{M}_{st}}+{{M}_{cr}}\]

To convert this biomass into carbon, coefficients of carbon content in the corresponding fraction are utilized. From the literature, Lakyda et al. (2012) have assumed values of 0.5 and 0.45 for wood and leaves respectively. Thus, the equations for calculating carbon stock in the forest stand are given as follows: 5 Mcrc=0.5Mbr+0.45Mlv \[M_{cr}^{c}=0.5{{M}_{br}}+0.45{{M}_{lv}}\] 6 Mtrc=0.5(M|st+Mbr)+0.45Mlv \[M_{tr}^{c}=0.5(M|st+{{M}_{br}})+0.45{{M}_{lv}}\]

Optical remote sensing of forests is based on developing relationships between spectral signatures and measured parameters. Multispectral satellite images provide the spatial distribution of land cover reflectance in the visible, near-infrared and infrared spectral ranges. These reflectance data can be used directly as well as for VI calculations. Depending on the complexity, such indices can be divided into simple ratios (SR), normalized differences (ND) and complex vegetation indices (Lu et al. 2004). A wide set of VIs has already been proposed to date; some of these are listed in Table 1. The theoretical basis for empirical VIs is derived from an examination of the typical spectral reflectance signatures of leaves. The reflected radiation in visible bands is very low as a result of high absorption by photosynthetically active pigments, with maximum absorption in blue (470 nm) and red (670 nm) wavelengths. Near-infrared radiation (NIR) is scattered (reflected and transmitted) with very little absorption in a manner dependent upon the structural properties of a canopy (leaf area index (LAI), leaf angle distribution and leaf morphology). As a result, the contrast between red and near-infrared responses is a sensitive measure of vegetation amount, with RED/NIR maximum over a dense canopy and minimal contrast over rare or no vegetation. A major finding on atmospheric effect minimization is the use of the difference in blue and red reflectances as an estimator of the atmosphere influence level. This concept is based on the wavelength dependency of aerosol scattering cross sections. In general, the scattering cross section in the blue band is larger than that in the red band. When the aerosol concentration is higher, the difference in the two bands becomes larger. This information is used to stabilize the index value against variations in aerosol concentration levels.

The enhanced vegetation index (EVI) incorporates this atmospheric resistance concept in the atmospheric resistant index (ARVI), along with the removal of soilbrightness-induced variations in VI as in the soil adjusted vegetation index (SAVI). The EVI additionally decouples the soil and atmospheric influences from the vegetation signal by including a feedback term for simultaneous correction.

Lu et al. (2004) showed that the relationships between TM reflectance and forest stand parameters vary depending on the characteristics of the study areas. However, not all vegetation indices are significantly related to forest stand parameters. What is crucial is the selection of suitable TM band (s) and VIs for relevant biophysical parameter estimation.

Unlike optical remote sensing, which allows for estimating only the top of the forest cover and does not provide any information about forest structure (due to impermeability of the optical waves), radar remote sensing uses the microwave portion of the electromagnetic spectrum. Canopy penetration varies with different wavelengths. Shorter wavelengths (e.g. X-band imagery at 3 cm) are reflected from the top of the canopy, while longer wavelengths (e.g. L-band imagery at 24 cm) go down to the ground and are reflected from there. Using these properties of different wavelengths makes it possible to discern information about canopy structure of a forested area from a multi-wavelength image and thus estimate the different component of AGB.

One more useful feature of radar data for studying forest structure is polarization. Transmitted and received radar signals are propagated in a certain plane. The propagation planes are usually horizontal (H) and vertical (V). Vertically polarized waves interact with the vertical stems of the forest cover, while horizontally polarized waves penetrate through the canopy. Thus, the combination of images with different channels can provide additional information about forest structure.

As with optical remote sensing data, two sets of texture features can be calculated using SAR data. The first set of calculations can be derived from the backscatter (sigma nought, σ⁰) distribution and the second one from the grey-level co-occurrence matrix. The intensity scenes of SAR images are converted in their corresponding backscattering coefficient (σ°) values using the following equation (Attarchi et al. 2014; Shimada et al. 2009): 7 σ0=10×log10(I2+Q2)+CF−32.0 \[{{\sigma }^{0}}=10\times {{\log }_{10}}({{I}^{2}}+{{Q}^{2}})+CF-32.0\] where:

CF – calibration factor = −83 dB,

To take into account the effects of relief, the topographic normalized backscattering coefficient (the corrected backscatter in gamma-nought γ°) can be obtained from the sigma-nought σ° value according to Ulander (1996), Castel et al. (2001) and Santoro et al. (2011): 8 γ0=σ0×AflatAslope×(cosθrefcosθloc)n where:

σ⁰ – the radar backscattering coefficient,

The exponent n is the optical canopy depth and ranges between 0 and 1. It is a site-specific factor and difficult to obtain in practice; therefore, it is set to 1 (Thiel et al. 2009; Kim 2012; Santoro et al. 2011; Attarchi et al. 2014).

Previous studies have shown that the forest backscatter can be described as a function of the growing stock volume, V (Pulliainen et al. 1994; Santoro et al. 2011): 9 σforo=σgro×e−βV+σvego(1−e−βV) \[\sigma _{for}^{o}=\sigma _{gr}^{o}\times {{e}^{-\beta V}}+\sigma _{veg}^{o}(1-{{e}^{-\beta V}})\]

The backscatter model in Eq. (9) contains three unknowns that need to be estimated: σgro$\sigma _{gr}^{o}$, σvego$\sigma _{veg}^{o}$ and β. These can be determined by means of least-squares regression, using a dataset of reference forest growing stock volume (GSV) measurements (Pulliainen et al. 1994; Santoro et al. 2002, 2013). After backscatter coefficients and texture features have been calculated, we can estimate the forest GSV.

Data fusion is a general multi-discipline approach. It combines data from multiple sources to improve the potential values and interpretation performances of the source data and to produce a high-quality visible representation of the data.

In general, remote sensing fusion techniques can be classified into three different levels: the pixel/data level, the feature level and the decision level (Pohl and van Genderen 1998; Zhang 2010). Pixel-level fusion is the combination of raw data from multiple sources into single resolution data, which is expected to be more accurate than either of the individual input data or it may reveal the changes between data sets acquired at different times (Zhang 2010).

Feature-level fusion extracts various features, e.g. edges, corners, lines and texture parameters, from different data sources and then combines them into one or more feature maps that may be used instead of the original data for further processing. This is particularly important when the number of available spectral bands becomes so large that it is impossible to analyse each band separately. Typically, in image processing, such fusion requires a precise (pixel-level) registration of the available images. Feature maps thus obtained are then used as inputs to pre-processing for image segmentation or change detection (Zhang 2010).

Decision-level fusion combines the results from multiple algorithms to yield a final fused decision. When the results from different algorithms are expressed as confidences (or scores) rather than decisions, it is called soft fusion; otherwise, it is called hard fusion. Methods of decision fusion include voting methods, statistical methods and fuzzy logic-based methods (Zhang 2010).

There are a set of multi-source remote sensing data fusion methods, which are based on different techniques: Markov random field (MRF) (Solberg et al. 1996), support vector machines (SVM) (Waske and Benediktsson 2007) and the decision fusion approach for multi-temporal classification (Jeon and Landgrebe 1999). However, the different nature and content of remote sensing imagery and GIS data prevent a direct comparison. Therefore, the integration of data from different applications must address the differences in the object model and the semantics of the objects themselves (Zhang 2010). Images are usually composed of a raster of pixels representing the intensities, whereas GIS data contain artificial objects (points, lines and polygons) with label forms, representing the objects or region affiliations. To combine segmented objects or primitives from remote sensing images and GIS data at the feature level or decision level, traditional pattern recognition methods can be used and have demonstrated their potential capabilities, e.g. knowledge-based techniques (Amarsaikhan and Douglas 2004), neural network and statistical approaches (Benediktsson and Kanellopoulos 1999), fuzzy set theories (Fauvel et al. 2006), Bayesian techniques and Dempster-Shafer-based methods (Zhang 2010).

HPF is one of the first developed image merge methods (Schowengerdt 1980, 2007a, b; Chavez et al. 1991; Pohl and van Genderen 2016), which has been used for image fusion for more than 30 years. It was primarily designed to improve the resolution of multispectral images by transferring spatial details derived from a higher resolution PAN or SAR image to a lower resolution multispectral image. This method is performed in three stages (Pohl and van Genderen 2016):

High-pass filtering the high-resolution SAR image

Thus, this method extracts high-frequency information from the SAR image, which is then added to each spectral channel of the multispectral image.

Figure 1a shows the result of merging the PALSAR image (HH mode) with RapidEye multispectral image using HPF method.

Figure 1.

IHS transform is one of the most popular methods of merging remote sensing images. It uses a mathematical colour model based on a cylindrical or spherical coordinate system. This method effectively separates spatial (I) and spectral (H, S) information from a standard RGB image (Pohl and van Genderen 2016).

There are two ways to use IHS transformation to merge images: direct and substitutional. The first way is to directly convert the three spectral channels of the image into I, H and S. The second method involves converting the three spectral channels from the RGB into IHS colour model, where colour aspects are separated by its average brightness. The hue and saturation in this case are related to the surface reflectivity or composition. The SAR image then replaces one of the components, and reverse transformation from IHS to RGB converts the data into their original colour model to produce a new synthesized image.

Figure 1b shows the result of merging the PALSAR image (HH mode) with a multispectral RapidEye image using IHS method.

Wavelet transform is another powerful mathematical tool used for signal processing. It is used to decompose a function or signal into several components (coefficients). When using this method to merge images, the main idea is to decompose the original images into a number of fragments using direct wavelet transform. The information is merged on the basis of the obtained coefficients of these fragments, and the inverse wavelet transform is applied to synthesize new image. This method is suitable for merging images from different sources with varied physical backgrounds (such as radar and optical images) because it decomposes images into different types of coefficients, preserving the source information. Based on these coefficients, a new image can be synthesized using inverse wavelet transform (Pajares and Cruz 2004).

Figure 1c shows the result of merging the PALSAR image (HH mode) with RapidEye multispectral image using WT method.

Hybrid data merging methods are widely used to compromise between spatial and spectral optimization. WHIS combines the IHS transformation approach with methods used to merge data with different spatial resolution, such as WT (Chibani and Houacine 2002; Gonzalez-Audicana et al. 2004; Zhang and Hong 2005). During wavelet-IHS transformation (WIHS), the multispectral data are converted to IHS colour model. Then, the intensity component I is decomposed using WT method. The SAR image is matched to I and then also decomposed by the WT method. The decomposed component I is replaced by the SAR decompositions, and the inverse IHS transformation is performed (Pohl and van Genderen 2016).

Figure 1d shows the result of merging the PALSAR image (HH mode) with RapidEye multispectral image using WHIS method.

The field data set of forest AGB was divided into two samples with the proportion of 70/30. The first sample (70% of data) was used to build models, which accuracy was further assessed using the second sample (30% of data). Regression analysis methods were used to find a set of relationships between forest AGB and varied combinations of remote sensing data for each approach. As a final result, one model that has the highest accuracy was selected for each approach.

In case of using HPF approach, the best result was achieved for merged HH PALSAR band with 4^th RapidEye spectral band. It is described by the following equation (Fig. 2): 10 AGB=9877.69×e−47.08×HPF(HH+B4) where:

HPF(HH + B4) – the signal from fused HH PALSAR band with 4^th RapidEye spectral band using HPF approach.

Figure 2.

In case of using IHS approach, the best result was achieved for merged HH PALSAR band with 4^th RapidEye spectral band. It is described by the following equation (Fig. 3): 11 AGB=1.3701e+16×IHS(HH+B4)11.15 where:

IHS(HH+B4) – the signal from fused HH PALSAR band with 4^th RapidEye spectral band using IHS approach.

Figure 3.

In case of using WT approach, the best result was achieved for merged HH PALSAR band with 4^th RapidEye spectral band. It is described by the following equation (Fig. 4): 12 AGB=35110+e−59.62×WT(HH+B4) where:

WT (HH+B4) – the signal from fused HH PALSAR band with 4^th RapidEye spectral band using WT approach.

Figure 4.

In case of using hybrid WIHS approach, the best result was achieved for merged HH PALSAR band with 4^th RapidEye spectral band. It is described by the following equation (Fig. 5): 13 AGB=19133.31×e−53.01×WIHS(HH+B4) where:

WIHS (HH+B4) – the signal from fused HH PALSAR band with 4^th RapidEye spectral band using WIHS approach.

Figure 5.

In case of using MLR approach, the best result was achieved for merged HH PALSAR band with 4^th RapidEye spectral band. It is described by the following equation: 14 AGB=308.53−2530.68B4+211.55HH where:

B4 is the reflection coefficient of the 4^th spectral band of the RapidEye image and HH is the amplitude of the reflected radar signal in the HH mode.

The accuracy of the obtained models was assessed with a test data set, which was not used for the models building. The modelled data were compared with the data of field measurements, and the correlations between measured and modelled data were calculated for each model (Tab. 2).

The correlation shows the relationship between the measured and modelled values. The index lays in a range from −1 to 1. Higher correlation means better model performance. A correlation of 0 indicates no relationship between the measured and modelled values.

Also, to assess the accuracy of each model, the min–max accuracy (MMA) was calculated: 15 MMA=mean(min(actuals,predicteds)max(actuals,predicteds))

This parameter shows how close the predicted values are to the actual ones. The MMA value ranges from 0 to 1, where a value of 1 indicates a perfect match between actual and predicted values. So, higher MMA score means better model performance.

Among the different methods of data fusion used in the research, the HPF method (correlation 0.804, MMA 0.652) showed the best result with a significant advantage over other methods.