Estimation of Paddy Rice Variables with a Modified Water Cloud Model and Improved Polarimetric Decomposition Using Multi-Temporal RADARSAT-2 Images

Abstract

Rice growth monitoring is very important as rice is one of the staple crops of the world. Rice variables as quantitative indicators of rice growth are critical for farming management and yield estimation, and synthetic aperture radar (SAR) has great advantages for monitoring rice variables due to its all-weather observation capability. In this study, eight temporal RADARSAT-2 full-polarimetric SAR images were acquired during rice growth cycle and a modified water cloud model (MWCM) was proposed, in which the heterogeneity of the rice canopy in the horizontal direction and its phenological changes were considered when the double-bounce scattering between the rice canopy and the underlying surface was firstly considered as well. Then, three scattering components from an improved polarimetric decomposition were coupled with the MWCM, instead of the backscattering coefficients. Using a genetic algorithm, eight rice variables were estimated, such as the leaf area index (LAI), rice height (h), and the fresh and dry biomass of ears (F_e and D_e). The accuracy validation showed the MWCM was suitable for the estimation of rice variables during the whole growth season. The validation results showed that the MWCM could predict the temporal behaviors of the rice variables well during the growth cycle (R² > 0.8). Compared with the original water cloud model (WCM), the relative errors of rice variables with the MWCM were much smaller, especially in the vegetation phase (approximately 15% smaller). Finally, it was discussed that the MWCM could be used, theoretically, for extensive applications since the empirical coefficients in the MWCM were determined in general cases, but more applications of the MWCM are necessary in future work.

1. Introduction

Rice is the main staple food in Asian countries, which feeds about 3.5 billion people worldwide and accounts for 90% of the global rice supply [1]. The estimation of rice variables arouses people’s interest [2,3] because rice variables are quantitative indicators for rice growth and productivity, and provide reliable information for agricultural management applications. Considering rice always grows in the rainy season in almost all Asian countries, the all-weather capability of synthetic aperture radar (SAR) makes it suitable for monitoring rice paddy fields when it is also sensitive to the structure of the rice canopy [4,5,6,7,8,9].

In order to relate SAR signatures to the rice variables, numerous scattering models have been developed, including empirical [9,10], semi-empirical [11,12], and physical models [13,14]. Among the three methods, the empirical model is the simplest way to estimate the rice variables, however, it lacks a theoretical basis and is only suitable for the specific study site [15]. The physical model consists of complicated equations describing the interaction between radar signals and targets accurately, but it is too complex to be applied in practice and insufficient input parameters could cause ambiguities [3]. The semi-empirical model as good compromise between the empirical and physical model, considering both the physical background and the simplicity [15]. Therefore, the semi-empirical model is the most popular method for rice variable estimation [5,12,16].

The water cloud model (WCM) is one of the typical semi-empirical scattering models, which was firstly proposed by Attema and Ulaby [17] and was then applied to the assessment of soil moisture and vegetation variables extensively since it is simple and practical [18,19,20,21]. However, the WCM simplicity is often accompanied by relatively low accuracies because many assumptions and simplifications are contained in the model [22]. In order to improve the accuracy of WCM, some researchers improved the scattering mechanisms of the WCM. For example, Ulaby et al. [23] developed the multi-layered vegetation model and separately considered the backscatter from different parts of the vegetation (leaves, stems, and ears). Kweon and Oh [22] modified the WCM by adding the leaf angle distribution function to accurately depict the canopy with angular scattering particles. Moreover, many other studies focused on the calibration of the model coefficients of WCM for specific area, vegetation type or SAR data. For instance, Prevot et al. [24] adjusted the input parameters of the WCM with the field data and expressed the WCM briefly with four empirical coefficients. Van Leeuwen and Clevers [25] undertook a study attempting to simulate the LAI values for sugar beet using both optical and SAR data. The WCM was used for the SAR data and a canopy parameter of the WCM related to attenuation was modified. Xu et al. [26] recommended a modified WCM by calibrating the model coefficients using ERS-1 SAR data and ground measurements, which made the root mean squared error (RMSE) of the leaf area index (LAI) less than 0.4. In addition, the WCM was also improved through considering a different optimization algorithm for estimating the model coefficients. Generally, the model coefficients of the WCM are estimated by nonlinear least squares regression analysis by fitting the model against experimental datasets [24,26,27]. For example, Prevot et al. [24] fitted the model coefficients of the WCM by using a quasi-Newton minimization algorithm, a first-order method based on Newtonian mechanics. Xu et al. [26] fitted the model coefficients on in situ data using a combined Gauss-Newton and modified Newton minimizing algorithm. Bolten et al. [28] used a quasi-Newtonian method to estimate A and B. The reliability of parameter estimates by least squares optimization depends greatly upon the quality of the experimental data and the nature of the objective function. Errors in the experimental data are inevitable. During the experiment, measurement errors generally occur due to human mistakes, the malfunction of instruments or a lack of precision. Recently, Kumar et al. [12] estimated the coefficients of the WCM using a genetic algorithm optimization and found that the root mean squared error obtained from the GA was about half of that obtained by the conventional gradient-based sequential unconstrained minimization technique.

According to the analysis above, it is seen that many improvements on the WCM were considered, including the calibration of the model coefficients using the in situ data and the optimization algorithms for accurate estimation results. However, there are still some flaws that are ignored in the WCM. First, the distribution of crop leaves in the horizontal direction is not random, especially in the early period of the growth season, which conflicts with the assumptions of the WCM. Second, the backscattering coefficients were always used directly in the WCM to build the relationship with the vegetation variables. The correlation between the different polarimetric channels could cause the equations of the WCM to be mutually related to some extent, which could produce problems [29]. In addition, the double-bounce was not considered in the WCM because the second-order scatterings between the canopy particles and the soil surface were ignored. However, the double-bounce actually makes significant contributions, especially in the early period of the rice growth cycle [30]. Thus the polarimetric scattering mechanism assumed in the WCM is not particularly well adapted to the scattering from crops, especially for rice.

