Effects of Mixture Mode on the Canopy Bidirectional Reflectance of Coniferous–Broadleaved Mixed Plantations

: One of the main initiatives for China to achieve the goal of being carbon neutral before 2060 is transforming monocultures into mixed plantations in subtropical China, because mixed forests possess a higher quality than monocultures in various ways. Very high spatial resolution (VHR) satellite imagery is very promising to precisely monitor the transformation process under the premise of clarifying the canopy reflectance anisotropy of mixed plantations. However, it is almost impossible to understand the canopy reflectance anisotropy of mixed plantations with real satellite data due to the extreme lack of multiangular VHR satellite images. In this study, the effects of the mixture mode on the canopy bidirectional reflectance factor (BRF) were comprehensively analyzed with simulated VHR images. The three ‐ dimensional (3D) Discrete Anisotropic Radiative Transfer model (DART) was used to construct a pure coniferous scene, a pure broadleaved scene, and 27 coniferous–broadleaved mixed plantation scenes containing 3 mixture patterns (i.e., mixed by single trees, mixed by stripes, and mixed by patches) and 9 mixing proportions (i.e., from 10% to 90% with the interval of 10%), and to simulate red (R) and near ‐ infrared (NIR) VHR images for these 3D scenes at both the solar principal plane (SPP) and perpendicular plane (PP) under dif ‐ ferent solar ‐ viewing geometries. Negative correlations were generally found between the canopy BRF and the ratio of conifers in a mixed stand. The anisotropy of conifer dominated plantations is more prominent than broadleaf dominated plantations, especially for the single tree mixture. Alt ‐ hough the level of anisotropy is much lower for PP than SPP, it should not be ignored, especially for the R band. Observations under large viewing zenith angles at PP are more preferred to study the effect of mixing proportions, followed by forward observations at SPP. The R band image has higher potential to distinguish mixture patterns for broadleaf ‐ dominated situations, while the NIR band image has a higher potential for conifer ‐ dominated situations. Furthermore, the canopy BRF generally increases with the solar zenith angle, and one meter can be considered as the optimal spatial resolution for the optical monitoring of the mixture mode. The findings of the current study add some valuable theoretical knowledge for the accurate monitoring of coniferous–broadleaved mixed plantations with VHR imagery.


