Dry Matter Yield and Nitrogen Content Estimation in Grassland Using Hyperspectral Sensor

Abstract

Estimation of Dry Matter Yield (DMY) and Nitrogen Content (NC) in forage is a big concern for growers. In this study, an estimation model of DMY and NC using Visible and Near Infrared (V-NIR) spectroscopy was developed. An adequate number of grass samples (5078) of perennial ryegrass (Lolium perenne), collected from Dutch grassland in 2019 and 2020 were sensed with a hyperspectral sensor, while grass height was recorded in situ by an ultrasonic sensor mounted on a tractor. The samples were treated with Artificial Intelligence (AI) techniques. PCA based feature selection was applied first, revealing that visible green wavelength (around 500 nm) and red edge wavelength (around 700 nm) were enough to express the overall variability of the dataset. Then, Feature Importance analysis of Random Forest Regressor showed that NIR wavelengths (around 910, 960 and 990nm) were the most sensitive in DMY estimation, while red edge (around 710 nm) and visible orange wavelengths (around 610 nm) were the most related to NC estimation. Finally, SHAP (SHapley Additive exPlanations) analysis was applied to the Random Forest estimation models, resulting in the visualization of wavelength selection, thus assisting in the interpretation of the results and the intermediate processes. Overall, this method can lead to the reduction of the number of wavelengths to be measured in the field and thus, to the possible development of a low cost hyperspectral sensor for the above purposes.

1. Introduction

Grassland research is one of the indispensable studies in Europe, given that approximately 34% of agricultural land in the European Union is regarded as permanent grassland [1]. In the Netherlands, United Kingdom, and Ireland, grassland covers 54%, 72%, and 92% of land, respectively [2]. Grassland is not only considered as a need for animal products, but also plays a vital role in the environment, such as minimizing the damage of erosion and flood [3]. Although the multifunctional traits of grassland are well-known in the context of the ecosystem, the area of grassland has declined in Europe over the decades due to land use change, which is related to political transformation and environmental concerns (livestock, methane emission), and intensification [4]. In this context, it is important to find a strategy to improve productivity within the limited grassland area to provide animal products in a more efficient way [3,5].

According to previous studies [6,7], nitrogen use efficiency in grassland is very limited, because more than 50% of the N applied is not assimilated by plants. Furthermore, environmental awareness of high nitrogen input through manure and chemical fertilizers is very important since it could potentially cause deterioration of groundwater and surface water [7].

Another aspect is that in grassland farming assessing Dry Matter yield (DMY) and nitrogen content (NC) is a big concern for farmers. DMY is not only directly related to the amount of pasture but also utilized for post processes such as tedding, raking, baling, or making silage in the form of DMY within a harvested pasture. NC is strongly related to crude protein content within the pasture, and crude protein is a key factor for the total intake and digestion of cattle.

In order to monitor DMY or NC, several approaches have been developed. In general, grass samples are collected in the field and sent to a laboratory to be analyzed by chemical methods. However, this is usually costly and takes time.

There have been studies to estimate DMY and NC using near infra-red spectroscopy. In this area of research, it is important to identify the wavelengths that contribute to the estimation, which will lead to the development of low-cost sensors [8,9,10,11]. One common method for this purpose is principal component analysis (PCA) [12,13,14]. For each feature, the contribution value to the principal components is calculated by PCA. Based on this contribution value, dominant features, i.e., bands in this case, can be selected. Another method for band selection is random forest regressor [15]. Band selection is performed according to the importance of the features calculated by the structure of random forest regressor. Other methods, genetic algorithm (GA), phased regression (PHR), or principal component regressor are also used [16,17]. However, with the above methods, feature importance is expressed as an absolute value, no matter if the correlation is positive or negative. If the transition of the importance of wavelength can be analyzed more precisely, such as when the importance becomes high, low, or even negative, the wavelength for sensors can be selected more reasonably. In addition, it could be useful for a deeper understanding of measured values or results of prediction models, and also for data pre-processing and designing models themselves.

