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

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 (Fe and De). 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 (R2 > 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.


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 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 (Fe and De), the water content of leaves, stalks, ears, and the total rice canopy (mv_l, mv_s, mv_e, and mv). The estimated results of the rice variables and the comparison between the MWCM and the WCM are discussed.

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 cover s a total area of 40 × 30 km 2 . 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 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 2 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 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 2 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.
where Fl, Fs, Fe, and Fp are the fresh biomass of leaves, stems, ears, and total rice canopy, respectively. Dl, Ds, De, and Dp are the dry biomass of leaves, stems, ears, and total rice canopy, respectively. Np is the planting density. hl, hs, he, and h denote the height of leaves, stem, ear layers, and total rice canopy, respectively.

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).

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). 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.

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., 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.

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).
where Fl, Fs, Fe, and Fp are the fresh biomass of leaves, stems, ears, and total rice canopy, respectively. Dl, Ds, De, and Dp are the dry biomass of leaves, stems, ears, and total rice canopy, respectively. Np is the planting density. hl, hs, he, and h denote the height of leaves, stem, ear layers, and total rice canopy, respectively.

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).   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.  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 (LL) and the width of the cuboid is equal to the row spacing (LR) 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 (DL and DS) 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 (DL), whose width is the short axis of the rice axis (DS), 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 − (DL × Ds)/(LL × LR), 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 W1 and that in the space part is W2, then W2 = n1 × W1, where coefficient n1 denotes the ratio of the VWC of the leaves in the space part (W2) to that in the rice part (W1), 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 n1 to be very small. From the tillering to the booting stage, rice leaves became larger, denser, and more inclined. Thus, the difference between W1 and W2 became smaller, resulting in a large value of n1. 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., W1-W2. Thus, n1 was close to one. From the milk to the mature stage, although the VWC in the rice part (W1) and the space part (W2) both decreased due 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 1 and that in the space part is W 2 , then W 2 = n 1 × W 1 , where coefficient n 1 denotes the ratio of the VWC of the leaves in the space part (W 2 ) to that in the rice part (W 1 ), 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 1 to be very small. From the tillering to the booting stage, rice leaves became larger, denser, and more inclined. Thus, the difference between W 1 and W 2 became smaller, resulting in a large value of n 1 . 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 1 -W 2 . Thus, n 1 was close to one. From the milk to the mature stage, although the VWC in the rice part (W 1 ) and the space part (W 2 ) both decreased due to the withered rice canopy, there was still no big difference between W 1 and W 2 because the leaves distributed randomly, implying that the coefficient n 1 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 3 and W 4 respectively, then W 4 = n 2 × W 3 , where coefficient n 2 denotes the ratio of the VWC of ears in the space part (W 4 ) to that in the rice part (W 3 ), 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 2 ≈ 0. From the milk to the mature stage, the ears became inclined due to the increase of biomass and, consequently, increased n 2 . Based on the analysis above, the value intervals of F, n 1 , and n 2 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.

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 2 ) and the rice part (W 1 ) 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.

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.

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 crosspolarized 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).

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., <SHHSHV*> and <SHV*SVV> [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:

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.

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).

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 HH S 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: 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).

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 | 2 /|S HH | 2 ) ≤ 2 dB; Yamaguchi et al. [35] added the horizontal and vertical dipole volume models when 10log (|S VV | 2 /|S HH | 2 ) ≤ −2 dB and 2 dB ≤ 10log (|S VV | 2 /|S HH | 2 ), respectively. In this paper, the generalized volume scattering model (GVSM) [40] is adopted: where γ = |S HH | 2 /|S VV | 2 , 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 Sections 3.2.1-3.2.3, the improved decomposition method is explained as follows: 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.
where f s , f d , f c , α, and β are as same as the parameters in Freeman and Yamaguchi decomposition.
where fs, fd, fc, α, and β are as same as the parameters in Freeman and Yamaguchi decomposition. γ=|SHH| 2 /|SVV| 2 , 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 Pd, Pv, and Ps from both the improved and the Freeman decomposition are similar. However, the 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).
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). It can be seen that the negative pixels with the improved decomposition are much fewer than those of the Freeman decomposition.

Scenario for Estimation of Rice Variables
Through applying Ps, Pd, and Pv from the improved polarimetric decomposition into the MWCM, three equations relating the three scattering components to the rice variables are built as follows: where the scattering components of V, D, and S with subscripts could be found in Figure 5.

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. 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). It can be seen that the negative pixels with the improved decomposition are much fewer than those of the Freeman decomposition.

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: where the scattering components of V, D, and S with subscripts could be found in Figure 5. (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].
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). 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) 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).

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 1 , and n 2 , 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.

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.

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.
Remote Sens. 2016 , 8, 878  14 of 21   2  2  2  ,  ,  ,  ,  ,  ,  , • 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.

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.

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 2 of the "y = x" line is calculated. In this study, the R 2 values for all rice variables, except mv_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 2 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

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.
with the mv, mv_l, and mv_e, the R 2 of mv_s is the minimum (0.7587) and the RMSE of mv_s is the largest (0.85 kg/m 3 ) (Figure 9g) because the estimated errors of mv_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 2 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 mv_s after the heading stage need more consideration due to the large scattering attenuation from leaves and ears.  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 2 of the "y = x" line is calculated. In this study, the R 2 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 2 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,  (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 2 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.

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.

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 AM and RM, respectively, and the absolute and relative errors with the WCM are abbreviated as AW and RW, respectively. Then Δa and Δr are attained by calculating (AW − AM) and (RW − RM), respectively ( Figure 10). It is noted that Δa of LAI, mv, and Me 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 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.

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. especially for LAI, h, mv, mv_l, and mv_s in the vegetation phase of the rice growth season and for Fe, De, and mv_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.

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.  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 1 , and n 2 are added. In this study, F, n 1 , and n 2 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 2 [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.

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 2 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 2 of m v_s is smaller (0.7585) and the RMSE is larger (0.85 kg/m 3 ). 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 1 , and n 2 in the MWCM are determined in general cases, but more applications of the MWCM are necessary in future work.