Rice Phenology Retrieval Based on Growth Curve Simulation and Multi-Temporal Sentinel-1 Data

Abstract

The accurate estimation and monitoring of phenology is necessary for modern agricultural industries. For crops with short phenology occurrence times, such as rice, Sentinel-1 can be used to effectively monitor the growth status in different phenology periods within a short time interval. Therefore, this study proposes a method to monitor rice phenology based on growth curve simulation by constructing a polarized growth index (PGI) and obtaining a polarized growth curve. A recursive neural network is used to realize the classification of phenology and use it as prior knowledge of rice phenology to divide and extract the phenological interval and date of rice in 2021. The experimental results show that the average accuracy of neural network phenological interval division reaches 93.5%, and the average error between the extracted and measured phenological date is 3.08 days, which proves the application potential of the method. This study will contribute to the technical development of planning, management and maintenance of renewable energy infrastructure related to phenology.

1. Introduction

Global agriculture is currently facing problems such as land and water scarcity, increasing pollution from surface sources, declining ecosystem stability and threats to food security, which seriously limit the achievement of the UN Sustainable Development goals. Therefore, modern agriculture needs technology and new energy sources as the driving force to advance agriculture. As a new type of land use, photovoltaic agriculture is a product of the development of a sophisticated combination of modern agriculture and clean energy, which can give full play to the advantages of zero pollution in photovoltaic power generation [1]. It enables the conversion of solar photovoltaic power to agricultural production activities: supporting greenhouse irrigation systems, supplementing light to crops, meeting the heating needs of greenhouse crops, etc. However, the combination of photovoltaic and agriculture still lacks theoretical research [2]. In order to arrange the photovoltaic cells without affecting the light of crops, the laying out of photovoltaic cells must take into account the demand for light intensity of crops in agricultural greenhouses. In addition, different crops have different temperature requirements. It is necessary to combine local land resources and environmental characteristics to design different sizes of photovoltaic power plant agricultural greenhouses required for different crops.

Crop phenology refers to the rhythm of growth and development of crops in response to periodic changes in light, precipitation, temperature and other conditions [3], which can reflect the characteristics of crop light and temperature requirements. Therefore, a study of the phenological characteristics is useful as a guide for the development of the photovoltaic agriculture industry. Rice is a crop that humans rely on for survival. The real-time grasp of rice phenology has important theoretical value and practical significance in greenhouse management, planning and decision-making [4,5,6]. In addition, rice phenology is one of the necessary inputs in rice growth and yield prediction models, ecosystem productivity models and surface process models [7,8,9] to facilitate the timely deployment of field management activities [10] such as fertilization, irrigation, weeding and pest control and to provide technical support for large-scale precision agriculture. As part of the terrestrial ecosystem, crop phenological study provides a foundation for a variety of research directions, including monitoring the crop interannual variability of surface features and long-term trends, according to the season for the phase spectrum response to classify land cover types, analysis of the dynamic phenology and climate variable correlation function between the measures of the response to climate change, evaluating the terrestrial carbon balance during the year and the interannual fluctuation [11].

Crop monitoring necessitates a high level of timeliness due to the quick growth of field crops and the severe dynamic changes that occur over time. Traditional field crop monitoring provides accurate biophysical parameter measurements and yield estimates. However, it is time-consuming, labor-intensive and expensive [12,13]. Satellite remote sensing is used as an alternative, providing global coverage data for large-scale crop phenology studies. The time-series remote sensing vegetation index is widely used in crop phenology detection [6,14,15,16]. High-resolution optical remote sensing is not always available in agricultural applications owing to interference from haze, clouds and rain. Satellites equipped with synthetic aperture radar (SAR) can operate 24/7 and are not affected by day/night. Thus, they play an important role in rice recognition and monitoring in cloudy and foggy areas [17,18]. Lopez-Sanchez et al. [19] explored the potential of TerraSAR-X data in rice phenology recognition in a threshold classification approach. Vicente-Guijalba et al. [20] used a Kalman filter to detect the phenology of barley, wheat and oats based on C-band Radarsat-2. Küçük et al. [21] retrieved rice phenology based on TerraSAR-X data of the same pole over multiple periods of time with machine learning and explored the potential of machine learning algorithms in phenological phase recognition. Some studies [21,22,23] accurately detected the correct phenological stage at a certain point in time. Owing to the relatively low time resolution of the SAR sensor, it is worth noting that none of them utilized time-series metrics to detect the dates of the phenological stages, as is common in optical remote sensing.

Finer temporal resolution data can capture more accurate vegetation phenological states. However, high temporal and spatial resolution cannot be taken into account at the same data source. Thus, how to balance the demand for high spatial and temporal resolution is critical. The launch of Sentinel-1 A and B solved the above problems. Sentinel-1, with C-band high-frequency, dual-polarization channels, spatial resolution and six-day short revisit time, has the same acquisition characteristics due to the constellation of two satellites, Sentinel 1A (S1A) and 1B (S1B), and is a potential SAR dataset for rice crop monitoring. Several studies [20,24,25,26] have demonstrated that C-band backscatter is particularly sensitive to the phenological cycle of crops based on their sensitivity to crop parameters. The short revisit time of Sentinel-1 allows the detection of instant morphological variations in crop growth, as observable phonological changes take a short time. This increased availability of open SAR data has facilitated the research of dense time-series at a spatial resolution of 10–20 m. Previous research focused on full-polarization data, and there were few studies on the application of the dual-polarization mode. Compared to full-band acquisition, the dual-band modes have the advantage of larger sweep width and less data volume but at the cost of limited polarization measurement information, offering some benefits for agencies in ongoing operational activities [27,28,29,30,31,32]. Consequently, Sentinel-1 images with fine spatial and temporal scales are used for crop monitoring and tracking the dynamics of crop growth [21,28,32,33,34,35,36].