For further utilization of near-infrared spectroscopy in grassland management, this research addresses the improvement of the estimation accuracy of DMY and NC and develops a feature selection method to determine the dominant wavelength. Combining hyperspectral absorbance data with other field measurements, specifically the plant height sensor, focuses on improving the accuracy of DMY and NC estimation. In addition, for more detailed wavelength analysis, several techniques to identify wavelengths which contribute for estimating DMY and NC of pasture are examined. Reducing the number of wavelengths to be measured by hyperspectral sensor might lead to the development of lower cost sensors. Each specific methodology is discussed in following sections.

2. Materials and Methods

2.1. Study Site

An experiment of grassland fields was conducted in the 2019 and 2020 crop seasons from May to October. The dominant grass vegetation was Lolium perenne, one of the perennial ryegrasses widely grown as pasture for livestock not only in the Netherlands but also around the world. The study sites were in Wageningen (51.964°N–51.966°N, 5.632°E–5.634°E) and De Marke (52.040°N–52.039°N, 6.355°E–6.338°E), both a part of The Netherlands. Both fields have 64 plots and each plot has an area of 10 m × 3 m, as shown in Figure 1. The soil type of the field in Wageningen is clay soil, whereas that in De Marke is sandy soil. Annual precipitation (mm/h), minimum temperature (°C), maximum temperature (°C), and average temperature (°C) in Wageningen throughout 2019 and 2020 are 2.3, 7.4, 15.2, and 11.1, respectively, while their statistics in De Marke are 2.2, 7.4, 15.0 and 11.0, respectively.

In Wageningen, harvesting and data acquisition were performed five times in 2019 (30th April, 17th June, 29th July, 3rd September and 29th October) and four times in 2020 (6th May, 22nd June, 11th September and 15th October). In De Marke, harvesting and data acquisition were performed four times in 2019 (13th May, 24th June, 2nd September and 21st October), and five times in 2020 (6th May, 8th June, 17th July, 16th September and 26th October)

In 2019, 32 out of 64 plots were unfertilized, and the remaining 32 fields were fertilized by applying nitrogen. Plots in Wageningen were fertilized in amounts of 140 kg/ha, 100 kg/ha, 60 kg/ha and 40 kg/ha from the 1st cut to the 4th cut, for a total of 340 kg/ha. Plots in De Marke were fertilized in amounts of 140 kg/ha, 90 kg/ha, 81 kg/ha, 54 kg/ha, and 40 kg/ha from the 1st cut to the 5th cut, for a total of 405 kg/ha. On the other hand, different fertilizing methods were employed in order to artificially create different growing conditions in 2020. Specifically, 8 out of 64 plots were unfertilized, 28 plots were equally fertilized, which means the same way as the fertilizing scheme in 2019, and the remaining 28 plots were fertilized in different ways in order to artificially create different growing conditions. For the remaining 28 plots, based on the NC from grass in 2019, more fertilizer was applied to the plots that had more nitrogen.

2.2. Data Collection

In this field experiment, the grass was cut four to five times a year. All plots per location were harvested on the same day using a forage harvester from Haldrup GmbH “https://www.haldrup.net/ (accessed on 9th January 2023)” specially designed for this experiment. The harvesting area was 10 × 1.5 m inside each plot.

First, the weight of harvested fresh grass on each plot was measured. After assessing the fresh weight per plot, approximately 300 to 500 g of grass was randomly sampled in order to determine the DMY. These samples were dried for 48 h at 70 °C and homogenized by grinding in a 0.5 mm mill. Samples were taken from these materials to assess the remaining moisture content by drying at 105 °C. Samples were taken from the remaining materials to assess the total NC by the method of “digestion H₂SO₄-H₂O₂-Se; SFA-Nt/Pt”.