Introduction
According to the latest national forest resource inventory results of China, the area of China's plantation is 79.54 million hectares, continuously ranking first in the world [1], with subtropical plantations occupying more than half of China's total plantation area [2]. However, the majority of plantations in subtropical China are coniferous monocultures, which was considered as one of the main causes for the low quality of China's plantations with weak ecosystem stability as well as low stock volume [3,4]. Sufficient water and heat resources in subtropical China are particularly favorable natural condi-tions for evergreen broadleaved tree species, which have been harvested and replaced by fast-growing high-yield coniferous monocultures since the 1950s [5,6]. Unfortunately, short rotations (~25 years each) of these coniferous monocultures for the past few decades substantially deteriorated the biodiversity and soil conditions in this area [5]. In order to solve this problem, the Chinese government has strongly campaigned for the forest quality improvement initiative since 2016, in which transforming the even-aged coniferous monocultures into coniferous-broadleaved mixed plantations is one of the preferred silvicultural solutions [2]. The reason is that, compared to monocultures, mixed plantations are considered to have better soil and habitat conditions, higher biodiversity, faster growth rates, more recreational values, lower damages to focal tree species from grazers, pest insects, fungal pathogens, wind and fire, as well as a wider variety of wood and other forest products [7][8][9].
For the sake of quantitatively evaluating the transformation process, it is necessary to observe structural parameters, such as mixing proportions and mixture patterns. Remote sensing overcomes many disadvantages of field investigations and laboratory analysis, and has become an important acquisition technology for forest parameters [10]. Among all kinds of remote sensing data, very high spatial resolution (VHR) satellite images being acquired by many commercial sensors (i.e., IKONOS, QuickBird, Ge-oEye, WorldView, Pléiades, KOMPSAT, SkySat, TripleSat, GaoFen-2, and Deimos) with a spatial resolution higher than 1 m, can be very attractive for accurately monitoring mixed forests, because of their high geometric accuracy and more affordable purchase prices than other VHR images, such as airborne data and UAV (unmanned aerial vehicle) data [11,12]. However, VHR satellites are not usually tasked to perform multiangular observations, which becomes one of the bottleneck constraints for VHR satellites to accurately monitor a forest canopy, because the spectral reflectance of a forest canopy is different under various solar-viewing geometries, which is described as canopy reflectance anisotropy [13,14]. The retrieval of forest parameters using spectral reflectance and vegetation indices should pay attention to canopy reflectance anisotropy. Therefore, understanding the relationships between canopy reflectance anisotropy and target forest parameters is the premise for accurately monitoring forests with VHR optical images. Moreover, the relationships may also be varied with solar angles [15].
Canopy reflectance models provide the possibility to remedy the extreme scarcity of real multiangular VHR satellite imagery. These physical models, which mainly include geometric optical (GO) models, radiative transfer (RT) models, and GO-RT models, were developed extensively over the past few decades. Although these physical models require intensive computations and are limited by the availability of various canopy structural parameters, it is more suitable for the simulation of complex and continuous vegetation canopies [16]. Compared to one-dimensional RT models, three-dimensional (3D) RT models can make greater use of the information that is provided by multi-source remote sensing data through creating and simulating more realistic natural scenes. They are capable of simulating the transfer processes of solar radiation coming from any direction into a natural land surface with complicated 3D structures, as well as simulating remote sensing imagery with any resolution at the top of atmosphere from any direction [17], and therefore are very precise for the estimation of vegetation parameters [18][19][20].
Many 3D RT models have been developed over the past two decades. For instance, RGM (the radiosity-graphics combined model) [21] and RAPID (radiosity applicable to porous individual objects model) [22] are based on radiance algorithms, Raytran [23] and LESS (large-scale remote sensing data and image simulation framework) [24] are based on ray tracing, and DART (Discrete Anisotropic Radiative Transfer model) [25] is based on the discrete coordinate method. DART has been extensively adopted and successfully validated by many scientific researches, either in applications or in modeling [26][27][28][29]. In addition to its accuracy, DART offers the advantage of easily creating any landscape, which is very suitable for the simulation of the directional reflectance of plantations with different mixture modes in this study.
The bidirectional reflectance factor (BRF) of the vegetation canopy in visible and near-infrared (NIR) bands are mainly influenced by three kinds of characteristics, which include the (1) optical properties of the ground and vegetation elements that depend on the canopy's physiological conditions (e.g., leaf water and pigment contents); (2) 3D architecture of individual plant elements, including the 3D distribution and orientation of the leaves; and (3) 2D distribution of the plant elements (e.g., the distance between plants) [29][30][31][32]. For instance, the reflectance of the Norway spruce canopy was simulated, and the biochemical parameters and their relationships with the vegetation indices were studied based on the DART model [31]. The effects of the canopy structures on the simulated top-of-canopy reflectance were evaluated with DART for 3D spruce representations, with different levels of complexity from geometrically simple to highly detailed architectures [32]. Hardly any previous researches have been conducted on the canopy reflectance anisotropic characteristics of mixed plantations. Therefore, in the current study, the relationships between mixture mode and canopy anisotropic reflectance were explored using the DART model, which can simulate tree structures and continuous vegetation components, in order to guide the application of VHR optical imagery over mixed forests.
The objective of this study is to analyze the influence of mixture modes of mixed plantations on canopy anisotropic reflectance. The first step was the construction of 3D forest scenes with various mixing proportions and mixture patterns, followed by the simulation of VHR images with different solar-viewing geometries. Then, the relationships between the three factors (mixing proportions, mixture patterns, and solar-viewing geometries) and canopy BRF of those simulated VHR images were analyzed. Meanwhile, the spatial scale effect analysis was also performed with the intention to obtain an optimal spatial resolution for VHR imagery applications on the mixed plantations.

Study Area
The study area (117°05'~117°40' E, 26°26'~27°04' N) is in the state-owned forest farm of Jiangle county, which is located in the northwest of the Fujian Province. It belongs to the subtropical monsoon climate zone. Based on the monthly temperature and precipitation data from the local meteorological station (26°44′ N, 117°28′ E, and 173.9 m above sea level), we calculated the annual lowest temperature (−2.8~7.5 °C), the annual highest temperature (32.2~38.5 °C), and the annual average precipitation (840~2050 mm) over the period of 1960 to 2012. The main species in this region include Cunninghamia lanceolata, Pinus massoniana, Michelia maccurrei, and Schima superba, with Cunninghamia lanceolata monocultures holding the highest proportions. In the forest quality improvement initiative, Cunninghamia lanceolata monocultures are planned to be transformed into mixed plantations with some native broadleaved species, such as Michelia maccurrei. Our previous research found that many archived even-aged pure coniferous plantations were in fact mixed with some non-coniferous species with different proportions (8.05~16.02%), and these non-coniferous species are the natural regeneration of native species with spatial aggregation patterns over very rugged terrain [33]. Therefore, the accurate monitoring of coniferous-broadleaved mixed plantations using VHR imagery is of practical significance for this region, which further requires the knowledge of reflectance anisotropy of mixed plantations.