Most methods of extracting phenological dates are based on the mathematical theory of feature point detection [22,23,37,38]. Raw time-series remote sensing vegetation index datasets often have noise errors, and smoothing is used to minimize the noise and reconstruct more representative indexes. However, these methods are often more sensitive to local fluctuations and data noise. In contrast, methods based on the entire annual time-series or low-order harmonic series usually yield smoother curves but may deviate from the actual growth trajectory and suffer from overcorrection. It is difficult to determine the best method in some cases. These research methods have no physiological basis, which limits the accuracy of the detection of specific phenolic indicators. The time-series analysis of the SAR data analysis is still in the exploratory stage.

The objective of this study was to construct a phenological change curve for rice using Sentinel-1 data in order to estimate the date of each phenological stage. This will provide a basis for the planning and design of photovoltaic agricultural projects that adapt to changes in climatic phenomena and environmental conditions. The current study [8,18,19,21] mainly applied classification methods to PolSAR features to classify the phenological intervals. These studies do not consider the temporal evolution of the crop and are unable to detect the complete date of the phenological stage. Moreover, the current research is limited to the use of data intensity features and band combination features for SAR data only, with few effective characterization indices that can describe the actual growth trajectory of rice.

For the above reasons, this study uses Sentinel-1 in a rice phenology monitoring method based on growth curve simulation. Combining rice growth morphology and key phenological periods, this study divided the Sentinel-1 time-series data to establish a simulated growth curve of rice, so that there is a good distinction between the phenological stages, and use it as prior knowledge to divide and extract the phenological interval and data, as follows: (1) This study analyzed and sorted the collected data in the experimental area and established the relationship between the phenological scale and the day of the year. Combined with the key phenological period of rice and the transit time of satellite photography, the Sentinel-1 time-series data in 2018 were divided into intervals of rice phenology. (2) Based on the 2018 Sentinel-1 time-series data, with the rice growth curve as the target, the objective programming function is used with the weighted fusion of the backscattering coefficients and H/α polarization decomposition parameters to construct a polarized growth index (PGI), so that the fusion feature changes can approach the growth pattern of crops. (3) Recurrent neural networks are used to train various features and their growth curve simulation features and are compared with the phenological phase classification results. In this study, the features with the best classification effects were selected as prior knowledge to screen the 2021 Sentinel-1 time-series data, divide the corresponding phenological interval, verify the reliability of the phenological division through the neural network and calculate the error between the extracted and measured phenological dates. This study can promote the combination of the photovoltaic industry and agricultural industry and contribute to the realization of agricultural science and technology and agricultural industrialization.

2. Study Area and Datasets

2.1. Study Area

The rice test areas were located in Jinhu (33°1′ N, 118°52′ E) and Liyang (31°30′ N, 119°12′ E), Jiangsu Province, China (Figure 1). This is in the Yangtze River Delta, the main rice-producing area in Eastern China, and has a subtropical monsoon climate, abundant rainfall, sufficient sunshine and a long frost-free period. The climatic and environmental conditions in the experimental area can meet the basic growing needs of the crops, and the long hours of sunshine, the flat topography of the farms and the fertile soil provide favorable conditions for the growth and development of local crops.

2.2. Ground Data

To ensure that the obtained Sentinel-1 data can accurately reflect the growth status of rice in the study area, it is necessary to conduct a field survey of rice at the time of satellite transit, as follows:

• Arrange several sampling sites evenly in the rice planting area and use GPS to record their locations.
• The field survey time interval is 3 to 4 days, and the time needs to cover the satellite transit time to guarantee that ground data collecting and radar satellite imaging are synchronized.
• According to the date of rice planting and the current status of the data acquisition channels, conduct field investigations to collect rice samples and obtain the growth parameters. Collect the ground data of rice samples at the sampling points in the rice-growing area. The contents of the collection include rice species, rice plant height, growth period, planting methods, collection time and weather conditions.

In this study, 80 sampling sites were set up in 2018 and 30 sampling sites in 2021 to determine the plant height of the rice crop during the growing season. Sampling was conducted at eight stages of the rice crop after transplanting: tillering, elongation, booting, heading, flowering, milk, dough and ripening. Plant height refers to the distance from the base of the plant to the uppermost heart leaves, measured with a straightedge. Based on the fieldwork results, the two regions are mainly planted in indica rice fields. The growth cycle is generally from mid-June to mid-October, about 120 days. The phenological period of the crop is influenced by meteorological conditions and regional differences (Figure 2). It can be expected that the variations of the backscattering coefficient were different over the analyzed years.

The phenological period of rice is recorded according to the BBCH (Biologische Bundesanstalt, Bundessortenamt und CHemische Industrie) standard, an internationally accepted standard for crop phenology [39], using consecutive numbers from 0 to 99 to characterize the entire rice growth period, as shown in