Grass height was measured with a Pasture Reader “http://pasturereader.com.au/ (accessed on 9th January 2023)” which is placed in front of the forage harvester. The Pasture Reader is a sensor that measures the height and the density of the grasses ultrasonic waves. Plant height was measured weekly for each plot.

Before each harvest, crop reflectance was measured by Tec5 HandySpec Field equipment “https://tec5.com/en/products/nir-systems/ (accessed on 9th January 2023)”. It is a spectroscopic sensor system that measures the light absorption of plants or soil in specific wavelength bands based on near-infrared spectroscopy. Measurements are independent of atmospheric influence because of permanent referencing of the ambient light (upper sensor). The spectral data can be used to evaluate biomass, leaf area and nitrogen status. The instrument measures incoming and reflected light simultaneously. The covered wavelength is 400–1000 nm with a spectral resolution of 10 nm. The sensor is held manually at a height of approx. 1.8 m above the top of the canopy and pointed directly downward. The field of view of the sensor is ±25°, which results in a measurement of roughly 0.5 m² of the crop. The sensors used for data collection are shown in Figure 2. The methods and parameters used to measure grass samples from each harvesting area are summarized in

Table 1. Measured Parameters of Grass.

Measured Parameter	Technological Device or Method for Data Acquisition	During Harvest	Lab Analysis
Plant height (cm)	Pasture Reader	x
Fresh Grass Weight (t/ha)	Weight scale	x
Hyper Spectrum reflectance	Tec5 HandySpec Field equipment	x
DMY (t/ha)	Absolutely dry method		x
NC (g N/kg DM)	Digestion H2SO4-H2O2-Se; SFA-Nt/Pt		x

.

Measurements were taken only when the crop was dry. Measurements were taken between 10:00 am and 4:00 pm in order to avoid low solar angles. The operator faced the sun to prevent shading the object. Measurements were made on six different spots within the harvest area of each plot.

Throughout two years of data collection, 2434 and 2644 measurements were collected during harvest in 2019 and 2020, respectively. Histograms of measured data are shown in Figure 3.

DMY and grass height in 2020 were higher than in 2019. There was a lack of precipitation during the crop season in 2020, but optimization of water supply in 2020 with a dripping irrigation system, which could have caused a higher yield. Accumulated rainfall in the pasture growing season in 2019 and 2020 was 474 mm and 411 mm, respectively, in Wageningen, and 480 mm and 425 mm, respectively, in De Marke, i.e., 13% and 11% less, respectively, in 2020. In Wageningen, no irrigation was applied in 2019, while the amount of irrigation in 2020 was 12 mm on 20th May, 15 mm on 27th May, and 15 mm on 10th August, totaling 42 mm. In De Marke, the amount of irrigation in 2019 was 27 mm on 11th July and 34 mm on 31st July, totaling 61 mm, while irrigation in 2020 was 17 mm on 17th April, 18 mm on 13th May, 21 mm on 4th June, 17 mm on 1st August, 22 mm on 7th August, 21 mm on 15th August and 12 mm on 20th September, totaling 129 mm. Here, the growing season was assumed to be from 1st April to 31st October.

2.3. Data Analysis

Using the measurement data collected as the described in previous section, estimation model of DMY and NC is developed. In this section, the methods of the model development and its analysis are explained.

2.3.1. Data Pre-Process

In general, measurement data of living organisms such as plants contain a lot of noise. Therefore, pre-processing of the data to create clean data is indispensable for the application of multivariate analysis, artificial intelligence (AI), etc. to obtain robust and reliable results. For this purpose, the data pre-process procedure shown in Figure 4 was applied.

Details of each procedure are explained in Appendix A. By pre-processing, noise and outliers are removed, the difference in the absolute value of each variable is smoothed out and the useful variables are emphasized, which in turn leads to good data analysis results.

2.3.2. Random Forest Regressor