Data Acquisitions
The field experiment was carried out in August 2019. At that time, the existing coniferous-broadleaved mixed plantations were conifer-dominated stands with irregular mixture patterns ( Figure 1, and detailed data are shown in Table S1). The plot size is 30 × 30 m, and each plot was divided evenly into 9 small squares to implement accurate tree location measurements and tree structural parameters. The total tree height, height of the first branch, diameter at breast height, tree locations, crown height, and crown radius were measured for each tree. In addition, the species and number of trees were also recorded. The tree height and height of the first branch were measured using TruPluse laser hypsometer. The crown height was determined by subtracting the height of the first branch from the total height, while the crown radius was measured roughly in four directions (north, south, east, and west) according to the crown projections. It should be noted that the main practical use of the field measurements in the current simulation study provides some references for the construction of forest scenes. More specific and extensive applications of these plot data appear in our follow-up researches.

Methods
The overall workflow of the study is illustrated in Figure 2. The main input parameters of DART model were configured according to the empirical values measured from the field investigations and simulated data from the 3-PGmix model (Physiological Principles in Predicting Growth for mixed stands). Pure and mixed plantation scenes were constructed, and then both the red (R) and NIR VHR images of these forest scenes were simulated at a solar principal plane (SPP) and perpendicular plane (PP). SPP has a relative azimuth of 0°/180° to the sun. PP is perpendicular to SPP and the relative azimuth is therefore 90°/270° to the sun. Almost no visible shaded elements can be observed from the backward scattering direction at SPP when the sun is behind the observer, which creates a typical phenomenon called "hotspot". On the contrary, viewing from the forward direction at the SPP that is opposite to the sun presents the largest amount of shadows, which is the so called "darkspot". Observations at SPP are helpful to assess the canopy structure of forests, since the differences in the forward and backward directions are dependent on the horizontal and vertical structure of the observed objects [34]. The anisotropic effect of the canopy BRF is generally smaller at PP, since neither the brightest "hotspot" nor the darkest "darkspot" exists at this plane, which, however, has its own advantages for VHR imagery applications on the forest canopy.

The DART Model
DART is one of the most accurate and comprehensive 3D RT models developed in CESBIO (Center for the Study of the Biosphere from Space, UPS-CNRS-CNES-IRD, France) since 1992. Its standard mode, called DART-FT, can simulate the BRF of any 3D landscape (for instance, forest stands, crops, and urban) with any spatial and spectral resolution under any solar-viewing geometry. To date, DART has been successfully employed in various scientific applications, and it has been successfully validated against field measurements and through Radiation Model Intercomparison (RAMI) [35][36][37].
A DART scene can range from simple structures, such as turbid layers, to highly detailed 3D objects. The scene can contain triangles and/or turbid media. A turbid medium that represents the leaves inside a crown is defined by its leaf angle distribution (LAD), leaf surface density (m 2 /m 3 ), and individual leaf optical properties (reflectance and transmittance spectra). Individual trees can be described simply as elementary geometrical shapes (conic, spherical, or ellipsoidal) filled with a turbid medium containing a mix of foliage and woody elements, combined with triangles to represent the trunk. Thus, the parameters defining a single tree include the structural parameters (i.e., size and shape of the crown, and distribution of the leaves within the crown) and optical parameters (e.g., optical properties of leaves and branches). Generally, a DART scene is defined by a voxel grid with a predefined size (e.g., 0.1 m in this study).
The new Monte Carlo mode, called DART-Lux, was developed in 2018. Its accuracy for simulating images, fluorescence, and LiDAR signals is similar to DART-FT, but with an efficiency that is a hundred times greater, since it uses a "Ray-Carlo" method (a combination of ray-tracing and Monte Carlo methods) in a multi-threaded algorithm to reduce the computation time. It adapts a Bidirectional Path Tracing (BPT) algorithm from LuxCoreRender, that solves the light transport equation with Monte Carlo integration techniques. The BPT algorithm constructs paths that start from the camera at one end, and from a light source at the other end, and connects them [26]. The version 5.7.11 of DART-Lux was adopted by our study. Further details about the standard mode and the Monte Carlo mode of the DART model can be found in [26,[38][39][40]].

