Forest canopy scattering properties with signal of opportunity reflectometry: theoretical simulations

In recent years, signal of opportunity reflectometry (SoOp-R) has become a promising remote sensing technique. This emerging technique employs the reflected signals from existing Global Navigation Satellite System (GNSS) or communication satellites to estimate geophysical parameters for Earth observation, such as wind speed, altimetry, significant wave height, soil moisture, etc. While its application for forest canopy monitoring is still in the initial stage, there are still many unknown relations between vegetation parameters and actual observations, and a proper theoretical basis needs to be established for simulation and analysis of the different observation geometries. In this paper, we develop a bistatic scattering model with various polarizations at different frequency bands. Our improved model is based on the first-order radiative transfer equation, and is developed based on the wave synthesis technique, after which it can be used for circular polarization signals in bistatic radar systems, i.e. the typical configuration of SoOp-R. We analyze the simulations of the P (0.25–0.5 GHz), L (0.5–1.5 GHz), C (4–8 GHz), and X (8–12 GHz) bands at the backscattering, specular cone, bistatic scattering, and perpendicular planes. The contributions of the different components to the total scattering are also analyzed. The results show that the coherent scattering at the specular cone is larger than the non-coherent scattering, while trunk-dominated forest canopy has strong scattering at the aforementioned different directions. Variations of canopy parameters such as trunk and branch diameters, tree density, and vegetation water content are also simulated at the specular cone plane, showing strong dependence on the final bistatic scattering observation. The simulation results show that the SoOp-R technique has a great potential for monitoring of canopy parameters.


Introduction
Signal of opportunity reflectometry (SoOp-R) is an innovative Earth observation technology that employs the existing satellite signals for terrestrial remote sensing. A clear example of this technology, recently under development, is the emerging Global Navigation Satellite System-reflectometry (GNSS-R) (Zavorotny and Gleason 2014;Cardellach et al. 2016;Edokossi and Calabia 2020). The basis of this technique aims to employ the GNSS signals reflected by the Earth's surface and received by a dedicated antenna. This new technology offers numerous advantages for surface geophysical parameters detection, including low cost, low power consumption, wide-coverage, and a high spatial and temporal resolution. Compared to the traditional monostatic microwave remote sensing with polar-orbiting satellites, the GNSS-R technology can provide continuous observations over the Earth's surface, and surely will become a powerful complement to the traditional microwave techniques.
GNSS-R has been used in ocean remote sensing due to the relative uniformity of the ocean medium, which allows neglecting the polarization characteristics Mayers and Ruf 2019). However, the changing geometry characteristics of GNSS-R observation over the complex land surface led to the initial study of land surface parameters, such as soil moisture and vegetation remote sensing (Li et al. 2017;Calabia et al. 2020;Jia and Savi 2016).
At present, the coherent signal of the direct and reflected signals observed by ground-based GNSS receivers is used to extract the signal-to-noise (SNR) ratio, so that the phase and effective reflectometer height can be used to determine soil moisture, snow thickness, and vegetation water content (Chew and Small 2013;Rodriguez-Alvarez and Camps 2010;Jin et al. 2016). In recent years, the GNSS-R receiver has been deployed onboard satellites, such as the Cyclone Global Navigation Satellite System (CYGNSS) (Larson 2016) satellites launched by the National Aeronautics and Space Administration (NASA) in December 2016, which have provided the most valuable dataset for GNSS-R research and applications. Although the initial objective of these satellite missions was to detect the wind speed for tropical cyclones, recent studies have exhibited the capabilities for sensing land surface attributes, including flood inundation (Ruf et al. 2013), soil moisture (Chew et al. 2017), and wetland extent (Yan and Huang 2020).
SoOp-R can employ numerous sources of radio signals, including the P (0.25-0.5 GHz), L (0.5-1.5 GHz), C (4-8 GHz), or X (8-12 GHz) bands. For instance, the P-band was employed to study land surface soil moisture and snow characteristics (Morris et al. 2019;Garrison 2019;Yueh et al. 2018). Currently, there are relatively few studies on vegetation monitoring using microwave SoOp remote sensing, and the numerous SoOp from communication satellites provide an unprecedented opportunity for microwave remote sensing reflectometry (Shah et al. 2019;Kurum et al. 2019). Previous studies on vegetation monitoring through SoOp-R were focused on correlation coefficients between the reflected signals and vegetation parameters. Recently, some researchers have employed TDS-1 and CYGNSS data to study the potential of GNSS-R for evaluating forest biomass Carrenoluengo et al. 2020). However, most of the existing works only focus on experimental analyses (such as ground-based and space-borne data), and less attention has been paid to the scattering mechanisms (Shah et al. 2019;Kurum et al. 2019;Eroglu et al. 2019;Carrenoluengo et al. 2020;Santi et al. 2020;Ferrazzoli et al. 2011). The experimental research without a physical basis is difficult to promote, and it greatly hinders the applications of SoOp-R over land surfaces. Besides, scattering characteristics have not been deeply explored yet, and the analysis of different observation geometries can provide new insights for future research. Therefore, it is urgent to reveal the physical laws that govern SoOp-R for land surfaces and clarify the mechanisms at different observation geometries.
This article focuses on vegetation monitoring through SoOp-R for modeling the scattering characteristics of forest canopies using the Michigan Microwave Canopy Scattering (MIMICS) model (Ulaby et al. 1988). The results provide a theoretical basis for the development of experimental inversion algorithms to estimate land surface parameters and design new sensors. This will contribute to experimental sensor design, model simulation, and data acquisition, interpretation, analysis, and model assimilation. This paper is organized as follows: Sect. "Theory and methods" introduces the theoretical basis of SoOp-R for forest canopy studies. In Sect. "Simulation results and analysis", the simulation and analyses of the different scattering characteristics are presented, and Sect. "Conclusions" summarizes the results, conclusions, and future research.

