Leaf area index estimation from the time-series SAR data using the AIEM-MWCM model

ABSTRACT

The leaf area index (LAI) inversion accuracy of the water cloud model (WCM) based on the SAR images is low. The main reason is that the WCM model assumes two independent contributions from plant and soil surfaces, without taking into account their interactions. Based on this, this study proposes an efficient LAI estimation method that combines the advanced integral equation model (AIEM) and modified water cloud model (MWCM). Images collected by GaoFen-3 synthetic aperture radar (GF-3 SAR) in the Xiangfu area in the east of Kaifeng City, Henan Province, are processed with a modified SAR-BM3D method and used as test data of the proposed method. The proposed AIEM-MWCM method that combines the AIEM and MWCM models is employed to perform the remote sensing inversion of winter wheat LAI throughout the growth cycle. The results show that the fitting accuracy of winter wheat LAI in the five growth stages achieved by the proposed AIEM-MWCM method inversion is better than that of the Dubois-MWCM model. The R² value of the proposed method is higher than 0.8, and its root mean squared error (RMSE) is lower than 0.3.

KEYWORDS:

• Leaf area index
• synthetic aperture radar
• advanced integral equation model
• modified water cloud model
• AIEM-MWCM model

1. Introduction

Leaf area index (LAI) is an important parameter that characterizes the energy transmission and exchange between soil, crops, and the atmosphere in farmland systems (Anderson et al. 2004; Jarlan et al. 2008). It is also a key index for measuring the biology and biochemistry of the vegetation canopy (Inge et al. 2004). In this paper, the LAI is defined as half of the sum of the total leaf area per unit surface area (Ren et al. 2018). Achieving rapid and accurate inversion of LAI is of great research significance for agricultural monitoring and biogeochemical cycles.

Traditional LAI ground survey techniques, including hemispherical photography, using optical sensors with the LAI-2000/2200, and terrestrial light detection and ranging (LIDAR) scanning, are typically time-consuming, labor intensive, and often destructive to vegetation, which limits their application and makes them unsuitable for high temporal and spatial frequency measurements (Inge et al. 2004; Yan et al. 2019). Moreover, capturing high spatial variability, which is often related to plant growth, via ground measurement is challenging. Remote sensing technology represents an effective method for regional or global crop LAI monitoring due to its advantages of a large coverage area and high efficiency. Based on optical remote-sensing images, these sensors are susceptible to meteorological conditions, such as clouds, rain, and fog, which not only reduces the inversion accuracy of crop LAI but also increases the uncertainty of the inversion result.

Since the wavelengths in the microwave region of the electromagnetic spectrum are less affected by the atmosphere, synthetic aperture radar (SAR) sensors, such as RADARSAT-1, Envisat-ASAR, and GF-3, can successfully image the Earth even when clouds are present. The main methods for LAI remote sensing inversion from the SAR remote sensing images include the empirical and water cloud models. The empirical model uses the backscattering coefficients of the horizontal transmitted and horizontally received (HH), horizontally transmitted and vertically received (HV), and vertically transmitted and vertically received (VV) polarization modes to establish a correlation with the measured LAI (Gao, Niu, and Wu 2010; Kurosu, Fujita, and Chibae 1995; Mansaray et al. 2017; Toan et al. 1989). Several studies have examined the relationship between the empirical model’s sensitivity and LAI for different frequencies and polarizations (Sonia et al. 2019; Pauline et al. 2014; Mandal et al. 2019). Although empirical models are easy to use, the transferability of the resulting statistical relationships usually faces certain constraints when data are collected under different conditions or using sensors with different characteristics.

The water cloud model (WCM) is a semi-empirical parameter inversion model first proposed by Attema and Ulaby in 1978 (Attema and Ulaby 1978). Recent research has examined the wide application of the WCM in numerous subsequent inversion studies on soil moisture and vegetation growth parameters (Emilie et al. 2015; Mehdi et al. 2015; Tao et al. 2016; Yadav, Prasad, and Bala 2021; Yang et al. 2016). The results have suggested that the WCM approach can be suitable for the operational monitoring of crops. However, the WCM involves many model parameters that are difficult to obtain. The WCM attempts to use straight forward vegetation parameters to represent the canopy scattering and attenuation components in the model (Shubham et al. 2022; Verma et al. 2022), which leads to an low accuracy problem of model inversion. Furthermore, the WCM ignores the mutually influencing interaction between plant and soil surfaces, which has a substantial influence on the total scattering and attenuation properties of the land surface.

In summary, The LAI inversion accuracy of the WCM based on the SAR images is low. The main reason is that the WCM model assumes two independent contributions from plant and soil surfaces, without taking into account their interactions. When the plant density is considerable, this might cause a large departure from the true scattering and attenuation properties of the ground surface. By using the GF-3 SAR images processed by the modified SAR-BM3D method, this study combines the advanced integral equation model (AIEM) and modified water cloud model (MWCM) and proposes an AIEM-MWCM inversion method for LAI in five growth stages of winter wheat. The proposed method combines the AIEM and MWCM models and makes full use of their respective advantages to invert surface parameters, which can improve inversion accuracy.

2. Materials and methods

2.1. Materials

2.1.1. Study area

The study area used in this work was the Xiangfu area (34°30′–34°56′ N, 114°07′–114°43′ E) in the eastern part of Kaifeng City, Henan Province, which belongs to the Yellow River alluvial plain. This is a typical winter wheat growing area and an important food-producing region. It has a size of 5 km × 5 km, an average altitude of 62.5–89.3 m, an average annual temperature of 14°C, an average annual precipitation of 628 mm, and an average frost-free period of 214 days. A temperate continental monsoon climate provides sufficient sunshine, suitable temperatures, and abundant precipitation, which are conducive to crop growth. The main crops in the study area include winter wheat, rice, and corn, but according to the growing period of crops, the main crop is winter wheat. The location of the study area and distribution of sampling points are shown in Figure 1.

