Dynamic Structure Learning through Graph Neural Network for Forecasting Soil Moisture in Precision Agriculture

1 Introduction

Precision agriculture [Zhang et al., 2002] is the science of observing, assessing and controlling agricultural practices such as monitoring soil, crop and climate in a field; detection and prevention of pest and disease attacks; providing a decision support system. It can help optimizing the natural resources needed for agricultural activities and thus ensures an efficient, sustainable and environment friendly development of agriculture. Remote sensing and IoT technology [Liaghat et al., 2010] can be used as effective tools in precision agriculture by providing high resolution satellite imagery data with rich information about crop status, crop and water stresses, and ground truth agricultural information from local sensors and agricultural drones. Thus, availability of historical and realtime data has significantly improved the scope of applying artificial intelligence for precision agricultural practices [Jha et al., 2019].

Soil moisture is an important component of precision agriculture [Zhang et al., 2002] for crop health and stress management, irrigation scheduling, food quality and supply chain. Soil moisture measures the amount of water stored in various depths of soil. Forecasting high resolution and accurate soil moisture well ahead of time helps to save natural resources such as water. Soil moisture can be measured by soil sensors in the field [Hummel et al., 2001] or by physics based land surface models [Rui and Beaudoing, 2011]. However, deploying soil moisture sensors on a vast region is expensive. They also cannot provide forecasts. Physics based soil moisture estimation can be quite accurate, but they need extremely rich set of input features such as different soil properties, landscape information, crop information which are difficult to obtain. These physical models are also computationally heavy, making them infeasible to scale over a larger region.

∗The work was done when Anoushka was an intern at IBM Research and Sambaran was affiliated to IBM Research prior to joining Amazon

†Both Anoushka and Sambaran contributed equally

Soil moisture is also directly related to weather parameters such as temperature, precipitation, and vegetation condition of the field such as Normalized Difference Vegetation Index (NDVI) [Pettorelli, 2013]. There are data driven and machine learning approaches present for soil moisture estimation. Conventional techniques such as Bayesian estimation, random forest and support vector regression are commonly used in industry for their simplicity to model soil moisture with respect to these flexible set of input features [Dasgupta et al., 2019; Ahmad et al., 2010]. However, these approaches essentially treat different locations independently while modeling soil moisture and fail to exploit the rich spatial dependency that soil moisture values over a region exhibit. For example, if there is a rainfall, soil moisture values over the whole region increases. In an agricultural region, soil moisture values exhibits similar patterns with varying crop cycles. There exist methods in machine learning to capture such spatial dependencies through latent variable models [Castro et al., 2012] and supervised deep learning architectures [Zhang et al., 2017]. But they are complex in nature (with large number of parameters) and difficult to apply for soil moisture estimation as getting ground truth soil moisture data is expensive because of the cost of deploying physical sensors.

In this paper, we propose soil moisture forecasting as a semi-supervised learning problem on dynamic graphs (also known as temporal graphs or sequence of graphs) where the structure of the graph and node attributes change over time to capture the spatio-temporal variation of soil moisture. In contrast to existing supervised techniques, please note that semisupervised nature of our solution enables us to train the network on both labeled (data points with ground truth soil moisture) and unlabeled (data points with missing soil moisture) data. This is important because (i) obtaining ground truth soil moisture is expensive and (ii) there can be large amount of missing ground truth soil moisture values because of device and communication failures. A graph can be visualized as a connected (by edges) set of entities (nodes). Intuitively, more the correlation between the soil moisture values in two locations more should be the chance of them to be connected in the graph. Recently, graph representation learning, or more specifically graph neural networks (GNNs) [Wu et al., 2019; Kipf and Welling, 2017; Bandyopadhyay and Peter, 2021] are able to achieve state-of-the-art performance on several learning tasks on graph. GNNs are also proposed to operate on dynamic graphs for traffic forecasting [Li et al., 2018; Chen et al., 2020] and social networks [Pareja et al., 2020]. However, these GNN based approaches cannot be used directly for modeling soil moisture. Unlike social networks where the graph structure is explicitly given [Karrer and Newman, 2011], there is no ground truth graph structure given as an input for the problem of spatio-temporal data modeling. As discussed before, it is possible that the correlation between soil moisture values in two nearby locations is negligible due to multiple hidden factors. Thus, we also need to learn the graph structure of the problem along with the prediction of soil moisture. Recently, research is conducted on graph structure learning using GNNs [Zhu et al., 2021]. However, those approaches are often limited to graphs without any temporal aspect and parameterized over the entries of the discrete adjacency matrix which make the optimization computationally expensive. As discussed in Section 3.3, we tackle these problems in a simple and effective way by introducing multiple regularizers. Following are the novel contributions we make in this paper:

• We pose the problem of soil moisture forecasting as a semi-supervised learning problem on dynamic graphs to jointly capture the varying degrees of spatial and temporal dependencies. As no ground truth graph structure is given as an input, we learn and update the sequence of graph structures in an end-to-end fashion.
• We propose a novel dynamic graph neural network, which is referred as DGLR (Dynamic Graph Learning through Recurrent graph neural network).
• We curate two real world soil moisture datasets. The source code of DGLR and the datasets are available at https://github.com/AnoushkaVyas/DGLR.

2 ProblemFormulation

We are given a set of N locations indexed by the set [N] = {1, 2, · · · ,N} from a geographic region. There are [T] = {1, 2, · · · , T} time steps. For each time step t, ground truth soil moisture (SM) values are given for a subset of locations N t tr ⊆ [N]. Please note that N t1 tr may be different fromN t2 tr for t1 ̸= t2 because of the presence of missing SM values.

We use sti ∈ R to denote the soil moisture at location i and time step t. For each location i and time step t, there are some D input features (such as temperature, relative humidity, precipitation, NDVI, etc.) available which might be useful to predict soil moisture. Let us denote those feature by an attribute vector xti ∈ RD , ∀i ∈ [N], t ∈ [T]. We also assume that the geographic distance between any two locations are given. Let us use dij to denote the distance between two locations i and j. Given all these information, our goal is to predict (forecast) the soil moisture at time step T +1 for each location i ∈ [N], such that the forecasting model considers both past soil moisture and input feature values, and also able to learn and exploit the relation between different locations in their soil moisture values.

In practice, retraining the model for forecasting soil moisture for every time step in future is expensive. So for our experiments in Section 4, we train (including validation) the models roughly on the first 70% to 80% of the time interval and predict soil moisture on the rest.

3 Solution Approach

There are multiple components of our integrated solution DGLR. Please note that all the parameters are trained in an end-to-end fashion. We explain each component below.

3.1 Initial Dynamic Graph Formulation

As explained in Section 1, there is no explicit graph structure given for the problem of soil moisture forecasting. To apply graph neural network in the first iteration, we use the following heuristic to form an initial graph structure. First, for each location i ∈ [N], we form a node (indexed by i) in the graph. Intuitively, if two locations are very close by, there soil moisture values can be more correlated compared to the points which are far away (there are exceptions to this and our learning algorithm is going to handle those). So, we connect any two nodes by an undirected and unweighted edge in the graph if the distance between the two corresponding locations dij is less than some pre-defined threshold θ. We also create a self-loop for every node in the graph. We assign attribute vector xti ∈ RF (as discussed in Section 2) to node i at time t. According to this construction, the link structure of the graph is same across different time steps, but node features are changing. Thus, the constructed graph is dynamic in nature 1 . Let us denote the set of graphs as {G1 , · · · ,GT }, where Gt = (V,Et ,Xt ) is the graph 2 at time step t and V = [N]. We denote the adjacency matrix of Gt as At . We also row-normalize the adjacency matrix by dividing each row with the degree of the corresponding node.

3.2 Temporal Graph Neural Network