Bistatic scattering geometry
The scattering geometry of the bistatic radar system is shown in Fig. 1, where θ and ϕ are the zenith and azimuth angles of the incoming signal, and the subscripts i and s represent the incident and scattering components, respectively. Note that the SoOP-R forms a bistatic radar system, and the BackScatter Alignment (BSA) convention is the standard for radar polarimetry. In order to obtain the different polarization combinations using the wave synthesis technique, we transform the polarization coordinate systems accordingly.

MIMICS model
The MIMICS model (Ulaby et al. 1988) is a very popular model used for backscattering systems, usually for monostatic radars. Unfortunately, this model cannot be directly used for SoOP-R, since the transmitters and receivers in SoOP-R follow a typical bistatic radar system (Fig. 1). Therefore, here we modify the backscattering model so that it can be used for the bistatic scattering systems. The method is based on adding the scattering geometry in the phase and extinction matrices implemented in the MIM-ICS model. The bistatic radar scattering model for the forest canopy (Bi-MIMICS model) uses an iterative algorithm to solve the radiation transfer equation (Ferrazzoli et al. 2011;Ulaby et al. 1988). The following equations are available in the original MIMICS handbook (Ulaby et al. 1988). However, the model is only provided in the backscattering mode. Here, we modify the model and obtain the bistatic radar form (Ferrazzoli et al. 2011;Ulaby et al. 1988), which is the typical form for SoOP-R. The development of the model is shown from Eq. 6 to Eq. 20. Note also that we have included both scattering angles zenith and azimuth.
The model simplifies the forest stands into 3 layers: the crown layer, the trunk layer, and the ground layer. As shown in Fig. 2, the crown layer is modeled in terms of the distribution of dielectric cylinders and disks, while  Scattering mechanisms in the first-order MIMICS model, including the GCG, CG, DC, GC, GT, DG, and TG terms. The SG term is not shown. The crown layer depth is Z 1 = d and the trunk layer depth Z 2 = Ht the trunk layer is treated as cylinders of uniform diameter. To simplify the calculation, we assume that the incident azimuth is ϕ ii = 0°. Then the angular relationship of backscattering is θ s = θ i , ϕ s = 180°; the angular relationship of mirroring scattering is θ s = θ i , ϕ s = 0°; and the angular relationship of forward scattering is θ s = 180°θ i , ϕ s = 0°. After the incident energy is scattered by the particles, the intensity of the scattered energy I s = (θ s , ϕ s ) and the intensity of the incident energy I i = (θ i , ϕ i ) can be related through the modified Mueller matrix L m : In this equation, (θ k , ϕ k ) is the particle orientation, r is the distance between the incident energy and the particle, and the modified Mueller matrix L m is defined by the electric field scattering matrix S as: In this matrix, the subscripts v and h indicate the vertical and horizontal polarizations. The superscript * is thetions. The superscript * is the conjunction, and ℜ and ℑ are the real and imaginary parts of the complex value. Finally, η is the intrinsic impedance. The first-order bistatic transformation matrix that connects the incident intensity and the scattered intensity is as follows: where the transformation matrix T is represented by the phase matrix and the extinction matrix, µ 0 and ϕ 0 are related to the incident angles. Both are calculated by the average modified Mueller matrix. The phase matrix is as follows: In this equation, the size, orientation, and distribution function of the scatterer are s k , (θ k , ϕ k ) , and f (s k ; θ k , ϕ k ) , respectively. N k is the density of the scatterer. The extinction matrix can be expressed as: (2) In this equation, is the average scattering amplitude coefficient, K indicates the type of the scatterer, pq is the polarization, and k 0 is the free space wavenumber.