Figure 1. Location of the study region and distribution of sampling points. (a) Location of Xiangfu area in Henan Province; (b) location of leaf area index sampling points.

2.1.2. Satellite data collection and processing

In this study, the GF-3 SAR data collected in the study area were downloaded from the China Centre for Resources Satellite Data and Application (CCRSDA) and used to develop the proposed model. The GF-3 SAR data were in the C band, the incidence angles ranged from 20° to 41°, and the transit days included March 2, March 21, April 7, May 28, and June 9, 2021. The operating frequency was 5.4 GHz, and the spatial resolution was 8 m. The images were processed by performing full polarization introduction, radiometric calibration, amplitude enhancement, and speckle filtering, and three scattering components were extracted using the Freeman-Durden polarimetric decomposition of the processed data. By using the TanDEM-X elevation model data (DEM), the GF-3 SAR data were calibrated to sigma nought using the calibration factors and sine of the local incidence angle obtained from the TanDEM-X DEM. To reduce speckle noise, this study performed speckle filtering using a BoxCar 5 × 5 kernel filter. The sigma-nought radiometric calibration was performed using the absolute calibration coefficients downloaded from CCRSDA (http://www.cresda.com/CN/).

The MODIS LAI product (MOD15A2) was downloaded from the United States Geological Survey (http://lpdaac.usgs.gov/). The orbital number was h26v05. The spatial resolution of the product was 500 m, and the temporal resolution was 8 days. The downloaded data were trans-projected and subjected to clipping to obtain the LAI images of the study area.

2.1.3. Field data observation

The aim of the experiment was to prevent changes in weather conditions from affecting the accuracy of the measured LAI. The LAI inversion of winter wheat in five growth stages was analyzed; the selected acquisition dates were March 1, March 24, April 20, May 27, and June 03, 2021, and they were almost synchronous with the satellite transit time. The daily acquisition interval was from 07:30 am to 09:30 am Fifty winter wheat quadrats with a size of 16 m × 16 m were used in the experiment, and three spots with the same growth status were selected. All quadrats were at least 30 m away from trees, buildings, and roads. During data sampling, sunlight was shielded, and LAI-2000 was used for measurement. In particular, a 180° covering cap was used to ensure that there were no measurement errors. In addition, the distance between the measured ground and the instrument was adjusted to 5 cm according to the standard to ensure that there was no error between the winter wheat canopy and the viewing angle range of the measurement instrument. Further, to avoid LAI measurement errors, the LAI of each sample site was measured four times, and the average value was used as the end value. The measured LAI data were divided according to the ratio of 3:2 into training and verification sets, i.e. 60 training datasets were selected for model inversion and the remaining 40 verification datasets were used for model verification. Furthermore, the longitude and latitude values of each winter wheat sampling point were recorded in real-time using a U.S. Global Positioning System (GPS) locator. The measurement process of winter wheat LAI and sampling point locations are presented in Figure 1.

2.2. Methods

The proposed method for LAI estimation uses the AIEM-MWCM model that is developed using the GF-3 satellite data, as shown in Figure 2. First, the modified SAR-BM3D was used to process the GF-3 satellite data by performing fully polarized guidance, radiometric calibration, amplitude enhancement, and speckle filtering to improve image quality. Then, Dubois and AIEM were used to calculate the soil backscattering coefficient and the radar vegetation index, respectively. Finally, the Dubois-MWCM and AIEM-MWCM models were constructed to retrieve the LAI values of winter wheat during the entire growth period.

Figure 2. Flowchart of leaf area index inversion by AIEM-MWCM.

2.2.1 Modified SAR-BM3D

The block-matching stage of the SAR-BM3D method aims to find the image block that is most similar to the reference block (Giovanni et al. 2017). The similarity measure based on the actual noise distribution is relatively single, and the filtered images can miss some details and generate unnecessary false information. Based on the SAR-BM3D method, the Pearson’s correlation coefficient is introduced to improve the similarity measure of the reference and matching blocks to improve the accuracy of image block matching and grouping. After the similarity measure is improved, the coordinated filtering and aggregation processing are performed to obtain the basic and final estimations, and the SAR image after noise reduction is obtained.

By matching and grouping a SAR image i with noise, the reference block with a size of $n\times n$ is selected, and the similarity measure $d_1(a,b)$ is calculated as follows: (1) $d_1(a,b)=(2L-1)\sum _{i=1}^{n^2}\log \left({\displaystyle \frac{a(i)}{b(i)}+{\displaystyle \frac{b(i)}{a(i)}}}\right).$(1) where $L$ is the apparent parameter.

The introduced Pearson’s correlation coefficient reflects the correlation degree between variables, and $|\rho |\le 1$; the Pearson’s correlation coefficient is calculated as follows: (2) $\rho ={\displaystyle \frac{\mathrm{cov}(a,b)}{\sqrt{{\delta }_a}\times \sqrt{{\delta }_b}}.}$(2) where $\mathrm{cov}(a,b)$ represents the covariance of the reference block a and matching block b; ${\delta }_a$ and${\delta }_b$ represent the variances of image blocks a and b, respectively.

When searching for similar blocks using the SAR-BM3D method, the smaller the similarity measure$d_1(a,b)$ is, the higher the similarity between the two image blocks is. Also, the larger the absolute value of $\rho$ is, the higher the degree of positive (negative) correlation will be, indicating a higher similarity between the two images. Therefore, the similarity measure is redefined as follows: (3) $d_1^{\mathrm{final}}(a,b)=(1-{\rho }^2)d_1(a,b).$(3) According to the preset similarity threshold $\tau$, M similar image blocks are searched, and a set of similar blocks is formed as follows: (4) $B_a=\{i\in I:d_1^{final}(a,b)\le \tau \}.$(4)

The similar blocks are integrated into a three-dimensional (3D) array $G_a$, and the tasks of 3D transformation, hard threshold processing, and 3D inverse transformation are performed to obtain the estimated value of the set of similar blocks as follows: (5) $R_a=T_{3D}^{-1}(X_{\gamma }(T_{3D}(G_a))).$(5) where $X_{\gamma }$ represents the wavelet hard threshold, and $T_{3\mathrm{D}}$ and $T_{3\mathrm{D}}^{-1}$ denote the 3D transformation and inverse transformation of a 3D array, respectively.

A pixel i can appear in several different combinations of similar blocks. The basic estimation image is obtained from the weighted average of the overlapping image blocks, given by (6) $x_{\mathrm{basic}}(i)={\displaystyle \frac{{\sum }_{i\in I}{\sum }_{i_m\in M}{W}_{\mathrm{basic}}{R}_{i_m}(i)}{{\sum }_{i\in I}{\sum }_{i_m\in M}{W}_{\mathrm{basic}}{P}_{i_m}(i)}.}$(6) where $R_{i_m}(i)$ represents the estimated value of an image block $i_m$ containing pixel i obtained after the 3D inverse transformation; $P_{i_m}(i)$ represents the characteristic function of a similar block; $W_{\mathrm{basic}}$ represents the basic estimated weight, and $W_{\mathrm{basic}}=\{\begin{array}{c}{({\sigma }^{2}N)}^{-1},N\ge 1\\ 1,N<1\end{array}$; $N$ denotes the number of all non-zeros retained after hard threshold filtering.

Before performing similar block matching on the original image, it is necessary to consider data obtained in the basic estimation stage and use the position information to perform block matching on the noisy image. According to the improved similarity measure method described in Equations (1)–(6) above, the improved similarity measure $d_2^{\mathrm{final}}(a,b,c,d)$ is expressed as follows: (7) $d_2(a,b,c,d)=\sum _{i=1}^{n^2}(2L-1)\log \left({\displaystyle \frac{a(i)}{b(i)}+{\displaystyle \frac{b(i)}{a(i)}}}\right)+L{\displaystyle \frac{{|c(i)-d(i)|}^2}{c(i)d(i)}}$(7) (8) $d_2^{\mathrm{final}}(a,b,c,d)=(1-{\rho }^2)d_2(a,b,c,d)$(8)

The similarity block grouping result of an image is composed of a 3D array $G$ formed by the basic estimation and a 3D array formed by the original noisy image. The two 3D arrays are transformed into the 3D domain, Wiener filtering is performed on the 3D array$G_W$, and the estimated value of the set of similar blocks is obtained by 3D inverse transformation, which can be expressed as follows: (9) $R_W=T_{3\mathrm{D}}^{-1}(W(T_{3\mathrm{D}}(G_W)))$(9) where $W(\cdot )$ represents the Wiener filtering function, and $W={\displaystyle \frac{{|T_{3\mathrm{D}}G|}^2}{{|T_{3\mathrm{D}}G|}^2+{\sigma }^2}{T}_{3\mathrm{D}}}$ and $T_{3\mathrm{D}}^{-1}$represent the 3D transformation and 3D inverse transformation, respectively.

Using the same method as for the basic estimation, the final estimation $Q_{\mathrm{final}}(i)$ is obtained by the weighted average, and then, the original position is restored to obtain the denoised image as follows: (10) $x_{\mathrm{final}}(i)={\displaystyle \frac{{\sum }_{i\in I}{\sum }_{i_m\in M_W}{W}_{wie}{R}_{i_m}^W(i)}{{\sum }_{i\in I}{\sum }_{i_m\in M_W}{W}_{wie}{P}_{i_m}^W(i)}}$(10) where $M_W$ represents the number of image blocks in the similar block set; $R_{i_m}^W(i)$ is the estimated value of the image block $i_m$ containing pixel i obtained after the 3D inverse transformation; $P_{i_m}^W(i)$ denotes the characteristic function of a similar block; $W_{wie}$ represents the final estimated weight, and$W_{wie}={\sigma }^{-2}{W}_2^{-2}$; the weighted weight depends on the Wiener filter coefficient$W$.

Based on the methods presented above, the equivalent number of loops (ENL) and edge preservation index (EPI) were used in the experiment to evaluate the effect of SAR image noise reduction. The equivalent visual number was used to evaluate the degree of suppression of speckle noise in SAR images. Generally, the larger the value of the equivalent visual number was, the better the suppression effect of speckle noise was. The equivalent apparent number was calculated as follows: (11) $ENL=\alpha ({\mu }^2/{\sigma }^2)$(11) where $\mu$ and $\sigma$ represent the mean and standard deviation, respectively; $\alpha$ is related to the SAR image format, and for the intensity format, $\alpha =1$, whereas for the amplitude format, $\alpha =4/\pi -1$.

In addition, the edge retention index was used to evaluate the edge retention ability of filtered images. Normally, the closer its value is to one, the better is the edge information retention ability of a method. The edge retention index was calculated as follows: (12) $EPI={\displaystyle \frac{\sum |p_s(i,j)-p_s(i-1,j+1)|}{\sum |p_o(i,j)-p_o(i-1,j+1)|}}$(12) where $p_o(i,j)$ and $p_s(i,j)$ represent the image pixels before and after filtering, respectively.

2.2.2. Bare soil backscatter modeling by Dubois

Dubois, Jakobl, and Engman (1995) established a semi-empirical SM retrieval model for a bare soil area. The Dubois empirical model was used on the bare surface, and f∈[1.5, 11 GHz], $k_s$∈[0.34, 2.5], and $\theta$∈[30°, 60°]; the retrieval effect of the root mean square height and soil water content was the best. The Dubois model can be conceptually described as follows: (13) ${\sigma }_{HH}^o={10}^{-2.75}{\displaystyle \frac{{\cos }^{1.5}\theta }{{\sin }^5\theta }{10}^{{0.028}_{{\epsilon }_r\tan (\theta )}}(k_s\sin \theta {)}^{1.4}{\delta }^{0.7}}$(13) (14) ${\sigma }_{VV}^o={10}^{-2.35}{\displaystyle \frac{{\cos }^3\theta }{{\sin }^3\theta }{10}^{{0.046}_{{\epsilon }_r\tan (\theta )}}(k_s\sin \theta {)}^{0.7}{\delta }^{0.7}}$(14) where ${\sigma }_{VV}^o$ is the VV polarization backscattering coefficient; ${\sigma }_{HH}^o$ is the HH polarization backscattering coefficient; $f$ is the GF-3 satellite frequency, and in this work, it has a value of 5.4 GHz; $\theta$ is the scattering angle; $k$ is the wavenumber; and $\delta$ is the wavelength.

2.2.3. Bare soil backscatter modeling by AIEM

The surface backscattering coefficient values obtained by the integral equation model (IEM) were not completely consistent with the measured values. Wu et al. (2001) proposed an improved IEM in 2001 and optimized the AIEM in 2004. Because of the complexity of the AIEM, it cannot be directly applied to the automatic operation of soil water content inversion. Therefore, based on the complex theoretical formula, the statistical analysis of the factors affecting the radar backscattering coefficient in the theoretical model, combined with the measured data, was used to construct a semi-empirical model of soil water content; radar backscattering was used to extract soil water content information. Finally, the AIEM was used to generate a bare soil backscatter database for the GF-3 configuration.

The (double-basis) scattering coefficient of the AIEM is expressed as follows: (15) ${\sigma }_{qp}^0={\sigma }_{qp}^k+{\sigma }_{qp}^c+{\sigma }_{qp}^{kc}$(15) where p and q represent the polarization modes of transmission and reception, respectively; ${\sigma }_{qp}^0$ is the scattering coefficient; ${\sigma }_{qp}^k$ is the Kirchhoff field term; ${\sigma }_{qp}^c$ is the compensating field term; and ${\sigma }_{qp}^{kc}$ is the cross term.

By substituting the calculation formulas of Kirchhoff's surface field term, compensation field term, and cross term into (15), the scattering coefficient of the AIEM can be obtained as follows: (16) ${\sigma }_{qp}^0={\displaystyle \frac{k_1^2}{2}{e}^{(-{\sigma }^2(k_z^2+k_{sz}^2))}\times \sum _{n=1}^{\mathrm{\infty }}|I_{qp}^n{|}^2{W}^{(n)}(k_{sx}-k_x,k_{sy}-k_y)}$(16) (17) $I_{qp}^n=(k_{sz}+k_z{)}^2{f}_{pq}{e}^{-s^2{k}_z{k}_{sz}}+(k_{sz}{)}^n{F}_{pq}(-k_x,-k_y)+(k_z{)}^n{F}_{pq}(-k_{sz},-k_{sy})$(17) where $k_1$ is the free space wave number of the transmission medium; $s$ is the root-mean-square height of the soil surface; $W^{(n)}(k_{sx}-k_x,k_{sy}-k_y)$ is the nth-order Fourier transform of the surface correlation function, where $k_z=k\cos {\theta }_i$, $k_{sz}=k\cos {\theta }_s$, $k_x=k\cos {\theta }_x$, $k_{sx}=k\cos {\theta }_s\cos {\phi }_s$, $k_y=k\cos {\theta }_i\cos \phi$, and $k_{sy}=k\cos {\theta }_s\cos {\phi }_s$; $\phi$ is the incidence azimuth; ${\phi }_s$ is the scattering azimuth; and $F_{pq}$ and $f_{pq}$ are functions related to Fresnel reflectance $R_{VV/HH}$. (18) $R_{VV}={\displaystyle \frac{{\epsilon }_r\cos \theta -\sqrt{{\epsilon }_r-{\sin }^2\theta }}{{\epsilon }_r\cos \theta +\sqrt{{\epsilon }_r-{\sin }^2\theta }}}$(18) (19) $R_{HH}={\displaystyle \frac{\cos \theta -\sqrt{{\epsilon }_r-{\sin }^2\theta }}{\cos \theta +\sqrt{{\epsilon }_r-{\sin }^2\theta }}}$(19) where ${\epsilon }_r$ is the dielectric constant of the soil; $l$ is a rough surface with the correlated length; and the roughness power spectrum can be expressed as follows: (20) $W^{(n)}(k_{sx}-k_x,k_{sy}-k_y)={\displaystyle \frac{l^2}{2n}{e}^{{{\displaystyle \frac{2kl\sin \theta }{4n}}}^2}}$(20) where $l$ is the correlation length.

Therefore, the AIEM can be conceptually described as follows: (21) ${\sigma }_{soil}=({\sigma }_{HH}^o,{\sigma }_{VV}^o)=AIEM(f,\theta ,PP,SM,s,l,ACF),$(21) where $\theta$ is the scattering angle; $PP$ is the polarization state (HH and VV polarizations); $SM$ is the input soil moisture of the Dobson model; $s$ is the root-mean-square height; $l$ is the correlation length; and $ACF(\cdot )$ is the adopted exponential autocorrelation function. The ranges (minimum and maximum) of each of the three main model parameters are summarized in Table 1.

Table 1. Value ranges of the advanced integral equation model (AIEM) parameters.

Parameter	Unit	Lower boundary	Upper boundary
$SM$	m^3m^−3	0.02	0.4
$s$	cm	−11	−10
$l$	cm	10	30

2.2.4. Dubois–MWCM and AIEW-MWCM

The MWCM is an important model used to study the relationship between radar backscattering and crop parameters (Zhang et al. 2018), and it is defined as follows: (22) ${\sigma }_0={\sigma }_{veg}(\theta )+{\tau }^2(\theta ){\sigma }_{soil}(\theta ),$(22) (23) ${\sigma }_{veg}(\theta )=A\times M\times \cos (\theta )\times [1-{\tau }^2(\theta )],$(23) (24) ${\tau }^2(\theta )=\exp \left[-2B\times {\displaystyle \frac{M}{\cos (\theta )}}\right],$(24) where ${\sigma }_0$ is the total backscattering coefficient of the surface canopy covered by random polarimetric vegetation; ${\sigma }_{veg}(\theta )$ is the backscattering coefficient of the vegetation layer; ${\sigma }_{soil}(\theta )$ is the direct surface backscattering coefficient; ${\tau }^2(\theta )$ is the double-layer attenuation factor of the radar wave penetrating the crop layer; $\theta$ is the incidence angle of the radar waves; M is the vegetation canopy parameter (i.e. LAI); B is the parameter that is defined by the vegetation type and incident electromagnetic wave frequency.

The model was built on the basis of the radiation transfer model, which considers only the volume scattering term from vegetation reflection and backscattering of the ground after double-layer attenuation of vegetation while ignoring other forms of scattering in vegetation and soil. The radar backscattering coefficient is sensitive to the surface dielectric constant and is closely related to soil moisture. The value of ${\sigma }_{soil}(\theta )$ is obtained by using (13) and (14) and is then substituted into (22)–(24), and the radar vegetation index ($f_{VI}$) is used to construct the Dubois-MWCM model, as shown in (25). Further, ${\sigma }_{soil}(\theta )$ is obtained by using (21) and is then substituted into (22)–(24); $f_{VI}$ is used to construct the AIEW-MWCM model, as shown in (26). (25) $\{\begin{array}{c}{\sigma }_0=A\times M\times \cos (\theta )\times [1-{\tau }^2(\theta )]+{\sigma }_{soil}(\theta )\times \exp [-2B\times M/\cos (\theta )]\\ {\sigma }_{HH}^o={10}^{-2.75}{\displaystyle \frac{{\cos }^{1.5}\theta }{{\sin }^5\theta }{10}^{{0.028}_{{\epsilon }_r\tan (\theta )}}{(k_s\sin \theta )}^{1.4}{\delta }^{0.7}}\\ {\sigma }_{VV}^o={10}^{-2.35}{\displaystyle \frac{{\cos }^3\theta }{{\sin }^3\theta }{10}^{{0.046}_{{\epsilon }_r\tan (\theta )}}{(k_s\sin \theta )}^{0.7}{\delta }^{0.7}}\\ f_{VI}=RV{I}_{Freeman}={\displaystyle \frac{P_V}{P_S+P_V+P_D}}\end{array}$(25) (26) $\{\begin{array}{c}{\sigma }_0=A\times M\times \cos (\theta )\times [1-{\tau }^2(\theta )]+{\sigma }_{soil}(\theta )\times \exp [-2B\times M/\cos (\theta )]\\ {\sigma }_{soil}=({\sigma }_{HH}^o,{\sigma }_{VV}^o)=AIEM(f,\theta ,PP,SM,s,l,ACF)\\ f_{VI}=RV{I}_{Freeman}={\displaystyle \frac{P_V}{P_S+P_V+P_D}}\\ \end{array}$(26) where $P_V$, $P_D$, and $P_S$ are the scattering components of the covariance matrix of the bulk scattering, even scattering, and odd scattering models, respectively; A and B are the Dubois-MWCM and AIEW-MWCM parameters determined by Levenberg–Marquardt nonlinear least square optimization.

The Dubois-MWCM and AIEW-MWCM denote the functions of a number of input parameters that quantify crop growth and soil surface parameters. However, SAR signals are also affected by the crop growth parameters, such as LAI, vegetation water content (VWC), and planting soil properties ($M_S$). The backscattering coefficient of crops is expressed as follows: (27) ${\sigma }_{mod}=f(LAI,f_{VI}).$(27) The polarization backscattering coefficient is mainly related to vegetation scattering. Therefore, the establishment of a lookup table (LUT) under the HH polarization can effectively invert the LAI obtained by the Dubois-MWCM and AIEW-MWCM models. The LUT method is, by nature, a direct comparison of the radiative transfer model backscattering value calculated from the GF-3 SAR images under the HV polarization against the backscattering value simulated by the Dubois-MWCM and AIEW-WCM models. It is a simple and widely used inversion method in the field of vegetation parameter retrieval from remote-sensing data. The LUT consists of two main steps: LUT generation and LUT searching.

In the LUT generation step, 90,000 parameter combinations were randomly generated following uniform distributions. The value ranges (i.e. the minimum and maximum values) of the five main model parameters are presented in Table 2.

Table 2. The range of all parameters in generating look up table.

Name of parameter	Unit	Lower limit	Upper limit	Sources
$\mathrm{LAI}$	m^2/m^2	0.1	7	Ground measurement
${\sigma }_{\mathrm{mod}}$	dB	−11	−9	AIEW-MWCM acquisition
${\sigma }_{\mathrm{mod}}$	dB	−12	−10	Dubois-MWCM acquisition
${\sigma }_{\mathrm{obs}}$	dB	−12	−10	Image acquisition
${\mathrm{f}}_{\mathrm{VI}}$	No unit	0.15	0.45	Image acquisition

In the LUT searching step, the mean square error Z was used to represent the difference between the observed and estimated backscattering values, and it was calculated as follows: (28) $Z=\min \sqrt{{\displaystyle \frac{1}{N}\sum _{i=1}^n{({\sigma }_{obs}-{\sigma }_{mod})}^2}}$(28) where ${\sigma }_{obs}$ is the backscattering value calculated from a GF-3 SAR image under the HV polarization, ${\sigma }_{mod}$ is the backscattering value simulated by the AIEW-MWCM model, and N is the number of pixel points in the image.

2.2.5. Evaluation methodology

The fitting between the simulated and measured LAI values of winter wheat was evaluated through accuracy evaluation to obtain the optimal inversion model. Two common evaluation metrics, R² (Wang et al. 2020) and RMSE (Amin et al. 2021; Jafari and Keshavarz 2021), were selected to evaluate the fitting accuracy, and they are respectively defined as follows: (29) $R^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y-y_i)}^2}{{\sum }_{c=1}^n{(\mathrm{y}-y_c)}^2},}$(29) (30) $\mathrm{RMSE}=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^n{(y_e-y_d)}^2}{n}}}.$(30) where R² is the determination coefficient; y is the measured ecological parameter value; $y_i$ is the estimated ecological parameter value; $y_c$ is the average value of the measured ecological parameter; RMSE is the root mean square error;$y_e$ is the measured LAI value; $y_d$ is the estimated LAI value; and n is the number of samples.

