# Modeling Soil Temperature for Different Days Using Novel Quadruplet Loss-Guided LSTM

Authors: Xuezhi Wang; Wenhui Li; Qingliang Li; Xiaoning Li
Journal: Computational Intelligence and Neuroscience
 (2022)

---

## Abstract

Soil temperature (Ts), a key variable in geosciences study, has generated growing interest among researchers. There are many factors affecting the spatiotemporal variation of Ts, which poses immense challenges for the Ts estimation. To enrich processing information on loss function and achieve better performance in estimation, the paper designed a new long short-term memory model using quadruplet loss function as an intelligence tool for data processing (QL-LSTM). The model in this paper combined the traditional squared-error loss function with distance metric learning between the sample features. It can zoom analyze the samples accurately to optimize the estimation accuracy. We applied the meteorological data from Laegern and Fluehli stations at 5, 10, and 15 cm depth on the 1st, 5th, and 15th day separately to verify the performance of the proposed soil temperature estimation model. Meanwhile, this paper inputs the variables into the proposed model including radiation, air temperature, vapor pressure deficit, wind speed, air pressure, and past Ts data. The performance of the model was tested by several error evaluation indices, including root mean square error (RMSE), mean absolute error (MAE), Nash-Sutcliffe model efficiency coefficient (NS), Willmott Index of Agreement (WI), and Legates and McCabe index (LMI). As the test results at different soil depths show, our model generally outperformed the four existing advanced estimation models, namely, backpropagation neural networks, extreme learning machines, support vector regression, and LSTM. Furthermore, as experiments show, the proposed model achieved the best performance at the 15 cm depth of soil on the 1st day at Laegern station, which achieved higher WI (0.998), NS (0.995), and LMI (0.938) values, and got lower RMSE (0.312) and MAE (0.239) values. Consequently, the QL-LSTM model is recommended to estimate daily Ts profiles estimation on the 1st, 5th, and 15th days.


---

## Body

## 1. Introduction

Soil temperature (Ts) is a main physical variable of the land surface, which has a direct influence on the atmosphere [1]. Relevant fields including geoscience and forestry application aspects have drawn attention from researchers [2, 3]. In principle, all the interactions in terrestrial ecosystems are companied by Ts variations since they involve energy exchanges. Ts is an essential factor for growing crops that can facilitate the development of the root system by impacting microbial activity, soil decomposition, and fluidity of soil water [4]. In addition, the death of animals and plants produces plenty of carbon substrates and a high volume of greenhouse gases in the soil. Consequently, it results in an increase in Ts, thus expediting carbon dioxide emission to the atmosphere [5]. Therefore, accurate Ts monitoring is crucial for agricultural management and atmosphere environment forecast. However, Ts data in most areas is still measured by using traditional sensors, and the Ts data cannot be collected at different depths [6].Therefore, it can be used to solve some problems in different fields for the study ofTs estimation. The essential environmental factors have a great influence on the accuracy of Ts estimation. At present, Ts is mainly predicted by methods based on physical models and data-driven methods. The physical method is based on the heat conduction model to estimate Ts [7]. Meanwhile, the method is greatly affected by related physical parameters and the scale problem [8].The data-driven methods can explore the internal relationship betweenTs and the surrounding environmental factors for Ts estimation. At present, several predictive models based on machine learning methods are used for estimating Ts [9–14]. For example, ANN is composed of a complex network structure that imitates the structure and function of the brain’s neural network, and it has powerful data processing capabilities. Bilgili applied the multilayer perceptron (MLP) model to adequately describe Ts distribution at a monthly temporal scale from meteorological data [15]. Kisi et al. used three machine learning models to estimate monthly Ts at the soil depth of 5 cm and 10 cm, respectively, and verified the predictive performance of radial basis neural networks performed better than that of generalized regression neural networks and multilayer perceptron models [9]. But generalized regression neural networks had the better performance for deeper depth (50 cm and 100 cm). Kisi et al. applied ANN-based models to predict long-term Ts at a monthly temporal scale, and they found that genetic programming generated the best performance with the meteorological data [16]. Zeynoddin et al. applied a multilayer perceptron (MLP) model to describe daily Ts distribution at three soil depths (5, 10, and 20 cm) from past measurements of Ts [17]. Samadianfard et al. processed the meteorological data such as Ta, W, RH, Rs, Sunshine hours (Sh), and air pressure (Ap) and integrated ANN-based models separately to predict Ts at a monthly temporal scale [18]. Mehdizadeh et al. noted that machine learning models combined with time series models performed better performance than the predictive models based on the single machine learning method or single time series method for predicting Ts at a daily temporal scale [19]. Moazenzadeh et al. proposed SVR with krill herd algorithm (SVR-KHA) method in modeling Ts estimation at different depths (5, 10, 20, 30, 50, and 100 cm), which achieved the best performance, compared to SVR and SVR with firefly algorithm (SVR-FA) [11]. Delbari et al. proposed an SVR-based model to compute daily Ts at three depths (5, 30, and 100 cm) in Iran [12]. The ELM network featured by a single hidden layer could improve the learning speed and accuracy of the algorithm and can model the accurate Ts. Nahvi et al. used the improved ELM model on daily Ts estimation based on the self-adaptive evolutionary method and verified the improved predictive model can estimate the adequate Ts [20]. Sanikhani et al. tested the data from the Mersin station, and the results that show the performance of ELM has the best predictive performance than other predictive models. [14]. Feng et al. tested the loess plateau data with ELM and random forests (RF) and ANN-based models showed that ELM had the better performance for estimation Ts of half-hourly at different soil depths [13].As a time loop neural network containing complex neural network modules, LSTM is used in this paper to solve long-term dependence problems, which can effectively alleviate gradient vanishing through the extraction of required features by the gate control unit. LSTM network [21] can learn long-term and short-term behaviors, and it has seen application in vast areas. By integrating LSTM and SVR, Guo et al. significantly improved the prediction accuracy of abnormal passenger flow fluctuations [22]. In hydrology, Zhang et al. designed a novel LSTM model with the dropout scheme to estimate the depth of the water table [23]. In the atmosphere field, Qing et al. estimated the solar irradiance based on the LSTM network [24]; the results showed the method could avoid the overfitting of the model. Li et al. designed a new GANs-LSTM model and noted that it could serve as an alternative method to estimate Ts [25].This article focuses on the following issues. First, we select the environmental factors which will affectTs estimation. Ts memory can help the predictive model “remember” a warm or cold condition when the anomaly is forgotten by the atmosphere forcing. In addition, recent literature reviews have revealed that the input for prediction models is either the past measurements of Ts or other meteorological information (Ta, W, RH, Rs, Sh, and Ap). Assume that the prediction models are constructed using input combinations of past Ts and other meteorological information; how does the prediction model performance? The second question is about the construction of loss function in LSTM. The predictive model for Ts estimation is a regression predictive modeling problem that involves predicting a real-valued quantity. The loss function is crucial for optimizing the predictive model which could express the degree of difference between predicted and observed Ts; meanwhile, it can optimize the predictive model by updating the weights. Recently, most previous studies in loss function of regression predictive model mainly focused on the distance metric learning between predicted values and real values [26–29]. However, the distance metric learning between the sample features (environmental factors) is usually ignored which has already been successfully applied to image processing [30–32]. To enrich information processing in loss function and further improve the estimation performance, how can we construct a novel loss function by combinations of distance metric learning? The last question is about timescale evaluating for Ts estimation. In previous studies, any evaluation at short-term Ts estimation (half-hourly, hourly, daily timescales) does not consider the timeliness of long-term Ts estimation. However, any evaluation at long-term Ts estimation (monthly timescale) does not include the information of Ts in a small timescale. An ideal decision-support tool for Ts estimation should provide a multifarious decision-making basis. How can we design a prediction scheme at the same timescale evaluation that provides not only the short-term decision-making basis but also the long-term decision-making basis?This paper proposed a novel quadruplet loss function based on the LSTM network that combines traditional squared-error loss function with distance metric learning between the sample features, called QL-LSTM. The traditional squared-error loss function is usually applied to the predictive task with great accuracy. The current limitation of this loss function, however, involves the special variation on Ts based on different predictors. As shown in Figure1, we have made labels according to Ts values. The Ts data which are in the same range are made the same label (the Ts data which are in the range of 8–12°C are labeled as “1”, the Ts data which are in the range of 12–16°C are labeled as “2”, and the Ts data which are in the range of 16–22°C are labeled as “3”). Meanwhile, Ta data are labeled as the same as Ts data. In Figure 1(a), we noticed that the Ts data with the same label are almost within a stable range. However, in the red ellipse, Figure 1(b), we observed that similar Ta values may have different labels (Ts data with similar Ta values may vary considerably). The data with this feature will prevent predictive models from accurately exploring the internal relationship between Ts and the surrounding environmental factors discovering. To address this problem, the idea of triplet loss [33] is considered in this paper. Triplet loss optimization allows the anchor and positive points to accumulate and therefore prevent the negative points and realize the similarity calculation of samples. This approach can enrich processing information of loss function and overcome the disadvantage of the traditional squared-error loss function and further improve the estimation performance.Figure 1
Variations of the daily air temperature (a) and soil temperature (b) at Laegern station (located in Switzerland) during 1st, January 2003–9th April 2003 (100 days).
(a)(b)The main three contributions of this research paper are summarized as follows:(1)
As we know, the proposed method that combined traditional squared-error loss function with distance metric learning between the sample features is a new approach to be used forTs estimation.(2)
Daily-scale prediction scheme was designed to provide the multifarious decision-making basis and was used to estimate theTs on the next 1st, 5th, and 15th day. To achieve this end, we input the meteorological and past Ts data to the estimation model.(3)
Results showed that our QL-LSTM method outperformed the existing advanced methods in most cases.


