Estimation of Postfire Reforestation with SAR Polarimetry and NDVI Time Series

Abstract

This communication is devoted to the methodology of remote complex analysis of forest restoration after strong wildfires. It is proposed to quantify the projective leaf/needles area index by multispectral optical images. The increase in dimensions of trunks and branches commensurate with a radar wavelength of 24 cm is estimated using radar polarimetric data. It is shown that the growth’s potential of aboveground biomass in different spots of test site ranges from 35 to 70% in the case under consideration. Such a new approach will make it possible to further consider more accurately the role of boreal forests as one of the largest carbon stocks.

1. Introduction

Modern climate changes [1] are largely associated with the human-induced problem of carbon balance. This problem resulted from emissions associated with energy use, urbanization, and land use changes. Today there has been progress in monitoring and accounting for carbon sources and carbon sinks [2,3,4] that shows a key role of forests in the global carbon cycle. Unfortunately, massive deforestation caused by long-standing wildfires and cuttings reduces the volume of carbon sinks. In recent years, this is especially true for the Siberian boreal forest [5,6,7]. One of the scenarios of significant carbon sinks increasing is massive reforestation on the areas affected by strong wildfires, as the young growth can fix large amounts of carbon in their woody biomass. Accordingly, a large-scale inventory of the current state of the forest is needed with retrospective estimates of the spatial distribution of the species composition and biomass of the forest before wildfires.

The fastest and most effective method for such assessments is remote sensing using optical sensors [5,6,7,8] and polarimetric radars [3,9], as well as comprehensive remote optical and microwave measurements of biophysical forest parameters [10]. At the same time, regularly increasing volumes of satellite data can be processed using standard calculation algorithms on remote servers, such as Google Earth Engine (GEE) [11,12,13], and presented as standard geoinformation raster layers. The GEE is a remote computing cluster with a multi-petabyte archive of public aerospace imagery covering the entire world that enables high-performance computations via interactive interface. It allows users to employ well-known geospatial data processing methods as well as to creating their own algorithms and software solutions. Users can request and analyze data from a shared catalog or upload their own data. A technical overview of the platform is presented by Gorelick et al. [11]. However, at present, there are no such platforms for processing and analyzing radar interferometric and polarimetric data. Therefore, specialized software should be used.

This study proposes a methodology for assessing reforestation processes based on the integrated use of polarimetric decomposition of L-band radar data and cloud computing of normalized difference vegetation index (NDVI) values using Google Earth Engine.

2. Satellite Data and the Test Site

For our study, we used quad-polarimetric data of L-band ALOS-2 PALSAR-2 from 5 April 2017 (Scene ID: ALOS2154891030-170405) and 30 August 2019 (Scene ID: ALOS2284561036-190830) as well as 32-Day NDVI composites based on Landsat-5 (1 January 1995–16 October 2011) and Landsat-8 (7 April 2013–6 March 2021) data that are available in the GEE data catalog.

The test site (107.64° N, 51.90° E; Figure 1) is a typical area of Siberian boreal forest that is located in the region of the Baikal Lake, near the city of Ulan-Ude. The area is mountainous with 680–800 m elevations and slopes up to 30°. After the forest fire that occurred in the spring of 2003, reforestation work was carried out, and pine seedlings (Pínus sylvéstris) were planted. There is an intensive reforestation at present (see Figure 2a).

Red polygon in the false color composite of SAR polarimetric data (see Figure 1b) and in the panchromatic optical image of Resurs-P satellite (2 m spatial resolution) (see Figure 1c) of 18 September 2020 shows the position of the reforestation test site with an area of 68.8 ha. The false color composite in Figure 1b gives a qualitative characteristic of a predominant radar backscattering mechanism. Such a representation is generally accepted in the literature on satellite radar polarimetry [14,15,16,17,18], as it has a fairly simple interpretation: shades of red are buildings and various natural objects that form a kind of dihedral scatterer with the surface; shades of green are volume discrete inhomogeneities, forest as usual; dark blue and blue are forestless rough surfaces; black is water and smooth surfaces.

For comparison, sections of mature unburned forest adjacent to the reforestation test site were selected. The control section of the pine forest outlined by a white polygon in Figure 1b,c is similar to the test site before the wildfire. Preliminary studies have shown that correlation between NDVI time series for control and reforestation sites is equal to 0.97 in 1995–2002 (Landsat-5 data). The second control section of mixed forest with a predominance of aspen and birch (yellow polygon in Figure 1b,c) was chosen considering the possibility of replacing a coniferous forest with a deciduous one after a fire. A photograph of the forest within the white polygon is shown in Figure 2b.