3. Results and discussion

3.1. Improved filtering analysis of SAR-BM3D method

To verify the effectiveness of the SAR-BM3D method, the Frost, Kuan, SAR-BM3D, and improved SAR-BM3D filters were applied to the SAR images of the jointing stage in the measured area, as shown in Figure 3 and Table 3.

Figure 3. Image filtering results obtained by (a) Frost filter; (b) Kuan filter; (c) SAR-BM3D filter; and (d) improved SAR-BM3D filter.

Table 3. Objective and quantitative comparison of filtering results of different methods.

Methods	The evaluation index
	ENL	EPI
Frost filter	1.3756	0.4395
Kuan filter	1.5210	0.3758
SAR-BM3D	1.3244	0.5133
Modified SAR-BM3D	1.7856	0.6988

As shown in Figure 3, the noise reduction effect after filtering by the Frost and Kuan filters was good, but their use made the image too smooth and blurred, resulting in many edge details being lost. The denoising effect of the SAR-BM3D filter was better than that of the Frost and Kuan filters, but the preserving effect of the image edge and texture information was incomplete. The improved SAR-BM3D filter has a good noise reduction effect; the filtered image was relatively clear, the visual effect in the homogenous area was good, and the image edge details were well-maintained.

Further, as presented in Table 3, the equivalent visual numbers of Frost and Kuan filtering were large, but the edge-preserving index was relatively small. The edge-preserving index of the wavelet threshold and SAR-BM3D approach was close to one, but the equivalent visual number was small. The improved SAR-BM3D filtering method was superior to the other filtering algorithms regarding the equivalent sight number and edge retention index. Based on the analysis results of subjective vision and objective quantitative indicators, the improved SAR-BM3D method exhibited good denoising performance and could effectively suppress the noise in SAR images and preserve the image edge details. However, there were still a few artifacts in the noisy area of SAR images with complex details.