## 2. Data and Methods

### 2.1. The Framework of Soil Temperature Estimation

The corresponding meteorological dataTs as the input of our QL-LSTM model are obtained from FLUXNET at first. In the meantime, several other advanced models based on data-driven technology (SVR, BPNN, ELM, and LSTM) were considered in Ts estimation. Traditional squared-error loss function and distance metric learning between the sample features were integrated into our model for accurate exploration of the internal relationship between Ts and the surrounding environmental factors. Finally, the comparison of model performance is reflected by five evaluating indicators (RMSE, MAE, NS, WI, and LMI). Figure 2 denoted the flow chart of soil temperature estimation.Figure 2
The flow chart of soil temperature estimation.


### 2.2. Long Short-Term Memory (LSTM) Network

LSTM can process and learn long-term dependence problems. Due to the characteristics of the LSTM network, we use it to explore the internal relationship betweenTs and the surrounding predictors. LSTM controls the transmission state through the gating state, remembers what needs to be remembered, and forgets unimportant information. Figure 3 shows the internal structure of an LSTM cell, and the calculation formula of the LSTM is as follows:(1)it=σWihht−1+Wixxt+Wixct−1+bi,ft=σWfhht−1+Wfxxt+Wfxct−1+bf,ct=ft⊗ct−1+it⊗tanhWchht−1+Wcxxt+bc,ot=σWohht−1+Woxxt+Wocct+bo,ht=ot⊗tanhct,y^t=Wyhht+bo,where xt is the input data, and y^t is the output data; it, ft, and ot denote the input gate, forget gate, and output gate; ct represents the unit status at the current moment; ht is the current output value; σ⋅ and tanh⋅ are the activation functions; W and b denote the weight matrix and bias term.Figure 3
The internal structure of an LSTM cell.


### 2.3. Triplet Loss

Triplet loss is a significant “learning criterion” for optimizing the predictive models, which is applied for adjusting the weight parameters of predictive models, including anchor (Anchor) example, positive (Positive) example, and negative (Negative) example. The similarity calculation of the samples is realized through triplet loss learning, which makes the anchor-to-positive distance smaller than the anchor-to-negative distance. And Figure4 denoted the visual representation of triplet loss.Figure 4
A visual representation of triplet loss.Equation (2) expresses the objective function of triplet loss as follows:(2)∑iNfxia−fxip22−fxia−fxin22+α+,where fxia, fxip, and fxin are the corresponding feature expression obtained by training a parameter in the triplet; α represents the minimum interval between the anchor-to-positive distance and the anchor-to-negative distance; the value of ⋅+ defines the degree of loss.


### 2.4. QL-LSTM Model

Previous analysis shows that LSTM with traditional squared-error loss function could not accurately discover the special relationship betweenTs and surrounding predictors. To address this problem, inspired by the study of triplet loss, we combined a predictive model with distance metric learning between the sample features. As far as we know, the method based on distance metric learning between the sample features has not been used to estimate Ts ever. It must be noted that the distance metric learning between the sample features is first proposed in the field of image processing. However, there is no description of the similarity of samples for Ts estimation. In this paper, the clustering method is used to label samples; thus, distance metric learning between the sample features could be further applied in Ts estimation.The framework of our QL-LSTM is shown in Figure5. Firstly, for its ability to cluster data efficiently and scalability, the Ts data were quantized by the clustering method (called Birch) [34]. In the quantization step, any Ts data quantized to the same label will be defined as similar samples (positive). In contrast, any Ts data quantized to different labels will be defined as the dissimilar samples (negative). It is worthwhile to observe that the number of labels should be neither too large nor too small [35]. Hence, the Calinski Harabasz Score (CH) and Y_Silhouette_score (S) are used to evaluate the quality of the cluster [36]. The larger value of CH or S, the better quality of the clustering results. Second, the labeled data are input into the predictive model (LSTM network). Finally, the weights of the predictive model are updated to reduce the loss based on our quadruplet loss function.Figure 5
The framework of QL-LSTM.We setX=xi,lii=1N as the input data, where li represents Ts labeled as “i” and xi represents the labeled environmental factors. Assume C is the total number of labels, where li∈1,2,3…,C. Then, we project an instance xi onto the estimate Ts by fLSTM.;θ:Rd⟶S1, where fLSTM is an LSTM network parameterized by θ. Let Xici=1Nc be the environmental factors in the i -th labeled samples. Nc represents the total number of samples. We evaluate the similarity between samples through cluster analysis and expect the output of the model closer to the true value.


#### 2.4.1. Hard Sample Mining

Hard sample mining generally refers to hard negative mining. Adding negative sample sets to participate in model training can improve the effectiveness of learning and training and mine hard negatives as much as possible [37, 38]. For each fixed picture, the farthest sample picture and the nearest negative sample picture in a training batch are applied to train the network to enhance the generalization ability of the network, so that the network can learn better representations.Inspired by TriHard loss, we first definexic as the test sample: Pc,iPc,i=xjc|j≠i,Pc,i=Nc−1 is a collection, which includes the samples with the same label; Nc,iNc,i=xjk|k≠c,Nc,i=∑k≠cNk represents the other samples’ collection. xic,yic,Pc,i∗,Nc,i∗ is the quadruplet data set we defined. Pc,i∗ is the positive set; Nc,i∗ represents the negative set, Pc,i∗ and Nc,i∗ represent positive and negative sample pairs, and these tuples form the training sample pairs. The query sample is represented as xic; when Sij+ satisfies the formula (3), xi,xj is the pair that we need.(3)Sij+>minxk∈Pc,iSik+μ,where Sij+fLSTMxi;θ,fLSTMxj;θ represents the similarity between two samples, where ⋅,⋅ represents the calculation of an n×n similarity matrix. Sij is the element in S at xi,xj, and μ as a hyperparameter impacts the quadruplet that can control the number of hard positive samples. The condition for selecting a hard and negative pair is the same as above:(4)Sij−<maxxk∈Nc,iSik−μ.


#### 2.4.2. Optimization Objective

For each test samplexic, we use the margin m to make it as close to the positive set Pc,i as possible, and as far away from the negative set Nc,i as possible. All the nontrivial positive points in Pc,i are pulled together by minimizing:(5)Lpxic;fLSTM=12∑xic∈Pc,i∗fLSTMxic−fLSTMxjc−r−m+2,where fLSTMxic and fLSTMxjc denote the estimated Ts of samples xic and xjc, respectively, and fLSTMxic−fLSTMxjc is the Euclidean distance between fLSTMxic and fLSTMxjc. Similarly, all nontrivial negative points in Nc,i need to push out of the boundary τ, by minimizing:(6)LNxic;fLSTM=12∑xjk∈Nc,i∗τ−fLSTMxic−fLSTMxjk+2.Meanwhile, we applied the squared-error loss function to the LSTM model forTs estimation, as follows:(7)LMSExic;fLSTM=12∑xic∈NcfLSTMxic−yic2.In the QL-LSTM, three minimization objectives were put into the model, and they are optimized at the same time:(8)LQL-LSTMxic;fLSTM=LPxic;fLSTM+LNxic;fLSTM+LMSExic;fLSTM,We incorporate stochastic gradient descent and minibatch into the QL-LSTM to optimize the estimation model.xic is a sample of the minibatch, which is obtained by sampling the labels of the training samples randomly, and serves as an anchor. We represent the QL-LSTM of each minibatch as(9)LQL−LSTMX;f=1N∑∀c∀iLQL−LSTMxic;f,where N denotes the batch size. Figure 6 represented the learning procedure of our QL-LSTM model.Figure 6
The learning procedure of QL-LSTM.