Table 1. Each growth stage of rice corresponds to the BBCH scale.

Major Stage	BBCH Scale	Description
Vegetative	00–09	Germination
	10–19	Leaf development
	20–29	Tillering
	30–39	Stem elongation
Reproductive	40–49	Booting
	50–59	Inflorescence emergence
	60–69	Flowering, anthesis
Ripening	70–79	Development of fruit
	80–89	Ripening
	90–99	Senescence

. In the actual field measurement process, because growth environments differ by area, the growth of rice in a given field is completely inconsistent. Therefore, when more than half of the rice in the same paddy field enters a certain phenological stage, that stage is recorded for the entire paddy field.

In practical agricultural applications, the monitoring of changes in crop phenological phases only needs to be specific to a few main phenological phases, which can not only ensure the overall grasp of the crop growth status but also reduce information redundancy and detection costs. The growth morphology of rice is also the standard for dividing the phenological period of this experiment. The phenological period is divided according to the growth status of rice stems, leaves and ears, which is beneficial to the modeling of rice parameter inversion. The phenological growth stages and rice morphology recorded during the ground survey are shown in Figure 3. The phenological evolution of rice during the field survey is shown in Figure 4.

2.3. Sentinel-1 Data

The radar data adopted the C-band Sentinel-1 dual-polarization (vertical–horizontal polarization, VH; vertical–vertical polarization, VV) SLC product, the acquisition mode was interferometric wide swath (IW) and the imaging mode was terrain observation by progressive scans (TOPS). According to the rice growth cycle in the experimental area and the satellite revisit cycle, the data collection time was from June to October 2018 and June to October 2021. The Sentinel-1 mission provides data with a repeat cycle of 12 days using one satellite and 6 days using the two Sentinel-1 satellites. In 2018, only Sentinel-1A transited the Jinhu rice experimental area, and the satellite revisit period was 12 days. In 2021, Sentinel-1A and Sentinel-1B both passed through the Liyang rice experimental area, and the satellite revisit period was 6 days. The Sentinel-1 data parameters are shown in

Table 2. Sentinel-1 data parameters [20,36].

Satellite	Sentinel-1
Swath width	250 km
Revisit period	12 d or 6 d
Spatial resolution	5 m × 20 m
Relative orbit number	69
Acquisition date	15 June 2018 to 13 October 2018; 11 June 2021 to 15 October 2021
Incident angles	32–34°; 38–40°
Polarization scheme	VH, VV

.

The ENVI SARscape platform was used to preprocess the Sentinel-1 radar image. The ENVI SARscape platform is based on the ENVI remote sensing image processing software of the Swiss SARMAP company, with the addition of a module ENVI SARscape that specializes in processing SAR radar images. The key operations are as follows:

• Multiview processing: average the azimuth or distance of the SLC data to improve the data intensity;
• Speckle filtering: remove the inherent speckle noise of radar images;
• Geocoding and radiometric calibration: combined with the SRTM data of the image coverage area, complete the geocoding and radiometric calibration and extract the backscatter coefficient of the test area in the remote sensing image.

Based on the field phenological observation results, the experiment selected seven satellite transit time points to divide the phenological phase into six stages: seedling to early-tillering (10–24), mid- to late-tillering (25–29), elongation to early-booting (30–44), mid-booting to early-flowering (45–64), mid-flowering to milk (65–79) and maturity (80–99). The time-series of backscattering coefficients of rice in the Jinhu area in 2018 is shown in Figure 5.

Unlike the full polarization, each resolution cell at each time point of dual-polarization SAR is defined by a 2 × 2 covariance matrix (C₂):

$$C_2=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]=\left[\begin{array}{cc}\langle S_{VV}{S}_{VV}^{*}\rangle & \langle S_{VV}{S}_{VH}^{*}\rangle \\ \langle S_{VH}{S}_{VV}^{*}\rangle & \langle S_{VH}{S}_{VH}^{*}\rangle \end{array}\right]$$

(1)

where C₁₁ is the VV backscattering coefficient, C₂₂ is the VH backscattering coefficient and C₂₁ and C₂₂ are complex numbers. * Indicates complex conjugation.

C₂ is formed by the second-order statistics of the scattering matrix, which provide information about the features of a deterministic target. In this study, the eigenvector and eigenvalue-based H/α decomposition using the Sentinel Application Platform (SNAP), an open-source software of the European Space Agency (ESA), and the related basic operations are consistent with ENVI SARscape platform preprocessing, and the key operations are as follows:

• Calibration: For polarized SAR processing, its radiometric calibration is a complex calibration, which is performed separately for the real and imaginary parts of the complex numbers.
• Polarimetric Metrices: Since the Sentinel-1 satellite has, at most, two polarization channels, only the covariance matrix (C₂) can be generated.
• Polarimetric Decomposition: For dual-polarization decomposition, only “H-Alpha Dual Pol Decomposition” can be performed.

This approach used the eigenvalues and the eigenvector of the C₂ to calculate the parameters entropy (H) and alpha angle (α). The H and α are defined as

$$\alpha ={\displaystyle \sum }_{i=1}^2{P}_{i}co{s}^{-1}\left(u_i\right) with P_i=\frac{{\lambda }_i}{{\lambda }_1+{\lambda }_2}$$

(2)

$$H={\displaystyle \sum }_{i=1}^2 - P_{i}lo{g}_2{P}_i$$