3.2. LAI inversion by the RVI_Freeman model

A regression model was constructed between the RVI_Freeman and winter wheat LAI for the entire growth cycle using five types of RVI_Freeman models, and the accuracy was evaluated. The evaluation results of the RVI_Freeman model are presented in Figure 4.

Figure 4. LAI Inversion by the RVI_Freeman model. (a–e) are from the jointing, booting, heading, filling, and maturing stages, respectively.

As shown in Figure 4, the inversion accuracy of the LAI of winter wheat during the entire growth stage obtained by the RVI_Freeman model was low and could not meet the inversion requirements. The inversion accuracy of the mature stage was the highest among all the stages, with R² = 0.5629 and RMSE = 0.6980. The second-best stage regarding the inversion accuracy was the booting stage, with R² = 0.5484 and RMSE = 0.7521. The jointing stage had the least inversion accuracy, with R² = 0.5295 and RMSE = 0.6988. Therefore, in the high-LAI area, the main component of the total scattering power was the volume scattering component, while secondary scattering and surface scattering were weak. However, in the low-LAI area, the radar scattering effect was relatively complex, and it was easy to form strong secondary scattering components, which reduced the sensitivity of the volume scattering component ratio value to the LAI value. In summary, for the radar vegetation index, the spectral saturation problem of winter wheat could be overcome by the penetration of radar to a certain extent.