### 2.5. Model Training and Testing

The input of our model is the corresponding meteorological data (Ta, W, Ap, Rs, VPD, and Ts) from Laegern and Fluehli stations in Switzerland. And we downloaded the data at https://fluxnet.fluxdata.org/ on FLUXNET with a total of 3,287 patterns from 2006 to 2014. Training datasets had 2465 patterns, and the rest as testing datasets.Comparing our QL-LSTM model with the other advanced methods (SVR, BPNN, ELM, and LSTM), meanwhile, we calculate several evaluation criteria to analyze the model performance, including model fitting degree and the accuracy of the estimation model, as follows:(10)RMSE=∑n=1Nyi−y^i2N,WI=1−∑n=1Nyi−y^i2∑n=1Ny^i−y¯+yi−y¯2,NS=1−∑n=1Nyi−y^i2∑n=1Nyi−y¯2,MAE=∑n=1Nyi−y^iN,LMI=1−∑n=1Nyi−y^i∑n=1Nyi−y¯,where N is the number of the whole data, yi denotes the observed value, y^i is the predicted value, and y¯ is the average of the true values.


### 2.6. Experiments

The data within half an hour is obtained from two meteorological stations in an ecological nature reserve, located in Switzerland, namely, Legern and Fluley. The corresponding meteorological data (Ta, W, Ap, Rs, VPD, and Ts) and past Ts data were input into the models. Meanwhile, the input variables are normalized to eliminate the dimensional influence between indicators. And the formula is as follows:(11)xnorm=x−xminxmax−xmin,where the minimum value of the sample data is represented by xmin, and the maximum value is represented by xmax. Moreover, we have conducted research on the influence of the surrounding environmental factors on the model prediction. And we found the value of Rs in the data of the two stations is low, which is close to the normal distribution.We conducted a statistical analysis of the data from the two stations. Table1 listed the details of variables (minimum value xmin, maximum value xmax, average value xmean, standard deviation zsd, skewness zs, and variation coefficient zv). We used the daily data to verify the performance of the model with every half an hour data. The results showed in Table 1 that Ap had the highest negative skewness and presents a normal distribution at 5 cm depth, which presented similar characteristics in both stations. Meanwhile, zv showed the biggest difference between the two stations. Ts at the 5 cm, 10 cm, and 15 cm depths range −1.888–26.876°C, −0.181–22.193°C, and 0.16–20.826°C, respectively. In summary, results showed that the values of zsd, zs, and zv change very slightly.Table 1
Statistical results of the applied data for Laegern and Fluehli stations.
StationVariablexminxmaxxmeanzsdzszvLaegernTa (°C)−14.50923.6467.8577.084−0.1180.901Rs (W/m2)175.545379.458305.08135.164−0.4180.115VPD (hpa)0.5415.9373.2712.3341.4240.713W (m/s)0.6688.0252.2371.0051.4110.449Ap (kpa)89.87695.16393.2370.714-0.6300.007Ts −5 cm (°C)−1.88826.87610.1046.0610.1030.599Ts −10 cm (°C)−0.18122.1939.7265.435−0.0310.558Ts −15 cm (°C)0.1619.3949.0105.025−0.0680.557FluehliTa (°C)−14.44822.8777.7086.906−0.1000.895Rs (W/m2)194.734377.444306.01532.805−0.3270.107VPD (hpa)0.4169.6622.1291.5431.4100.724W (m/s)0.3424.6361.4760.6190.8940.419Ap (kpa)82.30887.16485.4930.711−0.8430.008Ts −5 cm (°C)−0.3521.8228.7296.3380.0750.726Ts −10 cm (°C)−0.04421.7278.8366.2420.0710.706Ts −15 cm (°C)0.43220.8268.8136.0230.0620.683


## 2.1. The Framework of Soil Temperature Estimation

The corresponding meteorological dataTs as the input of our QL-LSTM model are obtained from FLUXNET at first. In the meantime, several other advanced models based on data-driven technology (SVR, BPNN, ELM, and LSTM) were considered in Ts estimation. Traditional squared-error loss function and distance metric learning between the sample features were integrated into our model for accurate exploration of the internal relationship between Ts and the surrounding environmental factors. Finally, the comparison of model performance is reflected by five evaluating indicators (RMSE, MAE, NS, WI, and LMI). Figure 2 denoted the flow chart of soil temperature estimation.Figure 2
The flow chart of soil temperature estimation.


## 2.2. Long Short-Term Memory (LSTM) Network

LSTM can process and learn long-term dependence problems. Due to the characteristics of the LSTM network, we use it to explore the internal relationship betweenTs and the surrounding predictors. LSTM controls the transmission state through the gating state, remembers what needs to be remembered, and forgets unimportant information. Figure 3 shows the internal structure of an LSTM cell, and the calculation formula of the LSTM is as follows:(1)it=σWihht−1+Wixxt+Wixct−1+bi,ft=σWfhht−1+Wfxxt+Wfxct−1+bf,ct=ft⊗ct−1+it⊗tanhWchht−1+Wcxxt+bc,ot=σWohht−1+Woxxt+Wocct+bo,ht=ot⊗tanhct,y^t=Wyhht+bo,where xt is the input data, and y^t is the output data; it, ft, and ot denote the input gate, forget gate, and output gate; ct represents the unit status at the current moment; ht is the current output value; σ⋅ and tanh⋅ are the activation functions; W and b denote the weight matrix and bias term.Figure 3
The internal structure of an LSTM cell.


## 2.3. Triplet Loss

Triplet loss is a significant “learning criterion” for optimizing the predictive models, which is applied for adjusting the weight parameters of predictive models, including anchor (Anchor) example, positive (Positive) example, and negative (Negative) example. The similarity calculation of the samples is realized through triplet loss learning, which makes the anchor-to-positive distance smaller than the anchor-to-negative distance. And Figure4 denoted the visual representation of triplet loss.Figure 4
A visual representation of triplet loss.Equation (2) expresses the objective function of triplet loss as follows:(2)∑iNfxia−fxip22−fxia−fxin22+α+,where fxia, fxip, and fxin are the corresponding feature expression obtained by training a parameter in the triplet; α represents the minimum interval between the anchor-to-positive distance and the anchor-to-negative distance; the value of ⋅+ defines the degree of loss.


## 2.4. QL-LSTM Model

Previous analysis shows that LSTM with traditional squared-error loss function could not accurately discover the special relationship betweenTs and surrounding predictors. To address this problem, inspired by the study of triplet loss, we combined a predictive model with distance metric learning between the sample features. As far as we know, the method based on distance metric learning between the sample features has not been used to estimate Ts ever. It must be noted that the distance metric learning between the sample features is first proposed in the field of image processing. However, there is no description of the similarity of samples for Ts estimation. In this paper, the clustering method is used to label samples; thus, distance metric learning between the sample features could be further applied in Ts estimation.The framework of our QL-LSTM is shown in Figure5. Firstly, for its ability to cluster data efficiently and scalability, the Ts data were quantized by the clustering method (called Birch) [34]. In the quantization step, any Ts data quantized to the same label will be defined as similar samples (positive). In contrast, any Ts data quantized to different labels will be defined as the dissimilar samples (negative). It is worthwhile to observe that the number of labels should be neither too large nor too small [35]. Hence, the Calinski Harabasz Score (CH) and Y_Silhouette_score (S) are used to evaluate the quality of the cluster [36]. The larger value of CH or S, the better quality of the clustering results. Second, the labeled data are input into the predictive model (LSTM network). Finally, the weights of the predictive model are updated to reduce the loss based on our quadruplet loss function.Figure 5
The framework of QL-LSTM.We setX=xi,lii=1N as the input data, where li represents Ts labeled as “i” and xi represents the labeled environmental factors. Assume C is the total number of labels, where li∈1,2,3…,C. Then, we project an instance xi onto the estimate Ts by fLSTM.;θ:Rd⟶S1, where fLSTM is an LSTM network parameterized by θ. Let Xici=1Nc be the environmental factors in the i -th labeled samples. Nc represents the total number of samples. We evaluate the similarity between samples through cluster analysis and expect the output of the model closer to the true value.


