The construction of the model can be divided into three parts: (1) The MOSUM method is used to detect mutation points in time series data ( Figures 2A, B ). (2) A model based on encoder decoder was constructed, where both the encoder and decoder are BiLSTM models ( Figures 2C, E ). In order to solve the problem of feature disappearance and information loss caused by covariate feature compression in long sequence data, we introduced an attention mechanism in the model to capture the long-term dependency relationship between input decoder covariates and feature variables ( Figure 2D ). The fully connected layer outputs the predicted results ( Figure 2F ). (3) During the training process, training batches containing mutation points are eliminated through pre-defined functions. This section focuses on the construction process of the TMEAD-BiLSTM model, the overall structure of which is shown in Figure 2 .

Time series data usually contains mutation points, and there is currently no widely recognized method to handle these mutation data points. The currently known processing methods include removing abnormal data or using linear interpolation for filling (Liguori et al., 2021). However, these methods may disrupt the time series properties of the data and affect the accuracy of the analysis. Time series prediction has made remarkable progress in recent years, evolving from traditional statistical methods and machine learning to the latest deep learning techniques, thereby advancing the field of time series analysis (Masini et al., 2023). However, there are still some basic problems in practice, and one of the main problems is the existence of lag difference. Lag differences are deviations or delays that occur in the predicted sequence, which can affect the accuracy and robustness of the model (Samanta et al., 2020). Time series data usually contain abrupt change, and the presence of abrupt change greatly reduces the prediction performance of the model when there is a lag in the prediction process, as shown in Figure 3 . As a forage crop, alfalfa shows a positive correlation between leaf area index (LAI), and final yield throughout its growth cycle. During mowing, alfalfa LAI was generally reaches its maximum value, and LAI mutations occur afterward, significantly impacting the prediction performance of deep learning models.

The MOSUM method is a technique used to detect structural changes in time series, commonly used for detecting outliers (Eichinger and Kirch, 2018). The core principle is to monitor the evolution of data by calculating the cumulative sum within a sliding window. As the window slides, calculate the magnitude of the cumulative sum change within the current window and set a predetermined threshold. When the cumulative sum changes beyond this threshold, the corresponding window position is marked as an outlier. The MOSUM method utilizes the bandwidth parameter (G value) to adjust the size of the sliding window to adapt to different degrees of variation. The mathematical formula of the MOSUM method is as follows:

Where, S_t represents the cumulative sum of the sliding window in the time series; R_t represents the amount of change in the time series at time t; X_i represents the data point in the time series; k represents the size of the sliding window; and G represents the bandwidth parameter.

As a forage crop aimed at obtaining plants, alfalfa has a positive correlation between leaf area index (LAI) and crop growth and final yield during its growth cycle. The maximum LAI of alfalfa during the growth period is generally observed during cutting, which leads to LAI mutations and significantly affects the predictive performance of deep learning models. Therefore, we consider using the MOSUM method to detect mutation points in alfalfa LAI data, as shown in Figure 4 .

The encoder network is a BiLSTM that sequentially transforms the input time series {x₁,…, x_i,…, x_t}. The mapping from x_i to e_h,i and e_b,i can be obtained, as shown in Formula 6–8:

Where, e_h,i and e_b,i are the forward and backward hidden layer outputs of encoder i at time; W _f and W _b represent the weight matrices of forward and backward units, respectively; f₁ and f₂ are units of LSTM; encoder_outputs_i represents the output of the encoder at time step i.

Using weighted attention to compute the weighted average of encoder outputs to generate context vectors:

Where, attn_weights_i represents the attention weights; context represents the context vectors.

The decoder network is also a BiLSTM, d_h,i and d_b,i are the hidden layer outputs of the decoder at time i. As shown in Formula 11, 12, d_h,i and d_b,i can be obtained:

Where, W _f and W _b indicates the weight matrices of forward and backward units, respectively; f₃ and f₄ are units of LSTM; d_h,i and d_b,i respectively represent the forward and backward hidden layer outputs of the decoder at time step i.