3. Polarimetric Data Decomposition and Classification

The most general relationship between the incident and received fields may be expressed in the form of a matrix equation:

$$\left[\begin{array}{c}{E}_H^r\\ E_V^r\end{array}\right]=\frac{e^{ikR}}{R}\left[\begin{array}{cc}{S}_{HH}& S_{HV}\\ S_{VH}& S_{VV}\end{array}\right]\left[\begin{array}{c}{E}_H^i\\ E_V^i\end{array}\right]$$

(1)

or

$$E^r=\frac{e^{ikR}}{R}\left[S\right]E^i$$

(2)

where $E^r$ and $E^i$ are the vectors of received and incident field, respectively, [S] is referred to as the scattering matrix or Sinclair matrix of the scatterer, R is a wave propagation path, and k is a phase constant. It should be pointed out that subscripts H and V of matrix elements mean horizontal (H) and vertical (V) polarizations. The first subscript on each of the elements refers to the polarization of the scattered wave, whereas the second subscript refers to the polarization of the incident wave.

Direct interpretation of scattering matrix is problematic; thus, decomposition methods were developed to represent it as a sum of known radar backscattering mechanisms with more clear interpretation [15,16,17,18].

One of the most common methods is the Cloude–Pottier decomposition [16], also named H-α-classification. The Cloude–Pottier decomposition is based on an analysis of coherence matrix eigenvalues. This matrix is defined by the outer product of a three-dimensional scattering vector written in the Pauli basis:

$$k_P=\frac{1}{\sqrt{2}}\left[\begin{array}{c}{S}_{HH}+S_{VV}\\ S_{HH}-S_{VV}\\ 2S_{HV}\end{array}\right]$$

(3)

Here, the S_HV = S_VH condition is supposed to be fulfilled. The coherence matrix [T] off-diagonal elements demonstrate the degree of k_p vector elements correlation:

$$\langle \left[T\right]\rangle =\frac{1}{2}\langle \left[\begin{array}{ccc}\left(S_{HH}+S_{VV}\right){(S_{HH}+S_{VV})}^{*}& \left(S_{HH}+S_{VV}\right){(S_{HH}-S_{VV})}^{*}& 2\left(S_{HH}+S_{VV}\right)S_{HV}^{*}\\ \left(S_{HH}-S_{VV}\right){(S_{HH}+S_{VV})}^{*}& \left(S_{HH}-S_{VV}\right){(S_{HH}-S_{VV})}^{*}& 2\left(S_{HH}-S_{VV}\right)S_{HV}^{*}\\ 2S_{HV}{(S_{HH}+S_{VV})}^{*}& 2S_{HV}{(S_{HH}-S_{VV})}^{*}& 4S_{HV}{S}_{HV}^{*}\end{array}\right]\rangle$$

(4)

As the matrix [T] is by definition Hermitian positive semi-definite, it can always be transformed to a diagonal form:

$$\left[T\right]=\left[U_3\right]\left[\begin{array}{ccc}{\mathsf{\lambda }}_1& 0& 0\\ 0& {\mathsf{\lambda }}_2& 0\\ 0& 0& {\mathsf{\lambda }}_3\end{array}\right]\left[U_3^{*T}\right]$$

(5)

where λ₁ ≥ λ₂ ≥ λ₃ are the eigenvalues of the coherence matrix that are real non-negative numbers, and the $\left[U_3\right]$ unitary matrix consists of its eigenvectors. As the dominant scattering mechanism in the model, the mechanism corresponding to the eigenvector with the maximum eigenvalue λ₁ is taken; the remaining two vectors describe the second and third most frequent scattering mechanisms.

The degree of scattering randomness is described by the entropy

$$H={\displaystyle \sum }_{j=1}^3-P_j{\log }_3{P}_j$$

(6)

where P_j is the probability of each of the three mechanisms:

$$P_j=\frac{{\lambda }_j}{{\lambda }_1+{\lambda }_2+{\lambda }_3}$$

(7)

If only one of the eigenvalues of the coherence matrix is non-zero, a single scattering mechanism is observed with probability one (deterministic case), and H = 0. If there is no allocated scattering mechanism, and λ₁ = λ₂ = λ₃, then H = 1. However, both these extreme cases never occur in reality, and in-between situations are commonly observed.

We can consider ${\lambda }_1\gg {\lambda }_2\cong {\lambda }_3$, as some approximation to H = 0. A typical example is the Bragg scattering from a rough surface. An approximation to H = 1 is the scattering on a cloud of chaotically oriented dipoles that is a model for dense forest crowns. High entropy values indicate the presence of vegetation.