### 2.4.1. Hard Sample Mining

Hard sample mining generally refers to hard negative mining. Adding negative sample sets to participate in model training can improve the effectiveness of learning and training and mine hard negatives as much as possible [37, 38]. For each fixed picture, the farthest sample picture and the nearest negative sample picture in a training batch are applied to train the network to enhance the generalization ability of the network, so that the network can learn better representations.Inspired by TriHard loss, we first definexic as the test sample: Pc,iPc,i=xjc|j≠i,Pc,i=Nc−1 is a collection, which includes the samples with the same label; Nc,iNc,i=xjk|k≠c,Nc,i=∑k≠cNk represents the other samples’ collection. xic,yic,Pc,i∗,Nc,i∗ is the quadruplet data set we defined. Pc,i∗ is the positive set; Nc,i∗ represents the negative set, Pc,i∗ and Nc,i∗ represent positive and negative sample pairs, and these tuples form the training sample pairs. The query sample is represented as xic; when Sij+ satisfies the formula (3), xi,xj is the pair that we need.(3)Sij+>minxk∈Pc,iSik+μ,where Sij+fLSTMxi;θ,fLSTMxj;θ represents the similarity between two samples, where ⋅,⋅ represents the calculation of an n×n similarity matrix. Sij is the element in S at xi,xj, and μ as a hyperparameter impacts the quadruplet that can control the number of hard positive samples. The condition for selecting a hard and negative pair is the same as above:(4)Sij−<maxxk∈Nc,iSik−μ.


### 2.4.2. Optimization Objective

For each test samplexic, we use the margin m to make it as close to the positive set Pc,i as possible, and as far away from the negative set Nc,i as possible. All the nontrivial positive points in Pc,i are pulled together by minimizing:(5)Lpxic;fLSTM=12∑xic∈Pc,i∗fLSTMxic−fLSTMxjc−r−m+2,where fLSTMxic and fLSTMxjc denote the estimated Ts of samples xic and xjc, respectively, and fLSTMxic−fLSTMxjc is the Euclidean distance between fLSTMxic and fLSTMxjc. Similarly, all nontrivial negative points in Nc,i need to push out of the boundary τ, by minimizing:(6)LNxic;fLSTM=12∑xjk∈Nc,i∗τ−fLSTMxic−fLSTMxjk+2.Meanwhile, we applied the squared-error loss function to the LSTM model forTs estimation, as follows:(7)LMSExic;fLSTM=12∑xic∈NcfLSTMxic−yic2.In the QL-LSTM, three minimization objectives were put into the model, and they are optimized at the same time:(8)LQL-LSTMxic;fLSTM=LPxic;fLSTM+LNxic;fLSTM+LMSExic;fLSTM,We incorporate stochastic gradient descent and minibatch into the QL-LSTM to optimize the estimation model.xic is a sample of the minibatch, which is obtained by sampling the labels of the training samples randomly, and serves as an anchor. We represent the QL-LSTM of each minibatch as(9)LQL−LSTMX;f=1N∑∀c∀iLQL−LSTMxic;f,where N denotes the batch size. Figure 6 represented the learning procedure of our QL-LSTM model.Figure 6
The learning procedure of QL-LSTM.


## 2.4.1. Hard Sample Mining

Hard sample mining generally refers to hard negative mining. Adding negative sample sets to participate in model training can improve the effectiveness of learning and training and mine hard negatives as much as possible [37, 38]. For each fixed picture, the farthest sample picture and the nearest negative sample picture in a training batch are applied to train the network to enhance the generalization ability of the network, so that the network can learn better representations.Inspired by TriHard loss, we first definexic as the test sample: Pc,iPc,i=xjc|j≠i,Pc,i=Nc−1 is a collection, which includes the samples with the same label; Nc,iNc,i=xjk|k≠c,Nc,i=∑k≠cNk represents the other samples’ collection. xic,yic,Pc,i∗,Nc,i∗ is the quadruplet data set we defined. Pc,i∗ is the positive set; Nc,i∗ represents the negative set, Pc,i∗ and Nc,i∗ represent positive and negative sample pairs, and these tuples form the training sample pairs. The query sample is represented as xic; when Sij+ satisfies the formula (3), xi,xj is the pair that we need.(3)Sij+>minxk∈Pc,iSik+μ,where Sij+fLSTMxi;θ,fLSTMxj;θ represents the similarity between two samples, where ⋅,⋅ represents the calculation of an n×n similarity matrix. Sij is the element in S at xi,xj, and μ as a hyperparameter impacts the quadruplet that can control the number of hard positive samples. The condition for selecting a hard and negative pair is the same as above:(4)Sij−<maxxk∈Nc,iSik−μ.


## 2.4.2. Optimization Objective

For each test samplexic, we use the margin m to make it as close to the positive set Pc,i as possible, and as far away from the negative set Nc,i as possible. All the nontrivial positive points in Pc,i are pulled together by minimizing:(5)Lpxic;fLSTM=12∑xic∈Pc,i∗fLSTMxic−fLSTMxjc−r−m+2,where fLSTMxic and fLSTMxjc denote the estimated Ts of samples xic and xjc, respectively, and fLSTMxic−fLSTMxjc is the Euclidean distance between fLSTMxic and fLSTMxjc. Similarly, all nontrivial negative points in Nc,i need to push out of the boundary τ, by minimizing:(6)LNxic;fLSTM=12∑xjk∈Nc,i∗τ−fLSTMxic−fLSTMxjk+2.Meanwhile, we applied the squared-error loss function to the LSTM model forTs estimation, as follows:(7)LMSExic;fLSTM=12∑xic∈NcfLSTMxic−yic2.In the QL-LSTM, three minimization objectives were put into the model, and they are optimized at the same time:(8)LQL-LSTMxic;fLSTM=LPxic;fLSTM+LNxic;fLSTM+LMSExic;fLSTM,We incorporate stochastic gradient descent and minibatch into the QL-LSTM to optimize the estimation model.xic is a sample of the minibatch, which is obtained by sampling the labels of the training samples randomly, and serves as an anchor. We represent the QL-LSTM of each minibatch as(9)LQL−LSTMX;f=1N∑∀c∀iLQL−LSTMxic;f,where N denotes the batch size. Figure 6 represented the learning procedure of our QL-LSTM model.Figure 6
The learning procedure of QL-LSTM.


## 2.5. Model Training and Testing

The input of our model is the corresponding meteorological data (Ta, W, Ap, Rs, VPD, and Ts) from Laegern and Fluehli stations in Switzerland. And we downloaded the data at https://fluxnet.fluxdata.org/ on FLUXNET with a total of 3,287 patterns from 2006 to 2014. Training datasets had 2465 patterns, and the rest as testing datasets.Comparing our QL-LSTM model with the other advanced methods (SVR, BPNN, ELM, and LSTM), meanwhile, we calculate several evaluation criteria to analyze the model performance, including model fitting degree and the accuracy of the estimation model, as follows:(10)RMSE=∑n=1Nyi−y^i2N,WI=1−∑n=1Nyi−y^i2∑n=1Ny^i−y¯+yi−y¯2,NS=1−∑n=1Nyi−y^i2∑n=1Nyi−y¯2,MAE=∑n=1Nyi−y^iN,LMI=1−∑n=1Nyi−y^i∑n=1Nyi−y¯,where N is the number of the whole data, yi denotes the observed value, y^i is the predicted value, and y¯ is the average of the true values.


## 2.6. Experiments