3.3. LAI inversion by the backscattering coefficient of Dubois-MWCM and AEIM-MWCM

The Dubois-MWCM and AEIM-MWCM models were employed to extract the backscattering values of the HH and VV polarization modes from the GF-3 SAR images, respectively. The LAI observations and predicted values were then compared, as shown in Tables 4 and 5.

Table 4. Accuracy evaluation of the backscattering coefficient model of Dubois-MWCM.

Growth Stage	Model	R^2	RMSE	Growth Stage	Model	R^2	RMSE
Jointing	HH	0.1863	0.8896	Booting	HH	0.1763	0.9324
	VV	0.1563	0.9324		VV	0.1432	0.9423
Heading	HH	0.1768	0.9236	Filling	HH	0.1769	0.9155
	VV	0.1325	0.9411		VV	0.1399	0.9356
Maturing	HH	0.1832	0.9153
	VV	0.1455	0.9367

Table 5. Accuracy evaluation of the backscattering coefficient model of AEIM-MWCM.

Growth Stage	Model	R^2	RMSE	Growth Stage	Model	R^2	RMSE
Jointing	HH	0.1986	0.8769	Booting	HH	0.2156	0.8632
	VV	0.1758	0.9125		VV	0.1866	0.9125
Heading	HH	0.1892	0.8963	Filling	HH	0.1856	0.9011
	VV	0.1639	0.9243		VV	0.1745	0.9586