(3)

where ${\lambda }_i$ ($i$ = 1, 2) is the eigenvalue, and $u_i$ ($i$ = 1, 2) is the corresponding eigenvector.

The value range of α is [0°, 90°]; values close to 0° indicate surface scattering, values close to 45° refer to volume scattering and values close to 90° indicate dihedral scattering. The value range of H is [0, 1], and its change of value from low to high describes the randomness of the target scattering process. That is, a low value of H indicates that only one scattering mechanism is dominant, and a higher value indicates more than two main scattering mechanisms. Multiple scattering mechanisms of a target should be considered comprehensively.

3. Methodology

A growth curve simulation is the construction of a curve of the polarized growth index (PGI) that approximates the growth law of rice to distinguish its various phenological stages, using the phenological classification results in 2018 as prior knowledge to retrieve the key phenological stages of rice in 2021. Its methodology has four parts: (i) rice plant height growth curve fitting, (ii) Gaussian fitting to satellite data, (iii) data fusion based on the rice canopy analysis and (iv) the polarized growth index obtained by goal programming.

3.1. Rice Growth Curve

As a growth curve, the logistic curve is often used in simulation studies of plant growth processes, including dry matter growth and leaf area dynamics. Crop growth is determined by crop, climate and soil factors. Plant height changes and dry matter accumulation during the growth period of crops will also show changes in compliance with the logistic curve [40], and the phenological information will have a corresponding stage. According to the available field data, the growth pattern of crops has a strong correlation with their height. The growth process of rice presents a change law that conforms to the logistic curve, i.e., growth is slower in the early stage, suddenly accelerates in the middle stage and becomes slower in the later stage. To objectively describe the growth process of rice, the plant height measured in the field at different growth periods was used as a benchmark. Plant height data are filtered and normalized, and the logistic curve model is used to fit the plant height data to obtain the growth curve of the target crop, as shown in Figure 6.

3.2. Gaussian Fitting

Since the units and ranges of each feature are completely different, it is necessary to separately count the normalized histograms of $S_{VH}$, $S_{VV}$, $H$ and $\alpha$ of rice to normalize the data. The black dots in Figure 7 show the normalized histograms of $S_{VH}$, $S_{VV}$, $H$, and $\alpha$ of rice on 15 June 2018. The abscissa represents the size of the normalized data, and the ordinate is the amount of data. Outliers can be found, which will lead to inaccurate analysis of the growth law of crops when directly using $S_{VH}$, $S_{VV}$, $H$ and $\alpha$. Therefore, to further analyze various features, it is necessary to use Gaussian fitting functions to fit the data of $S_{VH}$, $S_{VV}$, $H$ and $\alpha$ and select points of interest in each category.

Gaussian fitting can express parameters with clear physical meaning, such as peak shape, peak height and peak position in the data distribution. On the basis of keeping the data information valid, the original information is unified into a small number of uniquely determined Gaussian eigenvalue parameters to achieve the extraction and simplification of the original information. Figure 7 shows the fitting results of rice. The peak points of the fitted normalized histogram curve are the interest values $S_{VH}$, $S_{VV}$, $H$ and $\alpha$, and these will be used for the subsequent feature analysis of rice.

3.3. Data Fusion

During its entire growth period, as rice matures, the backscattering mechanism tends to be more complicated, and the changing law of different polarizations with growth provides a basis for the monitoring of growth. However, the backscatter characteristic changes with polarization tend to be more complicated, and the trend of the radar backscatter coefficient does not simply increase or decrease. The corresponding scattering parameters of H/α decomposition change regularly as the crop grows [41]. Specifically, the scattering mechanisms in the early stage of crop growth (e.g., seedling) are relatively simple and consist mainly of surface scattering, with scattering parameters corresponding to low-entropy scattering. As the morphology changes after emerging (e.g., stem elongation), the scattering parameters correspond to medium-entropy scattering. As time goes on, the planting density and crop height increase (e.g., flowering fruiting and maturing), which leads to a gradual increase in the stochasticity and scattering angle. Then, the scattering parameters represent high-entropy multiple scattering or high-entropy vegetation scattering, such as the scattering mechanisms of forest canopies and vegetation with random highly anisotropic scattering elements.

The mechanism model relies on the physical action process of the vegetation canopy structure, such as leaves and stems on electromagnetic waves, to simulate backscattering. The water cloud model proposed by Attema et al. [42] simplified the vegetation canopy as a uniform dielectric body composed of water droplets of equal size. When only single scattering is considered, it is considered that the backscattering coefficient of the vegetation canopy consists of two parts. One is canopy scattering from the body, and the other is soil scattering attenuated by the canopy, where the former is the main component. Then, the backscattering coefficient can be simulated as

$${\sigma }^0=A\ast cos\theta \ast V_1\ast \left[1-exp\left(-2B\ast V_2/cos\theta \right)\right]+{\sigma }_{soil}^0\ast exp\left(-2B\ast V_2/cos\theta \right)$$

(4)

where ${\sigma }_{soil}^0$ is the backscattering coefficients of soil, $V_1$ and $V_2$ can be described as the vegetation canopy structure, $\theta$ is the angle of incidence, ${\tau }^2$ is the two-way attenuation coefficient and $A$ and $B$ are the model coefficients.

$V_1$ is usually assumed to be constant, and $V_2$ describes the average plant height of rice. Finding the reverse solution to the equation, the average plant height of rice is