The data within half an hour is obtained from two meteorological stations in an ecological nature reserve, located in Switzerland, namely, Legern and Fluley. The corresponding meteorological data (Ta, W, Ap, Rs, VPD, and Ts) and past Ts data were input into the models. Meanwhile, the input variables are normalized to eliminate the dimensional influence between indicators. And the formula is as follows:(11)xnorm=x−xminxmax−xmin,where the minimum value of the sample data is represented by xmin, and the maximum value is represented by xmax. Moreover, we have conducted research on the influence of the surrounding environmental factors on the model prediction. And we found the value of Rs in the data of the two stations is low, which is close to the normal distribution.We conducted a statistical analysis of the data from the two stations. Table1 listed the details of variables (minimum value xmin, maximum value xmax, average value xmean, standard deviation zsd, skewness zs, and variation coefficient zv). We used the daily data to verify the performance of the model with every half an hour data. The results showed in Table 1 that Ap had the highest negative skewness and presents a normal distribution at 5 cm depth, which presented similar characteristics in both stations. Meanwhile, zv showed the biggest difference between the two stations. Ts at the 5 cm, 10 cm, and 15 cm depths range −1.888–26.876°C, −0.181–22.193°C, and 0.16–20.826°C, respectively. In summary, results showed that the values of zsd, zs, and zv change very slightly.Table 1
Statistical results of the applied data for Laegern and Fluehli stations.
StationVariablexminxmaxxmeanzsdzszvLaegernTa (°C)−14.50923.6467.8577.084−0.1180.901Rs (W/m2)175.545379.458305.08135.164−0.4180.115VPD (hpa)0.5415.9373.2712.3341.4240.713W (m/s)0.6688.0252.2371.0051.4110.449Ap (kpa)89.87695.16393.2370.714-0.6300.007Ts −5 cm (°C)−1.88826.87610.1046.0610.1030.599Ts −10 cm (°C)−0.18122.1939.7265.435−0.0310.558Ts −15 cm (°C)0.1619.3949.0105.025−0.0680.557FluehliTa (°C)−14.44822.8777.7086.906−0.1000.895Rs (W/m2)194.734377.444306.01532.805−0.3270.107VPD (hpa)0.4169.6622.1291.5431.4100.724W (m/s)0.3424.6361.4760.6190.8940.419Ap (kpa)82.30887.16485.4930.711−0.8430.008Ts −5 cm (°C)−0.3521.8228.7296.3380.0750.726Ts −10 cm (°C)−0.04421.7278.8366.2420.0710.706Ts −15 cm (°C)0.43220.8268.8136.0230.0620.683


## 3. Results and Discussion

For testing the superiority of our QL-LSTM model performance forTs estimation using scikit-learn, we compared our test results with those of other advanced models (SVR, BPNN, ELM, and LSTM).We choose default parameters for the SVR model. For the BPNN model, the square error is used as the loss function, and the optimization is Adam. The number of samples selected for the model is 400, the iteration is set to 500, the learning rate is set to 5.0e-4, and the size of the nodes is set to 128. The elm function was used to model the ELM model, the sigmoid served to activate the function in the hidden layer, and we set the same size of the nodes to BPNN. Furthermore, we set the hyperparameters of the LSTM to be the same as that of QL-LSTM. As can be seen from Table2 and 3, the different values of the hyperparameter can generate the different predictive results. When the number of samples selected for the model is set to 400, the iteration to 500, the numQL−LSTM to 128, and set the learning rate to 1.0e-3, the QL-LSTM model has the best performance.Table 2
Predictive performance with differentnumQL−LSTM and learning rates at Laegern station.
Learning ratenumQL−LSTMRMSEMAENSWILMI1.0e-5164.7354.0590.17990.3250.093322.9702.5080.6770.8450.439641.3461.0610.9330.9800.7631281.3651.0730.9310.9820.7602561.3411.0490.9340.9830.765161.2300.9620.9440.9860.7851.0e-4321.1750.9150.9490.9870.795641.1670.9030.9500.9870.7981281.1590.8840.9500.9870.8022561.1480.8710.9510.9870.805161.1290.8550.9530.9880.8085.0e-4321.0760.8090.9570.9890.819641.0010.7570.9630.9900.8301280.8170.6290.9750.9930.8592560.8520.6570.9730.9930.853161.0690.8020.9580.9890.8201.0e-3320.9500.7160.9660.9910.839640.8800.6760.9710.9920.8481280.8170.6250.9750.9930.8602560.8250.6340.9750.9930.858160.8710.6750.9720.9930.8495.0e-3320.8580.6610.9730.9930.852640.8680.6610.9720.9930.8521280.8520.6480.9730.9930.8552560.8540.6490.9730.9930.855Table 3
Predictive performance with different numbers of batch size and iterations at Laegern station.
BatchIterationRMSEMAENSWILMI1001000.8600.6660.9720.9930.8512000.8730.6760.9720.9930.8485000.8800.6670.9710.9920.8508001.0730.8250.9570.9890.8151000.8490.6540.9730.9930.8532002000.8170.6250.9750.9930.8605000.8430.6420.9730.9930.8568000.9160.6920.9690.9920.8451000.9730.7390.9650.9910.8343002000.8310.6440.9740.9930.8565000.8170.6260.9750.9930.8608000.8690.6640.9720.9930.8511001.0380.7840.9600.9900.8244002000.8510.6540.9730.9930.8535000.8090.6220.9760.9930.8608000.8170.6280.9750.9930.8591001.0760.8120.9570.9890.8185002000.8680.6660.9720.9930.8515000.8090.6230.9760.9930.8608000.8110.6260.9750.9930.860


### 3.1. Evaluation for the Hyperparameters in Quadruplet Loss Function

The quadruplet loss function has five main hyperparameters, which are the total number of labelsC, hyperparameter μ in equations (3) and (4), and τ and m in equations (5) and (6). When we evaluate the above hyperparameters in the quadruplet loss function, we set the parameters numQL−LSTM to 128, the learning rate to 1.0e-3, the iteration time to 500, and the batch size to 400. We first select the best C based on the Calinski Harabasz Score and Y_Silhouette_score. Figure 7 denotes the Calinski Harabasz Score and Y_Silhouette_score with different numbers of labels. It is observed that both scores achieve the best result when C is 25. Then, we gradually tune the hyperparameters, τ and m. Figure 8 represented the results of the estimation model with different μ, τ, and m in Laegern meteorological station. We can see that when we set μ to be 5.0e-3, τ to be 1.0e-3, and m to be 5.0e-5, and our QL-LSTM model could achieve the best estimation performance (RMSE = 0.789, MAE = 0.605, NS = 0.977, WI = 0.994, and LMI = 0.865). It is probably because the smaller hyperparameters we set, the less hard samples would be computed. Meanwhile, when we set the larger hyperparameters, the more redundant samples would be computed.Figure 7
The Calinski Harabasz Score and Y_Silhouette_score with different numbers of labels.
(a)Figure 8
The estimation results with differentμ, τ, and m at Laegern meteorological station.


### 3.2. The Impact of Different Inputs on the Performance of the Predictive Model

