Application of a Three-Dimensional Radiative Transfer Model to Retrieve the Species Composition of a Mixed Forest Stand from Canopy Reflected Radiation

Abstract

The paper introduces a three-dimensional model to derive the spatial patterns of photosynthetically active radiation (PAR) reflected and absorbed by a non-uniform forest canopy with a multi-species structure, as well as a model algorithm application to retrieve forest canopy composition from reflected PAR measured along some trajectory above the forest stand. This radiative transfer model is based on steady-state transport equations, initially suggested by Ross, and considers the radiative transfer as a function of the structure of individual trees and forest canopy, optical properties of photosynthesizing and non-photosynthesizing parts of the different tree species, soil reflection, and the ratio of incoming direct and diffuse solar radiation. Numerical experiments showed that reflected solar radiation of a typical mixed forest stand consisting of coniferous and deciduous tree species was strongly governed by canopy structure, soil properties and sun elevation. The suggested algorithm based on the developed model allows for retrieving the proportion of different tree species in a mixed forest stand from measured canopy reflection coefficients. The method accuracy strictly depends on the number of points for canopy reflection measurements.

1. Introduction

Solar radiation is both a direct and indirect driver for most biophysical and biochemical processes occurring in plant ecosystems. It influences canopy microclimate as well as CO₂ and H₂O exchange processes between plants and the atmosphere [1,2], thereby affecting not only the function and growth of the plant community, but also the concentration of greenhouse gases in the atmosphere and consequently, Earth’s climate system [3,4]. Plant canopy reflection and absorption of solar radiation are determined by factors including the structure of vegetation cover, the optical properties of photosynthesizing and non-photosynthesizing parts of plants, surface topography, and soil [1,5]. During recent decades, numerous studies have focused on the theory of solar radiation transfer in a plant canopy [1,6,7,8,9,10], leading to the development and application of mathematical models with different degrees of complexity [11,12,13,14,15,16,17]. Most simple approaches consider the plant canopy as a horizontally uniform turbid medium [13,18,19]. A vertically structured forest canopy can also be divided into many independent sub-layers [10]. The effects of local plant canopy heterogeneity on radiative fluxes in these one-dimensional (1D) approaches are usually ignored. A limited number of input parameters allow these models to be used for a broad spectrum of scientific tasks, including solving direct and inverse problems related to remote sensing [10,15,20].

The effects of local-scale heterogeneity on radiation fluxes at a regional scale are usually assumed to be very small and therefore negligible, but they can be relatively large on an ecosystem level and can cause large biases in the estimation of solar radiation that is reflected, absorbed, and transmitted by a plant canopy [21,22,23]. Complex forest canopy architecture is common along the boundaries between different forest types and this complexity can have significant effects on absorbed and scattered solar radiation. Scientific studies of radiative transfer in such non-uniform plant canopies require more sophisticated and process-based modeling approaches that consider three-dimensional (3D) canopy structure, including the individual structure of different tree species.

There are several detailed reviews of available 3D models for the description of radiative transfer within spatially distributed and non-uniform vegetation [7,10,20]. Many of these works were devoted to the ability of 1D and 3D models to reproduce the mean characteristics of canopy absorption, transmission, and reflection by allowing for mosaic vegetation structure [16,17,24]. Most studies show that 3D models of radiative transfer can be a universal tool to describe the 3D radiation absorption and scattering patterns for a spatially heterogeneous plant canopy. These models are presently used for various ecological and meteorological tasks involving satellite remote sensing (e.g., Landsat, MODIS, DSCOVR EPIC) [15,25,26,27,28,29]. Key parameters that can be determined using such models are total absorbed and reflected solar radiation by a plant canopy and their variance, the fractions of sunlit and shaded leaves, leaf area index of sunlit and shaded leaves, and soil surface reflection/absorption, etc. In considering the reflection of solar radiation from the forest canopy with a multi-specific structure, it is important to also estimate the proportion of various tree species within a forest stand. Such information can be derived from other biophysical properties of the forest canopy and soil, such as leaf area index and the amount of photosynthesizing plant biomass.

In this study we developed and applied a 3D model of radiative transfer in a non-uniform plant canopy to derive the possible influence of multi-specific forest structure on reflected photosynthetically active radiation (PAR). For our numerical experiments, we selected a mixed forest stand consisting of a coniferous and a deciduous tree species with different structure of individual trees, vertical and horizontal biomass distribution, and leaf and shoot optical properties. Mixed deciduous–coniferous forest is a dominant forest type at the southern boundary of the boreal forest and effects of forest structural and optical properties on radiative transfer (reflection, absorption, transmission) have not been sufficiently investigated in the area. Moreover, in this study we proposed and tested an algorithm to retrieve forest canopy composition (proportion of deciduous and coniferous tree species in a forest stand) from the intensity of reflected radiation measured by remote sensing at several points above the forest canopy.

2. Materials and Methods

2.1. The 3D Model Description

The 3D radiative transfer model is based on the main equations and assumptions suggested by Ross [1] and further developed and described by Myneni [7,30] and Knyazikhin [22,31]. According to this approach, the function characterizing the radiative field at each spatial point $r=\left\{x,y,z\right\}$ within the vegetation canopy is the radiance $I_{\lambda }\left(r,\Omega \right)$, which depends on wavelength λ and is calculated as the sum of the direct and diffuse radiation over all directions defined by vectors $\Omega =\left\{\mu ,\phi \right\}=\left\{\cos \phi \sin \theta ,\sin \phi \sin \theta ,\cos \theta \right\}$, ($\phi \in \left[0,2\pi \right]$, $\theta \in \left[0,\pi \right]$).

The spatial pattern of direct radiation within the plant canopy depends on the probability that a sunbeam incident on the upper canopy along the direction Ω₀ reaches some point r within the vegetation without being reflected or scattered by vegetation elements [7,30]. The intensity of direct solar radiation within a plant canopy can be expressed as $Q_{0,\lambda }\left(r\right)=T_{m,\lambda }\exp \left(-{\displaystyle \underset{0}{\overset{l_{r,{\Omega }_0}}{\int }}\sigma \left(r-s{\Omega }_0,{\Omega }_0\right)ds}\right)$.

Here, $T_{m,\lambda }$(W/m²) is the intensity of direct solar radiation at the upper boundary of the vegetation canopy, $l_{r,{\Omega }_0}$ is the distance between the point r and the upper boundary of vegetation cover along the direction Ω₀ to the sun, and $\sigma \left(r,\Omega \right)$ is the total cross-section of the interaction (scattering and absorption) of solar radiation with vegetation elements (see Appendix A).

The intensity of diffuse radiation at some point r within the vegetation over the direction defined by vector Ω is determined by the equation:

$$\Omega \cdot \nabla I_{\lambda ,d}\left(r,\Omega \right)+\sigma \left(r,\Omega \right)I_{\lambda ,d}\left(r,\Omega \right)={\displaystyle \underset{4\pi }{\int }{I}_{\lambda ,d}\left(r,{\Omega }^{\prime }\right){\sigma }_{s\lambda }\left(r,{\Omega }^{\prime }\to \Omega \right)\mathrm{d}{\Omega }^{\prime }}+{\sigma }_{s\lambda }\left(r,{\Omega }_0\to \Omega \right)Q_{0,\lambda }\left(r\right).$$

(1)

Here ${\sigma }_{s\lambda }\left(r,{\Omega }^{\prime }\to \Omega \right)$ is the differential cross-section for the scattering of the sun rays falling in the direction and scattered in the solid angle dΩ for corresponding r and Ω (see Appendix A).