Given a dynamic graph, we develop a temporal graph neural network which can generate embedding of a node in each time step and also use that to forecast soil moisture. Each layer of the proposed temporal graph neural network has two major components. There is a self-attention based GNN (with shared parameters) similar to [Veliˇckovi´c et al., 2018] which works on each graph, thus capturing the spatial dependency between the nodes. The updated node embeddings from the GNN is fed to a RNN which connects graphs over different time steps to capture the temporal dependency of soil moisture. The layer is formally discussed below.

1 A dynamic or temporal graph is a sequence of (correlated) graphs over time.

2 Initially, the edge set is same for all the graphs, but they will change during the course of learning.

For a graph Gt , we use a trainable parameter matrix W (shared over graphs ∀t ∈ [T]) to transform the initial feature vector xti as Wxti . As different neighboring locations may have significantly different impact on the soil moisture of a location, we use a trainable attention vector a ∈ R2K to learn the importance αij for any two neighboring nodes in a graph as follows:

where NGt (i) is the neighboring nodes of i (including i itself) in the graph Gt , · represents dot product between the two vectors and || is the vector concatenation. These normalized importance parameters are used to update the node features as:

where Atij is the (i, j)the element of At , i.e., weight of the edge (i, j) in Gt .

The sequence of updated node embeddings hti of a node i over different time steps are fed to a Gated Recurrent Unit (GRU) [Chung et al., 2014], which is a popularly used recurrent neural network. The GRU unit for ith node at time step t takes the GNN output hti (from Equation 2) and the GRU output ht-1 i to update hti , as shown below 3 .

where are update and reset gates of GRU at time t, respectively. W iU ,W iR,W iH , P iU , P iR, P iH ,BiU ,BiR,BiH are trainable parameters of GRU for the node i. Please note that we have not shared the parameters of GRUs for different nodes. This helps to boost the performance as the temporal trend of soil moisture values are quite different in different locations. The above completes the description of a GNN+GRU layer of our temporal graph neural network. The output of the GRU of a layer is fed as input to the next layer of GNN. Two such layers of GNN-GRU pair are used for all our experiments. We denote Ht = [ht1 · · · htN ]Trans , ∀t ∈ [T ].

For soil moisture prediction at time t, node embeddings from the final layer GRU are used. The node embeddings from t - w to t - 1 are used where w is a window of time steps. We concatenate ht-w i , · · · , ht-1 i and pass the vector to a fully connected layer with ReLU activation to predict the soil moisture of ith node at time t, which is denoted as ˆ sti . The window w is a sliding window with stride 1. The loss of soil moisture prediction is calculated as:

3 We use the same notation hti as the output of both GNN and GRU. But they update it sequentially.

3.3 Updating the Dynamic Graph Structures

In Section 3.1, we use a simple heuristic to form the initial graph which may not always reflect the actual spatial dependency of soil moisture values. In this section, we try to learn the link structure between different nodes to improve the quality of the graph such that it facilitates soil moisture prediction. Graph structure learning has received attention in recent literature. But a wide varieties of approaches are only limited to static graphs without temporal aspects and parameterize all the entries of the adjacency matrix (i.e., learning edges between the nodes) which makes the optimization combinatorial in nature [Fatemi et al., 2021; Selvan et al., 2020; Shang et al., 2021]. Our approach is quite feasible as we learn the temporal graph structures via node embeddings by introducing simple regularizers as discussed below.

First, we reconstruct the adjacency matrix ˆ At of the graph Gt by the similarity of the node embeddings. Formally, we obtain ˆ At from the node embedding matrix Ht as:

Thus, Atij (weight of an edge) is the dot product of the two corresponding node embeddings hti = Ht i: and htj = Ht j: . The element-wise use of activation function ReLU(·) ensures that reconstructed edge weights are non-negative and thus can be used in the message passing framework of GNN. Similar to At , we also row-normalize ˆ At by degrees of the nodes. However, the reconstructed graph can be noisy if the initially obtained node embeddings are erroneous. To avoid that, we use following regularization terms.