In this part, we analyzed the environmental factors that may affect our QL-LSTM model forTs estimation. Considering that the interaction between different environmental factors would have an impact on the Ts estimation, we combine the meteorological variables accordingly and input them into the submodels we set as follows:(i)
InputI1: Ta (d − 1)(ii)
InputI2: Ta (d − 1) + Rs(d − 1)(iii)
InputI3: Ta (d − 1) + Rs(d − 1) + VPD(d − 1)(iv)
InputI4: Ta (d − 1) + Rs(d − 1) + VPD(d − 1) + W(d − 1)(v)
InputI5: Ta (d − 1) + Rs(d − 1) + VPD(d − 1) + W(d − 1) + Ap(d − 1)(vi)
Output:Ts (d)Then, we consider that the pastTs will continue to have an impact on the future Ts estimation, so we have carried out lag processing for the past Ts on different days, as follows:(i)
InputI6: Ts (d − 1)(ii)
InputI7: Ts (d − 1) + Ts (d − 2)(iii)
InputI8: Ts (d − 1) + Ts (d − 2) + Ts (d − 3)(iv)
InputI9: Ts (d − 1) + Ts (d − 2) +Ts (d − 3) + Ts (d − 4)(v)
InputI10: Ts (d − 1) + Ts (d − 2) +Ts (d − 3) + Ts (d − 4) + Ts (d − 5)(vi)
Output:Ts (d)We input what we specified above into QL-LSTM to predict theTs (d) at the 5 cm depth of the Laegern station. For our model, we first selected the hyperparameters μ as 5.0e-3, τ as 1.0e-3, m as 5.0e-5, C as 25, numQL−LSTM as 128, learning rate as 1.0e-3, iteration time as 500, and batch size as 400, and the results are presented in Table 4. Obviously, the methods of QL-LSTM(I3) and QL-LSTM(I8) are better than the others, respectively. Meanwhile, we could conclude that W (d − 1), Ap (d − 1), Ts (d − 4), and Ts (d − 5) all have an influence on the performance of the predictive model. In addition, by comparing the estimation results between meteorological variables input and past Ts input, we found that our model with past Ts could achieve greater accuracy in modeling than the one with meteorological variables. The reason may be that the predictive model with past Ts input has stronger memory for Ts variables. The estimation of the future Ts should make the best use of its continuity; in this way, we can make a reliable Ts estimation, which not only continues its historical tendency but also conforms to its actual performance. Hence, we construct the predictive model (QL-LSTM(I11)) by combining the environmental factors (Ta (d − 1), Rs(d − 1), VPD(d − 1)) with past Ts (Ts (d − 1), Ts (d − 2), Ts (d − 3)), which is also considered in estimating the Ts (d) at the 5 cm depth of the Laegern station. Experiment results prove that it could achieve the best estimation performance (RMSE = 0.789, LMI = 0.865, WI = 0.994, NS = 0.977, and MAE = 0.605). Hence, the final input for the predictive models is the environmental factors (Ta(d − 1), Rs(d − 1), VPD(d − 1)) and the past Ts (Ts (d − 1), Ts (d − 2), Ts (d − 3)).Table 4
Predictive performance of QL-LSTM at 5 cm depth for the Laegern station.
MethodRMSEMAENSWILMIQL-LSTM(I1)1.4691.1430.9210.9800.744QL-LSTM(I2)1.4441.1130.9230.9810.751QL-LSTM(I3)1.2210.9370.9450.9860.790QL-LSTM(I4)1.3961.0780.9280.9820.759QL-LSTM(I5)1.4541.1430.9220.9800.744QL-LSTM(I6)1.1430.8660.9520.9870.806QL-LSTM(I7)1.0950.8230.9560.9880.815QL-LSTM(I8)1.0770.8150.9570.9890.817QL-LSTM(I9)1.0960.8340.9560.9880.813QL-LSTM(I10)1.1150.8420.9540.9880.811The three methods (QL-LSTM(I3), QL-LSTM(I8), and QL-LSTM(I11)) are used to test the data of the Laegern station. Figure 9 shows the linear relationship between the predicted value and the observed value. The QL-LSTM(I11) model gets the best predictive performance with y = 0.9899x + 0.3022 and the higher R2 (0.9792) compared with the others. In the frequency diagram (Figure 10) of the models (QL-LSTM(I3), QL-LSTM(I8), and QL-LSTM(I11)), the QL-LSTM(I11) also has a higher frequency (91%) compared to the others. Therefore, we can draw a conclusion that the predictive model (a combination of the environmental factors and the past Ts) normally outperformed the other two (by either past measurements of Ts or other meteorological information) in the Ts estimation.Figure 9
The scatterplots of the predictive model testing results (the values of estimated and observed) for the Laegern station. (a) QL-LSTM(I11), (b) QL-LSTM(I8) model, and (c) QL-LSTM(I3) model.
(a)(b)(c)Figure 10
The frequency plot of the predictive models (absolute estimation error) for the Laegern station. (a) QL-LSTM(I11), (b) QL-LSTM(I8) model, and (c) QL-LSTM(I3) model.
(a)(b)(c)


### 3.3. Comparison with Different Models

In this part, our QL-LSTM model was compared with several advanced models, including SVR, BPNN, ELM, and LSTM. The data ofTa, Rs, and VPD on day “d − 1”, and Ts data on different days were acted as input data to different predictive models, and the output was the predicted value of Ts on days “d”, “d + 5”, and “d + 15”. Time steps were in days.The testing results of five different models at 5, 10, and 15 cm depth on the 1st, 5th, and 15th days of the Laegern station were shown in Table5. And we can see that our QL-LSTM model performs better than the existing advanced models at the 5 cm depth on the 1st day. Specifically, the value of RMSE is 0.789, which is reduced relative to 13% (LSTM), 22% (ELM), 28% (BPNN), and 22% (SVR), respectively. The MAE values amount to 0.605 (QL-LSTM), and the others are 0.813 (SVR), 0.872 (BPNN), 0.824 (ELM), and 0.821 (LSTM). Meanwhile, the QL-LSTM model achieved a higher value of NS, WI, and LMI. Hence, it is obvious that our model had the best performance in this case. For the results of 5 cm depth on the 15th day, the LSTM achieved a higher WI (0.892) than estimation from other models on the 15th day, but it is similar to the WI (0.891) of our model. For 10 and 15 cm depth results on the 1st, 5th, and 15th days, the performance of our QL-LSTM model remains stable, although our QL-LSTM model has the lower values of WI (0.952 and 0.933) than the values of WI (0.954 and 0.934) on the LSTM model in individual cases. It can be found that the predictive performance will get better as the soil depth decreases (from 5 cm to 15 cm), but it will decrease as time goes on (from 1st to 15th days). The systematic errors caused this phenomenon for long-term estimation [39].Table 5
The predictive performance with different models at the Laegern station.
Depth (cm)DayMethodRMSEMAENSWILMI5dSVR1.0070.813−1315.50.412−882.6BPNN1.1010.8720.9300.9820.838ELM1.0160.8240.9250.9830.849LSTM0.9100.8210.9260.9830.850QL-LSTM0.7890.6050.9770.9940.865d + 5SVR2.6652.133-1205.30.419-845.1BPNN2.8122.2300.7000.9190.535ELM2.8322.2610.6930.9160.528LSTM2.6432.1130.7300.9260.560QL-LSTM2.4361.9080.7820.9370.573d + 15SVR3.3222.680−1215.20.422−846.9BPNN3.2712.6570.6010.8860.438ELM3.4032.7230.5720.8810.424LSTM3.2782.6520.6010.8920.439QL-LSTM3.0782.4280.6510.8910.45510dSVR1.0110.850−1291.50.414−870.2BPNN0.9860.8310.9230.9830.837ELM1.0900.8620.9160.9810.830LSTM0.9800.8240.9230.9830.839QL-LSTM0.7610.6050.9750.9930.854d + 5SVR2.1731.752−1211.90.417−842.1BPNN2.1881.7510.7740.9350.608ELM2.2911.8330.7620.9380.596LSTM2.1811.7450.7810.9540.615QL-LSTM1.9731.5450.8330.9520.628d + 15SVR2.7942.220−1244.50.431−850.8BPNN2.8312.2630.6520.9090.493ELM2.9202.3000.6320.8990.484LSTM2.7632.1310.6670.9070.509QL-LSTM2.5381.9750.7240.9170.52415dSVR0.5310.450−1304.50.424−873.5BPNN0.5250.4420.9420.9870.921ELM0.6370.4880.9400.9870.913LSTM0.5120.4360.9440.9880.927QL-LSTM0.3120.2390.9950.9980.938d + 5SVR1.7611.400−1262.30.426−857.9BPNN1.7731.4410.8260.9570.667ELM1.9211.5320.8000.9500.643LSTM1.7421.4080.8310.9580.675QL-LSTM1.5331.2030.8820.9680.687d + 15SVR2.3961.937−1266.30.431−857.4BPNN2.4061.9260.7070.9330.537ELM2.5082.0080.6740.9250.518LSTM2.4011.9180.7070.9340.539QL-LSTM2.1891.7190.7600.9330.553The same strategy was applied in the Fluehli station to further verify the performance of the models, with the results shown in Table6. And our QL-LSTM model performs better compared with others. However, for 5 cm depth on the 15th day and 15 cm depth on the 5th day, the BPNN model performs better with the results of RMSE = 2.081, LMI = 0.584, WI = 0.937, NS = 0.769, and MAE = 2.099, and RMSE = 1.832, LMI = 0.726, WI = 0.973, NS = 0.897, and MAE = 1.352 and LMI = 0.726. Our method does not perform well in some cases probably because the weights of the LSTM model are randomly selected to generate the nonoptimal solution. Meanwhile, our novel loss function (quadruplet loss) is applied based on the LSTM model; it only improved estimation performance to a certain extent against the LSTM model. All in all, the results of our testing on the data of different regions show that the performance of our QL-LSTM model is usually better for Ts prediction with different depths and days.Table 6
The predictive performance with different models at the Fluehli station.
Depth (cm)DayMethodRMSEMAENSWILMI5dSVR0.6910.549−1311.50.429−868.5BPNN0.7180.5500.9420.9180.916ELM0.6930.5380.9420.9880.918LSTM0.7230.5490.9410.9870.916QL-LSTM0.4920.3520.9920.9980.930d + 5SVR2.3161.745−1266.60.431−851.7BPNN2.2971.7160.8200.9550.687ELM2.4631.8800.7990.9490.654LSTM2.2911.7180.8210.9560.686QL-LSTM2.0841.5260.8720.9660.698d + 15SVR3.0232.316−1257.50.433−847.4BPNN2.8012.0990.7690.9370.584ELM3.1432.4000.6940.9180.534LSTM2.9712.3530.7230.9220.557QL-LSTM2.8152.1650.7660.9330.57010dSVR0.6280.518−1310.10.428−870.4BPNN0.6350.5080.9440.9880.926ELM0.6210.5100.9440.9880.926LSTM0.6100.5000.9440.9880.928QL-LSTM0.3810.2860.9950.9980.941d + 5SVR2.0431.529−1274.50.430−856.9BPNN2.0211.5100.8460.9630.721ELM2.1531.6000.8300.9580.697LSTM1.7771.5060.9010.9730.732QL-LSTM1.7941.3070.8990.9730.732d + 15SVR2.8132.143−1271.70.432−854.43BPNN2.7922.1120.7380.9330.594ELM2.9262.3000.7170.9260.558LSTM2.8022.1080.7360.9320.597QL-LSTM2.6011.9060.7890.9480.61215dSVR0.6310.506−1291.90.430−862.9BPNN0.6420.5000.9430.9880.926ELM0.4070.4680.9440.9880.930LSTM0.6180.4650.9440.9880.931QL-LSTM0.4090.2900.9940.9980.941d + 5SVR2.0281.551−1269.30.429−853.7BPNN1.8321.3520.8970.9730.726ELM2.1731.6380.8300.9580.691LSTM2.0421.5550.8460.9620.715QL-LSTM1.8351.3540.8960.9720.725d + 15SVR2.7712.156−1257.70.433−848.4BPNN2.7622.1340.7480.9350.595ELM2.8022.2000.7410.9310.579LSTM2.7952.1120.7420.9330.597QL-LSTM2.5521.9310.8050.9490.608