In order to solve the mentioned problems above, the present study aims to propose a modified water cloud model (MWCM), in which the heterogeneity of the rice canopy in the horizontal direction and the double-bounce scattering are added. The double-bounce, volume, and surface scattering components from an improved polarimetric decomposition are related to the rice variables with the MWCM, rather than the backscattering coefficients. Using a genetic algorithm (GA) to solve the MWCM, eight of the rice variables were estimated, including the leaf area index (LAI), rice height (h), the fresh and dry biomass of ears (F_e and D_e), the water content of leaves, stalks, ears, and the total rice canopy (m_{v_l}, m_{v_s}, m_{v_e}, and m_v). The estimated results of the rice variables and the comparison between the MWCM and the WCM are discussed.

2. Materials and Methods

2.1. Study Site

Our study was conducted in Jinhu, Jiangsu Province (33°07′05″N, 118°59′55.14″E), a major rice production region in Eastern China (Figure 1). The test site is located in a plain southeast of Hung-tse Lake and covers a total area of 40 × 30 km². It lies at a geographic elevation of approximately 20 m above sea level with a gentle and flat topography, where the annual average air temperature and total precipitation are 14.6 °C and 1085 mm, respectively.

In the test site, the transplanted hybrid rice fields were dominant and several varieties were included, such as “LIANGYOU-898” and “XIEYOU-9308”. All rice fields were cultivated once a year from late June to late October. The whole rice growth cycle was divided into three major development phases, which were subdivided into ten phenological stages according to the BBCH scale (from Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie) [31]: (1) Vegetative phase (from the seedling to the elongation stage): in the seedling stage, the bundles of 5–10 seedlings were transplanted with a spacing of 30 × 15 cm² under flooded fields. Then the seedlings split into several tillers and began to develop a root system in the tillering stage. The rice stems grew much longer in the elongation stage; (2) Reproductive phase (from the booting to the flowering stage): the upper part of stems thickened, and the flag leaves unfolded in the booting stage. The heading stage started with the emergence of the panicles. After the completion of the panicles, the anthers appeared and fully developed in the flowering stage; and (3) Maturity phase (from the dough to the mature stage): ears ripened and were ready to harvest (Figure 2). Normally, plant growth within each field was almost uniform when between-field variability was relatively large due to the difference in the soil conditions and the agronomic management practices.

2.2. Acquisition and Processing of RADARSAT-2 Images

In order to collect the rice growth conditions in different phenological stages, a stack of eight RADARSAT-2 full-polarimetric single-look complex (SLC) C-band SAR images were acquired from 27 June to 15 October 2012 (

Table 1. Available RADARSAT-2 images and corresponding rice phenology stages. Date: yy-mm-dd.

Date	Mode	Incidence Angle(°)	Pixel Spacing (A × R, m)	Phenology
2012-06-27	FQ20W 1	38–41	5.2 × 7.6	Seedling
2012-07-11	FQ9W	27–30	5.2 × 7.6	Tillering
2012-07-21	FQ20W	38–41	5.2 × 7.6	Elongation
2012-08-04	FQ9W	27–30	5.2 × 7.6	Booting
2012-08-28	FQ9W	27–30	5.2 × 7.6	Heading
2012-09-07	FQ20W	38–41	5.2 × 7.6	Flowering
2012-09-21	FQ9W	27–30	5.2 × 7.6	Dough
2012-10-15	FQ9W	27–30	5.2 × 7.6	Mature

¹ FQ**W stands for wide fine-quad mode.

). The eight SAR images were in two beam positions, FQ9W and FQ20W, which results in different incidence angles. The spatial resolution was higher than 10 m, which is suitable for estimating rice variables in field scale. The Next ESA SAR Toolbox (NEST) 5.1 (http://nest.array.ca/web/nest) was utilized for pre-processing these SAR images. The radiometric correction was applied to all image data based on radiometric parameters provided in the header files. Then a 7 × 7 enhanced Lee filter was applied to each image to reduce speckle noise, and all of the images were georectified based on the digital elevation map (DEM) using the Range-Doppler terrain model.

2.3. Ground Truth Measurements of Rice Variables

Detailed ground truth measurements were collected at 32 rice sample fields concurrently with the SAR data acquisitions. All sample fields in the ground campaign were 120 m × 120 m or larger, to ensure a sufficient number of pixels were included for analysis. For each rice field, the rice phenological stage was recorded based on the BBCH scale (Table 1), and the typical rice canopies were photographed (Figure 2). Then five or more representative rice plants were chosen randomly in each rice field, on which various rice variables were measured, such as the line spacing (L_L), row spacing (L_R), density of plants in the rice field (N_p), leaf angles, major and minor diameters of the rice cluster (D_L and D_S), etc. In detail, the LAI was measured using LAI-2200 and the height was measured by tapeline, whose precision is to the millimeter. The height between the lowest leaf node of the stem and the underlying surface was defined as the height of stems (h_s). The height between the highest leaf node of the stem and the vertex of the ears was defined as the height of the ear layer (h_e). The height between the lowest and highest leaf node was defined as the height of the leaf layer (h_l). The height of total rice canopy (h) was defined as the sum of h_l, h_s, and h_e. All of the rice variables were measured more than five times and then averaged to get the final results for each rice field. The wet and dry biomass of leaves, stems, ears and whole plants were obtained through destructive measurements of the sampled plants based on normal gravimetric processes. The volumetric water content (VWC) of leaves, stems, ears, and the total rice canopy (m_{v_l}, m_{v_s}, m_{v_e}, and m_v) were calculated by the fresh and dry biomass from different parts of the rice plants, respectively, i.e.,

$$\begin{array}{c}{m}_{v\_l}=\left(F_l-D_l\right)N_p/h_l\\ m_{v\_s}=\left(F_s-D_s\right)N_p/h_s\\ m_{v\_e}=\left(F_e-D_e\right)N_p/h_e\\ m_v=\left(F_p-D_p\right)N_p/h\end{array}\},$$

