Modeling of Environmental Impacts on Aerial Hyperspectral Images for Corn Plant Phenotyping

Aerial imaging technologies have been widely applied in agricultural plant remote sensing. However, an as yet unexplored challenge with field imaging is that the environmental conditions, such as sun angle, cloud coverage, temperature, and so on, can significantly alter plant appearance and thus affect the imaging sensor’s accuracy toward extracting plant feature measurements. These image alterations result from the complicated interaction between the real-time environments and plants. Analysis of these impacts requires continuous monitoring of the changes through various environmental conditions, which has been difficult with current aerial remote sensing systems. This paper aimed to propose a modeling method to comprehensively understand and model the environmental influences on hyperspectral imaging data. In 2019, a fixed hyperspectral imaging gantry was constructed in Purdue University’s research farm, and over 8000 repetitive images of the same corn field were taken with a 2.5 min interval for 31 days. Time-tagged local environment data, including solar zenith angle, solar irradiation, temperature, wind speed, and so on, were also recorded during the imaging time. The images were processed for phenotyping data, and the time series decomposition method was applied to extract the phenotyping data variation caused by the changing environments. An artificial neural network (ANN) was then built to model the relationship between the phenotyping data variation and environmental changes. The ANN model was able to accurately predict the environmental effects in remote sensing results, and thus could be used to effectively eliminate the environment-induced variation in the phenotyping features. The test of the normalized difference vegetation index (NDVI) calculated from the hyperspectral images showed that variance in NDVI was reduced by 79%. A similar performance was confirmed with the relative water content (RWC) predictions. Therefore, this modeling method shows great potential for application in aerial remote sensing applications in agriculture, to significantly improve the imaging quality by effectively eliminating the effects from the changing environmental conditions.