## 3.1. Evaluation for the Hyperparameters in Quadruplet Loss Function

The quadruplet loss function has five main hyperparameters, which are the total number of labelsC, hyperparameter μ in equations (3) and (4), and τ and m in equations (5) and (6). When we evaluate the above hyperparameters in the quadruplet loss function, we set the parameters numQL−LSTM to 128, the learning rate to 1.0e-3, the iteration time to 500, and the batch size to 400. We first select the best C based on the Calinski Harabasz Score and Y_Silhouette_score. Figure 7 denotes the Calinski Harabasz Score and Y_Silhouette_score with different numbers of labels. It is observed that both scores achieve the best result when C is 25. Then, we gradually tune the hyperparameters, τ and m. Figure 8 represented the results of the estimation model with different μ, τ, and m in Laegern meteorological station. We can see that when we set μ to be 5.0e-3, τ to be 1.0e-3, and m to be 5.0e-5, and our QL-LSTM model could achieve the best estimation performance (RMSE = 0.789, MAE = 0.605, NS = 0.977, WI = 0.994, and LMI = 0.865). It is probably because the smaller hyperparameters we set, the less hard samples would be computed. Meanwhile, when we set the larger hyperparameters, the more redundant samples would be computed.Figure 7
The Calinski Harabasz Score and Y_Silhouette_score with different numbers of labels.
(a)Figure 8
The estimation results with differentμ, τ, and m at Laegern meteorological station.


## 3.2. The Impact of Different Inputs on the Performance of the Predictive Model

In this part, we analyzed the environmental factors that may affect our QL-LSTM model forTs estimation. Considering that the interaction between different environmental factors would have an impact on the Ts estimation, we combine the meteorological variables accordingly and input them into the submodels we set as follows:(i)
InputI1: Ta (d − 1)(ii)
InputI2: Ta (d − 1) + Rs(d − 1)(iii)
InputI3: Ta (d − 1) + Rs(d − 1) + VPD(d − 1)(iv)
InputI4: Ta (d − 1) + Rs(d − 1) + VPD(d − 1) + W(d − 1)(v)
InputI5: Ta (d − 1) + Rs(d − 1) + VPD(d − 1) + W(d − 1) + Ap(d − 1)(vi)
Output:Ts (d)Then, we consider that the pastTs will continue to have an impact on the future Ts estimation, so we have carried out lag processing for the past Ts on different days, as follows:(i)
InputI6: Ts (d − 1)(ii)
InputI7: Ts (d − 1) + Ts (d − 2)(iii)
InputI8: Ts (d − 1) + Ts (d − 2) + Ts (d − 3)(iv)
InputI9: Ts (d − 1) + Ts (d − 2) +Ts (d − 3) + Ts (d − 4)(v)
InputI10: Ts (d − 1) + Ts (d − 2) +Ts (d − 3) + Ts (d − 4) + Ts (d − 5)(vi)
Output:Ts (d)We input what we specified above into QL-LSTM to predict theTs (d) at the 5 cm depth of the Laegern station. For our model, we first selected the hyperparameters μ as 5.0e-3, τ as 1.0e-3, m as 5.0e-5, C as 25, numQL−LSTM as 128, learning rate as 1.0e-3, iteration time as 500, and batch size as 400, and the results are presented in Table 4. Obviously, the methods of QL-LSTM(I3) and QL-LSTM(I8) are better than the others, respectively. Meanwhile, we could conclude that W (d − 1), Ap (d − 1), Ts (d − 4), and Ts (d − 5) all have an influence on the performance of the predictive model. In addition, by comparing the estimation results between meteorological variables input and past Ts input, we found that our model with past Ts could achieve greater accuracy in modeling than the one with meteorological variables. The reason may be that the predictive model with past Ts input has stronger memory for Ts variables. The estimation of the future Ts should make the best use of its continuity; in this way, we can make a reliable Ts estimation, which not only continues its historical tendency but also conforms to its actual performance. Hence, we construct the predictive model (QL-LSTM(I11)) by combining the environmental factors (Ta (d − 1), Rs(d − 1), VPD(d − 1)) with past Ts (Ts (d − 1), Ts (d − 2), Ts (d − 3)), which is also considered in estimating the Ts (d) at the 5 cm depth of the Laegern station. Experiment results prove that it could achieve the best estimation performance (RMSE = 0.789, LMI = 0.865, WI = 0.994, NS = 0.977, and MAE = 0.605). Hence, the final input for the predictive models is the environmental factors (Ta(d − 1), Rs(d − 1), VPD(d − 1)) and the past Ts (Ts (d − 1), Ts (d − 2), Ts (d − 3)).Table 4
Predictive performance of QL-LSTM at 5 cm depth for the Laegern station.
MethodRMSEMAENSWILMIQL-LSTM(I1)1.4691.1430.9210.9800.744QL-LSTM(I2)1.4441.1130.9230.9810.751QL-LSTM(I3)1.2210.9370.9450.9860.790QL-LSTM(I4)1.3961.0780.9280.9820.759QL-LSTM(I5)1.4541.1430.9220.9800.744QL-LSTM(I6)1.1430.8660.9520.9870.806QL-LSTM(I7)1.0950.8230.9560.9880.815QL-LSTM(I8)1.0770.8150.9570.9890.817QL-LSTM(I9)1.0960.8340.9560.9880.813QL-LSTM(I10)1.1150.8420.9540.9880.811The three methods (QL-LSTM(I3), QL-LSTM(I8), and QL-LSTM(I11)) are used to test the data of the Laegern station. Figure 9 shows the linear relationship between the predicted value and the observed value. The QL-LSTM(I11) model gets the best predictive performance with y = 0.9899x + 0.3022 and the higher R2 (0.9792) compared with the others. In the frequency diagram (Figure 10) of the models (QL-LSTM(I3), QL-LSTM(I8), and QL-LSTM(I11)), the QL-LSTM(I11) also has a higher frequency (91%) compared to the others. Therefore, we can draw a conclusion that the predictive model (a combination of the environmental factors and the past Ts) normally outperformed the other two (by either past measurements of Ts or other meteorological information) in the Ts estimation.Figure 9
The scatterplots of the predictive model testing results (the values of estimated and observed) for the Laegern station. (a) QL-LSTM(I11), (b) QL-LSTM(I8) model, and (c) QL-LSTM(I3) model.
(a)(b)(c)Figure 10
The frequency plot of the predictive models (absolute estimation error) for the Laegern station. (a) QL-LSTM(I11), (b) QL-LSTM(I8) model, and (c) QL-LSTM(I3) model.
(a)(b)(c)


## 3.3. Comparison with Different Models

