Domain-guided Machine Learning for Remotely Sensed In-Season Crop Growth Estimation

Advanced machine learning techniques have been used in remote sensing (RS) applications such as crop mapping and yield prediction, but remain under-utilized for tracking crop progress. In this study, we demonstrate the use of agronomic knowledge of crop growth drivers in a Long Short-Term Memory-based, domain-guided neural network (DgNN) for in-season crop progress estimation. The DgNN uses a branched structure and attention to separate independent crop growth drivers and capture their varying importance throughout the growing season. The DgNN is implemented for corn, using RS data in Iowa for the period 2003-2019, with USDA crop progress reports used as ground truth. State-wide DgNN performance shows significant improvement over sequential and dense-only NN structures, and a widely-used Hidden Markov Model method. The DgNN had a 4.0% higher Nash-Sutfliffe efficiency over all growth stages and 39% more weeks with highest cosine similarity than the next best NN during test years. The DgNN and Sequential NN were more robust during periods of abnormal crop progress, though estimating the Silking-Grainfill transition was difficult for all methods. Finally, Uniform Manifold Approximation and Projection visualizations of layer activations showed how LSTM-based NNs separate crop growth time-series differently from a dense-only structure. Results from this study exhibit both the viability of NNs in crop growth stage estimation (CGSE) and the benefits of using domain knowledge. The DgNN methodology presented here can be extended to provide near-real time CGSE of other crops.


Introduction
The increase in synchrony of global crop production and frequency of climate changedriven abnormal weather events is leading to higher variance in crop yields [1][2][3]. Most staple food crops are more vulnerable to yield loss in specific stages of growth and as such, accurate crop growth stage estimation (CGSE) is vital to track crop growth at different spatial scaleslocal, regional, and national -and anticipate and mitigate the effects of variable harvest. Highresolution Remote Sensing (RS) data have been successfully employed to track crop growth at regional scales, however current methods for CGSE utilize curve-fitting and simplistic Machine Learning (ML) models that cannot describe the more complex relationships between crop growth drivers and crop growth stage progress [4][5][6][7][8]. Many of these methods require full-season data and do not provide in-season CGSE information.
Advanced ML models have found success in applications such as crop-cover mapping/classification and yield estimation [9][10][11][12][13], but these models have yet to been applied for in-season CGSE. Whereas methods such as Neural Networks (NNs) have been used in crop mapping, for which researchers can utilize crop cover maps [14] to retrieve millions of crop cover examples per year, field-level crop growth stage (CGS) data is not publicly available and producing field scale ground truth data via field studies is prohibitively expensive. As such, arXiv:2106.13323v2 [cs.LG] 22 Sep 2021 large scale CGSE research relies on local and regional level crop progress data for ground truth. Even for the longest continually running sub-weekly temporal resolution remote sensing (RS) sensors (e.g. MODIS), there are only 21 full growing seasons of crop growth data. Constructing accurate, in-season ML approaches from such limited data is difficult, particularly with few example seasons of abnormal weather. In addition, many crop growth studies estimate events such as 'start of season', 'peak of season', 'end of season' etc. [6] [7], even though these events do not really describe phenological progress and knowledge of their timing may not be actionable.
Recently, domain knowledge has been used to improve the performance of ML techniques in applied research using techniques collectively known as Theory-guide Machine Learning (TgML) [15] [16]. TgML techniques include the use of physical models outputs [17], the integration of known domain limits into ML loss functions [18], and the designing of NN structures that reflect how variables interact within a real physical system [19]. These techniques have been shown to reduce the amount of data required to reach a given level of performance [15]. These Theory-guided Neural Networks have begun to significantly improve upon current state-of-the-art methods in applications such as [20] [21]. Whereas NNs have shown great promise in agricultural RS studies (e.g. [11] [12]), TgML methods have yet to utilize significant disciplinary advances in agriculture over the last two to three decades. The goal of this study is to understand the impact of incorporating domain knowledge into NN design for in-season CGSE at regional scales. Specifically, the objective of the study is to develop a Domain-guided NN (DgNN) that separates independent growth drivers and compare its performance to sequential NN structures of equivalent complexity. The TgML approach in this study is demonstrated for regional CGSE in field corn, which is one of the most cultivated crops in the world [22]. The methodology here, when paired with adequate crop mapping techniques, can be extended to track in-season growth of other crops.