Random forest regressor [14] is used to build a prediction model of DMY and NC of grass from hyperspectral measurement data. Random forest regressor is a kind of supervised learning algorithm for regression. It uses an ensemble learning method, which combines the output (prediction) from multiple decision trees.

Random forest regressor has less computational costs and it is easier to understand the output since it is based on a decision tree algorithm. The structure of the random forest regressor is shown in Figure 5.

For the purpose of building an estimation model, PLS (partial least squares regression) [15,16,17], MPLS (modified partial least squares regression) [18,19], PCR (principal component regression), and GA [13,15] are generally used in related research. Compared to them, fewer parameters of random forest regressor (number of estimators and max depth of each tree) are also one of the advantages.

“Number of estimators” and “max depth” are set to “100” and “none”, respectively. When “max depth” is “none,” nodes are expanded until all leaves become pure. Since the random forest-based approach is known to be parameter robust [20], the default value of 100 for the number of trees parameter of the software library “scikit-learn” described in the following section were applied.

By using random forest regressor, feature importance which indicates the magnitude of the effect of each feature (i.e., wavelength) on the estimation, can be calculated as the mean and standard deviation of accumulation of the impurity decrease within each tree.

2.3.3. X-Loading Analysis

Loadings in principal component analysis (PCA) are interpreted as the coefficients of the linear combination of the initial variables from which the principal components (PCs) are constructed. Loadings are also known as the values that represent the contribution of each variable (i.e., the importance of each variable) [7,21]. Therefore, X-Loadings analysis can reveal the important variables of principal components (in this case, the wavelength band), which explains the variance of the measured grass samples well.

2.3.4. SHAP Analysis

SHAP (shapley additive explanations) analysis is a method to calculate the contribution rate of each feature of each input value to the output of recognition or regression models [22]. The SHAP value is calculated by the following formula and expressed by each input value, that is, the absorbance magnitude of each wavelength band in this document. Hence, it is possible to explain why the prediction was made for the predicted value. For example, SHAP can visualize not only the importance of features that can be visualized by conventional methods such as PCA, but also the contribution of each data at the level of each variable and the positive or negative effect between features and output. More importantly, SHAP analysis is applicable to all AI models besides random forest regressor. Therefore, it is a means to solve the disadvantage of AI, which is considered to be one of the biggest disadvantages of being a black box.

When given the input $x$ and trained model $f$, the model $f$ is approximated with model $g$ to explain the contribution rate of each input variable. Model $g$ is defined as the following Equation (1):

$$g\left(z^{\prime }\right)={\varphi }_0+{{\displaystyle \sum }}_{i=1}^M{\varphi }_i{z}_i^{\prime }$$

(1)

where ${\varphi }_i$ is the contribution of $i$-th variable of input $x$ and $z_i^{\prime }$ indicates 1 when $i$-th variable is observed, otherwise 0. More details are explained in [22].

For the purpose of investigating the importance of features, feature selection methods such as PCA, random forest regressor, GA or phased regression (PHR) [12,13] are commonly used. GA is an algorithm that uses the evolutionary mechanism of organisms to evolve into highly superior individuals through crossover or mutation, and PHR is a model that performs stepwise regression. However, these methodologies only show the importance of features as an absolute value, no matter if the correlation is positive or negative. By using SHAP analysis, both strong positive and negative values can be calculated as positive or negative SHAP values. Besides, how each feature of each input value affects the output values can be calculated and visualized.

2.4. Computational Environment

All the computations in this paper were performed in the virtual environment of Python 3.7.1, which is built on Cluster Computing Environment of Wageningen University & Research. The operation system is Scientific Linux 7.9. Assigned CPU and memory to the virtual environment were one core and 8 GB, respectively.

For the purpose of statistical analysis of measurement data, the following packages were used:

• RandomForestRegressor (including feature selection),
• sklearn.ensemble (0.24.2)—“https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html (accessed on 9th January 2023)”,
• PCA (including X-loadings),
• sklearn.decomposition (0.24.2)—“https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html (accessed on 9th January 2023)”,
• SHAP,
• shap (0.40.0)—“https://pypi.org/project/shap/ (accessed on 9th January 2023)”.

3. Results & Discussion

3.1. Estimation of DMY and NC

In this paper, random forest regressor was used as one of the AI models, as is simple and easy to understand its structure and behavior. Using this model, the regression model was trained to estimate DMY and NC. As input features, reflection values from a hyperspectral sensor were used. Then, a sensor fusion approach was carried out to improve the estimation performance of the AI model. Data acquired by other sensors, specifically the plant height sensor, were integrated into hyperspectral absorbance data as input to the model. If performance could be improved by the sensor fusion approach, it would lead to the realization of a more accurate model or a more feasible system by combining it with lower cost sensor. Datasets were randomly split into 80% as training data and 20% as test data.

Table 2. Performance of estimation by random forest regressor using statistical parameters of r- squared (r²), root-mean-square deviation (RMSE) and mean absolute error (MAE).

	DMY	DMY (with Height)	NC	NC (with Height)
r2	0.94	0.97	0.88	0.90
RMSE	0.35	0.17	2.59	2.35
MAE	0.23	0.25	1.88	1.68

shows the results of the trained regression model evaluated by test data. R- squared (r²), root-mean-square deviation (RMSE) and mean absolute error (MAE) are shown here. The r² of DMY and NC are 0.94 and 0.88, respectively, and are considered acceptable. By adding height sensor data to input, r² values of DMY and NC are slightly improved by 0.03 and 0.02, respectively, and RMSE and MAE of DMY and NC values are decreased by 30% and 10%, respectively.

Here, the performance of DMY and NC (including crude protein content, which is strongly related to NC) calibration is shown in

Table 3. Performance of calibrations for DMY, NC or crude protein (CP) in related research.

Study	Country	Analyte	Parameters	Sample	r2		RMSE
					DMY	CP, NC	DMY	CP, NC
[13]	Austria, Netherlands	Fresh grass	CP	231	-	0.81	-	85.5 kg CP/ha
[15]	Japan	Fresh grass	CP	100	-	0.85	-	6.46 g/DM kg
[25]	Ireland	Fresh grass	DMY, CP	49	0.86	0.84	9.46 g/kg	20.38 g/DM kg
[26]	Germany	Dried, milled grass	Moisture, CP	1812	0.91	0.84	0.45	0.47
[16]	Chile	Fresh grass	DMY, CP	915	0.93	0.84	11.3 g/kg	22.2 g/DM kg
[17]	Italy	Fresh grass	DMY, CP	100	0.87	0.88	2.75 g/kg	2.14 g/DM kg
[27]	France	Fresh grass	CP	103	-	0.93	-	1.55 g/DM kg
[18]	Chile	Fresh grass	DMY, CP	107	0.99	0.91	6.55 g/kg	18.4 g/DM kg
[28]	USA	Fresh grass	NC	31	-	0.88	-	6 g/DM kg
[19]	Ireland	Fresh grass silage	DM, NC	136	0.85	0.78	-	4.8 g/DM kg
[23]	Ireland	Dried, milled grass	CP	2076	-	0.98	-	-
[24]	Ireland	Dried, milled grass	CP	153	-	0.96	-	-

. From this table, the performance of this study regarding DMY (with/without height) estimation is comparable to or better than the related research, which is from 0.85 at the minimum and 0.9 or higher. Since the sensors used can be easily mounted to tractors and can measure from a distance of 1 m, this performance is considered acceptable. Regarding NC, compared to related research, the result of statistical parameters of r² and RMSE are considered better, especially when combined with height sensor information, except compared to research that used dried and milled samples [23,24]. Most of the other research shows from 0.8 and up to around 0.9 of r². Figure 6 shows the pair plots between ground truth and predicted values by random forest regressor of DMY and NC.