(1)

where F_l, F_s, F_e, and F_p are the fresh biomass of leaves, stems, ears, and total rice canopy, respectively. D_l, D_s, D_e, and D_p are the dry biomass of leaves, stems, ears, and total rice canopy, respectively. N_p is the planting density. h_l, h_s, h_e, and h denote the height of leaves, stem, ear layers, and total rice canopy, respectively.

3. Methodology

3.1. Modified Water Cloud Model (MWCM)

From Figure 2, it is found that the distribution of leaves and ears in the horizontal direction is not random, thus the water content of the rice canopy in the horizontal direction is not uniform, which conflicts with the assumptions of the WCM [17,18]. In order to solve this problem, a “scattering cell” is defined, including two parts, i.e., the rice and the space parts (Figure 3).

3.1.1. Scattering Cell

From Figure 4a, it is seen that the rice field consists of many scattering cells and Figure 4b shows that a scattering cell is regarded as a cuboid including the leaves, stalks and air. The length of the cuboid is equal to the line spacing (L_L) and the width of the cuboid is equal to the row spacing (L_R) when the height of the cuboid is equal to the height of the rice cluster (h). In Figure 4c, the lowest leaf nodes of all stems of the rice cluster are illustrated when the leaves become invisible. It is assumed that the nodes make up an ellipse and the long and short diameter axes (D_L and D_S) of the rice cluster are measured in the ground truth measurement, respectively. Then the rice cluster in each scattering cell could be expressed as a cuboid, whose length is the long axis of the rice cluster (D_L), whose width is the short axis of the rice axis (D_S), and whose height is equal to the rice height (Figure 4d). Thus, a volume fraction coefficient F is defined as the ratio of the space volume to the scattering cell, i.e., F = 1 − (D_L × D_s)/(L_L × L_R), 0 < F < 1. In the seedling stage, the diameter of rice cluster is small due to the small rice plants. The space between the adjacent rice clusters accounts for a large portion of the scattering cell, implying F is large. In the tillering and elongation stage, many tillers are generated and the stems develop much longer. Thus, the long and short diameter axes of the rice cluster become large and the proportion of the rice cluster increases, which decreases F. From the heading to the mature stage, no more tillers are generated and the length of the stems remains stable; however, the inclination angles of stems and leaves become slightly larger. Thus, the diameter axes of rice cluster become larger, which causes F to slightly decrease.

In addition, it is noted that more of the surrounding air and less leaves are included in the space, compared with the rice cluster. Considering the VWC of air is much lower than that of the rice plant, the VWC of the rice part is larger than that of the space part. Assuming the VWC of the leaves in the rice part is W₁ and that in the space part is W₂, then W₂ = n₁ × W₁, where coefficient n₁ denotes the ratio of the VWC of the leaves in the space part (W₂) to that in the rice part (W₁), ranging from 0 to 1. In the seedling stage, rice plants and leaves were small. Almost all leaves were included in the rice parts and few leaves ran out to the space part, causing n₁ to be very small. From the tillering to the booting stage, rice leaves became larger, denser, and more inclined. Thus, the difference between W₁ and W₂ became smaller, resulting in a large value of n₁. From the heading to the flowering stage, the underlying surface of the rice field was almost completely covered by leaves. The density of leaves in the scattering cell was nearly homogenous, i.e., W₁–W₂. Thus, n₁ was close to one. From the milk to the mature stage, although the VWC in the rice part (W₁) and the space part (W₂) both decreased due to the withered rice canopy, there was still no big difference between W₁ and W₂ because the leaves distributed randomly, implying that the coefficient n₁ was nearly one as well.

Similar with the analysis above, assuming the VWC of ears in the rice part and that in the space are W₃ and W₄ respectively, then W₄ = n₂ × W₃, where coefficient n₂ denotes the ratio of the VWC of ears in the space part (W₄) to that in the rice part (W₃), ranging from 0 to 1. From the heading to the flowering stage, the ears were almost straight with small inclination angles (up to 5° approximately) because of their small biomass. Thus it was assumed that the ears were only distributed in the rice part during this period, i.e., n₂ ≈ 0. From the milk to the mature stage, the ears became inclined due to the increase of biomass and, consequently, increased n₂. Based on the analysis above, the value intervals of F, n₁, and n₂ were assumed empirically according to the ground measurement, which would be applied to the genetic algorithm (GA) as the initial conditions in different phenological stages.

3.1.2. Scattering Components from Multi-Layered Rice Canopy

Based on the scattering cell, the multi-layered rice canopy proposed by Ulaby et al. [23] is considered and Figure 5 shows the variation of the multi-layers at different phenological stages. In the seedling stage, the rice plants are short and mainly consisted of small leaves. The rice cluster of the scattering cell is considered as a layer consisting of identical water particles that are uniformly distributed. In addition, as almost no leaves were distributed in the space of the scattering cell, the water content in the space is negligible in the seedling stage. Therefore, four scattering mechanisms are considered in the seedling stage, including the volume scattering from the rice cluster (V_{f_r}), the underlying surface scattering through the rice cluster (S_{g_r}), the underlying surface scattering through the space (S_{g_s}), and the double-bounce between underlying surface and rice cluster (D_{g_f}) (Figure 5a). From the tillering to the booting stage, the volumetric water content in the space increases as the rice plant and leaves develop. The stem layer is also considered in the rice cluster since the rice stems grow longer. Thus, three more scattering interactions are considered in the MWCM, including the volume scattering from the leaf layer in the space (V_{f_s}), surface scattering from the stem layer (S_t), and the double-bounce between the underlying surface and the stem (D_{g_t}) (Figure 5b). After the booting stage, the ears emerge and the ear layer is considered in the MWCM (Figure 5c,d). From the heading to the flowering stage, rice plants were fully developed and completely covered the underlying surface, which made the VWC of the space (W₂) and the rice part (W₁) being almost equal. Consequently, the double-bounce scattering between the underlying surface and the leaf layer (D_{g_f}) was ignored. Meanwhile, rice ears came out with small inclination angles and were almost embodied in the rice part. Accordingly, the volumetric scattering from the ear layer in the rice part (V_{e_r}) and the double-bounce between the underlying surface and the ear layer (D_{g_e}) were added (Figure 5c). From the milk to the mature stage, rice ears leaned over and ran out to the space part due to the increasing biomass. Accordingly, the volumetric scattering from the ear layer in the space part (V_{e_s}) was considered (Figure 5d).