The additional boundary condition for solar radiation above the vegetation canopy for Equation (1) can be written as $I_{d,\lambda }\left(r_t,\Omega \right)=T_{d,\lambda }/2\pi$, where $T_{d,\lambda }$(W/m²) is the intensity of diffuse radiation incoming into the point $r_t=\left\{x,y,z_t\right\}$ at the upper boundary of the vegetation canopy.

The upward radiation in the point $r_b=\left\{x,y,z_b\right\}$ at the soil surface below tree crowns can be expressed as the sum of direct and diffuse solar radiation that is reflected off the soil surface. Direct radiation reaches the soil surface without interacting with vegetation elements, whereas diffuse solar radiation is scattered by vegetation elements before reflecting off the soil surface:

$$I_{d,\lambda }\left(r_b\right)={\rho }_{\lambda }{\displaystyle \underset{\left({\Omega }^{\prime }{n}_b\right)>0}{\int }{I}_{d,\lambda }\left(r_b,{\Omega }^{\prime }\right)\left|\left({\Omega }^{\prime },n_b\right)\right|d{\Omega }^{\prime }}+{\rho }_{\lambda }\left|\left({\Omega }_0,n_b\right)\right|Q_{0,\lambda }\left(r_b\right)$$

Here, $n_b$ is outward normal to the soil surface and ${\rho }_{\lambda }$ is the soil surface reflection coefficient. We assume that the ground surface is horizontally uniform and that reflected radiation is uniformly distributed in all directions of the upper hemisphere. The equation for total reflected radiation at soil surface $I_{d,\lambda }\left(r_b,\Omega \right)$ can be therefore written as:

$$I_{d,\lambda }\left(r_b,\Omega \right)=\frac{1}{2\pi }{\rho }_{\lambda }\left({\displaystyle \underset{\left({\Omega }^{\prime }{n}_b\right)>0}{\int }{I}_{d,\lambda }\left(r_b,{\Omega }^{\prime }\right)\left|\left({\Omega }^{\prime },n_b\right)\right|d{\Omega }^{\prime }}+\cos {\theta }_0\cdot Q_{0,\lambda }\left(r_b\right)\right),\left({\Omega }^{\prime }{n}_b\right)<0$$

We also assume that the forest plot under consideration is surrounded by a forest stand with similar leaf area indexes (LAI) and characterized by the same species composition. These assumptions allow us to parameterize solar radiation at the lateral boundary as ${I_{d,\lambda }\left(r_l,\Omega \right)|}_{\left(\Omega ,n_l\right)<0}={I_{d,\lambda }\left(r_l,\Omega \right)|}_{\left(\Omega ,n_l\right)>0}$. Here, $r_l=\left\{x_l,y_l,z\right\}$ is a point on the lateral boundary, and $n_l$ is its outward normal.

The reflected radiation is calculated in accordance with [23]. Let us assume that $I\left(r_M,\Omega \right)$ is the intensity of upward radiation reflected from some point $r_M=\left\{x_M,y_M,z_0\right\}$ at the upper boundary of the forest stand in the direction $\Omega =\left\{\xi ,\eta ,\mu \right\}$ ($\mu =\cos \theta$, θ is the solar zenith angle) (Figure 1). The intensity of the reflected radiation for some point at a reference height h placed above the point $r=\left\{x,y,z_0\right\}$ can be estimated as a total intensity of reflected radiation coming out from all points with coordinates $r_M=\left\{x_M,y_M,z_0\right\}$. Coordinates of vector $r_M$ can be expressed in general form as $r_M=r\{x,y,z_0\}-l$, where l is calculated as

$$l=\left\{h\mathrm{tg}\theta \frac{\xi }{\sqrt{{\xi }^2+{\eta }^2}},h\mathrm{tg}\theta \frac{\eta }{\sqrt{{\xi }^2+{\eta }^2}},z_0\right\}=\left\{\frac{h}{\mu }\sqrt{1-{\mu }^2}\frac{\xi }{\sqrt{1-{\mu }^2}},\frac{h}{\mu }\sqrt{1-{\mu }^2}\frac{\eta }{\sqrt{1-{\mu }^2}},z_0\right\}=\left\{\frac{h}{\mu }\xi ,\frac{h}{\mu }\eta ,z_0\right\}$$

2.2. Scenarios of Numerical Experiments

To derive the possible effects of multispecies forest structure on PAR albedo (the ratio of PAR irradiance reflected to the irradiance received by a surface) and angle distribution of reflected PAR, we considered a mixed forest stand consisting of deciduous (Silver birch; Betula pendula Roth.) and coniferous (Scots pine; Pinus sylvestris L.) tree species. Such species composition is very typical for the Southern-European taiga.

All trees of different species were randomly distributed within the 100 m × 100 m modeling domain (Figure 2). The proportion of each tree species based on aboveground biomass within the modeling domain of our numerical experiments was assumed to vary between 0 and 1. Thus, the model scenarios can imitate both monospecific (pine or birch) and mixed (pine and birch trees in different proportions) forest stands. All numerical experiments were conducted for two model scenarios imitating forest stands with high (LAI = 4 m² m⁻²) and low (LAI = 1 m² m⁻²) stock densities. We assumed tree species have the same height (10 m) but different vertical biomass distributions.

The spatial distributions of the trees of different species were simulated using a random number generator. It was assumed that the forest stock density in case of LAI = 4 m² m⁻² was 600 trees per hectare, and in case of LAI = 1 m² m⁻² was 155 trees per hectare. To specify the locations of the trees within the modeling domain, we generated in the first step 4 arrays of 400 random numbers via a built-in compiler function. We then took half of the numbers from each array as x coordinates of a point and the other half as y coordinates. This procedure defined 800 points. After that, we generated another 800 points by interchanging the x and y coordinates. We thus obtained 1600 points for possible tree locations. As 1600 trees is too large an amount for 1 hectare we excluded some trees from our consideration. For that we generated an array of 1600 markers. The value of 600 markers (for LAI = 4 m² m⁻²) related to the points of tree location were set to 1, while other markers were set to 0. Tree location points were chosen on the principle that the distance between trunks was greater than 2 m. Finally, we integrated the tree biomass distribution functions just for the tree trunk locations at points with numbers corresponding to markers equaling to 1.

The spatial distribution of biomass for individual trees was parameterized as (Figure 3):