$$h=\frac{cos\theta }{2B}ln\left(A\ast cos\theta \ast V_1-{\sigma }_{soil}^0\right)-\frac{cos\theta }{2B}ln\left(A\ast cos\theta \ast V_1-{\sigma }^0\right)$$

(5)

Let $a=ln\left(A\ast cos\theta \ast V_1-{\sigma }_{soil}^0\right)$, $b=\frac{cos\theta }{2B}$ and $c=A\ast cos\theta \ast V_1$. Then, the relationship between rice plant height and the backscatter coefficient can be described as

$$h=a\ast b-b\ast ln\left(c-{\sigma }^0\right)$$

(6)

The backscattering coefficient and polarization decomposition parameters of rice are mainly affected by the rice canopy. Therefore, this research model focuses on the rice canopy analysis. Referring to the above water cloud model, the relationship between the rice plant height and backscattering coefficient is logarithmic, and $T_{VV}+T_{VH}$ characterizes the rice canopy changes. According to the characteristics that H and α gradually increase [43], $T_H\times T_{\alpha }$ is used to characterize the change of differential characteristics of the H/α distribution to realize data fusion:

$$\mu =T_H\times T_{\alpha }+ln\left(T_{VV}+T_{VH}\right)$$

(7)

which can transform into

$$\eta =e^{T_H\times T_{\alpha }}\times \left(T_{VV}+T_{VH}\right)$$

(8)

where $T_{VH}$, $T_{VV}$, $T_H$ and $T_{\alpha }$ are the interest values of the VV backscattering coefficient, VH backscattering coefficient, polarization scattering entropy $H$ and average scattering angle $\alpha$, respectively, and $\mu$ and $\eta$ are the fusion parameters.

The normalized backscattering and polarization decomposition parameters are fused, and the obtained preliminary fusion parameter timing diagram curve is shown in Figure 8. Since the same crops will tend to have the same scattering mechanism in the same period, the introduction of polarization parameters can constrain the boundary range of various phenological periods and, to a certain extent, eliminate the misclassification phenomenon caused by only using backscattering coefficients for phenological division.

Comparing the curve in Figure 8 with the growth curve in Figure 8, it can be found that the fusion parameter timing curve of rice is very different from the ideal growth curve. The error of the fusion timing curve is large, and it is difficult to characterize the growth status of the entire growth cycle of rice, so the differences between the values of adjacent time points are not obvious. Therefore, this study planned the values of each period of the model with the goal of approximating the growth curve.

3.4. Polarized Growth Index

During its entire growth period, as the rice matures, the backscattering mechanism tends to be more complicated, and the changing law of different polarizations with growth provides a basis for the monitoring of rice growth.

S_VH, S_VV, H and α play different roles in the characterization of rice growth characteristics. The backscattering coefficient is very sensitive to the geometric structure and water content of rice, and H and α are related to the changes in rice scattering. Therefore, to assign different weights to each feature during fusion can bring the polarization time-series curve of rice closer to the ideal growth curve. In this study, the goal programming method is used to calculate the weights of S_VH, S_VV, H, α and ρ represents the polarized growth index:

$$\rho =W_1\times e^{T_H\times T_{\alpha }}\times T_{VV}+W_2\times e^{T_H\times T_{\alpha }}\times T_{VH}+C$$

(9)

where W₁ and W₂ ensure that the error between the polarized growth curve and the ideal growth curve is minimized, and C is the regulation parameter for fitting the polarized growth curve.

The goal constraint P is

$$\rho =W_1\times e^{T_H\times T_{\alpha }}\times T_{VV}+W_2\times e^{T_H\times T_{\alpha }}\times T_{VH}+C+d^{+}+d^{-}$$

(10)

where $d^{+}$ and $d^{-}$ are the positive and negative deviations, respectively, of the goal programming constraints, i.e., the part on the right of Equation (9) of the decision value, which exceeds or does not reach the target value $\tau$. P is a function with weights A, B and C as variables.

The goal function is

$$Z=Min\left(\Sigma P\left(\left|d^{+}\right|+\left|d^{-}\right|\right)\right)$$

(11)

When both the positive and negative variances are as small as possible, $\left|d^{+}\right|+\left|d^{-}\right|$ ensures that the target constraint is established. At this time, the goal function Z is the smallest, and the polarization growth characteristic is consistent with those in the ideal growth curve. The best values of the optimal weights in different periods are obtained through goal programming. Substituting the weights in Equation (9), the characteristic values of the polarized growth in different periods are obtained. Figure 9 shows the polarization growth curve of rice generated by the PGI. By comparing the polarization growth curve (Figure 9) and the ideal growth curve of rice (Figure 6), the law of the polarization growth curve and the ideal growth curve at each stage are seen to be basically consistent.

4. Results

4.1. 2018 Phenology Classification

This study used the Long-Short Term Memory (LSTM) unit model of the recurrent neural network (RNN) for sample training. LSTM can effectively retain the contextual semantic temporal relationship when processing long text and long sequence data [44]. The experiment selected two farmlands with the same growing conditions in the 2018 Jinhu area as training and test data with the number of pixels as 2527 and 317, respectively, to verify the application potential of the proposed growth curve simulation method in dividing different rice phenological intervals. The experiment took the product of H and α as the single characteristic parameter of H/α and combined it with backscattering coefficients S_VH and S_VV, the fusion parameter (FP) and various parameters simulated by the growth curve to train the samples of adjacent periods in the Sentinel-1 data and obtain the corresponding classifiers to realize the division of the phenological intervals.