In this part, our QL-LSTM model was compared with several advanced models, including SVR, BPNN, ELM, and LSTM. The data ofTa, Rs, and VPD on day “d − 1”, and Ts data on different days were acted as input data to different predictive models, and the output was the predicted value of Ts on days “d”, “d + 5”, and “d + 15”. Time steps were in days.The testing results of five different models at 5, 10, and 15 cm depth on the 1st, 5th, and 15th days of the Laegern station were shown in Table5. And we can see that our QL-LSTM model performs better than the existing advanced models at the 5 cm depth on the 1st day. Specifically, the value of RMSE is 0.789, which is reduced relative to 13% (LSTM), 22% (ELM), 28% (BPNN), and 22% (SVR), respectively. The MAE values amount to 0.605 (QL-LSTM), and the others are 0.813 (SVR), 0.872 (BPNN), 0.824 (ELM), and 0.821 (LSTM). Meanwhile, the QL-LSTM model achieved a higher value of NS, WI, and LMI. Hence, it is obvious that our model had the best performance in this case. For the results of 5 cm depth on the 15th day, the LSTM achieved a higher WI (0.892) than estimation from other models on the 15th day, but it is similar to the WI (0.891) of our model. For 10 and 15 cm depth results on the 1st, 5th, and 15th days, the performance of our QL-LSTM model remains stable, although our QL-LSTM model has the lower values of WI (0.952 and 0.933) than the values of WI (0.954 and 0.934) on the LSTM model in individual cases. It can be found that the predictive performance will get better as the soil depth decreases (from 5 cm to 15 cm), but it will decrease as time goes on (from 1st to 15th days). The systematic errors caused this phenomenon for long-term estimation [39].Table 5
The predictive performance with different models at the Laegern station.
Depth (cm)DayMethodRMSEMAENSWILMI5dSVR1.0070.813−1315.50.412−882.6BPNN1.1010.8720.9300.9820.838ELM1.0160.8240.9250.9830.849LSTM0.9100.8210.9260.9830.850QL-LSTM0.7890.6050.9770.9940.865d + 5SVR2.6652.133-1205.30.419-845.1BPNN2.8122.2300.7000.9190.535ELM2.8322.2610.6930.9160.528LSTM2.6432.1130.7300.9260.560QL-LSTM2.4361.9080.7820.9370.573d + 15SVR3.3222.680−1215.20.422−846.9BPNN3.2712.6570.6010.8860.438ELM3.4032.7230.5720.8810.424LSTM3.2782.6520.6010.8920.439QL-LSTM3.0782.4280.6510.8910.45510dSVR1.0110.850−1291.50.414−870.2BPNN0.9860.8310.9230.9830.837ELM1.0900.8620.9160.9810.830LSTM0.9800.8240.9230.9830.839QL-LSTM0.7610.6050.9750.9930.854d + 5SVR2.1731.752−1211.90.417−842.1BPNN2.1881.7510.7740.9350.608ELM2.2911.8330.7620.9380.596LSTM2.1811.7450.7810.9540.615QL-LSTM1.9731.5450.8330.9520.628d + 15SVR2.7942.220−1244.50.431−850.8BPNN2.8312.2630.6520.9090.493ELM2.9202.3000.6320.8990.484LSTM2.7632.1310.6670.9070.509QL-LSTM2.5381.9750.7240.9170.52415dSVR0.5310.450−1304.50.424−873.5BPNN0.5250.4420.9420.9870.921ELM0.6370.4880.9400.9870.913LSTM0.5120.4360.9440.9880.927QL-LSTM0.3120.2390.9950.9980.938d + 5SVR1.7611.400−1262.30.426−857.9BPNN1.7731.4410.8260.9570.667ELM1.9211.5320.8000.9500.643LSTM1.7421.4080.8310.9580.675QL-LSTM1.5331.2030.8820.9680.687d + 15SVR2.3961.937−1266.30.431−857.4BPNN2.4061.9260.7070.9330.537ELM2.5082.0080.6740.9250.518LSTM2.4011.9180.7070.9340.539QL-LSTM2.1891.7190.7600.9330.553The same strategy was applied in the Fluehli station to further verify the performance of the models, with the results shown in Table6. And our QL-LSTM model performs better compared with others. However, for 5 cm depth on the 15th day and 15 cm depth on the 5th day, the BPNN model performs better with the results of RMSE = 2.081, LMI = 0.584, WI = 0.937, NS = 0.769, and MAE = 2.099, and RMSE = 1.832, LMI = 0.726, WI = 0.973, NS = 0.897, and MAE = 1.352 and LMI = 0.726. Our method does not perform well in some cases probably because the weights of the LSTM model are randomly selected to generate the nonoptimal solution. Meanwhile, our novel loss function (quadruplet loss) is applied based on the LSTM model; it only improved estimation performance to a certain extent against the LSTM model. All in all, the results of our testing on the data of different regions show that the performance of our QL-LSTM model is usually better for Ts prediction with different depths and days.Table 6
The predictive performance with different models at the Fluehli station.
Depth (cm)DayMethodRMSEMAENSWILMI5dSVR0.6910.549−1311.50.429−868.5BPNN0.7180.5500.9420.9180.916ELM0.6930.5380.9420.9880.918LSTM0.7230.5490.9410.9870.916QL-LSTM0.4920.3520.9920.9980.930d + 5SVR2.3161.745−1266.60.431−851.7BPNN2.2971.7160.8200.9550.687ELM2.4631.8800.7990.9490.654LSTM2.2911.7180.8210.9560.686QL-LSTM2.0841.5260.8720.9660.698d + 15SVR3.0232.316−1257.50.433−847.4BPNN2.8012.0990.7690.9370.584ELM3.1432.4000.6940.9180.534LSTM2.9712.3530.7230.9220.557QL-LSTM2.8152.1650.7660.9330.57010dSVR0.6280.518−1310.10.428−870.4BPNN0.6350.5080.9440.9880.926ELM0.6210.5100.9440.9880.926LSTM0.6100.5000.9440.9880.928QL-LSTM0.3810.2860.9950.9980.941d + 5SVR2.0431.529−1274.50.430−856.9BPNN2.0211.5100.8460.9630.721ELM2.1531.6000.8300.9580.697LSTM1.7771.5060.9010.9730.732QL-LSTM1.7941.3070.8990.9730.732d + 15SVR2.8132.143−1271.70.432−854.43BPNN2.7922.1120.7380.9330.594ELM2.9262.3000.7170.9260.558LSTM2.8022.1080.7360.9320.597QL-LSTM2.6011.9060.7890.9480.61215dSVR0.6310.506−1291.90.430−862.9BPNN0.6420.5000.9430.9880.926ELM0.4070.4680.9440.9880.930LSTM0.6180.4650.9440.9880.931QL-LSTM0.4090.2900.9940.9980.941d + 5SVR2.0281.551−1269.30.429−853.7BPNN1.8321.3520.8970.9730.726ELM2.1731.6380.8300.9580.691LSTM2.0421.5550.8460.9620.715QL-LSTM1.8351.3540.8960.9720.725d + 15SVR2.7712.156−1257.70.433−848.4BPNN2.7622.1340.7480.9350.595ELM2.8022.2000.7410.9310.579LSTM2.7952.1120.7420.9330.597QL-LSTM2.5521.9310.8050.9490.608


## 4. Conclusions

Soil temperature (Ts) is a main physical variable of the land surface, which has an impact on many aspects, such as the growth and yield of crops. Therefore, how to predict Ts accurately is very important. This paper proposed the QL-LSTM model and compared it with the state-of-the-art predictive models to use the meteorological data and past Ts of the Laegern and Fluehli stations (Switzerland) for daily Ts estimation at 5, 10, and 15 cm depth on the 1st, 5th, and 15th days. The experiment results showed that the QL-LSTM model performed better than the existing advanced models for Ts estimation in multifarious cases.In addition, to enrich processing information in loss function and further improve estimation performance, we attempt to design the novel quadruplet loss function that combines the traditional squared-error loss function with distance metric learning between the sample features. Similar samples can be zoomed and the dissimilar samples can be pushed.The distance metric learning between the sample features is combined with the squared-error loss function, which could improve the estimation performance to a certain extent. However, the many hyperparameters in our method may cause sensitivity issues in estimation, which may lead to poor generalization ability of other estimations. In the future, the parametric adaptive method will be explored for a new loss function in the follow-up study.



---

*Source: 9016823-2022-02-17.xml*