It is noted that the underlying surface of the rice field changes along with the phenological stages. During the field experiments, it was found that the underlying surfaces of all rice fields were covered by water in the seedling stage. From the tillering to the booting stage, the underlying surfaces of most rice fields were flooded, except for some that were covered by mud. From the heading to the mature stage, nearly all of the underlying surfaces of rice fields were covered by wet soil. For convenience, it is assumed that rice fields are flooded from the seedling to the booting stage and then change to wet soil from the heading stage to the mature stage.

3.2. Improved Polarimetric Decomposition

The model-based decomposition has been extensively applied for its simplicity and easy implementation [32,33,34]. However, this decomposition methods has some deficiencies. Due to uncertainty and unsuitability of the volume scattering model, it usually leads to the overestimation of volume scattering contribution [34,35]. Additionally, for the flaw of the traditional estimation method for scattering model parameters, such as the assumption of reflection symmetry, there exists the problem of negative powers [36]. In order to solve the mentioned problems above and to enhance the accuracy of decomposition, a new decomposition method is proposed, including three improvements.

3.2.1. Deorientation

First, the deorientation [37,38] is applied. The purpose of the deorientation process is to remove the fluctuant influence of randomly-distributed target orientation angles on polarimetric scattering, and it can make identical targets with different orientation angles yield the same polarimetric matrix, which will finally lead to the same decomposition results.

In the deorientation process, the coherency matrix is rotated by a certain angle about the line of sight so that the cross-polarized scattering power is minimized. The minimization of the cross-polarized scattering means that the power is concentrated to the co-polarized channels. More detail information could be found in [37]. The deorientation process is denoted as C_θ = DEORIENT (C).

3.2.2. Distinguish the Reflection Symmetry

The difference between the reflection symmetry and non-reflection symmetry targets can be described by the product of the cross-polarization and co-polarization component in the Sinclair matrix, i.e., <S_HHS_HV*> and <S_HV*S_VV> [39]. In terms of the reflection symmetry target, the product is nearly zero. Thus, the averaged cross-polarization coefficient (ρ) can be generated to distinguish the reflection and non-reflection symmetry target:

$$\rho =\frac{1}{2}\left|\frac{\langle S_{HH}{S}_{HV}^{*}\rangle }{\sqrt{\langle {\left|S_{HH}\right|}^2\rangle \langle {\left|S_{HV}\right|}^2\rangle }}+\frac{\langle S_{HV}{S}_{VV}^{*}\rangle }{\sqrt{\langle {\left|S_{HV}\right|}^2\rangle \langle {\left|S_{VV}\right|}^2\rangle }}\right|,$$

(2)

where S_HH, S_HV, and S_VV are the parameters of the Sinclair matrix. The subscript * denotes the conjugate.

The range of ρ is [0,1]. When the value of ρ is close to 0, the target is inclined to reflection symmetry; and the target is inclined to non-reflection symmetry when ρ is larger than an appointed threshold (such as 0.1).

3.2.3. Generalized Volume Scattering Model

Modifying the volume scattering model is a good method to reduce the overestimation of the volume scattering power. Freeman [32] considered the forest canopy with randomly-oriented dipoles when −2 dB ≤ 10log (|S_VV|²/|S_HH|²) ≤ 2 dB; Yamaguchi et al. [35] added the horizontal and vertical dipole volume models when 10log (|S_VV|²/|S_HH|²) ≤ −2 dB and 2 dB ≤ 10log (|S_VV|²/|S_HH|²), respectively. In this paper, the generalized volume scattering model (GVSM) [40] is adopted:

$$C_{volume}=\frac{1}{\frac{3}{2}(1+\gamma )-\frac{\sqrt{\gamma }}{3}}\left[\begin{array}{ccc}\gamma & 0& \frac{\sqrt{\gamma }}{3}\\ 0& \frac{1+\gamma }{2}-\frac{\sqrt{\gamma }}{3}& 0\\ \frac{\sqrt{\gamma }}{3}& 0& 1\end{array}\right],$$

(3)

where γ = |S_HH|²/|S_VV|², denotes the HH and VV ratio.

From Equation (3), the GVSM can be used to describe the volume scattering continuously with the variation of γ [40]. When γ = 1, the GVSM is equivalent to the volume model in the Freeman decomposition; and when γ = 8/3 and 3/8, the GVSM is equivalent to the vertical and horizontal dipole volume models in the Yamaguchi decomposition.

Based on Section 3.2.1, Section 3.2.2 and Section 3.2.3, the improved decomposition method is explained as follows:

$$\begin{array}{c}\langle \left[C\right]\rangle =f_s\langle \left[C_{surface}\right]\rangle +f_d\langle \left[C_{double}\right]\rangle +f_v\langle \left[C_{volume}\right]\rangle \\ +\{\begin{array}{cc}\hfill f_c\langle \left[C_{helix}\right]\rangle ,& if\rho \ge \mathrm{threshold}\left(\mathrm{such}\mathrm{as}0.1\right)\hfill \\ \hfill 0,& if\rho \le \mathrm{threshold}\left(\mathrm{such}\mathrm{as}0.1\right)\hfill \end{array}\end{array},$$

(4)

where f_s, f_d, f_v, and f_c are coefficients to be determined. <[C_surface]>, <[C_double]>, <[C_volume]>, and <[C_helix]> are matrices corresponding to the surface, double-bounce, volume, and helix scattering mechanisms, respectively. Except <[C_volume]>, the other scattering mechanisms and matrices are adopted as same as that of the Yamaguchi decomposition.

Expanding the covariance matrices in Equation (4), Equation (5) is attained, from which the coefficients f_s, f_d, f_v, and f_c are calculated by using the similar method in the previous studies [34,35]. The flow chart of the improved decomposition method is shown in Figure 6.

$$\begin{array}{c}\left[\begin{array}{ccc}{C}_{11}& C_{12}& C_{13}\\ C_{21}& C_{22}& C_{23}\\ C_{31}& C_{32}& C_{33}\end{array}\right]=f_s\left[\begin{array}{ccc}{\left|\beta \right|}^2& 0& \beta \\ 0& 0& 0\\ {\beta }^{*}& 0& 1\end{array}\right]+f_d\left[\begin{array}{ccc}{\left|\alpha \right|}^2& 0& \alpha \\ 0& 0& 0\\ {\alpha }^{*}& 0& 1\end{array}\right]+\frac{f_v}{\frac{3}{2}(1+\gamma )-\frac{\sqrt{\gamma }}{3}}\left[\begin{array}{ccc}\gamma & 0& \frac{\sqrt{\gamma }}{3}\\ 0& \frac{1+\gamma }{2}-\frac{\sqrt{\gamma }}{3}& 0\\ \frac{\sqrt{\gamma }}{3}& 0& 1\end{array}\right]\\ +\{\begin{array}{ccc}\frac{f_c}{4}\left[\begin{array}{ccc}1& \mp j\sqrt{2}& -1\\ \pm j\sqrt{2}& 2& \mp j\sqrt{2}\\ {\beta }^{*}& \pm j\sqrt{2}& 1\end{array}\right]& ,& if\rho \ge \mathrm{threshold}(\mathrm{such}\mathrm{as}0.1)\\ 0& ,& if\rho \le \mathrm{threshold}(\mathrm{such}\mathrm{as}0.1)\end{array}\end{array},$$

(5)

where f_s, f_d, f_c, α, and β are as same as the parameters in Freeman and Yamaguchi decomposition. γ=|S_HH|²/|S_VV|², denotes the HH and VV ratio.

The results of the polarimetric scattering components could be attained during the whole growth season (Figure 7a). It is noted that the decomposition components have an unexpected decrease on 21 July and 7 September because the incident angles of SAR images on these dates are much larger than those on other dates (as seen in Table 1).

Compared with the results of the Freeman decomposition (Figure 7), the temporal behaviors of the P_d, P_v, and P_s from both the improved and the Freeman decomposition are similar. However, the P_v is significantly overestimated and the two other scattering components are underestimated with the Freeman decomposition. For example, the double-bounce scattering is dominant with the improved decomposition in the seedling stage, instead of the volume scattering. The possible reason is that the rice canopy is very small, leading to small volume scattering when the interaction between the specular reflection from water and the rice plants make the double-bounce dominant. In addition, the pixels denoting the negative scattering power are also considered. We subset some rice fields and then counted the negative pixels (

Table 2. Numbers of negative pixels in six rice fields.

No. of Rice Fields	Total Number of Pixels	Number of Negative Pixels		Percentage of Negative Pixels (%)
		Freeman Decomposition	Improved Decomposition	Freeman Decomposition	Improved Decomposition
(1)	408	31	10	7.60	2.45
(2)	2238	106	65	4.69	2.90
(3)	1516	240	80	15.83	5.28
(4)	518	71	32	13.90	6.18
(5)	213	19	8	9.39	3.76
(6)	556	22	13	4.14	2.34

). It can be seen that the negative pixels with the improved decomposition are much fewer than those of the Freeman decomposition.

3.3. Scenario for Estimation of Rice Variables

Through applying P_s, P_d, and P_v from the improved polarimetric decomposition into the MWCM, three equations relating the three scattering components to the rice variables are built as follows:

$$\begin{array}{l}{P}_v=V_{e\_r}+V_{e\_s}+V_{f\_r}+V_{f\_s}\\ P_d=D_{g-e}+D_{g-t}+D_{g-f}\\ P_s=S_{g\_r}+S_{g\_s}+S_t\end{array}\},$$

(6)

where the scattering components of V, D, and S with subscripts could be found in Figure 5.

3.3.1. Different Expressions of Equation (6)

Based on the previous studies [23,24,25,26,27], the scattering components of V, D, and S with subscripts in the right part of Equation (6) could be expressed as the function of different rice variables.

● LAI, h, m_{v_s}, and D_e

The scattering components in the right part of Equation (6) can be expressed as the function of LAI, h, m_{v_s}, and D_e (Equations (7)–(21)), based on the theory proposed by Ulaby et al. [23].

$$V_{e\_r}=(1-F)\times A_{e1}(\theta )\times D_e,$$

(7)

$$V_{e\_s}=F\times A_{e1}(\theta )\times n_2{D}_e,$$

(8)

$$V_{f\_r}=(1-F)\times A_{f1}\times (1-\exp (-B_{f1}L/h))\cos (\theta )\times (1-{\gamma }_{f\_r}^2(\theta )){\gamma }_{e\_r}^2(\theta ),$$

(9)

$$V_{f\_s}=F\times A_{f1}\times (1-\exp (-B_{f1}{n}_{1}L/h))\cos (\theta )\times (1-{\gamma }_{f\_s}^2(\theta )){\gamma }_{e\_s}^2(\theta ),$$

(10)

$$S_t=(1-F)\times A_{t1}{m}_{v}h\times {\gamma }_{f\_r}^2(\theta )\times {\gamma }_{e\_r}^2(\theta ),$$

(11)

$$S_{g\_r}=(1-F)\times C_{g1}(\theta )m_s\times {\gamma }_{f\_r}^2(\theta ){\gamma }_{e\_r}^2(\theta ){\gamma }_t^2(\theta ),$$

(12)

$$S_{g\_s}=F\times C_{g1}(\theta )m_s\times {\gamma }_{f\_s}^2(\theta ){\gamma }_{e\_t}^2(\theta ),$$

(13)

$$D_{g-f}=F\times C_{g2}(\theta )m_s\times A_{f2}\times (1-\exp (-B_{f2}L/h))\times {\gamma }_{f\_s}^2(\theta ),$$

(14)

$$D_{g-e}=F\times C_{g2}(\theta )m_s\times A_{e2}(\theta )\times D_e{\gamma }_{f\_s}^2(\theta ){\gamma }_{e\_s}^2(\theta ),$$

(15)

$$D_{g-t}=F\times C_{g2}(\theta )m_s\times A_{t2}{m}_{v}h\times {\gamma }_{f\_s}^2(\theta ){\gamma }_{e\_s}^2(\theta ),$$

(16)

where θ is the incident angle of the SAR image and:

$${\gamma }_{f\_r}^2(\theta )=\exp (-2{\alpha }_{f}L\sec (\theta )),$$

(17)

$${\gamma }_{f\_s}^2(\theta )=\exp (-2{\alpha }_f{n}_{1}L\sec (\theta )),$$

(18)

$${\gamma }_t^2(\theta )=\exp (-2{\alpha }_t{m}_{v}h\sec (\theta )),$$

(19)

$${\gamma }_{e\_r}^2(\theta )=\exp (-2{\alpha }_e{D}_e\sec (\theta )),$$

(20)

$${\gamma }_{e\_s}^2(\theta )=\exp (-2{\alpha }_e{n}_2{D}_e\sec (\theta )),$$

(21)

From Equations (7)–(21), there are twenty unknowns, including sixteen model coefficients and five rice variables. The sixteen model coefficients include F, n₁, n₂, A_e1(θ), A_e2(θ), A_f1, B_f1, A_f2, B_f2, A_t1, A_t2, C_g1(θ), C_g2(θ), α_f, α_t, and α_e. The four rice variables to estimate are D_e, L, h, and m_{v_s}.

Before the heading stage, there is no ear layer and the underlying surface of the rice field is water. Thus, D_e is not considered and the backscattering coefficient of the underlying surface related to m_s is replaced by the Fresnel reflection coefficient of water. Therefore, three rice variables, LAI, h, and m_{v_s}, are easily estimated.

From the heading to the mature stage, the ear is booting and the underlying surface of the rice field is wet soil. Thus, five rice variables need to be considered, including D_e, LAI, h, m_{v_s}, and m_s. Two assumptions are adopted here. First, we assume the backscattering coefficient of water is constant, i.e., C_g × m_s = C^’_g, which can eliminate the variable m_s. Second, the stem layer is ignored after the heading stage considering that the attenuation from ear and leaf layers are large and the scattering contribution from the stem layer is small, which results in the variable m_{v_s} being eliminated. According to the two assumptions, the three remaining rice variables D_e, LAI, and h are easily estimated by Equation (6).

● m_{v_l}, h, m_{v_s}, and D_e

The scattering components in the right part of Equation (6) could also be expressed as the function of m_{v_l}, h, m_{v_s}, and D_e. Considering the quantity m_{v_l}h for the top (leaf) layer is related to the wet and dry biomasses of the leaves, which, in turn, are related to the green LAI [23], we can replace the LAI in Equations (8)–(17) by the m_{v_l}h. Then the m_{v_l}, h, m_{v_s}, and D_e could be estimated similarly, i.e., the variables m_{v_l}, h and m_{v_s} are estimated before the heading stage and the variables D_e, m_{v_l}, and h are estimated after the heading stage.

● LAI, h, m_{v_s} and F_e

Similarly, assuming that D_e is linearly correlated to F_e, i.e., F_e = k × D_e, F_e could be estimated instead of D_e through replacing D_e in Equations (7), (14), and (19) by F_e.

● LAI, h, m_{v_s}, and m_{v_e}

Similar with the consideration that the quantity m_{v_l}h for the top (leaf) layer is related to the wet and dry biomasses of the leaves, which, in turn, are related to the green LAI [23], the m_{v_e}h could be estimated when the D_e in Equations (7), (14), and (19) by m_{v_e}h.

Except the four situations above, the scattering components could be expressed as different cases, such as “m_v, h, m_{v_s}, and D_e” and “LAI, h, and m_{v_s}”. Thus, we can estimate different rice variables through expressing the scattering components from the rice canopy as different combinations of rice variables in the right part of Equation (6).

3.3.2. Genetic Algorithm (GA)

A genetic algorithm [41] is used to solve Equation (6) for estimating rice variables in eight phenological stages. Nine rice fields, distributed uniformly in the study area, are used as training data for the model training and the remaining 23 rice fields are used for validation. Usually, a simple GA for optimization comprises six major components, namely genetic representation, reproduction, cross-over, mutation, a fitness function, and termination criteria.

● Initialization and genetic representation

In the genetic representation of the GA, each coefficient above was encoded by using binary coding as the gene substring. For F, n₁, and n₂, their initial values and ranges were determined by the ground data measurements and empirical assumption (Section 3.1.1). In terms of the other coefficients, their initial values and ranges were initialized by solving the MWCM with a set of training data. For example, the attenuation coefficients α_f and α_e were initialized to [0,1] and [0,4], respectively.

Based on their ranges, the gene substring lengths of the coefficients were determined when the desired solution accuracy to four decimal places was considered. Then the length of the chromosome (total of the gene substring) was determined, i.e., 161 in the seeding stage (27 June), 175 from the tillering (11 July) to the booting stage (4 August), 135 from the heading (28 August) to the flowering stage (7 September), and 149 from the dough (21 September) to the mature stage (15 October). The population of the coefficients was generated using a random generator and the population size of 100 was initialized.

● Reproduction

The chromosomes generated in the initial population were then chosen for participation in the reproduction process based on their fitness values [16]. In this study, the fitness value was calculated as Equations (22) and (23) and proportionate fitness selection was used. A chromosome with a higher fitness value had a higher probability of being copied into the cross-over pool.

$$f=\frac{1}{1+SSE},$$

(22)

$$SSE={\displaystyle \sum _{i=1}^n\left[{(P_{v,obs,i}-P_{v,est,i})}^2+{(P_{d,obs,i}-P_{d,est,i})}^2+{(P_{s,obs,i}-P_{s,est,i})}^2\right]}$$

(23)

● Cross-over

Cross-over is a recombinant operator that selects two chromosomes from the cross-over pool at random and cuts them into bits at a randomly chosen position [41]. The number of strings participating in mating depended on the cross-over probability. In this study, the cross-over probability was assumed to 0.8 and one-point cross-over was selected.

● Mutation

Mutation is an important process that permits new genetic material to be introduced to a population. A mutation probability is specified that permits random mutations to be made to individual genes (e.g., changing 1 to 0, and vice versa, for binary GAs). The mutation probability of 0.02 was selected in this study.

● Termination criteria

Finally, the termination criterion of the GA process was determined. The GA process could be stopped when the fitness criterion was satisfied, or the maximum number of generations was exceeded. In this study, the maximum number of generations was 5000.

3.3.3. Flow Chart for Rice Variable Estimation

Based on the MWCM, the improved decomposition components, and the genetic algorithm, the flow chart for rice variable estimation is shown (Figure 8), by which the coefficients of the MWCM and the rice variables in different phenological stages are estimated in Section 4.

4. Results and Discussion

4.1. Estimated Results

Based on the scenario for the estimation of rice variables (Section 3.3), nine rice fields, distributed uniformly in the study area, are used as training data for the GA and the remaining twenty-three rice fields are used for validation. The estimated results of rice variables in different phenological stages are shown in Figure 9.

For each variable, most estimated results from different fields (denoted as different dots in Figure 9) are around the ground truth data. Several statistical indicators are used to validate the performance of the MWCM in the multi-temporal prediction of rice variables. First, the determination coefficient R² of the “y = x” line is calculated. In this study, the R² values for all rice variables, except m_{v_s}, are larger than 0.8, implying that the “y = x” line can express the relationship between the estimated and ground truth data well. The R² of h is the largest (0.8914) (Figure 9d), denoting that the difference between the estimated results and ground truth of h is the smallest. The possible reason is that h is related to all scattering and attenuation from all layers in the rice plants (leaves, stalks, and ears). Then we attain that the RMSE of LAI, h, F_e, D_e, m_v, m_{v_l}, m_{v_s}, and m_{v_e} are 0.48, 10.37 cm, 0.27 kg/m², 0.22 kg/m², 0.57 kg/m³, 0.35 kg/m³, 0.85 kg/m³, and 0.46 kg/m³, respectively. Compared with the m_v, m_{v_l}, and m_{v_e}, the R² of m_{v_s} is the minimum (0.7587) and the RMSE of m_{v_s} is the largest (0.85 kg/m³) (Figure 9g) because the estimated errors of m_{v_s} are large, especially from the heading to the mature stage, considering the scattering attenuation from leaves and ears is very large. Therefore, the large R² and the small RMSE of all rice variables suggest the good performance of the MWCM in the multi-temporal prediction of rice variables when the estimated results of m_{v_s} after the heading stage need more consideration due to the large scattering attenuation from leaves and ears.

4.2. Comparison with the WCM

In order to prove the capability of the MWCM for improving the accuracy of the rice variables, the estimated results with the MWCM should be compared to the estimated results with the WCM proposed by Ulaby et al. [23].

For each rice variable, the absolute and relative errors with the MWCM are abbreviated as A_M and R_M, respectively, and the absolute and relative errors with the WCM are abbreviated as A_W and R_W, respectively. Then Δ_a and Δ_r are attained by calculating (A_W − A_M) and (R_W − R_M), respectively (Figure 10). It is noted that Δ_a of LAI, m_v, and M_e are multiplied by 100 for display in Figure 10, but the coefficient of 100 is ignored when Δ_a is analyzed in the paper.

From Figure 10, it is seen that Δ_a and Δ_r of all rice variables in all phenological stages are larger than zero, implying the absolute and relative errors with the WCM are larger than those with the MWCM in all phenological stages. For the LAI (Figure 10a), Δ_r is about 15% in the seedling stage (27 June), which is the largest in all phenological stages. With rice growth, Δ_r decrease. It is found that Δ_r from the seedling to the booting stage are clearly larger than those from the heading (28 August) to the mature stage (15 October). From the seedling to the booting stage (4 August), Δ_r is larger than 8%; from the heading to the mature stage (15 October), Δ_r is about 5%. Therefore, the MWCM improves the estimated accuracy of LAI efficiently, especially in the vegetative phase of the rice growth season (from the seedling to the booting stage). The possible reason is that the WCM assumes that the water particles are distributed uniformly in the rice canopy, which is not consistent with the fact that the water particles in the rice canopy are actually distributed heterogeneously in the vegetative phase of the rice growth season (Figure 5a,b). The MWCM greatly improves the estimated accuracy by considering the heterogeneity of water particles in this period. For h, m_v, m_{v_l}, and m_{v_s}, the temporal behavior of Δ_r is very similar with those of the LAI (Figure 10b,e–g). In addition, the largest Δ_r of F_e, D_e, and m_{v_e} are about 13% in the heading stage, and then decreases with rice growth (Figure 10c,d,h. The possible reason is the WCM ignores that the water particles are distributed heterogeneously in the heading stage when the ears just form and grow straight. The MWCM analyzes the variation of water content in the rice canopy during the whole rice growth season (Figure 5), which makes the estimated accuracy of rice variables greatly improved. In addition, the polarimetric scattering components, rather than backscattering coefficients, are used in the MWCM, which reduces the correlation of the model and is beneficial for the estimation of the rice variables. From all of the above, it can be found that the MWCM improves the estimated accuracy of rice variables efficiently, especially for LAI, h, m_v, m_{v_l}, and m_{v_s} in the vegetation phase of the rice growth season and for F_e, D_e, and m_{v_e} in the heading stage. Δ_a could be analyzed similarly considering that Δ_a is related to the absolute value of each rice variable in each phenological stage.

4.3. Discussion for the Key Improvements of the MWCM

Upon the analyses above, two main improvements in the MWCM result in the superiority of the MWCM (Figure 10), i.e., the consideration of the rice canopy heterogeneity in the horizontal direction and the application of the polarimetric decomposition components instead of the backscattering coefficients. In this section, the effects caused by the two improvements were quantified, respectively, and compared with each other.

In Figure 11, three different cases are set. In CASE I, the heterogeneity of the rice canopy in the horizontal direction is considered, when the backscattering coefficients are used in the MWCM. In CASE II, the scattering components from the improved polarimetric decomposition are used in the MWCM, while the heterogeneity of rice canopy in the horizontal direction is ignored. In CASE III, both the heterogeneity of the rice canopy and the decomposition components are considered in the MWCM. Then the relative error of LAI, h, and m_v are obtained for the three cases, which are compared with those with the original WCM, similar with Section 4.2 (Figure 11). It is noted that the LAI, h, and m_v are shown here as an example, and the similar results could be obtained for the other variables.

During the whole rice growth cycle, Δ_r for CASE III was the largest, implying that the estimation accuracy is the best by considering both the heterogeneity of rice canopy and the decomposition components. From the seedling to the booting stage (from 27 June to 4 August), Δ_r for CASE I (larger than 10%) is much larger than that for CASE II (2%–3%), implying that the estimation accuracy is improved considerably by considering the heterogeneity of rice canopy during the this period, whereas the estimation accuracy is just improved approximately 3% by applying the decomposition components. From the heading to the mature stage (from 28 August to 15 October), Δ_r for CASE II (approximately 5%) is larger than that for CASE I (less than 2%), implying the application of decomposition components in the MWCM could improve the estimation accuracy more than the consideration of the rice canopy heterogeneity during the reproductive phase of rice growth cycle. In sum, the heterogeneity of the rice canopy is a very essential consideration in the MWCM, especially during the vegetative phase. As the rice canopy becomes dense and uniform, the heterogeneity of rice canopy decreases. The application of decomposition components in the MWCM could improve the estimation accuracy with 3%–5% during the rice growth cycle.

Furthermore, it is important to consider whether or not the MWCM can be applied in general cases. The MWCM proposed in the study is mainly aimed at transplanting rice fields. The distribution of the rice plants in the direct-sown rice fields is more random than that in the transplanting rice fields, and the distribution of scattering elements (leaves, stalks, and ears) in the horizontal direction in the direct-sown rice fields is more uniform than that in the transplanting fields. Thus, the heterogeneity of water content of the rice canopy in the horizontal direction in the transplanting fields is more obvious than that of the direct-sown rice fields, which is considered in the MWCM (Section 3.1). In addition, the model coefficients of the MWCM are similar with those of the WCM except that the empirical parameters F, n₁, and n₂ are added. In this study, F, n₁, and n₂ in the MWCM are determined statistically based on the ground measurement (Section 2), which are mainly related to the transplanting condition and the rice phenological stage. The transplanting method most commonly used in China and Japan is considered, i.e., a bundle of 5–10 seedlings of 15–30 cm long are transplanted at a spacing of 30 × 15 cm² [3,42]. Eight phenological stages are taken into account, i.e., the seedling, tillering, elongation, booting, heading, flowering, dough, and mature stages. Accordingly, the MWCM with the statistical intervals of F and n in this study could be used in general cases, theoretically, but more applications of the MWCM are necessary in future work.

5. Conclusions

In this study, the heterogeneity of the rice canopy in the horizontal direction was considered and the double-bounce scattering was added in the MWCM. Then, improved polarimetric decomposition was adopted, which reduced the overestimation of volume scattering and the negative pixels from the Freeman decomposition. Afterwards, the scattering components P_d, P_v, and P_s from the improved decomposition were applied into the MWCM, instead of the polarimetric backscattering coefficients. By using a GA, eight rice variables were estimated in this study, including LAI, h, F_e, D_e, m_{v_l}, m_{v_s}, m_{v_e}, and m_v. The accuracy validation showed that the R² of all rice variables, except m_{v_s}, were larger than 0.8, and their RMSE were relatively small during the whole rice growth season, suggesting the good performance of the MWCM in the multi-temporal prediction of rice variables. Compared with m_{v_l}, m_{v_e}, and m_v, the R² of m_{v_s} is smaller (0.7585) and the RMSE is larger (0.85 kg/m³). The possible reason was that the large attenuation from leaves and ears decreased the scattering from the stalk layer, resulting in the large errors of m_{v_s}. The large scattering attenuation from ears impacted the scattering from the leaves layer, adding the errors of m_{v_l}. The scattering from the ear layer was not attenuated by any rice layers; thus, the relative errors of m_{v_e} were the smallest.

In addition, compared with the WCM, the relative errors of rice variables with the MWCM were much smaller than those with the WCM, especially for LAI, h, m_v, m_{v_l}, and m_{v_s} in the vegetation phase of the rice growth season (approximately 15% smaller) and for F_e, D_e, and m_{v_e} in the heading stage (approximately 13% smaller). Thus, the performance of the MWCM for the estimation of rice variables was validated and it was significant to consider the heterogeneity of the rice canopy in the horizontal direction for estimation of rice variables.

Finally, it was considered that the MWCM could be used for extensive applications, at least theoretically, since the empirical parameters F, n₁, and n₂ in the MWCM are determined in general cases, but more applications of the MWCM are necessary in future work.