Graph Closeness. Initially constructed graph with adjacency structure At , ∀t ∈ [T], though not perfect, is easy to explain and intuitive as soil moisture values in close by regions tend to be correlated. So we want the reconstructed graph ˆ At not to move very far away fromAt . To ensure that, we use the following graph closeness regularization:

One can use any distance function between the matrices such as KL divergence or any matrix norms. We use binary cross entropy between At and ˆ At as that produces better results.

Feature Smoothness. As explained in Section 1, features like weather parameters, vegetation type of the field have significant impact on soil moisture. So, locations having similar features often exhibit similar soil moisture values. It is not necessary that these locations have to be close to each other through geographic distance. Thus, we aim to connect such locations in the link structure of the graphs by using the feature smoothness regularizer as shown below:

Minimizing above ensures that for any time step t ∈ [T], if feature vectors for two nodes i and j are similar, i.e., ||xti -xtj ||22 is less, the optimization would try to assign higher values to ˆ Atij . Similarly, if ||xti - xtj ||22 is high for two locations i and j at time step t, feature smoothness would try to lower the weight of the edge (i, j) even if their geographic distance dij is low.

Target Smoothness. Finally, if two nodes have similar target variable (which is soil moisture for this work), message passing between them improves the performance of graph neural network, as observed in [Hou et al., 2020]. For example in a graph dataset, if nodes which are directly connected by an edge also share same node labels, it helps the performance of node classification via message passing GNNs. As soil moisture is continuous in nature, we use the square of the difference of two soil moisture values to find their similarity. We calculate this regularization term only on the training set N t tr for each time step t, where ground truth soil moisture is available.

Again, minimizing target smoothness ensures that less edge weights are assigned to a node pair where nodes have very different soil moisture values. Thus, during the application of GNN, it avoids mixing features which lead to different types of soil moisture values.

3.4 Joint Optimization and Training

There are multiple loss components in the overall solution of DGLR. So, we form the final loss function of DGLR by taking a linear combination of all as shown below.

where α1, α2, α3, α4 ∈ R+ are hyperparameters. WG and ΘR contain the trainable parameters of GNN and RNN layers respectively. In contrast to existing works which solve graph structure learning as an expensive approximation of a combinatorial optimization problem [Jin et al., 2020; Selvan et al., 2020], the variables in the optimization problem in Equation 12 are only the parameters of GNN and RNNs. Please note that we do not minimize Equations 9-11 directly with respect to the reconstructed adjacency matrix ˆ A, rather it is always calculated as ˆ At = ReLU(HtHtTrans) (Equation 8) (followed by row normalisation). Thus, instead of solving task learning (which in this case is predicting soil moisture) and graph structure learning as alternating minimization, we jointly solve both in the framework of neural networks using ADAM [Kingma and Ba, 2014].

Once we obtain the reconstructed adjacency matrix ˆ At after optimizing Equation 12, we use ˆ At , along with node feature matrices Xt , ∀t ∈ [T] for the next iteration of DGLR. The pseudo code of DGLR is presented in Algorithm 1.

Please note that DGLR is semi-supervised in nature. For example, there can be examples in the training where soil moisture value is not present due to some loss of information or faulty sensors, but other feature information is present. In a conventional supervised setup, such unlabeled examples (without ground truth) are ignored in the training. But in DGLR, both labeled and unlabeled examples are used to propagate information through the graph neural networks (Equation 2). Unlabeled examples can also contribute in the training loss function components such as Equations 9 and 10. Because of the semi-supervised nature, DGLR can work well even with lesser amount of labeled data or missing data, as shown in Section 4.3.