3.2. Wavelength Analysis

In this section, by using several methodologies, wavelengths that are important to describe the characteristics or composition of grass were investigated. One approach is based on PCA, which is one of the most common approaches in an unsupervised method for analyzing multivariate analysis. The other approach is based on the regression model, which is trained to estimate specific objective values.

There are related papers that analyzed the hyperspectral wavelength in grassland. The well-known wavelengths relevant to certain grass compositions and those revealed to be important in this section are listed in Appendix B. The similarity or differences between them are discussed here.

3.2.1. PCA Based Approach

Figure 7 shows the contribution rate and the cumulative contribution rate of components from PCA. From this figure, 90% of information is covered by three PCs, and 95% is covered by four PCs. Then, X-loadings were calculated of three principal components with the highest contribution rate.

Figure 8 shows the X-loadings value of the top three PCs, which covers 90% of the overall information. There are three peaks at around 530 nm, 700 nm, and over 900 nm, which are visible green, red edge, and NIR wavelengths, respectively. The list of bands known as absorption wavelengths of pasture components in the visible and near-infrared regions is shown in

Table A1. Contributed wavelength to estimation model, well-known wavelength and their chemical/biochemical characteristics.

Wavelength (nm)	Electron Transition/Bond Vibration/Red Edge	Biochemical Component	Experiment in This Research	References
430	Electron transition	Chlorophyll a	-	[32]
460	Electron transition	Chlorophyll b	-	[32]
523	Electron transition	-	3.2.1. PC1/2	-
530	Electron transition	-	3.2.1 X-loadings	-
532	Electron transition	-	3.2.1 PC1/2	-
535	Electron transition	Crude protein	-	[13]
539	Electron transition	-	3.2.1 PC1/2	-
542	Electron transition	-	3.2.1 PC1/2	-
545	Electron transition	Crude protein	-	[13]
609	Electron transition	Chlorophyll	-	-
612	Electron transition	-	3.2.2. NC estimation	-
680	Red edge	Chlorophyll Neutral detergent fiber	-	[32,33]
698	Red edge	-	3.2.1 PC1/2	-
699	Red edge	-	3.2.1 PC1/2	-
705	Red edge	Neutral detergent fiber	-	[32]
707	Red edge	Nitrogen	-	[12]
710	Red edge	-	3.2.2 NC estimation
711	Red edge	-	3.2.2 NC estimation
721	Red edge	Nitrogen	-	[12]
734	Red edge	-	3.2.1 PC1/2	-
736	Red edge	-	3.2.1 PC1/2	-
895	-	-	3.2.2 NC estimation
910	C-H stretch, 3rd overtone	Protein	-	[32]
912	C-H stretch, 3rd overtone	-	3.2.2 DMY estimation
930	C-H stretch, 3rd overtone	Lipid	-	[32]
940	C-H stretch, 3rd overtone	-	3.2.2 NC estimation X-loadings (PC1)
960	O-H bend, 1st overtone	-	3.2.2 DMY estimation
970	O-H bend, 1st overtone	Water, starch	-	[32]
990	O-H stretch, 2nd overtone	Starch	-	[32]
991	O-H stretch, 2nd overtone	-	3.2.2 DMY estimation	-

. Around 530 nm, where PC1, 2, and 3 have some peaks, there are crude protein wavelength bands at 535 nm and 545 nm. The wavelength band near the red edge 700–740 nm, where the peaks are seen in PC1 and 2, contains wavelength bands of major pasture components such as chlorophyll, neutral detergent fiber, and nitrogen and is commonly used in vegetation indices such as the normalized difference red edge index (NDRE) as an index for monitoring the vigor level of plants. The band of 900 nm and higher, which has many small peaks, mainly for PC2 and 3, has wavelength bands for proteins, lipids, and sugars, in addition to being highly absorbed by water.