$$F_{tree}\left(R,z\right)=\{\begin{array}{cc}\begin{array}{c}\max \left(F_{leaves}\left(R,z\right),k_B{F}_{branches}\left(R,z\right)\right),\\ k_B,\\ 0,\\ 0,\end{array}& \begin{array}{c}{R}_T\le R\le R_C\left(z\right),L_T\le z\le z_0;\\ R\le R_T,0\le z\le L_T;\\ R>R_C\left(z\right);\\ z>z_0,\end{array}\end{array}$$

where $R=\sqrt{{\left(x-x_T\right)}^2+{\left(y-y_T\right)}^2}$, $\left(x_T,y_T\right)$ are the coordinates of each tree trunk location, z is the height above the ground, $z_0=10$ m is the tree height, $R_T=0.05$ m is the tree trunk radius similar for all considered tree species, L_T is the height of crown base counting from the ground, factor $k_B=2/3$ shows the difference in the reflection and transmission coefficients of leaves and branches for PAR, $R_C\left(z\right)$ is a function specifying the crown radius at height z, F_leaves(R,z), F_branches(R,z) are the function of leaves and branches distribution within the tree crown, respectively. The expressions for these functions and the values of parameters for both tree species are detailed in Appendix B.

The optical properties of the leaves for Silver birch and Scots pine were taken from a pilot study [32]: the leaf (needle) reflection coefficient of Scots pine was 0.07, and the transmission coefficient was 0.04. The corresponding values for reflection and transmission coefficients of the Silver birch leaves were both 0.08.

The PAR albedo and the angle distribution of reflected radiation for forest stands under consideration was simulated under sunny weather conditions for different sun elevations (h₀) that can be observed in the mid latitudes in the northern hemisphere in summer.

The ratio between direct and total solar radiation was derived as a function of sun elevation (h₀) using results of long-term solar radiation measurements conducted at the meteorological station of Moscow State University in Moscow, Russia [33]. These measurements showed the ratio between direct (S) and diffuse (D) solar radiation can be very well described as a function of sin(h₀) using a second order polynomial function (r² = 0.99, p < 0.01): $S/D=0.159\cdot {\sin }^2(h_0)+3.48\cdot \sin (h_0)-0.122$.

To derive the effects of soil optical properties on reflected and absorbed PAR, all modeling experiments were conducted for both dark (ρ_s = 0.01) and light (ρ_s = 0.20) soils.

We also assumed in our numerical experiments that the atmospheric absorption term is negligible [22].

2.3. The Numerical Scheme

The numerical scheme suggested by Myneni et al. [34] was used to solve the integro-differential Equation (1). The equation was solved using discrete ordinates method. It is based on splitting the angular variable Ω into a smaller number of directions. Here, the discretization was performed using the Gauss–Markov quadrature formulae of the 17th order [35]. Equation (1) was solved for each of the discrete ordinates ${\Omega }_m,m=0\dots 109$ via the iteration procedure [34].

We selected the following boundaries for our modeling domain considering the horizontal and vertical dimensions of the selected forest plot (100 m × 100 m, tree height—10 m): $-54\mathrm{m}\le x,y\le 54\mathrm{m},0\mathrm{m}\le z\le 11.5\mathrm{m}$. The choice of domain dimensions was governed by the need to step back from the edges of the hectare domain under consideration to reduce the influence of the boundary conditions. The choice of domain was also constrained by computer resources. The entire modeling domain was also divided into uniform grid cells: $\left(x_i,y_j,z_k\right)$, $i=\overline{0,N_x-1}$, $j=\overline{0,N_y-1}$, $k=\overline{0,N_z-1}$ each of which was characterized by a specific density of photosynthesizing and non-photosynthesizing phyto-elements (branches, stem) for each tree species. The OZ-axes was directed from the top of our modeling domain to its bottom.

We use the notation $I{x}_{i-\frac{1}{2},j,k}^{\left(n+1\right)}\left({\Omega }_m\right)$, $i=\overline{0,N_x}$, $j=\overline{0,N_y-1}$ and $k=\overline{0,N_z-1}$ to specify the radiance spreading along the direction Ω_m at the face of (i,j,k) grid cell perpendicular to OX axis at (n + 1)th iteration. Similarly, notations $I{y}_{i,j-\frac{1}{2},k}^{\left(n+1\right)}\left({\Omega }_m\right)$, $i=\overline{0,N_x-1}$, $j=\overline{0,N_y}$ and $k=\overline{0,N_z-1}$ for radiance at the grid cell face perpendicular to OY axis; $I{z}_{i,j,k-\frac{1}{2}}^{\left(n+1\right)}\left({\Omega }_m\right)$, $i=\overline{0,N_x-1}$, $j=\overline{0,N_y-1}$, and $k=\overline{0,N_z}$ at the cell face perpendicular to OZ axis; and $I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)$ for the central-cell value.

Equation (1) can be rewritten in finite difference form as:

$$\begin{array}{l}-\frac{{\mu }_m}{\Delta z}\left(I{z}_{i,j,k+\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)-I{z}_{i,j,k-\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)\right)+\frac{{\eta }_m}{\Delta y}\left(I{y}_{i,j+\frac{1}{2},k}^{(n+1)}\left({\Omega }_m\right)-I{y}_{i,j-\frac{1}{2},k}^{(n+1)}\left({\Omega }_m\right)\right)+\\ +\frac{{\xi }_m}{\Delta x}\left(I{x}_{i+\frac{1}{2},j,k}^{(n+1)}\left({\Omega }_m\right)-I{x}_{i,j-\frac{1}{2},k}^{(n+1)}\left({\Omega }_m\right)\right)+\sigma \left(r_{i,j,k},{\Omega }_m\right)I_{i,j,k}^{(n+1)}=J_{i,j,k}^{(n)}\left({\Omega }_m\right),\end{array}$$

(2)

where $r_{i,j,k}=\left\{x_i,y_j,z_k\right\}$,

$$J_{i,j,k}^{(n)}\left({\Omega }_m\right)={\displaystyle \underset{4\pi }{\int }{I}_{i,j,k}^{(n)}\left(r_{i,j,k},{\Omega }^{\prime }\right){\sigma }_{s\lambda }\left(r_{i,j,k},{\Omega }^{\prime }\to {\Omega }_m\right)\mathrm{d}{\Omega }^{\prime }}+{\sigma }_{s\lambda }\left(r_{i,j,k},{\Omega }_0\to {\Omega }_m\right)Q_{0,\lambda }\left(r_{i,j,k}\right)$$

The value of the integral on the right-hand side was obtained numerically as:

$${\displaystyle \underset{4\pi }{\int }{I}_{i,j,k}^{(n)}\left(r_{i,j,k},{\Omega }^{\prime }\right){\sigma }_{s\lambda }\left(r_{i,j,k},{\Omega }^{\prime }\to {\Omega }_m\right)\mathrm{d}{\Omega }^{\prime }}=4\pi {\displaystyle \sum _{m^{\prime }=0}^{109}{I}_{i,j,k}^{(n)}\left(r_{i,j,k},{\Omega }_{m^{\prime }}\right){\sigma }_{s\lambda }\left(r_{i,j,k},{\Omega }_{m^{\prime }}\to {\Omega }_m\right)W\left({\Omega }_{m^{\prime }}\right)}$$

with discrete ordinates Ω_m and weights $W\left({\Omega }_m\right)$ taken according to Gauss–Markov quadrature formula of 17th order [35].

In order to link unknown variables in (2) we used the equalities called the diamond difference relations (DD method):

$$\begin{array}{l}{I}_{i,j,k}^{(n+1)}\left({\Omega }_m\right)=\frac{1}{2}\left(I{x}_{i-\frac{1}{2},j,k}^{\left(n+1\right)}\left({\Omega }_m\right)+I{x}_{i+\frac{1}{2},j,k}^{\left(n+1\right)}\left({\Omega }_m\right)\right),I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)=\frac{1}{2}\left(I{y}_{i,j-\frac{1}{2},k}^{\left(n+1\right)}\left({\Omega }_m\right)+I{y}_{i,j+\frac{1}{2},k}^{\left(n+1\right)}\left({\Omega }_m\right)\right),\\ I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)=\frac{1}{2}\left(I{z}_{i,j,k-\frac{1}{2}}^{\left(n+1\right)}\left({\Omega }_m\right)+I{z}_{i,j,k+\frac{1}{2}}^{\left(n+1\right)}\left({\Omega }_m\right)\right).\end{array}$$

(3)