The α parameter in the Cloude–Pottier decomposition is calculated as a weighted mean

$$\overline{\alpha }=P_1{\alpha }_1+P_2{\alpha }_2+P_3{\alpha }_3$$

(8)

where ${\alpha }_j$ is one of the angles that parametrize the $e_j$ eigenvector:

$$\left[e_j\right]=e^{i{\varphi }_j}{\left[\begin{array}{ccc}\mathit{cos}{\alpha }_j& \mathit{sin}{\alpha }_j\mathit{cos}{\beta }_j{e}^{i{\delta }_j}& \mathit{sin}{\alpha }_j\mathit{sin}{\beta }_j{e}^{i{\gamma }_j}\end{array}\right]}^T$$

(9)

The $\alpha$ angle takes on values from 0° (corresponds to the surface scattering mechanism) to 90° (corresponds to the double reflection from ideal conductive surface of a dihedral scatterer). An intermediate value of 45° corresponds to dipole scattering.

Using this method, sixteen different classes can be distinguished that are components of the three main physical mechanisms of the radar backscattering (see Figure 3):

• Single bounce (surface) scattering (no tree trunks and branches comparable with L-band wavelength, thus there is scattering from the surface only);
• Double bounce (dihedral) scattering (e.g., from the surface and further from the tree trunk or in reverse order;
• Volume scattering by a set of randomly located branches, tree trunks, having the dimensions comparable or greater than the L-band wavelength.

4. Results of Polarimetric Data Analysis

Let us consider the possibility of using this decomposition to assess forest recovery after a fire using the analysis of multi-temporal L-band polarimetric data. In this case, for the sake of simplicity, we evaluate changes in the three main physical scattering mechanisms, as was shown by Bondur et al. [17]. It should also be noted that the change in scattering mechanisms would be associated with the dynamics of the size of heterogeneities in forest regrowth. For example, individual trunks and branches, as well as a combination of trunks and branches, become comparable with a wavelength of 24 cm.

A classification using the Cloude–Pottier method was carried out according to radar polarimetry data for the years of 2017 and 2019. The comparison results are given in Figure 4. The colors of test site polygons are similar to polygons from Figure 1. Areas of green are reforestation areas, where the predominant mechanism changed from surface to volume scattering (see Figure 3).

Figure 5a shows a fragment of H-α classification result according to 2019 data, where classes with dominance of surface scattering are colored in shades of blue, classes with dominance of volume scattering are in green, and classes with dominance of double scattering are in red.

For further analysis, the test reforestation area will be represented as 15 separate sites (Figure 5b) with different regrowth density. Sites 1–7 have high density, 8–11 have medium density, and 12–15 have low density. A binary vegetation mask obtained from a high-resolution optical Resurs-P winter image of 4 March 2018 was used for site selection. The territory covered by young pine stock is green (NDVI ≥ 0.2, mask value is 1), without vegetation is white (NDVI < 0.2, mask value is 0).

Let us consider the statistics of the results of polarimetric decompositions of the ALOS-2 PALSAR-2 radar data for 2017 and 2019, grouped according to the type of dominant scattering mechanism (see

Table 1. H–α classification statistics for reforestation test sites and control pine and mixed forest sections.

Test Sites	Projective Leaf Area Index	Surface, %		Volume, %		Double, %		Field Measurements
		2017	2019	2017	2019	2017	2019	Trees/m2	Mean Trunk Diameter, m	Average Height, m
Site 1	0.73	83.71	75.93	16.29	24.07	0	0	2.35	0.11	3.3
Site 2	0.71	86.20	74.32	13.80	25.68	0	0	2.20	0.09	3.1
Site 3	0.71	98.51	88.85	1.49	11.15	0	0	2.20	0.04	2.3
Site 4	0.82	82.74	65.79	17.26	34.21	0	0	2.50	0.08	3.2
Site 5	0.81	93.07	83.73	4.27	16.27	2.67	0	2.50	0.10	3.8
Site 6	0.72	70.41	73.44	29.50	25.85	0.09	0.71	2.25	0.11	3.3
Site 7	0.67	93.61	92.01	5.95	7.72	0.44	0.27	2.00	0.08	3.0
Site 8	0.18	83.23	71.26	16.53	28.20	0.24	0.54	1.05	0.12	4.5
Site 9	0.24	86.11	85.32	13.69	14.68	0.20	0	1.35	0.06	2.9
Site 10	0.24	85.67	78.81	14.33	21.19	0	0	1.25	0.07	3.1
Site 11	0.36	86.96	85.44	13.04	14.56	0	0	1.55	0.05	2.3
Site 12	0.02	99.22	99.74	0.78	0.26	0	0	0.45	0.07	2.5
Site 13	0.03	87.73	72.86	12.27	25.84	0	1.30	0.55	0.09	4.2
Site 14	0.10	91.46	74.84	8.54	25.16	0	0	0.25	0.08	3.0
Site 15	0.02	84.38	89.79	12.91	10.21	2.70	0	0.45	0.07	2.7
Pine forest	0.98	3.88	3.28	70.77	77.92	25.35	18.79	0.13	0.30	19
Mixed forest	1	5.85	11.03	87.40	85.48	6.75	3.49	0.14	0.18	16

). Test and reference sites in this table are ranked according to the density of the forest canopy. The density was calculated as an average value over the territory of each test area on the binary image of the mask and represented the proportion of the territory occupied by forest vegetation. This density varied from 0—in the absence of vegetation to 1—for continuous dense forest. Data for 2017 and 2019 are given in the table. Changes in the values in the table show the dynamics over the 2.5-year period between surveys.

In general, throughout the study area of reforestation (red polygon), there is a decrease in the surface scattering component from 84.43% to 79.51%, and the dominance of the volume type of scattering increases from 12.44% to 20.38%. The dominant double scattering area is not large and almost does not change (from 0.13% to 0.11%), which is probably due to the insufficient thickness and weight of trunks of 10–15-year-old pine regrowth in the L-band.

As we can see from Table 1, the dynamics is the same for most reforestation sites. At the same time, it should be noted that there is a significant difference between the values corresponding to the volume scattering of the reference pine and mixed forest (70–87%) and those for the considered reforestation site (34% max). This indicates the potential for further tree growth.

5. Reforestation Assessment Based on the Results of Cloud Computations

To assess the reforestation process, the results of cloud computations (using GEE) of NDVI time series were used, based on the data of the Landsat-5 and Landsat-8 satellite sensors averaged over 32 days [19]. The relative regrowth index (RRI) was calculated using these values, according to the following equation [20]:

RRI = NDVI_test − NDVI_pine

(10)

where NDVI_test is the index averaged over the territory of the test site with forest regrowth, and NDVI_pine is the index of the reference pine forest site.

The result is given in Figure 6a. NDVI values of these two sites before the fire differ slightly; the RRI is close to zero. The correlation coefficient for two NDVI series for 1995–2002 is equal to 0.97.

After the fire in 2003, seasonal RRI fluctuations increased. This was related to intensive growth of grass within the test site (continuous projective green plant cover) in spring and summer and snow cover in winter (discrete pine regrowth cover). Thus, to distinguish conifer undergrowth and grass cover, let us consider RRI averaged over winter and summer seasons.

The averaged curves in Figure 6b demonstrate different rates of recovery to the level of the reference forest, i.e., level close to zero. The linear trends demonstrate that a coniferous forest reaches an equal, with a combination of dense grass, shrub vegetation, and deciduous trees, projective plant cover after 10 years or more.

Overall, the considered approach to the use of cloud computing for long-term series of the vegetation index will make it possible to predict the restoration of the forest with the replacement by more valuable coniferous tree species.

6. Conclusions

The technique of integrated remote sensing of reforestation processes after strong wildfires is proposed. The methodology is based on joint use of the following:

Fully polarimetric L-band ALOS-2 PALSAR-2 satellite radar data to assess the temporal dynamics of the physical mechanisms of radar backscatter in the areas of forest regrowth. A decrease in the surface scattering component of the radar backscattering and an increase in the scattering component by the volume inhomogeneity of the forest regrowth means an increase in the size and number of inhomogeneities (trunk and branches) of the forest regrowth, comparable with the radar wavelength.

The results of cloud computations of the vegetation index long-time series according to Landsat-5, 8 multispectral imagery and subsequent assessment of the relative regrowth index (projective plant cover increase) by coniferous regrowth on the basis of the comparison of winter and summer indices of the test and reference sites.

The proposed combination made it possible to quantitatively and qualitatively analyze the process of reforestation with a comparative assessment of the retrospective state of the forest before the fire, as well as from the point of view of seasonal changes in the projective cover and the geometric dimensions of forest regrowth. Long-time NDVI series have allowed us to reveal the forest recovery trend and to forecast further reforestation. Despite the lack of radar polarimetric data time series, the capability of assessing the processes of reforestation in terms of growth in the size of trunks, branches, i.e., aboveground tree biomass, has been demonstrated. In the case under consideration, it is shown that the growth potential of aboveground biomass in different areas ranges from 35% to 70%.

The described approach is novel and will further make it possible to more accurately account for the role of boreal forests, including forest regrowth, as one of the largest carbon stocks.