Wavelength bands known to be sensitive to major grass constituents, as described above, were also detected as important in the X-loading analysis, and this result is considered reasonable since the grass samples collected throughout the two-year experiment should vary in terms of these compositions.

Table 4. Top 10 wavelength which have high loadings for PC1 and PC2.

Wavelength [nm]	Loading
736	0.485
539	0.396
734	0.384
542	0.371
699	0.262
532	0.217
523	0.214
698	0.213
531	0.197
534	0.187

shows the top 10 wavelengths with high loading values for PC1 and PC2, which explain 85% of all the information. They are divided into wavelengths around 530 nm and 700 nm, which are the range of red edge and visible green, respectively. In the vicinity of 736 nm and 734 nm, which showed high loadings, there is a wavelength band of 707 nm, which is known to be sensitive to nitrogen. As mentioned above, it is also a red-edge wavelength band widely used as a vegetation value. Around the neighborhood of 539 nm and 542 nm, which have high loads, there are wavelength bands of 535 nm and 545 nm that respond with crude protein. From the viewpoint of PCA data separation, components such as nitrogen, chlorophyll, and crude protein content vary greatly depending on the grass samples contributing to the distribution of the data. Therefore, the result is considered to be understandable.

3.2.2. AI Model-Based Approach

By analyzing the structure of the random forest regressor model trained to estimate DMY and NC, the important wavelength for this purpose is revealed.

Figure 9 shows the importance of features to estimate DMY with/without height sensor data. For this purpose, the NIR wavelength range above 900 nm is judged as important both with and without a height sensor. Especially observed wavelength bands of 954 nm, 961 nm, 962 nm, 990 nm, and 991 nm are around 960 nm and 990 nm, which are reported to be related to starch and moisture. This is a reasonable result since starch is an important index which is related to grass growth, and water content obviously has a strong relation to DMY. Height sensor data also showed some amount of importance, and this is intuitively understandable.

Figure 10 represents the importance of features to estimate NC with/without height sensor data. In contrast to DMY estimating model, the red edge range strongly contributes to the estimation. NIR wavelengths and visible orange range around 600 nm also show high importance. This is a persuasive result since the wavelength bands determined to be important here are quite close to 705 nm and 721 nm, which are reported to be directly related to nitrogen, and wavelengths around 609 nm, 680 nm, and 705 nm, which are highly related to chlorophyll and neutral detergent fibers. When combined with height sensor data, this measurement data showed the greatest importance over the other wavelengths. The difference in importance between plant height and absorbance data is very large, but it is not an incorrect result simply because the amount of grass should be proportional to the amount of NC.

The top 20 SHAP values of the random forest regressor model, which was trained to estimate DMY (without height sensor data) are shown in Figure 11. In DMY estimation, NIR range is more relevant than the other wavelength range. The wavelengths chosen here are equivalent to the wavelengths chosen based on the feature importance of the random forest regressor. However, it also implies that a higher reflection value around 960 nm that is known to be relevant to water content, which means less water content, tends to result in less DMY. In contrast, a lower reflection value of 916 nm which is assigned to a higher amount of protein content, shows a negative relation to DMY. A higher reflection value of 916 nm, which means less protein content, can be both positive and negative in DM estimation, since red dots are observed for both positive and negative SHAP ranges.

Thus, unlike other methods such as random forest regressor, the biggest advantage of introducing SHAP is that it is possible to clearly indicate how “each feature value” of “each measurement data”, that is, the absorbance at each wavelength, affects the estimated target value.

Figure 12 shows the SHAP value transition over the feature value of random forest regressor to estimate DMY. Here, the features of 962 nm and 916 nm, which have high SHAP value contributions while having different trends, are chosen. From figure (a) 962 nm, the negative feature value tends to be higher in relation to the SHAP value, while the positive feature value tends to be a constant negative SHAP value around −0.3. In contrast, from Figure 12b, related to 916 nm, the positive feature value tends to be higher in relation to the SHAP value, while the negative feature value tends to lower the negative SHAP value. In addition, the range of the X-axis is −0.3 to −0.4 in the case of (a), while it is −0.06 to 0.08 in (b). This means that SHAP value varies considerably with little change in feature. Therefore, (b) is considered more sensitive to SHAP value.