Excluding $I{z}_{i,j,k+\frac{1}{2}}^{(n+1)}$, $I{y}_{i,j+\frac{1}{2},k}^{(n+1)}$, and $I{x}_{i+\frac{1}{2},j,k}^{\left(n+1\right)}$ from (2) and (3) we came to the following expression:

$$I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)=\frac{J_{i,j,k}^{(n)}\left({\Omega }_m\right)+2\frac{\left|{\mu }_m\right|}{\Delta z}I{z}_{i,j,k-\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)+2\frac{\left|{\eta }_m\right|}{\Delta x}I{y}_{i,j-\frac{1}{2},k}^{(n+1)}\left({\Omega }_m\right)+2\frac{\left|{\xi }_m\right|}{\Delta x}I{x}_{i-\frac{1}{2},j,k}^{(n+1)}\left({\Omega }_m\right)}{\sigma \left(r_{i,j,k},{\Omega }_m\right)+2\frac{\left|{\mu }_m\right|}{\Delta z}+2\frac{\left|{\xi }_m\right|}{\Delta y}+2\frac{\left|{\eta }_m\right|}{\Delta x}}.$$

(4)

Then we completed the following steps:

(1) In the first step we determined the radiance spreading along the directions ${\Omega }_m=\left\{{\xi }_m,{\eta }_m,{\mu }_m\right\},$ ${\xi }_m^2+{\eta }_m^2+{\mu }_m^2=1$ with ${\xi }_m>0$, ${\eta }_m\ge 0$, ${\mu }_m\le 0$ (the first octant). In this case the values $I{x}_{-\frac{1}{2},0,0}^{(n+1)}\left({\Omega }_m\right)$, $I{y}_{0,-\frac{1}{2},0}^{(n+1)}\left({\Omega }_m\right)$, $I{z}_{0,0,-\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)$ were known from the boundary conditions. Then the value of the radiance for $i=j=k=0$ was obtained from relation (4).

(2) Taking into account the known value $I_{0,0,0}^{(n+1)}\left({\Omega }_m\right)$ from (3) we expressed $I{z}_{0,0,\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)$, $I{y}_{0,\frac{1}{2},0}^{(n+1)}\left({\Omega }_m\right)$, $I{x}_{\frac{1}{2},0,0}^{\left(n+1\right)}\left({\Omega }_m\right)$—the values for all indices i, j, and k fulfilling the condition i + j + k = 1.

(3) Now the central-cell values $I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)$ for the indices with $i+j+k=1$ can be obtained from (4). And so on: taking into account the known value $I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)$ from (3) we expressed $I{z}_{i,j,k+\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)$, $I{y}_{i,j+\frac{1}{2},k}^{(n+1)}\left({\Omega }_m\right)$ and $I{x}_{i+\frac{1}{2},j,k}^{\left(n+1\right)}\left({\Omega }_m\right)$ for all indices i, j, and k satisfying the relation $i+j+k=c$, where c consistently takes the values $1,2,\dots ,N_x+N_y+N_z-3$.

The calculation procedure is schematically shown in Figure 4a.

(5) The same procedure was used for calculating the second octant ordinates: ${\xi }_m\le 0,{\eta }_m>0,{\mu }_m\le 0$ with the values $I{x}_{N_x-\frac{1}{2},0,0}^{(n+1)}\left({\Omega }_m\right)$, $I{y}_{N_x-1,-\frac{1}{2},0}^{(n+1)}\left({\Omega }_m\right)$, $I{z}_{N_x-1,0,-\frac{1}{2}}^{(n+1)}\left({\Omega }_m\right)$ taken from the boundary conditions. We obtain $I_{i,j,k}^{(n+1)}\left({\Omega }_m\right)$ from (4) at each step and then express $I{z}_{i,j,k+\frac{1}{2}}^{(n+1)}$, $I{y}_{i,j+\frac{1}{2},k}^{(n+1)}$, $I{x}_{i+\frac{1}{2},j,k}^{\left(n+1\right)}$ from (3) for all indices i, j, k fulfilling the conditions $N_x-1-i+j+k=c$, where $c=1,2,\dots ,N_x+N_y+N_z-3$. (Figure 4b).

(6) The calculations were similarly carried out for other octants (see Figure 4c–h).

(7) The flux of solar radiation at each point $\left(x_i,y_j,z_k\right)$ was calculated as the total intensity of scattered radiation in each direction ${\Omega }_m$ and summed with the intensity of direct solar radiation.

2.4. The Inverse Problem Statement to Retrieve the Forest Species Composition

Let us assume that $\left(x,y\right)\in D$ is the modeling domain with randomly distributed trees of different species. The functions ${\rho }_1\left(x,y\right)$ and ${\rho }_2\left(x,y\right)$ determine the spatial distributions of different tree species (trees per hectare) (here the subscripts “1” and “2” indicate deciduous and coniferous tree species, respectively) within the modeling domain. Let us also assume that we have known functions $I_1\left(\phi ,\theta \right)$ and $I_2\left(\phi ,\theta \right)$ describing the radiation intensity that is reflected from some average tree of each species in the direction determined by spherical coordinates $\left(\phi ,\theta \right)$.

Thus, for some very small pixel ds specified by coordinates $\left(x,y\right)$ on which both types of trees species are grown, we can calculate the radiation intensity that is measured by a PAR sensor located at some point with coordinates $\left(x_s,y_s,z_s\right)$ along the corresponding solid angle dΩ under which the surface element ds is visible: $d\Omega \equiv d\Omega (x,y,x_s,y_s,z_s,ds)\equiv k(x,y,x_s,y_s,z_s)ds$.

The total radiation intensity I reflected from the surface element ds is a function of angle coordinates $\phi (x,y,x_s,y_s)\equiv \arctan \frac{y_s-y}{x_s-x}$ and $\theta (x,y,x_s,y_s,z_s)\equiv \arctan \frac{\sqrt{{\left(x_s-x\right)}^2+{\left(y_s-y\right)}^2}}{z_s}$:

$$\begin{array}{l}I(\phi (x,y,x_s,y_s),\theta (x,y,x_s,y_s,z_s))d\Omega (x,y,x_s,y_s,z_s,ds)=\\ =I_1(\phi (x,y,x_s,y_s),\theta (x,y,x_s,y_s,z_s)){\rho }_1(x,y)ds+I_2(\phi (x,y,x_s,y_s),\theta (x,y,x_s,y_s,z_s)){\rho }_2(x,y)ds.\end{array}$$

In coordinates (x,y,z) this relation can be rewritten as:

$$I_1(x,y,x_s,y_s,z_s){\rho }_1(x,y)ds+I_2(x,y,x_s,y_s,z_s){\rho }_2(x,y)ds=I(x,y,x_s,y_s,z_s)k(x,y,x_s,y_s,z_s)ds.$$

(5)

The inverse problem can be posed as follows: it is necessary to determine the functions ${\rho }_1(x,y),{\rho }_2(x,y)\in W_2^2(D)$ from measured values of the function $I(\phi ,\theta ,x_s,y_s,z_s)\in L_2(D)$ along some trajectory above a forest canopy ($x_s=x_s\left(p\right)$, $y_s=y_s\left(p\right)$, and $z_s=z_s\left(p\right)$, $p\in \left[p_1,p_2\right]$).

2.5. The Model Algorithm for Inverse Problem Solving

The model algorithm can be divided into seven sequential steps:

1. In the first step we define the domain D on which model functions ${\rho }_1\left(x,y\right)$ and ${\rho }_2\left(x,y\right)$ are retrieved as a rectangle with following dimensions:$D\equiv \left\{\left(x,y\right):L_x\le x\le R_x,L_y\le y\le R_y\right\}$.