The 3-PGmix Model
3-PG (Physiological Principles in Predicting Growth) model is a simple process-based model to predict the growth and development of even-aged stands [41]. The core of this model is to dynamically simulate the gradual decline in solar radiation, carbon quantity balance, and water balance through a series of equations, and then to simulate the dynamic change of the forest division system and the external environment [42]. The 3-PGmix was developed from 3-PG to simulate the growth of mixed species forests with a new light sub-model, different versions of diameter distributions, and water balance calculations. It requires only readily available site and climatic data as the inputs, and predicts the time-course of stand development on a monthly basis in a form familiar to forest managers [43,44].
The 3-PGmix was used to simulate some forest growth parameters of a coniferous-broadleaved mixed plantation for the 3D forest scenes configuration of DART model (Table 1; Figure 3). Climatic data in our study area, including the monthly maximum and minimum temperatures, and the maximum and minimum precipitation levels, were obtained from China Meteorological Data Network (http://data.cma.cn/ (accessed on 15 October 2021). Other parameters were mainly determined through either the default values and statistical fittings, or according to a very representative reference [43][44][45][46]. The simulation of this study set 2000 as the initial year and 2100 as the end year. The time scale was set as 1 year to correspond with the temporal sampling of the climatic data. Since the maturation of Cunninghamia lanceolata is around 20 years [47], the tree structure parameters at the 20th year were retrieved from the results of the 3-PGmix as inputs for the DART model.

Parameterization for Simulation
Spectral heterogeneity also has important effects on the application of VHR imagery [48]. Typical spectral reflectance of different components, such as ground, leaf, and trunk for coniferous and broadleaved species in the spectral library of DART, were chosen for the VHR imagery simulations (detailed spectral reflectance of the different components are shown in Figure S1). In fact, 0.66 μm and 0.83 μm were set as the central wavelengths of the R and NIR bands for the simulated VHR images, respectively, to be consistent with the band settings of the most common VHR satellites, such as QuickBird and Worldview-2/3. The spectral reflectance data of the these two wavelength for different components are given in Table S2.
In our study, the DART model was used to construct mixed plantation scenes with three different mixture patterns (mixed by single trees, mixed by stripes, and mixed by patches), which are typical for mixed plantation management. Each mixture pattern was set to have 9 different mixing proportions, from 10% to 90%, with the interval of 10%. Among these mixing proportions, the percentage represents the ratio of conifers to the total number of trees in the scene. For example, 10% means that 10% of the trees in the scene are conifers and the other 90% are broadleaves. Therefore, 27 mixed plantations scenes, as well as a pure coniferous scene and a pure broadleaved scene were involved in the simulation. Figure 4 illustrates several typical spatial configurations of plantation scenes for the study. In order to smoothly present the variation of canopy BRF over the 2 observation planes, the viewing zenith angle (VZA) was set from −85° to 85° with the interval of 5° for the current study. Three different sets of solar zenith angle (SZA) and solar azimuth angle (SAA) were involved in the simulation, which were calculated for three local time (i.e., 9:30 a.m., 12:30 p.m., and 15:30 p.m. on 2nd July 2012). Therefore, a total of 6090 pairs of R and NIR VHR images were simulated to study the mixture mode effects on the reflectance anisotropic characteristics of the mixed plantations. Table 1 shows the specific settings of the input parameters for scene constructions and image simulations with DART, in which the tree parameters for broadleaf trees and conifer trees were generated from 3-PGmix (Figure 3).

Relationship Analysis
The relationship between canopy BRF and the mixture mode was performed from three perspectives ,which include (1) the effect of mixing proportions on canopy BRFs, (2) the effect of mixture patterns on canopy BRFs, and (3) the effects of solar-viewing geometries on canopy BRFs.
In order to focus on the effect of mixing proportions, one SZA (i.e., 24.4°) and one mixture pattern (i.e., mixed by single trees) were fixed for the analysis. In order to directly compare the level of anisotropy, we calculated the Normalized Difference between Extreme Values (NDEV) for each BRF curve as: which originated from Leblanc et al. [49] and Chen et al. [50]. A higher NDEV corresponds to higher levels of anisotropy. With NDEV, it is easy to identify the relative uncertainties caused by neglecting reflectance anisotropy in different situations.
When it comes to the effect of the mixture patterns, one SZA (24.4°) and one VZA (25°) were fixed for the analysis. In fact, this is the hotspot area with fully illuminated canopy surfaces. The coefficient of variation (CV) of the three mixture patterns was adopted to compare the variation of canopy BRFs among mixture patterns with different mixing proportions. The CV is the ratio of the standard deviation to the mean value, which can be used to reflect the relative variation of distribution independent of the unit of measurement. In addition, a SZA (24.4°) and a mixing proportion (50%) were fixed to analyze the variation of reflectance anisotropic characteristics with the mixture patterns, which can not only show the effect of the mixture patterns on the anisotropic characteristics, but can also indicate that certain viewing angles are more effective than others to differentiate the various mixture patterns with the same mixing proportions.
The effect of solar-viewing geometries on canopy BRFs was analyzed under three sets of solar angles at SPP. Under each solar-viewing geometry, the average NDEV, mean values, and the standard deviation of canopy BRFs of the scenes under 29 scenarios (27 mixed ones and 2 pure ones) were calculated and analyzed.

Spatial Scale Analysis
A higher spatial resolution of remote sensing imagery does not always guarantee better accuracy in either qualitative identification or quantitative inversion [51]. An optimal spatial resolution should therefore be determined before practical applications. In other words, different applications should have different selection criteria for the optimal spatial resolution, which not only meet the needs of certain research purposes, but also have the least data redundancy.
Optimal scale research methods, such as the average local variance (ALV), variogram, and least spectral difference, were widely used for the spatial scale analysis [52][53][54][55]. ALV is a texture analysis method, which is suitable for analyzing small scale images [54]. Therefore, it was preferred by the current study. First, a moving window was set on each image, the variance of all pixel values in the moving window was calculated through the whole image. Second, all the variance values were averaged to obtain the ALV of the whole image. Finally, the local variance map was drawn with the spatial resolution as the x-coordinate and the average local variance as the y-coordinate. Additionally, the optimal resolution was determined when the parabola reached the highest point for the ALV curves [55]. Moreover, the CV was used as an auxiliary means for the spatial scale analysis. Theoretically, the CV decreases with the decrease in the spatial resolution, until it reaches the threshold of optimal spatial resolution. In this study, the CV was used to depict the level of heterogeneity among the pixels.
The range of spatial resolutions for the spatial scale analysis was set from 0.1 m to 10 m (i.e., 0.1 m, 0.3 m, 0.5 m, 1 m, 5 m, 10 m), in order to cover the airborne remote sensing, VHR commercial satellites, as well as some high spatial resolution multi-spectral optical satellites, such as SPOT series and ALOS-AVNIR2. It should be noted that the range of spatial resolutions was limited by the predefined plot size, which was reset to 120 × 120 m for the spatial scale analysis in the study. Scene parameters for the 5 extra simulation experiments, except 0.1 m, were generally the same as in Table 1, among which SZA, the mixing proportion, and mixture pattern were set as 24.4°, 50%, and single trees, respectively.

Simulations of Sample Plots
The canopy BRFs of the six sample plots have different levels of anisotropic characteristics ( Figure 5). For instance, the R band BRFs of plot 2 are much higher than that of the other plots (Figure 5a,b), probably because of the higher exposure of soil background as well as the lower probability of mutual shadowing in plot 2 (Figure 1). At SPP, greater BRF differences among the six plots are around the hotspot area, rather than other forward and backward viewing angles for the R band (Figures 5a and 6a), while the reverse can be seen for the NIR band (Figures 5c and 6b). This suggests that the hotspot area is relatively more appropriate to monitor the variation caused by mixture modes than the other viewing directions at SPP for the R band rather than the NIR band. Bell-shaped patterns for the the R band ( Figure 5b) and bowl-shaped patterns for the NIR band (Figure 5d) are shown at PP. At PP, greater BRF differences among the six plots are around the nadir observation directions for the R band (Figures 5b and 6a), and are shown at larger VZAs for the NIR band (Figures 5d and 6b).
Generally, canopy BRF variations among the six plots are higher at SPP than PP, especially around the hotspot area for the R band (Figure 6a), while higher at PP than SPP for the NIR band, especially in the backward-facing observation directions ( Figure  6b). The similar reverse trend of canopy BRF variations is shown between the R band and the NIR band at PP, in which the greatest BRF differences for R band and the lowest BRF differences for the NIR band are all around the nadir observations ( Figure 6).
However, it is hard to ascribe all the BRF variations among the six plots to the effect of mixture mode on the canopy reflectance anisotropic characteristics, because of the irregular spatial distribution and complicated tree structures in the field plots. Specific knowledge about the theoretical effect of mixture mode on the canopy anisotropic reflectance of mixed plantations should be derived from hypothetical simulations.

The Effects of Mixing Proportions on Canopy BRFs
Generally, the BRF shapes under different mixing proportions are similar for the same observation planes, and canopy BRFs of conifer-dominated plantations are generally lower than broadleaf-dominated plantations. Obvious hotspots and darkspots can be seen at SPP (Figure 7a,c), and symmetrical patterns around the nadir are shown at PP, with the combination of general "bowl-shaped" patterns and local "bell-shaped" patterns around the nadir views (Figure 7b,d). Relatively, the "bell-shaped" patterns around the nadir views are more obvious for conifer-dominated plantations, and are also more obvious in the R band than in the NIR band (Figure 7b,d).
The maximum canopy BRF of the R band is above 0.03 at SPP, which is around the hotspot for the pure coniferous plantation, while the minimum is as low as 0.01 around the darkspot for the mixed plantation with 80% of conifers (Figure 7 a). It can be seen that the canopy BRF difference between the hotspot and darkspot is larger for conifer-dominated plantations than broadleaf-dominated plantations (Figure 7e,f). In another word, the anisotropy of conifer-dominated plantations is more prominent than broadleaf-dominated plantations. NDEV being 0.5 means that the maximum value of the BRF curve is three times higher than the minimum; being 0.3 means twice as high. The highest NDEV (~0.52) appears for mixed plantations with about 80% of conifers at SPP. All the NDEVs of the R band under all the mixture modes are higher than 0.32 at SPP. Except the NDEV of the R band at SPP, all the other three groups of NDEVs are lower than 0.3. In fact, for SPP, the maximum BRF is at a hotspot and the minimum is at a darkspot. The level of anisotropy is of course much higher for SPP than PP, and is much higher for the R band than the NIR band. At SPP, the level of anisotropy first increases and then decreases with the increase in the broadleaf proportions (Figure 7e,f). The CV was used to describe the difference of canopy BRFs among different mixing proportions, and to compare the intra differences between the two bands. A larger CV represents a greater difference among the different mixing proportions under the same set of solar-viewing geometries, which also means relative greater effects of mixing proportion on the canopy BRF. Large viewing angles are obviously better than small ones in highlighting the differences among different mixing proportions at PP for both bands. The CVs of the R band were generally higher than that of the NIR band, which means R band imagery has a relatively higher potential to monitor the variation of mixed plantations than the NIR band (Figure 8).

The Effects of Mixture Patterns on Canopy BRFs
Plantations with different mixture patterns have different spatial heterogeneity, which leads to different levels of canopy BRF anisotropy. At SPP, the level of anisotropy is the highest for plantations that contained mixed single trees, followed by plantations that contained a mixture of patches and stripes for the conifer-dominated plantations (Figure 7e,f). For broadleaf-dominated situations, plantations mixed by patches seem to have the highest level of anisotropy. The differences in the levels of anisotropy among the three mixture patterns are obviously smaller for broadleaf-dominated plantations than conifer-dominated situations.
Canopy BRF differences among the three typical mixture patterns for different mixing proportions were further illustrated in Figure 9. The general trend of canopy BRF variations with mixing proportions is similar for the three mixture patterns. In the R band, all the canopy BRFs of the three mixture patterns slightly decrease first and then increase sharply with mixing proportions, with the lowest at around 50% for both patterns that were mixed by single trees and stripes, and at around 70% for patches ( Figure  9a). In the NIR band, the canopy BRFs start to increase rapidly almost at the beginning of the increase in broadleaf proportions, and then tend to stabilize or even slightly decline after 50% (Figure 9b). The CV of canopy BRFs among these three mixture patterns at different mixing proportions and different VZAs were used to more straightforwardly compare their differences ( Figure 10). Overall, the CV gradually rises with the increase in broadleaf proportions for the R band, with the lowest at 20% of broadleaf and the highest at 20% of conifers (Figure 10a). It means that the effect of mixture patterns on canopy BRFs in the R band is relatively higher for broadleaf-dominated plantations and relatively lower for conifer-dominated plantations. On the contrary, the effect of the mixture patterns on canopy BRFs is relatively higher for conifer-dominated plantations than broadleaf-dominated plantations in the NIR band (Figure 10a). The results suggest that it should be relatively easier for the R band VHR images to distinguish the mixture patterns for broadleaf-dominated situations and relatively easier for the NIR band VHR images to distinguish the mixture patterns for conifer-dominated situations. The differences among the three mixture patterns appear to be much smaller when it comes to VZAs (Figure 10b). The differences of BRFs among the three patterns are the most obvious around the hotspot areas in the R band, while more profound differences can be observed with forward viewing directions than backward in the NIR band (Figure 10b).

The Effects of Solar-Viewing Geometries on Canopy BRFs
The mean values and the standard deviations of canopy BRFs at SPP for all 29 simulated forest scenes (27 mixed stands that contain 3 mixture patterns and 9 mixing proportions, and 2 pure stands) were used to analyze the effect of solar-viewing geometries on canopy BRFs (Figure 11). The canopy BRFs of both the R and NIR bands increase substantially with SZA in oblique observations (i.e., large VZA) for both backward and forward directions, while they decrease with SZA for small VZAs (−35° < VZA < 35° for the current study), although the differences are not as obvious as those in the large angle oblique observations (Figure 11a,b). The standard deviations are larger for large angle oblique observations, which indicate that canopy BRF differences among these 29 mixture modes can be observed more clearly under larger VZAs at SPP (Figure 11c,d). The range of canopy BRF at SPP is wider for larger SZA for both bands (Figure 11e,f), which means canopy BRFs under larger SZAs contain more information of mixture mode. The maximum BRFs clearly increase with SZA for both the R and NIR bands, but the minimum BRFs for different SZAs are almost the same. The mean NDEV of the 29 scenes is 0.542 for the R band when SZA is 54.3°, which means that the maximum R band BRFs are more than 3 times higher than the minimum BRFs. The maximum BRFs are about twice as high as the minimum BRFs when SZA is 22.4° for the R band, and when SZA is 54.3° for the NIR band.

Optimal Spatial Resolution
Two spatial scale analysis methods were implemented with the intention to retrieve an optimal spatial resolution for optical images to observe the mixture modes of mixed plantations. Resolutions from 0.1 m to 10 m were analyzed. The ALV increased steadily with the decrease in spatial resolution, and it remained unchanged when the spatial resolution decreased to 1 m. The CV sharply decreased with the decrease in the spatial resolution, with an obvious inflection point around the spatial resolution of 1 m. Both measures indicate that the optimal spatial resolution for the observation of mixed plantations is around 1 m (Figure 12).

Influence of Mixture Modes on Simulated Canopy BRFs
In this study, the mixture mode was studied from the perspective of mixing proportions and mixture patterns. The influence of mixing proportions on canopy BRF is more obvious than the influence of mixture patterns on canopy BRF. In other words, it should be easier for VHR satellites to monitor the mixing proportions under the same mixture pattern than to distinguish the different mixture patterns under the same mixing proportion. Generally, a positive correlation was found between canopy BRFs and broadleaf proportions, which is consistent with existing researches, such as [32] and [56]. For instance, Kuusk et al. [56] found that the reflectance of the NIR band increases with the proportion of broadleaved species for the dense coniferous-broadleaved mixed plantations. In fact, the relationship between canopy BRF and broadleaf proportion for each band depends on the spectral reflectance of the different components (i.e., ground, tree trunks, and leaves) and the proportion of each component in a single scene. The reason for the general positive correlation between the canopy BRF and broadleaf proportions is due in part to the higher reflectance of broadleaf canopies (Table S2).
The CVs of the R band BRFs among the different mixing proportions and among different mixture patterns were generally higher than that of the NIR band. It means that the R band imagery has a relatively higher potential to monitor the variation of mixed plantations. According to Ren [57], more differences were indeed found among the canopy spectra of coniferous, broadleaved, and mixed forests in the R band, while more spectral overlaps in the NIR band were found between the mixed and coniferous forests. The distinctive contribution of our study is to provide more detailed instructions for monitoring mixed plantations, such as selecting more preferred solar-viewing geometries.
To the best of our knowledge, the effects of the mixture mode on the level of reflectance anisotropy of mixed plantations have not been previously studied. Theoretically, the level of anisotropy is of course much higher for SPP than PP, and is much higher for the R band than the NIR band. In our study, in addition to confirming and quantifying these theories, we further obtained the explicit knowledge about the relationship between mixing proportion and the level of anisotropy under different mixture patterns. The knowledge should be very instructive for the accurate monitoring of mixed plantations with VHR imagery.

Effects of Different Solar-Viewing Geometries on Canopy BRFs
Understanding the effect of solar-viewing geometries on canopy BRFs is very important for the application of remote sensing images. Canopy BRFs observed at hotspots can be more than three times higher than those observed at darkspots in our study, which indicates the high level of canopy reflectance anisotropic characteristics. The level of anisotropy is different for various mixture modes and different SZAs (Figures 7e,f and 11e,f), which could provide valuable guides for the serious treatment of VHR imagery over the application of mixed plantation monitoring. Observations under some solar-viewing geometries were found to contain relatively more information about mixture mode than other solar-viewing geometries (Figures 8 and 10b), such as the forward observations at SPP, oblique observations with very large VZA, and larger SZAs. These are reasonable because the mixture mode can be seen as one of the measures of forest structures. With more shadows being cast under larger SZAs, more mutual shadows being observed in the forward directions at SPP, and a longer path length with larger VZAs, the range of canopy BRFs is significantly widened for the whole forest scene, which can highlight the differences among the different mixed plantations.

Optimal Spatial Resolution
In order to better distinguish the mixed plantations with different mixture modes, remote sensing images with enough detailed information should be preferred. The spatial resolution analysis in this study demonstrated that 1 m can be considered as the optimal spatial resolution for the observation of mixed plantations. In fact, the range of the spatial resolutions in our study was only set from 0.1 m to 10 m with the consideration of intensive computation, therefore, 1 m might be the local optimal rather than the global optimal spatial resolution. Moreover, it should also relate to the tree size parameters that are listed in Table 1. More spatial scale analysis experiments should be performed over a wider range of spatial resolutions and different settings of tree parameters with different dimensions of 3D scenes.

Suggestions for Future Research
Although some very instructive knowledge about the relationship between the mixture mode and canopy reflectance anisotropy has been summarized in the current study, the many limitations should be mentioned. For instance, only two sets of tree parameters were used to build those 3D mixed plantation scenes with three hypothetical mixture patterns. Only one possibility of the spatial distribution of the two species was chosen for each pair of mixing proportion and mixture pattern. Cone and ellipsoid were chosen, respectively, to represent the conifer and broadleaf tree crowns in a set of triangles with different optical properties.
Therefore, in the future, the mixture mode effects on the canopy BRFs could be studied using more realistic scenarios over larger spatial scales. For example, more than two species can be involved in the simulation. Even though only two species were considered, trees can still have different settings of structural parameters to approximate uneven-aged mixed plantations. More possibilities of the spatial arrangement of the tree locations for each mixture mode should be managed to enhance the credibility of the conclusion of the current study. The 3D objects of each species should be created or imported into DART to build more vivid forest scenes. Moreover, detailed field measurements of the mixed plantations, precise tree structural parameters, and the accurate spectra of each component (e.g., needles, leaves, trunks, and soil background), would further improve the simulation of the DART model for mixed plantations.

Conclusions
The precise monitoring mixed plantations with VHR imagery requires an understanding of the effects of mixture mode on canopy anisotropic reflectance. We analyzed the relationship between mixture modes and canopy anisotropic reflectance using the DART model. The main results indicate that (1) the mean canopy BRFs of the mixed plantations generally increase with broadleaved species proportions; (2) the anisotropy of conifer-dominated plantations is more prominent than broadleaf-dominated plantations, especially for the mixture pattern of single trees; (3) solar-viewing geometries have obvious influences on canopy BRFs, the canopy BRF generally increases with SZA, observations under large VZAs at PP are more preferable to study the effect of mixing proportions, followed by forward observations at SPP; (4) the R band image has a greater potential to distinguish the mixture patterns for broadleaf-dominated situations, while the NIR band image has a higher potential for conifer-dominated situations; and (5) 1 m can be considered as the optimal spatial resolution for VHR images to be applied to mixed plantations. Hopefully, this research can be treated as the reference for the selection and accurate application of VHR images in monitoring coniferous-broadleaved mixed plantations.

Supplementary
Materials: The following are available online at www.mdpi.com/article/10.3390/f13020235/s1, Figure S1: Spectral reflectance of each component of coniferous-broadleaved mixed plantations, Table S1: Six measured plots data, Table S2: Spectral reflectance data of the selected two wavelengths for different components. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data will be made available on Figshare upon paper acceptance.

Acknowledgments:
We gratefully appreciate Huaguo Huang, Jianbo Qi and Linyuan Li from the Beijing Forestry University for their valuable suggestions about this work. We also thank Xiaolong Yu from the Jiangle State-Owned Forest Farm for sharing valuable information about mixed plantations, as well as his assistance with field data collections in the study area. We would also like to thank Yingjie Wang from CESBIO for the technical assistance about the DART model. The most special gratitude goes to the academic editor and two anonymous reviewers for their constructive comments to give this work an overturning chance to be much more instructive for real applications.

Conflicts of Interest:
The authors declare no conflicts of interest.