Scattering mechanisms
There are 8 scattering mechanisms in the MIMICS model (Ferrazzoli et al. 2011;Ulaby et al. 1988). These include the direct specular scattering term (SG); the random rough surface term (DG); the direct crown bistatic scattering term (DC); the ground reflection-crown scattering-ground reflection term (GCG); the crown scattering-ground reflection term (CG); the ground reflection-crown scattering term (GC); the ground reflection-trunk scattering term (GT); and the trunk scattering-ground reflection term (TG). The graphic description of these terms is shown in (Fig. 2; Ferrazzoli et al. 2011;Ulaby et al. 1988), and the combined term is as follows: In the SG term, the direct signals are propagated through the crown and trunk layers and then reflected by the ground layer in the specular direction. Then, the energy is propagated upwards through the trunk and the crown layers. The whole path can be formulated as follows: In the CGC term, the incident energy is firstly propagated through the crown and trunk layers, then scattered by the ground layer, and then through the crown layer, where the volumetric scattering occurs. Part of the energy is scattered by the ground layer and then the upward signal is propagated through the trunk and crown layers. This process can be formulated with the following equation: In the CG term, the incident energy is first propagated by the crown and trunk layers, scattered by the ground, and then the signals are scattered by the crown layer. It can be formulated as follows: In the GC term, the energy is first scattered in the crown layer and then propagated by the trunk layer. Then the energy is scattered by the ground layer and then propagated by the trunk and crown layers. This process can be expressed as follows: In the DC term, the energy does not propagate in the boundary layer, but it is only scattered by the crown layer: In the TG term, the energy first propagates through the crown and trunk layers and it is scattered at the ground. Then a bistatic scattering occurs at the trunk layer, finally, the energy is propagated upwards through the trunk and crown layers: The GT term refers to the same path as the TG term, but in an inverse direction: Finally, the RG term represents the energy first propagated through the crown and trunk layers, at the ground layer bistatic scattering occurs, and the energy finally is propagated upwards through the trunk and the crown layers: For the above terms, k is the extinction matrix, ± indicate upward/downward directions, c and t indicate the crown and trunk layers, respectively, R is the reflectivity matrix of the specular surface, G is the rough surface scattering matrix, and A represents the scattering that occurs in the crown and trunk layers, which are calculated by the phase and the extinction matrices. For more details of the scattering models, please see the corresponding reference (Ulaby et al. 1988).

Wave synthesis technique
In microwave SoOp-R, the reflected signal on the ground surface is relatively weaker than the direct signal, and the receivers need specific features. These include accounting for the different polarizations, so that the reception of the reflected signal is strong enough for measurements. Unfortunately, the existing microwave scattering models are only developed for linear polarization, because these are usually used for the traditional radiometry and monostatic radars. However, the new emerging SoOP-R remote sensing technique uses the circular polarization, since it needs to overcome the ionospheric effects. Therefore, by employing the wave synthesis technique, here we present the required modifications to the existing scattering models, so that we obtain a model capable of estimating the bistatic scattering at various polarization combinations. The scattering characteristics for different polarizations are obtained according to the method of wave synthesis (Liang and Pierce 2005) as follows: In this equation, Mm is the modified Mueller matrix, and Ym is the modified stokes vector, where the upper subscript r and t indicate the polarization of transmitted and received signals. These depend on the variables of orientation angle ψ and ellipticity angle χ, respectively. After changing the orientation angles and ellipticity angle, we can get the bistatic scattering properties at various polarizations. (14)