Computation of the GNNs takes O(MTDK) time, where M is the average number of edges in the initial and reconstructed graphs over time. GRU takes another O(NTK) time to generate embeddings. Finally, the loss components collectively takes O(N 2TK) time. Hence, each iteration of DGLR takes O((N 2 + MD)TK) time. One can reduce the runtime of graph update regularizers by adding restriction on two nodes to get connected by an edge if their physical distance is more than some threshold, particularly for larger datasets. Besides, forming a graph across locations which are very far from each other does not make sense because of the significant change of input features, landscape, soil types, etc. Hence, DGLR can be used for real-world soil moisture prediction.

4 Experimental Section

4.1 Datasets and Experimental Setup

In this section, we summarize the datasets, baselines and the setup used for experiments 4 . Obtaining ground truth soil moisture is expensive due to the cost of physical sensors and most of the existing works from agriculture use their private and often interpolated datasets. So, we have curated two soil moisture datasets and the associated features from multiple sources. We refer them as Spain and USA. We use 6 input features to each location for each time step for Spain and 15 input features for USA. For Spain, they are NDVI obtained from MODIS, SAR back scattering coefficients VV and VH from Sentinel 1, weather parameters consisting of temperature, relative humidity and precipitation. Locally sensed soil moisture is collected from REMEDHUS [Sanchez et al., 2012] for Spain. For USA, the features (soil temperature, weather data, etc.) and soil moisture data are obtained from the SCAN network.

We have used a diverse set of baselines algorithms like SVR [Drucker et al., 1997] and SVR-Shared (for all the Shared models the parameters being shared for all the stations

4 More details can be found in the longer version of the work at https://arxiv.org/abs/2012.03506

Algorithm 1 DGLR

Input: Soil moisture locations [N] and distances dij, ∀i, j ∈ [N] between them. Temporal attribute matrix Xt = {xti | ∀i ∈ [N]} for each time step t ∈ [T +L]. Ground truth soil moisture values for (a subset of) locations from time 1, · · · , T. Window length w.

Output: Predict soil moisture values of all the locations from T +1, · · · , T +L.

Training

• 1: Construct an initial attributed graph Gt with normalized adjacency matrix At and node attribute matrix Xt for each time t ∈ [T ] where each node represents one location (Section 3.1).
• 2: Set ˜ At = At , ∀t ∈ [T ]
• 3: for pre-defined number of iterations do
• 4: Pass ˜ At and Xt to two layers of GNN+RNN pairs (Section 3.2) to generate all the node embeddings, ∀t ∈ [T ].
• 5: Concatenate ht-w i , · · · , ht-1 i and pass it to a fully connected layer to predict soil moisture ˆ sti, ∀i and ∀t ∈ {w + 1, · · · , T}
• 6: Update the parameters of GNN and RNNs by minimizing the losses at Equation 12.
• 7: Obtain ˆ At by Equation 8, followed by row normalization and set ˜ At = ˆ At , ∀t ∈ [T ]
• 8: end for Soil Moisture Prediction
• 9: Pass the reconstructed graph from the last time step of training, along with attribute matrices Xt-w · · · Xt-1 to the trained GNN+RNN model to predict soil moisture at time t, t = T +1, · · · , T +L.

in a dataset), Spatial-SVR where spatial information in SVR is used by concatenating a node's feature with the features of its nearest k neighbors, ARIMA [Hannan and Rissanen, 1982], RNN [Chung et al., 2014] and RNN-Shared, popular spatio-temporal GNN algorithms like DCRNN [Li et al., 2018] and EvolveGCN [Pareja et al., 2020].

Hyperparameters are tuned by using the validation set and test results are reported at the best validation epoch. The embedding dimension is set to 10 for Spain and 20 for USA. We run the algorithms for 2000 epochs on Spain and 500 epochs on USA. We used three different metrics during the test period of each time series. They are Root Mean Square Error (RMSE), Symmetric Mean Absolute Percentage Error (SMAPE) and Correlation Coefficient.

4.2 Results and Analysis