BiLSTM can extract complex features of time series, and attention layer output and hidden state information d_f,i-1 、d_b,i+1 can improve the prediction performance. Afterwards, the predicted values are output through the fully connected layer:

Where, output_i represents the prediction result; fc denotes the fully connected layer.

Before the start of training, the MOSUM method is used to detect potential mutation points in the data. When selecting a training batch, predefined functions are used to determine whether to include mutation points in the upcoming training batch. If a mutation point exists, the batch will be skipped until a batch without the mutation point is found. The specific process is shown in Figure 5 .

In detail, the batch size s is one of the most important control parameters of our method. We have to choose s in such a way that it could cover at most one change point. If s is larger than the distance between two change points c_i , c_j , the part of the time series between time index i and j would never be learned. The batches containing this part of the time series would always contain at least one change point and would, therefore, not be allowed. To avoid this, we develop a method that automatically selects the maximum batch size based on the change points. In this method we calculate the distance between all detected change points. We then set s_max to be half of the smallest distance between the change points. Taking half of this distance increases the probability of selecting the allowable batch between the change points, and the resulting number s_max is the maximum batch size that should be selected. See Algorithm 1 for details.

In the TMEAD-BiLSTM model, we pass a batch size s ≤ s_max and the detected change points c₀ ,…, c_k−1 to the algorithm. Our goal is to find the start and end points of the batches that do not contain any change points. To do this, we select indices that lie between 0 and n-s-1. For each change point, we check whether the index of the respective change point lies between the previously selected start index and end index of the batch. If the change point lies in the batch, we repeat the method and select a new start index, and perform the same process for this, as shown in Algorithm 2 . For example, a set of change points has been found through the MOSUM method:

Select batch size s ≤ s_max in Algorithm 1 . Select a random starting index start, that is in the range [0, n-s-1] of the time series. This gives us a batch with start point i_start and end point:

Assume that one of the change points lies in the batch, thus without loss of generality ∃j ∈ {0,…, k-1} such that c_j ∈[i_start , i_end ]. Then again a random start index i_start’ is chosen and:

If now ∀j∈{0,…, k − 1} holds:

In the following batch [i_start’ ,i_end’ ] is the valid batch.

The relationship between the relative leaf area index and the number of growing days is illustrated in Figure 6 . The Logistic, Richards, and Gompertz models were employed to fit the changes in the relative leaf area index over the growing days shown in Figure 6 , respectively. The fitting results of the logistic (Formula 22), Richards (Formula 23), and Gompertz (Formula 24) models are as follows:

Table 2 illustrates the fitting and validation results of three non-linear models for the relative leaf area index. It can be observed that all three models exhibit acceptable fits (R²>0.78), and the models were validated using test data. It is evident that there is no significant difference in the predictive performance of the three models for the alfalfa leaf area index.

Due to differences in alfalfa phenology under different meteorological conditions and soil moisture, growth days, and soil moisture at different depths (>0-10 cm, >10-20 cm, >20-30 cm) corresponding to LAI at time t were also used as input features. To ensure consistency in input features, the features were normalized and scaled to 0-1. Parallel control experiments were conducted between the TMEAD-BiLSTM model and four models: LSTM, BiLSTM, MLSTM and MBiLSTM. The prediction results of each model are compared using RMSE、MAE and R² measures the degree of fitting of the model. The prediction results of each model are shown in Table 3 .

As shown in Figure 7 , the model has a high accuracy in predicting 1.5 ≤ LAI ≤ 5.5, with both predicted and labeled values near the 1:1 line. When LAI < 1.5 and LAI > 5.5, the prediction accuracy is slightly lower compared to the previous situation. This may be because under different growth conditions, the time for alfalfa to turn green and set pods is different, and there are significant differences in LAI changes during the corresponding stages of LAI increase and decrease.