2. We divide region D into $N_x\times N_y$ grid cells, at the center of each we put the node. As a result, we obtain a $X_{N_x}\times Y_{N_y}$ mesh with a step $h_x=\left(R_x-L_x\right)/N_x$ with respect to the variable x and the step $h_y=\left(R_y-L_y\right)/N_y$ with respect to y:

$$X_{N_x}\times Y_{N_y}\equiv \left\{\begin{array}{l}\left(x_{n_x},y_{n_y}\right),1\le n_x\le N_x,1\le n_y\le N_y:x_{n_x}=L_x+\frac{h_x}{2}+\left(n_x-1\right)h_x,\\ y_{n_y}=L_y+\frac{h_y}{2}+\left(n_y-1\right)h_y\end{array}\right\}$$

3. We introduce a uniform grid on the part of the trajectory along which the measurements of functions $I_1$ and $I_2$ were made. For simplicity, we assume that $p\in \left[0,1\right],$ $x_s=x_s\left(p\right)\equiv L_s+\left(R_s-L_s\right)p,$ $y_s=y_s\left(p\right)\equiv l_y=const,$ $z_s=z_s\left(p\right)\equiv l_z=const.$ Thus, the corresponding mesh is uniform and co-oriented in space parallel to the coordinate axis Ox. The related mesh $S_{N_p}$ with $N_p-1$ intervals of length $h_s=\left(R_s-L_s\right)/\left(N_p-1\right)$ takes the form:$S_{N_p}\equiv \left\{s_{n_p},1\le n_p\le N_p:s_{n_p}=L_s+\left(n_p-1\right)h_s\right\}$.

4. We calculate the values of the functions $I_1(\phi ,\theta )$ and $I_2(\phi ,\theta )$:

$$\begin{array}{l}{I}_1(\phi (x_{n_x},y_{n_y},s_{n_p},l_y),\theta (x_{n_x},y_{n_y},s_{n_p},l_y,l_z))\equiv I_1{}_{n_x,n_y,n_p},\\ I_2(\phi (x_{n_x},y_{n_y},s_{n_p},l_y),\theta (x_{n_x},y_{n_y},s_{n_p},l_y,l_z))\equiv I_2{}_{n_x,n_y,n_p},\end{array}$$