Maturing	HH	0.1936	0.8936
	VV	0.1745	0.9314

Because the interaction mechanism between radar and vegetation was complex and affected by surface conditions, there was no strong correlation between the two backscattering coefficients and LAI, as shown in Tables 4 and 5. The correlation between the HH polarization backscattering coefficient and the LAI was higher than that of the VV polarization backscattering coefficient because the HH polarization could better reflect the secondary and surface scattering of the surface.

3.4. LAI inversion by Dubois-MWCM

Because the polarization backscattering coefficient was mainly related to vegetation scattering, the parameters of the Dubois-MWCM were determined by the nonlinear least squares method under the HH cross-polarization, as shown in Table 6. The relationship between the measured and predicted LAI values obtained by the Dubois-MWCM under the HH polarization at different stages is presented in Figure 5.

Figure 5. Correlation analysis of Dubois-MWCM and LAI. (a–e) are from the jointing, booting, heading, filling, and maturing stages, respectively.

Table 6. Model parameters at HH polarizations.

Polarization	Growth Stages	Model Parameters
		A	B
HH	Jointing	0.046	0.198
	Booting	0.037	0.210
	Heading	0.057	0.186
	Filling	0.048	0.132
	Maturing	0.036	0.185

As shown in Figure 5, under the HH polarization, the LAI inversion value obtained by the Dubois-MWCM had a strong correlation with the measured LAI. Among the five growth stages, the fitting accuracy between the jointing stage and measured LAI was the highest, having R² = 0.7865 and RMSE = 0.4589, and it was followed by the maturation (R² = 0.7850, RMSE = 0.4534), filling (R² = 0.7754, RMSE = 0.4589), heading (R² = 0.7731, RMSE = 0.4486), and booting stages (R² = 0.7584, RMSE = 0.4651). In the Dubois-MWCM model, the SAR backscattering coefficient was directly used in the winter wheat LAI inversion, and the impact of ground surface reflection was neglected, which led to low LAI inversion accuracies.