Study Area and Data
This study was conducted in the state of Iowa, US, from 2003 to 2019. The state consists of nine separate Agricultural Statistical Districts (ASDs) (see Figure 1) and had an average of 13.4 million acres of corn under cultivation across the study period [23].  [23] Location of corn field within the study region were obtained from the Corn-Soy Data Layer (CSDL) [24] for 2003-2007 and from the USDA Crop Data Layer (CDL) from 2008-2019 [25]. In Iowa, corn is typically planted in mid April / early May (week of year (WOY) [15][16][17][18][19][20][21][22][23][24], reaches its reproductive stages around late June (WOY 27 onward), and is harvested from early September through late November (WOY [36][37][38][39][40][41][42][43][44][45][46][47][48]. Weekly USDA-NASS Crop Progress Reports (CPRs), generated from grower and crop assessor surveys, were used as ground truth. CPR progress stages include Planted, Emerged, Silking, Dough, Dent, Mature, and Harvested. In this study, the Planted stage was replaced with Pre-Emergence, a placeholder progress stage that represents all crop/field states prior to emergence, and the Dough and Dent stages were combined as Grainfill. CGSE requires both canopy growth information and meteorological data. This study used ASD-wide means and standard deviations of field-level RS and other data shown in Table 1. Fields in each ASD were selected from the CSDL and CDL based on size and boundary criteria (see Figure 2).   Micro-meteorological observations obtained from DayMet were used to compute fieldlevel accumulated growing degree days (AGDD), which is a measure of accumulated temperature required for crop growth. AGDD is used to model progression through different corn growth stages, both in remote sensing studies and in mechanistic models [4][8] [26]. The total number of growing degree days (GDD) for a single 24 hour period is calculated using the function: where T max is lower value of daily maximum temperature and 34 • C, T min is the minimum recorded daily temperature, and T base is the minimum temperature above which GDD is accumulated, set to 8 • C [27]. AGDD is a running total of daily GDD values, and in this study it is calculated from April 8th of a given year, which is the date prior to first planting during the study period. Solar radiation inputs were converted from W/m 2 to MJ/m 2 /week using day length taken from DayMet to incorporate photoperiod information. Saturated hydraulic  [30] To measure canopy growth, 4-day MODIS Fraction of absorbed Photosynthetically Ac-1 tivate Radiation (FPAR) values for each field within an ASD were filtered to produce daily 2 time series data following the Savitsky-Golay (SG) filter method for NDVI used in [31]. In week vary slightly as the season progresses and more data is included within the long-term 8 filtering window. Noise filter adaptions to the existing SG method included rejecting points 9 with an absolute gradient of > 0.3 from the previous value, prior to September, the earliest 10 harvest over the study period. These adaptions prevented noisy, phenologically unrealistic 11 data from being included in the moving filter window. It should be noted that while this 12 filtering system is effective for a uni-modal crop such as corn, it may not be effective for crops 13 with more complex seasonal FPAR patterns such as winter wheat, where higher polynomial 14 filter parameters may be required. 15

16
Since CPRs are released every Monday from data collected during the prior week, weekly 17 meteorological data were obtained by aggregating field-scale daily data from Monday-Sunday. 18 The dataset spanned 38 weeks, WOY 13-51, encompassing the earliest planting and latest 19 harvest reported during the study period. To simulate in-season monitoring, one time series 20 was produced from pre-emergence to the 'current' week per field, totaling 39 time series. 21 Field-level time series for each input were then aggregated to ASD-level by calculating the 22 mean and standard deviation of the values across each district (median was used for rainfall), 23 with 12 total inputs (Table 1). These ASD-wide means and standard deviations formed the 24 un-scaled data for in the study. Solar radiation and rainfall were standardized using Z-score 25 scaling and AGDD and FPAR were standardized using MinMax scaling. To standardize the 26 length of each time series to 39 weekly values, all Z-scored in-season time series were zero 27 padded, while MinMax-scaled inputs were padded with 0.5. ASD location within Iowa was 28 represented using a one-hot location vector of length 9, with each bit representing an ASD. 29 The complete 17-year dataset consisted of 5967 time series, each of dimension (39 X 12) with 30 accompanying location vector. 31

32
The NNs in this study were based upon Long Short-Term Memory (LSTM) layers [32], 33 that are widely used in sequence identification / classification problems, such as speech 34 recognition, translation, and time series prediction [33]. LSTM has also found success in RS 35 studies, including crop classification [11] and yield prediction [34]. In this study, LSTM was 36 used because of its ability to handle time series data with variable length gaps between key 37 events [35], such as variable in-season crop growth data. Three NN implementations were 38 investigated, two reference structures and a third NN that incorporated domain knowledge. 39 The two reference structures included a dense NN, with traditional dense layers of decreasing 40 size ( Figure 3a) and a sequential NN in which LSTM layers were linearly chained ( Figure 3b). 41 In designing the third NN, interactions between the 12 different inputs (see Table 1) and 42 their effect on crop growth were considered. Domain knowledge was incorporated into the 43 NN by separating inputs into a branched, structure based on their relationship to crop growth. 44 TgML studies suggest that organizing NN inputs to reflect their real world interactions may 45 improve performance [15]. For example, Khandelwal et al. [16] [37], and excess solar radiation or photoperiod is 53 used in mechanistic crop models to determine growth stage timing (e.g. [26]). In addition, 54 soil moisture stress, due to low rainfall, in juvenile corn has been found to delay growth 55 progression and reduce final plant size [38] [39]. Typically, the effects of these two drivers 56 on canopy growth and crop progress are modeled separately [26,40,41]. In this study, solar 57 radiation and soil moisture-related inputs were separated from FPAR and AGDD using a

73
In this study, Kullback-Leibler Divergence (D KL ) was used as the loss function for all NNs. D KL is a measure of the difference between two probability distributions. D KL is often used in NN regression problems with targets that are distributions. Given two distributions P(x)a nd Q(x), D KL is calculated as: Here D KL is used as it provides a measure of the difference between the predicted and actual   ASD-level CGS estimates for the NNs and HMM were aggregated to state-level estimates for comparison via weighted sum, with ASD weights calculated based on the number of corn fields in each ASD that passed the processing criteria, explained in Section 2.1 and shown in Figure 2. Performance of the three NN structures and the HMM were evaluated against state-wide USDA CPRs using two metrics. The first, Nash-Sutcliffe efficiency (NSE), is a measure commonly used in hydrology and crop modeling to measure how well a model describes an observed time series versus the mean value of that time series, and is defined as: given growth stage over time.

110
The second metric, cosine similarity (CS), is a measure of the angle between two vectors in a multi-dimension space. CS between two vectors A and B is calculated as: were visualized for the layers in each NN feeding into the 128 node dense and softmax layers, 123 these being common to all three NNs (see Figure 3). Color representations of crop progress 124 were formed by reducing the six crop stages to three RGB channels using a UMAP reduction 125 to 3 dimensions, with 15 neighbors.         [18][19]. In addition, as seen in Table 6, the DgNN produced 186 the best estimates across all growth stages for the greatest number of weeks in each of the test 187 years, producing the highest CS value for 39% more weeks than the next best NN.         stages. This is expected given the large deviation in crop progress timing from the study 201 mean during that year (see Table 3). With the exception of the Emerged stage, the DgNN best 202 described each of the growth stages, as measured by NSE. The week-to-week CS performance 203 degraded less for the DgNN than the other NNs during the fast-moving WOY [33][34][35][36][37][38][39][40][41][42][43] Mature progression, as seen in Figure 6. progress to rates not seen since 1987 [48]. As shown in Table 3, growth stage time to 50% was  Figure 9.

214
In 2014, Silking NSE for the NNs was low, even though crop progress that year was 215 relative typical of the study period. All three NN structures were late in estimating the onset 216 of the Grainfill stage, as shown in Figure 9 and also reflected in the decline in CS between 217 WOY 30 and 34 (see Figure 7). Cumulative progress estimates, as seen in Figure 9, exhibit stage (see Figure 9). In the absence of measurable canopy change, AGDD is a useful proxy 244 for estimating progress from Silking to Grainfill. AGDD, however, is cultivar specific, and 245 cultivars are selected based on different factors, such as planting timing and drought risk.

246
Cultivar-specific variation in required AGDD for progress through Silking and the short 247 duration of that stage may reduce the effectiveness of AGDD as a proxy. In addition, AGDD

Conflicts of Interest:
The authors declare no conflict of interest.