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

Wuetal. Geosci. Lett. (2021) 8:25 https://doi.org/10.1186/s40562-021-00195-7 RESEARCH LETTER Forest canopy scattering properties withsignal ofopportunity reflectometry: theoretical simulations Xuerui Wu1,2,6, Andrés Calabia3,6, Jin Xu4,6* , Weihua Bai5,6 and Peng Guo1,2,6 Abstract 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. Keywords: Bistatic scattering, Forest canopy, Radiative transfer equation model, GNSS-reflectometry (GNSS-R), Polarization, Signal of opportunity-reflectometry (SoOP-R) © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. 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 GNSSR technology can provide continuous observations over the Earth’s surface, and surely will become a powerful complement to the traditional microwave techniques. Open Access *Correspondence: [email protected].cn 4 Maritime College, Guangdong Ocean University, Zhanjiang 524088, China Full list of author information is available at the end of the article Page 2 of 15 Wuetal. Geosci. Lett. (2021) 8:25 GNSS-R has been used in ocean remote sensing due to the relative uniformity of the ocean medium, which allows neglecting the polarization characteristics (Ruf etal. 2019; 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 (Eroglu etal. 2019; 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. Theory andmethods 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 MIMICS 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 Page 3 of 15 Wuetal. Geosci. Lett. (2021) 8:25 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 Fig. 1 Bistatic scattering geometry in SoOp-R. A GNSS satellite on the left and a CYGNSS satellite on the right Fig. 2 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 Z1 = d and the trunk layer depth Z2 = Ht Page 4 of 15 Wuetal. Geosci. Lett. (2021) 8:25 the trunk layer is treated as cylinders of uniform diameter. To simplify the calculation, we assume that the incident azimuth is ϕi i = 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 Is = (θs, ϕs ) and the intensity of the incident energy Ii = (θi, ϕi ) can be related through the modified Mueller matrix Lm: 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 Lm 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 sk,(θk,ϕk) , and f(sk;θk,ϕk) , respectively. Nk is the density of the scatterer. The extinction matrix can be expressed as: (1) I s(θs,ϕs)= 1 r 2Lm(θs,ϕs;θi,ϕi;θk,ϕk)Ii(θi,ϕi) . (2) L m=          |Svv| 2 |Svh| 2 ℜ�S ∗ vhSvv�ℑ�S ∗ vhSvv� |Shv|2|Shh|2ℜ�S∗ hhShv�ℑ�S∗ hvShh� 2ℜ�SvvS∗ hv�2ℜ�SvhS∗ hh�ℜ�SvvS∗ hh +SvhS∗ hv�−ℑ�SvvS∗ hh −SvhS∗ hv� 2ℑ � SvvS∗ hv� 2ℑ � SvhS∗ hh� ℑ � SvvS∗ hh +SvhS∗ hv� ℜ � SvvS∗ hh −SvhS∗ hv�          � η (3) Is(µ,ϕs)=T(µ,ϕs)I0(µ0,ϕ0), (4) P (θs,ϕs;θi,ϕi)=Nk  f(sk;θk,ϕk)Lm(θi,ϕi;θs,ϕs;θk,ϕk ) dskdθkdϕk . In this equation,  Spqk (θi,ϕi;θs,ϕs;θk,ϕk) k is the average scattering amplitude coefficient, K indicates the type of the scatterer, pq is the polarization, and k0 is the free space wavenumber. (5) k =            − 2 ℜ( M vv) 0 −ℜ( M vh)−ℑ( M vh) 0−2ℜ(Mhh)−ℜ(Mhv)ℑ(Mhv) −2ℜ(Mhv)−2ℜ(Mvh)−ℜ(Mvv +Mhh)ℑ(Mvv −Mhh) 2ℑ(M hv )−2ℑ(M vh )−ℑ(Mvv −M hh )−ℜ(Mvv +M hh )            , (6) Mpq =  K k=1 i2πN k k0 Spqk (θi,ϕi;θs,ϕs;θk,ϕk) k . 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: (7) T(µ,ϕ s ) = T SG + T GCG + T CG + T GC + T DC + T TG + T GT + T RG. Page 5 of 15 Wuetal. Geosci. Lett. (2021) 8:25 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: (8) T SG(µ,φs)=e−k + cd/µe−k + tH/µR(µ) e−k− tH/µe−k− cd/µ δ(µ − µi)δ(φ − φi). (9) T CGC (µ,φs)= 1 µ e−k+ cd/µe−k+ tH/µR(µ)e−k− tH/µACGC (−µ,φ,µ0,φ0 ) e−k+ tH/µ0R (µ0) e−k− tH/µ0e−k− cd/µ0. (10) T CG(µ,φs)= 1 µ e−k+ cd/µe−k+ tH/µ R (µ) e−k− tH/µACG ( − µ , φ , µ0 , φ0). (11) T GC (µ,φs)= 1 µ AGC (µ,φ,µ0,φ0)e−k+ tH/µ0R(µ0 ) e −k− tH/µ0 e −k− cd/µ0 . (12) T DC (µ,φs)= 1 µ ADC (µ,φ,−µ0,φ0) . (13) T TG (µ,φs)= 1 µ e−k+ cd/µe−k+ tH/µ R (µ) A TG ( − µ , φ ,− µ0 , φ0) e−k− cd/µ0 δ(µ − µ0). 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) T GT (µ,φs)= 1 µ e−k+ cd/µATG (−µ,φ,−µ0,φ0)R(µ0)e−k− tH/µ 0 e−k− cd/µ0 δ(µ − µ0) . (15) T RG(µ,φs)=e−k+ cd/µe−k + tH/µG(µ,φ,−µ0,φ0 ) e −k− tH/µ0 e −k− cd/µ0 . (16) σrt (ψ r ,χ r ,ψ t ,χ t )=4π ˜ Y r m M m Y t m Page 6 of 15 Wuetal. Geosci. Lett. (2021) 8:25 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). However, since the bistatic form of the MIMICS model has not been developed yet, here we upgrade the MIMICS model to the bistatic form and validate the results with the references (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 Page 7 of 15 Wuetal. Geosci. Lett. (2021) 8:25 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, Fig. 4 Canopy bistatic scattering versus scattering angles at P, L, C, and X bands at a backscattering plane, b θi = 30º and ϕs = 120 º, c specular scattering angle, and d incident and scattering angle, θi = θs, ϕs = 90º Page 8 of 15 Wuetal. Geosci. Lett. (2021) 8:25 Fig. 5 a P and b X-band canopy scattering contributions versus scattering angle from aspen at the backscattering plane for various polarizations Fig. 6 a P and b X-band canopy scattering contributions versus scattering angle from aspen at θi = 30°, ϕs = 120° for various polarizations Page 9 of 15 Wuetal. Geosci. Lett. (2021) 8:25 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. Fig. 7 a P and b X bands canopy scattering contributions versus scattering angle at the specular plane for various polarizations