Among the five growth stages of winter wheat, the jointing stage had the highest regional LAI inversion accuracy. When the winter wheat jointed and completely covered the ground surface, the accuracy of winter wheat LAI inversion was the highest because the microwave scattering process was simple, and the Dubois-MWCM contained fewer parameters at this growth stage.

3.5. LAI inversion by AEIM-MWCM

The parameters of the AEIM-MWCM model were determined by the nonlinear least squares method under the HH cross-polarization, as shown in Table 7. The AEIM-MWCM model was used to invert LAI by employing the least-squares algorithm and LUT method. The relationship between the LAI observations and the values predicted by the AEIM-MWCM model under the HH polarization is presented in Figure 6.

Figure 6. Correlation analysis of AEIM-MWCM and LAI. (a–e) are from the jointing, booting, heading, filling, and maturing stages, respectively.

Table 7. Model parameters at HH polarizations.

Polarization	Growth Stages	Model Parameters
		A	B
HH	Jointing	0.027	0.143
	Booting	0.038	0.135
	Heading	0.024	0.187
	Filling	0.049	0.203
	Maturing	0.036	0.189

As shown in Figure 6, the AEIM-MWCM model achieved the best inversion of winter wheat LAI at the maturing stage, having an R² of 0.8867 and an RMSE of 0.2989. The inversion accuracy of the stages decreased in the following order: jointing (R² = 0.8720 and RMSE = 0.3524), filling (R² = 0.8695 and RMSE = 0.3726), booting (R² = 0.8535 and RMSE = 0.4213), and heading stage (R² = 0.8519 and RMSE = 0.3395). In summary, the differences between the measured LAI and LAI value estimated from the GF-3 data were small, and the LAI retrieval results were satisfactory. These results indicate that the proposed LAI method for GF-3 data was reliable and that GF-3 data could provide acceptable performance for LAI estimation.

The C-band microwave signal penetrated the wheat layer to the ground surface, and it was reflected by the ground surface and received by the sensor. The microwave energy reflected by the ground surface accounted for a large proportion of the total microwave energy. Unlike the AEIM-MWCM model, the Dubois-MWCM model could effectively use the energy reflected from the crop canopy and ground surface. Therefore, the inversion accuracy of the AEIM-MWCM model was higher than that of the Dubois-MWCM model. Among the five growth stages of winter wheat, the jointing stage had the highest regional LAI inversion accuracy, and the heading stage had the lowest LAI inversion accuracy. When the winter wheat jointed and completely covered the ground surface, the accuracy of winter wheat LAI inversion was the highest because the microwave scattering process was simple, and the submodel contained fewer parameters at this growth stage. In contrast, at the heading stage, winter wheat ears and stem leaves showed obvious stratification, the microwave scattering process was complex, and the submodel contained a large number of parameters; thus, the LAI inversion accuracy of winter wheat was the lowest in the heading stage.

3.6. AEIM-MWCM verification

To verify the reliability of the LAI inversion model, the LAI value obtained by the AEIM-MWCM inversion was studied, and the measured LAI was taken as a true value. The regression analysis was conducted between the remaining 40 true values, other than the modeling samples and the LAI data obtained by the AEIM-MWCM model to obtain the correlation between the LAI inversion value and the true value in each growth stage, as shown in Figure 7.

Figure 7. The AEIM-MWCM model and LAI verification analysis. (a–e) are from the jointing, booting, heading, filling, and maturing stages, respectively.

As shown in Figure 7, the R² values between the LAI of the five growth stages of winter wheat (the jointing, booting, heading, filling, and maturing stages) inverted by the AEIM-MWCM model and the measured LAI were 0.8518, 0.8551, 0.8643, 0.8531, and 0.8609, and the RMSE values were 0.3924, 0.3762, 0.3869, 0.3989, and 0.3864, respectively. The comprehensive analysis results in Figure 7 showed that although the correlation of each growth stage was different, the R² was over 0.85, and the maximum RMSE was 0.3989, indicating that the LAI inversion of winter wheat in the five growth stages performed by the AEIM-MWCM reflected the growth and changes in winter wheat with preferable inversion results. The proposed model that combined the AEIM and MWCM models performed well and was feasible and applicable for regional winter wheat LAI inversion.