In Figure 13, the top 20 SHAP values of the random forest regressor model for estimating NC (without height sensor data) are shown. The features identified as having high relevance by SHAP are also close to the wavelength selected based on the feature importance of the random forest regressor. From this figure, red edge wavelength is found to be relevant in addition to NIR range. Less reflection value around 710 nm and 906 nm that are known to be related to NC and protein content, respectively, in related research, which means a high amount of NC or protein content, also shows the high NC in this analysis. A higher reflection value around 940 nm which is known to be related to lipid content, which means less amount of lipid, shows the relation to higher NC.

Same as Figure 12, Figure 14 exhibits the SHAP value transition over the feature value of Random Forest Regressor to estimate NC. Here, the features of 710nm and 940nm, which have high SHAP value contributions while having different trends, are chosen. From the figure of (a) 710nm, the lower feature value than –0.2 shows a very high relation and sensitivity to the positive SHAP value, while a higher feature value than –0.2 tends to be a constant negative SHAP value around –0.1. In contrast, from the figure of (b) 940nm, positive/negative relation between feature value and SHAP value can be observed. Still, a constant SHAP value over any feature value around zero is also observed. This implies that there are cases where 940nm directly affects SHAP values, and cases where other features rather than 940nm are related to SHAP values.

As stated in the discussion for each result, most wavelengths revealed that strongly contribute to the estimation in this paper can be linked to well-known wavelengths in the related papers shown in Table A1. In addition, some negative relations, such as lipid content against NC are also revealed. More importantly, this SHAP analysis revealed: “how the feature value changes affected the output value of the regressor, which is represented as SHAP value.” This is clearly outstanding from previous research and may encourage the use of AI models that have been black boxes until now. It also provides more detailed analysis and insight into the field of multivariate analysis, such as spectroscopy, insisting on the possibility of finding previously missed relationships between measurement data and estimated values.

4. Conclusions

In this research, using data collected with a hyperspectral sensor from grassland under different fertilization schemes, a regression model was constructed to estimate the DMY and NC of grass. The wavelength bands useful for regression analysis were clarified by X-Loading analysis, feature quantity analysis, and SHAP analysis. From PCA based approach, the visible green wavelength (around 500 nm) and red edge wavelength (around 700 nm) are revealed to represent the overall variability. According to the analysis of the feature importance of the Random Forest Regressor, NIR wavelengths (around 960, 910, and 990 nm) are more sensitive to DMY estimation, while red edge (about 710 nm) and visible orange wavelengths (about 610 nm) are more related to NC estimation.

Moreover, by introducing SHAP analysis, not only the positive effect of features but also negative relations such as protein content (around 906 nm) against DMY estimation or lipid content (940 nm) against NC estimation can be revealed, which is capable of being introduced to other composition analysis as well in the future. This method offers new options for approaches such as spectroscopic analysis and sensor fusion using AI estimation models, which have high estimation accuracy but lack convincing results due to the difficulty of logically explaining their behavior. And this approach might give opportunities for detailed data analysis, which leads to finding previously missed data relationships.

In addition, estimation models of DMY and NC, in which plant height data was added to the input in addition to the hyperspectral sensor data, showed that the estimation accuracy was improved by 30% and 10%, respectively.

In future works, it would be interesting to compare the results which are measured in different environments, such as grasslands where legumes are the main species. If there are differences in the results, investigating the detailed reasons using the results of SHAP analysis will give new insights. In addition, by adding optimal data besides plant height, such as weather data and soil data, it will be possible to build a more accurate or robust measurement model.