where $\phi (x_{n_x},y_{n_y},s_{n_p},l_y)=\{\begin{array}{ll}\arctan \frac{l_y-y_{n_y}}{s_{n_p}-x_{n_x}},\hfill & \mathrm{if}{l}_y-y_{n_y}\ge 0\mathrm{and}{s}_{n_p}-x_{n_x}\ge 0,\hfill \\ \arctan \frac{l_y-y_{n_y}}{s_{n_p}-x_{n_x}}+2\pi ,\hfill & \mathrm{if}{l}_y-y_{n_y}<0\mathrm{and}{s}_{n_p}-x_{n_x}\ge 0,\hfill \\ \arctan \frac{l_y-y_{n_y}}{s_{n_p}-x_{n_x}}+\pi ,\hfill & \mathrm{if}{s}_{n_p}-x_{n_x}<0,\hfill \end{array}$

$\theta (x_{n_x},y_{n_y},s_{n_p},l_y,l_z)=\arctan \frac{\sqrt{{\left(s_{n_p}-x_{n_x}\right)}^2+{\left(l_y-y_{n_y}\right)}^2}}{l_z}$ for all $n_x=\overline{1,N_x}$, $n_y=\overline{1,N_y}$, $n_p=\overline{0,N_p}$.

5. We calculate the values of model functions ${\rho }_1(x,y)$ and ${\rho }_2(x,y):$ ${\rho }_1(x_{n_x},y_{n_y})\equiv {\rho }_1{}_{n_x,n_y},{\rho }_2(x_{n_x},y_{n_y})\equiv {\rho }_2{}_{n_x,n_y},$ for all $n_x=\overline{1,N_x}$, $n_y=\overline{1,N_y}$.

6. The main equation (5) can be rewritten as:

$$I_1{}_{n_x,n_y,n_p}{\rho }_1{}_{n_x,n_y}{h}_x{h}_y+I_2{}_{n_x,n_y,n_p}{\rho }_2{}_{n_x,n_y}{h}_x{h}_y={\tilde{I}}_{n_x,n_y,n_p},$$

and for all $n_x=\overline{1,N_x}$, $n_y=\overline{1,N_y}$, $n_p=\overline{0,N_p}$ that provides the system of linear algebraic equations (SLAE) in the form

$$AX=Y$$

(6)

with matrix of operator A having a dimension $N_x{N}_y{N}_p\times 2N_x{N}_y$ and vector X of length $2N_x{N}_y$, with non-zero components as follows:

$$\begin{array}{r}{A}_{(n_p-1)N_x{N}_y+(n_x-1)N_y+(n_y-1)+1,(n_x-1)N_y+(n_y-1)+1}=I_1{}_{n_x,n_y,n_p}{h}_x{h}_y,\\ \mathrm{for}{n}_x=\overline{1,N_x},n_y=\overline{1,N_y}\mathrm{and}{n}_p=\overline{1,N_p},\\ A_{(n_p-1)N_x{N}_y+(n_x-1)N_y+(n_y-1)+1,N_x{N}_y+(n_x-1)N_y+(n_y-1)+1}=I_2{}_{n_x,n_y,n_p}{h}_x{h}_y,\\ \mathrm{for}{n}_x=\overline{1,N_x},n_y=\overline{1,N_y}\mathrm{and}{n}_p=\overline{1,N_p},\end{array}$$

$$X_n=\{\begin{array}{ll}{\rho }_1{}_{[(n-1)/N_y]+1,\{(n-1)/N_y\}+1},\hfill & \mathrm{if}n=\overline{1,N_x{N}_y},\hfill \\ {\rho }_2{}_{[(n-1-N_x{N}_y)/N_y]+1,\{(n-1-N_x{N}_y)/N_y\}+1},\hfill & \mathrm{if}n=\overline{N_x{N}_y+1,2N_x{N}_y}.\hfill \end{array}$$

Square brackets mean a whole part of the expression and figured brackets mean the decimal part.

As a result, via (6) using the model values of functions ${\rho }_1(x,y)$ and ${\rho }_2(x,y)$ on the introduced grid $X_{N_x}\times Y_{N_y}$ as the components of vector $X_{model}$, we simulate a model vector Y of input data having the length $N_x{N}_y{N}_p$ according to the following algorithm:

(a) Error simulation of operator A elements (that means a model error):

$$A_{hm,n}=A_{m,n}\left(1+\frac{\xi h}{100\%}\right)\mathrm{for}\mathrm{all}m=\overline{1,N_x{N}_y{N}_p},n=\overline{1,2N_x{N}_y}$$

where ξ is a uniformly distributed on the segment [−1,1] random value, h is a given error in percent.

(b) Calculation of vector Y from the ratio:

$$A_h{X}_{model}=Y$$

(c) Error simulation of Y components (that means the error of experimental data):

$$Y_{\delta }{}_m=Y_m\left(1+\frac{\xi \delta }{100\%}\right)\mathrm{for}m=\overline{1,N_x{N}_y{N}_p},$$

where ξ is a uniformly distributed random value on the segment [-1,1], δ is a given error in percent.

7. We seek for the solution of (6) with simulated vector $Y_{\delta }$ by the least square method:

$$X={(A^{T}A)}^{-1}{A}^T{Y}_{\delta }$$

If the solution of this equation is unstable, the regularizing algorithm can be used. It can be done, for instance, by using the Tikhonov regularizing algorithm based on the minimization of smoothing functional [36]. In this case the solution of (6) can be found as a minimum of Tikhonov functional $M^{\alpha }[x]$ (smoothing functional):

$$M^{\alpha }[X]=\Vert AX-Y_{\delta }{\Vert }^2+\alpha \Vert RX{\Vert }^2,$$

where α is the regularization parameter, R is the “smoothing operator” (for simplicity while solving the model problem, we shall use R = E, where E is a unit operator). The extremal of this functional can be found as the solution of SLAE:

$$(A^{T}A+\alpha R^{T}R)X=A^T{Y}_{\delta }-1A^T{Y}_{\delta }$$

(7)

In case the solution of (7) is stable with respect to the model errors (h) and the errors of input data (δ) we have α = 0. Otherwise if the solution of (7) is unstable with respect to the mentioned errors, then the choice of the regularization parameter α must be consistent with the errors δ of the input data Y_δ and h of the operator A. Such a choice can be made using, e.g., the generalized discrepancy principle [36]:

$$\rho (\alpha )=0,\mathrm{where}\rho (\alpha )=\Vert A{X}^{\alpha }-Y_{\delta }{\Vert }^2-{\left(\Delta +H\Vert X^{\alpha }\Vert \right)}^2-{\mu }^2$$

Here $X^{\alpha }$ is the extremal of the functional $M^{\alpha }\left[X\right]$ ($M^{\alpha }\left[X^{\alpha }\right]=\inf M^{\alpha }\left[X\right]$); $\Delta >0$ is the root-mean-square error of the specification ($Y_{\delta }$) in right-hand side of (6):

$${\Delta }^2={\displaystyle \sum _{m=1}^{N_x{N}_y{N}_p}{\left(Y_m\frac{\xi \delta }{100\%}\right)}^2};$$

$H>0$ is the root-mean-square error of operator A specification:

$$H^2={\displaystyle \sum _{m=1}^{N_x{N}_y{N}_p}{\displaystyle \sum _{n=1}^{2N_x{N}_y}{\left(A_{m,n}\frac{\xi h}{100\%}\right)}^2}};$$

where μ is the measure of incompatibility of SLAE (6), i.e., $\mu =\inf \Vert A_{h}X-Y_{\delta }\Vert =\{\Vert A_h{X}^{\alpha }-Y_{\delta }\Vert \alpha =0\}$.

2.6. Modeling Design for Inverse Problem Solution

To retrieve the forest species composition within our modeling domain D (Figure 2) with spatial dimensions (in meters) $L_x=-50$, $R_x=50$, $L_y=-50$, $R_y=50$, let us select the segment of the straight line defined by the parameters $L_s=-60$, $R_s=60$, $l_y=10$, $l_z=200$. This could correspond to the flight line of a helicopter, aircraft or drone carrying a PAR sensor 200 m above the forest (Figure 5).

The distribution of each tree species in a mixed forest stand is usually discrete. In our case, in order to test the suggested algorithm, we assumed the smoothed discrete distributions of different tree species within the modeling domain and described them by functions ${\rho }_1\left(x,y\right)$ and ${\rho }_2\left(x,y\right)$ specifying the density (trees per hectare) of each type of tree species (Figure 6).

The values of functions $I_1\left(\phi ,\theta \right)$ and $I_2\left(\phi ,\theta \right)$ where $I_{1,2}\left(\phi ,\theta \right)\equiv I_{1,2}(\phi (x_{n_x},y_{n_y},s_{n_p},l_y),\theta (x_{n_x},y_{n_y},s_{n_p},l_y,l_z))\equiv I_{1,2}{}_{n_x,n_y,n_p}$ for each sensor position $(s_{n_p},l_y,l_z)$ along the selected trajectory and for each node of the mesh $\left(n_x,n_y\right)$ were determined as radiation scattered by single trees along the direction $\left(\phi ,\theta \right)$, $0\le \phi <360°$, $0\le \theta <90°$ into the upper hemisphere. The reflected radiation was calculated using the 3D model, based on the assumption that the sun elevation was 30° ($h_0=30°$) and the sun rays were incident from the south.

3. Results and Discussion

3.1. Reflection of PAR for Mixed Forest Stand

Numerous experimental and modeling studies during recent decades showed that PAR albedo measured by different remote sensing equipment or by devices installed at meteorological towers can be broadly used to estimate canopy structure and architecture, surface net radiation, gross and net primary production, etc. [25,37,38,39,40]. Let us consider the possible effects of the change of forest species composition on variation of reflected PAR.

The numerical experiments provided using the 3D radiative transfer model showed a clear dependence of PAR albedo of both monospecific and mixed forest stands on canopy density (LAI), sun elevation, and soil reflection properties (Figure 7). The canopy albedo for sparse forest canopy and dark soils reached maximal values under low sun elevation, and minimal values under high sun elevation. This is because sunbeams travel longer distances through the forest canopy before interacting with the soil surface. This results in higher scattered direct solar radiation within the plant canopy; therefore, there is a higher vegetation contribution to the total forest reflection. The overall canopy reflection also increased with a higher proportion of birch trees in the forest canopy. This is mainly due to higher birch leaf reflection and transmission coefficients compared with pine. In the case of light soil, the dependence of forest albedo on sun elevation has the opposite effect. Simulated local minimums of forest albedo at h₀ = 20° can be explained by a reduced contribution of reflected direct solar radiation from light soils to the total forest albedo, joint effects of canopy architecture, and the optical properties of the soil and plant phyto-elements (leaves, branches, stem).

In the case of a dense forest canopy (LAI = 4.0 m² m⁻²), the albedo gradually increases with increasing sun elevation. The rate of this increase is slightly higher for dark soils. The clear dependence of albedo on species composition mainly manifested under high sun elevations. The forest canopy with a higher proportion of birch trees was characterized by higher albedo. Under low sun elevation (h₀ < 20°), the albedo was less dependent on species composition. A local minimum of canopy reflection under a sun height of ~45° can be caused by the influence of forest canopy architecture (spatial leaf area distribution; orientation of leaf blades) and the optical properties of the tree species and soil surface.

Summarizing all obtained results, it is possible to emphasize a high sensitivity of PAR albedo to forest canopy density, soil properties, and sun conditions. The dependence of forest albedo on the proportion of deciduous and coniferous species in a forest stand, except at low sun elevations (h₀ < 25–30°), was also shown. The results of our numerical experiments agreed well with available experimental and modeling data obtained by different scientific teams for various plant canopy types. In particular, Knyazikhin et al. [22] and Gravenhorst et al. [23] showed a clear impact of sun elevation on the PAR albedo of monospecific old-grown spruce forest in Solling, Germany. Maximal albedo was observed in the morning and evening, and the minimal values were close to noon. The soil had a minor impact on albedo because of a very dense overstory and grassy understory that almost completely covered the soil surface. Similar results were also obtained during the field observations provided by Pinker et al. (1980) [41] in a tropical dry evergreen forest in Southeast Asia and by Giambelluca et al. (1997) [42] at several forest and shrub locations in Amazonia.

Lower albedo at higher sun elevations can be explained by deeper penetration of solar photons into the vegetation cover under a lower solar zenith angle, more frequent interaction of solar photons with plant components, and higher probability to be absorbed by plant elements and soil surface before escaping back into the atmosphere [43,44]. Another factor contributing to the change in surface albedo with solar elevation angle is the soil surface fraction illuminated by direct solar beams. At the highest sun positions (i.e., highest sun elevation angle or lowest solar zenith angle), the contribution of directly illuminated soil surface into the total reflected solar radiation reach maximum values. As solar zenith angle and shadow fraction increase, the proportion of directly illuminated soil surface decrease and the contribution of the soil surface to surface albedo also decrease.

In discussing the possible effects of forest canopy architecture on radiative transfer characteristics, it is also necessary to consider the influence of leaf inclination on light reflection. In particular, Simioni et al. (2013) [45] showed that the light reflection and absorption of mixed deciduous and coniferous forests is highly sensitive to leaf blade inclination of different tree species, and erroneous assumptions about leaf inclination lead to larger uncertainties when modeling heterogeneous mixed canopies. In our numerical experiments we assumed that both our model trees (birch and pine) have the same spherical angular orientation of phyto-elements. Thus, possible uncertainties in estimating the canopy albedo, due to differences in angular orientation of phytoelements between different species in our numerical experiments, can be ignored.

Another factor that can highly influence the forest canopy albedo is the surface topography. Kranigk et al. (1994) [46] and Knyazikhin et al. (1996, 2005) [21,31] performed experimental and modeling studies on radiation regime in a spruce forest stand growing in the Harz Mountains in Central Germany. It was shown that complex topography can result in strong redistribution of reflected and penetrated solar radiation depending on the surface slope and aspect characteristics, and in turn, can influence forest net primary production and tree growth. In our numerical experiments, we considered the modeling scenario assuming only a flat topography. Thus, any uncertainties in radiation reflection characteristics, due to complex topography, can be ignored.

3.2. Retrieving Species Composition of a Mixed Forest Stand from Canopy Reflectance Properties

The results of numerical experiments to retrieve the density of deciduous and coniferous tree species (${\rho }_1\left(x,y\right)$ and ${\rho }_2\left(x,y\right)$) in a mixed forest stand from canopy reflection properties are shown in Figure 8, Figure 9 and Figure 10. As can be seen, with an increase of the error δ and h with which the canopy reflection measurements that can be carried out, the retrieving accuracy of the unknown functions drastically decreases. Figure 8 and Figure 9 represent examples of species composition retrieved for measurement error levels of $\left(\delta ;h\right)=\left(1;1\right)\%$ and $\left(\delta ;h\right)=\left(5;6\right)\%$ (the last is consistent with the model error mentioned in [23]) with a similar number of measurement points $N_p$). However, even with a sufficiently large error level in the measured signal and model error, it is possible to retrieve the unknown functions when measurements are carried out at a larger number of points. Figure 9 and Figure 10 show a comparison of the results for the same error levels of $\left(\delta ;h\right)=\left(5;6\right)\%$ in the experimental data with $N_p=10$ and $N_p=500$ measurement points.

The statements made above are supported by the following results. Figure 11 shows the dependence of the root-mean-square error ($s=\frac{\Vert X_{inv}-X_{\mathrm{mod}el}\Vert }{\sqrt{2N_x{N}_y}}$) for the components of the reconstructed vector $X_{inv}={\left(A_{}^{T}A+\alpha R^{T}R\right)}^{-1}{A}_{}^T{Y}_{\delta }$ in case $\alpha =0$ on the error levels δ of the input data and h of the model (for a fixed number of measurement points $N_p$). From this graph it can be seen that a decrease in the error δ with which the canopy reflection measurements are carried out, and the model error h, resulted in increase of the accuracies of vector $X_{inv}$ reconstruction, and hence, the grid values of functions ${\rho }_1\left(x,y\right)$ and ${\rho }_2\left(x,y\right)$. This means that $X_{inv}\to X_{model}$ when $\left(\delta ;h\right)\to 0.$ So, the suggested inverse problem statement is stable and the corresponding algorithm is therefore regularizing [36].

Figure 12 shows the dependence of the root-mean-square error in the components of the reconstructed vector $X_{inv}$ on the number of measurement points $N_p$ (for a fixed error level δ in the input data and h in the model). From this graph it is clearly seen that an increase of the input data volume (determined by the number $N_p$ of the measurement points) resulted in an increase of accuracies of vector $X_{inv}$ reconstruction and hence, the grid values of functions ${\rho }_1\left(x,y\right)$ and ${\rho }_2\left(x,y\right)$. These results support the notion of increasing the number of measuring points as much as possible.

3.3. Other Possible Ways to Reconstruct the Proportions of Different Tree Species in a Mixed Forest from Canopy Reflection Using the 3D Model

During recent decades, many algorithms have been suggested to quantify structural and optical aspects of vegetation canopies from remote sensing [25,28,45,46,47,48]. The problem of mapping forest composition in mixed forest stands using remote sensing data is still open. Several studies used hyperspectral and/or LiDAR data from LiCHy (Hyperspectral sensors, LiDAR and CCD (Charge-Coupled Device) camera) airborne systems to classify tree species composition in mixed forests [49] or high-resolution imagery for mapping forest cover types [50,51,52,53].

As shown in the previous section, the task to retrieve forest species composition from LAI, canopy, and soil reflection is not trivial, due to the multivariate relationships between solar radiation reflection, scattering coefficients, and tree species composition (different proportions of deciduous and coniferous trees in mixed forest stand) under various LAI and soil reflection. The mean canopy albedo is relatively insensitive to changes in the proportion of different species in a forest stand—especially in the case of dense forest and low sun elevations (Figure 7). However, the main conclusion that can be drawn from Section 3.2 is that the forest stand species composition can be reconstructed from information on canopy reflection properties.

Let us consider some other possibilities to solve this inverse problem using this 3D model and find effective numerical algorithms that would allow us to reconstruct forest species composition from other structural and optical properties of a plant canopy and soil (LAI, leaf reflection and transmission, soil reflection, etc.).

1. The functions ${\rho }_1(x,y),{\rho }_2(x,y)\in H(D)$ are known from measurements of intensity $I(\phi ,\theta ,x_s,y_s,z_s)\in L_2(D)$ at a single fixed point $(x_s,y_s,z_s)$.

This formulation of the problem is physically undetermined because it restores two functions ${\rho }_1$ and ${\rho }_2$ of two variables $\left(x,y\right)$ with one function I of two variables $\left(\phi ,\theta \right)$. In other words, for a fixed set of parameters $\left(x,y,x_s,y_s,z_s\right)$, we can write only one equation with two unknown values. Thus, for all possible $\left(x,y\right)\in D$ determining all viewing angles $\left(\phi ,\theta \right)$ after digitization of the corresponding domains we obtain a system of equations with twice as many unknown values as equations. Consequently, this formulation will be incorrect according to Hadamard, regarding the lack of uniqueness of the solution.

However, in case of the available and sufficiently precise a priori information about the functions ${\rho }_1$ and ${\rho }_2$ belonging to a certain space (class of functions) H, we can reformulate the problem as determining the functions ${\rho }_1(x,y),{\rho }_2(x,y)\in H(D)$ that realize minimum of the functional $F[{\rho }_1,{\rho }_2]={\displaystyle \underset{D}{\int }{(I_1{\rho }_1+I_2{\rho }_2-Ik)}^{2}ds}$ with the assumption that $H\left(D\right)$ is some compact set of functions, e.g., monotonic functions, piecewise-convex functions, functions of bounded variation, etc. Examples of analogous problem statements using a priori information on the desired functions are available in some studies [36,54,55,56].

Thus, the task to find ${\rho }_1(x,y),{\rho }_2(x,y)\in H(D)$ can be rewritten in the following form:

$$\{{\rho }_1,{\rho }_2\}=\underset{{\rho }_1,{\rho }_2\in H(D)}{\arg \min }F[{\rho }_1,{\rho }_2]\equiv \underset{{\rho }_1,{\rho }_2\in H(D)}{\arg \min }{\displaystyle \underset{D}{\int }{(I_1{\rho }_1+I_2{\rho }_2-Ik)}^{2}ds}.$$

The advantage of this approach is that (i) it is possible to find a solution to the problem having the function of radiation intensity distribution measured at a single point, and (ii) the formulation of the problem can be correct both in terms of uniqueness and in terms of stability of the solution. The disadvantages of this approach is that (i) for the case of each compact set $H\left(D\right)$ it is necessary to conduct theoretical analysis of the corresponding problem correctness with respect to the uniqueness and stability of its solution; (ii) the a priori information on the form of the set $H\left(D\right)$ may not be available.

This approach is not yet physically undetermined because we obtain basic additional information from the assumption that the recoverable functions have certain physical properties.

2. The functions ${\rho }_1(x,y),{\rho }_2(x,y)$ are known from measured values of function $I(\phi ,\theta ,x_s,y_s,z_s)\in L_2(D)$ at the two fixed points ${(x_s,y_s,z_s)}^1$ and ${(x_s,y_s,z_s)}^2$.

We have in this case a physically determined problem because it demands recovery of two functions ${\rho }_1$ and ${\rho }_2$ of two variables $\left(x,y\right)$ from two known functions of $I^1$ and $I^2$ of two variables $\left(\phi ,\theta \right)$: $I^n$ is the set of radiation intensities values obtained from point ${(x_s,y_s,z_s)}^n,n=1,2$.

Within this approach, it is no longer necessary to consider faithful a priori information about the form of the set $H\left(D\right)$. It is sufficient to assume that ${\rho }_1(x,y),{\rho }_2(x,y)\in W_2^2(D)$ and the problem to determine ${\rho }_1$ and ${\rho }_1$ can be solved in a classical way (e.g., the least square method):

$$\{{\rho }_1,{\rho }_2\}=\underset{{\rho }_1,{\rho }_2}{\arg \min }F[{\rho }_1,{\rho }_2],$$

where

$$\begin{array}{l}F[{\rho }_1,{\rho }_2]={\displaystyle \underset{D}{\int }}((I_1(x,y,x_s^1,y_s^1,z_s^1){\rho }_1(x,y)+I_2(x,y,x_s^1,y_s^1,z_s^1){\rho }_2(x,y)-\\ -I^1(x,y,x_s^1,y_s^1,z_s^1)k(x,y,x_s^1,y_s^1,z_s^1){)}^{2}ds+(I_2(x,y,x_s^2,y_s^2,z_s^2){\rho }_2(x,y)+\\ +I_2(x,y,x_s^2,y_s^2,z_s^2){\rho }_2(x,y)-I^2(x,y,x_s^2,y_s^2,z_s^2)k(x,y,x_s^2,y_s^2,z_s^2){)}^{2}ds.\end{array}$$

However, if the problem is ill-posed in terms of instability of the solution (see, for example [36]), it is necessary to build a regularizing algorithm to solve it, for example an algorithm based on Tikhonov functional minimization:

$$\{{\rho }_1,{\rho }_2\}=\underset{{\rho }_1,{\rho }_2}{\arg \min }{M}^{\alpha }[{\rho }_1,{\rho }_2],$$

where $M^{\alpha }[{\rho }_1,{\rho }_2]=F[{\rho }_1,{\rho }_2]+\alpha ({\Vert {\rho }_1\Vert }_{W_2^2}^2+{\Vert {\rho }_2\Vert }_{W_2^2}^2).$

For any regularization parameter $\alpha >0$, there is a unique extreme of Tikhonov functional realizing the minimum $M^{\alpha }[{\rho }_1,{\rho }_2]$. To select the regularization parameter, one can use the generalized discrepancy principle algorithm (see [36]). We then have $\gamma (\alpha )\equiv F[{\rho }_{{}_1}^{\alpha },{\rho }_{{}_2}^{\alpha }]-{\delta }^2=0$, where δ is an input error for the values $I^1$, $I^2$ and $\left\{{\rho }_1^{\alpha },{\rho }_2^{\alpha }\right\}$ tend to exact solution in the norm of $W_2^2$, and hence it is uniform in D when $\delta \to 0$.

3. The functions ${\rho }_1(x,y),{\rho }_2(x,y)\in W_2^2(D)$ are known from measured values of the function $I(\phi ,\theta ,x_s,y_s,z_s)\in L_2(D)$ at some surface (for example, $x_s=x_s(p_1,p_2),$ $y_s=y_s(p_1,p_2),$ $z_s=z_s(p_1,p_2),$ $\{p_1,p_2\}\in P$).

We have in this case a more over-determined problem than in all of the cases mentioned above, because it requires retrieving two functions ${\rho }_1$ and ${\rho }_2$ of two variables $\left(x,y\right)$ from one known function I of four variables $\left(\phi ,\theta ,p_1,p_2\right)$.

4. The functions ${\rho }_1(x,y),{\rho }_2(x,y)\in W_2^2(D)$ are known from measured values of the function $I(\phi ,\theta ,x_s,y_s,z_s)\in L_2(D)$ in a certain volume $\left(x_s,y_s,z_s\right)\in V$.

This retrieval problem is much more over-determined than the previous problem, because it demands the recovery of two functions ${\rho }_1$ and ${\rho }_2$ of two variables $\left(x,y\right)$ from one known function I of five variables $\left(\phi ,\theta ,x_s,y_s,z_s\right)$.

4. Conclusions

The results of numerical experiments using the 3D radiative transfer model showed that albedo of a mixed forest stand consisting of deciduous (birch) and coniferous (pine) was highly sensitivity to LAI, optical soil properties, and sun elevation. The highest sensitivity of forest albedo to species composition of the forest stand occurred under high sun elevations (higher than >30°) whereas at low solar elevations (<25–30°) this dependence was less pronounced. It was also shown that albedo of the dense forest stand (LAI = 4) was less sensitive to species composition changes than the sparser forest stand (LAI = 1). This dependence was pronounced for both light and dark soil types. These findings agreed well with results of prior experimental and modeling studies of radiative regimes in various types of forest canopies.

Our proposed algorithm provided a feasible approach to determine species composition of a mixed forest stand from PAR reflectance properties. This approach allowed us to quantify the density of at least two forest types from series the measurements from a theoretical helicopter, aircraft or drone with installed equipment for measuring reflected solar radiation. The accuracy of assessing vegetation species using this method depended on the number of points for canopy reflection measurements. An increase of the input data volume (determined by the number of the measurement points) resulted in higher accuracies for prediction of the proportions of different tree species. One limitation with this approach is that tree species with similar optical properties (reflection, transmission and absorption coefficients) cannot be readily differentiated. In this case, the method for retrieving forest species composition cannot be correctly applied.

The study introduced one way to solve the inverse problem and to quantify the proportions of different tree species in a forest stand using information on forest reflection properties, known LAI, and the optical properties of different tree species. Additional studies are needed to develop stable numerical algorithms for retrieving forest species composition. Moreover, it is crucial that each method be tested and compared with experimental field data. If successful, these algorithms might be a very useful tool for remote description of forest canopy characteristics and can be applied to various tasks of modern forest ecology and micrometeorology.