3.7. AEIM-MWCM cross validation

To further verify the accuracy of the inversion model, LAI images of the five typical growth stages of winter wheat obtained by inversion were resampled to 500 km, and 100 points were randomly selected. The LAI (GF-3 LAI) of the five growth stages of winter wheat obtained from the GF-3 remote sensing images and the MODIS LAI products (MOD15A2) of the corresponding phases were regressed and analyzed to obtain the correlation of the five growth stages of winter wheat, as shown in Figure 8.

Figure 8. Correlation analysis between GF-3 LAI and MODIS LAI. (a–e) are from the jointing, booting, heading, filling, and maturing stages, respectively.

As shown in Figure 8, the MODIS LAI and GF-3 LAI were strongly correlated in the five typical growth stages of winter wheat. Their correlation was the strongest at the jointing stage, with R² = 0.8658 and RMSE = 0.3910, and it was followed by those of the heading (R² = 0.8595 and RMSE = 0.3998), booting (R² = 0.8479 and RMSE = 0.4318), filling (R² = 0.8393 and RMSE = 0.4211), and maturing stages (R² = 0.8374 and RMSE = 0.4236). According to the statistics of the mean values of 100 sampling points, the mean values of GF-3 LAI and MODIS LAI at the jointing stage were 1.45 and 1.25, respectively. At the booting stage, the mean value of the GF-3 LAI was 2.68, and the mean value of the MODIS LAI was 2.54. The mean LAI values of the GF-3 and MODIS at the heading stage were 4.43 and 4.21, respectively. The mean value of the GF-3 LAI was 2.65, and the mean value of the MODIS LAI was 2.59 at the filling stage. The mean LAI values of the GF-3 and MODIS at the maturing stage were 1.31 and 1.14, respectively. The LAI of winter wheat at the five growth stages obtained by the MODIS was lower than the GF-3 LAI, which was related to the low spatial resolution of the MODIS and serious pixel mixing. The phenomenon of over- or under-estimation found in this study has been confirmed in other studies (Li et al. 2021).

3.8. Mapping of LAI inversion by AEIM-MWCM

Based on the verified reliability of each inversion model, the AEIM-MWCM was employed to generate winter wheat LAI spatial distribution maps for the five growth stages, as shown in Figure 9.

Figure 9. LAI images inverted by AEIM-MWCM. (a–e) are from the jointing, booting, heading, filling, and maturing stages, respectively.

The investigation results revealed that winter wheat in the test region was mainly distributed in the central areas. As shown in Figures 9 (a)–9(e), the non-vegetated area was red, and the vegetated area was colored. The LAI of winter wheat at the jointing stage was between 1 and 2, with an average value of 1.67. The larger LAI values of winter wheat were mainly concentrated in the central area, which was caused by the earlier sowing of winter wheat in this region. The LAI values of winter wheat at the booting, heading, and filling stages were 2.3–2.5, 2.9–4.3, and 2.5–3, with average values of 2.39, 3.94, and 2.52, respectively. Meanwhile, the larger LAI was mainly concentrated in the central area, whereas the LAI of winter wheat in the southwest area was generally smaller. The LAI value of winter wheat in the maturation stage was between 1 and 1.5, with an average value of 1.12. The LAI of winter wheat was mainly concentrated in the south. In summary, the LAI results of winter wheat in the five growth stages obtained by the AEIM-MWCM model were consistent with the local winter wheat growth status and situation.

However, there are several limitations in the application of the proposed model. For instance, in this study, the parameter calibration process is performed manually, so the obtained parameters might not be optimal. In the future, the automatic parameter optimization effect could be strengthened through a combination of different global optimization algorithms.

4. Conclusions

In this study, the AIEM-MWCM model is used to invert the LAI of winter wheat in the Xiangfu area east of Kaifeng City, Henan Province, during the entire growth cycle. The results indicate that the improved SAR-BM3D filtering method reduces the loss of detail information caused by the excessive smoothing of SAR images. The results also indicate that in the RVI_Freeman model, the problem of vegetation index saturation can be overcome to a certain extent by the penetration of radar. The results also show that the MWCM has the potential to provide improvements in the LAI estimation and discrimination of the scattering response from vegetation and bare soil. In addition, it is shown that the modeled backscattering coefficient calculated by the MWCM is better than that obtained by the WCM. Therefore, the MWCM can be successfully used to estimate the LAI of wheat crops. The AIEM-MWCM model has higher fitting accuracy than the Dubois-MWCM model in inverting the LAI of winter wheat during the entire growth cycle, with R² higher than 0.8 and RMSE lower than 0.22, which corresponds to the actual growth status of local winter wheat. The experiments in this paper are based on limited field measurement data in a single study area, resulting in some limitations in the model and conclusions. It is necessary to apply the model to different regions, expand the sample data, and conduct repeated verification experiments in subsequent studies.

Acknowledgement

We appreciate the detailed suggestions and comments of anonymous reviewers. The authors thank Guosheng Cai, Yushi Zhou, and Qinggang Zhang for assistance with flight tests and the other members of the Institute of Remote Sensing and Surveying and Mapping for assistance with LAI measurements.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The code used in this study are available by contacting the corresponding author.

Additional information

Funding

This research was funded by the 2016 National Key Research and Development Plan (grant number 2016YFC0803103), Research on Key Technology of Agricultural Remote Sensing Monitoring (grant number 12210243), and the Henan Provincial University Innovation Team Support Plan (grant number 14IRTSTHN026).