The LSTM network performs well in temporal data processing due to its structural features, but the prediction results are not ideal under the influence of multiple feature factors. In contrast, BiLSTM networks have the ability to perceive forward and backward data, enabling stronger causality between forward and backward data. Based on the data presented in Table 3 , it can be seen that the BiLSTM model has a 27.90% reduction in MAE compared to LSTM. In addition, after adding the MOSUM method, the prediction accuracy of all models was significantly improved, verifying that the MOSUM method can better preserve time series features and effectively reduce the impact of mutation points on model prediction accuracy. To further improve the model’s capacity to mine feature factors and time, this study combines the MOSUM approach with the BiLSTM encoder decoder network model based on attention mechanism. The prediction results show that the R² and RMSE of the TMEAD-BiLSTM model are 0.9986 and 0.0662, respectively.

Due to the influence of growth environment temperature, the growth rates of alfalfa vary between different cuttings (Lloveras et al., 1998). It is challenging to simulate the leaf area index changes for multiple cuttings within a year using nonlinear models. Therefore, we modeled the leaf area index of alfalfa for different cuttings separately. Although the variation trend of alfalfa leaf area index with growth days is basically consistent for the same cutting, there are significant differences in LAI values among the experimental plots. In order to analyze its intrinsic mechanism, the relative leaf area index was used to analyze its common growth characteristics.

Figure 8 shows the relationship between the relative leaf area index (RLAI) and growth days. The RLAI of alfalfa was fitted using the Logistic, Richards and Gompertz models. Tables 4 , 5 presents the fitting effects and fitting coefficients of the models, respectively. The R², RMSE and MAE values of three models were acceptable. Additionally, the maximum LAI of the fourth cutting in the observed region was around 5, and the LAI of the fourth cutting exhibited a declining trend in the later growth period ( Figure 8 ). The Logistic and Gompertz models effectively capture the decline process of the leaf area index of the fourth cutting of alfalfa (R² > 0.9). We utilized the Logistic, Richards, and Gompertz models to analyze the rate of change between the alfalfa leaf area index and growth days, computing the first derivative of the fitting curves for different cuttings. Setting growth days to 15, 20, 25, 30, 35, and 40, we averaged the slopes obtained from the three models. The results indicate that the fastest growth rates of alfalfa in different cutting occur at approximately 32, 23, 20, and 18 days after sowing for different cuttings, which is consistent with the results of Yang et al. (2024a)’s study.

The model obtained was validated using the test dataset, and the validation results are shown in Table 6 . The simulated values of the leaf area index (LAI) based on the Logistic and Gompertz models for alfalfa exhibit a high degree of agreement with the measured values. The R² values for the simulated curves of different cutting times are all greater than 0.9, indicating strong goodness-of-fit. Compared to the Logistic and Gompertz models, the Richards model shows no significant advantage in fitting the data from the first three cuttings, and it performs poorly in fitting the data from the fourth cutting. In summary, the fitting of the Logistic and Gompertz models is satisfactory.

Alfalfa is harvested multiple times annually, and the publicly published dataset used in this study includes data from 3-4 harvests per year. During the growth process of alfalfa, there are differences in the changes in LAI of alfalfa in different cutting. Therefore, there are differences in the accuracy of LAI prediction models for alfalfa in different cutting. As shown in Figure 9 , the alfalfa LAI dataset was divided into four batches, namely the first, second, third, and fourth batches. The accuracy analysis of these four batches was performed using the TMEAD-BiLSTM prediction model.

As shown in Table 7 , the prediction accuracy of the TMEAD-BiLSTM model was higher than that of other experimental models in four different crop cycles. Due to the significant impact of temperature on the growth of alfalfa, the growth periods of the second and third batches of alfalfa are from May to September each year, and the LAI growth rate is relatively high. The model’s ability to capture rapid changes in LAI is insufficient, resulting in a slight decrease in prediction accuracy.