The training and classification results are shown in Figure 10 and

Table 3. Accuracy results of rice phenology recognition.

Class	0–24	25–29	30–44	45–64	65–79	80–99	Total
LSTMVH	1.57%	5.36%	26.81%	37.54%	60.88%	5.68%	22.97%
LSTMVV	72.56%	0.63%	2.52%	42.90%	46.69%	58.99%	37.38%
LSTMH/α	86.43%	80.44%	87.07%	15.14%	24.29%	81.07%	62.41%
LSTMFP	58.68%	13.25%	0.63%	47.00%	53.00%	30.60%	33.86%
LSTMVH_GC	100%	97.16%	77.92%	74.13%	13.56%	64.67%	71.24%
LSTMVV_GC	100%	89.59%	0%	34.70%	82.01%	92.43%	66.46%
LSTMH/α_GC	100%	83.60%	79.50%	62.46%	41.96%	94.32%	76.97%
LSTMPGI	96.21%	100%	99.05%	100%	93.38%	93.06%	96.95%

, respectively, from which it can be seen that the classification accuracy of the LSTM neural network using the VV backscatter coefficient is slightly higher than that of the VH backscatter coefficient. This is because VV polarization is more sensitive than VH polarization to the growth of rice. However, the overall classification accuracy using only the backscatter coefficient as the input feature is relatively low. The overall recognition accuracy of the LSTM neural network classification using the H/α single feature is higher than that of the backscatter coefficient single feature method, reaching 62.41%. This is because, in the process of rice growth, the H/α polarization parameter changes from low-entropy to high-entropy scattering and has a stable changing trend. The classification accuracy of the fusion parameter is only 33.86%, which shows that the equal-weight fusion of scattering and polarization features is not conducive to phenology classification. The classification accuracy of all types of parameters simulated by growth curves is improved compared to the other parameters. Among them, the overall classification accuracy of the PGI can reach 96%, and the accuracy of each category is generally higher than that of the other input features, which shows the superiority of the proposed polarization growth characteristic in the division of phenological intervals.

4.2. 2021 Phenology Retrieval

Based on the above experiments, the PGI with the highest classification accuracy was selected as the prior knowledge of rice phenology, and the phenology retrieval experiment of rice in the Liyang area in 2021 was conducted. The growth cycle of rice in the Liyang and Jinhu areas is from mid-June to mid-October, but there are differences in the specific planting and harvesting times and, also, the satellite transit time. Sentinel-1A and Sentinel-1B both transited the Liyang area in 2021, the satellite revisit period was 6 days and the data were more abundant. The time-series diagram of the backscatter coefficient of rice in the Liyang area in 2021 is shown in Figure 11.

Since the goal planning method adopted in this article weighted and integrated the rice data at seven time points in 2018, the PGI data of rice have different numerical ranges on the polarization growth curve for each period. In this study, the data for different time periods in 2021 were substituted into the model for calculation, and when the data for a certain period do not match the model parameters, the PGI is not in the correct range of values or even has a large error. Therefore, the data for the time point corresponding to 2018 can be obtained based on the idea of a stepwise regression analysis and minimum distance clustering. The screening results are shown in Figure 12, where dots in different colors represent the results of coupling different model parameters with the data in 2021, and black dots represent the data points of the polarization growth curve in 2018.

This study sorted the data screening results as a test set and input them in the trained LSTM classifier, with the classification results as shown in Figure 13. It can be seen from the experimental results that the proposed growth curve simulation method has good stability for crop phenology retrieval in different regions and times and has a good classification ability for the phenological intervals. The average accuracy of various types reaches 93.5%, which verifies the correctness of the data screening experiment results. The classification accuracy for mid-flowering to milk (65–79) is only 73%, because during data screening, the data points near the dough stage in 2021 differ too much from the target planning values.

According to the above classification results, rice phenology obtained from satellite data can be determined according to the BBCH criteria, and the rice phenology evolution curve can be obtained by relying on the satellite transit time. Figure 14 compares the field survey phenology evolution curves, from which it can be seen that the evolution curve of rice phenology obtained through satellite data has the same trend as that of field survey phenology. The phenological evolution curve is smoothed, the difference between the key phenology periods is calculated, and the accuracy of the key phenology extraction results is shown in

Table 4. Accuracy of the key phenology extraction results in 2021.

BBCH	Field Observation Time/d	Retrieval Estimate Time/d	Error/d
10	159	162	3
20	176	178	2
25	189	186	3
30	206	204	2
40	221	224	3
45	228	234	6
50	236	238.5	2.5
60	246	247.5	1.5
65	249	252	3
70	253	256	3
80	261	264	3
90	270	276	6
100	286	288	2
Total			3.08

, from which it can be seen that the maximum error between the key phenology DOY obtained by this method and the actual phenology DOY is six days, and the average error is 3.08 days.

5. Discussion