Model validation and future prospects
The MIMICS model in its backscattering form is available to the scientific community and it has been validated with in situ measurements (Ulaby et al. 1988 (Ferrazzoli et al. 2011;Ulaby et al. 1988). In SoOP-R the transmitted signals are in right hand circular polarization to overcome the ionospheric effects. Since the existing forms of the scattering models are in linear polarization (Ferrazzoli et al. 2011;Ulaby et al. 1988), here we employ the wave synthesis technique to obtain the various polarization combinations. We modify the Stokes vectors to the linear polarization form and find out that the results are similar to that of the linear form of the model. Future works are addressed to validate the model with in situ measurements.

Simulation results and analysis
In this section, we simulate and analyze the scattering characteristics of various polarizations under different observation geometries, as well as the response of the different canopy parameters to the bistatic scattering characteristics. It is worth mentioning that once we have simulated some properties of the vegetation in the paper (Wu et al. (Under review)), where give the bistatic scattering of forest canopy at both circular and linear polarization. However, more detailed information of the different canopy parameters' effects on the final bistatic scattering are simulated here and in this paper, we only concentrate on the circular polarization.

Input parameter settings
The signal frequency range of the MIMICS model (Ulaby et al. 1988) is 0.5-10 GHz, and the incident angle θ has to be greater than 10°. The input information of the canopy and surface layers are shown in Table 2, and the dielectric constants of the soil, trunk, and branches are shown in Table 3.

Bistatic scattering under different observation geometries
In this section, we simulate the direction of the backscattering plane, specular scattering plane, specular cone direction, and perpendicular plane. In these directions, the trunk layer has the highest scattering influence. In other directions, the trunk layer attenuates most of the incident energy. The geometry of the specular cone is shown in Fig. 3. Figure 4 simulates the different frequencies on the P (0.25-0.5 GHz), L (0.5-1.5 GHz), C (4-8 GHz), and X (8-12 GHz) bands and the total scattering with different polarizations (RR, LR, VR, and HR) (Wu and Jin 2014). The incident angle used in Fig. 4a is equal to the scattering angle θ i = θ s , ranging from 10° to 85°, and the scattering azimuth angle is equal to 180°, which is backscattering. In Fig. 4b, θ i = 30°, the range for θ s is from 10° to 85°, ϕ s = 120º, with a bistatic scattering. In Fig. 5c, the incident angle is equal to the scattering angle, and the range is from 10° to 85°, ϕ s = 0º, which corresponds to specular scattering. In Fig. 5d, θ i = θ s , the range is from 10° to 85°, ϕ s = 90º, which corresponds to the perpendicular plane.
In Fig. 4b, when θ i = θ s = 30°, a scattering peak appears, which is caused by the scattering of the tree trunk. At the other angles, the scattering energy is relatively low, which is mainly caused by the scattering of the tree crown and the ground layer. In Fig. 4a, the total scattering energy in the X-band is higher at a large scattering angle, and lower at a small incident angle. In the LR, VR, and HR polarizations, the relationship between the total scattered energy, and the frequency is not obvious. This is mainly due to the different contributions of the canopy and ground layers for different frequencies. The intensity of the volume-scattering inside the canopy shows also different. In Fig. 4b, when θ i = θ s = 30° and the trunk layer Fig. 3 The geometry of the specular cone. While the incidence angle is given by Ki, the scattering azimuth angles given by Ks range from 0° to 360°, thus forming a specular cone attenuates the scattering, the main contribution comes from the canopy and the ground layers. In this figure, the C-band scattering is strongest for all polarizations, while the P-band has the weakest scattering, although for low incidence angles (θ s < 30°), the X-band seems to be the lowest at the VR and HR polarizations. Moreover, we can observe a relationship between L-band and X-band for the scattering to the incident angle. In Fig. 4c, the scattering of the different bands is more obvious, showing that the scattering value of X-band is the largest, followed by P, L, and C bands. Figure 4d shows the scattering properties at the perpendicular plane. From the simulations, we can see that when the scattering angle is lower than 75°, the difference between the X-band and other bands is very obvious. This is mainly due to different frequencies and penetration depth, resulting in different scattering energy levels in the canopy, trunk, and ground layers. Figure 5 shows a comparison of the contributions of various scattering mechanisms to the total scattered energy in the P and X bands (Wu et al. (Under review)). This comparison can better reveal the reasons for the various phenomena in Fig. 4. Figure 5a and b shows the contributions to the total scattering at the backscattering plane. At lower frequencies (e.g., P-band), both the trunk layer and the specular-ground component dominate the total scattering. However, for the X-band, Fig. 5b shows that the influence of the specular-ground component decrease especially for large backscattering angles, where the total scattering is dominated by the trunk layer. In this figure, the crown layer and the direct-ground component are very low and do not contribute to the total scattering. Figure 6a and b shows the contributions at the bistatic scattering geometry (θ i = 30°, ϕ s = 120°) for the P and X bands, at scattering angles from 10° to 85° (Wu et al. (Under review)). In the specular direction, the trunk layer dominates the total scattering both at P and X bands. For the other scattering geometries, not only the crown layer, but also the direct-ground component dominates the total scattering due to the longer wavelength penetration. For the X-band, only the volume-scattering from the crown layer contributes to the total scattering, while the influence of ground and the trunk layers disappear, for the scattering geometry but the specular direction. Figure 7a and b shows the contributions to the total scattering at the specular scattering plane (Wu et al. (Under review)). For the band P, the direct-ground component and the trunk layer dominate the total scattering. In this figure, the contribution of the trunk layer is a lager, while the influence of the crown layer volume-scattering is smaller and does not dominate the contribution to the total scattering. As for the X-band, the components that dominate depend on the scattering geometry. For smaller specular scattering angles, the direct-ground component is the largest contribution, but for larger specular scattering angles, the trunk layer dominates the total scattering. The volume-scattering from the crown layer and the specular-ground component shows a smaller magnitude. Figure 8a and b shows the contributions to the total scattering in the perpendicular direction (Wu et al. (Under review)), where the trunk layer scattering is very large and dominates the total scattering. At the P-band, the specular-ground influence is also very large and is one of the dominant components. However, for the X-band, the specular-ground component becomes weak and is not the dominant contribution. At the P-band, the volume-scattering contribution from the crown is larger than that in the X-band, which is larger due to the penetration properties.  Coherent and non-coherent parts of the scattering for different polarizations. The geometry for the coherent scattering is θ i = θ s =30°, and ϕ s ranges from 0° to 180°. For the non-coherent scattering, the geometry is θ i = 10°, ϕ s = 50°

Comparison between coherent and non-coherent components
Although it has always been thought that the energy received by GNSS-R receivers is mostly from the coherent scattering (especially in the above-mentioned CYG-NSS studies), as the microwave SoOp passes through the vegetation canopy, the energy scattering not only composes from the coherent component, but also from the non-coherent component. Therefore, it is worth to study the scattering for vegetation canopies. Figure 9 shows the comparison between coherent and non-coherent scattering for RR、LR、VR and HR polarizations. In this figure, the geometry for coherent scattering is θ i = θ s = 30°, and ϕ s ranges from 0° to 180°. Concerning the non-coherent scattering, the geometry is θ i = 10°, ϕ s = 50°, the range for ϕ s is the same as that for the coherent scattering. In this figure, we can see that for all the polarization, the coherent scattering is much larger than the non-coherent scattering. In Fig. 10a, the coherent scattering is shown for the different polarizations. The polarization LR is the largest from ϕ s = 0 to ϕ s = 75°. When the scattering azimuth angle is larger than 75, the HR polarization shows the largest scattering. The scattering from the RR polarization is the smallest when the scattering azimuth angle is below 75°. After this value, the VR polarization shows the smallest scattering. Figure 10b shows the noncoherent scattering, where the scattering geometry is θ i = 10°, ϕ s = 50°, and the range ϕ s goes from 0° to 180°. In this figure, we can see that the LR polarization shows the largest scattering, and the RR polarization has the smallest scattering. We can observe that the scattering values between VR and HR polarizations highly depend on the scattering azimuth angle.

Canopy parameters' effects on the bistatic scattering
In this section, we show the effects of different trunk and branch diameters, tree density, and vegetation water content to the total scattering for the different polarizations.

Trunk and branch diameter
In this experiment, we simulate 4 stands with different trunk and branch diameters. We set the parameters given in Tables 1 and 2. In this way, we configure 4 aspen stands with different biomass. We simulate the corresponding bistatic scattering at the band P for the different polarizations. Table 3 lists the diameters for the 4 aspen stands. Figure 11 shows the bistatic scattering at the band P at the different polarizations for the 4 stands shown in Table 3. The scattering geometry is θ i = θ s , and the scattering azimuth angles range from 0° to 180°. The direction ϕ s = 0° corresponds to the specular direction and ϕ s = 180° to the backscattering direction. From our simulation, we can see that different diameters of stands provide different bistatic scattering at these two directions. Therefore, we can clearly distinguish the stands' biomass, since different trunk diameters result in different bistatic scattering. Note that the variable branch diameter has a small effect on the bistatic scattering. This indicates that trunk    diameter dominates the variable contribution to the total scattering, together with the ground layer interaction. From our simulation, we see that larger biomass does not result in higher scattering. For some scattering azimuth angles, the diameters difference cannot be distinguished.

Tree density
In this experiment, we simulate 4 forest scenarios with different tree density. Table 4 shows the 4 different aspen stands with different tree density. For Aspen 1, Aspen 2, Aspen 3, and Aspen 4 stand, the tree densities are 20,000 trees/ha, 10,000 trees/ha, 700 trees/ha, and 500 trees/ha, respectively. Denser tree density corresponds to higher biomass. Figure 12 shows the bistatic scattering at the different polarizations in the band P for the 4 densities. The scattering geometry is the specular cone surface, θ i = θ s = 45°. In this figure, we can see that higher biomass density does not result in larger bistatic scattering. At the RR polarization we can see small differences for different tree densities, while for the LR, VR, and HR polarizations, high tree density strongly varies the the bistatic scattering. We can also observe small differences at the LR and VR polarizations for larger ϕ s . Figure 13 shows canopy scattering contributions of the band P along with the scattering azimuth angles at the different polarizations for the 4 tree densities. In this figure, the total scattering is dominated by the trunk layer. For all the polarization, the scattering azimuth angle increases inversely to trunk layer scattering.

Vegetation water content effects
In this experiment, we investigate the bistatic scattering at the P and L bands which results from 3 scenarios set up with different vegetation water content. When the vegetation water content of a canopy changes, the  permittivity of the scattered energy also changes. The reason is that the dielectric properties of the forest depend on the vegetation water content, and this results in the variation of the bistatic scattering. Table 5 lists the 3 scenarios with details on the different vegetation water content and permittivity of branches and trunks. These 3 different aspen stands are used for the inputs in the experiment. Figure 14 shows the bistatic scattering of various polarizations at the P and L bands for the 3 different scenarios of vegetation water content listed in Table 5. The scattering geometry for the simulations set θ i = θ s = 45°, and ϕ s ranges from 0° to 180°. For the RR polarization at the P-band, the differences for the 3 scenarios are more obvious at large scattering azimuth angles ( ϕ s > 100°). We can observe very small differences in the LR, VR, and HR polarizations. At the HR polarization, some can be seen for large scattering azimuth angles. The range of change is very small, not exceeding 3 dB, showing that the specular cone surface is not suitable for monitoring forest water content.

Conclusions
In SoOp-R, the transmitter and the receiver configuration belongs to the typical configuration of bistatic radar, which is very different from the traditional monostatic radiometer or monostatic radar. In addition, SoOp-R can provide a large number of angle measurements at different observation geometries. In this manuscript, we have presented a theoretical and simulation study on forest canopy monitoring using the MIMICS model (Ulaby et al. 1988;Wu et al. (Under review); Wu and Jin 2014;Liang and Pierce 2005), and revealed the relationship between the canopy parameters and the receiver observations under different observation geometries. Considering the circular polarization characteristics of microwave SoOp remote sensing, we have employed a polarization synthesis technique to calculate various polarization combinations (fully polarized bistatic scattering model). Then, we have simulated and analyzed the corresponding bistatic radar scattering for different polarization conditions. We also have compared and analyzed the contribution of various scattering mechanisms to the total scattered energy, and have revealed differences in scattering characteristics at different frequencies and for various polarizations. The results show that the coherent scattering at the specular cone is larger than the non-coherent scattering, while trunk-dominated forest canopy has strong scattering at the aforementioned different directions. Moreover, we have simulated and analyzed the influence of different canopy parameters to the total scattering energy, and compared and analyzed the difference of the corresponding bistatic scattering characteristics for different (1) trunk and branch diameters, (2) forest canopy density of trees, and (3) vegetation water content, showing strong dependence on the final bistatic scattering observation. It should be noted that the polarization GNSS-R is a potential aspects that should be dug in the future development since the scattering properties at different polarizations are simulated in detail. This will give a theoretical tool for the future GNSS-R vegetation remote sensing study. The theoretical simulations in this paper provide a solid theoretical support for the development of Microwave SoOp remote sensing. Future research is addressed but not restricted to evaluate and verify the corresponding model experiments with measurements data.