We run each algorithm 10 times on each dataset and reported the average metrics over all the locations. Table 1 shows the performance of all the baselines, along with DGLR and its variants. First, we can see that algorithms with non-shared parameters across the locations in a dataset performs better than their counterparts with shared parameters. This is because soil moisture patterns in a dataset varies significantly from one location to other. Next, the performance of dynamic GNN based baselines like DCRNN and EvolveGCN are not always better than a deep sequence modeling technique like RNN. This shows that using spatial information (via graph) may not always lead to a better performance, especially when the initially constructed graph is noisy. Finally, we find that

• (a) Feature Smoothness
• (b) Target Smoothness

Figure 1: 1a and 1b show loss components to depict the evolution of the graph structure of Spain over the iterations of DGLR (best seen in color).

DGLR is able to achieve better performance than all the baseline algorithms considered.

Comparing the performance of DGLR with its variants, we can see the importance of different components of DGLR such as non-shared parameters of the GRU layers across locations, the importance of learning the graph structure and the role of feature and target smoothness regularizers.

4.3 Sensitivity toMissing Soil Moisture Values

The datasets that we used here do not have any missing soil moisture values. But for real-world applications, it is possible to have missing soil moisture values in between due to problem in the local sensors or in the communication networks. DGLR, being a semi-supervised approach due to use of graph neural networks, can be trained on labeled and unlabeled data points together. To see the impact of missing ground truth values in the training, we randomly remove p% soil moisture values from the training data, where p ∈ {10, 20, 30}. Table 2 shows the test set performance of DGLR with different proportions of missing soil moisture values in the training. Please note that we still use the input features from those locations with missing soil moisture values and they do participate in the message passing framework of GNN. From Table 2, the performance of DGLR drops very slowly with more amount of missing soil moisture values in the training. In fact, its performance with 20% missing values is still better than most of the baselines in Table 1 with no missing values. So, it is clear that DGLR is quite robust in nature, and thus can be applicable for real-world soil moisture prediction.

4.4 Sensitivity to Initially Constructed Graph Structure

We use a simple heuristic to construct an initial graph in Section 3.1 on which we apply the proposed dynamic GNN to forecast soil moisture and update the graph structures over the iterations. In this section, we check the sensitivity of DGLR on different initially constructed graph. To do this, we vary the threshold parameter θ (Section 3.1). As the value of θ increases, nodes with larger distance will get connected in the initial graph, increasing the overall density. Interestingly, Table 3 shows that the performance of DGLR does not change much with different values of θ. This happens because DGLR also learns the optimal graph structure for soil moisture prediction. Thus, it is less sensitive to the heuristic-based initially constructed graph.

4.5 Evolution of Graph Structure

One key component of DGLR is to learn and update the graph structure as it predicts the soil moisture values (Section 3.3). Feature Smoothness (FS) loss (Eq. 10) measures if nodes having dissimilar features are connected. Target Smoothness (TS) loss (Eq. 11) measures if nodes having dissimilar soil moisture values are connected. For Spain, these losses over different iterations (Step 3 of Alg. 1) of DGLR are shown in Figure 1. Both FS loss and TS loss decrease significantly after the initial iterations and stabilize during the end of training. Thus, DGLR is able to connect nodes with similar features and similar soil moisture values by edges over the training such a way that it also helps to predict the soil moisture better.

5 Discussion and Conclusion

In this work, we have addressed the problem of soil moisture modeling and prediction which is of immense practical importance. We propose a semi-supervised dynamic graph neural network which also learns the temporal graph structures iteratively to predict soil moisture. Our solution is robust in nature with respect to missing ground truth soil moisture values and also able to achieve state-of-the-art performance for data driven soil moisture prediction on real-world soil moisture datasets. We hope that our solution along with the publicly available source code and datasets would help further AI research in this direction. In future, we want to apply this solution for other types of spatio-temporal data with different types of domain specific challenges. We will also check the impact of our solution in broader sustainability aspects such as improving production and quality of crops, preserving natural resources, etc.