The radar signal strength (backscatter) is affected by the dielectric characteristics of the target, geometric structure and radar system parameters. Using the experimental data of Sentinel-1 in different phases in 2018 and 2021, the changes in the backscattering characteristics of rice and its polarization characteristics during the complete rice growth cycle can be obtained. The phenological time-series of the VV backscattering coefficient and VH scattering coefficient of rice in 2018 and 2021 are shown in Figure 15, and the phenological time-series of the H/α single characteristic parameter is shown in Figure 16.

As can be seen in Figure 15, the change trend of backscattering in rice first rises, then falls and, finally, stabilizes. The backscattering law shown in these time-series is related to the growth law of rice and its growth environment. Since VV polarization is easily affected by the vertical structure, it is sensitive to changes in the rice morphology. In a short period of time after rice transplanting, seedlings are sparsely exposed to the water surface, causing a large amount of interactions between the water surface and the rice canopy. During seedling to early-tillering (10–24), as the rice grows, the canopy becomes denser, the surface scattering decreases and the volume scattering increases, because the rice has an obvious vertical structure in the early stage of growth, and the backscattering coefficient increases. In mid-tillering to early-booting (25–44), the rice canopy gradually increases, the leaves gradually change from a vertical to a horizontal structure and VV polarization is greatly disturbed by the vegetation, which makes the attenuation of the backscattering coefficient obvious. From mid-booting to maturity (45–99), owing to the appearance of rice ears, the structure of rice is very different from before, the scattering changes caused by the ear density dominate and the VV polarization shows an increasing trend. VH polarization is mainly affected by multiple scattering between different targets. Since rice is sparse in the early growth stage, its backscattering coefficient increasing. When the leaf density increases to a certain extent, a uniform and dense surface is formed on the surface of the rice area, and the backscattering coefficient remains relatively stable, but it is still affected by changes in the geometric structure of the rice.

Jiao Guo et al. [43] verified that the scattering parameters of dual H/α decomposition vary regularly as the crop grows before harvest (H and α gradually increase). Therefore, the H/α single feature approximately increases with the growth of the rice. It can be seen from Figure 16 that the time-series features of the H/α single-feature parameter show an increasing variation law. In seedling to maturity (10–99), the biological morphology of rice changes, the plant height continues to increase, the canopy gradually increases, the scattering mechanism is no longer single and the H/α single characteristic parameter steadily rises. Tillering (20–29) is dominated by vegetative organs such as growing roots, leaves and tillers, and the number of emerging leaves accounts for about three-fourths of the total number of leaves, so the H/α single characteristic parameter changes greatly at this stage. In elongation to milk (30–79), rice changes from vegetative growth (such as of roots, stems and leaves) to reproductive growth (flowering and fruiting). The emergence of rice ears complicates the scattering mechanism, and the time-series features of the H/α single-feature parameter show an increasing variation law. In maturity (80–99), the ears of rice already protrude from the sheath of the flag leaf. After the rice pollen is scattered, the grains begin to fill, firm and become mature. The H/α single characteristic parameter changes more slowly than in the previous stage.

To further verify the superiority of the polarization growth characteristics, the standard deviations of various scattering characteristics (VH backscattering coefficient, VV backscattering coefficient, H/α single feature and polarization growth characteristic) of rice were compared. The standard deviation can reflect the degree of dispersion of a dataset. The larger the standard deviation, the more obvious the discrimination. The results are shown in Figure 17. During the growth cycle of the crop, the VH backscattering coefficient, VV backscattering coefficient, H/α single feature and polarization growth characteristic all have certain changes in the standard deviation. On the whole, the parameter value of H/α is greater than that of the two types of backscattering coefficients, and the standard deviation of the polarization growth characteristic is the largest. This is because, in the process of crop growth, H/α has a steady growth trend in the time-series. Since the distribution variation of the polarization growth characteristic of rice is similar to the overall change trend of the rice growth cycle, the standard deviation of the polarization growth parameters is the highest for different stages. As the size of the standard deviation can reflect the distinguishing degree of each phenological period, the PGI proposed in this paper is more suitable for describing the differences between phenological periods than the other characteristics.

The above experiments show that the fusion parameter simulated by the growth curve achieves a good classification effect. To study the fusion parameter (FP) and polarization growth characteristic (PGI) time-series changes in more detail, six sets of data were randomly selected from the rice dataset to construct fusion parameter time-series diagrams and polarization growth characteristic time-series diagrams, as shown in Figure 18, where colors represent different sets of data.

The fusion parameter is further processed through goal planning to obtain the polarization growth feature timing diagram (Figure 18b). Compared with Figure 18a, it can be found that the overall trend of the polarization growth characteristics keeps rising, but the changes of the polarization growth characteristics are different in each group of curves. There are two main reasons for this. First, the polarized growth characteristics is based on the goal planning method to measure the optimal weight. Therefore, the time-series of the polarized growth index of rice is still slightly different from its ideal growth curve trend. Second, although the backscattering coefficients of the same category have similar changes in size, the specific sizes are not the same, resulting in different values of the polarization growth characteristics after target planning. In addition, because the training set is data from different regions, this may cause the sowing and maturity times to not be strictly consistent, so the scattering methods of the crops are not exactly the same, which will bring certain deviations to the data calculation process.