Introduction
Recent years have seen the rapid growth of remote sensing applications in the field of agriculture [1][2][3]. The advent and advances of low-weight and low-cost imaging platforms, and smart imaging devices resulted in the improved capability of the agricultural data collection. With various sensors such as red-green-blue (RGB), hyperspectral, and thermal cameras carried by these platforms, plant phenotypic properties are captured in images that largely facilitate the process of crop analyses, such as accessing plant biomass, nutrient level, diseases stresses, etc. [1,[4][5][6][7][8][9][10]. However, the changing environmental conditions have been reported to significantly impact the imaging results [11]. The intensity of remotely sensed images changes greatly based on when and where the image was captured [12][13][14].
One source of the variation arises from the complicated interactions between the camera's sensitivity, camera's view angle, plant canopy geometry, solar zenith angle, solar azimuth angle, and shadows [15][16][17][18]. Another source of variation results from plants' endogenous stress responses to the environmental conditions, with complicated interactions between their genetic backgrounds, external environments, and treatments [15,19]. All of these, collectively regarded as the environment-induced variation in phenotyping features, affect plants' final reflectance characteristics.
To reduce the impacts from environment variation, a relatively simple method involves standardizing or fixing the sampling time of the day and restricting imaging to clear conditions without cloud coverage [3,20,21]. Bellvert and Girona [22] suggested that the field phenotypic data should be collected around noon under clear weather conditions. Similarly, unmanned aerial vehicles (UAV) imaging is preferably performed at midday to ensure consistent data collection and analysis [22]. These restrictions could reduce the environmental impacts on the data, but they also significantly inhibit the imaging window and flexibilities. Practically speaking, performing the imaging at a fixed time is difficult, as many procedures, such as equipment setup, need to be completed before imaging and the field environment is naturally uncontrollable and unpredictable. These challenges often result in the collection of remote sensing data at different times of the day under varying environmental conditions. Therefore, the correction of the impact of different imaging time and varying environmental conditions to improve the quality of agricultural remote sensing is critically important to study.
Most aerial remote sensing systems require imaging white reference panels beside the target plants so that the sensing data are calibrated against reference values to remove the illumination variation between images [23]. White reference calibration is effective in compensating different lighting conditions, but these reference images do not precisely reflect the bidirectional reflectance (BRDF) on the leaf surface. Variations from the changes in leaf angles and the plant's endogenous responses remain. An existing calibration method, the combined PROSPECT leaf optical properties model and SAIL canopy bidirectional reflectance model (PROSAIL) [24], enables the prediction of the plant canopy spectral reflectance changes caused by the changing environmental conditions [17,24]. However, the model usually does not meet the accuracy requirement in plant phenotyping remote sensing [25]. Furthermore, the PROSAIL prediction theoretically requires three input variables, including leaf structure parameter, photosynthetic pigment concertation, and water content, which are difficult and costly to measure in remote sensing practices [26].
Another potential solution arises from the use of different regression methods to predict and compensate for the environmental effects. For example, a correction model with the polynomial regression method was developed to predict the crop reflectance as a function of solar zenith angle, time of day, and instantaneous clearness index (ICI). The capability of the model in reducing the diurnal variation with green normalized difference vegetation index (GNDVI) and some individual bands [27] was tested. However, that model only calibrates the imaging time and ambient lighting factors, while many more environmental condition changes, such as temperature and wind speed, also impact the imaging result. Moreover, the plant data was collected on a small portion of the leaf by a handheld radiometer with four bands, which may not properly simulate airborne remote sensing platforms carrying hyperspectral or multispectral cameras. The simple polynomial regression model could successfully describe the changes in data over three consecutive days. However, it may fail to represent the general pattern on other days when the plants are at different stages of their growth cycle. Therefore, a comprehensive environmental impact analysis for general aerial remote sensing images is still critically needed. This analysis requires the continuous collection of crop images at various plant stages through different environmental conditions, a task that has proven challenging with existing airborne remote sensing systems.
On the modeling method side, the artificial neural network (ANN) models, as opposed to conventional regression models, has received considerable attention because of its ability Remote Sens. 2021, 13, 2520 3 of 17 to learn the features directly from the raw data without prior knowledge and human effort in feature design [28]. Due to their better data utilization capacity, ANN models have outperformed conventional methods for solving regression problems in many ways [29]. For example, researchers have achieved high accuracies and efficiencies on modeling multivariable and time-series datasets [30][31][32]. Given the previous successful applications, an ANN model can prove a reliable and efficient alternative for modeling the environmentinduced effects in remote sensing data.
This article introduced the research work of correcting the aerial remote sensing results by modeling and analyzing the effects from the changing field environmental conditions, such as sun radiation, solar zenith angle, humidity, temperature, and wind speed. There are three major objectives in the work of this article:

1.
Collect time-series hyperspectral images of two varieties of corn plants with three nitrogen treatments from V4 to R1 every 2.5 min throughout the whole growing season, along with synchronized environmental condition data.

2.
Build a prediction model for the environment-induced variation in each of the measured phenotyping features (e.g., NDVI and RWC) with time-series decomposition and ANN method.

3.
Evaluate the performance of the trained ANN models and their effects in removing the environmental noise by comparing the variances in the phenotyping features (e.g., NDVI and RWC) before and after the model correction.

Experiment Design and Data Collection
To analyze the environmentally induced variation in phenotyping data, hyperspectral images of the crops and environmental data were collected from a corn field in the Purdue University Agronomy Center for Research and Education (ACRE). Two genotypes (B73 × Mo17 and P1105AM) of corn (Zea mays L.) were grown in the summer of 2019. Each genotype was treated with three different nitrogen (N) solutions: high N with 56 kg/ha (32 mL, 28-0-0 in 1L water), medium N with 28 kg/ha (16 mL, 28-0-0 in 1L water), and low N with 0 kg/ha (water). In total, six experimental plots existed, with one of two corn genotypes treated with one of three nitrogen levels; each plot had around 25 plant replicates. The abbreviation for each experimental plot is listed in Table 1. Table 1. Abbreviations of plant plots with different nitrogen treatments and genotypes.

Plant Groups
Genotypes N Treatments Abbrev P1105AM (Genotype 2) Low G2L 1-6 combined All combined All combined All Hyperspectral images of corn plants were continuously acquired using the Purdue field VNIR hyperspectral imaging gantry system [33]. To capture the instant environmental effects on the images, imaging frequency was set at 2.5 min. Starting from the vegetative growth stage, V4, the continuous imaging occurred for 31 consecutive days until the plants reached the reproductive stage, R1 ( Figure 1). Every day, imaging started at 8:00 am and ended at 7:30 pm. In total, we collected 8631 hyperspectral images of the same crop field ( Table 2) for this study. After data collection, the hyperspectral images were further processed to measure the plant phenotyping features of interest, including the reflectance spectrum, NDVI, and predicted RWC for each experimental plot. The reflectance spectrum was obtained from the averaged data of plant tissues using the image segmentation algorithm highlighted in [34]. The NDVI was calculated from the spectrum by following Equation (1) [35,36]. The plant's RWC was predicted with the pretrained partial least squares regression (PLSR) model [33].
Environmental data 31 8631 Air temperature (°C) Sun radiation (W/m 2 ) Wind speed (m/s) Solar zenith angle (degree) Humidity (%) Diurnal time (min) In addition to the hyperspectral imaging data, a local miniature weather station (Ambient Weather, Chandler, AZ, USA) was installed within the experimental plots to collect real-time time-tagged environmental data. The environmental data included air temperature (℃), solar radiation (W/m 2 ), wind speed (m/s), sun zenith angle (degree), humidity (%), and diurnal time (min) ( Table 2).

Time Series Decomposition for Environment-Induced Variation
The phenotyping data from 31 days were collected to provide enough images under various environmental conditions. However, besides the instant environmental effects,  In addition to the hyperspectral imaging data, a local miniature weather station (Ambient Weather, Chandler, AZ, USA) was installed within the experimental plots to collect real-time time-tagged environmental data. The environmental data included air temperature (°C), solar radiation (W/m 2 ), wind speed (m/s), sun zenith angle (degree), humidity (%), and diurnal time (min) ( Table 2).

Time Series Decomposition for Environment-Induced Variation
The phenotyping data from 31 days were collected to provide enough images under various environmental conditions. However, besides the instant environmental effects, the plant growth change and other day-by-day gradual weather changes also contributed to the variation among the images. These different components of variation need to be clearly separated before we can focus on modeling the instant environmental effects. As most of the field environment factors fluctuate over the course of a single day [37,38], we hypothesized that the higher frequency environment-induced variation could be identified by removing the lower frequency variation as the day-to-day trend. Thus, a time series decomposition method was applied, decomposing the original time-series phenotyping signal into two Remote Sens. 2021, 13, 2520 5 of 17 major parts: day-to-day trend (T t ) and daily instant changes (D t ) (Equation (2)). More specifically, T t is calculated with the locally estimated scatterplot smoothing (LOESS) method [39]. By fitting a non-parametric regression curve on the scattered plot of the data, the day-to-day changing trend can be extracted from the raw signal [40]. This trend mostly reflects the changes in the plant growth stage and general weather conditions over the 31 days of imaging. The daily instant changes (D t ) were calculated by subtracting the day-to-day trend (T t ) from the raw signal. D t contains the higher frequency variation components mostly caused by the plant's circadian behavior and environmental condition changes such as sun angle, solar radiation, and temperature changes during the day. In this study, D t was used as the output of the proposed model.

Environmental Data Transformation and Selection
To generate more discriminatory power in higher-dimensional feature spaces besides the original environmental variables (temperature, solar zenith angle, wind speed, etc.) for improved model accuracy, a feature transformation was performed by taking the square and square root of the measurements of the original environment factors [41]. These new non-linear variables have proven more effective in modeling environment variation [27]. Finally, the transformed variables were combined with the original variables for further processing.
After transformation, the features were selected to remove the irrelevant input of some environmental variables to reduce overfitting. This also lowered the difficulty of future applications, with fewer measurements required. A single-factor correlation analysis was performed. Each of the original or transformed environment variables was fitted with the calculated environment-induced variation in phenotyping data (daily instant changes D t in Equation (2)) in a linear regression model. The adjusted R 2 , which indicates the relevance of each feature to the estimated environmental variation, was calculated. A higher adjusted R 2 meant a variable was more relevant [42]. By comparing the adjusted R 2 of each variable, we determined the final list of input environmental variables for the model.

Data Quality Check
Training data quality is critically important for a supervised machine learning model. The data quality was checked before training the model, and the outlier data was removed [43]. For each phenotyping feature (NDVI and predicted RWC, etc.), the daily measurements between the upper inner fence (Q3+1.5IQR) and lower inner fence (Q1-1.5IQR) were kept [43]. IQR is the interquartile range, which equals the difference between the 75th (Q3) and 25th (Q1) percentiles. Meanwhile, image data before 10:00 am and after 5:30 pm was also removed, as it demonstrated extreme variance and noise [33]. Using the NDVI as an example, we employed the training data sizes shown in Table 3 to train the ANN model. The architecture of the proposed model is based on a feed-forward multi-layer perceptron (MLP) network, a class of ANN models ( Figure 2). Due to their adjustable architecture, MLP models are particularly flexible. This flexibility increases the suitability of MLP for regression prediction problems where a real-valued quantity is predicted, given a set of inputs such as time-series data [44]. In this study, after some speed-accuracy tradeoff pretests on model performance, the proposed ANN architecture was configured with a four-layer model containing an input layer, two hidden layers, and an output layer. After each hidden layer, the Leaky ReLU activation was performed to add non-linear properties to the function [45]. The selected environmental variables served as input for the model, whereas the D t environment-induced variation of selected phenotyping features was the output. To accelerate learning and lead to a faster convergence, both input and output data were normalized for modeling purposes [46], while the final prediction results were denormalized back to the original scale of the phenotyping feature. The architecture of the proposed model is based on a feed-forward multi-layer perceptron (MLP) network, a class of ANN models ( Figure 2). Due to their adjustable architecture, MLP models are particularly flexible. This flexibility increases the suitability of MLP for regression prediction problems where a real-valued quantity is predicted, given a set of inputs such as time-series data [44]. In this study, after some speed-accuracy tradeoff pretests on model performance, the proposed ANN architecture was configured with a four-layer model containing an input layer, two hidden layers, and an output layer. After each hidden layer, the Leaky ReLU activation was performed to add non-linear properties to the function [45]. The selected environmental variables served as input for the model, whereas the environment-induced variation of selected phenotyping features was the output. To accelerate learning and lead to a faster convergence, both input and output data were normalized for modeling purposes [46], while the final prediction results were denormalized back to the original scale of the phenotyping feature.

Training and Optimization
To train the network with minimum overfitting, the training process followed a fivefold cross-validation scenario. We randomly divided the whole dataset into five roughly equal subsets. In the first iteration, the first subset was used to test the model and the rest aided in training the model. This process was repeated until each subset has been used as

Training and Optimization
To train the network with minimum overfitting, the training process followed a fivefold cross-validation scenario. We randomly divided the whole dataset into five roughly equal subsets. In the first iteration, the first subset was used to test the model and the rest aided in training the model. This process was repeated until each subset has been used as the testing set. During training, the loss function was Mean Square Error (MSE). The network was initialized with the Kaiming weights [47]. All the ANN models were trained using the Adam optimizer [48].
The accuracy of the model was optimized by adjusting the learning rate, batch size, and the number of epochs. The learning rate controls the rate or speed at which the model learns [49]. The batch size establishes the accuracy of the error gradient estimate when training neural networks [50,51]. The number of epochs impacts the ability of the model to be generalized by determining how many times the model trains on the same data. Finally, by comparing the accuracy (R 2 and RMSE) of models with different combinations of the learning rate, batch size, and number of epochs, the model parameters with a batch size of 600 for 120 epochs with learning rate at 1e-3 were chosen for this study.

Evaluation Metrics
The performance of all the developed models was evaluated and compared with the coefficient of determination (R 2 ) and root mean square error (RMSE) between the prediction results and the original measurements. Meanwhile, we also compared the daily variances of the selected phenotyping features (e.g., NDVI) before and after the model correction. A two-sample t-test was performed to check if daily variance in features fell significantly.

Multi-Model Comparison Analysis across Genotypes and Nitrogen Treatments
The impacts of nitrogen treatments and genotypes on the ANN modeling were also investigated to determine whether a separate ANN model was needed for each case or if one general ANN model could fit the different treatments and genotypes. The ANN model of each experimental plot was trained separately (Table 3) and was tested on the other treatments and genotypes. We also built a general ANN model containing the entire sample data to provide a unified and general "all-in-one" correction approach. With the group-to-group cross-validation on each of the datasets, the R 2 and RMSE of each model's performance in the other datasets were examined to evaluate the generalization of models across genotypes and nitrogen treatments. For example, if the ANN model (ALL) resulted in similar outcomes as the individual plot models (G1H, G1M, G1L, G2H, G2M, G2L) for each plant plot, this unified model would be adopted. Otherwise, different models should be adopted separately for each different case. The aim was to find the most appropriate correction solution as the best balance between ease of use and accuracy.

Phenotyping Features for Testing the Model and Workflow
To demonstrate the detailed modeling procedures and performance evaluation, NDVI was chosen as the first example as it represents one of the most common plant features in remote sensing [52]. We also then tested the same ANN architecture and workflow on the RWC to validate the generalization of the proposed method.

Software and Computation
All the imaging processing work was implemented with Matlab R2018a software [53]. The modeling work was performed in the Python version 3.7.2 software environment [54]. The ANN model was implemented in PyTorch 0.4.1 [55]. The time-series data was analyzed and manipulated using Pandas [56] and Numpy [57]. The figures were drawn with Seaborn [58] and Matplotlib [59]. The Matlab and Python computations were all executed on a ThinkPad workstation P300 (Lenovo PC international, Morrisville, Morrisville, NC, USA) equipped with 16-gigabytes (GB) of random-access memory (RAM), a 3.70 GHz Intel ® Xeon™ E1270 processor, and Nvidia GTX 1070 GPU.

Time Series Decomposition Result
The time-series data of raw NDVI was successfully decomposed into the day-to-day trend (T t ) and daily instant changes (D t ) (Figure 3). The raw NDVI plot (row 1 in Figure 3) captured the variation in NDVI over the daytime period, with gaps indicating the time between 5:30 in the afternoon until 10:00 next day without imaging data. The raw NDVI plot showed a clear and repetitive V-shaped pattern for each day, which was caused by environment variation during imaging. The day-to-day trend, T t (row 2 in Figure 3), represented the changes of plant growth stage and plant health conditions. As plants mature, the NDVI was expected to increase until the reproductive stage. Meanwhile, the two big dips along the T t curve precisely captured the impacts from two severe temperature drops in the West Lafayette area. This kind of long-term environmental impact was not included in the proposed analysis.
periodical changes with V shapes. Due to the extreme weather conditions, parts of the data were missing on DAP 36, 37, 58, 59, and 60. Overall, remained substantially consistent through 31 days, while amplitude and minor skewness differences existed among the from different imaging days. For example, the of DAP 42 demonstrates a smaller amplitude than that of DAP 56. These differences were caused by different environmental conditions, which would be addressed by the environmental correction model in this study. Figure 3. The NDVI of the sample dataset (G1H) from the V4 stage to the R1 stage. The raw NDVI plot was decomposed into the day-to-day trend ( ) and the periodic change ( ). The red boxes are the days with incomplete data measurements due to the extreme weather conditions, which are DAP 36, 37, 58, 59, and 60.

Environmental Feature Selection and Range
The results of the single-factor correlation analysis for NDVI are shown in Figure 4. The environmental variables were all correlated with the environment-induced variation Severe temperature drops Figure 3. The NDVI of the sample dataset (G1H) from the V4 stage to the R1 stage. The raw NDVI plot was decomposed into the day-to-day trend (T t ) and the periodic change (D t ). The red boxes are the days with incomplete data measurements due to the extreme weather conditions, which are DAP 36, 37, 58, 59, and 60.
On the other hand, the daily instant changes (D t ) (row 3 in Figure 3) showed clear periodical changes with V shapes. Due to the extreme weather conditions, parts of the data were missing on DAP 36, 37, 58, 59, and 60. Overall, D t remained substantially consistent through 31 days, while amplitude and minor skewness differences existed among the D t from different imaging days. For example, the D t of DAP 42 demonstrates a smaller amplitude than that of DAP 56. These differences were caused by different environmental conditions, which would be addressed by the environmental correction model in this study.

Environmental Feature Selection and Range
The results of the single-factor correlation analysis for NDVI are shown in Figure 4. The environmental variables were all correlated with the environment-induced variation in NDVI, except for humidity and its derivatives. The adjusted R 2 values for humidity were almost 0, indicating no correlation found between humidity and NDVI changes. This confirmed the findings from the previous literature that while plant-sensing data were strongly impacted by environment factors such as air temperature, solar radiation, sun zenith angle, and diurnal time [17,27,60], humidity was rarely reported to demonstrate a similar impact. Thus, we removed humidity and its derivatives in the model. Finally, the input feature for each modeling sample was a 1-by-15 vector consisting of air temperature, solar radiation, wind speed, solar zenith angle, diurnal time, and their square or square root values. Besides, the range in environmental conditions experienced by the modeling data was shown in Table 4. These ranges illustrate the appropriate environmental condition to apply the model. zenith angle, and diurnal time [17,27,60], humidity was rarely reported to demonstrate a similar impact. Thus, we removed humidity and its derivatives in the model. Finally, the input feature for each modeling sample was a 1-by-15 vector consisting of air temperature, solar radiation, wind speed, solar zenith angle, diurnal time, and their square or square root values. Besides, the range in environmental conditions experienced by the modeling data was shown in Table 4. These ranges illustrate the appropriate environmental condition to apply the model.

Overall Performance
The R 2 and RMSE measure the precision of the predicted environment-induced variation in NDVI. The environment-induced variation predicted by the ANN model for the sample dataset (G1H) showed a fairly accurate linear relationship with the coefficient of determination (R 2 ) equal to 0.823 ( Figure 5). The RMSE also demonstrated a relatively low value of 0.00611. The prediction result was five-fold cross-validated.

Overall Performance
The R 2 and RMSE measure the precision of the predicted environment-induced variation in NDVI. The environment-induced variation predicted by the ANN model for the sample dataset (G1H) showed a fairly accurate linear relationship with the coefficient of determination (R 2 ) equal to 0.823 ( Figure 5). The RMSE also demonstrated a relatively low value of 0.00611. The prediction result was five-fold cross-validated.
The predicted environmentally induced variation was further used to correct the noise caused by environmental effects in the raw NDVI signal. Figure 6 shows the NDVI corrected by subtracting the predicted variation ( Figure 5) from the raw NDVI. In Figure 6b, each box represents the NDVI changes within a day. The trained ANN model largely eliminated the daily variance in the NDVI, so the boxes of the corrected NDVI (Figure 6b) were much more condensed compared with the original NDVI (Figure 6a). To facilitate the comparison, we compared the variances of NDVI before and after model correction with a two-sample t-test (Figure 7 and Table 5). The result confirmed that the daily variances in NDVI were significantly reduced (p-value < 0.01) by 79% on average, thereby confirming the ability of the proposed ANN model to effectively eliminate the environmentally induced effects on the raw signal. The predicted environmentally induced variation was further used to correct the noise caused by environmental effects in the raw NDVI signal. Figure 6 shows the NDVI corrected by subtracting the predicted variation ( Figure 5) from the raw NDVI. In Figure  6b, each box represents the NDVI changes within a day. The trained ANN model largely eliminated the daily variance in the NDVI, so the boxes of the corrected NDVI (Figure 6b) were much more condensed compared with the original NDVI (Figure 6a). To facilitate the comparison, we compared the variances of NDVI before and after model correction with a two-sample t-test (Figure 7 and Table 5). The result confirmed that the daily variances in NDVI were significantly reduced (p-value < 0.01) by 79% on average, thereby confirming the ability of the proposed ANN model to effectively eliminate the environmentally induced effects on the raw signal.    The predicted environmentally induced variation was further used to correct the noise caused by environmental effects in the raw NDVI signal. Figure 6 shows the NDVI corrected by subtracting the predicted variation ( Figure 5) from the raw NDVI. In Figure  6b, each box represents the NDVI changes within a day. The trained ANN model largely eliminated the daily variance in the NDVI, so the boxes of the corrected NDVI (Figure 6b) were much more condensed compared with the original NDVI (Figure 6a). To facilitate the comparison, we compared the variances of NDVI before and after model correction with a two-sample t-test (Figure 7 and Table 5). The result confirmed that the daily variances in NDVI were significantly reduced (p-value < 0.01) by 79% on average, thereby confirming the ability of the proposed ANN model to effectively eliminate the environmentally induced effects on the raw signal.   . Two-sample t-test between the variance of daily NDVI for the sample dataset (G1H) before and after ANN model correction.

Multi-Model Comparison Analysis across Genotypes and Nitrogen Treatments
The ANN models built for each dataset were tested on the other datasets to evaluate the drifts between different genotypes and nutrient treatments. For datasets from a different genotype or treatment, the ANN model demonstrated a weaker prediction performance compared to the results on dataset it had been trained with ( Figure 8). Notably, the predictions were least accurate when the ANN models trained with nitrogen-stressed plots (G1M, G2M, G1L, and G2L) were applied to the high-nitrogen groups (G1H and G2H), as shown within the red boxes in the Figure 8. The results of the multi-model comparison indicate that the nitrogen stress levels on plants should be considered when modeling the environment-induced variation in phenotyping features. Compared to the nitrogen treatments, the genotype difference demonstrated a minor impact. The R 2 between plots with the same nitrogen treatment but different genotypes were between 0.59-0.79 with RMSE between 0.009-0.013. The general ANN model (ALL) trained with the entire sample data performed well across the different genotypes and treatments with substantially high R 2 (0.617-0.843) and low RMSE values (0.008-0.0010). This allowed us to apply the same one ANN model (ALL) for diverse corn stages (already included in the modeling), genotypes, and treatments (validated in Figure 8).

Multi-Model Comparison Analysis across Genotypes and Nitrogen Treatments
The ANN models built for each dataset were tested on the other datasets to evaluate the drifts between different genotypes and nutrient treatments. For datasets from a different genotype or treatment, the ANN model demonstrated a weaker prediction performance compared to the results on dataset it had been trained with ( Figure 8). Notably, the predictions were least accurate when the ANN models trained with nitrogen-stressed plots (G1M, G2M, G1L, and G2L) were applied to the high-nitrogen groups (G1H and G2H), as shown within the red boxes in the Figure 8. The results of the multi-model comparison indicate that the nitrogen stress levels on plants should be considered when modeling the environment-induced variation in phenotyping features. Compared to the nitrogen treatments, the genotype difference demonstrated a minor impact. The R 2 between plots with the same nitrogen treatment but different genotypes were between 0.59-0.79 with RMSE between 0.009-0.013. The general ANN model (ALL) trained with the entire sample data performed well across the different genotypes and treatments with substantially high R 2 (0.617-0.843) and low RMSE values (0.008-0.0010). This allowed us to apply the same one ANN model (ALL) for diverse corn stages (already included in the modeling), genotypes, and treatments (validated in Figure 8

Modeling of Environmentally Induced Variation in Predicted RWC
Besides NDVI, the environmentally induced variation in the predicted RWC was also modeled and predicted. In Figure 9, the predicted and the measured variation were strongly correlated, with a R 2 of 0.791 and RMSE of 0.722%. With this model, the variance of the corrected RWC was significantly reduced (p-value < 0.01) by 72% on average compared to the raw predicted RWC data ( Figure 10). The successful application of the same proposed ANN architecture and decomposition method on the predicted RWC and NDVI indicated that this method has the potential to be generally applied on other phenotyping features that can be further explored. Moreover, the corrected predicted RWC plot ( Figure  10b) demonstrates a more obvious day-to-day trend than the raw predicted RWC ( Figure  10a). Therefore, with the environmental effects removed, plant remote-sensing researchers can more accurately track the plant growth signals.

Modeling of Environmentally Induced Variation in Predicted RWC
Besides NDVI, the environmentally induced variation in the predicted RWC was also modeled and predicted. In Figure 9, the predicted and the measured variation were strongly correlated, with a R 2 of 0.791 and RMSE of 0.722%. With this model, the variance of the corrected RWC was significantly reduced (p-value < 0.01) by 72% on average compared to the raw predicted RWC data ( Figure 10). The successful application of the same proposed ANN architecture and decomposition method on the predicted RWC and NDVI indicated that this method has the potential to be generally applied on other phenotyping features that can be further explored. Moreover, the corrected predicted RWC plot (Figure 10b) demonstrates a more obvious day-to-day trend than the raw predicted RWC (Figure 10a). Therefore, with the environmental effects removed, plant remote-sensing researchers can more accurately track the plant growth signals. Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 17

Discussion
Environmental change impacts in crop aerial remote sensing images were quantitatively investigated with the proposed ANN modeling approach. Over 8000 hyperspectral

Discussion
Environmental change impacts in crop aerial remote sensing images were quantitatively investigated with the proposed ANN modeling approach. Over 8000 hyperspectral

Discussion
Environmental change impacts in crop aerial remote sensing images were quantitatively investigated with the proposed ANN modeling approach. Over 8000 hyperspectral images of two varieties of corn with three nitrogen treatments were taken by the field imaging gantry at Purdue University over the 2019 growing season. The imaging covered the plant stages from V4 to R1. Crop phenotyping features such as NDVI, RWC, and two individual spectral bands (Red and NIR) were calculated from the imaging data. The proposed ANN model was trained with synchronized hyperspectral imaging data and environmental data (including sun radiation, solar zenith angle, diurnal time, temperature, and wind speed) to understand the correlation of the variations between the two datasets. A time-series decomposition method was applied to extract the phenotyping data variation caused by the changing environments. By learning the relationship between the phenotyping data variation and environmental changes, the developed ANN model was able to precisely predict the environmental effects on remote sensing results (i.e., 82.2% for NDVI), and thus could be used to effectively eliminate the environment-induced variation in the phenotyping features. The two-sample t-tests on the NDVI and predicted RWC of corn plants showed that the daily variances in NDVI and predicted RWC were significantly reduced by 79% and 72%, respectively. Thus, the ANN method can be used by remote sensing professionals to adjust or correct raw imaging data for changing environments to improve plant characterization.
The proposed ANN-based method showed a promising performance in modeling the environment-induced variations in different plant phenotyping features. However, the data used in the model was drawn from one single field test whose imaging data was collected from Purdue's field gantry system, which might induce systematic bias in the model. External validation data from the other remote sensing platforms, such as UAVs, are needed. In the next growing season, the proposed method will be validated with images from a field UAV system as well as RGBN camera-based imaging sensor (Ncam) [61]. Furthermore, this modeling method was developed based on corn images, which might limit the scope of application. It is necessary to conduct more tests on more diverse plant species (e.g., soybean, wheat, and rice). This will help further validate the developed models, as well as improve the robustness of the prospective models.
In the future, we will also continue exploring training models for all the single spectral bands to adjust/correct the whole spectrum data considering environmental variations. Remote sensing users could benefit from the spectrum calibration model to correct the prediction results from any plant feature prediction models.

Conclusions
In this paper, a new modeling method was successfully proposed to precisely predict the environmental effects on the hyperspectral imaging results (such as NDVI and predicted RWC) in aerial crop remote sensing. Over 8000 hyperspectral images, together with synchronized environment data, were collected over 31 days for field corn plants with different nitrogen treatments and genotypes. Experimental results demonstrated that the proposed ANN method could accurately predict the environment-induced variations in the selected phenotyping features. For example, the trained model for NDVI achieved promising predictive results for the sample dataset with an R 2 of 0.822 and an RMSE of 0.00611 compared with the measured variation. The predicted values were used to correct the raw phenotyping data, and the daily variance of NDVI was significantly reduced by 79%. The proposed method also achieved satisfactory results when tested on predicted RWC (daily variance reduced by 72%). The applicability of the proposed method on two different features highlighted its potential to correct the other phenotyping features of interest. Based on these results, this proposed modeling method can help agricultural remote sensing researchers to effectively eliminate the signal drifts caused by the environment variation, which will drastically increase the accuracy of field plant sensing.
Author Contributions: D.M. carried out the data collection, data analysis, and drafted the manuscript; J.J. contributed the research idea, led the designing, fund-raising, and construction of the field imaging facility, designed the experiments, and revised the manuscript; H.M. and M.R.T. guided the experimental design; T.U.R. and L.Z. made substantial contributions to system construction and data collection. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Sumitomo Chemical, grant number 16121941. The APC was funded by Purdue University.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to confidentiality.