Through the above experimental analysis, it can be seen that, compared with the other input characteristics, the polarization growth curve obtained by using the polarized growth index (PGI) in this paper can have a good classification effect on the rice phenological interval, and the phenology classification accuracy could reach 93% in 2021. This is because this study adopted the strategy with the least errors when using the goal plan to determine the weight. At the same time, the setting of the weight in the feature fusion process can better characterize its growth law. However, the polarization growth characteristic has certain limitations. For example, the relationship between $S_{VH}$, $S_{VV}$, $H$, $\alpha$ and the weight is not stable. This results from the large fluctuations in the backscattering coefficient at different stages of rice growth. For example, the backscattering coefficient is greatly affected by soil vegetation. Taking the experimental data as an example, the values of the backscattering coefficients of the two rice test areas are inconsistent, but they have similar changing trends.

To study the usefulness of the rice phenological information obtained by remote sensing has a certain guiding role for the development of the photovoltaic agriculture industry, the measured and estimated insolation and temperature conditions for crop growth were compared, as shown in Figure 19. From the macroscopic point of view, light is the main source of energy for photosynthesis in crops and plays a decisive role in promoting or preventing the reproductive growth of plants. The heat requirements of the crop are reflected in the accumulated temperature, which only accumulates to a certain level before the crop is able to grow properly. From Figure 19a,b, it can be seen that the change curves of the light and temperature accumulation data required by rice inferred from satellite data have the same trend as the actual change curves, and the estimated data are generally higher than the actual data. This is mainly due to the small average error between the rice phenological information obtained based on the proposed method and the phenological information from the field survey. More accurate phenological information can help photovoltaic agricultural practitioners to design different sizes of greenhouses for different crops.

In addition, the rice phenology evolution curve obtained based on the proposed method has the same trend as the field survey phenology evolution curve, with an average error of 3.08 days. This is mainly due to the high time resolution of the Sentinel-1 satellite and the superiority of the proposed PGI in phenological distinction. Since our extracted phenology results rely on satellite data, and the time resolution of the SAR sensor is relatively low, the satellite transit time cannot accurately correspond to the rice phenological transition time. Hence, there is bound to be an error between the obtained and real results. Küçük [21] obtained the results for the identification of six phenological periods using machine learning based on multi-temporal data from TerraSAR-X, with an overall accuracy of about 85%. Ze He [18] applied a decision tree classifier based on RADARSAT-2 to distinguish four phenological phases, with an overall accuracy level of 86.2%. De Bernardis [29] estimated that the key dates used a dual-polarization SAR time series to estimate the key dates and stages of the rice crop, verifying that the prediction error of the critical dates was within 5 days for about 50% of the fields and within 10 days for about 85% of the fields. Compared with the results of the above-mentioned studies, it is generally seen that the method proposed in this paper has achieved satisfactory results with less errors in identifying rice phenological periods.

From another point of view, this study needs to be further developed. To verify the stability of this method for crop classification in different regions, the test and training samples should have different positions. The proposed method is effective when the positions of the training and test samples are not significantly different. However, this study lacked sufficiently long time series experimental observations and was conducted only for specific crops at the study site. If the training and verification areas are geographically different, the application effect of the method will be less. This is because the differences in climate and environment across regions cause the same crop to have different backscatter levels, and there may be mismatches in goal planning. Crop phenology studies at the regional and higher scales still need to be further developed.

On the other hand, this study mainly focused on the relationship between temperature and light in response to phenology. Crop phenology is actually influenced by a combination of multiple factors such as variety, climate and agronomy [45,46,47]. If only the influence of a single climatic factor is considered, and the interaction between the climatic factors and the role of management factors such as agronomy are ignored, the uncertainty of the study results will be increased. Future research will consider further detailing the response of crop phenology to changes in the key climate factors, exploring the drivers of the phenological changes and their influence mechanisms. This will help to understand the response mechanism of crop phenology to the changes in the climate and management measures, providing a scientific basis for the planning and design of photovoltaic agricultural projects. It is of practical significance for understanding the impact of global changes on crop phenology and guiding regional agricultural production.

6. Conclusions

This study proposed a phenological period retrieval method based on growth curve simulations. Based on the data of Jinhu and Liyang of Sentinel-1, this study referred to the growth curve of rice to construct a polarized growth index (PGI) that could approach the growth of rice, and used it to divide and extract its phenological interval and date. The experimental results showed that, compared with the backscattering coefficient characteristics, polarization decomposition characteristics and other parameters of growth curve simulation, the polarized growth index (PGI) had the highest classification accuracy through the weighted fusion of each characteristic. The average accuracy of each period was about 90%, and the overall accuracy of rice phenology reached 96%. The average errors between the date of the key phenological period of rice extraction and that obtained by field ground investigation was 3.08 days.

Since the ideal growth trend of the same crop is similar, the method and the proposed parameters proposed can be applied to the phenological period retrieval of the same crop in different regions. In addition, an appropriate amount of multitemporal data and representative training samples are the basis for establishing a reliable growth model. As the number of images increases, the selection of information for rice growth scenes becomes easier. Therefore, the proposed method requires multitemporal images to cover the crop growth cycle as much as possible.

With the advent of the new energy revolution, the photovoltaic power generation industry has now become an international trend. Therefore, the growth curve simulation method in this paper achieved satisfactory results in the division and extraction of crop phenology, providing a basis for the planning and design of photovoltaic agricultural projects adapted to climatic phenomena and changing environmental conditions. It can greatly promote the development of photovoltaic agriculture, thus contributing to the realization of agricultural science and technology and agricultural industrialization, with more obvious economic, social and environmental benefits.